text:
function macWithCarry(a::ArbDigit, b::ArbDigit, c::ArbDigit, carry::DoubleArbDigit)
    @assert carry < (one(DoubleArbDigit) << BITS)
    carry += DoubleArbDigit(a) + DoubleArbDigit(b) * DoubleArbDigit(c)
    return carry % ArbDigit, carry >> BITS
end

function mulWithCarry(a::ArbDigit, b::ArbDigit, cin::DoubleArbDigit)
    cout = DoubleArbDigit(a) * DoubleArbDigit(b) + cin
    return cout % ArbDigit, cout >> BITS
end

function macDigit!(acc::AbstractVector{ArbDigit}, b::AbstractVector{ArbDigit}, c::ArbDigit)
    iszero(c) && return
    # we need to make sure that we have enough space
    # meaning length(b) +1 for carry
    @assert length(acc) > length(b) "accumulator needs extra space for carry!"

    carry = zero(DoubleArbDigit)
    a_lo = @view acc[begin:length(b)]
    a_hi = @view acc[length(b)+1:end]

    # we want to add `b*c+carry` into `acc`,
    # so we add elementwise into the lower
    # length(b) digits of `acc`.
    # whatever is left over is our true carry.
    for (a,(i,b)) in zip(a_lo, enumerate(b))
        a_lo[i], carry = macWithCarry(a, b, c, carry)
    end

    carry_hi, carry_lo = fromDoubleArbDigit(carry)

    # if we have carry from multiplying, we add it here
    final_carry = if iszero(carry_hi)
        _add2!(a_hi, carry_lo)
    else
        _add2!(a_hi, [carry_hi, carry_lo])
    end

    # there should never be anything left over
    # i.e. if this assert fires, the one above
    # should also have fired already
    @assert iszero(final_carry) "carry overflow during multiplication!"
end

function long_mul!(acc::AbstractVector{ArbDigit}, b::AbstractVector{ArbDigit}, c::AbstractVector{ArbDigit})
    # naive long multiplication
    # multiply each digit of `c` with `b` and
    # write to `acc`
    for (i,ci) in enumerate(c)
        macDigit!(@view(acc[i:end]), b, ci)
    end
end

function karatsuba!(acc::AbstractVector{ArbDigit}, a::AbstractVector{ArbDigit}, b::AbstractVector{ArbDigit})
    #= karatsuba multiplication
      exponents are just for numbering, not power function
      power is denoted by ^

      write x*y in two parts:
      a = a¹ * sh + a²
      b = b¹ * sh + b²

      now this becomes
      a*b = (a¹ * sh + a²) * (b¹ * sh + b²)
          =  a¹ * b¹            * sh^2
          + (a² * b¹ + a¹ * b²) * sh
          +  a² * b²

      now set p¹ = a¹ * b¹ and p³ = a² * b²:
      a*b =  p¹                 * sh^2
          + (a² * b¹ + a¹ * b²) * sh
          +  p³

      here comes karatsuba's trick:
              a¹ * b² + a² * b¹
          =   a¹ * b² + a² * b¹ - p¹ + p¹ - p³ + p³
          =   a¹ * b² + a² * b¹ - a¹ * b¹ - a² * b²  + p¹ + p³
          = -(a¹ * b¹ - a¹ * b² - a² * b¹ + a² * b²) + p¹ + p³
          = -((a² - a¹) * (b² - b¹))                 + p¹ + p³

      so we now have p² = (a² - a¹) * (b² - b¹) and thus
      a*b = p¹             * sh^2
          + (p¹ + p³ - p²) * sh
          + p³

      so our intermediate products are
      p¹  = a¹ * b¹
      p²  = (a² - a¹) * (b² - b¹)
      p³  = a² * b²

      which we can use to rearrange the above like so
      a*b = p¹ * sh^2
          + p¹ * sh
          + p³ * sh
          - p² * sh
          + p³
          = p³
          + p³ * sh
          + p¹ * sh
          + p¹ * sh^2
          - p² * sh
      with evaluation order from top to bottom. This means,
      by offsetting our index instead of really multiplying
      by `sh`, we can only calculate p¹, p² and p³ once each
      and reuse the result when adding. Subtracting p²*sh at
      the end is done to ensure we don't go negative. An
      informal argument for that is that we're multiplying
      positive integers - that can't ever give a negative
      result and since all transformations above are valid,
      the result has to be positive as well. However, to ensure
      we don't go negative during an intermediate step, we
      make sure to subtract as the last step.
    =#

    short, long = length(a) < length(b) ? (a,b) : (b,a)
    @assert length(short) <= length(long)
    @assert length(acc) >= length(short)+length(long)

    short_half_len = div(length(short), 2)
    short² = @view short[begin:short_half_len]
    short¹ = @view short[short_half_len+1:end]
    long²  = @view  long[begin:short_half_len]
    long¹  = @view  long[short_half_len+1:end]

    # uppers are longer, so we need at most that space
    p_len = length(short¹) + length(long¹)
    p = Vector{ArbDigit}(undef, p_len)
    p_num = ArbUInt(p)

    resize!(p, p_len)
    fill!(p, zero(ArbDigit))

    # p³ = a² * b²
    mac3!(p, short², long²)
    normalize!(p_num) # remove excess zeros
    # acc += p³
    add2!(acc, p)

    # acc += p³*sh
    add2!(@view(acc[short_half_len+1:end]), p)

    # zero buffer
    resize!(p, p_len)
    fill!(p, zero(ArbDigit))

    # p¹ = a¹ * b¹
    mac3!(p, short¹, long¹)
    normalize!(p_num) # remove excess zeros
    # acc += p¹*sh
    add2!(@view(acc[short_half_len+1:end]), p)

    # acc += p¹*sh^2
    add2!(@view(acc[2*short_half_len+1:end]), p)

    # now for p²
    # p²  = (a² - a¹) * (b² - b¹)
    # tmp¹  = a² - a¹
    # tmp²  = b² - b¹
    # this could be negative, so we have to take care
    sign¹, tmp¹ = sub_sign(short², short¹)
    sign², tmp² = sub_sign(long², long¹)

    if sign¹ * sign² > 0
        # both subtractions are positive, so we
        # have to _subtract_ from acc
        resize!(p, p_len)
        fill!(p, zero(ArbDigit))

        mac3!(p, tmp¹, tmp²)
        normalize!(p_num) # remove excess zeros

        # acc -= p²*sh
        sub2!(@view(acc[short_half_len+1:end]), p)

    elseif sign¹ * sign² < 0
        # one subtraction is negative, so we
        # have to _add_ to acc
        mac3!(@view(acc[short_half_len+1:end]), tmp¹, tmp²)

    else
        # one of the subtractions results in 0
        # we don't have to do anything!
    end

    # we're a mutating function - so return nothing
    return nothing
end

function mac3!(acc::AbstractVector{ArbDigit}, b::AbstractVector{ArbDigit}, c::AbstractVector{ArbDigit})
    if !iszero(length(b)) && iszero(first(b))
        nz = findfirst(!iszero, b)
        nz === nothing && return # only zeros!
        b = @view b[nz:end]
        acc = @view acc[nz:end]
    end
    if !iszero(length(c)) && iszero(first(c))
        nz = findfirst(!iszero, c)
        nz === nothing && return # only zeros!
        c = @view c[nz:end]
        acc = @view acc[nz:end]
    end

    short, long = length(b) < length(c) ? (b,c) : (c,b)

    if length(short) <= 32
        long_mul!(acc, long, short)
    else #if length(short) <= 256
        karatsuba!(acc, long, short)
    #else
        # TODO: implement Toom-3  multiplication for numbers larger than 256 digits
        #throw(ArgumentError("Multiplication for numbers with more than 256 digits has not been implemented yet."))
    end
end

function sub_sign(a::AbstractVector{ArbDigit}, b::AbstractVector{ArbDigit})
    if !isempty(a) && iszero(last(a))
        idx = something(findlast(!iszero, a), 0)
        a = @view a[begin:idx]
    end
    if !isempty(b) && iszero(last(b))
        idx = something(findlast(!iszero, b), 0)
        b = @view b[begin:idx]
    end

    ret = if a == b
        (0, ArbDigit[])
    elseif _less(<, a, b)
        retArr = copy(b)
        sub2!(retArr, a)
        (-1, retArr)
    else
        retArr = copy(a)
        sub2!(retArr, b)
        (1, retArr)
    end
    return ret
end

function mul3(x::AbstractVector{ArbDigit}, y::AbstractVector{ArbDigit})
    len = length(x) + length(y)
    prod = zeros(ArbDigit, len)

    mac3!(prod, x, y)
    ArbUInt(prod)
end

function scalar_mul!(a::ArbUInt, b::ArbDigit)
    if iszero(b)
        set_zero!(a)
    elseif isone(b)
        return
    else
        if ispow2(b)
            # this will break with too large numbers
            shl!(a, trailing_zeros(b))
        else
            carry = zero(DoubleArbDigit)
            for i in eachindex(a.data)
                a[i], carry = mulWithCarry(a[i], b, carry)
            end
            !iszero(carry) && push!(a.data, carry % ArbDigit)
        end
    end
end

Base.:(*)(a::ArbUInt, b::ArbUInt) = mul3(a.data, b.data)
Base.:(*)(a::ArbUInt, b::Unsigned) = scalar_mul!(deepcopy(a), ArbDigit(b))
Base.:(*)(b::Unsigned, a::ArbUInt) = scalar_mul!(deepcopy(a), ArbDigit(b))


text:
using DocStringExtensions

"""
# plotbounds

Plot regression line and PI interval with alpha=0.11.

$(SIGNATURES)

### Required arguments
```julia
* `df::DataFrame`                      : DataFrame with observed variables and scaled variables
* `xvar::Symbol`                       : X variable in df
* `yvar::Symbol`                       : Y variable in df
* `nt::NamedTuple`                     : NamedTuple with Stan samples
* `linkvars::Vector{Symbol}`           : Initial 2 Symbols are regression coefficients,
                                         3rd - only for :predict bounds indicates σ.
```
### Optional arguments
```julia
* `fnc = link`                         : Link function
* `fig::AbstractString=""`             : File to store the plot. If "", a plot is returned
* `stepsize::Float64=0.01`             : Stepsize for boundary accuracy
* `rescale_axis=true`                  : Display using un-standardized scale
* `title::AbstractString=""`           : Title for plot
* `lab::AbstractString=""`             : X axis variable label
```

This method is primarily intended to the PI region around the mean line.

The 2nd symbol in `linkvars` needs to match the matrix
```julia
linkvars = [:a, :bM]
```

For other options, :quantile and :hpdi, two parameters suffice (typically the itercept
and slope parameters).

"""
function pbounds(
    df::DataFrame,
    xvar::Symbol,
    yvar::Symbol,
    nt::NamedTuple,
    linkvars::Vector{Symbol};
    fnc::Function=link,
    fig::AbstractString="",
    stepsize=0.01,
    rescale_axis=true,
    lab::AbstractString="",
    title::AbstractString="")

    xbar = mean(df[:, xvar])
    xstd = std(df[:, xvar])
    ybar = mean(df[:, yvar])
    ystd = std(df[:, yvar])

    xvar_s = Symbol(String(xvar)*"_s")
    yvar_s = Symbol(String(yvar)*"_s")

    minx = minimum(df[!, xvar_s])
    maxx = maximum(df[!, xvar_s])
    x_s = minx:stepsize:maxx

    k = size(nt[Symbol(linkvars[2])], 1)
    m = fnc(collect(x_s), k)
    y_s = m * nt[Symbol(linkvars[2])]
    y_s .+= nt[Symbol(linkvars[1])]'
    mul = meanlowerupper(y_s)

    p = plot(;title)
    if rescale_axis
        xrange = rescale(x_s, xbar, xstd)
        mm = rescale(mul.mean, ybar, ystd)
        lm = rescale(mul.lower, ybar, ystd)
        um = rescale(mul.upper, ybar, ystd)
        plot!(p, xrange, mm;
            xlab=String(xvar), ylab=String(yvar),
            ribbon=(mm-lm, um-mm),
            lab)
    else
        plot!(p, xrange, mul.mean;
            xlab=String(xvar), ylab=String(yvar),
            ribbon=(mul.mean-mul.lower, mul.upper-mul.mean),
            lab)
    end

    if fig == ""
        return(p)
    else
        savefig(p, fig)
    end
end

nt = nt7_2s
linkvars = [:a, :bA, :sigma]
fnc=create_observation_matrix
lab="$(size(nt[Symbol(vars[2])], 1))-th degree polynomial"
title="R^2 = $(r2_is_bad(nt, df))"

display(pbounds(df, :mass, :brain, nt, linkvars; stepsize=0.05, fnc, lab, title))

ext: jl


text:
module window

end

max_stars_repo_name: tmaone/xsh

max_issues_repo_name: turing-complete/xsh

max_forks_repo_name: turing-complete/xsh


text:
################################################
############  PYTHON ENVIRONMENT  ##############
################################################

import Pkg
Pkg.activate(dirname(Base.current_project()))
Pkg.precompile()

## Set up Python environment
# Choose own Python environment with OGGM's installation
# Use same path as "which python" in shell
ENV["PYTHON"] = "/home/jovyan/.conda/envs/oggm_env/bin/python3.9" # path in JupyterHub
Pkg.build("PyCall")

using PyCall
using PyPlot # needed for Matplotlib plots

# Import OGGM sub-libraries in Julia
cfg = pyimport("oggm.cfg")
utils = pyimport("oggm.utils")
workflow = pyimport("oggm.workflow")
tasks = pyimport("oggm.tasks")
graphics = pyimport("oggm.graphics")
bedtopo = pyimport("oggm.shop.bedtopo")
# MBsandbox = pyimport("MBsandbox.mbmod_daily_oneflowline") # TODO: fix issue with Python version in Gemini HPC

