text:
"""
    HorizontalStructure

Horizontal memory structure. Decomposes stochastic program by scenario.

"""
struct HorizontalStructure{N, M, SP <: NTuple{M, AbstractScenarioProblems}} <: AbstractBlockStructure{N}
    decisions::NTuple{M, Decisions}
    proxy::JuMP.Model
    scenarioproblems::SP

    function HorizontalStructure(scenarioproblems::NTuple{M,AbstractScenarioProblems}) where M
        N = M + 1
        decisions = ntuple(Val(M)) do i
            Decisions()
        end
        SP = typeof(scenarioproblems)
        return new{N,M,SP}(decisions, Model(), scenarioproblems)
    end
end

function StochasticStructure(scenario_types::ScenarioTypes{M}, instantiation::Union{Horizontal, DistributedHorizontal}) where M
    scenarioproblems = ntuple(Val(M)) do i
        ScenarioProblems(scenario_types[i], instantiation)
    end
    return HorizontalStructure(scenarioproblems)
end

function StochasticStructure(scenarios::NTuple{M, Vector{<:AbstractScenario}}, instantiation::Union{Horizontal, DistributedHorizontal}) where M
    scenarioproblems = ntuple(Val(M)) do i
        ScenarioProblems(scenarios[i], instantiation)
    end
    return HorizontalStructure(scenarioproblems)
end

# Base overloads #
# ========================== #
function Base.print(io::IO, structure::HorizontalStructure{2})
    print(io, "Horizontal scenario problems \n")
    print(io, "============== \n")
    for (id, subproblem) in enumerate(subproblems(structure))
        @printf(io, "Subproblem %d (p = %.2f):\n", id, probability(scenario(structure, id)))
        print(io, subproblem)
        print(io, "\n")
    end
end

# MOI #
# ========================== #
function MOI.get(structure::HorizontalStructure, attr::MOI.AbstractModelAttribute)
    return MOI.get(backend(structure.proxy), attr)
end
function MOI.get(structure::HorizontalStructure, attr::MOI.AbstractVariableAttribute, index::MOI.VariableIndex)
    return MOI.get(backend(structure.proxy), attr, index)
end
function MOI.get(structure::HorizontalStructure, attr::Type{MOI.VariableIndex}, name::String)
    return MOI.get(backend(structure.proxy), attr, name)
end
function MOI.get(structure::HorizontalStructure, attr::MOI.AbstractConstraintAttribute, cindex::MOI.ConstraintIndex)
    return MOI.get(backend(structure.proxy), attr, cindex)
end

function MOI.set(structure::HorizontalStructure, attr::Union{MOI.AbstractModelAttribute, MOI.Silent}, value)
    MOI.set(scenarioproblems(structure), attr, value)
    MOI.set(backend(structure.proxy), attr, value)
end
function MOI.set(structure::HorizontalStructure, attr::MOI.AbstractVariableAttribute,
                 index::MOI.VariableIndex, value)
    MOI.set(scenarioproblems(structure), attr, index, value)
    MOI.set(backend(structure.proxy), attr, index, value)
    return nothing
end
function MOI.set(structure::HorizontalStructure, attr::MOI.AbstractConstraintAttribute,
                 cindex::MOI.ConstraintIndex, value)
    MOI.set(scenarioproblems(structure), attr, index, value)
    MOI.set(backend(structure.proxy), attr, cindex, value)
    return nothing
end

function MOI.is_valid(structure::HorizontalStructure, index::MOI.VariableIndex)
    return MOI.is_valid(backend(structure.proxy), index)
end

function MOI.add_constraint(structure::HorizontalStructure, f::MOI.AbstractFunction, s::MOI.AbstractSet)
    MOI.add_constraint(scenarioproblems(structure), f, s)
    return MOI.add_constraint(backend(structure.proxy), f, s)
end

function MOI.delete(structure::HorizontalStructure, index::MOI.Index)
    # TODO: more to do if index is decision
    MOI.delete(scenarioproblems(structure), index)
    MOI.delete(backend(structure.proxy), index)
    return nothing
end

# Getters #
# ========================== #
function structure_name(structure::HorizontalStructure)
    return "Horizontal"
end

# Setters #
# ========================== #
function untake_decisions!(structure::HorizontalStructure{2,1,NTuple{1,SP}}) where SP <: ScenarioProblems
    if untake_decisions!(structure.decisions[1])
        update_decisions!(scenarioproblems(structure), DecisionsStateChange())
    end
    return nothing
end
function untake_decisions!(structure::HorizontalStructure{2,1,NTuple{1,SP}}) where SP <: DistributedScenarioProblems
    sp = scenarioproblems(structure)
    @sync begin
        for (i,w) in enumerate(workers())
            @async remotecall_fetch(
                w, sp[w-1], sp.decisions[w-1]) do sp, d
                    if untake_decisions!(fetch(d))
                        update_decisions!(fetch(sp), DecisionsStateChange())
                    end
                end
        end
    end
    return nothing
end


text:
module Unwrapping

using FFTW

"""
    twopi(T)

yields the value of `2π` with floating-point type `T`.

"""
twopi(T::Type{<:AbstractFloat}) = 2*T(pi)

"""
    arc(x)

yields `x` wrapped in the range `[-π,+π]`.

"""
arc(x::Real) = arc(Float64(x))
arc(x::T) where {T<:AbstractFloat} = wrap(x, twopi(T))
# NOTE: Using `wrap(x, 2π)` is faster than `mod2pi(x + π) - π`.

"""
    wrap(x, p)

yields `x` wrapped in the range `[-p/2,+p/2]`, argument `p` is the period.  The
result is of floating-point type.

"""
wrap(x::Real, p::Real) = wrap(promote(x, p)...)
wrap(x::T, p::T) where {T<:Real} = wrap(Float64(x), Float64(p))
wrap(x::T, p::T) where {T<:AbstractFloat} = x - p*round(x/p)

"""
    unwrap(psi, p=2π; debiased=false) -> phi

yields the unwrapped phase `phi` that best fits, in a least squares sense, the
wrapped phase `psi`.  Optional argument `p` is the period of the wrapping.
The result has an undetermined constant offset.

The method is based on algorithm *"Robust two-dimensional weighted and
unweighted phase unwrapping that uses fast transforms and iterative methods"*
by Dennis C. Ghiglia and Louis A. Romero (J. Opt. Soc. Am. A, Vol. 11,
pp. 107-117, 1994).

Keyword `debiased` (true by default) specifies whether a constant bias (modulo
the period) should be avoided.  if `debiased` is false (the default), a
zero-mean solution is returned (up to rounding errors) as in the original
algorithm.  Estimating the bias has a cost, so do not set `debiased = true` if
this bias is irrelevant for your needs.

"""
function unwrap(psi::AbstractMatrix{T},
                period::Real; kwds...) where {T<:AbstractFloat}
    unwrap(psi, T(period); kwds...)
end

function unwrap(psi::AbstractMatrix{T},
                period::T = twopi(T);
                debiased::Bool = false) where {T<:AbstractFloat}
    Base.has_offset_axes(psi) && error("input array has non-standard axes")
    n1, n2 = size(psi)

    # We need 2 workspace arrays.
    wrk1 = Array{T}(undef, n1, n2)
    wrk2 = Array{T}(undef, n1, n2)

    # Compute wrapped phase Laplacian by finite differences.  We first compute
    # wrapped phase gradient along a given dimension by finite differences and
    # according to Eq. (2) or Eq. (3).  We then compute the curvature along
    # each dimension by finite differences of the gradient (in-place).  We
    # finally adds the curvatures to form the Laplacian.
    d1 = wrk1
    @inbounds for j2 in 1:n2
        @simd for j1 in 1:n1-1
            d1[j1,j2] = wrap(psi[j1+1,j2] - psi[j1,j2], period)
        end
        let j1 = n1
            d1[j1,j2] = 0
        end
        @simd for j1 in n1:-1:2
            d1[j1,j2] -= d1[j1-1,j2]
        end
    end
    d2 = wrk2
    @inbounds for j2 in 1:n2-1
        @simd for j1 in 1:n1
            d2[j1,j2] = wrap(psi[j1,j2+1] - psi[j1,j2], period)
        end
    end
    let j2 = n2
        @inbounds @simd for j1 in 1:n1
            d2[j1,j2] = 0
        end
    end
    @inbounds for j2 in n2:-1:2
        @simd for j1 in 1:n1
            d2[j1,j2] -= d2[j1,j2-1]
        end
    end
    rho = wrk1 # we re-use one of the workspaces to spare allocations
    @inbounds @simd for j1 in eachindex(wrk1, wrk2)
        rho[j1] = d1[j1] + d2[j1]
    end

    # Compute Q, the inverse of the DCT filter, that is the denominator in
    # Eq. (13) but set Q to 1 at the zero-th frequency to avoid division by
    # zero.  The zero-th frequency being undetermined is a direct consequence
    # of the phase being reconstructed from phase diffeences, it known up to a
    # constant offset.
    q = wrk2 # we re-use the other workspace to spare allocations
    q1 = Array{T}(undef, n1) # array of cosine along first dimension
    a1 = T(π)/n1
    @inbounds @simd for j1 in 1:n1
        q1[j1] = two_cos_minus_two(a1*(j1 - 1))
    end
    if n1 == n2
        @inbounds for j2 in 1:n2
            q2 = q1[j2]
            @simd for j1 in 1:n1
                q[j1,j2] = q1[j1] + q2
            end
        end
    else
        a2 = T(π)/n2
        @inbounds for j2 in 1:n2
            q2 = two_cos_minus_two(a2*(j2 - 1))
            @simd for j1 in 1:n1
                q[j1,j2] = q1[j1] + q2
            end
        end
    end
    q[1,1] = 1 # to avoid division by zero

    # Compute the solution (Julia's dct is of type II).
    dct!(rho)
    rho[1,1] = 0 # should be exactly zero without rounding errors
    phi = rho # we re-use the same workspace to spare allocations
    @inbounds @simd for j in eachindex(phi, rho, q)
        phi[j] = rho[j]/q[j]
    end
    idct!(phi)

    # Fix the piston (modulo the period) in the computed solution.
    debiased && fix_bias!(phi, psi, period)

    return phi
end

"""
    fft_unwrap(psi, p=2π; debiased=false) -> phi

yields the same result as `unwrap` but using a Fast Fourier Transform (FFT)
instead of a Cosine Fourier Transform (DCT).  This version is slower then the
one using the DCT and is provided as an example of how to use the FFT instead
of the DCT.

"""
function fft_unwrap(psi::AbstractMatrix{T},
                    period::T = twopi(T);
                    debiased::Bool = false) where {T<:AbstractFloat}
    Base.has_offset_axes(psi) && error("input array has non-standard axes")
    n1, n2 = size(psi)

    # We need 2 workspace arrays.
    wrk1 = Array{T}(undef, n1, n2)
    wrk2 = Array{T}(undef, n1, n2)

    # The mean level and mean slope (piston and tip-tilt) are undetermined when
    # the FFT is used instead of the DCT.  We therefore compute the discrete
    # wrapped gradient with the mean slope subtracted and then compute the
    # discrete curvature as before to estimate the discrete Laplacian.  After
    # reconstruction, we add the mean slope and the estimated mean level
    # (modulo the period).

    # Fist pass to estimate the mean wrapped slope along each dimension.
    s1 = zero(T)
    if n1 ≥ 2 && n2 ≥ 1
        @inbounds for j2 in 1:n2
            @simd for j1 in 1:n1-1
                s1 += wrap(psi[j1+1,j2] - psi[j1,j2], period)
            end
        end
        s1 /= (n1 - 1)*n2
    end
    s2 = zero(T)
    if n1 ≥ 1 && n2 ≥ 2
        @inbounds for j2 in 1:n2-1
            @simd for j1 in 1:n1
                s2 += wrap(psi[j1,j2+1] - psi[j1,j2], period)
            end
        end
        s2 /= n1*(n2 - 1)
    end

    # Compute the wrapped gradient (with mean slope subtracted) and the
    # curvature along each dimension.
    d1 = wrk1
    @inbounds for j2 in 1:n2
        @simd for j1 in 1:n1-1
            d1[j1,j2] = wrap(psi[j1+1,j2] - psi[j1,j2] - s1, period)
        end
        let j1 = n1
            d1[j1,j2] = 0
        end
        @simd for j1 in n1:-1:2
            d1[j1,j2] -= d1[j1-1,j2]
        end
    end
    d2 = wrk2
    @inbounds for j2 in 1:n2-1
        @simd for j1 in 1:n1
            d2[j1,j2] = wrap(psi[j1,j2+1] - psi[j1,j2] - s2, period)
        end
    end
    let j2 = n2
        @inbounds @simd for j1 in 1:n1
            d2[j1,j2] = 0
        end
    end
    @inbounds for j2 in n2:-1:2
        @simd for j1 in 1:n1
            d2[j1,j2] -= d2[j1,j2-1]
        end
    end

    # Compute the discrete Laplacian.
    rho = wrk1 # we re-use one of the workspaces to spare allocations
    @inbounds @simd for j1 in eachindex(wrk1, wrk2)
        rho[j1] = d1[j1] + d2[j1]
    end

    # Compute Q, the FFT filter implementing the discrete Laplacian.
    q = wrk2 # we re-use the other workspace to spare allocations
    q1 = Array{T}(undef, n1) # array of cosine along first dimension
    a1 = twopi(T)/n1
    @inbounds @simd for j1 in 1:n1
        q1[j1] = two_cos_minus_two(a1*(j1 - 1))
    end
    if n1 == n2
        @inbounds for j2 in 1:n2
            q2 = q1[j2]
            @simd for j1 in 1:n1
                q[j1,j2] = q1[j1] + q2
            end
        end
    else
        a2 = twopi(T)/n2
        @inbounds for j2 in 1:n2
            q2 = two_cos_minus_two(a2*(j2 - 1))
            @simd for j1 in 1:n1
                q[j1,j2] = q1[j1] + q2
            end
        end
    end
    q[1,1] = 1 # to avoid division by zero

    # Compute the solution.
    z = fft(rho)
    z[1,1] = 0 # should be exactly zero without rounding errors
    @inbounds @simd for j in eachindex(z, q)
        z[j] /= q[j]
    end
    ifft!(z)

    # Extract real part and restore the (centered) mean slopes.
    phi = wrk1
    c1 = T(n1 + 1)/2 # offset to have the slope centered
    c2 = T(n2 + 1)/2 # offset to have the slope centered
    @inbounds for j2 in 1:n2
        r2 = (j2 - c2)*s2
        @simd for j1 in 1:n1
            r1 = (j1 - c1)*s1
            phi[j1,j2] = z[j1,j2].re + (r1 + r2)
        end
    end

    # Fix the piston (modulo the period) in the computed solution.
    debiased && fix_bias!(phi, psi, period)

    return phi
end

