text:
using CUSPARSE
using CUDArt
using Base.Test

m = 25
n = 35
k = 10
blockdim = 5

##############
# test_csrsv #
##############

for (func,func2,typ) in ((:test_csrsv,:test_csrsv2,:CudaSparseMatrixCSR),
                         (:test_cscsv,:test_cscsv2,:CudaSparseMatrixCSC))
    @eval begin
        function $func(elty)
            A = rand(elty,m,m)
            A = triu(A)
            x = rand(elty,m)
            alpha = rand(elty)
            d_x = CudaArray(x)
            d_A = $typ(sparse(A))
            d_y = CUSPARSE.sv('N','T','U',alpha,d_A,d_x,'O')
            h_y = to_host(d_y)
            y = A\(alpha * x)
            @test_approx_eq(y,h_y)
            x = rand(elty,n)
            d_x = CudaArray(x)
            info = CUSPARSE.sv_analysis('N','T','U',d_A,'O')
            @test_throws(DimensionMismatch, CUSPARSE.sv_solve('N','U',alpha,d_A,d_x,info,'O'))
            A = sparse(rand(elty,m,n))
            d_A = $typ(A)
            @test_throws(DimensionMismatch, CUSPARSE.sv_analysis('T','T','U',d_A,'O'))
            CUSPARSE.cusparseDestroySolveAnalysisInfo(info)
        end

        function $func2(elty)
            A = rand(elty,m,m)
            A = triu(A)
            X = rand(elty,m)
            alpha = rand(elty)
            d_X = CudaArray(X)
            d_A = $typ(sparse(A))
            d_Y = CUSPARSE.sv2('N','U',alpha,d_A,d_X,'O')
            h_Y = to_host(d_Y)
            Y = A\(alpha * X)
            @test_approx_eq(Y,h_Y)
            d_y = UpperTriangular(d_A)\d_X
            h_y = to_host(d_y)
            y = A\X
            @test_approx_eq(y,h_y)
            A = sparse(rand(elty,m,n))
            d_A = $typ(A)
            @test_throws(DimensionMismatch, CUSPARSE.sv2('N','U',alpha,d_A,d_X,'O'))
        end
    end
end

types = [Float32,Float64,Complex64,Complex128]
for elty in types
    tic()
    test_csrsv(elty)
    test_cscsv(elty)
    println("cssv took ", toq(), " for ", elty)
    tic()
    test_csrsv2(elty)
    test_cscsv2(elty)
    println("cssv2 took ", toq(), " for ", elty)
end


text:
function plot_qubittimeevo(time_vec,excited_prob)
    p = scatter(;x=time_vec*1E9,y=excited_prob,mode="line")

    l = Layout(;title="Qubit Excited State Probability",
    xaxis=attr(title="Time (ns)", showgrid=false, zeroline=false),
    yaxis=attr(title="Excited Probability", range=[0.0, 1.0]))

    plot(p,l)
end

function plot_photonnumbers(photons)
    p = bar(;x=collect(0:length(photons)-1),y=photons)

    l = Layout(;title="Photon Number Population",
    xaxis=attr(title="n", showgrid=false, zeroline=false),
    yaxis=attr(title="Amplitude"))

    plot(p,l)
end

function plot_densitymatrix(densitymatrix)
    kets = 0
    for m = size(densitymatrix,1):-1:1
        if sum(abs(densitymatrix[m,:]))>0.0001
            kets = collect(0:m-1)
            break
        end
    end

    densitymatrix = densitymatrix[kets+1,kets+1]

    p = heatmap(;z=abs(densitymatrix),
    colorscale="Viridis")

    l = Layout(;title="Magnitude of Density Matrix Elements |ρ_{mn}|",
    xaxis_title="n",
    xaxis_ticktext=kets,
    xaxis_tickvals=kets,
    xaxis_side="top",
    yaxis_title="m",
    yaxis_ticktext=kets,
    yaxis_tickvals=kets,
    yaxis_autorange="reversed")

    plot(p,l)
end

function plot_wignerfunction(matrix;maxdisp=MAXDISP)
    x = linspace(-maxdisp,maxdisp,size(matrix,1))

    p = heatmap(;x=x,y=x,z=matrix,
    zsmooth="best",
    colorscale="Viridis")

    l = Layout(;title="Wigner Function")

    plot(p,l)
end

function plot_wignerfunction3d(matrix;maxdisp=MAXDISP)
    x = linspace(-maxdisp,maxdisp,size(matrix,1))

    p = surface(;x=x,y=x,z=matrix,
    zsmooth="best",
    colorscale="Viridis")

    l = Layout(;title="Wigner Function")

    plot(p,l)
end

ext: jl


text:
# """
# function that trimmes 3 polyA tail of a given fastq entry and outputs trimmed read with number of polyA indicated in read's name
# Arguments:
#     minimum_not_polyA - minimum length of not polyA strech in a read
#     minimum_polyA_length - minimum length of polyA strech
#     minimum_poly_A_between - minimum length of polyA strech that might occure between non A symbols forming a fragment to be trimmed off (starting from the very 3' end)
# Keywords arguments:
#     debug - turns out additional output
#     window_length - windows size for polyA trimming
#     max_nonA_in_window - tollerated nonA number per window
#
# The algorithms works like:
#     1. Trimmes offs remains of adapters at the 3'
#     2. Search for a window wich exceeds nonA limit and has no other high A content window towards 5' end
#     3. Ajusts the A position in the window
# """
# function trim_polyA_from_fastq_record(fq::FASTQ.Record,
#     minimum_not_polyA::Int64,
#     minimum_polyA_length::Int64,
#     minimum_poly_A_between::Int64;
#     debug=false,
#     window_length=10,
#     max_nonA_in_window=1
#     )::Tuple{FASTQ.Record,Bool}

for (seq,trimmed,len) in [
    ("GAATGTATGGTAGGAATGTATTCTCTTGTAGGAATGTAAATCTGTATTAAAAGGGGGTCCAAGCCAGGCCCCCAGGTCTTCAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAACAAAAAAAA","GAATGTATGGTAGGAATGTATTCTCTTGTAGGAATGTAAATCTGTATTAAAAGGGGGTCCAAGCCAGGCCCCCAGGTCTTC",67),
    ("AAATACATTTTCTGACATTGTAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA","AAATACATTTTCTGACATTGT",129),
    ("ATCCTTGGCTTTCTTCCGGGCCGGTTGAAAAAAAAAAAAAAAAAAAAAAAAAAA","ATCCTTGGCTTTCTTCCGGGCCGGTTG",27),
    ("GAATGTATGGTAGGAATGTATTCTCTTGTAGGAATGTAAATCTGTATTAAAAGGGGGTCCAAGCCAGGCCCCCAGGTCTTCAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAACAAAAAAAA","GAATGTATGGTAGGAATGTATTCTCTTGTAGGAATGTAAATCTGTATTAAAAGGGGGTCCAAGCCAGGCCCCCAGGTCTTC",67),
    ("GAATGTATGGTAGGAATGTATTCTCTTGTAGGAATGTAAATCTGTATTAAAAGGGGGTCCAAGCCAGGCCCCCAGGTCTTCAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAACAAAAAAAAGGTTATTTTCTTAAAT","GAATGTATGGTAGGAATGTATTCTCTTGTAGGAATGTAAATCTGTATTAAAAGGGGGTCCAAGCCAGGCCCCCAGGTCTTC",67),
    ("GAATGTATGGTAGGAATGTATTCTCTTGTAGGAATGTAAATCTGTATTAAAAGGGGGTCCAAGCCAGGCCCCCAGGTCTTCAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAGCAAAAAAAAGGTTATTTTCTTAAAT","GAATGTATGGTAGGAATGTATTCTCTTGTAGGAATGTAAATCTGTATTAAAAGGGGGTCCAAGCCAGGCCCCCAGGTCTTC",67),
    ("GAATGTATGGTAGGAATGTATTCTCTTGTAGGAATGTAAATCTGTATAAAAAAAAAAAAAAAAAAATAAAAGGGGGTCCAAGCCAGGCCCCCAGGTCTTC","GAATGTATGGTAGGAATGTATTCTCTTGTAGGAATGTAAATCTGTAT",19)
    ]
    description="blabla"
    q=fill(20,length(seq))
    q_trimmed=fill(20,length(trimmed))
    fq=FASTQ.Record("testing_seq",description,seq,q)
    fq_trimmed=FASTQ.Record("$len:A:testing_seq",description,trimmed,q_trimmed)
    output_fastq = trim_polyA_from_fastq_record(fq,20,10,4,debug=false)
    output_fastq_record=output_fastq[1]
    #println("KuKu",output_fastq[2]," ",output_fastq[3]," ")
    @test FASTQ.sequence(output_fastq_record)  == FASTQ.sequence(fq_trimmed)
    #test length calculation
    @test FASTQ.identifier(output_fastq_record)  == FASTQ.identifier(fq_trimmed)
end

# AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA
# CTTCTCATCTCCTTAATAAAATTTCATTAAAATCCCCTTGAACTCCCAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA
# AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAGAAAAAAAAAAAAAA
# CACTTTGTCTGTACATACTGGCCTCTGTGATTACATAGATCAGCCATTAAAATAAAACAAGCCTTAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA
# AGGCCCGCAAAGGCAAATAAATCCTTGTTTTGTCTTCACCCATGTAATAAAGGTGTTTATTGTTTTGTTCCCACAAAACAAAAAAAAAAAAAAACAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA
# TAAAACATACAAAATAAACATATTTTAAAACAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA
# AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA
# AAAATAAATACCTTCTAAAAATGTGCAGTAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA
# AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAGAATCGGTAAAGGGCCGGGTGTTACTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTT
# AAAAACTAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA
# TCTGTAGTGCTGTATTATTTTCAATAAATCTGGGACAACAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAATCAATTTTTTATTTTTTTTTTTTTTTTTTTTTTTTTTTTTTT
# AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA
# AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA
# CCTGGACCTGTGACATTCTGGACTATTTCTGTGTTTATTTGTGGCCGAGTGTAACAACCATATAATAAATCACCTCTTCCGCTGTTTTAGCTGAAGAATTAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA
# AGTTCGAGACCTGCCTGGGCAATATAGCGAGACCCCGTTCTCCAGAAAAAGGAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAGCAA
# ACAGTGACAACATTAAACTTAGAAAAGTTTTAAAACTTAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAACAAAAAACAAAACAAAAAAA
# CCAGTCCAAGAATCTGTTTAAAGTTCAGACTTCAAATAGTGGCAAATAAAAAGTGCTATTTGTGAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA
# AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA
# AGTCAAATCCATACCACCACTGAGATCTCATTTATTGCCACAGATGCACAAAATAAATAACCCAAAATCACAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA
# AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAATAAAAAAAAAACATAAATAAAAAATAATAAA
# CAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAACAAAACAAAACTCCCGGGCTGACCCGTAGTGGGCTCGCGCCCGTGAATCGCCACTGCACTCCAGCCTGGGCAACATAGCGAGACCCCGTC
# AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA
# AAACACAACGTTTCTCTAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA
# AAATACATTTTCTGACATTGTAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA
# AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAGAAAAAAAAAAAAAAAAA
# GGTGTCTATAAGAAGCTCACGGGCAAGGATGTTAATTTTGAATTCCCAGAGTTTCAATTGTAAACAAAAATGACTAAATAAAAAGTATATATTCACAGTAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA
# AAACACTCTTCAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAGATCCTGAGCTGGAGGCCCAGGGAGATG
# AAATAAAATAAAAAACAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAACAAAAAAAAAAAAAAAAAAAAACAAAAACACACACAAAAAAA
# AGTGAGCCAAGATCACACCACTGCACTCCAGCACTCCAGGTTGGATGACCGAGCACGACTCCATCTCAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA
# GCTGATCGTCAGATCAAATAAAGTTATAAAATTGCAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAC
# TAATTAATTTACCAAGATAATGTAATTACATTTGGTTTTGTAAGGTATACAGCAGTAATCTCCTATTTTGGTGTCAGTTTTTCAATAAAGTTTTGATTATGGGCAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA
# AAAGCAACTCAAGTTTCCTTTAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA
# TTGCAAAACAAGGAAATAAAGGCTTACTTCTTTAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA
# GAAGAGTAGCACATTAAAGTGATTTTATTGTTAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA
# AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA
# AGGACATGGTGCTGGTCCCTCCCTTAGCCCAAGGGAGAGGCAATAAAGAACACAAAGCTGTAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA
# TAAGCTGTTCTTGCATAGGCTCTTAAGCAGCATGGAAAAATGGTTGATGGAAAATAAACATCAGTTTCTAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA
# TTAGGCCCCAGCAACACGTCATTGTGTAAGGAAAAATAAAAGTGCTGCCGTAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA
# AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA
# AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA
# CAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAATAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA
# TGCTGAGAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA
# AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA
# AACCAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAACACAACAAACAAACAAAAGAAAAATA
# AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA
# TAAATTATCTTTTTCAGTCTCCTCTATAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA
# AATGAAAAACCATGATAATTCTTTGTATATAAAATAAACATTTGAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA
# AAATCCAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAACAAAAAACACCACCCCCCCACCCCCCCCCCCCCCCC
# TGAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA
# GATTGGCCAAGGCCGATGGCATCGTCTCAAAGAACTTTTGACTGGAGAGAATCACAGATGTGGAATATTTGTCATAAATAAATAATGAAAACCTAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA
# TTTTCCATGTGTTTATATTGTAAATATTTTTGATTTCATCAAATTATTTATTCATTAAAAGAAATTTTTGTGAAGCACAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA
# AAGGCGCACATGACTCCAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA
# AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA
# AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAGAAAAAAAAA
# AAATGGAAATAAAGCTTTTTGATGCAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAACAACAAAAACCAAAAACAACA

max_stars_repo_name: laurabiggins/polyAterminus

max_issues_repo_name: laurabiggins/polyAterminus

max_forks_repo_name: laurabiggins/polyAterminus


text:
using OrdinaryDiffEq
using Trixi

###############################################################################
# semidiscretization of the linear advection equation

