text:
## Agents consume, agents save, agents get expenditure shocks
## If agents default they are excluded from financial markets

# z ... income shock 1
# η ... income shock 2
# κ ... expenditure shock

using QuantEcon, Parameters, Optim, Interpolations
include("../src/ExogenousStates.jl")
using .ExogenousStates
using Aiyagari

## Incomes z, avg income = 1
z_grid = [0.5; 1.0; 1.5]
z_prob = [0.7 0.15 0.15;
          0.15 0.7 0.15;
          0.15 0.15 0.7]
z_MC = MarkovChain(z_prob, z_grid)

## Expenditure shock
κ_grid = [0.0, 0.26, 0.82]
κ_prob = [0.925 0.07 0.005;
          0.925 0.07 0.005;
          0.925 0.07 0.005]
κ_MC = MarkovChain(κ_prob, κ_grid)

# n_exo_states = n_states(exo_MC)
# n_endo_states = length(a_grid)
#
# V_default = zeros(n_endo_states, n_exo_states)
# V_no_default = zeros(n_endo_states, n_exo_states)

# expectation

u(c) = c > 0 ? log(c) : c * 1e5 - 1e3

q(a,z,agg_state) = (1-agg_state.θ(a,z))/(1+agg_state.r)

function Aiyagari.iterate_bellman!(value_new, value_old, policy, policies_full, a_grid, z_mc, converged, agg_state)
  β, γ = 0.95, 0.5
  a_min, a_max = extrema(a_grid)

  for (i_z, nt) in enumerate(z_mc.state_values)
    @unpack z, κ = nt
    # Create interpolated expected value function
    exp_value = value_old * z_mc.p[i_z,:]

    itp_exp_value = interpolate(exp_value, BSpline(Cubic(Line(OnGrid()))))
    sitp_exp_value = scale(itp_exp_value, a_grid)

    for (i_a, a) in enumerate(a_grid)

      obj = a_next -> begin
        c = a + z - κ - q(a_next,z,agg_state) * a_next
        u(c) + β * sitp_exp_value(a_next)
      end

      res = optimize(a_next -> - obj(a_next), a_min, a_max)

      a_no_def = Optim.minimizer(res)
      c_no_def = a + z - κ - a_no_def * q(a_no_def,z,agg_state)
      v_no_def = - Optim.minimum(res)

      # default
      c_default = (1-γ) * z

      a_default = (a + γ * z - κ) * (1+agg_state.r)
      v_default = u(c_default) + β * extrapolate(sitp_exp_value, Interpolations.Line())(a_default)

      if v_no_def > v_default
        policy[i_a, i_z] = a_no_def
        policies_full[i_a, i_z]= (default = false, c = c_no_def)
        value_new[i_a, i_z]    = v_no_def
        converged[i_a, i_z]    = Optim.converged(res)
      else
        policy[i_a, i_z] = a_default
        policies_full[i_a, i_z]= (default = true, c = c_default)
        value_new[i_a, i_z]    = v_default
        converged[i_a, i_z]    = Optim.converged(res)
      end
    end
  end
end

using StructArrays
using StatsBase: Weights
using DataFrames

mutable struct AggStateNoFreshStart{T1, T2, T3}
  r::T1
  default_choice::T2
  θ::T3 # default_probability
end

function AggStateNoFreshStart(r, a_grid::AbstractVector, exo_mc::MarkovChain)
  size = (length(a_grid), length(exo_mc.state_values))

  default_choice = falses(size)

  AggStateNoFreshStart(r, default_choice, default_probabilities(a_grid, exo_mc))
end

function update_probabilities!(agg_state, a_grid, exo_mc, policies_full, dist)
  agg_state.default_choice .=  default_choice(policies_full)
  agg_state.θ = default_probabilities(a_grid, exo_mc, agg_state.default_choice, dist)
end

function default_choice(policies_full)
  policies_SoA = StructArray(policies_full)
  default_choice = policies_SoA.default
end

function default_probabilities(a_grid, exo_mc, default_choice, dist)

  statespace = DataFrame([(a=a, exo...) for exo in exo_MC.state_values for a in a_grid])

  statespace[!,:pmf] = vec(dist)
  statespace[!,:default] = vec(default_choice)

  default_df = by(statespace, [:a, :z]) do df
    (default = mean(df.default, Weights(df.pmf)),)
  end

  default_df.default

  z_grid = unique(default_df.z)
  out = reshape(default_df.default, (length(a_grid), length(z_grid)))

  itp = interpolate((a_grid, z_grid), out, (Gridded(Linear()), Gridded(Constant())))
end

function default_probabilities(a_grid, exo_mc)
  size = (length(a_grid), length(exo_mc.state_values))

  default_probabilities(a_grid, exo_mc, falses(size), ones(size)/prod(size))
end

exo_MC = MarkovChain(z_MC, κ_MC, (:z, :κ))
exo_MC.state_values[1]

## Endogenous state
a_grid = LinRange(-10.0, 20.0, 300)

## Finding an equilibrium

# step 0: θ ≡ 0, r = r₀

using LinearAlgebra: norm, dot

function iterate_defaults(a_grid, exo_MC, agg_state; max_iter_Bellman=400) # solve for 0 profits of intermediary
  # Solve given "interest schedule"
  @unpack value, policy, policies_full = solve_bellman(a_grid, exo_MC, (default=true, c=1.5), agg_state, maxiter=max_iter_Bellman )
  dist = stationary_distribution(exo_MC, a_grid, policy)

  # Update interest schedule
  update_probabilities!(agg_state, a_grid, exo_MC, policies_full, dist)
  dist
end

function solve_defaults(a_grid, exo_MC, agg_state; max_iter_Bellman=400, max_iter=30)
  def_old = copy(agg_state.default_choice)
  θ_old = [agg_state.θ(a,z) for (a,z) in Iterators.product(a_grid, z_grid)]

  for i in 1:max_iter
    dist = iterate_defaults(a_grid, exo_MC, agg_state, max_iter_Bellman=max_iter_Bellman)

    def_new = agg_state.default_choice
    θ_new = [agg_state.θ(a,z) for (a,z) in Iterators.product(a_grid, z_grid)]

    #norm(θ_new - θ_old, Inf)

    def_switch = sum(def_new .!= def_old)

    mass = sum(dist[def_new .!= def_old])
    if def_switch/length(def_new) < 0.001 || mass < 1e-5
      break
    end
    if i%1 == 0
      @info "(#: $def_switch %: $(def_switch/length(def_new)), mass = $mass)"
    end

    if i == max_iter
      @warn "Default choices didn't converge (#: $def_switch %: $(def_switch/length(def_new)), mass = $mass)"
    end
    def_old .= def_new
    θ_old .= θ_new
  end
end

agg_state = AggStateNoFreshStart(0.05, a_grid, exo_MC)

solve_defaults(a_grid, exo_MC, agg_state)

agg_state.default_choice

begin
  plt_θ = plot(a_grid, def_prob)
  plt_dist = plot(a_grid, dist)
  display(plot(plt_θ, plt_dist))
end

## profits are zero automatically
## asset markets still need to clear
dist = stationary_distribution(exo_MC, a_grid, policy)

excess_demand(agg_state, a_grid, dist) = dot(a_grid, sum(dist, dims=2))

excess_demand(agg_state, a_grid, dist)

ζ(r, )

using StatsPlots

policies_SoA = StructArray(policies_full)

plot(a_grid, policies_SoA.default)

plot(agg_state.default_choice)

plot(a_grid, value)
plot(a_grid, policy)

@df default_df plot(:a, :default, group=:z)


text:
abstract type RosenbrockMutableCache <: OrdinaryDiffEqMutableCache end
################################################################################

# Shampine's Low-order Rosenbrocks

@cache mutable struct Rosenbrock23Cache{uType,rateType,JType,WType,TabType,TFType,UFType,F,JCType,GCType} <: RosenbrockMutableCache
  u::uType
  uprev::uType
  k₁::rateType
  k₂::rateType
  k₃::rateType
  du1::rateType
  du2::rateType
  f₁::rateType
  fsalfirst::rateType
  fsallast::rateType
  dT::rateType
  J::JType
  W::WType
  tmp::rateType
  tab::TabType
  tf::TFType
  uf::UFType
  linsolve_tmp::rateType
  linsolve::F
  jac_config::JCType
  grad_config::GCType
end

@cache mutable struct Rosenbrock32Cache{uType,rateType,JType,WType,TabType,TFType,UFType,F,JCType,GCType} <: RosenbrockMutableCache
  u::uType
  uprev::uType
  k₁::rateType
  k₂::rateType
  k₃::rateType
  du1::rateType
  du2::rateType
  f₁::rateType
  fsalfirst::rateType
  fsallast::rateType
  dT::rateType
  J::JType
  W::WType
  tmp::rateType
  tab::TabType
  tf::TFType
  uf::UFType
  linsolve_tmp::rateType
  linsolve::F
  jac_config::JCType
  grad_config::GCType
end

function alg_cache(alg::Rosenbrock23,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  k₁ = zero(rate_prototype)
  k₂ = zero(rate_prototype)
  k₃ = zero(rate_prototype)
  du1 = zero(rate_prototype)
  du2 = zero(rate_prototype)
  # f₀ = similar(u) fsalfirst
  f₁ = zero(rate_prototype)
  fsalfirst = zero(rate_prototype)
  fsallast = zero(rate_prototype)
  dT = zero(rate_prototype)
  if DiffEqBase.has_jac(f) && !DiffEqBase.has_Wfact(f) && f.jac_prototype !== nothing
    W = WOperator(f, dt, true)
    J = nothing # is J = W.J better?
  else
    J = false .* vec(rate_prototype) .* vec(rate_prototype)' # uEltype?
    W = similar(J)
  end
  tmp = zero(rate_prototype)
  tab = Rosenbrock23ConstantCache(real(uBottomEltypeNoUnits),identity,identity)
  tf = DiffEqDiffTools.TimeGradientWrapper(f,uprev,p)
  uf = DiffEqDiffTools.UJacobianWrapper(f,t,p)
  linsolve_tmp = zero(rate_prototype)
  linsolve = alg.linsolve(Val{:init},uf,u)

  grad_config = build_grad_config(alg,f,tf,du1,t)
  jac_config = build_jac_config(alg,f,uf,du1,uprev,u,tmp,du2)

  Rosenbrock23Cache(u,uprev,k₁,k₂,k₃,du1,du2,f₁,
                    fsalfirst,fsallast,dT,J,W,tmp,tab,tf,uf,linsolve_tmp,
                    linsolve,jac_config,grad_config)
end

function alg_cache(alg::Rosenbrock32,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  k₁ = zero(rate_prototype)
  k₂ = zero(rate_prototype)
  k₃ = zero(rate_prototype)
  du1 = zero(rate_prototype)
  du2 = zero(rate_prototype)
  # f₀ = similar(u) fsalfirst
  f₁ = zero(rate_prototype)
  fsalfirst = zero(rate_prototype)
  fsallast = zero(rate_prototype)
  dT = zero(rate_prototype)
  if DiffEqBase.has_jac(f) && !DiffEqBase.has_Wfact(f) && f.jac_prototype !== nothing
    W = WOperator(f, dt, true)
    J = nothing # is J = W.J better?
  else
    J = false .* vec(rate_prototype) .* vec(rate_prototype)' # uEltype?
    W = similar(J)
  end
  tmp = zero(rate_prototype)
  tab = Rosenbrock32ConstantCache(real(uBottomEltypeNoUnits),identity,identity)

  tf = DiffEqDiffTools.TimeGradientWrapper(f,uprev,p)
  uf = DiffEqDiffTools.UJacobianWrapper(f,t,p)
  linsolve_tmp = zero(rate_prototype)
  linsolve = alg.linsolve(Val{:init},uf,u)
  grad_config = build_grad_config(alg,f,tf,du1,t)
  jac_config = build_jac_config(alg,f,uf,du1,uprev,u,tmp,du2)
  Rosenbrock32Cache(u,uprev,k₁,k₂,k₃,du1,du2,f₁,fsalfirst,fsallast,dT,J,W,tmp,tab,tf,uf,linsolve_tmp,linsolve,jac_config,grad_config)
end

struct Rosenbrock23ConstantCache{T,TF,UF} <: OrdinaryDiffEqConstantCache
  c₃₂::T
  d::T
  tf::TF
  uf::UF
end

function Rosenbrock23ConstantCache(T::Type,tf,uf)
  c₃₂ = convert(T,6 + sqrt(2))
  d = convert(T,1/(2+sqrt(2)))
  Rosenbrock23ConstantCache(c₃₂,d,tf,uf)
end

function alg_cache(alg::Rosenbrock23,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}})
  tf = DiffEqDiffTools.TimeDerivativeWrapper(f,u,p)
  uf = DiffEqDiffTools.UDerivativeWrapper(f,t,p)
  Rosenbrock23ConstantCache(real(uBottomEltypeNoUnits),tf,uf)
end

struct Rosenbrock32ConstantCache{T,TF,UF} <: OrdinaryDiffEqConstantCache
  c₃₂::T
  d::T
  tf::TF
  uf::UF
end

function Rosenbrock32ConstantCache(T::Type,tf,uf)
  c₃₂ = convert(T,6 + sqrt(2))
  d = convert(T,1/(2+sqrt(2)))
  Rosenbrock32ConstantCache(c₃₂,d,tf,uf)
end

function alg_cache(alg::Rosenbrock32,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}})
  tf = DiffEqDiffTools.TimeDerivativeWrapper(f,u,p)
  uf = DiffEqDiffTools.UDerivativeWrapper(f,t,p)
  Rosenbrock32ConstantCache(real(uBottomEltypeNoUnits),tf,uf)
end

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

### 3rd order specialized Rosenbrocks

struct Rosenbrock33ConstantCache{TF,UF,Tab} <: OrdinaryDiffEqConstantCache
  tf::TF
  uf::UF
  tab::Tab
end

@cache mutable struct Rosenbrock33Cache{uType,rateType,JType,WType,TabType,TFType,UFType,F,JCType,GCType} <: RosenbrockMutableCache
  u::uType
  uprev::uType
  du::rateType
  du1::rateType
  du2::rateType
  k1::rateType
  k2::rateType
  k3::rateType
  k4::rateType
  fsalfirst::rateType
  fsallast::rateType
  dT::rateType
  J::JType
  W::WType
  tmp::rateType
  tab::TabType
  tf::TFType
  uf::UFType
  linsolve_tmp::rateType
  linsolve::F
  jac_config::JCType
  grad_config::GCType
end

