text:
module UnitSystems

#   This file is part of UnitSystems.jl
#   It is licensed under the MIT license
#   UnitSystems Copyright (C) 2020 Michael Reed
#       _           _                         _
#      | |         | |                       | |
#   ___| |__   __ _| | ___ __ __ ___   ____ _| | __ _
#  / __| '_ \ / _` | |/ / '__/ _` \ \ / / _` | |/ _` |
# | (__| | | | (_| |   <| | | (_| |\ V / (_| | | (_| |
#  \___|_| |_|\__,_|_|\_\_|  \__,_| \_/ \__,_|_|\__,_|
#
#   https://github.com/chakravala
#   https://crucialflow.com
#  _  _ _  _ _ ___ ____ _   _ ____ ___ ____ _  _ ____
#  |  | |\ | |  |  [__   \_/  [__   |  |___ |\/| [__
#  |__| | \| |  |  ___]   |   ___]  |  |___ |  | ___]

import Base: @pure, length, time, angle, rem

const Systems = (:Metric,:SI2019,:SI1976,:CODATA,:Conventional,:International,:InternationalMean,:MetricEngineering,:SI2019Engineering,:GravitationalMetric,:GravitationalSI2019,:MTS,:EMU,:ESU,:Gauss,:LorentzHeaviside,:Kennelly,:FPS,:IPS,:British,:English,:Survey,:FFF,:MPH,:KKH,:Nautical,:Meridian,:MeridianEngineering,:GravitationalMeridian,:IAU☉,:IAUE,:IAUJ,:Hubble,:Cosmological,:CosmologicalQuantum,:Planck,:PlanckGauss,:Stoney,:Hartree,:Rydberg,:Schrodinger,:Electronic,:Natural,:NaturalGauss,:QCD,:QCDGauss,:QCDoriginal)
const Dimensionless = (:coupling,:finestructure,:electronunit,:protonunit,:protonelectron,:darkenergydensity)
const Constants = (:lightspeed,:planck,:planckreduced,:electronmass,:molarmass,:boltzmann,:permeability,:rationalization,:lorentz,:luminousefficacy,:gravity) #angle
const Physics = (:turn,:spat,:atomicmass,:protonmass,:planckmass,:gravitation,:gaussgravitation,:einstein,:hartree,:rydberg,:bohr,:electronradius,:avogadro,:molargas,:stefan,:radiationdensity,:vacuumpermeability,:vacuumpermittivity,:electrostatic,:magnetostatic,:biotsavart,:elementarycharge,:faraday,:vacuumimpedance,:conductancequantum,:klitzing,:josephson,:magneticfluxquantum,:magneton)
const Derived = (:hyperfine,:loschmidt,:wienwavelength,:wienfrequency,:mechanicalheat,:solarmass,:jupitermass,:earthmass,:lunarmass,:earthradius,:greatcircle,:radarmile,:hubble,:cosmological,
    :radian,:steradian,:degree,:gradian,:arcminute,:arcsecond,
    :second,:minute,:hour,:day,:gaussianmonth,:siderealmonth,:synodicmonth,:year,:gaussianyear,:siderealyear,:jovianyear,
    :angstrom,:inch,:foot,:surveyfoot,:yard,:meter,:earthmeter,:mile,:statutemile,:meridianmile,:admiraltymile,:nauticalmile,:lunardistance,:astronomicalunit,:jupiterdistance,:lightyear,:parsec,
    :barn,:hectare,:acre,:surveyacre,
    :liter,:gallon,:quart,:pint,:cup,:fluidounce,:teaspoon,:tablespoon,
    :bubnoff,:ips,:fps,:fpm,:ms,:kmh,:mph,:knot,:mps,
    :grain,:gram,:earthgram,:kilogram,:tonne,:ton,:pound,:ounce,:slug,:slinch,:hyl,
    :dyne,:newton,:poundal,:poundforce,:kilopond,
    :psi,:pascal,:bar,:barye,:technicalatmosphere,:atmosphere,:inchmercury,:torr,
    :electronvolt,:erg,:joule,:footpound,:calorie,:kilocalorie,:meancalorie,:earthcalorie,:thermalunit,:gasgallon,:tontnt,
    :watt,:horsepower,:horsepowerwatt,:horsepowermetric,:electricalhorsepower,:tonsrefrigeration,:boilerhorsepower,
    :coulomb,:earthcoulomb,:ampere,:volt,:henry,:ohm,:siemens,:farad,:weber,:tesla,
    :abcoulomb,:abampere,:abvolt,:abhenry,:abohm,:abmho,:abfarad,:maxwell,:gauss,:oersted,:gilbert,
    :statcoulomb,:statampere,:statvolt,:stathenry,:statohm,:statmho,:statfarad,:statweber,:stattesla,
    :kelvin,:rankine,:celsius,:fahrenheit,:sealevel,:boiling,:mole,:earthmole,:poundmole,:slugmole,:slinchmole,:katal,:amagat,
    :lumen,:candela,:lux,:phot,:footcandle,:nit,:apostilb,:stilb,:lambert,:footlambert,:bril,
    :neper,:bel,:decibel,:hertz,:apm,:rpm,
    :kayser,:diopter,:gforce,:galileo,:eotvos,:darcy,:poise,:reyn,:stokes,:rayl,
    :mpge,:langley,:jansky,:solarflux,:curie,:sievert,:roentgen,:rem)
const Kinematic = (:solidangle,:time,:length,:area,:volume,:wavenumber,:angularwavenumber,:fuelefficiency,:numberdensity,:frequency,:angularfrequency,:frequencydrift,:speed,:acceleration,:jerk,:snap,:crackle,:pop,:volumeflow) #angle
const Mechanical = (:inertia,:mass,:massflow,:lineardensity,:areadensity,:density,:specificweight,:specificvolume,:force,:specificforce,:gravityforce,:pressure,:compressibility,:viscosity,:diffusivity,:rotationalinertia,:impulse,:momentum,:angularmomentum,:yank,:energy,:specificenergy,:action,:fluence,:power,:powerdensity,:intensity,:spectralflux,:soundexposure,:impedance,:specificimpedance,:admittance,:compliance,:inertance)
const Electromagnetic = (:charge,:chargedensity,:linearchargedensity,:exposure,:mobility,:current,:currentdensity,:resistance,:conductance,:resistivity,:conductivity,:capacitance,:inductance,:reluctance,:permeance,:permittivity,:permeability,:susceptibility,:specificsusceptibility,:demagnetizingfactor,:vectorpotential,:electricpotential,:magneticpotential,:electricfield,:magneticfield,:electricflux,:magneticflux,:electricfluxdensity,:magneticfluxdensity,:electricdipolemoment,:magneticdipolemoment,:electricpolarizability,:magneticpolarizability,:magneticmoment,:specificmagnetization,:polestrength)
const Thermodynamic = (:temperature,:entropy,:specificentropy,:volumeheatcapacity,:thermalconductivity,:thermalconductance,:thermalresistivity,:thermalresistance,:thermalexpansion,:lapserate)
const Molar = (:molarmass,:molality,:molaramount,:molarity,:molarvolume,:molarentropy,:molarenergy,:molarconductivity,:molarsusceptibility,:catalysis,:specificity)
const Photometric = (:luminousflux,:luminousintensity,:luminance,:illuminance,:luminousenergy,:luminousexposure,:luminousefficacy)
const Mechanics = [Kinematic...,Mechanical...]
const Convert = [Mechanics...,Electromagnetic...,Thermodynamic...,Molar...,Photometric...]

listext(x) = join(x,"`, `")

"""
    UnitSystems.similitude() = haskey(ENV,"SIMILITUDE")

An optional environment variable `ENV["SIMILITUDE"]` induces `UnitSystems.similitude()` to return `true`, giving flexibility for building dependencies whenever it is desirable to toggle usage between `UnitSystems` (default) and `Similitude` (requires environment variable specification). For example, in `MeasureSystems` and `Geophysics` this option is used to increase flexibility with variety in local compilation workflow.
"""
similitude() = haskey(ENV,"SIMILITUDE")

include("constant.jl")
Constant(D::Number) = Constant{D}()
@pure Constant(D::Float64) = Constant{D}()
@pure Constant(D::Int) = Constant{D}()
@pure Constant(D::Constant) = D
@pure constant(::Constant{D}) where D = D
Base.Int(::Constant{D}) where D = Constant(Int(D))
Base.show(io::IO,x::Constant{D}) where D  = show(io,D)
unit(x,::Constant{y}) where y = unit(x,y)
unit(::Constant{x},::Constant{y}) where {x,y} = Constant{unit(x,y)}()
unit(::Constant{x},y=1) where x = Constant{unit(x,y)}()

@pure measure(x) = x
@pure cache(x) = x
Constant(x) = x
Quantity(x) = x
Quantity(D,U,x) = x
const 𝟙,F,M,L,T,Q,Θ,N,J,A,Λ,C = 1,1,1,1,1,1,1,1,1,1,1,1

logdb(x) = 10log10(x)
expdb(x) = exp10(0.1)^x
const dB, Db = logdb, expdb

# universe

"""
    Coupling{αG,α,μₑᵤ,μₚᵤ,ΩΛ}

Specification of `Universe` with the dimensionless `Coupling` constants `coupling`, `finestructure`, `electronunit`, `protonunit`, `protonelectron`, and `darkenergydensity`.
Alterations to these values can be facilitated and quantified using parametric polymorphism.
Due to the `Coupling` interoperability, the `MeasureSystems` package is made possible to support calculations with `Measurements` having error standard deviations.
"""
struct Coupling{αG,α,μₑᵤ,μₚᵤ,ΩΛ} end
Coupling{αG,α,μₑᵤ,μₚᵤ}() where {αG,α,μₑᵤ,μₚᵤ} = Coupling(αG,α,μₑᵤ,μₚᵤ,ΩΛ)
Coupling(αG,α,μₑᵤ,μₚᵤ,ΩΛ) = Coupling{cache(αG),cache(α),cache(μₑᵤ),cache(μₚᵤ),cache(ΩΛ)}()
@pure coupling(U::Coupling{αG}) where αG = measure(αG)
@pure finestructure(U::Coupling{αG,α}) where {αG,α} = measure(α)
@pure electronunit(U::Coupling{αG,α,μₑᵤ}) where {αG,α,μₑᵤ} = measure(μₑᵤ)
@pure protonunit(U::Coupling{αG,α,μₑᵤ,μₚᵤ}) where {αG,α,μₑᵤ,μₚᵤ} = measure(μₚᵤ)
@pure protonelectron(U::Coupling) = protonunit(U)/electronunit(U)
@pure darkenergydensity(U::Coupling{αG,α,μₑᵤ,μₚᵤ,ΩΛ}) where {αG,α,μₑᵤ,μₚᵤ,ΩΛ} = measure(ΩΛ)
Base.display(U::Coupling) = println("Coupling{αG = $(coupling(U)), α = $(finestructure(U)), μₑᵤ = $(electronunit(U)), μₚᵤ = $(protonunit(U)), ΩΛ = $(darkenergydensity(U))}")

# unit systems

"""
    UnitSystem{kB, ħ, 𝘤, μ₀, mₑ, Mᵤ, (Kcd, θ, λ, αL, g₀, ...)}

Fundamental constants of physics are: `kB` Boltzmann's constant, `ħ` reduced Planck's constant, `𝘤` speed of light, `μ₀` vacuum permeability, `mₑ` electron rest mass, `Mᵤ` molar mass, `Kcd` luminous efficacy, `θ` angle measure, `λ` Gauss rationalization, `αL` Lorentz's constant, and `g₀` gravitational force reference.
Primarily the `Metric` SI unit system is used in addition to the historic `English` engineering unit system.
These constants induce derived values for `avogadro`, `boltzmann`, `molargas`, `planck`, `planckreduced`, `lightspeed`, `planckmass`, `atomicmass`, `protonmass`, `electronmass`, `newton`, `einstein`, `vacuumpermeability`, `vacuumpermittivity`, `electrostatic`, and
additional constants `molarmass`, `luminousefficacy`, `gravity`, `angle`, `turn`, `spat`, `stefan`, `radiationdensity`, `magnetostatic`, `lorentz`, `biotsavart`, `rationalization`, `vacuumimpedance`, `elementarycharge`, `magneton`, `conductancequantum`, `faraday`, `magneticfluxquantum`, `josephson`, `klitzing`, `hartree`, `rydberg`, `bohr`.

Standardized unit/derived quantities are `$(listext(Derived))`.

Additional reference `UnitSystem` variants: `EMU`, `ESU`, `Gauss`, `LorentzHeaviside`, `SI2019`, `SI1976`, `CODATA`, `Conventional`, `International`, `InternationalMean`, `MetricEngineering`, `GravitationalMetric`, `IAU`, `IAUE`, `IAUJ`, `FPS`, `IPS`, `British`, `Survey`, `Hubble`, `Cosmological`, `CosmologicalQuantum`, `Meridian`, `Nautical`, `MPH`, `KKH`, `MTS`, `FFF`; and natural atomic units based on gravitational `coupling` and `finestructure` constant (`Planck`, `PlanckGauss`, `Stoney`, `Hartree`, `Rydberg`, `Schrodinger`, `Electronic`, `Natural`, `NaturalGauss`, `QCD`, `QCDGauss`, and `QCDoriginal`).

**Dimensional unit conversions:**

Mechanics: `angle`, `$(listext(Kinematic))`, `$(listext(Mechanical))`;
Electromagnetics: `$(listext(Electromagnetic))`;
Thermodynamics: `$(listext(Thermodynamic))`,
`$(listext(Molar))`, `$(listext(Photometric))`.
""" #`Rᵤ,mᵤ,σ,ħ,μ₀,ε₀,kₑ,𝘦,𝔉,RK,Z₀,G₀`
struct UnitSystem{kB,ħ,𝘤,μ₀,mₑ,Mᵤ,extra}
    @pure UnitSystem{kB,ħ,𝘤,μ₀,mₑ,Mᵤ,extra}() where {kB,ħ,𝘤,μ₀,mₑ,Mᵤ,extra} = new{kB,ħ,𝘤,μ₀,mₑ,Mᵤ,extra}()