advection_velocity = (0.2, -0.7, 0.5)
equations = LinearScalarAdvectionEquation3D(advection_velocity)

# Create DG solver with polynomial degree = 3 and (local) Lax-Friedrichs/Rusanov flux as surface flux
solver = DGSEM(polydeg=3, surface_flux=flux_lax_friedrichs)

coordinates_min = (-1.0, -1.0, -1.0) # minimum coordinates (min(x), min(y), min(z))
coordinates_max = ( 1.0,  1.0,  1.0) # maximum coordinates (max(x), max(y), max(z))
trees_per_dimension = (1, 1, 1)

# Note that it is not necessary to use mesh polydeg lower than the solver polydeg
# on a Cartesian mesh.
# See https://doi.org/10.1007/s10915-018-00897-9, Section 6.
mesh = P4estMesh(trees_per_dimension, polydeg=3,
                 coordinates_min=coordinates_min, coordinates_max=coordinates_max)

# Refine bottom left quadrant of each tree to level 4
function refine_fn(p8est, which_tree, quadrant)
  if quadrant.level < 2 ||
      (quadrant.x == 0 && quadrant.y == 0 && quadrant.z == 0 && quadrant.level < 4)
    # return true (refine)
    return Cint(1)
  else
    # return false (don't refine)
    return Cint(0)
  end
end

# Refine recursively until each bottom left quadrant of a tree has level 3
# The mesh will be rebalanced before the simulation starts
refine_fn_c = @cfunction(refine_fn, Cint, (Ptr{Trixi.p8est_t}, Ptr{Trixi.p4est_topidx_t}, Ptr{Trixi.p8est_quadrant_t}))
Trixi.refine_p4est!(mesh.p4est, true, refine_fn_c, C_NULL)
Trixi.balance!(mesh)

# Refine everything again for convergence tests
refine_fn2(p8est, which_tree, quadrant) = Cint(1)
refine_fn_c2 = @cfunction(refine_fn2, Cint, (Ptr{Trixi.p8est_t}, Ptr{Trixi.p4est_topidx_t}, Ptr{Trixi.p8est_quadrant_t}))

# Abuse the variable `initial_refinement_level` for this because this will be increased
# during convergence tests
initial_refinement_level = 0
for i in 1:initial_refinement_level
  Trixi.refine_p4est!(mesh.p4est, false, refine_fn_c2, C_NULL)
end

# A semidiscretization collects data structures and functions for the spatial discretization
semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition_convergence_test, solver)

###############################################################################
# ODE solvers, callbacks etc.

# Create ODE problem with time span from 0.0 to 1.0
tspan = (0.0, 1.0)
ode = semidiscretize(semi, tspan)

# At the beginning of the main loop, the SummaryCallback prints a summary of the simulation setup
# and resets the timers
summary_callback = SummaryCallback()

# The AnalysisCallback allows to analyse the solution in regular intervals and prints the results
analysis_callback = AnalysisCallback(semi, interval=100)

# The SaveSolutionCallback allows to save the solution to a file in regular intervals
save_solution = SaveSolutionCallback(interval=100,
                                     solution_variables=cons2prim)

# The StepsizeCallback handles the re-calculcation of the maximum Δt after each time step
stepsize_callback = StepsizeCallback(cfl=1.6)

# Create a CallbackSet to collect all callbacks such that they can be passed to the ODE solver
callbacks = CallbackSet(summary_callback, analysis_callback, save_solution, stepsize_callback)

###############################################################################
# run the simulation

# OrdinaryDiffEq's `solve` method evolves the solution in time and executes the passed callbacks
sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false),
            dt=1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
            save_everystep=false, callback=callbacks);

# Print the timer summary
summary_callback()


text:
# Cutoff strategies for long-range interactions

export
    NoCutoff,
    DistanceCutoff,
    ShiftedPotentialCutoff,
    ShiftedForceCutoff

"""
    NoCutoff()

Placeholder cutoff that does not alter forces or potentials.
"""
struct NoCutoff end

cutoff_points(::Type{NoCutoff}) = 0

"""
    DistanceCutoff(dist_cutoff)

Cutoff that sets the potential and force to be zero past a specified cutoff point.
"""
struct DistanceCutoff{D, S, I}
    dist_cutoff::D
    sqdist_cutoff::S
    inv_sqdist_cutoff::I
end

function DistanceCutoff(dist_cutoff)
    return DistanceCutoff(dist_cutoff, dist_cutoff ^ 2, inv(dist_cutoff ^ 2))
end

cutoff_points(::Type{DistanceCutoff{D, S, I}}) where {D, S, I} = 1

"""
    ShiftedPotentialCutoff(dist_cutoff)

Cutoff that shifts the potential to be continuous at a specified cutoff point.
"""
struct ShiftedPotentialCutoff{D, S, I}
    dist_cutoff::D
    sqdist_cutoff::S
    inv_sqdist_cutoff::I
end

function ShiftedPotentialCutoff(dist_cutoff)
    return ShiftedPotentialCutoff(dist_cutoff, dist_cutoff ^ 2, inv(dist_cutoff ^ 2))
end

cutoff_points(::Type{ShiftedPotentialCutoff{D, S, I}}) where {D, S, I} = 1

"""
    ShiftedForceCutoff(dist_cutoff)

Cutoff that shifts the force to be continuous at a specified cutoff point.
"""
struct ShiftedForceCutoff{D, S, I}
    dist_cutoff::D
    sqdist_cutoff::S
    inv_sqdist_cutoff::I
end

function ShiftedForceCutoff(dist_cutoff)
    return ShiftedForceCutoff(dist_cutoff, dist_cutoff ^ 2, inv(dist_cutoff ^ 2))
end

cutoff_points(::Type{ShiftedForceCutoff{D, S, I}}) where {D, S, I} = 1

force_divr_cutoff(::NoCutoff, r2, inter, params) = force_divr_nocutoff(inter, r2, inv(r2), params)
potential_cutoff(::NoCutoff, r2, inter, params) = potential(inter, r2, inv(r2), params)

force_divr_cutoff(::DistanceCutoff, r2, inter, params) = force_divr_nocutoff(inter, r2, inv(r2), params)
potential_cutoff(::DistanceCutoff, r2, inter, params) = potential(inter, r2, inv(r2), params)

force_divr_cutoff(::ShiftedPotentialCutoff, r2, inter, params) = force_divr_nocutoff(inter, r2, inv(r2), params)

function potential_cutoff(cutoff::ShiftedPotentialCutoff, r2, inter, params)
    potential(inter, r2, inv(r2), params) - potential(inter, cutoff.sqdist_cutoff, cutoff.inv_sqdist_cutoff, params)
end

function force_divr_cutoff(cutoff::ShiftedForceCutoff, r2, inter, params)
    return force_divr_nocutoff(inter, r2, inv(r2), params) - force_divr_nocutoff(
                                inter, cutoff.sqdist_cutoff, cutoff.inv_sqdist_cutoff, params)
end

@fastmath function potential_cutoff(cutoff::ShiftedForceCutoff, r2, inter, params)
    invr2 = inv(r2)
    r = √r2
    rc = cutoff.dist_cutoff
    fc = force_divr_nocutoff(inter, cutoff.sqdist_cutoff, cutoff.inv_sqdist_cutoff, params) * r

    potential(inter, r2, invr2, params) - (r - rc) * fc -
        potential(inter, cutoff.sqdist_cutoff, cutoff.inv_sqdist_cutoff, params)
end


text:
using TableOperations, Tables, Test

ctable = (A=[1, missing, 3], B=[1.0, 2.0, 3.0], C=["hey", "there", "sailor"])
rtable = Tables.rowtable(ctable)
rtable2 = Iterators.filter(i -> i.a % 2 == 0, [(a=x, b=y) for (x, y) in zip(1:20, 21:40)])

struct ReallyWideTable
end
Tables.istable(::Type{ReallyWideTable}) = true
Tables.columnaccess(::Type{ReallyWideTable}) = true
Tables.columns(x::ReallyWideTable) = x
Tables.columnnames(x::ReallyWideTable) = [Symbol(:x, i) for i = 1:100_000]
Tables.getcolumn(x::ReallyWideTable, i::Int) = rand(10)
Tables.getcolumn(x::ReallyWideTable, nm::Symbol) = rand(10)
Tables.schema(x::ReallyWideTable) = Tables.Schema(Tables.columnnames(x), [Float64 for _ = 1:100_000])

@testset "TableOperations.transform" begin

tran = ctable |> TableOperations.transform(C=Symbol)
@test Tables.istable(typeof(tran))
@test !Tables.rowaccess(typeof(tran))
@test Tables.columnaccess(typeof(tran))
@test Tables.columns(tran) === tran
@test isequal(Tables.getcolumn(tran, :A), [1,missing,3])
@test isequal(Tables.getcolumn(tran, 1), [1,missing,3])

tran2 = rtable |> TableOperations.transform(C=Symbol)
@test Tables.istable(typeof(tran2))
@test Tables.rowaccess(typeof(tran2))
@test !Tables.columnaccess(typeof(tran2))
@test Tables.rows(tran2) === tran2
@test Base.IteratorSize(typeof(tran2)) == Base.HasShape{1}()
@test length(tran2) == 3
@test eltype(tran2) == TableOperations.TransformsRow{NamedTuple{(:A, :B, :C),Tuple{Union{Missing, Int},Float64,String}},NamedTuple{(:C,),Tuple{DataType}}}
trow = first(tran2)
@test trow.A === 1
@test trow.B === 1.0
@test trow.C == :hey
@test Tables.getcolumn(trow, 1) == 1
@test Tables.getcolumn(trow, :A) == 1
ctable2 = Tables.columntable(tran2)
@test isequal(ctable2.A, ctable.A)
@test ctable2.C == map(Symbol, ctable.C)

# test various ways of inputting TableOperations.transform functions
table = TableOperations.transform(ctable, Dict{String, Base.Callable}("C" => Symbol)) |> Tables.columntable
@test table.C == [:hey, :there, :sailor]

table = ctable |> TableOperations.transform(C=Symbol) |> Tables.columntable
@test table.C == [:hey, :there, :sailor]

table = TableOperations.transform(ctable, Dict{Symbol, Base.Callable}(:C => Symbol)) |> Tables.columntable
@test table.C == [:hey, :there, :sailor]

table = TableOperations.transform(ctable, Dict{Int, Base.Callable}(3 => Symbol)) |> Tables.columntable
@test table.C == [:hey, :there, :sailor]

# test simple TableOperations.transforms + return types
table = ctable |> TableOperations.transform(Dict("A"=>x->x+1)) |> Tables.columntable
@test isequal(table.A, [2, missing, 4])
@test typeof(table.A) == Vector{Union{Missing, Int}}

table = ctable |> TableOperations.transform(Dict("A"=>x->coalesce(x+1, 0))) |> Tables.columntable
@test table.A == [2, 0, 4]

table = ctable |> TableOperations.transform(Dict("A"=>x->coalesce(x+1, 0.0))) |> Tables.columntable
@test table.A == [2, 0.0, 4]

table = ctable |> TableOperations.transform(Dict(2=>x->x==2.0 ? missing : x)) |> Tables.columntable
@test isequal(table.B, [1.0, missing, 3.0])
@test typeof(table.B) == Vector{Union{Float64, Missing}}

# test row sinks
# test various ways of inputting TableOperations.transform functions
table = TableOperations.transform(ctable, Dict{String, Base.Callable}("C" => Symbol)) |> Tables.rowtable
@test table[1].C == :hey

table = ctable |> TableOperations.transform(C=Symbol) |> Tables.rowtable
@test table[1].C == :hey

table = TableOperations.transform(ctable, Dict{Symbol, Base.Callable}(:C => Symbol)) |> Tables.rowtable
@test table[1].C == :hey

table = TableOperations.transform(ctable, Dict{Int, Base.Callable}(3 => Symbol)) |> Tables.rowtable
@test table[1].C == :hey

# test simple transforms + return types
table = ctable |> TableOperations.transform(Dict("A"=>x->x+1)) |> Tables.rowtable
@test isequal(map(x->x.A, table), [2, missing, 4])
@test typeof(map(x->x.A, table)) == Vector{Union{Missing, Int}}

table = ctable |> TableOperations.transform(Dict("A"=>x->coalesce(x+1, 0))) |> Tables.rowtable
@test map(x->x.A, table) == [2, 0, 4]

table = ctable |> TableOperations.transform(Dict("A"=>x->coalesce(x+1, 0.0))) |> Tables.rowtable
@test map(x->x.A, table) == [2, 0.0, 4]

table = ctable |> TableOperations.transform(Dict(2=>x->x==2.0 ? missing : x)) |> Tables.rowtable
@test isequal(map(x->x.B, table), [1.0, missing, 3.0])
@test typeof(map(x->x.B, table)) == Vector{Union{Float64, Missing}}

end

@testset "TableOperations.select" begin

# 20
x = ReallyWideTable()
sel = TableOperations.select(x, :x1, :x2)
sch = Tables.schema(sel)
@test sch.names == (:x1, :x2)
@test sch.types == (Float64, Float64)
tt = Tables.columntable(sel)
@test tt.x1 isa Vector{Float64}

sel = TableOperations.select(x, 1, 2)
sch = Tables.schema(sel)
@test sch.names == (:x1, :x2)
@test sch.types == (Float64, Float64)
tt = Tables.columntable(sel)
@test tt.x1 isa Vector{Float64}

# 117
sel = TableOperations.select(ctable)
@test Tables.istable(typeof(sel))
@test Tables.schema(sel) == Tables.Schema((), ())
@test Tables.columnaccess(typeof(sel))
@test Tables.columns(sel) === sel
@test propertynames(sel) == ()
@test isequal(Tables.getcolumn(sel, 1), [1, missing, 3])
@test isequal(Tables.getcolumn(sel, :A), [1, missing, 3])
@test Tables.columntable(sel) == NamedTuple()
@test Tables.rowtable(sel) == NamedTuple{(), Tuple{}}[]

sel = ctable |> TableOperations.select(:A)
@test Tables.istable(typeof(sel))
@test Tables.schema(sel) == Tables.Schema((:A,), (Union{Int, Missing},))
@test Tables.columnaccess(typeof(sel))
@test Tables.columns(sel) === sel
@test propertynames(sel) == (:A,)

