text:
### BasicContMuvParameter

mutable struct BasicContMuvParameter <: Parameter{Continuous, Multivariate}
  key::Symbol
  index::Integer
  pdf::Union{ContinuousMultivariateDistribution, Void}
  prior::Union{ContinuousMultivariateDistribution, Void}
  setpdf::Union{Function, Void}
  setprior::Union{Function, Void}
  loglikelihood!::Union{Function, Void}
  logprior!::Union{Function, Void}
  logtarget!::Union{Function, Void}
  gradloglikelihood!::Union{Function, Void}
  gradlogprior!::Union{Function, Void}
  gradlogtarget!::Union{Function, Void}
  tensorloglikelihood!::Union{Function, Void}
  tensorlogprior!::Union{Function, Void}
  tensorlogtarget!::Union{Function, Void}
  dtensorloglikelihood!::Union{Function, Void}
  dtensorlogprior!::Union{Function, Void}
  dtensorlogtarget!::Union{Function, Void}
  uptogradlogtarget!::Union{Function, Void}
  uptotensorlogtarget!::Union{Function, Void}
  uptodtensorlogtarget!::Union{Function, Void}
  diffmethods::Union{DiffMethods, Void}
  diffopts::Union{DiffOptions, Void}
  states::VariableStateVector
  state::ParameterState{Continuous, Multivariate}

  function BasicContMuvParameter(
    key::Symbol,
    index::Integer,
    pdf::Union{ContinuousMultivariateDistribution, Void},
    prior::Union{ContinuousMultivariateDistribution, Void},
    setpdf::Union{Function, Void},
    setprior::Union{Function, Void},
    ll::Union{Function, Void},
    lp::Union{Function, Void},
    lt::Union{Function, Void},
    gll::Union{Function, Void},
    glp::Union{Function, Void},
    glt::Union{Function, Void},
    tll::Union{Function, Void},
    tlp::Union{Function, Void},
    tlt::Union{Function, Void},
    dtll::Union{Function, Void},
    dtlp::Union{Function, Void},
    dtlt::Union{Function, Void},
    uptoglt::Union{Function, Void},
    uptotlt::Union{Function, Void},
    uptodtlt::Union{Function, Void},
    diffmethods::Union{DiffMethods, Void},
    diffopts::Union{DiffOptions, Void},
    states::VariableStateVector,
    state::ParameterState{Continuous, Multivariate}
  )
    args = (setpdf, setprior, ll, lp, lt, gll, glp, glt, tll, tlp, tlt, dtll, dtlp, dtlt, uptoglt, uptotlt, uptodtlt)
    fnames = fieldnames(BasicContMuvParameter)[5:21]

    # Check that all generic functions have correct signature
    for i in 1:17
      if isa(args[i], Function) &&
        isa(args[i], Function) &&
        !method_exists(args[i], (BasicContMuvParameterState, VariableStateVector))
        error("$(fnames[i]) has wrong signature")
      end
    end

    new(
      key,
      index,
      pdf,
      prior,
      setpdf,
      setprior,
      ll,
      lp,
      lt,
      gll,
      glp,
      glt,
      tll,
      tlp,
      tlt,
      dtll,
      dtlp,
      dtlt,
      uptoglt,
      uptotlt,
      uptodtlt,
      diffmethods,
      diffopts,
      states,
      state
    )
  end
end

function BasicContMuvParameter!(
  parameter::BasicContMuvParameter,
  setpdf::Union{Function, Void},
  setprior::Union{Function, Void},
  ll::Union{Function, Void},
  lp::Union{Function, Void},
  lt::Union{Function, Void},
  gll::Union{Function, Void},
  glp::Union{Function, Void},
  glt::Union{Function, Void},
  tll::Union{Function, Void},
  tlp::Union{Function, Void},
  tlt::Union{Function, Void},
  dtll::Union{Function, Void},
  dtlp::Union{Function, Void},
  dtlt::Union{Function, Void},
  uptoglt::Union{Function, Void},
  uptotlt::Union{Function, Void},
  uptodtlt::Union{Function, Void}
)
  args = (setpdf, setprior, ll, lp, lt, gll, glp, glt, tll, tlp, tlt, dtll, dtlp, dtlt, uptoglt, uptotlt, uptodtlt)

  # Define setpdf (i = 1) and setprior (i = 2)
  for (i, setter, distribution) in ((1, :setpdf, :pdf), (2, :setprior, :prior))
    setfield!(
      parameter,
      setter,
      if isa(args[i], Function)
        (_state::BasicContMuvParameterState, _states::VariableStateVector) ->
          setfield!(parameter, distribution, args[i](_state, _states))
      else
        nothing
      end
    )
  end

  # Define loglikelihood! (i = 3) and gradloglikelihood! (i = 6)
  # plfield stands for parameter likelihood-related field respectively
  for (i, plfield) in ((3, :loglikelihood!), (6, :gradloglikelihood!))
    setfield!(
      parameter,
      plfield,
      if isa(args[i], Function)
        _state::BasicContMuvParameterState -> args[i](_state, parameter.states)
      else
        nothing
      end
    )
  end

  # Define logprior! (i = 4) and gradlogprior! (i = 7)
  # ppfield and spfield stand for parameter prior-related field and state prior-related field repsectively
  for (i , ppfield, spfield, f) in ((4, :logprior!, :logprior, logpdf), (7, :gradlogprior!, :gradlogprior, gradlogpdf))
    setfield!(
      parameter,
      ppfield,
      if isa(args[i], Function)
        _state::BasicContMuvParameterState -> args[i](_state, parameter.states)
      else
        if (
            isa(parameter.prior, ContinuousMultivariateDistribution) &&
            method_exists(f, (typeof(parameter.prior), Vector{eltype(parameter.prior)}))
          ) ||
          isa(args[2], Function)
          _state::BasicContMuvParameterState -> setfield!(_state, spfield, f(getfield(parameter, :prior), _state.value))
        else
          nothing
        end
      end
    )
  end

  # Define logtarget! (i = 5) and gradlogtarget! (i = 8)
  # ptfield, plfield and ppfield stand for parameter target, likelihood and prior-related field respectively
  # stfield, slfield and spfield stand for state target, likelihood and prior-related field respectively
  for (i , ptfield, plfield, ppfield, stfield, slfield, spfield, f) in (
    (5, :logtarget!, :loglikelihood!, :logprior!, :logtarget, :loglikelihood, :logprior, logpdf),
    (8, :gradlogtarget!, :gradloglikelihood!, :gradlogprior!, :gradlogtarget, :gradloglikelihood, :gradlogprior, gradlogpdf)
  )
    setfield!(
      parameter,
      ptfield,
      if isa(args[i], Function)
        _state::BasicContMuvParameterState -> args[i](_state, parameter.states)
      else
        if isa(args[i-2], Function) && isa(getfield(parameter, ppfield), Function)
          function (_state::BasicContMuvParameterState)
            getfield(parameter, plfield)(_state)
            getfield(parameter, ppfield)(_state)
            setfield!(_state, stfield, getfield(_state, slfield)+getfield(_state, spfield))
          end
        elseif (
            isa(parameter.pdf, ContinuousMultivariateDistribution) &&
            method_exists(f, (typeof(parameter.pdf), Vector{eltype(parameter.pdf)}))
          ) ||
          isa(args[1], Function)
          _state::BasicContMuvParameterState -> setfield!(_state, stfield, f(getfield(parameter, :pdf), _state.value))
        else
          nothing
        end
      end
    )
  end

  # Define tensorloglikelihood! (i = 9) and dtensorloglikelihood! (i = 12)
  # plfield stands for parameter likelihood-related field respectively
  for (i, plfield) in ((9, :tensorloglikelihood!), (12, :dtensorloglikelihood!))
    setfield!(
      parameter,
      plfield,
      if isa(args[i], Function)
        _state::BasicContMuvParameterState -> args[i](_state, parameter.states)
      else
        nothing
      end
    )
  end

  # Define tensorlogprior! (i = 10) and dtensorlogprior! (i = 13)
  # ppfield stands for parameter prior-related field respectively
  for (i, ppfield) in ((10, :tensorlogprior!), (13, :dtensorlogprior!))
    setfield!(
      parameter,
      ppfield,
      if isa(args[i], Function)
        _state::BasicContMuvParameterState -> args[i](_state, parameter.states)
      else
        nothing
      end
    )
  end

  # Define tensorlogtarget! (i = 11) and dtensorlogtarget! (i = 14)
  for (i , ptfield, plfield, ppfield, stfield, slfield, spfield) in (
    (
      11,
      :tensorlogtarget!, :tensorloglikelihood!, :tensorlogprior!,
      :tensorlogtarget, :tensorloglikelihood, :tensorlogprior
    ),
    (
      14,
      :dtensorlogtarget!, :dtensorloglikelihood!, :dtensorlogprior!,
      :dtensorlogtarget, :dtensorloglikelihood, :dtensorlogprior
    )
  )
    setfield!(
      parameter,
      ptfield,
      if isa(args[i], Function)
        _state::BasicContMuvParameterState -> args[i](_state, parameter.states)
      else
        if isa(args[i-2], Function) && isa(args[i-1], Function)
          function (_state::BasicContMuvParameterState)
            getfield(parameter, plfield)(_state)
            getfield(parameter, ppfield)(_state)
            setfield!(_state, stfield, getfield(_state, slfield)+getfield(_state, spfield))
          end
        else
          nothing
        end
      end
    )
  end

  # Define uptogradlogtarget!
  setfield!(
    parameter,
    :uptogradlogtarget!,
    if isa(args[15], Function)
      _state::BasicContMuvParameterState -> args[15](_state, parameter.states)
    else
      if isa(parameter.logtarget!, Function) && isa(parameter.gradlogtarget!, Function)
        function (state::BasicContMuvParameterState)
          parameter.logtarget!(state)
          parameter.gradlogtarget!(state)
        end
      else
        nothing
      end
    end
  )

  # Define uptotensorlogtarget!
  setfield!(
    parameter,
    :uptotensorlogtarget!,
    if isa(args[16], Function)
      _state::BasicContMuvParameterState -> args[16](_state, parameter.states)
    else
      if isa(parameter.logtarget!, Function) &&
        isa(parameter.gradlogtarget!, Function) &&
        isa(parameter.tensorlogtarget!, Function)
        function (state::BasicContMuvParameterState)
          parameter.logtarget!(state)
          parameter.gradlogtarget!(state)
          parameter.tensorlogtarget!(state)
        end
      else
        nothing
      end
    end
  )

  # Define uptodtensorlogtarget!
  setfield!(
    parameter,
    :uptodtensorlogtarget!,
    if isa(args[17], Function)
      _state::BasicContMuvParameterState -> args[17](_state, parameter.states)
    else
      if isa(parameter.logtarget!, Function) &&
        isa(parameter.gradlogtarget!, Function) &&
        isa(parameter.tensorlogtarget!, Function) &&
        isa(parameter.dtensorlogtarget!, Function)
        function (state::BasicContMuvParameterState)
          parameter.logtarget!(state)
          parameter.gradlogtarget!(state)
          parameter.tensorlogtarget!(state)
          parameter.dtensorlogtarget!(state)
        end
      else
        nothing
      end
    end
  )
end

BasicContMuvParameter(key::Symbol, index::Integer=0; signature::Symbol=:high, args...) =
  BasicContMuvParameter(key, Val{signature}, index; args...)

function BasicContMuvParameter(
  key::Symbol,
  ::Type{Val{:low}},
  index::Integer=0;
  pdf::Union{ContinuousMultivariateDistribution, Void}=nothing,
  prior::Union{ContinuousMultivariateDistribution, Void}=nothing,
  setpdf::Union{Function, Void}=nothing,
  setprior::Union{Function, Void}=nothing,
  loglikelihood::Union{Function, Void}=nothing,
  logprior::Union{Function, Void}=nothing,
  logtarget::Union{Function, Void}=nothing,
  gradloglikelihood::Union{Function, Void}=nothing,
  gradlogprior::Union{Function, Void}=nothing,
  gradlogtarget::Union{Function, Void}=nothing,
  tensorloglikelihood::Union{Function, Void}=nothing,
  tensorlogprior::Union{Function, Void}=nothing,
  tensorlogtarget::Union{Function, Void}=nothing,
  dtensorloglikelihood::Union{Function, Void}=nothing,
  dtensorlogprior::Union{Function, Void}=nothing,
  dtensorlogtarget::Union{Function, Void}=nothing,
  uptogradlogtarget::Union{Function, Void}=nothing,
  uptotensorlogtarget::Union{Function, Void}=nothing,
  uptodtensorlogtarget::Union{Function, Void}=nothing,
  diffmethods::Union{DiffMethods, Void}=nothing,
  diffopts::Union{DiffOptions, Void}=nothing,
  states::VariableStateVector=VariableState[],
  state::ParameterState{Continuous, Multivariate}=BasicContMuvParameterState(0)
)
  parameter = BasicContMuvParameter(key, index, pdf, prior, fill(nothing, 17)..., diffmethods, diffopts, states, state)

  BasicContMuvParameter!(
    parameter,
    setpdf,
    setprior,
    loglikelihood,
    logprior,
    logtarget,
    gradloglikelihood,
    gradlogprior,
    gradlogtarget,
    tensorloglikelihood,
    tensorlogprior,
    tensorlogtarget,
    dtensorloglikelihood,
    dtensorlogprior,
    dtensorlogtarget,
    uptogradlogtarget,
    uptotensorlogtarget,
    uptodtensorlogtarget
  )

  parameter
end

