text:
const Maybe{T} = Union{Nothing, T}

"""
    NoDefault

Type describes a field should have no default value.
"""
struct NoDefault end

"""
    const no_default = NoDefault()

Constant instance for [`NoDefault`](@ref) that
describes a field should have no default value.
"""
const no_default = NoDefault()

abstract type JLExpr end

"""
    JLFunction <: JLExpr

Type describes a Julia function declaration expression.

# Example

Construct a function expression

```julia
julia> JLFunction(;name=:foo, args=[:(x::T)], body= quote 1+1 end, head=:function, whereparams=[:T])
function foo(x::T) where {T}
    #= REPL[25]:1 =#
    1 + 1
end
```

Decompose a function expression

```julia
julia> ex = :(function foo(x::T) where {T}
           #= REPL[25]:1 =#
           1 + 1
       end)
:(function foo(x::T) where T
      #= REPL[26]:1 =#
      #= REPL[26]:3 =#
      1 + 1
  end)

julia> jl = JLFunction(ex)
function foo(x::T) where {T}
    #= REPL[26]:1 =#
    #= REPL[26]:3 =#
    1 + 1
end
```

Generate `Expr` from `JLFunction`

```julia
julia> codegen_ast(jl)
:(function foo(x::T) where T
      #= REPL[26]:1 =#
      #= REPL[26]:3 =#
      1 + 1
  end)
```
"""
mutable struct JLFunction <: JLExpr
    head::Symbol  # function def must have a head
    name::Any  # name can be nothing, Symbol, Expr
    args::Vector{Any}
    kwargs::Maybe{Vector{Any}}
    whereparams::Maybe{Vector{Any}}
    body::Any
    line::Maybe{LineNumberNode}
    doc::Maybe{String}
end

function JLFunction(;
        head=:function, name=nothing,
        args=[], kwargs=nothing,
        whereparams=nothing, body=Expr(:block),
        line=nothing, doc=nothing
    )
    JLFunction(head, name, args, kwargs, whereparams, body, line, doc)
end

"""
    JLField <: JLExpr
    JLField(name, type, line)

Type describes a Julia field in a Julia struct.
"""
mutable struct JLField <: JLExpr
    name::Symbol
    type::Any
    doc::Maybe{String}
    line::Maybe{LineNumberNode}
end

function JLField(;name, type=Any, doc=nothing, line=nothing)
    JLField(name, type, doc, line)
end

"""
    JLKwField <: JLExpr
    JLKwField(name, type, line, default=no_default)

Type describes a Julia field that can have a default value in a Julia struct.
"""
mutable struct JLKwField <: JLExpr
    name::Symbol
    type::Any
    doc::Maybe{String}
    line::Maybe{LineNumberNode}
    default::Any
end

function JLKwField(;name, type=Any, doc=nothing, line=nothing, default=no_default)
    JLKwField(name, type, doc, line, default)
end

"""
    JLStruct <: JLExpr

Type describes a Julia struct.

# Example

Construct a Julia struct.

```julia
julia> JLStruct(;name=:Foo, typevars=[:T], fields=[JLField(;name=:x, type=Int)])
struct Foo{T}
    x::Int64
end
```

Decompose a Julia struct expression

```julia
julia> ex = :(struct Foo{T}
           x::Int64
       end)
:(struct Foo{T}
      #= REPL[31]:2 =#
      x::Int64
  end)

julia> jl = JLStruct(ex)
struct Foo{T}
    #= REPL[31]:2 =#
    x::Int64
end
```

Generate a Julia struct expression

```julia
julia> codegen_ast(jl)
:(struct Foo{T}
      #= REPL[31]:2 =#
      x::Int64
  end)
```
"""
mutable struct JLStruct <: JLExpr
    name::Symbol
    ismutable::Bool
    typevars::Vector{Any}
    supertype::Any
    fields::Vector{JLField}
    constructors::Vector{JLFunction}
    line::Maybe{LineNumberNode}
    doc::Maybe{String}
    misc::Any
end

function JLStruct(;
    name, ismutable=false,
    typevars=[], supertype=nothing,
    fields=JLField[], constructors=JLFunction[],
    line=nothing, doc=nothing, misc=nothing)
    JLStruct(name, ismutable, typevars, supertype, fields, constructors, line, doc, misc)
end

"""
    JLKwStruct <: JLExpr

Type describes a Julia struct that allows keyword definition of defaults.
"""
mutable struct JLKwStruct <: JLExpr
    name::Symbol
    typealias::Maybe{String}
    ismutable::Bool
    typevars::Vector{Any}
    supertype::Any
    fields::Vector{JLKwField}
    constructors::Vector{JLFunction}
    line::Maybe{LineNumberNode}
    doc::Maybe{String}
    misc::Any
end

function JLKwStruct(;name, typealias=nothing,
    ismutable=false, typevars=[], supertype=nothing,
    fields=JLField[], constructors=JLFunction[],
    line=nothing, doc=nothing, misc=nothing)
    JLKwStruct(name, typealias, ismutable, typevars, supertype, fields, constructors, line, doc, misc)
end

"""
    JLIfElse <: JLExpr

`JLIfElse` describes a Julia `if ... elseif ... else ... end` expression. It allows one to easily construct
such expression by inserting condition and code block via a map.

# Example

### Construct JLIfElse object

One can construct an `ifelse` as following

```julia
julia> jl = JLIfElse()
nothing

julia> jl.map[:(foo(x))] = :(x = 1 + 1)
:(x = 1 + 1)

julia> jl.map[:(goo(x))] = :(y = 1 + 2)
:(y = 1 + 2)

julia> jl.otherwise = :(error("abc"))
:(error("abc"))

julia> jl
if foo(x)
    x = 1 + 1
elseif goo(x)
    y = 1 + 2
else
    error("abc")
end
```

### Generate the Julia `Expr` object

to generate the corresponding `Expr` object, one can call [`codegen_ast`](@ref).

```julia
julia> codegen_ast(jl)
:(if foo(x)
      x = 1 + 1
  elseif goo(x)
      y = 1 + 2
  else
      error("abc")
  end)
```
"""
mutable struct JLIfElse <: JLExpr
    map::OrderedDict{Any, Any}
    otherwise::Any
end

JLIfElse() = JLIfElse(OrderedDict(), nothing)

"""
    JLMatch <: JLExpr

`JLMatch` describes a Julia pattern match expression defined by
[`MLStyle`](https://github.com/thautwarm/MLStyle.jl). It allows
one to construct such expression by simply assign each code block
to the corresponding pattern expression.

# Example

One can construct a `MLStyle` pattern matching expression
easily by assigning the corresponding pattern and its result
to the `map` field.

```julia
julia> jl = JLMatch(:x)
#= line 0 =#
nothing

julia> jl = JLMatch(:x)
#= line 0 =#
nothing

julia> jl.map[1] = true
true

julia> jl.map[2] = :(sin(x))
:(sin(x))

julia> jl
#= line 0 =#
@match x begin
    1 => true
    2 => sin(x)
    _ =>     nothing
end
```

to generate the corresponding Julia `Expr` object, one can call [`codegen_ast`](@ref).

```julia
julia> codegen_ast(jl)
:(let
      true
      var"##return#263" = nothing
      var"##265" = x
      if var"##265" isa Int64
          #= line 0 =#
          if var"##265" === 1
              var"##return#263" = let
                      true
                  end
              #= unused:1 =# @goto var"####final#264#266"
          end
          #= line 0 =#
          if var"##265" === 2
              var"##return#263" = let
                      sin(x)
                  end
              #= unused:1 =# @goto var"####final#264#266"
          end
      end
      #= line 0 =#
      begin
          var"##return#263" = let
                  nothing
              end
          #= unused:1 =# @goto var"####final#264#266"
      end
      (error)("matching non-exhaustive, at #= line 0 =#")
      #= unused:1 =# @label var"####final#264#266"
      var"##return#263"
  end)
```
"""
struct JLMatch <: JLExpr
    item::Any
    map::OrderedDict{Any, Any}
    fallthrough::Any
    mod::Module
    line::LineNumberNode
end

"""
    JLMatch(item)

Generate an empty `JLMatch` object with given item expression.
`item` can be a `Symbol` or an `Expr`.
"""
JLMatch(item) = JLMatch(item, OrderedDict(), nothing, Main, LineNumberNode(0))

"""
    JLMatch(;item, map=OrderedDict(), fallthrough=nothing, mod=Main, line=LineNumberNode(0))

Create a `JLMatch` object from keyword arguments.

# Kwargs

- `item`: item to match
- `map`: the pattern=>result map, should be an `OrderedDict`.
- `fallthrough`: the result of fallthrough pattern `_`.
- `mod`: module to evaluate the expression.
- `line`: line number `LineNumberNode`.
"""
JLMatch(;item, map=OrderedDict(), fallthrough=nothing, mod=Main, line=LineNumberNode(0)) =
    JLMatch(item, map, fallthrough, mod, line)

struct JLFor <: JLExpr
    vars::Vector{Any}
    iterators::Vector{Any}
    kernel::Any
end

JLFor(;vars=[], iterators=[], kernel=nothing) = JLFor(vars, iterators, kernel)


text:
using Documenter
using LearningHorse.NeuralNetwork
makedocs(;
    sitename="LearningHorse",
    pages = [
        "Home" => "index.md",
        "Manual" => [
            "Regression" => "Regression.md",
            "Classification" => "Classification.md",
            "Preprocessing" => "Preprocessing.md",
            "LossFunction" => "LossFunction.md",
            "NeuralNetwork" => "neuralnetwork.md",
        ],
        "日本語版マニュアル" => "Japanese.md"
    ]
)

ext: jl


text:
export PoissonSparseIterative,solve!
using IterativeSolvers

struct PoissonSparseIterative
    spLxy::SparseMatrixCSC{Float64,Int64}

    function PoissonSparseIterative(sp)
        new(sp)
    end
end

function optimal_w(x1D)
    h=1/sqrt(length(x1D))
    w=2/(1+sin(π*h))
end

solve_sor!(x,a,b,nstep) = IterativeSolvers.sor!(x,a,b,optimal_w(x),maxiter=nstep)
solve_ssor!(x,a,b,nstep) = IterativeSolvers.ssor!(x,a,b,optimal_w(x),maxiter=nstep)
solve_gauss_seidel!(x,a,b,nstep) = IterativeSolvers.gauss_seidel!(x,a,b,maxiter=nstep)
solve_jacobi!(x,a,b,nstep) = IterativeSolvers.jacobi!(x,a,b,maxiter=nstep)

function solve!(Xxy,Bxy,sp::PoissonSparseIterative,method_name::String="cg")
    b=view(Bxy,1:length(Bxy))
    x=view(Xxy,1:length(Xxy))
    A=sp.spLxy

    method_name=="cg" && return IterativeSolvers.cg!(x,A,b;log=true,tol = 1e-10,maxiter=1000)
    # method_name="jacobi" && solve_cg!(Xxy,Bxy,sp)
    # method_name="gauss_seidel" && solve_cg!(Xxy,Bxy,sp)

    # method_name=="gauss_seidel" && return IterativeSolvers.gauss_seidel!(x,A,b,maxiter=1000)
    method_name=="gauss_seidel" && return solve_stationary!(Xxy,Bxy,sp::PoissonSparseIterative,solve_gauss_seidel!)
    method_name=="jacobi" && return solve_stationary!(Xxy,Bxy,sp::PoissonSparseIterative,solve_jacobi!)
    method_name=="sor" && return solve_stationary!(Xxy,Bxy,sp::PoissonSparseIterative,solve_sor!)
    method_name=="ssor" && return solve_stationary!(Xxy,Bxy,sp::PoissonSparseIterative,solve_ssor!)

    # method_name="ssor" && solve_cg!(Xxy,Bxy,sp)

end

function solve_stationary!(Xxy,Bxy,sp::PoissonSparseIterative,iterative_algorithm)
    b=view(Bxy,1:length(Bxy))
    x=view(Xxy,1:length(Xxy))
    A=sp.spLxy

    #LP: stationary methods looks to have problem with views...
    x1D=zeros(length(Xxy))
    b1D=zeros(length(Xxy))
    copyto!(b1D,b)
    copyto!(x1D,x)

    residual=Vector{Float64}()
    nstep=10
    # h=1/size(Bxy,1)
    # w=2/(1+sin(π*h))
    # @show w
    res=typemax(Float64)
    iter=1
    while res>1.e-8 && iter<2000
        iterative_algorithm(x1D,A,b1D,nstep)
        res=norm(A*x1D-b1D)
        for i=1:nstep
            push!(residual,res)
        end
        iter+=nstep
    end
    @show iterative_algorithm,iter
    copyto!(x,x1D)
    x,residual
end


text:
using Proj4
using Test

@testset "Proj4" begin

println("""
C library version: $(Proj4.version)  [\"$(Proj4._get_release())\"]
geodesic support: $(Proj4.has_geodesic_support)
""")

# Some very basic sanity checking
wgs84 = Projection("+proj=longlat +datum=WGS84 +no_defs")
utm56 = Projection("+proj=utm +zone=56 +south +datum=WGS84 +units=m +no_defs")
nad83 = Projection(Proj4.epsg[4269])

# Reference data computed using GeographicLib's GeoConvert tool
@test transform(wgs84, utm56, [150.0, -27.0, 0.0]) ≈
    [202273.912995055, 7010024.033113679, 0] atol=1e-6

dup_wgs84 = Projection("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs")
@test transform(wgs84, dup_wgs84, [150 -27 0]) ≈
    [150 -27 0] atol=1e-6
@test transform(dup_wgs84, wgs84, [150 -27 0]) ≈
    [150 -27 0] atol=1e-6

# Taken from https://github.com/proj4js/proj4js/blob/master/test/test.js
sweref99tm = Projection("+proj=utm +zone=33 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs")
rt90 = Projection("+lon_0=15.808277777799999 +lat_0=0.0 +k=1.0 +x_0=1500000.0 +y_0=0.0 +proj=tmerc +ellps=bessel +units=m +towgs84=414.1,41.3,603.1,-0.855,2.141,-7.023,0 +no_defs")
@test transform(sweref99tm, rt90, [319180 6399862]) ≈
    [1271137.927154 6404230.291456] atol=1e-2
@test transform(sweref99tm, rt90, [319480 6397862;
                                   329200 6599800]) ≈
    [1271414.272854 6402225.564811;
     1283568.895883 6604160.216834] atol=1e-2
@test transform(sweref99tm, rt90, [319480 6397862;
                                   329200 6599800;
                                   319480 6397862;
                                   329200 6599800]) ≈
    [1271414.272854 6402225.564811;
     1283568.895883 6604160.216834;
     1271414.272854 6402225.564811;
     1283568.895883 6604160.216834] atol=1e-2
@test transform(rt90, sweref99tm, [1271414.272854 6402225.564811;
                                   1283568.895883 6604160.216834;
                                   1271414.272854 6402225.564811;
                                   1283568.895883 6604160.216834]) ≈
    [319480 6397862;
     329200 6599800;
     319480 6397862;
     329200 6599800] atol=1e-2

