text:
int_rules_7_1_5 = @theory begin

    #= ::Subsection::Closed:: =#
    #= 7.1.5*u*(a+b*arcsinh(c*x))^n =#
    @apply_utils Antiderivative((~(a') + ~(b') * asinh(~(c') * ~x)) ^ ~(n') / (~(d') + ~(e') * ~x), ~x) => Subst(Antiderivative(((~a + ~b * ~x) ^ ~n * cosh(~x)) / (~c * ~d + ~e * sinh(~x)), ~x), ~x, asinh(~c * ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x) && IGtQ(~n, 0)
    @apply_utils Antiderivative((~(d') + ~(e') * ~x) ^ ~(m') * (~(a') + ~(b') * asinh(~(c') * ~x)) ^ ~(n'), ~x) => ((~d + ~e * ~x) ^ (~m + 1) * (~a + ~b * asinh(~c * ~x)) ^ ~n) / (~e * (~m + 1)) - ((~b * ~c * ~n) / (~e * (~m + 1))) * Antiderivative(((~d + ~e * ~x) ^ (~m + 1) * (~a + ~b * asinh(~c * ~x)) ^ (~n - 1)) / sqrt(1 + (~c) ^ 2 * (~x) ^ 2), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~m], ~x) && (IGtQ(~n, 0) && NeQ(~m, -1))
    @apply_utils Antiderivative((~d + ~(e') * ~x) ^ ~(m') * (~(a') + ~(b') * asinh(~(c') * ~x)) ^ ~n, ~x) => Antiderivative(ExpandIntegrand((~d + ~e * ~x) ^ ~m * (~a + ~b * asinh(~c * ~x)) ^ ~n, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x) && (IGtQ(~m, 0) && LtQ(~n, -1))
    @apply_utils Antiderivative((~(d') + ~(e') * ~x) ^ ~(m') * (~(a') + ~(b') * asinh(~(c') * ~x)) ^ ~(n'), ~x) => (1 / (~c) ^ (~m + 1)) * Subst(Antiderivative((~a + ~b * ~x) ^ ~n * cosh(~x) * (~c * ~d + ~e * sinh(~x)) ^ ~m, ~x), ~x, asinh(~c * ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~n], ~x) && IGtQ(~m, 0)
    @apply_utils Antiderivative(~Px * (~(a') + ~(b') * asinh(~(c') * ~x)), ~x) => With([u = IntHide(ExpandExpression(~Px, ~x), ~x)], Dist(~a + ~b * asinh(~c * ~x), u, ~x) - ~b * ~c * Antiderivative(SimplifyIntegrand(u / sqrt(1 + (~c) ^ 2 * (~x) ^ 2), ~x), ~x)) <-- FreeQ([~a, ~b, ~c], ~x) && PolynomialQ(~Px, ~x)
    #= Antiderivative((~Px)*((~a')+(~b')*asinh((~c')*(~x)))^(~n'),~x) :=  With([u=IntHide(Px,x)], Dist((a+b*asinh(c*x))^n,u,x) -  b*c*n*Antiderivative(SimplifyIntegrand(u*(a+b*asinh(c*x))^(n-1)/sqrt(1+c^2*x^2) ,x),x)) <-- FreeQ([a,b,c],x) && PolynomialQ(Px,x) && IGtQ(n,0) =#
    @apply_utils Antiderivative(~Px * (~(a') + ~(b') * asinh(~(c') * ~x)) ^ ~n, ~x) => Antiderivative(ExpandIntegrand(~Px * (~a + ~b * asinh(~c * ~x)) ^ ~n, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~n], ~x) && PolynomialQ(~Px, ~x)
    @apply_utils Antiderivative(~Px * (~(d') + ~(e') * ~x) ^ ~(m') * (~(a') + ~(b') * asinh(~(c') * ~x)), ~x) => With([u = IntHide(~Px * (~d + ~e * ~x) ^ ~m, ~x)], Dist(~a + ~b * asinh(~c * ~x), u, ~x) - ~b * ~c * Antiderivative(SimplifyIntegrand(u / sqrt(1 + (~c) ^ 2 * (~x) ^ 2), ~x), ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~m], ~x) && PolynomialQ(~Px, ~x)
    @apply_utils Antiderivative((~(f') + ~(g') * ~x) ^ ~(p') * (~d + ~(e') * ~x) ^ ~m * (~(a') + ~(b') * asinh(~(c') * ~x)) ^ ~n, ~x) => With([u = IntHide((~f + ~g * ~x) ^ ~p * (~d + ~e * ~x) ^ ~m, ~x)], Dist((~a + ~b * asinh(~c * ~x)) ^ ~n, u, ~x) - ~b * ~c * ~n * Antiderivative(SimplifyIntegrand((u * (~a + ~b * asinh(~c * ~x)) ^ (~n - 1)) / sqrt(1 + (~c) ^ 2 * (~x) ^ 2), ~x), ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g], ~x) && (IGtQ(~n, 0) && (IGtQ(~p, 0) && (ILtQ(~m, 0) && LtQ(~m + ~p + 1, 0))))
    @apply_utils Antiderivative(((~(f') + ~(g') * ~x + ~(h') * (~x) ^ 2) ^ ~(p') * (~(a') + ~(b') * asinh(~(c') * ~x)) ^ ~n) / (~d + ~(e') * ~x) ^ 2, ~x) => With([u = IntHide((~f + ~g * ~x + ~h * (~x) ^ 2) ^ ~p / (~d + ~e * ~x) ^ 2, ~x)], Dist((~a + ~b * asinh(~c * ~x)) ^ ~n, u, ~x) - ~b * ~c * ~n * Antiderivative(SimplifyIntegrand((u * (~a + ~b * asinh(~c * ~x)) ^ (~n - 1)) / sqrt(1 + (~c) ^ 2 * (~x) ^ 2), ~x), ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~h], ~x) && (IGtQ(~n, 0) && (IGtQ(~p, 0) && EqQ(~e * ~g - 2 * ~d * ~h, 0)))
    @apply_utils Antiderivative(~Px * (~d + ~(e') * ~x) ^ ~(m') * (~(a') + ~(b') * asinh(~(c') * ~x)) ^ ~(n'), ~x) => Antiderivative(ExpandIntegrand(~Px * (~d + ~e * ~x) ^ ~m * (~a + ~b * asinh(~c * ~x)) ^ ~n, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x) && (PolynomialQ(~Px, ~x) && (IGtQ(~n, 0) && IntegerQ(~m)))
    @apply_utils Antiderivative((~f + ~(g') * ~x) ^ ~(m') * (~d + ~(e') * (~x) ^ 2) ^ ~p * (~(a') + ~(b') * asinh(~(c') * ~x)), ~x) => With([u = IntHide((~f + ~g * ~x) ^ ~m * (~d + ~e * (~x) ^ 2) ^ ~p, ~x)], Dist(~a + ~b * asinh(~c * ~x), u, ~x) - ~b * ~c * Antiderivative(Dist(1 / sqrt(1 + (~c) ^ 2 * (~x) ^ 2), u, ~x), ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g], ~x) && (EqQ(~e, (~c) ^ 2 * ~d) && (IGtQ(~m, 0) && (ILtQ(~p + 1 / 2, 0) && (GtQ(~d, 0) && (LtQ(~m, -2 * ~p - 1) || GtQ(~m, 3))))))
    @apply_utils Antiderivative((~f + ~(g') * ~x) ^ ~(m') * (~d + ~(e') * (~x) ^ 2) ^ ~p * (~(a') + ~(b') * asinh(~(c') * ~x)) ^ ~(n'), ~x) => Antiderivative(ExpandIntegrand((~d + ~e * (~x) ^ 2) ^ ~p * (~a + ~b * asinh(~c * ~x)) ^ ~n, (~f + ~g * ~x) ^ ~m, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g], ~x) && (EqQ(~e, (~c) ^ 2 * ~d) && (IGtQ(~m, 0) && (IntegerQ(~p + 1 / 2) && (GtQ(~d, 0) && (IGtQ(~n, 0) && (EqQ(~n, 1) && GtQ(~p, -1) || (GtQ(~p, 0) || (EqQ(~m, 1) || EqQ(~m, 2) && LtQ(~p, -2)))))))))
    @apply_utils Antiderivative((~(f') + ~(g') * ~x) ^ ~m * sqrt(~d + ~(e') * (~x) ^ 2) * (~(a') + ~(b') * asinh(~(c') * ~x)) ^ ~(n'), ~x) => ((~f + ~g * ~x) ^ ~m * (~d + ~e * (~x) ^ 2) * (~a + ~b * asinh(~c * ~x)) ^ (~n + 1)) / (~b * ~c * sqrt(~d) * (~n + 1)) - (1 / (~b * ~c * sqrt(~d) * (~n + 1))) * Antiderivative((~d * ~g * ~m + 2 * ~e * ~f * ~x + ~e * ~g * (~m + 2) * (~x) ^ 2) * (~f + ~g * ~x) ^ (~m - 1) * (~a + ~b * asinh(~c * ~x)) ^ (~n + 1), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g], ~x) && (EqQ(~e, (~c) ^ 2 * ~d) && (ILtQ(~m, 0) && (GtQ(~d, 0) && IGtQ(~n, 0))))
    @apply_utils Antiderivative((~f + ~(g') * ~x) ^ ~(m') * (~d + ~(e') * (~x) ^ 2) ^ ~p * (~(a') + ~(b') * asinh(~(c') * ~x)) ^ ~(n'), ~x) => Antiderivative(ExpandIntegrand(sqrt(~d + ~e * (~x) ^ 2) * (~a + ~b * asinh(~c * ~x)) ^ ~n, (~f + ~g * ~x) ^ ~m * (~d + ~e * (~x) ^ 2) ^ (~p - 1 / 2), ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g], ~x) && (EqQ(~e, (~c) ^ 2 * ~d) && (IntegerQ(~m) && (IGtQ(~p + 1 / 2, 0) && (GtQ(~d, 0) && IGtQ(~n, 0)))))
    @apply_utils Antiderivative((~f + ~(g') * ~x) ^ ~(m') * (~d + ~(e') * (~x) ^ 2) ^ ~p * (~(a') + ~(b') * asinh(~(c') * ~x)) ^ ~(n'), ~x) => ((~f + ~g * ~x) ^ ~m * (~d + ~e * (~x) ^ 2) ^ (~p + 1 / 2) * (~a + ~b * asinh(~c * ~x)) ^ (~n + 1)) / (~b * ~c * sqrt(~d) * (~n + 1)) - (1 / (~b * ~c * sqrt(~d) * (~n + 1))) * Antiderivative(ExpandIntegrand((~f + ~g * ~x) ^ (~m - 1) * (~a + ~b * asinh(~c * ~x)) ^ (~n + 1), (~d * ~g * ~m + ~e * ~f * (2 * ~p + 1) * ~x + ~e * ~g * (~m + 2 * ~p + 1) * (~x) ^ 2) * (~d + ~e * (~x) ^ 2) ^ (~p - 1 / 2), ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g], ~x) && (EqQ(~e, (~c) ^ 2 * ~d) && (ILtQ(~m, 0) && (IGtQ(~p - 1 / 2, 0) && (GtQ(~d, 0) && IGtQ(~n, 0)))))
    @apply_utils Antiderivative(((~f + ~(g') * ~x) ^ ~(m') * (~(a') + ~(b') * asinh(~(c') * ~x)) ^ ~n) / sqrt(~d + ~(e') * (~x) ^ 2), ~x) => ((~f + ~g * ~x) ^ ~m * (~a + ~b * asinh(~c * ~x)) ^ (~n + 1)) / (~b * ~c * sqrt(~d) * (~n + 1)) - ((~g * ~m) / (~b * ~c * sqrt(~d) * (~n + 1))) * Antiderivative((~f + ~g * ~x) ^ (~m - 1) * (~a + ~b * asinh(~c * ~x)) ^ (~n + 1), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g], ~x) && (EqQ(~e, (~c) ^ 2 * ~d) && (IGtQ(~m, 0) && (GtQ(~d, 0) && LtQ(~n, -1))))
    @apply_utils Antiderivative(((~f + ~(g') * ~x) ^ ~(m') * (~(a') + ~(b') * asinh(~(c') * ~x)) ^ ~(n')) / sqrt(~d + ~(e') * (~x) ^ 2), ~x) => (1 / ((~c) ^ (~m + 1) * sqrt(~d))) * Subst(Antiderivative((~a + ~b * ~x) ^ ~n * (~c * ~f + ~g * sinh(~x)) ^ ~m, ~x), ~x, asinh(~c * ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~n], ~x) && (EqQ(~e, (~c) ^ 2 * ~d) && (IntegerQ(~m) && (GtQ(~d, 0) && (GtQ(~m, 0) || IGtQ(~n, 0)))))
    @apply_utils Antiderivative((~f + ~(g') * ~x) ^ ~(m') * (~d + ~(e') * (~x) ^ 2) ^ ~p * (~(a') + ~(b') * asinh(~(c') * ~x)) ^ ~(n'), ~x) => Antiderivative(ExpandIntegrand((~a + ~b * asinh(~c * ~x)) ^ ~n / sqrt(~d + ~e * (~x) ^ 2), (~f + ~g * ~x) ^ ~m * (~d + ~e * (~x) ^ 2) ^ (~p + 1 / 2), ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g], ~x) && (EqQ(~e, (~c) ^ 2 * ~d) && (IntegerQ(~m) && (ILtQ(~p + 1 / 2, 0) && (GtQ(~d, 0) && IGtQ(~n, 0)))))
    @apply_utils Antiderivative((~f + ~(g') * ~x) ^ ~(m') * (~d + ~(e') * (~x) ^ 2) ^ ~p * (~(a') + ~(b') * asinh(~(c') * ~x)) ^ ~(n'), ~x) => Simp((~d + ~e * (~x) ^ 2) ^ ~p / (1 + (~c) ^ 2 * (~x) ^ 2) ^ ~p) * Antiderivative((~f + ~g * ~x) ^ ~m * (1 + (~c) ^ 2 * (~x) ^ 2) ^ ~p * (~a + ~b * asinh(~c * ~x)) ^ ~n, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~n], ~x) && (EqQ(~e, (~c) ^ 2 * ~d) && (IntegerQ(~m) && (IntegerQ(~p - 1 / 2) && Not(GtQ(~d, 0)))))
    @apply_utils Antiderivative((log(~(h') * (~(f') + ~(g') * ~x) ^ ~(m')) * (~(a') + ~(b') * asinh(~(c') * ~x)) ^ ~(n')) / sqrt(~d + ~(e') * (~x) ^ 2), ~x) => (log(~h * (~f + ~g * ~x) ^ ~m) * (~a + ~b * asinh(~c * ~x)) ^ (~n + 1)) / (~b * ~c * sqrt(~d) * (~n + 1)) - ((~g * ~m) / (~b * ~c * sqrt(~d) * (~n + 1))) * Antiderivative((~a + ~b * asinh(~c * ~x)) ^ (~n + 1) / (~f + ~g * ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~h, ~m], ~x) && (EqQ(~e, (~c) ^ 2 * ~d) && (GtQ(~d, 0) && IGtQ(~n, 0)))
    @apply_utils Antiderivative(log(~(h') * (~(f') + ~(g') * ~x) ^ ~(m')) * (~d + ~(e') * (~x) ^ 2) ^ ~p * (~(a') + ~(b') * asinh(~(c') * ~x)) ^ ~(n'), ~x) => Simp((~d + ~e * (~x) ^ 2) ^ ~p / (1 + (~c) ^ 2 * (~x) ^ 2) ^ ~p) * Antiderivative(log(~h * (~f + ~g * ~x) ^ ~m) * (1 + (~c) ^ 2 * (~x) ^ 2) ^ ~p * (~a + ~b * asinh(~c * ~x)) ^ ~n, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~h, ~m, ~n], ~x) && (EqQ(~e, (~c) ^ 2 * ~d) && (IntegerQ(~p - 1 / 2) && Not(GtQ(~d, 0))))
    @apply_utils Antiderivative((~d + ~(e') * ~x) ^ ~m * (~f + ~(g') * ~x) ^ ~m * (~(a') + ~(b') * asinh(~(c') * ~x)), ~x) => With([u = IntHide((~d + ~e * ~x) ^ ~m * (~f + ~g * ~x) ^ ~m, ~x)], Dist(~a + ~b * asinh(~c * ~x), u, ~x) - ~b * ~c * Antiderivative(Dist(1 / sqrt(1 + (~c) ^ 2 * (~x) ^ 2), u, ~x), ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g], ~x) && ILtQ(~m + 1 / 2, 0)
    @apply_utils Antiderivative((~d + ~(e') * ~x) ^ ~(m') * (~f + ~(g') * ~x) ^ ~(m') * (~(a') + ~(b') * asinh(~(c') * ~x)) ^ ~(n'), ~x) => Antiderivative(ExpandIntegrand((~a + ~b * asinh(~c * ~x)) ^ ~n, (~d + ~e * ~x) ^ ~m * (~f + ~g * ~x) ^ ~m, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~n], ~x) && IntegerQ(~m)
    @apply_utils Antiderivative(~u * (~(a') + ~(b') * asinh(~(c') * ~x)), ~x) => With([v = IntHide(~u, ~x)], Dist(~a + ~b * asinh(~c * ~x), v, ~x) - ~b * ~c * Antiderivative(SimplifyIntegrand(v / sqrt(1 + (~c) ^ 2 * (~x) ^ 2), ~x), ~x) <-- InverseFunctionFreeQ(v, ~x)) <-- FreeQ([~a, ~b, ~c], ~x)
    @apply_utils Antiderivative(~Px * (~d + ~(e') * (~x) ^ 2) ^ ~p * (~(a') + ~(b') * asinh(~(c') * ~x)) ^ ~n, ~x) => With([u = ExpandIntegrand(~Px * (~d + ~e * (~x) ^ 2) ^ ~p * (~a + ~b * asinh(~c * ~x)) ^ ~n, ~x)], Antiderivative(u, ~x) <-- SumQ(u)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~n], ~x) && (PolynomialQ(~Px, ~x) && (EqQ(~e, (~c) ^ 2 * ~d) && IntegerQ(~p - 1 / 2)))
    @apply_utils Antiderivative(~(Px') * (~f + ~(g') * (~d + ~(e') * (~x) ^ 2) ^ ~p) ^ ~(m') * (~(a') + ~(b') * asinh(~(c') * ~x)) ^ ~(n'), ~x) => With([u = ExpandIntegrand(~Px * (~f + ~g * (~d + ~e * (~x) ^ 2) ^ ~p) ^ ~m * (~a + ~b * asinh(~c * ~x)) ^ ~n, ~x)], Antiderivative(u, ~x) <-- SumQ(u)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g], ~x) && (PolynomialQ(~Px, ~x) && (EqQ(~e, (~c) ^ 2 * ~d) && (IGtQ(~p + 1 / 2, 0) && IntegersQ(~m, ~n))))
    @apply_utils Antiderivative(~RFx * asinh(~(c') * ~x) ^ ~(n'), ~x) => With([u = ExpandIntegrand(asinh(~c * ~x) ^ ~n, ~RFx, ~x)], Antiderivative(u, ~x) <-- SumQ(u)) <-- FreeQ(~c, ~x) && (RationalFunctionQ(~RFx, ~x) && IGtQ(~n, 0))
    @apply_utils Antiderivative(~RFx * (~a + ~(b') * asinh(~(c') * ~x)) ^ ~(n'), ~x) => Antiderivative(ExpandIntegrand(~RFx * (~a + ~b * asinh(~c * ~x)) ^ ~n, ~x), ~x) <-- FreeQ([~a, ~b, ~c], ~x) && (RationalFunctionQ(~RFx, ~x) && IGtQ(~n, 0))
    @apply_utils Antiderivative(~RFx * (~d + ~(e') * (~x) ^ 2) ^ ~p * asinh(~(c') * ~x) ^ ~(n'), ~x) => With([u = ExpandIntegrand((~d + ~e * (~x) ^ 2) ^ ~p * asinh(~c * ~x) ^ ~n, ~RFx, ~x)], Antiderivative(u, ~x) <-- SumQ(u)) <-- FreeQ([~c, ~d, ~e], ~x) && (RationalFunctionQ(~RFx, ~x) && (IGtQ(~n, 0) && (EqQ(~e, (~c) ^ 2 * ~d) && IntegerQ(~p - 1 / 2))))
    @apply_utils Antiderivative(~RFx * (~d + ~(e') * (~x) ^ 2) ^ ~p * (~a + ~(b') * asinh(~(c') * ~x)) ^ ~(n'), ~x) => Antiderivative(ExpandIntegrand((~d + ~e * (~x) ^ 2) ^ ~p, ~RFx * (~a + ~b * asinh(~c * ~x)) ^ ~n, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x) && (RationalFunctionQ(~RFx, ~x) && (IGtQ(~n, 0) && (EqQ(~e, (~c) ^ 2 * ~d) && IntegerQ(~p - 1 / 2))))
    @apply_utils Antiderivative(~(u') * (~(a') + ~(b') * asinh(~(c') * ~x)) ^ ~(n'), ~x) => Unintegrable(~u * (~a + ~b * asinh(~c * ~x)) ^ ~n, ~x) <-- FreeQ([~a, ~b, ~c, ~n], ~x)
end

max_stars_repo_path: src/rules/7 Inverse hyperbolic functions/7.1 Inverse hyperbolic sine/7.1.5 u (a+b arcsinh(c x))^n.jl


text:
mutable struct NDMLUT{VT,VCT} <: NNLookUp
    θλ_σ::VT
    θμ_σ::VT
    θλ_σp::VT
    θμ_σp::VT
    _Π::VCT
end

lookup(net::NDM) = NDMLUT(net)
NDMLUT(net::NDM) = NDMLUT(similar(net.h_μ),
                          similar(net.h_μ),
                          similar(net.h_μ),
                          similar(net.h_μ),
                          similar(net.d_λ, complex(eltype(net.d_λ))))

function set_lookup!(lt::NDMLUT, W::NDM, c::NDMCache, σr, σc)
    ∑σ      = c.∑σ
    Δσ      = c.Δσ
    _Π_tmp  = c._Π_tmp
    T       = eltype(c.θλ_σ)

    # θλ_σ .= W.h_λ + W.w_λ*σr
    LinearAlgebra.BLAS.blascopy!(length(W.h_λ), W.h_λ, 1, lt.θλ_σ, 1)
    LinearAlgebra.BLAS.gemv!('N', one(T), W.w_λ, σr, one(T), lt.θλ_σ)

    # θμ_σ .= W.h_μ + W.w_μ*σr
    LinearAlgebra.BLAS.blascopy!(length(W.h_μ), W.h_μ, 1, lt.θμ_σ, 1)
    LinearAlgebra.BLAS.gemv!('N', one(T), W.w_μ, σr, one(T), lt.θμ_σ)

    #θλ_σp .= W.h_λ + W.w_λ*σc
    LinearAlgebra.BLAS.blascopy!(length(W.h_λ), W.h_λ, 1, lt.θλ_σp, 1)
    LinearAlgebra.BLAS.gemv!('N', T(1.0), W.w_λ, σc, T(1.0), lt.θλ_σp)

    #θμ_σp .= W.h_μ + W.w_μ*σc
    LinearAlgebra.BLAS.blascopy!(length(W.h_μ), W.h_μ, 1, lt.θμ_σp, 1)
    LinearAlgebra.BLAS.gemv!('N', T(1.0), W.w_μ, σc, T(1.0), lt.θμ_σp)

    ∑σ .= σr .+ σc
    Δσ .= σr .- σc

    #_Π         = (T(0.5)  * W.u_λ*∑σ
    #               + T(0.5)im * W.u_μ*Δσ .+ W.d_λ)
    LinearAlgebra.BLAS.blascopy!(length(W.d_λ), W.d_λ, 1, _Π_tmp, 1)
    LinearAlgebra.BLAS.gemv!('N', T(0.5), W.u_λ, ∑σ, one(T), _Π_tmp)
    lt._Π .= _Π_tmp
    LinearAlgebra.BLAS.gemv!('N', T(0.5), W.u_μ, Δσ, T(0.0), _Π_tmp)
    lt._Π .+= T(1.0)im.* _Π_tmp

    return lt
end

function update_lookup!(lt::NDMLUT, W::NDM, c::NDMCache, σr, σc,
                         changes_r, changes_c)
    θλ_σ    = lt.θλ_σ
    θμ_σ    = lt.θμ_σ
    θλ_σp   = lt.θλ_σp
    θμ_σp   = lt.θμ_σp
    Π       = lt._Π
    T       = eltype(c.θλ_σ)

    for (i, nv) = changes_r
        θλ_σ .+= W.w_λ[:,i] .* (nv - σr[i])
        θμ_σ .+= W.w_μ[:,i] .* (nv - σr[i])
        Π    .+= T(0.5)   .* W.u_λ[:,i] .* (nv - σr[i]) .+
                 T(0.5)im .* W.u_μ[:,i] .* (nv - σr[i])
    end
    for (i, nv) = changes_c
        θλ_σp .+= W.w_λ[:,i] .* (nv - σc[i])
        θμ_σp .+= W.w_μ[:,i] .* (nv - σc[i])
        Π     .+= T(0.5)   .* W.u_λ[:,i] .* (nv - σc[i]) .-
                  T(0.5)im .* W.u_μ[:,i] .* (nv - σc[i])
    end

    return lt
end

function logψ(W::NDM, c::NDMCache, lt::NDMLUT, σr, σc,
              changes_r, changes_c)
    θλ_σ    = c.θλ_σ;   θλ_σ  .= lt.θλ_σ
    θμ_σ    = c.θμ_σ;   θμ_σ  .= lt.θμ_σ
    θλ_σp   = c.θλ_σp;  θλ_σp .= lt.θλ_σp
    θμ_σp   = c.θμ_σp;  θμ_σp .= lt.θμ_σp
    Π       = c._Π;     Π     .= lt._Π
    Π_tmp   = c._Π2
    ∑σ      = c.∑σ
    Δσ      = c.Δσ
    T       = eltype(θλ_σ)

    ∑σ .= σr .+ σc
    Δσ .= σr .- σc

    for (i, nv) = changes_r
        θλ_σ .+= W.w_λ[:,i] .* (nv - σr[i])
        θμ_σ .+= W.w_μ[:,i] .* (nv - σr[i])
        Π    .+= T(0.5)   .* W.u_λ[:,i] .* (nv - σr[i]) .+
                 T(0.5)im .* W.u_μ[:,i] .* (nv - σr[i])
        ∑σ[i] += (nv - σr[i])
        Δσ[i] += (nv - σr[i])
    end
    for (i, nv) = changes_c
        θλ_σp .+= W.w_λ[:,i] .* (nv - σc[i])
        θμ_σp .+= W.w_μ[:,i] .* (nv - σc[i])
        Π     .+= T(0.5)   .* W.u_λ[:,i] .* (nv - σc[i]) .-
                  T(0.5)im .* W.u_μ[:,i] .* (nv - σc[i])
        ∑σ[i] += (nv - σc[i])
        Δσ[i] -= (nv - σc[i])
    end

    c.θλ_σ_tmp .= logℒ.(θλ_σ)
    ∑logℒ_λ_σ = sum(c.θλ_σ_tmp)

    c.θμ_σ_tmp .= logℒ.(θμ_σ)
    ∑logℒ_μ_σ = sum(c.θμ_σ_tmp)

    c.θλ_σp_tmp .= logℒ.(θλ_σp)
    ∑logℒ_λ_σp = sum(c.θλ_σp_tmp)

    c.θμ_σp_tmp .= logℒ.(θμ_σp)
    ∑logℒ_μ_σp = sum(c.θμ_σp_tmp)

    Π_tmp     .= logℒ.(Π)
    ∑Π        = sum(Π_tmp)

    Γ_λ = T(0.5) * (∑logℒ_λ_σ + ∑logℒ_λ_σp + ∑σ ⋅ W.b_λ)
    Γ_μ = T(0.5) * (∑logℒ_μ_σ - ∑logℒ_μ_σp + Δσ ⋅ W.b_μ)
    logψ = Γ_λ + T(1.0)im * Γ_μ + ∑Π

    return logψ
end

function logψ_and_∇logψ!(∇logψ, W::NDM, c::NDMCache, lt::NDMLUT,
                         σr_r, σc_r, changes_r, changes_c)
    θλ_σ    = c.θλ_σ;   θλ_σ  .= lt.θλ_σ
    θμ_σ    = c.θμ_σ;   θμ_σ  .= lt.θμ_σ
    θλ_σp   = c.θλ_σp;  θλ_σp .= lt.θλ_σp
    θμ_σp   = c.θμ_σp;  θμ_σp .= lt.θμ_σp
    Π       = c._Π;     Π     .= lt._Π
    Π_tmp   = c._Π2
    ∑σ      = c.∑σ
    Δσ      = c.Δσ
    T       = eltype(θλ_σ)

    σr = c.σr; σr .= σr_r
    σc = c.σc; σc .= σc_r
    ∑σ .= σr .+ σc
    Δσ .= σr .- σc

    for (i, nv) = changes_r
        θλ_σ .+= W.w_λ[:,i] .* (nv - σr[i])
        θμ_σ .+= W.w_μ[:,i] .* (nv - σr[i])
        Π    .+= T(0.5)   .* W.u_λ[:,i] .* (nv - σr[i]) .+
                 T(0.5)im .* W.u_μ[:,i] .* (nv - σr[i])
        ∑σ[i] += (nv - σr[i])
        Δσ[i] += (nv - σr[i])
        σr[i] = nv
    end
    for (i, nv) = changes_c
        θλ_σp .+= W.w_λ[:,i] .* (nv - σc[i])
        θμ_σp .+= W.w_μ[:,i] .* (nv - σc[i])
        Π     .+= T(0.5)   .* W.u_λ[:,i] .* (nv - σc[i]) .-
                  T(0.5)im .* W.u_μ[:,i] .* (nv - σc[i])
        ∑σ[i] += (nv - σc[i])
        Δσ[i] -= (nv - σc[i])
        σc[i] = nv
    end

    c.θλ_σ_tmp .= logℒ.(θλ_σ)
    ∑logℒ_λ_σ = sum(c.θλ_σ_tmp)

    c.θμ_σ_tmp .= logℒ.(θμ_σ)
    ∑logℒ_μ_σ = sum(c.θμ_σ_tmp)

    c.θλ_σp_tmp .= logℒ.(θλ_σp)
    ∑logℒ_λ_σp = sum(c.θλ_σp_tmp)

    c.θμ_σp_tmp .= logℒ.(θμ_σp)
    ∑logℒ_μ_σp = sum(c.θμ_σp_tmp)

    Π_tmp     .= logℒ.(Π)
    ∑Π        = sum(Π_tmp)

    Γ_λ = T(0.5) * (∑logℒ_λ_σ + ∑logℒ_λ_σp + ∑σ ⋅ W.b_λ)
    Γ_μ = T(0.5) * (∑logℒ_μ_σ - ∑logℒ_μ_σp + Δσ ⋅ W.b_μ)
    logψ = Γ_λ + T(1.0)im * Γ_μ + ∑Π

    # --- End common terms with computation of ψ --- #

    # Compute additional terms for derivatives
    ∂logℒ_λ_σ  = c.∂logℒ_λ_σ;  ∂logℒ_λ_σ  .= ∂logℒ.(θλ_σ)
    ∂logℒ_μ_σ  = c.∂logℒ_μ_σ;  ∂logℒ_μ_σ  .= ∂logℒ.(θμ_σ)
    ∂logℒ_λ_σp = c.∂logℒ_λ_σp; ∂logℒ_λ_σp .= ∂logℒ.(θλ_σp)
    ∂logℒ_μ_σp = c.∂logℒ_μ_σp; ∂logℒ_μ_σp .= ∂logℒ.(θμ_σp)
    ∂logℒ_Π    = c.∂logℒ_Π;    ∂logℒ_Π    .= ∂logℒ.(Π)

    # Store the derivatives
    ∇logψ.b_λ .= T(0.5)   .* ∑σ
    ∇logψ.b_μ .= T(0.5)im .* Δσ

    ∇logψ.h_λ .= T(0.5)   .* (∂logℒ_λ_σ .+ ∂logℒ_λ_σp)
    ∇logψ.h_μ .= T(0.5)im .* (∂logℒ_μ_σ .- ∂logℒ_μ_σp)

    ∇logψ.w_λ .= T(0.5)   .* (∂logℒ_λ_σ.*transpose(σr) .+ ∂logℒ_λ_σp.*transpose(σc))
    ∇logψ.w_μ .= T(0.5)im .* (∂logℒ_μ_σ.*transpose(σr) .- ∂logℒ_μ_σp.*transpose(σc))

    ∇logψ.d_λ .= ∂logℒ_Π
    ∇logψ.u_λ .= T(0.5) .* ∂logℒ_Π .* transpose(∑σ)
    ∇logψ.u_μ .= T(0.5)im .*  ∂logℒ_Π .* transpose(Δσ)

    return logψ
end

function Δ_logψ(W::NDM, c::NDMCache, lt::NDMLUT, σr, σc,
                     changes_r, changes_c)
    θλ_σ    = c.θλ_σ;   θλ_σ  .= lt.θλ_σ
    θμ_σ    = c.θμ_σ;   θμ_σ  .= lt.θμ_σ
    θλ_σp   = c.θλ_σp;  θλ_σp .= lt.θλ_σp
    θμ_σp   = c.θμ_σp;  θμ_σp .= lt.θμ_σp
    Π       = c._Π;     Π     .= lt._Π
    Π_tmp   = c._Π2;
    ∑σ      = c.∑σ
    Δσ      = c.Δσ
    T       = eltype(θλ_σ)

    ∑σ .= σr .+ σc
    Δσ .= σr .- σc

    logvaldiff = zero(eltype(Π))

    if !isempty(changes_r) || !isempty(changes_c)
        c.θλ_σ_tmp  .= logℒ.(θλ_σ)
        c.θμ_σ_tmp  .= logℒ.(θμ_σ)
        c.θλ_σp_tmp .= logℒ.(θλ_σp)
        c.θμ_σp_tmp .= logℒ.(θμ_σp)
        Π_tmp       .= logℒ.(Π)

        ∑logℒ_λ_σ  = sum(c.θλ_σ_tmp)
        ∑logℒ_μ_σ  = sum(c.θμ_σ_tmp)
        ∑logℒ_λ_σp = sum(c.θλ_σp_tmp)
        ∑logℒ_μ_σp = sum(c.θμ_σp_tmp)
        ∑Π_old     = sum(Π_tmp)

        Γ_λ_old = T(0.5) * (∑logℒ_λ_σ + ∑logℒ_λ_σp )
        Γ_μ_old = T(0.5) * (∑logℒ_μ_σ - ∑logℒ_μ_σp )
        logvaldiff -=  Γ_λ_old + T(1.0)im * Γ_μ_old + ∑Π_old

        for (i, nv) = changes_r
            θλ_σ .+= W.w_λ[:,i] .* (nv - σr[i])
            θμ_σ .+= W.w_μ[:,i] .* (nv - σr[i])
            Π    .+= T(0.5)   .* W.u_λ[:,i] .* (nv - σr[i]) .+
                     T(0.5)im .* W.u_μ[:,i] .* (nv - σr[i])
            logvaldiff += T(0.5)   .* W.b_λ[i] * (nv - σr[i])
            logvaldiff += T(0.5)im .* W.b_μ[i] * (nv - σr[i])
        end

        for (i, nv) = changes_c
            θλ_σp .+= W.w_λ[:,i] .* (nv - σc[i])
            θμ_σp .+= W.w_μ[:,i] .* (nv - σc[i])
            Π     .+= T(0.5)   .* W.u_λ[:,i] .* (nv - σc[i]) .-
                      T(0.5)im .* W.u_μ[:,i] .* (nv - σc[i])
            logvaldiff += T(0.5)   .* W.b_λ[i] * (nv - σc[i])
            logvaldiff -= T(0.5)im .* W.b_μ[i] * (nv - σc[i])
        end

        c.θλ_σ_tmp  .= logℒ.(θλ_σ)
        c.θμ_σ_tmp  .= logℒ.(θμ_σ)
        c.θλ_σp_tmp .= logℒ.(θλ_σp)
        c.θμ_σp_tmp .= logℒ.(θμ_σp)
        Π_tmp       .= logℒ.(Π)

        ∑logℒ_λ_σ  = sum(c.θλ_σ_tmp)
        ∑logℒ_μ_σ  = sum(c.θμ_σ_tmp)
        ∑logℒ_λ_σp = sum(c.θλ_σp_tmp)
        ∑logℒ_μ_σp = sum(c.θμ_σp_tmp)
        ∑Π         = sum(Π_tmp)

        Γ_λ = T(0.5) * (∑logℒ_λ_σ + ∑logℒ_λ_σp)
        Γ_μ = T(0.5) * (∑logℒ_μ_σ - ∑logℒ_μ_σp)
        logvaldiff +=  Γ_λ + T(1.0)im * Γ_μ + ∑Π
    end

    return logvaldiff
end

function Δ_logψ_and_∇logψ!(∇logψ, W::NDM, c::NDMCache, lt::NDMLUT,
                         σr_r, σc_r, changes_r, changes_c)
    θλ_σ    = c.θλ_σ;   θλ_σ  .= lt.θλ_σ
    θμ_σ    = c.θμ_σ;   θμ_σ  .= lt.θμ_σ
    θλ_σp   = c.θλ_σp;  θλ_σp .= lt.θλ_σp
    θμ_σp   = c.θμ_σp;  θμ_σp .= lt.θμ_σp
    Π       = c._Π;     Π     .= lt._Π
    Π_tmp   = c._Π2;
    ∑σ      = c.∑σ
    Δσ      = c.Δσ
    T       = eltype(θλ_σ)

    σr = c.σr; σr .= σr_r
    σc = c.σc; σc .= σc_r
    ∑σ .= σr .+ σc
    Δσ .= σr .- σc

    logvaldiff = zero(eltype(Π))

    if !isempty(changes_r) || !isempty(changes_c)
        c.θλ_σ_tmp  .= logℒ.(θλ_σ)
        c.θμ_σ_tmp  .= logℒ.(θμ_σ)
        c.θλ_σp_tmp .= logℒ.(θλ_σp)
        c.θμ_σp_tmp .= logℒ.(θμ_σp)
        Π_tmp       .= logℒ.(Π)

        ∑logℒ_λ_σ  = sum(c.θλ_σ_tmp)
        ∑logℒ_μ_σ  = sum(c.θμ_σ_tmp)
        ∑logℒ_λ_σp = sum(c.θλ_σp_tmp)
        ∑logℒ_μ_σp = sum(c.θμ_σp_tmp)
        ∑Π_old     = sum(Π_tmp)

        Γ_λ_old = T(0.5) * (∑logℒ_λ_σ + ∑logℒ_λ_σp )
        Γ_μ_old = T(0.5) * (∑logℒ_μ_σ - ∑logℒ_μ_σp )
        logvaldiff -=  Γ_λ_old + T(1.0)im * Γ_μ_old + ∑Π_old

        for (i, nv) = changes_r
            θλ_σ .+= W.w_λ[:,i] .* (nv - σr[i])
            θμ_σ .+= W.w_μ[:,i] .* (nv - σr[i])
            Π    .+= T(0.5)   .* W.u_λ[:,i] .* (nv - σr[i]) .+
                     T(0.5)im .* W.u_μ[:,i] .* (nv - σr[i])
            logvaldiff += T(0.5)   .* W.b_λ[i] * (nv - σr[i])
            logvaldiff += T(0.5)im .* W.b_μ[i] * (nv - σr[i])
            ∑σ[i] += (nv - σr[i])
            Δσ[i] += (nv - σr[i])
            σr[i] = nv
        end

        for (i, nv) = changes_c
            θλ_σp .+= W.w_λ[:,i] .* (nv - σc[i])
            θμ_σp .+= W.w_μ[:,i] .* (nv - σc[i])
            Π     .+= T(0.5)   .* W.u_λ[:,i] .* (nv - σc[i]) .-
                      T(0.5)im .* W.u_μ[:,i] .* (nv - σc[i])
            logvaldiff += T(0.5)   .* W.b_λ[i] * (nv - σc[i])
            logvaldiff -= T(0.5)im .* W.b_μ[i] * (nv - σc[i])
            ∑σ[i] += (nv - σc[i])
            Δσ[i] -= (nv - σc[i])
            σc[i] = nv
        end

        c.θλ_σ_tmp  .= logℒ.(θλ_σ)
        c.θμ_σ_tmp  .= logℒ.(θμ_σ)
        c.θλ_σp_tmp .= logℒ.(θλ_σp)
        c.θμ_σp_tmp .= logℒ.(θμ_σp)
        Π_tmp       .= logℒ.(Π)

        ∑logℒ_λ_σ  = sum(c.θλ_σ_tmp)
        ∑logℒ_μ_σ  = sum(c.θμ_σ_tmp)
        ∑logℒ_λ_σp = sum(c.θλ_σp_tmp)
        ∑logℒ_μ_σp = sum(c.θμ_σp_tmp)
        ∑Π         = sum(Π_tmp)

        Γ_λ = T(0.5) * (∑logℒ_λ_σ + ∑logℒ_λ_σp)
        Γ_μ = T(0.5) * (∑logℒ_μ_σ - ∑logℒ_μ_σp)
        logvaldiff +=  Γ_λ + T(1.0)im * Γ_μ + ∑Π
    end

    # Compute additional terms for derivatives
    ∂logℒ_λ_σ  = c.∂logℒ_λ_σ;  ∂logℒ_λ_σ  .= ∂logℒ.(θλ_σ)
    ∂logℒ_μ_σ  = c.∂logℒ_μ_σ;  ∂logℒ_μ_σ  .= ∂logℒ.(θμ_σ)
    ∂logℒ_λ_σp = c.∂logℒ_λ_σp; ∂logℒ_λ_σp .= ∂logℒ.(θλ_σp)
    ∂logℒ_μ_σp = c.∂logℒ_μ_σp; ∂logℒ_μ_σp .= ∂logℒ.(θμ_σp)
    ∂logℒ_Π    = c.∂logℒ_Π;    ∂logℒ_Π    .= ∂logℒ.(Π)

    # Store the derivatives
    ∇logψ.b_λ .= T(0.5)   .* ∑σ
    ∇logψ.b_μ .= T(0.5)im .* Δσ

    ∇logψ.h_λ .= T(0.5)   .* (∂logℒ_λ_σ .+ ∂logℒ_λ_σp)
    ∇logψ.h_μ .= T(0.5)im .* (∂logℒ_μ_σ .- ∂logℒ_μ_σp)

    ∇logψ.w_λ .= T(0.5)   .* (∂logℒ_λ_σ.*transpose(σr) .+ ∂logℒ_λ_σp.*transpose(σc))
    ∇logψ.w_μ .= T(0.5)im .* (∂logℒ_μ_σ.*transpose(σr) .- ∂logℒ_μ_σp.*transpose(σc))

    ∇logψ.d_λ .= ∂logℒ_Π
    ∇logψ.u_λ .= T(0.5) .* ∂logℒ_Π .* transpose(∑σ)
    ∇logψ.u_μ .= T(0.5)im .*  ∂logℒ_Π .* transpose(Δσ)

    return logvaldiff
end

ext: jl


text:
# Generic datum shifting function in ECEF coordinates

"""
    datum_shift_ECEF(dest_datum, source_datum)

Return a transformation object which transforms ECEF points from `source_datum`
to `dest_datum`.  The function should attempt to supply the current publicly
accepted best estimate.

Note that the best known version of a datum transformation will inherently
improve with time (see below), so we *cannot simultaneously* guarantee that:

1. We return the best publicly accepted version of a datum shift.
2. We return the same thing in future versions of Geodesy.jl.

This function will attempt to satisfy condition 1. rather than 2.; if you want a
stable version of a transformation you should use one of the lower level
functions, for example, `GDA94_from_ITRF_Dawson2010()`.

### Important note about accuracy

If you care about accuracy, you should *always* store long term archival data in
the source datum where possible, along with metadata defining the full datum and
coordinate system in use.  The time to do a datum shift is when you want to
compare information from two different datums.

Why all this bother about inaccuracy? The errors are twofold: First, the
parameters of a datum shift come from a physical measurement process. Typically
this involves measuring the coordinates of some physical locations in *both*
datums, and a physical measurement procedure is always subject to some inaccuracy.
Second, these **tie points** are used to infer a compact representation of the
datum shift, with as few numerical parameters as possible. A small number of
parameters will result in an overly smooth representation; this modelling error
is a second source of inaccuracy.  Both these errors can be reduced if you
choose to measure more tie points or improve the complexity of the numerical
model at a future date.
"""
datum_shift_ECEF(dest_datum, source_datum) = error("No datum shift implemented for $dest_datum from $source_datum")

# Convenience function for datums types without parameters
datum_shift_ECEF(::Type{D1}, ::Type{D2}) where {D1,D2} = datum_shift_ECEF(D1(), D2())
datum_shift_ECEF(::Type{D1}, d2) where {D1} = datum_shift_ECEF(D1(), d2)
datum_shift_ECEF(d1, ::Type{D2}) where {D2} = datum_shift_ECEF(d1, D2())

#-------------------------------------------------------------------------------
# Functions to compute transformations between known datums

"""
    GDA94_from_ITRF_Dawson2010(ITRF_year, epoch)

Return a `Transformation` converting ECEF points from a given ITRF to GDA94.
Datum shift parameters are taken from [1], supporting `ITRF_year` 2008, 2005,
2000, 1997, 1996.  `epoch` is the `Date` (or `DateTime`) of interest at which
the input `ECEF` coordinates to the transformation were measured in ITRF.

[1] J. Dawson and A. Woods, "ITRF to GDA94 coordinate transforms",
    Journal of Applied Geodesy, 4, p. 189 (2010).

TODO: We don't yet support `epoch` varying per input point, but there should be
a `Transformation` object for this at some stage.
"""
function GDA94_from_ITRF_Dawson2010(ITRF_year, epoch)
    # ITRF transformation parameters from
    # J. Dawson and A. Woods, "ITRF to GDA94 coordinate transforms",
    # Journal of Applied Geodesy, 4, p. 189 (2010):
    #
    # These are intentionally kept exactly the same as the paper for
    # consistency and ease of checking.
    #
    # (Unit key: ppb = parts per billion; mas = milli arc seconds; yr = year)
    #
    #          Values at 1994 reference epoch                                 Derivatives
    #  Name    tx      ty      tz     sc      rx       ry      rz             tx      ty    tz     sc     rx      ry      rz
    #  Units   mm      mm      mm     ppb     mas      mas     mas            mm/yr  mm/yr  mm/yr  ppb/yr mas/yr  mas/yr  mas/yr
    #
    table = Dict(
        2008 => ([-84.68, -19.42,  32.01, 9.710, -0.4254,  2.2578, 2.4015], [  1.42,  1.34,  0.90, 0.109, 1.5461, 1.1820, 1.1551]),
        2005 => ([-79.73,  -6.86,  38.03, 6.636, -0.0351,  2.1211, 2.1411], [  2.25, -0.62, -0.56, 0.294, 1.4707, 1.1443, 1.1701]),
        2000 => ([-45.91, -29.85, -20.37, 7.070, -1.6705,  0.4594, 1.9356], [ -4.66,  3.55, 11.24, 0.249, 1.7454, 1.4868, 1.2240]),
        1997 => ([-14.63, -27.62, -25.32, 6.695, -1.7893, -0.6047, 0.9962], [ -8.60,  0.36, 11.25, 0.007, 1.6394, 1.5198, 1.3801]),
        1996 => ([ 24.54, -36.43, -68.12, 6.901, -2.7359, -2.0431, 0.3731], [-21.80,  4.71, 26.27, 0.388, 2.0203, 2.1735, 1.6290]),
    )
    reference_epoch = DateTime(1994,1,1)

    if !haskey(table, ITRF_year)
        throw(ErrorException("No ITRF to GDA94 parameters for ITRF year $ITRF_year"))
    end
    params, rates = table[ITRF_year]

    # Fractional years since reference epoch.  Seems a bit unnecessarily cumbersome!
    millisecs_per_year = 365.25 * Dates.value(Dates.Millisecond(Dates.Day(1)))
    dt = Dates.value(Dates.Millisecond(convert(DateTime, epoch) - reference_epoch)) / millisecs_per_year

    # Convert units to meters and radians
    mas2rad = deg2rad(1e-3/(60*60))
    unitconv = [1e-3, 1e-3, 1e-3, 1e-9, mas2rad, mas2rad, mas2rad]
    Tx,Ty,Tz, Sc, Rx,Ry,Rz = unitconv .* (params .+ rates*dt)

    M = (1+Sc) * @SMatrix [1.0  Rz  -Ry;
                          -Rz  1.0   Rx;
                           Ry  -Rx  1.0]

    AffineMap(M, SVector(Tx,Ty,Tz))
end

datum_shift_ECEF(::GDA94, itrf::ITRF{Y}) where {Y} = GDA94_from_ITRF_Dawson2010(Y, itrf.epoch)
datum_shift_ECEF(itrf::ITRF{Y}, ::GDA94) where {Y} = inv(GDA94_from_ITRF_Dawson2010(Y, itrf.epoch))

# TODO - time-based transformation!
# make_datum_transform_ECEF{Y}(::ITRF{Y,Nothing}, ::GDA94) =


text:
#include("ps4000blockready.jl")
include("ps4000closeunit.jl")
#include("ps4000dataready.jl")
include("ps4000enumerateunits.jl")
#include("ps4000flashled.jl")
#include("ps4000getchannelinformation.jl")
#include("ps4000getmaxdownsampleratio.jl")
#include("ps4000getstreaminglatestvalues.jl")
#include("ps4000gettimebase.jl")
#include("ps4000gettimebase2.jl")
#include("ps4000gettriggerchanneltimeoffset.jl")
#include("ps4000gettriggerchanneltimeoffset64.jl")
#include("ps4000gettriggertimeoffset.jl")
#include("ps4000gettriggertimeoffset64.jl")
#include("ps4000getunitinfo.jl")
#include("ps4000getvalues.jl")
#include("ps4000getvaluesasync.jl")
#include("ps4000getvaluesbulk.jl")
#include("ps4000getvaluestriggerchanneltimeoffsetbulk.jl")
#include("ps4000getvaluestriggerchanneltimeoffsetbulk64.jl")
#include("ps4000getvaluestriggertimeoffsetbulk.jl")
#include("ps4000getvaluestriggertimeoffsetbulk64.jl")
#include("ps4000holdoff.jl")
#include("ps4000isledflashing.jl")
#include("ps4000isready.jl")
#include("ps4000istriggerorpulsewidthqualifierenabled.jl")
#include("ps4000memorysegments.jl")
#include("ps4000noofstreamingvalues.jl")
include("ps4000openunit.jl")
include("ps4000openunitasync.jl")
include("ps4000openunitasyncex.jl")
include("ps4000openunitex.jl")
include("ps4000openunitprogress.jl")
#include("ps4000runblock.jl")
#include("ps4000runstreaming.jl")
#include("ps4000runstreamingex.jl")
#include("ps4000setchannel.jl")
#include("ps4000setdatabuffer.jl")
#include("ps4000setdatabufferbulk.jl")
#include("ps4000setdatabuffers.jl")
#include("ps4000setdatabufferswithmode.jl")
#include("ps4000setdatabufferwithmode.jl")
#include("ps4000setets.jl")
#include("ps4000setetstimebuffer.jl")
#include("ps4000setetstimebuffers.jl")
#include("ps4000setexttriggerrange.jl")
#include("ps4000setpulsewidthqualifier.jl")
#include("ps4000setsiggenarbitrary.jl")
#include("ps4000setsiggenbuiltin.jl")
#include("ps4000setsimpletrigger.jl")
#include("ps4000settriggerchannelconditions.jl")
#include("ps4000settriggerchanneldirections.jl")
#include("ps4000settriggerchannelproperties.jl")
#include("ps4000settriggerdelay.jl")
#include("ps4000siggensoftwarecontrol.jl")
#include("ps4000stop.jl")
#include("ps4000streamingready.jl")


text:
const Curve_P_192 = ECP(
    p = 6277101735386680763835789423207666416083908700390324961279,
    n = 6277101735386680763835789423176059013767194773182842284081,
    #SEED = "3045ae6f c8422f64 ed579528 d38120ea e12196d5",
    #c = "3099d2bb bfcb2538 542dcd5f b078b6ef 5f3d6fe2 c745de65",
    b = "64210519 e59c80e7 0fa7e9ab 72243049 feb8deec c146b9b1",
    Gx = "188da80e b03090f6 7cbf20eb 43a18800 f4ff0afd 82ff1012",
    Gy = "07192b95 ffc8da78 631011ed 6b24cdd5 73f977a1 1e794811"
)

const Curve_P_244 = ECP(
    p = 26959946667150639794667015087019630673557916260026308143510066298881,
    n = 26959946667150639794667015087019625940457807714424391721682722368061,
    # SEED = "bd713447 99d5c7fc dc45b59f a3b9ab8f 6a948bc5",
    # c = "5b056c7e 11dd68f4 0469ee7f 3c7a7d74 f7d12111 6506d031 218291fb",
    b = "b4050a85 0c04b3ab f5413256 5044b0b7 d7bfd8ba 270b3943 2355ffb4",
    Gx = "b70e0cbd 6bb4bf7f 321390b9 4a03c1d3 56c21122 343280d6 115c1d21",
    Gy = "bd376388 b5f723fb 4c22dfe6 cd4375a0 5a074764 44d58199 85007e34"
)

const Curve_P_256 = ECP(
    p = 115792089210356248762697446949407573530086143415290314195533631308867097853951,
    n = 115792089210356248762697446949407573529996955224135760342422259061068512044369,
    #SEED = "c49d3608 86e70493 6a6678e1 139d26b7 819f7e90",
    #c = "7efba166 2985be94 03cb055c 75d4f7e0 ce8d84a9 c5114abc af317768 0104fa0d",
    b = "5ac635d8 aa3a93e7 b3ebbd55 769886bc 651d06b0 cc53b0f6 3bce3c3e 27d2604b",
    Gx = "6b17d1f2 e12c4247 f8bce6e5 63a440f2 77037d81 2deb33a0 f4a13945 d898c296",
    Gy = "4fe342e2 fe1a7f9b 8ee7eb4a 7c0f9e16 2bce3357 6b315ece cbb64068 37bf51f5"
)

const Curve_P_384 = ECP(
    p = 39402006196394479212279040100143613805079739270465446667948293404245721771496870329047266088258938001861606973112319,
    n = 39402006196394479212279040100143613805079739270465446667946905279627659399113263569398956308152294913554433653942643,
    #SEED = "a335926a a319a27a 1d00896a 6773a482 7acdac73",
    #c = "79d1e655 f868f02f ff48dcde e14151dd b80643c1 406d0ca1 0dfe6fc5 2009540a 495e8042 ea5f744f 6e184667 cc722483",
    b = "b3312fa7 e23ee7e4 988e056b e3f82d19 181d9c6e fe814112 0314088f 5013875a c656398d 8a2ed19d 2a85c8ed d3ec2aef",
    Gx = "aa87ca22 be8b0537 8eb1c71e f320ad74 6e1d3b62 8ba79b98 59f741e0 82542a38 5502f25d bf55296c 3a545e38 72760ab7",
    Gy = "3617de4a 96262c6f 5d9e98bf 9292dc29 f8f41dbd 289a147c e9da3113 b5f0b8c0 0a60b1ce 1d7e819d 7a431d7c 90ea0e5f"
)

const Curve_P_521 = ECP(
    p = 6864797660130609714981900799081393217269435300143305409394463459185543183397656052122559640661454554977296311391480858037121987999716643812574028291115057151,
    n = 6864797660130609714981900799081393217269435300143305409394463459185543183397655394245057746333217197532963996371363321113864768612440380340372808892707005449,
    #SEED = "d09e8800 291cb853 96cc6717 393284aa a0da64ba",
    #c = "0b4 8bfa5f42 0a349495 39d2bdfc 264eeeeb 077688e4 4fbf0ad8 f6d0edb3 7bd6b533 28100051 8e19f1b9 ffbe0fe9 ed8a3c22 00b8f875 e523868c 70c1e5bf 55bad637",
    b = "051 953eb961 8e1c9a1f 929a21a0 b68540ee a2da725b 99b315f3 b8b48991 8ef109e1 56193951 ec7e937b 1652c0bd 3bb1bf07 3573df88 3d2c34f1 ef451fd4 6b503f00",
    Gx = "c6 858e06b7 0404e9cd 9e3ecb66 2395b442 9c648139 053fb521 f828af60 6b4d3dba a14b5e77 efe75928 fe1dc127 a2ffa8de 3348b3c1 856a429b f97e7e31 c2e5bd66",
    Gy = "118 39296a78 9a3bc004 5c8a5fb4 2c7d1bd9 98f54449 579b4468 17afbd17 273e662c 97ee7299 5ef42640 c550b901 3fad0761 353c7086 a272c240 88be9476 9fd16650"
)

############### Koblitz curves ################

const Curve_K_163_PB = Koblitz{PB}(
    f = [163, 7, 6, 3, 0],
    a = 1,
    n = 5846006549323611672814741753598448348329118574063,
    Gx = "2 fe13c053 7bbc11ac aa07d793 de4e6d5e 5c94eee8",
    Gy = "2 89070fb0 5d38ff58 321f2e80 0536d538 ccdaa3d9"
)

const Curve_K_163_GNB = Koblitz{GNB}(
    n = 5846006549323611672814741753598448348329118574063,
    a = 1,
    m = 163,
    T = 4,
    Gx = "0 5679b353 caa46825 fea2d371 3ba450da 0c2a4541",
    Gy = "2 35b7c671 00506899 06bac3d9 dec76a83 5591edb2"
)

const Curve_K_233_PB = Koblitz{PB}(
    f = [233, 74, 0],
    a = 0,
    n = 3450873173395281893717377931138512760570940988862252126328087024741343,
    Gx = "172 32ba853a 7e731af1 29f22ff4 149563a4 19c26bf5 0a4c9d6e efad6126",
    Gy = "1db 537dece8 19b7f70f 555a67c4 27a8cd9b f18aeb9b 56e0c110 56fae6a3"
)

const Curve_K_233_GNB = Koblitz{GNB}(
    m = 233,
    T = 2,
    a = 0,
    n = 3450873173395281893717377931138512760570940988862252126328087024741343,
    Gx = "0fd e76d9dcd 26e643ac 26f1aa90 1aa12978 4b71fc07 22b2d056 14d650b3",
    Gy = "064 3e317633 155c9e04 47ba8020 a3c43177 450ee036 d6335014 34cac978"
)

const Curve_K_283_PB = Koblitz{PB}(
    f = [283, 12, 7, 5, 0],
    a = 0,
    n = 3885337784451458141838923813647037813284811733793061324295874997529815829704422603873,
    Gx = "503213f 78ca4488 3f1a3b81 62f188e5 53cd265f 23c1567a 16876913 b0c2ac24 58492836",
    Gy = "1ccda38 0f1c9e31 8d90f95d 07e5426f e87e45c0 e8184698 e4596236 4e341161 77dd2259"
)

const Curve_K_283_GNB = Koblitz{GNB}(
    m = 283,
    T = 6,
    a = 0,
    n = 3885337784451458141838923813647037813284811733793061324295874997529815829704422603873,
    Gx = "3ab9593 f8db09fc 188f1d7c 4ac9fcc3 e57fcd3b db15024b 212c7022 9de5fcd9 2eb0ea60",
    Gy = "2118c47 55e7345c d8f603ef 93b98b10 6fe8854f feb9a3b3 04634cc8 3a0e759f 0c2686b1"
)

const Curve_K_409_PB = Koblitz{PB}(
    f = [409, 87, 0],
    a = 0,
     n = 330527984395124299475957654016385519914202341482140609642324395022880711289249191050673258457777458014096366590617731358671,
    Gx = "060f05f 658f49c1 ad3ab189 0f718421 0efd0987 e307c84c 27accfb8 f9f67cc2 c460189e b5aaaa62 ee222eb1 b35540cf e9023746",
    Gy = "1e36905 0b7c4e42 acba1dac bf04299c 3460782f 918ea427 e6325165 e9ea10e3 da5f6c42 e9c55215 aa9ca27a 5863ec48 d8e0286b"
)

const Curve_K_409_GNB = Koblitz{GNB}(
    m = 409,
    T = 4,
    a = 0,
    n = 330527984395124299475957654016385519914202341482140609642324395022880711289249191050673258457777458014096366590617731358671,
    Gx = "1b559c7 cba2422e 3affe133 43e808b5 5e012d72 6ca0b7e6 a63aeafb c1e3a98e 10ca0fcf 98350c3b 7f89a975 4a8e1dc0 713cec4a",
    Gy = "16d8c42 052f07e7 713e7490 eff318ba 1abd6fef 8a5433c8 94b24f5c 817aeb79 852496fb ee803a47 bc8a2038 78ebf1c4 99afd7d6"
)

const Curve_K_571_PB = Koblitz{PB}(
    f = [571, 10, 5, 2, 0],
    n = 1932268761508629172347675945465993672149463664853217499328617625725759571144780212268133978522706711834706712800825351461273674974066617311929682421617092503555733685276673,
    a = 0,
    Gx = "26eb7a8 59923fbc 82189631 f8103fe4 ac9ca297 0012d5d4 60248048 01841ca4 43709584 93b205e6 47da304d b4ceb08c bbd1ba39 494776fb 988b4717 4dca88c7 e2945283 a01c8972",
    Gy = "349dc80 7f4fbf37 4f4aeade 3bca9531 4dd58cec 9f307a54 ffc61efc 006d8a2c 9d4979c0 ac44aea7 4fbebbb9 f772aedc b620b01a 7ba7af1b 320430c8 591984f6 01cd4c14 3ef1c7a3"
)

const Curve_K_571_GNB = Koblitz{GNB}(
    m = 571,
    T = 10,
    a = 0,
    n = 1932268761508629172347675945465993672149463664853217499328617625725759571144780212268133978522706711834706712800825351461273674974066617311929682421617092503555733685276673,
    Gx = "04bb2db a418d0db 107adae0 03427e5d 7cc139ac b465e593 4f0bea2a b2f3622b c29b3d5b 9aa7a1fd fd5d8be6 6057c100 8e71e484 bcd98f22 bf847642 37673674 29ef2ec5 bc3ebcf7",
    Gy = "44cbb57 de20788d 2c952d7b 56cf39bd 3e89b189 84bd124e 751ceff4 369dd8da c6a59e6e 745df44d 8220ce22 aa2c852c fcbbef49 ebaa98bd 2483e331 80e04286 feaa2530 50caff60"
)

##################### Binary curves ############################

const Curve_B_163_PB = EC2N{PB}(
    f = [163, 7, 6, 3, 0],
    a = 1,
    b = "2 0a601907 b8c953ca 1481eb10 512f7874 4a3205fd",
    n = 5846006549323611672814742442876390689256843201587,
    Gx = "3 f0eba162 86a2d57e a0991168 d4994637 e8343e36",
    Gy = "0 d51fbc6c 71a0094f a2cdd545 b11c5c0c 797324f1"
)

const Curve_B_163_GNB = EC2N{GNB}(
    m = 163,
    T = 4,
    a = 1,
    b = "6 645f3cac f1638e13 9c6cd13e f61734fb c9e3d9fb",
    n = 5846006549323611672814742442876390689256843201587,
    Gx = "0 311103c1 7167564a ce77ccb0 9c681f88 6ba54ee8",
    Gy = "3 33ac13c6 447f2e67 613bf700 9daf98c8 7bb50c7f"
)

const Curve_B_233_PB = EC2N{PB}(
    f = [233, 74, 0],
    a = 1,
    b = "066 647ede6c 332c7f8c 0923bb58 213b333b 20e9ce42 81fe115f 7d8f90ad",
    n = 6901746346790563787434755862277025555839812737345013555379383634485463,
    Gx = "0fa c9dfcbac 8313bb21 39f1bb75 5fef65bc 391f8b36 f8f8eb73 71fd558b",
    Gy = "100 6a08a419 03350678 e58528be bf8a0bef f867a7ca 36716f7e 01f81052"
)

const Curve_B_233_GNB = EC2N{GNB}(
    m = 233,
    T = 2,
    a = 1,
    b = "1a0 03e0962d 4f9a8e40 7c904a95 38163adb 82521260 0c7752ad 52233279",
    n = 6901746346790563787434755862277025555839812737345013555379383634485463,
    Gx = "18b 863524b3 cdfefb94 f2784e0b 116faac5 4404bc91 62a363ba b84a14c5",
    Gy = "049 25df77bd 8b8ff1a5 ff519417 822bfedf 2bbd7526 44292c98 c7af6e02"
)

const Curve_B_283_PB = EC2N{PB}(
    f = [283, 12, 7, 5, 0],
    n = 770675568902916283677847627294075626569625924376904889109196526770044277787378692871,
    a = 1,
    b = "27b680a c8b8596d a5a4af8a 19a0303f ca97fd76 45309fa2 a581485a f6263e31 3b79a2f5",
    Gx = "5f93925 8db7dd90 e1934f8c 70b0dfec 2eed25b8 557eac9c 80e2e198 f8cdbecd 86b12053",
    Gy = "3676854 fe24141c b98fe6d4 b20d02b4 516ff702 350eddb0 826779c8 13f0df45 be8112f4"
)

const Curve_B_283_GNB = EC2N{GNB}(
    m = 283,
    T = 6,
    a = 1,
    b = "157261b 894739fb 5a13503f 55f0b3f1 0c560116 66331022 01138cc1 80c0206b dafbc951",
    n = 7770675568902916283677847627294075626569625924376904889109196526770044277787378692871,
    Gx = "749468e 464ee468 634b21f7 f61cb700 701817e6 bc36a236 4cb8906e 940948ea a463c35d",
    Gy = "62968bd 3b489ac5 c9b859da 68475c31 5bafcdc4 ccd0dc90 5b70f624 46f49c05 2f49c08c"
)

const Curve_B_409_PB = EC2N{PB}(
    f = [409, 87, 0],
    n = 661055968790248598951915308032771039828404682964281219284648798304157774827374805208143723762179110965979867288366567526771,
    a = 1,
    b = "021a5c2 c8ee9feb 5c4b9a75 3b7b476b 7fd6422e f1f3dd67 4761fa99 d6ac27c8 a9a197b2 72822f6c d57a55aa 4f50ae31 7b13545f",
    Gx = "15d4860 d088ddb3 496b0c60 64756260 441cde4a f1771d4d b01ffe5b 34e59703 dc255a86 8a118051 5603aeab 60794e54 bb7996a7",
    Gy = "061b1cf ab6be5f3 2bbfa783 24ed106a 7636b9c5 a7bd198d 0158aa4f 5488d08f 38514f1f df4b4f40 d2181b36 81c364ba 0273c706"
)

const Curve_B_409_GNB = EC2N{GNB}(
    m = 409,
    T = 4,
    a = 1,
    b = "124d065 1c3d3772 f7f5a1fe 6e715559 e2129bdf a04d52f7 b6ac7c53 2cf0ed06 f610072d 88ad2fdc c50c6fde 72843670 f8b3742a",
    n = 661055968790248598951915308032771039828404682964281219284648798304157774827374805208143723762179110965979867288366567526771,
    Gx = "0ceacbc 9f475767 d8e69f3b 5dfab398 13685262 bcacf22b 84c7b6dd 981899e7 318c96f0 761f77c6 02c016ce d7c548de 830d708f",
    Gy = "199d64b a8f089c6 db0e0b61 e80bb959 34afd0ca f2e8be76 d1c5e9af fc7476df 49142691 ad303902 88aa09bc c59c1573 aa3c009a"
)

const Curve_B_571_PB = EC2N{PB}(
    f = [571, 10, 5, 2, 0],
    n = 3864537523017258344695351890931987344298927329706434998657235251451519142289560424536143999389415773083133881121926944486246872462816813070234528288303332411393191105285703,
    a = 1,
    b = "2f40e7e 2221f295 de297117 b7f3d62f 5c6a97ff cb8ceff1 cd6ba8ce 4a9a18ad 84ffabbd 8efa5933 2be7ad67 56a66e29 4afd185a 78ff12aa 520e4de7 39baca0c 7ffeff7f 2955727a",
    Gx = "303001d 34b85629 6c16c0d4 0d3cd775 0a93d1d2 955fa80a a5f40fc8 db7b2abd bde53950 f4c0d293 cdd711a3 5b67fb14 99ae6003 8614f139 4abfa3b4 c850d927 e1e7769c 8eec2d19",
    Gy = "37bf273 42da639b 6dccfffe b73d69d7 8c6c27a6 009cbbca 1980f853 3921e8a6 84423e43 bab08a57 6291af8f 461bb2a8 b3531d2f 0485c19b 16e2f151 6e23dd3c 1a4827af 1b8ac15b"
)

const Curve_B_571_GNB = EC2N{GNB}(
    m = 571,
    T = 10,
    a = 1,
    b = "3762d0d 47116006 179da356 88eeaccf 591a5cde a7500011 8d9608c5 9132d434 26101a1d fb377411 5f586623 f75f0000 1ce61198 3c1275fa 31f5bc9f 4be1a0f4 67f01ca8 85c74777",
    n = 3864537523017258344695351890931987344298927329706434998657235251451519142289560424536143999389415773083133881121926944486246872462816813070234528288303332411393191105285703,
    Gx = "0735e03 5def5925 cc33173e b2a8ce77 67522b46 6d278b65 0a291612 7dfea9d2 d361089f 0a7a0247 a184e1c7 0d417866 e0fe0feb 0ff8f2f3 f9176418 f97d117e 624e2015 df1662a8",
    Gy = "04a3642 0572616c df7e606f ccadaecf c3b76dab 0eb1248d d03fbdfc 9cd3242c 4726be57 9855e812 de7ec5c5 00b4576a 24628048 b6a72d88 0062eed0 dd34b109 6d3acbb6 b01a4a97"
)


text:
@info "Starting ModelNet tests..."
@testset "ModelNet10 PointCloud dataset" begin

    mn10 = ModelNet10(mode="point_cloud")
    for (split, train) in [
            ("train", true),
            ("test", false)
        ]

        dset = mn10(train=train, npoints=1024)

        @test dset isa Flux3D.Dataset.AbstractDataset
        @test dset.root == normpath(@__DIR__, "..", "datasets")
        @test dset.path == normpath(dset.root, "modelnet40_normal_resampled")
        @test dset.train == train
        @test dset.npoints == 1024
        # @test dset.sampling   #TODO add appropriate test
        # @test dset.transform   #TODO add appropriate test
        @test dset.datapaths isa AbstractArray

        if train
            @test dset.length == 3991
            @test size(dset) == (3991,)
            @test length(dset) == 3991
        else
            @test dset.length == 908
            @test size(dset) == (908,)
            @test length(dset) == 908
        end

        dpoint = dset[1]
        @test dpoint isa Dataset.AbstractDataPoint
        @test dpoint.data isa PointCloud
    end
end

@testset "ModelNet40 PointCloud dataset" begin

    mn40 = ModelNet40(mode="point_cloud")
    for (split, train) in [
            ("train", true),
            ("test", false)
        ]

        dset = mn40(train=train, npoints=1024)

        @test dset isa Flux3D.Dataset.AbstractDataset
        @test dset.root == normpath(@__DIR__, "..", "datasets")
        @test dset.path == normpath(dset.root, "modelnet40_normal_resampled")
        @test dset.train == train
        @test dset.npoints == 1024
        # @test dset.sampling   #TODO add appropriate test
        # @test dset.transform   #TODO add appropriate test
        @test dset.datapaths isa AbstractArray

        if train
            @test dset.length == 9843
            @test size(dset) == (9843,)
            @test length(dset) == 9843
        else
            @test dset.length == 2468
            @test size(dset) == (2468,)
            @test length(dset) == 2468
        end

        dpoint = dset[1]
        @test dpoint isa Dataset.AbstractDataPoint
        @test dpoint.data isa PointCloud
    end
end


text:
include("./shared.jl")
include("../scripts-part2/1D_SWE/shallow_water_1D_LF_xpu.jl")
include("../scripts-part2/1D_SWE/shallow_water_1D_LF.jl")

@testset "Height H between baseline and 1D xpu implementation" begin
      H_xpu = shallow_water_1D_xpu(; nx = n)
      H_cpu = shallow_water_1D(; nx = n)
      @test isapprox(H_xpu, H_cpu)
end;


text:
using JuLIP

import JuLIPMaterials
reload("JuLIPMaterials")
JM = JuLIPMaterials

at = bulk(:Fe)
calc = SitePotential(LennardJones(r0 = rnn(:Fe)) * C2Shift(2.5 * rnn(:Fe)))

lGF = JM.LGFs.LGF(calc, at)

# using JuLIPMaterials: Vec3, Mat3
# R = rand(SVec{Float64}, 10)
# H = rand(SMat{Float64}, 10)


text:
using Point2DTools
using SVG

function conv_svg(fname, points, upper, lower)
    cs = map(p->Circle(p.x, p.y, 2), points)
    hull = [Polyline(points[upper]; style=Style(;strokecolor="green")), Polyline(points[lower]; style=Style(;strokecolor="red"))]
    write("$fname.html", SVG.Svg(vcat(cs, hull)), 500, 500; inhtml = true)
end

function rand_svg()
    points(n=16) = 500 * rand(Point2D, n)
    ps = collect(points())
    conv_svg("rand", ps, convexperm(ps)...)
end

function rpoints_svg()
    rpoints = Point2D[[408, 22], [424, 65], [484, 416], [283, 161], [98, 358], [126, 288], [309, 56], [63, 410], [221, 8] ,[71, 378] ,[73, 138], [305, 105], [397, 475], [376, 212], [292, 304], [498, 303]])
    conv_svg("rpoints", rpoints, convexperm(rpoints)...)
end


text:
using Test
using LinearAlgebra
using ProximalOperators: NormL1
using ProximalAlgorithms

@testset "Autodiff ($T)" for T in [Float32, Float64, ComplexF32, ComplexF64]
    R = real(T)
    A = T[
        1.0 -2.0 3.0 -4.0 5.0
        2.0 -1.0 0.0 -1.0 3.0
        -1.0 0.0 4.0 -3.0 2.0
        -1.0 -1.0 -1.0 1.0 3.0
    ]
    b = T[1.0, 2.0, 3.0, 4.0]
    f(x) = R(1/2) * norm(A * x - b, 2)^2
    Lf = opnorm(A)^2
    m, n = size(A)

    @testset "Gradient" begin
        x = randn(T, n)
        gradfx, fx = ProximalAlgorithms.gradient(f, x)
        @test eltype(gradfx) == T
        @test typeof(fx) == R
        @test gradfx ≈ A' * (A * x - b)
    end

    @testset "Algorithms" begin
        lam = R(0.1) * norm(A' * b, Inf)
        @test typeof(lam) == R
        g = NormL1(lam)
        x_star = T[-3.877278911564627e-01, 0, 0, 2.174149659863943e-02, 6.168435374149660e-01]
        TOL = R(1e-4)
        solver = ProximalAlgorithms.FastForwardBackward(tol = TOL)
        x, it = solver(x0 = zeros(T, n), f = f, g = g, Lf = Lf)
        @test eltype(x) == T
        @test norm(x - x_star, Inf) <= TOL
        @test it < 100
    end
end


text:
# With node.

mutable struct WithNode <: TabularNode
    over::Union{SQLNode, Nothing}
    args::Vector{SQLNode}
    label_map::OrderedDict{Symbol, Int}

    function WithNode(; over = nothing, args, label_map = nothing)
        if label_map !== nothing
            return new(over, args, label_map)
        end
        n = new(over, args, OrderedDict{Symbol, Int}())
        for (i, arg) in enumerate(n.args)
            name = label(arg)
            if name in keys(n.label_map)
                err = DuplicateLabelError(name, path = [arg, n])
                throw(err)
            end
            n.label_map[name] = i
        end
        n
    end
end

WithNode(args...; over = nothing) =
    WithNode(over = over, args = SQLNode[args...])

"""
    With(; over = nothing, args)
    With(args...; over = nothing)

`With` assigns a name to a temporary dataset.  This dataset could be
referred to by name in the `over` query.

```
WITH \$args...
SELECT ...
FROM \$over
```

# Examples

```jldoctest
julia> person = SQLTable(:person, columns = [:person_id, :year_of_birth]);

julia> condition_occurrence =
           SQLTable(:condition_occurrence, columns = [:condition_occurrence_id,
                                                      :person_id,
                                                      :condition_concept_id]);

julia> q = From(person) |>
           Where(Fun.in(Get.person_id, From(:essential_hypertension) |>
                                       Select(Get.person_id))) |>
           With(:essential_hypertension =>
                    From(condition_occurrence) |>
                    Where(Get.condition_concept_id .== 320128));

julia> print(render(q))
WITH "essential_hypertension_1" ("person_id") AS (
  SELECT "condition_occurrence_1"."person_id"
  FROM "condition_occurrence" AS "condition_occurrence_1"
  WHERE ("condition_occurrence_1"."condition_concept_id" = 320128)
)
SELECT
  "person_1"."person_id",
  "person_1"."year_of_birth"
FROM "person" AS "person_1"
WHERE ("person_1"."person_id" IN (
  SELECT "essential_hypertension_1"."person_id"
  FROM "essential_hypertension_1"
))
```
"""
With(args...; kws...) =
    WithNode(args...; kws...) |> SQLNode

dissect(scr::Symbol, ::typeof(With), pats::Vector{Any}) =
    dissect(scr, WithNode, pats)

function PrettyPrinting.quoteof(n::WithNode, ctx::QuoteContext)
    ex = Expr(:call, nameof(With))
    if isempty(n.args)
        push!(ex.args, Expr(:kw, :args, Expr(:vect)))
    else
        append!(ex.args, quoteof(n.args, ctx))
    end
    if n.over !== nothing
        ex = Expr(:call, :|>, quoteof(n.over, ctx), ex)
    end
    ex
end

label(n::WithNode) =
    label(n.over)

rebase(n::WithNode, n′) =
    WithNode(over = rebase(n.over, n′), args = n.args)


text:
using Manifolds, Manopt, Test, ManifoldsBase, Dates

struct TestProblem <: Problem{AllocatingEvaluation} end
struct TestOptions <: Options end

@testset "StoppingCriteria" begin
    struct myStoppingCriteriaSet <: StoppingCriterionSet end
    @test_throws ErrorException get_stopping_criteria(myStoppingCriteriaSet())

    s = StopWhenAll(StopAfterIteration(10), StopWhenChangeLess(0.1))
    s2 = StopWhenAll([StopAfterIteration(10), StopWhenChangeLess(0.1)])
    @test get_stopping_criteria(s)[1].maxIter == get_stopping_criteria(s2)[1].maxIter

    s3 = StopWhenCostLess(0.1)
    p = GradientProblem(Euclidean(), (M, x) -> x^2, x -> 2x)
    o = GradientDescentOptions(Euclidean(), 1.0)
    @test !s3(p, o, 1)
    @test length(s3.reason) == 0
    o.x = 0.3
    @test s3(p, o, 2)
    @test length(s3.reason) > 0
    # repack
    sn = StopWhenAny(StopAfterIteration(10), s3)
    sn2 = StopWhenAny([StopAfterIteration(10), s3])
    @test get_stopping_criteria(sn)[1].maxIter == get_stopping_criteria(sn2)[1].maxIter
    @test get_stopping_criteria(sn)[2].threshold == get_stopping_criteria(sn2)[2].threshold
    # then s3 is the only active one
    @test get_active_stopping_criteria(sn) == [s3]
    @test get_active_stopping_criteria(s3) == [s3]
    @test get_active_stopping_criteria(StopAfterIteration(1)) == []
    sm = StopWhenAll(StopAfterIteration(10), s3)
    @test !sm(p, o, 9)
    @test sm(p, o, 11)
    an = sm.reason
    m = match(r"^((.*)\n)+", an)
    @test length(m.captures) == 2 # both have to be active
    update_stopping_criterion!(s3, :MinCost, 1e-2)
    @test s3.threshold == 1e-2
end

@testset "Test StopAfter" begin
    p = TestProblem()
    o = TestOptions()
    s = StopAfter(Second(1))
    s(p, o, 0) # Start
    @test s(p, o, 1) == false
    sleep(1.02)
    @test s(p, o, 2) == true
    @test_throws ErrorException StopAfter(Second(-1))
    @test_throws ErrorException update_stopping_criterion!(s, :MaxTime, Second(-1))
    update_stopping_criterion!(s, :MaxTime, Second(2))
    @test s.threshold == Second(2)
end

@testset "Stopping Criterion &/| operators" begin
    a = StopAfterIteration(200)
    b = StopWhenChangeLess(1e-6)
    c = StopWhenGradientNormLess(1e-6)
    d = StopWhenAll(a, b, c)
    @test typeof(d) === typeof(a & b & c)
    @test typeof(d) === typeof(a & (b & c))
    @test typeof(d) === typeof((a & b) & c)
    update_stopping_criterion!(d, :MinIterateChange, 1e-8)
    @test d.criteria[2].threshold == 1e-8
    e = StopWhenAny(a, b, c)
    @test typeof(e) === typeof(a | b | c)
    @test typeof(e) === typeof(a | (b | c))
    @test typeof(e) === typeof((a | b) | c)
    update_stopping_criterion!(e, :MinGradNorm, 1e-9)
    @test d.criteria[3].threshold == 1e-9
end

@testset "TCG stopping criteria" begin
    # create dummy criterion
    p = HessianProblem(Euclidean(), x -> x, (M, x) -> x, (M, x) -> x, x -> x)
    o = TruncatedConjugateGradientOptions(
        Euclidean(), 1.0, 0.0; trust_region_radius=2.0, randomize=false
    )
    o.new_model_value = 2.0
    o.model_value = 1.0
    s = StopWhenModelIncreased()
    @test !s(p, o, 0)
    @test s.reason == ""
    @test s(p, o, 1)
    @test length(s.reason) > 0
    s2 = StopWhenCurvatureIsNegative()
    o.δHδ = -1.0
    @test !s2(p, o, 0)
    @test s2.reason == ""
    @test s2(p, o, 1)
    @test length(s2.reason) > 0
    s3 = StopIfResidualIsReducedByPower(1.0)
    update_stopping_criterion!(s3, :ResidualPower, 0.5)
    @test s3.θ == 0.5
end

@testset "Stop with step size" begin
    p = GradientProblem(Euclidean(), (M, x) -> x^2, x -> 2x)
    o = GradientDescentOptions(Euclidean(), 1.0)
    s1 = StopWhenStepsizeLess(0.5)
    @test !s1(p, o, 1)
    @test s1.reason == ""
    o.stepsize = ConstantStepsize(; stepsize=0.25)
    @test s1(p, o, 2)
    @test length(s1.reason) > 0
    update_stopping_criterion!(o, :MaxIteration, 200)
    @test o.stop.maxIter == 200
    update_stopping_criterion!(s1, :MinStepsize, 1e-1)
    @test s1.threshold == 1e-1
end


text:
function get_anxious_index(outfile = "temp.xlsx")
  url = "https://www.philadelphiafed.org/-/media/research-and-data/real-time-center/survey-of-professional-forecasters/anxious-index/anxious_index_chart.xlsx?la=en"

  urldownload(url, parser = identity, save_raw = outfile)
  df = DataFrame(XLSX.readtable(outfile, "Data")...)
  outfile == "temp.xlsx" ? rm(outfile) : nothing

  return df
end


text:
# This file is a part of TypeDBClient.  License is MIT: https://github.com/Humans-of-Julia/TypeDBClient.jl/blob/main/LICENSE

const PULSE_INTERVAL_MILLIS = 5000

mutable struct  CoreSession <: AbstractCoreSession
    client::CoreClient
    database::CoreDatabase
    sessionID::Bytes
    transactions::Dict{UUID,T} where {T<:Optional{<:AbstractCoreTransaction}}
    type::EnumType
    accessLock::ReentrantLock
    options::TypeDBOptions
    isOpen::Bool
    networkLatencyMillis::Int
    timer::Optional{Controller}
    request_timeout::Float64
    error_break_time::Float64
end

Base.show(io::IO, session::AbstractCoreSession) = print(io, "Session(ID: $(bytes2hex(session.sessionID)))")

function CoreSession(client::AbstractCoreClient,
                     database::AbstractString,
                     type::Int32,
                     options::TypeDBOptions = TypeDBOptions();
                     request_timeout::Float64=Inf,
                     error_time::Float64 = 5.0)
    try
        #building open_request
        open_req = SessionRequestBuilder.open_req(
            database, type , copy_to_proto(options, Proto.Options)
        )
        # open the session
        startTime = now()
        grpc_controller = gRPCController()
        req_result, status  = Proto.session_open(client.core_stub.blockingStub, grpc_controller, open_req)
        res_id = grpc_result_or_error(req_result, status, result->result.session_id)
        endTime = now()

        database = CoreDatabase(database)
        networkLatencyMillis = (endTime - startTime).value
        session_id =  res_id
        transactions = Dict{UUID,AbstractCoreTransaction}()
        is_open = true

        t = Controller(true, 0.5)

        result = CoreSession(client,
                    database,
                    session_id,
                    transactions,
                    type,
                    ReentrantLock(),
                    options,
                    is_open,
                    networkLatencyMillis,
                    t,
                    request_timeout,
                    error_time)

        # initialize the ongoing pulse request
        make_pulse_request(result, t)

        # add the session to the client
        client.sessions[session_id] = result

        # the following transaction is useless for the construction of the the Session
        # and is only to initialize the first compilation of the ProtoBuf machinary
        # underneath the transaction. In the environment here it is possible to control
        # the tick of the pulse request a close as can be. If somebody find a way to
        # avoid the long starting time of the transaction this can be removed.
        trans = nothing
        try
            trans = transaction(result, Proto.Transaction_Type[:READ])
        catch ex
            @debug "First attempt for transaction done and not successful"
        end
        trans !== nothing && close(trans)

        t.duration_in_seconds =  (PULSE_INTERVAL_MILLIS / 1000) - 0.5

        return result
    catch ex
        throw(TypeDBClientException("Error construct a CoreSession",ex))
    end
end

"""
    make_pulse_request(session::T) where {T<:AbstractCoreSession}

This function make a pulse request to keep the session alive.
"""
function make_pulse_request(session::AbstractCoreSession, controller::Controller)
    @async begin
        while controller.running
            try
                pulsreq = SessionRequestBuilder.pulse_req(session.sessionID)
                req_result, status = Proto.session_pulse(session.client.core_stub.blockingStub, gRPCController() , pulsreq)
                result = grpc_result_or_error(req_result,status, result->result)
                @debug "Time: $(Dates.now())"

                if result.alive === false
                    close(session, false)
                    @debug "$session is closed"
                end
            catch ex
                close(session)
                throw(ex)
            end
            sleep(controller.duration_in_seconds)
        end
    end
end

transaction(session::AbstractCoreSession,
            type::EnumType;
            error_break_time::Real = session.error_break_time) =
                transaction(session,
                    type,
                    typedb_options_core(),
                    true,
                    error_break_time = error_break_time)

transaction(session::AbstractCoreSession,
    type::EnumType,
    blocking::Bool) = transaction(session,
                        type,
                        typedb_options_core(),
                        blocking)

function transaction(session::AbstractCoreSession,
            type::EnumType,
            options::TypeDBOptions,
            blocking::Bool = true;
            error_break_time::Real = session.error_break_time)

    lock(session.accessLock) do
        if !session.isOpen
            throw(TypeDBClientException(CLIENT_SESSION_CLOSED, bytes2hex(session.sessionID)))
        end

        transactionRPC = CoreTransaction(session,
                            session.sessionID,
                            type,
                            options,
                            use_blocking_stub = blocking,
                            error_break_time = error_break_time)

        session.transactions[transactionRPC.transaction_id] = transactionRPC

        return transactionRPC
    end
end

function close(session::AbstractCoreSession, session_alive::Bool=true)

        lock(session.accessLock) do
            if session.isOpen
                for (_,trans) in session.transactions
                    safe_close(trans)
                end
                remove_session(session.client, session)
                safe_close(session.timer)

                if session_alive
                    req = SessionRequestBuilder.close_req(session.sessionID)
                    stub = session.client.core_stub.blockingStub
                    Proto.session_close(stub, gRPCController(), req )
                end

                session.isOpen = false
            end
        end

    return true
end

function Base.delete!(session::AbstractCoreSession, trans_id::UUID)

        lock(session.accessLock) do
            delete!(session.transactions, trans_id)
        end

end

function is_open(session::AbstractCoreSession)
    return session.isOpen
end


text:
# Next, you should verify the life support rating, which can be determined by
# multiplying the oxygen generator rating by the CO2 scrubber rating.
#
# Both the oxygen generator rating and the CO2 scrubber rating are values that
# can be found in your diagnostic report - finding them is the tricky part.
# Both values are located using a similar process that involves filtering out
# values until only one remains. Before searching for either rating value,
# start with the full list of binary numbers from your diagnostic report and
# consider just the first bit of those numbers. Then:
#
# Keep only numbers selected by the bit criteria for the type of rating value
# for which you are searching. Discard numbers which do not match the bit
# criteria.
#
# - If you only have one number left, stop; this is the rating value for which
#   you are searching.
# - Otherwise, repeat the process, considering the next bit to the right.
# - The bit criteria depends on which type of rating value you want to find:
#
# To find oxygen generator rating, determine the most common value (0 or 1) in
# the current bit position, and keep only numbers with that bit in that
# position. If 0 and 1 are equally common, keep values with a 1 in the position
# being considered.
#
# To find CO2 scrubber rating, determine the least common value (0 or 1) in the
# current bit position, and keep only numbers with that bit in that position.
# If 0 and 1 are equally common, keep values with a 0 in the position being
# considered.
#
# For example, to determine the oxygen generator rating value using the same
# example diagnostic report from above:
#
# Start with all 12 numbers and consider only the first bit of each number.
# There are more 1 bits (7) than 0 bits (5), so keep only the 7 numbers with a
# 1 in the first position: 11110, 10110, 10111, 10101, 11100, 10000, and 11001.
#
# Then, consider the second bit of the 7 remaining numbers: there are more 0
# bits (4) than 1 bits (3), so keep only the 4 numbers with a 0 in the second
# position: 10110, 10111, 10101, and 10000.
#
# In the third position, three of the four numbers have a 1, so keep those
# three: 10110, 10111, and 10101.
#
# In the fourth position, two of the three numbers have a 1, so keep those two:
# 10110 and 10111.
#
# In the fifth position, there are an equal number of 0 bits and 1 bits (one
# each). So, to find the oxygen generator rating, keep the number with a 1 in
# that position: 10111.
#
# As there is only one number left, stop; the oxygen generator rating is 10111,
# or 23 in decimal.
#
# Then, to determine the CO2 scrubber rating value from the same example above:
#
# Start again with all 12 numbers and consider only the first bit of each
# number. There are fewer 0 bits (5) than 1 bits (7), so keep only the 5
# numbers with a 0 in the first position: 00100, 01111, 00111, 00010, and
# 01010.
#
# Then, consider the second bit of the 5 remaining numbers: there are fewer 1
# bits (2) than 0 bits (3), so keep only the 2 numbers with a 1 in the second
# position: 01111 and 01010.
#
# In the third position, there are an equal number of 0 bits and 1 bits (one
# each). So, to find the CO2 scrubber rating, keep the number with a 0 in that
# position: 01010.
#
# As there is only one number left, stop; the CO2 scrubber rating is 01010, or
# 10 in decimal.
#
# Finally, to find the life support rating, multiply the oxygen generator
# rating (23) by the CO2 scrubber rating (10) to get 230.
#
# Use the binary numbers in your diagnostic report to calculate the oxygen
# generator rating and CO2 scrubber rating, then multiply them together. What
# is the life support rating of the submarine? (Be sure to represent your
# answer in decimal, not binary.)

# used with searchsortedfirst below - because we're using by=, it needs to be
# the same length as a line of input
const srchterm = "111111111111"

# read and sort the lines - the sort makes it easy to find the most common bit
# in a column, and which values have that
lines = sort!(readlines("input.txt"))

"""
    dig(v, mc, colidx=1)

Recursively digs through data to find the datum according to the algo above.
`v` is the data we're considering, `mc` is true to search for the most common
bit or false otherwise, and `colidx` is which digit we're considering (1 being
the first digit and recursively increasing to 12 as the algo digs deeper).
"""
@views function dig(v, mc, colidx=1)
  len = length(v)
  if len == 1
    return v[1]
  end

  mid = len >>> 1 + 1
  first1idx = searchsortedfirst(v, srchterm; by=l -> l[colidx])

  # If mid and first1idx are equal, and the length is even, that means there's
  # an equal number of 1's and 0's, so we pick the 1's if mc is true and the
  # 0's otherwise. If that's not the case, then if first1idx <= mid there are
  # more 1's than 0's and we take them if mc is true, or take the zeros
  # otherwise.
  takehi = first1idx == mid && iseven(len) ? mc : (first1idx <= mid) == mc
  if takehi
    dig(v[first1idx:len], colidx + 1, mc)
  else
    dig(v[1:first1idx - 1], colidx + 1, mc)
  end
end

oxygen_rating = dig(lines, true)
co2_rating = dig(lines, false)

println("Oxygen: ", oxygen_rating)
println("CO2: ", co2_rating)

oxygen_rating = parse(Int, oxygen_rating; base=2)
co2_rating = parse(Int, co2_rating; base=2)

println("Result: ", oxygen_rating * co2_rating)


text:
using Random

Random.seed!(1)
len=15
xss_base = [randstring(len), rand(1:100, len), [randstring(rand(1:10)) for ii=1:len]]
const xss =  [xss_base; Base.Iterators.flatten.(xss_base)]

@testset "consolidate" begin
    for xs in xss
        for ii in 1:len
            @test consolidate(xs)[ii] == collect(xs)[ii]
        end
    end
end


text:
module Interpol

	include("inplaceadd.jl")

	include("Interpolation.jl")

	include("interpolations/linear.jl")
	include("interpolations/nearest.jl")

	include("Boundaries.jl")

	include("boundary_conditions/constant.jl")

	include("Compound.jl")

end

max_stars_repo_name: JKrehl/Interpol

max_issues_repo_name: JKrehl/Interpol

max_forks_repo_name: JKrehl/Interpol


text:
using Base.Test
using ThreePhotons

function test_random_rotations()
  p1 = [0,1,0]
  p2 = rand(3)#[0,0,1]#/sqrt(2)
  s1 = 0
  s2 = 0
  for i = 1:10000
      rot = random_rotation(3)
      s1 += rot*p1
      s2 += rot*p2
  end
  println(s1/10000)
  println(s2/10000)
end

function test_random_rotations2()
  val = sum([det(random_rotation(3)) for i = 1:10000])
  println(val)
end

#Tests if random_rotation for D=3 samples all rotations uniformly
function test_complettness_of_rotations(random_rotation_generator=random_rotation)

  angles = [get_euler_angles(random_rotation_generator(3)) for l = 1:100000]
  a = Base.hist(convert(Array{Float64},map((x)->x[1], angles)), 100)
  b = Base.hist(convert(Array{Float64},map((x)->x[2], angles)), 100)
  c = Base.hist(convert(Array{Float64},map((x)->x[3], angles)), 100)

  x1 = collect(a[1])
  x2 = collect(b[1])
  x3 = collect(c[1])
  plot(x1[1:length(x1)-1], a[2]/sum(abs, a[2]), label="phi")
  plot(x2[1:length(x2)-1], b[2]/sum(abs, b[2]), label="theta")
  plot(x3[1:length(x3)-1], c[2]/sum(abs, c[2]), label="gamma")
  xlabel("Angle [rad]")
  ylabel("Distribution of Angles")
  legend()
  # [a[1];]

end

function test_random_rotations3(func=random_rotation)
    global veclist = Float64[];
    init = [1,0,0]
    append!(veclist, init)
    num = int(1e3)
    for i = 1:num
        append!(veclist, func(3)*init)
    end
    # global z = zeros(3,num)
    global v = reshape(veclist, 3, num+1)
    scatter3D( v[1,:], v[2,:], v[3,:])
    xlim(-1, 1)
    ylim(-1,1)
    zlim(-1,1)
end

function test_random_rotations4(func=random_rotation)
    global veclist = Float64[];
    init = [1,0]
    append!(veclist, init)
    num = int(1e3)
    for i = 1:num
        append!(veclist, func(2)*init)
    end
    # global z = zeros(3,num)
    global v = reshape(veclist, 2, num+1)
    scatter( v[1,:], v[2,:])
    xlim(-1, 1)
    ylim(-1,1)
    # zlim(-1,1)
end

function test_random_roation_stepper()
    global veclist = Float64[];
    init = [1,0,0]
    num = int(1e4)
    for i = 1:num
        init = random_rotation_step(3, pi/8)*init
        append!(veclist, init)
    end
    # global z = zeros(3,num)
    global v = reshape(veclist, 3, num)
    scatter3D( v[1,:], v[2,:], v[3,:])
    xlim(-1, 1)
    ylim(-1,1)
    zlim(-1,1)
end

"""Tests SON_parametrized by throwing random angles and plotting the point cloud in 3D. This should give a seamlessly but ununiform sampling of the 3D sphere."""
function test_SON_paramtrized()
  r = 0:0.1:pi
  list = []
  for i = 1:1000
    push!(list, [0 0 1] * SON_parametrized(3, [2*pi*rand(), 2*pi*rand(), 2*pi*rand()]))
  end
  x = [list[i][1] for i=1:length(list)]
  y = [list[i][2] for i=1:length(list)]
  z = [list[i][3] for i=1:length(list)]
  scatter3D(x,y,z)
end


text:
"""
Abstract type for defining the strategy for injectimg memory noise.
Implemented strategies:

  * [`NoMemory`](@ref)
  * [`DeltaMemory`](@ref)
  * [`ShiftMemory`](@ref)
"""
abstract type MemoryStrategy end

"""
    NoMemory <: MemoryStrategy
Default strategy for [`MemoryStrategy`](@ref) indicating that
no memory noise will be used.
"""
struct NoMemory <: MemoryStrategy end

"""
    PurgeNegatives <: MemoryStrategy
Purge all negative sign walkers.
"""
struct PurgeNegatives <: MemoryStrategy end

"""
    DeltaMemory(Δ::Int) <: MemoryStrategy
Before updating the shift, memory noise with a memory length of `Δ` is applied,
where `Δ = 1` means no memory noise.

```
r̃ = (pnorm - tnorm)/(dτ*pnorm) + shift
r = r̃ - <r̃>
```
"""
mutable struct DeltaMemory <: MemoryStrategy
    Δ::Int # length of memory noise buffer
    pnorm::Float64 # previous norm
    noiseBuffer::DataStructures.CircularBuffer{Float64} # buffer for memory noise
end
DeltaMemory(Δ::Int) = DeltaMemory(Δ, NaN, DataStructures.CircularBuffer{Float64}(Δ))

"""
    DeltaMemory2(Δ::Int) <: MemoryStrategy
Before updating the shift, memory noise with a memory length of `Δ` is applied,
where `Δ = 1` means no memory noise.

```
r̃ = pnorm - tnorm + shift*dτ*pnorm
r = (r̃ - <r̃>)/(dτ*pnorm)
```
The long-term average of `r` is not guaranteed to be zero.
"""
mutable struct DeltaMemory2 <: MemoryStrategy
    Δ::Int # length of memory noise buffer
    pnorm::Float64 # previous norm
    noiseBuffer::DataStructures.CircularBuffer{Float64} # buffer for memory noise
end
DeltaMemory2(Δ::Int) = DeltaMemory2(Δ, NaN, DataStructures.CircularBuffer{Float64}(Δ))

"""
    DeltaMemory3(Δ::Int, level::Float64) <: MemoryStrategy
Before updating the shift, apply multiplicative memory noise with a
memory length of `Δ` at level `level`,
where `Δ = 1` means no memory noise.

```
r̃ = (pnorm - tnorm)/pnorm + dτ*shift
r = r̃ - <r̃>
w .*= 1 + level*r
```
"""
mutable struct DeltaMemory3 <: MemoryStrategy
    Δ::Int # length of memory noise buffer
    level::Float64 # previous norm
    noiseBuffer::DataStructures.CircularBuffer{Float64} # buffer for memory noise
end
DeltaMemory3(Δ::Int,level::Float64) = DeltaMemory3(Δ, level, DataStructures.CircularBuffer{Float64}(Δ))

"""
    ShiftMemory(Δ::Int) <: MemoryStrategy
Effectively replaces the fluctuating `shift` update procedure for the
coefficient vector by an averaged `shift` over `Δ` timesteps,
where `Δ = 1` means no averaging.
"""
struct ShiftMemory <: MemoryStrategy
    Δ::Int # length of memory noise buffer
    noiseBuffer::DataStructures.CircularBuffer{Float64} # buffer for memory noise
end
ShiftMemory(Δ::Int) = ShiftMemory(Δ, DataStructures.CircularBuffer{Float64}(Δ))

"""
    ProjectedMemory(Δ::Int, projector, pp::Number) <: MemoryStrategy
    ProjectedMemory(Δ::Int, projector, v::AbstractDVec)
Before updating the shift, apply memory noise to minimize the fluctuations
of the overlap of the coefficient vector with `projector`.
Averaging over `Δ` time steps is applied, where `Δ = 1` means no memory noise
is applied. Use `pp` to initialise the value of the projection or pass `v` in
order to initialise the projection with `pp = projector.v`.
```
r̃ = (projector⋅v - projector⋅w)/projector⋅v + dτ*shift
r = r̃ - <r̃>
```
where `v` is the coefficient vector before and `w` after applying a regular
FCIQMC step.
"""
mutable struct ProjectedMemory{D} <: MemoryStrategy
    Δ::Int # length of memory noise buffer
    pp::Float64 # previous projection
    projector::D # projector
    noiseBuffer::DataStructures.CircularBuffer{Float64} # buffer for memory noise
end
ProjectedMemory(Δ::Int, projector, pp::Number) = ProjectedMemory(Δ, pp, projector, DataStructures.CircularBuffer{Float64}(Δ))
function ProjectedMemory(Δ::Int, projector, v::AbstractDVec)
    pp = projector⋅v
    ProjectedMemory(Δ, pp, projector, DataStructures.CircularBuffer{Float64}(Δ))
end

mutable struct ProjectedMemory2{D} <: MemoryStrategy
    Δ::Int # length of memory noise buffer
    pp::Float64 # previous projection
    projector::D # projector
    noiseBuffer::DataStructures.CircularBuffer{Float64} # buffer for memory noise
end
ProjectedMemory2(Δ::Int, projector, pp::Number) = ProjectedMemory2(Δ, pp, projector, DataStructures.CircularBuffer{Float64}(Δ))
function ProjectedMemory2(Δ::Int, projector, v::AbstractDVec)
    pp = projector⋅v
    ProjectedMemory2(Δ, pp, projector, DataStructures.CircularBuffer{Float64}(Δ))
end

mutable struct ProjectedMemory3{D} <: MemoryStrategy
    Δ::Int # length of memory noise buffer
    pp::Float64 # previous projection
    projector::D # projector
    noiseBuffer::DataStructures.CircularBuffer{Float64} # buffer for memory noise
end
ProjectedMemory3(Δ::Int, projector, pp::Number) = ProjectedMemory3(Δ, pp, projector, DataStructures.CircularBuffer{Float64}(Δ))
function ProjectedMemory3(Δ::Int, projector, v::AbstractDVec)
    pp = projector⋅v
    ProjectedMemory3(Δ, pp, projector, DataStructures.CircularBuffer{Float64}(Δ))
end
mutable struct ProjectedMemory4{D} <: MemoryStrategy
    Δ::Int # length of memory noise buffer
    pp::Float64 # previous projection
    projector::D # projector
    noiseBuffer::DataStructures.CircularBuffer{Float64} # buffer for memory noise
end
ProjectedMemory4(Δ::Int, projector, pp::Number) = ProjectedMemory4(Δ, pp, projector, DataStructures.CircularBuffer{Float64}(Δ))
function ProjectedMemory4(Δ::Int, projector, v::AbstractDVec)
    pp = projector⋅v
    ProjectedMemory4(Δ, pp, projector, DataStructures.CircularBuffer{Float64}(Δ))
end


text:
module GCPTests

CredFilename = expanduser("~/secrets/tst-cld.json")

ts = [
  "credential"
  ,"session"
  ,"project"
  ,"bigquery"
]
for t in ts
  include("$(t).jl")
end

end # module GCPTests

max_stars_repo_name: ProfLeao/mais

max_issues_repo_name: ProfLeao/mais

max_forks_repo_name: ProfLeao/mais


text:
#define abstract planetary system
abstract type abstractPlanetarySystem end

#define structure that contains all (actively selected) plantery bodies in the system
# some validation is included to determine no unphysical combinations are selected! (only the Sun and the Moon, for example)
struct PlanetarySystem <: abstractPlanetarySystem
    Bodies::Vector{T} where {T <: abstractCelestialBody}
    StarsSelected::Int
    PlanetsSelected::Int
    MoonsSelected::Int
    InertialFrame_Body::Type{<:abstractCelestialBody}

    #Define constructor and validate entries
    function PlanetarySystem(Bodies::Vector{B}) where {B <: abstractCelestialBody}
        resultOK, StarsSelected, PlanetsSelected, MoonsSelected  = SystemValidation(Bodies)
        if resultOK
            if_id = GetInertialFrame(Bodies, StarsSelected, PlanetsSelected, MoonsSelected)
            return new(Bodies, StarsSelected, PlanetsSelected, MoonsSelected, if_id)
        else
            error("Unphysical system was proposed!")
        end
    end #constructor
end #struct

#Check if the system is physically reasonable
function SystemValidation(Bodies::Vector{T}) where {T <: abstractCelestialBody}
    #If true then more than one planet is selected without the Sun!
    StarSelected = 0
    PlanetSelected = 0
    MoonSelected = 0
    MoonValidated = 0
    result = false
    for body in Bodies
        if typeof(body) <: abstractStar #if currently selected body is a Star
            StarSelected += 1
        elseif typeof(body) <: abstractPlanet  #if currently selected body is a Planet
            PlanetSelected += 1
        elseif typeof(body) <: abstractNaturalSatellite  #if currently selected body is a Moon
            MoonSelected += 1
        end
    end

    #its okay when only the star, one planet, or a moon (with or without its main body) is selected
    if StarSelected == 1 #if star is selected
        if PlanetSelected == 0 && MoonSelected == 0 #only star selected
            return true, StarSelected, PlanetSelected, MoonSelected
        elseif PlanetSelected >= 1 && MoonSelected == 0 #only planets selected
            return true, StarSelected, PlanetSelected, MoonSelected
        elseif PlanetSelected >= 1 && MoonSelected >= 1 #Planets and moons selected -> must check if moon have their main bodies (because sun is present)
            for body in Bodies #loop over all bodies to find all moons
                if typeof(body) <: abstractNaturalSatellite #check if currently selected object is a moon
                    for body2 in Bodies #loop over all bodies to find the planet of the moon
                        if typeof(body2) <: abstractPlanet #check if currently selected body is a planet
                            if parent(body) == parent(body2)  #check if both share the same parent (the common barycenter of the planet-moon system)
                                MoonValidated += 1
                            end #validate
                        end #check if planet
                    end #loop over planets
                end #check if moon
            end #loop over moons
            if MoonSelected == MoonValidated #if the number of selected moons equals the number of validated moons, than OK
                return true, StarSelected, PlanetSelected, MoonSelected
            else
                return false, StarSelected, PlanetSelected, MoonSelected
            end
        else
            return false, StarSelected, PlanetSelected, MoonSelected
        end

    elseif StarSelected == 0 #star not selected
        if PlanetSelected == 0 && MoonSelected == 0 #nothing selected, technically valid
            return true, StarSelected, PlanetSelected, MoonSelected
        elseif PlanetSelected == 1 && MoonSelected == 0 #Single planet selected
            return true, StarSelected, PlanetSelected, MoonSelected
        elseif PlanetSelected >= 2 #Two or more planets selected without star, not valid
            return false, StarSelected, PlanetSelected, MoonSelected
        elseif PlanetSelected == 0 && MoonSelected == 1 #Single Moon Selected, is valid
            return true, StarSelected, PlanetSelected, MoonSelected
        elseif PlanetSelected ==0 && MoonSelected >= 2 #Two or more Moons selected, not valid
            return false, StarSelected, PlanetSelected, MoonSelected
        elseif PlanetSelected == 1 && MoonSelected >= 1 #Can be valid, if the single planet shares the barycenter with all moons (for example, jupiter system)
            for mooncand in Bodies
                if typeof(mooncand) <: abstractNaturalSatellite
                    for planetcand in Bodies
                        if typeof(planetcand) <: abstractPlanet
                            if parent(mooncand) == parent(planetcand)
                                MoonValidated += 1
                            end
                        end
                    end
                end
            end
            #Only valid if number of moons selected equal the number of moons validated
            (MoonSelected == MoonValidated)  ? result = true : result = false
            return result, StarSelected, PlanetSelected, MoonSelected
        else
            return false, StarSelected, PlanetSelected, MoonSelected
        end
    end

    #If the codes makes it to here, the validation has failed
    return false, StarSelected, PlanetSelected, MoonSelected

end #Function TwoPlanetsWithoutSunSelected

#Determine inertial reference frame from the list of Bodies
function GetInertialFrame(Bodies::Vector{T}, StarsSelected::Int ,PlanetsSelected::Int,
        MoonsSelected::Int) where {T <: abstractCelestialBody}
    if StarsSelected == 1
        return SSB
    else
        if PlanetsSelected == 0 && MoonsSelected == 0 #Nothing selected (free space)
            error("Nothing Selected!")
        elseif PlanetsSelected == 1 && MoonsSelected == 0 #Only a single planet
            return Bodies[1].Name
        elseif PlanetsSelected >= 2 #more than one planet without sun, should have been found by validation
            error("Validation has failed! More than one Planet selected without central star")
        elseif PlanetsSelected == 1 && MoonsSelected == 1 #Planet with moon (Validated Already)
            return parent(Bodies[1]) #Both have common barycenter
        elseif PlanetsSelected == 1 && MoonsSelected >= 1 #single planet with multiple moons (Validated already)
            return parent(Bodies[1]) #all have common barycenter
        end
    end
end


text:
# JuliaText TextAnalysis.jl Utility Functions

# The jackknife is a resampling technique especially useful for variance and bias estimation.
#Currently being used for averaging in ROUGE scores in evaluate.jl
# :param scores: List of integers to average
#:type scores: Array{Int64,1}

function jackknife_avg(scores)

    if length(collect(Set(scores))) == 1

        #= In case the elements of the array are all equal=#
        return scores[1]

    else

        #=store the maximum scores
        from the m different sets of m-1 scores.
        such that m is the len(score_list)=#

        average = []

        for i in scores
            # dummy : list a particular combo of m-1 scores
            dummy = [j for j in scores if i != j]
            append!(average, max(dummy...))
            end

        return sum(average)/length(average)

    end
end

# Returns an array of ngram tokens from the ngram doc.
# param ngram_doc : Output from the function ngrams(StringDocument(text), n)
function listify_ngrams(ngram_doc)
    flattened = []
    for item in ngram_doc
        for i in 1:item[2]
            push!(flattened, item[1])
        end
    end
    return flattened
end

#This function returns the longest common subsequence
#of two strings using the dynamic programming algorithm.
# param X : first string in tokenized form
#type (X) : Array{SubString{String},1}
# param Y : second string in tokenized form
#type (Y) : Array{SubString{String},1}
# param weighted : Weighted LCS is done if weighted is True (default)
#type (weighted) : Boolean
# param return_string : Function returns weighted LCS length when set to False (default).
#                       Function returns longest common substring when set to True.
#type (return_string) : Boolean
# param f: weighting function. The weighting function f must have the property
#          that f(x+y) > f(x) + f(y) for any positive integers x and y.
#type (f) : generic function which takes a float as an input and returns a float.

function weighted_lcs(X, Y, weighted = true, return_string = false, f = sqrt)

    m, n = length(X), length(Y)
    c_table = [zeros(n+1) for i in 1:m+1]
    w_table = [zeros(n+1) for i in 1:m+1]

    for i in 1:(m+1)

        for j in 1:(n+1)

            if i == 1 || j == 1
                continue

            elseif X[i-1] == Y[j-1]

                k = w_table[i-1][j-1]
                if weighted == true
                    increment = (f(k+1)) - (f(k))
                else
                    increment = 1
                end
                c_table[i][j] = c_table[i-1][j-1] + increment
                w_table[i][j] = k + 1

            else

                if c_table[i-1][j] > c_table[i][j-1]
                    c_table[i][j] = c_table[i-1][j]
                    w_table[i][j] = 0  # no match at i,j
                else
                    c_table[i][j] = c_table[i][j-1]
                    w_table[i][j] = 0  # no match at i,j
                end

            end

        end

    end

    lcs_length = (c_table[m+1][n+1])
    if return_string == false
        return lcs_length
    end

    if weighted == true
        lcs_length = c_table[m][n]^(2)
    end

    lcs_length = round(Int64, lcs_length)
    lcs_length = convert(Int64, lcs_length)
    lcs = ["" for i in  1:(lcs_length+1)]
    lcs[lcs_length+1] = ""
    i = m + 1
    j = n + 1

    while i>1 && j>1
        if X[i-1] == Y[j-1]
            lcs[lcs_length+1] = X[i-1]
            i -= 1
            j -= 1
            lcs_length -= 1

        elseif c_table[i-1][j] > c_table[i][j-1]
            i -= 1
        else
            j -= 1
        end
    end

    return (join(lcs, " "))  # the lcs string
end

# F-measure based on WLCS
# param beta : user defined parameter
#type (beta) : float

# param r_lcs : recall factor
#type (r_lcs) : float

# param p_lcs : precision factor
#type (p_lcs) : float

#score : f measure score between a candidate
#    	    and a reference

function fmeasure_lcs(RLCS, PLCS, beta=1)

    try
        return ((1+beta^2)*RLCS*PLCS)/(RLCS+(beta^2)*PLCS)
    catch ex
        if ex isa DivideError
            return 0
        else
            rethrow(ex)
        end
    end
end

# Calculates Word CoOccurence matrix.
# The input to the function can be a single string with a window size
# based on which the co-occurence will be evaluated. A list of strings
# of variable window size can also be taken as input.
# Returns a DataFrame which can be used to calculate the CoOccurence of a pair
# of words easily.

# param inputDoc : string or array of strings
# param integer  : sliding window size. The number of words before and after
#                  a token to be considered in context.
# param stripStopwords : Whether to remove stopwords or not.

# output : A NamedArray arr. arr[Word1, Word2] returns the co-occurrence value of
#          a pair of words. arr[Word, :] returns cooccurence value of Word with all
#          other words

function word_cooccurrence_matrix(inputDoc, window = 2, stripStopwords = false)
    if typeof(inputDoc) == String
        sd = StringDocument(lowercase(inputDoc))
        prepare!(sd, strip_punctuation)
        if stripStopwords == true
            prepare!(sd, strip_stopwords)
            sd = StringDocument(join([i for i in split(text(sd), " ") if i!= ""], " "))
        end
        sd = split(text(sd), " ")
        if length(sd) < 2*window + 1
                        throw(ArgumentError(string("Window Size is greater than the length of string provided\nPlease make sure that the value of 2*window+1 is lesser than the length of string provided.")))
        end
        candidate_list = []
        for item in enumerate(sd)
            if item[1] < length(sd) - 2 * window + 1
                temp = []
                for i in range(1 , 2 * window + 1)
                    push!(temp,sd[item[1] + i - 1])
                end
                push!(candidate_list,temp)
            end
        end
    else
        if typeof(inputDoc) == Array{String,1}

            if window != 0
                candidate_list = []
                for doc in enumerate(inputDoc)
                    sd = StringDocument(lowercase(doc[2]))
                    prepare!(sd, strip_punctuation)
                    if stripStopwords == true
                        prepare!(sd, strip_stopwords)
                        sd = StringDocument(join([i for i in split(text(sd), " ") if i!= ""], " "))
                    end
                    sd = split(text(sd), " ")
                    if length(sd) < 2*window + 1
                        throw(ArgumentError(string("Window Size is greater than the length of document provided at index ", doc[1], "\nPlease make sure that the value of 2*window+1 is lesser than the minimum length of documents provided.")))
                    end
                    for item in enumerate(sd)
                        if item[1] < length(sd) - 2 * window + 1
                            temp = []
                            for i in range(1 , 2 * window + 1)
                                push!(temp,sd[item[1] + i - 1])
                            end
                            push!(candidate_list,temp)
                        end
                    end
                end
            else
                candidate_list = []
                for doc in inputDoc
                    sd = StringDocument(lowercase(doc))
                    prepare!(sd, strip_punctuation)
                    if stripStopwords == true
                        prepare!(sd, strip_stopwords)
                        sd = StringDocument(join([i for i in split(text(sd), " ") if i!= ""], " "))
                    end
                    sd = split(text(sd), " ")
                    push!(candidate_list, sd)
                end
            end

        else
            throw(ArgumentError("Argument must be a String or a List of Strings i.e Array{String,1}"))
        end
    end
    wordlist = sort(unique(split(join([join(i, " ") for i in candidate_list], " "), " ")))
    word_matrix = []
    for i in wordlist
        temp = [0 for i in range(1, length(wordlist))]
        matches = Dict()
        if window == 0
            temp = [0 for i in range(1, length(wordlist))]
            matches = Dict()
            for keywords in candidate_list
                if i in keywords
                    for token in keywords
                        if token in keys(matches)
                            matches[token] += 1
                        else
                            matches[token] = 1
                        end
                    end
                end
            end
        else
            mid = window+1
            for keywords in candidate_list
                if i == keywords[mid]
                    for token in keywords
                        if token in keys(matches)
                            matches[token] += 1
                        else
                            matches[token] = 1
                        end
                    end
                else
                    if candidate_list[length(candidate_list)] == keywords
                        for keyword in enumerate(keywords)
                            if keyword[2] == i && keyword[1] > mid
                                for token in keywords[keyword[1]-mid+1:length(keywords)]
                                    if token in keys(matches)
                                        matches[token] += 1
                                    else
                                        matches[token] = 1
                                    end
                                end
                            end
                        end
                    end
                    if candidate_list[1] == keywords
                        for keyword in enumerate(keywords)
                            if keyword[2] == i && keyword[1] < mid
                                for token in keywords[1:keyword[1]+mid-1]
                                    if token in keys(matches)
                                        matches[token] += 1
                                    else
                                        matches[token] = 1
                                    end
                                end
                            end
                        end
                    end
                end
            end
        end
        for key in keys(matches)
            indexes = [i[1] for i in enumerate(wordlist) if i[2] == key]
            for index in indexes
                temp[index]  = matches[wordlist[index]]
            end
        end
        if length(word_matrix) == 0
            word_matrix = temp
        else
            word_matrix = hcat(word_matrix, temp)
        end
    end

    wordlist = convert(Array{String,1}, wordlist)
    coom = NamedArray(word_matrix, (wordlist, wordlist), ("Rows", "Cols"))

    return(coom)
end