function BasicContMuvParameter(
  key::Symbol,
  ::Type{Val{:high}},
  index::Integer=0;
  pdf::Union{ContinuousMultivariateDistribution, Void}=nothing,
  prior::Union{ContinuousMultivariateDistribution, Void}=nothing,
  setpdf::Union{Function, Void}=nothing,
  setprior::Union{Function, Void}=nothing,
  loglikelihood::Union{Function, Void}=nothing,
  logprior::Union{Function, Void}=nothing,
  logtarget::Union{Function, Void}=nothing,
  gradloglikelihood::Union{Function, Void}=nothing,
  gradlogprior::Union{Function, Void}=nothing,
  gradlogtarget::Union{Function, Void}=nothing,
  tensorloglikelihood::Union{Function, Void}=nothing,
  tensorlogprior::Union{Function, Void}=nothing,
  tensorlogtarget::Union{Function, Void}=nothing,
  dtensorloglikelihood::Union{Function, Void}=nothing,
  dtensorlogprior::Union{Function, Void}=nothing,
  dtensorlogtarget::Union{Function, Void}=nothing,
  uptogradlogtarget::Union{Function, Void}=nothing,
  uptotensorlogtarget::Union{Function, Void}=nothing,
  uptodtensorlogtarget::Union{Function, Void}=nothing,
  nkeys::Integer=0,
  vfarg::Bool=false,
  diffopts::Union{DiffOptions, Void}=nothing,
  states::VariableStateVector=VariableState[],
  state::ParameterState{Continuous, Multivariate}=BasicContMuvParameterState(0)
)
  inargs = (
    setpdf,
    setprior,
    loglikelihood,
    logprior,
    logtarget,
    gradloglikelihood,
    gradlogprior,
    gradlogtarget,
    tensorloglikelihood,
    tensorlogprior,
    tensorlogtarget,
    dtensorloglikelihood,
    dtensorlogprior,
    dtensorlogtarget,
    uptogradlogtarget,
    uptotensorlogtarget,
    uptodtensorlogtarget
  )

  fnames = fieldnames(BasicContMuvParameterState)[2:13]

  if vfarg
    if nkeys > 0
      if diffopts != nothing
        error("In the case of autodiff, if nkeys is not 0, then vfarg must be false")
      end
    else
      error("If vfarg=$vfarg, nkeys must be positive, got $nkeys")
    end
  else
    @assert nkeys >= 0 "If vfarg=$vfarg, nkeys must be non-negative, got $nkeys"
  end

  if diffopts != nothing
    for i in 3:5
      if isa(inargs[i], Function) && (!isa(inargs[i+3], Function) || (diffopts.order == 2 && !isa(inargs[i+6], Function)))
        diffopts.targets[i-2] = true
      end
    end
  end

  parameter = BasicContMuvParameter(key, index, pdf, prior, fill(nothing, 17)..., DiffMethods(), diffopts, states, state)

  outargs = Union{Function, Void}[nothing for i in 1:17]

  if vfarg
    for i in 1:2
      if isa(inargs[i], Function)
        outargs[i] = (_state::BasicContMuvParameterState, _states::VariableStateVector) ->
          inargs[i](Any[s.value for s in _states])
      end
    end

    for i in 3:14
      if isa(inargs[i], Function)
        outargs[i] = (_state::BasicContMuvParameterState, _states::VariableStateVector) ->
          setfield!(_state, fnames[i-2], inargs[i](Any[s.value for s in _states]))
      end
    end

    if isa(inargs[15], Function)
      outargs[15] = (_state::BasicContMuvParameterState, _states::VariableStateVector) ->
        (_state.logtarget, _state.gradlogtarget) = inargs[15](Any[s.value for s in _states])
    end

    if isa(inargs[16], Function)
      outargs[16] = (_state::BasicContMuvParameterState, _states::VariableStateVector) ->
        (_state.logtarget, _state.gradlogtarget, _state.tensorlogtarget) = inargs[16](Any[s.value for s in _states])
    end

    if isa(inargs[17], Function)
      outargs[17] = (_state::BasicContMuvParameterState, _states::VariableStateVector) ->
        (_state.logtarget, _state.gradlogtarget, _state.tensorlogtarget, _state.dtensorlogtarget) =
          inargs[17](Any[s.value for s in _states])
    end
  else
    if nkeys > 0
      for i in 1:2
        if isa(inargs[i], Function)
          outargs[i] = (_state::BasicContMuvParameterState, _states::VariableStateVector) ->
            inargs[i](_state.value, Any[s.value for s in _states])
        end
      end

      for i in 3:14
        if isa(inargs[i], Function)
          outargs[i] = (_state::BasicContMuvParameterState, _states::VariableStateVector) ->
            setfield!(_state, fnames[i-2], inargs[i](_state.value, Any[s.value for s in _states]))
        end
      end

      if isa(inargs[15], Function)
        outargs[15] = (_state::BasicContMuvParameterState, _states::VariableStateVector) ->
          (_state.logtarget, _state.gradlogtarget) = inargs[15](_state.value, Any[s.value for s in _states])
      end

      if isa(inargs[16], Function)
        outargs[16] = (_state::BasicContMuvParameterState, _states::VariableStateVector) ->
          (_state.logtarget, _state.gradlogtarget, _state.tensorlogtarget) =
            inargs[16](_state.value, Any[s.value for s in _states])
      end

      if isa(inargs[17], Function)
        outargs[17] = (_state::BasicContMuvParameterState, _states::VariableStateVector) ->
          (_state.logtarget, _state.gradlogtarget, _state.tensorlogtarget, _state.dtensorlogtarget) =
            inargs[17](_state.value, Any[s.value for s in _states])
      end
    else
      for i in 1:2
        if isa(inargs[i], Function)
          outargs[i] = (_state::BasicContMuvParameterState, _states::VariableStateVector) -> inargs[i](_state.value)
        end
      end

      for i in 3:14
        if isa(inargs[i], Function)
          outargs[i] = (_state::BasicContMuvParameterState, _states::VariableStateVector) ->
            setfield!(_state, fnames[i-2], inargs[i](_state.value))
        end
      end

      if isa(inargs[15], Function)
        outargs[15] = (_state::BasicContMuvParameterState, _states::VariableStateVector) ->
          (_state.logtarget, _state.gradlogtarget) = inargs[15](_state.value)
      end

      if isa(inargs[16], Function)
        outargs[16] = (_state::BasicContMuvParameterState, _states::VariableStateVector) ->
          (_state.logtarget, _state.gradlogtarget, _state.tensorlogtarget) = inargs[16](_state.value)
      end

      if isa(inargs[17], Function)
        outargs[17] = (_state::BasicContMuvParameterState, _states::VariableStateVector) ->
          (_state.logtarget, _state.gradlogtarget, _state.tensorlogtarget, _state.dtensorlogtarget) =
            inargs[17](_state.value)
      end
    end
  end

  for (i, j, field, distribution, f) in ((4, 2, :logprior, :prior, logpdf), (5, 1, :logtarget, :pdf, logpdf))
    if !isa(inargs[i], Function) &&
      (
        (
          isa(getfield(parameter, distribution), ContinuousMultivariateDistribution) &&
          method_exists(f, (typeof(getfield(parameter, distribution)), Vector{eltype(getfield(parameter, distribution))}))
        ) ||
        isa(inargs[j], Function)
      )
      outargs[i] = (_state::BasicContMuvParameterState, _states::VariableStateVector) ->
        setfield!(_state, field, f(getfield(parameter, distribution), _state.value))
    end
  end

  if diffopts != nothing
    if diffopts.mode == :reverse
      for (i, field) in ((3, :closurell), (4, :closurelp), (5, :closurelt))
        if isa(inargs[i], Function)
          setfield!(
            parameter.diffmethods,
            field,
            nkeys == 0 ? inargs[i] : _x -> inargs[i](_x, Any[s.value for s in parameter.states])
          )
        end
      end

      for (i, returnname, diffresult, diffmethod) in (
        (6, :gradloglikelihood, :resultll, :tapegll),
        (7, :gradlogprior, :resultlp, :tapeglp),
        (8, :gradlogtarget, :resultlt, :tapetlt)
      )
        if !isa(inargs[i], Function) && isa(inargs[i-3], Function)
          outargs[i] = (_state::BasicContMuvParameterState, _states::VariableStateVector) ->
            setfield!(
              _state,
              returnname,
              reverse_autodiff_gradient(
                getfield(_state.diffstate, diffresult), getfield(_state.diffmethods, diffmethod), _state.value
              )
            )
        end
      end

      if !isa(inargs[15], Function) && isa(inargs[5], Function)
        outargs[15] = (_state::BasicContMuvParameterState, _states::VariableStateVector) ->
          (_state.logtarget, _state.gradlogtarget) =
            reverse_autodiff_upto_gradient(_state.diffstate.resultlt, _state.diffmethods.tapeglt, _state.value)
      end

      if diffopts.order == 2
        for (i, returnname, diffresult, diffmethod) in (
          (9, :tensorloglikelihood, :resultll, :tapetll),
          (10, :tensorlogprior, :resultlp, :tapetlp),
          (11, :tensorlogtarget, :resultlt, :tapetlt)
        )
          if !isa(inargs[i], Function) && isa(inargs[i-6], Function)
            outargs[i] = (_state::BasicContMuvParameterState, _states::VariableStateVector) ->
              setfield!(
                _state,
                returnname,
                reverse_autodiff_negative_hessian(
                  getfield(_state.diffstate, diffresult),
                  getfield(_state.diffmethods, diffmethod),
                  _state.value
                )
              )
          end
        end

        if !isa(inargs[16], Function) && isa(inargs[5], Function)
          outargs[16] = (_state::BasicContMuvParameterState, _states::VariableStateVector) ->
            (_state.logtarget, _state.gradlogtarget, _state.tensorlogtarget) =
              reverse_autodiff_upto_negative_hessian(_state.diffstate.resultlt, _state.diffmethods.tapetlt, _state.value)
        end
      end
    elseif diffopts.mode == :forward
      for (i, field) in ((3, :closurell), (4, :closurelp), (5, :closurelt))
        if isa(inargs[i], Function)
          setfield!(
            parameter.diffmethods,
            field,
            nkeys == 0 ? inargs[i] : _x::Vector -> inargs[i](_x, Any[s.value for s in parameter.states])
          )
        end
      end

      for (i, returnname, diffresult, diffmethod, diffconfig) in (
        (6, :gradloglikelihood, :resultll, :closurell, :cfggll),
        (7, :gradlogprior, :resultlp, :closurelp, :cfgglp),
        (8, :gradlogtarget, :resultlt, :closurelt, :cfgglt)
      )
        if !isa(inargs[i], Function) && isa(inargs[i-3], Function)
          outargs[i] = (_state::BasicContMuvParameterState, _states::VariableStateVector) ->
            setfield!(
              _state,
              returnname,
              forward_autodiff_gradient(
                getfield(_state.diffstate, diffresult),
                getfield(_state.diffmethods, diffmethod),
                _state.value,
                getfield(_state.diffstate, diffconfig)
              )
            )
        end
      end

      if !isa(inargs[15], Function) && isa(inargs[5], Function)
        outargs[15] = (_state::BasicContMuvParameterState, _states::VariableStateVector) ->
          (_state.logtarget, _state.gradlogtarget) =
            forward_autodiff_upto_gradient(
              _state.diffstate.resultlt, _state.diffmethods.closurelt, _state.value, _state.diffstate.cfgglt
            )
      end

      if diffopts.order == 2
        for (i, returnname, diffresult, diffmethod, diffconfig) in (
          (9, :tensorloglikelihood, :resultll, :closurell, :cfgtll),
          (10, :tensorlogprior, :resultlp, :closurelp, :cfgtlp),
          (11, :tensorlogtarget, :resultlt, :closurelt, :cfgtlt)
        )
          if !isa(inargs[i], Function) && isa(inargs[i-6], Function)
            outargs[i] = (_state::BasicContMuvParameterState, _states::VariableStateVector) ->
              setfield!(
                _state,
                returnname,
                forward_autodiff_negative_hessian(
                  getfield(_state.diffstate, diffresult),
                  getfield(_state.diffmethods, diffmethod),
                  _state.value,
                  getfield(_state.diffstate, diffconfig)
                )
              )
          end
        end

        if !isa(inargs[16], Function) && isa(inargs[5], Function)
          outargs[16] = (_state::BasicContMuvParameterState, _states::VariableStateVector) ->
            (_state.logtarget, _state.gradlogtarget, _state.tensorlogtarget) =
              forward_autodiff_upto_negative_hessian(
                _state.diffstate.resultlt, _state.diffmethods.closurelt, _state.value, _state.diffstate.cfgtlt
              )
        end
      end
    end
  end

  BasicContMuvParameter!(parameter, outargs...)

  parameter
end

function set_tapes!(parameter::BasicContMuvParameter, pstate::ParameterState{Continuous, Multivariate})
  diffgtapes = (:tapegll, :tapeglp, :tapeglt)
  diffttapes = (:tapetll, :tapetlp, :tapetlt)

  for (i, diffclosure) in ((1, :closurell), (2, :closurelp), (3, :closurelt))
    if parameter.diffopts.targets[i]
      setfield!(
        parameter.diffmethods,
        diffgtapes[i],
        ReverseDiff.GradientTape(getfield(parameter.diffmethods, diffclosure), pstate.value)
      )

      if parameter.diffopts.compiled
        setfield!(
          parameter.diffmethods,
          diffgtapes[i],
          ReverseDiff.compile(getfield(parameter.diffmethods, diffgtapes[i]))
        )
      end

      if parameter.diffopts.order == 2
        difftapes = (:tapegll, :tapeglp, :tapeglt)

        setfield!(
          parameter.diffmethods,
          diffttapes[i],
          ReverseDiff.HessianTape(getfield(parameter.diffmethods, diffclosure), pstate.value)
        )

        if parameter.diffopts.compiled
          setfield!(
            parameter.diffmethods,
            diffttapes[i],
            ReverseDiff.compile(getfield(parameter.diffmethods, diffttapes[i]))
          )
        end
      end
    end
  end
end

