text:
zero(::Type{DoubleFloat{T}}) where {T<:AbstractFloat} = DoubleFloat(zero(T), zero(T))
one(::Type{DoubleFloat{T}}) where {T<:AbstractFloat} = DoubleFloat(one(T), zero(T))

nan(::Type{DoubleFloat{T}}) where {T<:AbstractFloat} = DoubleFloat(T(NaN), T(NaN))
inf(::Type{DoubleFloat{T}}) where {T<:AbstractFloat} = DoubleFloat(T(Inf), T(NaN))
posinf(::Type{DoubleFloat{T}}) where {T<:AbstractFloat} = DoubleFloat(T(Inf), T(NaN))
neginf(::Type{DoubleFloat{T}}) where {T<:AbstractFloat} = DoubleFloat(T(-Inf), T(NaN))

typemax(::Type{DoubleFloat{T}}) where {T<:AbstractFloat} = DoubleFloat{T}(typemax(T))
typemin(::Type{DoubleFloat{T}}) where {T<:AbstractFloat} = DoubleFloat{T}(typemin(T))

floatmax(::Type{DoubleFloat{T}}) where {T<:AbstractFloat} = DoubleFloat{T}(floatmax(T))
floatmin(::Type{DoubleFloat{T}}) where {T<:AbstractFloat} = DoubleFloat{T}(floatmin(T))

maxintfloat(::Type{DoubleFloat{Float64}}) = DoubleFloat{Float64}(ldexp(0.5, 107), ldexp(0.5, 54))
maxintfloat(::Type{DoubleFloat{Float32}}) = DoubleFloat{Float32}(ldexp(0.5f0, 49), ldexp(0.5f0, 25))
maxintfloat(::Type{DoubleFloat{Float16}}) = DoubleFloat{Float16}(65504.0) # floatmax(Double16)

isnan(x::DoubleFloat{T}) where {T<:IEEEFloat} = isnan(HI(x))
isinf(x::DoubleFloat{T}) where {T<:IEEEFloat} = (HI(x) === T(Inf))
isposinf(x::DoubleFloat{T}) where {T<:IEEEFloat} = (HI(x) === T(Inf))
isneginf(x::DoubleFloat{T}) where {T<:IEEEFloat} = (HI(x) === T(-Inf))

zero(::Type{DoubleFloat{DoubleFloat{T}}}) where {T<:AbstractFloat} = DoubleFloat(DoubleFloat(zero(T)), DoubleFloat(zero(T)))
one(::Type{DoubleFloat{DoubleFloat{T}}}) where {T<:AbstractFloat} = DoubleFloat(DoubleFloat(one(T)), DoubleFloat(zero(T)))

nan(::Type{DoubleFloat{DoubleFloat{T}}}) where {T<:AbstractFloat} = DoubleFloat(nan(DoubleFloat{T}), nan(DoubleFloat{T}))
inf(::Type{DoubleFloat{DoubleFloat{T}}}) where {T<:AbstractFloat} = DoubleFloat(inf(DoubleFloat{T}), zero(DoubleFloat{T}))
posinf(::Type{DoubleFloat{DoubleFloat{T}}}) where {T<:AbstractFloat} = DoubleFloat(posinf(DoubleFloat{T}), zero(DoubleFloat{T}))
neginf(::Type{DoubleFloat{DoubleFloat{T}}}) where {T<:AbstractFloat} = DoubleFloat(neginf(DoubleFloat{T}), zero(DoubleFloat{T}))

isnan(x::DoubleFloat{DoubleFloat{T}}) where {T<:IEEEFloat} = (HI(HI(x)) === T(NaN))
isinf(x::DoubleFloat{DoubleFloat{T}}) where {T<:IEEEFloat} = (HI(HI(x)) === T(Inf))
isposinf(x::DoubleFloat{DoubleFloat{T}}) where {T<:IEEEFloat} = (HI(HI(x)) === T(Inf))
isneginf(x::DoubleFloat{DoubleFloat{T}}) where {T<:IEEEFloat} = (HI(HI(x)) === T(-Inf))


text:
using ApproxFun, SpectralMeasures, Plots; Plots.gr(legend=false,linewidth=2,xlims=(-2,2),ylims=(0,2))

## 1by1 perturbation of Toeplitz
# J1 = [α/2 .5            ]
#      [ .5  0 .5         ]
#      [    .5  0 .5      ]
#      [       .5  0 .5   ]
#      [          ....... ]
#
#  We have several points of interest
#   α = 0 : Chebyshev of the 2nd Kind (semicircle)
#   α = -1 : Chebyshev of the 3rd kind
#   α = 1 : Chebyshev of the 4th kind
#  |α|≦ 1 : purely continuous spectrum
#  |α|> 1 : A single eigenvalue
#
#  In LightTable, "drag" k to see how the perturbation affects the spectral measure

k = 3
  plot(spectralmeasure([k/20],[.5]),title="\\alpha = $(k/10)")

## Basic 2by2 perturbation of Toeplitz
# J2 = [  0  β/⎷2            ]
#      [β/⎷2  0  .5          ]
#      [     .5   0  .5      ]
#      [         .5   0 .5   ]
#      [            ........ ]
#
#  We have several points of interest
#   β = 1/⎷2 : Chebyshev of the 2nd Kind
#   β = 1 : Chebyshev of the 1st kind
#  |β|≦ 1 : purely continuous spectrum
#  |β|> 1 : 2 eigenvalues
#
#  In LightTable, "drag" k to see how the perturbation affects the spectral measure

k = 5
  plot(spectralmeasure([0.],[.5+k/20]/sqrt(2)),title="\\beta=$(.5+k/20)")

####################
## Use the following for a closer look at what is going on with these operators

Δ = freejacobioperator()

a = 1.0
J1 = SymTriPertToeplitz([a],[.5],0.0,0.5)
C1 = connectioncoeffsoperator(J1)
Δ*C1-C1*J1
10.0*I
Q,L = ql(J1-10.0*I)

d,U = eigen(J1)

b = .5 + 0.1
J2 = SymTriPertToeplitz([0.],[b],0.0,0.5)

C2 = connectioncoeffsoperator(J2)
Δ*C2-C2*J2

ext: jl


text:
using MakieLayout
using Makie
using KernelDensity

function kdepoly!(la::LAxis, vec, reverse=false; kwargs...)
    kderesult = kde(vec; npoints=32)

    x = kderesult.x
    y = kderesult.density .* 100 # otherwise numbers are super small

    if reverse
        poly!(la, Point2.(y, x); kwargs...)
    else
        poly!(la, Point2.(x, y); kwargs...)
    end
end

begin
    scene = Scene(resolution = (1000, 1000));
    screen = display(scene)
    campixel!(scene);

    la1 = LAxis(scene)
    la2 = LAxis(scene)
    la3 = LAxis(scene)

    linkxaxes!(la1, la2)
    linkyaxes!(la1, la3)

    maingl = GridLayout(scene, 1, 1, alignmode=Outside(40), colsizes=Relative(1))

    gl2 = maingl[1, 1] = GridLayout(
        2, 2,
        rowsizes = [Relative(0.2), Relative(0.8)],
        colsizes = [Auto(), Aspect(1, 1)],
        addedrowgaps = [Fixed(10)],
        addedcolgaps = [Fixed(10)])

    gl2[2, 1] = la1
    la1.attributes.titlevisible[] = false

    gl2[1, 1] = la2
    la2.attributes.xlabelvisible[] = false
    la2.attributes.xticklabelsvisible[] = false
    la2.attributes.xticksvisible[] = false
    la2.attributes.titlevisible[] = false
    la2.attributes.ypanlock[] = true
    la2.attributes.yzoomlock[] = true

    gl2[2, 2] = la3
    la3.attributes.ylabelvisible[] = false
    la3.attributes.yticklabelsvisible[] = false
    la3.attributes.yticksvisible[] = false
    la3.attributes.titlevisible[] = false
    la3.attributes.xpanlock[] = true
    la3.attributes.xzoomlock[] = true

    maingl[0, 1] = LText(scene, text="Auto Limits", textsize=50)

    sleep(3)
    linkeddata = randn(200, 2) .* 15 .+ 50
    green = RGBAf0(0.05, 0.8, 0.3, 0.6)
    scatter!(la1, linkeddata, markersize=3, color=green, show_axis=false)
    sleep(1)
    kdepoly!(la2, linkeddata[:, 1], false, color=green, linewidth=2, show_axis=false)
    sleep(1)
    kdepoly!(la3, linkeddata[:, 2], true, color=green, linewidth=2, show_axis=false)
    sleep(1)

    linkeddata2 = randn(200, 2) .* 20 .+ 70
    red = RGBAf0(0.9, 0.1, 0.05, 0.6)
    scatter!(la1, linkeddata2, markersize=3, color=red, show_axis=false)
    sleep(1)
    kdepoly!(la2, linkeddata2[:, 1], false, color=red, linewidth=2, show_axis=false)
    sleep(1)
    kdepoly!(la3, linkeddata2[:, 2], true, color=red, linewidth=2, show_axis=false)
    sleep(1)

    linkeddata3 = randn(200, 2) .* 25 .+ 100
    blue = RGBAf0(0.05, 0.1, 0.9, 0.6)
    scatter!(la1, linkeddata3, markersize=3, color=blue, show_axis=false)
    sleep(1)
    kdepoly!(la2, linkeddata3[:, 1], false, color=blue, linewidth=2, show_axis=false)
    sleep(1)
    kdepoly!(la3, linkeddata3[:, 2], true, color=blue, linewidth=2, show_axis=false)
end


text:
import Pkg; Pkg.activate("../")
using Turing: Bijectors

b = Bijectors.SimplexBijector()
ib = inv(b)
x = [.3, .5, .2]
y = b(x)

Bijectors.logabsdetjac(b, x)
Bijectors.logabsdetjac(ib, y)

max_stars_repo_name: luiarthur/CytofDensityEstimation

max_issues_repo_name: luiarthur/CytofDensityEstimation

max_forks_repo_name: luiarthur/CytofDensityEstimation


text:
# Formulas from:
# http://www.fao.org/docrep/x0490e/x0490e07.htm

module Evapotranspiration

export Location, Angle, Degree, Radian,
       Kelvin, Celcius, Fahrenheit,
       RH, kPa, mm_day, MJ_m2hr, MJ_m2day, Meter,
       ave

abstract Angle
immutable Degree <: Angle
  value::Float64
end
immutable Radian <: Angle
  value::Float64
end
value(a::Angle) = a.value
Base.convert(::Type{Radian}, x::Degree) = π/180*value(x)
Base.promote_rule(::Type{Radian}, ::Type{Degree}) = Radian
Base.full{T<:Angle}(x::T) = x
Base.isnan{T<:Angle}(x::T) = isnan(x.value)
Base.isless{T<:Angle}(x::T, y::T) = isless(x,y)
Base.length{T<:Angle}(x::T) = 1
Base.zero{T<:Angle}(::Type{T}) = T(0.0)
Base.eltype{T<:Angle}(x::T) = eltype(value(x))
Base.getindex{T<:Angle}(x::T, i) = getindex(value(x), i)
Base.abs{T<:Angle}(x::T) = abs(value(x))
Base.sin(x::Radian) = sin(value(x))
Base.sin(x::Degree) = sind(value(x))
Base.cos(x::Radian) = cos(value(x))
Base.cos(x::Degree) = cosd(value(x))
Base.tan(x::Radian) = tan(value(x))
Base.tan(x::Degree) = tand(value(x))
for op in (:+, :-)
  @eval ($op){T<:Angle}(x::T, y::T) = T(($op)(x.value, y.value))
end
for op in (:*,)
  @eval ($op){T<:Radian,N<:Number}(x::N, y::T) = ($op)(x, y.value)
  @eval ($op){T<:Radian,N<:Number}(x::T, y::N) = ($op)(x.value, y)
end

abstract Temperature
immutable Celcius <: Temperature
  value::Float64
end
immutable Kelvin <: Temperature
  value::Float64
end
immutable Fahrenheit <: Temperature
  value::Float64
end
value(a::Temperature) = a.value
Base.convert(::Type{Kelvin},  x::Celcius)    = Kelvin(x.value + 273.15)
#Base.convert(::Type{Celcius}, x::Fahrenheit) = Celcius((x.value-32)/1.8)
Base.convert(::Type{Kelvin},  x::Fahrenheit) = Kelvin(1.8x.value - 457.67)
Kelvin(T::Celcius) = convert(Kelvin, T)
Kelvin(T::Fahrenheit) = convert(Kelvin, T)
Base.promote_rule(::Type{Kelvin}, ::Type{Celcius}) = Kelvin
Base.promote_rule(::Type{Kelvin}, ::Type{Fahrenheit}) = Kelvin
Base.promote_rule(::Type{Kelvin}, ::Type{Float64}) = Float64
Base.full{T<:Temperature}(x::T) = x
Base.isnan{T<:Temperature}(x::T) = isnan(x.value)
Base.isless{T<:Temperature}(x::T, y::T) = isless(x,y)
Base.length{T<:Temperature}(x::T) = 1
Base.zero{T<:Temperature}(::Type{T}) = T(0.0)
Base.eltype{T<:Temperature}(x::T) = eltype(value(x))
Base.getindex{T<:Temperature}(x::T, i) = getindex(value(x), i)
for op in (:+, :-)
  @eval ($op){T<:Temperature}(x::T, y::T) = T(($op)(x.value, y.value))
end
for op in (:*,)
  @eval ($op){T<:Kelvin}(x::T, y::T) = ($op)(x.value, y.value)
end
for op in (:+, :-, :*, :/)
  @eval ($op){T<:Temperature}(x::Number, y::T) = ($op)(x, y.value)
  @eval ($op){T<:Temperature}(x::T, y::Number) = ($op)(x.value, y)
end

immutable RH
  value::Float64
end
value(a::RH) = a.value
for op in (:+, :-)
  @eval ($op){T<:RH}(x::T, y::T) = T(($op)(x.value, y.value))
end
for op in (:+, :-, :*, :/)
  @eval ($op){T<:RH}(x::Number, y::T) = ($op)(x, y.value)
  @eval ($op){T<:RH}(x::T, y::Number) = ($op)(x.value, y)
end

immutable MJ_m2day  # MJ/(m^2*day)  (Radiation)
  value::Float64
end
value(a::MJ_m2day) = a.value
for op in (:+, :-)
  @eval ($op){T<:MJ_m2day}(x::T, y::T) = T(($op)(x.value, y.value))
end
for op in (:/,)
  @eval ($op){T<:MJ_m2day}(x::T, y::T) = ($op)(x.value, y.value)
end
for op in (:+, :-, :*, :/)
  @eval ($op){T<:MJ_m2day}(x::Number, y::T) = ($op)(x, y.value)
  @eval ($op){T<:MJ_m2day}(x::T, y::Number) = ($op)(x.value, y)
end
Base.full{T<:MJ_m2day}(x::T) = x
Base.isnan{T<:MJ_m2day}(x::T) = isnan(x.value)
Base.isless{T<:MJ_m2day}(x::T, y::T) = isless(x,y)
Base.length{T<:MJ_m2day}(x::T) = 1
Base.zero{T<:MJ_m2day}(::Type{T}) = T(0.0)
Base.eltype{T<:MJ_m2day}(x::T) = eltype(value(x))
Base.getindex{T<:MJ_m2day}(x::T, i) = getindex(value(x), i)