sel = ctable |> TableOperations.select(1)
@test Tables.istable(typeof(sel))
@test Tables.schema(sel) == Tables.Schema((:A,), (Union{Int, Missing},))
@test Tables.columnaccess(typeof(sel))
@test Tables.columns(sel) === sel
@test propertynames(sel) == (:A,)

sel = TableOperations.select(rtable)
@test Tables.rowaccess(typeof(sel))
@test Tables.rows(sel) === sel
@test Tables.schema(sel) == Tables.Schema((), ())
@test Base.IteratorSize(typeof(sel)) == Base.HasShape{1}()
@test length(sel) == 3
@test Base.IteratorEltype(typeof(sel)) == Base.HasEltype()
@test eltype(sel) == TableOperations.SelectRow{NamedTuple{(:A, :B, :C),Tuple{Union{Missing, Int},Float64,String}},()}
@test Tables.columntable(sel) == NamedTuple()
@test Tables.rowtable(sel) == [NamedTuple(), NamedTuple(), NamedTuple()]
srow = first(sel)
@test propertynames(srow) == ()

sel = rtable |> TableOperations.select(:A)
@test Tables.rowaccess(typeof(sel))
@test Tables.rows(sel) === sel
@test Tables.schema(sel) == Tables.Schema((:A,), (Union{Int, Missing},))
@test Base.IteratorSize(typeof(sel)) == Base.HasShape{1}()
@test length(sel) == 3
@test Base.IteratorEltype(typeof(sel)) == Base.HasEltype()
@test eltype(sel) == TableOperations.SelectRow{NamedTuple{(:A, :B, :C),Tuple{Union{Missing, Int},Float64,String}},(:A,)}
@test isequal(Tables.columntable(sel), (A = [1, missing, 3],))
@test isequal(Tables.rowtable(sel), [(A=1,), (A=missing,), (A=3,)])
srow = first(sel)
@test propertynames(srow) == (:A,)

# Testing issue where we always select the first column values, but using the correct name.
# NOTE: We don't use rtable here because mixed types produce TypeErrors which hide the
# underlying problem.
rtable2 = [(A = 1.0, B = 2.0), (A = 2.0, B = 4.0), (A = 3.0, B = 6.0)]
sel = rtable2 |> TableOperations.select(:B)
@test Tables.rowaccess(typeof(sel))
@test Tables.rows(sel) === sel
@test Tables.schema(sel) == Tables.Schema((:B,), (Float64,))
@test Base.IteratorSize(typeof(sel)) == Base.HasShape{1}()
@test length(sel) == 3
@test Base.IteratorEltype(typeof(sel)) == Base.HasEltype()
@test eltype(sel) == TableOperations.SelectRow{NamedTuple{(:A, :B,),Tuple{Float64,Float64}},(:B,)}
@test isequal(Tables.columntable(sel), (B = [2.0, 4.0, 6.0],))
@test isequal(Tables.rowtable(sel), [(B=2.0,), (B=4.0,), (B=6.0,)])
@test isequal(Tables.columntable(sel), (B = [2.0, 4.0, 6.0],))
@test isequal(Tables.rowtable(sel), [(B=2.0,), (B=4.0,), (B=6.0,)])
srow = first(sel)
@test propertynames(srow) == (:B,)
@test srow.B == 2.0 # What we expect

sel = rtable |> TableOperations.select(1)
@test Tables.rowaccess(typeof(sel))
@test Tables.rows(sel) === sel
@test Tables.schema(sel) == Tables.Schema((:A,), (Union{Int, Missing},))
@test Base.IteratorSize(typeof(sel)) == Base.HasShape{1}()
@test length(sel) == 3
@test Base.IteratorEltype(typeof(sel)) == Base.HasEltype()
@test eltype(sel) == TableOperations.SelectRow{NamedTuple{(:A, :B, :C),Tuple{Union{Missing, Int},Float64,String}},(1,)}
@test isequal(Tables.columntable(sel), (A = [1, missing, 3],))
@test isequal(Tables.rowtable(sel), [(A=1,), (A=missing,), (A=3,)])
srow = first(sel)
@test propertynames(srow) == (:A,)
@test Tables.getcolumn(srow, 1) == 1
@test Tables.getcolumn(srow, :A) == 1

table = ctable |> TableOperations.select(:A) |> Tables.columntable
@test length(table) == 1
@test isequal(table.A, [1, missing, 3])

table = ctable |> TableOperations.select(1) |> Tables.columntable
@test length(table) == 1
@test isequal(table.A, [1, missing, 3])

table = ctable |> TableOperations.select("A") |> Tables.columntable
@test length(table) == 1
@test isequal(table.A, [1, missing, 3])

# column re-ordering
table = ctable |> TableOperations.select(:A, :C) |> Tables.columntable
@test length(table) == 2
@test isequal(table.A, [1, missing, 3])
@test isequal(table[2], ["hey", "there", "sailor"])

table = ctable |> TableOperations.select(1, 3) |> Tables.columntable
@test length(table) == 2
@test isequal(table.A, [1, missing, 3])
@test isequal(table[2], ["hey", "there", "sailor"])

table = ctable |> TableOperations.select(:C, :A) |> Tables.columntable
@test isequal(ctable.A, table.A)
@test isequal(ctable[1], table[2])

table = ctable |> TableOperations.select(3, 1) |> Tables.columntable
@test isequal(ctable.A, table.A)
@test isequal(ctable[1], table[2])

# row sink
table = ctable |> TableOperations.select(:A) |> Tables.rowtable
@test length(table[1]) == 1
@test isequal(map(x->x.A, table), [1, missing, 3])

table = ctable |> TableOperations.select(1) |> Tables.rowtable
@test length(table[1]) == 1
@test isequal(map(x->x.A, table), [1, missing, 3])

table = ctable |> TableOperations.select("A") |> Tables.rowtable
@test length(table[1]) == 1
@test isequal(map(x->x.A, table), [1, missing, 3])

# column re-ordering
table = ctable |> TableOperations.select(:A, :C) |> Tables.rowtable
@test length(table[1]) == 2
@test isequal(map(x->x.A, table), [1, missing, 3])
@test isequal(map(x->x[2], table), ["hey", "there", "sailor"])

table = ctable |> TableOperations.select(1, 3) |> Tables.rowtable
@test length(table[1]) == 2
@test isequal(map(x->x.A, table), [1, missing, 3])
@test isequal(map(x->x[2], table), ["hey", "there", "sailor"])

table = ctable |> TableOperations.select(:C, :A) |> Tables.rowtable
@test isequal(ctable.A, map(x->x.A, table))
@test isequal(ctable[1], map(x->x[2], table))

table = ctable |> TableOperations.select(3, 1) |> Tables.rowtable
@test isequal(ctable.A, map(x->x.A, table))
@test isequal(ctable[1], map(x->x[2], table))

end

@testset "TableOperations.filter" begin

f = TableOperations.filter(x->x.B == 2.0, ctable)
@test Tables.istable(f)
@test Tables.rowaccess(f)
@test Tables.rows(f) === f
@test Tables.schema(f) == Tables.schema(f)
@test Base.IteratorSize(typeof(f)) == Base.SizeUnknown()
@test Base.IteratorEltype(typeof(f)) == Base.HasEltype()
@test eltype(f) == eltype(Tables.rows(ctable))
@test isequal(Tables.columntable(f), Tables.columntable(ctable |> TableOperations.filter(x->x.B == 2.0)))
@test length((TableOperations.filter(x->x.B == 2.0, ctable) |> Tables.columntable).B) == 1
@test length((TableOperations.filter(x->x.B == 2.0, rtable) |> Tables.columntable).B) == 1
@test length(TableOperations.filter(x->x.B == 2.0, ctable) |> Tables.rowtable) == 1
@test length(TableOperations.filter(x->x.B == 2.0, rtable) |> Tables.rowtable) == 1

end

@testset "TableOperations.map" begin

m = TableOperations.map(x->(A=x.A, C=x.C, B=x.B * 2), ctable)
@test Tables.istable(m)
@test Tables.rowaccess(m)
@test Tables.rows(m) === m
@test Tables.schema(m) === nothing
@test Base.IteratorSize(typeof(m)) == Base.HasLength()
@test Base.IteratorEltype(typeof(m)) == Base.EltypeUnknown()
@test isequal(Tables.columntable(m), Tables.columntable(ctable |> TableOperations.map(x->(A=x.A, C=x.C, B=x.B * 2))))
@test (TableOperations.map(x->(A=x.A, C=x.C, B=x.B * 2), ctable) |> Tables.columntable).B == [2.0, 4.0, 6.0]
@test (TableOperations.map(x->(A=x.A, C=x.C, B=x.B * 2), rtable) |> Tables.columntable).B == [2.0, 4.0, 6.0]
@test length(TableOperations.map(x->(A=x.A, C=x.C, B=x.B * 2), ctable) |> Tables.rowtable) == 3
@test length(TableOperations.map(x->(A=x.A, C=x.C, B=x.B * 2), rtable) |> Tables.rowtable) == 3

end

@testset "TableOperations.joinpartitions" begin

p = Tables.partitioner((ctable, ctable))
j = TableOperations.joinpartitions(p)
@test Tables.istable(j)
@test Tables.columnaccess(j)
@test Tables.schema(j) === Tables.schema(ctable)
@test Tables.columnnames(j) == Tables.columnnames(ctable)
@test isequal(Tables.getcolumn(j, 1), vcat(Tables.getcolumn(ctable, 1), Tables.getcolumn(ctable, 1)))
@test isequal(Tables.getcolumn(j, :A), vcat(Tables.getcolumn(ctable, :A), Tables.getcolumn(ctable, :A)))

# Test joinpartitions with promotion
t1 = (A=[1, 2, 3], B=[1, 2, 3], C=["hey", "there", "sailor"])
t2 = (A=[1, missing, 3], B=[1.0, 2.0, 3.0], C=["trim", "the", "sail"])
p = Tables.partitioner((t1, t2))
# Throws a method error trying to convert `missing` to `Int64`
@test_throws MethodError TableOperations.joinpartitions(p)
j = TableOperations.joinpartitions(p; promote=true)
@test Tables.istable(j)
@test Tables.columnaccess(j)
@test Tables.schema(j) !== Tables.schema(t1)
@test Tables.schema(j) === Tables.schema(t2)
@test Tables.columnnames(j) == Tables.columnnames(t1) == Tables.columnnames(t2)
@test isequal(Tables.getcolumn(j, 1), vcat(Tables.getcolumn(t1, 1), Tables.getcolumn(t2, 1)))
@test isequal(Tables.getcolumn(j, :A), vcat(Tables.getcolumn(t1, :A), Tables.getcolumn(t2, :A)))
end

@testset "TableOperations.makepartitions" begin

# columns
@test_throws ArgumentError TableOperations.makepartitions(ctable, 0)
parts = collect.(collect(Tables.partitions(TableOperations.makepartitions(ctable, 2))))
@test length(parts[1]) == 2
@test length(parts[2]) == 1
@test parts[2][1].A == 3

# rows
parts = collect(Tables.partitions(TableOperations.makepartitions(rtable, 2)))
@test length(parts[1]) == 2
@test length(parts[2]) == 1
@test parts[2][1].A == 3

# forward-only row iterator
parts = collect(Tables.partitions(TableOperations.makepartitions(rtable2, 3)))
@test length(parts) == 4
@test length(parts[1]) == 3
@test length(parts[end]) == 1
@test parts[end][1].a == 20

end

@testset "TableOperations.narrowtypes" begin

ctable_type_any = (A=Any[1, missing, 3], B=Any[1.0, 2.0, 3.0], C=Any["hey", "there", "sailor"])

nt = TableOperations.narrowtypes(ctable_type_any)
@test Tables.istable(nt)
@test Tables.columnaccess(nt)
@test Tables.schema(nt) == Tables.schema(ctable)
@test Tables.columnnames(nt) == Tables.columnnames(ctable)

end

@testset "TableOperations.dropmissing" begin

table = ctable |> TableOperations.dropmissing() |> Tables.columntable
@test isequal(table, Tables.columntable(TableOperations.dropmissing(ctable)))
@test length(table |> Tables.columntable) == 3
@test length(table |> Tables.rowtable) == 2

end


text:
# ProximalOperators.jl - library of commonly used functions in optimization, and associated proximal mappings and gradients

module ProximalOperators

using LinearAlgebra

const RealOrComplex{R <: Real} = Union{R, Complex{R}}
const HermOrSym{T, S} = Union{Hermitian{T, S}, Symmetric{T, S}}
const RealBasedArray{R} = AbstractArray{C, N} where {C <: RealOrComplex{R}, N}
const TupleOfArrays{R} = Tuple{RealBasedArray{R}, Vararg{RealBasedArray{R}}}
const ArrayOrTuple{R} = Union{RealBasedArray{R}, TupleOfArrays{R}}
const TransposeOrAdjoint{M} = Union{Transpose{C,M} where C, Adjoint{C,M} where C}
const Maybe{T} = Union{T, Nothing}

export ProximableFunction
export prox, prox!, gradient, gradient!

abstract type ProximableFunction end

# Utilities

include("utilities/approx_inequality.jl")
include("utilities/tuples.jl")
include("utilities/linops.jl")
include("utilities/symmetricpacked.jl")
include("utilities/uniformarrays.jl")
include("utilities/normdiff.jl")

# Basic functions