end # UnitSystem{kB,ħ,𝘤,μ₀,mₑ,Mᵤ,(Kcd,θ,λ,αL,g,C,τ,𝟐,𝟑,𝟓,𝟕,𝟏𝟏,𝟏𝟗,𝟒𝟑)}
function UnitSystem(kB,ħ,𝘤,μ₀,mₑ,Mᵤ=1,Kcd=1,θ=1,λ=1,αL=1,g=1,C=Universe,τ=2π,x=2,y=3,z=5,w=7,u=11,v=19,q=43)
    UnitSystem{cache(kB),cache(ħ),cache(𝘤),cache(μ₀),cache(mₑ),cache(Mᵤ),(cache(Kcd),cache(θ),cache(λ),cache(αL),cache(g),C,cache(τ),cache(x),cache(y),cache(z),cache(w),cache(u),cache(v),cache(q))}()
end
@pure boltzmann(::UnitSystem{k}) where k = measure(k)
@pure planckreduced(::UnitSystem{k,ħ}) where {k,ħ} = measure(ħ)
@pure lightspeed(::UnitSystem{k,ħ,𝘤}) where {k,ħ,𝘤} = measure(𝘤)
@pure permeability(::UnitSystem{k,ħ,𝘤,μ}) where {k,ħ,𝘤,μ} = measure(μ)
@pure vacuumpermeability(U::UnitSystem) = permeability(U)
@pure electronmass(::UnitSystem{k,ħ,𝘤,μ,m}) where {k,ħ,𝘤,μ,m} = measure(m)
@pure molarmass(::UnitSystem{k,ħ,𝘤,μ,m,M}) where {k,ħ,𝘤,μ,m,M} = measure(M)
@pure luminousefficacy(::UnitSystem{k,ħ,𝘤,μ,m,M,e}) where {k,ħ,𝘤,μ,m,M,e} = measure(e[1])
@pure angle(::UnitSystem{k,ħ,𝘤,μ,m,M,e}) where {k,ħ,𝘤,μ,m,M,e} = measure(e[2])
@pure rationalization(::UnitSystem{k,ħ,𝘤,μ,m,M,e}) where {k,ħ,𝘤,μ,m,M,e} = measure(e[3])
@pure lorentz(::UnitSystem{k,ħ,𝘤,μ,m,M,e}) where {k,ħ,𝘤,μ,m,M,e} = measure(e[4])
@pure gravity(::UnitSystem{k,ħ,𝘤,μ,m,M,e}) where {k,ħ,𝘤,μ,m,M,e} = measure(e[5])
@pure universe(::UnitSystem{k,ħ,𝘤,μ,m,M,e}) where {k,ħ,𝘤,μ,m,M,e} = measure(e[6])
@pure twopi(::UnitSystem{k,ħ,𝘤,μ,m,M,e}) where {k,ħ,𝘤,μ,m,M,e} = measure(e[7])
@pure two(::UnitSystem{k,ħ,𝘤,μ,m,M,e}) where {k,ħ,𝘤,μ,m,M,e} = measure(e[8])
@pure three(::UnitSystem{k,ħ,𝘤,μ,m,M,e}) where {k,ħ,𝘤,μ,m,M,e} = measure(e[9])
@pure five(::UnitSystem{k,ħ,𝘤,μ,m,M,e}) where {k,ħ,𝘤,μ,m,M,e} = measure(e[10])
@pure seven(::UnitSystem{k,ħ,𝘤,μ,m,M,e}) where {k,ħ,𝘤,μ,m,M,e} = measure(e[11])
@pure eleven(::UnitSystem{k,ħ,𝘤,μ,m,M,e}) where {k,ħ,𝘤,μ,m,M,e} = measure(e[12])
@pure nineteen(::UnitSystem{k,ħ,𝘤,μ,m,M,e}) where {k,ħ,𝘤,μ,m,M,e} = measure(e[13])
@pure fourtythree(::UnitSystem{k,ħ,𝘤,μ,m,M,e}) where {k,ħ,𝘤,μ,m,M,e} = measure(e[14])
# ΔνCs:s⁻¹, c:m⋅s⁻¹, h:kg⋅m²⋅s⁻¹, kB:kg⋅m²⋅s⁻²⋅K⁻¹, NA:mol⁻¹, Kcd: cd⋅sr⋅s³⋅kg⁻¹⋅m⁻²

function evaldim end
@pure isquantity(U) = false
@pure isquantity(A,B) = isquantity(A) && isquantity(B)
@pure isquantity(U::UnitSystem) = isquantity(boltzmann(U))
@pure isrationalized(U::UnitSystem) = rationalization(U) ≠ sphere(U)

normal(x) = x
@pure normal(x::Float64) = x
@pure normal(x::Int) = x
@pure normal(::UnitSystem{kB,ħ,𝘤,μ0,me,Mu,e}) where {kB,ħ,𝘤,μ0,me,Mu,e} = UnitSystem(normal(kB),normal(ħ),normal(𝘤),normal(μ0),normal(me),normal(Mu),normal(e[1]),normal(e[2]),normal(e[3]),normal(e[4]),normal(e[5]),e[6],normal(e[7]),normal(e[8]),normal(e[9]),normal(e[10]),normal(e[11]),normal(e[12]),normal(e[13]),normal(e[14]))

@pure unitname(::UnitSystem) = "Unknown"
Base.show(io::IO,U::UnitSystem) = print(io,unitname(normal(U)))
function Base.display(U::UnitSystem)
    println("UnitSystem: ", unitname(normal(U)))
    println("  kB  = $(boltzmann(U))")
    println("  ħ   = $(planckreduced(U))")
    println("  𝘤   = $(lightspeed(U))")
    println("  μ₀  = $(vacuumpermeability(U))")
    println("  mₑ  = $(electronmass(U))")
    println("  Mᵤ  = $(molarmass(U))")
    println("  Kcd = $(luminousefficacy(U))")
    println("  θ   = $(angle(U))")
    println("  λ   = $(rationalization(U)≠4π ? rationalization(U) : "4π")")
    println("  αL  = $(lorentz(U))")
    println("  g₀  = $(gravity(U))")
end

function (U::UnitSystem)(JK=1,Js=1,ms=1,Hm=1,kg=1)
    kB = boltzmann(U)*JK
    ħ = planckreduced(U)*Js
    𝘤 = lightspeed(U)*ms
    μ₀ = vacuumpermeability(U)*Hm
    mₑ = electronmass(U)*kg
    Mᵤ = molarmass(U)#kg#
    Kcd = luminousefficacy(U)#
    A = angle(U)
    λ = rationalization(U)
    αL = lorentz(U)
    g₀ = gravity(U)
    C = universe(U)
    τ = twopi(U)
    (x,y,z,u,v,w,p) = (two(U),three(U),five(U),seven(U),eleven(U),nineteen(U),fourtythree(U))
    UnitSystem(kB,ħ,𝘤,μ₀,mₑ,Mᵤ,Kcd,A,λ,isone(αL) ? αL : αL/ms,g₀,C,τ,x,y,z,u,v,w,p)
end

function DimensionSystem(U::UnitSystem,L,M,T,I,Θ,N,J,A,Λ,G,C)
    kB = boltzmann(U)*(M*L*L/(T*T)/Θ/G)
    ħ = reducedplanck(U)*(M*(L*L)/T/G/A)
    c = lightpseed(U)*(L/T)
    μ0 = vacuumpermeability(U)*(L*M/((T*T)*(I*I)*(C*C)*(A*A)*Λ))
    me = electronmass(U)*M
    Mu = molarmass(U)*(M/N)
    Kcd = luminousefficacy(U)*(J*(T*T*T)/(L*L*M)*G)
    λ = rationalization(U)*(Λ*(A*A))
    αL = lorentz(U)*C
    g₀ = gravity(U)*G
    τ = twopi(U)
    (x,y,z,u,v,w,p) = (two(U),three(U),five(U),seven(U),eleven(U),fourtythree(U))
    UnitSystem(kB,ħ,c,μ0,me,Mu,Kcd,A,λ,αL,g₀,universe(U),τ,x,y,z,u,v,w,p)
end

"""
    AstronomicalSystem(U::UnitSystem) = EntropySystem(U,𝟏,𝟏,𝟏/gravitation(U))

Constructs new `UnitSystem` from `U` with `mass` rescaled by `newton` gravitational constant. An example of this type is `AstronomicalSystem(Metric)`.
"""
AstronomicalSystem(u) = EntropySystem(u,one(u),one(u),inv(gravitation(u)))

"""
    ElectricSystem(U::UnitSystem,Ω,V) = EntropySystem(U,𝟏,𝟏,V^2/Ω,𝟏,vacuumpermeability(U)/Ω)

Constructs new `UnitSystem` from `U` with `mass` rescaled by `electricpotential` and `resistance`. In the `International` system, `Ωᵢₜ` and `Vᵢₜ` are used as definitions from the more recent United States results, while in `InternationalMean` an earlier estimate based on other nations was used.
"""
ElectricSystem(u,Ω,V) = EntropySystem(u,one(u),one(u),V*V/Ω,one(u),vacuumpermeability(u)/Ω)

"""
    GaussSystem(U::UnitSystem,μ0,λ,αL=𝟏,l=centi,m=milli,g0=gravity(U))

Constructs new `UnitSystem` from `U` rescaled for `CGS` with electromagnetic options. The first three options are to set the values for `vacuumpermeability`, `rationalization`, and `lorentz` constants. The following two parameters are scaling for `length` and `mass`, while the last is an option to change the `gravity` reference.

Examples include `EMU`, `ESU`, `Gauss`, `LorentzHeaviside`, and `Kennelly`.
"""
function GaussSystem(u,μ0,λ,αL=one(u),l=inv((two(u)*five(u))^2),m=inv((two(u)*five(u))^3),g0=gravity(u))
    EntropySystem(u,one(u),l,m,one(u),μ0,m==1/1000 ? one(u) : molarmass(u)/m,g0,m*(l*l),λ,αL)
end

"""
    EntropySystem(U::UnitSystem,t,l,m,θ=𝟏)
    EntropySystem(U::UnitSystem,t,l,m,θ,μ0,Mu=molarmass(U)/m,g0=gravity(U))

Constructs new `UnitSystem` from `U` rescaled along `time`, `length`, `mass`, and `temperature` by the first four parameters. Additional optional parameters allow for customization of the `vacuumpermeability`, `molarmass`, and `gravity` constants.

Examples of this type include `Nautical`, `Meridian`, `MeridianEngineering`, `GravitatonalMeridian`, `GravitationalMetric`, `GravitationalSI2019`, `MTS`, `KKH`, `MPH`, `IAU☉`, `IAUE`, `IAUJ`, `Hubble`, `Cosmological`, `CosmologicalQuantum`.
However, most other constructors for `UnitSystem` derivations are based on internally calling `EntropySystem`, such as `AstronomicalSystem`, `ElectricSystem`, `GaussSystem`, and `RankineSystem`.
This means `EntropySystem` also constructs the examples listed there.
"""
function EntropySystem(u,t,l,m,θ=one(u))
    EntropySystem(u,t,l,m,θ,permeability(u)/(m*l),molarmass(u)/m,gravity(u),m*l*l/(t*t))
