text:
using IntervalArithmetic
include("utils.jl")
import ..AbstractUncertainIndexValueDataset

export
sequence_exists,
sequence_exists!

"""
    sequence_exists(x, c::SequentialSamplingConstraint)

Does a point-by-point sequence through the uncertain dataset `x` exist that satisfies the criteria `c`?

If `x` is an `UncertainIndexValueDataset`, then check for a sequence through the indices only.

Before the check is performed, the distributions in `x` are truncated to the quantiles provided
by `c` to ensure they have finite supports.

## Example

```julia
# Create a set of time indices
# We construct this is such a way that we *know* an increasing sequence exists.
t = [UncertainValue(Normal, i, 2) for i in 1:N];
sequence_exists(t, StrictlyIncreasing(StartToEnd()))
```
"""
function sequence_exists end

# If data has uncertainties both in indices and values, check only for indices.
sequence_exists(x::AbstractUncertainIndexValueDataset, c) = sequence_exists(x.indices, c)

function sequence_exists(x, c::SequentialSamplingConstraint)
    lqs, uqs = get_quantiles(x, c)

    return sequence_exists(lqs, uqs, c), lqs, uqs
end

function sequence_exists!(lqs, uqs, x::AbstractUncertainValueDataset, c)
    get_quantiles!(lqs, uqs, x, c)

    return sequence_exists(lqs, uqs, c)
end

sequence_exists(x::AbstractUncertainIndexValueDataset, c::StrictlyIncreasing{StartToEnd}) =
    sequence_exists(x.indices, c)

###########################
# Concrete implementations
###########################
"""
    sequence_exists(udata::AbstractUncertainValueDataset, c::StrictlyDecreasing{StartToEnd})

Does a strictly decreasing sequence through the dataset exist? I.e,  check that a strictly
decreasing sequence can be found after first
constraining each distribution to the provided quantile range (this
is necessary because some distributions may have infinite support).
"""
function sequence_exists(lqs, uqs, c::StrictlyIncreasing{StartToEnd})
    L = length(lqs)
    if any(lqs .> uqs) # ties are allowed, because we have `CertainValue`s
        error("Not all `lqs[i]` are lower than uqs[i]. Quantile calculations are not meaningful.")
        return false
    end

    for i = 1:L-1
        if lqs[i] >= minimum(uqs[i+1:end])
            return false
        end
    end
    return true
end

function sequence_exists(lqs, uqs, c::StrictlyDecreasing{StartToEnd})
    L = length(lqs)
    if any(lqs .> uqs) # ties are allowed, because we have `CertainValue`s
        error("Not all `lqs[i]` are lower than uqs[i]. Quantile calculations are not meaningful.")
        return false
    end

    for i = 1:L-1
        if uqs[i] < maximum(lqs[i+1:end])
            return false
        end
    end
    return true
end


text: run(`ls`)

ext: jl

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
for (elty, relty, ext) in ((:Float32, :Float32, :s),
                           (:Float64, :Float64, :d),
                           (:ComplexF32, :Float32, :c),
                           (:ComplexF64, :Float64, :z))

    for (mat, sym) in ((:Matrix, "_"),
                       (:DistMatrix, "Dist_"))

        @eval begin

            function gemm(orientationOfA::Orientation, orientationOfB::Orientation, α::$elty, A::$mat{$elty}, B::$mat{$elty}, β::$elty, C::$mat{$elty})
                ElError(ccall(($(string("ElGemm", sym, ext)), libEl), Cuint,
                    (Orientation, Orientation, $elty, Ptr{Cvoid}, Ptr{Cvoid}, $elty, Ptr{Cvoid}),
                     orientationOfA, orientationOfB, α, A.obj, B.obj, β, C.obj))
                return C
            end

            function trsm(side::LeftOrRight, uplo::UpperOrLower, orientation::Orientation, diag::UnitOrNonUnit, α::$elty, A::$mat{$elty}, B::$mat{$elty})
                ElError(ccall(($(string("ElTrsm", sym, ext)), libEl), Cuint,
                    (LeftOrRight, UpperOrLower, Orientation, UnitOrNonUnit,
                     $elty, Ptr{Cvoid}, Ptr{Cvoid}),
                     side, uplo, orientation, diag,
                     α, A.obj, B.obj))
                return B
            end
        end
    end
end


text:
module TokamakSmallCylindrical

    import ElectromagneticFields.AxisymmetricTokamakCylindrical

    export charged_particle_3d_pode, charged_particle_3d_iode,
           hamiltonian, toroidal_momentum

    AxisymmetricTokamakCylindrical.@code() # inject magnetic field code

    include("charged_particle_3d_canonical.jl")

    const xᵢ = [1.05, 0.0, 0.0]
    const qᵢ = from_cartesian(0, xᵢ)
    const vᵢ = DF̄(0, xᵢ) * [2.1E-3, 4.3E-4, 0.0]
    const pᵢ = charged_particle_3d_pᵢ(qᵢ, vᵢ)

end


text:
function get_init_values_for_comparison(sim::Simulation)
    bus_size = PSID.get_bus_count(sim.simulation_inputs)
    system = PSID.get_system(sim)
    V_R = Vector{Float64}(undef, bus_size)
    V_I = Vector{Float64}(undef, bus_size)
    Vm = Vector{Float64}(undef, bus_size)
    θ = Vector{Float64}(undef, bus_size)
    for bus in PSY.get_components(PSY.Bus, system)
        bus_n = PSY.get_number(bus)
        bus_ix = PSID.get_lookup(sim.simulation_inputs)[bus_n]
        V_R[bus_ix] = sim.x0_init[bus_ix]
        V_I[bus_ix] = sim.x0_init[bus_ix + bus_size]
        Vm[bus_ix] = sqrt(V_R[bus_ix]^2 + V_I[bus_ix]^2)
        θ[bus_ix] = angle(V_R[bus_ix] + V_I[bus_ix] * 1im)
    end
    results =
        Dict{String, Vector{Float64}}("V_R" => V_R, "V_I" => V_I, "Vm" => Vm, "θ" => θ)
    for device in PSY.get_components(PSY.DynamicInjection, system)
        states = PSY.get_states(device)
        name = PSY.get_name(device)
        global_index = PSID.get_global_index(sim.simulation_inputs)[name]
        x0_device = Vector{Float64}(undef, length(states))
        for (i, s) in enumerate(states)
            x0_device[i] = sim.x0_init[global_index[s]]
        end
        results[name] = x0_device
    end
    dyn_branches = PSY.get_components(PSY.DynamicBranch, system)
    if !isempty(dyn_branches)
        for br in dyn_branches
            states = PSY.get_states(br)
            name = PSY.get_name(br)
            global_index = PSID.get_global_index(sim.simulation_inputs)[name]
            x0_br = Vector{Float64}(undef, length(states))
            for (i, s) in enumerate(states)
                x0_br[i] = sim.x0_init[global_index[s]]
            end
            printed_name = "Line " * name
            results[printed_name] = x0_br
        end
    end
    return results
end

function clean_extra_timestep!(t::Vector{Float64}, δ::Vector{Float64})
    idx = unique(i -> t[i], 1:length(t))
    return t[idx], δ[idx]
end

function get_csv_delta(str::AbstractString)
    M = readdlm(str, ',')
    return clean_extra_timestep!(M[:, 1], M[:, 2])
end

function get_csv_data(str::AbstractString)
    M = readdlm(str, ',')
    return M
end


text:
export extend, NfToNfMor, automorphisms, automorphism_group

struct NfMorSet{T}
  field::T
end

function show(io::IO, S::NfMorSet{T}) where {T}
  print(io, "Set of automorphisms of ", S.field)
end

parent(f::NfToNfMor) = NfMorSet(domain(f))

function parent(f::NumFieldMor)
  @assert domain(f) == codomain(f)
  return NfMorSet(domain(f))
end

function image(f::NumFieldMor, a::FacElem{S, T}) where {S <: NumFieldElem, T <: NumField}
  D = Dict{elem_type(codomain(f)), fmpz}(f(b) => e for (b, e) in a)
  return FacElem(D)
end

################################################################################
#
#  Some basic properties of NfToNfMor
#
################################################################################

is_injective(m::NumFieldMor) = true

is_surjective(m::NumFieldMor) = absolute_degree(domain(m) == absolute_degree(codomain(m)))

is_bijective(m::NumFieldMor) = is_surjective(m)

################################################################################
#
#  Generic groups to set of homomorphisms
#
################################################################################

mutable struct GrpGenToNfMorSet{S, T} <: Map{GrpGen, NfMorSet{T}, HeckeMap, GrpGenToNfMorSet{S, T}}
  G::GrpGen
  aut::Vector{S}
  header::MapHeader{GrpGen, NfMorSet{T}}

  function GrpGenToNfMorSet(aut::Vector{S}, G::GrpGen, s::NfMorSet{T}) where {S, T}
    z = new{S, T}()
    z.header = MapHeader(G, s)
    z.aut = aut
    z.G = G
    return z
  end
end

function GrpGenToNfMorSet(G::GrpGen, K::NumField)
  return GrpGenToNfMorSet(automorphisms(K), G, NfMorSet(K))
end

function GrpGenToNfMorSet(G::GrpGen, aut::Vector{S}, K::NumField) where S <: NumFieldMor
  return GrpGenToNfMorSet(aut, G, NfMorSet(K))
end

function image(f::GrpGenToNfMorSet, g::GrpGenElem)
  @assert parent(g) == f.G
  K = codomain(f).field
  return f.aut[g[]]
end

function (f::GrpGenToNfMorSet)(g::GrpGenElem)
  return image(f, g)
end

function preimage(f::GrpGenToNfMorSet{S, T}, a::S) where {S, T}
  K = codomain(f).field
  aut = f.aut
  for i in 1:length(aut)
    if a == aut[i]
      return domain(f)[i]
    end
  end
  error("something wrong")
end

function evaluate(f::fmpq_poly, a::nf_elem)
  #Base.show_backtrace(stdout, Base.stacktrace())
  R = parent(a)
  if iszero(f)
    return zero(R)
  end
  if a == gen(R) && parent(f) == parent(parent(a).pol)
    return R(f)
  end
  l = length(f) - 1
  s = R(coeff(f, l))
  for i in l-1:-1:0
    #s = s*a + R(coeff(f, i))
    mul!(s, s, a)
    add!(s, s, coeff(f, i))
  end
  return s
end

Base.copy(f::NfToNfMor) = f

################################################################################
#
#  is normal
#
################################################################################

@doc Markdown.doc"""
     is_normal(K::AnticNumberField) -> Bool

Returns true if $K$ is a normal extension of $\mathbb Q$, false otherwise.
"""
function is_normal(K::AnticNumberField)
  #Before computing the automorphisms, I split a few primes and check if the
  #splitting behaviour is fine
  c = get_attribute(K, :is_normal)
  if c isa Bool
    return c::Bool
  end
  fl = is_normal_easy(K)
  if !fl
    return false
  end
  if length(automorphisms(K, copy = false)) != degree(K)
    set_attribute!(K, :is_normal => false)
    return false
  else
    set_attribute!(K, :is_normal => true)
    return true
  end
end

function is_normal_easy(K::AnticNumberField)
  E = any_order(K)
  p = 1000
  ind = 0
  while ind < 15
    p = next_prime(p)
    F = GF(p, cached = false)
    Fx = PolynomialRing(F, cached = false)[1]
    fF = Fx(K.pol)
    if degree(fF) != degree(K) || iszero(discriminant(fF))
      continue
    end
    ind += 1
    dt = prime_decomposition_type(E, p)
    if !divisible(degree(K), length(dt))
      set_attribute!(K, :is_normal => false)
      return false
    end
    f = dt[1][1]
    for i = 2:length(dt)
      if f != dt[i][1]
        set_attribute!(K, :is_normal => false)
        return false
      end
    end
  end
  return true
end

is_normal(K::NumField) = length(automorphisms(K)) == degree(K)

################################################################################
#
#  IsCMfield
#
################################################################################
@doc Markdown.doc"""
    is_cm_field(K::AnticNumberField) -> Bool, NfToNfMor

Given a number field $K$, this function returns true and the complex conjugation
if the field is CM, false and the identity otherwise.
"""
function is_cm_field(K::NumField)
  c = get_attribute(K, :cm_field)
  if c !== nothing
    return true, c::morphism_type(K)
  end
  if isodd(degree(K)) || !is_totally_complex(K)
    return false, id_hom(K)
  end
  if is_automorphisms_known(K)
    auts = automorphisms(K, copy = false)
    return _find_complex_conj(auts)
  end
  if !is_cm_field_easy(K)
    return false, id_hom(K)
  end
  auts = _automorphisms_center(K)
  return _find_complex_conj(auts)
end

function _automorphisms_center(K::NumField)
  return automorphisms(K)
end

function is_cm_field_known(K::NumField)
  c = get_attribute(K, :cm_field)
  return c !== nothing
end

function _find_complex_conj(auts::Vector{NfToNfMor})
  K = domain(auts[1])
  for x in auts
    if !is_involution(x)
      continue
    end
    if is_complex_conjugation(x)
      set_attribute!(K, :cm_field => x)
      return true, x
    end
  end
  return false, id_hom(K)
end

function is_cm_field_easy(K::AnticNumberField)
  E = any_order(K)
  if is_maximal_order_known(K)
    E = maximal_order(K)
  end
  n = degree(E)
  g = zero_matrix(FlintZZ, n, n)
  B = basis(E, nf(E))
  prec = 32
  imgs = Vector{Vector{arb}}(undef, n)
  for i = 1:n
    imgs[i] = minkowski_map(B[i], prec)
  end
  i = 1
  t = arb()
  while i <= n
    j = i
    while j <= n
      el = imgs[i][1]*imgs[j][1]
      for k = 2:n
        mul!(t, imgs[i][k], imgs[j][k])
        add!(el, el, t)
      end
      if radius(el) > 1//16
        prec *= 2
        for k = i:n
          imgs[k] = minkowski_map(B[k], prec)
        end
        continue
      end
      fl, r = unique_integer(el)
      if !fl
        return false
      end
      j += 1
    end
    i += 1
  end
  return true

end

################################################################################
#
#  Induced image
#
################################################################################

function _evaluate_mod(f::fmpq_poly, a::nf_elem, d::fmpz)
  #Base.show_backtrace(stdout, Base.stacktrace())
  R = parent(a)
  if iszero(f)
    return zero(R)
  end
  l = length(f) - 1
  s = R(coeff(f, l))
  for i in l-1:-1:0
    #s = s*a + R(coeff(f, i))
    mul!(s, s, a)
    # TODO (easy): Once fmpq_poly_add_fmpq is improved in flint, remove the R(..)
    add!(s, s, R(coeff(f, i)))
    s = mod(s, d)
  end
  return s
end

(f::NfToNfMor)(x::NfOrdIdl) = induce_image(f, x)