# Just to verify that we get different Projections from different (projection) definitions
wgs84 = Projection(Proj4.epsg[4326]) # World
svy21 = Projection(Proj4.epsg[3414]) # Singapore
proj = Projection(Proj4.epsg[3906]) # based on the Australian National Spheroid (bessel ellipse)

wgs84_a, wgs84_es = spheroid_params(wgs84)
svy21_a, svy21_es = spheroid_params(svy21)
proj_a, proj_es = spheroid_params(proj)
@test wgs84_es  ≈ svy21_es atol=1e-6 # same ellipsoid
@test svy21_es  ≈ 0.0066943799901413165 atol=1e-6
@test proj_es   ≈ 0.006674372231802145 atol=1e-6
@test abs(proj_es - wgs84_es) > 1e-6
@test abs(proj_a - wgs84_a) > 500.0

@test !compare_datums(svy21, wgs84)
@test compare_datums(svy21, svy21)
@test is_latlong(wgs84)
@test !is_latlong(utm56)
@test !is_latlong(sweref99tm)
@test !is_latlong(svy21)
@test is_geocent(Projection(Proj4.epsg[4328])) # WGS 84 (geocentric)
@test !is_geocent(wgs84) # WGS 84 (geocentric)

#=
# TODO(cjf): It's unclear exatly what the below was meant to achieve.
# Commented out for now since the tests are broken.
function test_point(point)
    rev = reverse(vec(point))
    rad = deg2rad(point)
    rev_rad = deg2rad(point)

    @test lonlat2xy!(rev, wgs84) ≈ lonlat2xy(reverse(vec(point)), wgs84) atol=1e-6 # inplace
    @test rev ≈ rev_rad atol=1e-6
    @test point ≈ rad2deg(rad) atol=1e-6
end

for point in ([-73.78, 40.64], # JFK
              [-73.78  40.64],
              [103.99,  1.36], # SIN
              [103.99   1.36])
    test_point(point)
end
=#

p1, p2 = [-73.78, 40.64],[103.99, 30.36] # p1, p2 in degrees
proj1, proj2 = wgs84, proj # chosen to have different ellipses
q1, q2 = lonlat2xy(p1, proj1), lonlat2xy(p2, proj2)

if Proj4.has_geodesic_support
    r2 = transform(proj2, proj1, q2)
    dist_q1r2, azi1, azi2 = geod_inverse(q1, r2, proj1)
    dest, azi = geod_direct(q1, azi1, dist_q1r2, proj1)
    @test dest ≈ r2 atol=1e-6
    @test azi ≈ azi2 atol=1e-6
    @test geod_destination(q1, azi1, dist_q1r2, proj1) ≈ r2 atol=1e-6

    # # It is not necessarily symmetric:
    # # You can go q1(azi1) -[dist1]-> r2,
    # # but not    r2(azi2) -[dist1]-> q1
    # # So the following statements are false in general:
    # dest, azi = geod_direct(r2, azi2, dist1, proj1)
    # @test dest ≈ q1 atol=1e-6
    # @test azi ≈ azi1 atol=1e-6
    # @test geod_destination(r2, azi2, dist1, proj1) ≈ q1 atol=1e-6

    # # To get the reverse azimuth to move from r2 back to q1,
    # # make another call to geod_inverse:
    dist_r2q1, azi1, azi2 = geod_inverse(r2, q1, proj1)
    dest, azi = geod_direct(r2, azi1, dist_r2q1, proj1)
    @test dest ≈ q1 atol=1e-6
    @test azi ≈ azi2 atol=1e-6
    @test geod_destination(r2, azi1, dist_r2q1, proj1) ≈ q1 atol=1e-6

    # The distances from both calls to geod_inverse should be the same, i.e.
    @test dist_r2q1 ≈ dist_q1r2 atol=1e-6

    # Doublecheck when we perform the operations in the other projection:
    r1 = transform(proj1, proj2, q1)
    dist_r1q2, azi1, azi2 = geod_inverse(r1, q2, proj2)
    dest, azi = geod_direct(r1, azi1, dist_r1q2, proj2)
    @test dest ≈ q2 atol=1e-6
    @test azi ≈ azi2 atol=1e-6
    @test geod_destination(r1, azi1, dist_r1q2, proj2) ≈ q2 atol=1e-6
    dist_q2r1, azi1, azi2 = geod_inverse(q2, r1, proj1)
    dest, azi = geod_direct(q2, azi1, dist_q2r1, proj1)
    @test dest ≈ r1 atol=1e-6
    @test azi ≈ azi2 atol=1e-6
    @test geod_destination(q2, azi1, dist_q2r1, proj1) ≈ r1 atol=1e-6

    @test geod_distance(r1, q2, proj2) ≈ dist_r1q2 atol=1e-6
    @test geod_distance(q2, r1, proj1) ≈ dist_q2r1 atol=1e-6
    @test geod_distance(q1, r2, proj1) ≈ dist_q1r2 atol=1e-6
    @test geod_distance(r2, q1, proj1) ≈ dist_r2q1 atol=1e-6

    # The distances computed in both projections can be significantly different if the ellipsoids are different
    @test (dist_r1q2 - dist_q1r2) > 5e4
    @test (dist_q2r1 - dist_r2q1) > 5e4

    # So be consistent with your choice of projections when computing distances,
    # and deviate from WGS84 only if you know what you're doing
    @test geod_distance(r1, q2, proj2) ≈ geod_distance(q2, r1, proj2) atol=1e-6
    @test geod_distance(q1, r2, proj1) ≈ geod_distance(r2, q1, proj1) atol=1e-6
    @test geod_distance(r1, q2, proj2) ≈ geod_distance(q2, r1, proj2) atol=1e-6
    @test geod_distance(q1, r2, proj1) ≈ geod_distance(r2, q1, proj1) atol=1e-6
end

epsg_error = Int[]
for epsg_code in keys(Proj4.epsg)
    try
        proj = Projection(Proj4.epsg[epsg_code])
        proj_string1 = string(proj)
        proj1 = Projection(proj_string1)
        proj_string2 = string(proj1)
        proj2 = Projection(proj_string2)
        proj_string3 = string(proj2)
        @test proj_string1 == proj_string2
        @test proj_string2 == proj_string3
        @test proj_string3 == proj_string1
    catch
        push!(epsg_error, epsg_code)
    end
end

esri_error = Int[]
for esri_code in keys(Proj4.esri)
    try
        proj_string = Proj4.esri[esri_code]
        proj = Projection(proj_string)
        proj_string1 = string(proj)
        proj1 = Projection(proj_string1)
        proj_string2 = string(proj1)
        proj2 = Projection(proj_string2)
        proj_string3 = string(proj2)
        @test proj_string1 == proj_string2
        @test proj_string2 == proj_string3
        @test proj_string3 == proj_string1
    catch
        push!(esri_error, esri_code)
    end
end

if length(epsg_error) > 0 || length(esri_error) > 0
    errorFraction = (length(epsg_error) + length(esri_error)) /
                    (length(Proj4.epsg) + length(Proj4.esri))
    # Some errors are ok (due to old libproj versions), but a good fraction of
    # the strings should parse - if not something *really* wrong has probably
    # occurred.
    @test errorFraction < 0.1

    println(
    """
    The following projection strings could not be parsed by your version of
    libproj Note that this isn't necessarily a problem, but you won't be able
    to use the projections in question.

    total errors: $(round(100*errorFraction; digits=2))%
    """)
    for epsg_code in sort(epsg_error)
        println("[EPSG:$epsg_code] \"$(Proj4.epsg[epsg_code])\"")
    end
    for esri_code in sort(esri_error)
        println("[ESRI:$esri_code] \"$(Proj4.esri[esri_code])\"")
    end
end

end # testset "Proj4"


text:
using Test
import GDAL
import ArchGDAL; const AG = ArchGDAL

@testset "Test methods for raster dataset" begin
    AG.read("data/utmsmall.tif") do dataset
        @testset "Method 1" begin
            AG.copy(dataset, filename = "/vsimem/utmcopy.tif") do copydataset
                @test AG.ngcp(copydataset) == 0
                AG.getband(copydataset,1) do band
                    @test AG.noverview(band) == 0
                    AG.buildoverviews!(copydataset, Cint[2,4,8])
                    @test AG.noverview(band) == 3
                    AG.copywholeraster(dataset, copydataset,
                                       progressfunc=GDAL.gdaltermprogress)
                end
            end
        end
        @testset "Method 2" begin
            copydataset = AG.copy(dataset, filename = "/vsimem/utmcopy.tif")
            @test AG.ngcp(copydataset) == 0
            @test AG.noverview(AG.getband(copydataset,1)) == 0
            AG.buildoverviews!(copydataset, Cint[2,4,8])
            @test AG.noverview(AG.getband(copydataset,1)) == 3
            AG.copywholeraster(dataset, copydataset,
                               progressfunc=GDAL.gdaltermprogress)
        end
        AG.copyfiles("GTiff", "/vsimem/utmcopy2.tif", "/vsimem/utmcopy.tif")
        AG.update("/vsimem/utmcopy2.tif") do copydataset
            @test AG.ngcp(copydataset) == 0
            @test AG.noverview(AG.getband(copydataset,1)) == 3
            AG.copywholeraster(dataset, copydataset, options = ["COMPRESS=LZW"])
        end
    end
end

@testset "Test methods for vector dataset" begin
    dataset1 = AG.read("data/point.geojson")
    @test AG.nlayer(dataset1) == 1
    layer1 = AG.getlayer(dataset1, 0)
    @test AG.nfeature(layer1) == 4
    AG.getlayer(dataset1, 0) do layer1
        @test AG.nfeature(layer1) == 4
    end

    dataset2 = AG.create(AG.getdriver("Memory"))
    @test AG.nlayer(dataset2) == 0
    layer2 = AG.copy(layer1, dataset = dataset2, name = "copy")
    @test AG.nlayer(dataset2) == 1
    @test AG.nfeature(layer2) == 4
    @test AG.getname(layer2) == "copy"

    layer3a = AG.getlayer(dataset2, "copy")
    @test AG.nlayer(dataset2) == 1
    @test AG.nfeature(layer3a) == 4
    @test AG.getname(layer3a) == "copy"
    AG.getlayer(dataset2, "copy") do layer3a
        @test AG.nfeature(layer3a) == 4
        @test AG.getname(layer3a) == "copy"
    end

    layer3b = AG.copy(layer3a)
    @test AG.nlayer(dataset2) == 1 # layer3b is not associated with dataset2
    @test AG.getname(layer3b) == "copy(copy)"

    AG.copy(layer3b, dataset = dataset2)
    @test AG.nlayer(dataset2) == 2
    AG.deletelayer!(dataset2, 1)
    @test AG.nlayer(dataset2) == 1

    dataset4 = AG.create(tempname(), driver = AG.getdriver("KML"))
    @test AG.nlayer(dataset4) == 0
    layer4 = AG.createlayer(
        name = "layer4",
        dataset = dataset4,
        geom = GDAL.wkbLineString
    )
    @test AG.nlayer(dataset4) == 1

    AG.create(tempname(), driver = AG.getdriver("KML")) do dataset5
        @test AG.nlayer(dataset5) == 0
        AG.createlayer(
            name = "layer5",
            dataset = dataset5,
            geom = GDAL.wkbLineString
        ) do layer5
            @test AG.nfeature(layer5) == 0
        end
        @test AG.nlayer(dataset5) == 1
    end

    layer5 = AG.createlayer()
    @test AG.getname(layer5) == ""
    @test AG.nfeature(layer5) == 0
    @test AG.nfield(layer5) == 0
    @test AG.ngeom(layer5) == 1
end


text:
using Documenter
using LogarithmicNumbers

makedocs(
    sitename = "LogarithmicNumbers",
    format = Documenter.HTML(),
    modules = [LogarithmicNumbers],
    pages = ["Home" => "index.md", "Guide" => "guide.md", "Interoperability" => "interop.md"]
)

# Documenter can also automatically deploy documentation to gh-pages.
# See "Hosting Documentation" and deploydocs() in the Documenter manual
# for more information.
#=deploydocs(
    repo = "<repository url>"
)=#


text:
module GCoptimization

using Cxx

const libPath = joinpath(dirname(@__FILE__), "..", "deps", "usr", "lib")
const includePath = joinpath(dirname(@__FILE__), "..", "deps", "usr", "include")

global const libgco = Libdl.find_library(["libGCoptimization", "libGCoptimization-x86", "libGCoptimization-x64"], [libPath])
if isempty(libgco)
    error("GCoptimization not properly installed. Please run Pkg.build(\"GCoptimization\")")
end

function __init__()
    addHeaderDir(libPath, kind=C_System)
    Libdl.dlopen(libgco, Libdl.RTLD_GLOBAL)
    cxxinclude(includePath*"/GCoptimization.h")
end

# types
export GCObool, GCOnode_id, GCOLabelID, GCOEnergyTermType, GCOVarID, GCOSiteID
export GCOSparseDataCost

# low-level APIs
export GCoptimizationGridGraph, setSmoothCostVH,
       GCoptimizationGeneralGraph, setNeighbors

export expansion, alpha_expansion, swap, alpha_beta_swap

export setDataCost, setSmoothCost, setLabelCost, setLabelSubsetCost, whatLabel,
       setLabel, setLabelOrder, compute_energy, giveDataEnergy, giveSmoothEnergy,
       giveLabelEnergy, numSites, numLabels, setVerbosity

# high-level APIs
export gco_αexpansion, gco_αβswap

export gco_create, gco_setneighbors, gco_setdatacost, gco_setsmoothcost,
       gco_setlabelcost, gco_expansion, gco_swap, gco_energy, gco_getlabeling,
       gco_setlabeling, gco_setlabelorder

# types
include("types.jl")

# low-level APIs
include("lowlevel.jl")

# core functionalities
include("core.jl")

# high-level APIs
include("highlevel.jl")

end # module


text:
abstract type Chronology end

struct Synchronize <: Chronology
    from_horizon::Int64 #number of time periods to grab data from
    to_steps::Int64 #number of times to run using the same data
end

struct RecedingHorizon <: Chronology end

struct Sequential <: Chronology end


text:
using MaxMinFilters
using Test
using Random

len = 1000000
Random.seed!(0)
T = Float64
vec = rand(T, len)
window1 = 5
window2 = 200

mx_ref = similar(vec)
mn_ref = similar(vec)

mx = similar(vec)
mn = similar(vec)
using BenchmarkTools
deck1 = @benchmark movmaxmin!($mx, $mn, $vec, $window1)
cbuf1 = @benchmark movmaxmin_!($mx, $mn, $vec, $window1)
inli1 = @benchmark _movmaxmin!($mx,$mn, $vec, $window1)
deck2 = @benchmark movmaxmin!($mx, $mn, $vec, $window1)
cbuf2 = @benchmark movmaxmin_!($mx,$mn, $vec, $window1)
inli2 = @benchmark _movmaxmin!($mx,$mn, $vec, $window1)

t21 = @benchmark    _movmaxmin!($mx,$mn, $vec, $window2)
t22 = @benchmark  ie_movmaxmin!($mx,$mn, $vec, $window2)
t22 = @benchmark ien_movmaxmin!($mx,$mn, $vec, $window2)
t21 = @benchmark    _movmaxmin!($mx,$mn, $vec, $window2)
t22 = @benchmark  ie_movmaxmin!($mx,$mn, $vec, $window2)
t22 = @benchmark ien_movmaxmin!($mx,$mn, $vec, $window2)

