text:
abstract type NUTSState{F<:VariateForm} <: HMCSamplerState{F} end

mutable struct UnvNUTSState <: NUTSState{Univariate}
  pstateplus::ParameterState{Continuous, Univariate}
  pstateminus::ParameterState{Continuous, Univariate}
  pstateprime::ParameterState{Continuous, Univariate}
  pstatedprime::ParameterState{Continuous, Univariate}
  tune::MCTunerState
  ratio::Real
  momentum::Real
  momentumplus::Real
  momentumminus::Real
  momentumprime::Real
  oldhamiltonian::Real
  newhamiltonian::Real
  u::Real
  v::Integer
  j::Integer
  n::Integer
  nprime::Integer
  ndprime::Integer
  s::Bool
  sprime::Bool
  sdprime::Bool
  update::Bool
  count::Integer
  diagnosticindices::Dict{Symbol, Integer}

  function UnvNUTSState(
    pstateplus::ParameterState{Continuous, Univariate},
    pstateminus::ParameterState{Continuous, Univariate},
    pstateprime::ParameterState{Continuous, Univariate},
    pstatedprime::ParameterState{Continuous, Univariate},
    tune::MCTunerState,
    ratio::Real,
    momentum::Real,
    momentumplus::Real,
    momentumminus::Real,
    momentumprime::Real,
    oldhamiltonian::Real,
    newhamiltonian::Real,
    u::Real,
    v::Integer,
    j::Integer,
    n::Integer,
    nprime::Integer,
    ndprime::Integer,
    s::Bool,
    sprime::Bool,
    sdprime::Bool,
    update::Bool,
    count::Integer,
    diagnosticindices::Dict{Symbol, Integer}
  )
    if !isnan(ratio)
      @assert ratio > 0 "Acceptance ratio should be positive"
    end
    @assert count >= 0 "Number of iterations (count) should be non-negative"
    new(
      pstateplus,
      pstateminus,
      pstateprime,
      pstatedprime,
      tune,
      ratio,
      momentum,
      momentumplus,
      momentumminus,
      momentumprime,
      oldhamiltonian,
      newhamiltonian,
      u,
      v,
      j,
      n,
      nprime,
      ndprime,
      s,
      sprime,
      sdprime,
      update,
      count,
      diagnosticindices
    )
  end
end

UnvNUTSState(pstate::ParameterState{Continuous, Univariate}, tune::MCTunerState=BasicMCTune()) =
  UnvNUTSState(
    pstate,
    pstate,
    pstate,
    pstate,
    tune,
    NaN,
    NaN,
    NaN,
    NaN,
    NaN,
    NaN,
    NaN,
    NaN,
    0,
    0,
    0,
    0,
    0,
    true,
    true,
    true,
    true,
    0,
    Dict{Symbol, Integer}()
  )

mutable struct MuvNUTSState <: NUTSState{Multivariate}
  pstateplus::ParameterState{Continuous, Multivariate}
  pstateminus::ParameterState{Continuous, Multivariate}
  pstateprime::ParameterState{Continuous, Multivariate}
  pstatedprime::ParameterState{Continuous, Multivariate}
  tune::MCTunerState
  ratio::Real
  momentum::RealVector
  momentumplus::RealVector
  momentumminus::RealVector
  momentumprime::RealVector
  oldhamiltonian::Real
  newhamiltonian::Real
  u::Real
  v::Integer
  j::Integer
  n::Integer
  nprime::Integer
  ndprime::Integer
  s::Bool
  sprime::Bool
  sdprime::Bool
  update::Bool
  count::Integer
  diagnosticindices::Dict{Symbol, Integer}

  function MuvNUTSState(
    pstateplus::ParameterState{Continuous, Multivariate},
    pstateminus::ParameterState{Continuous, Multivariate},
    pstateprime::ParameterState{Continuous, Multivariate},
    pstatedprime::ParameterState{Continuous, Multivariate},
    tune::MCTunerState,
    ratio::Real,
    momentum::RealVector,
    momentumplus::RealVector,
    momentumminus::RealVector,
    momentumprime::RealVector,
    oldhamiltonian::Real,
    newhamiltonian::Real,
    u::Real,
    v::Integer,
    j::Integer,
    n::Integer,
    nprime::Integer,
    ndprime::Integer,
    s::Bool,
    sprime::Bool,
    sdprime::Bool,
    update::Bool,
    count::Integer,
    diagnosticindices::Dict{Symbol, Integer}
  )
    if !isnan(ratio)
      @assert ratio > 0 "Acceptance ratio should be positive"
    end
    @assert count >= 0 "Number of iterations (count) should be non-negative"
    new(
      pstateplus,
      pstateminus,
      pstateprime,
      pstatedprime,
      tune,
      ratio,
      momentum,
      momentumplus,
      momentumminus,
      momentumprime,
      oldhamiltonian,
      newhamiltonian,
      u,
      v,
      j,
      n,
      nprime,
      ndprime,
      s,
      sprime,
      sdprime,
      update,
      count,
      diagnosticindices
    )
  end
end

MuvNUTSState(pstate::ParameterState{Continuous, Multivariate}, tune::MCTunerState=BasicMCTune()) =
  MuvNUTSState(
    pstate,
    pstate,
    pstate,
    pstate,
    tune,
    NaN,
    Array{eltype(pstate)}(pstate.size),
    Array{eltype(pstate)}(pstate.size),
    Array{eltype(pstate)}(pstate.size),
    Array{eltype(pstate)}(pstate.size),
    NaN,
    NaN,
    NaN,
    0,
    0,
    0,
    0,
    0,
    true,
    true,
    true,
    true,
    0,
    Dict{Symbol, Integer}()
  )

mutable struct NUTS <: HMCSampler
  leapstep::Real
  maxδ::Integer
  maxndoublings::Integer

  function NUTS(leapstep::Real, maxδ::Integer, maxndoublings::Integer)
    @assert leapstep > 0 "Leapfrog step is not positive"
    @assert maxδ > 0 "maxδ is not positive"
    @assert maxndoublings > 0 "Maximum number of doublings is not positive"
    new(leapstep, maxδ, maxndoublings)
  end
end

NUTS(leapstep::Real=0.1; maxδ::Integer=1000, maxndoublings::Integer=5) = NUTS(leapstep, maxδ, maxndoublings)

function initialize!(
  pstate::ParameterState{Continuous, F},
  parameter::Parameter{Continuous, F},
  sampler::NUTS,
  outopts::Dict
) where F<:VariateForm
  parameter.uptogradlogtarget!(pstate)
  @assert isfinite(pstate.logtarget) "Log-target not finite: initial value out of support"
  @assert all(isfinite.(pstate.gradlogtarget)) "Gradient of log-target not finite: initial values out of support"

  if !isempty(outopts[:diagnostics])
    pstate.diagnostickeys = copy(outopts[:diagnostics])
    pstate.diagnosticvalues = Array{Any}(length(pstate.diagnostickeys))
  end
end

tuner_state(parameter::Parameter, sampler::NUTS, tuner::DualAveragingMCTuner) =
  DualAveragingMCTune(
  step=sampler.leapstep,
  λ=NaN,
  εbar=tuner.ε0bar,
  hbar=tuner.h0bar,
  accepted=0,
  proposed=0,
  totproposed=tuner.period
)

function sampler_state(
  parameter::Parameter{Continuous, Univariate},
  sampler::NUTS,
  tuner::VanillaMCTuner,
  pstate::ParameterState{Continuous, Univariate},
  vstate::VariableStateVector,
  diagnostickeys::Vector{Symbol}
)
  sstate = UnvNUTSState(
    generate_empty(pstate, parameter.diffmethods, parameter.diffopts), tuner_state(parameter, sampler, tuner)
  )
  set_diagnosticindices!(sstate, [:accept, :ndoublings], diagnostickeys)
  sstate
end

function sampler_state(
  parameter::Parameter{Continuous, Multivariate},
  sampler::NUTS,
  tuner::VanillaMCTuner,
  pstate::ParameterState{Continuous, Multivariate},
  vstate::VariableStateVector,
  diagnostickeys::Vector{Symbol}
)
  sstate = MuvNUTSState(
    generate_empty(pstate, parameter.diffmethods, parameter.diffopts), tuner_state(parameter, sampler, tuner)
  )
  set_diagnosticindices!(sstate, [:accept, :ndoublings], diagnostickeys)
  sstate
end

function sampler_state(
  parameter::Parameter{Continuous, Univariate},
  sampler::NUTS,
  tuner::DualAveragingMCTuner,
  pstate::ParameterState{Continuous, Univariate},
  vstate::VariableStateVector,
  diagnostickeys::Vector{Symbol}
)
  sstate = UnvNUTSState(
    generate_empty(pstate, parameter.diffmethods, parameter.diffopts), tuner_state(parameter, sampler, tuner)
  )

  sstate.tune.step = initialize_step!(
    sstate.pstateplus, pstate, randn(), sstate.tune.step, parameter.gradlogtarget!, typeof(tuner)
  )
  sstate.tune.μ = log(10*sstate.tune.step)
  set_diagnosticindices!(sstate, [:accept, :ndoublings, :a, :na], diagnostickeys)

  sstate
end

function sampler_state(
  parameter::Parameter{Continuous, Multivariate},
  sampler::NUTS,
  tuner::DualAveragingMCTuner,
  pstate::ParameterState{Continuous, Multivariate},
  vstate::VariableStateVector,
  diagnostickeys::Vector{Symbol}
)
  sstate = MuvNUTSState(
    generate_empty(pstate, parameter.diffmethods, parameter.diffopts), tuner_state(parameter, sampler, tuner)
  )

  sstate.tune.step = initialize_step!(
    sstate.pstateplus,
    sstate.momentum,
    pstate,
    randn(pstate.size),
    sstate.tune.step,
    parameter.gradlogtarget!,
    typeof(tuner)
  )
  sstate.tune.μ = log(10*sstate.tune.step)
  set_diagnosticindices!(sstate, [:accept, :ndoublings, :a, :na], diagnostickeys)

  sstate
end

function reset!(tune::DualAveragingMCTune, sampler::NUTS, tuner::DualAveragingMCTuner)
  tune.step = 1
  tune.εbar = tuner.ε0bar
  tune.hbar = tuner.h0bar
  (tune.accepted, tune.proposed, tune.totproposed, tune.rate) = (0, 0, tuner.period, NaN)
end

function reset!(
  sstate::MuvHMCState,
  pstate::ParameterState{Continuous, Univariate},
  parameter::Parameter{Continuous, Univariate},
  sampler::NUTS,
  tuner::DualAveragingMCTuner
)
  reset!(sstate.tune, sampler, tuner)
  sstate.tune.step = initialize_step!(
    sstate.pstateplus, pstate, randn(), sstate.tune.step, parameter.gradlogtarget!, typeof(tuner)
  )
  sstate.tune.μ = log(10*sstate.tune.step)
  sstate.count = 0
end

function reset!(
  sstate::MuvHMCState,
  pstate::ParameterState{Continuous, Multivariate},
  parameter::Parameter{Continuous, Multivariate},
  sampler::NUTS,
  tuner::DualAveragingMCTuner
)
  reset!(sstate.tune, sampler, tuner)
  sstate.tune.step = initialize_step!(
    sstate.pstateplus,
    sstate.momentum,
    pstate,
    randn(pstate.size),
    sstate.tune.step,
    parameter.gradlogtarget!,
    typeof(tuner)
  )
  sstate.tune.μ = log(10*sstate.tune.step)
  sstate.count = 0
end

uturn(positionplus::Real, positionminus::Real, momentumplus::Real, momentumminus::Real) =
  (positionplus-positionminus)*momentumplus < 0. || (positionplus-positionminus)*momentumminus < 0.

uturn(positionplus::RealVector, positionminus::RealVector, momentumplus::RealVector, momentumminus::RealVector) =
  dot(positionplus-positionminus, momentumplus) < 0. || dot(positionplus-positionminus, momentumminus) < 0.

