text:
### A Pluto.jl notebook ###
# v0.11.8

using Markdown
using InteractiveUtils

# ╔═╡ 25700c98-eebd-11ea-17eb-4398f75594e0
using Plots

# ╔═╡ 8666d19e-eebd-11ea-11ad-93097ff088f1
abstract type Tracker end

# ╔═╡ 8bb6e5e4-eebd-11ea-2019-49875c096a68

# ╔═╡ c637a44c-eebd-11ea-3d4f-9539b52a7b88

# ╔═╡ e4b6996e-eebf-11ea-3f8d-a13ddf340bbf
md"# Neighborhoods"

# ╔═╡ df0980b2-eec0-11ea-1775-e185922750ea
begin
	abstract type Neighborhood end

	struct OneFlip <: Neighborhood end

	oneflip = OneFlip()

	struct OneSwap <: Neighborhood end

	oneswap = OneSwap()
end

# ╔═╡ 47d523dc-eec1-11ea-2d55-e9c02f14cf4f
begin
	flip(s, i) = [j==i ? !sj : sj for (j, sj) in enumerate(s)]

	all_neighbors(s, ::OneFlip) = (flip(s, i) for i in 1:length(s))

	function swap(s, i, j)
		n = copy(s)
		n[i], n[j] = n[j], n[i]
		return n
	end

	all_neighbors(s, ::OneSwap) = (swap(s, i, j) for i in 1:length(s), j in 1:length(s) if s[i] != s[j])
end

# ╔═╡ 019334ee-eec2-11ea-2383-19f70f997887
all_neighbors([true, false, true, false], oneflip) |> collect

# ╔═╡ 9e624e5e-eec2-11ea-1aec-81cbeaee84c4
all_neighbors([true, false, true, false], oneswap) |> collect

# ╔═╡ 3f909330-eec3-11ea-1bd9-d7bb41cc052d
begin
	rand_neighbor(s, ::OneFlip) = flip(s, rand(1:length(s)))

	function rand_neighbor(s, ::OneSwap)
		i = rand(1:length(s))
		j = i
		while s[j] == s[i]
			j = rand(1:length(s))
		end
		return swap(s, i, j)
	end
end

# ╔═╡ c2abc1ea-eec3-11ea-3809-45bf92f5d133
rand_neighbor([true, false, true, false], oneflip)

# ╔═╡ caaf2bf2-eec3-11ea-1b01-9164c2d7d161
rand_neighbor([true, false, true, false], oneswap)

# ╔═╡ fcef48fe-eec3-11ea-22c2-0d17b1ee9af5
md"# Selectors"

# ╔═╡ 0492484a-eec4-11ea-091f-5146073f6824
abstract type Selector end

# ╔═╡ 1966a554-eec4-11ea-20d6-254b7898cf8d
begin

	struct BestNeighbor <: Selector end

	struct FirstNeighbor <: Selector end

	struct TwoStage <: Selector end

	struct RandomImprovement <: Selector end

	struct Metropolis <: Selector
		t
	end
end

# ╔═╡ b358ff22-eec4-11ea-2bbf-f12dbbf08c3a
begin
	function choose_next(f, s, N::Neighborhood, BestNeighbor)
		for i in

end

# ╔═╡ e6c4a97a-eebe-11ea-259c-d3433ddd579d
begin
	struct NoTracking <: Tracker end

	struct SolStoring{T} <: Tracker
		solutions::Vector{T}
	end

	notrack = NoTracking()
end

# ╔═╡ 25968390-eebe-11ea-2c8b-8d472bca08c8
begin
	# don't track anything
	track(::NoTracking, s) = nothing

	# store solutons
	track(tracker::SolStoring, s) = push!(tracker.solutions, s)
end

# ╔═╡ 98eb82b0-eebd-11ea-1403-a56e2c9baed5
begin
	choose_next(f, s, N::Neighborhood, S::Selector) = s

end

# ╔═╡ 3cd3fbba-eebd-11ea-164f-b39f495dea6e
function local_search(f, s₀, N::Neighborhood,
						S::Selector,
						tracker::Tracker=notrack;
						max_iter=1000)
	s = s₀
	track(f, s)
	for i in 1:max_iter
		s = choose_next(f, s, N, S)
		track(f, s)
	end
	return s
end

# ╔═╡ 5870cf14-eebe-11ea-07cd-d51bcd1fa702

# ╔═╡ Cell order:
# ╠═25700c98-eebd-11ea-17eb-4398f75594e0
# ╠═8666d19e-eebd-11ea-11ad-93097ff088f1
# ╠═8bb6e5e4-eebd-11ea-2019-49875c096a68
# ╠═c637a44c-eebd-11ea-3d4f-9539b52a7b88
# ╠═3cd3fbba-eebd-11ea-164f-b39f495dea6e
# ╠═e4b6996e-eebf-11ea-3f8d-a13ddf340bbf
# ╠═df0980b2-eec0-11ea-1775-e185922750ea
# ╠═47d523dc-eec1-11ea-2d55-e9c02f14cf4f
# ╠═019334ee-eec2-11ea-2383-19f70f997887
# ╠═9e624e5e-eec2-11ea-1aec-81cbeaee84c4
# ╠═3f909330-eec3-11ea-1bd9-d7bb41cc052d
# ╠═c2abc1ea-eec3-11ea-3809-45bf92f5d133
# ╠═caaf2bf2-eec3-11ea-1b01-9164c2d7d161
# ╠═fcef48fe-eec3-11ea-22c2-0d17b1ee9af5
# ╠═0492484a-eec4-11ea-091f-5146073f6824
# ╠═1966a554-eec4-11ea-20d6-254b7898cf8d
# ╠═b358ff22-eec4-11ea-2bbf-f12dbbf08c3a
# ╠═25968390-eebe-11ea-2c8b-8d472bca08c8
# ╠═e6c4a97a-eebe-11ea-259c-d3433ddd579d
# ╠═98eb82b0-eebd-11ea-1403-a56e2c9baed5
# ╠═5870cf14-eebe-11ea-07cd-d51bcd1fa702

max_stars_repo_name: mathematiguy/STMO

max_issues_repo_name: mathematiguy/STMO

max_forks_repo_name: mathematiguy/STMO


text:
using NLsolve

# Precondition: ∃ J ∈ ℤ: length(f) = 2^J
# Returns: ((a0, b0, b1, ...), length(a0))
function mallatdecomp(f, L, H, p)
    out = []
    a = f
    while length(a) >= 2p
        out = [H(a); out]
        a = L(a)
    end
    out = [a; out]
    return (out, length(a))
end

# Precondition: ∃ J ∈ ℤ: length(f) = 2^J,
#   length(a0) = length(b0),
#   length(b_{i+1}) = 2 length(b_i)
# Input: Wf = (a0, b0, b1, ...), n = length(a0)
function mallatrecon(Wf, n, iL, iH)
    a = Wf[1:n]
    i = n
    while i + n <= length(Wf)
        b = Wf[(i+1):(i+n)]
        a = iL(a) + iH(b)
        i += n
        n = 2n
    end
    return a
end

# Returns: out[i+1] = c_i for i ∈ {0, 2p-1} for DAUB(2p) wavelets
function daubgen(p)
    function f!(F,c)
        for m = 1:(p-1)
            k = collect((2m):(2p-1))
            F[m] = sum(c[k .+ 1] .* c[k .- 2m .+ 1])
        end
        F[p] = sum(c.^2) - 2
        k = collect(0:(2p - 1))
        for m = 1:p
            F[p+m] = sum((-1).^k .* k.^(m-1) .* c)
        end
    end

    function j!(J, c)
        for m = 1:(p-1)
            J[m, (2p-2m+1):(2p)] = zeros(1, 2m)
            for i = 0:(2p - 2m - 1)
                J[m,i+1] = c[i+1+2m]
            end
            for i = (2m):(2p-1)
                J[m,i+1] += c[i+1-2m]
            end
        end
        J[p,:] = 2 * c
        k = collect(0:(2p - 1))
        for m = 1:p
            J[p+m, :] = (-1).^k .* k.^(m-1)
        end
    end

    r = nlsolve(f!, j!, ones(2p))
    if converged(r)
        return r.zero
    end
    throw(Exception)
end

function filtergen(p)
    c = daubgen(p)

    function L(f)
        n = length(f)
        out = zeros(n >> 1)
        for i = 1:(n >> 1)
            s = 0
            for j = (2(i-p) + 1):(2i)
                s += f[(n+j-1)%n+1] * c[(2i - j + 1)]
            end
            out[i] = s / 2
        end
        return out
    end

    function iL(g)
        nhalf = length(g)
        n = 2 * nhalf
        out = zeros(n)
        for j = 1:n
            s = 0
            for i = ((j+1)>>1):(p-1 + (j+1)>>1)
                s += g[(i-1)%nhalf+1] * c[(2i - j + 1)]
            end
            out[j] = s
        end
        return out
    end

    function H(f)
        n = length(f)
        out = zeros(n >> 1)
        for i = 1:(n >> 1)
            s = 0
            for j = (2i-1):(2(p+i-1))
                s += f[(j-1)%n+1] * (-1)^(j+1) * c[(j -2i + 2)]
            end
            out[i] = s / 2
        end
        return out
    end

    function iH(g)
        nhalf = length(g)
        n = 2 * nhalf
        out = zeros(n)
        for j = 1:n
            s = 0
            for i = ((j+1)>>1 + 1 - p):((j+1)>>1)
                s += g[(nhalf+i-1)%nhalf+1] * (-1)^(j+1) * c[(j - 2i + 2)]
            end
            out[j] = s
        end
        return out
    end

    return (L, H, iL, iH)
end

ext: jl


text:
abstract type Initializers end

struct Uniform <: Initializers
    min::AbstractFloat
    max::AbstractFloat
    ortho::Bool
    function Uniform(;min=-1.0,max=1.0,ortho=false)
        @assert(min<max)
        new(min,max,ortho)
    end
end

function orthonormalize(X::Union{AbstractArray{T},AbstractArray{Complex{T}}}) where {T<:AbstractFloat}
    U,S,V = svd(X)
    return U*transpose(V)
end

function orthonormalize(X::AbstractArray{Quaternion{T}}) where {T<:AbstractFloat}
    Xc = Quaternions.equiv(X)
    U,S,V = svd(Xc)
    return Quaternions.equiv(U*transpose(V))
end

function rand_arr(T::Type{S}, dim::Tuple, init::Uniform) where {S<:AbstractFloat}
    M = rand(T, dim)*S(init.max-init.min) .+ T(init.min)
    return init.ortho ? orthonormalize(M) : M
end

function rand_arr(T::Type{Complex{S}}, dim::Tuple, init::Uniform) where {S<:AbstractFloat}
    M = rand(T, dim)*S(init.max-init.min) .+ T(init.min, init.min)
    return init.ortho ? orthonormalize(M) : M
end

function rand_arr(T::Type{Quaternion{S}}, dim::Tuple, init::Uniform) where {S<:AbstractFloat}
    M = rand(T, dim)*S(init.max-init.min) .+ T(init.min, init.min, init.min, init.min)
    return init.ortho ? orthonormalize(M) : M
end


text:
using Test
using QuantLib
using Dates

# build a fixed rate bond to test cash flows
settlement_date = Date(2008, 9, 18)
set_eval_date!(settings, settlement_date - Dates.Day(3))
bbCurveRate = 0.055
bbDayCount = QuantLib.Time.ActualActualBond()
bbIR = InterestRate(bbCurveRate, bbDayCount, CompoundedCompounding(), QuantLib.Time.Semiannual())
ts = FlatForwardTermStructure(settlement_date, bbIR.rate, bbIR.dc, bbIR.comp, bbIR.freq)
settlement_days = 3
face_amount = 100.0
fixed_schedule = QuantLib.Time.Schedule(Date(2007, 5, 15), Date(2017, 5, 15),
                       QuantLib.Time.TenorPeriod(QuantLib.Time.Semiannual()), QuantLib.Time.Unadjusted(),
                       QuantLib.Time.Unadjusted(), QuantLib.Time.DateGenerationBackwards(), false,
                       QuantLib.Time.USGovernmentBondCalendar())

pe = DiscountingBondEngine(ts)
fixedrate_bond = FixedRateBond(settlement_days, face_amount, fixed_schedule, 0.045,
                        QuantLib.Time.ISMAActualActual(), QuantLib.Time.ModifiedFollowing(), 100.0,
                        Date(2007, 5, 15), fixed_schedule.cal, pe)

# generate a floating rate bond
floating_cal = QuantLib.Time.TargetCalendar()
settlement_days = 3
cap_vol = ConstantOptionVolatility(3, floating_cal, QuantLib.Time.ModifiedFollowing(), 0.0, QuantLib.Time.Actual365())

fb_issue_date = Date(2005, 10, 21)
fb_bond_engine = DiscountingBondEngine(ts)
fb_dc = QuantLib.Time.Actual360()
conv = QuantLib.Time.ModifiedFollowing()

fb_schedule = QuantLib.Time.Schedule(Date(2005, 10, 21), Date(2010, 10, 21), QuantLib.Time.TenorPeriod(QuantLib.Time.Quarterly()),
                                        QuantLib.Time.Unadjusted(), QuantLib.Time.Unadjusted(), QuantLib.Time.DateGenerationBackwards(), false,
                                        QuantLib.Time.USNYSECalendar())
fixing_days = 2
in_arrears = true
gearings = ones(length(fb_schedule.dates) - 1)
spreads = fill(0.001, length(fb_schedule.dates) - 1)
libor_3m = usd_libor_index(QuantLib.Time.TenorPeriod(Dates.Month(3)), ts)
floating_bond = FloatingRateBond(settlement_days, face_amount, fb_schedule, libor_3m, fb_dc, conv, fixing_days, fb_issue_date, fb_bond_engine, in_arrears, 100.0, gearings, spreads, cap_vol=cap_vol)

# test
@test accrual_start_date(fixedrate_bond.cashflows.coupons[1]) == Date(2007, 5, 15)
@test accrual_end_date(fixedrate_bond.cashflows.coupons[1]) == Date(2007, 11, 15)
@test ref_period_start(fixedrate_bond.cashflows.coupons[1]) == Date(2007, 5, 15)
@test ref_period_end(fixedrate_bond.cashflows.coupons[1]) == Date(2007, 11, 15)
@test QuantLib.accrual_period(fixedrate_bond.cashflows.coupons[1]) == 0.5
@test QuantLib.get_latest_coupon(fixedrate_bond.cashflows) == fixedrate_bond.cashflows.coupons[end]

# test NPV calculations
@test npv(fixedrate_bond.cashflows, ts, ts.referenceDate, ts.referenceDate) == 94.67245111530346
@test npv(fixedrate_bond.cashflows, bbIR, true, ts.referenceDate, ts.referenceDate) == 94.71618240316593
@test QuantLib.npvbps(fixedrate_bond.cashflows, ts, ts.referenceDate, ts.referenceDate) == (94.67245111530346,0.07152580522996536)

# duration and yield
@test duration(ModifiedDuration(), fixedrate_bond.cashflows, bbIR, bbIR.dc, true, ts.referenceDate, ts.referenceDate) == 6.899980120600781
@test QuantLib.yield(fixedrate_bond.cashflows, 94.67245111530346, bbIR.dc, bbIR.comp, bbIR.freq, true, ts.referenceDate, ts.referenceDate, 1.0e-10, 100, 0.05) == 0.05506693320721388