value_support(::Type{BasicContMuvParameter}) = Continuous
value_support(::BasicContMuvParameter) = Continuous

variate_form(::Type{BasicContMuvParameter}) = Multivariate
variate_form(::BasicContMuvParameter) = Multivariate

default_state_type(::BasicContMuvParameter) = BasicContMuvParameterState

default_state(variable::BasicContMuvParameter, value::Vector{N}, outopts::Dict) where {N<:Real} =
  BasicContMuvParameterState(
    value,
    [getfield(variable, fieldnames(BasicContMuvParameter)[i]) == nothing ? false : true for i in 10:18],
    (haskey(outopts, :diagnostics) && in(:accept, outopts[:diagnostics])) ? [:accept] : Symbol[],
    variable.diffmethods,
    variable.diffopts
  )

show(io::IO, ::Type{BasicContMuvParameter}) = print(io, "BasicContMuvParameter")


text:
# These are the objects encountered in the caching experiments
@enum Food begin
    Mealworm = 1
    Waxworm
    Peanut
    SuetPellet
    Pinenut
    Kibble
    Cricket
    Pineapple
    Salami
    Stone
end
struct MaintenanceDiet end
const N_FOODTYPE = length(instances(Food))
mutable struct FoodItem
    id::Food
    n::Int
    freshness::Float64
    eatable::Bool
    cacheable::Bool
end
function Base.:(==)(f1::FoodItem, f2::FoodItem)
    f1.id == f2.id &&
    f1.n == f2.n &&
    f1.freshness == f2.freshness &&
    f1.eatable == f2.eatable &&
    f1.cacheable == f2.cacheable
end
(id::Food)(n = 1, cacheable = true, eatable = id != Stone, freshness = 1.) = FoodItem(id, n, freshness, eatable, cacheable)
for foodtype in instances(Food)
    powderedform = Symbol(:Powdered, foodtype)
    @eval @__MODULE__() begin
        $powderedform(n = 1) = $foodtype(n, false)
    end
end
isnut(f) = f == Peanut || f == Pinenut
isnut(f::FoodItem) = isnut(f.id)
isworm(f) = f == Mealworm || f == Waxworm
isworm(f::FoodItem) = isworm(f.id)
ismeat(f) = isworm(f) || f == Salami || f == Cricket || f == SuetPellet || f == Kibble
ismeat(f::FoodItem) = ismeat(f.id)
isfruit(f) = f == Pineapple
isfruit(f::FoodItem) = isfruit(f.id)
foodindex(f) = Int(f)
foodindex(f::FoodItem) = foodindex(f.id)

mutable struct Tray
	appearance::Int
    position::Int
	closed::Bool
	eatableitems::Array{FoodItem, 1}
end
function Base.:(==)(x::Tray, y::Tray)
    x.appearance == y.appearance &&
    x.position == y.position &&
    x.closed == y.closed &&
    x.eatableitems == y.eatableitems
end

Tray(p = 0, a = rand(Int); closed = false) = Tray(a, p, closed, FoodItem[])
degrade!(tray::Tray) = degrade!(tray, 0)
function degrade!(tray::Tray, foodtype)
    for o in tray.eatableitems
        (foodtype == 0 || o.id == foodtype) && (o.freshness = 0.)
    end
end
pilfer!(tray::Tray) = empty!(tray.eatableitems)
function pilfer!(tray::Tray, foodtype)
    tray.eatableitems = filter(x -> x.id != foodtype, tray.eatableitems)
end
countcache(tray::Tray) = countcache(tray, Any)
countcache(tray::Tray, typ) = countfooditems(tray.eatableitems, typ)
function countfooditems(array, typ)
	n = 0
	for o in array
		if o.id == typ || typ == Any
			n += o.n
		end
	end
	n
end
countfooditems(tray::Tray, typ) = countfooditems(tray.eatableitems, typ)
add!(model, typ, n) = add!(model, typ(n))
add!(model, objects::AbstractVector) = add!.(Ref(model), objects)
add!(list::Vector{FoodItem}, typ, n) = add!(list, typ(n))
function add!(list::Vector{FoodItem}, new_item::FoodItem; only_push = false)
    for item in list
        if item.id == new_item.id &&
           item.freshness == new_item.freshness &&
           item.eatable == new_item.eatable &&
           item.cacheable == new_item.cacheable
            only_push || (item.n += new_item.n)
            return list
        end
    end
    push!(list, new_item)
end
add!(tray::Tray, typ, n) = add!(tray.eatableitems, typ, n)
cover!(tray::Tray) = tray.closed = true
uncover!(tray::Tray) = tray.closed = false

struct InspectionObserver
	trayappearances::Array{Int, 1}
end
InspectionObserver() = InspectionObserver(Tray[])
countinspections(o) = length(o.trayappearances)
function countinspections(o::InspectionObserver, tray::Tray, N = typemax(Int))
	n = 0
    for appearance in Iterators.take(o.trayappearances, N)
		if appearance == tray.appearance
			n += 1
		end
	end
	n
end
firstinspection(o, tray) = countinspections(o, tray, 1)
observe!(o, tray) = push!(o.trayappearances, tray.appearance)

# Models need to implement the following functions
function add! end
function remove! end
function wait! end
function countfooditems end

ext: jl


text:
# new file for plotVariable related functionality

export plotSLAM2DSolveKeys

# function plotSolveKey(dfg::AbstractDFG,
#                       refSym::Symbol,
#                       refKey::Symbol,
#                       tstSym::Symbol,
#                       tstKey::Symbol;
#                       bw::AbstractVector{<:Real}=[0.001;],
#                       plotVarHack::Function=plotPose  )
#   #
#   pts, fctT = deconvSolveKey(dfg, refSym, refKey, tstSym, tstKey)
#   Xref = manikde!(pts[2],fctT)
#   Xtst = manikde!(pts[1],fctT)
#   mmdDist = mmd(Xref, Xtst, fctT, bw=bw)

#   # FIXME not all variables are Poses, so this is really hacky
#   plotVarHack(getManifold(fctT)(), [Xref, Xtst], "$fctT\n$(refSym)_$(refKey) <-> $(tstSym)_$tstKey\nmmd=$(round(mmdDist,digits=6))", c=["red","blue"], legend=["Identity";"Delta"])
# end

"""
    $SIGNATURES

Analyzing repeat solves is easier if you can animate repeat solves.

Notes
- Repeat solves can be stored using `solveTree!(fg, storeOld=true)`.
- Experimental, see DFG #641

Example

```julia
using Images, Caesar, RoMEPlotting

fg = generateCanonicalFG_Hexagonal()

for i in 1:10
  solveTree!(fg, storeOld=true)
end

# generate all the plots in RAM
plts = plotSLAM2DSolveKeys(fg)

# export all images to a video
imgs = convert.(Matrix{RGB},plts)
Caesar.writevideo("/tmp/test.avi", imgs)
@async run(`totem /tmp/test.avi`)
```

Related

writevideo, plotSolveKey, listSolveKeys
"""
function plotSLAM2DSolveKeys( dfg::AbstractDFG,
                              pattern::Regex=r"default_\d+";
                              xmin=nothing, xmax=nothing, ymin=nothing, ymax=nothing,
                              contour=false )
  #
  kys = listSolveKeys(dfg, filterSolveKeys=pattern) |> collect |> sortDFG
  kys .|> x -> plotSLAM2D(dfg, solveKey=x, contour=contour, xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax)
end

#


text:
int main() {
   int if;
   return 1;
}

max_stars_repo_name: carlostome/Javalette-Compiler

max_issues_repo_name: carlostome/Javalette-Compiler

max_forks_repo_name: carlostome/Javalette-Compiler


text:
###########################################################

"""
$(TYPEDSIGNATURES)

Real from ASCII files all observed/measured data from a set of experiments.
It can read data for one or more proteins.
The entalphy values are *scaled* by a factor given by the argument `scalfactor`,
 which defaults to 0.004184 (Cal/mol to kJ/mol).

# Arguments

- `inpdir`: directory containing the input data
- `proteinnames`: array of strings containing the names of proteins
- `scalfactor`=0.004184: scaling factor for enthalpy, defaults
                        to 0.004184 (Cal/mol to kJ/mol)
- `discninitrows`=0: number of initial rows of the data set to discard. This is used
                    to remove some initial data often affected by strong
                    instrument noise which could bias the fitting process.

"""
function readallexperiments(inpdir::String,proteinnames::Vector{String} ;
                            scalfactor::Float64=0.004184 ,
                            discninitrows::Integer=0)
    ## swith units, from Cal/mol to kJ/mol by multiplying by 0.004184
    @assert discninitrows>=0
    starow = discninitrows+1
    #@show discninitrows

    println("Reading data from directory: $inpdir")
    println(" Scaling enthalpy values by a factor $scalfactor")

    numprot = length(proteinnames)
    data = Dict()

    for p=1:numprot

        data[proteinnames[p]] = Dict()
        data[proteinnames[p]]["sdscon"] = Float64[]
        data[proteinnames[p]]["procon"] = Float64[]
        data[proteinnames[p]]["enout"] = Float64[]
        data[proteinnames[p]]["idxdata"] = UnitRange{Int64}[]

        sdscon = data[proteinnames[p]]["sdscon"]
        procon = data[proteinnames[p]]["procon"]
        enout = data[proteinnames[p]]["enout"]
        idxdata = data[proteinnames[p]]["idxdata"]

        # get a list of files in inpdir
        fllist = readdir(joinpath(inpdir,proteinnames[p]))
        # @show fllist
        # @show occursin(proteinnames[p],fllist[1])
        # @show endswith(fllist[1],".DAT")
        # select only files starting with
        checkfile(x) = (occursin(proteinnames[p],x) && ( endswith(x,".DAT") || endswith(x,".dat") ) )
        fllist = filter(checkfile,fllist)
        println(" Input file list: ")

        startind = 1
        for fl in fllist
            println("  $fl")

            seflname = joinpath(inpdir,proteinnames[p]*"/"*fl)
            csds1,cpro1,ceout1 = readsingleexperiment(seflname)

            # discard some initial rows (strongly affected by instrument error)
            csds1 = csds1[starow:end]
            cpro1 = cpro1[starow:end]
            ceout1 = ceout1[starow:end]

            # add to dictionary
            lendata = length(ceout1)
            push!(sdscon,csds1...)
            push!(procon,cpro1...)
            push!(enout,ceout1...)
            # indices of individual experiments
            push!(idxdata,startind:(startind+lendata-1))
            startind += lendata
        end

        ## scale enthalpy
        enout .*= scalfactor

    end

    return data
end

###########################################################

"""
$(TYPEDSIGNATURES)

Real from ASCII file the observed/measured data from an experiment.
"""
function readsingleexperiment(singlefl::String)
    ## Only column 3 (Xt), 4 (Mt) and 6 (NDH) are relevant
    #
    # Xt is the concentration of SDS in the cell (in mM). This
    # is initially 0, and increases during the experiment.
    #
    # Mt is the protein concentration. This is initially the
    # value stated in the header of the file (in mM), and
    # decreases slightly during the experiment (as it is
    # diluted with the SDS).
    #
    # NDH is the energy output of that single injection in J/mol..
    #

    # REMARK:
    # For DH, INJV, XMt and NDH an initial row with a zero should be added.
    # Xt and Mt are initially one row longer than the other columns, and
    #   this is corrected for by the additional zeros.

    data,header = readdlm(singlefl,header=true)
    header = vec(header)

    idxXt = findall(x->x=="Xt",header)[1]
    idxMt = findall(x->x=="Mt",header)[1]
    @assert idxXt<idxMt

    ## DH	INJV	Xt	Mt	XMt	NDH
    # the last row contains only Xt and Mt at the beginning...
    data2 = data[1:end-1,:]
    lasttXtMt = data[end,:]

    idx=findall(x->x=="Xt",header)[1]
    # Xt in the last row of data is in the second column
    sdscon = convert.(Float64,[data2[:,idx]; lasttXtMt[2]])

    idx=findall(x->x=="Mt",header)[1]
    # Xt in the last row of data is in the third column
    procon = convert.(Float64,[data2[:,idx]; lasttXtMt[3]])

    idx=findall(x->x=="NDH",header)[1]
    enout = convert.(Float64,[0.0; data2[:,idx]])

    return sdscon,procon,enout
end

###########################################################


text:
export Atoms

"""
    Atoms{T<:AbstractFloat}

Basic atomic parameters: element symbols, numbers and masses

Masses are converted to atomic units.
Constructed using either element symbols or masses.

```jldoctest
julia> Atoms(:H)
Atoms{Float64}([:H], [1], [1837.4715941070515])

julia> Atoms([:H, :H, :H, :C])
Atoms{Float64}([:H, :H, :H, :C], [1, 1, 1, 6], [1837.4715941070515, 1837.4715941070515, 1837.4715941070515, 21894.713607956142])

julia> Atoms([100, 200])
Atoms{Float64}([:X, :X], [0, 0], [100.0, 200.0])
```
"""
struct Atoms{T<:AbstractFloat}
    types::Vector{Symbol}
    numbers::Vector{Int}
    masses::Vector{T}
end

function Atoms{T}(atom_types::AbstractVector{Symbol}) where {T}
    types = Vector{Symbol}(atom_types)
    numbers = Vector{Int}([element.number for element in elements[atom_types]])
    masses = Vector{T}([austrip(element.atomic_mass) for element in elements[atom_types]])
    Atoms{T}(types, numbers, masses)
end

function Atoms{T}(masses::AbstractVector) where {T}
    natoms = length(masses)
    types = Vector{Symbol}(fill(:X, natoms))
    numbers = Vector{Int}(zeros(natoms))
    masses = Vector{T}(austrip.(masses))
    Atoms{T}(types, numbers, masses)
end

Atoms{T}(atom_type) where {T} = Atoms{T}([atom_type])
Atoms(atom_types) = Atoms{Float64}(atom_types)