# Essential Python libraries
np = pyimport("numpy")
xr = pyimport("xarray")
# matplotlib = pyimport("matplotlib")
# matplotlib.use("Qt5Agg")

###############################################
############  JULIA ENVIRONMENT  ##############
###############################################

## Environment and packages
using Distributed
using ProgressMeter
const processes = 10

if nprocs() < processes
    addprocs(processes - nprocs(); exeflags="--project")
end

println("Number of cores: ", nprocs())
println("Number of workers: ", nworkers())

@everywhere begin
    import Pkg
    Pkg.activate(dirname(Base.current_project()))
    Pkg.precompile()
end

@everywhere begin
using Statistics
using LinearAlgebra
using Random
using HDF5
using JLD2
using OrdinaryDiffEq
using DiffEqFlux
using Flux
using Tullio
using RecursiveArrayTools
using Infiltrator
using Plots
using ProgressMeter
using Dates # to provide correct Julian time slices
using PyCall
using ParallelDataTransfer
using BenchmarkTools

###############################################
#############    PARAMETERS     ###############
###############################################

include("helpers/parameters.jl")

###############################################
#############  ODINN LIBRARIES  ###############
###############################################

cd(@__DIR__)
root_dir = dirname(Base.current_project())

### Climate data processing  ###
include("helpers/climate.jl")
end # @everywhere

### Iceflow forward model  ###
# (includes utils.jl as well)
include("helpers/benchmark_iceflow.jl")

###############################
####  OGGM configuration  #####
###############################
cfg.initialize() # initialize OGGM configuration

PATHS = PyDict(cfg."PATHS")  # OGGM PATHS
home_dir = cd(pwd, "../../../../..")
PATHS["working_dir"] = joinpath(home_dir, "Python/OGGM_data")  # Choose own custom path for the OGGM data
PARAMS = PyDict(cfg."PARAMS")

# Multiprocessing
PARAMS["prcp_scaling_factor"], PARAMS["ice_density"], PARAMS["continue_on_error"]
PARAMS["use_multiprocessing"] = true # Let's use multiprocessing for OGGM

###############################################################
###########################  MAIN #############################
###############################################################

# Defining glaciers to be modelled with RGI IDs
# RGI60-11.03638 # Argentière glacier
# RGI60-11.01450 # Aletsch glacier
# RGI60-11.03646 # Bossons glacier
rgi_ids = ["RGI60-11.03638", "RGI60-11.01450", "RGI60-11.03646"]

### Initialize glacier directory to obtain DEM and ice thickness inversion  ###
# Where to fetch the pre-processed directories

# TODO: make this constant
base_url = ("https://cluster.klima.uni-bremen.de/~oggm/gdirs/oggm_v1.4/L1-L2_files/elev_bands")

# Where to fetch the pre-processed directories

# TODO: change to Lilian's version in notebook (ODINN_MB.ipynb)
# use elevation band  flowlines
# gdirs = workflow.init_glacier_directories(rgi_ids, from_prepro_level=2,
#                                           prepro_border=10,
#                                           prepro_base_url=base_url,
#                                           prepro_rgi_version="62")

(@isdefined gdirs) || (gdirs = workflow.init_glacier_directories(rgi_ids, from_prepro_level=3, prepro_border=40))
# gdirs = workflow.init_glacier_directories(rgi_ids, from_prepro_level=3, prepro_border=40)

glacier_filter = 1
gdir = gdirs[glacier_filter]
rgi_id = rgi_ids[glacier_filter]
println("Path to the DEM:", gdir.get_filepath("dem"))

# Obtain ice thickness inversion
if !@isdefined glacier_gd
    list_talks = [
        # tasks.glacier_masks,
        # tasks.compute_centerlines,
        # tasks.initialize_flowlines,
        # tasks.compute_downstream_line,
        tasks.gridded_attributes,
        tasks.gridded_mb_attributes,
        # tasks.prepare_for_inversion,  # This is a preprocessing task
        # tasks.mass_conservation_inversion,  # This gdirsdoes the actual job
        # tasks.filter_inversion_output,  # This smoothes the thicknesses at the tongue a little
        # tasks.distribute_thickness_per_altitude,
        bedtopo.add_consensus_thickness   # Use consensus ice thicknesses from Farinotti et al. (2019)
    ]
    for task in list_talks
        # The order matters!
        workflow.execute_entity_task(task, gdirs)
    end
end

# Plot glacier domain<
# graphics.plot_domain(gdirs)
# graphics.plot_distributed_thickness(gdir)

(@isdefined glacier_gd) || (glacier_gd = xr.open_dataset(gdir.get_filepath("gridded_data")))

# Broadcast necessary variables to all workers
sendto(workers(), glacier_gd=glacier_gd)
sendto(workers(), gdir=gdir)

@everywhere begin
nx = glacier_gd.y.size # glacier extent
ny = glacier_gd.x.size # really weird, but this is inversed
Δx = gdir.grid.dx
Δy = gdir.grid.dy
end # @everywhere

MBsandbox = pyimport("MBsandbox.mbmod_daily_oneflowline")

### Generate fake annual long-term temperature time series  ###
# This represents the long-term average air temperature, which will be used to
# drive changes in the `A` value of the SIA
(@isdefined temp_series) || (const temp_series, norm_temp_series = fake_temp_series(t₁))
# A_series = []
# for temps in temp_series
#     push!(A_series, A_fake.(temps))
# end
# display(Plots.plot(temp_series, xaxis="Years", yaxis="Long-term average air temperature", title="Fake air temperature time series"))
# display(Plots.plot(A_series, xaxis="Years", yaxis="A", title="Fake A reference time series"))

# Determine initial conditions
(@isdefined H₀) || (const H₀ = glacier_gd.consensus_ice_thickness.data) # initial ice thickness conditions for forward model
fillNaN!(H₀) # Fill NaNs with 0s to have real boundary conditions
smooth!(H₀)  # Smooth initial ice thickness to help the solver
(@isdefined B) || (const B = glacier_gd.topo.data .- H₀) # bedrock

# Run forward model for selected glaciers
if create_ref_dataset
    println("Generating reference dataset for training...")

    # Compute reference dataset in parallel
    @everywhere solver = Ralston()
    H_refs, V̄x_refs, V̄y_refs = generate_ref_dataset(temp_series, H₀)

    println("Saving reference data")
    jldsave(joinpath(root_dir, "data/PDE_refs_$rgi_id.jld2"); H_refs, V̄x_refs, V̄y_refs)
end

# Load stored PDE reference datasets
PDE_refs = load(joinpath(root_dir, "data/PDE_refs_$rgi_id.jld2"))

#######################################################################################################
#############################             Train UDE            ########################################
#######################################################################################################

UA = FastChain(
        FastDense(1,3, x->softplus.(x)),
        FastDense(3,10, x->softplus.(x)),
        FastDense(10,3, x->softplus.(x)),
        FastDense(3,1, sigmoid_A)
    )

# θ = initial_params(UA)

θ = load(joinpath(root_dir, "data/benchmark_weights.jld"))["θ_benchmark"]

current_epoch = 1
batch_size = length(temp_series)

multibatch = false

cd(@__DIR__)
const root_plots = cd(pwd, "../../plots")
# Train iceflow UDE in parallel
# First train with ADAM to move the parameters into a favourable space
@everywhere solver = ROCK4()
# @everywhere solver = TRBDF2()

train_settings = (ADAM(0.05), 20) # optimizer, epochs

iceflow_trained = @btime train_iceflow_UDE(H₀, UA, θ, train_settings, PDE_refs, temp_series)


text:
function PrintDivider()
println("************************************************************")
end


text:
using OrdinaryDiffEq, RecursiveArrayTools, DiffEqBase, Base.Test#, ParameterizedFunctions

f = function (t,u)
  - u + sin(t)
end

prob = ODEProblem(f,1.0,(0.0,10.0))

condition= function (t,u,integrator) # Event when event_f(t,u,k) == 0
  t - 2.95
end

affect! = function (integrator)
  integrator.u = integrator.u + 2
end

callback = ContinuousCallback(condition,affect!)

sol = solve(prob,Tsit5(),callback=callback)

f = function (t,u,du)
  du[1] = - u[1] + sin(t)
end

prob = ODEProblem(f,[1.0],(0.0,10.0))

condtion= function (t,u,integrator) # Event when event_f(t,u,k) == 0
  t - 2.95
end

affect! = function (integrator)
  integrator.u = integrator.u + 2
end

callback = ContinuousCallback(condtion,affect!)

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

#=
f = @ode_def BallBounce begin
  dy =  v
  dv = -g
end g=9.81
=#

f = function (t,u,du)
  du[1] = u[2]
  du[2] = -9.81
end

condtion= function (t,u,integrator) # Event when event_f(t,u,k) == 0
  u[1]
end

affect! = nothing
affect_neg! = function (integrator)
  integrator.u[2] = -integrator.u[2]
end

callback = ContinuousCallback(condtion,affect!,affect_neg!,interp_points=100)

u0 = [50.0,0.0]
tspan = (0.0,15.0)
prob = ODEProblem(f,u0,tspan)

sol = solve(prob,Tsit5(),callback=callback,adaptive=false,dt=1/4)

condtion_single = function (t,u,integrator) # Event when event_f(t,u,k) == 0
  u
end

affect! = nothing
affect_neg! = function (integrator)
  integrator.u[2] = -integrator.u[2]
end

callback_single = ContinuousCallback(condtion_single,affect!,affect_neg!,interp_points=100,idxs=1)

u0 = [50.0,0.0]
tspan = (0.0,15.0)
prob = ODEProblem(f,u0,tspan)

sol = solve(prob,Tsit5(),callback=callback_single,adaptive=false,dt=1/4)

#plot(sol,denseplot=true)

sol = solve(prob,Vern6(),callback=callback)
#plot(sol,denseplot=true)
sol = solve(prob,BS3(),callback=callback)

sol33 = solve(prob,Vern7(),callback=callback)

bounced = ODEProblem(f,sol[8],(0.0,1.0))
sol_bounced = solve(bounced,Vern6(),callback=callback,dt=sol.t[9]-sol.t[8])
#plot(sol_bounced,denseplot=true)
sol_bounced(0.04) # Complete density
@test maximum(maximum.(map((i)->sol.k[9][i]-sol_bounced.k[2][i],1:length(sol.k[9])))) == 0

sol2= solve(prob,Vern6(),callback=callback,adaptive=false,dt=1/2^4)
#plot(sol2)

sol2= solve(prob,Vern6())

sol3= solve(prob,Vern6(),saveat=[.5])

## Saving callback

condtion = function (t,u,integrator)
  true
end
affect! = function (integrator) end

save_positions = (true,false)
saving_callback = DiscreteCallback(condtion,affect!,save_positions=save_positions)

sol4 = solve(prob,Tsit5(),callback=saving_callback)

@test sol2(3) ≈ sol(3)

affect! = function (integrator)
  u_modified!(integrator,false)
end
saving_callback2 = DiscreteCallback(condtion,affect!,save_positions=save_positions)
sol4 = solve(prob,Tsit5(),callback=saving_callback2)

cbs = CallbackSet(saving_callback,saving_callback2)
sol4_extra = solve(prob,Tsit5(),callback=cbs)

@test length(sol4_extra) == 2length(sol4) - 1

condtion= function (t,u,integrator)
  u[1]
end

affect! = function (integrator)
  terminate!(integrator)
end

terminate_callback = ContinuousCallback(condtion,affect!)

tspan2 = (0.0,Inf)
prob2 = ODEProblem(f,u0,tspan2)

sol5 = solve(prob2,Tsit5(),callback=terminate_callback)

@test sol5[end][1] < 3e-12
@test sol5.t[end] ≈ sqrt(50*2/9.81)

affect2! = function (integrator)
  if integrator.t >= 3.5
    terminate!(integrator)
  else
    integrator.u[2] = -integrator.u[2]
  end
end
terminate_callback2 = ContinuousCallback(condtion,nothing,affect2!,interp_points=100)

sol5 = solve(prob2,Vern7(),callback=terminate_callback2)

@test sol5[end][1] < 1e-12
@test sol5.t[end] ≈ 3*sqrt(50*2/9.81)

condtion= function (t,u,integrator) # Event when event_f(t,u,k) == 0
  t-4
end

affect! = function (integrator)
  terminate!(integrator)
end

terminate_callback3 = ContinuousCallback(condtion,affect!,interp_points=100)

bounce_then_exit = CallbackSet(callback,terminate_callback3)

sol6 = solve(prob2,Vern7(),callback=bounce_then_exit)

@test sol6[end][1] > 0
@test sol6.t[end] ≈ 4


text:
# syntax: proto2
using Compat
using ProtoBuf
import ProtoBuf.meta
import ..hadoop

