text:
@testset "BestResponsePolicy" begin
    env = OpenSpielEnv("kuhn_poker";default_state_style=RLBase.Information{String}(), is_chance_agent_required=true)
    p = TabularRandomPolicy()
    bp0 = BestResponsePolicy(p, env, 0)
    bp1 = BestResponsePolicy(p, env, 1)

    @testset "Uniform Best Response" begin
        # taken from https://github.com/deepmind/open_spiel/blob/de2c9ef57038829c2709becb0104e3b2d76725cc/open_spiel/algorithms/best_response_test.cc#L307-L315
        expected_response_0 = Dict{String, Int}("0"=>1, "0pb"=>0, "1"=>1, "1pb"=>1, "2"=>0, "2pb"=>1)
        expected_response_1 = Dict{String, Int}("0b"=>0, "0p"=>1, "1b"=>1, "1p"=>1, "2b"=>1, "2p"=>1)
        # note that for "2", the values are the same, but the lower action will be chosen since we're using `>` instead of `>=` for selecting action with a larger value

        walk(env) do x
            if get_current_player(x) == 0
                @test bp0(x) == expected_response_0[get_state(x)]
            elseif get_current_player(x) == 1
                @test bp1(x) == expected_response_1[get_state(x)]
            end
        end
    end

    @testset "corner case test" begin
        # https://github.com/deepmind/open_spiel/blob/de2c9ef57038829c2709becb0104e3b2d76725cc/open_spiel/algorithms/best_response_test.cc#L337
        kuhn_exploitability_descent_iter4_policy = TabularRandomPolicy(table=Dict(
            "0"   => [0.567034158868, 0.432965841132],
            "0b"  => [0.602000197743, 0.397999802257],
            "0p"  => [0.520821285373, 0.479178714627],
            "0pb" => [0.621126761233, 0.378873238767],
            "1"   => [0.505160629764, 0.494839370236],
            "1b"  => [0.360357968472, 0.639642031528],
            "1p"  => [0.520821285373, 0.479178714627],
            "1pb" => [0.378873238767, 0.621126761233],
            "2"   => [0.419580194883, 0.580419805117],
            "2b"  => [0.202838286881, 0.797161713119],
            "2p"  => [0.5, 0.5],
            "2pb" => [0.202838286881, 0.797161713119]
        ))

        bp1 = BestResponsePolicy(kuhn_exploitability_descent_iter4_policy, env, 1)
        expected_response_1 = Dict("0b"=>0, "0p"=>0, "1b"=> 1, "1p"=> 1, "2b"=> 1, "2p"=> 1)

        walk(env) do x
            if get_current_player(x) == 1
                @test bp1(x) == expected_response_1[get_state(x)]
            end
        end
    end

    @testset "best response for optimal policy" begin

        bp0 = BestResponsePolicy(get_optimal_kuhn_policy(), env, 0)
        bp1 = BestResponsePolicy(get_optimal_kuhn_policy(), env, 1)
        expected_response_0 = Dict("0"=>0, "0pb"=>0, "1"=>0, "1pb"=>0, "2"=>0, "2pb"=>1)
        expected_response_1 = Dict("0b"=>0, "0p"=>0, "1p"=>0, "1b"=>0, "2p"=>1, "2b"=>1)

        walk(env) do x
            if get_current_player(x) == 0
                @test bp0(x) == expected_response_0[get_state(x)]
            elseif get_current_player(x) == 1
                @test bp1(x) == expected_response_1[get_state(x)]
            end
        end
    end

    @testset "check best response value" begin

        history_and_probs_0 = Dict(
            "2"       => 1.5,  "2 1 bb"  => 2.0,  "2 1 bp"  => 1.0,  "2 1 pbp" => -1.0,
            "2 1 pb"  => 2.0,  "2 1 pp"  => 1.0,  "2 0 bb"  => 2.0,  "2 1 p"   => 1.5,
            "2 0 pp"  => 1.0,  "2 0 pbb" => 2.0,  "2 0 p"   => 1.5,  "2 1 b"   => 1.5,
            "2 0 bp"  => 1.0,  "2 1 pbb" => 2.0,  "2 0"     => 1.5,  "2 1"     => 1.5,
            "2 0 pb"  => 2.0,  "2 0 b"   => 1.5,  "2 0 pbp" => -1.0, "1 0"     => 1.5,
            "0"       => -0.5, "1 2"     => -0.5, "0 2 p"   => -1.0, ""        => 0.5,
            "0 1"     => -0.5, "0 2"     => -0.5, "1 2 pp"  => -1.0, "0 1 p"   => -1.0,
            "1"       => 0.5,  "0 2 b"   => -0.5, "1 2 pb"  => -2.0, "0 1 b"   => -0.5,
            "1 2 pbb" => -2.0, "0 2 bb"  => -2.0, "1 2 b"   => -0.5, "0 1 pb"  => -1.0,
            "1 0 bp"  => 1.0,  "0 2 pp"  => -1.0, "1 2 p"   => -1.5, "0 2 bp"  => 1.0,
            "0 1 pp"  => -1.0, "1 0 bb"  => 2.0,  "1 2 bp"  => 1.0,  "1 0 pbp" => -1.0,
            "0 1 bp"  => 1.0,  "1 0 p"   => 1.5,  "1 2 pbp" => -1.0, "0 1 pbp" => -1.0,
            "1 0 pbb" => 2.0,  "1 2 bb"  => -2.0, "0 1 bb"  => -2.0, "1 0 b"   => 1.5,
            "0 1 pbb" => -2.0, "0 2 pb"  => -1.0, "1 0 pp"  => 1.0,  "0 2 pbb" => -2.0,
            "1 0 pb"  => 2.0,  "0 2 pbp" => -1.0)

        history_and_probs_1 = Dict(
            ""        => 0.416666666667, "0"       => 1.75,  "0 1"     => 1.75,
            "0 1 b"   => 2.0,            "0 1 bb"  => 2.0,   "0 1 bp"  => -1.0,
            "0 1 p"   => 1.5,            "0 1 pb"  => 1.5,   "0 1 pbb" => 2.0,
            "0 1 pbp" => 1.0,            "0 1 pp"  => 1.0,   "0 2"     => 1.75,
            "0 2 b"   => 2.0,            "0 2 bb"  => 2.0,   "0 2 bp"  => -1.0,
            "0 2 p"   => 1.5,            "0 2 pb"  => 1.5,   "0 2 pbb" => 2.0,
            "0 2 pbp" => 1.0,            "0 2 pp"  => 1.0,   "1"       => 0.5,
            "1 0"     => -0.75,          "1 0 b"   => -1.0,  "1 0 bb"  => -2.0,
            "1 0 bp"  => -1.0,           "1 0 p"   => -0.5,  "1 0 pb"  => -0.5,
            "1 0 pbb" => -2.0,           "1 0 pbp" => 1.0,   "1 0 pp"  => -1.0,
            "1 2"     => 1.75,           "1 2 b"   => 2.0,   "1 2 bb"  => 2.0,
            "1 2 bp"  => -1.0,           "1 2 p"   => 1.5,   "1 2 pb"  => 1.5,
            "1 2 pbb" => 2.0,            "1 2 pbp" => 1.0,   "1 2 pp"  => 1.0,
            "2"       => -1.0,           "2 0"     => -0.75, "2 0 b"   => -1.0,
            "2 0 bb"  => -2.0,           "2 0 bp"  => -1.0,  "2 0 p"   => -0.5,
            "2 0 pb"  => -0.5,           "2 0 pbb" => -2.0,  "2 0 pbp" => 1.0,
            "2 0 pp"  => -1.0,           "2 1"     => -1.25, "2 1 b"   => -2.0,
            "2 1 bb"  => -2.0,           "2 1 bp"  => -1.0,  "2 1 p"   => -0.5,
            "2 1 pb"  => -0.5,           "2 1 pbb" => -2.0,  "2 1 pbp" => 1.0,
            "2 1 pp"  => -1.0)

        p = TabularRandomPolicy()
        bp0 = BestResponsePolicy(p, env, 0)
        bp1 = BestResponsePolicy(p, env, 1)

        walk(env) do x
            @test RLZoo.best_response_value(bp0, x) ≈ history_and_probs_0[string(x.state)]
            @test RLZoo.best_response_value(bp1, x) ≈ history_and_probs_1[string(x.state)]
        end
    end

    @testset "Optimal value test" begin
        history_and_probs_0 = Dict(
            ""        => -0.05555555555555558,
            "1 2 pb"  => -1.0,
            "1 2 b"   => -2.0,
            "0 2 pp"  => -1.0,
            "0 1 bp"  => 1.0,
            "2 1 bp"  => 1.0,
            "2 0 pb"  => 2.0,
            "1 2 pp"  => -1.0,
            "2 0 b"   => 1.0,
            "0 1 bb"  => -2.0,
            "2 0 pp"  => 1.0,
            "2 0 p"   => 1.3333333333333333,
            "1 0"     => 0.3333333333333333,
            "1 0 bb"  => 2.0,
            "1 0 pbp" => -1.0,
            "1 2 bp"  => 1.0,
            "2 0 bp"  => 1.0,
            "0 1"     => -1.0,
            "0 2"     => -1.0,
            "1 0 pbb" => 2.0,
            "1 0 bp"  => 1.0,
            "2 0 bb"  => 2.0,
            "1 2 bb"  => -2.0,
            "2 1"     => 1.0,
            "2 1 bb"  => 2.0,
            "2 0 pbp" => -1.0,
            "1 2 p"   => -1.0,
            "0 2 bb"  => -2.0,
            "1 0 pp"  => 1.0,
            "0 2 b"   => -2.0,
            "2 1 pb"  => 2.0,
            "1 2 pbb" => -2.0,
            "1 2"     => -1.0,
            "0 1 pb"  => -1.0,
            "0 2 p"   => -1.0,
            "0 2 bp"  => 1.0,
            "1 0 pb"  => -1.0,
            "1 2 pbp" => -1.0,
            "2 1 pp"  => 1.0,
            "0 1 pp"  => -1.0,
            "2 1 pbb" => 2.0,
            "2 0"     => 1.3333333333333333,
            "1 0 b"   => 1.0,
            "0 2 pbp" => -1.0,
            "2 0 pbb" => 2.0,
            "0 1 pbp" => -1.0,
            "0 1 b"   => 0.0,
            "2 1 b"   => 1.3333333333333333,
            "2 1 pbp" => -1.0,
            "2"       => 1.1666666666666665,
            "1"       => -0.33333333333333337,
            "0"       => -1.0,
            "0 1 p"   => -1.0,
            "1 0 p"   => 0.3333333333333333,
            "0 2 pbb" => -2.0,
            "0 1 pbb" => -2.0,
            "2 1 p"   => 1.0,
            "0 2 pb"  => -1.0)

        history_and_probs_1 = Dict(
            ""        => 0.0555555555556,
            "0"       => 0.9,
            "0 1"     => 0.6,
            "0 1 b"   => -1.0,
            "0 1 bb"  => 2.0,
            "0 1 bp"  => -1.0,
            "0 1 p"   => 1.0,
            "0 1 pb"  => 1.0,
            "0 1 pbb" => 2.0,
            "0 1 pbp" => 1.0,
            "0 1 pp"  => 1.0,
            "0 2"     => 1.2,
            "0 2 b"   => 2.0,
            "0 2 bb"  => 2.0,
            "0 2 bp"  => -1.0,
            "0 2 p"   => 1.0,
            "0 2 pb"  => 1.0,
            "0 2 pbb" => 2.0,
            "0 2 pbp" => 1.0,
            "0 2 pp"  => 1.0,
            "1"       => 0.266666666667,
            "1 0"     => -1.0,
            "1 0 b"   => -1.0,
            "1 0 bb"  => -2.0,
            "1 0 bp"  => -1.0,
            "1 0 p"   => -1.0,
            "1 0 pb"  => -0.6,
            "1 0 pbb" => -2.0,
            "1 0 pbp" => 1.0,
            "1 0 pp"  => -1.0,
            "1 2"     => 1.53333333333,
            "1 2 b"   => 2.0,
            "1 2 bb"  => 2.0,
            "1 2 bp"  => -1.0,
            "1 2 p"   => 1.53333333333,
            "1 2 pb"  => 1.53333333333,
            "1 2 pbb" => 2.0,
            "1 2 pbp" => 1.0,
            "1 2 pp"  => 1.0,
            "2"       => -1.0,
            "2 0"     => -1.0,
            "2 0 b"   => -1.0,
            "2 0 bb"  => -2.0,
            "2 0 bp"  => -1.0,
            "2 0 p"   => -1.0,
            "2 0 pb"  => -2.0,
            "2 0 pbb" => -2.0,
            "2 0 pbp" => 1.0,
            "2 0 pp"  => -1.0,
            "2 1"     => -1.0,
            "2 1 b"   => -1.0,
            "2 1 bb"  => -2.0,
            "2 1 bp"  => -1.0,
            "2 1 p"   => -1.0,
            "2 1 pb"  => -2.0,
            "2 1 pbb" => -2.0,
            "2 1 pbp" => 1.0,
            "2 1 pp"  => -1.0)

        bp0 = BestResponsePolicy(get_optimal_kuhn_policy(), env, 0)
        bp1 = BestResponsePolicy(get_optimal_kuhn_policy(), env, 1)

        walk(env) do x
            @test RLZoo.best_response_value(bp0, x) ≈ history_and_probs_0[string(x.state)]
            @test RLZoo.best_response_value(bp1, x) ≈ history_and_probs_1[string(x.state)]
        end
    end
end


text:
# Run tests on functions in src/by_reference_doc/episodes.jl

using Spotify.Episodes

@testset verbose = true "GET-request endpoints for episodes" begin

    # Input argument from composite types defined in src/types.jl
    episode_id = SpEpisodeId()

    # Cycle through different input keywords for testing
    markets = ["US", "NL", "DE"]
    offsets = [7, 14]

    @testset "For market = $(market_id)" for market_id in markets

        @test ~isempty(episodes_get_single(episode_id, market = market_id)[1])
        @test ~isempty(episodes_get_multiple("$(episode_id),0Q86acNRm6V9GYx55SXKwf",
                                              market = market_id)[1])
        @test ~isempty(episodes_get_saved(market = market_id)[1])

    end

    @testset "For offset = $(num_offset)" for num_offset in offsets

        @test ~isempty(episodes_get_saved(offset = num_offset)[1])

    end

    @test ~isempty(episodes_get_contains("$(episode_id),0Q86acNRm6V9GYx55SXKwf")[1])