Base.length(atoms::Atoms) = length(atoms.types)
Base.range(atoms::Atoms) = range(1; length=length(atoms))

Base.IndexStyle(::Atoms) = IndexLinear()
Base.getindex(A::Atoms{T}, i::Int) where {T} = Atoms{T}(A.types[i])
Base.getindex(A::Atoms{T}, i::AbstractRange) where {T} = Atoms{T}(A.types[i])

function Base.:(==)(a::Atoms, b::Atoms)
    (a.types == b.types) && (a.numbers == b.numbers) && (a.masses ≈ b.masses)
end


text:
include("flight_enums.jl")

using Dates
using CSV
using Printf
import SQLite

mutable struct Flight <: Flighty
	id::UInt
	nr::UInt
	callsign::String
	aircraftType::String
	pilot::String
	copilot::Union{Missing, String}
	launch::Launch
	date::Date
	departureTime::Time
	arrivalTime::Time
	numberOfLandings::UInt
	departureLocation::String
	arrivalLocation::String

	comments::Union{Missing, String}

	dual::Bool
	instr::Bool
end

function Flight(row::NamedTuple)
	Flight(
		row.id,
		0,
		row.callsign,
		row.aircraftType,
		row.pilot,
		row.copilot,
		row.launch,
		Date(row.date, DateFormat("d.m.y")),
		Time(row.departureTime),
		Time(row.arrivalTime),
		row.numberOfLandings,
		row.departureLocation,
		row.arrivalLocation,
		row.comments,
		row.dual,
		row.instr
	)
end

function Flight(row::CSV.Row)
	inMyFlightbook = true
	!ismissing(row.Datum)		|| throw(ArgumentError("Datum nicht gegeben."))
	!ismissing(row.Lfz_) 		|| throw(ArgumentError("Kennzeichen nicht gegeben."))
	!ismissing(row.Lfz_Muster) 	|| throw(ArgumentError("Flugzeugtyp nicht gegeben."))
	!ismissing(row.Pilot)		|| throw(ArgumentError("Pilot nicht gegeben."))
	!ismissing(row.S_Art)		|| throw(ArgumentError("Startart nicht gegeben."))
	!ismissing(row.Start)		|| throw(ArgumentError("Startzeit nicht gegeben."))
	!ismissing(row.Landung)		|| throw(ArgumentError("Landezeit nicht gegeben."))
	!ismissing(row.Landungen)	|| throw(ArgumentError("Anzahl Landungen nicht gegeben."))
	!ismissing(row.Startort)	|| throw(ArgumentError("Startort nicht gegeben."))
	!ismissing(row.Landeort)	|| throw(ArgumentError("Landeort nicht gegeben."))
	!ismissing(row.Flugart)		|| throw(ArgumentError("Flugart nicht gegeben."))

	flightType = 	pm(strip, row.Flugart)
	pilotName = 	strip(row.Pilot)
	copilotName =	pm(strip, row.Begleiter_FI)

	typeNormalizations = [
		("GROB G 103 \"TWIN II\"", "Twin II"),
		("GROB G 103 C \"TWIN III\"", "Twin III"),
		("TWIN ASTIR TRAINER", "Twin I")
	]
	locNormalizations = [
		("Kammermark", "KM")
	]

	aircraftType = 	applyNormalizations(strip(row.Lfz_Muster), typeNormalizations)
	departureLoc = 	applyNormalizations(strip(row.Startort), locNormalizations)
	arrivalLoc =	applyNormalizations(strip(row.Landeort), locNormalizations)

	if lowercase(pilotName) in lowercase.(selfNames)
		pilot = true
	elseif !ismissing(copilotName) && lowercase(copilotName) in lowercase.(selfNames)
		pilot = false
	else
		throw(ArgumentError("Eigener Name weder im Feld 'Pilot' noch im Feld 'Copilot' gefunden."))
	end

	!ismissing(flightType) || throw(ArgumentError("Flugart nicht gegeben."))
	if !in(flightType, ["S", "N", "X", "G"])
		throw(ArgumentError("Ungültige Flugart '$flightType' vorgefunden."))
	end

	if !pilot
		# If copilot and training flight, then I was flight instructor
		if flightType == "S"
			instr = true
			dual = false
		# If copilot and no training flight, then not for my flightbook
		else
			instr = false
			dual = false
			inMyFlightbook = false
		end
	end

	if pilot
		# If pilot, then I was PIC or, if it is a training flight,
		# dual pilot (being instructed)
		instr = false
		dual = flightType == "S"
	end

	# If I was flight instructor or I was dual pilot (being instructed),
	# then pilot an copilot have to be swapped (compared to the usual
	# records made for the airfield)
	if (dual && pilot) || instr
		pilotName, copilotName = copilotName, pilotName
	end

	(
		Flight(
			0,
			0,
			row.Lfz_,
			aircraftType,
			pilotName,
			copilotName,
			Launch(row.S_Art),
			Date(row.Datum, DateFormat("d.m.y")),
			Time(row.Start),
			Time(row.Landung),
			row.Landungen,
			departureLoc,
			arrivalLoc,
			row.Bemerkung,
			dual,
			instr
		), inMyFlightbook
	)

end

function overlap(a::Flight, b::Flight)
	a.date == b.date &&
	a.arrivalTime > b.departureTime &&
	b.arrivalTime > a.departureTime
end

function fetchAllFlights(db::SQLite.DB)
	q = SQLite.Query(db, "SELECT * FROM flights")
	flights = [Flight(row) for row in q]

	sort!(flights)
	for i in 1:length(flights)
		flights[i].nr = i
	end

	return flights
end

# function printFlightHeader(aa::AutoAlign; showLocation=true, showComment=true)
# 	print(aa, "Nr.")
# 	print(aa, " ", "Datum")
# 	print(aa, "  ", "Reg.")
# 	print(aa, " ", "Typ")
# 	print(aa, " ", "S")
#
# 	print(aa, "  ", "Pilot")
# 	print(aa, "  ", "Co")
#
# 	if showLocation
# 		print(aa, "  ", "Von")
# 		print(aa, " ", "Nach")
# 	end
#
# 	print(aa, "  ", "Start")
# 	print(aa, "\u2E17", "Land.")
#
# 	print(aa, "  ", "Zeit")
#
# 	print(aa, " ", " ")
# 	print(aa, " ", "PIC")
# 	print(aa, " ", "DUAL")
# 	print(aa, " ", "INSTR")
#
# 	!showComment || print(aa, "  ", "Bemerkungen")
# 	println(aa)
# end

# function Base.print(aa::AutoAlign, f::Flight; showLocation=true, showComment=true)
# 	print(aa, f.nr)
# 	print(aa, " ", Dates.format(f.date, "dd.mm.yyyy"))
# 	print(aa, " ", f.callsign)
# 	print(aa, " ", f.aircraftType)
# 	print(aa, " ", convert(String, f.launch))
#
# 	print(aa, "  ", f.pilot)
# 	if ismissing(f.copilot)
# 		print(aa, "  ", "-")
# 	else
# 		print(aa, "  ", f.copilot)
# 	end
#
# 	if showLocation
# 		print(aa, "  ", f.departureLocation)
# 		print(aa, " ", f.arrivalLocation)
# 	end
#
# 	print(aa, "  ", Dates.format(f.departureTime, "HH:MM"))
# 	print(aa, "\u2E17", Dates.format(f.arrivalTime, "HH:MM"))
#
# 	timeString = asString(minutes(f))
# 	print(aa, "  ", timeString)
#
# 	flags = f.dual ? "D" : "P"
# 	flags *= f.instr ? "I" : ""
# 	print(aa, " ", flags)
#
# 	if f.dual
# 		print(aa, " ", " ")
# 		print(aa, " ", timeString)
# 	else
# 		print(aa, " ", timeString)
# 		print(aa, " ", " ")
# 	end
#
# 	if (f.instr)
# 		print(aa, " ", timeString)
# 	else
# 		print(aa, " ", " ")
# 	end
#
# 	if showComment
# 		print(aa, "  ", !ismissing(f.comments) ? f.comments : "")
# 	end
#
# 	println(aa)
# end

max_stars_repo_name: julianschick/flugbuch

max_issues_repo_name: julianschick/flugbuch

max_forks_repo_name: julianschick/flugbuch


text:
function predict_DNN_model_zscore(dataset::AbstractVector{Tuple{V1, V1, V2}}, model, ycol::Symbol,
    μ::AbstractFloat, σ::AbstractFloat,
    _eltype::DataType, output_size::Integer, output_dir::String) where
        {V1<:AbstractVector{<:Real}, V2<:AbstractVector{<:DateTime}}

    dnn_res_df = DataFrame(date = DateTime[], offset = Int[], y = _eltype[], yhat = _eltype[])

    for (x, y, dates) in dataset
        ŷ = model(x |> gpu)

        cpu_y = y |> cpu
        cpu_ŷ = ŷ |> cpu

        # 24 hour data
        org_y = unzscore(cpu_y, μ, σ)
        # 24 hour prediction
        org_ŷ = unzscore(cpu_ŷ, μ, σ)

        dnn_res_df_tmp = DataFrame(date = dates,
            offset = 1:output_size,
            y = org_y,
            yhat = org_ŷ)
        append!(dnn_res_df, dnn_res_df_tmp)
    end

    dnn_res_df
end

function predict_DNN_model_zscore(dataset::AbstractVector{Tuple{M1, M1, M2}}, model, ycol::Symbol,
    μ::AbstractFloat, σ::AbstractFloat,
    _eltype::DataType, output_size::Integer, output_dir::String) where
        {M1<:AbstractMatrix{<:Real}, M2<:AbstractMatrix{<:DateTime}}

    dnn_res_df = DataFrame(date = DateTime[], offset = Int[], y = _eltype[], yhat = _eltype[])

    for (x, y, dates) in dataset
        ŷ = model(x |> gpu)

        cpu_y = y |> cpu
        cpu_ŷ = ŷ |> cpu

        # 24 hour data
        org_y = unzscore(cpu_y, μ, σ)
        # 24 hour prediction
        org_ŷ = unzscore(cpu_ŷ, μ, σ)

        # size(x, 2) == size(y, 2) == size(dates, 2) == batch_size
        for i = 1:size(y, 2)
            dnn_res_df_tmp = DataFrame(date = dates,
                offset = 1:output_size,
                y = org_y[:, i],
                yhat = org_ŷ[:, i])
            append!(dnn_res_df, dnn_res_df_tmp)
        end
    end

    dnn_res_df
end

function predict_DNN_model_zscore_season(dataset::AbstractVector{T}, model, ycol::Symbol,
    μ::AbstractFloat, σ::AbstractFloat, season_table::AbstractNDSparse,
    _eltype::DataType, output_size::Integer, output_dir::String) where T <: Tuple

    dnn_res_df = DataFrame(date = DateTime[], offset = Int[], y = _eltype[], yhat = _eltype[])

    for (x, y, dates) in dataset
        ŷ = model(x |> gpu)

        cpu_y = y |> cpu
        cpu_ŷ = ŷ |> cpu

        # 24 hour data
        org_y = compose_seasonality(collect(DateTime.(dates, Local)),
            unzscore(cpu_y, μ, σ),
            season_table)

        # 24 hour prediction
        org_ŷ = compose_seasonality(collect(DateTime.(dates, Local)),
            unzscore(cpu_ŷ, μ, σ),
            season_table)

        dnn_res_df_tmp = DataFrame(date = DateTime.(dates, Local),
            offset = 1:output_size,
            y = org_y,
            yhat = org_ŷ)
        append!(dnn_res_df, dnn_res_df_tmp)
    end

    dnn_res_df
end

function predict_DNN_model_minmax(dataset::AbstractVector{T}, model, ycol::Symbol,
    _min::AbstractFloat, _max::AbstractFloat, a::AbstractFloat, b::AbstractFloat,
    _eltype::DataType, output_size::Integer, output_dir::String) where T <: Tuple

    dnn_res_df = DataFrame(date = DateTime[], offset = Int[], y = _eltype[], yhat = _eltype[])

    # https://en.wikipedia.org/wiki/Feature_scaling#Rescaling_(min-max_normalization)
    c = ((_max - _min) / (b - a))
    for (x, y, dates) in dataset
        ŷ = model(x |> gpu)

        cpu_y = y |> cpu
        cpu_ŷ = ŷ |> cpu

        # 24 hour data
        org_y = unminmax_scaling(cpu_y, _min, _max, a, b)
        # 24 hour prediction
        org_ŷ = unminmax_scaling(cpu_ŷ, _min, _max, a, b)

        dnn_res_df_tmp = DataFrame(date = DateTime.(dates, Local),
            offset = 1:output_size,
            y = org_y,
            yhat = org_ŷ)
        append!(dnn_res_df, dnn_res_df_tmp)
    end

    dnn_res_df
end

function predict_DNN_model_logzscore(dataset::AbstractVector{T}, model, ycol::Symbol,
    μ::AbstractFloat, σ::AbstractFloat,
    _eltype::DataType, output_size::Integer, output_dir::String) where T <: Tuple

    dnn_res_df = DataFrame(date = DateTime[], offset = Int[], y = _eltype[], yhat = _eltype[])

    for (x, y, dates) in dataset
        ŷ = model(x |> gpu)

        cpu_y = y |> cpu
        cpu_ŷ = ŷ |> cpu

        # 24 hour data
        org_y = exp.(unzscore(cpu_y, μ, σ)) .- 10.0
        # 24 hour prediction
        org_ŷ = exp.(unzscore(cpu_ŷ, μ, σ)) .- 10.0

        dnn_res_df_tmp = DataFrame(date = DateTime.(dates, Local),
            offset = 1:output_size,
            y = org_y,
            yhat = org_ŷ)
        append!(dnn_res_df, dnn_res_df_tmp)
    end

    dnn_res_df