end
function EntropySystem(u,t,l,m,θ,μ0,Mu=molarmass(u)/m,g0=gravity(u),e=m*l*l/(t*t),λ=one(u),αL=one(u))
    UnitSystem(
        boltzmann(u)*θ/e/g0,
        planckreduced(u)/t/e/g0,
        lightspeed(u)*t/l,
        μ0,
        electronmass(u)/m,
        Mu,
        luminousefficacy(u)*e/t*g0,
        angle(u),λ,αL,g0,universe(u),twopi(u),two(u),three(u),five(u),seven(u),eleven(u),nineteen(u),fourtythree(u))
end

@pure unit(x,y=1) = isapprox(y,x,rtol=eps()^0.9) ? y : x
@pure Base.one(U::UnitSystem) = unit(two(U)/two(U))
@pure Base.zero(U::UnitSystem) = one(U)-one(U)
@pure turn(U::UnitSystem) = twopi(U)/angle(U)
@pure solidangle(U::UnitSystem,S::UnitSystem) = unit(angle(U,S)^2)
@pure spat(U::UnitSystem) = two(U)*turn(U)/angle(U)*unit(turn(U)/normal(turn(U)))
@pure mass(U::UnitSystem,S::UnitSystem) = electronmass(U,S)
@pure electronmass(𝘩::Number,C::Coupling) = inv(finestructure(C))^2*R∞*2𝘩/𝘤
@pure planckmass(U::UnitSystem,C::Coupling=universe(U)) = electronmass(U,C)/√coupling(C)
@pure planck(U::UnitSystem,C::Coupling=universe(U)) = turn(U)*planckreduced(U,C)
@pure gravitation(U::UnitSystem,C::Coupling=universe(U)) = lightspeed(U,C)*planck(U,C)/normal(twopi(U))/planckmass(U,C)^2
@pure elementarycharge(U::UnitSystem,C::Coupling=universe(U)) = sqrt(two(U)*planck(U)/(vacuumpermeability(U)/finestructure(C))/(lightspeed(U)*rationalization(U)*lorentz(U)^2))

for unit ∈ (:coupling,:finestructure,:electronunit,:protonunit,:protonelectron,:darkenergydensity)
    @eval @pure $unit(U::UnitSystem) = $unit(universe(U))
end
for unit ∈ (:boltzmann,:planckreduced,:lightspeed,:vacuumpermeability,:permeability,:electronmass,:molarmass)
    @eval @pure $unit(U::UnitSystem,C::Coupling) = $unit(U)
end
for unit ∈ (Constants...,:angle,:vacuumpermeability)
    @eval @pure $unit(U::UnitSystem,S::UnitSystem) = isquantity(U,S) ? evaldim($unit)(U,S) : unit($unit(S)/$unit(U))
end
for unit ∈ (Convert...,:angle)
    @eval begin
        @pure @inline $unit(v::Real,U::UnitSystem) = isquantity(U) ? evaldim($unit)(v,U) : $unit(v,U,Metric)
        @pure @inline $unit(v::Real,U::UnitSystem,S::UnitSystem) = isquantity(U,S) ? evaldim($unit)(v,U,S) : (u=$unit(U,S);isone(u) ? v : v/u)
        @pure @inline $unit(v::Real,U::UnitSystem{kB,ħ,𝘤,μ₀,mₑ},S::UnitSystem{kB,ħ,𝘤,μ₀,mₑ}) where {kB,ħ,𝘤,μ₀,mₑ} = v
    end
    if unit ∉ (Constants...,:angle,:permeability)
        @eval @pure @inline $unit(U::UnitSystem) = isquantity(U) ? evaldim($unit)(U) : $unit(Natural,U)
    end
end
for unit ∈ (Systems...,Dimensionless...,Constants...,Physics...,Derived...,Convert...)
    @eval export $unit
end

# fundamental constants, αinv = (34259-1/4366.8123)/250 # 137.036 exactly?

const g₀,atm,T₀ = Constant(9.80665),Constant(101325.0),Constant(273.15)
const ft,ftUS,lb = Constant(0.3048),Constant(1200/3937),Constant(0.45359237)
const inHg,Ωᵢₜ,Vᵢₜ = Constant(1/3386.389),Constant(1.000495),Constant(1.00033)
const ΔνCs,Kcd,mP = Constant(9192631770.0),Constant(683*555.016/555),Constant(2.176434e-8)
const NA,kB,𝘩 = Constant(6.02214076e23),Constant(1.380649e-23),Constant(6.62607015e-34)
const 𝘤,𝘦 = Constant(299792458.),Constant(1.602176634e-19)
const μₑᵤ,μₚᵤ,μE☾ = Constant(1/1822.888486209),Constant(1.007276466621),Constant(81.300568)
const αinv,R∞ = Constant(137.035999084),Constant(10973731.5681601)
const RK1990,KJ1990,Rᵤ2014 = Constant(25812.807),Constant(4.835979e14),Constant(8.3144598)
const RK2014,KJ2014 = Constant(25812.8074555),Constant(4.835978525e14)
const GME,GMJ = Constant(398600441.8e6),Constant(1.26686534e17)
const kG,H0,ΩΛ = Constant(3548.18761),Constant(67.66),Constant(0.6889)
const aⱼ,au = Constant(365.25),Constant(149597870.7e3)
const LD,JD = Constant(384399e3),Constant(778479e6)
const zetta,zepto = Constant(1e21),Constant(1e-21)
const yotta,yocto = Constant(1e24),Constant(1e-24)
const 𝟏,𝟐,𝟑,𝟓,𝟕,𝟏𝟎,𝟏𝟏,𝟏𝟗,𝟒𝟑,τ,α = Constant(1),Constant(2),Constant(3),Constant(5),Constant(7),Constant(10),Constant(11),Constant(19),Constant(43),Constant(2π),inv(αinv)

include("initdata.jl")

const slug,lbm,lbmUS,rankine,kelvin = lb*g₀/ft,g₀/ft,g₀/ftUS,°R,K
const ħ1990,ħ2014 = planckreduced(Conventional),planckreduced(CODATA)
const mₑ1990,mₑ2014 = electronmass(Conventional),electronmass(CODATA)
const δμ₀,ly,mₛ,GG = μ₀-4π*1e-7,aⱼ*𝘤*DAY,GM☉/G,G
const units, temp = UnitSystem, temperature

# common conversion factors

const kcalₜₕ,kcal₄,kcal₁₀,kcal₂₀,kcalₘ,kcalᵢₜ = 4184,4204,4185.5,4182,4190,4186.8
const calₜₕ,cal₄,cal₁₀,cal₂₀,calₘ,calᵢₜ = (kcalₜₕ,kcal₄,kcal₁₀,kcal₂₀,kcalₘ,kcalᵢₜ)./1e3

@pure deka(U::UnitSystem) = two(U)*five(U)
@pure hecto(U::UnitSystem) = deka(U)^2
@pure kilo(U::UnitSystem) = deka(U)^3
@pure mega(U::UnitSystem) = (Constant(1.0)*kilo(U))^2
@pure giga(U::UnitSystem) = (Constant(1.0)*kilo(U))^3
@pure tera(U::UnitSystem) = (Constant(1.0)*kilo(U))^4
@pure peta(U::UnitSystem) = (Constant(1.0)*kilo(U))^5
@pure exa(U::UnitSystem) = (Constant(1.0)*kilo(U))^6
@pure zetta(U::UnitSystem) = (Constant(1.0)*kilo(U))^7
@pure yotta(U::UnitSystem) = (Constant(1.0)*kilo(U))^8
@pure deci(U::UnitSystem) = inv(deka(U))
@pure centi(U::UnitSystem) = inv(hecto(U))
@pure milli(U::UnitSystem) = inv(kilo(U))
@pure micro(U::UnitSystem) = inv(mega(U))
@pure nano(U::UnitSystem) = inv(giga(U))
@pure pico(U::UnitSystem) = inv(tera(U))
@pure femto(U::UnitSystem) = inv(peta(U))
@pure atto(U::UnitSystem) = inv(exa(U))
@pure zepto(U::UnitSystem) = inv(zetta(U))
@pure yocto(U::UnitSystem) = inv(yotta(U))
@pure kibi(U::UnitSystem) = two(U)^10
@pure mebi(U::UnitSystem) = two(U)^20
@pure gibi(U::UnitSystem) = two(U)^30
@pure tebi(U::UnitSystem) = two(U)^40
@pure pebi(U::UnitSystem) = two(U)^50
@pure exbi(U::UnitSystem) = two(U)^60
@pure zebi(U::UnitSystem) = (Constant(1.0)*two(U))^70
@pure yobi(U::UnitSystem) = (Constant(1.0)*two(U))^80

# physical constants

@pure electronmass(U::typeof(Planck),C::Coupling) = sqrt(spat(U)*coupling(C))
@pure electronmass(U::typeof(PlanckGauss),C::Coupling) = sqrt(coupling(C))
@pure electronmass(U::UnitSystem{kB,ħ,𝘤,μ₀,√(αG*αinv)},C::Coupling) where {kB,ħ,𝘤,μ₀} = sqrt(coupling(C)/finestructure(C))
@pure electronmass(U::UnitSystem{kB,ħ,𝘤,μ₀,1/μₚₑ},C::Coupling) where {kB,ħ,𝘤,μ₀} = 1/protonelectron(C)
@pure vacuumpermeability(U::UnitSystem{kB,ħ,𝘤,4π/αinv^2},C::Coupling) where {kB,ħ,𝘤} = spat(U)*finestructure(C)^2
@pure vacuumpermeability(U::UnitSystem{kB,ħ,𝘤,π/αinv^2},C::Coupling) where {kB,ħ,𝘤} = spat(U)/two(U)^2*finestructure(C)^2
@pure lightspeed(U::UnitSystem{kB,ħ,αinv},C::Coupling) where {kB,ħ} = inv(finestructure(C))
@pure lightspeed(U::UnitSystem{kB,ħ,2αinv},C::Coupling) where {kB,ħ} = two(U)/finestructure(C)
@pure planckreduced(U::UnitSystem{kB,αinv},C::Coupling) where kB = inv(finestructure(C))

@pure electronmass(U::UnitSystem{kB,ħ,𝘤,μ₀,mₑ},C::Coupling) where {kB,ħ,μ₀} = electronmass(U)#electronmass(planck(U),C)
@pure electronmass(U::UnitSystem{kB,ħ,100𝘤,μ₀,1000mₑ},C::Coupling) where {kB,ħ,μ₀} = electronmass(SI,C)*(two(U)*five(U))^3
@pure electronmass(U::UnitSystem{kB,ħ,𝘤,μ₀,mₑ/1000},C::Coupling) where {kB,ħ,μ₀} = electronmass(SI,C)/(two(U)*five(U))^3
@pure electronmass(U::UnitSystem{kB,ħ,𝘤,μ₀,electronmass(CODATA)},C::Coupling) where {kB,ħ,μ₀} = electronmass(planck(U),C)
@pure electronmass(U::UnitSystem{kB,ħ,𝘤,μ₀,electronmass(Conventional)},C::Coupling) where {kB,ħ,μ₀} = electronmass(planck(U),C)
@pure electronmass(U::UnitSystem{kB,ħ,𝘤/ftUS,μ₀,mₑ*ft/lb/g₀},C::Coupling) where {kB,ħ,μ₀} = electronmass(SI,C)*ft/lb/g₀
@pure vacuumpermeability(U::UnitSystem{kB,ħ,𝘤,μ₀},C::Coupling) where {kB,ħ,𝘤} = finestructure(C)*2𝘩/𝘤/𝘦^2
#@pure vacuumpermeability(U::typeof(ESU2019),C::Coupling) = 1e3*vacuumpermeability(SI,C)/𝘤^2
#@pure vacuumpermeability(U::typeof(EMU2019),C::Coupling) = 1e7*vacuumpermeability(SI,C)
@pure vacuumpermeability(U::typeof(CODATA),C::Coupling) = 2RK2014*finestructure(C)/𝘤
@pure vacuumpermeability(U::typeof(Conventional),C::Coupling) = 2RK1990*finestructure(C)/𝘤

convertext(unit,fun) = """
```Julia
$unit(U::UnitSystem,S::UnitSystem) = $fun
$unit(v::Real,U::UnitSystem,S::UnitSystem) = v/$unit(U,S)
```
"""

unitext(unit,text) = """
```Julia
$unit(U::UnitSystem) = $text
```
"""

include("kinematic.jl")
include("electromagnetic.jl")
include("thermodynamic.jl")
include("physics.jl")
include("systems.jl")