# Other calculations
@test accrued_amount(fixedrate_bond.cashflows, ts.referenceDate) == 1.5407608695652275
@test previous_cashflow_date(fixedrate_bond.cashflows, ts.referenceDate) == Date(2008, 5, 15)
@test accrual_days(fixedrate_bond.cashflows, bbIR.dc, ts.referenceDate) == 126
@test next_cashflow(fixedrate_bond.cashflows, ts.referenceDate) == 3
@test next_coupon_rate(fixedrate_bond.cashflows, ts.referenceDate) == 0.045
@test has_occurred(fixedrate_bond.cashflows.coupons[end-1], ts.referenceDate) == false
@test maturity_date(fixedrate_bond.cashflows) == Date(2017, 5, 15)
@test amount(fixedrate_bond.cashflows.coupons[1]) == 2.2499999999999964
@test QuantLib.calc_rate(fixedrate_bond.cashflows.coupons[1]) == 0.045

# Floating coupon calculations
@test QuantLib.calc_rate(floating_bond.cashflows[12]) == 0.05443944335694899
@test amount(floating_bond.cashflows[12]) == 1.3912302191220296
@test QuantLib.index_fixing(floating_bond.cashflows[12]) == 0.05343944335694899
@test accrued_amount(floating_bond.cashflows[12], ts.referenceDate) == 0.8922019883499973
@test accrued_amount(floating_bond.cashflows[1], ts.referenceDate) == 0.0


text:
struct HyperMeshData
    α::Vector{Float64}
    β::Vector{Float64}
    γ::Vector{Float64}
    λ::Vector{Float64}
    λi::Vector{Float64}
    Area0::Vector{Float64}
    Area1::Vector{Float64}
    rξ::Vector{MVector{2,Float64}}
    rξr::Vector{MVector{2,Float64}}
    rη::Vector{MVector{2,Float64}}
    S::Vector{MVector{2,Float64}}
    A::Vector{MMatrix{2,2,Float64,4}}
    B::Vector{MMatrix{2,2,Float64,4}}
    C::Vector{MMatrix{2,2,Float64,4}}
    function HyperMeshData(NI::Int)
        rξ = [MVector{2,Float64}(undef) for _ in 1:NI]
        rξr = [MVector{2,Float64}(undef) for _ in 1:NI]
        rη = [MVector{2,Float64}(undef) for _ in 1:NI]
        α = Vector{Float64}(undef, NI)
        β = Vector{Float64}(undef, NI)
        γ = Vector{Float64}(undef, NI)
        λ = Vector{Float64}(undef, NI)
        λi = Vector{Float64}(undef, NI)
        Area0 = Vector{Float64}(undef, NI)
        Area1 = Vector{Float64}(undef, NI)
        A = [MMatrix{2,2,Float64}(undef) for _ in 1:NI]
        B = [MMatrix{2,2,Float64}(undef) for _ in 1:NI]
        C = [MMatrix{2,2,Float64}(undef) for _ in 1:NI]
        S = [MVector{2,Float64}(undef) for _ in 1:NI]
        return new(α, β, γ, λ, λi, Area0, Area1, rξ, rξr, rη, S, A, B, C)
    end
end

function gen_mesh(
    mshr::HyperbolicMesher,
    p0::AbstractVector{MVector{2,Float64}},
    height::Real;
    mg_level::Int=0,
    timer::Bool=true,
)::AirfoilMesh
    @timeit to "gen_mesh" begin
        @assert height > 0
        @assert mg_level >= 0
        check_points(mshr.topo, p0)
        np_wake = 0
        if mshr.topo == :C && mshr.wake.len > 0
            wake = create_wake(mshr.wake, p0; mg_level=mg_level)
            np_wake = length(wake) + 1
            p0 = join_wake(p0, wake)
        end
        steps = gen_steps(mshr.step, convert(Float64, height), mg_level)
        @assert length(steps) % 2^mg_level == 0
        NJ = length(steps) + 1
        NI = length(p0)
        grid = Matrix{MVector{2,Float64}}(undef, NI, NJ)
        data = HyperMeshData(NI)
        grid[:, 1] .= p0

        if mshr.topo == :O # BCs are ignored in :O mesh
            bc0 = mshr.bc0
            bc1 = mshr.bc1
        else
            bc0 = if isa(mshr.bc0, SymmetryBoundary)
                dx, dy = p0[2] - p0[1]
                LineBoundary(-dy, dx)
            else
                mshr.bc0
            end
            bc1 = if isa(mshr.bc1, SymmetryBoundary)
                dx, dy = p0[end] - p0[end - 1]
                LineBoundary(-dy, dx)
            else
                mshr.bc1
            end
        end
        for j in 2:NJ
            for i in 1:NI
                grid[i, j] = MVector{2,Float64}(undef)
            end
            h0 = mshr.step(j - 2)
            h1 = mshr.step(j - 1)
            @timeit to "extrusion" hyp_run!(
                mshr,
                data,
                mshr.f_dissi,
                (1 + (1 - 0.4mshr.f_blend)^(j - 2)) / 2,
                h0,
                h1,
                bc0,
                bc1,
                view(grid, :, j - 1),
                view(grid, :, j),
            )
            if j == 2
                for i in 1:NI
                    rη = normalize(data.rη[i])
                    dη = grid[i, j] - grid[i, j - 1]
                    hc = norm(dη)
                    dc = dot(dη, rη)
                    if dc < 0 # negative direction
                        grid[i, j] .= rη * 0.5 * h1 + grid[i, j - 1]
                    elseif hc < h1
                        grid[i, j] .= (h1 + hc) / (2hc) * dη + grid[i, j - 1]
                    end
                end
            elseif (mshr.n_smooth > 0 && mshr.f_smooth > 0)
                α0 = mshr.f_smooth / 2
                α = min(α0, α0 * (j - 2) / NJ)
                if α > 0
                    β = 1 - α
                    for _ in 1:(mshr.n_smooth)
                        @timeit to "smooth sweept" for i in 2:(NI - 1)
                            dη = grid[i, j] - grid[i, j - 1]
                            hc = norm(dη)
                            if hc < 0.5h1
                                continue
                            end
                            p₋ = grid[i - 1, j]
                            p₊ = grid[i + 1, j]
                            p = grid[i, j]
                            lm = norm(p - p₋)
                            lp = norm(p₊ - p)
                            lmp = lm + lp
                            lm /= lmp
                            lp /= lmp
                            dss = lm * p₊ + lp * p₋
                            grid[i, j] .= β * p + α * dss
                        end
                    end
                end
            end
        end
    end
    if timer
        @info "$NI×$NJ Grid Extruded"
        show(to)
        println()
    end
    reset_timer!(to)
    nconn = if mshr.topo == :C
        np_wake
    elseif mshr.topo == :O
        NJ
    else
        0
    end
    return AirfoilMesh(mshr.topo, grid, nconn)
end

function hyp_run!(
    mshr::HyperbolicMesher,
    dat::HyperMeshData,
    fdis::Float64,
    fbld::Float64,
    h0::Float64,
    h1::Float64,
    bc0::BoundaryCondition,
    bc1::BoundaryCondition,
    r0::L0,
    r1::L1,
) where {L0<:AbstractVector{MVector{2,Float64}},L1<:AbstractVector{MVector{2,Float64}}}
    nξ = length(r0)
    # rξ
    @inbounds for ξ in 2:(nξ - 1)
        dat.rξ[ξ] .= (r0[ξ + 1] - r0[ξ - 1]) / 2
    end
    if mshr.topo == :O
        dat.rξ[1] .= (r0[2] - r0[nξ - 1]) / 2
        dat.rξ[nξ] .= dat.rξ[1]
    else
        dat.rξ[1] .= r0[2] - r0[1]
        dat.rξ[nξ] .= r0[nξ] - r0[nξ - 1]
    end
    # rξr, γ, Area0, Area1
    for ξ in 1:nξ
        xξ, yξ = dat.rξ[ξ]
        dat.rξr[ξ] .= (-yξ, xξ)
        dat.γ[ξ] = γ = xξ^2 + yξ^2
        γ_sqr = sqrt(γ)
        dat.Area0[ξ] = A0 = γ_sqr * h0
        dat.Area1[ξ] = γ_sqr * h1
        dat.rη[ξ] .= dat.rξr[ξ] * (A0 / γ)
    end
    # extrude ends
    if mshr.topo == :O
        r1[end] = r1[1] = r0[1] + dat.rξr[1] / sqrt(dat.γ[1]) * h1
    else
        r1[1] = r0[1] + one_step(bc0, h1, dat.rξr[1])
        r1[nξ] = r0[nξ] + one_step(bc1, h1, dat.rξr[nξ])
    end
    # vol_blending
    if fbld != 1
        AM1 = mean(dat.Area1)
        @. dat.Area1 = fbld * dat.Area1 + (1 - fbld) * AM1
    end
    # α, β, λ
    for ξ in 1:nξ
        dxξ0, dyξ0 = dat.rξ[ξ]
        dxη0, dyη0 = dat.rη[ξ]
        α = dxξ0 .* dxη0 - dyξ0 .* dyη0
        β = dxξ0 .* dyη0 + dxη0 .* dyξ0
        γ = dat.γ[ξ]
        dat.λ[ξ] = sqrt(α^2 + β^2) / γ * fdis
        dat.C[ξ][1, 1] = α
        dat.C[ξ][2, 1] = β
        dat.C[ξ][1, 2] = β
        dat.C[ξ][2, 2] = -α
        dat.C[ξ] ./= 2γ
        dat.A[ξ] .= -dat.C[ξ]
        dat.B[ξ] .= I2
        AS = (dat.Area0[ξ] + dat.Area1[ξ]) / γ
        dat.S[ξ] .= dat.rξr[ξ] * AS + r0[ξ]
    end
    # update dissipation at interface
    if fdis > 0
        for ξ in 1:(nξ - 1)
            dat.λi[ξ] = sqrt(dat.λ[ξ] * dat.λ[ξ + 1])
        end
        # apply dissipation
        for ξ in 2:(nξ - 1)
            λ₋ = 0.5(dat.λi[ξ - 1])
            λ₊ = 0.5(dat.λi[ξ])
            diag_sub(dat.A[ξ], λ₋)
            diag_sub(dat.C[ξ], λ₊)
            diag_add(dat.B[ξ], λ₋ + λ₊)
        end
    end
    # solve
    dat.S[2] -= dat.A[2] * r1[1]
    dat.S[nξ - 1] -= dat.C[nξ - 1] * r1[nξ]
    solve_block_tridiag!(
        view(dat.A, 2:(nξ - 1)),
        view(dat.B, 2:(nξ - 1)),
        view(dat.C, 2:(nξ - 1)),
        view(dat.S, 2:(nξ - 1)),
    )
    for ξ in 2:(nξ - 1)
        r1[ξ] .= dat.S[ξ]
    end
    return nothing
end


text:
# This file was generated, do not modify it. # hide
MLJ.save("tuned_iterated_pipe.jls", mach_tuned_iterated_pipe)

max_stars_repo_name: alan-turing-institute/MLJTutorials

max_issues_repo_name: alan-turing-institute/MLJTutorials

max_forks_repo_name: alan-turing-institute/MLJTutorials


text:
using FFTW

export Poisson

"""
    Poisson( grid )

Derived type to solve the Poisson equation on 2d regular
cartesian mesh with periodic boundary conditions on both sides
- kx   : wave number in x
- ky   : wave number in y
- k2   : ``k_x^2 + k_y^2``
- nc_x : cells number in x
- nc_y : cells number in y
- dx   : x step size
- dy   : y step size
- rht  : fft(rho)
"""
struct Poisson

    grid::Grid
    kx::Array{Float64,2}
    ky::Array{Float64,2}
    k2::Array{Float64,2}
    rht::Array{ComplexF64,2}

    function Poisson(grid::Grid)

        nx = grid.nx
        ny = grid.ny

        rht = zeros(ComplexF64, (div(nx, 2) + 1, ny))

        kx = zeros(nx ÷ 2 + 1, ny)
        ky = zeros(nx ÷ 2 + 1, ny)
        k2 = zeros(nx ÷ 2 + 1, ny)

        kx0 = 2π / grid.lx
        ky0 = 2π / grid.ly

        for ik = 1:nx÷2+1
            kx1 = (ik - 1) * kx0
            for jk = 1:ny÷2
                kx[ik, jk] = kx1
                ky[ik, jk] = (jk - 1) * ky0
            end
            for jk = ny÷2+1:ny
                kx[ik, jk] = kx1
                ky[ik, jk] = (jk - 1 - ny) * ky0
            end
        end

        kx[1, 1] = 1.0
        k2 .= kx .* kx .+ ky .* ky
        kx .= kx ./ k2
        ky .= ky ./ k2

        new(grid, kx, ky, k2, rht)

    end

end

export compute_psi!

"""
    compute_psi!( poisson, psi, omega )

computes ``\\psi`` from ``\\omega``

```math
-\\Delta \\psi(x,y) = \\omega(x,y)
```

"""
function compute_psi!(psi, poisson::Poisson, omega)

    poisson.rht .= rfft(omega)
    poisson.rht ./= poisson.k2
    psi .= irfft(poisson.rht, poisson.grid.nx)

end

export compute_velocities!

"""
    compute_velocities!( ux, uy, poisson, omega )

solves Poisson equation to compute velocity fields

```math
\\Delta \\psi(x,y) = - \\omega(x,y) \\\\
u(x,y) = \\nabla \\times \\psi(x,y)
```

"""
function compute_velocities!(ux, uy, poisson::Poisson, omega)

    poisson.rht .= rfft(omega)

    poisson.rht[1, 1] = 0.0
    poisson.rht .*= -1im .* poisson.ky

    ux .= irfft(poisson.rht, poisson.grid.nx)

    poisson.rht .= rfft(omega)

    poisson.rht[1, 1] = 0.0
    poisson.rht .*= 1im .* poisson.kx

    uy .= irfft(poisson.rht, poisson.grid.nx)

end


text:
# Convert array of particular type to array of Dual
function seed_duals(x::AbstractArray{V},::Type{T},
                    N::Integer) where {V,T}
  seeds = ForwardDiff.construct_seeds(ForwardDiff.Partials{N,V})
  duals = [ForwardDiff.Dual{T}(x[i],seeds[i]) for i in eachindex(x)]
end

# Relu function is used as a constrain for γ's values
relu(x::Real) = max(zero(x), x)

# update method
function update!(opt, ps::Array, ps̄::Array, N)
  len = length(ps)
  for i in 1:len
    p_value = [element.value for element in ps[i]]
    p_value .-= apply!(opt, p_value, ps̄[i])
    new_p = seed_duals(p_value, eltype(ps̄[i]), N)
    ps[i] .= new_p
  end
end

# partition the gradient into sub-arrays for each parameter
function grads(θ̄::Vector{T}, indexes) where T
  ll = length(indexes)
  partition_θ̄ = Array{Array{T, 1}}(undef, ll)
  for i in 1:ll
    partition_θ̄[i] = θ̄[indexes[i]]
  end
  return partition_θ̄
end

max_stars_repo_name: HamletWantToCode/DiffGPFlux

max_issues_repo_name: HamletWantToCode/DiffGPFlux