"""
    fix_bias!(x, y, p=2π) -> x

fixes a constant bias modulo the period `p` in the array `x` which is an
unwrapped but biased version of `y`.  The operation is done in-place and `x` is
returned.

"""
function fix_bias!(x::AbstractArray{T,N},
                   y::AbstractArray{T,N},
                   period::T = twopi(T)) where {T<:AbstractFloat,N}
    a = twopi(T)/period
    sum_cos = zero(T)
    sum_sin = zero(T)
    @inbounds @simd for j in eachindex(x, y)
        sin_j, cos_j = sincos(a*(y[j]- x[j]))
        sum_sin += sin_j
        sum_cos += cos_j
    end
    bias = atan(sum_sin, sum_cos)/a
    @inbounds @simd for j in eachindex(x)
        x[j] += bias
    end
    return x
end

"""
    two_cos_minus_two(x) -> 2*cos(x) - 2

computes accurately `2*cos(x) - 2` given the floating-point value `x`.

"""
two_cos_minus_two(x::AbstractFloat) = -4*sin(x/2)^2

end # module

ext: jl


text:
function event_driven(cp::ContactProcess{SIR}, x₀; tmax=100.0)
    N = nv(cp.g)
    β = cp.dynamics.β
    γ = cp.dynamics.δ
    events = PriorityQueue{EpidemicEvent,Float64}()
    # maximum 2 events per node, so 2*N events max
    ts = Vector{Float64}(undef, 2*N)
    S = Vector{Int32}(undef, 2*N)
    I = Vector{Int32}(undef, 2*N)
    R = Vector{Int32}(undef, 2*N)
    k = 1
    ts[k] = 0.0
    S[k] = N
    I[k] = 0
    R[k] = 0
    state = copy(x₀)
    for u in 1:N
        if x₀[u] == 2
            events[InfectionEvent(u)] = 0.0
        end
    end
    while !isempty(events)
        (e,t) = dequeue_pair!(events)
        if t != ts[k]
            k += 1
            ts[k] = t
            S[k] = S[k-1]
            I[k] = I[k-1]
            R[k] = R[k-1]
        end
        if typeof(e) == InfectionEvent
            S[k] -= 1
            I[k] += 1
            u = e.node
            state[u] = 2
            if γ != 0
                τᵣ = log(1.0/rand())/γ
                events[RecoveryEvent(u)] = t+τᵣ
            else
                τᵣ = Inf
            end
            for v in filter(x->state[x]==1, neighbors(cp.g, u))
                τᵢ = log(1.0/rand())/β
                if τᵢ >= τᵣ
                    continue
                end
                if haskey(events, InfectionEvent(v))
                    if t+τᵢ >= events[InfectionEvent(v)]
                        continue
                    end
                end
                events[InfectionEvent(v)] = t + τᵢ
            end
        else # recovery
            I[k] -= 1
            R[k] += 1
        end
    end
end


text:
include("holtslag_psim.jl")
include("holtslag_psis.jl")
include("constants.jl")

function VelocityScale(z,h,L,ustar,θv=nothing,wθv=nothing)
  xi=z/L
  psis=holtslag_psis(xi)
  psim=holtslag_psim(xi)
  if xi>0
      wt=ustar/psis
      wm=wt
  else
      if z/h<0.1
          wt=ustar/psis
          wm=ustar/psim
      else
          c1=0.6
          a=7.2
          wθv<0 ? println("wθv",wθv) : nothing
          #println("θv",θv)
          wstar=((gravity_mks/θv)*wθv*h)^(1/3)
          wm=(ustar^3+c1*wstar^3)^(1/3)
          xi1=0.1*z/L
          psis1=holtslag_psis(xi1)
          psim1=holtslag_psim(xi1)
          Pr=psis1/psim1*abs(xi1)+a*karman*(0.1*(wstar/wm))
          wt=wm/Pr
      end
  end
  return wt,wm
end


text:
# values of specification fields
@enum Alignment::UInt8 ALIGN_LEFT ALIGN_RIGHT ALIGN_CENTER
@enum Sign::UInt8 SIGN_PLUS SIGN_MINUS SIGN_SPACE SIGN_NONE
@enum Grouping::UInt8 GROUPING_COMMA GROUPING_UNDERSCORE

const FILL_DEFAULT = ' '
const SIGN_DEFAULT = SIGN_MINUS
const ALTFORM_DEFAULT = false
const ZERO_DEFAULT = false

const ALIGN_UNSPECIFIED = nothing
const WIDTH_UNSPECIFIED = nothing
const GROUPING_UNSPECIFIED = nothing
const PRECISION_UNSPECIFIED = nothing
const TYPE_UNSPECIFIED = nothing

isspecified(val) = val !== nothing
default(val, defval) = isspecified(val) ? val : defval

# format specification
struct Spec
    fill::Union{Char, Argument}
    align::Union{Alignment, Nothing}
    sign::Sign
    altform::Bool
    zero::Bool
    width::Union{Int, Nothing, Argument}
    grouping::Union{Grouping, Nothing}
    precision::Union{Int, Nothing, Argument}
    type::Union{Char, Nothing}
end

function Spec(
        spec::Spec;
        fill = spec.fill,
        align = spec.align,
        sign = spec.sign,
        altform = spec.altform,
        zero = spec.zero,
        width = spec.width,
        grouping = spec.grouping,
        precision = spec.precision,
        type = spec.type)
    return Spec(fill, align, sign, altform, zero, width, grouping, precision, type)
end

const SPEC_DEFAULT = Spec(
    FILL_DEFAULT,
    ALIGN_UNSPECIFIED,
    SIGN_DEFAULT,
    ALTFORM_DEFAULT,
    ZERO_DEFAULT,
    WIDTH_UNSPECIFIED,
    GROUPING_UNSPECIFIED,
    PRECISION_UNSPECIFIED,
    TYPE_UNSPECIFIED,
)


text:
@testset "basic operations" begin

# *, +, sin
x = Diffable(15)
y = Diffable(π/2)
z = x * y + sin(x)

@test backwards!(z) == 1.0
@test backwards!(y) == x.value
@test backwards!(x) == y.value + cos(x.value)

# negation
x = Diffable(4.0)
y = -x

@test backwards!(x) == -1.0

# division

# multiple instances of one variable
x = Diffable(3.0)
y = Diffable(2.0)
z = y * (x * x)
@test backwards!(z) == 1.0
@test backwards!(y) == x.value * x.value
@test backwards!(x) == y.value * 2 * x.value

# cos
x = Diffable(2.0)
z = Diffable(2.0) * cos(x)
@test backwards!(x) == -2.0 * sin(x.value)

# using a function
f(x) = Diffable(3.0) * x
x_1 = Diffable(2.0)
z = f(x_1)
@test backwards!(x_1) == 3.0

# nested functions

# identity

# power

# exponent

# relu

end #basic operations

@testset "interop with base types" begin

x = Diffable(15)
y = π/2
z = x * y + sin(x)

@test backwards!(x) == y + cos(x.value)

x = Diffable(15)
y = π/2
z = x * y + sin(x+2)

@test backwards!(x) == y + cos(x.value+2)

end #interop with base types

@testset "powers & natural log" begin

# x^n alone
a = Diffable(4.0)
y = a^3

@test backwards!(a) == 3 * 4.0^2

# b^x alone

# x^n chained

# b^x chained

# ℯ^x

x = Diffable(8.9)
y = 3 * exp(x)

@test backwards!(x) == 3 * exp(8.9)

x = Diffable(5.1)
y = exp(x^3)

@test backwards!(x) ≈ exp(5.1^3) * 3 * 5.1^2

# ln
x = Diffable(3.2)
y = log(x)

@test backwards!(x) == 1/3.2

x = Diffable(1.0)
y = exp(log(x))

@test backwards!(x) == 1.0

x = Diffable(-3.2)
@test_throws DomainError log(x)

end #powers & natural log


text:
"""
    FEMMDeforLinearESNICEModule

Formulation for the small displacement, small strain deformation
model for Nodally-Integrated Continuum Elements (NICE).

The approximation is  originally from Dohrmann et al IJNME 47 (2000).
The formulation was subsequently developed in Krysl, P. and Zhu, B.
Locking-free continuum displacement finite elements with nodal
integration, International Journal for Numerical Methods in Engineering,
76,7,1020-1043,2008.

The stabilization scheme comes from papers on energy-sampling stabilization
for mean-strain elements (Krysl).
"""
module FEMMDeforLinearESNICEModule

using FinEtools.FTypesModule: FInt, FFlt, FCplxFlt, FFltVec, FIntVec, FFltMat, FIntMat, FMat, FVec, FDataDict
import FinEtools.FENodeSetModule: FENodeSet
import FinEtools.FESetModule: FESet, FESetH8, FESetT4, manifdim, nodesperelem, gradN!
import FinEtools.IntegDomainModule: IntegDomain, integrationdata, Jacobianvolume
import FinEtools.FEMMDeforLinearBaseModule: FEMMDeforLinearAbstract
import FinEtools.DeforModelRedModule: DeforModelRed, DeforModelRed3D
import FinEtools.MatDeforModule: MatDefor
import FinEtools.MatDeforElastIsoModule: MatDeforElastIso
import FinEtools.FieldModule: ndofs, gatherdofnums!, gatherfixedvalues_asvec!, gathervalues_asvec!, gathervalues_asmat!
import FinEtools.NodalFieldModule: NodalField, nnodes
import FinEtools.CSysModule: CSys, updatecsmat!
import FinEtools.FENodeToFEMapModule: FENodeToFEMap
import FinEtools.DeforModelRedModule: nstressstrain, nthermstrain, Blmat!
import FinEtools.AssemblyModule: SysvecAssemblerBase, SysmatAssemblerBase, SysmatAssemblerSparseSymm, startassembly!, assemble!, makematrix!, makevector!, SysvecAssembler
using FinEtools.MatrixUtilityModule: add_btdb_ut_only!, complete_lt!, add_btv!, loc!, jac!, locjac!, adjugate3!
import FinEtools.FEMMDeforLinearBaseModule: stiffness, nzebcloadsstiffness, mass, thermalstrainloads, inspectintegpoints
import FinEtools.FEMMBaseModule: associategeometry!
import FinEtools.MatDeforModule: rotstressvec
import LinearAlgebra: mul!, Transpose, UpperTriangular, eigvals
At_mul_B!(C, A, B) = mul!(C, Transpose(A), B)
A_mul_B!(C, A, B) = mul!(C, A, B)
import LinearAlgebra: norm, qr, diag, dot, cond, I, cross
import Statistics: mean
import StatsBase: geomean

abstract type FEMMDeforLinearAbstractNICE <: FEMMDeforLinearAbstract end

mutable struct _NodalBasisFunctionGradients
    gradN::FFltMat
    patchconn::FIntVec
    Vpatch::FFlt
end

function _make_stabilization_material(material::M) where {M}
    ns = fieldnames(typeof(material))
    E = 0.0; nu = 0.0
    if :E in ns
        E = material.E
        if material.nu < 0.3
            nu = material.nu
        else
            nu = 0.3 + (material.nu - 0.3) / 2.0
        end
    else
        if :E1 in ns
            E = mean([material.E1, material.E2, material.E3])
            nu = min(material.nu12, material.nu13, material.nu23)
        else
            error("No clues on how to construct the stabilization material")
        end
    end
    return  MatDeforElastIso(material.mr, 0.0, E, nu, 0.0)
end

mutable struct FEMMDeforLinearESNICET4{MR<:DeforModelRed, S<:FESetT4, F<:Function, M<:MatDefor} <: FEMMDeforLinearAbstractNICE
    mr::Type{MR}
    integdomain::IntegDomain{S, F} # geometry data
    mcsys::CSys # updater of the material orientation matrix
    material::M # material object
    stabilization_material::MatDeforElastIso
    nodalbasisfunctiongrad::Vector{_NodalBasisFunctionGradients}
    ephis::Vector{FFlt}
    nphis::Vector{FFlt}
end

function FEMMDeforLinearESNICET4(mr::Type{MR}, integdomain::IntegDomain{S, F}, mcsys::CSys, material::M) where {MR<:DeforModelRed,  S<:FESetT4, F<:Function, M<:MatDefor}
    @assert mr == material.mr "Model reduction is mismatched"
    @assert (mr == DeforModelRed3D) "3D model required"
    stabilization_material = _make_stabilization_material(material)
    return FEMMDeforLinearESNICET4(mr, integdomain, mcsys, material, stabilization_material, _NodalBasisFunctionGradients[], fill(zero(FFlt), 1), fill(zero(FFlt), 1))
end

function FEMMDeforLinearESNICET4(mr::Type{MR}, integdomain::IntegDomain{S, F}, material::M) where {MR<:DeforModelRed,  S<:FESetT4, F<:Function, M<:MatDefor}
    @assert mr == material.mr "Model reduction is mismatched"
    @assert (mr == DeforModelRed3D) "3D model required"
    stabilization_material = _make_stabilization_material(material)
    return FEMMDeforLinearESNICET4(mr, integdomain, CSys(manifdim(integdomain.fes)), material, stabilization_material, _NodalBasisFunctionGradients[], fill(zero(FFlt), 1), fill(zero(FFlt), 1))
end

function centroid!(self::F, loc, X::FFltMat, conn::C) where {F<:FEMMDeforLinearESNICET4, C}
    weights = [0.250
                0.250
                0.250
                0.250]
    return loc!(loc, X, conn, reshape(weights, 4, 1))
end

function buffers1(self::FEMMDeforLinearAbstractNICE, geom::NodalField)
    fes = self.integdomain.fes
    nne = nodesperelem(fes); # number of nodes for element
    sdim = ndofs(geom);            # number of space dimensions
    mdim = manifdim(fes); # manifold dimension of the element
    # Prepare buffers
    loc = fill(zero(FFlt), 1, sdim); # quadrature point location -- buffer
    J = fill(zero(FFlt), sdim, mdim); # Jacobian matrix -- buffer
    adjJ = fill(zero(FFlt), sdim, mdim); # Jacobian matrix -- buffer
    csmatTJ = fill(zero(FFlt), mdim, mdim); # intermediate result -- buffer
    gradN = fill(zero(FFlt), nne, mdim);
    xl = fill(zero(FFlt), nne, mdim);
    lconn = collect(1:nne)
    return loc, J, adjJ, csmatTJ, gradN, xl, lconn
end

function buffers2(self::FEMMDeforLinearAbstractNICE, geom::NodalField, u::NodalField, npts::FInt)
    fes = self.integdomain.fes
    ndn = ndofs(u); # number of degrees of freedom per node
    nne = nodesperelem(fes); # number of nodes for element
    sdim = ndofs(geom);            # number of space dimensions
    mdim = manifdim(fes); # manifold dimension of the element
    nstrs = nstressstrain(self.mr);  # number of stresses
    elmatdim = ndn*nne;             # dimension of the element matrix
    # Prepare buffers
    elmat = fill(zero(FFlt), elmatdim, elmatdim);      # element matrix -- buffer
    B = fill(zero(FFlt), nstrs, elmatdim); # strain-displacement matrix -- buffer
    loc = fill(zero(FFlt), 1, sdim); # quadrature point location -- buffer
    J = fill(zero(FFlt), sdim, mdim); # Jacobian matrix -- buffer
    csmatTJ = fill(zero(FFlt), mdim, mdim); # intermediate result -- buffer
    Jac = fill(zero(FFlt), npts);
    D = fill(zero(FFlt), nstrs, nstrs); # material stiffness matrix -- buffer
    Dstab = fill(zero(FFlt), nstrs, nstrs); # material stiffness matrix -- buffer
    return dofnums, loc, J, csmatTJ, Jac, D, Dstab, elmat, B