const RK = klitzing(SI2019) #
const KJ = josephson(SI2019) #

end # module


text:
module GMMbox

using Distances
using Distributions
using Clustering
import Base.rand
import Base.rand!

export gaussian, logprob!, logprob, gmm, rand, rand!

mutable struct gaussian{T<:AbstractFloat}
  ndim::Int64
  mean::AbstractArray{T,2}
  cov::AbstractArray{T,2}
end

gaussian(ndim::Integer) = gaussian(ndim,zeros(ndim,1),eye(ndim))
gaussian(data::Matrix{T} where T<:AbstractFloat)  = gaussian(size(data,1),mean(data,2),cov(data,2))

mutable struct gmm{T<:AbstractFloat}
  nmix::Int64
  wts::AbstractArray{T,2}
  mix::AbstractArray{gaussian{T}}
end

function gmm(nmix::Integer, data::Matrix{T} where T <: AbstractFloat; method = "rand")
  if(method=="rand")
    ndim = size(data,1);
    dcov = cov(data,2);
    dmvn = MvNormal(vec(mean(data,2)),dcov);
    mix = Array{gaussian{eltype(data)}}(nmix);
    for i=1:nmix
      mixmean = convert(Array{eltype(data)},rand(dmvn));
      mix[i] = gaussian(ndim,reshape(mixmean,ndim,1),dcov);
    end
  elseif(method=="kmeans")
    km = kmeans(data,nmix);
    mix = Array{gaussian{eltype(data)}}(nmix);
    dcov = cov(data,2);
    for i=1:nmix
      mix[i] = gaussian(data[:,find(km.assignments.==i)]);
      mix[i].cov = dcov;
    end
  end

  return gmm(nmix,convert(Array{eltype(data)},(1/nmix)*ones(nmix,1)),mix);
end

rand!(g::gaussian{T},data::Matrix{T}) where T <: AbstractFloat = rand!(MvNormal(vec(g.mean),g.cov),data);
rand(g::gaussian{T},nvec::Integer) where T <: AbstractFloat = rand(MvNormal(vec(g.mean),g.cov),nvec);

function rand!(gm::gmm{T},data::Matrix{T}) where T <: AbstractFloat
  ftype = eltype(data);
  nvec = size(data,2);
  pmix = rand(ftype,nvec);
  nvecmix = zeros(Int64,gm.nmix);
  thr = zeros(ftype,1);
  for i=1:gm.nmix
    nvecmix[i] = sum((pmix.>=thr).&(pmix.<thr+gm.wts[i]))
    thr += gm.wts[i];
  end
  if(sum(nvecmix)<nvec)
    nvecmix[end] += nvec-sum(nvecmix);
    println("sum of number of samples per mixture not equal to requested number of samples");
  end
  sindx = 1;
  eindx = 0;
  for i=1:gm.nmix
    eindx += nvecmix[i];
    data[:,sindx:eindx] = rand(gm.mix[i],nvecmix[i]);
    sindx += nvecmix[i];
  end
  return nothing
end

function rand(gm::gmm{T},nvec::Integer) where T <: AbstractFloat
  ftype = eltype(gm.wts);
  pmix = rand(ftype,nvec);
  nvecmix = zeros(Int64,gm.nmix);
  thr = zeros(ftype,1);
  for i=1:gm.nmix
    nvecmix[i] = sum((pmix.>=thr).&(pmix.<thr+gm.wts[i]))
    thr += gm.wts[i];
  end
  if(sum(nvecmix)<nvec)
    nvecmix[end] += nvec-sum(nvecmix);
    println("sum of number of samples per mixture not equal to requested number of samples");
  end
  data = zeros(ftype,gm.mix[1].ndim,sum(nvec));
  sindx = 1;
  eindx = 0;
  for i=1:gm.nmix
    eindx += nvecmix[i];
    data[:,sindx:eindx] = rand(gm.mix[i],nvecmix[i]);
    sindx += nvecmix[i];
  end
  return data
end

logprob!(p::Array{T},g::gaussian{T},data::Matrix{T}) where T <: AbstractFloat = logpdf!(p,MvNormal(vec(g.mean),g.cov),data)
logprob(g::gaussian{T},data::Matrix{T}) where T <: AbstractFloat = logpdf(MvNormal(vec(g.mean),g.cov),data)

function logprob!(p::Array{T},gm::gmm{T},data::Matrix{T}) where T <: AbstractFloat
  mixprob = zeros(eltype(data),size(data,2),gm.nmix);
  for i=1:gm.nmix
    mixprob[:,i] = logprob(gm.mix[i],data);
  end
  maxprob = maximum(mixprob,2);
  broadcast!(-,mixprob,mixprob,maxprob);
  p .= vec(log.((exp.(mixprob))*gm.wts) .+ maxprob);
  return nothing
end

function logprob(gm::gmm{T},data::Matrix{T}) where T <: AbstractFloat
  mixprob = zeros(eltype(data),size(data,2),gm.nmix);
  for i=1:gm.nmix
    mixprob[:,i] = logprob(gm.mix[i],data);
  end
  maxprob = maximum(mixprob,2);
  broadcast!(-,mixprob,mixprob,maxprob);
  return vec(log.((exp.(mixprob))*gm.wts) .+ maxprob);
end

#=
function logprob!(p::Matrix{T},g::gaussian{T},data::Matrix{T}) where T <: AbstractFloat
  logpdf!(reshape(p,length(p)),MvNormal(vec(g.mean),g.cov),data);
  reshape(p,length(p),1);
end
=#

#=
#logpdf and logpdf! using Distributions.jl is faster than this version which uses Distances.jl for Float32 - not using this version for CPU code
function logprob!(p::Matrix{T},g::gaussian{T},data::Matrix{T}) where T <: AbstractFloat
  ndim = size(data,1);
  icov = inv(g.cov);
  lconst = ndim*log(2*pi) + log(det(g.cov));
  pairwise!(p,SqMahalanobis(icov),data,g.mean);
  p .+= lconst;
  p .*= -0.5;
  return nothing
end

function logprob(g::gaussian{T},data::Matrix{T}) where T <: AbstractFloat
  p = zeros(eltype(data),size(data,2),1);
  ndim = size(data,1);
  icov = inv(g.cov);
  lconst = ndim*log(2*pi) + log(det(g.cov));
  pairwise!(p,SqMahalanobis(icov),data,g.mean);
  p .+= lconst;
  p .*= -0.5;
  return p
end
=#

end

ext: jl


text:
# https://plot.ly/javascript/getting-started

is_subplot_supported(::PlotlyBackend) = true
# is_string_supported(::PlotlyBackend) = true
function _plotly_framestyle(style::Symbol)
    if style in (:box, :axes, :zerolines, :grid, :none)
        return style
    else
        default_style = get((semi = :box, origin = :zerolines), style, :axes)
        @warn("Framestyle :$style is not supported by Plotly and PlotlyJS. :$default_style was cosen instead.")
        default_style
    end
end

# --------------------------------------------------------------------------------------

using UUIDs

# ----------------------------------------------------------------

plotly_legend_pos(pos::Symbol) = get(
    (
        right = [1.0, 0.5],
        left = [0.0, 0.5],
        top = [0.5, 1.0],
        bottom = [0.5, 0.0],
        bottomleft = [0.0, 0.0],
        bottomright = [1.0, 0.0],
        topright = [1.0, 1.0],
        topleft = [0.0, 1.0],
    ),
    pos,
    [1.0, 1.0],
)
plotly_legend_pos(v::Tuple{S,T}) where {S<:Real, T<:Real} = v

function plotly_font(font::Font, color = font.color)
    KW(
        :family => font.family,
        :size   => round(Int, font.pointsize*1.4),
        :color  => rgba_string(color),
    )
end

function plotly_annotation_dict(x, y, val; xref="paper", yref="paper")
    KW(
        :text => val,
        :xref => xref,
        :x => x,
        :yref => yref,
        :y => y,
        :showarrow => false,
    )
end

function plotly_annotation_dict(x, y, ptxt::PlotText; xref="paper", yref="paper")
    merge(plotly_annotation_dict(x, y, ptxt.str; xref=xref, yref=yref), KW(
        :font => plotly_font(ptxt.font),
        :xanchor => ptxt.font.halign == :hcenter ? :center : ptxt.font.halign,
        :yanchor => ptxt.font.valign == :vcenter ? :middle : ptxt.font.valign,
        :rotation => -ptxt.font.rotation,
    ))
end

# function get_annotation_dict_for_arrow(plotattributes::KW, xyprev::Tuple, xy::Tuple, a::Arrow)
#     xdiff = xyprev[1] - xy[1]
#     ydiff = xyprev[2] - xy[2]
#     dist = sqrt(xdiff^2 + ydiff^2)
#     KW(
#         :showarrow => true,
#         :x => xy[1],
#         :y => xy[2],
#         # :ax => xyprev[1] - xy[1],
#         # :ay => xy[2] - xyprev[2],
#         # :ax => 0,
#         # :ay => -40,
#         :ax => 10xdiff / dist,
#         :ay => -10ydiff / dist,
#         :arrowcolor => rgba_string(plotattributes[:linecolor]),
#         :xref => "x",
#         :yref => "y",
#         :arrowsize => 10a.headwidth,
#         # :arrowwidth => a.headlength,
#         :arrowwidth => 0.1,
#     )
# end

function plotly_scale(scale::Symbol)
    if scale == :log10
        "log"
    else
        "-"
    end
end

function shrink_by(lo, sz, ratio)
    amt = 0.5 * (1.0 - ratio) * sz
    lo + amt, sz - 2amt
end

function plotly_apply_aspect_ratio(sp::Subplot, plotarea, pcts)
    aspect_ratio = get_aspect_ratio(sp)
    if aspect_ratio != :none
        if aspect_ratio == :equal
            aspect_ratio = 1.0
        end
        xmin,xmax = axis_limits(sp, :x)
        ymin,ymax = axis_limits(sp, :y)
        want_ratio = ((xmax-xmin) / (ymax-ymin)) / aspect_ratio
        parea_ratio = width(plotarea) / height(plotarea)
        if want_ratio > parea_ratio
            # need to shrink y
            ratio = parea_ratio / want_ratio
            pcts[2], pcts[4] = shrink_by(pcts[2], pcts[4], ratio)
        elseif want_ratio < parea_ratio
            # need to shrink x
            ratio = want_ratio / parea_ratio
            pcts[1], pcts[3] = shrink_by(pcts[1], pcts[3], ratio)
        end
        pcts
    end
    pcts
end

# this method gets the start/end in percentage of the canvas for this axis direction
function plotly_domain(sp::Subplot, letter)
    figw, figh = sp.plt[:size]
    pcts = bbox_to_pcts(sp.plotarea, figw*px, figh*px)
    pcts = plotly_apply_aspect_ratio(sp, sp.plotarea, pcts)
    i1,i2 = (letter == :x ? (1,3) : (2,4))
    [pcts[i1], pcts[i1]+pcts[i2]]
end

function plotly_axis(plt::Plot, axis::Axis, sp::Subplot)
    letter = axis[:letter]
    framestyle = sp[:framestyle]
    ax = KW(
        :visible    => framestyle != :none,
        :title      => axis[:guide],
        :showgrid   => axis[:grid],
        :gridcolor  => rgba_string(plot_color(axis[:foreground_color_grid], axis[:gridalpha])),
        :gridwidth  => axis[:gridlinewidth],
        :zeroline   => framestyle == :zerolines,
        :zerolinecolor => rgba_string(axis[:foreground_color_axis]),
        :showline   => framestyle in (:box, :axes) && axis[:showaxis],
        :linecolor  => rgba_string(plot_color(axis[:foreground_color_axis])),
        :ticks      => axis[:tick_direction] == :out ? "outside" : "inside",
        :mirror     => framestyle == :box,
        :showticklabels => axis[:showaxis],
    )

    if letter in (:x,:y)
        ax[:domain] = plotly_domain(sp, letter)
        if RecipesPipeline.is3d(sp)
            # don't link 3d axes for synchronized interactivity
            x_idx = y_idx = sp[:subplot_index]
        else
            x_idx, y_idx = plotly_link_indicies(plt, sp)
        end
        ax[:anchor] = "$(letter==:x ? "y$(y_idx)" : "x$(x_idx)")"
    end

    ax[:tickangle] = -axis[:rotation]
    ax[:type] = plotly_scale(axis[:scale])
    lims = axis_limits(sp, letter)

    if axis[:ticks] != :native || axis[:lims] != :auto
        ax[:range] = map(RecipesPipeline.scale_func(axis[:scale]), lims)
    end

    if !(axis[:ticks] in (nothing, :none, false))
        ax[:titlefont] = plotly_font(guidefont(axis))
        ax[:tickfont] = plotly_font(tickfont(axis))
        ax[:tickcolor] = framestyle in (:zerolines, :grid) || !axis[:showaxis] ? rgba_string(invisible()) : rgb_string(axis[:foreground_color_axis])
        ax[:linecolor] = rgba_string(axis[:foreground_color_axis])

        # ticks
        if axis[:ticks] != :native
            ticks = get_ticks(sp, axis)
            ttype = ticksType(ticks)
            if ttype == :ticks
                ax[:tickmode] = "array"
                ax[:tickvals] = ticks
            elseif ttype == :ticks_and_labels
                ax[:tickmode] = "array"
                ax[:tickvals], ax[:ticktext] = ticks
            end
        end
    else
        ax[:showticklabels] = false
        ax[:showgrid] = false
    end

    # flip
    if axis[:flip]
        ax[:range] = reverse(ax[:range])
    end

    ax