text:
### A Pluto.jl notebook ###
# v0.12.10

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ 8491344a-7c39-11eb-327e-df48b337d147
begin
	import Pkg

	Pkg.add("Plots")
	Pkg.add("PlutoUI")
	Pkg.add("Measures")

	using Plots
	using PlutoUI
	using Measures
end

# ╔═╡ 0904e59a-7c3b-11eb-307d-6fd2049ea8bc
md"To use Eulers Method on an ODE with parameters you can either define global variables and use them within the function or to increase runtime pass all the parameters in a variable P to the vector field."

# ╔═╡ 8f2743fe-7c39-11eb-27d6-d527b21f46b8
begin
	function euler(f,t₀,tₑ,Δt,x₀,P)
		T = t₀:Δt:tₑ
		F = Vector{typeof(x₀)}(undef,length(T))

		F[1] = x₀

		for (n,tₙ) ∈ enumerate(@view T[2:end])
			F[n+1] = F[n] .+ f(tₙ,F[n],P) .* Δt
		end
		return F
	end

	#reshape result for plotting
	reshape_result(F) = [[f[i] for f ∈ F] for i ∈ 1:length(F[1])]
end;

# ╔═╡ 925558e8-7cc7-11eb-37a8-c98fb8a33355
md"""
---
### One Locus, Two Alleles

We describe a population where every individual is characterized by one locus only. There are two possible alleles 0 and 1, where 0 is the wildtype and 1 the mutated allele. Individuals who carry the mutated allel on both copies of the allele are thought to have a severe diseas, causing them to be excluded from the birth process. Individuals give birth and die at constant rates $b, d$. The carrying capacity $K_{t}$ is time dependent and jumps at a certain point in time. The competition preassure between individuals is uniform $(b-d)/K_{t}$. Mutations appear at every birth at a constant rate $\mu$ from 0 to 1 only. The deterministic system can be described by the following system of non linear first order ODEs:

${
\frac{dX_{00}}{dt} = b
	\left(
		(1-\mu) \frac{(X_{00}+\frac{1}{2}X_{01})^2}{X_{00}+X_{01}}
	\right)
	- (d + c\Sigma)X_{00}
		\phantom{
			+ \mu  \frac{(X_{00}+\frac{1}{2}X_{01})X_{01}}{X_{00}+X_{01}}
			+ \mu^2  \frac{(X_{00}+\frac{1}{2}X_{01})^2}{X_{00}+X_{01}}
			}
	\phantom{\Bigg(}

}$

${
\frac{dX_{01}}{dt} = b
	\left(
		(1-\mu) \frac{(X_{00}+\frac{1}{2}X_{01})X_{01}}{X_{00}+X_{01}}
		+ \mu  \frac{(X_{00}+\frac{1}{2}X_{01})^2}{X_{00}+X_{01}}
	\right)
	- (d + c\Sigma)X_{01}
	\phantom{
			+ \mu^2  \frac{(X_{00}+\frac{1}{2}X_{01})^2}{X_{00}+X_{01}}
		}
}$

${
\frac{dX_{11}}{dt} = b
	\left(
		(1-\mu) \frac{\frac{1}{4}X_{01}^2}{X_{00}+X_{01}}
		+ \mu  \frac{(X_{00}+\frac{1}{2}X_{01})X_{01}}{X_{00}+X_{01}}
		+ \mu^2  \frac{(X_{00}+\frac{1}{2}X_{01})^2}{X_{00}+X_{01}}
	\right)
	- (d + c\Sigma)X_{11}
}$

where $\Sigma = X_{00} + X_{01} + X_{11}$ is the total population size.
"""

# ╔═╡ 62311b66-7c3b-11eb-335c-c53b566fcb9f
begin
	function OneLocus(t,X,P)
		X₀₀, X₀₁, X₁₁ = (x for x ∈ X)
		b, d, K, μ = (p for p ∈ P)

		Σ = sum(X)
		M = Σ - X₁₁
		D = d+((b-d)/K(t))*Σ

		return [
				(b * (1-μ)*(X₀₀+0.5*X₀₁)^2/M - D*X₀₀),
				(b * ((1-μ)*(X₀₀+0.5*X₀₁)*X₀₁/M + μ*(X₀₀+0.5*X₀₁)^2/M) - D*X₀₁),
				(b * ((1-μ)*0.25*(X₀₁)^2/M + μ*(X₀₀+0.5*X₀₁)*X₀₁/M + μ^2*(X₀₀+0.5*X₀₁)^2/M) - D*X₁₁)
				]
		end
end

# ╔═╡ 348cb7ba-7c40-11eb-1a80-4768b0c6b4d4
md"""
birth rate: $(@bind b Slider(0.0:0.1:1.0, show_value=true, default=1.0)) |
death rate: $(@bind d Slider(0.0:0.1:1.0, show_value=true, default=0.9)) |

initial population: $(@bind Kstart NumberField(1:10^7, default=500.0)) |
carrying capacity: $(@bind Kend NumberField(1:10^7, default=10_000.0)) |
burn in time: $(@bind tburn Slider(0:300,show_value=true,default=100))

mutation rate: $(@bind μ NumberField(0.0:0.0001:0.2, default=0.001))
"""

# ╔═╡ a19430e8-7c3d-11eb-0790-4b883ac9cd39
begin
	t_start = 0
	t_end = 500
	step = 0.01

	OL_0 = [Kstart,0.0,0.0]

	T = t_start:step:t_end
	K(t) = t <= tburn ? Kstart : Kend ;

	F = euler(OneLocus,t_start,t_end,step,OL_0,[b,d,K,μ])
	rF = reshape_result(F)
end;

# ╔═╡ d728c114-7c40-11eb-1fcb-a1ecbb31ae54
let
	p = plot(framestyle=:zerolines)

	plot!(T,rF,label=["(0,0)" "(0,1)" "(1,1)"])

	p
end

# ╔═╡ 20a4f4dc-7c41-11eb-26a8-21e34fcfd8e4
begin
	mutation_load(X) = (X[2] + 2*X[3])/sum(X)
	ill_individual(X) = X[3]/sum(X)
end;

# ╔═╡ 619eab70-7c41-11eb-1fc0-fb3ec4e64220
let
	p_left = plot(framestyle=:zerolines,legend=:bottomleft,rightmargin=12mm,size=(680,400))

	plot!(p_left,T,mutation_load.(F),label="mutation load",color=:red)

	p_right = twinx()
	plot!(p_right,framestyle=:zerolines,legend=:bottomright,grid=false)
	plot!(p_right,T,ill_individual.(F),label="ill_individual",color=:orange)

	vline!(p_right,[T[end]],color=:black,label="")

end

# ╔═╡ Cell order:
# ╟─8491344a-7c39-11eb-327e-df48b337d147
# ╟─0904e59a-7c3b-11eb-307d-6fd2049ea8bc
# ╠═8f2743fe-7c39-11eb-27d6-d527b21f46b8
# ╟─925558e8-7cc7-11eb-37a8-c98fb8a33355
# ╠═62311b66-7c3b-11eb-335c-c53b566fcb9f
# ╟─348cb7ba-7c40-11eb-1a80-4768b0c6b4d4
# ╟─d728c114-7c40-11eb-1fcb-a1ecbb31ae54
# ╟─619eab70-7c41-11eb-1fc0-fb3ec4e64220
# ╠═a19430e8-7c3d-11eb-0790-4b883ac9cd39
# ╟─20a4f4dc-7c41-11eb-26a8-21e34fcfd8e4

max_stars_repo_name: roccminton/PlutoNotebooks

max_issues_repo_name: roccminton/PlutoNotebooks

max_forks_repo_name: roccminton/PlutoNotebooks


text:
struct GageCard
    handle::CSHANDLE
    function GageCard()
        CsInitialize()
        hndl = CsGetSystem()
        hndl == 0 ? error("System not found or already in use.") : new(hndl)
    end
end

function free(sys::GageCard)
    CsFreeSystem(sys.handle)
end

function getconfig_acq(sys::GageCard)
    aquisitionconfig = CsGetAcq(sys.handle, CS_CURRENT_CONFIGURATION)
    return aquisitionconfig
end

function getconfig_trig(sys::GageCard)
    triggercfg = CsGetTrig(sys.handle, CS_CURRENT_CONFIGURATION)
    return triggercfg
end

function getconfig_chan(sys::GageCard)
    channelcfg = CsGetChan(sys.handle, CS_CURRENT_CONFIGURATION)
    return channelcfg
end

sys = GageCard()

acqconfig = getconfig_chan(sys)
trigger = getconfig_trig(sys)
channel = getconfig_chan(sys)

free(sys)

println(trigger.i32Source)


text:
module TestConstrainingLens

include("preamble.jl")

@testset begin
    obj = (a = 1, b = 1)
    constraint = constraining() do obj
        @set obj.b = obj.a
    end
    lens = constraint ∘ @lens _.a

    @test get(obj, lens) === 1
    @test set(obj, lens, 2) === (a = 2, b = 2)
end

end  # module


text:
# !! auto-generated code, do not edit !!
module visualization_msgs
using RobotOSData.Messages
using RobotOSData.StdMsgs
using RobotOSData.CommonMsgs
struct Marker <: Readable
    header::Header
    ns::String
    id::Int32
    type_::Int32
    action::Int32
    pose::geometry_msgs.Pose
    scale::geometry_msgs.Vector3
    color::std_msgs.ColorRGBA
    lifetime::ROSDuration
    frame_locked::Bool
    points::Vector{geometry_msgs.Point}
    colors::Vector{std_msgs.ColorRGBA}
    text::String
    mesh_resource::String
    mesh_use_embedded_materials::Bool
end
struct InteractiveMarkerControl <: Readable
    name::String
    orientation::geometry_msgs.Quaternion
    orientation_mode::UInt8
    interaction_mode::UInt8
    always_visible::Bool
    markers::Vector{Marker}
    independent_marker_orientation::Bool
    description::String
end
struct MenuEntry <: Readable
    id::UInt32
    parent_id::UInt32
    title::String
    command::String
    command_type::UInt8
end
struct InteractiveMarker <: Readable
    header::Header
    pose::geometry_msgs.Pose
    name::String
    description::String
    scale::Float32
    menu_entries::Vector{MenuEntry}
    controls::Vector{InteractiveMarkerControl}
end
struct InteractiveMarkerInit <: Readable
    server_id::String
    seq_num::UInt64
    markers::Vector{InteractiveMarker}
end
struct InteractiveMarkerPose <: Readable
    header::Header
    pose::geometry_msgs.Pose
    name::String
end
struct InteractiveMarkerUpdate <: Readable
    server_id::String
    seq_num::UInt64
    type_::UInt8
    markers::Vector{InteractiveMarker}
    poses::Vector{InteractiveMarkerPose}
    erases::Vector{String}
end
struct MarkerArray <: Readable
    markers::Vector{Marker}
end
struct ImageMarker <: Readable
    header::Header
    ns::String
    id::Int32
    type_::Int32
    action::Int32
    position::geometry_msgs.Point
    scale::Float32
    outline_color::std_msgs.ColorRGBA
    filled::UInt8
    fill_color::std_msgs.ColorRGBA
    lifetime::ROSDuration
    points::Vector{geometry_msgs.Point}
    outline_colors::Vector{std_msgs.ColorRGBA}
end
struct InteractiveMarkerFeedback <: Readable
    header::Header
    client_id::String
    marker_name::String
    control_name::String
    event_type::UInt8
    pose::geometry_msgs.Pose
    menu_entry_id::UInt32
    mouse_point::geometry_msgs.Point
    mouse_point_valid::Bool
end
end


text:
using CachedArrays
using Test

#construct 2d and 3d cached slice series
#these should generalize to N-dimensions, but will just test with 4 and 5
data4 = rand(2,3,2,2)
data5 = rand(2,3,2,2,2)
cs24 = CachedArray(data4, 2)
cs25 = CachedArray(data5, 2)

@test parent(cs24) == data4
@test parent(cs25) == data5

@test cached_axes(cs24) == axes(cs24)[1:2]
@test cached_axes(cs25) == axes(cs25)[1:2]

@test noncached_axes(cs24) == axes(cs24)[3:end]
@test noncached_axes(cs25) == axes(cs25)[3:end]

#2d caching, 4d array
#now test various getindex implementations
#(probably overkill to test all of these because this package only has one custom getindex method)

@test all(cs24[1:2, 1:1, 2, 1] .== data4[1:2, 1:1, 2, 1])
@test all(cs24[10] .== data4[10])
@test all(cs24[CartesianIndex(10)] .== data4[CartesianIndex(10)])
@test all(cs24[1:2, 1:1, 1:2, :] .== data4[1:2, 1:1, 1:2, :])
@test all(cs24[1:2, 1:1, 1, :] .== data4[1:2, 1:1, 1, :])

cs24[1:2,1,1,1] .= -1.0
@test all(cs24[1:2,1,1,1] .== -1.0)
@test all(cs24[1:2,1,1,1] .== data4[1:2,1,1,1])

#2d caching, 5d array
@test all(cs25[1:2, 1:1, 2, 1, 2] .== data5[1:2, 1:1, 2, 1, 2])
@test all(cs25[10] .== data5[10])
@test all(cs25[CartesianIndex(10)] .== data5[CartesianIndex(10)])
@test all(cs25[1:2, 1:1, 1:2, :, :] .== data5[1:2, 1:1, 1:2, :, :])
@test all(cs25[1:2, 1:1, 1, :, :] .== data5[1:2, 1:1, 1, :, :])

cs25[1:2,1,1,1,2] .= -1.0
@test all(cs25[1:2,1,1,1,2] .== -1.0)
@test all(cs25[1:2,1,1,1,2] .== data5[1:2,1,1,1,2])

#@code_warntype cs25[1:2, 1:1, 2, 1, 2]

#3D caching, 4d array
data4 = rand(2,3,2,2)
cs34 = CachedArray(data4, 3)
@test all(cs34[1:2, 1:1, 2, 1] .== data4[1:2, 1:1, 2, 1])
@test all(cs34[10] .== data4[10])
@test all(cs34[CartesianIndex(10)] .== data4[CartesianIndex(10)])
@test all(cs34[1:2, 1:1, 1:2, :] .== data4[1:2, 1:1, 1:2, :])
@test all(cs34[1:2, 1:1, 1, :] .== data4[1:2, 1:1, 1, :])

cs34[1:2,1,1,1] .= -1.0
@test all(cs34[1:2,1,1,1] .== -1.0)
@test all(cs34[1:2,1,1,1] .== data4[1:2,1,1,1])

#3d caching, 5d array
data5 = rand(2,3,2,2,2)
cs35 = CachedArray(data5, 3)
@test all(cs35[1:2, 1:1, 2, 1, 2] .== data5[1:2, 1:1, 2, 1, 2])
@test all(cs35[10] .== data5[10])
@test all(cs35[CartesianIndex(10)] .== data5[CartesianIndex(10)])
@test all(cs35[1:2, 1:1, 1:2, :, :] .== data5[1:2, 1:1, 1:2, :, :])
@test all(cs35[1:2, 1:1, 1, :, :] .== data5[1:2, 1:1, 1, :, :])

cs35[1:2,1,1,1,2] .= -1.0
@test all(cs35[1:2,1,1,1,2] .== -1.0)
@test all(cs35[1:2,1,1,1,2] .== data5[1:2,1,1,1,2])

#@code_warntype cs35[1:2, 1:1, 2, 1]


text:
module Kiltro

using LinearAlgebra
using Printf

