text:
# Atmospheric model which sits on top of the interior structure model

import ForwardDiff: derivative
import Optim: optimize, OnceDifferentiable, Fminbox, minimizer

# New types

"This system includes both temperature and optical depth details"
struct WithAtmosphere <: ModelComplexity end

"An opacity in units of m^2 / kg"
abstract type MassOpacity end

"A power law opacity of the form κ = C P^α T^β (in SI units)"
struct PowerLawOpacity{Op<:Opacity} <: MassOpacity
    C::Op
    α::Float64
    β::Float64
end
ConstantOpacity(C) = PowerLawOpacity(C, 0.0, 0.0)

abstract type OpacityRatio end

struct ConstantOpacityRatio <: OpacityRatio
    γ::Float64
end
(c::ConstantOpacityRatio)(P, T) = c.γ

struct VariableOpacityRatio <: OpacityRatio
    γ::Function
end
(v::VariableOpacityRatio)(P, T) = v.γ(P, T)

"The optical depth gradient w.r.t. mass, dτ/dm = -κ/4πr²"
struct OpticalDepthGradient <: StructureEquation
    κ::PowerLawOpacity
    is_radiative::Function
end

"Temperature gradient for a radiative two-stream atmosphere"
mutable struct TwoStreamTemperatureGradient{T1<:Temperature, T2<:Temperature} <: StructureEquation
    κ::PowerLawOpacity
    Tint::T1
    Tirr::T2
    γ::OpacityRatio
end

"Temperature gradient for a combined atmosphere and interior model"
struct CombinedTemperatureGradient <: StructureEquation
    atmosphere_gradient::TwoStreamTemperatureGradient
    interior_gradient::TemperatureGradient
    is_radiative::Function
end

"A planet with temperature dependence and an atmosphere"
mutable struct AtmospherePlanet <: PlanetSystem{WithAtmosphere}
    M::Float64
    structure_equations::EquationSet
    boundary_values::BoundaryValues{WithAtmosphere}
    solution_grid::Vector{Float64}
    radius_search_bracket::Vector{Float64}
    refine_surface!::Nullable{Function}
end

"Planetary structure that holds temperature and opacity info"
mutable struct AtmospherePlanetStructure <: PlanetStructure{WithAtmosphere}
    data::Matrix{Float64}

    function AtmospherePlanetStructure(data::Matrix{Float64})
        if size(data)[1] == 5
            new(data)
        else
            error("Expected a size (5,n) array: got $(size(data))")
        end
    end