end

function buffers3(self::FEMMDeforLinearAbstractNICE, geom::NodalField, u::NodalField)
    fes = self.integdomain.fes
    ndn = ndofs(u); # number of degrees of freedom per node
    nne = nodesperelem(fes); # number of nodes for element
    sdim = ndofs(geom);            # number of space dimensions
    mdim = manifdim(fes); # manifold dimension of the element
    nstrs = nstressstrain(self.mr);  # number of stresses
    elmatdim = ndn*nne;             # dimension of the element matrix
    # Prepare buffers
    elmat = fill(zero(FFlt), elmatdim, elmatdim);      # element matrix -- buffer
    dofnums = zeros(FInt, elmatdim); # degree of freedom array -- buffer
    B = fill(zero(FFlt), nstrs, elmatdim); # strain-displacement matrix -- buffer
    DB = fill(zero(FFlt), nstrs, elmatdim); # strain-displacement matrix -- buffer
    elvecfix = fill(zero(FFlt), elmatdim); # vector of prescribed displ. -- buffer
    elvec = fill(zero(FFlt), elmatdim); # element vector -- buffer
    gradN = fill(zero(FFlt), nne, mdim); # intermediate result -- buffer
    return dofnums, B, DB, elmat, elvec, elvecfix, gradN
end

function patchconn(fes, gl, thisnn)
    # Generate patch connectivity for a given node (thisnn)
    # from the connectivities of the finite elements attached to it.
    return vcat(collect(setdiff(Set([i for j=1:length(gl) for i in fes.conn[gl[j]]]), thisnn)), [thisnn])
end

function computenodalbfungrads(self, geom)
    # # Compute the nodal basis function gradients.
    # # Return the cell array of structures with attributes
    # %      bfun_gradients{nix}.Nspd= basis function gradient matrix
    # #        bfun_gradients{nix}.Vpatch= nodal patch volume
    # #        bfun_gradients{nix}.patchconn= nodal patch connectivity

    fes = self.integdomain.fes
    npts,  Ns,  gradNparams,  w,  pc = integrationdata(self.integdomain);
    loc, J, adjJ, csmatTJ, gradN, xl, lconn = buffers1(self, geom)

    # Get the inverse map from finite element nodes to geometric cells
    fen2fe = FENodeToFEMap(fes.conn, nnodes(geom));
    # Initialize the nodal gradients, nodal patch, and patch connectivity
    bfungrads = fill(_NodalBasisFunctionGradients(fill(0.0, 0, 0), fill(0, 0), 0.0), nnodes(geom));
    # Now loop over all finite element nodes in the map
    lnmap = fill(0, length(fen2fe.map)); # Local node map: buffer to speed up operations
    for nix = 1:length(fen2fe.map)
        gl = fen2fe.map[nix];
        thisnn = nix; # We are at this node
        if !isempty(gl) # This node has an element patch in this block
            # establish local numbering of all nodes of the patch @ node thisnn
            p = patchconn(fes, gl, thisnn);
            np = length(p);
            lnmap[p] .= 1:np;# now store the local numbers
            c = reshape(geom.values[thisnn, :], 1, ndofs(geom))
            updatecsmat!(self.mcsys, c, J, 0);
            gradNavg = fill(0.0, np, ndofs(geom));# preallocate strain-displacement matrix
            Vpatch = 0.0;
            for k = 1:length(gl)
                i = gl[k]
                kconn = collect(fes.conn[i]);
                pci = findfirst(cx -> cx == thisnn, kconn);# at which node in the element are we with this quadrature point?
                @assert 1 <= pci <= nodesperelem(fes)
                # centered coordinates of nodes in the material coordinate system
                for cn = 1:length(kconn)
                    xl[cn, :] = (reshape(geom.values[kconn[cn], :], 1, ndofs(geom)) - c) * self.mcsys.csmat
                end
                jac!(J, xl, lconn, gradNparams[pci])
                At_mul_B!(csmatTJ, self.mcsys.csmat, J); # local Jacobian matrix
                Jac = Jacobianvolume(self.integdomain, J, c, fes.conn[i], Ns[pci]);
                Vpatch += Jac * w[pci];
                sgradN = gradNparams[pci] * adjugate3!(adjJ, J);
                gradNavg[lnmap[kconn],:] += (w[pci] .* sgradN);
            end
            @assert Vpatch != 0
            gradNavg ./= Vpatch;
            bfungrads[nix] = _NodalBasisFunctionGradients(gradNavg, p, Vpatch);
            lnmap[p] .= 0; # Restore the buffer to pristine condition
        end
    end
    self.nodalbasisfunctiongrad = bfungrads
    return self
end

function aspectratio(X)
    edge1 = vec(X[2, :] - X[1, :])
    edge2 = vec(X[3, :] - X[1, :])
    edge3 = vec(X[4, :] - X[1, :])
    edge4 = vec(X[3, :] - X[2, :])
    edge5 = vec(X[4, :] - X[3, :])
    edge6 = vec(X[4, :] - X[2, :])
    V = dot(edge3, cross(edge1, edge2))
    A1 = norm(cross(edge1, edge2)) # This is twice the area of the triangle
    A2 = norm(cross(edge2, edge3))
    A3 = norm(cross(edge3, edge1))
    A4 = norm(cross(edge4, edge6))
    h1, h2, h3, h4 = V/A1, V/A2, V/A3, V/A4
    L1, L2, L3, L4, L5, L6 = norm(edge1), norm(edge2), norm(edge3), norm(edge4), norm(edge5), norm(edge6)
    f = maximum
    return h1/f([L1, L2, L4]), h2/f([L3, L2, L5]), h3/f([L1, L3, L6]), h4/f([L6, L5, L4]), V/6
end

"""
    associategeometry!(self::FEMMAbstractBase,  geom::NodalField{FFlt})

Associate geometry field with the FEMM.

Compute the  correction factors to account for  the shape of the  elements.
"""
function associategeometry!(self::F,  geom::NodalField{FFlt}) where {F<:FEMMDeforLinearESNICET4}
    # The coefficient set below was obtained by fitting the ratio of energies true/approximate
    # for the finite element model of six tetrahedra arranged into a rectangular block
    # and subject to pure bending
    # (a, b) = (1.9726538699841933, 0.23880554398506101)
    (a, b) = (2.101588423297799,  1.311321055432958) # Fitting for a small aspect-ratio range (1.0 to 10)
    fes = self.integdomain.fes
    self.ephis = fill(zero(FFlt), count(fes))
    evols = fill(zero(FFlt), count(fes))
    self.nphis = fill(zero(FFlt), nnodes(geom))
    nvols = fill(zero(FFlt), nnodes(geom))
    for i = 1:count(fes) # Loop over elements
        ar1, ar2, ar3, ar4, V = aspectratio(geom.values[collect(fes.conn[i]), :])
        evols[i] = V;
        ar = sort([ar1, ar2, ar3, ar4])
        self.ephis[i] = (1.0 / (b * minimum(ar) ^a) + 1.0) ^(-1)
        # Accumulate: the stabilization factor at the node is the weighted mean of the stabilization factors of the elements at that node
        for k = 1:nodesperelem(fes)
            nvols[fes.conn[i][k]] += evols[i]
            self.nphis[fes.conn[i][k]] += self.ephis[i] * evols[i]
        end
    end # Loop over elements
    # Now scale the values at the nodes with the nodal volumes
    for k = 1:length(nvols)
        self.nphis[k] /= nvols[k]
    end
    # Now calculate the nodal basis function gradients
    return computenodalbfungrads(self, geom)
end

"""
    stiffness(self::FEMMDeforLinearAbstractNICE, assembler::A,
      geom::NodalField{FFlt},
      u::NodalField{T}) where {A<:SysmatAssemblerBase, T<:Number}

Compute and assemble  stiffness matrix.
"""
function stiffness(self::FEMMDeforLinearAbstractNICE, assembler::A, geom::NodalField{FFlt}, u::NodalField{T}) where {A<:SysmatAssemblerBase, T<:Number}
    fes = self.integdomain.fes
    npts,  Ns,  gradNparams,  w,  pc = integrationdata(self.integdomain);
    dofnums, loc, J, csmatTJ, Jac, D, Dstab = buffers2(self, geom, u, npts)
    self.material.tangentmoduli!(self.material, D, 0.0, 0.0, loc, 0)
    self.stabilization_material.tangentmoduli!(self.stabilization_material, Dstab, 0.0, 0.0, loc, 0)
    elmatsizeguess = 4*nodesperelem(fes)*ndofs(u)
    startassembly!(assembler, elmatsizeguess, elmatsizeguess, nnodes(u) + count(fes), u.nfreedofs, u.nfreedofs);
    for nix = 1:length(self.nodalbasisfunctiongrad)
        gradN = self.nodalbasisfunctiongrad[nix].gradN
        patchconn = self.nodalbasisfunctiongrad[nix].patchconn
        Vpatch = self.nodalbasisfunctiongrad[nix].Vpatch
        c = reshape(geom.values[nix, :], 1, ndofs(geom))
        updatecsmat!(self.mcsys, c, J, 0);
        nd = length(patchconn) * ndofs(u)
        Bnodal = fill(0.0, size(D, 1), nd)
        Blmat!(self.mr, Bnodal, Ns[1], gradN, c, self.mcsys.csmat);
        elmat = fill(0.0, nd, nd) # Can we SPEED it UP?
        DB = fill(0.0, size(D, 1), nd)
        add_btdb_ut_only!(elmat, Bnodal, Vpatch, D, DB)
        add_btdb_ut_only!(elmat, Bnodal, -self.nphis[nix]*Vpatch, Dstab, DB)
        complete_lt!(elmat)
        dofnums = fill(0, nd)
        gatherdofnums!(u, dofnums, patchconn); # retrieve degrees of freedom
        assemble!(assembler, elmat, dofnums, dofnums); # assemble symmetric matrix
    end # Loop over elements
    dofnums, B, DB, elmat, elvec, elvecfix, gradN = buffers3(self, geom, u)
    # OPTIMIZATION: switch to a single-point quadrature rule here
    for i = 1:count(fes) # Loop over elements
        fill!(elmat,  0.0); # Initialize element matrix
        for j = 1:npts # Loop over quadrature points
            locjac!(loc, J, geom.values, fes.conn[i], Ns[j], gradNparams[j])
            Jac = Jacobianvolume(self.integdomain, J, loc, fes.conn[i], Ns[j]);
            updatecsmat!(self.mcsys, loc, J, fes.label[i]);
            At_mul_B!(csmatTJ, self.mcsys.csmat, J); # local Jacobian matrix
            gradN!(fes, gradN, gradNparams[j], csmatTJ);
            Blmat!(self.mr, B, Ns[j], gradN, loc, self.mcsys.csmat);
            add_btdb_ut_only!(elmat, B, self.ephis[i]*Jac*w[j], Dstab, DB)
        end # Loop over quadrature points
        complete_lt!(elmat)
        gatherdofnums!(u, dofnums, fes.conn[i]); # retrieve degrees of freedom
        assemble!(assembler, elmat, dofnums, dofnums); # assemble symmetric matrix
    end # Loop over elements
    return makematrix!(assembler);
end

function stiffness(self::FEMMDeforLinearAbstractNICE, geom::NodalField{FFlt},  u::NodalField{T}) where {T<:Number}
    assembler = SysmatAssemblerSparseSymm();
    return stiffness(self, assembler, geom, u);
end

"""
nzebcloadsstiffness(self::FEMMDeforLinearAbstract,  assembler::A,
  geom::NodalField{FFlt},
  u::NodalField{T}) where {A<:SysvecAssemblerBase, T<:Number}

Compute load vector for nonzero EBC for fixed displacement.
"""
function nzebcloadsstiffness(self::FEMMDeforLinearAbstractNICE,  assembler::A, geom::NodalField{FFlt}, u::NodalField{T}) where {A<:SysvecAssemblerBase, T<:Number}
    error("Not implemented yet")
end

function nzebcloadsstiffness(self::FEMMDeforLinearAbstractNICE, geom::NodalField{FFlt}, u::NodalField{T}) where {T<:Number}
    assembler = SysvecAssembler()
    return  nzebcloadsstiffness(self, assembler, geom, u);
end

end


text:
## Wetted-area method
#======================================================#

# Raymer form factor for wing, pylon, nacelle, fuselage? (Eq. 12.30)
form_factor_wing(x_c, t_c, Λ_m, M) = (1 + 0.6t_c / x_c + 100t_c^4) * (1.34M^0.18 * cos(Λ_m)^0.28)

function form_factor(wing :: HalfWing, M)
	# (x/c)_max, (t/c)_max, (t/c)_max sweep angles
	xcs, tcs, sweeps = max_tbyc_sweeps(wing, 60)
	Kf 				 = form_factor_wing.(xcs, tcs, sweeps, M)
end

# Schlichting averaged skin-friction coefficients
cf_lam(Re_c, k_lam = 1.) = 1.328 / √(Re_c * k_lam)
cf_turb(Re_c, M) = 0.455 / log10(Re_c)^2.58 / (1 + 0.144M^2)^0.65
cf_schlichting(Re, Re_xtr, M) = max(cf_lam(Re), cf_turb(Re, M) - (Re_xtr / 320 - 39) / Re)

function wetted_area_drag(mean_chords, S_wets, K_fs, x_tr, V, ρ, M, μ)
	# Skin-friction coefficients
	Re_c 		= reynolds_number.(ρ, V, mean_chords, μ)
	Re_xtr 		= reynolds_number.(ρ, V, mean_chords .* x_tr, μ)
	cfs 		= cf_schlichting.(Re_c, Re_xtr, M)

	# Profile drag
	Dp_by_q	 	= sum(@. cfs * S_wets * K_fs)
end

function wetted_area_drag(wing :: HalfWing, x_tr, V, ρ, a_ref = 330., μ = 1.5e-5)
	# Chord processing
	mean_chords = (fwdsum ∘ chords)(wing) / 2

	# Wetted areas
	S_wets = @. mean_chords * wing.spans / cos(wing.dihedrals)

	# Form factors
	M 	 = V / a_ref
	K_fs = form_factor(wing, M)

	wetted_area_drag(mean_chords, S_wets, K_fs, x_tr, V, ρ, M, μ)
end

profile_drag_coefficient(wing :: HalfWing, x_tr, V, rho_ref, a_ref, area_ref, μ) = wetted_area_drag(wing, x_tr, V, rho_ref, a_ref, μ) / area_ref
profile_drag_coefficient(wing :: Wing, x_tr, V, rho_ref, a_ref, area_ref, μ) = profile_drag_coefficient(left(wing), x_tr, V, rho_ref, a_ref, area_ref, μ) + profile_drag_coefficient(right(wing), x_tr, V, rho_ref, a_ref, area_ref, μ)

## Local-friction and local-dissipation method

max_forks_repo_name: GodotMisogi/AeroMDAO