# service methods for ApplicationClientProtocolService
const _ApplicationClientProtocolService_methods = MethodDescriptor[
        MethodDescriptor("getNewApplication", 1, hadoop.yarn.GetNewApplicationRequestProto, hadoop.yarn.GetNewApplicationResponseProto),
        MethodDescriptor("getApplicationReport", 2, hadoop.yarn.GetApplicationReportRequestProto, hadoop.yarn.GetApplicationReportResponseProto),
        MethodDescriptor("submitApplication", 3, hadoop.yarn.SubmitApplicationRequestProto, hadoop.yarn.SubmitApplicationResponseProto),
        MethodDescriptor("forceKillApplication", 4, hadoop.yarn.KillApplicationRequestProto, hadoop.yarn.KillApplicationResponseProto),
        MethodDescriptor("getClusterMetrics", 5, hadoop.yarn.GetClusterMetricsRequestProto, hadoop.yarn.GetClusterMetricsResponseProto),
        MethodDescriptor("getApplications", 6, hadoop.yarn.GetApplicationsRequestProto, hadoop.yarn.GetApplicationsResponseProto),
        MethodDescriptor("getClusterNodes", 7, hadoop.yarn.GetClusterNodesRequestProto, hadoop.yarn.GetClusterNodesResponseProto),
        MethodDescriptor("getQueueInfo", 8, hadoop.yarn.GetQueueInfoRequestProto, hadoop.yarn.GetQueueInfoResponseProto),
        MethodDescriptor("getQueueUserAcls", 9, hadoop.yarn.GetQueueUserAclsInfoRequestProto, hadoop.yarn.GetQueueUserAclsInfoResponseProto),
        MethodDescriptor("getDelegationToken", 10, hadoop.common.GetDelegationTokenRequestProto, hadoop.common.GetDelegationTokenResponseProto),
        MethodDescriptor("renewDelegationToken", 11, hadoop.common.RenewDelegationTokenRequestProto, hadoop.common.RenewDelegationTokenResponseProto),
        MethodDescriptor("cancelDelegationToken", 12, hadoop.common.CancelDelegationTokenRequestProto, hadoop.common.CancelDelegationTokenResponseProto),
        MethodDescriptor("moveApplicationAcrossQueues", 13, hadoop.yarn.MoveApplicationAcrossQueuesRequestProto, hadoop.yarn.MoveApplicationAcrossQueuesResponseProto),
        MethodDescriptor("getApplicationAttemptReport", 14, hadoop.yarn.GetApplicationAttemptReportRequestProto, hadoop.yarn.GetApplicationAttemptReportResponseProto),
        MethodDescriptor("getApplicationAttempts", 15, hadoop.yarn.GetApplicationAttemptsRequestProto, hadoop.yarn.GetApplicationAttemptsResponseProto),
        MethodDescriptor("getContainerReport", 16, hadoop.yarn.GetContainerReportRequestProto, hadoop.yarn.GetContainerReportResponseProto),
        MethodDescriptor("getContainers", 17, hadoop.yarn.GetContainersRequestProto, hadoop.yarn.GetContainersResponseProto),
        MethodDescriptor("submitReservation", 18, hadoop.yarn.ReservationSubmissionRequestProto, hadoop.yarn.ReservationSubmissionResponseProto),
        MethodDescriptor("updateReservation", 19, hadoop.yarn.ReservationUpdateRequestProto, hadoop.yarn.ReservationUpdateResponseProto),
        MethodDescriptor("deleteReservation", 20, hadoop.yarn.ReservationDeleteRequestProto, hadoop.yarn.ReservationDeleteResponseProto),
        MethodDescriptor("getNodeToLabels", 21, hadoop.yarn.GetNodesToLabelsRequestProto, hadoop.yarn.GetNodesToLabelsResponseProto),
        MethodDescriptor("getLabelsToNodes", 22, hadoop.yarn.GetLabelsToNodesRequestProto, hadoop.yarn.GetLabelsToNodesResponseProto),
        MethodDescriptor("getClusterNodeLabels", 23, hadoop.yarn.GetClusterNodeLabelsRequestProto, hadoop.yarn.GetClusterNodeLabelsResponseProto)
    ] # const _ApplicationClientProtocolService_methods
const _ApplicationClientProtocolService_desc = ServiceDescriptor("hadoop.yarn.ApplicationClientProtocolService", 1, _ApplicationClientProtocolService_methods)

ApplicationClientProtocolService(impl::Module) = ProtoService(_ApplicationClientProtocolService_desc, impl)

mutable struct ApplicationClientProtocolServiceStub <: AbstractProtoServiceStub{false}
    impl::ProtoServiceStub
    ApplicationClientProtocolServiceStub(channel::ProtoRpcChannel) = new(ProtoServiceStub(_ApplicationClientProtocolService_desc, channel))
end # mutable struct ApplicationClientProtocolServiceStub

mutable struct ApplicationClientProtocolServiceBlockingStub <: AbstractProtoServiceStub{true}
    impl::ProtoServiceBlockingStub
    ApplicationClientProtocolServiceBlockingStub(channel::ProtoRpcChannel) = new(ProtoServiceBlockingStub(_ApplicationClientProtocolService_desc, channel))
end # mutable struct ApplicationClientProtocolServiceBlockingStub

getNewApplication(stub::ApplicationClientProtocolServiceStub, controller::ProtoRpcController, inp::hadoop.yarn.GetNewApplicationRequestProto, done::Function) = call_method(stub.impl, _ApplicationClientProtocolService_methods[1], controller, inp, done)
getNewApplication(stub::ApplicationClientProtocolServiceBlockingStub, controller::ProtoRpcController, inp::hadoop.yarn.GetNewApplicationRequestProto) = call_method(stub.impl, _ApplicationClientProtocolService_methods[1], controller, inp)

getApplicationReport(stub::ApplicationClientProtocolServiceStub, controller::ProtoRpcController, inp::hadoop.yarn.GetApplicationReportRequestProto, done::Function) = call_method(stub.impl, _ApplicationClientProtocolService_methods[2], controller, inp, done)
getApplicationReport(stub::ApplicationClientProtocolServiceBlockingStub, controller::ProtoRpcController, inp::hadoop.yarn.GetApplicationReportRequestProto) = call_method(stub.impl, _ApplicationClientProtocolService_methods[2], controller, inp)

submitApplication(stub::ApplicationClientProtocolServiceStub, controller::ProtoRpcController, inp::hadoop.yarn.SubmitApplicationRequestProto, done::Function) = call_method(stub.impl, _ApplicationClientProtocolService_methods[3], controller, inp, done)
submitApplication(stub::ApplicationClientProtocolServiceBlockingStub, controller::ProtoRpcController, inp::hadoop.yarn.SubmitApplicationRequestProto) = call_method(stub.impl, _ApplicationClientProtocolService_methods[3], controller, inp)

forceKillApplication(stub::ApplicationClientProtocolServiceStub, controller::ProtoRpcController, inp::hadoop.yarn.KillApplicationRequestProto, done::Function) = call_method(stub.impl, _ApplicationClientProtocolService_methods[4], controller, inp, done)
forceKillApplication(stub::ApplicationClientProtocolServiceBlockingStub, controller::ProtoRpcController, inp::hadoop.yarn.KillApplicationRequestProto) = call_method(stub.impl, _ApplicationClientProtocolService_methods[4], controller, inp)

getClusterMetrics(stub::ApplicationClientProtocolServiceStub, controller::ProtoRpcController, inp::hadoop.yarn.GetClusterMetricsRequestProto, done::Function) = call_method(stub.impl, _ApplicationClientProtocolService_methods[5], controller, inp, done)
getClusterMetrics(stub::ApplicationClientProtocolServiceBlockingStub, controller::ProtoRpcController, inp::hadoop.yarn.GetClusterMetricsRequestProto) = call_method(stub.impl, _ApplicationClientProtocolService_methods[5], controller, inp)

getApplications(stub::ApplicationClientProtocolServiceStub, controller::ProtoRpcController, inp::hadoop.yarn.GetApplicationsRequestProto, done::Function) = call_method(stub.impl, _ApplicationClientProtocolService_methods[6], controller, inp, done)
getApplications(stub::ApplicationClientProtocolServiceBlockingStub, controller::ProtoRpcController, inp::hadoop.yarn.GetApplicationsRequestProto) = call_method(stub.impl, _ApplicationClientProtocolService_methods[6], controller, inp)

getClusterNodes(stub::ApplicationClientProtocolServiceStub, controller::ProtoRpcController, inp::hadoop.yarn.GetClusterNodesRequestProto, done::Function) = call_method(stub.impl, _ApplicationClientProtocolService_methods[7], controller, inp, done)
getClusterNodes(stub::ApplicationClientProtocolServiceBlockingStub, controller::ProtoRpcController, inp::hadoop.yarn.GetClusterNodesRequestProto) = call_method(stub.impl, _ApplicationClientProtocolService_methods[7], controller, inp)

getQueueInfo(stub::ApplicationClientProtocolServiceStub, controller::ProtoRpcController, inp::hadoop.yarn.GetQueueInfoRequestProto, done::Function) = call_method(stub.impl, _ApplicationClientProtocolService_methods[8], controller, inp, done)
getQueueInfo(stub::ApplicationClientProtocolServiceBlockingStub, controller::ProtoRpcController, inp::hadoop.yarn.GetQueueInfoRequestProto) = call_method(stub.impl, _ApplicationClientProtocolService_methods[8], controller, inp)

getQueueUserAcls(stub::ApplicationClientProtocolServiceStub, controller::ProtoRpcController, inp::hadoop.yarn.GetQueueUserAclsInfoRequestProto, done::Function) = call_method(stub.impl, _ApplicationClientProtocolService_methods[9], controller, inp, done)
getQueueUserAcls(stub::ApplicationClientProtocolServiceBlockingStub, controller::ProtoRpcController, inp::hadoop.yarn.GetQueueUserAclsInfoRequestProto) = call_method(stub.impl, _ApplicationClientProtocolService_methods[9], controller, inp)

getDelegationToken(stub::ApplicationClientProtocolServiceStub, controller::ProtoRpcController, inp::hadoop.common.GetDelegationTokenRequestProto, done::Function) = call_method(stub.impl, _ApplicationClientProtocolService_methods[10], controller, inp, done)
getDelegationToken(stub::ApplicationClientProtocolServiceBlockingStub, controller::ProtoRpcController, inp::hadoop.common.GetDelegationTokenRequestProto) = call_method(stub.impl, _ApplicationClientProtocolService_methods[10], controller, inp)

renewDelegationToken(stub::ApplicationClientProtocolServiceStub, controller::ProtoRpcController, inp::hadoop.common.RenewDelegationTokenRequestProto, done::Function) = call_method(stub.impl, _ApplicationClientProtocolService_methods[11], controller, inp, done)
renewDelegationToken(stub::ApplicationClientProtocolServiceBlockingStub, controller::ProtoRpcController, inp::hadoop.common.RenewDelegationTokenRequestProto) = call_method(stub.impl, _ApplicationClientProtocolService_methods[11], controller, inp)

cancelDelegationToken(stub::ApplicationClientProtocolServiceStub, controller::ProtoRpcController, inp::hadoop.common.CancelDelegationTokenRequestProto, done::Function) = call_method(stub.impl, _ApplicationClientProtocolService_methods[12], controller, inp, done)
cancelDelegationToken(stub::ApplicationClientProtocolServiceBlockingStub, controller::ProtoRpcController, inp::hadoop.common.CancelDelegationTokenRequestProto) = call_method(stub.impl, _ApplicationClientProtocolService_methods[12], controller, inp)

moveApplicationAcrossQueues(stub::ApplicationClientProtocolServiceStub, controller::ProtoRpcController, inp::hadoop.yarn.MoveApplicationAcrossQueuesRequestProto, done::Function) = call_method(stub.impl, _ApplicationClientProtocolService_methods[13], controller, inp, done)
moveApplicationAcrossQueues(stub::ApplicationClientProtocolServiceBlockingStub, controller::ProtoRpcController, inp::hadoop.yarn.MoveApplicationAcrossQueuesRequestProto) = call_method(stub.impl, _ApplicationClientProtocolService_methods[13], controller, inp)

getApplicationAttemptReport(stub::ApplicationClientProtocolServiceStub, controller::ProtoRpcController, inp::hadoop.yarn.GetApplicationAttemptReportRequestProto, done::Function) = call_method(stub.impl, _ApplicationClientProtocolService_methods[14], controller, inp, done)
getApplicationAttemptReport(stub::ApplicationClientProtocolServiceBlockingStub, controller::ProtoRpcController, inp::hadoop.yarn.GetApplicationAttemptReportRequestProto) = call_method(stub.impl, _ApplicationClientProtocolService_methods[14], controller, inp)

getApplicationAttempts(stub::ApplicationClientProtocolServiceStub, controller::ProtoRpcController, inp::hadoop.yarn.GetApplicationAttemptsRequestProto, done::Function) = call_method(stub.impl, _ApplicationClientProtocolService_methods[15], controller, inp, done)
getApplicationAttempts(stub::ApplicationClientProtocolServiceBlockingStub, controller::ProtoRpcController, inp::hadoop.yarn.GetApplicationAttemptsRequestProto) = call_method(stub.impl, _ApplicationClientProtocolService_methods[15], controller, inp)

getContainerReport(stub::ApplicationClientProtocolServiceStub, controller::ProtoRpcController, inp::hadoop.yarn.GetContainerReportRequestProto, done::Function) = call_method(stub.impl, _ApplicationClientProtocolService_methods[16], controller, inp, done)
getContainerReport(stub::ApplicationClientProtocolServiceBlockingStub, controller::ProtoRpcController, inp::hadoop.yarn.GetContainerReportRequestProto) = call_method(stub.impl, _ApplicationClientProtocolService_methods[16], controller, inp)

getContainers(stub::ApplicationClientProtocolServiceStub, controller::ProtoRpcController, inp::hadoop.yarn.GetContainersRequestProto, done::Function) = call_method(stub.impl, _ApplicationClientProtocolService_methods[17], controller, inp, done)
getContainers(stub::ApplicationClientProtocolServiceBlockingStub, controller::ProtoRpcController, inp::hadoop.yarn.GetContainersRequestProto) = call_method(stub.impl, _ApplicationClientProtocolService_methods[17], controller, inp)

submitReservation(stub::ApplicationClientProtocolServiceStub, controller::ProtoRpcController, inp::hadoop.yarn.ReservationSubmissionRequestProto, done::Function) = call_method(stub.impl, _ApplicationClientProtocolService_methods[18], controller, inp, done)
submitReservation(stub::ApplicationClientProtocolServiceBlockingStub, controller::ProtoRpcController, inp::hadoop.yarn.ReservationSubmissionRequestProto) = call_method(stub.impl, _ApplicationClientProtocolService_methods[18], controller, inp)