t21 = @benchmark    _movmaxmin!($mx,$mn, $vec, $window2)
t22 = @benchmark  ie_movmaxmin!($mx,$mn, $vec, $window2)
t22 = @benchmark ien_movmaxmin!($mx,$mn, $vec, $window2)
t21 = @benchmark    _movmaxmin!($mx,$mn, $vec, $window2)
t22 = @benchmark  ie_movmaxmin!($mx,$mn, $vec, $window2)
t22 = @benchmark ien_movmaxmin!($mx,$mn, $vec, $window2)

e = 1


text:
function haberman(m::D, s::AbstractMatrix{T}, f::AbstractVector{D},
                  β0::AbstractVector{T} = zeros(T, size(s)[2]*(m-1)),
                  method::Optim.SecondOrderOptimizer = Optim.Newton(),
                  options=Optim.Options()) where {T <: Real, D <: Integer}

  Q = size(s)[2]
  δ = zeros(T, m, Q)
  γ = zeros(T, m+1, m+1)
  η = zeros(T, m, m, m+1)

  pr = zeros(T, m, m+1)

  πs = zeros(T, m, m+1)
  πd = zeros(T, m, m, m+1)

  tmp = zeros(T, m-1)
  last_β = rand(T, Q*(m-1))
  function calculate_common!(x, last_x)
    if x != last_x
      copyto!(last_x,x)
      @inbounds for q in 1:Q
        δ[1,q] = 1
        δ[2:end,q] .= exp.(x[((q-1)*(m-1)+1):(q*(m-1))])
      end

      @inbounds for r in 0:m
        pr[1,r+1] = 0.5
        tmp .= 0.0
        for q in 1:Q
          @views tmp .+= r^(q-1) .* x[((q-1)*(m-1)+1):(q*(m-1))]
        end
        @views pr[2:end,r+1] .= 1 ./ (1 .+ exp.(-tmp))
        @views poisbin_sum_taub!(γ[:,r+1], pr[:,r+1])
        @views poisbin_sum_taub_dervs_2!(η, pr[:,r+1])

        for i in 1:m
          πs[i,r+1] = r == 0 ? 0.0 : pr[i,r+1] * η[i,i,r] / γ[r+1,r+1]
          for j in 1:m
            if r == 0
              πd[i,j,r+1] = 0.0
            elseif r == 1
              πd[i,j,r+1] = -pr[i,r+1] * pr[j,r+1] * f[r] * η[i,i,r] * η[j,j,r] / γ[r+1,r+1]^2
            else
              πd[i,j,r+1] = pr[i,r+1] * pr[j,r+1] * f[r] * (η[i,j,r-1] / γ[r+1,r+1] - η[i,i,r] * η[j,j,r] / γ[r+1,r+1]^2)
            end
          end
        end
      end
    end
  end

  function neglogLC(β)
    calculate_common!(β, last_β)
    result = 0.0
    @inbounds for q in 1:Q
      @views result -= s[:,q]'log.(δ[:,q])
    end
    @inbounds for r in 1:m
      result += f[r] * (log(γ[r+1,r+1]) - log(γ[1,r+1]))
    end
    return result
  end

  function neg_g!(storage, β)
    calculate_common!(β, last_β)
    @inbounds for i in 1:(m-1)
      for q in 1:Q
        storage[i + (q-1)*(m-1)] = -s[i+1,q]
      end
      for r in 1:m
        A = f[r] * πs[i+1,r+1]
        for q in 1:Q
          storage[i + (q-1)*(m-1)] += r^(q-1) * A
        end
      end
    end
  end

  function neg_h!(storage, β)
    calculate_common!(β, last_β)
    storage .= 0.0
    @inbounds for r in 1:m
      for i in 1:(m-1)
        A = f[r] * (πs[i+1,r+1] - πs[i+1,r+1]^2)
        for q1 in 1:Q
          for q2 in 1:Q
            storage[i + (q1-1)*(m-1),i + (q2-1)*(m-1)] += r^(q1-1 + q2-1) * A
          end
        end
        for j in (i+1):(m-1)
          B = πd[i+1,j+1,r+1]
          for q1 in 1:Q
            for q2 in 1:Q
              storage[i + (q1-1)*(m-1),j + (q2-1)*(m-1)] += r^(q1-1 + q2-1) * B
              storage[j + (q2-1)*(m-1),i + (q1-1)*(m-1)] += r^(q1-1 + q2-1) * B
            end
          end
        end
      end
    end
  end

  df = TwiceDifferentiable(neglogLC, neg_g!, neg_h!, β0)
  res = Optim.optimize(df, β0, method, options)
  delta_hat = zeros(m, Q)

  @inbounds for q in 1:Q
    delta_hat[1,q] = 1.0
    delta_hat[2:m,q] .= exp.(res.minimizer[((q-1)*(m-1)+1):(q*(m-1))])
  end

  return (res, delta_hat, πs)
end


text:
module MultiCellArrays

using Gridap
using Gridap.CachedArrays
using Gridap.CellArrayApply: CellArrayFromKernel

export MultiCellArray
export MultiCellMatrix
export MultiCellVector
import Base: iterate
import Base: length
import Gridap: reindex

struct MultiCellArray{T,N}
  blocks::Vector{<:CellValue{CachedArray{T,N,Array{T,N}}}}
  fieldids::Vector{NTuple{N,Int}}
end

const MultiCellMatrix{T} = MultiCellArray{T,2}

const MultiCellVector{T} = MultiCellArray{T,1}

function MultiCellArray(cellarrays::Vector{<:CellArray{T,N}},fieldids::Vector{NTuple{N,Int}}) where {T,N}
  @assert length(cellarrays) > 0
  @assert all( [ length(ca) == length(cellarrays[1]) for ca in cellarrays ])
  _cellarrays = [ _prepare(ca) for ca in cellarrays]
  MultiCellArray{T,N}(_cellarrays,fieldids)
end

function MultiCellMatrix(cellmatrices::Vector{<:CellMatrix{T}},fieldids::Vector{NTuple{2,Int}}) where T
  MultiCellArray(cellmatrices,fieldids)
end

function MultiCellVector(cellvectors::Vector{<:CellVector{T}},fieldids::Vector{NTuple{1,Int}}) where T
  MultiCellArray(cellvectors,fieldids)
end

function _prepare(ca::CellValue{CachedArray{T,N,Array{T,N}}}) where {T,N}
  ca
end

function _prepare(ca::CellArray{T,N}) where {T,N}
  k = ArrayKernelFromBroadcastedFunction(+)
  apply(k,ca)
end

function _prepare(ca::ConstantCellArray{T,N}) where {T,N}
  _ca = CachedArray(ca.value)
  ConstantCellValue(_ca,length(ca))
end

length(mca::MultiCellArray) = length(mca.blocks[1])

@inline function iterate(mca::MultiCellArray)
  zipped = zip(mca.blocks...)
  znext = iterate(zipped)
  if znext === nothing; return nothing end
  arrays, zstate = znext
  state = (zstate, zipped)
  (arrays, state)
end

@inline function iterate(mca::MultiCellArray,state)
  zstate, zipped = state
  znext = iterate(zipped, zstate)
  if znext === nothing; return nothing end
  arrays, zstate = znext
  state = (zstate, zipped)
  (arrays, state)
end

function reindex(mca::MultiCellArray,cn::CellNumber)
  blocks = [ reindex(block,cn)  for block in mca.blocks ]
  MultiCellArray(blocks,mca.fieldids)
end

end # module


text:
using NLOptControl, JuMP, Parameters, VehicleModels
include("plots.jl");include("autonomousControl.jl");

pa=Vpara(u0_= 15.); @unpack_Vpara pa # other vehicles parameters

c=defineCase("c3");

# initialize
n1 = NLOpt();s1=Settings(;MPC=false,save=true,simulate=false,format=:png);

# define
pa1 = VparaKB();  @unpack_VparaKB pa1 # vehicle parameters  TODO do I need these?
X0 = [x0_,y0_,psi0_,u0_]; #TODO pull these states in
XF = [NaN,NaN,NaN,NaN];

#TODO when running in parallel run several with various constraints on x_min
XL = [x_min,y_min,psi_min,u_min];
XU = [x_max,y_max,psi_max,u_max];
CL = [sa_min,ax_min];
CU = [sa_max,ax_max];
n1 = define(n1,stateEquations=KinematicBicycle,numStates=4,numControls=2,X0=X0,XF=XF,XL=XL,XU=XU,CL=CL,CU=CU)

# build
n1 = configure(n1,Ni=2,Nck=[15,10];(:integrationMethod => :ps),(:integrationScheme => :lgrExplicit),(:finalTimeDV => true))
defineSolver(n1,solver=:IPOPT)
mdl1 = build(n1);

# addtional information
names = [:x,:y,:psi,:ux];
descriptions = ["X (m)","Y (m)","Yaw Angle (rad)","Longitudinal Velocity (m/s)"];
stateNames(n1,names,descriptions)
names = [:sr,:jx];
descriptions = ["Steering Angle (rad)","Longitudinal Acceleration (m/s^2)"];
controlNames(n1,names,descriptions);

#mXL=Any[false,false,false,false];mXU=Any[false,false,false,-1];  # set to false if you don't want to taper that side
#linearStateTolerances(n;mXL=mXL,mXU=mXU);

# setup OCP
params1 = [pa1];   # vehicle parameters
n1,r1 = OCPdef(mdl1,n1,s1,params1);
x_ref = 200; y_ref = 100; # define target TODO get from other model
@NLobjective(mdl1, Min,  n1.tf + (r1.x[end,1]-x_ref)^2 + (r1.x[end,2]-y_ref)^2);

# obstacles
Q = size(c.A)[1]; # number of obstacles
@NLparameter(mdl1, a[i=1:Q] == c.A[i]);
@NLparameter(mdl1, b[i=1:Q] == c.B[i]);
@NLparameter(mdl1, X_0[i=1:Q] == c.X0_obs[i]);
@NLparameter(mdl1, Y_0[i=1:Q] == c.Y0_obs[i]);
@NLparameter(mdl1, speed_x[i=1:Q] == c.s_x[i]);
@NLparameter(mdl1, speed_y[i=1:Q] == c.s_y[i]);

# obstacle postion after the intial postion
X_obs=@NLexpression(mdl1, [j=1:Q,i=1:n1.numStatePoints], X_0[j] + speed_x[j]*n1.tV[i]);
Y_obs=@NLexpression(mdl1, [j=1:Q,i=1:n1.numStatePoints], Y_0[j] + speed_y[j]*n1.tV[i]);

# constraint on position
obs_con=@NLconstraint(mdl1, [j=1:Q,i=1:n1.numStatePoints-1], 1 <= ((r1.x[(i+1),1]-X_obs[j,i])^2)/((a[j]+sm)^2) + ((r1.x[(i+1),2]-Y_obs[j,i])^2)/((b[j]+sm)^2));
newConstraint(r1,obs_con,:obs_con);

# solve
optimize(mdl1,n1,r1,s1)

# post process
r1.results_dir = string(main_dir,"/results_1","/test_1_",n1.solver,"/")
resultsDir(r1.results_dir);
using PrettyPlots, Plots
pyplot();
allPlots(n1,r1,s1,1)

idx=1;
pp=statePlot(n1,r1,s1,idx,1,2;(:lims=>false));
pp=obstaclePlot(n1,r1,s1,c,idx,pp;(:append=>true)); # add obstacles
pp=vehiclePlot(n1,r1,s1,c,idx,pp;(:append=>true)); # add the vehicle


text:
"""
    Quaternion{T<:Real} <: Number

Quaternionic number type with elements of type `T`.

`QuaternionF16`, `QuaternionF32` and `QuaternionF64` are aliases for `Quaternion{Float16}`,
`Quaternion{Float32}` and `Quaternion{Float64}` respectively.  See also [`Rotor`](@ref) and
[`QuatVec`](@ref).

The functions

    Quaternion(w, x, y, z)
    Quaternion(x, y, z)
    Quaternion(w)
    Quaternion{T}(w, x, y, z)

create a new quaternion with the given components.  The argument `w` is the
scalar component, and `x`, `y`, and `z` are the corresponding "vector"
components.  If any of these arguments is missing, it will be set to zero.  The
type of the returned quaternion will be inferred from the input arguments, or
can be specified, by passing the type parameter `T` as above.

Note that the constants [`imx`](@ref), [`imy`](@ref), and [`imz`](@ref) can
also be used like the complex `im` to create new `Quaternion` object.

# Examples

```jldoctest
julia> Quaternion(1, 2, 3, 4)
1 + 2𝐢 + 3𝐣 + 4𝐤
julia> Quaternion{Float64}(1, 2, 3, 4)
1.0 + 2.0𝐢 + 3.0𝐣 + 4.0𝐤
julia> Quaternion(1.0, 2.0, 3.0, 4.0)
1.0 + 2.0𝐢 + 3.0𝐣 + 4.0𝐤
julia> Quaternion(2, 3, 4)
0 + 2𝐢 + 3𝐣 + 4𝐤
julia> Quaternion(1)
1 + 0𝐢 + 0𝐣 + 0𝐤
```
"""
struct Quaternion{T<:Real} <: AbstractQuaternion{T}
    components::SVector{4, T}
    Quaternion{T}(a::SVector{4, T}) where {T<:Real} = new{T}(a)
    Quaternion{T}(a::A) where {T<:Real, A<:AbstractArray} = new{T}(SVector{4, T}(a))
end