function build_tree!(
  sstate::UnvNUTSState,
  pstate::ParameterState{Continuous, Univariate},
  momentum::Real,
  u::Real,
  v::Real,
  j::Integer,
  step::Real,
  parameter::Parameter{Continuous, Univariate},
  maxδ::Integer,
  ::Type{NUTS},
  ::Type{VanillaMCTuner},
  outopts::Dict
)
  if j == 0
    local hamiltonianprime::Real

    sstate.momentumprime = leapfrog!(sstate.pstateprime, pstate, momentum, v*sstate.tune.step, parameter.gradlogtarget!)

    parameter.logtarget!(sstate.pstateprime)
    hamiltonianprime = hamiltonian(sstate.pstateprime.logtarget, sstate.momentumprime)

    sstate.nprime = Int(u <= hamiltonianprime)
    sstate.sprime = u < maxδ+hamiltonianprime

    return sstate.pstateprime,
      sstate.momentumprime,
      sstate.pstateprime,
      sstate.momentumprime,
      sstate.pstateprime,
      sstate.nprime,
      sstate.sprime
  else
    (
      sstate.pstateminus,
      sstate.momentumminus,
      sstate.pstateplus,
      sstate.momentumplus,
      sstate.pstateprime,
      sstate.nprime,
      sstate.sprime
    ) =
      build_tree!(sstate, pstate, momentum, u, v, j-1, step, parameter, maxδ, NUTS, VanillaMCTuner, outopts)

    if sstate.sprime
      if v == -1
        (sstate.pstateminus, sstate.momentumminus, _, _, sstate.pstatedprime, sstate.ndprime, sstate.sdprime) =
          build_tree!(
            sstate,
            sstate.pstateminus,
            sstate.momentumminus,
            u,
            v,
            j-1,
            step,
            parameter,
            maxδ,
            NUTS,
            VanillaMCTuner,
            outopts
          )
      else
        (_, _, sstate.pstateplus, sstate.momentumplus, sstate.pstatedprime, sstate.ndprime, sstate.sdprime) =
          build_tree!(
            sstate,
            sstate.pstateplus,
            sstate.momentumplus,
            u,
            v,
            j-1,
            step,
            parameter,
            maxδ,
            NUTS,
            VanillaMCTuner,
            outopts
          )
      end

      if (rand() <= sstate.ndprime/(sstate.ndprime+sstate.nprime))
        sstate.pstateprime.value = sstate.pstatedprime.value

        sstate.pstateprime.gradlogtarget = sstate.pstatedprime.gradlogtarget
        if in(:gradloglikelihood, outopts[:monitor]) && parameter.gradloglikelihood! != nothing
          sstate.pstateprime.gradloglikelihood = sstate.pstatedprime.gradloglikelihood
        end
        if in(:gradlogprior, outopts[:monitor]) && parameter.gradlogprior! != nothing
          sstate.pstateprime.gradlogprior = sstate.pstatedprime.gradlogprior
        end

        sstate.pstateprime.logtarget = sstate.pstatedprime.logtarget
        if in(:loglikelihood, outopts[:monitor]) && parameter.loglikelihood! != nothing
          sstate.pstateprime.loglikelihood = sstate.pstatedprime.loglikelihood
        end
        if in(:logprior, outopts[:monitor]) && parameter.logprior! != nothing
          sstate.pstateprime.logprior = sstate.pstatedprime.logprior
        end
      end

      sstate.nprime += sstate.ndprime

      sstate.sprime =
        sstate.sdprime &&
        !uturn(sstate.pstateplus.value, sstate.pstateminus.value, sstate.momentumplus, sstate.momentumminus)
    end

    return sstate.pstateminus,
      sstate.momentumminus,
      sstate.pstateplus,
      sstate.momentumplus,
      sstate.pstateprime,
      sstate.nprime,
      sstate.sprime
  end
end

function build_tree!(
  sstate::MuvNUTSState,
  pstate::ParameterState{Continuous, Multivariate},
  momentum::RealVector,
  u::Real,
  v::Real,
  j::Integer,
  step::Real,
  parameter::Parameter{Continuous, Multivariate},
  maxδ::Integer,
  ::Type{NUTS},
  ::Type{VanillaMCTuner},
  outopts::Dict
)
  if j == 0
    local hamiltonianprime::Real

    leapfrog!(sstate.pstateprime, sstate.momentumprime, pstate, momentum, v*sstate.tune.step, parameter.gradlogtarget!)

    parameter.logtarget!(sstate.pstateprime)
    hamiltonianprime = hamiltonian(sstate.pstateprime.logtarget, sstate.momentumprime)

    sstate.nprime = Int(u <= hamiltonianprime)
    sstate.sprime = u < maxδ+hamiltonianprime

    return sstate.pstateprime,
      sstate.momentumprime,
      sstate.pstateprime,
      sstate.momentumprime,
      sstate.pstateprime,
      sstate.nprime,
      sstate.sprime
  else
    (
      sstate.pstateminus,
      sstate.momentumminus,
      sstate.pstateplus,
      sstate.momentumplus,
      sstate.pstateprime,
      sstate.nprime,
      sstate.sprime
    ) =
      build_tree!(sstate, pstate, momentum, u, v, j-1, step, parameter, maxδ, NUTS, VanillaMCTuner, outopts)

    if sstate.sprime
      if v == -1
        (sstate.pstateminus, sstate.momentumminus, _, _, sstate.pstatedprime, sstate.ndprime, sstate.sdprime) =
          build_tree!(
            sstate,
            sstate.pstateminus,
            sstate.momentumminus,
            u,
            v,
            j-1,
            step,
            parameter,
            maxδ,
            NUTS,
            VanillaMCTuner,
            outopts
          )
      else
        (_, _, sstate.pstateplus, sstate.momentumplus, sstate.pstatedprime, sstate.ndprime, sstate.sdprime) =
          build_tree!(
            sstate,
            sstate.pstateplus,
            sstate.momentumplus,
            u,
            v,
            j-1,
            step,
            parameter,
            maxδ,
            NUTS,
            VanillaMCTuner,
            outopts
          )
      end

      if (rand() <= sstate.ndprime/(sstate.ndprime+sstate.nprime))
        sstate.pstateprime.value = copy(sstate.pstatedprime.value)

        sstate.pstateprime.gradlogtarget = copy(sstate.pstatedprime.gradlogtarget)
        if in(:gradloglikelihood, outopts[:monitor]) && parameter.gradloglikelihood! != nothing
          sstate.pstateprime.gradloglikelihood = copy(sstate.pstatedprime.gradloglikelihood)
        end
        if in(:gradlogprior, outopts[:monitor]) && parameter.gradlogprior! != nothing
          sstate.pstateprime.gradlogprior = copy(sstate.pstatedprime.gradlogprior)
        end

        sstate.pstateprime.logtarget = sstate.pstatedprime.logtarget
        if in(:loglikelihood, outopts[:monitor]) && parameter.loglikelihood! != nothing
          sstate.pstateprime.loglikelihood = sstate.pstatedprime.loglikelihood
        end
        if in(:logprior, outopts[:monitor]) && parameter.logprior! != nothing
          sstate.pstateprime.logprior = sstate.pstatedprime.logprior
        end
      end

      sstate.nprime += sstate.ndprime

      sstate.sprime =
        sstate.sdprime &&
        !uturn(sstate.pstateplus.value, sstate.pstateminus.value, sstate.momentumplus, sstate.momentumminus)
    end

    return sstate.pstateminus,
      sstate.momentumminus,
      sstate.pstateplus,
      sstate.momentumplus,
      sstate.pstateprime,
      sstate.nprime,
      sstate.sprime
  end
end

function build_tree!(
  sstate::UnvNUTSState,
  pstate::ParameterState{Continuous, Univariate},
  momentum::Real,
  oldhamiltonian::Real,
  u::Real,
  v::Real,
  j::Integer,
  step::Real,
  parameter::Parameter{Continuous, Univariate},
  maxδ::Integer,
  ::Type{NUTS},
  ::Type{DualAveragingMCTuner},
  outopts::Dict
)
  local aprime::Real
  local adprime::Real
  local naprime::Integer
  local nadprime::Integer

  if j == 0
    local hamiltonianprime::Real

    sstate.momentumprime = leapfrog!(sstate.pstateprime, pstate, momentum, v*sstate.tune.step, parameter.gradlogtarget!)

    parameter.logtarget!(sstate.pstateprime)
    hamiltonianprime = hamiltonian(sstate.pstateprime.logtarget, sstate.momentumprime)

    sstate.nprime = Int(u <= hamiltonianprime)
    sstate.sprime = u < maxδ+hamiltonianprime

    return sstate.pstateprime,
      sstate.momentumprime,
      sstate.pstateprime,
      sstate.momentumprime,
      sstate.pstateprime,
      sstate.nprime,
      sstate.sprime,
      min(1, exp(hamiltonianprime-oldhamiltonian)),
      1
  else
    (
      sstate.pstateminus,
      sstate.momentumminus,
      sstate.pstateplus,
      sstate.momentumplus,
      sstate.pstateprime,
      sstate.nprime,
      sstate.sprime,
      aprime,
      naprime
    ) =
      build_tree!(
        sstate, pstate, momentum, oldhamiltonian, u, v, j-1, step, parameter, maxδ, NUTS, DualAveragingMCTuner, outopts
      )

    if sstate.sprime
      if v == -1
        (
          sstate.pstateminus,
          sstate.momentumminus,
          _,
          _,
          sstate.pstatedprime,
          sstate.ndprime,
          sstate.sdprime,
          adprime,
          nadprime
        ) =
          build_tree!(
            sstate,
            sstate.pstateminus,
            sstate.momentumminus,
            oldhamiltonian,
            u,
            v,
            j-1,
            step,
            parameter,
            maxδ,
            NUTS,
            DualAveragingMCTuner,
            outopts
          )
      else
        (
          _,
          _,
          sstate.pstateplus,
          sstate.momentumplus,
          sstate.pstatedprime,
          sstate.ndprime,
          sstate.sdprime,
          adprime,
          nadprime
        ) =
          build_tree!(
            sstate,
            sstate.pstateplus,
            sstate.momentumplus,
            oldhamiltonian,
            u,
            v,
            j-1,
            step,
            parameter,
            maxδ,
            NUTS,
            DualAveragingMCTuner,
            outopts
          )
      end

      if (rand() <= sstate.ndprime/(sstate.ndprime+sstate.nprime))
        sstate.pstateprime.value = sstate.pstatedprime.value

        sstate.pstateprime.gradlogtarget = sstate.pstatedprime.gradlogtarget
        if in(:gradloglikelihood, outopts[:monitor]) && parameter.gradloglikelihood! != nothing
          sstate.pstateprime.gradloglikelihood = sstate.pstatedprime.gradloglikelihood
        end
        if in(:gradlogprior, outopts[:monitor]) && parameter.gradlogprior! != nothing
          sstate.pstateprime.gradlogprior = sstate.pstatedprime.gradlogprior
        end

        sstate.pstateprime.logtarget = sstate.pstatedprime.logtarget
        if in(:loglikelihood, outopts[:monitor]) && parameter.loglikelihood! != nothing
          sstate.pstateprime.loglikelihood = sstate.pstatedprime.loglikelihood
        end
        if in(:logprior, outopts[:monitor]) && parameter.logprior! != nothing
          sstate.pstateprime.logprior = sstate.pstatedprime.logprior
        end
      end

      sstate.nprime += sstate.ndprime

      sstate.sprime =
        sstate.sdprime &&
        !uturn(sstate.pstateplus.value, sstate.pstateminus.value, sstate.momentumplus, sstate.momentumminus)

      aprime += adprime
      naprime += nadprime
    end

    return sstate.pstateminus,
      sstate.momentumminus,
      sstate.pstateplus,
      sstate.momentumplus,
      sstate.pstateprime,
      sstate.nprime,
      sstate.sprime,
      aprime,
      naprime
  end
end

function build_tree!(
  sstate::MuvNUTSState,
  pstate::ParameterState{Continuous, Multivariate},
  momentum::RealVector,
  oldhamiltonian::Real,
  u::Real,
  v::Real,
  j::Integer,
  step::Real,
  parameter::Parameter{Continuous, Multivariate},
  maxδ::Integer,
  ::Type{NUTS},
  ::Type{DualAveragingMCTuner},
  outopts::Dict
)
  local aprime::Real
  local adprime::Real
  local naprime::Integer
  local nadprime::Integer

  if j == 0
    local hamiltonianprime::Real

    leapfrog!(sstate.pstateprime, sstate.momentumprime, pstate, momentum, v*sstate.tune.step, parameter.gradlogtarget!)

    parameter.logtarget!(sstate.pstateprime)
    hamiltonianprime = hamiltonian(sstate.pstateprime.logtarget, sstate.momentumprime)

    sstate.nprime = Int(u <= hamiltonianprime)
    sstate.sprime = u < maxδ+hamiltonianprime

    return sstate.pstateprime,
      sstate.momentumprime,
      sstate.pstateprime,
      sstate.momentumprime,
      sstate.pstateprime,
      sstate.nprime,
      sstate.sprime,
      min(1, exp(hamiltonianprime-oldhamiltonian)),
      1
  else
    (
      sstate.pstateminus,
      sstate.momentumminus,
      sstate.pstateplus,
      sstate.momentumplus,
      sstate.pstateprime,
      sstate.nprime,
      sstate.sprime,
      aprime,
      naprime
    ) =
      build_tree!(
        sstate, pstate, momentum, oldhamiltonian, u, v, j-1, step, parameter, maxδ, NUTS, DualAveragingMCTuner, outopts
      )

    if sstate.sprime
      if v == -1
        (
          sstate.pstateminus,
          sstate.momentumminus,
          _,
          _,
          sstate.pstatedprime,
          sstate.ndprime,
          sstate.sdprime,
          adprime,
          nadprime
        ) =
          build_tree!(
            sstate,
            sstate.pstateminus,
            sstate.momentumminus,
            oldhamiltonian,
            u,
            v,
            j-1,
            step,
            parameter,
            maxδ,
            NUTS,
            DualAveragingMCTuner,
            outopts
          )
      else
        (
          _,
          _,
          sstate.pstateplus,
          sstate.momentumplus,
          sstate.pstatedprime,
          sstate.ndprime,
          sstate.sdprime,
          adprime,
          nadprime
        ) =
          build_tree!(
            sstate,
            sstate.pstateplus,
            sstate.momentumplus,
            oldhamiltonian,
            u,
            v,
            j-1,
            step,
            parameter,
            maxδ,
            NUTS,
            DualAveragingMCTuner,
            outopts
          )
      end

      if (rand() <= sstate.ndprime/(sstate.ndprime+sstate.nprime))
        sstate.pstateprime.value = copy(sstate.pstatedprime.value)

        sstate.pstateprime.gradlogtarget = copy(sstate.pstatedprime.gradlogtarget)
        if in(:gradloglikelihood, outopts[:monitor]) && parameter.gradloglikelihood! != nothing
          sstate.pstateprime.gradloglikelihood = copy(sstate.pstatedprime.gradloglikelihood)
        end
        if in(:gradlogprior, outopts[:monitor]) && parameter.gradlogprior! != nothing
          sstate.pstateprime.gradlogprior = copy(sstate.pstatedprime.gradlogprior)
        end

        sstate.pstateprime.logtarget = sstate.pstatedprime.logtarget
        if in(:loglikelihood, outopts[:monitor]) && parameter.loglikelihood! != nothing
          sstate.pstateprime.loglikelihood = sstate.pstatedprime.loglikelihood
        end
        if in(:logprior, outopts[:monitor]) && parameter.logprior! != nothing
          sstate.pstateprime.logprior = sstate.pstatedprime.logprior
        end
      end

      sstate.nprime += sstate.ndprime

      sstate.sprime =
        sstate.sdprime &&
        !uturn(sstate.pstateplus.value, sstate.pstateminus.value, sstate.momentumplus, sstate.momentumminus)

      aprime += adprime
      naprime += nadprime
    end

    return sstate.pstateminus,
      sstate.momentumminus,
      sstate.pstateplus,
      sstate.momentumplus,
      sstate.pstateprime,
      sstate.nprime,
      sstate.sprime,
      aprime,
      naprime
  end