include("functions/cubeNormL2.jl")
include("functions/elasticNet.jl")
include("functions/huberLoss.jl")
include("functions/indAffine.jl")
include("functions/indBallL0.jl")
include("functions/indBallL1.jl")
include("functions/indBallL2.jl")
include("functions/indBallRank.jl")
include("functions/indBinary.jl")
include("functions/indBox.jl")
include("functions/indFree.jl")
include("functions/indGraph.jl")
include("functions/indHalfspace.jl")
include("functions/indHyperslab.jl")
include("functions/indNonnegative.jl")
include("functions/indNonpositive.jl")
include("functions/indPoint.jl")
include("functions/indPolyhedral.jl")
include("functions/indPSD.jl")
include("functions/indSimplex.jl")
include("functions/indSOC.jl")
include("functions/indSphereL2.jl")
include("functions/indStiefel.jl")
include("functions/indZero.jl")
include("functions/leastSquares.jl")
include("functions/linear.jl")
include("functions/logBarrier.jl")
include("functions/logisticLoss.jl")
include("functions/normL0.jl")
include("functions/normL1.jl")
include("functions/normL2.jl")
include("functions/normL21.jl")
include("functions/nuclearNorm.jl")
include("functions/quadratic.jl")
include("functions/sqrNormL2.jl")
include("functions/sumPositive.jl")
include("functions/sqrHingeLoss.jl")
include("functions/crossEntropy.jl")
include("functions/TotalVariation1D.jl")

# Calculus rules

include("calculus/conjugate.jl")
include("calculus/epicompose.jl")
include("calculus/distL2.jl")
include("calculus/moreauEnvelope.jl")
include("calculus/postcompose.jl")
include("calculus/precompose.jl")
include("calculus/precomposeDiagonal.jl")
include("calculus/regularize.jl")
include("calculus/separableSum.jl")
include("calculus/slicedSeparableSum.jl")
include("calculus/sqrDistL2.jl")
include("calculus/tilt.jl")
include("calculus/translate.jl")
include("calculus/sum.jl")
include("calculus/pointwiseMinimum.jl")

# Functions obtained from basic (as special cases or using calculus rules)

include("functions/hingeLoss.jl")
include("functions/indExp.jl")
include("functions/maximum.jl")
include("functions/normLinf.jl")
include("functions/sumLargest.jl")

function Base.show(io::IO, f::ProximableFunction)
    println(io, "description : ", fun_name(f))
    println(io, "domain      : ", fun_dom(f))
    println(io, "expression  : ", fun_expr(f))
    print(  io, "parameters  : ", fun_params(f))
end

fun_name(  f) = "n/a"
fun_dom(   f) = "n/a"
fun_expr(  f) = "n/a"
fun_params(f) = "n/a"

is_prox_accurate(f::ProximableFunction) = true
is_separable(f::ProximableFunction) = false
is_convex(f::ProximableFunction) = false
is_concave(f::ProximableFunction) = false
is_singleton(f::ProximableFunction) = false
is_cone(f::ProximableFunction) = false
is_affine(f::ProximableFunction) = is_singleton(f)
is_set(f::ProximableFunction) = is_cone(f) || is_affine(f)
is_positively_homogeneous(f::ProximableFunction) = is_cone(f)
is_support(f::ProximableFunction) = is_convex(f) && is_positively_homogeneous(f)
is_smooth(f::ProximableFunction) = false
is_quadratic(f::ProximableFunction) = false
is_generalized_quadratic(f::ProximableFunction) = is_quadratic(f) || is_affine(f)
is_strongly_convex(f::ProximableFunction) = false

"""
**Proximal mapping**

    y, fy = prox(f, x, γ=1.0)

Computes
```math
y = \\mathrm{prox}_{\\gamma f}(x) = \\arg\\min_z \\left\\{ f(z) + \\tfrac{1}{2\\gamma}\\|z-x\\|^2 \\right\\}.
```
Return values:
* `y`: the proximal point ``y``
* `fy`: the value ``f(y)``
"""
function prox(f::ProximableFunction, x::ArrayOrTuple{R}, gamma=R(1)) where R
    y = similar(x)
    fy = prox!(y, f, x, gamma)
    return y, fy
end

"""
**Proximal mapping (in-place)**

    fy = prox!(y, f, x, γ=1.0)

Computes
```math
y = \\mathrm{prox}_{\\gamma f}(x) = \\arg\\min_z \\left\\{ f(z) + \\tfrac{1}{2\\gamma}\\|z-x\\|^2 \\right\\}.
```
The resulting point ``y`` is written to the (pre-allocated) array `y`, which must have the same shape/size as `x`.

Return values:
* `fy`: the value ``f(y)``
"""
prox!
# TODO: should we put the following here instead? And remove the default value for gamma in each subtype
# prox!(y::ArrayOrTuple{R}, f::ProximableFunction, x::ArrayOrTuple{R}) = prox!(y, f, x, R(1))

"""
**Gradient mapping**

    gradfx, fx = gradient(f, x)

Computes the gradient (and value) of ``f`` at ``x``. If ``f`` is only *subdifferentiable* at ``x``, then return a subgradient instead.

Return values:
* `gradfx`: the (sub)gradient of ``f`` at ``x``
* `fx`: the value ``f(x)``
"""
function gradient(f::ProximableFunction, x)
    y = similar(x)
    fx = gradient!(y, f, x)
    return y, fx
end

"""
**Gradient mapping (in-place)**

    gradient!(gradfx, f, x)

Writes ``\\nabla f(x)`` to `gradfx`, which must be pre-allocated and have the same shape/size as `x`. If ``f`` is only *subdifferentiable* at ``x``, then writes a subgradient instead.

Return values:
* `fx`: the value ``f(x)``
"""
gradient!

end


text:
#
# colormaps.jl --
#
# Implements management of colors and colormaps for using with the PGPlot
# library.
#

module Colormaps

export
    RGBVec,
    palette

using Colors
using PGPlot.Bindings
import PGPlot.Bindings: pgqcr, pgscr

const DATA_DIR = normpath(joinpath(@__DIR__, "..", "data"))

"""

`RGBVec{T}(r,g,b)` represents an RGB color whose components have type `T`.

When `T` is a floating-point, the color triplet can be manipulated as a
3-element *vector*, that is `α*v + β*v` for any reals `α` and `β` and
colors `u` and `v` of type `RGBVec{<:AbstractFloat}` yields a color
of type `RGBVec{<:AbstractFloat}`.  This is useful to interpolate colors
and build colormaps.

Having `T = UInt8` can be used to parse/convert colors.

"""
struct RGBVec{T}
    r::T
    g::T
    b::T
end

Base.eltype(::RGBVec{T}) where {T} = T
RGBVec{T}(col::RGBVec{T}) where {T} = col

function RGBVec{T}(col::RGBVec{UInt8}) where {T<:AbstractFloat}
    a = one(T)/T(255)
    return RGBVec{T}(a*col.r, a*col.g, a*col.b)
end

function RGBVec{UInt8}(col::RGBVec{T}) where {T<:AbstractFloat}
    a = T(255)
    return RGBVec{T}(round(UInt8, a*clamp(col.r, zero(T), one(T))),
                     round(UInt8, a*clamp(col.g, zero(T), one(T))),
                     round(UInt8, a*clamp(col.b, zero(T), one(T))))
end

RGBVec{T}(rgb::NTuple{3,Real}) where {T} = RGB{T}(rgb...)
RGBVec(rgb::NTuple{3,Real}) = RGB(rgb...)

RGBVec{T}(col::Colorant) where {T} = RGBVec{T}(RGB(col))
RGBVec(col::Colorant) = RGBVec(RGB(col))

RGBVec{T}(col::RGB) where {T} = RGBVec{T}(col.r, col.g, col.b)
RGBVec(col::RGB) = RGBVec(col.r, col.g, col.b)

RGBVec{T}(col::Symbol) where {T} = RGBVec{T}(String(col))
RGBVec(col::Symbol) = RGBVec(String(col))

RGBVec{T}(col::AbstractString) where {T} = RGBVec{T}(parse(RGB, col))
RGBVec(col::AbstractString) = RGBVec(parse(RGB, col))

Base.convert(::Type{T}, arg::T) where {T<:RGBVec} = arg
Base.convert(::Type{T}, arg::Colorant) where {T<:RGBVec} = T(arg)
Base.convert(::Type{T}, arg::NTuple{3,Real}) where {T<:RGBVec} = T(arg)
Base.convert(::Type{T}, arg::RGBVec) where {T<:RGB} = T(arg)

black(::Type{RGBVec{T}}) where {T} = RGBVec{T}(zero(T), zero(T), zero(T))
white(::Type{RGBVec{T}}) where {T<:AbstractFloat} =
    RGBVec{T}(one(T), one(T), one(T))
white(::Type{RGBVec{T}}) where {T<:Unsigned} =
    RGBVec{T}(typemax(T), typemax(T), typemax(T))
gray(::Type{RGBVec{T}}, f::T) where {T<:AbstractFloat} = RGBVec{T}(f, f, f)
gray(::Type{RGBVec{T}}, f::Real) where {T<:AbstractFloat} =
    gray(RGBVec{T}, T(f))
gray(::Type{RGBVec{T}}, f::Real) where {T<:Unsigned} =
    (g = round(T, typemax(T)*clamp(f, oftype(f, 0), oftype(f, 1)));
     RGBVec{T}(g, g, g))

background(::Type{T}) where {T<:RGBVec} = pgqcr(T, 0)
foreground(::Type{T}) where {T<:RGBVec} = pgqcr(T, 1)

# Extend arithmetic operators to allow for simple color operations such as
# linear interpolation.  This is restricted to floating-point colorants.
Base.:( + )(c1::RGBVec{T1}, c2::RGBVec{T2}) where {T1<:Real,T2<:Real} =
    (T = float(promote_type(T1, T2)); RGBVec{T}(c1) + RGBVec{T}(c2))
Base.:( + )(c1::RGBVec{<:AbstractFloat}, c2::RGBVec{<:AbstractFloat}) =
    RGBVec(c1.r + c2.r, c1.g + c2.g, c1.b + c2.b)
Base.:( - )(c1::RGBVec{T1}, c2::RGBVec{T2}) where {T1<:Real,T2<:Real} =
    (T = float(promote_type(T1, T2)); RGBVec{T}(c1) - RGBVec{T}(c2))
Base.:( - )(c1::RGBVec{<:AbstractFloat}, c2::RGBVec{<:AbstractFloat}) =
    RGBVec(c1.r - c2.r, c1.g - c2.g, c1.b - c2.b)
Base.:( * )(col::RGBVec, α::Real) = α*col
Base.:( * )(α::Real, col::RGBVec{T}) where {T<:Real} = α*RGB{float(T)}(col)
Base.:( * )(α::Real, col::RGBVec{<:AbstractFloat}) =
    RGBVec(α*col.r, α*col.g, α*col.b)
Base.:( / )(col::RGBVec{T}, α::Real) where {T<:Real} = RGB{float(T)}(col)/α
Base.:( / )(col::RGBVec{<:AbstractFloat}, α::Real) =
    RGBVec(col.r/α, col.g/α, col.b/α)
Base.:( \ )(α::Real, col::RGBVec) = col/α

Base.Tuple(col::RGBVec) = (col.r, col.g, col.b)
Base.clamp(col::RGBVec{T}) where {T<:AbstractFloat} =
    RGBVec{T}(clamp(col.r, zero(T), one(T)),
              clamp(col.g, zero(T), one(T)),
              clamp(col.b, zero(T), one(T)))

function Base.tryparse(::Type{RGBVec{T}}, str::AbstractString) where {T}
    cols = split(str, ' ', keepempty=false)
    length(cols) == 3 || return nothing
    r = tryparse(T, cols[1])
    r === nothing && return nothing
    g = tryparse(T, cols[2])
    g === nothing && return nothing
    b = tryparse(T, cols[3])
    b === nothing && return nothing
    return RGBVec(r,g,b)
end

# Query indexed color as an RGBVec structure.
pgqcr(::Type{RGBVec}, ci::Integer) where {T<:AbstractFloat} =
    RGBVec(pgqcr(ci)...)
pgqcr(::Type{RGBVec{T}}, ci::Integer) where {T<:AbstractFloat} =
    RGBVec{T}(pgqcr(ci)...)
pgqcr(::Type{RGBVec{UInt8}}, ci::Integer) where {T<:AbstractFloat} =
    RGBVec{UInt8}(pgqcr(RGBVec, ci))

# Set indexed color with an RGBVec structure.
pgscr(ci::Integer, col::RGBVec{UInt8}) = pgscr(PGInt(ci), RGBVec{PGFloat}(col))
pgscr(ci::Integer, col::RGBVec{<:AbstractFloat}) =
    pgscr(PGInt(ci), col.r, col.g, col.b)

"""

```julia
find_file(name) -> path
```

yields the path to a readable graphics file.

"""
find_file(name::AbstractString) = find_file(convert(String, name))
function find_file(name::String)
    if isfile(name)
        return name
    else
        path = joinpath(DATA_DIR, name)
        if isfile(path)
            return path
        else
            throw_file_not_found(name)
        end
    end
end

@noinline throw_file_not_found(name::AbstractString) =
    throw(ArgumentError(string("file \"", name, "\" not found")))

"""

```julia
load_gist(name) -> lut
```

yields the colormap read in Gist file `name`.

"""
function load_gist(name::AbstractString)
    path = find_file(name)
    lut = Vector{RGBVec{UInt8}}(undef, 0)
    open(path, "r") do io
        load_gist!(lut, io)
    end
    return lut
end

load_gist(io::IO) = load_gist!(Vector{RGBVec{UInt8}}(undef, 0), io::IO)

"""

```julia
load_gist!(lut, name) -> lut
```

overwrites the contents of the colormap `lut` with the contents read in Gist
file `name`.

"""
function load_gist!(lut::AbstractVector{T}, io::IO) where {T<:RGBVec}
    resize!(lut, 0)
    while !eof(io)
        line = readline(io)
        rgb = tryparse(T, line)
        if rgb !== nothing
            push!(lut, rgb)
        end
    end
    return lut
end

"""

```julia
palette(cmap)
```

installs the colormap `cmap` (a name or a look-up table) in the current
plotting device.

The color index range may be specified:

```julia
palette(cmap, cmin, cmax)
```

If `cmin:cmax` is larger than the current index range, an attempt is made to
enlarge it.

Also see [`set_color_ramp`](@ref).

"""
function palette(ident::Union{AbstractString,AbstractVector{RGBVec{UInt8}}},
                 cmin::Union{Nothing,Integer} = nothing,
                 cmax::Union{Nothing,Integer} = nothing)
    palette(ident, get_color_index_range(cmin, cmax)...)
end

