text:
"""
    Saturation function for quadratic saturation models for machines
        Se(x) = B * (x - A)^2 / x
"""
function saturation_function(
    machine::Union{PSY.RoundRotorQuadratic, PSY.SalientPoleQuadratic},
    x::Number,
)
    Sat_A, Sat_B = PSY.get_saturation_coeffs(machine)
    return Sat_B * (x - Sat_A)^2 / x
end

"""
    Saturation function for exponential saturation models for machines
        Se(x) = B * x^A
"""
function saturation_function(
    machine::Union{PSY.RoundRotorExponential, PSY.SalientPoleExponential},
    x::Number,
)
    Sat_A, Sat_B = PSY.get_saturation_coeffs(machine)
    return Sat_B * x^Sat_A
end

function rectifier_function(I::Float64)
    if I <= 0.0
        return 1.0
    elseif I <= 0.433
        return 1.0 - 0.577 * I
    elseif I < 0.75
        return sqrt(0.75 - I^2)
    elseif I <= 1.0
        return 1.732 * (1.0 - I)
    else
        return 0.0
    end
end

function saturation_function(avr::PSY.ESAC1A, x::Number)
    Sat_A, Sat_B = PSY.get_saturation_coeffs(avr)
    return Sat_B * (x - Sat_A)^2 / x
end

function rectifier_function(I::T) where {T <: Number}
    if I <= 0.0
        return one(T)
    elseif I <= 0.433
        return 1.0 - 0.577 * I
    elseif I < 0.75
        return sqrt(0.75 - I^2)
    elseif I <= 1.0
        return 1.732 * (1.0 - I)
    else
        return zero(T)
    end
end


text:
immutable BSONOID
    _wrap_::Ptr{Void}
    _ref_::Any

    BSONOID() = begin
        buffer = Array(UInt8, 12)
        ccall(
            (:bson_oid_init, libbson),
            Void, (Ptr{UInt8}, Ptr{Void}),
            buffer,
            C_NULL
            )
            r = Compat.unsafe_convert(Ptr{UInt8}, buffer)
        new(r, buffer)
    end

    BSONOID(str::AbstractString) = begin
        cstr = string(str)

        isValid = ccall(
            (:bson_oid_is_valid, libbson),
            Bool, (Ptr{UInt8}, Csize_t),
            cstr,
            length(cstr)
            )
        isValid || error("'" * str * "': not a valid BSONOID string")

        buffer = Array(UInt8, 12)
        ccall(
            (:bson_oid_init_from_string, libbson),
            Void, (Ptr{UInt8}, Ptr{UInt8}),
            buffer,
            cstr
            )
            r = Compat.unsafe_convert(Ptr{UInt8}, buffer)
        new(r, buffer)
    end

    BSONOID(_ref_::Any) = new(pointer(_ref_), _ref_)
end
export BSONOID

import Base.==
==(lhs::BSONOID, rhs::BSONOID) = ccall(
    (:bson_oid_equal, libbson),
    Bool, (Ptr{Void}, Ptr{Void}),
    lhs._wrap_, rhs._wrap_
    )
export ==

hash(oid::BSONOID, h::UInt) = hash(
    ccall(
        (:bson_oid_hash, libbson),
        UInt32, (Ptr{UInt8},),
        oid._wrap_
        ),
    h
    )
export hash

function convert(::Type{AbstractString}, oid::BSONOID)
    cstr = Array(UInt8, 25)
    ccall(
        (:bson_oid_to_string, libbson),
        Void, (Ptr{UInt8}, Ptr{UInt8}),
        oid._wrap_,
        cstr
        )
    return String(unsafe_string(Compat.unsafe_convert(Ptr{UInt8}, cstr)))
end
export convert

string(oid::BSONOID) = convert(AbstractString, oid)
export string

show(io::IO, oid::BSONOID) = print(io, "BSONOID($(convert(AbstractString, oid)))")
export show

ext: jl


text:
using RobustAndOptimalControl
using ControlSystems

G = ssrand(3,4,5)

G = partition(G, 1, 1)

@test G isa ExtendedStateSpace
@test G == G
@test G ≈ G
@test ControlSystems.noutputs(G) == 3
@test ControlSystems.ninputs(G) == 4
@test ndims(G) == 2
@test size(G) == (3,4)
@test eltype(G) == typeof(G)
@test ControlSystems.numeric_type(G) == Float64
@test_throws ErrorException G[1]

@show partition(ssrand(2,2,1), 1, 1)

@test noise_mapping(G) == ss(G)[2:end, 1]


text:
### PHYSICAL FUNCTIONS ###
include("RK4_functions.jl")

get_elastic_stiffness_tensor(G,λ) = SymmetricTensor{4, 3}( (i,j,k,l) -> Dᵉ_func(i,j,k,l,G,λ))
get_elastic_stiffness_tensor(e::Elasticity) = SymmetricTensor{4, 3}( (i,j,k,l) -> Dᵉ_func(i,j,k,l,e.G,e.λ))
get_elastic_compliance_tensor(e::Elasticity) = SymmetricTensor{4, 3}( (i,j,k,l) -> Cᵉ_func(i,j,k,l,e.G,e.λ))

get_elastic_stiffness_tensor(r::Rheology) = get_elastic_stiffness_tensor(r.elasticity)
get_elastic_compliance_tensor(r::Rheology) = get_elastic_compliance_tensor(r.elasticity)
# Second invariant of the deviatoric stress tensor :
get_τ(s ; plas = :DruckerPrager) = sqrt(0.5 * s ⊡ s)
get_τ(s,::DruckerPrager) = sqrt(0.5 * s ⊡ s)
get_τ(s,::Damage) = sqrt(0.5 * s ⊡ s)
get_τ(s,::VonMises) = sqrt(3/2 * s ⊡ s)

get_τ(s,r::Rheology) = get_τ(s::AbstractTensor,r.plasticity)

### Damage functions ###

free_energy_convexity(r,D,A₁,B₁) = 1/compute_Γ(r,A₁,B₁) > 0 ? true : false

# eq 16 Bhat2012 & 2016 & notes (because in Bhat2011 c2 isn't the same form as in Harsha's notes) :
#compute_c1(d::Damage,D) = sqrt(1-cos(d.ψ)^2)/(π*cos(d.ψ)^(3/2)*((D/d.D₀)^(1/3) - 1 + d.β/cos(d.ψ))^(3/2))
function compute_c1(d::Damage,D)
    α = cosd(d.ψ)
    @assert α > 0
    if d.D₀>0
        @assert (D/d.D₀) >= 1
        sqrt(1-α^2)/(π*α^(3/2)*((D/d.D₀)^(1/3) - 1 + d.β/α)^(3/2))
    elseif (d.D₀==0) & (D==0)
        return sqrt(1-α^2)/(π*α^(3/2)*(d.β/α)^(3/2)) # lim D->0 D/D0 = 1
      end
end
# Perol&Bhat2016 : 1/α  or  Harsha's notes : 1/α^2 ???
function compute_c2(d::Damage,D)
    α = cosd(d.ψ)
    (d.D₀==0) && (return 0.0)
    return (sqrt(1 - α^2)/α^2) * (d.D₀^(2/3)/(1 - D^(2/3))) ###
  end

function compute_c3(d::Damage,D)
    α = cosd(d.ψ)
    @assert α > 0
    (d.D₀ == D) && (return 0.0)
    @assert (D/d.D₀) >= 1
    (2sqrt(α)/π)*((D/d.D₀)^(1/3) - 1)^(1/2)
end

compute_c1(r::Rheology,D) = compute_c1(r.damage::Damage,D)
compute_c2(r::Rheology,D) = compute_c2(r.damage::Damage,D)
compute_c3(r::Rheology,D) = compute_c3(r.damage::Damage,D)

function compute_c1c2c3(d::Damage,D)
    c1 = compute_c1(d,D)
    c2 = compute_c2(d,D)
    c3 = compute_c3(d,D)
    return c1, c2, c3
end

compute_c1c2c3(r::Rheology,D) = compute_c1c2c3(r.damage::Damage,D)

# eq 15 Bhat2012 (A1 : *c2*c3), Perol&Bhat2016 (A1 : ...*c2)*c3):
# Perol&Bhat2016 is the corrected version, and the one implemented
compute_A(r::Rheology,c1,c2,c3) = r.damage.μ*c1 + (1.0 + r.damage.μ*c2)*c3
compute_A(d::Damage,c1,c2,c3) = d.μ*c1 + (1.0 + d.μ*c2)*c3
compute_B(c1,c2,c3) = c1 + c2*c3

function compute_AB(d::Damage,c1,c2,c3)
    A = compute_A(d,c1,c2,c3)
    B = compute_B(c1,c2,c3)
    return A, B
end
compute_AB(r::Rheology,c1,c2,c3) = compute_AB(r.damage,c1,c2,c3)
compute_AB(r::Rheology,D) = compute_AB(r,compute_c1c2c3(r,D)...)

# eq 11 in Harsha's notes :
compute_A1(r::Rheology,A) = A * sqrt((π*r.damage.D₀*(1 - r.elasticity.ν))/cosd(r.damage.ψ)^3)
compute_B1(r::Rheology,B) = B * sqrt((π*r.damage.D₀*(1 - r.elasticity.ν))/cosd(r.damage.ψ)^3)

function compute_A1B1(r::Rheology,A,B)
    A1 = compute_A1(r,A)
    B1 = compute_B1(r,B)
    return A1, B1
end

compute_A1B1(r::Rheology,D) = compute_A1B1(r,compute_AB(r,D)...)

function compute_a1(B1,Γ)
    return (1/Γ)*(1 + B1^2/2)
end

function compute_b1(A1,B1,Γ)
    return -(1/Γ)*((A1*B1)/2)
end

function compute_b2(r::Rheology,A1,Γ)
    return (1/Γ)*(A1^2/2 + (3*(1-2r.elasticity.ν))/(2*(1+r.elasticity.ν)))
end

function compute_Γ(r::Rheology,A₁,B₁)
    ν = r.elasticity.ν
    return (3*(1-2ν))/(2*(1+ν)) + (3*(1-2ν)*B₁^2)/(4*(1+ν)) + A₁^2/2
end

function compute_Γ(r::Rheology,D)
    c1 = compute_c1(r,D)
    c2 = compute_c2(r,D)
    c3 = compute_c3(r,D)
    A = compute_A(r,c1,c2,c3)
    B = compute_B(c1,c2,c3)
    A₁ = compute_A1(r,A)
    B₁ = compute_B1(r,B)
    ν = r.elasticity.ν
    return (3*(1-2ν))/(2*(1+ν)) + (3*(1-2ν)*B₁^2)/(4*(1+ν)) + A₁^2/2
end

function compute_dc1dD(r::Rheology,D)
    d = r.damage
    α = cosd(d.ψ)
    term1 = (-sqrt(1-α^2))/(2*π*α^(3/2)*D^(2/3)*d.D₀^(1/3))
    term2 = ((D/d.D₀)^(1/3) - 1 + (d.β/α))^(-5/2)
    return term1 * term2
end

function compute_dc2dD(r::Rheology,D)
    d = r.damage
    α = cosd(d.ψ)
    term1 = (2*sqrt(1-α^2)*d.D₀^(2/3))/(3*α^2*D^(1/3))
    term2 = (1 - D^(2/3))^(-2)
    return term1 * term2
end

function compute_dc3dD(r::Rheology,D)
    d = r.damage
    α = cosd(d.ψ)
    term1 = (sqrt(α))/(3*π*D^(2/3)*d.D₀^(1/3))
    term2 = ((D/d.D₀)^(1/3) - 1)^(-1/2)

    # ensure term2 isn't Infinity
    (term2 == Inf) && (term2 = 1e15)
    return term1 * term2
end

function compute_dA1dD(r::Rheology,dc1dD,dc2dD,dc3dD,c2,c3)
    d = r.damage
    α = cosd(d.ψ)
    term1 = sqrt((π*d.D₀*(1 - r.elasticity.ν))/(α^3))
    term2 = d.μ*dc1dD + dc3dD + d.μ*c2*dc3dD + d.μ*c3*dc2dD
    return term1 * term2
end

function compute_dB1dD(r::Rheology,dc1dD,dc2dD,dc3dD,c2,c3)
    d = r.damage
    α = cosd(d.ψ)
    term1 = sqrt((π*d.D₀*(1 - r.elasticity.ν))/(α^3))
    term2 = dc1dD + c2*dc3dD + c3*dc2dD
    return term1 * term2
end

function compute_dΓdD(r::Rheology,A1,B1,dA1dD,dB1dD)
    return ((3*(1-2r.elasticity.ν))/(2*(1+r.elasticity.ν)))*B1*dB1dD + A1*dA1dD
end

function compute_dA1dt(r::Rheology,dc1dt,dc2dt,dc3dt,c2,c3)
    d = r.damage
    α = cosd(d.ψ)
    term1 = sqrt((π*d.D₀*(1 - r.elasticity.ν))/(α^3))
    term2 = d.μ*dc1dt + dc3dt + d.μ*c2*dc3dt + d.μ*c3*dc2dt
    return term1 * term2
end

function compute_dB1dt(r::Rheology,dc1dt,dc2dt,dc3dt,c2,c3)
    d = r.damage
    α = cosd(d.ψ)
    term1 = sqrt((π*d.D₀*(1 - r.elasticity.ν))/(α^3))
    term2 = dc1dt + c2*dc3dt + c3*dc2dt
    return term1 * term2
end

function compute_dΓdt(r::Rheology,A1,B1,dA1dt,dB1dt)
    return ((3*(1-2r.elasticity.ν))/(2*(1+r.elasticity.ν)))*B1*dB1dt + A1*dA1dt
end

# function compute_dΓdt(p::RockParams,A1,B1,dA1dt,dB1dt)
#     return ((3*(1-2p.ν))/(2*(1+p.ν)))*B1*dB1dt + A1*dA1dt
# end

function compute_da1dt(B1,Γ,dB1dt,dΓdt)
    return -(dΓdt/Γ^2)*(1 + B1^2/2) + (B1*dB1dt)/Γ
end

function compute_db1dt(A1,B1,Γ,dA1dt,dB1dt,dΓdt)
    return (dΓdt/Γ^2)*((A1*B1)/2) - (1/2Γ)*(dA1dt*B1 + A1*dB1dt)
end

function compute_db2dt(r::Rheology,A1,Γ,dA1dt,dΓdt)
    return -(dΓdt/Γ^2)*(A1^2/2 + (3*(1-2r.elasticity.ν))/(2*(1+r.elasticity.ν))) + (2/Γ)*A1*dA1dt
end

function compute_dσdt(r::Rheology,a1,b1,da1dt,db1dt,ϵ,γ,dϵdt,dγdt)
    return r.elasticity.G * (da1dt*ϵ + a1*dϵdt + db1dt*γ + b1*dγdt)
end

function compute_dτdt(r::Rheology,b1,b2,db1dt,db2dt,ϵ,γ,dϵdt,dγdt)
    return r.elasticity.G * (db1dt*ϵ + b1*dϵdt + db2dt*γ + b2*dγdt)
end

# Damage growth and KI

function compute_KI(r::Rheology,σ,τ,A,B)
    return (A*σ + B*τ) * sqrt(π*r.damage.a)
end

function compute_KI(r::Rheology,σ,τ,D)
    c1, c2, c3 = compute_c1c2c3(r,D)
    A, B = compute_AB(r,c1,c2,c3)
    # println("c1 : ",c1)
    # println("c2 : ",c2)
    # println("c3 : ",c3)
    # println("σ : ",σ)
    # println("τ : ",τ)
    return (A*σ + B*τ) * sqrt(π*r.damage.a)
end

function compute_KI(d::Damage,σij::SymmetricTensor,D)
    c1, c2, c3 = compute_c1c2c3(d,D)
    A, B = compute_AB(d,c1,c2,c3)
    p = 1/3 * tr(σij) # trial pressure, negative in compression
    sij = dev(σij) # trial deviatoric stress
    τ = get_τ(sij,d)
    return (A*p + B*τ) * sqrt(π*d.a)
end
compute_KI(r::Rheology,σij,D) = compute_KI(r.damage,σij,D)

