text:
#* I'm not sure what the tests should be right now, currently there are no actual tests, just
#* print statements and comparisons from the model to experimental data... - Tyler 07/06/2020
using MotorPower
using Test

#MOTOR

#maxon 305013 EC-4pole 30 Ø30 mm, brushless, 200 Watt
threshold = 0.0001
Vm = 24
I0 = 0.728
Imax = 7.58
eta_max = 0.89
Kv = 700 #rpm/volt
Kv = Kv * pi / 30

Q, omega, Po, eta_M, mass_M, Ro = MotorPower.motorVI(Vm, Imax, Kv, I0; numMotors=1, Case=0, mtype = "hobby")
RPM = omega * 30 / pi

# Nominalvoltage=24V
# Noloadspeed=16700rpm
# Noloadcurrent=728mA
Nominalspeed = 16100 #rpm
Nominaltorquemaxcontinuoustorque = 0.0946 #Nm
# Nominalcurrent(max.continuouscurrent)=7.58A
# Stalltorque=3220mNm
# Stallcurrent=236A
Maxefficiency = 0.89
# Characteristics
Terminalresistance = 0.102#Ω
# Terminalinductance=0.016mH
# Torqueconstant=13.6mNm/A
# Speedconstant=700rpm/V
# Speed/torquegradient=5.21rpm/mNm
# Mechanicaltimeconstant=1.82ms
# Rotorinertia=33.3gcm²
# Thermaldata
# Thermalresistancehousing-ambient=7.4K/W
# Thermalresistancewinding-housing=0.209K/W
# Thermaltimeconstantwinding=2.11s
# Thermaltimeconstantmotor=1180s
# Ambienttemperature=-20...+100°C
# Max.windingtemperature=+155°C
# Mechanicaldata
# Bearingtype=ballbearings
# Max.speed=25000rpm
# Axialplay=0-0.14mm
# Maxaxialload(dynamic)=5.5N
# Maxforceforpressfits(static)=73N
# (static,shaftsupported)=1300N
# Max.radialload=25N,5mmfromflange
# Otherspecifications
# Numberofpolepairs=2
# Numberofphases=3
# Numberofautoclavecycles=0
# Product
Weight=300 #g
println("\nHOBBY")
V2, Im2, Pshaft2, eta_M2, mass_M2, Ro2 = MotorPower.motorVI(Nominaltorquemaxcontinuoustorque, Nominalspeed * pi / 30, Kv, I0; numMotors=1, Case=1, mtype = "hobby")
println("Omega1: $omega Q1: $Q V1: $Vm Im1: $Imax Pshaft1: $Po eta_M1: $eta_M mass_M: $mass_M Ro1: $Ro")
println("Omega2: $(Nominalspeed * pi / 30) Q2: $Nominaltorquemaxcontinuoustorque V2: $V2 Im2: $Im2 Pshaft2: $Pshaft2 eta_M2: $eta_M2 mass_M2: $mass_M2 Ro2: $Ro2 ")

println()
println("Motor Model Validation Against Maxon 305013 Brushless Motor:")
println("Motor Model Case: V and I in")
println("Motor Efficiency Differs by: $((eta_M - Maxefficiency) / Maxefficiency * 100)%")
println("Motor Torque Differs by: $((Q-Nominaltorquemaxcontinuoustorque)/Nominaltorquemaxcontinuoustorque * 100)%")
println("Motor RPM Differs by: $((RPM - Nominalspeed) / Nominalspeed * 100)%")
println("Motor Mass Differs by: $((mass_M * 1000 - Weight) / Weight * 100)%")
println("Motor Ro Differs by: $((Ro * 1000 - Terminalresistance) / Terminalresistance * 100)%")

println()
println("Motor Model Case: Q and Ω in")
println("Motor Efficiency Differs by: $((eta_M2 - Maxefficiency) / Maxefficiency * 100)%")
println("Motor Voltage Differs by: $((V2 - Vm) / Vm * 100)%")
println("Motor Current Differs by: $((Im2 - Imax) / Imax * 100)%")
println("Motor Mass Differs by: $((mass_M2 * 1000 - Weight) / Weight * 100)%")
println("Motor Ro Differs by: $((Ro2 * 1000 - Terminalresistance) / Terminalresistance * 100)%")

println("\nASTROFLIGHT")
Q, omega, Po, eta_M, mass_M, Ro = MotorPower.motorVI(Vm, Imax, Kv, I0; numMotors=1, Case=0, mtype = "astroflight")
V2, Im2, Pshaft2, eta_M2, mass_M2, Ro2 = MotorPower.motorVI(Nominaltorquemaxcontinuoustorque, Nominalspeed * pi / 30, Kv, I0, numMotors=1, Case=1, mtype ="astroflight")
println("Omega1: $omega Q1: $Q V1: $Vm Im1: $Imax Pshaft1: $Po eta_M1: $eta_M mass_M: $mass_M Ro1: $Ro")
println("Omega2: $(Nominalspeed * pi / 30) Q2: $Nominaltorquemaxcontinuoustorque V2: $V2 Im2: $Im2 Pshaft2: $Pshaft2 eta_M2: $eta_M2 mass_M2: $mass_M2 Ro2: $Ro2 ")

println()
println("Motor Model Validation Against Maxon 305013 Brushless Motor:")
println("Motor Model Case: V and I in")
println("Motor Efficiency Differs by: $((eta_M - Maxefficiency) / Maxefficiency * 100)%")
println("Motor Torque Differs by: $((Q - Nominaltorquemaxcontinuoustorque) / Nominaltorquemaxcontinuoustorque * 100)%")
println("Motor RPM Differs by: $((RPM - Nominalspeed) / Nominalspeed * 100)%")
println("Motor Mass Differs by: $((mass_M * 1000 - Weight) / Weight * 100)%")
println("Motor Ro Differs by: $((Ro * 1000 - Terminalresistance) / Terminalresistance * 100)%")

println()
println("Motor Model Case: Q and Ω in")
println("Motor Efficiency Differs by: $((eta_M2 - Maxefficiency) / Maxefficiency * 100)%")
println("Motor Voltage Differs by: $((V2 - Vm) / Vm * 100)%")
println("Motor Current Differs by: $((Im2 - Imax) / Imax * 100)%")
println("Motor Mass Differs by: $((mass_M2 * 1000 - Weight) / Weight * 100)%")
println("Motor Ro Differs by: $((Ro2 * 1000 - Terminalresistance) / Terminalresistance * 100)%")

#GEARBOX

#maxon 305013 EC-4pole 30 Ø30 mm, brushless, 200 Watt
threshold = 0.0001
Vm = 24
I0 = 0.728
Imax = 7.58
eta_max = 0.89
Kv = 700 #rpm/volt
Kv = Kv * pi / 30

Q, omega, Po, eta_M, mass_M, Ro = MotorPower.motorVI(Vm, Imax, Kv, I0; numMotors=1, Case=0, mtype = "hobby")
RPM = omega * 30 / pi

# Nominalvoltage=24V
# Noloadspeed=16700rpm
# Noloadcurrent=728mA
Nominalspeed = 16100 #rpm
Nominaltorquemaxcontinuoustorque = 0.0946 #Nm
# Nominalcurrent(max.continuouscurrent)=7.58A
# Stalltorque=3220mNm
# Stallcurrent=236A
Maxefficiency = 0.89
# Characteristics
Terminalresistance = 0.102#Ω
# Terminalinductance=0.016mH
# Torqueconstant=13.6mNm/A
# Speedconstant=700rpm/V
# Speed/torquegradient=5.21rpm/mNm
# Mechanicaltimeconstant=1.82ms
# Rotorinertia=33.3gcm²
# Thermaldata
# Thermalresistancehousing-ambient=7.4K/W
# Thermalresistancewinding-housing=0.209K/W
# Thermaltimeconstantwinding=2.11s
# Thermaltimeconstantmotor=1180s
# Ambienttemperature=-20...+100°C
# Max.windingtemperature=+155°C
# Mechanicaldata
# Bearingtype=ballbearings
# Max.speed=25000rpm
# Axialplay=0-0.14mm
# Maxaxialload(dynamic)=5.5N
# Maxforceforpressfits(static)=73N
# (static,shaftsupported)=1300N
# Max.radialload=25N,5mmfromflange
# Otherspecifications
# Numberofpolepairs=2
# Numberofphases=3
# Numberofautoclavecycles=0
# Product
Weight = 300 #g
println("\nHOBBY")
V2, Im2, Pshaft2, eta_M2, mass_M2, Ro2 = MotorPower.motorVI(Nominaltorquemaxcontinuoustorque, Nominalspeed * pi / 30, Kv, I0; numMotors=1, Case=1, mtype = "hobby")
println("Omega1: $omega Q1: $Q V1: $Vm Im1: $Imax Pshaft1: $Po eta_M1: $eta_M mass_M: $mass_M Ro1: $Ro")
println("Omega2: $(Nominalspeed * pi / 30) Q2: $Nominaltorquemaxcontinuoustorque V2: $V2 Im2: $Im2 Pshaft2: $Pshaft2 eta_M2: $eta_M2 mass_M2: $mass_M2 Ro2: $Ro2 ")

println()
println("Motor Model Validation Against Maxon 305013 Brushless Motor:")
println("Motor Model Case: V and I in")
println("Motor Efficiency Differs by: $((eta_M - Maxefficiency) / Maxefficiency * 100)%")
println("Motor Torque Differs by: $((Q - Nominaltorquemaxcontinuoustorque) / Nominaltorquemaxcontinuoustorque * 100)%")
println("Motor RPM Differs by: $((RPM - Nominalspeed) / Nominalspeed * 100)%")
println("Motor Mass Differs by: $((mass_M * 1000 - Weight) / Weight * 100)%")
println("Motor Ro Differs by: $((Ro * 1000 - Terminalresistance) / Terminalresistance * 100)%")

println()
println("Motor Model Case: Q and Ω in")
println("Motor Efficiency Differs by: $((eta_M2 - Maxefficiency) / Maxefficiency * 100)%")
println("Motor Voltage Differs by: $((V2 - Vm) / Vm * 100)%")
println("Motor Current Differs by: $((Im2 - Imax) / Imax * 100)%")
println("Motor Mass Differs by: $((mass_M2 * 1000 - Weight) / Weight * 100)%")
println("Motor Ro Differs by: $((Ro2 * 1000 - Terminalresistance) / Terminalresistance * 100)%")

println("\nASTROFLIGHT")
Q, omega, Po, eta_M, mass_M, Ro = MotorPower.motorVI(Vm, Imax, Kv, I0; numMotors=1, Case=0, mtype = "astroflight")
V2, Im2, Pshaft2, eta_M2, mass_M2, Ro2 = MotorPower.motorVI(Nominaltorquemaxcontinuoustorque, Nominalspeed * pi / 30, Kv, I0, numMotors=1, Case=1, mtype ="astroflight")
println("Omega1: $omega Q1: $Q V1: $Vm Im1: $Imax Pshaft1: $Po eta_M1: $eta_M mass_M: $mass_M Ro1: $Ro")
println("Omega2: $(Nominalspeed * pi / 30) Q2: $Nominaltorquemaxcontinuoustorque V2: $V2 Im2: $Im2 Pshaft2: $Pshaft2 eta_M2: $eta_M2 mass_M2: $mass_M2 Ro2: $Ro2 ")

println()
println("Motor Model Validation Against Maxon 305013 Brushless Motor:")
println("Motor Model Case: V and I in")
println("Motor Efficiency Differs by: $((eta_M - Maxefficiency) / Maxefficiency * 100)%")
println("Motor Torque Differs by: $((Q - Nominaltorquemaxcontinuoustorque) / Nominaltorquemaxcontinuoustorque * 100)%")
println("Motor RPM Differs by: $((RPM - Nominalspeed) / Nominalspeed * 100)%")
println("Motor Mass Differs by: $((mass_M * 1000 - Weight) / Weight * 100)%")
println("Motor Ro Differs by: $((Ro * 1000 - Terminalresistance) / Terminalresistance * 100)%")

println()
println("Motor Model Case: Q and Ω in")
println("Motor Efficiency Differs by: $((eta_M2 - Maxefficiency) / Maxefficiency * 100)%")
println("Motor Voltage Differs by: $((V2 - Vm) / Vm * 100)%")
println("Motor Current Differs by: $((Im2 - Imax) / Imax * 100)%")
println("Motor Mass Differs by: $((mass_M2 * 1000 - Weight) / Weight * 100)%")
println("Motor Ro Differs by: $((Ro2 * 1000 - Terminalresistance) / Terminalresistance * 100)%")


text:
struct ModalC0BasisFunction <: Field end

struct ModalC0Basis{D,T,V} <: AbstractVector{ModalC0BasisFunction}
  orders::NTuple{D,Int}
  terms::Vector{CartesianIndex{D}}
  a::Vector{Point{D,V}}
  b::Vector{Point{D,V}}
  function ModalC0Basis{D}(
    ::Type{T},
    orders::NTuple{D,Int},
    terms::Vector{CartesianIndex{D}},
    a::Vector{Point{D,V}},
    b::Vector{Point{D,V}}) where {D,T,V}
    new{D,T,V}(orders,terms,a,b)
  end
end

@inline Base.size(a::ModalC0Basis{D,T,V}) where {D,T,V} = (length(a.terms)*num_components(T),)
@inline Base.getindex(a::ModalC0Basis,i::Integer) = ModalC0BasisFunction()
@inline Base.IndexStyle(::ModalC0Basis) = IndexLinear()

function ModalC0Basis{D}(
  ::Type{T},
  orders::NTuple{D,Int},
  a::Vector{Point{D,V}},
  b::Vector{Point{D,V}};
  filter::Function=_q_filter,
  sort!::Function=_sort_by_nfaces!) where {D,T,V}

  terms = _define_terms_mc0(filter, sort!, orders)
  ModalC0Basis{D}(T,orders,terms,a,b)
end

function ModalC0Basis{D}(
  ::Type{T},
  orders::NTuple{D,Int},
  sa::Point{D,V},
  sb::Point{D,V};
  filter::Function=_q_filter,
  sort!::Function=_sort_by_nfaces!) where {D,T,V}

  terms = _define_terms_mc0(filter, sort!, orders)
  a = fill(sa,length(terms))
  b = fill(sb,length(terms))
  ModalC0Basis{D}(T,orders,terms,a,b)