function palette(name::AbstractString, cmin::Int, cmax::Int)
    if endswith(name, ".gp")
        lut = load_gist(name)
        palette(lut, cmin, cmax)
    elseif name == "gray" || name == "+gray"
        set_color_ramp(cmin, cmax, 0)
    elseif name == "-gray"
        set_color_ramp(cmin, cmax, 1)
    elseif name == "bg-fg"
        set_color_ramp(cmin, cmax, 2)
    elseif name == "fg-bg"
        set_color_ramp(cmin, cmax, 3)
    elseif name == "red" || name == "+red"
        set_color_ramp(cmin, cmax, black(RGBVec{PGFloat}),
                       RGBVec{PGFloat}(1,0,0))
    elseif name == "-red"
        set_color_ramp(cmin, cmax, RGBVec{PGFloat}(1,0,0),
                       black(RGBVec{PGFloat}))
    elseif name == "green" || name == "+green"
        set_color_ramp(cmin, cmax, black(RGBVec{PGFloat}),
                       RGBVec{PGFloat}(0,1,0))
    elseif name == "-green"
        set_color_ramp(cmin, cmax, RGBVec{PGFloat}(0,1,0),
                       black(RGBVec{PGFloat}))
    elseif name == "blue" || name == "+blue"
        set_color_ramp(cmin, cmax, black(RGBVec{PGFloat}),
                       RGBVec{PGFloat}(0,0,1))
    elseif name == "-blue"
        set_color_ramp(cmin, cmax, RGBVec{PGFloat}(0,0,1),
                       black(RGBVec{PGFloat}))
    else
        throw_unknown_colormap(name)
    end
end

function palette(lut::AbstractVector{RGBVec{UInt8}}, cmin::Int, cmax::Int)
    length(lut) > 0 || error("no colors!")
    f = 1/255
    I = axes(lut, 1)
    imin, imax = Int(first(I)), Int(last(I))
    if cmin != cmax
        a = (imax - imin)/(cmax - cmin)
        for c in min(cmin,cmax):max(cmin,cmax)
            t = (c - cmin)*a + imin
            i0 = floor(Int, t)
            i1 = min(i0 + 1, imax)
            a1 = t - i0
            a0 = one(a1) - a1
            r = a0*lut[i0].r + a1*lut[i1].r
            g = a0*lut[i0].g + a1*lut[i1].g
            b = a0*lut[i0].b + a1*lut[i1].b
            pgscr(c, f*r, f*g, f*b)
        end
    else
        i = ((imax + imin + 1) >> 1)
        pgscr(c, f*lut[i].r, f*lut[i].g, f*lut[i].b)
    end
end

@noinline throw_unknown_colormap(name::AbstractString) =
    throw(ArgumentError(string("unknown colormap \"", name, "\"")))

get_color_index_range(::Nothing, ::Nothing) = get_color_index_range()

function get_color_index_range()
    cmin, cmax = pgqcir()
    return (Int(cmin), Int(cmax))
end

function get_color_index_range(cmin::Union{Nothing,Integer},
                               cmax::Union{Nothing,Integer})
    qmin, qmax = get_color_index_range()
    rmin, rmax = qmin, qmax
    if cmin !== nothing
        rmin = oftype(rmin, cmin)
    end
    if cmax !== nothing
        rmax = oftype(rmax, cmax)
    end
    if min(rmin, rmax) < qmin || max(rmin, rmax) > qmax
        pgscir(min(rmin, rmax), max(rmin, rmax))
    end
    return (rmin, rmax)
end

"""

```julia
set_color_ramp([cmin::Integer, cmax::Integer,]
               [flag=0 | lo::RGBVec, hi::RGBVec])
```

sets the current colormap with a linear ramp of shades of grays or of colors
interpolated between the background and the foreground color or between two
given colors.

Optional arguments `cmin` and `cmax` are to specify the range for the color
indices to set.  If unspecified, the full range of indices used for images
(cmap1) is modified.  Note that `cmin > cmax` is allowed to reverse the order
of colors.

Optional argument `flag` is an integer.  If the least significant bit of `flag`
is set.  Then the colors are reversed; if the second least significant bit of
`flag` is set, background and the foreground colors are interpolated;
otherwise, black and white colors are interpolated.

Two RGBVec colors, `lo` and `hi`, can be specified instead of `flag` to
interpolate between these two colors.

"""
set_color_ramp(flag::Integer = 0) = set_color_ramp(pgqcir()..., flag)

set_color_ramp(cmin::Integer, cmax::Integer, flag::Integer = 0) =
    set_color_ramp(Int(cmin), Int(cmax), Int(flag))

function set_color_ramp(cmin::Int, cmax::Int, flag::Int = 0)
    if (flag&2) == 1
        # Use background and foreground colors.
        col0 = background(RGBVec{PGFloat})
        col1 = foreground(RGBVec{PGFloat})
    else
        # Force black and white.
        col0 = black(RGBVec{PGFloat})
        col1 = white(RGBVec{PGFloat})
    end
    if (flag&1) == 1
        # Reverse colors.
        col0, col1 = col1, col0
    end
    set_color_ramp(col0, col1, cmin, cmax)
end

set_color_ramp(lo::RGBVec, hi::RGBVec) = set_color_ramp(pgqcir()..., lo, hi)

set_color_ramp(lo::RGBVec, hi::RGBVec, cmin::Integer, cmax::Integer) =
    set_color_ramp(cmin, cmax, lo, hi)

set_color_ramp(cmin::Integer, cmax::Integer, lo::RGBVec, hi::RGBVec) =
    set_color_ramp(Int(cmin), Int(cmax),
                   RGBVec{PGFloat}(lo), RGBVec{PGFloat}(hi))

function set_color_ramp(cmin::Int, cmax::Int,
                        lo::RGBVec{PGFloat}, hi::RGBVec{PGFloat})
    lo = clamp(lo)
    hi = clamp(hi)
    if cmin == cmax
        # Set all color indices to the mean level.
        col = (lo + hi)/2
        for c in min(cmin,cmax):max(cmin,cmax)
            pgscr(c, col)
        end
    else
        # Interpolate the given colors.
        f = one(PGFloat)/PGFloat(cmax - cmin)
        for c in min(cmin,cmax):max(cmin,cmax)
            a1 = (c - cmin)*f
            a0 = one(a1) - a1
            pgscr(c, a0*lo + a1*hi)
        end
    end
end

function set_standard_colors()
    pgscr(0, 0.0,0.0,0.0)
    pgscr(1, 1.0,1.0,1.0)
    pgscr(2, 1.0,0.0,0.0)
    pgscr(3, 0.0,1.0,0.0)
    pgscr(4, 0.0,0.0,1.0)
    pgscr(5, 0.0,1.0,1.0)
    pgscr(6, 1.0,0.0,1.0)
    pgscr(7, 1.0,1.0,0.0)
    pgscr(8, 1.0,0.5,0.0)
    pgscr(9, 0.5,1.0,0.0)
    pgscr(10, 0.0,1.0,0.5)
    pgscr(11, 0.0,0.5,1.0)
    pgscr(12, 0.5,0.0,1.0)
    pgscr(13, 1.0,0.0,0.5)
    pgscr(14, 0.333,0.333,0.333)
    pgscr(15, 0.667,0.667,0.667)
end

end # module


text:
module TestAqua

import Aqua
import Folds
using Test

# Default `Aqua.test_all(Folds)` does not work due to ambiguities
# in upstream packages.

@testset "Method ambiguity" begin
    Aqua.test_ambiguities(Folds)
end

Aqua.test_all(Folds; ambiguities = false)

@testset "Compare test/Project.toml and test/environments/main/Project.toml" begin
    @test Text(read(joinpath(@__DIR__, "Project.toml"), String)) ==
          Text(read(joinpath(@__DIR__, "environments", "main", "Project.toml"), String))
end

end  # module


text:
#= PRETTY-PRINTED TYPES

Magic is used to pretty-print the metric signature as it appears in a type:
`MixedMultivector{(1, 1, 1), S}` displays as `MixedMultivector{⟨+++⟩, S} where S`
=#

"""
	ShowWrapper

Wrapper around an all-bits object which displays as if the object was printed directly.

Examples
===
```jldoctest
julia> GeometricAlgebra.ShowWrapper("⟨+++⟩")
⟨+++⟩

julia> typeof(ans)
GeometricAlgebra.ShowWrapper{Symbol("⟨+++⟩")}
```
"""
struct ShowWrapper{object} end
ShowWrapper(a) = ShowWrapper{Symbol(a)}()
ShowWrapper(a::Symbol) = ShowWrapper(repr(a))

Base.show(io::IO, sw::ShowWrapper{object}) where object = print(io, object)

"""
	pretty_print_type_parameters(io, T, fns...)

Print a `Type` (`DataType` or `UnionAll`) with each parameter pretty-printed
if it is specified. The `i`th type parameter is displayed as what `fns[i]`
returns when called with the parameter.

Examples
===
```jldoctest; setup = :( using GeometricAlgebra: pretty_print_type_parameters )
julia> pretty_print_type_parameters(stdout, Array{T,2} where T,
       	identity,
       	n -> "\$n-dimensional")
Array{T, 2-dimensional} where T
```
"""
function pretty_print_type_parameters(io::IO, T::Type, fns...)
	typevars = TypeVar[]
	while T isa UnionAll # unwrap UnionAll and remember TypeVars
		push!(typevars, T.var)
		T = T.body
	end
	parameters = [p isa TypeVar ? p : ShowWrapper(fn(p)) for (p, fn) ∈ zip(T.parameters, fns)]
	S = T.name.wrapper{parameters...}
	for tv in reverse(typevars) # rewrap into UnionAll
		S = UnionAll(tv, S)
	end
	invoke(show, Tuple{IO,Type}, io, S) # show with default method
end

# pretty-print signature and display bits parameter as short binary string (instead of long hexadecimal number)
function Base.show(io::IO, T::Type{<:Blade})
	dim = isconcretetype(T) ? dimension(T) : 0
	pretty_print_type_parameters(io, T,
		show_signature,
		identity,
		function(bits)
			bits isa Unsigned || return bits
			"0b"*string(bits, pad=dim, base=2)
		end,
		identity
	)
end
Base.show(io::IO, T::Type{<:Multivector}) = pretty_print_type_parameters(io, T, show_signature, identity, identity)
Base.show(io::IO, T::Type{<:MixedMultivector}) = pretty_print_type_parameters(io, T, show_signature, identity)

# when type is displayed on its own, be explicit about pretty-printing to avoid confusing the user
function Base.show(io::IO, ::MIME"text/plain", T::Type{<:AbstractMultivector})
	show(io, T)
	Base.print_without_params(T) && return
	println(io)
	printstyled(io, "(pretty-printed ")
	invoke(show, Tuple{IO,Type}, io, T) # show with default method
	printstyled(io, ")")
end

# when signatures are displayed on their own, show both the shorthand and the full expression
function Base.show(io::IO, ::MIME"text/plain", sig::MetricSignature)
	println(io, show_signature(sig))
	printstyled(io, "(pretty-printed ")
	show(io, sig) # show with default method
	printstyled(io, ")")
end

#= BLADE DISPLAY METHODS =#

# show basis blade only without scalar coefficient
function show_basis_blade(io::IO, sig, bits)
	# blades with higher grade than dimension are always zero,
	#  and do not need to have any basis printed
	0 < grade(bits) <= dimension(sig) || return
	label = basis_blade_label(sig, bits_to_indices(bits))
	printstyled(io, label; bold=true)
end
show_basis_blade(io::IO, b::Blade) = show_basis_blade(io, signature(b), bitsof(b))

plaintext_repr(io, a) = repr(a; context=IOContext(io, :color=>false))

# make things line up right when printing in arrays
function Base.alignment(io::IO, b::Blade)
	(l, r) = Base.alignment(io, b.coeff)
	(l, length(plaintext_repr(io, b)) - l)
end
function Base.alignment(io::IO, b::CompositeMultivector)
	(0, length(plaintext_repr(io, b)))
end

"""
Display blade with parentheses surrounding coefficient if necessary.

```jldoctest
julia> GeometricAlgebra.show_blade(stdout, Blade{(x=1,)}(1 + im, 0b1))
(1 + 1im) x
```
"""
function show_blade(io::IO, b::Blade; compact=false)
	Base.show_unquoted(io, b.coeff, 0, Base.operator_precedence(:*))
	grade(b) == 0 && return
	compact || print(io, " ") # coefficient--basis separator
	(iszero(b) && compact) && print(io, '*')
	show_basis_blade(io, b)
end

#= COMPOSITE MULTIVECTOR DISPLAY METHODS =#

"""
Display a multivector as a column of blades, with coefficients aligned using
the native alignment mechanism, and blades basis aligned.

```jldoctest; setup = :( (x, y, z) = basis((x=1, y=1, z=1)) )
julia> GeometricAlgebra.show_multivector(stdout, 1e3x + y + 1e-3z)
1000.0   x
   1.0   y
   0.001 z
```
"""
function show_multivector(io::IO, a::Multivector; indent=0)
	iszero(a) && return print(io, " "^indent, zero(eltype(a)))

	comps = [(bitsof(u), u.coeff) for u ∈ blades(a)]
	alignments = [Base.alignment(io, coeff) for (_, coeff) ∈ comps]
	L = maximum(first.(alignments))
	R = maximum(last.(alignments))
	isfirst = true
	for ((bits, coeff), (l, r)) ∈ zip(comps, alignments)
		isfirst ? isfirst = false : println(io)
		print(io, " "^(L - l + indent))
		Base.show_unquoted(io, coeff, 0, Base.operator_precedence(:*))
		print(io, " "^(R - r), " ")
		show_basis_blade(io, signature(a), bits)
	end
end

function show_multivector_inline(io::IO, a::Multivector; compact=false, showzeros=false)
	if iszero(a)
		print(io, zero(eltype(a)))
		return
	end
	isfirst = true
	for b ∈ blades(a)
		(!showzeros || compact) && iszero(b) && continue
		isfirst ? isfirst = false : print(io, " + ")
		show_blade(io, b; compact)
	end
end