end

ext: jl


text:
using QuantumOptics: projector, tensor, SparseOperator, DenseOperator, basisstate, Ket
using LinearAlgebra: diagm
import QuantumOptics: displace, thermalstate, coherentthermalstate, fockstate

export create, destroy, number, displace, coherentstate, coherentthermalstate, fockstate,
       thermalstate, sigma, ionprojector, ionstate

#############################################################################################
# VibrationalMode operators
#############################################################################################

"""
    create(v::VibrationalMode)
returns the creation operator for `v` such that: `create(v) * v[i] = √(i+1) * v[i+1]`.
"""
create(v::VibrationalMode) = SparseOperator(v, diagm(-1 => sqrt.(1:v.N)))

"""
    destroy(v::VibrationalMode)
Returns the destruction operator for `v` such that: `destroy(v) * v[i] = √i * v[i-1]`.
"""
destroy(v::VibrationalMode) = create(v)'

"""
    number(v::VibrationalMode)
Returns the number operator for `v` such that:  `number(v) * v[i] = i * v[i]`.
"""
number(v::VibrationalMode) = SparseOperator(v, diagm(0 => 0:v.N))

"""
    displace(v::VibrationalMode, α::Number; method="truncated")
Returns the displacement operator ``D(α)`` corresponding to `v`.

If `method="truncated"` (default), the matrix elements are computed according to
``D(α) = exp[αa^† - α^*a]`` where ``a`` and ``a^†`` live in a truncated Hilbert space of
dimension `v.N+1`.
Otherwise if `method="analytic"`, the matrix elements are computed assuming an
infinite-dimension Hilbert space. In general, this option will not return a unitary operator.
"""
function displace(v::VibrationalMode, α::Number; method="truncated")
    @assert v.N ≥ abs(α) "`α` must be less than `v.N`"
    @assert method in ["truncated", "analytic"] "method ∉ [truncated, analytic]"
    D = zeros(ComplexF64, v.N+1, v.N+1)
    if α == 0
        return one(v)
    elseif method ≡ "analytic"
        @inbounds begin
            @simd for n in 1:v.N+1
                @simd for m in 1:v.N+1
                    D[n, m] = _Dnm(α, n, m)
                end
            end
        end
        return DenseOperator(v, D)
    elseif method ≡ "truncated"
        return exp(dense(α * create(v) - conj(α) * destroy(v)))
    end
end

"""
    thermalstate(v::VibrationalMode, n̄::Real; method="truncated")
Returns a thermal density matrix with ``⟨a^†a⟩ ≈ n̄``. Note: approximate because we are
dealing with a finite dimensional Hilbert space that must be normalized.

`method` can be set to either `"truncated"` (default) or `"analytic"`. In the former case,
the thermal density matrix is generated according to the formula:
``ρ_{th} = exp(-νa^†a/T) / Tr [exp(-νa^†a/T)]``. In the later case, the analytic formula,
assuming an infinite-dimensional Hilbert space, is used:
``[ρ_{th}]_{ij} = δ_{ij} \\frac{nⁱ}{(n+1)^{i+1}}.``
"""
function thermalstate(v::VibrationalMode, n̄::Real; method="truncated")
    @assert v.N ≥ n̄ "`n̄` must be less than `v.N`"
    @assert method in ["truncated", "analytic"] "method ∉ [truncated, analytic]"
    if n̄ == 0
        return v[0] ⊗ v[0]'
    elseif method ≡ "truncated"
        d = [(n̄ / (n̄+1))^i for i in 0:v.N]
        return DenseOperator(v, diagm(0 => d) ./ sum(d))
    elseif method ≡ "analytic"
        return DenseOperator(v, diagm(0 => [(n̄ / (n̄+1))^i / (n̄+1) for i in 0:v.N]))
    end
end

"""
    coherentstate(v::VibrationalMode, α::Number)
Returns a coherent state on `v` with complex amplitude ``α``.
"""
function coherentstate(v::VibrationalMode, α::Number)
    # this implementation is the same as in QuantumOptics.jl, but there the function is
    # restricted to v::FockBasis, so we must reimplement here
    @assert v.N ≥ abs(α) "`α` must be less than `v.N`"
    k = zeros(ComplexF64, v.N+1)
    k[1] = exp(-abs2(α) / 2)
    @inbounds for n=1:v.N
        k[n+1] = k[n] * α / √n
    end
    Ket(v, k)
end

"""
    coherentthermalstate(v::VibrationalMode, n̄::Real, α::Number; method="truncated)
Returns a displaced thermal state for `v`, which is created by applying a displacement
operation to a thermal state. The mean occupation of the thermal state is `n̄` and `α` is the
complex amplitude of the displacement.

`method` can be either `"truncated"` or `"analytic"` and this argument determines how the
displacement operator is computed (see: [`displace`](@ref)) .
"""
function coherentthermalstate(v::VibrationalMode, n̄::Real, α::Number; method="truncated")
    @assert (v.N ≥ n̄ && v.N ≥ abs(α)) "`n̄`, `α` must be less than `v.N`"
    @assert method in ["truncated", "analytic"] "method ∉ [truncated, analytic]"
    if method ≡ "truncated"
        d = displace(v, α)
    elseif method ≡ "analytic"
        d = displace(v, α, method="analytic")
    end
    d * thermalstate(v, n̄) * d'
end

"""
    fockstate(v::VibrationalMode, N::Int)
Returns the fockstate ``|N⟩`` on `v`.
"""
fockstate(v::VibrationalMode, N::Int) = v[N]

#############################################################################################
# Ion operators
#############################################################################################

"""
    ionstate(object, index)
For `object<:Ion` and `index<:String`, this returns the ket corresponding to the `Ion` being
in the state ``|index⟩``. The object can also be an `IonConfiguration` or `Trap` instance, in
which case ``N`` arguments should be given in place of `index`, where ``N`` equals the number
of ions in the `IonConfiguration` or `Trap`. This will return the state
``|index₁⟩⊗|index₂⟩⊗...⊗|index\\_N⟩``.

Rather than `index<:String`, one may also specify `index<:Int`. If `object<:Ion`, this will
return the ket given by ``|index⟩ = (0 ... 1 ... 0)ᵀ`` where the nonzero element in the
column vector is located at `index`.
"""
function ionstate(I::Ion, state::String)
    s = I.selected_level_structure.keys
    @assert state in s "index not in selected_level_structure: $s"
    i = findall(s .≡ state)[1]
    basisstate(I, i)
end

function ionstate(IC::IonConfiguration, states::Union{String,Int}...)
    ions = IC.ions
    L = length(ions)
    @assert L ≡ length(states) "wrong number of states"
    tensor([ionstate(ions[i], states[i]) for i in 1:L]...)
end

ionstate(T::Trap, states::Union{String,Int}...) = ionstate(T.configuration, states...)
ionstate(I::Ion, level::Int) = basisstate(I, level)

"""
    sigma(ion::Ion, ψ1::Union{String,Int}[, ψ2::Union{String,Int}])
Returns ``|ψ1\\rangle\\langle ψ2|``, where ``|ψ_i\\rangle`` corresponds to the state
returned by `ion[ψᵢ]`.

If ψ2 is not given, then ``|ψ1\\rangle\\langle ψ1|`` is returned.
"""
sigma(ion::Ion, ψ1::T, ψ2::T) where {T<:Union{String,Int}} = sparse(projector(ion[ψ1], dagger(ion[ψ2])))
sigma(ion::Ion, ψ1::Union{String,Int}) = sigma(ion, ψ1, ψ1)

"""
    ionprojector(obj, states::Union{String,Int}...; only_ions=false)

If `obj<:IonConfiguration` this will return ``|ψ₁⟩⟨ψ₁|⊗...⊗|ψ\\_N⟩⟨ψ\\_N|⊗𝟙``
where ``|ψᵢ⟩`` = `obj.ions[i][states[i]]` and the identity operator ``𝟙`` is over all of the
COM modes considered in `obj`.

If `only_ions=true`, then the projector is defined only over the ion subspace.

If instead `obj<:Trap`, then this is the same as `obj = Trap.configuration`.
"""
function ionprojector(IC::IonConfiguration, states::Union{String,Int}...; only_ions=false)
    ions = IC.ions
    L = length(ions)
    @assert L ≡ length(states) "wrong number of states"
    modes = get_vibrational_modes(IC)
    observable = tensor([projector(ions[i][states[i]]) for i in 1:L]...)
    if !only_ions
        for mode in modes
            observable = observable ⊗ one(mode)
        end
    end
    observable
end

function ionprojector(T::Trap, states::Union{String,Int}...; only_ions=false)
    ionprojector(T.configuration, states..., only_ions=only_ions)
end

#############################################################################################
# internal functions
#############################################################################################

# computes iⁿ(-i)ᵐ * (s! / ((s+1) * √(m!n!)))
function _pf(s::Int, n::Int, m::Int)
    n -= 1; m -= 1; s -= 1;
    @assert n<=s && m<=s
    val = 1. / (s+1)
    for i in 0:s-2
        if (m-i > 0) && (n-i > 0)
            val *= (s-i) / (√((m-i) * (n-i)))
        elseif m-i > 0
            val *= (s-i) / (√(m-i))
        elseif n-i > 0
            val *= (s-i) / (√(n-i))
        else
            val *= (s-i)
        end
    end
    (-1im)^n * 1im^m * val
end

# computes the coefficients for the 'probabilist's' Hermite polynomial of order n
function _He(n::Int)
    a = zeros(Float64, n+2, n+2)
    a[1, 1] = 1; a[2, 1] = 0; a[2, 2] = 1
    for i in 2:n+1, j in 1:n+1
        if j ≡ 1
            a[i+1, j] = -(i-1) * a[i-1, j]
        else
            a[i+1, j] = a[i, j-1] - (i-1) * a[i - 1, j]
        end
    end
    [a[n+1, k+1] for k in 0:n]
end

# computes He_n(x) (nth order Hermite polynomial)
function _fHe(x::Real, n::Int)
    n -= 1
    He = 1.0, x
    if n < 2
        return He[n+1]
    end
    for i in 2:n
        He = He[2], x * He[2] - (i-1) * He[1]
    end
    He[2]
end

# computes the matrix elements ⟨m|Dˢ(α)|n⟩ for the truncated displacement operator Dˢ(α)
# which exists in a Hilbert space of dimension s
function _Dtrunc(Ω, Δ, η, ν, rs, s, n, prefactor, timescale, L, t)
    d = complex(1, 0)
    for i in 1:L
        val = 0.
        Δn = n[1][i] - n[2][i]
        for r in rs[i]
            val += exp(im * r * abs(η[i])) * _fHe(r, n[2][i]) * _fHe(r, n[1][i]) / _fHe(r, s[i])^2
        end
        d *= (exp(im * Δn * (2π * ν[i] * timescale * t + π/2 + π * (sign(η[i] < 0))))
              * val * prefactor[i])
    end
    g = Ω * exp(-1im * t * Δ)
    g * d, g * conj(d)
end

# associated Laguerre polynomial
function _alaguerre(x::Real, n::Int, k::Int)
    L = 1.0, -x + k + 1
    if n < 2
        return L[n+1]
    end
    for i in 2:n
        L = L[2], ((k + 2i - 1 - x) * L[2] - (k + i - 1) * L[1]) / i
    end
    L[2]
end

# matrix elements of the displacement operator in the Fock Basis, assuming an
# infinite-dimensional Hilbert space. https://doi.org/10.1103/PhysRev.177.1857
function _Dnm(ξ::Number, n::Int, m::Int)
    if n < m
        return (-1)^isodd(abs(n-m)) * conj(_Dnm(ξ, m, n))
    end
    n -= 1; m -= 1
    s = 1.0
    for i in m+1:n
        s *= i
    end
    ret = sqrt(1 / s) * ξ^(n-m) * exp(-abs2(ξ) / 2.0) * _alaguerre(abs2(ξ), m, n-m)
    if isnan(ret)
        return 1.0 * (n==m)
    end
    ret
end


text:
module_path = splitdir(@__FILE__)[1]

# example 5 - difference between two flow fields

using JLD
using LinearAlgebra
# Load WFVisual module
import WFVisual
using Statistics
wfv=WFVisual

# Load GeometricTools: https://github.com/byuflowlab/GeometricTools.jl
import GeometricTools
gt=GeometricTools
                                              # Case of Julia v0.7
data_path = joinpath(module_path, "../datav07/")

# Create save path where to store vtks
save_path = "example_5/"         # Save path of this example
gt.create_path(save_path, true)

# --------------------- WIND FARM LAYOUT ---------------------------------------
turbine_x = [  0.0, 1000.0, -1000.0 ]

turbine_y = [  0.0, 0.0, 0.0]
turbine_z = [ 0.,  0.,  0.]

hub_height = [ 100.,  100.,  100.]
rotor_diameter = [ 100.,  100.,  100.]
nBlades = [3, 3, 3]

nturbs = length(turbine_x)
hub_height = zeros(nturbs) .+ 100.0
rotor_diameter = zeros(nturbs) .+ 100.0
nBlades = zeros(Int8,nturbs) .+ 3

wind_direction = 270.0

yaw = [ 0.,  0.,  0.] .- wind_direction

# --------------------- PERIMETER AND FLUID DOMAIN -----------------------------
NDIVSx = 200              # Cells in the parametric x-direction
NDIVSy = 200            # Cells in the parametric y-direction
NDIVSz = 50              # Cells in the geometric z-direction

# Dummy perimeter
Rper = 1500.0
perimeter_points = Rper.*[ [cos(a), sin(a), 0] for a in range(0, 2*pi, 179)]

# --------------------- GENERATE FLOW FIELD ------------------------------------

_zmin = 0
_zmax = maximum(hub_height) + 1.25*maximum(rotor_diameter)/2
z_off = 0.0
spl_s=0.01
fdom = wfv.generate_perimetergrid(perimeter_points,
                                  NDIVSx, NDIVSy, NDIVSz;
                                  z_min=_zmin+z_off, z_max=_zmax+z_off,
                                  verify_spline=false,
                                  spl_s=spl_s,
                                  save_path=nothing
                                )
# Guassian wake function
magVinf = 12.0                          # (m/s) wind speed
Vinf_Oaxis = gt.rotation_matrix2(0, 0, wind_direction) # Free-stream coordinate system
println(Vinf_Oaxis)
invVinf_Oaxis = Vinf_Oaxis'             # Inverse transformation
refH = mean(hub_height)                 # (m) reference height for shear layer

