text:
"""
Water flooding with BL assumptions
"""
function water_flood_numeric(core_props, fluids, rel_perms, core_flood; Nx = 30)
    pv_inj = core_flood.injected_pore_volume

    # Nx = 100 # number of cells in x direction
    H = core_props.length # [m] length of the domain in y direction
    D_core = core_props.diameter
    A_core=π*D_core^2/4
    V_core=H*A_core
    m = createMesh1D(Nx, H)

    ## define the physical parametrs
    # all the Corey-type relperm parametrs are defined for the oil-wet and water-wet
    # cases
    k0 = core_props.permeability # [m^2] average reservoir permeability
    phi0 = core_props.porosity # average porosity
    mu_oil = fluids.oil_viscosity # [Pa.s] oil viscosity
    mu_water = fluids.water_viscosity # [Pa.s] water viscosity
    krw0 = rel_perms.krw0
    kro0 = rel_perms.kro0
    nw = rel_perms.nw
    no = rel_perms.no
    sor=rel_perms.sor
    swc=rel_perms.swc

    p0 = core_flood.back_pressure # [bar] pressure
    u_in= core_flood.injection_velocity # [m/s] equal to 1 m/day
    sw0 = core_flood.initial_water_saturation

    sw_in = core_flood.injected_water_saturation

    eps1=1e-5
    perm_val= k0 #field2d(Nx,Ny,k0,V_dp,clx,cly)

    k=createCellVariable(m, perm_val)
    phi=createCellVariable(m, phi0)

    lw = geometricMean(k)/mu_water
    lo = geometricMean(k)/mu_oil

    ## Define the boundaries: all fixed Sw=1, fixed pressure everywhere(?)
    BCp = createBC(m) # Neumann BC for pressure
    BCs = createBC(m) # Neumann BC for saturation
    BCc = createBC(m) # Neumann BC for salinity

    BCp.left.a[:].=-k0/mu_water
    BCp.left.b[:].=0
    BCp.left.c[:].=u_in

    BCp.right.a[:].=0.0
    BCp.right.b[:].=1.0
    BCp.right.c[:].=p0

    BCs.left.a[:].=0.0
    BCs.left.b[:].=1.0
    BCs.left.c[:].=sw_in

    BCc.left.a[:].=0.0
    BCc.left.b[:].=1.0
    BCc.left.c[:].=1.0 # injection concentration

    ## define the time step and solver properties
    t_end = pv_inj*V_core*phi0/(u_in*A_core) # [s] final time
    dt=H/Nx/u_in/10
    SW0 = sw0
    eps_p = 1e-3 # pressure accuracy
    eps_sw = 1e-6 # saturation accuracy
    ## define the variables
    sw_old = createCellVariable(m, SW0, BCs)
    sw_new=createCellVariable(m, SW0, BCs)
    p_old = createCellVariable(m, p0, BCp)
    p_new= createCellVariable(m, p0, BCp)
    c_old = createCellVariable(m, 0.0, BCc)
    c_face=arithmeticMean(c_old) # just to initialize

    # define the rel perm functions
    KRW = sw -> CF.krw(sw, krw0, sor, swc, nw)
    KRO = sw -> CF.kro(sw, kro0, sor, swc, no)
    dKRW = sw -> CF.dkrwdsw(sw, krw0, sor, swc, nw)
    dKRO = sw -> CF.dkrodsw(sw, kro0, sor, swc, no)

    sw = copyCell(sw_old)
    oil_init=domainInt(1-sw_old)
    p = copyCell(p_old)
    uw = -gradientTerm(p_old) # an estimation of the water velocity
    ## start the main loop
    # generate intial pressure profile (necessary to initialize the fully
    # implicit solver)
    dp_core = zeros(1)
    dp_core[1]=H*u_in/(k0*(CF.kro(sw0, kro0, sor, swc, no)/mu_oil+
        CF.krw(sw0, krw0, sor, swc, nw)/mu_water))
    rec_fact=zeros(1)
    c_out_sal=zeros(1)
    c_out_sal[1]=c_face.xvalue[end]
    t_sec=zeros(1)
    t = 0.0
    dt0=dt
    dsw_alwd= 0.005
    dp_alwd= 100.0 # Pa
    dc_alwd= 0.005
    labdaw=createFaceVariable(m,1.0) # initialize
    FL=fluxLimiter("SUPERBEE")

    prog_1=Progress(100, 1)

    while (t<t_end)
        # Implicit loop
        error_p = 1e5
        error_sw = 1e5
        # implicit solver
        loop_count=0
        while ((error_sw>eps_sw))
            loop_count+=1
            if loop_count>10
                break
            end
            # calculate parameters
            pgrad = gradientTerm(p)
            uw = -labdaw*pgrad
            sw_face = upwindMean(sw, uw) # average value of water saturation
            # sw_ave=arithmeticMean(sw)
            # solve for pressure at known Sw
            labdao = lo*faceEval(KRO, sw_face)
            labdaw = lw*faceEval(KRW, sw_face)
            dlabdaodsw = lo*faceEval(dKRO, sw_face)
            dlabdawdsw = lw*faceEval(dKRW, sw_face)
            labda = labdao+labdaw
            dlabdadsw = dlabdaodsw+dlabdawdsw
            # compute [Jacobian] matrices
            Mdiffp1 = diffusionTerm(-labda)
            Mdiffp2 = diffusionTerm(-labdaw)
            # Mconvsw1 = convectionUpwindTerm(-dlabdadsw.pgrad, uw)
            # Mconvsw2 = convectionUpwindTerm(-dlabdawdsw*pgrad, uw)
            Mconvsw1 = convectionUpwindTerm(-dlabdadsw*pgrad, uw)
            Mconvsw2 = convectionUpwindTerm(-dlabdawdsw*pgrad, uw)
            Mtranssw2, RHStrans2 = transientTerm(sw_old, dt, phi)
            # Compute RHS values
            RHS1 = Mconvsw1*sw.value # divergenceTerm(-dlabdadsw*sw_face*pgrad)
            RHS2 = Mconvsw2*sw.value # divergenceTerm(-dlabdawdsw*sw_face*pgrad)
            # include boundary conditions
            Mbcp, RHSbcp = boundaryConditionTerm(BCp)
            Mbcsw, RHSbcsw = boundaryConditionTerm(BCs)
            # Couple the equations; BC goes into the block on the main diagonal
            M = [Mdiffp1+Mbcp Mconvsw1; Mdiffp2 Mconvsw2+Mtranssw2+Mbcsw]
            RHS = [RHS1+RHSbcp; RHS2+RHStrans2+RHSbcsw]
            x = M\RHS
            p_new.value[:] = (x[1:(Nx+2)])
            sw_new.value[:] = (x[(Nx+2)+1:end])

            error_p = maximum(abs.((internalCells(p_new)-internalCells(p))./internalCells(p_new)))
            error_sw = maximum(abs.(internalCells(sw_new)-internalCells(sw)))
            # dt_new=dt*min(dp_alwd/error_p, dsw_alwd/error_sw)
            # print("sw_error = $error_sw \n")
            # assign new values of p and sw
            p.value[:]=p_new.value[:]
            sw.value[:]=sw_new.value[:]
        end # end of implicit loop
        if loop_count>10
            p=copyCell(p_old)
            sw=copyCell(sw_old)
            dt=dt/5
            continue
        end
        # solve the advection equation
        # ϕ d/dt(c sw)=[sw dc/dt + c dsw/dt]
        dswdt=(sw_new-sw_old)/dt
        Msc=linearSourceTerm(phi*dswdt)
        sw_face = upwindMean(sw, -labdaw*pgrad) # average value of water saturation
        pgrad = gradientTerm(p_new)
        labdaw = lw*faceEval(KRW, sw_face)
        uw=-labdaw*pgrad
        Mtc, RHStc=transientTerm(c_old, dt, phi*sw)
        Mbcc, RHSbcc=boundaryConditionTerm(BCc)
        c_temp=copyCell(c_old)
        for j in 1:4
            Mconvc, RHSconv=convectionTvdTerm(uw, c_temp, FL)
            c_temp=solveLinearPDE(m, Msc+Mtc+Mconvc+Mbcc, RHStc+RHSbcc+RHSconv)
        end
        c_new=copyCell(c_temp)
        c_face=tvdMean(c_new, uw, FL)

        # adaptive time step
        dp = maximum(abs.((internalCells(p_new)-internalCells(p_old))./internalCells(p_new)))
        dsw = maximum(abs.(internalCells(sw_new)-internalCells(sw_old)))
        dc = maximum(abs.(internalCells(c_new)-internalCells(c_old)))

        # update the variables
        t+=dt
        ProgressMeter.update!(prog_1, floor(Int, t/t_end*100))
        c_old=copyCell(c_new)
        p_old = copyCell(p)
        sw_old = copyCell(sw)

        dt=min(dt*(dsw_alwd/dsw), 15*dt, t_end-t-eps())

        # calculate the recovery factor
        push!(rec_fact, (oil_init-domainInt(1-sw))/oil_init)
        push!(t_sec, t)
        push!(c_out_sal, c_face.xvalue[end])
        push!(dp_core, 0.5*sum(p.value[1:2])-p0)

        # print(t)
        #GR.plot(sw.value[2:end-1])
        # GR.plot(c_new.value[2:end-1])
        #GR.plot(SF.xvalue)
        # GR.imshow(sw.value[2:end-1,2:end-1])
        # visualizeCells(1-sw)
        # GR.plot(t_sec/3600/24, rec_fact)
        # xlabel('time [day]')
        # ylabel('recovery factor')
        # title([num2str(t/3600/24) ' day']) drawnow
    end
    pv = t_sec*u_in/(H*phi0)
    sw_face = linearMean(sw)
    c_face  = linearMean(c_old)
    x = m.facecenters.x
    t_sec, pv, rec_fact, dp_core, x, sw_face.xvalue, c_face.xvalue, c_out_sal
end # imb_impes

max_stars_repo_name: simulkade/peteng

max_issues_repo_name: simulkade/peteng

max_forks_repo_name: simulkade/peteng


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule ElasticFusionLib_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("ElasticFusionLib")
JLLWrappers.@generate_main_file("ElasticFusionLib", UUID("36eef7af-e89c-5840-8060-ddffb8b08188"))
end  # module ElasticFusionLib_jll

ext: jl


text:
using BlackBoxOptim

function fitness_for_opt(problem, numDimensions, populationSize, numSteps, method)

  problem = BlackBoxOptim.as_fixed_dim_problem(problem, numDimensions)

  println("\n$(problem.name), n = $(numdims(problem)), optimizer = $(string(method))")

  best, fitness = bboptimize(problem; method = method, parameters = {
    :NumDimensions => numDimensions,
    :PopulationSize => populationSize,
    :MaxSteps => numSteps
    })

  fitness
end


text:
"""
A structure containig the informations about the resulot of a ClaDS run. Contains the following fields :

- `tree::Tree`: the phylogeny on which the analysis was performed.
- `chains::Array{Array{Array{Float64,1},1},1}` : the mcmc chains
- `rtt_chains::Array{Array{Array{Float64,1},1},1}` : the mcmc chains with the rate through time
- `σ_map::Float64` : the σ parameter estimate
- `α_map::Float64` : the α parameter estimate
- `ε_map::Float64` : the ε parameter estimate
- `λ0_map::Float64` : the initial speciation rate estimate
- `λi_map::Array{Float64,1}` : the estimates of the branh specific speciation rates
- `λtip_map::Array{Float64,1}` : the estimates of the tip speciation rates
- `DTT_mean::Array{Float64,1}` : the diversity through time estimates
- `RTT_map::Array{Float64,1}` : the rate through time estimates
- `time_points::Array{Float64,1}` : the times at which `DTT_mean` and `RTT_map` are computed
- `enhanced_trees::Array{Tree,1}` : a sample of the complete tree distribution
- `gelm::Tuple{Int64,Float64}` : the gelman parameter
- `current_state` : other variables, by `infer_ClaDS` used to continue the run
"""
struct CladsOutput
    tree::Tree
    chains::Array{Array{Array{Float64,1},1},1}
    rtt_chains::Array{Array{Array{Float64,1},1},1}
    σ_map::Float64
    α_map::Float64
    ε_map::Float64
    λ0_map::Float64
    λi_map::Array{Float64,1}
    λtip_map::Array{Float64,1}
    DTT_mean::Array{Float64,1}
    RTT_map::Array{Float64,1}
    time_points::Array{Float64,1}
    enhanced_trees::Array{Tree,1}
    gelm::Tuple{Int64,Float64}
    current_state # all the rest that is mising to continue the run
end

function CladsOutput()
    return CladsOutput(
        Tree(),
        Array{Array{Array{Float64,1},1},1}(undef,0),
        Array{Array{Array{Float64,1},1},1}(undef,0),
        0.,0.,0.,0.,[0.],[0.],[0.],[0.],
        [0.],[Tree()],(0,10.),0.
    )
end

function print_CladsOutput(co::CladsOutput)
    n = n_tip(co.tree)
    nit = length(co.chains[1][1])
    to_print = string("ClaDS run for a $n tipped tree","\n",
        "$nit iterations, the gelman statistics in $(co.gelm[2])","\n",
        "inferred parameters : σ = $(co.σ_map), α = $(co.α_map), ε = $(co.ε_map), λ0 = $(co.λ0_map)"
        )
    println(to_print)
end