end

function ModalC0Basis{D}(
  ::Type{T},
  orders::NTuple{D,Int};
  filter::Function=_q_filter,
  sort!::Function=_sort_by_nfaces!) where {D,T}

  sa = Point{D,eltype(T)}(tfill(zero(eltype(T)),Val{D}()))
  sb = Point{D,eltype(T)}(tfill(one(eltype(T)),Val{D}()))
  ModalC0Basis{D}(T,orders,sa,sb,filter=filter,sort! = sort!)
end

function ModalC0Basis{D}(
  ::Type{T},
  order::Int,
  a::Vector{Point{D,V}},
  b::Vector{Point{D,V}};
  filter::Function=_q_filter,
  sort!::Function=_sort_by_nfaces!) where {D,T,V}

  orders = tfill(order,Val{D}())
  ModalC0Basis{D}(T,orders,a,b,filter=filter,sort! = sort!)
end

function ModalC0Basis{D}(
  ::Type{T},
  order::Int;
  filter::Function=_q_filter,
  sort!::Function=_sort_by_nfaces!) where {D,T}

  orders = tfill(order,Val{D}())
  ModalC0Basis{D}(T,orders,filter=filter,sort! = sort!)
end

# API

"""
    get_order(b::ModalC0Basis)
"""
function get_order(b::ModalC0Basis)
  maximum(b.orders)
end

"""
    get_orders(b::ModalC0Basis)
"""
function get_orders(b::ModalC0Basis)
  b.orders
end

return_type(::ModalC0Basis{D,T,V}) where {D,T,V} = T

# Field implementation

function return_cache(f::ModalC0Basis{D,T,V},x::AbstractVector{<:Point}) where {D,T,V}
  @assert D == length(eltype(x)) "Incorrect number of point components"
  np = length(x)
  ndof = length(f.terms)*num_components(T)
  n = 1 + _maximum(f.orders)
  r = CachedArray(zeros(T,(np,ndof)))
  v = CachedArray(zeros(T,(ndof,)))
  c = CachedArray(zeros(eltype(T),(D,n)))
  (r, v, c)
end

function evaluate!(cache,f::ModalC0Basis{D,T,V},x::AbstractVector{<:Point}) where {D,T,V}
  r, v, c = cache
  np = length(x)
  ndof = length(f.terms)*num_components(T)
  n = 1 + _maximum(f.orders)
  setsize!(r,(np,ndof))
  setsize!(v,(ndof,))
  setsize!(c,(D,n))
  for i in 1:np
    @inbounds xi = x[i]
    _evaluate_nd_mc0!(v,xi,f.a,f.b,f.orders,f.terms,c)
    for j in 1:ndof
      @inbounds r[i,j] = v[j]
    end
  end
  r.array
end

function return_cache(
  fg::FieldGradientArray{1,ModalC0Basis{D,V,W}},
  x::AbstractVector{<:Point}) where {D,V,W}

  f = fg.fa
  @assert D == length(eltype(x)) "Incorrect number of point components"
  np = length(x)
  ndof = length(f.terms)*num_components(V)
  xi = testitem(x)
  T = gradient_type(V,xi)
  n = 1 + _maximum(f.orders)
  r = CachedArray(zeros(T,(np,ndof)))
  v = CachedArray(zeros(T,(ndof,)))
  c = CachedArray(zeros(eltype(T),(D,n)))
  g = CachedArray(zeros(eltype(T),(D,n)))
  (r, v, c, g)
end

function evaluate!(
  cache,
  fg::FieldGradientArray{1,ModalC0Basis{D,T,V}},
  x::AbstractVector{<:Point}) where {D,T,V}

  f = fg.fa
  r, v, c, g = cache
  np = length(x)
  ndof = length(f.terms) * num_components(T)
  n = 1 + _maximum(f.orders)
  setsize!(r,(np,ndof))
  setsize!(v,(ndof,))
  setsize!(c,(D,n))
  setsize!(g,(D,n))
  for i in 1:np
    @inbounds xi = x[i]
    _gradient_nd_mc0!(v,xi,f.a,f.b,f.orders,f.terms,c,g,T)
    for j in 1:ndof
      @inbounds r[i,j] = v[j]
    end
  end
  r.array
end

function return_cache(
  fg::FieldGradientArray{2,ModalC0Basis{D,V,W}},
  x::AbstractVector{<:Point}) where {D,V,W}

  f = fg.fa
  @assert D == length(eltype(x)) "Incorrect number of point components"
  np = length(x)
  ndof = length(f.terms)*num_components(V)
  xi = testitem(x)
  T = gradient_type(gradient_type(V,xi),xi)
  n = 1 + _maximum(f.orders)
  r = CachedArray(zeros(T,(np,ndof)))
  v = CachedArray(zeros(T,(ndof,)))
  c = CachedArray(zeros(eltype(T),(D,n)))
  g = CachedArray(zeros(eltype(T),(D,n)))
  h = CachedArray(zeros(eltype(T),(D,n)))
  (r, v, c, g, h)
end

function evaluate!(
  cache,
  fg::FieldGradientArray{2,ModalC0Basis{D,T,V}},
  x::AbstractVector{<:Point}) where {D,T,V}

  f = fg.fa
  r, v, c, g, h = cache
  np = length(x)
  ndof = length(f.terms) * num_components(T)
  n = 1 + _maximum(f.orders)
  setsize!(r,(np,ndof))
  setsize!(v,(ndof,))
  setsize!(c,(D,n))
  setsize!(g,(D,n))
  setsize!(h,(D,n))
  for i in 1:np
    @inbounds xi = x[i]
    _hessian_nd_mc0!(v,xi,f.a,f.b,f.orders,f.terms,c,g,h,T)
    for j in 1:ndof
      @inbounds r[i,j] = v[j]
    end
  end
  r.array
end

# Optimizing evaluation at a single point

function return_cache(f::AbstractVector{ModalC0BasisFunction},x::Point)
  xs = [x]
  cf = return_cache(f,xs)
  v = evaluate!(cf,f,xs)
  r = CachedArray(zeros(eltype(v),(size(v,2),)))
  r, cf, xs
end

function evaluate!(cache,f::AbstractVector{ModalC0BasisFunction},x::Point)
  r, cf, xs = cache
  xs[1] = x
  v = evaluate!(cf,f,xs)
  ndof = size(v,2)
  setsize!(r,(ndof,))
  a = r.array
  copyto!(a,v)
  a
end

function return_cache(
  f::FieldGradientArray{N,<:AbstractVector{ModalC0BasisFunction}}, x::Point) where {N}
  xs = [x]
  cf = return_cache(f,xs)
  v = evaluate!(cf,f,xs)
  r = CachedArray(zeros(eltype(v),(size(v,2),)))
  r, cf, xs
end

function evaluate!(
  cache, f::FieldGradientArray{N,<:AbstractVector{ModalC0BasisFunction}}, x::Point) where {N}
  r, cf, xs = cache
  xs[1] = x
  v = evaluate!(cf,f,xs)
  ndof = size(v,2)
  setsize!(r,(ndof,))
  a = r.array
  copyto!(a,v)
  a
end

# Helpers

_s_filter_mc0(e,o) = ( sum( [ i for i in e if i>1 ] ) <= o )

_sort_by_tensor_prod!(terms,orders) = terms

function _sort_by_nfaces!(terms::Vector{CartesianIndex{D}},orders) where D

  # Generate indices of n-faces and order s.t.
  # (1) dimension-increasing (2) lexicographic
  bin_rang_nfaces = tfill(0:1,Val{D}())
  bin_ids_nfaces = collect(Iterators.product(bin_rang_nfaces...))
  sum_bin_ids_nfaces = [sum(bin_ids_nfaces[i]) for i in eachindex(bin_ids_nfaces)]
  bin_ids_nfaces = permute!(bin_ids_nfaces,sortperm(sum_bin_ids_nfaces))

  # Generate LIs of basis funs s.t. order by n-faces
  lids_b = LinearIndices(Tuple([orders[i]+1 for i=1:D]))

  eet = eltype(eltype(bin_ids_nfaces))
  f(x) = Tuple( x[i] == one(eet) ? (0:0) : (1:2) for i in 1:length(x) )
  g(x) = Tuple( x[i] == one(eet) ? (3:orders[i]+1) : (0:0) for i in 1:length(x) )
  rang_nfaces = map(f,bin_ids_nfaces)
  rang_own_dofs = map(g,bin_ids_nfaces)

  P = Int64[]
  for i = 1:length(bin_ids_nfaces)
    cis_nfaces = CartesianIndices(rang_nfaces[i])
    cis_own_dofs = CartesianIndices(rang_own_dofs[i])
    for ci in cis_nfaces
      ci = ci .+ cis_own_dofs
      P = vcat(P,reshape(lids_b[ci],length(ci)))
    end
  end

  permute!(terms,P)
end

function _compute_filter_mask(terms,filter,orders)
  g = (0 .* orders) .+ 1
  to = CartesianIndex(g)
  maxorder = _maximum(orders)
  term_to_is_fterm = lazy_map(t->filter(Int[Tuple(t-to)...],maxorder),terms)
  findall(term_to_is_fterm)
end

function _define_terms_mc0(filter,sort!,orders)
  terms = _define_terms(_q_filter,orders)
  sort!(terms,orders)
  mask = _compute_filter_mask(terms,filter,orders)
  collect(lazy_map(Reindex(terms),mask))
end

function _evaluate_1d_mc0!(v::AbstractMatrix{T},x,a,b,order,d) where T
  @assert order > 0
  n = order + 1
  z = one(T)
  @inbounds v[d,1] = z - x[d]
  @inbounds v[d,2] = x[d]
  if n > 2
    ξ = ( 2*x[d] - ( a[d] + b[d] ) ) / ( b[d] - a[d] )
    for i in 3:n
      @inbounds v[d,i] = -sqrt(2*i-3)*v[d,1]*v[d,2]*jacobi(ξ,i-3,1,1)/(i-2)
    end
  end
end

function _gradient_1d_mc0!(v::AbstractMatrix{T},x,a,b,order,d) where T
  @assert order > 0
  n = order + 1
  z = one(T)
  @inbounds v[d,1] = -z
  @inbounds v[d,2] = z
  if n > 2
    ξ = ( 2*x[d] - ( a[d] + b[d] ) ) / ( b[d] - a[d] )
    v1 = z - x[d]
    v2 = x[d]
    for i in 3:n
      j, dj = jacobi_and_derivative(ξ,i-3,1,1)
      @inbounds v[d,i] = -sqrt(2*i-3)*(v[d,1]*v2*j+v1*v[d,2]*j+v1*v2*(2/(b[d]-a[d]))*dj)/(i-2)
    end
  end
end

function _hessian_1d_mc0!(v::AbstractMatrix{T},x,a,b,order,d) where T
  @assert order > 0
  n = order + 1
  y = zero(T)
  z = one(T)
  @inbounds v[d,1] = y
  @inbounds v[d,2] = y
  if n > 2
    ξ = ( 2*x[d] - ( a[d] + b[d] ) ) / ( b[d] - a[d] )
    v1 = z - x[d]
    v2 = x[d]
    dv1 = -z
    dv2 = z
    for i in 3:n
      j, dj = jacobi_and_derivative(ξ,i-3,1,1)
      _, d2j = jacobi_and_derivative(ξ,i-4,2,2)
      @inbounds v[d,i] = -sqrt(2*i-3)*(2*dv1*dv2*j+2*(dv1*v2+v1*dv2)*(2/(b[d]-a[d]))*dj+v1*v2*d2j*2*i*((b[d]-a[d])^2))/(i-2)
    end
  end
end

function _evaluate_nd_mc0!(
  v::AbstractVector{V},
  x,
  a::Vector{Point{D,T}},
  b::Vector{Point{D,T}},
  orders,
  terms::AbstractVector{CartesianIndex{D}},
  c::AbstractMatrix{T}) where {V,T,D}

  dim = D
  o = one(T)
  k = 1
  l = length(terms)

  for (i,ci) in enumerate(terms)

    for d in 1:dim
      _evaluate_1d_mc0!(c,x,a[i],b[i],orders[d],d)
    end

    s = o
    for d in 1:dim
      @inbounds s *= c[d,ci[d]]
    end

    k = _set_value_mc0!(v,s,k,l)

  end

end

@inline function _set_value_mc0!(v::AbstractVector{V},s::T,k,l) where {V,T}
  m = zero(Mutable(V))
  z = zero(T)
  js = eachindex(m)
  for j in js
    for i in js
      @inbounds m[i] = z
    end
    @inbounds m[j] = s
    i = k+l*(j-1)
    @inbounds v[i] = m
  end
  k+1
end

@inline function _set_value_mc0!(v::AbstractVector{<:Real},s,k,l)
  @inbounds v[k] = s
  k+1
end

function _gradient_nd_mc0!(
  v::AbstractVector{G},
  x,
  a::Vector{Point{D,T}},
  b::Vector{Point{D,T}},
  orders,
  terms::AbstractVector{CartesianIndex{D}},
  c::AbstractMatrix{T},
  g::AbstractMatrix{T},
  ::Type{V}) where {G,T,D,V}

  dim = D
  z = zero(Mutable(VectorValue{D,T}))
  o = one(T)
  k = 1
  l = length(terms)

  for (i,ci) in enumerate(terms)

    for d in 1:dim
      _evaluate_1d_mc0!(c,x,a[i],b[i],orders[d],d)
      _gradient_1d_mc0!(g,x,a[i],b[i],orders[d],d)
    end

    s = z
    for i in eachindex(s)
      @inbounds s[i] = o
    end
    for q in 1:dim
      for d in 1:dim
        if d != q
          @inbounds s[q] *= c[d,ci[d]]
        else
          @inbounds s[q] *= g[d,ci[d]]
        end
      end
    end

    k = _set_gradient_mc0!(v,s,k,l,V)

  end

end

@inline function _set_gradient_mc0!(
  v::AbstractVector{G},s,k,l,::Type{<:Real}) where G

  @inbounds v[k] = s
  k+1
end