text:
function test_binary_accuracy_layer(backend::Backend, tensor_dim, T, threshold, eps)
  println("-- Testing BinaryAccuracyLayer on $(typeof(backend)){$T}...")

  ############################################################
  # Prepare Data for Testing
  ############################################################
  dims = tuple(rand(6:11, tensor_dim)...)
  println("    > $dims")
  if 0 == threshold
    preds = rand(T, dims)*4-2
    labels = round.(rand(T, dims))*2-1
  elseif 0.5 == threshold
    preds = round.(rand(T, dims))
    labels = rand(T, dims)
  else
    error("Threshold must be 0 or 0.5; was $threshold")
  end
  errs_mask = (preds.>threshold) .≠ (labels.>threshold)
  n = prod(dims)
  n_wrong = countnz(errs_mask)
  n_right = n - n_wrong

  ############################################################
  # Setup
  ############################################################
  layer = BinaryAccuracyLayer(; bottoms=[:predictions, :labels], threshold=threshold)
  pred_blob  = make_blob(backend, T, dims)
  label_blob = make_blob(backend, T, dims)
  copy!(pred_blob, preds)
  copy!(label_blob, labels)
  inputs = Blob[pred_blob, label_blob]

  state  = setup(backend, layer, inputs, Blob[])

  ###########################################################
  # Forward
  ###########################################################
  println("    > Forward")
  forward(backend, state, inputs)

  @test state.n_accum == n

  expected_acc = convert(T, n_right / n)
  @test abs(state.accuracy - expected_acc) < eps

  println("    > Forward Again")
  forward(backend, state, Blob[label_blob, label_blob])
  @test state.n_accum == 2n
  @test abs(state.accuracy - (expected_acc+1)/2) < eps

  println("    > Forward Again and Again")
  reset_statistics(state)
  forward(backend, state, inputs)
  @test state.n_accum == n
  @test abs(state.accuracy - expected_acc) < eps

  shutdown(backend, state)
end
function test_binary_accuracy_layer(backend::Backend, T, eps)
  for i in [2,4,5]
    for threshold in [0, 0.5]
      test_binary_accuracy_layer(backend, i, T, threshold, eps)
    end
  end
end
function test_binary_accuracy_layer(backend::Backend)
  test_binary_accuracy_layer(backend, Float32, 1e-4)
  test_binary_accuracy_layer(backend, Float64, 1e-8)
end

if test_cpu
  test_binary_accuracy_layer(backend_cpu)
end
if test_gpu
  test_binary_accuracy_layer(backend_gpu)
end


text:
"""
--------------------------------------------------------------------------------------\n
Author: Huckleberry Febbo, Graduate Student, University of Michigan
Date Create: 04/26/2018 Last Modified: 04/26/2018 \n
--------------------------------------------------------------------------------------\n
"""
function mpcPlot(n,idx,kwargs...)
  if !isdir(n.r.resultsDir); resultsDir!(n); end
  stp = [statePlot(n,idx,st;kwargs...) for st in 1:n.ocp.state.num]
  ctp = [controlPlot(n,idx,ctr;kwargs...) for ctr in 1:n.ocp.control.num]

  tp = tPlot(n,idx)

  if n.s.ocp.evalCostates && n.s.ocp.integrationMethod == :ps && n.s.ocp.evalConstraints
    csp = [costatePlot(n,idx,st;kwargs...) for st in 1:n.ocp.state.num]
    all = [stp;ctp;tp;csp]
  else
    all = [stp;ctp;tp]
  end

  h = plot(all...,size=_pretty_defaults[:size])
  if !_pretty_defaults[:simulate]; savefig(string(n.r.resultsDir,"mpc.",_pretty_defaults[:format])) end
  return h
end
"""
tp=tPlot(n,r,idx)
tp=tPlot(n,r,idx,tp;(:append=>true))
# plot the optimization times
--------------------------------------------------------------------------------------\n
Author: Huckleberry Febbo, Graduate Student, University of Michigan
Date Create: 3/11/2017, Last Modified: 12/02/2019 \n
Last Modified By: Alexander Buck
--------------------------------------------------------------------------------------\n
"""
function tPlot(n,idx::Int64,args...;kwargs...);
  r = n.r.ocp

  kw = Dict(kwargs);
  # check to see if user would like to add to an existing plot
  if !haskey(kw,:append); append=false;
  else; append = get(kw,:append,0);
  end
  if !append; tp=plot(0,leg=:false); else tp=args[1]; end

  if idx > length(r.dfsOpt[:tSolve])
      @warn "Cannot plot idx = ", idx, " because length(r.dfsOpt[:tSolve]) = ", length(r.dfsOpt[:tSolve]), ". \n
            reducing idx in tPlot!()."
      idx = length(r.dfsOpt[:tSolve])
  end

  # check to see if user would like to label legend
  if !haskey(kw,:legend);legend_string="";
  else; legend_string = get(kw,:legend,0);
  end

  # to avoid a bunch of jumping around in the simulation
	idx_max = length(r.dfsOpt[:tSolve])
	if (idx_max<10); idx_max=10 end

  if idx > 1
    scatter!(1:idx-1,r.dfsOpt[:tSolve][1:idx-1],marker=_pretty_defaults[:opt_marker][3],label=string(legend_string,"Previous Times"))
  end
  if isequal(r.dfsOpt[:status][idx],:Optimal)
    scatter!((idx,r.dfsOpt[:tSolve][idx]),marker=_pretty_defaults[:opt_marker][1],label=string(legend_string,"Optimal"))
  else
    scatter!((idx,r.dfsOpt[:tSolve][idx]),marker=_pretty_defaults[:opt_marker][2],label=string(legend_string,string(r.dfsOpt[:status][idx])) )
  end

  plot!(1:length(r.dfsOpt[:tSolve]),n.mpc.v.tex*ones(length(r.dfsOpt[:tSolve])),line=_pretty_defaults[:limit_lines][2],label="real-time threshhold",leg=:topright)

	ylims!(0,max(n.mpc.v.tex*1.2, maximum(r.dfsOpt[:tSolve])))
  xlims!(1,length(r.dfsOpt[:tSolve]))
	yaxis!("Optimization Time (s)")
	xaxis!("Evaluation Number")
  plot!(size=_pretty_defaults[:size]);
	if !_pretty_defaults[:simulate] savefig(string(n.r.resultsDir,"tplot.",_pretty_defaults[:format])) end

	return tp
end

"""
--------------------------------------------------------\n
Author: Huckleberry Febbo, Graduate Student, University of Michigan
Date Create: 7/04/2017, Last Modified: 7/04/2017 \n
--------------------------------------------------------------------------------------\n
"""
function optPlot(n)
  L = length(n.r.ocp.dfsOpt[:tSolve])
  temp = [n.r.ocp.dfsOpt[:objVal][jj] for jj in 1:L]
  val = [idx for tempM in temp for idx=tempM]
  opt = plot(1:L,val)
  yaxis!("Objective Function Values")
	xaxis!("Evaluation Number")
  savefig(string(n.r.resultsDir,"optPlot.",_pretty_defaults[:format]))
  return opt
end


text:
# Color grid so that all elements containing the i-th node have a different color.
# During threaded assembly each thread acts upon a single color. In this way, we
# can multithread along each color avoiding synchronisations between threads. It
# could be extended amongst different cpus, where each cpu takes different colors
function color_mesh(e2n)
    ## ---- Make list of elements within i-th node
    els = view(e2n,1:3,:)
    nel = size(els,2)
    els_in_node = elements_in_node(els)

    ## ---- Incidence matrix
    I, J, V = Int[], Int[], Bool[]
    for cells in els_in_node
        for cell1 in cells # All these cells have a neighboring node
            for cell2 in cells
                if cell1 != cell2
                    push!(I, cell1)
                    push!(J, cell2)
                    push!(V, true)
                end
            end
        end
    end
    incidence_matrix = sparse(I, J, V)

    ##  Fill elemental color matrix
    els_colors = fill(0,nel)
    max_color = 0
    @inbounds for iel in 1:nel
        occupied_color = Int64[]

        for r in nzrange(incidence_matrix, iel)
            el_neighbour = incidence_matrix.rowval[r]
            color = els_colors[el_neighbour]
            if color !=0
                push!(occupied_color, color)
            end
        end

        free_color = 0
        for color in 1:max_color
            if color ∉ occupied_color
                free_color = color
                break
            end
        end

        if free_color == 0 # no free color found, need to bump max colors
            max_color += 1
            free_color = max_color
        end

        els_colors[iel] = free_color
    end

    ## List of elements per color (this is what is needed for threading the assembler)
    # max_color -= 1
    color_list = [Int32[] for _ in 1:max_color]
    for (iel, color) in enumerate(els_colors)
        @inbounds push!(color_list[color],iel)
    end

    return els_colors, sort!(color_list)
end

## Make list of elements sharing i-th node
function elements_in_node(els)
    nodes = 1:maximum(els)
    els_in_node = [Int64[] for _ in axes(nodes,1)]
    for id in axes(els,2), v in view(els,:,id)
        @inbounds push!(els_in_node[v],id)
    end
    els_in_node
end

# Make vtk file of colored grid (WriteVTK pagkage)
function colored_grid_vtk(colors, gr)
    cells = [MeshCell(VTKCellTypes.VTK_TRIANGLE, gr.e2n[1:3,i]) for i in axes(gr.e2n,2)]
    vtkfile = vtk_grid("colored_grid", gr.x, gr.z, cells)
    vtkfile["color", VTKCellData()] = colors
    vtk_save(vtkfile)
end

# Make vtk file of colored grid
function colored_grid_vtk(colors, x, z, e2n)
    cells = [MeshCell(VTKCellTypes.VTK_TRIANGLE, e2n[1:3,i]) for i in axes(e2n,2)]
    vtkfile = vtk_grid("colored_grid", x, z, cells)
    vtkfile["color", VTKCellData()] = colors
    vtk_save(vtkfile)
end

#=
    check_colored_mesh(color_list)

    Checks whether coloring was succesfull
=#
function check_colored_mesh(color_list)
    ncolors = length(color_list)

    for (nc, node_list) in enumerate(color_list)
        for color in 1:ncolors
            if color!=nc
                for node in node_list
                    if node ∈ color_list[color]
                        println("node ", node, " found in color ", color)
                    end
                end
            end
        end
    end

end

max_stars_repo_name: albert-de-montserrat/Persephone

max_issues_repo_name: albert-de-montserrat/Persephone

max_forks_repo_name: albert-de-montserrat/Persephone


text:
#
# blas_get_num_threads()
# Get the number of BLAS threads
# This can be replaced by BLAS.get_num_threads() in Julia v1.6
#

function guess_vendor()
  # like determine_vendor, but guesses blas in some cases
  # where determine_vendor returns :unknown
  ret = BLAS.vendor()
  if Sys.isapple() && (ret == :unknown)
    ret = :osxblas
  end
  return ret
end

_tryparse_env_int(key) = tryparse(Int, get(ENV, key, ""))

blas_get_num_threads()::Union{Int,Nothing} = _get_num_threads()

function _get_num_threads(; _blas=guess_vendor())::Union{Int,Nothing}
  if _blas === :openblas || _blas === :openblas64
    return Int(ccall((BLAS.@blasfunc(openblas_get_num_threads), BLAS.libblas), Cint, ()))
  elseif _blas === :mkl
    return Int(ccall((:mkl_get_max_threads, BLAS.libblas), Cint, ()))
  elseif _blas === :osxblas
    key = "VECLIB_MAXIMUM_THREADS"
    nt = _tryparse_env_int(key)
    if nt === nothing
      @warn "Failed to read environment variable $key" maxlog = 1
    else
      return nt
    end
  else
    @assert _blas === :unknown
  end
  @warn "Could not get number of BLAS threads. Returning `nothing` instead." maxlog = 1
  return nothing
end


text:
t = Template(; user=me)
pkg_dir = joinpath(t.dir, test_pkg)

@testset "CirrusCI" begin
    @testset "Plugin creation" begin
        p = CirrusCI()
        @test isempty(p.gitignore)
        @test p.src == joinpath(DEFAULTS_DIR, "cirrus.yml")
        @test p.dest == ".cirrus.yml"
        @test p.badges == [
            Badge(
                "Build Status",
                "https://api.cirrus-ci.com/github/{{USER}}/{{PKGNAME}}.jl.svg",
                "https://cirrus-ci.com/github/{{USER}}/{{PKGNAME}}.jl",
            ),
        ]
        @test isempty(p.view)
        p = CirrusCI(; config_file=nothing)
        @test p.src === nothing
        p = CirrusCI(; config_file=test_file)
        @test p.src == test_file
        @test_throws ArgumentError CirrusCI(; config_file=fake_path)
    end

    @testset "Badge generation" begin
        p = CirrusCI()
        @test badges(p, me, test_pkg) == ["[![Build Status](https://api.cirrus-ci.com/github/$me/$test_pkg.jl.svg)](https://cirrus-ci.com/github/$me/$test_pkg.jl)"]
    end

    @testset "File generation" begin
        # Without a coverage plugin in the template, there should be no coverage step.
        p = CirrusCI()
        @test gen_plugin(p, t, test_pkg) == [".cirrus.yml"]
        @test isfile(joinpath(pkg_dir, ".cirrus.yml"))
        cirrus = read(joinpath(pkg_dir, ".cirrus.yml"), String)
        @test !occursin("coverage_script", cirrus)
        rm(joinpath(pkg_dir, ".cirrus.yml"))

        # Generating the plugin with Codecov in the template should create a post-test step.
        t.plugins[Codecov] = Codecov()
        gen_plugin(p, t, test_pkg)
        delete!(t.plugins, Codecov)
        cirrus = read(joinpath(pkg_dir, ".cirrus.yml"), String)
        @test occursin("coverage_script", cirrus)
        @test occursin("cirrusjl coverage", cirrus)
        rm(joinpath(pkg_dir, ".cirrus.yml"))

        p = CirrusCI(; config_file=nothing)
        @test isempty(gen_plugin(p, t, test_pkg))
        @test !isfile(joinpath(pkg_dir, ".cirrus.yml"))
    end
end

rm(pkg_dir; recursive=true)


text:
module IsingMCMC
##export E, MCMC

function genFunctor( getJ, getB)
    getJ = getJ
    getB = getB

    # E(S) = - 1/2\sum_{i,j, i not j} J_{i,j} si sj - \sum_i B_i si
    function E(S, N)
        Et = 0
        for i in 1:N
            J = Main.getNeighbors(i)
            Et += -0.5*S[i]*sum( [getJ(i,j) * S[j] for j in J] )
            Et += -    S[i]*getB(i)
        end
        Et
    end

    # differential energy
    function dE(S, k, N)
        dEt = 0
        J = Main.getNeighbors(k)
        dEt +=   S[k]*sum( [(getJ(k,j) + getJ(j,k)) * S[j] for j in J] )
        dEt += 2*S[k]*getB(k)
        dEt
    end

    # probability for dE
    PrE(dE, T) = exp(-dE/T)

    # flip spin at x
    function flipx!(list, x)
        list[x] *= -1
    end

    # MCMC
    function MCMC( T, N, trial)
        # initialize
        sim = []
        S = ones(N)
        for i in 1:N
            if rand() < 0.5
                flipx!(S,i) # random flip at first
            end
        end
        # Gibbs sampling position
        gpos = [rand(1:N) for i in 1:trial]
        # random value in [0,1]
        rn = [rand() for i in 1:trial]

        # MCMC trial
        for t in 1:trial
            k = gpos[t]
            de = dE(S, k, N)
            if rn[t] < PrE(de, T) # MH criteria
                flipx!(S, k) # change k
            end
            push!(sim, copy(S))
        end

        sim
    end

    return (E, MCMC)