function alg_cache(alg::ROS3P,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  du = zero(rate_prototype)
  du1 = zero(rate_prototype)
  du2 = zero(rate_prototype)
  k1 = zero(rate_prototype)
  k2 = zero(rate_prototype)
  k3 = zero(rate_prototype)
  k4 = zero(rate_prototype)
  fsalfirst = zero(rate_prototype)
  fsallast = zero(rate_prototype)
  dT = zero(rate_prototype)
  if DiffEqBase.has_jac(f) && !DiffEqBase.has_Wfact(f) && f.jac_prototype !== nothing
    W = WOperator(f, dt, true)
    J = nothing # is J = W.J better?
  else
    J = false .* vec(rate_prototype) .* vec(rate_prototype)' # uEltype?
    W = similar(J)
  end
  tmp = zero(rate_prototype)
  tab = ROS3PConstantCache(constvalue(uBottomEltypeNoUnits),constvalue(tTypeNoUnits))
  tf = DiffEqDiffTools.TimeGradientWrapper(f,uprev,p)
  uf = DiffEqDiffTools.UJacobianWrapper(f,t,p)
  linsolve_tmp = zero(rate_prototype)
  linsolve = alg.linsolve(Val{:init},uf,u)
  grad_config = build_grad_config(alg,f,tf,du1,t)
  jac_config = build_jac_config(alg,f,uf,du1,uprev,u,tmp,du2)
  Rosenbrock33Cache(u,uprev,du,du1,du2,k1,k2,k3,k4,
                    fsalfirst,fsallast,dT,J,W,tmp,tab,tf,uf,linsolve_tmp,
                    linsolve,jac_config,grad_config)
end

function alg_cache(alg::ROS3P,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}})
  tf = DiffEqDiffTools.TimeDerivativeWrapper(f,u,p)
  uf = DiffEqDiffTools.UDerivativeWrapper(f,t,p)
  Rosenbrock33ConstantCache(tf,uf,ROS3PConstantCache(constvalue(uBottomEltypeNoUnits),constvalue(tTypeNoUnits)))
end

@cache mutable struct Rosenbrock34Cache{uType,rateType,JType,WType,TabType,TFType,UFType,F,JCType,GCType} <: RosenbrockMutableCache
  u::uType
  uprev::uType
  du::rateType
  du1::rateType
  du2::rateType
  k1::rateType
  k2::rateType
  k3::rateType
  k4::rateType
  fsalfirst::rateType
  fsallast::rateType
  dT::rateType
  J::JType
  W::WType
  tmp::rateType
  tab::TabType
  tf::TFType
  uf::UFType
  linsolve_tmp::rateType
  linsolve::F
  jac_config::JCType
  grad_config::GCType
end

function alg_cache(alg::Rodas3,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  du = zero(rate_prototype)
  du1 = zero(rate_prototype)
  du2 = zero(rate_prototype)
  k1 = zero(rate_prototype)
  k2 = zero(rate_prototype)
  k3 = zero(rate_prototype)
  k4 = zero(rate_prototype)
  fsalfirst = zero(rate_prototype)
  fsallast = zero(rate_prototype)
  dT = zero(rate_prototype)
  if DiffEqBase.has_jac(f) && !DiffEqBase.has_Wfact(f) && f.jac_prototype !== nothing
    W = WOperator(f, dt, true)
    J = nothing # is J = W.J better?
  else
    J = false .* vec(rate_prototype) .* vec(rate_prototype)' # uEltype?
    W = similar(J)
  end
  tmp = zero(rate_prototype)
  tab = Rodas3ConstantCache(constvalue(uBottomEltypeNoUnits),constvalue(tTypeNoUnits))

  tf = DiffEqDiffTools.TimeGradientWrapper(f,uprev,p)
  uf = DiffEqDiffTools.UJacobianWrapper(f,t,p)
  linsolve_tmp = zero(rate_prototype)
  linsolve = alg.linsolve(Val{:init},uf,u)
  grad_config = build_grad_config(alg,f,tf,du1,t)
  jac_config = build_jac_config(alg,f,uf,du1,uprev,u,tmp,du2)
  Rosenbrock34Cache(u,uprev,du,du1,du2,k1,k2,k3,k4,
                    fsalfirst,fsallast,dT,J,W,tmp,tab,tf,uf,linsolve_tmp,
                    linsolve,jac_config,grad_config)
end

struct Rosenbrock34ConstantCache{TF,UF,Tab} <: OrdinaryDiffEqConstantCache
  tf::TF
  uf::UF
  tab::Tab
end

function alg_cache(alg::Rodas3,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}})
  tf = DiffEqDiffTools.TimeDerivativeWrapper(f,u,p)
  uf = DiffEqDiffTools.UDerivativeWrapper(f,t,p)
  Rosenbrock34ConstantCache(tf,uf,Rodas3ConstantCache(constvalue(uBottomEltypeNoUnits),constvalue(tTypeNoUnits)))
end

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

### ROS34PW methods

@ROS34PW(:cache)

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

### ROS4 methods

@Rosenbrock4(:cache)
jac_cache(c::Rosenbrock4Cache) = (c.J,c.W)

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

### Rodas methods

struct Rodas4ConstantCache{TF,UF,Tab} <: OrdinaryDiffEqConstantCache
  tf::TF
  uf::UF
  tab::Tab
end

@cache mutable struct Rodas4Cache{uType,rateType,JType,WType,TabType,TFType,UFType,F,JCType,GCType} <: RosenbrockMutableCache
  u::uType
  uprev::uType
  dense1::rateType
  dense2::rateType
  du::rateType
  du1::rateType
  du2::rateType
  k1::rateType
  k2::rateType
  k3::rateType
  k4::rateType
  k5::rateType
  k6::rateType
  fsalfirst::rateType
  fsallast::rateType
  dT::rateType
  J::JType
  W::WType
  tmp::rateType
  tab::TabType
  tf::TFType
  uf::UFType
  linsolve_tmp::rateType
  linsolve::F
  jac_config::JCType
  grad_config::GCType
end

function alg_cache(alg::Rodas4,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  dense1 = zero(rate_prototype)
  dense2 = zero(rate_prototype)
  du = zero(rate_prototype)
  du1 = zero(rate_prototype)
  du2 = zero(rate_prototype)
  k1 = zero(rate_prototype)
  k2 = zero(rate_prototype)
  k3 = zero(rate_prototype)
  k4 = zero(rate_prototype)
  k5 = zero(rate_prototype)
  k6 = zero(rate_prototype)
  fsalfirst = zero(rate_prototype)
  fsallast = zero(rate_prototype)
  dT = zero(rate_prototype)
  if DiffEqBase.has_jac(f) && !DiffEqBase.has_Wfact(f) && f.jac_prototype !== nothing
    W = WOperator(f, dt, true)
    J = nothing # is J = W.J better?
  else
    J = false .* vec(rate_prototype) .* vec(rate_prototype)' # uEltype?
    W = similar(J)
  end
  tmp = zero(rate_prototype)
  tab = Rodas4ConstantCache(constvalue(uBottomEltypeNoUnits),constvalue(tTypeNoUnits))

  tf = DiffEqDiffTools.TimeGradientWrapper(f,uprev,p)
  uf = DiffEqDiffTools.UJacobianWrapper(f,t,p)
  linsolve_tmp = zero(rate_prototype)
  linsolve = alg.linsolve(Val{:init},uf,u)
  grad_config = build_grad_config(alg,f,tf,du1,t)
  jac_config = build_jac_config(alg,f,uf,du1,uprev,u,tmp,du2)
  Rodas4Cache(u,uprev,dense1,dense2,du,du1,du2,k1,k2,k3,k4,
                    k5,k6,
                    fsalfirst,fsallast,dT,J,W,tmp,tab,tf,uf,linsolve_tmp,
                    linsolve,jac_config,grad_config)
end

function alg_cache(alg::Rodas4,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}})
  tf = DiffEqDiffTools.TimeDerivativeWrapper(f,u,p)
  uf = DiffEqDiffTools.UDerivativeWrapper(f,t,p)
  Rodas4ConstantCache(tf,uf,Rodas4ConstantCache(constvalue(uBottomEltypeNoUnits),constvalue(tTypeNoUnits)))
end

function alg_cache(alg::Rodas42,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  dense1 = zero(rate_prototype)
  dense2 = zero(rate_prototype)
  du = zero(rate_prototype)
  du1 = zero(rate_prototype)
  du2 = zero(rate_prototype)
  k1 = zero(rate_prototype)
  k2 = zero(rate_prototype)
  k3 = zero(rate_prototype)
  k4 = zero(rate_prototype)
  k5 = zero(rate_prototype)
  k6 = zero(rate_prototype)
  fsalfirst = zero(rate_prototype)
  fsallast = zero(rate_prototype)
  dT = zero(rate_prototype)
  if DiffEqBase.has_jac(f) && !DiffEqBase.has_Wfact(f) && f.jac_prototype !== nothing
    W = WOperator(f, dt, true)
    J = nothing # is J = W.J better?
  else
    J = false .* vec(rate_prototype) .* vec(rate_prototype)' # uEltype?
    W = similar(J)
  end
  tmp = zero(rate_prototype)
  tab = Rodas42ConstantCache(constvalue(uBottomEltypeNoUnits),constvalue(tTypeNoUnits))

  tf = DiffEqDiffTools.TimeGradientWrapper(f,uprev,p)
  uf = DiffEqDiffTools.UJacobianWrapper(f,t,p)
  linsolve_tmp = zero(rate_prototype)
  linsolve = alg.linsolve(Val{:init},uf,u)
  grad_config = build_grad_config(alg,f,tf,du1,t)
  jac_config = build_jac_config(alg,f,uf,du1,uprev,u,tmp,du2)
  Rodas4Cache(u,uprev,dense1,dense2,du,du1,du2,k1,k2,k3,k4,
                    k5,k6,
                    fsalfirst,fsallast,dT,J,W,tmp,tab,tf,uf,linsolve_tmp,
                    linsolve,jac_config,grad_config)
end

function alg_cache(alg::Rodas42,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}})
  tf = DiffEqDiffTools.TimeDerivativeWrapper(f,u,p)
  uf = DiffEqDiffTools.UDerivativeWrapper(f,t,p)
  Rodas4ConstantCache(tf,uf,Rodas42ConstantCache(constvalue(uBottomEltypeNoUnits),constvalue(tTypeNoUnits)))
end

function alg_cache(alg::Rodas4P,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  dense1 = zero(rate_prototype)
  dense2 = zero(rate_prototype)
  du = zero(rate_prototype)
  du1 = zero(rate_prototype)
  du2 = zero(rate_prototype)
  k1 = zero(rate_prototype)
  k2 = zero(rate_prototype)
  k3 = zero(rate_prototype)
  k4 = zero(rate_prototype)
  k5 = zero(rate_prototype)
  k6 = zero(rate_prototype)
  fsalfirst = zero(rate_prototype)
  fsallast = zero(rate_prototype)
  dT = zero(rate_prototype)
  if DiffEqBase.has_jac(f) && !DiffEqBase.has_Wfact(f) && f.jac_prototype !== nothing
    W = WOperator(f, dt, true)
    J = nothing # is J = W.J better?
  else
    J = false .* vec(rate_prototype) .* vec(rate_prototype)' # uEltype?
    W = similar(J)
  end
  tmp = zero(rate_prototype)
  tab = Rodas4PConstantCache(constvalue(uBottomEltypeNoUnits),constvalue(tTypeNoUnits))

  tf = DiffEqDiffTools.TimeGradientWrapper(f,uprev,p)
  uf = DiffEqDiffTools.UJacobianWrapper(f,t,p)
  linsolve_tmp = zero(rate_prototype)
  linsolve = alg.linsolve(Val{:init},uf,u)
  grad_config = build_grad_config(alg,f,tf,du1,t)
  jac_config = build_jac_config(alg,f,uf,du1,uprev,u,tmp,du2)
  Rodas4Cache(u,uprev,dense1,dense2,du,du1,du2,k1,k2,k3,k4,
                    k5,k6,
                    fsalfirst,fsallast,dT,J,W,tmp,tab,tf,uf,linsolve_tmp,
                    linsolve,jac_config,grad_config)
end

function alg_cache(alg::Rodas4P,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}})
  tf = DiffEqDiffTools.TimeDerivativeWrapper(f,u,p)
  uf = DiffEqDiffTools.UDerivativeWrapper(f,t,p)
  Rodas4ConstantCache(tf,uf,Rodas4PConstantCache(constvalue(uBottomEltypeNoUnits),constvalue(tTypeNoUnits)))
end

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

### Rosenbrock5

struct Rosenbrock5ConstantCache{TF,UF,Tab} <: OrdinaryDiffEqConstantCache
  tf::TF
  uf::UF
  tab::Tab
end

@cache mutable struct Rosenbrock5Cache{uType,rateType,JType,WType,TabType,TFType,UFType,F,JCType,GCType} <: RosenbrockMutableCache
  u::uType
  uprev::uType
  dense1::rateType
  dense2::rateType
  du::rateType
  du1::rateType
  du2::rateType
  k1::rateType
  k2::rateType
  k3::rateType
  k4::rateType
  k5::rateType
  k6::rateType
  k7::rateType
  k8::rateType
  fsalfirst::rateType
  fsallast::rateType
  dT::rateType
  J::JType
  W::WType
  tmp::rateType
  tab::TabType
  tf::TFType
  uf::UFType
  linsolve_tmp::rateType
  linsolve::F
  jac_config::JCType
  grad_config::GCType
end

function alg_cache(alg::Rodas5,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  dense1 = zero(rate_prototype)
  dense2 = zero(rate_prototype)
  du = zero(rate_prototype)
  du1 = zero(rate_prototype)
  du2 = zero(rate_prototype)
  k1 = zero(rate_prototype)
  k2 = zero(rate_prototype)
  k3 = zero(rate_prototype)
  k4 = zero(rate_prototype)
  k5 = zero(rate_prototype)
  k6 = zero(rate_prototype)
  k7 = zero(rate_prototype)
  k8 = zero(rate_prototype)
  fsalfirst = zero(rate_prototype)
  fsallast = zero(rate_prototype)
  dT = zero(rate_prototype)
  if DiffEqBase.has_jac(f) && !DiffEqBase.has_Wfact(f) && f.jac_prototype !== nothing
    W = WOperator(f, dt, true)
    J = nothing # is J = W.J better?
  else
    J = false .* vec(rate_prototype) .* vec(rate_prototype)' # uEltype?
    W = similar(J)
  end
  tmp = zero(rate_prototype)
  tab = Rodas5ConstantCache(constvalue(uBottomEltypeNoUnits),constvalue(tTypeNoUnits))

  tf = DiffEqDiffTools.TimeGradientWrapper(f,uprev,p)
  uf = DiffEqDiffTools.UJacobianWrapper(f,t,p)
  linsolve_tmp = zero(rate_prototype)
  linsolve = alg.linsolve(Val{:init},uf,u)
  grad_config = build_grad_config(alg,f,tf,du1,t)
  jac_config = build_jac_config(alg,f,uf,du1,uprev,u,tmp,du2)
  Rosenbrock5Cache(u,uprev,dense1,dense2,du,du1,du2,k1,k2,k3,k4,
                    k5,k6,k7,k8,
                    fsalfirst,fsallast,dT,J,W,tmp,tab,tf,uf,linsolve_tmp,
                    linsolve,jac_config,grad_config)
end

function alg_cache(alg::Rodas5,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}})
  tf = DiffEqDiffTools.TimeDerivativeWrapper(f,u,p)
  uf = DiffEqDiffTools.UDerivativeWrapper(f,t,p)
  Rosenbrock5ConstantCache(tf,uf,Rodas5ConstantCache(constvalue(uBottomEltypeNoUnits),constvalue(tTypeNoUnits)))
