text:
export divides

function Base.round(p::APL; args...)
    # round(0.1) is zero so we cannot use `mapcoefficientsnz`
    return mapcoefficients(p) do term
        round(term; args...)
    end
end

function Base.div(p::APL, α::Number, args...)
    return mapcoefficients(p) do term
        div(term, α, args...)
    end
end

"""
    divides(t1::AbstractTermLike, t2::AbstractTermLike)

Returns whether the monomial of t1 divides the monomial of t2.

### Examples

Calling `divides(2x^2y, 3xy)` should return false because `x^2y` does not divide `xy` since `x` has a degree 2 in `x^2y` which is greater than the degree of `x` on `xy`.
However, calling `divides(3xy, 2x^2y)` should return true.
"""
function divides(t1::AbstractTermLike, t2::AbstractTermLike)
    divides(monomial(t1), monomial(t2))
end
divides(t1::AbstractVariable, t2::AbstractVariable) = t1 == t2

Base.gcd(m1::AbstractMonomialLike, m2::AbstractMonomialLike) = mapexponents(min, m1, m2)
Base.lcm(m1::AbstractMonomialLike, m2::AbstractMonomialLike) = mapexponents(max, m1, m2)

struct Field end
struct UniqueFactorizationDomain end
const UFD = UniqueFactorizationDomain

algebraic_structure(::Type{<:Union{Rational, AbstractFloat}}, ::Type{<:Union{Rational, AbstractFloat}}) = Field()
algebraic_structure(::Type{<:Union{Rational, AbstractFloat}}, ::Type) = Field()
algebraic_structure(::Type, ::Type{<:Union{Rational, AbstractFloat}}) = Field()
algebraic_structure(::Type, ::Type) = UFD()

# _div(a, b) assumes that b divides a
_div(::Field, a, b) = a / b
_div(::UFD, a, b) = div(a, b)
_div(a, b) = _div(algebraic_structure(typeof(a), typeof(b)), a, b)
_div(m1::AbstractMonomialLike, m2::AbstractMonomialLike) = mapexponents(-, m1, m2)
function _div(t::AbstractTerm, m::AbstractMonomial)
    term(coefficient(t), _div(monomial(t), m))
end
function _div(t1::AbstractTermLike, t2::AbstractTermLike)
    term(_div(coefficient(t1), coefficient(t2)), _div(monomial(t1), monomial(t2)))
end
function _div(f::APL, g::APL)
    lt = leadingterm(g)
    rf = MA.copy_if_mutable(f)
    rg = removeleadingterm(g)
    q = zero(rf)
    while !iszero(rf)
        ltf = leadingterm(rf)
        if !divides(lt, ltf)
            # In floating point arithmetics, it may happen
            # that `rf` is not zero even if it cannot be reduced further.
            # As `_div` assumes that `g` divides `f`, we know that
            # `rf` is approximately zero anyway.
            break
        end
        qt = _div(ltf, lt)
        q = MA.add!!(q, qt)
        rf = MA.operate!!(removeleadingterm, rf)
        rf = MA.operate!!(MA.sub_mul, rf, qt, rg)
    end
    return q
end

Base.div(f::APL, g::Union{APL, AbstractVector{<:APL}}; kwargs...) = divrem(f, g; kwargs...)[1]
Base.rem(f::APL, g::Union{APL, AbstractVector{<:APL}}; kwargs...) = divrem(f, g; kwargs...)[2]

function pseudo_divrem(f::APL{S}, g::APL{T}, algo) where {S,T}
    return _pseudo_divrem(algebraic_structure(S, T), f, g, algo)
end

function _pseudo_divrem(::Field, f::APL, g::APL, algo)
    q, r = divrem(f, g)
    return one(q), q, r
end

function _pseudo_divrem(::UFD, f::APL, g::APL, algo)
    ltg = leadingterm(g)
    rg = removeleadingterm(g)
    ltf = leadingterm(f)
    if iszero(f) || !divides(monomial(ltg), ltf)
        return one(f), zero(f), zero(f)
    else
        st = constantterm(coefficient(ltg), f)
        new_f = st * removeleadingterm(f)
        qt = term(coefficient(ltf), _div(monomial(ltf), monomial(ltg)))
        new_g = qt * rg
        # Check with `::` that we don't have any type unstability on this variable.
        return convert(typeof(f), st), convert(typeof(f), qt), (new_f - new_g)::typeof(f)
    end
end

function pseudo_rem(f::APL{S}, g::APL{T}, algo) where {S,T}
    return _pseudo_rem(algebraic_structure(S, T), f, g, algo)
end

function _pseudo_rem(::Field, f::APL, g::APL, algo)
    return true, rem(f, g)
end

function _pseudo_rem(::UFD, f::APL, g::APL, algo)
    ltg = leadingterm(g)
    rg = removeleadingterm(g)
    ltf = leadingterm(f)
    if !divides(monomial(ltg), ltf)
        return false, f
    end
    while !iszero(f) && divides(monomial(ltg), ltf)
        new_f = constantterm(coefficient(ltg), f) * removeleadingterm(f)
        new_g = term(coefficient(ltf), _div(monomial(ltf), monomial(ltg))) * rg
        # Check with `::` that we don't have any type unstability on this variable.
        f = (new_f - new_g)::typeof(f)
        if algo.primitive_rem
            f = primitive_part(f, algo)::typeof(f)
        end
        if algo.skip_last && maxdegree(f) == maxdegree(g)
            break
        end
        ltf = leadingterm(f)
    end
    return true, f
end

function MA.promote_operation(
    ::typeof(pseudo_rem),
    ::Type{P},
    ::Type{Q},
) where {T,S,P<:APL{T},Q<:APL{S}}
    return _promote_operation(algebraic_structure(T, S), pseudo_rem, P, Q)
end
function _promote_operation(
    ::Field,
    ::typeof(pseudo_rem),
    ::Type{P},
    ::Type{Q},
) where {P<:APL,Q<:APL}
    return MA.promote_operation(rem, P, Q)
end
function _promote_operation(
    ::UFD,
    ::typeof(pseudo_rem),
    ::Type{P},
    ::Type{Q},
) where {T,S,P<:APL{T},Q<:APL{S}}
    U1 = MA.promote_operation(*, S, T)
    U2 = MA.promote_operation(*, T, S)
    # `promote_type(P, Q)` is needed for TypedPolynomials in case they use different variables
    return polynomialtype(promote_type(P, Q), MA.promote_operation(-, U1, U2))
end

function MA.promote_operation(::Union{typeof(div), typeof(rem)}, ::Type{P}, g::Type{Q}) where {T,S,P<:APL{T},Q<:APL{S}}
    U = MA.promote_operation(/, T, S)
    # `promote_type(P, Q)` is needed for TypedPolynomials in case they use different variables
    return polynomialtype(promote_type(P, Q), MA.promote_operation(-, U, U))
