text:
#=
  Module To Fit GLM
=#

function glm(::RegularData, x::Array{T, 2}, y::Array{T}, distrib::AbstractDistribution, link::AbstractLink,
              solver::AbstractSolver; inverse::AbstractInverse = GETRIInverse(),
              offset::Array{T, 1} = Array{T, 1}(undef, 0),
              weights = Array{T, 1}(undef, 0), control::Control{T} = Control{T}()) where {T <: AbstractFloat}

  y, mu, weights = init!(distrib, y, weights)
  eta = linkfun(link, mu)

  coef = zeros(T, size(x)[2])
  coefold = zeros(T, size(x)[2])

  absErr::T = T(Inf)
  relErr::T = T(Inf)

  residuals::Array{T, 1} = zeros(T, size(y))
  dev::T = T(Inf)
  devold::T = T(Inf)

  iter::Int64 = 1

  n, p = size(x)
  converged::Bool = false
  badBreak::Bool = false

  doOffset = false
  if length(offset) != 0
    doOffset = true
  end
  doWeights = false
  if length(weights) != 0
    doWeights = true
  end

  Cov::Array{T, 2} = zeros(T, (p, p))
  xwx::Array{T, 2} = zeros(T, (p, p))
  R::Array{T, 2} = zeros(T, (p, p))
  xw::Array{T, 2} = zeros(T, (0, 0))
  while (relErr > control.epsilon)
    if control.printError
      println("Iteration: $iter")
    end
    z::Array{T, 1} = Z(link, y, mu, eta)
    if doOffset
      z .-= offset
    end

    w::Array{T, 1} = W(solver, distrib, link, mu, eta)
    # w::Array{T, 1} = W(distrib, link, mu, eta)
    if doWeights
      w .*= weights
    end

    # xw::Array{T, 2} = copy(x)
    # for j in 1:p
    #   for i in 1:n
    #     xw[i, j] = xw[i, j] * w[i]
    #   end
    # end
    # z .*= w

    # Overwrites xw in the QR case!
    # coef = solve(solver, xw, z)

    # coef = solve!(solver, Ref(R), Ref(xwx), Ref(xw),
    #               Ref(x), Ref(z), Ref(w), Ref(coef))

    R, xwx, xw, x, z, w, coef = solve!(solver, R, xwx, xw,
                                    x, z, w, coef)

    # println("xwx:\n", xwx)
    # println("coef:\n", coef)

    if(control.printCoef)
      println(coef)
    end

    eta = x * coef
    if doOffset
      eta .+= offset
    end
    mu = linkinv(link, eta)

    # println("Iteration: $iter");

    if length(weights) == 0
      residuals = devianceResiduals(distrib, mu, y)
    else
      residuals = devianceResiduals(distrib, mu, y, weights)
    end

    dev = sum(residuals)

    absErr = absoluteError(dev, devold)
    relErr = relativeError(dev, devold)
    # absErr = absoluteError(coef, coefold)
    # relErr = relativeError(coef, coefold)

    frac::T = T(1)
    coefdiff = coef .- coefold
    # Step control
    while (dev > (devold + control.epsilon*dev))

      if control.printError
        println("\tStep control")
        println("\tFraction: $frac")
        println("\tDeviance: $dev")
        println("\tAbsolute Error: $absErr")
        println("\tRelative Error: $relErr")
      end

      frac *= 0.5
      coef = coefold .+ (coefdiff .* frac)

      if(control.printCoef)
        println(coef)
      end

      # Abstract this block away into a function so
      # that it doesn't need to be repeated
      eta = x * coef
      if doOffset
        eta .+= offset
      end
      mu = linkinv(link, eta)

      if length(weights) == 0
        residuals = devianceResiduals(distrib, mu, y)
      else
        residuals = devianceResiduals(distrib, mu, y, weights)
      end

      dev = sum(residuals)

      absErr = absoluteError(dev, devold)
      relErr = relativeError(dev, devold)

      if frac < control.minstep
        error("Step control exceeded.")
      end
    end

    devold = dev
    coefold = coef

    # Cov = xw' * xw
    Cov = cov(solver, inverse, R, xwx, xw)

    if control.printError
      # println("Iteration: $iter")
      println("Deviance: $dev")
      println("Absolute Error: $absErr")
      println("Relative Error: $relErr")
    end

    if iter >= control.maxit
      println("Maximum number of iterations " *
                string(control.maxit) * " has been reached.")
      badBreak = true
      break
    end

    iter += 1

  end

  if badBreak
    converged = false
  else
    converged = true
  end

  phi::T = T(1)
  if (typeof(distrib) != BernoulliDistribution) | (typeof(distrib) != BinomialDistribution) | (typeof(distrib) != PoissonDistribution)
    phi = dev/(n - p)
    Cov .*= phi
  end

  return GLM(link, distrib, phi, coef, Cov, iter, relErr, absErr, converged,
             dev, residuals)
end

#============================== BLOCK GLM FUNCTION ==============================#
#=
  Valid only for solvers: GESVSolver, POSVSolver, SYSVSolver
=#
function glm(::Block1D, x::Array{Array{T, 2}, 1}, y::Array{Array{T, 2}, 1}, distrib::AbstractDistribution, link::AbstractLink,
              solver::AbstractSolver; inverse::AbstractInverse = GETRIInverse(),
              offset::Array{Array{T, 1}, 1} = Array{Array{T, 1}, 1}(undef, 0),
              weights = Array{Array{T, 1}, 1}(undef, 0), control::Control{T} = Control{T}()) where {T <: AbstractFloat}

  y, mu, weights = init!(distrib, y, weights)
  eta = linkfun(link, mu)

  coef = zeros(T, size(x[1])[2])
  coefold = zeros(T, size(x[1])[2])

  absErr::T = T(Inf)
  relErr::T = T(Inf)

  nBlocks::Int64 = length(mu)

  residuals::Array{Array{T, 1}, 1} = Array{Array{T, 1}, 1}(undef, 0)
  dev::T = T(Inf)
  devold::T = T(Inf)

  iter::Int64 = 1

  n = sum(map((y) -> size(y)[1], x))
  p = size(x[1])[2]
  converged::Bool = false
  badBreak::Bool = false

  doOffset = false
  if length(offset) != 0
    doOffset = true
  end
  doWeights = false
  if length(weights) != 0
    doWeights = true
  end

  Cov::Array{T, 2} = zeros(T, (p, p))
  xwx::Array{T, 2} = zeros(T, (p, p))
  R::Array{T, 2} = zeros(T, (p, p))
  xw::Array{Array{T, 2}, 1} = Array{Array{T, 2}, 1}(undef, 0)
  while (relErr > control.epsilon)
    if control.printError
      println("Iteration: $iter")
    end
    z::Array{Array{T, 1}, 1} = Z(link, y, mu, eta)
    if doOffset
      for i in 1:nBlocks
        z[i] .-= offset[i]
      end
    end

    w::Array{Array{T, 1}, 1} = W(solver, distrib, link, mu, eta)
    if doWeights
      for i in 1:nBlocks
        w[i] .*= weights[i]
      end
    end

    # println("Coefficient: ", coef)

    # println("Typeof coeff before: ", typeof(coef))
    R, xwx, xw, x, z, w, coef = solve!(solver, R, xwx, xw,
                                    x, z, w, coef)
    coef = reshape(coef, (p,))

    # println("xwx:\n", xwx)
    # println("coef:\n", coef)

    if(control.printCoef)
      println(coef)
    end

    eta = [ (x[i] * coef)[:, 1] for i in 1:nBlocks]
    if doOffset
      for i in 1:nBlocks
        eta[i] .+= offset[i]
      end
    end
    mu = linkinv(link, eta)

    if length(weights) == 0
      residuals = devianceResiduals(distrib, mu, y)
    else
      residuals = devianceResiduals(distrib, mu, y, weights)
    end

    dev = T(0)
    for i in 1:nBlocks
      dev += sum(residuals[i])
    end

    absErr = absoluteError(dev, devold)
    relErr = relativeError(dev, devold)
    # absErr = absoluteError(coef, coefold)
    # relErr = relativeError(coef, coefold)

    frac::T = T(1)
    coefdiff = coef .- coefold
    # Step control
    while (dev > (devold + control.epsilon*dev))

      if control.printError
        println("\tStep control")
        println("\tFraction: $frac")
        println("\tDeviance: $dev")
        println("\tAbsolute Error: $absErr")
        println("\tRelative Error: $relErr")
      end

      frac *= 0.5
      coef = coefold .+ (coefdiff .* frac)

      if(control.printCoef)
        println(coef)
      end

      # Abstract this block away into a function so
      # that it doesn't need to be repeated
      eta = [ (x[i] * coef)[:, 1] for i in 1:nBlocks]
      if doOffset
        for i in 1:nBlocks
          eta[i] .+= offset[i]
        end
      end
      mu = linkinv(link, eta)

      if length(weights) == 0
        residuals = devianceResiduals(distrib, mu, y)
      else
        residuals = devianceResiduals(distrib, mu, y, weights)
      end

      # dev = sum(residuals)
      dev = T(0)
      for i in 1:nBlocks
        dev += sum(residuals[i])
      end

      absErr = absoluteError(dev, devold)
      relErr = relativeError(dev, devold)

      if frac < control.minstep
        error("Step control exceeded.")
      end
    end

    devold = dev
    coefold = coef

    # Cov = xw' * xw
    Cov = cov(solver, inverse, R, xwx)

    if control.printError
      # println("Iteration: $iter")
      println("Deviance: $dev")
      println("Absolute Error: $absErr")
      println("Relative Error: $relErr")
    end

    if iter >= control.maxit
      println("Maximum number of iterations " *
                string(control.maxit) * " has been reached.")
      badBreak = true
      break
    end

    iter += 1

  end

  if badBreak
    converged = false
  else
    converged = true
  end

  phi::T = T(1)
  if (typeof(distrib) != BernoulliDistribution) | (typeof(distrib) != BinomialDistribution) | (typeof(distrib) != PoissonDistribution)
    phi = dev/(n - p)
    Cov .*= phi
  end

  return GLMBlock1D(link, distrib, phi, coef, Cov, iter, relErr, absErr, converged,
             dev, residuals)
end

#============================== BLOCK PARALLEL GLM FUNCTION ==============================#
#=
  Valid only for solvers: GESVSolver, POSVSolver, SYSVSolver
=#
function glm(matrixType::Block1DParallel, x::Array{Array{T, 2}, 1}, y::Array{Array{T, 2}, 1}, distrib::AbstractDistribution, link::AbstractLink,
              solver::AbstractSolver; inverse::AbstractInverse = GETRIInverse(),
              offset::Array{Array{T, 1}, 1} = Array{Array{T, 1}, 1}(undef, 0),
              weights = Array{Array{T, 1}, 1}(undef, 0), control::Control{T} = Control{T}()) where {T <: AbstractFloat}

  #= Set BLAS threads =#
  set_num_threads(1)
  y, mu, weights = init!(matrixType, distrib, y, weights)
  eta = linkfun(matrixType, link, mu)

  coef = zeros(T, size(x[1])[2])
  coefold = zeros(T, size(x[1])[2])

  absErr::T = T(Inf)
  relErr::T = T(Inf)

  nBlocks::Int64 = length(mu)

  residuals::Array{Array{T, 1}, 1} = Array{Array{T, 1}, 1}(undef, 0)
  dev::T = T(Inf)
  devold::T = T(Inf)

  iter::Int64 = 1

  n = sum(map((y) -> size(y)[1], x))
  p = size(x[1])[2]
  converged::Bool = false
  badBreak::Bool = false

  doOffset = false
  if length(offset) != 0
    doOffset = true
  end
  doWeights = false
  if length(weights) != 0
    doWeights = true
  end

  Cov::Array{T, 2} = zeros(T, (p, p))
  xwx::Array{T, 2} = zeros(T, (p, p))
  R::Array{T, 2} = zeros(T, (p, p))
  xw::Array{Array{T, 2}, 1} = Array{Array{T, 2}, 1}(undef, 0)

  while (relErr > control.epsilon)

    if control.printError
      println("Iteration: $iter")
    end
    z::Array{Array{T, 1}, 1} = Z(matrixType, link, y, mu, eta)
    if doOffset
      @threads for i in 1:nBlocks
        z[i] .-= offset[i]
      end
    end

    w::Array{Array{T, 1}, 1} = W(matrixType, solver, distrib, link, mu, eta)
    if doWeights
      @threads for i in 1:nBlocks
        w[i] .*= weights[i]
      end
    end

    # println("Coefficient: ", coef)

    # println("Typeof coeff before: ", typeof(coef))
    R, xwx, xw, x, z, w, coef = solve!(matrixType, solver, R,
                                    xwx, xw, x, z, w, coef)
    coef = reshape(coef, (p,))

    # println("xwx:\n", xwx)
    # println("coef:\n", coef)

    if(control.printCoef)
      println(coef)
    end

    # eta = [ (x[i] * coef)[:, 1] for i in 1:nBlocks]
    @threads for i in 1:nBlocks
      eta[i] = (x[i] * coef)[:, 1]
    end
    if doOffset
      @threads for i in 1:nBlocks
        eta[i] .+= offset[i]
      end
    end
    mu = linkinv(matrixType, link, eta)

    if length(weights) == 0
      residuals = devianceResiduals(matrixType, distrib, mu, y)
    else
      residuals = devianceResiduals(matrixType, distrib, mu, y, weights)
    end

    dev = T(0)
    # Parallel Reduction Required
    for i in 1:nBlocks
      dev += sum(residuals[i])
    end

    absErr = absoluteError(dev, devold)
    relErr = relativeError(dev, devold)
    # absErr = absoluteError(coef, coefold)
    # relErr = relativeError(coef, coefold)

    frac::T = T(1)
    coefdiff = coef .- coefold
    # Step control
    while (dev > (devold + control.epsilon*dev))

      if control.printError
        println("\tStep control")
        println("\tFraction: $frac")
        println("\tDeviance: $dev")
        println("\tAbsolute Error: $absErr")
        println("\tRelative Error: $relErr")
      end

      frac *= 0.5
      coef = coefold .+ (coefdiff .* frac)

      if(control.printCoef)
        println(coef)
      end

      # Abstract this block away into a function so
      # that it doesn't need to be repeated

      # eta = [ (x[i] * coef)[:, 1] for i in 1:nBlocks]
      @threads for i in 1:nBlocks
        eta[i] = (x[i] * coef)[:, 1]
      end
      if doOffset
        for i in 1:nBlocks
          eta[i] .+= offset[i]
        end
      end
      mu = linkinv(matrixType, link, eta)

      if length(weights) == 0
        residuals = devianceResiduals(matrixType, distrib, mu, y)
      else
        residuals = devianceResiduals(matrixType, distrib, mu, y, weights)
      end

      # dev = sum(residuals)
      dev = T(0)
      # Parallel Reduction Required
      for i in 1:nBlocks
        dev += sum(residuals[i])
      end

      absErr = absoluteError(dev, devold)
      relErr = relativeError(dev, devold)

      if frac < control.minstep
        error("Step control exceeded.")
      end
    end

    devold = dev
    coefold = coef

    # Cov = xw' * xw
    Cov = cov(solver, inverse, R, xwx)

    if control.printError
      # println("Iteration: $iter")
      println("Deviance: $dev")
      println("Absolute Error: $absErr")
      println("Relative Error: $relErr")
    end

    if iter >= control.maxit
      println("Maximum number of iterations " *
                string(control.maxit) * " has been reached.")
      badBreak = true
      break
    end

    iter += 1

  end

  if badBreak
    converged = false
  else
    converged = true
  end

  phi::T = T(1)
  if (typeof(distrib) != BernoulliDistribution) | (typeof(distrib) != BinomialDistribution) | (typeof(distrib) != PoissonDistribution)
    phi = dev/(n - p)
    Cov .*= phi
  end

  set_num_threads(nthreads())

  return GLMBlock1D(link, distrib, phi, coef, Cov, iter, relErr, absErr, converged,
             dev, residuals)
