text:
################################################################################
#
#	CONCRETE TYPE
#
#   Lattice{S,B,U} <: AbstractLattice{S,B,U}
#   --> S is the site type (<: AbstractSite{LS,D})
#       --> D is the dimension of embedding space
#       --> LS is the label type of sites
#   --> B is the bond type (<: AbstractBond{LB,N})
#       --> N is the dimension of the Bravais lattice that the bond is located in
#       --> LB is the label type of bonds
#   --> U is the unitcell type (<: AbstractUnitcell{SU,BU})
#       --> SU is the site type of the unitcell (<: AbstractSite{LUS,DU})
#           --> DU is the dimension of embedding space (of the UC)
#           --> LUS is the label type of sites (of the UC)
#       --> BU is the bond type of the unitcell (<: AbstractBond{LUB,NU})
#           --> NU is the dimension of the Bravais lattice that the bond is located in (of the UC)
#           --> LUB is the label type of bonds (of the UC)
#
#   FILE CONTAINS
#       - concrete struct definition
#       - interface implementation
#
################################################################################

################################################################################
#
#   CONCRETE STRUCT DEFINITION
#
################################################################################
mutable struct Lattice{S,B,U} <: AbstractLattice{S,B,U}

    # basis vectors of the Bravais lattice
    lattice_vectors	:: Vector{Vector{Float64}}

    # all sites within the lattice
    sites			:: Vector{S}

    # list of bonds
    bonds			:: Vector{B}

    # unitcell
    unitcell        :: U

end

# export the concrete type
export Lattice

################################################################################
#
#	INTERFACING / ACCESSING LATTICES
#	(functions have to be overwritten by concrete types)
#
################################################################################

# default constructor interface
# used for creation of new lattices
function newLattice(
            ::Type{Lattice{S,B,U}},
            lattice_vectors :: Vector{<:Vector{<:Real}},
            sites           :: Vector{S},
            bonds           :: Vector{B},
            unitcell        :: U
        ) :: Lattice{S,B,U} where {D,N,LS,LB,U,S<:AbstractSite{LS,D},B<:AbstractBond{LB,N}}

    # return a newly created object
    Lattice{S,B,U}(lattice_vectors, sites, bonds, unitcell)
end

# accessing a list of lattice vectors
function latticeVectors(
            lattice :: Lattice{S,B,U}
        ) :: Vector{Vector{Float64}} where {D,N,LS,LB,U,S<:AbstractSite{LS,D},B<:AbstractBond{LB,N}}

    # return the list of lattice vectors
    return lattice.lattice_vectors
end
# setting a list of lattice vectors
function latticeVectors!(
            lattice         :: Lattice{S,B,U},
            lattice_vectors :: Vector{<:Vector{<:Real}}
        ) where {D,N,LS,LB,U,S<:AbstractSite{LS,D},B<:AbstractBond{LB,N}}

    # set the list of lattice vectors
    lattice.lattice_vectors = lattice_vectors
end

# accessing a list of sites
function sites(
            lattice :: Lattice{S,B,U}
        ) :: Vector{S} where {D,N,LS,LB,U,S<:AbstractSite{LS,D},B<:AbstractBond{LB,N}}

    # return the list of sites
    return lattice.sites
end
# setting a list of sites
function sites!(
            lattice :: Lattice{S,B,U},
            sites   :: Vector{S}
        ) where {D,N,LS,LB,U,S<:AbstractSite{LS,D},B<:AbstractBond{LB,N}}

    # set the list of sites
    lattice.sites = sites
end

# accessing a list of bonds
function bonds(
            lattice :: Lattice{S,B,U}
        ) :: Vector{B} where {D,N,LS,LB,U,S<:AbstractSite{LS,D},B<:AbstractBond{LB,N}}

    # return the list of bonds
    return lattice.bonds
end
# setting a list of bonds
function bonds!(
            lattice :: Lattice{S,B,U},
            bonds   :: Vector{B}
        ) where {D,N,LS,LB,U,S<:AbstractSite{LS,D},B<:AbstractBond{LB,N}}

    # set the list of bonds
    lattice.bonds = bonds
end

# accessing the unitcell
function unitcell(
            lattice :: Lattice{S,B,U}
        ) :: U where {D,N,LS,LB,U,S<:AbstractSite{LS,D},B<:AbstractBond{LB,N}}

    # return the unitcell
    return lattice.unitcell
end
# setting the unitcell
function unitcell!(
            lattice  :: Lattice{S,B,U},
            unitcell :: U
        ) :: U where {D,N,LS,LB,U,S<:AbstractSite{LS,D},B<:AbstractBond{LB,N}}

    # set the unitcell
    lattice.unitcell = unitcell
end

function saveLattice(
        lt :: L,
        fn :: AbstractString,
        group :: AbstractString = "lattice"
        ;
        append :: Bool = false
    ) where {LS,LB,D,N,S<:AbstractSite{LS,D},B<:AbstractBond{LB,N},SU,BU,U<:AbstractUnitcell{SU,BU},L<:Lattice{S,B,U}}

    # determine the mode based on if one wants to append stuff
    if append
        mode = "r+"
    else
        mode = "w"
    end

    # determine the site and bond group names
    group_sites    = group*"/sites"
    group_bonds    = group*"/bonds"
    group_unitcell = group*"/unitcell"

    # open the file in mode
    h5open(fn, mode) do ltfile
        # create the top level group
        group_lt = g_create(ltfile, group)
        # save into the attributes in which groups the sites and bonds are stored
        attrs(group_lt)["sites"]    = group_sites
        attrs(group_lt)["bonds"]    = group_bonds
        attrs(group_lt)["unitcell"] = group_unitcell
        # save types into the attributes
        attrs(group_lt)["site_type"]     = string(S)
        attrs(group_lt)["bond_type"]     = string(B)
        attrs(group_lt)["unitcell_type"] = string(U)
        # number of bravais lattice dimensions
        attrs(group_lt)["N"] = Int64(N)
        # save further attributes

        # write bravais lattice vectors
        for i in 1:N
            group_lt["a"*string(i)] = latticeVectors(lt)[i]
        end
    end

    # save sites, bonds and unitcell into the file
    saveSites(sites(lt), fn, group_sites, append=true)
    saveBonds(bonds(lt), fn, group_bonds, append=true)
    saveUnitcell(unitcell(lt), fn, group_unitcell, append=true)

    # return nothing
    return nothing
end

function loadLattice(
        ::Type{LT},
        fn :: AbstractString,
        group :: AbstractString = "lattice"
    ) where {SLT<:AbstractSite,BLT<:AbstractBond,ULT<:AbstractUnitcell,LT<:Union{Lattice{SLT,BLT,ULT},Lattice}}

    # read the attribute data
    attr_data = h5readattr(fn, group)

    # determine the bravais lattice dimension
    N = attr_data["N"]

    # determine the site and bond group names
    group_sites    = attr_data["sites"]
    group_bonds    = attr_data["bonds"]
    group_unitcell = attr_data["unitcell"]

    # load bravais lattice vectors
    lattice_vectors = Vector{Float64}[]
    for i in 1:N
        push!(lattice_vectors, h5read(fn, group*"/a"*string(i)))
    end

    # load the sites
    site_list = loadSites(fn, group_sites)
    # load the bonds
    bond_list = loadBonds(fn, group_bonds)
    # load the unitcell
    unitcell_loaded = loadUnitcell(fn, group_unitcell)

    # determine site and bond type
    S = typeof(site_list[1])
    B = typeof(bond_list[1])
    U = typeof(unitcell_loaded)

    # return new unitcell
    return newLattice(
        Lattice{S,B,U},
        lattice_vectors,
        site_list,
        bond_list,
        unitcell_loaded
    )
end

# convinience function for standard type
function loadLattice(
        fn :: AbstractString,
        group :: AbstractString = "lattice"
    )

    return loadLattice(Lattice, fn, group)
end


text:
using Distributions: Normal, MixtureModel, pdf
using HCubature
using Makie
using AbstractPlotting: textslider
const d = Normal(0, 1)

𝓝(x) = pdf.(d, x)
prior(x, a, b) = (1 - a) * 𝓝(x) + a * 𝓝(x - b)
ϕ(a, b) = 1. / 10.

domain_integral(f) = hcubature(x->f(x...), (0., -5.), (1., 5.))[1]

∏(v) = prod(v)
post(a, b, Xᴺ, β, Z) = ϕ(a, b) * ∏([prior(Xᵢ, a, b)^β for Xᵢ in Xᴺ])/Z

Ew(f, Xᴺ, β, Z) = domain_integral((a, b) -> f(a, b)*post(a, b, Xᴺ, β, Z))
pred(x, Xᴺ, β, Z) = Ew((a, b) -> prior(x, a, b), Xᴺ, β, Z)

function calc_partition(Xᴺ, β)
    domain_integral((a, b) -> ϕ(a, b) * ∏([prior(Xᵢ, a, b)^β for Xᵢ in Xᴺ]))
end

function main_dynamic()
    s_a₀, obs_a₀ = textslider(0. :0.01:1., "a₀", start=0.5)
    s_b₀, obs_b₀ = textslider(-5.:0.01:5., "b₀", start=3.0)
    n=100
    β = 1.
    alspace = 0.:0.01:1.
    blspace = -5.:0.01:5.
    function update!(z::Array{Float32, 2}, a₀, b₀)
        mixture = MixtureModel(
            Normal[
                Normal(0., 1.0),
                Normal(b₀, 1.0)
            ],
            [1 - a₀, a₀]
        )

        Xᴺ = rand(mixture, n)
        Z = calc_partition(Xᴺ, β)
        z = Float32[post(a, b, Xᴺ, β, Z)
                        for a in alspace, b in blspace]
    end

    args_n = (obs_a₀, obs_b₀)
    pos = lift(update!, Node(zeros(Float32, (n,n))), args_n...)
    sur = surface(alspace, blspace, pos)
    axis = sur[Axis]
    axis[:names][:axisnames] = ("a", "b", "post distribution")

    RecordEvents(
        hbox(
            sur,
            vbox(s_a₀, s_b₀),
            parent=Scene()
        ),
        "output"
    )
end

function main_static(a₀, b₀, n)
    β = 1.

    alspace = 0.:0.01:1.
    blspace = -5.:0.01:5.
    mixture = MixtureModel(
        Normal[
            Normal(0., 1.0),
            Normal(b₀, 1.0)
        ],
        [1 - a₀, a₀]
    )
    Xᴺ = rand(mixture, n)
    Z = calc_partition(Xᴺ, β)

    z = [post(a, b, Xᴺ, β, Z) for a in alspace, b in blspace]
    @show domain_integral((a,b) -> post(a, b, Xᴺ, β, Z))
    sur = surface(alspace, blspace, z)
    axis = sur[Axis]
    axis[:names][:axisnames] = ("a", "b", "post distribution")
    sur
end

ext: jl

max_stars_repo_name: terasakisatoshi/juliaExer

max_issues_repo_name: terasakisatoshi/juliaExer

max_forks_repo_name: terasakisatoshi/juliaExer


text:
using Adapt

struct Field{LX, LY, LZ, A, D, G, T, B} <: AbstractDataField{LX, LY, LZ, A, G, T, 3}
                   data :: D
           architecture :: A
                   grid :: G
    boundary_conditions :: B

    function Field{LX, LY, LZ}(data::D, arch::A, grid::G, bcs::B) where {LX, LY, LZ, D, A, G, B}
        T = eltype(grid)
        return new{LX, LY, LZ, A, D, G, T, B}(data, arch, grid, bcs)
    end
end

"""
    Field(LX, LY, LZ, [arch = CPU()], grid,
          [ bcs = FieldBoundaryConditions(grid, (LX, LY, LZ)),
           data = new_data(eltype(grid), arch, grid, (LX, LY, LZ))])

Construct a `Field` on `grid` with `data` on architecture `arch` with
boundary conditions `bcs`. Each of `(LX, LY, LZ)` is either `Center` or `Face` and determines
the field's location in `(x, y, z)`.

Example
=======

```jldoctest
julia> using Oceananigans

julia> ω = Field(Face, Face, Center, CPU(), RectilinearGrid(size=(1, 1, 1), extent=(1, 1, 1)))
Field located at (Face, Face, Center)
├── data: OffsetArrays.OffsetArray{Float64, 3, Array{Float64, 3}}, size: (1, 1, 1)
├── grid: RectilinearGrid{Float64, Periodic, Periodic, Bounded}(Nx=1, Ny=1, Nz=1)
└── boundary conditions: west=Periodic, east=Periodic, south=Periodic, north=Periodic, bottom=ZeroFlux, top=ZeroFlux, immersed=ZeroFlux
```
"""
function Field(LX, LY, LZ,
               arch::AbstractArchitecture,
               grid::AbstractGrid,
               bcs = FieldBoundaryConditions(grid, (LX, LY, LZ)),
               data = new_data(eltype(grid), arch, grid, (LX, LY, LZ)))

    # Make sure architecture is the same as the grid architecture
    # if hasproperty(grid, :architecture)
    #     arch = grid.architecture
    # end

    validate_field_data(LX, LY, LZ, data, grid)

    return Field{LX, LY, LZ}(data, arch, grid, bcs)
end

# Default CPU architecture
Field(LX, LY, LZ, grid::AbstractGrid, args...) = Field(LX, LY, LZ, CPU(), grid, args...)

# Canonical `similar` for AbstractField (doesn't transfer boundary conditions)
Base.similar(f::AbstractField{LX, LY, LZ, Arch}) where {LX, LY, LZ, Arch} = Field(LX, LY, LZ, Arch(), f.grid)

# Type "destantiation": convert Face() to Face and Center() to Center if needed.
destantiate(X) = typeof(X)
destantiate(X::DataType) = X

"""
    Field(L::Tuple, arch, grid, data, bcs)

Construct a `Field` at the location defined by the 3-tuple `L`,
whose elements are `Center` or `Face`.
"""
Field(L::Tuple, args...) = Field(destantiate.(L)..., args...)

#####
##### Special constructors for tracers and velocity fields
#####

"""
    CenterField([arch=CPU()], grid, args...)

Returns `Field{Center, Center, Center}` on `arch`itecture and `grid`.
Additional arguments are passed to the `Field` constructor.
"""
CenterField(arch::AbstractArchitecture, grid, args...) = Field(Center, Center, Center, arch, grid, args...)

"""
    XFaceField([arch=CPU()], grid, args...)

Returns `Field{Face, Center, Center}` on `arch`itecture and `grid`.
Additional arguments are passed to the `Field` constructor.
"""
XFaceField(arch::AbstractArchitecture, args...) = Field(Face, Center, Center, arch, args...)

"""
    YFaceField([arch=CPU()], grid, args...)

Returns `Field{Center, Face, Center}` on `arch`itecture and `grid`.
Additional arguments are passed to the `Field` constructor.
"""
YFaceField(arch::AbstractArchitecture, args...) = Field(Center, Face, Center, arch, args...)