end # of functor

end # of module


text:
"""
    population_data()

Returns a dictionary `D` with the US population data used in section 9.2.
The items in the dictionary are three vectors of length 100.

* `D["population"]`:  2010 population in millions for ages 0,...,99
* `D["birth_rate"]`:  birth rate
* `D["death_rate"]`:  death rate.
"""
function population_data()
population = [
   3.94415;
   3.97807;
   4.09693;
   4.11904;
   4.06317;
   4.05686;
   4.06638;
   4.03058;
   4.04649;
   4.14835;
   4.17254;
   4.11442;
   4.10624;
   4.11801;
   4.16598;
   4.24282;
   4.31614;
   4.39529;
   4.50085;
   4.58523;
   4.51913;
   4.35429;
   4.26464;
   4.19857;
   4.24936;
   4.26235;
   4.15231;
   4.24887;
   4.21525;
   4.22308;
   4.28567;
   3.97022;
   3.98685;
   3.88015;
   3.83922;
   3.95643;
   3.80209;
   3.93445;
   4.12188;
   4.36480;
   4.38327;
   4.11498;
   4.07610;
   4.10511;
   4.21150;
   4.50887;
   4.51976;
   4.53526;
   4.53880;
   4.60590;
   4.66029;
   4.46463;
   4.50085;
   4.38035;
   4.29200;
   4.25471;
   4.03751;
   3.93639;
   3.79493;
   3.64127;
   3.62113;
   3.49260;
   3.56318;
   3.48388;
   2.65713;
   2.68076;
   2.63914;
   2.64936;
   2.32367;
   2.14232;
   2.04312;
   1.94932;
   1.86427;
   1.73696;
   1.68449;
   1.62008;
   1.47107;
   1.45533;
   1.40012;
   1.37119;
   1.30851;
   1.21287;
   1.16142;
   1.07481;
   0.98572;
   0.91472;
   0.81421;
   0.71291;
   0.64062;
   0.53800;
   0.43556;
   0.34499;
   0.28139;
   0.21698;
   0.16944;
   0.12972;
   0.09522;
   0.06814;
   0.04590;
   0.03227;
];
d = [
   0.00623;
   0.00044;
   0.00027;
   0.00020;
   0.00016;
   0.00012;
   0.00011;
   0.00011;
   0.00012;
   0.00011;
   0.00010;
   0.00013;
   0.00013;
   0.00015;
   0.00020;
   0.00025;
   0.00037;
   0.00047;
   0.00064;
   0.00071;
   0.00076;
   0.00087;
   0.00087;
   0.00088;
   0.00094;
   0.00092;
   0.00095;
   0.00093;
   0.00099;
   0.00101;
   0.00103;
   0.00109;
   0.00110;
   0.00114;
   0.00115;
   0.00120;
   0.00131;
   0.00137;
   0.00146;
   0.00156;
   0.00162;
   0.00185;
   0.00201;
   0.00216;
   0.00243;
   0.00258;
   0.00298;
   0.00325;
   0.00351;
   0.00387;
   0.00413;
   0.00454;
   0.00494;
   0.00533;
   0.00571;
   0.00602;
   0.00670;
   0.00710;
   0.00769;
   0.00828;
   0.00860;
   0.00932;
   0.00998;
   0.01101;
   0.01250;
   0.01282;
   0.01404;
   0.01515;
   0.01687;
   0.01830;
   0.01967;
   0.02133;
   0.02347;
   0.02562;
   0.02800;
   0.03083;
   0.03441;
   0.03711;
   0.04126;
   0.04448;
   0.04964;
   0.05539;
   0.06149;
   0.06803;
   0.07673;
   0.08561;
   0.09540;
   0.10636;
   0.11802;
   0.13385;
   0.15250;
   0.16491;
   0.18738;
   0.20757;
   0.22688;
   0.25196;
   0.27422;
   0.29239;
   0.32560;
   0.34157;
];
b = [
   0.00000;
   0.00000;
   0.00000;
   0.00000;
   0.00000;
   0.00000;
   0.00000;
   0.00000;
   0.00000;
   0.00000;
   0.00020;
   0.00020;
   0.00020;
   0.00020;
   0.00020;
   0.01710;
   0.01710;
   0.01710;
   0.01710;
   0.01710;
   0.04500;
   0.04500;
   0.04500;
   0.04500;
   0.04500;
   0.05415;
   0.05415;
   0.05415;
   0.05415;
   0.05415;
   0.04825;
   0.04825;
   0.04825;
   0.04825;
   0.04825;
   0.02250;
   0.02250;
   0.02250;
   0.02250;
   0.02250;
   0.00510;
   0.00510;
   0.00510;
   0.00510;
   0.00510;
   0.00035;
   0.00035;
   0.00035;
   0.00035;
   0.00035;
   0.00000;
   0.00000;
   0.00000;
   0.00000;
   0.00000;
   0.00000;
   0.00000;
   0.00000;
   0.00000;
   0.00000;
   0.00000;
   0.00000;
   0.00000;
   0.00000;
   0.00000;
   0.00000;
   0.00000;
   0.00000;
   0.00000;
   0.00000;
   0.00000;
   0.00000;
   0.00000;
   0.00000;
   0.00000;
   0.00000;
   0.00000;
   0.00000;
   0.00000;
   0.00000;
   0.00000;
   0.00000;
   0.00000;
   0.00000;
   0.00000;
   0.00000;
   0.00000;
   0.00000;
   0.00000;
   0.00000;
   0.00000;
   0.00000;
   0.00000;
   0.00000;
   0.00000;
   0.00000;
   0.00000;
   0.00000;
   0.00000;
   0.00000;
];

    return Dict([("population", population),  ("birth_rate", b),
        ("death_rate", d)])
end

max_issues_repo_name: lynnlangit/julia-linear-algebra


text:
#include("../../src/ttv.jl")
#include("/Users/ericagol/Computer/Julia/regress.jl")

using NbodyGradient, PyPlot, DelimitedFiles
include("linreg.jl")

# This routine takes derivative of transit times with respect
# to the initial orbital elements.
n = 3
# Try a larger value of t0:
#t0 =  -200.0
t0 =  0.0
h  = 0.06
tmax = 400.0

# Read in initial conditions:
elements = readdlm("elements.txt",',')[1:n,:]

ntt = zeros(Int64,n)

# Make an array, tt,  to hold transit times:
# First, though, make sure it is large enough:
for i=2:n
  ntt[i] = ceil(Int64,tmax/elements[i,2])+3
end

nsub = 7 # How many times to divide the timestep (/2^(nsub-1)).
tt_save = zeros(nsub,n,maximum(ntt))
# Call the ttv function:

ic = ElementsIC(t0, n, elements[1:n,:])
intr = Integrator(h, t0, tmax)
s = State(ic)
tt = TransitTiming(intr.tmax, ic)
@time intr(s,tt)
tt_save[1,:,:] .= tt.tt

# Now, compute derivatives (with respect to initial cartesian positions/masses):
intr = Integrator(h/2, t0, tmax)
s = State(ic)
tt = TransitTiming(intr.tmax, ic)
@time intr(s,tt)
tt_save[2,:,:] .= tt.tt

intr = Integrator(h/4, t0, tmax)
s = State(ic)
tt = TransitTiming(intr.tmax, ic)
@time intr(s,tt)
tt_save[3,:,:] .= tt.tt

intr = Integrator(h/8, t0, tmax)
s = State(ic)
tt = TransitTiming(intr.tmax, ic)
@time intr(s,tt)
tt_save[4,:,:] .= tt.tt

intr = Integrator(h/16, t0, tmax)
s = State(ic)
tt = TransitTiming(intr.tmax, ic)
@time intr(s,tt)
tt_save[5,:,:] .= tt.tt

intr = Integrator(h/32, t0, tmax)
s = State(ic)
tt = TransitTiming(intr.tmax, ic)
@time intr(s,tt)
tt_save[6,:,:] .= tt.tt

intr = Integrator(h/128, t0, tmax)
s = State(ic)
tt = TransitTiming(intr.tmax, ic)
@time intr(s,tt)
tt_save[7,:,:] .= tt.tt

# Make a plot of transit time errors versus stepsize:
clf()
sigt = zeros(n-1,nsub-1)
tab = 0
#h_list = [h,h/2.,h/4.,h/8.]
h_list = [h,h/2.,h/4.,h/8,h/16,h/32]
#hlabel = ["h-h/64","h/2-h/64","h/4-h/64","h/8-h/64","h/16-h/64","h/32-h/64"]
hlabel = ["h-h/128","h/2-h/128","h/4-h/128","h/8-h/128","h/16-h/128","h/32-h/128"]
ch = ["black","red","green","blue","orange","magenta"]
for i=2:n
  for j=1:nsub-1
    tti1 = tt_save[j,i,1:tt.count[i]]
    tti_max = tt_save[nsub,i,1:tt.count[i]]
    diff = (tti1 .- tti_max)*24*3600 # seconds
    sigt[i-1,j]=std(diff)
    if i == n
#      plot(tti1,diff ./ h_list[j]^4,linestyle="dashed",c=ch[j])
#      plot(tti1,diff ,linestyle="dashed",c=ch[j])
      semilogy(tti1,abs.(diff) ,linestyle="dashed",c=ch[j])
    else
#      plot(tti1,diff ./ h_list[j]^4,label=hlabel[j],c=ch[j])
#      plot(tti1,diff,label=hlabel[j],c=ch[j])
      plot(tti1,abs.(diff),label=hlabel[j],c=ch[j])
    end
  end
end
xlabel("Elapsed time")
ylabel("Difference in transit time [sec]")
plot(tti1,diff .* h^4,".")
plot([50,50,45,50,55],[1e-2,1e-2/2^4,1e-2/2^4*1.5,1e-2/2^4,1e-2/2^4*1.5])
text(55,3e-3,L"$2^4$")
plot(tti1,tti1 .* 1e-3, linestyle="dotted")
text(100,0.2,L"$10^{-3}t$")
axis([-10,410,0.5e-10,1.0])
legend(ncol=2)

#read(stdin,Char)
#
## Make a plot of timing errors versus stepsize:
#clf()
#sigt = zeros(n-1,nsub-1)
#tab = 0
#for i=2:n
#  fn = zeros(Float64,2,count[i])
#  sig = ones(count[i])
#  tti_max = tt_save[nsub,i,1:count[i]]
#  for j=1:count[i]
#    fn[1,j] = 1.0
#    fn[2,j] = round(Int64,(tti_max[j]-elements[i,3])/elements[i,2])
#  end
#  for j=1:nsub-1
#    tti1 = tt_save[j,i,1:count[i]]
#    #coeff,cov = regress(fn,tti1-tti_max,sig)
#    p = Regression(fn[2,:],tti1 .- tti_max)
#    coeff = [p.α, p.β]
#    diff = tti1 .- tti_max .- coeff[1] .- (coeff[2] * fn[2,:])
#    sigt[i-1,j]=std(diff)
##    print("planet: ",i-1," stepsize: ",j," sig: ",sigt[i-1,j])
#    if i == n
#      plot(tti1,diff ./ h_list[j]^4,linestyle="dashed",c=ch[j])
#    else
#      plot(tti1,diff ./ h_list[j]^4,label=hlabel[j],c=ch[j])
#    end
#  end
#end
#ylabel("TTV difference / h^4")
#xlabel("Transit time")
#legend(loc="lower left")
#
# #read(stdin,Char)
# clf()
#
#
# # Make a plot of some TTVs:
# loglog(h_list,sigt[1,:] .* 24.0 .* 3600.0,".",markersize=15,label="planet b")
# loglog(h_list,sigt[1,1] .* 24.0 .* 3600.0 .* (h_list ./ h[1]).^4,label=L"$\propto h^4$")
# loglog(h_list,sigt[2,:] .* 24.0 .* 3600.0,".",markersize=15,label="planet c")
# loglog(h_list,sigt[2,1] .* 24.0 .* 3600.0 .* (h_list ./ h[1]).^4,label=L"$\propto h^4$")
# legend(loc = "upper left")
# ylabel("RMS TTV error [sec]")
# xlabel("Step size [day]")
#
## PyPlot.savefig("timing_error_vs_h_kep88_no_outer.pdf",bbox_inches="tight")

max_forks_repo_name: langfzac/nbg-papers


text:
include("common.jl")

mmesh = Mesh("rectangle.stl")
c1=(1.0, 0.0)
c2=(1.0, 1.0)
fault_idx, fault_coords, fault_normal_vec = get_fault_info(mmesh, c1, c2)

mmesh, conn = fix_mesh(mmesh, fault_idx, rectangle_redblack)
close("all")
visualize_mesh(mmesh)
v = 0.05*(rand(mmesh.nnode).-0.5)
u = 0.05*(rand(mmesh.nnode).-0.5)
scatter(mmesh.nodes[:,1] + u, mmesh.nodes[:,2] + v, s = 2)
# plot(fault_coords[1,:], fault_coords[2,:], "r.", label="fault")
# legend()
savefig("mesh.png")

bdnode = bcnode((x, y)->(x<1e-5 || x>2-1e-5), mmesh)

Conn = compute_slip_boundary_condition( conn, mmesh)
NT = 100
Δt = 0.004
λ = 1.0
μ = 1.0
invη = ones(get_ngauss(mmesh))
sn = size(conn,1)
slipvec = [zeros(NT, sn) ones(NT, sn)] # acceleration
gravity = 0.0

U, Sigma, Varepsilon = subduction_solver(invη, gravity, Conn, slipvec, bdnode, mmesh)

sess = Session(); init(sess)
Uval, Stress, Strain = run(sess, [U, Sigma, Varepsilon])
p = visualize_displacement(Uval[1:5:end,:], mmesh)
saveanim(p, "subduction.gif")

close("all")
u = Uval[end,:]
# u = (@. sqrt(u[1:mmesh.ndof]^2 + u[mmesh.ndof+1:2mmesh.ndof]^2))
visualize_scalar_on_fem_points(u[mmesh.ndof+1:2mmesh.ndof],mmesh)
savefig("stress_yy.png")
# node = [0.0 0.0
#         1.0 0.0
#         0.0 1.0
#         5.0 0.0
#         6.0 0.0
#         5.0 1.0]

# elem = [1 2 3
#         4 5 6]

# mmesh = Mesh(node, elem)
# close("all")
# visualize_mesh(mmesh)
# savefig("mesh.png")


text:
import Base: isfinite, isnan, precision, iszero, eps,
    typemin, typemax, floatmin, floatmax,
    sign_mask, exponent_mask, exponent_one, exponent_half,
    significand_mask, round, Int32, Int64,
    +, -, *, /, ^, ==, <, <=, >=, >, !=, inv,
    abs, sqrt, exp, log, log2, log10, sin, cos, tan, asin,
    acos, atan, sinh, cosh, tanh, asinh, acosh, atan,
    exponent,
    bitstring