end

#============================== GRADIENT DESCENT ==============================#
#==============================================================================#

#============================== REGULAR DATA ==============================#
function glm(matrixType::RegularData, x::Array{T, 2}, y::Array{T},
              distrib::AbstractDistribution, link::AbstractLink,
              solver::AbstractGradientDescentSolver;
              inverse::AbstractInverse = GETRIInverse(),
              control::Control{T} = Control{T}(),
              calculateCovariance::Bool = true,
              doStepControl::Bool = true,
              offset::Array{T, 1} = Array{T, 1}(undef, 0),
              weights = Array{T, 1}(undef, 0)) where {T <: AbstractFloat}

  y, weights = init!(solver, distrib, y, weights)

  coef = zeros(T, size(x)[2])
  coefold = zeros(T, size(x)[2])

  eta = x * coef
  doOffset = false
  if length(offset) != 0
    doOffset = true
  end

  if doOffset
    eta .+= offset
  end
  mu = linkinv(link, eta)

  absErr::T = T(Inf)
  relErr::T = T(Inf)

  residuals::Array{T, 1} = zeros(T, size(y))
  dev::T = T(Inf)
  devold::T = T(Inf)

  iter::Int64 = 1

  n, p = size(x)
  converged::Bool = false
  badBreak::Bool = false

  doWeights = false
  if length(weights) != 0
    doWeights = true
  end

  vcov::Array{T, 2} = fill(T(1), (p, p))
  xwx::Array{T, 2} = zeros(T, (p, p))
  w::Array{T, 1} = zeros(T, (0,))

  while (relErr > control.epsilon)

    if control.printError
      println("Iteration: $iter")
    end

    solver = iteration(solver, iter)

    if typeof(solver) <: NesterovSolver
      coef .= NesterovModifier(solver, coef);

      eta = x * coef
      if doOffset
        eta .+= offset
      end
      mu = linkinv(link, eta)

      coef .= NesterovUnModifier(solver, coef);
    end

    solver, coef = solve!(matrixType, solver, distrib, link, y, x, mu, eta, coef)
    # println("Coef gradient descent regular data: ", coef, "\n")

    if(control.printCoef)
      println(coef)
    end

    eta = x * coef
    if doOffset
      eta .+= offset
    end
    mu = linkinv(link, eta)

    if length(weights) == 0
      residuals = devianceResiduals(distrib, mu, y)
    else
      residuals = devianceResiduals(distrib, mu, y, weights)
    end

    dev = sum(residuals)

    absErr = absoluteError(dev, devold)
    relErr = relativeError(dev, devold)

    frac::T = T(1)
    coefdiff = coef .- coefold
    # Step control
    while (doStepControl && (dev > (devold + control.epsilon*dev)))

      if control.printError
        println("\tStep control")
        println("\tFraction: $frac")
        println("\tDeviance: $dev")
        println("\tAbsolute Error: $absErr")
        println("\tRelative Error: $relErr")
      end

      frac *= 0.5
      coef = coefold .+ (coefdiff .* frac)

      if(control.printCoef)
        println(coef)
      end

      # Abstract this block away into a function so
      # that it doesn't need to be repeated
      eta = x * coef
      if doOffset
        eta .+= offset
      end
      mu = linkinv(link, eta)

      if length(weights) == 0
        residuals = devianceResiduals(distrib, mu, y)
      else
        residuals = devianceResiduals(distrib, mu, y, weights)
      end

      dev = sum(residuals)

      absErr = absoluteError(dev, devold)
      relErr = relativeError(dev, devold)

      if frac < control.minstep
        error("Step control exceeded.")
      end
    end

    devold = dev
    coefold = coef

    if control.printError
      println("Deviance: $dev")
      println("Absolute Error: $absErr")
      println("Relative Error: $relErr")
    end

    if iter >= control.maxit
      println("Maximum number of iterations " *
                string(control.maxit) * " has been reached.")
      badBreak = true
      break
    end

    iter += 1

  end

  if badBreak
    converged = false
  else
    converged = true
  end

  phi::T = T(1)
  if calculateCovariance
    w = W(solver, distrib, link, mu, eta)
    if doWeights
      w .*= weights
    end
    xwx = XWX(x, w)
    vcov = cov(solver, inverse, xwx)
    if (typeof(distrib) != BernoulliDistribution) | (typeof(distrib) != BinomialDistribution) | (typeof(distrib) != PoissonDistribution)
      phi = dev/(n - p)
      vcov .*= phi
    end
  end

  return GLM(link, distrib, phi, coef, vcov, iter, relErr, absErr, converged,
             dev, residuals)
end

#============================== BLOCK DATA GRADIENT DESCENT ==============================#

function glm(matrixType::Block1D, x::Array{Array{T, 2}, 1}, y::Array{Array{T, 2}, 1},
              distrib::AbstractDistribution, link::AbstractLink,
              solver::AbstractGradientDescentSolver;
              inverse::AbstractInverse = GETRIInverse(),
              control::Control{T} = Control{T}(),
              calculateCovariance::Bool = true,
              doStepControl::Bool = true,
              offset::Array{Array{T, 1}, 1} = Array{Array{T, 1}, 1}(undef, 0),
              weights = Array{Array{T, 1}, 1}(undef, 0)) where {T <: AbstractFloat}

  y, weights = init!(solver, distrib, y, weights)

  coef = zeros(T, size(x[1])[2])
  coefold = zeros(T, size(x[1])[2])

  nBlocks::Int64 = length(y)
  eta = [ (x[i] * coef)[:, 1] for i in 1:nBlocks]
  doOffset = false
  if length(offset) != 0
    doOffset = true
  end

  if doOffset
    for i in 1:nBlocks
      eta[i] .+= offset[i]
    end
  end
  mu = linkinv(link, eta)

  absErr::T = T(Inf)
  relErr::T = T(Inf)

  residuals::Array{Array{T, 1}, 1} = Array{Array{T, 1}, 1}(undef, 0)
  dev::T = T(Inf)
  devold::T = T(Inf)

  iter::Int64 = 1

  n = sum(map((y) -> size(y)[1], x))
  p = size(x[1])[2]
  converged::Bool = false
  badBreak::Bool = false

  doWeights = false
  if length(weights) != 0
    doWeights = true
  end

  vcov::Array{T, 2} = fill(T(1), (p, p))
  xwx::Array{T, 2} = zeros(T, (p, p))
  w::Array{Array{T, 1}, 1} = Array{Array{T, 1}, 1}(undef, 0)

  while (relErr > control.epsilon)

    if control.printError
      println("Iteration: $iter")
    end

    solver = iteration(solver, iter)

    if typeof(solver) <: NesterovSolver
      coef .= NesterovModifier(solver, coef);

      eta = [ (x[i] * coef)[:, 1] for i in 1:nBlocks]
      if doOffset
        for i in 1:nBlocks
          eta[i] .+= offset[i]
        end
      end
      mu = linkinv(link, eta)

      coef .= NesterovUnModifier(solver, coef);
    end

    solver, coef = solve!(matrixType, solver, distrib, link, y, x, mu, eta, coef)
    coef = reshape(coef, (p,))

    if(control.printCoef)
      println(coef)
    end

    eta = [ (x[i] * coef)[:, 1] for i in 1:nBlocks]
    if doOffset
      for i in 1:nBlocks
        eta[i] .+= offset[i]
      end
    end
    mu = linkinv(link, eta)

    if length(weights) == 0
      residuals = devianceResiduals(distrib, mu, y)
    else
      residuals = devianceResiduals(distrib, mu, y, weights)
    end

    dev = T(0)
    for i in 1:nBlocks
      dev += sum(residuals[i])
    end

    absErr = absoluteError(dev, devold)
    relErr = relativeError(dev, devold)
    # absErr = absoluteError(coef, coefold)
    # relErr = relativeError(coef, coefold)

    frac::T = T(1)
    coefdiff = coef .- coefold
    # Step control
    while (doStepControl && (dev > (devold + control.epsilon*dev)))

      if control.printError
        println("\tStep control")
        println("\tFraction: $frac")
        println("\tDeviance: $dev")
        println("\tAbsolute Error: $absErr")
        println("\tRelative Error: $relErr")
      end

      frac *= 0.5
      coef = coefold .+ (coefdiff .* frac)

      if(control.printCoef)
        println(coef)
      end

      # Abstract this block away into a function so
      # that it doesn't need to be repeated
      eta = [ (x[i] * coef)[:, 1] for i in 1:nBlocks]
      if doOffset
        for i in 1:nBlocks
          eta[i] .+= offset[i]
        end
      end
      mu = linkinv(link, eta)

      if length(weights) == 0
        residuals = devianceResiduals(distrib, mu, y)
      else
        residuals = devianceResiduals(distrib, mu, y, weights)
      end

      dev = T(0)
      for i in 1:nBlocks
        dev += sum(residuals[i])
      end

      absErr = absoluteError(dev, devold)
      relErr = relativeError(dev, devold)

      if frac < control.minstep
        error("Step control exceeded.")
      end
    end

    devold = dev
    coefold = coef

    if control.printError
      # println("Iteration: $iter")
      println("Deviance: $dev")
      println("Absolute Error: $absErr")
      println("Relative Error: $relErr")
    end

    if iter >= control.maxit
      println("Maximum number of iterations " *
                string(control.maxit) * " has been reached.")
      badBreak = true
      break
    end

    iter += 1

  end

  if badBreak
    converged = false
  else
    converged = true
  end

  phi::T = T(1)
  if calculateCovariance
    w = W(solver, distrib, link, mu, eta)
    if doWeights
      for i in 1:nBlocks
        w[i] .*= weights[i]
      end
    end
    xwx = XWX(x, w)
    vcov = cov(solver, inverse, xwx)
    if (typeof(distrib) != BernoulliDistribution) | (typeof(distrib) != BinomialDistribution) | (typeof(distrib) != PoissonDistribution)
      phi = dev/(n - p)
      vcov .*= phi
    end
  end

  return GLMBlock1D(link, distrib, phi, coef, vcov, iter, relErr, absErr, converged,
             dev, residuals)
end

#============================== BLOCK PARALLEL GRADIENT DESCENT ==============================#
#=
  Valid only for solvers: GESVSolver, POSVSolver, SYSVSolver