@doc raw"""
    Rotor{T<:Real} <: Number

Quaternion of unit magnitude with elements of type `T`.  These objects can be
significantly faster *and* more accurate in certain operations representing
rotations.

A rotor is typically considered to be an element of the group
``\mathrm{Spin}(3) ≃ \mathrm{SU}(2)``, which can be thought of as the subgroup
of quaternions with norm 1.  They are particularly useful as representations of
rotations because a rotor ``R`` acts on a vector ``\vec{v}`` by "conjugation" as

```math
\vec{v}' = R\, \vec{v}\, R^{-1}.
```

This preserves the inner product between any two vectors conjugated in this
way, and so is a rotation.  Note that, because there are two factors of ``R``
here, the sign of ``R`` does not affect the result.  Therefore,
``\mathrm{Spin}(3)`` forms a *double* cover of the rotation group
``\mathrm{SO}(3)``.  For this reason, it will occasionally be useful to
disregard or arbitrarily change the sign of a `Rotor` (as in [`distance`](@ref)
functions) — though this is not generally the default, and may cause problems
if the input rotors change sign when the corresponding rotations are not so
different (cf. [`unflip`](@ref)).

`RotorF16`, `RotorF32` and `RotorF64` are aliases for `Rotor{Float16}`,
`Rotor{Float32}` and `Rotor{Float64}` respectively.  See also
[`Quaternion`](@ref) and [`QuatVec`](@ref).

The functions

    Rotor(w, x, y, z)
    Rotor(x, y, z)
    Rotor(w)
    Rotor{T}(w, x, y, z)

create a new rotor with the given components (where the components are as
described in [`Quaternion`](@ref)), automatically normalizing them on input.
If you would like to bypass this normalization step, you can call

    Rotor{T}(v)

where `v<:AbstractArray`, and can be converted to an `SVector{4, T}`.

However, once a `Rotor` is created, its norm will always be assumed to be
precisely 1.  So if its true norm is significantly different, you will like see
weird results — including vectors with very different lengths after "rotation"
by a non-unit `Rotor`.

Note that simply creating a `Quaternion` that happens to have norm 1 does not
make it a `Rotor`.  However, you can pass such a `Quaternion` to the `Rotor`
function and get the desired result.

# Examples

```jldoctest
julia> Rotor(1, 2, 3, 4)
0.18257418583505536 + 0.3651483716701107𝐢 + 0.5477225575051661𝐣 + 0.7302967433402214𝐤
julia> Rotor(Quaternion(1, 2, 3, 4))
0.18257418583505536 + 0.3651483716701107𝐢 + 0.5477225575051661𝐣 + 0.7302967433402214𝐤
julia> Rotor{Float16}(1, 2, 3, 4)
0.1826 + 0.3652𝐢 + 0.548𝐣 + 0.7305𝐤
julia> Rotor(2, 3, 4)
0.0 + 0.3713906763541037𝐢 + 0.5570860145311556𝐣 + 0.7427813527082074𝐤
julia> Rotor(1)
1 + 0𝐢 + 0𝐣 + 0𝐤
```
"""
struct Rotor{T<:Real} <: AbstractQuaternion{T}
    components::SVector{4, T}
    Rotor{T}(a::SVector{4, T}) where {T<:Real} = new{T}(a)
    Rotor{T}(a::A) where {T<:Real, A<:AbstractArray} = new{T}(SVector{4, T}(a))
end

"""
    QuatVec{T<:Real} <: Number

Pure-vector quaternion with elements of type `T`.  These objects can be
significantly faster *and* more accurate in certain operations than general
`Quaternion`s.

`QuatVecF16`, `QuatVecF32` and `QuatVecF64` are aliases for `QuatVec{Float16}`,
`QuatVec{Float32}` and `QuatVec{Float64}` respectively.  See also
[`Quaternion`](@ref) and [`Rotor`](@ref).

The functions

    QuatVec(w, x, y, z)
    QuatVec(x, y, z)
    QuatVec(w)
    QuatVec{T}(w, x, y, z)

create a new rotor with the given components (where the components are as
described in [`Quaternion`](@ref)), except that the scalar argument `w` is
always set to 0.

# Examples

```jldoctest
julia> QuatVec(1, 2, 3, 4)
0 + 2𝐢 + 3𝐣 + 4𝐤
julia> QuatVec(Quaternion(1, 2, 3, 4))
0 + 2𝐢 + 3𝐣 + 4𝐤
julia> QuatVec(2, 3, 4)
0 + 2𝐢 + 3𝐣 + 4𝐤
julia> QuatVec(1)
0 + 0𝐢 + 0𝐣 + 0𝐤
```
"""
struct QuatVec{T<:Real} <: AbstractQuaternion{T}
    components::SVector{4, T}
    QuatVec{T}(a::SVector{4, T}) where {T<:Real} = new{T}(a)
    QuatVec{T}(a::A) where {T<:Real, A<:AbstractArray} = new{T}(SVector{4, T}(a))
end

# Untyped constructor from SVector
# (::Type{QT})(v::SVector{4, T}) where {T<:Real, QT<:AbstractQuaternion} = QT{eltype(v)}(v)
Quaternion(v::SVector{4, T}) where {T<:Real} = Quaternion{eltype(v)}(v)
Rotor(v::SVector{4, T}) where {T<:Real} = Rotor{eltype(v)}(v)
QuatVec(v::SVector{4, T}) where {T<:Real} = QuatVec{eltype(v)}(v)

# Constructor from all 4 components
(::Type{QT})(w, x, y, z) where {QT<:AbstractQuaternion} = (v=SVector{4}(w, x, y, z); QT{eltype(v)}(v))
(::Type{QT})(w, x, y, z) where {T<:Real, QT<:AbstractQuaternion{T}} = QT(SVector{4, T}(w, x, y, z))
Rotor(w, x, y, z) = (n=√(w^2+x^2+y^2+z^2); Rotor(SVector{4}(w/n, x/n, y/n, z/n)))
Rotor{T}(w, x, y, z) where {T<:Real} = (n=√T(w^2+x^2+y^2+z^2); Rotor{T}(SVector{4, T}(w/n, x/n, y/n, z/n)))
QuatVec(w, x, y, z) = QuatVec(SVector{4}(oftype(w, false), x, y, z))
QuatVec{T}(w, x, y, z) where {T<:Real} = QuatVec{T}(SVector{4, T}(oftype(w, false), x, y, z))

# Constructor from vector components
(::Type{QT})(x, y, z) where {QT<:AbstractQuaternion} = (v=SVector{4}(false, x, y, z); QT{eltype(v)}(v))
(::Type{QT})(x, y, z) where {T<:Real, QT<:AbstractQuaternion{T}} = QT(SVector{4, T}(false, x, y, z))
Rotor(x, y, z) = (n=√(x^2+y^2+z^2); Rotor(SVector{4}(false, x/n, y/n, z/n)))
Rotor{T}(x, y, z) where {T<:Real} = (n=√T(x^2+y^2+z^2); Rotor{T}(SVector{4, T}(false, x/n, y/n, z/n)))
QuatVec(x, y, z) = QuatVec(SVector{4}(false, x, y, z))
QuatVec{T}(x, y, z) where {T<:Real} = QuatVec{T}(SVector{4, T}(false, x, y, z))

# Constructor from scalar component
# (::Type{QT})(w::Real) where {QT<:AbstractQuaternion} = (v=SVector{4}(w, false, false, false); QT{eltype(v)}(v))
# (::Type{QT})(w::Real) where {T<:Real, QT<:AbstractQuaternion{T}} = QT(SVector{4, T}(w, false, false, false))
Quaternion(w::Real) = Quaternion(SVector{4}(w, false, false, false))
Quaternion(w::Symbolics.Num) = Quaternion(SVector{4}(w, false, false, false))
Quaternion{T}(w::Real) where {T<:Real} = Quaternion{T}(SVector{4, T}(w, false, false, false))
Rotor(w::Real) = Rotor(SVector{4}(one(w), false, false, false))
Rotor(w::Symbolics.Num) = Rotor(SVector{4}(one(w), false, false, false))
Rotor{T}(w::Real) where {T<:Real} = Rotor{T}(SVector{4, T}(one(T), false, false, false))
QuatVec(w::Real) = QuatVec(SVector{4, typeof(w)}(false, false, false, false))
QuatVec(w::Symbolics.Num) = QuatVec(SVector{4, typeof(w)}(false, false, false, false))
QuatVec{T}(w::Real) where {T<:Real} = QuatVec{T}(SVector{4, T}(false, false, false, false))

# Copy constructor
# (::Type{QT})(q::AbstractQuaternion) where {QT<:AbstractQuaternion} = QT{eltype(q.components)}(q.components)
# (::Type{QT})(q::AbstractQuaternion{S}) where {T<:Real, S<:Real, QT<:AbstractQuaternion{T}} = QT(SVector{4, T}(q.components))
# (::Type{QT})(q::AbstractQuaternion{T}) where {T<:Real, QT<:AbstractQuaternion{T}} = QT(SVector{4, T}(q.components))
Quaternion(q::AbstractQuaternion{T}) where {T<:Real} = Quaternion(q.components...)
Quaternion{T}(q::AbstractQuaternion{S}) where {T<:Real, S<:Real} = Quaternion{T}(q.components...)
Rotor(q::QT) where {T<:Real, QT<:AbstractQuaternion{T}} = Rotor(q.components...)
Rotor{T}(q::AbstractQuaternion{S}) where {T<:Real, S<:Real} = Rotor{T}(q.components...)
QuatVec(q::AbstractQuaternion{T}) where {T<:Real} = QuatVec(q.components...)
QuatVec{T}(q::AbstractQuaternion{S}) where {T<:Real, S<:Real} = QuatVec{T}(q.components...)

# Type constructors
(::Type{QT})(::Type{T}) where {T<:Real, QT<:AbstractQuaternion} = QT{T}
(::Type{QT})(::Type{<:AbstractQuaternion{T}}) where {T<:Real, QT<:AbstractQuaternion} = QT{T}

# Handy aliases like `ComplexF64`, etc.
const QuaternionF64 = Quaternion{Float64}
const QuaternionF32 = Quaternion{Float32}
const QuaternionF16 = Quaternion{Float16}
const RotorF64 = Rotor{Float64}
const RotorF32 = Rotor{Float32}
const RotorF16 = Rotor{Float16}
const QuatVecF64 = QuatVec{Float64}
const QuatVecF32 = QuatVec{Float32}
const QuatVecF16 = QuatVec{Float16}

# Handy constants like `im`
"""
    imx

The quaternionic unit associated with rotation about the `x` axis.  Can also be entered as Unicode
bold: `𝐢`.

# Examples
```jldoctest
julia> imx * imx
-1 + 0𝐢 + 0𝐣 + 0𝐤
julia> 1.2imx
0.0 + 1.2𝐢 + 0.0𝐣 + 0.0𝐤
```
"""
const imx = QuatVec(false, true, false, false)
const 𝐢 = imx

"""
    imy

The quaternionic unit associated with rotation about the `y` axis.  Can also be entered as Unicode
bold: `𝐣`.

# Examples
```jldoctest
julia> imy * imy
-1 + 0𝐢 + 0𝐣 + 0𝐤
julia> 1.2imy
0.0 + 0.0𝐢 + 1.2𝐣 + 0.0𝐤
```
"""
const imy = QuatVec(false, false, true, false)
const 𝐣 = imy

"""
    imz

The quaternionic unit associated with rotation about the `z` axis.  Can also be entered as Unicode
bold: `𝐤`.

# Examples
```jldoctest
julia> imz * imz
-1 + 0𝐢 + 0𝐣 + 0𝐤
julia> 1.2imz
0.0 + 0.0𝐢 + 0.0𝐣 + 1.2𝐤
```
"""
const imz = QuatVec(false, false, false, true)
const 𝐤 = imz

# Essential constructors
Base.zero(::Type{QT}) where {T<:Real, QT<:AbstractQuaternion{T}} = QT(false, false, false, false)
Base.zero(q::QT) where {T<:Real, QT<:AbstractQuaternion{T}} = Base.zero(QT)
Base.zero(::Type{Rotor}) = throw(DomainError("Rotor", "Zero is not a possible rotor."))
Base.zero(::Type{Rotor{T}}) where T = throw(DomainError("Rotor", "Zero is not a possible rotor."))

Base.one(::Type{QT}) where {T<:Real, QT<:AbstractQuaternion{T}} = QT(true, false, false, false)
Base.one(q::QT) where {T<:Real, QT<:AbstractQuaternion{T}} = Base.one(QT)
Base.one(::Type{QuatVec}) = throw(DomainError("QuatVec", "One is not a possible 3-vector."))
Base.one(::Type{QuatVec{T}}) where T = throw(DomainError("QuatVec", "One is not a possible 3-vector."))

# Getting pieces of quaternions
@inline function Base.getproperty(q::AbstractQuaternion, sym::Symbol)
    @inbounds begin
        if sym === :w
            return q.components[1]
        elseif sym === :x
            return q.components[2]
        elseif sym === :y
            return q.components[3]
        elseif sym === :z
            return q.components[4]
        elseif sym === :re
            return q.components[1]
        elseif sym === :im
            return q.components[2:4]
        elseif sym === :vec
            return q.components[2:4]
        else # fallback to getfield
            return getfield(q, sym)
        end
    end
end
@inline Base.getindex(q::AbstractQuaternion, i::Int) = (@boundscheck checkbounds(q.components,i); q.components[i])
Base.@propagate_inbounds Base.getindex(q::AbstractQuaternion, I) = [q[i] for i in I]
Base.real(::Type{QT}) where {T<:Real, QT<:AbstractQuaternion{T}} = real(T)
Base.real(q::AbstractQuaternion{T}) where {T<:Real} = q.re
Base.imag(q::AbstractQuaternion{T}) where {T<:Real} = q.im

# Type games
wrapper(::T) where {T} = wrapper(T)
wrapper(T::UnionAll) = T
wrapper(T::Type{Q}) where {S<:Real, Q<:AbstractQuaternion{S}} = T.name.wrapper
wrapper(::Type{T}, ::Type{T}) where {T<:AbstractQuaternion} = wrapper(T)
wrapper(::Type{<:AbstractQuaternion}, ::Type{<:AbstractQuaternion}) = Quaternion

wrapper(::Type{<:AbstractQuaternion}, ::Val{OP}, ::Type{<:AbstractQuaternion}) where {OP} = Quaternion
wrapper(::Type{<:AbstractQuaternion}, ::Val{OP}, ::Type{<:Real}) where {OP} = Quaternion
wrapper(::Type{<:Real}, ::Val{OP}, ::Type{<:AbstractQuaternion}) where {OP} = Quaternion
wrapper(::Type{<:AbstractQuaternion}, ::Val{OP}, ::Type{<:Symbolics.Num}) where {OP} = Quaternion
wrapper(::Type{<:Symbolics.Num}, ::Val{OP}, ::Type{<:AbstractQuaternion}) where {OP} = Quaternion

wrapper(::Type{<:Rotor}, ::Val{*}, ::Type{<:Rotor}) = Rotor
wrapper(::Type{<:Rotor}, ::Val{/}, ::Type{<:Rotor}) = Rotor

wrapper(::Type{<:QuatVec}, ::Val{+}, ::Type{<:QuatVec}) = QuatVec
wrapper(::Type{<:QuatVec}, ::Val{-}, ::Type{<:QuatVec}) = QuatVec
wrapper(::Type{<:QuatVec}, ::Val{*}, ::Type{<:Real}) = QuatVec
wrapper(::Type{<:QuatVec}, ::Val{/}, ::Type{<:Real}) = QuatVec
wrapper(::Type{<:Real}, ::Val{*}, ::Type{<:QuatVec}) = QuatVec
wrapper(::Type{<:Real}, ::Val{/}, ::Type{<:QuatVec}) = QuatVec
wrapper(::Type{<:QuatVec}, ::Val{*}, ::Type{<:Symbolics.Num}) = QuatVec
wrapper(::Type{<:QuatVec}, ::Val{/}, ::Type{<:Symbolics.Num}) = QuatVec
wrapper(::Type{<:Symbolics.Num}, ::Val{*}, ::Type{<:QuatVec}) = QuatVec
wrapper(::Type{<:Symbolics.Num}, ::Val{/}, ::Type{<:QuatVec}) = QuatVec