export QuadratureRule, FunctionSpace, Solve

include("Utils.jl")

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

function NodeArrangementHex(C)
  """It manages the order of nodes and faces in an Hexaedral element."""
  linear_bases_idx = [1,(C+2),(C+2)^2,(C+2)^2-(C+1)]
  quad_aranger = [linear_bases_idx; delete1d(collect(range(1,stop=(C+2)^2)),linear_bases_idx)]
  element_numbering = copy(quad_aranger)
  for i in 0:C
    faces_z = quad_aranger.+maximum(element_numbering)
    element_numbering = [element_numbering; faces_z]
    if i==C
      element_numbering = element_numbering[.~in1d(element_numbering,faces_z[1:4])]
      element_numbering = [element_numbering[1:4];faces_z[1:4];element_numbering[5:end]]
    end
  end
  traversed_edge_numbering_hex = nothing

  # GET FACE NUMBERING ORDER FROM TETRAHEDRAL ELEMENT
  face_0,face_1,face_2,face_3 = [],[],[],[]
  if C==0
    face_0 = [0,1,2,3]  # constant Z =-1 plane
    face_1 = [4,5,6,7]  # constant Z = 1 plane
    face_2 = [0,1,5,4]  # constant Y =-1 plane
    face_3 = [3,2,6,7]  # constant Y = 1 plane
    face_4 = [0,3,7,4]  # constant X =-1 plane
    face_5 = [1,2,6,5]  # constant X = 1 plane

    face_numbering = [face_0; face_1; face_2; face_3; face_4; face_5] .+ 1
    face_numbering = transpose(reshape(face_numbering,4,6))
  else
    # THIS IS A FLOATING POINT BASED ALGORITHM
    error("Polynomial degreee not implemented yet")
  end

  return face_numbering, traversed_edge_numbering_hex, element_numbering
end

function Material(mu,kappa,F)
  """It computes the stress and elasticity of the material from the deformation gradient"""
  I = [1.0 0.0 0.0; 0.0 1.0 0.0; 0.0 0.0 1.0]
  J = det(F)
  b = F*transpose(F)
  trb = tr(b)

  energy = 0.5*mu*(trb -3.0)

  stress = mu*J^(-5.0/3.0)*(b - 1.0/3.0*trb*I) + kappa*(J-1.0)*I

  II_ijkl = einsum("ijkl",I,I)
  II_ikjl = einsum("ikjl",I,I)
  II_iljk = einsum("iljk",I,I)
  bI_ijkl = einsum("ijkl",b,I)
  Ib_ijkl = einsum("ijkl",I,b)
  elasticity = 2*mu*J^(-5.0/3.0)*(1.0/9.0*trb*II_ijkl - 1.0/3.0*(bI_ijkl + Ib_ijkl) + 1.0/6.0*trb*(II_ikjl + II_iljk) )
  elasticity += kappa*((2.0*J-1.0)*II_ijkl - (J-1.0)*(II_ikjl + II_iljk))
  H_Voigt = Voigt(elasticity)

  return stress, H_Voigt
end

function QuadratureRule(N,a,b)
  """Return integration points and weights from Gauss-Legendre polynomials.
     This is a gaussian quadrature rule."""
  deps = eps(1.0)
  x = zeros(N)
  w = zeros(N)
  M = div((N+1),2)
  xm = 0.5*(b+a)
  xl = 0.5*(b-a)
  for i in 1:M
    z = cos(pi*(i-0.25)/(N+0.5))
    z1 = 2.0
    while abs(z-z1) > deps
      global pp
      p1 = 1.0
      p2 = 0.0
      for j in 1:N
        p3 = p2
        p2 = p1
        p1 = ((2.0*j-1.0)*z*p2-(j-1.0)*p3)/j
      end
      pp = N*(z*p1-p2)/(z*z-1.0)
      z1 = z
      z = z1-p1/pp
    end
    x[i] = xm-xl*z
    x[N+1-i] = xm+xl*z
    w[i] = 2.0*xl/((1.0-z*z)*pp*pp)
    w[N+1-i] = w[i]
  end
  if size(x,1) == size(w,1)
    N = size(x,1)
    x_flatten = zeros(Int(N*N*N),3)
    w_flatten = zeros(Int(N*N*N))
    counter = 1
    for i in 1:N
      for j in 1:N
        for k in 1:N
          w_flatten[counter] = w[i]*w[j]*w[k]
          x_flatten[counter,1] = x[i]
          x_flatten[counter,2] = x[j]
          x_flatten[counter,3] = x[k]
          counter += 1
        end
      end
    end
  end
  return x_flatten, w_flatten
end

function OneDLagrange(C,xi)
  """Returns a Lagrange interpolation (Base) and its derivatives given the
     interpolation order C and the point xi."""
  n = C+2
  ranger = collect(range(0,stop=(n-1)))
  eps = collect(range(-1.0,1.0,length=n))

  A = zeros(n,n)
  A[:,1] .= 1.0
  for i in 2:n
    A[:,i] = eps.^(i-1)
  end

  RHS = zeros(n,n)
  for i in 1:n
    RHS[i,i] = 1
  end
  coeff = inv(A)*RHS
  xis = ones(n).*xi.^ranger
  N = transpose(coeff)*xis
  dN = transpose(coeff[2:end,:]) * xis[1:end-1].*(1 .+ ranger[1:end-1])

  return N, dN, eps
end

function HexLagrange(C,zeta,eta,beta)
  """Computes C order Lagrangian bases with equally spaced points
     from 1D Lagrange interpolation"""

  Neta = zeros(C+2); Nzeta = zeros(C+2); Nbeta = zeros(C+2)

  Nzeta .= OneDLagrange(C,zeta)[1]
  Neta .=  OneDLagrange(C,eta)[1]
  Nbeta .=  OneDLagrange(C,beta)[1]

  # Ternsorial product
  node_arranger = NodeArrangementHex(C)[3]
  Bases = zeros((C+2)*(C+2)*(C+2))
  for i in 1:(C+2)
    for j in 1:(C+2)
      for k in 1:(C+2)
        Bases[(i-1)*(C+2)*(C+2)+(j-1)*(C+2)+k] = Nbeta[i]*Neta[j]*Nzeta[k]
      end
    end
  end
  Bases = Bases[(node_arranger)]

  return Bases
end

function HexGradLagrange(C,zeta,eta,beta)
  """Computes gradient of C order Lagrangian bases with equally spaced points
     from 1D Lagrange interpolation"""

  gBases = zeros((C+2)^3,3)
  Nzeta = zeros(C+2); Neta = zeros(C+2); Nbeta = zeros(C+2)
  gNzeta = zeros(C+2); gNeta = zeros(C+2);  gNbeta = zeros(C+2)
  # Compute each from one-dimensional bases
  Nzeta = OneDLagrange(C,zeta)[1]
  Neta = OneDLagrange(C,eta)[1]
  Nbeta = OneDLagrange(C,beta)[1]
  gNzeta = OneDLagrange(C,zeta)[2]
  gNeta = OneDLagrange(C,eta)[2]
  gNbeta = OneDLagrange(C,beta)[2]

  # Ternsorial product
  node_arranger = NodeArrangementHex(C)[3]
  g0 = zeros((C+2)*(C+2)*(C+2))
  g1 = zeros((C+2)*(C+2)*(C+2))
  g2 = zeros((C+2)*(C+2)*(C+2))
  for i in 1:(C+2)
    for j in 1:(C+2)
      for k in 1:(C+2)
        g0[(i-1)*(C+2)*(C+2)+(j-1)*(C+2)+k] = Nbeta[i]*Neta[j]*gNzeta[k]
        g1[(i-1)*(C+2)*(C+2)+(j-1)*(C+2)+k] = Nbeta[i]*gNeta[j]*Nzeta[k]
        g2[(i-1)*(C+2)*(C+2)+(j-1)*(C+2)+k] = gNbeta[i]*Neta[j]*Nzeta[k]
      end
    end
  end
  gBases[:,1] = g0[node_arranger]
  gBases[:,2] = g1[node_arranger]
  gBases[:,3] = g2[node_arranger]

  return gBases
end

function FunctionSpace(z,w)
  """Collects the Bases and its gradient from quadrature rule to build
     the function space for the element."""
  C = 0
  ndim = 3
  ns = Int((C+2)^ndim)
  Basis = zeros(ns,size(w,1))
  gBasisx = zeros(ns,size(w,1))
  gBasisy = zeros(ns,size(w,1))
  gBasisz = zeros(ns,size(w,1))

  for i in 1:size(w,1)
    ndummy = HexLagrange(C,z[i,1],z[i,2],z[i,3])
    dummy = HexGradLagrange(C,z[i,1],z[i,2],z[i,3])
    Basis[:,i] = ndummy
    gBasisx[:,i] = dummy[:,1]
    gBasisy[:,i] = dummy[:,2]
    gBasisz[:,i] = dummy[:,3]
  end
  Jm = zeros(ndim,size(Basis,1),size(w,1))
  AllGauss = zeros(size(w,1))
  for counter in 1:size(w,1)
    # GRADIENT TENSOR IN PARENT ELEMENT [\nabla_\varepsilon (N)]
    Jm[1,:,counter] = gBasisx[:,counter]
    Jm[2,:,counter] = gBasisy[:,counter]
    Jm[3,:,counter] = gBasisz[:,counter]

    AllGauss[counter] = w[counter]
  end
  return AllGauss, Basis, Jm
end

function KinematicMeasures(AllGauss,Jm,LagrangeElemCoords,EulerElemCoords)
  """Given the function space, Lagrangian and Eulerian coordinates returns
     the kinematics measures (deformation measures) such as the deformation gradient."""

  #ParentGradientX = np.einsum('ijk,jl->kil', Jm, LagrangeElemCoords)
  ParentGradientX = zeros(size(Jm,3),size(Jm,1),size(LagrangeElemCoords,2))
  for i in 1:size(Jm,1)
    for j in 1:size(Jm,2)
      for k in 1:size(Jm,3)
        for l in 1:size(LagrangeElemCoords,2)
          ParentGradientX[k,i,l] += Jm[i,j,k]*LagrangeElemCoords[j,l]
        end
      end
    end
  end
  #MaterialGradient = np.einsum('ijk,kli->ijl', inv(ParentGradientX), Jm)
  invParentGradientX = zeros(size(ParentGradientX,1),size(ParentGradientX,2),size(ParentGradientX,3))
  for i in 1:size(ParentGradientX,1)
    invParentGradientX[i,:,:] = inv(ParentGradientX[i,:,:])
  end
  MaterialGradient = zeros(size(ParentGradientX,1),size(ParentGradientX,2),size(Jm,2))
  for i in 1:size(invParentGradientX,1)
    for j in 1:size(invParentGradientX,2)
      for k in 1:size(Jm,1)
        for l in 1:size(Jm,2)
          MaterialGradient[i,j,l] += invParentGradientX[i,j,k]*Jm[k,l,i]
        end
      end
    end
  end
  #F = np.einsum('ij,kli->kjl', EulerElemCoords, MaterialGradient)
  F = zeros(size(MaterialGradient,1),size(EulerElemCoords,2),size(MaterialGradient,2))
  for i in 1:size(EulerElemCoords,1)
    for j in 1:size(EulerElemCoords,2)
      for k in 1:size(MaterialGradient,1)
        for l in 1:size(MaterialGradient,2)
          F[k,j,l] += EulerElemCoords[i,j]*MaterialGradient[k,l,i]
        end
      end
    end
  end

  #ParentGradientx = np.einsum('ijk,jl->kil',Jm, EulerElemCoords)
  ParentGradientx = zeros(size(Jm,3),size(Jm,1),size(EulerElemCoords,2))
  for i in 1:size(Jm,1)
    for j in 1:size(Jm,2)
      for k in 1:size(Jm,3)
        for l in 1:size(EulerElemCoords,2)
          ParentGradientx[k,i,l] += Jm[i,j,k]*EulerElemCoords[j,l]
        end
      end
    end
  end
  #SpatialGradient = np.einsum('ijk,kli->ilj',inv(ParentGradientx), Jm)
  invParentGradientx = zeros(size(ParentGradientx,1),size(ParentGradientx,2),size(ParentGradientx,3))
  for i in 1:size(ParentGradientx,1)
    invParentGradientx[i,:,:] = inv(ParentGradientx[i,:,:])
  end
  SpatialGradient = zeros(size(ParentGradientx,1),size(Jm,2),size(ParentGradientx,2))
  for i in 1:size(invParentGradientx,1)
    for j in 1:size(invParentGradientx,2)
      for k in 1:size(Jm,1)
        for l in 1:size(Jm,2)
          SpatialGradient[i,l,j] += invParentGradientx[i,j,k]*Jm[k,l,i]
        end
      end
    end
  end
  #detJ = np.einsum('i,i,i->i',AllGauss[:,0],np.abs(det(ParentGradientX)),np.abs(StrainTensors['J']))
  detJ = zeros(size(AllGauss,1))
  for i in 1:size(AllGauss,1)
    detParentGradientX = abs(det(ParentGradientX[i,:,:]))
    J = det(F[i,:,:])
    detJ[i] = AllGauss[i]*detParentGradientX*J
  end
  return SpatialGradient,F,detJ
end

function FillConstitutiveB(SpatialGradient,ndim,nvar,VoigtSize)
  """It fills the B-matrix (or N-gradient) to multiply the the elasticity tensor for
     its integration in the element."""

  nodeperelem = size(SpatialGradient,2)
  B = zeros(nodeperelem*nvar,VoigtSize)

  for node in 1:nodeperelem
    B[nvar*node-2,1] = SpatialGradient[1,node]
    B[nvar*node-1,2] = SpatialGradient[2,node]
    B[nvar*node,3] = SpatialGradient[3,node]

    B[nvar*node-1,6] = SpatialGradient[3,node]
    B[nvar*node,6] = SpatialGradient[2,node]

    B[nvar*node-2,5] = SpatialGradient[3,node]
    B[nvar*node,5] = SpatialGradient[1,node]

    B[nvar*node-2,4] = SpatialGradient[2,node]
    B[nvar*node-1,4] = SpatialGradient[1,node]
  end
  return B
end

function FillGeometricB(SpatialGradient,CauchyStress,ndim,nvar)
  """It fills the geometric B-matrix to be multiply with the stress
     to get the geometric stiffness matrix from integration."""

  nodeperelem = size(SpatialGradient,2)
  B = zeros(nvar*nodeperelem,ndim*ndim)
  S = zeros(ndim*ndim,ndim*ndim)

  for node in 1:nodeperelem
    B[nvar*node-2,1] = SpatialGradient[1,node]
    B[nvar*node-2,2] = SpatialGradient[2,node]
    B[nvar*node-2,3] = SpatialGradient[3,node]
    B[nvar*node-1,4] = SpatialGradient[1,node]
    B[nvar*node-1,5] = SpatialGradient[2,node]
    B[nvar*node-1,6] = SpatialGradient[3,node]
    B[nvar*node,7] = SpatialGradient[1,node]
    B[nvar*node,8] = SpatialGradient[2,node]
    B[nvar*node,9] = SpatialGradient[3,node]
  end

  S[1,1] = CauchyStress[1,1]
  S[1,2] = CauchyStress[1,2]
  S[1,3] = CauchyStress[1,3]
  S[2,1] = CauchyStress[2,1]
  S[2,2] = CauchyStress[2,2]
  S[2,3] = CauchyStress[2,3]
  S[3,1] = CauchyStress[3,1]
  S[3,2] = CauchyStress[3,2]
  S[3,3] = CauchyStress[3,3]

  S[4,4] = CauchyStress[1,1]
  S[4,5] = CauchyStress[1,2]
  S[4,6] = CauchyStress[1,3]
  S[5,4] = CauchyStress[2,1]
  S[5,5] = CauchyStress[2,2]
  S[5,6] = CauchyStress[2,3]
  S[6,4] = CauchyStress[3,1]
  S[6,5] = CauchyStress[3,2]
  S[6,6] = CauchyStress[3,3]

  S[7,7] = CauchyStress[1,1]
  S[7,8] = CauchyStress[1,2]
  S[7,9] = CauchyStress[1,3]
  S[8,7] = CauchyStress[2,1]
  S[8,8] = CauchyStress[2,2]
  S[8,9] = CauchyStress[2,3]
  S[9,7] = CauchyStress[3,1]
  S[9,8] = CauchyStress[3,2]
  S[9,9] = CauchyStress[3,3]

  return B,S