=#
function glm(matrixType::Block1DParallel, x::Array{Array{T, 2}, 1},
              y::Array{Array{T, 2}, 1}, distrib::AbstractDistribution,
              link::AbstractLink, solver::AbstractGradientDescentSolver;
              inverse::AbstractInverse = GETRIInverse(),
              control::Control{T} = Control{T}(),
              calculateCovariance::Bool = true,
              doStepControl::Bool = true,
              offset::Array{Array{T, 1}, 1} = Array{Array{T, 1}, 1}(undef, 0),
              weights = Array{Array{T, 1}, 1}(undef, 0)) where {T <: AbstractFloat}

  #= Set BLAS threads =#
  set_num_threads(1)

  y, weights = init!(solver, matrixType, distrib, y, weights)

  coef = zeros(T, size(x[1])[2])
  coefold = zeros(T, size(x[1])[2])

  nBlocks::Int64 = length(y)
  eta::Array{Array{T, 1}, 1} = Array{Array{T, 1}, 1}(undef, nBlocks)
  @threads for i in 1:nBlocks
    eta[i] = (x[i] * coef)[:, 1]
  end

  doOffset = false
  if length(offset) != 0
    doOffset = true
  end
  if doOffset
    @threads for i in 1:nBlocks
      eta[i] .+= offset[i]
    end
  end
  mu = linkinv(matrixType, link, eta)

  absErr::T = T(Inf)
  relErr::T = T(Inf)

  residuals::Array{Array{T, 1}, 1} = Array{Array{T, 1}, 1}(undef, 0)
  dev::T = T(Inf)
  devold::T = T(Inf)

  iter::Int64 = 1

  n = sum(map((y) -> size(y)[1], x))
  p = size(x[1])[2]
  converged::Bool = false
  badBreak::Bool = false

  doWeights = false
  if length(weights) != 0
    doWeights = true
  end

  vcov::Array{T, 2} = fill(T(1), (p, p))
  xwx::Array{T, 2} = zeros(T, (p, p))
  w::Array{Array{T, 1}, 1} = Array{Array{T, 1}, 1}(undef, 0)

  while (relErr > control.epsilon)

    if control.printError
      println("Iteration: $iter")
    end

    solver = iteration(solver, iter)

    if typeof(solver) <: NesterovSolver
      coef .= NesterovModifier(solver, coef);

      @threads for i in 1:nBlocks
        eta[i] = (x[i] * coef)[:, 1]
      end
      if doOffset
        @threads for i in 1:nBlocks
          eta[i] .+= offset[i]
        end
      end
      mu = linkinv(matrixType, link, eta)

      coef .= NesterovUnModifier(solver, coef);
    end

    # println("Coefficient: ", coef)

    solver, coef = solve!(matrixType, solver, distrib, link, y, x, mu, eta, coef)
    coef = reshape(coef, (p,))

    if(control.printCoef)
      println(coef)
    end

    @threads for i in 1:nBlocks
      eta[i] = (x[i] * coef)[:, 1]
    end
    if doOffset
      @threads for i in 1:nBlocks
        eta[i] .+= offset[i]
      end
    end
    mu = linkinv(matrixType, link, eta)

    if length(weights) == 0
      residuals = devianceResiduals(matrixType, distrib, mu, y)
    else
      residuals = devianceResiduals(matrixType, distrib, mu, y, weights)
    end

    dev = T(0)
    # Parallel Reduction Required
    for i in 1:nBlocks
      dev += sum(residuals[i])
    end

    absErr = absoluteError(dev, devold)
    relErr = relativeError(dev, devold)
    # absErr = absoluteError(coef, coefold)
    # relErr = relativeError(coef, coefold)

    frac::T = T(1)
    coefdiff = coef .- coefold
    # Step control
    while (doStepControl && (dev > (devold + control.epsilon*dev)))

      if control.printError
        println("\tStep control")
        println("\tFraction: $frac")
        println("\tDeviance: $dev")
        println("\tAbsolute Error: $absErr")
        println("\tRelative Error: $relErr")
      end

      frac *= 0.5
      coef = coefold .+ (coefdiff .* frac)

      if(control.printCoef)
        println(coef)
      end

      # Abstract this block away into a function so
      # that it doesn't need to be repeated

      # eta = [ (x[i] * coef)[:, 1] for i in 1:nBlocks]
      @threads for i in 1:nBlocks
        eta[i] = (x[i] * coef)[:, 1]
      end
      if doOffset
        for i in 1:nBlocks
          eta[i] .+= offset[i]
        end
      end
      mu = linkinv(matrixType, link, eta)

      if length(weights) == 0
        residuals = devianceResiduals(matrixType, distrib, mu, y)
      else
        residuals = devianceResiduals(matrixType, distrib, mu, y, weights)
      end

      dev = T(0)
      # Parallel Reduction Required
      for i in 1:nBlocks
        dev += sum(residuals[i])
      end

      absErr = absoluteError(dev, devold)
      relErr = relativeError(dev, devold)

      if frac < control.minstep
        error("Step control exceeded.")
      end
    end

    devold = dev
    coefold = coef

    if control.printError
      # println("Iteration: $iter")
      println("Deviance: $dev")
      println("Absolute Error: $absErr")
      println("Relative Error: $relErr")
    end

    if iter >= control.maxit
      println("Maximum number of iterations " *
                string(control.maxit) * " has been reached.")
      badBreak = true
      break
    end

    iter += 1

  end

  if badBreak
    converged = false
  else
    converged = true
  end

  phi::T = T(1)
  if calculateCovariance
    w = W(matrixType, solver, distrib, link, mu, eta)
    if doWeights
      @threads for i in 1:nBlocks
        w[i] .*= weights[i]
      end
    end
    xwx = XWX(x, w)
    vcov = cov(solver, inverse, xwx)
    if (typeof(distrib) != BernoulliDistribution) | (typeof(distrib) != BinomialDistribution) | (typeof(distrib) != PoissonDistribution)
      phi = dev/(n - p)
      vcov .*= phi
    end
  end

  set_num_threads(nthreads())

  return GLMBlock1D(link, distrib, phi, coef, vcov, iter, relErr, absErr, converged,
             dev, residuals)
end

max_stars_repo_name: dataPulverizer/glmSolver

max_issues_repo_name: dataPulverizer/glmSolver

max_forks_repo_name: dataPulverizer/glmSolver


text:
"""
`m = maxfinite(A)` calculates the maximum value in `A`, ignoring any values that are not finite (Inf or NaN).
"""
function maxfinite(A::AbstractArray{T}) where T
    ret = sentinel_max(T)
    for a in A
        ret = maxfinite_scalar(a, ret)
    end
    ret
end

function maxfinite(f, A::AbstractArray)
    ret = sentinel_max(typeof(f(first(A))))
    for a in A
        ret = maxfinite_scalar(f(a), ret)
    end
    ret
end

"""
`m = maxabsfinite(A)` calculates the maximum absolute value in `A`, ignoring any values that are not finite (Inf or NaN).
"""
function maxabsfinite(A::AbstractArray{T}) where T
    ret = sentinel_min(typeof(abs(A[1])))
    for a in A
        ret = maxfinite_scalar(abs(a), ret)
    end
    ret
end

minfinite_scalar(a::T, b::T) where {T} = isfinite(a) ? (b < a ? b : a) : b
maxfinite_scalar(a::T, b::T) where {T} = isfinite(a) ? (b > a ? b : a) : b
minfinite_scalar(a::T, b::T) where {T <: Union{Integer,FixedPoint}} = b < a ? b : a
maxfinite_scalar(a::T, b::T) where {T <: Union{Integer,FixedPoint}} = b > a ? b : a
minfinite_scalar(a, b) = minfinite_scalar(promote(a, b)...)
maxfinite_scalar(a, b) = maxfinite_scalar(promote(a, b)...)

function minfinite_scalar(c1::C, c2::C) where C <: AbstractRGB
    C(minfinite_scalar(c1.r, c2.r),
      minfinite_scalar(c1.g, c2.g),
      minfinite_scalar(c1.b, c2.b))
end
function maxfinite_scalar(c1::C, c2::C) where C <: AbstractRGB
    C(maxfinite_scalar(c1.r, c2.r),
      maxfinite_scalar(c1.g, c2.g),
      maxfinite_scalar(c1.b, c2.b))
end

sentinel_min(::Type{T}) where {T <: Union{Integer,FixedPoint}} = typemax(T)
sentinel_max(::Type{T}) where {T <: Union{Integer,FixedPoint}} = typemin(T)
sentinel_min(::Type{T}) where {T <: AbstractFloat} = convert(T, NaN)
sentinel_max(::Type{T}) where {T <: AbstractFloat} = convert(T, NaN)
sentinel_min(::Type{C}) where {C <: AbstractRGB} = _sentinel_min(C, eltype(C))
_sentinel_min(::Type{C},::Type{T}) where {C <: AbstractRGB,T} = (s = sentinel_min(T); C(s, s, s))
sentinel_max(::Type{C}) where {C <: AbstractRGB} = _sentinel_max(C, eltype(C))
_sentinel_max(::Type{C},::Type{T}) where {C <: AbstractRGB,T} = (s = sentinel_max(T); C(s, s, s))

difftype(::Type{T}) where {T <: Integer} = Int
difftype(::Type{T}) where {T <: Real} = Float32
difftype(::Type{Float64}) = Float64
difftype(::Type{CV}) where {CV <: Colorant} = difftype(CV, eltype(CV))
difftype(::Type{CV}, ::Type{T}) where {CV <: RGBA,T <: Real} = RGBA{Float32}
difftype(::Type{CV}, ::Type{Float64}) where {CV <: RGBA} = RGBA{Float64}
difftype(::Type{CV}, ::Type{Float64}) where {CV <: Gray} = Gray{Float64}
difftype(::Type{CV}, ::Type{T}) where {CV <: BGRA,T <: Real} = BGRA{Float32}
difftype(::Type{CV}, ::Type{Float64}) where {CV <: BGRA} = BGRA{Float64}
difftype(::Type{CV}, ::Type{T}) where {CV <: AbstractRGB,T <: Real} = RGB{Float32}
difftype(::Type{CV}, ::Type{Float64}) where {CV <: AbstractRGB} = RGB{Float64}

accum(::Type{T}) where {T <: Integer} = Int
accum(::Type{Float32})    = Float32
accum(::Type{T}) where {T <: Real} = Float64
accum(::Type{C}) where {C <: Colorant} = base_colorant_type(C){accum(eltype(C))}
function sumdiff(f, A::AbstractArray, B::AbstractArray)
    axes(A) == axes(B) || throw(DimensionMismatch("A and B must have the same axes"))
    T = promote_type(difftype(eltype(A)), difftype(eltype(B)))
    s = zero(accum(eltype(T)))
    for (a, b) in zip(A, B)
        x = convert(T, a) - convert(T, b)
        s += f(x)
    end
    s
end

"`s = ssd(A, B)` computes the sum-of-squared differences over arrays/images A and B"
ssd(A::AbstractArray, B::AbstractArray) = sumdiff(abs2, A, B)

"`s = sad(A, B)` computes the sum-of-absolute differences over arrays/images A and B"
sad(A::AbstractArray, B::AbstractArray) = sumdiff(abs, A, B)

function approx_difference(
        A::AbstractArray, B::AbstractArray,
        sigma::AbstractVector{T}=ones(ndims(A)),
        eps::AbstractFloat=1e-2
    ) where T <: Real

    if length(sigma) != ndims(A)
        error("Invalid sigma in test_approx_eq_sigma_eps. Should be ndims(A)-length vector of the number of pixels to blur.  Got: $sigma")
    end
    kern = KernelFactors.IIRGaussian(sigma)
    Af = imfilter(A, kern, NA())
    Bf = imfilter(B, kern, NA())
    diffscale = max(maxabsfinite(Af), maxabsfinite(Bf))
    d = sad(Af, Bf)
    return d / (length(Af) * diffscale)
end

ext: jl


text:
include("lpc.jl")

infile = "in.wav"
if length(ARGS) >= 1
    infile = ARGS[1]
end

pcm, fs = wavread(infile)
xs = pcm[:, 1]

out, n, frameTypes = lpc_encode_decode(xs, fs; excitationFile="excitation.wav")

outfile = splitext(infile)[1] * "_daftpunk.wav"
dump_voicing_labels(out, fs, n, frameTypes, outfile)


text:
using Mimi
using Distributions
using Query
using DataFrames
using IterTools
using DelimitedFiles

using Test

N = 100

mcs = @defmcs begin
    # Define random variables. The rv() is required to disambiguate an
    # RV definition name = Dist(args...) from application of a distribution
    # to an external parameter. This makes the (less common) naming of an
    # RV slightly more burdensome, but it's only required when defining
    # correlations or sharing an RV across parameters.
    rv(name1) = Normal(1, 0.2)
    rv(name2) = Uniform(0.75, 1.25)
    rv(name3) = LogNormal(20, 4)

    # define correlations
    name1:name2 = 0.7
    name1:name3 = 0.5

    # assign RVs to model Parameters
    share = Uniform(0.2, 0.8)
    sigma[:, Region1] *= name2
    sigma[2020:5:2050, (Region2, Region3)] *= Uniform(0.8, 1.2)

    depk = [Region1 => Uniform(0.08, 0.14),
            Region2 => Uniform(0.10, 1.50),
            Region3 => Uniform(0.10, 0.20)]

    # indicate which parameters to save for each model run. Specify
    # a parameter name or [later] some slice of its data, similar to the
    # assignment of RVs, above.
    save(grosseconomy.K, grosseconomy.YGROSS, emissions.E, emissions.E_Global)
end

Mimi.reset_compdefs()
include("../../examples/tutorial/02-two-region-model/main.jl")

m = model

# Optionally, user functions can be called just before or after a trial is run
function print_result(m::Model, mcs::MonteCarloSimulation, trialnum::Int)
    ci = Mimi.compinstance(m.mi, :emissions)
    value = Mimi.get_variable_value(ci, :E_Global)
    println("$(ci.comp_id).E_Global: $value")
end

output_dir = joinpath(tempdir(), "mcs")

generate_trials!(mcs, N, sampling=LHS, filename=joinpath(output_dir, "trialdata.csv"))

# Test that the proper number of trials were saved
d = readdlm(joinpath(output_dir, "trialdata.csv"), ',')
@test size(d)[1] == N+1 # extra row for column names

# Run trials 1:N, and save results to the indicated directory

Mimi.set_model!(mcs, m)

run_mcs(mcs, N, output_dir=output_dir)

# From MCS discussion 5/23/2018
# generate_trials(mcs, samples=load("foo.csv"))
#
# run_mcs(mcs, [:foo=>m1,:bar=>m2], output_vars=[:foo=>[:grosseconomy=>[:bar,:bar2,:bar3], :comp2=>:var2], :bar=>[]], N, output_dir="/tmp/Mimi")
# run_mcs(mcs, m1, output_vars=[:grosseconomy=>:asf, :foo=>:bar], N, output_dir="/tmp/Mimi")
# run_mcs(mm, output_vars=[(:base,:compname,:varname), (:)], N, output_dir="/tmp/Mimi")
# run_mcs(mcs, mcs, mm, output_vars=[:grosseconomy=>:asf, :foo=>:bar], N, output_dir="/tmp/Mimi")
# run_mcs(mcs, m, output_vars=[(:base,:compname,:varname), (:)], N, output_dir="/tmp/Mimi")

# run_mcs(mcs, m, N, post_trial_func=print_result, output_dir="/tmp/Mimi")

function show_E_Global(year::Int; bins=40)
    df = @from i in E_Global begin
             @where i.time == year
             @select i
             @collect DataFrame
        end
    histogram(df[:E_Global], bins=bins,
              title="Distribution of global emissions in $year",
              xlabel="Emissions")
end

#
# Test scenario loop capability
#
global loop_counter = 0

function outer_loop_func(mcs::MonteCarloSimulation, tup)
    global loop_counter
    loop_counter += 1

    # unpack tuple (better to use NT here?)
    (scen, rate) = tup
    @debug "outer loop: scen:$scen, rate:$rate"
end

function inner_loop_func(mcs::MonteCarloSimulation, tup)
    global loop_counter
    loop_counter += 1

    # unpack tuple (better to use NT here?)
    (scen, rate) = tup
    @debug "inner loop: scen:$scen, rate:$rate"