function wake(X)

  k = 0.0325 # we can let this vary
  CT = 8.0/9.0
  loss1 = 0
  alpha = 0.12

  loss2 = 0
  alpha_jensen = 0.1

  # Calculate loss in freestream direction
  for i = 1:length(turbine_x)     # Iterate over turbines

      # Distance from this turbine in global coordinates
      DX = X - [turbine_x[i], turbine_y[i], turbine_z[i]+hub_height[i]]

      # Distance in free-stream coordinates
      dX = invVinf_Oaxis*DX
      dx = dX[1]
      dy = dX[2]
      dz = dX[3]
      if dx > 0.
          sigma = k*dx + rotor_diameter[i]/sqrt(8.0)
          this_loss = (1.0-sqrt(1.0-CT/(8*sigma^2/(rotor_diameter[i]^2))))*
                          exp(-1.0/2.0*(dy/sigma)^2)*exp(-1.0/2.0*(dz/sigma)^2)
          loss1 += this_loss^2
      end

      r = rotor_diameter[i]/2.0 + alpha_jensen*dx
      if dx > 0. && abs(dy) <= r && abs(dz) <= r
        this_loss = 2.0/3.0 * (rotor_diameter[i]/2.0 / r)^2
        loss2 += this_loss^2
      end
  end
  loss1 = sqrt(loss1)
  loss2 = sqrt(loss2)

  windspeed= magVinf
  # Uncomment this line adds the shear layer
  # windspeed *= (X[3]/refH)^alpha

  # Return velocity vector in global coordinate system
  # return Vinf_Oaxis*[ windspeed*(1-loss), 0, 0 ]
  speed1 = windspeed*(1-loss1)
  speed2 = windspeed*(1-loss2)

  return Vinf_Oaxis*[ speed1 - speed2, 0, 0 ]
  # return Vinf_Oaxis*[ speed2, 0, 0 ]
end

gt.calculate_field(fdom, wake, "wake", "vector", "node")

# --------------------- GENERATE WIND FARM -------------------------------------
wfv.generate_windfarm(rotor_diameter, hub_height, nBlades,
                              turbine_x, turbine_y, turbine_z,
                              yaw,
                              perimeter_points, fdom;
                              NDIVSx=NDIVSx, NDIVSy=NDIVSy, NDIVSz=NDIVSz,
                              save_path=save_path, spl_s=0.01,
                              data_path=data_path, paraview=false);
nothing


text:
using ArgParse
using CUDAapi
using Dates
using LinearAlgebra
using Logging
using MPI
using Printf
using CLIMAParameters

using ..Atmos
using ..Callbacks
using ..ColumnwiseLUSolver
using ..ConfigTypes
using ..Diagnostics
using ..DGmethods
using ..DGmethods:
    vars_state_conservative, vars_state_auxiliary, update_auxiliary_state!
using ..DGmethods.NumericalFluxes
using ..HydrostaticBoussinesq
using ..Mesh.Grids
using ..Mesh.Topologies
using ..MoistThermodynamics
using ..MPIStateArrays
using ..ODESolvers
using ..TicToc
using ..VariableTemplates

using CuArrays, CuArrays.CUDAdrv, CuArrays.CUDAnative

function _init_array(::Type{CuArray})
    comm = MPI.COMM_WORLD
    # allocate GPUs among MPI ranks
    local_comm =
        MPI.Comm_split_type(comm, MPI.MPI_COMM_TYPE_SHARED, MPI.Comm_rank(comm))
    # we intentionally oversubscribe GPUs for testing: may want to disable this for production
    CUDAnative.device!(MPI.Comm_rank(local_comm) % length(devices()))
    CuArrays.allowscalar(false)
    return nothing
end

_init_array(::Type{Array}) = nothing

const cuarray_pkgid =
    Base.PkgId(Base.UUID("3a865a2d-5b23-5a0f-bc46-62713ec82fae"), "CuArrays")
function gpu_allowscalar(val)
    if haskey(Base.loaded_modules, CLIMA.cuarray_pkgid)
        Base.loaded_modules[CLIMA.cuarray_pkgid].allowscalar(val)
    end
    return nothing
end

# Note that the initial values specified here are overwritten by the
# command line argument defaults in `parse_commandline()`.
Base.@kwdef mutable struct CLIMA_Settings
    disable_gpu::Bool = false
    show_updates::String = "60secs"
    diagnostics::String = "never"
    vtk::String = "never"
    monitor_timestep_duration::String = "never"
    monitor_courant_numbers::String = "never"
    checkpoint::String = "never"
    checkpoint_keep_one::Bool = true
    checkpoint_at_end::Bool = false
    checkpoint_dir::String = "checkpoint"
    restart_from_num::Int = -1
    log_level::String = "INFO"
    output_dir::String = "output"
    integration_testing::Bool = false
    array_type
end
const Settings = CLIMA_Settings(array_type = Array)

"""
    parse_commandline()
"""
function parse_commandline(custom_settings)
    exc_handler =
        isinteractive() ? ArgParse.debug_handler : ArgParse.default_handler
    s = ArgParseSettings(
        prog = PROGRAM_FILE,
        description = "Climate Machine: an Earth System Model that automatically learns from data\n",
        preformatted_description = true,
        epilog = """
            Any <interval> unless otherwise stated may be specified as:
                - 2hours or 10mins or 30secs => wall-clock time
                - 9.5smonths or 3.3sdays or 1.5shours => simulation time
                - 1000steps => simulation steps
                - never => disable
                - default => use experiment specified interval (only for diagnostics at present)
            """,
        preformatted_epilog = true,
        version = string(CLIMA_VERSION),
        exc_handler = exc_handler,
    )
    add_arg_group!(s, "CLIMA")

    @add_arg_table! s begin
        "--disable-gpu"
        help = "do not use the GPU"
        action = :store_true
        "--show-updates"
        help = "interval at which to show simulation updates"
        metavar = "<interval>"
        arg_type = String
        default = "60secs"
        "--diagnostics"
        help = "interval at which to collect diagnostics"
        metavar = "<interval>"
        arg_type = String
        default = "never"
        "--vtk"
        help = "interval at which to output VTK"
        metavar = "<interval>"
        arg_type = String
        default = "never"
        "--monitor-timestep-duration"
        help = "interval in time-steps at which to output wall-clock time per time-step"
        metavar = "<interval>"
        arg_type = String
        default = "never"
        "--monitor-courant-numbers"
        help = "interval at which to output acoustic, advective, and diffusive Courant numbers"
        metavar = "<interval>"
        arg_type = String
        default = "never"
        "--checkpoint"
        help = "interval at which to create a checkpoint"
        metavar = "<interval>"
        arg_type = String
        default = "never"
        "--checkpoint-keep-all"
        help = "keep all checkpoints (instead of just the most recent)"
        action = :store_true
        "--checkpoint-at-end"
        help = "create a checkpoint at the end of the simulation"
        action = :store_true
        "--checkpoint-dir"
        help = "the directory in which to store checkpoints"
        metavar = "<path>"
        arg_type = String
        default = "checkpoint"
        "--restart-from-num"
        help = "checkpoint number from which to restart (in <checkpoint-dir>)"
        metavar = "<number>"
        arg_type = Int
        default = -1
        "--log-level"
        help = "set the log level to one of debug/info/warn/error"
        metavar = "<level>"
        arg_type = String
        default = "info"
        "--output-dir"
        help = "directory for output data"
        metavar = "<path>"
        arg_type = String
        default = "output"
        "--integration-testing"
        help = "enable integration testing"
        action = :store_true
    end

    if custom_settings !== nothing
        import_settings!(s, custom_settings)
    end

    return parse_args(s)
end

"""
    CLIMA.array_type()

Return the array type used by CLIMA. This defaults to (CPU-based) `Array`
and is only correctly set (based on choice from the command line, from
an environment variable, or from experiment code) after `CLIMA.init()`
is called.
"""
array_type() = Settings.array_type

"""
    CLIMA.init(
        ;
        disable_gpu = false,
        arg_settings = nothing,
    )

Perform necessary initializations for CLIMA:
- Initialize MPI.
- Parse command line arguments. To support experiment-specific arguments,
`arg_settings` may be specified (it is an `ArgParse.ArgParseSettings`);
it will be imported into CLIMA's settings.
- Determine whether GPU(s) is available and should be used (pass
`disable-gpu = true` if not) and set the CLIMA array type appropriately.
- Set up the global logger.

Returns a `Dict` containing non-CLIMA command-line arguments.
"""
function init(; disable_gpu = false, arg_settings = nothing)
    # set up timing mechanism
    tictoc()

    # parse command line arguments
    parsed_args = nothing
    try
        parsed_args = parse_commandline(arg_settings)
        Settings.disable_gpu = disable_gpu || parsed_args["disable-gpu"]
        delete!(parsed_args, "disable-gpu")
        Settings.show_updates = parsed_args["show-updates"]
        delete!(parsed_args, "show-updates")
        Settings.diagnostics = parsed_args["diagnostics"]
        delete!(parsed_args, "diagnostics")
        Settings.vtk = parsed_args["vtk"]
        delete!(parsed_args, "vtk")
        Settings.monitor_timestep_duration =
            parsed_args["monitor-timestep-duration"]
        delete!(parsed_args, "monitor-timestep-duration")
        Settings.monitor_courant_numbers =
            parsed_args["monitor-courant-numbers"]
        delete!(parsed_args, "monitor-courant-numbers")
        Settings.log_level = uppercase(parsed_args["log-level"])
        delete!(parsed_args, "log-level")
        Settings.checkpoint = parsed_args["checkpoint"]
        delete!(parsed_args, "checkpoint")
        Settings.checkpoint_keep_one = !parsed_args["checkpoint-keep-all"]
        delete!(parsed_args, "checkpoint-keep-all")
        Settings.checkpoint_at_end = parsed_args["checkpoint-at-end"]
        delete!(parsed_args, "checkpoint-at-end")
        Settings.checkpoint_dir = parsed_args["checkpoint-dir"]
        delete!(parsed_args, "checkpoint-dir")
        Settings.restart_from_num = parsed_args["restart-from-num"]
        delete!(parsed_args, "restart-from-num")
        Settings.output_dir = parsed_args["output-dir"]
        delete!(parsed_args, "output-dir")
        Settings.integration_testing = parsed_args["integration-testing"]
        delete!(parsed_args, "integration-testing")
    catch
        Settings.disable_gpu = disable_gpu
    end

    # initialize MPI
    if !MPI.Initialized()
        MPI.Init()
    end

    # set up the array type appropriately depending on whether we're using GPUs
    if get(ENV, "CLIMA_GPU", "") != "false" &&
       !Settings.disable_gpu &&
       CUDAapi.has_cuda_gpu()
        atyp = CuArrays.CuArray
    else
        atyp = Array
    end
    _init_array(atyp)
    Settings.array_type = atyp

    # create the output directory if needed
    if Settings.diagnostics !== "never" || Settings.vtk !== "never"
        mkpath(Settings.output_dir)
    end
    if Settings.checkpoint !== "never" || Settings.checkpoint_at_end
        mkpath(Settings.checkpoint_dir)
    end

    # set up logging
    loglevel = Settings.log_level == "DEBUG" ? Logging.Debug :
        Settings.log_level == "WARN" ? Logging.Warn :
        Settings.log_level == "ERROR" ? Logging.Error : Logging.Info
    # TODO: write a better MPI logging back-end and also integrate Dlog for large scale
    logger_stream = MPI.Comm_rank(MPI.COMM_WORLD) == 0 ? stderr : devnull
    global_logger(ConsoleLogger(logger_stream, loglevel))

    return parsed_args
end

include("driver_configs.jl")
include("solver_configs.jl")
include("diagnostics_configs.jl")

"""
    CLIMA.invoke!(
        solver_config::SolverConfiguration;
        diagnostics_config = nothing,
        user_callbacks = (),
        check_euclidean_distance = false,
        adjustfinalstep = false,
        user_info_callback = (init) -> nothing,
    )

Run the simulation defined by `solver_config`.

Keyword Arguments:

The `user_callbacks` are passed to the ODE solver as callback functions;
see [`ODESolvers.solve!]@ref().

If `check_euclidean_distance` is `true, then the Euclidean distance
between the final solution and initial condition function evaluated with
`solver_config.timeend` is reported.

The value of 'adjustfinalstep` is passed to the ODE solver; see
[`ODESolvers.solve!]@ref().

The function `user_info_callback` is called after the default info
callback (which is called every `Settings.show_updates` interval). The
single input argument `init` is `true` when the callback is called for
initialization (before time stepping begins) and `false` when called
during the actual ODE solve; see [`GenericCallbacks`](@ref) and
[`ODESolvers.solve!]@ref().
"""
function invoke!(
    solver_config::SolverConfiguration;
    diagnostics_config = nothing,
    user_callbacks = (),
    check_euclidean_distance = false,
    adjustfinalstep = false,
    user_info_callback = (init) -> nothing,
)
    mpicomm = solver_config.mpicomm
    dg = solver_config.dg
    bl = dg.balancelaw
    Q = solver_config.Q
    FT = eltype(Q)
    timeend = solver_config.timeend
    init_on_cpu = solver_config.init_on_cpu
    init_args = solver_config.init_args
    solver = solver_config.solver

    # set up callbacks
    callbacks = ()

    # info callback
    cb_updates = Callbacks.show_updates(
        Settings.show_updates,
        solver_config,
        user_info_callback,
    )
    if cb_updates !== nothing
        callbacks = (callbacks..., cb_updates)
    end

    # diagnostics callback(s)
    if Settings.diagnostics !== "never" && diagnostics_config !== nothing
        dgn_starttime = replace(string(now()), ":" => ".")
        Diagnostics.init(mpicomm, dg, Q, dgn_starttime, Settings.output_dir)

        dgncbs = Callbacks.diagnostics(
            Settings.diagnostics,
            solver_config,
            dgn_starttime,
            diagnostics_config,
        )
        callbacks = (callbacks..., dgncbs...)
    end

    # vtk callback
    cb_vtk = Callbacks.vtk(Settings.vtk, solver_config, Settings.output_dir)
    if cb_vtk !== nothing
        callbacks = (callbacks..., cb_vtk)
    end

    # timestep duration monitor
    cb_mtd = Callbacks.monitor_timestep_duration(
        Settings.monitor_timestep_duration,
        Settings.array_type,
        mpicomm,
    )
    if cb_mtd !== nothing
        callbacks = (callbacks..., cb_mtd)
    end

    # Courant number monitor
    cb_mcn = Callbacks.monitor_courant_numbers(
        Settings.monitor_courant_numbers,
        solver_config,
    )
    if cb_mcn !== nothing
        callbacks = (callbacks..., cb_mcn)
    end

    # checkpointing callback
    cb_checkpoint = Callbacks.checkpoint(
        Settings.checkpoint,
        Settings.checkpoint_keep_one,
        solver_config,
        Settings.checkpoint_dir,
    )
    if cb_checkpoint !== nothing
        callbacks = (callbacks..., cb_checkpoint)
    end

    # user callbacks
    callbacks = (user_callbacks..., callbacks...)

    # initial condition norm
    eng0 = norm(Q)
    @info @sprintf(
        """
Starting %s
    dt              = %.5e
    timeend         = %8.2f
    number of steps = %d
    norm(Q)         = %.16e""",
        solver_config.name,
        solver_config.dt,
        solver_config.timeend,
        solver_config.numberofsteps,
        eng0
    )

    # run the simulation
    @tic solve!
    solve!(
        Q,
        solver;
        timeend = timeend,
        callbacks = callbacks,
        adjustfinalstep = adjustfinalstep,
    )
    @toc solve!

    # write end checkpoint if requested
    if Settings.checkpoint_at_end
        Callbacks.write_checkpoint(
            solver_config,
            Settings.checkpoint_dir,
            solver_config.name,
            mpicomm,
            solver_config.numberofsteps,
        )
    end

    # fini diagnostics groups
    if Settings.diagnostics !== "never" && diagnostics_config !== nothing
        currtime = ODESolvers.gettime(solver)
        for dgngrp in diagnostics_config.groups
            dgngrp(currtime, fini = true)
        end
    end

    engf = norm(Q)
    @info @sprintf(
        """
Finished
    norm(Q)            = %.16e
    norm(Q) / norm(Q₀) = %.16e
    norm(Q) - norm(Q₀) = %.16e""",
        engf,
        engf / eng0,
        engf - eng0
    )

    if check_euclidean_distance
        Qe =
            init_ode_state(dg, timeend, init_args...; init_on_cpu = init_on_cpu)
        engfe = norm(Qe)
        errf = euclidean_distance(Q, Qe)
        @info @sprintf(
            """
Euclidean distance
    norm(Q - Qe)            = %.16e
    norm(Q - Qe) / norm(Qe) = %.16e""",
            errf,
            errf / engfe
        )
    end

    return engf / eng0