updateReservation(stub::ApplicationClientProtocolServiceStub, controller::ProtoRpcController, inp::hadoop.yarn.ReservationUpdateRequestProto, done::Function) = call_method(stub.impl, _ApplicationClientProtocolService_methods[19], controller, inp, done)
updateReservation(stub::ApplicationClientProtocolServiceBlockingStub, controller::ProtoRpcController, inp::hadoop.yarn.ReservationUpdateRequestProto) = call_method(stub.impl, _ApplicationClientProtocolService_methods[19], controller, inp)

deleteReservation(stub::ApplicationClientProtocolServiceStub, controller::ProtoRpcController, inp::hadoop.yarn.ReservationDeleteRequestProto, done::Function) = call_method(stub.impl, _ApplicationClientProtocolService_methods[20], controller, inp, done)
deleteReservation(stub::ApplicationClientProtocolServiceBlockingStub, controller::ProtoRpcController, inp::hadoop.yarn.ReservationDeleteRequestProto) = call_method(stub.impl, _ApplicationClientProtocolService_methods[20], controller, inp)

getNodeToLabels(stub::ApplicationClientProtocolServiceStub, controller::ProtoRpcController, inp::hadoop.yarn.GetNodesToLabelsRequestProto, done::Function) = call_method(stub.impl, _ApplicationClientProtocolService_methods[21], controller, inp, done)
getNodeToLabels(stub::ApplicationClientProtocolServiceBlockingStub, controller::ProtoRpcController, inp::hadoop.yarn.GetNodesToLabelsRequestProto) = call_method(stub.impl, _ApplicationClientProtocolService_methods[21], controller, inp)

getLabelsToNodes(stub::ApplicationClientProtocolServiceStub, controller::ProtoRpcController, inp::hadoop.yarn.GetLabelsToNodesRequestProto, done::Function) = call_method(stub.impl, _ApplicationClientProtocolService_methods[22], controller, inp, done)
getLabelsToNodes(stub::ApplicationClientProtocolServiceBlockingStub, controller::ProtoRpcController, inp::hadoop.yarn.GetLabelsToNodesRequestProto) = call_method(stub.impl, _ApplicationClientProtocolService_methods[22], controller, inp)

getClusterNodeLabels(stub::ApplicationClientProtocolServiceStub, controller::ProtoRpcController, inp::hadoop.yarn.GetClusterNodeLabelsRequestProto, done::Function) = call_method(stub.impl, _ApplicationClientProtocolService_methods[23], controller, inp, done)
getClusterNodeLabels(stub::ApplicationClientProtocolServiceBlockingStub, controller::ProtoRpcController, inp::hadoop.yarn.GetClusterNodeLabelsRequestProto) = call_method(stub.impl, _ApplicationClientProtocolService_methods[23], controller, inp)

export ApplicationClientProtocolService, ApplicationClientProtocolServiceStub, ApplicationClientProtocolServiceBlockingStub, getNewApplication, getApplicationReport, submitApplication, forceKillApplication, getClusterMetrics, getApplications, getClusterNodes, getQueueInfo, getQueueUserAcls, getDelegationToken, renewDelegationToken, cancelDelegationToken, moveApplicationAcrossQueues, getApplicationAttemptReport, getApplicationAttempts, getContainerReport, getContainers, submitReservation, updateReservation, deleteReservation, getNodeToLabels, getLabelsToNodes, getClusterNodeLabels


text:
#=
  Package: FishEBM
  File: spawn.jl
  Devin Rose
  Functions for simulating one week of a spawn season
  Created: June 2016
=#

"""
  Description: Generates a brood size for each spawning location based on number
    of adults in the spawning location, and their age specific fecundity. Brood
    size is affected by a compensation factor based on the total adult
    population in the environment.

  Returns: Operates directly on agent_db

  Last Update: September 2016
"""
function spawn!(agent_db::Vector{EnviroAgent}, adult_a::AdultAssumptions,
  age_assumpt::AgentAssumptions, enviro_a::EnvironmentAssumptions, week::Int64,
  carryingcapacity::Float64, sdf::DataFrame)

  adult_pop = getStagePopulation(4, week, agent_db, age_assumpt)
  totalPop = 0
  for stage = 1:4
    totalPop += getStagePopulation(stage, week, agent_db, age_assumpt)
  end

  if isnan(adult_a.fecunditycompensation)
    compensation_factor_a = 1
  else
    compensation_factor_a = 2*(1-cdf(Normal(carryingcapacity, carryingcapacity/adult_a.fecunditycompensation), totalPop))
  end

  @assert(0.0 <= compensation_factor_a <= 2.0, "Population regulation has failed in spawn (part a), compensation expected:0.0-2.0, actual:$compensation_factor_a")

  if isnan(adult_a.maturitycompensation)
    compensation_factor_b = 1
  else
    compensation_factor_b = 2*(1-cdf(Normal(carryingcapacity, carryingcapacity/adult_a.maturitycompensation), totalPop))
  end

  @assert(0.0 <= compensation_factor_b <= 2.0, "Population regulation has failed in spawn (part b), compensation expected:0.0-2.0, actual:$compensation_factor_a")

  newClass = length(agent_db[1].weekNum)

  #Check if currently in same spawning season as most recent class. If so, continue adding brood to that class.
  #If not, add a new class to each agent and begin adding brood to the new class.
  if week - agent_db[1].weekNum[newClass] > 12
    for i = 1:length(agent_db)
      push!((agent_db[i]).alive, 0)
      push!((agent_db[i]).weekNum, week)
    end
    newClass += 1
  end

  brood_size = fill(0, size(adult_a.broodsize))

  for i = 1:length(enviro_a.spawningHash)
    if isEmpty(agent_db[enviro_a.spawningHash[i]]) == false
      for age = 2:8
        ageSpecificPop = getAgeSpecificPop(age, week, agent_db[enviro_a.spawningHash[i]].alive, agent_db[enviro_a.spawningHash[i]].weekNum, age_assumpt)
        if ageSpecificPop > 1
          #ageSpecificPop/2 implies a 50% male/female ratio, this will be number of females able to lay eggs
          numSpawningAdults = rand(1:ageSpecificPop/2)
          for j = 1:numSpawningAdults
            brood = rand(Poisson(compensation_factor_a*adult_a.broodsize[age - 1]), rand(Binomial(ageSpecificPop, cdf(Binomial(length(adult_a.broodsize)+2, min(1, compensation_factor_b*adult_a.halfmature/(length(adult_a.broodsize)+2))), age)*0.5)))
            for k = 1:length(brood)
              agent_db[enviro_a.spawningHash[i]].alive[newClass] += brood[k]
              brood_size[age - 1] += brood[k]
            end #for k=1:brood
          end #for j = 1:numSpawningAdults
        end #if ageSpecificPop
      end #for ages
    end #if isEmpty
  end #for i=1:length spawningHash

  #Push to weekly spawn number data frame
  push!(sdf, (vcat(week, brood_size..., sum(brood_size))))

end


text:
# Following code can be plugged into cellbox.jl

begin
    function  loss_all(p)
        loss_epoch = zeros(Float32, n_exp);
        for i_exp in 1:n_exp
            loss_epoch[i_exp] = loss_neuralode(p, i_exp)
        end
        loss_total = mean(loss_epoch);
    end
    p_grad = Zygote.gradient(x -> loss_neuralode(x), p)[1]

    backend(:pyplot)

    xs = vcat(["α"], [string("W i=", i) for i = 1:ns])
    ys = [string("W j=", i) for i = 1:ns]
    plt = heatmap(xs, ys, abs.(p_grad) ./ maximum(abs.(p_grad)), aspect_ratio = 1);
    title!(plt, "Scaled sensitivity map p(α | W)");

    plt_p = heatmap(xs, ys, abs.(p), aspect_ratio = 1);
    title!(plt_p, "map p");

    plt_p_gold = heatmap(xs, ys, abs.(p_gold), aspect_ratio = 1);
    title!(plt_p_gold, "map p_gold");

    plt_p_diff = heatmap(xs, ys, abs.(p_gold - p), aspect_ratio = 1);
    title!(plt_p_diff, "map p_gold - p");

    plt_all = plot(plt, plt_p, plt_p_gold, plt_p_diff, layout=(2, 2), size=(1600, 1600));
    png(plt_all, string(fig_path, "/sens"))

    backend(:gr)
end


text:
function getmetaboliteintervals(tag::String,
    molecule_names,
    ΩS_ppm::Vector{Vector{Vector{T}}},
    Δsys_border::Vector{Vector{T}};
    min_dist::T = 0.15) where T <: Real

    n_select = findfirst(xx->(xx==tag), molecule_names)
    C0 = ΩS_ppm[n_select]
    N_groups = length(C0)

    Δcompound_border = Δsys_border[n_select]
    @assert N_groups == length(Δcompound_border)

    st_c = Vector{Vector{T}}(undef, N_groups)
    fin_c = Vector{Vector{T}}(undef, N_groups)

    for i = 1:N_groups
        st_c[i], fin_c[i] = segmentmoleculefreqs2(C0[i],
            Δcompound_border[i], min_dist = min_dist)
    end

    return MetaboliteRegionType(st_c, fin_c)
end

# A region needs to contain all interval C0[l] +/- r such that
#   all intervals are within min_dist from each other.
# r and min_dst are in ppm.
function segmentmoleculefreqs2(C0::Vector{T}, border::T;
    min_dist::T = 0.15) where T <: Real

    sort_inds = sortperm(C0)
    C = C0[sort_inds]

    st_c = Vector{T}(undef, 0)
    fin_c = Vector{T}(undef, 0)

    push!(st_c, C[1])

    for l = 2:length(C)
        left = C[l-1] + border
        right = C[l] -border

        if (right-left) > min_dist
            push!(fin_c, C[l-1])
            push!(st_c, C[l])
        end
    end
    push!(fin_c, C[end])

    return st_c, fin_c
end

function constructgraph(sts::Vector{T}, fins::Vector{T}) where T

    # construct graph.
    N_intervals = length(sts)
    g = LightGraphs.SimpleGraph(N_intervals)

    for i = 1:N_intervals
        flags = isthereoverlap(sts[i], fins[i], sts, fins)

        for j = 1:length(flags)
            if flags[j]
                #add edge i to j, j to i.
                LightGraphs.add_edge!(g, i, j)
                LightGraphs.add_edge!(g, j, i)
            end
        end
    end

    return g
end

function isthereoverlap(st::T, fin::T, sts, fins) where T

    flags = collect( isoverlapped(st, fin, sts[i], fins[i]) for i = 1:length(fins))

    return flags
end

# assumes st_a < fin_a, st_b < fin_b.
function isoverlapped(st_a::T, fin_a::T, st_b::T, fin_b::T) where T

    if st_a < st_b && st_b < fin_a < fin_b
        return true

    elseif st_b < st_a && st_a < fin_b < fin_a
        return true

    elseif st_b < st_a < fin_b
        return true

    elseif st_a < st_b < fin_a
        return true

    end

    return false
end


text:
@testset "E522" begin
    msgs = lintstr("""
    s = I
    println(s[1])
    """)

    @test messageset(msgs) == Set([:E522])
    @test msgs[1].variable == "s[1]"
    @test occursin("indexing UniformScaling", msgs[1].message)
    @test occursin("with index types $Int is not supported", msgs[1].message)

    # dicts
    @test messageset(lintstr("""
    s = keys(Dict(1 => 2))
    println(s["foo"])
    """)) == Set([:E522])
    @test messageset(lintstr("""
    s = Dict(1 => 2)
    println(s["foo"])
    """)) == Set([:E522])
    @test isempty(messageset(lintstr("""
    s = Dict(1 => 2)
    println(s[1])
    """)))

    msgs = lintstr("""
    function f()
        a = 1
        d = Dict{Symbol,Int}(:a=>1, :b=>2)
        x = d[a]
        return x
    end
    """)
    @test messageset(msgs) == Set([:E522, :E539])
    @test msgs[1].variable == "d[a]"
    @test occursin("indexing Dict", msgs[1].message)
    @test occursin("Int", msgs[1].message)

    # strings
    msgs = lintstr("""
    function f()
        b = repeat(" ", 10)
        b[:start]
    end
    """)
    @test messageset(msgs) == Set([:E522])
    @test msgs[1].variable == "b[:start]"
    @test occursin("indexing String", msgs[1].message)
    @test occursin("with index types Symbol is not supported", msgs[1].message)

    # zero-dimensional indexing
    msgs = lintstr("""
    d = Dict()
    x = d[]
    """)
    @test messageset(msgs) == Set([:E522, :E539])
    @test occursin("indexing Dict", msgs[1].message)
    @test occursin("with no indices is not supported", msgs[1].message)

    msgs = lintstr("""
    a = ""
    a[]
    """)
    @test messageset(msgs) == Set([:E522])
    @test occursin("indexing String with no indices", msgs[1].message)

    # issue 196
    @test messageset(lintstr("""
    s = Dict(:b => 2)
    println(keys(s)[1])
    """)) == Set([:E522])
end


text:
# This file contains algorithms needed for the solution of
# the parametric Risch differential equation from chapter 7 of the book
#
#    Manuel Bronstein, Symbolic Integration I: Transcendental Functions,
#    2nd ed., Springer 2005.
#