immutable MJ_m2hr  # MJ/(m^2*hr)  (Radiation)
  value::Float64
end
value(a::MJ_m2hr) = a.value
for op in (:+, :-)
  @eval ($op){T<:MJ_m2hr}(x::T, y::T) = T(($op)(x.value, y.value))
end
for op in (:/,)
  @eval ($op){T<:MJ_m2hr}(x::T, y::T) = ($op)(x.value, y.value)
end
for op in (:+, :-, :*, :/)
  @eval ($op){T<:MJ_m2hr}(x::Number, y::T) = ($op)(x, y.value)
  @eval ($op){T<:MJ_m2hr}(x::T, y::Number) = ($op)(x.value, y)
end
Base.full{T<:MJ_m2hr}(x::T) = x
Base.isnan{T<:MJ_m2hr}(x::T) = isnan(x.value)
Base.isless{T<:MJ_m2hr}(x::T, y::T) = isless(x,y)
Base.length{T<:MJ_m2hr}(x::T) = 1
Base.zero{T<:MJ_m2hr}(::Type{T}) = T(0.0)
Base.eltype{T<:MJ_m2hr}(x::T) = eltype(value(x))
Base.getindex{T<:MJ_m2hr}(x::T, i) = getindex(value(x), i)

immutable mm_day  # mm/day (Evapotranspiration)
  value::Float64
end
value(a::mm_day) = a.value
for op in (:+, :-)
  @eval ($op){T<:mm_day}(x::T, y::T) = T(($op)(x.value, y.value))
end
for op in (:+, :-, :*, :/)
  @eval ($op){T<:mm_day}(x::Number, y::T) = ($op)(x, y.value)
  @eval ($op){T<:mm_day}(x::T, y::Number) = ($op)(x.value, y)
end
Base.full{T<:mm_day}(x::T) = x
Base.isnan{T<:mm_day}(x::T) = isnan(x.value)
Base.isless{T<:mm_day}(x::T, y::T) = isless(x,y)
Base.length{T<:mm_day}(x::T) = 1
Base.zero{T<:mm_day}(::Type{T}) = T(0.0)
Base.eltype{T<:mm_day}(x::T) = eltype(value(x))
Base.getindex{T<:mm_day}(x::T, i) = getindex(value(x), i)

immutable kPa  # kPa (1000 Pascals)
  value::Float64
end
value(a::kPa) = a.value
for op in (:+, :-)
  @eval ($op){T<:kPa}(x::T, y::T) = T(($op)(x.value, y.value))
end
for op in (:+, :-, :*, :/)
  @eval ($op){T<:kPa}(x::Number, y::T) = ($op)(x, y.value)
  @eval ($op){T<:kPa}(x::T, y::Number) = ($op)(x.value, y)
end
Base.full{T<:kPa}(x::T) = x
Base.isnan{T<:kPa}(x::T) = isnan(x.value)
Base.isless{T<:kPa}(x::T, y::T) = isless(x,y)
Base.length{T<:kPa}(x::T) = 1
Base.zero{T<:kPa}(::Type{T}) = T(0.0)
Base.eltype{T<:kPa}(x::T) = eltype(value(x))
Base.getindex{T<:kPa}(x::T, i) = getindex(value(x), i)

immutable Meter
  value::Float64
end
value(a::Meter) = a.value
for op in (:+, :-)
  @eval ($op){T<:Meter}(x::T, y::T) = T(($op)(x.value, y.value))
end
for op in (:+, :-, :*, :/)
  @eval ($op){T<:Meter}(x::Number, y::T) = ($op)(x, y.value)
  @eval ($op){T<:Meter}(x::T, y::Number) = ($op)(x.value, y)
end
Base.full{T<:Meter}(x::T) = x
Base.isnan{T<:Meter}(x::T) = isnan(x.value)
Base.isless{T<:Meter}(x::T, y::T) = isless(x,y)
Base.length{T<:Meter}(x::T) = 1
Base.zero{T<:Meter}(::Type{T}) = T(0.0)
Base.eltype{T<:Meter}(x::T) = eltype(value(x))
Base.getindex{T<:Meter}(x::T, i) = getindex(value(x), i)

for op in (:*,)
  @eval ($op)(x::kPa, y::RH) = kPa(($op)(x.value, y.value))
  @eval ($op)(x::RH, y::kPa) = kPa(($op)(x.value, y.value))
end

immutable Location
  latitude::Radian  # above (+ve) or below (-ve) the equator
  Lz::Degree # Lz = longitude of the centre of the local time zone [degrees west of Greenwich]
  Lm::Degree # Lm = longitude of the measurement site [degrees west of Greenwich]
  z::Meter   # elevation above sea level
  # Angstrom's values:
  as::Float64 # regression constant
  bs::Float64
  albedo::Float64
  function Location(lat, Lz, Lm, z=Meter(100), as=0.25, bs=0.5, albedo=0.23)
     new(lat, Lz, Lm, z, as, bs, albedo)
  end
end
function Location(latitude::Degree, Lz::Degree, Lm::Degree, z::Meter=Meter(100), as=0.25, bs=0.5, albedo=0.23)
  Location(radians(latitude), Lz, Lm, z, as, bs, albedo)
end
latitude(loc::Location) = loc.latitude
Lz(loc::Location) = loc.Lz
Lm(loc::Location) = loc.Lm
Lzm(loc::Location) = value(Lz(loc) - Lm(loc))
elevation(loc::Location) = loc.z
as(loc::Location) = loc.as
bs(loc::Location) = loc.bs
albedo(loc::Location) = loc.albedo

using Dates

if false
if !haskey(ENV, "FORECAST_API_KEY")
  error("tests expect a FORECAST_API_KEY to be set in your environment variables")
end

options = {"APIKey" => ENV["FORECAST_API_KEY"],
           "units" => "ca"}
end

ave(x::Number, y::Number) = (x+y)/2
ave{T}(x::T, y::T) = T((x+y)/2)
# Midpoint of hours
function ave(::Type{Hour}, hr1, hr2)
    if hr2 < hr1  # time wrapped around
        hr2 += 24
    end
    (hr1 + hr2)/2
end

# Atmospheric pressure (P)
# P =  atmospheric pressure [kPa],
# z = elevation above sea level [m],
# T = temperature [°C]
P(z::Meter, T::Kelvin) = kPa(101.3((T-0.0065z)/T)^5.26)   # (7)
P(z::Meter, T::Temperature=Celcius(20)) = P(z, Kelvin(T))

## psychrometric constant, gamma, kPa/C
# γ  = psychrometric constant [kPa/°C],
# p  = atmospheric pressure [kPa],
# λ  = latent heat of vaporization, 2.45 [MJ/kg],
# cp = specific heat at constant pressure, 1.013e-3 [MJ/kg/°C],
# e  = ratio molecular weight of water vapour/dry air = 0.622.
γ(z::Meter, T::Kelvin, e=0.622, cp=1.013e-3, λ=2.54) = cp*P(z,T)/(e*λ)  # (8)
γ(p::kPa, e=0.622, cp=1.013e-3, λ=2.54) = cp*p/(e*λ)  # (8)
γ(z::Meter, T::Temperature=Celcius(20), args...) = γ(z, convert(Kelvin, T), args...)

## Relativive humidity, Rh
# ea = actual vapour pressure [kPa],
# e° = saturation vapour pressure at temperature T [kPa],
RH(ea::kPa,T=20) = 100ea/e°(T)   # (10)

## Saturation vapour pressure, e°, in kPa:
# e° = saturation vapour pressure at temperature T [kPa],
# T  = air temperature [°C],
e°(T::Celcius) = kPa(0.6108*exp(17.27T/(T+237.3)))  # (11)

## Mean saturation vapour pressure, Es, in kPa:
# Julia 0.3 bug: would like to use eₛ (instead of Es) but it doesn't work.
Es(Tmin::Celcius, Tmax::Celcius) = ave(e°(Tmax), e°(Tmin))   # (12)

# Slope of saturation vapour pressure curve (Δ) [kPa/C]
Δ(T::Celcius) = 4098*(0.6108*exp(17.27T/(T+237.3)))/(T+237.3)^2 # (13)
Δ(Tmin::Celcius, Tmax::Celcius) = Δ(ave(Tmin, Tmax))

# Actual vapour pressure (Ea) derived from dewpoint temperature
Ea(Tdew::Celcius) = e°(Tdew) # (14)
Ea(Tdew::Temperature) = eᵒ(convert(Celcius, Tdew))

# daily mean RH comes from forecast.io
# This is preferred over Ea(Tdew)
Ea(Tmin::Celcius, Tmax::Celcius, RHmean::RH) = RHmean*ave(e°(Tmax),e°(Tmin)) # (19)
Ea(Tmin::Temperature, Tmax::Temperature, RHmean::RH) = Ea(Celcius(Tmin), Celcius(Tmax), RHmean)

# This is preferred over other methods
Ea(Tmin::Celcius, Tmax::Celcius, RHmin::RH, RHmax::RH) = ave(RHmin*e°(Tmax), RHmax*e°(Tmin)) # (17)

## Vapor pressure deficit, Es - Ea, kPa:
Esa(Tmin::Celcius, Tmax::Celcius; Tdew::Celcius=error("No Tdew")) = Es(Tmin, Tmax) - Ea(Tdew)
Esa(Tmin::Celcius, Tmax::Celcius, RHmean::RH) = Es(Tmin, Tmax) - Ea(Tmin, Tmax, RHmean)
Esa(Tmin::Celcius, Tmax::Celcius, RHmin::RH, RHmax::RH) = Es(Tmin, Tmax) - Ea(Tmin, Tmax, RHmin, RHmax)

######################
# Radiation:
######################

# Extraterrestial radiation

MJ2mm(Ra::MJ_m2day) = MJ_m2day(value(Ra)/2.45)  # (20)  # converts MJ to mm (over whatever time period (eg Day or Hour))
MJ2mm(Ra::MJ_m2hr) = MJ_m2day(value(Ra)/2.45)

# Maximum extraterrestrial radiation for daily periods, Ra
# Ra  = extraterrestrial radiation [MJ/m^2/day],
# ϕ   = latitude [rad] (Equation 22),
# dᵣ  = inverse relative distance Earth-Sun (Equation 23),
# δ   = solar decimation (Equation 24) [rad].
# ωs  = sunset hour angle (Equation 25) [rad],
const Gsc = 0.0820 # Solar constant, Gsc, in MJ/(m^2*min):
function _Ra(ϕ::Radian, dᵣ, δ::Radian, ωs::Radian)
  MJ_m2day(24*60/π * Gsc * dᵣ * (ωs*sin(ϕ)*sin(δ) + cos(ϕ)*cos(δ)*sin(ωs)))  # (21)
end
# latitude = in [deg]
function Ra(::Type{Day}, J, loc::Location)  # (21)
  ϕ = latitude(loc)
  d = δ(J)
  _Ra(ϕ, dᵣ(J), d, ωs(ϕ, d))
end

radians(deg::Degree) = Radian(π/180*value(deg))   # (22)
radians(loc::Location) = radians(latitude(loc))
# dᵣ inverse relative distance Earth-Sun, dᵣ, unitless
dᵣ(J::Day) = 1 + 0.033cos(2π/365*Dates.value(J))  # (23)
# solar decimation, δ
δ(J::Day) = Radian(0.409*sin(2π/365*Dates.value(J)-1.39))  # (24)
# sunset hour angle, ws, in radians (zero is solar noon)
ωs(ϕ::Radian,δ::Radian) = Radian(acos(-tan(ϕ)*tan(δ)))   # (25)
# sun angle at midpoint of period shorter than a day, ws, in radians (zero is solar noon)
# t = standard clock time at the midpoint of the period [hour].
#   - For example for a period between 14.00 and 15.00 hours, t = 14.5
# Lz = longitude of the centre of the local time zone [degrees west of Greenwich]
# Lm = longitude of the measurement site [degrees west of Greenwich]
ωs(::Type{Hour}, t, J::Day, loc::Location) = π/12*((t+0.06667*Lzm(loc) + Sc(J))-12)  # (31)
# Sc seasonal correction for solar time [hour]:
function Sc(J::Day)
    b = 2π*(Dates.value(J)-81)/364   # (33)
    0.1645*sin(2b)-0.1255*cos(b)-0.025*sin(b)  # (32)
end

# Max daylight hours
daylighthours(J::Day, loc::Location) = 24/π*ωs(latitude(loc), δ(J))  # (34)
# n = average hours of sun per day
cloudcover(n, J::Day, loc::Location) = 1-n/daylighthours(J,loc)

function duration(::Type{Hour}, hr1, hr2)
    if hr2 < hr1  # time wrapped around
        hr2 += 24
    end
    hr2 - hr1
end
# hr is clock hrs (at midpoint of time period)
function issunset(T::Type{Hour}, hr, J::Day, loc::Location)
    wsun = ωs(T, hr, J, loc)  # 0 rad is noon
    ϕ = latitude(loc)
    d = δ(J)
    # sunset hour angle
    wsunset = ωs(ϕ, d)
    abs(wsun) > abs(wsunset)
end
# Extraterrestrial radiation for hourly periods, Ra, in MJ/(m^2*hour)
function Ra(T::Type{Hour}, J::Day, hr1, hr2, loc::Location) # (28)
    # angle above (below) equator
    ϕ = latitude(loc)
    # solar decimation, d, radians
    d = δ(J)
    # midpiont of hours, eg if 14h and 15h then t=14.5
    hr = ave(T, hr1, hr2)
    # sunset hour angle, in radians
    if issunset(T, hr, J, loc)
        # sun is below horizon:
        return ra = MJ_m2hr(0)
    end
    # The solar time angles at the beginning and end of the period are given by:
    t1 = duration(Hour, hr1, hr2)
    wsun = ωs(T, hr, J, loc)  # sun angle, 0 rad is noon
    w1 = wsun - π*t1/24   # (29)
    w2 = wsun + π*t1/24   # (30)
    ra = 12*60/π*Gsc*dᵣ(J)*((w2-w1)*sin(ϕ)*sin(d)+cos(ϕ)*cos(d)*(sin(w2)-sin(w1)))  # (28)
    MJ_m2hr(ra)
end

### Solar or Shortwave Radiation, Rs

# The amount of radiation that penetrates the atmosphere, Rs, units same as Ra
# Rs solar or shortwave radiation
# cloud_cover = percentage of cloud coverage, 0.5 is half covered
# Ra = extraterrestrial radiation (same units as Rs)
# as = regression constant, expressing the fraction of extraterrestrial radiation reaching the earth on overcast days (n = 0),
# as+bs = fraction of extraterrestrial radiation reaching the earth on clear days (n ≡ N (clear sky)).
# n = hours of sun per day
# N = max hours of sun per day
function Rs{T}(cloud_cover, J::Day, Ra::T, loc::Location)  # (35)
    N = daylighthours(J, loc)
    sun_pct = (1 - cloud_cover)
    coef = as(loc) + bs(loc)*sun_pct
    rs = T(coef*Ra)