"""
Display inhomogeneous `MixedMultivector` with each grade on a new line.
"""
function show_mixedmultivector(io::IO, a::MixedMultivector; inline, indent=0)
	firstline = true
	if iszero(a)
		print(io, " "^indent, zero(eltype(a)))
		return
	end
	for k ∈ 0:dimension(a)
		ak = grade(a, k)
		if iszero(ak) continue end
		if firstline
			firstline = false
		else
			print(io, inline ? " + " : "\n")
		end
		print(io, " "^indent)
		showparens = inline && (0 < k < dimension(a))
		showparens && print(io, "(")
		show_multivector_inline(io, ak; compact=inline)
		showparens && print(io, ")")
	end
end

show_header(io::IO, a::MixedMultivector) = println(io, "$(typeof(a)):")
show_header(io::IO, a::HomogeneousMultivector{sig,k}) where {sig,k} = println(io, "Grade-$k $(typeof(a)):")

Base.show(io::IO, a::Blade) = show_blade(io, a; compact=true)
function Base.show(io::IO, ::MIME"text/plain", a::Blade)
	show_header(io, a)
	print(io, " ")
	show_blade(io, a)
end

Base.show(io::IO, a::Multivector) = show_multivector_inline(io, a; compact=true)
function Base.show(io::IO, ::MIME"text/plain", a::Multivector)
	show_header(io, a)
	show_multivector(io, a; indent=1)
end

Base.show(io::IO, a::MixedMultivector) = show_mixedmultivector(io, a; inline=true)
function Base.show(io::IO, ::MIME"text/plain", a::MixedMultivector)
	show_header(io, a)
	show_mixedmultivector(io, a; inline=false, indent=1)
end


text:
"""
    ConfidenceInterval

An estimated interval around an unknown quantity at a given confidence level
"""
struct ConfidenceInterval{T <: AbstractFloat}
    lower::T
    upper::T
    confidence::T
end
lower(interval::ConfidenceInterval) = interval.lower
upper(interval::ConfidenceInterval) = interval.upper
confidence(interval::ConfidenceInterval) = interval.confidence
length(interval::ConfidenceInterval) = upper(interval) - lower(interval)
function in(x::T, interval::ConfidenceInterval) where T <: Real
    isinf(x) && return false
    return interval.lower <= x <= interval.upper
end
function issubset(a::ConfidenceInterval, b::ConfidenceInterval)
    return a.lower >= b.lower && a.upper <= b.upper
end
function show(io::IO, interval::ConfidenceInterval)
    intervalstr = @sprintf "Confidence interval (p = %.0f%%): [%.2f − %.2f]" 100*interval.confidence interval.lower interval.upper
    print(io, intervalstr)
    return io
end
"""
    StochasticSolution

Approximate solution to a stochastic model to a given level of confidence.
"""
struct StochasticSolution{T <: AbstractFloat, A <: AbstractVector}
    x::A
    Q::T
    confidence::ConfidenceInterval{T}

    function StochasticSolution(x::AbstractVector, Q::Real, interval::ConfidenceInterval{T}) where T <: AbstractFloat
        A = typeof(x)
        return new{T, A}(x, Q, interval)
    end
end

function EmptySolution()
    return StochasticSolution([], 0.0, ConfidenceInterval(-Inf, Inf, 1.0))
end
"""
    decision(solution::StochasticSolution)

Return the optimal first stage decision of the approximate `solution` to some stochastic model.
"""
decision(solution::StochasticSolution) = solution.x
"""
    optimal_value(solution::StochasticSolution)

Return the optimal value of the approximate `solution` to some stochastic model
"""
optimal_value(solution::StochasticSolution) = solution.Q
"""
    confidence_interval(solution::StochasticSolution)

Return a confidence interval around the optimal value of the approximate `solution` to some stochastic model
"""
confidence_interval(solution::StochasticSolution) = solution.confidence
function show(io::IO, solution::StochasticSolution)
    println(io, "Stochastic solution")
    println(io, "Optimal value: $(solution.Q)")
    print(io, solution.confidence)
    return io
end


text:
import Main.Distributions:cdf,pdf,quantile
typealias Dist Main.Distributions.Distribution{Main.Distributions.Univariate,Main.Distributions.Continuous}

# quash ambiguity
cdf{S<:Sculpt, Q<:Qualia, C<:Clay}(d::Main.Distributions.Triweight, x::Flex{S,Q,C}) = ErrorException("Not Implemented")
pdf{S<:Sculpt, Q<:Qualia, C<:Clay}(d::Main.Distributions.AbstractMixtureModel{Main.Distributions.Univariate, Main.Distributions.Continuous}, x::FlexFloat.Flex{S,Q,C})= ErrorException("Not Implemented")
pdf{S<:Sculpt, Q<:Qualia, C<:Clay}(d::Main.Distributions.Triweight, x::Flex{S,Q,C}) = ErrorException("Not Implemented")

@inline negabs(x::Clay) = -abs(x)

   for (fn) in (:cdf,) # increasing functions
       @eval begin
           function ($fn){S<:Sculpt, Q<:Qualia, C<:Clay}(d::Dist, x::Flex{S,Q,C})
               loIsOpen, hiIsOpen = boundries(S)
               if loIsOpen
                   with_rounding(C, RoundDown) do
                       lo = ($fn)(d, x.lo)
                   end
               else
                   lo = ($fn)(d, x.lo)
               end
               if hiIsOpen
                   with_rounding(C, RoundUp) do
                       hi = ($fn)(d, x.hi)
                   end
               else
                   hi = ($fn)(d, x.hi)
               end

               if lo>hi
                 loIsOpen, hiIsOpen = hiIsOpen, loIsOpen
                 # swap x.lo,x.hi inline
                 if loIsOpen
                     with_rounding(C, RoundDown) do
                         lo = ($fn)(d, x.hi)
                     end
                 else
                     lo = ($fn)(d, x.hi)
                 end
                 if hiIsOpen
                     with_rounding(C, RoundUp) do
                         hi = ($fn)(d, x.lo)
                     end
                 else
                     hi = ($fn)(d, x.lo)
                 end
               end

               Flex{CLCL,Q,C}(minmax(lo,hi)...)
           end
       end
   end

   for (fn) in (:pdf,) # increasing functions (pdf(-x) == pdf(x), pdf(-x) increases)
       @eval begin
           function ($fn){S<:Sculpt, Q<:Qualia, C<:Clay}(d::Dist, x::Flex{S,Q,C})
               loIsOpen, hiIsOpen = boundries(S)
               if loIsOpen
                   with_rounding(C, RoundDown) do
                       lo = ($fn)(d, negabs(x.lo))
                   end
               else
                   lo = ($fn)(d, negabs(x.lo))
               end
               if hiIsOpen
                   with_rounding(C, RoundUp) do
                       hi = ($fn)(d, negabs(x.hi))
                   end
               else
                   hi = ($fn)(d, negabs(x.hi))
               end

               if lo>hi
                 loIsOpen, hiIsOpen = hiIsOpen, loIsOpen
                 # swap x.lo,x.hi inline
                 if loIsOpen
                     with_rounding(C, RoundDown) do
                         lo = ($fn)(d, x.hi)
                     end
                 else
                     lo = ($fn)(d, x.hi)
                 end
                 if hiIsOpen
                     with_rounding(C, RoundUp) do
                         hi = ($fn)(d, x.lo)
                     end
                 else
                     hi = ($fn)(d, x.lo)
                 end
               end

               Flex{CLCL,Q,C}(minmax(lo,hi)...)
           end
       end
   end

# the implementation of quantile is not ok with directed rounding on 01-Dec-2015
#      roundingDn(quantile(Normal(),0.95) > roundingNearest(Normal(),0.95)
#  and roundingUp(quantile(Normal(),0.95) > roundingNearest(Normal(),0.95)

 for (fn) in (:quantile,) # increasing functions
       @eval begin
           function ($fn){S<:Sculpt, Q<:Qualia, C<:Clay}(d::Dist, x::Flex{S,Q,C})
               loIsOpen, hiIsOpen = boundries(S)
               if loIsOpen
                   loNearest = ($fn)(d, x.lo)
                   with_rounding(C, RoundDown) do
                       lo = min(($fn)(d, x.lo), loNearest)
                   end
               else
                   lo = ($fn)(d, x.lo)
               end
               if hiIsOpen
                   hiNearest = ($fn)(d, x.hi)
                   with_rounding(C, RoundUp) do
                       hi = max(($fn)(d, x.hi), hiNearest)
                   end
               else
                   hi = ($fn)(d, x.hi)
               end

               if lo>hi
                 loIsOpen, hiIsOpen = hiIsOpen, loIsOpen
                 # swap x.lo,x.hi inline
                 if loIsOpen
                     with_rounding(C, RoundDown) do
                         lo = ($fn)(d, x.hi)
                     end
                 else
                     lo = ($fn)(d, x.hi)
                 end
                 if hiIsOpen
                     with_rounding(C, RoundUp) do
                         hi = ($fn)(d, x.lo)
                     end
                 else
                     hi = ($fn)(d, x.lo)
                 end
               end

               Flex{CLCL,Q,C}(minmax(lo,hi)...)
           end
       end
   end


text:
@testset "Filters" begin

@testset "Box filter" begin
    # Box filter always returns 1.0
    B = GRAYTR.BoxFilter(1.0, 2.0)
    @test GRAYTR.evaluate(B, 0.0, 1.0) ≈ 1.0
    @test GRAYTR.evaluate(B, 2.0, -1.0) ≈ 1.0
end

@testset "Triangle filter" begin
    T = GRAYTR.TriangleFilter(1.0, 2.0)
    @test GRAYTR.evaluate(T, 0.0, 0.0) ≈ 2.0
    @test GRAYTR.evaluate(T, 0.5, 0.0) ≈ 1.0
    @test GRAYTR.evaluate(T, 0.5, 1.0) ≈ 0.5
    @test GRAYTR.evaluate(T, 10.0, 0.0) ≈ 0.0
end

end

max_stars_repo_name: dronir/GRAYTR

max_issues_repo_name: dronir/GRAYTR

max_forks_repo_name: dronir/GRAYTR


text:
module Yota

export
    grad,
    gradtape

include("core.jl")

end


text:
s = Dict("output" => Dict("branch_flows" => true), "conv_losses_mp" => true)
@testset "test ac polar opf" begin
    @testset "3-bus case" begin
        result = run_acdcopf("../test/data/case3.m", ACPPowerModel, ipopt_solver; setting = s)

        @test result["termination_status"] == LOCALLY_SOLVED
        @test isapprox(result["objective"], 5907; atol = 1e0)
    end
    @testset "5-bus ac dc case" begin
        result = run_acdcopf("../test/data/case5_acdc.m", ACPPowerModel, ipopt_solver; setting = s)

        @test result["termination_status"] == LOCALLY_SOLVED
        @test isapprox(result["objective"], 194.14; atol = 1e0)
    end
    @testset "5-bus ac dc case with 2 seperate ac grids" begin
        result = run_acdcopf("../test/data/case5_2grids.m", ACPPowerModel, ipopt_solver; setting = s)

        @test result["termination_status"] == LOCALLY_SOLVED
        @test isapprox(result["objective"], 397.36; atol = 1e0)
    end
    @testset "24-bus rts ac dc case with three zones" begin
        result = run_acdcopf("../test/data/case24_3zones_acdc.m", ACPPowerModel, ipopt_solver; setting = s)

    @test isapprox(result["objective"], 150228.15; atol = 1e0)
    end
end

@testset "test dc opf" begin
    @testset "3-bus case" begin
        result = run_acdcopf("../test/data/case3.m", DCPPowerModel, ipopt_solver; setting = s)

        @test result["termination_status"] == LOCALLY_SOLVED
        @test isapprox(result["objective"], 5782.03; atol = 1e0)
    end
    @testset "5-bus ac dc case" begin
        result = run_acdcopf("../test/data/case5_acdc.m", DCPPowerModel, ipopt_solver; setting = s)

        @test result["termination_status"] == LOCALLY_SOLVED
        @test isapprox(result["objective"], 178.314; atol = 1e0)
    end
    @testset "5-bus ac dc case with 2 seperate ac grids" begin
        result = run_acdcopf("../test/data/case5_2grids.m", DCPPowerModel, ipopt_solver; setting = s)

        @test result["termination_status"] == LOCALLY_SOLVED
        @test isapprox(result["objective"], 379.842; atol = 1e0)
    end

    @testset "24-bus rts ac dc case with three zones" begin
       result = run_acdcopf("../test/data/case24_3zones_acdc.m", DCPPowerModel, ipopt_solver; setting = s)

       @test result["termination_status"] == LOCALLY_SOLVED
       @test isapprox(result["objective"], 144791.0; atol = 1e0)
    end
end

@testset "test soc (BIM) opf" begin
    @testset "3-bus case" begin
        result = run_acdcopf("../test/data/case3.m", SOCWRPowerModel, ipopt_solver; setting = s)

        @test result["termination_status"] == LOCALLY_SOLVED
        @test isapprox(result["objective"], 5746.7; atol = 1e0)
    end
    @testset "5-bus ac dc case" begin
        result = run_acdcopf("../test/data/case5_acdc.m", SOCWRPowerModel, ipopt_solver; setting = s)

        @test result["termination_status"] == LOCALLY_SOLVED
        @test isapprox(result["objective"], 183.76; atol = 1e0)
    end
    @testset "5-bus ac dc case with 2 seperate ac grids" begin
        result = run_acdcopf("../test/data/case5_2grids.m", SOCWRPowerModel, ipopt_solver; setting = s)

        @test result["termination_status"] == LOCALLY_SOLVED
        @test isapprox(result["objective"], 363.50; atol = 1e0)
    end
    @testset "24-bus rts ac dc case with three zones" begin
        result = run_acdcopf("../test/data/case24_3zones_acdc.m", SOCWRPowerModel, ipopt_solver; setting = s)

        @test result["termination_status"] == LOCALLY_SOLVED
        @test isapprox(result["objective"], 150156.24; atol = 1e0)
    end
end