"""
    plot_CladsOutput(co::CladsOutput ; method = "tree", ...)

Plots various aspects of the output of ClaDS

# Arguments
- `co::CladsOutput` : A `CladsOutput` object, the output of a ClaDS run.

# Keyword arguments
- `method::String` : A `String` indicating what aspect of the output should be plotted, see details.
"""
function plot_CladsOutput(co::CladsOutput ; method = "tree",
    nplot = 50, alpha_col = 0.05, options = "", id_par = 1, force_ylim = Array{Float64,1}(undef,0))

    if method == "tree"
        plot_ClaDS(co.tree,co.λi_map, options = options)
    elseif method == "DTT"
        plot_DTT(co, n_ltt = nplot, alpha_col = alpha_col, force_ylim = force_ylim)
    elseif method == "RTT"
        plot_RTT(co, nplot = nplot, alpha_col = alpha_col, force_ylim = force_ylim)
    elseif method == "density"
        plot_density(co, id_par)
    elseif method == "chain"
        plot_chain(co, id_par)
    end
end

function tip_rate(co::CladsOutput, sp_name::String)
    number = 0
    dist = Inf
    tip_labs = tip_labels(co.tree)
    for itl in 1:length(tip_labs)
        d = Levenshtein()(sp_name, tip_labs[itl])
        if d < dist
            dist = d
            number = itl
            if d == 0
                break
            end
        end
    end
    return co.λtip_map[number]
end

max_stars_repo_name: hmorlon/JPANDA

max_issues_repo_name: hmorlon/JPANDA

max_forks_repo_name: hmorlon/JPANDA


text:
using CausalTrees
using Test

@testset "CausalTrees.jl" begin
    @test CausalTrees.update_average(2.5,4.0,2) == 3.0
    @test CausalTrees.update_average(3.0,4.0,3,true) == 2.5
end


text:
@defstruct TiedInnerProductLayer Layer (
  (name :: AbstractString = "", !isempty(name)),
  param_key :: AbstractString = "",
  (tied_param_key :: AbstractString = "", !isempty(tied_param_key)),
  (tops :: Vector{Symbol} = Symbol[], length(tops) > 0),
  (bottoms :: Vector{Symbol} = Symbol[], length(bottoms) == length(tops)),
  bias_init :: Initializer = ConstantInitializer(0),
  bias_regu :: Regularizer = NoRegu(),
  bias_cons :: Constraint = NoCons(),
  bias_lr :: AbstractFloat = 2.0,
  neuron :: ActivationFunction = Neurons.Identity()
)
@characterize_layer(TiedInnerProductLayer,
  can_do_bp  => true,
  has_param  => true,
  has_neuron => true
)

type TiedInnerProductLayerState <: LayerState
  layer      :: TiedInnerProductLayer
  blobs      :: Vector{Blob}
  blobs_diff :: Vector{Blob}

  parameters :: Vector{Parameter}

  W  :: Blob # reference to weights of tied layer, no learning for W here

  b  :: Blob
  ∇b :: Blob

  # a all-1 vector used in gemm to help bias calculation
  bias_multipliers :: Vector{Blob}

  frozen :: Bool

  TiedInnerProductLayerState(backend::Backend, layer::TiedInnerProductLayer, shared_params, inputs::Vector{Blob}) = begin
    fea_size  = get_fea_size(inputs[1])
    data_type = eltype(inputs[1])
    for i = 2:length(inputs)
      @assert get_fea_size(inputs[i]) == fea_size
      @assert eltype(inputs[i]) == data_type
    end

    tied_params = registry_get(backend, layer.tied_param_key)
    @assert length(tied_params) == 2
    @assert tied_params[1].name == "weight"
    W = tied_params[1].blob
    @assert size(W, 2) == fea_size
    out_dim = size(W, 1)

    blobs = Array{Blob}(length(inputs))
    blobs_diff = Array{Blob}(length(inputs))
    bias_multipliers = Array{Blob}(length(inputs))

    for i = 1:length(inputs)
      nums = get_num(inputs[i])
      blobs[i] = make_blob(backend, data_type, out_dim, nums)
      blobs_diff[i] = make_blob(backend, data_type, out_dim, nums)
      bias_multipliers[i] = make_blob(backend, ones(data_type, 1, nums))
    end

    state = new(layer, blobs, blobs_diff)
    state.bias_multipliers = bias_multipliers

    if shared_params != nothing
      @assert length(shared_params) == 1
      @assert shared_params[1].name == "bias"
      @assert size(shared_params[1].blob) == (out_dim, 1)
      m_debug("TiedInnerProductLayer($(layer.name)): sharing bias")

      params = [share_parameter(backend, shared_params[1])]
    else
      params = [make_parameter(backend, "bias", data_type, (out_dim,1),
          layer.bias_init, layer.bias_regu, layer.bias_cons, layer.bias_lr)]
    end

    state.W = W
    state.b = params[1].blob
    state.∇b = params[1].gradient
    state.parameters = params
    state.frozen = false

    return state
  end
end

function freeze!(state::TiedInnerProductLayerState)
  state.frozen = true
end
function unfreeze!(state::TiedInnerProductLayerState)
  state.frozen = false
end
function is_frozen(state::TiedInnerProductLayerState)
  state.frozen
end

function setup(backend::Backend, layer::TiedInnerProductLayer, shared_state, inputs::Vector{Blob}, diffs::Vector{Blob})
  TiedInnerProductLayerState(backend, layer, shared_state, inputs)
end

function setup(backend::Backend, layer::TiedInnerProductLayer, inputs::Vector{Blob}, diffs::Vector{Blob})
  setup(backend, layer, nothing, inputs, diffs)
end

function shutdown(backend::Backend, state::TiedInnerProductLayerState)
  map(destroy, state.blobs)
  map(destroy, state.blobs_diff)
  map(destroy, state.bias_multipliers)
  map(destroy, state.parameters)
end

function forward(backend::CPUBackend, state::TiedInnerProductLayerState, inputs::Vector{Blob})
  # note state.W is literally from the tied layer, we should use its transpose
  recon_dim  = size(state.W, 1)
  hidden_dim = size(state.W, 2)

  dtype = eltype(state.W)
  for i = 1:length(inputs)
    input = inputs[i]
    N = get_num(input)
    output = state.blobs[i]
    # output = (W^T)^T * X
    BLAS.gemm!('N', 'N', one(dtype), state.W.data,
                reshape(input.data, (hidden_dim,N)), zero(dtype), output.data)
    # output += bias
    BLAS.gemm!('N', 'N', one(dtype), state.b.data,
                state.bias_multipliers[i].data, one(dtype), output.data)
  end
end

function backward(backend::CPUBackend, state::TiedInnerProductLayerState, inputs::Vector{Blob}, diffs::Vector{Blob})
  recon_dim = size(state.W, 1)
  hidden_dim = size(state.W, 2)
  data_type  = eltype(state.W)

  # used in BLAS, at first it is zero, indicating overwriting the data
  # then it becomes one, indicating adding to the data
  zero_and_then_one = zero(data_type)

  for i = 1:length(inputs)
    input = inputs[i]
    batch_size = get_num(input)
    ∂f_∂o = state.blobs_diff[i]

    if !state.frozen
      # ∂f/∂b = sum(∂f/∂o, 2)
      BLAS.gemm!('N', 'N', one(data_type), ∂f_∂o.data,
                 reshape(state.bias_multipliers[i].data, (batch_size, 1)), zero_and_then_one,
                 state.∇b.data)
    end

    zero_and_then_one = one(data_type)

    # if back propagate down
    if isa(diffs[i], CPUBlob)
      # ∂f/∂x = W^T * [∂f/∂o]
      BLAS.gemm!('T', 'N', one(data_type), state.W.data,
                 ∂f_∂o.data, zero(data_type),
                 reshape(diffs[i].data, (hidden_dim, batch_size)))
    end
  end
end


text:
module GasChromatographySimulator

using Reexport
@reexport using Interpolations
using QuadGK
@reexport using DifferentialEquations
using ForwardDiff
@reexport using DataFrames
@reexport using CSV
using Plots
using HypertextLiteral
using PlutoUI

# some constants
Tst = 273.15            # K
R = 8.31446261815324    # J mol⁻¹ K⁻¹
Tn = 25.0 + Tst         # K
pn = 101300             # Pa

# ---Begin-Structures---
"""
    Column(L, d, a_d, df, a_df, sp, gas)

Structure describing the GC Column.

# Arguments

* `L`: Length of the capillary measured in m (meter)
* `d`: A function `d(x, a_d)` of `x`, the position along the capillary, describing the diameter in m (meter).
* `a_d`: Parameters of the diameter function.
* `d_f`: A function `d_f(x, a_df)` of `x`, describing the film thickness in m (meter).
* `sp`: The name of the stationary phase.
* `gas`: The name of the mobile phase. Allowed values: He, H2 or N2.
"""
struct Column{Fd<:Function, Fdf<:Function}
    L                       # column length in m
    d::Fd                   # column internal diameter in m as function of x
    a_d::Array{Float64,1}   # parameters of the diameters function d(x)
    df::Fdf                 # column film thickness in m as function of x
    a_df::Array{Float64}    # parameters of the film thickness function df(x)
    sp::String              # stationary phase of the column
    gas::String             # gas of the mobile phase ["He", "H2", "N2"]
    Column(L, d, a_d, df, a_df, sp, gas) = (gas!="He" && gas!="H2" && gas!="N2") ? error("Wrong selection for 'gas'. Choose 'He', 'H2' or 'N2'.") : new{typeof(d), typeof(df)}(L, d, a_d, df, a_df, sp, gas)
end

"""
    Program(time_steps, temp_steps, Fpin_steps, pout_steps, gf, a_gf, T_itp, Fpin_itp, pout_itp)

Structure to describe the temperature and flow/pressure program of a GC Column. The function `gf` describes an optional thermal gradient.

# Arguments
* `time_steps`: Time steps in s (seconds).
* `temp_steps`: Temperature steps in °C (degree celsius).
* `Fpin_steps`: Flow steps in m³/s resp. inlet pressure steps in Pa(a) (pascal, absolute).
* `pout_steps`: Outlet pressure steps in Pa(a) (pascal, absolute).
* `gf`: Gradient function `gf(x, a_gf)`, describes the thermal gradient.
* `a_gf`: Parameters of the gradient function.
* `T_itp`: Interpolated (linear) temperature `T(x,t)`, constructed from `time_steps`, `temp_steps` and `gf`.
* `Fpin_itp`: Interpolated (linear) flow/inlet pressure `Fpin(t)`, constructed from `time_steps` and `Fpin_steps`.
* `pout_itp`: Interpolated (linear) outlet pressure `pout(t)`, constructed from `time_steps` and `pout_steps`.

Note: The length of the arrays `time_steps`, `temp_steps`, `pin_steps` and `a_gf`
have to be the same.
"""
struct Program{Fgf<:Function}
    time_steps::Array{<:Real, 1}            # vector time steps for the temperature program
    temp_steps::Array{<:Real, 1}            # vector temperature steps for the temperature program
    Fpin_steps::Array{<:Real, 1}            # vector flow resp. inlet pressure steps for the pressure program
    pout_steps::Array{<:Real, 1}            # vector outlet pressure steps for the pressure program
    gf::Fgf                                 # function of x of the gradient form
    a_gf::Array{<:Real}                     # parameters of the gradient function gf(x)
  	T_itp::Interpolations.Extrapolation     # interpolation function of T(x,t)
    Fpin_itp::Interpolations.Extrapolation  # interpolation function of Fpin(t)
    pout_itp::Interpolations.Extrapolation	# interpolation function of pout(t)
    # add inner constructor to check the lengths of the Arrays and of the result of gf
    Program(ts, Ts, Fpis, pos, gf, agf, T, Fpin, po) = (length(ts)!=length(Ts) || length(ts)!=length(gf(0.0)) || length(ts)!=length(Fpis) || length(ts)!=length(pos)) || length(ts)!=size(agf)[1] ? error("Mismatch between length(time_steps) = $(length(ts)), length(temp_steps) = $(length(Ts)), length(Fpin_steps) = $(length(Fpis)), length(pout_steps) = $(length(pos)), length(gf(0.0)) = $(length(gf(0.0))) and size(a_gf)[1] = $(size(agf)[1]).") : new{typeof(gf)}(ts, Ts, Fpis, pos, gf, agf, T, Fpin, po)
end

"""
    Substance(name, CAS, Tchar, θchar, ΔCp, φ₀, ann, Dag, t₀, τ₀)

Structure to describe the properties of a solute, which migrates through the GC Column. These datas are in most cases read from a database with the function `load_solute_database()`.

# Arguments
* `name`: Name of the solute.
* `CAS`: CAS number of the solute.
* `Tchar`: Characterisic temperature (in K). One of the three distribution-centric thermodynamic parameters describing the retention of this solute on the given stationary phase.
* `θchar`: Characterisic parameters (in °C). One of the three distribution-centric thermodynamic parameters describing the retention of this solute on the given stationary phase.
* `ΔCp`: Change of the isobaric heat capacity moving from the mobile to the stationary phase (in J mol⁻¹ K⁻¹). One of the three distribution-centric thermodynamic parameters describing the retention of this solute on the given stationary phase.
* `φ₀`: Dimensionless film thickness (φ ≈ df/d) of the column for which the thermodynamic parameters (Tchar, θchar, ΔCp) were estimated.
* `ann`: Annotations. In most cases the source of the data is noted here.
* `Dag`: The diffusitivity of the solute `a` in the mobile phase `g` (in...). It is calculated by the function `diffusitivity()`.
* `t₀`: Initial time of the solute (in s) at the start of the simulation.
* `τ₀`: Initial peak width of the solute (in s) at the start of the simulation.

See also: [`load_solute_database`](@ref)
"""
struct Substance
    name::String        # name of solute
    CAS::String         # CAS registry number
    Tchar               # characteristic temperature in K
    θchar               # characteristic thermal constant in °C
    ΔCp                 # 3rd parameter
    φ₀                  # dimless film thickness for which Tchar, θchar and ΔCp were estimated
    ann::String         # annotations, e.g. the source of the data from which Tchar, θchar and ΔCp were estimated
    Dag                 # diffusion coefficient of analyt in a gas, calculate from structure (or from measurements)
    t₀                  # initial time in s
    τ₀                  # initial peak width in s