primitive type BFloat16 <: AbstractFloat 16 end

# Floating point property queries
for f in (:sign_mask, :exponent_mask, :exponent_one,
            :exponent_half, :significand_mask)
    @eval $(f)(::Type{BFloat16}) = UInt16($(f)(Float32) >> 16)
end

iszero(x::BFloat16) = reinterpret(UInt16, x) & ~sign_mask(BFloat16) == 0x0000
isfinite(x::BFloat16) = (reinterpret(UInt16,x) & exponent_mask(BFloat16)) != exponent_mask(BFloat16)
isnan(x::BFloat16) = (reinterpret(UInt16,x) & ~sign_mask(BFloat16)) > exponent_mask(BFloat16)
precision(::Type{BFloat16}) = 8
eps(::Type{BFloat16}) = Base.bitcast(BFloat16, 0x3c00)

round(x::BFloat16, r::RoundingMode{:Up}) = BFloat16(ceil(Float32(x)))
round(x::BFloat16, r::RoundingMode{:Down}) = BFloat16(floor(Float32(x)))
round(x::BFloat16, r::RoundingMode{:Nearest}) = BFloat16(round(Float32(x)))

Int64(x::BFloat16) = Int64(Float32(x))
Int32(x::BFloat16) = Int32(Float32(x))

## floating point traits ##
"""
    InfB16
Positive infinity of type [`BFloat16`](@ref).
"""
const InfB16 = reinterpret(BFloat16, 0x7f80)

"""
    NaNB16
A not-a-number value of type [`BFloat16`](@ref).
"""
const NaNB16 = reinterpret(BFloat16, 0x7fc0)

# More floating point property queries
typemin(::Type{BFloat16}) = -InfB16
typemax(::Type{BFloat16}) = InfB16
floatmax(::Type{BFloat16}) = reinterpret(BFloat16, 0x7f7f)
floatmin(::Type{BFloat16}) = reinterpret(BFloat16, 0x0080)

# Truncation from Float32
Base.uinttype(::Type{BFloat16}) = UInt16
Base.trunc(::Type{BFloat16}, x::Float32) = reinterpret(BFloat16,
        (reinterpret(UInt32, x) >> 16) % UInt16
    )

# Conversion from Float32
function BFloat16(x::Float32)
    isnan(x) && return NaNB16
    # Round to nearest even (matches TensorFlow and our convention for
    # rounding to lower precision floating point types).
    h = reinterpret(UInt32, x)
    h += 0x7fff + ((h >> 16) & 1)
    return reinterpret(BFloat16, (h >> 16) % UInt16)
end

# Conversion from Float64
function BFloat16(x::Float64)
	BFloat16(Float32(x))
end

# Conversion from Integer
function BFloat16(x::Integer)
	convert(BFloat16, convert(Float32, x))
end

# Expansion to Float32
function Base.Float32(x::BFloat16)
    reinterpret(Float32, UInt32(reinterpret(UInt16, x)) << 16)
end

# Expansion to Float64
function Base.Float64(x::BFloat16)
    Float64(Float32(x))
end

# Truncation to integer types
Base.unsafe_trunc(T::Type{<:Integer}, x::BFloat16) = unsafe_trunc(T, Float32(x))

# Basic arithmetic
for f in (:+, :-, :*, :/, :^)
    @eval ($f)(x::BFloat16, y::BFloat16) = BFloat16($(f)(Float32(x), Float32(y)))
end
-(x::BFloat16) = reinterpret(BFloat16, reinterpret(UInt16, x) ⊻ sign_mask(BFloat16))
^(x::BFloat16, y::Integer) = BFloat16(^(Float32(x), y))

for F in (:abs, :sqrt, :exp, :log, :log2, :log10,
          :sin, :cos, :tan, :asin, :acos, :atan,
          :sinh, :cosh, :tanh, :asinh, :acosh, :atanh)
  @eval begin
    $F(x::BFloat16) = BFloat16($F(Float32(x)))
  end
end

const ZeroBFloat16 = BFloat16(0.0f0)
const OneBFloat16 = BFloat16(1.0f0)
Base.zero(::Type{BFloat16}) = ZeroBFloat16
Base.one(::Type{BFloat16}) = OneBFloat16

inv(x::BFloat16) = one(BFloat16) / x

# Floating point comparison
function ==(x::BFloat16, y::BFloat16)
    ix = reinterpret(UInt16, x)
    iy = reinterpret(UInt16, y)
    # NaNs (isnan(x) || isnan(y))
    if (ix|iy)&~sign_mask(BFloat16) > exponent_mask(BFloat16)
        return false
    end
    # Signed zeros
    if (ix|iy)&~sign_mask(BFloat16) == 0
        return true
    end
    return ix == iy
end

for op in (:<, :<=, :>, :>=, :!=)
    @eval ($op)(a::BFloat16, b::BFloat16) = ($op)(Float32(a), Float32(b))
end

Base.widen(::Type{BFloat16}) = Float32
Base.promote_rule(::Type{Float32}, ::Type{BFloat16}) = Float32
Base.promote_rule(::Type{Float64}, ::Type{BFloat16}) = Float64
for t in (Int8, Int16, Int32, Int64, Int128, UInt8, UInt16, UInt32, UInt64, UInt128)
    @eval Base.promote_rule(::Type{BFloat16}, ::Type{$t}) = BFloat16
end

# Wide multiplication
Base.widemul(x::BFloat16, y::BFloat16) = Float32(x) * Float32(y)

# Showing
function Base.show(io::IO, x::BFloat16)
    hastypeinfo = BFloat16 === get(io, :typeinfo, Any)
    if isinf(x)
        print(io, x < 0 ? "-InfB16" : "InfB16")
    elseif isnan(x)
        print(io, "NaNB16")
    else
        hastypeinfo || print(io, "BFloat16(")
        show(IOContext(io, :typeinfo=>Float32), Float32(x))
        hastypeinfo || print(io, ")")
    end
end

# Random
import Random: rand, randn, randexp, AbstractRNG, Sampler
rand(rng::AbstractRNG, ::Sampler{BFloat16}) = convert(BFloat16, rand(rng))
randn(rng::AbstractRNG, ::Type{BFloat16}) = convert(BFloat16, randn(rng))
randexp(rng::AbstractRNG, ::Type{BFloat16}) = convert(BFloat16, randexp(rng))

# Exponent
exponent(x::BFloat16) = exponent(Float32(x))

# Bitstring
bitstring(x::BFloat16) = bitstring(reinterpret(UInt16, x))


text:
function mul!(z::acb, x::acb, y::arb)
  ccall((:acb_mul_arb, libarb), Nothing, (Ref{acb}, Ref{acb}, Ref{arb}, Int),
        z, x, y, parent(z).prec)
  return z
end

################################################################################
#
#  Small interface to characters
#
################################################################################

function Base.show(io::IO, C::RCFCharacter)
  println(IOContext(io, :compact => true), "Character of $(C.C)")
end

function character(C::ClassField, x::GrpAbFinGenElem, mGhat::Map)
  return RCFCharacter(C, x, mGhat)
end

function _conjugate(chi::RCFCharacter)
  return character(chi.C, -chi.x, chi.mGhat)
end

function iszero(chi::RCFCharacter)
  return iszero(chi.x)
end

function ==(x::RCFCharacter, y::RCFCharacter)
  if x.C === y.C
    return x.x == y.x
  end
  error("Not yet implemented!")
end

function conductor(chi::RCFCharacter)
  if isdefined(chi, :conductor)
    return chi.conductor
  end
  C = chi.C
  x = chi.x
  mGhat = chi.mGhat
  mp = mGhat(x) #The character as a map
  k, mk = kernel(mp)
  q, mq = cokernel(mk, false)
  C1 = ray_class_field(C.rayclassgroupmap, C.quotientmap*mq)
  chi.conductor, chi.conductor_inf_plc = conductor(C1)
  return chi.conductor
end

(chi::RCFCharacter)(I::NfOrdIdl, prec::Int) = image(chi, I, prec)
(chi::RCFCharacter)(I::GrpAbFinGenElem, prec::Int) = image(chi, I, prec)

function image(chi::RCFCharacter{MapClassGrp, T}, I::NfOrdIdl, prec::Int) where T
  CC = AcbField(prec)
  x = chi.x
  mGhat = chi.mGhat
  mpc = mGhat(x)
  C = chi.C
  mR = C.rayclassgroupmap
  mQ = C.quotientmap
  img = lift(mpc(mQ(mR\I)))
  if iszero(img)
    return one(CC)
  end
  return cispi(2*CC(img))
end

function image(chi::RCFCharacter{MapRayClassGrp, T}, I::NfOrdIdl, prec::Int) where T
  CC = AcbField(prec)
  x = chi.x
  mGhat = chi.mGhat
  mpc = mGhat(x)
  c = conductor(chi)
  dc = defining_modulus(chi.C)[1]
  if c == dc
    if !is_coprime(I, c)
      return zero(CC)
    end
    C = chi.C
    mR = C.rayclassgroupmap
    if !isdefined(mR, :prime_ideal_preimage_cache)
      mR.prime_ideal_preimage_cache = Dict{NfOrdIdl, GrpAbFinGenElem}()
    end
    if haskey(mR.prime_ideal_preimage_cache, I)
      pim = mR.prime_ideal_preimage_cache[I]
    else
      pim = mR\I
      mR.prime_ideal_preimage_cache[I] = pim
    end
    mQ = C.quotientmap
    img = lift(mpc(mQ(pim)))
    if iszero(img)
      return one(CC)
    end
    return cispi(2*CC(img))
  end
  assure_with_conductor(chi)
  if !is_coprime(I, conductor(chi))
    return zero(CC)
  end
  mR = chi.mrcond
  mp = chi.mp_cond
  mQ = chi.C.quotientmap
  if !isdefined(mR, :prime_ideal_preimage_cache)
    mR.prime_ideal_preimage_cache = Dict{NfOrdIdl, GrpAbFinGenElem}()
  end
  if haskey(mR.prime_ideal_preimage_cache, I)
    el = mR.prime_ideal_preimage_cache[I]
  else
    el = mR\I
    mR.prime_ideal_preimage_cache[I] = el
  end
  el = mp\el
  el = mQ(el)
  el = mpc(el)
  img = lift(el)
  if iszero(img)
    return one(CC)
  end
  return cispi(2*CC(img))
end

function image(chi::RCFCharacter, x::GrpAbFinGenElem, prec::Int)
  CC = AcbField(prec)
  mp = chi.mGhat(chi.x)
  img = lift(mp(x))
  if iszero(img)
    return one(CC)
  end
  return cispi(CC(2*img))
end

function assure_with_conductor(chi::RCFCharacter)
  if isdefined(chi, :mrcond)
    return nothing
  end
  C = chi.C
  mR = C.rayclassgroupmap
  c = conductor(chi)
  r, mr = ray_class_group(c, chi.conductor_inf_plc, n_quo = degree(C))
  lp, sR = find_gens(mR)
  imgs = GrpAbFinGenElem[mr\x for x in lp]
  mpR = hom(sR, imgs)
  chi.mrcond = mr
  chi.mp_cond = mpR
  return nothing
end

################################################################################
#
#  RCF using Stark units
#
################################################################################

function rcf_using_stark_units(C::T; cached::Bool = true) where T <: ClassField_pp
  K = base_field(C)
  @hassert :ClassField 1 is_totally_real(K)
  c, inf_plc = conductor(C)
  @hassert :ClassField 1 isempty(inf_plc)
  C1, mp, v = _find_suitable_quadratic_extension(C)
  kmp, mkmp = kernel(mp)
  comp = mkmp*C1.quotientmap
  imgc, mimgc = image(comp)
  @hassert :ClassField 1 order(imgc) == 2
  simgs, msimgs = snf(imgc)
  y = mimgc(msimgs(simgs[1]))
  Ghat, mGhat = dual(codomain(C1.quotientmap))
  p = 64
  #I don't need the character with value 1 on the generator of the
  #quadratic extension
  chars = RCFCharacter{MapRayClassGrp, GrpAbFinGenMap}[character(C1, x, mGhat) for x in Ghat if !iszero(lift(mGhat(x)(y)))]
  @hassert :ClassField 1 length(chars) == degree(C)
  #First, we get an approximation of the defpoly at v
  fl = false
  nb = 1
  while !fl
    approximations_derivative_Artin_L_functions = approximate_derivative_Artin_L_function(chars, p)
    el = approximate_artin_zeta_derivative_at_0(C1, mp, approximations_derivative_Artin_L_functions)
    fl, nb = approximate_defpoly(K, el)
    if !fl
      p *= 2
      @vprint :ClassField 1 "Guess not useful. Increasing precision to $(2*p) \n"
    elseif nb < p
      f = find_defining_polynomial(K, el, v)
      if degree(f) != degree(C)
        continue
      end
      mR = C.rayclassgroupmap
      mQ = C.quotientmap
      ng, mng = norm_group(f, mR, false, cached = false, check = false, of_closure = false)
      if iszero(mng*mQ)
        C.A = number_field(f, cached = false, check = false)[1]
        return nothing
      end
    end
  end
  p = max(64+p, nb)
  @vprint :ClassField 1 "Guess: we need precision $(p) \n"
  while true
    approximations_derivative_Artin_L_functions = approximate_derivative_Artin_L_function(chars, p)
    el = approximate_artin_zeta_derivative_at_0(C1, mp, approximations_derivative_Artin_L_functions)
    f = find_defining_polynomial(K, el, v)
    if degree(f) != degree(C)
      p = 2*p
      @vprint :ClassField 1 "Wrong guess, setting precision to $(p) \n"
      continue
    end
    mR = C.rayclassgroupmap
    mQ = C.quotientmap
    ng, mng = norm_group(f, mR, false, cached = false, check = false, of_closure = false)
    if iszero(mng*mQ)
      C.A = number_field(f, cached = false, check = false)[1]
      return nothing
    end
    p *= 2
    @vprint :ClassField 1 "The reconstructed polynomial is wrong, setting precision to $(p) \n"
  end
end

function power_sums(v::Vector{T}) where T <: FieldElem
  res = Vector{T}(undef, length(v))
  pows = T[x for x in v]
  res[1] = sum(pows)
  for i = 2:length(v)
    for j = 1:length(pows)
      pows[j] *= v[j]
    end
    res[i] = sum(pows)
  end
  return res
end

function approximate_defpoly(K::AnticNumberField, el::Vector{acb})
  rts = arb[2*real(cosh(x)) for x in el]
  ps = power_sums(rts)
  pol = power_sums_to_polynomial(ps)
  @vprint :ClassField 1 "Approximation: $pol \n"
  l2norm = fmpz(0)
  for i = 0:degree(pol)
    c = coeff(pol, i)
    l2norm += upper_bound(fmpz, abs(c)^2)
  end
  return true, nbits(l2norm)
end

function find_defining_polynomial(K::AnticNumberField, el::Vector{acb}, v::InfPlc)
  OK = maximal_order(K)
  rts = arb[2*real(cosh(x)) for x in el]
  ps = power_sums(rts)
  pol = power_sums_to_polynomial(ps)
  attempt = _find_coeffs(K, pol, v)
  if degree(attempt) == length(el)
    return attempt
  end
  rt1 = arb[x^2-2 for x in rts]
  ps = power_sums(rt1)
  pol = power_sums_to_polynomial(ps)
  return _find_coeffs(K, pol, v)
