text:
module DiffEqDevTools

using DiffEqBase, RecipesBase, RecursiveArrayTools

import Base: length

const TIMESERIES_ERRORS = Set([:l2,:l∞,:L2,:L∞])
const DENSE_ERRORS = Set([:L2,:L∞])

include("benchmark.jl")
include("convergence.jl")
include("plotrecipes.jl")
include("test_solution.jl")
include("ode_tableaus.jl")
include("tableau_info.jl")

export ConvergenceSimulation, Shootout, ShootoutSet, TestSolution

#Benchmark Functions
export ode_shootout, ode_shootoutset, ode_workprecision, ode_workprecision_set

export test_convergence, appxtrue!, appxtrue

#Plot Functions
export stability_region

#Tableus
export constructEuler, constructKutta3, constructRK4, constructRK438Rule,
       constructImplicitEuler, constructMidpointRule, constructTrapezoidalRule,
       constructLobattoIIIA4, constructLobattoIIIB2, constructLobattoIIIB4,
       constructLobattoIIIC2, constructLobattoIIIC4, constructLobattoIIICStar2,
       constructLobattoIIICStar4, constructLobattoIIID2, constructLobattoIIID4,
       constructRadauIA3, constructRadauIA5, constructRadauIIA3, constructRadauIIA5,
       constructRalston, constructHeun, constructRKF5, constructBogakiShampine3,
       constructCashKarp, constructRKF8, constructDormandPrince8,
       constructMSRI1,constructFeagin10, constructFeagin12, constructFeagin14,
       constructDormandPrince8_64bit, constructRKF5, constructRungeFirst5,
       constructCassity5, constructLawson5, constructLutherKonen5, constructLutherKonen52,
       constructLutherKonen53, constructPapakostasPapaGeorgiou5, constructPapakostasPapaGeorgiou52,
       constructTsitouras5, constructBogakiShampine5, constructSharpSmart5,
       constructButcher6, constructButcher7, constructDverk, constructClassicVerner6,
       constructClassicVerner7, constructClassicVerner8, constructClassicVerner92,
       constructVernerRobust7, constructEnrightVerner7, constructTanakaYamashitaStable7,
       constructTanakaYamashitaEfficient7, constructSharpSmart7, constructSharpVerner7,
       constructVernerEfficient7, constructCooperVerner8, constructCooperVerner82,
       constructTsitourasPapakostas8, constructdverk78, constructEnrightVerner8,
       constructCurtis8, constructVernerRobust9, constructVernerEfficient9,
       constructSharp9, constructTsitouras9, constructTsitouras92,constructFeagin14Tableau,
       constructFeagin12Tableau, constructOno12, constructCurtis10, constructOno10, constructFeagin10Tableau,
       constructCurtis10, constructBaker10, constructHairer10, constructButcher63,
       constructButcher6, constructButcher62, constructVerner6, constructDormandPrince6,
       constructSharpVerner6, constructVerner9162, constructVerner916, constructVernerRobust6,
       constructVernerEfficient6, constructPapakostas6, constructLawson6,
       constructTsitourasPapakostas6, constructDormandLockyerMcCorriganPrince6,
       constructTanakaKasugaYamashitaYazaki6D, constructTanakaKasugaYamashitaYazaki6C,
       constructTanakaKasugaYamashitaYazaki6B, constructTanakaKasugaYamashitaYazaki6A,
       constructMikkawyEisa, constructChummund6, constructChummund62,
       constructHuta62, constructHuta6, constructRKF4, constructVerner7, constructVerner8,
       constructVerner9, constructVerner6

end # module

max_stars_repo_name: davidmuaddib/DiffEqDevTools

max_issues_repo_name: davidmuaddib/DiffEqDevTools

max_forks_repo_name: davidmuaddib/DiffEqDevTools


text:
"""
    mutable struct Replica{AS<:AbstractSampler}

Replica type for replica exchange Monte Carlo. A Sunny sampler type must
be provided during construction.
"""
mutable struct Replica{AS <: AbstractSampler}
    rank::Int64             # MPI rank of sampler is ∈ [0, N_ranks-1]
    N_ranks::Int64          # Total number of MPI ranks (samplers, replicas)
    N_rex_up::Int64         # Count of total number of attempted exchanges btw. replicas (i, i+1)
    rex_dir::Int64          # Replica's current exchange direction ("up"=2, "down"=1)
    α::Float64              # The current value for control parameter
    nn_ranks::Vector{Int64} # MPI ranks of nearest-neighbor samplers, -1 signals no neighbor
    label::Int64            # Label for replica control parameter tracking. Starts as N_ranks+rank+1,
                            #   and updates as swaps occur. Label is < 0 when last α is most recently
                            #   visited and > 0 when first α most recently visited
    N_up::Int64             # Number of exchanges with replica most recently visiting first α
    N_down::Int64           # Number of exchanges with replica most recently visiting last α
    sampler::AS             # Sampler containing/updating this replica's system
end

"""
Initialize MPI and return process rank, number of ranks
"""
function init_MPI()
    if !MPI.Initialized()
        MPI.Init()
    end
    return (MPI.Comm_rank(MPI.COMM_WORLD), MPI.Comm_size(MPI.COMM_WORLD))
end

"""
Constructor which initializes MPI communicator and sets sampler. -> Use on
single communicator with no groups for now
"""
function Replica(sampler::AS, α::Float64) where {AS <: AbstractSampler}
    # MPI must be initialized by user after importing MPI
    if !MPI.Initialized()
        println("""MPI not initialized. Please add "MPI.Init()" at start of script""")
        exit(-1)
    end

    rank = MPI.Comm_rank(MPI.COMM_WORLD)
    N_ranks = MPI.Comm_size(MPI.COMM_WORLD)

    Random.seed!(round(Int64, time()*1000))

    # Even rank exch. down when rex_dir==1, up when rex_dir==2
    nn_ranks = [rank-1, rank+1]

    # First and last replicas only exch. in one dir.
    for i in 1:2
        if (nn_ranks[i] < 0) || (nn_ranks[i] > N_ranks-1)
            nn_ranks[i] = -1
        end
    end

    # Start extremal ranks with up/down trip direction
    if rank == 0
        label = 1
    elseif rank == N_ranks-1
        label = -N_ranks
    else
        label = N_ranks + rank+1
    end

    # Start even ranks exchanging down
    rex_dir = (rank % 2 == 0) ? 1 : 2

    return Replica(rank, N_ranks, 0, rex_dir, α, nn_ranks, label, 0, 0, 0, sampler)
end

"""
Propose and accept/reject a replica exchange. -> Swap spin configurations
instead of the α's for now (no bookkeeping or master proc. needed). -> Use
deterministic even/odd scheme.
"""
function replica_exchange!(replica::Replica)
    # First and last replicas only exch. in one dir.
    if replica.nn_ranks[replica.rex_dir] < 0
        return false
    end

    if replica.rex_dir == 2
        replica.N_rex_up += 1
    end

    # Rank of exchange partner
    rex_rank = replica.nn_ranks[replica.rex_dir]

    # Action of this rank's configuration
    S_curr = running_energy(replica.sampler) / get_temp(replica.sampler)

    # Backup current configuration
    backup_sites = deepcopy(replica.sampler.system.sites)

    # Swap trial configuration with partner
    MPI.Sendrecv!(
                        backup_sites, rex_rank, 1,
        replica.sampler.system.sites, rex_rank, 1,
        MPI.COMM_WORLD
    )

    # This rank's action with partner's configuration
    S_exch = energy(replica.sampler.system) / get_temp(replica.sampler)

    # change in action due to exchange
    ΔS = S_curr - S_exch

    # Calculate exchange probability
    rex_accept = false
    if replica.rank < rex_rank
        rex_ΔS = MPI.Recv(Float64, rex_rank, 2, MPI.COMM_WORLD)[1]
        ln_P = ΔS + rex_ΔS

        # Replica exchange acceptance criterion
        if (ln_P >= 0) || (rand() < exp(ln_P))
            rex_accept = true
        end

        MPI.Send(rex_accept, rex_rank, 3, MPI.COMM_WORLD)
    else
        MPI.Send(ΔS, rex_rank, 2, MPI.COMM_WORLD)

        rex_accept = MPI.Recv(Bool, rex_rank, 3, MPI.COMM_WORLD)[1]
    end

    # Reject exchange
    if !rex_accept
        replica.sampler.system.sites .= backup_sites
        return false
    end

    # Recalculate energy and magnetization of new state
    reset_running_energy!(replica.sampler)
    reset_running_mag!(replica.sampler)

    return true
end

"""
Exchange labels between two replicas.

Each replica has a label that starts as (N_ranks + rank+1) until the replica
reaches an extremal temperature, after which the label ∈ [1, N_ranks]. Labels
are > 0 after most recently visiting first α and < 0 after most recently
visiting lasst α

The label trip direction (±) is used to calculate replica flow for FBOPT and
used to track all replica trajectories.

TODO: use to calculate average round-trip time for replicas
"""
function label_exchange!(replica::Replica, rex_rank::Int64)
    # Exchange labels
    rex_label = [0]
    MPI.Sendrecv!(
        [replica.label], rex_rank, 4,
        rex_label,       rex_rank, 4,
        MPI.COMM_WORLD
    )
    replica.label = rex_label[1]

    # Set label of replica as initialized (has trip direction)
    if (replica.rank == 0) || (replica.rank == replica.N_ranks-1)
        if replica.label > replica.N_ranks
            replica.label -= replica.N_ranks
        end
    end

    # Change the trip direction when α extrema are reached
    if (replica.rank == 0) && (replica.label < 0)
        replica.label *= -1
    elseif (replica.rank == replica.N_ranks-1) && (replica.label > 0)
        replica.label *= -1
    end

    return nothing
end

"""
Gather data from each process and print to specified filename
-> only handle 1D data buffer right now
"""
function gather_print(replica::Replica, data, fname::String)
    # Have root (rank 0) gather data from all processes
    data_all = MPI.Gather(data, 0, MPI.COMM_WORLD)

    # Write gathered data to file w/ column formatting
    if replica.rank == 0
        data_all = reshape(data_all, :, replica.N_ranks)

        f = open(fname, "w")
        for data_c in eachcol(data_all)
            for val in data_c
                print(f, val, "\t")
            end
            print(f, "\n")
        end
        close(f)
    end

    return nothing
end

"""
Print xyz formatted (Lx, Ly, Lz, Sx, Sy, Sz) configurations to file
"""
function xyz_to_file(sys::SpinSystem, output::IOStream)
    sites = reinterpret(reshape, Float64, sys.lattice)
    spins = reinterpret(reshape, Float64, sys.sites)
    xyz = vcat(sites, spins)
    xyz = reshape(xyz, 6, :)

    for c in eachcol(xyz)
        @printf(output, "%f\t%f\t%f\t%f\t%f\t%f\n", c...)
    end
    println(output, "")

    return nothing
end

"""
Update α schedule using feedback-optimization scheme from
Katzgraber et. al. (https://arxiv.org/pdf/cond-mat/0602085v3.pdf).

Note: use f(α) that goes from fist α → last α instead of original paper
Note: use (k-1)/(M-1) as target for integral in Eq. 11 instead of k/M
"""
function feedback_update!(replica::Replica, set_α!::F; w::Float64=0.0, dα′::Float64=1e-5) where {F <: Function}
    # Storage for collecting updated α
    _α = [replica.α]

    # Replica flow
    N_labeled = (replica.N_up + replica.N_down)
    f = (N_labeled > 0) ? replica.N_down/N_labeled : 0.0

    # Gather α and replica flow 'f' from all replicas
    all_α_and_f = MPI.Gather([replica.α, f], 0, MPI.COMM_WORLD)

    if replica.rank == 0
        all_α_and_f = reshape(all_α_and_f, :, replica.N_ranks)
        α = all_α_and_f[1, :]
        f = all_α_and_f[2, :]

        M = replica.N_ranks
        L = collect(range(0, 1.0, length=M))

        # Use 'smoothed replica flow' from Hamze et. al. (https://arxiv.org/pdf/1004.2840.pdf)
        fs = (1-w).*f .+ (w .* L)

        # Use crude linear approximation for df/dα
        Δf = fs[2:end] .- fs[1:end-1]
        Δf[Δf .< 0] .= 0.0

        Δα = α[2:end] .- α[1:end-1]

        # Calculated normalized density for adjusted α's
        η = sqrt.( (Δf ./ Δα) ./ Δα )
        C = sum(η .* Δα)
        η ./= C

        # New α's have fixed end points
        α_new = zeros(M)
        α_new[1] = α[1]
        α_new[M] = α[M]

        α′= α[1]
        cηf = 0.0
        pos = 2

        # find new α's using Eq. 11 from Katzgraber et. al.
        for i in 1:M-1
            for j in 1:Int(round(Δα[i]/dα′))
                cηf += η[i]*dα′
                α′+= dα′

                if (cηf >= L[pos]) && (pos < M)
                    α_new[pos] = α′
                    pos += 1
                end
            end
        end

        # Send out updated α to all replicas
        MPI.Scatter!(MPI.UBuffer(α_new, 1), _α, 0, MPI.COMM_WORLD)
    else
        MPI.Scatter!(nothing, _α, 0, MPI.COMM_WORLD)
    end

    # Assign new value of α and set in replica sampler
    # -> corresponds to update of kT or other Hamiltonian parameter
    set_α!(replica, _α[1])

    return nothing
end