end

function predict_DNN_model_invzscore(dataset::AbstractVector{T}, model, ycol::Symbol,
    μ::AbstractFloat, σ::AbstractFloat,
    _eltype::DataType, output_size::Integer, output_dir::String) where T <: Tuple

    dnn_res_df = DataFrame(date = DateTime[], offset = Int[], y = _eltype[], yhat = _eltype[])

    for (x, y, dates) in dataset
        ŷ = model(x |> gpu)

        cpu_y = y |> cpu
        cpu_ŷ = ŷ |> cpu

        # 24 hour data
        org_y = 1.0 ./ (unzscore(cpu_y, μ, σ)) .- 10.0
        # 24 hour prediction
        org_ŷ = 1.0 ./ (unzscore(cpu_ŷ, μ, σ)) .- 10.0

        dnn_res_df_tmp = DataFrame(date = DateTime.(dates, Local),
            offset = 1:output_size,
            y = org_y,
            yhat = org_ŷ)
        append!(dnn_res_df, dnn_res_df_tmp)
    end

    dnn_res_df
end

function predict_DNN_model_logminmax(dataset::AbstractVector{T}, model, ycol::Symbol,
    _min::AbstractFloat, _max::AbstractFloat, a::AbstractFloat, b::AbstractFloat,
     _eltype::DataType, output_size::Integer, output_dir::String) where T <: Tuple

    dnn_res_df = DataFrame(date = DateTime[], offset = Int[], y = _eltype[], yhat = _eltype[])

    # https://en.wikipedia.org/wiki/Feature_scaling#Rescaling_(min-max_normalization)
    c = ((_max - _min) / (b - a))
    for (x, y, dates) in dataset
        ŷ = model(x |> gpu)

        cpu_y = y |> cpu
        cpu_ŷ = ŷ |> cpu

        # 24 hour data
        org_y = exp.(unminmax_scaling(cpu_y, _min, _max, a, b)) .- 10.0
        # 24 hour prediction
        org_ŷ = exp.(unminmax_scaling(cpu_ŷ, _min, _max, a, b)) .- 10.0

        dnn_res_df_tmp = DataFrame(date = DateTime.(dates, Local),
            offset = 1:output_size,
            y = org_y,
            yhat = org_ŷ)
        append!(dnn_res_df, dnn_res_df_tmp)
    end

    dnn_res_df
end

function predict_RNN_model_zscore(dataset::AbstractVector{Tuple{F4, AF2, F2, D1}}, model, ycol::Symbol,
    μ::AbstractFloat, σ::AbstractFloat,
     _eltype::DataType, output_size::Integer, output_dir::String) where
        {F4<:AbstractArray{<:Real, 4}, F2<:AbstractArray{<:Real, 2},
        AF2<:AbstractVector{<:AbstractArray{<:Real, 2}}, D1<:AbstractVector{<:DateTime}}
    # pair test set (xe, xd, y) is same as batch with batch_size == 1, but `dates` are Vector
    rnn_res_df = DataFrame(date = DateTime[], offset = Int[], y = _eltype[], yhat = _eltype[])

    for (xe, xd, y, dates) in dataset
        ŷ = model(xe |> gpu, xd |> gpu)

        cpu_y = y[:, 1] |> cpu
        cpu_ŷ = ŷ[:, 1] |> cpu

        # 24 hour data
        org_y = unzscore(cpu_y, μ, σ)
        # 24 hour prediction
        org_ŷ = unzscore(cpu_ŷ, μ, σ)

        rnn_res_df_tmp = DataFrame(date = dates,
            offset = 1:output_size,
            y = org_y,
            yhat = org_ŷ)
        append!(rnn_res_df, rnn_res_df_tmp)
    end

    rnn_res_df
end

function predict_RNN_model_zscore(dataset::AbstractVector{Tuple{F4, AF2, F2, D2}}, model, ycol::Symbol,
    μ::AbstractFloat, σ::AbstractFloat,
     _eltype::DataType, output_size::Integer, output_dir::String) where
        {F4<:AbstractArray{<:Real, 4}, F2<:AbstractArray{<:Real, 2},
        AF2<:AbstractVector{<:AbstractArray{<:Real, 2}}, D2<:AbstractMatrix{<:DateTime}}

    # batch test set, only difference with pair test set is type of `dates`
    rnn_res_df = DataFrame(date = DateTime[], offset = Int[], y = _eltype[], yhat = _eltype[])

    for (xe, xd, y, dates) in dataset
        ŷ = model(xe |> gpu, xd |> gpu)

        cpu_y = y |> cpu
        cpu_ŷ = ŷ |> cpu

        # size(x, 2) == size(y, 2) == size(dates, 2) == batch_size
        for i = 1:size(y, 2)
            # 24 hour data
            org_y = unzscore(cpu_y, μ, σ)
            # 24 hour prediction
            org_ŷ = unzscore(cpu_ŷ, μ, σ)

            rnn_res_df_tmp = DataFrame(date = dates[:, i],
                offset = 1:output_size,
                y = org_y[:, i],
                yhat = org_ŷ[:, i])
            append!(rnn_res_df, rnn_res_df_tmp)
        end
    end

    rnn_res_df
end

function predict_RNN_model_minmax(dataset::AbstractVector{Tuple{V1, V1, V1, V2}}, model, ycol::Symbol,
    _min::AbstractFloat, _max::AbstractFloat, a::AbstractFloat, b::AbstractFloat,
    _eltype::DataType, output_size::Integer, output_dir::String) where
        {V1<:AbstractVector{<:Real}, V2<:AbstractVector{<:DateTime}}

    rnn_res_df = DataFrame(date = DateTime[], offset = Int[], y = _eltype[], yhat = _eltype[])

    # https://en.wikipedia.org/wiki/Feature_scaling#Rescaling_(min-max_normalization)
    c = ((_max - _min) / (b - a))
    for (xe, xd, y, dates) in dataset
        ŷ = model(xe |> gpu, xd |> gpu)

        cpu_y = y[:, 1] |> cpu
        cpu_ŷ = ŷ[:, 1] |> cpu

        # 24 hour data
        org_y = unminmax_scaling(cpu_y, _min, _max, a, b)
        # 24 hour prediction
        org_ŷ = unminmax_scaling(cpu_ŷ, _min, _max, a, b)

        rnn_res_df_tmp = DataFrame(date = dates,
            offset = 1:output_size,
            y = org_y,
            yhat = org_ŷ)
        append!(rnn_res_df, rnn_res_df_tmp)
    end

    rnn_res_df
end

function predict_RNN_model_minmax(dataset::AbstractVector{Tuple{M1, M1, M1, M2}}, model, ycol::Symbol,
    _min::AbstractFloat, _max::AbstractFloat, a::AbstractFloat, b::AbstractFloat,
    _eltype::DataType, output_size::Integer, output_dir::String) where
        {M1<:AbstractMatrix{<:Real}, M2<:AbstractMatrix{<:DateTime}}

    rnn_res_df = DataFrame(date = DateTime[], offset = Int[], y = _eltype[], yhat = _eltype[])

    # https://en.wikipedia.org/wiki/Feature_scaling#Rescaling_(min-max_normalization)
    c = ((_max - _min) / (b - a))
    for (xe, xd, y, dates) in dataset
        ŷ = model(xe |> gpu, xd |> gpu)

        cpu_y = y[:, 1] |> cpu
        cpu_ŷ = ŷ[:, 1] |> cpu

        # 24 hour data
        org_y = unminmax_scaling(cpu_y, _min, _max, a, b)
        # 24 hour prediction
        org_ŷ = unminmax_scaling(cpu_ŷ, _min, _max, a, b)

        # size(x, 2) == size(y, 2) == size(dates, 2) == batch_size
        for i = 1:size(x, 2)
            rnn_res_df_tmp = DataFrame(date = dates[:, i],
                offset = 1:output_size,
                y = org_y[:, i],
                yhat = org_ŷ[:, i])
            append!(rnn_res_df, rnn_res_df_tmp)
        end
    end

    rnn_res_df
end

"""
    export_CSV

Basic NDSparse
"""
function export_CSV(dates::Array{DateTime, 1}, dnn_table::Array{IndexedTable, 1},
    ycol::Symbol, output_size::Integer,
    _eltype::DataType, output_dir::String, output_prefix::String)

    dfs = Array{DataFrame}(undef, output_size)

    for i = 1:output_size
        i_pad = lpad(i, 2, '0')
        Base.Filesystem.mkpath(output_dir * "$(i_pad)/")
        plottable_path::String = output_dir * "$(i_pad)/" * "$(output_prefix)_plottable_$(i_pad)h.csv"

        dates_h = dates .+ Dates.Hour(i)
        y = JuliaDB.select(dnn_table[i], :y)
        ŷ = JuliaDB.select(dnn_table[i], :ŷ)

        len_data = min(length(dates), length(y), length(ŷ))

        df = DataFrame(
            date = dates_h[1:len_data], y = y[1:len_data], yhat = ŷ[1:len_data])

        dfs[i] = df

        CSV.write(plottable_path, df)
    end

    dfs
end

"""
    export_CSV

with seasonality decomposition
"""
function export_CSV(dates::Array{DateTime, 1}, table::AbstractNDSparse, season_table::AbstractNDSparse,
    ycol::Symbol, output_size::Integer, output_dir::String, output_prefix::String)

    dfs = Array{DataFrame}(undef, output_size)

    for i = 1:output_size
        i_pad = lpad(i, 2, '0')
        Base.Filesystem.mkpath(output_dir * "$(i_pad)/")
        plotdf_path::String = output_dir * "$(i_pad)/" * "$(output_prefix)_plottable_$(i_pad)h.csv"

        dates_h = dates .+ Dates.Hour(i)
        y = JuliaDB.select(table[i], :y)
        ŷ = compose_seasonality(dates_h, JuliaDB.select(table[i], :ŷ), season_table)

        len_data = min(length(dates), length(y), length(ŷ))

        df = DataFrame(
            date = dates_h[1:len_data], y = y[1:len_data], yhat = ŷ[1:len_data])

        sort!(df, (:date))
        dfs[i] = df

        CSV.write(plotdf_path, df)
    end

    dfs
end

"""
    export_CSV

new DataFrame
"""
function export_CSV(dates::Array{DateTime, 1}, df::DataFrame,
    ycol::Symbol, output_size::Integer, output_dir::String, output_prefix::String)

    offset_dfs = Array{DataFrame}(undef, output_size)

    for i = 1:output_size
        i_pad = lpad(i, 2, '0')
        Base.Filesystem.mkpath(output_dir * "$(i_pad)/")
        plotdf_path::String = output_dir * "$(i_pad)/" * "$(output_prefix)_plottable_$(i_pad)h.csv"

        dates_h = dates .+ Dates.Hour(i)
        offset_df = filter(row -> row[:offset] == i, df)

        sort!(offset_df, (:date))
        offset_dfs[i] = offset_df

        CSV.write(plotdf_path, offset_df)
    end

    offset_dfs
end

"""
    export_CSV

new DataFrame + seasonality
"""
function export_CSV(dates::Array{DateTime, 1}, df::DataFrame, season_table::AbstractNDSparse,
    ycol::Symbol, output_size::Integer, output_dir::String, output_prefix::String)

    offset_dfs = Array{DataFrame}(undef, output_size)

    for i = 1:output_size
        i_pad = lpad(i, 2, '0')
        Base.Filesystem.mkpath(output_dir * "$(i_pad)/")
        plotdf_path::String = output_dir * "$(i_pad)/" * "$(output_prefix)_plottable_$(i_pad)h.csv"

        dates_h = dates .+ Dates.Hour(i)
        offset_df = filter(row -> row[:offset] == i, df)

        sort!(offset_df, (:date))
        offset_dfs[i] = offset_df

        CSV.write(plotdf_path, offset_df)
    end

    offset_dfs
end

function compute_corr(dnn_table::Array{IndexedTable, 1},
    output_size::Integer, output_dir::String, output_prefix::String)

    corr_path = output_dir * output_prefix * "_corr.csv"
    corr = zeros(output_size)

    for i = 1:output_size
        corr[i] = Statistics.cor(JuliaDB.select(dnn_table[i], :y), JuliaDB.select(dnn_table[i], :ŷ))
    end

    df = DataFrame(hour = collect(1:output_size), corr = corr)

    CSV.write(corr_path, df)

    df
end

function compute_corr(dfs::Array{DataFrame, 1},
    output_size::Integer, output_dir::String, output_prefix::String)

    corr_path = output_dir * output_prefix * "_corr.csv"
    corr = zeros(output_size)

    for i = 1:output_size
        _df = dfs[i]
        corr[i] = Statistics.cor(_df[!, :y], _df[!, :yhat])
    end

    df = DataFrame(hour = collect(1:output_size), corr = corr)

    CSV.write(corr_path, df)

    df
end


text:
using tinygrad

Base.broadcastable(m::Diffable) = Ref(m) #treat diffable as scalar if Diffable not <: Number

function descend(n = 100, α = 1e-1, epochs = 100)

    f(x) = 5x + 4 #true
    y(x) = f(x) + randn() * 0.5 #noised

    ŷ(m,b,x) = m*x + b # estimate

    # MSE loss
    loss(y,ŷ) = abs.(y .- ŷ) .^ 2 |> x-> sum(x) / n # extra abs resolves sign issues on extraneous Diffable promotion

    # make the training data (with pre-normalized inputs)
    xs = range(-1.0; stop = 1.0, length = n) |> collect
    ys = y.(xs)

    m = abs(randn()) |> Diffable
    b = 0 |> Diffable

    ms = Vector{Float64}([m.value])
    bs = Vector{Float64}([b.value])

    for epoch in 1:epochs
        ŷ(x) = ŷ(m,b,x)

        ŷs = ŷ.(xs)

        currentloss = loss(ys, ŷs)
        backwards!(m)
        backwards!(b)

        m -= m.grad * α; clear!(m)
        b -= b.grad * α; clear!(b)

        push!(ms, m.value)
        push!(bs, b.value)

        epoch % floor(epochs/10) == 0 ? println("Loss on iteration $epoch: $(currentloss.value)") : nothing
    end

    @info "m: $(m.value); b: $(b.value)"
    return ms, bs, xs, ys
