text:
using UPPE

using PhysicalConstants.CODATA2018
const C0 = SpeedOfLightInVacuum.val
const EPS0 = VacuumElectricPermittivity.val

cd(@__DIR__)

include("air.jl")
UPPE.permittivity(w) = permittivity(w)
UPPE.permeability(w) = permeability(w)

# ******************************************************************************
grid = GridRT(
    GPU();
    ru = 1e-3,   # [m] unit of space in r direction
    rmax = 10,   # [ru] area in spatial domain
    Nr = 500,   # number of points in spatial domain
    rguard = 1,   # [ru] the width of the lossy slab at the end of r grid
    kguard = 45,   # [degrees] the cut-off angle for wave vectors
    tu = 1e-15,   # [s] unit of time
    tmin = -200,   # [tu] left boundary of t grid
    tmax = 200,   # [tu] right boundary of t grid
    Nt = 2^11,   # number of points on t grid
    tguard = 20,   # [tu] the width of the lossy slab at the ends of t grid
    wguard = 1e16,   # [1/s] the cut-off angular frequency
)

# ******************************************************************************
lam0 = 0.8e-6   # [m] central wavelength
a0 = 2e-3 / (2 * sqrt(log(2)))  # [m] initial beam radius
tau0 = 35e-15 / (2 * sqrt(log(2)))  # [s] initial pulse duration
W = 2e-3   # [J] initial pulse energy

w0 = 2 * pi * C0 / lam0   # central frequency
I0 = W / (pi^1.5 * tau0 * a0^2)   # initial pulse intensity
Iu = I0   # [W/m**2] unit of intensity
n0 = refractive_index(w0)
Eu = intensity2field(Iu, n0)

E = zeros((grid.Nr, grid.Nt))
for it=1:grid.Nt
for ir=1:grid.Nr
    E[ir,it] = sqrt(I0 / Iu) *
               exp(-0.5 * (grid.r[ir] * grid.ru)^2 / a0^2) *
               exp(-0.5 * (grid.t[it] * grid.tu)^2 / tau0^2) *
               cos(w0 * grid.t[it] * grid.tu)
end
end
field = Field(grid, w0, E; Eu, Iu)

# ******************************************************************************
function kerr_response_func!(ne, Pnl, J, E, p, z)
    Rk, = p
    @. Pnl = Rk * real(E)^3
    return nothing
end

function kerr_response(grid::Grid{T}, field, n2) where T
    (; arch) = grid
    (; Eu, w0, E) = field

    chi3 = chi3_func(w0, n2)
    Rk = EPS0 * chi3 * Eu^3
    p_resp = (Rk, )

    # response arrays:
    ne = adapt_array(arch, zeros(size(E)))
    Pnl = zero(E)
    J = zero(E)

    TF = typeof(kerr_response_func!)
    TP = typeof(p_resp)
    return Response{T, TF, TP}(ne, Pnl, J, kerr_response_func!, p_resp)
end

response =  kerr_response(grid, field, n2)

model = Model(
    grid, field, response;
    zu = 1,   # [m] unit of space in z direction
    z = 0,   # [m] initial propagation distance
    zmax = 5,   # [zu] propagation distance
    kparaxial = true,   # paraxial approximation for the linear term
    qparaxial = true,   # paraxial approximation for the nonlinear term
    nonlinearity = true,   # switch for the nonlinear term
)

simulation = Simulation(
    model, n2, N0;
    prefix = "results/",
    dz0 = 0.01,   # [zu] initial z step
    dzout = 0.5,   # [zu] z step for writing field into output file
    phimax = pi/200,  # maximum nonlinear phase increment for adaptive z step
    Icrit = 10,   # [Iu] maixmum allowed intensity (stop if exceeded)
    alg = RK4(),   # Algorithm of ODE solver for nonlinear part
)

run!(simulation)


text:
export leastsq,∂leastsq
"least square"
function leastsq(obj,pre)
  cost=sum((obj-pre).^2)/size(pre)[end]
end
function ∂leastsq(obj,pre)
  cost=pre-obj
end

ext: jl


text:
using DataFrames
using JLD2
using Base.Test

opt = get_options()

info("Load metadata from excel table stations_metadata.xlsx")

df_sel = load_metadata(opt["stat_sel"])

info("Clean the metadata table and save to stat_selected.csv")

df_nice = clean_metadata(df_sel)

writetable(Pkg.dir("IcanProj", "data", "stat_selected.csv"), df_nice)

info("Get data and save in WatershedData struct to wsh_info.jld2")

wsh_info = get_watershed_data(df_sel)

@save Pkg.dir("IcanProj", "data", "wsh_info.jld2") wsh_info

# Run some basic tests on wsh_info

for wsh_single in wsh_info

    @show wsh_single.name

    @test all(wsh_single.elev .>= 0)
    @test all(wsh_single.elev .<= 3000)
    @test size(wsh_single.gridcel) == (50, 50)

end


text:
"""
Datasets can load and transform data and can be used in a pipeline.
Typically they would have length and index methods implemented and need to return
only a single sample. The MiniBatch transformer would collect these single samples
and put them into a mini batch.

```julia
data = SomeDataset() |> SomeTransformer() |> MiniBatch()
train!(workout, data)
```

Please note that in case all your data would fit in memory, you can feed it
directly to the train! function (so a dataset is not required)

```julia
X = [randn(10,10,16) for i in 1:100]
Y = [randn(1,16) for i in 1:100]
train!(workout, zip(X,Y))
```

"""
abstract type Dataset end

"""
A dataset that contains random generated samples. Ideal for quick testing
of a model. The random values will be drawn from a normal distribution.

You need to provide the shape of X, Y and the number of samples that the
dataset should contain. Optionally you can specify a sleep value to simulate
IO blocking.

# Usage

```julia
xshape = (28,28,1)
yshape = (10,)
ds = TestDataset(xshape, yshape, 60000)

ds = TestDataset((100,), (1,), 1000, sleep=0.1)
```
"""
struct TestDataset <: Dataset
    x::Vector{Array}
    y::Vector{Array}
	sleep::Float64

	function TestDataset(shapeX::Tuple, shapeY::Tuple, samples::Int; sleep=0, dtype=Float32)
		x = [randn(dtype, shapeX...) for _ in 1:samples]
		y = [randn(dtype, shapeY...) for _ in 1:samples]
		new(x,y,sleep)
	end
end

Base.length(ds::TestDataset) = length(ds.x)

function Base.getindex(ds::TestDataset, idx)
	ds.sleep > 0 && sleep(ds.sleep)
	return (ds.x[idx], ds.y[idx])
end

"""
Dataset that contains two vectors, one for the input data and one for the labels.
"""
struct VectorDataset{A,B} <: Dataset
    x::A
    y::B
end

Base.length(ds::VectorDataset) = length(ds.x)
Base.getindex(ds::VectorDataset, idx) = return (ds.x[idx], ds.y[idx])

"""
Dataset that loads data from a JLD(2) file. The format
that is stored is expected to be in the shape of: (X,Y)

Note: Tested with JLD2 library only.

Example:
========

```julia
using JLD2

jldopen("example.jld2", "w") do file
    file["image1"] = (randn(Float32,28,28,1), rand(0:9,1))
    file["image2"] = (randn(Float32,28,28,1), rand(0:9,1))
end

f = jldopen("example.jld2", "r")
ds = JLD2Dataset(f)
ds[1]
```
"""
struct JLDDataset <: Dataset
	f
    keys

	JLDDataset(f) = new(f,keys(f))
end

Base.length(ds::JLDDataset) = length(ds.keys)

function Base.getindex(ds::JLDDataset, idx)
	key = ds.keys[idx]
	ds.f[key]
end

"""
Dataset that loads an single image from a file and optionally resizes the image.
The labels are passed as is.

# Usage

```julia
ds = ImageDataset(filenames, labels, resize=(200,200))
```
"""
struct ImageDataset <: Dataset
    filenames::Vector{String}
    labels::Vector
	resize::Union{Nothing,Tuple}

	function ImageDataset(filenames, labels; resize=nothing)
		@assert length(filenames) == length(labels)
		try
			@eval import Images
			@eval import ImageMagick
		catch
			@warn "Package Images or ImageMagick not installed"
		end
		new(filenames, labels, resize)
	end

end

Base.length(ds::ImageDataset) = length(ds.filenames)

function Base.getindex(ds::ImageDataset, idx)
	filename = ds.filenames[idx]
	# img = Images.load(filename) has multi-thread issue
    img = ImageMagick.load(filename)
	img = Images.channelview(img)
	img = permutedims(img, [2,3,1])
	img = convert(Array{Float32}, img)
	if ds.resize !== nothing
		img = Images.imresize(img, ds.resize...)
	end
	return (img, ds.labels[idx])
end

"""
Dataset that retrieves is data from a dataframe. The provided column names
for X and Y can be either a single Symbol or a Vector of Symbols.

# Usage

```julia
df = DataFrame(randn(4,20))
ds = DFDataset(df, :x1, :x2)
ds = DFDataset(df, [:x1, :x3], [:x2,:x4])
```
"""
struct DFDataset <: Dataset
	df
    X::Vector{Symbol}
	Y::Vector{Symbol}

	DFDataset(df,X,Y) = new(df, makeArray(X), makeArray(Y))
end

Base.length(ds::DFDataset) = size(ds.df,1)

function Base.getindex(ds::DFDataset, idx)
	df = ds.df
	(Vector(df[idx, ds.X]), Vector(df[idx, ds.Y]))
end


text:
struct ThreadTask
  p::Ptr{UInt}
end
Base.pointer(tt::ThreadTask) = tt.p

@inline taskpointer(tid::T) where {T} = THREADPOOLPTR[] + tid*(THREADBUFFERSIZE%T)

@inline function _call(p::Ptr{UInt})
  fptr = load(p + sizeof(UInt), Ptr{Cvoid})
  assume(fptr ≠ C_NULL)
  ccall(fptr, Cvoid, (Ptr{UInt},), p)
end
@inline function launch(f::F, tid::Integer, args::Vararg{Any,K}) where {F,K}
  p = taskpointer(tid)
  f(p, args...)
  state = _atomic_xchg!(p, TASK) # exchange must happen atomically, to prevent it from switching to `WAIT` after reading
  state == WAIT && wake_thread!(tid)
  return nothing
end

function (tt::ThreadTask)()
  p = pointer(tt)
  max_wait = one(UInt32) << 20
  wait_counter = max_wait
  GC.@preserve THREADPOOL begin
    while true
      if _atomic_state(p) == TASK
        # if _atomic_cas_cmp!(p, TASK, EXEC)
        _call(p)
        # store!(p, SPIN)
        _atomic_store!(p, SPIN)
        wait_counter = zero(UInt32)
        continue
      end
      pause()
      if (wait_counter += one(UInt32)) > max_wait
        wait_counter = zero(UInt32)
        _atomic_cas_cmp!(p, SPIN, WAIT) && Base.wait()
      end
    end
  end
end

# 1-based tid, pushes into task 2-nthreads()
@noinline function wake_thread!(_tid::T) where {T <: Integer}
  tid = _tid % Int
  tidp1 = tid + one(tid)
  assume(unsigned(length(Base.Workqueues)) > unsigned(tid))
  assume(unsigned(length(TASKS)) > unsigned(tidp1))
  @inbounds push!(Base.Workqueues[tidp1], TASKS[tid])
  ccall(:jl_wakeup_thread, Cvoid, (Int16,), tid % Int16)
end

@noinline function checktask(tid)
  t = TASKS[tid]
  if istaskfailed(t)
    display(t)
    dump(t)
    println()
    initialize_task(tid)
    return true
  end
  yield()
  false
end
# 1-based tid
@inline wait(tid::Integer) = wait(taskpointer(tid), tid)
@inline wait(p::Ptr{UInt}) = wait(p, (p - THREADPOOLPTR[]) ÷ (THREADBUFFERSIZE))
@inline function wait(p::Ptr{UInt}, tid)
  counter = 0x00000000
  while _atomic_state(p) == TASK
    pause()
    if ((counter += 0x00000001) > 0x00010000)
      checktask(tid) && return true
    end
  end
  false
end


text:
function compile_pricer(
            provider::MarketData.AbstractMarketDataProvider,
            pricing_date::Date,
            model::PricingModel,
            c::Contract,
            attr::ContractAttributes;
            compiler::Symbol=:interpreter
        )

    Compiler.compile_pricer(model, provider, attr, pricing_date, c, AbstractPricer, compiler=compiler)
end

function compile_cashflow_pricer(
            provider::MarketData.AbstractMarketDataProvider,
            pricing_date::Date,
            model::PricingModel,
            c::Contract,
            attr::ContractAttributes;
            compiler::Symbol=:interpreter
        )

	Compiler.compile_pricer(model, provider, attr, pricing_date, c, AbstractCashflowPricer, compiler=compiler)
end

function compile_pricing_program(
            provider::MarketData.AbstractMarketDataProvider,
            pricing_date::Date,
            model::PricingModel,
            c::Contract,
            attr::ContractAttributes;
        )

    Compiler.compile_pricing_program(model, provider, attr, pricing_date, c, AbstractPricer)
end

function compile_cashflow_pricing_program(
            provider::MarketData.AbstractMarketDataProvider,
            pricing_date::Date,
            model::PricingModel,
            c::Contract,
            attr::ContractAttributes;
        )

    Compiler.compile_pricing_program(model, provider, attr, pricing_date, c, AbstractCashflowPricer)
end

function create_pricer(
        compiler_result;
        compiler::Symbol=:interpreter,
        memory_buffer::Union{Nothing, Vector{Any}}=nothing,
        input_values_buffer::Union{Nothing, Vector{Any}}=nothing,
        auto_resize_buffers::Bool=true
    )

    return Compiler.create_pricer(compiler_result, compiler=compiler, memory_buffer=memory_buffer, input_values_buffer=input_values_buffer, auto_resize_buffers=auto_resize_buffers)
end

"""
    price(p::AbstractPricer, scenario::Scenario) :: Real

Returns the price for the contract.
"""
function price end

"""
    exposures(method::AbstractExposuresMethod, p::AbstractPricer, scenario::Scenario) :: ExposureResult

Risk factors mapping.
"""
function exposures end

"""
    riskfactors(pricer::AbstractPricer) -> itr

Returns an iterator for the risk factors of the contract.
"""
function riskfactors end

"""
	eachcashflow(p::AbstractCashflowPricer, scenario::Scenario) -> itr

Returns an iterator for the cashflow of the contract.
"""
function eachcashflow end


text:
f = open("input07.txt", "r")
    l = readlines(f)
close(f)

positions = [parse(Int, n) for n in split(l[1], ",")]

options = collect(minimum(positions):maximum(positions))
consumption = zeros(length(options))
for (i, k) in enumerate(options)
    diff = positions - ones(length(positions)) * k
    fuel = sum([abs(j) for j in diff])
    consumption[i] = fuel
end