"""
Run a feedback-optimized parallel tempering simulation and return an optimized
temperature set. No measurements are made, but this is used to diagnose replica
flow in PT simulations.

# Arguments
- `replica::Replica`: Replica type that contains MPI info and system data

- `set_α!::Function`: User-provided function for setting control parameters
  in the replica's sampler or system. Needs form: set_α!(::Replica, ::Float64)

- `therm_mcs::Int64`: Number of initial thermalization MC sweeps

- `rex_interval::Int64`: Number of MC sweeps between replica exchange attempts

- `max_mcs_opt::Int64`: Maximum number of MC sweeps to allow during optimization

- `update_interval::Int64`: Number of MC sweeps between feedback updates

- `w::Float64`: Damping factor (w ∈ [0,1]) for feedback updates: 0 gives
  aggressive updates, 1 gives leaves set of α unchanged

- `dα′::Float64`: Step size for integration during feedback update

- `print_ranks_trajectory::Vector{Float64}`: If an MPI rank is in this array, the
  timeseries for the replica starting with this rank will printed to file

- `trajectory_interval::Int64`: Number of replica exchange attempts between printed
  timeseries points

"""
function run_FBO!(
    replica::Replica,
    set_α!::F;
    therm_mcs::Int64=1000,
    rex_interval::Int64=1,
    max_mcs::Int64=500_000,
    update_interval::Int64=100_000,
    w::Float64=0.0,
    dα′::Float64=1e-5,
    print_ranks_trajectory::Vector{Int64}=Int64[],
    trajectory_interval::Int64=500_000
) where {F <: Function}

    # Initialize energy and equilibrate replica to it's distribution
    reset_running_energy!(replica.sampler)
    reset_running_mag!(replica.sampler)

    thermalize!(replica.sampler, therm_mcs)

    replica.sampler.nsweeps = 1

    rex_accepts = zeros(Int64, 2)
    N_updates = 0
    N_rex = 0

    # Start feedback optimization simulation
    for mcs in 1:max_mcs
        sample!(replica.sampler)

        if (replica.rank == 0) && (mcs % 1000 == 0)
            println("FBO ", mcs, " mcs")
        end

        # Attempt replica exchange
        if mcs % rex_interval == 0
            N_rex += 1

            if replica_exchange!(replica)
                rex_accepts[replica.rex_dir] += 1

                # Exchange labels
                label_exchange!(replica, replica.nn_ranks[replica.rex_dir])
            end

            # Record replica flow
            if replica.nn_ranks[replica.rex_dir] != -1
                if replica.label < 0
                    replica.N_down += 1
                elseif replica.label < replica.N_ranks+1
                    replica.N_up += 1
                end
            end

            # Print out replica exchange timeseries -- inefficient IO
            if (N_rex % trajectory_interval == 0) && (abs(replica.label)-1 in print_ranks_trajectory)
                if replica.label <= replica.N_ranks
                    fname = @sprintf("FBO-update%03d_trajectory%03d.dat", N_updates, abs(replica.label)-1)
                    output = open(fname, "a")
                    println(output, replica.rank)
                    close(output)
                end
            end

            # Alternate up/down pairs of replicas for exchanges
            replica.rex_dir = 3 - replica.rex_dir
        end

        if mcs % update_interval == 0
            # Replica flow 'f' and smoothed replica flow 'fs'
            N_labeled = (replica.N_up + replica.N_down)
            f = (N_labeled > 0) ? replica.N_down/N_labeled : 0.0
            fs = (1-w)*f + w*replica.rank/(replica.N_ranks-1)

            # Acceptance rate between replicas i and i+1
            A = (replica.rank == replica.N_ranks-1) ? 0.0 : rex_accepts[2]/replica.N_rex_up

            # Print out α, replica flow, and acceptance rates for each update
            fname = @sprintf("FBO-update%03d.dat", N_updates)
            gather_print(
                replica,
                [replica.α, fs, A],
                fname
            )

            # Perform feedback optimization for control parameters
            feedback_update!(replica, set_α!; w=w, dα′=dα′)

            N_updates += 1
            rex_accepts[1] = rex_accepts[2] = 0
            reset_stats!(replica)
        end
    end

    rex_accepts[1] = rex_accepts[2] = 0
    reset_stats!(replica)

    # Return a vector of optimized α's
    return MPI.Allgather([replica.α], MPI.COMM_WORLD)
end

"""
Reset the exchange statistics and labels stored in replica
"""
function reset_stats!(replica::Replica)
    replica.N_up = replica.N_down = 0
    replica.N_rex_up = 0

    # Reset label
    if replica.rank == 0
        replica.label = 1
    elseif replica.rank == replica.N_ranks-1
        replica.label = -replica.N_ranks
    else
        replica.label = replica.N_ranks + replica.rank+1
    end
    return nothing
end

"""
Start a parallel tempering (PT) simulation that saves measured averages,
histograms, minimal energies, and configurations.

Run w/ the command "mpiexec -n [n_procs] julia --project [julia_script.jl]".

# Arguments
- `replica::Replica`: Replica type that contains MPI info and system data

- `therm_mcs::Int64`: Number of initial thermalization MC sweeps

- `measure_interval::Int64`: Number of MC sweeps between measurements

- `rex_interval::Int64`: Number of MC sweeps between replica exchange attempts

- `max_mcs::Int64`: Maximum number of MC sweeps to allow during simulation

- `bin_size::Float64`: Width in energy space for each histogram bin

- `print_hist::Bool`: Print energy histograms to file for each rank if true

- `print_xyz_ranks::Vector{Float64}`: If an MPI rank is in this array, the xyz
  coordinates for new minimal energies are printed to file

- `print_ranks_trajectory::Vector{Float64}`: If an MPI rank is in this array, the
  timeseries for the replica starting with this rank will printed to file

- `trajectory_interval::Int64`: Number of replica exchange attempts between printed
  timeseries points

- `measure_replica_flow::Bool`: If true, then the flux of downward replica flow in
  α space is measured and printed after the simulation

"""
function run_REMC!(
    replica::Replica;
    therm_mcs::Int64=1000,
    measure_interval::Int64=10,
    rex_interval::Int64=1,
    max_mcs::Int64=1_000_000,
    bin_size::Float64=1.0,
    print_hist::Bool=false,
    print_xyz_ranks::Vector{Int64}=Int64[],
    print_ranks_trajectory::Vector{Int64}=Int64[],
    trajectory_interval::Int64=1_000_000,
    measure_replica_flow::Bool=false
)
    # Initialize energy and equilibrate replica to it's distribution
    reset_running_energy!(replica.sampler)
    reset_running_mag!(replica.sampler)

    thermalize!(replica.sampler, therm_mcs)

    replica.sampler.nsweeps = 1

    # Manually include these basic measurements for now
    U  = 0.0
    U² = 0.0
    M⃗  = [0.0, 0.0, 0.0]
    M  = 0.0
    M² = 0.0
    C  = 0.0
    X  = 0.0

    system_size = length(replica.sampler.system)
    Emin = typemax(Float64)
    E_hist = BinnedArray{Float64, Int64}(bin_size=bin_size)
    m_hist = BinnedArray{Float64, Int64}(bin_size=1.0)
    rex_accepts = zeros(Int64, 2)
    N_measure = 0
    N_rex = 0

    # Start PT with finite length
    for mcs in 1:max_mcs
        sample!(replica.sampler)

        # Crude progress printing to stdout
        if (replica.rank == 0) && (mcs % 1000 == 0)
            println("REMC ", mcs, " mcs")
        end

        # Attempt replica exchange
        if mcs % rex_interval == 0
            N_rex += 1

            if replica_exchange!(replica)
                rex_accepts[replica.rex_dir] += 1

                # Exchange labels
                if measure_replica_flow
                    label_exchange!(replica, replica.nn_ranks[replica.rex_dir])
                end
            end

            # Record replica flow (for quality evaluation)
            if measure_replica_flow && replica.nn_ranks[replica.rex_dir] != -1
                if replica.label < 0
                    replica.N_down += 1
                elseif replica.label < replica.N_ranks+1
                    replica.N_up += 1
                end
            end

            # Print out replica exchange timeseries -- inefficient IO
            if (N_rex % trajectory_interval == 0) && (abs(replica.label)-1 in print_ranks_trajectory)
                if replica.label <= replica.N_ranks
                    fname = @sprintf("trajectory%03d.dat", abs(replica.label)-1)
                    rex_output = open(fname, "a")
                    println(rex_output, replica.rank)
                    close(rex_output)
                end
            end

            # Alternate up/down pairs of replicas for exchanges
            replica.rex_dir = 3 - replica.rex_dir
        end

        # Measurements
        if mcs % measure_interval == 0
            E = running_energy(replica.sampler)
            m⃗ = running_mag(replica.sampler)
            m = norm(m⃗)

            # New minimum energy measured
            if E < Emin
                Emin = E

                # Print out evolution of Emin w.r.t MC time
                f = open(@sprintf("Emin%03d.dat", replica.rank), "a")
                println(f, mcs," ",Emin)
                close(f)

                # Overwrite Emin coordinates in saved file
                if replica.rank in print_xyz_ranks
                    xyz_output = open(@sprintf("Emin%03d.xyz", replica.rank), "w")
                    xyz_to_file(replica.sampler.system, xyz_output)
                    close(xyz_output)
                end
            end

            U  += E
            U² += E^2
            M⃗ .+= m⃗
            M  += m
            M² += m^2
            N_measure += 1

            # Histograms for energy and magnetization
            E_hist[E] += 1
            m_hist[m] += 1
        end
    end

    # Normalize averages
    kT = get_temp(replica.sampler)
    U  /= N_measure
    U² /= N_measure
    M⃗ ./= N_measure
    M  /= N_measure
    M² /= N_measure
    C = 1.0/kT^2 * (U² - U^2)
    X = 1.0/kT   * (M² - M^2)

    # Acceptance rate between replicas i and i+1
    A = (replica.rank == replica.N_ranks-1) ? 0.0 : rex_accepts[2]/replica.N_rex_up

    # Gather and print: rank, "down" exch. accepts, "up" exch. accepts, acceptance ratio
    gather_print(
        replica,
        [replica.α, A],
        "replica_exchanges.dat"
    )

    # Gather and print thermodynamic measurements
    gather_print(
        replica,
        [replica.α, (M⃗ ./ system_size)..., M/system_size, X/system_size, U/system_size, C/system_size],
        "measurements.dat"
    )

    # Write histogram to file for each process if specified
    if print_hist
        f = open(@sprintf("E-hist%03d.dat", replica.rank), "w")
        println(f, E_hist)
        close(f)
        f = open(@sprintf("m-hist%03d.dat", replica.rank), "w")
        println(f, m_hist)
        close(f)
    end

    # Print replica flow 'f' to file
    if measure_replica_flow
        N_labeled = (replica.N_up + replica.N_down)
        f = (N_labeled > 0) ? replica.N_down/N_labeled : 0.0

        gather_print(
            replica,
            [replica.α, f],
            "replica_flow.dat"
        )
    end

    return :SUCCESS
end

ext: jl


text:
function generateCorrRandom(rhoz)

        e = randn(size(rhoz, 2))

        C = transpose(chol(rhoz))

        *(C,e)

end


text:
struct Foo
    X::Vector{Int}
end

object = Foo([7, 7, 7])

object.X

Y = object.X
(object.X)[1] = 8

object.X

Y

mutable struct Foo
    x::Vector{Float64}
end

function bar(object::Foo)

    object.x[1]=2.3
    mu = mean(object.x)

    return object, mu
end

function bar2(object::Foo)
    x = object.x

    x[1]=2.3
    mu = mean(x)

    object.x = x

    return object, mean(x)
end

object = Foo(randn(10^6))

@allocated bar(object)
@allocated bar2(object)


text:
using TestDatasets
using Test

@testset "TestDatasets.jl" begin
    # Write your tests here.
end

max_stars_repo_name: dvillacis/TestDatasets

max_issues_repo_name: dvillacis/TestDatasets

max_forks_repo_name: dvillacis/TestDatasets


text:
##############################################################################################################
# Optimization Run for NICE With Revenue Recycling
##############################################################################################################

#---------------------------------------------------------------------
# Set up first-stage global optimization for NICE + Revenue Recycling.
#---------------------------------------------------------------------

println("Beginning global optimization for NICE + Revenue Recycling.")

# Create NICE + Revenue Recycling objective function with user-specified settings.
recycle_objective = construct_nice_recycle_objective(nice, revenue_recycling = true)

# Create an NLopt optimization object.
recycle_opt_global = Opt(global_opt_algorithm, n_objectives)

# Set upper and lower bounds for global CO₂ tax (upper bound = maximum RICE backstop prices).
upper_bound = maximum(rice_backstop_prices, dims=2)[2:(n_objectives+1)]

lower_bounds!(recycle_opt_global, zeros(n_objectives))
upper_bounds!(recycle_opt_global, upper_bound)

# Set maximum run time.
maxtime!(recycle_opt_global, global_stop_time)

# Set convergence tolerance.
ftol_rel!(recycle_opt_global, global_tolerance)

# Set objective function.
max_objective!(recycle_opt_global, (x, grad) -> recycle_objective(x))

# Optimize model using global optimization algorithm, initialized at a fraction of the upper bound.
recycle_max_welfare_global, recycle_opt_tax_global, recycle_convergence_global = optimize(recycle_opt_global, upper_bound ./ 100)

println("Convergence result = ", recycle_convergence_global)

# Calculate full tax time series (just to have to cross-check with local solution).
recycle_full_global_opt_tax, recycle_global_opt_regional_mitigation = mu_from_tax(recycle_opt_tax_global, rice_backstop_prices)

#----------------------------------------------------------------------
# Set up second-stage local optimization for NICE + Revenue Recycling.
#----------------------------------------------------------------------

println("Beginning local optimization for NICE + Revenue Recycling.")

# Create an NLopt optimization object for second-stage local optimization.
recycle_opt = Opt(local_opt_algorithm, n_objectives)

# Set upper and lower bounds for global CO₂ tax (upper bound = maximum RICE backstop prices).
lower_bounds!(recycle_opt, zeros(n_objectives))
upper_bounds!(recycle_opt, upper_bound)