end

# Clear-sky (n ≡ N) solar radiation, Rso
# When as and bs are available:
Rso(Ra,as,bs) = (as+bs)*Ra  # (36)
# When as and bs are not available:
# z = station elevation above sea level [m]
Rso(Ra,z) = (0.75 + 2e-5*z)*Ra  # (37)

# Rs/Rso: Ratio of solar radiation that reaches the earth (compared to a clear sky day)
function Rsso(cloud_cover, J::Day, Ra, loc)
    r1 = Rs(cloud_cover, J, Ra, loc)
    r2 = Rs(0,           J, Ra, loc)
    r1/r2
end

# Net solar radiation
# a = albedo or canopy reflection coefficient
# For example a=0.23 for the hypothetical grass reference crop [dimensionless],
Rns(Rs, loc::Location) = MJ_m2day((1-albedo(loc))Rs)  # (38)

# Net long wave radiation, Rnl [MJ/(m^*day)]
# Rnl = net outgoing longwave radiation [MJ m-2 day-1],
# σ = Stefan-Boltzmann constant [4.903 10-9 MJ K-4 m-2 day-1],
# Tmax = maximum temperature during the 24-hour period [°C],
# Tmin = minimum temperature during the 24-hour period [°C],
# Ea   = actual vapour pressure [kPa],
# Rs/Rso = relative shortwave radiation, 0 <= Rsso <= 1
# Rs measured or calculated. (Equation 35) solar radiation [MJ m-2 day-1],
# Rso calculated (Equation 36 or 37) clear-sky radiation [MJ m-2 day-1].
function Rnl(Rsso, Tmin::Kelvin, Tmax::Kelvin, Ea::kPa)  # (39)
    σ = 4.903e-9
    blackbody = σ * ave(Tmax^4, Tmin^4)
    air_humidity_correction = 0.34 - 0.14*sqrt(value(Ea))
    effect_of_cloudiness = 1.35 * Rsso - 0.351
    MJ_m2day(blackbody * air_humidity_correction * effect_of_cloudiness)
end
function Rnl(Rsso, Tmin::Temperature, Tmax::Temperature, Ea::kPa)
    Rnl(Rsso, convert(Kelvin, Tmin), convert(Kelvin, Tmax), Ea::kPa)
end

# Daily radiation, Rn [MJ/(m^*day)]:
Rn(Rns, Rnl) = Rns - Rnl   # (40)

function net_radiation(T::Type{Day}, cloud_cover, J::Day, Tmin::Temperature, Tmax::Temperature, RHmin::RH, RHmax::RH, loc::Location)
  ra   = Ra(T, J, loc)
  rs   = Rs(cloud_cover, J, ra, loc)
  rsso = Rsso(cloud_cover, J, ra, loc)
  ea   = Ea(Tmin, Tmax, RHmin, RHmax)
  rnl  = Rnl(rsso, Tmin, Tmax, ea)
  rns  = Rns(rs, loc)
  rn   = Rn(rns, rnl)
end
function net_radiation(T::Type{Day}, cloud_cover, J::Day, Tmin::Temperature, Tmax::Temperature, RHmean::RH, loc::Location)
  ra   = Ra(T, J, loc)
  rs   = Rs(cloud_cover, J, ra, loc)
  rsso = Rsso(cloud_cover, J, ra, loc)
  ea   = Ea(Tmin, Tmax, RHmean)
  rnl  = Rnl(rsso, Tmin, Tmax, ea)
  rns  = Rns(rs, loc)
  rn   = Rn(rns, rnl)
end

# Net radiation [MJ/(m^2*hour)]
function net_radiation(T::Type{Hour}, cloud_cover, J::Day, Tmin, Tmax, RHmin::RH, RHmax::RH, hr1, hr2, loc::Location)
  ra   = Ra(T, J, hr1, hr2, loc) # (28) [MJ/(m^2*hour)]
  rs   = Rs(cloud_cover, J, ra, loc)
  rsso = Rsso(cloud_cover, J, ra, loc)
  ea   = Ea(Tmin, Tmax, RHmin, RHmax)
  rnl  = Rnl(rsso, Tmin, Tmax, ea)
  rns  = Rns(rs, loc)
  rn   = Rn(rns, rnl)
end
function net_radiation(T::Type{Hour}, cloud_cover, J::Day, Tmin, Tmax, RHmean::RH, hr1, hr2, loc::Location)
  ra   = Ra(T, J, hr1, hr2, loc) # (28) [MJ/(m^2*hour)]
  rs   = Rs(cloud_cover, J, ra, loc)
  rsso = Rsso(cloud_cover, J, ra, loc)
  ea   = Ea(Tmin, Tmax, RHmean)
  rnl  = Rnl(rsso, Tmin, Tmax, ea)
  rns  = Rns(rs, loc)
  rn   = Rn(rns, rnl)
end

# Soil heat flux (inaccurate method for over many days)
# G  = soil heat flux [MJ/(m^2*day)],
# cs = soil heat capacity [MJ/(m^3*°C)],
# T1 = air temperature at time 1 [°C],
# T2 = air temperature at time 2 [°C],
# Δt = length of time interval between T1 and T2 [day],
# Δz = effective soil depth [m].
# For a day Δz is between 0.1 and 0.2, 1 month is about 1, 4 months can be 2.0.
x = [1, 30, 4*30]
y = [0.1, 1, 2]
const Δz_slope = (x\y)[1]  # line fit
const Δz_b = y[1]
Δz(::Type{Day}, days) = Δz_slope * days + Δz_b
function soil_heat_flux{Temp<:Union(Celcius,Kelvin)}(T::Type{Day}, Δt, T1::Temp, T2::Temp, Δz=Δz(T,Δt), soil_heat_capacity=2.1)
    0.408*soil_heat_capacity * (T2-T1) * Δz / Δt   # (41)
end

# For monthly periods
# T1 = mean air temperature of previous month
# T2 = mean air temperature of current month
# T3 = mean air temperature of next month
# Δt = length of time interval between T1 and T2 [month],
function soil_heat_flux(T::Type{Month}, Δt, T1::Celcius, T2::Celcius)  # (43) and (44)
    warming = sign(value(T2-T1))
    G = warming*0.14/Δt*(T2-T1)
end

# For hourly or shorter periods
function soil_heat_flux(T::Type{Hour}, hr1, hr2, Rn, latitude, J, Lz, Lm, T1, T2=T1+1)  # (45)
    warming = sign(T2-T1)
    hr = ave(T, hr1, hr2)
    coef = ifelse(issunset(T, hr, latitude, J, Lz, Lm), 0.5, 0.1)
    G = warming*coef*Rn
end

# Windspeed height compensation (for weather stations)
# u2 = wind speed at 2 m above ground surface [m/s],
# uz = measured wind speed at z m above ground surface [m/s],
# z  = height of measurement above ground surface [m].
function u2(windspeed, z=10)
    # speeds are in m/s or kmph
    return windspeed*4.87/log(67.8z-5.42)
end

# ETo = reference evapotranspiration [mm/day],
# Δ   = slope vapour pressure curve [kPa/°C],
# Rn  = net radiation at the crop surface [MJ/(m^2*day)],
# G   = soil heat flux density [MJ/(m^2*day)],
# γ   = psychrometric constant [kPa/°C].
# T   = mean daily air temperature at 2 m height [K],
# u₂  = wind speed at 2 m height [m/s],
# Es  = saturation vapour pressure [kPa],
# Ea  = actual vapour pressure [kPa],
# Esa = Es - Ea = saturation vapour pressure deficit [kPa],
function ETo(Δ, Rn, G, γ, T::Kelvin, u₂, Esa)
   (0.408*Δ*(Rn-G) + 900γ/value(T)*u₂*Esa)/(Δ+γ*(1+0.34u₂))  # (6)
end
function ETo(Td::Type{Day}, Δt, Tmin::Celcius, Tmax::Celcius, RHmin::RH, RHmax::RH, P_kPa::kPa, u10_kph, cloud_cover, J::Day, loc::Location)
    D = Δ(Tmin,Tmax)
    rn = net_radiation(Td, cloud_cover, J, Tmin, Tmax, RHmin, RHmax, loc)
    G = soil_heat_flux(Td, Δt, Tmin, Tmax)
    g = γ(P_kPa)
    T = Kelvin(ave(Tmin,Tmax))
    u₂ = u2(u10_kph/3.6, 10)
    esa = Esa(Tmin, Tmax, RHmin, RHmax)
    ETo(D, rn, G, g, T, u₂, esa)
end
function ETo(Td::Type{Day}, Δt, Tmin::Celcius, Tmax::Celcius, RHmean::RH, P_kPa::kPa, u10_kph, cloud_cover, J::Day, loc::Location)
    D = Δ(Tmin,Tmax)
    rn = net_radiation(Td, cloud_cover, J, Tmin, Tmax, RHmean, loc)
    G = soil_heat_flux(Td, Δt, Tmin, Tmax)
    g = γ(P_kPa)
    T = Kelvin(ave(Tmin,Tmax))
    u₂ = u2(u10_kph/3.6, 10)
    esa = Esa(Tmin, Tmax, RHmean)
    ETo(D, rn, G, g, T, u₂, esa)
end

saturation_vapor_pressure(ave_temp_celcius) = 10^(8.07131-1730.63/(ave_temp_celcius+233.426))

function saturation_vapor_pressure_slope(ave_temp_celcius, dt=0.1)
    p1 = saturation_vapor_pressure(ave_temp_celcius)
    p2 = saturation_vapor_pressure(ave_temp_celcius+dt*5/9)
    return (p2-p1)/dt
end

temperature_coef(ave_temp_celcius, scaling=0.03567) = scaling*5.67e-8*(273+ave_temp_celcius)^4

function evaporation_mm_per_day(wind_speed_kmph, ave_temp_celcius, dew_point_celcius)
    mi_per_km = 1/1.61
    wind_speed_mi_per_day = wind_speed_kmph * 24 * mi_per_km
    ea = saturation_vapor_pressure(ave_temp_celcius)
    ed = saturation_vapor_pressure(dew_point_celcius)
    0.35*(ea-ed)*(1+0.0098*wind_speed_mi_per_day)
end

function degree_days_mm_per_day(ave_temp_celcius, dew_point_celcius, solar_radiation, sunshine_pct)
    mi_per_km = 1/1.61
    Ra = solar_radiation
    Re = 0.15 # reflexivity of forrest
    T = ave_temp_celcius
    Tcoef = 0.03567*5.67e-8*(273+T)^4
    ed = saturation_vapor_pressure(dew_point_celcius)
    S = sunshine_pct
    (Ra*(1-Re)*(0.18+0.55*S))-(Tcoef*(0.56-0.092*ed^0.5)*(0.1+0.9*S))
end

function potential_evaporation(ave_temp_celcius, dew_point_celcius, wind_speed_kmph, solar_radiation, sunshine_pct)
    Slope = saturation_vapor_pressure_slope(ave_temp_celcius)
    E = evaporation_mm_per_day(wind_speed_kmph, ave_temp_celcius, dew_point_celcius)
    H = degree_days_mm_per_day(ave_temp_celcius, dew_point_celcius, solar_radiation, sunshine_pct)
    ((Slope*H)+(0.27*E))/(Slope+0.27)
end

end # module

#function print_evaporation(forecast_data) {
#    vec = forecast_data
#    totalDewPoint = 0
#    totalRH = 0
#    numDewPoints = 0
#    for (i=0; i < vec.length; i++) {
#        v = vec[i]
#	v.date = Date(v.time*1000)
#	v.dom = v.date.getDate()
#	v.hr = v.date.getHours()
#        ave_temp_c = NaN
#        if ("temperature" in v) {
#            ave_temp_c = v.temperature
#        } else if (("temperatureMax" in v) && ("temperatureMin" in v)) {
#            ave_temp_c = v.temperatureMin + 0.75*(v.temperatureMax - v.temperatureMin)
#        }
#        dew_temp_c = v.dewPoint
#        solar_radiation = 6
#        sunshine = 1-v.cloudCover; # account for night having no sun.
#        wind_speed_kmph = v.windSpeed
#        E = evaporation_mm_per_day(wind_speed_kmph, ave_temp_c, dew_temp_c)
#        H = degree_days_mm_per_day(ave_temp_c, dew_temp_c, 6, sunshine)
#                U = potential_evaporation(ave_temp_c, dew_temp_c, wind_speed_kmph, solar_radiation, sunshine)
#        precipIntensity = NaN
#        if ("percipIntensityMax" in v) {
#            precipIntensity = v.percipIntensityMax
#        }
#        if ("precipIntensity" in v) {
#            precipIntensity = v.precipIntensity
#        }
#        percip_mm_day = v.precipIntensity; # * vec[i].precipProbability
#        if (v.dom == 18) {
#		console.log(v.date+': '+v.humidity+', '+v.dewPoint)
#		totalDewPoint += v.dewPoint
#		totalRH += v.humidity
#		numDewPoints += 1
#	}
#        if (typeof percip_mm_day == 'undefined') {
#            percip_mm_day = 0
#        }
#        AM = percip_mm_day - U
#        t = Date(1000*vec[i].time)
#        #console.log(vec[i].time+": "+t + ": E="+ E.toFixed(2) + ", H=" + H.toFixed(2)+", ETo="+U.toFixed(2)+", percip="+percip_mm_day.toFixed(2)+", AM="+AM.toFixed(2))
#
#	# New method:
#	Tmin = v.temperatureMin
#	Tmax = v.temperatureMax
#	Tdew = v.dewPoint
#	RHmean = v.humidity
#	P_kPa = v.pressure
#	# TODO: scale down wind because of shelter from trees, maybe 0.25
#	windspeed_10m_kmph = v.windSpeed
#	cloud_cover = v.cloudCover
#	J = J(v.time)
#        eto_day = ETo_day(Tmin, Tmax, RHmean, P_kPa, windspeed_10m_kmph, cloud_cover, latitude, J)
#	eto_diff = U - eto_day
#        #console.log("ETo day="+eto_day.toFixed(2)+" [mm/day]   diff="+eto_diff)
#    }
#    console.log("Total dew: "+ totalDewPoint+" RH: "+totalRH)
#    aveDewPoint = totalDewPoint/numDewPoints
#    aveRH = totalRH/numDewPoints
#    console.log("Ave: Dew Point: "+ aveDewPoint+", RH: "+ aveRH)
#}
#
#
#function print_daily_forecast(data) {
#  daily = data.daily.data
#  for (i=0; i<daily.length; i++) {
#      di = daily[i]
#      console.log(Date(1000*di.time) + ":")
#      console.log("  - sunrise: " + Date(1000*di.sunriseTime))
#      console.log("  - sunset:  " + Date(1000*di.sunsetTime))
#      console.log("  - TmaxTime: " + Date(1000*di.temperatureMaxTime))
#      console.log("  - TminTime:  " + Date(1000*di.temperatureMinTime))
#      console.log(di)
#  }
#}
#
#forecast.get(latitude, longitude, options, function (err, res, data) {
#  if (err) throw err
#  #console.log('res: ' + util.inspect(res))
#  #console.log('data: ' + util.inspect(data, true, 5))
#  #console.log('daily data: ' + util.inspect(data, true, 5))
#  #print_daily_forecast(data)
#  print_evaporation(data.daily.data)
#  print_evaporation(data.hourly.data)
#})
#
#