end

loop_counter = 0

run_mcs(mcs, N;
        output_dir=output_dir,
        scenario_args=[:scen => [:low, :high],
                       :rate => [0.015, 0.03, 0.05]],
        scenario_func=outer_loop_func,
        scenario_placement=Mimi.OUTER)

@test loop_counter == 6

loop_counter = 0

run_mcs(mcs, N;
        output_dir=output_dir,
        scenario_args=[:scen => [:low, :high],
                       :rate => [0.015, 0.03, 0.05]],
        scenario_func=inner_loop_func,
        scenario_placement=Mimi.INNER)

@test loop_counter == N * 6

function other_loop_func(mcs::MonteCarloSimulation, tup)
    global loop_counter
    loop_counter += 10
end

function pre_trial(mcs::MonteCarloSimulation, trialnum::Int, ntimesteps::Int, tup::Tuple)
    global loop_counter
    loop_counter += 1
end

loop_counter = 0

run_mcs(mcs, N;
        output_dir=output_dir,
        pre_trial_func=pre_trial,
        scenario_func=other_loop_func,
        scenario_args=[:scen => [:low, :high],
                       :rate => [0.015, 0.03, 0.05]])

@test loop_counter == 6 * N + 60

function post_trial(mcs::MonteCarloSimulation, trialnum::Int, ntimesteps::Int, tup::Union{Nothing,Tuple})
    global loop_counter
    loop_counter += 1

    m = mcs.models[1]
    # println("grosseconomy.share: $(m[:grosseconomy, :share])")
end

loop_counter = 0

N = 10
run_mcs(mcs, N;
        output_dir=output_dir,
        post_trial_func=post_trial)

@test loop_counter == N

N = 1000

# Test new values generated for RANDOM sampling

generate_trials!(mcs, N, sampling=RANDOM)
trial1 = copy(mcs.rvdict[:name1].dist.values)

generate_trials!(mcs, N, sampling=RANDOM)
trial2 = copy(mcs.rvdict[:name1].dist.values)

@test length(trial1) == length(trial2)
@test trial1 != trial2

# Test new values generated for LHS sampling

generate_trials!(mcs, N, sampling=LHS)
trial1 = copy(mcs.rvdict[:name1].dist.values)

generate_trials!(mcs, N, sampling=LHS)
trial2 = copy(mcs.rvdict[:name1].dist.values)

@test length(trial1) == length(trial2)
@test trial1 != trial2


text:
include("src/document.jl")

import Base.Filesystem

build_assets()

run(`pdflatex fibonacci_prime_factors_scomposition.tex`)

max_stars_repo_name: manuelelucchi/Fibonacci

max_issues_repo_name: manuelelucchi/Fibonacci

max_forks_repo_name: manuelelucchi/Fibonacci


text:
using JLD2
using Plots

include("show.jl")
include("errors.jl")
include("interpolate.jl")

include("../deflector_def.jl")

paramplots = Vector{Plots.Plot}(undef, 0)
fieldBIEplots = Vector{Plots.Plot}(undef, 0)
fieldRCWAplots = Vector{Plots.Plot}(undef, 0)
fieldDeltaRCWAplots = Vector{Plots.Plot}(undef, 0)
fielderrorplots = Vector{Plots.Plot}(undef, 0)
errorbydistanceplots = Vector{Plots.Plot}(undef, 0)
L2errors = Vector{Any}(undef, 0)
RCWAfielderrorplots = Vector{Plots.Plot}(undef, 0)
RCWAerrorbydistanceplots = Vector{Plots.Plot}(undef, 0)
RCWAL2errors = Vector{Any}(undef, 0)

jldopen("scripts/deflector.jld2", "r") do file
    for e in exponents
        N = 2^e
        prefix = "deflector/k$(k)/n$(N)/i$(i)/"
        BIEsol = file[prefix * "BIE"]
        _RCWAsol = file[prefix * "RCWA"]
        _DeltaRCWAsol = file[prefix * "DeltaRCWA"]
        x, y₁, y₂ = BIEsol.x, BIEsol.y₁, BIEsol.y₂
        RCWAsol = interpolate(x, y₁, y₂, _RCWAsol)
        DeltaRCWAsol = interpolate(x, y₁, y₂, _DeltaRCWAsol)
        push!(paramplots, show_params(η, μ, RCWAsol.x))
        push!(fieldBIEplots, show_fields(BIEsol))
        push!(fieldRCWAplots, show_fields(RCWAsol))
        push!(fieldDeltaRCWAplots, show_fields(DeltaRCWAsol))
        BIEerrors = get_errors(BIEsol, DeltaRCWAsol)
        RCWAerrors = get_errors(RCWAsol, DeltaRCWAsol)
        push!(fielderrorplots, show_fields(BIEerrors))
        push!(errorbydistanceplots, show_errorbydistance(BIEerrors))
        push!(L2errors, errorL2(BIEerrors))
        push!(RCWAfielderrorplots, show_fields(RCWAerrors))
        push!(RCWAerrorbydistanceplots, show_errorbydistance(RCWAerrors))
        push!(RCWAL2errors, errorL2(RCWAerrors))
    end
end

for (j, p) in enumerate(paramplots)
    savefig(p, "scripts/analyses/paramplot$(j).png")
end
for (j, p) in enumerate(fieldBIEplots)
    savefig(p, "scripts/analyses/fieldBIEplot$(j).png")
end
for (j, p) in enumerate(fieldRCWAplots)
    savefig(p, "scripts/analyses/fieldRCWAplot$(j).png")
end
for (j, p) in enumerate(fieldDeltaRCWAplots)
    savefig(p, "scripts/analyses/fieldDeltaRCWAplot$(j).png")
end
for (j, p) in enumerate(fielderrorplots)
    savefig(p, "scripts/analyses/fielderrorplot$(j).png")
end
for (j, p) in enumerate(errorbydistanceplots)
    savefig(p, "scripts/analyses/errorbydistanceplot$(j).png")
end
for (j, p) in enumerate(RCWAfielderrorplots)
    savefig(p, "scripts/analyses/RCWAfielderrorplot$(j).png")
end
for (j, p) in enumerate(RCWAerrorbydistanceplots)
    savefig(p, "scripts/analyses/RCWAerrorbydistanceplot$(j).png")
end
savefig(show_errorL2(exponents, L2errors), "scripts/analyses/L2error.png")
savefig(show_errorL2(exponents, RCWAL2errors), "scripts/analyses/RCWAL2error.png")


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

#Period types
"""
    Dates.value(x::Period) -> Int64

For a given period, return the value associated with that period.  For example,
`value(Millisecond(10))` returns 10 as an integer.
"""
value(x::Period) = x.value

# The default constructors for Periods work well in almost all cases
# P(x) = new((convert(Int64,x))
# The following definitions are for Period-specific safety
for period in (:Year, :Quarter, :Month, :Week, :Day, :Hour, :Minute, :Second, :Millisecond, :Microsecond, :Nanosecond)
    period_str = string(period)
    accessor_str = lowercase(period_str)
    # Convenience method for show()
    @eval _units(x::$period) = " " * $accessor_str * (abs(value(x)) == 1 ? "" : "s")
    # periodisless
    @eval periodisless(x::$period, y::$period) = value(x) < value(y)
    # AbstractString parsing (mainly for IO code)
    @eval $period(x::AbstractString) = $period(Base.parse(Int64, x))
    # The period type is printed when output, thus it already implies its own typeinfo
    @eval Base.typeinfo_implicit(::Type{$period}) = true
    # Period accessors
    typs = period in (:Microsecond, :Nanosecond) ? ["Time"] :
           period in (:Hour, :Minute, :Second, :Millisecond) ? ["Time", "DateTime"] : ["Date", "DateTime"]
    reference = period === :Week ? " For details see [`$accessor_str(::Union{Date, DateTime})`](@ref)." : ""
    for typ_str in typs
        @eval begin
            @doc """
                $($period_str)(dt::$($typ_str)) -> $($period_str)

            The $($accessor_str) part of a $($typ_str) as a `$($period_str)`.$($reference)
            """ $period(dt::$(Symbol(typ_str))) = $period($(Symbol(accessor_str))(dt))
        end
    end
    @eval begin
        @doc """
            $($period_str)(v)

        Construct a `$($period_str)` object with the given `v` value. Input must be
        losslessly convertible to an [`Int64`](@ref).
        """ $period(v)
    end
end

#Print/show/traits
Base.print(io::IO, x::Period) = print(io, value(x), _units(x))
Base.show(io::IO, ::MIME"text/plain", x::Period) = print(io, x)
Base.show(io::IO, p::P) where {P<:Period} = print(io, P, '(', value(p), ')')
Base.zero(::Union{Type{P},P}) where {P<:Period} = P(0)
Base.one(::Union{Type{P},P}) where {P<:Period} = 1  # see #16116
Base.typemin(::Type{P}) where {P<:Period} = P(typemin(Int64))
Base.typemax(::Type{P}) where {P<:Period} = P(typemax(Int64))
Base.isfinite(::Union{Type{P}, P}) where {P<:Period} = true

# Default values (as used by TimeTypes)
"""
    default(p::Period) -> Period

Returns a sensible "default" value for the input Period by returning `T(1)` for Year,
Month, and Day, and `T(0)` for Hour, Minute, Second, and Millisecond.
"""
function default end

default(p::Union{T,Type{T}}) where {T<:DatePeriod} = T(1)
default(p::Union{T,Type{T}}) where {T<:TimePeriod} = T(0)

(-)(x::P) where {P<:Period} = P(-value(x))
==(x::P, y::P) where {P<:Period} = value(x) == value(y)
==(x::Period, y::Period) = (==)(promote(x, y)...)
Base.isless(x::P, y::P) where {P<:Period} = isless(value(x), value(y))
Base.isless(x::Period, y::Period) = isless(promote(x, y)...)

# Period Arithmetic, grouped by dimensionality:
for op in (:+, :-, :lcm, :gcd)
    @eval ($op)(x::P, y::P) where {P<:Period} = P(($op)(value(x), value(y)))
end

/(x::P, y::P) where {P<:Period} = /(value(x), value(y))
/(x::P, y::Real) where {P<:Period} = P(/(value(x), y))
div(x::P, y::P, r::RoundingMode) where {P<:Period} = div(value(x), value(y), r)
div(x::P, y::Real, r::RoundingMode) where {P<:Period} = P(div(value(x), Int64(y), r))

for op in (:rem, :mod)
    @eval begin
        ($op)(x::P, y::P) where {P<:Period} = P(($op)(value(x), value(y)))
        ($op)(x::P, y::Real) where {P<:Period} = P(($op)(value(x), Int64(y)))
    end
end

(*)(x::P, y::Real) where {P<:Period} = P(value(x) * y)
(*)(y::Real, x::Period) = x * y

(*)(A::Period, B::AbstractArray) = Broadcast.broadcast_preserving_zero_d(*, A, B)
(*)(A::AbstractArray, B::Period) = Broadcast.broadcast_preserving_zero_d(*, A, B)

# intfuncs
Base.gcdx(a::T, b::T) where {T<:Period} = ((g, x, y) = gcdx(value(a), value(b)); return T(g), x, y)
Base.abs(a::T) where {T<:Period} = T(abs(value(a)))
Base.sign(x::Period) = sign(value(x))

periodisless(::Period,::Year)        = true
periodisless(::Period,::Quarter)     = true
periodisless(::Year,::Quarter)       = false
periodisless(::Period,::Month)       = true
periodisless(::Year,::Month)         = false
periodisless(::Quarter,::Month)      = false
periodisless(::Period,::Week)        = true
periodisless(::Year,::Week)          = false
periodisless(::Quarter,::Week)       = false
periodisless(::Month,::Week)         = false
periodisless(::Period,::Day)         = true
periodisless(::Year,::Day)           = false
periodisless(::Quarter,::Day)        = false
periodisless(::Month,::Day)          = false
periodisless(::Week,::Day)           = false
periodisless(::Period,::Hour)        = false
periodisless(::Minute,::Hour)        = true
periodisless(::Second,::Hour)        = true
periodisless(::Millisecond,::Hour)   = true
periodisless(::Microsecond,::Hour)   = true
periodisless(::Nanosecond,::Hour)    = true
periodisless(::Period,::Minute)      = false
periodisless(::Second,::Minute)      = true
periodisless(::Millisecond,::Minute) = true
periodisless(::Microsecond,::Minute) = true
periodisless(::Nanosecond,::Minute)  = true
periodisless(::Period,::Second)      = false
periodisless(::Millisecond,::Second) = true
periodisless(::Microsecond,::Second) = true
periodisless(::Nanosecond,::Second)  = true
periodisless(::Period,::Millisecond) = false
periodisless(::Microsecond,::Millisecond) = true
periodisless(::Nanosecond,::Millisecond)  = true
periodisless(::Period,::Microsecond)      = false
periodisless(::Nanosecond,::Microsecond)  = true
periodisless(::Period,::Nanosecond)       = false

# return (next coarser period, conversion factor):
coarserperiod(::Type{P}) where {P<:Period} = (P, 1)
coarserperiod(::Type{Nanosecond})  = (Microsecond, 1000)
coarserperiod(::Type{Microsecond}) = (Millisecond, 1000)
coarserperiod(::Type{Millisecond}) = (Second, 1000)
coarserperiod(::Type{Second}) = (Minute, 60)
coarserperiod(::Type{Minute}) = (Hour, 60)
coarserperiod(::Type{Hour}) = (Day, 24)
coarserperiod(::Type{Day}) = (Week, 7)
coarserperiod(::Type{Month}) = (Year, 12)