best_fuel = minimum(consumption)
println("Result 1: $best_fuel")

consumption = zeros(length(options))
for (i, k) in enumerate(options)
    diff = positions - ones(length(positions)) * k
    eachfuel = [sum(collect(1:abs(j))) for j in diff]
    fuel = sum([abs(j) for j in eachfuel])
    consumption[i] = fuel
end

new_best_fuel = minimum(consumption)
println("Result 2: $new_best_fuel")

max_stars_repo_name: jung-benjamin/advent-of-code

max_issues_repo_name: jung-benjamin/advent-of-code

max_forks_repo_name: jung-benjamin/advent-of-code


text:
module Namtso

using Dates,JSON,Random

export EChart,series!,public_render!

include("Base.jl")
include("ECharts/ECharts.jl")

end # module


text:
# License for this file: MIT (expat)
# Copyright 2017-2018, DLR Institute of System Dynamics and Control

"""
    module ContactDetectionMPR

Detects contact between solid objects3D with the improved MPR algorithm.

For further information see [A. Neumayr and M. Otter, Collision Handling with Variable-Step Integrators, EOOLT'17]

# Main developers
Andrea Neumayr and Martin Otter, [DLR - Institute of System Dynamics and Control](https://www.dlr.de/sr/en)
"""
module ContactDetectionMPR

using StaticArrays

import Modia3D
import Modia3D.Basics
import Modia3D.Composition
import Modia3D.Shapes
import Printf

using StaticArrays
using OrderedCollections

include("analyzeMPR.jl")
include("utilitiesMPR.jl")
include("utilitiesPairID.jl")
include("mpr.jl")
include("handler.jl")

end


text:
import Base: getindex, setindex!, (+), (*), (-), length, setindex!, show

export RV, RV_types, E, Var, validate!
export vals, probs, report
export Uniform_RV, Binomial_RV, Bernoulli_RV

"""
`RV` represents a discrete random variable with finite support.
"""
mutable struct RV{S<:Number, T<:Real}
  data::Dict{S,T}
  valid::Bool
  function RV{S,T}() where {S,T}
    d = Dict{S,T}()
    new(d,false)
  end
end

RV_types(X::RV{S,T}) where {S,T} = (S,T)

"""
`length(X::RV)` returns the number of values in the random variable `X`.
"""
length(X::RV) = length(X.data)

function show(io::IO, X::RV)
  S,T = RV_types(X)
  print(io, "RV{$S,$T} with $(length(X.data)) values")
end

"""
`vals(X::RV)` returns an iterator of the values this random variable
can take. Use `X[v]` to get the associate probability of the
value `v`.
"""
vals(X::RV) = keys(X.data)

"""
`probs(X::RV)` returns an iterator of the probabilities associated with
the values in `X`.
"""
function probs(X::RV)
  validate!(X)
  return values(X.data)
end

"""
`validate!(X)` ensures that the probabilies of the values in `X`
sum to one. If not, they are rescaled.
"""
function validate!(X::RV)
  if !X.valid
    u = sum(values(X.data))
    for x in keys(X.data)
      X.data[x]/= u
    end
    X.valid = true
  end
  nothing
end

"""
`E(X)` is the expected value of `X`.
"""
function E(X::RV{S,T}) where {S,T}
  @assert length(X)>0 "Cannot compute the expected value: no values!"
  validate!(X)
  return  sum( k*X.data[k] for k in keys(X.data))
end

"""
`Var(Y)` is the variance of `Y`.
"""
function Var(X::RV{S,T}) where {S,T}
  @assert length(X)>0 "Cannot compute the variance: no values!"
  validate!(X)
  exex = E(X)^2
  exx  = sum( k*k*X.data[k] for k in keys(X.data))
  return exx - exex
end

"""
`Bernoulli(p)` makes a single coin flip RV.
"""
function Bernoulli_RV(p::T) where {T}
  @assert 0<=p && p<=1 "p must be in [0,1]"
  X = RV{Int,T}()
  X[1] = p
  X[0] = 1-p
  X.valid = true
  return X
end

"""
`Binomial_RV(n,p)` returns a binomial random variable.
"""
function Binomial_RV(n::S, p::T) where {S<:Integer,T}
  @assert n>=0 "n must be nonnegative"
  @assert 0<=p && p<=1 "probability must be in [0,1]"
  X = RV{S,T}()
  for k=0:n
    X[k] = binomial(n,k)*(p^k)*(1-p)^(n-k)
  end
  return X
end

"""
`Uniform_RV(n)` returns the uniform distribution on `{1,2,...,n}`.
"""
function Uniform_RV(n::Int)
  X = RV{Int,Rational{Int}}()
  for k=1:n
    X[k] = 1//n
  end
  return X
end

"""
`X[v]` returns the probability of `v` in the random variable `X`.
Note that we validate `X` (with `validate!`) before retrieving
the value.
"""
function getindex(X::RV{S,T}, k::S) where {S,T}
  validate!(X)
  try
    return X.data[k]
  catch
    return zero(T)
  end
end

function setindex!(X::RV{S,T}, p::Real, k::S) where {S,T}
  @assert p>=0 "Probability must be nonnegative"
  X.data[k] = T(p)
  X.valid = false
  return p
end

"""
`X+Y`: sum of independent random variables.
"""
function (+)(X::RV, Y::RV)
  S = typeof( first(vals(X)) + first(vals(Y)) )
  T = typeof( first(probs(X)) + first(probs(Y)) )

  Z = RV{S,T}()
  for a in keys(X.data)
    for b in keys(Y.data)
      if !haskey(Z.data,a+b)
        Z.data[a+b] = 0
      end
      Z.data[a+b] += X.data[a]*Y.data[b]
    end
  end
  validate!(Z)
  return Z
end

"""
`X-Y`: difference of independent random variables.
"""
(-)(X::RV,Y::RV) = X + (-Y)

"""
`a*X`: scalar multiple of the random variable `X`.
"""
function (*)(a::Number, X::RV)
  S = typeof( first(vals(X)) * a)
  T = typeof( first(probs(X)) )
  aX = RV{typeof(a),T}()
  for k in vals(X)
    aX[a*k] = X[k]
  end
  return aX
end

"""
`-X`: negative of a random variable.
"""
function (-)(X::RV{S,T}) where {S,T}
  negone = - one(S)
  return negone * X
end

# This implementation is somewhat inefficient
"""
`random_choice(X)` for a random variable `X` returns a value of `X`
according to its probability distribution. That is, the probability
a value `v` is returned is `X[v]`.
"""
function random_choice(X::RV)
  validate!(X)
  return random_choice(X.data)
end

"""
`report(X)` prints out a list of the values of `X` and their
associated probabilities
"""
function report(X::RV)
  A = collect(vals(X))
  try
    sort!(A)
  finally
    for a in A
      println("$a\t$(X[a])")
    end
  end
  nothing
end


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

import UnitSystems: listext, Kinematic, Mechanical, Electromagnetic, Thermodynamic, Molar, Photometric

# dimension

struct Dimension{D} <: AbstractModule
    @pure Dimension{D}() where D = new{D}()
end

@pure Dimension(D::Group) = Dimension{D}()
@pure isunknown(x::Dimension{D}) where D = isunknown(D)
@pure isdimension(x::Dimension{D}) where D = typeof(D)<:Dimension
@pure dimension(::Dimension{D}) where D = D

value(::Dimension{D}) where D = value(D)

showgroup(io::IO,x::Dimension{D},u) where D = showgroup(io,D,u,'𝟙')
showgroup(io::IO,x,u,c) = showgroup(io,x,u)

Base.show(io::IO,x::Dimension{D}) where D  = show(io,D)

Base.:(==)(::Dimension{A},::Dimension{B}) where {A,B} = A == B

preferred(a,b) = a
function Base.:+(a::Dimension{A},b::Dimension{B}) where {A,B}
    if A.v==B.v
        preferred(a,b)
    elseif islog(A) && islog(B)
        Dimension{A+B}()
    else
        error("addition of Group $a + $b is not valid")
    end
end
function Base.:-(a::Dimension{A},b::Dimension{B}) where {A,B}
    islog(A) && islog(B) ? Dimension{A-B}() : a+b
end

Base.sqrt(::Dimension{D}) where D = Dimension{sqrt(D)}()
Base.cbrt(::Dimension{D}) where D = Dimension{cbrt(D)}()
Base.inv(::Dimension{D}) where D = Dimension{inv(D)}()
Base.log(x::Dimension{D}) where D = Dimension{log(D)}()
Base.log2(x::Dimension{D}) where D = Dimension{log2(D)}()
Base.log10(x::Dimension{D}) where D = Dimension{log10(D)}()
Base.log(b::Number,x::Dimension{D}) where D = Dimension{log(b,D)}()
Base.exp(x::Dimension{D}) where D = Dimension{exp(D)}()
Base.exp2(x::Dimension{D}) where D = Dimension{exp2(D)}()
Base.exp10(x::Dimension{D}) where D = Dimension{exp10(D)}()
Base.:^(::UnitSystems.Constant{B},x::Dimension{D}) where {B,D} = Dimension{B^D}()
Base.:^(b::T,x::Dimension{D}) where {T<:Number,D} = Dimension{b^D}()
Base.:^(::Dimension{D},b::Number) where D = Dimension{D^b}()
Base.:^(::Dimension{D},b::Integer) where D = Dimension{D^b}()
Base.:^(::Dimension{D},b::Rational) where D = Dimension{D^b}()
Base.:*(::Dimension{A},::Dimension{B}) where {A,B} = Dimension{A*B}()
Base.:/(::Dimension{A},::Dimension{B}) where {A,B} = Dimension{A/B}()
Base.:*(a::Number,b::Dimension{D}) where D = Dimension{a*D}()
Base.:*(a::Dimension{D},b::Number) where D = Dimension{D*b}()
Base.:/(a::Number,b::Dimension) = a*inv(b)
Base.:/(a::Dimension,b::Number) = a*inv(b)

const isq = Values('F','M','L','T','Q','Θ','N','J','A','R','C')
const dims = length(isq)
factorfind(x,k,i=0) = iszero(x) ? (x,0) : (r = x%k; iszero(r) ? factorfind(x÷k,k,i+1) : (x,i))

const 𝟙 = Dimension(valueat(0,dims))
for i ∈ 1:dims
    @eval const $(Symbol(isq[i])) = Dimension(valueat($i,dims))
end
const usq = Values(F,M,L,T,Q,Θ,N,J,A,R,C)

@pure function factorize(x::Float64,A::Group{Int,N},B::Group{Int,N},C::Group{Int,N},D::Group{Int,N},E::Group{Int,N},F::Group{Int,N},G::Group{Int,N},H::Group{Int,N}) where N
    if isinteger(x)
        try
            return factorize(Int(x),A,B,C,D,E,F,G,H)
        catch
        end
    end
    (x,a) = factorfind(x,2π)
    (x,b) = factorfind(x,π/200)
    (x,c) = factorfind(x,π/180)
    Group(zeros(Values{N,Int}),x)*A^a*(A/B^4/D^2)^b*(A/B^3/C^2/D)^c
end

@pure function factorize(x::Int,A::Group{Int,N},B::Group{Int,N},C::Group{Int,N},D::Group{Int,N},E::Group{Int,N},F::Group{Int,N},G::Group{Int,N},H::Group{Int,N}) where N
    (x,b) = factorfind(x,2)
    (x,c) = factorfind(x,3)
    (x,d) = factorfind(x,5)
    (x,e) = factorfind(x,7)
    (x,f) = factorfind(x,11)
    (x,g) = factorfind(x,19)
    (x,h) = factorfind(x,43)
    Group(zeros(Values{N,Int}),x)*B^b*C^c*D^d*E^e*F^f*G^g*H^h
end

const basis = Values("kB", "NA", "𝘩", "𝘤", "𝘦", "Kcd", "ΔνCs", "R∞", "α", "μₑᵤ", "μₚᵤ", "ΩΛ", "H0", "g₀", "aⱼ", "au", "ft", "ftUS", "lb", "T₀", "atm", "inHg", "RK90", "KJ90", "RK", "KJ", "Rᵤ2014", "Ωᵢₜ", "Vᵢₜ", "kG", "mP", "GME", "GMJ", "φ", "γ", "ℯ", "τ", "2", "3", "5", "7", "11", "19","43")
const vals = length(basis)
@pure function constant(d::Group,C::Coupling=UnitSystems.Universe); cs =
    UnitSystems.kB^makeint(d.v[1])*
    UnitSystems.NA^makeint(d.v[2])*
    UnitSystems.𝘩^makeint(d.v[3])*
    UnitSystems.𝘤^makeint(d.v[4])*
    UnitSystems.𝘦^makeint(d.v[5])*
    UnitSystems.Kcd^makeint(d.v[6])*
    UnitSystems.ΔνCs^makeint(d.v[7])*
    UnitSystems.R∞^makeint(d.v[8])*
    inv(UnitSystems.αinv)^makeint(d.v[9])*
    UnitSystems.μₑᵤ^makeint(d.v[10])*
    UnitSystems.μₚᵤ^makeint(d.v[11])*
    UnitSystems.ΩΛ^makeint(d.v[12])*
    UnitSystems.H0^makeint(d.v[13])*
    UnitSystems.g₀^makeint(d.v[14])*
    UnitSystems.aⱼ^makeint(d.v[15])*
    UnitSystems.au^makeint(d.v[16])*
    UnitSystems.ft^makeint(d.v[17])*
    UnitSystems.ftUS^makeint(d.v[18])*
    UnitSystems.lb^makeint(d.v[19])*
    UnitSystems.T₀^makeint(d.v[20])*
    UnitSystems.atm^makeint(d.v[21])*
    UnitSystems.inHg^makeint(d.v[22])*
    UnitSystems.RK1990^makeint(d.v[23])*
    UnitSystems.KJ1990^makeint(d.v[24])*
    UnitSystems.RK2014^makeint(d.v[25])*
    UnitSystems.KJ2014^makeint(d.v[26])*
    UnitSystems.Rᵤ2014^makeint(d.v[27])*
    UnitSystems.Ωᵢₜ^makeint(d.v[28])*
    UnitSystems.Vᵢₜ^makeint(d.v[29])*
    UnitSystems.kG^makeint(d.v[30])*
    UnitSystems.mP^makeint(d.v[31])*
    UnitSystems.GME^makeint(d.v[32])*
    UnitSystems.GMJ^makeint(d.v[33])*
    Base.MathConstants.φ^makeint(d.v[34])*
    Base.MathConstants.γ^makeint(d.v[35])*
    Base.MathConstants.ℯ^makeint(d.v[36])*
    (2π)^makeint(d.v[37]); is =
    2.0^makeint(d.v[38])*
    3.0^makeint(d.v[39])*
    5.0^makeint(d.v[40])*
    7.0^makeint(d.v[41])*
    11.0^makeint(d.v[42])*
    19.0^makeint(d.v[43])*
    43.0^makeint(d.v[44])
    return cs*(is*d.c)