end


text:
```
a(n) Abstract `seed` of idea
that has value

value
is a value fixed (or changing?)
what changes in a value

definitely:
    Any `thing` : requires time & effort

    effort: could be the production material the `thing` is made of

how is a value related to an idea?
(Realizing) the relation is where the bottle neck is

is there something else (you could think of ?)
```
struct Value <: idea

    time::DateTime #estimated time - how much time spent (is it kuantifiable?, monitizable?)
    effort
end

max_stars_repo_name: adamwillisXanax/accountingJulia

max_issues_repo_name: adamwillisXanax/accountingJulia

max_forks_repo_name: adamwillisXanax/accountingJulia


text:
function isapprox_zero(x::Symbolics.Num)
    return  x === zero(x)
end

# Cheat and don't convert to a Num
function Base.convert(::Type{Symbolics.Num}, z::Z4Group0)
    if isreal(z)
        return convert(Int, z)
    else
        return convert(Complex{Int}, z)
    end
end


text:
using Revise
module LeverageChanges
using Revise

using  CSV, DataFrames, Serialization, Dates, Statistics,
  Distributed, Finometrics, StatsBase,
  CodecLz4, CodecZlib,  Distributions, LinearAlgebra,
  CuArrays,#Plots,
  #Gadfly, Cairo, Fontconfig, Compose
  VegaLite

if  pwd() ∉ LOAD_PATH
  push!(LOAD_PATH,pwd())
end

import Base: eltype, getindex, iterate, length, isfinite, merge

#necessary due to exporting conflict
describe(df::DataFrame) = DataFrames.describe(df)

#############################CONSTANTS#####################
TEST_OUTPUT = true

const OUT_PATH = pwd() * "\\output"
const DATA_PATH = pwd() * "\\data"
const ARCHIVE_PATH = pwd() * "\\working"
const TEX_PATH = "C:\\Users\\Clinton\\Dropbox\\Apps\\Overleaf\\LeverageChangesRep"
const TABLE_PATH = "$TEX_PATH\\sub\\tab"
const FIGURE_PATH = "$TEX_PATH\\sub\\fig"
const DEFAULT_DECIMALS = 2

const CRSP_TYPE = :monthly
const CRSP_DATE_FORMAT = dateformat"yyyymmdd"

const CRSPM_NAME = "CRSP-M"
const CRSPM_PARTS = ["CRSP-M"]

const CRSPD_NAME = "CRSP-D"
const CRSPD_PARTS = ["CRSP-D-II", "CRSP-D-I"]

const COMP_NAME = "COMP"
const COMP_A_NAME = "COMP-A"
const COMP_Q_NAME = "COMP-Q"
const COMP_DATE_FORMAT = dateformat"yyyymmdd"

const CCM_NAME = "CCM"
const CCM_DATE_FORMAT = dateformat"yyyymmdd"
const LAST_DATE_STRING = "21001231"
const LAST_DATE = Date(LAST_DATE_STRING, CCM_DATE_FORMAT)

const MONTHS_2_STALE_LONG = 18
const MONTHS_2_STALE_SHORT = 6 #maybe mess with this?
const WINSORIZE_PROP = 0.001

const UNIV_NAME = "UNIV"

const LEV_NAME = "LEV"

const PORT_NAME = "PORT"
const PANEL_NAME = "PANEL"

const VALIDATE_MERGED = true

const USE_QUARTERLY_DATES = true

##constants related to outside data
const VWCRSP_NAME = "VWCRSP-D"
const RFCRSP_NAME = "RFCRSP-M"
const RETAINED_VWCRSP_COLS = [:date, :vwindd, :ewindd]
const RETAINED_RFCRSP_COLS = [:date, :t30ind]
const IDXCRSP_PATH = "$DATA_PATH\\crspidx"
const IDXCRSP_DATE_FORMAT = CRSP_DATE_FORMAT
const IDX_DATE_COL = :caldt

const FF5_NAME = "FF5"
const FF3M_NAME = "FF3M"
const FF_PATH = "$DATA_PATH\\french"
const FF_DATE_FORMAT = CRSP_DATE_FORMAT

const EVENT_NAME = "event"
const EVENT_PATH = "$DATA_PATH\\event"
const CRSP_DISTRIBUTIONS_NAME  = "crspdistributions"

#############################Multi-threading

#available lock if needed- use lock(spinlock) unlock(spinlock) etc
const spinlock = Threads.SpinLock()

#####################REFERENCE CONSTANTS
const PARALLEL = Ref{Bool}(true)
const REPLICATION_TYPE = Ref{Symbol}(:daily)
const QRTYPE = Ref{Type}(Matrix{Float64})

#NOTE: This macro allows for the parallelization routines to be shut down
#by flicking the above PARALLEL switch
macro mpar(expr)
  quote
    if PARALLEL[]
      :($(Threads.@threads($expr)))
    else
      :($($expr))
    end
  end
end

###########################IO Compression/Decompression routines
IN_CSV_STREAM(p::String) = LZ4DecompressorStream(open(p))
IN_CSV_STREAM(F::Function, p::String) = open(F, LZ4DecompressorStream, p)

IN_JLS_STREAM(p::String) = LZ4DecompressorStream(open(p))
IN_JLS_STREAM(F::Function, p::String) = open(F, LZ4DecompressorStream, p)

OUT_JLS_STREAM(p::String) = LZ4CompressorStream(open(p, "w"))
OUT_JLS_STREAM(F::Function, p::String) = open(F, LZ4CompressorStream, p, "w")

############################FILES
include("PanelFormation\\PreprocessCRSP.jl")
include("PanelFormation\\PreprocessCOMP.jl")
include("PanelFormation\\MergeCOMPCCM.jl")
include("PanelFormation\\MergeCRSPCOMP.jl")
include("Utilities.jl")
include("PanelFormation\\FormDataSeries.jl")
include("PortfolioFormation\\PortfolioConstruction.jl")
include("PortfolioFormation\\ConstructPortfolios.jl")
include("PanelFormation\\MergeOthers.jl")

include("Analysis\\AnalyzeFactors.jl")
include("Analysis\\FactorTests.jl")
include("Tables\\Table1.jl")
include("Tables\\Table2.jl")
include("Tables\\Table3.jl")
include("Tables\\Table4.jl")
include("Tables\\Table5.jl")
include("Tables\\Table6.jl")
include("Tables\\Table7.jl")
include("Tables\\Table8.jl")
include("Figures\\Figure1.jl")
include("Analysis\\Replicate.jl")

include("Event\\SEO.jl")
include("Event\\Distributions.jl")
include("Event\\FormEvent.jl")

end # module

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
"""
Write model configurations to options file (OPS.h).
"""
function write_modcfg(modcfg)

    info("Writing model configs (opts.h)")

    str = """
    /* Process options */
    #define ALBEDO $(modcfg["albedo"])  /* snow albedo: 0 - diagnostic, 1 - prognostic                          */
    #define CANMOD $(modcfg["canmod"])  /* forest canopy: 0 - zero layer, 1 - one layer                         */
    #define CONDCT $(modcfg["condct"])  /* snow thermal conductivity: 0 - constant, 1 - Yen (1981)              */
    #define DENSTY $(modcfg["densty"])  /* snow density: 0 - constant, 1 - Verseghy (1991), 2 - Anderson (1976) */
    #define EXCHNG $(modcfg["exchng"])  /* turbulent exchange: 0 - constant, 1 - Louis (1979)                   */
    #define HYDROL $(modcfg["hydrol"])  /* snow hydraulics: 0 - free draining, 1 - bucket                       */

    /* Driving data options */
    #define DRIV1D 0   /* 1D driving data format: 0 - FSM, 1 - ESM-SnowMIP                  */
    #define DOWNSC 0   /* 1D driving data downscaling: 0 - no, 1 - yes                      */
    #define DEMHDR 0   /* DEM header: 0 - none, 1 - ESRI format                             */
    #define SWPART 0   /* SW radiation: 0 - total, 1 - direct and diffuse calculated        */
    #define ZOFFST 1   /* Measurement height offset: 0 - above ground, 1 - above canopy top */
    """

    # Write file

    fid = open("OPTS.h", "w")

    write(fid, str)

    close(fid)

    return nothing

end

"""
Compile the model for current configuration and netcdf version.
"""
function compile_netcdf(path, modcfg, icfg)

    curr = pwd()

    cd(joinpath(path, "src"))

    write_modcfg(modcfg)

    icfg = string(icfg)

    dest = joinpath(path, "bin", "FSM_$(icfg)")

    mods = ["DATANC.F90", "MODULES.F90"]

    routines = ["CANOPY.F90", "CUMULATE.F90", "DRIVENC.F90", "DUMP.F90", "EBALFOR.F90",
                "EBALSRF.F90", "FSMNC.F90", "LUDCMP.F90", "OUTPUTNC.F90", "PHYSICS.F90", "QSAT.F90",
                "RADIATION.F90", "READMAPS.F90", "SETUPNC.F90", "SNOW.F90", "SOIL.F90", "SFEXCH.F90",
                "THERMAL.F90", "TRIDIAG.F90"]

    run(`gfortran -o $dest -O3 $mods $routines -g -O2 -I/usr/include -L/usr/lib -lnetcdff -lnetcdf`)

    rm.(filter!(f -> endswith(f, "mod"), readdir()))

    cd(curr)

end

"""
Compile the model for current configuration and standard version.
"""
function compile_standard(path, modcfg, icfg)

    curr = pwd()

    cd(joinpath(path, "src"))

    write_modcfg(modcfg)

    icfg = string(icfg)

    dest = joinpath(path, "bin", "FSM_$(icfg)")

    mods = ["MODULES.F90"]

    routines = ["CANOPY.F90", "CUMULATE.F90", "DRIVE.F90", "DUMP.F90", "EBALFOR.F90",
                "EBALOPN.F90", "FSM2.F90", "LUDCMP.F90", "OUTPUT.F90", "PHYSICS.F90",
                "QSAT.F90", "READMAPS.F90", "SETUP.F90", "SFEXCH.F90", "SNOW.F90", "SOIL.F90",
                "SWRAD.F90", "THERMAL.F90", "TRIDIAG.F90"]

    run(`gfortran -o $dest -O3 $mods $routines -g -O2`)

    rm.(filter!(f -> endswith(f, "mod"), readdir()))

    cd(curr)

end

"""
Generate dataframe with configurations.
"""
function cfg_table()

    df = DataFrame(albedo = [],
                   canmod = [],
                   condct = [],
                   densty = [],
                   exchng = [],
                   hydrol = [])

    for alb = 0:1, can = 0, con = 0:1, den = 0:1, exh = 0:1, hyd = 0:1
        push!(df, [alb, can, con, den, exh, hyd])
    end

    df[:cfg] = 1:32

    return df

end

"""
Compile all configuratins in df_cfg.
"""
function compile_all(path, df_cfg, version)

    icfg = 1

    for row in eachrow(df_cfg)

        @info "Compiling configuration $(icfg)"

        modcfg = Dict("albedo" => row[:albedo],
                      "canmod" => row[:canmod],
                      "condct" => row[:condct],
                      "densty" => row[:densty],
                      "exchng" => row[:exchng],
                      "hydrol" => row[:hydrol])

        if version == "netcdf"
            compile_netcdf(path, modcfg, icfg)
        end

        if version == "standard"
            compile_standard(path, modcfg, icfg)
        end

        icfg += 1

    end

end


text:
function pcontrib{T<:Integer}(p::T, a::T)
    n = one(T)
    pcon = one(T)
    for i in 1:a
        n *= p
        pcon += n
    end
    return pcon
end

function divisorsum{T<:Integer}(n::T)
    dsum = one(T)
    for (p, a) in factor(n)
        dsum *= pcontrib(p, a)
    end
    dsum -= n
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
using ScientificColourMaps
using PyPlot
using Test

cmaps = [ "acton"  ,"bamako" ,"batlow" ,"berlin" ,"bilbao" ,"broc"   ,"brocO"  ,"buda"   ,"cork",
"corkO"  ,"davos"  ,"devon"  ,"grayC"  ,"hawaii" ,"imola"  ,"lajolla","lapaz"  ,"lisbon" ,"nuuk",
"oleron" ,"oslo"   ,"roma"   ,"romaO"  ,"tofino" ,"tokyo"  ,"turku"  ,"vik"    ,"vikO"]