# Stores multiple periods in greatest to least order by type, not values,
# canonicalized to eliminate zero periods, merge equal period types,
# and convert more-precise periods to less-precise periods when possible
"""
    CompoundPeriod

A `CompoundPeriod` is useful for expressing time periods that are not a fixed multiple of
smaller periods. For example, \"a year and a  day\" is not a fixed number of days, but can
be expressed using a `CompoundPeriod`. In fact, a `CompoundPeriod` is automatically
generated by addition of different period types, e.g. `Year(1) + Day(1)` produces a
`CompoundPeriod` result.
"""
struct CompoundPeriod <: AbstractTime
    periods::Array{Period, 1}
    function CompoundPeriod(p::Vector{Period})
        n = length(p)
        if n > 1
            sort!(p, rev=true, lt=periodisless)
            # canonicalize p by merging equal period types and removing zeros
            i = j = 1
            while j <= n
                k = j + 1
                while k <= n
                    if typeof(p[j]) == typeof(p[k])
                        p[j] += p[k]
                        k += 1
                    else
                        break
                    end
                end
                if p[j] != zero(p[j])
                    p[i] = p[j]
                    i += 1
                end
                j = k
            end
            n = i - 1 # new length
            p  = resize!(p, n)
        elseif n == 1 && value(p[1]) == 0
            p = Period[]
        end

        return new(p)
    end
end

"""
    Dates.periods(::CompoundPeriod) -> Vector{Period}

Return the `Vector` of `Period`s that comprise the given `CompoundPeriod`.

!!! compat "Julia 1.7"
    This function requires Julia 1.7 or later.
"""
periods(x::CompoundPeriod) = x.periods

"""
    CompoundPeriod(periods) -> CompoundPeriod

Construct a `CompoundPeriod` from a `Vector` of `Period`s. All `Period`s of the same type
will be added together.

# Examples
```jldoctest
julia> Dates.CompoundPeriod(Dates.Hour(12), Dates.Hour(13))
25 hours

julia> Dates.CompoundPeriod(Dates.Hour(-1), Dates.Minute(1))
-1 hour, 1 minute

julia> Dates.CompoundPeriod(Dates.Month(1), Dates.Week(-2))
1 month, -2 weeks

julia> Dates.CompoundPeriod(Dates.Minute(50000))
50000 minutes
```
"""
CompoundPeriod(p::Vector{<:Period}) = CompoundPeriod(Vector{Period}(p))

CompoundPeriod(t::Time) = CompoundPeriod(Period[Hour(t), Minute(t), Second(t), Millisecond(t),
                                                Microsecond(t), Nanosecond(t)])

CompoundPeriod(p::Period...) = CompoundPeriod(Period[p...])

"""
    canonicalize(::CompoundPeriod) -> CompoundPeriod

Reduces the `CompoundPeriod` into its canonical form by applying the following rules:

* Any `Period` large enough be partially representable by a coarser `Period` will be broken
  into multiple `Period`s (eg. `Hour(30)` becomes `Day(1) + Hour(6)`)
* `Period`s with opposite signs will be combined when possible
  (eg. `Hour(1) - Day(1)` becomes `-Hour(23)`)

# Examples
```jldoctest
julia> Dates.canonicalize(Dates.CompoundPeriod(Dates.Hour(12), Dates.Hour(13)))
1 day, 1 hour

julia> Dates.canonicalize(Dates.CompoundPeriod(Dates.Hour(-1), Dates.Minute(1)))
-59 minutes

julia> Dates.canonicalize(Dates.CompoundPeriod(Dates.Month(1), Dates.Week(-2)))
1 month, -2 weeks

julia> Dates.canonicalize(Dates.CompoundPeriod(Dates.Minute(50000)))
4 weeks, 6 days, 17 hours, 20 minutes
```
"""
canonicalize(x::Period) = canonicalize(CompoundPeriod(x))
function canonicalize(x::CompoundPeriod)
    # canonicalize Periods by pushing "overflow" into a coarser period.
    p = x.periods
    n = length(p)
    if n > 0
        pc = sizehint!(Period[], n)
        P = typeof(p[n])
        v = value(p[n])
        i = n - 1
        while true
            Pc, f = coarserperiod(P)
            if i > 0 && typeof(p[i]) == P
                v += value(p[i])
                i -= 1
            end
            v0 = f == 1 ? v : rem(v, f)
            v0 != 0 && push!(pc, P(v0))
            if v != v0
                P = Pc
                v = div(v - v0, f)
            elseif i > 0
                P = typeof(p[i])
                v = value(p[i])
                i -= 1
            else
                break
            end
        end
        p = reverse!(pc)
        n = length(p)
    else
        return x
    end

    # reduce the amount of mixed positive/negative Periods.
    if n > 0
        pc = sizehint!(Period[], n)
        i = n
        while i > 0
            j = i

            # Determine sign of the largest period in this group which
            # can be converted into via coarserperiod.
            last = Union{}
            current = typeof(p[i])
            while i > 0 && current != last
                if typeof(p[i]) == current
                    i -= 1
                end
                last, current = current, coarserperiod(current)[1]
            end
            s = sign(value(p[i + 1]))

            # Adjust all the periods in the group based upon the
            # largest period sign.
            P = typeof(p[j])
            v = 0
            while j > i
                Pc, f = coarserperiod(P)
                if j > 0 && typeof(p[j]) == P
                    v += value(p[j])
                    j -= 1
                end
                v0 = f == 1 ? v : mod(v, f * s)
                v0 != 0 && push!(pc, P(v0))
                if v != v0
                    P = Pc
                    v = div(v - v0, f)
                elseif j > 0
                    P = typeof(p[j])
                    v = 0
                else
                    break
                end
            end
        end
        p = reverse!(pc)
    end

    return CompoundPeriod(p)
end

Base.convert(::Type{CompoundPeriod}, x::Period) = CompoundPeriod(Period[x])
function Base.string(x::CompoundPeriod)
    if isempty(x.periods)
        return "empty period"
    else
        s = ""
        for p in x.periods
            s *= ", " * string(p)
        end
        return s[3:end]
    end
end
Base.show(io::IO,x::CompoundPeriod) = print(io, string(x))

# E.g. Year(1) + Day(1)
(+)(x::Period,y::Period) = CompoundPeriod(Period[x, y])
(+)(x::CompoundPeriod, y::Period) = CompoundPeriod(vcat(x.periods, y))
(+)(y::Period, x::CompoundPeriod) = x + y
(+)(x::CompoundPeriod, y::CompoundPeriod) = CompoundPeriod(vcat(x.periods, y.periods))
# E.g. Year(1) - Month(1)
(-)(x::Period, y::Period) = CompoundPeriod(Period[x, -y])
(-)(x::CompoundPeriod, y::Period) = CompoundPeriod(vcat(x.periods, -y))
(-)(x::CompoundPeriod) = CompoundPeriod(-x.periods)
(-)(y::Union{Period, CompoundPeriod}, x::CompoundPeriod) = (-x) + y

GeneralPeriod = Union{Period, CompoundPeriod}
(+)(x::GeneralPeriod) = x

(==)(x::CompoundPeriod, y::Period) = x == CompoundPeriod(y)
(==)(x::Period, y::CompoundPeriod) = y == x
(==)(x::CompoundPeriod, y::CompoundPeriod) = canonicalize(x).periods == canonicalize(y).periods

Base.isequal(x::CompoundPeriod, y::Period) = isequal(x, CompoundPeriod(y))
Base.isequal(x::Period, y::CompoundPeriod) = isequal(y, x)
Base.isequal(x::CompoundPeriod, y::CompoundPeriod) = isequal(x.periods, y.periods)

# Capture TimeType+-Period methods
(+)(a::TimeType, b::Period, c::Period) = (+)(a, b + c)
(+)(a::TimeType, b::Period, c::Period, d::Period...) = (+)((+)(a, b + c), d...)

function (+)(x::TimeType, y::CompoundPeriod)
    for p in y.periods
        x += p
    end
    return x
end
(+)(x::CompoundPeriod, y::TimeType) = y + x

function (-)(x::TimeType, y::CompoundPeriod)
    for p in y.periods
        x -= p
    end
    return x
end

# Fixed-value Periods (periods corresponding to a well-defined time interval,
# as opposed to variable calendar intervals like Year).
const FixedPeriod = Union{Week, Day, Hour, Minute, Second, Millisecond, Microsecond, Nanosecond}

# like div but throw an error if remainder is nonzero
function divexact(x, y)
    q, r = divrem(x, y)
    r == 0 || throw(InexactError(:divexact, Int, x/y))
    return q
end

# TODO: this is needed to prevent undefined Period constructors from
# hitting the deprecated construct-to-convert fallback.
(::Type{T})(p::Period) where {T<:Period} = convert(T, p)::T

# FixedPeriod conversions and promotion rules
const fixedperiod_conversions = [(:Week, 7), (:Day, 24), (:Hour, 60), (:Minute, 60), (:Second, 1000),
                                 (:Millisecond, 1000), (:Microsecond, 1000), (:Nanosecond, 1)]
for i = 1:length(fixedperiod_conversions)
    T, n = fixedperiod_conversions[i]
    N = Int64(1)
    for j = (i - 1):-1:1 # less-precise periods
        Tc, nc = fixedperiod_conversions[j]
        N *= nc
        vmax = typemax(Int64) ÷ N
        vmin = typemin(Int64) ÷ N
        @eval function Base.convert(::Type{$T}, x::$Tc)
            $vmin ≤ value(x) ≤ $vmax || throw(InexactError(:convert, $T, x))
            return $T(value(x) * $N)
        end
    end
    N = n
    for j = (i + 1):length(fixedperiod_conversions) # more-precise periods
        Tc, nc = fixedperiod_conversions[j]
        @eval Base.convert(::Type{$T}, x::$Tc) = $T(divexact(value(x), $N))
        @eval Base.promote_rule(::Type{$T}, ::Type{$Tc}) = $Tc
        N *= nc
    end
end

# other periods with fixed conversions but which aren't fixed time periods
const OtherPeriod = Union{Month, Quarter, Year}
let vmax = typemax(Int64) ÷ 12, vmin = typemin(Int64) ÷ 12
    @eval function Base.convert(::Type{Month}, x::Year)
        $vmin ≤ value(x) ≤ $vmax || throw(InexactError(:convert, Month, x))
        Month(value(x) * 12)
    end
end
Base.convert(::Type{Year}, x::Month) = Year(divexact(value(x), 12))
Base.promote_rule(::Type{Year}, ::Type{Month}) = Month

let vmax = typemax(Int64) ÷ 4, vmin = typemin(Int64) ÷ 4
    @eval function Base.convert(::Type{Quarter}, x::Year)
        $vmin ≤ value(x) ≤ $vmax || throw(InexactError(:convert, Quarter, x))
        Quarter(value(x) * 4)
    end
end
Base.convert(::Type{Year}, x::Quarter) = Year(divexact(value(x), 4))
Base.promote_rule(::Type{Year}, ::Type{Quarter}) = Quarter

let vmax = typemax(Int64) ÷ 3, vmin = typemin(Int64) ÷ 3
    @eval function Base.convert(::Type{Month}, x::Quarter)
        $vmin ≤ value(x) ≤ $vmax || throw(InexactError(:convert, Month, x))
        Month(value(x) * 3)
    end
end
Base.convert(::Type{Quarter}, x::Month) = Quarter(divexact(value(x), 3))
Base.promote_rule(::Type{Quarter}, ::Type{Month}) = Month

# fixed is not comparable to other periods, except when both are zero (#37459)
(==)(x::FixedPeriod, y::OtherPeriod) = iszero(x) & iszero(y)
(==)(x::OtherPeriod, y::FixedPeriod) = y == x

const zero_or_fixedperiod_seed = UInt === UInt64 ? 0x5b7fc751bba97516 : 0xeae0fdcb
const nonzero_otherperiod_seed = UInt === UInt64 ? 0xe1837356ff2d2ac9 : 0x170d1b00
otherperiod_seed(x::OtherPeriod) = iszero(value(x)) ? zero_or_fixedperiod_seed : nonzero_otherperiod_seed
# tons() will overflow for periods longer than ~300,000 years, implying a hash collision
# which is relatively harmless given how infrequent such periods should appear
Base.hash(x::FixedPeriod, h::UInt) = hash(tons(x), h + zero_or_fixedperiod_seed)
# Overflow can also happen here for really long periods (~8e17 years)
Base.hash(x::Year, h::UInt) = hash(12 * value(x), h + otherperiod_seed(x))
Base.hash(x::Quarter, h::UInt) = hash(3 * value(x), h + otherperiod_seed(x))
Base.hash(x::Month, h::UInt) = hash(value(x), h + otherperiod_seed(x))

function Base.hash(x::CompoundPeriod, h::UInt)
    isempty(x.periods) && return hash(0, h + zero_or_fixedperiod_seed)
    for p in x.periods
        h = hash(p, h)
    end
    return h
end

Base.isless(x::FixedPeriod, y::OtherPeriod) = throw(MethodError(isless, (x, y)))
Base.isless(x::OtherPeriod, y::FixedPeriod) = throw(MethodError(isless, (x, y)))

Base.isless(x::Period, y::CompoundPeriod) = CompoundPeriod(x) < y
Base.isless(x::CompoundPeriod, y::Period) = x < CompoundPeriod(y)
Base.isless(x::CompoundPeriod, y::CompoundPeriod) = tons(x) < tons(y)
# truncating conversions to milliseconds, nanoseconds and days:
# overflow can happen for periods longer than ~300,000 years
toms(c::Nanosecond)  = div(value(c), 1000000)
toms(c::Microsecond) = div(value(c), 1000)
toms(c::Millisecond) = value(c)
toms(c::Second)      = 1000 * value(c)
toms(c::Minute)      = 60000 * value(c)
toms(c::Hour)        = 3600000 * value(c)
toms(c::Day)         = 86400000 * value(c)
toms(c::Week)        = 604800000 * value(c)
toms(c::Month)       = 86400000.0 * 30.436875 * value(c)
toms(c::Quarter)     = 86400000.0 * 91.310625 * value(c)
toms(c::Year)        = 86400000.0 * 365.2425 * value(c)
toms(c::CompoundPeriod) = isempty(c.periods) ? 0.0 : Float64(sum(toms, c.periods))
tons(x)              = toms(x) * 1000000
tons(x::Microsecond) = value(x) * 1000
tons(x::Nanosecond)  = value(x)
tons(c::CompoundPeriod) = isempty(c.periods) ? 0.0 : Float64(sum(tons, c.periods))
days(c::Millisecond) = div(value(c), 86400000)
days(c::Second)      = div(value(c), 86400)
days(c::Minute)      = div(value(c), 1440)
days(c::Hour)        = div(value(c), 24)
days(c::Day)         = value(c)
days(c::Week)        = 7 * value(c)
days(c::Year)        = 365.2425 * value(c)
days(c::Quarter)     = 91.310625 * value(c)
days(c::Month)       = 30.436875 * value(c)
days(c::CompoundPeriod) = isempty(c.periods) ? 0.0 : Float64(sum(days, c.periods))