end

# convertunit

struct ConvertUnit{D,U,S} <: AbstractModule
    @pure ConvertUnit{D,U,S}() where {D,U,S} = new{D,normal(U),normal(S)}()
end

convertdim(::ConvertUnit{D,U,S}) where {D,U,S} = convertdim(D,U,S)
convertdim(::Dimension{D},U,S) where D = convertdim(D,U,S)
dimconvert(x::T,d,u,s) where T = (isone(ratio(d,u,s)) ? zero(x) : x)::T
convertdim(d::Group{T},U,S) where T = Dimension{Group{T,11}(dimconvert.(d.v,usq,Ref(U),Ref(S)))}()

function Base.show(io::IO,::ConvertUnit{D,U,S}) where {D,U,S}
    d = convertdim(D,U,S)
    print(io, ratio(D,U,S), " [")
    showgroup(io,S(d),S)
    print(io, "]/[")
    showgroup(io,U(d),U)
    print(io, "] ", unitname(U), " -> ", unitname(S))
end

@pure Base.inv(::ConvertUnit{D,U,S}) where {D,U,S} = ConvertUnit{inv(D),U,S}()

(D::Dimension)(U::UnitSystem,S::UnitSystem) = ConvertUnit{D,U,S}()
(D::Dimension)(U::UnitSystem) = U(ratio(D,Natural,U),D)
(D::Dimension)(v::Real,U::UnitSystem,S::UnitSystem=Metric) = v/ratio(D,U,S) #U(_,D)
#(D::Dimension)(v::Number,U::UnitSystem,S::UnitSystem=Natural) = U(v/ratio(D,U,S),D)

Base.log(x::ConvertUnit{D,U,S}) where {D,U,S} = ConvertUnit{log(D),U,S}()
Base.log2(x::ConvertUnit{D,U,S}) where {D,U,S} = ConvertUnit{log2(D),U,S}()
Base.log10(x::ConvertUnit{D,U,S}) where {D,U,S} = ConvertUnit{log10(D),U,S}()
Base.log(b::Number,x::ConvertUnit{D,U,S}) where {D,U,S} = ConvertUnit{log(b,D),U,S}()
Base.exp(x::ConvertUnit{D,U,S}) where {D,U,S} = ConvertUnit{exp(D),U,S}()
Base.exp2(x::ConvertUnit{D,U,S}) where {D,U,S} = ConvertUnit{exp2(D),U,S}()
Base.exp10(x::ConvertUnit{D,U,S}) where {D,U,S} = ConvertUnit{exp10(D),U,S}()
Base.:^(a::Number,b::ConvertUnit{D,U,S}) where {D,U,S} = ConvertUnit{a^D,U,S}()
Base.:^(a::ConvertUnit{D,U,S},b::Integer) where {D,U,S} = Quantity{D^b,U,S}()
Base.:^(a::ConvertUnit{D,U,S},b::Rational{Int}) where {D,U,S} = Quantity{D^b,U,S}()
Base.:*(a::ConvertUnit{A,U,S},b::ConvertUnit{B,U,S}) where {A,B,U,S} = ConvertUnit{A*B,U,S}()
Base.:/(a::ConvertUnit{A,U,S},b::ConvertUnit{B,U,S}) where {A,B,U,S} = ConvertUnit{A/B,U,S}()

# quantity

struct Quantity{D,U,T} <: AbstractModule
    v::T
    Quantity{D,U,T}(v) where {T,D,U} = new{D,normal(U),T}(v)
    @pure Quantity{D,U,Int}(v) where {D,U} = new{D,normal(U),Int}(v)
    @pure Quantity{D,U,Float64}(v) where {D,U} = new{D,normal(U),Float64}(v)
end

Quantity{D,U}(v::Quantity) where {D,U} = Quantity{D,U}(v.v)
Quantity{D,U}(v::T) where {D,U,T} = Quantity{D,U,T}(v)
@pure Quantity{D,U}(v::Int) where {D,U} = Quantity{D,U,Int}(v)::Quantity{D,normal(U),Int}
@pure Quantity{D,U}(v::Float64) where {D,U} = Quantity{D,U,Float64}(v)::Quantity{D,normal(U),Float64}
@pure Quantity(D,U,v::Int) = Quantity{D,U}(v)
@pure Quantity(D,U,v::Float64) = Quantity{D,U}(v)
Quantity(D,U,v) = Quantity{D,U}(v)

Base.convert(::Type{Float64},x::Quantity) = convert(Float64,x.v)
Base.convert(::Type{Float64},x::Quantity{D,U,Float64}) where {D,U} = x.v
Base.convert(::Type{T},x::Quantity{D,U,T}) where {D,U,T<:Number} = x.v

import UnitSystems: evaldim
UnitSystems.isquantity(U::Quantity) = true

quantity(q) = q
quantity(q::Quantity) = q.v
@pure quantity(q::Quantity{D,U,Int} where {D,U}) = q.v
@pure quantity(q::Quantity{D,U,Float64} where {D,U}) = q.v
@pure dimensions(::Quantity{D}) where D = D
@pure Dimension(::Quantity{D}) where D = D
@pure unitsystem(::Quantity{D,U}) where {D,U} = dimension(U)
@pure unitsystem2(::Quantity{D,U}) where {D,U} = U

function Base.show(io::IO,x::Quantity{D,U}) where {D,U}
    print(io, x.v, " [")
    showgroup(io,U(D),U)
    print(io, "] ", unitname(U))
end

Base.log(x::Quantity{D,U}) where {D,U} = Quantity{log(D),U}(log(x.v))
Base.log2(x::Quantity{D,U}) where {D,U} = Quantity{log2(D),U}(log2(x.v))
Base.log10(x::Quantity{D,U}) where {D,U} = Quantity{log10(D),U}(log10(x.v))
Base.log(b::Number,x::Quantity{D,U}) where {D,U} = Quantity{log(b,D),U}(log(b,x.v))
Base.exp(x::Quantity{D,U}) where {D,U} = Quantity{exp(D),U}(exp(x.v))
Base.exp2(x::Quantity{D,U}) where {D,U} = Quantity{exp2(D),U}(exp2(x.v))
Base.exp10(x::Quantity{D,U}) where {D,U} = Quantity{exp10(D),U}(exp10(x.v))
#Base.:^(a::Constant,b::Quantity{D,U}) where {D,U} = Quantity{a^D,U}(a^b.v)
Base.:^(a::UnitSystems.Constant,b::Quantity{D,U}) where {D,U} = Quantity{a^D,U}(a^b.v)
Base.:^(a::Number,b::Quantity{D,U}) where {D,U} = Quantity{a^D,U}(a^b.v)
Base.:^(a::Quantity{D,U},b::Number) where {D,U} = Quantity{D^b,U}(a.v^b)
Base.:^(a::Quantity{D,U},b::Integer) where {D,U} = Quantity{D^b,U}(a.v^b)
Base.:^(a::Quantity{D,U},b::Rational{Int}) where {D,U} = Quantity{D^b,U}(a.v^b)
Base.:+(a::Number,b::Quantity{D,U}) where {D,U} = U(D)==𝟙 ? Quantity{D,U}(a+b.v) : throw(error("$(U(D)) ≠ 𝟙 "))
Base.:+(a::Quantity{D,U},b::Number) where {D,U} = U(D)==𝟙 ? Quantity{D,U}(a.v+b) : throw(error("$(U(D)) ≠ 𝟙 "))
Base.:-(a::Number,b::Quantity{D,U}) where {D,U} = U(D)==𝟙 ? Quantity{D,U}(a-b.v) : throw(error("$(U(D)) ≠ 𝟙 "))
Base.:-(a::Quantity{D,U},b::Number) where {D,U} = U(D)==𝟙 ? Quantity{D,U}(a.v-b) : throw(error("$(U(D)) ≠ 𝟙 "))
Base.:-(a::Quantity{D,U},b::Quantity{D,U}) where {D,U} = Quantity{D,U}(a.v-b.v)
Base.:+(a::Quantity{D,U},b::Quantity{D,U}) where {D,U} = Quantity{D,U}(a.v+b.v)
Base.:*(a::Real,b::Quantity{D,U}) where {D,U} = Quantity{D,U}(a*b.v)
Base.:*(a::Quantity{D,U},b::Real) where {D,U} = Quantity{D,U}(a.v*b)
Base.:*(a::Complex,b::Quantity{D,U}) where {D,U} = Quantity{D,U}(a*b.v)
Base.:*(a::Quantity{D,U},b::Complex) where {D,U} = Quantity{D,U}(a.v*b)
Base.:*(a::Quantity{A,U},b::Quantity{B,U}) where {A,B,U} = Quantity{A*B,U}(a.v*b.v)
Base.:/(a::Quantity{A,U},b::Quantity{B,U}) where {A,B,U} = Quantity{A/B,U}(a.v/b.v)
Base.:/(a::Quantity{D,U},b::Quantity{D,U}) where {D,U} = Quantity{D/D,U}(a.v/b.v)
Base.:/(a::Quantity{D,A},b::Quantity{D,B}) where {A,B,D} = ConvertUnit{(a.v/b.v)*D,A,B}()
Base.:/(a::Number,b::Quantity{D}) where D = a*inv(b)
Base.:/(a::Quantity{D},b::Number) where D = a*inv(b)
Base.:/(a::Quantity{D},b::UnitSystems.Constant) where D = a*inv(b)
Base.:/(a::UnitSystems.Constant,b::Quantity{D}) where D = a*inv(b)
Base.:/(a::Quantity{D,U},b::ConvertUnit{D,S,U}) where {D,U,S} = a*ConvertUnit{D,U,S}()
Base.:-(a::Quantity{D,U}) where {D,U} = Quantity{D,U}(-a.v)
Base.inv(a::Quantity{D,U}) where {D,U} = Quantity{inv(D),U}(inv(a.v))
Base.sqrt(a::Quantity{D,U}) where {D,U} = Quantity{sqrt(D),U}(sqrt(a.v))
Base.cbrt(a::Quantity{D,U}) where {D,U} = Quantity{cbrt(D),U}(cbrt(a.v))

Base.:*(a::Dimension{A},b::Quantity{D,U}) where {D,U,A} = Quantity{a,U}(10^-(isunknown(A) ? A.v.v[end] : A.v[end]))*b
Base.:*(a::Quantity{D,U},b::Dimension{B}) where {D,U,B} = a*Quantity{b,U}(10^-(isunknown(B) ? B.v.v[end] : B.v[end]))
Base.:*(a::Quantity{D,U},b::ConvertUnit{A,U,U}) where {D,U,A} = D==A || A==𝟙 ? a : error("ConvertUnit incompatible $D ≠ $A")
Base.:*(a::ConvertUnit{A,U,U},b::Quantity{D,U}) where {D,U,A} = D==A || A==𝟙 ? b : error("ConvertUnit incompatible $D ≠ $A")
Base.:*(a::ConvertUnit{A,S,U},b::Quantity{D,U}) where {D,U,A,S} = (A==D && S==U) ? b : inv(A)==D ? ConvertUnit{inv(A),U,S}()*b : error("ConvertUnit incompatible $(inv(A)) ≠ $D")
Base.:*(a::Quantity{D,U},b::ConvertUnit{A,S,U}) where {D,U,A,S} = (D==A && S==U) ? a : inv(A)==D ? a*ConvertUnit{inv(A),U,S}() : error("ConvertUnit incompatible $D ≠ $(inv(A))")
Base.:*(a::ConvertUnit{A,U,S},b::Quantity{D,U}) where {A,D,U,S} = A==D ? Quantity{D,S}(ratio(D,U,S)*b.v) : error("ConvertUnit incompatible $A ≠ $B")
Base.:*(a::Quantity{D,U},b::ConvertUnit{A,U,S}) where {A,D,U,S} = A==D ? Quantity{D,S}(a.v*ratio(D,U,S)) : error("ConvertUnit incompatible $A ≠ $B")

Base.:(==)(a::Quantity{A,U},b::Quantity{B,U}) where {A,B,U} = U(A) == U(B) && a.v == b.v
Base.:(==)(a::Number,b::Quantity{D}) where D = iszero(norm(first(value(D),dims-1))) && a == b.v*10^last(value(D))
Base.:(==)(a::Quantity{D},b::Number) where D = iszero(norm(first(value(D),dims-1))) && b == a.v*10^last(value(D))

@pure isunknown(x) = false

@pure sumabs(::Dimension{D}) where D = sum(abs.(D.v))
@pure add(d::Dimension{D},::Dimension{D},U) where D = islog(D) ? d+d : d
@pure function add(a::Dimension{A},b::Dimension{B},U) where {A,B}
    Ua,Ub = U(a),U(b)
    if Ua == Ub
        sa,sb = sumabs(a),sumabs(b)
        sb < sa || (sb == sa && Ub == b) ? b : a
    elseif islog(A) && islog(B)
        a + b
    else
        throw(error("$(Ua) ≠ $(Ub)"))
    end
end
@pure sub(d::Dimension{D},::Dimension{D},U) where D = islog(D) ? 𝟙 : d
@pure function sub(a::Dimension{A},b::Dimension{B},U) where {A,B}
    islog(A) && islog(B) ? a - b : add(a,b)
end

Base.:+(a::Quantity{A,U},b::Quantity{B,U}) where {A,B,U} = Quantity{add(A,B,U),U}(a.v+b.v)
Base.:-(a::Quantity{A,U},b::Quantity{B,U}) where {A,B,U} = Quantity{sub(A,B,U),U}(a.v-b.v)

# other

Base.isone(x::Dimension{D}) where D = isone(D)
Base.isone(x::Quantity{D}) where D = false
Base.iszero(x::Dimension{D}) where D = iszero(D)

# Quantities

struct Quantities{D,U,N,T} <: TupleVector{N,T}
    v::Values{N,T}
    Quantities{D,U,N,T}(v::Values{N,T}) where {D,U,N,T} = new{D,normal(U),N,T}(v)
end

Quantities{D,U,N,T}(v::Tuple) where {D,U,N,T} = Quantities{D,U,N,T}(Values{N,T}(v))
Quantities{D,U}(v::Values{N,T}) where {D,U,N,T} = Quantities{D,U,N,T}(v)
Quantities(D,U,v) = Quantities{D,U}(v)

(S::UnitSystem)(x::Quantities{D,U}) where {D,U} = Quantities{D,S}(D.(x.v,Ref(U),Ref(S)))

export Quantity, Quantities

Base.getindex(x::Quantities{D,U},i::Integer) where {D,U} = Quantity{D,U}(x.v[i])
Base.getindex(x::Quantities{D,U},i::Int) where {D,U} = Quantity{D,U}(x.v[i])

# functors
# 1,2,3,4, 5, 6, 7,  8,9,10,11
#kB,ħ,𝘤,μ₀,mₑ,Mᵤ,Kcd,A,λ,αL,g₀
# F,M,L,T, Q, Θ, N,  J,A,R, C