function compute_dDdl(r::Rheology,D)
    d = r.damage
    return (3*D^(2/3)*d.D₀^(1/3))/(cosd(d.ψ)*d.a)
end

function compute_subcrit_damage_rate(r::Rheology, KI, D)
    ((KI <= 0) | (D >= 1)) && (return 0.0)
    d = r.damage
    e = r.elasticity
    ρ = 2700 ##### TODO better
    Vs = sqrt(e.G/ρ)
    Vr = Vs * (0.862 + 1.14e.ν)/(1 + e.ν)

    dDdl = compute_dDdl(r,D) # damage derivative wrt crack length
    dldt = min(d.l̇₀*(KI/d.K₁c)^(d.n),Vr)  # cracks growth rate
    @assert dDdl * dldt >= 0
    return dDdl * dldt
end

function compute_subcrit_damage_rate(r::Rheology, σij::SymmetricTensor, D)
    KI = compute_KI(r,σij,D)
    ((KI <= 0) | (D >= 1)) && (return 0.0)
    d = r.damage
    e = r.elasticity
    ρ = 2700 ##### TODO better
    Vs = sqrt(e.G/ρ)
    Vr = Vs * (0.862 + 1.14e.ν)/(1 + e.ν)

    dDdl = compute_dDdl(r,D) # damage derivative wrt crack length
    dldt = min(d.l̇₀*(KI/d.K₁c)^(d.n),Vr)  # cracks growth rate
    @assert dDdl * dldt >= 0
    return dDdl * dldt
end

# dDdt_thresh=1e-10
# KIC_max=1e8
# ΔKI = d.K₁c*(dDdt_thresh/(1.5e3*0.24))^(1/d.n)
# KI = -1e6:100:3e6
# x = 2.5.*KI./ΔKI
# KIC = (1.0 .- tanh.(x)).*(KIC_max - d.K₁c).*0.5 .+ d.K₁c

function compute_subcrit_damage_rate_variable_KIC(r::Rheology, KI, D ; dDdt_thresh=1e-40, KIC_max=1e15 )
    #((KI <= 0) | (D >= 1)) && (return 0.0)
    d = r.damage
    e = r.elasticity
    ρ = 2700 ##### TODO better
    Vs = sqrt(e.G/ρ)
    Vr = Vs * (0.862 + 1.14e.ν)/(1 + e.ν)

    dDdl = compute_dDdl(r,D) # damage derivative wrt crack length

    # get KI length scale of transition from high KIC to real KIC
    #println("KI : ",KI)
    KIC = compute_KIC(r, KI, dDdt_thresh, KIC_max)

    dldt = min(d.l̇₀*(abs(KI)/KIC)^(d.n),Vr)  # cracks growth rate

    try @assert dDdl * dldt >= 0
    catch e
        println("KI : ", KI)
        println("KIC : ", KIC)
        throw(e)
    end

    return dDdl * dldt
end

function compute_KIC(r, KI, dDdt_thresh, KIC_max)
    d = r.damage
    ΔKI = d.K₁c*(dDdt_thresh/(1.5e3*d.l̇₀))^(1/d.n)
    x = 2.5*KI/ΔKI
    #println("x : ", x)
    #println("tanh(x) : ", tanh(x))
    KIC = (1.0 - tanh(x))*(KIC_max - d.K₁c)*0.5 + d.K₁c
    return KIC
end

function compute_subcrit_damage_rate(r::Rheology, σ, τ, D)
    KI = compute_KI(r,σ,τ,D)
    # return zero if KI is negative
    ((KI <= 0) | (D >= 1)) && (return 0.0)

    d = r.damage
    e = r.elasticity
    ρ = 2700 ##### TODO better
    Vs = sqrt(e.G/ρ)
    Vr = Vs * (0.862 + 1.14e.ν)/(1 + e.ν)

    dDdl = compute_dDdl(r,D) # damage derivative wrt crack length
    dldt = min(d.l̇₀*(KI/d.K₁c)^(d.n),Vr)  # cracks growth rate  speed
    # println("d.l̇₀ = ",d.l̇₀)
    # println("d.K₁c = ",d.K₁c)
    # println("d.n = ",d.n)
    # println("Vr = ",Vr)
    # println("dldt = ",dldt)
    @assert dDdl >= 0

    return dDdl * dldt
end

function get_damage_constrained_Δt(model,u,ΔD_max)
    nu,nD = getnbasefunctions.(model.cellvalues_tuple)
    nqp = getnquadpoints(model.cellvalues_tuple[2])
    Δt_max = 1e9
    dDdt_max = 0.0
    for cell in CellIterator(model.dofhandler)
        cellid = cell.current_cellid.x
        states = model.material_state[cellid]
        r = model.material_properties[cellid]

        cell_dofs = celldofs(cell)
        ue = u[cell_dofs]
        De = exp.(ue[nu+1:end])
        for qp in nqp
            state = states[qp]
            D = function_value(model.cellvalues_tuple[2],qp,De)
            KI = compute_KI(r,state.temp_σ,D)
            dDdt = compute_subcrit_damage_rate(r, KI, D)
            Δt_max = min(ΔD_max/dDdt,Δt_max)
            dDdt_max = max(dDdt,dDdt_max)
        end
    end
    #println("dDdt_max based on converged σ and D = ",dDdt_max)
    return Δt_max
end

# function get_damage_constrained_Δt(model,ΔD_max)
#     nqp = getnquadpoints(model.cellvalues_tuple[1])
#     KI = 0.0
#     for cell in CellIterator(model.dofhandler)
#         cellid = cell.current_cellid.x
#         r = model.material_properties[cellid]
#         for qp in nqp
#             state = model.material_state[cellid][qp]
#             KI = max(KI,compute_KI(r,state.temp_σ,D))
#         end
#     end
#     KI_max =
#     D = function_value(model.cellvalues_tuple[2],qp,De)
#     KI = compute_KI(r,state.temp_σ,D)
#     dDdt = compute_subcrit_damage_rate(r, KI, D)
#     Δt_max = min(ΔD_max/dDdt,Δt_max)
#     dDdt_max = max(dDdt,dDdt_max)
#     println("dDdt_max based on converged σ and D = ",dDdt_max)
#     return Δt_max
# end
# function compute_subcrit_damage_rate(r::Rheology, σ::T, D) where {T<:AbstractArray}
#     d = r.damage
#     e = r.elasticity
#     G = e.E / (2*(1 + e.ν))
#     Vs = sqrt(G/2700) ##### TODO better
#
#     KI = compute_KI(r,σ,D)
#     dDdl = compute_dDdl(r,D) # damage derivative wrt crack length
#     dldt = min(d.l̇₀*(KI/d.K₁c)^(d.n),Vs)  # cracks growth rate #TODO should be rayleigh wave speed
#     return dDdl * dldt
# end

# function compute_dσijdt(r,ϵij, D, dϵijdt, dDdt)
#
#     # strain decomposition
#     ϵ = tr(ϵij)
#     e = dev(ϵij)
#     γ = sqrt(2.0 * e ⊡ e)
#     dϵdt = tr(dϵijdt)
#     dedt = dev(dϵijdt)
#     dγdt = sqrt(2.0 * dedt ⊡ dedt)
#
#     # damage constants
#     c1, c2, c3 = compute_c1c2c3(r,D)
#     A, B = compute_AB(r,c1,c2,c3)
#     A1 = compute_A1(r,A)
#     B1 = compute_B1(r,B)
#     Γ = compute_Γ(r,A1,B1)
#
#     # derivatives
#     dc1dD = compute_dc1dD(r,D)
#     dc2dD = compute_dc2dD(r,D)
#     dc3dD = compute_dc3dD(r,D)
#     dA1dD = compute_dA1dD(r,dc1dD,dc2dD,dc3dD,c2,c3)
#     dB1dD = compute_dB1dD(r,dc1dD,dc2dD,dc3dD,c2,c3)
#     dΓdD = compute_dΓdD(r,A₁,B₁,dA1dD,dB1dD)
#
#     dA1dt = dA1dD * dDdt
#     dB1dt = dB1dD * dDdt
#     dΓdt = dΓdD * dDdt
#     # unpack
#     G = r.elasticity.G
#     ν = r.elasticity.ν
#
#     # Identity matrix
#     Id = SymmetricTensor{2,3}(δ)
#     # dσijdt, with σij = F*H
#     F = G / Γ
#     H = ( (3*(1-2ν))/(1+ν) + A1^2 - A1*B1*ϵ/γ )*ϵij +
#         ( (3ν/(1+ν) + B1^2/2 - A1^2/3 + A1*B1*ϵ/(3γ))*ϵ - A1*B1*γ/2 ) * Id
#
#     dFdt = -G/Γ^2 * dΓdt
#     dHdt = ( 2A1*dA1dt - (dA1dt*B1 + A1*dB1dt)*ϵ/γ - A1*B1*(dϵdt/γ - ϵ*dγdt/γ^2) )*ϵij +
#            ( (3*(1-2ν))/(1+ν) + A1^2 - A1*B1*ϵ/γ )*dϵijdt +
#            ( (B1*dB1dt - (2/3)*A1*dA1dt + (1/3)*((dA1dt*B1 + A1*dB1dt)*ϵ/γ - A1*B1*(dϵdt/γ - ϵ*dγdt/γ^2)))*ϵ +
#              (3ν/(1+ν) + B1^2/2 - A1^2/3 + A1*B1*ϵ/(3γ))*dϵdt -
#              0.5*((dA1dt*B1 + A1*dB1dt)*γ + A1*B1*dγdt) ) * Id
#
#     return dFdt*H + F*dHdt
# end

# function state_tensors_system!(du,u,p,t)
#     D, ϵ, σ = u
#     rheology, Δt, Δϵ = p
#     KI = compute_KI(rheology,σ,D)
#     du[1] = dD = (KI > 0) ? compute_subcrit_damage_rate(rheology, σ, D) : 0.0
#     du[2] = dϵ = Δϵ/Δt
#     du[3] = dσ = compute_dσijdt(rheology, ϵ, D, du[2], du[1])
# end
function compute_ϵ_oop(r::Rheology,A1,B1,σᵢⱼ::SymmetricTensor)

    σoop = σᵢⱼ[3,3]
    σ = 1/3*tr(σᵢⱼ)
    sᵢⱼ = dev(σᵢⱼ)
    τ = get_τ(sᵢⱼ,r.damage)
    G = r.elasticity.G
    ν = r.elasticity.ν

    # stress tensor calculation
    term1 = ( 1 + A1*B1*σ/(2*τ) + B1^2/2) * σoop
    term2 = (3ν/(1+ν) + A1*B1*σ/(2*τ) - A1^2/3 + B1^2/2) * σ
    term3 = A1*B1*τ/3
    return (1/(2G)) * (term1 - term2 + term3)
end
compute_ϵ_oop(r,D,σᵢⱼ)=compute_ϵ_oop(r,compute_A1B1(r,D)...,σᵢⱼ)

function compute_ϵ_oop(r::Rheology,A1,B1,σᵢⱼ::SymmetricTensor)

    σoop = σᵢⱼ[3,3]
    σ = 1/3*tr(σᵢⱼ)
    sᵢⱼ = dev(σᵢⱼ)
    τ = get_τ(sᵢⱼ,r.damage)
    G = r.elasticity.G
    ν = r.elasticity.ν

    # stress tensor calculation
    term1 = ( 1 + A1*B1*σ/(2*τ) + B1^2/2) * σoop
    term2 = (3ν/(1+ν) + A1*B1*σ/(2*τ) - A1^2/3 + B1^2/2) * σ
    term3 = A1*B1*τ/3
    return (1/(2G)) * (term1 - term2 + term3)
end

function compute_ϵij(r,A1,B1,σᵢⱼ::SymmetricTensor)

    if all(σᵢⱼ.==0)
        return σᵢⱼ
    end
    σ = 1/3*tr(σᵢⱼ)
    sᵢⱼ = dev(σᵢⱼ)
    τ = get_τ(sᵢⱼ,r.damage)
    G = r.elasticity.G
    ν = r.elasticity.ν
    Id = SymmetricTensor{2,3}(δ)
    # strain invariants
    term1 = ( 1 + A1*B1*σ/(2*τ) + B1^2/2) * σᵢⱼ
    term2 = (3ν/(1+ν) + A1*B1*σ/(2*τ) - A1^2/3 + B1^2/2) * σ
    term3 = A1*B1*τ/3
    return (1/(2G)) * (term1 + (term3 - term2)*Id)
end

compute_ϵij(r::Rheology,D,σᵢⱼ::SymmetricTensor) = (D==0) ? compute_ϵij(r,0.0,0.0,σᵢⱼ) : compute_ϵij(r,compute_A1B1(r,D)...,σᵢⱼ)

function compute_σij(r,A1,B1,Γ,ϵij)
    # TODO make a visco elastic version of this function

    if all(x->x==0,ϵij)
        return SymmetricTensor{2,3}(ϵij)
    end

    G = r.elasticity.G
    ν = r.elasticity.ν
    Id = SymmetricTensor{2,3}(δ)
    # strain invariants
    ϵ = tr(ϵij)
    e = dev(ϵij)
    γ = sqrt(2.0 * e ⊡ e)
    (γ == 0) && (γ = nextfloat(0.0))
    (γ == Inf) && (γ = prevfloat(Inf))

    # stress tensor calculation
    term1 = ( (3*(1-2ν))/(1+ν) + A1^2 - A1*B1*ϵ/γ ) * ϵij
    term2 = (3ν/(1+ν) + B1^2/2 - A1^2/3 + A1*B1*ϵ/(3γ)) * ϵ
    term3 = -A1*B1*γ/2
    return (G/Γ) * (term1 + (term2 + term3)*Id)
end

# CHECKED WITH ALL USED FUNCTIONS
function compute_σij(r,D,ϵij)
    # TODO make a visco elastic version of this function

    if all(x->x==0,ϵij)
        return SymmetricTensor{2,3}(ϵij)
    end
    # unpack
    G = r.elasticity.G
    ν = r.elasticity.ν

    # Create Identity 3*3 tensor
    Id = SymmetricTensor{2,3}(δ)

    # Damage constants
    c1, c2, c3 = compute_c1c2c3(r,D)
    A, B = compute_AB(r,c1,c2,c3)
    A1 = compute_A1(r,A)
    B1 = compute_B1(r,B)
    Γ = compute_Γ(r,A1,B1)

    # strain invariants
    ϵ = tr(ϵij)
    e = dev(ϵij)
    γ = sqrt(2.0 * e ⊡ e)
    (γ == 0) && (γ = nextfloat(0.0))
    (γ == Inf) && (γ = prevfloat(Inf))
    # stress tensor calculation
    term1 = ( (3*(1-2ν))/(1+ν) + A1^2 - A1*B1*ϵ/γ ) * ϵij
    term2 = (3ν/(1+ν) + B1^2/2 - A1^2/3 + A1*B1*ϵ/(3γ)) * ϵ
    term3 = -A1*B1*γ/2
    return (G/Γ) * (term1 + (term2 + term3)*Id)
end

function compute_σ̇ij(r,D,Ḋ,ϵij,ϵ̇ij)

    if all(x->x==0,ϵij)
        all(x->x==0,ϵ̇ij) && (return SymmetricTensor{2,3}(ϵij))
    end
    # unpack
    G = r.elasticity.G
    ν = r.elasticity.ν

    # Create Identity 3*3 tensor
    𝕀 = SymmetricTensor{2,3}(δ)

    # Damage constants
    c1, c2, c3 = compute_c1c2c3(r,D)
    A, B = compute_AB(r,c1,c2,c3)
    A1 = compute_A1(r,A)
    B1 = compute_B1(r,B)
    Γ = compute_Γ(r,A1,B1)
    dc1dD = compute_dc1dD(r,D)
    dc2dD = compute_dc2dD(r,D)
    dc3dD = compute_dc3dD(r,D)
    dA1dD = compute_dA1dD(r,dc1dD,dc2dD,dc3dD,c2,c3)
    dA1dD = compute_dB1dD(r,dc1dD,dc2dD,dc3dD,c2,c3)
    dΓdD = compute_dΓdD(r,A1,B1,dA1dD,dB1dD)

    # strain invariants
    ϵ = tr(ϵij)
    e = dev(ϵij)
    γ = sqrt(2.0 * e ⊡ e)
    (γ == 0) && (γ = nextfloat(0.0))
    (γ == Inf) && (γ = prevfloat(Inf))

    ϵ̇ = tr(ϵ̇ij)
    ė = dev(ϵ̇ij)
    γ̇ = e ⊡ ė / (2*γ)

    # stress tensor calculation
    f1 = f₁(ϵ,γ,A1,B1)
    f2 = f₂(ϵ,γ,A1,B1)
    f3 = f₃(A1,B1)
    df1dt = df₁dt(ϵ, γ, A1, B1, ϵ̇, γ̇, Ḋ, dA1dD, dB1dD)
    df2dt = df₂dt(ϵ, γ, A1, B1, ϵ̇, γ̇, Ḋ, dA1dD, dB1dD)
    df3dt = df₃dt(D,Ḋ,dA1dD, dB1dD)

    term1 = df1dt*ϵij + f1*ϵ̇ij + (df2dt*ϵ + f2*ϵ̇ + df3dt*γ + f3*γ̇)*𝕀
    term2 = f1*ϵij + (f2*ϵ + f3*γ)*𝕀
    return (G/Γ) * (term1 - (1/Γ)*dΓdD*Ḋ*term2)