max_stars_repo_name: udesou/julia

max_issues_repo_name: zhuangchaoxi/julia

max_forks_repo_name: zhuangchaoxi/julia


text:
# Hicks, G. A., & Ray, W. H. (1971).
# Approximation methods for optimal control synthesis.
# The Canadian Journal of Chemical Engineering, 49(4), 522-528.

function reference_tracking(; store_results=false::Bool)
    datadir = nothing
    if store_results
        datadir = generate_data_subdir(@__FILE__)
    end

    # irreversible reaction
    reaction = "irreversible"
    time = 20.0f0  # final time (θ in original paper)
    cf = 1.0f0
    Tf = 300.0f0
    Tc = 290.0f0
    J = 100.0f0
    α = 1.95f-4
    k10 = 300.0f0
    N = 25.2f0
    u_upper = 1500.0f0
    u_lower = 0.0f0
    y1s = 0.408126f0
    y2s = 3.29763f0
    us = 370.0f0

    # # reversible reaction
    # reaction = "reversible"
    # time = 1f0  # final time (θ in original paper)
    # cf = 1f0
    # Tf = 323f0
    # Tc = 326f0
    # J = 100f0
    # α = 1f0
    # k10 = 1.5f7
    # k20 = 1.5f10
    # N = 10f0
    # γ = 1.5f0
    # u_upper = 9.5f0
    # u_lower = 0f0
    # y1s = 0.433848f0
    # y2s = 0.659684f0
    # us = 3.234f0

    # adimensional constants
    yf = Tf / (J * cf)
    yc = Tc / (J * cf)

    # case 1
    # α1 = 1f6
    # α2 = 2f3
    # α3 = 1f-3

    # case 2
    # α1 = 1f6
    # α2 = 2f3
    # α3 = 0f0

    # case 3
    # α1 = 10f0
    # α2 = 1f0
    # α3 = 0.1f0

    # case 4
    # α1 = 10f0
    # α2 = 1f0
    # α3 = 0f0

    # custom case
    α1 = 1.0f0
    α2 = 1.0f1
    α3 = 1.0f-1

    # initial conditions and timepoints
    @show u0 = [1.0f0, yf]
    @show tspan = (0.0f0, time)
    Δt = 1f-2

    function system!(du, u, p, t, controller)

        # neural network outputs controls taken by the system
        y1, y2 = u
        c = controller(u, p)[1]  # controller output is a 1-element container

        # reaction rate
        r = k10 * y1 * exp(-N / y2)  # irreversible (case 1 & case 2)
        # r = k10 * y1 * exp(-N/y2) - k20 * exp(-γ*N/y2) * (1-y1)  # reversible (case 3 & case 4)

        # dynamics of the controlled system
        y1_prime = (1 - y1) / time - r
        y2_prime = (yf - y2) / time + r - α * c * (y2 - yc)

        # update in-place
        @inbounds begin
            du[1] = y1_prime
            du[2] = y2_prime
        end
    end

    # set arquitecture of neural network controller
    controller = FastChain(
        FastDense(2, 16, tanh_fast),
        FastDense(16, 16, tanh_fast),
        FastDense(16, 2),
        (x, p) -> [u_lower + (u_upper - u_lower) * sigmoid_fast(x[1])],  # controllers ∈ [u_lower, u_upper]
    )

    controlODE = ControlODE(controller, system!, u0, tspan; Δt)

    # model weights are destructured into a vector of parameters
    θ = initial_params(controller)

    # define objective function to optimize
    function loss(controlODE, params)

        # curious error with ROS3P()
        sol = solve(controlODE, params) |> Array # integrate ODE system

        sum_squares = 0.0f0
        for state in eachcol(sol)
            control = controller(state, params)
            sum_squares +=
                α1 * (state[1] - y1s)^2 + α2 * (state[2] - y2s)^2 + α3 * (control[1] - us)^2
        end
        return sum_squares * 0.01f0
    end

    # closures to comply with required interface
    loss(params) = loss(controlODE, params)
    function plotting_callback(params, loss)
        return plot_simulation(controlODE, params; only=:controls, show=loss)
    end

    plot_simulation(controlODE, θ; only=:controls, show=loss(θ))

    result = sciml_train(
        loss,
        θ,
        LBFGS(; linesearch=BackTracking());
        allow_f_increases=true,
        cb=plotting_callback,
    )

    plot_simulation(
        controlODE,
        result.minimizer;
        only=:states,
        vars=[1],
        show=loss(result.minimizer),
        yrefs=[y1s],
    )
    plot_simulation(
        controlODE,
        result.minimizer;
        only=:states,
        vars=[2],
        show=loss(result.minimizer),
        yrefs=[y2s],
    )
    plot_simulation(
        controlODE,
        result.minimizer;
        only=:controls,
        show=loss(result.minimizer),
        yrefs=[us],
    )

    metadata = Dict(
        :loss => loss(result.minimizer),
        :reaction => reaction,
        :α1 => α1,
        :α2 => α2,
        :α3 => α3,
    )
    return store_simulation(
        reaction, controlODE, result.minimizer; metadata, datadir
    )
end  # script wrapper


text:
using BHAPtfem, DataFrames, GR

ProjDir = dirname(@__FILE__)
ProjName = split(ProjDir, "/")[end]

#bha = BHA(ProjName, ProjDir)
#bha.ratio = 1.7
#bha.medium = :lightmud

segs2 = [
# Element type,  Material,    Length,     OD,         ID,        fc
  :bit,           :steel ,       0.00,    2.75,      12.25,     0.0,
  :collar ,       :monel ,      25.00,    2.75,       7.75,     0.0,
  :stabilizer,    :steel,        0.00,    2.75,      12.125,    0.0,
  :collar,        :steel,       30.00,    2.75,       7.75,     0.0,
  :stabilizer,    :steel,        0.00,    2.75,      12.125,    0.0,
  :collar,        :steel,       30.00,    2.75,       7.75,     0.0,
  :stabilizer,    :steel,        0.00,    2.75,      12.125,    0.0,
  :collar,        :steel,       90.00,    2.75,       7.75,     0.0,
  :stabilizer,    :steel,        0.00,    2.75,      12.125,    0.0,
  :pipe,          :steel,      100.00,    2.75,       7.75,     0.0
];

segs1 = reshape(segs2, (6, :));
syms = Array{Symbol, 2}(reshape(segs1[1:2, :], (2, :)))
segs = Array{Float64, 2}(reshape(segs1[3:6, :], (4, :)))

df = DataFrame([Symbol, Symbol, Float64], [:eltype, :material, :length], 0)
append!(df, DataFrame(eltype=:bit, material=:steel, length=0.0))

using GR
x = 0:pi/10:2*pi
y1 = sin.(x)
y2 = sin.(x .- 0.25)
y3 = sin.(x .- 0.5)

plot(x,y1,"g",x,y2,"--",x,y2,"g--o",x,y3,"-",x,y3,"c*")
savefig("$(ProjDir)/gr_plot_test.pdf")

df


text:
using RvSpectMLBase
using Experimental
using Test

@testset "project_flux_common_wavelengths" begin
    len = 2000
    nspec = 100
    lgx = range(log(5000.0),stop=log(5100),length=len)
    x = exp.(lgx)
    P = 0.5*log(maximum(x)/minimum(x))
    y1 = sin.(2π.*lgx./P)
    #dx = vcat(x[2]-x[1], 0.5*(x[3:end].-x[1:end-2]), x[end]-x[end-1])
    v1 = 0.001
    yo1 = y1 .+ v1.*randn(length(y1))
    yo2 = mapreduce(z->y1 .+ v1.*randn(length(y1)),hcat,1:nspec)
    v2 = v1*ones(size(yo2))
    myo2 = Experimental.calc_mean_spectrum(yo2,v2)
    @assert sum(abs2.((myo2 .- y1)./v1)) < 3*length(y1)/sqrt(length(y1))

end


text:
a = 2;
b = 3;
c = a + b;
print(c);

max_stars_repo_name: cadamsmith/julia-programming

max_issues_repo_name: cadamsmith/julia-programming

max_forks_repo_name: cadamsmith/julia-programming


text:
# Integration tests. Load the example game files, and verify that no exceptions are thrown
# No user input is inserted, so no actual gameplay is tested. Simply that the file loads
# Running these needs a graphical system -- so xvfb if running on a headless system.
@testset "Examples" begin
    @testset "basic" begin
        g = GameZero.initgame(joinpath("..","example","BasicGame","basic.jl"), true)
        @test g!=nothing
        sleep(0.2)
        GameZero.quitSDL(g)
    end

    @testset "basic2" begin
        g = GameZero.initgame(joinpath("..","example","BasicGame","basic2.jl"), true)
        @test g!=nothing
        sleep(0.2)
        GameZero.quitSDL(g)
    end
    end


text:
## Warning: type piracy

## Prior likelihood

function MeasureTheory.logdensity(prior::Dists.Dirichlet, dist::Dists.Categorical)
    return logdensity(prior, Dists.probs(dist))
end

## Sufficient statistics with prior

function Dists.suffstats(
    ::Type{<:Dists.Categorical},
    prior::Dists.Dirichlet,
    x::AbstractArray{<:Integer},
)
    return Dists.CategoricalStats(Dists.add_categorical_counts!(prior.alpha .- 1, x))
end

function Dists.suffstats(
    ::Type{<:Dists.Categorical},
    prior::Dists.Dirichlet,
    x::AbstractArray{<:Integer},
    w::AbstractArray{Float64},
)
    return Dists.CategoricalStats(Dists.add_categorical_counts!(prior.alpha .- 1, x, w))
end


text:
type ArrayContainer <: RoaringContainer
    card::UInt
    cap::UInt
    arr::Vector{UInt16}
end

ArrayContainer(size::Integer = ARRAY_DEFAULT_SIZE) = ArrayContainer(0, size, Vector{UInt16}(size))

cardinality(x::ArrayContainer) = x.card

capacity(x::ArrayContainer) = x.cap

@inline function capacity!(x::ArrayContainer, y::Integer)
    x.cap = y
end

@inline function Base.copy(x::ArrayContainer)
    return ArrayContainer(cardinality(x), capacity(x), copy(x.arr))
end

@inline function Base.copy!(dest::ArrayContainer, source::ArrayContainer)
    cardi = cardinality(source)
    if cardinality(dest) > capacity(dest)
        grow!(dest, cardi, typemax(UInt32), false)
    end
    cardinality!(dest, cardi)
    copy!(dest.arr, 1, source.arr, 1, length(source.arr))
end

@inline function shrink!(x::ArrayContainer)
    if cardinality(x) != capacity(x)
        cp = capacity(x)
        cd = cardinality(x)
        savings = cp - cd
        capacity!(x, cd)
        resize!(x.arr, cd)
    else
        savings = UInt(0)
    end
    return savings
end

@inline function growth_capacity(c::Integer)
    if c <= 0
        return ARRAY_DEFAULT_SIZE
    elseif c < 64
        return Int32(c * 2)
    elseif c < 1024
        return Int32(floor(c * 3 / 2))
    else
        return Int32(floor(c * 5 / 4))
    end
end
@inline function growth_capacity(x::ArrayContainer)
    return growth_capacity(capacity(x))
end

function grow!(x::ArrayContainer, min::Integer, max::Integer, preserve::Bool)
    ncap = clamp(growth_capacity(x), min, max)
    ncap = ncap > (max - max / 16) ? max : ncap
    capacity!(x, ncap)
    if preserve
        resize!(x.arr, ncap)
    else
        x.arr = Vector{UInt16}(ncap)
    end
    return x
end

#=
static inline bool array_container_full(const array_container_t *array) {
    return array->cardinality == array->capacity;
}
=#
isfull(x::ArrayContainer) = cardinality(x) == capacity(x)

#=
/* Append x to the set. Assumes that the value is larger than any preceding
 * values.  */
static void array_container_append(array_container_t *arr, uint16_t pos) {
    const int32_t capacity = arr->capacity;

    if (array_container_full(arr)) {
        array_container_grow(arr, capacity + 1, INT32_MAX, true);
    }

    arr->array[arr->cardinality++] = pos;
}
=#
"""
    append!(arr::ArrayContainer, pos::UInt16)

Append `pos` to the `ArrayContainer` `arr`, this assumes that the value being
appended is larger than any value currently in the container.
"""
function append!(arr::ArrayContainer, pos::UInt16)
    cap = capacity(arr)
    if isfull(arr)
        grow!(arr, cap + 1, typemax(Int32), true)
    end
    arr.card += 1
    arr.arr[arr.card] = pos
end

#=
Add all the values in [min,max) (included) at a distance k*step from min.
The container must have a size less or equal to ARRAY_DEFAULT_MAX_SIZE
after this addition.

void array_container_add_from_range(array_container_t *arr, uint32_t min,
                                    uint32_t max, uint16_t step) {
    for (uint32_t value = min; value < max; value += step) {
        array_container_append(arr, value);
    }
}
=#
function insert!{T<:Unsigned}(container::ArrayContainer, r::UnitRange{T})
    for i in r

    end
    return container
end


text:
export expand_project, expand_file, ExpandOptions

using Pkg
using TOML

function rm_using(names::Vector{Symbol}, ex)
    for name in names
        ex = rm_using(name, ex)
    end
    return ex
end

function rm_include(paths::Vector{String}, ex)
    for path in paths
        ex = rm_include(path, ex)
    end
    return ex
end

function substitute(f, ex, new)
    @match ex begin
        GuardBy(f) => new
        Expr(head, args...) => Expr(head, map(x->substitute(f, x, new), args)...)
        _ => ex
    end
end