# Set maximum run time.
maxtime!(recycle_opt, local_stop_time)

# Set convergence tolerance.
ftol_rel!(recycle_opt, local_tolerance)

# Set objective function.
max_objective!(recycle_opt, (x, grad) -> recycle_objective(x))

# Optimize model using local optimization algorithm, initialized at global optimization solution.
recycle_max_welfare, recycle_opt_tax, recycle_convergence = optimize(recycle_opt, recycle_opt_tax_global)

println("Convergence result = ", recycle_convergence)

#---------------------------------------------------
# Run Model Under Optimal Policy and Save Results
#---------------------------------------------------

# Get optimal regional CO₂ mitigation rates from final optimal tax.
recycle_full_opt_tax, recycle_opt_regional_mitigation = mu_from_tax(recycle_opt_tax, rice_backstop_prices)

# Run model under optimal policy.
update_param!(nice, :MIU, recycle_opt_regional_mitigation)
update_param!(nice, :global_carbon_tax, recycle_full_opt_tax)
run(nice)

# Save results.
if save_results == true
    save_nice_recycle_results(nice, bau_model, recycle_full_opt_tax, recycle_full_global_opt_tax, output_directory, revenue_recycling=true)
end

println("Recycle runs complete.")


text:
module Compute
include("GaussianSlab.jl")
include("MajorEllipticTorus.jl")
include("DoubleEllipticTorus.jl")
include("EllipticTube.jl")
include("Tube.jl")
include("Rod.jl")
include("Rods.jl")
include("Sphere.jl")
include("Plane.jl")

function get_fn(shape::AbstractString, params)
    shapes = Dict(
        :GaussianSlab => GaussianSlab.get_fn,
        :MajorEllipticTorus => MajorEllipticTorus.get_fn,
        :DoubleEllipticTorus => DoubleEllipticTorus.get_fn,
        :EllipticTube => EllipticTube.get_fn,
        :Tube => Tube.get_fn,
        :Rod => Rod.get_fn,
        :Rods => Rods.get_fn,
        :Sphere => Sphere.get_fn,
        :Plane => Plane.get_fn,
    )
    return shapes[Symbol(shape)](params)
end
export get_fn
end

max_stars_repo_name: yuhangwang/SimShape

max_issues_repo_name: yuhangwang/SimShape

max_forks_repo_name: yuhangwang/SimShape


text:
module FiniteVolumes

using ForwardDiff
using LinearAlgebra
using StaticArrays

# CORE

include("./flux_function.jl")
export LinearAdvectionFlux, AdvectionFlux, RotationFlux, ShallowWater

include("./scheme.jl")
export Upwind, Centered, NeumannBC

include("./courant.jl")

# EXAMPLE SUBMODULE
# Cartesian meshes
include("./cartesian_mesh/CartesianMeshes.jl")
using FiniteVolumes.CartesianMeshes
export CartesianMesh, PeriodicCartesianMesh, cell_centers

include("./splitting.jl")
export directional_splitting

# Time solver
include("./euler_explicit.jl")
export FixedCourant, TimeDependantFlux

# EXPERIMENTS
include("./diffusion.jl")
export DiffusionFlux

include("./experimental.jl")

end


text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     comment_magics: false
#     formats: ipynb,jl:light
#     text_representation:
#       extension: .jl
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.7.1
#   kernelspec:
#     display_name: Julia 1.5.3
#     language: julia
#     name: julia-1.5
# ---

# # Maxwell solver FEM 1D

using GEMPIC, LinearAlgebra, Plots

const mode = 2

eta1_min = .0
eta1_max = 2π
nc_eta1  = 128
Lx = eta1_max - eta1_min
delta_eta1 = Lx / nc_eta1
deg = 3

mesh = Mesh( eta1_min , eta1_max, nc_eta1)
maxwell_1d = Maxwell1DFEM(mesh, deg);

ex = zeros(Float64, nc_eta1)
ey = zeros(Float64, nc_eta1)
bz = zeros(Float64, nc_eta1);

bz_exact = zeros(Float64, nc_eta1)
ex_exact = zeros(Float64, nc_eta1)
ey_exact = zeros(Float64, nc_eta1)
rho      = zeros(Float64, nc_eta1)
sval     = zeros(Float64, nc_eta1);

cos_k(x) = cos(mode*2*pi*x/Lx)
sin_k(x) = sin(mode*2*pi*x/Lx)

# Test Poisson
#-------------
# Set exact solution
for i = 1:nc_eta1
   xi = eta1_min + (i-1)*delta_eta1
   ex_exact[i] = sin_k(xi)/(2.0*mode*pi/Lx)
end

x = LinRange(eta1_min, eta1_max, nc_eta1+1)[1:end-1]

compute_rhs_from_function!( rho, maxwell_1d, cos_k, deg)

compute_e_from_rho!( ex, maxwell_1d, rho )

plot(x, ex_exact)
sval = eval_uniform_periodic_spline_curve(deg-1, ex)
scatter!(x, sval)

dt = .5 * delta_eta1

for i = 1:nc_eta1
   xi = eta1_min + (i-1)*delta_eta1
   ex_exact[i] = - cos_k(xi)*dt
end

compute_rhs_from_function!(rho, maxwell_1d, cos_k, deg-1)
fill!(ex, 0.0)
compute_e_from_j!(ex, maxwell_1d, dt .* rho, 1 )
sval =  eval_uniform_periodic_spline_curve(deg-1, ex);

plot(x, ex_exact)
scatter!(x, sval)

# +
mode     = 2
eta1_min = .0
eta1_max = 2π
nc_eta1  = 256

Lx = eta1_max - eta1_min

delta_eta1 = Lx / nc_eta1

deg = 3

maxwell_1d = Maxwell1DFEM(mesh, deg)

ex = zeros(Float64, nc_eta1)
ey = zeros(Float64, nc_eta1)
bz = zeros(Float64, nc_eta1)

bz_exact = zeros(Float64, nc_eta1)
ex_exact = zeros(Float64, nc_eta1)
ey_exact = zeros(Float64, nc_eta1)
rho      = zeros(Float64, nc_eta1)
sval     = zeros(Float64, nc_eta1)

cos_k(x) = cos(mode*2*pi*x/Lx)
sin_k(x) = sin(mode*2*pi*x/Lx)

# +
# Test Maxwell on By and Ez
#--------------------------
# Set time stepping parameters

time  = 0.0
dt    = .5 * delta_eta1
nstep = 10

# Compute initial fields
ex = 0.0 # 0-form -> splines of degree deg
l2projection!(bz, maxwell_1d, cos_k, deg-1) # 0-form -> splines of degree deg-1

#plot(bz)

# +

compute_b_from_e!( bz, maxwell_1d, 0.5*dt, ey)

plot(bz)
# -

compute_e_from_b!( ey, maxwell_1d,     dt, bz)
plot(ey)

# +
compute_b_from_e!( bz, maxwell_1d, 0.5*dt, ey)

time = time + dt

for i = 1:nc_eta1
   xi = eta1_min + (i-1)*delta_eta1
   ey_exact[i] = sin(mode * 2π * xi/Lx) * sin(mode * 2π * time/Lx)
   bz_exact[i] = cos(mode * 2π * xi/Lx) * cos(mode * 2π * time/Lx)
end

sval = eval_uniform_periodic_spline_curve(deg, ey)
@show err_ey = norm(sval .- ey_exact)

sval = eval_uniform_periodic_spline_curve(deg-1, bz)
@show err_bz = norm(sval .- bz_exact)
# -


text:
export OSProc, Context, addprocs!, rmprocs!

"""
    Processor

An abstract type representing a processing device and associated memory, where
data can be stored and operated on. Subtypes should be immutable, and
instances should compare equal if they represent the same logical processing
device/memory. Subtype instances should be serializable between different
nodes. Subtype instances may contain a pointer to a "parent" `Processor` to
make it easy to transfer data to/from other types of `Processor` at runtime.
"""
abstract type Processor end

const PROCESSOR_CALLBACKS = []

"""
    execute!(proc::Processor, f, args...) -> Any

Executes the function `f` with arguments `args` on processor `proc`. This
function can be overloaded by `Processor` subtypes to allow executing function
calls differently than normal Julia.
"""
execute!

"""
    iscompatible(proc::Processor, opts, f, Targs...) -> Bool

Indicates whether `proc` can execute `f` over `Targs` given `opts`. `Processor`
subtypes should overload this function to return `true` if and only if it is
essentially guaranteed that `f(::Targs...)` is supported. Additionally,
`iscompatible_func` and `iscompatible_arg` can be overriden to determine
compatibility of `f` and `Targs` individually. The default implementation
returns `false`.
"""
iscompatible(proc::Processor, opts, f, Targs...) =
    iscompatible_func(proc, opts, f) &&
    all(x->iscompatible_arg(proc, opts, x), Targs)
iscompatible_func(proc::Processor, opts, f) = false
iscompatible_arg(proc::Processor, opts, x) = false

"""
    default_enabled(proc::Processor) -> Bool

Returns whether processor `proc` is enabled by default (opt-out). `Processor` subtypes can override this function to make themselves opt-in (default returns `false`).
"""
default_enabled(proc::Processor) = false

"""
    get_processors(proc::Processor) -> Vector{T} where T<:Processor

Returns the full list of processors contained in `proc`, if any. `Processor`
subtypes should overload this function if they can contain sub-processors. The
default method will return a `Vector` containing `proc` itself.
"""
get_processors(proc::Processor) = Processor[proc]

"""
    get_parent(proc::Processor) -> Processor

Returns the parent processor for `proc`. The ultimate parent processor is an
`OSProc`. `Processor` subtypes should overload this to return their most
direct parent.
"""
get_parent

"""
    move(from_proc::Processor, to_proc::Processor, x)

Moves and/or converts `x` such that it's available and suitable for usage on
the `to_proc` processor. This function can be overloaded by `Processor`
subtypes to transport arguments and convert them to an appropriate form before
being used for exection. Subtypes of `Processor` wishing to implement efficient
data movement should provide implementations where `x::Chunk`.
"""
move(from_proc::Processor, to_proc::Processor, x) = x

"""
    capacity(proc::Processor=OSProc()) -> Int

Returns the total processing capacity of `proc`.
"""
capacity(proc=OSProc()) = length(get_processors(proc))
capacity(proc, ::Type{T}) where T =
    length(filter(x->x isa T, get_processors(proc)))

"""
    OSProc <: Processor

Julia CPU (OS) process, identified by Distributed pid. Executes thunks when
threads or other "children" processors are not available, and/or the user has
not opted to use those processors.
"""
struct OSProc <: Processor
    pid::Int
    attrs::Dict{Symbol,Any}
    children::Vector{Processor}
    queue::Vector{Processor}
end
const OSPROC_CACHE = Dict{Int,OSProc}()
OSProc(pid::Int=myid()) = get!(OSPROC_CACHE, pid) do
    remotecall_fetch(get_osproc, pid, pid)
end
function get_osproc(pid::Int)
    proc = OSProc(pid, Dict{Symbol,Any}(), Processor[], Processor[])
    for cb in PROCESSOR_CALLBACKS
        try
            child = Base.invokelatest(cb, proc)
            child !== nothing && push!(proc.children, child)
        catch err
            @error "Error in processor callback" exception=(err,catch_backtrace())
        end
    end
    proc
end
function add_callback!(func)
    push!(Dagger.PROCESSOR_CALLBACKS, func)
    empty!(OSPROC_CACHE)
end
Base.:(==)(proc1::OSProc, proc2::OSProc) = proc1.pid == proc2.pid
iscompatible(proc::OSProc, opts, f, args...) =
    any(child->iscompatible(child, opts, f, args...), proc.children)
iscompatible_func(proc::OSProc, opts, f) =
    any(child->iscompatible_func(child, opts, f), proc.children)
iscompatible_arg(proc::OSProc, opts, args...) =
    any(child->
        all(arg->iscompatible_arg(child, opts, arg), args),
    proc.children)
get_processors(proc::OSProc) =
    vcat((get_processors(child) for child in proc.children)...,)
function choose_processor(options, f, Targs)
    osproc = OSProc()
    if isempty(osproc.queue)
        for child in osproc.children
            grandchildren = get_processors(child)
            append!(osproc.queue, grandchildren)
        end
    end
    @assert !isempty(osproc.queue)
    for i in 1:length(osproc.queue)
        proc = popfirst!(osproc.queue)
        push!(osproc.queue, proc)
        if !iscompatible(proc, options, f, Targs...)
            continue
        end
        if options.proclist === nothing
            default_enabled(proc) && return proc
        elseif options.proclist isa Function
            options.proclist(proc) && return proc
        elseif any(p->proc isa p, options.proclist)
            return proc
        end
    end
    throw(ProcessorSelectionException(options.proclist, osproc.queue, f, Targs))
end
struct ProcessorSelectionException <: Exception
    proclist
    procsavail::Vector{Processor}
    f
    args
end
function Base.show(io::IO, pex::ProcessorSelectionException)
    println(io, "(Worker $(myid())) Exhausted all available processor types!")
    println(io, "  Proclist: $(pex.proclist)")
    println(io, "  Procs Available: $(pex.procsavail)")
    println(io, "  Function: $(pex.f)")
    print(io, "  Arguments: $(pex.args)")
end

execute!(proc::OSProc, f, args...) = f(args...)
default_enabled(proc::OSProc) = true

"""
    ThreadProc <: Processor

Julia CPU (OS) thread, identified by Julia thread ID.
"""
struct ThreadProc <: Processor
    owner::Int
    tid::Int