end
f₁(ϵ,γ,A1,B1) = 3*(1-2ν)/(1+ν) + A1^2 - A1*B1*ϵ/γ
f₂(ϵ,γ,A1,B1) = 3*ν/(1+ν) + B1^2/2 - A1^2/3 + A1*B1*ϵ/(3*γ)
f₃(A1,B1) = -A1*B1/2

df₁dt(ϵ, γ, A1, B1, ϵ̇, γ̇, Ḋ, dA1dD, dB1dD) = -(A1*B1/γ)*ϵ̇ + (A1*B1*ϵ/γ^2)*γ̇ + (2*A1*dA1dD - (dA1dD*B1 + A1*dB1dD)*ϵ/γ)*Ḋ
df₂dt(ϵ, γ, A1, B1, ϵ̇, γ̇, Ḋ, dA1dD, dB1dD) = (A1*B1/(3*γ))*ϵ̇ - (A1*B1*ϵ/(3*γ^2))*γ̇ + (B1*dB1dD - 2/3*A1*dA1dD - (dA1dD*B1 + A1*dB1dD)*ϵ/(3γ))*Ḋ
df₃dt(D,Ḋ,dA1dD, dB1dD) = -0.5*(dA1dD*B1 + A1*dB1dD)*Ḋ

function state_system!(du,u,p,t)
    D, ϵ, γ, σ, τ = u
    r, dϵdt, dγdt = p

    # damage constants
    (D < r.damage.D₀) && (println("D = ",D); @warn "D < D0, something went wrong")
    (D == r.damage.D₀) && (D += 1e-9) # insure D > D0 to prevent singularity
    isnan(D) && println("D is NaN")

    c1, c2, c3 = compute_c1c2c3(r,D)
    A, B = compute_AB(r,c1,c2,c3)

    # TODO : Check KI sign to avoid unnecessary calculations

    A1 = compute_A1(r,A)
    B1 = compute_B1(r,B)
    Γ = compute_Γ(r,A1,B1)
    a1 = compute_a1(B1,Γ)
    b1 = compute_b1(A1,B1,Γ)
    b2 = compute_b2(r,A1,Γ)

    # derivatives
    # D
    dDdt = compute_subcrit_damage_rate(r, σ, τ, D)
    #println("D = ", D)
    #println("dDdt = ", dDdt, "\n")

    # c1, c2, c3
    dc1dt = compute_dc1dD(r,D) * dDdt
    dc2dt = compute_dc2dD(r,D) * dDdt
    dc3dt = compute_dc3dD(r,D) * dDdt

    # println("dDdt : ", dDdt)
    # println("dc1dt : ",dc1dt)
    # println("dc2dt : ",dc2dt)
    # println("dc3dt : ",dc3dt)
    # A1, B1
    dA1dt = compute_dA1dt(r,dc1dt,dc2dt,dc3dt,c2,c3)
    dB1dt = compute_dB1dt(r,dc1dt,dc2dt,dc3dt,c2,c3)

    # Γ
    dΓdt = compute_dΓdt(r,A1,B1,dA1dt,dB1dt)

    # a1, b1, b2
    da1dt = compute_da1dt(B1,Γ,dB1dt,dΓdt)
    db1dt = compute_db1dt(A1,B1,Γ,dA1dt,dB1dt,dΓdt)
    db2dt = compute_db2dt(r,A1,Γ,dA1dt,dΓdt)

    # println("Γ : ",Γ)
    # println("A1 : ",A1)
    # println("B1 : ",B1)
    # println("dΓdt : ",dΓdt) # issue here with Nan
    # println("dA1dt : ",dA1dt) # issue here with Nan
    # println("dB1dt : ",dB1dt) # issue here with Nan
    # println("a1,b1,da1dt,db1dt,ϵ,γ,dϵdt,dγdt : ")
    # println(a1)
    # println(b1)
    # println(da1dt) # issue here with Nan
    # println(db1dt) # issue here with Nan
    # println(ϵ)
    # println(γ)
    # println(dϵdt)
    # println(dγdt)
    du[1] = dD = dDdt
    du[2] = dϵ = dϵdt
    du[3] = dγ = dγdt
    du[4] = dσ = compute_dσdt(r,a1,b1,da1dt,db1dt,ϵ,γ,dϵdt,dγdt)
    du[5] = dτ = compute_dτdt(r,b1,b2,db1dt,db2dt,ϵ,γ,dϵdt,dγdt)

    # assertions
    @assert dD >= 0
end

#### TEST ####

# elasticity = Elasticity(E = 70e9,
#                                     ν = 0.3)
# plasticity = DruckerPrager(μ = 0.6,
#                                       C = 1e6)
# damage = BRSDamage( μ = 0.6, # Friction coef
#                     β = 0.1, # Correction factor
#                     K₁c = 1.74e6, # Critical stress intensity factor (Pa.m^(1/2))
#                     a = 1e-3, # Initial flaw size (m)
#                     ψ = atan(0.7), # crack angle to the principal stress (radians)
#                     D₀ = 0.3, #0.1 Initial flaw density
#                     n = 34.0, # Stress corrosion index
#                     l̇₀ = 0.24, # Ref. crack growth rate (m/s)
#                     H = 50e3, # Activation enthalpy (J/mol)
#                     A = 5.71 ) # Preexponential factor (m/s)
# r = Rheology(damage,nothing,elasticity,nothing)
#
# D = range(r.damage.D₀,1.0,length = 50)
# c1 = compute_c1.(Ref(r),D)
# c2 = compute_c2.(Ref(r),D)
# c3 = compute_c3.(Ref(r),D)
# A = compute_A.(Ref(r),c1,c2,c3)
# B = compute_B.(c1,c2,c3)
#
# A1 = compute_A1.(Ref(r),A)
# B1 = compute_B1.(Ref(r),B)
#Γ = compute_Γ.(Ref(r),A1,B1)
#Γ = compute_Γ.(Ref(r),D)
# convex_cond = 1 ./ Γ
#
# p = lineplot(D,convex_cond)
# lineplot!(p,D[1:end-1],A1[1:end-1])
# lineplot!(p,D[1:end-1],Γ[1:end-1])
# show(p)
# pyplot()
#
# σ = -10e6
# τ = 55e6
# D_vec = [0.1:0.001:0.999;]
# KI_vec = compute_KI.(Ref(r), Ref(σ), Ref(τ), D_vec)
# dDdt_vec = compute_subcrit_damage_rate.(Ref(r), KI_vec, D_vec)
# p = plot(D_vec,dDdt_vec)
#     ylims!(0,1)
# display(p)
# p2 = plot(D_vec,KI_vec./r.damage.K₁c)
#     ylims!(0,10)
# display(p2)


text:
using Test
using Printf: @sprintf
using LinearAlgebra: dot
using StaticArrays: SVector
using QuantumLattices.Essentials.Terms
using QuantumLattices.Essentials.Spatials: Point,PID,Bond,Bonds,Lattice,decompose,acrossbonds,zerothbonds
using QuantumLattices.Essentials.DegreesOfFreedom: IDFConfig,Table,IID,Index,Internal,FilteredAttributes,OID,Operator,Operators,coordabsent,coordpresent,Boundary
using QuantumLattices.Interfaces: rank,update!,expand!,id,kind,reset!
using QuantumLattices.Prerequisites: Float,decimaltostr
using QuantumLattices.Prerequisites.CompositeStructures: NamedContainer
using QuantumLattices.Mathematics.AlgebraOverFields: ID,SimpleID,idtype
import QuantumLattices.Interfaces: dimension,expand
import QuantumLattices.Essentials.DegreesOfFreedom: isHermitian,otype
import QuantumLattices.Essentials.Terms: couplingcenter,couplingcenters,abbr

struct TID <: IID nambu::Int end
Base.adjoint(sl::TID)=TID(3-sl.nambu)

struct TIndex{S} <: Index{PID{S},TID}
    scope::S
    site::Int
    nambu::Int
end
Base.fieldnames(::Type{<:TIndex})=(:scope,:site,:nambu)
Base.union(::Type{P},::Type{I}) where {P<:PID,I<:TID}=TIndex{fieldtype(P,:scope)}

struct TFock <: Internal{TID} end
dimension(f::TFock)=2
Base.getindex(f::TFock,i::Int)=TID(i)

struct TCID <: SimpleID
    center::Int
    nambu::Int
end
Base.fieldnames(::Type{<:TCID})=(:center,:nambu)

struct TCoupling{V<:Number,I<:ID{TCID}} <: Coupling{V,I}
    value::V
    id::I
end
Base.repr(tc::TCoupling)=@sprintf "%s ph(%s)@(%s)" decimaltostr(tc.value) join(tc.id.nambus,',') join(tc.id.centers,',')
function expand(tc::TCoupling,pids::NTuple{R,PID},focks::NTuple{R,TFock},kind::Union{Val{K},Nothing}=nothing) where {R,K}
    nambus=tc.id.nambus
    pids=NTuple{rank(tc),eltype(pids)}(pids[tc.id.centers[i]] for i=1:rank(tc))
    return ((tc.value,NTuple{rank(tc),TIndex{fieldtype(pids|>eltype,:scope)}}(TIndex(pids[i].scope,pids[i].site,nambus[i]) for i=1:rank(tc))),)
end

struct TOperator{V<:Number,I<:ID{OID}} <: Operator{V,I}
    value::V
    id::I
end

abbr(::Type{<:Term{ST,:TermMu} where ST})=:tmu
isHermitian(::Type{<:Term{ST,:TermMu} where ST})=true
otype(T::Type{<:Term{ST,:TermMu} where ST},I::Type{<:OID})=TOperator{T|>valtype,ID{I,T|>rank}}

abbr(::Type{<:Term{ST,:TermHopping} where ST})=:thp
isHermitian(::Type{<:Term{ST,:TermHopping} where ST})=false
otype(T::Type{<:Term{ST,:TermHopping} where ST},I::Type{<:OID})=TOperator{T|>valtype,ID{I,T|>rank}}

@testset "Subscript" begin
    sub=@subscript (x1,4,4,x2) with x1<x2
    @test sub==deepcopy(sub)
    @test isequal(sub,deepcopy(sub))
    @test rank(sub)==rank(typeof(sub))==2
    @test dimension(sub)==dimension(typeof(sub))==4
    @test string(sub)=="(x1,4,4,x2) with $(sub.identifier)"
    @test sub(Val('M'),1,2)==(1,4,4,2)
    @test sub(Val('C'),1,2)==true
    @test sub(Val('C'),2,1)==false

    sub=@subscript (x1,x2,x1,x2)
    @test sub(Val('M'),1,2)==(1,2,1,2)
    @test sub(Val('C'),1,2)==true

    sub=Subscript{4}()
    @test rank(sub)==1 && dimension(sub)==4
    @test string(sub)=="(*,*,*,*)"
    @test sub(Val('M'),2)==(2,2,2,2)
    @test sub(Val('C'),2)==true

    sub=Subscript((1,2,2,1))
    @test rank(sub)==0 && dimension(sub)==4
    @test string(sub)=="(1,2,2,1)"
    @test sub(Val('M'))==(1,2,2,1)
    @test sub(Val('C'))==true
end

@testset "Subscripts" begin
    sub1=@subscript (x1,2) with x1<2
    sub2=@subscript (y1,y2,y1,y2) with y1<y2
    subs=Subscripts(sub1,sub2)
    @test rank(subs)==rank(typeof(subs))==3
    @test rank(subs,1)==rank(typeof(subs),1)==1
    @test rank(subs,2)==rank(typeof(subs),2)==2
    @test dimension(subs)==dimension(typeof(subs))==6
    @test dimension(subs,1)==dimension(typeof(subs),1)==2
    @test dimension(subs,2)==dimension(typeof(subs),2)==4
    @test subs(Val('M'),(1,2,3))==(1,2,2,3,2,3)
    @test subs(Val('C'),(1,2,3))==true
    @test subs(Val('C'),(2,2,3))==false
    @test subs(Val('C'),(1,3,2))==false
    @test subs(Val('C'),(2,3,2))==false

    sub3=Subscript((6,6))
    @test subs==sub1*sub2
    @test subs*sub3==Subscripts(sub1,sub2,sub3)
    @test sub3*subs==Subscripts(sub3,sub1,sub2)
    @test subs*Subscripts(sub3)==Subscripts(sub1,sub2,sub3)
    @test expand(subs*sub3,(3,3,2,3,2,3,8,9))|>collect==[(1,2,1,2,1,2,6,6),(1,2,1,3,1,3,6,6),(1,2,2,3,2,3,6,6)]
end

@testset "centerrelated" begin
    @test couplingcenter(Coupling,2,2,Val(1))==1
    @test couplingcenter(Coupling,1,2,Val(1))==1
    @test couplingcenters(Coupling,('*','*'),Val(1))==(1,1)
    @test couplingcenters(Coupling,(1,2,3),Val(3))==(1,2,3)
    @test couplingcenters("(1,2,3,4)")==(1,2,3,4)
end

@testset "couplings" begin
    @test @couplings(TCoupling(1.0,ID(TCID(1,1),TCID(2,1))))==Couplings(TCoupling(1.0,ID(TCID(1,1),TCID(2,1))))
    @test @couplings(TCoupling(1.0,ID(TCID(1,1),TCID(2,1)))+TCoupling(1.0,ID(TCID(1,2),TCID(2,2))))==TCoupling(1.0,ID(TCID(1,1),TCID(2,1)))+TCoupling(1.0,ID(TCID(1,2),TCID(2,2)))
end

@testset "TermFunction" begin
    ta=TermAmplitude()
    @test ta()==1

    ta=TermAmplitude(x->x+3.0)
    @test ta(1.0)==4.0

    tcs=TCoupling(1.0,ID(TCID(1,1),TCID(2,1)))+TCoupling(2.0,ID(TCID(1,2),TCID(2,2)))
    termcouplings=TermCouplings(tcs)
    @test termcouplings==deepcopy(TermCouplings(tcs))
    @test isequal(termcouplings,deepcopy(TermCouplings(tcs)))
    @test termcouplings()==tcs

    tcs1=TCoupling(1.0,ID(TCID(1,1),TCID(2,1)))+TCoupling(1.0,ID(TCID(1,2),TCID(2,2)))
    tcs2=TCoupling(1.0,ID(TCID(1,2),TCID(2,1)))+TCoupling(1.0,ID(TCID(2,1),TCID(1,2)))
    termcouplings=TermCouplings(i->(tcs1,tcs2)[(i-1)%2+1])
    @test termcouplings(1)==tcs1
    @test termcouplings(2)==tcs2

    termmodulate=TermModulate(:t)
    @test termmodulate(t=1)==1
    @test termmodulate(mu=1)==nothing
    @test termmodulate()==nothing

    termmodulate=TermModulate(:t,t->t*2.0)
    @test termmodulate(2)==4
end