end

ms, bs, xs, ys = descend()

using Plots; gr()

# create a scatter plot for each set of values
function create_plot(m,b,iter)

    ŷs = [m * x + b for x in xs]
    plot(xs, ŷs, w = 3)
    scatter!(xs, ys, markerstrokewidth = 0, legend = nothing, title = "Iteration $iter")
end

animation = @animate for i in 1:40
    create_plot(ms[i], bs[i], i)
end

gif(animation, "linear_fit.gif", fps = 3)


text:
using Alea, FactCheck

facts("OrthPoly") do
  context("ChebyshevTPolynomials") do
    polys = ChebyshevTPolynomials()
    cheb4(x) = [1, x, 2x.^2 - 1, 4x.^3 - 3x, 8x.^4 - 8x.^2 + 1]
    @fact evaluate(polys, 4, 2) --> cheb4(2)
    @fact recurrence_coeff(polys, 3) --> (0, 2, 1)
    @pending norm(polys, 3) --> 1
  end
end

facts("GaussIntegration") do
  @fact gauss_rule(HermitePolynomials(), 1) --> ([0.0], [1.0])
  @fact gauss_rule(HermitePolynomials(), 2) --> ([-1.0, 1.0], [0.5, 0.5])
  @fact gauss_rule(HermitePolynomials(), 3) --> roughly(([-√3, 0, √3], [1/6, 4/6, 1/6]))

  @fact gauss_rule(LegendrePolynomials(), 2) --> ([-√(1/3), √(1/3)], [0.5, 0.5])
  @fact gauss_rule(LegendrePolynomials(), 3) --> roughly(([-√(3/5), 0, √(3/5)], [5, 8, 5]/18))

  @fact gauss_rule(LaguerrePolynomials(0.0), 2) --> roughly(([2-√2, 2+√2],
                                                           [1/2+1/√8, 1/2-1/√8]))
end


text:
export @metal

"""
    @metal [kwargs...] func(args...)

High-level interface for executing code on a GPU. The `@metal` macro should prefix a call,
with `func` a callable function or object that should return nothing. It will be compiled to
a Metal function upon first use, and to a certain extent arguments will be converted and
managed automatically using `mtlconvert`. Finally, a call to `mtlcall` is
performed, creating a command buffer in the current global command queue then committing it.

There is one supported keyword argument that influences the behavior of `@metal`.
- `launch`: whether to launch this kernel, defaults to `true`. If `false` the returned
  kernel object should be launched by calling it and passing arguments again.
"""
macro metal(ex...)
    call = ex[end]
    kwargs = ex[1:end-1]

    # destructure the kernel call
    Meta.isexpr(call, :call) || throw(ArgumentError("second argument to @metal should be a function call"))
    f = call.args[1]
    args = call.args[2:end]

    code = quote end
    vars, var_exprs = assign_args!(code, args)

    # group keyword argument
    macro_kwargs, compiler_kwargs, call_kwargs, other_kwargs =
        split_kwargs(kwargs, [:launch], [:name], [:grid, :threads])
    if !isempty(other_kwargs)
        key,val = first(other_kwargs).args
        throw(ArgumentError("Unsupported keyword argument '$key'"))
    end

    # handle keyword arguments that influence the macro's behavior
    launch = true
    for kwarg in macro_kwargs
        key,val = kwarg.args
        if key == :launch
            isa(val, Bool) || throw(ArgumentError("`launch` keyword argument to @metal should be a Bool"))
            launch = val::Bool
        else
            throw(ArgumentError("Unsupported keyword argument '$key'"))
        end
    end
    if !launch && !isempty(call_kwargs)
        error("@metal with launch=false does not support launch-time keyword arguments; use them when calling the kernel")
    end

    # FIXME: macro hygiene wrt. escaping kwarg values (this broke with 1.5)
    #        we esc() the whole thing now, necessitating gensyms...
    @gensym f_var kernel_f kernel_args kernel_tt kernel

    # convert the arguments, call the compiler and launch the kernel
    # while keeping the original arguments alive
    push!(code.args,
        quote
            $f_var = $f
            GC.@preserve $(vars...) $f_var begin
                local $kernel_f = $mtlconvert($f_var)
                local $kernel_args = map($mtlconvert, ($(var_exprs...),))
                local $kernel_tt = Tuple{map(Core.Typeof, $kernel_args)...}
                local $kernel = $mtlfunction($kernel_f, $kernel_tt; $(compiler_kwargs...))
                if $launch
                    $kernel($(var_exprs...); $(call_kwargs...))
                end
                $kernel
            end
         end)
    return esc(code)
end

## argument conversion

struct Adaptor end

# convert Metal Buffers Metal device pointers # should be generic
Adapt.adapt_storage(to::Adaptor, buf::MtlBuffer{T}) where {T} = reinterpret(Core.LLVMPtr{T,AS.Device}, buf.handle)

# Base.RefValue isn't GPU compatible, so provide a compatible alternative
struct MtlRefValue{T} <: Ref{T}
  x::T
end
Base.getindex(r::MtlRefValue) = r.x
Adapt.adapt_structure(to::Adaptor, r::Base.RefValue) = MtlRefValue(adapt(to, r[]))

"""
  mtlconvert(x)

This function is called for every argument to be passed to a kernel, allowing it to be
converted to a GPU-friendly format. By default, the function does nothing and returns the
input object `x` as-is.

Do not add methods to this function, but instead extend the underlying Adapt.jl package and
register methods for the the `Metal.Adaptor` type.
"""
mtlconvert(arg) = adapt(Adaptor(), arg)

## host-side kernel API

struct HostKernel{F,TT}
    f::F
    fun::MtlFunction
    arg_info::Any
end

"""
    mtlfunction(f, tt=Tuple{}; kwargs...)

Low-level interface to compile a function invocation for the currently-active GPU, returning
a callable kernel object. For a higher-level interface, use [`@metal`](@ref).

The output of this function is automatically cached, i.e. you can simply call `mtlfunction`
in a hot path without degrading performance. New code will be generated automatically when
the function changes, or when different types or keyword arguments are provided.
"""
function mtlfunction(f::F, tt::TT=Tuple{}; name=nothing, kwargs...) where {F,TT}
    dev = MtlDevice(1)
    cache = get!(()->Dict{UInt,Any}(), mtlfunction_cache, dev)
    source = FunctionSpec(f, tt, true, name)
    target = MetalCompilerTarget(macos=macos_version(); kwargs...)
    params = MetalCompilerParams()
    job = CompilerJob(target, source, params)
    fun, arguments = GPUCompiler.cached_compilation(cache, job,
                                                    mtlfunction_compile, mtlfunction_link)
    # compilation is cached on the function type, so we can only create a kernel object here
    # (as it captures the function _instance_). we may want to cache those objects.
    HostKernel{F,tt}(f, fun, arguments)
end

const mtlfunction_cache = Dict{Any,Any}()

function mtlfunction_compile(@nospecialize(job::CompilerJob))
    # TODO: on 1.9, this actually creates a context. cache those.
    JuliaContext() do ctx
        mtlfunction_compile(job, ctx)
    end
end
function mtlfunction_compile(@nospecialize(job::CompilerJob), ctx::Context)
    mi, mi_meta = GPUCompiler.emit_julia(job)
    ir, ir_meta = GPUCompiler.emit_llvm(job, mi; ctx)
    image, asm_meta = GPUCompiler.emit_asm(job, ir; strip=false, format=LLVM.API.LLVMObjectFile)
    # TODO: don't strip IR
    return (image, entry=LLVM.name(ir_meta.entry), arguments=job.meta[:arguments])
end

function mtlfunction_link(@nospecialize(job::CompilerJob), compiled)
    dev = current_device()
    lib = MtlLibraryFromData(dev, compiled.image)
    MtlFunction(lib, compiled.entry), compiled.arguments
end

## kernel launching and argument encoding

function (kernel::HostKernel)(args...; grid::MtlDim=1, threads::MtlDim=1)
    grid = MtlDim3(grid)
    threads = MtlDim3(threads)
    (grid.width>0 && grid.height>0 && grid.depth>0) ||
        throw(ArgumentError("Grid dimensions should be non-null"))
    (threads.width>0 && threads.height>0 && threads.depth>0) ||
        throw(ArgumentError("Threadgroup dimensions should be non-null"))

    pipeline_state = MtlComputePipelineState(kernel.fun.lib.device, kernel.fun)
    (threads.width * threads.height * threads.depth) > pipeline_state.maxTotalThreadsPerThreadgroup &&
        throw(ArgumentError("Max total threadgroup size should not exceed $(pipeline_state.maxTotalThreadsPerThreadgroup)"))

    cmdq = global_queue(kernel.fun.lib.device)
    cmdbuf = MtlCommandBuffer(cmdq)
    MtlComputeCommandEncoder(cmdbuf) do cce
        MTL.set_function!(cce, pipeline_state)

        # encode arguments
        idx = 1
        args = map(mtlconvert, (kernel.f, args...))
        @assert length(args) == length(kernel.arg_info)
        for (arg, arg_info) in zip(args, kernel.arg_info)
            arg_info.kind == GPUCompiler.GhostArgument && continue

            if arg_info.kind == GPUCompiler.BufferArgument
                if arg isa Core.LLVMPtr
                    buf = MtlBuffer{Nothing}(Base.bitcast(MTL.MTLBuffer, arg))
                    set_buffer!(cce, buf, 0, idx)
                else
                    @assert isbits(arg)
                    ref = Base.RefValue(arg)
                    GC.@preserve ref begin
                        ptr = Base.unsafe_convert(Ptr{Nothing}, ref)
                        set_bytes!(cce, ptr, sizeof(ref), idx)
                    end
                end
            elseif arg_info.kind == GPUCompiler.ArrayArgument
                @assert isbits(arg)
                ref = Base.RefValue(arg)
                GC.@preserve ref begin
                    ptr = Base.unsafe_convert(Ptr{Nothing}, ref)
                    set_bytes!(cce, ptr, sizeof(ref), idx)
                end
            elseif arg_info.kind == GPUCompiler.StructArgument
                @assert isbits(arg)
                ref = Base.RefValue(arg)
                GC.@preserve ref begin
                    ptr = Base.unsafe_convert(Ptr{Nothing}, ref)
                    set_bytes!(cce, ptr, sizeof(ref), idx)
                end
            elseif arg_info.kind == GPUCompiler.IndirectStructArgument
                # create an argument encoder
                arg_enc = MtlArgumentEncoder(kernel.fun, idx)
                arg_buf = alloc(Cchar, kernel.fun.lib.device, sizeof(arg_enc), storage=Shared)
                MTL.assign_argument_buffer!(arg_enc, arg_buf, 0)

                # encode fields
                function encode_field(arg, arg_info)
                    @assert fieldcount(typeof(arg)) == length(arg_info.fields)
                    for (field_name, field_info) in zip(fieldnames(typeof(arg)), arg_info.fields)
                        field_info.kind == GPUCompiler.GhostArgument && continue
                        field = getfield(arg, field_name)
                        if field_info.kind == GPUCompiler.BufferArgument
                            @assert field isa Core.LLVMPtr
                            buf = MtlBuffer{Nothing}(Base.bitcast(MTL.MTLBuffer, field))
                            set_buffer!(arg_enc, buf, 0, field_info.id)
                            MTL.use!(cce, buf, MTL.ReadWriteUsage)
                        elseif field_info.kind == GPUCompiler.ConstantArgument
                            set_constant!(arg_enc, field, field_info.id)
                        elseif field_info.kind == GPUCompiler.IndirectStructArgument
                            encode_field(field, field_info)
                        elseif field_info.kind == GPUCompiler.ArrayArgument
                            # XXX: this seems weird?
                            set_constant!(arg_enc, field, field_info.id)
                        else
                            error("Unknown struct field: $(field_info.kind) $(field)")
                        end
                    end
                end
                encode_field(arg, arg_info)

                # encode argument
                set_buffer!(cce, arg_buf, 0, idx)
            else
                error("Unknown argument kind: $(arg_info.kind)")
            end

            idx += 1
        end

        MTL.append_current_function!(cce, grid, threads)
    end
    commit!(cmdbuf)
end


text:
module Turnover

    using DataFrames
    using TumorGrowth: clones_by_mutations

    export mfreqs, reduced_μ, reduced_μ!, remove_mutations, remove_mutations!
    export W_orphaned, W_estranged
    export estranged, estranged_treeless, estranged_expected
    export orphaned_red, orphaned_green, orphaned_red_treeless, orphaned_green_treeless, orphaned_red_expected, orphaned_green_expected

    ##########################
    ### Neutral Simulation ###
    ##########################

    include("neutral_haplotype_growth_v2.jl")