end
iscompatible(proc::ThreadProc, opts, f, args...) = true
iscompatible_func(proc::ThreadProc, opts, f) = true
iscompatible_arg(proc::ThreadProc, opts, x) = true
@static if VERSION >= v"1.3.0-DEV.573"
    function execute!(proc::ThreadProc, f, args...)
        sch_handle = task_local_storage(:sch_handle)
        task = Threads.@spawn begin
            task_local_storage(:processor, proc)
            task_local_storage(:sch_handle, sch_handle)
            f(args...)
        end
        try
            fetch(task)
        catch err
            @static if VERSION >= v"1.1"
                stk = Base.catch_stack(task)
                err, frames = stk[1]
                rethrow(CapturedException(err, frames))
            else
                rethrow(task.result)
            end
        end
    end
else
    # TODO: Use Threads.@threads?
    function execute!(proc::ThreadProc, f, args...)
        sch_handle = task_local_storage(:sch_handle)
        task = @async begin
            task_local_storage(:processor, proc)
            task_local_storage(:sch_handle, sch_handle)
            f(args...)
        end
        try
            fetch(task)
        catch err
            @static if VERSION >= v"1.1"
                stk = Base.catch_stack(task)
                err, frames = stk[1]
                rethrow(CapturedException(err, frames))
            else
                rethrow(task.result)
            end
        end
    end
end
get_parent(proc::ThreadProc) = OSProc(proc.owner)
default_enabled(proc::ThreadProc) = true

# TODO: ThreadGroupProc?

"A context represents a set of processors to use for an operation."
mutable struct Context
    procs::Vector{Processor}
    proc_lock::ReentrantLock
    log_sink::Any
    log_file::Union{String,Nothing}
    profile::Bool
    options
end

"""
    Context(;nthreads=Threads.nthreads()) -> Context
    Context(xs::Vector{OSProc}) -> Context
    Context(xs::Vector{Int}) -> Context

Create a Context, by default adding each available worker once from
every available thread. Use the `nthreads` keyword to use a different
number of threads.

It is also possible to create a Context from a vector of [`OSProc`](@ref),
or equivalently the underlying process ids can also be passed directly
as a `Vector{Int}`.

Special fields include:
- 'log_sink': A log sink object to use, if any.
- `log_file::Union{String,Nothing}`: Path to logfile. If specified, at
scheduler termination logs will be collected, combined with input thunks, and
written out in DOT format to this location.
- `profile::Bool`: Whether or not to perform profiling with Profile stdlib.
"""
Context(procs::Vector{P}=Processor[OSProc(w) for w in workers()];
        proc_lock=ReentrantLock(), log_sink=NoOpLog(), log_file=nothing,
        profile=false, options=nothing) where {P<:Processor} =
    Context(procs, proc_lock, log_sink, log_file, profile, options)
Context(xs::Vector{Int}) = Context(map(OSProc, xs))
procs(ctx::Context) = lock(ctx) do
    copy(ctx.procs)
end

"""
    write_event(ctx::Context, event::Event)

Write a log event
"""
function write_event(ctx::Context, event::Event)
    write_event(ctx.log_sink, event)
end

"""
    lock(f, ctx::Context)

Acquire `ctx.proc_lock`, execute `f` with the lock held, and release the lock when `f` returns.
"""
Base.lock(f, ctx::Context) = lock(f, ctx.proc_lock)

"""
    addprocs!(ctx::Context, xs)

Add new workers `xs` to `ctx`.

Workers will typically be assigned new tasks in the next scheduling iteration if scheduling is ongoing.

Workers can be either `Processor`s or the underlying process IDs as `Integer`s.
"""
addprocs!(ctx::Context, xs::AbstractVector{<:Integer}) = addprocs!(ctx, map(OSProc, xs))
addprocs!(ctx::Context, xs::AbstractVector{<:OSProc}) = lock(ctx) do
    append!(ctx.procs, xs)
end

"""
    rmprocs!(ctx::Context, xs)

Remove the specified workers `xs` from `ctx`.

Workers will typically finish all their assigned tasks if scheduling is ongoing but will not be assigned new tasks after removal.

Workers can be either `Processor`s or the underlying process IDs as `Integer`s.
"""
rmprocs!(ctx::Context, xs::AbstractVector{<:Integer}) = rmprocs!(ctx, map(OSProc, xs))
rmprocs!(ctx::Context, xs::AbstractVector{<:OSProc}) = lock(ctx) do
    filter!(p -> (p ∉ xs), ctx.procs)
end

"Gets the current processor executing the current thunk."
thunk_processor() = task_local_storage(:processor)::Processor


text:
function abs(::Type{FastMath}, arg0::RealFieldElement)
    return jcall(FastMath, "abs", RealFieldElement, (RealFieldElement,), arg0)
end

function abs(::Type{FastMath}, arg0::jdouble)
    return jcall(FastMath, "abs", jdouble, (jdouble,), arg0)
end

function abs(::Type{FastMath}, arg0::jfloat)
    return jcall(FastMath, "abs", jfloat, (jfloat,), arg0)
end

function abs(::Type{FastMath}, arg0::jint)
    return jcall(FastMath, "abs", jint, (jint,), arg0)
end

function abs(::Type{FastMath}, arg0::jlong)
    return jcall(FastMath, "abs", jlong, (jlong,), arg0)
end

function acos(::Type{FastMath}, arg0::RealFieldElement)
    return jcall(FastMath, "acos", RealFieldElement, (RealFieldElement,), arg0)
end

function acos(::Type{FastMath}, arg0::jdouble)
    return jcall(FastMath, "acos", jdouble, (jdouble,), arg0)
end

function acosh(::Type{FastMath}, arg0::RealFieldElement)
    return jcall(FastMath, "acosh", RealFieldElement, (RealFieldElement,), arg0)
end

function acosh(::Type{FastMath}, arg0::jdouble)
    return jcall(FastMath, "acosh", jdouble, (jdouble,), arg0)
end

function add_exact(::Type{FastMath}, arg0::jint, arg1::jint)
    return jcall(FastMath, "addExact", jint, (jint, jint), arg0, arg1)
end

function add_exact(::Type{FastMath}, arg0::jlong, arg1::jlong)
    return jcall(FastMath, "addExact", jlong, (jlong, jlong), arg0, arg1)
end

function asin(::Type{FastMath}, arg0::RealFieldElement)
    return jcall(FastMath, "asin", RealFieldElement, (RealFieldElement,), arg0)
end

function asin(::Type{FastMath}, arg0::jdouble)
    return jcall(FastMath, "asin", jdouble, (jdouble,), arg0)
end

function asinh(::Type{FastMath}, arg0::RealFieldElement)
    return jcall(FastMath, "asinh", RealFieldElement, (RealFieldElement,), arg0)
end

function asinh(::Type{FastMath}, arg0::jdouble)
    return jcall(FastMath, "asinh", jdouble, (jdouble,), arg0)
end

function atan(::Type{FastMath}, arg0::RealFieldElement)
    return jcall(FastMath, "atan", RealFieldElement, (RealFieldElement,), arg0)
end

function atan(::Type{FastMath}, arg0::jdouble)
    return jcall(FastMath, "atan", jdouble, (jdouble,), arg0)
end

function atan2(::Type{FastMath}, arg0::RealFieldElement, arg1::RealFieldElement)
    return jcall(FastMath, "atan2", RealFieldElement, (RealFieldElement, RealFieldElement), arg0, arg1)
end

function atan2(::Type{FastMath}, arg0::jdouble, arg1::jdouble)
    return jcall(FastMath, "atan2", jdouble, (jdouble, jdouble), arg0, arg1)
end

function atanh(::Type{FastMath}, arg0::RealFieldElement)
    return jcall(FastMath, "atanh", RealFieldElement, (RealFieldElement,), arg0)
end

function atanh(::Type{FastMath}, arg0::jdouble)
    return jcall(FastMath, "atanh", jdouble, (jdouble,), arg0)
end

function cbrt(::Type{FastMath}, arg0::RealFieldElement)
    return jcall(FastMath, "cbrt", RealFieldElement, (RealFieldElement,), arg0)
end

function cbrt(::Type{FastMath}, arg0::jdouble)
    return jcall(FastMath, "cbrt", jdouble, (jdouble,), arg0)
end

function ceil(::Type{FastMath}, arg0::RealFieldElement)
    return jcall(FastMath, "ceil", RealFieldElement, (RealFieldElement,), arg0)
end

function ceil(::Type{FastMath}, arg0::jdouble)
    return jcall(FastMath, "ceil", jdouble, (jdouble,), arg0)
end

function copy_sign(::Type{FastMath}, arg0::RealFieldElement, arg1::RealFieldElement)
    return jcall(FastMath, "copySign", RealFieldElement, (RealFieldElement, RealFieldElement), arg0, arg1)
end

function copy_sign(::Type{FastMath}, arg0::RealFieldElement, arg1::jdouble)
    return jcall(FastMath, "copySign", RealFieldElement, (RealFieldElement, jdouble), arg0, arg1)
end

function copy_sign(::Type{FastMath}, arg0::jdouble, arg1::jdouble)
    return jcall(FastMath, "copySign", jdouble, (jdouble, jdouble), arg0, arg1)
end

function copy_sign(::Type{FastMath}, arg0::jfloat, arg1::jfloat)
    return jcall(FastMath, "copySign", jfloat, (jfloat, jfloat), arg0, arg1)
end

function cos(::Type{FastMath}, arg0::RealFieldElement)
    return jcall(FastMath, "cos", RealFieldElement, (RealFieldElement,), arg0)
end

function cos(::Type{FastMath}, arg0::jdouble)
    return jcall(FastMath, "cos", jdouble, (jdouble,), arg0)
end

function cosh(::Type{FastMath}, arg0::RealFieldElement)
    return jcall(FastMath, "cosh", RealFieldElement, (RealFieldElement,), arg0)
end

function cosh(::Type{FastMath}, arg0::jdouble)
    return jcall(FastMath, "cosh", jdouble, (jdouble,), arg0)
end

function decrement_exact(::Type{FastMath}, arg0::jint)
    return jcall(FastMath, "decrementExact", jint, (jint,), arg0)
end

function decrement_exact(::Type{FastMath}, arg0::jlong)
    return jcall(FastMath, "decrementExact", jlong, (jlong,), arg0)
end

function exp(::Type{FastMath}, arg0::RealFieldElement)
    return jcall(FastMath, "exp", RealFieldElement, (RealFieldElement,), arg0)
end

function exp(::Type{FastMath}, arg0::jdouble)
    return jcall(FastMath, "exp", jdouble, (jdouble,), arg0)
end

function expm1(::Type{FastMath}, arg0::RealFieldElement)
    return jcall(FastMath, "expm1", RealFieldElement, (RealFieldElement,), arg0)
end

function expm1(::Type{FastMath}, arg0::jdouble)
    return jcall(FastMath, "expm1", jdouble, (jdouble,), arg0)
end

function floor(::Type{FastMath}, arg0::RealFieldElement)
    return jcall(FastMath, "floor", RealFieldElement, (RealFieldElement,), arg0)
end

function floor(::Type{FastMath}, arg0::jdouble)
    return jcall(FastMath, "floor", jdouble, (jdouble,), arg0)
end

function floor_div(::Type{FastMath}, arg0::jint, arg1::jint)
    return jcall(FastMath, "floorDiv", jint, (jint, jint), arg0, arg1)
end

function floor_div(::Type{FastMath}, arg0::jlong, arg1::jint)
    return jcall(FastMath, "floorDiv", jlong, (jlong, jint), arg0, arg1)
end

function floor_div(::Type{FastMath}, arg0::jlong, arg1::jlong)
    return jcall(FastMath, "floorDiv", jlong, (jlong, jlong), arg0, arg1)
end

function floor_mod(::Type{FastMath}, arg0::jint, arg1::jint)
    return jcall(FastMath, "floorMod", jint, (jint, jint), arg0, arg1)
end

function floor_mod(::Type{FastMath}, arg0::jlong, arg1::jint)
    return jcall(FastMath, "floorMod", jint, (jlong, jint), arg0, arg1)
end

function floor_mod(::Type{FastMath}, arg0::jlong, arg1::jlong)
    return jcall(FastMath, "floorMod", jlong, (jlong, jlong), arg0, arg1)
end

function fma(::Type{FastMath}, arg0::jdouble, arg1::jdouble, arg2::jdouble)
    return jcall(FastMath, "fma", jdouble, (jdouble, jdouble, jdouble), arg0, arg1, arg2)
end

function fma(::Type{FastMath}, arg0::jfloat, arg1::jfloat, arg2::jfloat)
    return jcall(FastMath, "fma", jfloat, (jfloat, jfloat, jfloat), arg0, arg1, arg2)
end

function get_exponent(::Type{FastMath}, arg0::jdouble)
    return jcall(FastMath, "getExponent", jint, (jdouble,), arg0)
end

function get_exponent(::Type{FastMath}, arg0::jfloat)
    return jcall(FastMath, "getExponent", jint, (jfloat,), arg0)
end

function hypot(::Type{FastMath}, arg0::RealFieldElement, arg1::RealFieldElement)
    return jcall(FastMath, "hypot", RealFieldElement, (RealFieldElement, RealFieldElement), arg0, arg1)
end

function hypot(::Type{FastMath}, arg0::jdouble, arg1::jdouble)
    return jcall(FastMath, "hypot", jdouble, (jdouble, jdouble), arg0, arg1)
end

function i_e_e_eremainder(::Type{FastMath}, arg0::RealFieldElement, arg1::RealFieldElement)
    return jcall(FastMath, "IEEEremainder", RealFieldElement, (RealFieldElement, RealFieldElement), arg0, arg1)
end

function i_e_e_eremainder(::Type{FastMath}, arg0::RealFieldElement, arg1::jdouble)
    return jcall(FastMath, "IEEEremainder", RealFieldElement, (RealFieldElement, jdouble), arg0, arg1)
end