@inline function _set_gradient_mc0!(
  v::AbstractVector{G},s,k,l,::Type{V}) where {V,G}

  T = eltype(s)
  m = zero(Mutable(G))
  w = zero(V)
  z = zero(T)
  for (ij,j) in enumerate(CartesianIndices(w))
    for i in CartesianIndices(m)
      @inbounds m[i] = z
    end
    for i in CartesianIndices(s)
      @inbounds m[i,j] = s[i]
    end
    i = k+l*(ij-1)
    @inbounds v[i] = m
  end
  k+1
end

function _hessian_nd_mc0!(
  v::AbstractVector{G},
  x,
  a::Vector{Point{D,T}},
  b::Vector{Point{D,T}},
  orders,
  terms::AbstractVector{CartesianIndex{D}},
  c::AbstractMatrix{T},
  g::AbstractMatrix{T},
  h::AbstractMatrix{T},
  ::Type{V}) where {G,T,D,V}

  dim = D
  z = zero(Mutable(TensorValue{D,D,T}))
  o = one(T)
  k = 1
  l = length(terms)

  for (i,ci) in enumerate(terms)

    for d in 1:dim
      _evaluate_1d_mc0!(c,x,a[i],b[i],orders[d],d)
      _gradient_1d_mc0!(g,x,a[i],b[i],orders[d],d)
      _hessian_1d_mc0!(h,x,a[i],b[i],orders[d],d)
    end

    s = z
    for i in eachindex(s)
      @inbounds s[i] = o
    end
    for r in 1:dim
      for q in 1:dim
        for d in 1:dim
          if d != q && d != r
            @inbounds s[r,q] *= c[d,ci[d]]
          elseif d == q && d ==r
            @inbounds s[r,q] *= h[d,ci[d]]
          else
            @inbounds s[r,q] *= g[d,ci[d]]
          end
        end
      end
    end

    k = _set_gradient_mc0!(v,s,k,l,V)

  end

end

ext: jl


text:
#=
Stock dynamics the a lake model.

=#

using Plots
pyplot()
include("lake_model.jl")

lm = LakeModel()
N_0 = 150      # Population
e_0 = 0.92     # Initial employment rate
u_0 = 1 - e_0  # Initial unemployment rate
T = 50         # Simulation length

E_0 = e_0 * N_0
U_0 = u_0 * N_0
X_0 = [E_0; U_0]

X_path = simulate_stock_path(lm, X_0, T)

titles = ["Employment" "Unemployment" "Labor Force"]
dates = collect(1:T)

x1 = squeeze(X_path[1, :], 1)
x2 = squeeze(X_path[2, :], 1)
x3 = squeeze(sum(X_path, 1), 1)

plot(dates, Vector[x1, x2, x3], title=titles, layout=(3, 1), legend=:none)

max_stars_repo_name: mwhchen/quantecon

max_issues_repo_name: mwhchen/quantecon

max_forks_repo_name: mwhchen/quantecon


text:
using OrdinaryDiffEq, DifferentialEquations

tspan = (0.0f0, 50.0f0)

"""
Simple SIR disease model from:
    https://en.wikipedia.org/wiki/Compartmental_models_in_epidemiology

Note: this does not include birth and death rates
"""
function sir_system(du, u, p, t;
                     U_func_time=U_func_time_trivial,
                     U_func_space=U_func_space_trivial)
    Ft = U_func_time(t)
    Fs = U_func_space(u)
    beta, gamma = p
    S, I, R = u

    N = sum([S,I,R])
    tmp = beta*I*S/N
    du[1] = -tmp
    du[2] = tmp - gamma*I
    du[3] = gamma*I

    # This forcing function must keep the overall population constant
    du[1] += -Ft[1]
    du[3] += Ft[1]
    # du .+= Ft .+ Fs
end

# Generate data, with a time component
p = [0.5, 0.2]
u0 = [999.0, 1.0, 0.0]
ts = range(tspan[1], tspan[2], length=5001)

function solve_sir_system(;U_func_time=U_func_time_trivial,
                              U_func_space=U_func_space_trivial)
    prob = ODEProblem((du, u, p, t)->sir_system(du, u, p, t,
                                U_func_time=U_func_time,
                                U_func_space=U_func_space),
                    u0, tspan, p)
    sol = solve(prob, Tsit5(), saveat=ts);#, adaptive=true,
        #dt=1e-5);#, abstol=1e-12, reltol=1e-5);
    # sol = solve(prob, AB4(), saveat=ts, dt=1e-5);
    return sol
end

#####
##### Turing function for Bayesian parameter estimation
#####
# Unforced version for use with Turing.jl
function sir_system(u, p, t)
    u = convert.(eltype(p),u)
    beta, gamma = p
    S, I, R = u

    N = sum([S,I,R])
    tmp = beta*I*S/N
    du = [-tmp; tmp - gamma*I; gamma*I]
    return du
end

#####
##### True model in SINDy syntax
#####
a, b = p
a = a / sum(u0)
#      S  I  R  c SS SI SR II IR RR
 A = [[0  0  0  0  0 -a  0  0  0  0];
      [0 -b  0  0  0  a  0  0  0  0];
      [0  b  0  0  0  0  0  0  0  0]]
n = size(A, 1)
sindy_library = Dict(
    "cross_terms"=>2,
    "constant"=>nothing
);
core_dyn_true = sindyModel(ts, A,
                            convert_string2function(sindy_library),
                            ["S", "I", "R"])

# Actually export
export ts, solve_sir_system, sir_system, core_dyn_true
    # ,lorenz_grad_residual


text:
"""Julia program to implement Linear Search algorithm.
Given an array, Linear Search searches for an element by traversing the entire array once
and comparing each element with the element to be found out.
"""

function linear_search(arr, n, ele)
    # Traverse the entire array once
    for i in 1:n
        if(arr[i] == ele)
            return true
        end
    end
    return false
end

print("How many numbers are present in the array? ")
n = readline()
n = parse(Int, n)
if (n <= 0)
    println("Array is Empty!!!")
    exit()
end
arr = Int[]
print("Enter the numbers: ")
arr = [parse(Int, num) for num in split(readline())]
print("Which number do you want to search in the array? ")
ele = readline()
ele = parse(Int, ele)
res = linear_search(arr, n, ele)
if res
    print("The number $ele is present in the array.")
else
    print("The number $ele is not present in the array")
end

"""
Time Complexity - O(n), where 'n' is the size of the array
Space Complexity - O(n)

SAMPLE INPUT AND OUTPUT

SAMPLE I

How many numbers are present in the array? 5
Enter the numbers: 1 2 3 4 5
Which number do you want to search in the array? 6
The number 6 is not present in the array

SAMPLE II

How many numbers are present in the array? 3
Enter the numbers: 3 1 2
Which number do you want to search in the array? 2
The number 2 is present in the array.

"""


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

"""
    Ray(p, v)

A ray originating at point `p`, pointed in direction `v`.
It can be called as `r(t)` with `t > 0` to cast it at
`p + t * v`.
"""
struct Ray{Dim,T} <: Primitive{Dim,T}
  p::Point{Dim,T}
  v::Vec{Dim,T}
end

Ray(p::Tuple, v::Tuple) = Ray(Point(p), Vec(v))

paramdim(::Type{<:Ray}) = 1

isconvex(::Type{<:Ray}) = true

function (r::Ray)(t)
  if t < 0
    throw(DomainError(t, "r(t) is not defined for t < 0."))
  end
  r.p + t * r.v
end


text:
"""
Graph Laplacian of the line graph. With keyword `boundary = false`,
a graph Laplacian of the circle.
"""
linelaplacian(n; kwargs...) = linelaplacian(Float64,n; kwargs...)
function linelaplacian(T,n; boundary=true)
    A = SymTridiagonal(2ones(T,n), -ones(T,n-1))
    if boundary
        A[1,1] = A[end,end] = true
    end
    A
end

"""
Graph Laplacian of a `m×n` lattice.
"""
function gridlaplacian(T, m, n)
    S = sparse(T(0.0)I, n*m, n*m)
    linear = LinearIndices((1:m, 1:n))
    for i in 1:m
        for j in 1:n
            for (i2, j2) in ((i + 1, j), (i, j + 1))
                if i2 <= m && j2 <= n
                    S[linear[i, j], linear[i2, j2]] -= 1
                    S[linear[i2, j2], linear[i, j]] -= 1

                    S[linear[i, j], linear[i, j]] += 1
                    S[linear[i2, j2], linear[i2, j2]] += 1
                end
            end
        end
    end
    S
end

function gridderiv(T, m, n)
    S1 = sparse(T(0.0)*I, n*m, n*m)
    S2 = sparse(T(0.0)*I, n*m, n*m)
    S3 = sparse(T(0.0)*I, n*m, n*m)
    S4 = sparse(T(0.0)*I, n*m, n*m)

    linear = LinearIndices((1:m, 1:n))
    for i in 1:m
        for j in 1:n
            if i > 1
                S1[linear[i, j], linear[i, j]] -= 1
                S1[linear[i, j], linear[i - 1, j]] += 1
            end
            if j > 1
                S2[linear[i, j], linear[i, j]] -= 1
                S2[linear[i, j], linear[i, j - 1]] += +1
            end
            if i < m
                S3[linear[i, j], linear[i, j]] -= 1
                S3[linear[i, j], linear[i + 1, j]] += 1
            end
            if j < n
                S4[linear[i, j], linear[i, j]] -= 1
                S4[linear[i, j], linear[i, j + 1]] += +1
            end
        end
    end
    S1, S2, S3, S4
end

function boundary(A)
    m, n = size(A)
    B = zero(A)
    for i in 1:m
        for j in 1:n
            v = A[i,j]
            for (i2, j2) in ((i + 1, j), (i + 1, j +1 ), (i, j + 1))
                if 1 <= i2 <= m && 1<= j2 <= n
                    if v != A[i2, j2]
                        B[i,j] = 1.0
                        break
                    end
                end
            end
        end
    end
    B
end

function downop(T, m, n)
    S = sparse(T(0.0)I, (n÷2)*(m÷2), n*m)
    linearj = LinearIndices((1:m, 1:n))
    lineari = LinearIndices((1:m÷2, 1:n÷2))

    for i in 2:2:m
        for j in 2:2:n
            S[lineari[i÷2,j÷2], linearj[i, j]] = 1/4
            S[lineari[i÷2,j÷2], linearj[i-1, j]] = 1/4
            S[lineari[i÷2,j÷2], linearj[i, j-1]] = 1/4
            S[lineari[i÷2,j÷2], linearj[i-1, j-1]] = 1/4
        end
    end
    S