"""
    ZFaceField([arch=CPU()], grid, args...)

Returns `Field{Center, Center, Face}` on `arch`itecture and `grid`.
Additional arguments are passed to the `Field` constructor.
"""
ZFaceField(arch::AbstractArchitecture, args...) = Field(Center, Center, Face, arch, args...)

# Default CPU architectures...
CenterField(grid::AbstractGrid, args...) = CenterField(CPU(), grid, args...)
 XFaceField(grid::AbstractGrid, args...) =  XFaceField(CPU(), grid, args...)
 YFaceField(grid::AbstractGrid, args...) =  YFaceField(CPU(), grid, args...)
 ZFaceField(grid::AbstractGrid, args...) =  ZFaceField(CPU(), grid, args...)

Adapt.adapt_structure(to, field::Field) = Adapt.adapt(to, field.data)


text:
@test Benitor.bubbleSort!( [3, 4, 7, 5, 8, 9, 6, 1, 2] ) == [1, 2, 3, 4, 5, 6, 7, 8, 9]
@test Benitor.bubbleSort!( [3.4, 5.2, 4.3, 1.0, 0.0]  ) == [0.0, 1.0, 3.4, 4.3, 5.2]
@test Benitor.bubbleSort!( [3.4, 5.2, 4.3, 1.0, 0.0, -2.3, -10.3] ) == [-10.3, -2.3, 0.0, 1.0, 3.4, 4.3, 5.2]
@test Benitor.bubbleSort!( [π, 3.4, 4.3, -1.0, -10.3] ) == [-10.3, -1.0, 3.141592653589793, 3.4, 4.3]
@test collect(skipmissing(Benitor.bubbleSort!( [missing, 3.12, -6.66, 7.5] ))) == collect(skipmissing( [missing, -6.66, 3.12, 7.5] ))
@test Benitor.bubbleSort!( [1//3, 1//6, -1//3, 1//0] ) == [-1//3, 1//6, 1//3, 1//0]


text:
# MIT license
# Copyright (c) Microsoft Corporation. All rights reserved.
# See LICENSE in the project root for full license information.

@testset "Visualization" begin
    if get(ENV, "CI", nothing) == "true" && Sys.iswindows() return #OpenGL is unreliable on headless Windows VMs on github. This fails unpredictably and prevents pull requests from being approved. These tests are not essential, so turn them off when running on windows.
    else
        # empty Vis.draw() call to clear Makie @info message:
        # "Info: Makie/Makie is caching fonts, this may take a while. Needed only on first run!"
        Vis.draw()

        # test that this all at least runs
        surf1 = AcceleratedParametricSurface(TestData.beziersurface(), 15)
        surf2 = AcceleratedParametricSurface(TestData.upsidedownbeziersurface(), 15)

        m = (
            leaf(surf1, translation(-0.5, -0.5, 0.0)) ∩
            Cylinder(0.3, 5.0) ∩
            leaf(surf1, Transform(0.0, Float64(π), 0.0, 0.5, -0.5, 0.0))
        )()
        Vis.draw(m)
        @test_nowarn Vis.draw(m)

        m = (
            leaf(surf1, translation(-0.5, -0.5, 0.0)) ∩
            Cylinder(0.3, 5.0) ∩
            leaf(surf2, translation(-0.5, -0.5, 0.0))
        )()
        @test_nowarn Vis.draw(m)

        m = (Cylinder(0.5, 3.0) ∪ Sphere(1.0))()
        @test_nowarn Vis.draw(m)

        m = (Cylinder(0.5, 3.0) ∩ Sphere(1.0))()
        @test_nowarn Vis.draw(m)

        m = (Cylinder(0.5, 3.0) - Sphere(1.0))()
        @test_nowarn Vis.draw(m)

        @test_nowarn Vis.draw(Examples.hex3RGB(),[0 1 0;0 0 1])
        @test_nowarn Examples.drawhexrect()
        @test_nowarn Examples.drawhexneighbors()
    end
end # testset Visualization


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

# To distinguish between upstream changes and changes to the JLL, we use:
#     version = 100 * upstream version + offset
# See C/Coin-OR/coin-or-common.jl for more details.
version = v"0.600.200"

# Collection of sources required to complete build
sources = [
    GitSource(
        "https://github.com/osqp/osqp.git",
        "f9fc23d3436e4b17dd2cb95f70cfa1f37d122c24",
    ),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/osqp
git submodule update --init --recursive
mkdir build
cd build/
cmake \
    -DCMAKE_INSTALL_PREFIX=${prefix} \
    -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} \
    -DCMAKE_BUILD_TYPE=Release \
    ..
make
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms(;experimental=true)

# The products that we will ensure are always built
products = [
    LibraryProduct("libqdldl", :qdldl),
    LibraryProduct("libosqp", :osqp)
]

# Dependencies that must be installed before this package can be built
dependencies = Dependency[
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(
    ARGS,
    "OSQP",
    version,
    sources,
    script,
    platforms,
    products,
    dependencies;
    julia_compat = "1.6",
)

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
#####
##### 'Tupled closure' implementation: 1-tuple, 2-tuple, and then n-tuple by induction
#####

#####
##### Stress divergences
#####

for stress_div in (:∂ⱼ_τ₁ⱼ, :∂ⱼ_τ₂ⱼ, :∂ⱼ_τ₃ⱼ)
    @eval begin
        @inline $stress_div(i, j, k, grid::AbstractGrid, closures::Tuple{<:Any}, clock, U, Ks, args...) =
                    $stress_div(i, j, k, grid, closures[1], clock, U, Ks[1], args...)

        @inline $stress_div(i, j, k, grid::AbstractGrid, closures::Tuple{<:Any, <:Any}, clock, U, Ks, args...) = (
                    $stress_div(i, j, k, grid, closures[1], clock, U, Ks[1], args...)
                  + $stress_div(i, j, k, grid, closures[2], clock, U, Ks[2], args...))

        @inline $stress_div(i, j, k, grid::AbstractGrid, closures::Tuple{<:Any, <:Any, <:Any}, clock, U, Ks, args...) = (
                    $stress_div(i, j, k, grid, closures[1], clock, U, Ks[1], args...)
                  + $stress_div(i, j, k, grid, closures[2], clock, U, Ks[2], args...)
                  + $stress_div(i, j, k, grid, closures[3], clock, U, Ks[3], args...))

        @inline $stress_div(i, j, k, grid::AbstractGrid, closures::Tuple, clock, U, Ks, args...) = (
                    $stress_div(i, j, k, grid, closures[1:2], clock, U, Ks[1:2], args...)
                  + $stress_div(i, j, k, grid, closures[3:end], clock, U, Ks[3:end], args...))
    end
end

#####
##### Tracer flux divergences
#####

@inline ∇_dot_qᶜ(i, j, k, grid::AbstractGrid, closures::Tuple{<:Any}, c, iᶜ, clock, Ks, args...) =
        ∇_dot_qᶜ(i, j, k, grid, closures[1], c, iᶜ, clock, Ks[1], args...)

@inline ∇_dot_qᶜ(i, j, k, grid::AbstractGrid, closures::Tuple{<:Any, <:Any}, c, iᶜ, clock, Ks, args...) = (
        ∇_dot_qᶜ(i, j, k, grid, closures[1], c, iᶜ, clock, Ks[1], args...)
      + ∇_dot_qᶜ(i, j, k, grid, closures[2], c, iᶜ, clock, Ks[2], args...))

@inline ∇_dot_qᶜ(i, j, k, grid::AbstractGrid, closures::Tuple{<:Any, <:Any, <:Any}, c, iᶜ, clock, Ks, args...) = (
        ∇_dot_qᶜ(i, j, k, grid, closures[1], c, iᶜ, clock, Ks[1], args...)
      + ∇_dot_qᶜ(i, j, k, grid, closures[2], c, iᶜ, clock, Ks[2], args...)
      + ∇_dot_qᶜ(i, j, k, grid, closures[3], c, iᶜ, clock, Ks[3], args...))

@inline ∇_dot_qᶜ(i, j, k, grid::AbstractGrid, closures::Tuple, c, iᶜ, clock, Ks, args...) = (
        ∇_dot_qᶜ(i, j, k, grid, closures[1:2], c, iᶜ, clock, Ks[1:2], args...)
      + ∇_dot_qᶜ(i, j, k, grid, closures[3:end], c, iᶜ, clock, Ks[3:end], args...))

#####
##### Utilities
#####

with_tracers(tracers, closure_tuple::Tuple) =
    Tuple(with_tracers(tracers, closure) for closure in closure_tuple)

function calculate_diffusivities!(diffusivity_fields_tuple, closure_tuple::Tuple, args...)
    for (α, closure) in enumerate(closure_tuple)
        @inbounds diffusivity_fields = diffusivity_fields_tuple[α]
        calculate_diffusivities!(diffusivity_fields, closure, args...)
    end
    return nothing
end

#####
##### Support for VerticallyImplicitTimeDiscretization
#####

const EC = AbstractTurbulenceClosure{<:ExplicitTimeDiscretization}
const VIC = AbstractTurbulenceClosure{<:VerticallyImplicitTimeDiscretization}

# Filter explicitly-discretized closures.
@inline z_diffusivity(clo::Tuple{<:EC},        iᶜ, Ks, args...) = tuple(0)
@inline z_diffusivity(clo::Tuple{<:VIC},       iᶜ, Ks, args...) = tuple(z_diffusivity(clo[1], iᶜ, Ks[1], args...))
@inline z_diffusivity(clo::Tuple{<:VIC, <:EC}, iᶜ, Ks, args...) = tuple(z_diffusivity(clo[1], iᶜ, Ks[1], args...))
@inline z_diffusivity(clo::Tuple{<:EC, <:VIC}, iᶜ, Ks, args...) = tuple(z_diffusivity(clo[2], iᶜ, Ks[2], args...))

@inline z_diffusivity(clo::Tuple{<:VIC, <:VIC}, iᶜ, Ks, args...) = tuple(z_diffusivity(clo[1], iᶜ, Ks[1], args...),
                                                                         z_diffusivity(clo[2], iᶜ, Ks[2], args...))

@inline z_diffusivity(clo::Tuple, iᶜ, Ks, args...) = tuple(z_diffusivity(clo[1:2],   iᶜ, Ks[1:2],   args...)...,
                                                           z_diffusivity(clo[3:end], iᶜ, Ks[3:end], args...)...)

@inline z_viscosity(clo::Tuple{<:EC},         Ks, args...) = tuple(0)
@inline z_viscosity(clo::Tuple{<:VIC},        Ks, args...) = tuple(z_viscosity(clo[1], Ks[1], args...))
@inline z_viscosity(clo::Tuple{<:VIC, <:EC},  Ks, args...) = tuple(z_viscosity(clo[1], Ks[1], args...))
@inline z_viscosity(clo::Tuple{<:EC, <:VIC},  Ks, args...) = tuple(z_viscosity(clo[2], Ks[2], args...))

@inline z_viscosity(clo::Tuple{<:VIC, <:VIC}, Ks, args...) = tuple(z_viscosity(clo[1], Ks[1], args...),
                                                                   z_viscosity(clo[2], Ks[2], args...))

@inline z_viscosity(clo::Tuple, Ks, args...) = tuple(z_viscosity(clo[1:2],   Ks[1:2], args...)...,
                                                     z_viscosity(clo[3:end], Ks[3:end], args...)...)

for coeff in (:νᶜᶜᶜ, :νᶠᶠᶜ, :νᶠᶜᶠ, :νᶜᶠᶠ, :κᶜᶜᶠ, :κᶜᶠᶜ, :κᶠᶜᶜ)
    @eval begin
        @inline $coeff(i, j, k, grid, clock, ν::Tuple{C1})     where C1       = $coeff(i, j, k, grid, clock, ν[1])
        @inline $coeff(i, j, k, grid, clock, ν::Tuple{C1, C2}) where {C1, C2} = $coeff(i, j, k, grid, clock, ν[1])   + $coeff(i, j, k, grid, clock, ν[2])
        @inline $coeff(i, j, k, grid, clock, ν::Tuple)                        = $coeff(i, j, k, grid, clock, ν[1:2]) + $coeff(i, j, k, grid, clock, ν[3:end])
    end
end

const ImplicitClosure = AbstractTurbulenceClosure{TD} where TD <: VerticallyImplicitTimeDiscretization
const ExplicitOrNothing = Union{ExplicitTimeDiscretization, Nothing}

@inline combine_time_discretizations(disc) = disc

@inline combine_time_discretizations(::ExplicitTimeDiscretization, ::VerticallyImplicitTimeDiscretization)           = VerticallyImplicitTimeDiscretization()
@inline combine_time_discretizations(::VerticallyImplicitTimeDiscretization, ::ExplicitTimeDiscretization)           = VerticallyImplicitTimeDiscretization()
@inline combine_time_discretizations(::VerticallyImplicitTimeDiscretization, ::VerticallyImplicitTimeDiscretization) = VerticallyImplicitTimeDiscretization()
@inline combine_time_discretizations(::ExplicitTimeDiscretization, ::ExplicitTimeDiscretization)                     = ExplicitTimeDiscretization()

@inline combine_time_discretizations(disc1, disc2, other_discs...) =
    combine_time_discretizations(combine_time_discretizations(disc1, disc2), other_discs...)

@inline time_discretization(closures::Tuple) = combine_time_discretizations(time_discretization.(closures)...)


text:
"""
Generates the integrand for use in the Sauter-Schwab quadrature formulas.

        generate_integrand(kernel, test_local_space, trial_local_space,
                test_chart, trial_chart)

`kernel` is a function that takes two neighborhoods `x` and `y` and returns a
scalar or dyadic value. `kernel` is the part of the integrand that is most
easily described in term of cartesian coordinates, that is non-separable in the
two arguments, and that does not depend on which of the local shape functions
is considered. In boundary element methods it is the Green function (fundamental
solution) of the underlying partial differential equations.
"""
function generate_integrand_uv(kernel, testref, trialref, testel, trialel)

    function k3(u,v)
        out = @SMatrix zeros(3,3)

        x = neighborhood(testel,u)
        y = neighborhood(trialel,v)

        kernelval = kernel(x,y)
        f = testref(x)
        g = trialref(y)

        jx = jacobian(x)
        jy = jacobian(y)
        ds = jx*jy

        return  SMatrix{3,3}([dot(f[i][1], kernelval*g[j][1])*ds for i=1:3, j=1:3])
    end

    return k3
end


text:
using Dorothy
using Formats
using FormatStreams

datapath = joinpath(@__DIR__, "..", "data")

include("trr.jl")
include("xtc.jl")


text:
"""
    node_hash(node, key2hash)

Calculates and returns the hash corresponding to a `Dispatcher` task graph
node i.e. `DispatchNode` using the hashes of its dependencies, input arguments
and source code of the function associated to the `node`. Any available hashes
are taken from `key2hash`.
"""
function node_hash(node::DispatchNode, key2hash::Dict{T,String}) where T
    hash_code = source_hash(node)
    hash_arguments = arg_hash(node)
    hash_dependencies = dep_hash(node, key2hash)

    node_hash = __hash(join(hash_code, hash_arguments, hash_dependencies))
    subgraph_hash = Dict("code" => hash_code,
                         "args" => hash_arguments,
                         "deps" => hash_dependencies)
    return node_hash, subgraph_hash
end

"""
    source_hash(node)

Hashes the lowered representation of the source code of the function
associated with `node`. Useful for `Op` nodes, the other node types
do not have any associated source code.

# Examples
```julia
julia> using DispatcherCache: source_hash
       f(x) = x + 1
       g(x) = begin
               #comment
               x + 1
              end
       node_f = @op f(1)
       node_g = @op g(10)
       # Test
	   source_hash(node_f) == source_hash(node_g)
# true
```
"""
source_hash(node::Op) = begin
    f = node.func
    local code
    try
        code = join(code_lowered(f)[1].code, "\n")
    catch
        code = get_label(node)
        @warn "Cannot hash code for node $(code) (using label)."
    end
    return __hash(code)
end

source_hash(node::DispatchNode) = __hash(nothing)

"""
    arg_hash(node)

Hash the data arguments (in certain cases configuration fields) of the
dispatch `node`.

# Examples
```julia
julia> using DispatcherCache: arg_hash, __hash
       f(x) = println("\$x")
       arg = "argument"
       node = @op f(arg)
       arg_hash(node)
# "d482b7b1b5357c33"

julia> arg_hash(node) == __hash(hash(nothing) + hash(arg) + hash(typeof(arg)))
# true
```
"""
arg_hash(node::Op) = begin
    h = hash(nothing)
    arguments = (arg for arg in node.args if !(arg isa DispatchNode))
    if !isempty(arguments)
        h += mapreduce(+, arguments) do x
                hash(x) + hash(typeof(x))
            end
    end
    kwarguments = ((k,v) for (k,v) in node.kwargs if !(v isa DispatchNode))
    if !isempty(kwarguments)
        h += mapreduce(+, kwarguments) do x
                k, v = x
                hash(k) + hash(v) + hash(typeof(v))
            end
    end
    return __hash(h)
end

arg_hash(node::DataNode) = __hash(node.data)

arg_hash(node::IndexNode) = __hash(node.index)

arg_hash(node::DispatchNode) = __hash(nothing)

"""
    dep_hash(node, key2hash)

Hash the dispatch node dependencies of `node` using their existing hashes if possible.
"""
dep_hash(node::DispatchNode, key2hash) = begin
    h = __hash(nothing)
    nodes = dependencies(node)
    if isempty(nodes)
        return __hash(h)
    else
        for node in nodes
            h *= get(key2hash, node, node_hash(node, key2hash)[1])
        end
        return __hash(h)
    end
end

"""
    __hash(something)

Return a hexadecimal string corresponding to the hash of sum
of the hashes of the value and type of `something`.

# Examples
```julia
julia> using DispatcherCache: __hash
       __hash([1,2,3])
# "f00429a0d65eb7cb"
```
"""
function __hash(something)
    h = hash(hash(typeof(something)) + hash(something))
    return string(h, base=16)
end


text:
Xc2, yc2 = gen_classif(classes=["A", "B", "C"])

@testset "LDA" begin
    m, f = simple_test_classif_prob(BayesianLDA(), Xc2, yc2)
    fp = fitted_params(m, f)
    @test keys(fp) == (:coef, :intercept, :covariance, :means,
                       :priors, :scalings, :xbar, :classes,
                       :explained_variance_ratio)
    infos = info_dict(m)
    @test infos[:input_scitype] == Table(Continuous)
    @test infos[:target_scitype] == AbstractVector{<:Finite}
    @test !isempty(infos[:docstring])
end

@testset "QDA" begin
    m, f = simple_test_classif_prob(BayesianQDA(), Xc2, yc2)
    fp = fitted_params(m, f)
    @test keys(fp) == (:covariance, :means, :priors,
                       :rotations, :scalings)
    infos = info_dict(m)
    @test infos[:input_scitype] == Table(Continuous)
    @test infos[:target_scitype] == AbstractVector{<:Finite}
    @test !isempty(infos[:docstring])
end


text:
"""
    struct ZScore

Represents a z-scoring scheme, akin to `StatsModels.AbstractContrasts`.  Pass as
value in `Dict` as hints to `schema` (or as `contrasts` kwarg for `fit`).

## Examples

Can specify the center and scale values to use:

```
julia> schema((x=collect(1:10), ), Dict(:x => ZScore(; center=5, scale=3)))
StatsModels.Schema with 1 entry:
  x => x(centered: 5 scaled: 3)
```

Or scale will be automatically computed if left out:

```
julia> schema((x=collect(1:10), ), Dict(:x => ZScore()))
StatsModels.Schema with 1 entry:
  x => x(centered: 5.5 scaled: 3.03)
```
"""
struct ZScore
    center::Any
    scale::Any
end

ZScore(; center=nothing, scale=nothing) = ZScore(center, scale)

function ZScore(xs::AbstractArray, zs::ZScore)
    center = _standard(xs, zs.center)
    scale = _standard(xs, zs.scale)
    return ZScore(center, scale)
end

"""
    struct ZScoredTerm{T,C,S} <: AbstractTerm

A lazily z-scored term. A wrapper around an `T<:AbstractTerm` which will
produce scaled values with `modelcols` by subtracting `center` from each
element and then dividing by `scale`.

## Fields

- `term::T`: The wrapped term.
- `center::C`: The center value which is subtracted from the resulting `modelcols`.
- `scale::S`: The scale value which the resulting `modelcols` are divided by.

## Examples

Directly construct with given scale:

```
julia> d = (x=collect(1:10), );

julia> t = concrete_term(term(:x), d)
x(continuous)

julia> ts = ZScoredTerm(t, 3, 5)
x(centered: 3 scaled: 5)

julia> hcat(modelcols(t + ts, d)...)
10×2 Matrix{Float64}:
  1.0  -0.4
  2.0  -0.2
  3.0   0.0
  4.0   0.2
  5.0   0.4
  6.0   0.6
  7.0   0.8
  8.0   1.0
  9.0   1.2
 10.0   1.4
```

Construct with lazy scaling via [`ZScore`](@ref)

```
julia> ts = concrete_term(term(:x), d, ZScore())
x(centered: 5.5 scaled: 3.03)

julia> hcat(modelcols(t + ts, d)...)
10×2 Matrix{Float64}:
  1.0  -1.4863
  2.0  -1.15601
  3.0  -0.825723
  4.0  -0.495434
  5.0  -0.165145
  6.0   0.165145
  7.0   0.495434
  8.0   0.825723
  9.0   1.15601
 10.0   1.4863
```

Or similarly via schema hints:

```
julia> sch = schema(d, Dict(:x => ZScore()))
StatsModels.Schema with 1 entry:
  x => x(centered: 5.5 scaled: 3.03)
```

"""
struct ZScoredTerm{T,C,S} <: AbstractTerm
    term::T
    center::C
    scale::S
end

function StatsModels.concrete_term(t::Term, xs::AbstractArray, z::ZScore)
    return zscore(StatsModels.concrete_term(t, xs, nothing), ZScore(xs, z))
end

# run-time constructors:
function StatsBase.zscore(t::ContinuousTerm, z::ZScore)
    return ZScoredTerm(t, something(z.center, t.mean), something(z.scale, sqrt(t.var)))
end
function StatsBase.zscore(t::ContinuousTerm; center=nothing, scale=nothing)
    return ZScoredTerm(t, center, scale)
end
function StatsBase.zscore(t::AbstractTerm)
    throw(ArgumentError("can only compute z-score for ContinuousTerm; must provide scale value via zscore(t; center, scale)"))
end

function StatsBase.zscore(t::AbstractTerm, z::ZScore)
    z.scale !== nothing && z.center !== nothing ||
        throw(ArgumentError("can only compute z-score for ContinuousTerm; must provide scale via zscore(t; center, scale)"))
    return ZScoredTerm(t, z.center, z.scale)
end

function StatsModels.modelcols(t::ZScoredTerm, d::NamedTuple)
    return zscore(modelcols(t.term, d), t.center, t.scale)
end

function StatsBase.coefnames(t::ZScoredTerm)
    if StatsModels.width(t.term) == 1
        return "$(coefnames(t.term))(centered: $(_round(t.center)) scaled: $(_round(t.scale)))"
    elseif length(t.scale) > 1
        return string.(vec(coefnames(t.term)), "(centered: ", _round.(vec(t.center)),
                       " scaled: ", _round.(vec(t.scale)), ")")
    else
        return string.(coefnames(t.term), "(centered: ", _round(t.center), " scaled: ",
                       _round(t.scale), ")")
    end
end
# coef table: "x(scaled: 5.5)"
function Base.show(io::IO, t::ZScoredTerm)
    return print(io, "$(t.term)(centered: $(_round(t.center)) scaled: $(_round(t.scale)))")
end
# regular show: "x(scaled: 5.5)", used in displaying schema dicts
function Base.show(io::IO, ::MIME"text/plain", t::ZScoredTerm)
    return print(io, "$(t.term)(centered: $(_round(t.center)) scaled: $(_round(t.scale)))")
end
# long show: "x(scaled: 5.5)"

# statsmodels glue code:
StatsModels.width(t::ZScoredTerm) = StatsModels.width(t.term)
# don't generate schema entries for terms which are already scaled
StatsModels.needs_schema(::ZScoredTerm) = false
StatsModels.termsyms(t::ZScoredTerm) = StatsModels.termsyms(t.term)
StatsModels.termvars(t::ZScoredTerm) = StatsModels.termvars(t.term)


text:
module JuMPTutorials

end


text:
reddit_init() = register(DataDep(
    "Reddit",
    """
    The Reddit dataset from the `"Inductive Representation Learning on
    Large Graphs" <https://arxiv.org/abs/1706.02216>`_ paper, containing
    Reddit posts belonging to different communities.
    """,
    "https://data.dgl.ai/dataset/reddit.zip",
    "9a16353c28f8ddd07148fc5ac9b57b818d7911ea0fbe9052d66d49fc32b372bf";
    post_fetch_method=preprocess_reddit,
))

function preprocess_reddit(local_path)
    unzip_zip(local_path)

    graph_file = datadep"Reddit/reddit_graph.npz"
    data_file = datadep"Reddit/reddit_data.npz"

    graphfile = replace(local_path, "reddit.zip"=>"reddit.graph.jld2")
    allfile = replace(local_path, "reddit.zip"=>"reddit.all.jld2")
    rawfile = replace(local_path, "reddit.zip"=>"reddit.raw.jld2")
    metadatafile = replace(local_path, "reddit.zip"=>"reddit.metadata.jld2")

    graph, X, y = to_reddit_rawfile(graph_file, data_file, rawfile)

    sg = to_simplegraph(graph)
    all_X = Matrix(X')
    all_y = Matrix{UInt16}(y')
    meta = (graph=(num_V=nv(sg), num_E=ne(sg)),
            all=(features_dim=size(all_X), labels_dim=size(all_y))
            )

    @save graphfile sg
    @save allfile all_X all_y
    @save metadatafile meta
end

function to_reddit_rawfile(graph_file, data_file, rawfile)
    py"""
    import numpy as np
    import scipy.sparse as sp
    graph = np.load($graph_file, allow_pickle=True)
    data = np.load($data_file, allow_pickle=True)
    """

    graph = sparse(Vector(py"graph['row']") .+ 1,
                   Vector(py"graph['col']") .+ 1,
                   Vector{Int32}(py"graph['data']"))
    X = Matrix{Float32}(py"data['feature']")
    y = Vector{Int32}(py"data['label']")
    ids = Vector{Int32}(py"data['node_ids']")
    types = Vector{UInt8}(py"data['node_types']")
    raw = Dict(:graph=>graph, :X=>X, :y=>y, :ids=>ids, :types=>types)

    @save rawfile raw

    graph, X, y
end

struct Reddit <: Dataset
end

function graphdata(::Reddit)
    file = datadep"Reddit/reddit.graph.jld2"
    @load file sg
    sg
end

function alldata(::Reddit)
    file = datadep"Reddit/reddit.all.jld2"
    @load file all_X all_y
    all_X, all_y
end

function rawdata(::Reddit)
    file = datadep"Reddit/reddit.raw.jld2"
    @load file raw
    raw
end

function metadata(::Reddit)
    file = datadep"Reddit/reddit.metadata.jld2"
    @load file meta
    meta
end


text:
"""
Nonnegativity-constrained proximal operator of nnz()
    nnz_y = proxl0nonneg!( x, a, y )
solves the following problem
    min_y    a nnz(y) + 1/2 ||y-x||^2    s.t.    y ≥ 0
where `a` is a given nonnegative scalar.
The analytical expression is available:
    y = { x  if  x ≥ √2a,  0  otherwise
"""
function proxl0nonneg!( x::Vector{R}, a::R, y::Vector{R} ) where {R <: Real}
    @assert a >= R(0)
    y .= x
    if a > R(0)
        y[.!(x .> sqrt( 2 * a ))] .= R(0)
    else
        y .= max.( R(0), x )
    end
    nnz_y = sum(y .> R(0))
    return nnz_y
end

"""
Simplex-constrained proximal operator of nnz()
    nnz_y = proxl0simplex!( x, a, b, y )
solves the following problem
    min_y    a nnz(y) + 1/2 ||y-x||^2    s.t.    y ≥ 0,    1 ⋅ y = b
where `a` and `b` are given nonnegative scalars.
"""
function proxl0simplex!( x::Vector{R}, a::R, b::R, y::Vector{R} ) where {R <: Real}
    @assert a >= R(0)
    @assert b >= R(0)
    if a == R(0)
        projsimplex!(x, b, y)
        nnz_y = sum(y .> R(0))
        return nnz_y
    end
    n = length(x)
    y .= x # max.( x, 0.0 ) # project
    prm = sortperm(y) # sort
    y .= y[prm]
    mv = Vector(0:n-1) # number of zeros
    lv = (b .- reverse(cumsum(reverse(y)))) ./ (n .- mv)
    feas = (y .+ lv .> R(0)) # feasibility
    if !any(feas)
        m = n - 1
        l = b - x[n]
        #c = a + 0.5 * sum( x[1:m].^2 ) + 0.5 * (l ^ 2) # cost
    else
        sv = cumsum([R(0); y[1:n-1]] .^ 2)
        mv = mv[feas]
        lv = lv[feas]
        sv = sv[feas]
        cv = a .* (n .- mv) .+ 0.5 .* sv .+ 0.5 .* (n .- mv) .* (lv .^ 2) # cost
        c, i = findmin(cv) # find minimum
        m = mv[i]
        l = lv[i]
    end
    y[1:m] .= R(0)
    y[m+1:n] .+= l
    nnz_y = n - m
    y .= y[invperm(prm)] # unsort
    return nnz_y
end

"""
Projection onto the simplex
    projsimplex!(x, b, y)
corresponds to solving the following problem
        min_y    1/2 ||y-x||^2    s.t.    y ≥ 0,    1 ⋅ y = b
where `b` is a given nonnegative scalar.
See arxiv.org/abs/1101.6081
"""
function projsimplex!(x::Vector{R}, b::R, y::Vector{R}) where {R <: Real}
    @assert b >= R(0)
    if b == R(0)
        y .= R(0)
        return nothing
    end
    n = length(x)
    y .= sort(x, rev=true) # sort
    flag = false
    tmpsum = R(0)
    for i in 1:n-1
        tmpsum += y[i] # cumulative sum
        tmpmax = (tmpsum - b) / i
        if tmpmax >= y[i+1]
            flag = true
            break
        end
    end
    if !flag
        tmpmax = (tmpsum + y[n] - b) / n
    end
    y .= max.(x .- tmpmax, R(0)) # shift, project
    return nothing
end


text:
# test default solvers

include("linprog.jl")
linprogtest()

include("mixintprog.jl")
mixintprogtest()

include("quadprog.jl")
quadprogtest()

# Test conic fallback for LPs
include("conicinterface.jl")
coniclineartest(MathProgBase.defaultLPsolver, duals=true)

# Test LP fallback for conics
include("linproginterface.jl")
linprogsolvertest(MathProgBase.defaultConicsolver)


text:
module HighlyAdaptiveLasso

using MLJModelInterface
using MLJModelInterface: @mlj_model
using RCall
using DataFrames

################################
########## Structure ###########
################################

DESCR =
"""
    HALRegressor/HALClassifier(;formula=nothing,
                 max_degree=nothing,
                 smoothness_orders=nothing,
                 num_knots=nothing,
                 reduce_basis=nothing,
                 lambda=nothing,
                 cv_select=nothing,
                 n_folds=nothing)

An interface respecting the MMI for a deterministic/probabilistic regression/classification model that wraps the
HAL implementation from [hal9001](https://github.com/tlverse/hal9001/)m version 0.4.1.

# Arguments

Only some hyper parameters are exposed for now. Refer to the original repository for more details on possible values.

- formula: A String, eg ~h(x1) + h(x2, x1) + h(x3)
- max_degree: Integer
- smoothness_orders: Integer
- num_knots: Integer, Vector, ...
- reduce_basis: Float64
- lambda: Float64 or Vector of such if running in cross validation mode
- cv_select: Bool
- n_folds: Integer

"""

"$DESCR"
@mlj_model mutable struct HALRegressor <: MLJModelInterface.Deterministic
    formula::Union{String, Nothing} = nothing
    max_degree::Union{Int, Nothing} = nothing
    smoothness_orders::Union{Int, Nothing} = nothing
    num_knots = nothing
    reduce_basis::Union{Float64, Nothing} = nothing
    family::String = "gaussian"
    lambda = nothing

    cv_select::Union{Bool, Nothing} = nothing
    n_folds::Union{Int, Nothing} = nothing
end

"$DESCR"
@mlj_model mutable struct HALClassifier <: MLJModelInterface.Probabilistic
    formula = nothing
    max_degree::Union{Int, Nothing} = nothing
    smoothness_orders::Union{Int, Nothing} = nothing
    num_knots = nothing
    reduce_basis::Union{Float64, Nothing} = nothing
    family::String = "binomial"
    lambda = nothing

    cv_select::Union{Bool, Nothing} = nothing
    n_folds::Union{Int, Nothing} = nothing
end

HAL = Union{HALRegressor, HALClassifier}

make_fitresult(m::HAL, fitresult, y) = (fitresult, y[1])

adapt_y(m::HALRegressor, y) = y
adapt_y(m::HALClassifier, y) = int(y) .- 1

"""
    fit(m::HAL, verbosity::Int, Xt, y)

"""
function MLJModelInterface.fit(m::HAL, verbosity::Int, Xt, yt)
    R"library(hal9001)"

    X = DataFrame(Xt)
    @rput X

    y = adapt_y(m, yt)
    @rput y

    # I can't make dynamical assignement of variables work
    # So I rely on this ugly iteraction over variables
    # to build the string for RCall
    fitstring = "fit_hal(X, y, family = '$(m.family)'"
    if m.formula !== nothing
        R"formula <- $(m.formula)"
        fitstring *= ", formula = formula"
    end
    if m.max_degree !== nothing
        R"max_degree <- $(m.max_degree)"
        fitstring *= ", max_degree = max_degree"
    end
    if m.smoothness_orders !== nothing
        R"smoothness_orders <- $(m.smoothness_orders)"
        fitstring *= ", smoothness_orders = smoothness_orders"
    end
    if m.num_knots !== nothing
        R"num_knots <- $(m.num_knots)"
        fitstring *= ", num_knots = num_knots"
    end
    if m.reduce_basis !== nothing
        R"reduce_basis <- $(m.reduce_basis)"
        fitstring *= ", reduce_basis = reduce_basis"
    end
    if m.lambda !== nothing
        R"lambda <- $(m.lambda)"
        fitstring *= ", lambda = lambda"
    end
    # Fit control string
    if m.cv_select !== nothing || m.n_folds !== nothing
        fitstring *= ", fit_control = list("
        if m.cv_select !== nothing
            R"cv_select <- $(m.cv_select)"
            fitstring *= "cv_select = cv_select, "
        end
        if m.n_folds !== nothing
            R"n_folds <- $(m.n_folds)"
            fitstring *= "n_folds = n_folds, "
        end
        fitstring = fitstring[1:end - 2] * ")"
    end
    fitstring *= ")"
    # Rcall

    fitresult = reval(fitstring)

    cache = nothing
    report = fitstring
    return (make_fitresult(m, fitresult, yt), cache, report)
end

function MLJModelInterface.fitted_params(m::HALRegressor, fitresult)
    return rcopy(fitresult)
end

function MLJModelInterface.predict(m::HALClassifier, (fitresult, decode), Xnew)
    Xnew = DataFrame(Xnew)
    @rput Xnew
    ypred = R"predict($fitresult, new_data = $Xnew)"
    μ = rcopy(ypred)
    return UnivariateFinite(classes(decode), μ, augment=true)
end

function MLJModelInterface.predict(m::HALRegressor, (fitresult, _), Xnew)
    Xnew = DataFrame(Xnew)
    @rput Xnew
    ypred = R"predict($fitresult, new_data = $Xnew)"
    return rcopy(ypred)
end

################################
########### METADATA ###########
################################

MLJModelInterface.metadata_pkg.((HALRegressor, HALClassifier),
    name       = "HighlyAdaptiveLasso",
    uuid       = "c5dac772-1445-43c4-b698-9440de7877f6",
    url        = "https://github.com/olivierlabayle/HighlyAdaptiveLasso.jl",
    julia      = true,
    license    = "MIT",
    is_wrapper = false,
)

MLJModelInterface.metadata_model(HALRegressor,
    input_scitype    = MLJModelInterface.Table(MLJModelInterface.Continuous),
    target_scitype   = AbstractVector{MLJModelInterface.Continuous},
    supports_weights = false,
    descr            = DESCR,
	load_path        = "HighlyAdaptiveLasso.HALRegressor"
)

MLJModelInterface.metadata_model(HALClassifier,
    input_scitype    = MLJModelInterface.Table(MLJModelInterface.Continuous),
    target_scitype   = AbstractVector{<:Finite},
    supports_weights = false,
    descr            = DESCR,
	load_path        = "HighlyAdaptiveLasso.HALClassifier"
)

################################
########### Exports ############
################################

export HALRegressor, HALClassifier
export fit, fitted_params, predict

end


text:
__precompile__()
module RegTools
using GLM, DataFrames,Gadfly, Distributions
export avPlot,vif,rsquared,adjrsquared,rstudent,jackknife,halfnorm,cooksdistance

include("diagnostics.jl")
include("misc.jl")
include("modsel.jl")

end # module

max_stars_repo_name: wakakusa/RegTools

max_issues_repo_name: wakakusa/RegTools

max_forks_repo_name: wakakusa/RegTools


text:
#=
The data argument provides a collection of data to train with
(usually a set of inputs x and target outputs y)
For example, here's a dummy data set with only one data point:
=#

x = rand(28*28)
y = rand(10)
data = [(x,y)]
@show data

data = [(x,y),(x,y),(x,y)]
# or equivalently
data = Iterators.repeated((x,y),3)

xs = [rand(784),rand(784),rand(784)]
ys = [rand( 10),rand( 10),rand( 10)]

data = zip(xs,ys)

using Flux: @epochs

@epochs 2 println("Hello")


text:
# --------10--------20--------30--------40--------50--------60--------70--------80--------90--------100-------110-------120--------130
struct Time{TF <: AbstractFloat,TI <: Integer}
    nstep_cyc::TI
    dt::TF   # real time propagation
    dh::Complex{TF}   # imaginary time-propagation
    nstep_tot::TI # num. of time step
    ncyc_tot::TI # num. of optical cycle
    tcyc::Vector{TF} # time step in optical cycle of the first laser
    t::Vector{TF} # time step in atomic unit
end
function Time(λ, nstep_cyc, ncyc_tot)
    T = 2π / wlen2au(λ)
    nstep_cyc = nstep_cyc
    dt = T / nstep_cyc
    dh = dt / im
    nstep_tot = nstep_cyc * ncyc_tot
    tcyc = collect(range(0.0, ncyc_tot, length=nstep_tot))
    t = tcyc * T
    return Time(nstep_cyc, dt, dh, nstep_tot, ncyc_tot, tcyc, t)
end
function Time(pulse::Pulse, nstep_cyc, ncyc_tot)
    @unpack λ = pulse
    return Time(λ, nstep_cyc, ncyc_tot)
end


text:
using Donquicote, Test, BenchmarkTools
using Roots, ForwardDiff

f(x) = exp(x) - x^4
∂f(x) = exp(x) - 3*x^4
∂2f(x) = exp(x) - 2*x^3
D(f) = x -> ForwardDiff.derivative(f,float(x))

println("Bisection (built-in):")
@btime find_zero(f, (1,3), Bisection(); xatol = 1e-7)
println("FalsePosition (build-in):")
@btime find_zero(f, (1,3), FalsePosition(), xatol = 1e-7)
println("Newton (built-in):")
@btime find_zero((f, D(f)), 2, Roots.Newton(); xatol = 1e-7)
println("Newton (custom):")
@btime newton(f, D(f), 2)
println("Halley (built-in):")
@btime find_zero((f, D(f), D(D(f))), 2, Roots.Halley(); xatol = 1e-7)
println("Halley (custom):")
@btime haley(f, D(f), D(D(f)), 2)
println("Schroeder (built-in):")
@btime find_zero((f, D(f), D(D(f))), 2, Roots.Schroeder(); xatol = 1e-7)
println("Schroeder (custom):")
@btime schroeder(f, D(f), D(D(f)), 2)

##

V = ones(3)
V0 = zeros(V)


text:
# p6u.jl - variable coefficient wave equation - UNSTABLE VARIANT

# Grid, variable coefficient, and initial data:
N = 128; h = 2*pi/N; x = h*(1:N);
c = @. .2 + sin(x-1)^2;
t = 0; dt = 1.9/N;
v = @. exp(-100*(x-1)^2);
vold = @. exp(-100*(x-.2*dt-1)^2);

# Time-stepping by leap frog formula:
tmax = 8; tplot = .15; clf();
plotgap = round(Int,tplot/dt);
nplots = round(Int,tmax/tplot);
data = [v zeros(N,nplots)]; tdata = t;
for i = 1:nplots
    for n = 1:plotgap
        t = t+dt;
        v_hat = fft(v);
        w_hat = 1im*[0:N/2-1; 0; -N/2+1:-1] .* v_hat;
        w = real.(ifft(w_hat));
        vnew = vold - 2*dt*c.*w;        # leap frog formula
        vold = v; v = vnew;
    end
    data[:,i+1] = v; tdata = [tdata; t];
    if norm(v,Inf)>2.5
        data = data[:,1:i+1];
        break
    end
end

# Plot results:
mesh(x,tdata,data',ccount=0);
xlim(0,2*pi); ylim(0,tmax); zlim(-3,3);
xlabel("x"); ylabel("t"); zlabel("u"); gca()[:view_init](70,10-90);


text:
using Test
using DFTK
import DFTK: mpi_sum

include("testcases.jl")
using Random
Random.seed!(0)

function test_consistency_term(term; rtol=1e-3, atol=1e-8, ε=1e-8, kgrid=[1, 2, 3],
                               lattice=silicon.lattice, Ecut=10, spin_polarization=:none)
    sspol = spin_polarization != :none ? " ($spin_polarization)" : ""
    @testset "Hamiltonian consistency $(typeof(term)) $sspol" begin
        Si = ElementPsp(14, psp=load_psp(silicon.psp))
        atoms = [Si => silicon.positions]
        model = Model(lattice; n_electrons=silicon.n_electrons, atoms=atoms, terms=[term],
                      spin_polarization=spin_polarization)
        basis = PlaneWaveBasis(model, Ecut; kgrid=kgrid, use_symmetry=false)

        n_electrons = silicon.n_electrons
        n_bands = div(n_electrons, 2)
        filled_occ = DFTK.filled_occupation(model)

        ψ = [Matrix(qr(randn(ComplexF64, length(G_vectors(basis.kpoints[ik])), n_bands)).Q)
             for ik in 1:length(basis.kpoints)]
        occupation = [filled_occ * rand(n_bands) for ik in 1:length(basis.kpoints)]
        occ_scaling = n_electrons / sum(sum(occupation))
        occupation = [occ * occ_scaling for occ in occupation]
        ρ, ρspin = compute_density(basis, ψ, occupation)

        dψ = [randn(ComplexF64, size(ψ[ik])) for ik = 1:length(basis.kpoints)]
        ψ_trial = ψ .+ ε .* dψ
        ρ_trial, ρspin_trial = compute_density(basis, ψ_trial, occupation)

        @assert length(basis.terms) == 1
        E0, ham = energy_hamiltonian(basis, ψ, occupation; ρ=ρ, ρspin=ρspin)
        E1, _ = energy_hamiltonian(basis, ψ_trial, occupation; ρ=ρ_trial, ρspin=ρspin_trial)
        diff = (E1.total - E0.total)/ε

        diff_predicted = 0.0
        for (ik, kpt) in enumerate(basis.kpoints)
            Hψ = ham.blocks[ik]*ψ[ik]
            dψHψ = sum(occupation[ik][iband] * real(dot(dψ[ik][:, iband], Hψ[:, iband]))
                       for iband=1:n_bands)

            diff_predicted += 2 * basis.kweights[ik] * dψHψ
        end
        diff_predicted = mpi_sum(diff_predicted, basis.comm_kpts)

        err = abs(diff - diff_predicted)
        @test err < rtol * abs(E0.total) || err < atol
    end
end

test_consistency_term(Kinetic())
test_consistency_term(AtomicLocal())
test_consistency_term(AtomicNonlocal())
test_consistency_term(ExternalFromReal(X -> cos(X[1])))
test_consistency_term(ExternalFromFourier(X -> randn()))
test_consistency_term(PowerNonlinearity(1.0, 2.0))
test_consistency_term(Hartree())
test_consistency_term(Ewald())
test_consistency_term(PspCorrection())
test_consistency_term(Xc(:lda_xc_teter93))
test_consistency_term(Xc(:lda_xc_teter93), spin_polarization=:collinear)
test_consistency_term(Xc(:gga_x_pbe), spin_polarization=:collinear)

let
    a = 6
    pot(x, y, z) = (x - a/2)^2 + (y - a/2)^2
    Apot(x, y, z) = .2 * [y - a/2, -(x - a/2), 0]
    Apot(X) = Apot(X...)
    test_consistency_term(Magnetic(Apot);  kgrid=[1, 1, 1], lattice=[a 0 0; 0 a 0; 0 0 0], Ecut=20)
    test_consistency_term(DFTK.Anyonic(1); kgrid=[1, 1, 1], lattice=[a 0 0; 0 a 0; 0 0 0], Ecut=20)
end


text:
"""
    row_countlines(io::IO)

    Count number of rows in a file, also deteremines EOL character/padding from first line.
    Test last line of file to see if it is 0 or len bytes
"""
# To support bytes or characters this no longer does malformed testing
function row_countlines(io::IO)
    rows = 0
    # EOL detection, if we don't have an EoL doesn't matter
    start_pos = position(io)
    line = readline(io, keep=false)
    eolpad = (eof(io) || length(line) < 2 || (line[end-1] != '\r')) ? 0 : 1
    seek(io, start_pos)
    rows = mod_countlines(io)
    return (rows, eolpad)
end

# version of countlines() that checks last line for non-empty.
function mod_countlines(io::IO)
    b=[UInt8(0)]::Vector{UInt8}
    eof(io) && return 0
    l = countlines(io)
    readbytes!(skip(io, -1), b, 1)
    (b[1] == UInt8('\n')) ? l : l+1
end

"""
    FWF.readsplitline!(vals::Vector{Union{Missing,String}}, source::FWF.Source)
    FWF.readsplitline!(vals::Vector{Union{Missing,String}}, io::IO, columnwidths::Vector{UnitRange{Int}}, trim::Bool, errorlevel::Symbol)

Read next line from a `FWF.Source` or `IO` as a `Vector{String}` and
store the values in `vals`.
Fields are determined by the field widths stored in the `source` options or `columnwidths`
Fields will be trimed if `trim` is true
Row or rows will be skipped if there is an error found if `errorlevel` is `:skip`
and as much as possible is parsed if it is `:parse`.
The contents of `vals` are replaced.
Returns length of line read in bytes or characters following `unitbytes` keyword argument.
"""
# This function is pretty simple
# * Read line of input from the source
# * Ensure it meets specifications
# * Break it into chunks based on column widths

function readsplitline!(vals::Vector{Union{Missing,String}}, source::FWF.Source)
    line_len = readsplitline!(vals, source.io, source.options.columnrange,
                              source.options.trimstrings, source.options.unitbytes,
                              source.options.errorlevel)
    if source.line_len == -1
        source.line_len = line_len
    elseif source.line_len != line_len
        println(STDERR, "Inconsistent line lengths. Current $(line_len) and first valid length was $(source.line_len)")
    end
    return line_len
end

function readsplitline!(vals::Vector{Union{Missing,String}}, io::IO, columnwidths::Vector{UnitRange{Int}}, trim::Bool=true, unitbytes=true, errorlevel=:parse)
    empty!(vals)
    # Parameter validation
    ((columnwidths == nothing) || (isempty(columnwidths))) && throw(ArgumentError("No column widths provided"))
    eof(io) && (throw(ArgumentError("IO not available")))

    our_length = unitbytes ? sizeof : length

    rowlength = last(last(columnwidths))
    # Read a line and validate
    test = true
    line = ""
    line_len::Int = -1

    while test
        eof(io) && (throw(ArgumentError("Unable to find next valid line")))
        line = readline(io)
        line_len = our_length(line)
        # we do not care what is in the line beyond rowlength
        if line_len < rowlength
            !(errorlevel in (:parse, :skip)) && throw(ParsingException("Invalid length line: "*string(our_length(line))))
            println(STDERR, "Invalid length line($(our_length(line))):", line)
            if errorlevel == :parse
                test = false
            end
        else
            test = false
        end
    end

    chr = 0
    ind = 0
    # Break it up into chunks
    for range in columnwidths
        if errorlevel == :parse && first(range) > line_len
            push!(vals, missing)
            continue
        end
        if unitbytes || isascii(line)
            # Julia 0.7 is more strict about parsing string ranges
            # below will fail if first(range) is not a valid index into line
            # TODO: simplify when only Julia 0.7 or higher is supported
            str = line[first(range):prevind(line, last(range)+1)]
        else
            # TODO: this code will be simpler when only Julia 0.7 or higher is supported
            while chr < first(range)
                ind = nextind(line, ind)
                chr += 1
            end
            ind_start = ind
            while chr < last(range)
                ind = nextind(line, ind)
                chr += 1
            end
            if errorlevel == :parse
                # TODO: in Julia 1.0 we will have to replace endof
                str = line[ind_start:min(ind, endof(line))]
            else
                str = line[ind_start:ind]
            end
        end
        # strip returns SubString in Julia 0.7 or higher
        push!(vals, trim ? String(strip(str)) : str)
    end
    return line_len
end

"""
`FWF.read(fullpath::Union{AbstractString,IO}, columnwidths::Union{Vector{UnitRange{Int}}, Vector{Int}}, sink::Type{T}=DataFrame, args...; kwargs...)` => `typeof(sink)`

`FWF.read(fullpath::Union{AbstractString,IO}, columnwidths::Union{Vector{UnitRange{Int}}, Vector{Int}}, sink::Data.Sink; kwargs...)` => `Data.Sink`

parses a fixed width file into a Julia structure (a DataFrame by default, but any valid `Data.Sink` may be requested).

Positional arguments:

* `fullpath`; can be a file name (string) or other `IO` instance
* `columnwidths`; can be a vector of integers or consecutive unit ranges that represents the column widths
  examples: [4,4,8] or [1:4,5:8, 9:16]; ranges may be discontinous, e.g. [1:2, 5:10].
  Widths are interpreted either in bytes or characters as defined by `unitbytes` keyword argument
  Each valid line in the file must have its width at least equal to width implied by passed `columnwidths`
  (if it is wider it is silently ignored).
* `sink::Type{T}`; `DataFrame` by default, but may also be other `Data.Sink` types that support streaming via
  `Data.Field` interface; note that the method argument can be the *type* of `Data.Sink`, plus any required
  arguments the sink may need (`args...`) or an already constructed `sink` may be passed (2nd method above)

Keyword Arguments:

* `unitbytes::Bool`: whether to treat field ranges in bytes or characters; default = true (bytes)
* `usemissings::Bool`: whether to use missings, all fields will be unioned with Missing; default = true
                        if not set default values of 0, Date() and "" will be used for missing values
* `trimstrings::Bool`: trim whitespace from all strings; default = true
* `errorlevel`  : if `:parse` then as much as possible is parsed and missing data is replaced by `missing`;
                if `:skip` then malformed line is skipped on error;
                if `:error` then an exception is thrown on error; default `:parse`
* `use_mmap::Bool=true`: whether the underlying file will be mmapped or not while parsing; note that on Windows machines, the underlying file will not be "deletable" until Julia GC has run (can be run manually via `gc()`) due to the use of a finalizer when reading the file.
* `skip::Int`: number of rows at start of file to skip; default = 0
* `rows::Int`: maximum number of rows to read from file; default = 0 (whole file)
* `types`: a vector of how to parse each column. (String, Int, Float64, Missing) are valid types, missing will convert whole comumn to `missing`. Pass in the format for Date columns as DateFormat("")
            example: [String, Int, DateFormat("mmddyyyy")]
* `header`: column names can be provided as a Vector{String} or parameter can be set to `true` to use the first row as values or `false` to auto-generate names
* `missings`: a Vector{String} that represents all values that should be converted to missing; example: ["***", "NA", "NULL", "####"]
* `append::Bool=false`: if the `sink` argument provided is an existing table, `append=true` will append the source's data to the existing data instead of doing a full replace
* `transforms::Dict{Union{String,Int},Function}`: a Dict of transforms to apply to values as they are parsed. Note that a column can be specified by either number or column name.

Simple example using inline text:
```
julia> FWF.read(IOBuffer("abc12310102017\ndef45610112017\n"), [3,3,8], types=[String,Int,DateFormat("mmddyyyy")])
```
Example usage from test dataset:
```
# Setup data frame with control information:
julia> format = DataFrame(["PSEUDO-ID" 9 Int;
    "EMPLOYEE_NAME" 23 Missing;
    "FILE_DATE" 8  DateFormat("yyyymmdd");
    "AGENCY" 2 String;
    "SUB_AGENCY" 2 String;
    "DUTY_STATION" 9 String;
    "AGE" 6 String;
    "EDUCATION_LEVEL" 2 String;
    "PAY_PLAN" 2 String;
    "GRADE" 2 String;
    "LOS_LEVEL" 6 String;
    "OCCUPATION" 4 String;
    "OCCUPATIONAL_CATEGORY" 1 String;
    "ADJUSTED_BASIC_PAY" 6 Int;
    "SUPERVISORY_STATUS" 1 String;
    "TYPE_OF_APPOINTMENT" 2 String;
    "WORK_SCHEDULE" 1 String;
    "NSFTP_IND"  1 String])

#Setup NA values
julia> strrep(s::String, r::UnitRange) = [repeat(s, n) for n in r]
julia> naValues = vcat(strrep("*", 1:23), strrep("#", 1:23), "NAME WITHHELD BY AGENCY", "NAME WITHHELD BY OPM", "NAME UNKNOWN", "UNSP", "<NA>", "000000", "999999", "")

julia> dt = FWF.read("sal.txt", convert(Array{Int},format[:x2]),
        header=convert(Array{String}, format[:x1]), types=convert(Array{Union{Type, DateFormat}},format[:x3]),
        missings=naValues)
200000×18 DataFrames.DataFrame. Omitted printing of 8 columns
│ Row    │ PSEUDO-ID │ EMPLOYEE_NAME │ FILE_DATE  │ AGENCY │ SUB_AGENCY │ DUTY_STATION │ AGE   │ EDUCATION_LEVEL │ PAY_PLAN │ GRADE │
├────────┼───────────┼───────────────┼────────────┼────────┼────────────┼──────────────┼───────┼─────────────────┼──────────┼───────┤
│ 1      │ 278418    │ missing       │ 1981-03-31 │ AA     │ 00         │ 110010001    │ 35-39 │ 07              │ GS       │ 09    │
│ 2      │ 485025    │ missing       │ 1981-03-31 │ AA     │ 00         │ 110010001    │ 40-44 │ missing         │ GS       │ 13    │
│ 3      │ 990825    │ missing       │ 1981-03-31 │ AA     │ 00         │ 110010001    │ 25-29 │ 15              │ GS       │ 15    │
│ 4      │ 1220286   │ missing       │ 1981-03-31 │ AA     │ 00         │ 110010001    │ 35-39 │ missing         │ SR       │ 00    │
⋮
│ 199996 │ 6885054   │ missing       │ 1981-03-31 │ DJ     │ 02         │ missing      │ 20-24 │ 07              │ GS       │ 04    │
│ 199997 │ 6885056   │ missing       │ 1981-03-31 │ DJ     │ 02         │ missing      │ 30-34 │ 13              │ GS       │ 10    │
│ 199998 │ 6885057   │ missing       │ 1981-03-31 │ DJ     │ 02         │ missing      │ 30-34 │ 22              │ GS       │ 10    │
│ 199999 │ 6885058   │ missing       │ 1981-03-31 │ DJ     │ 02         │ missing      │ 30-34 │ 18              │ GS       │ 10    │
│ 200000 │ 6885059   │ missing       │ 1981-03-31 │ DJ     │ 02         │ missing      │ 30-34 │ 04              │ GS       │ 04    │
...
```
"""
function read end

function read(fullpath::Union{AbstractString,IO}, columnwidths::Union{Vector{UnitRange{Int}}, Vector{Int}}, sink::Type=DataFrame, args...; append::Bool=false, transforms::Dict=Dict{Int,Function}(), transpose::Bool=false, kwargs...)
    source =  Source(fullpath, columnwidths; kwargs...)
    sink = Data.stream!(source, sink, args...; append=append, transforms=transforms)
    return Data.close!(sink)
end

function read(fullpath::Union{AbstractString,IO}, columnwidths::Union{Vector{UnitRange{Int}}, Vector{Int}}, sink::T; append::Bool=false, transforms::Dict=Dict{Int,Function}(), transpose::Bool=false, kwargs...) where {T}
    source = Source(fullpath, columnwidths; kwargs...)
    sink = Data.stream!(source, sink; append=append, transforms=transforms)
    return Data.close!(sink)
end

#read(source::FWF.Source, columnwidths::Union{Vector{UnitRange{Int}}, sink=DataFrame, args...; append::Bool=false, transforms::Dict=Dict{Int,Function}()) = (sink = Data.stream!(source, sink, args...; append=append, transforms=transforms); return Data.close!(sink))
#read(source::FWF.Source, columnwidths::Union{Vector{UnitRange{Int}}, sink::T; append::Bool=false, transforms::Dict=Dict{Int,Function}()) where {T} = (sink = Data.stream!(source, sink; append=append, transforms=transforms); return Data.close!(sink))
#read(source::CSV.TransposedSource, sink=DataFrame, args...; append::Bool=false, transforms::Dict=Dict{Int,Function}()) = (sink = Data.stream!(source, sink, args...; append=append, transforms=transforms); return Data.close!(sink))
#read(source::CSV.TransposedSource, sink::T; append::Bool=false, transforms::Dict=Dict{Int,Function}()) where {T} = (sink = Data.stream!(source, sink; append=append, transforms=transforms); return Data.close!(sink))


text:
using SnpArrays, WiSER, TrajGWAS, Random, DataFrames, CSV, Statistics, Distributions, LinearAlgebra
using GeneticVariation, VCFTools, BGEN

# Simulate PLINK data
hapmap = SnpArray("hapmap3.bed")
famids = CSV.read("hapmap3.fam", header=false)[!, 1]
hapmap = convert(Matrix{Float64}, hapmap, model=ADDITIVE_MODEL, impute=true)
toplist = sortperm(vec(var(hapmap, dims=1))[1:12000], rev=true)[1:30]
# pick a few SNPs with not too much correlation
lastcor = 1
for i in toplist, j in toplist #find most independent SNPs
    curcor = abs(Statistics.cor(@view(hapmap[:, i]), @view(hapmap[:, j])))
    if curcor < lastcor
        lastcor = curcor
        println("$i and $j cor is $curcor")
    end
end

m, _ = size(hapmap)
ni = 10
Random.seed!(310)
sex = rand(0:1, m)
timeinvarX = [sex hapmap[:, 7828] hapmap[:, 8250] hapmap[:, 5] hapmap[:, 11793]]
p = 7
l = 7
q = 1

βtrue = [10.0; -3.2; 0.5; 0.5; 1.5; 0.25; 0.75]
τtrue = [0.2; -0.3; 0.5; -0.5; 0.5; 0.25; 0.15]
Σγ    = Matrix(Diagonal([2.0]))
δγω   = [0.0]
σω    = [0.0]

Σγω   = [Σγ δγω; δγω' σω]
Lγω   = cholesky(Symmetric(Σγω), check = false).L
Lγ    = Lγω[1:q, 1:q]
lγω   = Lγω[q + 1, 1:q]
lω    = Lγω[q + 1, q + 1]

γω = Vector{Float64}(undef, q + 1)
z  = similar(γω) # hold vector of iid std normal

df = DataFrame()
obsvec = Vector{WSVarLmmObs{Float64}}(undef, m)
for i in 1:m
    sex_i = timeinvarX[i, 1]
    snp1_i = timeinvarX[i, 2]
    snp2_i = timeinvarX[i, 3]
    snp3_i = timeinvarX[i, 4]
    snp4_i = timeinvarX[i, 5]
    # Xu from a normal dist 0 mean 2 std
    OnMeds = rand([0, 1], ni)
    # first column intercept, remaining entries iid std normal
    X = Matrix{Float64}(undef, ni, p)
    X[:, 1] .= 1
    X[:, 2] .= sex_i
    X[:, 3] = OnMeds
    X[:, 4] .= snp1_i
    X[:, 5] .= snp2_i
    X[:, 6] .= snp3_i
    X[:, 7] .= snp4_i

    # first column intercept, remaining entries iid std normal
    Z = Matrix{Float64}(undef, ni, q)
    Z[:, 1] .= 1
    # first column intercept, remaining entries iid std normal
    W = Matrix{Float64}(undef, ni, l)
    W[:, 1] .= 1
    W[:, 2] .= sex_i
    W[:, 3] = OnMeds
    W[:, 4] .= snp1_i
    W[:, 5] .= snp2_i
    W[:, 6] .= snp3_i
    W[:, 7] .= snp4_i
    # generate random effects: γω = Lγω * z
    mul!(γω, Lγω, randn!(z))
    # generate y
    μy = X * βtrue + Z * γω[1:q]
    @views ysd = exp.(0.5 .* (W * τtrue .+ dot(γω[1:q], lγω) .+ γω[end]))
    y = ysd .* randn(ni) .+ μy
    # form a VarLmmObs instance
    obsvec[i] = WSVarLmmObs(y, X, Z, W)
    tempdf = hcat(DataFrame(X[:, 2:end],
        [:sex, :onMeds, :snp1, :snp2, :snp3, :snp4]),
        DataFrame(y = y, id = fill(famids[i], ni))
    )
    df = vcat(df, tempdf)
end

fullmodel = WSVarLmmModel(@formula(y ~ 1 + sex + onMeds + snp1 + snp2 + snp3 + snp4),
              @formula(y ~ 1),
              @formula(y ~ 1 + sex + onMeds + snp1 + snp2 + snp3 + snp4),
              :id,
              df)
nullmodel = WSVarLmmModel(@formula(y ~ 1 + sex + onMeds),
                @formula(y ~ 1),
                @formula(y ~ 1 + sex + onMeds),
                :id,
                df)
WiSER.fit!(fullmodel)
WiSER.fit!(nullmodel)

CSV.write("trajgwas_plinkex.csv", df)

trajgwas(@formula(y ~ 1 + sex + onMeds),
@formula(y ~ 1),
@formula(y ~ 1 + sex + onMeds),
:id, "trajgwas_plinkex.csv", "hapmap3")

trajgwas(@formula(y ~ 1 + sex + onMeds),
       @formula(y ~ 1),
       @formula(y ~ 1 + sex + onMeds),
       :id,
       "trajgwas_plinkex.csv",
       "hapmap3",
       pvalfile = "snpset.window5.pval.txt",
       analysistype="snpset",
       snpset = 10)

# Simulate VCF data

D = convert_ds(Float64, "test_vcf.vcf.gz"; key="DS", impute=true, center=false, scale=false)

# pick a few SNPs with not too much correlation
using Statistics
toplist = sortperm(vec(var(D, dims=[1])), rev=true)[1:30]
lastcor = 1
for i in toplist, j in toplist #find most independent SNPs
    curcor = abs(Statistics.cor(D[:, i], D[:, j]))
    if curcor < lastcor
        lastcor = curcor
        println("$i and $j cor is $curcor")
    end
end

people, _ = VCFTools.nsamples("test_vcf.vcf.gz"), nrecords("test_vcf.vcf.gz")
reader = VCF.Reader(openvcf("test_vcf.vcf.gz", "r"))
famids = VCF.header(reader).sampleID
close(reader)
Random.seed!(310)
sex = rand(0:1, people)
m = people

ni = 10
Random.seed!(310)
sex = rand(0:1, m)
timeinvarX = [sex D[:, 86] D[:, 316] D[:, 656]]
p = 6
l = 6
q = 1

βtrue = [10.0; -3.2; 0.5; 0.5; 1.5; 0.25]
τtrue = [0.2; -0.3; 0.5; -0.5; 0.5; 0.25]
Σγ    = Matrix(Diagonal([2.0]))
δγω   = [0.0]
σω    = [0.0]

Σγω   = [Σγ δγω; δγω' σω]
Lγω   = cholesky(Symmetric(Σγω), check = false).L
Lγ    = Lγω[1:q, 1:q]
lγω   = Lγω[q + 1, 1:q]
lω    = Lγω[q + 1, q + 1]

γω = Vector{Float64}(undef, q + 1)
z  = similar(γω) # hold vector of iid std normal

df = DataFrame()
obsvec = Vector{WSVarLmmObs{Float64}}(undef, m)
for i in 1:m
    sex_i = timeinvarX[i, 1]
    snp1_i = timeinvarX[i, 2]
    snp2_i = timeinvarX[i, 3]
    snp3_i = timeinvarX[i, 4]
    # Xu from a normal dist 0 mean 2 std
    OnMeds = rand([0, 1], ni)
    # first column intercept, remaining entries iid std normal
    X = Matrix{Float64}(undef, ni, p)
    X[:, 1] .= 1
    X[:, 2] .= sex_i
    X[:, 3] = OnMeds
    X[:, 4] .= snp1_i
    X[:, 5] .= snp2_i
    X[:, 6] .= snp3_i

    # first column intercept, remaining entries iid std normal
    Z = Matrix{Float64}(undef, ni, q)
    Z[:, 1] .= 1
    # first column intercept, remaining entries iid std normal
    W = Matrix{Float64}(undef, ni, l)
    W[:, 1] .= 1
    W[:, 2] .= sex_i
    W[:, 3] = OnMeds
    W[:, 4] .= snp1_i
    W[:, 5] .= snp2_i
    W[:, 6] .= snp3_i
    # generate random effects: γω = Lγω * z
    mul!(γω, Lγω, randn!(z))
    # generate y
    μy = X * βtrue + Z * γω[1:q]
    @views ysd = exp.(0.5 .* (W * τtrue .+ dot(γω[1:q], lγω) .+ γω[end]))
    y = ysd .* randn(ni) .+ μy
    # form a VarLmmObs instance
    obsvec[i] = WSVarLmmObs(y, X, Z, W)
    tempdf = hcat(DataFrame(X[:, 2:end],
        [:sex, :onMeds, :snp1, :snp2, :snp3]),
        DataFrame(y = y, id = fill(famids[i], ni))
    )
    df = vcat(df, tempdf)
end

fullmodel = WSVarLmmModel(@formula(y ~ 1 + sex + onMeds + snp1 + snp2 + snp3),
              @formula(y ~ 1),
              @formula(y ~ 1 + sex + onMeds + snp1 + snp2 + snp3),
              :id,
              df)
nullmodel = WSVarLmmModel(@formula(y ~ 1 + sex + onMeds),
                @formula(y ~ 1),
                @formula(y ~ 1 + sex + onMeds),
                :id,
                df)
WiSER.fit!(fullmodel)
WiSER.fit!(nullmodel)

CSV.write("trajgwas_vcf_ex.csv", df)

trajgwas(@formula(y ~ 1 + sex + onMeds),
        @formula(y ~ 1),
        @formula(y ~ 1 + sex + onMeds),
        :id,
        "trajgwas_vcf_ex.csv",
        "test_vcf",
        geneticformat = "VCF",
        vcftype = :DS,
        pvalfile = "trajgwas_vcfex.pval.txt")

# Simulate BGEN data
b = Bgen(BGEN.datadir("example.8bits.bgen"))
people = n_samples(b)
sampleids = samples(b)
variants = parse_variants(b; from_bgen_start=false)
#sample 3 and 15 and 25 to be causal
v3 = minor_allele_dosage!(b, variants[3]; T=Float64, mean_impute=true)
v15 = minor_allele_dosage!(b, variants[15]; T=Float64, mean_impute=true)
v25 = minor_allele_dosage!(b, variants[25]; T=Float64, mean_impute=true)

ni = 10
m = people
Random.seed!(310)
sex = rand(0:1, m)
timeinvarX = [sex v3 v15 v25]
p = 6
l = 6
q = 1

βtrue = [10.0; -3.2; 0.5; 0.5; 1.5; 0.25]
τtrue = [0.2; -0.3; 0.5; -0.5; 0.5; 0.25]
Σγ    = Matrix(Diagonal([2.0]))
δγω   = [0.0]
σω    = [0.0]

Σγω   = [Σγ δγω; δγω' σω]
Lγω   = cholesky(Symmetric(Σγω), check = false).L
Lγ    = Lγω[1:q, 1:q]
lγω   = Lγω[q + 1, 1:q]
lω    = Lγω[q + 1, q + 1]

γω = Vector{Float64}(undef, q + 1)
z  = similar(γω) # hold vector of iid std normal

df = DataFrame()
obsvec = Vector{WSVarLmmObs{Float64}}(undef, m)
for i in 1:m
    sex_i = timeinvarX[i, 1]
    snp1_i = timeinvarX[i, 2]
    snp2_i = timeinvarX[i, 3]
    snp3_i = timeinvarX[i, 4]
    # Xu from a normal dist 0 mean 2 std
    OnMeds = rand([0, 1], ni)
    # first column intercept, remaining entries iid std normal
    X = Matrix{Float64}(undef, ni, p)
    X[:, 1] .= 1
    X[:, 2] .= sex_i
    X[:, 3] = OnMeds
    X[:, 4] .= snp1_i
    X[:, 5] .= snp2_i
    X[:, 6] .= snp3_i

    # first column intercept, remaining entries iid std normal
    Z = Matrix{Float64}(undef, ni, q)
    Z[:, 1] .= 1
    # first column intercept, remaining entries iid std normal
    W = Matrix{Float64}(undef, ni, l)
    W[:, 1] .= 1
    W[:, 2] .= sex_i
    W[:, 3] = OnMeds
    W[:, 4] .= snp1_i
    W[:, 5] .= snp2_i
    W[:, 6] .= snp3_i
    # generate random effects: γω = Lγω * z
    mul!(γω, Lγω, randn!(z))
    # generate y
    μy = X * βtrue + Z * γω[1:q]
    @views ysd = exp.(0.5 .* (W * τtrue .+ dot(γω[1:q], lγω) .+ γω[end]))
    y = ysd .* randn(ni) .+ μy
    # form a VarLmmObs instance
    obsvec[i] = WSVarLmmObs(y, X, Z, W)
    tempdf = hcat(DataFrame(X[:, 2:end],
        [:sex, :onMeds, :snp1, :snp2, :snp3]),
        DataFrame(y = y, id = fill(sampleids[i], ni))
    )
    df = vcat(df, tempdf)
end

fullmodel = WSVarLmmModel(@formula(y ~ 1 + sex + onMeds + snp1 + snp2 + snp3),
              @formula(y ~ 1),
              @formula(y ~ 1 + sex + onMeds + snp1 + snp2 + snp3),
              :id,
              df)
nullmodel = WSVarLmmModel(@formula(y ~ 1 + sex + onMeds),
                @formula(y ~ 1),
                @formula(y ~ 1 + sex + onMeds),
                :id,
                df)
WiSER.fit!(fullmodel)
WiSER.fit!(nullmodel)

CSV.write("trajgwas_bgen_ex.csv", df)

trajgwas(@formula(y ~ 1 + sex + onMeds),
        @formula(y ~ 1),
        @formula(y ~ 1 + sex + onMeds),
        :id,
        "trajgwas_bgen_ex.csv",
        "example.8bits",
        usespa = false,
        geneticformat = "BGEN",
        pvalfile = "trajgwas_bgen.pval.txt")

trajgwas_res = CSV.read("trajgwas_bgen.pval.txt", DataFrame)

trajgwas(@formula(y ~ 1 + sex + onMeds),
        @formula(y ~ 1),
        @formula(y ~ 1 + sex + onMeds),
        :id,
        "trajgwas_bgen_ex.csv",
        "example.8bits",
        usespa = true,
        geneticformat = "BGEN",
        pvalfile = "trajgwas_bgen_spa.pval.txt")

trajgwas_res_spa = CSV.read("trajgwas_bgen_spa.pval.txt", DataFrame)


text:
export play

using Colors

# coordinate transform for Makie.jl
transform(x::Int) = p -> CartesianIndex(p[2], x-p[1]+1)

using .Makie

function init_screen(w::Observable{<:AbstractGridWorld}; resolution=(1000,1000))
    scene = Scene(resolution = resolution, raw = true, camera = campixel!)

    area = scene.px_area
    grid_size = size(w[].world)[2:3]
    grid_inds = CartesianIndices(grid_size)
    tile_size = @lift((widths($area)[1] / size($w.world, 2), widths($area)[2] / size($w.world, 3)))
    T = transform(size(w[].world, 2))
    boxes(pos, s) = [FRect2D((T(p).I .- (1,1)) .* s, s) for p in pos]
    centers(pos, s) = [(T(p).I .- (0.5,0.5)) .* s for p in pos]

    # 1. paint background
    poly!(scene, area)

    # 2. paint each kind of object
    for o in get_object(w[])
        if o !== EMPTY
            scatter!(scene, @lift(centers(findall($w.world[o,:,:]), $tile_size)), color=get_color(o), marker=convert(Char, o), markersize=@lift(minimum($tile_size)))
        end
    end

    # 3. paint stroke
    poly!(scene, @lift(boxes(vec(grid_inds), $tile_size)), color=:transparent, strokecolor = :lightgray, strokewidth = 4)

    # 3. paint agent's view
    view_boxes = @lift boxes([p for p in get_agent_view_inds($w) if p ∈ grid_inds], $tile_size)
    poly!(scene, view_boxes, color="rgba(255,255,255,0.2)")

    # 4. paint agent
    agent = @lift(get_agent($w))
    agent_position = @lift((T(get_agent_pos($w)).I .- (0.5, 0.5)).* $tile_size)
    scatter!(scene, agent_position, color=@lift(get_color($agent)), marker=@lift(convert(Char, $agent)), markersize=@lift(minimum($tile_size)))

    display(scene)
    scene
end

function play(environment::AbstractGridWorld;file_name=nothing,frame_rate=24)
    print("""
    Key bindings:
    ←: TurnLeft
    →: TurnRight
    ↑: MoveForward
    q: Quit
    """)
    w = environment
    w_node = Node(w)
    scene = init_screen(w_node)
    is_quit = Ref(false)

    if !isnothing(file_name)
        vs = VideoStream(scene)
        recordframe!(vs)
    end

    on(scene.events.keyboardbuttons) do b
        if ispressed(b, Keyboard.left)
            w(TURN_LEFT)
            w_node[] = w
        elseif ispressed(b, Keyboard.right)
            w(TURN_RIGHT)
            w_node[] = w
        elseif ispressed(b, Keyboard.up)
            w(MOVE_FORWARD)
            w_node[] = w
        elseif ispressed(b, Keyboard.q)
            is_quit[] = true
        end
        isnothing(file_name) || recordframe!(vs)
    end

    try
        while !is_quit[]
            sleep(0.5)
        end
    catch
    end
    isnothing(file_name) || save(file_name, vs;framerate=frame_rate)
    Makie.GLMakie.destroy!(Makie.GLMakie.global_gl_screen())
end


text:
function rhs!(dX::Vector{T},X::Vector{T},α::T,F::T,s_inv::T) where {T<:AbstractFloat}

    n = length(dX)
    @boundscheck n == length(X) || throw(BoundsError())

    # periodic boundary conditions
    @inbounds dX[1] = (X[2]-X[n-1])*X[n]*s_inv - α*X[1] + F
    @inbounds dX[2] = (X[3]-X[n])*X[1]*s_inv - α*X[2] + F

    # in the middle of the domain
    for i in 3:n-1
        @inbounds dX[i] = (X[i+1] - X[i-2])*X[i-1]*s_inv - α*X[i] + F
    end

    # periodic boundary conditions
    @inbounds dX[n] = (X[1] - X[n-2])*X[n-1]*s_inv - α*X[n] + F
end


text:
using ApproxOperator, LinearAlgebra

aps = import_msh("./msh/cube_pfm.msh")
nₚ = length(aps["surface"][1].nodes)
set_integration_rule!(aps["surface"],Val(:PhaseFieldFriction))

E = 2.0e11
ν = 0.3
c = 40
𝜙 = 15/180*π
𝜙ᵣ = 15/180*π
θ = 0.
u₀ = 5.0e-3
nls = 50
k̄ = 2700
l = 0.01
tol = 1e-10
maxiter = 50

b = (x,y,z)->(0.0,0.0)
n = (x,y,z)->(1.0,0.0,1.0)
ū = (x,y,z)->(0.0,0.0,0.0)

ops = [
    NonlinearPlaneStress_C_Ω(b),
    NonlinearPlaneStress_Γᵍ_penalty(ū,n,E*1e7),
    SecondOrderPhaseField(k̄,l),
    Update_Friction_PhaseField_PlaneStress(E,ν,c,𝜙,𝜙ᵣ,1e-6),
    VTKExport(
                filename="figure/movie00.vtk",
                topic="phase field modeling fracture",
                parameters=Dict{String,Float64}("E"=>E,"ν"=>ν)
            )
]

d = zeros(2*nₚ)
k = zeros(2*nₚ,2*nₚ)
f = zeros(2*nₚ)
kᵥ = zeros(nₚ,nₚ)
fᵥ = zeros(nₚ)
dᵥ = zeros(nₚ)
dᵗ = zeros(nₚ)

ops[5](aps["surface"],d,dᵥ,Val(:PFM_PlaneStress))

for t in 1:nls
    # define essential BC
    u = (x,y,z)->(cos(θ)*u₀/nls*t,sin(θ)*u₀/nls*t)
    op_EBC = NonlinearPlaneStress_Γᵍ_penalty(u,n,E*1e7)
    iter = 0

    # Alternate minimization algorithm
    while iter ≤ maxiter
        dᵗ .= dᵥ
        iter += 1

        # update history field
        ops[4](aps["surface"],d,dᵥ)

        # loop of elasticity
        iteru = 0
        while iteru ≤ maxiter
            iteru += 1
            fill!(k,0.0)
            fill!(f,0.0)
            ops[1](aps["surface"],k,f)
            ops[2](aps["disp_fix"],k,f,d)
            op_EBC(aps["disp"],k,f,d)
            erru = norm(f)/E/1e7
            print("u: time step = $t, alternate iter = $iter, disp iter = $iteru, erru = $erru\n")
            if erru < tol;break;end
            d .+= k\f
        end

        # loop of phase field
        iterv = 0
        while iterv ≤ maxiter
            iterv += 1
            fill!(kᵥ,0.0)
            fill!(fᵥ,0.0)
            ops[3](aps["surface"],kᵥ,fᵥ,dᵥ)
            errv = norm(fᵥ)/k̄
            print("v: time step = $t, alternate iter = $iter, damage iter = $iterv, errv = $errv\n")
            if errv < tol;break;end
            dᵥ .+= kᵥ\fᵥ
        end

        err = norm(dᵗ .- dᵥ)/norm(dᵥ)
        print("err = $err\n")
        if err < tol;break;end
    end
    for ap in aps["surface"]
        for qw in ap.qw
            qw.ℋₙ = qw.ℋ
            qw.ℋₘ = max(qw.ℋₘ,qw.ℋₙ)
        end
    end
    # export to VTK file
    ops[5].filename = "figure/movie"*string(t,pad=2)*".vtk"
    ops[5](aps["surface"],d,dᵥ,Val(:PFM_PlaneStress))
end

max_stars_repo_name: jmwjc/ApproxOperator

max_issues_repo_name: jmwjc/ApproxOperator

max_forks_repo_name: jmwjc/ApproxOperator


text:
"""
atleast2d \\

For formatting the dimensions of the random variables to at least 2 dimensions  \\
Arguments:  random_variable \\
Returns:    A 2-d version of the variable \\
"""
function atleast2d(random_variable::Array)
  is_one_d = 1 == ndims(random_variable)
  return is_one_d ? reshape(random_variable, size(random_variable)[1], 1) : random_variable
end


text:
using Parameters, Test, Markdown, REPL

# misc
a8679 = @eval (a=1, b=2)
ra8679 = @eval (a=1, b=44)
@test ra8679 == reconstruct(a8679, b=44)
@test_throws ErrorException reconstruct(a8679, c=44)

a8679 = Dict(:a=>1, :b=>2)
@test Dict(:a=>1, :b=>44) == reconstruct(a8679, b=44)
@test a8679 == Dict(:a=>1, :b=>2)
@test_throws ErrorException reconstruct(a8679, c=44)

struct A8679
    a
    b
end
a8679 = A8679(1, 2)
@test A8679(1, 44) == reconstruct(a8679, b=44)
@test_throws ErrorException reconstruct(a8679, c=44)

##########
# @with_kw
##########

"Test documentation"
@with_kw struct MT1
    "Field r"
    r::Int = 4
    "A field"
    c = "sdaf"
end
@test MT1().r==4
@test MT1().c=="sdaf"
@test "Test documentation\n" == Markdown.plain(@doc MT1)
# https://github.com/JuliaLang/julia/issues/27092 means this does not work:
# @test "A field Default: sdaf\n" == Markdown.plain(@doc MT1.c)
@test "Field r Default: 4\n" == Markdown.plain(REPL.fielddoc(MT1, :r))
@test "A field Default: sdaf\n" == Markdown.plain(REPL.fielddoc(MT1, :c))

abstract type AMT1_2 end
"Test documentation with type-parameter"
@with_kw struct MT1_2{T} <: AMT1_2
    "Field r"
    r::Int = 4
    "A field"
    c = "sdaf"
end
@test MT1_2{Int}().r==4
@test MT1_2{Int}().c=="sdaf"
@test "Test documentation with type-parameter\n" == Markdown.plain(@doc MT1_2)
const TMT1_2 = MT1_2{Int} # Julia bug https://github.com/JuliaLang/julia/issues/27656
@test "Field r Default: 4\n" == Markdown.plain(REPL.fielddoc(TMT1_2, :r))
@test "A field Default: sdaf\n" == Markdown.plain(REPL.fielddoc(TMT1_2, :c))

"Test documentation with bound type-parameter"
@with_kw struct MT1_3{T} <: AMT1_2
    "Field r"
    r::Int = 4
    "A field"
    c::T = "sdaf"
end
@test MT1_3().r==4
@test MT1_3().c=="sdaf"
@test "Test documentation with bound type-parameter\n" == Markdown.plain(@doc MT1_3)
const TMT1_3 = MT1_3{Int} # Julia bug https://github.com/JuliaLang/julia/issues/27656
@test "Field r Default: 4\n" == Markdown.plain(REPL.fielddoc(TMT1_3, :r))
@test "A field Default: sdaf\n" == Markdown.plain(REPL.fielddoc(TMT1_3, :c))

# parameter-less
@with_kw_noshow mutable struct MT2
    r::Int
    c
    a::Float64
end
MT2(r=4, a=5., c=6)
MT2(r=4, a=5, c="asdf")
MT2(4, "dsaf", 5)
@test_throws ErrorException MT2(r=4)
@test_throws ErrorException MT2()

@with_kw mutable struct MT3
    r::Int=5
    a::Float64
end
MT3(r=4, a=5.)
MT3(r=4, a=5)
MT3(a=5)
MT3(4,5)
@test_throws ErrorException MT3(r=4)
@test_throws ErrorException MT3()

# parameter-less with supertype
abstract type MM1 end
abstract type MM2{T} end
@with_kw mutable struct MT3_1 <: MM1
    r::Int=5
    a::Float64
end
MT3_1(r=4, a=5.)
MT3_1(r=4, a=5)
MT3_1(a=5)
MT3_1(4,5)
@test_throws ErrorException MT3_1(r=4)
@test_throws ErrorException MT3_1()

@with_kw mutable struct MT3_2 <: MM2{Int}
    r::Int=5
    a::Float64
end
MT3_2(r=4, a=5.)
MT3_2(r=4, a=5)
MT3_2(a=5)
MT3_2(4,5)
@test_throws ErrorException MT3_2(r=4)
@test_throws ErrorException MT3_2()

# with type-parameters
@with_kw struct MT4{R,I}
    r::R=5
    a::I
end
@test MT4(r=4, a=5.0)==MT4{Int,Float64}(4,5.0)
MT4{Float32, Int}(r=4, a=5.)
MT4{Float32, Int}(a=5.)
MT4{Float32, Int}(5.4, 4)  # inner positional
mt4=MT4(5.4, 4) # outer positional
@test MT4(mt4)==mt4 # outer reconstruct
@test MT4(mt4; a=77)==MT4(5.4, 77)
@test_throws ErrorException MT4{Float32, Int}()
@test_throws InexactError MT4{Float32,Int}(a=5.5)
@test_throws InexactError MT4{Float32,Int}(5.5, 5.5)

# with type-parameters 2
abstract type AMT{R<:Real} end
@with_kw mutable struct MT5{R,I<:Integer} <: AMT{R}
    r::R=5
    a::I
end
@test_throws MethodError MT5(r=4, a=5.) # a has wrong type
MT5{Float32, Int}(r=4, a=5.)  # a gets converted
MT5{Float32, Int}(a=5.)
MT5{Float32, Int}(5.4, 4)  # inner positional
mt5=MT5(5.4, 4) # outer positional
@test MT5(mt5).r==mt5.r # outer reconstruct
@test MT5(mt5).a==mt5.a # outer reconstruct
@test MT5(mt5; a=77).a==MT5(5.4, 77).a
@test MT5(mt5; a=77).r==MT5(5.4, 77).r
@test_throws ErrorException MT5{Float32, Int}()
@test_throws InexactError MT5{Float32,Int}(a=5.5)
@test_throws InexactError MT5{Float32,Int}(5.5, 5.5)
@test_throws  MethodError MT5(5., "asdf")
@test_throws  TypeError MT5( "asdf", 5)
@test_throws  TypeError MT5{Float64, String}(5., "asdf")
@test_throws  TypeError MT5{String, Int}("asdf", 6)

# with type parameters and supertype
@with_kw mutable struct MT4_1{T} <: MM1
    r::Int=5
    a::T
end
MT4_1{Float64}(r=4, a=5.)
MT4_1{Float64}(r=4, a=5)
MT4_1{Float64}(a=5)
MT4_1{Float64}(4,5)
@test_throws ErrorException MT4_1{Float64}(r=4)
@test_throws ErrorException MT4_1{Float64}()

@with_kw mutable struct MT4_2{T} <: MM2{T}
    r::Int=5
    a::T
end
MT4_2{Float64}(r=4, a=5.)
MT4_2{Float64}(r=4, a=5)
MT4_2{Float64}(a=5)
MT4_2{Float64}(4,5)
@test_throws ErrorException MT4_2{Float64}(r=4)
@test_throws ErrorException MT4_2{Float64}()

# user defined inner positional constructor
@with_kw struct MT6{R,I<:Integer} <: AMT{R}
    r::R=5
    a::I
    MT6{R,I}(r) where {R,I} = new{R,I}(r,r)
    MT6{R,I}(r,a) where {R,I} = (@assert a>r; new{R,I}(r,a))
end
@test_throws MethodError MT6(r=4, a=5.) # need to specify type parameters
MT6{Float32, Int}(r=4, a=5.)
MT6{Float32, Int}(a=6.)
MT6{Float32, Int}(5.4, 6)  # inner positional
mt6=MT6(5.4, 6) # outer positional
@test MT6(mt6)==mt6 # outer reconstruct
@test MT6(mt6; a=77)==MT6(5.4, 77)
@test_throws ErrorException MT6{Float32, Int}()
@test_throws InexactError MT6{Float32,Int}(a=5.5)
@test_throws InexactError MT6{Float32,Int}(5.5, 6.5)
@test_throws  MethodError MT6(5., "asdf")
@test_throws  TypeError MT6( "asdf", 5)
@test_throws  TypeError MT6{Float64, String}(5., "asdf")
@test_throws  TypeError MT6{String, Int}("asdf", 6)

# user defined BAD inner positional constructor
@with_kw mutable struct MT7{R,I<:Integer} <: AMT{R}
    r::R=5
    a::I
    MT7{R,R}(r::R) where {R} = new{R,R}(r,r+8)
    # no MT7(r,a)
end
@test_throws MethodError MT7{Float32, Int}(r=4, a=5.)

# user defined BAD inner positional constructor
tmp = :(struct MT8{R,I<:Integer} <: AMT{R}
        r::R=5
        a::I
        MT8() = new{Int,Int}(5,6) # this would shadow the keyword constructor!
        MT8{R,I}(r,a) where {R,I} = new{R,I}(r,a)
        end)
@test_throws ErrorException Parameters.with_kw(tmp, @__MODULE__)
tmp = :(mutable struct MT8{R,I<:Integer} <: AMT{R}
        r::R=5
        a::I
        MT8(;a=7) = new{Int,Int}(5,a) # this would shadow the keyword constructor!
        MT8{R,I}(r,a) where {R,I} = new{R,I}(r,a)
        end)
@test_throws ErrorException Parameters.with_kw(tmp, @__MODULE__)

# default type annotation (adapted from MT6 test above)
@with_kw struct MT8{R,I<:Integer} <: AMT{R} @deftype R
    r=5
    a::I
    MT8{R,R}(r::R) where {R} = new{R,R}(r,r)
    MT8{R,I}(r,a) where {R,I} = (@assert a>r; new{R,I}(r,a))
end
@test_throws MethodError MT8(r=4, a=5.) # need to specify type parameters
MT8{Float32, Int}(r=4, a=5.)
MT8{Float32, Int}(a=6.)
MT8{Float32, Int}(5.4, 6)  # inner positional
mt6=MT8(5.4, 6) # outer positional
@test MT8(mt6)==mt6 # outer reconstruct
@test MT8(mt6; a=77)==MT8(5.4, 77)
@test_throws ErrorException MT8{Float32, Int}()
@test_throws InexactError MT8{Float32,Int}(a=5.5)
@test_throws InexactError MT8{Float32,Int}(5.5, 6.5)
@test_throws  MethodError MT8(5., "asdf")
@test_throws  TypeError MT8( "asdf", 5)
@test_throws  TypeError MT8{Float64, String}(5., "asdf")
@test_throws  TypeError MT8{String, Int}("asdf", 6)
@test MT8.var.name==:R
@test MT8.body.var.name==:I
@test MT8.body.var.ub==Integer
@test MT8{Float32,Int32}.types[1]==Float32
@test MT8{Float32,Int32}.types[2]==Int32

# parameter interdependence
@with_kw struct MT9{R<:Real}
    a::R = 5
    b::R
    c::R = a+b
end
@test MT9{Float64}(b=1).c==6
@test MT9{Float64}(b=1, c=1).c==1

@with_kw struct MT10{R<:Real}
    b::R = 6
    c::R = a+b
    a::R = 5
end
@test_throws UndefVarError MT10{Float64}() # defaults are evaluated in order
@test_throws UndefVarError MT10{Float64}(b=1).c
@test MT10{Float64}(b=1, c=1).c==1

# binding outside variables
a_ = 7
b_ = [1,2]
@with_kw struct MT11
    a::Int=a_  # a::Int=a is not possible as the outside a gets shadowed
    b::Vector{Int}=b_
end
m = MT11()
@test m.a===a_
@test m.b===b_
a_ = 5
b_[1] = 2
@test m.a===7
@test m.b[1]==2

## (Un)pack
@with_kw struct P1
    r::Int
    c
    a::Float64
end

let
    mt = P1(r=4, a=5, c=6)
    @unpack_P1 mt
    @test r===4
    @test c===6
    @test a===5.
    r = 1
    a = 2
    c = 3
    @test_throws LoadError eval(:(@pack!_P1 mt))
    mt2 = let P1 = error  # `P1` should not be resolved in this namespace
        @pack_P1
    end
    @test mt2 === P1(r=r, c=c, a=a)
end

@with_kw mutable struct P1m
    r::Int
    c
    a::Float64
end

let
    mt = P1m(r=4, a=5, c=6)
    @unpack_P1 mt
    @test r===4
    @test c===6
    @test a===5.
    r = 1
    a = 2
    c = 3
    mt2 = let P1m = error  # `P1m` should not be resolved in this namespace
        @pack_P1m
    end
    @test (mt2.r, mt2.c, mt2.a) == (r, c, a)
    @test (mt.r, mt.c, mt.a) === (4, 6, 5.0)  # mt is unmodified
    @pack_P1m! mt
    if Int==Int64
        @test string(mt) == "P1m\n  r: Int64 1\n  c: Int64 3\n  a: Float64 2.0\n"
    else
        @test string(mt) == "P1m\n  r: Int32 1\n  c: Int32 3\n  a: Float64 2.0\n"
    end
end

### Assertions
@with_kw struct MT12
    a=5; @assert a>=5
    b
    @assert b>a
end

@test_throws AssertionError MT12(b=2)
@test_throws AssertionError MT12(a=1,b=2)
@test MT12(b=6)==MT12(5,6)

# only asserts allowed if no inner constructors
@test_throws ErrorException Parameters.with_kw(:(struct MT13
                                                   a=5;
                                                   @assert a>=5
                                                   MT13(a) = new(8)
                                                 end),
                                               @__MODULE__)

# issue #29: assertions with parameterized types
@with_kw struct MT12a{R}
    a::Array{R,1}
    @assert 1 == length(a)
end
@test_throws AssertionError MT12a([1,2])
@test MT12a([1]).a==MT12a(a=[1]).a

####
# issue 10: infer type parameters from kw-args
@with_kw struct I10{T} @deftype Int
    a::T
    b = 10
    c::T="aaa"
end
@test_throws ErrorException I10()
@test I10(1,2,3)==I10{Int}(1,2,3)
@test_throws MethodError I10(a=10) # typeof(a)!=typeof(c)
a_ = I10(a="asd")
b_ = I10{String}("asd",10,"aaa")
@test a_==b_

@with_kw struct I10a{T}
    a::T
end
@test I10a([1]).a==I10a(a=[1]).a
@test I10a(1)==I10a(a=1)
@with_kw struct I10b{T}
    a::Vector{T}
end
@test I10b([1]).a==I10b(a=[1]).a
@test_throws MethodError I10b(a=1)

# issue #12: only one-liner inner constructors were parsed correctly.
@with_kw struct T9867
    r::Int
    function T9867(r)
        new(r)
    end
end
@test_throws ErrorException T9867()
@test T9867(r=2).r == 2

# issue #56: spliced-in modules
@eval @with_kw struct SplicedModule <: $Base.AbstractVector{Int}
    x
end
@test SplicedModule(12).x == 12

# NamedTuples
###
MyNT = @with_kw (a=1, b="test", w=:uu)
@test MyNT()==(a=1, b="test", w=:uu)
@test MyNT(b=1)==(a=1, b=1, w=:uu)
@test MyNT(1,2,"a")==(a=1, b=2, w="a")
@test_throws MethodError MyNT(c=1)
@test_throws LoadError @eval @with_kw (a::Int=1,) # no type annotations allowed
@test_throws LoadError @eval @with_kw(a=1,) # no space

MyNT2 = @with_kw (a=1, b="test", w)
@test_throws ErrorException MyNT2() # no default for w

MyNT3 = @with_kw (a=1, b=a)
@test MyNT3()==(a=1, b=1)
@test MyNT3(a=2)==(a=2, b=2)
@test MyNT3(b=2)==(a=1, b=2)

@test_throws ErrorException MyNT().z # Undefined field access
q = x -> x^3
scopingTest = @with_kw (q = q,)
@test scopingTest().q(2) == 8
x = [1, 2, 3]
scopingTest = @with_kw (x = x,)
@test scopingTest().x == [1, 2, 3]
x = [1, 2, 3]
immutabilityTest = @with_kw (f = x,)
obj = immutabilityTest()
x = [4, 5, 6]
@test obj.f == [1, 2, 3] # Test immutability
@test immutabilityTest().f == [4, 5, 6] # Test rebinding for new object

# Exotic input test
naughtyInputs = ["!@#\$%^&*()`~", :"()esc(;", :x, :(eval(:x)), true, esc(true), "-9223372036854775808/-1", :(0/0), "-9223372036854775808/-1", "0xabad1dea", :(@test)]
for input in naughtyInputs
    naughtyTest = @with_kw (x = input, y = 2)
    @test naughtyTest().x == input
end

# Other behavior checks
compoundTest = @with_kw (x = (2 > 1), y = 2)
@test compoundTest().x == true # Compound evaluation check
foo(x) = 3*x
anonymousTest = @with_kw (x = x -> foo(x), y = 2)
@test anonymousTest().x(4) == 12 # Anonymous functions/value clash check
symbolTest = @with_kw (x = :x, y = 2)
@test symbolTest().x == :x

###########################
# Packing and unpacking @unpack, @pack!
##########################
# Example with dict:
d = Dict{Symbol,Any}(:a=>5.0,:b=>2,:c=>"Hi!")
@unpack a, c = d
@test a == 5.0 #true
@test c == "Hi!" #true

d = Dict("a"=>5.0,"b"=>2,"c"=>"Hi!")
@unpack a, c = d
@test a == 5.0 #true
@test c == "Hi!" #true

# Example with named tuple
@eval d = (a=5.0, b=2, c="Hi!")
@unpack a, c = d
@test a == 5.0 #true
@test c == "Hi!" #true

mutable struct PropertyExample
    a
    last_set_property
end
Base.getproperty(::PropertyExample, name::Symbol) = String(name)
Base.setproperty!(d::PropertyExample, name::Symbol, value) =
    setfield!(d, :last_set_property, (name, value))
Base.propertynames(::PropertyExample) = (:A, :B, :C)

let d = PropertyExample(:should_be_ignored, nothing)
    @unpack a, b = d
    @test a == "a"
    @test b == "b"

    a = "a value"
    @pack! d = a
    @test getfield(d, :last_set_property) == (:a, "a value")

    @test type2dict(d) == Dict(:A => "A", :B => "B", :C => "C")
end

# TODO add test with non String string

# Example with type:
mutable struct A; a; b; c; end
d = A(4,7.0,"Hi!")
@unpack a, c = d
@test a == 4 #true
@test c == "Hi!" #true

## Packing

# Example with dict:
a = 5.0
c = "Hi!"
d = Dict{Symbol,Any}()
@pack! d = a, c
@test d==Dict{Symbol,Any}(:a=>5.0,:c=>"Hi!")

d = Dict{String,Any}()
@pack! d = a, c
@test d==Dict{String,Any}("a"=>5.0,"c"=>"Hi!")

# Example with type:
a = 99
c = "HaHa"
d = A(4,7.0,"Hi")
@pack! d = a, c
@test d.a == 99
@test d.c == "HaHa"

# older tests ported
mutable struct UP1
    aUP1
    bUP1
end
uu = UP1(1,2)
@test_throws ErrorException @unpack cUP1 = uu
@test_throws ErrorException @unpack aUP1, cUP1 = uu

aUP1, bUP1 = 0, 0
@unpack aUP1 = uu
@test aUP1==1
@test bUP1==0
aUP1, bUP1 = 0, 0
@unpack aUP1, bUP1 = uu
@test aUP1==1
@test bUP1==2

vv = uu
aUP1 = 99
@pack! uu = aUP1
@test uu==vv
@test uu.aUP1==99

struct UP2
    aUP2
    bUP2
end
uu = UP2(1,2)
@test_throws ErrorException @unpack cUP2 = uu
@test_throws ErrorException @unpack aUP2, cUP2 = uu

aUP1 = 99
@test_throws ErrorException @pack! uu = aUP1

# check that inference works
struct UP3
    a::Float64
    b::Int
end
function f(u::UP3)
    @unpack a,b = u
    a,b
end
@inferred f(UP3(1,2))

#
@with_kw struct UP4{T}
    g::T=9
    a::Float64=4
end
@test typeof(UP4())==UP4{Int}
@test UP4().g===9
@test UP4().a===4.0

# Issue 21
# A macro to create the same fields in several types:
macro def(name, definition)
    return quote
        macro $(esc(name))()
            esc($(Expr(:quote, definition)))
        end
    end
end
@def sharedparams1 begin
    a::Float64  = 1.0
    b::Int = 1
end

@def sharedparams2 begin
    c::Float64  = 1.0
    d::Int = 1
end

@with_kw struct MyType1
    @sharedparams1
    @sharedparams2
end

@test MyType1()==MyType1(1,1,1,1)

@def sharedparams3 begin
    e::Float64  = 1.0
    @assert x>0  # nested macros which shouldn't be expanded are not allowed
end

@test_throws AssertionError Parameters.with_kw(:(struct MyType2
    @sharedparams1
    @sharedparams2
    @sharedparams3
end), @__MODULE__)

@def sharedparams4 begin
    e::Float64  = 1.0
    @sharedparams1 # not allowed as it will lead to nested begin-end block
end
@test_throws ErrorException Parameters.with_kw(:(struct MyType2
    @sharedparams4
end), @__MODULE__)

### New 0.6 type system
@with_kw struct V06{T} @deftype Array{I,1} where I<:Integer
    a::T
    b = [10]
    c::Vector{S} where S<:AbstractString=["aaa"]
end
V06(a=88)
V06(a=88, b=[1], c=["a"])
V06(88, [1], ["a"])

### test escaping
module TestModule

using Parameters: @unpack, @pack!, @with_kw

@with_kw mutable struct TestStruct
    x::Int = 1
    y::Float64 = 42.0
end

function test_function(z::TestStruct)
    @unpack x, y = z
    y += x
    @pack! z = y
    z
end

end

z2 = TestModule.test_function(TestModule.TestStruct(; y = 9.0))
@test z2.x == 1 && z2.y == 10.0

# issue 91
@testset "i91" begin
    function foo()
        a = @with_kw (x = 1, y = 2)
        return a
    end
    @test foo()() == (x = 1, y = 2)
    @test foo()(y=4) == (x = 1, y = 4)
end

####
# @consts
####
@consts begin
    a34 = 1
    b34 = 2.0
    d34 = sin(56)
end
@test a34 == 1
@test b34 == 2.0
@test_throws ErrorException global b34 = 1
## this test is dependent on Julia version, leave it:
# @test_warn "WARNING: redefinition of constant b34. This may fail, cause incorrect answers, or produce other errors." global b34 = 4.0