function i_e_e_eremainder(::Type{FastMath}, arg0::jdouble, arg1::jdouble)
    return jcall(FastMath, "IEEEremainder", jdouble, (jdouble, jdouble), arg0, arg1)
end

function increment_exact(::Type{FastMath}, arg0::jint)
    return jcall(FastMath, "incrementExact", jint, (jint,), arg0)
end

function increment_exact(::Type{FastMath}, arg0::jlong)
    return jcall(FastMath, "incrementExact", jlong, (jlong,), arg0)
end

function log(::Type{FastMath}, arg0::RealFieldElement)
    return jcall(FastMath, "log", RealFieldElement, (RealFieldElement,), arg0)
end

function log(::Type{FastMath}, arg0::jdouble)
    return jcall(FastMath, "log", jdouble, (jdouble,), arg0)
end

function log(::Type{FastMath}, arg0::jdouble, arg1::jdouble)
    return jcall(FastMath, "log", jdouble, (jdouble, jdouble), arg0, arg1)
end

function log10(::Type{FastMath}, arg0::RealFieldElement)
    return jcall(FastMath, "log10", RealFieldElement, (RealFieldElement,), arg0)
end

function log10(::Type{FastMath}, arg0::jdouble)
    return jcall(FastMath, "log10", jdouble, (jdouble,), arg0)
end

function log1p(::Type{FastMath}, arg0::RealFieldElement)
    return jcall(FastMath, "log1p", RealFieldElement, (RealFieldElement,), arg0)
end

function log1p(::Type{FastMath}, arg0::jdouble)
    return jcall(FastMath, "log1p", jdouble, (jdouble,), arg0)
end

function main(::Type{FastMath}, arg0::Vector{JString})
    return jcall(FastMath, "main", void, (Vector{JString},), arg0)
end

function max(::Type{FastMath}, arg0::RealFieldElement, arg1::RealFieldElement)
    return jcall(FastMath, "max", RealFieldElement, (RealFieldElement, RealFieldElement), arg0, arg1)
end

function max(::Type{FastMath}, arg0::jdouble, arg1::jdouble)
    return jcall(FastMath, "max", jdouble, (jdouble, jdouble), arg0, arg1)
end

function max(::Type{FastMath}, arg0::jfloat, arg1::jfloat)
    return jcall(FastMath, "max", jfloat, (jfloat, jfloat), arg0, arg1)
end

function max(::Type{FastMath}, arg0::jint, arg1::jint)
    return jcall(FastMath, "max", jint, (jint, jint), arg0, arg1)
end

function max(::Type{FastMath}, arg0::jlong, arg1::jlong)
    return jcall(FastMath, "max", jlong, (jlong, jlong), arg0, arg1)
end

function min(::Type{FastMath}, arg0::RealFieldElement, arg1::RealFieldElement)
    return jcall(FastMath, "min", RealFieldElement, (RealFieldElement, RealFieldElement), arg0, arg1)
end

function min(::Type{FastMath}, arg0::jdouble, arg1::jdouble)
    return jcall(FastMath, "min", jdouble, (jdouble, jdouble), arg0, arg1)
end

function min(::Type{FastMath}, arg0::jfloat, arg1::jfloat)
    return jcall(FastMath, "min", jfloat, (jfloat, jfloat), arg0, arg1)
end

function min(::Type{FastMath}, arg0::jint, arg1::jint)
    return jcall(FastMath, "min", jint, (jint, jint), arg0, arg1)
end

function min(::Type{FastMath}, arg0::jlong, arg1::jlong)
    return jcall(FastMath, "min", jlong, (jlong, jlong), arg0, arg1)
end

function multiply_exact(::Type{FastMath}, arg0::jint, arg1::jint)
    return jcall(FastMath, "multiplyExact", jint, (jint, jint), arg0, arg1)
end

function multiply_exact(::Type{FastMath}, arg0::jlong, arg1::jint)
    return jcall(FastMath, "multiplyExact", jlong, (jlong, jint), arg0, arg1)
end

function multiply_exact(::Type{FastMath}, arg0::jlong, arg1::jlong)
    return jcall(FastMath, "multiplyExact", jlong, (jlong, jlong), arg0, arg1)
end

function multiply_full(::Type{FastMath}, arg0::jint, arg1::jint)
    return jcall(FastMath, "multiplyFull", jlong, (jint, jint), arg0, arg1)
end

function multiply_high(::Type{FastMath}, arg0::jlong, arg1::jlong)
    return jcall(FastMath, "multiplyHigh", jlong, (jlong, jlong), arg0, arg1)
end

function next_after(::Type{FastMath}, arg0::jdouble, arg1::jdouble)
    return jcall(FastMath, "nextAfter", jdouble, (jdouble, jdouble), arg0, arg1)
end

function next_after(::Type{FastMath}, arg0::jfloat, arg1::jdouble)
    return jcall(FastMath, "nextAfter", jfloat, (jfloat, jdouble), arg0, arg1)
end

function next_down(::Type{FastMath}, arg0::jdouble)
    return jcall(FastMath, "nextDown", jdouble, (jdouble,), arg0)
end

function next_down(::Type{FastMath}, arg0::jfloat)
    return jcall(FastMath, "nextDown", jfloat, (jfloat,), arg0)
end

function next_up(::Type{FastMath}, arg0::jdouble)
    return jcall(FastMath, "nextUp", jdouble, (jdouble,), arg0)
end

function next_up(::Type{FastMath}, arg0::jfloat)
    return jcall(FastMath, "nextUp", jfloat, (jfloat,), arg0)
end

function pow(::Type{FastMath}, arg0::RealFieldElement, arg1::RealFieldElement)
    return jcall(FastMath, "pow", RealFieldElement, (RealFieldElement, RealFieldElement), arg0, arg1)
end

function pow(::Type{FastMath}, arg0::RealFieldElement, arg1::jdouble)
    return jcall(FastMath, "pow", RealFieldElement, (RealFieldElement, jdouble), arg0, arg1)
end

function pow(::Type{FastMath}, arg0::RealFieldElement, arg1::jint)
    return jcall(FastMath, "pow", RealFieldElement, (RealFieldElement, jint), arg0, arg1)
end

function pow(::Type{FastMath}, arg0::jdouble, arg1::jdouble)
    return jcall(FastMath, "pow", jdouble, (jdouble, jdouble), arg0, arg1)
end

function pow(::Type{FastMath}, arg0::jdouble, arg1::jint)
    return jcall(FastMath, "pow", jdouble, (jdouble, jint), arg0, arg1)
end

function pow(::Type{FastMath}, arg0::jdouble, arg1::jlong)
    return jcall(FastMath, "pow", jdouble, (jdouble, jlong), arg0, arg1)
end

function random(::Type{FastMath})
    return jcall(FastMath, "random", jdouble, ())
end

function rint(::Type{FastMath}, arg0::RealFieldElement)
    return jcall(FastMath, "rint", RealFieldElement, (RealFieldElement,), arg0)
end

function rint(::Type{FastMath}, arg0::jdouble)
    return jcall(FastMath, "rint", jdouble, (jdouble,), arg0)
end

function round(::Type{FastMath}, arg0::RealFieldElement)
    return jcall(FastMath, "round", jlong, (RealFieldElement,), arg0)
end

function round(::Type{FastMath}, arg0::jdouble)
    return jcall(FastMath, "round", jlong, (jdouble,), arg0)
end

function round(::Type{FastMath}, arg0::jfloat)
    return jcall(FastMath, "round", jint, (jfloat,), arg0)
end

function scalb(::Type{FastMath}, arg0::RealFieldElement, arg1::jint)
    return jcall(FastMath, "scalb", RealFieldElement, (RealFieldElement, jint), arg0, arg1)
end

function scalb(::Type{FastMath}, arg0::jdouble, arg1::jint)
    return jcall(FastMath, "scalb", jdouble, (jdouble, jint), arg0, arg1)
end

function scalb(::Type{FastMath}, arg0::jfloat, arg1::jint)
    return jcall(FastMath, "scalb", jfloat, (jfloat, jint), arg0, arg1)
end

function signum(::Type{FastMath}, arg0::RealFieldElement)
    return jcall(FastMath, "signum", RealFieldElement, (RealFieldElement,), arg0)
end

function signum(::Type{FastMath}, arg0::jdouble)
    return jcall(FastMath, "signum", jdouble, (jdouble,), arg0)
end

function signum(::Type{FastMath}, arg0::jfloat)
    return jcall(FastMath, "signum", jfloat, (jfloat,), arg0)
end

function sin(::Type{FastMath}, arg0::RealFieldElement)
    return jcall(FastMath, "sin", RealFieldElement, (RealFieldElement,), arg0)
end

function sin(::Type{FastMath}, arg0::jdouble)
    return jcall(FastMath, "sin", jdouble, (jdouble,), arg0)
end

function sin_cos(::Type{FastMath}, arg0::RealFieldElement)
    return jcall(FastMath, "sinCos", FieldSinCos, (RealFieldElement,), arg0)
end

function sin_cos(::Type{FastMath}, arg0::jdouble)
    return jcall(FastMath, "sinCos", SinCos, (jdouble,), arg0)
end

function sinh(::Type{FastMath}, arg0::RealFieldElement)
    return jcall(FastMath, "sinh", RealFieldElement, (RealFieldElement,), arg0)
end

function sinh(::Type{FastMath}, arg0::jdouble)
    return jcall(FastMath, "sinh", jdouble, (jdouble,), arg0)
end

function sqrt(::Type{FastMath}, arg0::RealFieldElement)
    return jcall(FastMath, "sqrt", RealFieldElement, (RealFieldElement,), arg0)
end

function sqrt(::Type{FastMath}, arg0::jdouble)
    return jcall(FastMath, "sqrt", jdouble, (jdouble,), arg0)
end

function subtract_exact(::Type{FastMath}, arg0::jint, arg1::jint)
    return jcall(FastMath, "subtractExact", jint, (jint, jint), arg0, arg1)
end

function subtract_exact(::Type{FastMath}, arg0::jlong, arg1::jlong)
    return jcall(FastMath, "subtractExact", jlong, (jlong, jlong), arg0, arg1)
end

function tan(::Type{FastMath}, arg0::RealFieldElement)
    return jcall(FastMath, "tan", RealFieldElement, (RealFieldElement,), arg0)
end

function tan(::Type{FastMath}, arg0::jdouble)
    return jcall(FastMath, "tan", jdouble, (jdouble,), arg0)
end

function tanh(::Type{FastMath}, arg0::RealFieldElement)
    return jcall(FastMath, "tanh", RealFieldElement, (RealFieldElement,), arg0)
end

function tanh(::Type{FastMath}, arg0::jdouble)
    return jcall(FastMath, "tanh", jdouble, (jdouble,), arg0)
end

function to_degrees(::Type{FastMath}, arg0::RealFieldElement)
    return jcall(FastMath, "toDegrees", RealFieldElement, (RealFieldElement,), arg0)
end

function to_degrees(::Type{FastMath}, arg0::jdouble)
    return jcall(FastMath, "toDegrees", jdouble, (jdouble,), arg0)
end

function to_int_exact(::Type{FastMath}, arg0::jlong)
    return jcall(FastMath, "toIntExact", jint, (jlong,), arg0)
end

function to_radians(::Type{FastMath}, arg0::RealFieldElement)
    return jcall(FastMath, "toRadians", RealFieldElement, (RealFieldElement,), arg0)
end

function to_radians(::Type{FastMath}, arg0::jdouble)
    return jcall(FastMath, "toRadians", jdouble, (jdouble,), arg0)
end

function ulp(::Type{FastMath}, arg0::jdouble)
    return jcall(FastMath, "ulp", jdouble, (jdouble,), arg0)
end

function ulp(::Type{FastMath}, arg0::jfloat)
    return jcall(FastMath, "ulp", jfloat, (jfloat,), arg0)
end


text:
module DustExtinction

using FITSIO  # for SFD98 dust maps

import Base: show

export ccm89,
       od94,
       download_sfd98,
       SFD98Map,
       ebv_galactic

# -----------------------------------------------------------------------------
# Extinction Laws

# Optical coefficients
const ccm89_ca = [1., 0.17699, -0.50447, -0.02427, 0.72085, 0.01979, -0.77530,
                  0.32999]
const ccm89_cb = [0., 1.41338, 2.28305, 1.07233, -5.38434, -0.62251, 5.30260,
                  -2.09002]
const od94_ca = [1., 0.104, -0.609, 0.701, 1.137, -1.718, -0.827, 1.647,
                 -0.505]
const od94_cb = [0., 1.952, 2.908, -3.989, -7.985, 11.102, 5.491, -10.805,
                 3.347]

function ccm89like(w::Real, r_v, c_a, c_b)
    x = 1.e4 / w
    a = 0.
    b = 0.
    if x < 0.3
        error("wavelength out of range")
    elseif x < 1.1  # Near IR
        y = x^1.61
        a = 0.574y
        b = -0.527y
    elseif x < 3.3  # Optical
        y = x - 1.82
        yn = 1.
        a = c_a[1]
        b = c_b[1]
        for i = 2:length(c_a)
            yn *= y
            a += c_a[i] * yn
            b += c_b[i] * yn
        end
    elseif x < 8.  # UV
        a =  1.752 - 0.316x - (0.104 / ((x-4.67)^2 + 0.341))
        b = -3.090 + 1.825x + (1.206 / ((x-4.62)^2 + 0.263))
        if x > 5.9
            y = x - 5.9
            y2 = y * y
            y3 = y2 * y
            a += -0.04473y2 - 0.009779y3
            b += 0.2130y2 + 0.1207y3
        end
    elseif x < 10.
        y = x - 8.
        y2 = y * y
        y3 = y2 * y
        a = -0.070y3 + 0.137y2 - 0.628y - 1.073
        b = 0.374y3 - 0.420y2 + 4.257y + 13.670
    else
        error("wavelength out of range")
    end

    a + b/r_v
