text:
using Test

for testscen in 1:2
    valdir, scenario, use_permafrost, use_seaice = get_scenario(testscen)
    println(scenario)

    include("../src/main_model.jl")

    m = getpage(scenario, use_permafrost, use_seaice)
    run(m)

    if scenario == "NDCs"
        while m[:Discontinuity,:occurdis_occurrencedummy] != [0.,0.,0.,0.,0.,0.,1.,1.,1.,1.]
            println(m[:Discontinuity,:occurdis_occurrencedummy])
            run(m)
        end
    else
        while m[:Discontinuity,:occurdis_occurrencedummy] != [0.,0.,0.,0.,0.,0.,0.,0.,0.,0.]
            println(m[:Discontinuity,:occurdis_occurrencedummy])
            run(m)
        end
    end

    if updatetestdata
        include("../src/utils/save_parameters.jl")
        savepagedata(m, :Population, :pop_population, "test/validationdata/$valdir/pop_population.csv")
        savepagedata(m, :GDP, :cons_percap_consumption, "test/validationdata/$valdir/cons_percap_consumption.csv")
        savepagedata(m, :GDP, :gdp, "test/validationdata/$valdir/gdp.csv")
        savepagedata(m, :Discontinuity, :rcons_per_cap_DiscRemainConsumption, "test/validationdata/$valdir/rcons_per_cap_DiscRemainConsumption.csv")
        savepagedata(m, :TotalAbatementCosts, :tct_per_cap_totalcostspercap, "test/validationdata/$valdir/tct_per_cap_totalcostspercap.csv")
        savepagedata(m, :AdaptiveCostsEconomic, :ac_adaptivecosts, "test/validationdata/ac_adaptationcosts_economic.csv")
        savepagedata(m, :AdaptiveCostsNonEconomic, :ac_adaptivecosts, "test/validationdata/ac_adaptationcosts_noneconomic.csv")
        savepagedata(m, :AdaptiveCostsSeaLevel, :ac_adaptivecosts, "test/validationdata/$valdir/ac_adaptationcosts_sealevelrise.csv")
        savepagedata(m, :TotalAdaptationCosts, :act_adaptationcosts_total, "test/validationdata/$valdir/act_adaptationcosts_tot.csv")
        savepagedata(m, :TotalAdaptationCosts, :act_percap_adaptationcosts, "test/validationdata/$valdir/act_percap_adaptationcosts.csv")
        savepagedata(m, :SLRDamages, :rcons_per_cap_SLRRemainConsumption, "test/validationdata/$valdir/rcons_per_cap_SLRRemainConsumption.csv")
        savepagedata(m, :SLRDamages, :cons_percap_aftercosts, "test/validationdata/$valdir/cons_percap_aftercosts.csv")
        savepagedata(m, :MarketDamagesBurke, :rcons_per_cap_MarketRemainConsumption, "test/validationdata/$valdir/rcons_per_cap_MarketRemainConsumption.csv")
        savepagedata(m, :MarketDamagesBurke, :rgdp_per_cap_MarketRemainGDP, "test/validationdata/$valdir/rgdp_per_cap_MarketRemainGDP.csv")
        savepagedata(m, :NonMarketDamages, :rcons_per_cap_NonMarketRemainConsumption, "test/validationdata/$valdir/rcons_per_cap_NonMarketRemainConsumption.csv")
        savepagedata(m, :NonMarketDamages, :rgdp_per_cap_NonMarketRemainGDP, "test/validationdata/$valdir/rgdp_per_cap_NonMarketRemainGDP.csv")
        savepagedata(m, :SLRDamages, :rcons_per_cap_SLRRemainConsumption, "test/validationdata/$valdir/rcons_per_cap_SLRRemainConsumption.csv")
        savepagedata(m, :SLRDamages, :rgdp_per_cap_SLRRemainGDP, "test/validationdata/$valdir/rgdp_per_cap_SLRRemainGDP.csv")
        savepagedata(m, :SLRDamages, :igdp_ImpactatActualGDPperCapSLR, "test/validationdata/$valdir/igdp_ImpactatActualGDPperCap_sea.csv")
        savepagedata(m, :SLRDamages, :isat_ImpactinclSaturationandAdaptationSLR, "test/validationdata/$valdir/isat_ImpactinclSaturationandAdaptation_SLRise.csv")
        savepagedata(m, :EquityWeighting, :wtct_percap_weightedcosts, "test/validationdata/$valdir/wtct_percap_weightedcosts.csv")
        savepagedata(m, :EquityWeighting, :pct_percap_partiallyweighted, "test/validationdata/$valdir/pct_percap_partiallyweighted.csv")
        savepagedata(m, :EquityWeighting, :dr_discountrate, "test/validationdata/$valdir/dr_discountrate.csv")
        savepagedata(m, :EquityWeighting, :dfc_consumptiondiscountrate, "test/validationdata/$valdir/dfc_consumptiondiscountrate.csv")
        savepagedata(m, :EquityWeighting, :pct_partiallyweighted, "test/validationdata/$valdir/pct_partiallyweighted.csv")
        savepagedata(m, :EquityWeighting, :pcdt_partiallyweighted_discounted, "test/validationdata/$valdir/pcdt_partiallyweighted_discounted.csv")
        savepagedata(m, :EquityWeighting, :wacdt_partiallyweighted_discounted, "test/validationdata/$valdir/wacdt_partiallyweighted_discounted.csv")
        savepagedata(m, :EquityWeighting, :aact_equityweightedadaptation_discountedaggregated, "test/validationdata/$valdir/aact_equityweightedadaptation_discountedaggregated.csv")
        savepagedata(m, :EquityWeighting, :wit_equityweightedimpact, "test/validationdata/$valdir/wit_equityweightedimpact.csv")
        savepagedata(m, :EquityWeighting, :addt_equityweightedimpact_discountedaggregated, "test/validationdata/$valdir/addt_equityweightedimpact_discountedaggregated.csv")
    end

    # climate component
    temp = m[:ClimateTemperature,:rt_g_globaltemperature]
    temp_compare = readpagedata(m, "test/validationdata/$valdir/rt_g_globaltemperature.csv")
    @test temp ≈ temp_compare rtol = 1e-3

    slr = m[:SeaLevelRise,:s_sealevel]
    slr_compare = readpagedata(m, "test/validationdata/$valdir/s_sealevel.csv")
    @test slr ≈ slr_compare rtol = 1e-2

    # Socio-Economics
    gdp = m[:GDP,:gdp]
    gdp_compare = readpagedata(m, "test/validationdata/$valdir/gdp.csv")
    @test gdp ≈ gdp_compare rtol = 1

    pop = m[:Population,:pop_population]
    pop_compare = readpagedata(m, "test/validationdata/$valdir/pop_population.csv")
    @test pop ≈ pop_compare rtol = 0.001

    # Abatement Costs
    abatement = m[:TotalAbatementCosts,:tct_totalcosts]
    abatement_compare = readpagedata(m, "test/validationdata/$valdir/tct_totalcosts.csv")
    @test abatement ≈ abatement_compare rtol = 1e-2

    # Adaptation Costs
    adaptation = m[:TotalAdaptationCosts,:act_adaptationcosts_total]
    adaptation_compare = readpagedata(m, "test/validationdata/$valdir/act_adaptationcosts_tot.csv")
    @test adaptation ≈ adaptation_compare rtol = 1e-2

    # Damages
    damages = m[:Discontinuity,:rcons_per_cap_DiscRemainConsumption]
    damages_compare = readpagedata(m, "test/validationdata/$valdir/rcons_per_cap_DiscRemainConsumption.csv")
    @test damages ≈ damages_compare rtol = 10
    # SLR damages
    slrdamages = m[:SLRDamages,:rcons_per_cap_SLRRemainConsumption]
    slrdamages_compare = readpagedata(m, "test/validationdata/$valdir/rcons_per_cap_SLRRemainConsumption.csv")
    @test slrdamages ≈ slrdamages_compare rtol = 0.1
    # Market damages
    mdamages = m[:MarketDamages,:rcons_per_cap_MarketRemainConsumption]
    mdamages_compare = readpagedata(m, "test/validationdata/$valdir/rcons_per_cap_MarketRemainConsumption.csv")
    @test mdamages ≈ mdamages_compare rtol = 1
    # NonMarket Damages
    nmdamages = m[:NonMarketDamages,:rcons_per_cap_NonMarketRemainConsumption]
    nmdamages_compare = readpagedata(m, "test/validationdata/$valdir/rcons_per_cap_NonMarketRemainConsumption.csv")
    @test nmdamages ≈ nmdamages_compare rtol = 1

    te = m[:EquityWeighting, :te_totaleffect]
    if scenario == "NDCs"
        te_compare = 1.0320923880568126e9
    elseif scenario == "2 degC Target"
        te_compare = 5.0496576104580307e8
    end
    @test te ≈ te_compare rtol = 1e4
end


text:
import Base.size
import Base.isequal
import Base.==
import Base.hash
import Base.show

# MetNet
size(metnet::MetNet) = size(metnet.S)
size(metnet::MetNet, dim) = size(metnet.S, dim)

function ==(metnet1::MetNet, metnet2::MetNet)
    return metnet1.S == metnet2.S && metnet1.lb == metnet2.lb &&
    metnet1.ub == metnet2.ub && metnet1.mets == metnet2.mets &&
    metnet1.rxns == metnet2.rxns && metnet1.rxns == metnet2.rxns
end

isequal(metnet1::MetNet, metnet2::MetNet) = (metnet1 == metnet2)

hash(m::MetNet, h::Int = 0) = hash((:MetNet, m.S, m.b, m.lb, m.ub, h))

show(io::IO, m::MetNet) = summary(io, m)

# MetStates
show(io::IO, s::AbstractMetState) = _print_state_head(io, s)

ext: jl


text:
# type test
F = Graph()
@test isa(fadj(F), Dictionary{Any,Dictionary})
@test isa(badj(F), Dictionary{Any,Dictionary})
@test F.nn == nn(F)
@test F.na == na(F)
@test F.fadj == fadj(F)
@test F.badj == badj(F)

# add node if empty & if exist
nodelistt = [1,2,3,4,5]
for u in eachindex(nodelistt)
    @test nodelistt[u] == add_node!(F, nodelistt[u])
end
for u in eachindex(nodelistt)
    @test nothing == add_node!(F, nodelistt[u])
end

# add nodes from a vector
G = Graph()
add_nodes_from!(G, nodelistt)
@test F == G

# rem node if exist & if empty
for u in eachindex(nodelistt)
    @test rem_node!(G, u)
end
for u in eachindex(nodelistt)
    @test !rem_node!(G, u)
end

# remove nodes from a vector
rem_nodes_from!(F, nodelistt)
@test isempty(fadj(F))

# remove nodes while there's several not in list
add_nodes_from!(F, nodelistt)
res = pop!(nodelistt)
rem_nodes_from!(F, nodelistt)
push!(nodelistt, res)
@test length(fadj(F)) == 1

# add arc
@test add_arc!(F, 1, 2) == 1
@test Nexus.new_arc_key(F, 1, 2) == add_arc!(F, 1, 2)
@test length(fadj(F)[1][2]) == 2
@test add_arc!(F, 1, 2, "haha") == "haha"

# rem arc sampai habis
@test rem_arc!(F, 1, 2)
@test rem_arc!(F, 1, 2)
@test !rem_arc!(F, 1, 2, 2)
@test rem_arc!(F, 1, 2, 1)
@test !rem_arc!(F, 1, 2)

# add arcs from!
arclist = [(1, 2),(2, 3),(3, 4),(4, 5)]
@test add_arcs_from!(G, arclist) == [1,1,1,1]
@test empty!(G) == Graph()
@test add_arcs_from!(G, arclist) == [1,1,1,1]

# rem arcs from!
@test rem_arcs_from!(F, arclist) == 0 # no arc to remove
@test rem_arcs_from!(G, arclist) == length(arclist)


text:
#= BINARY AND DECIMAL PREFIXES:
2^10 ... kilo ... 10^3
2^20 ... mega ... 10^6
2^30 ... giga ... 10^9
...
2^(10*i) ... 10^(3*i)
=#

function simulate(multipliersOfInerest)
    i = 0
    for m in multipliersOfInerest
        while (big(2)^big(10 * i) < m*big(10)^big(3 * i))
            i +=1
        end
        println("Binary 2^$(10*i) is $m times larger than decimal 10^$(3*i) for i = $i")
    end
end

# Usage: simulate(multipliersOfInerest), where multipliersOfInerest ... range, array or something similar.
println("Usage example: simulate(union(2:10, 100, 1000))")

max_stars_repo_name: zznidar/binary-vs-decimal-prefixes

max_issues_repo_name: zznidar/binary-vs-decimal-prefixes

max_forks_repo_name: zznidar/binary-vs-decimal-prefixes


text:
using LinearMaps
using IterativeSolvers
using ProgressMeter