end

function GetTotalTraction(CauchyStress)
  """Returns the stress as a vector into Voigt notation."""
  TotalTraction = zeros(6)
  TotalTraction[1] = CauchyStress[1,1]
  TotalTraction[2] = CauchyStress[2,2]
  TotalTraction[3] = CauchyStress[3,3]
  TotalTraction[4] = CauchyStress[1,2]
  TotalTraction[5] = CauchyStress[1,3]
  TotalTraction[6] = CauchyStress[2,3]
  return TotalTraction
end

function ConstitutiveStiffnessIntegrand(SpatialGradient, CauchyStress, H_Voigt)
  """Applies to displacement based formulation, to integrate the constitutive stiffness"""
  ndim = 3
  nvar = 3
  VoigtSize = size(H_Voigt,1)
  # SpatialGradient(ndim x nodesperelem) and B(nodesperelem*nvar,H_VoigtSize)
  SpatialGradient = transpose(SpatialGradient)
  B = FillConstitutiveB(SpatialGradient,ndim,nvar,VoigtSize)

  BDB = B*(H_Voigt*(transpose(B)))

  t=zeros(size(B,1))
  TotalTraction = GetTotalTraction(CauchyStress)
  t = B*TotalTraction

  return BDB, t
end

function GeometricStiffnessIntegrand(SpatialGradient, CauchyStress,ndim,nvar)
    """Applies to displacement based, displacement potential based and all mixed
    formulations that involve static condensation"""

  SpatialGradient = transpose(SpatialGradient)

  B,S = FillGeometricB(SpatialGradient,CauchyStress,ndim,nvar)
  BDB = (B*S)*transpose(B)

  return BDB
end

function Assemble(AllGauss,Basis,Jm,points,elements,mu,kappa,Eulerx,ndim,nvar)
  """It assembles the stiffness and internal forces by element."""

  nodeperelem = size(Basis,1)
  ranger = collect(1:1:nodeperelem)

  stiffness_global = zeros(size(points,1)*nvar,size(points,1)*nvar)
  traction = zeros(size(points,1)*nvar)

  # loop on elements
  for elem in 1:size(elements,1)
    LagrangeElemCoords = points[elements[elem,:],:]
    EulerElemCoords = Eulerx[elements[elem,:],:]

    stiffness = zeros(nodeperelem*nvar,nodeperelem*nvar)
    tractionforce = zeros(nodeperelem*nvar)
    SpatialGradient,F,detJ = KinematicMeasures(AllGauss,Jm,LagrangeElemCoords,EulerElemCoords)

    for counter in 1:size(AllGauss,1)
      # COMPUTE CAUCHY STRESS AND HESSIAN FROM MATERIAL AT GAUSS POINT
      CauchyStress,H_Voigt = Material(mu,kappa,F[counter,:,:])
      # COMPUTE THE TANGENT STIFFNESS MATRIX
      BDB_1, t = ConstitutiveStiffnessIntegrand(SpatialGradient[counter,:,:], CauchyStress, H_Voigt)
      # COMPUTE GEOMETRIC STIFFNESS MATRIX
      BDB_1 += GeometricStiffnessIntegrand(SpatialGradient[counter,:,:],CauchyStress,ndim,nvar)
      # INTEGRATE TRACTION FORCE
      tractionforce += t*detJ[counter]
      # INTEGRATE STIFFNESS
      stiffness += BDB_1*detJ[counter]
    end

    for i in 0:(nvar-1)
      for j in 0:(nvar-1)
        stiffness_global[nvar.*elements[elem,:].-i,nvar.*elements[elem,:].-j] += stiffness[nvar.*ranger.-i,nvar.*ranger.-j]
      end
      traction[nvar.*elements[elem,:].-i] += tractionforce[nvar.*ranger.-i]
    end
  end

  return stiffness_global,traction
end

function NewtonRaphson(Increment,AllGauss,Basis,Jm,points,elements,mu,kappa,
                       norm_residual,NormForces,K,Residual,NodalForces,Eulerx,columns_in,ndim,nvar)
  """This is the Newton-Raphson loop of iterations to find a solution"""

  Tolerance = 1.0e-6
  Iter = 1

  while norm_residual > Tolerance

    # GET REDUCED FORCE VECTOR
    F_b = Residual[columns_in]
    # GET REDUCED STIFFNESS MATRIX
    K_b = K[columns_in,columns_in]

    sol = -inv(K_b)*F_b

    # UPDATE THE FREE DOFs
    # GET TOTAL SOLUTION
    TotalSol = zeros(size(K,1),1)
    TotalSol[columns_in] = sol
    # RE-ORDER SOLUTION COMPONENTS
    dU = transpose(reshape(TotalSol,nvar,Int(size(TotalSol,1)/nvar)))

    # UPDATE THE EULERIAN COMPONENTS, UPDATE THE GEOMETRY
    Eulerx += dU[:,1:ndim]

    # RE-ASSEMBLE - COMPUTE STIFFNESS AND INTERNAL TRACTION FORCES
    K,TractionForces = Assemble(AllGauss,Basis,Jm,points,elements,mu,kappa,Eulerx,ndim,nvar)

    # FIND THE RESIDUAL
    Residual[columns_in] = TractionForces[columns_in] - NodalForces[columns_in]

    # SAVE THE NORM
    abs_norm_residual = norm(Residual[columns_in])
    if Iter==1
      NormForces = norm(Residual[columns_in])
    end
    norm_residual = abs(norm(Residual[columns_in])/NormForces)

    @printf("Iteration %d for increment %d. Residual (abs) %e\t Residual (rel) %e\n",Iter,Increment,abs_norm_residual,norm_residual)
    #{0:>16.7g}

    # BREAK BASED ON RELATIVE NORM
    if abs(abs_norm_residual) < Tolerance
      break
    end

    # BREAK BASED ON INCREMENTAL SOLUTION - KEEP IT AFTER UPDATE
    if norm(dU) <=  Tolerance*0.1
      @printf("Incremental solution within tolerance i.e. norm(dU): %f\n",norm(dU))
      break
    end

    # UPDATE ITERATION NUMBER
    Iter += 1

    # BREAK BASED ON ITERATION LIMIT NUMBER
    if Iter==10
      break
    end

  end
  return Eulerx, K, Residual
end

function StaticSolver(LoadIncrements,AllGauss,Basis,Jm,points,elements,mu,kappa,
                      applied_dirichlet,NeumannForces,
                      K,Residual,NodalForces,TotalDisp,Eulerx,
                      columns_in,columns_out,ndim,nvar)
  """Static solver"""

  NormForces = 0.0
  LoadFactor = 1.0/LoadIncrements
  AppliedDirichletInc = zeros(size(applied_dirichlet,1))
  LoadFactorInc = 0.0

  for Increment in 1:LoadIncrements

    LoadFactorInc += LoadFactor

    DeltaF = LoadFactor*NeumannForces
    NodalForces += DeltaF
    # OBRTAIN INCREMENTAL RESIDUAL - CONTRIBUTION FROM BOTH NEUMANN AND DIRICHLET
    F = zeros(size(Residual))
    nnz_cols = .~isapprox.(applied_dirichlet,0.0)
    if size(columns_out[nnz_cols],1)==0
      F[columns_in] = F[columns_in]
    else
      F[columns_in] = F[columns_in] - (K[columns_in,columns_out[nnz_cols]]*applied_dirichlet[nnz_cols])
    end
    Residual = -F
    Residual -= DeltaF
    AppliedDirichletInc = LoadFactor*applied_dirichlet

    # LET NORM OF THE FIRST RESIDUAL BE THE NORM WITH RESPECT TO WHICH WE
    # HAVE TO CHECK THE CONVERGENCE OF NEWTON RAPHSON. TYPICALLY THIS IS
    # NORM OF NODAL FORCES
    if Increment==1
      NormForces = norm(Residual)
      # AVOID DIVISION BY ZERO
      if isapprox(NormForces,0.0)
        NormForces = 1e-14
      end
    end

    norm_residual = norm(Residual)/NormForces

    # APPLY INCREMENTAL DIRICHLET PER LOAD STEP (THIS IS INCREMENTAL NOT ACCUMULATIVE)
    # GET TOTAL SOLUTION
    TotalSol = zeros(size(K,1),1)
    TotalSol[columns_out] = AppliedDirichletInc
    # RE-ORDER SOLUTION COMPONENTS
    IncDirichlet = transpose(reshape(TotalSol,nvar,Int(size(TotalSol,1)/nvar)))
    # UPDATE EULERIAN COORDINATE
    Eulerx += IncDirichlet[:,1:ndim]

    #Newton-Raphson loop of iterations
    Eulerx,K,Residual = NewtonRaphson(Increment,AllGauss,Basis,Jm,points,elements,mu,kappa,norm_residual,NormForces,
                                      K,Residual,NodalForces,Eulerx,columns_in,ndim,nvar)

    # UPDATE DISPLACEMENTS FOR THE CURRENT LOAD INCREMENT
    TotalDisp[:,1:ndim,Increment] = Eulerx - points

    @printf("\nFinished Load increment %d\n",Increment)

    MaxDisp = [maximum(TotalDisp[:,1,Increment]),maximum(TotalDisp[:,2,Increment]),maximum(TotalDisp[:,3,Increment])]
    MinDisp = [minimum(TotalDisp[:,1,Increment]),minimum(TotalDisp[:,2,Increment]),minimum(TotalDisp[:,3,Increment])]
    #@printf("\nMinimum and maximum incremental solution values at increment are %e %e %e \n",MaxDisp)
    println("Minimum incremental solution values ",MinDisp)
    println("Maximum incremental solution values ",MaxDisp)
    println()

  end
  return TotalDisp
end

function Solve(AllGauss,Basis,Jm,points,elements,mu,kappa,dirichlet_flags,neumann_flags,ndim,nvar)
  """The heart of the solver."""

  # It runs the solver to seek the solution under a Newton-Raphson scheme of iteration.
  println("Pre-processing the information. Getting paths, solution parameters, mesh info, interpolation info etc...")
  println("Number of nodes is ",size(points,1)," and number of DoFs is ",size(points,1)*nvar)

  LoadIncrements = 3
  # INITIATE DATA FOR THE ANALYSIS
  NodalForces = zeros(size(points,1)*nvar)
  Residual = zeros(size(points,1)*nvar)
  TotalDisp = zeros(size(points,1),nvar,LoadIncrements)

  println("Assembling the system and acquiring neccessary information for the analysis...")

  # APPLY DIRICHELT BOUNDARY CONDITIONS AND GET DIRICHLET RELATED FORCES
  flat_dirich = vec(transpose(dirichlet_flags))
  columns_out = collect(1:size(flat_dirich,1))[.~isnan.(flat_dirich)] #Int64
  applied_dirichlet = flat_dirich[.~isnan.(flat_dirich)]
  # GENERAL PROCEDURE - GET REDUCED MATRICES FOR FINAL SOLUTION
  columns_in = delete1d(collect(1:nvar*size(points,1)),columns_out)

  # FIND PURE NEUMANN (EXTERNAL) NODAL FORCE VECTOR
  NeumannForces = zeros(size(points,1)*nvar)
  flat_neu = vec(transpose(neumann_flags))
  to_apply = collect(1:size(flat_neu,1))[.~isnan.(flat_neu)]
  applied_neumann = flat_neu[.~isnan.(flat_neu)]
  NeumannForces[to_apply] = applied_neumann

  # initialize euler coordinates
  Eulerx = points

  # assembly
  K,TractionForces = Assemble(AllGauss,Basis,Jm,points,elements,mu,kappa,Eulerx,ndim,nvar)

  TotalDsip = StaticSolver(LoadIncrements,AllGauss,Basis,Jm,points,elements,mu,kappa,
                           applied_dirichlet,NeumannForces,
                           K,Residual,NodalForces,TotalDisp,Eulerx,
                           columns_in,columns_out,ndim,nvar)

  return TotalDisp
end

end

max_stars_repo_name: jdlaubrie/Kiltro

max_issues_repo_name: jdlaubrie/Kiltro

max_forks_repo_name: jdlaubrie/Kiltro


text:
"""
Black-box system under test (SUT). In this case,
we use an MNIST neural network classifier.
"""
module SystemUnderTest

export accuracy, train, trainandsave, load, getdata

using Base.Iterators: repeated
using BSON: @save, @load
using Colors
using CUDA
using Images
using MLDatasets
using NNlib
using Parameters
using Flux
using Flux.Data: DataLoader
using Flux: @epochs, onehotbatch, onecold, logitcrossentropy, throttle
using Statistics

@with_kw mutable struct Args
    α::Float64 = 3e-4      # learning rate
    batchsize::Int = 1024  # batch size
    epochs::Int = 20       # number of epochs
    device::Function = gpu # set as gpu if available
    throttle::Int = 1      # throttle print every X seconds
end

function getdata(args=Args())
    # load dataset
    xtrain, ytrain = MLDatasets.MNIST.traindata(Float32)
    xtest, ytest = MLDatasets.MNIST.testdata(Float32)

    # reshape data to flatten each image into a vector
    xtrain = Flux.flatten(xtrain)
    xtest = Flux.flatten(xtest)

    # one-hot encode the labels
    ytrain, ytest = onehotbatch(ytrain, 0:9), onehotbatch(ytest, 0:9)

    # split into batches
    traindata = DataLoader(xtrain, ytrain, batchsize=args.batchsize, shuffle=true)
    testdata = DataLoader(xtest, ytest, batchsize=args.batchsize)

    return traindata, testdata
end

function buildmodel(; imgsize=(28,28,1), nclasses=10)
    return Chain(
            Dense(prod(imgsize), 32, relu),
            Dense(32, nclasses))
end

"""
Cost function:

\$\\mathcal{J}(\\mathbf{\\hat{y}}, \\mathbf{y}) = \\frac{1}{m}\\sum \\mathcal{L}(\\hat{y}, y)\$

with loss:

\$\\mathcal{L}(\\hat{y}, y) = -\\frac{1}{n}\\sum_{i=1}^n y \\left(\\hat{y} - \\log\\left(\\sum e^{\\hat{y}}\\right)\\right)\$
"""
cost(dataloader, model) = mean(logitcrossentropy(model(x), y) for (x,y) in dataloader)

function accuracy(dataloader, model)
    acc = 0
    for (x,y) in dataloader
        acc += sum(onecold(cpu(model(x))) .== onecold(cpu(y))) / size(x,2)
    end
    return acc/length(dataloader)