end

"""
    Options(alg, abstol, reltol, Tcontrol, odesys, ng, vis)

Structure describing some general options for the simulation.

# Arguments
* `alg`: The algorithm used for the ODE solver. The algorithms `OwrenZen3()`, `OwrenZen4()` and `OwrenZen5()` are recommended.
* `abstol`: The absolute tolerance for the ODE solver. Recommended value 1e-6 to 1e-8.
* `reltol`: The relative tolerance for the ODE solver. Recommended value 1e-3 to 1e-5.
* `Tcontrol`: Option defining at which point of the column the temperature program is calculated. The options are `inlet` (x=0) and `outlet` (x=L).
* `odesys`: Combine the ODEs for migration and peak-width into a system of ODEs (`odesys = true`) or solve the two ODEs separately (`odesys = false`).
* `ng`: Option to calculate the simulation without a gradient (`ng = true`) or with a gradient (`ng = false`). This distinction is made because of partly manuall differentiation (problem of automatic differentiation with integrals, e.g. in the `flow_restriction()` function. -> **TODO**: test package Quadrature.jl as alternative to QuadGK.jl for integration)
* `vis`: Used model of viscosity. `HP` is a second-order polynomial taken from the HP flow calculator. `Blumberg` is an emperical formula according to the book
    `Temperature-programmed Gas Chromatography` by Leonid M. Blumberg (2010, Wiley-VCH)
* `control`: Control of the "Flow" or of the "Pressure" (at column inlet) during the program

**TODO**: add option for the retention model ('ABC', 'K-centric')

For more informations about the arguments `alg`, `abstol` and `reltol` see the documentation of the DifferentialEquations.jl package.
"""
struct Options
    alg                 # algorithmen for the ODE solver
    abstol              # absolute tolerance for ODE solver
    reltol              # relative tolerance for ODE solver
    Tcontrol::String    # temperature control at 'inlet' (top) or 'outlet' (bottom) of the column
	odesys::Bool  		# calculate the two ODEs (migration and peak-width) separately (false) or
                        # combined as a system of ODEs (true)
    ng::Bool            # non-gradient calculation, ignores a defined spatial change of d, df or T
    vis::String         # viscosity model 'HP' or 'Blumberg'
    control::String     # control of the 'Flow' or of the inlet 'Pressure' during the program
end

"""
    Parameters(col, prog, sub, opt)

Structure describing all parameters for the simulation of a GC system.

# Arguments
* `col`: Structure `Column` describing the parameters of the GC column and
    mobile phase gas.
* `prog`: Structure `Program` describing the temperature and pressure
    program of a GC Column.
* `sub`: An array of the structure `Substance` describing the parameters of
    the solutes which are separated in the GC simulation.
* `opt`: Structure `Options` describing additional option parameters.
"""
struct Parameters
    col::Column
    prog::Program
    sub::Array{Substance,1}
    opt::Options
end
# ---End-Structures---

# ---Begin-Constructor-functions-for-Structures---

"""
    Column(L, d, df, sp, gas)

Construct the structure `Column` with given values for the case
of constant diameter `d` and film thickness `df`.

# Arguments
* `L`: Length of the capillary measured in m (meter).
* `d`: Diameter of the capillary measured in m (meter).
* `d_f`: Film thickness of the capillary measured in m (meter).
* `sp`: The name of the stationary phase.
* `gas`: The name of the mobile phase. Allowed values: He, H2 or N2.

# Examples
```julia
julia> Column(10.0, 0.1e-3, 0.1e-6, "DB5", "He")
	```
"""
function Column(L, d, df, sp, gas)
    # function to construct the Column structure
    # for the case of constant diameter and constant film thickness
    d_func(x) = gradient(x, [d])
    df_func(x) = gradient(x, [df])
    col = Column(L, d_func, [d], df_func, [df], sp, gas)
    return col
end

"""
    Program(time_steps, temp_steps, Fpin_steps, pout_steps, a_gf, Tcontrol, L)

Construct the structure `Program` with given values.

# Arguments
* `time_steps`: Time steps in s (seconds).
* `temp_steps`: Temperature steps in °C (degree celsius).
* `Fpin_steps`: Flow steps in m³/s resp. inlet pressure steps in Pa(a).
* `pout_steps`: Outlet pressure steps in Pa(a) (pascal, absolute).
* `a_gf`: Parameters of the gradient function.
* `Tcontrol`: Option defining at which point of the column the temperature
program is calculated. The options are `inlet` (x=0) and `outlet` (x=L).
* `L`: Length of the capillary measured in m (meter).

The length of the arrays `time_steps`, `temp_steps`, `pin_steps`, `pout_steps` and `a_gf`
have to be the same.

The arguments `Tcontrol` and `L` are used to construct the thermal gradient
function `gf(x)` and the temperature interpolation `T_itp(x,t)`.

# Examples
```julia
julia> Program([0.0, 60.0, 300.0, 120.0],
        [40.0, 40.0, 320.0, 320.0],
        300000.0.*ones(4),
        zeros(4),
        [[20.0, 20.0, 20.0, 20.0] zeros(4) 10.0.*ones(4) [0.0, -2.0, -5.0, -5.0]],
        "inlet",
        10.0)
```
"""
function Program(time_steps::Array{<:Real, 1}, temp_steps::Array{<:Real, 1}, Fpin_steps::Array{<:Real, 1}, pout_steps::Array{<:Real, 1}, a_gf::Array{<:Real, 2}, Tcontrol, L)
    # function to construct the Program structure
    # using as gradient function the exponential model 'gf_exp(x,a_gf,Tcontrol)'
    gf(x) = gradient(x, a_gf; Tcontrol=Tcontrol)
    T_itp = temperature_interpolation(time_steps, temp_steps, gf, L)
    Fpin_itp = steps_interpolation(time_steps, Fpin_steps)
    pout_itp = steps_interpolation(time_steps, pout_steps)
    prog = Program(time_steps, temp_steps, Fpin_steps, pout_steps, gf, a_gf, T_itp, Fpin_itp, pout_itp)
    return prog
end

"""
    Program(time_steps, temp_steps, Fpin_steps, pout_steps, ΔT_steps, x₀_steps, L₀_steps, α_steps, Tcontrol, L)

Construct the structure `Program` with given values.

# Arguments
* `time_steps`: Time steps in s (seconds).
* `temp_steps`: Temperature steps in °C (degree celsius).
* `Fpin_steps`: Flow steps in m³/s resp. inlet pressure steps in Pa(a).
* `pout_steps`: Outlet pressure steps in Pa(a) (pascal, absolute).
* `ΔT_steps`: Parameters of the gradient function. Temperature difference in °C.
* `x₀_steps`: Parameters of the gradient function. Spatial offset of the gradient in m.
* `L₀_steps`: Parameters of the gradient function. Distance over which the temperature difference in ΔT_steps is measured, in m.
* `α_steps`: Parameters of the gradient function. Factor describing the gradient profile.
* `Tcontrol`: Option defining at which point of the column the temperature
program is calculated. The options are `inlet` (x=0) and `outlet` (x=L).
* `L`: Length of the capillary measured in m (meter).

The length of the arrays `time_steps`, `temp_steps`, `Fpin_steps`, `pout_steps`, `ΔT_steps`, `x₀_steps`, `L₀_steps` and `α_steps`
have to be the same.

The arguments `Tcontrol` and `L` are used to construct the thermal gradient
function `gf(x)` and the temperature interpolation `T_itp(x,t)`.

# Examples
```julia
julia> Program([0.0, 60.0, 300.0, 120.0],
        [40.0, 40.0, 320.0, 320.0],
        300000.0.*ones(4),
        zeros(4),
        [20.0, 20.0, 20.0, 20.0],
        zeros(4),
        10.0.*ones(4),
        [0.0, -2.0, -5.0, -5.0],
        "inlet",
        10.0)
```
"""
function Program(time_steps::Array{<:Real, 1}, temp_steps::Array{<:Real, 1}, Fpin_steps::Array{<:Real, 1}, pout_steps::Array{<:Real, 1}, ΔT_steps::Array{<:Real, 1}, x₀_steps::Array{<:Real, 1}, L₀_steps::Array{<:Real, 1}, α_steps::Array{<:Real, 1}, Tcontrol, L)
    # function to construct the Program structure
    # using as gradient function the exponential model 'gf_exp(x,a_gf,Tcontrol)'
    a_gf = [ΔT_steps x₀_steps L₀_steps α_steps]
    gf(x) = gradient(x, a_gf; Tcontrol=Tcontrol)
    T_itp = temperature_interpolation(time_steps, temp_steps, gf, L)
    Fpin_itp = steps_interpolation(time_steps, Fpin_steps)
    pout_itp = steps_interpolation(time_steps, pout_steps)
    prog = Program(time_steps, temp_steps, Fpin_steps, pout_steps, gf, a_gf, T_itp, Fpin_itp, pout_itp)
    return prog
end

"""
    Program(TP, FpinP, poutP, ΔTP, x₀P, L₀P, αP, Tcontrol, L; time_unit="min")

Construct the structure `Program` with conventional formulation (see [`conventional_program`](@ref)) of programs.

# Arguments
* `TP`: conventional formulation of a temperature program.
* `FpinP`: conventional formulation of a Flow (in m³/s) resp. inlet pressure (in Pa(a)) program.
* `poutP`: conventional formulation of a outlet pressure (in Pa(a)) program.
* `ΔTP`: conventional formulation of a program for temperature gradient parameter ΔT.
* `x₀P`: conventional formulation of a program for temperature gradient parameter x₀.
* `L₀P`: conventional formulation of a program for temperature gradient parameter L₀.
* `αP`: conventional formulation of a program for temperature gradient parameter αP.
* `Tcontrol`: Option defining at which point of the column the temperature
program is calculated. The options are `inlet` (x=0) and `outlet` (x=L).
* `L`: Length of the capillary measured in m (meter).
* `time_unit`: unit of time in the programs, "min"` (default) times are measured in minutes, "s" times are measured in seconds.

The arguments `Tcontrol` and `L` are used to construct the thermal gradient
function `gf(x)` and the temperature interpolation `T_itp(x,t)`.

"""
function Program(TP::Array{<:Real, 1}, FpinP::Array{<:Real, 1}, poutP::Array{<:Real, 1}, ΔTP::Array{<:Real, 1}, x₀P::Array{<:Real, 1}, L₀P::Array{<:Real, 1}, αP::Array{<:Real, 1}, Tcontrol::String, L::Float64; time_unit="min")
    # using as gradient function the exponential model 'gf_exp(x,a_gf,Tcontrol)
    ts = Array{Array{Float64,1}}(undef, 7)
    ts[1], Ts = GasChromatographySimulator.conventional_program(TP; time_unit=time_unit)
    ts[2], Fps = GasChromatographySimulator.conventional_program(FpinP; time_unit=time_unit)
    ts[3], pouts = GasChromatographySimulator.conventional_program(poutP; time_unit=time_unit)
    ts[4], ΔTs = GasChromatographySimulator.conventional_program(ΔTP; time_unit=time_unit)
    ts[5], x₀s = GasChromatographySimulator.conventional_program(x₀P; time_unit=time_unit)
    ts[6], L₀s = GasChromatographySimulator.conventional_program(L₀P; time_unit=time_unit)
    ts[7], αs = GasChromatographySimulator.conventional_program(αP; time_unit=time_unit)
    time_steps = Float64[]
    for i=1:7
        time_steps = GasChromatographySimulator.common_time_steps(time_steps, ts[i])
    end
    temp_steps = GasChromatographySimulator.new_value_steps(Ts, ts[1], time_steps)
    Fpin_steps = GasChromatographySimulator.new_value_steps(Fps, ts[2], time_steps)
    pout_steps = GasChromatographySimulator.new_value_steps(pouts, ts[3], time_steps)
    ΔT_steps = GasChromatographySimulator.new_value_steps(ΔTs, ts[4], time_steps)
    x₀_steps = GasChromatographySimulator.new_value_steps(x₀s, ts[5], time_steps)
    L₀_steps = GasChromatographySimulator.new_value_steps(L₀s, ts[6], time_steps)
    α_steps = GasChromatographySimulator.new_value_steps(αs, ts[7], time_steps)

    a_gf = [ΔT_steps x₀_steps L₀_steps α_steps]
    gf(x) = gradient(x, a_gf; Tcontrol=Tcontrol)
    T_itp = temperature_interpolation(time_steps, temp_steps, gf, L)
    Fpin_itp = steps_interpolation(time_steps, Fpin_steps)
    pout_itp = steps_interpolation(time_steps, pout_steps)
    prog = Program(time_steps, temp_steps, Fpin_steps, pout_steps, gf, a_gf, T_itp, Fpin_itp, pout_itp)
    return prog
end