end

show(io::IO, sampler::NUTS) =
  print(
    io,
    "NUTS sampler: leap step = ",
    sampler.leapstep,
    ", doubling threshold = ",
    sampler.maxδ,
    ", maximum number of doublings = ",
    sampler.maxndoublings
  )


text:
meanVect = [27.155405405405407, 26.163835263835264]
covMat = [16.1253895583728 13.046961200891614; 13.046961200891614 12.367253570765161]

ext: jl

max_stars_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_issues_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_forks_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook


text:
using StructuredMatrices: nlogdet
using ReverseDiffExpressionsBase: tadd

"""
σ⁻²   : multiply with δ² / -2 to get kernel-component of density
ld    : multiply with -logdet_coef to get logdet component of density
∂k∂i  : multiply with δ²σ⁻² to get ∂kernel-component/∂input
∂ld∂i : multiply with -logdet_coef to get ∂logdet-component/∂input
"""
struct Precision{T,P,D}
    σ⁻²::T
    ld::D
    ∂k∂i::P
end
struct PrecisionArray{S,T,N,X,L,V<:AbstractFixedSizeArray{S,T,N,X,L},P,D<:AbstractFixedSizeArray{S,T,N,X,L}} <: AbstractFixedSizeArray{S,T,N,X,L}
    σ⁻²::V
    ld::D
    ∂k∂i::P
end
@inline Base.getindex(A::PrecisionArray, i...) = A.σ⁻²[i...]
const AbstractPrecision{T,P,D} = Union{Precision{T,P,D},PrecisionArray{<:Any,<:Any,<:Any,<:Any,<:Any,T,P,D}}

using ReverseDiffExpressionsBase: One, Zero

@inline logdet_coef(Y, args...) = Float64(size(Y, 1))
@inline function canonicalize_Σ(σ::Real)
    σ⁻¹ = Base.FastMath.inv_fast( VectorizationBase.extract_data(σ) )
    σ⁻² = Base.FastMath.abs2_fast( σ⁻¹ )
    Precision( σ⁻², log(σ), nothing )
end
@inline function ∂canonicalize_Σ(σ::Real)
    σ⁻¹ = Base.FastMath.inv_fast(VectorizationBase.extract_data(σ))
    σ⁻² = Base.FastMath.abs2_fast( σ⁻¹ )
    # nσ⁻¹ = Base.FastMath.sub_fast( σ⁻¹ )
    Precision( σ⁻², log(σ), σ⁻¹ )
end
@inline precision(v) = SIMDPirates.vabs2(SIMDPirates.vinv(v))
@inline function canonicalize_Σ(σ::AbstractFixedSizeArray{S,T,N,X,L}) where {S,T,N,X,L}
    σ⁻² = LazyMap( precision, σ )
    logσ = LazyMap( SLEEFPirates.log, σ )
    PrecisionArray{S,T,N,X,L,typeof(σ⁻²),Nothing,typeof(logσ)}( σ⁻², logσ, nothing )
end
@inline function ∂canonicalize_Σ(σ::AbstractFixedSizeArray{S,T,N,X,L}) where {S,T,N,X,L}
    # Trust compiler to elliminate redundant
    σ⁻¹ = LazyMap( vinv, σ )
    σ⁻² = LazyMap( precision, σ )
    logσ = LazyMap( SLEEFPirates.log, σ )
    PrecisionArray{S,T,N,X,L,typeof(σ⁻²),typeof(σ⁻¹),typeof(logσ)}( σ⁻², logσ, σ⁻¹ )
end
@inline function canonicalize_Σ(σ::S) where {S <: Integer}
    T = promote_type(S, Float64)
    canonicalize_Σ(convert(T, σ))
end
@inline canonicalize_Σ(σ²::UniformScaling{Bool}) = One()
@inline function canonicalize_Σ(σ²I::UniformScaling{T}) where {T <: Real}
    σ² = σ²I.λ
    σ⁻² = Base.FastMath.inv_fast(VectorizationBase.extract_data(σ²))
    Precision(
        σ⁻², Base.FastMat.div_fast(Base.log(σ⁻²),2), nothing
    )
end
@inline canonicalize_Σ(σ::AbstractPrecision) = σ
@inline ∂canonicalize_Σ(σ::AbstractPrecision) = σ
@inline Base.inv(σ::AbstractPrecision) = σ.σb⁻¹
@inline SIMDPirates.vinv(σ::AbstractPrecision) = σ.σ⁻¹
@inline loginvroot(σ::AbstractPrecision) = σ.logσ
@inline function Base.materialize(σ::Precision)
    Precision(
        Base.materialize(σ.σ⁻²),
        Base.materialize(σ.ld),
        Base.materialize(σ.∂k∂i)
    )
end
@inline function Base.materialize(σ::PrecisionArray{S,T,N,X,L,V,P,D}) where {S,T,N,X,L,V,P,D}
    PrecisionArray(
        Base.materialize(σ.σ⁻²),
        Base.materialize(σ.ld),
        Base.materialize(σ.∂k∂i)
    )
end
# @inline Base.inv(σ::Union{Precision,PrecisionArray}) = σ.σ⁻¹
# @inline SIMDPirates.vinv(σ::Union{Precision,PrecisionArray}) = σ.σ⁻¹
# @inline loginvroot(σ::Union{Precision,PrecisionArray}) = σ.logσ

@inline precision(λ::AbstractPrecision) = λ.σ⁻²
@inline precision(::One) = One()
@inline LinearAlgebra.logdet(λ::AbstractPrecision) = λ.ld
@inline StructuredMatrices.nlogdet(λ::AbstractPrecision) = Base.FastMath.sub_fast(λ.ld)
@inline LinearAlgebra.logdet(::One) = Zero()
@inline StructuredMatrices.nlogdet(::One) = Zero()
@inline ∂k∂i(λ::AbstractPrecision) = λ.∂k∂i
@inline ∂k∂i(λ::AbstractPrecision{T,Nothing}) where {T} = Base.FastMath.mul_fast(T(0.5), λ.σ⁻²) # Variance
@inline ∂ld∂i(λ::AbstractPrecision{T,Nothing}) where {T} = Base.FastMath.mul_fast(T(0.5), λ.σ⁻²) # Variance
@inline ∂ld∂i(λ::AbstractPrecision{T,P}) where {T,P<:Real} = λ.∂k∂i # St.Dev

"""
The normal distribution can be split into
Normal(args...) = Normal_kernel(args...) + logdet_coef(args...)*logdet(last(args))
"""
Normal(y) = Normal_kernel(y)
Normal(::Val{track}, y) where {track} = (first(track) == true ? Normal_kernel(y) : Zero())
∂Normal!(∂y, y) = ∂Normal_kernel!(∂y, y)
function Normal(::Val{track}, args...) where {track}
    fargs = Base.front(args)
    L = canonicalize_Σ(last(args))
    if last(track) == true
        tadd(vmul(logdet_coef(fargs..., L), nlogdet(L)), Normal_kernel(fargs...,L))
    else
        Normal_kernel(fargs...,L)
    end
end
@generated function ∂Normal!(args::Vararg{<:Any,N}) where {N}
    Nargs = N >>> 1
    # @assert N & 1 == 0 "Number of arguments is not odd." # separate ::StackPointer function
    if args[Nargs + (N & 1)] === Nothing
        quote
            @inbounds L = canonicalize_Σ(args[$N])
            @inbounds $(Expr(:call, :∂Normal_kernel!, [:(args[$n]) for n ∈ 1:N-1]..., L)) # if uninitialized, this initializes
        end
    else
        if isodd(N) # then there is a StackPointer
            quote
                @inbounds (sp,(∂Ldiv∂Σ, L)) = ∂canonicalize_Σ(first(args), args[$N])
                @inbounds (sp, kern) = $(Expr(:call, :∂Normal_kernel!, [:(args[$n]) for n ∈ 1:N-1]..., L)) # if uninitialized, this initializes
                @inbounds ∂TARGETdiv∂L = initialized(args[$(Nargs+(N&1))]) # so here we adjust type
                nlogdetL = ∂nlogdet!(∂TARGETdiv∂L, L, logdet_coef(args...))
                ∂TARGETdiv∂Σ = ∂TARGETdiv∂L # alias
                ReverseDiffExpressionsBase.RESERVED_INCREMENT_SEED_RESERVED!(∂TARGETdiv∂Σ, ∂Ldiv∂Σ, ∂TARGETdiv∂L) # adjust based on canonicalization
                (sp, tadd(nlogdetL, kern))
            end
        else # there is no StackPointer
            quote
                @inbounds (∂Ldiv∂Σ, L) = ∂canonicalize_Σ(args[$N])
                @inbounds kern = $(Expr(:call, :∂Normal_kernel!, [:(args[$n]) for n ∈ 1:N-1]..., L)) # if uninitialized, this initializes
                @inbounds ∂TARGETdiv∂L = initialized(args[$Nargs]) # so here we adjust type
                nlogdetL = ∂nlogdet!(∂TARGETdiv∂L, L, logdet_coef(args...))
                ∂TARGETdiv∂Σ = ∂TARGETdiv∂L # alias
                ReverseDiffExpressionsBase.RESERVED_INCREMENT_SEED_RESERVED!(∂TARGETdiv∂Σ, ∂Ldiv∂Σ, ∂TARGETdiv∂L) # adjust based on canonicalization
                tadd(nlogdetL, kern)
            end
        end
    end
end

push!(DISTRIBUTION_DIFF_RULES, :Normal)


text:
"""
	delaunayTriangulation(V::Matrix)::Array{Array{Int64,1},1}

Return highest level simplices of Delaunay triangulation.

# Examples 1D
```jldoctest

julia> V = [1. 2. 5. 6. 0. 7.];

julia> DT = AlphaStructures.delaunayTriangulation(V)
5-element Array{Array{Int64,1},1}:
 [1, 5]
 [1, 2]
 [2, 3]
 [3, 4]
 [4, 6]
```

# Examples 2D
```jldoctest

julia> V = [
 1. 2. 1. 2. ;
 0. 0. 1. 2.
];

julia> DT = AlphaStructures.delaunayTriangulation(V)
2-element Array{Array{Int64,1},1}:
 [1, 2, 3]
 [2, 3, 4]
```

# Examples 3D
```jldoctest

julia> V = [
 1. 2. 1. 2. ;
 0. 0. 1. 2. ;
 3. 1. 0. 2.
];

julia> DT = AlphaStructures.delaunayTriangulation(V)
1-element Array{Array{Int64,1},1}:
 [1, 2, 3, 4]
```
"""
function delaunayTriangulation(V::Matrix)::Array{Array{Int64,1},1}
	dim = size(V, 1)
	@assert dim > 0 "delaunayTriangulation: V do not contains points."
	@assert dim < 4 "delaunayTriangulation: Function not yet Programmed."

	if dim == 1
		vertices = vcat(V...)
		p = sortperm(vertices)
		upper_simplex = [[p[i],p[i+1]] for i=1:length(p)-1]
		return sort(sort.(upper_simplex))

	elseif dim == 2
		println("Delaunay triangulation of Triangulate package")
		upper_simplex = triangulate2D(V)
		return sort(upper_simplex)
	else
		println("Delaunay triangulation of Delaunay package")
		upper_simplex = DelaunayTriangulation(V)
		return sort(upper_simplex)
	end

end

"""
    delaunay_triangulation(points::Matrix) -> Array{Array{Int64,1},1}

Delaunay triangulation of points in d-dimensional Euclidean space.
Lar interface of Delaunay.jl.
"""
function DelaunayTriangulation(points::Matrix)

	centroid(points::Union{Matrix,Array{Float64,1}}) = (sum(points,dims=2)/size(points,2))[:,1]

    centroid = centroid(points)
    T = apply_matrix(traslation(-centroid...),points)
    mesh = delaunay(permutedims(T));
    DT = [mesh.simplices[c,:] for c in 1:size(mesh.simplices,1)]
    return sort.(DT)
end

"""
    delaunay_triangulation(points::Matrix) -> Array{Array{Int64,1},1}

Delaunay triangulation of points in d-dimensional Euclidean space.
Lar interface of Delaunay.jl.
"""
function triangulate2D(points::Matrix)
	triin = Triangulate.TriangulateIO()
	triin.pointlist = points
	(triout, vorout) = triangulate("Q", triin)
	trias = Array{Int64,1}[c[:] for c in eachcol(triout.trianglelist)]
	return sort.(trias)
end

# """
# 	delaunayMATLAB(V::Matrix)
#
# Delaunay triangulation algorithm in MATLAB.
# """
# function delaunayMATLAB(V::Matrix)
#
# 	dim = size(V,1)
# 	@assert dim <=3 "delaunayMATLAB: input points have invalid dimension."
#
# 	W = convert(Matrix,V')
# 	@mput W
# 	mat"DT = delaunay(W)"
# 	@mget DT
# 	DT = convert(Array{Int64,2},DT)
# 	DT = [DT[i,:] for i in 1:size(DT,1)]
#
# 	return DT
# end