end

"""
`ccm89(wave, r_v)`

Clayton, Cardelli and Mathis (1989) dust law. Returns the extinction
in magnitudes at the given wavelength(s) `wave` (in Angstroms),
relative to the extinction at 5494.5 Angstroms. The parameter `r_v`
changes the shape of the function.  A typical value for the Milky Way
is 3.1. An error is raised for wavelength values outside the range of
support, 1000. to 33333.33 Angstroms.
"""
ccm89(w::Real, r_v) = ccm89like(w, r_v, ccm89_ca, ccm89_cb)

"""
`od94(wave, r_v)`

O'Donnell (1994) dust law, which is identical to the Clayton, Cardelli
and Mathis (1989) dust law, except that different coefficients are
used in the optical (3030.3 to 9090.9 Angstroms). Returns the
extinction in magnitudes at the given wavelength(s) `wave` (in
Angstroms), relative to the extinction at 5494.5 Angstroms. The
parameter `r_v` changes the shape of the function.  A typical value
for the Milky Way is 3.1.  An error is raised for wavelength values
outside the range of support, 1000. to 33333.33 Angstroms.
"""
od94(w::Real, r_v) = ccm89like(w, r_v, od94_ca, od94_cb)

# Vectorized versions (vectorized on wavelength only)
for f = (:ccm89, :od94)
    @eval begin
        ($f)(w::AbstractArray{T,1}, r_v) where T<:Real =
            [ ($f)(w[i], r_v) for i=1:length(w) ]
        ($f)(w::AbstractArray{T,2}, r_v) where T<:Real =
            [ ($f)(w[i,j], r_v) for i=1:size(w,1), j=1:size(w,2) ]
        ($f)(w::AbstractArray{T}, r_v) where T<:Real  =
            reshape([ ($f)(w[i], r_v) for i=1:length(w) ], size(w))
    end
end

# -----------------------------------------------------------------------------
# SFD98 Dust Maps

# It would be good to find a more permanent remote location of these maps,
# but they're here for the time being.
const SFD98_BASEURL = "http://sncosmo.github.io/data/dust/"

# Download 4096x4096 E(B-V) maps to $SFD98_DIR directory.
"""
`download_sfd98([destdir])`

Download the Schlegel, Finkbeiner and Davis (1998) dust maps to the given
directory. If the directory is ommitted, the `SFD98_DIR` environment variable
is used as the destination directory.
"""
function download_sfd98(destdir::AbstractString)
    for fname in ["SFD_dust_4096_ngp.fits", "SFD_dust_4096_sgp.fits"]
        dest = joinpath(destdir, fname)
        if isfile(dest)
            info("$dest already exists, skipping download.")
        else
            download(SFD98_BASEURL * fname, dest)
        end
    end
end

function download_sfd98()
    haskey(ENV, "SFD98_DIR") || error("SFD98_DIR environment variable not set")
    destdir = ENV["SFD98_DIR"]
    download_sfd98(destdir)
end

"""
`SFD98Map([mapdir])`

Schlegel, Finkbeiner and Davis (1998) dust map. `mapdir` should be a
directory containing the two FITS files defining the map,
`SFD_dust_4096_[ngp,sgp].fits`. If `mapdir` is omitted, the
`SFD98_DIR` environment variable is used. Internally, this type keeps
the FITS files defining the map open, speeding up repeated queries
for E(B-V) values.
"""
mutable struct SFD98Map
    mapdir::String
    ngp::ImageHDU
    ngp_size::Tuple{Int, Int}
    ngp_crpix1::Float64
    ngp_crpix2::Float64
    ngp_lam_scal::Float64
    sgp::ImageHDU
    sgp_size::Tuple{Int, Int}
    sgp_crpix1::Float64
    sgp_crpix2::Float64
    sgp_lam_scal::Float64

    function SFD98Map(mapdir::AbstractString)
        ngp = FITS(joinpath(mapdir, "SFD_dust_4096_ngp.fits"))[1]
        ngp_size = size(ngp)
        ngp_crpix1 = read_key(ngp, "CRPIX1")[1]
        ngp_crpix2 = read_key(ngp, "CRPIX2")[1]
        ngp_lam_scal = read_key(ngp, "LAM_SCAL")[1]
        sgp = FITS(joinpath(mapdir, "SFD_dust_4096_sgp.fits"))[1]
        sgp_size = size(sgp)
        sgp_crpix1 = read_key(sgp, "CRPIX1")[1]
        sgp_crpix2 = read_key(sgp, "CRPIX2")[1]
        sgp_lam_scal = read_key(sgp, "LAM_SCAL")[1]
        new(mapdir,
            ngp, ngp_size, ngp_crpix1, ngp_crpix2, ngp_lam_scal,
            sgp, sgp_size, sgp_crpix1, sgp_crpix2, sgp_lam_scal)
    end
end

function SFD98Map()
    haskey(ENV, "SFD98_DIR") || error("SFD98_DIR environment variable not set")
    SFD98Map(ENV["SFD98_DIR"])
end

show(io::IO, map::SFD98Map) = print(io, "SFD98Map(\"$(map.mapdir)\")")

# Convert from galactic longitude/latitude to lambert pixels.
# See SFD 98 Appendix C. For the 4096x4096 maps, lam_scal = 2048,
# crpix1 = 2048.5, crpix2 = 2048.5.
function galactic_to_lambert(crpix1, crpix2, lam_scal, n, l, b)
    x = lam_scal * sqrt(1. - n*sin(b)) * cos(l) + crpix1
    y = -lam_scal * n * sqrt(1. - n*sin(b)) * sin(l) + crpix2
    return x, y
end

# l, b: galactic coordinates in radians
#
# uses bilinear interpolation
"""
`ebv_galactic(dustmap, l, b)`

Get E(B-V) value from a `SFD98Map` instance at galactic coordinates
(`l`, `b`), given in radians. `l` and `b` may be Vectors. Uses linear
interpolation between pixel values.
"""
function ebv_galactic(dustmap::SFD98Map, l::Real, b::Real)
    if b >= 0.
        hdu = dustmap.ngp
        crpix1 = dustmap.ngp_crpix1
        crpix2 = dustmap.ngp_crpix2
        lam_scal = dustmap.ngp_lam_scal
        xsize, ysize = dustmap.ngp_size
        n = 1.
    else
        hdu = dustmap.sgp
        crpix1 = dustmap.sgp_crpix1
        crpix2 = dustmap.sgp_crpix2
        lam_scal = dustmap.sgp_lam_scal
        xsize, ysize = dustmap.sgp_size
        n = -1.
    end

    x, y = galactic_to_lambert(crpix1, crpix2, lam_scal, n, l, b)

    # determine interger pixel locations and weights for bilinear interpolation
    xfloor = floor(x)
    xw = x - xfloor
    x0 = round(Int, xfloor)
    yfloor = floor(y)
    yw = y - yfloor
    y0 = round(Int, yfloor)

    # handle cases near l = [0, pi/2. pi, 3pi/2] where two pixels
    # are out of bounds. This is made simpler because we know from the
    # galactic_to_lambert() transform that only x or y will be near
    # the image bounds, but not both.
    if x0 == 0
        data = read(hdu, 1, y0:y0+1)
        val = (1 - yw) * data[1] + yw * data[2]
    elseif x0 == xsize
        data = read(hdu, xsize, y0:y0+1)
        val = (1 - yw) * data[1] + yw * data[2]
    elseif y0 == 0
        data = read(hdu, x0:x0+1, 1)
        val = (1 - xw) * data[1] + xw * data[2]
    elseif y0 == ysize
        data = read(hdu, x0:x0+1, xsize)
        val = (1 - xw) * data[1] + xw * data[2]
    else
        data = read(hdu, x0:x0+1, y0:y0+1)
        val = ((1 -xw) * (1 -yw) * data[1, 1] +
               xw      * (1 -yw) * data[2, 1] +
               (1 -xw) * yw      * data[1, 2] +
               xw      * yw      * data[2, 2])
    end

    return convert(Float64, val)
end

# array version
function ebv_galactic(dustmap::SFD98Map, l::Vector{T}, b::Vector{T}) where T <: Real
    m = length(l)
    length(b) == m || error("length of l and b must match")
    result = Array(Float64, m)
    for i=1:m
        result[i] = ebv_galactic(dustmap, l[i], b[i])
    end
    return result
end

end # module


text:
include("test_utils.jl")
using KoLesky
using Test
import KoLesky.distance

# mindist = minimum([distance(dofs[k],dofs[l]) for k in 1 : length(dofs) for l in 1 : length(dofs) if k != l])
# diam = maximum([distance(dofs[k],dofs[l]) for k in 1 : length(dofs) for l in 1 : length(dofs)])
# sort_distance = sort([distance(dofs[k],dofs[l]) for k in 1 : length(dofs) for l in 1 : length(dofs) if k < l])
# r = sqrt(minimum(sort_distance[2:end] ./ sort_distance[1 : (end - 1)]))
#
# n_bins = 2 * round(Int, log(r, diam / mindist))

# r_P, r_ℓ = r_maximin(dofs; r=r, n_bins=n_bins)
# P, ℓ = KoLesky.naive_maximin(dofs)

# r_P, r_ℓ = r_maximin(dofs; r=r, n_bins=n_bins)
# P, ℓ = KoLesky.naive_maximin(dofs)

@testset "KoLesky.jl" begin
    @testset "Maximin" begin
        dofs = mat2points(rand(3, 10))
        @test maximin(dofs,)[1] == KoLesky.naive_maximin(dofs)[1]
        @test maximin(dofs,)[2] ≈ KoLesky.naive_maximin(dofs)[2]
    end
    @testset "SuperNodes" begin
        P, ℓ = KoLesky.naive_maximin(dofs); λ = 1.5; ρ = 3.0
        @test KoLesky.is_complete(KoLesky.partition_into_supernodes(dofs, ℓ, P, λ, ρ))
    end
end


text:
using Test
using EllZeroTrendFiltering: add_quadratic!

global const TEST_PRINT_LEVEL = 0

Random.seed!(27182)
#---

# Given a matrix with coefficients for quadratic polynomials this function
# constructs the quadratic polynomials and inserts them into a piecewise
# quadratic object
function verify_piecewise_quadratic(c_mat, show_plot=false)

    poly_list = [QuadraticPolynomial(c_mat[k, :]) for k in 1:size(c_mat,1)]

    if show_plot
        xgrid = linspace(-2,2)
        plot(show=true)
        for p in poly_list
            plot!(xgrid, p.(xgrid))
        end
    end

    # Insert quadratics into piecewise quadfatic
    pwq = create_new_pwq(Float64)
    for p in poly_list
        TEST_PRINT_LEVEL > 0 && println("Inserting: ", p)
        add_quadratic!(pwq, p)

        TEST_PRINT_LEVEL > 0 && println("After Insertion: ")
        TEST_PRINT_LEVEL > 0 && println(pwq)
    end

    # Check that the intersections between the quadratics are in increasing order
    TEST_PRINT_LEVEL > 0 && println("Checking that intersections points are increasing: ")
    for λ in pwq
        @test λ.left_endpoint < get_right_endpoint(λ)
    end

    # Check for continuity of the piecewise quadratic
    TEST_PRINT_LEVEL > 0 && println("Checking continuity: ")
    for λ in pwq
        x = get_right_endpoint(λ)
        if x == Inf
            break
        end

        TEST_PRINT_LEVEL > 0 && println("Continutiy diff: ", λ.p(x) - λ.next.p(x))
        @test λ.p(x) ≈ λ.next.p(x)
    end

    # Check that the piecewise quadratic is smaller than all the quadratics
    # that it was formed by
    x_grid = -5:0.1:5
    y_pwq = pwq(x_grid)
    # Test both versions of evaluating piecewise-quadratic
    @test y_pwq ≈ pwq.(x_grid)  rtol=sqrt(eps())

    # Test get_vals
    x, y, x_all, y_all = get_vals(pwq)
    @test pwq(x) ≈ y            rtol = sqrt(eps())

    for p in poly_list
        TEST_PRINT_LEVEL > 0 && @printf("suboptimality diff: %2.3f\n", maximum(y_pwq - p.(x_grid)))
        #println("   ",  x_grid[ (y_pwq - p.(x_grid)) .> 0])
        @test (y_pwq .<= p.(x_grid)) == trues(length(x_grid))
    end

end

pwq = generate_PiecewiseQuadratic(([2.0, 2, 1], -Inf), ([2.0, -2, 1], 0.0))

@test length(pwq) == 2
add_quadratic!(pwq, QuadraticPolynomial([1, 0, 1.0]))

#Test print string
str = "PiecewiseQuadratic{Float64} with 3 elements:\n[  -∞ , -0.50]\t  :   1.00*x^2 + 2.00*x + 2.00   \n[-0.50, 0.50]\t  :   1.00*x^2 + 0.00*x + 1.00   \n[0.50, ∞  ]\t  :   1.00*x^2 + -2.00*x + 2.00   \n"
str2 = sprint(print, pwq)
@test str == str2

@test length(pwq) == 3

#---
# Tests with (2x^2 + 1) and (x^2 + 2) which have intersections in -1, +1
# and 5/2*x(x-1) + 1 = 2.5x^2 - 2.5x + 1

p1 = QuadraticPolynomial([2.0, 0, 1])
p2 = QuadraticPolynomial([1.0, 0, 2])
p3 = QuadraticPolynomial(2.5, -2.5, 1.0)

pwq1 = create_new_pwq(p1)
add_quadratic!(pwq1, p2)
@test length(pwq1) == 3
add_quadratic!(pwq1, p3)
@test length(pwq1) == 4

pwq2 = create_new_pwq(p2)
add_quadratic!(pwq2, p3)
@test length(pwq2) == 3
add_quadratic!(pwq2, p1)
@test length(pwq2) == 4