end
function Base.divrem(f::APL{T}, g::APL{S}; kwargs...) where {T, S}
    rf = convert(MA.promote_operation(div, typeof(f), typeof(g)), MA.copy_if_mutable(f))
    q = zero(rf)
    r = zero(rf)
    lt = leadingterm(g)
    rg = removeleadingterm(g)
    lm = monomial(lt)
    while !iszero(rf)
        ltf = leadingterm(rf)
        if isapproxzero(ltf; kwargs...)
            rf = MA.operate!!(removeleadingterm, rf)
        elseif divides(lm, ltf)
            qt = _div(ltf, lt)
            q = MA.add!!(q, qt)
            rf = MA.operate!!(removeleadingterm, rf)
            rf = MA.operate!!(MA.sub_mul, rf, qt, rg)
        elseif lm > monomial(ltf)
            # Since the monomials are sorted in decreasing order,
            # lm is larger than all of them hence it cannot divide any of them
            r = MA.add!!(r, rf)
            break
        else
            r = MA.add!!(r, ltf)
            rf = MA.operate!!(removeleadingterm, rf)
        end
    end
    q, r
end
function Base.divrem(f::APL{T}, g::AbstractVector{<:APL{S}}; kwargs...) where {T, S}
    rf = convert(MA.promote_operation(div, typeof(f), eltype(g)), MA.copy_if_mutable(f))
    r = zero(rf)
    q = similar(g, typeof(rf))
    for i in eachindex(q)
        q[i] = zero(rf)
    end
    lt = leadingterm.(g)
    rg = removeleadingterm.(g)
    lm = monomial.(lt)
    useful = BitSet(eachindex(g))
    while !iszero(rf)
        ltf = leadingterm(rf)
        if isapproxzero(ltf; kwargs...)
            rf = MA.operate!!(removeleadingterm, rf)
            continue
        end
        divisionoccured = false
        for i in useful
            if divides(lm[i], ltf)
                qt = _div(ltf, lt[i])
                q[i] = MA.add!!(q[i], qt)
                rf = MA.operate!!(removeleadingterm, rf)
                rf = MA.operate!!(MA.sub_mul, rf, qt, rg[i])
                divisionoccured = true
                break
            elseif lm[i] > monomial(ltf)
                # Since the monomials are sorted in decreasing order,
                # lm is larger than all of them hence it cannot divide any of them
                delete!(useful, i)
            end
        end
        if !divisionoccured
            if isempty(useful)
                r = MA.add!!(r, rf)
                break
            else
                r = MA.add!!(r, ltf)
                rf = MA.operate!!(removeleadingterm, rf)
            end
        end
    end
    q, r
end


text:
#using MinkowskiReduction
using Test
using Random
using LinearAlgebra

@testset "MinkowskiReduction.jl" begin
    U=[1, 2, 3];V=[-1, 2, 3];W=[3, 0, 4]
    @test all(minkReduce(U,V,W) .≈ ([-2.0, 0.0, 0.0], [0.0, -2.0, 1.0], [-1.0, 2.0, 3.0]))
    @test orthogonalityDefect(minkReduce(U,V,W)...)≈1.0458250331675945
    @test orthogonalityDefect(U,V,W)==4.375
    fcc = [[1,1,0],[1,0,1],[0,1,1]]
    @test orthogonalityDefect(fcc...)≈1.4142135623730954
    m = DeviousMat(26) # Largest size that doesn't overflow
    @test det(hcat(minkReduce(m[:,1],m[:,2],m[:,3])...))==1
    @test all(MinkowskiReduction.shortenW_in_UVW(U,V,W) .≈ ([-2.0, 0.0, 0.0], [1.0, 2.0, 3.0], [-2.0, -2.0, 1.0]))
    @test (m = DeviousMat(26); all(minkReduce(m[:,1],m[:,2],m[:,3],true) .≈ ([-1.0, 0.0, 0.0], [0.0, 0.0, -1.0], [0.0, -1.0, 0.0], 15)))
    @test (m = DeviousMat(20); all(minkReduce(m[:,1],m[:,2],m[:,3],true) .≈ ([0.0, 0.0, -1.0], [1.0, 0.0, 0.0], [0.0, 1.0, 0.0], 12)))
    @test (Random.seed!(1234); MinkowskiReduction.RandLowerTri(35)==[1 0;-28 1])
    @test (MinkowskiReduction.FibonacciMat(45)==[1836311903 2971215073; 1134903170 1836311903])
    @test (Random.seed!(1234);MinkowskiReduction.RandUnimodMat(4)==[-1 -1; 0 -1])
    @test (Random.seed!(1234);MinkowskiReduction.RandUnimodMat(2)==[0 1; -1 -1])
    @test !(U = [1, 2, 9]; V = [-3, 4, 3]; W = [3, 0, 4]; isMinkReduced(U,V,W))
    @test !(U = [1, 2, 3]; V = [-3, 4, 3]; W = [3, 0, 4]; isMinkReduced(U,V,W))
    @test !(U = [1, 0, 0]; V = [-1, 0, 1]; W = [3, 0, 4]; isMinkReduced(U,V,W))
    @test !(U = [1, 2, 3]; V = [-1, 2, 3]; W = [3, 0, 4]; isMinkReduced(U,V,W))
    @test !(U = [1, 0, 0]; V = [0,-4,-3]; W = [-3, 0, 4]; isMinkReduced(U,V,W))
    @test !(U = [1, 0, 0]; V = [0, 4, 3]; W = [3, 0, 4]; isMinkReduced(U,V,W))
    @test !(U = [0, 4, 3]; V = [-5, 0, 0]; W = [3, 0, 4]; isMinkReduced(U,V,W))
    @test !(U = [1,0,0]; V = [0,1,0]; W = [0,1,-2]; isMinkReduced(U,V,W))
    @test !(U = [1,0,-2]; V = [-1,1,-2]; W = [0,0,3]; isMinkReduced(U,V,W))
    @test !(δ=.50;U = [1,0,0]; V = [-.4,1,0]; W = [-δ,-δ,5]; isMinkReduced(U,V,W))
    @test !(U = [2,0,0]; V = [-.9,2,0]; W = [-1,-1,10]; isMinkReduced(U,V,W))
    @test !(U = [2,0,0]; V = -[-.9,2,0]; W = [-1,-1,10]; isMinkReduced(U,V,W))
    @test !(U = -[2,0,0]; V = -[-.9,2,0]; W = [-1,-1,10]; isMinkReduced(U,V,W))
    @test !(U = -[2,0,0]; V = [-.9,2,0]; W = [-1,-1,10]; isMinkReduced(U,V,W))
    @test (U = -[2,0,0]; V = [0,2,0]; W = [-1,-1,10]; isMinkReduced(U,V,W))
    U = -[2,0,0]; V = [0,2,0]; W=U+V # Linearly dependedent basis
    @test_throws ErrorException minkReduce(U,V,W)
    U = -[2.,0,0]; V = [0,2.,0]; W=U+V+[0,0,1e-170]
    @test_throws ErrorException minkReduce(U,V,W)
    @test (U = -[2.,0,0]; V = [0,2.,0]; W=U+V+[0,0,1e-150];minkReduce(U,V,W)==
    ([0.0, 0.0, 1.0e-150], [-2.0, -0.0, -0.0], [0.0, 2.0, 0.0]))
    U = -[2^63-1,0,0]; V = [0,2^63-1,0]; W=U+V+[0,0,1e-150];
    @test_throws InexactError minkReduce(U,V,W)