@testset "test soc distflow opf_bf" begin
    @testset "3-bus case" begin
        result = run_acdcopf_bf("../test/data/case3.m", SOCBFPowerModel, ipopt_solver; setting = s)

        @test result["termination_status"] == LOCALLY_SOLVED
        @test isapprox(result["objective"], 5746.7; atol = 1e0)
    end
    @testset "5-bus ac dc case" begin
        result = run_acdcopf_bf("../test/data/case5_acdc.m", SOCBFPowerModel, ipopt_solver; setting = s)

        @test result["termination_status"] == LOCALLY_SOLVED
        @test isapprox(result["objective"], 183.91; atol = 1e0)
    end
    @testset "5-bus ac dc case with 2 seperate ac grids" begin
        result = run_acdcopf_bf("../test/data/case5_2grids.m", SOCBFPowerModel, ipopt_solver; setting = s)

        @test result["termination_status"] == LOCALLY_SOLVED
        @test isapprox(result["objective"], 363.50; atol = 1e0)
    end
    @testset "24-bus rts ac dc case with three zones" begin
        result = run_acdcopf_bf("../test/data/case24_3zones_acdc.m", SOCBFPowerModel, ipopt_solver; setting = s)

        @test result["termination_status"] == LOCALLY_SOLVED
        @test isapprox(result["objective"], 150156.26; atol = 1e0)
    end
end

@testset "test qc opf" begin
    @testset "3-bus case" begin
        result = run_acdcopf("../test/data/case3.m", QCRMPowerModel, ipopt_solver; setting = s)

        @test result["termination_status"] == LOCALLY_SOLVED
        @test isapprox(result["objective"], 5780; atol = 1e0)
    end
    @testset "5-bus ac dc case" begin
        result = run_acdcopf("../test/data/case5_acdc.m", QCRMPowerModel, ipopt_solver; setting = s)

        @test result["termination_status"] == LOCALLY_SOLVED
        @test isapprox(result["objective"], 183.76; atol = 1e0)
    end
    @testset "24-bus rts ac dc case with three zones" begin
        result = run_acdcopf("../test/data/case24_3zones_acdc.m", QCRMPowerModel, ipopt_solver; setting = s)

        @test result["termination_status"] == LOCALLY_SOLVED
        @test isapprox(result["objective"], 150156.25; atol = 1e0)
    end
end

@testset "test qc opf with trilinear convexhull relaxation" begin
    @testset "3-bus case" begin
        result = run_acdcopf("../test/data/case3.m", QCRMPowerModel, ipopt_solver; setting = s)

        @test result["termination_status"] == LOCALLY_SOLVED
        @test isapprox(result["objective"], 5780; atol = 1e0)
    end
    @testset "5-bus ac dc case" begin
        result = run_acdcopf("../test/data/case5_acdc.m", QCRMPowerModel, ipopt_solver; setting = s)

        @test result["termination_status"] == LOCALLY_SOLVED
        @test isapprox(result["objective"], 183.76; atol = 1e0)
    end
    @testset "24-bus rts ac dc case with three zones" begin
        result = run_acdcopf("../test/data/case24_3zones_acdc.m", QCRMPowerModel, ipopt_solver; setting = s)

        @test result["termination_status"] == LOCALLY_SOLVED
        @test isapprox(result["objective"],  150156.25; atol = 1e0)
    end
end


text:
using FluxRC, KitBase, Plots, OrdinaryDiffEq, LinearAlgebra
using ProgressMeter: @showprogress

cd(@__DIR__)
ps = UnstructPSpace("square.msh")

ncell = size(ps.cellid, 1)
nface = size(ps.faceType, 1)

ax = -1.0
ay = -1.0

u = zeros(size(ps.cellid, 1))
for i in axes(u, 1)
    u[i] = max(
        exp(-300 * ((ps.cellCenter[i, 1] - 0.5)^2 + (ps.cellCenter[i, 2] - 0.5)^2)),
        1e-4,
    )
end
u0 = deepcopy(u)

cell_normal = zeros(ncell, 3, 2)
for i = 1:ncell
    pids = ps.cellid[i, :]

    cell_normal[i, 1, :] .= unit_normal(ps.points[pids[1], :], ps.points[pids[2], :])
    cell_normal[i, 2, :] .= unit_normal(ps.points[pids[2], :], ps.points[pids[3], :])
    cell_normal[i, 3, :] .= unit_normal(ps.points[pids[3], :], ps.points[pids[1], :])

    p =
        [
            ps.points[pids[1], :] .+ ps.points[pids[2], :],
            ps.points[pids[2], :] .+ ps.points[pids[3], :],
            ps.points[pids[3], :] .+ ps.points[pids[1], :],
        ] / 2

    for j = 1:3
        if dot(cell_normal[i, j, :], p[j][1:2] - ps.cellCenter[i, 1:2]) < 0
            cell_normal[i, j, :] .= -cell_normal[i, j, :]
        end
    end
end

cell_length = zeros(ncell, 3)
for i = 1:ncell
    pids = ps.cellid[i, :]
    cids = ps.cellNeighbors[i, :]

    for j = 1:3
        if cids[j] > 0
            nodes = intersect(pids, ps.cellid[cids[j], :])
            cell_length[i, j] = norm(ps.points[nodes[1], 1:2] - ps.points[nodes[2], 1:2])
        end
    end
end

function flux_ad(fL, fR, uL, uR)
    au = @. (fL - fR) / (uL - uR + 1e-6)
    return @. 0.5 * (fL + fR) #- 0.5 * abs(au) * (uL - uR) # 这个有问题
end

function dudt!(du, u, p, t)
    ps, ax, ay = p

    du .= 0.0

    f = zeros(ncell, 2)
    for i in axes(f, 1)
        if ps.cellType[i] == 0
            f[i, :] .= [ax * u[i], ay * u[i]]

            u1 = u[ps.cellNeighbors[i, 1]]
            u2 = u[ps.cellNeighbors[i, 2]]
            u3 = u[ps.cellNeighbors[i, 3]]

            f1 = [ax * u1, ay * u1]
            f2 = [ax * u2, ay * u2]
            f3 = [ax * u3, ay * u3]

            fa1 = flux_ad(f[i, :], f1, u[i], u1)
            fa2 = flux_ad(f[i, :], f2, u[i], u2)
            fa3 = flux_ad(f[i, :], f3, u[i], u3)

            du[i] =
                -(
                    dot(fa1, cell_normal[i, 1, :]) * cell_length[i, 1] +
                    dot(fa2, cell_normal[i, 2, :]) * cell_length[i, 2] +
                    dot(fa3, cell_normal[i, 3, :]) * cell_length[i, 3]
                )

            du[i] /= ps.cellArea[i]
        end
    end

    return nothing
end

du = zero(u)

p = (ps, ax, ay)
dudt!(du, u, p, 0.0)

tspan = (0.0, 1.0)
#p = ()
prob = ODEProblem(dudt!, u, tspan, p)

dt = 0.01
itg = init(prob, Euler(), save_everystep = false, adaptive = false, dt = dt)

@showprogress for iter = 1:10
    step!(itg)
end

write_vtk(ps.points, ps.cellid, itg.u)


text:
#
# Profile analysis
#

# Returns a separate mean for each ensemble member along the x-axis in the many-columns scenario
ensemble_mean(c::AbstractDataField) = mean(c, dims = (2,3))

column_mean(c::AbstractDataField) = mean(c, dims = 3)

"""
    struct ValueProfileAnalysis{D, A}

A type for doing analyses on a discrepancy profile located
at cell centers. Defaults to taking the mean square difference between
the model and data coarse-grained to the model grid.
"""
struct ValueProfileAnalysis{D, A}
    discrepancy :: D
       analysis :: A
end

ValueProfileAnalysis(grid; analysis=mean) = ValueProfileAnalysis(CenterField(grid), analysis)
ValueProfileAnalysis(; analysis=mean) = ValueProfileAnalysis(nothing, analysis)
on_grid(profile::ValueProfileAnalysis, grid) = ValueProfileAnalysis(grid; analysis=profile.analysis)

"""
    struct GradientProfileAnalysis{D, A}

A type for combining discreprancy between the fields and field gradients.
Defaults to taking the mean square difference between
the model and data coarse-grained to the model grid.
"""
mutable struct GradientProfileAnalysis{D, G, F, W, A}
     ϵ :: D
    ∇ϵ :: G
    ∇ϕ :: F
    gradient_weight :: W
    value_weight :: W
    analysis :: A
end

GradientProfileAnalysis(grid; analysis=mean, gradient_weight=1.0, value_weight=1.0) =
    GradientProfileAnalysis(CenterField(grid), FaceField(grid), FaceField(grid),
                            gradient_weight, value_weight, analysis)

GradientProfileAnalysis(; analysis=mean, gradient_weight=1.0, value_weight=1.0) =
    GradientProfileAnalysis(nothing, nothing, nothing, gradient_weight, value_weight, analysis)

function on_grid(profile::GradientProfileAnalysis, grid)
    return GradientProfileAnalysis(grid;
                                          analysis = profile.analysis,
                                   gradient_weight = profile.gradient_weight,
                                      value_weight = profile.value_weight)
end


text:
# Minesweeper:

mutable struct Field
    size::Tuple{Int, Int}
    numbers::Array{Int, 2}
    possible_mines::Array{Bool, 2}
    actual_mines::Array{Bool, 2}
    visible::Array{Bool, 2}
end

function Field(x, y)
    size = (x, y)
    actual_mines = convert(Array{Bool, 2}, rand(x, y) .< 0.15)
    possible_mines = zeros(Bool, x, y)
    numbers = zeros(Int, x, y)
    visible = zeros(Bool, x, y)
    for i = 1:x
        for j = 1:y
            n = 0
            for di = -1:1
                for dj = -1:1
                    n += (0 < di+i <= x && 0 < dj+j <= y) ?
                         (actual_mines[di+i, dj+j] ? 1 : 0) : 0
                end
            end
            numbers[i, j] = n
        end
    end
    return Field(size, numbers, possible_mines, actual_mines, visible)
end

function printfield(f::Field; showall = false)
    spaces = Int(floor(log(10, f.size[2])))

    str = " "^(4+spaces)
    for i in 1:f.size[1]
        str *= string(" ", i, " ")
    end
    str *= "\n" * " "^(4+spaces) * "___"^f.size[1] * "\n"
    for j = 1:f.size[2]
        str *= " " * string(j) * " "^(floor(log(10, j)) > 0 ? 1 : spaces+1) * "|"
        for i = 1:f.size[1]
            if showall
                str *= f.actual_mines[i, j] ? " * " : "   "
            else
                if f.visible[i, j]
                    str *= " " * string(f.numbers[i, j] > 0 ? f.numbers[i, j] : " ") * " "
                elseif f.possible_mines[i, j]
                    str *= " ? "
                else
                    str *= " . "
                end
            end
        end
        str *= "\r\n"
    end
    println("Found " * string(length(f.possible_mines[f.possible_mines.==true])) *
            " of " * string(length(f.actual_mines[f.actual_mines.==true])) * " mines.\n")
    print(str)
end

function parse_input(str::String)
    input = split(chomp(str), " ")
    mode = input[1]
    println(str)
    coords =  length(input) > 1 ? (parse(Int,input[2]), parse(Int,input[3])) : (0, 0)
    return mode, coords
end

function eval_input(f::Field, str::String)
    mode, coords = parse_input(str)
    (coords[1] > f.size[1] || coords[2] > f.size[2]) && return true
    if mode == "o"
        reveal(f, coords...) || return false
    elseif mode == "m" && f.visible[coords...] == false
        f.possible_mines[coords...] = !f.possible_mines[coords...]
    elseif mode == "close"
        error("You closed the game.")
    end
    return true
end

function reveal(f::Field, x::Int, y::Int)
    (x > f.size[1] || y > f.size[2]) && return true # check for index out of bounds
    f.actual_mines[x, y] && return false # check for mines
    f.visible[x, y] = true
    if f.numbers[x, y] == 0
        for di = -1:1
            for dj = -1:1
                if (0 < di+x <= f.size[1] && 0 < dj+y <= f.size[2]) &&
                        f.actual_mines[x+di, y+dj] == false &&
                        f.visible[x+di, y+dj] == false
                    reveal(f, x+di, y+dj)
                end
            end
        end
    end
    return true
end

function play()
    print("\nWelcome to Minesweeper\n\nEnter the gridsize x  y:\n")
    s = split(readline(), " ")
    f = Field(parse.(Int,s)...)
    won = false
    while true
        printfield(f)
        print("\nWhat do you do? (\"o x y\" to reveal a field; \"m x y\" to toggle a mine; close)\n")
        eval_input(f, readline()) || break
        print("_"^80 * "\n")
        (f.actual_mines == f.possible_mines ||
              f.visible == .!f.actual_mines) && (won = true; break)
    end
    println(won ? "You won the game!" : "You lost the game:\n")
    printfield(f, showall = true)
end

play()

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
using Documenter
using MultivariateTests

makedocs(modules=[MultivariateTests],
         clean=false,
         format=:html,
         sitename="MultivariateTests.jl",
         authors="Alex Arslan and other contributors",
         pages=Any["Home" => "index.md",
                   "Measures of Dispersion" => "dispersion.md",
                   "Partial Correlation" => "partialcor.md",
                   "Equality of Mean Vectors" => "hotelling.md",
                   "Equality of Covariance Matrices" => "covariance.md"])

deploydocs(repo="github.com/ararslan/MultivariateTests.jl.git",
           target="build",
           deps=nothing,
           make=nothing,
           julia="0.6")


text:
using Test
using ClusterManagers
using Distributed

@testset "ElasticManager" begin
    TIMEOUT = 10.

    em = ElasticManager(addr=:auto, port=0)

    # launch worker
    run(`sh -c $(ClusterManagers.get_connect_cmd(em))`, wait=false)

    # wait at most TIMEOUT seconds for it to connect
    @test :ok == timedwait(TIMEOUT) do
        length(em.active) == 1
    end

    wait(rmprocs(workers()))
end

@static if Sys.iswindows()
    windows_which(command) = `powershell.exe -Command Get-Command $command`
    is_lsf_installed() = success(windows_which("bsub.exe"))
    is_sge_installed() = success(windows_which("qsub.exe"))