@test pwq1[1].p === pwq2[1].p == p1
@test pwq1[2].left_endpoint == pwq2[2].left_endpoint == -1
@test pwq1[2].p === pwq2[2].p == p2
@test pwq1[3].left_endpoint == pwq2[3].left_endpoint == 0
@test pwq1[3].p === pwq2[3].p

#---
c_mat1 = [
    2.53106  1.59097   1.60448
    2.51349  0.681205  1.60553
    2.09364  0.714858  1.08607
    ]
verify_piecewise_quadratic(c_mat1)

#---
c_mat2 = [
    2.01532  1.59269   1.65765
    2.50071  1.56421   1.53899
    2.02767  0.706143  1.129
    ]
verify_piecewise_quadratic(c_mat2)

#---

c_mat3 = [
    2.42298  0.796953  1.02011
    2.45055  1.32235   1.4894
    2.1762   0.855014  1.0647
    ]
verify_piecewise_quadratic(c_mat3)
#---
# Random test case
N = 10
c_mat_rand = [2 .+ rand(N) 2*rand(N) 1 .+ rand(N)]

verify_piecewise_quadratic(c_mat_rand)


text:
# Phlia and I always talked about building up a list of algorithms under their types
abstract type algorithm end # capitalisation?
abstract type gradientBased <: algorithm end
abstract type gradientFree <: algorithm end


text:
macro genstruct(x)
    return :(mutable struct $x end)
end

casslib = "libcassandra"

@genstruct CassFuture
@genstruct CassCluster
@genstruct CassSession
@genstruct CassString
@genstruct CassStatement
@genstruct CassResult
@genstruct CassIterator
@genstruct CassRow
@genstruct CassValue
@genstruct CassPrepared
@genstruct CassBatch

function cql_cluster_set_concurrency(cluster::Ptr{CassCluster}, nthreads::Int64)
    val = ccall(
            (:cass_cluster_set_num_threads_io, casslib),
            UInt16,
            (Ptr{CassCluster}, UInt32),
            cluster, nthreads)
    return val::UInt16
end

function cql_cluster_set_connections_per_host(cluster::Ptr{CassCluster}, siz::Int64)
    val1 = ccall( (:cass_cluster_set_max_connections_per_host, "libcassandra"), UInt16, (Ptr{CassCluster}, UInt32), cluster, siz)
    val2 = ccall( (:cass_cluster_set_core_connections_per_host, "libcassandra"), UInt16, (Ptr{CassCluster}, UInt32), cluster, siz)
    val = val1 | val2
    return val::UInt16
end

function cql_cluster_set_write_bytes_high_water_mark(cluster::Ptr{CassCluster}, siz::Int64)
    val = ccall(
        (:cass_cluster_set_write_bytes_high_water_mark, casslib),
        UInt16,
        (Ptr{CassCluster}, UInt32),
        cluster, siz)
    return val::UInt16
end

function cql_cluster_set_pending_requests_high_water_mark(cluster::Ptr{CassCluster}, siz::Int64)
    val = ccall(
        (:cass_cluster_set_pending_requests_high_water_mark, casslib),
        UInt16,
        (Ptr{CassCluster}, UInt32),
        cluster, siz)
    return val::UInt16
end

function cql_cluster_set_queue_size(cluster::Ptr{CassCluster}, siz::Int64)
    val1 = ccall(
        (:cass_cluster_set_queue_size_io, casslib),
        UInt16,
        (Ptr{CassCluster}, UInt32),
        cluster, siz)
    val2 = ccall(
        (:cass_cluster_set_queue_size_event, casslib),
        UInt16,
        (Ptr{CassCluster}, UInt32),
        cluster, siz)
    val = val1 | val2
    return val::UInt16
end

function cql_future_error_code(future::Ptr{CassFuture})
    val = ccall( (:cass_future_error_code, "libcassandra"),UInt16, (Ptr{CassFuture},),future)
    return val::UInt16
end

function cql_future_error_message(future::Ptr{CassFuture}, strref::Ref{Ptr{UInt8}}, siz::Ref{Csize_t})
    ccall( (:cass_future_error_message, "libcassandra"), Nothing, (Ptr{CassFuture}, Ref{Ptr{UInt8}}, Ref{Csize_t}), future, strref, siz)
end

function cql_cluster_new()
    val =  ccall((:cass_cluster_new, "libcassandra"), Ptr{CassCluster},(),)
    return val::Ptr{CassCluster}
end

function cql_session_new()
    val =  ccall((:cass_session_new, "libcassandra"), Ptr{CassSession},(),)
    return val::Ptr{CassSession}
end

function cql_cluster_set_credentials(cluster::Ptr{CassCluster}, username::String, password::String)
    ccall(
        (:cass_cluster_set_credentials, casslib),
        Nothing,
        (Ptr{CassCluster}, Cstring, Cstring),
        cluster, username, password)
end

function cql_cluster_set_contact_points(cluster::Ptr{CassCluster}, hosts::String)
    ccall((:cass_cluster_set_contact_points, "libcassandra"), Nothing,(Ptr{CassCluster}, Cstring), cluster, hosts)
end

function cql_cluster_set_whitelist_filtering(cluster::Ptr{CassCluster}, hosts::String)
    ccall(
        (:cass_cluster_set_whitelist_filtering, casslib),
        Nothing,
        (Ptr{CassCluster}, Cstring),
        cluster, hosts)
end

function cql_cluster_set_blacklist_filtering(cluster::Ptr{CassCluster}, hosts::String)
    ccall(
        (:cass_cluster_set_blacklist_filtering,casslib),
        Nothing,
        (Ptr{CassCluster}, Cstring),
        cluster, hosts)
end

function cql_session_connect(session::Ptr{CassSession}, cluster::Ptr{CassCluster})
    val = ccall( (:cass_session_connect,"libcassandra"), Ptr{CassFuture}, (Ptr{CassSession}, Ptr{CassCluster}), session, cluster)
    return val::Ptr{CassFuture}
end

function cql_session_free(session::Ptr{CassSession})
    ccall(
        (:cass_session_free, casslib),
        Nothing,
        (Ptr{CassCluster},),
        session)
end

function cql_cluster_free(cluster::Ptr{CassCluster})
    ccall( (:cass_cluster_free, "libcassandra"), Nothing, (Ptr{CassCluster},), cluster)
end

function cql_result_row_count(result::Ptr{CassResult})
    val = ccall( (:cass_result_row_count, "libcassandra"), Int32, (Ptr{CassResult},), result)
    return val::Int32
end

function cql_result_column_count(result::Ptr{CassResult})
    val = ccall( (:cass_result_column_count, "libcassandra"), Int32, (Ptr{CassResult},), result)
    return val::Int32
end

function cql_iterator_next(iterator::Ptr{CassIterator})
    next = ccall( (:cass_iterator_next, "libcassandra"), UInt8, (Ptr{CassIterator},), iterator)
    val = ifelse(next == 0, false, true)
    return val::Bool
end

function cql_future_free(future::Ptr{CassFuture})
    ccall( (:cass_future_free, "libcassandra"), Nothing, (Ptr{CassFuture},), future)
end

function cql_result_column_type(result::Ptr{CassResult}, idx::Int64)
    val = ccall( (:cass_result_column_type, "libcassandra"), UInt16, (Ptr{CassResult}, UInt32), result, idx)
    return val::UInt16
end

function cql_value_get_int8(val::Ptr{CassValue}, out::Ref{Cshort})
    err = ccall( (:cass_value_get_int8, "libcassandra"), Cushort, (Ptr{CassValue}, Ref{Cshort}), val, out)
    return err::UInt16
end

function cql_value_get_int16(val::Ptr{CassValue}, out::Ref{Cshort})
    err = ccall( (:cass_value_get_int16, "libcassandra"), Cushort, (Ptr{CassValue}, Ref{Cshort}), val, out)
    return err::UInt16
end

function cql_value_get_int64(val::Ptr{CassValue}, out::Ref{Clonglong})
    err = ccall( (:cass_value_get_int64, "libcassandra"), Cushort, (Ptr{CassValue}, Ref{Clonglong}),val, out)
    return err::UInt16
end

function cql_value_get_int32(val::Ptr{CassValue}, out::Ref{Cint})
    err = ccall( (:cass_value_get_int32, "libcassandra"), Cushort, (Ptr{CassValue}, Ref{Cint}), val, out)
    return err::UInt16
end

function cql_result_column_name(val::Ptr{CassResult}, pos::Int, out::Ref{Ptr{UInt8}}, siz::Ref{Csize_t})
    err = ccall( (:cass_result_column_name, "libcassandra"), Cushort, (Ptr{CassResult}, Clonglong, Ref{Ptr{UInt8}}, Ref{Csize_t}), val, pos, out, siz)
    return err::UInt16
end

function cql_value_get_uint32(val::Ptr{CassValue}, out::Ref{Cuint})
    err = ccall( (:cass_value_get_uint32, "libcassandra"), Cushort, (Ptr{CassValue}, Ref{Cuint}), val, out)
    return err::UInt16
end

function cql_value_get_bool(val::Ptr{CassValue}, out::Ref{Cint})
    err = ccall( (:cass_value_get_bool, "libcassandra"), Cushort, (Ptr{CassValue}, Ref{Cint}), val, out)
    return err::UInt16
end

function cql_value_get_string(val::Ptr{CassValue}, out::Ref{Ptr{UInt8}}, siz::Ref{Csize_t})
    err = ccall( (:cass_value_get_string, "libcassandra"), Cushort, (Ptr{CassValue}, Ref{Ptr{UInt8}}, Ref{Csize_t}), val, out, siz)
    return err::UInt16
end

function cql_value_get_float(val::Ptr{CassValue}, out::Ref{Cfloat})
    err = ccall( (:cass_value_get_float, "libcassandra"), Cushort, (Ptr{CassValue}, Ref{Cfloat}), val, out)
    return err::UInt16
end

function cql_value_get_double(val::Ptr{CassValue}, out::Ref{Cdouble})
    err = ccall( (:cass_value_get_double, "libcassandra"), Cushort, (Ptr{CassValue}, Ref{Cdouble}), val, out)
    return err::UInt16
end

function cql_statement_free(statement::Ptr{CassStatement})
    ccall( (:cass_statement_free,  "libcassandra"), Nothing,(Ptr{CassStatement},), statement)
end

function cql_result_free(result::Ptr{CassResult})
    ccall( (:cass_result_free,  "libcassandra"), Nothing, (Ptr{CassResult},), result)
end

function cql_iterator_free(iterator::Ptr{CassIterator})
    ccall( (:cass_iterator_free,  "libcassandra"), Nothing, (Ptr{CassIterator},), iterator)
end

function cql_statement_new(query::String, params::Int)
    statement = ccall( (:cass_statement_new, "libcassandra"), Ptr{CassStatement}, (Cstring, Clonglong), query, params)
    return statement::Ptr{CassStatement}
end

function cql_statement_set_paging_size(statement::Ptr{CassStatement}, pgsize::Int)
    ccall( (:cass_statement_set_paging_size, "libcassandra"), Nothing, (Ptr{CassStatement}, Cint), statement, pgsize)
end

function cql_statement_set_request_timeout(statement::Ptr{CassStatement}, timeout::Int)
    err = ccall( (:cass_statement_set_request_timeout,  "libcassandra"), Cushort, (Ptr{CassStatement}, Clonglong), statement, timeout)
    return err::UInt16
end

function cql_set_cluster_timeout(cluster::Ptr{CassCluster}, millisecs::Int)
    ccall((:cass_cluster_set_request_timeout, "libcassandra"), Nothing, (Ptr{CassCluster}, Cint), cluster, millisecs)
end

function cql_set_cluster_heartbeat_interval(cluster::Ptr{CassCluster}, millisecs::Int)
    ccall((:cass_cluster_set_connection_heartbeat_interval, "libcassandra"), Nothing, (Ptr{CassCluster}, Cint), cluster, millisecs)
end

function cql_set_cluster_idle_timeout(cluster::Ptr{CassCluster}, secs::Int)
    ccall((:cass_cluster_set_connection_idle_timeout, "libcassandra"),Nothing, (Ptr{CassCluster}, Cint), cluster, secs )
end

function cql_set_cluster_reconnect_timeout(cluster::Ptr{CassCluster}, millisecs::Int)
    ccall((:cass_cluster_set_reconnect_wait_time, "libcassandra"), Nothing, (Ptr{CassCluster}, Cint), cluster, millisecs)
end

function cql_session_execute(session::Ptr{CassSession}, statement::Ptr{CassStatement})
    future = ccall( (:cass_session_execute, "libcassandra"), Ptr{CassFuture}, (Ptr{CassSession}, Ptr{CassStatement}), session, statement)
    return future::Ptr{CassFuture}
end

function cql_future_get_result(future::Ptr{CassFuture})
    result = ccall( (:cass_future_get_result,  "libcassandra"), Ptr{CassResult}, (Ptr{CassFuture},), future)
    return result::Ptr{CassResult}
end

function cql_iterator_from_result(result::Ptr{CassResult})
    iterator = ccall( (:cass_iterator_from_result,  "libcassandra"), Ptr{CassIterator}, (Ptr{CassResult},), result)
    return iterator::Ptr{CassIterator}
end

function cql_iterator_get_row(iterator::Ptr{CassIterator})
    row = ccall( (:cass_iterator_get_row,  "libcassandra"), Ptr{CassRow}, (Ptr{CassIterator},), iterator)
    return row::Ptr{CassRow}
end

function cql_row_get_column(row::Ptr{CassRow}, pos::Int64)
    val = ccall( (:cass_row_get_column,  "libcassandra"), Ptr{CassValue}, (Ptr{CassRow}, Clonglong), row, pos)
    return val::Ptr{CassValue}