end

function _find_coeffs(K, pol, v)
  Kt = PolynomialRing(K, "t", cached = false)[1]
  OK = maximal_order(K)
  B = basis(OK, K)
  bconjs = [real(evaluate(x, v, 2*precision(parent(pol)))) for x in B]
  coeffs = Vector{nf_elem}(undef, degree(pol)+1)
  for i = 1:degree(pol)
    c = coeff(pol, i-1)
    bn = 3*nbits(Hecke.upper_bound(fmpz, c))
    fl, comb = _approximate(c, bconjs, bn)
    if !fl
      add = 10
      while !fl && add < 100
        fl, comb = _approximate(c, bconjs, bn)
        add += 10
      end
      return Kt()
    end
    coeffs[i] = dot(comb, B)
  end
  coeffs[end] = one(K)
  return Kt(coeffs)
end

function _approximate(el::arb, A::Vector{arb}, bits::Int)
  n = length(A)
  V0 = floor(ldexp(el, bits) + 0.5)
  V = [floor(ldexp(s, bits) + 0.5) for s in A]
  M = zero_matrix(ZZ, n + 1, n + 2)
  for i = 1:n
    M[i, i] = 1
    flag, M[i, n + 2] = unique_integer(V[i])
    !flag && return false, fmpz[M[1, 1]]
  end
  M[n + 1, n + 1] = 1
  flag, M[n + 1, n + 2] = unique_integer(-V0)
  !flag && return false, fmpz[M[1, 1]]
  lll!(M, lll_ctx(0.99999, 0.5001))
  if !isone(abs(M[1, n + 1]))
    return false, fmpz[M[1, 1]]
  end
  if isone(M[1, n + 1])
    res = fmpz[M[1, i] for i = 1:n]
  else
    res = fmpz[-M[1, i] for i = 1:n]
  end
  return true, res
end

function _find_suitable_quadratic_extension(C::T) where T <: ClassField_pp
  c = factored_modulus(C)
  K = base_field(C)
  mR = C.rayclassgroupmap
  mQ = C.quotientmap
  R = codomain(mQ)
  OK = maximal_order(K)
  real_plc = real_places(K)
  @hassert :ClassField 1 length(real_plc) == degree(K)
  v = real_plc[end]
  w = real_plc[1:end-1]
  inf_plc_ram = Set(w)
  bound = 10
  ctx = rayclassgrp_ctx(OK, Int(exponent(C))*2)
  allow_cache!(ctx.class_group_map)
  lc = ideals_up_to(OK, bound, conductor(C)[1])
  cnt = 0
  while true
    @vprint :ClassField 1 "Batch of ideals with $(length(lc)) elements \n"
    for I in lc
      for i = 1:length(real_plc)
        v = real_plc[i]
        w = [real_plc[j] for j = 1:length(real_plc) if j != i]
        lf = factor(I)
        newc = merge(max, c, lf)
        r, mr = ray_class_group_quo(OK, newc, w, ctx)
        gens, group_gens = find_gens(mr)
        images = GrpAbFinGenElem[mQ(mR\J) for J in gens]
        mp = hom(group_gens, images, check = false)
        k, mk = kernel(mp)
        ls = subgroups(k, quotype = Int[2], fun = (x, y) -> sub(x, y, false)[2])
        for ms in ls
          ms::GrpAbFinGenMap
          q, mq = cokernel(ms*mk, false)
          Cnew = ray_class_field(mr, mq)
          cnew, inf_plcc = conductor(Cnew)
          if length(inf_plcc) != length(real_plc)-1
            continue
          end
          acceptable = true
          for (P, v) in c
            pdt1 = prime_decomposition_type(C, P)
            pdt2 = prime_decomposition_type(Cnew, P)
            if pdt1[3] != pdt2[3]
              acceptable = false
              break
            end
          end
          if acceptable
            return Cnew, mp, v
          end
        end
      end
    end
    #Need to consider more moduli.
    bound *= 2
    lc1 = ideals_up_to(OK, bound, conductor(C)[1])
    lc = setdiff(lc1, lc)
  end
end