"""
    Program(time_steps::Array{<:Real, 1}, temp_steps::Array{<:Real, 1}, Fpin_steps::Array{<:Real, 1}, pout_steps::Array{<:Real, 1}, L)

Construct the structure `Program` with given values for the case
without a thermal gradient.

# Arguments
* `time_steps::Array{<:Real, 1}`: Time steps in s (seconds).
* `temp_steps::Array{<:Real, 1}`: Temperature steps in °C (degree celsius).
* `Fpin_steps::Array{<:Real, 1}`: Flow steps in m³/s resp. inlet pressure steps in Pa(a).
* `pout_steps::Array{<:Real, 1}`: Outlet pressure steps in Pa(a) (pascal, absolute).
* `L`: Length of the capillary measured in m (meter).

The length of the arrays `time_steps`, `temp_steps`, `Fpin_steps` and `pout_steps`
have to be the same.

The argument `L` is used to construct the temperature interpolation `T_itp(x,t)`.

# Examples
```julia
julia> Program([0.0, 60.0, 300.0, 120.0],
        [40.0, 40.0, 320.0, 320.0],
        300000.0.*ones(4),
        zeros(4),
        10.0)
```
"""
function Program(time_steps::Array{<:Real, 1}, temp_steps::Array{<:Real, 1}, Fpin_steps::Array{<:Real, 1}, pout_steps::Array{<:Real, 1}, L)
    # function to construct the Program structure
    # without a thermal gradient
    # using as gradient function the exponential model 'gf_exp(x,a_gf,Tcontrol)'
    a_gf = [zeros(length(time_steps)) zeros(length(time_steps)) L.*ones(length(time_steps)) zeros(length(time_steps))]
    gf(x) = gradient(x, a_gf)
    T_itp = temperature_interpolation(time_steps, temp_steps, gf, L)
    Fpin_itp = steps_interpolation(time_steps, Fpin_steps)
    pout_itp = steps_interpolation(time_steps, pout_steps)
    prog = Program(time_steps, temp_steps, Fpin_steps, pout_steps, gf, a_gf, T_itp, Fpin_itp, pout_itp)
    return prog
end

"""
    Program(TP, FpinP, L; pout="vacuum", time_unit="min")

Construct the structure `Program` with conventional formulation (see [`conventional_program`](@ref)) of programs for the case
without a thermal gradient.

# Arguments
* `TP`: conventional formulation of a temperature program.
* `FpinP`: conventional formulation of a Flow (in m³/s) resp. inlet pressure (in Pa(a)) program.
* `L`: Length of the capillary measured in m (meter).
* `pout`: Outlet pressure, "vacuum" (default) or "atmosphere".
* `time_unit`: unit of time in the programs, "min"` (default) times are measured in minutes, "s" times are measured in seconds.

The argument `L` is used to construct the temperature interpolation `T_itp(x,t)`.

# Examples
```julia
julia> Program((40.0, 1.0, 5.0, 280.0, 2.0, 20.0, 320.0, 2.0),
                (400000.0, 10.0, 5000.0, 500000.0, 20.0),
                10.0)
```
"""
function Program(TP, FpinP, L; pout="vacuum", time_unit="min")
    ts1, Ts = GasChromatographySimulator.conventional_program(TP; time_unit=time_unit)
    ts2, Fps = GasChromatographySimulator.conventional_program(FpinP; time_unit=time_unit)
    time_steps = GasChromatographySimulator.common_time_steps(ts1, ts2)
    temp_steps = GasChromatographySimulator.new_value_steps(Ts, ts1, time_steps)
    Fpin_steps = GasChromatographySimulator.new_value_steps(Fps, ts2, time_steps)
    if pout == "vacuum"
        pout_steps = zeros(length(time_steps))
    else
        pout_steps = pn.*ones(length(time_steps))
    end
    a_gf = [zeros(length(time_steps)) zeros(length(time_steps)) L.*ones(length(time_steps)) zeros(length(time_steps))]
    gf(x) = gradient(x, a_gf)
    T_itp = temperature_interpolation(time_steps, temp_steps, gf, L)
    Fpin_itp = steps_interpolation(time_steps, Fpin_steps)
    pout_itp = steps_interpolation(time_steps, pout_steps)
    prog = Program(time_steps, temp_steps, Fpin_steps, pout_steps, gf, a_gf, T_itp, Fpin_itp, pout_itp)
    return prog
end

"""
    Options(;alg=OwrenZen5(), abstol=1e-6, reltol=1e-3, Tcontrol="inlet", odesys=true, ng=false, vis="Blumberg", control="Pressure")

Construct the structure `Options` with default values.

# Arguments
* `alg`: The algorithm used for the ODE solver. The algorithms
    `OwrenZen3()`, `OwrenZen4()` and `OwrenZen5()` are recommended.
* `abstol`: The absolute tolerance for the ODE solver. Recommended value
    1e-6 to 1e-8.
* `reltol`: The relative tolerance for the ODE solver. Recommended value
1e-3 to 1e-5.
* `Tcontrol`: Option defining at which point of the column the temperature
    program is calculated. The options are `inlet` (x=0) and `outlet` (x=L).
* `odesys`: Combine the ODEs for migration and peak-width into a system of
    ODEs (`odesys = true`) or solve the two ODEs separately (`odesys = false`).
* `ng`: Option to calculate the simulation without a gradient (`ng = true`)
    or with a gradient (`ng = false`). This distinction is made because of
    partly manuall differentiation (problem of automatic differentiation with
    integrals, e.g. in the `flow_restriction()` function. -> TODO: test
    package Quadrature.jl as alternative to QuadGK.jl for integration)
* `vis`: Used model of viscosity. `HP` is a second-order polynomial taken from the HP flow calculator. `Blumberg` is an emperical formula according to the book
    `Temperature-programmed Gas Chromatography` by Leonid M. Blumberg (2010, Wiley-VCH)
* `control`: Control of the "Flow" or of the "Pressure" (at column inlet) during the program

For more informations about the arguments `alg`, `abstol` and `reltol` see
the documentation of the DifferentialEquations.jl package.

# Examples
```julia
julia> Options()
```

```julia
julia> Options(abstol=1e-8, Tcontrol="outlet")
```
"""
function Options(;alg=OwrenZen5(), abstol=1e-6, reltol=1e-3, Tcontrol="inlet", odesys=true, ng=false, vis="Blumberg", control="Pressure")
    opt = Options(alg, abstol, reltol, Tcontrol, odesys, ng, vis, control)
    return opt
end

"""
    Options(alg, abstol, reltol, Tcontrol, odesys; ng=false, vis="Blumberg", control="Pressure")

Construct the structure `Options` with given values.

# Arguments
* `alg`: The algorithm used for the ODE solver. The algorithms
    `OwrenZen3()`, `OwrenZen4()` and `OwrenZen5()` are recommended.
* `abstol`: The absolute tolerance for the ODE solver. Recommended value
    1e-6 to 1e-8.
* `reltol`: The relative tolerance for the ODE solver. Recommended value
1e-3 to 1e-5.
* `Tcontrol`: Option defining at which point of the column the temperature
    program is calculated. The options are `inlet` (x=0) and `outlet` (x=L).
* `odesys`: Combine the ODEs for migration and peak-width into a system of
    ODEs (`odesys = true`) or solve the two ODEs separately (`odesys = false`).
* `ng`: Option to calculate the simulation without a gradient (`ng = true`)
    or with a gradient (`ng = false`). This distinction is made because of
    partly manuall differentiation (problem of automatic differentiation with
    integrals, e.g. in the `flow_restriction()` function. -> TODO: test
    package Quadrature.jl as alternative to QuadGK.jl for integration)
* `vis`: Used model of viscosity. `HP` is a second-order polynomial taken from the HP flow calculator. `Blumberg` is an emperical formula according to the book
    `Temperature-programmed Gas Chromatography` by Leonid M. Blumberg (2010, Wiley-VCH)
* `control`: Control of the "Flow" or of the "Pressure" (at column inlet) during the program

For more informations about the arguments `alg`, `abstol` and `reltol` see
the documentation of the DifferentialEquations.jl package.

# Examples
```julia
julia> Options(OwrenZen3(), 1e-7, 1e-4, "inlet", true)
```

```julia
julia> Options(OwrenZen3(), 1e-7, 1e-4, "inlet", true; ng=true, vis="HP", control="Flow")
```
"""
function Options(alg, abstol, reltol, Tcontrol, odesys; ng=false, vis="Blumberg", control="Pressure")
    opt = Options(alg, abstol, reltol, Tcontrol, odesys, ng, vis, control)
    return opt
end

# Aliases for compatibility
constructor_System = Column
constructor_Program = Program
constructor_Options = Options
#---End-Constructor-functions-for-Structures---

#---Begin-Functions-used-for-Parameter-construction
"""
    gradient(x, a; Tcontrol="inlet")

Defines the gradient of the column diameter, film thickness or
temperature along the GC column.

# Arguments
* `x`: Position along the GC column, in m.
* `a`: Parameters of the gradient function

The form of `a` decides the actual used function for the gradient. The
following options are realized:
* `a` is a single value (e.g. Float or Int): The gradient function is constant for all positions `x` with the value of `a`.
* `a` is a 1D-array of length = 4: The gradient function is a exponential function. The 4 values of `a`: `a = [f₀, x₀, L₀, α]` with
    * `f₀`: Start value at `x = x₀`.
    * `x₀`: Shift in `x` position.
    * `L₀`: Distance over which the value drops from f₀ to 0.
    * `α`: Factor describing the gradient profile.
    if `α <= 0`
        `f = f₀ .* (1 .- exp.(α.*(1 .- (x.-x₀) ./ L₀)) + (1 .- (x.-x₀) ./ L₀)
        .* exp.(α))`
    if `α > 0`
        `f = f₀ .* (exp.(-α.*(x.-x₀) ./ L₀) .- (x.-x₀) ./ L₀ .* exp.(-α))`
* `a` is a 2D-array with size = (n, 4): The gradient function is an
    exponential function. The 4 entrys have the same meaning as above, but their
    values can change over the times defined in the `time_steps` of the
    Program structure. At these `time_steps[i]` the gradient function is
    described by the corresponding parameters `a[i,:]`. Between the
    `time_steps[i]` and `time_steps[j]` the value of the gradient function at position `x` is linear
    interpolated from the gradien functions defined by `a[i,:]` and `a[j,:]`.

# Examples
```julia
julia> d(x) = gradient(x, 0.1e-3)
```

```julia
julia> gf(x) = gradient(x, [[20.0, 20.0, 20.0, 20.0] zeros(4) 10.0.*ones(4) [0.0, -2.0, -5.0, -5.0]])
```
"""
# gradient functions
function gradient(x, a; Tcontrol="inlet")
    if size(a) == (1,)
        # constant value, no gradient
        f = a[1]
    elseif size(a) > (1,)
        if  length(size(a))==1
            if length(a)==4
                # for diameter or film thickness, values of parameters 'a' are fixed
                # over time
                f₀ = a[1] # start value
                x₀ = a[2] # shift in x
                L₀ = a[3] # length of the linear segment
                α  = a[4] # exponential factor, α=0 -> linear
                if α<=0 # concave/linear, weak change at beginning and strong change at end of column
                    f = f₀ .* (1 .- exp.(α.*(1 .- (x.-x₀) ./ L₀)) + (1 .- (x.-x₀) ./ L₀) .* exp.(α))
                elseif α>0 # convex function, strong change at beginning and weak change at end of the column
                    f = f₀ .* (exp.(-α.*(x.-x₀) ./ L₀) .- (x.-x₀) ./ L₀ .* exp.(-α))
                end
            # other functions ...
            end
        elseif length(size(a)) == 2
            if a==zeros(size(a))
                f = zeros(size(a)[1])
                #if size(a)[2] == 1
                # for thermal gradient, no change in time  of values of the
                # parameters
            elseif size(a)[2] == 4
                # for thermal gradient, values of parameters 'a[i,:]' can change
                # over time
                f₀ = a[:,1] # start value
                x₀ = a[:,2] # shift in x (e.g. to correct for the real position of IR-sensors)
                L₀ = a[:,3] # length of the linear segment
                α  = a[:,4] # exponential factor, α=0 -> linear
                f = Array{Float64}(undef, length(α))
                for i=1:length(α)
                    if α[i]<=0 # concave function, weak change at beginning and strong change at end of column
                        if Tcontrol=="inlet"
                            f[i] = f₀[i] .* (.- exp.(α[i].*(1 .- (x.-x₀[i]) ./ L₀[i])) + (1 .- (x.-x₀[i]) ./ L₀[i]) .* exp.(α[i]))
                        elseif Tcontrol=="outlet"
                            f[i] = f₀[i] .* (1 .- exp.(α[i].*(1 .- (x.-x₀[i]) ./ L₀[i])) + (1 .- (x.-x₀[i]) ./ L₀[i]) .* exp.(α[i]))
                        end
                    elseif α[i]>0 # convex function, strong change at beginning and weak change at end of the column
                        if Tcontrol=="inlet"
                            f[i] = f₀[i] .* (exp.(-α[i].*(x.-x₀[i]) ./ L₀[i]) .- (x.-x₀[i]) ./ L₀[i] .* exp.(-α[i]) .- 1)
                        elseif Tcontrol=="outlet"
                            f[i] = f₀[i] .* (exp.(-α[i].*(x.-x₀[i]) ./ L₀[i]) .- (x.-x₀[i]) ./ L₀[i] .* exp.(-α[i]))
                        end
                    end
                end
            # other functions ...
            end
        end
    end
    return f
end

"""
    temperature_interpolation(time_steps, temp_steps, gradient_function, L)

Construct the temperature function depending on position `x` and
time `t`.

# Arguments
* `time_steps::Array{<:Real,1}`: Time steps in s (seconds).
* `temp_steps::Array{<:Real,1}`: Temperature steps in °C (degree celsius).
* `gf::Function`: Gradient function `gf(x, a_gf)`, describes the thermal gradient.
* `L::Float64`: Length of the capillary measured in m (meter).

For the spatial dependency of the interpolated temperature `T_ipt(x,t)` the
gradient function `gf` is calculated every 1e-3 m (1 mm). Positions
inbetween are linear interpolated. For the temporal dependency the
temperatures `temp_steps` defined at the `time_steps` are linear
interpolated over time `t`.

# Examples
```julia
julia> T_itp = temperature_interpolation([0.0, 60.0, 300.0, 120.0], [40.0, 40.0, 320.0, 320.0], gf, 10.0)
```
"""
function temperature_interpolation(time_steps::Array{<:Real,1}, temp_steps::Array{<:Real,1}, gradient_function::Function, L)
	T(x) = temp_steps .+ gradient_function(x)
	nx = 0.0:1e-3:L # mm exact
	nt = cumsum(time_steps)
	Tmat = Array{Float64}(undef, length(nx), length(nt))
	for j=1:length(nt)
		for i=1:length(nx)
			Tmat[i,j] = T(nx[i])[j] + 273.15
		end
	end
	T_itp = LinearInterpolation((nx, nt), Tmat, extrapolation_bc=Flat())
	return T_itp