text:
################################################################################
#  Copyright 2020, Tom Van Acker                                               #
################################################################################
# MultiStateSystems.jl                                                         #
# A Julia package to solve multi-state system models.                          #
# See http://github.com/timmyfaraday/MultiStateSystems.jl                      #
################################################################################

# network (abbr: ntw)
## structs
struct Network{I<:Int} <: AbstractNetwork{I}
    graph::_MG.DiMultigraph{I}

    props::PropDict

    cmp::Vector{PropDict}
    src::Vector{PropDict}
    usr::Vector{PropDict}

    clib::LibDict
    slib::LibDict
    ulib::LibDict
end

## constructors
################################################################################
# WARNING:  The empty constructor needs to be last in order to overwrite the   #
#           empty constructor created by other contructors, see: discourse -   #
#           keyword argument contructor breaks incomplete constructor.         #                                               #
################################################################################
"""
    Network()

An network constructor.

# Example
```julia-repl
julia> ntw = Network()
```
"""
function Network()
    graph = _MG.DiMultigraph(0)
    props = PropDict(:info => NetworkInfo(),
                     :msr  => Set{Symbol}())
    cmp, src, usr = PropDict[], PropDict[], PropDict[]
    clib, slib, ulib = LibDict(), LibDict(), LibDict()

    return Network(graph,props,cmp,src,usr,clib,slib,ulib)
end

## functions
is_directed(::Type{Network}) = true
is_directed(::Type{Network{Int}}) = true
is_directed(ntw::Network) = true

props(ntw::AbstractNetwork) = ntw.props

get_prop(ntw::AbstractNetwork, prop::Symbol) =
    haskey(props(ntw), prop) ? props(ntw)[prop] : ~ ;
has_prop(ntw::AbstractNetwork, prop::Symbol) = haskey(props(ntw), prop)
set_prop!(ntw::AbstractNetwork, prop::Symbol, value) =
    set_props!(ntw, PropDict(prop => value))
set_props!(ntw::AbstractNetwork, dict::PropDict) = merge!(ntw.props, dict)

nv(ntw::AbstractNetwork) = _MG.nv(ntw.graph)
add_vertex!(ntw::AbstractNetwork) = _MG.add_vertex!(ntw.graph)
add_vertex!(ntw::AbstractNetwork, node::Int) =
    !has_vertex(ntw,node) ? add_vertices!(ntw,node-nv(ntw)) : ~ ;
add_vertices!(ntw::AbstractNetwork, n::Int) = _MG.add_vertices!(ntw.graph, n)
has_vertex(ntw::AbstractNetwork, x...) = _MG.has_vertex(ntw.graph, x...)

ne(ntw::AbstractNetwork) = _MG.ne(ntw.graph)
add_edge!(ntw::AbstractNetwork, x...) = _MG.add_edge!(ntw.graph, x...)
has_edge(ntw::AbstractNetwork, x...) = _MG.has_edge(ntw.graph, x...)
mul_edge(ntw::AbstractNetwork,edge::Tuple{Int,Int}) =
    _MG.mul(ntw.graph, edge[1], edge[2])
mul_edge(graph::_MG.DiMultigraph,edge::Tuple{Int,Int}) =
    _MG.mul(graph, edge[1], edge[2])
update_lib!(type::Symbol,array::Array,lib::Dict) =
    haskey(lib,array[end][type]) ? push!(lib[array[end][type]],length(array)) :
                                   lib[array[end][type]] = [length(array)] ;

weights(ntw::AbstractNetwork) = _LG.weights(ntw.graph)
has_path(ntw::AbstractNetwork, s_node::Int, u_node::Int) =
    _LG.has_path(ntw.graph, s_node, u_node)
max_paths(graph::_MG.DiMultigraph) =
    _MG.nv(graph) + _MG.ne(graph, count_mul = true)
nodal_paths(graph::_MG.DiMultigraph, s_node::Int, u_node::Int) =
    sort(_LG.yen_k_shortest_paths(graph, s_node, u_node,
                             _LG.weights(graph), max_paths(graph)).paths)
mul_path(graph::_MG.DiMultigraph, npath::Array{Int,1}) =
    [1:mul_edge(graph,(npath[ni],npath[ni+1])) for ni in 1:length(npath)-1]

function get_extended_graph(ntw::AbstractNetwork, u_node::Int)
    graph = _MG.copy(ntw.graph)
    _MG.add_vertex!(graph)

    for src in ntw.src if has_path(ntw, src[:node], u_node)
        _MG.add_edge!(graph, _MG.nv(graph), src[:node])
    end end

    return graph, _MG.nv(graph)
end
function paths(ntw::AbstractNetwork, u_node::Int)
    npaths, cpaths = Vector{Int}[], Vector{Expr}[]

    graph, x_node = get_extended_graph(ntw, u_node)
    for npath in unique!(nodal_paths(graph, x_node, u_node))
        for mul in Iterators.product(mul_path(graph,npath)...)
            cpath = Vector{Expr}()
            # sources
            add_src_expr!(cpath, ntw, npath[2], mul[1])
            # components
            for nn in 2:length(npath)-1
                fr, to, ml = npath[nn], npath[nn+1], mul[nn]
                # fr-node cmp
                add_cmp_expr!(cpath, ntw, fr)
                # edge cmp
                edge = _MG.MultipleEdge(fr, to, ml)
                add_cmp_expr!(cpath, ntw, edge)
                # edge = _MG.MultipleEdge(to, fr, ml) # not necessary in a directed graph
                # add_cmp_expr!(cpath, ntw, edge)
            end
            # end-node cmp
            add_cmp_expr!(cpath, ntw, npath[end])

            push!(npaths,npath)
            push!(cpaths,cpath)
    end end
    return npaths, cpaths
end
function ntws(ntw::AbstractNetwork)
    cntr = 1
    ntws = Vector{AbstractNetwork}([ntw])
    while length(ntws) >= cntr
        push!(ntws,setdiff(get_ntw(ntws[cntr]),ntws)...)
        cntr += 1
    end
    for nn in ntws set_msr!(nn) end
    if any(x -> get_info(x,:dependent_sources),ntws)
        ni = findlast(x -> get_info(x,:dependent_sources),ntws)
        ntws[1].props[:source_ugf] = UGF(ntws[ni].props[:msr],ntws[ni].src[1][:std])
        for nn in ntws set_info!(nn,:dependent_sources,true) end
    end
    return reverse(ntws)
end

# info
## structs
mutable struct NetworkInfo{B<:Bool} <: AbstractInfo
    solved::B
    dependent_sources::B
    eval_dep::B

    # default constructor
    NetworkInfo() = new{Bool}(false,false,false)
end
mutable struct ComponentInfo{B<:Bool} <: AbstractInfo
    eval_dep::B
    eval_dep_id::B

    # default constructor
    ComponentInfo() = new{Bool}(false,false)
end
mutable struct UserInfo{B<:Bool} <: AbstractInfo
    eval_dep::B
    eval_dep_id::B

    # default constructor
    UserInfo() = new{Bool}(false,false)
end
mutable struct SourceInfo{B<:Bool} <: AbstractInfo
    eval_dep::B
    eval_dep_id::B

    # default constructor
    SourceInfo() = new{Bool}(false,false)
end

## functions
get_info(ntw::AbstractNetwork, info::Symbol) = getproperty(ntw.props[:info],info)
set_info!(ntw::AbstractNetwork, info::Symbol, value::Bool) =
    setproperty!(ntw.props[:info], info, value)
get_info(prt::PropDict, info::Symbol) = getproperty(prt[:info],info)
set_info!(prt::PropDict, info:: Symbol, value::Bool) =
    setproperty!(prt[:info], info, value)
set_info!(info::AbstractInfo, field::Symbol, value::Bool) =
    setproperty!(info, field, value)

# dependence
## source dependence
### functions
init_source_dep(ntw::AbstractNetwork, kwargs::Iterators.Pairs) =
    if haskey(kwargs, :dep)
        set_info!(ntw, :dependent_sources, kwargs[:dep])
    end

## evaluation dependence
### functions
init_eval_dep_src(ntw::AbstractNetwork, kwargs::Iterators.Pairs, Ns::Int) =
    if haskey(kwargs, :eval_dep)
        set_info!(ntw, :eval_dep, kwargs[:eval_dep])
        return length(ntw.src) .+ 1:Ns
    else
        return nothing
    end
init_eval_dep_cmp(ntw::AbstractNetwork, kwargs::Iterators.Pairs, Nc::Int) =
    if haskey(kwargs, :eval_dep)
        set_info!(ntw, :eval_dep, kwargs[:eval_dep])
        return length(ntw.cmp) .+ 1:Nc
    else
        return nothing
    end
function init_eval_dep_cmp(ntw::AbstractNetwork, Nc::Int)
    set_info!(ntw, :eval_dep, true)
    return length(ntw.cmp) .+ 1:Nc
end
init_eval_dep_usr(ntw::AbstractNetwork, kwargs::Iterators.Pairs, Nu::Int) =
    if haskey(kwargs, :eval_dep)
        set_info!(ntw, :eval_dep, kwargs[:eval_dep])
        return length(ntw.usr) .+ 1:Nu
    else
        return nothing
    end
set_eval_dep!(prop_dict::PropDict, ni::Int, eval_dep_ids) =
    if haskey(prop_dict, :eval_dep)
        prop_dict[:eval_dep_ids] = eval_dep_ids
        prop_dict[:eval_dep_id] = ni == 1 ? false : true ;
    end
set_eval_info!(info::AbstractInfo, prop_dict::PropDict) =
    if haskey(prop_dict, :eval_dep)
        for i_key in [:eval_dep,:eval_dep_id]
            set_info!(info, i_key, prop_dict[i_key])
            delete!(prop_dict, i_key)
        end
    end

# elements (abbr: elm)
## functions
elements(ntw::AbstractNetwork) = Iterators.flatten((ntw.cmp,ntw.src))
elm_range(elm::PropDict) =
    !get_info(elm,:eval_dep_id) ? length(elm[:ugf].val) : 1 ;
get_idx(ntw::AbstractNetwork) = [1:elm_range(elm) for elm in elements(ntw)]
get_idx_itr(ntw::AbstractNetwork) = Iterators.product(get_idx(ntw)...)
get_val(ntw::AbstractNetwork) = [elm[:ugf].val for elm in elements(ntw)]
get_prb(ntw::AbstractNetwork) = [elm[:ugf].prb for elm in elements(ntw)
                                               if  !get_info(elm,:eval_dep_id)]
get_ntw(ntw::AbstractNetwork) =
    [elm[:ntw][1] for elm in elements(ntw) if haskey(elm,:ntw)]

## component (abbr: cmp)
### functions
nc(ntw::AbstractNetwork) = length(ntw.cmp)
nc(ntw::AbstractNetwork, c_key::UIE) = length(ntw.clib[c_key])
cmp(ntw::AbstractNetwork) = ntw.cmp
cmp(ntw::AbstractNetwork, c_key::UIE) = ntw.cmp[ntw.clib[c_key][1]]
cmp_id(ntw::AbstractNetwork, c_key::UIE) = ntw.clib[c_key][1]
cmp_ix(ntw::AbstractNetwork, c_key::UIE) =
    get_info(cmp(ntw, c_key), :eval_dep) ?
        first(cmp(ntw, c_key)[:eval_dep_ids]) :
        cmp_id(ntw, c_key) ;
cmp_ids(ntw::AbstractNetwork) = 1:nc(ntw)
cmp_ids(ntw::AbstractNetwork, c_key::UIE) = ntw.clib[c_key]
cmp_keys(ntw::AbstractNetwork) = keys(ntw.clib)
cmp_expr(nc::Int) = :(val[$nc][idx[$nc]])
cmp_expr(nc::Int, ni::Int) = :(val[$nc][idx[$ni]])
add_cmp_expr!(cpath::Vector{Expr}, ntw::AbstractNetwork, c_key::UIE) =
    haskey(ntw.clib, c_key) ? push!(cpath, cmp_expr(cmp_id(ntw, c_key),
                                                    cmp_ix(ntw, c_key))) : ~ ;

"""
    add_components!(ntw::MultiStateSystems.AbstractNetwork; kwargs...)

Adds a single component to the network `ntw` and fills their corresponding
`PropDict` with the named arguments `kwargs`.

# Example
```julia-repl
julia> ntwᵖʷʳ = ntw()
julia> add_components!(ntwᵖʷʳ, edge = (1,2),
                               name = "cable 1",
                               std  = STD(power = [0u"MW",1500u"MW"],
                                          prob  = [0.2,0.8]))
```
"""
function add_component!(ntw::AbstractNetwork; kwargs...)
    (haskey(kwargs, :node) || haskey(kwargs, :edge)) || return false

    if haskey(kwargs,:node)
        add_vertex!(ntw, kwargs[:node])

        info = ComponentInfo()

        push!(ntw.cmp, PropDict(:info => info, kwargs...))

        update_lib!(:node, ntw.cmp, ntw.clib)
    end

    if haskey(kwargs,:edge)
        edge = kwargs[:edge]
        add_vertex!(ntw, maximum(edge))
        add_edge!(ntw, edge[1], edge[2])
        edge = _MG.MultipleEdge(edge[1], edge[2], mul_edge(ntw,edge))

        info = ComponentInfo()

        push!(ntw.cmp, Dict(:edge => edge, :info => info, reduce(kwargs, 1, excl=[:edge])...))

        update_lib!(:edge,ntw.cmp,ntw.clib)
    end

    return true