end

function train(; kwargs...)
    # initialize model parameters
    args = Args(; kwargs...)

    # load data
    traindata, testdata = getdata(args)

    # construct model
    m = buildmodel()
    traindata = args.device.(traindata)
    testdata = args.device.(testdata)
    m = args.device(m)
    loss(x,y) = logitcrossentropy(m(x), y)

    # callback, optimizer, and training
    callback = throttle(()->@show(cost(traindata, m)), args.throttle)
    opt = ADAM(args.α)

    @epochs args.epochs Flux.train!(loss, params(m), traindata, opt, cb=callback)

    @show accuracy(traindata, m)
    @show accuracy(testdata, m)

    return m, traindata, testdata
end

function predict(m, x)
    return onecold(cpu(m(x)))
end

function trainandsave(; kwargs...)
    m, traindata, testdata = train(; kwargs...)
    m = cpu(m)
    @save "models/sut.bson" m
end

"""
Load pre-trained neural network model.

```julia
using Flux
using CUDA
m = load(\"models/sut.bson\")
```
"""
function load(model_filename_bson)
    @load model_filename_bson m
    return m
end

end # module


text:
using TableSchema
using Test
using Dates

import TableSchema: read, validate, infer, commit, save

import DelimitedFiles: readdlm

include("read.jl")
include("schema.jl")
include("validate.jl")
include("infer.jl")
include("edit.jl")
include("save.jl")
# include("changes.jl")


text:
import Pkg; Pkg.add("PackageCompiler")
import PackageCompiler

PackageCompiler.create_app(
    ".",
    "pangraph",
    executables=["pangraph"=>"julia_main"],
    force=true,
    include_lazy_artifacts=false,
    precompile_execution_file="trace.jl"
)

max_stars_repo_name: neherlab/pangraph

max_issues_repo_name: neherlab/pangraph

max_forks_repo_name: neherlab/pangraph


text:
function next_time_period(TIME::Timer, QC::Array{QuayCrane, 1})
    TIME.period += 1
    l = Array{Int, 1}()
    #update all QC atributtes
    for qc in QC
        #update QC position
        if qc.status == "moving"
            if  qc.next_bay < qc.current_bay
                qc.current_bay -= CTS.tt
            elseif qc.current_bay < qc.next_bay
                qc.current_bay += CTS.tt
            end
        end
        #update QC time left and status
        if qc.time_left > 1
            qc.time_left -= 1
        elseif qc.time_left == 1
            if length(qc.task_buffer) > 0 && qc.status != "waiting to move"
                qc.time_left = 2*qc.task_buffer[1].t
                qc.status = "loading"
                deleteat!(qc.task_buffer, 1)
            elseif length(qc.task_buffer) > 0 && qc.status == "waiting to move"
                qc.time_left = travel_time(qc.current_bay, qc.task_buffer[1].b, CTS)
                qc.status = "moving"
            else
                qc.time_left = 0
                qc.status = "idle"
            end
        end
        if qc.status == "idle"
            push!(l, qc.q)
        end
    end
    TIME.available_cranes = l
    #println(TIME)
end

function travel_time(current_bay::Int, target_bay::Int, CTS::Constants)
    return(abs(current_bay - target_bay)*CTS.tt)
end

function update_quay_crane(TIME::Timer, QC::Array{QuayCrane, 1}, task::Task, q::Int, CTS::Constants)
    if task.b == QC[q].current_bay
        QC[q].time_left = 2*task.t
        QC[q].status = "loading"
        QC[q].next_bay = QC[q].current_bay
    else
        QC[q].time_left = travel_time(QC[q].current_bay, task.b, CTS)
        QC[q].status = "moving"
        QC[q].next_bay = task.b
        if task.t != 0 && task.p !=0
            push!(QC[q].task_buffer, task)
        end
    end
end

function update_load_seq(TIME::Timer, LS::LoadingSequence, QC::Array{QuayCrane, 1}, task::Task, q::Int, CTS::Constants)
    LS.len += 1
    LS.tasks_left -= 1
    start_time = TIME.period + travel_time(QC[q].current_bay, task.b, CTS)
    push!(LS.order, (task=task, start_time=start_time, qc=q))
    push!(LS.filled_pos, task.p)
    push!(LS.loaded_cont,task.c)
end

function update_wait_quay_crane(TIME::Timer, QC::Array{QuayCrane, 1}, task::Task, q::Int, move_q::Int, move_bay::Int, CTS::Constants)
    if task.b == QC[q].current_bay
        QC[q].time_left = travel_time(QC[move_q].current_bay, QC[move_q].next_bay, CTS)
        QC[q].status = "waiting to start"
        QC[q].next_bay = QC[q].current_bay
        push!(QC[q].task_buffer, task)
    elseif travel_time(QC[q].current_bay, task.b, CTS) >= travel_time(QC[move_q].current_bay, QC[move_q].next_bay, CTS)
        QC[q].time_left = travel_time(QC[q].current_bay, task.b, CTS)
        QC[q].status = "moving"
        QC[q].next_bay = task.b
        if task.t != 0 && task.p !=0
            push!(QC[q].task_buffer, task)
        end
    else
        QC[q].time_left = travel_time(QC[move_q].current_bay, QC[move_q].next_bay, CTS) - travel_time(QC[q].current_bay, task.b, CTS)
        QC[q].status = "waiting to move"
        QC[q].next_bay = task.b
        if task.t != 0 && task.p !=0
            push!(QC[q].task_buffer, task)
        end
    end
end

function update_wait_load_seq(TIME::Timer, LS::LoadingSequence, QC::Array{QuayCrane, 1}, task::Task, q::Int, move_q::Int, move_bay::Int, CTS::Constants)
    LS.len += 1
    LS.tasks_left -= 1
    start_time = TIME.period + max(travel_time(QC[q].current_bay, task.b, CTS), travel_time(QC[move_q].current_bay, QC[move_q].next_bay, CTS))
    push!(LS.order, (task=task, start_time=start_time, qc=q))
    push!(LS.filled_pos, task.p)
    push!(LS.loaded_cont,task.c)
end

function add_task(task::Task, q::Int, TIME::Timer, LS::LoadingSequence, QC::Array{QuayCrane, 1}, CTS::Constants)
    update_load_seq(TIME, LS, QC, task, q, CTS)
    update_quay_crane(TIME, QC, task, q, CTS)
    deleteat!(TIME.available_cranes, findall(x->x==q, TIME.available_cranes))
end

function add_task_move(task::Task, q::Int, move_q::Int, move_bay::Int, move_status::String, LS::LoadingSequence, TIME::Timer, QC::Array{QuayCrane, 1}, CTS::Constants)
    if move_status == "idle"
        update_quay_crane(TIME, QC, Task(0, move_bay, 0, 0), move_q, CTS)
        update_wait_quay_crane(TIME, QC, task, q, move_q, move_bay, CTS)
        update_wait_load_seq(TIME, LS, QC, task, q, move_q, move_bay, CTS)
        deleteat!(TIME.available_cranes, findall(x->x==q, TIME.available_cranes))
        deleteat!(TIME.available_cranes, findall(x->x==move_q, TIME.available_cranes))
    elseif move_status == "moving"
        update_wait_quay_crane(TIME, QC, task, q, move_q, move_bay, CTS)
        update_wait_load_seq(TIME, LS, QC, task, q, move_q, move_bay, CTS)
        deleteat!(TIME.available_cranes, findall(x->x==q, TIME.available_cranes))
        deleteat!(TIME.available_cranes, findall(x->x==move_q, TIME.available_cranes))
    end
end

max_stars_repo_licenses: Zlib MIT

max_issues_repo_licenses: Zlib MIT

max_forks_repo_licenses: Zlib MIT


text:
# Low Rank Approximation

mutable struct LRA{T}
	om::LinearModel{T}
	optm::LRAOpt{T}
	gx::X{T}
	sx::X{T}
	opG::LinearMaps.LinearMap{T}
	opS::LinearMaps.LinearMap{T}
	err::DataFrames.DataFrame
end

"""
Constructor for the blind deconvolution problem
"""
function LRA(nt, nr, nra;
	       dobs=nothing, gobs=nothing, sobs=nothing, verbose=false, attrib_inv=:g,
	       )
	# determine type of IBD
	if(!(dobs===nothing))
		(size(dobs)≠(nt,nr)) && error("dobs size error")
		T=eltype(dobs)
	else
		T1=eltype(sobs)
		T2=eltype(gobs)
		(T1≠T2) ? error("type difference") : (T=T1)
		(size(sobs)≠(nt,nra)) && error("sobs size error")
		(size(gobs)≠(nra,nr)) && error("gobs size error")
	end

	# store observed data
	om=LinearModel(nt,nr,nra,T, d=dobs, g=gobs, s=sobs)

	# create models depending on mode
	optm=LRAOpt(nt,nr,nra,T)

	# inversion variables allocation
	gx=X(length(optm.cal.g), T)
	sx=X(length(optm.cal.s), T)

	err=DataFrame(g=[], s=[], d=[])

	# dummy
	opG=LinearMap{T}(x->0.0, y->0.0,1,1, ismutating=true)
	opS=LinearMap{T}(x->0.0, y->0.0,1,1, ismutating=true)

	pa=LRA(om,optm,sx,gx,opG,opS,err)

	# update operators
	pa.opG=create_operator(pa, G())
	pa.opS=create_operator(pa, S())

	!(gobs===nothing) && copyto!(pa.optm.obs.g, gobs)
	!(sobs===nothing) && copyto!(pa.optm.obs.s, sobs)
	if(dobs===nothing)
		mul!(pa.optm.obs.d,pa.optm.obs.s,pa.optm.obs.g)
	else
		copyto!(pa.optm.obs.d, dobs)
	end

	initialize!(pa)

	return pa

end

function update_prepare!(pa::LRA, ::S)
end
function update_prepare!(pa::LRA, ::G)
end
function update_finalize!(pa::LRA, ::S)
end
function update_finalize!(pa::LRA, ::G)
end

function bd!(pa::LRA, io=stdout; tol=1e-6)

	if(io===nothing)
		logfilename=joinpath(pwd(),string("XBD",Dates.now(),".log"))
		io=open(logfilename, "a+")
	end

	update_all!(pa, io, max_reroundtrips=1, max_roundtrips=100000, roundtrip_tol=tol)

	err!(pa)
end

function F!(y, x, pa::LRA, ::S)
	for i in eachindex(x)
		pa.optm.cal.g[i]=x[i]
	end
	mul!(pa.optm.cal.d, pa.optm.cal.s, pa.optm.cal.g)
	copyto!(y,pa.optm.cal.d)
	return nothing
end

function F!(y, x, pa::LRA, ::G)
	for i in eachindex(x)
		pa.optm.cal.s[i]=x[i]
	end
	mul!(pa.optm.cal.d, pa.optm.cal.s, pa.optm.cal.g)
	copyto!(y,pa.optm.cal.d)
	return nothing
end

function F!(pa::LRA,x::AbstractVector, attrib::S)
	compute=(x!=pa.sx.last_x)
	if(compute)
		x_to_model!(x, pa, attrib) #
		copyto!(pa.sx.last_x, x)
		mul!(pa.optm.cal.d, pa.optm.cal.s, pa.optm.cal.g)
		return pa
	end
end

function F!(pa::LRA,x::AbstractVector, attrib::G)
	compute=(x!=pa.gx.last_x)
	if(compute)
		x_to_model!(x, pa, G())
		copyto!(pa.gx.last_x, x)
		mul!(pa.optm.cal.d, pa.optm.cal.s, pa.optm.cal.g)
		return pa
	end
end

function Fadj!(y, x, pa::LRA, ::S)
	copyto!(pa.optm.ddcal,x)
	mul!(pa.optm.cal.g,transpose(pa.optm.cal.s),pa.optm.ddcal)
	copyto!(y,pa.optm.cal.g)
	return nothing
end

function Fadj!(y, x, pa::LRA, ::G)
	copyto!(pa.optm.ddcal,x)
	mul!(pa.optm.cal.s,pa.optm.ddcal,transpose(pa.optm.cal.g))
	copyto!(y,pa.optm.cal.s)
	return nothing
end

function x_to_model!(x, pa::LRA, ::G)
	for i in eachindex(pa.optm.cal.g)
		pa.optm.cal.g[i]=x[i]*pa.gx.preconI[i]
	end
	return nothing
end
function x_to_model!(x, pa::LRA, ::S)
	for i in eachindex(pa.optm.cal.s)
		pa.optm.cal.s[i]=x[i]*pa.gx.preconI[i]
	end
	return nothing
end

function model_to_x!(x, pa::LRA, ::S)
	for i in eachindex(x)
		x[i]=pa.optm.cal.s[i]*pa.sx.precon[i]
	end
	return nothing
end
function model_to_x!(x, pa::LRA, ::G)
	for i in eachindex(x)
		x[i]=pa.optm.cal.g[i]*pa.gx.precon[i]
	end
	return nothing
end

function initialize!(pa::LRA)
	for i in eachindex(pa.optm.cal.g)
		x=(pa.gx.precon[i]≠0.0) ? randn() : 0.0
		pa.optm.cal.g[i]=x
	end

	for i in eachindex(pa.optm.cal.s)
		x=(pa.sx.precon[i]≠0.0) ? randn() : 0.0
		pa.optm.cal.s[i]=x
	end
end


text:
# Code derived from metabolicEP (https://github.com/anna-pa-m/Metabolic-EP)

# new Q_n tilted first moment
function newav(s, μ, av, va, siteflagave, siteflagvar, lb, ub, minvar, maxvar)
    sqrts = sqrt(s)
    xinf = (lb - μ) / sqrts
    xsup = (ub - μ) / sqrts
    scra1 , scra12 = compute_mom5d(xinf, xsup)
    avnew  = siteflagave ? μ + scra1 * sqrts : av # if is fixed use av
    varnew = siteflagvar ? max(minvar,s * (1.0 + scra12)) : va
    isnan(avnew) || isnan(varnew) && println("avnew = $avnew varnew = $varnew")
    return avnew, varnew
end


text:
using PtFEM
using Test

ProjDir = @__DIR__

#=
Compare formulas at:
http://www.awc.org/pdf/codes-standards/publications/design-aids/AWC-DA6-BeamFormulas-0710.pdf
=#