text:
#-------------------------------------------------------------------------------
# File: for-test.jl
# Description: This file contains bad complexity functions relatives to moves.
# Date: November 4, 2019
# Author: Jonathan Fontaine, Killian Fretaud, Rémi Garcia,
#         Boualem Lamraoui, Benoît Le Badezet, Benoit Loger
#-------------------------------------------------------------------------------

"""
    TEST_move_exchange!(solution::Solution, car_pos_a::Int, car_pos_b::Int, instance::Instance)

Bad complexity version for function move_exchange.
"""
function TEST_move_exchange!(solution::Solution, car_pos_a::Int, car_pos_b::Int, instance::Instance)
    if car_pos_a > car_pos_b
        return TEST_move_exchange!(solution, car_pos_b, car_pos_a, instance)
    end
    if car_pos_a == car_pos_b
        return solution
    end
    solution.sequence[car_pos_a], solution.sequence[car_pos_b] = solution.sequence[car_pos_b], solution.sequence[car_pos_a]
    update_matrices!(solution, instance)

    return solution
end

"""
    TEST_cost_move_exchange(solution::Solution, car_pos_a::Int, car_pos_b::Int,
                       instance::Instance, objective::Int)

Bad complexity version for function cost_move_exchange.
"""
function TEST_cost_move_exchange(solution::Solution, car_pos_a::Int, car_pos_b::Int,
                            instance::Instance, objective::Int)
    s = deepcopy(solution)
    s = TEST_move_exchange!(s, car_pos_a, car_pos_b, instance)

    return cost(s, instance, objective) - cost(solution, instance, objective)
end

"""
    TEST_move_insertion!(solution::Solution, old_index::Int, new_index::Int, instance::Instance)

Bad complexity version for function move_insertion.
"""
function TEST_move_insertion!(solution::Solution, old_index::Int, new_index::Int, instance::Instance)

    car_inserted = solution.sequence[old_index]
    if old_index < new_index
        for car_moved_pos in old_index:(new_index-1)
            solution.sequence[car_moved_pos] = solution.sequence[car_moved_pos+1]
        end
        solution.sequence[new_index] = car_inserted
    end
    if old_index > new_index
        for car_moved_pos in old_index:-1:(new_index+1)
            solution.sequence[car_moved_pos] = solution.sequence[car_moved_pos-1]
        end
        solution.sequence[new_index] = car_inserted
    end
    update_matrices!(solution, instance)

    return solution
end

"""
    TEST_cost_move_insertion(solution::Solution, car_pos_a::Int,
                        instance::Instance, objective::Int)

Bad complexity version for function cost_move_insertion.
"""
function TEST_cost_move_insertion(solution::Solution, car_pos_a::Int,
                             instance::Instance, objective::Int)

    cost_on_objective = zeros(Int, instance.nb_cars, 3)
    for i in (instance.nb_late_prec_day+1):instance.nb_cars
        s = deepcopy(solution)
        s = TEST_move_insertion!(s, car_pos_a, i, instance)
        c = cost(s, instance, objective) - cost(solution, instance, objective)
        cost_on_objective[i, :] = c
    end

    return cost_on_objective
end


text:
"""
Now if we pad to the start on the left hand side
and read to the right we see that the bit strings are all uniquely
identifiable. For example if I see 101, I know that represents
#9 and isn't just a subset for something else.

 0, 010
 4, 1100
 5, 1101
 6, 011
 7, 00
 8, 100
 9, 101
10, 1110
11, 11110
16, 111110
17, 1111110
18, 1111111

Note however that bits are never processed left to right, they are
processed right to left. So for processing of our data we need to reverse
the bits. Also note that this does not mean swapping the first and last bit
in a byte since the bit strings are of a given length.

So here is a side by side of the original values and what we need.

    Orig         New   Level & Index    New as # and hex
 0,     010 -     010   L3 - 1          2
 4,    1100 -    0011       L4 - 1      3
 5,    1101 -    1011       L4 - 2      11 b
 6,     011 -     110   L3 - 2          6
 7,      00 -      00 L2 - 1            0
 8,     100 -     001   L3 - 3          1
 9,     101 -     101   L3 - 4          5
10,    1110 -    0111       L4 - 3      7
11,   11110 -   01111 L5 - 1
16,  111110 -  011111 L6 - 1
17, 1111110 - 0111111 L7 - 1
18, 1111111 - 1111111 L7 - 2

"""


text:
using Plots
using DataFrames
using LinearAlgebra
using Distributions

"""
    vowelPlot(f1, f2, cats; [meansOnly=false, addLabels=true, ell=false, ellPercent=0.67, nEllPts=500, markersize=1, linewidth=2, kw...])

Create an F1-by-F2 vowel plot. The `f1` values are displayed along the x-axis, and the `f2` values are displayed along the y-axis, with each unique vowel class in `cats` being represented with a new color. The series labels in the legend will take on the unique values contained in `cats`. The alternate display whereby reversed F2 is on the x-axis and reversed F1 is on the y-axis can be created by passing the F2 values in for the `f1` argument and F1 values in for the `f2` argument, and then using the `:flip` magic argument provided by the `Plots` package.

If `meansOnly` is set to true, only the mean values for each vowel category are plotted. Using `ell=true` will plot a data ellipse that approximately encompases the percentage of data specified by `ellPercent`. The ellipse is represented by a number of points specified with `nEllPts`. Other arguments to `plot` are passed in through the splatted `kw` argument. Setting the `addLabels` argument to `true` will add the text label of the vowel category above and to the right of the mean.

Args
======

* `f1` The F1 values, or otherwise the values to plot on the x-axis
* `f2` The F2 values, or otherwise the values to plot on the y-axis
* `cats` The vowel categories associated with each F1, F2 pair
* `meansOnly` (keyword) Plot only mean value for each category
* `addLabels` (keyword) Add labels for each category to the plot near the mean
* `ell` (keyword) Whether to add data ellipses to the plot
* `ellPercent` (keyword) How much of the data distribution the ellipse should cover (approximately)
* `nEllPts` (keyword) How many points should be used when plotting the ellipse
* `markersize` (keyword) How large the markers should be; passed directly to `plot`
* `linewidth` (keyword) How wide the line for the llipses should be; passed directly to `plot`
"""
function vowelPlot(f1, f2, cats; meansOnly=false, addLabels=false, ell=false, ellPercent=0.67, nEllPts=500, markersize=3, linewidth=3, kw...)
  d = DataFrame(f1=f1, f2=f2, cat=cats)
  groups = groupby(d, :cat)

  if meansOnly
    p = Plots.scatter([mean(groups[1].f1)], [mean(groups[1].f2)], label=groups[1].cat[1],
      markersize=markersize; kw...)
  else
    p = Plots.scatter(groups[1].f1, groups[1].f2, label=groups[1].cat[1],
      markersize=markersize; kw...)
  end

  jt1 = abs(mean(f1) / 20)
  jt2 = abs(mean(f2) / 20)

  if addLabels
    p = annotate!(mean(groups[1].f1) + jt1, mean(groups[1].f2) + jt2, text(groups[1].cat[1], color=palette(:auto)[1]))
  end

  if ell
    e = ellipsePts(groups[1].f1, groups[1].f2, percent=ellPercent, nPoints=nEllPts)
    p = Plots.plot!(e[:,1], e[:,2], label="", color=1, linewidth=linewidth)
  end

  for (i, g) in enumerate(groups)
    if i == 1 continue end

    if meansOnly
      p = Plots.scatter!([mean(g.f1)], [mean(g.f2)], label=g.cat[1], color=i,
        markersize=markersize)
    else
      p = Plots.scatter!(g.f1, g.f2, label=g.cat[1], color=i, markersize=markersize)
    end

    if addLabels
      p = annotate!(mean(g.f1) + jt1, mean(g.f2) + jt2, text(g.cat[1], color=palette(:auto)[i]))
    end

    if ell
      e = ellipsePts(g.f1, g.f2, percent=ellPercent, nPoints=nEllPts)
      p = Plots.plot!(e[:,1], e[:,2], label="", color=i, linewidth=linewidth)
    end
  end
  return p
end