@testset "ColorMaps exist?" begin
    @testset for cmap in cmaps
        @eval begin
            @test typeof($(Symbol(cmap))) == ColorMap
        end
    end
end


text:
module TestDistance
using ParallelKMeans: pairwise!, SingleThread, MultiThread
using Test

@testset "naive singlethread pairwise" begin
    X = [1.0 2.0; 3.0 5.0; 4.0 6.0]
    y = [1.0 2.0; ]
    r = Array{Float64, 2}(undef, 3, 1)

    pairwise!(r, X, y)
    @test all(r .≈ [0.0, 13.0, 25.0])
end

@testset "multithread pairwise" begin
    X = [1.0 2.0; 3.0 5.0; 4.0 6.0]
    y = [1.0 2.0; ]
    r = Array{Float64, 2}(undef, 3, 1)

    pairwise!(r, X, y, MultiThread())
    @test all(r .≈ [0.0, 13.0, 25.0])
end

end # module


text:
replace_dots(field) = searchindex(field, ".") > 0 ? replace(field, ".", "_") : field

# _allowsnulls(A::Array) = false
# _allowsnulls(A::AbstractArray) = true
_hasnulls(A::Array) = false
_hasnulls(A::AbstractArray) = anynull(A)


text:
# BigBed Reader
# =============

struct Reader <: BioCore.IO.AbstractReader
    stream::IO
    header::BBI.Header
    zooms::Vector{BBI.Zoom}
    summary::BBI.Summary
    # chrom name => (ID, length)
    chroms::Dict{String,Tuple{UInt32,Int}}
    # chrom ID => name
    chrom_names::Dict{UInt8,String}
    # data index
    index::BBI.RTree
end

function Base.eltype(::Type{Reader})
    return Record
end

function BioCore.IO.stream(reader::Reader)
    return reader.stream
end

"""
    BigBed.Reader(input::IO)

Create a reader for bigBed file format.

Note that `input` must be seekable.
"""
function Reader(input::IO)
    # read header
    header = read(input, BBI.Header)
    if header.magic != BBI.BED_MAGIC
        error("invalid BigBed magic")
    elseif header.version < 3
        error("not a supported version of BigBed")
    end
    # read zoom objects
    zoom_headers = Vector{BBI.ZoomHeader}(undef, header.zoom_levels)
    read!(input, zoom_headers)
    zooms = [BBI.Zoom(input, h, header.uncompress_buf_size) for h in zoom_headers]
    sort!(zooms, by=z->z.header.reduction_level)
    # read summary, B tree, and R tree
    seek(input, header.total_summary_offset)
    summary = read(input, BBI.Summary)
    chromindex = BBI.BTree(input, header.chromosome_tree_offset)
    chroms = Dict(name => (id, Int(len)) for (name, id, len) in BBI.chromlist(chromindex))
    chrom_names = Dict(id => name for (name, (id, len)) in chroms)
    index = BBI.RTree(input, header.full_index_offset)
    return Reader(input, header, zooms, summary, chroms, chrom_names, index)
end

"""
    chromlist(reader::BigBed.Reader)::Vector{Tuple{String,Int}}

Get the `(name, length)` pairs of chromosomes/contigs.
"""
function chromlist(reader::Reader)::Vector{Tuple{String,Int}}
    return sort!([(name, Int(len)) for (name, (id, len)) in reader.chroms], by=x->x[1])
end

const data_machine = (function ()
    cat = Automa.RegExp.cat
    rep = Automa.RegExp.rep
    opt = Automa.RegExp.opt

    record = let
        chromid = re"...."
        chromid.actions[:exit] = [:record_chromid]

        chromstart = re"...."
        chromstart.actions[:exit] = [:record_chromstart]

        chromend = re"...."
        chromend.actions[:exit] = [:record_chromend]

        name = re"[ -~]*"
        name.actions[:enter] = [:mark]
        name.actions[:exit]  = [:record_name]

        score = re"[0-9]+"
        score.actions[:enter] = [:mark]
        score.actions[:exit]  = [:record_score]

        strand = re"[+\-.?]"
        strand.actions[:enter] = [:record_strand]

        thickstart = re"[0-9]+"
        thickstart.actions[:enter] = [:mark]
        thickstart.actions[:exit]  = [:record_thickstart]

        thickend = re"[0-9]+"
        thickend.actions[:enter] = [:mark]
        thickend.actions[:exit]  = [:record_thickend]

        itemrgb = cat(re"[0-9]+", opt(cat(',', re"[0-9]+", ',', re"[0-9]+")))
        itemrgb.actions[:enter] = [:mark]
        itemrgb.actions[:exit]  = [:record_itemrgb]

        blockcount = re"[0-9]+"
        blockcount.actions[:enter] = [:mark]
        blockcount.actions[:exit]  = [:record_blockcount]

        # comma-separated values
        csv(x) = cat(rep(cat(x, ',')), opt(x))

        blocksizes = let
            blocksize = re"[0-9]+"
            blocksize.actions[:enter] = [:mark]
            blocksize.actions[:exit]  = [:record_blocksizes_blocksize]

            csv(blocksize)
        end
        blocksizes.actions[:exit] = [:record_blocksizes]

        blockstarts = let
            blockstart = re"[0-9]+"
            blockstart.actions[:enter] = [:mark]
            blockstart.actions[:exit]  = [:record_blockstarts_blockstart]

            csv(blockstart)
        end
        blockstarts.actions[:exit] = [:record_blockstarts]

        cat(
            chromid, chromstart, chromend,
            opt(cat(      name,
            opt(cat('\t', score,
            opt(cat('\t', strand,
            opt(cat('\t', thickstart,
            opt(cat('\t', thickend,
            opt(cat('\t', itemrgb,
            opt(cat('\t', blockcount,
            opt(cat('\t', blocksizes,
            opt(cat('\t', blockstarts)))))))))))))))))))
    end
    record.actions[:enter] = [:anchor]
    record.actions[:exit]  = [:record]

    data = rep(cat(record, '\0'))

    return Automa.compile(data)
end)()

const actions = Dict(
    :record_chromid                => :(record.chromid    = unsafe_load(convert(Ptr{UInt32}, pointer(data, p-4))); record.ncols += 1),
    :record_chromstart             => :(record.chromstart = unsafe_load(convert(Ptr{UInt32}, pointer(data, p-4))); record.ncols += 1),
    :record_chromend               => :(record.chromend   = unsafe_load(convert(Ptr{UInt32}, pointer(data, p-4))); record.ncols += 1),
    :record_name                   => :(record.name       = (mark:p-1) .- offset; record.ncols += 1),
    :record_score                  => :(record.score      = (mark:p-1) .- offset; record.ncols += 1),
    :record_strand                 => :(record.strand     =       p    .- offset; record.ncols += 1),
    :record_thickstart             => :(record.thickstart = (mark:p-1) .- offset; record.ncols += 1),
    :record_thickend               => :(record.thickend   = (mark:p-1) .- offset; record.ncols += 1),
    :record_itemrgb                => :(record.itemrgb    = (mark:p-1) .- offset; record.ncols += 1),
    :record_blockcount             => :(record.blockcount = (mark:p-1) .- offset; record.ncols += 1),
    :record_blocksizes_blocksize   => :(push!(record.blocksizes, (mark:p-1) .- offset)),
    :record_blocksizes             => :(record.ncols += 1),
    :record_blockstarts_blockstart => :(push!(record.blockstarts, (mark:p-1) .- offset)),
    :record_blockstarts            => :(record.ncols += 1),
    :record => quote
        BioCore.ReaderHelper.resize_and_copy!(record.data, data, BioCore.ReaderHelper.upanchor!(stream):p-1)
        record.filled = (offset+1:p-1) .- offset
        found_record = true
        @escape
    end,
    :countrecord => :(),
    :mark => :(mark = p),
    :anchor => :(BioCore.ReaderHelper.anchor!(stream, p); offset = p - 1))

mutable struct Record
    chromid::UInt32
    chromstart::UInt32
    chromend::UInt32
    # data and filled range
    data::Vector{UInt8}
    filled::UnitRange{Int}
    # number of columns
    ncols::Int
    # indexes
    name::UnitRange{Int}
    score::UnitRange{Int}
    strand::Int
    thickstart::UnitRange{Int}
    thickend::UnitRange{Int}
    itemrgb::UnitRange{Int}
    blockcount::UnitRange{Int}
    blocksizes::Vector{UnitRange{Int}}
    blockstarts::Vector{UnitRange{Int}}
    # reader
    reader::Reader

    function Record(chromid, chromstart, chromend,
                    data, filled, ncols,
                    name, score, strand, thickstart, thickend,
                    itemrgb, blockcount, blocksizes, blockstarts)
        return new(chromid, chromstart, chromend,
                   data, filled, ncols,
                   name, score, strand, thickstart, thickend,
                   itemrgb, blockcount, blocksizes, blockstarts)
    end
end

eval(
    BioCore.ReaderHelper.generate_read_function(
        Reader,
        data_machine,
        :(offset = mark = 0),
        actions))

# Iterator
# --------

mutable struct IteratorState
    state::BioCore.Ragel.State
    record::Record
    n_records::UInt64
    current_record::UInt64
end

function Base.iterate(reader::Reader)
    seek(reader.stream, reader.header.full_data_offset)
    # this is defined as UInt32 in the spces but actually UInt64
    record_count = read(reader.stream, UInt64)
    datastream = Libz.ZlibInflateInputStream(reader.stream)
    parser_state = BioCore.Ragel.State(data_machine.start_state, datastream)
    state = IteratorState(parser_state, Record(), record_count, 0)
    return iterate(reader, state)
end

function Base.iterate(reader::Reader, state::IteratorState)
    if state.current_record < state.n_records
        _read!(reader, state.state, state.record)
        state.record.reader = reader
        state.current_record += 1
        return copy(state.record), state
    end
    return nothing
end


text:
#This file is to randomly select 5% individuals (i.e., 1055) from 21100 individuals in Christensenet al.(2021).
# genotypes("genotypes.txt"), omics("M.txt"), phenotypes("y.txt"), breeding values("ybv.txt") data in Christensenet al.(2021) were used
# above data can be found in Christensenet al.(2021) or below link:
# http://genoweb.toulouse.inra.fr/~alegarra/GOBLUP/

using DataFrames,CSV,Random,DelimitedFiles,Distributions
cd("/Users/tianjing/Box/Omicsdata/andrew1055")

Random.seed!(123)

nData=20  #create 20 replicates
#selet 5% animals per generation
for rep in 1:nData
    mkdir("data$rep")
    println("rep$rep------------")
    pct=0.05
    selectIndiv=sample(1:1100,Int(1100*pct),replace=false, ordered=true)  #1100 individuals in generation1 in Christensenet al.(2021)
    println("generation:",1," start:",1," ","end:",1100," ","all:",1100, " selected:",length(selectIndiv))
    for i in 2:11 #11 generation in total
        start_i = 1100+(i-2)*2000+1  #2000 individuals in generation1 in Christensenet al.(2021)
        end_i = start_i+2000-1
        selectIndiv_i = sample(start_i:end_i,Int(2000*pct),replace=false, ordered=true)
        selectIndiv = [selectIndiv ; selectIndiv_i]
        println("generation:",i," start:",start_i," ","end:",end_i," ","all:",end_i-start_i+1, " selected:",length(selectIndiv_i))

    end
    @show length(selectIndiv)
    println(selectIndiv[1:5])
    writedlm("data$rep/select_1055_index_data$rep.txt",selectIndiv)
end

#save y (phenotype), omics
y=vec(readdlm("/Users/tianjing/Box/Omicsdata/andrewlegarra_simulated_omics/raw/y.txt"))  #phenotype data in Christensenet al.(2021)
M=readdlm("/Users/tianjing/Box/Omicsdata/andrewlegarra_simulated_omics/raw/M.txt")       #omics data in Christensenet al.(2021)
for rep in 1:nData #20 reps
    selectIndiv=Int.(vec(readdlm("data$rep/select_1055_index_data$rep.txt")))
    # save y
    println("number of individuals is ", length(y))
    y_select=y[selectIndiv]
    println("number of individuals after selection is ", length(y_select))
    writedlm("data$rep/y.txt",y_select)
    y_df=DataFrame(ID=selectIndiv,y=y_select)
    CSV.write("data$rep/y.csv", y_df)

    # save omics
    println("number of individuals is ", size(M,1), ", number of genes is ", size(M,2))
    M_select=M[selectIndiv,:]
    println("after selection, number of individuals is ", size(M_select,1), ", number of genes is ", size(M_select,2))
    writedlm("data$rep/omics.txt",M_select)
    names=["gene$i" for i in 1:size(M_select,2)]
    M_df=DataFrame(M_select,names)
    insertcols!(M_df, 1, :ID => selectIndiv)
    CSV.write("data$rep/omics.csv", M_df)
end

#remove ar from y
# ar is polygenetic effect whose covariance matrix is defined by the pedigree or/and genotypes in Christensenet al.(2021)
# this part is ignored here for simplicity
using DataFrames,CSV,Random,DelimitedFiles,Distributions
cd("/Users/tianjing/Box/Omicsdata/andrew1055")
tbv=readdlm("/Users/tianjing/Box/Omicsdata/andrewlegarra_simulated_omics/raw/tbv.txt")  #breeding values in Christensenet al.(2021)
ar=tbv[:,2]
y=vec(readdlm("/Users/tianjing/Box/Omicsdata/andrewlegarra_simulated_omics/raw/y.txt")) #phenotype data in Christensenet al.(2021)
y_no_ar= y-ar

for rep in 1:20
    println("rep$rep------------")
    selectIndiv=Int.(vec(readdlm("data$rep/select_1055_index_data$rep.txt")))
    # save y_no_ar
    y_no_ar_select = y_no_ar[selectIndiv]
    writedlm("data$rep/y_no_ar.txt",y_no_ar_select)
    y_no_ar_df=DataFrame(ID=selectIndiv, y=y_no_ar_select)
    CSV.write("data$rep/y_no_ar.csv", y_no_ar_df)
end

# save breeding values
tbv=readdlm("/Users/tianjing/Box/Omicsdata/andrewlegarra_simulated_omics/raw/tbv.txt") #breeding values in Christensenet al.(2021)
tbv_m=tbv[:,1]  #breeding values from omics (i.e., no ar)
for rep in 1:20
    println("rep$rep------------")
    selectIndiv=Int.(vec(readdlm("data$rep/select_1055_index_data$rep.txt")))
    bv_select = tbv_m[selectIndiv]
    writedlm("data$rep/tbv_m.txt",bv_select)
    bv_select_df=DataFrame(ID=selectIndiv, tbv_m=bv_select)
    CSV.write("data$rep/tbv_m.csv", bv_select_df)