"""
    ParamRdeNormalDenominator(f, gs, D) -> (a, b, Gs, h)

Normal part of the denominator.

Given a field `k`, a derivation `D` on `k[t]`, `f` in `k(t)` which is weakly normalized with respect to `t`,
and `gs=[g₁,...,gₘ]` with `gᵢ` in `k(t)`, return  `a`, `h` in `k[t]`, `b` in `k⟨t⟩` and  `Gs=[G₁,...,Gₘ]`
with `G_i` in `k(t)` such that for any solution `c₁`,...,`cₘ` in `Const(k)` and `y` in `k(t)`
of `D(y)+f*y=∑ᵢcᵢ*gᵢ`, `q=y*h` in `k⟨t⟩` satisfies `a*D(q)+b*q=∑ᵢcᵢ*Gᵢ`.

(Here, `k⟨t⟩` denotes the elements of `k(t)` which are reduced w.r.t. `D`.)

See [Bronstein's book](https://link.springer.com/book/10.1007/b138171), Section 7.1, p. 219.
"""
function ParamRdeNormalDenominator(f::F, gs::Vector{F}, D::Derivation) where
    {T<:FieldElement, P<:PolyElem{T}, F<:FracElem{P}}
    iscompatible(f, D) && all(iscompatible(g, D) for g in gs) ||
        error("rational functions f and g_i must be in the domain of derivation D")
    # Note: f must be weakly normalized which we do not check. It is recommended
    # to use this function only for rational functions f which were computed by WeakNormalizer.
    (dn, ds) = SplitFactor(denominator(f), D)
    (en, es) = SplitFactor(lcm([denominator(g) for g in gs]), D)
    p = gcd(dn, en)
    h = divexact(gcd(en, derivative(en)), gcd(p, derivative(p)))
    dnh2 = dn*h^2
    dn*h, dn*h*f-dn*D(h), [dnh2*g for g in gs], h
end