end
function add_component!(ntw::AbstractNetwork, node::Int, dict::Dict=PropDict())
    add_vertex!(ntw, node)

    info = ComponentInfo()
    set_eval_info!(info, dict)

    push!(ntw.cmp, PropDict(:node => node, :info => info, dict...))

    update_lib!(:node, ntw.cmp, ntw.clib)
end
function add_component!(ntw::AbstractNetwork, edge::Tuple{Int,Int}, dict::Dict=PropDict())
    add_vertex!(ntw, maximum(edge))
    add_edge!(ntw, edge[1], edge[2])
    edge = _MG.MultipleEdge(edge[1], edge[2], mul_edge(ntw, edge))

    info = ComponentInfo()
    set_eval_info!(info, dict)

    push!(ntw.cmp,Dict(:edge => edge, :info => info, dict...))

    update_lib!(:edge, ntw.cmp, ntw.clib)
end
function add_bidirectional_component!(ntw::AbstractNetwork; kwargs...)
    haskey(kwargs, :edge) || return false

    eval_dep_ids = init_eval_dep_cmp(ntw, 2)

    for (ne,edge) in enumerate([kwargs[:edge], reverse(kwargs[:edge])])
        prop_dict = reduce(kwargs, 1, excl=[:edge])
        prop_dict[:eval_dep] = true
        set_eval_dep!(prop_dict, ne, eval_dep_ids)
        add_component!(ntw, edge, prop_dict)
    end

    return true
end
function add_bidirectional_component!(ntw::AbstractNetwork, edge::Tuple{Int,Int}, dict::Dict=PropDict())
    if !haskey(dict,:eval_dep_ids)
        eval_dep_ids = init_eval_dep_cmp(ntw, 2)
    end
    for (ne,edge) in enumerate([edge, reverse(edge)])
        if !haskey(dict,:eval_dep_ids)
            dict[:eval_dep] = true
            set_eval_dep!(dict, ne, eval_dep_ids)
        end
        add_component!(ntw, edge, dict)
    end
end
"""
    add_components!(ntw::MultiStateSystems.AbstractNetwork; kwargs...)

Adds multiple components to the network `ntw` and fills their corresponding
`PropDict` with the named arguments `kwargs`. Either an uniform arguments is
given which holds for all components or an array is given whith specific
argument for each component.

# Example
```julia-repl
julia> ntwᵖʷʳ = ntw()
julia> add_components!(ntwᵖʷʳ, edge = [(1,2),(1,2),(2,3)],
                               name = ["cable 1","cable 2","cable 3"],
                               std  = [STD(power = [0u"MW",1500u"MW"],
                                           prob  = [0.2,0.8]),
                                       STD(power = [0u"MW",2000u"MW"],
                                           prob  = [0.4,0.6]),
                                       STD(power = [0u"MW",1800u"MW",4000u"MW"],
                                           prob = [0.1,0.2,0.7])])
```
"""
function add_components!(ntw::AbstractNetwork; kwargs...)
    (test(kwargs) && (haskey(kwargs, :node) || haskey(kwargs, :edge))) || return false

    if haskey(kwargs, :node)
        node, Nn = kwargs[:node], length(kwargs[:node])

        eval_dep_ids = init_eval_dep_cmp(ntw, kwargs, Nn)

        for ni in 1:Nn
            prop_dict = reduce(kwargs, ni, excl=[:node])
            set_eval_dep!(prop_dict, ni, eval_dep_ids)

            add_component!(ntw, node[ni], prop_dict)
    end end
    if haskey(kwargs,:edge)
        edge, Ne = kwargs[:edge], length(kwargs[:edge])

        eval_dep_ids = init_eval_dep_cmp(ntw, kwargs, Ne)

        for ni in 1:Ne
            prop_dict = reduce(kwargs, ni, excl=[:edge])
            set_eval_dep!(prop_dict, ni, eval_dep_ids)

            add_component!(ntw, edge[ni], prop_dict)
    end end
    return true
end
function add_bidirectional_components!(ntw::AbstractNetwork; kwargs...)
    (test(kwargs) && haskey(kwargs, :edge)) || return false

    edge, Ne = kwargs[:edge], length(kwargs[:edge])

    eval_dep_ids = init_eval_dep_cmp(ntw, kwargs, Ne)

    for ne in 1:Ne
        prop_dict = reduce(kwargs, ne, excl=[:edge])
        set_eval_dep!(prop_dict, ne, eval_dep_ids)

        add_bidirectional_component!(ntw, edge[ne], prop_dict)
    end
    return true
end

## source (abbr: src)
### functions
ns(ntw::AbstractNetwork) = length(ntw.src)
ns(ntw::AbstractNetwork,s_node::Int) = length(ntw.slib[s_node])
src(ntw::AbstractNetwork) = ntw.src
src(ntw::AbstractNetwork,s_node::Int, mul::Int) = ntw.src[ntw.slib[s_node][mul]]
src_id(ntw::AbstractNetwork, s_node::Int, mul::Int) =
    nc(ntw) + ntw.slib[s_node][mul]
src_ix(ntw::AbstractNetwork, s_node::Int, mul::Int) =
    get_info(src(ntw, s_node, mul), :eval_dep) ?
        nc(ntw) + first(src(ntw, s_node, mul)[:eval_dep_ids]) :
        nc(ntw) + ntw.slib[s_node][mul] ;
src_ids(ntw::AbstractNetwork, s_node::Int) = ntw.slib[s_node]
src_ids(ntw::AbstractNetwork) = 1:ns(ntw)
src_expr(ns::Int) = :(val[$ns][idx[$ns]])
src_expr(ns::Int, ni::Int) = :(val[$ns][idx[$ni]])
src_nodes(ntw::AbstractNetwork) = keys(ntw.slib)
add_src_expr!(cpath::Vector{Expr}, ntw::AbstractNetwork, node::Int, mul::Int) =
    haskey(ntw.slib, node) ? push!(cpath, src_expr( src_id(ntw, node, mul),
                                                    src_ix(ntw, node, mul))) : ~ ;

"""
    add_source!(ntw::MultiStateSystems.AbstractNetwork; kwargs...)

Adds a single source to the network `ntw` and fills their corresponding
`PropDict` with the named arguments `kwargs`.

# Example
```julia-repl
julia> ntwᵖʷʳ = ntw()
julia> stdᵍᵉⁿ = solvedSTD(prob = [0.1,0.2,0.7],
                          flow = [0.0u"MW",0.5u"MW",2.0u"MW"])
julia> add_source!(ntwᵖʷʳ, node = 1,
                           name = "generator 1",
                           std  = stdᵍᵉⁿ)
```
"""
function add_source!(ntw::AbstractNetwork; kwargs...)
    haskey(kwargs, :node) || return false

    add_vertex!(ntw, kwargs[:node])

    info = UserInfo()
    init_source_dep(ntw, kwargs)

    push!(ntw.src, Dict(:info => info, kwargs...))

    update_lib!(:node, ntw.src, ntw.slib)

    return true
end
function add_source!(ntw::AbstractNetwork, node::Int, dict::Dict=PropDict())
    add_vertex!(ntw, node)

    info = UserInfo()
    set_eval_info!(info, dict)

    push!(ntw.src,Dict(:node => node, :info => info, dict...))

    update_lib!(:node, ntw.src, ntw.slib)
end
"""
    add_sources!(ntw::MultiStateSystems.AbstractNetwork; kwargs...)

Adds multiple sources to the network `ntw` and fills their corresponding
`PropDict` with the named arguments `kwargs`. Either an uniform arguments is
given which holds for all components or an array is given whith specific
argument for each component.

# Example
```julia-repl
julia> ntwᵖʷʳ = ntw()
julia> stdᵍᵉⁿ = solvedSTD(prob = [0.1,0.2,0.7],
                          flow = [0.0u"MW",0.5u"MW",2.0u"MW"])
julia> add_sources!(ntwᵖʷʳ, node = 1:5,
                            std  = stdᵍᵉⁿ,
                            dep  = true)
```
"""

function add_sources!(ntw::AbstractNetwork; kwargs...)
    (test(kwargs) && haskey(kwargs,:node)) || return false

    node, Nn = kwargs[:node], length(kwargs[:node])

    init_source_dep(ntw, kwargs)
    eval_dep_ids = init_eval_dep_src(ntw, kwargs, Nn)

    for ni in 1:Nn
        prop_dict = reduce(kwargs, ni, excl=[:node])
        set_eval_dep!(prop_dict, ni, eval_dep_ids)

        add_source!(ntw, node[ni], prop_dict)
    end

    return true
end

# user (abbr: usr)
## functions
nu(ntw::AbstractNetwork) = length(ntw.usr)
nu(ntw::AbstractNetwork,u_node::Int) = length(ntw.ulib[u_node])
usr(ntw::AbstractNetwork) = ntw.usr
usr(ntw::AbstractNetwork,u_node::Int) = ntw.usr[ntw.ulib[u_node]]
usr_ids(ntw::AbstractNetwork) = 1:nu(ntw)
usr_ids(ntw::AbstractNetwork,u_node::Int) = ntw.ulib[u_node]
usr_nodes(ntw::AbstractNetwork) = keys(ntw.ulib)
"""
    add_user!(ntw::MultiStateSystems.AbstractNetwork; kwargs...)

Adds a single user to the network `ntw` and fills their corresponding
`PropDict` with the named arguments `kwargs`.

# Example
```julia-repl
julia> ntwᵖʷʳ = ntw()
julia> add_source!(ntwᵖʷʳ, node = 1,
                           ind  = [:EENS])
```
"""
function add_user!(ntw::AbstractNetwork; kwargs...)
    haskey(kwargs,:node) || return false

    add_vertex!(ntw,kwargs[:node])

    info = UserInfo()

    push!(ntw.usr,PropDict(:info => info, kwargs...))
    update_lib!(:node,ntw.usr,ntw.ulib)
    return true
end
function add_user!(ntw::AbstractNetwork, node::Int, dict::Dict=PropDict())
    add_vertex!(ntw,node)

    info = UserInfo()
    set_eval_info!(info, dict)

    push!(ntw.usr,PropDict(:node => node, :info => info, dict...))
    update_lib!(:node,ntw.usr,ntw.ulib)
end
"""
    add_users!(ntw::MultiStateSystems.AbstractNetwork; kwargs...)

Adds multiple users to the network `ntw` and fills their corresponding
`PropDict` with the named arguments `kwargs`. Either an uniform arguments is
given which holds for all components or an array is given whith specific
argument for each component.

# Example
```julia-repl
julia> ntwᵖʷʳ = ntw()
julia> add_sources!(ntwᵖʷʳ, node = [1,5,8],
                            ind  = [:EENS])
```
"""
function add_users!(ntw::AbstractNetwork; kwargs...)
    (test(kwargs) && haskey(kwargs,:node)) || return false

    node, Nn = kwargs[:node], length(kwargs[:node])

    eval_dep_ids = init_eval_dep_usr(ntw, kwargs, Nn)

    for ni in 1:Nn
        prop_dict = reduce(kwargs, ni, excl=[:node])
        set_eval_dep!(prop_dict, ni, eval_dep_ids)

        add_user!(ntw, node[ni], prop_dict)
    end
    return true
end


text:
# hierarchical multi-objetive LP
#

using Gurobi, Test

env = Gurobi.Env()
setparam!(env, "OutputFlag", 0)

A = [1. 1.]
lb = [0.; 0.]
b = [1.]
f = [1.; 1.]
model = gurobi_model(env, sense = :maximize, A = A, b = b, lb = lb, f = [0., 0.])

Gurobi.set_multiobj_n!(model, 2)
Gurobi.update_model!(model)
Gurobi.set_multiobj!(model, 1, [1., 1.], 10, 1.0)

Gurobi.set_multiobj!(model, 2, [1., 0.], 1, 1.0)
Gurobi.update_model!(model)

optimize(model)
@test get_solution(model) == [1., 0.]

Gurobi.set_multiobj!(model, 2, [0., 1.], 1, 1.0)
Gurobi.update_model!(model)

optimize(model)
@test get_solution(model) == [0., 1.]


text:
## Assumptions
## http://docs.sympy.org/0.7.2/modules/assumptions/index.html

"""
refine: http://docs.sympy.org/dev/modules/assumptions/refine.html
"""
refine(ex, assumpts...) = sympy_meth(:refine, ex, assumpts...)
export refine

"""

`ask`. Returns true, false or nothing

Example:

```
ask(Qinteger(x*y), And(Qinteger(x), Qinteger(y)))
## really slow isprime:
filter(x -> ask(Qprime(x)), [1:1000])
```

"""
ask(x::Sym, args...) = sympy_meth(:ask, x, args...)
export ask

## should we support & and | for (sym,sym) pairs? Not sure
## dependso on what happens to | (x==0) and ex == x usage for ==
## for now, we can combine terms logically with And, Or, Not...
## these are in logic module

## simple methods (x, args) -> y (y coercion happens via PyCall)
logic_sympy_methods = (
                     :And, :Or, :Not,
                     :Xor, :Nand, :Nor, :Implies,
                     :Equivalent,
                     :satisfiable
                     )

## We make a module Q to hold the assumptions
## this follows this page http://docs.sympy.org/0.7.5/_modules/sympy/assumptions/ask.html
module Q
import SymPy
import PyCall