end

#save genotype
using DataFrames,CSV,Random,DelimitedFiles,Distributions
cd("/Users/tianjing/Box/Omicsdata/andrew1055")

Random.seed!(123)

geno=readdlm("/Users/tianjing/Box/Omicsdata/andrewlegarra_simulated_omics/raw/genotypes.txt")[:,2] #genotypes data in Christensenet al.(2021)
genptype = zeros(length(geno),15000)
for i in 1:21100
    indi=parse.(Int32,split(geno[i],""))
    genptype[i,:] = indi
end

for rep in 1:20
    selectIndiv=Int.(vec(readdlm("data$rep/select_1055_index_data$rep.txt")))

    # save genotype
    println("number of individuals is ", size(genptype,1), ", number of SNP is ", size(genptype,2))
    genptype_select=genptype[selectIndiv,:]
    println("after selection, number of individuals is ", size(genptype_select,1), ", number of SNP is ", size(genptype_select,2))
    writedlm("data$rep/geno.txt",genptype_select)
    names=["snp$i" for i in 1:size(genptype_select,2)]
    genptype_df=DataFrame(genptype_select,names)
    insertcols!(genptype_df, 1, :ID => selectIndiv)
    CSV.write("data$rep/geno.csv", genptype_df)
end

#save tbv (=tbv_m+ar)
using DataFrames,CSV,Random,DelimitedFiles,Distributions
cd("/Users/tianjing/Box/Omicsdata/andrew1055")
tbv=readdlm("/Users/tianjing/Box/Omicsdata/andrewlegarra_simulated_omics/raw/tbv.txt")
tbv=tbv[:,3]
#save tbv_m
for rep in 1:20
    println("rep$rep------------")
    selectIndiv=Int.(vec(readdlm("data$rep/select_1055_index_data$rep.txt")))
    tbv_select = tbv[selectIndiv]
    writedlm("data$rep/tbv.txt",tbv_select)
    tbv_select_df=DataFrame(ID=selectIndiv, tbv=tbv_select)
    CSV.write("data$rep/tbv.csv", tbv_select_df)
end

max_stars_repo_name: zhaotianjing/NN-LMM

max_issues_repo_name: zhaotianjing/NN-LMM

max_forks_repo_name: zhaotianjing/NN-LMM


text:
"""
    cfactor!(A::AbstractMatrix)

A slightly modified version of `chol!` from `Base`

Uses `inject!` (as opposed to `copy!`), `downdate!` (as opposed to `syrk!` or `gemm!`)
and recursive calls to `cfactor!`.

Note: The `cfactor!` method for dense matrices calls `LAPACK.potrf!` directly to avoid
errors being thrown when `A` is computationally singular
"""
function cfactor!(A::AbstractMatrix)
    n = Compat.LinAlg.checksquare(A)
    for k = 1:n
        Akk = A[k, k]
        for i in 1:(k - 1)
            downdate!(Akk, A[i, k])  # A[k,k] -= A[i,k]'A[i,k]
        end
        Akk = cfactor!(Akk)          # right Cholesky factor of A[k,k]
        for j in (k + 1):n
            for i in 1:(k - 1)
                downdate!(A[k, j], A[i, k], A[i, j]) # A[k,j] -= A[i,k]'A[i,j]
            end
            LinAlg.Ac_ldiv_B!(Akk, A[k, j])
        end
    end
    UpperTriangular(A)
end

function cfactor!{T <: AbstractFloat}(D::Diagonal{T})
    dd = D.diag
    for i in eachindex(dd)
        dd[i] = sqrt(dd[i])
    end
    D
end

cfactor!{T <: LinAlg.BlasFloat}(A::Matrix{T}) = UpperTriangular(LAPACK.potrf!('U', A)[1])

function cfactor!{T}(A::HBlkDiag{T})
    Aa = A.arr
    r, s, t = size(Aa)
    if r ≠ s
        throw(ArgumentError("HBlkDiag matrix A must be square"))
    end
    scm = Array(T, (r, r))
    for k in 1 : t  # FIXME: Lots of allocations in this loop
        for j in 1 : r, i in 1 : j
            scm[i, j] = Aa[i, j, k]
        end
        LAPACK.potrf!('U', scm)
        for j in 1 : r, i in 1 : j
            Aa[i, j, k] = scm[i, j]
        end
    end
    UpperTriangular(A)
end

"""
    downdate!(C::AbstractMatrix, A::AbstractMatrix)
    downdate!(C::AbstractMatrix, A::AbstractMatrix, B::AbstractMatrix)

Subtracts, in place, `A'A` or `A'B` from `C`
"""
function downdate! end

downdate!{T <: LinAlg.BlasFloat}(C::DenseMatrix{T}, A::DenseMatrix{T}) =
    BLAS.syrk!('U', 'T', -one(T), A, one(T), C)

downdate!{T <: LinAlg.BlasFloat}(C::DenseMatrix{T}, A::DenseMatrix{T}, B::DenseMatrix{T}) =
    BLAS.gemm!('T', 'N', -one(T), A, B, one(T), C)

function downdate!{T}(C::Diagonal{T}, A::SparseMatrixCSC{T})
    m, n = size(A)
    dd = C.diag
    if length(dd) ≠ n
        throw(DimensionMismatch("size(C,2) ≠ size(A,2)"))
    end
    nz = nonzeros(A)
    for j in eachindex(dd)
        for k in nzrange(A, j)
            @inbounds dd[j] -= abs2(nz[k])
        end
    end
    C
end

function downdate!{T}(C::Diagonal{T},A::Diagonal{T})
    if size(C) ≠ size(A)
        throw(DimensionMismatch("size(C) ≠ size(A)"))
    end
    c, a = C.diag, A.diag
    for i in eachindex(c)
        c[i] -= abs2(a[i])
    end
    C
end

function downdate!{T}(C::Diagonal{T}, A::Diagonal{T}, B::Diagonal{T})
    if !(size(C) == size(A) == size(B))
        throw(DimensionMismatch("need size(C) == size(A) == size(B)"))
    end
    c, a, b = C.diag, A.diag, B.diag
    for i in eachindex(c)
        c[i] -= a[i] * b[i]
    end
    C
end

function downdate!{T}(C::DenseMatrix{T}, A::Diagonal{T}, B::DenseMatrix{T})
    a = A.diag
    if ((m, n) = size(B)) ≠ size(C)
        throw(DimensionMismatch("size(B) ≠ size(C)"))
    end
    if length(a) ≠ m
        throw(DimensionMismatch("size(A,2) ≠ size(B,1)"))
    end
    for j in 1:n, i in 1:m
        C[i,j] -= a[i] * B[i,j]
    end
    C
end

function downdate!{T}(C::DenseMatrix{T}, A::SparseMatrixCSC{T}, B::DenseMatrix{T})
    m, n = size(A)
    r, s = size(C)
    if r ≠ n || s ≠ size(B,2) || m ≠ size(B,1)
        throw(DimensionMismatch("size(C,1) ≠ size(A,2) or size(C,2) ≠ size(B,2) or size(A,1) ≠ size(B,1)"))
    end
    nz = nonzeros(A)
    rv = rowvals(A)
    for jj in 1:s, j in 1:n, k in nzrange(A,j)
        C[j,jj] -= nz[k] * B[rv[k], jj]
    end
    C
end

function downdate!{T}(C::DenseMatrix{T}, A::SparseMatrixCSC{T}, B::SparseMatrixCSC{T})
    AtB = A'B
    if size(C) ≠ size(AtB)
        throw(DimensionMismatch("size(C) ≠ size(A'B)"))
    end
    atbrv = rowvals(AtB)
    atbnz = nonzeros(AtB)
    for j in 1:size(AtB, 2)
        for k in nzrange(AtB, j)
            C[atbrv[k], j] -= atbnz[k]
        end
    end
    C
end

function downdate!{T}(C::DenseMatrix{T}, A::SparseMatrixCSC{T})
    m, n = size(A)
    if n ≠ Compat.LinAlg.checksquare(C)
        throw(DimensionMismatch("C is not square or size(C,2) ≠ size(A,2)"))
    end
    # FIXME: avoid allocation by caching a transposed matrix and just fill in the new values
    # alternatively, work with the lower Cholesky factor L instead of R
    At = A'
    rv = rowvals(A)
    nz = nonzeros(A)
    rvt = rowvals(At)
    nzt = nonzeros(At)
    cp = A.colptr
    @inbounds for j in 1:n
        for jp in nzrange(A, j)
            nzB = nz[jp]
            k = rv[jp]
            for kp in nzrange(At, k)
                C[rvt[kp], j] -= nzt[kp] * nzB
            end
        end
    end
    C
end

## FIXME: Need a downdate! method for SparseMatrixCSC, SparseMatrixCSC, SparseMatrixCSC (3 or more nested)

"""
    inflate!(A::HblkDiag)
    inflate!(A::Diagonal)
    inflate!(A::StridedMatrix)

Adds an identity to `A` in place.
"""
function inflate! end

function inflate!(A::HBlkDiag)
    Aa = A.arr
    r, s, k = size(Aa)
    for j in 1 : k, i in 1 : min(r, s)
        Aa[i, i, j] += 1
    end
    A
end
function inflate!{T<:AbstractFloat}(A::StridedMatrix{T})
    n = Compat.LinAlg.checksquare(A)
    for i in 1 : n
        @inbounds A[i, i] += one(T)
    end
    A
end
function inflate!{T <: AbstractFloat}(D::Diagonal{T})
    d = D.diag
    for i in eachindex(d)
        d[i] += one(T)
    end
    D
end