function rm_using(name::Symbol, ex)
    @match ex begin
        :(using $(&name)) => nothing
        :(using $(&name).$(_...)) => nothing
        Expr(head, args...) => Expr(head, map(x->rm_using(name, x), args)...)
        _ => ex
    end
end

function rm_include(path::String, ex)
    @match ex begin
        :(include($(&path))) => nothing
        Expr(head, args...) => Expr(head, map(x->rm_include(path, x), args)...)
        _ => ex
    end
end

function insert_toplevel(ex, code)
    @switch ex begin
        @case Expr(:module, bare, name, Expr(:block, stmts...))
            stmts = map(x->insert_toplevel(x, code), stmts)
            return Expr(:module, bare, name, Expr(:block, code, stmts...))
        @case _
            return ex
    end
end

function expand_macro(m::Module, ex; macronames=[])
    @switch ex begin
        @case Expr(:macrocall, name, line, xs...)
            if name in macronames
                return expand_macro(m, macroexpand(m, ex); macronames=macronames)
            else
                xs = map(x->expand_macro(m, x; macronames=macronames), xs)
                return Expr(:macrocall, name, line, xs...)
            end
        @case Expr(head, args...)
            args = map(args) do x
                expand_macro(m, x; macronames=macronames)
            end
            return Expr(head, args...)
        @case _
            return ex
    end
end

"""
    ExpandOptions(;kw...)

An option type to specify available options for [`expand_project`](@ref).
See [`expand_project`](@ref) documentation for detailed kwargs.
"""
Base.@kwdef struct ExpandOptions
    mod::Module
    macronames::Vector{Any}=[]
    project::Symbol = nameof(mod)
    project_toml::String = "Project.toml"
    postfix::Symbol = :Lite
    build_dir::String = "build"
    uuid::String
    exclude_src::Vector{String} = String[] # src/*.jl to exclude
    exclude_paths::Vector{String} = String[] # misc path to exclude
    exclude_modules::Vector{Symbol} = Symbol[] # modules to remove
    src_dont_touch::Vector{String} = String[] # src/* not touch
end

# NOTE:
# this should not be an API
# it should be part of the implement detail
# of expand_file
function _replace_include(ex, options::ExpandOptions)
    @match ex begin
        :(include($path)) => begin
            if path in options.src_dont_touch
                return ex
            else
                return _xinclude_generated(path)
            end
        end
        ::Expr => Expr(ex.head, map(x->_replace_include(x, options), ex.args)...)
        _ => ex
    end
end

function _insert_include_generated(ex)
    include_generated_def = quote
        @static if !isdefined(@__MODULE__(), :include_generated)
            function _include_generated(_path::String)
                Base.@_noinline_meta
                mod = @__MODULE__()
                path, prev = Base._include_dependency(mod, _path)
                code = read(path, String)
                tls = task_local_storage()
                tls[:SOURCE_PATH] = path
                try
                    ex = include_string(mod, "quote $code end", path)
                    mod.eval(mod.eval(ex))
                    return
                finally
                    if prev === nothing
                        delete!(tls, :SOURCE_PATH)
                    else
                        tls[:SOURCE_PATH] = prev
                    end
                end
            end
        end
    end
    return insert_toplevel(ex, include_generated_def)
end

function _insert_var_str(ex)
    var_str = quote
        @static if VERSION < v"1.3"
            macro var_str(x)
                Symbol(x)
            end
        end
    end

    return insert_toplevel(ex, var_str)
end

function _xinclude_generated(path)
    return quote
        _include_generated($path)
    end
end

function expand_file(src, dst; kw...)
    expand_file(src, dst, ExpandOptions(;kw...))
end

function parse_file(src)
    raw = read(src, String)
    ex = Meta.parse("begin $raw end")
    return prettify(ex)
end

function expand_file(src, dst, options::ExpandOptions)
    ex = parse_file(src)
    ispath(dirname(dst)) || mkpath(dirname(dst))

    @info "substitute module identifier"
    old_mod = options.project
    new_mod = Symbol(options.project, options.postfix)
    ex = subtitute(ex, old_mod=>new_mod)
    ex = expand_macro(options.mod, ex; macronames=options.macronames)
    ex = rm_include(options.exclude_src, ex)
    ex = rm_include(options.exclude_paths, ex)
    ex = rm_using(options.exclude_modules, ex)
    ex = _insert_var_str(ex)
    ex = _insert_include_generated(ex)
    ex = _replace_include(ex, options)
    ex = prettify(ex)

    open(dst, "w+") do io
        println(io, ex)
    end
end

"""
    expand_project(; kw...)

Expand a list of macros for a given Julia project.

# Kwargs

- `mod::Module`: **require** package module to expand.
- `macronames::Vector{Any}`: macro names to expand, should be either a `Symbol` or `GlobalRef`.
- `project::Symbol`: Project name, default is `nameof(mod)`.
- `project_toml::String`: `Project.toml` name, should be either `Project.toml` or `JuliaProject.toml`.
- `postfix::Symbol`: the generated project postfix, default is `:Lite`.
- `build_dir::String`: the build directory to put generated proejct into.
- `uuid::String`: **require** generated project UUID.
- `exclude_src::Vector{String}`: list of file to exclude inside `<project>/src/`.
- `exclude_paths::Vector{String}`: misc paths to exclude inside `<project>`.
- `exclude_modules::Vector{Symbol}`: modules/packages to remove.
- `src_dont_touch::Vector{String}`: list files in `<project>/src/` that should not be manipulated but direct copied.
"""
function expand_project(; kw...)
    expand_project(ExpandOptions(;kw...))
end

"""
    expand_project(options::ExpandOptions)

Expand given macros in a Julia project with given `options`.
"""
function expand_project(options::ExpandOptions)
    if Sys.iswindows()
        error("expand_project does not support windows")
    end

    pathof(options.mod) === nothing && error("not a project module")
    project_dir = dirname(dirname(pathof(options.mod)))
    @info "expanding macro"
    src_dir = joinpath(project_dir, "src")
    for (root, dirs, files) in walkdir(src_dir)
        for file in files
            src = joinpath(root, file)
            relsrc = relpath(src, src_dir)
            relsrc in options.exclude_src && continue
            if file == string(options.project, ".jl") && root == src_dir
                dst = joinpath(options.build_dir, "src",
                    string(options.project, options.postfix, ".jl"))
            else
                dst = joinpath(options.build_dir, "src", relsrc)
            end

            if relsrc in options.src_dont_touch
                cp(src, dst; force=true)
            else
                @info "expanding..." src dst
                expand_file(src, dst, options)
            end
        end
    end

    @info "copying other files..."
    for each in readdir(project_dir)
        each == "src" && continue
        each == options.project_toml && continue
        each in options.exclude_paths && continue
        src = joinpath(project_dir, each)

        if isfile(src)
            dst = joinpath(options.build_dir, each)
            _cp(src, dst, options)
        else
            # check if there are files to exclude
            for (root, dirs, files) in walkdir(src)
                for file in files
                    path = joinpath(root, file)
                    relpath(path, project_dir) in options.exclude_paths && continue
                    dst = joinpath(options.build_dir, relpath(path, project_dir))
                    _cp(path, dst, options)
                end
            end
        end
    end

    # rm packages + rename project
    @info "generating new $(options.project_toml)"
    d = TOML.parsefile(joinpath(project_dir, options.project_toml))
    d["name"] = string(options.project, options.postfix)
    d["uuid"] = options.uuid

    # NOTE:
    # we don't need to care about test deps
    for each in options.exclude_modules
        package = string(each)
        delete!(d["deps"], package)
        delete!(d["compat"], package)
    end

    open(joinpath(options.build_dir, options.project_toml), "w+") do io
        TOML.print(io, d; sorted=true, by=key -> (Pkg.Types.project_key_order(key), key))
    end
    return
end

function _cp(src, dst, options::ExpandOptions)
    @info "copying..." src dst
    raw = read(src, String)
    # replace project module name
    old = string(options.mod)
    new = string(options.mod, options.postfix)
    raw = replace(raw, old=>new)

    for each in options.exclude_paths
        # use local path
        path = normpath(relpath(each, dirname(src)))
        raw = replace(raw, "include(\"$path\")"=>"nothing")
    end
    dst_dir = dirname(dst)
    ispath(dst_dir) || mkpath(dst_dir)
    write(dst, raw)
end


text:
isdefined(Base, :__precompile__) && __precompile__(false)

module RiskAdjustedLinearizations

import Base: show, getindex
import DiffEqBase: get_tmp
using ArrayInterface, FastGaussQuadrature, FiniteDiff, ForwardDiff, LinearAlgebra, Printf
using SparseArrays, SparseDiffTools, SparsityDetection, UnPack
using BandedMatrices: Zeros
using DiffEqBase: DiffCache, dualcache
using NLsolve: nlsolve, OnceDifferentiable

# Utilities
include("util.jl") # in principle, several of these utility functions could work for other packages

# RiskAdjustedLinearization
include("cache_types.jl")
include("risk_adjusted_linearization.jl")
include("sparse_jacobian_helpers.jl")

# Numerical Algorithms
include("numerical_algorithms/compute_psi.jl")
include("numerical_algorithms/blanchard_kahn.jl")
include("numerical_algorithms/relaxation.jl")
include("numerical_algorithms/homotopy.jl")
include("numerical_algorithms/solve.jl")

# Simulations and Simulation-based Diagnostics
include("simulation/simulate.jl")
include("simulation/impulse_responses.jl")
include("simulation/gausshermite_expectation.jl")
include("simulation/euler_equation_error.jl")

export
    # main functions
    RiskAdjustedLinearization, update!, nonlinear_system, linearized_system, solve!,
    steady_state_errors,

    # simulation related functions
    simulate, impulse_responses, gausshermite_expectation, euler_equation_error,
    dynamic_euler_equation_error,

    # sparsity helpers
    preallocate_jac_cache, compute_sparsity_pattern, update_sparsity_pattern!
end # module


text:
module Simkl

using HTTP
using JSON
using Dates: now
using Base64

const redirect_uri = Ref("urn:ietf:wg:oauth:2.0:oob")
const simkl_oauth = "https://simkl.com/oauth/authorize"
const simkl_pin = "https://api.simkl.com/oauth/pin"
const simkl_all_items = "https://api.simkl.com/sync/all-items"

const cache_path = Ref(get(ENV, "XDG_CACHE_HOME", "$(ENV["HOME"])/.cache"))
const creds_path = Ref(joinpath(cache_path[], "simkl_creds.json"))
const items_path = Ref(joinpath(cache_path[], "simkl_items.json"))
const creds = IdDict{String, String}()
merge!(creds, JSON.parse(read(creds_path[], String)))
const access_token = Ref(get(creds, "access_token", ""))
const headers = []

@inline function tv_url(ids)
    if "imdb" ∈ keys(ids)
        "https://www.imdb.com/title/" * ids["imdb"]
    elseif "tvdbslug" ∈ keys(ids)
        "https://www.thetvdb.com/series/" * ids["tvdbslug"]
    elseif "anidb" ∈ keys(ids)
        "https://anidb.net/anime/" * ids["anidb"]
    else
        ""
    end
end

@inline function simkl_url(type, show::String)
    "https://simkl.com/" * type * "/" * show
end

function get_simkl_pin()
    query = Dict("client_id" => creds["client_id"],
                 "redirect_uri" => redirect_uri[])
    res = HTTP.request("GET", simkl_pin, headers; query)
    body = JSON.parse(String(res.body))
    body["user_code"], body["verification_url"]
end

function get_simkl_token(code)
    query = Dict("client_id" => creds["client_id"])
    res = HTTP.request("GET", joinpath(simkl_pin, code), headers; query)
    body = JSON.parse(String(res.body))
    get(body, "access_token", "")
end

function set_headers!()
    empty!(headers)
    push!(headers, "Content-Type" => "application/json")
    push!(headers, "Authorization" => "Bearer $(access_token[])")
    push!(headers, "simkl-api-key" => creds["client_id"])
end

@doc "Checks if an ACCESS_TOKEN key is present in the file pointed by CREDS_PATH. If not present initiate
a pin verification procedure."
function simkl_auth()
    token = ""
    if isempty(get(creds, "access_token", ""))
        code, url = get_simkl_pin()
        display("Verify pin: $code at $url")
        sl = 1
        while true
            token = get_simkl_token(code)
            isempty(token) || break
            sleep(sl)
            sl += 1
        end
        creds["code"] = code
        creds["access_token"] = token
        access_token[] = token
        write(creds_path[], JSON.json(creds))
        display("Saved new access token to $(creds_path[])")
    else
        display("Access token already available.")
    end
    set_headers!()
end

function simkl_fetch_all_items(type="", status=""; reset=false)
    date_from = reset ? "" : get(creds, "date_from", "")
    query = Dict()
    prev_items = nothing
    isempty(date_from) || begin
	    query["date_from"] = date_from
        isfile(items_path[]) && begin
            prev_items = JSON.parse(read(items_path[], String))
        end
    end
    set_headers!()
    res = HTTP.request("GET", joinpath(simkl_all_items, type, status), headers; query)
    if isnothing(prev_items)
        items_str = String(res.body)
        write(items_path[], items_str)
        prev_items = JSON.parse(items_str)
    else
        items = JSON.parse(String(res.body))
        if !isnothing(items)
            merge!(prev_items, items)
            @show prev_items
            write(items_path[], JSON.json(prev_items))
        end
    end
    creds["date_from"] = string(now())
    write(creds_path[], JSON.json(creds))
    prev_items
end

function simkl_get_all_items(update=false; reset=false, kwargs...)
    if update || !isfile(items_path[])
        simkl_fetch_all_items(;reset, kwargs...)
    else
	    JSON.parse(read(items_path[], String))
    end
end