##http://docs.sympy.org/dev/_modules/sympy/assumptions/ask.html#ask
Q_predicates = (:antihermitian,
                :bounded, :finite, # bounded deprecated
                :commutative,
#                :complex,
                :composite,
                :even,
                :extended_real,
                :hermitian,
                :imaginary,
                :infinitesimal,
                :infinity, :infinite, # :infinity deprecated
                :integer,
                :irrational,
                :rational,
                :algebraic,
                :transcendental,
                :negative,
                :nonzero, :zero,
                :positive,
                :prime,
                :real,
                :odd,
                :is_true,
                :nonpositive,
                :nonnegative,
                :symmetric,
                :invertible,
                :singular,
                :orthogonal,
                :unitary,
                :normal,
                :positive_definite,
                :upper_triangular,
                :lower_triangular,
                :diagonal,
                :triangular,
                :unit_triangular,
                :fullrank,
                :square,
                :real_elements,
                :complex_elements,
                :integer_elements)

for meth in Q_predicates
        nm = string(meth)
        @eval begin
            @doc """
`$($nm)`: a SymPy function.
The SymPy documentation can be found through: http://docs.sympy.org/latest/search.html?q=$($nm)
""" ->
            ($meth)(x) = PyCall.pycall(SymPy.sympy["Q"][$nm], SymPy.Sym, x)::SymPy.Sym
        end
    end
end

export Q


text:
# This file is auto-generated by AWSMetadata.jl
using AWS
using AWS.AWSServices: sagemaker_a2i_runtime

using AWS.Compat
using AWS.UUIDs
"""
    DeleteHumanLoop()

Deletes the specified human loop for a flow definition.

# Required Parameters
- `HumanLoopName`: The name of the human loop that you want to delete.

"""

delete_human_loop(HumanLoopName; aws_config::AWSConfig=global_aws_config()) = sagemaker_a2i_runtime("DELETE", "/human-loops/$(HumanLoopName)"; aws_config=aws_config)
delete_human_loop(HumanLoopName, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = sagemaker_a2i_runtime("DELETE", "/human-loops/$(HumanLoopName)", args; aws_config=aws_config)

"""
    DescribeHumanLoop()

Returns information about the specified human loop.

# Required Parameters
- `HumanLoopName`: The name of the human loop that you want information about.

"""

describe_human_loop(HumanLoopName; aws_config::AWSConfig=global_aws_config()) = sagemaker_a2i_runtime("GET", "/human-loops/$(HumanLoopName)"; aws_config=aws_config)
describe_human_loop(HumanLoopName, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = sagemaker_a2i_runtime("GET", "/human-loops/$(HumanLoopName)", args; aws_config=aws_config)

"""
    ListHumanLoops()

Returns information about human loops, given the specified parameters. If a human loop was deleted, it will not be included.

# Required Parameters
- `FlowDefinitionArn`: The Amazon Resource Name (ARN) of a flow definition.

# Optional Parameters
- `CreationTimeAfter`: (Optional) The timestamp of the date when you want the human loops to begin in ISO 8601 format. For example, 2020-02-24.
- `CreationTimeBefore`: (Optional) The timestamp of the date before which you want the human loops to begin in ISO 8601 format. For example, 2020-02-24.
- `MaxResults`: The total number of items to return. If the total number of available items is more than the value specified in MaxResults, then a NextToken is returned in the output. You can use this token to display the next page of results.
- `NextToken`: A token to display the next page of results.
- `SortOrder`: Optional. The order for displaying results. Valid values: Ascending and Descending.
"""

list_human_loops(FlowDefinitionArn; aws_config::AWSConfig=global_aws_config()) = sagemaker_a2i_runtime("GET", "/human-loops", Dict{String, Any}("FlowDefinitionArn"=>FlowDefinitionArn); aws_config=aws_config)
list_human_loops(FlowDefinitionArn, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = sagemaker_a2i_runtime("GET", "/human-loops", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("FlowDefinitionArn"=>FlowDefinitionArn), args)); aws_config=aws_config)

"""
    StartHumanLoop()

Starts a human loop, provided that at least one activation condition is met.

# Required Parameters
- `FlowDefinitionArn`: The Amazon Resource Name (ARN) of the flow definition associated with this human loop.
- `HumanLoopInput`: An object that contains information about the human loop.
- `HumanLoopName`: The name of the human loop.

# Optional Parameters
- `DataAttributes`: Attributes of the specified data. Use DataAttributes to specify if your data is free of personally identifiable information and/or free of adult content.
"""

start_human_loop(FlowDefinitionArn, HumanLoopInput, HumanLoopName; aws_config::AWSConfig=global_aws_config()) = sagemaker_a2i_runtime("POST", "/human-loops", Dict{String, Any}("FlowDefinitionArn"=>FlowDefinitionArn, "HumanLoopInput"=>HumanLoopInput, "HumanLoopName"=>HumanLoopName); aws_config=aws_config)
start_human_loop(FlowDefinitionArn, HumanLoopInput, HumanLoopName, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = sagemaker_a2i_runtime("POST", "/human-loops", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("FlowDefinitionArn"=>FlowDefinitionArn, "HumanLoopInput"=>HumanLoopInput, "HumanLoopName"=>HumanLoopName), args)); aws_config=aws_config)

"""
    StopHumanLoop()

Stops the specified human loop.

# Required Parameters
- `HumanLoopName`: The name of the human loop that you want to stop.

"""

stop_human_loop(HumanLoopName; aws_config::AWSConfig=global_aws_config()) = sagemaker_a2i_runtime("POST", "/human-loops/stop", Dict{String, Any}("HumanLoopName"=>HumanLoopName); aws_config=aws_config)
stop_human_loop(HumanLoopName, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = sagemaker_a2i_runtime("POST", "/human-loops/stop", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("HumanLoopName"=>HumanLoopName), args)); aws_config=aws_config)


text:
using SOtoZulip

include("configuration.jl")

db = getdb(SODB)
SOtoZulip.create_so_tables(db)


text:
using Unitful
Unitful.preferunits(u"cm,g,s"...)

function main()

    a = 9.81*u"m/s^2"
    v0 = 5u"m/s"

    x(t) = v0*t + 1/2*a*t^2
    v(t) = v0 + a*t

    t1 = 5*u"s"

    # displaying units
    println(a)
    println(v0)
    println(x(t1))
    println(v(t1))

    # checking unit math
    # println(x(t1) + v(t1))

    # converting units
    println(uconvert(u"yr", t1))

    # checking dimensions
    h = 10*u"ħ"
    println(h)
    println(typeof(h))

    # looking at SI units
    println(upreferred(h))

end

main()

max_stars_repo_name: nunesmelo/djs-office-hours

max_issues_repo_name: nunesmelo/djs-office-hours

max_forks_repo_name: nunesmelo/djs-office-hours


text:
# linear algebra matrix types
import Base.*

const Untracked = Input

# the following types have immutable data, but aren't constants
# so we create Input op for them, but don't assign argid

# WARNING: highly experimental, only partially implemented!

# function *(A::UpperTriangular{T1,A}, B::TArray{T2,1}) where {T1,T2,A<:TArray}
#     println("oooo!")
#     tA = record!(B.tape, Untracked, to_device(B.tape.device, A))
#     return tA * B
# end

# UpperTriangular(x::TArray) = record!()

# when this conversion happens?
# 1. in record_struct!
# 2. as part of subroutine -> record!(tape, Call, UpperTriangular, (x,))

# tA = record!(tape, Input, UpperTriangular(rand(3,3)))
# tA * B
# ^ works fine

# UpperTriangular{T,TArray{T,N}} -> untrack -> wrap into UpperTriangular -> track result
# record!(tape, Call, *, (A, tB))


text:
#Day6
answers = readlines("./src/Day6/input.txt")

allgroups = Array{Set,1}()
activegroup = Set()

for an in answers
    if length(an) == 0
        push!(allgroups, copy(activegroup))
        activegroup = Set()
    else
        push!(activegroup, split(an, "")...)
    end
end
if length(activegroup) > 0  # last group might not have ended on a blank line
    push!(allgroups, copy(activegroup))
end

sum(length.(allgroups))

allgroups = Array{Dict,1}()
activegroup = Dict{Char,Int}(collect('a':'z') .=> 0)
activegroup['_'] = 0

for an in answers
    if length(an) == 0
        push!(allgroups, copy(activegroup))
        activegroup = Dict{Char,Int}(collect('a':'z') .=> 0)
        activegroup['_'] = 0
    else
        for r in an
            activegroup[r] += 1
        end
        activegroup['_'] += 1
    end
end
if length(activegroup) > 0  # last group might not have ended on a blank line
    push!(allgroups, copy(activegroup))
end

vals = 0
for i in allgroups
    respondents = i['_']
    val = count(==(respondents), values(i)) - 1
    vals += val
end
vals


text:
using PyPlot
using celerite

# Compute time for Generalized Rybicki-Press with J Lorentzian components
# for benchmarking:

function time_julia_nexp()
omega = 2pi/12.203317
alpha_final = [1.0428542, 0.30345984]
alpha_imag = [0.0, 0.1]
beta_real_final = [0.1229159,0.09086397]
beta_imag_final = [0.0,omega]
#beta_imag_final = [0.0,0.0]
p_final = length(alpha_final)
p0_final = 0
for i=1:p_final
  if beta_imag_final[i] == 0.0
    p0_final +=1
  end
end

w0 =  0.03027
#nt = [16,16,32,64,128,256,512,1024,2048,4096,8192,16384,32768,65536,131072,262144,524288]
nt = 1024
#nt = [16]
nj = [1,2,4,8,16,32,64,128,256]
nnt = length(nj)
time_compute_final = zeros(nnt)
time_likelihood_final = zeros(nnt)
nrep = 1
for it=1:nnt
  n = nt
  j = nj[it]
  t = collect(linspace(0,n-1,n))
# White noise component:

# Use white noise for doing the timing test:
  y=randn(n)*sqrt(w0)
# Now use Ambikarasan O(N) method:
#  nex = (2p+1)*n-2p
#  width = 2p+3
#  m1 = p+1
#  aex::Array{Complex{Float64},2} = zeros(Complex{Float64},width,nex)
#  al_small::Array{Complex{Float64},2} = zeros(Complex{Float64},m1,nex)
#  indx::Vector{Int64} = collect(1:nex)
#  tic()
# See if I can look into types of this routine.
#  @code_warntype logdeta = lorentz_likelihood_hermitian_band_init(alpha,beta,w0,t,nex,aex,al_small,indx)
#  logdeta = lorentz_likelihood_hermitian_band_init(alpha,beta,w0,t,nex,aex,al_small,indx)
#  aex_full,bex_full  = lorentz_likelihood_hermitian(alpha,beta,w0,t,y)
#  eigval_aex = eigvals(aex_full)
# Now use the new slimmed-down real version:
  p_final = j
  p0_final = 0
  nex_final = (4(p_final-p0_final)+2p0_final+1)*(n-1)+1
  m1_final = 2(p_final-p0_final)+p0_final+2
  if p0_final == p_final
    m1_final = p0_final + 1
  end
  width_final = 2m1_final+1
  aex_final = zeros(Float64,width_final,nex_final)
  al_small_final = zeros(Float64,m1_final,nex_final)
  indx_final= collect(1:nex_final)
  logdeta_final=0.0
  alpha_final = 10.^(rand(j))
  alpha_imag = 10.^(rand(j))
  beta_real_final = 10.^(rand(j))
  beta_imag_final = 10.^(rand(j))
  tic()
#  @code_warntype   compile_matrix_symm(alpha_final,beta_real_final,beta_imag_final,w0,t,nex_final,aex_final,al_small_final,indx_final)
  for irep=1:nrep
    logdeta_final= compile_matrix_symm(alpha_final,alpha_imag,beta_real_final,beta_imag_final,w0,t,nex_final,aex_final,al_small_final,indx_final)
  end
  time_compute_final[it] = toq();
  bex = zeros(Float64,nex_final)
  log_like_final = 0.0
#  @code_warntype compute_likelihood(p_final,y,aex_final,al_small_final,indx_final,logdeta_final,bex)
  tic()
  for irep=1:nrep
    log_like_final= compute_likelihood(p_final,p0_final,y,aex_final,al_small_final,indx_final,logdeta_final,bex)
  end
  time_likelihood_final[it] = toq();
  println(n," final:   ",time_compute_final[it]," ",time_likelihood_final[it])

#  time_compute[it] = toq();
#  tic()
#  log_like = lorentz_likelihood_hermitian_band_save(p,y,aex,al_small,indx,logdeta)
#  time_likelihood[it] = toq();
#  println(n," ",time_compute[it]," ",time_likelihood[it])
#  println("Eigenvalues: ",eigval_aex)
#  read(STDIN,Char)
end

using PyPlot
loglog(nj,time_compute_final+time_likelihood_final,label="Compute time")
loglog(nj,(time_compute_final[8]+time_likelihood_final[8]).*(nj./nj[8]).^3,label="Compute time")
xlabel("J")
ylabel("Time [sec]")
legend(loc="upper left")
return nt,time_compute_final,time_likelihood_final
end


text:
using Serialization
using JSON3
using Colors, TestImages, ImageIO, ImageShow, FileIO

function get_saved_img_bin()
    pth = raw"C:\_LabView_projects\ZMQ\LabView0mqJl.jl\src\Julia\this-and-that\sers.jsr"

    nbd = nbdd = nothing
    open(pth, "r") do io
        nbd=deserialize(io)
        nbdd=deserialize(io)
    end;

    return (; bin_data=nbd, bindata_descr=[nbdd])
end

function img_from_saved()
    bin2nums(;get_saved_img_bin()...)[:rgbimg]
end

# function fia(arr, seq)
#     arr = UInt8.(arr)
#     sec = UInt8.(seq)
#     ls = length(seq)
#     for i in 1:(lastindex(arr)-ls)
#         if arr[i:i+ls-1] == seq
#             return i
#         end
#     end
#     return nothing
# end