end

function plotly_polaraxis(sp::Subplot, axis::Axis)
    ax = KW(
        :visible => axis[:showaxis],
        :showline => axis[:grid],
    )

    if axis[:letter] == :x
        ax[:range] = rad2deg.(axis_limits(sp, :x))
    else
        ax[:range] = axis_limits(sp, :y)
        ax[:orientation] = -90
    end

    ax
end

function plotly_layout(plt::Plot)
    plotattributes_out = KW()

    w, h = plt[:size]
    plotattributes_out[:width], plotattributes_out[:height] = w, h
    plotattributes_out[:paper_bgcolor] = rgba_string(plt[:background_color_outside])
    plotattributes_out[:margin] = KW(:l=>0, :b=>20, :r=>0, :t=>20)

    plotattributes_out[:annotations] = KW[]

    multiple_subplots = length(plt.subplots) > 1

    for sp in plt.subplots
        spidx = multiple_subplots ? sp[:subplot_index] : ""
        x_idx, y_idx = multiple_subplots ? plotly_link_indicies(plt, sp) : ("", "")
        # add an annotation for the title... positioned horizontally relative to plotarea,
        # but vertically just below the top of the subplot bounding box
        if sp[:title] != ""
            bb = plotarea(sp)
            tpos = sp[:titlelocation]
            xmm = if tpos == :left
                left(bb)
            elseif tpos == :right
                right(bb)
            else
                0.5 * (left(bb) + right(bb))
            end
            titlex, titley = xy_mm_to_pcts(xmm, top(bbox(sp)), w*px, h*px)
            title_font = font(titlefont(sp), :top)
            push!(plotattributes_out[:annotations], plotly_annotation_dict(titlex, titley, text(sp[:title], title_font)))
        end

        plotattributes_out[:plot_bgcolor] = rgba_string(sp[:background_color_inside])

        # set to supported framestyle
        sp[:framestyle] = _plotly_framestyle(sp[:framestyle])

        # if any(RecipesPipeline.is3d, seriesargs)
        if RecipesPipeline.is3d(sp)
            azim = sp[:camera][1] - 90 #convert azimuthal to match GR behaviour
            theta = 90 - sp[:camera][2] #spherical coordinate angle from z axis
            plotattributes_out[:scene] = KW(
                Symbol("xaxis$(spidx)") => plotly_axis(plt, sp[:xaxis], sp),
                Symbol("yaxis$(spidx)") => plotly_axis(plt, sp[:yaxis], sp),
                Symbol("zaxis$(spidx)") => plotly_axis(plt, sp[:zaxis], sp),

                #2.6 multiplier set camera eye such that whole plot can be seen
                :camera => KW(
                    :eye => KW(
                        :x => cosd(azim)*sind(theta)*2.6,
                        :y => sind(azim)*sind(theta)*2.6,
                        :z => cosd(theta)*2.6,
                    ),
                ),
            )
        elseif ispolar(sp)
            plotattributes_out[Symbol("angularaxis$(spidx)")] = plotly_polaraxis(sp, sp[:xaxis])
            plotattributes_out[Symbol("radialaxis$(spidx)")] = plotly_polaraxis(sp, sp[:yaxis])
        else
            plotattributes_out[Symbol("xaxis$(x_idx)")] = plotly_axis(plt, sp[:xaxis], sp)
            # don't allow yaxis to be reupdated/reanchored in a linked subplot
            spidx == y_idx ? plotattributes_out[Symbol("yaxis$(y_idx)")] = plotly_axis(plt, sp[:yaxis], sp) : nothing
        end

        # legend
        plotattributes_out[:showlegend] = sp[:legend] != :none
        xpos,ypos = plotly_legend_pos(sp[:legend])
        if sp[:legend] != :none
            plotattributes_out[:legend] = KW(
                :bgcolor  => rgba_string(sp[:background_color_legend]),
                :bordercolor => rgba_string(sp[:foreground_color_legend]),
                :font     => plotly_font(legendfont(sp)),
                :tracegroupgap => 0,
                :x => xpos,
                :y => ypos
            )
        end

        # annotations
        for ann in sp[:annotations]
            append!(plotattributes_out[:annotations], KW[plotly_annotation_dict(locate_annotation(sp, ann...)...; xref = "x$(x_idx)", yref = "y$(y_idx)")])
        end
        # series_annotations
        for series in series_list(sp)
            anns = series[:series_annotations]
            for (xi,yi,str,fnt) in EachAnn(anns, series[:x], series[:y])
                push!(plotattributes_out[:annotations], plotly_annotation_dict(
                    xi,
                    yi,
                    PlotText(str,fnt); xref = "x$(x_idx)", yref = "y$(y_idx)")
                )
            end
        end

        # # arrows
        # for sargs in seriesargs
        #     a = sargs[:arrow]
        #     if sargs[:seriestype] in (:path, :line) && typeof(a) <: Arrow
        #         add_arrows(sargs[:x], sargs[:y]) do xyprev, xy
        #             push!(plotattributes_out[:annotations], get_annotation_dict_for_arrow(sargs, xyprev, xy, a))
        #         end
        #     end
        # end

        if ispolar(sp)
            plotattributes_out[:direction] = "counterclockwise"
        end

        plotattributes_out
    end

    # turn off hover if nothing's using it
    if all(series -> series.plotattributes[:hover] in (false,:none), plt.series_list)
        plotattributes_out[:hovermode] = "none"
    end

    plotattributes_out = recursive_merge(plotattributes_out, plt.attr[:extra_plot_kwargs])
end

function plotly_layout_json(plt::Plot)
    JSON.json(plotly_layout(plt), 4)
end

plotly_colorscale(cg::ColorGradient, α = nothing) =
    [[v, rgba_string(plot_color(cg.colors[v], α))] for v in cg.values]
function plotly_colorscale(c::AbstractVector{<:Colorant}, α = nothing)
    if length(c) == 1
        return [
            [0.0, rgba_string(plot_color(c[1], α))],
            [1.0, rgba_string(plot_color(c[1], α))],
        ]
    else
        vals = range(0.0, stop = 1.0, length = length(c))
        return [[vals[i], rgba_string(plot_color(c[i], α))] for i in eachindex(c)]
    end
end
function plotly_colorscale(cg::PlotUtils.CategoricalColorGradient, α = nothing)
    n = length(cg)
    cinds = repeat(1:n, inner = 2)
    vinds = vcat((i:(i + 1) for i in 1:n)...)
    return [
        [cg.values[vinds[i]], rgba_string(plot_color(color_list(cg)[cinds[i]], α))]
        for i in eachindex(cinds)
    ]
end
plotly_colorscale(c, α = nothing) = plotly_colorscale(_as_gradient(c), α)

get_plotly_marker(k, def) = get(
    (
        rect = "square",
        xcross = "x",
        x = "x",
        utriangle = "triangle-up",
        dtriangle = "triangle-down",
        star5 = "star-triangle-up",
        vline = "line-ns",
        hline = "line-ew",
    ),
    k,
    def,
)

# find indicies of axes to which the supblot links to
function plotly_link_indicies(plt::Plot, sp::Subplot)
    if plt[:link] in (:x, :y, :both)
        x_idx = sp[:xaxis].sps[1][:subplot_index]
        y_idx = sp[:yaxis].sps[1][:subplot_index]
    else
        x_idx = y_idx = sp[:subplot_index]
    end
    x_idx, y_idx
end

# the Shape contructor will automatically close the shape. since we need it closed,
# we split by NaNs and then construct/destruct the shapes to get the closed coords
function plotly_close_shapes(x, y)
    xs, ys = nansplit(x), nansplit(y)
    for i=eachindex(xs)
        shape = Shape(xs[i], ys[i])
        xs[i], ys[i] = coords(shape)
    end
    nanvcat(xs), nanvcat(ys)
end

function plotly_data(series::Series, letter::Symbol, data)
    axis = series[:subplot][Symbol(letter, :axis)]

    data = if axis[:ticks] == :native && data !== nothing
        plotly_native_data(axis, data)
    else
       data
    end

    if series[:seriestype] in (:heatmap, :contour, :surface, :wireframe)
        plotly_surface_data(series, data)
    else
        plotly_data(data)
    end
end
plotly_data(v) = v !== nothing ? collect(v) : v
plotly_data(v::AbstractArray) = v
plotly_data(surf::Surface) = surf.surf
plotly_data(v::AbstractArray{R}) where {R<:Rational} = float(v)

plotly_surface_data(series::Series, a::AbstractVector) = a
plotly_surface_data(series::Series, a::AbstractMatrix) = transpose_z(series, a, false)
plotly_surface_data(series::Series, a::Surface) = plotly_surface_data(series, a.surf)

function plotly_native_data(axis::Axis, data::AbstractArray)
    if !isempty(axis[:discrete_values])
        construct_categorical_data(data, axis)
    elseif axis[:formatter] in (datetimeformatter, dateformatter, timeformatter)
        plotly_convert_to_datetime(data, axis[:formatter])
    else
        data
    end
end
plotly_native_data(axis::Axis, a::Surface) = Surface(plotly_native_data(axis, a.surf))

function plotly_convert_to_datetime(x::AbstractArray, formatter::Function)
    if formatter == datetimeformatter
        map(xi -> replace(formatter(xi), "T" => " "), x)
    elseif formatter == dateformatter
        map(xi -> string(formatter(xi), " 00:00:00"), x)
    elseif formatter == timeformatter
        map(xi -> string(Dates.Date(Dates.now()), " ", formatter(xi)), x)
    else
        error("Invalid DateTime formatter. Expected Plots.datetime/date/time formatter but got $formatter")
    end
end
#ensures that a gradient is called if a single color is supplied where a gradient is needed (e.g. if a series recipe defines marker_z)
as_gradient(grad::ColorGradient, α) = grad
as_gradient(grad, α) = cgrad(alpha = α)