#     include("neutral_haplotype_growth_v2_removedead.jl")
    include("neutral_haplotype_growth_clones_v2.jl")

    """
        remove_mutations!(htumor, keep)
    Removes all mutations that are not in the collection `keep` and deletes their haplotypes from the haplotype tumor. Remaining identical haplotypes are recombined. Mutation indices will generally not coincide with row indices anymore!
    """
    function remove_mutations!(htumor, keep)
        removed = setdiff(1:nrow(htumor)-1, keep)
        for (i,muts) in enumerate(htumor.mutations)
            filter!(in(keep), muts)
            n = htumor.n[i]
            htumor.n[i] = 0
            if isempty(muts)
                htumor.n[1] += n
            else
                htumor.n[last(muts)+1] += n
            end
        end
        delete!(htumor, removed.+1)
        return htumor
    end
    remove_mutations(htumor, keep) = remove_mutations!(deepcopy(htumor), keep)

    """
        reduced_μ!(htumor, x)
    Artificially reduces mutation rate of a **complete (!) haplotype-tumor** by sampling from its set of mutations, where each mutation is selected with probability x.

    Remaining identical haplotypes are recombined and mutation indices mapped to `1:nrow(htumor)-1`.
    """
    function reduced_μ!(htumor, x)
        keep = findall(rand(nrow(htumor)-1).<=x)
        remove_mutations!(htumor, keep)
        rebrand = (keep .=> 1:nrow(htumor)-1)
        for muts in htumor.mutations
            replace!(muts, rebrand...)
        end
        return htumor
    end
    reduced_μ(htumor, x) = reduced_μ!(deepcopy(htumor), x)

    ##########################
    ### Theoretical result ###
    ##########################

    """
        W_orphaned(q; N)
    Analytical solution to (red) **orphaned turnover**. Under this definition the probability of being orphaned
    depends on mutation extinction probability `q = d/b` and (threshold-) size `N` only.
    """
    W_orphaned(q; N) = q / (2*log(N)) * (1-N^(-2))

    """
        W_estranged(q; N)
    Analytical solution to **estranged turnover**. A haplotype clone's probability of becoming estranged, i.e. losing it's haplotype parent,
    depends on birth rate `b`, death rate `d`, mutation rate `μ` and (threshold-) time `T`.

    Assumes **Branchingprocess `v2`**: independent mutation of offspring pair at division. `μ` is the probability of mutation per cell division.
    The mutation probability per cell division per cell is `μ/2`.
    """
    W_estranged(d; b, μ, T) = (μ/2)*(d + b*(μ/2)^2) / (1-(μ/2))^2 / (b*(1-μ/2)-d) * (1-exp(-2*(b-d-b*μ/2)*T)) / (1-exp(-2*b*μ/2*T))

    #########################
    ### METHODS WITH TREE ###
    #########################

    """
        mfreqs(htumor)

    Get mutation frequencies for a **complete (!) haplotype-tumor**,
    `complete` meaning that haplotypes with size `n` = 0 are listed.
    """
    function mfreqs(htumor)
        count = zeros(Int, nrow(htumor)-1)
        for (i, muts) in enumerate(htumor.mutations[2:end])
            count[muts] .+= htumor.n[i+1]
        end
        count ./ sum(htumor.n)
    end

    """
        estranged(htumor)

    Estranged turnover on **complete (!) haplotype-tumors**,
    *complete* meaning that haplotypes with size *n* = 0 are listed.
    The function reports for all offspring mutations wether their parental
    haplotype went extinct in form of

        DataFrame( mutation :: Vector{Int}, isestranged :: Vector{Bool})
    """
    function estranged(htumor)
        # freqs = mfreqs(htumor)
        # daughtertypes = htumor.mutations[ [false, (freqs .> 0)...] .& (length.(htumor.mutations) .> 1)  ]

        daughtertypes = htumor.mutations[ (htumor.n .> 0) .& (length.(htumor.mutations) .> 1)  ]

        parents_n = [htumor.n[muts[end-1]+1] for muts=daughtertypes]
        return DataFrame(  mutation = last.(daughtertypes), isestranged = iszero.(parents_n))
    end

    ### red-clones

    """
        orphaned_red(htumor)

    Orphaned turnover based on parent clones on **complete (!) haplotype-tumors**,
    `complete` meaning that haplotypes with size `n` = 0 are listed.
    The function returns the orphan-status of all offspring mutations of every clone as

        DataFrame( mutation :: Vector{Int}, isorphaned :: Vector{Bool})

    This result reduces to the method based on offspring clones by counting unique
    mutations in `mutation` to get total offspring count and by counting a mutation as
    orphaned if it is orphaned with respect to any of the parental clones

        total count: orphaned_red(htumor).mutation |> unique |> length
        orphan count: unique(orphaned_red(htumor)).orphaned |> count

    """
    function orphaned_red(htumor)
        mtypes = clones_by_mutations(htumor)[1]
        freqs = sum.(getproperty.(mtypes,:n))

        orphans = DataFrame(mutation = Int[], isorphaned = Bool[])
        for (m,redclone) in enumerate(mtypes)
            n_clone = freqs[m]
            m_green = last.(redclone.mutations[Not(1)])
            filter!( m-> !iszero(freqs[m]), m_green )
            n_green = freqs[m_green]
            isempty(n_green) && continue

            append!(orphans.mutation, m_green)
            append!(orphans.isorphaned, iszero.(n_clone .- n_green))
        end
        return orphans
    end

    # To recover "green-clones" method:
    #     - total count: orphaned_red(htumor).mutation |> unique |> length
    #     - orphan count: unique(orphaned_red(htumor)).orphaned |> count

    ### green-clones
    """
        orphaned_green(htumor)

    Orphaned turnover based on offspring clones on **complete (!) haplotype-tumors**,
    `complete` meaning that haplotypes with size `n` = 0 are listed.
    The function returns the orphan-status of all offspring mutations as

        DataFrame( mutation :: Vector{Int}, isorphaned :: Vector{Bool})
    """
    function orphaned_green(tumor)
        mtypes, mutations = clones_by_mutations(tumor)
        n = sum.(getproperty.(mtypes,:n))

        typeone = tumor.mutations[Not(1)] .|> first |> unique
        m_green = filter(m -> !iszero(n[m]), setdiff(mutations, typeone))
        m_red = [ first(mtypes[m].mutations)[end-1] for m=m_green]

        countorph = count(iszero, n[m_red] .- n[m_green])

        orphans = DataFrame(mutation = m_green, isorphaned = iszero.(n[m_red] .- n[m_green]) )
        return orphans
    end

    #########################
    ### TREE-LESS METHODS ###
    #########################

    """
        estranged_treeless(htumor)

    Tree-less version of the `estranged` algorithm on **incomplete (!) haplotype-tumors**. In an incomplete
    haplotype-tumor no haplotype has size `n` = 0 (those are to be removed beforehand).
    The function returns

        DataFrame( mutation :: Vector{Int}, isestranged :: Vector{Float64}, isgreen :: Vector{Float64} )

    Summing over `isestranged` gives the number of estranged mutations and summing over `isgreen` the total number
    of offspring haplotypes.
    """
    function estranged_treeless(htumor)
        # htumor = filter(h -> !iszero(h.n), tumor)

        mtypes, mutations = clones_by_mutations(htumor)

        truncs = [ intersect(Set.(mtype.mutations)...) for mtype in mtypes]
        htypes = Set.(htumor.mutations)

        mutdict = Dict(mutations .=> 1:length(mutations))
        estrangedmuts = DataFrame(  mutation = mutations,
                                    isestranged = zeros(length(mutations)),
                                    isgreen = zeros(length(mutations))
                                    )
        for (i,m) in enumerate(mutations)
            trunc = truncs[i]
            truncp = filter( mp -> truncs[mutdict[mp]] == trunc, trunc)
            t = length(truncp)

            lastparentestranged = !(setdiff(trunc, truncp) in htypes)
            lastmutalive = trunc in htypes

            if t > 1
                estrangedmuts.isestranged[i] = lastmutalive/t
                estrangedmuts.isgreen[i] = lastmutalive/t
            elseif isone(length(trunc))
                estrangedmuts.isestranged[i] = 0.
                estrangedmuts.isgreen[i] = 0.
            else
                estrangedmuts.isestranged[i] = lastparentestranged && lastmutalive
                estrangedmuts.isgreen[i] = lastmutalive
            end

        end
        return estrangedmuts
    end

    ### orphaned tree-less

    """
        orphaned_green_treeless(htumor)

    Tree-less version of the `orphaned_green` algorithm on **incomplete (!) haplotype-tumors**. In an incomplete
    haplotype-tumor no haplotype has size `n` = 0 (those are to be removed beforehand).
    The function returns

        DataFrame( mutation :: Vector{Int}, isorphaned :: Vector{Float64}, isgreen :: Vector{Float64} )

    Summing over `isorphaned` gives the number of orphans and summing over `isgreen` the total number
    of offspring mutations.
    """
    function orphaned_green_treeless(htumor)
        # htumor = filter(h -> !iszero(h.n), tumor)

        mtypes, mutations = clones_by_mutations(htumor)

        truncs = [ intersect(Set.(mtype.mutations)...) for mtype in mtypes]

        mutdict = Dict(mutations .=> 1:length(mutations))
        orphanedmuts = DataFrame(  mutation = mutations,
                                    isorphaned = zeros(length(mutations)),
                                    isgreen = zeros(length(mutations))
                                    )
        for (i,m) in enumerate(mutations)
            trunc = truncs[i]
            truncp = filter( mp -> truncs[mutdict[mp]] == trunc, trunc)
            t = length(truncp)
            orphanedmuts.isorphaned[i] = 1 - 1/t
            orphanedmuts.isgreen[i] = truncp == trunc ? 1 - 1/t : 1
        end
        return orphanedmuts
    end

    """
        orphaned_red_treeless(htumor)

    Tree-less version of the `orphaned_red` algorithm on **incomplete (!) haplotype-tumors**. In an incomplete
    haplotype-tumor no haplotype has size `n` = 0 (those are to be removed beforehand).
    The function returns

        DataFrame( mutation :: Vector{Int}, isorphaned :: Vector{Float64}, isgreen :: Vector{Float64} )

    Summing over `isorphaned` gives the number of orphans and summing over `isgreen` the total number
    of offspring mutations.
    """
    function orphaned_red_treeless(htumor)
        # htumor = filter(h -> !iszero(h.n), tumor)

        mtypes, mutations = clones_by_mutations(htumor)

        truncs = [ intersect(Set.(mtype.mutations)...) for mtype in mtypes]

        mutdict = Dict(mutations .=> 1:length(mutations))
        orphanedmuts = DataFrame(  mutation = mutations,
                                    isorphaned = zeros(length(mutations)),
                                    isgreen = zeros(length(mutations))
                                    )
        for (i,m) in enumerate(mutations)
            trunc = truncs[i]
            twigs = setdiff( union(mtypes[i].mutations...), trunc )
            truncp = filter( mp -> truncs[mutdict[mp]] == trunc, trunc)
            t = length(truncp)
            orphanedmuts.isorphaned[i] = (t-1)/2
            orphanedmuts.isgreen[i] = length(twigs) + (t-1)/2
        end
        return orphanedmuts
    end

    """
        orphaned_red_expected(q, clonesizes, nbirth)

    Expected orphaned turnover considering full ancestry of a mutant clone. Given extinction probability `q` and observed ancestral clone sizes at birth `nbirth` an offspring mutation `m` is expected to be orphaned with respect to an ancestor of size `n` at birth of `m` with probability `q^n`.
    The function returns the probability of being orphaned for all offspring mutations as

        DataFrame( mutation :: Vector{Int}, isorphaned :: Vector{Float64})
    """
    function orphaned_red_expected(q, nbirth)
        extant = filter(m -> !isempty(nbirth[m]), 1:length(nbirth))
        orphans = DataFrame(mutation = Int[], isorphaned = Float64[])
        for m_green in extant
            parentsizes = nbirth[m_green]
            append!( orphans.mutation, fill(m_green, length(parentsizes)) )
            append!( orphans.isorphaned, q .^ parentsizes)
        end
        return orphans
    end

    """
        orphaned_green_expected(htumor, q, nbirth)

    Expected orphaned turnover based on offspring mutations on **complete (!) haplotype-tumors**,
    `complete` meaning that haplotypes with size `n` = 0 are listed. Given extinction probability `q`
    and observed parent clone sizes at birth `nbirth` a (green) offspring mutation is expected to be orphaned
    with probability `q^n`.
    The function returns the probability of being orphaned for all offspring mutations as

        DataFrame( mutation :: Vector{Int}, isorphaned :: Vector{Float64})
    """
    function orphaned_green_expected(htumor, q, nbirth)
        mtypes = clones_by_mutations(htumor)[1]
        freqs = sum.(getproperty.(mtypes,:n))

        orphans = DataFrame(mutation = Int[], isorphaned = Float64[])
        for (m,redclone) in enumerate(mtypes)
            m_green = last.(redclone.mutations[Not(1)])
            filter!( m-> !iszero(freqs[m]), m_green )
            isempty(m_green) && continue

            append!(orphans.mutation, m_green)
            append!(orphans.isorphaned, q .^ nbirth[m_green] )
        end
        return orphans
    end

    """
        estranged(htumor, q, nbirth)

    Estranged turnover on **complete (!) haplotype-tumors**,
    *complete* meaning that haplotypes with size *n* = 0 are listed. Given extinction probability `q`
    and observed parent clone sizes at birth `nbirth`, an offspring mutation is expected to be estranged
    with probability `q^n`.
    The function returns the probability of being orphaned for all offspring mutations as

        DataFrame( mutation :: Vector{Int}, isestranged :: Vector{Float64})
    """
    function estranged_expected(htumor, q, nbirth)
        freqs = mfreqs(htumor)
        daughters = filter(h -> h.n >= 1 && length(h.mutations)>1, htumor)

        mutations = last.(daughters.mutations)
        return DataFrame(  mutation = mutations, isestranged = q .^ nbirth[mutations])
    end
end  # module Turnover

# function estranged_treeless_v2(htumor)
#     clones, mutations = clones_by_mutations(htumor)
#     m_indices = Dict(mutations .=> 1:length(mutations))
#     htypes = filter(h-> length(h)>1, htumor.mutations)
#     out = DataFrame(  clone = htypes, isestranged = Vector{Bool}(undef, length(htypes)) )
#
#     for (i,h) in enumerate(htypes)
#         n = length(h)
#         isestranged = true
#         for m in h
#             isestranged || break
#             for potential_parent in clones[m_indices[m]].mutations
#                 if length(potential_parent) == n-1 && isone(length(setdiff(h, potential_parent)))
#                     isestranged = false
#                     break
#                 end
#             end
#         end
#         out.isestranged[i] = isestranged
#     end
#     return out
# end