@testset "Term" begin
    point=Point(PID(1,1),(0.0,0.0),(0.0,0.0))
    config=IDFConfig{TFock}(pid->TFock(),[PID(1,1)])
    term=Term{'F',:TermMu,2}(:mu,1.5,0,couplings=TCoupling(1.0,ID(TCID(1,2),TCID(1,1))),amplitude=(bond->3),modulate=false)
    @test term|>statistics==term|>typeof|>statistics=='F'
    @test term|>kind==term|>typeof|>kind==:TermMu
    @test term|>id==term|>typeof|>id==:mu
    @test term|>valtype==term|>typeof|>valtype==Float
    @test term|>rank==term|>typeof|>rank==2
    @test term|>abbr==term|>typeof|>abbr==:tmu
    @test term|>isHermitian==term|>typeof|>isHermitian==true
    @test term==deepcopy(term)
    @test isequal(term,deepcopy(term))
    @test string(term)=="TermMu{2F}(id=mu,value=1.5,bondkind=0,factor=1.0)"
    @test repr(term,point,config)=="tmu: 4.5 ph(2,1)@(1,1)"
    @test +term==term
    @test -term==term*(-1)==replace(term,factor=-term.factor)
    @test 2*term==term*2==replace(term,factor=term.factor*2)
    @test term/2==term*(1/2)==replace(term,factor=term.factor/2)

    p1=Point(PID(1,1),(0.0,0.0),(0.0,0.0))
    p2=Point(PID(1,2),(0.0,0.0),(0.0,0.0))
    config=IDFConfig{TFock}(pid->TFock(),[PID(1,1),PID(1,2)])
    tcs1=Couplings(TCoupling(1.0,ID(TCID(1,2),TCID(1,2))))
    tcs2=Couplings(TCoupling(1.0,ID(TCID(1,1),TCID(1,1))))
    term=Term{'F',:TermMu,2}(:mu,1.5,0,couplings=bond->(tcs1,tcs2)[bond.pid.site%2+1],amplitude=bond->3,modulate=true)
    @test repr(term,p1,config)=="tmu: 4.5 ph(1,1)@(1,1)"
    @test repr(term,p2,config)=="tmu: 4.5 ph(2,2)@(1,1)"
    @test one(term)==replace(term,value=1.0)
    @test zero(term)==replace(term,value=0.0)
    @test term.modulate(mu=4.0)==4.0
    @test term.modulate(t=1.0)==nothing
    @test update!(term,mu=4.25)==replace(term,value=4.25)
    @test term.value==4.25
end

@testset "expand" begin
    point=Point(PID(1,1),(0.0,0.0),(0.0,0.0))
    config=IDFConfig{TFock}(pid->TFock(),[PID(1,1)])
    table=Table(config,filter(attr->attr≠:nambu,FilteredAttributes(TIndex)))
    term=Term{'F',:TermMu,2}(:mu,1.5,0,couplings=TCoupling(1.0,ID(TCID(1,2),TCID(1,1))),amplitude=bond->3.0,modulate=true)
    operators=Operators(TOperator(+2.25,(TIndex(1,1,2),TIndex(1,1,1)),seqs=(1,1)))
    @test expand(term,point,config,table,true,coordabsent)==operators
    @test expand(term,point,config,table,false,coordabsent)==operators*2

    bond=Bond(1,Point(PID('b',2),(1.5,1.5),(1.0,1.0)),Point(PID('a',1),(0.5,0.5),(0.0,0.0)))
    config=IDFConfig{TFock}(pid->TFock(),[PID('a',1),PID('b',2)])
    table=Table(config,filter(attr->attr≠:nambu,FilteredAttributes(TIndex)))
    term=Term{'F',:TermHopping,2}(:t,1.5,1,couplings=TCoupling(1.0,ID(TCID(1,2),TCID(2,1))),amplitude=bond->3.0,modulate=true)
    operators=Operators(TOperator(4.5,(TIndex('a',1,2),TIndex('b',2,1)),rcoords=(SVector(0.5,0.5),SVector(1.5,1.5)),icoords=(SVector(0.0,0.0),SVector(1.0,1.0)),seqs=(1,2)))
    @test expand(term,bond,config,table,true)==operators
    @test expand(term,bond,config,table,false)==operators+operators'

    lattice=Lattice("Tuanzi",[Point(PID(1,1),(0.0,0.0),(0.0,0.0))],vectors=[[1.0,0.0]],neighbors=1)
    bonds=Bonds(lattice)
    config=IDFConfig{TFock}(pid->TFock(),lattice.pids)
    table=Table(config,filter(attr->attr≠:nambu,FilteredAttributes(TIndex)))
    term=Term{'F',:TermMu,2}(:mu,1.5,0,couplings=TCoupling(1.0,ID(TCID(1,2),TCID(1,1))),amplitude=bond->3.0,modulate=true)
    operators=Operators(TOperator(+2.25,(TIndex(1,1,2),TIndex(1,1,1)),seqs=(1,1)))
    @test expand(term,bonds,config,table,true,coordabsent)==operators
    @test expand(term,bonds,config,table,false,coordabsent)==operators*2
end

@testset "Parameters" begin
    ps1=Parameters{(:t1,:t2,:U)}(1.0im,1.0,2.0)
    ps2=Parameters{(:t1,:U)}(1.0im,2.0)
    ps3=Parameters{(:t1,:U)}(1.0im,2.1)
    @test match(ps1,ps2)==true
    @test match(ps1,ps3)==false
end

@testset "Generator" begin
    lattice=Lattice("Tuanzi",[Point(PID(1,1),(0.0,0.0),(0.0,0.0)),Point(PID(1,2),(0.5,0.0),(0.0,0.0))],vectors=[[1.0,0.0]],neighbors=1)
    bonds=Bonds(lattice)
    config=IDFConfig{TFock}(pid->TFock(),lattice.pids)
    table=Table(config,filter(attr->attr≠:nambu,FilteredAttributes(TIndex)))
    boundary=Boundary()
    t=Term{'F',:TermHopping,2}(:t,2.0,1,couplings=TCoupling(1.0,ID(TCID(1,2),TCID(2,1))))
    μ=Term{'F',:TermMu,2}(:μ,1.0,0,couplings=TCoupling(1.0,ID(TCID(1,2),TCID(1,1))),modulate=true)
    tops1=expand(t,filter(acrossbonds,bonds,Val(:exclude)),config,table,true,coordpresent)
    tops2=expand(one(t),filter(acrossbonds,bonds,Val(:include)),config,table,true,coordpresent)
    μops=expand(one(μ),filter(zerothbonds,bonds,Val(:include)),config,table,true,coordpresent)

    optp=TOperator{Float,ID{OID{TIndex{Int},SVector{2,Float},SVector{2,Float},Int},2}}
    genops=GenOperators(tops1,NamedContainer{(:μ,)}((μops,)),NamedContainer{(:t,:μ)}((tops2,Operators{optp|>idtype,optp}())))
    @test genops==deepcopy(genops) && isequal(genops,deepcopy(genops))
    @test genops==GenOperators((t,μ),bonds,config,table,true,coordpresent)
    @test genops|>eltype==genops|>typeof|>eltype==optp
    @test genops|>idtype==genops|>typeof|>idtype==optp|>idtype
    @test expand!(Operators{idtype(optp),optp}(),genops,boundary,t=2.0,μ=1.5)==tops1+tops2*2.0+μops*1.5
    @test empty!(deepcopy(genops))==empty(genops)==GenOperators(empty(μops),NamedContainer{(:μ,)}((empty(μops),)),NamedContainer{(:t,:μ)}((empty(μops),empty(μops))))
    @test reset!(deepcopy(genops),(t,μ),bonds,config,table,true,coordpresent)==genops

    gen=Generator((t,μ),bonds,config,table,true,boundary)
    @test gen==deepcopy(gen) && isequal(gen,deepcopy(gen))
    @test Parameters(gen)==Parameters{(:t,:μ)}(2.0,1.0)
    @test expand!(Operators{idtype(optp),optp}(),gen)==expand(gen)==tops1+tops2*2.0+μops
    @test expand(gen,:t)==tops1+tops2*2.0
    @test expand(gen,:μ)==μops
    @test expand(gen,1)+expand(gen,2)+expand(gen,3)+expand(gen,4)==expand(gen)
    @test expand(gen,:μ,1)+expand(gen,:μ,2)==μops
    @test expand(gen,:t,3)==tops1
    @test expand(gen,:t,4)==tops2*2.0
    @test empty!(deepcopy(gen))==Generator{coordpresent}((t,μ),empty(bonds),empty(config),empty(table),true,boundary)==empty(gen)
    @test reset!(empty(gen),lattice)==gen
    @test update!(gen,μ=1.5)|>expand==tops1+tops2*2.0+μops*1.5
end


text:
using jInv.Mesh;
using FactoredEikonalFastMarching;
using Printf

#using PyPlot
#close("all")

include("runAccuracyExperiment.jl");
include("getWorkUnit.jl");

"""
A function for running the experiments in the paper:
Eran Treister and Eldad Haber, A fast marching algorithm for the factored eikonal equation, Journal of Computational Physics, 324, 210-225, 2016.
To acctually run the experiments in the paper: increase numOfRefinements
"""
function runPaperExperiments(numOfRefinements)
println("******************** 2D experiments ***************************");

I = [4,8];
n = zeros(Int64,2);

h0 = [0.05,0.05];
WU = zeros(numOfRefinements)
for k = 1:numOfRefinements
	h = (0.5^k)*h0;
	n[1] = round(Int64,(I[1]/h[1])+1);
	n[2] = round(Int64,(I[2]/h[2])+1);
	WU[k] = GetWorkunitForLoop(n,h);
end

println("******************** Gradient Slowness ***************************");
println("h\t\tn\t\terror(1st)[max,l2]\ttime(1st)\terror(2nd)[max,l2]\ttime(2nd)")

h0 = [0.05,0.05];
for k=1:numOfRefinements
	h = (0.5^k)*h0;
	n[1] = round(Int64,(I[1]/h[1])+1);
	n[2] = round(Int64,(I[2]/h[2])+1);
	(kappaSquared,src,T_exact) = getAnalyticalConstGrad2D(n,h);
	runExperimentAndWriteResults(kappaSquared,h,src,n,T_exact,WU[k]);
end

println("******************** Gradient Velocity ***************************");

println("h\t\tn\t\terror(1st)[max,l2]\ttime(1st)\terror(2nd)[max,l2]\ttime(2nd)")
h0 = [0.05,0.05];
for k=1:numOfRefinements
	h = (0.5^k)*h0;
	n[1] = round(Int64,(I[1]/h[1])+1);
	n[2] = round(Int64,(I[2]/h[2])+1);
	(kappaSquared,src,T_exact) = getAnalyticalConstGradInv2D(n,h);
	runExperimentAndWriteResults(kappaSquared,h,src,n,T_exact,WU[k]);
end

println("******************** Gaussian factor ***************************");

println("h\t\tn\t\terror(1st)[max,l2]\ttime(1st)\terror(2nd)[max,l2]\ttime(2nd)")
h0 = [0.05,0.05];
for k=1:numOfRefinements
	h = (0.5^k)*h0;
	n[1] = round(Int64,(I[1]/h[1])+1);
	n[2] = round(Int64,(I[2]/h[2])+1);
	(kappaSquared,src,T_exact) = getSmoothGaussianMedium(n,h);
	runExperimentAndWriteResults(kappaSquared,h,src,n,T_exact,WU[k]);
end

println("*******************************************************************************")
println("*******************************************************************************")
println("************************* 3-D EXPERIMENTS *************************************")
println("*******************************************************************************")
println("*******************************************************************************")

##########################################################################################
##########################################################################################
# 3D experiments:
##########################################################################################

I = [1.6,1.6,0.8];
n = zeros(Int64,3);

h0 = [0.1,0.1,0.1];
WU = zeros(numOfRefinements)
for k = 1:numOfRefinements
	h = (0.5^k)*h0;
	n[1] = round(Int64,(I[1]/h[1])+1);
	n[2] = round(Int64,(I[2]/h[2])+1);
	n[3] = round(Int64,(I[3]/h[3])+1);
	WU[k] = GetWorkunitForLoop(n,h);
end

println("******************** Gradient Slowness ***************************");
println("h\t\tn\t\terror(1st)[max,l2]\ttime(1st)\terror(2nd)[max,l2]\ttime(2nd)")

for k=1:numOfRefinements
	h = (0.5^k)*h0;
	n[1] = round(Int64,(I[1]/h[1])+1);
	n[2] = round(Int64,(I[2]/h[2])+1);
	n[3] = round(Int64,(I[3]/h[3])+1);
	(kappaSquared,src,T_exact) = getAnalyticalConstGrad3D(n,h);
	runExperimentAndWriteResults(kappaSquared,h,src,n,T_exact,WU[k]);
end

println("******************** Gradient Velocity ***************************");
println("h\t\tn\t\terror(1st)[max,l2]\ttime(1st)\terror(2nd)[max,l2]\ttime(2nd)")
for k=1:numOfRefinements
	h = (0.5^k)*h0;
	n[1] = round(Int64,(I[1]/h[1])+1);
	n[2] = round(Int64,(I[2]/h[2])+1);
	n[3] = round(Int64,(I[3]/h[3])+1);
	(kappaSquared,src,T_exact) = getAnalyticalConstGradInv3D(n,h);
	runExperimentAndWriteResults(kappaSquared,h,src,n,T_exact,WU[k]);
end

println("******************** Gaussian factor ***************************");
println("h\t\tn\t\terror(1st)[max,l2]\ttime(1st)\terror(2nd)[max,l2]\ttime(2nd)")
for k=1:numOfRefinements
	h = (0.5^k)*h0;
	n[1] = round(Int64,(I[1]/h[1])+1);
	n[2] = round(Int64,(I[2]/h[2])+1);
	n[3] = round(Int64,(I[3]/h[3])+1);
	(kappaSquared,src,T_exact) = getSmoothGaussianMedium(n,h);
	runExperimentAndWriteResults(kappaSquared,h,src,n,T_exact,WU[k]);
end
return;
end
############################################################################################

runPaperExperiments(2);


text:
function dv_cont(T::Float32, P::Float32)::Float32
    P_atm = P * 1.01325e-5 # Pa -> atm
    return 1e-4 * (0.211/P_atm) * ( (T/273.) ^ 1.94)
end