# get a dictionary representing the series params (plotattributes is the Plots-dict, plotattributes_out is the Plotly-dict)
function plotly_series(plt::Plot, series::Series)
    st = series[:seriestype]

    sp = series[:subplot]
    clims = get_clims(sp, series)

    if st == :shape
        return plotly_series_shapes(plt, series, clims)
    end

    plotattributes_out = KW()

    # these are the axes that the series should be mapped to
    x_idx, y_idx = plotly_link_indicies(plt, sp)
    plotattributes_out[:xaxis] = "x$(x_idx)"
    plotattributes_out[:yaxis] = "y$(y_idx)"
    plotattributes_out[:showlegend] = should_add_to_legend(series)

    if st == :straightline
        x, y = straightline_data(series, 100)
        z = series[:z]
    else
        x, y, z  = series[:x], series[:y], series[:z]
    end

    x, y, z = (plotly_data(series, letter, data)
        for (letter, data) in zip((:x, :y, :z), (x, y, z))
    )

    plotattributes_out[:name] = series[:label]

    isscatter = st in (:scatter, :scatter3d, :scattergl)
    hasmarker = isscatter || series[:markershape] != :none
    hasline = st in (:path, :path3d, :straightline)
    hasfillrange = st in (:path, :scatter, :scattergl, :straightline) &&
        (isa(series[:fillrange], AbstractVector) || isa(series[:fillrange], Tuple))

    plotattributes_out[:colorbar] = KW(:title => sp[:colorbar_title])

    if is_2tuple(clims)
        plotattributes_out[:zmin], plotattributes_out[:zmax] = clims
    end

    # set the "type"
    if st in (:path, :scatter, :scattergl, :straightline, :path3d, :scatter3d)
        return plotly_series_segments(series, plotattributes_out, x, y, z, clims)

    elseif st == :heatmap
        x = heatmap_edges(x, sp[:xaxis][:scale])
        y = heatmap_edges(y, sp[:yaxis][:scale])
        plotattributes_out[:type] = "heatmap"
        plotattributes_out[:x], plotattributes_out[:y], plotattributes_out[:z] = x, y, z
        plotattributes_out[:colorscale] = plotly_colorscale(series[:fillcolor], series[:fillalpha])
        plotattributes_out[:showscale] = hascolorbar(sp)

    elseif st == :contour
        filled = isfilledcontour(series)
        plotattributes_out[:type] = "contour"
        plotattributes_out[:x], plotattributes_out[:y], plotattributes_out[:z] = x, y, z
        plotattributes_out[:ncontours] = series[:levels] + 2
        plotattributes_out[:contours] = KW(:coloring => filled ? "fill" : "lines", :showlabels => series[:contour_labels] == true)
        plotattributes_out[:colorscale] = plotly_colorscale(series[:linecolor], series[:linealpha])
        plotattributes_out[:showscale] = hascolorbar(sp) && hascolorbar(series)

    elseif st in (:surface, :wireframe)
        plotattributes_out[:type] = "surface"
        plotattributes_out[:x], plotattributes_out[:y], plotattributes_out[:z] = x, y, z
        if st == :wireframe
            plotattributes_out[:hidesurface] = true
            wirelines = KW(
                :show => true,
                :color => rgba_string(plot_color(series[:linecolor], series[:linealpha])),
                :highlightwidth => series[:linewidth],
            )
            plotattributes_out[:contours] = KW(:x => wirelines, :y => wirelines, :z => wirelines)
            plotattributes_out[:showscale] = false
        else
            plotattributes_out[:colorscale] = plotly_colorscale(series[:fillcolor], series[:fillalpha])
            plotattributes_out[:opacity] = series[:fillalpha]
            if series[:fill_z] !== nothing
                plotattributes_out[:surfacecolor] = plotly_surface_data(series, series[:fill_z])
            end
            plotattributes_out[:showscale] = hascolorbar(sp)
        end

    else
        @warn("Plotly: seriestype $st isn't supported.")
        return KW()
    end

    # add "marker"
    if hasmarker
        inds = eachindex(x)
        plotattributes_out[:marker] = KW(
            :symbol => get_plotly_marker(series[:markershape], string(series[:markershape])),
            # :opacity => series[:markeralpha],
            :size => 2 * _cycle(series[:markersize], inds),
            :color => rgba_string.(plot_color.(get_markercolor.(series, inds), get_markeralpha.(series, inds))),
            :line => KW(
                :color => rgba_string.(plot_color.(get_markerstrokecolor.(series, inds), get_markerstrokealpha.(series, inds))),
                :width => _cycle(series[:markerstrokewidth], inds),
            ),
        )
    end

    plotly_polar!(plotattributes_out, series)
    plotly_hover!(plotattributes_out, series[:hover])

    return [plotattributes_out]
end

function plotly_series_shapes(plt::Plot, series::Series, clims)
    segments = iter_segments(series)
    plotattributes_outs = Vector{KW}(undef, length(segments))

    # TODO: create a plotattributes_out for each polygon
    # x, y = series[:x], series[:y]

    # these are the axes that the series should be mapped to
    x_idx, y_idx = plotly_link_indicies(plt, series[:subplot])
    plotattributes_base = KW(
        :xaxis => "x$(x_idx)",
        :yaxis => "y$(y_idx)",
        :name => series[:label],
        :legendgroup => series[:label],
    )

    x, y = (plotly_data(series, letter, data)
        for (letter, data) in zip((:x, :y), shape_data(series, 100))
    )

    for (i,rng) in enumerate(segments)
        length(rng) < 2 && continue

        # to draw polygons, we actually draw lines with fill
        plotattributes_out = merge(plotattributes_base, KW(
            :type => "scatter",
            :mode => "lines",
            :x => vcat(x[rng], x[rng[1]]),
            :y => vcat(y[rng], y[rng[1]]),
            :fill => "tozeroy",
            :fillcolor => rgba_string(plot_color(get_fillcolor(series, clims, i), get_fillalpha(series, i))),
        ))
        if series[:markerstrokewidth] > 0
            plotattributes_out[:line] = KW(
                :color => rgba_string(plot_color(get_linecolor(series, clims, i), get_linealpha(series, i))),
                :width => get_linewidth(series, i),
                :dash => string(get_linestyle(series, i)),
            )
        end
        plotattributes_out[:showlegend] = i==1 ? should_add_to_legend(series) : false
        plotly_polar!(plotattributes_out, series)
        plotly_hover!(plotattributes_out, _cycle(series[:hover], i))
        plotattributes_outs[i] = plotattributes_out
    end
    if series[:fill_z] !== nothing
        push!(plotattributes_outs, plotly_colorbar_hack(series, plotattributes_base, :fill))
    elseif series[:line_z] !== nothing
        push!(plotattributes_outs, plotly_colorbar_hack(series, plotattributes_base, :line))
    elseif series[:marker_z] !== nothing
        push!(plotattributes_outs, plotly_colorbar_hack(series, plotattributes_base, :marker))
    end
    plotattributes_outs
end

function plotly_series_segments(series::Series, plotattributes_base::KW, x, y, z, clims)
    st = series[:seriestype]
    sp = series[:subplot]
    isscatter = st in (:scatter, :scatter3d, :scattergl)
    hasmarker = isscatter || series[:markershape] != :none
    hasline = st in (:path, :path3d, :straightline)
    hasfillrange = st in (:path, :scatter, :scattergl, :straightline) &&
        (isa(series[:fillrange], AbstractVector) || isa(series[:fillrange], Tuple))

    segments = iter_segments(series)
    plotattributes_outs = fill(KW(), (hasfillrange ? 2 : 1 ) * length(segments))

    for (i,rng) in enumerate(segments)
        plotattributes_out = deepcopy(plotattributes_base)
        plotattributes_out[:showlegend] = i==1 ? should_add_to_legend(series) : false
        plotattributes_out[:legendgroup] = series[:label]

        # set the type
        if st in (:path, :scatter, :scattergl, :straightline)
            plotattributes_out[:type] = st==:scattergl ? "scattergl" : "scatter"
            plotattributes_out[:mode] = if hasmarker
                hasline ? "lines+markers" : "markers"
            else
                hasline ? "lines" : "none"
            end
            if series[:fillrange] == true || series[:fillrange] == 0 || isa(series[:fillrange], Tuple)
                plotattributes_out[:fill] = "tozeroy"
                plotattributes_out[:fillcolor] = rgba_string(plot_color(get_fillcolor(series, clims, i), get_fillalpha(series, i)))
            elseif typeof(series[:fillrange]) <: Union{AbstractVector{<:Real}, Real}
                plotattributes_out[:fill] = "tonexty"
                plotattributes_out[:fillcolor] = rgba_string(plot_color(get_fillcolor(series, clims, i), get_fillalpha(series, i)))
            elseif !(series[:fillrange] in (false, nothing))
                @warn("fillrange ignored... plotly only supports filling to zero and to a vector of values. fillrange: $(series[:fillrange])")
            end
            plotattributes_out[:x], plotattributes_out[:y] = x[rng], y[rng]

        elseif st in (:path3d, :scatter3d)
            plotattributes_out[:type] = "scatter3d"
            plotattributes_out[:mode] = if hasmarker
                hasline ? "lines+markers" : "markers"
            else
                hasline ? "lines" : "none"
            end
            plotattributes_out[:x], plotattributes_out[:y], plotattributes_out[:z] = x[rng], y[rng], z[rng]
        end

        # add "marker"
        if hasmarker
            plotattributes_out[:marker] = KW(
                :symbol => get_plotly_marker(_cycle(series[:markershape], i), string(_cycle(series[:markershape], i))),
                # :opacity => series[:markeralpha],
                :size => 2 * _cycle(series[:markersize], i),
                :color => rgba_string(plot_color(get_markercolor(series, clims, i), get_markeralpha(series, i))),
                :line => KW(
                    :color => rgba_string(plot_color(get_markerstrokecolor(series, i), get_markerstrokealpha(series, i))),
                    :width => _cycle(series[:markerstrokewidth], i),
                ),
            )
        end

        # add "line"
        if hasline
            plotattributes_out[:line] = KW(
                :color => rgba_string(plot_color(get_linecolor(series, clims, i), get_linealpha(series, i))),
                :width => get_linewidth(series, i),
                :shape => if st == :steppre
                    "vh"
                elseif st == :steppost
                    "hv"
                else
                    "linear"
                end,
                :dash => string(get_linestyle(series, i)),
            )
        end

        plotly_polar!(plotattributes_out, series)
        plotly_hover!(plotattributes_out, _cycle(series[:hover], rng))

        if hasfillrange
            # if hasfillrange is true, return two dictionaries (one for original
            # series, one for series being filled to) instead of one
            plotattributes_out_fillrange = deepcopy(plotattributes_out)
            plotattributes_out_fillrange[:showlegend] = false
            # if fillrange is provided as real or tuple of real, expand to array
            if typeof(series[:fillrange]) <: Real
                series[:fillrange] = fill(series[:fillrange], length(rng))
            elseif typeof(series[:fillrange]) <: Tuple
                f1 = typeof(series[:fillrange][1]) <: Real ? fill(series[:fillrange][1], length(rng)) : series[:fillrange][1][rng]
                f2 = typeof(series[:fillrange][2]) <: Real ? fill(series[:fillrange][2], length(rng)) : series[:fillrange][2][rng]
                series[:fillrange] = (f1, f2)
            end
            if isa(series[:fillrange], AbstractVector)
                plotattributes_out_fillrange[:y] = series[:fillrange][rng]
                delete!(plotattributes_out_fillrange, :fill)
                delete!(plotattributes_out_fillrange, :fillcolor)
            else
                # if fillrange is a tuple with upper and lower limit, plotattributes_out_fillrange
                # is the series that will do the filling
                fillrng = Tuple(series[:fillrange][i][rng] for i in 1:2)
                plotattributes_out_fillrange[:x], plotattributes_out_fillrange[:y] = concatenate_fillrange(x[rng], fillrng)
                plotattributes_out_fillrange[:line][:width] = 0
                delete!(plotattributes_out, :fill)
                delete!(plotattributes_out, :fillcolor)
            end

            plotattributes_outs[(2 * i - 1):(2 * i)] = [plotattributes_out_fillrange, plotattributes_out]
        else
            plotattributes_outs[i] = plotattributes_out
        end
    end

    if series[:line_z] !== nothing
        push!(plotattributes_outs, plotly_colorbar_hack(series, plotattributes_base, :line))
    elseif series[:fill_z] !== nothing
        push!(plotattributes_outs, plotly_colorbar_hack(series, plotattributes_base, :fill))
    elseif series[:marker_z] !== nothing
        push!(plotattributes_outs, plotly_colorbar_hack(series, plotattributes_base, :marker))
    end

    plotattributes_outs
end

function plotly_colorbar_hack(series::Series, plotattributes_base::KW, sym::Symbol)
    plotattributes_out = deepcopy(plotattributes_base)
    cmin, cmax = get_clims(series[:subplot])
    plotattributes_out[:showlegend] = false
    plotattributes_out[:type] = RecipesPipeline.is3d(series) ? :scatter3d : :scatter
    plotattributes_out[:hoverinfo] = :none
    plotattributes_out[:mode] = :markers
    plotattributes_out[:x], plotattributes_out[:y] = [series[:x][1]], [series[:y][1]]
    if RecipesPipeline.is3d(series)
        plotattributes_out[:z] = [series[:z][1]]
    end
    # zrange = zmax == zmin ? 1 : zmax - zmin # if all marker_z values are the same, plot all markers same color (avoids division by zero in next line)
    plotattributes_out[:marker] = KW(
        :size => 0,
        :opacity => 0,
        :color => [0.5],
        :cmin => cmin,
        :cmax => cmax,
        :colorscale => plotly_colorscale(series[Symbol("$(sym)color")], 1),
        :showscale => hascolorbar(series[:subplot]),
    )
    return plotattributes_out
end

function plotly_polar!(plotattributes_out::KW, series::Series)
    # convert polar plots x/y to theta/radius
    if ispolar(series[:subplot])
        theta, r = pop!(plotattributes_out, :x), pop!(plotattributes_out, :y)
        plotattributes_out[:t] = rad2deg.(theta)
        plotattributes_out[:r] = r
    end
end

function plotly_hover!(plotattributes_out::KW, hover)
    # hover text
    if hover === nothing || all(in([:none, false]), hover)
        plotattributes_out[:hoverinfo] = "none"
    elseif any(!isnothing, hover)
        plotattributes_out[:hoverinfo] = "text"
        plotattributes_out[:text] = hover
    end