function induce_image(f::NfToNfMor, x::NfOrdIdl)
  K = domain(f)
  if K != codomain(f)
    OK = maximal_order(codomain(f))
    @assert is_maximal(order(x))
    assure_2_normal(x)
    I = ideal(OK, x.gen_one, OK(f(x.gen_two.elem_in_nf)))
    I.gens_normal = x.gens_normal
    return I
  end

  if isone(x)
    return x
  end

  OK = order(x)
  K = nf(OK)
 if has_2_elem(x) && is_maximal_known(OK) && is_maximal(OK)
    int_in_ideal = x.gen_one
    if has_minimum(x)
      int_in_ideal = minimum(x, copy = false)
    elseif has_norm(x)
      int_in_ideal = norm(x, copy = false)
    end
    if is_coprime(index(OK, copy = false), int_in_ideal) && fits(Int, int_in_ideal^2)
    #The conjugate of the prime will still be a prime over the minimum
    #I just need to apply the automorphism modularly
      return induce_image_easy(f, x)
    end
  end
  I = ideal(OK)
  if isdefined(x, :gen_two)
    new_gen_two = f(K(x.gen_two))
    if has_minimum(x)
      new_gen_two = mod(new_gen_two, minimum(x, copy = false)^2)
    end
    if is_maximal_known(OK) && is_maximal(OK)
      I.gen_two = OK(new_gen_two, false)
    else
      I.gen_two = OK(new_gen_two)
    end
  end
  if isdefined(x, :princ_gen)
     if is_maximal_known(OK) && is_maximal(OK)
      I.princ_gen = OK(f(K(x.princ_gen)), false)
    else
      I.princ_gen = OK(f(K(x.princ_gen)))
    end
  end
  for i in [:gen_one, :is_prime, :gens_normal, :gens_weakly_normal, :is_principal,
          :iszero, :minimum, :norm, :splitting_type]
    if isdefined(x, i)
      setfield!(I, i, getfield(x, i))
    end
  end
  if !has_2_elem(I)
    #I need to translate the basis matrix
    bb = Vector{NfOrdElem}(undef, degree(K))
    B = basis(x, copy = false)
    for i = 1:length(bb)
      bb[i] = OK(f(K(B[i])))
    end
    I.basis = bb
    M = zero_matrix(FlintZZ, degree(K), degree(K))
    for i = 1:degree(K)
      el = coordinates(I.basis[i])
      for j = 1:degree(K)
        M[i, j] = el[j]
      end
    end
    I.basis_matrix = M
  end
  return I
end

function induce_image_easy(f::NfToNfMor, P::NfOrdIdl)
  OK = order(P)
  K = nf(OK)
  R = ResidueRing(FlintZZ, Int(minimum(P, copy = false))^2, cached = false)
  Rx = PolynomialRing(R, "t", cached = false)[1]
  fmod = Rx(K.pol)
  prim_img = Rx(image_primitive_element(f))
  gen_two = Rx(P.gen_two.elem_in_nf)
  img = compose_mod(gen_two, prim_img, fmod)
  new_gen = OK(lift(K, img), false)
  res = ideal(OK, minimum(P), new_gen)
  if isdefined(P, :princ_gen)
    res.princ_gen = OK(f(K(P.princ_gen)), false)
  end
  for i in [:is_prime, :gens_normal, :gens_weakly_normal, :is_principal,
          :minimum, :norm, :splitting_type]
    if isdefined(P, i)
      setfield!(res, i, getfield(P, i))
    end
  end
  return res
end

################################################################################
#
#  Maps to algebras
#
################################################################################

# Embedding of a number field into an algebra over Q.
mutable struct NfAbsToAbsAlgAssMor{S} <: Map{AnticNumberField, S, HeckeMap, NfAbsToAbsAlgAssMor}
  header::MapHeader{AnticNumberField, S}
  mat::fmpq_mat
  t::fmpq_mat

  function NfAbsToAbsAlgAssMor{S}(K::AnticNumberField, A::S, M::fmpq_mat) where { S <: AbsAlgAss{fmpq} }
    z = new{S}()
    z.mat = M
    z.t = zero_matrix(FlintQQ, 1, degree(K))

    function _image(x::nf_elem)
      for i = 1:degree(K)
        z.t[1, i] = coeff(x, i - 1)
      end
      s = z.t*z.mat
      return A([ s[1, i] for i = 1:dim(A) ])
    end

    z.header = MapHeader{AnticNumberField, S}(K, A, _image)
    return z
  end
end

function NfAbsToAbsAlgAssMor(K::AnticNumberField, A::S, M::fmpq_mat) where { S <: AbsAlgAss{fmpq} }
  return NfAbsToAbsAlgAssMor{S}(K, A, M)
end

function haspreimage(m::NfAbsToAbsAlgAssMor, a::AbsAlgAssElem)
  A = parent(a)
  t = matrix(FlintQQ, 1, dim(A), coefficients(a))
  b, p = can_solve_with_solution(m.mat, t, side = :left)
  if b
    return true, domain(m)([ p[1, i] for i = 1:nrows(m.mat) ])
  else
    return false, zero(domain(m))
  end
end

################################################################################
#
#  Order of an automorphism in the automorphisms group
#
################################################################################
@doc Markdown.doc"""
    is_involution(f::NfToNfMor) -> Bool

Returns true if $f$ is an involution, i.e. if $f^2$ is the identity, false otherwise.
"""
function is_involution(f::NfToNfMor)
  K = domain(f)
  @assert K == codomain(f)
  if image_primitive_element(f) == gen(K)
    return false
  end
  p = 2
  R = ResidueRing(FlintZZ, p, cached = false)
  Rt = PolynomialRing(R, "t", cached = false)[1]
  fmod = Rt(K.pol)
  while iszero(discriminant(fmod))
    p = next_prime(p)
    R = ResidueRing(FlintZZ, p, cached = false)
    Rt = PolynomialRing(R, "t", cached = false)[1]
    fmod = Rt(K.pol)
  end
  i = 2
  ap = Rt(image_primitive_element(f))
  fp = compose_mod(ap, ap, fmod)
  return fp == gen(Rt)
end

#@doc Markdown.doc"""
#    _order(f::NfToNfMor) -> Int
#
#If $f$ is an automorphism of a field $K$, it returns the order of $f$ in the automorphism group of $K$.
#"""
function _order(f::NfToNfMor)
  K = domain(f)
  @req K === codomain(f) "The morphism must be an automorphism"
  if image_primitive_element(f) == gen(K)
    return 1
  end
  p = 2
  R = ResidueRing(FlintZZ, p, cached = false)
  Rt = PolynomialRing(R, "t", cached = false)[1]
  fmod = Rt(K.pol)
  while iszero(discriminant(fmod))
    p = next_prime(p)
    R = ResidueRing(FlintZZ, p, cached = false)
    Rt = PolynomialRing(R, "t", cached = false)[1]
    fmod = Rt(K.pol)
  end
  i = 2
  ap = Rt(image_primitive_element(f))
  fp = compose_mod(ap, ap, fmod)
  while fp != gen(Rt)
    i += 1
    fp = compose_mod(ap, fp, fmod)
  end
  return i
end

function small_generating_set(G::Vector{NfToNfMor})

  if length(G) == 1
    return G
  end

  firsttry = 10
  secondtry = 20
  thirdtry = 30

	K = domain(G[1])
	p = 2
  R = GF(p, cached = false)
	Rx = PolynomialRing(R, "x", cached = false)[1]
	while iszero(discriminant(Rx(K.pol)))
		p = next_prime(p)
	  R = GF(p, cached = false)
		Rx = PolynomialRing(R, "x", cached = false)[1]
	end

  given_gens = gfp_poly[Rx(image_primitive_element(x)) for x in G]
	orderG = length(closure(given_gens, (x, y) -> Hecke.compose_mod(x, y, Rx(K.pol)), gen(Rx)))
  # First try one element

  for i in 1:firsttry
    trygen = _non_trivial_randelem(G, id_hom(K))
    if length(closure(gfp_poly[Rx(image_primitive_element(trygen))], (x, y) -> Hecke.compose_mod(x, y, Rx(K.pol)), gen(Rx))) == orderG
      return NfToNfMor[trygen]
    end
  end

  for i in 1:secondtry
    gens = NfToNfMor[_non_trivial_randelem(G, id_hom(K)) for i in 1:2]
    gens_mod = gfp_poly[Rx(image_primitive_element(x)) for x in gens]
    if length(closure(gens_mod, (x, y) -> Hecke.compose_mod(x, y, Rx(K.pol)), gen(Rx))) == orderG
      return unique(gens)
    end
  end

  for i in 1:thirdtry
    gens = NfToNfMor[_non_trivial_randelem(G, id_hom(K)) for i in 1:3]
    gens_mod = gfp_poly[Rx(image_primitive_element(x)) for x in gens]
    if length(closure(gens_mod, (x, y) -> Hecke.compose_mod(x, y, Rx(K.pol)), gen(Rx))) == orderG
      return unique(gens)
    end
  end

  # Now use that unconditionally log_2(|G|) elements generate G

  b = ceil(Int, log(2, orderG))
  @assert orderG <= 2^b

  j = 0
  while true
    if j > 2^20
      error("Something wrong with generator search")
    end
    j = j + 1
    gens = NfToNfMor[_non_trivial_randelem(G, id_hom(K)) for i in 1:b]
    gens_mod = gfp_poly[Rx(image_primitive_element(x)) for x in gens]
    if length(closure(gens_mod, (x, y) -> Hecke.compose_mod(x, y, Rx(K.pol)), gen(Rx))) == orderG
      return unique(gens)
    end
  end
end

function _order(G::Vector{NfToNfMor})
  K = domain(G[1])
	p = 2
  R = GF(p, cached = false)
	Rx = PolynomialRing(R, "x", cached = false)[1]
	while iszero(discriminant(Rx(K.pol)))
		p = next_prime(p)
	  R = GF(p, cached = false)
		Rx = PolynomialRing(R, "x", cached = false)[1]
	end
  given_gens = gfp_poly[Rx(image_primitive_element(x)) for x in G]
	return length(closure(given_gens, (x, y) -> Hecke.compose_mod(x, y, Rx(K.pol)), gen(Rx)))
end

################################################################################
#
#  Frobenius automorphism
#
################################################################################

function frobenius_automorphism(P::NfOrdIdl)
  @assert is_prime(P)
  OK = order(P)
  K = nf(OK)
  @assert is_maximal_known_and_maximal(OK)
  @assert ramification_index(P) == 1
  @assert is_normal(K)
  K = nf(OK)
  auts = decomposition_group(P)
  F, mF = ResidueField(OK, P)
  p = minimum(P, copy = false)
  genF = elem_in_nf(mF\gen(F))
  powgen = gen(F)^p
  for i = 1:length(auts)
    img = auts[i](genF)
    if mF(OK(img, false)) == powgen
      return auts[i]
    end
  end
  error("Something went wrong")
end


text:
import Base: isidentifier, is_id_start_char, is_id_char

const RESERVED_WORDS = Set(["begin", "while", "if", "for", "try",
    "return", "break", "continue", "function", "macro", "quote", "let",
    "local", "global", "const", "type",
    "immutable", "do", "module", "baremodule", "using", "import", "struct",
    "export", "importall", "end", "else", "elseif", "catch", "finally"])

if VERSION < v"0.7.0-DEV.3220"
    push!(RESERVED_WORDS, "bitstype", "typealias", "abstract")
end

function identifier(s::AbstractString)
    s = Unicode.normalize(s)
    if !isidentifier(s)
        s = makeidentifier(s)
    end
    Symbol(in(s, RESERVED_WORDS) ? "_"*s : s)
end

function makeidentifier(s::AbstractString)
    i = start(s)
    done(s, i) && return "x"

    res = IOBuffer(zeros(UInt8, sizeof(s)+1), write=true)

    (c, i) = next(s, i)
    under = if is_id_start_char(c)
        write(res, c)
        c == '_'
    elseif is_id_char(c)
        write(res, 'x', c)
        false
    else
        write(res, '_')
        true
    end

    while !done(s, i)
        (c, i) = next(s, i)
        if c != '_' && is_id_char(c)
            write(res, c)
            under = false
        elseif !under
            write(res, '_')
            under = true
        end
    end

    return String(take!(res))
end

function make_unique(names::Vector{Symbol}; makeunique::Bool=false)
    seen = Set{Symbol}()
    names = copy(names)
    dups = Int[]
    for i in 1:length(names)
        name = names[i]
        in(name, seen) ? push!(dups, i) : push!(seen, name)
    end

    if length(dups) > 0
        if !makeunique
            Base.depwarn("Duplicate variable names are deprecated: pass makeunique=true to add a suffix automatically.", :make_unique)
            # TODO: uncomment the lines below after deprecation period
            # msg = """Duplicate variable names: $(u[dups]).
            #          Pass makeunique=true to make them unique using a suffix automatically."""
            # throw(ArgumentError(msg))
        end
    end

    for i in dups
        nm = names[i]
        k = 1
        while true
            newnm = Symbol("$(nm)_$k")
            if !in(newnm, seen)
                names[i] = newnm
                push!(seen, newnm)
                break
            end
            k += 1
        end
    end

    return names
end

"""
    gennames(n::Integer)

Generate standardized names for columns of a DataFrame. The first name will be `:x1`, the
second `:x2`, etc.
"""
function gennames(n::Integer)
    res = Array{Symbol}(undef, n)
    for i in 1:n
        res[i] = Symbol(@sprintf "x%d" i)
    end
    return res
end

"""
    countmissing(a::AbstractArray)

Count the number of `missing` values in an array.
"""
function countmissing(a::AbstractArray)
    res = 0
    for x in a
        res += ismissing(x)
    end
    return res
end

function countmissing(a::CategoricalArray)
    res = 0
    for x in a.refs
        res += x == 0
    end
    return res
end

# Gets the name of a function. Used in groupeDataFrame/grouping.jl
function _fnames(fs::Vector{T}) where T<:Function
    λcounter = 0
    names = map(fs) do f
        name = string(f)
        if name == "(anonymous function)" # Anonymous functions with Julia < 0.5
            λcounter += 1
            name = "λ$(λcounter)"
        end
        name
    end
    names
end


text:
#=
Copyright (c) 2015, Intel Corporation

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:

    * Redistributions of source code must retain the above copyright notice,
      this list of conditions and the following disclaimer.
    * Redistributions in binary form must reproduce the above copyright
      notice, this list of conditions and the following disclaimer in the
      documentation and/or other materials provided with the distribution.
    * Neither the name of Intel Corporation nor the names of its contributors
      may be used to endorse or promote products derived from this software
      without specific prior written permission.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE
FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
=#

export LSTMLayer

function LSTMLayer(name::Symbol, net::Net, input_ensemble::AbstractEnsemble, n_outputs::Int)
    @assert ndims(input_ensemble) == 1

    for ens in [:i, :C_sim, :f, :o]
        @eval $(symbol(ens, :h)) = FullyConnectedEnsemble($net, length($input_ensemble), $n_outputs)
        @eval $(symbol(ens, :x)) = InnerProductLayer($net, $input_ensemble, $n_outputs)
    end

    # i = σ(ih + ix)
    i = σ(net, +(net, ih, ix))
    # C_sim = tanh(C_simh + C_simx)
    C_sim = tanh(net, +(net, C_simh, C_simx))
    # f = σ(fh + fx)
    f = σ(net, +(net, fh, fx))

    # C_t = i * C_sim + f * C_{t-1}
    f_C = MulEnsemble(net, (n_outputs, ))
    add_connections(net, f, f_C, (i) -> (i,))
    C = +(net, *(net, i, C_sim), f_C)
    add_connections(net, C, f_C, (i) -> (i,); recurrent=true)

    # oC = V * C
    oC = InnerProductLayer(net, C, n_outputs)

    # o = sigmoid(oh + ox + oC)
    o = σ(net, +(net, oC, oh, ox))

    # h = o * tanh(C)
    h = *(net, o, tanh(net, C; copy=true))

    # Connect h back to each gate
    for ens in [ih, C_simh, fh, oh]
        add_connections(net, h, ens, (i) -> (1:n_outputs, ); recurrent=true)
    end