else
    is_lsf_installed() = success(`which bsub`)
    is_sge_installed() = success(`which qsub`)
end

if is_lsf_installed()

@testset "LSFManager" begin
    p = addprocs_lsf(1, bsub_flags=`-P scicompsoft`)
    @test nprocs() == 2
    @test workers() == p
    @test fetch(@spawnat :any myid()) == p[1]
    @test remotecall_fetch(+,p[1],1,1) == 2
    rmprocs(p)
    @test nprocs() == 1
    @test workers() == [1]
end

end

if is_sge_installed()
  @testset "SGEManager" begin
    p = addprocs_sge(1, queue=``)
    @test nprocs() == 2
    @test workers() == p
    @test fetch(@spawnat :any myid()) == p[1]
    @test remotecall_fetch(+,p[1],1,1) == 2
    rmprocs(p)
    @test nprocs() == 1
    @test workers() == [1]
  end
end


text:
function ueberpruefenormalenVektor(u,v,nx1,nx2,nx3,b3)
println("Eure Richtungsvektoren lauten:");
println("u = ",round(u[1];digits=4));
println("    ",round(u[2];digits=4));
println("    ",round(u[3];digits=4),"\n");
println("v = ",round(v[1];digits=4));
println("    ",round(v[2];digits=4));
println("    ",round(v[3];digits=4),"\n");

u1 = u[1];
u2 = u[2];
u3 = u[3];
v1 = v[1];
v2 = v[2];
v3 = v[3];

x1 = 1;
A = [u2 u3;
     v2 v3];
b = [-u1; -v1];
yz = A\b;
x2 = yz[1];
x3 = yz[2];
fac = nx1/x1;

# Berechnung des Punktes
    M = [nx1 u[1] v[1];
         nx2 u[2] v[2];
         nx3 u[3] v[3]];
    P = M\b3;

if isnan(nx1) && isnan(nx2) && isnan(nx3)
elseif (!isnan(nx1) || !isnan(nx2) || !isnan(nx3)) && (isnan(nx1) || isnan(nx2) || isnan(nx3))
    println("Ihr müsst erst alle NaNs ersetzen, bevor eure Lösung überprüft werden kann.")
elseif norm(fac - nx2/x2) < 0.01 && norm(fac - nx3/x3) < 0.01
    println("Ihr habt den Normalenvektor korrekt berechent! Euer Schnittpunkt P lautet:")
    println("P = ",round(P[1];digits=4));
    println("    ",round(P[2];digits=4));
    println("    ",round(P[3];digits=4),"\n");
else
    println("Überprüft erneut eure Rechnung und Eingabe.")
end

end


text:
function life_rule(old)
    m, n = size(old)
    new = similar(old, m-2, n-2)
    for j = 2:n-1
        for i = 2:m-1
            nc = +(old[i-1,j-1], old[i-1,j], old[i-1,j+1],
                   old[i  ,j-1],             old[i  ,j+1],
                   old[i+1,j-1], old[i+1,j], old[i+1,j+1])

            new[i-1,j-1] = (nc == 3 ? 1 :
                            nc == 2 ? old[i,j] :
                            0)
        end
    end
    new
end

function life_step(d)
    DArray(size(d),procs(d)) do I
        # fetch neighborhood - toroidal boundaries
        top   = mod(first(I[1])-2,size(d,1))+1
        bot   = mod( last(I[1])  ,size(d,1))+1
        left  = mod(first(I[2])-2,size(d,2))+1
        right = mod( last(I[2])  ,size(d,2))+1

        old = Array(Bool, length(I[1])+2, length(I[2])+2)
        @sync begin
            @async old[1      , 1      ] = d[top , left]   # left side
            @async old[2:end-1, 1      ] = d[I[1], left]
            @async old[end    , 1      ] = d[bot , left]
            @async old[1      , 2:end-1] = d[top , I[2]]
            @async old[2:end-1, 2:end-1] = d[I[1], I[2]]   # middle
            @async old[end    , 2:end-1] = d[bot , I[2]]
            @async old[1      , end    ] = d[top , right]  # right side
            @async old[2:end-1, end    ] = d[I[1], right]
            @async old[end    , end    ] = d[bot , right]
        end

        life_rule(old)
    end
end

function plife(m, n)
    w = Window("parallel life", n, m)
    c = Canvas(w)
    pack(c)
    done = false
    c.mouse.button1press = (c,x,y)->(done=true)
    cr = getgc(c)

    grid = DArray(I->convert(Array{Bool,2},randbool(map(length,I))), (m, n), workers())

    last = time(); f = 1
    while !done
        @async begin
            img = convert(Array{Uint32,2},grid) .* 0x00ffffff
            set_source_surface(cr, CairoRGBSurface(img), 0, 0)
            paint(cr)
            reveal(c)
        end
        t = time()
        if t-last > 2
            println("$(f/(t-last)) FPS")
            last = t; f = 0
        end
        grid = life_step(grid)
        f += 1
        sleep(0.01)
    end
end

max_stars_repo_name: jakebolewski/julia

max_issues_repo_name: jakebolewski/julia

max_forks_repo_name: jakebolewski/julia


text:
@snn_kw struct IZParameter{FT=Float32}
    a::FT = 0.01
    b::FT = 0.2
    c::FT = -65
    d::FT = 2
end

@snn_kw mutable struct IZ{VFT=Vector{Float32},VBT=Vector{Bool}}
    param::IZParameter = IZParameter()
    N::Int32 = 100
    v::VFT = fill(-65.0, N)
    u::VFT = param.b * v
    fire::VBT = zeros(Bool, N)
    I::VFT = zeros(N)
    records::Dict = Dict()
end

"""
[Izhikevich Neuron](https://www.izhikevich.org/publications/spikes.htm)
"""
IZ

function integrate!(p::IZ, param::IZParameter, dt::Float32)
    @unpack N, v, u, fire, I = p
    @unpack a, b, c, d = param
    @inbounds for i = 1:N
        v[i] += 0.5f0dt * (0.04f0v[i]^2 + 5f0v[i] + 140f0 - u[i] + I[i])
        v[i] += 0.5f0dt * (0.04f0v[i]^2 + 5f0v[i] + 140f0 - u[i] + I[i])
        u[i] += dt * (a * (b * v[i] - u[i]))
    end
    @inbounds for i = 1:N
        fire[i] = v[i] > 30f0
        v[i] = ifelse(fire[i], c, v[i])
        u[i] += ifelse(fire[i], d, 0f0)
    end
end


text:
module Figures

using Pages, JSON

export figure, Plotly, Style

include("packages/plotly/plotly.jl")

import .Plotly: newPlot

struct Display <: AbstractDisplay end

const current = Ref{String}("")

const Style = Dict{String,String}

function figure(io::IO,id,style=Style())
    id = string(id)
    if tryparse(Int,id) !== nothing
        id = "figure"*id
    end
    global current[] = id

    style = merge(Dict{String,String}(
        "position"      => "absolute",
        "border"        => "thin solid lightgrey",
        "left"          => "50px",
        "top"           => "50px",
        "width"         => "700px",
        "height"        => "450px",
        "box-shadow"    => "3px 3px 5px lightgrey",
        "background"    => "white"),style)

    print(io,"Figures.Figure('$(id)',$(json(style)));")
    return io
end

function closeall(io::IO)
    print(io,"Figures.closeall();")
    return io
end

function close(io::IO,id)
    id = string(id)
    if tryparse(Int,id) !== nothing
        id = "figure"*id
    end
    if id == current[]
        current[] = ""
    end
    print(io,"Figures.close('$(id)');")
    return io
end

for (m,argsin,kwargsin,argsout,kwargsout) in (
        (:figure, (:id, Expr(:kw,:style,:(Style()))), (), (:id,:style), ()),
        (:close, (:id,), (), (:id,), ()),
        (:closeall, (), (), (), ()),
        (:newPlot, (:id, :traces),(
                Expr(:kw,:layout,:(default["layout"])),
                Expr(:kw,:config,:(default["config"]))),
            (:id,:traces),(Expr(:kw,:layout,:layout),Expr(:kw,:config,:config))
        )
    )
    @eval begin
        function $m($(argsin...);$(kwargsin...))
            script = String(take!($m(IOBuffer(),$(argsout...),$(kwargsout...))))
            Pages.broadcast("script",script)
            return
        end
        function $m(route::String,client_id::String,$(argsin...);$(kwargsin...))
            script = String(take!($m(IOBuffer(),$(argsout...),$(kwargsout...))))
            Pages.message(route,client_id,"script",script)
            return
        end
    end
end

start(port=8000) = @async Pages.start(port)
syncstart(port=8000) = Pages.start(port)

include("../examples/examples.jl")

"Download file to a folder relative to /libs"
function download(url::String,path::String)
    Base.download(url,joinpath(@__DIR__,"..","libs",path))
end

function __init__()
    Endpoint("/") do request::HTTP.Request
        read(joinpath(@__DIR__,"index.html"),String)
    end

    Endpoint("/figures.js") do request::HTTP.Request
        read(joinpath(@__DIR__,"..","libs","figures.js"),String)
    end

    Endpoint("/libs/d3/5.9.2/d3.min.js") do request::HTTP.Request
        read(joinpath(@__DIR__,"..","libs","d3","5.9.2","d3.min.js"),String)
    end

    Endpoint("/libs/plotly/1.45.3/plotly.min.js") do request::HTTP.Request
        read(joinpath(@__DIR__,"..","libs","plotly","1.45.3","plotly.min.js"),String)
    end

    pushdisplay(Display())
end

end


text:
import AIBECS: units, @units
import AIBECS: prior, @prior
import AIBECS: limits, @limits
import AIBECS: initial_value, @initial_value
import AIBECS: flattenable, @flattenable
import AIBECS: description, @description
const ∞ = Inf
@initial_value @units @description @limits @flattenable struct TestParameters{T} <: AbstractParameters{T}
    xgeo::T | 2.17 | u"mmol/m^3" | "Geological mean P concentration"             |  (0,∞) | true
    τgeo::T |  1.0 | u"Myr"      | "Geological restoring timescale"              |  (0,∞) | false
    k::T    | 10.0 | u"μmol/m^3" | "Half-saturation constant (Michaelis-Menten)" |  (0,∞) | true
    w₀::T   |  1.0 | u"m/d"      | "Sinking velocity at surface"                 |  (0,∞) | true
    w′::T   | 0.22 | u"m/d/m"    | "Vertical gradient of sinking velocity"       | (-∞,∞) | true
    τDOP::T | 0.25 | u"yr"       | "Remineralization timescale (DOP to DIP)"     |  (0,∞) | true
    τPOP::T | 5.25 | u"d"        | "Dissolution timescale (POP to DOP)"          |  (0,∞) | true
    σ::T    |  1/3 | NoUnits     | "Fraction of quick local uptake recycling"    |  (0,1) | true
    τDIP::T | 30.0 | u"d"        | "Uptake maximum timescale (DIP to POP)"       |  (0,∞) | true
end
function prior(::Type{T}, s::Symbol) where {T<:AbstractParameters}
    if flattenable(T, s)
        lb, ub = limits(T, s)
        if (lb, ub) == (0,∞)
            μ = log(initial_value(T, s))
            return LogNormal(μ, 1.0)
        elseif (lb, ub) == (-∞,∞)
            μ = initial_value(T, s)
            σ = 10.0 # Assumes that a sensible unit is chosen (i.e., that within 10.0 * U)
            return Normal(μ, σ)
        else
            return LocationScale(lb, ub-lb, LogitNormal()) # <- The LogitNormal works well for Optim?
        end
    else
        return nothing
    end
end
prior(::T, s::Symbol) where {T<:AbstractParameters} = prior(T,s)
prior(::Type{T}) where {T<:AbstractParameters} = Tuple(prior(T,s) for s in AIBECS.symbols(T))
prior(::T) where {T<:AbstractParameters} = prior(T)
Distributions.gradlogpdf(::Uniform, x::T) where {T<:Real} = zero(T) # Required for these tests because undefined in Distributions...

p = TestParameters()

@testset "Parameters" begin
    println(p)
    @test p isa AIBECS.AbstractParameters{Float64}
    @testset "Adding a dual-valued vector" begin
        using DualNumbers
        dualp = p + ε * ones(sum(flattenable(p)))
        @test dualp isa AIBECS.AbstractParameters{Dual{Float64}}
        println(dualp)
    end
    @testset "Non-exported functions" begin
        @test AIBECS.symbols(p) isa NTuple{length(fieldnames(typeof(p))), Symbol}
        @test AIBECS.flattenable_symbols(p) isa NTuple{length(p), Symbol}
        @test length(values(p)) == length(fieldnames(typeof(p)))
        @test AIBECS.flattenable_values(p) isa Vector{Float64}
        @test length(AIBECS.flattenable_values(p)) == length(p)
    end
    @testset "Exported functions" begin
        @test values(p) isa Vector{Float64}
        @test vec(p) isa Vector{Float64}
        @test size(p) == size(vec(p))
        @test length(p) == length(vec(p))
        println(latex(p))
    end
    @testset "Unpacking" begin
        @unpack xgeo, τDIP = p
        @test xgeo == ustrip(upreferred(p.xgeo * units(p, :xgeo)))
        @test τDIP == ustrip(upreferred(p.τDIP * units(p, :τDIP)))
    end
    @testset "p <-> λ transformations" begin
        λ = p2λ(p)
        np = length(λ)
        λtest = randn(np)
        @test p2λ(λ2p(TestParameters, λ)) ≈ λ
        @test p2λ(λ2p(TestParameters, λtest)) ≈ λtest
        @testset "$k" for k in AIBECS.flattenable_symbols(p)
            D = prior(p, k)
            @test D isa ContinuousUnivariateDistribution
            lb, ub = limits(p, k)
            b = bijector(D)
            @test b(lb) == -Inf
            @test b(ub) == +Inf
            b⁻¹ = inv(b)
            @test b⁻¹(-Inf) == lb
            @test b⁻¹(+Inf) == ub
        end
    end
end


text:
# Autogenerated wrapper script for Eigen_jll for any
JLLWrappers.@generate_wrapper_header("Eigen")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@generate_init_footer()
end  # __init__()