# bd = LabView0mqJl.data2bin(img, "kwn", [1,3,2]); println(bd==bd0); println(Int.(bd0[1:15])); println(Int.(bd[1:15]))
# typeof(prw) = Array{UInt8,3}
# size(rwv) = (3, 256, 256)
# false
# [129, 42, 63, 170, 62, 102, 157, 60, 95, 153, 55, 92, 156, 55, 99]
# [129, 42, 63, 170, 62, 102, 157, 60, 95, 153, 55, 92, 156, 55, 99]

# using LabView0mqJl
#
# nbd, = LabView0mqJl.get_saved_img_bin();
# img = LabView0mqJl.img_from_saved()
# bd = LabView0mqJl.data2bin(img, "kwn")
# bd == nbd


text:
using Poly

# multiply arrays and double
instructions = [Instruction(:mult, :(out[i, j] += A[i, k] * B[k, j]), Set()),
                Instruction(:double, :(out[i, j] *= 2), Set())]

domains = [Domain(:i, 1, :(size(out, 1)), :(i += 1), Set(), []),
           Domain(:j, 1, :(size(out, 2)), :(j += 1), Set(), []),
           Domain(:k, 1, :(size(A, 2)), :(k += 1), Set(), [])]

kern = compile(LoopKernel(instructions, domains, [:out, :A, :B], [Array{Float64, 2}, Array{Float64, 2}, Array{Float64, 2}], []))

A = rand(10, 10)
B = rand(10, 10)
out = zeros(10, 10)

kern(out=out, A=A, B=B)

@show isapprox(out, A*B*2)

# add arrays with global offset
instructions = [Instruction(:offset, :(offset = 4), Set()),
                Instruction(:add, :(out[i, j] = A[i, j] + B[i, j] + offset), Set())]

domains = [Domain(:i, 1, :(size(out, 1)), :(i += 1), Set(), []),
           Domain(:j, 1, :(size(out, 2)), :(j += 1), Set(), [])]

kern = compile(LoopKernel(instructions, domains, [:out, :A, :B], [Array{Float64, 2}, Array{Float64, 2}, Array{Float64, 2}], []))

A = rand(10, 10)
B = rand(10, 10)
out = zeros(10, 10)

kern(out=out, A=A, B=B)

@show isapprox(out, A+B.+4)


text:
# functions to split data into train and test periods for the EMSx challenge

function write_site_file(filename, df; compressed=true)
    if compressed
        open(GzipCompressorStream, "$(filename).gz", "w") do stream
            CSV.write(stream, df)
        end
    else
        CSV.write("$(filename)", df)
    end
end

function read_site_file(file::String; kw...)
    return CSV.read(file, DataFrame)
end

function find_site_files(folder)
    files = readdir(folder)
    filter!(f -> occursin(r"^(\d+)\.csv(\.gz)?$",f), files)
    sort!(files)
    unique!(f -> match(r"(\d+)", f)[1], files) # since files is sorted, ".csv" will be chosen over ".csv.gz"
    return files
end

function initialize_data(path_to_data_folder::String,
                         path_to_test_periods_csv::String =
                             joinpath(DIR, "metadata", "test_periods.csv");
                         delete_files::Bool = true,
                         progress::Bool = true,
                         compressed = true)
    println()
    println("Splitting train and test datasets")
    mkpath(joinpath(path_to_data_folder, "test"))
    mkpath(joinpath(path_to_data_folder, "train"))
    full_site_files = find_site_files(path_to_data_folder)

    @showprogress for full_site_file in full_site_files
        train_test_split(full_site_file,
                         path_to_data_folder,
                         path_to_test_periods_csv;
                         compressed = compressed)
        delete_files && rm(joinpath(path_to_data_folder, full_site_file))
    end
end

function initialize_data_parallel(path_to_data_folder::String,
                                  path_to_test_periods_csv::String =
                                      joinpath(DIR, "metadata", "test_periods.csv");
                                  delete_files::Bool = true,
                                  progress::Bool = true,
                                  compressed = true)
    println()
    mkpath(joinpath(path_to_data_folder, "test"))
    mkpath(joinpath(path_to_data_folder, "train"))
    full_site_files = find_site_files(path_to_data_folder)

    if progress
        prog = ParallelProgress(length(full_site_files);
                                desc = "Splitting train and test datasets ")
        update!(prog, 0)
    end

    pmap(full_site_files) do full_site_file
        train_test_split(full_site_file,
                         path_to_data_folder,
                         path_to_test_periods_csv;
                         compressed = compressed)
        delete_files && rm(joinpath(path_to_data_folder, full_site_file))
        progress && next!(prog)
    end
    progress && finish!(prog)
    return
end

function train_test_split(site_file::String,
                          path_to_data_folder::String,
                          path_to_test_periods_csv::String;
                          compressed = true)

    test_periods = CSV.read(path_to_test_periods_csv, DataFrame)
    date_format = dateformat"y-m-dTH:M:S+z"
    site_id = parse(Int, match(r"(\d+)", site_file)[1])

    data = read_site_file(joinpath(path_to_data_folder, site_file))

    data.timestamp = DateTime.(data.timestamp, date_format)
    sort!(data, :timestamp)

    periods = test_periods[test_periods.site_id .== site_id, :test_periods][1]
    periods = [parse(Int, id) for id in split(periods[2:end-1], ",")]
    test_data = DataFrame()

    for period in periods
        df = data[data.period_id .== period, :]
        timestamp = df[1, :timestamp]
        history_span = timestamp-Day(1):Minute(15):timestamp-Minute(15)
        history = data[data.timestamp .∈ Ref(history_span), :]
        history[!, :period_id] = fill(period, 96)
        new_period = vcat(history, df)
        test_data = vcat(test_data, new_period)
    end

    test_file = joinpath(path_to_data_folder, "test", "$(site_id).csv")
    write_site_file(test_file, test_data; compressed=compressed)

    train_file = joinpath(path_to_data_folder, "train", "$(site_id).csv")
    train_data = data[data.period_id .∉ Ref(periods), :]
    write_site_file(train_file, train_data; compressed=compressed)

    return nothing
end


text:
# by lazarusA # HIDE
# using GLMakie # HIDE
using CairoMakie
CairoMakie.activate!() # HIDE
let
    U = LinRange(-pi, pi, 100) # 50
    V = LinRange(-pi, pi, 20)
    x1 = [cos(u) + .5 * cos(u) * cos(v)      for u in U, v in V]
    y1 = [sin(u) + .5 * sin(u) * cos(v)      for u in U, v in V]
    z1 = [.5 * sin(v)                        for u in U, v in V]
    x2 = [1 + cos(u) + .5 * cos(u) * cos(v)  for u in U, v in V]
    y2 = [.5 * sin(v)                        for u in U, v in V]
    z2 = [sin(u) + .5 * sin(u) * cos(v)      for u in U, v in V]
    fig = Figure(resolution =(900,600))
    ax = LScene(fig, scenekw = (camera = cam3d!, show_axis = true))
    wireframe!(x1,y1,z1, shading = false)
    axis = ax.scene[OldAxis]
    tstyle = axis[:names]
    tstyle[:textsize] = 16
    tstyle[:gap] = 8
    axis[:ticks][:textcolor] = :black
    axis[:ticks][:textsize] = 10
    fig[1,1] = ax
    fig
    save(joinpath(@__DIR__, "output", "wireTori.png"), fig, px_per_unit = 2.0) # HIDE
end
using Pkg # HIDE
Pkg.status(["CairoMakie"]) # HIDE

max_stars_repo_name: MisterBiggs/BeautifulMakie

max_issues_repo_name: MisterBiggs/BeautifulMakie

max_forks_repo_name: MisterBiggs/BeautifulMakie


text:
using Test
include(joinpath("..", "testhelpers.jl"))

@testset "MPIStateArrays reductions" begin
    runmpi(joinpath(@__DIR__, "basics.jl"))
    runmpi(joinpath(@__DIR__, "broadcasting.jl"))
    runmpi(joinpath(@__DIR__, "reductions.jl"))
    runmpi(joinpath(@__DIR__, "reductions.jl"), ntasks = 3)
    runmpi(joinpath(@__DIR__, "varsindex.jl"))
end


text:
# Raytracer.jl
# Raytracing for the generation of photorealistic images in Julia
# Copyright (c) 2021 Samuele Colombo, Paolo Galli

# Main unit test file

using Test

using Documenter, LinearAlgebra, StaticArrays

using Raytracer, Raytracer.Interpreter

##############
# Source code

@testset "Color" begin
    include("test_colors.jl")
end
@testset "HDR Image" begin
    include("test_hdrimage.jl")
end

@testset "PCG" begin
    include("test_pcg.jl")
end

@testset "Geometry" begin
    include("test_geometry.jl")
end
@testset "Transformations" begin
    include("test_transformations.jl")
end
@testset "Ray" begin
    include("test_ray.jl")
end
@testset "Cameras" begin
    include("test_cameras.jl")
end
@testset "Materials" begin
    include("test_materials.jl")
end
@testset "Shapes" begin
    include("test_shapes.jl")
end
@testset "World" begin
    include("test_world.jl")
end
@testset "Renderers" begin
    include("test_renderers.jl")
end
@testset "ImageTracer" begin
    include("test_imagetracer.jl")
end

################
# Documentation

DocMeta.setdocmeta!(Raytracer, :DocTestSetup, :(using Raytracer); recursive=true)
@testset "Docs" begin
    doctest(Raytracer, manual=false)
end


text:
using LightXML
using CSV
using DataFrames

"""
    dict_to_array(dict::Dict)

Given a dictionary, returns a tuple of arrays with the keys and values.
"""
function dict_to_array(dict::Dict)
    keys = []
    vals = []

    for (key, val) in dict
        push!(keys, key)
        push!(vals, val)
    end

    return keys, vals
end

"""
    strip_newline(val::String)

Replaces new line characters with spaces.
"""
function strip_newline(val::String)
    res = replace(val, "\n" => " ")
    res = replace(res, "\r" => " ")
    return res
end
strip_newline(::Missing) = missing

"""
    parse_MedlineDate(ml_dt::String)

Parses the contents of the MedlineDate element and returns a tuple of the year and month.
"""
function parse_MedlineDate(ml_dt::String)
    year = ""
    month = ""
    matches = split(ml_dt, " ", limit = 2)
    try
        year = matches[1]
        month = matches[2]
    catch
        @debug "Couldn't fully parse date: $ml_dt"
    end

    return year, month
end

"""
    parse_year(yr::String)

Parses the string year and returns an integer with the first year in range.
"""
function parse_year(yr::AbstractString)
    try
        Base.parse(Int64, yr[1:4])
    catch
        missing
    end
end

"""
    parse_month(mon::String)

Parses the string month (month or season) and returns an integer with the first month in range.
"""
function parse_month(mon::AbstractString)
    try
        trim_mon = lowercase(mon[1:3])

        if trim_mon == "jan" || trim_mon == "win"
            1
        elseif trim_mon == "feb"
            2
        elseif trim_mon == "mar"
            3
        elseif trim_mon == "apr" || trim_mon == "spr"
            4
        elseif trim_mon == "may"
            5
        elseif trim_mon == "jun"
            6
        elseif trim_mon == "jul" || trim_mon == "sum"
            7
        elseif trim_mon == "aug"
            8
        elseif trim_mon == "sep"
            9
        elseif trim_mon == "oct" || trim_mon == "fal"
            10
        elseif trim_mon == "nov"
            11
        elseif trim_mon == "dec"
            12
        else
            missing
        end

    catch
        missing
    end

end

"""
    parse_orcid(raw_orc::String)

Takes a string containing an ORC ID (url, 16 digit string) and returns a formatted ID (0000-1111-2222-3333).
"""
function parse_orcid(raw_orc::String)
    if occursin(r"^[0-9]{16}$", raw_orc)
        return "$(raw_orc[1:4])-$(raw_orc[5:8])-$(raw_orc[9:12])-$(raw_orc[13:16])"
    else
        reg = match(r"^.*([0-9]{4}-[0-9]{4}-[0-9]{4}-[0-9a-zA-Z]{4}).*$", raw_orc)

        return reg == nothing ? "PARSE_ERROR" : reg.captures[1]
    end
end

"""
    parse_author

Takes xml for author, and returns parsed elements
"""
function parse_author(xml::LightXML.XMLElement)

    first_name = missing :: Union{Missing,String}
    initials = missing :: Union{Missing,String}
    last_name = missing :: Union{Missing,String}
    suffix = missing :: Union{Missing,String}
    orcid = missing :: Union{Missing,String}
    collective = missing :: Union{Missing,String}
    affs = ""

    for names in child_elements(xml)
        names_name = name(names)
        if names_name == "LastName"
            last_name = content(names)
        elseif names_name == "ForeName"
            first_name = content(names)
        elseif names_name == "Initials"
            initials = content(names)
        elseif names_name == "Suffix"
            suffix = content(names)
        elseif names_name == "Identifier" && attribute(names, "Source") == "ORCID"
            orcid = parse_orcid(content(names))
        elseif names_name == "CollectiveName"
            collective = content(names)
        elseif names_name == "AffiliationInfo"
            for affil_info in child_elements(names)
                if name(affil_info) == "Affiliation"
                    affs *= content(affil_info) * "; "
                end
            end
        end
    end

    affil_str = (affs == "" ? missing : affs[1:end-2])

    return last_name, first_name, initials, suffix, orcid, collective, affil_str
end