end

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

### RosenbrockW6S4O

@RosenbrockW6S4OS(:cache)

ext: jl


text:
module CharLM
using Knet, ArgParse, JLD

function main(args=ARGS)
    global net, vocab, text, data
    s = ArgParseSettings()
    s.description="charlm.jl (c) Deniz Yuret, 2016. This is the character-level language model from a popular blog post 'The Unreasonable Effectiveness of Recurrent Neural Networks' by Andrej Karpathy, 2015. (http://karpathy.github.io/2015/05/21/rnn-effectiveness)"
    s.exc_handler=ArgParse.debug_handler
    @add_arg_table s begin
        ("--datafiles"; nargs='+'; help="if provided, use first file for training, second for dev, others for test.")
        ("--loadfile"; help="initialize model from file")
        ("--savefile"; help="save final model to file")
        ("--bestfile"; help="save best model to file")
        ("--generate"; arg_type=Int; default=0; help="if non-zero generate given number of bytes.")
        ("--nlayer"; arg_type=Int; default=1)
        ("--hidden"; arg_type=Int; default=256)
        ("--embedding"; arg_type=Int; default=256)
        ("--epochs"; arg_type=Int; default=1)
        ("--batchsize"; arg_type=Int; default=128)
        ("--seqlength"; arg_type=Int; default=100)
        ("--decay"; arg_type=Float64; default=0.9)
        ("--lr"; arg_type=Float64; default=1.0)
        ("--gclip"; arg_type=Float64; default=5.0)
        ("--dropout"; arg_type=Float64; default=0.0)
        ("--seed"; arg_type=Int; default=42)
    end
    isa(args, AbstractString) && (args=split(args))
    o = parse_args(args, s; as_symbols=true); println(o)
    o[:seed] > 0 && setseed(o[:seed])
    text = !isempty(o[:datafiles]) ? map(f->readall(f), o[:datafiles]) : cell(0)
    !isempty(text) && info("Chars read: $((map(length,text)...))")
    vocab = o[:loadfile]!=nothing ? load(o[:loadfile], "vocab") : Dict{Char,Int}()
    if !isempty(text)
        if isempty(vocab)
            for t in text, c in t; get!(vocab, c, 1+length(vocab)); end
        else
            for t in text, c in t; haskey(vocab, c) || error("Unknown char $c"); end
        end
    end
    o[:vocabsize] = length(vocab); info("$(o[:vocabsize]) unique chars")
    net = o[:loadfile]!=nothing ? load(o[:loadfile], "net") : compile(:charlm; o...)
    !isempty(text) && bigtrain(net, text, vocab, o)
    o[:generate] > 0 && generate(net, vocab, o)
end

function generate(f, vocab, o)
    alpha = Array(Char, length(vocab))
    for (c,i) in vocab; alpha[i]=c; end
    x=zeros(Float32, length(alpha), 1)
    y=zeros(Float32, length(alpha), 1)
    xi = 1
    for i=1:o[:generate]
        copy!(y, forw(f,x))
        x[xi] = 0
        r = rand(Float32)
        p = 0
        for c=1:length(y)
            p += y[c]
            r <= p && (xi=c; break)
        end
        x[xi] = 1
        print(alpha[xi])
    end
    println()
end

function bigtrain(net, text, vocab, o)
    data = map(t->minibatch(t, vocab; o...), text)
    setp(net, lr=o[:lr])
    loss = zeros(length(data))
    lastloss = bestloss = Inf
    for epoch=1:o[:epochs]
        for d=1:length(data)
            loss[d] = (d==1 ? train(net, data[d], softloss; o...) : test(net, data[d], softloss))
        end
        println((epoch, o[:lr], loss...)); flush(STDOUT)
        if length(data) > 1
            loss[2] > lastloss && (o[:lr] *= o[:decay]; setp(net, lr=o[:lr]))
            lastloss = loss[2]
            if lastloss < bestloss
                bestloss = lastloss
                o[:bestfile]!=nothing && save(o[:bestfile], "net", clean(net), "vocab", vocab)
            end
        end
    end
    o[:savefile]!=nothing && save(o[:savefile], "net", clean(net), "vocab", vocab)
    return bestloss
end

@knet function charlm(x; embedding=0, nlayer=0, vocabsize=0, o...)
    y = wdot(x; out=embedding)
    z = repeat(y; o..., frepeat=:charlayer, nrepeat=nlayer)
    return wbf(z; out=vocabsize, f=:soft)
end

@knet function charlayer(x; dropout=0, hidden=0, o...)
    y = lstm(x; out=hidden)
    return drop(y; pdrop=dropout)
end

function test(f, data, loss)
    sumloss = 0.0
    reset!(f)
    T = length(data)-1
    for t=1:T
        x = data[t]
        ygold = data[t+1]
        ypred = forw(f,x)
        sumloss += loss(ypred, ygold)
    end
    sumloss/T
end

function train(f, data, loss; gclip=0, seqlength=100, dropout=0, o...)
    ystack = cell(0)
    sumloss = 0.0
    reset!(f, keepstate=false)
    T = length(data)-1
    for t=1:T
        x = data[t]
        ygold = data[t+1]
        ypred = sforw(f,x; dropout=(dropout>0))
        sumloss += loss(ypred, ygold)
        push!(ystack,ygold)
        if (t%seqlength == 0 || t==T)
            while !isempty(ystack)
                ygold = pop!(ystack)
                sback(f,ygold,loss)
            end
            g = (gclip > 0 ? gnorm(f) : 0)
            update!(f; gscale=(g > gclip > 0 ? gclip/g : 1))
            reset!(f, keepstate=true)
            # @show (t, sumloss/t)
        end
    end
    sumloss/T
end

# Create minibatches.  Each minibatch is a (vocabsize, batchsize)
# matrix with one-hot columns.  The vocab dictionary maps characters
# to integer indices.  The data array returned will have
# T=|chars|/batchsize minibatches.  The columns of minibatch t refer
# to characters t, t+T, t+2T etc.  During training if x=data[t], then
# y=data[t+1].

function minibatch(chars, vocab; vocabsize=0, batchsize=0, o...)
    data = cell(0)
    T = div(length(chars), batchsize)
    for t=1:T
        d=zeros(Float32, vocabsize, batchsize)
        for b=1:batchsize
            c = vocab[chars[t + (b-1) * T]]
            d[c,b] = 1
        end
        push!(data, d)
    end
    return data
end

# temp workaround: prevents error in running finalizer: ErrorException("auto_unbox: unable to determine argument type")
@gpu atexit(()->(for r in net.reg; r.out0!=nothing && Main.CUDArt.free(r.out0); end))

!isinteractive() && !isdefined(Core.Main, :load_only) && main(ARGS)

end # module

# Complete works of Shakespeare:
# http://www.gutenberg.org/files/100/100.txt

# Karpathy blog:
# http://karpathy.github.io/2015/05/21/rnn-effectiveness/
#  we will encode each character into a vector using 1-of-k encoding.
#  The RNN is trained with mini-batch Stochastic Gradient Descent. I like to use RMSProp (a per-parameter adaptive learning rate) to stablilize the updates.
#  At test time, we feed a character into the RNN and get a distribution over what characters are likely to come next. We sample from this distribution, and feed it right back in to get the next letter.
#  PaulGraham:  Lets train a 2-layer LSTM with 512 hidden nodes (approx. 3.5 million parameters), and with dropout of 0.5 after each layer. We'll train with batches of 100 examples and truncated backpropagation through time of length 100 characters.
#  Shakespeare: 3-layer RNN with 512 hidden nodes on each layer
#  temperature? you can also generate an infinite amount of your own samples at different temperatures with the provided code.

# const path = Pkg.dir("Knet/data/100.txt")
# const url = "http://www.gutenberg.org/files/100/100.txt"
# isfile(path) || (info("Downloading $url"); save(get(url), path))

# Parameters to play with:
# lr=0.001
# gclip=0
# network depth=1
# network width=512
# fbias=1 # 0 did not help
# winit=Xavier # Gaussian did not help
# embedding matrix=none
# dropout=none
# keepstate=no
# rmsprop=no

# start score:  	(3.2801345440473475,131202.10162734985,39999)
# keepstate=true: 	(3.2360347558469806,129438.15419912338,39999)
# lr=0.002		(3.217129004007876,128681.94303131104,39999)
# lr=0.005		(3.16123230188065,126446.13084292412,39999)
# lr=0.01		(2.9789575188841657,119155.32179784775,39999)
# lr=0.02		(12.209143830274368,488353.5440671444,39999) blows up need gradient clip
# lr=0.02 gclip=10	(3.253291426058063,130128.40375089645,39999)
# lr=2.0 gclip=1.0	(2.7846112851396305,111381.66679430008,39999)
# lr=2.0 gclip=5.0	(2.842096986512893,113681.0373635292,39999)
# lr=1.0 gclip=5.0	(2.685244952175879,107407.11284208298,39999)
# lr=0.5 gclip=5.0	(2.726828870917801,109070.42800784111,39999)
# lr=1.0 gclip=10.0	(2.8446781882471064,113784.28285169601,39999)
# lr=1.0 gclip=2.0	(2.7840308177804656,111358.44868040085,39999)
# fixing at lr=1.0 gclip=5.0
# H=1024		(2.988540398555683,119538.62740182877,39999)
# H=256			(2.5366065015989308,101461.72345745564,39999) but this is for one epoch not convergence. and train not dev.

# with dev set: lr=1.0 gclip=5.0 (ep,lr,trn,dev)
# H=512	(1,1.0,2.6913729234082586,2.0592177872112076)
# H=256	(1,1.0,2.5404179560464306,1.9690401755821478)
# H=128 (1,1.0,2.452683317178178,1.955850279233761)

# fix H=128 play with embedding E
# E=16 (1,1.0,2.260463515062112,1.8129167987150456)
# E=32 (1,1.0,2.2403026809662343,1.8023061217339227)
# E=64 (1,1.0,2.2276982673749544,1.7939553815922231)
# E=128 (1,1.0,2.21924334593858,1.7907175187843698)

# H=256 E=128 (1,1.0,2.2865155415581584,1.7820057602604338)
# H=E=512 (1,1.0,2.3814222320723157,1.7827530246014816)

# fix H=E=256
# H=E=256 (1,1.0,2.274747170884358,1.7639083232801682)
# 2-layer: (1,1.0,2.814780586749161,2.1773407803894065)


text:
using Documenter, TerraDG
#https://juliadocs.github.io/Documenter.jl/stable/man/guide/index.html
# Run this by julia --project=. --color=yes docs/make.jl in the root
# of this repository.
makedocs(
    modules = [TerraDG],
    sitename="TerraDG Documentation")


text:
filter_ext(path, ext) = filter(filename -> splitext(filename)[2] == ext, readdir(path))

maybewindow2indices(A, window::Tuple) =
    maybewindow2indices(A, dims(A), window::Tuple)
maybewindow2indices(A, dims::Tuple, window::Tuple) =
    window == () ? () : to_indices(A, dims2indices(dims, window))

# Read from the paraent dataset, using the window indices if they exist
# We try to load as little data from disk as possible.
readwindowed(A::AbstractGeoArray, window::Tuple{}) = GeoArray(A)
readwindowed(A, window::Tuple{}) = Array(A)
readwindowed(A, window::Tuple{}, I...) = A[I...]
readwindowed(A, window::Tuple, I...) = A[Base.reindex(window, I)...]
readwindowed(A, window::Tuple) = readwindowed(A, window...)
readwindowed(A, i, I...) = A[i, I...]
readwindowed(A) = Array(A)

# Get a metadata field
getmeta(A::AbstractGeoArray, key, fallback) = getmeta(metadata(A), key, fallback)
getmeta(m::Metadata, key, fallback) = get(val(m), key, fallback)
getmeta(m::Nothing, key, fallback) = fallback

# Check that arrayu order amtches expectation
checkarrayorder(A, order::Order) = map(d -> checkarrayorder(d, order), dims(A))
checkarrayorder(A, order::Tuple) = map(checkarrayorder, dims(A), order)
checkarrayorder(dim::Dimension, order::Order) = begin
    ao = arrayorder(dim)
    ao == order || @warn "Array order for `$(DD.basetypeof(order))` is `$ao`, usualy `$order`"
end

cleankeys(keys) = Tuple(Symbol.(keys))

#=
Shift the index from the current loci to the new loci. We only actually
shift Regular Intervals, and do this my multiplying the offset of
-1, -0.5, 0, 0.5 or 1 by the absolute value of the span.

TODO: move this to DimensionalData.jl
=#
shiftindexloci(locus::Locus, dim::Dimension) = shiftindexloci(mode(dim), locus, dim)
shiftindexloci(::IndexMode, ::Locus, dim::Dimension) = dim
shiftindexloci(mode::AbstractSampled, locus::Locus, dim::Dimension) =
    shiftindexloci(span(mode), sampling(mode), locus, dim)
shiftindexloci(span::Span, sampling::Sampling, ::Locus, dim::Dimension) = dim
shiftindexloci(span::Regular, sampling::Intervals, destlocus::Locus, dim::Dimension) =
    rebuild(dim, val(dim) .+ (abs(step(span)) * offset(locus(sampling), destlocus)))

offset(::Start, ::Center) = 0.5
offset(::Start, ::End) = 1
offset(::Center, ::Start) = -0.5
offset(::Center, ::End) = 0.5
offset(::End, ::Start) = -1
offset(::End, ::Center) = -0.5
offset(::T, ::T) where T<:Locus = 0


text:
"""
    ReplitDatabaseCore

The core functions which access the database.
This implementation accesses the database using `curl` via shell calls.
The `ReplitDB` interface type depends on these functions.

## Exports

- `get(key; url)`
- `set!(key, value; url)`
- `list([prefix]; url)`
- `delete!(key; url)`
"""
module ReplitDatabaseCore
export delete!, get, list, set!

using HTTP
using URIs: escapeuri, unescapeuri

"""
    delete!(key::AbstractString; url::AbstractString=ENV["REPLIT_DB_URL"])

Deletes the value associated with the given key in the Repl.it database.
"""
function delete!(key::AbstractString; url::AbstractString=ENV["REPLIT_DB_URL"])
    HTTP.request("DELETE", "$(url)/$(escapeuri(key))")
end

"""
    get(key::AbstractString; url::AbstractString=ENV["REPLIT_DB_URL"])

Get the value associated with the given key in the Repl.it database.
"""
function get(key::AbstractString; url::AbstractString=ENV["REPLIT_DB_URL"])
    try
        String(HTTP.get("$(url)/$(escapeuri(key))").body)
    catch err
        if err.status == 404
            throw(KeyError(key))
        else
            throw(err)
        end
    end
end