Base.eltype(::Type{<:AbstractQuaternion{T}}) where {T} = T
Base.widen(::Type{Q}) where {Q<:AbstractQuaternion} = wrapper(Q){widen(eltype(Q))}
Base.float(::Type{Q}) where {Q<:AbstractQuaternion{<:AbstractFloat}} = Q
Base.float(::Type{Q}) where {Q<:AbstractQuaternion} = wrapper(Q){float(eltype(Q))}
Base.float(q::AbstractQuaternion{T}) where {T<:AbstractFloat} = q
Base.float(q::AbstractQuaternion{T}) where {T} = wrapper(q){float(T)}(float(q.components))

Base.big(::Type{Q}) where {Q<:AbstractQuaternion} = wrapper(Q){big(eltype(Q))}
Base.big(q::AbstractQuaternion{T}) where {T<:Real} = wrapper(q){big(T)}(q)

Base.promote_rule(::Type{Q}, ::Type{S}) where {Q<:AbstractQuaternion, S<:Real} =
    wrapper(Q){promote_type(eltype(Q),S)}
Base.promote_rule(::Type{QuatVec{T}}, ::Type{S}) where {T<:Real, S<:Real} =
    Quaternion{promote_type(T,S)}
Base.promote_rule(::Type{Q1}, ::Type{Q2}) where {Q1<:AbstractQuaternion, Q2<:AbstractQuaternion} =
    wrapper(wrapper(Q1), wrapper(Q2)){promote_type(eltype(Q1),eltype(Q2))}


text:
@doc raw"""
    f_diff_weight(k, j)

Weight coefficient
```math
c_{j}^{k}=(-1)^{j}\binom{k}{j},
```
of the ``k^{th}``-order finite difference operator ``\nabla^k`` and corresponding to the function value ``f[n-j]``.
#### Example:
```
k = 5; j = 3
f_diff_weight(k, j)
 -10
```
"""
f_diff_weight(k::Int, j::Int) = Base.iseven(j) ? Base.binomial(k,j) : -Base.binomial(k,j)

# ==============================================================================

@doc raw"""
    f_diff_weights(k)

Weight vector ``[c_k^k,\ \ldots,\ c_0^k]`` defining the ``k^{th}``-order finite difference operator,
```math
\nabla^{k}f[n]	=[c_{k}^{k},\thinspace c_{k-1}^{k},\thinspace\ldots,c_{0}^{k}]\left[\begin{array}{c}
f[n-k]\\
\vdots\\
f[n]
\end{array}\right]=\sum_{j=0}^{k}c^k[j]f[n-k+j],
```
where ``c^k[j] ≡ c_{k-j}^{k}``.
This form is designed for use with *analytic* functions, ``f``, tabulated in *forward* order, ``f[n-k], ...,f[n]``.
#### Example:
```
k = 3
f_diff_weights(k)
4-element Vector{Int64}:
  1
 -3
  3
 -1
```
"""
f_diff_weights(k::Int) = [CamiXon.f_diff_weight(k, k-j) for j=0:k]

# ==============================================================================

@doc raw"""
    f_diff_weights_array(kmax)

Collection of weight vectors, ``c^0,\ \ldots,\ c^k``, defining the finite difference operators ``\nabla^0,\ \ldots,\ \nabla^k``
under the convention ``c^k ≡ [c_k^k,\ c_1^k,\ldots,\ c_0^k])``.
#### Example:
```
kmax = 3
∇ = f_diff_weights_array(kmax)
4-element Vector{Vector{Int64}}:
 [1]
 [-1, 1]
 [1, -2, 1]
 [-1, 3, -3, 1]
```
"""
f_diff_weights_array(kmax::Int) = [CamiXon.f_diff_weights(k)  for k=0:kmax]

# ==============================================================================

@doc raw"""
    f_diff_expansion_weights(a, ∇)

Weight vector ``b^k ≡ [b_k^k,\ ,\ldots,\ b_0^k]`` corresponding to the expansion coefficients
``[a_0^k,\ ,\ldots,\ a_k^k]`` of the ``k^{th}``-order finite-difference expansion,

```math
\sum_{p=0}^{k}a_{p}\nabla^{p}f[n]=\sum_{j=0}^{k}b^k[j]f[n-k+j],
```

where ``b^k[j] \equiv b_{k-j}^k`` and ``f[n-k], ...,f[n]`` are elements of the analytic function ``f`` tabulated in *forward* order.
Note the difference in ordering between the finite-difference expansion *coefficients*,
``a_{0},\ \ldots,\ a_{k}``, and the finite-difference expansion *weights*, ``b_k^{k},\ \ldots,\ b_0^{k}``.
Note further the difference in ``k`` dependence: the *weights*,
``b_j^k``, are ``k``*-dependent*, whereas the *coefficients*, ``a_j``, are not.
#### Example:
```
k=5
∇ = f_diff_weights_array(k)
a = UnitRange(0,k)
b = f_diff_expansion_weights(a, ∇)
6-element Vector{Int64}:
  -5
  29
 -69
  85
 -55
  15
```
"""
function f_diff_expansion_weights(coeffs, ∇)
# ======================================================================================
#   function weights of finite-difference summation
# ======================================================================================
    k = Base.length(coeffs)-1
    #return [sum([coeffs[1+p] * ∇[1+p][1+p-j] for p=j:k]) for j=0:k]
    return [Base.sum([coeffs[1+p] * ∇[1+p][1+p-j] for p=j:k]) for j=k:-1:0]
end

# ==============================================================================

@doc raw"""
    f_diff_expansion_coeffs_lagrange(k::Int, x::T) where T<:Real

Finite-difference expansion coefficient vector ``[l_0(x),\ \ldots,\ l_p(x)]`` defining
``k^{th}``-order lagrangian interpolation of the tabulated analytic function ``f(n+x)``
at offset position ``x`` with respect to position ``n``,
```math
f[n+x] = (1 - \nabla)^{-x} f[n] \equiv \sum_{p=0}^{\infty}l_p(x)∇^p f[n],
```
where ``l_0\equiv 1`` and ``l_p(x) = x(x+1)(x+2)\cdots(x+p-1)/p!``. Interpolation corresponds to the interval
``-k\le\ x\le 0``; extrapolation to ``x\ge 0``.
#### Examples:
```
k = 5; x = 1
l = f_diff_expansion_coeffs_lagrange(k,x); println(l)
 [1, 1, 1, 1, 1, 1]
```
"""
function f_diff_expansion_coeffs_lagrange(k::Int, x::T) where T<:Real
# ======================================================================================
#   f_difference expansion coefficients for lagrange interpolation for position n+x
# ======================================================================================
    l = Base.ones(T,k+1)
    x == 1 ? l : x ==0 ? (for i=2:k+1; l[i] = 0 end) : (for i=1:k; l[i+1] = l[i]*(x+i-1)/i end)
    return l
end

# ==============================================================================

@doc raw"""
    summation_range(n, i, k, m)

Summation range for interpolation position ``0\le i/m \le 1`` used in ``k^{th}``-order
lagrangian interpolation of the anaytic function ``f`` tabulated in forward
order on a uniform grid of ``n`` points, ``f[1],\ \ldots,\ f[n]``;
``m`` is the multiplier defining the interpolation grid size.

#### Examples:
```
n = 7; k = 2; m = 1
o = [summation_range(n,i,k,m) for i=0:(n-1)*m]; println(o)
 UnitRange{Int64}[1:3, 2:4, 3:5, 4:6, 5:7, 5:7, 5:7]
```
"""
function summation_range(n::Int, i::Int, k::Int, m::Int)
# ================================================================================================
#   summation range for point position i lagrangian interpolation
# ================================================================================================
      0 ≤ i ≤ n*m || error("Error: position index i outside index range 0 ≤ i ≤ n⋅m")
     return i < (n-1-k)*m  ? UnitRange(i÷m+1,i÷m+k+1) : UnitRange(n-k,n)
end

# ==============================================================================

@doc raw"""
    f_diff_function_sequences(f, k::Int, m=1)

Finite-difference summation sequences of function values given in forward order
for use in ``k^{th}``-order lagrangian interpolation of the anaytic function
``f`` tabulated in forward order on a regular grid of ``n`` points, ``f[1], ...,f[n]``;
``m`` is the multiplier defining the interpolation grid size. Each sequence consists of ``k⋅m+1`` function values.
#### Example:
```
f = [0,1,2,3,4,5,6]
k = 2
o = f_diff_function_sequences(f, k); println(o)
 [[0, 1, 2], [1, 2, 3], [2, 3, 4], [3, 4, 5], [4, 5, 6], [4, 5, 6], [4, 5, 6]]
```
"""
function f_diff_function_sequences(f, k::Int, m=1)
# ================================================================================================
#   finite-difference function values for interpolation range of lagrangian interpolation
# ================================================================================================
    n = Base.length(f)
    return [f[CamiXon.summation_range(n,i,k,m)] for i=0:(n-1)*m]
end

# ==============================================================================

@doc raw"""
    lagrangian_interpolation(f::Vector{Float64}, domain::ClosedInterval{Float64}; k=1, m=1)

``k^{th}``-order lagrangian *interpolation* of the analytic function ``f`` tabulated
in forward order on a regular grid of ``n`` points, ``f[1],\ \ldots,\ f[n]``;
``m`` is the multiplier defining the interpolation grid size.
#### Example:
```
f = [0.0,1,2,3,4,5,6,7]
domain = 0.0..1.0
(X,Y) = lagrangian_interpolation(f, domain; k=2, m=2); println((X,Y))
 (0.0:0.07142857142857142:1.0, [0.0, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0, 5.5, 6.0, 6.5, 7.0])
```
"""
function lagrange_interpolation(f::Vector{Float64}, domain::ClosedInterval{Float64}; k=3, m=1)
# ======================================================================================
#   lagrangian (k+1)-point interpolation at i interpolation points
# ======================================================================================
    n = length(f)

    ∇ = CamiXon.f_diff_weights_array(k)
    l = [CamiXon.f_diff_expansion_coeffs_lagrange(k, x) for x=-k:1/m:0]
    w = [CamiXon.f_diff_expansion_weights(l[i], ∇) for i ∈ eachindex(l)]
    w1 = Base.append!(Base.repeat(w[1:m],n-k-1),w)
    w2 = CamiXon.f_diff_function_sequences(f, k, m)

    X = Base.range(domain.left, domain.right, length=(n-1)*m+1)
    Y = [w1[i] ⋅ w2[i] for i ∈ Base.eachindex(w1)]

    return X, Y

end

@doc raw"""
    lagrangian_extrapolation(f::Vector{Float64}, domain::ClosedInterval{Float64}; k=1, e=1, m=1)

``k^{th}``-order lagrangian *extrapolation* up to position ``n+e`` of the analytic function
``f`` tabulated in forward order at ``n`` points, ``f[1],\ \ldots,\ f[n]``;
``m`` is the multiplier defining the interpolation grid size.
#### Example:
```
f = [0.0,1,2,3,4,5,6,7]
domain = 0.0..1.0
(X,Y) = lagrangian_extrapolation(f, domain; k=2, e=1, m=2); println((X,Y))
 (0.0:0.07142857142857142:1.0, [0.0, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0, 5.5, 6.0, 6.5, 7.0])
```
"""
function lagrange_extrapolation(f::Vector{Float64}, domain::ClosedInterval{Float64}; k=1, e=1, m=1)
# ======================================================================================
#   lagrangian (k+1)-point interpolation at μ interpolation points
# ======================================================================================
    n = Base.length(f)

    ∇ = CamiXon.f_diff_weights_array(k)
    l = [CamiXon.f_diff_expansion_coeffs_lagrange(k, x) for x=0:1/m:e]
    w1 = [CamiXon.f_diff_expansion_weights(l[i], ∇) for i ∈ Base.eachindex(l)]
    w2 = CamiXon.f_diff_function_sequences(f, k, m)[end]

    ΔX = (domain.right - domain.left)/((n-1)*m)
    X = Base.range(domain.right, domain.right + ΔX * m*e, length=m*e+1)
    Y = [w1[i] ⋅ w2 for i ∈ Base.eachindex(w1)]

    return X, Y

end

# ===================================== f_diff_expansion_coeffs_differentiation(k, x) ====