end

function cql_result_has_more_pages(result::Ptr{CassResult})
    hasmore = ccall( (:cass_result_has_more_pages, "libcassandra"), Cint, (Ptr{CassResult},), result)
    out = convert(Bool, hasmore)
    return out::Bool
end

function cql_statement_set_paging_state(statement::Ptr{CassStatement}, result::Ptr{CassResult})
    ccall( (:cass_statement_set_paging_state, "libcassandra"), Nothing, (Ptr{CassStatement}, Ptr{CassResult}), statement, result)
end

function cql_future_wait(future::Ptr{CassFuture})
    ccall( (:cass_future_wait, "libcassandra"), Nothing, (Ptr{CassFuture},), future)
end

function cql_session_prepare(session::Ptr{CassSession}, query::String)
    future = ccall( (:cass_session_prepare,  "libcassandra"), Ptr{CassFuture}, (Ptr{CassSession}, Cstring), session, query)
    return future::Ptr{CassFuture}
end

function cql_batch_new(batch_type::UInt8)
    #=
    CASS_BATCH_TYPE_LOGGED = 0x00
    CASS_BATCH_TYPE_UNLOGGED = 0x01
    CASS_BATCH_TYPE_COUNTER = 0x02
    =#
    batch = ccall(
                (:cass_batch_new, casslib),
                Ptr{CassBatch},
                (Cuchar,),
                batch_type)
    return batch::Ptr{CassBatch}
end

function cql_future_get_prepared(future::Ptr{CassFuture})
    prep = ccall(
            (:cass_future_get_prepared,casslib),
            Ptr{CassPrepared},
            (Ptr{CassFuture},),
            future)
    return prep::Ptr{CassPrepared}
end

function cql_prepared_bind(prep::Ptr{CassPrepared})
    statement = ccall(
                    (:cass_prepared_bind, casslib),
                    Ptr{CassStatement},
                    (Ptr{CassPrepared},),
                    prep)
    return statement::Ptr{CassStatement}
end

function cql_statement_bind_string(statement::Ptr{CassStatement}, pos::Int, data::String)
    ccall(
        (:cass_statement_bind_string, casslib),
        Nothing,
        (Ptr{CassStatement}, Cint, Cstring),
        statement, pos, data)
end

function cql_statement_bind_int8(statement::Ptr{CassStatement}, pos::Int, data::Int8)
    ccall(
        (:cass_statement_bind_int8, casslib),
        Nothing,
        (Ptr{CassStatement}, Cint, Cshort),
        statement, pos, data)
end

function cql_statement_bind_int16(statement::Ptr{CassStatement}, pos::Int, data::Int16)
    ccall(
        (:cass_statement_bind_int16, casslib),
        Nothing,
        (Ptr{CassStatement}, Cint, Cshort),
        statement, pos, data)
end

function cql_statement_bind_int32(statement::Ptr{CassStatement}, pos::Int, data::Int32)
    ccall(
        (:cass_statement_bind_int32, casslib),
        Nothing,
        (Ptr{CassStatement}, Cint, Cint),
        statement, pos, data)
end

function cql_statement_bind_int64(statement::Ptr{CassStatement}, pos::Int, data::Int64)
    ccall(
        (:cass_statement_bind_int64, casslib),
        Nothing,
        (Ptr{CassStatement}, Cint, Clonglong),
        statement, pos, data)
end

function cql_statement_bind_bool(statement::Ptr{CassStatement}, pos::Int, data::Bool)
    ccall(
        (:cass_statement_bind_bool, casslib),
        Nothing,
        (Ptr{CassStatement}, Cint, Cint),
        statement, pos, data)
end

function cql_statement_bind_uint32(statement::Ptr{CassStatement}, pos::Int, data::UInt32)
    ccall(
        (:cass_statement_bind_uint32, casslib),
        Nothing,
        (Ptr{CassStatement}, Cint, Cuint),
        statement, pos, data)
end

function cql_statement_bind_double(statement::Ptr{CassStatement}, pos::Int, data::Float64)
    ccall(
        (:cass_statement_bind_double, casslib),
        Nothing,
        (Ptr{CassStatement}, Cint, Cdouble),
        statement, pos, data)
end

function cql_statement_bind_float(statement::Ptr{CassStatement}, pos::Int, data::Float32)
    ccall(
        (:cass_statement_bind_float, casslib),
        Nothing,
        (Ptr{CassStatement}, Cint, Cfloat),
        statement, pos, data)
end

function cql_statement_bind_null(statement::Ptr{CassStatement}, pos::Int, ::Missing)
    ccall(
        (:cass_statement_bind_null, casslib),
        Nothing,
        (Ptr{CassStatement}, Cint),
        statement, pos)
end

function cql_batch_add_statement(batch::Ptr{CassBatch}, statement::Ptr{CassStatement})
    ccall(
        (:cass_batch_add_statement, casslib),
        Nothing,
        (Ptr{CassBatch}, Ptr{CassStatement}),
        batch, statement)
end

function cql_session_execute_batch(session::Ptr{CassSession}, batch::Ptr{CassBatch})
    future = ccall(
                (:cass_session_execute_batch, casslib),
                Ptr{CassFuture},
                (Ptr{CassSession}, Ptr{CassBatch}),
                session, batch)
    return future::Ptr{CassFuture}
end

function cql_batch_free(batch::Ptr{CassBatch})
    ccall(
        (:cass_batch_free, casslib),
        Nothing,
        (Ptr{CassBatch},),
        batch)
end

function cql_prepared_free(prep::Ptr{CassPrepared})
    ccall(
        (:cass_prepared_free, casslib),
        Nothing,
        (Ptr{CassPrepared},),
        prep)
end


text:
# this test file is called directly by hijack, could be a runtests.jl file
# we test that included submodules have Revise working

using Test
include("submodule.jl")

RUN = []
@testset "submodules_tests" begin
    push!(RUN, 1)
end

max_issues_repo_name: rfourquet/InlineTest

max_forks_repo_name: rfourquet/InlineTest


text:
"""
    TypeA()

Lorem ipsum dolor sit amet.
"""
struct TypeA

end


text:
mutable struct NumFieldEmbNfNS{S, U} <: NumFieldEmb{S}
  field::S             # Number field
  base_field_emb::U    # Embedding of base field
  data::Vector{acb}    # For each L = K[x]/(g_i) component a root p(g_i)
  absolute_index::Int  # Bookkeeping
  isreal::Bool         # Whether the embedding is real
  conjugate::Int       # The absolute index of the conjugate

  function NumFieldEmbNfNS{S, U}(field::S,
                                 base_field_emb::U,
                                 data::Vector{acb},
                                 absolute_index::Int,
                                 isreal::Bool,
                                 conjugate::Int) where {S,  U}
    z = new{S, U}(field, base_field_emb, data,
                  absolute_index, isreal, conjugate)
  end
end

function embedding_type(::Type{NfRelNS{T}}) where {T}
  return NumFieldEmbNfNS{NfRelNS{T}, embedding_type(parent_type(T))}
end

embedding_type(L::NfRelNS{T}) where {T} = embedding_type(NfRelNS{T})

isreal(P::NumFieldEmbNfNS) = P.isreal

isimaginary(P::NumFieldEmbNfNS) = !P.isreal

number_field(f::NumFieldEmbNfNS) = f.field

_absolute_index(P::NumFieldEmbNfNS) = P.absolute_index

function conj(f::NumFieldEmbNfNS)
  return complex_embeddings(number_field(f))[f.conjugate]
end

function Base.show(io::IO, ::MIME"text/plain", f::NumFieldEmbNfNS)
  print(io, "Embedding of\n")
  println(io, number_field(f))
  print(io, "extending the \n", f.base_field_emb, "\n")
  print(io, "with roots ≈ ")
  print(io, "[ ")
  for i in 1:length(f.data)
    _print_acb_neatly(io, f.data[i])
    if i < length(f.data)
      print(io, ", ")
    end
  end
  print(io, "]")
end

function Base.show(io::IO, f::NumFieldEmbNfNS)
  print(io, "Embedding corresponding to (", f.base_field_emb, ") and ≈ ")
  print(io, "[ ")
  for i in 1:length(f.data)
    _print_acb_neatly(io, f.data[i])
    if i < length(f.data)
      print(io, ", ")
    end
  end
  print(io, "]")
end

function complex_embeddings(L::NfRelNS{T}; conjugates::Bool = true) where {T}
  res = get_attribute!(L, :complex_embeddings) do
    return _complex_embeddings(L)
  end::Vector{embedding_type(L)}
  if conjugates
    return res
  else
    r, s = signature(L)
    return res[1:r + s]
  end
end

function _complex_embeddings(L::NfRelNS{T}) where {T}
  r, s = signature(L)
  K = base_field(L)
  S = embedding_type(L)
  data = _conjugates_data_new(L, 32)
  ind = 1
  res = Vector{embedding_type(L)}(undef, r + 2*s)
  for (p, rts) in data
    res[ind] = S(L, p, rts, ind, ind <= r, ind > r ? ind + s : ind)
    if ind > r
      res[ind + s] = S(L, conj(p), conj.(rts), ind + s, false, ind)
    end
    ind += 1
  end
  return res
end

function (g::NumFieldEmbNfNS)(a::NfRelNSElem, prec::Int = 32)
  # This is very slow.
  @req number_field(g) === parent(a) "Parent mismatch"
  f = data(a)
  wprec = prec
  L = parent(a)
  r, s = signature(L)
  # We only store data fo non-conjugated complex_embeddings
  if _absolute_index(g) > r + s
    return conj(conj(g)(a, prec))
  end
  K = base_field(L)
  plcK = complex_embeddings(K)
  pols = Vector{Generic.MPoly{acb}}(undef, length(plcK))
  r, s = signature(L)

  while true
    data = _conjugates_data_new(L, wprec)[_absolute_index(g)]
    p = data[1]
    prec1 = precision(parent(data[2][1]))
    for j = 1:length(data[2])
      prec1 = max(prec1, precision(parent(data[2][j])))
    end
    CC = AcbField(prec1, cached = false)
    CCy, y = PolynomialRing(CC, ngens(L), cached = false)
    fatp = map_coefficients(let wprec = wprec; x -> evaluate(x, data[1], wprec) end, f, parent = CCy)
    pt = data[2]
    for c in fatp.coeffs
      c.parent = CC
    end
    for i in 1:ngens(L)
      pt[i].parent = CC
    end

    o = evaluate(fatp, pt)
    if radiuslttwopower(o, -prec)
      return expand!(o, -prec)
    end
    wprec = 2 * wprec
  end
end

evaluate(a::NfRelNSElem, g::NumFieldEmbNfNS, prec::Int = 32) = g(a, prec)

################################################################################
#
#  Conjugates data
#
################################################################################

function _conjugates_data_new(L::NfRelNS{T}, p::Int) where T
  cd = get_attribute(L, :conjugates_data_new)
  if cd === nothing
    D = Dict{Int, Vector{Tuple{embedding_type(base_field(L)), Vector{acb}}}}()
    res = __conjugates_data_new(L, p)
    D[p] = res
    set_attribute!(L, :conjugates_data_new => D)
    return res
  end
  cd::Dict{Int, Vector{Tuple{embedding_type(base_field(L)), Vector{acb}}}}
  if haskey(cd, p)
    res = cd[p]::Vector{Tuple{embedding_type(base_field(L)), Vector{acb}}}
    return res
  end
  res = __conjugates_data_new(L, p)
  cd[p] = res
  return res::Vector{Tuple{embedding_type(base_field(L)), Vector{acb}}}
end

function __conjugates_data_new(L::NfRelNS{T}, p::Int) where T
  data = [_conjugates_data_new(component(L, j)[1], p) for j = 1:ngens(L)]
  plcs = complex_embeddings(base_field(L), conjugates = false)
  r, s = signature(L)
  res = Vector{Tuple{embedding_type(base_field(L)), Vector{acb}}}(undef, r+s)
  r_cnt = 0
  c_cnt = 0
  for P in plcs
    datas = [x for y in data for x in y if x[1] == P]
    if isreal(P)
      ind_real, ind_complex = enumerate_conj_prim_rel(datas)
      for y in ind_real
        r_cnt += 1
        res[r_cnt] = (P, acb[datas[j][2][y[j]] for j = 1:length(y)])
      end
      for y in ind_complex
        c_cnt += 1
        res[r + c_cnt] = (P, acb[datas[j][2][y[j]] for j = 1:length(y)])
      end
    else
      it = cartesian_product_iterator([1:length(x[2]) for x in datas], inplace = true)
      for y in it
        c_cnt += 1
        res[r + c_cnt] = (P, acb[datas[j][2][y[j]] for j = 1:length(y)])
      end
    end
  end
  return res
end

################################################################################
#
#  Restrict
#
################################################################################

function restrict(f::NumFieldEmbNfNS, K::NumField)
  if K === number_field(f)
    return f
  end
  @req _appears_as_base_field(K, number_field(f)) "Number field is not a base field"
  if K === base_field(number_field(f))
    return f.base_field_emb
  else
    return restrict(f.base_field_emb, K)
  end
end

################################################################################
#
#  Restrict from bigger field
#
################################################################################

function restrict(e::NumFieldEmb, f::NumFieldMor{<: NfRelNS, <: Any, <: Any})
  @req number_field(e) === codomain(f) "Number fields do not match"
  L = domain(f)
  emb = complex_embeddings(L)
  # I need to find the embedding of the base_field of L
  K = base_field(L)
  # This is the natural map K -> L
  KtoL = hom(K, L)
  res_base_field = restrict(e, KtoL * f)
  cn = [ all(overlaps.(ee.data, e.(f.(gens(L))))) && ee.base_field_emb === res_base_field for ee in emb]
  @assert count(cn) == 1
  i = findfirst(cn)
  return emb[i]
end