end

ext: jl


text:
function vblast{T}(H::Array{T,2},mu=Inf)
# (W,P,B) = vblast(H)
#  Returns a VBLAST matrix decomposition of H such that
#  H = pinv(W)*B*P' or B = W*H*P.  Here P is a permutation matrix, B is
#  lower triangular with ones on the diagonal, and W has orthogonal rows.
# (W,P,B) = vblast(H,mu)
#  If an SNR argument mu is passed in, a regularized ("MMSE")
#  decomposition is done, with the result that W will no longer have
#  orthogonal rows and B is no longer lower triangular.
#
# Example: For decoding with N receivers
#  N=4;H = rand(N,N); include("vblast.jl");(W,P,B) = vblast(H)
# Example: For MMSE decoding with N receivers mu=rho/N
#  (W,P,B) = vblast(H,rho/size(H,1));

(K,N) = size(H);

W = zeros(T,K,K);
P = zeros(Int,N,N);
Hp = H;
k = zeros(Int,K,1);
indexH = [1:N;];
for i=1:N
    G = inv(Hp'*Hp + eye(N-i+1)/mu)*Hp'; # Do pseudoinverse
    gn = diag(G*G');
    (mx,mi) = findmin(abs(gn))
    k[i] = indexH[mi];
    W[i,:] = G[mi,:];
    P[k[i],i]=1;
    indexH = setdiff(indexH,k[i]);
    Hp = H[:,indexH];
end
B = W*H*P  # Calculating B may not be required in some real-world scenarios
return (W,P,B)
end


text:
using Documenter
include("../src/SugarKelp.jl")
Documenter.makedocs(
    root = "docs",
    source = "src",
    build = "build",
    clean = true,
    doctest = true,
    modules = Module[SugarKelp],
    repo = "https://github.com/jagoosw/SugarKelp.jl",
    highlightsig = true,
    sitename = "SugarKelp.jl Documentation",
    expandfirst = [],
    pages = [
        "Index" => "index.md",
        "Function Documentation" => "functions.md"
    ],
    format = Documenter.HTML(prettyurls = false)
)

deploydocs(
    repo = "github.com/jagoosw/SugarKelp.jl.git",
    devbranch = "main"
)


text:
using ForwardDiff
using Plots
using ReverseDiff
using DifferentialEquations
using Flux
using Statistics
# gr()
# pyplot()
plotly()

# Define SHO Hamiltonian Function
H(q, p) = 0.5*q^2 + 0.5*p^2

function dH(q, p)
    dHdq, dHdp = ForwardDiff.gradient(x->H(x[1], x[2]), [q, p])
    return dHdq, dHdp
end

function dH_analytic(q, p)
    return q, p
end

function H_vec_field(q, p)
    dHdq, dHdp = dH(q, p)
    q̇ = dHdp
    ṗ = - dHdq

    return q̇, ṗ
end

# Visualize the H, it's level curves, and the flow as
# as tangent vectors to level curves
# i.e. governing idea is H = const
qs = -1.5:0.1:1.5
ps = -1.5:0.1:1.5

Q = []
P = []
Q̇ = []
Ṗ = []

for q ∈ qs[1:2:end], p ∈ ps[1:2:end]
    q̇, ṗ = H_vec_field(q, p)
    push!(Q, q)
    push!(P, p)
    push!(Q̇, q̇)
    push!(Ṗ, ṗ)
end

p1 = contour(qs, ps, H, xlabel="q", ylabel="p", aspect_ratio=1.0, fill=true, colorbar_title="H", size=(800,800), c=:viridis)
α=0.1
quiver!(p1, Q, P, quiver=(α .* Q̇, α .* Ṗ), aspect_ratio=1.0, color=:white)
xlims!(p1, -1.5, 1.5)
ylims!(p1, -1.5, 1.5)

Ham(q,p, param) = q^2 + p^2
p0=1.0
q0=0.0
tspan=(0.0, 10.0)
prob = HamiltonianProblem(Ham, p0, q0, tspan, dt=0.01)
sol = solve(prob, SymplecticEuler())

plot!(p1, sol[1,:], sol[2, :], color=:green, label="q₀=$(q0), p₀=$(p0)")

p0=1.5
q0=0.0
tspan=(0.0, 10.0)
prob = HamiltonianProblem(Ham, p0, q0, tspan, dt=0.01)
sol = solve(prob, SymplecticEuler())

plot!(p1, sol[1,:], sol[2, :], color=:cyan, label="q₀=$(q0), p₀=$(p0)")

# Let's set up the Neural Network to approximate the Hamiltonian
# generate training data with H(q,p) <= 1
Npoints = 500
r_train = sqrt.(rand(1, Npoints))
φ_train = 2π.*rand(1, Npoints)
data = vcat(r_train .* cos.(φ_train), r_train .* sin.(φ_train))

# set up target data: here we use analytic result for time derivatives
target = zero(data)
# i.e. q̇= dHdp, ṗ = -dHdq
target[1,:] .= data[2, :]
target[2,:] .= -data[1, :]

p2 =  scatter!(p1, data[1,:], data[2,:], c=:green, label="training points")
p3 = plot(Q, P, H.(Q, P), st = :surface, xlabel = "q", ylabel = "p", zlabel = "H")

# -------------------------------------------------------------------------------------------------------
# create model for HNN

# define the struct
struct HNN{M, R, P}
    model::M  # the internal NN
    re::R   # for recreating internal NN
    p::P  # for holding current params of internal NN

    # define the constructor
    function HNN(model)
        p, re = Flux.destructure(model)
        return new{typeof(model), typeof(re), typeof(p)}(model, re, p)
    end
end

# define the trainable paramaters
Flux.trainable(hnn::HNN) = (hnn.p,)

function _hamiltonian_flow(re, p, x)
    dHdX = Flux.gradient(x->sum(re(p)(x)), x)[1]  # note: re(p)(x) == Model(x)
    n = size(x,1) ÷ 2  # i.e. how many p's and q's
    return cat(dHdX[(n+1):2n, :], dHdX[1:n, :], dims=1)
end

# define how to call the HNN on data
(hnn::HNN)(X, p=hnn.p) = _hamiltonian_flow(hnn.re, p, X)

hnn = HNN(
    Chain(Dense(2, 200, relu), Dense(200, 200, relu), Dense(200, 1))
)

dataloader = Flux.Data.DataLoader((data, target), batchsize=500, shuffle=true)

p = hnn.p

opt = ADAM(10^(-3))

loss(x, y, p) = mean((hnn(x, p) .- y) .^ 2)
loss(data, target, p)

callback() = println("Loss Neural Hamiltonian DE = $(loss(data, target, p))")
callback()

test_gs = ReverseDiff.gradient(p -> loss(data, target, p), p)

epochs = 2000
for epoch in 1:epochs
    for (x, y) in dataloader
        gs = ReverseDiff.gradient(p -> loss(x, y, p), p)
        Flux.Optimise.update!(opt, p, gs)
    end
    if epoch % 100 == 1
        callback()
    end
end
callback()

#p3 = plot(Q, P, H.(Q, P), st = :surface, xlabel = "q", ylabel = "p", zlabel = "H")
#p3 = plot(Q, P, HNN.model.(Q, P), st = :surface, xlabel = "q", ylabel = "p", zlabel = "H")


text:
using Documenter, ImageCore

makedocs(modules  = [ImageCore],
         format   = Documenter.Formats.HTML,
         sitename = "ImageCore",
         pages    = ["index.md", "views.md", "traits.md", "reference.md"])

deploydocs(repo   = "github.com/JuliaImages/ImageCore.jl.git",
           julia  = "0.5",
           target = "build",
           deps   = nothing,
           make   = nothing)
#           deps   = Deps.pip("mkdocs", "python-markdown-math"))


text:
# using Distributed
# addprocs(5)

using Caesar, RoME
using RoMEPlotting, Gadfly, KernelDensityEstimatePlotting, Cairo, Fontconfig
using KernelDensityEstimate
using IncrementalInference
using DocStringExtensions
using DelimitedFiles, JLD, DistributedFactorGraphs

include(joinpath(@__DIR__,"slamUtils.jl"))
include(joinpath(@__DIR__,"plotSASUtils.jl"))
include(joinpath(@__DIR__,"expDataUtils.jl"))

function main(expID::String, datastart::Int, dataend::Int, fgap::Int, gps_gap::Int; saswindow::Int=8, trialID::Int=1)
    ## Default parameters
    println("Start Script\n")

    sas_gap = saswindow+fgap;
    tstart = 1_000_000;
    element = 2

    totalposes = dataend-datastart;
    datawindow = collect(datastart:dataend);
    allpaths = getPaths(expID,"dynsas", currloc = "local", trialID = trialID);
    scriptHeader = joinpath(allpaths[3],expID*"_fgap$(fgap)_gpsgap$(gps_gap)_swind$(saswindow)_f$(datastart)t$(dataend)_")

    poses = Dict{Int,Array}();
    nav = Dict{Int,Array}();
    errorinds = [];
    allsasframes = Dict{Int,Array}();

    igt = loadGT(datawindow,expID,allpaths);
    posData = loadNav(datawindow,allpaths)
    dposData = deepcopy(posData)
    i=1
    while i+gps_gap < size(dposData,1)
        cumulativeDrift!(@view(dposData[i:i+gps_gap-1,:]),[0.0;0],[0.2,0.2])
         i+=gps_gap
        if i+gps_gap > size(dposData,1)
            cumulativeDrift!(@view(dposData[i-gps_gap:end,:]),[0.0;0],[0.2,0.2])
        end
    end

    # posPl = Gadfly.plot(layer(x=posData[:,1],y=posData[:,2], Geom.path, Theme(default_color=colorant"green")), layer(x=dposData[:,1],y=dposData[:,2],Geom.path)); posPl |> PDF(scriptHeader*"posGT.pdf")

    chirpFile = joinpath(ENV["HOME"],"data","sas","chirp250.txt");
    cfgFile = joinpath(ENV["HOME"],"data","sas","SAS2D.yaml")
    cfgd=loadConfigFile(cfgFile)

    fg = initfg();
    beacon = :l1
    addVariable!(fg, beacon, Point2 )
    tree = emptyBayesTree();

    pose_counter = 1
    sas_counter = 1
    sas_gap_counter = 0;

    while pose_counter < totalposes-2
        current_symbol = Symbol("x$pose_counter")
        addVariable!(fg, current_symbol, DynPoint2(ut=tstart))
        println("Adding Variable x$pose_counter")
        sas_gap_counter += 1
        tstart += 1_000_000

        #Interpolate velocity from GPS ground truth every gps_gap poses
       if pose_counter == 1 || mod(pose_counter,gps_gap) == 0
           xdotp = posData[pose_counter+1,1] - posData[pose_counter,1];
           ydotp = posData[pose_counter+1,2] - posData[pose_counter,2];
            dpμ = [posData[pose_counter,1];posData[pose_counter,2];xdotp;ydotp];
            dpσ = Matrix(Diagonal([0.1;0.1;0.2;0.2].^2))

            pp = DynPoint2VelocityPrior(MvNormal(dpμ,dpσ))
            addFactor!(fg, [current_symbol;], pp, autoinit=false)
            # initManual!(fg,current_symbol,kde!(rand(MvNormal(dpμ,dpσ),100)))
       end

       #odo factors
       if pose_counter > 1
           xdotp = dposData[pose_counter-1,1] - dposData[pose_counter,1];
           ydotp = dposData[pose_counter-1,2] - dposData[pose_counter,2];
           xdotf = dposData[pose_counter+1,1] - dposData[pose_counter,1];
           ydotf = dposData[pose_counter+1,2] - dposData[pose_counter,2];
           dpμ = [xdotp;ydotp;xdotf-xdotp;xdotf-xdotp];
           dpσ = Matrix(Diagonal([0.2;0.2;0.2;0.2].^2))
           vp = VelPoint2VelPoint2(MvNormal(dpμ,dpσ))
           addFactor!(fg, [Symbol("x$(pose_counter-1)");current_symbol], vp, autoinit=false)
           # initManual!(fg,current_symbol,kde!(rand(MvNormal(dpμ,dpσ),100)))
       end

       if sas_gap_counter >= sas_gap
           println("Trying to add SAS at pose: $pose_counter")
           navframes = collect(pose_counter-saswindow+1:1:pose_counter)
           navtemp = posData[navframes,:];
           navchecked, errorind = sanitycheck_nav(navtemp);
           push!(errorinds,errorind);

           if navchecked
               poses[sas_counter] = [Symbol("x$i") for i in pose_counter-saswindow+1:pose_counter]
               currentstart = datastart+pose_counter-saswindow
               println("Adding SAS at dataframe: $currentstart")
               sasframes = collect(currentstart:1:currentstart+saswindow-1);
               allsasframes[sas_counter] = sasframes;
               waveformData = importdata_waveforms(sasframes,element,datadir=allpaths[1]);
               sas2d = prepareSAS2DFactor(saswindow, waveformData, rangemodel=:Correlator,cfgd=cfgd, chirpFile=chirpFile)
               addFactor!(fg, [beacon;poses[sas_counter]], sas2d, autoinit=false)

               getSolverParams(fg).drawtree = false
               getSolverParams(fg).showtree = false

               if sas_counter > 1
                   tree = solveTree!(fg,tree)
               else
                   tree = solveTree!(fg)
               end

               writeGraphPdf(fg,viewerapp="", engine="neato", filepath=scriptHeader*"fg.pdf")
               drawTree(tree, filepath=scriptHeader*"bt.pdf")

               plotSASDefault(fg,expID, posData,igt,datadir=allpaths[1],savedir=scriptHeader*"$sas_counter.pdf")

               if expID == "dock"
                   l1fit = getKDEMean(getBelief(getVariable(fg,:l1)))
                   meanerror = sqrt((igt[1]-l1fit[1])^2+(igt[2]-l1fit[2])^2)

                   l1max = getKDEMax(getBelief(getVariable(fg,:l1)))
                   maxerror = sqrt((igt[1]-l1max[1])^2+(igt[2]-l1max[2])^2)

                   mykde = getBelief(getVariable(fg,:l1))
                   mynorm = kde!(rand(fit(MvNormal,getVal(fg,:l1)),200))
                   kld = min(abs(KernelDensityEstimate.kld(mynorm,mykde)),abs(KernelDensityEstimate.kld(mykde,mynorm)))
               elseif expID == "drift"
                   l1fit = getKDEMean(getBelief(getVariable(fg,:l1)))
                   meanerror = sqrt((mean(igt[:,1])-l1fit[1])^2+(mean(igt[:,2])-l1fit[2])^2)

                   l1max = getKDEMax(getBelief(getVariable(fg,:l1)))
                   maxerror = sqrt((mean(igt[:,1])-l1max[1])^2+(mean(igt[:,2])-l1max[2])^2)

                   mykde = getBelief(getVariable(fg,:l1))
                   mynorm = kde!(rand(fit(MvNormal,getVal(fg,:l1)),200))
                   kld = min(abs(KernelDensityEstimate.kld(mynorm,mykde)),abs(KernelDensityEstimate.kld(mykde,mynorm)))
               end

               jldname2 = scriptHeader * "solve_$(sas_counter).jld"
               JLD.save(jldname2,"beacon",getVal(fg,:l1),"posData",posData,"dposData", dposData,"gps_gap", gps_gap, "poses",poses,"sasframes", allsasframes, "l1fit",l1fit, "meanerror",meanerror,"l1max",l1max,"maxerror",maxerror,"kld",kld)

               saveDFG(fg,scriptHeader * "fg$(sas_counter)")

               sas_counter +=1
               sas_gap_counter = 0
           end
       end
       pose_counter+=1

       println("Trying next Variable: $pose_counter")
       println("Trying next SAS: $sas_counter")

   end

   return fg
end

fg = main("drift",1550,1580,1,19,trialID=1)


text:
abstract type SaturationGoal end

reached(g::EGraph, goal::Nothing) = false
reached(g::EGraph, goal::SaturationGoal) = false

"""
This goal is reached when the `exprs` list of expressions are in the
same equivalence class.
"""
struct EqualityGoal <: SaturationGoal
    exprs::Vector{Any}
    ids::Vector{EClassId}
    function EqualityGoal(exprs, eclasses)
        @assert length(exprs) == length(eclasses) && length(exprs) != 0
        new(exprs, eclasses)
    end
end

function reached(g::EGraph, goal::EqualityGoal)
    all(x -> in_same_class(g, goal.ids[1], x), @view goal.ids[2:end])
end

"""
Boolean valued function as an arbitrary saturation goal.
User supplied function must take an [`EGraph`](@ref) as the only parameter.
"""
struct FunctionGoal <: SaturationGoal
    fun::Function
end

function reached(g::EGraph, goal::FunctionGoal)::Bool
    fun(g)
end

mutable struct Report
    reason::Union{Symbol, Nothing}
    egraph::EGraph
    iterations::Int
    to::TimerOutput
end

Report() = Report(nothing, EGraph(), 0, TimerOutput())
Report(g::EGraph) = Report(nothing, g, 0, TimerOutput())

# string representation of timedata
function Base.show(io::IO, x::Report)
    g = x.egraph
    println(io, "Equality Saturation Report")
    println(io, "=================")
    println(io, "\tStop Reason: $(x.reason)")
    println(io, "\tIterations: $(x.iterations)")
    # println(io, "\tRules applied: $(g.age)")
    println(io, "\tEGraph Size: $(g.numclasses) eclasses, $(length(g.memo)) nodes")
    print_timer(io, x.to)
end

"""
Configurable Parameters for the equality saturation process.
"""
@with_kw mutable struct SaturationParams
    timeout::Int = 8
    timelimit::Period = Second(-1)
    # default sizeout. TODO make this bytes
    # sizeout::Int = 2^14
    matchlimit::Int = 5000
    eclasslimit::Int = 5000
    enodelimit::Int = 15000
    goal::Union{Nothing, SaturationGoal} = nothing
    stopwhen::Function = ()->false
    scheduler::Type{<:AbstractScheduler} = BackoffScheduler
    schedulerparams::Tuple=()
    threaded::Bool = false
    timer::Bool = true
    printiter::Bool = false
    simterm::Function = similarterm
end

struct Match
    rule::AbstractRule
    # the rhs pattern to instantiate
    pat_to_inst
    # the substitution
    sub::Sub
    # the id the matched the lhs
    id::EClassId
end

const MatchesBuf = Vector{Match}

function cached_ids(g::EGraph, p::AbstractPat)# ::Vector{Int64}
    if isground(p)
        id = lookup_pat(g, p)
        !isnothing(id) && return [id]
    else
        return collect(keys(g.classes))
    end
    return []
end

function cached_ids(g::EGraph, p) # p is a literal
    id = lookup(g, ENodeLiteral(p))
    !isnothing(id) && return [id]
    return []
end

# FIXME
function cached_ids(g::EGraph, p::PatTerm)
    # println("pattern $p, $(p.head)")
    # println("all ids")
    # keys(g.classes) |> println
    # println("cached symbols")
    # cached = get(g.symcache, p.head, Set{Int64}())
    # println("symbols where $(p.head) appears")
    # appears = Set{Int64}()
    # for (id, class) ∈ g.classes
    #     for n ∈ class
    #         if n.head == p.head
    #             push!(appears, id)
    #         end
    #     end
    # end
    # # println(appears)
    # if !(cached == appears)
    #     @show cached
    #     @show appears
    # end

    collect(keys(g.classes))
    # cached
    # get(g.symcache, p.head, [])
end

# function cached_ids(g::EGraph, p::PatLiteral)
#     get(g.symcache, p.val, [])
# end

function (r::SymbolicRule)(g::EGraph, id::EClassId)
    ematch(g, r.ematch_program, id) .|> sub -> Match(r, r.right, sub, id)
end

function (r::DynamicRule)(g::EGraph, id::EClassId)
    ematch(g, r.ematch_program, id) .|> sub -> Match(r, nothing, sub, id)
end

function (r::BidirRule)(g::EGraph, id::EClassId)
    vcat(ematch(g, r.ematch_program_l, id) .|> sub -> Match(r, r.right, sub, id),
        ematch(g, r.ematch_program_r, id) .|> sub -> Match(r, r.left, sub, id))
end

"""
Returns an iterator of `Match`es.
"""
function eqsat_search!(egraph::EGraph, theory::Vector{<:AbstractRule},
    scheduler::AbstractScheduler, report; threaded=false)
    match_groups = Vector{Match}[]
    function pmap(f, xs)
        # const propagation should be able to optimze one of the branch away
        if threaded
            # # try to divide the work evenly between threads without adding much overhead
            # chunks = Threads.nthreads() * 10
            # basesize = max(length(xs) ÷ chunks, 1)
            # ThreadsX.mapi(f, xs; basesize=basesize)
            ThreadsX.map(f, xs)
        else
            map(f, xs)
        end
    end

    inequalities = filter(Base.Fix2(isa, UnequalRule), theory)
    # never skip contradiction checks
    append_time = TimerOutput()
    for rule ∈ inequalities
        @timeit report.to repr(rule) begin
            ids = cached_ids(egraph, rule.left)
            rule_matches = pmap(i -> rule(egraph, i), ids)
            @timeit append_time "appending matches" begin
                append!(match_groups, rule_matches)
            end
        end
    end

    other_rules = filter(theory) do rule
        !(rule isa UnequalRule)
    end
    for rule ∈ other_rules
        @timeit report.to repr(rule) begin
            # don't apply banned rules
            if !cansearch(scheduler, rule)
                # println("skipping banned rule $rule")
                continue
            end
            ids = cached_ids(egraph, rule.left)
            rule_matches = pmap(i -> rule(egraph, i), ids)

            n_matches = isempty(rule_matches) ? 0 : sum(length, rule_matches)
            # @show (rule, n_matches)
            can_yield = inform!(scheduler, rule, n_matches)
            if can_yield
                @timeit append_time "appending matches" begin
                    append!(match_groups, rule_matches)
                end
            end
        end
    end

    # @timeit append_time "appending matches" begin
    #     result = reduce(vcat, match_groups) # this should be more efficient than multiple appends
    # end
    merge!(report.to, append_time, tree_point=["Search"])

    return Iterators.flatten(match_groups)
    # return result
end

function (rule::UnequalRule)(g::EGraph, match::Match; simterm=similarterm)
    lc = match.id
    rinst = instantiate(g, match.pat_to_inst, match.sub, rule; simterm=simterm)
    rc, node = addexpr!(g, rinst)

    if find(g, lc) == find(g, rc)
        @log "Contradiction!" rule
        return :contradiction
    end
    return nothing
end

function (rule::SymbolicRule)(g::EGraph, match::Match; simterm=similarterm)
    rinst = instantiate(g, match.pat_to_inst, match.sub, rule; simterm=simterm)
    rc, node = addexpr!(g, rinst)
    merge!(g, match.id, rc.id)
    return nothing
end

function (rule::DynamicRule)(g::EGraph, match::Match; simterm=similarterm)
    f = rule.rhs_fun
    actual_params = [instantiate(g, PatVar(v, i, alwaystrue), match.sub, rule) for (i, v) in enumerate(rule.patvars)]
    r = f(g[match.id], match.sub, g, actual_params...)
    isnothing(r) && return nothing
    rc, node = addexpr!(g, r)
    merge!(g, match.id, rc.id)
    return nothing
end

function eqsat_apply!(g::EGraph, matches, rep::Report, params::SaturationParams)
    i = 0
    # println.(matches)
    for match ∈ matches
        i += 1

        # if params.eclasslimit > 0 && g.numclasses > params.eclasslimit
        #     @log "E-GRAPH SIZEOUT"
        #     rep.reason = :eclasslimit
        #     return
        # end

        if reached(g, params.goal)
            @log "Goal reached"
            rep.reason = :goalreached
            return
        end

        rule = match.rule
        # println("applying $rule")

        halt_reason = rule(g, match; simterm=params.simterm)
        if (halt_reason !== nothing)
            rep.reason = halt_reason
            return
        end

        # println(rule)
        # println(sub)
        # println(l); println(r)
        # display(egraph.classes); println()
    end
end

import ..@log

"""
Core algorithm of the library: the equality saturation step.
"""
function eqsat_step!(g::EGraph, theory::Vector{<:AbstractRule}, curr_iter,
        scheduler::AbstractScheduler, match_hist::MatchesBuf,
        params::SaturationParams, report)

    instcache = Dict{AbstractRule, Dict{Sub, EClassId}}()

    setiter!(scheduler, curr_iter)

    matches = @timeit report.to "Search" eqsat_search!(g, theory, scheduler, report; threaded=params.threaded)

    # matches = setdiff!(matches, match_hist)

    @timeit report.to "Apply" eqsat_apply!(g, matches, report, params)

    # union!(match_hist, matches)

    if report.reason === nothing && cansaturate(scheduler) && isempty(g.dirty)
        report.reason = :saturated
    end
    @timeit report.to "Rebuild" rebuild!(g)

    return report, g
end

"""
Given an [`EGraph`](@ref) and a collection of rewrite rules,
execute the equality saturation algorithm.
"""
function saturate!(g::EGraph, theory::Vector{<:AbstractRule}, params=SaturationParams())
    curr_iter = 0

    sched = params.scheduler(g, theory, params.schedulerparams...)
    match_hist = MatchesBuf()
    report = Report(g)

    start_time = Dates.now().instant

    !params.timer && disable_timer!(report.to)
    timelimit = params.timelimit > Second(0)

    while true
        curr_iter+=1

        params.printiter && @info("iteration ", curr_iter)

        report, egraph = eqsat_step!(g, theory, curr_iter, sched, match_hist, params, report)

        elapsed = Dates.now().instant - start_time

        if timelimit && params.timelimit <= elapsed
            report.reason = :timelimit
            break
        end

        # report.reason == :matchlimit && break
        if !(report.reason isa Nothing)
            break
        end

        if curr_iter >= params.timeout
            report.reason = :timeout
            break
        end

        if params.eclasslimit > 0 && g.numclasses > params.eclasslimit
            # println(params.eclasslimit)
            report.reason = :eclasslimit
            break
        end

        if reached(g, params.goal)
            report.reason = :goalreached
            break
        end
    end
    report.iterations = curr_iter
    @log report

    return report
end

function areequal(theory::Vector, exprs...; params=SaturationParams())
    g = EGraph(exprs[1])
    areequal(g, theory, exprs...; params=params)
end

function areequal(g::EGraph, t::Vector{<:AbstractRule}, exprs...; params=SaturationParams())
    @log "Checking equality for " exprs
    if length(exprs) == 1; return true end
    # rebuild!(G)

    @log "starting saturation"

    n = length(exprs)
    ids = Vector{EClassId}(undef, n)
    nodes = Vector{AbstractENode}(undef, n)
    for i ∈ 1:n
        ec, node = addexpr!(g, exprs[i])
        ids[i] = ec.id
        nodes[i] = node
    end

    goal = EqualityGoal(collect(exprs), ids)

    # alleq = () -> (all(x -> in_same_set(G.uf, ids[1], x), ids[2:end]))

    params.goal = goal
    # params.stopwhen = alleq

    report = saturate!(g, t, params)

    # display(g.classes); println()
    if !(report.reason === :saturated) && !reached(g, goal)
        return missing # failed to prove
    end
    return reached(g, goal)
end

macro areequal(theory, exprs...)
    esc(:(areequal($theory, $exprs...)))
end

macro areequalg(G, theory, exprs...)
    esc(:(areequal($G, $theory, $exprs...)))
end


text:
using MPSKit,MPSKitModels,TensorKit,Test

#we pass this function to dmrg, it gets executed once per iteration
function finalize(iter,state,ham,envs)
    final_bonddim = 12;

    upperbound = max_Ds(state);
    shouldincrease = reduce((a,i) -> a && dim(virtualspace(state,i)) >= final_bonddim || upperbound[i+1] == dim(virtualspace(state,i)),1:length(state),init=true);

    if (shouldincrease)
        (state,envs) = changebonds(state, ham, OptimalExpand(),envs);
    end

    return (state,envs)
end

#defining the hamiltonian
th = nonsym_ising_ham(lambda = 4.0);
szt = TensorMap([1 0;0 -1],ℂ^2,ℂ^2)

# ------------------------------------
# |     Drmg                         |
# ------------------------------------
ts = FiniteMPS(fill(TensorMap(rand,ComplexF64,ℂ^1*ℂ^2,ℂ^1),10));
(ts,envs,_) = find_groundstate(ts,th,Dmrg(finalize = finalize));

szval_finite = sum(expectation_value(ts,szt))/length(ts)
@test szval_finite ≈ 0 atol=1e-12

# ------------------------------------
# |     Drmg2                        |
# ------------------------------------
ts = FiniteMPS(fill(TensorMap(rand,ComplexF64,ℂ^1*ℂ^2,ℂ^1),10));
(ts,envs,_) = find_groundstate(ts,th,Dmrg2(trscheme = truncdim(15)));

szval_finite = sum(expectation_value(ts,szt))/length(ts)
@test szval_finite ≈ 0 atol=1e-12

# ------------------------------------
# |     Vumps                        |
# ------------------------------------
ts = InfiniteMPS([ℂ^2],[ℂ^50]);
(ts,envs,_) = find_groundstate(ts,th,Vumps(maxiter=400));

szval_infinite = sum(expectation_value(ts,szt))/length(ts)
@test szval_infinite ≈ 0 atol=1e-12

# ------------------------------------
# |     Gradient optimization        |
# ------------------------------------
ts = InfiniteMPS([ℂ^2], [ℂ^5]);
(ts, envs, _) = find_groundstate(ts, th, GradientGrassmann(maxiter=400));

szval_infinite = sum(expectation_value(ts,szt))/length(ts)
@test szval_infinite ≈ 0 atol=1e-12


text:
# test KernelEmbeddingDistribution

using BenchmarkTools
using Distributions
using ApproxManifoldProducts
import ApproxManifoldProducts: mmd!
# using StaticArrays
# using KernelDensityEstimate

function main()
  P = randn(100)
  Q = randn(100)
  # v = StaticArray(0.0,)
  v = Float64[0.0]

  @btime mmd!(v, P, Q)

  nothing
end

main()

function main2()

  P1_ = kde!(randn(100))
  Q1 = kde!(randn(100))

  @btime kld( P1_, Q1)
  nothing
end

main2()

#


text:
using FermiCG
using Printf
using Test
using LinearAlgebra
using Profile
using HDF5
using Random
using StatProfilerHTML

using PyCall

molecule = "
C     1.259020     0.711710     0.000000
C     2.503660     1.386310     0.000000
C     3.708550     0.701730     0.000000
C     3.708550    -0.701730     0.000000
C     2.503660    -1.386310     0.000000
C     1.259020    -0.711710     0.000000
C    -0.013150    -1.446190     0.000000
C    -0.051250    -2.861390     0.000000
C    -1.246560    -3.562570     0.000000
C    -2.461990    -2.860840     0.000000
C    -2.452410    -1.475080     0.000000
C    -1.245860    -0.734490     0.000000
C    -1.245860     0.734490     0.000000
C    -0.013150     1.446190     0.000000
C    -0.051250     2.861390     0.000000
C    -1.246560     3.562570     0.000000
C    -2.461990     2.860840     0.000000
C    -2.452410     1.475080     0.000000
H     2.529850     2.475260     0.000000
H     4.650090     1.255640     0.000000
H     4.650090    -1.255640     0.000000
H     2.529850    -2.475260     0.000000
H     0.878710    -3.428550     0.000000
H    -1.237620    -4.654920     0.000000
H    -3.412460    -3.399270     0.000000
H    -3.408560    -0.953290     0.000000
H     0.878710     3.428550     0.000000
H    -1.237620     4.654920     0.000000
H    -3.412460     3.399270     0.000000
H    -3.408560     0.953290     0.000000
"

atoms = []
#molecule = lstrip(molecule)
for (li,line) in enumerate(split(rstrip(lstrip(molecule)), "\n"))
    l = split(line)
    push!(atoms, Atom(li, l[1], parse.(Float64,l[2:4])))
end

cas_nel = 18
cas_norb = 18

clusters    = [(1:6), (7:12), (13:16)]
init_fspace = [(3, 3),(3, 3),(3, 3)]

na = sum([i[1] for i in init_fspace])
nb = sum([i[2] for i in init_fspace])

basis = "sto-3g"
mol     = Molecule(0,1,atoms,basis)

tot_na = (18*6 + 12) ÷ 2
tot_nb = (18*6 + 12) ÷ 2
tot_n_elec = tot_na + tot_nb

#np = pyimport("numpy")
#C = np.load("test/c18.orbs.npy")
# get integrals
mf = FermiCG.pyscf_do_scf(mol, verbose=2)
C = mf.mo_coeff
nbas = size(C)[1]

FermiCG.pyscf_write_molden(mol, C, filename="hf.molden")

#
# Get active space
act_space = [46, 51, 52, collect(55:68)..., 72]
n_frozen = 51
#n_frozen = 4*18 + 12
inactive_space = setdiff(1:nbas,act_space)

occ_space = inactive_space[1:n_frozen]
vir_space = inactive_space[n_frozen+1:end]
Cact = C[:,act_space]
Cocc = C[:,occ_space]
Cvir = C[:,vir_space]

ints = FermiCG.pyscf_build_ints(mol, C[:,act_space], 2.0*Cocc*Cocc');

#pyscf = pyimport("pyscf")
#mcscf = pyimport("pyscf.mcscf")
#ao2mo = pyimport("pyscf.ao2mo")
#cas_norb = length(act_space)
#mycas = pyscf.mcscf.CASSCF(mf, cas_norb, 18)
#h1e_cas, ecore = mycas.get_h1eff(mo_coeff = Cocc)  #core core orbs to form ecore and eff
#h2e_cas = ao2mo.kernel(mol, Cact, aosym="s4",compact=false).reshape(4 * ((cas_norb), ))

Cl = FermiCG.localize(Cact, "boys", mf)
#
# use fiedler vector to reorder orbitals
h,j,k = FermiCG.pyscf_get_jk(mol, C[:,1:tot_na] * C[:,1:tot_na]')
Cl = FermiCG.fiedler_sort(Cl,k)

FermiCG.pyscf_write_molden(mol, Cl, filename="boys.molden")
FermiCG.pyscf_write_molden(mol, Cocc, filename="occ.molden")

n_core = 18
core_space = inactive_space[1:n_core]
sig_space = inactive_space[n_core+1:end]
Csig = FermiCG.localize(C[:,sig_space], "boys", mf)
Csig = FermiCG.fiedler_sort(Csig,(h+j-.5*k))
FermiCG.pyscf_write_molden(mol, Csig, filename="sig.molden")

S = FermiCG.get_ovlp(mf)
U =  Cact' * S * Cl
println(" Rotate Integrals")
flush(stdout)
ints = FermiCG.orbital_rotation(ints,U)
println(" done.")
flush(stdout)

max_stars_repo_name: nmayhall-vt/FermiCG

max_issues_repo_name: nmayhall-vt/FermiCG

max_forks_repo_name: nmayhall-vt/FermiCG


text:
# CPII-E 32
# peso especifico do cimento kg/m3
gama_c = 3100.0

# AREIA
# modulo de finura da areia
MF = 2.6
# peso especifico kg/m3
gama_a = 2650.0
# massa unitaria kg/m3
massa_a = 1470.0

# BRITA
# peso especifico kg/cm3
gama_b = 2700.0
# massa unitaria kg/cm3
massa_b = 1500.0
Dmax = 25 #mm

# CONCRETO
fck = 25.0 #MPa
abat = 90 #mm
abat_tolerancia = 10 #+/-mm
sd = 5.5

# PROPORÇÃO DAS BRITAS
B1 = 0.8
B2 = 0.2

# 1 - calcular resistencia do corpo de prova MPa
@show fc28 = fck + 1.645 * sd

# 2 - encontrar relação água-cimento na curva de Abrams
@show relacao_ac = 0.475

# 3 - calcular modulo de elasticidade GPa
@show Eci = (5600.0 * sqrt(fck)) / 1000.0

# 4 - Encontrar consumo de água na tabela l/m3
@show consumo_agua = 200

# 5 - consumo de cimento kg/m3
@show consumo_cimento = consumo_agua / relacao_ac

# 6 - consumo de agregado graúdo m3
@show consumo_agr_graudo = 0.715 * massa_b

# 6.1 - consumo de agregado graudo B1 e B2
@show consumo_B1 = B1 * consumo_agr_graudo
@show consumo_B2 = B2 * consumo_agr_graudo

# 7 - consumo de agregado miudo m3
@show volume_areia = 1.0 - ((consumo_cimento / gama_c) + (consumo_agr_graudo / gama_b) + (consumo_agua / 1000.0))
@show consumo_areia = volume_areia * gama_a

# 8 - apresentacao do traço
_c = consumo_cimento / consumo_cimento
_a = consumo_areia / consumo_cimento
_b1 = consumo_B1 / consumo_cimento
_b2 = consumo_B2 / consumo_cimento
_ac = consumo_agua / consumo_cimento

println("\n\n")
println("TRAÇO: $(round(_c, digits=3)):$(round(_a, digits=3)):$(round(_b1, digits=3)):$(round(_b2, digits=3)):$(round(_ac, digits=3))")
println("\n--- RESUMO ---\n")
println("Consumo cimento (kg): $(round(consumo_cimento ,digits=1)) --> $(round(consumo_cimento / 50.0,digits=1)) sacos de 50kg")
println("Consumo areia (kg): $(round(consumo_areia ,digits=1))")
println("Consumo brita 1 (kg): $(round(consumo_B1 ,digits=1))")
println("Consumo brita 2 (kg): $(round(consumo_B2 ,digits=1))")
println("Consumo água (l): $(round(consumo_agua ,digits=1))")
println("\n\n")

max_stars_repo_name: rsd-consultorias/scripts-engenharia

max_issues_repo_name: rsd-consultorias/scripts-engenharia

max_forks_repo_name: rsd-consultorias/scripts-engenharia


text:
using Pkg
Pkg.add("Distributions")
Pkg.add("Plots")
Pkg.add("PyPlot")