"""
    list(; url::AbstractString=ENV["REPLIT_DB_URL"])
    list(prefix::AbstractString; url::AbstractString=ENV["REPLIT_DB_URL"])

Lists all of the keys in the Repl.it database, or all of the keys starting with `prefix` if specified.
"""
function list(prefix=""::AbstractString; url::AbstractString=ENV["REPLIT_DB_URL"])::Vector{String}
    escapedurl = "$(url)?prefix=$(escapeuri(prefix))&encode=true"
    result = String(HTTP.get(escapedurl).body)
    isempty(result) ? String[] : unescapeuri.(split(result, '\n'))
end

"""
    set!(key::AbstractString, value::AbstractString; url::AbstractString=ENV["REPLIT_DB_URL"])

Set the given key/value pair in the Repl.it database.
"""
function set!(key::AbstractString, value::AbstractString; url::AbstractString=ENV["REPLIT_DB_URL"])
    # success(`curl -s $(url) -d $(key * "=" * value)`)
    try
        response = HTTP.post(url, ["Content-Type" => "application/x-www-form-urlencoded"], "$(key)=$(value)")
    catch response
    finally
        response.status == 200
    end
end

end # module


text:
#############
# In this example, we sample a bivariate function:
#
#   f(x,y) = 1/(1+x^2+y^2),
#
# on the reference triangle with vertices (0,0), (0,1), and (1,0) and analyze it
# in a Proriol series. Then, we find Proriol series for each component of its
# gradient by term-by-term differentiation of our expansion, and we compare them
# with the true Proriol series by sampling an exact expression for the gradient.
#
# We analyze f(x,y) on an N×M mapped tensor product grid defined by:
#
#   x = (1+u)/2, and y = (1-u)*(1+v)/4, where:
#
#   uₙ = cos[(n+1/2)π/N], for 0 ≤ n < N, and
#
#   vₘ = cos[(m+1/2)π/M], for 0 ≤ m < M;
#
# we convert the function samples to mapped Chebyshev² coefficients using
# `plan_tri_analysis`; and finally, we transform the mapped Chebyshev²
# coefficients to Proriol coefficients using `plan_tri2cheb`.
#
# For the storage pattern of the arrays, please consult the documentation.
#############

using FastTransforms, LinearAlgebra

f = (x,y) -> 1/(1+x^2+y^2)
fx = (x,y) -> -2x/(1+x^2+y^2)^2
fy = (x,y) -> -2y/(1+x^2+y^2)^2

N = 10
M = N

α, β, γ = 0, 0, 0

u = [sinpi((N-2n-1)/(2N)) for n in 0:N-1]
v = [sinpi((M-2m-1)/(2M)) for m in 0:M-1]

# Instead of using the u, v grid, we use one with more accuracy near the origin.
x = [sinpi((2N-2n-1)/(4N))^2 for n in 0:N-1]
w = [sinpi((2M-2m-1)/(4M))^2 for m in 0:M-1]