import UnitSystems: UnitSystem
UnitSystem(::Dimension{D}) where D = Dimension{UnitSystem(D)}()
UnitSystem(g::LogGroup{B}) where B = LogGroup{B}(UnitSystem(value(g)))
UnitSystem(g::ExpGroup{B}) where B = ExpGroup{B}(UnitSystem(value(g)))
function UnitSystem(d::Group)
    Group(Values(
        -d.v[6],
        d.v[3]+d.v[4]+d.v[5]/2-d.v[1]-d.v[8],
        3(d.v[1])+2d.v[6]+4d.v[8]-d.v[3]-2d.v[4]-d.v[5]/2,
        d.v[5]/-2,
        d.v[2]+d.v[6]+d.v[7]+2(d.v[1]+d.v[8])-d.v[3]-d.v[4],
        -d.v[7],
        d.v[8],
        d.v[3]+d.v[4]+d.v[9]+d.v[5]/2-d.v[1]-d.v[8],
        d.v[10]-d.v[5]/2,
        -(d.v[5]+d.v[11]),
        d.v[3]+d.v[4]-d.v[6]-2(d.v[1]+d.v[8])))
end
function UnitSystem(d::Group{<:Integer})
    Group(Values(
        -d.v[6],
        d.v[3]+d.v[4]+d.v[5]//2-d.v[1]-d.v[8],
        3(d.v[1])+2d.v[6]+4d.v[8]-d.v[3]-2d.v[4]-d.v[5]//2,
        d.v[5]//-2,
        d.v[2]+d.v[6]+d.v[7]+2(d.v[1]+d.v[8])-d.v[3]-d.v[4],
        -d.v[7],
        d.v[8],
        d.v[3]+d.v[4]+d.v[9]+d.v[5]//2-d.v[1]-d.v[8],
        d.v[10]-d.v[5]//2,
        -(d.v[5]+d.v[11]),
        d.v[3]+d.v[4]-d.v[6]-2(d.v[1]+d.v[8])))
end

Quantity(::ConvertUnit{D,U,S}) where {D,U,S} = Quantity{D,S}(ratio(D,U,S))
(u::UnitSystem)(::ConvertUnit{D,U,S}) where {D,U,S} = Quantity{D,u}(ratio(D,U,S))
(u::UnitSystem)(a::Number, d::Dimension) = Quantity{d,u}(a)
(u::UnitSystem)(a::Number, d::AbelianGroup) = Quantity{Dimension(d),u}(a)
(u::UnitSystem)(::Dimension{D}) where D = Dimension{normal(u)(D)}()
(u::UnitSystem)(d::LogGroup{B}) where B = LogGroup{B}(u(d.v))
(u::UnitSystem)(d::ExpGroup{B}) where B = ExpGroup{B}(u(d.v))


text:
################################################################################
#
#   group_elements.jl : Interface for group elements
#
################################################################################
# Obligatory methods
################################################################################

@doc Markdown.doc"""
    parent(g::GroupElement)

Return the parent object of the group element.
"""
Base.parent(g::GroupElement) =
    throw(InterfaceNotImplemented(:Group, "Base.parent(::$(typeof(g)))"))

@doc Markdown.doc"""
    parent_type(::Type{<:GroupElement})

Return the type of parent object, based on the type of its elements.
"""
parent_type(::Type{GEl}) where {GEl <: GroupElement} =
    throw(InterfaceNotImplemented(
        :Group,
        "GroupsCore.parent_type(::Type{$GEl})"
       ))
parent_type(g::GroupElement) = parent_type(typeof(g))

@doc Markdown.doc"""
    ==(g::GEl, h::GEl) where {GEl <: GroupElement}

Return `true` if and only if the mathematical equality $g = h$ holds.

!!! note

    This function may not return, due to unsolvable word problem.
"""
Base.:(==)(g::GEl, h::GEl) where {GEl <: GroupElement} = throw(
    InterfaceNotImplemented(:Group, "Base.:(==)(::$GEl, ::$GEl)"),
)

Base.deepcopy_internal(g::GroupElement, stackdict::IdDict) = throw(
    InterfaceNotImplemented(
        :Group,
        "Base.deepcopy_internal(::$(typeof(g)), ::IdDict)",
    ),
)
# TODO: Technically, it is not necessary to implement `deepcopy_internal` method
# if `parent(g)` can be reconstructed exactly from `g` (i.e. either it's cached,
# or a singleton). However by defining this fallback we force everybody to
# implement it, except isbits group elements.

@doc Markdown.doc"""
    inv(g::GroupElement)

Return $g^{-1}$, the group inverse.
"""
Base.inv(g::GroupElement) =
    throw(InterfaceNotImplemented(:Group, "Base.inv(::$(typeof(g)))"))

@doc Markdown.doc"""
    *(g::GEl, h::GEl) where {GEl <: GroupElement}

Return $g h$, the result of group binary operation.
"""
Base.:(*)(g::GEl, h::GEl) where {GEl <: GroupElement} = throw(
    InterfaceNotImplemented(
        :Group,
        "Base.:(*)(::$(typeof(g)), ::$(typeof(g)))",
    ),
)

@doc Markdown.doc"""
    isfiniteorder(g::GroupElement)

Return `true` if $g$ is of finite order, possibly without computing it.
"""
function isfiniteorder(g::GroupElement)
    isfinite(parent(g)) && return true
    throw(
        InterfaceNotImplemented(:Group, "GroupsCore.isfiniteorder(::$(typeof(g)))"),
    )
end

################################################################################
# Default implementations
################################################################################

@doc Markdown.doc"""
    one(g::GroupElement)

Return the identity element in the group of $g$.
"""
Base.one(g::GroupElement) = one(parent(g))

@doc Markdown.doc"""
    order(::Type{T} = BigInt, g::GroupElement) where T

Return the order of $g$ as an instance of `I`. If $g$ is of infinite order
`GroupsCore.InfiniteOrder` exception will be thrown.

!!! warning

    Only arbitrary sized integers are required to return a mathematically
    correct answer.
"""
function order(::Type{T}, g::GroupElement) where T
    isfiniteorder(g) || throw(InfiniteOrder(g))
    isone(g) && return T(1)
    o = T(1)
    gg = deepcopy(g)
    out = similar(g)
    while !isone(gg)
        o += T(1)
        gg = mul!(out, gg, g)
    end
    return o
end
order(g::GroupElement) = order(BigInt, g)

@doc Markdown.doc"""
    conj(g::GEl, h::GEl) where {GEl <: GroupElement}

Return conjugation of $g$ by $h$, i.e. $h^{-1} g h$.
"""
Base.conj(g::GEl, h::GEl) where {GEl <: GroupElement} = conj!(similar(g), g, h)

@doc Markdown.doc"""
    ^(g::GEl, h::GEl) where {GEl <: GroupElement}

Alias for [`conj`](@ref GroupsCore.conj).
"""
Base.:(^)(g::GEl, h::GEl) where {GEl <: GroupElement} = conj(g, h)

@doc Markdown.doc"""
    commutator(g::GEl, h::GEl, k::GEl...) where {GEl <: GroupElement}

Return the left associative iterated commutator $[[g, h], ...]$, where
$[g, h] = g^{-1} h^{-1} g h$.
"""
function commutator(g::GEl, h::GEl, k::GEl...) where {GEl <: GroupElement}
    res = commutator!(similar(g), g, h)
    for l in k
        res = commutator!(res, res, l)
    end
    return res
end

Base.literal_pow(::typeof(^), g::GroupElement, ::Val{-1}) = inv(g)

Base.:(/)(g::GEl, h::GEl) where {GEl <: GroupElement} =
    div_right!(similar(g), g, h)

################################################################################
# Default implementations that (might) need performance modification
################################################################################

@doc Markdown.doc"""
    similar(g::GroupElement)

Return a group element sharing the parent with $g$. Might be arbitrary and
possibly uninitialized.
"""
Base.similar(g::GroupElement) = one(g)

@doc Markdown.doc"""
    isone(g::GroupElement)

Return true if $g$ is the identity element.
"""
Base.isone(g::GroupElement) = g == one(g)

@doc Markdown.doc"""
    isequal(g::GEl, h::GEl) where {GEl <: GroupElement}

Return the "best effort" equality for group elements.

The implication `isequal(g, h)` → $g = h$, must be always satisfied, i.e.
`isequal(g, h)` might return `false` even if `g == h` holds (i.e. $g$ and $h$
are mathematically equal).

For example, in a finitely presented group, `isequal` may return the equality
of words.
"""
Base.isequal(g::GEl, h::GEl) where {GEl <: GroupElement} = g == h

function Base.:(^)(g::GroupElement, n::Integer)
    n == 0 && return one(g)
    n < 0 && return inv(g)^-n
    return Base.power_by_squaring(g, n)
end

# NOTE: Modification RECOMMENDED for performance reasons
Base.hash(g::GroupElement, h::UInt) = hash(typeof(g), h)

################################################################################
# Mutable API where modifications are recommended for performance reasons
################################################################################

@doc Markdown.doc"""
    one!(g::GroupElement)

Return `one(g)`, possibly modifying `g`.
"""
one!(g::GroupElement) = one(parent(g))

@doc Markdown.doc"""
    inv!(out::GEl, g::GEl) where {GEl <: GroupElement}

Return `inv(g)`, possibly modifying `out`. Aliasing of `g` with `out` is
allowed.
"""
inv!(out::GEl, g::GEl) where {GEl <: GroupElement} = inv(g)

@doc Markdown.doc"""
    mul!(out::GEl, g::GEl, h::GEl) where {GEl <: GroupElement}

Return $g h$, possibly modifying `out`. Aliasing of `g` or `h` with `out` is
allowed.
"""
mul!(out::GEl, g::GEl, h::GEl) where {GEl <: GroupElement} = g * h

@doc Markdown.doc"""
    div_right!(out::GEl, g::GEl, h::GEl) where {GEl <: GroupElement}

Return $g h^{-1}$, possibly modifying `out`. Aliasing of `g` or `h` with `out`
is allowed.
"""
div_right!(out::GEl, g::GEl, h::GEl) where {GEl <: GroupElement} =
    mul!(out, g, inv(h))

@doc Markdown.doc"""
    div_left!(out::GEl, g::GEl, h::GEl) where {GEl <: GroupElement}

Return $h^{-1} g$, possibly modifying `out`. Aliasing of `g` or `h` with `out`
is allowed.
"""
function div_left!(out::GEl, g::GEl, h::GEl) where {GEl <: GroupElement}
    out = (out === g || out === h) ? inv(h) : inv!(out, h)
    return mul!(out, out, g)
end

@doc Markdown.doc"""
    conj!(out::GEl, g::GEl, h::GEl) where {GEl <: GroupElement}

Return $h^{-1} g h$, possibly modifying `out`. Aliasing of `g` or `h` with
`out` is allowed.
"""
function conj!(out::GEl, g::GEl, h::GEl) where {GEl <: GroupElement}
    out = (out === g || out === h) ? inv(h) : inv!(out, h)
    out = mul!(out, out, g)
    return mul!(out, out, h)
end

@doc Markdown.doc"""
    commutator!(out::GEl, g::GEl, h::GEl) where {GEl <: GroupElement}

Return $g^{-1} h^{-1} g h$, possibly modifying `out`. Aliasing of `g` or `h`
with `out` is allowed.
"""
function commutator!(out::GEl, g::GEl, h::GEl) where {GEl <: GroupElement}
    # TODO: can we make commutator! with 3 arguments without allocation??
    out = conj!(out, g, h)
    return div_left!(out, out, g)
end


text:
using CompScienceMeshes
using BEAST
using Base.Test

fn = joinpath(dirname(@__FILE__),"assets","rect1.in")

mesh = readmesh(fn)
edgs = skeleton(mesh,1)

charts = [chart(edgs,edg) for edg in cells(edgs)]
ctrs = [cartesian(center(cht)) for cht in charts]

ND = BEAST.nedelec(mesh, edgs)

@test numfunctions(ND) == numcells(edgs)
@test length(ND.fns[1]) == 1
@test length(ND.fns[2]) == 1
@test length(ND.fns[3]) == 2
@test length(ND.fns[4]) == 1
@test length(ND.fns[5]) == 1

# center of the diagonal
ctr = center(charts[3])
face_charts = [chart(mesh,fce) for fce in cells(mesh)]
nbd1 = neighborhood(face_charts[1], carttobary(face_charts[1], ctr))
nbd2 = neighborhood(face_charts[2], carttobary(face_charts[2], ctr))
t = tangents(ctr,1)
ut = t / norm(t)

ndlocal = refspace(ND)
fn = ND.fns[3]
v1 = dot(fn[1].coeff*ndlocal(nbd1)[1][1],ut)
v2 = dot(fn[2].coeff*ndlocal(nbd2)[2][1],ut)
@test v1 ≈ v2 ≈ -1/√2 ≈ -1/volume(charts[3])


text:
using JSON
using JLD
using StatsBase
using BenchmarkTools
using Glob
using HTTP
# using Libz
# using Base.LibGit2

# push!(LOAD_PATH, @__DIR__)
# using Check

const cachefolder = joinpath(@__DIR__, "dlcache")
const repobase = "https://raw.githubusercontent.com/JuliaCI/BaseBenchmarkReports/master"

function summarize(benchmarks)
    if isdefined(benchmarks, :data)
        return map(benchmarks.data) do entry
            string(entry.first) => summarize(entry.second)
        end
    elseif isdefined(benchmarks, :times)
        return summarystats(benchmarks.times)
    else
        return nothing
    end
end

datedirname(date) = string("daily_", Dates.year(date), "_", Dates.month(date), "_", Dates.day(date))

function get_daily_data(date)
    dailydir = joinpath(cachefolder, date)
    isdir(dailydir) || mkdir(dailydir)
    filegz = joinpath(dailydir, "data.tar.gz")
    if !isfile(filegz)
        uri = string("$repobase/$date/data.tar.gz")
        datagz = HTTP.get(uri)
        status = HTTP.status(datagz)
        if status != 200
            warn("$uri => $status")
            return nothing # assume there is no data for this day
        end
        write(filegz, HTTP.bytes(datagz))
    end
    run(`tar -xzf $filegz -C $dailydir`)
    files = readdir(glob"data/*_primary.jld", dailydir)
    if isempty(files)
        data = nothing # no data?
    else
        length(files) == 1 || warn("Too many primary.jld files for $date; taking one.")
        data = BenchmarkTools.load(files[1], "results")
    end
    rm(joinpath(dailydir, "data"), recursive = true)
    return data
end

function get_daily_summary(date)
    file = joinpath(cachefolder, "$date.json")
    if !isfile(file)
        data = summarize(get_daily_data(date))
        open(file, "w") do io
            JSON.print(io, data)
        end
    end
    return JSON.parse(readstring(file))
end