@doc raw"""
    f_diff_expansion_coeffs_differentiation(k::Int, x::T) where T<:Real

Finite-difference expansion coefficient vector ``[l_0^{\prime}(x),\ \ldots,\ l_p^{\prime}(x)]`` defining
``k^{th}``-order lagrangian *differentiation* of the tabulated analytic function ``f(n+x)``
at position ``x``,
```math
\frac{df}{dx}[n+x]=\sum_{p=0}^{k}l_{p}^{\prime}(x)\nabla^{p}f[n]
```
#### Example:
```
k = 2; x = 0
o = f_diff_expansion_coeffs_differentiation(k,x); println(o)
 [0.0, 1.0, -1.5]
```
"""
function f_diff_expansion_coeffs_differentiation(k::Int, x::T) where T<:Real
# ======================================================================================
#   finite difference expansion coeffs for differentiation in interval -k ≤ x ≤ 0
# ======================================================================================
    a = Base.prepend!([1//i for i=1:k],[0//1])
    b = CamiXon.f_diff_expansion_coeffs_lagrange(k, x)

    return CamiXon.polynom_product_expansion(a, b, k)

end

# =================== create_lagrange_differentiation_weights(k, x) ============

@doc raw"""
    create_lagrange_differentiation_weights(k::Int, x::T) where T<:Real

``k^{th}``-order Lagrange differentiation weights vector, ``s^k(x) ≡ [s_k^k(x),\ ,\ldots,\ s_0^k(x)]``, where ``x`` is the position relative point ``n``.

```math
\frac{df}{dx}[n+x]= \sum_{j=0}^{k}s_{k-j}^k(x)f[n-k+j],
```
where ``s^k_x[j] ≡ s_{k-j}(x)^k``.
#### Example:
```
k = 3
x = 0
ldw = create_lagrange_differentiation_weights(k,x); println(ldw)
```
"""
function create_lagrange_differentiation_weights(k::Int, x::T) where T<:Real

    ∇ = CamiXon.f_diff_weights_array(k)
    coeffs = CamiXon.f_diff_expansion_coeffs_differentiation(k,-k+x)

    return CamiXon.f_diff_expansion_weights(coeffs,∇)

end

# =============================== create_lagrange_differentiation_matrix(k) ====

@doc raw"""
    create_lagrange_differentiation_matrix(k::Int)

Lagrange differentiation matrix, ``m[i,j]=s_{k-j}^k(i)``, for ``k^{th}``-order lagrangian differentiation,
```math
\frac{dy}{dx}[i]= \sum_{j=0}^{k}m[i,j]y[j],
```
#### Example:
```
k = 3
create_lagrange_differentiation_matrix(k)
 4×4 Matrix{Rational{Int64}}:
  -11//6   3//1  -3//2   1//3
   -1//3  -1//2   1//1  -1//6
    1//6  -1//1   1//2   1//3
   -1//3   3//2  -3//1  11//6
```
"""
function create_lagrange_differentiation_matrix(k::Int)

    m = Matrix{Rational{Int}}(undef,k+1,k+1)

    ∇ = CamiXon.f_diff_weights_array(k)

    for i=0:k
        coeffs = CamiXon.f_diff_expansion_coeffs_differentiation(k,-k+i)
        m[1+i,1:k+1] = CamiXon.f_diff_expansion_weights(coeffs,∇)
    end

    return m

end

# ================================f_diff_expansion_coeffs_array_differentiation(k, m) ====

@doc raw"""
    lagrange_differentiation(f::Vector{Float64}, domain::ClosedInterval{Float64}; k=1, m=1)

``k^{th}``-order lagrangian *differentiation* of the analytic function ``f``, tabulated
in forward order on a regular grid of ``n`` points, ``f[1],\ \ldots,\ f[n]``;
``m`` is the multiplier for intermediate positions
#### Example:
```
f = [0.0,1,2,3,4,5]
domain = 0.0..5.0
X,Y = lagrangian_differentiation(f, domain; k=2, i = 0); println(X,Y)
 (0.0:1.0:5.0, [1.0, 1.0, 1.0, 1.0, 1.0, 1.0])
```
"""
function lagrange_differentiation(f::Vector{Float64}, domain::ClosedInterval{Float64}; k=1, m=1)
# ======================================================================================
#   lagrangian (k+1)-point differentiation at i interpolation points
# ======================================================================================
    n = Base.length(f)

    ∇ = CamiXon.f_diff_weights_array(k)
    l = [CamiXon.f_diff_expansion_coeffs_differentiation(k, x) for x=-k:1/m:0]
    w = [CamiXon.f_diff_expansion_weights(l[i], ∇) for i ∈ Base.eachindex(l)]
    w1 = Base.append!(repeat(w[1:m],n-k-1),w)
    w2 = CamiXon.f_diff_function_sequences(f, k, m)

    X = Base.range(domain.left, domain.right, length=(n-1)*m+1)
    Y = (n-1)/(domain.right-domain.left) .*  [w1[i] ⋅ w2[i] for i ∈ Base.eachindex(w1)]

    return X, Y

end

# ======================== trapezoidal_weights(k; rationalize=false, devisor=false) ===================================

@doc raw"""
    trapezoidal_weights(k::Int [; rationalize=false [, devisor=false]])

Weight coefficient vector ``a=[a_1,\cdots,\ a_k]`` of trapeziodal rule optimized for functions of polynomial form,
```math
    ∫_0^n f(x) dx = a_1 (f_0+f_n)+\cdots+a_k (f_{k-1}+f_{n-k+1}) + (f_k+\cdots+f_{n-k}),
```
where ``k`` is *odd*. The rule is exact for polynonials of degree ``d=0,\ 1,\cdots,\ k-1``.
For ``k=1`` the rule reduces to the ordinary trapezoidal rule. By default the output is in Float64,
optionally the output is rational, with or without specification of the gcd devisor.
#### Example::
```
[trapezoidal_weights(k; rationalize=true, devisor=true) for k=1:2:9]
5-element Vector{Tuple{Int64, Int64, Vector{Int64}}}:
 (1, 2, [1])
 (3, 24, [9, 28, 23])
 (5, 1440, [475, 1902, 1104, 1586, 1413])
 (7, 120960, [36799, 176648, 54851, 177984, 89437, 130936, 119585])
 (9, 7257600, [2082753, 11532470, 261166, 16263486, -1020160, 12489922, 5095890, 7783754, 7200319])
```
"""
function trapezoidal_weights(k::Int; rationalize=false, devisor=false)
# ==============================================================================
# trapezoidal_weights(k; rationalize=false)
# ==============================================================================
    Base.isodd(k) || return error("jwError: method requires odd k")

    l = k - 1
    σ = Base.Matrix{Int}(undef,k,l+1)

    for i=0:k-1
        σ[i+1,:] = CamiXon.polynom_power([i,-1],l)      # corresponds to coeffs a0,...,ak
        σ[i+1,1] = σ[i+1,1] + i^l                       # correction of coeff a0
    end

    F = CamiXon.faulhaber_polynom(l+1)
    s = [CamiXon.polynom_power([-k,1],p) for p=0:l+1] .* F
    s[1][1] = -CamiXon.faulhaber_summation(k-1,l)//1

    c = [Base.sum([s[p+1][i+1] for p=i:l+1]) for i=0:l+1][1:end-1]

    σ = Base.inv(Base.transpose(σ))

    o = -σ * c  # solving matrix equation results in trapezoidal_weights as real numbers

    if rationalize
        a = CamiXon.f_diff_expansion_coeffs_adams_moulton(k)
        D = Base.denominator(Base.gcd(a))       # trapezoidal_weights divisor == Adams-Moulton devisor
        o = devisor ? (k, D, Base.round.(Int, o* D)) : Base.round.(Int, o* D) // D              # convert to Rational{Int}
    end

    return o

end

# ======================== trapezoidal_integration(f, domain; k=5) =============

@doc raw"""
    trapezoidal_integration(f, domain, weights)

Integral of the tabulated function ``f=[f_0,\cdots,\ f_n]`` over the `domain` ``a..b`` using the optimized
trapezoidal rule with endpoint correction by the weightsvector `weights`,
```math
    ∫_0^n f(x) dx = a_1 (f_0+f_n)+\cdots+a_k (f_{k-1}+f_{n-k+1}) + (f_k+\cdots+f_{n-k}).
```
The rule is exact for polynonials of degree ``d=0,\ 1,\cdots,\ k-1``. For ``k=1`` the rule reduces to the
ordinary trapezoidal rule (weights = [1/2]).
#### Examples::
```
p = 3
c = [1 for i=0:p]
pol = ImmutablePolynomial(c,:z)
Ipol = integrate(pol)
n = 10

domain = 0.0..5.0
x = collect(range(domain, n))
f = pol.(x .-2.5)

w3 = trapezoidal_weights(3)
trapezoidal_integration(f, domain, w3)
 15.416666666666673

Ipol(2.5)-Ipol(-2.5)
 15.41666666666666
```
"""
function trapezoidal_integration(f, domain, weights)

    n = Base.length(f)
    k = Base.length(weights)
    s = (domain.right-domain.left)/(n-1)
    a = Base.ones(n); a[1:k] = weights; a[end-k+1:end] = Base.reverse(weights)
    o = (f ⋅ a) * s

    return o

end


text:
#=
  Monty Hall problem

  From the PyMC3 model in the talk
  "Carlo Revolution is Open Source: Probabilistic Programming with PyMC3?Austin Rochford"
  https://www.safaribooksonline.com/videos/open-data-science/9780135432792/9780135432792-ODSC_11
  Around time 9:00

  Which mean that if we selected d1 it will be 1/3 chance of being the price door.
  Changing to d3 would - however - give a 2/3 change of getting the price.

  Distributions of variable prize
  door3      =>    6685  (0.668500)
  door1      =>    3312  (0.331200)
  door2      =>       3  (0.000300)

  See ~/blog/monty_hall.blog
      ~/psi/monty_hall.psi
      ~/webppl/monty_hall.wppl

=#

# using Memoization
using Turing # , StatsPlots, DataFrames
# using ReverseDiff, Zygote, Tracker
# Turing.setadbackend(:reversediff)
# Turing.setadbackend(:zygote)
# Turing.setadbackend(:tracker)
include("jl_utils.jl")

@model function monty_hall()
    # There are 3 doors.
    d1 = 1
    d2 = 2
    d3 = 3

    # The prize can be behind any door 1..3.
    prize = DiscreteUniform(1,3)

    # Which door will Monty open?
    # Assumption (WLOG): We always select door 1.
    prize ~ DiscreteUniform(1,3)
    open  ~ DiscreteUniform(1,3);
    open ~ (prize == d1) ? Categorical([0.0, 0.5,0.5]) :
            (prize == d2 ? Categorical([0.0, 0.0,1.0]) : Categorical([0.0, 1.0,0.0]))

    # We see that Monty opens door 2.
    true ~ Dirac(open == d2)

    # What are the probabilities that the price is behind
    # - door d1 (the one we selected, i.e don't switch)
    # - or door d3 (i.e. switch door)

end

model = monty_hall()
num_chns = 4

# chns = sample(model, Prior(), 10_000)
chns = sample(model, MH(), 10_000)
# chns = sample(model, PG(20), 1_000)
# chns = sample(model, IS(), 10_000)
# chns = sample(model, SMC(), 10_000)

display(chns)

show_var_dist_pct(chns, :prize,["door1","door2","door3"])

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
# Type hierarchy to aid in splitting up of SHA2 algorithms
# as SHA224/256 are similar, and SHA-384/512 are similar
@compat abstract type SHA_CTX end
@compat abstract type SHA2_CTX <: SHA_CTX end
@compat abstract type SHA3_CTX <: SHA_CTX end
import Base: copy

# We derive SHA1_CTX straight from SHA_CTX since it doesn't have a
# family of types like SHA2 or SHA3 do
type SHA1_CTX <: SHA_CTX
    state::Array{UInt32,1}
    bytecount::UInt64
    buffer::Array{UInt8,1}
    W::Array{UInt32,1}
end

# SHA2 224/256/384/512-bit Context Structures
type SHA2_224_CTX <: SHA2_CTX
    state::Array{UInt32,1}
    bytecount::UInt64
    buffer::Array{UInt8,1}
end

type SHA2_256_CTX <: SHA2_CTX
    state::Array{UInt32,1}
    bytecount::UInt64
    buffer::Array{UInt8,1}
end

type SHA2_384_CTX <: SHA2_CTX
    state::Array{UInt64,1}
    bytecount::UInt128
    buffer::Array{UInt8,1}
end

type SHA2_512_CTX <: SHA2_CTX
    state::Array{UInt64,1}
    bytecount::UInt128
    buffer::Array{UInt8,1}
end

# Typealias common nicknames for SHA2 family of functions
const SHA224_CTX = SHA2_224_CTX
const SHA256_CTX = SHA2_256_CTX
const SHA384_CTX = SHA2_384_CTX
const SHA512_CTX = SHA2_512_CTX

# SHA3 224/256/384/512-bit context structures
type SHA3_224_CTX <: SHA3_CTX
    state::Array{UInt64,1}
    bytecount::UInt128
    buffer::Array{UInt8,1}
end
type SHA3_256_CTX <: SHA3_CTX
    state::Array{UInt64,1}
    bytecount::UInt128
    buffer::Array{UInt8,1}
end
type SHA3_384_CTX <: SHA3_CTX
    state::Array{UInt64,1}
    bytecount::UInt128
    buffer::Array{UInt8,1}
end
type SHA3_512_CTX <: SHA3_CTX
    state::Array{UInt64,1}
    bytecount::UInt128
    buffer::Array{UInt8,1}
end

# Define constants via functions so as not to bloat context objects.  Yay dispatch!

# Digest lengths for SHA1, SHA2 and SHA3.  This is easy to figure out from the typename
digestlen(::Type{SHA1_CTX}) = 20
digestlen(::Type{SHA2_224_CTX}) = 28
digestlen(::Type{SHA3_224_CTX}) = 28
digestlen(::Type{SHA2_256_CTX}) = 32
digestlen(::Type{SHA3_256_CTX}) = 32
digestlen(::Type{SHA2_384_CTX}) = 48
digestlen(::Type{SHA3_384_CTX}) = 48
digestlen(::Type{SHA2_512_CTX}) = 64
digestlen(::Type{SHA3_512_CTX}) = 64

# SHA1 and SHA2 have differing element types for the internal state objects
state_type(::Type{SHA1_CTX}) = UInt32
state_type(::Type{SHA2_224_CTX}) = UInt32
state_type(::Type{SHA2_256_CTX}) = UInt32
state_type(::Type{SHA2_384_CTX}) = UInt64
state_type(::Type{SHA2_512_CTX}) = UInt64

# blocklen is the number of bytes of data processed by the transform!() function at once
blocklen(::Type{SHA1_CTX}) = UInt64(64)
blocklen(::Type{SHA2_224_CTX}) = UInt64(64)
blocklen(::Type{SHA2_256_CTX}) = UInt64(64)
blocklen(::Type{SHA2_384_CTX}) = UInt64(128)
blocklen(::Type{SHA2_512_CTX}) = UInt64(128)

blocklen(::Type{SHA3_224_CTX}) = UInt64(25*8 - 2*digestlen(SHA3_224_CTX))
blocklen(::Type{SHA3_256_CTX}) = UInt64(25*8 - 2*digestlen(SHA3_256_CTX))
blocklen(::Type{SHA3_384_CTX}) = UInt64(25*8 - 2*digestlen(SHA3_384_CTX))
blocklen(::Type{SHA3_512_CTX}) = UInt64(25*8 - 2*digestlen(SHA3_512_CTX))

# short_blocklen is the size of a block minus the width of bytecount
short_blocklen{T<:Union{SHA1_CTX,SHA2_CTX}}(::Type{T}) = blocklen(T) - 2*sizeof(state_type(T))

# Once the "blocklen" methods are defined, we can define our outer constructors for SHA types:
SHA2_224_CTX() = SHA2_224_CTX(copy(SHA2_224_initial_hash_value), 0, zeros(UInt8, blocklen(SHA2_224_CTX)))
SHA2_256_CTX() = SHA2_256_CTX(copy(SHA2_256_initial_hash_value), 0, zeros(UInt8, blocklen(SHA2_256_CTX)))
SHA2_384_CTX() = SHA2_384_CTX(copy(SHA2_384_initial_hash_value), 0, zeros(UInt8, blocklen(SHA2_384_CTX)))
SHA2_512_CTX() = SHA2_512_CTX(copy(SHA2_512_initial_hash_value), 0, zeros(UInt8, blocklen(SHA2_512_CTX)))

SHA3_224_CTX() = SHA3_224_CTX(zeros(UInt64, 25), 0, zeros(UInt8, blocklen(SHA3_224_CTX)))
SHA3_256_CTX() = SHA3_256_CTX(zeros(UInt64, 25), 0, zeros(UInt8, blocklen(SHA3_256_CTX)))
SHA3_384_CTX() = SHA3_384_CTX(zeros(UInt64, 25), 0, zeros(UInt8, blocklen(SHA3_384_CTX)))
SHA3_512_CTX() = SHA3_512_CTX(zeros(UInt64, 25), 0, zeros(UInt8, blocklen(SHA3_512_CTX)))

# Nickname'd outer constructor methods for SHA2
SHA224_CTX = SHA2_224_CTX
SHA256_CTX = SHA2_256_CTX
SHA384_CTX = SHA2_384_CTX
SHA512_CTX = SHA2_512_CTX

# SHA1 is special; he needs extra workspace
SHA1_CTX() = SHA1_CTX(copy(SHA1_initial_hash_value), 0, zeros(UInt8, blocklen(SHA1_CTX)), Vector{UInt32}(80))

# Copy functions
copy{T <: SHA1_CTX}(ctx::T) = T(copy(ctx.state), ctx.bytecount, copy(ctx.buffer), copy(ctx.W))
copy{T <: SHA2_CTX}(ctx::T) = T(copy(ctx.state), ctx.bytecount, copy(ctx.buffer))
copy{T <: SHA3_CTX}(ctx::T) = T(copy(ctx.state), ctx.bytecount, copy(ctx.buffer))

# Make printing these types a little friendlier
import Base.show
show(io::IO, ::SHA1_CTX) = write(io, "SHA1 hash state")
show(io::IO, ::SHA2_224_CTX) = write(io, "SHA2 224-bit hash state")
show(io::IO, ::SHA2_256_CTX) = write(io, "SHA2 256-bit hash state")
show(io::IO, ::SHA2_384_CTX) = write(io, "SHA2 384-bit hash state")
show(io::IO, ::SHA2_512_CTX) = write(io, "SHA2 512-bit hash state")
show(io::IO, ::SHA3_224_CTX) = write(io, "SHA3 224-bit hash state")
show(io::IO, ::SHA3_256_CTX) = write(io, "SHA3 256-bit hash state")
show(io::IO, ::SHA3_384_CTX) = write(io, "SHA3 384-bit hash state")
show(io::IO, ::SHA3_512_CTX) = write(io, "SHA3 512-bit hash state")


text:
# reshape test
x = create_var(6, 4)
y = reshape(x, 8, 3)
z = reshape(x, 4, 6)
v = vec(x)
test_expr(y, (x,))
test_expr(z, (x,))
test_expr(v, (x,))

info("All reshape tests passed")


text:
# Rather than use a condirional tree to determine which peoperty is
# being queried by `getproperty`, we can method specialize on Val types.
# We should then be able to automate `propertynames` based on the
# defined methods.

# A struct might have fields, which are exposed as properties.  The
# `Val` methods will not shadow the method that implements that.

export @njl_getprop

"""
    propertynames_from_val_methods(t::Type, private::Bool)
Compute and return the `propertynames` for type `t` from
its fields and `Val` specialized methods.
"""
function propertynames_from_val_methods(t::Type, private::Bool)
    props = collect(fieldnames(t))
    for m in methods(Base.getproperty, [t, Any]).ms
        if m.sig isa DataType
            if m.sig.parameters[2] == t
                p = m.sig.parameters[3]
                if p <: Val
                    push!(props, p.parameters[1])
                end
            end
        end
    end
    props
end

"""
    @njl_getprop MyStruct
Define the methods necessary so that `Val `getproperties` of`MyStruct`
will find `Val` specialized properties.
"""
macro njl_getprop(structname)
    quote
        Base.getproperty(o::$(esc(structname)), prop::Symbol) =
            Base.getproperty(o, Val(prop))
        Base.getproperty(o::$(esc(structname)), prop::Val{T}) where {T} =
            getfield(o, T)
        Base.propertynames(o::$(esc(structname)), private::Bool=false) =
            propertynames_from_val_methods(typeof(o), private)
    end
end


text:
struct Context
    threadid::Int # So that the thread that the context is associated with calls the opengl operations in the queue
    queue::Vector{Function} # Closures to be executed when the context becomes current again
end

Context() = Context(threadid(), Function[])

Base.push!(c::Context, f::Function) = push!(c.queue, f)

const GLOBAL_CONTEXT = Base.RefValue{Union{Nothing,Context}}(nothing)
const GLOBAL_CONTEXTS = Dict{Any, Context}()

current_context() = GLOBAL_CONTEXT[]
is_current_context(x::Context) = x == GLOBAL_CONTEXT[]
is_current_context(x) =
    haskey(GLOBAL_CONTEXTS, x) && GLOBAL_CONTEXTS[x] === GLOBAL_CONTEXT[]

clear_context!() = GLOBAL_CONTEXT[] = Context()

function set_context!(x)
    if haskey(GLOBAL_CONTEXTS, x)
        c = GLOBAL_CONTEXTS[x]
        GLOBAL_CONTEXT[] = c
        @tspawnat c.threadid begin
            for f in c.queue
                f()
            end
            empty!(c.queue)
        end
    else
        c = Context()
        GLOBAL_CONTEXTS[x] = c
        GLOBAL_CONTEXT[] = c
    end
end

function exists_context()
    if current_context() === nothing
        @error "Couldn't find valid OpenGL Context. OpenGL Context active?"
    end
end

function context_command(c::Context, f::Function)
    if !is_current_context(c)
        push!(c, f)
    else
        f()
    end
end


text:
#-----------------------------------------------------------------------------------------------
# 1st run - packages installation

# ERROR: LoadError: ArgumentError: Package DataFrames not found in current path:
# - Run `import Pkg; Pkg.add("DataFrames")` to install the DataFrames package.
import Pkg; Pkg.add("DataFrames")
import Pkg;
Pkg.add("CSV")
Pkg.add("StatsBase")

#-----------------------------------------------------------------------------------------------
using DataFrames, CSV, Statistics, StatsBase
data = CSV.read("../externals/Core.Math.Data/data/Pejcic_318.csv", copycols = true)

# @show data.ATV
# println()

for col in eachcol(data);
    if  (
        eltype(col[2]) <: Real)
        println(mean(col[2])
        );
    end;
end;


text:
let
    #= line 1 =#
    var"|Prim|b_op_Getter"::Function = function op_Getter(x)
                #= none:1 =#
                #= none:4 =#
                function (a,)
                    #= none:4 =#
                    #= none:7 =#
                    x(a)
                end
            end::Function
    var"|Prim|c_unit"::Symbol = Symbol("{}")::Symbol
    var"|Prim|d_Vec"::Symbol = Symbol("Vec@10_|Prim|e")::Symbol
    var"|Prim|f_make"::Function = function (var"|Prim|g_t",)
                #= ../examples/prim.mlfs:14 =#
                var"|Prim|g_t"::Function
            end::Function
    var"|Prim|h_mkVec"::Function = (_->begin
                    #= none:1 =#
                    #= none:3 =#
                    []
                end)::Function
    var"|Prim|i_pushVec"::Any = (x->begin
                    #= none:1 =#
                    #= none:3 =#
                    function (a,)
                        #= none:3 =#
                        #= none:6 =#
                        push!(x, a)
                        #= none:8 =#
                        ()
                    end
                end)::Any
    var"|Prim|j_println"::Function = function (x,)
                #= none:1 =#
                #= none:4 =#
                println(x)
                #= none:6 =#
                ()
            end::Function
    var"|Prim|k"::Function = function (_,)
            #= ../examples/prim.mlfsf:1 =#
            var"|Prim|j_println"
        end
    var"|Prim|l"::Function = function (_,)
            #= ../examples/prim.mlfsf:1 =#
            var"|Prim|i_pushVec"
        end
    var"|Prim|m"::Function = function (_,)
            #= ../examples/prim.mlfsf:1 =#
            var"|Prim|h_mkVec"
        end
    var"|Prim|n"::Function = function (_,)
            #= ../examples/prim.mlfsf:1 =#
            var"|Prim|f_make"
        end
    var"|Prim|o"::Function = function (_,)
            #= ../examples/prim.mlfsf:1 =#
            Symbol("Vec@10_|Prim|e")
        end
    var"|Prim|p"::Function = function (_,)
            #= ../examples/prim.mlfsf:1 =#
            Symbol("{}")
        end
    var"|Prim|q"::Function = function (_,)
            #= ../examples/prim.mlfsf:1 =#
            var"|Prim|b_op_Getter"
        end
    var"|Prim|"::Any = Symbol("Module Prim")
    Prim::Any = Symbol("namespace {Module Prim, {{make, |Prim|f_make, forall t@13071667606244709930 a@17423321405210106963.implicit[(Type t@13071667606244709930)->t@13071667606244709930 a@17423321405210106963]->(Type t@13071667606244709930)->t@13071667606244709930 a@17423321405210106963}, {unit, |Prim|c_unit, Type {}}, {mkVec, |Prim|h_mkVec, implicit[forall a@10415488706922051249.(Type Vec@10_|Prim|e)->Vec@10_|Prim|e a@10415488706922051249]}, {println, |Prim|j_println, forall a@14273736158010134791.a@14273736158010134791->{}}, {op_Getter, |Prim|b_op_Getter, forall a@3960223772273614135 b@8222540082496598410 c@13760051442097314445.(field {a@3960223772273614135, b@8222540082496598410, c@13760051442097314445})->a@3960223772273614135->c@13760051442097314445}, {Vec, |Prim|d_Vec, Type Vec@10_|Prim|e}, {pushVec, |Prim|i_pushVec, implicit[forall a@13063604110054457948.field {Vec@10_|Prim|e a@13063604110054457948, push, a@13063604110054457948->{}}]}}}")
    (var"|Prim|j_println"::Function)(1::Int64)::Tuple
    (var"|Prim|j_println"::Function)(:i64::Symbol)::Tuple
    var"|main|r_a"::Any = (var"|Prim|f_make"(var"|Prim|h_mkVec"::Function)::Function)(Symbol("Vec@10_|Prim|e")::Symbol)::Any
    var"|main|s_c"::Tuple = ((op_Getter(var"|Prim|i_pushVec"::Any))(var"|main|r_a"::Any)::Function)(1::Int64)::Tuple
    var"|main|t_mk2"::Function = function (var"|main|u_x",)
                #= ../examples/main.mlfs:17 =#
                (var"|main|u_x"::Any, var"|main|u_x"::Any)::Tuple
            end::Function
    (var"|Prim|j_println"::Function)((var"|main|t_mk2"::Function)(var"|main|r_a"::Any)::Tuple)::Tuple
    var"|main|v_index"::Function = (x::Vector->begin
                    #= none:1 =#
                    #= none:3 =#
                    i::Int64->begin
                            #= none:3 =#
                            #= none:5 =#
                            x[i]
                        end
                end)::Function
    var"|main|w_CanAdd1"::Symbol = Symbol("CanAdd1@24_|main|x")::Symbol
    var"|main|y_int_add1"::Function = (_->begin
                    #= none:1 =#
                    #= none:3 =#
                    x::Int64->begin
                            #= none:3 =#
                            #= none:5 =#
                            x + 1
                        end
                end)::Function
    var"|main|z_float_add1"::Function = (_->begin
                    #= none:1 =#
                    #= none:3 =#
                    x::Float64->begin
                            #= none:3 =#
                            #= none:5 =#
                            x + 1.0
                        end
                end)::Function
    var"|main|ab_vec_head_add1"::Function = function (var"|main|bb_instance",)
                #= ../examples/main.mlfs:35 =#
                function (var"|main|cb_x",)
                        #= ../examples/main.mlfs:35 =#
                        ((var"|main|bb_instance"::Function)(Symbol("CanAdd1@24_|main|x")::Symbol)::Function)(((var"|main|v_index"::Function)(var"|main|cb_x"::Any)::Function)(1::Int64)::Any)::Any
                    end::Function
            end::Function
    (var"|Prim|j_println"::Function)((var"|main|ab_vec_head_add1"(var"|main|y_int_add1"::Function)::Function)(var"|main|r_a"::Any)::Int64)::Tuple
    var"|main|db"::Function = function (_,)
            #= ../examples/main.mlfsf:1 =#
            var"|main|ab_vec_head_add1"
        end
    var"|main|eb"::Function = function (_,)
            #= ../examples/main.mlfsf:1 =#
            var"|main|z_float_add1"
        end
    var"|main|fb"::Function = function (_,)
            #= ../examples/main.mlfsf:1 =#
            var"|main|y_int_add1"
        end
    var"|main|gb"::Function = function (_,)
            #= ../examples/main.mlfsf:1 =#
            Symbol("CanAdd1@24_|main|x")
        end
    var"|main|hb"::Function = function (_,)
            #= ../examples/main.mlfsf:1 =#
            var"|main|v_index"
        end
    var"|main|ib"::Function = function (_,)
            #= ../examples/main.mlfsf:1 =#
            var"|main|t_mk2"
        end
    var"|main|jb"::Function = function (_,)
            #= ../examples/main.mlfsf:1 =#
            var"|main|s_c"
        end
    var"|main|kb"::Function = function (_,)
            #= ../examples/main.mlfsf:1 =#
            var"|main|r_a"
        end
    var"|main|lb"::Function = function (_,)
            #= ../examples/main.mlfsf:1 =#
            var"|Prim|i_pushVec"
        end
    var"|main|mb"::Function = function (_,)
            #= ../examples/main.mlfsf:1 =#
            Symbol("Vec@10_|Prim|e")
        end
    var"|main|nb"::Function = function (_,)
            #= ../examples/main.mlfsf:1 =#
            var"|Prim|b_op_Getter"
        end
    var"|main|ob"::Function = function (_,)
            #= ../examples/main.mlfsf:1 =#
            var"|Prim|j_println"
        end
    var"|main|pb"::Function = function (_,)
            #= ../examples/main.mlfsf:1 =#
            var"|Prim|h_mkVec"
        end
    var"|main|qb"::Function = function (_,)
            #= ../examples/main.mlfsf:1 =#
            Symbol("{}")
        end
    var"|main|rb"::Function = function (_,)
            #= ../examples/main.mlfsf:1 =#
            var"|Prim|f_make"
        end
    var"|main|sb"::Function = function (_,)
            #= ../examples/main.mlfsf:1 =#
            var"|Prim|"
        end
    var"|main|"::Any = Symbol("Module main")
    main::Any = Symbol("namespace {Module main, {{Prim, |Prim|, Module Prim}, {make, |Prim|f_make, forall t@13071667606244709930 a@17423321405210106963.implicit[(Type t@13071667606244709930)->t@13071667606244709930 a@17423321405210106963]->(Type t@13071667606244709930)->t@13071667606244709930 a@17423321405210106963}, {c, |main|s_c, {}}, {mkVec, |Prim|h_mkVec, implicit[forall a@10415488706922051249.(Type Vec@10_|Prim|e)->Vec@10_|Prim|e a@10415488706922051249]}, {int_add1, |main|y_int_add1, implicit[(Type CanAdd1@24_|main|x)->i64->i64]}, {op_Getter, |Prim|b_op_Getter, forall a@3960223772273614135 b@8222540082496598410 c@13760051442097314445.(field {a@3960223772273614135, b@8222540082496598410, c@13760051442097314445})->a@3960223772273614135->c@13760051442097314445}, {float_add1, |main|z_float_add1, implicit[(Type CanAdd1@24_|main|x)->f64->f64]}, {CanAdd1, |main|w_CanAdd1, Type CanAdd1@24_|main|x}, {index, |main|v_index, forall a@12721974860504392992.(Vec@10_|Prim|e a@12721974860504392992)->i64->a@12721974860504392992}, {unit, |Prim|c_unit, Type {}}, {println, |Prim|j_println, forall a@14273736158010134791.a@14273736158010134791->{}}, {mk2, |main|t_mk2, forall a@11153002359981584258.a@11153002359981584258->{a@11153002359981584258, a@11153002359981584258}}, {a, |main|r_a, Vec@10_|Prim|e i64}, {vec_head_add1, |main|ab_vec_head_add1, forall a@3814065663426908098.implicit[(Type CanAdd1@24_|main|x)->a@3814065663426908098->a@3814065663426908098]->(Vec@10_|Prim|e a@3814065663426908098)->a@3814065663426908098}, {Vec, |Prim|d_Vec, Type Vec@10_|Prim|e}, {pushVec, |Prim|i_pushVec, implicit[forall a@13063604110054457948.field {Vec@10_|Prim|e a@13063604110054457948, push, a@13063604110054457948->{}}]}}}")
    0
end


text:
using MLJ

LinearRegressor = @load LinearRegressor pkg=MLJLinearModels

import RDatasets: dataset
import DataFrames: describe, select, Not, rename!
boston = dataset("MASS", "Boston")
first(boston, 3)

describe(boston, :mean, :std, :eltype)

data = coerce(boston, autotype(boston, :discrete_to_continuous));

y = data.MedV
X = select(data, Not(:MedV));

mdl = LinearRegressor()

X_uni = select(X, :LStat) # only a single feature
mach_uni = machine(mdl, X_uni, y)
fit!(mach_uni)

fp = fitted_params(mach_uni)
@show fp.coefs
@show fp.intercept

using PyPlot

figure(figsize=(8,6))
plot(X.LStat, y, ls="none", marker="o")
Xnew = (LStat = collect(range(extrema(X.LStat)..., length=100)),)
plot(Xnew.LStat, MLJ.predict(mach_uni, Xnew))

mach = machine(mdl, X, y)
fit!(mach)

fp = fitted_params(mach)
coefs = fp.coefs
intercept = fp.intercept
for (name, val) in coefs
    println("$(rpad(name, 8)):  $(round(val, sigdigits=3))")
end
println("Intercept: $(round(intercept, sigdigits=3))")

ŷ = MLJ.predict(mach, X)
round(rms(ŷ, y), sigdigits=4)

figure(figsize=(8,6))
res = ŷ .- y
stem(res)

figure(figsize=(8,6))
hist(res, density=true)

X2 = hcat(X, X.LStat .* X.Age);

rename!(X2, :x1 => :interaction);

mach = machine(mdl, X2, y)
fit!(mach)
ŷ = MLJ.predict(mach, X2)
round(rms(ŷ, y), sigdigits=4)

X3 = hcat(X.LStat, X.LStat.^2) |> MLJ.table
mach = machine(mdl, X3, y)
fit!(mach)
ŷ = MLJ.predict(mach, X3)
round(rms(ŷ, y), sigdigits=4)

Xnew = (LStat = Xnew.LStat, LStat2 = Xnew.LStat.^2)

figure(figsize=(8,6))
plot(X.LStat, y, ls="none", marker="o")
plot(Xnew.LStat, MLJ.predict(mach, Xnew))

# This file was generated using Literate.jl, https://github.com/fredrikekre/Literate.jl

max_issues_repo_name: alan-turing-institute/MLJTutorials

max_forks_repo_name: alan-turing-institute/MLJTutorials


text:
"""
    $(TYPEDEF)
Superclass of all key-value datastores.
"""
abstract type AbstractDataStore{T}
end


text:
export uncsolver

"""
    uncsolver(nlp; options...)

Template for an unconstrained solver using the first order conditions

    ‖∇f(xₖ)‖ ≤ ϵₐ + ϵᵣ‖∇f(x₀)‖,

where ϵₐ is an absolute tolerance and ϵᵣ is a relative tolerance.

Options:
- atol: absolute tolerance for the first order condition (default: 1e-6)
- rtol: relative tolerance for the first order condition (default: 1e-6)
- max_eval: maximum number of [functions] evaluations, use ≤ 0 for unlimited (default: 1000)
- max_iter: maximum number of iterations, use ≤ 0 for unlimited (default: 0)
- max_time: maximum elapsed time in seconds, use ≤ 0 for unlimited (default: 10)

Disclaimers for the developer:
- nlp should be the only mandatory argument
- these five options are the current default for other JSO-compliant solvers
- always return a GenericExecutionStats
"""
function uncsolver(
  nlp::AbstractNLPModel;
  atol::Real = 1e-6,
  rtol::Real = 1e-6,
  max_eval::Int = 1000,
  max_iter::Int = 0,
  max_time::Float64 = 10.0
)

  if !unconstrained(nlp)
    error("Problem is not unconstrained")
  end

  x = copy(nlp.meta.x0)

  f(x) = obj(nlp, x)
  ∇f(x) = grad(nlp, x)
  #= Re: Hessian
  Hessians are expensive, so `hess(nlp, x)` only returns the lower triangle of the symmetric Hessian at x.
  You can use it as a symmetric matrix with

      H(x) = Symmetric(hess(nlp, x), :L)

  However, most symmetric linear solvers only use the lower half.
  Furthermore, for large problems with sparse Hessians, you can obtain the triplet list `(i, j, aᵢⱼ)` using

      rows, cols = hess_structure(nlp)
      vals = hess_coord(nlp, x)
      # or
      hess_coord!(nlp, x, vals)

  Finally, if your method only uses matrix-vector products, you can use

      Hv = hprod(nlp, x, v)
      # or
      hprod!(nlp, x, v, Hv)

  Or potentially, a linear operator that automatically does that (read more on github.com/JuliaSmoothOptimizers/LinearOperators.jl).
  In the example below, the Hessian is not created explicitly (unless the model itself doesn't support it):

      H = hess_op(nlp, x)
      Hv = H * v

  Notice that you can't access elements of the Hessian in this case.
  =#

  fx = f(x)
  ∇fx = ∇f(x)

  ϵ = atol + rtol * norm(∇fx)
  t₀ = time()

  iter = 0
  Δt = time() - t₀
  solved = norm(∇fx) < ϵ # First order stationary
  tired = neval_obj(nlp) ≥ max_eval > 0|| iter ≥ max_iter > 0 || Δt ≥ max_time > 0 # Excess time, iteration, evaluations

  # status must be one of a few options found in SolverTools.show_statuses()
  # A good default value is :unknown.
  status = :unknown

  # log_header is up for some rewrite in the future. For now, it simply prints the column names with some spacing
  @info log_header(
    [:iter, :fx, :ngx, :nf, :Δt],
    [Int, Float64, Float64, Int, Float64],
    hdr_override=Dict(:fx => "f(x)", :ngx => "‖∇f(x)‖", :nf => "#f")
  )
  # log_row uses the type information of each value, thus we use `Any` here.
  @info log_row(
    Any[iter, fx, norm(∇fx), neval_obj(nlp), Δt]
  )

  # This template implements a simple steepest descent method without any hopes of working.
  # This is where most of your change will happen
  α = 1.0
  while !(solved || tired)
    x -= α * ∇fx
    α *= 0.99
    fx = f(x)
    ∇fx = ∇f(x)

    iter += 1
    Δt = time() - t₀
    solved = norm(∇fx) < ϵ # First order stationary
    tired = neval_obj(nlp) ≥ max_eval > 0|| iter ≥ max_iter > 0 || Δt ≥ max_time > 0 # Excess time, iteration, evaluations

    @info log_row(
      Any[iter, fx, norm(∇fx), neval_obj(nlp), Δt]
    )
  end

  if solved
    status = :first_order
  elseif tired
    if neval_obj(nlp) ≥ max_eval > 0
      status = :max_eval
    elseif iter ≥ max_iter > 0
      status = :max_iter
    elseif Δt ≥ max_time > 0
      status = :max_time
    end
  end

  return GenericExecutionStats(
    status,
    nlp,
    solution=x,
    objective=f(x),
    dual_feas=norm(∇fx),
    elapsed_time=Δt,
    iter=iter
  )
end


text:
function find_vtd_directory(error_on_fail::Bool=true)
    @static is_linux() ? let
                hd = homedir()
                for content in readdir(hd)
                    path = joinpath(hd, content)
                    if ismatch(r"VTD\.\d+\.\d+", content) && isdir(path)
                        return path
                    end
                end
                error_on_fail && error("VTD directory not found")
            end : let
                error_on_fail && error("non-linux os's not supported!")
            end

    return "" # return empty string if no error was thrown and path was not found
end
function start_vires_vtd_tasks()
    dir = joinpath(find_vtd_directory(), "Data", "Setups", "Current", "Bin")
    cd(dir) do
        run(`./startTasks`)
    end
end
function stop_vires_components()
    dir = joinpath(find_vtd_directory(), "Data", "Setups", "Current", "Bin")
    cd(dir) do
        run(`./stopTasks`)
    end
end
function rdbsniffer(; port::Integer=RDB_FEEDBACK_PORT)
    # run RDB sniffer on given port in debug mode

    @lintpragma("Ignore unused port")
    dir = joinpath(find_vtd_directory(), "Runtime", "Tools", "RDBSniffer")
    cd(dir) do
        run(`./rdbSniffer -c udp -p $port -d`)
    end
end

max_stars_repo_name: sisl/VirtualTestDrive

max_issues_repo_name: sisl/VirtualTestDrive

max_forks_repo_name: sisl/VirtualTestDrive


text:
using Mimi
include("../src/main_model.jl")

m = getpage()
run(m)
explore(m)

max_stars_repo_name: lrennels/PAGE-ICE

max_issues_repo_name: lrennels/PAGE-ICE

max_forks_repo_name: lrennels/PAGE-ICE


text:
export Linear
using Knet.Ops21: linear
using AutoGrad: Param

"""
    Linear(inputsize, outputsize; winit, binit, activation, dropout)
    Linear(w; bias, inputsize, outputsize, activation, dropout)

Return a function that generalizes matrix multiplication to possibly more than 2 dims with
reshapes:

    w(M...,N...) * x(N...,K...) => y(M...,K...)

where `inputsize=(N...)`, `outputsize=(M...)` each of which can be a single dimension or a
tuple of dimensions. Optionally dropout is applied to the input, a bias of size `(M...)` is
added to the output and an activation function `activation` is applied to the
output. No bias will be used unless `bias` or `binit` are specified.

The first form takes the sizes and initializes `w` and `bias` using the distributions given
by `winit` and `binit` with an array type that matches the first input. The second form
takes existing weights `w` and `bias` and performs no initialization.

Arguments:
* `w`: the linear transformation weights
* `bias=nothing`: optional bias weights
* `winit=𝑼(√(6/(prod(inputsize)+prod(outputsize))))`: distribution for weight initialization
* `binit=nothing`: distribution for bias initialization, `nothing` means no bias
* `inputsize=size(w)[end]`: size of the input tensor (excluding batch, beam etc. dimensions)
* `outputsize=size(w)[1:end-1]`: size of the output tensor (excluding batch, beam etc. dimensions)
* `activation=nothing`: apply activation function to output unless equal to `nothing`
* `dropout=0`: apply dropout with this probability to input if non-zero

References:
* torch.nn.Linear
* tf.keras.layers.Dense
"""
mutable struct Linear
    w
    bias
    winit
    binit
    inputsize
    outputsize
    activation
    dropout
end

function Linear(
    w;
    inputsize=size(w)[end],
    outputsize=(ndims(w)==2 ? size(w,1) : size(w)[1:end-1]),
    bias=nothing,
    activation=nothing,
    dropout=0,
)
    @assert size(w) == (outputsize..., inputsize...) "size(w) must be $((outputsize..., inputsize...))"
    @assert bias === nothing || bsimilar(size(bias), outputsize) "size(bias) must be $(outputsize) not $(size(bias))"
    w = (w isa Param ? w : Param(w))
    bias = (bias isa Nothing || bias isa Param ? bias : Param(bias))
    Linear(w, bias, nothing, nothing, inputsize, outputsize, activation, dropout)
end

function Linear(
    inputsize, outputsize;
    winit=𝑼(√(6/(prod(inputsize)+prod(outputsize)))),
    binit=nothing,
    activation=nothing,
    dropout=0,
)
    Linear(nothing, nothing, winit, binit, inputsize, outputsize, activation, dropout)
end

function (l::Linear)(x)
    initlinear(l, x)
    if l.dropout != 0; x = dropout(x, l.dropout); end
    y = linear(l.w, x)
    if l.bias !== nothing; y = y .+ l.bias; end
    if l.activation !== nothing; y = l.activation(y); end
    return y
end

function initlinear(l::Linear, x)
    if l.w === nothing
        wsize = (l.outputsize..., l.inputsize...)
        l.w = Param(copyto!(similar(x, wsize...), l.winit(eltype(x), wsize...)))
    end
    if l.bias === nothing && l.binit !== nothing
        bsize = l.outputsize
        l.bias = Param(copyto!(similar(x, bsize...), l.binit(eltype(x), bsize...)))
    end
end

function bsimilar(a,b)          # compare dimensions ignoring trailing ones
    dimpair(d) = (d isa Integer ? (d,()) : d === () ? (1,()) : (d[1],d[2:end]))
    a === () && b === () && return true
    a1,a2 = dimpair(a)
    b1,b2 = dimpair(b)
    a1 == b1 && bsimilar(a2,b2)
end

# function (l::Linear)(x::MaskedArray)
#     (a,m) = (x.array, x.mask)
#     @assert m===nothing || all(size(m,i) == 1 || size(m,i) == size(a,i) for i in 1:ndims(a))
#     if m === nothing
#         return MaskedArray(l(a), nothing)
#     elseif size(m,1) == 1   # avoid mask multiplication if possible
#         b = l(a)
#         if ndims(b) > ndims(m)
#             m = reshape(m, ntuple(i->1, ndims(b)-ndims(m))..., size(m)...)
#         end
#         return MaskedArray(b, m)
#     else
#         return MaskedArray(l(a .* oftype(a,m)), nothing)
#     end
# end

##
# In RNNs where input has shape (X,B,T) we want the result to be (Y,B,T). In CNNs where
# input has shape (X1,X2,C,B) we want the output to be (Y,B). So in general we want to be
# able to specify how many initial dims of x should be replaced. Does it ever get replaced
# with more than one dim? In my vaswani it does.


text:
using Test
using JDF
using DataFrames: DataFrame
using PooledArrays
using Missings: allowmissing

@testset "JDF.jl pooledarrays arrays" begin
    df = DataFrame([rand(1:10, 100) for i = 1:3])
    df[!, :x3] = df[:, :x1] |> allowmissing
    df[!, :x4] = string.(df[!, :x2]) |> allowmissing
    df[50, :x3] = missing
    df[66, :x4] = missing

    df[!, :x1] = PooledArray(df[!, :x1])
    df[!, :x2] = PooledArray(string.(df[!, :x2]))
    df[!, :x3] = PooledArray(df[!, :x3])
    df[!, :x4] = PooledArray(df[!, :x4])

    JDF.save(df, "a3pooled.jdf")
    df2 = JDF.load("a3pooled.jdf")
    @test size(df2, 2) == 4
    @test size(df2, 1) == 100
    @time df2[!, :x1] isa PooledVector{Int}
    @time df2[!, :x2] isa PooledVector{String}
    @time df2[!, :x3] isa PooledVector{Union{Missing,Int}}
    @time df2[!, :x4] isa PooledVector{Union{Missing,String}}

    for n in names(df)
        @test isequal(df2[n], df[n])
    end

    rm("a3pooled.jdf", force = true, recursive = true)
end

max_stars_repo_name: xiaodaigh/JDF

max_issues_repo_name: xiaodaigh/JDF

max_forks_repo_name: xiaodaigh/JDF


text:
using HorizonSideRobots
include("BackPath.jl")

function put_marker_in_corners!(r::Robot)
    back_path = BackPath(r)
    sides = [Nord, Ost, Sud, West]
    for side in sides
        while !isborder(r, side)
            move!(r, side)
        end
        putmarker!(r)
    end
    back!(r, back_path)
end

max_stars_repo_name: Mananhina/mirea-progs

max_issues_repo_name: Mananhina/mirea-progs

max_forks_repo_name: Mananhina/mirea-progs