end

"""
    step_interpolation(time_steps, steps)

Construct a linear interpolated function depending on time `t` of the `steps`-values over `time_steps`.

# Arguments
* `time_steps::Array{<:Real,1}`: Time steps in s (seconds).
* `steps::Array{<:Real,1}`: steps, e.g. pressure or flow.

# Examples
```julia
julia> pin_itp = steps_interpolation([0.0, 60.0, 300.0, 120.0],
                                    [300000.0, 300000.0, 400000.0, 400000.0])
```
"""
function steps_interpolation(time_steps::Array{<:Real,1}, steps::Array{<:Real,1})
    s_itp = LinearInterpolation((cumsum(time_steps), ), steps, extrapolation_bc=Flat())
    return s_itp
end

"""
    load_solute_database(db, sp, gas, solutes, t₀, τ₀)

Load the data of `solutes` for the stationary phase `sp` and the mobile
phase `gas` from the database `db` into an array of the structure `Substance`.

# Arguments
* `db::DataFrame`: DataFrame of the database.
* `sp::String`: Name of the stationary phase.
* `gas::String`: Name of the mobile phase.
* `solutes::Array{<:AbstractString,1}`: Name of the solutes.
* `t₀::Array{Float64,1}`: Initial start times of the solutes.
* `τ₀::Array{Float64,1}`: Initial peak widths of the solutes.

# Examples
```julia
julia> sub = load_solute_database(db, "DB5", "He", ["C10", "C11"], [0.0, 0.0], [0.5, 0.5])
```
"""
function load_solute_database(db::DataFrame, sp::String, gas::String, solutes::Array{<:AbstractString,1}, t₀::Array{Float64,1}, τ₀::Array{Float64,1})
	# load the information about a solute from a data base and collecting these informations in the
    # structure SubstanceGC2
    #
	# db ... Dataframe of the database
	# sp ... Name of the stationary phase
	# gas ... Name of the used mobile phase gas
	# solutes ... Names of the solutes for which the informations should be used
    # τ₀ ... initial values of the peak width for the solutes
    # t₀ ... initial values of the time for the solutes
	#
	# at the moment only the new database config is supported (new data is appended in new rows)
	# also, if for a solute no data is available in the database no error is given!!!
	if size(db)[2]>14
		# old database config (additional stationary phases are appended in new columns, only one row for a solute)
        error("Data format not supported. Use the appended database structure.")
	else
		# new database config (simpler, additional stationary phases are appended in new rows, multiple entrys for a solute)
		# 1. Filter the stationary phase
		phases_db = unique(db.Phase)
		if isa(findfirst(phases_db.==sp), Nothing) && sp!=""
			error("Unknown selction of stationary phase. Choose one of these: $phases_db")
		elseif sp=="" # no stationary phase is selected, like for transferlines
			db_filtered = db[!,1:8] # use only the data of the first 8 columns
			# 2. Filter the solutes, not using multiple entrys
			db_filtered_1=unique(db_filtered[in(solutes).(db_filtered.Name),:])
			# use placeholder values
			Tchar = ones(length(solutes))
			θchar = ones(length(solutes))
			ΔCp = ones(length(solutes))
			φ₀ = ones(length(solutes))
			Annotation = fill("no sp", length(solutes))
		else
			db_filtered = filter([:Phase] => x -> x==sp, db)
			# 2. Filter the solutes
			db_filtered_1=db_filtered[in(solutes).(db_filtered.Name),:]
			# values
			Tchar = db_filtered_1.Tchar.+Tst
			θchar = db_filtered_1.thetachar
			ΔCp = db_filtered_1.DeltaCp
			φ₀ = db_filtered_1.phi0
			Annotation = db_filtered_1.Annotation
		end
        sub = Array{Substance}(undef, length(solutes))
        for i=1:length(solutes)
			Dag = diffusivity(db_filtered_1.Molmass[i],
											db_filtered_1.Cnumber[i],
											db_filtered_1.Hnumber[i],
											db_filtered_1.Onumber[i],
											db_filtered_1.Nnumber[i],
											db_filtered_1.Ringnumber[i],
											gas)
			sub[i] = Substance(db_filtered_1.Name[i],
										db_filtered_1.CAS[i],
										Tchar[i],
										θchar[i],
										ΔCp[i],
										φ₀[i],
										Annotation[i],
										Dag,
										t₀[i],
										τ₀[i])
		end
	end
	return sub
end

"""
    load_solute_database(db_path, db, sp, gas, solutes, t₀, τ₀)

Load the data of `solutes` for the stationary phase `sp` and the mobile
phase `gas` from the database file `db` (located in `db_path`) into an array
of the structure `Substance`.

# Arguments
* `db_path::String`: Path to the database file.
* `db::String`: Name of the database file.
* `sp::String`: Name of the stationary phase.
* `gas::String`: Name of the mobile phase.
* `solutes::Array{<:AbstractString,1}`: Name of the solutes.
* `t₀::Array{Float64,1}`: Initial start times of the solutes.
* `τ₀::Array{Float64,1}`: Initial peak widths of the solutes.

# Examples
```julia
julia> sub = load_solute_database("path/to/the/file", "db.csv", "DB5", "He", ["C10", "C11"], [0.0, 0.0], [0.5, 0.5])
```
"""
function load_solute_database(db_path::String, db::String, sp::String, gas::String, solutes::Array{<:AbstractString,1}, t₀::Array{Float64,1}, τ₀::Array{Float64,1})
	# load the information about a solute from a data base and collecting these informations in the
    # structure Substance
    #
	# db_path ... Path to the database file
	# db ... Name of the database
	# sp ... Name of the stationary phase
	# gas ... Name of the used mobile phase gas
	# solutes ... Names of the solutes for which the informations should be used
    # τ₀ ... initial values of the peak width for the solutes
    # t₀ ... initial values of the time for the solutes
	db_dataframe = DataFrame(CSV.File(string(db_path,"/",db), header=1, silencewarnings=true))
    sub = load_solute_database(db_dataframe, sp, gas, solutes, t₀, τ₀)
	return sub
end

"""
    all_solutes(sp, db)

Extract the name of all solutes for which data in a database `db` and the
stationay phase `sp` is available.

# Arguments
* `sp`: Name of the stationary phase.
* `db`: DataFrame of the database.

# Examples
```julia
julia> all = all_solutes("DB5", db)
```
"""
function all_solutes(sp, db)
	db_filter = filter([:Phase] => x -> x==sp, db)
	solutes = string.(db_filter.Name)
	return solutes
end

"""
    conventional_program(CP; time_unit="min")

Translate the conventional program notation into a vector of time steps and value steps (temperature, pressure, flow) used in GasChromatographySimulator.Program

The conventional temperature program is defined as an array of the following form (for a temperature program):
`CP = [T₁, t₁, r₁, T₂, t₂, r₂, T₃, t₃, ...]` corresponding to the notation:
`T₁(t₁) - r₁ - T₂(t₂) - r₂ - T₃(t₃) - ...` which can be read as:
Starting temperature `T₁` is holded for time `t₁`. After the holding time the temperature increases linearly with the heating rate `r₁`, until temperature `T₂` is reached. This temperature is held for the time `t₂` after which the temperature increases linearly by the heating rate `r₂` until temperature `T₃` is reached, which is hold for the time `t₃`, and so on.

The option `time_unit` determines the unit of time in the program `CP`. For `time_unit = "min"` (default) the times are measured in minutes and the heating rates in °C/min. For `time_unit = "s"` the times are measured in seconds and the heating rate in °C/s.
"""
function conventional_program(CP; time_unit="min")
    if time_unit == "min"
        c = 60.0
    elseif time_unit == "s"
        c = 1.0
    end
    value_steps = Float64[]
    for i=1:3:length(CP) # every third CP-entry starting from 1 is a value step
        push!(value_steps, CP[i])
        push!(value_steps, CP[i])
    end
    hold_times = Float64[]
    for i=2:3:length(CP)
        push!(hold_times, CP[i])
    end
    heating_rates = Float64[]
    for i=3:3:length(CP)
        push!(heating_rates, CP[i])
    end
    time_steps = Array{Float64}(undef, length(value_steps))
    time_steps[1] = 0.0
    for i=2:2:length(value_steps)
        time_steps[i] = hold_times[Int(i/2)] * c
        if i<length(value_steps)
            time_steps[i+1] = (value_steps[i+1] - value_steps[i])/heating_rates[Int(i/2)] * c
        end
    end
    return time_steps, value_steps
end

"""
    temperature_program(time_steps, value_steps; time_unit="min")

Translate the vector of time steps and value steps (temperature, pressure, flow) into a conventional program notation.

The conventional temperature program is defined as an array of the following form (for a temperature program):
`CP = [T₁, t₁, r₁, T₂, t₂, r₂, T₃, t₃, ...]` corresponding to the notation:
`T₁(t₁) - r₁ - T₂(t₂) - r₂ - T₃(t₃) - ...` which can be read as:
Starting temperature `T₁` is holded for time `t₁`. After the holding time the temperature increases linearly with the heating rate `r₁`, until temperature `T₂` is reached. This temperature is held for the time `t₂` after which the temperature increases linearly by the heating rate `r₂` until temperature `T₃` is reached, which is hold for the time `t₃`, and so on.

The option `time_unit` determines the unit of time in the program `CP`. For `time_unit = "min"` (default) the times are measured in minutes and the heating rates in °C/min. For `time_unit = "s"` the times are measured in seconds and the heating rate in °C/s.
"""
function temperature_program(time_steps, value_steps; time_unit="min")
    if time_unit == "min"
        c = 60.0
    elseif time_unit == "s"
        c = 1.0
    end
    # identify temperature pairs (the same value of temperature at neigboring elements)
    index_pair = Int[]
    for i=1:(length(value_steps)-1)
        if value_steps[i] == value_steps[i+1]
            push!(index_pair, i)
        end
    end
    # identify single temperatures
    index_single = Int[]
    for i=1:length(value_steps)
        if (i in index_pair) == false && (i in index_pair.+1) == false
            push!(index_single, i)
        end
    end
    values = value_steps[sort([index_pair; index_single])]
    # every (1+(i-1)*3)th element of VP is a value element of `values`

    thold = Array{Float64}(undef, length(values))
    a = sort([index_pair.+1; index_single])
    for i=1:length(values)
        if a[i] in index_single # holding time is zero for single entrys
            thold[i] = 0.0
        else # holding times are for paired temperatures the time_steps[index_pair.+1]
            thold[i] = time_steps[a[i]] / c
        end
    end
    # every (2+(i-1)*3)th element of VP is a holding time

    rate = Array{Float64}(undef, length(values)-1)
    theat = time_steps[sort([index_pair; index_single])[2:end]]
    for i=1:(length(values)-1)
        rate[i] = (values[i+1] - values[i])/theat[i] * c
    end

    VP = Array{Float64}(undef, 2 + (length(values)-1)*3)
    for i=1:length(values)
        VP[1+(i-1)*3] = values[i]
    end
    for i=1:length(thold)
        VP[2+(i-1)*3] = thold[i]
    end
    for i=1:length(rate)
        VP[3+(i-1)*3] = rate[i]
    end
    return VP
end

"""
    common_time_steps(time_steps_1, time_steps_2)

Estimate a new set of time steps, which represents the combination of `time_steps_1` and `time_steps_2`.
"""
function common_time_steps(time_steps_1, time_steps_2)
	# constructs the new timesteps common for all moduls
	ctselements = Float64[]
    for j=1:length(time_steps_1)
        push!(ctselements, cumsum(time_steps_1)[j])
    end
    for j=1:length(time_steps_2)
        push!(ctselements, cumsum(time_steps_2)[j])
    end
	new_time_steps = [0.0; diff(sort(unique(ctselements)))]
	return new_time_steps
end

"""
    new_value_steps(value_steps, time_steps, new_time_steps)

Estimate the new value steps at the `new_time_steps` from the original set of `value_steps` over `time_steps`. The new values at new time steps are calculated from a linear change of the value between the original time steps.
"""
function new_value_steps(value_steps, time_steps, new_time_steps)
    new_values = Array{Float64}(undef, length(new_time_steps))
    index_calc = Int[]
    index_old = Int[]
    for i=1:length(new_time_steps)
        if cumsum(new_time_steps)[i] in cumsum(time_steps)
            j = findfirst(cumsum(new_time_steps)[i].==cumsum(time_steps))
            new_values[i] = value_steps[j]
            push!(index_old, i)
        else
            push!(index_calc, i)
        end
    end
    for i=1:length(index_calc)
        i1 = findlast(index_old.<index_calc[i])
        i2 = findfirst(index_old.>index_calc[i])
        if isnothing(i1) || isnothing(i2)
            new_values[index_calc[i]] = value_steps[end]
        else
            if time_steps[i2] == 0.0
                i2 = i2 + 1
            end
            v1 = value_steps[i1]
            v2 = value_steps[i2]
            t1 = cumsum(time_steps)[i1]
            t2 = cumsum(time_steps)[i2]
            rate = (v2 - v1)/(t2 - t1)
            if t1 == t2
                new_values[index_calc[i]] = v2
            else
                new_values[index_calc[i]] = v1 + rate * (cumsum(new_time_steps)[index_calc[i]] - t1)
            end
        end
    end
    return new_values
end

#---End-Functions-used-for-Parameter-construction---