# take a Vector{Dict{String, Float64}} and turn it into Dict{String, Vector{Float64}}
# but with nesting and nullable
function transpose_daily_summary(data)
    groups = Set{String}()
    for day in data
        isa(day, Dict) && union!(groups, keys(day))
    end

    if isempty(groups)
        output = Float64[ (day === nothing ? NaN : day) for day in data ]
    else
        output = Dict{String, Any}()
        sizehint!(output, length(groups))
        for group in groups
            values = Any[ (isa(day, Dict) ? get(day, group, nothing) : day) for day in data ]
            output[group] = transpose_daily_summary(values)
        end
    end
    return output
end

# take a Dict{String, Dict{String, ...{String, Vector{Float64}}}
# and flatten it into Dict{String, Vector{Float64}} by name-mangling the paths together
flatten_summary(data) = flatten_summary("", data)
flatten_summary(prefix::String, data::Vector{Float64}) = Dict(prefix => data)
function flatten_summary(prefix::String, data::Dict)
    output = Dict{String, Vector{Float64}}()
    for (k, v) in data
        merge!(output, flatten_summary(isempty(prefix) ? k : "$prefix.$k", v))
    end
    return output
end

if false
    today = trunc(now(), Dates.Day)
    until = today - Dates.Day(90)
    headers = String[]
    data = Any[]
    while today > until
        date = datedirname(today)
        push!(headers, date)
        push!(data, get_daily_summary(date))
        today -= Dates.Day(1)
    end
    missing = data .== nothing
    data = transpose_daily_summary(data)
    data = flatten_summary(data)
end

max_stars_repo_name: vtjnash/PerfPlots

max_issues_repo_name: vtjnash/PerfPlots

max_forks_repo_name: vtjnash/PerfPlots


text:
# Maximum feasibility of 0.4769921436588103 for [1,3,5,8]
# Mean feasibility of 0.3221046443268665

include("CenterOfMass.jl")
include("setup_parameters.jl")

mass_resolution = 0.1
masses = collect(min_mass:mass_resolution:max_mass)
resolution = 0.1

max_robots = 4
actuator_limit = 6.0

#########################
# Myopic feasibility rate
#########################
boundary_vec = flatten(circle_ps)
min_p = minimum(boundary_vec, 2)[1:2]
max_p = maximum(boundary_vec, 2)[1:2]
position_ranges = map(x->x[1]:resolution:x[2], zip(min_p, max_p))
ranges = (position_ranges..., masses)

robots_feasible(robots, theta) =
  lifting_feasibility(attachment_ws[robots],
                      theta[3]*g*point_to_wrench(theta[1:2]),
                      actuator_limit)

maximum_feasibility_rate = 0
maximum_feasibility_combination = zeros(Int64,0)

good_inds = zeros(Bool,map(length,ranges))
for (ii,x) = enumerate(ranges[1])
  for (jj,y) = enumerate(ranges[2])
    for (kk,m) = enumerate(ranges[3])
      Convex.clearmemory()
      if interior_q([x;y])
        gravity_w = m*g*point_to_wrench([x;y])
        feasible = check_feasible_configuration(Array{Array{Float64}}(0), attachment_ws,
          gravity_w, actuator_limit, max_robots)

        if feasible
          good_inds[ii,jj,kk] = true
        end
      end
    end
  end
end
num_good = sum(good_inds)

mean_feasibility = 0.0
cs = combinations(1:length(attachment_ws), max_robots)
for (ii, c) = enumerate(cs)
  num_feasible = 0
  for (jj,x) = enumerate(ranges[1])
    for (kk,y) = enumerate(ranges[2])
      for (ll,m) = enumerate(ranges[3])
        Convex.clearmemory()
        theta = [x,y,m]
        if good_inds[jj,kk,ll]
          num_feasible += robots_feasible(c, theta)
        end
      end
    end
  end
  feasibility = num_feasible / num_good
  mean_feasibility += feasibility
  if feasibility > maximum_feasibility_rate
    maximum_feasibility_combination = c
    maximum_feasibility_rate = feasibility
    println("New best: $(feasibility)")
  end
  println("Finished $(ii) of $(length(cs)), best: $(maximum_feasibility_rate)")
end
mean_feasibility /= length(cs)

println("Maximum feasibility of $(maximum_feasibility_rate) for $(maximum_feasibility_combination)")
println("Mean feasibility of $(mean_feasibility)")

max_stars_repo_name: mcorah/CenterOfMass

max_issues_repo_name: mcorah/CenterOfMass

max_forks_repo_name: mcorah/CenterOfMass


text:
using OrdinaryDiffEq
using Trixi

###############################################################################
# semidiscretization of the compressible Euler multicomponent equations

# 1) Air  2) Helium  3) Neon  4) Argon  5) Hydrogen  6) Methane  7) Butan
equations           = CompressibleEulerMulticomponentEquations2D(gammas        = (1.4, 1.667, 1.667, 1.667, 1.4, 1.3, 1.091),
                                                                 gas_constants = (0.287, 2.0769, 0.4119, 0.2081, 4.124, 0.5182, 0.1433))

initial_condition   = initial_condition_shock_bubble

surface_flux        = flux_lax_friedrichs
volume_flux         = flux_chandrashekar
basis               = LobattoLegendreBasis(3)
indicator_sc        = IndicatorHennemannGassner(equations, basis,
                                                alpha_max=0.6,
                                                alpha_min=0.001,
                                                alpha_smooth=true,
                                                variable=density_pressure)
volume_integral     = VolumeIntegralShockCapturingHG(indicator_sc;
                                                     volume_flux_dg=volume_flux,
                                                     volume_flux_fv=surface_flux)
solver              = DGSEM(basis, surface_flux, volume_integral)

coordinates_min     = (-3.0, -3.0)
coordinates_max     = ( 3.0,  3.0)
mesh                = TreeMesh(coordinates_min, coordinates_max,
                               initial_refinement_level=3,
                               n_cells_max=1_000_000)

semi                = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)

###############################################################################
# ODE solvers, callbacks etc.

tspan               = (0.0, 0.01)
ode                 = semidiscretize(semi, tspan)

summary_callback    = SummaryCallback()

analysis_interval   = 300 # change to 200+ for whole calc.
analysis_callback   = AnalysisCallback(semi, interval=analysis_interval,
                                       extra_analysis_integrals=(Trixi.density,))

alive_callback      = AliveCallback(analysis_interval=analysis_interval)

save_solution       = SaveSolutionCallback(interval=300,     # 40 or change to 200+ for whole calc.
                                           save_initial_solution=true,
                                           save_final_solution=true,
                                           solution_variables=cons2prim)

stepsize_callback   = StepsizeCallback(cfl=0.3)

callbacks           = CallbackSet(summary_callback,
                                  analysis_callback,
                                  alive_callback,
                                  save_solution,
                                  stepsize_callback)

###############################################################################
# run the simulation
sol                 = solve(ode, CarpenterKennedy2N54(williamson_condition=false),
                                                      dt=1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
                                                      save_everystep=false,
                                                      callback=callbacks,
                                                      maxiters=1e5);
summary_callback() # print the timer summary


text:
export CaesarConfig, VisualizationConfig

struct VisualizationConfig
    drawDepth::Bool
    cleanSlamInDb::Bool
    drawEdges::Bool
end

struct CaesarConfig
    numParticles::Int64
    multiSession::Vector{AbstractString}
    visualizationConfig::VisualizationConfig
end


text:
using Documenter
using SimpleFWA

makedocs(
    sitename = "SimpleFWA",
    modules = [SimpleFWA]
)

deploydocs(
    repo = "git@github.com:hondoRandale/SimpleFWA.jl.git",
)

deploydocs(
    repo   = "github.com/USER_NAME/PACKAGE_NAME.jl.git",
    deps   = Deps.pip("mkdocs", "pygments", "python-markdown-math"),
    make   = () -> run(`mkdocs build`),
    target = "site"
)


text:
using Test
using QuantumOptics
using FFTW, LinearAlgebra, Random

@testset "particle" begin

Random.seed!(0)

D(op1::AbstractOperator, op2::AbstractOperator) = abs(tracedistance_nh(dense(op1), dense(op2)))

N = 200
xmin = -32.5
xmax = 20.1

basis_position = PositionBasis(xmin, xmax, N)
basis_momentum = MomentumBasis(basis_position)

b2 = PositionBasis(basis_momentum)
@test basis_position.xmax - basis_position.xmin ≈ b2.xmax - b2.xmin
@test basis_position.N == b2.N

# Test Gaussian wave-packet in both bases
x0 = 5.1
p0 = -3.2
sigma = 1.
sigma_x = sigma/sqrt(2)
sigma_p = 1.0/(sigma*sqrt(2))

psi0_bx = gaussianstate(basis_position, x0, p0, sigma)
psi0_bp = gaussianstate(basis_momentum, x0, p0, sigma)

@test 1 ≈ norm(psi0_bx)
@test 1 ≈ norm(psi0_bp)

p_bx = momentum(basis_position)
x_bx = position(basis_position)

@test 1e-10 > D(p_bx, transform(basis_position, basis_momentum)*dense(momentum(basis_momentum))*transform(basis_momentum, basis_position))

p_bp = momentum(basis_momentum)
x_bp = position(basis_momentum)

@test x0 ≈ expect(x_bx, psi0_bx)
@test 0.1 > abs(p0 - expect(p_bx, psi0_bx))
@test p0 ≈ expect(p_bp, psi0_bp)
@test 0.1 > abs(x0 - expect(x_bp, psi0_bp))

@test 1e-13 > abs(variance(x_bx, psi0_bx) - sigma^2/2)
@test 1e-13 > abs(variance(x_bp, psi0_bp) - sigma^2/2)
@test 1e-13 > abs(variance(p_bx, psi0_bx) - 1/(2*sigma^2))
@test 1e-13 > abs(variance(p_bp, psi0_bp) - 1/(2*sigma^2))

# Test potentialoperator
V(x) = x^2
V_bx = potentialoperator(basis_position, V)
V_bp = potentialoperator(basis_momentum, V)

@test expect(V_bp, psi0_bp) ≈ expect(V_bx, psi0_bx)

# Test FFT transformation
function transformation(b1::MomentumBasis, b2::PositionBasis, psi::Ket)
    Lp = (b1.pmax - b1.pmin)
    dx = particle.spacing(b2)
    if b1.N != b2.N || abs(2*pi/dx - Lp)/Lp > 1e-12
        throw(IncompatibleBases())
    end
    N = b1.N
    psi_shifted = exp.(1im*b2.xmin*(particle.samplepoints(b1) .- b1.pmin)).*psi.data
    psi_fft = exp.(1im*b1.pmin*particle.samplepoints(b2)).*ifft(psi_shifted)*sqrt(N)
    return Ket(b2, psi_fft)
end

function transformation(b1::PositionBasis, b2::MomentumBasis, psi::Ket)
    Lx = (b1.xmax - b1.xmin)
    dp = particle.spacing(b2)
    if b1.N != b2.N || abs(2*pi/dp - Lx)/Lx > 1e-12
        throw(IncompatibleBases())
    end
    N = b1.N
    psi_shifted = exp.(-1im*b2.pmin*(particle.samplepoints(b1) .- b1.xmin)).*psi.data
    psi_fft = exp.(-1im*b1.xmin*particle.samplepoints(b2)).*fft(psi_shifted)/sqrt(N)
    return Ket(b2, psi_fft)
end

psi0_bx_fft = transformation(basis_position, basis_momentum, psi0_bx)
psi0_bp_fft = transformation(basis_momentum, basis_position, psi0_bp)

@test 0.1 > abs(x0 - expect(x_bp, psi0_bx_fft))
@test p0 ≈ expect(p_bp, psi0_bx_fft)
@test 0.1 > abs(p0 - expect(p_bx, psi0_bp_fft))
@test x0 ≈ expect(x_bx, psi0_bp_fft)

@test 1e-12 > norm(psi0_bx_fft - psi0_bp)
@test 1e-12 > norm(psi0_bp_fft - psi0_bx)

Tpx = transform(basis_momentum, basis_position)
Txp = transform(basis_position, basis_momentum)
psi0_bx_fft = Tpx*psi0_bx
psi0_bp_fft = Txp*psi0_bp

@test 0.1 > abs(x0 - expect(x_bp, psi0_bx_fft))
@test p0 ≈ expect(p_bp, psi0_bx_fft)
@test 0.1 > abs(p0 - expect(p_bx, psi0_bp_fft))
@test x0 ≈ expect(x_bx, psi0_bp_fft)

@test 1e-12 > norm(psi0_bx_fft - psi0_bp)
@test 1e-12 > norm(psi0_bp_fft - psi0_bx)

@test 1e-12 > norm(dagger(Tpx*psi0_bx) - dagger(psi0_bx)*dagger(Tpx))
@test 1e-12 > norm(dagger(Txp*psi0_bp) - dagger(psi0_bp)*dagger(Txp))

# Test gemv!
psi_ = deepcopy(psi0_bp)
operators.gemv!(Complex(1.), Tpx, psi0_bx, Complex(0.), psi_)
@test 1e-12 > norm(psi_ - psi0_bp)
operators.gemv!(Complex(1.), Tpx, psi0_bx, Complex(1.), psi_)
@test 1e-12 > norm(psi_ - 2*psi0_bp)

psi_ = deepcopy(psi0_bx)
operators.gemv!(Complex(1.), Txp, psi0_bp, Complex(0.), psi_)
@test 1e-12 > norm(psi_ - psi0_bx)
operators.gemv!(Complex(1.), Txp, psi0_bp, Complex(1.), psi_)
@test 1e-12 > norm(psi_ - 2*psi0_bx)

alpha = complex(3.2)
beta = complex(-1.2)
randdata1 = rand(ComplexF64, N)
randdata2 = rand(ComplexF64, N)

state = Ket(basis_position, randdata1)
result_ = Ket(basis_momentum, copy(randdata2))
result0 = alpha*dense(Tpx)*state + beta*result_
operators.gemv!(alpha, Tpx, state, beta, result_)
@test 1e-11 > norm(result0 - result_)

state = Bra(basis_position, randdata1)
result_ = Bra(basis_momentum, copy(randdata2))
result0 = alpha*state*dense(Txp) + beta*result_
operators.gemv!(alpha, state, Txp, beta, result_)
@test 1e-11 > norm(result0 - result_)

state = Ket(basis_momentum, randdata1)
result_ = Ket(basis_position, copy(randdata2))
result0 = alpha*dense(Txp)*state + beta*result_
operators.gemv!(alpha, Txp, state, beta, result_)
@test 1e-11 > norm(result0 - result_)

state = Bra(basis_momentum, randdata1)
result_ = Bra(basis_position, copy(randdata2))
result0 = alpha*state*dense(Tpx) + beta*result_
operators.gemv!(alpha, state, Tpx, beta, result_)
@test 1e-11 > norm(result0 - result_)