end

# get a list of dictionaries, each representing the series params
function plotly_series(plt::Plot)
    slist = []
    for series in plt.series_list
        append!(slist, plotly_series(plt, series))
    end
    slist
end

# get json string for a list of dictionaries, each representing the series params
plotly_series_json(plt::Plot) = JSON.json(plotly_series(plt), 4)

# ----------------------------------------------------------------

html_head(plt::Plot{PlotlyBackend}) = plotly_html_head(plt)
html_body(plt::Plot{PlotlyBackend}) = plotly_html_body(plt)

const ijulia_initialized = Ref(false)

function plotly_html_head(plt::Plot)
    local_file = ("file:///" * plotly_local_file_path)
    plotly =
        use_local_dependencies[] ? local_file : "https://cdn.plot.ly/plotly-latest.min.js"

    include_mathjax = get(plt[:extra_plot_kwargs], :include_mathjax, "")
    mathjax_file = include_mathjax != "cdn" ? ("file://" * include_mathjax) : "https://cdnjs.cloudflare.com/ajax/libs/mathjax/2.7.7/MathJax.js?config=TeX-MML-AM_CHTML"
    mathjax_head = include_mathjax == "" ? "" : "<script src=\"$mathjax_file\"></script>\n\t\t"

    if isijulia() && !ijulia_initialized[]
        # using requirejs seems to be key to load a js depency in IJulia!
        # https://requirejs.org/docs/start.html
        # https://github.com/JuliaLang/IJulia.jl/issues/345
        display("text/html", """
            <script type="text/javascript">
                requirejs([$(repr(plotly))], function(p) {
                    window.Plotly = p
                });
            </script>
        """)
        ijulia_initialized[] = true
    end
    return "$mathjax_head<script src=$(repr(plotly))></script>"
end

function plotly_html_body(plt, style = nothing)
    if style === nothing
        w, h = plt[:size]
        style = "width:$(w)px;height:$(h)px;"
    end
    uuid = UUIDs.uuid4()
    html = """
        <div id=\"$(uuid)\" style=\"$(style)\"></div>
        <script>
        PLOT = document.getElementById('$(uuid)');
        Plotly.plot(PLOT, $(plotly_series_json(plt)), $(plotly_layout_json(plt)));
        </script>
    """
    html
end

function js_body(plt::Plot, uuid)
    js = """
        PLOT = document.getElementById('$(uuid)');
        Plotly.plot(PLOT, $(plotly_series_json(plt)), $(plotly_layout_json(plt)));
    """
end

function plotly_show_js(io::IO, plot::Plot)
    data = []
    for series in plot.series_list
        append!(data, plotly_series(plot, series))
    end
    layout = plotly_layout(plot)
    JSON.print(io, Dict(:data => data, :layout => layout))
end

# ----------------------------------------------------------------

Base.showable(::MIME"application/prs.juno.plotpane+html", plt::Plot{PlotlyBackend}) = true

function _show(io::IO, ::MIME"application/vnd.plotly.v1+json", plot::Plot{PlotlyBackend})
    plotly_show_js(io, plot)
end

function _show(io::IO, ::MIME"text/html", plt::Plot{PlotlyBackend})
    write(io, standalone_html(plt))
end

function _display(plt::Plot{PlotlyBackend})
    standalone_html_window(plt)
end


text:
using NKLandscapes
using FactCheck

srand(1)

facts("NKLandscapes.jl") do
  landscapes = [
    NKLandscape(10, 0),
    NKLandscape(10, 1),
    NKLandscape(10, 2),
    NKLandscape(10, 3),
    NKLandscape(10, 4),
    NKLandscape(10, 5),
    NKLandscape(10, 6),
    NKLandscape(10, 7),
    NKLandscape(10, 8),
    NKLandscape(10, 9),
    NKLandscape(10, 4; near=true),
    NKqLandscape(10, 1, 2),
    NKpLandscape(10, 1, 0.90),
  ]

  for l = landscapes
    context("$(l)") do
      context("Enumeration") do
        g = rand(Genotype, l)
        f = fitness(g, l)

        context("gtoi and itog") do
          ig = gtoi(g, l)
          gg = itog(ig, l)
          @fact gg --> g
        end

        context("fitness") do
          ig = gtoi(g, l)
          @fact fitness(ig, l) --> fitness(g, l)
        end
      end

      context("Landscapes") do
        g = rand(Genotype, l)
        f = fitness(g, l)

        context("Neighbors should differ at one locus") do
          function test_neighbors(genotype, landscape)
            nbrs = all_neighbors(genotype, landscape)
            for i = 1:number_neighbors(genotype, landscape)
              nbr = nbrs[:, i]
              @fact (genotype - nbr) |> sum |> abs --> 1
            end
          end

          test_neighbors(g, l)
        end

        context("Random neighbor should be a neighbor") do
          g0 = random_neighbor(g, l)
          @fact (g - g0) |> sum |> abs --> 1
        end

        if typeof(l) != NKLandscape
          context("Neutral neighbors should have the same fitness") do
            function test_neutral_neighbors(genotype, landscape)
              nbrs = neutral_neighbors(genotype, landscape)
              score = fitness(genotype, landscape)
              @fact size(nbrs)[1] --> landscape.n
              @fact size(nbrs)[2] --> greater_than_or_equal(0)
              for j = 1:size(nbrs)[2]
                @fact fitness(nbrs[:,j], landscape) --> score
              end
            end

            test_neutral_neighbors(g, l)
          end
        end

        context("Fitter neighbors should all be fitter") do
          function test_fitter_neighbors(nbrs, landscape, score)
            @fact size(nbrs)[1] --> landscape.n
            @fact size(nbrs)[2] --> greater_than_or_equal(0)
            for j = 1:size(nbrs)[2]
              @fact fitness(nbrs[:,j], landscape) --> greater_than(score)
            end
          end

          fn = fitter_neighbors(g, l)
          test_fitter_neighbors(fn, l, f)
        end

        context("Fittest n neighbors should be all neighbors") do
          fnn = fittest_neighbors(g, l, l.n)
          @fact size(fnn)[1] --> l.n
          @fact size(fnn)[2] --> number_neighbors(g, l)
        end

        context("Fittest 1 neighbor should be the fittest neighbor") do
          nbrs = fitter_neighbors(g, l)
          if length(nbrs) > 0
            nbr1 = fittest_neighbor(g, l)
            @fact nbr1 --> nbrs[:,end]
          end
        end

        context("Adaptive walks") do
          function test_adaptive_walk(walk_function, genotype, landscape)
            walk = walk_function(genotype, landscape)
            @fact walk.length --> greater_than_or_equal(0)
            for i = 2:walk.length
              @fact fitness(walk.history_list[:,i], landscape) --> greater_than(fitness(genotype, landscape))
            end
          end

          context("Random adaptive walk should terminate and move uphill") do
            test_adaptive_walk(random_adaptive_walk, g, l)
          end

          context("Greedy adaptive walk should terminate and move uphill") do
            test_adaptive_walk(greedy_adaptive_walk, g, l)
          end

          context("Reluctant adaptive walk should terminate and move uphill") do
            test_adaptive_walk(reluctant_adaptive_walk, g, l)
          end
        end

        context("Neutral walks") do
          function test_neutral_walk(walk_function, genotype, landscape)
            walk = walk_function(genotype, landscape)
            @fact walk.length --> greater_than_or_equal(0)
            f = fitness(genotype, landscape)
            for i = 1:walk.length
              f0 = fitness(walk.history_list[:,i], landscape)
              @fact f0 --> roughly(f)
            end
          end

          context("Neutral walk should terminate and move sideways") do
            test_neutral_walk(neutral_walk, g, l)
          end
        end
      end

      context("Populations") do
        context("Should be the correct size") do
          function test_population_size(p::Population, l::Landscape, n)
            @fact popsize(p) --> n
            @fact size(p)[1] --> l.n
            @fact size(p)[2] --> n
          end

          for n = [1, 10]
            rp = rand(Population, l, n)
            test_population_size(rp, l, n)
            zp = zeros(Population, l, n)
            test_population_size(zp, l, n)
          end
        end

        context("Should compute fitnesses") do
          function test_population_fitnesses(p::Population, l::Landscape)
            fs = popfits(p, l)
            for i = 1:popsize(p)
              @fact fs[i] --> fitness(p[:,i], l)
            end
          end

          for n = [1, 10]
            rp = rand(Population, l, n)
            test_population_fitnesses(rp, l)
            zp = zeros(Population, l, n)
            test_population_fitnesses(zp, l)
          end
        end
      end

      context("Meta populations") do
        context("Should be the correct size") do
          function test_meta_population_size(p::MetaPopulation, l::Landscape, n, k)
            @fact popsize(p) --> n
            @fact popct(p) --> k
            @fact size(p)[1] --> l.n
            @fact size(p)[2] --> n
            @fact size(p)[3] --> k
          end

          for n = [1, 10]
            for k = [1, 10]
              rp = rand(MetaPopulation, l, n, k)
              test_meta_population_size(rp, l, n, k)
              zp = zeros(MetaPopulation, l, n, k)
              test_meta_population_size(zp, l, n, k)
            end
          end
        end

        context("Should compute fitness") do
          function test_meta_population_fitnesses(p::MetaPopulation, l::Landscape)
            fs = popfits(p, l)
            for ip = 1:popct(p)
              for ig = 1:popsize(p)
                @fact fs[ig,ip] --> fitness(p[:,ig,ip], l)
              end
            end
          end

          for n = [1, 10]
            for k = [1, 10]
              rp = rand(MetaPopulation, l, n, k)
              test_meta_population_fitnesses(rp, l)
              zp = zeros(MetaPopulation, l, n, k)
              test_meta_population_fitnesses(zp, l)
            end
          end
        end
      end

      for n = [1, 10, 100]
        context("Selection") do
          context("Should result in a population of the same size") do
            rp = rand(Population, l, n)
            np = propsel(rp, l)
            @fact popsize(np) --> n
          end

          context("Should contain only members of the original population") do
            rp = rand(Population, l, n)
            np = propsel(rp, l)
            for i_n = 1:n
              founds = map(1:n) do i_o
                np[:,i_n] == rp[:,i_o]
              end
              @fact any(founds) --> true
            end
          end

          context("Should create a new population when using propsel") do
            rp = rand(Population, l, n)
            np = propsel(rp, l)
            @fact is(rp, np) --> false
          end
        end

        context("Mutation") do
          context("Should result in a population of the same size") do
            rp = rand(Population, l, n)
            np = bwmutate(rp, l, 1.0)
            @fact popsize(np) --> n
          end

          context("Should mutate population members") do
            rp = rand(Population, l, n)
            np = bwmutate(rp, l, 1.0)
            @fact np --> not(rp)
          end
        end

        context("Migration") do
          context("Should move genotypes between populations") do
            rp = rand(MetaPopulation, l, n, 4)
            np = migrate(rp, 1.0, n / 2.0 |> ceil |> Int)
            @fact np --> not(rp)
          end
        end

        context("Linear migration") do
          context("Should move genotypes between populations") do
            rp = rand(MetaPopulation, l, n, 4)
            np = linmigrate(rp, 1.0, n / 2.0 |> ceil |> Int)
            @fact np --> not(rp)
          end
        end
      end
    end
  end

  context("Enumeration") do
    l = NKLandscape(2, 0)

    function testlevel(f, b)
      if f < 0.5
        @fact b --> 0
      else
        @fact b --> 1
      end
    end

    function testfitness(f, g)
      fg = fitness(g, l)
      @fact f --> roughly(fg)
    end

    context("lsfits") do
      fits = lsfits(l)

      testfitness(fits[1], [1, 1])
      testfitness(fits[2], [1, 2])
      testfitness(fits[3], [2, 1])
      testfitness(fits[4], [2, 2])
    end

    context("fitlevs") do
      levs = fitlevs(l, 2)

      f1 = fitness([1, 1], l)
      testlevel(f1, levs[1])

      f2 = fitness([1, 2], l)
      testlevel(f2, levs[2])

      f3 = fitness([2, 1], l)
      testlevel(f3, levs[3])

      f4 = fitness([2, 2], l)
      testlevel(f4, levs[4])
    end

    context("levcounts") do
      cts = levcounts(l, 2)

      gfits = [
        fitness([1, 1], l),
        fitness([1, 2], l),
        fitness([2, 1], l),
        fitness([2, 2], l)
      ]
      low = 0
      high = 0
      for gf = gfits
        if gf < 0.5
          low += 1
        else
          high += 1
        end
      end

      @fact cts[1] --> low
      @fact cts[2] --> high
    end

    context("neighbors") do
      la = NKLandscape(2, 0, a=3)
      @fact_throws neighbors(0, la)

      nbrs = neighbors(IntGenotype(0), l)
      @fact length(nbrs) --> 2
      @fact IntGenotype(1) --> anyof(nbrs...)
      @fact IntGenotype(2) --> anyof(nbrs...)
    end
  end