"""
    ellipsePts(f1, f2; percent=0.95, nPoints=500)

Calculates `nPoints` points of the perimeter of a data ellipse for `f1` and `f2` with approximately the percent of the data spcified by `percent` contained within the ellipse. Points are returned in counter-clockwise order as the polar angle of rotation moves from 0 to 2π.

See Friendly, Monette, and Fox (2013, Elliptical insights: Understanding statistical methods through elliptical geometry, Statistical science 28(1), 1-39) for more information on the calculation process.

Args
======

* `f1` The F1 values or otherwise x-axis values
* `f2` The F2 values or otherwise y-axis values
* `percent` (keyword) Percent of the data distribution the ellipse should approximately cover
* `nPoints` (keyword) How many points to use when drawing the ellipse
"""
function ellipsePts(f1, f2; percent=0.95, nPoints=500)
  Σ = cov(hcat(f1, f2))
  L = cholesky(Σ).L

  # we have two dimensions (F1 and F2), so the chisq distribution has
  # two degrees of freedom
  r = sqrt(quantile(Chisq(2), percent))
  t = range(0, 2π, length=nPoints)
  S = hcat(cos.(t), sin.(t)) .* r

  pts = L * S'
  pts[1,:] .+= mean(f1)
  pts[2,:] .+= mean(f2)
  return Array(pts')
end


text:
p_rules = Rules([0.05], ["significant", "not significant"], true)

@test interpret(0.04, p_rules) == "significant"

@test format_p(0.12, stars=true) == "p > .1"
@test format_p(0.06, stars=true) == "p = .06°"
@test format_p(0.04, stars=true) == "p < .05*"
@test format_p(0.005, stars=true) == "p < .01**"
@test format_p(0.0001, stars=true) == "p < .001***"


text:
"""
    get_variables_str(client::Client,
                      args::AbstractDict,
                      arg_names::AbstractDict,
                      query::String;
                      typedict=client.query_to_args_map)

Returns the variable string, not bracketed, for the `args` of a given `query`.
For each variable, the string is:

`\$var_name: var_type`

where `var_name` is specified in `arg_names`. This allows multiple args of the same type
to be parameterised with different names. For queries and mutations, `arg_names` is
calculated by `get_query_args_str`.

# Arguments
- `client::Client`: client.
- `args::AbstractDict`: contains name of type and value of argument as key value pairs.
    For nested arguments, the value should be a `AbstractDict` or `Vector{<:AbstractDict})`.
    See examples below.
- `arg_names::AbstractDict`: contains name of type and name to be used in variable string
    as key value pairs. The structure and keys of `arg_names` should match that of `args`.
- `query::String`: query/mutation/subscription name
- `typedict=client.query_to_args_map`: dictionary to look up type in. For top level variables,
    `query_to_args_map` is used. For lower level variables, `type_to_fields_map` is used.

# Examples
```jldoctest; setup=:(using GraphQLClient: get_variables_str, prettify_query, Client)
julia> client = Client("https://countries.trevorblades.com");

julia> args = Dict("code" => "BR");

julia> arg_names = Dict("code" => "code1");

julia> str = get_variables_str(client, args, arg_names, "country")
\"\\\$code1: ID!,\"
```
"""
function get_variables_str(client::Client,
                           args::AbstractDict,
                           arg_names::AbstractDict,
                           query::String;
                           typedict=client.query_to_args_map)

    # Query might be appended with !, so remove
    query = replace(query, r"!" => "")

    # Query might be a list containing a string, so remove outer square brackets
    query = replace(query, r"^\[|\]$" => "")

    str = ""
    for (key, val) in args
        !haskey(typedict[query], key) && throw(GraphQLClientException("Cannot query field \"$key\" on type \"$query\""))
        if isa(val, AbstractDict)
            str *= get_variables_str(client, val, arg_names[key], typedict[query][key], typedict=client.input_object_fields_to_type_map)
        elseif isa(val, Vector{<:AbstractDict})
            str *= join(get_variables_str.(Ref(client), val, arg_names[key], Ref(typedict[query][key]), typedict=client.input_object_fields_to_type_map))
        else
            str *= "\$" * arg_names[key] * ": " * typedict[query][key]
            str *= ","
        end
    end
    return str
end

"""
    get_query_variables(args, arg_names)

Returns single level dictionary with variable name and value as key val pairs.

`args` and `arg_names` must be dictionaries with the same structure and keys. `args` contains
the original variable name and its value, `arg_names` contains the original variable name and
the name being used in the query, which may be different from the original variable name, for
example if multiple of the same original name are used.

# Examples
```jldoctest; setup=:(using GraphQLClient: get_query_args_str, get_query_variables, prettify_query)
julia> args = Dict("int" => 1, "vector" => [1, 2, 3], "dict" => Dict("float" => 1.0), "vec_dict" => [Dict("bool" => true, "int" => 2)]);

julia> _, arg_names = get_query_args_str(args);

julia> arg_names
Dict{String, Any} with 4 entries:
  "int"      => "int"
  "dict"     => Dict{String, Any}("float"=>"float")
  "vec_dict" => Dict{String, Any}[Dict("int"=>"int__2", "bool"=>"bool__2")]
  "vector"   => "vector__1"

julia> get_query_variables(args, arg_names)
Dict{String, Any} with 5 entries:
  "int"       => 1
  "int__2"    => 2
  "bool__2"   => true
  "float"     => 1.0
  "vector__1" => [1, 2, 3]
```
"""
function get_query_variables(args, arg_names; name_tracking=[])
    variables = Dict{String, Any}()
    for (key, val) in args
        if val isa AbstractDict
            merge!(variables, get_query_variables(val, arg_names[key], name_tracking=name_tracking))
        elseif val isa Vector{<:AbstractDict}
            for i in eachindex(val)
                merge!(variables, get_query_variables(val[i], arg_names[key][i], name_tracking=name_tracking))
            end
        else
            # Check name is unique and add to tracking if so
            arg_names[key] in name_tracking && throw(GraphQLClientException("Duplicate name \"$(arg_names[key])\" found when creating variable dictionary for HTTP query."))
            push!(name_tracking, arg_names[key])
            variables[arg_names[key]] = val
        end
    end
    return variables
end


text:
export get_entanglement_entropy

"""
    do_schmidt_decomposition(L::Int,
                             N::Int,
                             i::Int,
                             state::Vector{<:Number})

Do Schmidt decomposition of `state` at site `i`.
"""
function do_schmidt_decomposition(L::Int,
                                  N::Int,
                                  i::Int,
                                  state::Vector{<:Number})
    # Get basis of states.
    basis = make_LN_basis(L, N)
    return do_schmidt_decomposition(L, i, basis, state)
end

"""
    do_schmidt_decomposition(L::Int,
                             i::Int,
                             basis::Vector{Int},
                             state::Vector{<:Number})

Do Schmidt decomposition of `state` at site `i`.
"""
function do_schmidt_decomposition(L::Int,
                                  i::Int,
                                  basis::Vector{Int},
                                  state::Vector{<:Number})
    # Matrix to store the decomposed state in the A and B subspaces.
    dim_A = 1<<(L-i)
    dim_B = 1<<i
    M = zeros(eltype(state), (dim_A, dim_B))

    mask = 1<<i - 1
    # s, sA, and sB refer to the indices of states in the whole space, A, and B.
    for s=1:length(basis)
        sA = basis[s]>>i
        sB = mask&basis[s]
        M[sA+1, sB+1] = state[s]
    end

    svals = svdvals(M)
    return svals
end

"""
    get_entanglement_entropy(L::Int,
                             N::Int,
                             i::Int,
                             state::Vector{<:Number})

Compute the entanglement entropy of `state` at site `i`.
"""
function get_entanglement_entropy(L::Int,
                                  N::Int,
                                  i::Int,
                                  state::Vector{<:Number})
    # Get basis of states.
    basis = make_LN_basis(L, N)
    return get_entanglement_entropy(L, i, basis, state)
end

"""
    get_entanglement_entropy(L::Int,
                             i::Int,
                             basis::Vector{Int},
                             state::Vector{<:Number})

Compute the entanglement entropy of `state` at site `i`.
"""
function get_entanglement_entropy(L::Int,
                                  i::Int,
                                  basis::Vector{Int},
                                  state::Vector{<:Number})
    svals = do_schmidt_decomposition(L, i, basis, state)
    svals2 = svals.^2
    S = 0.
    for i=1:findlast(svals2 .> 1e-14)
        S -= svals2[i]*log(svals2[i])
    end
    return S
end


text:
import TrajectoryOptimization: gen_usrfun_ipopt
using Random, LinearAlgebra, Plots
using BenchmarkTools, Statistics
Random.seed!(7)
include("N_plots.jl")

# Solver Options
dt = 0.01
integration = :rk3
method = :hermite_simpson

function convergence_plot(stat_i,stat_d;kwargs...)
    plot(log.(abs.(stat_i["cost"])),width=2,label="iLQR")
    plot!(log.(abs.(stat_d["cost"])), ylabel="Cost (abs log)",xlabel="iterations",width=2,label="DIRCOL"; kwargs...)
end

###################
## Parallel Park ##
###################

opts = SolverOptions()
opts.verbose = false
opts.cost_tolerance = 1e-6
opts.cost_tolerance_intermediate = 1e-5
opts.constraint_tolerance = 1e-5
opts.resolve_feasible = false
opts.outer_loop_update_type = :default
opts.use_nesterov = true
opts.penalty_scaling = 200
opts.penalty_initial = 1
opts.R_infeasible = 10

dircol_options = Dict("tol"=>opts.cost_tolerance,"constr_viol_tol"=>opts.constraint_tolerance)

# Set up model, objective, and solver
model, = TrajectoryOptimization.Dynamics.dubinscar
n, m = model.n,model.m

x0 = [0.0;0.0;0.]
xf = [0.0;1.0;0.]
tf =  3.
Qf = 100.0*Diagonal(I,n)
Q = (1e-3)*Diagonal(I,n)
R = (1e-2)*Diagonal(I,m)

obj = LQRObjective(Q, R, Qf, tf, x0, xf)

x_min = [-0.25; -0.001; -Inf]
x_max = [0.25; 1.001; Inf]

obj_con_box = TrajectoryOptimization.ConstrainedObjective(obj,x_min=x_min,x_max=x_max)

solver_uncon  = Solver(model, obj, integration=integration, dt=dt, opts=opts)
solver_con_box = Solver(model, obj_con_box, integration=integration, dt=dt, opts=opts)
solver_nn = Solver(model, obj_con_box, integration=integration, dt=dt, opts=opts)
solver_nn.opts.use_nesterov = false

U0 = rand(solver_uncon.model.m,solver_uncon.N)
X0 = line_trajectory(solver_con_box)
X0_rollout = rollout(solver_uncon, U0)

@time results_uncon, stats_uncon = TrajectoryOptimization.solve(solver_uncon,U0)

solver_con = Solver(model, obj_con_box, integration=integration, dt=dt, opts=opts)
solver_con.opts.use_nesterov = false
solver_con.opts.square_root = false
solver_con.opts.penalty_scaling = 100
solver_con.opts.penalty_initial = 0.01
solver_con.opts.outer_loop_update_type = :feedback
solver_con.opts.cost_tolerance_intermediate = 1e-3
solver_con.opts.constraint_tolerance_intermediate = 1
solver_con.opts.use_penalty_burnin = false
solver_con.opts.verbose = false
@time results_con_box, stats_con_box = solve(solver_con,U0)
stats_con_box["c_max"][end]

solver_con_n = Solver(solver_con)
solver_con_n.opts.use_nesterov = true
solver_con_n.opts.square_root = true
results_n, stats_n = solve(solver_con_n,U0)

plot(stats_con_box["c_max"],yscale=:log10)
plot!(stats_n["c_max"])

solver_inf = Solver(model, obj_con_box, integration=integration, dt=dt, opts=opts)
solver_inf.opts.use_nesterov = false
solver_inf.opts.penalty_scaling = 100
solver_inf.opts.penalty_initial = 0.001
solver_inf.opts.outer_loop_update_type = :feedback
solver_inf.opts.cost_tolerance_intermediate = 1e-3
solver_inf.opts.constraint_tolerance_intermediate = 1
solver_inf.opts.use_penalty_burnin = false
solver_inf.opts.verbose = false
solver_inf.opts.R_infeasible = 10
@time res_inf, stats_inf = TrajectoryOptimization.solve(solver_inf,X0,U0)
stats_inf["iterations"]

@time res_d, stats_d = TrajectoryOptimization.solve_dircol(solver_con_box, X0_rollout, U0, options=dircol_options)

# Parallel Park (boxed)
plt = plot(title="Parallel Park")#,aspect_ratio=:equal)
plot!(x_min[1]*ones(1000),collect(range(x_min[2],stop=x_max[2],length=1000)),color=:red,width=2,label="")
plot!(x_max[1]*ones(1000),collect(range(x_min[2],stop=x_max[2],length=1000)),color=:red,width=2,label="")
plot!(collect(range(x_min[1],stop=x_max[1],length=1000)),x_min[2]*ones(1000),color=:red,width=2,label="")
plot!(collect(range(x_min[1],stop=x_max[1],length=1000)),x_max[2]*ones(1000),color=:red,width=2,label="")
plot_trajectory!(to_array(results_uncon.X),width=2,color=:blue,label="Unconstrained")
plot_trajectory!(to_array(results_con_box.X),width=2,color=:green,label="Constrained",legend=:bottomright)
plot_trajectory!(to_array(res_inf.X),width=2,color=:yellow,label="Constrained (infeasible)",legend=:bottomright)
plot_trajectory!(res_d.X,width=2,color=:black,linestyle=:dash,label="DIRCOL")

# Stats comparison
eval_f, = gen_usrfun_ipopt(solver_con_box,method)
res_i = DircolVars(results_con_box)
stats_con_box["runtime"]
stats_d["runtime"]
stats_con_box["iterations"]
stats_d["iterations"]
eval_f(res_i.Z)
eval_f(res_d.Z)
TrajectoryOptimization._cost(solver_con_box,results_con_box)
plot(stats_uncon["gradient_norm"][10:end])
stats_uncon["gradient_norm"]

# Convergence Behavior
convergence_plot(stats_con_box,stats_d)
plot!(log.(abs.(stats_inf["cost"])),linewidth=2,label="iLQR (infeasible)")
convergence_rate(stats_uncon,tail=0.5,plot_fit=true)
convergence_rate(stats_con_box,tail=0.5,plot_fit=true)
convergence_rate(stats_inf,tail=0.5,plot_fit=true)

p = plot(stats_con_box["gradient_norm"],yscale=:log10,label="no nesterov",color=:blue,xlabel="iterations",ylabel="gradient",title="Constrained Parallel Park")
ylim = collect(ylims(p))
plot_vertical_lines!(stats_con_box["outer_updates"],ylim,linecolor=:blue,linestyle=:dash,label="")
plot!(stats_n["gradient_norm"],yscale=:log10,label="nesterov",linecolor=:red,legend=:bottomleft)
plot_vertical_lines!(stats_n["outer_updates"],ylim,linecolor=:red,linestyle=:dash,label="")

# STEP SIZE COMPARISONS #
# Unconstrained
Ns = [21,31,51,81,101,201,301]
disable_logging(Logging.Debug)
group = "dubinscar/parallelpark/unconstrained"
run_step_size_comparison(model, obj, U0, group, Ns, integrations=[:midpoint,:rk3],dt_truth=1e-3,opts=opts,benchmark=true)
plot_stat("runtime",group,legend=:bottomright,title="Unconstrained Parallel Park")
plot_stat("iterations",group,legend=:bottom,title="Unconstrained Parallel Park")
plot_stat("error",group,yscale=:log10,legend=:right,title="Unconstrained Parallel Park")

# Constrained
Ns = [21,31,51,81,151,201,301]
disable_logging(Logging.Debug)
group = "dubinscar/parallelpark/constrained"
run_step_size_comparison(model, obj_con_box, U0, group, Ns, integrations=[:midpoint,:rk3],dt_truth=1e-3,opts=solver_con.opts)
plot_stat("runtime",group,legend=:topleft,title="Constrained Parallel Park")
plot_stat("iterations",group,legend=:topright,title="Constrained Parallel Park")
plot_stat("error",group,yscale=:log10,legend=:topright,title="Constrained Parallel Park")
plot_stat("c_max",group,yscale=:log10,legend=:right,title="Constrained Parallel Park")

# Infeasible
Ns = [21,31,51,81,151,201,301]
disable_logging(Logging.Debug)
group = "dubinscar/parallelpark/infeasible"
run_step_size_comparison(model, obj_con_box, U0, group, Ns, integrations=[:midpoint,:rk3],dt_truth=1e-3,opts=solver_inf.opts, infeasible=true)
plot_stat("runtime",group,legend=:topleft,title="Constrained Parallel Park (infeasible)")
plot_stat("iterations",group,legend=:topright,title="Constrained Parallel Park (infeasible)")
plot_stat("error",group,yscale=:log10,legend=:topright,title="Constrained Parallel Park (infeasible)")
plot_stat("c_max",group,yscale=:log10,legend=:right,title="Constrained Parallel Park (infeasible)")

dt_truth = 1e-3
solver_truth = Solver(model, obj, dt=dt_truth)
X_truth, U_truth = get_dircol_truth(solver_truth,X0,U0,group)[2:3]
interp(t) = TrajectoryOptimization.interpolate_trajectory(solver_truth, X_truth, U_truth, t)

opts.use_nesterov = false
err, err_final, stats = run_Ns(model, obj_con_box, X0, U0, Ns, interp, opts=opts)
opts.use_nesterov = true
err_n, err_final_n, stats_n = run_Ns(model, obj_con_box, X0, U0, Ns, interp, opts=opts)
iters = [[stat["iterations"] for stat in sts] for sts in [stats,stats_n]]
runtime = [[stat["runtime"] for stat in sts] for sts in [stats,stats_n]]
plot_vals(Ns,iters,["default","nesterov"],"iterations")
plot_vals(Ns,runtime,["default","nesterov"],"runtime")

########################
## Obstacle Avoidance ##
########################

opts.penalty_initial = 0.1

x0 = [0.0;0.0;0.]
xf = [10.0;10.0;0.]
tf =  3.0
Qf = 100.0*Diagonal(I,n)
Q = (1e-3)*Diagonal(I,n)
R = (1e-2)*Diagonal(I,m)

obj = LQRObjective(Q, R, Qf, tf, x0, xf)

u_min = [-100.; -100.]
u_max = [100.; 100.]
x_min = [-Inf; -Inf; -Inf]
x_max = [Inf; Inf; Inf]

n_circles = 3
r = 1.0
circles = ((2.5,2.5,r),(5.,5.,r),(7.5,7.5,r))

function cI(c,x,u)
    for i = 1:n_circles
        c[i] = circle_constraint(x,circles[i][1],circles[i][2],circles[i][3])
    end
    c
end

obj_con_obstacles = TrajectoryOptimization.ConstrainedObjective(obj,x_min=x_min,x_max=x_max,cI=cI)
obj_con_obstacles_control = TrajectoryOptimization.ConstrainedObjective(obj,u_min=u_min, u_max=u_max,x_min=x_min,x_max=x_max,cI=cI)

solver_uncon_obstacles = Solver(model, obj, integration=integration, dt=dt, opts=opts)
solver_con_obstacles = Solver(model, obj_con_obstacles, integration=integration, dt=dt, opts=opts)
solver_con_obstacles.opts.penalty_scaling = 10
solver_con_obstacles.opts.use_nesterov = false

# -Initial state and control trajectories
U0 = rand(solver_uncon.model.m,solver_uncon_obstacles.N)
X0 = line_trajectory(solver_con_obstacles)
X0_rollout = rollout(solver_con_obstacles,U0)

@time results_uncon_obstacles, stats_uncon_obstacles = TrajectoryOptimization.solve(solver_uncon_obstacles,U0)
@time results_con_obstacles, stats_con_obstacles = TrajectoryOptimization.solve(solver_con_obstacles,U0)
results_inf, stats_inf = TrajectoryOptimization.solve(solver_con_obstacles,X0,U0)
res_d, stats_d = solve_dircol(solver_con_obstacles,X0_rollout,U0,options=dircol_options)

# Obstacle Avoidance
plt = plot(title="Obstacle Avoidance")
plot_obstacles(circles)
plot_trajectory!(to_array(results_uncon_obstacles.X),width=2,color=:blue,label="Unconstrained",aspect_ratio=:equal,xlim=[-1,11],ylim=[-1,11])
plot_trajectory!(to_array(results_con_obstacles.X),width=2,color=:green,label="Constrained")
plot_trajectory!(to_array(results_con_obstacles.X),width=2,color=:orange,label="Infeasible")
plot_trajectory!(res_d.X,width=2,color=:black,linestyle=:dash,label="DIRCOL",legend=:bottomright)
display(plt)

group
solver_truth = Solver(model,obj_con_obstacles,dt=1e-3)
_,X_truth,U_truth = get_dircol_truth(solver_truth,X0,U0,group)
plot_trajectory!(X_truth,width=2,color=:grey,linestyle=:dash,label="DIRCOL (truth)",legend=:bottomright)
res_truth = DircolVars(X_truth,U_truth)

# Stats comparison
eval_f_truth, = gen_usrfun_ipopt(solver_truth,method)
eval_f, = gen_usrfun_ipopt(solver_con_obstacles,method)
res_i = DircolVars(to_array(results_con_obstacles.X),to_array(results_con_obstacles.U))
res_inf = DircolVars(results_inf)
stats_con_obstacles["runtime"]
stats_d["runtime"]
stats_con_obstacles["iterations"]
stats_d["iterations"]
convergence_plot(stats_con_obstacles,stats_d)
plot!(stats_inf["cost"],width=2,label="Infeasible")
eval_f(res_i.Z)
eval_f(res_d.Z)
eval_f(res_inf.Z)
eval_f_truth(res_truth.Z)
TrajectoryOptimization._cost(solver_con_obstacles,results_con_obstacles)
TrajectoryOptimization._cost(solver_con_obstacles,results_inf)

# KNOT POINT COMPARISON
Ns = [21,41,51,81,101,201,301]
opts.penalty_scaling = 10
opts.use_nesterov = false
disable_logging(Logging.Debug)
group = "dubinscar/obstacles/constrained"
run_step_size_comparison(model, obj_con_obstacles, U0, group, Ns, integrations=[:midpoint,:rk3],dt_truth=1e-3,opts=opts)
plot_stat("runtime",group,legend=:topleft,title="Obstacle Avoidance")
plot_stat("iterations",group,legend=:left,title="Obstacle Avoidance")
plot_stat("error",group,yscale=:log10,legend=:topright,title=" Obstacle Avoidance")
plot_stat("c_max",group,yscale=:log10,legend=:right,title="Obstacle Avoidance")

Ns = [21,41,51,81,101,201,301]
opts.penalty_scaling = 10
opts.use_nesterov = false
disable_logging(Logging.Debug)
group = "dubinscar/obstacles/infeasible"
run_step_size_comparison(model, obj_con_obstacles, U0, group, Ns, integrations=[:midpoint,:rk3],dt_truth=1e-3,opts=opts,infeasible=true)
plot_stat("runtime",group,legend=:topleft,title="Obstacle Avoidance")
plot_stat("iterations",group,legend=:left,title="Obstacle Avoidance")
plot_stat("error",group,yscale=:log10,legend=:topright,title=" Obstacle Avoidance")
plot_stat("c_max",group,yscale=:log10,legend=:right,title="Obstacle Avoidance")

# Nesterov comparison
opts.use_nesterov = false
err, err_final, stats = run_Ns(model, obj_con_obstacles, X0, U0, Ns, interp, opts=opts, infeasible=true)
opts.use_nesterov = true
err_n, err_final_n, stats_n = run_Ns(model, obj_con_obstacles, X0, U0, Ns, interp, opts=opts, infeasible=true)
iters = [[stat["iterations"] for stat in sts] for sts in [stats,stats_n]]
runtime = [[stat["runtime"] for stat in sts] for sts in [stats,stats_n]]
plot_vals(Ns,iters,["default","nesterov"],"iterations")
plot_vals(Ns,runtime,["default","nesterov"],"runtime")

############
## Escape ##
############

x0 = [2.5;2.5;0.]
xf = [7.5;2.5;0.]
tf =  3.0
Qf = 100.0*Diagonal(I,n)
Q = (1e-3)*Diagonal(I,n)
R = (1e-2)*Diagonal(I,m)

obj = LQRObjective(Q, R, Qf, tf, x0, xf)

r = 0.5
s1 = 30; s2 = 50; s3 = 15
circles = []
for i in range(0,stop=5,length=s1)
    push!(circles,(0.,i,r))
end
for i in range(0,stop=5,length=s1)
    push!(circles,(5.,i,r))
end
for i in range(0,stop=5,length=s1)
    push!(circles,(10.,i,r))
end
for i in range(0,stop=10,length=s2)
    push!(circles,(i,0.,r))
end
for i in range(0,stop=3,length=s3)
    push!(circles,(i,5.,r))
end
for i in range(5,stop=8,length=s3)
    push!(circles,(i,5.,r))
end

n_circles = 3*s1 + s2 + 2*s3

# plt = plot(title="Escape",aspect_ratio=:equal,xlim=[-1,11],ylim=[-1,11])
# plot_obstacles(circles,:red)
# display(plt)

function cI(c,x,u)
    for i = 1:n_circles
        c[i] = circle_constraint(x,circles[i][1],circles[i][2],circles[i][3])
    end
    c
end

obj_con_obstacles = TrajectoryOptimization.ConstrainedObjective(obj,cI=cI,u_min=-10,u_max=10)

solver_uncon_obstacles = Solver(model, obj, dt=dt, opts=opts)
solver_con_obstacles = Solver(model, obj_con_obstacles, dt=dt, opts=opts)
solver_con_obstacles.opts.R_infeasible = 1.0
n,m,N = get_sizes(solver_con_obstacles)

# -Initial state and control trajectories
X_guess = [2.5 2.5 0.;4. 5. .785;5. 6.25 0.;7.5 6.25 -.261;9 5. -1.57;7.5 2.5 0.]
X0 = TrajectoryOptimization.interp_rows(N,tf,Array(X_guess'))
# X0 = line_trajectory(solver_uncon)
U0 = ones(m,N)

# @time results_uncon_obstacles, stats_uncon_obstacles = TrajectoryOptimization.solve(solver_uncon_obstacles,U0)
# @time results_con_obstacles, stats_con_obstacles = TrajectoryOptimization.solve(solver_con_obstacles,U0)
solver = Solver(model, obj_con_obstacles, N=101, opts=SolverOptions())
solver.opts.R_infeasible = 1
solver.opts.square_root
solver.opts.resolve_feasible = false
solver.opts.cost_tolerance = 1e-6
solver.opts.cost_tolerance_intermediate = 1e-3
solver.opts.constraint_tolerance = 1e-5
solver.opts.constraint_tolerance_intermediate = 0.01
solver.opts.penalty_scaling = 50
solver.opts.penalty_initial = 10
solver.opts.outer_loop_update_type = :default
solver.opts.use_penalty_burnin = false
solver.opts.use_nesterov = false
solver.opts.live_plotting = true
solver.opts.verbose = true
@time res_inf, stats_inf = solve(solver,X0,U0)
stats_inf["iterations"]

res_d,stats_d = solve_dircol(solver,X0,U0,options=dircol_options)

# Escape
plt = plot(title="Escape",aspect_ratio=:equal)
plot_obstacles(circles)
# plot_trajectory!(to_array(results_uncon_obstacles.X),width=2,color=:blue,label="Unconstrained",aspect_ratio=:equal,xlim=[-1,11],ylim=[-1,11])
# plot_trajectory!(to_array(results_con_obstacles.X),width=2,color=:green,label="Constrained",aspect_ratio=:equal,xlim=[-1,11],ylim=[-1,11])
plot_trajectory!(to_array(res_inf.X),width=2,color=:purple,label="Infeasible",aspect_ratio=:equal,xlim=[-1,11],ylim=[-1,11])
plot_trajectory!(res_d.X,width=2,color=:yellow,label="DIRCOL",aspect_ratio=:equal,xlim=[-1,11],ylim=[-1,11])
plot!(X0[1,:],X0[2,:],label="Infeasible Initialization",width=1,color=:purple,linestyle=:dash)
display(plt)

# Dircol Truth
group = "dubinscar/escape/infeasible"
solver_truth = Solver(model,obj_con_obstacles,N=1001)
solve_dircol(solver_truth,Array(res_d.X), Array(res_d.U))
run_dircol_truth(solver_truth, Array(res_d.X), Array(res_d.U), group::String)

Ns = [101,151,201,251,301]
disable_logging(Logging.Debug)
run_step_size_comparison(model, obj_con_obstacles, U0, group, Ns, integrations=[:midpoint,:rk3],dt_truth=solver_truth.dt,opts=solver.opts,infeasible=true,X0=X0)
plot_stat("runtime",group,legend=:topleft,title="Obstacle Avoidance",yscale=:log10)
plot_stat("iterations",group,legend=:right,title="Obstacle Avoidance",yscale=:log10)
plot_stat("error",group,yscale=:log10,legend=:topright,title=" Obstacle Avoidance")
plot_stat("c_max",group,yscale=:log10,legend=:right,title="Obstacle Avoidance")

X_truth, U_truth = get_dircol_truth(solver_truth,X0,U0,group)[2:3]
interp(t) = TrajectoryOptimization.interpolate_trajectory(solver_truth, X_truth, U_truth, t)

opts_no_nesterov = copy(solver.opts)
opts_no_nesterov.use_nesterov = false
opts_no_nesterov.outer_loop_update_type = :feedback
err, err_final, stats = run_Ns(model, obj_con_box, X0, U0, Ns, interp, opts=opts_no_nesterov)
opts_no_nesterov.use_nesterov = true
err_n, err_final_n, stats_n = run_Ns(model, obj_con_box, X0, U0, Ns, interp, opts=opts_no_nesterov)
iters = [[stat["iterations"] for stat in sts] for sts in [stats,stats_n]]
runtime = [[stat["runtime"] for stat in sts] for sts in [stats,stats_n]]
plot_vals(Ns,iters,["default","nesterov"],"iterations")
plot_vals(Ns,runtime,["default","nesterov"],"runtime")

plot(stats[2]["gradient_norm"],yscale=:log10)
plot!(stats_n[2]["gradient_norm"])


text:
"""
    AbstractTimeSpan

A type repesenting a continuous, inclusive span between two points in time.

All subtypes of `AbstractTimeSpan` must implement:

- `first(::AbstractTimeSpan)::Nanosecond`: return the first nanosecond contained in `span`
- `last(::AbstractTimeSpan)::Nanosecond`: return the last nanosecond contained in `span`

For convenience, many Onda functions that accept `AbstractTimeSpan` values also accept
`Dates.Period` values.

See also: [`TimeSpan`](@ref)
"""
abstract type AbstractTimeSpan end

Base.first(t::Period) = convert(Nanosecond, t)

Base.last(t::Period) = convert(Nanosecond, t)

function _validate_timespan(first::Nanosecond, last::Nanosecond)
    if first > last
        throw(ArgumentError("start of time span should precede end, got $first and $last"))
    end
    return nothing
end

#####
##### `TimeSpan`
#####

"""
    TimeSpan(first, last)

Return `TimeSpan(Nanosecond(first), Nanosecond(last))::AbstractTimeSpan`.

See also: [`AbstractTimeSpan`](@ref)
"""
struct TimeSpan <: AbstractTimeSpan
    first::Nanosecond
    last::Nanosecond
    function TimeSpan(first::Nanosecond, last::Nanosecond)
        _validate_timespan(first, last)
        return new(first, last)
    end
    TimeSpan(first, last) = TimeSpan(Nanosecond(first), Nanosecond(last))
end

"""
    TimeSpan(x)

Return `TimeSpan(first(x), last(x))`.

See also: [`AbstractTimeSpan`](@ref)
"""
TimeSpan(x) = TimeSpan(first(x), last(x))

Base.first(span::TimeSpan) = span.first

Base.last(span::TimeSpan) = span.last

#####
##### `AbstractTimeSpan` Utilities
#####

"""
    contains(a, b)

Return `true` if the timespan `b` lies entirely within the timespan `a`, return `false` otherwise.
"""
function contains(a, b)
    a, b = TimeSpan(a), TimeSpan(b)
    return first(a) <= first(b) && last(a) >= last(b)
end

"""
    overlaps(a, b)

Return `true` if the timespan `a` and the timespan `b` overlap, return `false` otherwise.
"""
function overlaps(a, b)
    a, b = TimeSpan(a), TimeSpan(b)
    starts_earlier, starts_later = ifelse(first(b) > first(a), (a, b), (b, a))
    return last(starts_earlier) >= first(starts_later)
end

"""
    shortest_timespan_containing(spans)

Return the shortest possible `TimeSpan` containing all timespans in `spans`.

`spans` is assumed to be an iterable of timespans.
"""
function shortest_timespan_containing(spans)
    isempty(spans) && throw(ArgumentError("input iterator must be nonempty"))
    lo, hi = Nanosecond(typemax(Int64)), Nanosecond(0)
    for span in spans
        lo = min(first(span), lo)
        hi = max(last(span), hi)
    end
    return TimeSpan(lo, hi)
end

"""
    duration(span)

Return the duration of `span` as a `Period`.

For `span::AbstractTimeSpan`, this is equivalent to `last(span) - first(span)`.

For `span::Period`, this function is the identity.
"""
duration(t::AbstractTimeSpan) = last(t) - first(t)
duration(t::Period) = t

nanoseconds_per_sample(sample_rate) = inv(sample_rate) * 1_000_000_000

"""
    index_from_time(sample_rate, sample_time)

Given `sample_rate` in Hz, return the integer index of the most recent sample
taken at `sample_time`. Note that `sample_time` must be non-negative and support
`convert(Nanosecond, sample_time)`.

Examples:

```
julia> index_from_time(1, Second(0))
1

julia> index_from_time(1, Second(1))
2

julia> index_from_time(100, Millisecond(999))
100

julia> index_from_time(100, Millisecond(1000))
101
```
"""
function index_from_time(sample_rate, sample_time)
    time_in_nanoseconds = convert(Nanosecond, sample_time).value
    time_in_nanoseconds >= 0 || throw(ArgumentError("`sample_time` must be >= 0 nanoseconds"))
    ns_per_sample = nanoseconds_per_sample(sample_rate)
    return floor(Int, time_in_nanoseconds / ns_per_sample) + 1
end

"""
    index_from_time(sample_rate, span::AbstractTimeSpan)

Return the `UnitRange` of indices corresponding to `span` given `sample_rate` in Hz:

```
julia> index_from_time(100, TimeSpan(Second(0), Second(1)))
1:100

julia> index_from_time(100, TimeSpan(Second(1), Second(1)))
101:101

julia> index_from_time(100, TimeSpan(Second(3), Second(6)))
301:600
```
"""
function index_from_time(sample_rate, span::AbstractTimeSpan)
    i = index_from_time(sample_rate, first(span))
    j = index_from_time(sample_rate, last(span))
    j = i == j ? j : (j - 1)
    return i:j
end

"""
    time_from_index(sample_rate, sample_index)

Given `sample_rate` in Hz and assuming `sample_index > 0`, return the earliest
`Nanosecond` containing `sample_index`.

Examples:

```
julia> time_from_index(1, 1)
0 nanoseconds

julia> time_from_index(1, 2)
1000000000 nanoseconds

julia> time_from_index(100, 100)
990000000 nanoseconds

julia> time_from_index(100, 101)
1000000000 nanoseconds
```
"""
function time_from_index(sample_rate, index)
    index > 0 || throw(ArgumentError("`index` must be > 0"))
    return Nanosecond(ceil(Int, (index - 1) * nanoseconds_per_sample(sample_rate)))
end

"""
    time_from_index(sample_rate, sample_range::AbstractUnitRange)

Return the `TimeSpan` corresponding to `sample_range` given `sample_rate` in Hz:

```
julia> time_from_index(100, 1:100)
TimeSpan(0 nanoseconds, 1000000000 nanoseconds)

julia> time_from_index(100, 101:101)
TimeSpan(1000000000 nanoseconds, 1000000000 nanoseconds)

julia> time_from_index(100, 301:600)
TimeSpan(3000000000 nanoseconds, 6000000000 nanoseconds)
```
"""
function time_from_index(sample_rate, sample_range::AbstractUnitRange)
    i, j = first(sample_range), last(sample_range)
    j = j == i ? j : j + 1
    return TimeSpan(time_from_index(sample_rate, i),
                    time_from_index(sample_rate, j))
end


text:
# Hello World してみる
println("Hello World")

max_stars_repo_name: tetsugps/julia

max_issues_repo_name: tetsugps/julia

max_forks_repo_name: tetsugps/julia


text:
module api

include("types.jl")

const paths = Sys.isapple() ? String["/System/Library/Frameworks/OpenCL.framework"] : String[]

import Libdl

const libopencl = Libdl.find_library(["libOpenCL", "OpenCL"], paths)
@assert libopencl != ""

function _ocl_func(func, ret_type, arg_types)
    local args_in = Symbol[Symbol("arg$i")
                           for (i, T) in enumerate(arg_types.args)]

    esc(quote
        function $func($(args_in...))
            ccall(($(string(func)), libopencl),
                   $ret_type,
                   $arg_types,
                   $(args_in...))
        end
    end)
end

macro ocl_func(func, ret_type, arg_types)
    _ocl_func(func, ret_type, arg_types)
end

const CL_callback  = Ptr{Nothing}

abstract type CL_user_data_tag end
const CL_user_data = Ptr{CL_user_data_tag}

Base.cconvert(::Type{Ptr{CL_user_data_tag}}, obj::T) where {T} = Ref{T}(obj)

Base.unsafe_convert(P::Type{Ptr{CL_user_data_tag}}, ptr::Ref) = P(Base.unsafe_convert(Ptr{Cvoid}, ptr))
Base.unsafe_convert(P::Type{Ptr{CL_user_data_tag}}, ptr::Ptr) = P(Base.unsafe_convert(Ptr{Cvoid}, ptr))

include("api/opencl_1.0.0.jl")
include("api/opencl_1.1.0.jl")
include("api/opencl_1.2.0.jl")
include("api/opencl_2.0.0.jl")

function parse_version(version_string)
    mg = match(r"^OpenCL ([0-9]+)\.([0-9]+) .*$", version_string)
    if mg === nothing
        error("Non conforming version string: $(ver)")
    end
    return VersionNumber(parse(Int, mg.captures[1]),
                                 parse(Int, mg.captures[2]))
end

end


text:
#                        .-'''-.                               _..._
#                       '   _    \         _______          .-'_..._''.
#  /|                 /   /` '.   \        \  ___ `'.     .' .'      '.\
#  ||                .   |     \  '         ' |--.\  \   / .'
#  ||        .-,.--. |   '      |  '        | |    \  ' . '                                     .|
#  ||  __    |  .-. |\    \     / /  __     | |     |  '| |                 __                .' |_
#  ||/'__ '. | |  | | `.   ` ..' /.:--.'.   | |     |  || |              .:--.'.         _  .'     |
#  |:/`  '. '| |  | |    '-...-'`/ |   \ |  | |     ' .'. '             / |   \ |      .' |'--.  .-'
#  ||     | || |  '-             `" __ | |  | |___.' /'  \ '.          .`" __ | |     .   | / |  |
#  ||\    / '| |                  .'.''| | /_______.'/    '. `._____.-'/ .'.''| |   .'.'| |// |  |
#  |/\'..' / | |                 / /   | |_\_______|/       `-.______ / / /   | |_.'.'.-'  /  |  '.'
#  '  `'-'`  |_|                 \ \._,\ '/                          `  \ \._,\ '/.'   \_.'   |   /
#                                 `--'  `"                               `--'  `"             `'-'

using Base.Broadcast
using Base.Broadcast: Broadcasted, AbstractArrayStyle, DefaultArrayStyle, broadcasted,
  instantiate, materialize, flatten, combine_eltypes, _broadcast_getindex

# Structural utilities

using Base: tail

tcat(x) = x
tcat(x, y, z...) = tcat((x..., y...), z...)

broadcast_args(x) = (x,)
broadcast_args(bc::Broadcasted) = tcat(map(broadcast_args, bc.args)...)

_unflatten(x, xs) = first(xs), tail(xs)

_unflatten(x::Tuple{}, xs) = (), xs

function _unflatten(x::Tuple, xs)
  t1, xs1 = _unflatten(first(x), xs)
  t2, xs2 = _unflatten(tail(x), xs1)
  (t1, t2...), xs2
end

function _unflatten(bc::Broadcasted, xs)
  t, xs′ = _unflatten(bc.args, xs)
  (args=t,f=nothing,axes=nothing), xs′
end

unflatten(x, xs) = _unflatten(x, xs)[1]

unflatten(x, xs::Nothing) = nothing

trim(x, Δ) = reshape(Δ, ntuple(i -> size(Δ, i), Val(ndims(x))))

unbroadcast(x::AbstractArray, Δ) =
  size(x) == size(Δ) ? Δ :
  length(x) == length(Δ) ? trim(x, Δ) :
    trim(x, sum(Δ, dims = ntuple(i -> size(x, i) == 1 ? i : ndims(Δ)+1, Val(ndims(Δ)))))

unbroadcast(x::Union{Number,Ref}, Δ) = sum(Δ)

# Reverse Mode

Jtrivial(f, a...) = nothing
Jtrivial(::typeof(+), a...) = a
Jtrivial(::typeof(-), a, b) = (a..., .-b...)

trivia(_) = (1,)
function trivia(bc::Broadcasted)
  t = map(trivia, bc.args)
  any(t -> t === nothing, t) && return
  Jtrivial(bc.f, t...)
end

Joutput(f, a...) = nothing
Joutput(::typeof(exp), x) = map(t -> y -> y*t, x)

function Jbroadcast(bc::Broadcasted)
  t = map(trivia, bc.args)
  any(t -> t === nothing, t) && return
  Joutput(bc.f, t...)
end

@inline function unbroadcast_r(x, y, ȳ, j::J) where J
  trim(x, j.(y).*ȳ)
end

@inline function unbroadcast_r(x::Number, y, ȳ, j::J) where J
  x̄ = zero(float(x))
  @simd for I in eachindex(y)
    @inbounds x̄ += j(y[I])*ȳ[I]
  end
  return x̄
end

function ∇broadcast_r(bc::Broadcasted, J)
  y = copy(instantiate(bc))
  back(ȳ) = map(unbroadcast_r, broadcast_args(bc), map(_ -> y, J), map(_ -> ȳ, J), J)
  return y, back
end

# Mixed Mode

import ForwardDiff
using ForwardDiff: Dual

dual(x, p) = x
dual(x::Real, p) = Dual(x, p)

dualtype(::Type{Dual{G,T,P}}) where {G,T,P} = T
dualtype(T) = T

function dual_function(f::F) where F
  function (args::Vararg{Any,N}) where N
    ds = map(args, ntuple(identity,Val(N))) do x, i
      dual(x, ntuple(j -> i==j, Val(N)))
    end
    return f(ds...)
  end
end

dualify(bc::Broadcasted{S}) where S = Broadcasted{S}(dual_function(bc.f), bc.args, bc.axes)

@inline function broadcast_gradient!(bc::Broadcasted, dest::AbstractArray, grads...)
  @simd for I in eachindex(bc)
    @inbounds begin
      out = bc[I]
      dest[I] = ForwardDiff.value(out)
      Δs = out isa Dual ? out.partials.values : map(_ -> false, grads)
      map((g, p) -> g[I] = p, grads, Δs)
    end
  end
end

function broadcast_gradient(bc::Broadcasted, ::Type{T}) where T
  dest = similar(bc, T)
  grads = map(_ -> similar(bc, promote_type(T,Bool)), bc.args)
  broadcast_gradient!(bc, dest, grads...)
  return dest, grads
end

@inline function ∇broadcast_f(bc′::Broadcasted)
  bc = dualify(instantiate(flatten(bc′)))
  T = combine_eltypes(bc.f, bc.args)
  T <: Bool && return copy(bc′), _ -> nothing
  y, gs = broadcast_gradient(bc, dualtype(T))
  back(Δ) = map((x, d) -> unbroadcast(x, Δ.*d), bc.args, gs)
  return y, back
end

function ∇broadcast_f(bc::Broadcasted{<:AbstractArrayStyle{0}})
  out = dualify(instantiate(flatten(bc)))[]
  return out.value, Δ -> map(x -> x*Δ, out.partials.values)
end

∇broadcast(bc::Broadcasted, ::Nothing) = ∇broadcast_f(bc)
∇broadcast(bc::Broadcasted, J) = ∇broadcast_r(bc, J)
∇broadcast(bc::Broadcasted) = ∇broadcast(bc, Jbroadcast(bc))

@adjoint function broadcasted(f, args...)
  broadcasted(f, args...), Δ -> (nothing, Δ.args...)
end

@adjoint function materialize(bc::Broadcasted{<:DefaultArrayStyle})
  y, back = ∇broadcast(bc)
  y, Δ -> (unflatten(bc, back(Δ)),)
end


text:
__precompile__(true)

# TODO x/0 = 0
# TODO ϕ * 2, ϕ1 + 5, etc...
# TODO ϕ .< 3, ϕ1 == ϕ2, etc... (how would this work?)
# TODO assignment tuple indexing, e.g. ϕ["a", 1]
# TODO push(ϕ, dims...), push(ϕ, Pair{Symbol}...)
# TODO better / faster broadcast
# TODO custom error type for NegativeError
# TODO also, better error type of others (e.g. Bounds error, DimensionMismatch, ...)
# TODO broadcast_reduce tag-team
# TODO map_reduce as well
# TODO add @boundscheck and @inbounds where applicable

module Factors

import Base: .==, .!=, .<, .<=, .>, .>=, in, ==, *, /, +, -, hash
import DataFrames

export
    # Dimension stuff
    Dimension,
    RangeDimension,

    name,
    eltype,
    support,
    spttype,
    update,

    Assignment,

    Factor,
    # basic access
    scope,
    potential,
    lengths,
    # interact with dimensions
    getdim,
    push,
    pattern,
    pattern_states,
    getdim,
    reducedim!,
    Z,
    marginalize,
    # assignment conversions
    at2sub,
    sub2at,
    at2a,
    a2at,
    a2sub,
    sub2a,

    # negatives stuff
    set_negative_mode,
    NegativeMode,
    NegativeIgnore,
    NegativeWarn,
    NegativeError

include("dimensions.jl")   # define Dimension
include("factors_main.jl") # define Factors
include("auxiliary.jl")    # define Assignment

include("factors_broadcast.jl")
include("factors_index.jl")
include("factors_iter.jl")
include("factors_join.jl")
include("factors_map.jl")
include("factors_reduce.jl")

include("dataframes.jl")
include("errors.jl")
include("io.jl")
include("negatives.jl")

end # module


text:
using Documenter, WebDAV

makedocs(;
         modules=[WebDAV],
         format=Documenter.HTML(),
         pages=[
             "index.md",
         ],
         repo="https://github.com/Alexander-Barth/WebDAV.jl/blob/{commit}{path}#L{line}",
         sitename="WebDAV",
         authors="Alexander Barth",
)

deploydocs(;
    repo="github.com/Alexander-Barth/WebDAV.jl",
)


text:
struct KSH <: AbstractShell
    path::AbstractString
    prompt::AbstractString
end

KSH(path::AbstractString) = KSH(path, "(playground)> ")
KSH() = KSH(strip(readstring(`which ksh`)))

function Base.run(shell::KSH, env::Environment)
    init(env)
    prompt = getprompt(shell, env)
    ENV["PS1"] = prompt

    usr_rc = join(Path(home()), ".kshrc")
    pg_rc = join(parent(Path(ENV["JULIA_PKGDIR"])), ".kshrc")

    if !exists(pg_rc)
        exists(usr_rc) ? cp(usr_rc, pg_rc, follow_symlinks=true) : touch(pg_rc)

        content = string(
            "\nexport PATH=\"", ENV["PATH"], "\"\n",
            "export PS1=\"$prompt\"\n",
            "export JULIA_PKGDIR=\"", ENV["JULIA_PKGDIR"], "\"\n",
        )

        if haskey(ENV, "HISTFILE")
            content = string(content, "export HISTFILE=\"", ENV["HISTFILE"], "\"\n")
        end

        write(pg_rc, content, "a")
    end

    ENV["ENV"] = pg_rc
    @mock run(`$(shell.path) -i`)
end


text:
## Credit to Simon Danisch for most of the following code

function AbstractPlotting.to_vertices(cells::AbstractVector{<: JuAFEM.Node{N, T}}) where {N, T}
    convert_arguments(nothing, cells)[1]
end

function AbstractPlotting.to_gl_indices(cells::AbstractVector{<: JuAFEM.Cell})
    tris = GLTriangle[]
    for cell in cells
        to_triangle(tris, cell)
    end
    tris
end

function to_triangle(tris, cell::Union{JuAFEM.Hexahedron, JuAFEM.QuadraticHexahedron})
    nodes = cell.nodes
    push!(tris, GLTriangle(nodes[1], nodes[2], nodes[5]))
    push!(tris, GLTriangle(nodes[5], nodes[2], nodes[6]))

    push!(tris, GLTriangle(nodes[6], nodes[2], nodes[3]))
    push!(tris, GLTriangle(nodes[3], nodes[6], nodes[7]))

    push!(tris, GLTriangle(nodes[7], nodes[8], nodes[3]))
    push!(tris, GLTriangle(nodes[3], nodes[8], nodes[4]))

    push!(tris, GLTriangle(nodes[4], nodes[8], nodes[5]))
    push!(tris, GLTriangle(nodes[5], nodes[4], nodes[1]))

    push!(tris, GLTriangle(nodes[1], nodes[2], nodes[3]))
    push!(tris, GLTriangle(nodes[3], nodes[1], nodes[4]))
end

function to_triangle(tris, cell::Union{JuAFEM.Tetrahedron, JuAFEM.QuadraticTetrahedron})
    nodes = cell.nodes
    push!(tris, GLTriangle(nodes[1], nodes[3], nodes[2]))
    push!(tris, GLTriangle(nodes[3], nodes[4], nodes[2]))
    push!(tris, GLTriangle(nodes[4], nodes[3], nodes[1]))
    push!(tris, GLTriangle(nodes[4], nodes[1], nodes[2]))
end

function to_triangle(tris, cell::Union{JuAFEM.Quadrilateral, JuAFEM.QuadraticQuadrilateral})
    nodes = cell.nodes
    push!(tris, GLTriangle(nodes[1], nodes[2], nodes[3]))
    push!(tris, GLTriangle(nodes[3], nodes[4], nodes[1]))
end

function to_triangle(tris, cell::Union{JuAFEM.Triangle, JuAFEM.QuadraticTriangle})
    nodes = cell.nodes
    push!(tris, GLTriangle(nodes[1], nodes[2], nodes[3]))
end

function AbstractPlotting.convert_arguments(P, x::AbstractVector{<: JuAFEM.Node{N, T}}) where {N, T}
    convert_arguments(P, reinterpret(Point{N, T}, x))
end

function visualize(mesh::JuAFEM.Grid{dim}, u) where {dim}
    T = eltype(u)
    nnodes = length(mesh.nodes)
    #TODO make this work without creating a Node
    if dim == 2
        nodes = broadcast(1:nnodes, mesh.nodes) do i, node
            JuAFEM.Node(ntuple(Val{3}) do j
                if j < 3
                    node.x[j]
                else
                    zero(T)
                end
            end)
        end
        u = [u; zeros(T, 1, nnodes)]
    else
        nodes = mesh.nodes
    end

    cnode = AbstractPlotting.Node(zeros(Float32, length(mesh.nodes)))
    scene = Makie.mesh(nodes, mesh.cells, color = cnode, colorrange = (0.0, 33.0), shading = false);
    mplot = scene[end]
    displacevec = reinterpret(GeometryTypes.Vec{3, Float64}, u, (size(u, 2),))
    displace = norm.(displacevec)
    new_nodes = broadcast(1:length(mesh.nodes), nodes) do i, node
        JuAFEM.Node(ntuple(Val{3}) do j
            node.x[j] + u[j, i]
        end)
    end
    mesh!(nodes, mesh.cells, color = (:gray, 0.4))

    scatter!(Point3f0.(getfield.(new_nodes, :x)), markersize = 0.1);
    # TODO make mplot[1] = new_nodes work
    mplot.input_args[1][] = new_nodes
    # TODO make mplot[:color] = displace work
    push!(cnode, displace)
    points = reinterpret(Point{3, Float64}, nodes)
    #arrows!(points, displacevec, linecolor = (:black, 0.3))

    scene
end

function visualize(problem::LinearElasticityProblem{dim, T}, u) where {dim, T}
    mesh = problem.ch.dh.grid
    node_dofs = problem.metadata.node_dofs
    nnodes = JuAFEM.getnnodes(mesh)
    node_displacements = reshape(u[node_dofs], dim, nnodes)
    visualize(mesh, node_displacements)
end

function visualize(problem::LinearElasticityProblem{dim, T}) where {dim, T}
    mesh = problem.ch.dh.grid
    nnodes = JuAFEM.getnnodes(mesh)
    node_displacements = zeros(T, dim, nnodes)
    visualize(mesh, node_displacements)
end

function visualize(mesh::JuAFEM.Grid{dim, N, T}) where {dim, N, T}
    nnodes = JuAFEM.getnnodes(mesh)
    node_displacements = zeros(T, dim, nnodes)
    visualize(mesh, node_displacements)
end


text:
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:light
#     text_representation:
#       extension: .jl
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.11.3
#   kernelspec:
#     display_name: Julia 1.6.0
#     language: julia
#     name: julia-1.6
# ---

# + [markdown] slideshow={"slide_type": "slide"}
# # [JuliaCon 2021](https://juliacon.org/2021/) Workshop
#
# **Title:** Modeling Marine Ecosystems At Multiple Scales Using Julia
#
# **Speakers:** [Gael Forget](https://github.com/gaelforget), [Benoit Pasquier](https://github.com/briochemc), [Zhen Wu](https://github.com/zhenwu0728)
#
# - material : <https://github.com/JuliaOcean/MarineEcosystemsJuliaCon2021.jl>
# - streaming : <https://www.youtube.com/watch?v=UCIRrXz2ZS0>
# - Q & A : <https://pigeonhole.at/MARINE>

# + [markdown] slideshow={"slide_type": "subslide"}
# Life in the oceans is strongly connected to our climate. In this workshop, you will learn to use packages from the [JuliaOcean](https://github.com/JuliaOcean) and [JuliaClimate](https://github.com/JuliaClimate) organizations that provide a foundation for studying marine ecosystems across a wide range of scales. We will run agent-based models to explore individual microbes and processes that drive species interactions. On the other end of the model hierarchy, we will simulate planetary-scale transports that control ocean biogeography and climate change.

# + [markdown] slideshow={"slide_type": "slide"}
# ## Schedule
#
# - Introduction of the topics covered, presenters, installation, and workshop roadmap (15 minutes).
#
# - [AIBECS.jl](https://github.com/JuliaOcean/AIBECS.jl): concept, implementation, tutorial workthough (30 minutes + 10' for questions; [this notebook URL](https://raw.githubusercontent.com/JuliaOcean/MarineEcosystemsJuliaCon2021.jl/main/src/AIBECSExample.jl))
#
# - [PlanktonIndividuals.jl](https://github.com/JuliaOcean/PlanktonIndividuals.jl): concept, implementation, tutorial workthough (30 minutes + 10' for questions; [this notebook URL](https://raw.githubusercontent.com/JuliaOcean/MarineEcosystemsJuliaCon2021.jl/main/src/PlanktonIndividualExample.jl))
#
# - [ClimateModels.jl](https://github.com/gaelforget/ClimateModels.jl) and [MITgcmTools.jl](https://github.com/gaelforget/MITgcmTools.jl): concept, implementation, tutorial workthough (20 minutes + 10' for questions; [this notebook URL](https://raw.githubusercontent.com/JuliaOcean/MarineEcosystemsJuliaCon2021.jl/main/src/MITgcm_tutorial_global_oce_biogeo.jl))
#
# - [ClimateModels.jl](https://github.com/gaelforget/ClimateModels.jl) and [IndividualDisplacements.jl](https://github.com/JuliaClimate/IndividualDisplacements.jl): concept, implementation, tutorial workthough (20 minutes + 10' for questions; [this notebook URL](https://raw.githubusercontent.com/JuliaOcean/MarineEcosystemsJuliaCon2021.jl/main/src/IndividualDisplacementsExample.jl))
#
# - Q&A, tutorials, etc wrap-up
#

# + [markdown] slideshow={"slide_type": "slide"}
# ## Context
#
# - marine ecosystems and climate:
#     - from microbes O($10^{-6}m$) to whole Earth O($10^6$m); from seconds to millenia
#     - life in a moving fluid; interactions between physics, chemistry, and biology
#     - many fascinating and crucial questions w.r.t. climate change (carbon, food, biodiversity, ...)
# - hierarchy of models :
#     - single ODE to massively parallel HPC simulations
#     - various languages; julia as the future and to unify

# + [markdown] slideshow={"slide_type": "-"} cell_style="split"
#
#
#

# + [markdown] slideshow={"slide_type": "-"} cell_style="split"
#
#
#

# + [markdown] slideshow={"slide_type": "subslide"}
# - Github Organizations :
#     - JuliaOcean
#     - JuliaClimate
#     - ...
# - Funding / Support :
#     - NASA (PO, MAP, IDS)
#     - Simons Foundation (CBIOMES, Scope, Gradients)
#     - ...

# + [markdown] slideshow={"slide_type": "slide"}
# ## How To
#
# To run the notebooks of this workshop on your machine, you need to:
#
# 1. **Install julia** (v1.6.0 or later, <https://julialang.org/>)
#
# 1. **Start julia**
#
# 1. **Add Pluto.jl** (v0.15.0 or later, https://github.com/fonsp/Pluto.jl)
#
#     ```julia
#     import Pkg
#     Pkg.add("Pluto")
#     ```
#
# 1. **Use Pluto to run a notebook**
#
#     ```julia
#     using Pluto
#     Pluto.run(notebook="https://raw.githubusercontent.com/JuliaOcean/MarineEcosystemsJuliaCon2021.jl/main/src/AIBECSExample.jl")
#     ```

# + [markdown] slideshow={"slide_type": "subslide"}
# **Alternatively**, instead of your own computer, you can just launch a Pluto instance in the cloud using [JuliaHub.com](https://juliahub.com/ui/Home), paste a notebook URL in the [Pluto start page](https://github.com/fonsp/Pluto.jl), and click open.

# + [markdown] slideshow={"slide_type": "slide"}
# ## Schedule
#
# - Introduction of the topics covered, presenters, installation, and workshop roadmap (15 minutes).
#
# - [AIBECS.jl](https://github.com/JuliaOcean/AIBECS.jl): concept, implementation, tutorial workthough (30 minutes + 10' for questions; [this notebook URL](https://raw.githubusercontent.com/JuliaOcean/MarineEcosystemsJuliaCon2021.jl/main/src/AIBECSExample.jl))
#
# - [PlanktonIndividuals.jl](https://github.com/JuliaOcean/PlanktonIndividuals.jl): concept, implementation, tutorial workthough (30 minutes + 10' for questions; [this notebook URL](https://raw.githubusercontent.com/JuliaOcean/MarineEcosystemsJuliaCon2021.jl/main/src/PlanktonIndividualExample.jl))
#
# - [ClimateModels.jl](https://github.com/gaelforget/ClimateModels.jl) and [MITgcmTools.jl](https://github.com/gaelforget/MITgcmTools.jl): concept, implementation, tutorial workthough (20 minutes + 10' for questions; [this notebook URL](https://raw.githubusercontent.com/JuliaOcean/MarineEcosystemsJuliaCon2021.jl/main/src/MITgcm_tutorial_global_oce_biogeo.jl))
#
# - [ClimateModels.jl](https://github.com/gaelforget/ClimateModels.jl) and [IndividualDisplacements.jl](https://github.com/JuliaClimate/IndividualDisplacements.jl): concept, implementation, tutorial workthough (20 minutes + 10' for questions; [this notebook URL](https://raw.githubusercontent.com/JuliaOcean/MarineEcosystemsJuliaCon2021.jl/main/src/IndividualDisplacementsExample.jl))
#
# - Q&A, tutorials, etc wrap-up
#


text:
###############################################################################
#
# Set X and Y limits for the axis
#
###############################################################################
"""
    set_xlims!(axs::Array, xlims::Array)
    set_xlims!(axs::Array,
               xlims::Union{Array{Int,1},Array{Float32,1},Array{Float64,1}}
    )

Set the X-axis limits, given
- `axs` Array of axis
- `xlims` Array of X-axis limits
"""
function set_xlims!(axs::Array, xlims::Array)
    # set X-axis limit for each panel
    for indx in eachindex(axs)
        axs[indx].set_xlim(xlims[indx]);
    end

    return nothing
end

function set_xlims!(
            axs::Array,
            xlims::Union{Array{Int,1},Array{Float32,1},Array{Float64,1}}
)
    # set X-axis limit for each panel
    for indx in eachindex(axs)
        axs[indx].set_xlim(xlims);
    end

    return nothing
end

"""
    set_ylims!(axs::Array, ylims::Array)
    set_ylims!(axs::Array,
               ylims::Union{Array{Int,1},Array{Float32,1},Array{Float64,1}}
    )

Set the Y-axis limits, given
- `axs` Array of axis
- `ylims` Array of Y-axis limits
"""
function set_ylims!(axs::Array, ylims::Array)
    # set Y-axis limit for each panel
    for indx in eachindex(axs)
        axs[indx].set_ylim(ylims[indx]);
    end

    return nothing
end

function set_ylims!(
            axs::Array,
            ylims::Union{Array{Int,1},Array{Float32,1},Array{Float64,1}}
)
    # set X-axis limit for each panel
    for indx in eachindex(axs)
        axs[indx].set_ylim(ylims);
    end

    return nothing
end

"""
    set_xylims!(axs::Array, xlims::Array, ylims::Array)

Set the X-axis limits, given
- `axs` Array of axis
- `xlims` Array of X-axis limits
- `ylims` Array of Y-axis limits
"""
function set_xylims!(axs::Array, xlims::Array, ylims::Array)
    set_xlims!(axs, xlims);
    set_ylims!(axs, ylims);

    return nothing
end


text:
# Evaluates dual objective function and its gradient
function dualobj!(ncm, G, H, L, τ, α, method;
    computeV::Bool = false,
    scaleX::Bool = true
    )

    n = ncm.n
    d = ncm.d
    g = ncm.g

    M = ncm.M
    R = ncm.R
    H2 = ncm.H2
    Y = ncm.Y
    ∇fY = ncm.∇fY
    Γ = ncm.Γ
    V = ncm.V
    X = ncm.X
    Z = ncm.Z
    Rd = ncm.Rd

    proj = ncm.proj
    res = ncm.res

    Xnew = res.X
    y = res.y
    Λ = res.Λ

    rpRef = res.rpRef
    rdRef = res.rdRef

    fvals = res.fvals
    resvals = res.resvals
    distvals = res.distvals

    res.fgcountRef[] += 1
    fgcount = res.fgcountRef[]

    if method == :IR
        Lbar = L / τ
    elseif method == :IER
        Lbar = 1/α + L
    else
        Lbar = L
    end
    invLbar = 1/Lbar

    # ∇fY = H2.*(Y - G)
    # M = Y - (τ/L)*(∇f(Y) - Diag(y))
    # X = M
    @inbounds for j = 1:n
        for i = 1:j
            ∇fY.data[i, j] = H2.data[i, j] * (Y.data[i, j] - G.data[i, j])
            M.data[i, j] = Y.data[i, j] - invLbar * ∇fY.data[i, j]
        end
        M.data[j, j] += invLbar * y[j]
    end

    proj(M, X, Λ)

    if scaleX
        # Ensure that diag(Xnew).==1 exactly
        @inbounds for j = 1:n
            if X.data[j, j] > 0.0
                d[j] = 1.0 / sqrt(X.data[j, j])
            else
                d[j] = 1.0
            end
        end
    end

    # Λ  = Lbar*(X - M)        # Λ is psd
    # Γ  = -Diag(y) - Λ        # Γ is an ε-subgradient
    # Z  = Xnew - Y            # used to compute ||Xnew - Y||
    # V  = ∇f(Y) + Lbar*Z + Γ  # used in the IR and IER updates
    # R  = H.*(Xnew - G)       # used to compute primal objective
    # Rd = H.*R + Γ            # used to compute dual feasibility

    @inbounds for j = 1:n
        for i = 1:j
            if scaleX
                Xnew.data[i, j] = d[i] * d[j] * X.data[i, j]
            else
                Xnew.data[i, j] = X.data[i, j]
            end
            #Λ.data[i,j] = Lbar*(X.data[i,j] - M.data[i,j])
            Λ.data[i, j] *= Lbar
            Γ.data[i, j] = -Λ.data[i, j]
            Z.data[i, j] = Xnew.data[i, j] - Y.data[i, j]
            if computeV
                #V.data[i, j] = ∇fY.data[i, j] + Lbar * Z.data[i, j] + Γ.data[i, j]
                V.data[i,j] = Lbar * (Xnew.data[i,j] - X.data[i,j])
            end
            R.data[i, j] = H.data[i, j] * (Xnew.data[i, j] - G.data[i, j])
            Rd.data[i, j] = H.data[i, j] * R.data[i, j] + Γ.data[i, j]
        end
        Γ.data[j, j] -= y[j]
        #if computeV
        #    V.data[j, j] -= y[j]
        #end
        Rd.data[j, j] -= y[j]
    end

    # Compute primal and dual feasibility
    @inbounds for j = 1:n
        d[j] = Xnew.data[j, j] - 1.0
    end
    rpRef[] = norm(d)
    rdRef[] = fronorm(Rd, proj.work)

    fvals[fgcount] = 0.5 * fronorm(R, proj.work)^2
    resvals[fgcount] = max(rpRef[], rdRef[])
    distvals[fgcount] = fronorm(Z, proj.work)

    # Compute the gradient of the dual function
    @inbounds for j = 1:n
        g[j] = X.data[j, j] - 1.0
    end

    # Compute the value of the dual function
    λ = proj.w
    dualobjval = 0.0
    @inbounds for j = 1:n
        tmp = λ[j]
        if tmp > 0.0
            dualobjval += tmp^2
        end
    end
    dualobjval *= 0.5 * Lbar
    dualobjval -= sum(y)

    return dualobjval
end

max_stars_repo_name: krislock/IR-FISTA

max_issues_repo_name: krislock/IR-FISTA

max_forks_repo_name: krislock/IR-FISTA