"""
    ParamRdeSpecialDenomExp(a, b, gs, D) -> (A, B, Gs, h)

Special part of the denominator - hyperexponential case.

Given a field `k`, a derivation `D` on `k[t]` with `D(t)/t` in `k`, `a` in `k[t]`, `b` in `k⟨t⟩`, `gs=[g₁,...,gₘ]`
with `gᵢ` in `k(t)` , `a≠0` and `gcd(a,t)=1`, return `A`, `B`, `h` in `k[t]`,  `Gs=[G₁,...,Gₘ]` with `Gᵢ` in `k(t)`
such that for any solution `c₁`,...,`cₘ` in `Const(k)` and `q` in `k⟨t⟩`
of `a*D(q)+b*q=∑ᵢcᵢ*gᵢ`, `r=q*h` in `k[t]` satisfies `A*D(r)+B*r=∑ᵢcᵢ*Gᵢ`.

(Here, `k⟨t⟩` denotes the elements of `k(t)` which are reduced w.r.t. `D`.)

See [Bronstein's book](https://link.springer.com/book/10.1007/b138171), Section 7.1, p. 221.
"""
function ParamRdeSpecialDenomExp(a::P, b::F, gs::Vector{F}, D::Derivation) where
    {T<:FieldElement, P<:PolyElem{T}, F<:FracElem{P}}
    ishyperexponential(D) ||
        error("monomial of derivation D must be hyperexponential")
    iscompatible(a, D) && iscompatible(b, D) && all([iscompatible(g, D) for g in gs]) ||
        error("polynomial a and rational functions b and g_i must be in the domain of derivation D")
    !iszero(a) || error("a must be != 0")
    isreduced(b, D) ||
        error("rational function b must be reduced with respect to derivation D")
    t = gen(parent(a))
    degree(gcd(t, a))==0 || error("gcd(a, t) must be == 1")
    p = t
    nb = iszero(b) ? typemax(Int) : valuation(b, p)
    nc = minimum([iszero(g) ? typemax(Int) : valuation(g, p) for g in gs])
    n = min(0, nc - min(0, nb))
    if nb==0
        α = constant_coefficient(Remainder(-b//a, p))
        w = coeff(MonomialDerivative(D), 1)
        n0, s, z, ρ = ParametricLogarithmicDerivative(α, w, BaseDerivation(D))
        if ρ<0
            error("ParametricLogarithmicDerivative failed")
        end
        if  ρ>0 && n0==1 && !iszero(z)
            n = min(n, s)
        end
    end
    N = max(0, -nb)
    p_power_N = p^N
    p_power_N_minus_n = p^(N-n)
    b1 = (b+n*a*divexact(D(p), p))*p_power_N
    @assert isone(denominator(b1))
    a*p_power_N, numerator(b1), [g*p_power_N_minus_n for g in gs], p^(-n)
end

"""
    ParamRdeSpecialDenomTan(a, b, gs, D) -> (A, B, Gs, h)

Special part of the denominator - hypertangent case.

Given a field `k` not containing `√-1`, a derivation `D` on `k[t]` with `D(t)/(t^2+1)` in `k`,
`a≠0` in `k[t]` with `gcd(a,t^2+1)=1`, `b` in `k⟨t⟩`, and `gs=[g₁,...,gₘ]` with `gᵢ` in `k(t)`,
return  `A`, `B`, `h` in `k[t]` and `Gs=[G₁,...,Gₘ]` with `Gᵢ` in `k(t)`
such that for any solution `c₁`,...,`cₘ` in `Const(k)` and `q` in `k⟨t⟩`
of `a*D(q)+b*q=∑ᵢcᵢ*gᵢ`, `r=q*h` in `k[t]` satisfies `A*D(r)+B*r=∑ᵢcᵢ*Gᵢ`.

(Here, `k⟨t⟩` denotes the elements of `k(t)` which are reduced w.r.t. `D`.)

See [Bronstein's book](https://link.springer.com/book/10.1007/b138171), Section 7.1, p. 222.
"""
function ParamRdeSpecialDenomTan(a::P, b::F, gs::Vector{F}, D::Derivation) where
    {T<:FieldElement, P<:PolyElem{T}, F<:FracElem{P}}
    !iszero(a) || error("a must be != 0")
    ishypertangent(D) ||
        error("monomial of derivation D must be hypertangent")
    iscompatible(a, D) && iscompatible(b, D) && all([iscompatible(g, D) for g in gs]) ||
        error("polynomial a and rational functions b and g_i must be in the domain of derivation D")
    isreduced(b, D) ||
        error("rational function b must be reduced with respect to derivation D")
    t = gen(parent(a))
    p = t^2+1
    degree(gcd(a, p))==0 || error("gcd(a, t^2+1) must be == 1")
    nb = iszero(b) ? typemax(Int) : valuation(b, p)
    nc = minimum([iszero(g) ? typemax(Int) : valuation(g, p) for g in gs])
    n = min(0, nc - min(0, nb))
    if nb==0
        αI_plus_β = Remainder(-b//a, p)
        α = coeff(αI_plus_β, 1)
        β = coeff(αI_plus_β, 0)
        η = divexact(MonomialDerivative(D), p)
        v, ρ = InFieldDerivative(2*β, D0)
        D0 = BaseDerivation(D)
        if ρ>0 && !iszero(v)
            _, I, D0I = Complexify(parent(η), D0)
            n0, m, v, ρ = ParametricLogarithmicDerivative(α*I+β, 2*η*I, D0I)
            if  ρ>0 && n0==1 && !iszero(v)
                n = min(n, m)
            end
        end
    end
    N = max(0, -nb)
    p_power_N = p^N
    b1 = (b+n*a*divexact(D(p), p))*p_power_N
    @assert isone(denominator(b1))
    a*p_power_N, numerator(b1), [g*p_power_N_minus_n for g in gs], p^(-n)
end

"""
    ParamRdeSpecialDenominator(a, b, gs, D) -> (A, B, Gs, h)

Special part of the denominator.

Given a field `k`, a derivation `D` on `k[t]`, `a≠0` in `k[t]`, `b` in `k⟨t⟩`, `gs=[g₁,...,gₘ]` with
`g₁`,...,`gₘ` in `k(t)`, return  `A`, `B`, `h` in `k[t]` and `Gs=[G₁,...,Gₘ]` with `Gᵢ` in `k(t)`
such that for any solution `c₁`,...,`cₘ` in `Const(k)` and `q` in `k⟨t⟩`
of `a*D(q)+b*q=∑ᵢcᵢ*gᵢ`, `r=q*h` in `k[t]` satisfies `A*D(r)+B*r=∑ᵢcᵢ*Gᵢ`.

(Here, `k⟨t⟩` denotes the elements of `k(t)` which are reduced w.r.t. `D`.)

See [Bronstein's book](https://link.springer.com/book/10.1007/b138171), Section 7.1, p. 221.
"""
function ParamRdeSpecialDenominator(a::P, b::F, gs::Vector{F}, D::Derivation) where
    {T<:FieldElement, P<:PolyElem{T}, F<:FracElem{P}}
    iscompatible(a, D) && iscompatible(b, D) && all([iscompatible(g, D) for g in gs]) ||
        error("polynomial a and rational functions b and g_i must be in the domain of derivation D")
    !iszero(a) || error("a must be != 0")
    isreduced(b, D) ||
        error("rational function b must be reduced with respect to derivation D")
    if isprimitive(D)
        @assert isone(denominator(b))
        return a, numerator(b), gs, one(b)
    elseif ishyperexponential(D)
        t = gen(parent(a))
        d = gcd(a, t)
        a = divexact(a, d)
        b = b//d
        gs = [g//d for g in gs]
        return ParamRdeSpecialDenomExp(a, b, gs, D)
    elseif ishypertangent(D)
        t = gen(parent(a))
        p = t^2 + 1
        d = gcd(a, p)
        a = divexact(a, d)
        b = b//d
        gs = [g//d for g in gs]
        return ParamRdeSpecialDenomTan(a, b, gs, D)
    else
        H = MonomialDerivative(D)
        throw(NotImplementedError("ParamRdeSpecialDenominator: monomial derivative $H\n@ $(@__FILE__):$(@__LINE__)"))
    end
end

"""
    LinearConstraints(a, b, gs, D) -> (qs, M)

Generate linear constraints on the constants.

Given a field `k`, a derivation `D` on `k[t]`, `a`, `b` in `k[t]` and `gs=[g₁,...,gₘ]` with
`gᵢ` in `k(t)`, return `qs=[q₁,...,qₘ]` with `qᵢ` in `k[t]`
and a matrix `M` with entries in `k` sucht that for any solution
`c₁`,...,`cₘ` in `Const(k)` and `p` in `k[t]` of  `a*D(p)+b*p=∑ᵢcᵢ*gᵢ`,
`x=[c₁,...,cₘ]` is a solution of `M*x=0`, and `p` and the `cᵢ` satisfy
 `a*D(p)+b*p=∑ᵢcᵢ*qᵢ`.

See [Bronstein's book](https://link.springer.com/book/10.1007/b138171), Section 7.1, p. 223.
"""
function LinearConstraints(a::P, b::P, gs::Vector{F}, D::Derivation) where
    {T<:FieldElement, P<:PolyElem{T}, F<:FracElem{P}}
    iscompatible(a, D) && iscompatible(b, D) && all([iscompatible(g, D) for g in gs]) ||
        error("polynomials a and b and rational functions g_i must be in the domain of derivation D")
    d = lcm([denominator(g) for g in gs])
    qrs = [divrem(numerator(d*g), d) for g in gs]
    qs = [q for (q, r) in qrs]
    rs = [r for (q, r) in qrs]
    if all([iszero(r) for r in rs])
        N = -1
    else
        N = maximum([degree(r) for r in rs])
    end
    M = [coeff(r, i) for i=0:N, r in rs]
    return qs, M
end

RowEchelon(A::Matrix{T}; cut_zero_rows::Bool=true) where T<:FieldElement = RowEchelon(A, T[], cut_zero_rows=cut_zero_rows)

function RowEchelon(A::Matrix{T}, u::Vector{T}; cut_zero_rows::Bool=true) where T<:FieldElement
    # based on https://github.com/blegat/RowEchelon.jl/blob/master/src/RowEchelon.jl
    # Note: input arguments A, u are changed on output
    nr, nc = size(A)
    uzero = length(u)==0
    i = 1
    j = 1
    while i <= nr && j <= nc
        p = findfirst(x->!iszero(x), A[i:nr,j])
        if p===nothing
            j += 1
        else
            p = p + i - 1
            for k=j:nc
                A[i, k], A[p, k] = A[p, k], A[i, k]
            end
            if !uzero
                u[i], u[p] = u[p], u[i]
            end
            d = A[i,j]
            for k = j:nc
                A[i,k] //= d
            end
            if !uzero
                u[i] //= d
            end
            for k = 1:nr
                if k != i
                    d = A[k,j]
                    for l = j:nc
                        A[k,l] -= d*A[i,l]
                    end
                    if !uzero
                        u[k] -= d*u[i]
                    end
                end
            end
            i += 1
            j += 1
        end
    end
    if cut_zero_rows
        A = A[1:i-1,:]
    end
    if uzero
        return A
    else
        A, u
    end
end

"""
    ConstantSystem(A[, u], D) -> (B[, v])

Generate a system for the constant solutions

Given a field `K`, a derivation `D` on `K` with constant field `C`, a matrix `A`
with entries in `K` and optionally a vector `u` with coefficients in `K`,
return a matrix `B` with entries in `C` and, in the case that `u` is provided,
a vector `v` with coefficients in `K`. If no `u` is provided, or all entries
of the returned `v` are actually elements of `C`, then the solutions in `C`
of `A*x=u` (resp. `A*x=0` if no `u` is provided) are exactly all the solutions
of `B*x=v` (resp. `B*x=0`). Otherwise, if `v` has a nonconstant coefficient,
`A*x=u` has no constant solution.

See [Bronstein's book](https://link.springer.com/book/10.1007/b138171), Section 7.1, p. 225.
"""
ConstantSystem(A::Matrix{T}, D::Derivation) where T<:FieldElement =
    ConstantSystem(A, T[], D)

function ConstantSystem(A::Matrix{T}, u::Vector{T}, D::Derivation) where T<:FieldElement
    # check compatibility
    # check dimensions
    # Note: input arguments A, u are changed on output
    uzero = length(u)==0
    if uzero
        A = RowEchelon(A)
    else
        A, u = RowEchelon(A, u)
    end
    m, n = size(A)
    j = 1
    while j<=n
        i = findfirst(x->!isconstant(x, D), A[:,j])
        if i===nothing
            j += 1
            continue
        end
        # enlarge A by one row and u by one entry
        A = vcat(A, [zero(A[1,1]) for l=1:1, k=1:n])
        if !uzero
            push!(u, zero(A[1,1]))
        end
        Daij = D(A[i,j])
        for k=1:n
            A[m+1, k] =  D(A[i, k])//Daij
        end
        if !uzero
            u[m+1] = D(u[i])//Daij
        end
        for s=1:m
            asj = A[s,j]
            for k=1:n
                A[s, k] -= asj*A[m+1, k]
            end
            if !uzero
                u[s] -= asj*u[m+1]
            end
        end
        j += 1
        m += 1
    end
    # return only nonzero rows
    nz = [i for i=1:size(A,1) if !(all([iszero(A[i,j]) for j=1:size(A,2)]) && (uzero || iszero(u[i])))]
    CT = typeof(zero(constant_field(D))) # ensure B has eltype of constant field, even if it is empty
    B = CT[constantize(A[i,j], D) for i in nz, j=1:n]
    if uzero
        return B
    else
        B, u[nz]
    end
end

"""
    ParamRdeBoundDegreePrim(a, b, qs, D) -> n

Bound on polynomial solutions - primitive case.

Given a field `k`, a derivation `D` on `k[t]` with `D(t)` in `k`, `a`, `b` in `k[t]` with `a≠0`,
and `qs=[q₁,...,qₘ]` with `qᵢ` in `k[t]`, return integer `n`
such that `degree(q)≤n` for any solution  `c₁`,...,`cₘ` in `Const(k)`
and `q` in `k[t]` of  `a*D(q)+b*q=∑ᵢcᵢ*qᵢ`.

See [Bronstein's book](https://link.springer.com/book/10.1007/b138171), Section 7.1, p. 228.
"""
function ParamRdeBoundDegreePrim(a::P, b::P, qs::Vector{P}, D::Derivation) where
    {T<:FieldElement, P<:PolyElem{T}}
    isprimitive(D) ||
        error("monomial of derivation D must be primitive")
    iscompatible(a, D) && iscompatible(b, D) && all([iscompatible(q, D) for q in qs]) ||
        error("polynomial a and rational functions b and q_i must be in the domain of derivation D")
    !iszero(a) || error("polynomial a must be nonzero")
    da = degree(a)
    db = degree(b)
    dc = maximum([degree(q) for q in qs])
    if db>da
        n = max(0, dc-db)
    else
        n = max(0, dc-da+1)
    end
    D0 = BaseDerivation(D)
    if db==da-1
        α = -leading_coefficient(b)//leading_coefficient(a)
        z, s0, ρ = LimitedIntegrate(α, leading_coefficient(D), D0)
        if ρ>0 && isrational(s0)
            s = rationalize_over_Int(s0)
            if denominator(s)==1
                n = max(n, numerator(s))
            end
        end
    end
    if db==da
        α = -leading_coefficient(b)//leading_coefficient(a)
        z, ρ = InFieldDerivative(α, D0)
        if ρ>0 && !iszero(z)
            β = -leading_coefficient(a*D0(z)+b*z)//(z*leading_coefficient(a))
            w, s0, ρ = LimitedIntegrate(β, leading_coefficient(D), D0)
            if ρ>0 && isrational(s0)
                m = rationalize_over_Int(s0)
                if denominator(s)==1
                    n = max(n, numerator(s))
                end
            end
        end
    end
    n
end

"""
    ParamRdeBoundDegreeBase(a, b, qs) -> n

Bound on polynomial solutions - base case.

Given a field `k`, `a`, `b` in `k[t]` with `a≠0`, and `qs=[q₁,...,qₘ]` with
`qᵢ` in `k[t]`,  return integer `n`
such that `degree(q)≤n` for any solution  `c₁`,...,`cₘ` in `Const(k)`
and `q` in `k[t]` of `a*(d/dt)(q)+b*q=∑ᵢcᵢ*qᵢ`.

See [Bronstein's book](https://link.springer.com/book/10.1007/b138171), Section 7.1, p. 228.
"""
function ParamRdeBoundDegreeBase(a::P, b::P, qs::Vector{P}) where
    {T<:FieldElement, P<:PolyElem{T}}
    !iszero(a) || error("polynomial a must be nonzero")
    da = degree(a)
    db = degree(b)
    dc = maximum([degree(q) for q in qs])
    n = max(0, dc - max(db, da - 1))
    if db==da-1
        s0 = -leading_coefficient(b)//leading_coefficient(a)
        if isrational(s0)
            s = rationalize_over_Int(s0)
            if isone(denominator(s))
                n = max(0, numerator(s), dc - db)
            end
        end
    end
    return n
end

"""
    ParamRdeBoundDegreeExp(a, b, qs, D) -> n

Bound on polynomial solutions - hyperexponential case.

Given a field `k`, a derivation `D` on `k[t]` with `D(t)/t` in `k`, `a`, `b` in `k[t]`,
with `a≠0`, and `qs=[q₁,...,qₘ]` with `qᵢ` in `k[t]` return integer `n`
such that `degree(q)≤n` for any solution  `c₁`,...,`cₘ` in `Const(k)`
and `q` in `k[t]` of  `a*D(q)+b*q=∑ᵢcᵢ*qᵢ`.

See [Bronstein's book](https://link.springer.com/book/10.1007/b138171), Section 7.1, p. 229.
"""
function ParamRdeBoundDegreeExp(a::P, b::P, qs::Vector{P}, D::Derivation) where
    {T<:FieldElement, P<:PolyElem{T}}
    ishyperexponential(D) ||
        error("monomial of derivation D must be hyperexponential")
    iscompatible(a, D) && iscompatible(b, D) && all([iscompatible(q, D) for q in qs]) ||
        error("polynomial a and rational functions b and q_i must be in the domain of derivation D")
    !iszero(a) || error("polynomial a must be nonzero")
    da = degree(a)
    db = degree(b)
    dc = maximum([degree(q) for q in qs])
    n = max(0, dc - max(db, da))
    if da==db
        α = -leading_coefficient(b)//leading_coefficient(a)
        w = coeff(MonomialDerivative(D), 1)
        n0, s, z, ρ = ParametricLogarithmicDerivative(α, w, BaseDerivation(D))
        if ρ<0
            error("ParametricLogarithmicDerivative failed")
        end
        if  ρ>0 && n0==1 && !iszero(z)
            n = max(n, s)
        end
    end
    n
end

"""
    ParamRdeBoundDegreeNonLinear(a, b, qs, D) -> n

Bound on polynomial solutions - nonlinear case.

Given a  field `k`, a derivation `D` on `k[t]` with `degree(D(t))≥2` , `a`, `b` in `k[t]` with `a≠0`,
and `qs=[q₁,...,qₘ]` with `qᵢ` in `k[t]`, return integer `n`
such that `degree(q)≤n` for any solution `c₁`,...,`cₘ` in `Const(k)`
and `q` in `k[t]` of `a*D(q)+b*q=∑ᵢcᵢ*qᵢ`.

See [Bronstein's book](https://link.springer.com/book/10.1007/b138171), Section 7.1, p. 230.
"""
function ParamRdeBoundDegreeNonLinear(a::P, b::P, qs::Vector{P}, D::Derivation) where
    {T<:FieldElement, P<:PolyElem{T}}
    isnonlinear(D) ||
        error("monomial of derivation D must be nonlinear")
    iscompatible(a, D) && iscompatible(b, D) && all([iscompatible(q, D) for q in qs]) ||
        error("polynomial a and rational functions b and q_i must be in the domain of derivation D")
    !iszero(a) || error("polynomial a must be nonzero")
    da = degree(a)
    db = degree(b)
    dc = maximum([degree(q) for q in qs])
    δ = degree(D)
    λ = leading_coefficient(D)
    n = max(0, dc - max(da + δ - 1, db))
    if db==da+δ-1
        s0 = -leading_coefficient(b)/(λ*leading_coefficient(a))
        if isrational(s0)
            s = rationalize_over_Int(s0)
            if isone(denominator(s))
                n = max(0, numerator(s), dc - db)
            end
        end
    end
    n
end

"""
    ParamRdeBoundDegree(a, b, qs, D) -> n

Bound on polynomial solutions.

Given a  field `k`, a derivation `D` on `k[t]`, `a`, `b` in `k[t]` with `a≠0`,
 and `qs=[q₁,...,qₘ]` with `qᵢ` in `k[t]`, return integer `n`
such that `degree(q)≤n` for any solution `c₁`,...,`cₘ` in `Const(k)`
and `q` in `k[t]` of  `a*D(q)+b*q=∑ᵢcᵢ*qᵢ`.

See [Bronstein's book](https://link.springer.com/book/10.1007/b138171), Section 7.1, p. 227.
"""
function ParamRdeBoundDegree(a::P, b::P, qs::Vector{P}, D::Derivation) where
    {T<:FieldElement, P<:PolyElem{T}}
    iscompatible(a, D) && iscompatible(b, D) && all([iscompatible(q, D) for q in qs]) ||
        error("polynomial a and rational functions b and q_i must be in the domain of derivation D")
    !iszero(a) || error("polynomial a must be nonzero")
    if isbasic(D)
        return ParamRdeBoundDegreeBase(a, b, qs)
    elseif isprimitive(D)
        return ParamRdeBoundDegreePrim(a, b, qs, D)
    elseif ishyperexponential(D)
        return ParamRdeBoundDegreeExp(a, b, qs, D)
    elseif isnonlinear(D)
        return ParamRdeBoundDegreeNonLinear(a, b, qs, D)
    else
        H = MonomialDerivative(D)
        throw(NotImplementedError("ParamRdeBoundDegree: monomial derivative $H\n@ $(@__FILE__):$(@__LINE__)"))
    end
end

"""
    ParSPDE(a, b, qs, D, n) -> (A, B, Qs, rs, N)

Parametric Special Polynomial Differential Equation algorithm.

Given a field `k`, a derivation `D` on `k[t]`, `a`, `b` in `k[t]` with `degree(a)>0` and `gcd(a,b)=1`,
`qs=[q₁,...,qₘ]` with `qᵢ` in `k[t]`, and an integer `n` , return `A`, `B` in `k[t]`,
`Qs=[Q₁,...,Qₘ]` with `Qᵢ` in `k[t]`, `rs=[r₁,...,rₘ]` with `rᵢ` in `k[t]`,
and an integer `N` such that for any solution `c₁`,...,`cₘ` in `Const(k)`
and `q` in `k[t]` of degree at most `n` of `a*D(q)+b*q=∑ᵢcᵢ*qᵢ`, `p=(q-∑ᵢcᵢ*rᵢ)/a`
has degree at most `N` and satisfies `A*D(p)+B*p=∑ᵢcᵢ*Qᵢ`.

See [Bronstein's book](https://link.springer.com/book/10.1007/b138171), Section 7.1, p. 231.
"""
function ParSPDE(a::P, b::P, qs::Vector{P}, D::Derivation, n::Int) where
    {T<:FieldElement, P<:PolyElem{T}}
    iscompatible(a, D) && iscompatible(b, D) && all([iscompatible(q, D) for q in qs]) ||
        error("polynomial a and rational functions b and q_i must be in the domain of derivation D")
    degree(a)>0 || error("degree(a) must be > 0")
    degree(gcd(a, b))<=0 || error("gcd(a,b) must be 1")
    rzs = [gcdx(b, a, q) for q in qs]
    a, b + D(a), [z-D(r) for (r, z) in rzs], [r for (r, z) in rzs], n-degree(a)
end

"""
    ParamPolyRischDENoCancel1(b, qs, D, n) -> (hs, A)

Parametric polynomial Risch differential equation - no cancellation.

Given a field `k`, a derivation `D` on `k[t]`, `b≠0` in `k[t]`, `qs=[q₁,...,qₘ]` with
`qᵢ` in `k[t]`, and an integer `n` such that either
`D=d/dt` or `degree(b)>max(0, degree(D(t))-1)`, return `hs=[h₁,...,hᵣ]` with
`hᵢ` in `k[t]` and  a matrix `A` with entries in `Const(k)` such that
if `c₁`,...,`cₘ` in `Const(k)` and `q` in `k[t]` satisfy `degree(q)<=n` and
`D(q)+b*q=∑ᵢcᵢ*qᵢ` then `q=∑ⱼdⱼ*hⱼ` where `d₁`,...,`dᵣ` are in `Const(k)` and
`A*[c₁,...,cₘ,d₁,...,dᵣ]=0`.

See [Bronstein's book](https://link.springer.com/book/10.1007/b138171), Section 7.1, p. 234.
"""
function ParamPolyRischDENoCancel1(b::P, qs::Vector{P}, D::Derivation, n::Int) where
    {T<:FieldElement, P<:PolyElem{T}}
    # Note: this implementation changes the input parameters qs!
    iscompatible(b, D) && all([iscompatible(q, D) for q in qs]) ||
        error("polynomials a and q_i must be in the domain of derivation D")
    !iszero(b) || error("polynomial b must be nonzero")
    isbasic(D) || degree(b)>max(0, degree(D)-1) ||
        error("either derivation D must be basic or degree(b)>max(0, degree(D)-1)")
    t = gen(parent(b))
    db = degree(b)
    bd = leading_coefficient(b)
    m = length(qs)
    hs = [zero(b) for i=1:m]
    ss = [zero(bd) for i=1:m]
    while n>=0
        for i=1:m
            ss[i] = coeff(qs[i], n + db)//bd
            sitn = ss[i]*t^n
            hs[i] += sitn
            qs[i] -= D(sitn) + b*sitn
        end
        n -= 1
    end
    if all([iszero(q) for q in qs])
        dc = -1
    else
        dc = maximum([degree(q) for q in qs])
    end
    M = [coeff(q, i) for i=0:dc, q in qs]
    A = ConstantSystem(M, BaseDerivation(D))
    C = constant_field(D)
    neq = size(A, 1)
    A = vcat(hcat(A, zeros(C, neq, m)), zeros(C, m, 2*m))
    for i=1:m
        A[i+neq, i] = one(C)
        A[i+neq, m+i] = -one(C)
    end
    hs, A
end

"""
    ParamPolyRischDENoCancel2(b, qs, D, n) -> (hs, A)

Parametric polynomial Risch differential equation - no cancellation.

Given a field `k`, a derivation `D` on `k[t]`, `b` in `k[t]`, `qs=[q₁,...,qₘ]` with
`qᵢ` in `k[t]`, and an integer `n` such that `degree(b)<degree(D(t)-1)` and either
`D=d/dt` or `degree(D(t))≥2`, return `hs=[h₁,...,hᵣ]` with
`hᵢ` in `k[t]` and  a matrix `A` with entries in `Const(k)` such that
if `c₁`,...,`cₘ` in `Const(k)` and `q` in `k[t]` satisfy `degree(q)≤n` and
`D(q)+b*q=∑ᵢcᵢ*qᵢ` then `q=∑ⱼdⱼ*hⱼ` where `d₁`,...,`dᵣ` are in `Const(k)` and
`A*[c₁,...,cₘ,d₁,...,dᵣ]=0`.

See [Bronstein's book](https://link.springer.com/book/10.1007/b138171), Section 7.1, p. 238.
"""
function ParamPolyRischDENoCancel2(b::P,  qs::Vector{P}, D::Derivation, n::Int) where
    {T<:FieldElement, P<:PolyElem{T}}
    iscompatible(b, D) && all([iscompatible(q, D) for q in qs]) ||
        error("polynomials a and q_i must be in the domain of derivation D")
    δ = degree(D)
    iszero(b) || degree(b)<δ-1 || error("degree(b)<degree(D)-1 must hold")
    # Note: In Bronstein's book the convention degree(0)=-infinity is used,
    # but in AbstractAlgebra we have degree(0)=-1, so we have to explicitly
    # include the clause iszero(b) in the above condition.
    isbasic(D) || δ>=2 ||
        error("either derivation D must be basic or degree(D)>=2")
    t = gen(parent(b))
    λ = leading_coefficient(D)
    Z  = zero(b)
    Z0 = zero(λ)
    m = length(qs)
    hs = [zero(b) for i=1:m]
    ss = [zero(λ) for i=1:m]
    while n>0
        nλ = n*λ
        for i=1:m
            ss[i] = coeff(qs[i], n + δ -1)//(nλ)
            sitn = ss[i]*t^n
            hs[i] += sitn
            qs[i] -= D(sitn) + b*sitn
        end
        n -= 1
    end
    C = constant_field(D)
    if degree(b)>0
        for i=1:m
            ss[i] = coeff(qs[i], degree(b))//leading_coefficient(b)
            hs[i] += ss[i]
            qs[i] -= D(ss[i]) + b*ss[i]
        end
        if all([iszero(q) for q in qs])
            dc = -1
        else
            dc = maximum([degree(q) for q in qs])
        end
        M = [coeff(q, i) for i=0:dc, q in qs]
        A = ConstantSystem(M, BaseDerivation(D))
        neq = size(A, 1)
        A = vcat(hcat(A, zeros(C, neq, m)), zeros(C, m, 2*m))
        for i=1:m
            A[i+neq, i] = one(C)
            A[i+neq, m+i] = -one(C)
        end
        return hs, A
    else
        D0 = BaseDerivation(D)
        b0 = constant_coefficient(b) # b \in k
        fs, B = ParamRischDE(b0, [constant_coefficient(q) for q in qs], D0)
        if all([iszero(q) for q in qs])
            if all([iszero(D0(f)+b0*f) for f in fs])
                dc = -1
            else
                dc = 0
            end
        else
            dc = maximum([degree(q) for q in qs])
        end
        r = length(fs)
        M = zeros(parent(b0), dc+1, m+r)
        for i=0:dc
            for j = 1:m
                M[i+1,j] = coeff(qs[j], i)
            end
        end
        if dc>=0
            for j=1:r
                M[1,j+m] = -D0(fs[j]) - b0*fs[j]
            end
        end
        A = ConstantSystem(M, D0)
        A = vcat(A, B)
        neq = size(A, 1)
        A = vcat(hcat(A, zeros(C, neq, m)), zeros(C, m, 2*m+r))
        for i=1:m
            A[i+neq, i] = one(C)
            A[i+neq, m+r+i] = -one(C)
        end
        return vcat(fs .+ zero(b), hs), A
    end
end

"""
    ParamPolyRischDECancelLiouvillian(b, qs, D, n) -> (hs, A)

Parametric polynomial Risch differential equation - no cancellation.

Given a field `k`, a derivation `D` on `k[t]` with `D≠d/dt` and `D(t)` in `k`
or `D(t)/t` in `k`, `b` in `k`, `qs=[q₁,...,qₘ]` with
`qᵢ` in `k[t]`, and an integer `n`,   return `hs=[h₁,...,hᵣ]` with
`hᵢ` in `k[t]` and  a matrix `A` with entries in `Const(k)` such that
if `c₁`,...,`cₘ` in `Const(k)` and `q` in `k[t]` satisfy `degree(q)≤n` and
`D(q)+b*q=∑ᵢcᵢ*qᵢ` then `q=∑ⱼdⱼ*hⱼ` where `d₁`,...,`dᵣ` are in `Const(k)` and
`A*[c₁,...,cₘ,d₁,...,dᵣ]=0`.

See [Bronstein's book](https://link.springer.com/book/10.1007/b138171), Section 7.1, p. 241.
"""
function ParamPolyRischDECancelLiouvillian(b::T,  qs::Vector{P}, D::Derivation, n::Int) where
    {T<:FieldElement, P<:PolyElem{T}}
    isprimitive(D) || ishyperexponential(D) ||
    error("monomial of derivation D must be primitive or hyperexponential")
    D0 = BaseDerivation(D)
    iscompatible(b, D0) ||
        error("coefficient b must be in the domain of the base derivation of D")
    all([iscompatible(q, D) for q in qs]) ||
        error("polynomials q_i must be in the domain of derivation D")
    m = length(qs)
    C = constant_field(D)
    if iszero(b) && all([iszero(q) for q in qs])
        return P[], zeros(C, 0, m)
    end
    H = MonomialDerivative(D)
    if ishyperexponential(D)
        w = coeff(H, 1)
    else
        w = zero(parent(b))
    end
    t = gen(parent(H))
    first = true
    qns = copy(qs)
    hs = P[]
    while n>=0
        fns, An = ParamRischDE(b + n*w, [coeff(q, n) for q in qns], D0)
        hns = [f*t^n for f in fns]
        qns = vcat(qns, [-(D(h) + b*h) for h in hns]) # 2nd displayed eq. on p. 242 (sign of b wrong there !?!?)
        if first
            A = An
            hs = hns
            first = false
        else
            A = vcat( hcat(A, zeros(C, size(A, 1),   size(An,2) - size(A, 2) )),
                           An)
            hs = vcat(hs, hns)
        end
        n -= 1
    end
    r = length(hs)
    Dhbhs = [D(h) + b*h for h in hs]  # sign of b wrong on  p. 242 !?!?!
    if all([iszero(q) for q in qs]) && all([iszero(Dhbh) for Dhbh in Dhbhs])
        dc = -1
    else
        dc = max( maximum([degree(q) for q in qs]), maximum([degree(Dhbh) for Dhbh in Dhbhs]))
    end
    M = zeros(parent(b), dc+1, m+r)
    for i=0:dc
        for j = 1:m
            M[i+1,j] = coeff(qs[j], i)
        end
        for j=1:r
            M[i+1,j+m] = -coeff(Dhbhs[j], i)
        end
    end
    A0 = ConstantSystem(M, D0)
    A = vcat(A0, A)
    hs, A
end

"""
    ParamPolyRischDE(b, qs, D, n) -> (hs, A)

Parametric polynomial Risch differential equation.

Given a field `k`, a derivation `D` on `k[t]`, `b` in `k[t]`, `qs=[q₁,...,qₘ]` with
`qᵢ` in `k[t]` and an integer `n`, return `hs=[h₁,...,hᵣ]` with
`hᵢ` in `k[t]` and  a matrix `A` with entries in `Const(k)` such that
if `c₁`,...,`cₘ` in `Const(k)` and `q` in `k[t]` satisfy `degree(q)≤n` and
`D(q)+b*q=∑ᵢcᵢ*qᵢ` then `q=∑ⱼdⱼ*hⱼ` where `d₁`,...,`dᵣ` are in `Const(k)` and
`A*[c₁,...,cₘ,d₁,...,dᵣ]=0`.

See [Bronstein's book](https://link.springer.com/book/10.1007/b138171), Section 7.1, p. 241.
"""
function ParamPolyRischDE(b::P,  qs::Vector{P}, D::Derivation, n::Int) where
    {T<:FieldElement, P<:PolyElem{T}}
    iscompatible(b, D) || all([iscompatible(q, D) for q in qs]) ||
        error("coefficient b and polynomials q_i must be in the domain of derivation D")
    δ = degree(D)
    if !iszero(b) && (isbasic(D) || degree(b)>max(0, δ-1))
        return ParamPolyRischDENoCancel1(b, qs, D, n)
    elseif (iszero(b) || degree(b)<δ-1) && (isbasic(D) || δ>=2)
        return ParamPolyRischDENoCancel2(b, qs, D, n)
    elseif δ>=2 && degree(b)==δ-1
        if ishypertangent(D)
            throw(NotImplementedError("ParamPolyRischDE: no cancellation, degree(b)==δ-1, hypertangent case\n@ $(@__FILE__):$(@__LINE__)"))
        else
            H = MonomialDerivative(D)
            throw(NotImplementedError("ParamPolyRischDE: no cancellation, degree(b)==δ-1, monomial derivative $H\n@ $(@__FILE__):$(@__LINE__)"))
        end
    elseif isprimitive(D) || ishyperexponential(D)
        @assert δ<=1 && degree(b)<=0 && !isbasic(D)
        return ParamPolyRischDECancelLiouvillian(constant_coefficient(b), qs, D, n)
    elseif ishypertangent(D)
        throw(NotImplementedError("ParamPolyRischDE: hypertangent, cancellation case\n@ $(@__FILE__):$(@__LINE__)"))
    else
        H = MonomialDerivative(D)
        throw(NotImplementedError("ParamPolyRischDE: cancelling case, monomial derivative $H\n@ $(@__FILE__):$(@__LINE__)"))
    end
end

"""
    ParamPolyCoeffsRischDE(a, b, gs, D) -> (hs, A)

Parametric Risch differential equation with polynomial coefficients.

Given a field `k`, a derivation `D` on `k[t]`, `a`, `b` in `k[t]` with `a≠0`, and `gs=[g₁,...,gₘ]` with
`gᵢ` in `k(t)`, return `hs=[h₁,...,hᵣ]` with
`hᵢ` in `k(t)` and  a matrix `A` in reduced row echelon form with entries in `Const(k)` such that
`c₁`,...,`cₘ` in `Const(k)` and `q` in `k[t]` is a solution of
`a*D(q)+b*q=∑ᵢcᵢ*gᵢ` if and only if `q=∑ⱼdⱼ*hⱼ` for `d₁`,...,`dᵣ` in `Const(k)` with
`A*[c₁,...,cₘ,d₁,...,dᵣ]=0`.

See [Bronstein's book](https://link.springer.com/book/10.1007/b138171), Section 7.1, p. 222.
"""
function ParamPolyCoeffsRischDE(a::P, b::P, gs::Vector{F}, D::Derivation;
                                n::Int=typemin(Int), exit_if_small_nullspace::Bool=false) where
                                {T<:FieldElement, P<:PolyElem{T}, F<:FracElem{P}}
    iscompatible(a, D) && iscompatible(b, D) && all([iscompatible(g, D) for g in gs]) ||
        error("polynomials a and b and rational functions g_i must be in the domain of derivation D")
    m = length(gs)
    d = gcd(a, b)
    a = divexact(a, d)
    b = divexact(b, d)
    gs = [g//d for g in gs]
    qs, M = LinearConstraints(a, b, gs, D)
    A = ConstantSystem(M, BaseDerivation(D))
    if exit_if_small_nullspace
        A = RowEchelon(A)
        if size(A, 2) - size(A, 1) <= 1
            # rank(nullspace) <= 1
            return F[], A, qs, true
        end
    end
    if n==typemin(Int)
        n = ParamRdeBoundDegree(a, b, qs, D)
    end
    aprod = one(parent(a))
    Rs = zeros(parent(a//one(a)), m)
    while n>=0 && degree(a)>0
        a0 = a
        a, b, qs, rs, n = ParSPDE(a, b, qs, D, n)
        aprod *= a0
        for i=1:m
            Rs[i] = (Rs[i] + rs[i])//a0
        end
        d = gcd(a, b)
        a = divexact(a, d)
        b = divexact(b, d)
        gs = [q//d for q in qs]
        qs, M = LinearConstraints(a, b, gs, D)
        A = vcat(A, ConstantSystem(M, BaseDerivation(D)))
    end
    C = constant_field(D)
    if n<0
        hs = F[]
        A = vcat(A, ConstantSystem(RowEchelon([q//1 for i=1:1, q in qs]), D))
    else
        a = constant_coefficient(a)
        b = divexact(b, a)
        qs = [divexact(q, a) for q in qs]
        hs, A1 = ParamPolyRischDE(b, qs, D, n)
        A = vcat(hcat(A, zeros(C, size(A, 1), length(hs))), A1)
    end
    r = length(hs)
    s = sum([1 for R in Rs if !iszero(R)])
    hs = vcat([aprod*h//1 for h in hs], [aprod*R for R in Rs if !iszero(R)])
    neq = size(A, 1)
    A = vcat(hcat(A, zeros(C, neq, s)), zeros(C, s, m+r+s))
    i1 = 0
    for i=1:m
        if !iszero(Rs[i])
            i1 += 1
            A[neq+i1, i] = one(C)
            A[neq+i1, m+r+i1] = -one(C)
        end
    end
    A = RowEchelon(A)
    if exit_if_small_nullspace
        return hs, A, P[], false
    else
        hs, A
    end
end

"""
    ParamRischDE(f, gs, D) -> (hs, A)

Parametric Risch differential equation.

Given a field `K`, a derivation `D` on `K`, `f` in `K`, `gs=[g₁,...,gₘ]` with
`gᵢ` in `K`, return `hs=[h₁,...,hᵣ]` with
`hᵢ` in `K` and  a matrix `A` reduced row echelon form with entries in `Const(K)` such that
`c₁`,...,`cₘ` in `Const(K)` and `y` in `K` is a solution of
`D(y)+f*y=∑ᵢcᵢ*gᵢ` if and only if `y=∑ⱼdⱼ*hⱼ` for `d₁`,...,`dᵣ` in `Const(K)` with
`A*[c₁,...,cₘ,d₁,...,dᵣ]=0`.

See [Bronstein's book](https://link.springer.com/book/10.1007/b138171), Section 7.1, p. 217.
"""
function ParamRischDE(f::F, gs::Vector{F}, D::Derivation) where F<:FieldElement
    iszero(D) || error("base case only for null derivations")
    #base case => pure (linear) algebra problem ...
    iscompatible(f, D) && all(iscompatible(g, D) for g in gs) ||
        error("rational functions f and g_i must be in the domain of derivation D")
    m = length(gs)
    C = parent(f)
    if iszero(f)
        hs = [one(parent(f))]
        A = reshape(vcat(gs, zero(parent(f))), (1, m+1))
    else
        A = zeros(C, m , 2*m)
        for i=1:m
            A[i, i] = one(C)
            A[i, m+i] = -one(C)
        end
        hs = [gs[i]//f for i=1:m]
    end
    hs, RowEchelon(A)
end

function ParamRischDE(f::F, gs::Vector{F}, D::Derivation) where
    {T<:FieldElement, P<:PolyElem{T}, F<:FracElem{P}}
    iscompatible(f, D) && all(iscompatible(g, D) for g in gs) ||
        error("rational functions f and g_i must be in the domain of derivation D")

    h0 = WeakNormalizer(f, D)

    f1 = f - D(h0)//h0
    gs = [h0*g for g in gs]
    a, b, gs, h1 = ParamRdeNormalDenominator(f1, gs, D)

    a, b, gs, h2 =  ParamRdeSpecialDenominator(a, b, gs, D)

    hs, A = ParamPolyCoeffsRischDE(a, b, gs, D) # A in reduced echelon form

    h012 = h0*h1*h2
    [h//h012 for h in hs], A
end

"""
    LimitedIntegrateReduce(f, ws, D) -> (a, b, h, N, g, vs)

Reduction to a polynomial problem.

Given a field `k`, a derivation `D` on `k[t]`, `f` in `k(t)` and `ws=[w₁,...,wₘ]` with
`wᵢ` in `k(t)`, return `a`, `b`, `h` in `k[t]`, `g` in `k(t)`, `vs=[v₁,...,vₘ]` with
`vᵢ` in `k(t)` and an integer `N` such that
for any solution `v` in `k(t)`, `c_1`,...,`cₘ` in Const(k) of `f = D(v)+∑ᵢcᵢ*wᵢ`,
`p=v*h` is in `k⟨t⟩`, and `p` and the `c_i` satisfy `a*D(p)+b*p=g+∑ᵢcᵢ*vᵢ`.

(Here, `k⟨t⟩` denotes the elements of `k(t)` which are reduced w.r.t. `D`.)

Furthermore, if `S₁ⁱʳʳ==Sⁱʳʳ`, which is indicated by `is_Sirr1_eq_Sirr(D)==true`,
then `p` is in `k[t]`, and if `t` is nonlinear or Liouvillian over `k`, then `degree(p)≤N`.

See [Bronstein's book](https://link.springer.com/book/10.1007/b138171), Section 7.2, p. 248.
"""
function LimitedIntegrateReduce(f::F, ws::Vector{F}, D::Derivation) where
    {T<:FieldElement, P<:PolyElem{T}, F<:FracElem{P}}
    # Note: LimitedIntegrateReduce seems to be the only algorithm in Bronstein's book,
    # where he messed something up. Already in equation (7.32) of Corollary 7.2.1 (p.247)
    # the coefficient of D(p) contains a factor hs too many. This then reproduces in the
    # algoritm.
    iscompatible(f, D) && all(iscompatible(w, D) for w in ws) ||
        error("rational functions f and w_i must be in the domain of derivation D")
    dn, ds = SplitFactor(denominator(f), D)
    eness = [SplitFactor(denominator(w), D) for w in ws]
    c = lcm(vcat(dn, [en for (en, es) in eness]))
    hn = gcd(c, derivative(c))
    a = hn
    b = -D(hn)
    N = 0
    if is_Sirr1_eq_Sirr(D)
       hs = lcm(vcat(ds, [es for (en, es) in eness]))
       h = hn*hs # in Bronsteins's book, wrongly a = hn*hs
       b = -D(hn)-divexact(hn*D(hs),hs)
       μ = min(iszero(f) ? typemax(Int) : valuation_infinity(f), minimum([iszero(w) ? typemax(Int) : valuation_infinity(w) for w in ws]))
       N = degree(hn) + degree(hs) + max(0, 1 - degree(D) - μ)
    end
    hhn = h*hn
    # in Bronsteins'book, wrongly always the factor a*hn instead of h*hn
    a, b, h, N, hhn*f, [-hhn*w for w in ws]
end

function solve_x1_eq_1(A::Matrix{T}) where T<:FieldElement
    # return solution x of A*x = 0 with x[1]==1 if it exists and x=0 otherwise.
    # A must be in  reduced row echelon form
    nr, nc = size(A)
    @assert nr>0 && nc>0
    Z = zero(A[1,1])
    E = one(A[1,1])
    x = fill(Z, nc)
    if iszero(A[1,1])
        x[1] = E
        return x
    elseif (isone(A[1,1]) && (all([iszero(a) for a in A[1,2:nc]])))
        return x
    else
        # A[1,1]==1 && not all other elements in the first row are ==0
        k = findfirst(a->!iszero(a), A[1, 2:nc])+1
        p = [findfirst(a->isone(a), A[i,:]) for i=1:nr]
        x[k] = E
        x[p] = -A*x
        @assert !iszero(x[1])
        x .//= x[1]
        return x
    end
end

"""
    LimitedIntegrate(f, ws, D) -> (v, cs, ρ)

Limited integration problem.

Given a field `k`, a derivation `D` on `k[t]`, `f` in `k(t)`, `ws=[w₁,...,wₘ]` with
`wᵢ` in `k(t)`, return either `ρ=0`, in which case the equation
`f = D(v)+∑ᵢcᵢ*wᵢ` has no solution `v` in `k(t)`, `cs=[c₁,...,cₘ]`, `cᵢ` in `Const(k)`,
or `ρ=1` and a solution `v`, `cs` of that equation.

See [Bronstein's book](https://link.springer.com/book/10.1007/b138171), Section 7.2, p. 245.
"""
function LimitedIntegrate(f::F, w::F, D::Derivation) where F <: FieldElement
    # this is a wrapper for the case length(ws)==1 for convenience
    y, cs, ρ = LimitedIntegrate(f, [w], D)
    if ρ <=0
        return zero(f), zero(constant_field(D)), ρ
    else
        return y, cs[1], 1
    end
end

function LimitedIntegrate(f::F, ws::Vector{F}, D::Derivation) where F <: FieldElement
    # base case f, ws[i] \in constant field, D must be the null derivation
    iscompatible(f, D) && all(iscompatible(w, D) for w in ws) ||
        error("rational functions f and w_i must be in the domain of derivation D")
    iszero(D) || error("base case only for null derivations")
    Z = zero(f)
    m = length(ws)
    A = reshape(vcat(-f, ws), (1, m+1))
    A = RowEchelon(A)
    cs = solve_x1_eq_1(A)
    if isone(cs[1])
        return Z, cs[2:end], 1
    else
        return Z, fill(Z, m), 0
    end
end

function LimitedIntegrate(f::F, ws::Vector{F}, D::Derivation) where
    {T<:FieldElement, P<:PolyElem{T}, F<:FracElem{P}}
    iscompatible(f, D) && all(iscompatible(w, D) for w in ws) ||
        error("rational functions f and w_i must be in the domain of derivation D")
    Z = zero(f)
    no_solution = (Z, F[], 0)
    (a, b, h, n, g, vs) = LimitedIntegrateReduce(f, ws, D)
    if is_Sirr1_eq_Sirr(D)
        gs = vcat(g, vs)
        # Note:  ParamPolyCoeffsRischDE does the following division by gcd(a,b) anyway.
        # But in case of small_nullspace==true we need these *reduced* a, b, so we compute them here.
        d = gcd(a, b)
        a = divexact(a, d)
        b = divexact(b, d)
        gs = [g//d for g in gs]
        hs, A, qs, small_nullspace = ParamPolyCoeffsRischDE(a, b, gs, D, n=n, exit_if_small_nullspace=true)
        # A in reduced echelon form
        if !small_nullspace # regular exit
            cds = solve_x1_eq_1(A) # solution of A*(cs, ds)=0 with cs[1]==1
            if iszero(cds[1])
                return no_solution
            end
            @assert isone(cds[1])
            m = length(gs)
            r = length(hs)
            cs = cds[1:m]
            ds = cds[m+1:m+r]
            y = Z
            for j=1:length(hs)
                y +=  ds[j]*hs[j]
            end
            y = y//h
            return y, cs[2:m], 1
        end
        # premature exit
        rg = size(A, 2)-size(A, 1)
        if rg==0
            return no_solution
        end
        @assert rg==1
        cs = solve_x1_eq_1(A) # solution of A*cs=0 with cs[1]==1
        if iszero(cs[1])
            return no_solution
        end
        m = length(qs)
        c = sum([cs[i]*qs[i] for i=1:m])
        b, c, n, α, β, ρ = SPDE(a, b, c, D, n)
        ρ>=1 || return Z, F[], ρ
        z, ρ = PolyRischDE(b, c, D, n)
        ρ>=1 || return Z, F[], ρ
        return (α*z+β)//h, cs[2:m], 1
    else
        throw(NotImplementedError("LimitedIntegrate: !is_Sirr1_eq_Sirr(D)\n@ $(@__FILE__):$(@__LINE__)"))
    end
end


text:
import Base.show

# fallback for Abstract
function show(io::IO, a::AbstractEdge)
    print(io, "$(a.edge)")
end

function show(io::IO, a::AbstractNode)
    print(io, "$(a.node)")
end

function show(io::IO, a::SimpleNode)
    print(io, "$(a.node)")
end

function show(io::IO, a::ModifyingNode)
    print(io, "$(a.node) * $(a.modifiers)")
end

function show(io::IO, a::DirectedEdge)
    print(io, "$(a.src) → $(a.dst)")
end

function show(io::IO, a::UndirectedEdge)
    print(io, "$(a.src) ↔ $(a.dst)")
end

function show(io::IO, a::Vector{E}) where {E <: AbstractEdge}
    for e in a
       print(io, "$(e)\n")
    end
end

function show(io::IO, ::MIME"text/plain", a::Vector{E}) where {E <: AbstractEdge}
    show(io, a)
end

function show(io::IO, a::Arrow)
    print(io, "$(a.edges)")
end

function show(io::IO, a::ModifiedEdge)
    print(io, "$(a.edge) * $(a.modifiers)")
end

function show(io::IO, a::ModifiedNode)
    print(io, "$(a.node)^$(a.modifiers)")
end

function show(io::IO, a::Dict{Symbol, M}) where {M <: Modifier}
    if length(a) > 0
        keyset = [k for k in keys(a)]
        if length(keyset) > 1
            print(io, "[")
            for i in 1:(length(keyset) - 1)
                k = keyset[i]
                print(io, "$(a[k]), ")
            end
            print(io, "$(a[last(keyset)])]")
        else
            print(io, "$(a[last(keyset)])")
        end
    else
        print(io, nothing)
    end
end

function show(io::IO, ::MIME"text/plain", a::Dict{Symbol, M}) where {M <: Modifier}
    show(io, a)
end


text:
"""
    getchannelhistory(token::String)
Takes in the Slack App Auth token and channel ID.

Returns the 800 most recent messags (this may not be the case always). Could be time based...

Tokens can be accessed at "https://api.slack.com/apps/random-letters-and-numbers/oauth?"
"""
function getchannelhistory(token::String, channelid::String)

    HTTP.setuseragent!("HTTP.jl")
    url = "https://slack.com/api/conversations.history?token=$(token)&channel=$(channelid)"
    r = HTTP.get(url, ["Content-Type" => "application/json"])

    json_obj = JSON.parse(String(r.body))
    message_holder = []

    #TODO: Need to figure out a parameter to constrain the time period in which messages are returned.

    for (k,v) in json_obj
        if occursin("messages", k)
            for item in v

                # Declare the holder variables
                edited = false
                ts = ""
                thread_ts = ""
                user = ""
                text = ""

                # println("Starting to parse a new message") #DEBUG
                for (key,value) in item
                    if key == "edited"
                        edited = true

                    # TODO: also convert this to Date Time and add that for print out.
                    elseif key == "ts"
                        ts = value

                    elseif key == "text"
                        text = value

                    elseif key == "user"
                        user = value

                    # It's a thread not a message...
                    elseif key == "thread_ts"
                        thread_ts = value
                    end
                end

                # If ts and thread_ts are the same, then this is a thread so call the thread function
                if ts == thread_ts && ts != ""

                    try
                        thread_msg = getthreads(token, channelid, ts)
                        push!(message_holder, thread_msg)
                    catch
                        println("Failed to create a threaded message object")
                    end

                # Else it's a normal message so create it here and push it into an array.
                else
                    msg = message(edited, ts, user, text)
                    push!(message_holder, msg)
                end
            end
        end
    end

    return(message_holder)

end