# Test gemm!
rho0_xx = tensor(psi0_bx, dagger(psi0_bx))
rho0_xp = tensor(psi0_bx, dagger(psi0_bp))
rho0_px = tensor(psi0_bp, dagger(psi0_bx))
rho0_pp = tensor(psi0_bp, dagger(psi0_bp))

rho_ = DenseOperator(basis_momentum, basis_position)
operators.gemm!(Complex(1.), Tpx, rho0_xx, Complex(0.), rho_)
@test 1e-12 > D(rho_, rho0_px)
@test 1e-12 > D(Tpx*rho0_xx, rho0_px)

rho_ = DenseOperator(basis_position, basis_momentum)
operators.gemm!(Complex(1.), rho0_xx, Txp, Complex(0.), rho_)
@test 1e-12 > D(rho_, rho0_xp)
@test 1e-12 > D(rho0_xx*Txp, rho0_xp)

rho_ = DenseOperator(basis_momentum, basis_momentum)
operators.gemm!(Complex(1.), Tpx, rho0_xp, Complex(0.), rho_)
@test 1e-12 > D(rho_, rho0_pp)
@test 1e-12 > D(Tpx*rho0_xx*Txp, rho0_pp)

rho_ = DenseOperator(basis_momentum, basis_momentum)
operators.gemm!(Complex(1.), rho0_px, Txp, Complex(0.), rho_)
@test 1e-12 > D(rho_, rho0_pp)
@test 1e-12 > D(Txp*rho0_pp*Tpx, rho0_xx)

alpha = complex(3.2)
beta = complex(-1.2)
randdata1 = rand(ComplexF64, N, N)
randdata2 = rand(ComplexF64, N, N)

op = DenseOperator(basis_position, basis_position, randdata1)
result_ = DenseOperator(basis_momentum, basis_position, copy(randdata2))
result0 = alpha*dense(Tpx)*op + beta*result_
operators.gemm!(alpha, Tpx, op, beta, result_)
@test 1e-11 > D(result0, result_)

result_ = DenseOperator(basis_position, basis_momentum, copy(randdata2))
result0 = alpha*op*dense(Txp) + beta*result_
operators.gemm!(alpha, op, Txp, beta, result_)
@test 1e-11 > D(result0, result_)

op = DenseOperator(basis_momentum, basis_momentum, randdata1)
result_ = DenseOperator(basis_position, basis_momentum, copy(randdata2))
result0 = alpha*dense(Txp)*op + beta*result_
operators.gemm!(alpha, Txp, op, beta, result_)
@test 1e-11 > D(result0, result_)

result_ = DenseOperator(basis_momentum, basis_position, copy(randdata2))
result0 = alpha*op*dense(Tpx) + beta*result_
operators.gemm!(alpha, op, Tpx, beta, result_)
@test 1e-11 > D(result0, result_)

# Test FFT with lazy product
psi_ = deepcopy(psi0_bx)
operators.gemv!(Complex(1.), LazyProduct(Txp, Tpx), psi0_bx, Complex(0.), psi_)
@test 1e-12 > norm(psi_ - psi0_bx)
@test 1e-12 > norm(Txp*(Tpx*psi0_bx) - psi0_bx)

psi_ = deepcopy(psi0_bx)
id = dense(identityoperator(basis_momentum))
operators.gemv!(Complex(1.), LazyProduct(Txp, id, Tpx), psi0_bx, Complex(0.), psi_)
@test 1e-12 > norm(psi_ - psi0_bx)
@test 1e-12 > norm(Txp*id*(Tpx*psi0_bx) - psi0_bx)

# Test dense FFT operator
Txp_dense = DenseOperator(Txp)
Tpx_dense = DenseOperator(Tpx)
@test isa(Txp_dense, DenseOperator)
@test isa(Tpx_dense, DenseOperator)
@test 1e-5 > D(Txp_dense*rho0_pp*Tpx_dense, rho0_xx)

# Test FFT in 2D
N = [21, 18]
xmin = [-32.5, -10π]
xmax = [24.1, 9π]

basis_position = [PositionBasis(xmin[i], xmax[i], N[i]) for i=1:2]
basis_momentum = MomentumBasis.(basis_position)

x0 = [5.1, -0.2]
p0 = [-3.2, 1.33]
sigma = [1., 0.9]
sigma_x = sigma./sqrt(2)
sigma_p = 1.0/(sigma.*sqrt(2))

Txp = transform(tensor(basis_position...), tensor(basis_momentum...))
Tpx = transform(tensor(basis_momentum...), tensor(basis_position...))

Txp_sub = [transform(basis_position[i], basis_momentum[i]) for i=1:2]
Tpx_sub = dagger.(Txp_sub)
Txp_dense = dense.(Txp_sub)
Txp_comp = tensor(Txp_dense...)

difference = (dense(Txp) - Txp_comp).data
@test isapprox(difference, zero(difference); atol=1e-12)

psi0_x = [gaussianstate(basis_position[i], x0[i], p0[i], sigma_x[i]) for i=1:2]
psi0_p = [gaussianstate(basis_momentum[i], x0[i], p0[i], sigma_p[i]) for i=1:2]

psi_p_fft = Tpx*tensor(psi0_x...)
psi_p_fft2 = tensor((Tpx_sub.*psi0_x)...)
@test norm(psi_p_fft - psi_p_fft2) < 1e-15

psi_x_fft = Txp*tensor(psi0_p...)
psi_x_fft2 = tensor((Txp_sub.*psi0_p)...)
@test norm(psi_p_fft - psi_p_fft2) < 1e-15

psi_p_fft = dagger(tensor(psi0_x...))*Txp
psi_p_fft2 = tensor((dagger.(psi0_x).*Txp_sub)...)
@test norm(psi_p_fft - psi_p_fft2) < 1e-15

psi_x_fft = dagger(tensor(psi0_p...))*Tpx
psi_x_fft2 = tensor((dagger.(psi0_p).*Tpx_sub)...)
@test norm(psi_p_fft - psi_p_fft2) < 1e-15