function approximate_artin_zeta_derivative_at_0(C::ClassField, mp::GrpAbFinGenMap, D::Dict{S, T}) where {S, T}
  Dzeta = Vector{acb}()
  ks = keys(D)
  CC = parent(first(values(D)))
  R = codomain(C.quotientmap)
  mp1 = hom(R, codomain(mp), GrpAbFinGenElem[mp(C.quotientmap\x) for x in gens(R)])
  k1, mk1 = kernel(mp1)
  ck, mck = cokernel(mk1)
  for x1 in ck
    x = mck\x1
    v = zero(CC)
    for k in ks
      v += D[k]*conj(k(x, CC.prec))
    end
    push!(Dzeta, v//(degree(C)))
  end
  return Dzeta
end

function approximate_derivative_Artin_L_function(chars::Vector{RCFCharacter{MapRayClassGrp, GrpAbFinGenMap}}, prec::Int)
  chars1 = Vector{RCFCharacter{MapRayClassGrp, GrpAbFinGenMap}}()
  idxs = Tuple{Int, Bool}[]
  for i = 1:length(chars)
    idx = 0
    cx = _conjugate(chars[i])
    for j = 1:length(chars1)
      if cx == chars1[j]
        idx = j
        break
      end
    end
    if iszero(idx)
      #The conjugate is not there.
      push!(chars1, chars[i])
      push!(idxs, (length(chars1), false))
    else
      push!(idxs, (idx, true))
    end
  end
  coeffs1 = _approximate_derivative_Artin_L_function(chars1, prec)
  if length(chars1) == length(chars)
    return Dict{RCFCharacter{MapRayClassGrp, GrpAbFinGenMap}, acb}(zip(chars, coeffs1))
  end
  coeffs = Dict{RCFCharacter{MapRayClassGrp, GrpAbFinGenMap}, acb}()
  for i = 1:length(chars)
    if idxs[i][2]
      coeffs[chars[i]] = conj(coeffs1[idxs[i][1]])
    else
      coeffs[chars[i]] = coeffs1[idxs[i][1]]
    end
  end
  return coeffs
end

function _limx(K::AnticNumberField, target_prec::Int)
  R = ArbField(3*target_prec)
  r1, r2 = signature(K)
  n = degree(K)
  r = r1 + r2
  c_1 = n*R(2)^fmpq(-2*r2, n)
  p1 = (2*const_pi(R))^(r-1)
  p2 = n^r*p1//2^r2
  c_0 = sqrt(p2//(c_1^(r+1)))
  A0 = log(c_0*(fmpz(2)^target_prec))
  p2 = A0//c_1
  p1 = n*(r+1)*log(p2)//(2*(r+1)+4*A0)
  return p2^(fmpq(n, 2))//(p1+1)
end

function _find_N0(K::AnticNumberField, target_prec::Int, C::arb)
  N0 = C*_limx(K, target_prec)
  return N0
end

function _B(K::AnticNumberField, target_prec::Int)
  R = ArbField(3*target_prec)
  return degree(K)*R(2)^(target_prec)//sqrt(const_pi(R)^3*_limx(K, target_prec))
end

function _find_i0(K::AnticNumberField, target_prec::Int)
  B = _B(K, target_prec)
  limx = _limx(K, target_prec)
  i0 = 300
  while limx^i0*(factorial(fmpz(div(i0, 2)))) < B
    i0 *= 2
  end
  imax = i0
  imin = div(i0, 2)
  while imax - imin > 2
    attempt = div(imax + imin, 2)
    if limx^attempt*(factorial(fmpz(div(attempt, 2)))) < B
      imin = attempt
    else
      imax = attempt
    end
  end
  return imax
end

function _approximate_derivative_Artin_L_function(chars::Vector, target_prec::Int)
  K = base_field(chars[1].C)
  n = degree(K)
  if n == 2
   return  _approximate_derivative_Artin_L_function_quadratic(chars, target_prec)
  end
  prec = min(10, div(degree(chars[1].C), 2))*target_prec
  RR = ArbField(prec)
  maxC = (root(norm(conductor(chars[1].C)[1])*abs(discriminant(maximal_order(K))), 2)+1)//(sqrt(const_pi(RR))^degree(K))
  nterms = Int(Hecke.upper_bound(fmpz, target_prec*maxC//2))
  i0 = _find_i0(K, target_prec)
  Acoeffs = _compute_A_coeffs(n, i0, prec)
  factorials = Vector{fmpz}(undef, n)
  factorials[1] = fmpz(1)
  for i = 2:length(factorials)
    factorials[i] = factorials[i-1]*(i-1)
  end
  coeffs_0 = Vector{arb}[_Aij_at_0(i, n, Acoeffs[i+1], factorials) for i = 0:i0]
  coeffs_1 = Vector{arb}[_Aij_at_1(i, n, Acoeffs[i+1], factorials) for i = 0:i0]
  den = 2*sqrt(const_pi(RR)^(n-1))
  res = Vector{acb}()
  for i = 1:length(chars)
    x = chars[i]
    A = _A_function(x, prec)
    lambda = _lambda_and_artin(x, target_prec, coeffs_0, coeffs_1, compute_coeffs_L_function(x, nterms, prec))
    num = A*lambda
    push!(res, num/den)
  end
  return res
end

function _approximate_derivative_Artin_L_function_quadratic(chars::Vector, target_prec::Int)
  prec = min(10, div(degree(chars[1].C), 2))*target_prec
  RR = ArbField(prec)
  res = Vector{acb}()
  den = 2*sqrt(const_pi(RR))
  @vtime :ClassField 1 vf = compute_values_f_quadratic(chars, target_prec)
  for i = 1:length(chars)
    x = chars[i]
    A = _A_function(x, prec)
    lambda = _lambda_and_artin_quadratic(x, target_prec, vf)
    num = A*lambda
    push!(res, num/den)
  end
  return res
end

function compute_values_f_quadratic(chars::Vector, target_prec::Int)
  prec = min(10, div(degree(chars[1].C), 2))*target_prec
  res = Dict{NfOrdIdl, Vector{Tuple{arb, arb}}}()
  for x in chars
    c = conductor(x)
    if haskey(res, c)
      continue
    end
    cx = _C(x, prec)
    nterms = Int(Hecke.upper_bound(fmpz, (target_prec*cx)//2))
    v = Vector{Tuple{arb, arb}}(undef, nterms)
    v0 = _evaluate_f_x_0(cx, prec, 2*target_prec, nterms)
    for i = 1:length(v)
      val = cx//i
      v1 = _evaluate_f_x_1(val, prec)
      v[i] = (v0[i], v1)
    end
    res[c] = v
  end
  return res
end

################################################################################
#
#  Coefficients L-function
#
################################################################################

function compute_coeffs_L_function(chi::T, n::Int, prec::Int) where T <: RCFCharacter
  C = chi.C
  OK = base_ring(C)
  lp = prime_ideals_up_to(OK, n)
  sort!(lp, by = x -> norm(x, copy = false))
  CC = AcbField(prec)
  coeffs_old = Vector{acb}(undef, n)
  coeffs_old[1] = one(CC)
  for j = 2:n
    coeffs_old[j] =  zero(CC)
  end
  coeffs_new = Vector{acb}(undef, n)
  for h = 1:length(lp)
    P = lp[h]
    r = chi(P, prec)
    np = Int(norm(P, copy = false))
    j = np
    while j <= n
      v = valuation(j, np)
      res = zero(CC)
      if iszero(r)
        coeffs_new[j] = res
        j += np
        continue
      end
      r1 = one(CC)
      aux = CC()
      for k = 0:v
        kk = divexact(j, np^k)
        ckk = coeffs_old[kk]
        mul!(aux, ckk, r1)
        add!(res, res, aux)
        #res += ckk*r1
        mul!(r1, r1, r)
      end
      coeffs_new[j] = res
      j += np
    end
    #Now, I need to shift the column...
    j = np
    while j <= n
      coeffs_old[j] = coeffs_new[j]
      j += np
    end
  end
  return coeffs_old
end

function ideals_up_to(OK::NfOrd, n::Int, coprime_to::NfOrdIdl = ideal(OK, 1))

  lp = prime_ideals_up_to(OK, n)
  filter!(x -> is_coprime(x, coprime_to), lp)
  lI = NfOrdIdl[ideal(OK, 1)]
  for i = 1:length(lp)
    lnew = NfOrdIdl[]
    P = lp[i]
    nP = Int(norm(P, copy = false))
    @assert nP <= n
    expon = Int(flog(fmpz(n), nP))
    for j = 1:length(lI)
      I = lI[j]
      if norm(I, copy = false)*nP > n
        break
      end
      push!(lnew, I*P)
      for s = 2:expon
        if nP^s*norm(I, copy = false) > n
          break
        end
        push!(lnew, I*P^s)
      end
    end
    append!(lI, lnew)
    sort!(lI, by = x -> norm(x, copy = false))
  end
  return lI
end

################################################################################
#
#  Correction term
#
################################################################################

function _A_function(chi::RCFCharacter, prec::Int)
  R = AcbField(prec)
  if iszero(chi)
    #We have the trivial character sending everything to 1.
    #we return 0
    return zero(R)
  end
  C = chi.C
  cC = conductor(C)[1]
  cchi = conductor(chi)
  if cchi == cC
    return one(R)
  end
  fcC = factor(cC)
  fcchi = factor(cchi)
  res = one(R)
  for (p, v) in fcC
    if haskey(fcchi, p)
      continue
    end
    res = mul!(res, res, 1 - chi(p, prec))
  end
  return res
end

################################################################################
#
#  Artin root number
#
################################################################################

function artin_root_number(chi::RCFCharacter, prec::Int)
  R = AcbField(prec)
  C = chi.C
  c = conductor(chi)
  OK = base_ring(C)
  D = different(OK)
  J = D*c
  lfJ = factor(J)
  lambda = OK(approximate(collect(values(lfJ)), collect(keys(lfJ))))
  lambda = make_positive(lambda, minimum(J)^2)
  @hassert :ClassField 1 is_totally_positive(lambda)
  @hassert :ClassField 1 lambda in J
  g = numerator(simplify(ideal(OK, lambda) * inv(J)))
  @hassert :ClassField 1 is_coprime(g, c)
  u = idempotents(g, c)[1]
  u = make_positive(u, minimum(g^2))
  @hassert :ClassField 1 is_totally_positive(u)
  @hassert :ClassField 1 u in g
  h = numerator(simplify(ideal(OK, u) * inv(g)))
  @hassert :ClassField 1 is_coprime(h, c)
  Qc, mQc = quo(OK, c)
  G, mG = multiplicative_group(Qc)
  reps = NfOrdElem[make_positive(lift(mG(x)), minimum(c)^2) for x in G]
  for (j, x) in enumerate(G)
    @hassert :ClassField 1 x == mG\Qc(reps[j])
    @hassert :ClassField 1 is_totally_positive(reps[j])
  end
  Gsum = R(0)
  for i = 1:length(reps)
    el = reps[i].elem_in_nf*u.elem_in_nf//lambda.elem_in_nf
    trel = 2*tr(el)
    newtrel = fmpq(mod(numerator(trel), 2*denominator(trel)), denominator(trel))
    expi = cispi(R(newtrel))
    Gsum += chi(ideal(OK, reps[i]), prec) * expi
  end
  Gsum = mul!(Gsum, Gsum, chi(h, prec))
  res = (-onei(R))^length(chi.conductor_inf_plc)*Gsum/sqrt(R(norm(c)))
  return res
end

################################################################################
#
#  Lambda evaluation for quadratic fields
#
################################################################################

function _lambda_and_artin_quadratic(chi::RCFCharacter, target_prec::Int, vf::Dict{NfOrdIdl, Vector{Tuple{arb, arb}}})
  prec = min(10, div(degree(chi.C), 2))*target_prec
  Wchi = artin_root_number(chi, prec)
  CC = AcbField(prec)
  c = conductor(chi)
  vals = vf[c]
  nterms_chi = length(vals)
  @vprint :ClassField 1 "Computing $(nterms_chi) coefficients Artin L function \n"
  @vtime :ClassField 1 val_chi = compute_coeffs_L_function(chi, nterms_chi, prec)
  res1 = zero(CC)
  res2 = zero(CC)
  aux = CC()
  for i = 1:nterms_chi
    if iszero(val_chi[i])
      continue
    end
    mul!(aux, val_chi[i], vals[i][1])
    add!(res1, res1, aux)
    #res1 += val_chi[i]*vals[i][1]
    mul!(aux, conj(val_chi[i]), vals[i][2])
    add!(res2, res2, aux)
    #res2 += conj(val_chi[i])*vals[i][2]
  end
  return res1 + Wchi*res2
end

function _evaluate_f_x_1(x::arb, prec::Int)
  RR = parent(x)
  return x*sqrt(const_pi(RR))*exp(-2//x)
end

function _evaluate_f_x_0(x::arb, prec::Int)
  RR = parent(x)
  CC = AcbField(precision(RR))
  return 2*sqrt(const_pi(RR))*real(exp_integral_e(one(CC), CC(2//x)))
end

function _evaluate_f_x_0(x::arb, prec::Int, tolerance::Int, N::Int)
  CC = AcbField(prec)
  RR = ArbField(prec)
  res = Vector{arb}(undef, N)
  A = 2//x
  res[N] = real(exp_integral_e(one(CC), CC(N*A)))
  nstop = Int(upper_bound(fmpz, ceil(4//A)))
  n = N
  e0 = exp(A)
  e1 = exp(-N*A)
  th = fmpq(1, 2)^tolerance
  while n > nstop
    first = true
    res[n-1] = zero(RR)
    f0 = e1
    f1 = -f0//n
    m = 1
    d = fmpq(-1)
    s = res[n]
    while abs(s) > th
      add!(res[n-1], res[n-1], s)
      #res[n-1] += s
      s = d*f1
      f0 = -A*f0
      f1 = -(m*f1+f0)//n
      m += 1
      d = -d//m
    end
    n = n-1
    e1 = e0*e1
  end
  for i = 1:nstop
    res[i] = real(exp_integral_e(one(CC), CC(2*i//x)))
  end
  cc = 2*sqrt(const_pi(RR))
  for i = 1:N
    mul!(res[i], res[i], cc)
    #res[i] = cc*res[i]
  end
  return res
end

################################################################################
#
#  Lambda evaluation - general case
#
################################################################################

function _lambda_and_artin(chi::RCFCharacter, target_prec::Int, coeffs_0, coeffs_1, coeffs_chi)
  prec = min(10, div(degree(chi.C), 2))*target_prec
  K = base_field(chi.C)
  Wchi = artin_root_number(chi, prec)
  CC = AcbField(prec)
  res1 = zero(CC)
  res2 = zero(CC)
  cchi = _C(chi, prec)
  n = degree(K)
  nterms_chi = Int(Hecke.upper_bound(fmpz, target_prec*cchi//2))
  res1 = zero(CC)
  res2 = zero(CC)
  aux = CC()
  for i = 1:nterms_chi
    if iszero(coeffs_chi[i])
      continue
    end
    evpoint = cchi//i
    ev0, ev1 = _evaluate_f_x_0_1(evpoint, n, target_prec, coeffs_0, coeffs_1)
    mul!(aux, coeffs_chi[i], ev0)
    add!(res1, res1, aux)
    mul!(aux, conj(coeffs_chi[i]), ev1)
    add!(res2, res2, aux)
  end
  return res1 + Wchi*res2
end

function _C(chi::RCFCharacter, prec::Int)
  RR = ArbField(prec)
  c = conductor(chi)
  OK = order(c)
  nc = norm(c)
  p = const_pi(RR)^degree(OK)
  d = sqrt(RR(abs(discriminant(OK))))*sqrt(RR(nc))
  return d//sqrt(p)
end

function _evaluate_f_x_0_1(x::arb, n::Int, target_prec::Int, coeffs_0::Vector{Vector{arb}}, coeffs_1::Vector{Vector{arb}})
  RR = ArbField(3*target_prec)
  res0 = zero(RR)
  res1 = zero(RR)
  lnx = log(x)
  powslogx = powers(lnx, n-1)
  ix = inv(x)
  th = RR(fmpq(1, 2)^(3+target_prec))
  ixpow = RR(1)
  aux = RR()
  for i = 0:length(coeffs_0)-1
    if iseven(i)
      m = 1
    else
      m = n-1
    end
    aij1 = coeffs_1[i+1]
    auxres1 = zero(RR)
    for j = 2:m+1
      mul!(aux, aij1[j], powslogx[j-1])
      add!(auxres1, auxres1, aux)
      #auxres1 += aij1[j]*powslogx[j-1]
    end
    auxres1 = mul!(auxres1, auxres1, ixpow)
    res1 = add!(res1, res1, auxres1)

    aij0 = coeffs_0[i+1]
    auxres0 = zero(RR)
    for j = 2:length(aij0)
      mul!(aux, aij0[j], powslogx[j-1])
      add!(auxres0, auxres0, aux)
      #auxres0 += aij0[j]*powslogx[j-1]
    end
    auxres0 = mul!(auxres0, auxres0, ixpow)
    res0 = add!(res0, res0, auxres0)
    if iszero(mod(i, 15)) && abs(auxres0) < th && abs(auxres1) < th
      break
    end
    #res0 += ixpow*auxres0
    ixpow = mul!(ixpow, ixpow, ix)
  end
  res1 += x*gamma(fmpq(1, 2), RR)*gamma(fmpz(1), RR)^(n-1)

  #CC = AcbField(prec)
  #res0int = (one(CC)/(2*const_pi(CC)*onei(CC)))*Nemo.integrate(CC, y -> x^y * gamma(y/2)*gamma((y+1)/2)^(n-1)/y, 1.1 - (nterms+0.1)*onei(CC), 1.1 + (nterms+0.1) * onei(CC))
  #res1int = (one(CC)/(2*const_pi(CC)*onei(CC)))*Nemo.integrate(CC, y -> x^y * gamma(y/2)*gamma((y+1)/2)^(n-1)/(y-1), 1.1 - (nterms+0.1)*onei(CC), 1.1 + (nterms+0.1) * onei(CC))
  #@assert overlaps(real(res0int), res0)
  #@assert overlaps(real(res1int), res1)
  #@show res0
  #@show res0int
  #return real(res0int), real(res1int)
  return res0, res1
end

function _Aij_at_0(i::Int, n::Int, aij::Vector{arb}, factorials::Vector{fmpz})
  #aij starts with ai0 and finishes with ain
  CC = parent(aij[1])
  if iseven(i)
    m = 1
  else
    m = n-1
  end
  if iszero(i)
    D = aij[1:m+2]
  else
    D = Vector{arb}(undef, m+1)
    D[m+1] = -aij[m+1]/i
    for j = m:-1:1
      D[j] = (D[j+1] - aij[j])/i
    end
  end
  for i = 2:length(D)
    D[i] = D[i]//factorials[i-1]
  end
  #=
  ev = CC(-i)+CC(0.001)
  @show val1 = (gamma(ev/2)*gamma((ev+1)/2)^(n-1))/ev
  @show val2 = sum(D[j+1]/((ev+i)^j) for j = 0:m)
  @hassert :ClassField 1 radiuslttwopower(abs(val1-val2), 16)
  =#
  return D
end

function _Aij_at_1(i::Int, n::Int, aij::Vector{arb}, factorials::Vector{fmpz})
  #aij starts with ai0 and finishes with ain
  if iseven(i)
    m = 1
  else
    m = n-1
  end
  if iszero(i)
    aij = aij[2:end]
  end
  CC = parent(aij[1])
  D = Vector{arb}(undef, m+1)
  D[m+1] = -aij[m+1]/(i+1)
  for j = m:-1:1
    D[j] = (D[j+1] - aij[j])/(i+1)
  end
  for i = 2:length(D)
    D[i] = D[i]//factorials[i-1]
  end
  #=
  ev = CC(-i)+CC(0.001)
  val1 = (gamma(ev/2)*gamma((ev+1)/2)^(n-1))/(ev-1)
  val2 = sum(D[j+1]/((ev+i)^j) for j = 0:m)
  @hassert :ClassField 1 radiuslttwopower(abs(val1-val2), 8)
  =#
  return D
end

function _compute_A_coeffs(n::Int, nterms::Int, prec::Int)
  RR = ArbField(prec)
  #res[j] contains A_{i,(j-1)}
  res_final = Vector{Vector{arb}}(undef, nterms+1)
  spi = sqrt(const_pi(RR))
  coeffs_exp_even = _coeffs_exp_even(n, nterms, RR)
  coeffs_exp_odd = _coeffs_exp_odd(n, nterms, RR)
  spipow = spi^(n-1)
  for i = 0:nterms
    if iszero(i)
      #In this case, I need the i-1st coefficient
      res = Vector{arb}(undef, 3)
      r = spipow*_coeff_0_even(n, 0)
      r1 = _coeff_exp_0(n, RR)
      #Careful: in this case the indices are shifted.
      mul!(r1[3], r1[3], r)
      res[1] = r1[3]
      mul!(r1[2], r1[2], r)
      res[2] = r1[2]
      mul!(r1[1], r1[1], r)
      res[3] = r1[1]
    elseif iseven(i)
      res = Vector{arb}(undef, 2)
      q = divexact(i, 2)
      r = spipow*_coeff_0_even(n, q)
      res[2] = r
      r1 = coeffs_exp_even[q]
      res[1] = r*r1
    else
      res = Vector{arb}(undef, n+1)
      q = divexact(i-1, 2)
      r0 = spi*_coeff_0_odd(n, q)
      vg = coeffs_exp_odd[q+1]
      res[n+1] = zero(RR)
      for j = 1:n
        mul!(vg[n-j+1], vg[n-j+1], r0)
        res[j] = vg[n-j+1]
      end
    end
    res_final[i+1] = res
  end
  return res_final
end

function _coeff_0_odd(n::Int, q::Int)
  exc_num = fmpz(2)^(n+2*q)
  exc_den = (2*q+1)*factorial(fmpz(2*q))*factorial(fmpz(q))^(n-2)
  r = fmpq(exc_num, exc_den)
  if iseven(q) || iseven(n)
    return -r
  else
    return r
  end
end

function _coeff_exp_odd(n::Int, q::Int, RR::ArbField)
  res = Vector{arb}(undef, n-1)
  res[1] = _sum_pow_inv_odd(q, 1) + (n-1)*_sum_pow_inv_even(q, 1) - log(RR(2)) - fmpq(n, 2)*const_euler(RR)
  for k = 2:n-1
    res[k] = (-1)^k*zeta(k, RR)*(1+fmpq(n-2, fmpz(2)^k)) + _sum_pow_inv_odd(q, k) + (n-1)*_sum_pow_inv_even(q, k)
    res[k] = res[k]/k
  end
  RRx = PowerSeriesRing(RR, n, "x", cached = false)[1]
  g = RRx(res, length(res), n, 1)
  gexp = exp(g)
  return arb[coeff(gexp, i) for i = 0:n-1]
end

function _coeffs_exp_odd(n::Int, nterms::Int, RR::ArbField)
  nt = div(nterms, 2)
  if isodd(nterms)
    nt += 1
  end
  RRx = PowerSeriesRing(RR, n, "x", cached = false)[1]
  res = Vector{Vector{arb}}(undef, nt)
  sum_pow_inv_odd = Vector{fmpq}(undef, n-1)
  sum_pow_inv_even = Vector{fmpq}(undef, n-1)
  for i = 1:n-1
    sum_pow_inv_odd[i] = fmpq(0)
    sum_pow_inv_even[i] = fmpq(0)
  end
  inexc =  log(RR(2)) + fmpq(n, 2)*const_euler(RR)
  zeta_k = Vector{arb}(undef, n-2)
  for i = 1:n-2
    if iseven(i)
      zeta_k[i] = -zeta(i+1, RR)
    else
      zeta_k[i] = zeta(i, RR)
    end
  end
  for q = 0:nt-1
    for j = 1:n-1
      sum_pow_inv_odd[j] += fmpq(1, (2*q+1)^j)
      if !iszero(q)
        sum_pow_inv_even[j] += (n-1)*fmpq(1, (2*q)^j)
      end
    end
    res_q = Vector{arb}(undef, n-1)
    res_q[1] = sum_pow_inv_odd[1] + sum_pow_inv_even[1] - inexc
    for k = 2:n-1
      res_q[k] = zeta_k[k-1]*(1+fmpq(n-2, fmpz(2)^k)) + sum_pow_inv_odd[k] + sum_pow_inv_even[k]
      res_q[k] = res_q[k]/k
    end
    g = RRx(res_q, length(res_q), n, 1)
    gexp = exp(g)
    res[q+1] = arb[coeff(gexp, i) for i = 0:n-1]
  end
  return res
end

function _coeffs_exp_even(n::Int, nterms::Int, RR::ArbField)
  inexc = fmpq(n, 2)*const_euler(RR)+(n-1)*log(RR(2))
  res = Vector{arb}(undef, div(nterms, 2))
  sum_inv_odd = fmpq(0)
  sum_inv_even = fmpq(1, 2)
  res[1] = sum_inv_even - inexc
  for q = 2:length(res)
    sum_inv_odd += (n-1)*fmpq(1, (2*q-1))
    sum_inv_even += fmpq(1, 2*q)
    res[q] = sum_inv_odd + sum_inv_even - inexc
  end
  return res
end

function _coeff_exp_0(n::Int, RR::ArbField)
  c0 = -fmpq(n, 2)*const_euler(RR)-(n-1)*log(RR(2))
  c1 = zeta(2, RR)*fmpq(3*n-2, 8)
  RRx = PowerSeriesRing(RR, 3, "x", cached = false)[1]
  g = RRx([c0, c1], 2, 3, 1)
  expg = exp(g)
  return arb[coeff(expg, i) for i = 0:2]
end

function _coeff_0_even(n::Int, q::Int)
  num = 2 * fmpz(4)^(q*(n-1))*factorial(fmpz(q))^(n-2)
  den = factorial(fmpz(2*q))^(n-1)
  r = fmpq(num, den)
  if iseven(q) || iseven(n)
    return r
  else
    return -r
  end
end