@doc raw"""
Compute the independent-particle susceptibility. Will blow up for large systems.
For non-spin-polarized calculations the matrix dimension is
`prod(basis.fft_size)` × `prod(basis.fft_size)` and
for collinear spin-polarized cases it is
`2prod(basis.fft_size)` × `2prod(basis.fft_size)`.
In this case the matrix has effectively 4 blocks, which are:
```math
\left(\begin{array}{cc}
    (χ_0)_{αα}  & (χ_0)_{αβ} \\
    (χ_0)_{βα}  & (χ_0)_{ββ}
\end{array}\right)
```
"""
function compute_χ0(ham; temperature=ham.basis.model.temperature)
    # We're after χ0(r,r') such that δρ = ∫ χ0(r,r') δV(r') dr'
    # where (up to normalizations)
    # ρ = ∑_nk f(εnk - εF) |ψnk|^2
    # ∑_nk f(εnk - εF) = N_el
    # Everything is summed on k so we omit it for notational simplicity

    # We differentiate wrt a variation δV of the external potential
    # δρ = ∑_n (f'n δεn |ψn|^2 + 2Re fn ψn* δψn - f'n δεF |ψn|^2
    # with fn = f(εnk - εF), f'n = f'(εnk - εF)
    # δN_el = 0 = ∑_n f'n (δεn - δεF)

    # Now we use from first order perturbation theory
    # δεn = <ψn|δV|ψn>
    # δψn = ∑_{m != n} <ψm|δV|ψn> |ψm> / (εn-εm)

    # for δεF we get, with DOS = -∑_n f'n and LDOS = -∑_n f'n |ψn|^2
    # δεF = 1/DOS ∫ δV(r) LDOS(r)dr

    # for δρ we note ρnm = ψn* ψm, and we get
    # δρ = LDOS δεF + ∑_n f'n <ρn|δV> ρn + ∑_{n,m != n} 2Re fn ρnm <ρmn|δV> / (εn-εm)
    # δρ = LDOS δεF + ∑_n f'n <ρn|δV> ρn + ∑_{n,m != n} (fn-fm)/(εn-εm) ρnm <ρnm|δV>
    # The last two terms merge with the convention that (f(x)-f(x))/(x-x) = f'(x) into
    # δρ = LDOS δεF + ∑_{n,m} (fn-fm) ρnm <ρmn|δV> / (εn-εm)
    # Therefore the kernel is LDOS(r) LDOS(r') / DOS + ∑_{n,m} (fn-fm)/(εn-εm) ρnm(r) ρmn(r')
    basis = ham.basis
    model = basis.model
    filled_occ = filled_occupation(model)
    n_spin   = basis.model.n_spin_components
    fft_size = basis.fft_size
    n_fft    = prod(fft_size)

    length(model.symmetries) == 1 || error("Disable symmetries completely for computing χ0")

    EVs = [eigen(Hermitian(Array(Hk))) for Hk in ham.blocks]
    Es = [EV.values for EV in EVs]
    Vs = [EV.vectors for EV in EVs]
    occ, εF = compute_occupation(basis, Es, temperature=temperature)

    χ0 = zeros(eltype(basis), n_spin * n_fft, n_spin * n_fft)
    for (ik, kpt) in enumerate(basis.kpoints)
        # The sum-over-states terms of χ0 are diagonal in the spin blocks (no αβ / βα terms)
        # so the spin of the kpt selects the block we are in
        spinrange = kpt.spin == 1 ? (1:n_fft) : (n_fft+1:2n_fft)
        χ0σσ = @view χ0[spinrange, spinrange]

        N = length(G_vectors(basis, basis.kpoints[ik]))
        @assert N < 10_000
        E = Es[ik]
        V = Vs[ik]
        Vr = cat(G_to_r.(Ref(basis), Ref(kpt), eachcol(V))..., dims=4)
        Vr = reshape(Vr, n_fft, N)
        @showprogress "Computing χ0 for k-point $ik/$(length(basis.kpoints)) ..." for m = 1:N, n = 1:N
            enred = (E[n] - εF) / temperature
            @assert occ[ik][n] ≈ filled_occ * Smearing.occupation(model.smearing, enred)
            ddiff = Smearing.occupation_divided_difference
            ratio = filled_occ * ddiff(model.smearing, E[m], E[n], εF, temperature)
            # dvol because inner products have a dvol in them
            # so that the dual gets one : |f> -> <dvol f|
            # can take the real part here because the nm term is complex conjugate of mn
            # TODO optimize this a bit... use symmetry nm, reduce allocs, etc.
            factor = basis.kweights[ik] * ratio * basis.dvol

            @views χ0σσ .+= factor .* real(conj((Vr[:, m] .* Vr[:, m]'))
                                           .*   (Vr[:, n] .* Vr[:, n]'))
        end
    end
    mpi_sum!(χ0, basis.comm_kpts)

    # Add variation wrt εF (which is not diagonal wrt. spin)
    if temperature > 0
        dos  = compute_dos(εF, basis, Es)
        ldos = compute_ldos(εF, basis, Es, Vs)
        χ0 .+= vec(ldos) .* vec(ldos)' .* basis.dvol ./ sum(dos)
    end
    χ0
end

# make ldiv! act as a given function
struct FunctionPreconditioner
    precondition!  # precondition!(y, x) applies f to x and puts it into y
end
LinearAlgebra.ldiv!(y::T, P::FunctionPreconditioner, x) where {T} = P.precondition!(y, x)::T
LinearAlgebra.ldiv!(P::FunctionPreconditioner, x) = (x .= P.precondition!(similar(x), x))
precondprep!(P::FunctionPreconditioner, ::Any) = P

# Solves (1-P) (H-εn) (1-P) δψn = - (1-P) rhs
# where 1-P is the projector on the orthogonal of ψk
function sternheimer_solver(Hk, ψk, ψnk, εnk, rhs;
                            ψk_extra=zeros(size(ψk,1), 0), εk_extra=zeros(0),
                            tol_cg=1e-12, verbose=false)
    basis = Hk.basis
    kpoint = Hk.kpoint
    model = basis.model
    temperature = model.temperature

    # We use a Schur decomposition of the orthogonal of the occupied states
    # into a part where we have the partially converged, non-occupied bands
    # (which are Rayleigh-Ritz wrt to Hk) and the rest.
    #
    # /!\ This is only implemented for insulators at the moment, WIP for
    # metals, in which case we use empty arrays and all computations are
    # equivalent to invert the Sternheimer system in Ran(1-P) without the Schur
    # trick.

    # TODO
    # finite temperature not supported yet => remove extra bands
    if !iszero(temperature)
        ψk_extra = zeros(size(ψk, 1), 0)
        εk_extra = zeros(0)
    end

    # Projectors:
    # projector onto the computed and converged states
    P(ϕ) = ψk * (ψk' * ϕ)
    # projector onto the computed but nonconverged states
    P_extra(ϕ) = ψk_extra * (ψk_extra' * ϕ)
    # projector onto the computed (converged and unconverged) states
    P_computed(ϕ) = P(ϕ) + P_extra(ϕ)
    # Q = 1-P is the projector onto the orthogonal of converged states
    Q(ϕ) = ϕ - P(ϕ)
    # R = 1-P_computed is the projector onto the orthogonal of computed states
    R(ϕ) = ϕ - P_computed(ϕ)

    # We put things into the form
    # δψkn = ψk_extra * αkn + δψknᴿ ∈ Ran(Q)
    # where δψknᴿ ∈ Ran(R).
    # Note that, if ψk_extra = [], then 1-P = 1-P_computed and
    # δψkn = δψknᴿ is obtained by inverting the full Sternheimer
    # equations in Ran(Q) = Ran(R)
    #
    # This can be summarized as the following:
    #
    # <---- P ----><------------ Q = 1-P -----------------
    #              <-- P_extra -->
    # <--------P_computed -------><-- R = 1-P_computed ---
    # |-----------|--------------|------------------------
    # 1     N_occupied  N_occupied + N_extra

    # ψk_extra are not converged but have been Rayleigh-Ritzed (they are NOT
    # eigenvectors of H) so H(ψk_extra) = ψk_extra' (Hk-εn) ψk_extra should be a
    # real diagonal matrix.
    H(ϕ) = Hk * ϕ - εnk * ϕ
    ψk_exHψk_ex = Diagonal(real.(εk_extra .- εnk))

    # 1) solve for δψknᴿ
    # ----------------------------
    # writing αkn as a function of δψknᴿ, we get that δψknᴿ
    # solves the system (in Ran(1-P_computed))
    #
    # R * (H - εn) * (1 - M * (H - εn)) * R * δψknᴿ = R * (1 - M) * b
    #
    # where M = ψk_extra * (ψk_extra' (H-εn) ψk_extra)^{-1} * ψk_extra'
    # is defined above and b is the projection of -rhs onto Ran(Q).
    #
    b = -Q(rhs)
    bb = R(b -  H(ψk_extra * (ψk_exHψk_ex \ ψk_extra'b)))
    function RAR(ϕ)
        Rϕ = R(ϕ)
        # A denotes here the Schur complement of (1-P) (H-εn) (1-P)
        # with the splitting Ran(1-P) = Ran(P_extra) ⊕ Ran(R)
        ARϕ = Rϕ - ψk_extra * (ψk_exHψk_ex \ ψk_extra'H(Rϕ))
        R(H(ARϕ))
    end
    precon = PreconditionerTPA(basis, kpoint)
    precondprep!(precon, ψnk)
    function R_ldiv!(x, y)
        x .= R(precon \ R(y))
    end
    J = LinearMap{eltype(ψk)}(RAR, size(Hk, 1))
    δψknᴿ = cg(J, bb, Pl=FunctionPreconditioner(R_ldiv!),
               reltol=0, abstol=tol_cg, verbose=verbose)

    # 2) solve for αkn now that we know δψknᴿ
    # Note that αkn is an empty array if there is no extra bands.
    αkn = ψk_exHψk_ex \ ψk_extra' * (b - H(δψknᴿ))

    δψkn = ψk_extra * αkn + δψknᴿ
end

# Apply the four-point polarizability operator χ0_4P = -Ω^-1
# Returns δψ corresponding to a change in *total* Hamiltonian δH
# We start from
# P = f(H-εF), tr(P) = N
# where P is the density matrix, f the occupation function
# δεn = <ψn|δV|ψn>
# 0 = ∑_n fn' (δεn - δεF) determines δεF
# where fn' = f'((εn-εF)/T)/T

# Then <ψm|δP|ψn> = (fm-fn)/(εm-εn) <ψm|δH|ψn>
# Except for the diagonal which is
# <ψn|δP|ψn> = (fn'-δεF) <ψn|δH|ψn>

# We want to represent this with a δψ. We do *not* impose that
# δψ is orthogonal at finite temperature.
# We get
# δP = ∑_k fk (|δψk><ψk| + |ψk><δψk|)
# Identifying with <ψm|δP|ψn> we get for the diagonal terms
# <ψn|δψn> fn = fn'<ψn|δH-δεF|ψn>
# and for the off-diagonal
# (fm-fn)/(εm-εn) <ψm|δH|ψn> = fm <δψm|ψn> + fn <ψm|δψn>

# We split the computation of δψn in two contributions:
# for the already-computed states, we add an explicit contribution
# for the empty states, we solve a Sternheimer equation
# (H-εn) δψn = - P_{ψ^⟂} δH ψn

# The off-diagonal explicit term needs a careful consideration of stability.
# Let <ψm|δψn> = αmn <ψm|δH|ψn>. αmn has to satisfy
# fn αmn + fm αnm = ratio = (fn-fm)/(εn-εm)   (*)
# The usual way is to impose orthogonality (=> αmn=-αnm),
# but this means that αmn = 1/(εm-εn), which is unstable
# Instead, we minimize αmn^2 + αnm^2 under the linear constraint (*), which leads to
# αmn = ratio * fn / (fn^2 + fm^2)
# fn αmn = ratio * fn^2 / (fn^2 + fm^2)

# This formula is very nice
# - It gives a vanishing contribution fn αmn for empty states
#   (note that α itself blows up, but it's compensated by fn)
# - In the case where fn=1/0 or fm=0 we recover the same formulas
#   as the ones with orthogonality
# - When n=m it gives the correct contribution
# - It does not blow up for degenerate states
function compute_αmn(fm, fn, ratio)
    ratio == 0 && return ratio
    ratio * fn / (fn^2 + fm^2)
end

@views @timing function apply_χ0_4P(ham, ψ, occ, εF, eigenvalues, δHψ;
                                    ψ_extra=[zeros(size(ψk, 1), 0) for ψk in ψ],
                                    ε_extra=[zeros(0) for _ in eigenvalues],
                                    kwargs_sternheimer...)
    basis  = ham.basis
    model = basis.model
    temperature = model.temperature
    filled_occ = filled_occupation(model)
    T      = eltype(basis)
    Nk = length(basis.kpoints)

    # First compute δεF
    δεF = zero(T)
    δocc = [zero(occ[ik]) for ik = 1:Nk]  # = fn' * (δεn - δεF)
    if temperature > 0
        # First compute δocc without self-consistent Fermi δεF
        D = zero(T)
        for ik = 1:Nk
            for (n, εnk) in enumerate(eigenvalues[ik])
                enred = (εnk - εF) / temperature
                δεnk = real(dot(ψ[ik][:, n], δHψ[ik][:, n]))
                fpnk = (filled_occ
                        * Smearing.occupation_derivative(model.smearing, enred)
                        / temperature)
                δocc[ik][n] = δεnk * fpnk
                D += fpnk * basis.kweights[ik]
            end
        end
        # compute δεF
        D = mpi_sum(D, basis.comm_kpts)  # equal to minus the total DOS
        δocc_tot = mpi_sum(sum(basis.kweights .* sum.(δocc)), basis.comm_kpts)
        δεF = δocc_tot / D
    end

    # compute δψnk band per band
    δψ = zero.(ψ)
    for ik = 1:Nk
        ψk = ψ[ik]
        δψk = δψ[ik]

        εk = eigenvalues[ik]
        for n = 1:length(εk)
            fnk = filled_occ * Smearing.occupation(model.smearing, (εk[n]-εF) / temperature)

            # explicit contributions (nonzero only for temperature > 0)
            for m = 1:length(εk)
                fmk = filled_occ * Smearing.occupation(model.smearing, (εk[m]-εF) / temperature)
                ddiff = Smearing.occupation_divided_difference
                ratio = filled_occ * ddiff(model.smearing, εk[m], εk[n], εF, temperature)
                αmn = compute_αmn(fmk, fnk, ratio)  # fnk * αmn + fmk * αnm = ratio
                δψk[:, n] .+= ψk[:, m] .* αmn .* (dot(ψk[:, m], δHψ[ik][:, n]) - (n == m) * δεF)
            end

            # Sternheimer contribution
            # TODO unoccupied orbitals may have a very bad sternheimer solve, be careful about this
            δψk[:, n] .+= sternheimer_solver(ham.blocks[ik], ψk, ψk[:, n], εk[n], δHψ[ik][:, n];
                                             ψk_extra=ψ_extra[ik], εk_extra=ε_extra[ik],
                                             kwargs_sternheimer...)
        end
    end
    δψ
end

"""
Get the density variation δρ corresponding to a total potential variation δV.

Note: This function assumes that all bands contained in `ψ` and `eigenvalues` are
sufficiently converged. By default the `self_consistent_field` routine of `DFTK`
returns `3` extra bands, which are not converged by the eigensolver
(see `n_ep_extra` parameter). These should be discarded from `ψ` before using this
function. However, one can still use additional information from nonconverged bands
stored in `ψ_extra` (typically, the `3` extra bands returned by default in SCF
routine).
"""
function apply_χ0(ham, ψ, εF, eigenvalues, δV;
                  ψ_extra=[zeros(size(ψk, 1), 0) for ψk in ψ],
                  ε_extra=[zeros(0) for _ in eigenvalues],
                  kwargs_sternheimer...)

    basis = ham.basis
    model = basis.model
    occ = [filled_occupation(model) *
           Smearing.occupation.(model.smearing, (eigenvalues[ik] .- εF) ./ model.temperature)
           for ik = 1:length(basis.kpoints)]

    # Make δV respect the basis symmetry group, since we won't be able
    # to compute perturbations that don't anyway
    δV = symmetrize_ρ(basis, δV)

    # Normalize δV to avoid numerical trouble; theoretically should
    # not be necessary, but it simplifies the interaction with the
    # Sternheimer linear solver (it makes the rhs be order 1 even if
    # δV is small)
    normδV = norm(δV)
    normδV < eps(typeof(εF)) && return zero(δV)
    δV ./= normδV

    δHψ = [DFTK.RealSpaceMultiplication(basis, kpt, @views δV[:, :, :, kpt.spin]) * ψ[ik]
           for (ik, kpt) in enumerate(basis.kpoints)]
    δψ = apply_χ0_4P(ham, ψ, occ, εF, eigenvalues, δHψ;
                     ψ_extra, ε_extra, kwargs_sternheimer...)
    δρ = DFTK.compute_δρ(basis, ψ, δψ, occ)
    δρ * normδV
end

"""
We use here the full scfres, with a distinction between converged bands and
nonconverged extra bands used in the SCF.
"""
function apply_χ0(scfres, δV; kwargs_sternheimer...)
    n_ep_extra = scfres.n_ep_extra

    ψ       = [@view ψk[:, 1:end-n_ep_extra]     for ψk in scfres.ψ]
    ψ_extra = [@view ψk[:, end-n_ep_extra+1:end] for ψk in scfres.ψ]

    ε       = [εk[1:end-n_ep_extra]     for εk in scfres.eigenvalues]
    ε_extra = [εk[end-n_ep_extra+1:end] for εk in scfres.eigenvalues]

    apply_χ0(scfres.ham, ψ, scfres.εF, ε, δV; ψ_extra, ε_extra,
             kwargs_sternheimer...)
end


text:
module Jvuetify
include("utilfuns.jl")
export jhtml, htmltable, jggplot, plotlygraph
end # module


text:
# block is for if there are multiple blocks of parameters being produce in succession
# then it is meant to indicate the ordering of blocks for visual comparison.
# I.e. if param_range 1:10 corresponds to the first block, and 11:20 corresponds to the
# second block then the files will be aptly named ...block=1.pdf, ...block=2.pdf
function plot_posterior_intervals(m::AbstractDSGEModel; cloud::ParticleCloud = ParticleCloud(m, 0),
                                  df::DataFrame = DataFrame(),
                                  plotroot::String = figurespath(m, "estimate"),
                                  verbose::Symbol = :low, label::String = "", title::String
                                  = isempty(label) ? "" : "$label Posterior Intervals",
                                  block::Int64 = 0, include_fixed::Bool = false,
                                  param_range::UnitRange = include_fixed ? UnitRange(1, n_parameters(m)) : UnitRange(1, n_parameters_free(m)),
                                  excl_list::Vector{Symbol} = Vector{Symbol}(undef, 0))

    # If the moments to be plotted are not pre-provided
    if isempty(df)
        df = load_posterior_moments(m; cloud = cloud, load_bands = true, include_fixed =
                                    include_fixed, excl_list = excl_list)
    end

    # Indexing out a block of the parameters
    df = df[param_range, :]
    df[:param_inds] = param_range
    parameter_labels = df[:param]

    # Plot non-outliers
    p = @df df plot(:param_inds, :post_ub, label = "", linealpha = 0, marker =
                    :hline, markerstrokecolor = :black, title = title, series_annotations =
                    parameter_labels)
    @df df plot!(p, :param_inds, :post_lb, label = "", linealpha = 0, marker =
                 :hline, markerstrokecolor = :black)
    @df df plot!(:param_inds, :post_mean, label = "", linealpha = 0, marker =
                 :hline, markerstrokecolor = :black)

    rescale_annotations_fontsize!(p, 6)
    relocate_annotations!(p, :right, :top)

    if !isempty(plotroot)
        vint = data_vintage(m)
        output_file = joinpath(plotroot, "posterior_intervals_" * "vint=$vint")
        if !iszero(block)
            output_file = output_file * "_block=$block"
        end
        output_file = output_file * ".pdf"
        save_plot(p, output_file, verbose = verbose)
    else
        return p
    end
end

# Calculate the posterior intervals averaged over the estimations in the
# clouds vector
function plot_posterior_intervals(m::AbstractDSGEModel, clouds::Vector{ParticleCloud};
                                  plotroot::String = figurespath(m, "estimate"),
                                  verbose::Symbol = :low, label::String = "",
                                  title::String = isempty(label) ? "" : "$label Posterior Intervals",
                                  block::Int64 = 0, include_fixed::Bool = false,
                                  param_range::UnitRange = include_fixed ? UnitRange(1, n_parameters(m)) : UnitRange(1, n_parameters_free(m)),
                                  excl_list::Vector{Symbol} = Vector{Symbol}(undef, 0))

    df = load_posterior_moments(m, clouds, load_bands = true,
                                include_fixed = include_fixed, excl_list = excl_list)

    plot_posterior_intervals(m; df = df, plotroot = plotroot, verbose = verbose,
                             label = label, title = title, param_range = param_range,
                             block = block, include_fixed = include_fixed, excl_list = excl_list)
end

# block is for if there are multiple blocks of parameters being produce in succession
# then it is meant to indicate the ordering of blocks for visual comparison.
# I.e. if param_range 1:35 corresponds to the first block, and 36:68 corresponds to the
# second block then the files will be aptly named ...b1.pdf, ...b2.pdf
function plot_posterior_interval_comparison(m_baseline::AbstractDSGEModel,
                                            m_comparison::AbstractDSGEModel;
                                            cloud_baseline::ParticleCloud = ParticleCloud(m_baseline, 0),
                                            cloud_comparisons::Vector{ParticleCloud} = [ParticleCloud(m_comparison, 0)],
                                            df_baseline::DataFrame = DataFrame(),
                                            df_comparisons::Vector{DataFrame} = [DataFrame()],
                                            plotroot::String = figurespath(m_baseline, "estimate"),
                                            verbose::Symbol = :low,
                                            in_deviations::Bool = false,
                                            scale_by_std::Bool = false,
                                            base_label::String = "",
                                            comp_labels::Vector{String} = [""],
                                            title::String = in_deviations ? "$comp_label deviations from $base_label Interval Comparisons" : "$base_label and $comp_label Interval Comparisons",
                                            block::Int64 = 0, include_fixed::Bool = false,
                                            param_range::UnitRange = include_fixed ? UnitRange(1, n_parameters(m_baseline)) : UnitRange(1, n_parameters_free(m_baseline)),
                                            excl_list::Vector{Symbol} = Vector{Symbol}(undef, 0),
                                            filename_tag::String = "")

    if isempty(df_baseline)
        df_baseline = load_posterior_moments(m_baseline; cloud = cloud_baseline, load_bands = true,
                                             include_fixed = include_fixed, excl_list = excl_list)
    end

    if isempty(df_comparisons[1])
        for i=1:length(cloud_comparisons)
            df_comparisons[i] = load_posterior_moments(m_comparison; cloud = cloud_comparisons[i],
                                                      load_bands = true, include_fixed = include_fixed,
                                                      excl_list = excl_list)
        end
    end
    for i=1:length(df_comparisons)
        df_comparisons[i] = df_comparisons[i][param_range, :]
        df_comparisons[i][:param_inds] = param_range
    end

    # Indexing out a block of the parameters
    df_baseline = df_baseline[param_range, :]
  #  df_comparison = df_comparison[param_range, :]
    df_baseline[:param_inds] = param_range
  #  df_comparison[:param_inds] = param_range

    parameter_labels = df_baseline[:param]

    if in_deviations
        df_deviations = DataFrame()
        df_deviations[:param] = df_baseline[:param]
        df_deviations[:param_inds] = df_baseline[:param_inds]

        # Baseline model mean normalized to 0
        df_deviations[:base_mean] = zeros(size(df_baseline, 1))

        # Comparison model mean and bands normalized to deviations from baseline mean and
        # bands (hence the mean deviation may be higher than either of the bands' deviations
        # in the resulting plot)
        df_deviations[:comp_mean] = df_comparison[:post_mean] - df_baseline[:post_mean]
        df_deviations[:comp_ub] = df_comparison[:post_ub] - df_baseline[:post_ub]
        df_deviations[:comp_lb] = df_comparison[:post_lb] - df_baseline[:post_lb]

        # Scale the comparison mean and lower and upper bounds deviations by the std. of the
        # parameter under the baseline
        if scale_by_std
            df_deviations[:comp_mean] = df_deviations[:comp_mean]./df_baseline[:post_std]
            df_deviations[:comp_ub] = df_deviations[:comp_ub]./df_baseline[:post_std]
            df_deviations[:comp_lb] = df_deviations[:comp_lb]./df_baseline[:post_std]
        end

        p = @df df_deviations plot(:param_inds, :comp_ub, label = "", linealpha = 0, marker =
                                :hline, markerstrokecolor = :blue, title = "",
                                series_annotations = parameter_labels, legend=:bottomright)
        @df df_deviations plot!(p, :param_inds, :comp_lb, label = "", linealpha = 0, marker =
                                :hline, markerstrokecolor = :blue, legend=:bottomright)
        @df df_deviations plot!(p, :param_inds, :base_mean, label = "", marker = :hline,
                                markerstrokecolor = :black, linecolor = :black, legend=:bottomright)
        @df df_deviations plot!(p, :param_inds, :comp_mean, label = "", linealpha = 0, marker =
                                :hline, markerstrokecolor = :red, legend=:bottomright)

        rescale_annotations_fontsize!(p, 6)
        relocate_annotations!(p, :right, :top)
    else
        p = @df df_baseline plot(:param_inds, :post_ub, label = base_label, linealpha = 0, marker =
                              :hline, markerstrokecolor = :black, title = "", series_annotations =
                              parameter_labels, legend=:bottomright)
        @df df_baseline plot!(p, :param_inds, :post_lb, label = "", linealpha = 0, marker =
                              :hline, markerstrokecolor = :black, legend=:bottomright)
        @df df_baseline plot!(p, :param_inds, :post_mean, label = "", linealpha = 0, marker =
                              :hline, markerstrokecolor = :black, legend=:bottomright)
        colors = [:red, :blue, :green, :orange, :yellow]
        for i = 1:length(df_comparisons)
            @df (df_comparisons[i]) plot!(p, :param_inds, :post_ub, label = comp_labels[i], linealpha = 0, marker =
                                    :hline, markerstrokecolor = colors[i], legend=:bottomright)
            @df (df_comparisons[i]) plot!(p, :param_inds, :post_lb, label = "", linealpha = 0, marker =
                                    :hline, markerstrokecolor = colors[i], legend=:bottomright)
            @df (df_comparisons[i]) plot!(p, :param_inds, :post_mean, label = "", linealpha = 0, marker =
                                    :hline, markerstrokecolor = colors[i], legend=:bottomright)
        end
        rescale_annotations_fontsize!(p, 6)
        relocate_annotations!(p, :right, :top)
    end

    if !isempty(plotroot)
        vint = data_vintage(m_baseline)
        output_file = joinpath(plotroot, "posterior_comparison_intervals_" *
                               "deviat=$(in_deviations)_" * "vint=$vint" * "_" *filename_tag)
        if !iszero(block)
            output_file = output_file * "_block=$block"
        end
        output_file = output_file * ".pdf"
        save_plot(p, output_file, verbose = verbose)
    else
        return p
    end
end

# Calculate the posterior intervals averaged over the estimations in the
# clouds vectors
function plot_posterior_interval_comparison(m_baseline::AbstractDSGEModel,
                                            m_comparison::AbstractDSGEModel,
                                            clouds_baseline::Vector{ParticleCloud},
                                            clouds_comparisons::Vector{Vector{ParticleCloud}};
                                            plotroot::String = figurespath(m_baseline, "estimate"),
                                            verbose::Symbol = :low,
                                            in_deviations::Bool = false,
                                            scale_by_std::Bool = false,
                                            base_label::String = "",
                                            comp_labels::Vector{String} = [""],
                                            title::String = in_deviations ? "$(comp_labels[1]) deviations from $base_label Interval Comparisons" : "$base_label and $(comp_labels[1]) Interval Comparisons",
                                            block::Int64 = 0, include_fixed::Bool = false,
                                            param_range::UnitRange = include_fixed ? UnitRange(1, n_parameters(m_baseline)) : UnitRange(1, n_parameters_free(m_baseline)),
                                            excl_list::Vector{Symbol} = Vector{Symbol}(undef, 0),
                                            filename_tag::String = "")

    df_baseline = load_posterior_moments(m_baseline, clouds_baseline, load_bands = true,
                                         include_fixed = include_fixed, excl_list = excl_list)
    df_comparisons = Vector{DataFrame}(undef, length(clouds_comparisons))
    for i=1:length(clouds_comparisons)
        df_comparisons[i] = load_posterior_moments(m_comparison, clouds_comparisons[i], load_bands = true,
                                         include_fixed = include_fixed, excl_list = excl_list)
    end

    plot_posterior_interval_comparison(m_baseline, m_comparison, df_baseline = df_baseline, df_comparisons = df_comparisons,
                                       plotroot = plotroot, verbose = verbose, in_deviations = in_deviations, scale_by_std = scale_by_std,
                                       base_label = base_label, comp_labels = comp_labels, title = "", param_range = param_range, block =
                                       block, excl_list = excl_list, filename_tag = filename_tag)
end


text:
module ReactiveMPTest

## https://discourse.julialang.org/t/generation-of-documentation-fails-qt-qpa-xcb-could-not-connect-to-display/60988
## https://gr-framework.org/workstations.html#no-output
ENV["GKSwstype"] = "100"

using Test, Documenter, ReactiveMP, Distributions
using TestSetExtensions
using Aqua

# DocMeta.setdocmeta!(ReactiveMP, :DocTestSetup, :(using ReactiveMP, Distributions); recursive=true)

# Example usage of a reduced testset
# julia --project --color=yes -e 'import Pkg; Pkg.test(test_args = [ "distributions:normal_mean_variance" ])'

function addtests(filename)
    key = filename_to_key(filename)
    if isempty(enabled_tests) || key in enabled_tests
        include(filename)
    end
end

function key_to_filename(key)
    splitted = split(key, ":")
    return length(splitted) === 1 ? string("test_", first(splitted), ".jl") : string(join(splitted[1:end - 1], "/"), "/test_", splitted[end], ".jl")
end

function filename_to_key(filename)
    splitted   = split(filename, "/")
    if length(splitted) === 1
        return replace(replace(first(splitted), ".jl" => ""), "test_" => "")
    else
        path, name = splitted[1:end - 1], splitted[end]
        return string(join(path, ":"), ":", replace(replace(name, ".jl" => ""), "test_" => ""))
    end
end

enabled_tests = lowercase.(ARGS)

if isempty(enabled_tests)
    println("Running all tests...")
    Aqua.test_all(ReactiveMP; ambiguities=false)
    # doctest(ReactiveMP)
else
    println("Running specific tests: $enabled_tests")
end

@testset ExtendedTestSet "ReactiveMP" begin

    function key_to_filename(key)
        splitted = split(key, ":")
        return length(splitted) === 1 ? string("test_", first(splitted), ".jl") : string(join(splitted[1:end - 1], "/"), "/test_", splitted[end], ".jl")
    end

    function filename_to_key(filename)
        splitted   = split(filename, "/")
        if length(splitted) === 1
            return replace(replace(first(splitted), ".jl" => ""), "test_" => "")
        else
            path, name = splitted[1:end - 1], splitted[end]
            return string(join(path, ":"), ":", replace(replace(name, ".jl" => ""), "test_" => ""))
        end
    end

    function addtests(filename)
        key = filename_to_key(filename)
        if isempty(enabled_tests) || key in enabled_tests
            include(filename)
        end
    end

    @testset "Testset helpers" begin
        @test key_to_filename(filename_to_key("distributions/test_normal_mean_variance.jl")) == "distributions/test_normal_mean_variance.jl"
        @test filename_to_key(key_to_filename("distributions:normal_mean_variance")) == "distributions:normal_mean_variance"
        @test key_to_filename(filename_to_key("test_message.jl")) == "test_message.jl"
        @test filename_to_key(key_to_filename("message")) == "message"
    end

    addtests("algebra/test_correction.jl")
    addtests("algebra/test_permutation_matrix.jl")

    addtests("test_math.jl")
    addtests("test_helpers.jl")

    addtests("constraints/prod/test_prod_final.jl")

    addtests("test_distributions.jl")
    addtests("distributions/test_common.jl")
    addtests("distributions/test_bernoulli.jl")
    addtests("distributions/test_beta.jl")
    addtests("distributions/test_categorical.jl")
    addtests("distributions/test_contingency.jl")
    addtests("distributions/test_exp_linear_quadratic.jl")
    addtests("distributions/test_dirichlet_matrix.jl")
    addtests("distributions/test_dirichlet.jl")
    addtests("distributions/test_gamma.jl")
    addtests("distributions/test_mv_normal_mean_covariance.jl")
    addtests("distributions/test_mv_normal_mean_precision.jl")
    addtests("distributions/test_mv_normal_weighted_mean_precision.jl")
    addtests("distributions/test_normal_mean_variance.jl")
    addtests("distributions/test_normal_mean_precision.jl")
    addtests("distributions/test_normal_weighted_mean_precision.jl")
    addtests("distributions/test_normal.jl")
    addtests("distributions/test_pointmass.jl")
    addtests("distributions/test_wishart.jl")
    addtests("distributions/test_sample_list.jl")

    addtests("test_message.jl")

    addtests("test_variable.jl")

    addtests("test_node.jl")
    addtests("nodes/flow/test_flow.jl")
    addtests("nodes/test_addition.jl")
    addtests("nodes/test_bifm.jl")
    addtests("nodes/test_bifm_helper.jl")
    addtests("nodes/test_subtraction.jl")
    addtests("nodes/test_probit.jl")
    addtests("nodes/test_autoregressive.jl")

    addtests("rules/flow/test_marginals.jl")
    addtests("rules/flow/test_in.jl")
    addtests("rules/flow/test_out.jl")

    addtests("rules/addition/test_marginals.jl")
    addtests("rules/addition/test_in1.jl")
    addtests("rules/addition/test_in2.jl")
    addtests("rules/addition/test_out.jl")

    addtests("rules/bifm/test_marginals.jl")
    addtests("rules/bifm/test_in.jl")
    addtests("rules/bifm/test_out.jl")
    addtests("rules/bifm/test_zprev.jl")
    addtests("rules/bifm/test_znext.jl")

    addtests("rules/bifm_helper/test_in.jl")
    addtests("rules/bifm_helper/test_out.jl")

    addtests("rules/normal_mixture/test_out.jl")

    addtests("rules/subtraction/test_marginals.jl")
    addtests("rules/subtraction/test_in1.jl")
    addtests("rules/subtraction/test_in2.jl")
    addtests("rules/subtraction/test_out.jl")

    addtests("rules/bernoulli/test_out.jl")
    addtests("rules/bernoulli/test_p.jl")
    addtests("rules/bernoulli/test_marginals.jl")

    addtests("rules/beta/test_out.jl")
    addtests("rules/beta/test_marginals.jl")

    addtests("rules/dot_product/test_out.jl")
    addtests("rules/dot_product/test_in1.jl")
    addtests("rules/dot_product/test_in2.jl")
    addtests("rules/dot_product/test_marginals.jl")

    addtests("rules/normal_mean_variance/test_out.jl")
    addtests("rules/normal_mean_variance/test_mean.jl")
    addtests("rules/normal_mean_variance/test_var.jl")

    addtests("rules/normal_mean_precision/test_out.jl")
    addtests("rules/normal_mean_precision/test_mean.jl")
    addtests("rules/normal_mean_precision/test_precision.jl")

    addtests("rules/mv_normal_mean_covariance/test_out.jl")
    addtests("rules/mv_normal_mean_covariance/test_mean.jl")

    addtests("rules/mv_normal_mean_precision/test_out.jl")
    addtests("rules/mv_normal_mean_precision/test_mean.jl")
    addtests("rules/mv_normal_mean_precision/test_precision.jl")

    addtests("rules/probit/test_out.jl")
    addtests("rules/probit/test_in.jl")

    addtests("rules/wishart/test_marginals.jl")
    addtests("rules/wishart/test_out.jl")

    addtests("models/test_lgssm.jl")
    addtests("models/test_hgf.jl")
    addtests("models/test_ar.jl")
    addtests("models/test_gmm.jl")
    addtests("models/test_hmm.jl")

end

end


text:
# This file was generated, do not modify it. # hide
xgb.num_round = 200;


text:
function compute_RMSE(U, V, Y, r, d1, d2, rows_t, vals_t, cols_t)
	res = 0.0
	n = 0.0
	for i = 1:d1
		tmp = nzrange(Y, i)
		d2_bar = rows_t[tmp];
		vals_d2_bar = vals_t[tmp];
		ui = U[:, i]
		len = size(d2_bar)[1]
		for j = 1:len
			J = d2_bar[j];
			vj = V[:, J]
			res += (vals_d2_bar[j] - dot(ui,vj))^2
			n += 1.0
		end
	end
	return (res / n)^0.5
end

function compute_RMSE_train(U, V, X, r, d1, d2, rows, vals, cols)
	res = 0.0
	n = 0.0
	for i = 1:d1
		tmp = nzrange(X, i)
		d2_bar = rows[tmp];
		vals_d2_bar = vals[tmp];
		ui = U[:, i]
		len = size(d2_bar)[1]
		for j = 1:len
			J = d2_bar[j];
			vj = V[:, J]
			res += (vals_d2_bar[j] - dot(ui,vj))^2
			n += 1.0
		end
	end
	return (res / n)^0.5
end

function compute_RMSE(U, V, Y, r, d1, d2, rows_t, vals_t, cols_t, threshold)
    res = 0.0
    n = 0.0

    u_bar = zeros(r)
    for i = 1:d1
        u_bar += U[:, i]
    end
    u_bar /= d1

    for i = 1:d1
        tmp = nzrange(Y, i)
        d2_bar = rows_t[tmp];
        vals_d2_bar = vals_t[tmp];
        #ui = U[:, i]
        #ui = U[:, i] * threshold + (1 - threshold) * u_bar
        if rand(1)[1] > threshold
            i = rand(1:d1)
        end
        ui = U[:, i]
        len = size(d2_bar)[1]
        for j = 1:len
            J = d2_bar[j];
            vj = V[:, J]
            res += (vals_d2_bar[j] - dot(ui,vj))^2
            n += 1.0
        end
    end
    return (res / n)^0.5
end

function compute_RMSE_train(U, V, X, r, d1, d2, rows, vals, cols, threshold)
    res = 0.0
    n = 0.0
    u_bar = zeros(r)
    for i = 1:d1
        u_bar += U[:, i]
    end
    u_bar /= d1
    for i = 1:d1
        tmp = nzrange(X, i)
        d2_bar = rows[tmp];
        vals_d2_bar = vals[tmp];
        #ui = U[:, i]
        #ui = U[:, i] * threshold + (1 - threshold) * u_bar
        if rand(1)[1] > threshold
            i = rand(1:d1)
        end
        ui = U[:, i]
        len = size(d2_bar)[1]
        for j = 1:len
            J = d2_bar[j];
            vj = V[:, J]
            res += (vals_d2_bar[j] - dot(ui,vj))^2
            n += 1.0
        end
    end
    return (res / n)^0.5
end

function objective(U, V, X, d1, lambda, rows, vals)
	res = 0.0
	res = lambda * (vecnorm(U) ^ 2 + vecnorm(V) ^ 2)
	for i in 1:d1
		tmp = nzrange(X, i)
		d2_bar = rows[tmp];
		vals_d2_bar = vals[tmp];
		len = size(d2_bar)[1];
		for k in 1:len
			j = d2_bar[k]
			res += (vals_d2_bar[k] - dot(U[:,i], V[:,j]))^2
		end
	end
	return res
end

function compute_pairwise_error_ndcg(U, V, Y, r, d1, d2, rows_t, vals_t, cols_t, ndcg_k)
	sum_error = 0.; ndcg_sum = 0.;
	for i = 1:d1
		tmp = nzrange(Y, i)
		d2_bar = rows_t[tmp];
		vals_d2_bar = vals_t[tmp];
		ui = U[:, i]
		len = size(d2_bar)[1]
		score = zeros(len)
		for j = 1:len
			J = d2_bar[j];
			vj = V[:, J]
			score[j] = dot(ui,vj)
		end
		error_this = 0
		n_comps_this = 0
		for j in 1:(len - 1)
			jval = vals_d2_bar[j]
			for k in (j + 1):len
				kval = vals_d2_bar[k]
				if score[j] >= score[k] && jval < kval
					error_this += 1
				end
				if score[j] <= score[k] && jval > kval
					error_this += 1
				end
				n_comps_this += 1
			end
		end
		sum_error += error_this / n_comps_this

		p1 = sortperm(score, rev = true)
		p1 = p1[1:ndcg_k]
		M1 = vals_d2_bar[p1]
		p2 = sortperm(vals_d2_bar, rev = true)
		p2 = p2[1:ndcg_k]
		M2 = vals_d2_bar[p2]
		dcg = 0.; dcg_max = 0.
		for k = 1:ndcg_k
			dcg += (2 ^ M1[k] - 1) / log2(k + 1)
			dcg_max += (2 ^ M2[k] - 1) / log2(k + 1)
		end
		ndcg_sum += dcg / dcg_max
	end
	return sum_error / d1, ndcg_sum / d1
end

function update(U, V, X, r, d1, d2, lambda, rows, vals, stepsize, cols, threshold, p)
	l = rand(1:length(rows))
	j = rows[l]
	i = cols[l]
	#i = rand(1:d1)
	#j = rand(1:d2)
	#eij = X[j, i] - dot(U[:,i], V[:,j])

	if rand(1)[1] > threshold
		i = rand(1:d1)
	end
	if rand(1)[1] > threshold
		j = rand(1:d2)
	end

	eij = vals[l] - dot(U[:,i], V[:,j])

	for k in 1:r
		if rand(1)[1] <= p
			U[k, i] +=  stepsize * (eij * V[k,j] - lambda * U[k,i]) / p
		end
	end
	for k in 1:r
		if rand(1)[1] <= p
			V[k,j] += stepsize * (eij * U[k,i] - lambda * V[k,j]) / p
		end
	end
	#U[:,i] = ui
	#V[:,j] = vj
	return U, V
end

# command to run julia program after include this file
# main("data/ml1m_train_ratings.csv", "data/ml1m_test_ratings.csv", 100, 0.1, 0.99)
# p is probability of keeping embeddings

function main(train, test, r, lambda, threshold, p)
	X = readdlm(train, ',' , Int64);
	x = vec(X[:,1]);
	y = vec(X[:,2]);
	v = vec(X[:,3]);
	Y = readdlm(test, ',' , Int64);
	xx = vec(Y[:,1]);
	yy = vec(Y[:,2]);
	vv = vec(Y[:,3]);
	# userid; movieid
	# n = 6040; msize = 3952;
	# depending on the size of X, read n_users and n_items from python output
	n = max(maximum(x), maximum(xx)); msize = max(maximum(y), maximum(yy));
	println("Training dataset ", train, " and test dataset ", test, " are loaded. \n There are ", n, " users and ", msize, " items in the dataset.")
	#n = 1496; msize = 3952;
	#n = 12851; msize = 65134
	#n = 221004; msize = 17771
	X = sparse(x, y, v, n, msize); # userid by movieid
	Y = sparse(xx, yy, vv, n, msize);
	# julia column major
	# now moveid by userid
	X = X';
	Y = Y';

	# too large to debug the algorithm, subset a small set: 500 by 750
	#X = X[1:500, 1:750];
	#X = X[1:1000, 1:2000];
	rows = rowvals(X);
	vals = nonzeros(X);
	cols = zeros(Int, size(vals)[1]);

	d2, d1 = size(X);
	cc = 0;
	for i = 1:d1
		tmp = nzrange(X, i);
		nowlen = size(tmp)[1];
		for j = 1:nowlen
			cc += 1
			cols[cc] = i
		end
	end

	rows_t = rowvals(Y);
	vals_t = nonzeros(Y);
	cols_t = zeros(Int, size(vals_t)[1]);
	cc = 0;
	for i = 1:d1
		tmp = nzrange(Y, i);
		nowlen = size(tmp)[1];
		for j = 1:nowlen
			cc += 1
			cols_t[cc] = i
		end
	end

	#r = 100;
	#lambda = 0.1
	#lambda = 7000 # works better for movielens10m data
	#lambda = 10000; # works better for netflix data
	ndcg_k = 10;
	# initialize U, V

	U = randn(r, d1); V = randn(r, d2);
	# U = 0.01*randn(r, d1); V = 0.01*randn(r, d2); # works better for netflix data

	#stepsize = 0.0001
	stepsize = 0.005
        totaltime = 0.00000;
	#println("iter time objective_function pairwise_error NDCG RMSE")
	println("iter time objective_function train_RMSE test_RMSE")

	#pairwise_error, ndcg = compute_pairwise_error_ndcg(U, V, Y, r, d1, d2, rows_t, vals_t, cols_t, ndcg_k)

	nowobj = objective(U, V, X, d1, lambda, rows, vals)
	rmse = compute_RMSE(U, V, Y, r, d1, d2, rows_t, vals_t, cols_t)
	rmse_tr = compute_RMSE_train(U, V, X, r, d1, d2, rows, vals, cols)
	#println("[", 0, ", ", totaltime, ", ", nowobj, ", ", pairwise_error, ", ", ndcg, ", ", rmse, ", ", rmse_tr, "],")
	println("[", 0, ", ", totaltime, ", ", nowobj, ", ", rmse, ", ", rmse_tr, "],")

	#for iter in 1:150000000
	for iter in 1:3000000000
		tic();

 		#V, m, nowobj  = update_V(U, V, X, r, d1, d2, lambda, rows, vals, stepsize, cols)

 		#U, nowobj = update_U(U, V, X, r, d1, d2, lambda, rows, vals, stepsize, m)

		U, V = update(U, V, X, r, d1, d2, lambda, rows, vals, stepsize, cols, threshold, p)

		totaltime += toq();

		# need to add codes for computing pairwise error and NDCG

		#pairwise_error = compute_pairwise_error(U, V, Y, r, d1, d2, rows_t, vals_t, cols_t)
		#ndcg = compute_NDCG(U, V, Y, r, d1, d2, rows_t, vals_t, cols_t, ndcg_k)
	 	#if iter % 3000000 == 0
	 	if iter % 10000000 == 0
        	#if iter % 100000000 == 0
            		#threshold += (1 - threshold) / 3
	 		nowobj = objective(U, V, X, d1, lambda, rows, vals)
	 		#pairwise_error, ndcg = compute_pairwise_error_ndcg(U, V, Y, r, d1, d2, rows_t, vals_t, cols_t, ndcg_k)
	 		rmse = compute_RMSE(U, V, Y, r, d1, d2, rows_t, vals_t, cols_t)
	 		rmse_tr = compute_RMSE_train(U, V, X, r, d1, d2, rows, vals, cols)
			#rmse = compute_RMSE(U, V, Y, r, d1, d2, rows_t, vals_t, cols_t, threshold)
			#rmse_tr = compute_RMSE_train(U, V, X, r, d1, d2, rows, vals, cols, threshold)
			#println("[", iter, ", ", totaltime, ", ", nowobj, ", ", pairwise_error, ", ", ndcg, ", ", rmse_tr, ", ", rmse, "],")
			println("[", iter, ", ", totaltime, ", ", nowobj, ", ", rmse_tr, ", ", rmse, "],")
		end
	end
	return V, U
end


text:
"""
    tpors(xi, eta, ra, dec)

In the tangent plane projection, given the rectangular coordinates
of a star and its spherical coordinates, determine the spherical
coordinates of the tangent point.

### Given ###

- `xi`, `eta`: rectangular coordinates of star image (Note 2)
- `a`, `b`: star's spherical coordinates (Note 3)

### Returned ###

- `status`: number of solutions:
    - 0 = no solutions returned (Note 5)
    - 1 = only the first solution is useful (Note 6)
    - 2 = both solutions are useful (Note 6)
- `a01`, `b01`: tangent point's spherical coordinates, Soln. 1
- `a02`, `b02`: tangent point's spherical coordinates, Soln. 2

### Notes ###

1. The tangent plane projection is also called the "gnomonic
   projection" and the "central projection".

2. The eta axis points due north in the adopted coordinate system.
   If the spherical coordinates are observed (RA,Dec), the tangent
   plane coordinates (xi,eta) are conventionally called the
   "standard coordinates".  If the spherical coordinates are with
   respect to a right-handed triad, (xi,eta) are also right-handed.
   The units of (xi,eta) are, effectively, radians at the tangent
   point.

3. All angular arguments are in radians.

4. The angles a01 and a02 are returned in the range 0-2pi.  The
   angles b01 and b02 are returned in the range +/-pi, but in the
   usual, non-pole-crossing, case, the range is +/-pi/2.

5. Cases where there is no solution can arise only near the poles.
   For example, it is clearly impossible for a star at the pole
   itself to have a non-zero xi value, and hence it is meaningless
   to ask where the tangent point would have to be to bring about
   this combination of xi and dec.

6. Also near the poles, cases can arise where there are two useful
   solutions.  The return value indicates whether the second of the
   two solutions returned is useful;  1 indicates only one useful
   solution, the usual case.

7. The basis of the algorithm is to solve the spherical triangle PSC,
   where P is the north celestial pole, S is the star and C is the
   tangent point.  The spherical coordinates of the tangent point are
   [a0,b0];  writing rho^2 = (xi^2+eta^2) and r^2 = (1+rho^2), side c
   is then (pi/2-b), side p is sqrt(xi^2+eta^2) and side s (to be
   found) is (pi/2-b0).  Angle C is given by sin(C) = xi/rho and
   cos(C) = eta/rho.  Angle P (to be found) is the longitude
   difference between star and tangent point (a-a0).

8. This function is a member of the following set:

   | spherical         |   vector        | solve for |
   |:------------------|:----------------|:----------|
   |  [`tpxes`](@ref)  | [`tpxev`](@ref) | xi,eta    |
   |  [`tpsts`](@ref)  | [`tpstv`](@ref) | star      |
   | *[`tpors`](@ref)* | [`tporv`](@ref) | origin    |

### Called ###

- [`anp`](@ref): normalize angle into range 0 to 2pi

### References ###

- Calabretta M.R. & Greisen, E.W., 2002, "Representations of
  celestial coordinates in FITS", Astron.Astrophys. 395, 1077

- Green, R.M., "Spherical Astronomy", Cambridge University Press,
  1987, Chapter 13.
"""
function tpors(xi, eta, ra, dec)
    az1 = Ref{Cdouble}()
    az2 = Ref{Cdouble}()
    bz1 = Ref{Cdouble}()
    bz2 = Ref{Cdouble}()
    status = ccall((:eraTpors, liberfa), Cint,
                   (Cdouble, Cdouble, Cdouble, Cdouble,
                    Ref{Cdouble}, Ref{Cdouble}, Ref{Cdouble}, Ref{Cdouble}),
                   xi, eta, ra, dec, az1, az2, bz1, bz2)
    return status, az1[], az2[], bz1[], bz2[]
end

"""
    tporv(xi, eta, v)

In the tangent plane projection, given the rectangular coordinates
of a star and its direction cosines, determine the direction
cosines of the tangent point.

### Given ###

- `xi`, `eta`: rectangular coordinates of star image (Note 2)
- `v`: star's direction cosines (Note 3)

### Returned ###

- `status`: number of solutions:
    - 0 = no solutions returned (Note 4)
    - 1 = only the first solution is useful (Note 5)
    - 2 = both solutions are useful (Note 5)
- `v01`: tangent point's direction cosines, Solution 1
- `v02`: tangent point's direction cosines, Solution 2

### Notes ###

1. The tangent plane projection is also called the "gnomonic
   projection" and the "central projection".

2. The eta axis points due north in the adopted coordinate system.
   If the direction cosines represent observed (RA,Dec., the tangent
   plane coordinates (xi,eta. are conventionally called the
   "standard coordinates".  If the direction cosines are with
   respect to a right-handed triad, (xi,eta. are also right-handed.
   The units of (xi,eta. are, effectively, radians at the tangent
   point.

3. The vector v must be of unit length or the result will be wrong.

4. Cases where there is no solution can arise only near the poles.
   For example, it is clearly impossible for a star at the pole
   itself to have a non-zero xi value, and hence it is meaningless
   to ask where the tangent point would have to be.

5. Also near the poles, cases can arise where there are two useful
   solutions.  The return value indicates whether the second of the
   two solutions returned is useful;  1 indicates only one useful
   solution, the usual case.

6. The basis of the algorithm is to solve the spherical triangle
   PSC, where P is the north celestial pole, S is the star and C is
   the tangent point.  Calling the celestial spherical coordinates
   of the star and tangent point (a,b. and (a0,b0) respectively, and
   writing rho^2 = (xi^2+eta^2. and r^2 = (1+rho^2), and
   transforming the vector v into (a,b. in the normal way, side c is
   then (pi/2-b., side p is sqrt(xi^2+eta^2) and side s (to be
   found. is (pi/2-b0), while angle C is given by sin(C) = xi/rho
   and cos(C. = eta/rho;  angle P (to be found) is (a-a0).  After
   solving the spherical triangle, the result (a0,b0. can be
   expressed in vector form as v0.

7. This function is a member of the following set:

   | spherical       | vector            | solve for |
   |:----------------|:------------------|:----------|
   | [`tpxes`](@ref) |  [`tpxev`](@ref)  | xi,eta    |
   | [`tpsts`](@ref) |  [`tpstv`](@ref)  | star      |
   | [`tpors`](@ref) | *[`tporv`](@ref)* | origin    |

### References ###

- Calabretta M.R. & Greisen, E.W., 2002, "Representations of
  celestial coordinates in FITS", Astron.Astrophys. 395, 1077

- Green, R.M., "Spherical Astronomy", Cambridge University Press,
  1987, Chapter 13.
"""
function tporv(xi, eta, v)
    vz1 = Array{Cdouble}(undef, 3)
    vz2 = Array{Cdouble}(undef, 3)
    status = ccall((:eraTporv, liberfa), Cint,
                   (Cdouble, Cdouble, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}),
                   xi, eta, v, vz1, vz2)
    return status, vz1, vz2
end

"""
    tpsts(xi, eta, raz, decz)

In the tangent plane projection, given the star's rectangular
coordinates and the spherical coordinates of the tangent point,
solve for the spherical coordinates of the star.

### Given ###

- `xi`, `eta`: rectangular coordinates of star image (Note 2)
- `a0`, `b0`: tangent point's spherical coordinates

### Returned ###

- `a`, `b`: star's spherical coordinates

### Notes ###

1. The tangent plane projection is also called the "gnomonic
   projection" and the "central projection".

2. The eta axis points due north in the adopted coordinate system.
   If the spherical coordinates are observed (RA,Dec), the tangent
   plane coordinates (xi,eta) are conventionally called the
   "standard coordinates".  If the spherical coordinates are with
   respect to a right-handed triad, (xi,eta) are also right-handed.
   The units of (xi,eta) are, effectively, radians at the tangent
   point.

3. All angular arguments are in radians.

4. This function is a member of the following set:

   |   spherical       | vector           | solve for |
   |:------------------|:-----------------|:----------|
   |  [`tpxes`](@ref)  | [`tpxev`](@ref)  | xi,eta    |
   | *[`tpsts`](@ref)* | [`tpstv`](@ref)  | star      |
   |  [`tpors`](@ref)  | [`tporv`](@ref)  | origin    |

### Called ###

- [`anp`](@ref) normalize angle into range 0 to 2pi

### References ###

- Calabretta M.R. & Greisen, E.W., 2002, "Representations of
  celestial coordinates in FITS", Astron.Astrophys. 395, 1077

- Green, R.M., "Spherical Astronomy", Cambridge University Press,
  1987, Chapter 13.
"""
function tpsts(xi, eta, raz, decz)
    ra = Ref{Cdouble}()
    dec = Ref{Cdouble}()
    ccall((:eraTpsts, liberfa), Cvoid,
          (Cdouble, Cdouble, Cdouble, Cdouble, Ref{Cdouble}, Ref{Cdouble}),
          xi, eta, raz, decz, ra, dec)
    return ra[], dec[]
end

"""
    tpstv(xi, eta, vz)

In the tangent plane projection, given the star's rectangular
coordinates and the direction cosines of the tangent point, solve
for the direction cosines of the star.

### Given ###

- `xi`, `eta`: rectangular coordinates of star image (Note 2)
- `v0`: tangent point's direction cosines

### Returned ###

- `v`: star's direction cosines

### Notes ###

1. The tangent plane projection is also called the "gnomonic
   projection" and the "central projection".

2. The eta axis points due north in the adopted coordinate system.
   If the direction cosines represent observed (RA,Dec), the tangent
   plane coordinates (xi,eta) are conventionally called the
   "standard coordinates".  If the direction cosines are with
   respect to a right-handed triad, (xi,eta) are also right-handed.
   The units of (xi,eta) are, effectively, radians at the tangent
   point.

3. The method used is to complete the star vector in the (xi,eta)
   based triad and normalize it, then rotate the triad to put the
   tangent point at the pole with the x-axis aligned to zero
   longitude.  Writing (a0,b0) for the celestial spherical
   coordinates of the tangent point, the sequence of rotations is
   (b-pi/2) around the x-axis followed by (-a-pi/2) around the
   z-axis.

4. If vector v0 is not of unit length, the returned vector v will
   be wrong.

5. If vector v0 points at a pole, the returned vector v will be
   based on the arbitrary assumption that the longitude coordinate
   of the tangent point is zero.

6. This function is a member of the following set:

   | spherical       | vector            | solve for |
   |:----------------|:------------------|:----------|
   | [`tpxes`](@ref) |  [`tpxev`](@ref)  | xi,eta    |
   | [`tpsts`](@ref) | *[`tpstv`](@ref)* | star      |
   | [`tpors`](@ref) |  [`tporv`](@ref)  | origin    |

### References ###

- Calabretta M.R. & Greisen, E.W., 2002, "Representations of
  celestial coordinates in FITS", Astron.Astrophys. 395, 1077

- Green, R.M., "Spherical Astronomy", Cambridge University Press,
  1987, Chapter 13.
"""
function tpstv(xi, eta, vz)
    v = Array{Cdouble}(undef, 3)
    ccall((:eraTpstv, liberfa), Cvoid,
          (Cdouble, Cdouble, Ptr{Cdouble}, Ptr{Cdouble}),
          xi, eta, vz, v)
    return v
end

"""
    tpxes(ra, dec, raz, decz)

In the tangent plane projection, given celestial spherical
coordinates for a star and the tangent point, solve for the star's
rectangular coordinates in the tangent plane.

### Given ###

- `a`, `b`: star's spherical coordinates
- `a0`, `b0`: tangent point's spherical coordinates

### Returned ###

- `status`:
    - 0 = OK
    - 1 = star too far from axis
    - 2 = antistar on tangent plane
    - 3 = antistar too far from axis
- `xi`, `eta`: rectangular coordinates of star image (Note 2)

### Notes ###

1. The tangent plane projection is also called the "gnomonic
   projection" and the "central projection".

2. The eta axis points due north in the adopted coordinate system.
   If the spherical coordinates are observed (RA,Dec), the tangent
   plane coordinates (xi,eta) are conventionally called the
   "standard coordinates".  For right-handed spherical coordinates,
   (xi,eta) are also right-handed.  The units of (xi,eta) are,
   effectively, radians at the tangent point.

3. All angular arguments are in radians.

4. This function is a member of the following set:

   | spherical         | vector          | solve for |
   |:------------------|:----------------|:----------|
   | *[`tpxes`](@ref)* | [`tpxev`](@ref) | xi,eta    |
   |  [`tpsts`](@ref)  | [`tpstv`](@ref) | star      |
   |  [`tpors`](@ref)  | [`tporv`](@ref) | origin    |

### References ###

- Calabretta M.R. & Greisen, E.W., 2002, "Representations of
  celestial coordinates in FITS", Astron.Astrophys. 395, 1077

- Green, R.M., "Spherical Astronomy", Cambridge University Press,
  1987, Chapter 13.
"""
function tpxes(ra, dec, raz, decz)
    xi = Ref{Cdouble}()
    eta = Ref{Cdouble}()
    status = ccall((:eraTpxes, liberfa), Cint,
                   (Cdouble, Cdouble, Cdouble, Cdouble, Ref{Cdouble}, Ref{Cdouble}),
                   ra, dec, raz, decz, xi, eta)
    return status, xi[], eta[]
end

"""
    tpxev(v, vz)

In the tangent plane projection, given celestial direction cosines
for a star and the tangent point, solve for the star's rectangular
coordinates in the tangent plane.

### Given ###

- `v`: direction cosines of star (Note 4)
- `v0`: direction cosines of tangent point (Note 4)

### Returned ###

- `status`:
    - 0 = OK
    - 1 = star too far from axis
    - 2 = antistar on tangent plane
    - 3 = antistar too far from axis
- `xi`, `eta`: tangent plane coordinates of star

### Notes ###

1. The tangent plane projection is also called the "gnomonic
   projection" and the "central projection".

2. The eta axis points due north in the adopted coordinate system.
   If the direction cosines represent observed (RA,Dec), the tangent
   plane coordinates (xi,eta) are conventionally called the
   "standard coordinates".  If the direction cosines are with
   respect to a right-handed triad, (xi,eta) are also right-handed.
   The units of (xi,eta) are, effectively, radians at the tangent
   point.

3. The method used is to extend the star vector to the tangent
   plane and then rotate the triad so that (x,y) becomes (xi,eta).
   Writing (a,b) for the celestial spherical coordinates of the
   star, the sequence of rotations is (a+pi/2) around the z-axis
   followed by (pi/2-b) around the x-axis.

4. If vector v0 is not of unit length, or if vector v is of zero
   length, the results will be wrong.

5. If v0 points at a pole, the returned (xi,eta) will be based on
   the arbitrary assumption that the longitude coordinate of the
   tangent point is zero.

6. This function is a member of the following set:

   | spherical       | vector            | solve for |
   |:----------------|:------------------|:----------|
   | [`tpxes`](@ref) | *[`tpxev`](@ref)* | xi,eta    |
   | [`tpsts`](@ref) |  [`tpstv`](@ref)  | star      |
   | [`tpors`](@ref) |  [`tporv`](@ref)  | origin    |

### References ###

- Calabretta M.R. & Greisen, E.W., 2002, "Representations of
  celestial coordinates in FITS", Astron.Astrophys. 395, 1077

- Green, R.M., "Spherical Astronomy", Cambridge University Press,
  1987, Chapter 13.
"""
function tpxev(v, vz)
    xi = Ref{Cdouble}()
    eta = Ref{Cdouble}()
    status = ccall((:eraTpxev, liberfa), Cint,
                   (Ptr{Cdouble}, Ptr{Cdouble}, Ref{Cdouble}, Ref{Cdouble}),
                   v, vz, xi, eta)
    return status, xi[], eta[]
end

"""
    tr(r)

Transpose an r-matrix.

!!! warning "Deprecated"
    Use `r'` instead.

### Given ###

- `r`: R-matrix

### Returned ###

- `rt`: Transpose

### Note ###

   It is permissible for r and rt to be the same array.

### Called ###

- [`cr`](@ref): copy r-matrix

"""
tr

function _tr(r)
    @checkdims 3 3 r
    rt = zeros(Cdouble, 3, 3)
    ccall((:eraTr, liberfa), Cvoid,
          (Ptr{Cdouble}, Ptr{Cdouble}),
          permutedims(r), rt)
    return permutedims(rt)
end

@deprecate tr(r) r'

"""
    trxpv(r, pv)

Multiply a pv-vector by the transpose of an r-matrix.

!!! warning "Deprecated"
    Use `[r' * pv[1], r' * pv[2]]` instead.

### Given ###

- `r`: R-matrix
- `pv`: Pv-vector

### Returned ###

- `trpv`: R * pv

### Note ###

   It is permissible for pv and trpv to be the same array.

### Called ###

- [`tr`](@ref): transpose r-matrix
- [`rxpv`](@ref): product of r-matrix and pv-vector

"""
trxpv

function _trxpv(r, pv)
    @checkdims 3 3 r
    _pv = array_to_cmatrix(pv; n=3)
    rp = zeros(Cdouble, 3, 2)
    ccall((:eraTrxpv, liberfa), Cvoid,
            (Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}),
            permutedims(r), _pv, rp)
    return cmatrix_to_array(rp)
end

@deprecate trxpv(r, pv) [r' * pv[1], r' * pv[2]]

"""
    trxp(r, p)

Multiply a p-vector by the transpose of an r-matrix.

!!! warning "Deprecated"
    Use `r' * p` instead.

### Given ###

- `r`: R-matrix
- `p`: P-vector

### Returned ###

- `trp`: R * p

### Note ###

   It is permissible for p and trp to be the same array.

### Called ###

- [`tr`](@ref): transpose r-matrix
- [`rxp`](@ref): product of r-matrix and p-vector

"""
trxp

function _trxp(r, p)
    @checkdims 3 3 r
    @checkdims 3 p
    rp = zeros(Cdouble, 3)
    ccall((:eraTrxp, liberfa), Cvoid,
            (Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}),
            permutedims(r), p, rp)
    return rp
end

@deprecate trxp(r, p) r' * p

"""
    taiut1(tai1, tai2, dta)

Time scale transformation:  International Atomic Time, TAI, to
Universal Time, UT1.

### Given ###

- `tai1`, `tai2`: TAI as a 2-part Julian Date
- `dta`: UT1-TAI in seconds

### Returned ###

- `ut11`, `ut12`: UT1 as a 2-part Julian Date

### Notes ###

1. tai1+tai2 is Julian Date, apportioned in any convenient way
   between the two arguments, for example where tai1 is the Julian
   Day Number and tai2 is the fraction of a day.  The returned
   UT11,UT12 follow suit.

2. The argument dta, i.e. UT1-TAI, is an observed quantity, and is
   available from IERS tabulations.

### Reference ###

- Explanatory Supplement to the Astronomical Almanac,
    P. Kenneth Seidelmann (ed), University Science Books (1992)

"""
taiut1

"""
    tdbtt(tdb1, tdb2, dtr)

Time scale transformation:  Barycentric Dynamical Time, TDB, to
Terrestrial Time, TT.

### Given ###

- `tdb1`, `tdb2`: TDB as a 2-part Julian Date
- `dtr`: TDB-TT in seconds

### Returned ###

- `tt1`, `tt2`: TT as a 2-part Julian Date

### Notes ###

1. tdb1+tdb2 is Julian Date, apportioned in any convenient way
   between the two arguments, for example where tdb1 is the Julian
   Day Number and tdb2 is the fraction of a day.  The returned
   tt1,tt2 follow suit.

2. The argument dtr represents the quasi-periodic component of the
   GR transformation between TT and TCB.  It is dependent upon the
   adopted solar-system ephemeris, and can be obtained by numerical
   integration, by interrogating a precomputed time ephemeris or by
   evaluating a model such as that implemented in the ERFA function
   [`dtdb`](@ref).   The quantity is dominated by an annual term of 1.7 ms
   amplitude.

3. TDB is essentially the same as Teph, the time argument for the
   JPL solar system ephemerides.

### References ###

- McCarthy, D. D., Petit, G. (eds.), IERS Conventions (2003),
    IERS Technical Note No. 32, BKG (2004)

- IAU 2006 Resolution 3

"""
tdbtt

"""
    tttdb(tt1, tt2, dtr)

Time scale transformation:  Terrestrial Time, TT, to Barycentric
Dynamical Time, TDB.

### Given ###

- `tt1`, `tt2`: TT as a 2-part Julian Date
- `dtr`: TDB-TT in seconds

### Returned ###

- `tdb1`, `tdb2`: TDB as a 2-part Julian Date

### Notes ###

1. tt1+tt2 is Julian Date, apportioned in any convenient way between
   the two arguments, for example where tt1 is the Julian Day Number
   and tt2 is the fraction of a day.  The returned tdb1,tdb2 follow
   suit.

2. The argument dtr represents the quasi-periodic component of the
   GR transformation between TT and TCB.  It is dependent upon the
   adopted solar-system ephemeris, and can be obtained by numerical
   integration, by interrogating a precomputed time ephemeris or by
   evaluating a model such as that implemented in the ERFA function
   [`dtdb`](@ref).   The quantity is dominated by an annual term of 1.7 ms
   amplitude.

3. TDB is essentially the same as Teph, the time argument for the JPL
   solar system ephemerides.

### References ###

- McCarthy, D. D., Petit, G. (eds.), IERS Conventions (2003),
    IERS Technical Note No. 32, BKG (2004)

- IAU 2006 Resolution 3

"""
tttdb

"""
    ttut1(tt1, tt2, dt)

Time scale transformation:  Terrestrial Time, TT, to Universal Time,
UT1.

### Given ###

- `tt1`, `tt2`: TT as a 2-part Julian Date
- `dt`: TT-UT1 in seconds

### Returned ###

- `ut11`, `ut12`: UT1 as a 2-part Julian Date

### Notes ###

1. tt1+tt2 is Julian Date, apportioned in any convenient way between
   the two arguments, for example where tt1 is the Julian Day Number
   and tt2 is the fraction of a day.  The returned ut11,ut12 follow
   suit.

2. The argument dt is classical Delta T.

### Reference ###

- Explanatory Supplement to the Astronomical Almanac,
    P. Kenneth Seidelmann (ed), University Science Books (1992)

"""
ttut1

for name in ("taiut1",
             "tdbtt",
             "tttdb",
             "ttut1")
    f = Symbol(name)
    fc = "era" * uppercasefirst(name)
    @eval begin
        function ($f)(a, b, c)
            r1 = Ref{Cdouble}()
            r2 = Ref{Cdouble}()
            i = ccall(($fc, liberfa), Cint,
                      (Cdouble, Cdouble, Cdouble, Ref{Cdouble}, Ref{Cdouble}),
                      a, b, c, r1, r2)
            @assert i == 0
            r1[], r2[]
        end
    end
end

"""
    taitt(tai1, tai2)

Time scale transformation:  International Atomic Time, TAI, to
Terrestrial Time, TT.

### Given ###

- `tai1`, `tai2`: TAI as a 2-part Julian Date

### Returned ###

- `tt1`, `tt2`: TT as a 2-part Julian Date

### Note ###

   tai1+tai2 is Julian Date, apportioned in any convenient way
   between the two arguments, for example where tai1 is the Julian
   Day Number and tai2 is the fraction of a day.  The returned
   tt1,tt2 follow suit.

### References ###

- McCarthy, D. D., Petit, G. (eds.), IERS Conventions (2003),
    IERS Technical Note No. 32, BKG (2004)

- Explanatory Supplement to the Astronomical Almanac,
    P. Kenneth Seidelmann (ed), University Science Books (1992)

"""
taitt

"""
    taiutc(tai1, tai2)

Time scale transformation:  International Atomic Time, TAI, to
Coordinated Universal Time, UTC.

### Given ###

- `tai1`, `tai2`: TAI as a 2-part Julian Date (Note 1)

### Returned ###

- `utc1`, `utc2`: UTC as a 2-part quasi Julian Date (Notes 1-3)

### Notes ###

1. tai1+tai2 is Julian Date, apportioned in any convenient way
   between the two arguments, for example where tai1 is the Julian
   Day Number and tai2 is the fraction of a day.  The returned utc1
   and utc2 form an analogous pair, except that a special convention
   is used, to deal with the problem of leap seconds - see the next
   note.

2. JD cannot unambiguously represent UTC during a leap second unless
   special measures are taken.  The convention in the present
   function is that the JD day represents UTC days whether the
   length is 86399, 86400 or 86401 SI seconds.  In the 1960-1972 era
   there were smaller jumps (in either direction) each time the
   linear UTC(TAI) expression was changed, and these "mini-leaps"
   are also included in the ERFA convention.

3. The function [`d2dtf`](@ref) can be used to transform the UTC quasi-JD
   into calendar date and clock time, including UTC leap second
   handling.

4. The warning status "dubious year" flags UTCs that predate the
   introduction of the time scale or that are too far in the future
   to be trusted.  See [`dat`](@ref) for further details.

### Called ###

- [`utctai`](@ref): UTC to TAI

### References ###

- McCarthy, D. D., Petit, G. (eds.), IERS Conventions (2003),
    IERS Technical Note No. 32, BKG (2004)

- Explanatory Supplement to the Astronomical Almanac,
    P. Kenneth Seidelmann (ed), University Science Books (1992)

"""
taiutc

"""
    tcbtdb(tcb1, tcb2)

Time scale transformation:  Barycentric Coordinate Time, TCB, to
Barycentric Dynamical Time, TDB.

### Given ###

- `tcb1`, `tcb2`: TCB as a 2-part Julian Date

### Returned ###

- `tdb1`, `tdb2`: TDB as a 2-part Julian Date

### Notes ###

1. tcb1+tcb2 is Julian Date, apportioned in any convenient way
   between the two arguments, for example where tcb1 is the Julian
   Day Number and tcb2 is the fraction of a day.  The returned
   tdb1,tdb2 follow suit.

2. The 2006 IAU General Assembly introduced a conventional linear
   transformation between TDB and TCB.  This transformation
   compensates for the drift between TCB and terrestrial time TT,
   and keeps TDB approximately centered on TT.  Because the
   relationship between TT and TCB depends on the adopted solar
   system ephemeris, the degree of alignment between TDB and TT over
   long intervals will vary according to which ephemeris is used.
   Former definitions of TDB attempted to avoid this problem by
   stipulating that TDB and TT should differ only by periodic
   effects.  This is a good description of the nature of the
   relationship but eluded precise mathematical formulation.  The
   conventional linear relationship adopted in 2006 sidestepped
   these difficulties whilst delivering a TDB that in practice was
   consistent with values before that date.

3. TDB is essentially the same as Teph, the time argument for the
   JPL solar system ephemerides.

### Reference ###

- IAU 2006 Resolution B3

"""
tcbtdb

"""
    tcgtt(tcg1, tcg2)

Time scale transformation:  Geocentric Coordinate Time, TCG, to
Terrestrial Time, TT.

### Given ###

- `tcg1`, `tcg2`: TCG as a 2-part Julian Date

### Returned ###

- `tt1`, `tt2`: TT as a 2-part Julian Date

### Note ###

   tcg1+tcg2 is Julian Date, apportioned in any convenient way
   between the two arguments, for example where tcg1 is the Julian
   Day Number and tcg22 is the fraction of a day.  The returned
   tt1,tt2 follow suit.

### References ###

- McCarthy, D. D., Petit, G. (eds.), IERS Conventions (2003),.
    IERS Technical Note No. 32, BKG (2004)

- IAU 2000 Resolution B1.9

"""
tcgtt

"""
    tdbtcb(tdb1, tdb2)

Time scale transformation:  Barycentric Dynamical Time, TDB, to
Barycentric Coordinate Time, TCB.

### Given ###

- `tdb1`, `tdb2`: TDB as a 2-part Julian Date

### Returned ###

- `tcb1`, `tcb2`: TCB as a 2-part Julian Date

### Notes ###

1. tdb1+tdb2 is Julian Date, apportioned in any convenient way
   between the two arguments, for example where tdb1 is the Julian
   Day Number and tdb2 is the fraction of a day.  The returned
   tcb1,tcb2 follow suit.

2. The 2006 IAU General Assembly introduced a conventional linear
   transformation between TDB and TCB.  This transformation
   compensates for the drift between TCB and terrestrial time TT,
   and keeps TDB approximately centered on TT.  Because the
   relationship between TT and TCB depends on the adopted solar
   system ephemeris, the degree of alignment between TDB and TT over
   long intervals will vary according to which ephemeris is used.
   Former definitions of TDB attempted to avoid this problem by
   stipulating that TDB and TT should differ only by periodic
   effects.  This is a good description of the nature of the
   relationship but eluded precise mathematical formulation.  The
   conventional linear relationship adopted in 2006 sidestepped
   these difficulties whilst delivering a TDB that in practice was
   consistent with values before that date.

3. TDB is essentially the same as Teph, the time argument for the
   JPL solar system ephemerides.

### Reference ###

- IAU 2006 Resolution B3

"""
tdbtcb

"""
    tttai(tt1, tt2)

Time scale transformation:  Terrestrial Time, TT, to International
Atomic Time, TAI.

### Given ###

- `tt1`, `tt2`: TT as a 2-part Julian Date

### Returned ###

- `tai1`, `tai2`: TAI as a 2-part Julian Date

### Note ###

   tt1+tt2 is Julian Date, apportioned in any convenient way between
   the two arguments, for example where tt1 is the Julian Day Number
   and tt2 is the fraction of a day.  The returned tai1,tai2 follow
   suit.

### References ###

- McCarthy, D. D., Petit, G. (eds.), IERS Conventions (2003),
    IERS Technical Note No. 32, BKG (2004)

- Explanatory Supplement to the Astronomical Almanac,
    P. Kenneth Seidelmann (ed), University Science Books (1992)

"""
tttai

"""
    tttcg(tt1, tt2)

Time scale transformation:  Terrestrial Time, TT, to Geocentric
Coordinate Time, TCG.

### Given ###

- `tt1`, `tt2`: TT as a 2-part Julian Date

### Returned ###

- `tcg1`, `tcg2`: TCG as a 2-part Julian Date

### Note ###

   tt1+tt2 is Julian Date, apportioned in any convenient way between
   the two arguments, for example where tt1 is the Julian Day Number
   and tt2 is the fraction of a day.  The returned tcg1,tcg2 follow
   suit.

### References ###

- McCarthy, D. D., Petit, G. (eds.), IERS Conventions (2003),
    IERS Technical Note No. 32, BKG (2004)

- IAU 2000 Resolution B1.9

"""
tttcg

for name in ("taitt",
             "taiutc",
             "tcbtdb",
             "tcgtt",
             "tdbtcb",
             "tttai",
             "tttcg")
    f = Symbol(name)
    fc = "era" * uppercasefirst(name)
    @eval begin
        function ($f)(a, b)
            r1 = Ref{Cdouble}()
            r2 = Ref{Cdouble}()
            i = ccall(($fc, liberfa), Cint,
                      (Cdouble, Cdouble, Ref{Cdouble}, Ref{Cdouble}),
                      a, b, r1, r2)
            @assert i == 0
            return r1[], r2[]
        end
    end
end

"""
    tf2a(s, ihour, imin, sec)

Convert hours, minutes, seconds to radians.

### Given ###

- `s`: Sign:  '-' = negative, otherwise positive
- `ihour`: Hours
- `imin`: Minutes
- `sec`: Seconds

### Returned ###

- `rad`: Angle in radians

### Notes ###

1.  The result is computed even if any of the range checks fail.

2.  Negative ihour, imin and/or sec produce a warning status, but
    the absolute value is used in the conversion.

3.  If there are multiple errors, the status value reflects only the
    first, the smallest taking precedence.

"""
tf2a

"""
    tf2d(s, ihour, imin, sec)

Convert hours, minutes, seconds to days.

### Given ###

- `s`: Sign:  '-' = negative, otherwise positive
- `ihour`: Hours
- `imin`: Minutes
- `sec`: Seconds

### Returned ###

- `days`: Interval in days

### Notes ###

1.  The result is computed even if any of the range checks fail.

2.  Negative ihour, imin and/or sec produce a warning status, but
    the absolute value is used in the conversion.

3.  If there are multiple errors, the status value reflects only the
    first, the smallest taking precedence.

"""
tf2d

for name in ("tf2a",
             "tf2d")
    f = Symbol(name)
    fc = "era" * uppercasefirst(name)
    @eval begin
        function ($f)(s, ideg, iamin, asec)
            rad = Ref{Cdouble}()
            i = ccall(($fc, liberfa), Cint,
                      (Cchar, Cint, Cint, Cdouble, Ref{Cdouble}),
                       s, ideg, iamin, asec, rad)
            if i == 1
                throw(ERFAException("ihour outside range 0-23"))
            elseif i == 2
                throw(ERFAException("imin outside range 0-59"))
            elseif i == 3
                throw(ERFAException("sec outside range 0-59.999..."))
            end
            return rad[]
        end
    end
end


text:
export iCycle, iPath, iGrid, iKnight, iCube, iShift

"""
`iCycle(n::Int, simple::Bool=true)` creates an implicit graph that is
a cycle graph with `n` vertices `1` through `n`. When `simple`
is `true`, the graph is undirected; when `false` it's a directed
cycle `1 → 2 → 3 → ⋯ → n → 1`.
"""
function iCycle(n::Int, simple::Bool = true)::ImplicitGraph{Int}
    if n < 3
        error("Number of vertices must be at least 3")
    end

    function N(v::Int)
        a = mod1(v + 1, n)
        b = mod1(v - 1, n)
        if simple
            return [b, a]
        end
        return [a]
    end

    has_vertex(v::Int) = 1 <= v <= n

    return ImplicitGraph{Int}(has_vertex, N)
end

"""
`iPath(simple::Bool=true)` creates an implicit graph that is
a path graph on the integers. If `simple` is `true` vertex `v` is
has an edge to `v-1` and `v+1`; otherwise, `v` has an edge only to
`v+1`.
"""
function iPath(simple::Bool = true)::ImplicitGraph{Int}
    yes(v::Int)::Bool = true

    function N(v::Int)::Vector{Int}
        if simple
            return [v - 1, v + 1]
        else
            return [v + 1]
        end
    end
    return ImplicitGraph{Int}(yes, N)
end

"""
`iGrid()` returns an infinite two-dimensional grid graph. Vertices are
of type `Tuple{Int,Int}`.
"""
function iGrid()::ImplicitGraph{Tuple{Int,Int}}
    yes(v::Tuple{Int,Int})::Bool = true
    function N(v::Tuple{Int,Int})::Vector{Tuple{Int,Int}}
        a, b = v
        return [(a, b - 1), (a, b + 1), (a - 1, b), (a + 1, b)]
    end
    return ImplicitGraph{Tuple{Int,Int}}(yes, N)
end

"""
`iKnight()` returns the Knight's move graph on an infinite
chessboard.
"""
function iKnight()::ImplicitGraph{Tuple{Int,Int}}
    yes(v::Tuple{Int,Int})::Bool = true
    function N(v::Tuple{Int,Int})::Vector{Tuple{Int,Int}}
        a, b = v
        neigh = [
            (a + 1, b + 2),
            (a + 1, b - 2),
            (a + 2, b + 1),
            (a + 2, b - 1),
            (a - 1, b + 2),
            (a - 1, b - 2),
            (a - 2, b + 1),
            (a - 2, b - 1),
        ]
        return neigh
    end
    return ImplicitGraph{Tuple{Int,Int}}(yes, N)
end

"""
`iCube(d::Int)` creates an (implict) `d`-dimensional cube graph.
"""
function iCube(d::Int)::ImplicitGraph{String}
    if d < 1
        error("Dimension must be positive")
    end

    function dvec_check(s::String)::Bool
        if length(s) != d
            return false
        end
        for i = 1:d
            if s[i] ∉ "01"
                return false
            end
        end
        return true
    end

    function N(v::String)
        result = Vector{String}(undef, d)
        for i = 1:d
            head = v[1:i-1]
            c = v[i] == '0' ? "1" : "0"
            tail = v[i+1:end]
            result[i] = head * c * tail
        end
        return result
    end

    return ImplicitGraph{String}(dvec_check, N)
end

using IterTools
"""
`iShift(alphabet,n)` creates an implicit shift digraph whose vertices
are `n`-tuples of elements of `alphabet`.
"""
function iShift(alphabet, n::Int)::ImplicitGraph
    elts = collect(distinct(alphabet))
    T = eltype(elts)

    function has_vertex(v)::Bool
        for i = 1:n
            if v[i] ∉ elts
                return false
            end
        end
        return true
    end

    function N(v)
        base = v[2:end]
        result = [(base..., j) for j in elts]
        return result
    end

    return ImplicitGraph{NTuple{n,T}}(has_vertex, N)
end


text:
"""
## Visualization tool for the clawpack output
VisClaw.jl is a Julia package for plotting simulation results of the clawpack.\n
https://github.com/hydrocoast/VisClaw.jl

### Examples
- $(dirname(dirname(pathof(VisClaw))))/Examples_using_Plots.ipynb
- $(dirname(dirname(pathof(VisClaw))))/Examples_using_GMT.ipynb

### Author
Takuya Miyashita (miyashita@hydrocoast.jp)\n
Doctoral student, Kyoto University, 2018\n
"""
module VisClaw

using Statistics: mean
using DelimitedFiles: readdlm
using Interpolations: Interpolations
using Printf
using Dates
using NetCDF: NetCDF

using GeometricalPredicates: GeometricalPredicates
using Plots: Plots

## define CLAW path from shell
include("clawpath.jl")
export CLAW

## define structs and basic functions
const KWARG = Dict{Symbol,Any}
const emptyF = Array{Float64}(undef, 0, 0)
const timedict = Dict(:second => 1.0, :minute => 60.0, :hour => 3600.0, :day => 24*3600.0)
const varnameset(D,k,v) = haskey(D,k) ? k : v

include("structclaw.jl")
include("amrutils.jl")
include("replaceunit.jl")
include("converttodatetime.jl")
include("getvarname_nctopo.jl")

## load
include("loaddata.jl")
include("loadtrack.jl")
include("loadtopo.jl")
include("loadfgmaxdata.jl")
include("loadfgmax.jl")
include("loadfort.jl")
include("loadgauge.jl")

## print
include("printtopo.jl")

## convert data
include("gaugemax.jl")
include("gaugeinterp.jl")
include("coarsegridmask.jl")

## setup
include("plotsargs.jl")
include("plotstools.jl")

## plot (using Plots)
include("plots2d.jl")
include("plotscheck.jl")
include("plotstopo.jl")
include("plotsgaugewaveform.jl")
include("plotsgaugevelocity.jl")
include("plotsgaugelocation.jl")
include("plotsfgmax.jl")
include("plotstrack.jl")

## plot (using GMT)
using GMT:GMT
include("gmttools.jl")
include("gmttopo.jl")
include("gmtgauge.jl")
include("gmtsurface.jl")
include("gmtarrows.jl")
include("gmttrack.jl")

## general functions
export geodata, amrdata, surgedata, gaugedata, fgmaxdata, regiondata
export topodata, dtopodata
export loadfgmax
export loadtopo, loaddeform, loaddtopo
export loadgauge
export loadtrack
export loadsurface, loadcurrent, loadstorm
export printtopoESRI, printtopo, printdtopo
export coarsegridmask!
export axesratio
export replaceunit!, converttodatetime!
export gaugemax, gaugeinterp

## functions with Plots.jl
export plotsamr
export plotscheck
export gridnumber!, tilebound!
export plotscoastline, plotscoastline!
export plotsfgmax, plotsfgmax!
export plotstopo, plotstopo!
export plotsdtopo, plotsdtopo!
export plotstoporange, plotstoporange!
export plotsgaugelocation, plotsgaugelocation!
export plotsgaugewaveform, plotsgaugewaveform!
export plotsgaugevelocity, plotsgaugevelocity!
export plotstrack, plotstrack!
export plotsgif, plotssavefig

## functions with GMT.jl
export getR, getR_tile, getJ, geogrd
export landmask_asc, landmask_grd
export tilegrd_xyz, tilegrd, tilegrd_mask
export arrowgrd, arrowscalegrd
export gmttopo
export gmtgaugewaveform, gmtgaugewaveform!
export gmtgaugevelocity, gmtgaugevelocity!
export gmtgaugelocation, gmtgaugelocation!
export gmtgaugeannotation!
export gmttoporange!
export gmtcoastline, gmtcoastline!
export gmttrack, gmttrack!

## uniform-grid interpolation
using PyCall: PyCall
include("scipyinterp.jl")

export interpsurface

end


text:
# Conventional Beamformer (CBF)
# Construct filter before beamforming

mutable struct CBFFilterConfig
  fFloor::Float64
  fCeil::Float64
  dataLen::Int
  nPhones::Int
  azimuths::Array{Float64}
  soundSpeed::Float64
  FFTfreqs::Array{Float64}
end

function compare(a::CBFFilterConfig,b::CBFFilterConfig)
  TP = true
  TP = TP && a.fFloor == b.fFloor       # Float64
  TP = TP && a.fCeil == b.fCeil       # Float64
  TP = TP && a.dataLen == b.dataLen       # Int
  TP = TP && a.nPhones == b.nPhones       # Int
  TP = TP && a.azimuths == b.azimuths       # Array{Float64}
  TP = TP && a.soundSpeed == b.soundSpeed       # Float64
  TP = TP && a.FFTfreqs == b.FFTfreqs
  TP
end

function getCBFFilter2Dsize(thisCfg::CBFFilterConfig)
    return (length(thisCfg.azimuths),thisCfg.dataLen,thisCfg.nPhones)
end

function constructCBFFilter2D!(thisCfg::CBFFilterConfig,arrayPos::Array,filterOut::Array,sourceXY::Array, dataHolder::Array)

    # combination of fCeil and fSampling limit information content extracted by thirsCfg.azimuths. Sample faster to get higher res BF
    for iter in eachindex(thisCfg.azimuths)
        @inbounds sourceXY[1] = -cos(thisCfg.azimuths[iter])
        @inbounds sourceXY[2] = -sin(thisCfg.azimuths[iter])
        #tDelays[:] = (arrayPos * sourceXY)/thisCfg.soundSpeed

        for piter in 1:thisCfg.nPhones
            @inbounds dt = (arrayPos[piter,1]*sourceXY[1] + arrayPos[piter,2]*sourceXY[2])/thisCfg.soundSpeed

            copyto!(view(dataHolder,piter,:),thisCfg.FFTfreqs)
            lmul!(-2im*pi*dt,view(dataHolder,piter,:))

            # for fftiter in 1:thisCfg.dataLen
            #     @inbounds dataHolder[piter,fftiter] = exp(-2im*pi*dt*thisCfg.FFTfreqs[fftiter])
            # end
        end
        # @time filterOut[iter,:,:] = transpose(conj!(exp.(-2im*pi*tDelays*FFTfreqs'))) #conjugate transpose in place
        @inbounds adjoint!(view(filterOut,iter,:,:),exp.(dataHolder))
    end
    nothing
end

function CBF2D_DelaySum!(thisCfg::CBFFilterConfig,dataIn::Array,dataOut::Array,temp1::Array,temp2::Array,thisFilter::Array; dochecks::Bool=false)
    if dochecks
      if size(dataIn,1) !=thisCfg.dataLen || size(dataIn,2)!=thisCfg.nPhones
          error("Data Array wrong dimensions")
      end
      if size(thisFilter)!=(length(thisCfg.azimuths),thisCfg.dataLen,thisCfg.nPhones)
          error("Provided Filter does not match provided config.")
      end
    end
    # minimal syntax
#    for iter in eachindex(thisCfg.azimuths)
#        dataOut[iter] = sum(norm.(sum(dataIn.*thisFilter[iter,:,:],2),2),1)[1]
#    end

    for iter in eachindex(thisCfg.azimuths)
        # for diter in 1:thisCfg.dataLen
        #      for phiter in 1:thisCfg.nPhones
        #          temp2[diter,phiter] = dataIn[diter,phiter]
        #      end
        #  end
        copyto!(temp2,dataIn)
        #temp2 = deepcopy(dataIn) # think somethign not zeroing properly?
        @fastmath temp2 .*= @view thisFilter[iter,:,:]
        @fastmath @inbounds sum!(temp1,temp2)

        @fastmath @inbounds @simd for ditr in 1:thisCfg.dataLen
            temp1[ditr] = norm(temp1[ditr])
        end

        @inbounds dataOut[iter] = sum(temp1)
        #dataOut[iter] = sum(norm.(temp1,2))
    end

    # dataOut = deepcopy(dataOut)

    # brute force in-place syntax
    # fill!(dataOut, 0.0)
    # phones_combined_ = zeros(thisCfg.dataLen) .+ im*0.0
    # for iter in eachindex(thisCfg.azimuths)
    #   fill!(phones_combined_, 0.0)
    #   for diter in 1:thisCfg.dataLen
    #     for phiter in 1:thisCfg.nPhones
    #       phones_combined_[diter] += dataIn[diter,phiter]*thisFilter[iter,diter,phiter]
    #     end
    #     dataOut[iter] += norm(phones_combined_[diter],2)
    #   end
    # end
    # @assert norm(dataOut - dataOutTest) < 0.00001

    # TODO - some hybrid clean in-place syntax
    nothing
end

"""
    $(SIGNATURES)

Phase shifting ztransWave contents of leave one out element of `elemPositions` (positions of array elements)
based on dx,dy position perturbation.  The ztransWave is chirp z-transform of LOO elements's waveform data,
and only considering subset of frequencies in ztransWave (i.e. smaller than size of raw waveform data)
"""
function phaseShiftSingle!(sourceXY::Vector{Float64},
                           thisCfg::CBFFilterConfig,
                           azi::R,
                           arrayElemPos::Array,
                           ztransWave::Array{Complex{R2}}  ) where {R <: Real, R2 <: Real}

  #sourceXY[1] = -cos(azi) # unit vector in azimuth direction
  #sourceXY[2] = -sin(azi)
  # eliminating data placeholder/ no impact on mallocs
  dt = (arrayElemPos[1] * -cos(azi) + arrayElemPos[2] * -sin(azi))/thisCfg.soundSpeed
  #dt = (positions' * sourceXY)/thisCfg.soundSpeed
  ztransWave .*= conj(exp.(-2im*pi*dt*thisCfg.FFTfreqs))

  # @assert thisCfg.dataLen == length(ztransWave)
  # for j in 1:thisCfg.dataLen
  #    ztransWave[j] *= conj(exp(-2im*pi*dt*thisCfg.FFTfreqs[j]))
  # end
  nothing
end

"""
    $(SIGNATURES)

Generate the incoming waveform (into `bfOutLIEl`) at the desired look angle `aziInd`.
"""
function liebf!(bfOutLIEl, cztWaveLIEl, bfFilterLIEl, aziIndl, temp; normalize::Bool=false)

  # perform conventional beam forming with single look angle (on leave-ins)
  #bfOutLIEl[:] = sum(cztWaveLIEl.*bfFilterLIEl[aziIndl,:,:],2)

  copy!(temp,cztWaveLIEl)
  temp .*= @view bfFilterLIEl[aziIndl,:,:]
  sum!(bfOutLIEl,temp)

  # OPTIONAL normalize according to number of leave in elements
  # REQUIRED when doing residual rather than corr of waveform (measured - predict)
  if normalize
    bfOutLIEl ./= size(cztWaveLIEl,2) # corrWaveLIEl
  end

  # predict time series of beam formed output of leave-in elements
  # ifft!(bfOutLIE)
  nothing
end

#