end
PlanetStructure(m, r, P, T, τ) = AtmospherePlanetStructure(hcat(m, r, P, T, τ)')

"Holds physical values of mass, radius, pressure, temperature, opacity"
mutable struct AtmosphereValues{Ma<:Mass,Ra<:Radius,
                      Pr<:Pressure,Te<:Temperature,
                      Op<:Opacity} <: ValueSet{WithAtmosphere}
    m::Ma
    r::Ra
    P::Pr
    T::Te
    τ::Op
end
ValueSet(m, r, P, T, τ) = AtmosphereValues(m, r, P, T, τ)

# Temperature profiles

"Angle-dependent atmospheric temperature profile"
function temperature_profile(τ::OpticalDepth,
                             T_int::Temperature, T_irr::Temperature,
                             γ::Dimensionless, μ::Dimensionless=μ_isotropic)
    internal_term = 3/4 * T_int^4 * (2/3 + τ)
    irradiation_term = 3/4 * T_irr^4 * μ * (2/3 + μ/γ + (γ/3μ - μ/γ)*exp(-γ*τ/μ))

    T⁴ = internal_term + irradiation_term
    return (T⁴ ^ (1/4))::Temperature
end
"Derivative of the angle-dependent atmospheric temperature profile"
function temperature_profile_deriv(τ::OpticalDepth,
                                   T_int::Temperature, T_irr::Temperature,
                                   γ::Dimensionless, μ::Dimensionless=μ_isotropic)
    T = τ -> temperature_profile(τ, T_int, T_irr, γ, μ)
    dTdτ = derivative(T, τ)
end

# Evaluating structure

function (κ::PowerLawOpacity)(P::Pressure, T::Temperature)
    C = κ.C
    α = κ.α
    β = κ.β

    return (C * (P/1Pa)^α * (T/1K)^β)::Opacity
end

"The optical-depth--radius relation for a radiative layer"
function dτdm(κ::Dimensionless, r::Distance)
    return -κ/(4π*r^2)
end

function (odg::OpticalDepthGradient)(vs::ValueSet)
    if isphysical(vs) && odg.is_radiative(vs)
        r = radius(vs)
        P = pressure(vs)
        T = temperature(vs)
        κ = odg.κ(P, T)

        return dτdm(κ, r)
    else
        return zero(Float64)
    end
end

function (tst::TwoStreamTemperatureGradient)(vs::ValueSet)
    dTdm = zero(Float64)

    if isphysical(vs)
        r = radius(vs)
        P = pressure(vs)
        T = temperature(vs)
        τ = opticaldepth(vs)
        κ = tst.κ(P, T)
        Tint = tst.Tint
        Tirr = tst.Tirr
        γ = tst.γ(P, T)
        dTdτ = temperature_profile_deriv(τ, Tint, Tirr, γ)
        dTdm = dTdτ * dτdm(κ, r)
    end

    return dTdm
end

function (ctg::CombinedTemperatureGradient)(vs::ValueSet)
    dTdm = zero(Float64)
    atmosphere = ctg.atmosphere_gradient
    interior = ctg.interior_gradient
    is_radiative = ctg.is_radiative

    dTdm = is_radiative(vs) ? atmosphere(vs) : interior(vs)

    return dTdm
end

# FIXME: this is a workaround for julia issue #14919
# (the macro approach ran into segfault issues)
for structEq in (MassContinuity, PressureBalance, TemperatureGradient,
                 CombinedTemperatureGradient, TwoStreamTemperatureGradient,
                 OpticalDepthGradient)
    @eval begin
        let MRP = Ogre.MassRadiusPressure, PV = Ogre.PhysicalValues, AV=Ogre.AtmosphereValues
            (eq::$structEq)(m, y::Vector) = eq(m, y...)
            (eq::$structEq)(m, r, P) = eq(MRP(m, r, P))
            (eq::$structEq)(m, r, P, T) = eq(PV(m, r, P, T))
            (eq::$structEq)(m, r, P, T, τ) = eq(AV(m, r, P, T, τ))
        end
    end
end

function (tg::TemperatureGradient)(av::AtmosphereValues)
    m = mass(av)
    r = radius(av)
    P = pressure(av)
    T = temperature(av)

    tg(m, r, P, T)
end
(eos::MassPiecewiseEOS)(vs::AtmosphereValues) = extracteos(eos, vs)(vs)

# Other features of the atmospheric model

"""Make a power law opacity from the functions in Kurosaki et al.
These functions have the form κ = D * (P / 1 bar)^α * (T / 1000 K)^β cm^2/g.
This function returns a power law opacity κ(P, T), accepting P/Pa and T/K."""
function kurosaki_opacity(D, α, β)
    C = D * cm^2/g * (1/1bar)^α * (1/1000K)^β
    PowerLawOpacity(C, α, β)
end

const water_opacity_r_vis = kurosaki_opacity(2.20, 1.0, -0.4)
const water_opacity_r_th = kurosaki_opacity(3.07e2, 0.9, -4.0)
const water_opacity_p_vis = kurosaki_opacity(1.94e4, 0.01, 1.0)
const water_opacity_p_th = kurosaki_opacity(4.15e5, 0.01, -1.1)

const water_opacity = water_opacity_r_th
const α = water_opacity_r_th.α
const β = water_opacity_r_th.β
const C = water_opacity_r_th.C * m^2 / kg
const water_opacity_ratio = VariableOpacityRatio((P,T) -> water_opacity_r_vis(P, T) / water_opacity_r_th(P, T))

is_radiative_default(vs::ValueSet) = (pressure(vs) < P_rad_max)

function scale_height(T, g)
    H = (R_H2O * T / g)
    return H
end

function surface_opticaldepth(M, R, T, γ)
    g = surface_gravity(M, R)
    H = scale_height(T, g)
    τ = (1/γ) * sqrt(H / (2π*(α+1)*R))
    return τ
end

function surface_pressure(M, R, T, τ)
    # this function is written in dimensionless form so we divide out units
    let G = G/(m^3/kg/s^2), M = M/kg, R = R/m, C = C/(m^2/kg), T = T/K
        num = G * M * (α+1) * τ
        den = R^2 * C * T^β
        P = (num/den)^(1/(α+1))
        return P * Pa
    end
end

function surface_temperature(τsurf, Tint, Tirr, γ)
    temperature_profile(τsurf, Tint/K, Tirr/K, γ) * K
end

function surface_T_and_τ(Tint, Tirr, M, R, γ)
    τguess = 0.01
    Tguess = Tirr
    τf(T) = surface_opticaldepth(M, R, T, γ)
    Tf(τ) = surface_temperature(τ, Tint, Tirr, γ)
    lsq(T, τ) = (Tf(τ)/K - T)^2 + (τf(T*K) - τ)^2
    lsq(Tτ) = lsq(Tτ...)
    lbounds = [0., 0.]
    ubounds = [Inf, Inf]
    result = optimize(OnceDifferentiable(lsq), [Tguess/K, τguess],
                      lbounds, ubounds, Fminbox())
    T = minimizer(result)[1] * K
    τ = minimizer(result)[2]
    return (T, τ)
end

function T_at_P(ps::PlanetStructure{WithAtmosphere}, P)
    Ps = pressure(ps)
    i = findfirst(x -> x >= P, Ps)
    return temperature(ps)[i]
end

# Other additions to Ogre methods

function blank_structure(sys::PlanetSystem{WithAtmosphere})
    n = npoints(sys)
    m = nvars(WithAtmosphere)
    AtmospherePlanetStructure(fill(NaN, m, n))
end

atm_gradient(sys::AtmospherePlanet) = sys.structure_equations[3].atmosphere_gradient

function refine_surface!(system::AtmospherePlanet)
    if !isnull(system.refine_surface!)
        do_refine! = get(system.refine_surface!)
        do_refine!(system)
    end
end

nvars(::Type{WithAtmosphere}) = 5 # M, R, P, T, τ
opticaldepth(av::AtmosphereValues) = av.τ
temperature(av::AtmosphereValues) = av.T
opticaldepth(ps::PlanetStructure{WithAtmosphere}) = ps.data[5, :]
temperature(ps::PlanetStructure{WithAtmosphere}) = ps.data[4, :]
nonmass(ps::PlanetStructure{WithAtmosphere}) = ps.data[2:5, :]
nonmass(av::AtmosphereValues) = [radius(av), pressure(av), temperature(av),
                                 opticaldepth(av)]
setnonmass!(ps::PlanetStructure{WithAtmosphere}, i, rPTτ) = (ps.data[2:5, i] = rPTτ; nothing)
function density(ps::PlanetStructure{WithAtmosphere}, T, sys::PlanetSystem)
    eos = sys.structure_equations[1].eos
    map(eos, pressure(ps), temperature(ps))
end
centre(ps::PlanetStructure{WithAtmosphere}) = AtmosphereValues(ps.data[:, end]...)
surface(ps::PlanetStructure{WithAtmosphere}) = AtmosphereValues(ps.data[:, 1]...)

function Base.show(io::IO, atm::AtmosphereValues)
    m, r, P, T, τ = atm.m, atm.r, atm.P, atm.T, atm.τ
    show(io, "$(m/M_earth) M⊕, $(r/R_earth) R⊕, $(P/Pa) Pa, $(T/K) K, τ=$τ")
end

function isphysical(vs::AtmosphereValues)
    (mass(vs) > 0 && radius(vs) > 0 && pressure(vs) > 0 &&
     temperature(vs) > 0 && opticaldepth(vs) > 0)
end

# FIXME: this is a workaround for julia issue #14919
macro addEOSCall(eos)
    esc(
        quote
            let P = Ogre.pressure, T = Ogre.temperature,
                MRP = Ogre.MassRadiusPressure, PV = Ogre.PhysicalValues,
                AV = Ogre.AtmosphereValues
                (e::$eos)(vs::MRP) = e(P(vs))
                (e::$eos)(vs::PV) = e(P(vs), T(vs))
                (e::$eos)(vs::AV) = e(P(vs), T(vs))
            end
        end
    )
end

for eos in (TFD, BME3, BME4, Vinet, PolytropicEOS, WaterData.OutOfDomainEOS, BoundedEOS, IAPWS, MGDPressureEOS, PressurePiecewiseEOS, LineEOS, GridEOS, StitchedEOS)
    @addEOSCall eos
end

"Do-everything constructor for a watery planet with a watery atmosphere"
function AtmospherePlanet(M::Mass, eos::EOS, Cₚ::HeatCapacity,
                          bvs::BoundaryValues{WithAtmosphere},
                          Tint_initial::Temperature, Tirr::Temperature,
                          grid=linspace(M, 0, defaults.total_points),
                          r_bracket=defaults.R_bracket, κ_const=nothing,
                          γ_const=nothing, Psurf=P_rad_max, refine_surface=nothing)

    # set up structural equations
    masscontinuity = MassContinuity(eos)
    adiabatic_gradient = TemperatureGradient(eos, Cₚ)

    is_variable_opacity = κ_const == nothing
    is_variable_opacity_ratio = γ_const == nothing
    κ = (is_variable_opacity ? water_opacity : ConstantOpacity(κ_const))
    γ = (is_variable_opacity_ratio ? water_opacity_ratio : ConstantOpacityRatio(γ_const))
    changed_Pradmax = Psurf != P_rad_max
    is_radiative = if changed_Pradmax
        vs -> pressure(vs) < Psurf
    else
        is_radiative_default
    end
    wateratm_gradient = TwoStreamTemperatureGradient(κ, Tint_initial, Tirr, γ)
    opticaldepth_gradient = OpticalDepthGradient(κ, is_radiative)
    combined_temperature = CombinedTemperatureGradient(wateratm_gradient,
                                                       adiabatic_gradient,
                                                       is_radiative)

    structure = EquationSet([masscontinuity, pressurebalance,
                             combined_temperature, opticaldepth_gradient])

    AtmospherePlanet(M, structure, bvs, grid, r_bracket, refine_surface)
end

function mass_grid(M, f, Npoints)
    minusM = [0.]
    append!(minusM, logspace(-12, log10(f), Npoints÷2))
    append!(minusM, linspace(f, 1, Npoints÷2))
    grid = M * (1 - minusM)
end

function mass_grid_linear(M, f, Npoints)
    grid = linspace(M, 0, Npoints)
end


text:
"""
Test of the nonlinear FEM machine
"""
module m1testop1
using FinEtools
using FinEtoolsDeforLinear.DeforModelRedModule: DeforModelRed3D
using FinEtoolsDeforNonlinear
using FinEtoolsDeforNonlinear.MatDeforNeohookeanModule: MatDeforNeohookean
using FinEtoolsDeforNonlinear.FEMMDeforNonlinearModule: FEMMDeforNonlinear
using FinEtoolsDeforNonlinear.FEMMDeforNonlinearBaseModule: stiffness
using LinearAlgebra: norm
using SparseArrays
using DelimitedFiles
using Test
function test()
    mr = DeforModelRed3D
    E, nu = 7.0*phun("MPa"), 0.3
    m = MatDeforNeohookean(mr, E, nu)
    L= 6/2*phun("mm");
    H = 2/2*phun("mm");
    W = 2/2*phun("mm");
    fens, fes = H8block(L, W, H, 1, 1, 1)
    femm = FEMMDeforNonlinear(mr, IntegDomain(fes, GaussRule(3, 2)), m)

    geom = NodalField(fens.xyz)
    un1 = NodalField(zeros(size(fens.xyz, 1), 3)) # displacement field
    un = NodalField(zeros(size(fens.xyz, 1), 3)) # displacement field
    numberdofs!(un1)
    tn, dtn = 0.0, 1.0
    K = stiffness(femm, geom, un1, un, tn, dtn)
    ijs = readdlm("m1testop1.csv", ',', Float64)
    matK  = sparse(Int.(ijs[:, 1]), Int.(ijs[:, 2]), ijs[:, 3], size(K, 1), size(K, 2))
    # @show norm(K - matK) / maximum(matK[:])
    @test norm(K - matK) / maximum(matK[:]) < 1.0e-4
end
end
using .m1testop1
m1testop1.test()

module m1testop2
using FinEtools
using FinEtoolsDeforLinear.DeforModelRedModule: DeforModelRed3D
using FinEtoolsDeforNonlinear
using FinEtoolsDeforNonlinear.MatDeforNeohookeanModule: MatDeforNeohookean
using FinEtoolsDeforNonlinear.FEMMDeforNonlinearModule: FEMMDeforNonlinear
using FinEtoolsDeforNonlinear.FEMMDeforNonlinearBaseModule: stiffness
using LinearAlgebra: norm
using SparseArrays
using DelimitedFiles
using Test
function test()
    mr = DeforModelRed3D
    E, nu = 7.0*phun("MPa"), 0.3
    m = MatDeforNeohookean(mr, E, nu)
    L= 6/2*phun("mm");
    H = 2/2*phun("mm");
    W = 2/2*phun("mm");
    fens, fes = H8hexahedron([0 0 0; 1.2 0 -0.1; 1.0 0.9 0.1; -0.05 1.1 0; 0 0 1.03; 1.2 0 1.1; 1.0 0.95 0.81; -0.05 1.01 0.9], 1, 1, 1)
    femm = FEMMDeforNonlinear(mr, IntegDomain(fes, GaussRule(3, 2)), m)

    geom = NodalField(fens.xyz)
    un1 = NodalField(zeros(size(fens.xyz, 1), 3)) # displacement field
    un = NodalField(zeros(size(fens.xyz, 1), 3)) # displacement field
    numberdofs!(un1)
    tn, dtn = 0.0, 1.0
    K = stiffness(femm, geom, un1, un, tn, dtn)
    ijs = readdlm("m1testop2.csv", ',', Float64)
    matK  = sparse(Int.(ijs[:, 1]), Int.(ijs[:, 2]), ijs[:, 3], size(K, 1), size(K, 2))
    # @show norm(K - matK) / maximum(matK[:])
    @test norm(K - matK) / maximum(matK[:]) < 1.0e-4
end
end
using .m1testop2
m1testop2.test()

module m1testop3
using FinEtools
using FinEtoolsDeforLinear.DeforModelRedModule: DeforModelRed3D
using FinEtoolsDeforNonlinear
using FinEtoolsDeforNonlinear.MatDeforNeohookeanModule: MatDeforNeohookean
using FinEtoolsDeforNonlinear.FEMMDeforNonlinearModule: FEMMDeforNonlinear
using FinEtoolsDeforNonlinear.FEMMDeforNonlinearBaseModule: stiffness
using LinearAlgebra: norm
using SparseArrays
using DelimitedFiles
using Test
function test()
    mr = DeforModelRed3D
    E, nu = 7.0*phun("MPa"), 0.3
    m = MatDeforNeohookean(mr, E, nu)
    L= 6/2*phun("mm");
    H = 2/2*phun("mm");
    W = 2/2*phun("mm");
    fens, fes = H8hexahedron([0 0 0; 1.2 0 -0.1; 1.0 0.9 0.1; -0.05 1.1 0; 0 0 1.03; 1.2 0 1.1; 1.0 0.95 0.81; -0.05 1.01 0.9], 1, 1, 1)
    Rmout = fill(0.0, 3, 3)
    rv = vec([1 1.1 0.3])
    rotmat3!(Rmout, rv)
    mcsys = CSys(Rmout)
    femm = FEMMDeforNonlinear(mr, IntegDomain(fes, GaussRule(3, 2)), mcsys, m)

    geom = NodalField(fens.xyz)
    un1 = NodalField(zeros(size(fens.xyz, 1), 3)) # displacement field
    un = NodalField(zeros(size(fens.xyz, 1), 3)) # displacement field
    numberdofs!(un1)
    tn, dtn = 0.0, 1.0
    K = stiffness(femm, geom, un1, un, tn, dtn)
    ijs = readdlm("m1testop3.csv", ',', Float64)
    matK  = sparse(Int.(ijs[:, 1]), Int.(ijs[:, 2]), ijs[:, 3], size(K, 1), size(K, 2))
    # @show norm(K - matK) / maximum(matK[:])
    @test norm(K - matK) / maximum(matK[:]) < 1.0e-4
end
end
using .m1testop3
m1testop3.test()

module m1testop4
using FinEtools
using FinEtoolsDeforLinear.DeforModelRedModule: DeforModelRed3D
using FinEtoolsDeforNonlinear
using FinEtoolsDeforNonlinear.MatDeforNeohookeanModule: MatDeforNeohookean
using FinEtoolsDeforNonlinear.FEMMDeforNonlinearModule: FEMMDeforNonlinear
using FinEtoolsDeforNonlinear.FEMMDeforNonlinearBaseModule: stiffness
using LinearAlgebra: norm, cross
using SparseArrays
using DelimitedFiles
# using Debugger
using Test
function test()
    mr = DeforModelRed3D
    E, nu = 7.0*phun("MPa"), 0.3
    m = MatDeforNeohookean(mr, E, nu)
    L= 6/2*phun("mm");
    H = 2/2*phun("mm");
    W = 2/2*phun("mm");
    fens, fes = H8hexahedron([0 0 0; 1.2 0 -0.1; 1.0 0.9 0.1; -0.05 1.1 0; 0 0 1.03; 1.2 0 1.1; 1.0 0.95 0.81; -0.05 1.01 0.9], 1, 1, 1)
    function  mcsys(csmatout::FFltMat, XYZ, tangents, fe_label)
        # @bp
        n1 = vec(XYZ); n1[1] = 0; n1  = n1 / norm(n1);
        n3 = vec([1, 0, 0]);
        n2 = cross(n3, n1);
        csmatout .= hcat(n1, n2, n3);
    end
    femm = FEMMDeforNonlinear(mr, IntegDomain(fes, GaussRule(3, 2)), CSys(3, 3, mcsys), m)

    geom = NodalField(fens.xyz)
    un1 = NodalField(zeros(size(fens.xyz, 1), 3)) # displacement field
    un = NodalField(zeros(size(fens.xyz, 1), 3)) # displacement field
    numberdofs!(un1)
    tn, dtn = 0.0, 1.0
    K = stiffness(femm, geom, un1, un, tn, dtn)
    ijs = readdlm("m1testop4.csv", ',', Float64)
    matK  = sparse(Int.(ijs[:, 1]), Int.(ijs[:, 2]), ijs[:, 3], size(K, 1), size(K, 2))
    # @show norm(K - matK) / maximum(matK[:])
    @test norm(K - matK) / maximum(matK[:]) < 1.0e-4
end
end
using .m1testop4
m1testop4.test()

module m1testop5
using FinEtools
using FinEtoolsDeforLinear.DeforModelRedModule: DeforModelRed3D
using FinEtoolsDeforNonlinear
using FinEtoolsDeforNonlinear.MatDeforNeohookeanModule: MatDeforNeohookean
using FinEtoolsDeforNonlinear.FEMMDeforNonlinearModule: FEMMDeforNonlinear
using FinEtoolsDeforNonlinear.FEMMDeforNonlinearBaseModule: stiffness
using LinearAlgebra: norm
using SparseArrays
using DelimitedFiles
using Test
function test()
    mr = DeforModelRed3D
    E, nu = 7.0*phun("MPa"), 0.3
    m = MatDeforNeohookean(mr, E, nu)
    L= 6/2*phun("mm");
    H = 2/2*phun("mm");
    W = 2/2*phun("mm");
    fens, fes = H8hexahedron([0 0 0; 1.2 0 -0.1; 1.0 0.9 0.1; -0.05 1.1 0; 0 0 1.03; 1.2 0 1.1; 1.0 0.95 0.81; -0.05 1.01 0.9], 1, 1, 1)
    femm = FEMMDeforNonlinear(mr, IntegDomain(fes, GaussRule(3, 2)), m)

    geom = NodalField(fens.xyz)
    un1 = NodalField(zeros(size(fens.xyz, 1), 3)) # displacement field
    un = NodalField(zeros(size(fens.xyz, 1), 3)) # displacement field
    numberdofs!(un1)
    d(XYZ) = 0.1 * [XYZ[1]^2*XYZ[3], XYZ[2]^2*XYZ[3], XYZ[2]*XYZ[3]-XYZ[3]^3];
    for i = 1:count(fens)
        un1.values[i, :]  .+= d(fens.xyz[i, :]);
    end
    #println("un1.values = $(un1.values)")
    tn, dtn = 0.0, 1.0
    K = stiffness(femm, geom, un1, un, tn, dtn)
    ijs = readdlm("m1testop5.csv", ',', Float64)
    matK  = sparse(Int.(ijs[:, 1]), Int.(ijs[:, 2]), ijs[:, 3], size(K, 1), size(K, 2))
    # @show norm(K - matK) / maximum(matK[:])
    @test norm(K - matK) / maximum(matK[:]) < 1.0e-4
end
end
using .m1testop5
m1testop5.test()

module m1testop6
using FinEtools
using FinEtoolsDeforLinear.DeforModelRedModule: DeforModelRed3D
using FinEtoolsDeforNonlinear
using FinEtoolsDeforNonlinear.MatDeforNeohookeanModule: MatDeforNeohookean
using FinEtoolsDeforNonlinear.FEMMDeforNonlinearModule: FEMMDeforNonlinear
using FinEtoolsDeforNonlinear.FEMMDeforNonlinearBaseModule: stiffness
using LinearAlgebra: norm, cross
using SparseArrays
using DelimitedFiles
using Test
function test()
    mr = DeforModelRed3D
    E, nu = 7.0*phun("MPa"), 0.3
    m = MatDeforNeohookean(mr, E, nu)
    L= 6/2*phun("mm");
    H = 2/2*phun("mm");
    W = 2/2*phun("mm");
    fens, fes = H8hexahedron([0 0 0; 1.2 0 -0.1; 1.0 0.9 0.1; -0.05 1.1 0; 0 0 1.03; 1.2 0 1.1; 1.0 0.95 0.81; -0.05 1.01 0.9], 1, 1, 1)
    function  mcsys(csmatout::FFltMat, XYZ, tangents, fe_label)
        # @bp
        n1 = vec(XYZ); n1[1] = 0; n1  = n1 / norm(n1);
        n3 = vec([1, 0, 0]);
        n2 = cross(n3, n1);
        csmatout .= hcat(n1, n2, n3);
    end
    femm = FEMMDeforNonlinear(mr, IntegDomain(fes, GaussRule(3, 2)), CSys(3, 3, mcsys), m)

    geom = NodalField(fens.xyz)
    un1 = NodalField(zeros(size(fens.xyz, 1), 3)) # displacement field
    un = NodalField(zeros(size(fens.xyz, 1), 3)) # displacement field
    numberdofs!(un1)
    d(XYZ) = 0.1 * [XYZ[1]^2*XYZ[3], XYZ[2]^2*XYZ[3], XYZ[2]*XYZ[3]-XYZ[3]^3];
    for i = 1:count(fens)
        un1.values[i, :]  .+= d(fens.xyz[i, :]);
    end
    #println("un1.values = $(un1.values)")
    tn, dtn = 0.0, 1.0
    K = stiffness(femm, geom, un1, un, tn, dtn)
    ijs = readdlm("m1testop5.csv", ',', Float64)
    matK  = sparse(Int.(ijs[:, 1]), Int.(ijs[:, 2]), ijs[:, 3], size(K, 1), size(K, 2))
    # @show norm(K - matK) / maximum(matK[:])
    @test norm(K - matK) / maximum(matK[:]) < 1.0e-4
end
end
using .m1testop6
m1testop6.test()

module m1testop7
using FinEtools
using FinEtoolsDeforLinear.DeforModelRedModule: DeforModelRed3D
using FinEtoolsDeforNonlinear
using FinEtoolsDeforNonlinear.MatDeforNeohookeanModule: MatDeforNeohookean
using FinEtoolsDeforNonlinear.FEMMDeforNonlinearModule: FEMMDeforNonlinear
using FinEtoolsDeforNonlinear.FEMMDeforNonlinearBaseModule: stiffness
using LinearAlgebra: norm, cross
using SparseArrays
using DelimitedFiles
using Test
function test()
    mr = DeforModelRed3D
    E, nu = 7.0*phun("MPa"), 0.3
    m = MatDeforNeohookean(mr, E, nu)
    L= 6/2*phun("mm");
    H = 2/2*phun("mm");
    W = 2/2*phun("mm");
    fens, fes = H8hexahedron([0 0 0; 1.2 0 -0.1; 1.0 0.9 0.1; -0.05 1.1 0; 0 0 1.03; 1.2 0 1.1; 1.0 0.95 0.81; -0.05 1.01 0.9], 1, 1, 1)
    Rmout = fill(0.0, 3, 3)
    rv = vec([-0.56 -0.1361 0.35])
    rotmat3!(Rmout, rv)
    mcsys = CSys(Rmout)
    femm = FEMMDeforNonlinear(mr, IntegDomain(fes, GaussRule(3, 2)), mcsys, m)

    geom = NodalField(fens.xyz)
    un1 = NodalField(zeros(size(fens.xyz, 1), 3)) # displacement field
    un = NodalField(zeros(size(fens.xyz, 1), 3)) # displacement field
    numberdofs!(un1)
    d(XYZ) = 0.1 * [XYZ[1]^2*XYZ[3], XYZ[2]^2*XYZ[3], XYZ[2]*XYZ[3]-XYZ[3]^3];
    for i = 1:count(fens)
        un1.values[i, :]  .+= d(fens.xyz[i, :]);
    end
    #println("un1.values = $(un1.values)")
    tn, dtn = 0.0, 1.0
    K = stiffness(femm, geom, un1, un, tn, dtn)
    ijs = readdlm("m1testop5.csv", ',', Float64)
    matK  = sparse(Int.(ijs[:, 1]), Int.(ijs[:, 2]), ijs[:, 3], size(K, 1), size(K, 2))
    # @show norm(K - matK) / maximum(matK[:])
    @test norm(K - matK) / maximum(matK[:]) < 1.0e-4
end
end
using .m1testop7
m1testop7.test()

module m1testop7a
using FinEtools
using FinEtoolsDeforLinear.DeforModelRedModule: DeforModelRed3D
using FinEtoolsDeforNonlinear
using FinEtoolsDeforNonlinear.MatDeforNeohookeanModule: MatDeforNeohookean
using FinEtoolsDeforNonlinear.FEMMDeforNonlinearModule: FEMMDeforNonlinear
using FinEtoolsDeforNonlinear.FEMMDeforNonlinearBaseModule: restoringforce
using LinearAlgebra: norm, cross
using SparseArrays
using DelimitedFiles
using Test
function test()
    mr = DeforModelRed3D
    E, nu = 7.0*phun("MPa"), 0.3
    m = MatDeforNeohookean(mr, E, nu)
    L= 6/2*phun("mm");
    H = 2/2*phun("mm");
    W = 2/2*phun("mm");
    fens, fes = H8hexahedron([0 0 0; 1.2 0 -0.1; 1.0 0.9 0.1; -0.05 1.1 0; 0 0 1.03; 1.2 0 1.1; 1.0 0.95 0.81; -0.05 1.01 0.9], 1, 1, 1)
    Rmout = fill(0.0, 3, 3)
    rv = vec([-0.56 -0.1361 0.35])
    rotmat3!(Rmout, rv)
    mcsys = CSys(Rmout)
    femm = FEMMDeforNonlinear(mr, IntegDomain(fes, GaussRule(3, 2)), mcsys, m)

    geom = NodalField(fens.xyz)
    un1 = NodalField(zeros(size(fens.xyz, 1), 3)) # displacement field
    un = NodalField(zeros(size(fens.xyz, 1), 3)) # displacement field
    numberdofs!(un1)
    d(XYZ) = 0.1 * [XYZ[1]^2*XYZ[3], XYZ[2]^2*XYZ[3], XYZ[2]*XYZ[3]-XYZ[3]^3];
    for i = 1:count(fens)
        un1.values[i, :]  .+= d(fens.xyz[i, :]);
    end
    # println("un1.values = $(un1.values)")
    tn, dtn = 0.0, 1.0
    F = restoringforce(femm, geom, un1, un, tn, dtn)
    matF = readdlm("m1testop7.csv", ',', Float64)
    # @show norm(F - matF) / maximum(matF)
    @test norm(F - matF) / maximum(matF) < 1.0e-4
end
end
using .m1testop7a
m1testop7a.test()

module m1testop8
using FinEtools
using FinEtoolsDeforLinear.DeforModelRedModule: DeforModelRed3D
using FinEtoolsDeforNonlinear
using FinEtoolsDeforNonlinear.MatDeforNeohookeanModule: MatDeforNeohookean
using FinEtoolsDeforNonlinear.FEMMDeforNonlinearModule: FEMMDeforNonlinear
using FinEtoolsDeforNonlinear.FEMMDeforNonlinearBaseModule: restoringforce
using LinearAlgebra: norm, cross
using SparseArrays
using DelimitedFiles
using Test
function test()
    mr = DeforModelRed3D
    E, nu = 7.0*phun("MPa"), 0.3
    m = MatDeforNeohookean(mr, E, nu)
    L= 6/2*phun("mm");
    H = 2/2*phun("mm");
    W = 2/2*phun("mm");
    fens, fes = H8hexahedron([0 0 0; 1.2 0 -0.1; 1.0 0.9 0.1; -0.05 1.1 0; 0 0 1.03; 1.2 0 1.1; 1.0 0.95 0.81; -0.05 1.01 0.9], 1, 1, 1)
    function  mcsys(csmatout::FFltMat, XYZ, tangents, fe_label)
        # @bp
        n1 = vec(XYZ); n1[1] = 0; n1  = n1 / norm(n1);
        n3 = vec([1, 0, 0]);
        n2 = cross(n3, n1);
        csmatout .= hcat(n1, n2, n3);
    end
    femm = FEMMDeforNonlinear(mr, IntegDomain(fes, GaussRule(3, 2)), CSys(3, 3, mcsys), m)

    geom = NodalField(fens.xyz)
    un1 = NodalField(zeros(size(fens.xyz, 1), 3)) # displacement field
    un = NodalField(zeros(size(fens.xyz, 1), 3)) # displacement field
    numberdofs!(un1)
    d(XYZ) = 0.1 * [XYZ[1]^2*XYZ[3], XYZ[2]^2*XYZ[3], XYZ[2]*XYZ[3]-XYZ[3]^3];
    for i = 1:count(fens)
        un1.values[i, :]  .+= d(fens.xyz[i, :]);
    end
    # println("un1.values = $(un1.values)")
    tn, dtn = 0.0, 1.0
    F = restoringforce(femm, geom, un1, un, tn, dtn)
    matF = readdlm("m1testop7.csv", ',', Float64)
    # @show norm(F - matF) / maximum(matF)
    @test norm(F - matF) / maximum(matF) < 1.0e-4
end
end
using .m1testop8
m1testop8.test()

module m1testop9
using FinEtools
using FinEtoolsDeforLinear.DeforModelRedModule: DeforModelRed3D
using FinEtoolsDeforNonlinear
using FinEtoolsDeforNonlinear.MatDeforNeohookeanModule: MatDeforNeohookean
using FinEtoolsDeforNonlinear.FEMMDeforNonlinearModule: FEMMDeforNonlinear
using FinEtoolsDeforNonlinear.FEMMDeforNonlinearBaseModule: restoringforce
using LinearAlgebra: norm, cross
using SparseArrays
using DelimitedFiles
using Test
function test()
    mr = DeforModelRed3D
    E, nu = 7.0*phun("MPa"), 0.3
    m = MatDeforNeohookean(mr, E, nu)
    L= 6/2*phun("mm");
    H = 2/2*phun("mm");
    W = 2/2*phun("mm");
    fens, fes = H8hexahedron([0 0 0; 1.2 0 -0.1; 1.0 0.9 0.1; -0.05 1.1 0; 0 0 1.03; 1.2 0 1.1; 1.0 0.95 0.81; -0.05 1.01 0.9], 1, 1, 1)
    femm = FEMMDeforNonlinear(mr, IntegDomain(fes, GaussRule(3, 2)), m)

    geom = NodalField(fens.xyz)
    un1 = NodalField(zeros(size(fens.xyz, 1), 3)) # displacement field
    un = NodalField(zeros(size(fens.xyz, 1), 3)) # displacement field
    numberdofs!(un1)
    d(XYZ) = 0.1 * [XYZ[1]^2*XYZ[3], XYZ[2]^2*XYZ[3], XYZ[2]*XYZ[3]-XYZ[3]^3];
    for i = 1:count(fens)
        un1.values[i, :]  .+= d(fens.xyz[i, :]);
    end
    # println("un1.values = $(un1.values)")
    tn, dtn = 0.0, 1.0
    F = restoringforce(femm, geom, un1, un, tn, dtn)
    matF = readdlm("m1testop7.csv", ',', Float64)
    # @show norm(F - matF) / maximum(matF)
    @test norm(F - matF) / maximum(matF) < 1.0e-4
end
end
using .m1testop9
m1testop9.test()

module m1testop10
using FinEtools
using FinEtoolsDeforLinear.DeforModelRedModule: DeforModelRed3D
using FinEtoolsDeforNonlinear
using FinEtoolsDeforNonlinear.MatDeforNeohookeanModule: MatDeforNeohookean
using FinEtoolsDeforNonlinear.FEMMDeforNonlinearModule: FEMMDeforNonlinear
using FinEtoolsDeforNonlinear.FEMMDeforNonlinearBaseModule: nzebcloads
using LinearAlgebra: norm
using SparseArrays
using DelimitedFiles
using Test
function test()
    mr = DeforModelRed3D
    E, nu = 7.0*phun("MPa"), 0.3
    m = MatDeforNeohookean(mr, E, nu)
    L= 6/2*phun("mm");
    H = 2/2*phun("mm");
    W = 2/2*phun("mm");
    fens, fes = H8block(L, W, H, 1, 1, 1)
    femm = FEMMDeforNonlinear(mr, IntegDomain(fes, GaussRule(3, 2)), m)

    geom = NodalField(fens.xyz)
    un1 = NodalField(zeros(size(fens.xyz, 1), 3)) # displacement field
    un = NodalField(zeros(size(fens.xyz, 1), 3)) # displacement field
    for i in [3, 7, 8]
        for j in 1:3
            setebc!(un1, [i], true, j, 0.0)
        end
    end
    numberdofs!(un1)
    du = deepcopy(un1)
    tn, dtn = 0.0, 1.0
    d(XYZ) = 100.1 * [XYZ[1]^2*XYZ[3], XYZ[2]^2*XYZ[3], XYZ[2]*XYZ[3]-XYZ[3]^3];
    for i = 1:count(fens)
        du.fixed_values[i, :]  .+= d(fens.xyz[i, :]);
    end
    # @show du
    F = nzebcloads(femm, geom, un1, un, du, tn, dtn)
    matF = readdlm("m1testop10.csv", ',', Float64)
    # @show norm(K - matK) / maximum(matK[:])
    @test norm(F - matF) / maximum(matF[:]) < 1.0e-4
end
end
using .m1testop10
m1testop10.test()

module m1testop11
using FinEtools
using FinEtoolsDeforLinear.DeforModelRedModule: DeforModelRed3D
using FinEtoolsDeforNonlinear
using FinEtoolsDeforNonlinear.MatDeforNeohookeanModule: MatDeforNeohookean
using FinEtoolsDeforNonlinear.FEMMDeforNonlinearModule: FEMMDeforNonlinear
using FinEtoolsDeforNonlinear.FEMMDeforNonlinearBaseModule: geostiffness
using LinearAlgebra: norm
using SparseArrays
using DelimitedFiles
using Test
function test()
    mr = DeforModelRed3D
    E, nu = 7.0*phun("MPa"), 0.3
    m = MatDeforNeohookean(mr, E, nu)
    L= 6/2*phun("mm");
    H = 2/2*phun("mm");
    W = 2/2*phun("mm");
    fens, fes = H8hexahedron([0 0 0; 1.2 0 -0.1; 1.0 0.9 0.1; -0.05 1.1 0; 0 0 1.03; 1.2 0 1.1; 1.0 0.95 0.81; -0.05 1.01 0.9], 1, 1, 1)
    Rmout = fill(0.0, 3, 3)
    rv = vec([-0.56 -0.1361 0.35])
    rotmat3!(Rmout, rv)
    mcsys = CSys(Rmout)
    femm = FEMMDeforNonlinear(mr, IntegDomain(fes, GaussRule(3, 2)), mcsys, m)

    geom = NodalField(fens.xyz)
    un1 = NodalField(zeros(size(fens.xyz, 1), 3)) # displacement field
    un = NodalField(zeros(size(fens.xyz, 1), 3)) # displacement field
    numberdofs!(un1)
    d(XYZ) = 0.1 * [XYZ[1]^2*XYZ[3], XYZ[2]^2*XYZ[3], XYZ[2]*XYZ[3]-XYZ[3]^3];
    for i = 1:count(fens)
        un1.values[i, :]  .+= d(fens.xyz[i, :]);
    end
    #println("un1.values = $(un1.values)")
    tn, dtn = 0.0, 1.0
    Kgeo = geostiffness(femm, geom, un1, un, tn, dtn)
    # @show Kgeo[5,5]
    ijs = readdlm("m1testop11.csv", ',', Float64)
    matK  = sparse(Int.(ijs[:, 1]), Int.(ijs[:, 2]), ijs[:, 3], size(Kgeo, 1), size(Kgeo, 2))
    # @show norm(K - matK) / maximum(matK[:])
    @test norm(Kgeo - matK) / maximum(matK[:]) < 1.0e-4
end
end
using .m1testop11
m1testop11.test()

module m1testop12
using FinEtools
using FinEtoolsDeforLinear.DeforModelRedModule: DeforModelRed3D
using FinEtoolsDeforNonlinear
using FinEtoolsDeforNonlinear.MatDeforNeohookeanModule: MatDeforNeohookean
using FinEtoolsDeforNonlinear.FEMMDeforNonlinearModule: FEMMDeforNonlinear
using FinEtoolsDeforNonlinear.FEMMDeforNonlinearBaseModule: geostiffness
using LinearAlgebra: norm
using SparseArrays
using DelimitedFiles
using Test
function test()
    mr = DeforModelRed3D
    E, nu = 7.0*phun("MPa"), 0.3
    m = MatDeforNeohookean(mr, E, nu)
    L= 6/2*phun("mm");
    H = 2/2*phun("mm");
    W = 2/2*phun("mm");
    fens, fes = H8hexahedron([0 0 0; 1.2 0 -0.1; 1.0 0.9 0.1; -0.05 1.1 0; 0 0 1.03; 1.2 0 1.1; 1.0 0.95 0.81; -0.05 1.01 0.9], 1, 1, 1)
    # Rmout = fill(0.0, 3, 3)
    # rv = vec([-0.56 -0.1361 0.35])
    # rotmat3!(Rmout, rv)
    # mcsys = CSys(Rmout)
    # femm = FEMMDeforNonlinear(mr, IntegDomain(fes, GaussRule(3, 2)), mcsys, m)
    femm = FEMMDeforNonlinear(mr, IntegDomain(fes, GaussRule(3, 2)), m)

    geom = NodalField(fens.xyz)
    un1 = NodalField(zeros(size(fens.xyz, 1), 3)) # displacement field
    un = NodalField(zeros(size(fens.xyz, 1), 3)) # displacement field
    numberdofs!(un1)
    d(XYZ) = 0.1 * [XYZ[1]^2*XYZ[3], XYZ[2]^2*XYZ[3], XYZ[2]*XYZ[3]-XYZ[3]^3];
    for i = 1:count(fens)
        un1.values[i, :]  .+= d(fens.xyz[i, :]);
    end
    #println("un1.values = $(un1.values)")
    tn, dtn = 0.0, 1.0
    Kgeo = geostiffness(femm, geom, un1, un, tn, dtn)
    # @show Kgeo[5,5]
    ijs = readdlm("m1testop11.csv", ',', Float64)
    matK  = sparse(Int.(ijs[:, 1]), Int.(ijs[:, 2]), ijs[:, 3], size(Kgeo, 1), size(Kgeo, 2))
    # @show norm(K - matK) / maximum(matK[:])
    @test norm(Kgeo - matK) / maximum(matK[:]) < 1.0e-4
end
end
using .m1testop12
m1testop12.test()

ext: jl


text:
using DiffEqProblemLibrary, DiffEqBase, DASSL
using DiffEqProblemLibrary.DAEProblemLibrary: importdaeproblems; importdaeproblems()
using DiffEqProblemLibrary.DAEProblemLibrary: prob_dae_resrob

prob = prob_dae_resrob

sol = solve(prob,dassl())

sol = solve(prob,dassl(),abstol = 1e-1, reltol = 1e-2)


text:
# Example taken from Golightly and Wilkinson 2005 `Bayesian Inference for Stochastic
# Kinetic Models Using a Diffusion Approximation`

using Bridge
using StaticArrays, LinearAlgebra
import Bridge: b, σ, B, β, a, constdiff

struct Prokaryote{T} <: ContinuousTimeProcess{ℝ{4,T}}
    c₁::T
    c₂::T
    c₃::T
    c₄::T
    c₅::T
    c₆::T
    c₇::T
    c₈::T
    K::Float64

    function Prokaryote(c₁::T, c₂::T, c₃::T, c₄::T, c₅::T, c₆::T, c₇::T,
                        c₈::T, K) where T
        new{T}(c₁, c₂, c₃, c₄, c₅, c₆, c₇, c₈, K)
    end
end

struct ProkaryoteAux{O,R,S1,S2,S3} <: ContinuousTimeProcess{ℝ{4,R}}
    c₁::R
    c₂::R
    c₃::R
    c₄::R
    c₅::R
    c₆::R
    c₇::R
    c₈::R
    K::Float64
    t::Float64
    u::S1
    T::Float64
    v::S2
    aux::SArray{Tuple{4},S3,1,4}

    function ProkaryoteAux(c₁::R, c₂::R, c₃::R, c₄::R, c₅::R, c₆::R, c₇::R,
                           c₈::R, K, t, u::S1, T, v::S2, ::O,
                           aux::SArray{Tuple{4},S3,1,4}=@SVector[1.0,2.0,3.0,4.0]
                           ) where {O,R,S1,S2,S3}
        new{O,R,S1,S2,S3}(c₁, c₂, c₃, c₄, c₅, c₆, c₇, c₈, K, t, u, T, v, aux)
    end

    function ProkaryoteAux(P::ProkaryoteAux{O,R,S1,S2,S3},
                           aux::SArray{Tuple{4},S3,1,4}) where {O,R,S1,S2,S3}
        new{O,R,S1,S2,S3}(P.c₁, P.c₂, P.c₃, P.c₄, P.c₅, P.c₆, P.c₇, P.c₈, P.K,
                          P.t, P.u, P.T, P.v, aux)
    end
end

const ProkaryoteTypes = Union{Prokaryote, ProkaryoteAux}
# x <-> (RNA, P, P₂, DNA)
@inline _aux₁(x, P::ProkaryoteTypes) = P.c₅*x[2]*(x[2]-1)
@inline _aux₂(x, P::ProkaryoteTypes) = P.c₆*x[3]
@inline _aux₃(x, P::ProkaryoteTypes) = P.c₂*(P.K - x[4])
@inline _aux₄(x, P::ProkaryoteTypes) = P.c₁*x[3]*x[4]
@inline _aux₁₂(x, P::ProkaryoteTypes) = 2.0*_aux₂(x, P) - _aux₁(x, P)
@inline _aux₃₄(x, P::ProkaryoteTypes) = _aux₃(x, P) - _aux₄(x, P)

function b(t, x, P::Prokaryote{T}) where T
    k₁ = _aux₁₂(x, P)
    k₂ = _aux₃₄(x, P)
    ℝ{4}(P.c₃*x[4] - P.c₇*x[1], P.c₄*x[1] + k₁ - P.c₈*x[2], k₂ - 0.5 * k₁, k₂)
end

function _σ_prokaryote(x, P)
    k₁ = sqrt(0.5*_aux₁(x, P))
    k₂ = sqrt(_aux₂(x, P))
    k₃ = sqrt(_aux₃(x, P))
    k₄ = sqrt(_aux₄(x, P))
    _O = zero(x[1])

    @SMatrix [_O  _O  sqrt(P.c₃*x[4]) _O _O _O -sqrt(P.c₇*x[1]) _O;
              _O _O _O sqrt(P.c₄*x[1]) -2.0*k₁ 2.0*k₂ _O -sqrt(P.c₈*x[2]);
              -k₄ k₃ _O _O k₁ -k₂ _O _O;
              -k₄ k₃ _O _O _O _O _O _O]
end

function _a_prokaryote(x, P)
    k₁ = _aux₁(x, P) + 2.0*_aux₂(x, P)
    k₂ = _aux₃(x, P) + _aux₄(x, P)
    s₁ = P.c₃*x[4]+P.c₇*x[1]
    s₂ = P.c₄*x[1]+P.c₈*x[2]
    _O = zero(x[1])

    @SMatrix [ s₁ _O _O _O;
              _O s₂+2*k₁ -k₁ _O;
              _O -k₁ 0.5*k₁+k₂ k₂;
              _O _O k₂ k₂]
end

σ(t, x, P::Prokaryote) = _σ_prokaryote(x, P)
a(t, x, P::Prokaryote) = _a_prokaryote(x, P)
domain(P::Prokaryote) = BoundedDomain( LowerBoundedDomain((0.0, 1.0, 0.0, 0.0),
                                                          (1,2,3,4)),
                                       UpperBoundedDomain((P.K,), (4,)) )
constdiff(::Prokaryote) = false
clone(P::Prokaryote, θ) = Prokaryote(θ..., P.K)
params(P::Prokaryote) = [P.c₁, P.c₂, P.c₃, P.c₄, P.c₅, P.c₆, P.c₇, P.c₈]

# x <-> (RNA, P, P₂, DNA)
phi(::Val{0}, t, x, P::Prokaryote) = (zero(x[1]), zero(x[1]), zero(x[1]), zero(x[1]))
phi(::Val{1}, t, x, P::Prokaryote) = (zero(x[1]), zero(x[1]), -x[3]*x[4], -x[3]*x[4])
phi(::Val{2}, t, x, P::Prokaryote) = (zero(x[1]), zero(x[1]), P.K-x[4], P.K-x[4])
phi(::Val{3}, t, x, P::Prokaryote) = (x[4], zero(x[1]), zero(x[1]), zero(x[1]))
phi(::Val{4}, t, x, P::Prokaryote) = (zero(x[1]), x[1], zero(x[1]), zero(x[1]))
phi(::Val{5}, t, x, P::Prokaryote) = (zero(x[1]), -x[2]*(x[2]-1), 0.5*x[2]*(x[2]-1), zero(x[1]))
phi(::Val{6}, t, x, P::Prokaryote) = (zero(x[1]), 2.0*x[3], zero(x[1]), zero(x[1]))
phi(::Val{7}, t, x, P::Prokaryote) = (-x[1], zero(x[1]), zero(x[1]), zero(x[1]))
phi(::Val{8}, t, x, P::Prokaryote) = (zero(x[1]), -x[2], zero(x[1]), zero(x[1]))

nonhypo(P::Prokaryote, x) = x
num_non_hypo(P::Type{<:Prokaryote}) = 4

observables(::ProkaryoteAux{O}) where O = O()

function B(t, P::ProkaryoteAux{Val{(true,true,true,true)}})
    k₁ = P.c₅ - 2*P.c₅*P.v[2]
    k₂ = -P.c₁*P.v[4]
    k₃ = -P.c₂-P.c₁*P.v[3]
    @SMatrix [-P.c₇  0.0  0.0  P.c₃;
              P.c₄  k₁-P.c₈  2.0*P.c₆  0.0;
              0.0  -0.5*k₁  k₂-P.c₆  k₃;
              0.0  0.0  k₂ k₃]
end

function β(t, P::ProkaryoteAux{Val{(true,true,true,true)}})
    k₁ = P.c₂*P.K + P.c₁*P.v[3]*P.v[4]
    @SVector [0.0,
              P.c₅*P.v[2]^2,
              k₁ - 0.5*P.c₅*P.v[2]^2,
              k₁]
end

function B(t, P::ProkaryoteAux{Val{:custom}})
    k₁ = P.c₅ - 2*P.c₅*P.aux[2]
    k₂ = -P.c₁*P.aux[4]
    k₃ = -P.c₂-P.c₁*P.aux[3]
    @SMatrix [-P.c₇  0.0  0.0  P.c₃;
              P.c₄  k₁-P.c₈  2.0*P.c₆  0.0;
              0.0  -0.5*k₁  k₂-P.c₆  k₃;
              0.0  0.0  k₂ k₃]
end

function β(t, P::ProkaryoteAux{Val{:custom}})
    k₁ = P.c₂*P.K + P.c₁*P.aux[3]*P.aux[4]
    @SVector [0.0,
              P.c₅*P.aux[2]^2,
              k₁ - 0.5*P.c₅*P.aux[2]^2,
              k₁]
end

#=
function B(t, P::ProkaryoteAux{Val{:custom}})
    @SMatrix [0.0  0.0  0.0  0.0;
              0.0  0.0  0.0  0.0;
              0.0  0.0  0.0  0.0;
              0.0  0.0  0.0 0.0]
end

function β(t, P::ProkaryoteAux{Val{:custom}})
    @SVector [0.0,
              0.0,
              0.0,
              0.0]
end
=#

σ(t, P::ProkaryoteAux{Val{(true,true,true,true)}}) = _σ_prokaryote(P.v, P)
a(t, P::ProkaryoteAux{Val{(true,true,true,true)}}) = _a_prokaryote(P.v, P)
σ(t, P::ProkaryoteAux{Val{:custom}}) = _σ_prokaryote(P.aux, P)
a(t, P::ProkaryoteAux{Val{:custom}}) = _a_prokaryote(P.aux, P)
a(t, x, P::ProkaryoteAux) = a(t, P)
σ(t, x, P::ProkaryoteAux) = σ(t, P)
b(t, x, P::ProkaryoteAux) = B(t, P)*x + β(t, P)

constdiff(::ProkaryoteAux) = true
clone(P::ProkaryoteAux, θ) = ProkaryoteAux(θ..., P.K, P.t, P.u, P.T, P.v, observables(P), P.aux)
clone(P::ProkaryoteAux, θ, v) = ProkaryoteAux(θ..., P.K, P.t, zero(v), P.T, v, observables(P), P.aux)
params(P::ProkaryoteAux) = [P.c₁, P.c₂, P.c₃, P.c₄, P.c₅, P.c₆, P.c₇, P.c₈]
depends_on_params(::ProkaryoteAux) = (1,2,3,4,5,6,7,8)

recentre(P::ProkaryoteAux{Val{:custom}}, tt, X̄) = ProkaryoteAux(P, X̄[end])
update_λ!(P::ProkaryoteAux{Val{:custom}}, λ) = nothing


text:
function V0inv(m::LinearMixedModel)

    Zt=MixedModels.zt(m)
    Lt=SelectMixedModels.Λ0(m)
    X=m.X.m
    #Vinv = eye(size(m)[1]) - Zt'*inv(Zt*Zt' + inv(Lt*Lt'))*Zt
    Vinv=vinv(m)

    return(Vinv - Vinv*X*inv(X'*Vinv*X)*X'*Vinv)
end

function sumAY(m::LinearMixedModel)
     Vi = V0inv(m)
     return( (m.y' * Vi, Vi*m.y, m.y' * Vi *m.y))
end

function MakeBlkDeriv(m::LinearMixedModel)

    vv = SparseMatrixCSC{Float64,Int}[]
    for i in 1:length(m.λ)

        nl = length(m.facs[i].pool)
        ll = m.λ[i]
        #p = size(ll,1)
        p = size(ll)[1]
        if p == 1
            isa(ll,MixedModels.PDScalF) || error("1×1 λ section should be a PDScalF type")
            push!(vv,m.Ztblks[i]'*speye(nl)*m.Ztblks[i])
        else
            if isa(ll, MixedModels.PDDiagF)
                for k=1:p
                    sfll = sparse(zeros(p,p))
                    sfll[k,k]=1.0
                    push!(vv,m.Ztblks[i]'*blkdiag([sfll for _ in 1:nl]...)*m.Ztblks[i])
                end
            else
                for k=1:p
                    for l=k:p
                        sfll = sparse(zeros(p,p))
                        if k==l
                            sfll[l,k]=1.0
                            push!(vv,m.Ztblks[i]'*blkdiag([sfll for _ in 1:nl]...)*m.Ztblks[i])
                        else
                            sfll[l,k]=1.0
                            sfll[k,l]=1.0
                            push!(vv,m.Ztblks[i]'*blkdiag([sfll for _ in 1:nl]...)*m.Ztblks[i])
                        end
                    end
                end
            end
        end
    end
    #Triangular(blkdiag(vv...),:L,false)
    vv
end

function yAWAy(m::LinearMixedModel, W::Array{SparseMatrixCSC{Float64,Int64},1} )
    p=length(W)
    e = m.y-m.μ
    ewe=zeros(p)

    for i in 1:p
       ewe[i] =  (e' * W[i] *e)[1]
    end

    return(ewe)
end

function conAICbc(mm::LinearMixedModel)

    llk = condll(mm)

    oldf=mm.f
    newf=checkzerovar(mm)

    if oldf ==newf
        m=mm
    else
        retrms = filter(x->Meta.isexpr(x,:call) && x.args[1] == :|, ModelFrame(newf,mm.mf.df).terms.terms)
        length(retrms) > 0 || return(SelectMixedModels.IC( 0.0 , Inf, "VBc", [0.0,0.0], Inf ))
        m=fit(lmm(newf,mm.mf.df))
    end

    Zt=MixedModels.zt(m)
    Lt=SelectMixedModels.Λ0(m)
    X=m.X.m

    try
        Vinv = eye(size(m)[1]) - Zt'*inv(Zt*Zt' + inv(Lt*Lt'))*Zt
    catch
        Vinv = eye(size(m)[1])
    end

    if countnz(Vinv) == size(m)[1]
        return(SelectMixedModels.IC( 0.0 , Inf, "VBc", [0.0,0.0], Inf ))
    else
        Vinv = eye(size(m)[1]) - Zt'*inv(Zt*Zt' + inv(Lt*Lt'))*Zt
    end

    A= Vinv - Vinv*X*inv(X'*Vinv*X)*X'*Vinv
    rho =size(m)[1]-trace(A)

    (nn,pp,qq,tt)=size(m)

    W=MakeBlkDeriv(m)
    eWe=yAWAy(m, W)
    e=m.y-m.μ
    tye= (m.y'*e)[1]
    theta=MixedModels.θ(m)
    s=length(theta)

    B=zeros(s,s)

    for j in 1:s
        for l in 1:j
            B[j,l] =  -tye^2 *trace(W[j]*Vinv*W[l]*Vinv)/nn
            B[j,l] -=  eWe[j] * eWe[l] +2.0* (e' * W[l] *A* W[j] *e)[1]*tye

            if j != l
                B[j,l] =B[l,j]
            end
         end
    end

    Binv=inv(B)

    G= zeros(s,nn)

    for j in 1:s
      G[j,:] = 2.0* (tye * (e' * W[j] *A - eWe[j] * e')[1] )
    end

    bc=0.0

    for j in 1:s
      bc += sum ( Binv[j,:]*G*A*W[j]*e )
    end

    df = rho+bc+1.0
    val = -2.0*llk+2.0*df

    cr=[1.0,1.0]

    SelectMixedModels.IC(llk, df, "VBc", cr, val )
end

max_stars_repo_name: ilovealldata/SelectMixedModels

max_issues_repo_name: ilovealldata/SelectMixedModels

max_forks_repo_name: ilovealldata/SelectMixedModels


text:
# Test if we can transform part of a Hessenberg matrix to Schur form.
# Here we look at some edge cases.

using Test, LinearAlgebra
using ArnoldiMethod: eigenvalues, local_schurfact!, is_offdiagonal_small,
                     NotWanted

include("utils.jl")

@testset "Schur factorization" begin
    let
        # 2-by-2 matrix with distinct eigenvalues while H[2,1] != 0
        H = [1.0 2.0; 3.0 4.0]
        H′ = copy(H)
        Q = Matrix{Float64}(I, 2, 2)

        @test local_schurfact!(H′, 1, 2, Q, eps(), 2)
        @test norm(H * Q - Q * H′) < 10eps()
        @test sort!(eigenvalues(H′), by=realimag) ≈ sort!(eigvals(H′), by=realimag)
        @test sort!(eigenvalues(H′), by=realimag) ≈ sort!(eigvals(H), by=realimag)
        @test iszero(H′[2,1])
    end

    let
        # 2-by-2 matrix with distinct eigenvalues while H[2,1] = 0
        H = [1.0 2.0; 0.0 4.0]
        H′ = copy(H)
        Q = Matrix{Float64}(I, 2, 2)

        @test local_schurfact!(H′, 1, 2, Q, eps(), 2)
        @test norm(H * Q - Q * H′) < 10eps()
        @test sort!(eigenvalues(H′), by=realimag) ≈ sort!(eigvals(H′), by=realimag)
        @test sort!(eigenvalues(H′), by=realimag) ≈ sort!(eigvals(H), by=realimag)
        @test iszero(H′[2,1])
    end

    let
        # 2-by-2 matrix with conjugate eigenvalues
        H = [1.0 4.0; -5.0 3.0]
        H′ = copy(H)
        Q = Matrix{Float64}(I, 2, 2)

        @test local_schurfact!(H′, 1, 2, Q, eps(), 2)
        @test norm(H * Q - Q * H′) < 10eps()
        @test sort!(eigenvalues(H′), by=realimag) ≈ sort!(eigvals(H′), by=realimag)
        @test sort!(eigenvalues(H′), by=realimag) ≈ sort!(eigvals(H), by=realimag)
    end

    # Larger real matrix.

    let
        n = 10
        # Transforming a 1+i by n-i block of the matrix H into upper triangular form
        for i = 0 : 4
            Q = Matrix{Float64}(I, n, n)
            H = triu(randn(Float64, n, n))
            H[1+i:n-i,1+i:n-i] = normal_hessenberg_matrix(Float64, i+1:n-i)
            H′ = copy(H)

            # Test that the procedure has converged
            @test local_schurfact!(H′, 1+i, n-i, Q)

            for j = 1+i:n-i-1
                t = H′[j,j] + H′[j+1,j+1]
                d = H′[j,j] * H′[j+1,j+1] - H′[j+1,j] * H′[j,j+1]

                # Test if subdiagonal is small. If not, check if conjugate eigenvalues.
                @test is_offdiagonal_small(H′, j) || t^2 < 4d
            end

            # Test that the elements below the subdiagonal are 0
            @test is_hessenberg(H′)

            # Test that the partial Schur decomposition relation holds
            @test norm(H * Q - Q * H′) < 1000eps()

            # Test that the eigenvalues of H are the same before and after transformation
            @test sort!(eigvals(H), by=realimag) ≈ sort!(eigvals(H′), by=realimag)
        end

    end

    # COMPLEX ARITHMETIC

    let
        n = 10
        # Transforming a 1+i by 10-i block of the matrix H into upper triangular form
        for i = 0 : 4
            Q = Matrix{ComplexF64}(I, n, n)
            H = triu(randn(ComplexF64, n, n))
            H[i+1:n-i,i+1:n-i] = normal_hessenberg_matrix(ComplexF64, (i+1:n-i) .* (1 + im))
            H′ = copy(H)

            # Test that the procedure has converged
            @test local_schurfact!(H′, 1+i, n-i, Q)

            # Test if subdiagonal is small.
            for j = 1+i:n-i-1
                @test iszero(H′[j+1,j])
            end

            # Test that the elements below the subdiagonal are 0
            @test is_hessenberg(H′)

            # Test that the partial Schur decomposition relation holds
            @test norm(H * Q - Q * H′) < 1000eps()

            # Test that the eigenvalues of H are the same before and after transformation
            @test sort!(eigvals(H), by=realimag) ≈ sort!(eigvals(H′), by=realimag)
        end
    end
end

@testset "Schur with nearly repeated eigenvalues" begin
    # Matrix with nearly repeated eigenpair could converge very slowly or
    # stagnate completely when there's a tiny perturbation in the computed
    # shift.
    mat(ε) = [2   0    0  ;
              5ε  1-ε  2ε ;
              0   3ε   1+ε]

    @test local_schurfact!(mat(eps()))
end


text:
#exprGenerator.jl

using Calculus, SymbolicUtils, SymbolicUtils.Code
import Base.*, Base.show

@syms x f(x) df(x)

struct tabularFunction
	vals::Matrix{Int64}
	x0::Int64
	baseSym::Char
	tabularFunction(x,y)=new(rand(UnitRange(-10,10),(4,3)),x,y)
end

function Base.println(t::tabularFunction)
	println("x  $(t.baseSym)(x) $(t.baseSym)'(x) $(t.baseSym)''(x)")
	for i in 1:4
		println("$((t.x0)+i)  $(t.vals[i,:])")
	end
end

function charVector(s::String)
	v = Char[]
	for c in s
		push!(v,c)
	end
	v
end

function *(v::Vector{Char})
	s = ""
	for i in v
		s=s*i
	end
	s
end

function compose(s::String,f::String)
	v,w = charVector(s),charVector(f)
	final = Char[]

	index = 1

	for c in v
		if c == 'x'
			push!(final,'(')
			for d in w
				push!(final,d)
			end
			push!(final,')')
		else
		push!(final,c)
		end
		index += 1
	end
	*(final)
end


text:
push!(LOAD_PATH,"../src/") # load module from local directory

include("../src/PolarisableIonTightBinding.jl")
using Plots

# Plot pair potential; O-H bond integrals
plot([PITB.OOPP, PITB.OHVssσ, PITB.OHVpsσ],1.0,4.0)


text:
@doc raw"""
Spin-``S`` XXZ model represented as the following Hamiltonian,

``
\mathcal{H} = \sum_{i,j} \left[ J_{ij}^z S_i^z S_j^z
            + \frac{J_{ij}^{xy}}{2} (S_i^+ S_j^- + S_i^-S_j^+) \right]
            - \sum_i \Gamma_i S_i^x,
``

where ``S^x, S^y, S^z`` are ``x, y`` and ``z`` component of spin operator with length ``S``,
and ``S^\pm \equiv S^x \pm iS^y`` are ladder operator.
A state is represented by a product state (spins at ``\tau=0``) of local ``S^z`` diagonal basis and an operator string (perturbations).
A local spin with length ``S`` is represented by a symmetrical summation of ``2S`` sub spins with length ``1/2``.
"""
mutable struct QuantumXXZ <: Model
    lat :: Lattice
    S2 :: Int
    spins :: Matrix{Int}
    ops :: Vector{LocalLoopOperator}
    rng :: Random.MersenneTwister

    function QuantumXXZ(lat::Lattice, S::Real, rng::Random.AbstractRNG)
        if round(2S) != 2S
            error("`S` should be integer or half-integer")
        end
        S2 = round(Int,2S)
        model = new()
        model.lat = lat
        model.rng = rng
        model.S2 = S2
        model.spins = rand(model.rng,[1,-1], 1, numsites(lat)*S2)
        model.ops = LocalLoopOperator[]
        return model
    end
end
QuantumXXZ(lat::Lattice, S::Real) = QuantumXXZ(lat, S, Random.seed!(Random.MersenneTwister(0)))
QuantumXXZ(lat::Lattice, S::Real, seed) = QuantumXXZ(lat, S, Random.seed!(Random.MersenneTwister(0), seed...))

@doc doc"""
    QuantumXXZ(param)

Generates `QuantumXXZ` using `param["Lattice"]`, `param["S"]` and `param["Seed"]` (if defined).
Each subspin will be initialized independently and randomly.
"""
function QuantumXXZ(param::Parameter)
    lat = generatelattice(param)
    S = param["S"]
    if "Seed" in keys(param)
        return QuantumXXZ(lat,S,param["Seed"])
    else
        return QuantumXXZ(lat,S)
    end
end

@inline site2subspin(site::Integer, ss::Integer, S2::Integer) = (site-1)*S2+ss
@inline function subspin2site(subspin::Integer, S2::Integer)
    return ceil(Int, subspin/S2), mod1(subspin,S2)
end

include("update.jl")
include("estimator.jl")
include("postproc.jl")


text: @functor Base.RefValue


text:
module OpenStreetMapIO

include("map_types.jl")
include("io_pbf.jl")
include("io_xml.jl")

export readpbf, readosm, queryoverpass
export Map, Node, Way, Relation, BBox, LatLon

end


text:
module HWfuncapp

using FastGaussQuadrature # to get chebyshevnodes

# you dont' have to use PyPlot. I did much of it in PyPlot, hence
# you will see me often qualify code with `PyPlot.plot` or so.
# using PyPlot
import ApproXD: getBasis, BSpline
using Distributions
using BasisMatrices
using LinearAlgebra, SpecialFunctions
using ApproxFun
using Plots

export q1, q2, q3, q7

ChebyT(x,deg) = cos(acos(x)*deg)
unitmap(x,lb,ub) = 2 .* (x .- lb) ./ (ub .- lb) .- 1	#[a,b] -> [-1,1]
abmap(x,lb,ub) = 0.5 .* (ub .+ lb) .+ 0.5 .* (ub .- lb) .* x	# [-1,1] -> [a,b]
function chebpol(x, n)
	Φ = zeros(size(x,1),n) # chebyshev polynomials basis
	for i in 1:size(x,1)
		for j in 1:n
			Φ[i,j] = ChebyT(x[i],j-1)
			# Φ[i,j] = abmap(Φ[i,j], lb, ub)
		end
	end
	return Φ
end

function q1(n=15)
	lb = -3
	ub = 3
	x = range(lb, stop = ub ,length = n)
	# generate Chebyshev nodes
	S, y = gausschebyshev(n)
	# Scale Chebyshev nodes on -3 to 3 domain
	ϕ = abmap(S, lb, ub) # Chebyshev nodes on -3, 3
	# define function
	f(x) = x .+ 2x.^2 - exp.(-x)
	# Evaluate function at chebyshev nodes
	Y = f(ϕ)

	# Build Function grid
	Φ = chebpol(S, n)

	# Estimate parameters of approximation
	c = Φ\Y

	# Evaluate approximation
	Yhat = Φ*c

	# Test accuracy on larger sample
	n_new = range(lb, stop = ub ,length = 100)
	# Evaluate function
	Y_new = f(n_new)
	# Build Polynomial grid
	Φ_new = chebpol(unitmap(n_new,lb,ub), n)

	# Evaluate using fitted coefficients
	Yhat_new = Φ_new*c

	# Plot
	p = Any[]
	push!(p,Plots.plot(n_new, [Y_new, Yhat_new], label=["True value" "Approximation"], marker=([:none :diamond])))
	err = Y_new .- Yhat_new
	push!(p, Plots.plot(n_new, err, label="Approximation error"))
	Plots.plot(p...)
	# without using PyPlot, just erase the `PyPlot.` part
	Plots.savefig(Plots.plot(p...), joinpath(dirname(@__FILE__),"..","q1.png"))
	return Dict("error"=>maximum(abs,err))
end

function q2(b=4)
	@assert b > 0
	n = 15
	# use ApproxFun.jl to do the same:
	deg = n - 1
	ub = b
	lb = -b
	S = Chebyshev(lb..ub)
	p = range(lb,stop=ub,length=n)  # a non-default grid
	# define function
	f(x) = x .+ 2x.^2 - exp.(-x)
	# evaluate function
	v = f(p)           # values at the non-default grid
	V = Array{Float64}(undef,n,deg + 1) # Create a Vandermonde matrix by evaluating the basis at the grid

	for k = 1:deg+1
	    V[:,k] = Fun(S,[zeros(k-1);1]).(p) # Evaluate Chebyshev basis at p
		# [zeros(k-1);1] is the identity matrix, all observations have same k degree Chebyshev basis function, evaluated at a different point
	end
	V
	g = Fun(S,V\v);
	@show g(1.1)
	@show f(1.1)

	n_new = range(lb, stop = ub, length = 100)
	p = Any[]
	push!(p,Plots.plot(n_new, [f(n_new), g.(n_new)], label=["True value" "Approximation"], marker=([:none :diamond])))
	err = f(n_new) .- g.(n_new)
	push!(p, Plots.plot(n_new, err, label="Approximation error"))
	fig = Plots.plot(p...)

	Plots.savefig(fig,joinpath(dirname(@__FILE__),"..","q2.png"))
end

function q3(b=10)

	x = Fun(identity,-b..b)
	f = sin(x^2)
	g = cos(x)
	h = f - g
	r = roots(h)

	Plots.plot(h, label="h(x)")
	Plots.scatter!(r,h.(r), label="Roots")
	Plots.savefig(joinpath(dirname(@__FILE__),"..","q3.png"))
	# xbis = Fun(identity,-b..0)
	g = cumsum(h) # indefinite  integral
	g = g + h(-b) # definite integral with constant of integration
	integral = norm(g(0) - g(-b)) # definite integral from -b to 0
	# p is your plot
	return (integral)
end

# optinal
function q4()

	return fig
end

# I found having those useful for q5
mutable struct ChebyType
	f::Function # fuction to approximate
	nodes::Union{Vector,LinRange} # evaluation points
	basis::Matrix # basis evaluated at nodes
	coefs::Vector # estimated coefficients

	deg::Int 	# degree of chebypolynomial
	lb::Float64 # bounds
	ub::Float64

	# constructor
	function ChebyType(_nodes::Union{Vector,LinRange},_deg,_lb,_ub,_f::Function)
		n = length(_nodes)
		y = _f(_nodes)
		_basis = Float64[ChebyT(unitmap(_nodes[i],_lb,_ub),j) for i=1:n,j=0:_deg]
		_coefs = _basis \ y  # type `?\` to find out more about the backslash operator. depending the args given, it performs a different operation
		# create a ChebyComparer with those values
		new(_f,_nodes,_basis,_coefs,_deg,_lb,_ub)
	end
end

# function to predict points using info stored in ChebyType
function predict(Ch::ChebyType,x_new)

	true_new = Ch.f(x_new)
	basis_new = Float64[ChebyT(unitmap(x_new[i],Ch.lb,Ch.ub),j) for i=1:length(x_new),j=0:Ch.deg]
	basis_nodes = Float64[ChebyT(unitmap(Ch.nodes[i],Ch.lb,Ch.ub),j) for i=1:length(Ch.nodes),j=0:Ch.deg]
	preds = basis_new * Ch.coefs
	preds_nodes = basis_nodes * Ch.coefs

	return Dict("x"=> x_new,"truth"=>true_new, "preds"=>preds, "preds_nodes" => preds_nodes)
end

function q5(deg=(5,9,15),lb=-1.0,ub=1.0)

	runge(x) = 1.0 ./ (1 .+ 25 .* x.^2)

	PyPlot.savefig(joinpath(dirname(@__FILE__),"..","q5.png"))

end

function q6()

	# compare 2 knot vectors with runge's function

	PyPlot.savefig(joinpath(dirname(@__FILE__),"..","q6.png"))

end

function q7(n=13)
	f(x) = abs.(x).^0.5

	# Regular grid
	bs = BSpline(13,3,-1,1) #13 knots, degree 3 in [0,1], no multiplicity
	# 3 knots grid
	multiknots = vcat(range(-1,stop = -0.1,length = 5),0,0,0, range(0.1,stop = 1,length =5))
	bs2 = BSpline(multiknots, 3)

	# Evaluate function of interest
	x = range(-1,stop =1.0, length = 65)
	Y = f(x)

	# Evaluate basis on grid
	B = Array(getBasis(collect(x),bs))
	B2 = Array(getBasis(collect(x),bs2))

	# Solve for parameters
	c = B\Y
	c2 = B2\Y

	# Compute fitted values
	Yhat = B*c
	Yhatbis = B2*c2

	# Compute fitting errors
	err = Y .- Yhat
	errbis = Y .- Yhatbis

	# Plot
	p = Any[]
	push!(p,Plots.plot(x, Y, title = "True function"))
	push!(p,Plots.plot(x, [Yhat, Yhatbis], label = ["uniform" "multiplicity"], title="Approximations"))
	push!(p,Plots.plot(x, [err, errbis], label = ["uniform" "multiplicity"], title="Errors"))
	return Plots.plot(p...)
	Plots.savefig(joinpath(dirname(@__FILE__),"..","q7.png"))
end

	# function to run all questions
function runall()
	@info("running all questions of HW-funcapprox:")
	q1(15)
	q2(3)
	q3(10)
	q4()
	q5()
	q6()
	q7()
end

end # module


text:
module BMTest

using DelimitedFiles
using LinearAlgebra
using Random
using Test
using Statistics

import BoltzmannMachines
const BMs = BoltzmannMachines

macro check(expr)
   msg = string(expr) # TODO: print vales
   quote
      if !($(esc(expr)))
         @warn "Check failed: " $msg
         return false
      else
         true # no return
      end
   end
end

function createsamples(nsamples::Int, nvariables::Int, samerate::Float64 = 0.7)

   x = round.(rand(nsamples, nvariables))
   samerows = ((x[:,1] .== x[:,2]) .& (x[:,2] .== x[:,3]))
   sameratestart = count(samerows) / nsamples

   sameratediff = samerate - sameratestart

   if sameratediff > 0
      addedsamerange = findall(.!samerows)
      shuffle!(addedsamerange)
      nsameadded = round(Int, sameratediff * nsamples)
      addedsamerange = addedsamerange[1:nsameadded]
      x[addedsamerange, 3] = x[addedsamerange, 2] = x[addedsamerange, 1]
   else
      subtractedsamerange = find(samerows)
      shuffle!(subtractedsamerange)
      nsamesubtracted = round(Int, - sameratediff * nsamples)
      subtractedsamerange = subtractedsamerange[1:nsamesubtracted]
      for row in subtractedsamerange
         while ((x[row,1] == x[row,2]) & (x[row,2] == x[row,3]))
            xrow = view(x, row, 1:3)
            rand!(xrow)
            map!(round, xrow, xrow)
         end
      end
   end
   x
end

function randrbm(nvisible, nhidden, factorw = 1.0, factora = 1.0, factorb = 1.0)
   w = factorw*randn(nvisible, nhidden)
   a = factora*rand(nvisible)
   b = factorb*(0.5 .- rand(nhidden))
   BMs.BernoulliRBM(w, a, b)
end

function randgbrbm(nvisible, nhidden, factorw = 1.0, factora = 1.0, factorb = 1.0, factorsd = 1.0)
   w = factorw*randn(nvisible, nhidden)
   a = factora*rand(nvisible)
   b = factorb*ones(nhidden)
   sd = factorsd*ones(nvisible)
   BMs.GaussianBernoulliRBM(w, a, b, sd)
end

function randbgrbm(nvisible, nhidden, factorw = 1.0, factora = 1.0, factorb = 1.0)
   w = factorw*randn(nvisible, nhidden)
   a = factora*rand(nvisible)
   b = factorb*ones(nhidden)
   BMs.BernoulliGaussianRBM(w, a, b)
end

function randdbm(nunits)
   nrbms = length(nunits) - 1
   dbm = BMs.BasicDBM(undef, nrbms)
   for i = 1:nrbms
      dbm[i] = randrbm(nunits[i], nunits[i+1])
   end
   dbm
end

"""
Tests the functions for computing the activation potentials
"""
function test_potentials()
   nvisible = 3
   nhidden = 2
   nsamples = 5
   vv = rand(nsamples, nvisible)

   # Test Test activation potential of hidden nodes for BernoulliRBM
   hh = rand(nsamples, nhidden)
   rbm = BMTest.randrbm(nvisible, nhidden)
   BMs.hiddenpotential!(hh, rbm, vv)
   @test isapprox(sum(abs.(hh - BMs.hiddenpotential(rbm, vv))), 0)

   # Test that activation potential of hidden nodes of a
   # Binomial2BernoulliRBM is the same as that of a BernoulliRBM
   b2brbm = BMs.Binomial2BernoulliRBM(rbm.weights, rbm.visbias, rbm.hidbias)
   hh2 = BMs.hiddenpotential(b2brbm, vv)
   @test sum(abs.(hh - hh2)) == 0
   BMs.hiddenpotential!(hh, b2brbm, vv)
   @test isapprox(sum(abs.(hh - hh2)), 0)

   # Test activation potential of visible nodes for BernoulliRBM
   BMs.visiblepotential!(vv, rbm, hh)
   @test isapprox(sum(abs.(vv - BMs.visiblepotential(rbm, hh))), 0)

   # Test activation potential of the visible nodes of a Binomial2BernoulliRBM
   BMs.visiblepotential!(vv, b2brbm, hh)
   @test isapprox(sum(abs.(vv- BMs.visiblepotential(b2brbm, hh))), 0)

   # Test activation potential of hidden nodes for GBRBM
   hh = rand(nsamples, nhidden)
   gbrbm = BMTest.randgbrbm(nvisible, nhidden)
   BMs.hiddenpotential!(hh, gbrbm, vv)
   @test isapprox(sum(abs.(hh - BMs.hiddenpotential(gbrbm, vv))), 0)

   # Test activation potential of visible nodes for GBRBM
   BMs.visiblepotential!(vv, gbrbm, hh)
   @test isapprox(sum(abs.(vv - BMs.visiblepotential(gbrbm, hh))), 0)

   # Test activation potential of hidden nodes for PartitionedRBM.
   # For a PartitionedRBM consisting of BernoulliRBMs, the
   # activation potential has to be the same as the activation potential
   # of the BernoulliRBM resulting from joining the RBMs.
   nvisible2 = 7
   nhidden2 = 4
   vv = rand(nsamples, nvisible + nvisible2)
   hhpartitioned = rand(nsamples, nhidden + nhidden2)
   rbm2 = BMTest.randrbm(nvisible2, nhidden2)
   prbm = BMs.PartitionedRBM{BMs.BernoulliRBM}([rbm; rbm2])
   joinedrbm = BMs.joinrbms([rbm; rbm2])
   BMs.hiddenpotential!(hhpartitioned, joinedrbm, vv)
   @test isapprox(sum(abs.(hhpartitioned - BMs.hiddenpotential(joinedrbm, vv))), 0)

   # Test activation potential of visible nodes for PartitionedRBM
   BMs.visiblepotential!(vv, prbm, hhpartitioned)
   @test isapprox(sum(abs.(vv - BMs.visiblepotential(joinedrbm, hhpartitioned))), 0)
end

function rbmexactloglikelihoodvsbaserate(x::Matrix{Float64}, nhidden::Int)
   a = BMs.logit.(vec(mean(x, dims = 1)))
   nvisible = length(a)
   rbm = BMs.BernoulliRBM(zeros(nvisible, nhidden), a, ones(nhidden))
   baserate = BMs.bernoulliloglikelihoodbaserate(x)
   exactloglik = BMs.loglikelihood(rbm, x, BMs.exactlogpartitionfunction(rbm))
   baserate - exactloglik
end

function bgrbmexactloglikelihoodvsbaserate(x::Matrix{Float64}, nhidden::Int)
   a = BMs.logit.(vec(mean(x, dims = 1)))
   nvisible = length(a)
   bgrbm = BMs.BernoulliGaussianRBM(zeros(nvisible, nhidden), a, ones(nhidden))
   baserate = BMs.bernoulliloglikelihoodbaserate(x)
   exactloglik = BMs.loglikelihood(bgrbm, x, BMs.exactlogpartitionfunction(bgrbm))
   baserate - exactloglik
end

function gbrbmexactloglikelihoodvsbaserate(x::Matrix{Float64}, nhidden::Int)
   a = vec(mean(x, dims = 1))
   nvisible = length(a)
   sd = vec(std(x, dims = 1))
   gbrbm = BMs.GaussianBernoulliRBM(zeros(nvisible, nhidden), a, ones(nhidden), sd)
   baserate = BMs.gaussianloglikelihoodbaserate(x)
   exactloglik = BMs.loglikelihood(gbrbm, x, BMs.exactlogpartitionfunction(gbrbm))
   baserate - exactloglik
end

function test_likelihoodbaserates()
   @test isapprox(BMs.bernoulliloglikelihoodbaserate(10),
         BMs.bernoulliloglikelihoodbaserate(rand([0.0 1.0], 10000, 10)),
         atol = 0.002)

   # Test exact loglikelihood of a baserate RBM
   # vs the calculated loglikelihoodbaserate
   x = rand([0.0 0.0 0.0 1.0 1.0], 100, 10);
   @test isapprox(BMTest.bgrbmexactloglikelihoodvsbaserate(x, 10), 0, atol = 1e-10)
   @test isapprox(BMTest.rbmexactloglikelihoodvsbaserate(x, 10), 0, atol = 1e-10)
   x = rand(100, 10) + randn(100, 10);
   @test isapprox(BMTest.gbrbmexactloglikelihoodvsbaserate(x, 10), 0, atol = 1e-10)

   @test isapprox(BMs.bernoulliloglikelihoodbaserate(10),
         BMs.bernoulliloglikelihoodbaserate(rand([0.0 1.0], 10000, 10)),
         atol = 0.002)

   # Test exact loglikelihood of a baserate RBM
   # vs the calculated loglikelihoodbaserate
   x = rand([0.0 0.0 0.0 1.0 1.0], 100, 10);
   @test isapprox(BMTest.bgrbmexactloglikelihoodvsbaserate(x, 10), 0, atol = 1e-10)
   @test isapprox(BMTest.rbmexactloglikelihoodvsbaserate(x, 10), 0, atol = 1e-10)
   x = rand(100, 10) + randn(100, 10);
   @test isapprox(BMTest.gbrbmexactloglikelihoodvsbaserate(x, 10), 0, atol = 1e-10)

   nvisible = 10
   nhidden = 10
   bgrbm = BMs.BernoulliGaussianRBM(zeros(nvisible, nhidden), rand(nvisible), ones(nhidden));
   @test isapprox(BMs.logpartitionfunction(bgrbm, 0.0),
         BMs.exactlogpartitionfunction(bgrbm))
end

function test_stackrbms_preparetrainlayers()
   x = Matrix{Float64}(undef, 22, 10);
   epochs = 17
   learningrate = 0.007
   nhiddens = [7;6;5]
   learningrate1 = 0.006
   learningrate2 = 0.017
   epochs1 = 27
   epochs2 = 17
   batchsize = 5
   batchsize1 = 1
   batchsize2 = 4
   optimizer = BMs.loglikelihoodoptimizer(learningrate = learningrate)
   optimizer1 = BMs.loglikelihoodoptimizer(learningrate = learningrate1)
   optimizer2 = BMs.loglikelihoodoptimizer(learningrate = learningrate2)

   # no layer-specific instructions
   trainlayers = Vector{BMs.TrainLayer}()
   trainlayers = BMs.stackrbms_preparetrainlayers(trainlayers, x, epochs,
         learningrate, nhiddens, batchsize, optimizer)
   @test length(trainlayers) == length(nhiddens)
   @test all(map(t -> t.learningrate == learningrate, trainlayers))
   @test all(map(t -> t.epochs == epochs, trainlayers))
   @test map(t-> t.nhidden, trainlayers) == nhiddens
   @test all(map(t-> t.batchsize == batchsize, trainlayers))
   @test all(map(t-> t.optimizer.learningrate == learningrate, trainlayers))

   # layer-specific instructions without partitioning
   trainlayers = [
         BMs.TrainLayer(learningrate = learningrate1, epochs = epochs1,
               batchsize = batchsize1, nhidden = nhiddens[1],
               optimizer = optimizer1);
         BMs.TrainLayer(learningrate = learningrate2, epochs = epochs2,
               batchsize = batchsize2, nhidden = nhiddens[2],
               optimizer = optimizer2);
         BMs.TrainLayer(nhidden = nhiddens[3])
   ]
   trainlayers = BMs.stackrbms_preparetrainlayers(trainlayers, x, epochs,
         learningrate, Vector{Int}(), batchsize, optimizer)
   @test length(trainlayers) == length(nhiddens)
   @test trainlayers[1].learningrate == learningrate1
   @test trainlayers[2].learningrate == learningrate2
   @test trainlayers[3].learningrate == learningrate
   @test trainlayers[1].epochs == epochs1
   @test trainlayers[2].epochs == epochs2
   @test trainlayers[3].epochs == epochs
   @test map(t-> t.nhidden, trainlayers) == nhiddens
   @test trainlayers[1].batchsize == batchsize1
   @test trainlayers[2].batchsize == batchsize2
   @test trainlayers[3].batchsize == batchsize
   @test trainlayers[1].optimizer.learningrate == learningrate1
   @test trainlayers[2].optimizer.learningrate == learningrate2
   @test trainlayers[3].optimizer.learningrate == learningrate

   # partitioning
   trainlayers = [
         BMs.TrainPartitionedLayer([
            BMs.TrainLayer(learningrate = learningrate1, epochs = epochs1,
                  rbmtype = BMs.GaussianBernoulliRBM, nhidden = 3,
                  nvisible = 3);
            BMs.TrainLayer(nhidden = 4, nvisible = 7)
         ]);
         BMs.TrainLayer(nhidden = 6);
         BMs.TrainLayer(nhidden = 5)
   ]
   trainlayers = BMs.stackrbms_preparetrainlayers(trainlayers, x, epochs,
         learningrate, Vector{Int}(), batchsize, BMs.NoOptimizer())
   @test length(trainlayers) == length(nhiddens)
   @test trainlayers[1].parts[1].learningrate == learningrate1
   @test trainlayers[1].parts[2].learningrate == learningrate
   @test trainlayers[3].learningrate == learningrate
   @test trainlayers[1].parts[1].epochs == epochs1
   @test trainlayers[1].parts[2].epochs == epochs
   @test trainlayers[2].epochs == epochs
   @test trainlayers[3].epochs == epochs
   @test [trainlayers[1].parts[1].nhidden + trainlayers[1].parts[2].nhidden;
         trainlayers[2].nhidden; trainlayers[3].nhidden] == nhiddens

   # incorrect partitioning must not be allowed
   trainlayers[1].parts[1].nvisible = 20
   @test_throws ErrorException BMs.stackrbms_preparetrainlayers(trainlayers,
         x, epochs, learningrate, Vector{Int}(), batchsize, BMs.NoOptimizer())

   nothing
end

function test_summing_out()
   # Test exact computation of log partition function of DBM:
   # Compare inefficient simple implementation with more efficient one that
   # sums out layers analytically.
   nhiddensvecs = Array[[5;4;6;3], [5;6;7], [4;4;2;4;3], [4;3;2;2;3;2]];
   for i = 1:length(nhiddensvecs)
      dbm = BMTest.randdbm(nhiddensvecs[i]);
      @test isapprox(BMs.exactlogpartitionfunction(dbm),
            BMTest.exactlogpartitionfunctionwithoutsummingout(dbm))
   end
end

"""
    testaisvsexact(bm, percentalloweddiff)
Tests whether the exact log partition function is approximated by the value
estimated by AIS for the given Boltzmann Machine.
The test is successful, if the difference in percent of the log of the exact value
between the exact log partition function and AIS-estimated one is less than
`percentalloweddiff`.
"""
function checkaisvsexact(bm::BMs.AbstractBM, percentalloweddiff::Float64;
      ntemperatures = 100, nparticles = 100)

   exact = BMs.exactlogpartitionfunction(bm)
   estimated = BMs.logpartitionfunction(bm;
         ntemperatures = ntemperatures, nparticles = nparticles)

   @check abs((exact - estimated)/exact) < percentalloweddiff / 100
end

function exactlogpartitionfunctionwithoutsummingout(dbm::BMs.BasicDBM)
   nlayers = length(dbm) + 1
   u = BMs.initcombination(dbm)
   z = 0.0
   while true
      z += exp(-BMs.energy(dbm, u))
      BMs.next!(u) || break
   end
   log(z)
end

function testsummingoutforexactloglikelihood(nunits::Vector{Int})
   x = BMTest.createsamples(1000, nunits[1]);
   dbm = BMs.stackrbms(x, nhiddens = nunits[2:end], batchsize = 5,
         epochs = 50, predbm = true, learningrate = 0.001);
   dbm = BMs.traindbm!(dbm, x,
         learningrates = [0.02*ones(10); 0.01*ones(10); 0.001*ones(10)],
         epochs = 30);
   logz = BMs.exactlogpartitionfunction(dbm)
   @test isapprox(BMs.exactloglikelihood(dbm, x, logz),
         BMTest.exactloglikelihoodwithoutsummingout(dbm, x, logz))
end

function exactloglikelihoodwithoutsummingout(dbm::BMs.BasicDBM, x::Array{Float64,2},
      logz = BMs.exactlogpartitionfunction(dbm))

   nsamples = size(x,1)
   nlayers = length(dbm) + 1

   u = BMs.initcombination(dbm)
   logp = 0.0
   for j = 1:nsamples
      u[1] = vec(x[j,:])

      p = 0.0
      while true
         p += exp(-BMs.energy(dbm, u))

         # next combination of hidden nodes' activations
         BMs.next!(u[2:end]) || break
      end

      logp += log(p)
   end

   logp /= nsamples
   logp -= logz
   logp
end

function test_dbmjoining()
   dbm1 = BMTest.randdbm([5;4;3])
   dbm2 = BMTest.randdbm([4;5;2])
   dbm3 = BMTest.randdbm([6;7;8])

   exactlogpartitionfunction1 = BMs.exactlogpartitionfunction(dbm1)
   exactlogpartitionfunction2 = BMs.exactlogpartitionfunction(dbm2)
   exactlogpartitionfunction3 = BMs.exactlogpartitionfunction(dbm3)

   jointdbm1 = BMs.joindbms(BMs.BasicDBM[dbm1, dbm2])
   # Test use of visibleindexes
   indexes = randperm(15)

   jointdbm2 = BMs.joindbms(BMs.BasicDBM[dbm1, dbm2, dbm3],
            [indexes[1:5], indexes[6:9], indexes[10:15]])

   @test isapprox(
               exactlogpartitionfunction1 + exactlogpartitionfunction2,
               BMs.exactlogpartitionfunction(jointdbm1))

   @test isapprox(
            exactlogpartitionfunction1 + exactlogpartitionfunction2 +
                  exactlogpartitionfunction3,
            BMs.exactlogpartitionfunction(jointdbm2))
end

"""
Tests whether the empirical loglikelihood from samples that are generated from
the exact distribution of Boltzmann machines is equal to the exact loglikelihood.
"""
function check_exactsampling()
   nsamples = 200
   x = BMs.barsandstripes(nsamples, 9)[:, 1:6]

   dbm1 = BMs.fitdbm(x; epochspretraining = 20, epochs = 15, nhiddens = [8; 5])
   BMTest.check_exactsampling(dbm1, x)

   dbm1 = BMs.fitdbm(x; epochspretraining = 20, epochs = 15, nhiddens = [6; 5; 4])
   BMTest.check_exactsampling(dbm1, x)

   rbm = BMs.fitrbm(x; epochs = 100, nhidden = 20)
   BMTest.check_exactsampling(rbm, x)
end

function check_exactsampling(bm::BMs.AbstractBM, x::Matrix{Float64};
      percentalloweddiff = 2.0)

   mbdist = BMs.MultivariateBernoulliDistribution(bm)
   exactsamples = BMs.samples(mbdist, 35000)

   emploglikexactsamples = BMs.empiricalloglikelihood(x, exactsamples)
   exactloglik = BMs.exactloglikelihood(bm, x)
   @check abs((emploglikexactsamples - exactloglik) / exactloglik) <
         percentalloweddiff / 100
end

"""
Fits a larger RBM and a DBM and tests whether the AIS estimation
yields approximately the same results, if performed parallel or not.
"""
function check_likelihoodconsistency()
   nsamples = 1000
   nvariables = 64
   x = BMs.barsandstripes(nsamples, nvariables)
   dbm = BMs.fitdbm(x; epochspretraining = 20, epochs = 15,
         nhiddens = [50; 25; 20])

   logp1 = BMs.loglikelihood(dbm, x[1:10,:];
         parallelized = false, nparticles = 300, ntemperatures = 200)
   logp2 = BMs.loglikelihood(dbm, x[1:10,:];
         parallelized = true, nparticles = 300, ntemperatures = 200)
   @test isapprox(logp1, logp2, rtol = 0.055)

   rbm1 = BMs.fitrbm(x; epochs = 20)
   logz1 = BMs.logpartitionfunction(rbm1; parallelized = false)
   logz1_2 = BMs.logpartitionfunction(rbm1; parallelized = true)
   @test isapprox(logz1, logz1_2, rtol = 0.02)

   logp1 = BMs.loglikelihood(rbm1, x, logz1)
   rbm2 = BMs.fitrbm(x; epochs = 20, nhidden = 30)
   logz2 = BMs.logpartitionfunction(rbm2)
   logp2 = BMs.loglikelihood(rbm2, x, logz2)

   samplemeanrbm =
         BMs.BernoulliRBM(fill(0.0, nvariables, 1), vec(mean(x, dims = 1)), fill(0.0, 1))
   samplemeanrbm2 =
         BMs.BernoulliRBM(fill(0.0, size(rbm1.weights)), vec(mean(x, dims = 1)),
               fill(0.0, length(rbm1.hidbias)))

   # Annealing between two RBMs with zero weights
   @check isapprox(
         BMs.logmeanexp(BMs.aislogimpweights(samplemeanrbm, samplemeanrbm2)),
         (BMs.logpartitionfunctionzeroweights(samplemeanrbm2) -
               BMs.logpartitionfunctionzeroweights(samplemeanrbm)))

   # Annealing from RBM to RBM with sample mean and zero weights
   @check abs(logz1 - BMs.logpartitionfunctionzeroweights(samplemeanrbm) -
         BMs.logmeanexp(BMs.aislogimpweights(samplemeanrbm, rbm1))) <
               logz1 * 0.013

   # Annealing from one RBM to another RBM: Compare partition functions ...
   logimpweights = BMs.aislogimpweights(rbm1, rbm2;
         ntemperatures = 1000, nparticles = 200, burnin = 10)
   @check abs(BMs.logmeanexp(logimpweights) - (logz2 - logz1)) <
         max(logz2, logz1) * 0.03
   # ... and loglikelihood
   lldiff2 = BMs.loglikelihooddiff(rbm1, rbm2, x, logimpweights)
   lldiff1 = logp1 - logp2
   # very weak test, TODO: check and improve
   @check sign(lldiff1) == sign(lldiff2)
end

"""
Tests the likelihood estimation approaches for GaussianBernoulliRBMs
and alternative GaussianBernoulliRBMs.
"""
function check_likelihoodconsistency_gaussian(gbrbmtype::Type{GBRBM}
      ) where GBRBM <: Union{BMs.GaussianBernoulliRBM, BMs.GaussianBernoulliRBM2}

   x = BMs.curvebundles(nperbundle = 300, nvariables = 10,
         nbundles = 4)

   gbrbm1 = BMs.fitrbm(x; nhidden = 30, rbmtype = gbrbmtype,
         cdsteps = 10,
         pcd = false,
         learningrate = 0.001, epochs = 50,
         sdlearningrate = 0.00001)

   datainitrbm = GBRBM(zeros(size(x,2), 1),
         vec(mean(x, dims = 1)), [0.0], vec(std(x, dims = 1)))

   # Compare estimated difference of log partition functions to difference
   # between estimation and exact log partition function of GBRBM with zero weights
   logz1 = BMs.logpartitionfunction(gbrbm1)
   @check abs(logz1 - BMs.logpartitionfunctionzeroweights(datainitrbm) -
         BMs.logmeanexp(BMs.aislogimpweights(datainitrbm, gbrbm1,
               nparticles = 200, ntemperatures = 200))) < 0.03 * logz1

   gbrbm2 = BMs.fitrbm(x; nhidden = 15, rbmtype = gbrbmtype,
         cdsteps = 10,
         learningrate = 0.001, epochs = 100,
         sdlearningrates = [fill(0.0, 30); fill(0.00001, 300)])

   logz2 = BMs.logpartitionfunction(gbrbm2)

   logdiffimpweights = BMs.aislogimpweights(gbrbm1, gbrbm2)
   logzdiff = BMs.logmeanexp(logdiffimpweights)

   # Annealing from one RBM to another RBM: Compare partition functions ...
   @check abs(logzdiff - (logz2 - logz1)) < max(logz2, logz1) * 0.03
   # ... and loglikelihood
   logp1 = BMs.loglikelihood(gbrbm1, x, logz1)
   logp2 = BMs.loglikelihood(gbrbm2, x, logz2)
   lldiff2 = BMs.loglikelihooddiff(gbrbm1, gbrbm2, x, logdiffimpweights)
   lldiff1 = logp1 - logp2
   # very weak test, TODO: check and improve
   @check sign(lldiff1) == sign(lldiff2)
end

"
Tests whether the log-likelihood of Binomial2BernoulliRBMs is approximately
equal to the empirical loglikelihood of data generated by the models, and whether
the partition function estimated by AIS is near the exact value.
"
function check_b2brbm()
   x = BMTest.createsamples(100, 4) + BMTest.createsamples(100, 4)
   b2brbm = BMs.fitrbm(x, rbmtype = BMs.Binomial2BernoulliRBM, epochs = 150,
         nhidden = 4, learningrate = 0.001)
   @check checklikelihoodempirically(b2brbm, x; percentalloweddiff = 1.0)
end

function check_rbm()
   x = BMTest.createsamples(100, 4)
   rbm = BMs.fitrbm(x, epochs = 30,
         nhidden = 4, learningrate = 0.001)
   @check checklikelihoodempirically(rbm, x)
end

function check_gbrbm(gbrbmtype::Type{GBRBM}
      ) where GBRBM <:Union{BMs.GaussianBernoulliRBM, BMs.GaussianBernoulliRBM2}

   x = BMs.curvebundles(nperbundle = 300, nvariables = 10,
         nbundles = 4)
   gbrbm = BMs.fitrbm(x; nhidden = 15, rbmtype = gbrbmtype,
         cdsteps = 10,
         pcd = false,
         learningrate = 0.001, epochs = 50,
         sdlearningrate = 0.00001)

   logz = BMs.logpartitionfunction(gbrbm; parallelized = true,
         ntemperatures = 500, nparticles = 500)
   estloglik = BMs.loglikelihood(gbrbm, x, logz)
   exactloglik = BMs.exactloglikelihood(gbrbm, x)
   @check abs((exactloglik - estloglik) / exactloglik) < 2.5 / 100
end

function test_rbm_monitoring(gbrbmtype::Type{GBRBM}
      ) where GBRBM <:Union{BMs.GaussianBernoulliRBM, BMs.GaussianBernoulliRBM2}

   x = BMs.curvebundles(nperbundle = 20, nvariables = 5,
      nbundles = 4)

   # Test some more options and monitoring
   monitor = BMs.Monitor()
   datadict = BMs.DataDict("Small subset" => x[1:3, :])
   gbrbm = BMs.fitrbm(x; nhidden = 5, rbmtype = gbrbmtype,
         learningrate = 0.000001, epochs = 20,
         batchsize = 10,
         pcd = false,
         sdlearningrate = 0.000000001,
         monitoring = (rbm, epoch) -> begin
            if epoch == 10 || epoch == 20
               BoltzmannMachines.monitorexactloglikelihood!(monitor, rbm, epoch, datadict)
               BoltzmannMachines.monitorreconstructionerror!(monitor, rbm, epoch, datadict)
            end
         end)
   nothing
end

function check_mdbm_rbm_b2brbm()
   nsamples = 100
   nvariables = 4

   x1 = BMTest.createsamples(nsamples, nvariables) +
         BMTest.createsamples(nsamples, nvariables);
   x2 = BMTest.createsamples(nsamples, nvariables);
   x = hcat(x1, x2);
   trainlayers1 = [
         BMs.TrainPartitionedLayer([
            BMs.TrainLayer(rbmtype = BMs.Binomial2BernoulliRBM,
                  nvisible = nvariables, nhidden = nvariables - 1,
                  learningrate = 0.0015);
            BMs.TrainLayer(rbmtype = BMs.BernoulliRBM,
                  nvisible = nvariables, nhidden = nvariables)
         ]);
         BMs.TrainLayer(nhidden = nvariables);
         BMs.TrainLayer(nhidden = 3)]

   dbm1 = BMs.fitdbm(x, epochs = 15,
         epochspretraining = 15,
         learningratepretraining = 0.001,
         pretraining = trainlayers1)

   @check BMTest.checklikelihoodempirically(dbm1, x; percentalloweddiff = 5.5,
         ntemperatures = 300, nparticles = 300)

   # partitioned second layer
   trainlayers2 = [
         trainlayers1[1],
         BMs.TrainPartitionedLayer([
            BMs.TrainLayer(nhidden = 3);
            BMs.TrainLayer(nhidden = 2)
         ]),
         BMs.TrainLayer(nhidden = 4),
         BMs.TrainLayer(nhidden = 3)
   ]

   dbm2 = BMs.fitdbm(x, pretraining = trainlayers2)
   @check BMTest.checklikelihoodempirically(dbm2, x; percentalloweddiff = 6.3,
         ntemperatures = 300, nparticles = 300)
end

function checklikelihoodempirically(rbm::BMs.AbstractRBM, x::Matrix{Float64};
      percentalloweddiff = 0.5, ntemperatures::Int = 100, nparticles::Int = 100)

   logz = BMs.logpartitionfunction(rbm; parallelized = true,
      nparticles = nparticles, ntemperatures = ntemperatures)
   estloglik = BMs.loglikelihood(rbm, x, logz)
   exactloglik = BMs.exactloglikelihood(rbm, x)
   @check abs((exactloglik - estloglik) / exactloglik) < percentalloweddiff / 100

   emploglik = BMs.empiricalloglikelihood(rbm, x, 1000000)
   @check abs((exactloglik - emploglik) / exactloglik) < percentalloweddiff / 100
end

function checklikelihoodempirically(dbm::BMs.MultimodalDBM, x::Matrix{Float64};
      percentalloweddiff = 0.5, ntemperatures::Int = 100, nparticles::Int = 100)

   emploglik = BMs.empiricalloglikelihood(dbm, x, 1000000)
   estloglik = BMs.loglikelihood(dbm, x; parallelized = true,
         ntemperatures = ntemperatures, nparticles = nparticles)
   exactloglik = BMs.exactloglikelihood(dbm, x)
   @check abs((exactloglik - emploglik) / exactloglik) < percentalloweddiff / 100
   @check abs((exactloglik - estloglik) / exactloglik) < percentalloweddiff / 100
end

"""
Test DBMs with Gaussian visible nodes.
"""
function check_mdbm_gaussianvisibles()

   x = convert(Matrix{Float64}, readdlm(
         joinpath(dirname(pathof(BoltzmannMachines)), "..", "test/data/iris.csv"), ',',
         header = true)[1][:,1:4]);

   datadict = BMs.DataDict("x" => x)

   monitor1 = BMs.Monitor()
   trainlayers = [
         BMs.TrainLayer(rbmtype = BMs.GaussianBernoulliRBM,
               nhidden = 4,
               sdlearningrate = 0.000001,
               monitoring = (rbm, epoch) -> begin
                  BMs.monitorexactloglikelihood!(monitor1, rbm, epoch, datadict)
               end);
         BMs.TrainLayer(nhidden = 4);
         BMs.TrainLayer(nhidden = 4)]

   learningrates = [0.02*ones(10); 0.01*ones(10); 0.001*ones(10)]

   seed = round(Int, rand()*typemax(Int), RoundDown)
   Random.seed!(seed)
   dbm1 = BMs.stackrbms(x, epochs = 20, predbm = true, learningrate = 0.001,
         trainlayers = trainlayers)
   # BoltzmannMachinesPlots.plotevaluation(monitor1, BMs.monitorexactloglikelihood)
   dbm1 = BMs.traindbm!(dbm1, x,
         learningrates = learningrates,
         epochs = 30);

   Random.seed!(seed)
   dbm2 = BMs.fitdbm(x, epochs = 30,
         epochspretraining = 20,
         learningratepretraining = 0.001,
         pretraining = trainlayers,
         learningrates = learningrates)

   # first and second dbm must be equal
   @test length(dbm1) == length(dbm2)
   @test isapprox(dbm1[1].sd, dbm2[1].sd)
   for i in 1:length(dbm1)
      @test isapprox(dbm1[i].weights, dbm2[i].weights)
      @test isapprox(dbm1[i].visbias, dbm2[i].visbias)
      @test isapprox(dbm1[i].hidbias, dbm2[i].hidbias)
   end

   # Test AIS
   @check BMTest.checkaisvsexact(dbm1, 0.6)

   # Test exact likelihood vs estimated likelihood
   estloglik = BMs.loglikelihood(dbm1, x; nparticles = 400, ntemperatures = 200)
   exactloglik = BMs.exactloglikelihood(dbm1, x)
   @check abs((exactloglik - estloglik) / exactloglik) < 3.5 / 100
end

function test_mdbm_architectures()
   nsamples = 5
   data = hcat(float.(rand(nsamples, 3) .< 0.5), randn(nsamples, 3))

   BMs.fitdbm(data; epochs = 2, pretraining = [
         BMs.TrainPartitionedLayer([
               BMs.TrainLayer(nvisible = 3, nhidden = 3);
               BMs.TrainLayer(nvisible = 3, nhidden = 3,
                     rbmtype = BMs.GaussianBernoulliRBM)]);
               BMs.TrainLayer(nhidden = 4);
               BMs.TrainLayer(nhidden = 3)])

   BMs.fitdbm(data; epochs = 2, pretraining = [
         BMs.TrainPartitionedLayer([
               BMs.TrainLayer(nvisible = 3, nhidden = 3);
               BMs.TrainLayer(nvisible = 3, nhidden = 3,
                     rbmtype = BMs.GaussianBernoulliRBM2)]);
         BMs.TrainPartitionedLayer([
               BMs.TrainLayer(nhidden = 2)
               BMs.TrainLayer(nhidden = 2)]);
         BMs.TrainLayer(nhidden = 4)])

   data2 = hcat(data,
         float.(rand(nsamples, 3) .< 0.5) .+ float.(rand(nsamples, 3) .< 0.5))
   BMs.fitdbm(data2; epochs = 2, pretraining = [
         BMs.TrainPartitionedLayer([
               BMs.TrainLayer(nvisible = 3, nhidden = 3);
               BMs.TrainLayer(nvisible = 3, nhidden = 3,
                     rbmtype = BMs.GaussianBernoulliRBM2);
               BMs.TrainLayer(nvisible = 3, nhidden = 3,
                     rbmtype = BMs.Binomial2BernoulliRBM)]);
         BMs.TrainPartitionedLayer([
               BMs.TrainLayer(nvisible = 3, nhidden = 3)
               BMs.TrainLayer(nvisible = 6, nhidden = 3)]);
         BMs.TrainLayer(nhidden = 3)])
   nothing
end

end # of module BMTest


text:
"""
    InwardsPopulationDispersal <: NeighborhoodRule

    InwardsPopulationDispersal(; kw...)
    InwardsPopulationDispersal{R}(; kw...)
    InwardsPopulationDispersal{R,W}(; kw...)

Implements deterministic dispersal from populations in neighboring cells to the current
cell.

The result should be identical to those obtained substituting `InwardsDispersal` for
[`OutwardsDispersal`](@ref) but will perform better when populations are spread across the
grid.

# Keywords

- `neighborhood`: Any DynamicGrids.jl `Neighborhood`, or an
    already constructed [`DispersalKernel`](@ref). Using this keyword means `radius` is
    ignored, and for a `DispersalKernel`, all other keywords are ignored.
- `neighborhood`: `Neighborhood` object specifying the range from the origin of the
    discretised dispersal kernal. Defaults to `Window(radius)`.
- `formulation`: kernel formulation object holding the exact form of the kernal.
    Default [`ExponentialKernel`](@ref).
- `cellsize`: the cell size of the discretised kernal (i.e. simulation grid size).
    Default is 1.0.
- `distancemethod`: [`DistanceMethod`](@ref) object for calculating distance between cells.
    The default is [`CentroidToCentroid`](@ref).

Pass grid `Symbol`s to `R` or both `R` and `W` type parameters to use to specific grids.
"""
struct InwardsDispersal{R,W,N<:AbstractKernelNeighborhood} <: NeighborhoodRule{R,W}
    neighborhood::N
end
function InwardsDispersal{R,W}(; kw...) where {R,W}
    InwardsDispersal{R,W}(DispersalKernel(; kw...))
end

@inline applyrule(data, rule::InwardsDispersal, N, I) = kernelproduct(rule)


text:
g = incidence_list(6)
@test length(nodes(g)) == 6
add_edge!(g, 1, 3, 1)
add_edge!(g, 2, 4, 2)
add_edge!(g, 3, 4, 4)
es = edges(g)
@test es[1] == (1,3,1)
@test es[2] == (2,4,2)
@test es[3] == (3,4,4)
@test num_timestamps(g) == 4