(1 .+ u)./2 ≈ x
(1 .- u).*(1 .+ v')/4 ≈ reverse(x).*w'

# On the mapped tensor product grid, our function samples are:
F = [f(x[n+1], x[N-n]*w[m+1]) for n in 0:N-1, m in 0:M-1]

P = plan_tri2cheb(F, α, β, γ)
PA = plan_tri_analysis(F)

# Its Proriol-(α,β,γ) coefficients are:
U = P\(PA*F)

# Similarly, our function's gradient samples are:
Fx = [fx(x[n+1], x[N-n]*w[m+1]) for n in 0:N-1, m in 0:M-1]
Fy = [fy(x[n+1], x[N-n]*w[m+1]) for n in 0:N-1, m in 0:M-1]

# For the partial derivative with respect to x, Olver et al.
# derive simple expressions for the representation of this component
# using a Proriol-(α+1,β,γ+1) series. For the partial derivative with respect
# to y, the analogous formulae result in a Proriol-(α,β+1,γ+1) series.
# These expressions are adapted from https://arxiv.org/abs/1902.04863.
Gx = zeros(Float64, N, M)
for m = 0:M-2
    for n = 0:N-2
        cf1 = m == 0 ? sqrt((n+1)*(n+2m+α+β+γ+3)/(2m+β+γ+2)*(m+γ+1)*8) : sqrt((n+1)*(n+2m+α+β+γ+3)/(2m+β+γ+1)*(m+β+γ+1)/(2m+β+γ+2)*(m+γ+1)*8)
        cf2 = sqrt((n+α+1)*(m+1)/(2m+β+γ+2)*(m+β+1)/(2m+β+γ+3)*(n+2m+β+γ+3)*8)
        Gx[n+1, m+1] = cf1*U[n+2, m+1] + cf2*U[n+1, m+2]
    end
end
Px = plan_tri2cheb(Fx, α+1, β, γ+1)
Ux = Px\(PA*Fx)

Gy = zeros(Float64, N, M)
for m = 0:M-2
    for n = 0:N-2
        Gy[n+1, m+1] = 4*sqrt((m+1)*(m+β+γ+2))*U[n+1, m+2]
    end
end
Py = plan_tri2cheb(Fy, α, β+1, γ+1)
Uy = Py\(PA*Fy)

# The 2-norm relative error in differentiating the Proriol series
# for f(x,y) term-by-term and its sampled gradient is:
hypot(norm(Ux-Gx), norm(Uy-Gy))/hypot(norm(Ux), norm(Uy))

# This error can be improved upon by increasing N and M.


text:
using SpecialFunctions
# using LinearAlgebra

l2norm(x) = sqrt(sum(abs2, x))

function λ(κ,m, μ,ψ; n=1)
    cosψ = cos(ψ)
    tanψ = tan(ψ)
    if n != 1
        numerator = 2^n * π^((n-1)/2) * ((μ / cosψ) - m*im ) * gamma((n+1)/2) * ((tanψ)^(n-1))
        denominator = (((μ / cosψ) - m*im )^2 + (l2norm(κ) * tanψ)^2)^((n+1)/2) * cosψ
    elseif n == 1
        numerator = 2 * ((μ / cosψ) - m*im )
        denominator = (((μ / cosψ) - m*im )^2 + (l2norm(κ) * tanψ)^2) * cosψ
    end
    return numerator / denominator
end

function Φ(κ,m, x,z)
    return exp((κ'x + (m*z)) * im)
end


text:
# more types for algorhythm
# abstract type capitalaction end
# abstract type active <: capitalaction end
# abstract type passive <: capitalaction end
# # only use capital action for multiplechoicevar. will try to code without

"""ItpArray allows the array to be interpolated. Remove when Interpolations.jl
supports array indexing with `[]`."""
struct ItpArray{T,N,Tarr<:Array{T,N},
		Titp<:AbstractInterpolation{T,N}} <: AbstractArray{T,N}
	arr::Tarr
	itp::Titp
end
# overload array functions
Base.size(A::ItpArray) = size(A.arr)
Base.getindex(A::ItpArray, args...) = getindex(A.arr,args...)
Base.setindex!(A::ItpArray, args...) = setindex!(A.arr, args...)
# overload interpolator
function (A::ItpArray)(x::Vararg{Number,N}) where {T,N}
	return A.itp(x...) #A.itp()(x) # size(A.itp)
end
function itparray(A::AbstractArray{T,N}, knots::NTuple{N,AbstractVector{T}}) where {T,N}
	itp = interpolate(knots, A, Gridded(Linear()))
	itp = extrapolate(itp, Flat())
	ItpArray{T,N,typeof(A),typeof(itp)}(A,itp)
end

abstract type AbstractDDPSolution{NS,NC} end

"""Solution object.
Parameters:
    - NS = number of states variables
    - NC = number of choice variables"""
struct DDPSolution{NS,NC,valueType,policyType} <: AbstractDDPSolution{NS,NC}

    prob::DDP
    value::valueType
    policy::NTuple{NC,policyType}

end

"""Includes solution for exact initialization.
Parameters:
    - NE0 = number of exogenous states variables at t=0
    - NC0 = number of choice variables at t=0"""
struct DDPSolutionZero{NS,NC,NE0,NC0,valueType,policyType,
		value0Type,policy0Type} <: AbstractDDPSolution{NS,NC}

	prob::DDP
	value::valueType
	policy::NTuple{NC,policyType}
	value0::value0Type
	policy0::NTuple{NC0,policy0Type}

end

function createsolution(p::DDP, meshValFun::Array{T,NS},
                                tmeshPolFun::NTuple{NC,Array{T,NS}}) where {T,NS,NC}

	value = itparray(meshValFun, p.tStateVectors)
  	policy = [itparray(pol, p.tStateVectors) for pol in tmeshPolFun]
  	policy = tuple(policy...)

   if p.options.initialize == nothing
      return DDPSolution{NS,NC,typeof(value),typeof(policy[1])}(p, value, policy)
   else
      meshValFunZero, tmeshPolFunZero = initialendogstatevars(p, meshValFun)
		tExogStateVectorsZero = getnonchoicevarszero(p)
		NE0 = length(tExogStateVectorsZero)
		NC0 = length(p.options.initialize.tChoiceVectorsZero)
		value0 = itparray(meshValFunZero, tExogStateVectorsZero)
		policy0 = [itparray(pol, tExogStateVectorsZero) for pol in tmeshPolFunZero]
		policy0 = tuple(policy0...)
      return DDPSolutionZero{NS,NC,NE0,NC0,typeof(value),typeof(policy[1]),
			typeof(value0),typeof(policy0[1])}(p, value, policy, value0, policy0)
   end
end

value(sol::AbstractDDPSolution) = sol.value
policy(sol::AbstractDDPSolution) = sol.policy
policy(sol::AbstractDDPSolution{NS,1}) where NS = sol.policy[1]

value0(sol::DDPSolutionZero) = sol.value0
policy0(sol::DDPSolutionZero) = sol.policy0
policy0(sol::DDPSolutionZero{NS,1}) where NS  = sol.policy0[1]

"""Compare two DDP solutions. Returns true when they are approximately the same,
otherwise false."""
function isapprox(sol1::AbstractDDPSolution,
        sol2::AbstractDDPSolution; kwargs...)

	test = []
    for field in fieldnames(typeof(sol1))
		field1 = getfield(sol1,field)
		field2 = getfield(sol2,field)
        if typeof(field1) <: Tuple
            append!(test,[isapprox(field1[i], field2[i]; kwargs...) for i = 1 : length(field1)])
        elseif typeof(field1) <: AbstractArray
            append!(test,isapprox(field1, field2; kwargs...))
        end
    end
	return .!any(.!test)
end

# function compare(sol1::AbstractDDPSolution,
#         sol2::AbstractDDPSolution;
# 		str::String = "compare solutions",
#         tol::Real = 1e-4)
#     @testset "$str" begin
#         for field in fieldnames(typeof(sol1))
# 			field1 = getfield(sol1,field)
# 			field2 = getfield(sol2,field)
#             if typeof(field1) <: Tuple
#                 [@test isapprox(field1[i], field2[i], rtol=tol) for i = 1 : length(field1)]
#             elseif typeof(field1) <: Array
#                 @test isapprox(field1, field2, rtol=tol)
#             end
#         end
#     end
# end


text:
using LinearAlgebra;

datadir=joinpath("data","exp");
graph_m6_SID=import_compgraph(joinpath(datadir,"exp_sid_m6_opt_rho2_22.cgr"));
graph_m7_mono=import_compgraph(joinpath(datadir,"exp_mono_m7_opt_rho6_4.cgr"));
graph_m8_mono=import_compgraph(joinpath(datadir,"exp_mono_m8_opt_rho13_5.cgr"));

# To create the corresponding graph_exp_native
AA=ones(Float64,1,1)
(graph_native,_)=graph_exp_native_jl(AA*2.5);
(graph_native2,_)=graph_exp_native_jl(AA*5.5);
(graph_native3,_)=graph_exp_native_jl(AA*13.5);

runs=["expm_matlab"             :expm_matlab        2.5
      "expm_matlab"             :expm_matlab        6.0
      "expm_matlab"             :expm_matlab        13.5
      "exp_julia"               :exp_julia          2.5
      "exp_julia"               :exp_julia          6.0
      "exp_julia"               :exp_julia          13.5
      "expmpoly_matlab"         :expmpoly_matlab    2.5
      "expmpoly_matlab"         :expmpoly_matlab    6.0
      "expmpoly_matlab"         :expmpoly_matlab    13.5
      "sid_m6_opt_rho2_7"       graph_m6_SID        2.5
      "mono_m7_opt_rho6_4"      graph_m7_mono       6.0
      "mono_m8_opt_rho13_5"     graph_m8_mono       2.5
      "native_73_jl"            graph_native        2.5
      "native_83_jl"            graph_native2       6.0
      "native_93_jl"            graph_native3       13.5];

names=runs[:,1];
graphs=runs[:,2];
normvalues=runs[:,3];
for (i,g)=enumerate(graphs)
    if (g isa Compgraph)
        if (count(values(g.operations) .== :ldiv) == 0)
            degopt=Degopt(g);
            normalize!(degopt);
            (g,_)=graph_degopt(degopt);
        end
        x=get_coeffs(g)
        if (norm(imag.(x))<eps()*100)
            set_coeffs!(g,real.(x));
        else
            error("non-real coeffs");
        end
        g=Compgraph(Float64,g);
        compress_graph!(g);
        graphs[i]=g;
    end
end

# To reduce the memory footprint
priohelp=Dict();
priohelp[:B2]=-1e8;
priohelp[:Ba3_2]=-90;
priohelp[:Ba3]=-90;
priohelp[:Bb3_2]=-90;
priohelp[:Bb3]=-90;
priohelp[:B3]=-90;

priohelp[:Ba4_2]=-80;
priohelp[:Ba4_3]=-80;
priohelp[:Ba4]=-80;
priohelp[:Bb4_2]=-80;
priohelp[:Bb4_3]=-80;
priohelp[:Bb4]=-80;
priohelp[:B4]=-80;

priohelp[:Ba5_2]=-70;
priohelp[:Ba5_3]=-70;
priohelp[:Ba5_4]=-70;
priohelp[:Ba5]=-70;
priohelp[:Bb5_2]=-70;
priohelp[:Bb5_3]=-70;
priohelp[:Bb5_4]=-70;
priohelp[:Bb5]=-70;
priohelp[:B5]=-70;

priohelp[:Ba6_2]=-60;
priohelp[:Ba6_3]=-60;
priohelp[:Ba6_4]=-60;
priohelp[:Bb6_2]=-60;
priohelp[:Bb6_3]=-60;
priohelp[:Bb6_4]=-60;
priohelp[:Ba6]  =-60;
priohelp[:Bb6]  =-60;
priohelp[:B6]  =-60;

priohelp[:Ba7_2]=0;
priohelp[:Ba7_3]=0;
priohelp[:Bb7_2]=0;
priohelp[:Bb7_3]=0;

priohelp[:T2k2]=1000;
priohelp[:T2k3]=1000;
priohelp[:T2k4]=1000;
priohelp[:T2k5]=1000;

#
#addit=-500;
#factor=10;
#for (i,s)=enumerate(get_topo_order_degopt(6))
#    priohelp[s]=addit+factor*i;
#end
#

#priohelp=Dict();
#addit=-10;
#priohelp[:Ba2]=-10+addit;
#priohelp[:Bb2]=-10+addit;
#priohelp[:B2]=-10+addit;
#priohelp[:Ba3]=-9+addit;
#priohelp[:Ba3_2]=-9+addit;
#priohelp[:Bb3]=-8+addit;
#priohelp[:Bb3_2]=-8+addit;
#priohelp[:B3]=-8+addit;
#
#priohelp[:Ba4]=-7+addit;
#priohelp[:Ba4_2]=-7+addit;
#priohelp[:Ba4_3]=-7+addit;
#priohelp[:Bb4]=-6+addit;
#priohelp[:Bb4_2]=-6+addit;
#priohelp[:Bb4_3]=-6+addit;
#priohelp[:B4]=-6+addit;
#
#priohelp[:Ba5]=-5+addit;
#priohelp[:Ba5_2]=-5
#priohelp[:Ba5_3]=-5+addit;
#priohelp[:Ba5_4]=-5+addit;
#priohelp[:Bb5]=-4+addit;
#priohelp[:Bb5_2]=-4+addit;
#priohelp[:Bb5_3]=-4+addit;
#priohelp[:Bb5_4]=-4+addit;
#priohelp[:B5]=-4+addit;
#
#
#
#priohelp[:Ba6]=-3+addit;
#priohelp[:Ba6_2]=-3
#priohelp[:Ba6_3]=-3+addit;
#priohelp[:Ba6_4]=-3+addit;
#priohelp[:Ba6_5]=-3+addit;
#priohelp[:Bb6]=-2+addit;
#priohelp[:Bb6_2]=-2+addit;
#priohelp[:Bb6_3]=-2+addit;
#priohelp[:Bb6_4]=-2+addit;
#priohelp[:Bb6_5]=-2+addit;
#priohelp[:B6]=-2+addit;
#
#
#priohelp[:Ba7]=-1.5+addit;
#priohelp[:Ba7_2]=-1.5
#priohelp[:Ba7_3]=-1.5+addit;
#priohelp[:Ba7_4]=-1.5+addit;
#priohelp[:Ba7_5]=-1.5+addit;
#priohelp[:Ba7_6]=-1.5+addit;
#priohelp[:Bb7]=-1.1+addit;
#priohelp[:Bb7_2]=-1.1+addit;
#priohelp[:Bb7_3]=-1.1+addit;
#priohelp[:Bb7_4]=-1.1+addit;
#priohelp[:Bb7_5]=-1.1+addit;
#priohelp[:Bb7_6]=-1.1+addit;
#priohelp[:B7]=-1.1+addit;
#

max_stars_repo_name: matrixfunctions/GraphMatFunData

max_issues_repo_name: matrixfunctions/GraphMatFunData

max_forks_repo_name: matrixfunctions/GraphMatFunData


text:
GMSH_FOUND = true
gmsh = Sys.which("gmsh")
gmsh_jl = ""

if gmsh === nothing
  s = """
  No gmsh installation found in the system.
  Install gmsh and make sure that the gmsh executable
  can be found via the PATH environment variable.
  """
  @warn s
  GMSH_FOUND = false
else
  gmsh_root = dirname(gmsh)[1:end-4]
  gmsh_jl = joinpath(gmsh_root,"lib","gmsh.jl")
  @info "Using the gmsh Julia API found in $gmsh_jl"
end

if ! isfile(gmsh_jl)
  @warn "gmsh Julia API not found"
  GMSH_FOUND = false
end

open("deps.jl","w") do f
  println(f, "# This file is automatically generated")
  println(f, "# Do not edit")
  println(f)
  println(f, :(const GMSH_FOUND = $GMSH_FOUND))
  println(f, :(const gmsh_jl    = $gmsh_jl))
end


text:
module CGDycore

export strcmp

#=
Look for

# TODO: check translation with Oswald
=#

# using StatsBase
using LinearAlgebra
using SparseArrays
include("matlab_intrinsics.jl")

include("Grid/Node/Node.jl")
include("Grid/Edge/Edge.jl")
include("Grid/Face/Face.jl")
# include("Grid/Cell/Cell.jl")
include("Grid/GridStruct.jl")

# include("Grid/Annulus.jl")
# include("Grid/AnnulusGrid.jl")
# include("Grid/CFPToGrid.jl")
# include("Grid/Cart.jl")
include("Grid/CartGrid.jl")
include("Grid/CubedGrid.jl")
# include("Grid/DGNodes.jl")
# include("Grid/FacesInEdges.jl")
include("Grid/FacesInNodes.jl")
# # include("Grid/GaussLobattoQuad.jl") # delete?
# include("Grid/Grid2DCFP.jl")
# include("Grid/GridToKiteGrid.jl")
# include("Grid/GridToKiteGridNeu.jl")
# include("Grid/HexagonGrid.jl")
# include("Grid/HexagonGridNeu.jl")
# include("Grid/InputGrid.jl")
# include("Grid/JacobiAnnulus.jl")
# include("Grid/JacobiAnnulusCyl.jl")
# include("Grid/JacobiCart.jl")
# include("Grid/JacobiCart1.jl")
# include("Grid/JacobiCartV1.jl")
# include("Grid/JacobiDG.jl")
# include("Grid/JacobiDG1.jl")
# include("Grid/JacobiDG2.jl")
include("Grid/JacobiDG3.jl")
# include("Grid/JacobiShallowSphere3.jl")
# include("Grid/JacobiSphere.jl")
# include("Grid/JacobiSphere2.jl")
# include("Grid/JacobiSphere2TT.jl")
include("Grid/JacobiSphere3.jl")
# include("Grid/JacobiSphereDG.jl")
# include("Grid/JacobiTri.jl")
# include("Grid/MetricGuba.jl")
# include("Grid/MetricGubaOld.jl")
# include("Grid/MetricGubaOld2.jl")
include("Grid/OrientFaceCart.jl")
include("Grid/OrientFaceSphere.jl")
include("Grid/Orientation.jl")
# include("Grid/Plot2.jl")
# include("Grid/Plot2BBB.jl")
# include("Grid/PlotFaceGrid.jl")
# include("Grid/PlotFaceKiteGrid.jl")
# include("Grid/PolygonGrid.jl")
# include("Grid/PolygonToGrid.jl")
include("Grid/Renumbering.jl")
# include("Grid/SingleHexagonGrid.jl")
include("Grid/Topo.jl")
# include("Grid/Topo2.jl")
include("Grid/TransCart.jl")
include("Grid/TransSphere.jl")
# include("Grid/TriGrid.jl")
# include("Grid/cart2Radial.jl")
include("Grid/cart2sphere.jl")
include("Grid/hS.jl")
# include("Grid/sphere2cart.jl")
include("Grid/vtkWriteHex.jl")
# include("Grid/vtkWriteQuad.jl")

include("DG/DLagrange.jl")
# include("DG/DerivativeMatrix.jl")
include("DG/DerivativeMatrixSingle.jl")
# include("DG/GaussLegendreQuad.jl")
include("DG/GaussLobattoQuad.jl")
# include("DG/GetQuadMeth.jl")
# include("DG/GradGrad.jl")
# include("DG/InitMethodDG.jl")
include("DG/Lagrange.jl")
include("DyCore/Average.jl")
include("DyCore/AverageFB.jl")
include("DyCore/BoundaryW.jl")
# include("DyCore/BoundaryWOutput.jl")
include("DyCore/DampingKoeff.jl")
include("DyCore/Discretization.jl")
# include("DyCore/Energy.jl")
include("DyCore/FCurlNon3Vec.jl")
# include("DyCore/FDiv3UpwindVec.jl")
include("DyCore/FDiv3Vec.jl")
# include("DyCore/FDivGrad2Vec.jl")
include("DyCore/FDivGrad2VecDSS.jl")
include("DyCore/FDivRhoGrad2Vec.jl")
include("DyCore/FGrad3Vec.jl")
include("DyCore/FGradDiv2Vec.jl")
include("DyCore/FGradDiv2VecDSS.jl")
include("DyCore/FRotCurl2Vec.jl")
include("DyCore/FRotCurl2VecDSS.jl")
include("DyCore/FVort2VecDSS.jl")
# include("DyCore/FcnDiff.jl")
include("DyCore/FcnNHCurlVec.jl")
# include("DyCore/FcnVec.jl")
include("DyCore/HyperDiffusionVec.jl")
# include("DyCore/Jac.jl")
include("DyCore/JacSchur.jl")
include("DyCore/MassCG.jl")
include("DyCore/NumberingFemCG.jl")
# include("DyCore/PlotCG.jl")
# include("DyCore/PlotDG.jl")
include("DyCore/Project.jl")
include("DyCore/ProjectVec.jl")
include("DyCore/Source.jl")
# include("DyCore/SphereGrid.jl")
include("DyCore/simpson.jl")
include("DyCore/vtkCG.jl")
include("DyCore/vtkCGGrid.jl")
include("DyCore/vtkOutput.jl")
# include("IntegrationMethods/JacVec.jl")
# include("IntegrationMethods/MethodRKIMEX.jl")
# include("IntegrationMethods/MethodRos.jl")
# include("IntegrationMethods/Rosenbrock.jl")
include("IntegrationMethods/RosenbrockMethod.jl")
include("IntegrationMethods/RosenbrockSchur.jl")
# include("IntegrationMethods/RosenbrockSchur1.jl")
include("IntegrationMethods/RungeKuttaExplicit.jl")
# include("IntegrationMethods/RungeKuttaHEVI.jl")
include("IntegrationMethods/RungeKuttaMethod.jl")
include("IntegrationMethods/SchurSolve.jl")
include("IntegrationMethods/JacStruc.jl")
# include("Model/EnergyToPot.jl")
# include("Model/EtaFromZ.jl")
# include("Model/PhiBaroWave.jl")
include("Model/PhysParameters.jl")
# include("Model/PotToEnergy.jl")
# include("Model/Pres.jl")
include("Model/Pressure.jl")
# include("Model/SymJac.jl")
# include("Model/TBaroWave.jl")
# include("Model/TFromThetaRho.jl")
include("Model/dPresdTh.jl")
include("Model/fRho.jl")
include("Model/fRhoBGrd.jl")
include("Model/fT.jl")
include("Model/fTBGrd.jl")
include("Model/fTheta.jl")
include("Model/fThetaBGrd.jl")
# include("Model/fTr.jl")
include("Model/fVel.jl")
# include("Model/fVelW.jl")
include("Model/fpBGrd.jl")
# include("Model/pFromTheta.jl")

end # module


text:
using Revise, Pkg

#if  pwd() ∉ LOAD_PATH
#  push!(LOAD_PATH,pwd())
#end

Pkg.develop(PackageSpec(path=pwd()))

import QFactors

#Pkg.activate(pwd())

#revise(QFactors)

#Pkg.develop(PackageSpec(path=pwd()))

#QFactors.recompressqfactorfiles()
#QFactors.recompressqfactorfile("COMP-A")

@time begin
  QFactors.prepcrsp(refreshcrsp=false)
  #QFactors.prepcomp(refreshcomp=false)
end

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
function calc_W_from_B(Tk, B, P, EPS=1e-8, MAXITER=100)

    pws = Pws(B)
    efac = efactor(B,P)
    xsv = efac * pws / P
    w2 = Mv/Ma * xsv/(1-xsv)

    w = (enthalpyair(B,P) - enthalpyair(Tk,P) - w2*(enthalpywi(B) - enthalpyvapor(B))) /
        (enthalpyvapor(Tk) - enthalpywi(B))

    for iter = 1:MAXITER
        f = aux_WB(w, Tk, B, P, pws, efac)
        df = (aux_WB(w+1e-4*w2, Tk, B, P, pws, efac) - f) / (1e-4*w2)
        dw = -f/df
        w = w + dw

        if abs(dw) < EPS*w2
            return w
        end

    end
    return w
end

function aux_WB(w, Tk, B, P, pws, efac)
    xv1 = w / (Mv/Ma+w)
    xv2 = efac * pws / P
    w2 = Mv / Ma * xv2 / (1.0-xv2)
    #(1.0+w)*enthalpymoist(Tk,P,xv1) + (w2-w)*enthalpywi(B) - (1.0+w2)*enthalpymoist(B,P,xv2)
    enthalpymoist(Tk,P,xv1) + (w2-w)*enthalpywi(B) - enthalpymoist(B,P,xv2)
end

function aux_WB(w, Tk, B, P)
    pws = Pws(B)

    efac = efactor(B,P)

    xv1 = w / (Mv/Ma+w)
    xv2 = efac * pws / P
    w2 = Mv / Ma * xv2 / (1.0-xv2)
    #(1.0+w)*enthalpymoist(Tk,P,xv1) + (w2-w)*enthalpywi(B) - (1.0+w2)*enthalpymoist(B,P,xv2)
    enthalpymoist(Tk,P,xv1) + (w2-w)*enthalpywi(B) - enthalpymoist(B,P,xv2)
end

function calcwetbulb(Tk, P, xv, EPS=1e-8, MAXITER=200)

    w = humrat(xv)

    B = Tk - 1.0 # Initial guess
    h = 1e-7
    for i = 1:MAXITER
        f = aux_WB(w, Tk, B, P)
        df = (aux_WB(w, Tk, B + h, P) - f) / h
        dB = -f/df
        B = B + dB

        if abs(dB) < EPS
            return B
        end
    end
    throw(ConvergenceError("Wet Bulb temperature failed to converge!", B, MAXITER, dB))

end


text:
module ZooTests

using Photon, Test
using Photon.Zoo: VGG16, AlexNet, DenseNet169, DenseNet121

# include("../src/zoo/zoo.jl")
# using .Zoo

"""
############################################
#### Start of Testing of Specific Models ###
############################################
"""

getimages(s=224) = KorA(randn(Float32,s,s,3,4))

function test_vgg()
    model = VGG16()
    images = getimages()
    pred = model(images)
    @test size(pred) == (1000,size(images,4))
end

@testset "VGG" begin
    test_vgg()
end

function test_alexnet()
    model = AlexNet(classes=500)
    images = getimages()
    pred = model(images)
    @test size(pred) == (500,size(images,4))
end

@testset "AlexNet" begin
    test_alexnet()
end

function test_densenet()
    model = DenseNet169()
    images = getimages()
    pred = model(images)
    @test size(pred) == (1000, size(images,4))

    model = DenseNet121(classes=500)
    images = getimages()
    pred = model(images)
    @test size(pred) == (500, size(images,4))

end

@testset "DenseNet" begin
    test_densenet()
end

end


text:
@testset "Exceptions" begin
  # number of simulations is negative
  @test_throws ErrorException SimpleBirthDeathProcess.rand_continuous(-1, 1, 1, [1, 1])
  @test_throws ErrorException SimpleBirthDeathProcess.rand_discrete(-1, 1, 1, 1, [1, 1])

  # number of simulations is zero
  @test_throws ErrorException SimpleBirthDeathProcess.rand_continuous(0, 1, 1, [1, 1])
  @test_throws ErrorException SimpleBirthDeathProcess.rand_discrete(0, 1, 1, 1, [1, 1])

  # initial population size is zero or negative
  @test_throws ErrorException SimpleBirthDeathProcess.rand_continuous(1, -1, 1, [1, 1])
  @test_throws ErrorException SimpleBirthDeathProcess.rand_continuous(1,  0, 1, [1, 1])

  @test_throws ErrorException SimpleBirthDeathProcess.rand_continuous(-1, 1, [1, 1])
  @test_throws ErrorException SimpleBirthDeathProcess.rand_continuous( 0, 1, [1, 1])

  @test_throws ErrorException SimpleBirthDeathProcess.rand_discrete(1, -1, 1, 1, [1, 1])
  @test_throws ErrorException SimpleBirthDeathProcess.rand_discrete(1,  0, 1, 1, [1, 1])

  @test_throws ErrorException SimpleBirthDeathProcess.rand_discrete(-1, 1, 1, [1, 1])
  @test_throws ErrorException SimpleBirthDeathProcess.rand_discrete( 0, 1, 1, [1, 1])

  # observation time is negative (continuous)
  @test_throws ErrorException SimpleBirthDeathProcess.rand_continuous(1, 1, -1, [1, 1])
  @test_throws ErrorException SimpleBirthDeathProcess.rand_continuous(1, -1, [1, 1])

  # observation time is negative (discrete)
  @test_throws ErrorException SimpleBirthDeathProcess.rand_discrete(1, 1, -1, 1, [1, 1])
  @test_throws ErrorException SimpleBirthDeathProcess.rand_discrete(1, -1, 1, [1, 1])

  # time lag is negative or zero (discrete)
  @test_throws ErrorException SimpleBirthDeathProcess.rand_discrete(1, 1, 1, -1, [1, 1])
  @test_throws ErrorException SimpleBirthDeathProcess.rand_discrete(1, 1, -1, [1, 1])

  @test_throws ErrorException SimpleBirthDeathProcess.rand_discrete(1, 1, 1, 0, [1, 1])
  @test_throws ErrorException SimpleBirthDeathProcess.rand_discrete(1, 1, 0, [1, 1])

  # birth rate is negative
  @test_throws ErrorException SimpleBirthDeathProcess.rand_continuous(1, 1, 1, [-1, 1])
  @test_throws ErrorException SimpleBirthDeathProcess.rand_continuous(1, 1, [-1, 1])

  @test_throws ErrorException SimpleBirthDeathProcess.rand_discrete(1, 1, 1, 1, [-1, 1])
  @test_throws ErrorException SimpleBirthDeathProcess.rand_discrete(1, 1, 1, [-1, 1])

  # death rate is negative
  @test_throws ErrorException SimpleBirthDeathProcess.rand_continuous(1, 1, 1, [1, -1])
  @test_throws ErrorException SimpleBirthDeathProcess.rand_continuous(1, 1, [1, -1])

  @test_throws ErrorException SimpleBirthDeathProcess.rand_discrete(1, 1, 1, 1, [1, -1])
  @test_throws ErrorException SimpleBirthDeathProcess.rand_discrete(1, 1, 1, [1, -1])

  # birth rate and death rate both equal to zero
  @test_throws ErrorException SimpleBirthDeathProcess.rand_continuous(1, 1, 1, [0, 0])
  @test_throws ErrorException SimpleBirthDeathProcess.rand_continuous(1, 1, [0, 0])

  @test_throws ErrorException SimpleBirthDeathProcess.rand_discrete(1, 1, 1, 1, [0, 0])
  @test_throws ErrorException SimpleBirthDeathProcess.rand_discrete(1, 1, 1, [0, 0])
end

@testset "Type stability" begin
  @testset for I = [Int8, Int16, Int32, Int64, Int128], F = [Float32, Float64, BigFloat]
    # all operations happen with the same type passed to the function...
    @test isa(SimpleBirthDeathProcess.rand_continuous(I(2), I(100), I(1), [F(1), F(1)]), Vector{SimpleBirthDeathProcess.ObservationContinuousTime{F}})
    @test isa(SimpleBirthDeathProcess.rand_continuous(I(2), I(100), F(1), [F(1), F(1)]), Vector{SimpleBirthDeathProcess.ObservationContinuousTime{F}})
    @test isa(SimpleBirthDeathProcess.rand_continuous(I(100), I(1), [F(1), F(1)]), SimpleBirthDeathProcess.ObservationContinuousTime{F})
    @test isa(SimpleBirthDeathProcess.rand_continuous(I(100), F(1), [F(1), F(1)]), SimpleBirthDeathProcess.ObservationContinuousTime{F})

    @test isa(SimpleBirthDeathProcess.rand_discrete(I(2), I(100), I(1), I(1), [F(1), F(1)]), SimpleBirthDeathProcess.ObservationDiscreteTimeEven{F, I})
    @test isa(SimpleBirthDeathProcess.rand_discrete(I(2), I(100), I(1), F(1), [F(1), F(1)]), SimpleBirthDeathProcess.ObservationDiscreteTimeEven{F, I})
    @test isa(SimpleBirthDeathProcess.rand_discrete(I(100), I(1), I(1), [F(1), F(1)]), SimpleBirthDeathProcess.ObservationDiscreteTimeEven{F, I})
    @test isa(SimpleBirthDeathProcess.rand_discrete(I(100), I(1), F(1), [F(1), F(1)]), SimpleBirthDeathProcess.ObservationDiscreteTimeEven{F, I})

    # ...but if parameters are passed as integers, operations use Float64
    @test isa(SimpleBirthDeathProcess.rand_continuous(I(2), I(100), I(1), [I(1), I(1)]), Vector{SimpleBirthDeathProcess.ObservationContinuousTime{Float64}})
    @test isa(SimpleBirthDeathProcess.rand_continuous(I(2), I(100), I(1), [F(1), I(1)]), Vector{SimpleBirthDeathProcess.ObservationContinuousTime{F}})
    @test isa(SimpleBirthDeathProcess.rand_continuous(I(2), I(100), I(1), [I(1), F(1)]), Vector{SimpleBirthDeathProcess.ObservationContinuousTime{F}})
    @test isa(SimpleBirthDeathProcess.rand_continuous(I(100), I(1), [I(1), I(1)]), SimpleBirthDeathProcess.ObservationContinuousTime{Float64})
    @test isa(SimpleBirthDeathProcess.rand_continuous(I(100), I(1), [F(1), I(1)]), SimpleBirthDeathProcess.ObservationContinuousTime{F})
    @test isa(SimpleBirthDeathProcess.rand_continuous(I(100), I(1), [I(1), F(1)]), SimpleBirthDeathProcess.ObservationContinuousTime{F})

    # simulations happen as Float64 but we only see the realizations as integer counts
    @test isa(SimpleBirthDeathProcess.rand_discrete(I(2), I(100), I(1), I(1), [I(1), I(1)]), SimpleBirthDeathProcess.ObservationDiscreteTimeEven{Float64, I})
    @test isa(SimpleBirthDeathProcess.rand_discrete(I(100), I(1), I(1), [F(1), I(1)]), SimpleBirthDeathProcess.ObservationDiscreteTimeEven{F, I})
    @test isa(SimpleBirthDeathProcess.rand_discrete(I(100), I(1), I(1), [I(1), F(1)]), SimpleBirthDeathProcess.ObservationDiscreteTimeEven{F, I})

    @test isa(SimpleBirthDeathProcess.rand_discrete(I(100), I(1), F(1), [I(1), I(1)]), SimpleBirthDeathProcess.ObservationDiscreteTimeEven{F, I})
    @test isa(SimpleBirthDeathProcess.rand_discrete(I(100), I(1), F(1), [F(1), I(1)]), SimpleBirthDeathProcess.ObservationDiscreteTimeEven{F, I})
    @test isa(SimpleBirthDeathProcess.rand_discrete(I(100), I(1), F(1), [I(1), F(1)]), SimpleBirthDeathProcess.ObservationDiscreteTimeEven{F, I})
  end
end

@testset "Observation time is zero" begin
  I = Int32
  F = Float32

  x = SimpleBirthDeathProcess.rand_continuous(I(2), I(100), I(0), [I(1), I(1)])
  @test isa(x, Vector{SimpleBirthDeathProcess.ObservationContinuousTime{Float64}})
  @test length(x) == 2
  @test x[1].tot_births === zero(Int)
  @test x[1].tot_deaths === zero(Int)
  @test x[1].integrated_jump === zero(Float64)
  @test x[1].sum_log_n === zero(Float64)
  @test x[1].len === zero(Int)
  @test x[1].waiting_time == Float64[]
  @test x[1].initial_population_size == Int(100)
  @test x[1].increment == Int[]
  @test x[2].tot_births === zero(Int)
  @test x[2].tot_deaths === zero(Int)
  @test x[2].integrated_jump === zero(Float64)
  @test x[2].sum_log_n === zero(Float64)
  @test x[2].len === zero(Int)
  @test x[2].waiting_time == Float64[]
  @test x[2].initial_population_size == Int(100)
  @test x[2].increment == Int[]

  x = SimpleBirthDeathProcess.rand_continuous(I(2), I(100), F(0), [F(1), F(1)])
  @test isa(x, Vector{SimpleBirthDeathProcess.ObservationContinuousTime{F}})
  @test length(x) == 2
  @test x[1].tot_births === zero(Int)
  @test x[1].tot_deaths === zero(Int)
  @test x[1].integrated_jump === zero(F)
  @test x[1].sum_log_n === zero(Float64)
  @test x[1].len === zero(Int)
  @test x[1].waiting_time == F[]
  @test x[1].initial_population_size == Int(100)
  @test x[1].increment == Int[]
  @test x[2].tot_births === zero(Int)
  @test x[2].tot_deaths === zero(Int)
  @test x[2].integrated_jump === zero(F)
  @test x[2].sum_log_n === zero(Float64)
  @test x[2].len === zero(Int)
  @test x[2].waiting_time == F[]
  @test x[2].initial_population_size == Int(100)
  @test x[2].increment == Int[]

  x = SimpleBirthDeathProcess.rand_continuous(I(100), I(0), [I(1), I(1)])
  @test isa(x, SimpleBirthDeathProcess.ObservationContinuousTime{Float64})
  @test x.tot_births === zero(Int)
  @test x.tot_deaths === zero(Int)
  @test x.integrated_jump === zero(Float64)
  @test x.sum_log_n === zero(Float64)
  @test x.len === zero(Int)
  @test x.waiting_time == Float64[]
  @test x.initial_population_size == Int(100)
  @test x.increment == Int[]

  x = SimpleBirthDeathProcess.rand_continuous(I(100), F(0), [F(1), F(1)])
  @test isa(x, SimpleBirthDeathProcess.ObservationContinuousTime{F})
  @test x.tot_births === zero(Int)
  @test x.tot_deaths === zero(Int)
  @test x.integrated_jump === zero(F)
  @test x.sum_log_n === zero(Float64)
  @test x.len === zero(Int)
  @test x.waiting_time == F[]
  @test x.initial_population_size == Int(100)
  @test x.increment == Int[]

  x = SimpleBirthDeathProcess.rand_discrete(I(2), I(100), I(0), I(1), [I(1), I(1)])
  @test isa(x, SimpleBirthDeathProcess.ObservationDiscreteTimeEven{Float64, I})
  @test x.n == Int(2)
  @test x.k == zero(Int)
  @test x.u == one(Float64)
  @test size(x.state) == (1, 2)
  @test x.state[1, 1] == I(100)
  @test x.state[1, 2] == I(100)

  x = SimpleBirthDeathProcess.rand_discrete(I(2), I(100), I(0), F(1), [F(1), F(1)])
  @test isa(x, SimpleBirthDeathProcess.ObservationDiscreteTimeEven{F, I})
  @test x.n == Int(2)
  @test x.k == zero(Int)
  @test x.u == one(F)
  @test size(x.state) == (1, 2)
  @test x.state[1, 1] == I(100)
  @test x.state[1, 2] == I(100)

  x = SimpleBirthDeathProcess.rand_discrete(I(100), I(0), I(1), [I(1), I(1)])
  @test isa(x, SimpleBirthDeathProcess.ObservationDiscreteTimeEven{Float64, I})
  @test x.n == one(Int)
  @test x.k == zero(Int)
  @test x.u == one(Float64)
  @test size(x.state) == (1, 1)
  @test x.state[1, 1] == I(100)

  x = SimpleBirthDeathProcess.rand_discrete(I(100), I(0), F(1), [F(1), F(1)])
  @test isa(x, SimpleBirthDeathProcess.ObservationDiscreteTimeEven{F, I})
  @test x.n == one(Int)
  @test x.k == zero(Int)
  @test x.u == one(F)
  @test size(x.state) == (1, 1)
  @test x.state[1, 1] == I(100)
end

@testset "Correct output" begin
  # we check if the expected value and variance are close to the theoretical one
  n = 100_000
  i = 10_000
  t = 10.0
  k = 10
  u = 1.0
  λ = 0.001
  μ = 0.002

  θ = (λ - μ) * t
  e = exp(θ)
  v = (λ + μ) * e * expm1(θ) / (i * (λ - μ))

  x = SimpleBirthDeathProcess.rand_continuous(n, i, t, [λ, μ])
  @test isa(x, Vector{SimpleBirthDeathProcess.ObservationContinuousTime{Float64}})
  @test abs(mean([(1 + sum(x[s].increment) / i) for s = 1:n]) - e) < 1e-5
  @test abs(var([(1 + sum(x[s].increment) / i) for s = 1:n]) - v) < 1e-5

  x = SimpleBirthDeathProcess.rand_discrete(n, i, k, u, [λ, μ])
  @test isa(x, SimpleBirthDeathProcess.ObservationDiscreteTimeEven{Float64, Int64})
  @test size(x.state) == (k + 1, n)
  @test abs(mean(x.state[end, :] ./ i) - e) < 1e-5
  @test abs(var(x.state[end, :] ./ i) - v) < 1e-5
end


text:
#=--- TerminalGraphics / deps / build.jl ----------------------------------=#

#=-------------------------------------------------------------------------=#
#
#   Download an build a recent version of 'libsixel'
#
#=-------------------------------------------------------------------------=#

using BinDeps

validate_libsixel(name, handle) = begin
    # try to locate sixel_encoder_encode_bytes, which is only
    # part of the newest, yet un-'released' version, but strictly necessary.
    println("Checking library $name ($handle)")
    try
        Base.Libdl.dlsym(handle, :sixel_encoder_encode_bytes)
        return true
    catch
        #
    end
    false
end

@BinDeps.setup

libsixel = library_dependency( "libsixel",
                aliases = ["libsixel.so.1.0.6","libsixel.so.1","libsixel.so"],
                validate = validate_libsixel
                )

# build from source
provides(Sources, URI("https://github.com/saitoha/libsixel/archive/master.zip")
                , libsixel
                , unpacked_dir = "libsixel-master")
provides(BuildProcess,
                Autotools(libtarget="src/libsixel.la")
                , libsixel)

@BinDeps.install Dict(:libsixel => :libsixel)

#=--- end of file ---------------------------------------------------------=#


text:
# data/z-test.jl

using Test

@test ir_load_brainweb_t1_256() !== nothing


text:
@testset "transformations" begin

    p1 = Point([-75.69926351308823,45.43145021122502])
    poly = Polygon([[[0, 29], [3.5, 29], [2.5, 32], [0, 29]]])

    @test transform_rotate(geojson=p1, angle=80., pivot=Point([-75.6, 45.3])).coordinates ≈ [-75.433041, 45.391507]
    @test transform_rotate(geojson=poly, angle=200., pivot=Point([2.48291015625, 27.819644755099446])).coordinates ≈
        [[[4.331318,25.963562],[1.0811,27.016436],[0.878141,23.89646],[4.331318,25.963562]]] atol = 0.0001

    coll = FeatureCollection([Feature(Point([-75.69926351308823,45.43145021122502])),
            Feature(Polygon([[[0, 29], [3.5, 29], [2.5, 32], [0, 29]]]))])

	@test transform_scale(coll, 0.1).features[2].geometry.coordinates ≈ [[[
									1.3495,
									29.675
								],
								[
									1.700667,
									29.675
								],
								[
									1.598959,
									29.975
								],
								[
									1.3495,
									29.675]]]

    @test transform_scale(coll, 5., "sw").features[1].geometry.coordinates ≈ [-75.699264,45.43145]
    @test_throws ErrorException transform_scale(coll, 5., "qw")

	poly2 = Polygon([[[100, 0], [101, 0], [101, 1], [100, 1], [100, 0]]])
	feats = explode(poly2).features
	res = [Point([100, 0]), Point([101, 0]), Point([101, 1]), Point([100, 1]), Point([100, 0])]

	for i in eachindex(feats)
		@test feats[i].geometry.coordinates == res[i].coordinates
	end

	@test flip(Point([77.34374999999999,43.58039085560784,3000])).coordinates == [43.58039085560784,77.34374999999999,3000]

	p2 = Point([92.46093749999999,54.67383096593114])
	poly2 = Polygon([[[48.1640625,20.632784250388028],[76.640625,20.632784250388028],
            [
              76.640625,
              38.8225909761771
            ],
            [
              48.1640625,
              38.8225909761771
            ],
            [
              48.1640625,
              20.632784250388028]]])

	@test polygon_tangents(p2, poly2).features[1].geometry.coordinates ≈ [48.1640625,38.8225909761771]
	@test polygon_tangents(p2, poly2).features[2].geometry.coordinates ≈ [76.640625, 20.632784250388028]

	poly3 = Polygon([[[-2.275543, 53.464547 ],
		[-2.275543, 53.489271 ],
		[-2.215118, 53.489271 ],
		[-2.215118, 53.464547 ],
		[-2.275543, 53.464547 ]]])

	@test polygon_to_line(poly3).coordinates ≈ [
			[
				-2.275543,
				53.464547
			],
			[
				-2.275543,
				53.489271
			],
			[
				-2.215118,
				53.489271
			],
			[
				-2.215118,
				53.464547
			],
			[
				-2.275543,
				53.464547
			]]

	poly5 = Polygon([[[0, 29], [3.5, 29], [2.5, 32], [0, 29]]])
	poly6 = Polygon([[
        [
          121.46484375,
          -23.88583769986199,
          10
        ],
        [
          125.5078125,
          -27.488781168937983,
          11
        ],
        [
          129.638671875,
          -25.562265014427492,
          12
        ],
        [
          129.2431640625,
          -22.06527806776582,
          13
        ],
        [
          124.93652343749999,
          -20.385825381874263,
          14
        ],
        [
          121.46484375,
          -23.88583769986199,
          10
        ]]])

	@test transform_translate(poly5, 300, 70).coordinates ≈ [
					[
						[
							2.911836,
							29.922757
						],
						[
							6.411836,
							29.922757
						],
						[
							5.504792,
							32.922757
						],
						[
							2.911836,
							29.922757
						]
					]]

	@test transform_translate(poly6, 120, 165, 15, false, "miles").coordinates ≈ [
					[
						[
							121.959747,
							-25.563437,
							25
						],
						[
							126.018505,
							-29.166381,
							26
						],
						[
							130.140552,
							-27.239865,
							27
						],
						[
							129.731173,
							-23.742878,
							28
						],
						[
							125.418767,
							-22.063425,
							29
						],
						[
							121.959747,
							-25.563437,
							25
						]
					]]

	poly7 = Polygon([[[
					-20031393.380526427,
					-1860477.2684611906
				],
				[
					19998984.080533516,
					-1906645.233545436
				],
				[
					20014882.98241683,
					-1959233.9090056375
				],
				[
					-20002653.0578912,
					-1937831.5410857894
				],
				[
					-19990117.385252435,
					-1892275.0722278224
				],
				[
					-20031393.380526427,
					-1860477.2684611906
					]]])

	#@test convert_to(poly7, "wgs84").coordinates

	res = convert_to(Point([10, 40]), "mercator")
	@test convert_to(res, "wgs84") ≈ Point([10, 40]) atol=0.5

	pt1 = Point([50, 51])
	pt2 = MultiPoint([[100, 101], [101, 102]])

	fc1 = combine(FeatureCollection([Feature(pt1), Feature(pt2)]))

	@test fc1.features[1].geometry.coordinates ≈ [[50, 51], [100, 101], [101, 102]]

	l1 = LineString([[102.0,-10.0],[130.0,4.0]])
	l2 = LineString([[40.0,-20.0],[150.0,18.0]])

	fc2 = combine(FeatureCollection([Feature(l1), Feature(l2)]))

	@test fc2.features[1].geometry.coordinates == [[[102, -10], [130, 4], [40, -20], [150, 18]]]

	pl1 = Polygon([[
            [20.0, 0.0],
            [101.0, 0.0],
            [101.0, 1.0],
            [100.0, 1.0],
            [100.0, 0.0],
            [20.0, 0.0]]])

	pl2 = MultiPolygon([
		[[
            [30.0, 0.0],
            [102.0, 0.0],
            [103.0, 1.0],
            [30.0, 0.0]
        ]],
        [[
                [20.0, 5.0],
                [92.0, 5.0],
                [93.0, 6.0],
                [20.0, 5.0]
            ],
            [
                [25, 5],
                [30, 5],
                [30, 5.5],
                [25, 5]]]])

	fc3 = combine(FeatureCollection([Feature(pl1), Feature(pl2)]))
	@test fc3.features[1].geometry.coordinates == [[[[20, 0], [101.0, 0.0], [101.0, 1.0], [100.0, 1.0], [100.0, 0.0], [20, 0],[30.0, 0.0], [102.0, 0.0], [103.0, 1.0], [30.0, 0.0]]]]

	l3 = LineString([[102.0, -10.0],[130.0, 4.0]])
	ml = MultiLineString([[[40.0, -20.0],[150.0, 18.0], [50, -10],[160, 28]]])

	fc4 = combine(FeatureCollection([Feature(l3), Feature(ml)]))
	@test fc4.features[1].geometry.coordinates == [[[102., -10.], [130., 4.], [40., -20.], [150., 18.], [50, -10], [160, 28]]]

	fc4 = FeatureCollection([Feature(pt1)])
	fc5 = FeatureCollection([Feature(Polygon([[[49, 52], [51, 50.5], [49, 51], [49, 52], [49, 52]]]), Dict("poly" => true))])
	@test tag(fc4, fc5,"poly", "inside").features[1].properties == Dict("inside" => true)

	fcoll = GeoJSON.parsefile(string(@__DIR__) * "/simplify.geojson")
	fcoll2 = GeoJSON.parsefile(string(@__DIR__) * "/simplify2.geojson")

	@test simplify(fcoll, 100).features[1].geometry.coordinates ≈ [
			[
				[
					26.14843,
					-28.297552
				],
				[
					26.150354,
					-28.302606
				],
				[
					26.135463,
					-28.304283
				],
				[
					26.14843,
					-28.297552
				]
			]
		]

	@test_throws ErrorException simplify(fcoll2, -1)
	@test simplify(fcoll2, 0.005, true).features[1].geometry.coordinates ≈ [
			[
				[
					179.975281,
					-16.51477
				],
				[
					179.980431,
					-16.539127
				],
				[
					180.0103,
					-16.523328
				],
				[
					180.007553,
					-16.534848
				],
				[
					180.018196,
					-16.539127
				],
				[
					180.061455,
					-16.525632
				],
				[
					180.066605,
					-16.513124
				],
				[
					180.046349,
					-16.479547
				],
				[
					180.086861,
					-16.44761
				],
				[
					180.084114,
					-16.441354
				],
				[
					180.055618,
					-16.439707
				],
				[
					180.026093,
					-16.464732
				],
				[
					180.01442,
					-16.464073
				],
				[
					179.975281,
					-16.51477
				]
			]
		]
end


text:
function testidea1()

    rg = MersenneTwister(1)

    # define parameter distributions

    mu    = [3.57; 9.95; 3.87; 9.95; 28.22; 14.07]

    Sigma = [0.0148; 0.0706; 0.0187; 0.0266; 0.1006; 0.0238]

    # Since Sigma is a diagonal covariance matrix
    # we can sample parameters a and b independently
    # (in the future these parameters will be correlated)

    pdf_a = MvNormal(mu[1:3], Diagonal(Sigma[1:3]))

    pdf_b = MvNormal(mu[4:6], Diagonal(Sigma[4:6]))

    # Distribution of normed vectors (this needs to be refined)
    A = reduce(hcat, [(a=rand(rg,pdf_a); a/norm(a)) for i=1:10_000])'

    pdf_a_norm = MvNormal(vec(mean(A,dims=1)), cov(A))

    # draw samples for b
    B = [rand(pdf_b) for i=1:200]

    figure(1)
    cla()
    xrange = collect(LinRange(-2,2,100))

    for i=1:5

        a, b = rand(rg, pdf_a), rand(rg, pdf_b)

        for x in xrange
               local p = b + a*x
               plot3D(p[1], p[2], p[3], "r.", alpha=0.3)
        end
    end

    function test(x)

        l = zeros(length(B))

        for (i,b) in enumerate(B)

            local v = x - b

            l[i] = logpdf(pdf_a_norm, v / norm(v))

        end

        return logsumexp(l) - log(length(B))

    end

end


text:
module BasisSetModule
using StringParserPEG
export BasisSet, readBasisSetTX93
using ..BaseModule
using ..AtomModule
import Base.print
import Base.show
import Base.==
import Base.≈
using Printf

struct PrimitiveGaussianBasisFunctionDefinition
    exponent::Float64
    prefactor::Float64
end
print(io::IO,pgbf::PrimitiveGaussianBasisFunctionDefinition) = print(io,@sprintf("lin:%9.5f exp:%12.5f",pgbf.prefactor, pgbf.exponent))
show(io::IO,::MIME"text/plain",pgbf::PrimitiveGaussianBasisFunctionDefinition) = print(io,pgbf)
≈(pgbfdef1::PrimitiveGaussianBasisFunctionDefinition,pgbfdef2::PrimitiveGaussianBasisFunctionDefinition;kwargs...) = (≈(pgbfdef1.exponent,pgbfdef2.exponent;kwargs...) && ≈(pgbfdef1.prefactor,pgbfdef2.prefactor;kwargs...))
≈(vec1::Vector{PrimitiveGaussianBasisFunctionDefinition},vec2::Vector{PrimitiveGaussianBasisFunctionDefinition};kwargs...) = all(.≈(vec1,vec2;kwargs...))

struct ContractedGaussianBasisFunctionDefinition
    lQuantumNumber::LQuantumNumber
    primitives::Array{PrimitiveGaussianBasisFunctionDefinition,1}
end
function show(io::IO,::MIME"text/plain",cgbf::ContractedGaussianBasisFunctionDefinition)
  str  = cgbf.lQuantumNumber.symbol*"─ $(cgbf.primitives[1])\n"
  for prim in cgbf.primitives[2:end-1]
    str *= "│  $(prim)\n"
  end
  str *= "╰─ $(cgbf.primitives[end])\n"
  print(io,str)
end
function show(io::IO,mime::MIME"text/plain",vec::Vector{ContractedGaussianBasisFunctionDefinition})
  for cgbf in vec
    show(io,mime,cgbf)
  end
end
==(cgbfdef1::ContractedGaussianBasisFunctionDefinition,cgbfdef2::ContractedGaussianBasisFunctionDefinition) = (cgbfdef1.lQuantumNumber==cgbfdef2.lQuantumNumber && cgbfdef1.primitives == cgbfdef2.primitives)
==(vec1::Vector{ContractedGaussianBasisFunctionDefinition},vec2::Vector{ContractedGaussianBasisFunctionDefinition}) = all(vec1 .== vec2)
≈(cgbfdef1::ContractedGaussianBasisFunctionDefinition,cgbfdef2::ContractedGaussianBasisFunctionDefinition;kwargs...) = (cgbfdef1.lQuantumNumber==cgbfdef2.lQuantumNumber && all(.≈(cgbfdef1.primitives, cgbfdef2.primitives; kwargs...)))
≈(vec1::Vector{ContractedGaussianBasisFunctionDefinition},vec2::Vector{ContractedGaussianBasisFunctionDefinition};kwargs...) = all(.≈(vec1, vec2;kwargs...))

struct BasisSet
    definitions::Dict{Element,Array{ContractedGaussianBasisFunctionDefinition,1}}
end
function show(io::IO,mime::MIME"text/plain",basset::BasisSet)
  for (key,val) in basset.definitions
    println(io,"$key => ")
    show(io,mime,val)
  end
end
≈(bassset1::BasisSet,basset2::BasisSet;kwargs...) = keys(basset1.definitions)==keys(basset2.definitions) && all(≈(basset1.defintions[el],basset2.definitons[el];kwargs...) for el in keys(basset1.definitions))

tx93 = Grammar("""
	       start   => (basset & -(*(comment))) {liftchild}
	       basset  => *((-(*(comment)) & elementgroup) {liftchild}) {"basset"}
	       comment => ?(r(^!.*)r | space) & eol

	       elementgroup  => (elementheader & CGBFs) {"entry"}
	       elementheader => (-('FOR' & space) & element & -(?(space) & eol)) {liftchild}
	       element       => r([A-Za-z]*)r {"element"}
	       CGBFs	       => +(CGBF) {"cgbfs"}

	       lqn => ('SPD' | 'SP' | 'L' | 'S' | 'P' | 'D' | 'F' | 'G' | 'H') {nodevalue}
	       CGBF    => (lqn & +(PGBF)) {"cgbf"}
	       PGBF    => (-(space) & float & +(-(space) & float) & -(?(space) & eol)) {"pgbf"}
	       """,standardrules)

fromTX93(node,children,::MatchRule{:pgbf}) = [children[1],[n.children[1] for n in children[2].children]]
function QLcgbf(lins,exps,lqn,index=1)
  ContractedGaussianBasisFunctionDefinition(
     LQuantumNumber(lqn),
     [PrimitiveGaussianBasisFunctionDefinition(exp,pre[index]) for (pre,exp) in zip(lins,exps)]
  )
end
function fromTX93(node,children,::MatchRule{:cgbf})
  lqn  = children[1]
  exps = [p[1] for p in children[2].children]
  lins = [p[2] for p in children[2].children]
  if lqn == "SPD"
    @assert length(lins[1]) == 3
    sshell = QLcgbf(lins,exps,"S",1)
    pshell = QLcgbf(lins,exps,"P",2)
    dshell = QLcgbf(lins,exps,"D",3)
    return [sshell,pshell,dshell]
  elseif lqn == "SP" || lqn == "L"
    @assert length(lins[1]) == 2
    sshell = QLcgbf(lins,exps,"S",1)
    pshell = QLcgbf(lins,exps,"P",2)
    return [sshell,pshell]
  elseif lqn in ["S","P","D","F","G","H"]
    @assert length(lins[1]) == 1
    shell = QLcgbf(lins,exps,lqn,1)
    return [shell]
  else
    @error("Unexpected sharedsp lqn ($lqn)")
  end
end
fromTX93(node,children,::MatchRule{:cgbfs}) = vcat(children...)
fromTX93(node,children,::MatchRule{:element}) = Element(node.value)
fromTX93(node,children,::MatchRule{:entry}) =(children[1] => children[2])
fromTX93(node,children,::MatchRule{:basset}) = BasisSet(Dict(children...))

function BasisSetTX93(str::AbstractString)
  (ast,pos,err) = parse(tx93,str)
  if err != nothing
    @error err
  end
  return transform(fromTX93,ast)
end

function readBasisSetTX93(filename::AbstractString)
  fd = open(filename)
  basSet = BasisSetTX93(read(fd,String))
  close(fd)
  return basSet
end

end # module


text:
using Test
using BinaryBuilderBase
using BinaryBuilderBase: platform_dlext, platform_exeext

@testset "Wrappers utilities" begin
    @test nbits(Platform("i686", "linux")) == 32
    @test nbits(Platform("x86_64", "linux"; march="avx")) == 64
    @test nbits(Platform("armv7l", "linux")) == 32
    @test nbits(Platform("aarch64", "linux"; cuda="10.1")) == 64
    @test nbits(Platform("powerpc64le", "linux")) == 64
    @test nbits(AnyPlatform()) == 64

    @test proc_family(Platform("i686", "linux")) == "intel"
    @test proc_family(Platform("x86_64", "linux"; march="avx")) == "intel"
    @test proc_family(Platform("armv7l", "linux")) == "arm"
    @test proc_family(Platform("aarch64", "linux"; cuda="10.1")) == "arm"
    @test proc_family(Platform("powerpc64le", "linux")) == "power"
    @test proc_family(AnyPlatform()) == "any"

    @test platform_exeext(Platform("i686", "linux")) == ""
    @test platform_exeext(Platform("x86_64", "freebsd")) == ""
    @test platform_exeext(Platform("x86_64", "macos")) == ""
    @test platform_exeext(Platform("i686", "windows")) == ".exe"
    @test platform_exeext(Platform("x86_64", "linux"; march="avx512")) == ""

    @test aatriplet(Platform("x86_64", "linux"; libc="glibc", libgfortran_version=v"4.0.0", march="avx", cuda="9.2")) == "x86_64-linux-gnu"
end

# Are we using docker? If so, test that the docker runner works...
@testset "Runner utilities" begin
    # Test that is_ecryptfs works for something we're certain isn't encrypted
    if isdir("/proc")
        isecfs = (false, "/proc/")
        @test_logs (:info, "Checking to see if /proc/ is encrypted...") @test BinaryBuilderBase.is_ecryptfs("/proc"; verbose=true) == isecfs
        @test_logs (:info, "Checking to see if /proc/ is encrypted...") @test BinaryBuilderBase.is_ecryptfs("/proc/"; verbose=true) == isecfs
        @test_logs (:info, "Checking to see if /proc/not_a_file is encrypted...") @test BinaryBuilderBase.is_ecryptfs("/proc/not_a_file"; verbose=true) == isecfs
    else
        @test_logs (:info, "Checking to see if /proc/ is encrypted...") @test BinaryBuilderBase.is_ecryptfs("/proc"; verbose=true) == (false, "/proc")
        @test_logs (:info, "Checking to see if /proc/ is encrypted...") @test BinaryBuilderBase.is_ecryptfs("/proc/"; verbose=true) == (false, "/proc/")
        @test_logs (:info, "Checking to see if /proc/not_a_file is encrypted...") @test BinaryBuilderBase.is_ecryptfs("/proc/not_a_file"; verbose=true) == (false, "/proc/not_a_file")
    end

    if isa(preferred_runner(), BinaryBuilderBase.DockerRunner)
        @testset "Docker image importing" begin
            # First, delete the docker image, in case it already existed
            BinaryBuilderBase.delete_docker_image()

            # Next, import it and ensure that doesn't throw
            rootfs = first(BinaryBuilderBase.choose_shards(platform))
            mktempdir() do dir
                @test BinaryBuilderBase.import_docker_image(rootfs, dir; verbose=true) === nothing
            end

            # Test that deleting the docker image suceeds, now that we know
            # it exists
            @test BinaryBuilderBase.delete_docker_image()
        end
    end

    @testset "hello world" begin
        mktempdir() do dir
            ur = preferred_runner()(dir; platform=Platform("x86_64", "linux"; libc="musl"))
            iobuff = IOBuffer()
            @test run(ur, `/bin/bash -c "echo test"`, iobuff)
            seek(iobuff, 0)
            # Test that we get the output we expect (e.g. the second line is `test`)
            @test split(String(read(iobuff)), "\n")[2] == "test"
        end
    end

    if lowercase(get(ENV, "BINARYBUILDER_FULL_SHARD_TEST", "false")) == "true"
        @info("Beginning full shard test... (this can take a while)")
        platforms = supported_platforms()
    else
        platforms = (Platform("x86_64", "linux"; libc="musl"),)
    end

    # Checks that the wrappers provide the correct C++ string ABI
    @testset "Compilation - C++ string ABI" begin
        mktempdir() do dir
            # Host is x86_64-linux-musl-cxx11 and target is x86_64-linux-musl-cxx03
            ur = preferred_runner()(dir; platform=Platform("x86_64", "linux"; libc="musl", cxxstring_abi="cxx03"), preferred_gcc_version=v"5")
            iobuff = IOBuffer()
            test_script = raw"""
            set -e

            # Building for the target uses C++03 string ABI
            echo 'int main() {return 0;}' | SUPER_VERBOSE=1 ${CC} -x c - 2>&1 | grep -- "-D_GLIBCXX_USE_CXX11_ABI=0"

            # Building for the host uses C++11 string ABI
            echo 'int main() {return 0;}' | SUPER_VERBOSE=1 ${HOSTCC} -x c - 2>&1 | grep -v -- "-D_GLIBCXX_USE_CXX11_ABI=0"
            """
            @test run(ur, `/bin/bash -c "$(test_script)"`, iobuff; tee_stream=devnull)
        end
    end

    # This tests only that compilers for all platforms can build a simple C program
    @testset "Compilation - $(platform)" for platform in platforms
        mktempdir() do dir
            ur = preferred_runner()(dir; platform=platform)
            iobuff = IOBuffer()
            @test run(ur, `/bin/bash -c "echo 'int main() {return 0;}' | cc -x c -"`, iobuff; tee_stream=devnull)
            seekstart(iobuff)
            @test split(String(read(iobuff)), "\n")[2] == ""
        end
    end

    # This tests that compilers for all Intel Linux platforms can build a simple
    # C program that we can also run
    @testset "Compilation and running" begin
        mktempdir() do dir
            platforms = filter(p -> Sys.islinux(p) && proc_family(p) == "intel", supported_platforms())

            @testset "C - $(platform)" for platform in platforms
                ur = preferred_runner()(dir; platform=platform)
                iobuff = IOBuffer()
                test_c = """
                #include <stdio.h>
                int main() {
                    printf("Hello World!\\n");
                    return 0;
                }
                """
                cmd = `/bin/bash -c "echo '$(test_c)' > test.c && cc -o test test.c && ./test"`
                @test run(ur, cmd, iobuff; tee_stream=devnull)
                seekstart(iobuff)
                # Test that we get the output we expect
                @test endswith(readchomp(iobuff), "Hello World!")
            end

            # This tests that compilers for all Intel Linux platforms can build a simple
            # Fortran program that we can also run
            @testset "Fortran - $(platform)" for platform in filter(p -> Sys.islinux(p) && proc_family(p) == "intel", supported_platforms())
                ur = preferred_runner()(dir; platform=platform)
                iobuff = IOBuffer()
                test_f = """
                      program hello
                          print *, "Hello World!"
                      end program
                """
                cmd = `/bin/bash -c "echo '$(test_f)' > test.f && gfortran -o test test.f && ./test"`
                if arch(platform) == "i686" && libc(platform) == "musl"
                    @test_broken run(ur, cmd, iobuff; tee_stream=devnull)
                else
                    @test run(ur, cmd, iobuff; tee_stream=devnull)
                    seekstart(iobuff)
                    # Test that we get the output we expect
                    @test endswith(readchomp(iobuff), "Hello World!")
                end
            end
        end
    end

    @testset "Locking microarchitecture" begin
        mktempdir() do dir
            platform = Platform("x86_64", "linux"; libc="musl")
            cmd = `/bin/bash -c "echo 'int main() {return 0;}' | cc -x c -march=native -"`
            ur = preferred_runner()(dir; platform=platform, lock_microarchitecture=true)
            iobuff = IOBuffer()
            @test !run(ur, cmd, iobuff; tee_stream=devnull)
            seekstart(iobuff)
            @test split(String(read(iobuff)), "\n")[2] == "Cannot force an architecture"

            ur = preferred_runner()(dir; platform=platform, lock_microarchitecture=false)
            iobuff = IOBuffer()
            @test run(ur, cmd, iobuff; tee_stream=devnull)
            seekstart(iobuff)
            @test split(String(read(iobuff)), "\n")[2] == ""
        end
    end

    @testset "Unsafe flags" begin
        mktempdir() do dir
            platform = Platform("x86_64", "linux"; libc="musl")
            cmd = `/bin/bash -c "echo 'int main() {return 0;}' | cc -x c -Ofast -"`
            ur = preferred_runner()(dir; platform=platform, allow_unsafe_flags=false)
            iobuff = IOBuffer()
            @test !run(ur, cmd, iobuff; tee_stream=devnull)
            seekstart(iobuff)
            lines = split(String(read(iobuff)), "\n")
            @test lines[2] == "You used one or more of the unsafe flags: -Ofast, -ffast-math, -funsafe-math-optimizations"
            @test lines[3] == "Please repent."

            ur = preferred_runner()(dir; platform=platform, allow_unsafe_flags=true)
            iobuff = IOBuffer()
            @test run(ur, cmd, iobuff; tee_stream=devnull)
            seekstart(iobuff)
            @test split(String(read(iobuff)), "\n")[2] == ""
        end
    end
end


text:
function pascal(a::Integer, b::Integer, mid::Integer, top::Integer)
    yd = round((top - 4 * (a + b)) / 7)
    !isinteger(yd) && return 0, 0, 0
    y  = Int(yd)
    x  = mid - 2a - y
    return x, y, y - x
end

x, y, z = pascal(11, 4, 40, 151)
if !iszero(x)
    println("Solution: x = $x, y = $y, z = $z.")
else
    println("There is no solution.")
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder

name = "LibHSAKMT"
version = v"1.9.0"

# Collection of sources required to build LibHSAKMT
sources = [
    "https://github.com/RadeonOpenCompute/ROCT-Thunk-Interface/archive/roc-1.9.0.tar.gz" =>
    "00fe91f3dcb9f3246945a2421e3df3cf618bb392c1b22aa485cf74ebc816f7d0",

]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir
cd ROCT-Thunk-Interface-roc-1.9.0/
mkdir build
cd build
sed -i 's/^\(.*\)nodelete\(.*\)$/\1nodelete -L\/workspace\/destdir\/lib\/\2/g' ../CMakeLists.txt
cmake -DCMAKE_INSTALL_PREFIX=$prefix -DCMAKE_TOOLCHAIN_FILE=/opt/$target/$target.toolchain ..
make -j${nproc}
make install
rsync -a ../include/ $WORKSPACE/destdir/include/

"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    Linux(:x86_64, libc=:glibc),
    Linux(:aarch64, libc=:glibc),
    Linux(:powerpc64le, libc=:glibc),
    Linux(:aarch64, libc=:musl)
]

# The products that we will ensure are always built
products(prefix) = [
    LibraryProduct(prefix, "libhsakmt", Symbol(""))
]

# Dependencies that must be installed before this package can be built
dependencies = [
    "https://github.com/jpsamaroo/PCIUtilsBuilder/releases/download/v3.6.2/build_PCIUtils.v3.6.2.jl",
    "https://github.com/jpsamaroo/LibNUMABuilder/releases/download/v2.0.12/build_LibNUMA.v2.0.12.jl"
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: jpsamaroo/LibHSAKMTBuilder

max_issues_repo_name: jpsamaroo/LibHSAKMTBuilder

max_forks_repo_name: jpsamaroo/LibHSAKMTBuilder


text:
function longest_repeated_substring_T(text, nodes, root, minlen = nothing)

    if isnothing(minlen)
        minlen = 2
    end

    @assert minlen >= 1

    lenlongest = nothing
    ids = []

    strlongest = Dict()
    for nd in nodes

        if nd == root
            continue
        end

        repnode = length(nd.at) > 1
        if repnode
            @assert false
            mylen = nd.pflen + length(nd.s)
        else
            mylen = nd.pflen
        end

        if mylen < minlen
            continue
        end

        if isnothing(lenlongest) || mylen > lenlongest
            lenlongest = mylen
            ids = [nd.parent, nd]

            key = text[nd.sorg:nd.sorg+lenlongest-1]
            empty!(strlongest)
            strlongest[key] = union(Set(nd.parent.at), Set(nd.at))

        elseif mylen == lenlongest
            push!(ids, nd.parent)
            push!(ids, nd)

            key = text[nd.sorg:nd.sorg+lenlongest-1]

            if haskey(strlongest, key)
                union!(strlongest[key], Set(nd.parent.at))
                union!(strlongest[key], Set(nd.at))
            else
                strlongest[key] = union(Set(nd.parent.at), Set(nd.at))
            end

        end

    end

    if isnothing(lenlongest)
        return (nothing, strlongest)
    end

    return (lenlongest, strlongest)
end


text:
using Documenter, FoosAndBars

makedocs(
  sitename="FoosAndBars",
  doctest = false,
  strict = false,
  format = Documenter.HTML(prettyurls = get(ENV, "CI", nothing) == "true",),
  clean = false,
  pages = Any[
    "Home" => "index.md",
  ]
)

deploydocs(
           repo = "github.com/odunbar/FoosAndBars.git",
           target = "build",
          )

max_stars_repo_name: odunbar/FoosAndBars

max_issues_repo_name: odunbar/FoosAndBars

max_forks_repo_name: odunbar/FoosAndBars