max_stars_repo_name: aangaji/Turnover

max_issues_repo_name: aangaji/Turnover

max_forks_repo_name: aangaji/Turnover


text:
ProduceTestQuantity(n::Int) = n == -1 ? rand() : (n > 0 ? rand(n) : throw("Got invalid length value $n."))
ProduceTestQuantity(Tup::Tuple{Vararg{Int}}) = ((@assert all(x->x > 0, Tup) "Got invalid Tuple $Tup."); rand(Tup...))

GetOutLength(F::Function, In::Int=_GetArgLength(F)) = GetOutLength(F, ProduceTestQuantity(In))
"""
    GetOutLength(F::Function, input::Union{Number,AbstractVector{<:Number}})
Returns output dimensions of given `F`. If it outputs arrays of more than one dimension, a tuple is returned.
This can also be used to determine the approximate size of the input for mutating `F` which accept 2 arguments.

!!! note
    Discriminates between `Real` and `Vector{Real}` of length one, i.e.:
    `Real`↦`-1` and `x::AbstractVector{<:Real}`↦`length(x)`.
"""
function GetOutLength(F::Function, testinput::Union{Number,AbstractVector{<:Number}})
    if MaximalNumberOfArguments(F) == 1
        output = F(testinput)
        output isa Number ? -1 : (output isa AbstractVector ? length(output) : size(output))
    else
        _GetArgLengthInPlace(F)[1]
    end
end

function GetInOut(F::Function)
    if MaximalNumberOfArguments(F) == 1
        In = _GetArgLengthOutOfPlace(F);    Out = GetOutLength(F,In);    (In,Out)
    else
        # _GetArgLengthInPlace uses opposite order
        _GetArgLengthInPlace(F) |> reverse
    end
end

GetInOut(F::Function, input::Union{Number, AbstractVector{<:Number}}) = (In = (input isa Number ? -1 : length(input)); (In,GetOutLength(F,input)))

function _GetFirstDeriv(F::Function, InOut::Tuple{Int,Union{Int,Tuple}}=GetInOut(F); ADmode::Union{Val,Symbol}=Val(:Symbolic))
    if !(Out(InOut) isa Number)
        GetMatrixJac(ADmode,F,abs(In(InOut)),_MakeTuple(Out(InOut)))
    elseif In(InOut) == -1
        # Out is an Int but must be Tuple
        Out(InOut) == -1 ? GetDeriv(ADmode,F) : GetMatrixJac(ADmode,F,abs(In(InOut)),_MakeTuple(Out(InOut)))
    else
        Out(InOut) == -1 ? GetGrad(ADmode,F,abs(In(InOut))) : GetJac(ADmode,F,abs(In(InOut)))
    end
end
function _GetSecondDeriv(F::Function, dF::Function, InOut::Tuple{Int,Union{Int,Tuple}}=GetInOut(F); ADmode::Union{Val,Symbol}=Val(:Symbolic))
    if Out(InOut) isa Number && Out(InOut) == -1
        ## For scalar functions, Symbolics often produces false results.
        In(InOut) == -1 ? GetDeriv(ADmode, dF) : GetHess(Val(:ForwardDiff), F, abs(In(InOut)))
    else
        ## Already know info about dimensions, do not call GetDoubleJac.
        GetMatrixJac(ADmode,dF,abs(In(InOut)),(Out(InOut)...,abs(In(InOut))))
    end
end
function _GetSecondDeriv(F::Function, InOut::Tuple{Int,Union{Int,Tuple}}=GetInOut(F); ADmode::Union{Val,Symbol}=Val(:Symbolic))
    _GetSecondDeriv(F, _GetFirstDeriv(F,InOut;ADmode=ADmode), InOut; ADmode=ADmode)
end

"""
    DerivableFunction(F::Function; ADmode::Union{Val,Symbol}=Val(:Symbolic))
    DerivableFunction(F::Function, testinput::Union{Number,AbstractVector{<:Number}}; ADmode::Union{Val,Symbol}=Val(:Symbolic))
    DerivableFunction(F::Function, dF::Function; ADmode::Union{Val,Symbol}=Val(:Symbolic))
    DerivableFunction(F::Function, dF::Function, ddF::Function)
Stores derivatives of a given function (as well as input-output dimensions) for potentially faster computations when derivatives are known.
"""
mutable struct DerivableFunction <: Function
    F::Function
    dF::Function
    ddF::Function
    InOut::Tuple{Val,Val}
    function DerivableFunction(F::Function, testinput::Union{Number,AbstractVector{<:Number}}, InOut::Tuple{Int,Union{Int,Tuple}}=GetInOut(F,testinput); ADmode::Union{Val,Symbol}=Val(:Symbolic))
        DerivableFunction(F, _GetFirstDeriv(F,InOut;ADmode=ADmode), InOut; ADmode=ADmode)
    end
    function DerivableFunction(F::Function, InOut::Tuple{Int,Union{Int,Tuple}}=GetInOut(F); ADmode::Union{Val,Symbol}=Val(:Symbolic))
        DerivableFunction(F, _GetFirstDeriv(F,InOut;ADmode=ADmode), InOut; ADmode=ADmode)
    end
    function DerivableFunction(F::Function, dF::Function, InOut::Tuple{Int,Union{Int,Tuple}}=GetInOut(F); ADmode::Union{Val,Symbol}=Val(:Symbolic))
        DerivableFunction(F, dF, _GetSecondDeriv(F,dF,InOut;ADmode=ADmode), InOut)
    end
    function DerivableFunction(F::Function, dF::Function, ddF::Function, InOut::Tuple{Int,Union{Int,Tuple}}=GetInOut(F))
        DerivableFunction(F, dF, ddF, (Val(InOut[1]), Val(InOut[2])))
    end
    function DerivableFunction(F::Function, dF::Function, ddF::Function, InOut::Tuple{Val,Val})
        @assert In(InOut) ≥ -1 && (Out(InOut) isa Number ? (Out(InOut) ≥ -1) : all(x->x≥1, Out(InOut)))
        new(F, dF, ddF, InOut)
    end
    DerivableFunction(D::DerivableFunction, args...; kwargs...) = D
end
DFunction = DerivableFunction

(D::DFunction)(x) = EvalF(D, x)

InOut(D::DFunction) = (In(D), Out(D))
In(D::DFunction) = In(D.InOut);      In(InOut::Tuple) = _content(InOut[1])
Out(D::DFunction) = Out(D.InOut);    Out(InOut::Tuple) = _content(InOut[2])
_content(::Val{N}) where N = _content(N)
_content(N::Int) = N;   _content(tup::Tuple) = tup
_content(N) = throw("Unexpected content.")

EvalF(D::DFunction, x; kwargs...) = D.F(x; kwargs...)
EvalF(D::DFunction) = D.F
EvaldF(D::DFunction, x; kwargs...) = D.dF(x; kwargs...)
EvaldF(D::DFunction) = D.dF
EvalddF(D::DFunction, x; kwargs...) = D.ddF(x; kwargs...)
EvalddF(D::DFunction) = D.ddF

EvalF(F::Function, x; kwargs...) = F(x; kwargs...)
EvalF(F::Function) = F

# No inference is performed on the input function here, it is assumed to not output scalar values
EvaldF(F::Function, x; ADmode::Union{Symbol,Val}=Val(:ForwardDiff)) = GetMatrixJac(ADmode, F)(x)
EvaldF(F::Function; ADmode::Union{Symbol,Val}=Val(:ForwardDiff)) = GetMatrixJac(ADmode, F)
function EvalddF(F::Function, x; ADmode::Union{Symbol,Val}=Val(:ForwardDiff))
    m=GetArgLength(F);  f=_SizeTuple(F, m)
    GetMatrixJac(ADmode, GetMatrixJac(ADmode, F, m, f), m, (f..., m))(x)
end
function EvalddF(F::Function; ADmode::Union{Symbol,Val}=Val(:ForwardDiff))
    m=GetArgLength(F);  f=_SizeTuple(F, m)
    GetMatrixJac(ADmode, GetMatrixJac(ADmode, F, m, f), m, (f..., m))
end

_InputSpace(D::DFunction) = _SpaceWord(In(D))
_OutputSpace(D::DFunction) = _SpaceWord(Out(D))
_SpaceWord(n::Int) = n == -1 ? "R" : "R^$n"
_SpaceWord(Tup::Tuple) = "R^(" * prod([string(n)*"×" for n in Tup])[1:end-1] * ")"
Base.show(io::IO, D::DFunction) = println(io, "DerivableFunction : " * _InputSpace(D) * " ⟶  " * _OutputSpace(D))


text:
# ------------------------------------------------------------------
# Licensed under the ISC License. See LICENSE in the project root.
# ------------------------------------------------------------------

# ----------------
# IMPLEMENTATIONS
# ----------------
include("solutions/estimation.jl")
include("solutions/simulation.jl")


text:
module test_dataframes_df

using Test
using DataFrames

# result = Repo.query([SELECT * FROM Temp])
result = [(id=1, title="Hello"),
          (id=2, title="World"),
         ]

df = DataFrame(result, copycols=false)
@test df.id == [1, 2]
@test df.title == ["Hello", "World"]
@test size(df) == (2, 2)
@test result == NamedTuple.(eachrow(df))

end # module test_dataframes_df

using Jive
@If VERSION >= v"1.7" module test_dataframes_df_destructing

using Test
using DataFrames

result = [(id=1, title="Hello"),
          (id=2, title="World"),
         ]

df = DataFrame(result, copycols=false)

row = eachrow(df)[1]
(; id, title) = row
@test id == 1
@test title == "Hello"

titles = []
for (; title) in eachrow(df)
    push!(titles, title)
end
@test titles == ["Hello", "World"]

end # @If VERSION >= v"1.7" module test_dataframes_df_destructing


text:
using Test
using ClosedIntervals

A = ClosedInterval(1,5)
B = ClosedInterval(7,9)
C = ClosedInterval(9,1)

@test A+B == C
@test isempty(A*B)
@test A == ClosedInterval(right(A),left(A))
@test length(A) == 4

A = ClosedInterval(4.)
B = ClosedInterval(1.)
C = ClosedInterval(4.,1.)

@test A+B == C
@test length(A)==0
@test isempty(A*B)

A = ClosedInterval(1,6)
B = ClosedInterval(8,3)
C = ClosedInterval(3,6)

@test A*B == C
@test A ∧ B == C
@test A+B == ClosedInterval(1,8)
@test A ∨ B == ClosedInterval(1,8)

J = ClosedInterval(0,0)
K = EmptyInterval(Int)
@test J != K

KK = EmptyInterval(Float64)
@test K==KK

@test 5..3 == ClosedInterval(3,5)
@test 2±1  == ClosedInterval(1,3)


text:
immutable QuasiHyperbolicIntercept{T} <: Model
    β₀::T
    a::T
    β::T
    δ::T
end

function Base.show(io::IO, ::Type{QuasiHyperbolicIntercept})
    return print(io, "Intercept Quasi-Hyperbolic")
end

function Base.show(io::IO, m::QuasiHyperbolicIntercept)
    @printf(
        io,
        "Intercept Quasi-Hyperbolic (β₀ = %.3f, a = %.3f, β = %.3f, δ = %.3f)",
        m.β₀,
        m.a,
        m.β,
        m.δ,
    )
end

function initial_parameters(::Type{QuasiHyperbolicIntercept})
    return Float64[0.0, log(1.0), logit(0.5), logit(0.5)]
end

canonical_parameters(m::QuasiHyperbolicIntercept) = (m.β₀, m.a, m.β, m.δ)

function QuasiHyperbolicIntercept(parameters::Vector)
    β₀ = parameters[1]
    a = exp(parameters[2])
    β = invlogit(parameters[3])
    δ = invlogit(parameters[4])
    return QuasiHyperbolicIntercept(β₀, a, β, δ)
end

@inline function predict(
    m::QuasiHyperbolicIntercept,
    x2::Real,
    t2::Real,
    x1::Real,
    t1::Real,
)
    β₀, a, β, δ = m.β₀, m.a, m.β, m.δ
    u1 = x1 * ifelse(t1 == 0.0, 1.0, β) * δ^t1
    u2 = x2 * ifelse(t2 == 0.0, 1.0, β) * δ^t2
    z = β₀ + a * (u2 - u1)
    p = invlogit(z)
    return p
end

@inline function gradient_component(
    m::QuasiHyperbolicIntercept,
    x2::Real,
    t2::Real,
    x1::Real,
    t1::Real,
)
    β₀, a, β, δ = m.β₀, m.a, m.β, m.δ
    βt2 = ifelse(t2 == 0.0, 1.0, β)
    βt1 = ifelse(t1 == 0.0, 1.0, β)
    Δβt2 = ifelse(t2 == 0.0, 0.0, 1.0)
    Δβt1 = ifelse(t1 == 0.0, 0.0, 1.0)
    δt2 = δ^t2
    δt1 = δ^t1
    gr = (
        1.0,
        (x2 * βt2 * δt2 - x1 * βt1 * δt1) * a,
        a * (x2 * Δβt2 * δt2 - x1 * Δβt1 * δt1) * β * (1 - β),
        a * ((x2 * βt2 * t2 * δt2 - x1 * βt1 * t1 * δt1) / δ) * δ * (1 - δ),
    )
    return gr
end


text:
@test simplify(:x/:x,Div)==1

E=@equ E=sqrt(p^2*c^2+m^2*c^4)
vars=@equs p=sqrt(2)*1e6/c m=0.5e6/c^2
r=E&vars[2]
@test r.rhs==Equations.Sqrt(:p*:p*:c*:c+2.5e11)