@doc "Generates an HTML list of completed shows with links to simkl and imdb."
function simkl_completed_shows_list(types=["shows", "anime"], status=["completed", "watching"])
    io = IOBuffer()
	all_items = simkl_get_all_items()
    write(io, "
")
    for tp in types
        el_tp = tp === "movies" ? "movie" : "show"
        type_items = all_items[tp]
        if length(type_items) > 0
            write(io, "
")
            write(io, titlecase(tp))
            write(io, "

")
        for itm in type_items
            if itm["status"] ∈ status
                write(io, "
")
                ids = itm[el_tp]["ids"]
                write(io, "")
                write(io, itm[el_tp]["title"])
                write(io, "")
                write(io, "")
                write(io, "")
                write(io, "")
                write(io, "
")
            end
        end
            write(io, "
")
        end
    end
    write(io, "
")
    list = String(take!(io))
    close(io)
    list
end

function hfun_simkl_list()
	simkl_completed_shows_list()
end

function simkl_get_shows(status="watching"; types = ["shows", "anime"])
    all_items = simkl_get_all_items()
    shows = []
    for tp in types
        for el in all_items[tp]
            if el["status"] === status
                push!(shows, el)
            end
        end
    end
    shows
end

function simkl_get_show_by_title(title; status="watching")
    shows = simkl_get_shows(status)
    for s in shows
        s["show"]["title"] === title && return s
    end
end

export hfun_simkl_list

function get_show_id(show)
    ids = show["ids"]
    k = keys(show["ids"])
    # medusa doesn't want imdb as indexer
    # "imdb" ∈ k && return "imdb" => ids["imdb"]
    "tvdb" ∈ k && return "tvdb" => ids["tvdb"]
    "tmdb" ∈ k && return "tmdb" => ids["tmdb"]
    "anidb" ∈ k && return "imdb" => ids["anidb"]
end

end


text:
using WeatherStation
using Dates

basedir = expanduser("~/.julia/dev/WeatherStation/examples")
parameter = "temperature"
taverage = 1000 * 60 * 10

#=
curl 'http://192.168.1.61:8081/_dash-update-component' -H 'User-Agent: Mozilla/5.0 (X11; Ubuntu; Linux x86_64; rv:76.0) Gecko/20100101 Firefox/76.0' -H 'Accept: application/json' -H 'Accept-Language: en-US,de;q=0.7,en;q=0.3' --compressed -H 'Referer: http://192.168.1.61:8081/' -H 'Content-Type: application/json' -H 'X-CSRFToken: undefined' -H 'Origin: http://192.168.1.61:8081' -H 'DNT: 1' -H 'Connection: keep-alive' --data-raw '{"output":"graph.figure","changedPropIds":["date-picker-range.end_date"],"inputs":[{"id":"parameter","property":"value","value":"temperature"},{"id":"taverage","property":"value","value":600000},{"id":"date-picker-range","property":"start_date","value":"2020-07-01"},{"id":"date-picker-range","property":"end_date","value":"2020-08-04"}]}'
curl 'http://192.168.1.61:8081/_dash-update-component' -H 'User-Agent: Mozilla/5.0 (X11; Ubuntu; Linux x86_64; rv:76.0) Gecko/20100101 Firefox/76.0' -H 'Accept: application/json' -H 'Accept-Language: en-US,de;q=0.7,en;q=0.3' --compressed -H 'Referer: http://192.168.1.61:8081/' -H 'Content-Type: application/json' -H 'X-CSRFToken: undefined' -H 'Origin: http://192.168.1.61:8081' -H 'DNT: 1' -H 'Connection: keep-alive' --data-raw '{"output":"graph.figure","changedPropIds":["date-picker-range.end_date"],"inputs":[{"id":"parameter","property":"value","value":"temperature"},{"id":"taverage","property":"value","value":600000},{"id":"date-picker-range","property":"start_date","value":"2020-07-01"},{"id":"date-picker-range","property":"end_date","value":"2020-08-04"}]}'
=#

timerange = (typemin(DateTime),typemax(DateTime))

timerange = (DateTime(2020,07,06),DateTime(2020,07,08))

WeatherStation.loadavg(basedir,parameter,taverage;
	timerange = timerange)


text:
using StaticArrays, CUDAnative

function Base.sort!(a::MArray{Tuple{3}})
  # Use a (Bose-Nelson Algorithm based) sorting network from
  # <http://pages.ripco.net/~jgamble/nw.html>.
  a[2], a[3] = minmax(a[2], a[3])
  a[1], a[3] = minmax(a[1], a[3])
  a[1], a[2] = minmax(a[1], a[2])
end

function knl!()
  a = @MVector [7, 3, 4]
  sort!(a)
  @cuprintf("(%ld, %ld, %ld)\n", a[1], a[2], a[3])
  nothing
end

@cuda knl!()


text:
macro F_YF() #(V, U, Ax, R, SA)		# front tire force
  code = quote
    if length(Ax) != 1
      Fyf = zeros(length(Ax),1)
      for ii in eachindex(Ax)
        Fyf[ii] = 278.2789450759981 + ((PD2*(FzF0 - (Ax[ii] - V[ii]*R[ii])*KZX)^2 + PD1*(FzF0 - (Ax[ii] - V[ii]*R[ii])*KZX))*sin(PC1*atan((((PK1*sin(2*atan(PK2*(FzF0 - (Ax[ii] - V[ii]*R[ii])*KZX))))/(((PD2*PC1*(FzF0 - (Ax[ii] - V[ii]*R[ii])*KZX)^2 + PD1*PC1*(FzF0 - (Ax[ii] - V[ii]*R[ii])*KZX)) + ((PD2*PC1*(FzF0 - (Ax[ii] - V[ii]*R[ii])*KZX)^2 + PD1*PC1*(FzF0 - (Ax[ii] - V[ii]*R[ii])*KZX)))/(((PD2*PC1*(FzF0 - (Ax[ii] - V[ii]*R[ii])*KZX)^2 + PD1*PC1*(FzF0 - (Ax[ii] - V[ii]*R[ii])*KZX))^2 )^(0.5)))+EP))*((atan((V[ii] + la*R[ii])/(U[ii]+EP)) - SA[ii]) + PH2*(FzF0 - (Ax[ii] - V[ii]*R[ii])*KZX) + PH1)) - ((PE2*(FzF0 - (Ax[ii] - V[ii]*R[ii])*KZX) + PE1)*(1 - PE3)*(((atan((V[ii] + la*R[ii])/(U[ii]+EP)) - SA[ii]) + PH2*(FzF0 - (Ax[ii] - V[ii]*R[ii])*KZX) + PH1))/((((atan((V[ii] + la*R[ii])/(U[ii]+EP)) - SA[ii]) + PH2*(FzF0 - (Ax[ii] - V[ii]*R[ii])*KZX) + PH1)^2 )^(0.5)))*((((PK1*sin(2*atan(PK2*(FzF0 - (Ax[ii] - V[ii]*R[ii])*KZX))))/(((PD2*PC1*(FzF0 - (Ax[ii] - V[ii]*R[ii])*KZX)^2 + PD1*PC1*(FzF0 - (Ax[ii] - V[ii]*R[ii])*KZX)) + ((PD2*PC1*(FzF0 - (Ax[ii] - V[ii]*R[ii])*KZX)^2 + PD1*PC1*(FzF0 - (Ax[ii] - V[ii]*R[ii])*KZX)))/(((PD2*PC1*(FzF0 - (Ax[ii] - V[ii]*R[ii])*KZX)^2 + PD1*PC1*(FzF0 - (Ax[ii] - V[ii]*R[ii])*KZX))^2 )^(0.5)))+EP))*((atan((V[ii] + la*R[ii])/(U[ii]+EP)) - SA[ii]) + PH2*(FzF0 - (Ax[ii] - V[ii]*R[ii])*KZX) + PH1)) - atan((((PK1*sin(2*atan(PK2*(FzF0 - (Ax[ii] - V[ii]*R[ii])*KZX))))/(((PD2*PC1*(FzF0 - (Ax[ii] - V[ii]*R[ii])*KZX)^2 + PD1*PC1*(FzF0 - (Ax[ii] - V[ii]*R[ii])*KZX)) + ((PD2*PC1*(FzF0 - (Ax[ii] - V[ii]*R[ii])*KZX)^2 + PD1*PC1*(FzF0 - (Ax[ii] - V[ii]*R[ii])*KZX)))/(((PD2*PC1*(FzF0 - (Ax[ii] - V[ii]*R[ii])*KZX)^2 + PD1*PC1*(FzF0 - (Ax[ii] - V[ii]*R[ii])*KZX))^2 )^(0.5)))+EP))*((atan((V[ii] + la*R[ii])/(U[ii]+EP)) - SA[ii]) + PH2*(FzF0 - (Ax[ii] - V[ii]*R[ii])*KZX) + PH1)))))) + (PV2*(FzF0 - (Ax[ii] - V[ii]*R[ii])*KZX)^2 + PV1*(FzF0 - (Ax[ii] - V[ii]*R[ii])*KZX)));
      end
    else
      Fyf = Float64
      Fyf = 278.2789450759981 + ((PD2*(FzF0 - (Ax - V*R)*KZX)^2 + PD1*(FzF0 - (Ax - V*R)*KZX))*sin(PC1*atan((((PK1*sin(2*atan(PK2*(FzF0 - (Ax - V*R)*KZX))))/(((PD2*PC1*(FzF0 - (Ax - V*R)*KZX)^2 + PD1*PC1*(FzF0 - (Ax - V*R)*KZX)) + ((PD2*PC1*(FzF0 - (Ax - V*R)*KZX)^2 + PD1*PC1*(FzF0 - (Ax - V*R)*KZX)))/(((PD2*PC1*(FzF0 - (Ax - V*R)*KZX)^2 + PD1*PC1*(FzF0 - (Ax - V*R)*KZX))^2 )^(0.5)))+EP))*((atan((V + la*R)/(U+EP)) - SA) + PH2*(FzF0 - (Ax - V*R)*KZX) + PH1)) - ((PE2*(FzF0 - (Ax - V*R)*KZX) + PE1)*(1 - PE3)*(((atan((V + la*R)/(U+EP)) - SA) + PH2*(FzF0 - (Ax - V*R)*KZX) + PH1))/((((atan((V + la*R)/(U+EP)) - SA) + PH2*(FzF0 - (Ax - V*R)*KZX) + PH1)^2 )^(0.5)))*((((PK1*sin(2*atan(PK2*(FzF0 - (Ax - V*R)*KZX))))/(((PD2*PC1*(FzF0 - (Ax - V*R)*KZX)^2 + PD1*PC1*(FzF0 - (Ax - V*R)*KZX)) + ((PD2*PC1*(FzF0 - (Ax - V*R)*KZX)^2 + PD1*PC1*(FzF0 - (Ax - V*R)*KZX)))/(((PD2*PC1*(FzF0 - (Ax - V*R)*KZX)^2 + PD1*PC1*(FzF0 - (Ax - V*R)*KZX))^2 )^(0.5)))+EP))*((atan((V + la*R)/(U+EP)) - SA) + PH2*(FzF0 - (Ax - V*R)*KZX) + PH1)) - atan((((PK1*sin(2*atan(PK2*(FzF0 - (Ax - V*R)*KZX))))/(((PD2*PC1*(FzF0 - (Ax - V*R)*KZX)^2 + PD1*PC1*(FzF0 - (Ax - V*R)*KZX)) + ((PD2*PC1*(FzF0 - (Ax - V*R)*KZX)^2 + PD1*PC1*(FzF0 - (Ax - V*R)*KZX)))/(((PD2*PC1*(FzF0 - (Ax - V*R)*KZX)^2 + PD1*PC1*(FzF0 - (Ax - V*R)*KZX))^2 )^(0.5)))+EP))*((atan((V + la*R)/(U+EP)) - SA) + PH2*(FzF0 - (Ax - V*R)*KZX) + PH1)))))) + (PV2*(FzF0 - (Ax - V*R)*KZX)^2 + PV1*(FzF0 - (Ax - V*R)*KZX)));
    end
    Fyf
  end
  return esc(code)
end


text:
import Base.convert
import Base.string
import Base.print
import Base.println

# Bitcoin's variable length integer
# https://en.bitcoin.it/wiki/Protocol_specification#Variable_length_integer
#
# Only providing convenience functions to and from Array{Uint8} because varint
# seems like its primarily a transport structure
#
# VarInt are created as big-endian here, but are sent little-endian as per
# https://bitcoin.org/en/developer-reference#compactsize-unsigned-integers
function to_varint(x::Integer)
  result = Array(Uint8, 0)

  if x < 0
    error("Negative values for VarInt undefined.")
  elseif x < 0xfd
    append!(result, Crypto.int2oct(uint8(x)))
  elseif x <= 0xffff
    append!(result, [0xfd])
    append!(result, Crypto.int2oct(uint16(x)))
  elseif x <= 0xffffffff
    append!(result, [0xfe])
    append!(result, Crypto.int2oct(uint32(x)))
  else
    append!(result, [0xff])
    append!(result, Crypto.int2oct(uint64(x)))
  end

  return result
end

# Bitcoin's variable length string
# https://en.bitcoin.it/wiki/Protocol_specification#Variable_length_string
function to_varstring(x::String)
  return [to_varint(length(x)), x.data]
end

function reverse_endian(hex_string::String)
  return join([hex(x, 2) for x in reverse(hex2oct(hex_string))])
end

function reverse_endian(hex_data::Integer)
  data_type = typeof(hex_data)
  num_bytes = sizeof(data_type)
  counter = 0
  result = 0

  # Assumes BigInts are at least 16 bytes
  while counter < num_bytes || hex_data > 0
    result = convert(data_type, result << 8 + hex_data & 0xff)
    hex_data >>>= 8
    counter += 1
  end

  return result
end

function get_checksum(message::String; is_hex = false)
  message = Crypto.digest("SHA256", Crypto.digest("SHA256", message, is_hex = true))
  return Crypto.oct2hex(message[1:4])
end

function get_checksum(message::Array{Uint8})
  return Crypto.digest("SHA256", Crypto.digest("SHA256", message))[1:4]
end

# Convenience function for converting to Array{Uint8}
# Must be defined for given types
function bytearray(x)
  convert(Array{Uint8}, x)
end

function string(array::Array{Uint8})
  return string("[", join([string("0x", hex(x, 2)) for x in array], ", "), "]")
end

print(x::Array{Uint8}) = print(string(x))
println(x::Array{Uint8}) = print(x); print("\n")

function string(x::Unsigned)
  return string("0x", hex(x, div(sizeof(x), 2)))
end

print(x::Unsigned) = print(string(x))
println(x::Unsigned) = print(x); print("\n")