text:
using Frames
import MAT
function setupframes(path)
	#include("/home/fredrikb/robotlib/frames/Frames.jl")
	# path = "/home/fredrikb/work/flexifab/frames/"

	# Add frame names to the dictionary
	add_frame_name!("SEAM","Weld seam frame")
	add_frame_name!("TAB","Table frame")

	# Read matrices from file
	T_RB_Tm = MAT.matread(path*"T_RB_T.mat")["T_RB_T"]
	T_TF_TCPm = MAT.matread(path*"T_TF_TCP.mat")["T_TF_TCP"]
	T_T_TABm = MAT.matread(path*"T_T_Table.mat")["T_T_Table"]

	# Create frames from matrices
	T_RB_T = Frame(T_RB_Tm,"RB","T")
	T_S_D = Frame(T_TF_TCPm,"S","D")
	T_T_TAB = Frame(T_T_TABm,"T","TAB")

	# Read point clouds generated by nikon software from file
	cloud_seam = readcloud(path*"CloudSeam_edge.txt")
	plane_seam = readplane(path*"PlaneSeam_edge.txt")

	# Project points onto plane and fit a line
	cloud_seam_projected = project(plane_seam,cloud_seam)
	line_seam = fitline(cloud_seam_projected)

	# Create a frame from the measured features
	T_T_SEAM = framefromfeatures(("z+",line_seam),("y-",plane_seam),cloud_seam_projected[1],"SEAM")
	T_RB_SEAM = T_RB_T*T_T_SEAM
	T_RB_TAB = T_RB_T*T_T_TAB
	T_TAB_SEAM = inv(T_T_TAB)*T_T_SEAM

	# Write frames to file
	MAT.matwrite(path*"T_TAB_SEAM.mat",["T_TAB_SEAM" => T_TAB_SEAM.T])
	MAT.matwrite(path*"T_T_SEAM.mat",["T_T_SEAM" => T_T_SEAM.T])
	MAT.matwrite(path*"T_RB_TAB.mat",["T_RB_TAB" => T_RB_TAB.T])
	println("Wrote T_TAB_SEAM, T_T_SEAM, T_RB_TAB to files in $path")

	cloud_seam_RB = T_RB_T*cloud_seam
	cloud_seam_projected_RB = T_RB_T*cloud_seam_projected
	plane_seam_RB = T_RB_T*plane_seam
	line_seam_RB = T_RB_T*line_seam

	# Plot results
	plotframe(Frame(eye(4),"RB","U"),200, label=true)
	plotpoints(cloud_seam_RB,"bx")
	plotpoints(cloud_seam_projected_RB,"rx")
	plotline(line_seam_RB,"r",500,label="Line seam")
	plotplane(plane_seam_RB,"b",200,label="Plane seam")
	plotframe(T_RB_SEAM,200, label=true)
	plotframe(T_RB_TAB,200, label=true)

	PyPlot.xlabel("x")
	PyPlot.ylabel("y")
	PyPlot.zlabel("z")
	PyPlot.axis("scaled")

	if true
		add_frame_name!("SEAM2","Frame for second part of Weld seam")
		cloud_seam2 = readcloud(path*"CloudSeam2_edge.txt")
		cloud_seam2_projected = project(plane_seam,cloud_seam2)
		line_seam2 = fitline(cloud_seam2_projected)

		T_T_SEAM2 = framefromfeatures(("z+",line_seam2),("y-",plane_seam),cloud_seam2_projected[1],"SEAM2")
		T_RB_SEAM2 = T_RB_T*T_T_SEAM2
		T_TAB_SEAM2 = inv(T_T_TAB)*T_T_SEAM2;
		MAT.matwrite("T_TAB_SEAM2.mat",["T_TAB_SEAM2" => T_TAB_SEAM2.T])
		MAT.matwrite("path*"T_T_SEAM2.mat",["T_T_SEAM2" => T_T_SEAM2.T])
		println("Wrote T_TAB_SEAM2, T_T_SEAM2 to files in $path")

		cloud_seam2_RB = T_RB_T*cloud_seam2
		cloud_seam2_projected_RB = T_RB_T*cloud_seam2_projected
		plane_seam2_RB = T_RB_T*plane_seam
		line_seam2_RB = T_RB_T*line_seam2

		plotpoints(cloud_seam2_RB,"bx")
		plotpoints(cloud_seam2_projected_RB,"rx")
		plotline(line_seam2_RB,"r",500,label="Line seam2")
		plotplane(plane_seam2_RB,"b",200,label="Plane seam2")
		plotframe(T_RB_SEAM2,200, label=true)
	end

end


text:
#Evaluating the likelihood requires only filtering ("forward").

# WRAPPER FUNCTIONS

#data: several iid individuals
function loglikelihood{S<:Array}(model::DynamicDiscreteModel,data::Array{S,1})
	n=length(data)
	T=0.0
	llk=0.0
	for i=1:n
		Ti=float(length(data[i]))
		llk+=loglikelihood(model,data[i])*Ti
		T+=Ti
	end
	llk/T
end

#data: several iid individuals
function loglikelihood_jac{S<:Array}(model::DynamicDiscreteModel,data::Array{S,1})
	dtheta=size(model.rhojac)
	n=length(data)
	T=0.0
	llk=0.0
	llkjac=zeros(dtheta)
	for i=1:n
		Ti=float(length(data[i]))
		temp=loglikelihood_jac(model,data[i])
		llk+=temp[1]*Ti
		llkjac+=temp[2]*Ti
		T+=Ti
	end
	llk/T,llkjac/T
end

#simulate iid individuals with heterogeneous number of periods of observation
function rand(model::DynamicDiscreteModel,T::Array{Int,1})
	n=length(T)
	data=Array(Array,n)
	for i=1:n
		data[i]=rand(model,T[i])
	end
	data
end

# (T,n) -> n iid individuals with T periods of observation each
rand(model::DynamicDiscreteModel,T,n)=rand(model,fill(T,n))
# rand(model::StatisticalModel,Tn::Tuple{Int,Int})=rand(model,Tn[1],Tn[2])

# CORE NON-JACOBIAN METHODS

function filterstep!(model::DynamicDiscreteModel,iy,jy)
	dx=length(model.phi)

	#jx=1
	#phi is normalized from the previous step
	model.phi[1]=model.psi[1]/model.rho[1]
	model.psi[1]=model.m[1,iy,1,jy]*model.phi[1]
	for ix=2:dx
		model.phi[ix]=model.psi[ix]/model.rho[1]
		model.psi[1]+=model.m[ix,iy,1,jy]*model.phi[ix]
	end
	model.rho[1]=model.psi[1]

	for jx=2:dx
		model.psi[jx]=model.m[1,iy,jx,jy]*model.phi[1]
		for ix=2:dx
			model.psi[jx]+=model.m[ix,iy,jx,jy]*model.phi[ix]
		end
		model.rho[1]+=model.psi[jx]
	end
end

function loglikelihood(model::DynamicDiscreteModel,data::Array{Int,1})
	lambda=log(sum(model.mu[:,data[1]]))				#log-normalization factor for numerical stability

	model.psi[:]=model.mu[:,data[1]]/sum(model.mu[:,data[1]])		#actual filter
	model.rho[1]=1

	for t=2:length(data)
		filterstep!(model,data[t-1],data[t])
		lambda+=log(model.rho[1])
	end
	lambda/length(data)
end

# CORE JACOBIAN METHODS

function filterstep_jac!(model::DynamicDiscreteModel,iy,jy)
	dx,dtheta=size(model.phijac)

	#jx=1
	#phi is normalized from the previous step
	model.phi[1]=model.psi[1]/model.rho[1]
	model.psi[1]=model.m[1,iy,1,jy]*model.phi[1]
	for itheta=1:dtheta
		model.phijac[1,itheta]=(model.psijac[1,itheta]-model.phi[1]*model.rhojac[itheta])/model.rho[1]
		model.psijac[1,itheta]=model.mjac[1,iy,1,jy,itheta]*model.phi[1]+model.m[1,iy,1,jy]*model.phijac[1,itheta]
	end
	for ix=2:dx
		model.phi[ix]=model.psi[ix]/model.rho[1]
		model.psi[1]+=model.m[ix,iy,1,jy]*model.phi[ix]
		for itheta=1:dtheta
			model.phijac[ix,itheta]=(model.psijac[ix,itheta]-model.phi[ix]*model.rhojac[itheta])/model.rho[1]
			model.psijac[1,itheta]+=model.mjac[ix,iy,1,jy,itheta]*model.phi[ix]+model.m[ix,iy,1,jy]*model.phijac[ix,itheta]
		end
	end
	model.rho[1]=model.psi[1]
	for itheta=1:dtheta
		model.rhojac[itheta]=model.psijac[1,itheta]
	end

	for jx=2:dx
		model.psi[jx]=model.m[1,iy,jx,jy]*model.phi[1]
		for itheta=1:dtheta
			model.psijac[jx,itheta]=model.mjac[1,iy,jx,jy,itheta]*model.phi[1]+model.m[1,iy,jx,jy]*model.phijac[1,itheta]
		end
		for ix=2:dx
			model.psi[jx]+=model.m[ix,iy,jx,jy]*model.phi[ix]
			for itheta=1:dtheta
				model.psijac[jx,itheta]+=model.mjac[ix,iy,jx,jy,itheta]*model.phi[ix]+model.m[ix,iy,jx,jy]*model.phijac[ix,itheta]
			end
		end
		model.rho[1]+=model.psi[jx]
		for itheta=1:dtheta
			model.rhojac[itheta]+=model.psijac[jx,itheta]
		end
	end
end

function loglikelihood_jac(model::DynamicDiscreteModel,data::Array{Int,1})
	lambda=log(sum(model.mu[:,data[1]]))				#log-normalization factor for numerical stability
	dx,dtheta=size(model.phijac)
	lambdajac=zeros(dtheta)

	model.psi[:]=model.mu[:,data[1]]/sum(model.mu[:,data[1]])		#actual filter
	model.rho[1]=1
	model.psijac[:]=zeros(dx,dtheta)		#actual filter
	model.rhojac[:]=zeros(dtheta)

	for t=2:length(data)
		filterstep_jac!(model,data[t-1],data[t])
		lambda+=log(model.rho[1])
		for itheta=1:dtheta
			lambdajac[itheta]+=model.rhojac[itheta]/model.rho[1]
		end
	end
	lambda/length(data),lambdajac/length(data)
end

#non efficient log-likelihood + jacobian implementation. keep around for test purposes
function loglikelihood_jac_nonefficient(model::DynamicDiscreteModel,data::Array{Int,1})
	dx,dy=size(model.mu)
	dtheta=size(model.mjac)[5] 					 #dimension of the parameter
	rho=log(sum(model.mu[:,data[1]]))				#log-normalization factor for numerical stability
	rhojac=zeros(1,dtheta)
	phi=model.mu[:,data[1]]/sum(model.mu[:,data[1]])		#actual filter
	phijac=zeros(dx,dtheta)
	phihatjac=zeros(dx,dtheta)

	for t=2:length(data)
		mm=reshape(model.m[:,data[t-1],:,data[t]],(dx,dx))'
		phihat=mm*phi
		for itheta=1:dtheta
			mmjac=reshape(model.mjac[:,data[t-1],:,data[t],itheta],(dx,dx))'
			phihatjac[:,itheta]=mm*phijac[:,itheta]+mmjac*phi
		end
		lambda=sum(phihat)
		lambdajac=sum(phihatjac,1)
		phi=phihat/lambda
		rho+=log(lambda)
		for itheta=1:dtheta
			phijac[:,itheta]=(phihatjac[:,itheta]*lambda - phihat*lambdajac[itheta])/lambda^2
			rhojac[itheta]+=lambdajac[itheta]/lambda
		end
	end
	rho,rhojac
end


text:
struct Wilson_GeneralDirac_FermiAction{Dim,Dirac,fermion,gauge} <: Wilsontype_FermiAction{Dim,Dirac,fermion,gauge}
    hascovnet::Bool
    covneuralnet::Union{Nothing,CovNeuralnet{Dim}}
    diracoperator::Dirac
    _temporary_fermionfields::Vector{fermion}
    _temporary_gaugefields::Vector{gauge}

    function Wilson_GeneralDirac_FermiAction(D::Dirac_operator{Dim},hascovnet,covneuralnet,parameters_action) where Dim

        num = 8

        x = D._temporary_fermi[1]
        xtype = typeof(x)
        _temporary_fermionfields = Array{xtype,1}(undef,num)
        for i=1:num
            _temporary_fermionfields[i] = similar(x)
        end

        Utemp = D.U[1]
        Utype = typeof(Utemp)
        numU = 2
        _temporary_gaugefields = Array{Utype,1}(undef,numU)
        for i=1:numU
            _temporary_gaugefields[i] = similar(Utemp)
        end

        return new{Dim,typeof(D),xtype,Utype}(hascovnet,covneuralnet,D,
                        _temporary_fermionfields,
                        _temporary_gaugefields)

    end
end

function calc_dSfdU!(dSfdU::Vector{<: AbstractGaugefields},fermi_action::Wilson_GeneralDirac_FermiAction{Dim,Dirac,fermion,gauge} ,U::Vector{<: AbstractGaugefields},ϕ::AbstractFermionfields) where  {Dim,Dirac,fermion,gauge,hascloverterm}
    #println("------dd")
    W = fermi_action.diracoperator(U)
    WdagW = DdagD_Wilson_GeneralDirac_operator(W)
    X = fermi_action._temporary_fermionfields[5]
    Y = fermi_action._temporary_fermionfields[6]
    #X = (D^dag D)^(-1) ϕ
    #
    #println("Xd ",X[1,1,1,1,1,1])

    solve_DinvX!(X,WdagW,ϕ)
    #set_wing_fermion!(X)
    #println("X ",X[1,1,1,1,1,1])
    clear_U!(dSfdU)

    calc_dSfdU_fromX!(dSfdU,Y,fermi_action,U,X)
    #println("----aa--")
    set_wing_U!(dSfdU)
end

function calc_dSfdU_fromX!(dSfdU::Vector{<: AbstractGaugefields},Y,fermi_action::Wilson_GeneralDirac_FermiAction{Dim,Dirac,fermion,gauge} ,U,X;coeff = 1) where {Dim,Dirac,fermion,gauge}
    W = fermi_action.diracoperator(U)
    mul!(Y,W,X)
    #set_wing_fermion!(Y)

    temp0_f = fermi_action._temporary_fermionfields[3]
    temp1_f = fermi_action._temporary_fermionfields[4]
    temp0_g = fermi_action._temporary_gaugefields[end]
    #temp1_g = fermi_action._temporary_gaugefields[2]

    #set_wing_U!(U)

    numterms = get_numterms(W)

    #=

    for μ=1:Dim
        #!  Construct U(x,mu)*P1

        # U_{k,μ} X_{k+μ}
        Xplus = shift_fermion(X,μ)
        #@time mul!(temp0_f,U[μ],X)
        diracterm = W[1]
        coeff_i = get_coefficient(diracterm)

        mul!(temp1_f,rminusγ1[:,:,μ],Xplus) #Xplus Y^dag

        # κ (r-γ_μ) U_{k,μ} X_{k+μ}
        mul!(temp0_f,coeff_i,temp1_f)

        # κ ((r-γ_μ) U_{k,μ} X_{k+μ}) ⊗ Y_k
        cross!(temp0_g,Y',temp0_f)
        add_U!(dSfdU[μ],-coeff,temp0_g)

        Yplus = shift_fermion(Y,μ) #X   Yplus^dag U_{k,μ}^dag * rplusγ1[:,:,μ]
        mul!(temp1_f,rplusγ1[:,:,μ],Yplus)
        mul!(temp0_f,coeff_i,temp1_f)

        cross!(temp0_g,X',temp0_f)
        add_U!(dSfdU[μ],-coeff,temp0_g)
        #add_U!(dSfdU[μ],-coeff,temp0_g')

        #= memo
        -Ydag otimes rminusγ1*Xplus
        -Xdag otimes rplusγ1*Yplus

        =#
    end
    return

    =#

    diracterm = W[1]
    coeff_i = get_coefficient(diracterm)
    #clear_U!(dSfdU)
    #return

    #=

    for μ=1:Dim
        Xplus = shift_fermion(X,μ)
        #println(rplusγ1[:,:,μ])
        mul!(temp0_f,rminusγ1[:,:,μ],Xplus)
        #println(Xplus[1,1,1,1,1,1])
        #println(temp0_f[1,1,1,1,1,1])
        cross!(temp0_g,Y',temp0_f)

        #mul!(temp0_g,Y',temp0_f)
       # println(temp0_g[:,:,1,1,1,1])

        add_U!(dSfdU[μ],-coeff*coeff_i,temp0_g)
        println("c 1, ", temp0_g[1,1,1,1,1,1])

        Yplus = shift_fermion(Y,μ)
        mul!(temp0_f,rplusγ1[:,:,μ],Yplus)
        cross!(temp0_g,X',temp0_f)
        #mul!(temp0_g,X',temp0_f)
        add_U!(dSfdU[μ],-coeff*coeff_i,temp0_g)
        println("c 2, ",temp0_g[1,1,1,1,1,1])

    end
    =#

    #return

    #error("dd")

    for i=1:numterms
        diracterm = W[i]
        position  = collect(get_fermionposition(diracterm))
        Sdag = get_spinormatrix_dagger(diracterm)
        S = get_spinormatrix(diracterm)
        #γ5S = get_spinormatrix_γ5M(diracterm)
        coeff_i = get_coefficient(diracterm)
        #println("i = $i ",coeff_i)
        #if coeff_i < 1e-15
        #    continue
        #end
        for μ=1:Dim
            derivatives = get_derivatives(diracterm,μ)
            linkinfo_derivatives_leftdag = get_linkinfo_derivatives_leftdag(diracterm,μ)
            linkinfo_derivatives_right = get_linkinfo_derivatives_right(diracterm,μ)
            numderivatives = length(derivatives)

            for inum=1:numderivatives
                dFdμ_i = derivatives[inum]
                #dFdμ_i_position_left = collect(dFdμ_i.position)
                dFdμ_i_position_left = -collect(dFdμ_i.position)
                dFdμ_i_position_right = dFdμ_i_position_left .+ position

                #println(Tuple(dFdμ_i_position_left))
                #println(Tuple(dFdμ_i_position_right))
                #if coeff_i > 1e-15
                Xshifted = shift_fermion(X,Tuple(dFdμ_i_position_right))
                Yshifted = shift_fermion(Y,Tuple(dFdμ_i_position_left))

                #if coeff_i > 1e-15
                gaugefield_fermion_mul!(temp1_f,linkinfo_derivatives_right[inum],U,Xshifted,fermi_action._temporary_fermionfields)
                #end
                #set_wing_fermion!(temp1_f)
                mul!(temp1_f,S)
                #set_wing_fermion!(temp1_f)

                #if coeff_i > 1e-15
                gaugefield_fermion_mul!(temp0_f,linkinfo_derivatives_leftdag[inum],U,Yshifted,fermi_action._temporary_fermionfields)
                #end
                #set_wing_fermion!(temp0_f)
                cross!(temp0_g,temp0_f',temp1_f)
                #set_wing_U!(temp0_g)
                #println("1, ",coeff*coeff_i*temp0_g[:,:,1,1,1,1])
                #mul!(temp0_g,temp1_f,temp0_f')

                #if coeff_i > 1e-15
                add_U!(dSfdU[μ],-coeff*coeff_i,temp0_g)
                #end
                #println("1, μ = $μ ", dSfdU[μ][1,1,1,1,1,1])
                #set_wing_U!(dSfdU[μ])

                #if coeff_i > 1e-15

                Xshifted = shift_fermion(X,Tuple(dFdμ_i_position_left))
                Yshifted = shift_fermion(Y,Tuple(dFdμ_i_position_right))

                gaugefield_fermion_mul!(temp1_f,linkinfo_derivatives_right[inum],U,Yshifted,fermi_action._temporary_fermionfields)
                #set_wing_fermion!(temp1_f)
                mul!(temp1_f,Sdag)
                #set_wing_fermion!(temp1_f)

                gaugefield_fermion_mul!(temp0_f,linkinfo_derivatives_leftdag[inum],U,Xshifted,fermi_action._temporary_fermionfields)
                #set_wing_fermion!(temp0_f)
                cross!(temp0_g,temp0_f',temp1_f)
                #set_wing_U!(temp0_g)

                #println("1, ",coeff*coeff_i*temp0_g[:,:,1,1,1,1])
                #mul!(temp0_g,temp1_f,temp0_f')

                add_U!(dSfdU[μ],-coeff*coeff_i,temp0_g)
                #println("2, μ = $μ ", dSfdU[μ][1,1,0,1,1,1])
                #end

                #set_wing_U!(dSfdU[μ])

            end

            #=

            derivatives = get_derivatives_dag(diracterm,μ)
            linkinfo_derivatives_leftdag = get_linkinfo_derivatives_dag_leftdag(diracterm,μ)
            linkinfo_derivatives_right = get_linkinfo_derivatives_dag_right(diracterm,μ)
            numderivatives = length(derivatives)

            for inum=1:numderivatives
                dFdμ_i = derivatives[inum]
                dFdμ_i_position_left = collect(dFdμ_i.position)
                dFdμ_i_position_right = dFdμ_i_position_left .+ position

                Xshifted = shift_fermion(X,Tuple(dFdμ_i_position_left))
                Yshifted = shift_fermion(Y,Tuple(dFdμ_i_position_right))

                gaugefield_fermion_mul!(temp1_f,linkinfo_derivatives_right[inum],U,Yshifted,fermi_action._temporary_fermionfields)
                mul!(temp1_f,S')

                gaugefield_fermion_mul!(temp0_f,linkinfo_derivatives_leftdag[inum],U,Xshifted,fermi_action._temporary_fermionfields)
                cross!(temp0_g,temp0_f',temp1_f)
                #mul!(temp0_g,temp1_f,temp0_f')

                add_U!(dSfdU[μ],-coeff*coeff_i,temp0_g)
                println("2, ", temp0_g[1,1,1,1,1,1])

                #println("2, ",coeff*coeff_i*temp0_g[:,:,1,1,1,1])

            end
            =#
        end

        #show(diracterm)

        #D_positions = get_positions(diracterm)
        #D_directions = get_directions(diracterm)
        #D_isdagvectors = get_isdagvectors(diracterm)

    end

    return

    #=
        for μ=1:Dim
            derivatives = get_derivatives(diracterm,μ)
            linkinfo_derivatives_leftdag = get_linkinfo_derivatives_leftdag(diracterm,μ)
            linkinfo_derivatives_right = get_linkinfo_derivatives_right(diracterm,μ)
            numderivatives = length(derivatives)
            if numderivatives != 0
                println("i =$i μ = $μ, num = $numderivatives")
                show( derivatives)
            end

            for inum=1:numderivatives
                dFdμ_i = derivatives[inum]
                dFdμ_i_left = get_leftlinks(dFdμ_i)
                dFdμ_i_right = get_rightlinks(dFdμ_i)
                dFdμ_i_position_left = collect(dFdμ_i.position)
                dFdμ_i_position_right = dFdμ_i_position_left .+ position
                #show(dFdμ_i_left )
                #show(dFdμ_i_right )
                #println(dFdμ_i_position_right )
                #println(dFdμ_i_position_left )
                println(Tuple(dFdμ_i_position_left))
                println(Tuple(dFdμ_i_position_right))
                Xshifted = shift_fermion(X,Tuple(dFdμ_i_position_right))
                Yshifted = shift_fermion(Y,Tuple(dFdμ_i_position_left))

                gaugefield_fermion_mul!(temp1_f,linkinfo_derivatives_right[inum],U,Xshifted,fermi_action._temporary_fermionfields)
                mul!(temp1_f,S)

                gaugefield_fermion_mul!(temp0_f,linkinfo_derivatives_leftdag[inum],U,Yshifted,fermi_action._temporary_fermionfields)
                mul!(temp0_g,temp0_f',temp1_f)
                #mul!(temp0_g,temp1_f,temp0_f')

                add_U!(dSfdU[μ],coeff*coeff_i,temp0_g)
            end

            derivatives_dag = get_derivatives_dag(diracterm,μ)

            linkinfo_derivatives_dag_leftdag = get_linkinfo_derivatives_dag_leftdag(diracterm,μ)
            linkinfo_derivatives_dag_right = get_linkinfo_derivatives_dag_right(diracterm,μ)
            numderivatives = length(derivatives_dag)
            for inum=1:numderivatives
                println("position = ",derivatives_dag[inum].position)
            end

            if numderivatives != 0
                println("dag i =$i μ = $μ, num = $numderivatives")
                show( derivatives_dag)
            end

            for inum=1:numderivatives
                dFdμ_i = derivatives_dag[inum]
                dFdμ_i_left = get_leftlinks(dFdμ_i)
                dFdμ_i_right = get_rightlinks(dFdμ_i)
                dFdμ_i_position_left = collect(dFdμ_i.position)
                dFdμ_i_position_right = dFdμ_i_position_left .+ position
                println(Tuple(dFdμ_i_position_left))
                println(Tuple(dFdμ_i_position_right))
                Xshifted = shift_fermion(X,Tuple(-dFdμ_i_position_left))
                Yshifted = shift_fermion(Y,Tuple(-dFdμ_i_position_right))

                gaugefield_fermion_mul!(temp1_f,linkinfo_derivatives_dag_right[inum],U,Yshifted,fermi_action._temporary_fermionfields)

                mul!(temp1_f,S')

                gaugefield_fermion_mul!(temp0_f,linkinfo_derivatives_dag_leftdag[inum],U,Xshifted,fermi_action._temporary_fermionfields)
                mul!(temp0_g,temp0_f',temp1_f)
                #mul!(temp0_g,temp1_f,temp0_f')

                add_U!(dSfdU[μ],-coeff*coeff_i,temp0_g)
            end
        end

    end

    error("dte")

    =#
    return

    for μ=1:Dim
        #!  Construct U(x,mu)*P1

        # U_{k,μ} X_{k+μ}
        Xplus = shift_fermion(X,μ)

        #@time mul!(temp0_f,U[μ],X)

        mul!(temp0_f,U[μ],Xplus)

        # (r-γ_μ) U_{k,μ} X_{k+μ}
        mul!(temp1_f,view(W.rminusγ,:,:,μ),temp0_f)

        # κ (r-γ_μ) U_{k,μ} X_{k+μ}
        mul!(temp0_f,W.hopp[μ],temp1_f)

        # κ ((r-γ_μ) U_{k,μ} X_{k+μ}) ⊗ Y_k
        mul!(temp0_g,temp0_f,Y')

        add_U!(UdSfdU[μ],-coeff,temp0_g)

        #!  Construct P2*U_adj(x,mu)
        # Y_{k+μ}^dag U_{k,μ}^dag
        Yplus = shift_fermion(Y,μ)
        mul!(temp0_f,Yplus',U[μ]')

        # Y_{k+μ}^dag U_{k,μ}^dag*(r+γ_μ)
        mul!(temp1_f,temp0_f,view(W.rplusγ,:,:,μ))

        # κ Y_{k+μ}^dag U_{k,μ}^dag*(r+γ_μ)
        mul!(temp0_f,W.hopm[μ],temp1_f)

        # X_k ⊗ κ Y_{k+μ}^dag U_{k,μ}^dag*(r+γ_μ)
        mul!(temp0_g,X,temp0_f)

        add_U!(UdSfdU[μ],coeff,temp0_g)

    end

end


text:
################################################################################
#
#  Ramification
#
################################################################################

export schur_index, is_eichler

################################################################################
#
#  Is split
#
################################################################################

@doc Markdown.doc"""
    is_split(A::AbsAlgAss, p) -> Bool

Return whether the $\mathbf{Q}$-algebra $A$ is split at $p$. The object $p$ can be an integer or `inf`.
"""
function is_split(A::AbsAlgAss, p)
  return schur_index(A, p) == 1
end

@doc Markdown.doc"""
    is_split(A::AlgAss{fmpq}) -> Bool

Given a central $\mathbf{Q}$-algebra $A$, return `true` if $A$ splits.
"""
function is_split(A::AbsAlgAss{fmpq})
  i = schur_index(A, inf)
  if i == 2
    @vprint :AlgAssOrd 1 "Not split at the infinite prime\n"
    return false
  end
  O = any_order(A)
  fac = factor(root(abs(discriminant(O)),2))
  for (p,j) in fac
    O1 = pmaximal_overorder(O, Int(p))
    if valuation(discriminant(O1), Int(p)) != 0
      @vprint :AlgAssOrd 1 "Not split at $p\n"
      return false
    end
  end
  return true
end

function is_split(A::AbsAlgAss{nf_elem})
  K = base_ring(A)
  for p in infinite_places(K)
    if !is_split(A, p)
      return false
    end
  end
  O1 = maximal_order(A)
  return isone(discriminant(O1))
end

function is_split(A::AbsAlgAss, P::InfPlc)
  if iscomplex(P)
    return true
  end
  return schur_index(A, P) == 1
end

################################################################################
#
#  Ramified infinite places
#
################################################################################

function ramified_infinite_places(A::AlgAss{nf_elem})
  K = base_ring(A)
  inf_plc = Vector{InfPlc}()
  places = real_places(K)
  for p in places
    if !is_split(A, p)
      push!(inf_plc, p)
    end
  end

  return inf_plc
end

function ramified_infinite_places_of_center(A::AbsAlgAss)
  dec = decompose(A)
  C, = center(A)
  res = Vector{InfPlc}[]
  for i in 1:length(dec)
    K, = component(Field, C, i)
    B, = _as_algebra_over_center(dec[i][1])
    plcs = ramified_infinite_places(B)
    L = base_ring(B)
    fl, f = is_isomorphic_with_map(L, K)
    @assert fl
    push!(res, [induce_image(f, p) for p in plcs])
  end

  return res
end

################################################################################
#
#  Schur indices
#
################################################################################

# We follow
# Nebe--Steel, "Recognition of division algebras",
# https://doi.org/10.1016/j.jalgebra.2009.04.026

@doc Markdown.doc"""
   schur_index(A::AlgAss{fmpq}, p::Union{IntegerUnion, PosInf}) -> Int

Determine the Schur index of $A$ at $p$, where $p$ is either a prime or `inf`.
"""
schur_index(A::AbsAlgAss{fmpq}, ::Union{IntegerUnion, PosInf})

function schur_index(A::AbsAlgAss{fmpq}, ::PosInf)
  @req iscentral(A) "Algebra must be central"
  @req issimple(A) "Algebra must be simple"

  if dim(A) % 4 != 0
    return 1
  end

  x = trace_signature(A)
  n = root(dim(A), 2)
  if x[1] == divexact(n * (n + 1), 2)
    return 1
  else
    return 2
  end
end

function schur_index(A::AbsAlgAss{nf_elem}, P::InfPlc)
  @req iscentral(A) "Algebra must be central"
  @req issimple(A) "Algebra must be simple"

  if dim(A) % 4 != 0
    return 1
  end

  x = trace_signature(A, P)
  n = root(dim(A), 2)
  if x[1] == divexact(n * (n + 1),2)
    return 1
  else
    return 2
  end
end

#  Schur Index at p

function schur_index(A::AbsAlgAss, p::IntegerUnion)
  @req iscentral(A) "Algebra must be central"
  @req issimple(A) "Algebra must be simple"

  d = discriminant(maximal_order(A))
  v = valuation(d, p)
  if v == 0
    return 1
  end
  s = root(dim(A), 2)
  t = s - divexact(v, s)
  return divexact(s, t)
end

function schur_index(A::AbsAlgAss{<: NumFieldElem}, p::NumFieldOrdIdl)
  @req iscentral(A) "Algebra must be central"
  @req issimple(A) "Algebra must be simple"

  M = maximal_order(A)
  d = discriminant(maximal_order(A))
  @req order(p) === base_ring(M) "Prime ideal of wrong order"
  v = valuation(d, p)
  if v == 0
    return 1
  end
  s = root(dim(A), 2)
  t = s - divexact(v, s)
  return divexact(s, t)
end

function schur_index(A::AbsAlgAss{fmpq})
  e = schur_index(A, inf)
  for p in prime_divisors(discriminant(maximal_order(A)))
    e = lcm(e, schur_index(A, p))
  end
  return e
end

function schur_index(A::AbsAlgAss{<: NumFieldElem})
  rlp = infinite_places(base_ring(A))
  e = schur_index(A, rlp[1])
  for i in 2:length(rlp)
    e = lcm(e, schur_index(A, rlp[i]))
  end
  for (p, _) in factor(discriminant(maximal_order(A)))
    e = lcm(e, schur_index(A, p))
  end
  return e
end

################################################################################
#
#  Eichler condition
#
################################################################################

# Tests whether A fulfils the Eichler condition relative to the maximal Z-order
# of base_ring(A)
function is_eichler(A::AbsAlgAss{nf_elem})
  @assert issimple(A)
  @assert iscentral(A)

  if !istotally_real(base_ring(A))
    return true
  end

  if dim(A) != 4
    return true
  end

  K = base_ring(A)
  places = real_places(K)
  for p in places
    if is_split(A, p)
      return true
    end
  end
  return false
end

function is_eichler(A::AbsAlgAss{fmpq})
  @assert issimple(A)
  @assert iscentral(A)
  if dim(A) != 4
    return true
  end
  O = Order(A, basis(A))
  return schur_index(A, inf) == 1
end