#---Begin-Solving-Functions---
"""
    simulate(par::Parameters)

Simulate the GC system defined by the structure `par`.

Note: Based on the option for `odesys` the result is different. For `odesys =
true` the result is a dataframe (the peaklist) and the solution of the ODEs
as a system (solution structure from DifferentialEquations.jl). If `odesys =
false` the result is a dataframe (the peaklist) and the two solutions of the
ODEs for time ``t(z)`` and peak variance ``τ²(z)``.
"""
function simulate(par)
    if par.opt.odesys==true
        sol = solve_system_multithreads(par)
    	pl = GasChromatographySimulator.peaklist(sol, par)
        return pl, sol
	else
		sol, peak = solve_multithreads(par)
    	pl = GasChromatographySimulator.peaklist(sol, peak, par)
        return pl, sol, peak
	end
end

"""
    solve_system_multithreads(par::Parameters)

Simulate the GC system defined by the structure `par` by solving the
ODEs for ``t(z)`` and ``τ²(z)`` together as a system of ODEs using multiple
threads (parallel computing) for the simulation of different solutes.

Note: The result is an array of the solution structure from DifferentialEquations.jl.

# Examples

```julia
julia> sol = solve_system_multithreads(par)
```
"""
function solve_system_multithreads(par)
	n = length(par.sub)
	sol = Array{Any}(undef, n)
	Threads.@threads for i=1:n
		sol[i] = solving_odesystem_r(par.col, par.prog, par.sub[i], par.opt)
	end
	return sol
end

"""
    solve_multithreads(par::Parameters)

Simulate the GC system defined by the structure `par` by solving the
ODEs for ``t(z)`` and ``τ²(z)`` separatly (solving ``t(z)`` and using this result
to solve for ``τ²(z)``) using multiple threads (parallel computing) for the
simulation of different solutes.

Note: The result are two arrays of the solution structure from
DifferentialEquations.jl.

# Examples

```julia
julia> sol, peak = solve_multithreads(par)
```
"""
function solve_multithreads(par)
    n = length(par.sub)
    sol = Array{Any}(undef, n)
    peak = Array{Any}(undef, n)
    Threads.@threads for i=1:n
        sol[i] = solving_migration(par.col, par.prog, par.sub[i], par.opt)
        peak[i] = solving_peakvariance(sol[i], par.col, par.prog, par.sub[i], par.opt)
    end
    return sol, peak
end

"""
    solving_migration(col::Column, prog::Program, sub::Substance, opt::Options)

Solve for the migration ``t(z)`` of solute `sub` in the GC Column `col` with
the program `prog` and the options `opt`.

Note: The result is the solution structure from
DifferentialEquations.jl.
"""
function solving_migration(col::Column, prog::Program, sub::Substance, opt::Options)
	f_tz(t,p,z) = residency(z, t, prog.T_itp, prog.Fpin_itp, prog.pout_itp, col.L, col.d, col.df, col.gas, sub.Tchar, sub.θchar, sub.ΔCp, sub.φ₀; ng=opt.ng, vis=opt.vis, control=opt.control)
    t₀ = sub.t₀
    zspan = (0.0,col.L)
    prob_tz = ODEProblem(f_tz, t₀, zspan)
    solution_tz = solve(prob_tz, alg=opt.alg, abstol=opt.abstol,reltol=opt.reltol)
    return solution_tz
end

"""
    solving_peakvariance(solution_tz, col::Column, prog::Program, sub::Substance, opt::Options)

Solve for the development of the peak variance ``τ²(z)`` of solute `sub` in the GC Column `col` with
the program `prog` and the options `opt` during its migration defined by `solution_tz`.

Note: The result is the solution structure from
DifferentialEquations.jl.
"""
function solving_peakvariance(solution_tz, col, prog, sub, opt)
    t(z) = solution_tz(z)
    p = [col, prog, sub, opt]
    f_τ²z(τ²,p,z) = peakode(z, t(z), τ², col, prog, sub, opt)
    τ²₀ = sub.τ₀^2
    zspan = (0.0, col.L)
    prob_τ²z = ODEProblem(f_τ²z, τ²₀, zspan, p)
    solution_τ²z = solve(prob_τ²z, alg=opt.alg, abstol=opt.abstol,reltol=opt.reltol)
    return solution_τ²z
end

"""
    solving_odesystem_r(col::Column, prog::Program, sub::Substance, opt::Options)

Solve the migration ``t(z)`` and peak variance development ``τ²(z)`` of solute `sub` in the GC Column `col` with
the program `prog` and the options `opt` as a system of ODEs.

Note: The result is the solution structure from
DifferentialEquations.jl.

See also: [`odesystem_r!`](@ref)
"""
function solving_odesystem_r(col::Column, prog::Program, sub::Substance, opt::Options)
    t₀ = [sub.t₀; sub.τ₀^2]
    zspan = (0.0,col.L)
	p = [col, prog, sub, opt]
    prob = ODEProblem(odesystem_r!, t₀, zspan, p)

    solution = solve(prob, alg=opt.alg, abstol=opt.abstol,reltol=opt.reltol)

    if solution.t[end]<col.L
        solution = solve(prob, alg=opt.alg, abstol=opt.abstol,reltol=opt.reltol, dt=col.L/1000000)
    end
    return solution
end

"""
    odesystem_r!(dt, t, p, z)

The ODE system for migration ``t(z)`` and peak variance development
``τ²(z)``.

``\\frac{dt}{dz} = r(z, t(z))``

``\\frac{dτ^2}{dz} = H(z, t(z)) r(z, t(z)) + 2 τ^2(z, t(z))
\\frac{∂r}{∂t}(z,t(z))``

See also: [`solving_odesystem_r`](@ref), [`peakode`](@ref)
"""
function odesystem_r!(dt, t, p, z)
    # t[1] ... t time
    # t[2] ... τ² band variance
	col = p[1]
	prog = p[2]
	sub = p[3]
	opt = p[4]
    dt[1] = residency(z, t[1], prog.T_itp, prog.Fpin_itp, prog.pout_itp, col.L, col.d, col.df, col.gas, sub.Tchar, sub.θchar, sub.ΔCp, sub.φ₀; ng=opt.ng, vis=opt.vis, control=opt.control)
    dt[2] = peakode(z, t[1], t[2], col, prog, sub, opt)
end

"""
    peakode(z, t, τ², col, prog, sub, opt)

The second ODE function for the ODE system describing the peak variance
development ``τ²(z)``, using (in parts) automatic differentiation.

``\\frac{dτ^2}{dz} = H(z, t(z)) r(z, t(z)) + 2 τ^2(z, t(z))
\\frac{∂r}{∂t}(z,t(z))``

**TODO**: alternative to QuadGK.jl for integration which is available for
ForwardDiff.jl

See also: [`solving_odesystem_r`](@ref), [`odesystem_r!`](@ref)
"""
function peakode(z, t, τ², col, prog, sub, opt)
	# alternative function
    if opt.ng==true
        r_ng(zt) = residency(zt[1], zt[2], prog.T_itp, prog.Fpin_itp, prog.pout_itp, col.L, col.d, col.df, col.gas, sub.Tchar, sub.θchar, sub.ΔCp, sub.φ₀; ng=true, vis=opt.vis, control=opt.control)
        H_ng(z,t) = plate_height(z, t, prog.T_itp, prog.Fpin_itp, prog.pout_itp, col.L, col.d, col.df, col.gas, sub.Tchar, sub.θchar, sub.ΔCp, sub.φ₀, sub.Dag; ng=true, vis=opt.vis, control=opt.control)
        ∂r∂t_ng(z,t) = ForwardDiff.gradient(r_ng, [z, t])[2]
        return H_ng(z,t)*r_ng([z,t])^2 + 2*τ²*∂r∂t_ng(z,t)
    else
		d(z) = col.d(z)
		T(z,t) = prog.T_itp(z, t)
		η(z,t) = GasChromatographySimulator.viscosity(z, t, prog.T_itp, col.gas; vis=opt.vis)
		c(zt) = η(zt[1], zt[2])*T(zt[1], zt[2])
		∂c∂t(z,t) = ForwardDiff.gradient(c, [z,t])[2]
        #if opt.control == "Pressure"
		#    pi2(t) = prog.Fpin_itp(t)^2
        #elseif opt.control == "Flow"
            pi2(t) = inlet_pressure(t, prog.T_itp, prog.Fpin_itp, prog.pout_itp, col.L, col.d, col.gas; vis=opt.vis, control=opt.control)^2
        #end
		po2(t) = prog.pout_itp(t)^2
		∂pi2∂t(t) = ForwardDiff.derivative(pi2, t)
		∂po2∂t(t) = ForwardDiff.derivative(po2, t)
        κ(z,t) = flow_restriction(z, t, prog.T_itp, col.d, col.gas; vis=opt.vis)
		κL(t) = κ(col.L,t)
		∂κ∂t(z,t) = quadgk(y -> d(y)^-4*∂c∂t(y,t), 0.0, z)[1]
        ∂κL∂t(t) = ∂κ∂t(col.L,t)
		e(z,t) = κ(z,t)/κL(t)
		∂e∂t(z,t) = (∂κ∂t(z,t)*κL(t)-κ(z,t)*∂κL∂t(t))/κL(t)^2
		p(z,t) = pressure(z, t, prog.T_itp, prog.Fpin_itp, prog.pout_itp, col.L, col.d, col.gas; vis=opt.vis, control=opt.control)
		∂p∂t(z,t) = 1/(2*p(z,t))*(∂pi2∂t(t)-(∂e∂t(z,t)*(pi2(t)-po2(t))+e(z,t)*(∂pi2∂t(t)-∂po2∂t(t))))
		rM(z,t) = mobile_phase_residency(z, t, prog.T_itp, prog.Fpin_itp, prog.pout_itp, col.L, col.d, col.gas; vis=opt.vis, control=opt.control)
		a(z,t) = κL(t)*p(z,t)
		∂a∂t(z,t) = ∂κL∂t(t)*p(z,t)+κL(t)*∂p∂t(z,t)
		b(zt) = T(zt[1],zt[2])*(pi2(zt[2])-po2(zt[2]))
		∂b∂t(z,t) = ForwardDiff.gradient(b, [z,t])[2]
		∂rM∂t(z,t) = 64*d(z)^2*((∂a∂t(z,t)*b([z,t])-a(z,t)*∂b∂t(z,t))/b([z,t])^2)

		k(zt) = retention_factor(zt[1], zt[2], prog.T_itp, col.d, col.df, sub.Tchar, sub.θchar, sub.ΔCp, sub.φ₀)
        ∂k∂t(z,t) = ForwardDiff.gradient(k, [z, t])[2]

		r(z,t) = residency(z, t, prog.T_itp, prog.Fpin_itp, prog.pout_itp, col.L, col.d, col.df, col.gas, sub.Tchar, sub.θchar, sub.ΔCp, sub.φ₀; vis=opt.vis, control=opt.control)
        H(z,t) = plate_height(z, t, prog.T_itp, prog.Fpin_itp, prog.pout_itp, col.L, col.d, col.df, col.gas, sub.Tchar, sub.θchar, sub.ΔCp, sub.φ₀, sub.Dag; vis=opt.vis, control=opt.control)
        ∂r∂t(z,t) = ∂rM∂t(z,t)*(1+k([z,t])) + rM(z,t)*∂k∂t(z,t)
        # AutoDiff problems with the intgral in flow_restriction κ and κL makes
        # it nessecary to partly calculate the derivative manually
        # package Quadrature.jl makes it possible to differentiate for parameters
        # inside the integral, but not for the differentiation for upper or lower bound
        # of the integral (issue #56)
        # if this can be solved, it should be straigthforward with AutoDiff
        return H(z,t)*r(z,t)^2 + 2*τ²*∂r∂t(z,t)
    end
end
#---End-Solving-Functions---

#---Begin-Result-Functions---
"""
    peaklist(sol, par)

Construct a DataFrame with the peak list of the solution `sol` of the
simulation of the GC system defined by `par`.

# Output

The peaklist DataFrame consists of the entrys:
* `Name`: Name of the solute.
* `tR`: Retention time of the solute (in s).
* `τR`: Peak width of the solute (in s).
* `TR`: Temperature of the end of the column at the retention time (in °C).
* `σR`: Band width of the solute at retention time (in m).
* `uR`: Solute velocity at retention time (in m/s).
* `kR`: Retention factor of the solute at retention time.

# Examples

```julia
julia> pl = peaklist(sol, par)
...
```
"""
function peaklist(sol, par)
	n = length(par.sub)
    # sol is solution from ODE system
    Name = Array{String}(undef, n)
    tR = Array{Float64}(undef, n)
    TR = Array{Float64}(undef, n)
    σR = Array{Float64}(undef, n)
    uR = Array{Float64}(undef, n)
    τR = Array{Float64}(undef, n)
    kR = Array{Float64}(undef, n)
    Res = fill(NaN, n)
    Threads.@threads for i=1:n
        Name[i] = par.sub[i].name
        if sol[i].t[end]==par.col.L
            tR[i] = sol[i].u[end][1]
            TR[i] = par.prog.T_itp(par.col.L, tR[i]) - 273.15
            uR[i] = 1/residency(par.col.L, tR[i], par.prog.T_itp, par.prog.Fpin_itp, par.prog.pout_itp, par.col.L, par.col.d, par.col.df, par.col.gas, par.sub[i].Tchar, par.sub[i].θchar, par.sub[i].ΔCp, par.sub[i].φ₀; vis=par.opt.vis, control=par.opt.control)
            τR[i] = sqrt(sol[i].u[end][2])
            σR[i] = τR[i]*uR[i]
            kR[i] = retention_factor(par.col.L, tR[i], par.prog.T_itp, par.col.d, par.col.df, par.sub[i].Tchar, par.sub[i].θchar, par.sub[i].ΔCp, par.sub[i].φ₀)
        else
            tR[i] = NaN
            TR[i] = NaN
            uR[i] = NaN
            τR[i] = NaN
            σR[i] = NaN
            kR[i] = NaN
        end
    end
    df = sort!(DataFrame(Name = Name, tR = tR, τR = τR, TR=TR, σR = σR, uR = uR, kR = kR, ), [:tR])
    Threads.@threads for i=1:n-1
        Res[i] = (df.tR[i+1] - df.tR[i])/(2*(df.τR[i+1] + df.τR[i]))
    end
    df[!, :Res] = Res
    return df