# Note: If needed it could be further refactored to to that author, journal is a type
"""
    parse_articles(xml)

Parses a PubMedArticleSet that matches the NCBI-XML format
"""
function parse_articles(xml::LightXML.XMLElement)

    articles = get_elements_by_tagname(xml,"PubmedArticle")
    n_articles = length(articles)

    #basic
    pmid = Vector{Int64}(undef, n_articles)
    url = Vector{Union{Missing, String}}(undef, n_articles)
    title = Vector{String}(undef, n_articles)
    auth_cite = Vector{Union{Missing, String}}(undef, n_articles)
    pub_year = Vector{Union{Missing, Int64}}(undef, n_articles)
    pub_month = Vector{Union{Missing, Int64}}(undef, n_articles)
    pub_dt_desc = Vector{String}(undef, n_articles)
    journal_title = Vector{Union{Missing, String}}(undef, n_articles)
    journal_iso_abbrv = Vector{Union{Missing, String}}(undef, n_articles)
    journal_issn = Vector{Union{Missing, String}}(undef, n_articles)
    volume = Vector{Union{Missing, String}}(undef, n_articles)
    issue = Vector{Union{Missing, String}}(undef, n_articles)
    pages = Vector{Union{Missing, String}}(undef, n_articles)

    mesh_qual=Dict{Int64, String}() #qid, qdesc
    mesh_desc=Dict{Int64, String}() #did, ddesc

    #author_ref
    au_pmid = Vector{Int64}()
    au_lname = Vector{Union{String,Missing}}()
    au_fname = Vector{Union{String,Missing}}()
    au_inits = Vector{Union{String,Missing}}()
    au_suffix = Vector{Union{String,Missing}}()
    au_orcid = Vector{Union{String,Missing}}()
    au_coll = Vector{Union{String,Missing}}()
    au_affs = Vector{Union{String,Missing}}()

    #pub_type
    pt_pmid = Vector{Int64}()
    pt_uid = Vector{Int64}()
    pt_name = Vector{String}()

    #abstract_full
    af_pmid = Vector{Int64}()
    af_text = Vector{String}()

    #abstract_structured
    as_pmid = Vector{Int64}()
    as_nlm = Vector{Union{String,Missing}}()
    as_label = Vector{Union{String,Missing}}()
    as_text = Vector{String}()

    #mesh_heading
    mh_pmid = Vector{Int64}()
    mh_did = Vector{Int64}()
    mh_dmaj = Vector{Int8}()
    mh_qid = Vector{Union{Int64,Missing}}()
    mh_qmaj = Vector{Union{Int64,Missing}}()

    i = 1 ::Int64
    for article in articles

        # Initialize optional 1:1 article attributes
        this_issn = missing :: Union{Missing,String}
        this_volume = missing :: Union{Missing,String}
        this_issue = missing :: Union{Missing,String}
        this_journal_title = missing :: Union{Missing,String}
        this_iso_abbrv = missing :: Union{Missing,String}
        this_pages = missing :: Union{Missing,String}
        @inbounds auth_cite[i] = missing

        this_pmid = 0

        for tdat in child_elements(article)
            if name(tdat) == "MedlineCitation"
                for mc in child_elements(tdat)
                    if name(mc) == "PMID"
                        this_pmid = Base.parse(Int64, content(mc)) ::Int64
                        @inbounds url[i] = string("http://www.ncbi.nlm.nih.gov/pubmed/", this_pmid)
                        @inbounds pmid[i] = this_pmid
                    elseif name(mc) == "Article"
                        for a_info in child_elements(mc)

                            # Journal
                            if name(a_info) == "Journal"
                                for j_info in child_elements(a_info)

                                    # ISSN?
                                    if name(j_info) == "ISSN"
                                    this_issn = content(j_info)

                                    # JournalIssue
                                    elseif name(j_info) == "JournalIssue"
                                        for j_issue in child_elements(j_info)

                                            # Volume?
                                            if name(j_issue) == "Volume"
                                                this_volume = content(j_issue)

                                            # Issue?
                                            elseif name(j_issue) == "Issue"
                                                this_issue = content(j_issue)

                                            # PubDate
                                            elseif name(j_issue) == "PubDate"
                                                ystr = "" :: String
                                                mstr = "" :: String
                                                for j_pub_dt in child_elements(j_issue)
                                                    if name(j_pub_dt) == "MedlineDate"
                                                        ystr, mstr = parse_MedlineDate(content(j_pub_dt))
                                                    elseif name(j_pub_dt) == "Year"
                                                        ystr = content(j_pub_dt)
                                                    elseif name(j_pub_dt) == "Month" || name(j_pub_dt) == "Season"
                                                        mstr = content(j_pub_dt)
                                                    end
                                                end
                                                @inbounds pub_year[i] = parse_year(ystr)
                                                @inbounds pub_month[i] = parse_month(mstr)
                                                @inbounds pub_dt_desc[i] = ystr * (mstr == "" ? "" : " " * mstr)
                                            end # JournalIssue if
                                        end # JournalIssue For

                                    # Title?
                                    elseif name(j_info) == "Title"
                                        this_journal_title = strip_newline(content(j_info))

                                    # ISO Abbreviation?
                                    elseif name(j_info) == "ISOAbbreviation"
                                        this_iso_abbrv = content(j_info)
                                    end #Journal If
                                end # Journal For

                            # Article Title
                            elseif name(a_info) == "ArticleTitle"
                                @inbounds title[i] = strip_newline(content(a_info))

                            # Pagination?
                            elseif name(a_info) == "Pagination"
                                ml_pgn = missing :: Union{Missing,String}
                                start_page = missing :: Union{Missing,String}
                                end_page = missing :: Union{Missing,String}
                                for p_info in child_elements(a_info)
                                    if name(p_info) == "MedlinePgn"
                                        ml_pgn = content(p_info)
                                    elseif name(p_info) == "StartPage"
                                        start_page = content(p_info)
                                    elseif name(p_info) == "EndPage"
                                        end_page = content(p_info)
                                    end # p_info if
                                end # p_info for

                                if !ismissing(ml_pgn)
                                    this_pages = ml_pgn
                                else
                                    this_pages = (!ismissing(end_page) ? (start_page * "-" * end_page) : start_page)
                                end

                            # Abstract?
                            elseif name(a_info) == "Abstract"
                                # AbstractText+
                                # parse_abstracts!(a_info, abstract_structured, abstract_full, this_pmid)
                                abstract_full_text = "" :: String
                                abs_is_struct = false :: Bool
                                for abs in child_elements(a_info)
                                    if name(abs) == "AbstractText"
                                        if !abs_is_struct && has_attributes(abs)
                                            abs_is_struct = true
                                        end

                                        if abs_is_struct
                                            label = has_attribute(abs, "Label") ? attribute(abs,"Label") : missing :: Union{Missing,String}
                                            nlm_category = has_attribute(abs, "NlmCategory") ? attribute(abs, "NlmCategory") : missing :: Union{Missing,String}
                                            abs_text = strip_newline(content(abs))
                                            push!(as_pmid, this_pmid)
                                            push!(as_nlm, nlm_category)
                                            push!(as_label, label)
                                            push!(as_text, abs_text)
                                            abstract_full_text *= !ismissing(label) ? label * ": " * abs_text * " " : "NO LABEL: " * abs_text * " "
                                        else
                                            abstract_full_text *= strip_newline(content(abs)) * " "
                                        end
                                    end #abstract if
                                end #abstract for
                                push!(af_pmid, this_pmid)
                                push!(af_text, abstract_full_text[1:end-1])

                            # AuthorList?
                            elseif name(a_info) == "AuthorList"
                                # parse_authors!(a_info, authors, auth_cite, i, this_pmid)
                                this_auth_cite = "" :: String

                                for author in child_elements(a_info)
                                    last_name, first_name, initials, suffix, orcid, collective, affs = parse_author(author)

                                    this_auth_cite *= !ismissing(first_name) ? "$last_name, $first_name; " : (!ismissing(last_name) ? "$last_name; " : "")

                                    push!(au_pmid, this_pmid)
                                    push!(au_lname, last_name)
                                    push!(au_fname, first_name)
                                    push!(au_inits, initials)
                                    push!(au_suffix, suffix)
                                    push!(au_orcid, orcid)
                                    push!(au_coll, collective)
                                    push!(au_affs, strip_newline(affs))
                                end

                                @inbounds auth_cite[i] = (this_auth_cite == "" ? missing : this_auth_cite[1:end-2])

                            # PublicationTypeList
                            elseif name(a_info) == "PublicationTypeList"
                                #Publication Type+
                                # parse_pubtypes!(a_info, pub_type, this_pmid)
                                for pt in child_elements(a_info)
                                    desc = content(pt) :: String
                                    ui = attribute(pt, "UI") :: String
                                    uid = length(ui) > 1 ? Base.parse(Int64, ui[2:end]) : -1
                                    push!(pt_pmid, this_pmid)
                                    push!(pt_uid, uid)
                                    push!(pt_name, desc)
                                end
                            end # Article If
                        end # Article For

                    # MeshHeadingList?
                    elseif name(mc) == "MeshHeadingList"
                        for heading in child_elements(mc)
                        # parse_meshheadings!(mc, mesh_heading, mesh_desc, mesh_qual, this_pmid)
                            desc_uid = -1
                            desc_maj = -1
                            qual = missing :: Union{Missing,String}

                            for header in child_elements(heading)
                                header_name = name(header) :: String
                                if header_name == "DescriptorName"
                                    desc = content(header) :: String
                                    desc_maj = attribute(header, "MajorTopicYN") == "Y" ? 1 : 0
                                    desc_ui = attribute(header, "UI")
                                    desc_uid = length(desc_ui) > 1 ? Base.parse(Int, desc_ui[2:end]) : -1
                                    mesh_desc[desc_uid] = desc
                                elseif header_name == "QualifierName"
                                    qual = content(header)
                                    qual_maj = attribute(header, "MajorTopicYN") == "Y" ? 1 : 0 :: Int
                                    qual_ui = attribute(header, "UI")
                                    qual_uid = length(qual_ui) > 1 ? Base.parse(Int, qual_ui[2:end]) : -1

                                    mesh_qual[qual_uid] = qual

                                    push!(mh_pmid, this_pmid)
                                    push!(mh_did, desc_uid)
                                    push!(mh_dmaj, desc_maj)
                                    push!(mh_qid, qual_uid)
                                    push!(mh_qmaj, qual_maj)
                                end
                            end

                            if ismissing(qual)
                                qual_uid = missing
                                qual_maj = missing

                                push!(mh_pmid, this_pmid)
                                push!(mh_did, desc_uid)
                                push!(mh_dmaj, desc_maj)
                                push!(mh_qid, qual_uid)
                                push!(mh_qmaj, qual_maj)
                            end
                        end # heading for
                    end # MedlineCitation If
                end # medline citation if
            end # if nodename == MedlineCitation
        end # tdat For

        @inbounds journal_issn[i] = this_issn
        @inbounds volume[i] = this_volume
        @inbounds issue[i] = this_issue
        @inbounds journal_title[i] = this_journal_title
        @inbounds journal_iso_abbrv[i] = this_iso_abbrv
        @inbounds pages[i] = this_pages

        i += 1
    end # Document For

    dfs = Dict{String,DataFrame}()

    dfs["basic"] = DataFrame(pmid = pmid,
        pub_year = pub_year,
        pub_month = pub_month,
        pub_dt_desc = pub_dt_desc,
        title = title,
        authors = auth_cite,
        journal_title = journal_title,
        journal_ISSN = journal_issn,
        journal_volume = volume,
        journal_issue = issue,
        journal_pages = pages,
        journal_iso_abbreviation = journal_iso_abbrv,
        url = url)

    dfs["abstract_full"] = DataFrame(pmid = af_pmid,
        abstract_text = af_text)

    dfs["abstract_structured"] = DataFrame(pmid = as_pmid,
        nlm_category = as_nlm,
        label = as_label,
        abstract_text = as_text)

    dfs["author_ref"] = DataFrame(pmid = au_pmid,
        last_name = au_lname,
        first_name = au_fname,
        initials = au_inits,
        suffix = au_suffix,
        orcid = au_orcid,
        collective = au_coll,
        affiliation = au_affs)

    dfs["mesh_heading"] = DataFrame(pmid = mh_pmid,
        desc_uid = mh_did,
        desc_maj_status = mh_dmaj,
        qual_uid = mh_qid,
        qual_maj_status = mh_qmaj)

    did, dname = dict_to_array(mesh_desc)
    dfs["mesh_desc"] = DataFrame(uid = did,
        name = dname)

    qid, qname = dict_to_array(mesh_qual)
    dfs["mesh_qual"] = DataFrame(uid = qid,
        name = qname)

    dfs["pub_type"] = DataFrame(pmid = pt_pmid,
        uid = pt_uid,
        name = pt_name)

    return dfs
end

"""
    dfs_to_csv(dfs::Dict, path::String, [file_prefix::String])
Takes output of toDataFrames and writes to CSV files at the provided path and with the file prefix.
"""
function dfs_to_csv(dfs::Dict{String,DataFrame}, path::String, file_prefix::String="")
    [CSV.write(joinpath(path,"$file_prefix$k.csv"),v, missingstring = "NULL") for (k, v) in dfs]
    return nothing
end

"""
    remove_csvs(dfs, path, file_prefix)
Removes all of the CSV files associated with a dictionary of dataframes
"""
function remove_csvs(dfs::Dict{String,DataFrame}, path::String, file_prefix::String="")
    [rm(joinpath(path,"$(file_prefix)$k.csv"), force=true) for (k,v) in dfs]
    return nothing
end

"""
    remove_csvs(paths::Vector)
Removes all of the CSV files associated with an array of paths
"""
function remove_csvs(paths::Vector{String})
    [rm(path, force=true) for path in paths]
    return nothing
end