data = Dict(
  # Frame(nels, nn, ndim, nst, nip, finite_element(nod, nodof))
  :struc_el => Frame(20, 21, 3, 1, 1, Line(2, 3)),
  :properties => [1.0e3 1.0e3 1.0e3 3.0e5;],
  :x_coords => range(0, stop=4, length=21),
  :y_coords => zeros(21),
  :z_coords => zeros(21),
  :g_num => [
    collect(1:20)';
    collect(2:21)'],
  :support => [
    (1, [0 0 0 0 0 0]),
    (21, [0 0 0 0 0 0]),
    ],
  :loaded_nodes => [
      ( 1, [0.0 -1.0 0.0 0.0 0.0 0.0]),
      ( 2, [0.0 -1.0 0.0 0.0 0.0 0.0]),
      ( 3, [0.0 -1.0 0.0 0.0 0.0 0.0]),
      ( 4, [0.0 -1.0 0.0 0.0 0.0 0.0]),
      ( 5, [0.0 -1.0 0.0 0.0 0.0 0.0]),
      ( 6, [0.0 -1.0 0.0 0.0 0.0 0.0]),
      ( 7, [0.0 -1.0 0.0 0.0 0.0 0.0]),
      ( 8, [0.0 -1.0 0.0 0.0 0.0 0.0]),
      ( 9, [0.0 -1.0 0.0 0.0 0.0 0.0]),
      (10, [0.0 -1.0 0.0 0.0 0.0 0.0]),
      (11, [0.0 -1.0 0.0 0.0 0.0 0.0]),
      (12, [0.0 -1.0 0.0 0.0 0.0 0.0]),
      (13, [0.0 -1.0 0.0 0.0 0.0 0.0]),
      (14, [0.0 -1.0 0.0 0.0 0.0 0.0]),
      (15, [0.0 -1.0 0.0 0.0 0.0 0.0]),
      (16, [0.0 -1.0 0.0 0.0 0.0 0.0]),
      (17, [0.0 -1.0 0.0 0.0 0.0 0.0]),
      (18, [0.0 -1.0 0.0 0.0 0.0 0.0]),
      (19, [0.0 -1.0 0.0 0.0 0.0 0.0]),
      (20, [0.0 -1.0 0.0 0.0 0.0 0.0]),
      (21, [0.0 -1.0 0.0 0.0 0.0 0.0])
    ]
)

data |> display
println()

m, dis_df, fm_df = p44(data)

println("Displacements:")
m.displacements |> display
println()

println("Actions:")
m.actions |> display
println()

println("y displacements:")
m.displacements[2,:] |> display
println()

println("y moment actions:")
m.actions[12,:] |> display
println()

if VERSION.minor == 6
  using Plots
  gr(size=(400,600))

  p = Vector{Plots.Plot{Plots.GRBackend}}(undef, 3)
  p[1] = plot(m.x_coords, m.displacements[2,:], ylim=(-0.004, 0.001), lab="Displacement",
   xlabel="x [m]", ylabel="deflection [m]", color=:red)
  p[2] = plot(m.actions[2,:], lab="Shear force", ylim=(-10, 10), xlabel="element",
    ylabel="shear force [N]", palette=:greens,fill=(0,:auto),α=0.6)
  p[3] = plot(m.actions[12,:], lab="Moment", ylim=(-7, 7), xlabel="element",
    ylabel="moment [Nm]", palette=:grays,fill=(0,:auto),α=0.6)

  plot(p..., layout=(3, 1))
  savefig(ProjDir*"/figure-23.png")
  #=
  plot!()
  gui()
  =#
end

if VERSION.minor > 5

  # See figure 24 in above reference (Δmax):
  @eval @test m.displacements[2,11] ≈ -10000 * 4^3 / (192 * 1.0e6) atol=10.0*eps()

  # See figure 24 in above reference (Mmax):
  @eval @test m.actions[12,10] ≈ (10000 * 4 / 8) atol=10.0*eps()
end

max_stars_repo_path: examples/Nodal Formulas/figure-23.jl

max_issues_repo_path: examples/Nodal Formulas/figure-23.jl

max_forks_repo_path: examples/Nodal Formulas/figure-23.jl


text:
module TestManual

using Test
using Documenter
using Trixi2Img

DocMeta.setdocmeta!(Trixi2Img, :DocTestSetup, :(using Trixi2Img); recursive=true)
doctest(Trixi2Img, manual=false)

end # module


text:
"""
A module that makes expressions like `-π` behave like Irrational, rather
than Float64.

Generates IrrationalExpr objects when arithmetic operations are applied to
objects of type Irrational, and these operations are carried out in the
destination type when type conversion occurs.

For example, `BigFloat(π) + BigFloat(-π)` will evaluate to BigFloat(0),
rather than to approximately 1.2e-16.
"""

module IrrationalExpressions

import Base: +, -, *, /, convert, promote_rule, show

immutable IrrationalExpr{op, N} <: Real
  args::NTuple{N,Union{IrrationalExpr, Irrational, Rational, Integer}}
end

@generated convert{T<:AbstractFloat,op,N}(::Type{T}, x::IrrationalExpr{op,N}) =
  Expr(:call, op, [Expr(:call, :convert, T, :( x.args[$i] )) for i=1:N]...)

promote_rule{T1<:AbstractFloat, T2<:IrrationalExpr}(::Type{T1}, ::Type{T2}) = T1
promote_rule{T2<:IrrationalExpr}(::Type{BigFloat}, ::Type{T2}) = BigFloat

## Unary operators
(+)(x::IrrationalExpr) = x
(-)(x::IrrationalExpr) = IrrationalExpr{:(-),1}((x,))
(-)(x::Irrational) = IrrationalExpr{:(-),1}((x,))

## Binary operators
ops = (:(+), :(-), :(*), :(/))
types = (IrrationalExpr, Irrational, Rational, Integer)
for op in ops, i in eachindex(types), j in eachindex(types)
  if i<=2 || j<=2
    @eval $op(x::$(types[i]), y::$(types[j])) = IrrationalExpr{Symbol($op),2}((x,y))
  end
end

# We define expr() as a conveninent middle step to get to strings,
# but this also allows eval(expr(x)) as a roundabout way to get x.
expr(x) = x
expr{sym}(::Irrational{sym}) = sym
expr{op,N}(x::IrrationalExpr{op,N}) = Expr(:call, op, map(expr,x.args)...)

show(io::IO, x::IrrationalExpr) = print(io, string(expr(x)), " = ", string(Float64(x))[1:end-2], "...")

end # module


text:
module PhysicsFormulae

using LinearAlgebra
include("./Initial.jl")
using .Initial

const G = 6.67408e-11

# Calcualte the difference in velocities of two nodes due to gravitational force
function acceleration(node1, node2)
    # First we calculate the distance
    displacement = node1.position - node2.position
    distance = sum(square.(displacement))

    dv = normalize(displacement)
    dv = ((G * node2.mass) / square(distance)) .* dv

    return dv
end

# Calcualte the velocities of two nodes after collusion
function bounce(node1, node2)
    m₁, m₂ = node1.mass, node2.mass
    v₁, v₂ = node1.velocity, node2.velocity

    vf₁ = v₁ .* ((m₁ - m₂) / (m₂ + m₁))
    vf₁ += v₂ .* ((2 * m₂) / (m₂ + m₁))

    vf₂ = v₁ .* ((2 * m₁) / (m₂ + m₁))
    vf₂ += v₂ .* ((m₂ - m₁) / (m₂ + m₁))

    return vf₁, vf₂
end

# Calculate the combined velocity of two nodes after collusion
function combinedVelocity(node1, node2)
    m₁, m₂ = node1.mass, node2.mass
    v₁, v₂ = node1.velocity, node2.velocity

    vf = (m₁ .* v₁ + m₂ .* v₂) / (m₁ + m₂)

    return vf
end

function collideIfColliding(node1, node2)
    p₁, p₂ = node1.position, node2.position
    r₁, r₂ = node1.radius, node2.radius

    minDistanceForCollusion = r₁ + r₂
    distanceBetweenNodes = sum(square.(p₁ - p₂))

    if distanceBetweenNodes <= minDistanceForCollusion
        if onCollusion == "bounce"
            bounce(node1, node2)
        else
            combinedVelocity(node1, node2)
        end
    end
end

square(n) = n^2

export acceleration, bounce, combinedVelocity, collideIfColliding
end

max_stars_repo_name: PseudoDevelopers/Physics-Simulations

max_issues_repo_name: PseudoDevelopers/Physics-Simulations

max_forks_repo_name: PseudoDevelopers/Physics-Simulations


text:
"""
    ManyBodyBasis(b, occupations)

Basis for a many body system.

The basis has to know the associated one-body basis `b` and which occupation states
should be included. The occupations_hash is used to speed up checking if two
many-body bases are equal.
"""
struct ManyBodyBasis{S,B,H,UT} <: Basis
    shape::S
    onebodybasis::B
    occupations::Vector{S}
    occupations_hash::UT

    function ManyBodyBasis{S,B,H}(onebodybasis::B, occupations::Vector{S}) where {S,B,H}
        @assert isa(H, UInt)
        new{S,B,H,typeof(H)}([length(occupations)], onebodybasis, occupations, hash(hash.(occupations)))
    end
end
ManyBodyBasis(onebodybasis::B, occupations::Vector{S}) where {B,S} = ManyBodyBasis{S,B,hash(hash.(occupations))}(onebodybasis,occupations)

"""
    fermionstates(Nmodes, Nparticles)
    fermionstates(b, Nparticles)

Generate all fermionic occupation states for N-particles in M-modes.
`Nparticles` can be a vector to define a Hilbert space with variable
particle number.
"""
fermionstates(Nmodes::T, Nparticles::T) where T = _distribute_fermions(Nparticles, Nmodes, 1, zeros(Int, Nmodes), Vector{Int}[])
fermionstates(Nmodes::T, Nparticles::Vector{T}) where T = vcat([fermionstates(Nmodes, N) for N in Nparticles]...)
fermionstates(onebodybasis::Basis, Nparticles) = fermionstates(length(onebodybasis), Nparticles)

"""
    bosonstates(Nmodes, Nparticles)
    bosonstates(b, Nparticles)

Generate all bosonic occupation states for N-particles in M-modes.
`Nparticles` can be a vector to define a Hilbert space with variable
particle number.
"""
bosonstates(Nmodes::T, Nparticles::T) where T = _distribute_bosons(Nparticles, Nmodes, 1, zeros(Int, Nmodes), Vector{Int}[])
bosonstates(Nmodes::T, Nparticles::Vector{T}) where T = vcat([bosonstates(Nmodes, N) for N in Nparticles]...)
bosonstates(onebodybasis::Basis, Nparticles) = bosonstates(length(onebodybasis), Nparticles)

==(b1::ManyBodyBasis, b2::ManyBodyBasis) = b1.occupations_hash==b2.occupations_hash && b1.onebodybasis==b2.onebodybasis

"""
    basisstate(b::ManyBodyBasis, occupation::Vector)

Return a ket state where the system is in the state specified by the given
occupation numbers.
"""
function basisstate(basis::ManyBodyBasis, occupation::Vector)
    index = findfirst(isequal(occupation), basis.occupations)
    if isa(index, Nothing)
        throw(ArgumentError("Occupation not included in many-body basis."))
    end
    basisstate(basis, index)
end

function isnonzero(occ1, occ2, index)
    for i=1:length(occ1)
        if i == index
            if occ1[i] != occ2[i] + 1
                return false
            end
        else
            if occ1[i] != occ2[i]
                return false
            end
        end
    end
    true
end

"""
    create(b::ManyBodyBasis, index)

Creation operator for the i-th mode of the many-body basis `b`.
"""
function create(b::ManyBodyBasis, index)
    result = SparseOperator(b)
    # <{m}_i| at |{m}_j>
    for i=1:length(b)
        occ_i = b.occupations[i]
        if occ_i[index] == 0
            continue
        end
        for j=1:length(b)
            if isnonzero(occ_i, b.occupations[j], index)
                result.data[i, j] = sqrt(occ_i[index])
            end
        end
    end
    result
end

"""
    destroy(b::ManyBodyBasis, index)

Annihilation operator for the i-th mode of the many-body basis `b`.
"""
function destroy(b::ManyBodyBasis, index)
    result = SparseOperator(b)
    # <{m}_j| a |{m}_i>
    for i=1:length(b)
        occ_i = b.occupations[i]
        if occ_i[index] == 0
            continue
        end
        for j=1:length(b)
            if isnonzero(occ_i, b.occupations[j], index)
                result.data[j, i] = sqrt(occ_i[index])
            end
        end
    end
    result
end

"""
    number(b::ManyBodyBasis, index)

Particle number operator for the i-th mode of the many-body basis `b`.
"""
function number(b::ManyBodyBasis, index)
    result = SparseOperator(b)
    for i=1:length(b)
        result.data[i, i] = b.occupations[i][index]
    end
    result
end

"""
    number(b::ManyBodyBasis)

Total particle number operator.
"""
function number(b::ManyBodyBasis)
    result = SparseOperator(b)
    for i=1:length(b)
        result.data[i, i] = sum(b.occupations[i])
    end
    result
end

function isnonzero(occ1, occ2, index1, index2)
    for i=1:length(occ1)
        if i == index1 && i == index2
            if occ1[i] != occ2[i]
                return false
            end
        elseif i == index1
            if occ1[i] != occ2[i] + 1
                return false
            end
        elseif i == index2
            if occ1[i] != occ2[i] - 1
                return false
            end
        else
            if occ1[i] != occ2[i]
                return false
            end
        end
    end
    true
end

"""
    transition(b::ManyBodyBasis, to, from)

Operator ``|\\mathrm{to}⟩⟨\\mathrm{from}|`` transferring particles between modes.
"""
function transition(b::ManyBodyBasis, to, from)
    result = SparseOperator(b)
    # <{m}_j| at_to a_from |{m}_i>
    for i=1:length(b)
        occ_i = b.occupations[i]
        if occ_i[from] == 0
            continue
        end
        for j=1:length(b)
            occ_j = b.occupations[j]
            if isnonzero(occ_j, occ_i, to, from)
                result.data[j, i] = sqrt(occ_i[from])*sqrt(occ_j[to])
            end
        end
    end
    result
end

# Calculate many-Body operator from one-body operator
"""
    manybodyoperator(b::ManyBodyBasis, op)

Create the many-body operator from the given one-body operator `op`.

The given operator can either be a one-body operator or a
two-body interaction. Higher order interactions are at the
moment not implemented.

The mathematical formalism for the one-body case is described by

```math
X = \\sum_{ij} a_i^† a_j ⟨u_i| x | u_j⟩
```

and for the interaction case by

```math
X = \\sum_{ijkl} a_i^† a_j^† a_k a_l ⟨u_i|⟨u_j| x |u_k⟩|u_l⟩
```

where ``X`` is the N-particle operator, ``x`` is the one-body operator and
``|u⟩`` are the one-body states associated to the
different modes of the N-particle basis.
"""
function manybodyoperator(basis::ManyBodyBasis, op)
    @assert op.basis_l == op.basis_r
    if op.basis_l == basis.onebodybasis
        result =  manybodyoperator_1(basis, op)
    elseif op.basis_l == basis.onebodybasis ⊗ basis.onebodybasis
        result = manybodyoperator_2(basis, op)
    else
        throw(ArgumentError("The basis of the given operator has to either be equal to b or b ⊗ b where b is the 1st quantization basis associated to the nparticle basis."))
    end
    result
end

function manybodyoperator_1(basis::ManyBodyBasis, op::Operator)
    N = length(basis)
    S = length(basis.onebodybasis)
    result = DenseOperator(basis)
    @inbounds for n=1:N, m=1:N
        for j=1:S, i=1:S
            C = coefficient(basis.occupations[m], basis.occupations[n], [i], [j])
            if C != 0.
                result.data[m,n] += C*op.data[i,j]
            end
        end
    end
    return result
end
manybodyoperator_1(basis::ManyBodyBasis, op::AdjointOperator) = dagger(manybodyoperator_1(basis,dagger(op)))

function manybodyoperator_1(basis::ManyBodyBasis, op::SparseOpPureType)
    N = length(basis)
    S = length(basis.onebodybasis)
    result = SparseOperator(basis)
    M = op.data
    @inbounds for colindex = 1:M.n
        for i=M.colptr[colindex]:M.colptr[colindex+1]-1
            row = M.rowval[i]
            value = M.nzval[i]
            for m=1:N, n=1:N
                C = coefficient(basis.occupations[m], basis.occupations[n], [row], [colindex])
                if C != 0.
                    result.data[m, n] += C*value
                end
            end
        end
    end
    return result
end

function manybodyoperator_2(basis::ManyBodyBasis, op::Operator)
    N = length(basis)
    S = length(basis.onebodybasis)
    @assert S^2 == length(op.basis_l)
    @assert S^2 == length(op.basis_r)
    result = DenseOperator(basis)
    op_data = reshape(op.data, S, S, S, S)
    occupations = basis.occupations
    @inbounds for m=1:N, n=1:N
        for l=1:S, k=1:S, j=1:S, i=1:S
            C = coefficient(occupations[m], occupations[n], [i, j], [k, l])
            result.data[m,n] += C*op_data[i, j, k, l]
        end
    end
    return result
end

function manybodyoperator_2(basis::ManyBodyBasis, op::SparseOpType)
    N = length(basis)
    S = length(basis.onebodybasis)
    result = SparseOperator(basis)
    occupations = basis.occupations
    rows = rowvals(op.data)
    values = nonzeros(op.data)
    @inbounds for column=1:S^2, j in nzrange(op.data, column)
        row = rows[j]
        value = values[j]
        for m=1:N, n=1:N
            # println("row:", row, " column:"column, ind_left)
            index = Tuple(CartesianIndices((S, S, S, S))[(column-1)*S^2 + row])
            C = coefficient(occupations[m], occupations[n], index[1:2], index[3:4])
            if C!=0.
                result.data[m,n] += C*value
            end
        end
    end
    return result
end

# Calculate expectation value of one-body operator
"""
    onebodyexpect(op, state)

Expectation value of the one-body operator `op` in respect to the many-body `state`.
"""
function onebodyexpect(op::AbstractOperator, state::Ket)
    @assert isa(state.basis, ManyBodyBasis)
    @assert op.basis_l == op.basis_r
    if state.basis.onebodybasis == op.basis_l
        result = onebodyexpect_1(op, state)
    # Not yet implemented:
    # elseif state.basis.basis ⊗ state.basis.basis == op.basis_l
    #     result = onebodyexpect_2(op, state)
    else
        throw(ArgumentError("The basis of the given operator has to either be equal to b or b ⊗ b where b is the 1st quantization basis associated to the nparticle basis of the state."))
    end
    result
end

function onebodyexpect(op::AbstractOperator, state::AbstractOperator)
    @assert op.basis_l == op.basis_r
    @assert state.basis_l == state.basis_r
    @assert isa(state.basis_l, ManyBodyBasis)
    if state.basis_l.onebodybasis == op.basis_l
        result = onebodyexpect_1(op, state)
    # Not yet implemented
    # elseif state.basis.basis ⊗ state.basis.basis == op.basis_l
    #     result = onebodyexpect_2(op, state)
    else
        throw(ArgumentError("The basis of the given operator has to either be equal to b or b ⊗ b where b is the 1st quantization basis associated to the nparticle basis of the state."))
    end
    result
end
onebodyexpect(op::AbstractOperator, states::Vector) = [onebodyexpect(op, state) for state=states]

function onebodyexpect_1(op::Operator, state::Ket)
    N = length(state.basis)
    S = length(state.basis.onebodybasis)
    result = complex(0.)
    occupations = state.basis.occupations
    for m=1:N, n=1:N
        value = conj(state.data[m])*state.data[n]
        for i=1:S, j=1:S
            C = coefficient(occupations[m], occupations[n], [i], [j])
            if C != 0.
                result += C*op.data[i,j]*value
            end
        end
    end
    result
end

function onebodyexpect_1(op::Operator, state::Operator)
    N = length(state.basis_l)
    S = length(state.basis_l.onebodybasis)
    result = complex(zero(promote_type(eltype(op),eltype(state))))
    occupations = state.basis_l.occupations
    @inbounds for s=1:N, t=1:N
        value = state.data[t,s]
        for i=1:S, j=1:S
            C = coefficient(occupations[s], occupations[t], [i], [j])
            if !iszero(C)
                result += C*op.data[i,j]*value
            end
        end
    end
    result
end

function onebodyexpect_1(op::SparseOpPureType, state::Ket)
    N = length(state.basis)
    S = length(state.basis.onebodybasis)
    result = complex(0.)
    occupations = state.basis.occupations
    M = op.data
    @inbounds for colindex = 1:M.n
        for i=M.colptr[colindex]:M.colptr[colindex+1]-1
            row = M.rowval[i]
            value = M.nzval[i]
            for m=1:N, n=1:N
                C = coefficient(occupations[m], occupations[n], [row], [colindex])
                if C != 0.
                    result += C*value*conj(state.data[m])*state.data[n]
                end
            end
        end
    end
    result
end

function onebodyexpect_1(op::SparseOpPureType, state::Operator)
    N = length(state.basis_l)
    S = length(state.basis_l.onebodybasis)
    result = complex(0.)
    occupations = state.basis_l.occupations
    M = op.data
    @inbounds for colindex = 1:M.n
        for i=M.colptr[colindex]:M.colptr[colindex+1]-1
            row = M.rowval[i]
            value = M.nzval[i]
            for s=1:N, t=1:N
                C = coefficient(occupations[s], occupations[t], [row], [colindex])
                if C != 0.
                    result += C*value*state.data[t,s]
                end
            end
        end
    end
    result
end

"""
Calculate the matrix element <{m}|at_1...at_n a_1...a_n|{n}>.
"""
function coefficient(occ_m, occ_n, at_indices, a_indices)
    occ_m = copy(occ_m)
    occ_n = copy(occ_n)
    C = 1.
    for i=at_indices
        if occ_m[i] == 0
            return 0.
        end
        C *= sqrt(occ_m[i])
        occ_m[i] -= 1
    end
    for i=a_indices
        if occ_n[i] == 0
            return 0.
        end
        C *= sqrt(occ_n[i])
        occ_n[i] -= 1
    end
    if occ_m == occ_n
        return C
    else
        return 0.
    end
end

function _distribute_bosons(Nparticles, Nmodes, index, occupations, results)
    if index==Nmodes
        occupations[index] = Nparticles
        push!(results, copy(occupations))
    else
        for n=Nparticles:-1:0
            occupations[index] = n
            _distribute_bosons(Nparticles-n, Nmodes, index+1, occupations, results)
        end
    end
    return results
end

function _distribute_fermions(Nparticles, Nmodes, index, occupations, results)
    if (Nmodes-index)+1<Nparticles
        return results
    end
    if index==Nmodes
        occupations[index] = Nparticles
        push!(results, copy(occupations))
    else
        for n=min(1,Nparticles):-1:0
            occupations[index] = n
            _distribute_fermions(Nparticles-n, Nmodes, index+1, occupations, results)
        end
    end
    return results
end


text:
@testset "Interface" begin
    @test GreeDS().alg == PCA()
    @test GreeDS().threshold == 0
    @test GreeDS(7).alg == PCA(7)
end

@testset "Decomposition" begin
    data = 4 .* randn(rng, 30, 101, 101) .+ 10
    angles = sort!(90randn(rng, 30)) |> normalize_par_angles

    # get sizes correct for ncomps
    for N in [1, 3, 5]
        A, w = @inferred decompose(GreeDS(PCA(ncomps=N)), data, angles)
        @test size(A) == (N, 101 * 101)
        @test size(w) == (30, N)
    end
    @test_throws ErrorException decompose(GreeDS(PCA(40)), data, angles)
    A, w = decompose(PCA(30; pratio=0.5), data, angles)
    @test size(A, 1) < 30
    @test size(w, 2) < 30

    # default is to use whole cube
    S = reconstruct(GreeDS(PCA()), data, angles)
    @test size(S) == (30, 101, 101)
end

@testset "ADI Trivial" begin
    cube = ones(10, 100, 100)
    angles = zeros(10)

    reduced_5 = GreeDS(PCA(5))(cube, angles)
    reduced_10 = GreeDS(PCA(10))(cube, angles)

    @test size(reduced_5) == size(reduced_10) == (100, 100)
    @test all(x -> isapprox(x, 0, atol = 1e-9), reduced_5)
    @test all(x -> isapprox(x, 0, atol = 1e-9), reduced_10)
end

@testset "RDI Trivial" begin
    data = randn(rng, 30, 101, 101)
    angles = sort!(90rand(rng, 30)) |> normalize_par_angles

    # not supported for GreeDS alg
    S = data .- reconstruct(GreeDS(PCA(1)), data, angles, zeros(30, 101, 101))
    @test S ≈ data rtol=2e-1
end

@testset "Decomposition - TPCA" begin
    data = 4 .* randn(rng, 30, 101, 101) .+ 10
    angles = sort!(90randn(rng, 30)) |> normalize_par_angles

    # get sizes correct for ncomps
    for N in [1, 3, 5]
        # A, w = @inferred decompose(GreeDS(TPCA(N)), data, angles)
        A, w = decompose(GreeDS(TPCA(N)), data, angles)
        @test size(A) == (N, 101 * 101)
        @test size(w) == (30, N)
    end
    @test_throws ErrorException decompose(GreeDS(TPCA(40)), data, angles)

    # default is to use whole cube
    S = reconstruct(GreeDS(TPCA()), data, angles)
    @test size(S) == (30, 101, 101)
end

@testset "Decomposition - NMF" begin
    data = 4 .* randn(rng, 10, 101, 101) .+ 10
    angles = sort!(90randn(rng, 10)) |> normalize_par_angles

    # get sizes correct for ncomps
    N = 3
    A, w = decompose(GreeDS(NMF(N)), data, angles)
    @test size(A) == (N, 101 * 101)
    @test size(w) == (10, N)

    @test_throws ErrorException decompose(GreeDS(NMF(40)), data, angles)

    # default is to use whole cube
    S = reconstruct(GreeDS(NMF()), data, angles)
    @test size(S) == (10, 101, 101)
end


text:
using Distributed
export parallelize

parallelize(f, bd, t, n::Int) = parallelize(f, bd, t, [randominside(bd) for i in 1:n])
parallelize(f, bd, t, n::Int, ω) = parallelize(f, bd, t, [randominside(bd, ω) for i in 1:n])

"""
    parallelize(f, bd::Billiard, t, particles; partype = :threads)
Parallelize function `f` across the available particles. The parallelization type can
be `:threads` or `:pmap`, which use threads or a worker pool initialized with `addprocs`
_before_ `using DynamicalBilliards`.

`particles` can be:
* A `Vector` of particles.
* An integer `n` optionally followed by an angular velocity `ω`.
  This uses [`randominside`](@ref).

The functions usable here are:
* [`meancollisiontime`](@ref)
* [`escapetime`](@ref)
* [`lyapunovspectrum`](@ref) (returns only the maximal exponents)
* [`boundarymap`](@ref) (returns vector of vectors of 2-vectors _and_ `arcintervals`)
"""
function parallelize(f, bd::Billiard, t, particles::Vector{<:AbstractParticle};
    partype = :threads)
    if partype == :threads
        return threads_pl(f, bd, t, particles)
    elseif partype == :pmap
        return pmap_pl(f, bd, t, particles)
    end
end

function threads_pl(f, bd, t, particles)
    ret = _retinit(f, particles)
    Threads.@threads for i in 1:length(particles)
        @inbounds ret[i] = _getval(f, particles[i], bd, t)
    end
    return ret
end

function pmap_pl(f, bd, t, particles)
    g(p) = _getval(f, p, bd, t)
    ret = pmap(g, particles)
    return ret
end

_retinit(f, p::Vector{<:AbstractParticle{T}}) where {T} = zeros(T, length(p))
_retinit(::typeof(boundarymap), p::Vector{<:AbstractParticle{T}}) where {T} =
    Vector{Vector{SV{T}}}(undef, length(p))

_getval(f, p, bd, t) = f(p, bd, t)
_getval(f::typeof(lyapunovspectrum), p, bd, t) = @inbounds f(p, bd, t)[1]
_getval(f::typeof(lyapunovspectrum!), p, bd, t) = @inbounds f(p, bd, t)[1]

# Methods for boundary map are trickier because of the weird call signature
# and return signature
function threads_pl(f::typeof(boundarymap), bd, t, particles)
    intervals = arcintervals(bd)
    ret = _retinit(f, particles)
    Threads.@threads for i in 1:length(particles)
        @inbounds ret[i] = f(particles[i], bd, t, intervals)[1]
    end
    return ret, intervals
end
function pmap_pl(f::typeof(boundarymap), bd, t, particles)
    intervals = arcintervals(bd)
    g(p) = f(p, bd, t, intervals)
    ret = pmap(g, particles)
    return ret, intervals
end


text:
function is_installed(pkg)
    try
        return(Pkg.installed(pkg) === nothing ? false: true)
    catch
        return(false)
    end
end

if is_installed("NDSparseData")
    using NDSparseData

    source_ndsparsearray1 = NDSparse([fill("New York",3); fill("Boston",3)],
                            repmat(Date(2016,7,6):Date(2016,7,8), 2),
                            [91,89,91,95,83,76])

    q = @from i in source_ndsparsearray1 begin
        @where i.index[1]=="Boston"
        @select i.value
        @collect
    end

    @test isa(q, Array{Int,1})
    @test length(q)==3
    @test q==[95,83,76]

    source_ndsparsearray2 = NDSparse(Columns(city = [fill("New York",3); fill("Boston",3)],
                            date = repmat(Date(2016,7,6):Date(2016,7,8), 2)),
                            [91,89,91,95,83,76])

    q = @from i in source_ndsparsearray2 begin
        @where i.index.city=="New York"
        @select i.value
        @collect
    end

    @test isa(q, Array{Int,1})
    @test length(q)==3
    @test q==[91,89,91]
end


text:
## Exercise 8-5
## Encapsulate this code in a function named count, and generalize it so that it accepts the string and the letter as arguments.
println("Ans 1: ")

function count(word::AbstractString, to_count::Char)
    counter = 0
    for letter in word
        if letter == to_count
            counter = counter + 1
        end
    end
    counter
end

println(count("banana", 'a'))

## Then rewrite the function so that instead of traversing the string, it uses the three-parameter version of find from the previous section.
println("Ans 2: ")

function find(word, letter, start::Int = 1)
    word = word[start : end]
    index = firstindex(word)
    while index <= sizeof(word)
        if word[index] == letter
            return index + start - 1
        end
        index = nextind(word, index)
    end
    -1
end

function count_find(word::AbstractString, letter::Char)
    counter = 0
    last_ind = 1
    while true
        ind = find(word, letter, last_ind)
        if ind > 0
            counter += 1
            last_ind = ind + 1
        end

        if last_ind > length(word) || ind == -1
            break;
        end
    end
    counter
end

println(count_find("bananana", 'a'))

println("End.")


text:
using Myterial
using Documenter

DocMeta.setdocmeta!(Myterial, :DocTestSetup, :(using Myterial); recursive=true)

makedocs(;
    modules=[Myterial],
    authors="FedeClaudi <federicoclaudi@protonmail.com> and contributors",
    repo="https://github.com/fedeclaudi/Myterial.jl/blob/{commit}{path}#{line}",
    sitename="Myterial.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://fedeclaudi.github.io/Myterial.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/fedeclaudi/Myterial.jl",
    devbranch="main",
)