difference = (dense(Txp) - identityoperator(DenseOperator, Txp.basis_l)*Txp).data
@test isapprox(difference, zero(difference); atol=1e-12)
@test_throws AssertionError transform(tensor(basis_position...), tensor(basis_position...))
@test_throws particle.IncompatibleBases transform(SpinBasis(1//2)^2, SpinBasis(1//2)^2)

@test dense(Txp) == dense(Txp_sub[1] ⊗ Txp_sub[2])

# Test ket only FFTs
Txp = transform(tensor(basis_position...), tensor(basis_momentum...); ket_only=true)
Tpx = transform(tensor(basis_momentum...), tensor(basis_position...); ket_only=true)

Txp_sub = [transform(basis_position[i], basis_momentum[i]; ket_only=true) for i=1:2]
Tpx_sub = dagger.(Txp_sub)

psi_p_fft = Tpx*tensor(psi0_x...)
psi_p_fft2 = tensor((Tpx_sub.*psi0_x)...)
@test norm(psi_p_fft - psi_p_fft2) < 1e-15

psi_x_fft = Txp*tensor(psi0_p...)
psi_x_fft2 = tensor((Txp_sub.*psi0_p)...)
@test norm(psi_p_fft - psi_p_fft2) < 1e-15

psi_p_fft = dagger(tensor(psi0_x...))*Txp
psi_p_fft2 = tensor((dagger.(psi0_x).*Txp_sub)...)
@test norm(psi_p_fft - psi_p_fft2) < 1e-15

psi_x_fft = dagger(tensor(psi0_p...))*Tpx
psi_x_fft2 = tensor((dagger.(psi0_p).*Tpx_sub)...)
@test norm(psi_p_fft - psi_p_fft2) < 1e-15

psi_x_fft = Txp*tensor(psi0_p...)
psi_x_fft2 = tensor(Txp_sub...)*tensor(psi0_p...)
@test norm(psi_x_fft - psi_x_fft2) < 1e-15

# Test composite basis of mixed type
bc = FockBasis(2)
psi_fock = fockstate(FockBasis(2), 1)
psi1 = tensor(psi0_p[1], psi_fock, psi0_p[2])
psi2 = tensor(psi0_x[1], psi_fock, psi0_x[2])

basis_l = tensor(basis_position[1], bc, basis_position[2])
basis_r = tensor(basis_momentum[1], bc, basis_momentum[2])
Txp = transform(basis_l, basis_r; ket_only=true)
Tpx = transform(basis_r, basis_l; ket_only=true)

psi1_fft = Txp*psi1
psi1_fft2 = tensor(Txp_sub[1]*psi0_p[1], psi_fock, Txp_sub[2]*psi0_p[2])
@test norm(psi1_fft - psi1_fft2) < 1e-15

psi2_fft = Tpx*psi2
psi2_fft2 = tensor(Tpx_sub[1]*psi0_x[1], psi_fock, Tpx_sub[2]*psi0_x[2])
@test norm(psi2_fft - psi2_fft2) < 1e-15

Txp = transform(basis_l, basis_r)
Txp_sub = [transform(basis_position[i], basis_momentum[i]) for i=1:2]
difference = (dense(Txp) - tensor(dense(Txp_sub[1]), dense(one(bc)), dense(Txp_sub[2]))).data
@test isapprox(difference, zero(difference); atol=1e-12)

basis_l = tensor(bc, basis_position[1], basis_position[2])
basis_r = tensor(bc, basis_momentum[1], basis_momentum[2])
Txp2 = transform(basis_l, basis_r)
Tpx2 = transform(basis_r, basis_l)
difference = (dense(Txp) - permutesystems(dense(Txp2), [2, 1, 3])).data
@test isapprox(difference, zero(difference); atol=1e-13)
difference = (dense(dagger(Txp)) - permutesystems(dense(Tpx2), [2, 1, 3])).data
@test isapprox(difference, zero(difference); atol=1e-13)

# Test potentialoperator in more than 1D
N = [21, 18]
xmin = [-32.5, -10π]
xmax = [24.1, 9π]

basis_position = [PositionBasis(xmin[i], xmax[i], N[i]) for i=1:2]
basis_momentum = MomentumBasis.(basis_position)

bcomp_pos = tensor(basis_position...)
bcomp_mom = tensor(basis_momentum...)
V(x, y) = sin(x*y) + cos(x)
xsample, ysample = samplepoints.(basis_position)
V_op = diagonaloperator(bcomp_pos, [V(x, y) for y in ysample for x in xsample])
V_op2 = potentialoperator(bcomp_pos, V)
@test V_op == V_op2

basis_position = PositionBasis.(basis_momentum)
bcomp_pos = tensor(basis_position...)
Txp = transform(bcomp_pos, bcomp_mom)
Tpx = transform(bcomp_mom, bcomp_pos)
xsample, ysample = samplepoints.(basis_position)
V_op = Tpx*dense(diagonaloperator(bcomp_pos, [V(x, y) for y in ysample for x in xsample]))*Txp
V_op2 = potentialoperator(bcomp_mom, V)
@test V_op == V_op2

N = [17, 12, 9]
xmin = [-32.5, -10π, -0.1]
xmax = [24.1, 9π, 22.0]

basis_position = [PositionBasis(xmin[i], xmax[i], N[i]) for i=1:3]
basis_momentum = MomentumBasis.(basis_position)

bcomp_pos = tensor(basis_position...)
bcomp_mom = tensor(basis_momentum...)
V(x, y, z) = exp(-z^2) + sin(x*y) + cos(x)
xsample, ysample, zsample = samplepoints.(basis_position)
V_op = diagonaloperator(bcomp_pos, [V(x, y, z) for z in zsample for y in ysample for x in xsample])
V_op2 = potentialoperator(bcomp_pos, V)
@test V_op == V_op2

# Test error messages
b1 = PositionBasis(-1, 1, 50)
b2 = MomentumBasis(-1, 1, 30)
@test_throws bases.IncompatibleBases transform(b1, b2)
@test_throws bases.IncompatibleBases transform(b2, b1)

bc1 = b1 ⊗ bc
bc2 = b2 ⊗ bc
@test_throws bases.IncompatibleBases transform(bc1, bc2)
@test_throws bases.IncompatibleBases transform(bc2, bc1)

b1 = PositionBasis(-1, 1, 50)
b2 = MomentumBasis(-1, 1, 50)
bc1 = b1 ⊗ bc
bc2 = b2 ⊗ bc
@test_throws bases.IncompatibleBases transform(bc1, bc2)
@test_throws bases.IncompatibleBases transform(bc2, bc1)
@test_throws bases.IncompatibleBases transform(bc1, bc2; index=[2])

bc1 = b1 ⊗ b2
bc2 = b1 ⊗ b2
@test_throws bases.IncompatibleBases transform(bc1, bc2)
@test_throws bases.IncompatibleBases transform(bc2, bc1)

@test_throws bases.IncompatibleBases potentialoperator(bc ⊗ bc, V)

end # testset


text:
abstract type Visitor{T <: Integer} end

"""
    ncycles_n_i(n::Integer, i::Integer)

Compute the theoretical maximum number of cycles of size `i` in a directed graph of `n`
 vertices.
"""
ncycles_n_i(n::Integer, i::Integer) =
    binomial(big(n), big(n - i + 1)) * factorial(big(n - i))

"""
    maxsimplecycles(n::Integer)

Compute the theoretical maximum number of cycles in a directed graph of `n` vertices,
assuming there are no self-loops.

### References
- [Johnson](http://epubs.siam.org/doi/abs/10.1137/0204007).
"""
maxsimplecycles(n::Integer) = sum(x -> ncycles_n_i(n, x), 1:(n - 1))

"""
    maxsimplecycles(dg::::IsDirected, byscc::Bool = true)

Compute the theoretical maximum number of cycles in the directed graph `dg`.

The computation can be performed assuming the graph is complete or taking into account the
decomposition in strongly connected components (`byscc` parameter).

### Performance
A more efficient version is possible.

### References
- [Johnson](http://epubs.siam.org/doi/abs/10.1137/0204007)
"""
function maxsimplecycles end
@traitfn function maxsimplecycles(dg::::IsDirected, byscc::Bool=true)
    c = 0
    n = nv(dg)
    if !byscc
        c = maxsimplecycles(n)
    else
        for scc in strongly_connected_components(dg)
            if length(scc) > 1
                c += maxsimplecycles(length(scc))
            end
        end
    end
    return c
end

"""
```
type JohnsonVisitor{T<:Integer} <: Visitor{T}
    stack::Vector{T}
    blocked::BitArray
    blockedmap::Vector{Set{T}}
end

JohnsonVisitor(dg::::IsDirected)

```

Composite type that regroups the information needed for Johnson's algorithm.

`stack` is the stack of visited vertices. `blocked` is a boolean for each
vertex that tells whether it is blocked or not. `blockedmap` tells which
vertices to unblock if the key vertex is unblocked.

`JohnsonVisitor` may also be constructed directly from the directed graph.
"""
struct JohnsonVisitor{T <: Integer} <: Visitor{T}
    stack::Vector{T}
    blocked::BitArray
    blockedmap::Vector{Set{T}}
end

@traitfn function JohnsonVisitor(dg::::IsDirected)
    T = eltype(dg)
    return JohnsonVisitor(Vector{T}(), falses(nv(dg)), [Set{T}() for i in vertices(dg)])
end

"""
    unblock!{T<:Integer}(v::T, blocked::BitArray, B::Vector{Set{T}})

Unblock the vertices recursively.

`v` is the vertex to unblock, `blocked` tells whether a vertex is blocked or
not and `B` is the map that tells if the unblocking of one vertex should
unblock other vertices.
"""
function unblock!(v::T, blocked::BitArray, B::Vector{Set{T}}) where T <: Integer
    blocked[v] = false
    for w in B[v]
        delete!(B[v], w)
        if blocked[w]
            unblock!(w, blocked, B)
        end
    end
end

"""
    circuit{T<:Integer}(v::T, dg::::IsDirected, vis::JohnsonVisitor{T},
    allcycles::Vector{Vector{T}}, vmap::Vector{T}, startnode::T = v)

Return one step of the recursive version of simple cycle detection,
using a DFS algorithm.

* `v`: the vertex considered in this iteration of the DFS
* `dg`: the digraph from which cycles are computed
* `visitor`: Informations needed for the cycle computation, contains:
    * `stack`: the stack of parent vertices
    * `blocked`: tells whether a vertex has already been explored or not
    * `blockedmap`: mapping of the blocking / unblocking consequences
* `allcycles`: output containing the cycles already detected
* `vmap`: vector map containing the link from the old to the new nodes of the directed graph
* `startnode = v`: optional argument giving the starting node. In the first iteration,
 the same as v, otherwise it should be passed.

### Implementation Notes
Implements Johnson's CIRCUIT function. This is a recursive version.
Modifies the vector of cycles, when needed.

### References
- [Johnson](http://epubs.siam.org/doi/abs/10.1137/0204007)
"""
function circuit end
@traitfn function circuit(v::T, dg::::IsDirected, vis::JohnsonVisitor{T},
allcycles::Vector{Vector{T}}, vmap::Vector{T}, startnode::T=v) where T <: Integer
    done = false
    push!(vis.stack, v)
    vis.blocked[v] = true
    for w in outneighbors(dg, v)
        if w == startnode
            push!(allcycles, vmap[vis.stack])
            done = true
        elseif !vis.blocked[w]
            circuit(w, dg, vis, allcycles, vmap, startnode) && (done = true) #This is different from done = circuit(...). It keeps the previous value of done in the for loop
        end
    end
    if done
        unblock!(v, vis.blocked, vis.blockedmap)
    else
        for w in outneighbors(dg, v)
            if !in(vis.blockedmap[w], v)
                push!(vis.blockedmap[w], v)
            end
        end
    end
    pop!(vis.stack)
    return done
end

"""
    simplecycles(dg::::IsDirected)

Compute and return all cycles of the given directed graph using Johnson's algorithm.

### Performance
The number of cycles grows more than exponentially with the number of vertices,
you might want to use the algorithm with a ceiling -- `simplecycles_iter` -- on large directed graphs
(slightly slower). If you want to have an idea of the possible number of cycles,
look at function `maxsimplecycles(dg::DiGraph, byscc::Bool = true)`. If you only need
short cycles of a limited length, `simplecycles_limited_length` can be more efficient.

### References
- [Johnson](http://epubs.siam.org/doi/abs/10.1137/0204007)
"""
function simplecycles end
@traitfn function simplecycles(dg::::IsDirected)
    sccs = strongly_connected_components(dg)
    cycles = Vector{Vector{Int}}() # Pas très cohérent : devrait être du type de dg.
    for scc in sccs
        for i in 1:length(scc)
            wdg, vmap = induced_subgraph(dg, scc[i:end])
            visitor = JohnsonVisitor(wdg)
            circuit(1, wdg, visitor, cycles, vmap) # 1 is the startnode.
        end
    end
    return cycles
end

##########################################################
#### Iterative version, using Tasks, of the previous algorithms.
"""
    circuit_iter{T<:Integer}(v::T, dg::::IsDirected, vis::JohnsonVisitor{T},
    vmap::Vector{T}, cycle::Channel, startnode::T = v)

Execute one step of the recursive version of simple cycle detection, using a DFS algorithm.
Return `true` if a circuit has been found in the current exploration.

# Arguments
* v: the vertex considered in this iteration of the DFS
* dg: the digraph from which cycles are computed
* visitor: Informations needed for the cycle computation, contains:
    * stack: the stack of parent vertices
    * blocked: tells whether a vertex has already been explored or not
    * blockedmap: mapping of the blocking / unblocking consequences
* `vmap`: vector map containing the link from the old to the new nodes of the directed graph
* `cycle`: storage of the channel
* startnode = v: optional argument giving the starting node. In the first iteration,
the same as v, otherwise it should be passed.

### Implementation Notes
The CIRCUIT function from Johnson's algorithm, recursive and iterative version.
Produce a cycle when needed, can be used only inside a Channel.

### References
- [Johnson](http://epubs.siam.org/doi/abs/10.1137/0204007)
"""
function circuit_iter end
@traitfn function circuit_iter(v::T, dg::::IsDirected, vis::JohnsonVisitor{T},
vmap::Vector{T}, cycle::Channel, startnode::T=v) where T <: Integer
    done = false
    push!(vis.stack, v)
    vis.blocked[v] = true
    for w in outneighbors(dg, v)
        if w == startnode
            put!(cycle, vmap[vis.stack])
            done = true
        elseif !vis.blocked[w]
            circuit_iter(w, dg, vis, vmap, cycle, startnode) && (done = true) #This is different from done = circuit(...). It keeps the previous value of done in the for loop
        end
    end
    if done
        unblock!(v, vis.blocked, vis.blockedmap)
    else
        for w in outneighbors(dg, v)
            if !in(vis.blockedmap[w], v)
                push!(vis.blockedmap[w], v)
            end
        end
    end
    pop!(vis.stack)
    return done
end

"""
    itercycles(dg::::IsDirected, cycle::Channel)

Compute all cycles of the given directed graph, using Johnson's algorithm.

### Implementation Notes
Iterative version of the algorithm, using Channels to stop the exploration
after a given number of cycles.

### References
- [Johnson](http://epubs.siam.org/doi/abs/10.1137/0204007)
"""
function itercycles end
@traitfn function itercycles(dg::::IsDirected, cycle::Channel)
    sccs = strongly_connected_components(dg)
    for scc in sccs
        while length(scc) >= 1
            wdg, vmap = induced_subgraph(dg, scc)
            popfirst!(scc)
            visitor = JohnsonVisitor(wdg)
            circuit_iter(1, wdg, visitor, vmap, cycle)
        end
    end
end

"""
    simplecyclescount(dg::DiGraph, ceiling = 10^6)

Count the number of cycles in a directed graph, using Johnson's algorithm.
Return the minimum of the ceiling and the number of cycles.

### Implementation Notes
The `ceiling` is here to avoid memory overload if there are a lot of cycles in the graph.
Default value is 10^6, but it can be higher or lower. You can use the function
`maxsimplecycles(dg::DiGraph, byscc::Bool = true)` to get an idea of the
theoretical maximum number or cycles.

### References
- [Johnson](http://epubs.siam.org/doi/abs/10.1137/0204007)
"""
function simplecyclescount end
@traitfn function simplecyclescount(dg::::IsDirected, ceiling=10^6)
    len = 0
    for cycle in Iterators.take(Channel(c -> itercycles(dg, c)), ceiling)
        len += 1
    end
    return len
end

"""
    simplecycles_iter(dg::DiGraph, ceiling = 10^6)

Search all cycles of the given directed graph, using Johnson's algorithm,
up to the ceiling (to avoid memory overload).

### Implementation Notes
If the graph is small, the ceiling will not be reached and
`simplecycles(dg::DiGraph)` is more efficient. It avoids the overhead of the
counting and testing if the ceiling is reached. It returns all the cycles of the
directed graph if the `ceiling` is not reached, a subset of them otherwise.

To get an idea of the possible number of cycles, use function
`maxsimplecycles(dg::DiGraph, byscc::Bool = true) on the directed graph.

### References
- [Johnson](http://epubs.siam.org/doi/abs/10.1137/0204007)
"""
function simplecycles_iter end
@traitfn simplecycles_iter(dg::::IsDirected, ceiling=10^6) =
    collect(Iterators.take(Channel(c -> itercycles(dg, c)), ceiling))

"""
    simplecycleslength(dg::DiGraph, ceiling = 10^6)

Search all cycles of the given directed graph, using Johnson's algorithm,
and return a tuple representing the cycle length and the number of cycles.

### Implementation Notes
To get an idea of the possible number of cycles, using function
`maxsimplecycles(dg::DiGraph, byscc::Bool = true)` on the directed graph.

If the `ceiling` is reached (`ncycles = ceiling`), the output is only
a subset of the cycles lengths.

### References
- [Johnson](http://epubs.siam.org/doi/abs/10.1137/0204007)
"""
function simplecycleslength end
@traitfn function simplecycleslength(dg::::IsDirected, ceiling=10^6)
    ncycles = 0
    cyclelength = zeros(Int, nv(dg))
    for cycle in Iterators.take(Channel(c -> itercycles(dg, c)), ceiling)
        cyclelength[length(cycle)] += 1
        ncycles += 1
    end
    return cyclelength, ncycles
end


text:
#=
Tools for working with Markov Chains

@author : Spencer Lyon, Zac Cranko

@date: 07/10/2014

References
----------

http://quant-econ.net/jl/finite_markov.html
=#
import LightGraphs: DiGraph, period, attracting_components,
                    strongly_connected_components, is_strongly_connected

____solver_docs = """
solve x(P-I)=0 using either an eigendecomposition, lu factorization, or an
algorithm presented by Grassmann-Taksar-Heyman (GTH)

##### Arguments

- `p::Matrix` : valid stochastic matrix

##### Returns

- `x::Matrix`: A matrix whose columns contain stationary vectors of `p`

##### References

The following references were consulted for the GTH algorithm

- W. K. Grassmann, M. I. Taksar and D. P. Heyman, "Regenerative Analysis and
Steady State Distributions for Markov Chains, " Operations Research (1985),
1107-1116.
- W. J. Stewart, Probability, Markov Chains, Queues, and Simulation, Princeton
University Press, 2009.

"""

"""
Finite-state discrete-time Markov chain.

It stores useful information such as the stationary distributions, and
communication, recurrent, and cyclic classes, and allows simulation of state
transitions.

##### Fields

- `p::Matrix` The transition matrix. Must be square, all elements must be
positive, and all rows must sum to unity
"""
type MarkovChain{T<:Real}
    p::Matrix{T} # valid stochastic matrix

    function MarkovChain(p::Matrix)
        n, m = size(p)

        if n != m
            throw(ArgumentError("stochastic matrix must be square"))
        elseif any(p.<0)
            throw(ArgumentError("stochastic matrix must have nonnegative elements"))
        elseif any(x->!isapprox(sum(x), 1), [p[i, :] for i = 1:n])
            throw(ArgumentError("stochastic matrix rows must sum to 1"))
        end

        return new(p)
    end
end

# Provide constructor that infers T from eltype of matrix
MarkovChain{T<:Real}(p::Matrix{T}) = MarkovChain{T}(p)

"Number of states in the markov chain `mc`"
n_states(mc::MarkovChain) = size(mc.p, 1)

function Base.show(io::IO, mc::MarkovChain)
    println(io, "Discrete Markov Chain")
    println(io, "stochastic matrix:")
    println(io, mc.p)
end

# solve x(P-I)=0 by eigen decomposition
"$(____solver_docs)"
function eigen_solve{T}(p::Matrix{T})
    ef = eigfact(p')
    isunit = map(x->isapprox(x, 1), ef.values)
    x = real(ef.vectors[:, isunit])
    x ./= sum(x, 1) # normalisation
    for i = 1:length(x)
        x[i] = isapprox(x[i], zero(T)) ? zero(T) : x[i]
    end
    any(x .< 0) && warn("something has gone wrong with the eigen solve")
    x
end

# solve x(P-I)=0 by lu decomposition
"$(____solver_docs)"
function lu_solve{T}(p::Matrix{T})
    n, m = size(p)
    x   = vcat(Array(T, n-1), one(T))
    u   = lufact(p' - one(p))[:U]
    for i = n-1:-1:1 # backsubstitution
        x[i] = -sum([x[j]*u[i, j] for j=i:n])/u[i, i]
    end
    x ./= norm(x, 1) # normalisation
    for i = 1:length(x)
        x[i] = isapprox(x[i], zero(T)) ? zero(T) : x[i]
    end
    any(x .< 0) && warn("something has gone wrong with the lu solve")
    x
end

"$(____solver_docs)"
gth_solve{T<:Integer}(a::Matrix{T}) = gth_solve(convert(Array{Rational, 2}, a))

# solve x(P-I)=0 by the GTH method
function gth_solve{T<:Real}(original::Matrix{T})
    a = copy(original)
    n = size(a, 1)
    x = zeros(T, n)

    @inbounds for k in 1:n-1
        scale = sum(a[k, k+1:n])
        if scale <= zero(T)
            # There is one (and only one) recurrent class contained in
            # {1, ..., k};
            # compute the solution associated with that recurrent class.
            n = k
            break
        end
        a[k+1:n, k] /= scale

        for j in k+1:n, i in k+1:n
            a[i, j] += a[i, k] * a[k, j]
        end
    end

    # backsubstitution
    x[n] = 1
    @inbounds for k in n-1:-1:1, i in k+1:n
        x[k] += x[i] * a[i, k]
    end

    # normalisation
    x / sum(x)
end

"""
Find the recurrent classes of the `MarkovChain`

##### Arguments

- `mc::MarkovChain` : MarkovChain instance containing a valid stochastic matrix

##### Returns

- `x::Vector{Vector}`: A `Vector` containing `Vector{Int}`s that describe the
recurrent classes of the transition matrix for p

"""
recurrent_classes(mc::MarkovChain) = attracting_components(DiGraph(mc.p))

communication_classes(mc::MarkovChain) = strongly_connected_components(DiGraph(mc.p))

is_irreducible(mc::MarkovChain) =  is_strongly_connected(DiGraph(mc.p))

is_aperiodic(mc::MarkovChain) = period(mc) == 1

function period(mc::MarkovChain)
    g = DiGraph(mc.p)
    recurrent = attracting_components(g)
    periods   = Int[]

    for r in recurrent
        push!(periods, period(g[r]))
    end
    return lcm(periods)
end
"""
calculate the stationary distributions associated with a N-state markov chain

##### Arguments

- `mc::MarkovChain` : MarkovChain instance containing a valid stochastic matrix
- `;method::Symbol(:gth)`: One of `gth`, `lu`, and `eigen`; specifying which
of the three `_solve` methods to use.

##### Returns

- `dists::Matrix{Float64}`: N x M matrix where each column is a stationary
distribution of `mc.p`

"""
function mc_compute_stationary{T}(mc::MarkovChain{T}; method::Symbol=:gth)
    solvers = Dict(:gth => gth_solve, :lu => lu_solve, :eigen => eigen_solve)
    solve   = solvers[method]

    recurrent = recurrent_classes(mc)

    # unique recurrent class
    length(recurrent) == 1 && return solve(mc.p)

    # more than one recurrent classes
    stationary_dists = Array(T, n_states(mc), length(recurrent))
    for (i, class) in enumerate(recurrent)
        dist        = zeros(T, n_states(mc))
        dist[class] = solve(mc.p[class, class]) # recast to correct dimension
        stationary_dists[:, i] = dist
    end
    return stationary_dists
end
"""
Simulate a Markov chain starting from an initial state

##### Arguments

- `mc::MarkovChain` : MarkovChain instance containing a valid stochastic matrix
- `init::Int(rand(1:n_states(mc)))` : The index of the initial state. This should
be an integer between 1 and `n_states(mc)`
- `sample_size::Int(1000)`: The number of samples to collect
- `;burn::Int(0)`: The burn in length. Routine drops first `burn` of the
`sample_size` total samples collected

##### Returns

- `samples::Vector{Int}`: Vector of simulated states

"""
function mc_sample_path(mc::MarkovChain,
                        init::Int=rand(1:n_states(mc)),
                        sample_size::Int=1000;
                        burn::Int=0)
    samples = Array(Int, sample_size)
    samples[1] = init
    mc_sample_path!(mc, samples)
    samples[burn+1:end]
end

"""
Simulate a Markov chain starting from an initial distribution

##### Arguments

- `mc::MarkovChain` : MarkovChain instance containing a valid stochastic matrix
- `init::Vector` : A vector of length `n_state(mc)` specifying the number
probability of being in seach state in the initial period
- `sample_size::Int(1000)`: The number of samples to collect
- `;burn::Int(0)`: The burn in length. Routine drops first `burn` of the
`sample_size` total samples collected

##### Returns

- `samples::Vector{Int}`: Vector of simulated states

"""
function mc_sample_path(mc::MarkovChain,
                        init::Vector = vcat(1, zeros(n_states(mc)-1)),
                        sample_size::Int=1000; burn::Int=0)
    # ensure floating point input for Categorical()
    init = convert(Vector{Float64}, init)
    return mc_sample_path(mc, rand(Categorical(init)), sample_size, burn=burn)
end

# duplicate python functionality
function mc_sample_path(mc::MarkovChain,
                        state::Real,
                        sample_size::Int=1000; burn::Int=0)
    init = zeros(n_states(mc))
    init[state] = 1
    return mc_sample_path(mc, init, sample_size, burn=burn)
end

"""
Fill `samples` with samples from the Markov chain `mc`

##### Arguments

- `mc::MarkovChain` : MarkovChain instance containing a valid stochastic matrix
- `samples::Array{Int}` : Pre-allocated vector of integers to be filled with
samples from the markov chain `mc`. The first element will be used as the
initial state and all other elements will be over-written.

##### Returns

None modifies `samples` in place
"""
function mc_sample_path!(mc::MarkovChain, samples::Array)
    # ensure floating point input for Categorical()
    p       = convert(Matrix{Float64}, mc.p)
    dist    = [Categorical(vec(p[i, :])) for i=1:n_states(mc)]
    for t=2:length(samples)
        samples[t] = rand(dist[samples[t-1]])
    end
    Void
end

function simulate(mc::MarkovChain,
                  init::Vector = vcat(1, zeros(n_states(mc)-1)),
                  sample_size=1000)
    return mc_sample_path(mc, init, sample_size)
end


text:
import Base.cross
import WiltonInts84

export MWSingleLayer3D, MWHyperSingular, MWWeaklySingular
export MWDoubleLayer3D
export PlaneWaveMW
export TangTraceMW, CrossTraceMW

export curl

abstract type MaxwellOperator3D <: IntegralOperator end
abstract type MaxwellOperator3DReg <: MaxwellOperator3D end

immutable KernelValsMaxwell3D{T,U,P,Q}
    "gamma = im * wavenumber"
    gamma::U
    vect::P
    dist::T
    green::U
    gradgreen::Q
end

function kernelvals(biop::MaxwellOperator3D, p, q)

    γ = biop.gamma
    r = p.cart - q.cart
    R = norm(r)
    γR = γ*R

    expn = exp(-γR)
    green = expn / (4pi*R)
    gradgreen = -(γ + 1/R) * green / R * r

    KernelValsMaxwell3D(γ, r, R, green, gradgreen)
end

function kernelvals(kernel::MaxwellOperator3DReg, p, q)

    γ = kernel.gamma
    r = p.cart - q.cart
    R = norm(r)
    γR = γ*R

    Exp = exp(-γ*R)
    green = (Exp - 1 + γR - 0.5*γR^2) / (4pi*R)
    gradgreen = ( - (γR + 1)*Exp + (1 - 0.5*γR^2) ) * (r/R^3) / (4π)

    KernelValsMaxwell3D(γ, r, R, green, gradgreen)
end

immutable MWSingleLayer3D{T,U} <: MaxwellOperator3D
  gamma::T
  α::U
  β::U
end

MWSingleLayer3D(gamma)  = MWSingleLayer3D(gamma, -gamma, -1/(gamma))
MWWeaklySingular(gamma) = MWSingleLayer3D(gamma, 1, 0)
MWHyperSingular(gamma)  = MWSingleLayer3D(gamma, 0, 1)

module Maxwell3D

    using ..BEAST

    """
        singlelayer(;gamma, alpha, beta)

    Bilinear form given by:

    ```math
        α ∬_{Γ×Γ} j(x)⋅k(y) G_{γ}(x,y) + β ∬_{Γ×Γ} div j(x) div k(y) G_{γ}(x,y)
    ```

    with ``G_{γ} = e^{-γ|x-y|} / 4π|x-y|``.
    """
    singlelayer(;
        gamma=error("propagation constant required"),
        alpha=-gamma,
        beta=-1/gamma) = MWSingleLayer3D(gamma, alpha, beta)

    planewave(;
            direction    = error("missing arguement `direction`"),
            polarization = error("missing arguement `polarization`"),
            wavenumber   = error("missing arguement `wavenumber`"),
            amplitude    = one(real(typeof(wavenumber)))) =
        BEAST.PlaneWaveMW(direction, polarization, wavenumber, amplitude)
end

export Maxwell3D

struct MWSingleLayer3DReg{T,U} <: MaxwellOperator3DReg
    gamma::T
    α::U
    β::U
end

struct MWSingleLayer3DSng{T,U} <: MaxwellOperator3D
    gamma::T
    α::U
    β::U
end

scalartype(op::MaxwellOperator3D) = typeof(op.gamma)

regularpart(op::MWSingleLayer3D) = MWSingleLayer3DReg(op.gamma, op.α, op.β)
singularpart(op::MWSingleLayer3D) = MWSingleLayer3DSng(op.gamma, op.α, op.β)

function quaddata(op::MaxwellOperator3D, g::RefSpace, f::RefSpace, tels, bels)

    tqd = quadpoints(g, tels, (2,6))
    bqd = quadpoints(f, bels, (3,7))

    return QuadData(tqd, bqd)
end

# use Union type so this code can be shared between the operator
# and its regular part.
MWSL3DGen = Union{MWSingleLayer3D,MWSingleLayer3DReg}
function integrand(biop::MWSL3DGen, kerneldata, tvals, tgeo, bvals, bgeo)

  gx = tvals[1]
  fy = bvals[1]

  dgx = tvals[2]
  dfy = bvals[2]

  G = kerneldata.green
  γ = kerneldata.gamma

  α = biop.α
  β = biop.β

  t = (α * dot(gx, fy) + β * (dgx*dfy)) * G
end

immutable MWDoubleLayer3D{T} <: MaxwellOperator3D
  gamma::T
end

immutable MWDoubleLayer3DSng{T} <: MaxwellOperator3D
    gamma::T
end

immutable MWDoubleLayer3DReg{T} <: MaxwellOperator3DReg
    gamma::T
end

regularpart(op::MWDoubleLayer3D) = MWDoubleLayer3DReg(op.gamma)
singularpart(op::MWDoubleLayer3D) = MWDoubleLayer3DSng(op.gamma)

const MWDL3DGen = Union{MWDoubleLayer3D,MWDoubleLayer3DReg}
function integrand(biop::MWDL3DGen, kerneldata, tvals, tgeo, bvals, bgeo)
    g = tvals[1]
    f = bvals[1]
    ∇G = kerneldata.gradgreen
    g ⋅ (∇G × f)
end

function innerintegrals!(op::MWSingleLayer3DSng, p, g, f, t, s, z, strat::WiltonSEStrategy, dx)

    γ = op.gamma
    x = cartesian(p)
    n = cross(s[1]-s[3],s[2]-s[3])
    n /= norm(n)
    ρ = x - ((x-s[1]) ⋅ n) * n

    scal, vec = WiltonInts84.wiltonints(s[1], s[2], s[3], x, Val{1})

    # \int \frac{1}{4 \pi R}
    ∫G = (scal[2] - γ*scal[3] + 0.5*γ^2*scal[4]) / (4π)
    # \int \frac{y}{4 \pi R}
    ∫Gy = SVector((
        (vec[2][1] + scal[2]*ρ[1] - γ*(vec[3][1]+scal[3]*ρ[1]) + 0.5*γ^2*(vec[4][1]+scal[4]*ρ[1]))/(4π),
        (vec[2][2] + scal[2]*ρ[2] - γ*(vec[3][2]+scal[3]*ρ[2]) + 0.5*γ^2*(vec[4][2]+scal[4]*ρ[2]))/(4π),
        (vec[2][3] + scal[2]*ρ[3] - γ*(vec[3][3]+scal[3]*ρ[3]) + 0.5*γ^2*(vec[4][3]+scal[4]*ρ[3]))/(4π),
    ))

    c₁ = op.α
    c₂ = op.β

    α = 1 / volume(t) / volume(s) / 4
    for i in 1 : numfunctions(g)
        a = t[i]
        g = x - a
        dg = 2

        for j in 1 : numfunctions(f)
            b = s[j]

            ∫Gf = SVector(∫Gy[1]-∫G*b[1], ∫Gy[2]-∫G*b[2], ∫Gy[3]-∫G*b[3])
            ∫Gdf = 2 * ∫G
            dg∫Gf = g[1]*∫Gf[1] + g[2]*∫Gf[2] + g[3]*∫Gf[3]
            z[i,j] += ( α*c₁*dg∫Gf + α*c₂*dg*∫Gdf ) * dx

        end # next j
    end #

end

function innerintegrals!(op::MWDoubleLayer3DSng, p, g, f, t, s, z, strat::WiltonSEStrategy, dx)

    γ = op.gamma
    x = cartesian(p)
    n = cross(s[1]-s[3],s[2]-s[3])
    n /= norm(n)
    ρ = x - ((x-s[1]) ⋅ n) * n

    scal, vec, grad = WiltonInts84.wiltonints(s[1], s[2], s[3], x, Val{1})

    # \int \nabla G_s with G_s = \nabla (1/R + 0.5*γ^2*R) / (4\pi)
    ∫∇G = (-grad[1] - 0.5*γ^2*grad[3]) / (4π)

    α = 1 / volume(t) / volume(s) / 4
    for i in 1 : numfunctions(g)
        a = t[i]
        g = (x - a)

        for j in 1 : numfunctions(f)
            b = s[j]

            z[i,j] += ( α * ( (x-b) × g ) ⋅ ∫∇G ) * dx

        end # next j
    end #

end

function select_quadrule()
    try
        Pkg.installed("BogaertInts10")
        info("`BogaertInts10` detected: enhanced quadrature enabled.")
        @eval using BogaertInts10
        @eval include("bogaertints.jl")
        @eval quadrule(op::MaxwellOperator3D, g::RTRefSpace, f::RTRefSpace, i, τ, j, σ, qd) = qrib(op, g, f, i, τ, j, σ, qd)
    catch
        info("Cannot find package `BogaertInts10`. Default quadrature strategy used.")
        @eval quadrule(op::MaxwellOperator3D, g::RTRefSpace, f::RTRefSpace, i, τ, j, σ, qd) = qrdf(op, g, f, i, τ, j, σ, qd)
    end
end

select_quadrule()

function qrib(op::MaxwellOperator3D, g::RTRefSpace, f::RTRefSpace, i, τ, j, σ, qd)
  # defines coincidence of points
  dtol = 1.0e3 * eps(eltype(eltype(τ.vertices)))

  # decides on whether to use singularity extraction
  xtol = 0.2

  k = norm(op.gamma)

  hits = 0
  xmin = xtol
  for t in τ.vertices
    for s in σ.vertices
      d = norm(t-s)
      xmin = min(xmin, k*d)
      if d < dtol
        hits +=1
        break
      end
    end
  end

  hits == 3   && return BogaertSelfPatchStrategy(5)
  hits == 2   && return BogaertEdgePatchStrategy(8, 4)
  hits == 1   && return BogaertPointPatchStrategy(2, 3)
  xmin < xtol && return WiltonSEStrategy(
    qd.tpoints[1,i],
    DoubleQuadStrategy(
      qd.tpoints[2,i],
      qd.bpoints[2,j],
    ),
  )
  return DoubleQuadStrategy(
    qd.tpoints[1,i],
    qd.bpoints[1,j],
  )

end

function qrdf(op::MaxwellOperator3D, g::RTRefSpace, f::RTRefSpace, i, τ, j, σ, qd)
  # defines coincidence of points
  dtol = 1.0e3 * eps(eltype(eltype(τ.vertices)))

  # decides on whether to use singularity extraction
  xtol = 0.2

  k = norm(op.gamma)

  hits = 0
  xmin = xtol
  for t in τ.vertices
    for s in σ.vertices
      d = norm(t-s)
      xmin = min(xmin, k*d)
      if d < dtol
        hits +=1
        break
      end
    end
  end

  xmin < xtol && return WiltonSEStrategy(
    qd.tpoints[1,i],
    DoubleQuadStrategy(
      qd.tpoints[2,i],
      qd.bpoints[2,j],
    ),
  )
  return DoubleQuadStrategy(
    qd.tpoints[1,i],
    qd.bpoints[1,j],
  )

end


text:
module PlutoPlotly

using Reexport
@reexport using PlotlyBase

export PlutoPlot, get_plotly_version, change_plotly_version, check_plotly_version, force_pluto_mathjax_local, htl_js, add_plotly_listener!, add_class!, remove_class!, add_js_listener!
export plot, push_script!
include("../notebooks/wrapper.jl")

function __init__()
	if !is_inside_pluto()
		@warn "You loaded this package outside of Pluto, this is not the intended behavior and you should use either PlotlyBase or PlotlyJS directly"
	end
end

end