end

"""
    peaklist(sol, peak, par)

Construct a DataFrame with the peak list of the solution `sol` and `peak` of
the simulation of the GC system defined by `par`.

# Output

The peaklist DataFrame consists of the entrys:
* `Name`: Name of the solute.
* `tR`: Retention time of the solute (in s).
* `τR`: Peak width of the solute (in s).
* `TR`: Temperature of the end of the column at the retention time (in °C).
* `σR`: Band width of the solute at retention time (in m).
* `uR`: Solute velocity at retention time (in m/s).
* `kR`: Retention factor of the solute at retention time.

# Examples

```julia
julia> pl = peaklist(sol, peak, par)
...
```
"""
function peaklist(sol, peak, par)
	n = length(par.sub)
    Name = Array{String}(undef, n)
    tR = Array{Float64}(undef, n)
    TR = Array{Float64}(undef, n)
    σR = Array{Float64}(undef, n)
    uR = Array{Float64}(undef, n)
    τR = Array{Float64}(undef, n)
    kR = Array{Float64}(undef, n)
    # add resolution
    Res = fill(NaN, n)
    Threads.@threads for i=1:n
        Name[i] = par.sub[i].name
        if sol[i].t[end]==par.col.L
            tR[i] = sol[i].u[end]
            TR[i] = par.prog.T_itp(par.col.L, tR[i]) - 273.15
            uR[i] = 1/residency(par.col.L, tR[i], par.prog.T_itp, par.prog.Fpin_itp, par.prog.pout_itp, par.col.L, par.col.d, par.col.df, par.col.gas, par.sub[i].Tchar, par.sub[i].θchar, par.sub[i].ΔCp, par.sub[i].φ₀; vis=par.opt.vis, control=par.opt.control)
            τR[i] = sqrt(peak[i].u[end])
            σR[i] = τR[i]*uR[i]
            kR[i] = retention_factor(par.col.L, tR[i], par.prog.T_itp, par.col.d, par.col.df, par.sub[i].Tchar, par.sub[i].θchar, par.sub[i].ΔCp, par.sub[i].φ₀)
        else
            tR[i] = NaN
            TR[i] = NaN
            uR[i] = NaN
            τR[i] = NaN
            σR[i] = NaN
            kR[i] = NaN
        end
    end
    df = sort!(DataFrame(Name = Name, tR = tR, τR = τR, TR=TR, σR = σR, uR = uR, kR = kR, ), [:tR])
    Threads.@threads for i=1:n-1
        Res[i] = (df.tR[i+1] - df.tR[i])/(2*(df.τR[i+1] + df.τR[i]))
    end
    df[!, :Res] = Res
    return df
end

"""
    sol_extraction(sol, par)

Extract the points z=t, t=u1, τ²=u2 from the solution `sol` of
the ODE system of the GC system defined by `par` and exports them in a DataFrame.

# Examples

```julia
df_sol = sol_extraction(sol, par)
...
```
"""
function sol_extraction(sol, par)
    # extract the points z=t, t=u1, τ²=u2 from the solution of
    # the ODE system
	n = length(par.sub)
    sol_z = Array{Any}(undef, n)
    sol_t = Array{Any}(undef, n)
    sol_τ² = Array{Any}(undef, n)
    solutes = Array{String}(undef, n)
    for i=1:n
        sol_z[i] = sol[i].t
        temp_t = Array{Float64}(undef, length(sol[i].t))
        temp_τ² = Array{Float64}(undef, length(sol[i].t))
        for j=1:length(sol[i].t)
                temp_t[j] = sol[i].u[j][1]
                temp_τ²[j] = sol[i].u[j][2]
        end
        sol_t[i] = temp_t
        sol_τ²[i] = temp_τ²
        solutes[i] = par.sub[i].name
    end
    df_sol = DataFrame(name=solutes, z=sol_z, t=sol_t, τ²=sol_τ²)
    return df_sol
end

"""
    sol_extraction(sol, peak, par)

Extract the points z_t=sol.t, t=sol.u, z_τ²=peak.t and τ²=peak.u from the
solution `sol` and `peak` of the ODEs of the GC system defined by `par` and exports them in a DataFrame.

# Examples

```julia
df_sol = sol_extraction(sol, peak, par)
...
```
"""
function sol_extraction(sol, peak, par)
    # extract the points z=t, t=u, from the solution of
    # the first ODE (sol_tz)
    # and the points z=t, τ²=u, from the solution of
    # the second ODE (peak_τz)
	n = length(par.sub)
    sol_z = Array{Any}(undef, n)
    sol_t = Array{Any}(undef, n)
    peak_z = Array{Any}(undef, n)
    peak_τ² = Array{Any}(undef, n)
    solutes = Array{String}(undef, n)
    for i=1:n
        sol_z[i] = sol[i].t
        sol_t[i] = sol[i].u
        peak_z[i] = peak[i].t
        peak_τ²[i] = peak[i].u
		solutes[i] = par.sub[i].name
	end
    df_sol = DataFrame(name=solutes, z_t=sol_z, t=sol_t, z_τ²=peak_z, τ²=peak_τ²)
    return df_sol
end

#----begin-notebooks-functions----------------------------------------------------------------------
## Physical-model-functions
include("./Model.jl")
## UI-functions
include("./UI.jl")
## plot-functions
include("./Plot.jl")
## misc-functions
include("./Misc.jl")
##---end-notebooks-functions---------------------------------------------------------------------------

end # module


text:
using GPUAcceleratedTracking, DrWatson, Tracking, GNSSSignals, StructArrays, ProgressMeter
@quickactivate "GPUAcceleratedTracking"

allparams = Dict(
    "processor"   => ["CPU"],
    "GNSS"  => ["GPSL5"],
    "num_samples" => 2 .^ (15:18),
    "num_ants" => [1,4],
    "num_correlators" => [3],
    "algorithm" => [
        "1_4_cplx_multi_textmem",
        # "2_4_cplx_multi_textmem",
        # "3_4_cplx_multi_textmem",
        # "4_4_cplx_multi_textmem",
        # "5_4_cplx_multi_textmem"
    ]
)

dicts = dict_list(allparams)

@showprogress 1 "Benchmarking kernel algorithms" for (_, d) in enumerate(dicts)
    benchmark_results = run_kernel_benchmark(d)
    @tagsave(
        datadir("benchmarks/kernel/atbenchmark", savename("KernelBenchmark", d, "jld2")),
        benchmark_results
    )
end

max_stars_repo_name: mfkiwl/GPUAcceleratedTracking

max_issues_repo_name: mfkiwl/GPUAcceleratedTracking

max_forks_repo_name: mfkiwl/GPUAcceleratedTracking


text:
mutable struct Agent
  name::AbstractString
  capacity::Int64
  capacity_reduction::Int64
end
hash(a::Agent) = hash(a.name)
isequal(a::Agent, b::Agent) = a.name == b.name
Agent(name::AbstractString) = Agent(name, 1, 1)

@enum PreferenceRanking strictly_prefers = 2 weakly_prefers = 1 indifferent = 0 weakly_not_prefer = -1 strictly_not_prefer = -2

struct Preferences
  a::Agent
  b::Agent
  ranking::PreferenceRanking
end

struct RankedPreferences
  ranking::Dict{Agent, Vector{Agent}}
end
RankedPreferences() = RankedPreferences(Dict{Agent, Vector{Agent}}())

function add_rank!(rankings::RankedPreferences, a::Agent, b::Vector{Agent})
  rankings.ranking[a] = b
end
copy(r::RankedPreferences) = RankedPreferences(copy(r.ranking))

struct Match
  results::Dict{Agent, Agent}
end

function show(io::IO, m::Match)
  s = ""
  for (key, value) in m.results
    s = s * key.name * " matches with " * value.name * "\n"
  end
  print(io, s)
end

>(a::Agent, b::Agent) = Preference(a, b, PreferenceRanking(2))
>=(a::Agent, b::Agent) = Preference(a, b, PreferenceRanking(1))
~(a::Agent, b::Agent) = Preference(a, b, PreferenceRanking(0))
<=(a::Agent, b::Agent) = Preference(a, b, PreferenceRanking(-1))
<(a::Agent, b::Agent) = Preference(a, b, PreferenceRanking(-2))

function prefers(ranking::RankedPreferences, agent::Agent, a::Agent, b::Agent)
  ranked_list = ranking.ranking[agent]
  for i in 1:length(ranked_list)
    if ranked_list[i] == a
      return true
    elseif ranked_list[i] == b
      return false
    end
  end
end


text:
# Tests based on the Raster API tutorial found at http://gdal.org/gdal_tutorial.html

using GDAL

utmsmall = joinpath(@__DIR__, "data/utmsmall.tif")
utmsmall_copy = joinpath(@__DIR__, "tmp/utmsmall.tif")

# Opening the File
dataset = GDAL.gdalopen(utmsmall, GDAL.GA_ReadOnly) # file from GDAL SVN

# Getting Dataset Information
driver = GDAL.gdalgetdriverbyname("GTiff")
@test GDAL.gdalgetdrivershortname(driver) == "GTiff"
@test GDAL.gdalgetdriverlongname(driver) == "GeoTIFF"

@test GDAL.gdalgetrasterxsize(dataset) == 100
@test GDAL.gdalgetrasterysize(dataset) == 100
@test GDAL.gdalgetrastercount(dataset) == 1
nad27_prefix = "PROJCS[\"NAD27 / UTM zone 11N\",GEOGCS[\"NAD27\",DATUM[\"North_American_Datum_1927\","
@test startswith(GDAL.gdalgetprojectionref(dataset), nad27_prefix)

geotransform = fill(0.0, 6)
GDAL.gdalgetgeotransform(dataset, geotransform)
@test geotransform == [440720.0, 60.0, 0.0, 3.75132e6, 0.0, -60.0]

# Fetching a Raster Band
band = GDAL.gdalgetrasterband(dataset, 1)

blockxsize, blockysize = Ref(Cint(-1)), Ref(Cint(-1))
GDAL.gdalgetblocksize(band, blockxsize, blockysize)
@test blockxsize[] == 100
@test blockysize[] == 81

@test GDAL.gdalgetdatatypename(GDAL.gdalgetrasterdatatype(band)) == "Byte"
@test GDAL.gdalgetcolorinterpretationname(GDAL.gdalgetrastercolorinterpretation(band)) ==
      "Gray"

gotmin, gotmax = Ref(Cint(-1)), Ref(Cint(-1))
@test GDAL.gdalgetrasterminimum(band, gotmin) == 0.0
@test GDAL.gdalgetrastermaximum(band, gotmax) == 255.0
@test gotmin[] == gotmax[] == 0

@test GDAL.gdalgetoverviewcount(band) == 0
# note that if there is no raster color table, no error (C_None) is set by GDAL
# so there is no GDALError thrown here
@test GDAL.gdalgetrastercolortable(band) == C_NULL

# Reading Raster Data
xsize = GDAL.gdalgetrasterbandxsize(band)
@test xsize == 100
scanline = fill(0.0f0, xsize)
GDAL.gdalrasterio(
    band,
    GDAL.GF_Read,
    0,
    0,
    xsize,
    1,
    scanline,
    xsize,
    1,
    GDAL.GDT_Float32,
    0,
    0,
)
@test scanline == Float32[
    107.0f0,
    123.0f0,
    132.0f0,
    115.0f0,
    132.0f0,
    132.0f0,
    140.0f0,
    132.0f0,
    132.0f0,
    132.0f0,
    107.0f0,
    132.0f0,
    107.0f0,
    132.0f0,
    132.0f0,
    107.0f0,
    123.0f0,
    115.0f0,
    156.0f0,
    148.0f0,
    107.0f0,
    132.0f0,
    107.0f0,
    115.0f0,
    99.0f0,
    123.0f0,
    99.0f0,
    74.0f0,
    115.0f0,
    82.0f0,
    115.0f0,
    115.0f0,
    107.0f0,
    123.0f0,
    123.0f0,
    99.0f0,
    123.0f0,
    123.0f0,
    115.0f0,
    115.0f0,
    107.0f0,
    90.0f0,
    99.0f0,
    107.0f0,
    107.0f0,
    99.0f0,
    123.0f0,
    107.0f0,
    140.0f0,
    123.0f0,
    123.0f0,
    115.0f0,
    99.0f0,
    132.0f0,
    123.0f0,
    115.0f0,
    115.0f0,
    123.0f0,
    132.0f0,
    115.0f0,
    123.0f0,
    132.0f0,
    214.0f0,
    156.0f0,
    165.0f0,
    148.0f0,
    115.0f0,
    148.0f0,
    156.0f0,
    148.0f0,
    140.0f0,
    165.0f0,
    156.0f0,
    197.0f0,
    156.0f0,
    197.0f0,
    140.0f0,
    173.0f0,
    156.0f0,
    165.0f0,
    148.0f0,
    156.0f0,
    206.0f0,
    214.0f0,
    181.0f0,
    206.0f0,
    173.0f0,
    222.0f0,
    206.0f0,
    255.0f0,
    214.0f0,
    173.0f0,
    214.0f0,
    255.0f0,
    214.0f0,
    247.0f0,
    255.0f0,
    230.0f0,
    206.0f0,
    197.0f0,
]

# Closing the Dataset
GDAL.gdalclose(dataset)

# Techniques for Creating Files
@test GDAL.gdalgetmetadataitem(driver, "DCAP_CREATE", "") == "YES"
@test GDAL.gdalgetmetadataitem(driver, "DCAP_CREATECOPY", "") == "YES"