end


text:
"""
Copyright Marcus Greiff 2020

Tests for the matrix_transposition() function
"""

# Include modules and recompile if running test standalone
if !isdefined(Main, :runAllTests)
    using Test, AerialVehicleControl, LinearAlgebra
    recompile();
end

# Testing fixture
function fixture(mA,nA,mB,nB)
    A  = rand(Cdouble,mA,nA)
    Ac = deepcopy(A)
    Am = matrix_double_t(A)
    B  = rand(Cdouble,mB,nB)
    Bc = deepcopy(B)
    Bm = matrix_double_t(B)
    return A, Ac, Am, B, Bc, Bm
end

# Compact call for the matrix_double_addition function
function test_call_matrix_double_transposition(Am, Bm)
    status = ccall((:matrix_double_transposition, AerialVehicleControl.CONT_LIB_PATH),
        Cint,
        (Ref{matrix_double_t}, Ref{matrix_double_t}),
        Am,
        Bm)
    return status
end

@verbose(1, "Testing matrix transposition...")
@testset "Matrix transposition" begin
    @testset "Good inputs" begin
        for m = 1:10
            for n = 1:10
                # Fixture
                mA, nA, mB, nB = m, n, n, m
                A, Ac, Am, B, Bc, Bm = fixture(mA,nA,mB,nB)
                # Call
                status = test_call_matrix_double_transposition(Am, Bm)
                # Test
                @test isapprox(status, 1)
                @test isapprox(A, Ac)
                @test !isapprox(B, Bc)
                @test isapprox(B, A')
            end
        end
    end
end


text:
# Copyright (c) 2018: Matthew Wilhelm & Matthew Stuber.
# This code is licensed under MIT license (see LICENSE.md for full details)
#############################################################################
# McCormick.jl
# A McCormick operator library in Julia
# See https://github.com/PSORLab/McCormick.jl
#############################################################################
# src/forward_operators/arithmetic.jl
# Contains definitions of +, -, /, *, promotions, conversion, one, zero.
#############################################################################

# Defines functions required for linear algebra packages
@inline nan(::Type{MC{N,T}}) where {N, T <: RelaxTag} = MC{N,T}(NaN, NaN, Interval{Float64}(NaN),
                                                                fill(NaN, SVector{N,Float64}),
																fill(NaN, SVector{N,Float64}), true)
@inline nan(x::MC{N,T}) where {N, T <: RelaxTag} = MC{N,T}(NaN, NaN, Interval{Float64}(NaN),
                                                           fill(NaN, SVector{N,Float64}),
														   fill(NaN, SVector{N,Float64}), true)

@inline one(::Type{MC{N,T}}) where {N, T <: RelaxTag} = MC{N,T}(1.0, 1.0, one(Interval{Float64}),
                                                                zero(SVector{N,Float64}),
																zero(SVector{N,Float64}), true)
@inline one(x::MC{N,T}) where {N, T <: RelaxTag} = MC{N,T}(1.0, 1.0, one(Interval{Float64}),
                                                           zero(SVector{N,Float64}),
														   zero(SVector{N,Float64}), true)

@inline zero(::Type{MC{N,T}}) where {N, T <: RelaxTag} = MC{N,T}(0.0, 0.0, zero(Interval{Float64}),
                                                                 zero(SVector{N,Float64}),
																 zero(SVector{N,Float64}), true)
@inline zero(x::MC{N,T}) where {N, T <: RelaxTag} = zero(MC{N,T})

@inline real(x::MC) = x
@inline dist(x1::MC, x2::MC) = max(abs(x1.cc - x2.cc), abs(x1.cv - x2.cv))
@inline eps(x::MC) = max(eps(x.cc), eps(x.cv))
@inline mid(x::MC) = mid(x.Intv)

# Unsafe addition
@inline function plus_kernel(x::MC{N,T}, y::MC{N,T}, z::Interval{Float64}) where {N, T <: Union{NS, MV, Diff}}
	MC{N,T}(x.cv + y.cv, x.cc + y.cc, z, x.cv_grad + y.cv_grad, x.cc_grad + y.cc_grad, (x.cnst && y.cnst))
end
@inline +(x::MC, y::MC) = plus_kernel(x, y, x.Intv + y.Intv)
@inline plus_kernel(x::MC, y::Interval{Float64}) = x
@inline +(x::MC) = x

@inline minus_kernel(x::MC{N,T}, z::Interval{Float64}) where {N, T <: RelaxTag} = MC{N,T}(-x.cc, -x.cv, z, -x.cc_grad, -x.cv_grad, x.cnst)
@inline -(x::MC) = minus_kernel(x, -x.Intv)
@inline -(x::MC{N,T}, y::MC{N,T}) where {N, T <: RelaxTag} = minus_kernel(x, y, x.Intv - y.Intv)

# Unsafe subtraction
@inline function minus_kernel(x::MC{N,T}, y::MC{N,T}, z::Interval{Float64}) where {N, T <: Union{NS, MV, Diff}}
	MC{N,T}(x.cv - y.cc, x.cc - y.cv, z, x.cv_grad - y.cc_grad, x.cc_grad - y.cv_grad, (x.cnst && y.cnst))
end

################## CONVERT THROUGH BINARY DEFINITIONS #########################
# Unsafe scalar addition
@inline function plus_kernel(x::MC{N,T}, y::Float64, z::Interval{Float64}) where {N, T <: Union{NS, MV, Diff}}
	MC{N,T}(x.cv + y, x.cc + y, z, x.cv_grad, x.cc_grad, x.cnst)
end
@inline +(x::MC, y::Float64) = plus_kernel(x, y, x.Intv + y)
@inline +(y::Float64, x::MC) = plus_kernel(x, y, x.Intv + y)
@inline +(x::MC{N,T}, y::Interval{Float64}) where {N, T<:RelaxTag} = x + MC{N,T}(y)
@inline +(y::Interval{Float64}, x::MC{N,T}) where {N, T<:RelaxTag} = x + MC{N,T}(y)

@inline plus_kernel(x::MC, y::C, z::Interval{Float64}) where {C <: NumberNotRelax} = plus_kernel(x, convert(Float64, y), z)
@inline plus_kernel(x::C, y::MC, z::Interval{Float64}) where {C <: NumberNotRelax} = plus_kernel(y, convert(Float64, x), z)
@inline +(x::MC, y::C) where {C <: NumberNotRelax} = x + convert(Float64, y)
@inline +(y::C, x::MC) where {C <: NumberNotRelax} = x + convert(Float64, y)

# Unsafe scalar subtraction
@inline function minus_kernel(x::MC{N,T}, c::Float64, z::Interval{Float64}) where {N, T <: Union{NS, MV, Diff}}
	MC{N,T}(x.cv - c, x.cc - c, z, x.cv_grad, x.cc_grad, x.cnst)
end
@inline function minus_kernel(c::Float64, x::MC{N,T}, z::Interval{Float64}) where {N, T <: Union{NS, MV, Diff}}
	MC{N,T}(c - x.cc, c - x.cv, z, -x.cc_grad, -x.cv_grad, x.cnst)
end
@inline -(x::MC, c::Float64) = minus_kernel(x, c, x.Intv - c)
@inline -(c::Float64, x::MC) = minus_kernel(c, x, c - x.Intv)
@inline -(x::MC{N,T}, y::Interval{Float64}) where {N, T<:RelaxTag} = x - MC{N,T}(y)
@inline -(y::Interval{Float64}, x::MC{N,T}) where {N, T<:RelaxTag} = MC{N,T}(y) - x

@inline minus_kernel(x::MC, y::C, z::Interval{Float64}) where {C <: NumberNotRelax} = minus_kernel(x, convert(Float64, y), z)
@inline minus_kernel(y::C, x::MC, z::Interval{Float64}) where {C <: NumberNotRelax} = minus_kernel(convert(Float64, y), x, z)
@inline -(x::MC, c::C) where {C <: NumberNotRelax} = x - convert(Float64,c)
@inline -(c::C, x::MC) where {C <: NumberNotRelax} = convert(Float64,c) - x

# Unsafe Scalar Multiplication
@inline function mult_kernel(x::MC{N,T}, c::Float64, z::Interval{Float64}) where {N, T <: Union{NS, MV, Diff}}
	if c >= 0.0
		zMC = MC{N,T}(c*x.cv, c*x.cc, z, c*x.cv_grad, c*x.cc_grad, x.cnst)
	else
		zMC = MC{N,T}(c*x.cc, c*x.cv, z, c*x.cc_grad, c*x.cv_grad, x.cnst)
	end
	return zMC
end
@inline *(x::MC, c::Float64) = mult_kernel(x, c, c*x.Intv)
@inline *(c::Float64, x::MC) = mult_kernel(x, c, c*x.Intv)
@inline *(x::MC{N,T}, y::Interval{Float64}) where {N, T<:RelaxTag} = x*MC{N,T}(y)
@inline *(y::Interval{Float64}, x::MC{N,T}) where {N, T<:RelaxTag} = MC{N,T}(y)*x

@inline mult_kernel(x::MC, c::C, z::Interval{Float64}) where {C <: NumberNotRelax} = mult_kernel(x, convert(Float64, c), z)
@inline mult_kernel(c::C, x::MC, z::Interval{Float64}) where {C <: NumberNotRelax} =  mult_kernel(x, convert(Float64, c), z)
@inline *(c::C, x::MC) where {C <: NumberNotRelax}  = x*Float64(c)
@inline *(x::MC, c::C) where {C <: NumberNotRelax}  = x*Float64(c)

# Unsafe scalar division
@inline div_kernel(x::MC, y::Float64, z::Interval{Float64}) = mult_kernel(x, inv(y), z)
@inline div_kernel(x::Float64, y::MC, z::Interval{Float64}) = mult_kernel(inv(y), x, z)
@inline div_kernel(x::MC, y::C, z::Interval{Float64}) where {C <: NumberNotRelax}  = mult_kernel(x, inv(y), z)
@inline div_kernel(x::C, y::MC, z::Interval{Float64}) where {C <: NumberNotRelax}  = mult_kernel(inv(y), x, z)
@inline /(x::MC, y::Float64) = x*inv(y)
@inline /(x::Float64, y::MC) = x*inv(y)
@inline /(x::MC, y::C) where {C <: NumberNotRelax} = x*inv(convert(Float64,y))
@inline /(x::C, y::MC) where {C <: NumberNotRelax} = convert(Float64,x)*inv(y)
@inline /(x::MC{N,T}, y::Interval{Float64}) where {N, T<:RelaxTag} = x/MC{N,T}(y)
@inline /(y::Interval{Float64}, x::MC{N,T}) where {N, T<:RelaxTag} = MC{N,T}(y)/x

# Maximization
@inline max_kernel(c::Float64, x::MC, z::Interval{Float64}) = max_kernel(x, c, z)
@inline max_kernel(x::MC, c::C, z::Interval{Float64}) where {C <: NumberNotRelax} = max_kernel(x, convert(Float64, c), z)
@inline max_kernel(c::C, x::MC, z::Interval{Float64}) where {C <: NumberNotRelax} = max_kernel(x, convert(Float64, c), z)

@inline max(c::Float64, x::MC) = max_kernel(x, c, max(x.Intv, c))
@inline max(x::MC, c::C) where {C <: NumberNotRelax} = max_kernel(x, convert(Float64, c), max(x.Intv, c))
@inline max(c::C, x::MC) where {C <: NumberNotRelax} = max_kernel(x, convert(Float64, c), max(x.Intv, c))
@inline max(x::MC{N,T}, y::Interval{Float64}) where {N, T<:RelaxTag} = max(x, MC{N,T}(y))
@inline max(y::Interval{Float64}, x::MC{N,T}) where {N, T<:RelaxTag} = max(MC{N,T}(y), x)

# Minimization
@inline min_kernel(x::MC, c::C, z::Interval{Float64}) where {C <: NumberNotRelax} = min_kernel(x, convert(Float64, c), z)
@inline min_kernel(c::C, x::MC, z::Interval{Float64}) where {C <: NumberNotRelax} = min_kernel(x, convert(Float64, c), z)

@inline min(c::Float64, x::MC) = min_kernel(x, c, min(x.Intv, c))
@inline min(x::MC, c::C) where {C <: NumberNotRelax} = min_kernel(x, convert(Float64, c), min(x.Intv, c))
@inline min(c::C, x::MC) where {C <: NumberNotRelax} = min_kernel(x, convert(Float64, c), min(x.Intv, c))
@inline min(x::MC{N,T}, y::Interval{Float64}) where {N, T<:RelaxTag} = min(x, MC{N,T}(y))
@inline min(y::Interval{Float64}, x::MC{N,T}) where {N, T<:RelaxTag} = min(MC{N,T}(y), x)

# Add fma function
@inline fma(x::MC, y::Float64, z::Float64) = x*y + z
@inline fma(x::MC, y::MC, z::Float64) = x*y + z
@inline fma(x::MC, y::Float64, z::MC) = x*y + z
@inline fma(x::MC, y::MC, z::MC) = x*y + z
@inline fma(x::Float64, y::MC, z::Float64) = x*y + z
@inline fma(x::Float64, y::MC, z::MC) = x*y + z
@inline fma(x::Float64, y::Float64, z::MC) = x*y + z

@inline fma(x::MC, y::Float64, z::Float64, q::Interval{Float64}) = x*y + z
@inline fma(x::MC, y::MC, z::Float64, q::Interval{Float64}) = x*y + z
@inline fma(x::MC, y::Float64, z::MC, q::Interval{Float64}) = x*y + z
@inline fma(x::MC, y::MC, z::MC, q::Interval{Float64}) = x*y + z
@inline fma(x::Float64, y::MC, z::Float64, q::Interval{Float64}) = x*y + z
@inline fma(x::Float64, y::MC, z::MC, q::Interval{Float64}) = x*y + z
@inline fma(x::Float64, y::Float64, z::MC, q::Interval{Float64}) = x*y + z

# Promote and Convert
promote_rule(::Type{MC{N,T}}, ::Type{S}) where {S<:NumberNotRelax, N, T <: RelaxTag} = MC{N,T}
promote_rule(::Type{MC{N,T}}, ::Type{S}) where {S<:Real, N, T <: RelaxTag} = MC{N,T}

convert(::Type{MC{N,T}}, x::S) where {S<:NumberNotRelax, N, T <: RelaxTag} = MC{N,T}(Interval{Float64}(x))
convert(::Type{MC{N,T}}, x::S) where {S<:AbstractInterval, N, T <: RelaxTag} = MC{N,T}(Interval{Float64}(x.lo, x.hi))
Interval(x::MC{N,T}) where {N,T<:RelaxTag} = x.Intv
Interval{Float64}(x::MC{N,T}) where {N,T<:RelaxTag} = x.Intv

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
using ADCME
using PyCall
using LinearAlgebra
using PyPlot
using Random
using AdFem
Random.seed!(233)

function univariate_fem_stiffness(hmat,m,n,h)
    univariate_fem_stiffness_ = load_op_and_grad("./build/libUnivariateFemStiffness","univariate_fem_stiffness", multiple=true)
    hmat,m,n,h = convert_to_tensor([hmat,m,n,h], [Float64,Int32,Int32,Float64])
    univariate_fem_stiffness_(hmat,m,n,h)
end

m = 10
n = 5
h = 1.0
hmat = zeros(4*m*n, 2, 2)
for i = 1:4*m*n
    hmat[i,:,:] = diagm(0=>ones(2))
end

# TODO: specify your input parameters
u = univariate_fem_stiffness(hmat,m,n,h)
sess = Session(); init(sess)
@show run(sess, u)

# uncomment it for testing gradients
H = rand(2,2)
H = H+ H'
K = compute_fem_stiffness_matrix1(H, m, n, h)
u = univariate_fem_stiffness(H,m,n,h)
A = SparseTensor(u..., (m+1)*(n+1), (m+1)*(n+1))
@show maximum(abs.(run(sess, A)-K))
# error()

# TODO: change your test parameter to `m`
#       in the case of `multiple=true`, you also need to specify which component you are testings
# gradient check -- v
function scalar_function(hmat)
    return sum(univariate_fem_stiffness(hmat,m,n,h)[3]^2)
end

# TODO: change `m_` and `v_` to appropriate values
# m_ = constant(rand(4*m*n, 2, 2))
# v_ = rand(4*m*n, 2, 2)

m_ = constant(rand( 2, 2))
v_ = rand( 2, 2)
y_ = scalar_function(m_)
dy_ = gradients(y_, m_)
ms_ = Array{Any}(undef, 5)
ys_ = Array{Any}(undef, 5)
s_ = Array{Any}(undef, 5)
w_ = Array{Any}(undef, 5)
gs_ =  @. 1 / 10^(1:5)

for i = 1:5
    g_ = gs_[i]
    ms_[i] = m_ + g_*v_
    ys_[i] = scalar_function(ms_[i])
    s_[i] = ys_[i] - y_
    w_[i] = s_[i] - g_*sum(v_.*dy_)
end

sess = Session(); init(sess)
sval_ = run(sess, s_)
wval_ = run(sess, w_)
close("all")
loglog(gs_, abs.(sval_), "*-", label="finite difference")
loglog(gs_, abs.(wval_), "+-", label="automatic differentiation")
loglog(gs_, gs_.^2 * 0.5*abs(wval_[1])/gs_[1]^2, "--",label="\$\\mathcal{O}(\\gamma^2)\$")
loglog(gs_, gs_ * 0.5*abs(sval_[1])/gs_[1], "--",label="\$\\mathcal{O}(\\gamma)\$")

plt.gca().invert_xaxis()
legend()
xlabel("\$\\gamma\$")
ylabel("Error")

savefig("test.png")


text:
# test the DynPoint2 functionality

using RoME
using Test
# using RoMEPlotting,
# , IncrementalInference, Distributions
# using KernelDensityEstimate #, KernelDensityEstimatePlotting

@testset "test DynPoint2..." begin

# N = 75
global fg = initfg()

# add two point locations
global v0 = addVariable!(fg, :x0, DynPoint2; nanosecondtime=0)

global v1 = addVariable!(fg, :x1, DynPoint2; nanosecondtime=10^9)

# Prior factor as boundary condition
global pp0 = DynPoint2VelocityPrior(MvNormal([zeros(2);10*ones(2)], 0.1*Matrix{Float64}(LinearAlgebra.I, 4,4)))
global f0 = addFactor!(fg, [:x0;], pp0)
# conditional likelihood between Dynamic Point2
global dp2dp2 = DynPoint2DynPoint2(MvNormal([10*ones(2);zeros(2)], 0.1*Matrix{Float64}(LinearAlgebra.I, 4,4)))
global f1 = addFactor!(fg, [:x0;:x1], dp2dp2)

# Graphs.plot(fg.g)
ensureAllInitialized!(fg)

# global tree = wipeBuildNewTree!(fg)
# inferOverTree!(fg, tree)
tree, smt, hist = solveTree!(fg)

# X1 = getVal(fg, :x1)
@show global x0 = getKDEMax(getBelief(fg, :x0))

@test norm(x0[1:2] - [0.0;0.0]) < 1.0
@test norm(x0[3:4] - [10.0;10.0]) < 1.0

@show global x1 = getKDEMax(getBelief(fg, :x1))

@test norm(x1[1:2] - [20.0;20.0]) < 1.0
@test norm(x1[3:4] - [10.0;10.0]) < 1.0

end

@testset "test DynPoint2VelocityPrior..." begin

# N = 75
global fg = initfg()

# add two point locations
global v0 = addVariable!(fg, :x0, DynPoint2; nanosecondtime=0)
global v1 = addVariable!(fg, :x1, DynPoint2; nanosecondtime=10^9)
global v2 = addVariable!(fg, :x2, DynPoint2; nanosecondtime=2*10^9)

# Prior factor as boundary condition
global pp0 = DynPoint2VelocityPrior(MvNormal([zeros(2);10*ones(2)], 0.1*Matrix{Float64}(LinearAlgebra.I, 4,4)))
global f0 = addFactor!(fg, [:x0;], pp0)

# conditional likelihood between Dynamic Point2
global dp2dp2 = VelPoint2VelPoint2(MvNormal([10*ones(2);zeros(2)], 0.1*Matrix{Float64}(LinearAlgebra.I, 4,4)))
global f1 = addFactor!(fg, [:x0;:x1], dp2dp2)

# conditional likelihood between Dynamic Point2
global dp2dp2 = VelPoint2VelPoint2(MvNormal([10*ones(2);zeros(2)], 0.1*Matrix{Float64}(LinearAlgebra.I, 4,4)))
global f2 = addFactor!(fg, [:x1;:x2], dp2dp2)

# Graphs.plot(fg.g)
ensureAllInitialized!(fg)

@test isInitialized(fg, :x0)
@test isInitialized(fg, :x1)
@test isInitialized(fg, :x2)

# global tree = wipeBuildNewTree!(fg)
# inferOverTree!(fg, tree)
tree, smt, hist = solveTree!(fg)

# X1 = getVal(fg, :x1)
@show global x0 = getKDEMax(getBelief(fg, :x0))
@test norm(x0[1:2] - [0.0;0.0]) < 1.0
@test norm(x0[3:4] - [10.0;10.0]) < 1.0

@show global x1 = getKDEMax(getBelief(fg, :x1))
@test norm(x1[1:2] - [10.0;10.0]) < 1.0
@test norm(x1[3:4] - [10.0;10.0]) < 1.0

@show global x2 = getKDEMax(getBelief(fg, :x2))
@test norm(x2[1:2] - [20.0;20.0]) < 1.0
@test norm(x2[3:4] - [10.0;10.0]) < 1.0

end

@testset "test VelPoint2VelPoint2" begin

global N=100
global pμ = [0.0,0,10.0,0]
global pσ = Matrix(Diagonal([0.1;0.1;0.1;0.1].^2))

global fg = initfg();

addVariable!(fg, :x1, DynPoint2;  nanosecondtime=0)
global pp = DynPoint2VelocityPrior(MvNormal(pμ,pσ))
addFactor!(fg, [:x1;], pp, graphinit=false)

ensureAllInitialized!(fg)

addVariable!(fg, :x2, DynPoint2;  nanosecondtime=1_000_000_000)
global dpμ = [10.0;0;0;0];
global dpσ = Matrix(Diagonal([1.0;1;0.5;0.01].^2))
global vp = VelPoint2VelPoint2(MvNormal(dpμ,dpσ))
addFactor!(fg, [:x1,:x2], vp, graphinit=false)

# writeGraphPdf(fg)

# note, nothign is initialized yet...
@test isInitialized(fg, :x1)
# @test !isInitialized(fg, :x2)

# lets init the first
# IncrementalInference.doautoinit!(fg, :x1)

global X1 = getVal(fg, :x1)
@test 0.7*N < sum(-0.5 .< X1[1,:] .< 0.5)
@test 0.7*N < sum(-0.5 .< X1[2,:] .< 0.5)
@test 0.7*N < sum(9.0 .< X1[3,:] .< 11.0)
@test 0.7*N < sum(9.0 .< X1[3,:] .< 11.0)

# now test "projection" / probabilistic-convolution through VelPointVelPoint{D, N}
global pX2 = approxConv(fg, :x1x2f1, :x2)

@test 0.7*N < sum(9 .< pX2[1,:] .< 11)
@test 0.7*N < sum(-2 .< pX2[2,:] .< 2)
@test 0.7*N < sum(9.0 .< pX2[3,:] .< 11.0)
@test 0.7*N < sum(9.0 .< pX2[3,:] .< 11.0)

# using RoMEPlotting, Gadfly
# using KernelDensityEstimate
# using KernelDensityEstimatePlotting
# X1 = getVal(fg,:x1)
# plot(x=pX2[1,:], y=pX2[2,:], Geom.hexbin)
# pl = plotKDE(kde!(pX2[1:2,:]), dims=[1;2], levels=3, c=["blue"])
# pl = plotKDE(kde!(X1[1:2,:]), dims=[1;2], levels=3, c=["blue"])

end

#


text:
using Documenter
using Pkg
Pkg.activate(".")
Pkg.instantiate()
using ATA

makedocs(
    sitename = "ATA",
    format = Documenter.HTML(),
    modules = [ATA],
    doctest = true,
    pages = ["index.md", "model_build.md", "opt.md", "compact.md", "print_and_plot.md", "structs.md", "utils.md", "examples.md"],
)

deploydocs(repo = "github.com/giadasp/ATA.jl.git", devurl = "docs", devbranch = "release")


text:
struct VEMFunctionSpace{dim,T<:Real,FE<:AbstractVirtualElement,mtype <: AbstractPolytopalMesh} <: AbstractVEMFunctionSpace{dim,T,FE}
  element::FE
  mesh::mtype
  components::Int
end

getelement(fs::VEMFunctionSpace) = fs.element
getmesh(fs::VEMFunctionSpace) = fs.mesh

function Base.show(io::IO, fs::VEMFunctionSpace)
  if fs.components == 1
      println(io, "Discrete scalar Virtual function space")
  else
      println(io, "Discrete vectorial Virtual function space (components = $(fs.components))")
  end
  println(io, "Virtual Element: \n", getelement(fs))
  println(io, "Domain mesh: \n", getmesh(fs))
end

function VEMFunctionSpace(mesh::AbstractPolytopalMesh{dim,T}, elem::AbstractVirtualElement; components::Int=1) where {dim,T}
  VEMFunctionSpace{dim,T,typeof(elem),typeof(mesh)}(elem,mesh,components)
end

function getnlocaldofs(fs::VEMFunctionSpace, cell::Cell)
  getnlocaldofs(fs.element, cell)*fs.components
end

getncomponents(fs::VEMFunctionSpace) = fs.components

max_stars_repo_name: Paulms/jFEMToools

max_issues_repo_name: Paulms/jFEMToools

max_forks_repo_name: Paulms/jFEMToools


text:
# # Electromechanic break
#md # [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/models/EMBrake.ipynb)
#
#md # !!! note "Overview"
#md #     System type: Affine system\
#md #     State dimension: 48\
#md #     Application domain: Mechanical Engineering
#
# ## Model description
#
# The system is described by the linear differential equations:
#
# ```math
#   \begin{array}{lcl}
#   \dot{x}(t) &=& Ax(t) + Bu(t),\qquad u(t) \in \mathcal{U} \\
#   y(t) &=& C x(t)
#   \end{array}
# ```
#
# There are two versions of this benchmark:
#
# - *(time-varying inputs):* The inputs can change arbitrarily over time: $\forall t: u(t)\in \mathcal{U}$.
#
# - *(constant inputs):* The inputs are uncertain only in their initial value, and
#    constant over time: ``u(0)\in \mathcal{U}``, ``\dot u (t)= 0``.

using ReachabilityAnalysis, SparseArrays, JLD2 #!jl

# adapt default LazySets tolerances
set_rtol(Float64, 1e-12)
set_ztol(Float64, 1e-13)

## Fixed parameters

# The model without parameter variation is defined below.

function embrake_no_pv(; Tsample=1.E-4, ζ=1e-6, x0=0.05)
    ## model's constants
    L = 1.e-3
    KP = 10000.
    KI = 1000.
    R = 0.5
    K = 0.02
    drot = 0.1
    i = 113.1167

    ## state variables: [I, x, xe, xc]
    A = Matrix([-(R+K^2/drot)/L 0      KP/L  KI/L;
                          K/i/drot      0      0      0;
                          0             0      0      0;
                          0             0      0      0])

    EMbrake = @system(x' = Ax)

    ## initial conditions
    I₀  = Singleton([0.0])
    x₀  = Singleton([0.0])
    xe₀ = Singleton([0.0])
    xc₀ = Singleton([0.0])
    X₀ = I₀ × x₀ × xe₀ × xc₀

    ## reset map
    Ar = sparse([1, 2, 3, 4, 4], [1, 2, 2, 2, 4], [1., 1., -1., -Tsample, 1.], 4, 4)
    br = sparsevec([3, 4], [x0, Tsample*x0], 4)
    reset_map(X) = Ar * X + br

    ## hybrid system with clocked linear dynamics
    ha = HACLD1(EMbrake, reset_map, Tsample, ζ)
    return IVP(ha, X₀)
end

## Parameter variation

# The model with parameter variation described below consists of changing only one
# coefficient, which corresponds to the Flow* settings in [^SO15].
function embrake_pv_1(; Tsample=1.E-4, ζ=1e-6, Δ=3.0, x0=0.05)
    ## model's constants
    L = 1.e-3
    KP = 10000.
    KI = 1000.
    R = 0.5
    K = 0.02
    drot = 0.1
    i = 113.1167
    p = 504. + (-Δ .. Δ)

    ## state variables: [I, x, xe, xc]
    A = IntervalMatrix([-p            0      KP/L   KI/L;
                        K/i/drot      0      0      0;
                        0             0      0      0;
                        0             0      0      0])

    EMbrake = @system(x' = Ax)

    ## initial conditions
    I₀  = Singleton([0.0])
    x₀  = Singleton([0.0])
    xe₀ = Singleton([0.0])
    xc₀ = Singleton([0.0])
    X₀ = I₀ × x₀ × xe₀ × xc₀

    ## reset map
    Ar = sparse([1, 2, 3, 4, 4], [1, 2, 2, 2, 4], [1., 1., -1., -Tsample, 1.], 4, 4)
    br = sparsevec([3, 4], [x0, Tsample*x0], 4)
    reset_map(X) = Ar * X + br

    ## hybrid system with clocked linear dynamics
    ha = HACLD1(EMbrake, reset_map, Tsample, ζ)
    return IVP(ha, X₀)
end

## Extended parameter variation

# In the following model, considered in [...], we vary all constants by χ% with
# respect to their nominal values. The variation percentage defaults to 5%.

function embrake_pv_2(; Tsample=1.E-4, ζ=1e-6, x0=0.05, χ=5.0)
    ## model's constants
    Δ = -χ/100 .. χ/100
    L = 1.e-3 * (1 + Δ)
    KP = 10000. * (1 + Δ)
    KI = 1000. * (1 + Δ)
    R = 0.5 * (1 + Δ)
    K = 0.02 * (1 + Δ)
    drot = 0.1 * (1 + Δ)
    i = 113.1167 * (1 + Δ)

    ## state variables: [I, x, xe, xc]
    A = IntervalMatrix([-(R+K^2/drot)/L 0      KP/L  KI/L;
                          K/i/drot      0      0      0;
                          0             0      0      0;
                          0             0      0      0])

    EMbrake = @system(x' = Ax)

    ## initial conditions
    I₀  = Singleton([0.0])
    x₀  = Singleton([0.0])
    xe₀ = Singleton([0.0])
    xc₀ = Singleton([0.0])
    X₀ = I₀ × x₀ × xe₀ × xc₀

    ## reset map
    Ar = sparse([1, 2, 3, 4, 4], [1, 2, 2, 2, 4], [1., 1., -1., -Tsample, 1.], 4, 4)
    br = sparsevec([3, 4], [x0, Tsample*x0], 4)
    reset_map(X) = Ar * X + br

    ## hybrid system with clocked linear dynamics
    ha = HACLD1(EMbrake, reset_map, Tsample, ζ)
    return IVP(ha, X₀)
end

# ## References

# [^SO15]: Strathmann, Thomas, and Jens Oehlerking. *Verifying Properties of an Electro-Mechanical Braking System.* ARCH CPSWeek. 2015.

# []:


text:
#=
cuboid_im.jl
2019-07-01 Helena
=#

using Plots
#using MIRT

"""
`phantom = cuboid_im(ig, params;
    oversample=1, hu_scale=1, how=:auto, return_params=false)`

generate cuboid phantom image from parameters:
    `[x_center y_center z_center x_diameter y_diameter z_diameter
        xy_angle_degrees z_angle_degrees amplitude]`

in
* `ig`            `image_geom()` object
* `params`        `[9 N]`  cuboid parameters. if empty, use defualt
                           note: "diameter" not "radius"

options
* `oversample::Integer`   oversampling factor
* `how::Symbol`          `:sample` use samples
                        	`:lowmem1` one slice per time
                            `:exact` perfect partial volume if angle* = 0
                            default: `:exact` if non rotated, else `:sample`
* `return_params::Bool`    if true, return both phantom and params

out
* `phantom`         `[nx ny nz]` image
* `params`          `[9 N]` cuboid parameters (only return if return_params=true)
"""
function cuboid_im(ig::MIRT_image_geom,
    params::AbstractArray{<:Real,2};
    oversample::Integer=1,
    how::Symbol=:auto,
    return_params::Bool=false)

    size(params,2) != 9 && throw("bad cuboid parameter vector size")

    if oversample > 1
        ig = ig.over(oversample)
    end
    args = (ig.nx, ig.ny, ig.nz, params,
    ig.dx, ig.dy, ig.dz, ig.offset_x, ig.offset_y, ig.offset_z, oversample)

	if any(params[:,4:6] .< 0)
		throw("cuboid `diameters` must be nonnegative")
	end

	phantom = zeros(Float32, ig.nx, ig.ny, ig.nz)

	rotated = params[:, 7:8] .!= 0

	if how == :auto
		if any(rotated)
			how = :sample
		else
			how = :exact
		end
	end

    if how == :exact
        phantom += cuboid_im_exact(args...)
        if oversample != 1
            throw("ignoring oversample $oversample")
        end
    elseif how == :lowmem1
        phantom += cuboid_im_lowmem1(args...)
    elseif how == :sample
        phantom += cuboid_im_sample(args...)
    else
        throw("bad how $how")
    end

    if return_params
        return (phantom, params)
    end
    return phantom
end

"""
`cuboid_im_sample()``

:sample
"""
function cuboid_im_sample(nx, ny, nz, params, dx, dy, dz, offset_x, offset_y, offset_z, over)
    phantom = zeros(Float32, nx*over, ny*over, nz*over)

    wx = (nx*over-1)/2 + offset_x*over
    wy = (ny*over-1)/2 + offset_y*over
    wz = (nz*over-1)/2 + offset_z*over
    x = ((0:nx*over-1) .- wx) * dx / over
    y = ((0:ny*over-1) .- wy) * dy / over
    z = ((0:nz*over-1) .- wz) * dz / over
    (xx, yy, zz) = ndgrid(x, y ,z)

    # ticker reset
    ne = size(params)[1]
    for ie in 1:ne
        # ticker(mfile, ie, ne)

        par = params[ie, :]
        cx = par[1]
        cy = par[2]
        cz = par[3]
        rx = par[4]
        ry = par[5]
        rz = par[6]

        theta = par[7] * (pi/180)
        phi = par[8] * (pi/180)
        if phi != 0
            throw("phi rotation not done")
        end

        (x, y) = rotate2d(xx .- cx, yy .- cy, theta)
        z = zz .- cz

		# rx, ry, rz are 'diameters' not 'radius'
		tmp = (abs.(x / rx) .<= 1/2) .& (abs.(y / ry) .<= 1/2) .& (abs.(z / rz) .<= 1/2)
		value = Float32(par[9])
		phantom += value * tmp
	end
	phantom = downsample3(phantom, over)

    return phantom
end

"""
`cuboid_im_exact()``

:exact

non-rotated cuboid -- rotation not done
"""
function cuboid_im_exact(nx, ny, nz, params, dx, dy, dz, offset_x, offset_y, offset_z, over)
	phantom = zeros(Float32, nx, ny, nz)

	wx = (nx-1)/2 + offset_x
	wy = (ny-1)/2 + offset_y
	wz = (nz-1)/2 + offset_z
	x = ((0:nx-1) .- wx) * dx
	y = ((0:ny-1) .- wy) * dy
	z = ((0:nz-1) .- wz) * dz
	(xx, yy, zz) = ndgrid(x, y, z)

	# length of intersection of [a, b] with [c, d]
	fun = (a, b, c, d) ->
		max(min(d,b) - max(a,c), 0)
	fun2 = (x, h) ->
		fun(x - 0.5, x + 0.5, -h, h)

	adx = abs(dx)
	ady = abs(dy)
	adz = abs(dz)

	#ticker reset
	ne = size(params)[1]
	for ie in 1:ne
		#ticker(mfilename, ie, ne)

		par = params[ie, :]
		cx = par[1]
		cy = par[2]
		cz = par[3]
		rx = par[4]
		ry = par[5]
		rz = par[6]

		if par[7] != 0 || par[8] != 0
			throw("rotation not done")
		end

		# rx, ry, rz are "diameters" not "radius"
		x = (xx .- cx) / adx
		hx2 = rx / adx / 2 # half width in pixels
		y = (yy .- cy) / ady
		hy2 = ry / ady / 2
		z = (zz .- cz) / adz
		hz2 = rz / adz / 2

		fx = fun2.(x, hx2)
		fy = fun2.(y, hy2)
		fz = fun2.(z, hz2)
		tmp = (fx .* fy .* fz)
		value = Float32(par[9])
		phantom += value * tmp
	end

    return phantom
end

"""
`cuboid_im_lowmem1`

this version does one 'one slice at a time' to reduce memory
"""
function cuboid_im_lowmem1(nx, ny, nz, params, dx, dy, dz, offset_x, offset_y, offset_z, over)
	phantom = zeros(Float32, nx, ny, nz)
	for iz in 1:nz
		offset_z_new = (nz - 1)/2 + offset_z - (iz-1)
		phantom[:,:,iz] = cuboid_im_sample(nx, ny, 1, params,
	    	dx, dy, dz, offset_x, offset_y, offset_z_new, over)
	end
    return phantom
end

"""
`default_parameters()`
"""
function default_parameters(xfov, yfov, zfov)
	params = [
	0	0		0	60	60	2	0	0	1
	20	20		3	10	20	2	0	0	1
	1.4	-0.5	1	2	2.7	3.2	0	0	1
	]
	return params
end

"""
`rotated_parameters()`
"""
function rotated_parameters(xfov, yfov, zfov)
	params = [
	0	0	0	60	60	2	45	0	1
	0	0	0	50	50	3	45	0	1
	]
	return params
end

#=
"""
`phantom = cuboid_im(nx, dx, params; args...)`

specifying voxel size `dx` and cuboid `params`
"""
function cuboid_im(nx::Integer, dx::Real, params; args...)
	ig = image_geom(nx=nx, dx=1)
	return cuboid_im(ig, params; args...)
end

"""
`phantom = cuboid_im(nx::Integer, params; args...)`

voxel size `dx=1` and cuboid `params`
"""
function cuboid_im(nx::Integer, params; args...)
	return cuboid_im(nx, 1., params; args...)
end

"""
`phantom = cuboid_im(nx::Integer; ny::Integer=nx,dx::Real=1, args...)`

image of size 'nx' by 'ny' (default 'nx') with specified 'dx' (default 1),
defaults to `:default_cuboid`
"""
function cuboid_im(nx::Integer; ny::Integer=nx, dx::Real=1, args...)
	ig = image_geom(nx=nx, ny=ny, dx=dx)
	return cuboid_im(ig, :default; args...)
end

"""
`phantom = cuboid_im(nx::Integer, ny::integer; args...)`

`:default` of size `nx` by `ny`
"""
function cuboid_im(nx::Integer, ny::Integer; args...)
	return cuboid_im(nx, ny=ny, dx=1.; args...)
end
=#

"""
`phantom = cuboid_im(ig, code, args...)`

`code = :default | :rotate` # add more options
"""
function cuboid_im(ig::MIRT_image_geom, params::Symbol; args...)
	fov = ig.fovs
	if params == :default
		params = default_parameters(fov, fov, fov)
	elseif params == :rotate
		params = rotated_parameters(fov, fov, fov)
	else
		throw("bad phantom symbol $params")
	end
	return cuboid_im(ig, params; args...)
end

"""
`phantom = cuboid_im(ig; args...)`

`:default` default parameters
"""
function cuboid_im(ig::MIRT_image_geom; args...)
	return cuboid_im(ig, :default; args...)
end

"""
`cuboid_im()`

shows docstring(s)
"""
function cuboid_im()
	@doc cuboid_im
end

"""
`cuboid_im_show()`
"""
function cuboid_im_show()
	ig = image_geom(nx=2^3, ny=2^3, nz=2^3, dz=1, fov=240)

	x1 = cuboid_im(ig, :default, how=:exact)
	p1 = jim(x1, title="exact")

	x2 = cuboid_im(ig, :default, how=:sample)
	p2 = jim(x2, title="sample")

	x3 = cuboid_im(ig, :default, how=:lowmem1)
	p3 = jim(x3, title="lowmem")

	x4 = cuboid_im(ig, :rotate)
	p4 = jim(x4, title="rotate") # runs sample

	plot(p1, p2, p3, p4)
end

"""
`cuboid_im_test()`
"""
function cuboid_im_test()
	ig = image_geom(nx=2^3, ny=2^3, nz=2^3, dz=1, fov=240)
	fov = 100
	cuboid_im_show()

	diam = abs.([2*ig.dx 2.7*ig.dy 3.2*ig.dz])
	params = [1.4 -0.5 1 diam 0 0 1]
	phantom_exact = cuboid_im(ig, params, how=:exact)
	dxyz = abs.(ig.dx * ig.dy * ig.dz)
	vol_true = prod(diam)
	vol_phantom = sum(phantom_exact[:]) * dxyz
	isapprox(vol_true, vol_phantom)

	phantom_sample = cuboid_im(ig, :default, how=:sample)
	vol_phantom2 = sum(phantom_sample[:]) * dxyz
	isapprox(vol_phantom2, vol_true)

	# test
	#x4 = cuboid_im(ig, :rotate, how=:exact) exact does not support rotation
	x3 = cuboid_im(ig, :default, how=:lowmem1)
	x5 = cuboid_im(ig)
	true
end

"""
`cuboid_im(:test)`

`cuboid_im(:show)`

run tests
"""
function cuboid_im(test::Symbol)
	if test == :show
		return cuboid_im_show()
	end
	test != :test && throw(ArgumentError("test $test"))
	cuboid_im()
	cuboid_im(:show)
	cuboid_im_test()
	true
end

#cuboid_im(:show)
#cuboid_im(:test)


text:
using FractrixZeta

game = FractrixZeta.main()

max_stars_repo_name: zoemcc/fractrixzeta

max_issues_repo_name: zoemcc/fractrixzeta

max_forks_repo_name: zoemcc/fractrixzeta


text:
# Common backend methods.
#

# Generate single call, replace [xr]Type with typenames and assign letter.
macro blis_ccall_typed(funcname, ret, ctypename, rtypename, named_param_pairs...)
    # C and Julia function names.
    cfunc  = string("bli_", cblas_typechar[ eval(ctypename) ], funcname)
    jlfunc = Symbol("bli_", cblas_typechar[ eval(ctypename) ], funcname, "!")
    ccfunc = Expr(:string, cfunc)

    # Separate variable names and their types.
    names  = named_param_pairs[1:2:end]
    params = named_param_pairs[2:2:end]

    # For typed parameters and return value, replace [xr]Type with real typenames.
    typed_params = ( iType -> eval(:( ( (xType, rType) -> $iType )($ctypename, $rtypename) )) ).(params)
    typed_return =            eval(:( ( (xType, rType) -> $ret   )($ctypename, $rtypename) ))

    # Convert type specifiers to static symbol.
    typed_param_syms = Expr(:tuple, typed_params...)
    typed_return_sym = Symbol(typed_return)

    # Docstring.
    funcdoc = Expr(:string, """
        Autogenerated BLIS Typed-API function.
        For usage refer the
        [BLIS docs for `$cfunc`.](https://github.com/flame/blis/blob/0.8.0/docs/BLISTypedAPI.md#$funcname)
        """)

    # Define Julia function as ccall.
    return quote
        $(esc(jlfunc))( $(names...) ) = begin
            ccall(dlsym(libblis, $(esc(ccfunc))), $typed_return_sym, $typed_param_syms, $(names...))
        end

        @doc $funcdoc $jlfunc
    end
end

# Ganerate a famity of calls.
macro blis_ccall_group(funcname, ret, named_param_pairs...)
    return quote
        @blis_ccall_typed($funcname, $ret, Cfloat    , Cfloat , $(named_param_pairs...))
        @blis_ccall_typed($funcname, $ret, Cdouble   , Cdouble, $(named_param_pairs...))
        @blis_ccall_typed($funcname, $ret, ComplexF32, Cfloat , $(named_param_pairs...))
        @blis_ccall_typed($funcname, $ret, ComplexF64, Cdouble, $(named_param_pairs...))
    end
end


text:
##################### Without HetaSimulator ###############
using Parameters, NLopt, LikelihoodProfiler, CSV, DataFrames

pars = (
  K_t_mg = 101,       # constant of dissociation ATP-Mg {uM}  {1}
  K_th_mg = 10616,    # constant of dissociation ATPH-Mg {uM} {1}
  K_d_mg = 901,       # constant of dissociation ADP-Mg {uM}  {1}
  K_dh_mg = 38168,    # constant of dissociation ADPH-Mg {uM} {1}

  K_t_h = 0.323,     # constant of dissociation ATP-H {uM} {1}
  K_d_h = 0.445,     # constant of dissociation ADP-H {uM} {1}

  Mg_o = 0, # 1000,   # external Mg {uM} {2}
  Mg_i = 0, # 1000,   # internal Mg {uM} {2}
  pH_o = 7.2,        # external pH {--} {2}
  pH_i = 7.8,        # internal pH {--} {2}

  T_i = 5000, # 5000, # internal ATP {uM}
  T_o = 2000, # 3000, # external ATP {uM}
  D_i = 5000, # 5000, # internal ADP {uM}
  D_o = 50,          # external ADP {uM}

  fi = 6.9,          # {--} dimentionless potential
  # ! fitted parameters  see article PMID: 16239329
  del_D = 0,         # {--}
  del_T = 0.07,      # {--}
  k1_ANT = 35,       # {1/min}
  k2_ANT = 10.8,     # {1/min}
  k3_ANT = 21,       # {1/min}
  k4_ANT = 29,       # {1/min}
  K_D_o_ANT = 51,    # {uM}
  K_T_o_ANT = 57,    # {uM}
  A = 0.268, # -3.272670e-01,   # -a1 in article {--}
  B = -0.205, # 2.532980e-01,   # -a2 in article {--}
  C = 0.187 # 2.149946e-01,    # a3 in article {--}
)

function ant_model(p)

  _pars = merge(pars,p)
  @unpack K_t_mg, K_th_mg, K_d_mg, K_dh_mg, K_t_h, K_d_h, Mg_o, Mg_i, pH_o, pH_i, T_i, T_o, D_i, D_o, fi, del_D, del_T, k1_ANT, k2_ANT, k3_ANT, k4_ANT, K_D_o_ANT, K_T_o_ANT, A, B, C = _pars

  H_o = (exp(-pH_o*log(10))*1e6)
  H_i = (exp(-pH_i*log(10))*1e6)

  T_o_free = (T_o/(1+Mg_o/K_t_mg+H_o/K_t_h+Mg_o*H_o/(K_th_mg*K_t_h)))
  T_i_free = (T_i/(1+Mg_i/K_t_mg+H_i/K_t_h+Mg_i*H_i/(K_th_mg*K_t_h)))
  D_o_free = (D_o/(1+Mg_o/K_d_mg+H_o/K_d_h+Mg_o*H_o/(K_dh_mg*K_d_h)))
  D_i_free = (D_i/(1+Mg_i/K_d_mg+H_i/K_d_h+Mg_i*H_i/(K_dh_mg*K_d_h)))

  k1_ANT_fi = (k1_ANT*exp((A*(-4)+B*(-4)+C)*fi))
  k2_ANT_fi = (k2_ANT*exp((A*(-3)+B*(-4)+C)*fi))
  k3_ANT_fi = (k3_ANT*exp((A*(-4)+B*(-3)+C)*fi))
  k4_ANT_fi = (k4_ANT*exp((A*(-3)+B*(-3)+C)*fi))
  K_D_o_ANT_fi = (K_D_o_ANT*exp(3*del_D*fi))
  K_T_o_ANT_fi = (K_T_o_ANT*exp(4*del_T*fi))
  q_fi = (k3_ANT_fi*K_D_o_ANT_fi*exp(fi)/(k2_ANT_fi*K_T_o_ANT_fi))

  ch = (k1_ANT_fi*T_i_free*T_o_free*q_fi/K_T_o_ANT_fi+k2_ANT_fi*T_i_free*D_o_free*q_fi/K_D_o_ANT_fi+k3_ANT_fi*D_i_free*T_o_free/K_T_o_ANT_fi+k4_ANT_fi*D_i_free*D_o_free/K_D_o_ANT_fi)
  zn = ((T_i_free*q_fi+D_i_free)*(1+T_o_free/K_T_o_ANT_fi+D_o_free/K_D_o_ANT_fi))

  ch2 = (k2_ANT_fi*T_i_free*D_o_free*q_fi/K_D_o_ANT_fi+k4_ANT_fi*D_i_free*D_o_free/K_D_o_ANT_fi)

  x1 = (T_o/(1+Mg_o/K_t_mg+H_o/K_t_h+Mg_o*H_o/(K_th_mg*K_t_h)))
  x2 = (T_o/(1+Mg_o/K_d_mg+H_o/K_d_h+Mg_o*H_o/(K_dh_mg*K_d_h)))

  ch3 = (k1_ANT_fi*T_i_free*x1*q_fi/K_T_o_ANT_fi+k3_ANT_fi*D_i_free*x1/K_T_o_ANT_fi)
  zn3 = ((T_i_free*q_fi+D_i_free)*(1+x1/K_T_o_ANT_fi+x2/K_D_o_ANT_fi))

  Vm_T = (((k1_ANT_fi*T_i_free*q_fi+k3_ANT_fi*D_i_free)/(T_i_free*q_fi+D_i_free))*(1/(1+K_T_o_ANT_fi/K_D_o_ANT_fi)))
  ch4 = (k2_ANT_fi*T_i_free*x2*q_fi/K_D_o_ANT_fi+k4_ANT_fi*D_i_free*x2/K_D_o_ANT_fi)
  #ch_obm = (k2_ANT_fi*T_i_free*D_o_free*q_fi/K_D_o_ANT_fi-k3_ANT_fi*D_i_free*T_o_free/K_T_o_ANT_fi)
  #v_obm_ANT = (ch_obm/zn)
  #ff = (fi*26.5)
  #eff = (ch_obm/ch)

  return (
    v1_ANT = ch/zn,
    v2_ANT = ch2/zn,
    v3_ANT = ch3/zn3,
    v4_ANT = ch4/zn3,
    Vm_T = (((k1_ANT_fi*T_i_free*q_fi+k3_ANT_fi*D_i_free)/(T_i_free*q_fi+D_i_free))*(1/(1+K_T_o_ANT_fi/K_D_o_ANT_fi))),
    Vm_D = ((k2_ANT_fi*T_i_free*q_fi+k4_ANT_fi*D_i_free)/(T_i_free*q_fi+D_i_free))*(1/(1+K_D_o_ANT_fi/K_T_o_ANT_fi)),
  )
end

#sim_ant_model(par, fixed, out, fit=()) = ant_model(merge(fixed,fit))

p1 = [288.1, 162.8, 92.91, 48.34, 27.79]
data1 = [28.74, 25.03, 21.11, 16.47, 12.05]
sim1(t,p=()) = ant_model(merge((fi=0, D_o=0, D_i=0, T_i=10000, T_o=t),p))[:v1_ANT]

p2 = [298.6, 71.66, 32.63, 13.01]
data2 = [8.272, 5.063, 3.209, 1.64]
sim2(t,p=()) = ant_model(merge((fi=6.9, D_o=0, D_i=0, T_i=10000, T_o=t),p))[:v1_ANT]

p3 = [303.7, 159, 83.03, 49.1, 29.02]
data3 = [25.24, 22.46, 17.83, 13.26, 9.769]
sim3(d,p=()) = ant_model(merge((fi=0,D_i=10000,T_o=0,T_i=0,D_o=d),p))[:v1_ANT]

p4 = [271.3, 163.5, 93.89, 51.23, 30.23]
data4 = [23.67, 21.75, 18.47, 14.9, 11.77]
sim4(d,p=()) = ant_model(merge((fi=6.9, D_o=d, D_i=10000, T_i=0, T_o=0),p))[:v1_ANT]

p5 = [100.4, 63.86, 41.15, 25.4, 16.2, 10.21]
data5 = [26.95, 23.19, 19.19, 14.54, 10.83, 8.239]
sim5(d,p=()) = ant_model(merge((fi=6.9, D_o=d, D_i=5000, T_i=5000, T_o=0),p))[:v2_ANT]

p6 = [100.5, 63.14, 40.27, 24.84, 15.57, 9.76]
data6 = [24.72, 19.71, 15.42, 11.48, 8.357, 5.709]
sim6(d,p=()) = ant_model(merge((fi=6.9, D_o=d, D_i=5000, T_i=5000, T_o=100),p))[:v2_ANT]

p7 = [101.3, 62.48, 39.88, 24.64, 15.86, 9.581]
data7 = [21.19, 15.89, 11.89, 8.945, 6.297, 4.061]
sim7(d,p=()) = ant_model(merge((fi=6.9, D_o=d, D_i=5000, T_i=5000, T_o=400),p))[:v2_ANT]

p8  =[98.08, 62.69, 43.41, 25.27, 16.06, 10.26]
data8 = [10.95, 8.945, 7.297, 5.355, 4.061, 3.001]
sim8(d,p=()) = ant_model(merge((fi=0, D_o=d, D_i=5000, T_i=5000, T_o=0),p))[:v2_ANT]

p9 = [99.61, 64.95, 40.07, 24.52, 15.33, 9.467]
data9 = [9.71, 8.533, 5.944, 4.296, 3.119, 2.06]
sim9(d,p=()) = ant_model(merge((fi=0, D_o=d, D_i=5000, T_i=5000, T_o=20),p))[:v2_ANT]

p10 = [97.73, 63.2, 37.93, 25.35, 14.42, 8.128]
data10 = [7.18, 5.179, 3.59, 2.472, 1.589, 0.9416]
sim10(d,p=()) = ant_model(merge((fi=0, D_o=d, D_i=5000, T_i=5000, T_o=100),p))[:v2_ANT]

p11 = [300.5, 164.3, 91.74, 48.9, 29.3]
data11 = [34.55, 30.04, 25.67, 19.78, 14.69]
sim11(t,p=()) = ant_model(merge((fi=6.9, D_i=5000, T_i=5000, T_o=t),p))[:v4_ANT]

p12 = [149.5, 44.74, 22.63, 11.17]
data12 = [9.527, 7.127, 5.236, 3.636]
sim12(t,p=()) = ant_model(merge((fi=0, D_i=5000, T_i=5000, T_o=t),p))[:v4_ANT]

p13 = [180, 57.98, 27.78, 11.44]
data13 = [12.22, 9.236, 7.273, 4.073]
sim13(t,p=()) = ant_model(merge((fi=0, D_i=5000, T_i=5000, T_o=t),p))[:v3_ANT]

p14 = [173.2, 69.07, 30.23, 15.31]
data14 = [2.4, 1.818, 1.382, 1.018]
sim14(t,p=()) = ant_model(merge((fi=6.9, D_i=5000, T_i=5000, T_o=t),p))[:v3_ANT]

p15 = [0, 6.9]
data15 = [14.14, 2.53]
sim15(f,p=()) = ant_model(merge((fi=f, D_i=5000, T_i=5000, D_o=0, T_o=400),p))[:Vm_T]

p16 = [0, 6.9]
data16 = [10.57, 39.68]
sim16(f,p=()) = ant_model(merge((fi=f,),p))[:Vm_D]

p_vec = [p1,p2,p3,p4,p5,p6,p7,p8,p9,p10,p11,p12,p13,p14,p15,p16]
data_vec = [data1,data2,data3,data4,data5,data6,data7,data8,data9,data10,data11,data12,data13,data14,data15,data16]
sim_vec = [sim1,sim2,sim3,sim4,sim5,sim6,sim7,sim8,sim9,sim10,sim11,sim12,sim13,sim14,sim15,sim16]

function loss_component(pfit, p, sim, data)
  loss = 0.
    for i in eachindex(p)
      loss += (sim(p[i],pfit) - data[i])^2
    end
  return loss
end

function loss_func(pfit::NamedTuple)
  loss = 0.
  for i in eachindex(p_vec)
    loss += loss_component(pfit, p_vec[i], sim_vec[i], data_vec[i])
  end
  loss
end

function fit_params(
  p;
  ftol_abs = 0.0,
  ftol_rel = 1e-4,
  xtol_rel = 0.0,
  xtol_abs = 0.0,
  fit_alg = :LN_NELDERMEAD,
  maxeval = 10000
)
  obj(pfit::Vector, g::Vector) = loss_func(NamedTuple{Tuple(keys(p))}(pfit))

  opt = Opt(fit_alg, length(p))
  opt.min_objective = obj

  opt.ftol_rel = ftol_rel
  opt.ftol_abs = ftol_abs

  opt.xtol_rel = xtol_rel
  opt.xtol_abs = xtol_abs

  opt.maxeval = maxeval

  return NLopt.optimize(opt, collect(values(p)))
end

fitpars = (
  # ! fitted parameters  see article PMID: 16239329
  del_D = 0,         # {--}
  del_T = 0.07,      # {--}
  k1_ANT = 35,       # {1/min}
  k2_ANT = 10.8,     # {1/min}
  k3_ANT = 21,       # {1/min}
  k4_ANT = 29,       # {1/min}
  K_D_o_ANT = 51,    # {uM}
  K_T_o_ANT = 57,    # {uM}
  A = 0.268, # -3.272670e-01,   # -a1 in article {--}
  B = -0.205, # 2.532980e-01,   # -a2 in article {--}
  C = 0.187 # 2.149946e-01,    # a3 in article {--}
)

(minf, minx, ret) = fit_params(fitpars)

ci = [get_interval(
  minx,
  i,
  (p) -> loss_func(NamedTuple{Tuple(keys(fitpars))}(p)),
  :CICO_ONE_PASS,
  theta_bounds = [(-Inf, Inf) for v in minx], #fill((1e-10,1e10), length(p_names[:scn1])),
  scan_bounds=(-100.0, 1e6),
  #local_alg=:LD_MMA,
  scan_tol=1e-5,
  autodiff=false,
  scale =  fill(:direct, length(fitpars)),
  loss_crit = minf + 3.84
  ) for i in 1:length(fitpars)]

plb_1 = [iv.result[1].value for iv in ci]
pub_1 = [iv.result[2].value for iv in ci]

df = DataFrame(param = collect(keys(fitpars)), optim = minx, lb = plb_1, ub = pub_1)

max_stars_repo_name: insysbio/ant-mito-model

max_issues_repo_name: insysbio/ant-mito-model

max_forks_repo_name: insysbio/ant-mito-model


text:
"""
    PointCache{T} <: AbstractCache

An abstract cache subtype that contains a dictionary of points/costs and a vector
that stores the order of incumbent points.
"""
mutable struct PointCache{T} <: AbstractCache
    #TODO: Record the type of point: incumbent, violates unrelaxable constraints,
    #violates relaxable constraints for iteration i, etc.

    #TODO: Set maximum cache size

    #TODO: An alternative data structure would be much faster

    #Map a point to a cost value
    costs::Dict{Vector{T},T}
    #List the incumbent points in the order they are considered
    order::Vector{Vector{T}}

    function PointCache{T}() where T
        c = new()
        c.costs = Dict{Vector{T},T}()
        c.order = Vector{Vector{T}}()
        return c
    end
end

"""
    CachePush(p::AbstractProblem, x::Vector, cost)

Add point `x` and its cost `cost` to the cache of `p`.
"""
CachePush(p::AbstractProblem{T}, x::Vector{T}, cost::T) where T = CachePush(p.cache, x, cost)

"""
    CachePush(p::AbstractProblem)

Add the feasible and infeasible incumbent points (assuming neither are `nothing`)
to the cache.
"""
function CachePush(p::AbstractProblem)
    p.x != nothing && CachePush(p.cache, p.x, p.x_cost)
    p.i != nothing && CachePush(p.cache, p.i, p.i_cost)
end

function CachePush(c::PointCache{T}, x::Vector{T}, cost::T) where T
    push!(c.costs, x=>cost)
end

#TODO add equivilent for infeasible
"""
    CacheOrderPush(p::AbstractProblem{T}) where T

Add the feasible incumbent point to the order vector.
"""
CacheOrderPush(p::AbstractProblem{T}) where T = CacheOrderPush(p.cache, p.x)
function CacheOrderPush(c::PointCache{T},
                        x::Union{Vector{T},Nothing},
                       ) where T
    x === nothing && return

    if length(c.order) == 0 || x != c.order[end]
        push!(c.order, x)
    end
end

"""
    CacheQuery(p::AbstractProblem, x::Vector)

Query the cache of `p` to find if it has a cost value for point `x`. Alias
to `haskey`.
"""
CacheQuery(p::AbstractProblem, x::Vector) = CacheQuery(p.cache, x)
CacheQuery(c::PointCache{T}, x::Vector{T}) where T = haskey(c.costs, x)

"""
    CacheGet(p::AbstractProblem, x::Vector)

Return the cost of point `x` in the cache of `p`. Does not check if
`x` is in the cache, use `CacheQuery` to check.
"""
CacheGet(p::AbstractProblem, x::Vector) = CacheGet(p.cache, x)
(CacheGet(c::PointCache{T}, x::Vector{T})::T) where T = c.costs[x]

"""
    CacheRandomSample(p::AbstractProblem, n::Int)

Returns a uniformly sampled collection of `n` points from the cache. Points
can be repeated in the sample.
"""
CacheRandomSample(p::AbstractProblem, n::Int) = CacheRandomSample(p.cache, n)
(CacheRandomSample(c::PointCache{T}, n::Int)::Vector{Vector{T}}) where T =
    length(c.order) == 0 ? [] : rand(c.order, n)

"""
    CacheInitialPoint(p::AbstractProblem)

Return a tuple of the initial point added to the cache and its cost.
"""
CacheInitialPoint(p::AbstractProblem) = CacheInitialPoint(p.cache)
(CacheInitialPoint(c::PointCache{T})::Tuple{Vector,T}) where T = (c.order[1], c.costs[c.order[1]])

"""
    CacheGetRange(p::AbstractProblem, points::Vector)::Vector{Vector}

Return a vector of costs corresponding to the vector of points.
"""
CacheGetRange(p::AbstractProblem, points::Vector)::Vector = CacheGetRange(p.cache, points)
CacheGetRange(c::PointCache{T}, points::Vector) where T =
              CacheGetRange(c, convert(Vector{Vector{T}}, points))
(CacheGetRange(c::PointCache, points::Vector{Vector{T}})::Vector{T}) where T = map(p->CacheGet(c,p), points)

"""
    CacheFilter(p::AbstractProblem{T}, points::Vector{T})::Tuple{Vector{Vector{T}},Vector{Vector{T}}} where T

Return a tuple where the first entry is the set of input points in the cache and the
second is the set of input points not in the cache.
"""
(CacheFilter(p::AbstractProblem{T}, points
            )::Tuple{Vector{Vector{T}},Vector{Vector{T}}}) where T =
            CacheFilter(p.cache, convert(Vector{Vector{T}}, points))

(CacheFilter(p::PointCache{T}, points
            )::Tuple{Vector{Vector{T}},Vector{Vector{T}}}) where T =
            CacheFilter(p, convert(Vector{Vector{T}}, points))

function CacheFilter(c::PointCache{T}, points::Vector{Vector{T}}
                    )::Tuple{Vector{Vector{T}},Vector{Vector{T}}} where T
    qt = p -> CacheQuery(c, p)
    qf = p -> !CacheQuery(c, p)
    return filter(qt, points),filter(qf, points)
end


text:
function parse_h2_technologies!(dtr::DieterModel, df::DataFrame)
# function parse_h2_technologies!(dtr::DieterModel, path::AbstractString)
    # df = df::DataFrame.read(path)
    dtr.sets[:H2Technologies] = disallowmissing(unique(df[!,:H2Technologies]))

    dtr.sets[:P2G] = disallowmissing(unique([row[:H2Technologies] for row in eachrow(df)
        if row[:H2Type] == "P2G"]))

    dtr.sets[:G2P] = disallowmissing(unique([row[:H2Technologies] for row in eachrow(df)
        if row[:H2Type] == "G2P"]))

    dtr.sets[:GasStorages] = disallowmissing(unique([row[:H2Technologies] for row in eachrow(df)
        if row[:H2Type] == "GasStorage"]))

    params = map_idcol(df, [:Region, :H2Technologies], skip_cols=Symbol[])
    for (k,v) in params update_dict!(dtr.parameters, k, v) end

    rel_node_h2tech = create_relation(df,:Region,:H2Technologies,:Efficiency)
    dtr.data["relations"]["rel_node_h2tech"] = rel_node_h2tech

    Nodes = dtr.sets[:Nodes]
    H2Technologies = dtr.sets[:H2Technologies]
    dtr.sets[:Nodes_H2Tech] = Dieter.tuple2_filter(rel_node_h2tech, Nodes, dtr.sets[:H2Technologies])
    dtr.sets[:Nodes_P2G] = Dieter.tuple2_filter(rel_node_h2tech, Nodes, dtr.sets[:P2G])
    dtr.sets[:Nodes_G2P] = Dieter.tuple2_filter(rel_node_h2tech, Nodes, dtr.sets[:G2P])
    dtr.sets[:Nodes_GasStorages] = Dieter.tuple2_filter(rel_node_h2tech, Nodes, dtr.sets[:GasStorages])

    return nothing
end

function calc_inv_gas!(dtr::DieterModel)
    Nodes_H2Tech = dtr.sets[:Nodes_H2Tech]

    T = dtr.sets[:H2Technologies]
    oc = dtr.parameters[:H2OvernightCost]
    lt = dtr.parameters[:Lifetime]
    i = dtr.settings[:interest]

    dict = Dict((n,t) => oc[n,t]*annuity(i, lt[n,t]) for (n,t) in Nodes_H2Tech if t in T)
    update_dict!(dtr.parameters, :InvestmentCost, dict)

    fc = dtr.parameters[:FuelCost]      # FuelCost in currency/MWh-thermal-input
    eff = dtr.parameters[:Efficiency]   # Efficiency is unitless (MWh-output/MWh-thermal-input)
    vc = dtr.parameters[:VariableCost]

    marginalcost = Dict((n,t) => fc[n,t]/eff[n,t] + vc[n,t] for (n,t) in Nodes_H2Tech)

    update_dict!(dtr.parameters, :MarginalCost, marginalcost)

    return nothing
end

# # Hydrogen - code for testing sub-functions
# fileDict["h2_technologies"] = joinpath(datapath,"h2","h2_technologies.sql")
# dfDict["h2_technologies"] = parse_file(fileDict["h2_technologies"]; dataname=dataname)
#
# Dieter.parse_h2_technologies!(dtr, dfDict["h2_technologies"])
#
# rel_node_h2tech = create_relation(dfDict["h2_technologies"],:Region,:H2Technologies,:Efficiency)
# Nodes = dtr.sets[:Nodes]
# H2Technologies = dtr.sets[:H2Technologies]
# dtr.sets[:Nodes_H2Tech] = Dieter.tuple2_filter(rel_node_h2tech, Nodes, H2Technologies)
# Dieter.calc_inv_gas!(dtr)
#
# temp_h2_set = dtr.settings[:h2]
# dtr.settings[:h2] = missing
# Dieter.parse_extensions!(dtr,dataname=sql_db_path)
# dtr.settings[:h2] = temp_h2_set


text:
"""Struct holding the prognostic spectral variables."""
struct PrognosticVariables{NF<:AbstractFloat}
    vor     ::Array{Complex{NF},3}      # Vorticity of horizontal wind field
    div     ::Array{Complex{NF},3}      # Divergence of horizontal wind field
    Tabs    ::Array{Complex{NF},3}      # Absolute temperature [K]
    logp0   ::Array{Complex{NF},2}      # Log of surface pressure [log(Pa)]
    humid   ::Array{Complex{NF},3}      # Specific humidity [g/kg]
    geopot  ::Array{Complex{NF},3}      # Atmospheric geopotential [m²/s²]
end

"""Initialize prognostic variables from rest or restart from file."""
function initial_conditions(    P::Params,                      # Parameter struct
                                B::Boundaries{NF},              # Boundaries struct
                                G::GeoSpectral{NF}              # GeoSpectral struct
                                ) where {NF<:AbstractFloat}     # number format NF

    @unpack initial_conditions = P

    if initial_conditions == :rest
        ProgVars = initialize_from_rest(P,B,G)
    elseif initial_conditions == :restart
        ProgVars = initialize_from_file(P,B,G)
    else
        throw(error("Incorrect initialization option, $initial_conditions given."))
    end

    return ProgVars
end

"""Initialize a PrognosticVariables struct for an atmosphere at rest. No winds,
hence zero vorticity and divergence, but temperature, pressure and humidity are
initialised """
function initialize_from_rest(  P::Params,
                                B::Boundaries{NF},
                                G::GeoSpectral{NF}
                                ) where {NF<:AbstractFloat}

    @unpack nlev = G.geometry
    @unpack mx, nx = G.spectral
    @unpack ϕ0 = B

    # conversion to type NF later when creating a PrognosticVariables struct
    vor    = zeros(Complex{Float64}, mx, nx, nlev)  # Vorticity
    div    = zeros(Complex{Float64}, mx, nx, nlev)  # Divergence
    Tabs   = zeros(Complex{Float64}, mx, nx, nlev)  # Absolute Temperature
    logp0  = zeros(Complex{Float64}, mx, nx)        # logarithm of surface pressure
    humid  = zeros(Complex{Float64}, mx, nx, nlev)  # specific humidity
    geopot = zeros(Complex{Float64}, mx, nx, nlev)  # geopotential

    initialize_temperature!(Tabs,P,B,G)             # temperature from lapse rates
    logp0_grid = initialize_pressure!(logp0,P,B,G)  # pressure from temperature profile
    initialize_humidity!(humid,logp0_grid,P,G)      # specific humidity from pressure
    geopotential!(geopot,ϕ0spectral,Tabs,G)         # geopotential from surface geopotential

    # conversion to NF happens here implicitly
    return PrognosticVariables{NF}(vor, div, Tabs, logp0, humid, geopot)
end

"""Initialize spectral temperature from surface absolute temperature and constant
lapse rate (troposphere) and zero lapse rate (stratosphere)."""
function initialize_temperature!(   Tabs::AbstractArray{Complex{NF},3}, # spectral temperature in 3D
                                    P::Params,                          # Parameters struct
                                    B::Boundaries{NF},                  # Boundaries struct
                                    G::GeoSpectral{NF}                  # Geospectral struct
                                    ) where {NF<:AbstractFloat}         # number format NF

    # Compute spectral surface geopotential
    @unpack ϕ0 = B
    ϕ0_spectral = spectral(ϕ0,G)

    # Tabs_ref: Reference absolute T [K] at surface z = 0, constant lapse rate
    # Tabs_top: Reference absolute T in the stratosphere [K], lapse rate = 0
    @unpack Tabs_ref, Tabs_top, γ, g, R = P

    γ_g = γ/g/1000                  # Lapse rate [K/m] scaled by gravity
    T0 = -γ_g*ϕ0_spectral            # Surface air temperature
    T0[1,1] += Complex(√2*Tabs_ref) # adjust mean value (spectral coefficient 1,1) with Tabs_ref

    # Stratosphere, set the first spectral coefficient (=mean value)
    # in two uppermost levels (i.e. k=1,2) for lapse rate = 0
    #TODO why √2? (normalisation?)
    Tabs[1,1,1] = Complex(√2*Tabs_top)
    Tabs[1,1,2] = Complex(√2*Tabs_top)

    # Temperature at tropospheric levels
    @unpack σ_full = G.geometry

    for k in 3:nlev
        Tabs[:,:,k] = T0*σ_full[k]^(R*γ_g)
    end
end

"""Initialize the logarithm of surface pressure `logp0` consistent with temperature profile."""
function initiliaze_pressure!(  logp0::AbstractArray{Complex{NF},2},    # logarithm of surface pressure
                                P::Params,                              # Parameters struct
                                B::Boundaries{NF},                      # Boundaries struct
                                G::GeoSpectral{NF}                      # Geospectral struct
                                ) where {NF<:AbstractFloat}             # number format NF

    @unpack nlon, nlat, p0_ref = P
    @unpack Tabs_ref, Tabs_top, γ, g, R = P
    @unpack ϕ0 = B                  # load surface geopotential

    γ_g = γ/g/1000                  # Lapse rate [K/m] scaled by gravity
    logp0_ref = log(p0_ref)         # logarithm of reference surface pressure
    logp0_grid = zeros(nlon, nlat)  # logarithm of surface pressure by grid point

    for j in 1:nlat
        for i in 1:nlon
            logp0_grid[i,j] = logp0_ref + log(1 - γ_g*ϕ0[i,j]/Tabs_ref)/(R*γ_g)
        end
    end

    spectral!(logp0,logp0_grid,G)   # convert to spectral space
    spectral_truncation!(logp0,G)   # smooth in spectral space

    return logp0_grid
end

"""Initialize specific humidity in spectral space."""
function initialize_humidity!(  humid::AbstractArray{Complex{NF},3},    # specific humidity
                                logp0_grid::AbstractArray{NF,2},        # logarithm of surface pressure
                                P::Params,                              # Parameters struct
                                G::GeoSpectral{NF}                      # Geospectral struct
                                ) where {NF<:AbstractFloat}             # number format NF

    mx,nx,nlev = size(humid)
    @unpack nlon, nlat = P
    @boundscheck nlon, nlat == size(logp0_grid) || throw(BoundsError())

    @unpack es_ref, rh_ref = P      # reference saturation water vapour pressure [Pa] + relative humidity [1]
    @unpack hscale, hshum = P       # scale height [km], scale height for spec humidity [km]
    q_ref = rh_ref*0.622*es_ref     # reference specific humidity [Pa]
    q_exp = hscale/hshum            # ratio of scale heights [1]

    # Specific humidity at the surface (grid space)
    humid0_grid = q_ref*exp.(q_exp*logp0_grid)

    humid0 = spectral(humid0_grid,G)
    spectral_truncation!(humid0,G)

    # Specific humidity at tropospheric levels (zero in the stratosphere[?])
    for k in 3:nlev
        for j in 1:nx
            for i in 1:mx
                humid[i,j,k] = humid0[i,j]*σ_full[k]^q_exp
            end
        end
    end
end


text:
import Random.shuffle!
@testset "fsm_reorder.jl" begin
    DEFAULT_CAPACITY = 65000
    TIMEOUT_DFLT = 1000
    NREPS = 64
    MAX_PAYLOAD_SIZE = 1460
    @testset "non-overlapping out-of-order segments" begin
        for i in 1:NREPS
            rx_isn, tx_isn = WrappingInt32(rand(UInt32)), WrappingInt32(rand(UInt32))
            conn = fsm_in_established(tx_isn, rx_isn; cap=DEFAULT_CAPACITY)
            seq_size = Tuple{UInt, UInt}[]
            datalen = 0
            while datalen < DEFAULT_CAPACITY
                sz = min(rand(UInt) % (MAX_PAYLOAD_SIZE - 1) + 1, DEFAULT_CAPACITY - datalen)
                push!(seq_size, (datalen, sz))
                datalen += sz
            end
            shuffle!(seq_size)
            d = rand('a':'z', datalen) |> join
            min_expect_ackno = rx_isn + 1
            max_expect_ackno = rx_isn + 1
            for (off, sz) in seq_size
                send_data!(conn, rx_isn + 1 + off, tx_isn + 1, @view(d[1+off:off+sz]))
                (off == min_expect_ackno.val) && (min_expect_ackno = min_expect_ackno + sz)
                max_expect_ackno = max_expect_ackno + sz
                seg = expect_seg(conn; ack=true)
                seg_hdr = seg.header
                @test seg_hdr.ackno.val >= min_expect_ackno.val
                @test seg_hdr.ackno.val <= max_expect_ackno.val
            end
            tick!(conn, 1)
            expect_data!(conn, d)
        end
    end

    @testset "overlapping out-of-order segments" begin
        for i in 1:NREPS
            rx_isn, tx_isn = WrappingInt32(rand(UInt32)), WrappingInt32(rand(UInt32))
            conn = fsm_in_established(tx_isn, rx_isn; cap=DEFAULT_CAPACITY)
            seq_size = Tuple{UInt, UInt}[]
            datalen = 0
            while datalen < DEFAULT_CAPACITY
                sz = min(rand(UInt) % (MAX_PAYLOAD_SIZE - 1) + 1, DEFAULT_CAPACITY - datalen)
                rem = MAX_PAYLOAD_SIZE - sz
                offs = UInt(0)
                if rem == 0
                    offs = 0
                elseif rem == 1
                    offs = min(UInt(1), datalen)
                else
                    offs = min(datalen, rem, 1 + rand(UInt) % (rem - 1))
                end
                @test sz + offs <= MAX_PAYLOAD_SIZE
                push!(seq_size, (datalen - offs, sz + offs))
                datalen += sz
            end
            @test datalen <= DEFAULT_CAPACITY
            shuffle!(seq_size)
            d = rand('a':'z', datalen) |> join
            min_expect_ackno = rx_isn + 1
            max_expect_ackno = rx_isn + 1
            for (off, sz) in seq_size
                send_data!(conn, rx_isn + 1 + off, tx_isn + 1, @view(d[1+off:off+sz]))
                (off <= min_expect_ackno.val && off + sz > min_expect_ackno.val) &&
                    (min_expect_ackno = WrappingInt32(sz + off))
                max_expect_ackno = max_expect_ackno + sz

                seg = expect_seg(conn; ack=true)
                seg_hdr = seg.header
                @test seg_hdr.ackno.val >= min_expect_ackno.val
                @test seg_hdr.ackno.val <= max_expect_ackno.val
            end
            tick!(conn, 1)
            expect_data!(conn, d)
        end
    end

end

max_stars_repo_name: AquaIndigo/JLSponge

max_issues_repo_name: AquaIndigo/JLSponge

max_forks_repo_name: AquaIndigo/JLSponge


text:
"""
module jInvVisPyPlot

Visualization tools for jInv

"""
module jInvVisPyPlot
    using jInv.Mesh
    using PyPlot
    include("plotModel.jl");
end

max_stars_repo_name: JuliaInv/jInvVisPyPlot

max_issues_repo_name: JuliaInv/jInvVisPyPlot

max_forks_repo_name: JuliaInv/jInvVisPyPlot


text:
# Testing a pretrained normalizing flow for seismic interpolation
# Authors: Rajiv Kumar
#          Ali Siahkoohi, alisk@gatech.edu
# Date: February 2021
"""
Example:
julia test_hint_interpolation.jl
"""

using DrWatson
@quickactivate :WavefieldRecoveryUQ

using InvertibleNetworks
using JLD2
using HDF5
using Random
using Statistics
using ArgParse
using ProgressMeter
using PyPlot
using Seaborn
using LinearAlgebra
using Flux: gpu, cpu
using Flux.Data: DataLoader

set_style("whitegrid")
rc("font", family="serif", size=13)
font_prop = matplotlib.font_manager.FontProperties(
    family="serif",
    style="normal",
    size=13
)
sfmt=matplotlib.ticker.ScalarFormatter(useMathText=true)
sfmt.set_powerlimits((0, 0))
matplotlib.use("Agg")

s = ArgParseSettings()
@add_arg_table! s begin
    "--max_epoch"
        help = "Maximum # of epochs"
        arg_type = Int
        default = 12
    "--epoch"
        help = "Epoch to load net params"
        arg_type = Int
        default = -1
    "--lr"
        help = "starting learning rate"
        arg_type = Float32
        default = 1f-3
    "--lr_step"
        help = "lr scheduler applied at evey lr_step epoch"
        arg_type = Int
        default = 1
    "--batchsize"
        help = "batch size"
        arg_type = Int
        default = 8
    "--n_hidden"
        help = "# hidden channels"
        arg_type = Int
        default = 64
    "--depth"
        help = "depth of the network"
        arg_type = Int
        default = 8
    "--subsample_fac"
        help = "subsampling factor"
        arg_type = Float32
        default = 7.5f-1
    "--sim_name"
        help = "simulation name"
        arg_type = String
        default = "wavefield_recovery_and_UQ"
end
parsed_args = parse_args(s)

max_epoch = parsed_args["max_epoch"]
lr = parsed_args["lr"]
lr_step = parsed_args["lr_step"]
batchsize = parsed_args["batchsize"]
n_hidden = parsed_args["n_hidden"]
depth = parsed_args["depth"]
subsample_fac = parsed_args["subsample_fac"]
sim_name = parsed_args["sim_name"]

if parsed_args["epoch"] == -1
    parsed_args["epoch"] = parsed_args["max_epoch"]
end
epoch = parsed_args["epoch"]

# Define raw data directory
mkpath(datadir("dataset"))
data_path = datadir("dataset", "frequency_slices.h5")

# Download the dataset into the data directory if it does not exist
if isfile(data_path) == false
    run(`wget https://www.dropbox.com/s/wcqucrlye615fhi/'
        'frequency_slices.h5 -q -O $data_path`)
end

# Loading the experiment—only network weights and training loss
Params, fval, exp_path = load_experiment(parsed_args; return_path=true)

mask = wload(exp_path)["M.M"]
AN_params = wload(exp_path)["AN_params"]
train_idx = wload(exp_path)["train_idx"]
test_idx = wload(exp_path)["test_idx"]

# Load seismic images and create training and testing data
h5open(data_path, "r") do file
    X_orig = file["dataset"][:, :, :, 1:2:end]

    X_orig = permutedims(X_orig, (2, 3, 1, 4))
    global nx, ny, nc, nsamples = size(X_orig)

    # Whiten data
    global AN = ActNorm(nsamples)
    put_params!(AN, convert(Array{Any,1}, AN_params))
    X_orig = AN.forward(X_orig)

    # Forward operator
    global M = Mask(mask, subsample_fac)

    global idx = shuffle(test_idx)[1]

    global X_fixed = X_orig[:, :, :, idx:idx]
    global Y_fixed = M*X_orig[:, :, :, idx:idx]

end

# Create network
CH = NetworkConditionalHINT(
    nc, n_hidden, depth, k1=3, k2=3, p1=1, p2=1
)

# Loading the experiment—only network weights and training loss
put_params!(CH, convert(Array{Any,1}, Params))

# Now select single fixed sample from all Ys
Zy_fixed = CH.forward_Y(Y_fixed)

# Draw new Zx, while keeping Zy fixed
n_samples = 1000
X_post = zeros(Float32, nx, ny, nc, n_samples)
CH = CH |> gpu

test_loader = DataLoader(
    (randn(Float32, nx, ny, nc, n_samples), repeat(Zy_fixed, 1, 1, 1, n_samples)),
    batchsize=8, shuffle=false
)

p = Progress(length(test_loader))
for (itr, (X, Y)) in enumerate(test_loader)
    Base.flush(Base.stdout)
    counter = (itr - 1)*size(X)[4] + 1

    X = X |> gpu
    Y = Y |> gpu

    X_post[:, :, :, counter:(counter+size(X)[4]-1)] = (CH.inverse(X, Y)[1] |> cpu)
    ProgressMeter.next!(p)
end

X_post = AN.inverse(X_post)
X_fixed = AN.inverse(X_fixed)
Y_fixed = AN.inverse(Y_fixed)

# Some stats
X_post_mean = mean(X_post; dims=4)
X_post_std = std(X_post; dims=4)

min_std=1.5f0*minimum(X_post_std)
max_std=0.93f0*maximum(X_post_std)

X_post = (1 .- M.M) .* X_post .+ (M.M .* Y_fixed)
X_post_mean = (1 .- M.M) .* X_post_mean .+ (M.M .* Y_fixed)
X_post_std = (1 .- M.M) .* X_post_std

# Some stats
X_post_mean = mean(X_post; dims=4)
X_post_std = std(X_post; dims=4)

save_dict = @strdict max_epoch epoch lr lr_step batchsize n_hidden depth sim_name
save_path = plotsdir(sim_name, savename(save_dict; digits=6))

snr_list = []
signal_to_noise(x, xhat) = (-20.0 * log(norm(x - xhat)/norm(x))/log(10.0))

for j = 1:n_samples
    push!(snr_list, signal_to_noise(X_post[:, :, :, j], X_fixed[:, :, :, 1]))
end
X_post_mean_snr = signal_to_noise(X_post_mean[:, :, :, 1], X_fixed[:, :, :, 1])

# Training loss
fig = figure("training logs", figsize=(7, 2.5))
if epoch == parsed_args["max_epoch"]
    plot(range(0, epoch, length=length(fval)), fval, color="#d48955")
else
    plot(
        range(0, epoch, length=length(fval[1:findfirst(fval .== 0f0)-1])),
        fval[1:findfirst(fval .== 0f0)-1], color="#d48955"
    )
end
title("Negative log-likelihood")
ylabel("Training objective")
xlabel("Epochs")
safesave(joinpath(save_path, "log.png"), fig)
close(fig)

# Plot posterior samples, mean and standard deviation
fig = figure(figsize=(18, 8))
subplot(2,4,1)
imshow(
    X_fixed[:, :, 1, 1], vmin=-6.0, vmax=6.0, aspect=1, cmap="Greys", resample=true,
    interpolation="lanczos", filterrad=1
)
title(L"Fully sampled: $\mathbf{x} \sim p_x$")
colorbar(fraction=0.047, pad=0.01)
grid(false)

subplot(2,4,2)
imshow(
    X_post[:, :, 1, 1], vmin=-6.0, vmax=6.0, aspect=1, cmap="Greys", resample=true,
    interpolation="lanczos", filterrad=1
)
title(L"Conditional sample: $\widehat{\mathbf{x}} \sim \mathbf{x} \mid \mathbf{y}$")
colorbar(fraction=0.047, pad=0.01)
grid(false)

subplot(2,4,3)
imshow(
    X_post[:, :, 1, 2], vmin=-6.0, vmax=6.0, aspect=1, cmap="Greys", resample=true,
    interpolation="lanczos", filterrad=1
)
title(L"Conditional sample: $\widehat{\mathbf{x}} \sim \mathbf{x} \mid \mathbf{y}$")
colorbar(fraction=0.047, pad=0.01)
grid(false)

subplot(2,4,4)
imshow(
    X_post_mean[:, :, 1, 1], vmin=-6.0, vmax=6.0, aspect=1, cmap="Greys", resample=true,
    interpolation="lanczos", filterrad=1
)
title(L"Conditional mean: $\mathbb{E} \left [\mathbf{x} \mid \mathbf{y} \right ]$")
colorbar(fraction=0.047, pad=0.01)
grid(false)

subplot(2,4,5)
imshow(
    Y_fixed[:, :, 1, 1], vmin=-6.0, vmax=6.0, aspect=1, cmap="Greys", resample=true,
    interpolation="lanczos", filterrad=1
)
title(L"Subsampled $\mathbf{y} = \mathbf{M} \odot \mathbf{x}$ ")
colorbar(fraction=0.047, pad=0.01)
grid(false)

subplot(2,4,6)
imshow(
    X_post[:, :, 1, 4], vmin=-6.0, vmax=6.0, aspect=1, cmap="Greys", resample=true,
    interpolation="lanczos", filterrad=1
)
title(L"Conditional sample: $\widehat{\mathbf{x}} \sim \mathbf{x} \mid \mathbf{y}$")
colorbar(fraction=0.047, pad=0.01)
grid(false)

subplot(2,4,7)
imshow(
    X_post[:, :, 1, 5], vmin=-6.0, vmax=6.0, aspect=1, cmap="Greys", resample=true,
    interpolation="lanczos", filterrad=1
)
title(L"Conditional sample: $\widehat{\mathbf{x}} \sim \mathbf{x} \mid \mathbf{y}$")
colorbar(fraction=0.047, pad=0.01)
grid(false)

subplot(2,4,8)
imshow(
    X_post_std[:, :, 1,1], vmin=min_std, vmax=max_std,
    cmap="OrRd", resample=true, interpolation="lanczos", filterrad=1
)
title("Pointwise standard deviation");
colorbar(fraction=0.047, pad=0.01)
grid(false)
tight_layout()
safesave(joinpath(save_path, "posterior-"*string(test_idx[1])*".png"), fig)
close(fig)


text:
#!/usr/bin/env julia

using Luxor, Random

Random.seed!(42)

using Test

function test_circular_arrows_1(pos)
    gsave()
    froma = rescale(rand(1:100), 1, 100, 0, 2pi)
    toa =   rescale(rand(1:100), (1, 100), (0, 2pi))
    sethue("black")
    arrow(pos, 100, froma, toa, linewidth=rand(1:6), arrowheadlength=rand(10:30))
    text(string("from: ", round(rad2deg(froma), digits=1)), pos)
    text(string("to: ", round(rad2deg(toa), digits=1)), pos.x, pos.y+10)
    sethue("magenta")
    arrow(pos, 100, toa, froma; linewidth=rand(1:6), arrowheadlength=rand(10:30))
    text(string("from: ", round(rad2deg(toa), digits=1)), pos.x, pos.y+20)
    text(string("to: ", round(rad2deg(froma), digits=1)), pos.x, pos.y+30)
    grestore()
end

function test_circular_arrows_2(pos, w)
  sethue("black")
  for a in 50:10:w
    randomhue()
    starta = rand(0:pi/12:2pi)
    finisha = rand(0:pi/12:2pi)
    arrow(pos, a, starta, finisha, linewidth=rand(1:6), arrowheadlength=rand(1:20))
  end
end

function arrow_arc_test(fname)
  pagewidth, pageheight = 2000, 2000
  Drawing(pagewidth, pageheight, fname)
  origin() # move 0/0 to center
  background(1, 1, 0.9, 1)
  setopacity(0.5)
  setline(2)
  pagetiles = Tiler(pagewidth, pageheight, 4, 4, margin=50)
  for (pos, n) in pagetiles
    if isodd(n)
      test_circular_arrows_1(pos)
    else
      test_circular_arrows_2(pos, pagetiles.tilewidth/2)
    end
  end

  # test Bezier arrow
  setopacity(1.0)
  sethue("red")
  arrow(pagetiles[1][1], pagetiles[3][1], pagetiles[6][1], pagetiles[7][1],
    linewidth=5,
    arrowheadlength=50,
    arrowheadfill = false,
    arrowheadangle=π/4,
    startarrow=true)

  @test finish() == true
  println("...finished arrow-test: saved in $(fname)")
end

fname = "arrow-arctest.png"
arrow_arc_test(fname)


text: bibrefs() = _bibrefs(bibtex_paths())