function dv(T::Float32, r::Float32, P::Float32; accom=


text:
using Term.TermMarkdown
using Markdown

m1 = md"""
This is a test function for math syntax

```math
\\alpha^2 - \\sqrt9 \\in [1, 2, 3] \\cup [4, 5, 6]
```

whut
"""

m2 = md"""
# Markdown rendering in Term.jl
## twp
### three
#### four
##### five
###### six
"""

m3 = md"""
This is an example of markdown content rendered in Term.jl.
You can use markdown syntax to make words **bold** and *italic* or insert `literals`.

You markdown can include in-line latex ``\LaTeX  \frac{{1}}{{2}}`` and maths in a new line too:

```math
f(a) = \frac{1}{2\pi}\int_{0}^{2\pi} (\alpha+R\cos(\theta))d\theta
```

You can also have links: [Julia](http://www.julialang.org) and
footnotes [^1] for your content [^named].

And, of course, you can show some code too:

```julia
function say_hi(x)
    print("Hellow World")
end
```

---

You can use "quotes" to highlight a section:

> Multi-line quotes can be helpful to make a
> paragram stand out, so that users won't miss it!
> You can use **other inline syntax** in you `quotes` too.

but if you really need to grab someone's attention, use admonitions:

!!! note
    You can use different levels

!!! warning
    to send different messages

!!! danger
    to your reader

!!! tip "Wow!"
    Turns out that admonitions can be pretty useful!
    What will you use them for?

---

Of course you can have classic lists:
* item one
* item two
* And a sublist:
    + sub-item one
    + sub-item two

and ordered lists too:
1. item one
2. item two
3. item three

!!! note "Tables"
    You can use the [Markdown table syntax](https://www.markdownguide.org/extended-syntax/#tables)
    to insert tables - Term.jl will convert them to Table object!

| Term | handles | tables|
|:---------- | ---------- |:------------:|
| Row `1`    | Column `2` |              |
| *Row* 2    | **Row** 2  | Column ``3`` |

----

This is where you print the content of your foot notes:

[^1]: Numbered footnote text.

[^note]:
    Named footnote text containing several toplevel elements.

"""

for (i, m) in enumerate([m1, m2, m3])
    # tprintln(parse_md(m; width=60))
    # tprint(m)
    t = parse_md(m; width = 60)
    tofile(string(t), "./txtfiles/markdown_$i.txt")
end

@testset "Test Markdown Strings" begin
    for (i, m) in enumerate([m1, m2, m3])
        t = parse_md(m; width = 60)
        @test fromfile("./txtfiles/markdown_$i.txt") == t
    end
end


text:
# This works best in the REPL
using ObjCWrap

# Load AppKit which contains the NSSound class
framework("AppKit")

# Load NSSound class. In Julia terms, this means setting a variable
# named NSSound to point to the NSSound class object
@classes NSSound

function play(name::String)
  @objc begin
    sound = [NSSound soundNamed:name]
    if [sound isPlaying] |> Bool
      [sound stop]
    end
    [sound play]
  end
end

play("Purr")


text:
# Autogenerated wrapper script for XML2_jll for armv6l-linux-musleabihf
export libxml2, xmlcatalog, xmllint

using Zlib_jll
using Libiconv_jll
JLLWrappers.@generate_wrapper_header("XML2")
JLLWrappers.@declare_library_product(libxml2, "libxml2.so.2")
JLLWrappers.@declare_executable_product(xmlcatalog)
JLLWrappers.@declare_executable_product(xmllint)
function __init__()
    JLLWrappers.@generate_init_header(Zlib_jll, Libiconv_jll)
    JLLWrappers.@init_library_product(
        libxml2,
        "lib/libxml2.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_executable_product(
        xmlcatalog,
        "bin/xmlcatalog",
    )

    JLLWrappers.@init_executable_product(
        xmllint,
        "bin/xmllint",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
## ------------------------------------------------------------------
macro new_note!(ex...)
    label = _insert_label!("new_note!", ex, __source__)
    _add_secobj!(() -> RBNote(label), __source__)
end


text:
using LinearAlgebra

"""
Construct a statically-sized, mutable array of dimensions S (expressed as a Tuple{...}) using the data from a.
# Example
```
Matrix3(a0,a1,a2,a3,a4,a5,a6,a7,a8)= Matrix3([a0 a1 a2 ; a3 a4 a5 ; a6 a7 a8])
```
"""
const Matrix3=MMatrix{3, 3, Float64}
Matrix3(a0,a1,a2,a3,a4,a5,a6,a7,a8)= Matrix3([a0 a1 a2 ; a3 a4 a5 ; a6 a7 a8])
Matrix3() =  Matrix4(1,0,0, 0,1,0, 0,0,1)

"""
Construct a statically-sized, mutable array of dimensions S (expressed as a Tuple{...}) using the data from a.
# Example
```
Matrix4(a0,a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12,a13,a14,a15)= Matrix4([a0 a1 a2 a3 ; a4 a5 a6 a7; a8 a9 a10 a11; a12 a13 a14 a15])
```
"""
const Matrix4=MMatrix{4, 4, Float64}
Matrix4(a0,a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12,a13,a14,a15)= Matrix4([a0 a1 a2 a3 ; a4 a5 a6 a7; a8 a9 a10 a11; a12 a13 a14 a15])
Matrix4() =  Matrix4(1,0,0,0, 0,1,0,0, 0,0,1,0, 0,0,0,1)

 """
 	flatten(T::Matrix3)

Flatten `T`  by row: value of `Matrix3` into a 1D array,.
 # Example
 ```
 julia> T = convert(GL.Matrix3, rand(Float64, 3,3))
 3×3 StaticArrays.MArray{Tuple{3,3},Float64,2,9}:
  0.774754  0.681844  0.642835
  0.686775  0.448813  0.0144597
  0.824035  0.438451  0.0153561

 julia> GL.flatten(T::GL.Matrix3)
 9-element Array{Float32,1}:
  0.77475363
  0.6818439
	...
  0.4384512
  0.01535608
 ```
 """
function flatten(T::Matrix3)
	return Vector{Float32}([
		T[1,1],T[1,2],T[1,3],
		T[2,1],T[2,2],T[2,3],
		T[3,1],T[3,2],T[3,3]])
end

"""
	flatten(T::Matrix4)

Flatten `T` by row a value of `Matrix4` into a 1D array.

# Example

```
julia> T = convert(GL.Matrix4, rand(Float64, 4,4))
4×4 StaticArrays.MArray{Tuple{4,4},Float64,2,16}:
 0.23079   0.620179  0.362604   0.136963
	...		 ...	   ...		  ...
 0.328367  0.709699  0.170108   0.848778

julia> GL.flatten(T::GL.Matrix4)
16-element Array{Float32,1}:
 0.23079032
 0.6201795
	...
 0.17010835
 0.84877783
```
"""
function flatten(T::Matrix4)
	return Vector{Float32}([
		T[1,1],T[1,2],T[1,3],T[1,4],
		T[2,1],T[2,2],T[2,3],T[2,4],
		T[3,1],T[3,2],T[3,3],T[3,4],
		T[4,1],T[4,2],T[4,3],T[4,4]])
end

"""
	dropW(T::Matrix4)

Remove the last row and column of a `Matrix4`.
# Example
```
julia> typeof(T)
StaticArrays.MArray{Tuple{4,4},Float64,2,16}

julia> GL.dropW(T)
3×3 StaticArrays.MArray{Tuple{3,3},Float64,2,9}:
 0.23079   0.620179  0.362604
 0.650641  0.431527  0.0229223
 0.534276  0.190774  0.946666
```
"""
function dropW(T::Matrix4)
	return Matrix3(
		T[1,1],T[1,2],T[1,3],
		T[2,1],T[2,2],T[2,3],
		T[3,1],T[3,2],T[3,3])
end

"""
	translateMatrix(vt::Point3d)

Build a translation `Matrix4` from a translaction vector `Point3d`.
# Example
```
julia> GL.translateMatrix(GL.Point3d(1,2,3))
4×4 MArray{Tuple{4,4},Float64,2,16}:
 1.0  0.0  0.0  1.0
 0.0  1.0  0.0  2.0
 0.0  0.0  1.0  3.0
 0.0  0.0  0.0  1.0
```
"""
function translateMatrix(vt::Point3d)
	return Matrix4(
		1.0, 0.0, 0.0, vt[1],
		0.0, 1.0, 0.0, vt[2],
		0.0, 0.0, 1.0, vt[3],
		0.0, 0.0, 0.0, 1.0)
end

"""
	scaleMatrix(vs::Point3d)

Build a scaling `Matrix4` from a `vs` scaling vector `Point3d`.
# Example
```
julia> GL.scaleMatrix(GL.Point3d(1.1,1.2,1.3))
4×4 MArray{Tuple{4,4},Float64,2,16}:
 1.1  0.0  0.0  0.0
 0.0  1.2  0.0  0.0
 0.0  0.0  1.3  0.0
 0.0  0.0  0.0  1.0
```
"""
function scaleMatrix(vs::Point3d)
	return Matrix4(
		vs[1], 0.0, 0.0, 0.0,
		0.0, vs[2], 0.0, 0.0,
		0.0, 0.0, vs[3], 0.0,
		0.0, 0.0, 0.0, 1.0)
end

"""
	lookAtMatrix(eye::Point3d, center::Point3d, up::Point3d)

Build a *lookAtMatrix* `Matrix4`.

In order to set a camera position, all really needed is a point to set the *camera position* in world space, referred to as the `from` point, and a point that defines *where the camera is looking at*, referred as the `to` point. The camera orintation is specified by the `up` vector, whose projection is seen vertical in screen coordinates.

# Example
```
julia> GL.lookAtMatrix(GL.Point3d(2,2,2), GL.Point3d(0,0,0), GL.Point3d(0,0,1))
4×4 MArray{Tuple{4,4},Float64,2,16}:
 -0.707107   0.707107  0.0        0.0
 -0.408248  -0.408248  0.816497   0.0
  0.57735    0.57735   0.57735   -3.4641
  0.0        0.0       0.0        1.0
```
"""
function lookAtMatrix(eye::Point3d, center::Point3d, up::Point3d)
	forward = normalized(center-eye)
	side    = normalized(cross(forward,up))
	up      = cross(side,forward)
	m = Matrix4(
		side[1],up[1],-forward[1], 0.0,
		side[2],up[2],-forward[2], 0.0,
		side[3],up[3],-forward[3], 0.0,
		    0.0,  0.0,        0.0, 1.0
	)
	return transpose(m) * translateMatrix(-1*eye)
end

"""
	perspectiveMatrix(fovy::Float64, aspect::Float64, zNear::Float64, zFar::Float64)

Generation of the perspective matrix. Using 4 real parameters:
-	`fovy::Float64` is the ``y`` component of ``field-ov-view`` in radians
-	`aspect::Float64` is the ratio ``width/height`` of viewport in screen space
-	`zNear::Float64` is -1 in clip space.
-	`zFar::Float64` is a positive number in clip space.
"""
function perspectiveMatrix(fovy::Float64, aspect::Float64, zNear::Float64, zFar::Float64)
	radians =  deg2rad(fovy/2.0)
	cotangent = cos(radians) / sin(radians)
	m=Matrix4()
	m[1,1] = cotangent / aspect
	m[2,2] = cotangent
	m[3,3] = -(zFar + zNear) / (zFar - zNear)
	m[3,4] = -1.0
	m[4,3] = -2.0 * zNear * zFar / (zFar - zNear)
	m[4,4] =  0.0
	return transpose(m)
end

function orthoMatrix(left::Float64, right::Float64, bottom::Float64, top::Float64, nearZ::Float64, farZ::Float64)
	m=Matrix4()
	m[1,1] = 2 / (right-left); m[1,2] =                0; m[1,3] =                 0; m[1,4] = -(right+left) / (right-left)
	m[2,1] =                0; m[2,2] = 2 / (top-bottom); m[2,3] =                 0; m[2,4] = -(top+bottom) / (top-bottom)
	m[3,1] =                0; m[3,2] =                0; m[3,3] = -2 / (farZ-nearZ); m[3,4] = -(farZ+nearZ) / (farZ-nearZ)
	m[4,1] =                0; m[4,2] =                0; m[4,3] =                 0; m[4,4] = 1
	return m
end

"""
	getLookAt(T::Matrix4)

Return `pos`,`dir`,`vup` vectors in world space.
"""
function getLookAt(T::Matrix4)
	T=inv(T)
	pos=          (Point3d(  T[1,4], T[2,4], T[3,4]))
	dir=normalized(Point3d( -T[1,3],-T[2,3],-T[3,3]))
	vup=normalized(Point3d(  T[1,2], T[2,2], T[3,2]))
	pos,dir,vup
end


text:
# differentiated version of homotopy.jl
import PDESolver.evalHomotopyJacobian

function evalHomotopyJacobian(mesh::AbstractMesh, sbp::AbstractOperator,
                              eqn::EulerData{Tsol, Tres}, opts::Dict,
                              assembler::AssembleElementData, lambda::Number
                             ) where {Tsol, Tres}

  if opts["homotopy_function"] == "ViscousBO"
    sensor = ShockSensorBO{Tsol, Tres}(mesh, sbp, opts)
    sensor.alpha = lambda
    capture = eqn.shock_capturing
    applyShockCapturing_diff(mesh, sbp, eqn, opts, sensor, capture, assembler)

  elseif opts["homotopy_function"] == "FirstOrderDissipation"
    calcHomotopyDiss_jac(mesh, sbp, eqn, opts, assembler, lambda)
  end
end

function calcHomotopyDiss_jac(mesh::AbstractDGMesh{Tmsh}, sbp,
                          eqn::EulerData{Tsol, Tres}, opts, assembler, lambda) where {Tsol, Tres, Tmsh}

  # some checks for when parallelism is enabled
  @assert opts["parallel_data"] == PARALLEL_DATA_ELEMENT
  for i=1:mesh.npeers
    @assert eqn.shared_data[i].recv_waited
  end

  params = eqn.params
  #----------------------------------------------------------------------------
  # volume dissipation

  # compute the D operator in each direction
  D = zeros(mesh.numNodesPerElement, mesh.numNodesPerElement, mesh.dim)
  for d=1:mesh.dim
    D[:, :, d] = inv(diagm(sbp.w))*sbp.Q[:, :, d]
  end

  res_jac = eqn.params.calc_face_integrals_data.res_jacLL
  t2_dot = eqn.params.calc_face_integrals_data.res_jacLR  # work array
  t1 = zeros(Tsol, mesh.numDofPerNode, mesh.numNodesPerElement)
  lambda_dot = zeros(Tres, mesh.numDofPerNode)

  nrm = zeros(Tmsh, mesh.dim)
  for el=1:mesh.numEl
    q_el = sview(eqn.q, :, :, el)
    fill!(res_jac, 0.0)
    for d1=1:mesh.dim
      fill!(t1, 0.0)
      fill!(t2_dot, 0.0)

      differentiateElement!(sbp, d1, q_el, t1)

      # t2 = t1*lambda_p
      # contribution t2_dot = lambda_p*t1_dot
      for p=1:mesh.numNodesPerElement
        q_p = sview(eqn.q, :, p, el)

        # get vector in xi direction defined by dim
        for k=1:mesh.dim
          nrm[k] = mesh.dxidx[d1, k, p, el]
        end

        lambda_max = getLambdaMax(eqn.params, q_p, nrm)

        for q=1:mesh.numNodesPerElement
          for j=1:mesh.numDofPerNode
#            for i=1:mesh.numDofPerNode
              t2_dot[j, j, p, q] += lambda_max*D[p, q, d1]
#            end
          end
        end
      end   # end loop p

      # contribution t2_dot += t1*lambda_p_dot
      for p=1:mesh.numNodesPerElement
        q_p = sview(eqn.q, :, p, el)

        # get vector in xi direction defined by dim
        for k=1:mesh.dim
          nrm[k] = mesh.dxidx[d1, k, p, el]
        end

        getLambdaMax_diff(eqn.params, q_p, nrm, lambda_dot)

        for j=1:mesh.numDofPerNode
          for i=1:mesh.numDofPerNode
            t2_dot[i, j, p, p] += t1[i, p]*lambda_dot[j]
          end
        end
      end  # end loop p

      # apply Q^T
#      for d2=1:mesh.dim
        for p=1:mesh.numNodesPerElement
          for q=1:mesh.numNodesPerElement
            for c=1:mesh.numNodesPerElement
              for j=1:mesh.numDofPerNode
                for i=1:mesh.numDofPerNode
                  res_jac[i, j, p, q] += sbp.Q[c, p, d1]*t2_dot[i, j, c, q]
                end
              end
            end
          end
        end
#      end  # end loop d2

    end  # end loop d1

    # negate res_jac for consistency with physics module
    for i=1:length(res_jac)
      res_jac[i] = -lambda*res_jac[i]
    end

    assembleElement(assembler, mesh, el, res_jac)
  end  # end loop el

  fill!(t2_dot, 0.0)

  #----------------------------------------------------------------------------
  # interface terms

  @unpack params.calc_face_integrals_data q_faceL q_faceR flux_dotL flux_dotR res_jacLL res_jacLR res_jacRL res_jacRR

  #=
  q_faceL = zeros(Tsol, mesh.numDofPerNode, mesh.numNodesPerFace)
  q_faceR = zeros(q_faceL)
#  nrm2 = eqn.params.nrm2
  flux_dotL = zeros(Tres, mesh.numDofPerNode, mesh.numDofPerNode, mesh.numNodesPerFace)
  flux_dotR = zeros(flux_dotL)

  res_jacLL = params.res_jacLL
  res_jacLR = params.res_jacLR
  res_jacRL = params.res_jacRL
  res_jacRR = params.res_jacRR
  =#
  lambda_dotL = zeros(Tres, mesh.numDofPerNode)
  lambda_dotR = zeros(Tres, mesh.numDofPerNode)

  for i=1:mesh.numInterfaces
    iface_i = mesh.interfaces[i]
    qL = sview(eqn.q, :, :, iface_i.elementL)
    qR = sview(eqn.q, :, :, iface_i.elementR)
    fill!(res_jacLL, 0.0)
    fill!(res_jacLR, 0.0)
    fill!(res_jacRL, 0.0)
    fill!(res_jacRR, 0.0)

    interiorFaceInterpolate!(mesh.sbpface, iface_i, qL, qR, q_faceL, q_faceR)

    # calculate the flux jacobian at each face node
    for j=1:mesh.numNodesPerFace
      qL_j = sview(q_faceL, :, j)
      qR_j = sview(q_faceR, :, j)

      # get the face normal
      nrm2 = sview(mesh.nrm_face, :, j, i)

      lambda_max = getLambdaMaxSimple_diff(eqn.params, qL_j, qR_j, nrm2,
                                           lambda_dotL, lambda_dotR)

      for k=1:mesh.numDofPerNode
        # flux[k, j] = 0.5*lambda_max*(qL_j[k] - qR_j[k])
        for m=1:mesh.numDofPerNode
          flux_dotL[m, k, j] = 0.5*lambda_dotL[k]*(qL_j[m] - qR_j[m])
          flux_dotR[m, k, j] = 0.5*lambda_dotR[k]*(qL_j[m] - qR_j[m])
        end
        flux_dotL[k, k, j] += 0.5*lambda_max
        flux_dotR[k, k, j] -= 0.5*lambda_max
      end
    end  # end loop j

    # multiply by lambda here and it will get carried through
    # interiorFaceCombined_jac
    scale!(flux_dotL, lambda)
    scale!(flux_dotR, lambda)

    # compute dR/dq
    interiorFaceCombined_jac!(mesh.sbpface, iface_i, flux_dotL, flux_dotR,
                             res_jacLL, res_jacLR, res_jacRL, res_jacRR,
                             SummationByParts.Subtract())
    # assemble into the Jacobian
    assembleInterface(assembler, mesh.sbpface, mesh, iface_i, res_jacLL, res_jacLR,
                                                res_jacRL, res_jacRR)

  end  # end loop i

  #----------------------------------------------------------------------------
  # skipping boundary integrals
  # use nrm2, flux_jfacL from interface terms above
  if opts["homotopy_addBoundaryIntegrals"]
    qg = zeros(Complex128, mesh.numDofPerNode)
    q_faceLc = zeros(Complex128, mesh.numDofPerNode, mesh.numNodesPerFace)
    h = 1e-20
    pert = Complex128(0, h)
    for i=1:mesh.numBoundaryFaces
      bndry_i = mesh.bndryfaces[i]
      qL = sview(eqn.q, :, :, bndry_i.element)
#      resL = sview(res, :, :, bndry_i.element)
      fill!(q_faceLc, 0.0)
      fill!(res_jac, 0.0)

      boundaryFaceInterpolate!(mesh.sbpface, bndry_i.face, qL, q_faceLc)

      # compute flux jacobian at each node
      for j=1:mesh.numNodesPerFace
        q_j = sview(q_faceLc, :, j)
        for m=1:mesh.numDofPerNode
          q_j[m] += pert
    #      dxidx_j = sview(mesh.dxidx_bndry, :, :, j, i)

          # calculate boundary state
          coords = sview(mesh.coords_bndry, :, j, i)
          calcFreeStream(eqn.params_complex, coords, qg)

          # calculate face normal
          nrm2 = sview(mesh.nrm_bndry, :, j, i)

          # calculate lambda_max
          lambda_max = getLambdaMaxSimple(eqn.params_complex, q_j, qg, nrm2)

          # calculate dissipation
          for k=1:mesh.numDofPerNode
            flux_dotL[k, m, j] = lambda*imag(0.5*lambda_max*(q_j[k] - qg[k]))/h
          end

          q_j[m] -= pert
        end  # end loop m
      end  # end loop j

      boundaryFaceCombined_jac!(mesh.sbpface, bndry_i.face, flux_dotL, res_jac,
                               SummationByParts.Subtract())

      assembleBoundary(assembler, mesh.sbpface, mesh, bndry_i, res_jac)
    end  # end loop i
  end

  #----------------------------------------------------------------------------
  # shared face integrals
  # use q_faceL, q_faceR, lambda_dotL, lambda_dotR, flux_dotL, flux_dotR
  # from above

  workarr = zeros(q_faceR)
  for peer=1:mesh.npeers
    # get data for this peer
    interfaces_peer = mesh.shared_interfaces[peer]

    qR_peer = eqn.shared_data[peer].q_recv
#    dxidx_peer = mesh.dxidx_sharedface[peer]
    nrm_peer = mesh.nrm_sharedface[peer]
    start_elnum = mesh.shared_element_offsets[peer]

    for i=1:length(interfaces_peer)
      iface_i = interfaces_peer[i]
      qL_i = sview(eqn.q, :, :, iface_i.elementL)
      qR_i = sview(qR_peer, :, :, iface_i.elementR - start_elnum + 1)
      fill!(res_jacLL, 0.0)
      fill!(res_jacLR, 0.0)

      # interpolate to face
      interiorFaceInterpolate!(mesh.sbpface, iface_i, qL_i, qR_i, q_faceL, q_faceR)
      # compute flux at every face node
      for j=1:mesh.numNodesPerFace
        qL_j = sview(q_faceL, :, j)
        qR_j = sview(q_faceR, :, j)
        nrm2 = sview(nrm_peer, :, j, i)

        # get max wave speed
        lambda_max = getLambdaMaxSimple_diff(eqn.params, qL_j, qR_j, nrm2,
                                             lambda_dotL, lambda_dotR)

        # calculate flux
        for k=1:mesh.numDofPerNode
          # flux[k, j] = 0.5*lambda_max*(qL_j[k] - qR_j[k])
          for m=1:mesh.numDofPerNode
            flux_dotL[m, k, j] = 0.5*lambda_dotL[k]*(qL_j[m] - qR_j[m])
            flux_dotR[m, k, j] = 0.5*lambda_dotR[k]*(qL_j[m] - qR_j[m])
          end
          flux_dotL[k, k, j] += 0.5*lambda_max
          flux_dotR[k, k, j] -= 0.5*lambda_max
        end
      end  # end loop j

      # multiply by lambda here and it will get carried through
      # interiorFaceCombined_jac
      scale!(flux_dotL, lambda)
      scale!(flux_dotR, lambda)

      # compute dR/dq
      interiorFaceCombined_jac!(mesh.sbpface, iface_i, flux_dotL, flux_dotR,
                                res_jacLL, res_jacLR, res_jacRL, res_jacRR,
                                SummationByParts.Subtract())

     assembleSharedFace(assembler, mesh.sbpface, mesh, iface_i, res_jacLL, res_jacLR)
    end  # end loop i
  end  # end loop peer

  return nothing
end

"""
  Differentiated version of [`getLambdaMaxSimple`](@ref)

  **Inputs**

   * params
   * qL
   * qR
   * dir

  **Inputs/Outputs**

   * lambda_dotL: derivative of lambda wrt. qL
   * lambda_dotR: derivative of lambda wrt. qR
"""
function getLambdaMaxSimple_diff(params::ParamType{Tdim},
                      qL::AbstractVector{Tsol}, qR::AbstractVector{Tsol},
                      dir::AbstractVector{Tmsh},
                      lambda_dotL::AbstractVector{Tres},
                      lambda_dotR::AbstractVector{Tres}) where {Tsol, Tres, Tmsh, Tdim}

  q_avg = params.get_lambda_max_simple_data.q_avg

  for i=1:length(q_avg)
    q_avg[i] = 0.5*(qL[i] + qR[i])
  end

  lambda_max = getLambdaMax_diff(params, q_avg, dir, lambda_dotL)

  for i=1:length(lambda_dotL)
    lambda_dotL[i] *= 0.5
    lambda_dotR[i] = lambda_dotL[i]
  end

  return lambda_max
end


text:
function stepSparseCoding(ihtsc::Type{SparseCoders.IHT}, options, cb::DT, x::AbstractArray; shape::Shapes.AbstractShape=Shapes.Vec(size(x)), vlevel::Integer=0, kwargs...) where {DT<:LearningTarget}
    ts = createTransform(cb, shape)
    # number of non-zero coefficients
    iht = ihtsc(x, t->synthesize(ts, t), t->analyze(ts, t); options...)

    # initial sparse vector y0
    y0 = analyze(ts, x)
    y_opt, loss_iht = iht(y0, isverbose=(vlevel>=3))

    return (y_opt, loss_iht)
end

function stepSparseCoding(istasc::Type{TI}, options, cb::DT, x::AbstractArray; shape::Shapes.AbstractShape=Shapes.Vec(size(x)), vlevel::Integer=0, kwargs...) where {TI<:Union{SparseCoders.ISTA,SparseCoders.FISTA},DT<:LearningTarget}
    ts = createTransform(cb, shape)
    # number of non-zero coefficients
    ∇f = (_y) -> begin
        -real.(analyze(ts, x - synthesize(ts, _y)))
    end
    ista = istasc(∇f; options...)

    # initial sparse vector y0
    y0 = analyze(ts, x)
    y_opt, _dummy = ista(y0)

    loss_ista_x = norm(x - synthesize(ts, y_opt))^2/2
    return (y_opt, loss_ista_x)
end

max_stars_repo_name: nixir/MDCDL

max_issues_repo_name: nixir/MDCDL

max_forks_repo_name: nixir/MDCDL


text:
funname(x) = split(string(x), ".")[end]
asarray(x) = [x]
asarray(x::Tuple) = collect(x)
asarray(x::AbstractArray) = x

getselectionpairs(df, s...) = reduce_renames(postprocess(df, resolve_query(df, [s...])))
getselectionkeys(df, s...) = first.(getselectionpairs(df, s...))

function whichkeys(d, v)
    out = Vector{Symbol}()
    for (k,vals) in d
        if v in vals
            push!(out, k)
        end
    end
    length(out) == 0 ? nothing : out
end


text:
module LaTeXTabulars

using ArgCheck
using DocStringExtensions
using Parameters

export Rule, CMidRule, MultiColumn, Tabular, latex_tabular

# cells

"""
    $SIGNATURES

Print a the contents of `cell` to `io` as LaTeX.

!!! NOTE

    Methods are defined for some specific types, but if you want full control
    (eg rounding), use an `<: AbstractString`, eg `String` or `LaTeXString`.
"""
function latex_cell(io::IO, cell::T) where T
    info("Define a `latex_cell` for writing $T objects to LaTeX.")
    throw(MethodError(latex_cell, Tuple{IO, T}))
end

latex_cell(io::IO, x::Real) = print(io, string(x))

latex_cell(io::IO, s::AbstractString) = print(io, s)

"""
    MultiColumn(n, pos, cell)

For `\\multicolumn{n}{pos}{cell}`. Use the symbols `:l`, `:c`, `:r` for `pos`.
"""
struct MultiColumn
    n::Int
    pos::Symbol
    cell
end

function latex_cell(io::IO, mc::MultiColumn)
    @unpack n, pos, cell = mc
    @argcheck(pos ∈ (:l, :c, :r),
              "$pos is not a recognized position. Use :l, :c, :r.")
    print(io, "\\multicolumn{$(n)}{$(pos)}{")
    latex_cell(io, mc.cell)
    print(io, "}")
end

# non-cell-like objects

struct Rule{T} end

"""
    $SIGNATURES

Horizontal rule. The `kind` of the rule is specified by a symbol, which will
generally be printed as `\\KINDrule` for rules in `booktabs`, eg `Rule(:top)`
prints `\\toprule`. To obtain a `\\hline`, use `Rule{:h}`.
"""
Rule(kind::Symbol = :h) = Rule{kind}()

latex_line(io::IO, ::Rule{:top}) = println(io, "\\toprule ")

latex_line(io::IO, ::Rule{:mid}) = println(io, "\\midrule ")

latex_line(io::IO, ::Rule{:bottom}) = println(io, "\\bottomrule ")

latex_line(io::IO, ::Rule{:h}) = println(io, "\\hline ")

latex_line(io::IO, r::Rule) = error("Don't know how to print $(typeof(r)).")

"""
    CMidRule([wd], [trim], left, right)

Will be printed as `\\cmidrule[wd](trim)[left-right]`. When `wd` or `trim` is
`nothing`, it is omitted. Use with the `booktabs` LaTeX package.
"""
struct CMidRule
    wd::Union{Void, AbstractString}
    trim::Union{Void, AbstractString}
    left::Int
    right::Int
    function CMidRule(wd, trim, left, right)
        @argcheck 1 ≤ left ≤ right
        new(wd, trim, left, right)
    end
end

CMidRule(trim, left, right) = CMidRule(nothing, trim, left, right)

CMidRule(left, right) = CMidRule(nothing, left, right)

function latex_line(io::IO, rule::CMidRule)
    @unpack wd, trim, left, right = rule
    print(io, "\\cmidrule")
    wd ≠ nothing && print(io, "[$(wd)]")
    trim ≠ nothing && print(io, "($(trim))")
    print(io, "{$(left)-$(right)} ") # NOTE trailing space important
end

function latex_line(io::IO, cells)
    for (column, cell) in enumerate(cells)
        column == 1 || print(io, " & ")
        latex_cell(io, cell)
    end
    println(io, " \\\\")
end

function latex_line(io::IO, M::AbstractMatrix)
    for i in indices(M, 1)
        latex_line(io, M[i, :])
    end
end

function latex_line(io::IO, lines::Tuple)
    for line in lines
        latex_line(io, line)
    end
end

# tabular and similar environments

abstract type TabularLike end

"""
    Tabular(cols)

For the LaTeX environment `\\begin{tabular}{cols} ... \\end{tabular}`.
"""
struct Tabular <: TabularLike
    "A column specification, eg `\"llrr\"`."
    cols::AbstractString
end

latex_env_begin(io::IO, t::Tabular) = println(io, "\\begin{tabular}{$(t.cols)}")

latex_env_end(io::IO, t::Tabular) = println(io, "\\end{tabular}")

"""
    $SIGNATURES

Print `lines` to `io` as a LaTeX using the given environment.

Each `line` in `lines` can be

- a rule-like object, eg [`Rule`] or [`CMidRule`],

- an iterable (eg `AbstractVector`) of cells,

- a `Tuple`, which is treated as multiple lines (“splat” in place), which is
  useful for functions that generate lines with associated rules, or multiple
  `CMidRule`s,

- a matrix, each row of which is treated as a line.

See [`latex_cell`](@ref) for the kinds of cell supported (particularly
[`MultiColumn`](@ref), but for full formatting control, use an `String` or
`LaTeXString` for cells.
"""
function latex_tabular(io::IO, t::TabularLike, lines)
    latex_env_begin(io, t)
    for line in lines
        latex_line(io, line)
    end
    latex_env_end(io, t)
end

latex_tabular(io::IO, t::TabularLike, lines::AbstractMatrix) =
    latex_tabular(io, t, [lines])

"""
    $SIGNATURES

LaTeX output as a string.
"""
function latex_tabular(::Type{String}, t::TabularLike, lines)
    io = IOBuffer()
    latex_tabular(io, t, lines)
    String(take!(io))
end

"""
    $SIGNATURES

Write a `tabular`-like LaTeX environment to `filename`, which is **overwritten**
if it already exists.
"""
function latex_tabular(filename::AbstractString, t::TabularLike, lines)
    open(filename, "w") do io
        latex_tabular(io, t, lines)
    end
end

end # module


text:
random_seed = nothing

data_path = ""
data_prefix = "2d1ksample"

#Model Parameters
iterations = 99
hard_clustering = false
argmax_sample_stop = 10 #Change to hard assignment from soft at iterations - argmax_sample_stop
split_stop  = 10 #Stop split/merge moves at  iterations - split_stop

total_dim = 2

α = 100000.0

initial_clusters = 1

use_dict_for_global = false

hyper_params = niw_hyperparams(1.0,
    [0,0],
    5.0,
    Matrix{Float64}(I, 2, 2)*1.0)

#Saving specifics:
enable_saving = true
model_save_interval = 50
save_path = ""
overwrite_prec = false
save_file_prefix = "checkpoint_"


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

using Test

using Base.LinAlg: BlasComplex, BlasFloat, BlasReal, QRPivoted

let n = 10
    srand(1234321)

    Areal  = randn(n,n)/2
    Aimg   = randn(n,n)/2

    @testset for eltya in (Float32, Float64, Complex64, Complex128, BigFloat, Int)
        A = eltya == Int ?
                rand(1:7, n, n) :
                convert(Matrix{eltya}, eltya <: Complex ?
                    complex.(Areal, Aimg) :
                    Areal)

        if eltya != BigFloat
            H = hessfact(A)
            @test size(H[:Q], 1) == size(A, 1)
            @test size(H[:Q], 2) == size(A, 2)
            @test size(H[:Q]) == size(A)
            @test_throws KeyError H[:Z]
            @test convert(Array, H) ≈ A
            @test (H[:Q] * H[:H]) * H[:Q]' ≈ A
            @test (H[:Q]' *A) * H[:Q] ≈ H[:H]
            #getindex for HessenbergQ
            @test H[:Q][1,1] ≈ Array(H[:Q])[1,1]
        end
    end
end

max_stars_repo_name: xorJane/julia

max_issues_repo_name: JuliaLangEs/julieta

max_forks_repo_name: JuliaLangEs/julieta


text:
using CompScienceMeshes
using Test

#include(Pkg.dir("CompScienceMeshes","examples","matlab_patches.jl"))

sphere = readmesh(joinpath(dirname(@__FILE__),"assets","sphere2.in"))
index_in_sphere = mapper(sphere)
for (i,_c) in enumerate(cells(sphere))
    @eval @test index_in_sphere[$_c] == $i
end

isdownbelow(c) = center(chart(sphere,c))[3] < 0
southern_hemisphere = submesh(isdownbelow, sphere)

# build the restriction operator
R = spzeros(Int, numcells(southern_hemisphere), numcells(sphere))
for (i,_c) in enumerate(cells(southern_hemisphere))
    j = index_in_sphere[_c]
    R[i,j] = 1
end

@test nnz(R) == numcells(southern_hemisphere)

rect = meshrectangle(1.0, 1.0, 1.0, 3)
@test numcells(rect) == 2
tri = Mesh(rect.vertices, rect.faces[1:1])
@test numcells(tri) == 1

R = restriction(tri, rect)
@test R == [1 0]


text:
using CUSPARSE
using CUDArt
using Base.Test

m = 25
n = 35
k = 10
blockdim = 5

#############
# test_doti #
#############

function test_doti(elty)
    x = sparsevec(rand(1:m,k), rand(elty,k), m)
    y = rand(elty,m)
    d_x = CudaSparseVector(x)
    d_y = CudaArray(y)
    ddot = CUSPARSE.doti(d_x,d_y,'O')
    #compare
    dot = zero(elty)
    for i in 1:length(x.nzval)
        if VERSION >= v"0.5.0-dev+742"
            dot += x.nzval[i] * y[x.nzind[i]]
        else
            dot += x.nzval[i] * y[x.rowval[i]]
        end
    end
    @test ddot ≈ dot
end

##############
# test_dotci #
##############

function test_dotci(elty)
    x = sparsevec(rand(1:m,k), rand(elty,k), m)
    y = rand(elty,m)
    d_x = CudaSparseVector(x)
    d_y = CudaArray(y)
    ddot = CUSPARSE.dotci(d_x,d_y,'O')
    #compare
    dot = zero(elty)
    for i in 1:length(x.nzval)
        if VERSION >= v"0.5.0-dev+742"
            dot += conj(x.nzval[i]) * y[x.nzind[i]]
        else
            dot += conj(x.nzval[i]) * y[x.rowval[i]]
        end
    end
    @test ddot ≈ dot
end

types = [Float32,Float64,Complex64,Complex128]
for elty in types
    tic()
    test_doti(elty)
    println("doti took ", toq(), " for ", elty)
end
types = [Complex64,Complex128]
for elty in types
    tic()
    test_dotci(elty)
    println("dotci took ", toq(), " for ", elty)
end


text:
using Plots, GraphRecipes
gr(size=(400,400))

"""
## mesh

Plots displacements and directions

### Function
```julia
mesh(data, g_coord, g_num, disp, ampl, pdir)
```

### Arguments
```julia
* data::Dict                 : Input dictionary
* g_coord::Array{Float64, 2} : Coordinate array
* g_num::Array{Int, 2}       : Global node numbering array
* disp::DataFrame            : Displacements DataFrame
* ampl::Float64              : Amplification for derivatives
* pdir::AbstractString       : Nodal freedom matrix (updated)
```
"""
function mesh(data::Dict, g_coord::Array{Float64,2}, g_num::Array{Int, 2},
    disp, ampl, pdir)
  Plots.plot(leg=false, title="Nodes (blue), $(ampl)*displacements (red)")
  for xc in data[:x_coords]
    for yc in data[:y_coords]
      plot!([xc, xc], [minimum(data[:y_coords]), maximum(data[:y_coords])],
        line=(:solid), color=:lightgray)
      plot!([minimum(data[:x_coords]), maximum(data[:x_coords])],
        [yc, yc], line=(:solid), color=:lightgray)
    end
  end
  vx = ampl*convert(Array, disp[:,1])
  vy = ampl*convert(Array, disp[:,2])
  scatter!(g_coord[1,:], g_coord[2,:], marker=(:circle,2,0.4,stroke(1,:blue)))
  quiver!(g_coord[1,:], g_coord[2,:], quiver=(vx, vy))
  savefig(pdir*"/Ex61.1.jl_dis.pdf")
  plot!()
end


text:
Base.length(cd::CoolDict) = length(getfield(cd, :_data))
Base.iterate(cd::CoolDict, args...; kwargs...) = iterate(getfield(cd, :_data), args...; kwargs...)

Base.getindex(cd::CoolDict, key) = getfield(cd, :_data)[key]
Base.setindex!(cd::CoolDict, value, key) = (getfield(cd, :_data)[key] = value)

Base.propertynames(cd::CoolDict) = Tuple(keys(cd))

function Base.getproperty(cd::CoolDict, key::Symbol)
    if key == :_data
        error("Hey, the $key field is sorta private. If you want to access it, you can use getfield(mycooldict, :$key) though. Sorry.")
    elseif key in keys(cd)
        return getfield(cd, :_data)[key]
    else
        return error("type CoolDict has no field $key")
    end
end
Base.setproperty!(cd::CoolDict, key::Symbol, value) = (getfield(cd, :_data)[key] = value)

Base.keys(cd::CoolDict) = keys(getfield(cd, :_data))
Base.values(cd::CoolDict) = values(getfield(cd, :_data))

Base.convert(T::Type{<:Dict}, cd::CoolDict) = T(getfield(cd, :_data))
Base.convert(::Type{NamedTuple}, cd::CoolDict) = (; getfield(cd, :_data)...)

function Base.show(io::IO, ::MIME"text/plain", cd::CoolDict{V}) where {V}
    println(io, "CoolDict{$V} with $(length(cd)) entries:")
    for (key, val) in cd
        print(io, "  $(string(key)) = ")
        show(io, val)
        println(io,"")
    end
    return nothing
end

Base.show(io::IO, cd::CoolDict) = print(io, "CoolDict$(NamedTuple(cd))")


text:
"""
    Fermi.CoupledCluster.RCCD

Fermi struct that holds information about RCCD wavefunctions

_struct tree:_

**RCCD** <: AbstractCCWavefunction <: AbstractCorrelatedWavefunction <: AbstractWavefunction
"""
mutable struct RCCD{T} <: AbstractCCWavefunction
    CorrelationEnergy::T
    MP2Energy::T
    T2::_T where _T <: Fermi.AbstractTensor
end

# Most general function. It defines the precision and call a precision-specific function
function RCCD()
    precision_selector = Dict{Any,Any}("single" => Float32,
                                       "double" => Float64)
    prec = precision_selector[Fermi.CurrentOptions["precision"]]
    RCCD{prec}()
end

function RCCD{T}() where T <: AbstractFloat
    alg_selector = Dict{Any,Any}("dpd" => DPD(),
                                 "ctf" => CTF())
    int_selector = Dict{Any,Any}("df" => DF(),
                                 "conv" => Conventional(),
                                 "cd" => CD())
    alg = alg_selector[Fermi.CurrentOptions["cc_alg"]]
    int = int_selector[Fermi.CurrentOptions["ints_type"]]
    RCCD{T}(alg,int)
end

#implementations
include("RCCD_DPD_Conventional.jl")

"""
Docstring
"""
function RCCD{T}(ref::Fermi.HartreeFock.RHF.RHFWavefunction,alg::T2,int::T3) where { T <: AbstractFloat,
                                                                                T2 <: CCAlgorithm,
                                                                                T3 <: CCIntegrals }
end


text:
################################################################################
#
#          EllCrv/EllCrvDiv.jl : Divisors on Elliptic Curves and Logarithms of
#                                   Divisors
#
################################################################################

################################################################################
#
#  Submodules
#
################################################################################

include("EllCrv.jl")

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

export AbstractDivisor, EllCrvDivisor, EllCrvModel, Logarithm
export EllipticCurveDivisor, ord, assoc

################################################################################
#
#  Abstract Types/Structs
#
################################################################################

abstract type AbstractDivisor{EllCrvPt} end

mutable struct EllCrvDivisor{EllCrvPt} <: AbstractDivisor{EllCrvPt}
    degree::Int
    coeff::Array{Int, 1}
    points::Tuple{Vararg{EllCrvPt, N} where N}
    func_field_without_zero::AbstractAlgebra.Field
    rat_func::AbstractAlgebra.FieldElem
    principal::Bool

    # blowup(i::AbstractVarietyMap) --> AbstractVariety
    # if exceptionalDiv(blowup) --> EllCrvDivisor
    exceptional::Bool

    # if prod(coeff::Array{T, 1}) >= 0
    effective::Bool

    function EllCrvDivisor{EllCrvPt}(coeffs::Array{Int, 1},
        points::Tuple{Vararg{EllCrvPt, N} where N}, check::Bool = true) where
        EllCrvPt
        if check
            if assoc(rat_func, coeffs, points)
                ECD = new{EllCrvPt}()
                ECD = +(coeffs[1]*EllCrvPt...)
                ECD.principal = true
            elseif prod(coeffs) >= 0
                ECD = new{EllCrvPt}()
                ECD = +(coeffs[1]*EllCrvPt...)
                ECD.effective = true
            end
        else
            ECD = new{EllCrvPt}()
            ECD = +(coeffs[1]*EllCrvPt...)
        end
        return ECD
    end
end

mutable struct EllCrvModel
    curve::EllCrv
    point::EllCrvPt
    divisor::EllCrvDivisor
end

# Logarithms of Divisors
mutable struct Logarithm
    base_ec_model::EllCrvModel
    divisor::EllCrvDivisor
    #...
end

################################################################################
#
#  Constructors for users
#
################################################################################

function EllipticCurveDivisor(coeff::Array{Int, 1},
    points::Tuple{Vararg{EllCrvPt, N} where N}, check::Bool = true)
    if check
        ECD = EllCrvDivisor{EllCrvPt}(coeff, points, check)
        return ECD
    end
end

################################################################################
#
#  Divisorial Type Methods
#
################################################################################

function ord(rat_func::AbstractAlgebra.FieldElem,
    points::Tuple{Vararg{EllCrvPt, N} where N})

    if AbstractAlgebra.denominator(rat_func) == 0

    end

    return
end

function assoc(rat_func::AbstractAlgebra.FieldElem,
    coeff::Array{Int, 1}, points::Tuple{Vararg{EllCrvPt, N} where N})

    coeff = ord(rat_func, points)
    return EllipticCurveDivisor(coeff, points, true)
end


text:
using StructuresKit

#Test bimoment calculation for a single span beam.
#Example II-11, AISI D100-17, 9CS2.5x059, 25 ft. span, load is 10 lbs/ft
#AISC assumes vertical location of load is at center of web, so do that here also

#Ix Iy Ixy J Cw
sectionProperties = [(10.2911, 0.69678, 0.0, 0.001022, 11.9)]

#E  ν
materialProperties = [(29500, 0.30)]

#ax ay
loadLocation = [(-1.05, 0.0)]

#kx kϕ
springStiffness = [(0.0, 0.0)]

#member information
#L dL SectionProperties MaterialProperties LoadLocation SpringStiffness
memberDefinitions = [(25*12,3.0,1,1,1,1)]

#end boundary conditions
#type=1 u''=v''=ϕ''=0 (simply supported), type=2 u'=v'=ϕ'=0  (fixed), type=3 u''=v''=ϕ''=u'''=v'''=ϕ'''=0 (free end, e.g., a cantilever)
endBoundaryConditions = [1 1]

#supports
#location where u=v=ϕ=0
supports = [0.0 25.0*12]

#load  qx   qy
uniformLoad = (0.0, 10/12/1000)  #10 lbs/ft to kips/in.

z, u, v, ϕ, beamProperties = Beam.solve(memberDefinitions, sectionProperties, materialProperties, loadLocation, springStiffness, endBoundaryConditions, supports, uniformLoad)

B = InternalForces.bimoment(z, beamProperties.dm, ϕ, beamProperties.E, beamProperties.Cw)

Bmax = maximum(B)

B_D100 = 7.44   #from Bob Glauz's ballot example, see /testfiles/InternalForcesTest5/495A_Examples_v2.pdf

Error= abs((Bmax - B_D100)/ B_D100)

#accept 1% error from numerical solution
@test Error <= 0.01


text:
#!/usr/bin/env julia
#
# Jerzy Wroczyński (nr. indeksu 250075)
# Obliczenia Naukowe, Laboratorium
# Lista 2, Zadanie 4.
#

# imports
using Polynomials
using LinearAlgebra

# the polynomial coefficients
coeffs = [1.0, -210.0, 20615.0,-1256850.0, 53327946.0, -1672280820.0, 40171771630.0, -756111184500.0, 11310276995381.0, -135585182899530.0, 1307535010540395.0, -10142299865511450.0, 63030812099294896.0, -311333643161390640.0, 1206647803780373360.0, -3599979517947607200.0,
 8037811822645051776.0, -12870931245150988800.0,
 13803759753640704000.0, -8752948036761600000.0,
 2432902008176640000.0]

# if given any arguments run the program with a slight modification
if length(ARGS) > 0
    coeffs[2] = coeffs[2] - 2^(-23)
end

# generate the polynomial
reverse!(coeffs)
P = Polynomial(coeffs)

# this is the canonical form of the polynomial
function p(x)
    return (x - 20)*(x - 19)*(x - 18)*(x - 17)*(x - 16)*(x - 15)*(x - 14)*(x - 13)*(x - 12)*(x - 11)*(x - 10)*(x - 9)*(x - 8)*(x - 7)*(x - 6)*(x - 5)*(x - 4)*(x - 3)*(x - 2)*(x - 1)
end

# calculate the roots of the polynomial (given its natural form)
rs = roots(P)

# print out the results
for k in 1:20

    zₖ = rs[k]
    print(
        k, " & ",
        zₖ, " & ",
        norm(P(zₖ)), " & ",
        norm(p(zₖ)), " & ",
        norm(zₖ - k), "\\\\\n"
    )

end

max_stars_repo_name: jerry-sky/academic-notebook

max_issues_repo_name: jerry-sky/academic-notebook

max_forks_repo_name: jerry-sky/academic-notebook


text:
# ---
# title: 381. Insert Delete GetRandom O(1) - Duplicates allowed
# id: problem381
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: Array, Hash Table, Design
# link: <https://leetcode.com/problems/insert-delete-getrandom-o1-duplicates-allowed/description/>
# hidden: true
# ---
#
# Design a data structure that supports all following operations in _average_
# **O(1)** time.
#
# **Note: Duplicate elements are allowed.**
#
#   1. `insert(val)`: Inserts an item val to the collection.
#   2. `remove(val)`: Removes an item val from the collection if present.
#   3. `getRandom`: Returns a random element from current collection of elements. The probability of each element being returned is **linearly related** to the number of same value the collection contains.
#
# **Example:**
#
#
#
#     // Init an empty collection.
#     RandomizedCollection collection = new RandomizedCollection();
#
#     // Inserts 1 to the collection. Returns true as the collection did not contain 1.
#     collection.insert(1);
#
#     // Inserts another 1 to the collection. Returns false as the collection contained 1. Collection now contains [1,1].
#     collection.insert(1);
#
#     // Inserts 2 to the collection, returns true. Collection now contains [1,1,2].
#     collection.insert(2);
#
#     // getRandom should return 1 with the probability 2/3, and returns 2 with the probability 1/3.
#     collection.getRandom();
#
#     // Removes 1 from the collection, returns true. Collection now contains [1,2].
#     collection.remove(1);
#
#     // getRandom should return 1 and 2 both equally likely.
#     collection.getRandom();
#
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
import Base./, Base.%

function S_poly(A::Multinomial, B::Multinomial)
    x = A[end].indeterminate
    tmp = [max(A[end].exponent[i], B[end].exponent[i]) for i in 1:x.degree]
    return B[end].coefficient*x^(tmp-A[end].exponent)*A -
        A[end].coefficient*x^(tmp-B[end].exponent)*B
end

"""
Does monomial a divide monomial b?
"""
function divides(a::Monomial, b::Monomial)
   !isapprox(a.coefficient, 0) &&
       all(a.exponent .<= b.exponent) &&
       a.indeterminate == b.indeterminate
end

"""
Quotient of b and a. Assumes divisibility
"""
function /(b::Monomial, a::Monomial)
    Monomial(b.indeterminate,
             b.exponent-a.exponent,
             b.coefficient/a.coefficient)
end

""" function %(A::Multinomial, G::Array{Multinomial,1})

Divide A by the array of multinomials, G returning the remainder.
"""
function %(A::Multinomial, G::Array{Multinomial,1})
    h = deepcopy(A)
    i = 1
    while i <= length(G)
        g = G[i]
        # find first j, if any, such that g[end] divides h[end-j]
        j = 0
        while j < length(h)
            if divides(g[end], h[end-j])
                break
            else
                j += 1
            end
        end
        if j < length(h)
            h = h - (h[end-j]/g[end])*g
            i = 1
        else
            i += 1
        end
    end
    return h
end

""" function iszero(A::Multinomial)

Is multinomial A zero?
"""
function iszero(A::Multinomial)
    length(A)==0 || (length(A)==1 && isapprox(A[1].coefficient, 0))
end

"""
function buchberger(F::Array{Multinomial,1})

Buchberger's algorithm for finding a Groebner basis for the ideal generated by the system of multinomials, F.

This is a straightforward, unoptimized implementation of the algorithm as described by Buchberger himself in [Gröbner Bases: A Short Introduction for Systems Theorists](http://link.springer.com/chapter/10.1007/3-540-45654-6_1)[PDF](http://people.reed.edu/~davidp/pcmi/buchberger.pdf) as follows:

Start with G = F .
For any pair of polynomials f1 , f2 in G:
  Compute the S✁polynomial of f1, f2
  and reduce it to a reduced form h w.r.t. G.
  If h = 0, consider the next pair.
  If h is not 0, add h to G and iterate.

Note that G never shrinks. In the worst case it can become intractably large. Moreover, the final G may contain many redundant multinomials. Neither of these problems aare faced here.
"""
function buchberger(F::Array{Multinomial,1};
                    increment::Int=1000,
                    maxiterations::Int=1000,
                    # for re-entry:
                    p1::Int=1,
                    p2::Int=2)
    ## Initialization ##
    # Expandable storage for Groebner basis
    nG = max(length(F), increment)
    iG = length(F)
    G = Array(Multinomial, nG)
    # Initialize G = F
    for i in 1:iG G[i] = deepcopy(F[i]) end
    ## Main loop ##
    iteration = 1
    # Continue while new pairs are available and max iterations
    # are not exceeded.
    while p2 <= iG && iteration <= maxiterations
        f1 = G[p1]
        f2 = G[p2]
        # Compute their S-polynomial
        S = S_poly(f1,f2)
        # Find S mod G
        h = S%G[1:iG]
        if !iszero(h)
            iG += 1
            if iG > length(G)
                G = vcat(G,Array(Monomial, increment))
            end
            G[iG] = h
        end
        # bookkeeping
        iteration += 1
        # next pair: note p2 increments more slowly than p1
        # so that all pairs p1 < p2 <= k are covered before
        # pairs p1 < p2 = k+1 are reached.
        p1, p2 = p1+1 < p2 ? (p1+1,p2) : (1, p2+1)
    end
    # Return G, p1, p2
    return G, iG, p1, p2
end

""" normalize!(G!::Array{Multinomial,1})

Convert all leading coefficients in  Groebner basis, G, to 1 modififying G in place.
"""
function normalize!(G::Array{Multinomial,1})
    for g in G
        d = g[end].coefficient
        for m in g
            m.coefficient /= d
        end
    end
    nothing
end

function mark_nonredundant(G::Array{Multinomial,1})
    marks = trues(length(G))
    for i in eachindex(G)
        m = G[i][end]
        for j in eachindex(G)
            if j==i continue end
            if divides(G[j][end], m)
                marks[i] = false
                break
            end
        end
    end
    return marks
end

function reduce_Groebner(G::Array{Multinomial,1})
    Gtmp = deepcopy(G)
    normalize!(Gtmp)
    Gtmp = Gtmp[mark_nonredundant(Gtmp)]
    idx = trues(length(Gtmp))
    for i in eachindex(Gtmp)
        idx[i] = false
        Gtmp[i] = Gtmp[i]%Gtmp[idx]
        idx[i] = true
    end
    return Gtmp
end

function demo_mod()
    x = Indeterminate(2)
    P = 1 + x^[1,2]
    println("Example 1\n")
    @show P
    println()
    G = Array(Multinomial,2)
    G[1]= 1 + x^[1,1]
    G[2]= 1 + x^[0,1]
    @show G[1]
    @show G[2]
    println()
    @show P%G
    println("\nExample 2\n")
    P = x^[1,2] + x^[2,1] + x^[0,2]
    @show P
    println()
    G[1] = x^[1,1]-1
    G[2] = x^[0,2]-1
    @show G[1]
    @show G[2]
    println()
    @show P%G
    nothing
end

function demo_alg()
    @show x = Indeterminate(3)
    println()
    @show A = 1 + 2*x^[1,2,3] + 3*x^[3,2,1]
    println()
    @show B = 2 + 5* x^[1,3,2]
    println()
    @show S_poly(A,B)
    F = Array(Multinomial,2)
    F[1] = A; F[2] = B
    println()
    @show F
    println()
    G, iG, p1, p2 = buchberger(F, maxiterations=10)
    @show iG, p1, p2
    println()
    @show G[iG]
    nothing
end

function demo_bb()
    @show x = Indeterminate(2)
    println()
    F = Array(Multinomial,2)
    F[1] = x^[1,1] - x^[1,0]
    F[2] = x^[2,0] - x^[0,1]
    @show F[1]
    @show F[2]
    println("\nRunning buchberger on F. Result:\n")
    G, iG, p1, p2 = buchberger(F, maxiterations=50)
    @show iG, p1, p2
    println()
    for i in 1:iG @show G[i] end
    println("\nReducing the Groebner basis found above.\n")
    Gred = reduce_Groebner(G[1:iG])
    for i in eachindex(Gred) @show Gred[i] end
    nothing
end

function demo_reduce()
    @show x = Indeterminate(2)
    println()
    G = Array(Multinomial,5)
    G[1] = x^[3,0] - 2*x^[1,1]
    G[2] = x^[2,1] - 2*x^[0,2] + x^[1,0]
    G[3] = [x^[2,0]]
    G[4] = [x^[1,1]]
    G[5] = x^[0,2] - (1//2)*x^[1,0]
    @show G[1]
    @show G[2]
    @show G[3]
    @show G[4]
    @show G[5]
    println("\nReducing the Groebner basis shown above.\n")
    Gred = reduce_Groebner(G[1:5])
    for i in eachindex(Gred) @show Gred[i] end
    nothing
end

function bug1()
    x=Indeterminate(2)
    @show -x^[0,2]
    println()
    @show x^[0,4]-x^[0,2]
    println()
    @show -(x^[0,4],x^[0,2])
    nothing
end


text:
using Documenter
using DiffEqBase

makedocs(
    sitename = "DiffEqBase",
    format = :html,
    modules = [DiffEqBase],
    pages = [
        "index.md",
        "API" => [
            "Overview" => "api/overview.md",
            "api/functions.md",
            "api/problems.md",
        ],
    ],
)

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


text:
# Unit Tests of functions within src/SwCore.jl

using NCDatasets, Test

include("../src/SwCore.jl")
include("../src/Csw.jl")

using .sw_core_mod, .program

current_directory = @__DIR__

input_file_name = relpath("test/data/inputs/c_sw_12x24.nc", ".")

example_dataset = NCDataset(input_file_name)

state = State(example_dataset)

function subroutine_tests()
    @testset "c_sw subroutine tests" begin

        # @test c_sw!(state, 1)

        # @test divergence_corner!(state, 1)

        # @test d2a2c_vect!(state, 1)

        # @test edge_interpolate4!(state.ua, state.dxa)

        # @test fill2_4corners!(state.delp[:,:,1], state.pt[:,:,1], 1, true, true, true, true, state.npx, state.npy)

        # @test fill_4corners!(state.w[:,:,1], 1, true, true, true, true, state.npx, state.npy)

    end
end

function load_state_tests()

   # Test various variables, attributes and dimensions from .nc file are properly loaded into Julia Struct
    @test state.isd == -2

end

function write_and_print_tests()
    @testset "write_state test"

        # create a new .nc file, write some data to it, test that the data was writen correctly, delete the file

    @testset "print_state test"

        # create a new .test file, print some data to it, test that the data was writen correctly, delete the file

end

@testset "All tests" begin
    subroutine_tests()
    # load_state_tests()
    # write_print_tests()
end


text:
#
# Copyright (c) 2021 Tobias Thummerer, Lars Mikelsons, Josef Kircher
# Licensed under the MIT license. See LICENSE file in the project root for details.
#

tool = ENV["EXPORTINGTOOL"]
pathToFMU = joinpath(dirname(@__FILE__), "..", "model", tool, "SpringFrictionPendulum1D.fmu")

myFMU = fmiLoad(pathToFMU)

@test fmiGetVersion(myFMU) == "2.0"
@test fmiGetTypesPlatform(myFMU) == "default"

@test fmiGetModelName(myFMU) == "SpringFrictionPendulum1D"
@test fmiGetVariableNamingConvention(myFMU) == "structured"
@test fmiIsCoSimulation(myFMU) == true
@test fmiIsModelExchange(myFMU) == true

if tool == "Dymola/2020x"
    @test fmiGetGUID(myFMU) == "{a357de9f-764f-41f1-81c4-b2d344c98a37}"
    @test fmiGetGenerationTool(myFMU) == "Dymola Version 2020x (64-bit), 2019-10-10"
    @test fmiGetGenerationDateAndTime(myFMU) == "2021-05-06T08:23:56Z"
    @test fmiGetNumberOfEventIndicators(myFMU) == 24
    @test fmiCanGetSetState(myFMU) == true
    @test fmiCanSerializeFMUstate(myFMU) == true
    @test fmiProvidesDirectionalDerivative(myFMU) == true
elseif tool == "OpenModelica/v1.17.0"
    @test fmiGetGUID(myFMU) == "{8584aa5b-179e-44ed-9ba6-d557ed34541e}"
    @test fmiGetGenerationTool(myFMU) == "OpenModelica Compiler OMCompiler v1.17.0"
    @test fmiGetGenerationDateAndTime(myFMU) == "2021-06-21T11:48:49Z"
    @test fmiGetNumberOfEventIndicators(myFMU) == 14
    @test fmiCanGetSetState(myFMU) == false
    @test fmiCanSerializeFMUstate(myFMU) == false
    @test fmiProvidesDirectionalDerivative(myFMU) == false
else
    @warn "Unknown exporting tool `$tool`. Skipping model description tests."
end

fmiUnload(myFMU)


text:
# This file was generated, do not modify it. # hide
Xnew = (LStat = Xnew.LStat, LStat2 = Xnew.LStat.^2)

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

savefig("assets/literate/ISL-lab-3-lreg.svg") # hide


text:
#Tests des packages
include("../ALAS/include.jl")

using MPCCmod
using OutputRelaxationmod
#package pour plot
using PyPlot

using NLPModels
using JuMP

#Est-ce que l'on veut les figures :
plot=false

ex1=JuMP.Model()
JuMP.@variable(ex1,x[1:2],start=1.0)
JuMP.@NLobjective(ex1,Min,x[1]-x[2])
JuMP.@constraint(ex1,1-x[2]>=0)
ex1=MathProgNLPModel(ex1)
G=JuMP.Model()
JuMP.@variable(G,x[1:2],start=1.0)
JuMP.@constraint(G,x[1]>=0)
JuMP.@NLobjective(G,Min,0.0)
G=MathProgNLPModel(G)
H=JuMP.Model()
JuMP.@variable(H,x[1:2],start=1.0)
JuMP.@constraint(H,x[2]>=0)
JuMP.@NLobjective(H,Min,0.0)
H=MathProgNLPModel(H)

ex1= MPCCmod.MPCC(ex1,G,H)

ex2=JuMP.Model()
JuMP.@variable(ex2,x[1:2],start=-1.0)
JuMP.@NLobjective(ex2,Min,0.5*((x[1]-1)^2+(x[2]-1)^2))
ex2=MathProgNLPModel(ex2)
G=JuMP.Model()
JuMP.@variable(G,x[1:2],start=1.0)
JuMP.@constraint(G,x[1]>=0)
JuMP.@NLobjective(G,Min,0.0)
G=MathProgNLPModel(G)
H=JuMP.Model()
JuMP.@variable(H,x[1:2],start=1.0)
JuMP.@constraint(H,x[2]>=0)
JuMP.@NLobjective(H,Min,0.0)
H=MathProgNLPModel(H)

ex2= MPCCmod.MPCC(ex2,G,H)

# Exemple 3 :
# minimize x1+x2-x3
# s.t.     0<=x1 _|_ x2>=0
#          -4x1+x3<=0
#          -4x2+x3<=0

ex3=JuMP.Model()
JuMP.@variable(ex3,x[1:3],start=1.0)
JuMP.@NLobjective(ex3,Min,x[1]+x[2]-x[3])
JuMP.@NLconstraint(ex3,c1,-4*x[1]+x[3]<=0)
JuMP.@NLconstraint(ex3,c2,-4*x[2]+x[3]<=0)
ex3=MathProgNLPModel(ex3)
G=JuMP.Model()
JuMP.@variable(G,x[1:3],start=1.0)
JuMP.@constraint(G,x[1]>=0)
JuMP.@NLobjective(G,Min,0.0)
G=MathProgNLPModel(G)
H=JuMP.Model()
JuMP.@variable(H,x[1:3],start=1.0)
JuMP.@constraint(H,x[2]>=0)
JuMP.@NLobjective(H,Min,0.0)
H=MathProgNLPModel(H)

ex3= MPCCmod.MPCC(ex3,G,H)


text:
@testset "types" begin
    @testset "Mag" begin
        @test Mag() isa Mag
        @test Mag(Mag()) isa Mag
        @test Mag(Arf()) isa Mag
    end

    @testset "Arf" begin
        @test Arf() isa Arf
        @test precision(Arf(prec = 80)) == 80
        @test Arf(UInt64(1)) isa Arf
        @test precision(Arf(UInt64(1), prec = 80)) == 80
        @test Arf(1) isa Arf
        @test precision(Arf(1, prec = 80)) == 80
    end

    @testset "Arb" begin
        @test Arb() isa Arb
        @test precision(Arb(prec = 80)) == 80
    end

    @testset "Acb" begin
        @test Acb() isa Acb
        @test precision(Acb(prec = 80)) == 80
    end

    @testset "ArbPoly" begin
        @test ArbPoly() isa ArbPoly
        @test precision(ArbPoly(prec = 80)) == 80
    end

    @testset "ArbSeries" begin
        @test ArbSeries(3) isa ArbSeries
        @test precision(ArbSeries(3, prec = 80)) == 80
        @test ArbSeries(3).degree == 3
    end

    @testset "AcbPoly" begin
        @test AcbPoly() isa AcbPoly
        @test precision(AcbPoly(prec = 80)) == 80
    end

    @testset "AcbSeries" begin
        @test AcbSeries(3) isa AcbSeries
        @test precision(AcbSeries(3, prec = 80)) == 80
        @test AcbSeries(3).degree == 3
    end
end