# Using CreateCopy()
ds_src = GDAL.gdalopen(utmsmall, GDAL.GA_ReadOnly)

progressfunc = convert(Ptr{GDAL.GDALProgressFunc}, C_NULL)

ds_dst = GDAL.gdalcreatecopy(driver, utmsmall_copy, ds_src, 0, C_NULL, progressfunc, C_NULL)
GDAL.gdalclose(ds_dst)
GDAL.gdalclose(ds_src)

rm(utmsmall_copy)
# rest of the tutorial is skipped, perhaps we should wrap "cpl_string.h"
# to be able to build papszOptions, or write it in julia
# #include "cpl_string.h"
# ...
#     char **papszOptions = NULL;
#
#     papszOptions = CSLSetNameValue( papszOptions, "TILED", "YES" );
#     papszOptions = CSLSetNameValue( papszOptions, "COMPRESS", "PACKBITS" );
#     hDstDS = GDALCreateCopy( hDriver, pszDstFilename, hSrcDS, FALSE,
#                              papszOptions, GDALTermProgres, NULL );
#     /* Once we're done, close properly the dataset */
#     if( hDstDS != NULL )
#         GDALClose( hDstDS );
#     CSLDestroy( papszOptions );


text:
using Documenter, Weave, StatsBase
using ContinuousWavelet

# Convert jmd -> md
foreach(
    filename -> weave(
        joinpath(dirname(pathof(ContinuousWavelet)), "..", "docs", "jmd", filename),
        doctype = "github",
        fig_path = joinpath("docs", "fig"),
        fig_ext = ".svg",
        out_path = joinpath("docs", "src"),
    ),
    [
        "wavelets.jmd",
        "transform.jmd",
        "synthesis.jmd",
        "coherence.jmd",
        "utils.jmd",
        "examples.jmd"
    ]
)

# assemble docs from md sources
makedocs(
    root = joinpath(dirname(pathof(ContinuousWavelet)), "..", "docs"),
    sitename = "ContinuousWavelet",
    pages = [
        "index.md",
        "transform.md",
        "synthesis.md",
        "wavelets.md",
        "coherence.md",
        "utils.md",
        "examples.md"
    ],
)

# upload
deploydocs(repo = "github.com/hmatuschek/ContinuousWavelet.jl.git", push_preview = true)


text:
@testset "EUtils" begin

    @testset "einfo" begin
        sleep(0.3)
        res = einfo(db="pubmed")
        @test res.status == 200
        @test startswith(Dict(res.headers)["Content-Type"], "text/xml")
        body = parse_xml(String(res.body))
        @test isa(body, XMLDict.XMLDictElement)
        @test first(body)[1] != "ERROR"
    end

    @testset "esearch" begin
        sleep(0.3)
        res = esearch(db="pubmed", term="""(Asthma[MeSH Major Topic]) AND
                                        ("1/1/2018"[Date - Publication] :
                                        "3000"[Date - Publication])""")
        @test res.status == 200
        @test startswith(Dict(res.headers)["Content-Type"], "text/xml")
        body = parse_xml(String(res.body))
        @test isa(body, XMLDict.XMLDictElement)
        @test first(body)[1] != "ERROR"
    end

    @testset "epost" begin
        sleep(0.3)
        ctx = Dict()
        res = epost(ctx, db="protein", id="NP_005537")
        @test res.status == 200
        @test startswith(Dict(res.headers)["Content-Type"], "text/xml")
        body = parse_xml(String(res.body))
        @test isa(body, XMLDict.XMLDictElement)
        @test haskey(ctx, :WebEnv)
        @test haskey(ctx, :query_key)
    end

    @testset "esummary" begin
        # esummary doesn't seem to support accession numbers
        sleep(0.3)
        res = esummary(db="protein", id="15718680,157427902,119703751")
        @test res.status == 200
        @test startswith(Dict(res.headers)["Content-Type"], "text/xml")
        body = parse_xml(String(res.body))
        @test isa(body, XMLDict.XMLDictElement)
        @test first(body)[1] != "ERROR"

        sleep(0.3)
        res = esummary(db="protein", id=["15718680", "157427902", "119703751"])
        @test res.status == 200
        @test startswith(Dict(res.headers)["Content-Type"], "text/xml")
        body = parse_xml(String(res.body))
        @test isa(body, XMLDict.XMLDictElement)
        @test first(body)[1] != "ERROR"

        # esearch then esummary
        sleep(0.3)
        query = "asthma[mesh] AND leukotrienes[mesh] AND 2009[pdat]"
        ctx = Dict()
        res = esearch(ctx, db="pubmed", term=query, usehistory=true, retmode="xml")
        @test res.status == 200
        sleep(0.3)
        res = esummary(ctx, db="pubmed")
        @test res.status == 200

        sleep(0.3)
        ctx = Dict()
        res = esearch(ctx, db="pubmed", term=query, usehistory=true, retmode="json")
        @test res.status == 200
        sleep(0.3)
        res = esummary(ctx, db="pubmed")
        @test res.status == 200
    end

    @testset "efetch" begin
        sleep(0.3)
        res = efetch(db="nuccore", id="NM_001178.5", retmode="xml", idtype="acc")
        @test res.status == 200
        @test startswith(Dict(res.headers)["Content-Type"], "text/xml")
        @test isa(parse_xml(String(res.body)), XMLDict.XMLDictElement)

        # epost then efetch
        sleep(0.3)
        ctx = Dict()
        res = epost(ctx, db="protein", id="NP_005537")
        @test res.status == 200
        res = efetch(ctx, db="protein", retmode="xml")
        @test res.status == 200
    end

    @testset "elink" begin
        sleep(0.3)
        res = elink(dbfrom="protein", db="gene", id="NM_001178.5")
        @test res.status == 200
        @test startswith(Dict(res.headers)["Content-Type"], "text/xml")
        @test isa(parse_xml(String(res.body)), XMLDict.XMLDictElement)
    end

    @testset "egquery" begin
        sleep(0.3)
        res = egquery(term="""(Asthma[MeSH Major Topic]) AND
                              ("1/1/2018"[Date - Publication] :
                              "3000"[Date - Publication])""")
        @test res.status == 200
        @test startswith(Dict(res.headers)["Content-Type"], "text/xml")
        @test isa(parse_xml(String(res.body)), XMLDict.XMLDictElement)
    end

    @testset "espell" begin
        sleep(0.3)
        res = espell(db="pmc", term="fiberblast cell grwth")
        @test res.status == 200
        @test startswith(Dict(res.headers)["Content-Type"], "text/xml")
        body = parse_xml(String(res.body))
        @test isa(body, XMLDict.XMLDictElement)
        replaced_spell = body["SpelledQuery"]["Replaced"][2]
        @test replaced_spell == "growth"
    end

    @testset "ecitmatch" begin
        sleep(0.3)
        res = ecitmatch(
            db="pubmed",
            retmode="xml",
            bdata="proc natl acad sci u s a|1991|88|3248|mann bj|Art1|")
        @test res.status == 200
        # ECitMatch is mysterious: it returns a plain text data even though it
        # requires retmode="xml".
        # @show res.headers["Content-Type"]
    end
end


text:
export build_lsoptim_objective

function build_lsoptim_objective(prob::DiffEqBase.DEProblem,t,data,alg;
                                 prob_generator = STANDARD_PROB_GENERATOR,
                                 kwargs...)
  vec_data = vec(data)
  data_length = length(vec_data)
  cost_function = function (out,p)
  tmp_prob = prob_generator(prob,p)
  sol = solve(tmp_prob,alg;saveat=t,save_everystep=false,dense=false,kwargs...)
  y = vec(sol)
  y_excess = length(y) - data_length + 1
  y = y[y_excess:length(y)]
  out .= y .- vec_data
  end
end


text:
# Utility method for state gradient accumulation.
@inline fold_sum(::Nothing, ::Nothing) = nothing
@inline fold_sum(a::A, ::Nothing) where A = a
@inline fold_sum(::Nothing, a::A) where A = a
@inline fold_sum(a::A, b::A) where A = a + b

@inline _sum(::Nothing) = nothing
@inline _sum(x::Vector) = sum(x)

function choice_gradients(trace::VectorTrace{UnfoldType,T,U}, selection::Selection, retval_grad) where {T,U}
    kernel_has_grads = has_argument_grads(trace.gen_fn.kernel)
    if kernel_has_grads[1]
        error("Cannot differentiate with respect to index in unfold")
    end

    args = get_args(trace)
    (len, init_state, params) = unpack_args(args)
    state_has_grad = kernel_has_grads[2]
    params_has_grad = kernel_has_grads[3 : end]
    value_choices = choicemap()
    gradient_choices = choicemap()
    params_grad = Vector(undef, length(params_has_grad))
    for (i, has_grad) in enumerate(params_has_grad)
        if has_grad
            params_grad[i] = Vector(undef, len)
        else
            params_grad[i] = nothing
        end
    end

    local kernel_arg_grads :: Tuple = (nothing, nothing)
    for key = len : -1 : 1 # walks backward over chain
        subtrace = trace.subtraces[key]
        subselection = getindex(selection, key)
        kernel_retval_grad = (retval_grad == nothing) ? nothing : retval_grad[key]
        if state_has_grad
            kernel_retval_grad = fold_sum(kernel_retval_grad, kernel_arg_grads[2])
        end
        kernel_arg_grads, kernel_value_choices, kernel_gradient_choices = choice_gradients(subtrace, subselection, kernel_retval_grad)
        set_submap!(value_choices, key, kernel_value_choices)
        set_submap!(gradient_choices, key, kernel_gradient_choices)
        @assert kernel_arg_grads[1] == nothing
        state_has_grad || @assert kernel_arg_grads[2] == nothing
        for (i, (grad, has_grad)) in enumerate(zip(kernel_arg_grads[3:end], params_has_grad))
            if has_grad
                params_grad[i][key] = grad
            end
        end
    end
    ((nothing, kernel_arg_grads[2], map(_sum, params_grad)...), value_choices, gradient_choices)
end

function accumulate_param_gradients!(trace::VectorTrace{UnfoldType,T,U}, retval_grad) where {T,U}
    kernel_has_grads = has_argument_grads(trace.gen_fn.kernel)
    if kernel_has_grads[1]
        error("Cannot differentiate with respect to index in unfold")
    end

    args = get_args(trace)
    (len, init_state, params) = unpack_args(args)
    state_has_grad = kernel_has_grads[2]
    params_has_grad = kernel_has_grads[3:end]
    params_grad = Vector(undef, length(params_has_grad))
    for (i, has_grad) in enumerate(params_has_grad)
        if has_grad
            params_grad[i] = Vector(undef, len)
        else
            params_grad[i] = nothing
        end
    end

    local kernel_arg_grads::Tuple = (nothing, nothing) # initialize for (len index, state)
    for key = len : -1 : 1 # walks backward over chain
        subtrace = trace.subtraces[key]
        kernel_retval_grad = (retval_grad == nothing) ? nothing : retval_grad[key]
        if state_has_grad
            kernel_retval_grad = fold_sum(kernel_retval_grad, kernel_arg_grads[2])
        end
        kernel_arg_grads = accumulate_param_gradients!(subtrace, kernel_retval_grad)
        @assert kernel_arg_grads[1] == nothing
        state_has_grad || @assert kernel_arg_grads[2] == nothing
        for (i, (grad, has_grad)) in enumerate(zip(kernel_arg_grads[3:end], params_has_grad))
            if has_grad
                params_grad[i][key] = grad
            end
        end
    end
    (nothing, kernel_arg_grads[2], map(_sum, params_grad)...)
end


text:
# Code that doesn't parse under Julia 0.5
# This is part of fundef.jl

function gather_wheres(ex)
    if @capture(ex, (f_ where {params1__}))
        f2, params2 = gather_wheres(f)
        (f2, (params1..., params2...))
    else
        (ex, ())
    end
end

function longdef1_where(ex)
    if @capture(ex, (fcall_ = body_))
        fcall2, whereparams = gather_wheres(fcall)
        @match fcall2 begin
            (f_(args__)) =>
                @q function $f($(args...)) where {$(whereparams...)}
                    $body end
            (f_(args__)::rtype_) =>
                @q function ($f($(args...))::$rtype) where {$(whereparams...)}
                    $body end
            _ => ex
        end
    else
        ex
    end
end
function splitwhere(fdef)
    @assert(@capture(longdef1(fdef),
                     function (fcall_ | fcall_) body_ end),
            "Not a function definition: $fdef")
    fcall2, whereparams = gather_wheres(fcall)
    return fcall2, body, whereparams
end

"""
     combinedef(dict::Dict)

`combinedef` is the inverse of `splitdef`. It takes a splitdef-like dict
and returns a function definition. """
function combinedef(dict::Dict)
    rtype = get(dict, :rtype, :Any)
    params = get(dict, :params, [])
    # LightGraph.jl has outer-constructors with both normal-params and where-params
    # TODO: write test
    wparams = get(dict, :whereparams, [])
    name = dict[:name]
    name_p = isempty(params) ? name : :($name{$(params...)})
    if isempty(wparams)
        :(function $name_p($(dict[:args]...);
                           $(dict[:kwargs]...))::$rtype
          $(dict[:body])
          end)
    else
        :(function $name_p($(dict[:args]...);
                           $(dict[:kwargs]...))::$rtype where {$(wparams...)}
          $(dict[:body])
          end)
    end
end


text:
p = collect(permutations(1:6))
io = open("output.txt","w")
write(io,string(length(p))*"\n")
for i in p
	write(io,join(string.(i)," ")*"\n")
end
close(io)

max_stars_repo_name: Ved-Mahajan/JustForFun

max_issues_repo_name: Ved-Mahajan/JustForFun

max_forks_repo_name: Ved-Mahajan/JustForFun