end

include("fastfunc.jl")
include("fastnets.jl")


text:
### Independent recovery on vintage 5
## author: Ziyi Yin (ziyi.yin@gatech.edu)

using DrWatson
@quickactivate "CompassTimeLapseCCS"
using TimeProbeSeismic
using JLD2, FFTW, DSP, PyPlot, Images

JLD2.@load "../model/CompassTimeLapse.jld2"
JLD2.@load "../data/d_lin_cut.jld2" d_lin_cut q

nv = 5

isdir("../results") || mkdir("../results")

Random.seed!(nv)

srcGeometry = q.geometry
recGeometry_stack = [d_lin_cut[i].geometry for i = 1:num_vintage]

model0_stack = [Model(n,d,o,m0_stack[i]; nb=100) for i = 1:num_vintage]

opt = Options(isic=true, limit_m=true)

ntComp = get_computational_nt(srcGeometry, recGeometry_stack[1], model0_stack[1])
info = Info(prod(n), d_lin_cut[1].nsrc, ntComp)

F0_stack = [judiProjection(info, recGeometry_stack[i])*judiModeling(info, model0_stack[i]; options=opt)*adjoint(judiProjection(info, srcGeometry)) for i = 1:num_vintage]

# Preconditioner

function dip(x,n;k=20)
    image = reshape(x,n)
    image_ext = zeros(Float32,n[1],2*n[2])
    image_ext[:,1:n[2]] = image
    image_ext[:,n[2]+1:end] = reverse(image,dims=2)
    image_f = fftshift(fft(image_ext))
    mask = ones(Float32,n[1],2*n[2])
    for i = 1:n[1]
        for j = 1:2*n[2]
            if (i-(n[1]+1)/2-k*j+k*(2*n[2]+1)/2)*(i-(n[1]+1)/2+k*j-k*(2*n[2]+1)/2)>0
                mask[i,j] = 0f0
            end
        end
    end
    mask = convert(Array{Float32},imfilter(mask,Kernel.gaussian(20)))
    image_f1 = mask .* image_f
    image_out = (vec(real.(ifft(ifftshift(image_f1)))[:,1:n[2]])+vec(real.(ifft(ifftshift(image_f1)))[:,end:-1:n[2]+1]))/2f0
    return image_out
end

D = joLinearFunctionFwd_T(prod(n), prod(n),
                                 v -> dip(v,n;k=10),
                                 w -> dip(w,n;k=10),
                                 Float32,Float32,name="dip filter")
Tm = judiTopmute(model0_stack[1].n, idx_wb, 1)  # Mute water column
S = judiDepthScaling(model0_stack[1])  # depth scaling
Mr = S*Tm*D

# Linearized Bregman parameters
nn = prod(model0_stack[1].n)
x = zeros(Float32, nn)
z = zeros(Float32, nn)

batchsize = 8
niter = 23

# Soft thresholding functions and Curvelet transform
soft_thresholding(x::Array{Float64}, lambda) = sign.(x) .* max.(abs.(x) .- convert(Float64, lambda), 0.0)
soft_thresholding(x::Array{Float32}, lambda) = sign.(x) .* max.(abs.(x) .- convert(Float32, lambda), 0f0)

n = model0_stack[1].n
C0 = joCurvelet2D(n[1], 2*n[2]; zero_finest = false, DDT = Float32, RDT = Float64)

function C_fwd(im, C, n)
	im = hcat(reshape(im, n), reshape(im, n)[:, end:-1:1])
	coeffs = C*vec(im)/sqrt(2f0)
	return coeffs
end

function C_adj(coeffs, C, n)
	im = reshape(C'*coeffs, n[1], 2*n[2])
	return vec(im[:, 1:n[2]] .+ im[:, end:-1:n[2]+1])/sqrt(2f0)
end

C = joLinearFunctionFwd_T(size(C0, 1), n[1]*n[2],
                          x -> C_fwd(x, C0, n),
                          b -> C_adj(b, C0, n),
                          Float32,Float64, name="Cmirrorext")       # mirror-extended curvelet

src_list = collect(1:d_lin_cut[1].nsrc)
lambda = 0f0

ps = 32     # number of probing vectors

# Main loop
for j = 1:niter
    # Select batch and set up left-hand preconditioner
    length(src_list) < batchsize && (global src_list = collect(1:d_lin_cut[nv].nsrc))
    src_list = src_list[randperm(MersenneTwister(nv+2000*j),length(src_list))]
    if j == 1
        global i = [pop!(src_list) for b=1:2*batchsize]
    else
        global i = [pop!(src_list) for b=1:batchsize]
    end
    println("Vintage $(nv) LS-RTM Iteration: $(j), imaging sources $(i)")
    flush(Base.stdout)

    J_probe = judiJacobian(F0_stack[nv][i],q[i],ps,d_lin_cut[nv][i])

    residual = J_probe*Mr*x-d_lin_cut[nv][i]
    phi = 0.5 * norm(residual)^2
    g = Mr'*J_probe'*residual

    # Step size and update variable
    t = Float32.(2*phi/norm(g)^2)

    # Update variables and save snapshot
    global z -= t*g
    C_z = C*z
    (j==1) && (global lambda = quantile(abs.(C_z), .9))   # estimate thresholding parameter in 1st iteration
    global x = adjoint(C)*soft_thresholding(C_z, lambda)

    @printf("At iteration %d function value is %2.2e and step length is %2.2e \n", j, phi, t)
    @printf("Lambda is %2.2e \n", lambda)

    JLD2.@save "../result/IndpRecovVintage$(nv)Iter$(j).jld2" x z g lambda phi
end


text:
using AbstractPlotting
using AbstractPlotting.MakieLayout
using GLMakie

# GUI for spirograph
# https://en.wikipedia.org/wiki/Spirograph

function x(t; R=1, k=1/4, l=1/4)
    R*[(1-k)*cos(t) + l*k*cos((1-k)/k*t), (1-k)*sin(t) - l*k*sin((1-k)/k*t)]
end

# where we lay our scene:
scene, layout = layoutscene()

flyt = GridLayout()
flyt.halign[] = :left # fails?
flyt.valign[] = :top

layout[1,1] = flyt
p = layout[1,2] = LAxis(scene)
rowsize!(layout, 1, Relative(1))
colsize!(layout, 2, Relative(2/3))

flyt[1,1] = LText(scene, "t")
ts = flyt[1,2] = LSlider(scene, range = 2pi:pi/8:40pi)

flyt[2, 1] = LText(scene, "k = r/R")
k = flyt[2,2] = LSlider(scene, range = 0.01:0.01:1.0, startvalue=1/4)

flyt[3,1] = LText(scene, "l=ρ/r")
l = flyt[3,2] = LSlider(scene,  range = 0.01:0.01:1.0, startvalue=1/4)

data = lift(ts.value, k.value, l.value) do ts,k,l

    ts′ = range(0, ts, length=1000)
    xys = Point2f0.(x.(ts′, R=1, k=k, l=l))

end

lines!(p, data)
xlims!(p, (-1, 1))
ylims!(p, (-1, 1))

scene


text:
"""
$(TYPEDEF)

Basic convolutional module with batch normalization.

## Fields

$(TYPEDFIELDS)

## Notes

In [MobileNetV3], the activations are `relu` and `hardσ`, while in EfficientNet,
they are `swish` and `σ`.
"""
@with_kw struct ConvBlock <: ModuleSpec
    "kernel size"
    ksize::Int
    "number of input kernels"
    k_in::Int
    "number of output kernels"
    k_out::Int
    "non-linearity/activation function (applied after BN)"
    σ = relu
    "stride of the convolution"
    stride::Int = 1
end
function (cb::ConvBlock)()
    @unpack ksize, k_in, k_out, σ, stride = cb
    return Chain(
        Conv((ksize, ksize), k_in => k_out; pad = ksize ÷ 2, stride = stride),
        BatchNorm(k_out, σ),
    )

end

"""
$(TYPEDEF)

A squeeze-and-excitation block.

## Fields

$(TYPEDFIELDS)

## Notes

In [MobileNetV3], the activations are `relu` and `hardσ`, while in EfficientNet,
they are `swish` and `σ`.
"""
@with_kw struct SqueezeExcitation <: ModuleSpec
    "Number of input kernels and output kernels"
    k::Integer
    "Ratio to calculate `k_mid`"
    ratio::Integer = 4
    "Ratio to calculate `k_mid`"
    k_mid::Integer = k ÷ ratio
    "Activation function to apply after \"squeezing\""
    σ1 = relu
    "Activation function to apply after \"excitation\""
    σ2 = relu
end
function (se::SqueezeExcitation)()
    @unpack k, k_mid, σ1, σ2 = se
    return SkipConnection(
        Chain(
            GlobalMeanPool(),
            Conv((1, 1), k => k_mid),
            BatchNorm(k_mid, σ1),
            Conv((1, 1), k_mid => k),
            BatchNorm(k, σ2),
        ),
        bmul,
    )
end

"""
$(TYPEDEF)

A depthwise separable convolution block.

## Fields

$(TYPEDFIELDS)
"""
@with_kw struct DepthwiseSeparable <: ModuleSpec
    "Kernel size of convolution"
    ksize::Integer
    "Number of input kernels"
    k_in::Integer
    "Number of expanded kernels"
    k_out::Integer = 6k_in
    "Activation function to apply after convolution"
    σ = relu
    "Stride of convolution"
    stride::Integer = 1
end

function (dws::DepthwiseSeparable)()
    @unpack ksize, k_in, k_out, stride = dws
    return Chain(
        Conv((1, 1), k_in => k_out),
        BatchNorm(k_out, σ),

        DepthwiseConv(
            (ksize, ksize),
            k_out => k_out;
            pad = ksize ÷ 2,
            stride = stride
        ),
        BatchNorm(k_out, σ)
    )
end

"""
$(TYPEDEF)

Mobile Inverted Bottleneck block, as described in
[Searching for MobileNetV3](https://arxiv.org/abs/1905.02244).

Also used in [EfficientNet: Rethinking Model Scaling for Convolutional Neural Networks](https://arxiv.org/abs/1905.11946v3)

## Fields

$(TYPEDFIELDS)
"""
@with_kw struct MBConv <: ModuleSpec
    "Kernel size of expansion block"
    ksize::Int
    "Number of input kernels"
    k_in::Int
    "Number of kernels after expansion"
    k_exp::Int
    "Number of output kernels"
    k_out::Int
    "Non-linearity/activation function to use in expansion block"
    σ = relu
    "Stride of the expansion block"
    stride = 1
    "Whether to add a squeeze-excitation block between expansion and projection"
    has_se = true
    """
    Convolution block to use, either [`DepthwiseSeparable`](@ref) or
    [`ConvBlock`](@ref).
    Depthwise separable convolutions are used in the literature, but as of
    `Flux@0.10.4`, the `DepthwiseConv` layer does not run on GPU, see
    [Flux.jl Issue #459](https://github.com/FluxML/Flux.jl/issues/459).
    """
    convblock = ConvBlock
end
function (mbc::MBConv)()
    @unpack ksize, k_in, k_exp, k_out, σ, stride, has_se, convblock = mbc

    block = Chain(
        # expansion
        convblock(ksize = ksize, k_in = k_in, k_out = k_exp; σ = σ, stride = stride)(),
        has_se ? SqueezeExcitation(k = k_exp)() : identity,

        # projection
        Conv((1, 1), k_exp => k_out),
    )

    # add residual connection if input and output are the same size
    # as done in MobileNetV3
    if k_in == k_out && stride == 1
        block = Residual(block)
    end

    return block

end

# Utilities

"""
    bmul(x, weights) = x .* weights

Broadcasted multiplication.
"""
bmul(x, weights) = x .* weights


text:
module Birdsong
using STFT
using WAV
using FDA
using DataFrames
using ProgressMeter
using StatsBase
using ArgCheck

#Universal things that are package internals
include("constants.jl")
#include("typealiases.jl")

#Functions that are package internals
include("dtw/dtw.jl")

#Functions that are exported for the user
include("find_sound.jl")
include("readsongs.jl")
include("sonogram.jl")
include("syllableFinder.jl")
include("summaryStats.jl")
include("callTrace.jl")
include("tonesmooth.jl")
include("visualize3d.jl")

end # module