end
#=
m = n = 8
L = downop(F0, m, n); Λ = gridlaplacian(F0, m, n); Λ2 = gridlaplacian(F0, m÷2, n÷2)
heatmap([Matrix(8L*Λ*L') Matrix(Λ2) Matrix(8L*Λ*L' - Λ2)])
L = downop(F0, m, n); Λ = gridlaplacian(F0, m, n) + I/2; Λ2 = gridlaplacian(F0, m÷2, n÷2) + I
invm(x) = inv(Matrix(x))
F = Float64
L = downop(F, m, n); Λ = gridlaplacian(F, m, n)/2 + I/4; Λ2 = gridlaplacian(F, m÷2, n÷2) + I
heatmap([Matrix(L*invm(Λ)*L') Matrix(invm(Λ2)) Matrix(L*invm(Λ)*L' - invm(Λ2))])
mean(diag(L*invm(Λ)*L') - diag(invm(Λ2)))
.32007- .32008
=#


text:
# const TRAVIS_SKIP = VERSION.minor != 4 && !isnothing(get(ENV, "TRAVIS_BRANCH", nothing))

function test_fixed_size(M, K, N)
    A = @StrideArray rand(M,K);
    B = @StrideArray rand(K,N);
    At = (@StrideArray rand(K,M))';
    Bt = (@StrideArray rand(N,K))';
    Aa = Array(A); Ba = Array(B);
    Aat = Array(At); Bat = Array(Bt);
    time = @elapsed(@test Aa * Ba ≈ A * B)
    @show M, K, N, time
    time = @elapsed(@test Aa * Bat ≈ A * Bt)
    @show M, K, N, time
    time = @elapsed(@test Aat * Ba ≈ At * B)
    @show M, K, N, time
    time = @elapsed(@test Aat * Bat ≈ At * Bt)
    @show M, K, N, time
    nothing
end

function test_fixed_size(M)
    K = N = M
    A = @StrideArray rand(M,K);
    B = @StrideArray rand(K,N);
    Aa = Array(A); Ba = Array(B);
    time = @elapsed(@test Aa * Ba ≈ A * B)
    @show M, K, N, time
    time = @elapsed(@test Aa * Ba' ≈ A * B')
    @show M, K, N, time
    time = @elapsed(@test Aa' * Ba ≈ A' * B)
    @show M, K, N, time
    time = @elapsed(@test Aa' * Ba' ≈ A' * B')
    @show M, K, N, time
    nothing
end

# gmul(A, B) = LinearAlgebra.generic_matmatmul('N','N', A, B)

@testset "MatrixMultiply" begin
    @testset "matmul" begin
        @time for T in (Float32, Float64, Int32, Int64)
            @show T, @__LINE__
            logMmax, logKmax, logNmax = log.(PaddedMatrices.matmul_params(T))
            @time for logN ∈ range(0, min(logNmax + 0.7, 9.210340371976184), length = 7)
                N = round(Int, exp(logN))
                for logM ∈ range(0, logMmax + 0.7, length = 7)
                    M = round(Int, exp(logM))
                    C1 = Matrix{T}(undef, M, N);
                    C2 = similar(C1); C3 = similar(C1); C4 = similar(C1); C5 = similar(C1);
                    for logK ∈ range(0, logKmax + 0.7, length = 7)
                        K = round(Int, exp(logK))

                        gopc = 2e-9M*K*N
                        A = rand(T, M, K); At = copy(A');
                        B = rand(T, K, N); Bt = copy(B');

                        blastime = @elapsed mul!(C1, A, B)
                        # blastime = @elapsed (C1 = gmul(A, B))
                        bops = gopc / blastime
                        @show (M,K,N), blastime, bops, T
                        s_time_nn = @elapsed matmul_serial!(C2, A, B); s_ops_nn = gopc / s_time_nn
                        @test C1 ≈ C2
                        s_time_nt = @elapsed matmul_serial!(C3, A, Bt'); s_ops_nt = gopc / s_time_nt
                        @test C1 ≈ C3
                        @show (M,K,N), s_time_nn, s_ops_nn, s_time_nt, s_ops_nt
                        if T <: AbstractFloat
                            fill!(C2, NaN); fill!(C3, NaN)
                        else
                            fill!(C2, -99999); fill!(C3, -99999)
                        end
                        t_time_nn = @elapsed matmul!(C2, A, B); t_ops_nn = gopc / t_time_nn
                        @test C1 ≈ C2
                        t_time_nt = @elapsed matmul!(C3, A, Bt'); t_ops_nt = gopc / t_time_nt
                        @test C1 ≈ C3
                        @show (M,K,N), t_time_nn, t_ops_nn, t_time_nt, t_ops_nt

                        s_time_tn = @elapsed matmul_serial!(C4, At', B); s_ops_tn = gopc / s_time_tn
                        @test C1 ≈ C4
                        s_time_tt = @elapsed matmul_serial!(C5, At', Bt'); s_ops_tt = gopc / s_time_tt
                        @test C1 ≈ C5
                        @show (M,K,N), s_time_tn, s_ops_tn, s_time_tt, s_ops_tt
                        if T <: AbstractFloat
                            fill!(C4, NaN); fill!(C5, NaN)
                        else
                            fill!(C4, -99999); fill!(C5, -99999)
                        end
                        t_time_tn = @elapsed matmul!(C4, At', B); t_ops_tn = gopc / t_time_tn
                        @test C1 ≈ C4
                        t_time_tt = @elapsed matmul!(C5, At', Bt'); t_ops_tt = gopc / t_time_tt
                        @test C1 ≈ C5
                        @show (M,K,N), t_time_tn, t_ops_tn, t_time_tt, t_ops_tt
                    end
                end
            end
        end
    end
    @testset "FixedSize" begin
        r = 2:7
        for M ∈ r, K ∈ r, N ∈ r
            test_fixed_size(M, K, N)
        end
        r = 8:33
        for M ∈ r
            test_fixed_size(M)
        end
        M = K = N = 80
        A = @StrideArray rand(M,K);
        B = @StrideArray rand(K,N);
        C = StrideArray{Float64}(undef, (StaticInt(M),StaticInt(N)));
        M, K, N = 23, 37, 19
        @views begin
            Av = A[1:M, 1:K];
            Bv = B[1:K, 1:N];
            Cv = C[1:M, 1:N];
            Avsl = A[StaticInt(1):M, StaticInt(1):K];
            Bvsl = B[StaticInt(1):K, StaticInt(1):N];
            Cvsl = C[StaticInt(1):M, StaticInt(1):N];
            Avsr = A[StaticInt(1):StaticInt(M), StaticInt(1):StaticInt(K)];
            Bvsr = B[StaticInt(1):StaticInt(K), StaticInt(1):StaticInt(N)];
            Cvsr = C[StaticInt(1):StaticInt(M), StaticInt(1):StaticInt(N)];
        end
        Creference = Array(Av) * Array(Bvsl);
        time = @elapsed mul!(Cv, Av, Bv)
        @test Creference ≈ Cv
        @show M, K, N, time
        time = @elapsed mul!(Cvsl, Avsl, Bvsl)
        @test Creference ≈ Cv
        @show M, K, N, time
        time = @elapsed mul!(Cvsr, Avsr, Bvsr)
        @test Creference ≈ Cv
        @show M, K, N, time

        Av = A[1:M, 1:K];
        Bv = B[1:K, 1:N];
        Cv = C[1:M, 1:N];
        Avsl = A[StaticInt(1):M, StaticInt(1):K];
        Bvsl = B[StaticInt(1):K, StaticInt(1):N];
        Cvsl = C[StaticInt(1):M, StaticInt(1):N];
        Avsr = A[StaticInt(1):StaticInt(M), StaticInt(1):StaticInt(K)];
        Bvsr = B[StaticInt(1):StaticInt(K), StaticInt(1):StaticInt(N)];
        Cvsr = C[StaticInt(1):StaticInt(M), StaticInt(1):StaticInt(N)];
        Creference = Array(Av) * Array(Bvsl);
        time = @elapsed mul!(Cv, Av, Bv)
        @test Creference ≈ Cv
        @show M, K, N, time
        time = @elapsed mul!(Cvsl, Avsl, Bvsl)
        @test Creference ≈ Cv
        @show M, K, N, time
        time = @elapsed mul!(Cvsr, Avsr, Bvsr)
        @test Creference ≈ Cv
        @show M, K, N, time

    end
end


text:
@testset "NoneSelector" begin
    st = Date(2019, 1, 1)
    ed = Date(2019, 2, 1)
    date_range = st:Day(1):ed

    selector = NoneSelector()
    @test selector isa DateSelector

    validation, holdout = partition(date_range, selector)
    @test validation == date_range
    @test isempty(holdout)
end


text:
"""
Returns bonding rules involving R-group-tagged atoms
"""
function tagged_bonding_rules()::Array{BondingRule}
    newrules = []
    for rule ∈ rc[:bonding_rules]
        if rule.species_i != :*
            push!(newrules, BondingRule(Symbol("$(rule.species_i)!"), rule.species_j, rule.max_dist))
            push!(newrules, BondingRule(Symbol("$(rule.species_j)!"), rule.species_i, rule.max_dist))
            push!(newrules, BondingRule(Symbol("$(rule.species_i)!"), Symbol("$(rule.species_j)!"), rule.max_dist))
        end
    end
    return newrules
end

"""
Returns R group indices (whichever atoms have species symbols appended by '!')
"""
function r_group_indices(xtal::Crystal)::Array{Int}
    R = []
    for (idx, label) in enumerate(xtal.atoms.species) # loop over crystal atoms to find tags
        # if String representation of label Symbol ends in !, atom is in R
        tokens = split("$label", rc[:r_tag])
        if length(tokens) == 2 && tokens[2] == "" # other ! in symbol not tolerated.
            push!(R, idx)
        end
    end
    return R
end

"""
Un-tags R group atoms (removes '!' suffix)
"""
function untag_r_group!(xtal::Crystal)
    r = r_group_indices(xtal) # get indices of R group
    for i ∈ r
        xtal.atoms.species[i] = Symbol(split("$(xtal.atoms.species[i])", rc[:r_tag])[1])
    end
end

"""
Returns a copy of a crystal w/ R group atoms deleted
"""
function subtract_r_group(xtal::Crystal)::Crystal
    not_r = [i for i ∈ 1:length(xtal.atoms.species) if !(i ∈ r_group_indices(xtal))]
    coords = @view xtal.atoms.coords[not_r]
    species = @view xtal.atoms.species[not_r]
    return Crystal("no_r_$(xtal.name)", xtal.box, Atoms(species, coords), xtal.charges)
end

## moiety import function (exposed)
@doc raw"""

    q = moiety(xyz_filename)

Generates a moiety (`Crystal`) from an .xyz file found in `rc[:paths][:moieties]`.

Use `set_path_to_data` or set `rc[:paths][:moieties]` to change the path from which the XYZ file is read.

Atoms appended with '!' are tagged for replacement via `substructure_replace`.

Bonds are inferred automatically via `infer_bonds!`.

# Arguments
- `xyz_filename::Union{String,Nothing}` the moiety input file name, an `.xyz` file; if set to `nothing` the moiety is the null set.
- `bonding_rules::Union{Vector{BondingRule},Nothing}` (optional) a list of rules to use for inferring the bonding network of the atoms loaded from the XYZ file. If set to `nothing`, the default rules are used.
"""
function moiety(name::Union{String,Nothing}; bonding_rules::Union{Vector{BondingRule},Nothing}=nothing)::Crystal
    # make box (arbitrary unit cube)
    box = unit_cube()
    # handle deletion option (replace-with-nothing)
    if !isnothing(name)
        fx = Frac(read_xyz("$(rc[:paths][:moieties])/$name"), box)
    else
        name = "nothing"
        fx = Atoms{Frac}(0)
    end
    # generate Crystal from moiety XYZ coords
    charges = Charges{Frac}(0)
    moiety = Crystal(name, box, fx, charges)
    # ID R group
    R_group_indices = r_group_indices(moiety)
    # handle custom vs. default bonding rules
    if isnothing(bonding_rules)
        infer_bonds!(moiety, false)
    else
        infer_bonds!(moiety, false; bonding_rules=bonding_rules)
    end
    # sort by node degree
    order = sortperm(degree(moiety.bonds), rev=true)
    # ordered atoms
    if length(R_group_indices) > 0
        order_wo_R = order[
            [i for i ∈ 1:length(order) if !(order[i] ∈ R_group_indices)]]
    else
        order_wo_R = order
    end
    # append R-group to the end
    order = vcat(order_wo_R, R_group_indices)
    # rebuild Atoms
    atoms = Atoms(moiety.atoms.species[order], moiety.atoms.coords[order])
    # nodes are sorted by bond order, and R group is moved to end & tagged w/ !
    moiety = Crystal(name, box, atoms, charges)
    infer_bonds!(moiety, false)
    return moiety
end


text:
using Ratios, Test
using FixedPointNumbers

@testset "SimpleRatio" begin
    r = SimpleRatio(1,2)
    @test convert(Float64, r) == 0.5
    @test convert(Float32, r) == 0.5f0
    @test convert(BigFloat, r) == BigFloat(1)/2

    r2 = SimpleRatio(2,3)
    @test r*r2 == SimpleRatio(2,6) == SimpleRatio(1,3)
    @test r2*3 == 3*r2 == 2
    @test r*false == false*r == 0
    @test r/r2 == SimpleRatio(3,4)
    @test r/2 == SimpleRatio(1,4)
    @test 2/r == 4
    @test 4 == 2/r
    @test r+1 == 1+r == SimpleRatio(3,2)
    @test r-1 == SimpleRatio(-1,2)
    @test 1-r == r
    @test r+r2 == SimpleRatio(7,6)
    @test r-r2 == SimpleRatio(-1,6)
    @test r^2 == SimpleRatio(1,4)
    @test -r == SimpleRatio(-1,2)
    @test 0.2*r ≈ 0.1
    @test r == 0.5
    @test 0.5 == r

    r3 = SimpleRatio(7,3)
    @test r2+r3 === SimpleRatio(9,3)
    @test r2-r3 === SimpleRatio(-5,3)

    r8 = SimpleRatio{Int8}(1, 20)
    @test r8 + r8 == SimpleRatio(1, 10)

    @test_throws OverflowError -SimpleRatio(0x02,0x03)

    @test r + SimpleRatio(0x02,0x03) == SimpleRatio(7,6)

    @test SimpleRatio(11, 10) == 11//10
    @test 1//3 + SimpleRatio(1, 5) == 8//15

    @test isfinite(SimpleRatio(0,0)) == false
    @test isfinite(SimpleRatio(1,0)) == false
    @test isfinite(SimpleRatio(2,1)) == true

    @test SimpleRatio(5,3) * 0.035N0f8 == SimpleRatio{Int}(rationalize((5*0.035N0f8)/3))
    @test SimpleRatio(5,3) * 0.035N4f12 == SimpleRatio{Int}(rationalize((5*0.035N4f12)/3))
    @test SimpleRatio(5,3) * -0.03Q0f7 == SimpleRatio{Int}(rationalize((5.0*(-0.03Q0f7))/3))
    r = @inferred(SimpleRatio(0.75Q0f7))
    @test r == 3//4 && r isa SimpleRatio{Int16}

    # common_denominator
    @test common_denominator(SimpleRatio(2,7), SimpleRatio(3,11), SimpleRatio(-1,5)) ===
        (SimpleRatio(2*11*5,385), SimpleRatio(3*7*5,385), SimpleRatio(-1*7*11,385))
    @test common_denominator(SimpleRatio(5,12), SimpleRatio(4,15), SimpleRatio(-1,9)) ===
        (SimpleRatio(75,180), SimpleRatio(48,180), SimpleRatio(-20,180))
    @test_throws OverflowError common_denominator(SimpleRatio{Int8}(1, 20), SimpleRatio{Int8}(2, 21))
    @test common_denominator(SimpleRatio{Int8}(1, 20), SimpleRatio{Int8}(3, 20)) ===
        (SimpleRatio{Int8}(1, 20), SimpleRatio{Int8}(3, 20))
end


text:
struct Simple
    baz::Int
end

struct Composite
    composite::Simple
end


text:
name = "LLVM"
llvm_full_version = v"9.0.1+4"
libllvm_version = v"9.0.1+5"

# Include common LLVM stuff
include("../common.jl")
build_tarballs(ARGS, configure_extraction(ARGS, llvm_full_version, name, libllvm_version)...; skip_audit=true)

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
# TODO: provide outer constructor checking dimensions?
# TODO: think about associating it to a coordinate system
"""
    PointForcesMoments

Forces and moments applied on a point. This objects can be added, translated or
rotated to be expressed in other coordinate system.

# Constructor

    PointForcesMoments(point::Array{Number, 1}, forces::Array{Number, 1},
                       moments::Array{Number, 1})
"""
struct PointForcesMoments
    point::Array{Number, 1}
    forces::Array{Number, 1}
    moments::Array{Number, 1}
end

get_point(pfm::PointForcesMoments) = pfm.point
get_forces(pfm::PointForcesMoments) = pfm.forces
get_moments(pfm::PointForcesMoments) = pfm.moments

"""
    translate_forces_moments(pfm::PointForcesMoments, point::Array{T, 1}) where T<:Number

Calculate equivalent `PointForcesMoments` in point Q given the forces and
moments at point P. Forces are preserved and moments are calculated taking into
account the arm between the two points.

``M_{Q} = M_{P} + r \times F``
"""
function translate_forces_moments(pfm::PointForcesMoments, point::Array{T, 1}) where T<:Number
    r = pfm.point - point
    if isapprox(r, [0., 0., 0.])
        return pfm
    else
        f = pfm.forces
        m = pfm.moments + cross(r, f)
        return PointForcesMoments(point, f, m)
    end
end

"""
    +(pfm1::PointForcesMoments, pfm2::PointForcesMoments, point::Array{T, 1})
      where T<:Number

Add two PointForcesMoments objects expressing the result in point P.
"""
function +(pfm1::PointForcesMoments, pfm2::PointForcesMoments, point::Array{T, 1}) where T<:Number

    pfm1 = translate_forces_moments(pfm1, point)
    pfm2 = translate_forces_moments(pfm2, point)
    PointForcesMoments(point, pfm1.forces+pfm2.forces, pfm1.moments+pfm2.moments)
end

"""
Add two PointForcesMoments objects. The result will be expressed in the first
point.
"""
function +(pfm1::PointForcesMoments, pfm2::PointForcesMoments)
    +(pfm1, pfm2, pfm1.point)
end

"""
Multiply the forces and moments by a number.
"""
function *(pfm::PointForcesMoments, n::Number)
    PointForcesMoments(pfm.point, n*pfm.forces, n*pfm.moments)
end

*(n::Number, pfm::PointForcesMoments) = *(pfm, n)

"""
    -(pfm1::PointForcesMoments, pfm2::PointForcesMoments, point::Array{T, 1})
      where T<:Number

Substract two PointForcesMoments objects expressing the result in point P.
"""
function -(pfm1::PointForcesMoments, pfm2::PointForcesMoments, point::Array{T, 1}) where T<:Number
    +(pfm1, (-1) * pfm2, point)
end

"""
Substract two PointForcesMoments objects. The result will be expressed in the first
point.
"""
function -(pfm1::PointForcesMoments, pfm2::PointForcesMoments)
    -(pfm1, pfm2, pfm1.point)
end

"""
Check if two PointForcesMoments are equivalent when expressed at the same point.
"""
function isapprox(pfm1::PointForcesMoments, pfm2::PointForcesMoments; kwargs...)
    trans_pfm2 = translate_forces_moments(pfm2, pfm1.point)
    cond = isapprox(pfm1.forces, trans_pfm2.forces; kwargs...) &
           isapprox(pfm1.moments, trans_pfm2.moments; kwargs...)
    return cond
end

"""
    rotate(pfm::PointForcesMoments, psi, theta, phi)

Express the forces and moments in other coordinate system using the Euler angles
that orientate the original pfm with respect to the second.
"""
function rotate(pfm::PointForcesMoments, psi, theta, phi)
    pfm = PointForcesMoments(pfm.point,
                             body2hor(pfm.forces..., psi, theta, phi),
                             body2hor(pfm.moments..., psi, theta, phi)
                             )
end


text:
using GLMakie

function cubedshellwarp(a, b, c, R = max(abs(a), abs(b), abs(c)))

    function f(sR, ξ, η)
        X, Y = tan(π * ξ / 4), tan(π * η / 4)
        x1 = sR / sqrt(X^2 + Y^2 + 1)
        x2, x3 = X * x1, Y * x1
        x1, x2, x3
    end

    fdim = argmax(abs.((a, b, c)))
    if fdim == 1 && a < 0
        # (-R, *, *) : Face I from Ronchi, Iacono, Paolucci (1996)
        x1, x2, x3 = f(-R, b / a, c / a)
    elseif fdim == 2 && b < 0
        # ( *,-R, *) : Face II from Ronchi, Iacono, Paolucci (1996)
        x2, x1, x3 = f(-R, a / b, c / b)
    elseif fdim == 1 && a > 0
        # ( R, *, *) : Face III from Ronchi, Iacono, Paolucci (1996)
        x1, x2, x3 = f(R, b / a, c / a)
    elseif fdim == 2 && b > 0
        # ( *, R, *) : Face IV from Ronchi, Iacono, Paolucci (1996)
        x2, x1, x3 = f(R, a / b, c / b)
    elseif fdim == 3 && c > 0
        # ( *, *, R) : Face V from Ronchi, Iacono, Paolucci (1996)
        x3, x2, x1 = f(R, b / c, a / c)
    elseif fdim == 3 && c < 0
        # ( *, *,-R) : Face VI from Ronchi, Iacono, Paolucci (1996)
        x3, x2, x1 = f(-R, b / c, a / c)
    else
        error("invalid case for cubedshellwarp: $a, $b, $c")
    end

    return x1, x2, x3
end

n = 20

r = range(-1,1,length=n+1)

xa = zeros(n+1,n+1)
xb = zeros(n+1,n+1)
xc = zeros(n+1,n+1)

xa .= r
xb .= r'
xc .= 1

a = [cubedshellwarp(a,b,c)[1] for (a,b,c) in zip(xa,xb,xc)]
b = [cubedshellwarp(a,b,c)[2] for (a,b,c) in zip(xa,xb,xc)]
c = [cubedshellwarp(a,b,c)[3] for (a,b,c) in zip(xa,xb,xc)]

sc = Scene()

wireframe!(sc,
           a,b,c,
           show_axis=false,
           linewidth=1.2)
wireframe!(sc,
           a,b,-c,
           show_axis=false,
           linewidth=1.2)
wireframe!(sc,
           c,a,b,
           show_axis=false,
           linewidth=1.2)
wireframe!(sc,
           -c,a,b,
           show_axis=false,
           linewidth=1.2)
wireframe!(sc,
           b,c,a,
           show_axis=false,
           linewidth=1.2)
wireframe!(sc,
           b,-c,a,
           show_axis=false,
           linewidth=1.2)

sc


text:
"""
A handle to a cursor structure for navigating through a database.
"""
mutable struct Cursor
    handle::Ptr{Nothing}
    Cursor(cur::Ptr{Nothing}) = new(cur)
end

"Check if cursor is open"
isopen(cur::Cursor) = cur.handle != C_NULL

"Create a cursor"
function open(txn::Transaction, dbi::DBI)
    cur_ptr_ref = Ref{Ptr{Nothing}}(C_NULL)
    ret = ccall((:mdb_cursor_open, liblmdb), Cint,
                 (Ptr{Nothing}, Cuint, Ptr{Ptr{Nothing}}),
                  txn.handle, dbi.handle, cur_ptr_ref)
    (ret != 0) && throw(LMDBError(ret))
    return Cursor(cur_ptr_ref[])
end

"Wrapper of Cursor `open` for `do` construct"
function open(f::Function, txn::Transaction, dbi::DBI)
    cur = open(txn, dbi)
    try
        f(cur)
    finally
        close(cur)
    end
end

"Close a cursor"
function close(cur::Cursor)
    if cur.handle == C_NULL
        warn("Cursor is already closed")
    end
    ccall((:mdb_cursor_close, liblmdb), Nothing, (Ptr{Nothing},), cur.handle)
    cur.handle = C_NULL
    return
end

"Renew a cursor"
function renew(txn::Transaction, cur::Cursor)
    ret = ccall((:mdb_cursor_renew, liblmdb), Cint,
                 (Ptr{Nothing}, Ptr{Nothing}), txn.handle, cur.handle)
    (ret != 0) && throw(LMDBError(ret))
    return ret
end

"Return the cursor's transaction"
function transaction(cur::Cursor)
    txn_ptr = ccall((:mdb_cursor_txn, liblmdb), Ptr{Nothing}, (Ptr{Nothing},), cur.handle)
    (txn_ptr == C_NULL) && return nothing
    return Transaction(txn_ptr)
end

"Return the cursor's database"
function database(cur::Cursor)
    dbi = ccall((:mdb_cursor_dbi, liblmdb), Cuint, (Ptr{Nothing},), cur.handle)
    (dbi == 0) && return nothing
    return DBI(dbi, "")
end

"Type to implement the Iterator interface"
mutable struct KeyIterator
   cur::Cursor
   keytype::Type
end

"Iterate over keys"
function Base.iterate(iter::KeyIterator, first=true)
    # Setup parameters
    mdb_key_ref = Ref(MDBValue())
    mdb_val_ref = Ref(MDBValue())
    NOTFOUND::Cint = -30798

    cursor_op = first ? FIRST : NEXT
    ret = ccall( (:mdb_cursor_get, liblmdb), Cint,
               (Ptr{Nothing}, Ptr{MDBValue}, Ptr{MDBValue}, Cint),
                iter.cur.handle, mdb_key_ref, mdb_val_ref, Cint(cursor_op))

    if ret == 0
        # Convert to proper type
        return (convert(iter.keytype, mdb_key_ref), false)
    elseif ret == NOTFOUND
        return nothing
    else
        throw(LMDBError(ret))
    end
end

Base.IteratorSize(::KeyIterator) = Base.SizeUnknown()
Base.eltype(iter::KeyIterator) = iter.keytype

"Return iterator over keys of uniform, specified type"
function keys(cur::Cursor, keytype::Type{T}) where T
    return  KeyIterator(cur, keytype)
end

"""Retrieve by cursor.

This function retrieves key/data pairs from the database.
"""
function get(cur::Cursor, key, ::Type{T}, op::CursorOps=SET_KEY) where T
    # Setup parameters
    k = isbitstype(typeof(key)) ? [key] :  key
    mdb_key_ref = Ref(MDBValue(k))
    mdb_val_ref = Ref(MDBValue())

    # Get value
    ret = ccall( (:mdb_cursor_get, liblmdb), Cint,
                  (Ptr{Nothing}, Ptr{MDBValue}, Ptr{MDBValue}, Cint),
                   cur.handle, mdb_key_ref, mdb_val_ref, Cint(op))
    (ret != 0) && throw(LMDBError(ret))

    # Convert to proper type
    return convert(T, mdb_val_ref)
end

"""Store by cursor.

This function stores key/data pairs into the database. The cursor is positioned at the new item, or on failure usually near it.
"""
function put!(cur::Cursor, key, val; flags::Cuint = zero(Cuint))
    k = isbitstype(typeof(key)) ? [key] :  key
    mdb_key_ref = Ref(MDBValue(k))
    v = isbitstype(typeof(val)) ? [val] :  val
    mdb_val_ref = Ref(MDBValue(v))

    ret = ccall((:mdb_cursor_put, liblmdb), Cint,
                 (Ptr{Nothing}, Ptr{MDBValue}, Ptr{MDBValue}, Cuint),
                  cur.handle, mdb_key_ref, mdb_val_ref, flags)

    (ret != 0) && throw(LMDBError(ret))
    return ret
end

"Delete current key/data pair to which the cursor refers"
function delete!(cur::Cursor; flags::Cuint = zero(Cuint))
    ret = ccall((:mdb_cursor_del, liblmdb), Cint,
                 (Ptr{Nothing}, Cuint), cur.handle, flags)
    (ret != 0) && throw(LMDBError(ret))
    return ret
end

"Return count of duplicates for current key"
function count(cur::Cursor)
    countp = Csize_t[0]
    ret = ccall( (:mdb_cursor_count, liblmdb), Cint,
                  (Ptr{Nothing}, Csize_t), cur.handle, countp)
    (ret != 0) && throw(LMDBError(ret))
    return Int(countp[1])
end


text:
using Kinetic, Plots, Flux
using KitBase.JLD2
using KitBase.ProgressMeter: @showprogress
using BenchmarkTools

function split_regime!(regime, regime0, ks, ctr, nn)
    @inbounds Threads.@threads for i in eachindex(regime)
        if i == 0
            sw = (ctr[i+1].w .- ctr[i].w) / ks.ps.dx[i]
        else
            sw = (ctr[i].w .- ctr[i-1].w) / ks.ps.dx[i]
        end
        tau = vhs_collision_time(ctr[i].prim, ks.gas.μᵣ, ks.gas.ω)
        regime[i] = Int(round(nn([ctr[i].w; sw; tau])[1]))

        if regime[i] == 1 && regime0[i] == 0
            Mu, Mxi, _, _1 = gauss_moments(ctr[i].prim, ks.gas.K)
            a = pdf_slope(ctr[i].prim, sw, ks.gas.K)
            swt = -ctr[i].prim[1] .* moments_conserve_slope(a, Mu, Mxi, 1)
            A = pdf_slope(ctr[i].prim, swt, ks.gas.K)
            tau = vhs_collision_time(ctr[i].prim, ks.gas.μᵣ, ks.gas.ω)
            ctr[i].f .= chapman_enskog(ks.vs.u, ctr[i].prim, a, A, tau)
        end
    end

    regime0 .= regime

    return nothing
end

function up!(ks, ctr, dt, regime)
    res = zeros(3)
    avg = zeros(3)

    @inbounds Threads.@threads for i = 1:ks.ps.nx
        if regime[i] == 0
            KitBase.step!(
                face[i].fw,
                ctr[i].w,
                ctr[i].prim,
                face[i+1].fw,
                ks.gas.γ,
                ks.ps.dx[i],
                res,
                avg,
            )
        else
            KitBase.step!(
                face[i].fw,
                face[i].ff,
                ctr[i].w,
                ctr[i].prim,
                ctr[i].f,
                face[i+1].fw,
                face[i+1].ff,
                ks.vSpace.u,
                ks.vSpace.weights,
                ks.gas.γ,
                ks.gas.μᵣ,
                ks.gas.ω,
                ks.gas.Pr,
                ks.ps.dx[i],
                dt,
                res,
                avg,
                :bgk,
            )
        end
    end
end

cd(@__DIR__)
#@load "../nn_scalar.jld2" nn
@load "/home2/vavrines/Coding/Flowmachine/src/1d/sampler/nn_rif.jld2" nn

set = Setup(case = "sod", space = "1d1f1v", maxTime = 0.15)
ps = PSpace1D(0.0, 1.0, 200, 1)
vs = VSpace1D(-5.0, 5.0, 100)
gas = Gas(Kn = 1e-4, K = 0.0, γ = 3.0)
ib = IB1F(ib_sod(set, ps, vs, gas)...)
ks = SolverSet(set, ps, vs, gas, ib)
ctr, face = init_fvm(ks, ks.ps, :dynamic_array; structarray = true)
for i in eachindex(face)
    face[i].fw .= 0.0
    face[i].ff .= 0.0
end

t = 0.0
dt = timestep(ks, ctr, t)
nt = Int(ks.set.maxTime ÷ dt) + 1
res = zero(ctr[1].w)
avg = zero(res)
regime = ones(Int, axes(ks.ps.x))
regime0 = deepcopy(regime)

@time for iter = 1:nt
    #reconstruct!(ks, ctr)
    split_regime!(regime, regime0, ks, ctr, nn)

    for i in eachindex(face)
        if 1 ∉ (regime[i], regime[i-1])
            flux_gks!(
                face[i].fw,
                ctr[i-1].w .+ ctr[i-1].sw .* ks.ps.dx[i-1] / 2,
                ctr[i].w .- ctr[i].sw .* ks.ps.dx[i] / 2,
                ks.gas.K,
                ks.gas.γ,
                ks.gas.μᵣ,
                ks.gas.ω,
                dt,
                ks.ps.dx[i-1] / 2,
                ks.ps.dx[i] / 2,
                ctr[i-1].sw,
                ctr[i].sw,
            )
        else
            flux_kfvs!(
                face[i].fw,
                face[i].ff,
                ctr[i-1].f .+ ctr[i-1].sf .* ks.ps.dx[i-1] / 2,
                ctr[i].f .- ctr[i].sf .* ks.ps.dx[i] / 2,
                ks.vs.u,
                ks.vs.weights,
                dt,
                ctr[i-1].sf,
                ctr[i].sf,
            )
        end
    end

    up!(ks, ctr, dt, regime)

    global t += dt
end

plot(ks, ctr)

max_stars_repo_name: CSMMLab/Flowmachine

max_issues_repo_name: CSMMLab/Flowmachine

max_forks_repo_name: CSMMLab/Flowmachine


text:
using Oceananigans
using Oceananigans.Units

using Oceananigans.TurbulenceClosures.CATKEVerticalDiffusivities:
    CATKEVerticalDiffusivity,
    SurfaceTKEFlux,
    MixingLength

using OceanTurbulenceParameterEstimation.Observations: FieldTimeSeriesCollector
using JLD2
using GLMakie
using Printf
using Statistics

dir = "six_day_suite"
filename = "free_convection"
xy_filepath = joinpath(dir, filename * "_xy_slice.jld2")

#####
##### Metadata
#####

xy_file = jldopen(xy_filepath)

Nx = xy_file["grid/Nx"]
Ny = xy_file["grid/Ny"]
Nz = xy_file["grid/Nz"]
Hx = xy_file["grid/Hx"]
Hy = xy_file["grid/Hy"]
Hz = xy_file["grid/Hz"]
Lx = xy_file["grid/Lx"]
Ly = xy_file["grid/Ly"]
Lz = xy_file["grid/Lz"]

iterations = parse.(Int, keys(xy_file["timeseries/t"]))
times = [xy_file["timeseries/t/$i"] for i in iterations]
Nt = length(times)

Qᵀ = xy_file["parameters/temperature_flux"]
Qᵘ = xy_file["parameters/momentum_flux"]
f  = xy_file["parameters/coriolis_parameter"]
α  = xy_file["parameters/thermal_expansion_coefficient"]
f  = xy_file["parameters/coriolis_parameter"]
dTdz = xy_file["parameters/dθdz_deep"]

close(xy_file)

H = Hz
chop(a::AbstractArray{<:Any, 1}) = a[1+H:end-H]
chop(a::AbstractArray{<:Any, 2}) = a[1+H:end-H, 1+H:end-H]

function extract_slices(dir, filename; name="w")
    xz_filepath = joinpath(dir, filename * "_xz_slice.jld2")
    yz_filepath = joinpath(dir, filename * "_yz_slice.jld2")
    xy_filepath = joinpath(dir, filename * "_xy_slice.jld2")

    slices = []
    for (filepath, dims) in zip((yz_filepath, xz_filepath, xy_filepath), (1, 2, 3))
        file = jldopen(filepath)
        push!(slices, [chop(dropdims(file["timeseries/$name/$i"]; dims)) for i in iterations])
        close(file)
    end

    return (yz=slices[1], xz=slices[2], xy=slices[3])
end

#####
##### 3D visualization
#####

grid = RectilinearGrid(size=(Nx, Ny, Nz), halo=(Hx, Hy, Hz), x=(0, Lx), y=(0, Ly), z=(-Lz, 0))

e = extract_slices(dir, filename, name="e")
u = extract_slices(dir, filename, name="u")
w = extract_slices(dir, filename, name="w")
T = extract_slices(dir, filename, name="T")
ϵ = extract_slices(dir, filename, name="ϵ")

x, y, z = nodes((Center, Center, Center), grid)

x_xz = repeat(x, 1, Nz)
y_xz = 0.995 * Ly * ones(Nx, Nz)
z_xz = repeat(reshape(z, 1, Nz), Nx, 1)

x_yz = 0.995 * Lx * ones(Ny, Nz)
y_yz = repeat(y, 1, Nz)
z_yz = repeat(reshape(z, 1, Nz), grid.Ny, 1)

# Slight displacements to "stitch" the cube together
x_xy = x
y_xy = y
z_xy = - 0.001 * Lz * ones(grid.Nx, grid.Ny)

#####
##### Animate!
#####

fig = Figure(resolution=(1200, 1200))

ax_e = Axis3(fig[1, 1:4], aspect=:data,
             xlabel="x (m)", ylabel="y (m)", zlabel="z (m)",
             azimuth=1.26, elevation=0.53, perspectiveness=0.5)

hidedecorations!(ax_e)
hidespines!(ax_e)

slider = Slider(fig[1, 5], range=1:Nt, horizontal=false, startvalue=1)

n = slider.value
#n = Observable(1)

###
### Vertical velocity
###

w_yz = @lift w.yz[$n]
w_xz = @lift w.xz[$n]
w_xy = @lift w.xy[$n]

e_yz = @lift e.yz[$n]
e_xz = @lift e.xz[$n]
e_xy = @lift e.xy[$n]

ϵ_yz = @lift log10.(ϵ.yz[$n])
ϵ_xz = @lift log10.(ϵ.xz[$n])
ϵ_xy = @lift log10.(ϵ.xy[$n])

u_yz = @lift u.yz[$n]
u_xz = @lift u.xz[$n]
u_xy = @lift u.xy[$n]

T_yz = @lift T.yz[$n]
T_xz = @lift T.xz[$n]
T_xy = @lift T.xy[$n]

colormap_T = :oslo
colormap_u = :balance
colormap_w = :balance
colormap_e = :solar
colormap_ϵ = :viridis

w_max = maximum(maximum.(abs, w.xz))
colorrange_w = (-w_max/2, w_max/2)

u_max = maximum(maximum.(abs, u.xz))
colorrange_u = (-u_max/2, u_max/2)

e_max = maximum(maximum.(abs, e.xz))
colorrange_e = (0, e_max/6)

# Compute maximum temperature for colorrange with a moving average
ΔT = 0.04 # Absolute range
Δn = 2 # Half-width of moving average
Tf_max = maximum(T.xy[end])
T0_max = maximum(T.xy[1])

colorrange_T = @lift begin
    if $n < Δn + 1
        T_max = T0_max
    elseif $n < Nt - Δn
        T_max = mean([maximum(T.xy[nn]) for nn in $n-Δn:$n+Δn])
    else
        T_max = Tf_max
    end

    T_min = T_max - ΔT

    (T_min, T_max)
end

function box!(ax, yz, xz, xy; colormap, colorrange)
    sfc = surface!(ax, x_xz, y_xz, z_xz; color=xz, colormap, colorrange)
          surface!(ax, x_yz, y_yz, z_yz; color=yz, colormap, colorrange)
          surface!(ax, x_xy, y_xy, z_xy; color=xy, colormap, colorrange)
    return sfc
end

sfc_e = box!(ax_e, e_yz, e_xz, e_xy, colormap=colormap_e, colorrange=colorrange_e)

cp_e = fig[2, 2:3] = Colorbar(fig, sfc_e,
                              flipaxis = false,
                              vertical = false,
                              label = "Turbulent kinetic energy (m² s⁻²)",
                              textsize = 24)

title = @lift begin
    t_hours = times[$n] / hours
    @sprintf("Convecting ocean surface boundary layer after %.1f hours", t_hours)
end

lab = Label(fig[0, 2:3], title, textsize=28)

display(fig)

#=
record(fig, joinpath(dir, filename * "_LES_one_field.mp4"), 1:Nt; framerate=12) do nn
    @info "Drawing frame $nn of $Nt..."
    n[] = nn
end
=#

max_stars_repo_name: glwagner/LocalOceanClosureCalibration


text:
using JuliaFormatter

include("clima_formatter_options.jl")

format(@__FILE__; clima_formatter_options...)


text:
"""
  nn.jl File

Neural Network implementation (Module BetaML.Nn)

`?BetaML.Nn` for documentation

- [Importable source code (most up-to-date version)](https://github.com/sylvaticus/BetaML.jl/blob/master/src/Nn.jl) - [Julia Package](https://github.com/sylvaticus/BetaML.jl)
- [Demonstrative static notebook](https://github.com/sylvaticus/BetaML.jl/blob/master/notebooks/Nn.ipynb)
- [Demonstrative live notebook](https://mybinder.org/v2/gh/sylvaticus/BetaML.jl/master?filepath=notebooks%2FNn.ipynb) (temporary personal online computational environment on myBinder) - it can takes minutes to start with!
- Theory based on [MITx 6.86x - Machine Learning with Python: from Linear Models to Deep Learning](https://github.com/sylvaticus/MITx_6.86x) ([Unit 3](https://github.com/sylvaticus/MITx_6.86x/blob/master/Unit%2003%20-%20Neural%20networks/Unit%2003%20-%20Neural%20networks.md))
- New to Julia? [A concise Julia tutorial](https://github.com/sylvaticus/juliatutorial) - [Julia Quick Syntax Reference book](https://julia-book.com)

"""

"""
    BetaML.Nn module

Implement the functionality required to define an artificial Neural Network, train it with data, forecast data and assess its performances.

Common type of layers and optimisation algorithms are already provided, but you can define your own ones subclassing respectively the `AbstractLayer` and `OptimisationAlgorithm` abstract types.

The module provide the following type or functions. Use `?[type or function]` to access their full signature and detailed documentation:

# Model definition:

- `DenseLayer`: Classical feed-forward layer with user-defined activation function
- `DenseNoBiasLayer`: Classical layer without the bias parameter
- `VectorFunctionLayer`: Parameterless layer whose activation function run over the ensable of its nodes rather than on each one individually
- `buildNetwork`: Build the chained network and define a cost function
- `getParams(nn)`: Retrieve current weigthts
- `getGradient(nn)`: Retrieve the current gradient of the weights
- `setParams!(nn)`: Update the weigths of the network
- `show(nn)`: Print a representation of the Neural Network

Each layer can use a default activation function, one of the functions provided in the `Utils` module (`relu`, `tanh`, `softmax`,...) or you can specify your own function. The derivative of the activation function can be optionally be provided, in such case training will be quicker, altought this difference tends to vanish with bigger datasets.
You can alternativly implement your own layer defining a new type as subtype of the abstract type `AbstractLayer`. Each user-implemented layer must define the following methods:

- A suitable constructor
- `forward(layer,x)`
- `backward(layer,x,nextGradient)`
- `getParams(layer)`
- `getGradient(layer,x,nextGradient)`
- `setParams!(layer,w)`
- `size(layer)`

# Model training:

- `trainingInfo(nn)`: Default callback function during training
- `train!(nn)`:  Training function
- `singleUpdate!(θ,▽;optAlg)`: The parameter update made by the specific optimisation algorithm
- `SGD`: The default optimisation algorithm
- `ADAM`: A faster moment-based optimisation algorithm (added in v0.2.2)

To define your own optimisation algorithm define a subtype of `OptimisationAlgorithm` and implement the function `singleUpdate!(θ,▽;optAlg)` and eventually `initOptAlg(⋅)` specific for it.

# Model predictions and assessment:

- `predict(nn)`: Return the output given the data
- `loss(nn)`: Compute avg. network loss on a test set
- `Utils.accuracy(ŷ,y)`: Categorical output accuracy

While high-level functions operating on the dataset expect it to be in the standard format (nRecords × nDimensions matrices) it is custom to represent the chain of a neural network as a flow of column vectors, so all low-level operations (operating on a single datapoint) expect both the input and the output as a column vector.
"""
module Nn

#import Base.Threads.@spawn

using Random, Zygote, ProgressMeter, Reexport
import Distributions: Uniform

using ForceImport
@force using ..Api
@force using ..Utils

import Base.size
import Base: +, -, *, /, sum, sqrt

# module own functions
export AbstractLayer, forward, backward, getParams, getNParams, getGradient, setParams!, size, NN,
       buildNetwork, predict, loss, train!, getindex, initOptAlg!, singleUpdate!,
       DenseLayer, DenseNoBiasLayer, VectorFunctionLayer,
       Learnable,
       show

# for working on gradient as e.g [([1.0 2.0; 3.0 4.0], [1.0,2.0,3.0]),([1.0,2.0,3.0],1.0)]
"""
Learnable(data)

Structure representing the learnable parameters of a layer or its gradient.

The learnable parameters of a layers are given in the form of a N-tuple of Array{Float64,N2} where N2 can change (e.g. we can have a layer with the first parameter being a matrix, and the second one being a scalar).
We wrap the tuple on its own structure a bit for some efficiency gain, but above all to define standard mathematic operations on the gradients without doing "type pyracy" with respect to Base tuples.
"""
mutable struct Learnable
    data::Tuple{Vararg{Array{Float64,N} where N}}
    function Learnable(data)
        return new(data)
    end
end
function +(items::Learnable...)
  values = collect(items[1].data)
  N = length(values)
  @inbounds for item in items[2:end]
       @inbounds  @simd for n in 1:N # @inbounds  @simd
          values[n] += item.data[n]
      end
  end
  return Learnable(Tuple(values))
end
sum(items::Learnable...)  = +(items...)
function -(items::Learnable...)
  values = collect(items[1].data)
  N = length(values)
 @inbounds for item in items[2:end]
       @inbounds @simd for n in 1:N # @simd
          values[n] -= item.data[n]
      end
  end
  return Learnable(Tuple(values))
end
function *(items::Learnable...)
  values = collect(items[1].data)
  N = length(values)
  @inbounds for item in items[2:end]
      @inbounds @simd for n in 1:N # @simd
          values[n] = values[n] .* item.data[n]
      end
  end
  return Learnable(Tuple(values))
end
+(item::Learnable,sc::Number) =  Learnable(Tuple([item.data[i] .+ sc for i in 1:length(item.data)]))
+(sc::Number, item::Learnable) = +(item,sc)
-(item::Learnable,sc::Number) = Learnable(Tuple([item.data[i] .- sc for i in 1:length(item.data)]))
-(sc::Number, item::Learnable) = (-(item,sc)) * -1
*(item::Learnable,sc::Number) = Learnable(item.data .* sc)
*(sc::Number, item::Learnable) = Learnable(sc .* item.data)
/(item::Learnable,sc::Number) = Learnable(item.data ./ sc)
/(sc::Number,item::Learnable,) = Learnable(Tuple([sc ./ item.data[i] for i in 1:length(item.data)]))
sqrt(item::Learnable) = Learnable(Tuple([sqrt.(item.data[i]) for i in 1:length(item.data)]))
/(item1::Learnable,item2::Learnable) = Learnable(Tuple([item1.data[i] ./ item2.data[i] for i in 1:length(item1.data)]))

#=
# not needed ??
function Base.iterate(iter::Learnable, state=(iter.data[1], 1))
           element, count = state
           if count > length(iter)
               return nothing
           elseif count == length(iter)
               return (element, (iter.data[count], count + 1))
           end
           return (element, (iter.data[count+1], count + 1))
end
Base.length(iter::Learnable) = length(iter.data)
#Base.eltype(iter::Learnable) = Int
=#

## Sckeleton for the layer functionality.
# See nn_default_layers.jl for actual implementations

abstract type AbstractLayer end

include("Nn_default_layers.jl")

"""
   forward(layer,x)

Predict the output of the layer given the input

# Parameters:
* `layer`:  Worker layer
* `x`:      Input to the layer

# Return:
- An Array{T,1} of the prediction (even for a scalar)
"""
function forward(layer::AbstractLayer,x)
 error("Not implemented for this kind of layer. Please implement `forward(layer,x)`.")
end

"""
   backward(layer,x,nextGradient)

Compute backpropagation for this layer

# Parameters:
* `layer`:        Worker layer
* `x`:            Input to the layer
* `nextGradient`: Derivative of the overaall loss with respect to the input of the next layer (output of this layer)

# Return:
* The evaluated gradient of the loss with respect to this layer inputs

"""
function backward(layer::AbstractLayer,x,nextGradient)
    error("Not implemented for this kind of layer. Please implement `backward(layer,x,nextGradient)`.")
end

"""
   getParams(layer)

Get the layers current value of its trainable parameters

# Parameters:
* `layer`:  Worker layer

# Return:
* The current value of the layer's trainable parameters as tuple of matrices. It is up to you to decide how to organise this tuple, as long you are consistent with the `getGradient()` and `setParams()` functions. Note that starting from BetaML 0.2.2 this tuple needs to be wrapped in its `Learnable` type.
"""
function getParams(layer::AbstractLayer)
  error("Not implemented for this kind of layer. Please implement `getParams(layer)`.")
end

"""
   getGradient(layer,x,nextGradient)

Compute backpropagation for this layer

# Parameters:
* `layer`:        Worker layer
* `x`:            Input to the layer
* `nextGradient`: Derivative of the overaall loss with respect to the input of the next layer (output of this layer)

# Return:
* The evaluated gradient of the loss with respect to this layer's trainable parameters as tuple of matrices. It is up to you to decide how to organise this tuple, as long you are consistent with the `getParams()` and `setParams()` functions. Note that starting from BetaML 0.2.2 this tuple needs to be wrapped in its `Learnable` type.
"""
function getGradient(layer::AbstractLayer,x,nextGradient)
    error("Not implemented for this kind of layer. Please implement `getGradient(layer,x,nextGradient)`.")
  end

"""
     setParams!(layer,w)

Set the trainable parameters of the layer with the given values

# Parameters:
* `layer`: Worker layer
* `w`:     The new parameters to set (Learnable)

# Notes:
*  The format of the tuple wrapped by Learnable must be consistent with those of the `getParams()` and `getGradient()` functions.
"""
function setParams!(layer::AbstractLayer,w)
    error("Not implemented for this kind of layer. Please implement `setParams!(layer,w)`.")
end

"""
    size(layer)

Get the dimensions of the layers in terms of (dimensions in input , dimensions in output)

# Notes:
* You need to use `import Base.size` before defining this function for your layer
"""
function size(layer::AbstractLayer)
    error("Not implemented for this kind of layer. Please implement `size(layer)`.")
end

"""getNParams(layer)

Return the number of parameters of a layer.

It doesn't need to be implemented by each layer type, as it uses getParams().
"""
function getNParams(layer::AbstractLayer)
    pars = getParams(layer)
    nP = 0
    for p in pars.data
        nP += *(size(p)...)
    end
    return nP
end

# ------------------------------------------------------------------------------
# NN-related functions
"""
   NN

Representation of a Neural Network

# Fields:
* `layers`:  Array of layers objects
* `cf`:      Cost function
* `dcf`:     Derivative of the cost function
* `trained`: Control flag for trained networks
"""
mutable struct NN
    layers::Array{AbstractLayer,1}
    cf::Function
    dcf::Union{Function,Nothing}
    trained::Bool
    name::String
end

"""
   buildNetwork(layers,cf;dcf,name)

Instantiate a new Feedforward Neural Network

Parameters:
* `layers`: Array of layers objects
* `cf`:     Cost function
* `dcf`:    Derivative of the cost function [def: `nothing`]
* `name`:   Name of the network [def: "Neural Network"]

# Notes:
* Even if the network ends with a single output note, the cost function and its derivative should always expect y and ŷ as column vectors.
"""
function buildNetwork(layers,cf;dcf=nothing,name="Neural Network")
    return NN(layers,cf,dcf,false,name)
end

"""
   predict(nn,x)

Network predictions

# Parameters:
* `nn`:  Worker network
* `x`:   Input to the network (n × d)
"""
#=
function predict(nn::NN,x)
    makeColVector(x)
    values = x
    for l in nn.layers
        values = forward(l,values)
    end
    return values
end
=#

function predict(nn::NN,x)
    x = makeMatrix(x)
    # get the output dimensions
    n = size(x)[1]
    d = size(nn.layers[end])[2]
    out = zeros(n,d)
    for i in 1:size(x)[1]
        values = x[i,:]
        for l in nn.layers
            values = forward(l,values)
        end
        out[i,:] = values
    end
    return out
end

"""
   loss(fnn,x,y)

Compute avg. network loss on a test set (or a single (1 × d) data point)

# Parameters:
* `fnn`: Worker network
* `x`:   Input to the network (n) or (n x d)
* `y`:   Label input (n) or (n x d)
"""
function loss(nn::NN,x,y)
    x = makeMatrix(x)
    y = makeMatrix(y)
    (n,d) = size(x)
    #(nn.trained || n == 1) ? "" : @warn "Seems you are trying to test a neural network that has not been tested. Use first `train!(nn,x,y)`"
    ϵ = 0.0
    for i in 1:n
        ŷ = predict(nn,x[i,:]')[1,:]
        ϵ += nn.cf(ŷ,y[i,:])
    end
    return ϵ/n
end

"""
   getParams(nn)

Retrieve current weigthts

# Parameters:
* `nn`: Worker network

# Notes:
* The output is a vector of tuples of each layer's input weigths and bias weigths
"""
@inline function getParams(nn::NN)
  return [getParams(l) for l in nn.layers]
end

"""
   getGradient(nn,x,y)

Retrieve the current gradient of the weigthts (i.e. derivative of the cost with respect to the weigths)

# Parameters:
* `nn`: Worker network
* `x`:   Input to the network (d,1)
* `y`:   Label input (d,1)

#Notes:
* The output is a vector of tuples of each layer's input weigths and bias weigths
"""
function getGradient(nn::NN,x::Union{T,AbstractArray{T,1}},y::Union{T2,AbstractArray{T2,1}}) where { T <: Number, T2 <: Number}

  x = makeColVector(x)
  y = makeColVector(y)

  nLayers = length(nn.layers)

  # Stap 1: Forward pass
  forwardStack = Vector{Vector{Float64}}(undef,nLayers+1)

  forwardStack[1] = x
  @inbounds for (i,l) in enumerate(nn.layers)
      forwardStack[i+1] = forward(l,forwardStack[i])
  end

  # Step 2: Backpropagation pass
  backwardStack = Vector{Vector{Float64}}(undef,nLayers+1)
  if nn.dcf != nothing
    backwardStack[end] = nn.dcf(forwardStack[end],y) # adding dϵ_dHatY
  else
    backwardStack[end] = gradient(nn.cf,forwardStack[end],y)[1] # using AD from Zygote
  end
  @inbounds for lidx in nLayers:-1:1
     l = nn.layers[lidx]
     dϵ_do = backward(l,forwardStack[lidx],backwardStack[lidx+1])
     backwardStack[lidx] = dϵ_do
  end

  # Step 3: Computing gradient of weigths
  dWs = Array{Learnable,1}(undef,nLayers)
  @inbounds for lidx in 1:nLayers
     dWs[lidx] = getGradient(nn.layers[lidx],forwardStack[lidx],backwardStack[lidx+1])
  end

  return dWs
end

"""
   getGradient(nn,xbatch,ybatch)

Retrieve the current gradient of the weigthts (i.e. derivative of the cost with respect to the weigths)

# Parameters:
* `nn`:      Worker network
* `xbatch`:  Input to the network (n,d)
* `ybatch`:  Label input (n,d)

#Notes:
* The output is a vector of tuples of each layer's input weigths and bias weigths
"""
function getGradient(nn,xbatch::AbstractArray{T,2},ybatch::AbstractArray{T2,2}) where {T <: Number, T2 <: Number}
    #return [getGradient(nn,xbatch[j,:],ybatch[j,:]) for j in 1:size(xbatch,1)]
    bSize = size(xbatch,1)
    gradients = Array{Vector{Learnable},1}(undef,bSize)
    # Note: in Julia 1.6 somehow the multithreading is less efficient than in Julia 1.5
    # Using @inbounds @simd result faster than using 4 threads, so reverting to it.
    # But to keep following the evolution, as there seems to be some issues on performances
    # in Julia 1.6: https://discourse.julialang.org/t/drop-of-performances-with-julia-1-6-0-for-interpolationkernels/58085
    # Maybe when that's solved it will be again more convenient to use multi-threading
    #Threads.@threads
    @inbounds  for j in 1:bSize # @simd
       gradients[j] =  getGradient(nn,xbatch[j,:],ybatch[j,:])
    end
    return gradients
end

"""
   setParams!(nn,w)

Update weigths of the network

# Parameters:
* `nn`: Worker network
* `w`:  The new weights to set
"""
function setParams!(nn::NN,w)
    for lidx in 1:length(nn.layers)
        setParams!(nn.layers[lidx],w[lidx])
    end
end

"""
  show(nn)

Print a representation of the Neural Network (layers, dimensions..)

# Parameters:
* `nn`: Worker network
"""
function show(nn::NN)
  trainedString = nn.trained == true ? "trained" : "non trained"
  println("*** $(nn.name) ($(length(nn.layers)) layers, $(trainedString))\n")
  println("#\t # In \t # Out \t Type")
  for (i,l) in enumerate(nn.layers)
    shapes = size(l)
    println("$i \t $(shapes[1]) \t\t $(shapes[2]) \t\t $(typeof(l)) ")
  end
end

"getNParams(nn) - Return the number of trainable parameters of the neural network."
function getNParams(nn::NN)
    nP = 0
    for l in nn.layers
        nP += getNParams(l)
    end
    return nP
end

Base.getindex(n::NN, i::AbstractArray) = NN(n.layers[i]...)

# ------------------------------------------------------------------------------
# Optimisation-related functions

"""
    OptimisationAlgorithm

Abstract type representing an Optimisation algorithm.

Currently supported algorithms:
- `SGD` (Stochastic) Gradient Descent

See `?[Name OF THE ALGORITHM]` for their details

You can implement your own optimisation algorithm using a subtype of `OptimisationAlgorithm` and implementing its constructor and the update function `singleUpdate(⋅)` (type `?singleUpdate` for details).

"""
abstract type OptimisationAlgorithm end

include("Nn_default_optalgs.jl")

"""
   trainingInfo(nn,x,y;n,batchSize,epochs,verbosity,nEpoch,nBatch)

Default callback funtion to display information during training, depending on the verbosity level

# Parameters:
* `nn`: Worker network
* `x`:  Batch input to the network (batchSize,d)
* `y`:  Batch label input (batchSize,d)
* `n`: Size of the full training set
* `nBatches` : Number of baches per epoch
* `epochs`: Number of epochs defined for the training
* `verbosity`: Verbosity level defined for the training (NONE,LOW,STD,HIGH,FULL)
* `nEpoch`: Counter of the current epoch
* `nBatch`: Counter of the current batch

#Notes:
* Reporting of the error (loss of the network) is expensive. Use `verbosity=NONE` for better performances
"""
function trainingInfo(nn,x,y;n,nBatches,epochs,verbosity,nEpoch,nBatch)
   if verbosity == NONE
       return false # doesn't stop the training
   end

   nMsgDict = Dict(LOW => 0, STD => 10,HIGH => 100, FULL => n)
   nMsgs = nMsgDict[verbosity]
   batchSize = size(x,1)

   if verbosity == FULL || ( nBatch == nBatches && ( nEpoch == 1  || nEpoch % ceil(epochs/nMsgs) == 0))

      ϵ = loss(nn,x,y)
      println("Training.. \t avg ϵ on (Epoch $nEpoch Batch $nBatch): \t $(ϵ)")
   end
   return false
end

"""
   train!(nn,x,y;epochs,batchSize,sequential,optAlg,verbosity,cb)

Train a neural network with the given x,y data

# Parameters:
* `nn`:         Worker network
* `x`:          Training input to the network (records x dimensions)
* `y`:          Label input (records x dimensions)
* `epochs`:     Number of passages over the training set [def: `100`]
* `batchSize`:  Size of each individual batch [def: `min(size(x,1),32)`]
* `sequential`: Wether to run all data sequentially instead of random [def: `false`]
* `optAlg`:     The optimisation algorithm to update the gradient at each batch [def: `ADAM()`]
* `verbosity`:  A verbosity parameter for the trade off information / efficiency [def: `STD`]
* `cb`:         A callback to provide information. [def: `trainingInfo`]
* `rng`:        Random Number Generator (see [`FIXEDSEED`](@ref)) [deafult: `Random.GLOBAL_RNG`]

# Return:
- A named tuple with the following information
  - `epochs`: Number of epochs actually ran
  - `ϵ_epochs`: The average error on each epoch (if `verbosity > LOW`)
  - `θ_epochs`: The parameters at each epoch (if `verbosity > STD`)

# Notes:
- Currently supported algorithms:
    - `SGD`, the classical (Stochastic) Gradient Descent optimiser
    - `ADAM`,  an adaptive moment estimation optimiser
- Look at the individual optimisation algorithm (`?[Name OF THE ALGORITHM]`) for info on its parameter, e.g. [`?SGD`](@ref SGD) for the Stochastic Gradient Descent.
- You can implement your own optimisation algorithm using a subtype of `OptimisationAlgorithm` and implementing its constructor and the update function `singleUpdate!(⋅)` (type `?singleUpdate!` for details).
- You can implement your own callback function, altought the one provided by default is already pretty generic (its output depends on the `verbosity` parameter). See [`trainingInfo`](@ref) for informations on the cb parameters.
- Both the callback function and the [`singleUpdate!`](@ref) function of the optimisation algorithm can be used to stop the training algorithm, respectively returning `true` or `stop=true`.
- The verbosity can be set to any of `NONE`,`LOW`,`STD`,`HIGH`,`FULL`.
- The update is done computing the average gradient for each batch and then calling `singleUpdate!` to let the optimisation algorithm perform the parameters update
"""
function train!(nn::NN,x,y; epochs=100, batchSize=min(size(x,1),32), sequential=false, verbosity::Verbosity=STD, cb=trainingInfo, optAlg::OptimisationAlgorithm=ADAM(),rng = Random.GLOBAL_RNG)#,   η=t -> 1/(1+t), λ=1, rShuffle=true, nMsgs=10, tol=0optAlg::SD=SD())
    if verbosity > STD
        @codeLocation
    end
    x = makeMatrix(x)
    y = makeMatrix(y)
    (n,d)     = size(x)
    batchSize = min(size(x,1),batchSize)
    if verbosity > NONE # Note that are two "Verbosity type" objects. To compare with numbers use Int(NONE) > 1
        println("***\n*** Training $(nn.name) for $epochs epochs with algorithm $(typeof(optAlg)).")
    end
    ϵ_epoch_l = Inf
    θ_epoch_l = getParams(nn)
    ϵ_epoch   = loss(nn,x,y)
    θ_epoch   = getParams(nn)
    ϵ_epochs  = Float64[]
    θ_epochs  = []

    initOptAlg!(optAlg::OptimisationAlgorithm;θ=getParams(nn),batchSize=batchSize,x=x,y=y)

    timetoShowProgress = verbosity > NONE ? 1 : typemax(Int64)
    @showprogress timetoShowProgress "Training the Neural Network..."    for t in 1:epochs
       batches = batch(n,batchSize,sequential=sequential,rng=rng)
       nBatches = length(batches)
       if t == 1
           if (verbosity >= STD) push!(ϵ_epochs,ϵ_epoch); end
           if (verbosity > STD) push!(θ_epochs,θ_epoch); end
       end
       for (i,batch) in enumerate(batches)
           xbatch = x[batch, :]
           ybatch = y[batch, :]
           θ   = getParams(nn)
           # remove @spawn and fetch (on next row) to get single thread code
           # note that there is no random number issue here..
           #gradients   = @spawn getGradient(nn,xbatch,ybatch)
           #sumGradient = sum(fetch(gradients))
           gradients   = getGradient(nn,xbatch,ybatch)
           sumGradient = sum(gradients)

           ▽   = sumGradient / length(batch)
           #▽   = gradDiv.(gradSum([getGradient(nn,xbatch[j,:],ybatch[j,:]) for j in 1:batchSize]), batchSize)
           res = singleUpdate!(θ,▽;nEpoch=t,nBatch=i,nBatches=nBatches,xbatch=xbatch,ybatch=ybatch,optAlg=optAlg)
           setParams!(nn,res.θ)
           cbOut = cb(nn,xbatch,ybatch,n=d,nBatches=nBatches,epochs=epochs,verbosity=verbosity,nEpoch=t,nBatch=i)
           if(res.stop==true || cbOut==true)
               nn.trained = true
               return (epochs=t,ϵ_epochs=ϵ_epochs,θ_epochs=θ_epochs)
           end
       end
       if (verbosity >= STD)
           ϵ_epoch_l = ϵ_epoch
           ϵ_epoch = loss(nn,x,y)
           push!(ϵ_epochs,ϵ_epoch);
       end
       if (verbosity > STD)
           θ_epoch_l = θ_epoch
           θ_epoch = getParams(nn)
           push!(θ_epochs,θ_epoch); end
    end

    if (verbosity > NONE)
        if verbosity == LOW
            ϵ_epoch = loss(nn,x,y)
        end
        println("Training of $epochs epoch completed. Final epoch error: $(ϵ_epoch).");
     end
    nn.trained = true
    return (epochs=epochs,ϵ_epochs=ϵ_epochs,θ_epochs=θ_epochs)
end

"""
   singleUpdate!(θ,▽;nEpoch,nBatch,batchSize,xbatch,ybatch,optAlg)

Perform the parameters update based on the average batch gradient.

# Parameters:
- `θ`:         Current parameters
- `▽`:         Average gradient of the batch
- `nEpoch`:    Count of current epoch
- `nBatch`:    Count of current batch
- `nBatches`:  Number of batches per epoch
- `xbatch`:    Data associated to the current batch
- `ybatch`:    Labels associated to the current batch
- `optAlg`:    The Optimisation algorithm to use for the update

# Notes:
- This function is overridden so that each optimisation algorithm implement their
own version
- Most parameters are not used by any optimisation algorithm. They are provided
to support the largest possible class of optimisation algorithms
- Some optimisation algorithms may change their internal structure in this function
"""
function singleUpdate!(θ,▽;nEpoch,nBatch,nBatches,xbatch,ybatch,optAlg::OptimisationAlgorithm)
   return singleUpdate!(θ,▽,optAlg;nEpoch=nEpoch,nBatch=nBatch,nBatches=nBatches,xbatch=xbatch,ybatch=ybatch)
end

function singleUpdate!(θ,▽,optAlg::OptimisationAlgorithm;nEpoch,nBatch,nBatches,xbatch,ybatch)
    error("singleUpdate() not implemented for this optimisation algorithm")
end

"""
   initOptAlg!(optAlg;θ,batchSize,x,y)

Initialize the optimisation algorithm

# Parameters:
- `optAlg`:    The Optimisation algorithm to use
- `θ`:         Current parameters
- `batchSize`:    The size of the batch
- `x`:   The training (input) data
- `y`:   The training "labels" to match
* `rng`: Random Number Generator (see [`FIXEDSEED`](@ref)) [deafult: `Random.GLOBAL_RNG`]

# Notes:
- Only a few optimizers need this function and consequently ovverride it. By default it does nothing, so if you want write your own optimizer and don't need to initialise it, you don't have to override this method
"""
initOptAlg!(optAlg::OptimisationAlgorithm;θ,batchSize,x,y,rng = Random.GLOBAL_RNG) = nothing

#=
        if rShuffle
           # random shuffle x and y
           ridx = shuffle(1:size(x)[1])
           x = x[ridx, :]
           y = y[ridx , :]
        end
        ϵ = 0
        #η = dyn_η ? 1/(1+t) : η
        ηₜ = η(t)*λ
        for i in 1:size(x)[1]
            xᵢ = x[i,:]'
            yᵢ = y[i,:]'
            W  = getParams(nn)
            dW = getGradient(nn,xᵢ,yᵢ)
            newW = gradientDescentSingleUpdate(W,dW,ηₜ)
            setParams!(nn,newW)
            ϵ += loss(nn,xᵢ,yᵢ)
        end
        if nMsgs != 0 && (t % ceil(maxEpochs/nMsgs) == 0 || t == 1 || t == maxEpochs)
          println("Avg. error after epoch $t : $(ϵ/size(x)[1])")
        end

        if abs(ϵl/size(x)[1] - ϵ/size(x)[1]) < (tol * abs(ϵl/size(x)[1]))
            if nMsgs != 0
                println((tol * abs(ϵl/size(x)[1])))
                println("*** Avg. error after epoch $t : $(ϵ/size(x)[1]) (convergence reached")
            end
            converged = true
            break
        else
            ϵl = ϵ
        end
    end
    if nMsgs != 0 && converged == false
        println("*** Avg. error after epoch $maxEpochs : $(ϵ/size(x)[1]) (convergence not reached)")
    end
    nn.trained = true
end

 =#

end # end module


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, Pkg

name = "SCIP"
version = v"0.1.3"

# Collection of sources required to complete build
sources = [
    ArchiveSource("https://scip.zib.de/download/release/scipoptsuite-7.0.3.tgz", "5af5185a6e60cc62d1a89e3ac4fe22d32351a5158c2c04a95e180e76eb98cc07"),
]

# bliss script, ignored for now
script_bliss = raw"""
cd $WORKSPACE/srcdir
cd bliss
mkdir build
cd build
cmake -DCMAKE_INSTALL_PREFIX=$prefix -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} -DCMAKE_BUILD_TYPE=Release ..
make -j${nproc}
make install
cd $WORKSPACE/srcdir
"""

# Bash recipe for building across all platforms
script = raw"""
cd scipoptsuite*
mkdir build
cd build/
cmake -DCMAKE_INSTALL_PREFIX=$prefix\
  -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN}\
  -DCMAKE_BUILD_TYPE=Release\
  -DPAPILO=0\
  -DZIMPL=OFF\
  -DGCG=0\
  -DREADLINE=OFF\
  -DBOOST=off\
  -DSYM=bliss\
  -DIPOPT_DIR=${prefix} -DIPOPT_LIBRARIES=${libdir} ..
make -j${nproc}
make install

mkdir -p ${prefix}/share/licenses/SCIP
for dir in papilo scip soplex; do
    cp $WORKSPACE/srcdir/scipoptsuite*/${dir}/COPYING ${prefix}/share/licenses/SCIP/LICENSE_${dir}
done
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    Platform("i686", "linux"; libc="glibc"),
    Platform("x86_64", "linux"; libc="glibc"),
    Platform("aarch64", "linux"; libc="glibc"),
    Platform("armv7l", "linux"; libc="glibc"),
    Platform("powerpc64le", "linux"; libc="glibc"),
    Platform("i686", "linux"; libc="musl"),
    Platform("x86_64", "linux"; libc="musl"),
    Platform("aarch64", "linux"; libc="musl"),
    Platform("armv7l", "linux"; libc="musl"),
    Platform("x86_64", "macos"),
    Platform("x86_64", "freebsd"),
    Platform("i686", "windows"),
    Platform("x86_64", "windows"),
]

platforms = expand_cxxstring_abis(platforms)

# The products that we will ensure are always built
products = [
    LibraryProduct("libscip", :libscip),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency(PackageSpec(name="bliss_jll", uuid="508c9074-7a14-5c94-9582-3d4bc1871065")),
    Dependency(PackageSpec(name="GMP_jll", uuid="781609d7-10c4-51f6-84f2-b8444358ff6d"), v"6.1.2"),
    Dependency(PackageSpec(name="CompilerSupportLibraries_jll", uuid="e66e0078-7015-5450-92f7-15fbd957f2ae")),
    Dependency(PackageSpec(name="Ipopt_jll", uuid="9cc047cb-c261-5740-88fc-0cf96f7bdcc7")),
    Dependency(PackageSpec(name="Zlib_jll", uuid="83775a58-1f1d-513f-b197-d71354ab007a")),
]

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

max_stars_repo_name: jameskermode/Yggdrasil

max_issues_repo_name: jameskermode/Yggdrasil

max_forks_repo_name: jameskermode/Yggdrasil


text:
# ---
# title: 1621. Number of Sets of K Non-Overlapping Line Segments
# id: problem1621
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Dynamic Programming
# link: <https://leetcode.com/problems/number-of-sets-of-k-non-overlapping-line-segments/description/>
# hidden: true
# ---
#
# Given `n` points on a 1-D plane, where the `ith` point (from `0` to `n-1`) is
# at `x = i`, find the number of ways we can draw **exactly** `k` **non-
# overlapping** line segments such that each segment covers two or more points.
# The endpoints of each segment must have **integral coordinates**. The `k` line
# segments **do not** have to cover all `n` points, and they are **allowed** to
# share endpoints.
#
# Return _the number of ways we can draw_`k` _non-overlapping line segments_ _._
# Since this number can be huge, return it **modulo** `109 + 7`.
#
#
#
# **Example 1:**
#
# ![](https://assets.leetcode.com/uploads/2020/09/07/ex1.png)
#
#
#
#     Input: n = 4, k = 2
#     Output: 5
#     Explanation: The two line segments are shown in red and blue.
#     The image above shows the 5 different ways {(0,2),(2,3)}, {(0,1),(1,3)}, {(0,1),(2,3)}, {(1,2),(2,3)}, {(0,1),(1,2)}.
#
# **Example 2:**
#
#
#
#     Input: n = 3, k = 1
#     Output: 3
#     Explanation: The 3 ways are {(0,1)}, {(0,2)}, {(1,2)}.
#
#
# **Example 3:**
#
#
#
#     Input: n = 30, k = 7
#     Output: 796297179
#     Explanation: The total number of possible ways to draw 7 line segments is 3796297200. Taking this number modulo 109 + 7 gives us 796297179.
#
#
# **Example 4:**
#
#
#
#     Input: n = 5, k = 3
#     Output: 7
#
#
# **Example 5:**
#
#
#
#     Input: n = 3, k = 2
#     Output: 1
#
#
#
# **Constraints:**
#
#   * `2 <= n <= 1000`
#   * `1 <= k <= n-1`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end