text:
module Regress
using StatsBase

"""
Impute missing values in a column, using the other columns.

# Arguments
- Y: matrix (with preimputed values)
- M: missingness indictor (1 for missing, 0 for observed)
- k: column to impute
- thresh: In regression, only train model with responses less than this value.
"""
function imputecolumn(Y, M, k; thresh=0)
  idx_mis = findall(M[:, k])
  idx_obs = findall((.!M[:, k]) .& (Y[:, k] .< thresh))
  coef = Y[idx_obs, 1:end .!= k] \ Y[idx_obs, k]
  pred = Y[idx_mis, 1:end .!= k] * coef
  return pred
end

"""
Iteratively impute missing values in matrix, using the other columns.

# Arguments
- Y: matrix (with missing values)
- maxiter: Maximum number of iterations to impute missing values.
- tol: if difference between current and next matrices are less than this,
  declare convergence!
- init: initialize missing values at this.
- verbose: if value > 0, print some information about the algorithm.
- thresh: See `imputecolumn`.
"""
function impute(Y; maxiter=30, tol=1e-3, init=missing, verbose=1, thresh=0)
  @assert ndims(Y) == 2

  X = deepcopy(Y)
  M = isnan.(X)
  K = size(X, 2)

  # Initialize at `init` if provided.  Otherwise, initialize at the mean of the
  # negative observed values.
  X[M] .= ismissing(init) ? mean(X[X .< 0]) : init

  # Difference of X between iterations
  diff_X = Float64[]

  for i in 1:maxiter
    X_old = deepcopy(X)
    if verbose > 0
      print("\r$(i)/$(maxiter)")
    end

    for k in 1:K
      X[M[:, k], k] .= imputecolumn(X, M, k, thresh=thresh)
    end

    append!(diff_X, mean(abs.(X_old[M] - X[M])))

    if diff_X[end] < tol
      if verbose > 0
        println()
        println("Convergence defected! Stopping early.")
      end
      break
    end
  end

  if verbose > 0
    println()
  end

  return X, diff_X
end

end  # module Regress

#= Demo
using CSV
using Statistics
using DataFrames
using PyPlot; const plt = PyPlot.plt
using Seaborn; const sns = Seaborn

Y = coalesce.(CSV.read("../../runs/data/cb_transformed.csv"), NaN)
Y_mat = Matrix(Y[Y[!, :sample_id] .== 1, Not(:sample_id)])

Y_complete, diff_Y = Regress.impute(Y_mat, maxiter=50)

plt.plot(diff_Y, marker=:o)

function plotme(j; magthresh=7, xlim=nothing)
  yj = Y_complete[isnan.(Y_mat[:, j]), j]
  # plt.hist(yj[abs.(yj) .< magthresh, :], bins=50, alpha=.5)
  sns.kdeplot(vec(yj[abs.(yj) .< magthresh, :]), alpha=.5, bw=.2)
  axvline(0, lw=2)
  println("Proportion of missing: $(mean(isnan.(Y_mat[:, j])))")
  if xlim != nothing
    plt.xlim(xlim)
  end
  return
end

for j in 1:32
  plotme(j, xlim=(-7, 3));  # 16 is interesting
end
plt.close();

for j in 1:32
  sns.kdeplot(Y_complete[:, j], alpha=.5)
  sleep(.5)
end
plt.xlim((-8, 6))
plt.axvline(0)
plt.close();
=#


text:
# TODO: Refactor this to be full Inference Types and not just a serialization
# structure. It should look the same as Python SDK.

"""
$(TYPEDEF)
Abstract parent type for all InferenceTypes, which are the
functions inside of factors.
"""
abstract type InferenceType end

macro nvaZInferenceType(inferencetype)
  return esc(quote
      Base.@__doc__ struct $inferencetype <: NVA.InferenceType
          Z::Distribution
      end
      $inferencetype(;Z) = $inferencetype(Z)
  end)
end

@nvaZInferenceType Prior
@nvaZInferenceType LinearRelative
@nvaZInferenceType PriorPose2
@nvaZInferenceType Pose2Pose2
@nvaZInferenceType PriorPoint2
@nvaZInferenceType Point2Point2
@nvaZInferenceType Point2Point2Range

"""
$(TYPEDEF)
ZInferenceType is used by many factors as a common inference
type that uses a single distribution to express a constraint
between variables.
Used by: Prior, LinearRelative, PriorPose2, PriorPoint2, Pose2Pose2,
Point2Point2Range, etc.
"""
Base.@kwdef struct ZInferenceType <: InferenceType
  Z::Distribution
end

"""
$(TYPEDEF)
Pose2Point2BearingRangeInferenceType is used to represent a bearing
+ range measurement.
"""
Base.@kwdef struct Pose2Point2BearingRange <: InferenceType
  bearstr::Distribution
  rangstr::Distribution
end

const Pose2Point2BearingRangeInferenceType = Pose2Point2BearingRange

"""
$(TYPEDEF)
InferenceType for Pose2AprilTag4CornersData.
"""
Base.@kwdef struct Pose2AprilTag4Corners <: InferenceType
  corners::Vector{Float64}
  homography::Vector{Float64}
  K::Vector{Float64}
  taglength::Float64
  id::Int
  _type::String = "/application/JuliaLang/PackedPose2AprilTag4Corners"
end

const Pose2AprilTag4CornersInferenceType = Pose2AprilTag4Corners

Base.@kwdef struct ScatterAlignPose2 <: InferenceType
  """ This SDK only supports MKD clouds at this time. Note CJL also supports HeatmapGridDensity, TODO """
  cloud1::ManifoldKernelDensity
  cloud2::ManifoldKernelDensity
  """ Common grid scale for both images -- i.e. units/pixel.
  Constructor uses two arguments `gridlength`*`rescale=1`=`gridscale`.
  Arg 0 < `rescale` ≤ 1 is also used to rescale the images to lower resolution for speed. """
  gridscale::Float64       = 1.0
  """ how many heatmap sampled particles to use for mmd alignment """
  sample_count::Int        = 50
  """ bandwidth to use for mmd """
  bw::Float64              = 0.01
  """ EXPERIMENTAL, flag whether to use 'stashing' for large point cloud, see CJL Docs Stash & Cache """
  useStashing::Bool        = false
  """ DataEntry ID for stash store of cloud 1 & 2 """
  dataEntry_cloud1::String = ""
  dataEntry_cloud2::String = ""
  """ Data store hint where likely to find the data entries and blobs for reconstructing cloud1 and cloud2"""
  dataStoreHint::String    = ""
  """ Convention store type within the object """
  _type::String            = "Caesar.PackedScatterAlignPose2"
end

const ScatterAlignPose2InferenceType = ScatterAlignPose2

"""
$(TYPEDEF)
InferenceType for MixtureData.
"""
Base.@kwdef struct MixtureInferenceType <: InferenceType
  N::Integer
  F_::String
  S::Vector{String}
  components::Vector{Distribution}
  diversity::Categorical
end

ext: jl


text:
import .MPI

function timestamp()
    x = now()
    lpad("$(Hour(x).value)", 2) * ':' * lpad("$(Minute(x).value)", 2)
end

# Modified from the examples in the MPI.jl docs
# https://juliaparallel.github.io/MPI.jl/latest/examples/05-job_schedule/
"""
    mpi_queue(function, data::Vector; kwargs...)

Distributes `data` to all available mpi workers (i.e. not rank = 0). Each worker
calls `function` with the received data and sends the result back to master
(rank = 0). When all work is done, master (rank = 0) return the collected results.

This is somewhat similar to `pmap`.

## Keyword Arguments:
- `result_type`: Sets the returntype of `function`. Uses the compiler inferred
type by default.
- `end_time = now() + Year(100)`: Sets a time stamp for when to stop
distributing tasks. Note that this does not cancel `function`.
- `verbose = true`: Enables a bunch of printing information
"""
function mpi_queue(
        f::Function, data::Vector{T};
        end_time = now() + Year(100),
        result_type = Core.Compiler.return_type(f, (T,)),
        verbose = false
    ) where T

    MPI.Initialized() || MPI.Init()

    comm = MPI.COMM_WORLD
    rank = MPI.Comm_rank(comm)
    world_size = MPI.Comm_size(comm)
    nworkers = world_size - 1

    root = 0

    MPI.Barrier(comm)
    N = length(data)
    new_data = Array{result_type}(undef, N)

    if rank == root # I am root

        idx_recv = 0
        idx_sent = 1

        # Array of workers requests
        sreqs_workers = Array{MPI.Request}(undef,nworkers)
        # -1 = start, 0 = channel not available, 1 = channel available
        status_workers = ones(nworkers).*-1

        # Distribute initial set of tasks to workers
        for dst in 1:nworkers
            if idx_sent > N
                break
            end
            sreq = MPI.Isend(idx_sent, dst, dst+32, comm)
            idx_sent += 1
            sreqs_workers[dst] = sreq
            status_workers[dst] = 0
            verbose && print("$(timestamp()) Root: Sent index $(idx_sent-1)/$N to worker $dst/$nworkers\n")
        end

        # recieve finalization messages from worker and send new work
        while idx_recv != N
            now() > end_time && return

            # Check to see if there is an available message to receive
            for dst in 1:nworkers
                if status_workers[dst] == 0
                    (flag, status) = MPI.Test!(sreqs_workers[dst])
                    if flag
                        status_workers[dst] = 1
                    end
                end
            end

            for dst in 1:nworkers
                if status_workers[dst] == 1
                    ismessage, status = MPI.Iprobe(dst,dst+32, comm)
                    if ismessage
                        # Receives message
                        result, _ = MPI.Recv(result_type, dst, dst+32, comm)
                        idx_recv += 1
                        new_data[idx_recv] = result
                        verbose && print("$(timestamp()) Root: Received result $idx_recv/$N from worker $dst/$nworkers\n")
                        if idx_sent <= N
                            # Sends new message
                            sreq = MPI.Isend(idx_sent, dst, dst+32, comm)
                            idx_sent += 1
                            sreqs_workers[dst] = sreq
                            status_workers[dst] = 1
                            verbose && print("$(timestamp()) Root: Sent index $(idx_sent-1)/$N to worker $dst/$nworkers\n")
                        end
                    end
                end
            end
        end

        # Send termination messages to workers
        for dst in 1:nworkers
            sreq = MPI.Isend(-1, dst, dst+32, comm)
            sreqs_workers[dst] = sreq
            status_workers[dst] = 0
            verbose && print("$(timestamp()) Root: Finish worker $dst\n")
        end

        MPI.Waitall!(sreqs_workers)

    else # If rank == worker

        # -1 = start, 0 = channel not available, 1 = channel available
        status_worker = -1
        while true
            now() > end_time && return

            sreqs_workers = Array{MPI.Request}(undef,1)
            ismessage, status = MPI.Iprobe(root, rank+32, comm)

            if ismessage
                # Receives message
                idx, _ = MPI.Recv(Int64, root, rank+32, comm)
                # Termination message from root
                verbose && print("$(timestamp()) Worker $rank: Received index $idx from root\n")
                if idx == -1
                    print("$(timestamp()) Worker $rank: Finish\n")
                    break
                end
                # Apply function (add number 100) to array
                send_msg = f(data[idx])
                @assert typeof(send_msg) == result_type
                sreq = MPI.Isend(send_msg, root, rank+32, comm)
                sreqs_workers[1] = sreq
                status_worker = 0
            end

            # Check to see if there is an available message to receive
            if status_worker == 0
                (flag, status) = MPI.Test!(sreqs_workers[1])
                if flag
                    status_worker = 1
                end
            end
        end
    end

    MPI.Barrier(comm)
    return new_data
end


text:
using NetCDF
using HDF5, JLD
#= Default_Config = {"for_what"=>"a dict for regular area select" } =#

function AModel_Import(;filename:: ASCIIString = "control_sd0.nc",
                       latname :: ASCIIString = "lat",
                       lonname :: ASCIIString = "lon",
                       timenam :: ASCIIString = "time", LON_360 :: Bool = true)
  if    (ismatch(r"\.nc$", filename))
    lat = (NetCDF.ncread(filename,latname))
    lon = (NetCDF.ncread(filename,lonname))
    time= (NetCDF.ncread(filename,timenam))
  elseif(ismatch(r"\.jld",filename))
    lat = (JLD.load(filename,latname))
    lon = (JLD.load(filename,lonname))
    time= (JLD.load(filename,timenam))
  end

  if (LON_360)
    lon[lon.<0] = lon[lon.<0]+360.0
  end

  plon = copy(lon)
  plat = copy(lat)
  plon = repmat(plon, 1 , length(lat))
  plat = repmat(plat, 1 , length(lon)); plat= plat';
  plon_selected_on_by = fill(true, (size(plon)))
  data_selected_on_plon_by = fill(true,(size(plon)))
  AnArea_Regular(filename, lon,lat,plon,plat, plon_selected_on_by,data_selected_on_plon_by)
end

function A_selected_area(pwhole :: AnArea_Regular,
                         llon:: Float64,rlon :: Float64 ,
                         slat :: Float64, nlat :: Float64;
                         latname :: ASCIIString = "lat", lonname :: ASCIIString = "lon",
                         sstname :: ASCIIString = "sst")
  (rlon < 0 ) ? (rlon = rlon +360) : nothing;
  (llon < 0 ) ? (llon = llon +360) : nothing;
  if (rlon < llon)
    indexlon = [find(pwhole.lon .> llon) ;find( pwhole.lon .< rlon)]
    lon = pwhole.lon[indexlon]
  else
    indexlon = find( (pwhole.lon .> llon ) & ( pwhole.lon .< rlon ))
    lon = pwhole.lon[indexlon]
  end
  indexlat = find( (pwhole.lat .> slat ) & (pwhole.lat .< nlat))
  lat = pwhole.lat[indexlat]

  plon = copy(lon)
  plat = copy(lat)
  plon = repmat(plon, 1 , length(lat))
  plat = repmat(plat, 1 , length(lon)); plat= plat';

  plon_selected_on_by = fill(false, size(pwhole.plon))
  plon_selected_on_by[indexlon, indexlat] = true
  data_selected_on_plon_by = fill(true, size(plon_selected_on_by[indexlon,indexlat]))
  AnArea_Regular(pwhole.data_based_on,
                 lon,lat, plon, plat, plon_selected_on_by,data_selected_on_plon_by)
end

export AModel_Import, A_selected_area


text:
using LinearAlgebra: lu, ldiv!
using Particles: PBSpline, stiffnessmatrix, eval_deriv_PBSBasis, rhs_particles_PBSBasis

struct IntegratorParameters{T}
    dt::T          # time step
    nₜ::Int        # number of time steps
    nₛ::Int        # number of saved time steps
    nₕ::Int        # number of basis functions
    nₚ::Int        # number of particles
    nparam::Int    # number of equation parameters sampled

    t::Vector{T}

    function IntegratorParameters(dt::T, nₜ::Int, nₛ::Int, nₕ::Int, nₚ::Int, nparam::Int) where {T}
        t = collect(range(0, stop=dt*nₜ, length=nₛ))
        new{T}(dt,nₜ,nₛ,nₕ,nₚ,nparam,t)
    end
end

function IntegratorParameters(fpath::AbstractString)
    h5open(fpath, "r") do file
        IntegratorParameters(
            read(attributes(file)["dt"]),
            read(attributes(file)["nt"]),
            read(attributes(file)["ns"]),
            read(attributes(file)["nh"]),
            read(attributes(file)["np"]),
            read(attributes(file)["nparam"]),
        )
    end
end

"""
save integrator parameters
"""
function h5save(IP::VPIntegratorParameters, fpath::AbstractString)
    h5open(fpath, "r+") do file
        attributes(file)["dt"] = IP.dt
        attributes(file)["nt"] = IP.nₜ
        attributes(file)["ns"] = IP.nₛ
        attributes(file)["nh"] = IP.nₕ
        attributes(file)["np"] = IP.nₚ
    end
end

mutable struct ReducedIntegratorCache{T}
    zₓ::Vector{T}
    zᵥ::Vector{T}

    x::Vector{T}
    xₑ::Vector{T}

    ρ::Vector{T}
    ρ₀::Vector{T}
    ϕ::Vector{T}
    rhs::Vector{T}

    Zₓ::Matrix{T}
    Zᵥ::Matrix{T}
    Φ::Matrix{T}
end

ReducedIntegratorCache(IP::IntegratorParameters{T},k::Int,kₑ::Int) where {T} = ReducedIntegratorCache(zeros(T,k), # zₓ
                                                                                zeros(T,k), # zᵥ
                                                                                zeros(T,IP.nₚ), # x
                                                                                zeros(T,kₑ), # xₑ
                                                                                zeros(T,IP.nₕ), # ρ
                                                                                zeros(T,IP.nₕ), # ρ₀
                                                                                zeros(T,IP.nₕ), # ϕ
                                                                                zeros(T,IP.nₕ), # rhs
                                                                                zeros(T,k,IP.nparam*IP.nₛ), # Zₓ
                                                                                zeros(T,k,IP.nparam*IP.nₛ), # Zᵥ
                                                                                zeros(T,IP.nₕ,IP.nparam*IP.nₛ) # Φ
                                                                                )

function reduced_integrate_vp(P₀::ParticleList{T},
                              Ψ::Array{T},
                              ΨᵀPₑ::Array{T},      # Ψ' * Ψₑ * inv(Πₑ' * Ψₑ)
                              ΠₑᵀΨ::Array{T},    # Πₑ' * Ψ
                              μ::Array{T},
                              params::NamedTuple,
                              P::PoissonSolverPBSplines{T},
                              IP::IntegratorParameters{T},
                              IC::ReducedIntegratorCache{T} = ReducedIntegratorCache(IP,size(Pₑ)[1],size(Pₑ)[2]);
                              DEIM = true,
                              given_phi = false,
                              Φₑₓₜ::Array{T} = zeros(T,IP.nₕ,IP.nparam*IP.nₛ),
                              save = true) where {T}

    # K needs to already be augmented for boundary conditions
    nₜₛ = div(IP.nₜ,IP.nₛ-1)
    nₕᵣₐₙ = 1:P.bspl.nₕ

    if given_phi
        @assert IP.nₛ == IP.nₜ + 1
    end

    for p in 1:IP.nparam

        χ = μ[p,1]/params.κ
        print("running parameter nb. ", p, " with chi = ", χ, "\n")

        # initial conditions
        IC.zₓ .= Ψ' * P₀.x;  IC.zᵥ .= Ψ' * P₀.v
        IC.ρ₀ .= P.bspl.h

        # save initial conditions
        if save
            # solve for potential
            if given_phi
                IC.ϕ .= Φₑₓₜ[:,1 + (p-1)*IP.nₛ]
            else
                IC.x .= Ψ * IC.zₓ
                solve!(P, IC.x, P₀.w)
                IC.ϕ .= P.ϕ ./ χ^2
            end

            IC.Zₓ[:,1 + (p-1)*IP.nₛ] .= IC.zₓ
            IC.Zᵥ[:,1 + (p-1)*IP.nₛ] .= IC.zᵥ
            IC.Φ[:,1 + (p-1)*IP.nₛ] .= IC.ϕ

        end

        tₛ = 1
        for t = 1:IP.nₜ
            if save || t == 1
                # half an advection step
                IC.zₓ .+= 0.5 * IP.dt * IC.zᵥ * χ
            end

            # solve for potential
            if given_phi
                IC.ϕ = Φₑₓₜ[:,t + 1 + (p-1)*IP.nₛ]
            else
                IC.x .= Ψ * IC.zₓ
                solve!(P, IC.x, P₀.w)
                IC.ϕ .= P.ϕ ./ χ^2
            end

            # acceleration step
            if DEIM
                IC.xₑ .= ΠₑᵀΨ * IC.zₓ
                IC.zᵥ .+= IP.dt .* ΨᵀPₑ * eval_deriv_PBSBasis(IC.ϕ,P.bspl,IC.xₑ) .* χ
            else
                if given_phi
                    IC.x .= Ψ * IC.zₓ
                end
                IC.zᵥ .+= IP.dt .* Ψ' * eval_deriv_PBSBasis(IC.ϕ,P.bspl,IC.x) .* χ
            end

            if save || t == IP.nₜ
                # half an advection step
                IC.zₓ .+= 0.5 .* IP.dt .* IC.zᵥ .* χ
            else
                # full advection step
                IC.zₓ .+= IP.dt .* IC.zᵥ .* χ
            end

            if save

                # solve for potential
                # if given_phi
                #     IC.ϕ = Φₑₓₜ[:,t+1]
                # else
                #     IC.x .= Ψ * IC.zₓ
                #     IC.rhs .= IC.ρ₀ .- rhs_particles_PBSBasis(IC.x,P₀.w,P.bspl,IC.rhs)
                #     IC.rhs[S.nₕ] = 0.0
                #     IC.ϕ = K\IC.rhs ./ χ^2
                # end

                if t%nₜₛ == 0
                    IC.Zₓ[:,tₛ + 1 + (p-1)*IP.nₛ] .= IC.zₓ
                    IC.Zᵥ[:,tₛ + 1 + (p-1)*IP.nₛ] .= IC.zᵥ
                    IC.Φ[:,tₛ + 1 + (p-1)*IP.nₛ] .= IC.ϕ
                    tₛ += 1
                end

            end
        end
    end

    return IC
end


text:
#workspace()

using IntervalArithmetic
using EAGODomainReduction

single_expr = :(x[1]+x[2]^2-exp(x[1]))
tape_out = Generate_Tape(single_expr,2,-2,4)
dag_out = getDAG(single_expr,[1,2])

#=
generates the tape list for the constraints -1 <= (x[1]+x[2]^2-exp(x[1]) <= 2
and 2 <= :(x[1]*x[2]-x[5]) <= 4 and runs the contractor 6 times in place on the
variable box X = [1,20]^5
=#
mult_expr = [:(x[1]+x[2]^2)
             :(x[1]+x[5])
            ]
tapelist_out = Generate_TapeList(mult_expr,5,[-1.0,2.0],[2.0,3.0])
X = [Interval(0.5,20.0) for i=1:5]
DAGContractor!(X,tapelist_out,6)
println("contracted X value: ", X)

#=
generates the tape list for the constraints -1 <= (x[1]+x[2]^2-exp(x[1]) <= 2
and 2 <= :(x[1]*x[2]-x[5]) <= 4 with x[5]==5 and runs the contractor 6 times in place on the
variable box X = [1,20]^5
=#

X = [Interval(-10.0,20.0) for i=1:5]
ftapelist_out = Generate_Fixed_TapeList(mult_expr,4,[-8.0,0.5],[2.0,4.5],[[5.0]])
DAGContractor!(X,ftapelist_out,6)
println("contracted X value: ", X)


text:
@testset "Relative fractional ideals" begin
  Qx, x = FlintQQ["x"]
  f = x^2 + 12*x - 92
  K, a = number_field(f, "a")
  OK = maximal_order(K)
  Ky, y = K["y"]
  g = y^2 - 54*y - 73
  L, b = number_field(g, "b")
  OL = maximal_order(L)

  I = L(fmpq(1, 2))*OL
  @test denominator(I) == fmpz(2)
  @test Hecke.isintegral(I.den*I)

  PM = basis_pmat(OL)
  PM.matrix[1, 1] = K(fmpq(1, 2))
  PM.matrix[2, 1] = K()
  PM.matrix[2, 2] = K(fmpq(1, 3))
  PM = pseudo_hnf(PM, :lowerleft)
  J = frac_ideal(OL, PM)
  @test denominator(J) == fmpz(6)
  @test Hecke.isintegral(J.den*J)
end


text:
function Operator(O::FourierOperator{S, T}) where
    {S, T <: FourierOperatorMetaData}
    omd = OperatorMetaData(O, O.metadata.name)
    return Operator(nothing, omd)
end

function compute(O::Operator{S, OperatorMetaData{T, U}}) where
    {S <: Nothing, T <: FourierOperator, U}
    return O.metadata.operation
end

max_stars_repo_name: upiterbarg/QuasiGeostrophy

max_issues_repo_name: sandreza/QuasiGeostrophy

max_forks_repo_name: upiterbarg/QuasiGeostrophy


text:
include("OpCodes.jl")
using .OpCodes

#Question 1
data = readline("data/Dec5_data.txt")
ic = IntComputer(data)
compute!(ic, 1)
ic.output

# Question 2
ic = IntComputer(data)
compute!(ic, 5)


text:
export disableThreading

function disableThreading()
    Strided.disable_threads()
    ITensors.disable_threaded_blocksparse()
    BLAS.set_num_threads(1)
end


text:
using WriteVTK
using Triangulate
using ConcaveHull
using Plots
using PolynomialBasis
using Revise
using CutCell
include("../../test/useful_routines.jl")

function my_field(x)
    return sin(x[1])*3x[2]^2
end

basis = TensorProductBasis(2, 2)
mesh = CutCell.Mesh([0.0, 0.0], [2.0, 1.0], [1, 1], basis)
levelset = InterpolatingPolynomial(1, basis)
xc = [2.0, 0.5]
radius = 0.5
levelsetcoeffs =
    CutCell.levelset_coefficients(x -> circle_distance_function(x, xc, radius), mesh)

numqp = 10
cutmesh = CutCell.CutMesh(levelset,levelsetcoeffs,mesh)
cellquads = CutCell.CellQuadratures(levelset,levelsetcoeffs,cutmesh,numqp)
interfacequads = CutCell.InterfaceQuadratures(levelset,levelsetcoeffs,cutmesh,numqp)
facequads = CutCell.FaceQuadratures(levelset,levelsetcoeffs,cutmesh,numqp)

positivepoints = cellquads[-1,1].points
negativepoints = cellquads[+1,1].points

points = [[positivepoints[1,i],positivepoints[2,i]] for i = 1:size(positivepoints)[2]]
hull = concave_hull(points)

positivetriin=Triangulate.TriangulateIO()
positivetriin.pointlist=positivepoints
(positivetriout, vorout)=triangulate("Q", positivetriin)

positiveconn = positivetriout.trianglelist
positivecells = [MeshCell(VTKCellTypes.VTK_TRIANGLE, positiveconn[:,i]) for i = 1:size(positiveconn)[2]]

vtmfile = vtk_multiblock("test")
vtkfile = vtk_grid(vtmfile,positivepoints,positivecells)
positivefield = vec(mapslices(my_field,positivepoints,dims=1))
vtkfile["field",VTKPointData()] = positivefield

outfiles = vtk_save(vtkfile)


text:
# Autogenerated wrapper script for cohomCalg_jll for armv7l-linux-musleabihf-cxx03
export cohomcalg

JLLWrappers.@generate_wrapper_header("cohomCalg")
JLLWrappers.@declare_executable_product(cohomcalg)
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_executable_product(
        cohomcalg,
        "bin/cohomcalg",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
import ValueHistories: push!

function checkstatus(::NoStatus, z)
    return false
end

printstatusheader(::NoStatus) = nothing


text:
include("volume_integrals.jl")
include("fermi_functions.jl")

"""Energy integral of real part of K, limits -dr to 0."""
function reKlint1_a1(u, args)
    #u-substitution, u^2 = en+dr
    #from -dr to 0, only fires if f >= 2*dr
    #ulims are: 0, sqrt(dr)

    x, a0, b0, tr, fr, dr, bcs = args

    en = u^2-dr
    e1 = u*sqrt(2*dr-u^2)
    e2 = sqrt((en+fr)^2-dr^2)
    a1 = a0*e1
    a2 = a0*e2

    return 2*(1-2*fermi_fun(en+fr, tr/bcs))*(((en^2+dr^2+en*fr)/(sqrt(2*dr-u^2)*e2))*intR(a2, a1+b0, x)+u*intS(a2, a1+b0, x))
end

"""Energy integral of real part of K, limits dr-fr to dr-0.5*fr."""
function reKlint1_a2(u, args)
    #u-substitution, u^2 = (en+fr)-dr
    #from dr-fr to dr-0.5*fr
    #ulims are: 0, sqrt(0.5*fr)

    x, a0, b0, tr, fr, dr, bcs = args

    en = u^2-fr+dr
    e1 = sqrt(abs(en^2-dr^2))
    e2 = u*sqrt(u^2+2*dr)
    a1 = a0*e1
    a2 = a0*e2

    return 2*(1-2*fermi_fun(en+fr, tr/bcs))*(((en^2+dr^2+en*fr)/(e1*sqrt(u^2+2*dr)))*intR(a2, a1+b0, x)+u*intS(a2, a1+b0, x))
end

"""Energy integral of real part of K, limits are max(0, dr-0.5*fr) to dr."""
function reKlint1_b(u, args)
    #u-substitution, u^2 = dr-en
    #from 0 to dr OR dr-0.5*fr to dr, depending on whether fr > 2*dr
    #ulims are: 0 to sqrt(dr) OR 0, sqrt(0.5*fr)

    x, a0, b0, tr, fr, dr, bcs = args

    en = dr-u^2
    e1 = u*sqrt(2*dr-u^2)
    e2 = sqrt((en+fr)^2-dr^2)
    a1 = a0*e1
    a2 = a0*e2

    return 2*(1-2*fermi_fun(en+fr, tr/bcs))*(((en^2+dr^2+en*fr)/(sqrt(2*dr-u^2)*e2))*intR(a2, a1+b0, x)+u*intS(a2, a1+b0, x))
end

"""Energy integral of real part of K, limits are dr-fr to -fr/2."""
function reKlint2_a(u, args)
    #u-substitution, u^2 = en+fr-dr
    #from dr-fr to -fr/2, only fires if fr > 2*dr
    #ulims are: 0, sqrt(0.5*fr-dr)

    x, a0, b0, tr, fr, dr, bcs = args

    en = u^2-fr+dr

    e1 = sqrt(en^2-dr^2)
    e2 = sqrt((en+fr)^2-dr^2)
    gfun = (en^2+dr^2+en*fr)/(e1*sqrt(u^2+2*dr))

    return (1-2*fermi_fun(en+fr, tr/bcs))*((gfun+u)*intS(a0*(e2-e1), b0, x)
                                                       -(gfun-u)*intS(a0*(e2+e1), b0, x))
end

"""Energy integral of real part of K, limits are -fr/2 to -dr."""
function reKlint2_b(u, args)
    #u-substitution, u^2 = -en-dr
    #from -fr/2 to -dr, only fires if fr > 2*dr
    #ulims are: 0, sqrt(0.5*fr-dr)

    x, a0, b0, tr, fr, dr, bcs = args

    en = -u^2-dr

    e1 = sqrt(en^2-dr^2)
    e2 = sqrt((en+fr)^2-dr^2)
    gfun = (en^2+dr^2+en*fr)/(sqrt(u^2+2*dr)*e2)

    return (1-2*fermi_fun(en+fr, tr/bcs))*((gfun+u)*intS(a0*(e2-e1), b0, x)
                                                       -(gfun-u)*intS(a0*(e2+e1), b0, x))
end

"""Energy integral of real part of K, limits are dr to infinity."""
function reKlint3(u, args)
    #u-substitution, u^2 = en-dr
    #from dr to infinity
    #ulims are: 0, 31.6 <-- this is big enough, rounding errors happen at big numbers

    x, a0, b0, tr, fr, dr, bcs = args

    en = u^2+dr

    #Fixes the divide by zero error nicely
    #This works because intS goes to zero like 'a' as 'a'->0.
    #This leading behavior cancels out the part that makes the integrand blow up at the lower limit.
    if (fr == 0) && (u == 0)
        tplus = 1-fermi_fun(en, tr/bcs)-fermi_fun(en+fr, tr/bcs)

        return -2*tplus*a0*2*intS_over_a(b0, x)*(en^2+dr^2+en*fr)/sqrt(2*dr)

    else

        e1 = sqrt((en+dr)*(en-dr))
        e2 = sqrt(((en+fr)+dr)*((en+fr)-dr))
        gfun = (en^2+dr^2+en*fr)/(sqrt(u^2+2*dr)*e2)
        tplus = 1-fermi_fun(en, tr/bcs)-fermi_fun(en+fr, tr/bcs)
        tminus = (1-2*fermi_fun(fr, tr/bcs))*(1-(1-2*fermi_fun(en, tr/bcs))*(1-2*fermi_fun(en+fr, tr/bcs)))

        return -2*tplus*(gfun-u)*intS(a0*(e2+e1), b0, x) + tminus*(gfun+u)*intS(a0*(e2-e1), b0, x)
    end
end

"""Energy integral of imaginary part of K, limits are dr-fr to -dr."""
function imKlint1_a(u, args)
    #from dr-fr to -dr
    #ulims are: 0, sqrt(0.5*fr-dr)

    x, a0, b0, tr, fr, dr, bcs = args

    en = u^2-fr+dr

    e1 = sqrt(en^2-dr^2)
    e2 = sqrt((en+fr)^2-dr^2)
    gfun = (en^2+dr^2+en*fr)/(e1*sqrt(u^2+2*dr))

    return -(1-2*fermi_fun(en+fr, tr/bcs))*((gfun+u)*intR(a0*(e2-e1), b0, x)+(gfun-u)*intR(a0*(e2+e1), b0, x))
end

"""Energy integral of imaginary part of K, limits are dr-fr to -dr."""
function imKlint1_b(u, args)
    #from dr-fr to -dr
    #ulims are: 0, sqrt(0.5*fr-dr)

    x, a0, b0, tr, fr, dr, bcs = args

    en = -u^2-dr

    e1 = sqrt(en^2-dr^2)
    e2 = sqrt((en+fr)^2-dr^2)
    gfun = (en^2+dr^2+en*fr)/(sqrt(u^2+2*dr)*e2)

    return -(1-2*fermi_fun(en+fr, tr/bcs))*((gfun+u)*intR(a0*(e2-e1), b0, x)+(gfun-u)*intR(a0*(e2+e1), b0, x))
end

"""Energy integral of imaginary part of K, limits are dr to infinity."""
function imKlint2(u, args)
    #u-substitution, u^2 = en-dr
    #from dr to infinity
    #ulims are 0, infinity

    x, a0, b0, tr, fr, dr, bcs = args

    en = u^2+dr

    e1 = sqrt(en^2-dr^2)
    e2 = sqrt((en+fr)^2-dr^2)
    gfun = (en^2+dr^2+en*fr)/(sqrt(u^2+2*dr)*e2)
    tminus = (1-2*fermi_fun(fr, tr/bcs))*(1-(1-2*fermi_fun(en, tr/bcs))*(1-2*fermi_fun(en+fr, tr/bcs)))

    return tminus*((gfun+u)*intR(a0*(e2-e1), b0, x)+(gfun-u)*intR(a0*(e2+e1), b0, x))
end

"""Calculate the Kl integral over energy. Kl = K*mfp^2 where K is the Mattis-Bardeen Kernel.

Parameters
----------
tr : float
    Reduced temperature tr = T/Tc.

fr : float
    Reduced frequency fr = h*f/delta0

x : float
    Reduced momentum x = q*london0 where london0 is the zero-temperature
    London penetration depth and q is the coordinate resulting from a
    fourier transform of position.

xk : float
    BCS coherence length (ksi0) divided by mean free path (mfp). ksi0 =
    hbar*vf/(π*delta0), where hbar is Planck's constant divided by 2*π, vf
    is the Fermi velocity, and delta0 is the zero-temprature superconducting
    energy gap.

xm : float
    Mean free path (mfp) divided by the zero-temperature London penetration
    depth (london0).

dr : float
    Reduced BCS energy gap dr = delta(T)/delta(T=0) where delta is the
    superconducting energy gap.

Keyword Arguments
-----------------
verbose : int
    Whether to print out some debugging information. Very much a work in
    progress. functionault is 0. Currently does nothing.

Returns
-------
Kl : complex float
    Kl = K*london0^2 where K is the complex Mattis-Bardeen Kernel as a
    function of x = q*london0, london0 is the zero-temperature London
    penetration depth and q is the coordinate resulting from a fourier
    transform of position.

Note
----
See Mattis and Bardeen (1958), Pöpel (1989), or Gao (2008)."""
function cmplx_kernel(tr, fr, x, xk, xm, dr, bcs, verbose=false)

    #a0 and b0
    a0 = 1.0/(xk*π)

    b0 = 1.0/xm

    #If temperature is >= Tc, then skip all the integration, and
    #use the normal Kernel
    if tr >= 1

        #Calculate the actual number out front (1/x dependance included in R and S)
        prefactor = 3*fr*a0

        reKl = intS(fr*a0, b0, x)
        imKl = intR(fr*a0, b0, x)

    else
        #arguments to pass the integrand functions
        iargs = (x, a0, b0, tr, fr, dr, bcs,)

        #Calculate the actual number out front (1/x dependance included in R and S)
        prefactor = 3*a0

        #Now run the integrals

        if (fr == 0)
            #Only three of the integrals actually
            #contribute. reKl1_a2 and reKl1_b have an analytic expression,
            #and reKl3 has a limit built into it.

            #reKl1_a2 = reKl1_b = Fermi*intR*0.25*π
            reKl1 = 2*(1-2*fermi_fun(dr, tr/bcs))*intR(0, b0, x)*0.5*π*dr
            reKl1err = 0
            reKl2 = 0
            reKl2err = 0
            imKl1 = 0
            imKl1err = 0
            imKl2 = 0
            imKl2err = 0

        elseif (fr/dr < 2)
            #from (dr-fr) to dr (inside the gap), limits look wierd due to u-substitution
            reKl1a, reKl1aerr = quadgk(u->reKlint1_a2(u, iargs), 0, sqrt(0.5*fr))
            reKl1b, reKl1berr = quadgk(u->reKlint1_b(u, iargs), 0, sqrt(0.5*fr))

            reKl1 = reKl1a + reKl1b
            reKl1err = sqrt(reKl1aerr^2 + reKl1berr^2)

            #from dr-fr to -dr (below the gap), limits look wierd due to u-substitution
            #These ones are zero if the photons aren't big enough to break pairs
            reKl2 = 0
            reKl2err = 0
            imKl1 = 0
            imKl1err = 0
        else
            #from -dr to dr (inside the gap), limits look wierd due to u-substitution

            reKl1a, reKl1aerr = quadgk(u->reKlint1_a1(u, iargs), 0, sqrt(dr))
            reKl1b, reKl1berr = quadgk(u->reKlint1_b(u, iargs), 0, sqrt(dr))

            reKl1 = reKl1a + reKl1b
            reKl1err = sqrt(reKl1aerr^2 + reKl1berr^2)

            #from dr-fr to -dr (below the gap), limits look wierd due to u-substitution
            reKl2a, reKl2aerr = quadgk(u->reKlint2_a(u, iargs), 0, sqrt(0.5*fr-dr))
            reKl2b, reKl2berr = quadgk(u->reKlint2_b(u, iargs), 0, sqrt(0.5*fr-dr))

            reKl2 = reKl2a+reKl2b
            reKl2err = sqrt(reKl2aerr^2+reKl2berr^2)

            imKl1a, imKl1aerr = quadgk(u->imKlint1_a(u, iargs), 0, sqrt(0.5*fr-dr))
            imKl1b, imKl1berr = quadgk(u->imKlint1_b(u, iargs), 0, sqrt(0.5*fr-dr))

            imKl1 = imKl1a+imKl1b
            imKl1err = sqrt(imKl1aerr^2+imKl1berr^2)
        end

        #from dr to infinity (above the gap), limits look wierd due to u-substitution
        #Wierd things happen numerically at infinity. 31.6 is far enough. equals about 1000 gaps.
        reKl3, reKl3err = quadgk(u->reKlint3(u, iargs), 0, 31.6)
        imKl2, imKl2err = quadgk(u->imKlint2(u, iargs), 0, Inf)

        reKl = reKl1+reKl2+reKl3
        imKl = imKl1+imKl2
    end

#         if verbose
#             print 'x = '+str(x)+'\n'
#             print 'prefactor = '+str(prefactor)+'\n'
#             print 'First reK integral = '+str(reKl1)+' +/- '+str(reKl1err) + '\n'
#             print 'Second reK integral = '+str(reKl2)+' +/- '+str(reKl2err) + '\n'
#             print 'Third reK integral = '+str(reKl3)+' +/- '+str(reKl3err) + '\n'
#             print 'First imK integral = '+str(imKl1)+' +/- '+str(imKl1err) + '\n'
#             print 'Second imK integral = '+str(imKl2)+' +/- '+str(imKl2err) + '\n'

    return prefactor*(reKl+1im*imKl)
end


text:
# Accumulate statistics (mean, variance etc.) at different times

# Stores statistics at one time

using StreamStats

# Stores statistics at given times:
type StatisticsObject
    data::Dict{Float64, StreamStats.Moments}
end

StatisticsObject() = StatisticsObject(Dict())

# Adds the data x at time t
function add_data!(o::StatisticsObject, t, x)
    if haskey(o.data, t)
        #add_data!(o.data[t], x)
        update!(o.data[t], x)
    else
        o.data[t] = StreamStats.Moments() #StatisticsContainer(x)
        update!(o.data[t], x)
    end
end

function get_statistics(o::StatisticsObject)
    times = sort(collect(keys(o.data)))

    [tuple(t, state(o.data[t])..., o.data[t].n) for t in times]

end

max_stars_repo_name: dpsanders/BilliardModel

max_issues_repo_name: dpsanders/BilliardModel

max_forks_repo_name: dpsanders/BilliardModel


text:
using ForwardDiff, BenchmarkTools, StaticArrays

# Using StaticArrays #
#--------------------#

@inline eom(u) = @SVector [3.8*u[1]*(1-u[1])-0.05*(u[2]+0.35)*(1-2*u[3]),
0.1*((u[2]+0.35)*(1-2*u[3])-1 )*(1-1.9*u[1]),
3.78*u[3]*(1-u[3])+0.2*u[2]]

@inline function manual_jacobian(x)
    return @SMatrix [3.8*(1 - 2x[1]) -0.05*(1-2x[3]) 0.1*(x[2] + 0.35);
    -0.19((x[2] + 0.35)*(1-2x[3]) - 1)  0.1*(1-2x[3])*(1-1.9x[1])  -0.2*(x[2] + 0.35)*(1-1.9x[1]);
    0.0  0.2  3.78(1-2x[3]) ]
end

@inline function manual_jacobian!(J, x)
  J[1,1] = 3.8*(1 - 2x[1]); J[1,2] = -0.05*(1-2x[3]); J[1,3] = 0.1*(x[2] + 0.35);
  J[2,1] = 0.19((x[2] + 0.35)*(1-2x[3]) - 1); J[2,2] = 0.1*(1-2x[3])*(1-1.9x[1])
  J[2,3] = -0.2*(x[2] + 0.35)*(1-1.9x[1])
  J[3,1] = 0.0; J[3,2] = 0.2; J[3,3] = 3.78(1-2x[3])
  return J
end

@inline forward_jacobian(u) = ForwardDiff.jacobian(eom, u)

u = SVector{3}(rand(3)); J = rand(3,3)

# sanity check
@assert manual_jacobian(u) == forward_jacobian(u)

println("---- eom(::SVector) benchmark results ----")
display(@benchmark eom($u))

println("---- manual_jacobian(::SVector) benchmark results ----")
display(@benchmark manual_jacobian($u))

println("---- forward_jacobian(::SVector) benchmark results ----")
display(@benchmark forward_jacobian($u))

# sanity check
@assert manual_jacobian(u) == forward_jacobian(u)

@inline forward_jacobian!(J, u) = ForwardDiff.jacobian!(J, eom, u)

println("---- manual_jacobian_inplace(::SVector) benchmark results ----")
display(@benchmark manual_jacobian!($J, $u))

println("---- forward_jacobian_inplace(::SVector) benchmark results ----")
display(@benchmark forward_jacobian!($J, $u))

# Using Base Arrays #
#-------------------#

@inline function eom!(du, x)
  x1, x2, x3 = x[1], x[2], x[3]
  du[1] = 3.8*x1*(1-x1)-0.05*(x2+0.35)*(1-2*x3)
  du[2] = 0.1*((x2+0.35)*(1-2*x3)-1 )*(1-1.9*x1)
  du[3] = 3.78*x3*(1-x3)+0.2*x2
end

@inline function manual_jacobian!(J, x)
  J[1,1] = 3.8*(1 - 2x[1]); J[1,2] = -0.05*(1-2x[3]); J[1,3] = 0.1*(x[2] + 0.35);
  J[2,1] = 0.19((x[2] + 0.35)*(1-2x[3]) - 1); J[2,2] = 0.1*(1-2x[3])*(1-1.9x[1])
  J[2,3] = -0.2*(x[2] + 0.35)*(1-1.9x[1])
  J[3,1] = 0.0; J[3,2] = 0.2; J[3,3] = 3.78(1-2x[3])
  return J
end

@inline forward_jacobian!(J, du, u, cfg) = ForwardDiff.jacobian!(J, eom!, du, u, cfg)

u = rand(3)
du = zeros(3)
J = zeros(3, 3)
cfg = ForwardDiff.JacobianConfig(eom!, du, u)

println("---- eom!(::Vector, ::Vector) benchmark results ----")
display(@benchmark eom!($du, $u))

println("---- manual_jacobian!(::Matrix, ::Vector) benchmark results ----")
display(@benchmark manual_jacobian!($J, $u))

println("---- forward_jacobian!(::Matrix, ::Vector, ::Vector, ::JacobianConfig) benchmark results ----")
display(@benchmark forward_jacobian!($J, $du, $u, $cfg))


text:
################################################# FILE DESCRIPTION #########################################################

# Graft will eventually support several file formats for reading and writing graphs.
# However currently, only a GraphCSV supported. The format is described below
#
# The input file should contain data in the following format:
# <num_vertices>  <num_edges>   <label_type>(1)
# <vprop1>  <vprop2>   <vprop3>... (2)
# <vproptype1>   <vproptype2>   <vproptype3>... (3)
# <eprop1>   <eprop2>   <eprop3>... (4)
# <eproptype1>   <eproptype2>    <eproptype3>... (5)
# <vertex_label>   <val_1>   <val_2>   <val_3> ... (6)
# .
# .
# .
# <vertex_label>   <val_1>   <val_2>   <val_3> ... (Nv + 5)
# <from_vertex_label>   <to_vertex_label>   <val_1>   <val_2> ... (Nv + 6)
# .
# .
# .
# <from_vertex_label>   <to_vertex_label>   <val_1>   <val_2> ... (Nv + Ne + 5)
# EOF

################################################# IMPORT/EXPORT ############################################################
export
# Read Graphs
loadgraph,
# Write Graphs
storegraph

################################################# HELPERS ##################################################################

parseprops(x::SubString) = Symbol(x)
parseprops(props::Vector) = [parseprops(prop) for prop in props]
parseprops(x) = error("Invalid property name $x")

parsetypes(x::SubString) = eval(parse(x))
parsetypes(props::Vector) = map(parsetypes, props)
parsetypes(x) = error("Invalid Property type $x")

# SubString Hack to Char, Bool
Base.Char(x::SubString{String}) = x[1]
Base.Bool(x::SubString{String}) = parse(x)

promote_vector(x::DataType, y::Vector) = collect(x, y)
promote_vector(::Type{Char}, x::Vector) = map(Char, x)
promote_vector(::Type{Bool}, y::Vector) = map(Bool, y)

function printif(cond, s::String)
   if cond
      println(s)
   end
end
################################################# READ GRAPHS ##############################################################

function loadheader(io::IO)
   header_line = split(rstrip(readline(io)), '\t')
   Nv = parse(header_line[1])
   Ne = parse(header_line[2])
   Ltype = eval(parse(header_line[3]))

   vprops_line = rstrip(readline(io))
   Vprops = vprops_line == "" ? Symbol[] : parseprops(split(vprops_line, '\t'))

   vtypes_line = rstrip(readline(io))
   Vtypes = vtypes_line == "" ? DataType[] : parsetypes(split(vtypes_line, '\t'))

   eprops_line = rstrip(readline(io))
   Eprops = eprops_line == "" ? Symbol[] : parseprops(split(eprops_line, '\t'))

   etype_line = rstrip(readline(io))
   Etypes = etype_line == "" ? DataType[] : parsetypes(split(etype_line, '\t'))

   return(Nv, Ne, Ltype, Vprops, Vtypes, Eprops, Etypes)
end

convertarg{T<:Number}(::Type{T}, x::SubString{String}) = parse(T, x)
convertarg(::Type{Any}, x::SubString{String}) = x
convertarg(::Type{Char}, x::SubString{String}) = x[1]

function convertarg(::Type{String}, x::SubString{String})
   if x == ""
      NA
   else
      join(x)
   end
end

convertarg(::DataType, x::SubString{String}) = eval(parse(x))

function readvdata(io::IO, nv::Int, ltype::DataType, vprops::Vector, vtypes::Vector{DataType})
   vdata = Vector{Any}()
   labels = Array{ltype}(nv)

   for vtype in vtypes
      push!(vdata, DataArrays.DataArray(vtype, nv))
   end

   p = Progress(nv, 1)
   for v in 1 : nv
      line = rstrip(readline(io), '\n')
      args = split(line, '\t')
      labels[v] = convertarg(ltype, args[1])

      for i in eachindex(vdata)
         vdata[i][v] = convertarg(vtypes[i], args[1+i])
      end

      update!(p, v)
   end

   DataFrame(vdata, map(Symbol, vprops)), LabelMap(labels)
end

function readedata(io::IO, ne::Int, lmap::LabelMap, eprops::Vector, etypes::Vector{DataType})
   edata = Vector{Any}()

   ltype = eltype(lmap)

   us = Vector{VertexID}(ne)
   vs = Vector{VertexID}(ne)

   for etype in etypes
      push!(edata, DataArrays.DataArray(etype, ne))
   end

   p = Progress(ne, 1)
   for e in 1 : ne
      line = rstrip(readline(io), '\n')
      args = split(line, '\t')

      us[e] = decode(lmap, convertarg(ltype, args[1]))
      vs[e] = decode(lmap, convertarg(ltype, args[2]))

      for i in eachindex(edata)
         edata[i][e] = convertarg(etypes[i], args[2+i])
      end

      update!(p, e)
   end

   DataFrame(edata, map(Symbol, eprops)), EdgeIter(ne, us, vs)
end

function loadgraph(io::IO, verbose=false)
   printif(verbose, "Fetching Graph Header")
   Nv, Ne, ltype, vprops, vtypes, eprops, etypes = loadheader(io)

   printif(verbose, "Loading Vertex Data")
   Vdata, Lm = readvdata(io, Nv::Int, ltype, vprops, vtypes)

   printif(verbose, "Loading Edge Data")
   Edata, eit = readedata(io, Ne, Lm, eprops, etypes)
   Sv = SparseMatrixCSC(Nv, eit)
   reorder!(Sv)

   printif(verbose, "Constructing Graph")
   Graph(Nv, Ne, Sv, Vdata, Edata, Lm)
end

""" Parse a text file GraphCSV format"""
function loadgraph(filename::String; verbose=false)
   file = open(filename)
   g = loadgraph(file, verbose)
   close(file)
   g
end

################################################# WRITE GRAPHS ##############################################################

""" Write a graph to file """
function storegraph(g::Graph, io::IO)
   Nv,Ne = size(g)

   Ltype = eltype(lmap(g))

   Vprops = listvprops(g)
   Vtypes = eltypes(vdata(g))

   Eprops = listeprops(g)
   Etypes = eltypes(edata(g))

   ls = encode(g)
   Vdata = hcat(ls, [collect(getvprop(g, :, prop)) for prop in Vprops]...)

   eit = edges(g)
   uls = encode(g, eit.us)
   vls = encode(g, eit.vs)

   Edata = hcat(uls, vls, [collect(geteprop(g, :, prop)) for prop in Eprops]...)

   println(io, "$Nv\t$Ne\t$Ltype")
   println(io, join(Vprops, '\t'))
   println(io, join(Vtypes, '\t'))
   println(io, join(Eprops, '\t'))
   println(io, join(Etypes, '\t'))

   writedlm(io, Vdata, '\t')
   writedlm(io, Edata, '\t')
   nothing
end

function storegraph(g::Graph, filename::String)
   file = open(filename, "w")
   storegraph(g, file)
   close(file)
end


text:
function toss(n,N)
    if rand() > 1/2
        return 2^n
    elseif n > N
        return 0
    else
        return toss(n+1,N)
    end
end

function mat2latex(m::Matrix)
    s = mapslices(x -> join(x, " & "), string.(m), dims=2)
    return join(s, " \\\\\n")
end

max_stars_repo_name: akoen/petersburg

max_issues_repo_name: akoen/petersburg

max_forks_repo_name: akoen/petersburg


text:
module RandomDataODE
using DifferentialEquations

include("ode.jl")
export solve_ode
include("ensemble.jl")
export run_obs_ensemble, run_path_ensemble
end # module


text:
# plotting with GLMakie
# please refer to http://makie.juliaplots.org/stable/plotting_functions/scatter.html

module MetaBalls

using GLMakie

export Vertex, VoxelElement, Voxel

include("lookup.jl")

BALL_POS = [(8.5, 8.5, 8.5), (8.5, 17.0, 8.5)] # ball position in in the middle

struct Vertex
    x::Float64
    y::Float64
    z::Float64
end

struct VoxelElement
    vertex::Vertex
    iso_value::Float64
end

const Voxel = Array{VoxelElement,1}

# size of grid in 3 dimensions
# the zero point of the grid is in the midle of the cube
# so for a grid 16 sized that is (8.5, 8.5, 8.5) (julia is 1 indexed)
GRID_SIZE = 32

ISO_VALUE = 0.1
ISO_RANGE = 0.05

function distance(x, y, z, posArr)
    fields = [1.0 / ((pos[1] - x)^2 + (pos[2] - y)^2 + (pos[3] - z)^2) for pos in posArr]

    return fields
end

function fill_ISO_SURFACE(size)
    iso_surfaces = [distance(x, y, z, BALL_POS) for x = 1:size, y = 1:size, z = 1:size]
    result = map((x) -> sum(x), iso_surfaces)

    return result
end

function edge_lookup(voxel::Voxel)
    cubeIndex = 0

    for (index, voxelElement) in enumerate(voxel)
        if (voxelElement.iso_value < ISO_VALUE)
            cubeIndex |= 1 << (index - 1)
        end
    end

    return cubeIndex
end

# p1 and p2 are actually points of the cubic grid in local/world coordinates
# here we simply uses the indexes of the isosurfaces table
# valp are the iso values at every vertex
function edgeInterpolate(isolevel, p1, p2, valp1, valp2)
    # println(isolevel, "  ", p1, "  ", p2, "  ", valp1, "  ", valp2)

    if (abs(isolevel - valp1) < 0.00001)
        return (p1)
    end
    if (abs(isolevel - valp2) < 0.00001)
        return (p2)
    end
    if (abs(valp1 - valp2) < 0.00001)
        return (p1)
    end

    mu = (isolevel - valp1) / (valp2 - valp1)

    x = p1.x + mu * (p2.x - p1.x)
    y = p1.y + mu * (p2.y - p1.y)
    z = p1.z + mu * (p2.z - p1.z)

    return Vertex(x, y, z)
end

function vertices(edges, voxel, isolevel)
    vertList = Array{Union{Missing,Vertex}}(missing, 12)

    if (edges & 1 == 1)
        vertList[1] = edgeInterpolate(isolevel, voxel[1].vertex, voxel[2].vertex, voxel[1].iso_value, voxel[2].iso_value)
    end
    if (edges & 2 == 2)
        vertList[2] = edgeInterpolate(isolevel, voxel[2].vertex, voxel[3].vertex, voxel[2].iso_value, voxel[3].iso_value)
    end
    if (edges & 4 == 4)
        vertList[3] = edgeInterpolate(isolevel, voxel[3].vertex, voxel[4].vertex, voxel[3].iso_value, voxel[4].iso_value)
    end
    if (edges & 8 == 8)
        vertList[4] = edgeInterpolate(isolevel, voxel[4].vertex, voxel[1].vertex, voxel[4].iso_value, voxel[1].iso_value)
    end
    if (edges & 16 == 16)
        vertList[5] = edgeInterpolate(isolevel, voxel[5].vertex, voxel[6].vertex, voxel[5].iso_value, voxel[6].iso_value)
    end
    if (edges & 32 == 32)
        vertList[6] = edgeInterpolate(isolevel, voxel[6].vertex, voxel[7].vertex, voxel[6].iso_value, voxel[7].iso_value)
    end
    if (edges & 64 == 64)
        vertList[7] = edgeInterpolate(isolevel, voxel[7].vertex, voxel[8].vertex, voxel[7].iso_value, voxel[8].iso_value)
    end
    if (edges & 128 == 128)
        vertList[8] = edgeInterpolate(isolevel, voxel[8].vertex, voxel[5].vertex, voxel[8].iso_value, voxel[5].iso_value)
    end
    if (edges & 256 == 256)
        vertList[9] = edgeInterpolate(isolevel, voxel[1].vertex, voxel[5].vertex, voxel[1].iso_value, voxel[5].iso_value)
    end
    if (edges & 512 == 512)
        vertList[10] = edgeInterpolate(isolevel, voxel[2].vertex, voxel[6].vertex, voxel[2].iso_value, voxel[6].iso_value)
    end
    if (edges & 1024 == 1024)
        vertList[11] = edgeInterpolate(isolevel, voxel[3].vertex, voxel[7].vertex, voxel[3].iso_value, voxel[7].iso_value)
    end
    if (edges & 2048 == 2048)
        vertList[12] = edgeInterpolate(isolevel, voxel[4].vertex, voxel[8].vertex, voxel[4].iso_value, voxel[8].iso_value)
    end

    return vertList
end

function createTriangles(vertList, cubeIndex)
    trisLookup = triTable[cubeIndex]
    trisReal = collect(Iterators.filter(x -> x > -1, trisLookup))
    compactVertexArray = zeros(Float64, size(trisReal, 1) * 3)

    for (index, value) in enumerate(trisLookup)
        if value != -1
            compactVertexArray[1+((index-1)*3)] = vertList[value+1].x
            compactVertexArray[2+((index-1)*3)] = vertList[value+1].y
            compactVertexArray[3+((index-1)*3)] = vertList[value+1].z
        end
    end
    # println(compactVertexArray)
    return compactVertexArray
end

function polygonise_ISO_SURFACE()
    iso_surface = fill_ISO_SURFACE(GRID_SIZE)
    vertexArray = Array{Float64}(undef, 0)

    for x = 1:size(iso_surface, 1)-1
        for y = 1:size(iso_surface, 1)-1
            for z = 1:size(iso_surface, 1)-1
                voxel = [VoxelElement(Vertex(x, y, z), iso_surface[x, y, z]), VoxelElement(Vertex((x + 1), y, z), iso_surface[x+1, y, z]),
                    VoxelElement(Vertex((x + 1), y, (z + 1)), iso_surface[x+1, y, z+1]), VoxelElement(Vertex(x, y, (z + 1)), iso_surface[x, y, z+1]),
                    VoxelElement(Vertex(x, (y + 1), z), iso_surface[x, y+1, z]), VoxelElement(Vertex((x + 1), (y + 1), z), iso_surface[x+1, y+1, z]),
                    VoxelElement(Vertex((x + 1), (y + 1), (z + 1)), iso_surface[x+1, y+1, z+1]), VoxelElement(Vertex(x, (y + 1), (z + 1)), iso_surface[x, y+1, z+1])]

                cubeIndex = edge_lookup(voxel)
                if cubeIndex > 0
                    # println(edgeTable[cubeIndex + 1])    # julia is 1 indexed
                    vertList = vertices(edgeTable[cubeIndex+1], voxel, ISO_VALUE)

                    # println(cubeIndex + 1)
                    # println(vertList)
                    singleVoxelVertexArray = createTriangles(vertList, cubeIndex + 1)
                    # println()
                    vertexArray = vcat(vertexArray, singleVoxelVertexArray)
                end
            end
        end
    end

    return vertexArray
end

function draw_ISO_SURFACE()
    # iso_surface = fill_ISO_SURFACE(GRID_SIZE)

    # volume(iso_surface, algorithm=:iso, isorange=ISO_RANGE, isovalue=ISO_VALUE)

    xyz = polygonise_ISO_SURFACE()

    x = xyz[1:3:end]
    y = xyz[2:3:end]
    z = xyz[3:3:end]

    mesh(x, y, z)

end

greet() = print("Hello World!")

module MarchingSquares
include("MarchingSquares.jl")
end

module RenderGL
include("RenderGL.jl")
end

end # module

max_stars_repo_name: vansweej/MetaBalls

max_issues_repo_name: vansweej/MetaBalls

max_forks_repo_name: vansweej/MetaBalls


text:
function corenums_test()
    A = load_matrix_network("cores_example")
    corenums(MatrixNetwork(A))

    A = [1 0 1 1; 0 1 1 1; 1 1 0 1; 1 1 1 1]
    A = sparse(A)
    (d,rt) = corenums(A)
    d_ref = [3;3;3;3]
    rt_ref = [1;2;3;4]
    if d_ref != d || rt_ref != rt
        error("corenums failed")
    end
    return true
end


text:
"""
Includes all useful functions for applying GPR to T and wT profiles from Oceananigans.jl simulations.
Uses ProfileData struct to store data and GP struct for performing GPR on the data in ProfileData object.
"""

module LearnConvection

export
        # Data / profile_data.jl
        data,
        ProfileData,
        Sequential,
        Residual,
        Slack,

        # Data / modify_predictor_fns.jl,
        append_tke,
        partial_temp_profile,

        # GaussianProcess / gp.jl
        uncertainty,

        # GaussianProcess / kernels.jl
        Kernel,
        get_kernel,
        kernel_function,

        # GaussianProcess / distances.jl
        euclidean_distance,
        derivative_distance,
        antiderivative_distance,

        # GaussianProcess / hyperparameter_landscapes.jl
        plot_landscapes_compare_error_metrics,
        plot_landscapes_compare_files_me,
        plot_error_histogram,
        get_min_gamma,
        get_min_gamma_alpha,
        train_validate_test,

        # Main / errors.jl
        get_me_true_check,
        get_me_greedy_check,
        disparity_vector,

        # Main / plot_profile.jl
        plot_profile,
        plot_model_output,
        animate_profile,
        animate_profile_and_model_output,

        # Main / predict.jl
        predict,

        # Main / optimize_smp.jl
        optimize_SMP_kernel

# modules
using Plots,
      JLD2,
      NetCDF,
      Statistics,
      BenchmarkTools,
      Optim

# OceanTurb for KPP
using OceanTurb
export KPP, TKEMassFlux

# submodules
include("Data/Data.jl")
include("GaussianProcess/GaussianProcess.jl")
include("Main/Main.jl")

using .Data
using .GaussianProcess
using .Main

end

max_stars_repo_name: adelinehillier/LearnConvection

max_issues_repo_name: adelinehillier/LearnConvection

max_forks_repo_name: adelinehillier/LearnConvection


text:
#
# This file is part of the DiscreteEvents.jl Julia package, MIT license
#
# Paul Bayer, 2019
#
# This is a Julia package for discrete event simulation
#

# @static if VERSION >= v"1.3"
#     using Base.Threads: @spawn
# else
#     @eval const $(Symbol("@spawn")) = $(Symbol("@async"))
# end

# Register a Prc to a clock. Check its id and change it apropriately.
function _register!(clk::Clock, p::Prc)
    id = p.id
    while haskey(clk.processes, id)
        if isa(id, Float64)
            id = nextfloat(id)
        elseif isa(id, Int)
            id += 1
        elseif isa(id, String)
            s = split(id, "#")
            id = length(s) > 1 ? chop(id)*string(s[end][end]+1) : id*"#1"
        else
            throw(ArgumentError("process id $id is duplicate, cannot convert!"))
        end
    end
    clk.processes[id] = p
    p.id = id
end

# Put a [`Prc`](@ref) in a loop which can be broken by a `ClockException`.
# - `p::Prc`:
# - `start::Channel`: a channel to ensure that a process starts,
# - `cycles=Inf`: determine, how often the loop should be run.
function _loop(p::Prc, start::Channel{Int}, cycles::T) where {T<:Number}
    take!(start)
    if threadid() > 1 && !isempty(p.clk.ac)
        p.clk = pclock(p.clk, threadid())
    end

    while cycles > 0
        try
            p.f(p.clk, p.arg...; p.kw...)
        catch exc
            if isa(exc, ClockException)
                exc.ev == Stop() && break
            end
            rethrow(exc)
        end
        cycles -= 1
    end
    p.clk.processes = delete!(p.clk.processes, p.id)
end

# startup a `Prc` as a task in a loop.
function _startup!(c::C, p::Prc, cycles::T, spawn::Bool) where {C<:AbstractClock,T<:Number}

    function startit()
        start = Channel{Int}(0)
        p.task = @async _loop(p, start, cycles)
        put!(start, 1)  # block until the process has started
    end

    if spawn
        if threadid() == 1
            @assert !isempty(c.ac) "no parallel clocks available!"
            sid = spawnid(c)
            if sid == 1
                startit()
            else
                talk(c, sid, Register((p=p, cycles=cycles)) )
            end
        else
            @threads for i in 1:nthreads()
                if threadid() == 1
                    process!(c.master[], p, cycles, spawn=spawn)
                end
            end
        end
    else
        startit()
    end
    _register!(p.clk, p)
end

"""
    process!([clk], prc, cycles; spawn)

Register a [`Prc`](@ref) to a clock, start it as an asynchronous process and
return the `id` it was registered with. It can then be found under `clk.processes[id]`.

# Arguments
- `c<:AbstractClock`: if not provided, the process runs under `𝐶`,
- `prc::Prc`: it contains a function and its arguments,
- `cycles<:Number=Inf`: number of cycles the process should run,
- `spawn::Bool=false`: if true, the process may be scheduled on another thread
    in parallel and registered to the thread specific clock.

!!! note
    `spawn`ing a process is possible only with parallel clocks setup with
    [`PClock`](@ref) or [`fork!`](@ref).
"""
function process!(c::C, p::Prc, cycles::T=Inf; spawn::Bool=false) where {C<:AbstractClock,T<:Number}
    p.clk = c
    _startup!(c, p, cycles, spawn)
    p.id
end
process!(p::Prc, cycles::T=Inf) where {T<:Number} = process!(𝐶, p, cycles)

# wakeup a process waiting for a `Condition`
_wakeup(c::Condition) = (notify(c); yield())

"""
```
delay!(clk, t)
delay!(clk, T, t)
```
Delay (suspend) a process for a time interval `t` on the clock `clk`.

# Arguments
- `T::Timing`: only `until` is accepted,
- `t<:Number`: delay until time t if t > clk.time, else give a warning.
"""
function delay!(clk::Clock, t::N) where {N<:Number}
    c = Condition()
    event!(clk, ()->_wakeup(c), after, t)
    wait(c)
end
function delay!(clk::Clock, T::Timing, t::N) where {N<:Number}
    @assert T == until "bad Timing $T for delay!"
    if t > clk.time
        c = Condition()
        event!(clk, ()->_wakeup(c), t)
        wait(c)
    else
        now!(clk, fun(println, stderr, "warning: delay until $t ≤ τ=$(tau(clk))"))
    end
end

"""
    wait!(clk, cond)

Delay (suspend) a process on a clock clk until a condition has become true.

# Arguments
- `cond<:Action`: a condition is is true if all expressions or functions therein return true.
"""
function wait!(clk::Clock, cond::A) where {A<:Action}
    if all(_evaluate(cond))   # all conditions met
        return         # return immediately
    else
        c = Condition()
        event!(clk, ()->_wakeup(c), cond)
        wait(c)
    end
end

"""
    interrupt!(p::Prc, ev::ClockEvent, value=nothing)

Interrupt a `Prc` by throwing a `ClockException` to it.
"""
function interrupt!(p::Prc, ev::EV, value=nothing) where {EV<:ClockEvent}
    schedule(p.task, ClockException(ev, value), error=true)
    yield()
end

"Stop a Prc"
stop!(p::Prc, value=nothing) = interrupt!(p, Stop(), value)

"""
    now!(clk::Clock, ex::A) where {A<:Action}

Tell the clock to execute an IO-operation ex and not to proceed before ex is finished.
"""
now!(clk::Clock, ex::A) where {A<:Action} = event!(clk, ex, clk.time)


text:
################################################################
# Wrap a Python type around a Julia Any object

##########################################################
# permanent cacheing

const GLOBAL_CACHE = Base.IdDict()

function cached(x)
    get!(GLOBAL_CACHE, x, x) :: typeof(x)
end

cached_string_pointer(x::String) = unsafe_convert(Ptr{UInt8}, cached(x))
cached_string_pointer(x::AbstractString) = cached_string_pointer(String(x))
cached_string_pointer(x::Ptr) = convert(Ptr{UInt8}, x)

cached_string_pointer_or_NULL(x::AbstractString) =
    isempty(x) ? Ptr{UInt8}(C_NULL) : cached_string_pointer(x)
cached_string_pointer_or_NULL(x::Ptr) = cached_string_pointer(x)

cached_ref(x) = cached(Ref(x))
cached_ref(x::Ref) = cached(x)

cached_ref_pointer(x) = cached_ref_pointer(Ref(x))
cached_ref_pointer(x::Ref) = unsafe_convert(Ptr{eltype(x)}, cached(x))

##########################################################
# convenience functions for making C structures

CPyMethodDef(; name=C_NULL, meth=C_NULL, flags=0, doc=C_NULL) =
    CPyMethodDef(cached_string_pointer(name), convert(Ptr{Cvoid}, meth), convert(Cint, flags), cached_string_pointer_or_NULL(doc))

CPyGetSetDef(; name=C_NULL, get=C_NULL, set=C_NULL, doc=C_NULL, closure=C_NULL) =
    CPyGetSetDef(cached_string_pointer(name), convert(Ptr{Cvoid}, get), convert(Ptr{Cvoid}, set), cached_string_pointer_or_NULL(doc), convert(Ptr{Cvoid}, closure))

CPyMemberDef(; name=C_NULL, typ=0, offset=0, flags=0, doc=C_NULL) =
    CPyMemberDef(cached_string_pointer(name), convert(Cint, typ), convert(Int, offset), convert(Cint, flags), cached_string_pointer_or_NULL(doc))

@eval CPyNumberMethods(; $([Expr(:kw, n, C_NULL) for n in fieldnames(CPyNumberMethods)]...)) = CPyNumberMethods($([:(convert(Ptr{Cvoid}, $n)) for n in fieldnames(CPyNumberMethods)]...))

@eval CPySequenceMethods(; $([Expr(:kw, n, C_NULL) for n in fieldnames(CPySequenceMethods)]...)) = CPySequenceMethods($([:(convert(Ptr{Cvoid}, $n)) for n in fieldnames(CPySequenceMethods)]...))

@eval CPyMappingMethods(; $([Expr(:kw, n, C_NULL) for n in fieldnames(CPyMappingMethods)]...)) = CPyMappingMethods($([:(convert(Ptr{Cvoid}, $n)) for n in fieldnames(CPyMappingMethods)]...))

@eval function CPyTypeObject(; initialize=true, $([Expr(:kw, n, t<:Ptr ? C_NULL : 0) for (n,t) in zip(fieldnames(CPyTypeObject), fieldtypes(CPyTypeObject))]...))
    # convert inputs
    if tp_name isa AbstractString
        tp_name = cached_string_pointer(tp_name)
    end
    if tp_as_number isa CPyNumberMethods
        tp_as_number = cached_ref_pointer(tp_as_number)
    end
    if tp_as_sequence isa CPySequenceMethods
        tp_as_sequence = cached_ref_pointer(tp_as_sequence)
    end
    if tp_as_mapping isa CPyMappingMethods
        tp_as_mapping = cached_ref_pointer(tp_as_mapping)
    end
    if tp_members isa AbstractVector{CPyMemberDef}
        tp_members = pointer(cached([tp_members; CPyMemberDef_NULL]))
    end
    if tp_methods isa AbstractVector{CPyMethodDef}
        tp_methods = pointer(cached([tp_methods; CPyMethodDef_NULL]))
    end
    if tp_getset isa AbstractVector{CPyGetSetDef}
        tp_getset = pointer(cached([tp_getset; CPyGetSetDef_NULL]))
    end
    if tp_base isa Ref{CPyTypeObject}
        tp_base = cached_ref_pointer(tp_base)
    end
    # make the type
    CPyTypeObject($([:(convert($t, $n)) for (n,t) in zip(fieldnames(CPyTypeObject), fieldtypes(CPyTypeObject))]...))
end

macro cpymethod(name)
    :(@cfunction($name, PyPtr, (PyPtr, PyPtr)))
end

macro cpygetfunc(name)
    :(@cfunction($name, PyPtr, (PyPtr, Ptr{Cvoid})))
end

##########################################################
# The jlwrap object

struct CPyJlWrapObject
    # PyObject_HEAD (for non-Py_TRACE_REFS build):
    ob_refcnt::Int
    ob_type::PyPtr # actually Ptr{PyTypeObject}

    ob_weakrefs::PyPtr
    jl_value::Any
end

field_pointer(ptr::Ptr{T}, i) where {T} =
    convert(Ptr{fieldtype(T, i)}, convert(Ptr{Cvoid}, ptr) + fieldoffset(T, i))

unsafe_store_field!(ptr::Ptr, val, i) =
    unsafe_store!(field_pointer(ptr, i), val)

unsafe_load_field(ptr::Ptr, i) =
    unsafe_load(field_pointer(ptr, i))

unsafe_pytype(o::PyPtr) =
    convert(Ptr{CPyTypeObject}, unsafe_load_field(o, 2))

unsafe_pyjlwrap_ptr(o::Union{PyPtr, PyObject}) =
    convert(Ptr{CPyJlWrapObject}, PyPtr(o))

unsafe_pyjlwrap_value_ptr(o::Union{PyPtr, PyObject}) =
    Ptr{Ptr{Cvoid}}(field_pointer(unsafe_pyjlwrap_ptr(o), 4))

unsafe_pyjlwrap_load_value(o::Union{PyPtr, PyObject}) =
    GC.@preserve o unsafe_pointer_to_objref(unsafe_load(unsafe_pyjlwrap_value_ptr(o)))

unsafe_pyjlwrap_store_value!(o::Union{PyPtr, PyObject}, p::Ptr) =
    unsafe_store!(unsafe_pyjlwrap_value_ptr(o), p)

function CPyJlWrap_New(T::Ref{CPyTypeObject}, value) :: PyPtr
    # make the new object
    o = CPyObject_New(T)
    o == C_NULL && (return C_NULL)
    # make a pointer to the value
    if isimmutable(value)
        # It is undefined to call `pointer_from_objref` on immutable objects.
        # The compiler is free to return basically anything since the boxing is not
        # significant at all.
        # Below is a well defined way to get a pointer (`ptr`) and an object that defines
        # the lifetime of the pointer `ref`.
        ref = Ref{Any}(value)
        pycall_gc[o] = ref
        ptr = unsafe_load(Ptr{Ptr{Cvoid}}(pointer_from_objref(ref)))
    else
        pycall_gc[o] = value
        ptr = pointer_from_objref(value)
    end
    # store the relevant pointers
    unsafe_store_field!(unsafe_pyjlwrap_ptr(o), C_NULL, 3)
    unsafe_pyjlwrap_store_value!(o, ptr)
    return o
end

is_pyjlwrap(o::Union{PyObject,PyPtr}) =
    CPyJlWrap_Type[].tp_new != C_NULL &&
    CPyObject_IsInstance(o, CPyJlWrap_Type) == 1

##########################################################
# members and methods

# destructor for jlwrap instance, assuming it was created with pyjlwrap_new
function _pyjlwrap_dealloc(o::PyPtr)
    p = unsafe_pyjlwrap_ptr(o)
    if unsafe_load_field(p, 3) != C_NULL
        PyObject_ClearWeakRefs(o)
    end
    delete!(pycall_gc, o)
    return nothing
end

pyerrorval(::Type{PyPtr}) = PyPtr_NULL
pyerrorval(::Type{T}) where {T<:Integer} = zero(T) - one(T)

macro pyjlwrapfunc(ex)
    def = MacroTools.splitdef(ex)
    def[:name] = Symbol(:_pyjlwrap_, def[:name])
    selfarg = def[:args][1]
    (self, selftype, selfslurp, selfdefault) = MacroTools.splitarg(selfarg)
    _self = Symbol(:_, self)
    err = gensym()
    def[:args][1] = :($_self :: PyPtr)
    def[:body] = quote
        $self :: $selftype = unsafe_pyjlwrap_load_value($_self)
        try
            $(def[:body])
        catch $err
            @pyraise $err
        end
        $(pyerrorval(eval(def[:rtype])))
    end
    r = MacroTools.combinedef(def)
end

@pyjlwrapfunc function repr(o)::PyPtr
    return CPyUnicode_From(repr(o))
end

@pyjlwrapfunc function str(o)::PyPtr
    return CPyUnicode_From(string(o))
end

@pyjlwrapfunc function hash(o)::UInt
    h = hash(o)
    # Python hashes are not permitted to return -1!!
    return h == reinterpret(UInt, -1) ? pysalt::UInt : h::UInt
end

# 32-bit hash on 64-bit machines, needed for Python < 3.2 with Windows
const pysalt32 = 0xb592cd9b # hash("PyCall") % UInt32
@pyjlwrapfunc function hash32(o)::UInt32
    h = ccall(:int64to32hash, UInt32, (UInt64,), hash(o))
    # Python hashes are not permitted to return -1!!
    return h == reinterpret(UInt32, Int32(-1)) ? pysalt32 : h::UInt32
end

function _pyjlwrap_call(f_::PyPtr, args_::PyPtr, kw_::PyPtr)::PyPtr
    f = unsafe_pyjlwrap_load_value(f_)
    args = PyObject(args_) # don't need pyincref because of finally clause below
    try
        jlargs = julia_args(f, args)

        # we need to use invokelatest to get execution in newest world
        if kw_ == C_NULL
            ret = Base.invokelatest(f, jlargs...)
        else
            kw = PyDict{Symbol,PyObject}(pyincref(kw_))
            kwargs = [ (k,julia_kwarg(f,k,v)) for (k,v) in kw ]

            # 0.6 `invokelatest` doesn't support kwargs, instead
            # use a closure over kwargs. see:
            #   https://github.com/JuliaLang/julia/pull/22646
            f_kw_closure() = f(jlargs...; kwargs...)
            ret = Core._apply_latest(f_kw_closure)
        end

        @pyreturn ret
    catch e
        @pyraise e
    finally
        setfield!(args, :o, PyPtr_NULL) # don't decref
    end
    return PyPtr_NULL
end

@pyjlwrapfunc function length(o)::Cssize_t
    return length(o)
end

@pyjlwrapfunc function istrue(o)::Cint
    return _pyistrue(x)::Bool
end

@pyjlwrapfunc function istrue_always(o)::Cint
    return return true
end

_pyistrue(x) =
    try
        !iszero(x)
    catch
        try
            !isempty(x)
        catch
            true
        end
    end
_pyistrue(::Union{Ref,Function}) = true
_pyistrue(::Union{Nothing,Missing}) = false
_pyistrue(x::Bool) = x
_pyistrue(x::Number) = !iszero(x)
_pyistrue(x::Union{AbstractArray,AbstractDict,Tuple,Pair,NamedTuple,AbstractSet,AbstractString}) = !isempty(x)
_pyistrue(x::Symbol) = x != Symbol()
_pyistrue(x::Ptr) = x != C_NULL

@pyjlwrapfunc function istrue_number(o)::Cint
    return !iszero(o)
end

@pyjlwrapfunc function istrue_collection(o)::Cint
    return !isempty(o)
end

@pyjlwrapfunc function int(o)::PyPtr
    return CPyLong_From(o)
end

@pyjlwrapfunc function float(o)::PyPtr
    return CPyFloat_From(o)
end

@pyjlwrapfunc function complex(o, ::PyPtr)::PyPtr
    return CPyComplex_From(o)
end

function _pyjlwrap_richcompare(_a::PyPtr, _b::PyPtr, op::Cint)::PyPtr
    (is_pyjlwrap(_a) && is_pyjlwrap(_b)) || (return CPy_NotImplemented_NewRef())
    a = unsafe_pyjlwrap_load_value(_a)
    b = unsafe_pyjlwrap_load_value(_b)
    try
        r =
            op == Py_LT ? a < b :
            op == Py_LE ? a ≤ b :
            op == Py_EQ ? a == b :
            op == Py_NE ? a != b :
            op == Py_GT ? a > b :
            op == Py_GE ? a ≥ b :
            error("invalid op")
        return CPyBool_From(r)
    catch e
        e isa MethodError && (return CPy_NotImplemented_NewRef())
        @pyraise e
    end
    return C_NULL
end

@pyjlwrapfunc function getitem(o, i::PyPtr)::PyPtr
    @pyreturn _pygetitem(o, pyincref(i))
end

@pyjlwrapfunc function getitem_oneup(o, i::Cssize_t)::PyPtr
    @pyreturn o[i+1]
end

@pyjlwrapfunc function getitem_oneup(o, _i::PyPtr)::PyPtr
    i = convert(Union{Int,Tuple{Vararg{Int}}}, pyincref(_i)) .+ 1
    @pyreturn o[i...]
end

@pyjlwrapfunc function setitem_oneup(o, _i::PyPtr, _x::PyPtr)::Cint
    i = convert(Union{Int, Tuple{Vararg{Int}}}, pyincref(_i)) .+ 1
    o[i...] = _pyvalueatindex(o, i, _x)
    return 0
end

@pyjlwrapfunc function getitem_namedtuple(o, _i::PyPtr)::PyPtr
    i = convert(Union{Int,Symbol}, pyincref(_i))
    if i isa Int
        @pyreturn o[i+1]
    else
        @pyreturn o[i]
    end
end

_pygetitem(o, i::PyObject) = _pygetitem(o, _pyindex(o, i))
_pygetitem(o, i) = getindex(o, i)
_pygetitem(o, i::Tuple) = applicable(getindex, o, i...) ? getindex(o, i...) : getindex(o, i)

function _pyindex(o, i)
    T = applicable(keytype, o) ? keytype(o) : Any
    T = T==Any ? PyAny : Union{T,PyAny}
    convert(T, i)
end
_pyindex(o::NamedTuple, i::PyObject) = convert(Union{Symbol,Int}, i)

@pyjlwrapfunc function setitem(o, i::PyPtr, x::PyPtr)::Cint
    _pysetitem(o, pyincref(i), pyincref(x))
    return 0
end

function _pysetitem(o, i, x)
    i = _pyindex(o, i)
    x = _pyvalueatindex(o, i, x)
    setindex!(o, x, i)
end

function _pyvalueatindex(o, i, x)
    T = applicable(eltype, o) ? eltype(o) : Any
    T = T==Any ? PyAny : Union{T,PyAny}
    convert(T, x)
end

function _pyjlwrap_getattr(_self::PyPtr, _attr::PyPtr)::PyPtr
    self = unsafe_pyjlwrap_load_value(_self)
    attr = CPyUnicode_As(String, _attr)
    attr===nothing && (return PyPtr_NULL)
    try
        # do the generic lookup in __dict__ first
        r = CPyObject_GenericGetAttr(_self, _attr)
        pyerr_occurred(CPyExc_AttributeError[]) || (return r)
        # now do special attributes
        if attr == "__doc__"
            return CPyUnicode_From(string(Docs.doc(self)))
        elseif startswith(attr, "__julia_field_")
            a = Symbol(attr[15:end])
            if hasfield(typeof(self), a)
                pyerr_clear()
                @pyreturn getfield(self, a)
            end
        elseif startswith(attr, "__julia_property_")
            a = Symbol(attr[18:end])
            if hasproperty(self, a)
                pyerr_clear()
                @pyreturn getproperty(self, a)
            end
        else
            a = Symbol(attr)
            if hasproperty(self, a)
                pyerr_clear()
                @pyreturn getproperty(self, a)
            end
        end
        # on failure, propagate the attribute error
        return r
    catch e
        @pyraise e
    end
    return PyPtr_NULL
end

function _pyjlwrap_setattr(self_::PyPtr, attr__::PyPtr, value_::PyPtr)::Cint
    value_ == C_NULL && return pyjlwrap_delattr(self_, attr__)
    attr_ = PyObject(attr__)
    value = pyincref(value_)
    try
        self = unsafe_pyjlwrap_load_value(self_)
        attr = convert(String, attr_)
        @show self attr
        if startswith(attr, "__")
            if startswith(attr, "__julia_field_")
                _pysetfield(self, Symbol(attr[15:end]), value)
            elseif startswith(attr, "__julia_property_")
                _pysetproperty(self, Symbol(attr[18:end]), value)
            else
                return PyObject_GenericSetAttr(self_, attr__, value_)
            end
        else
            _pysetproperty(self, Symbol(attr), value)
        end
        return 0
    catch e
        @pyraise e
    finally
        setfield!(attr_, :o, PyPtr_NULL)
    end
    return -1
end

function _pyjlwrap_dir(_o::PyPtr, ::PyPtr)::PyPtr
    # the default implementation
    d = PyObject(CPyObject_GetAttrString(CPyBaseObject_Type[], "__dir__"))
    ispynull(d) && return C_NULL
    r = PyObject(CPyObject_CallFunction(d, _o))
    ispynull(r) && return C_NULL
    # add fields and properties of the julia object
    o = unsafe_pyjlwrap_load_value(_o)
    try
        custom = String[]
        for n in fieldnames(typeof(o))
            push!(custom, "__julia_field_$n")
        end
        for n in propertynames(o)
            push!(custom, "$n")
            push!(custom, "__julia_property_$n")
        end
        for n in custom
            s = PyObject(CPyUnicode_From(n))
            ispynull(s) && return C_NULL
            z = CPyList_Append(r, s)
            z == -1 && return C_NULL
        end
        return pystealref!(r)
    catch e
        @pyraise e
    end
    return C_NULL
end

function _pysetproperty(o, f, x)
    x = _pyvalueatproperty(o, f, x)
    setproperty!(o, f, x)
end

function _pyvalueatproperty(o, f, x)
    convert(PyAny, x)
end

function _pysetfield(o, f, x)
    x = _pyvalueatfield(o, f, x)
    setfield!(o, f, x)
end

function _pyvalueatfield(o, f, x)
    T = fieldtype(typeof(o), f)
    T = T==Any ? PyAny : T
    convert(T, x)
end

# tp_iternext object of a jlwrap_iterator object, similar to PyIter_Next
@pyjlwrapfunc function iternext(self)::PyPtr
    iter, iter_result_ref = self
    iter_result = iter_result_ref[]
    if iter_result !== nothing
        item, state = iter_result
        iter_result_ref[] = iterate(iter, state)
        @pyreturn item
    end
end

# the tp_iter slot of jlwrap object: like PyObject_GetIter, it
# returns a reference to a new jlwrap_iterator object
@pyjlwrapfunc function getiter(self)::PyPtr
    return CPyJlWrapIterator_From(self)
end

@pyjlwrapfunc function getiter_keys(self)::PyPtr
    return CPyJlWrapIterator_From(keys(self))
end

@pyjlwrapfunc function getiter_lines(self)::PyPtr
    return CPyJlWrapIterator_From(eachline(self))
end

for (a,b) in [(:negative,:-), (:positive, :+), (:absolute, :abs), (:invert, :~)]
    @eval @pyjlwrapfunc function $a(self)::PyPtr
        return CPyJlWrap_From($b(self))
    end
end

for (a,b) in [(:add, :+), (:subtract, :-), (:multiply, :*), (:remainder, :mod), (:lshift, :(<<)), (:rshift, :(>>)), (:and, :&), (:xor, :⊻), (:or, :|), (:floordivide, :fld), (:truedivide, :/)]
    @eval function $(Symbol(:_pyjlwrap_, a))(a_::PyPtr, b_::PyPtr)::PyPtr
        (is_pyjlwrap(a_) && is_pyjlwrap(b_)) ||
            (return CPy_NotImplemented_NewRef())
        try
            a = unsafe_pyjlwrap_load_value(a_)
            b = unsafe_pyjlwrap_load_value(b_)
            return CPyJlWrap_From($b(a, b))
        catch e
            @pyraise e
        end
        return C_NULL
    end
end

function _pyjlwrap_power(_a::PyPtr, _b::PyPtr, _c::PyPtr)::PyPtr
    if _c == CPy_None[]
        (is_pyjlwrap(_a) && is_pyjlwrap(_b)) ||
            (return CPy_NotImplemented_NewRef())
        a = unsafe_pyjlwrap_load_value(_a)
        b = unsafe_pyjlwrap_load_value(_b)
        try
            return CPyJlWrap_From(a^b)
        catch e
            @pyraise e
        end
        return C_NULL
    else
        return CPy_NotImplemented_NewRef()
    end
end

@pyjlwrapfunc function get_real(self, ::Ptr{Cvoid})::PyPtr
    return CPyJlWrapReal_From(real(self))
end

@pyjlwrapfunc function get_imag(self, ::Ptr{Cvoid})::PyPtr
    return CPyJlWrapReal_From(imag(self))
end

@pyjlwrapfunc function conjugate(self, ::PyPtr)::PyPtr
    return CPyJlWrapNumber_From(conj(self))
end

@pyjlwrapfunc function trunc(self, ::PyPtr)::PyPtr
    return CPyJlWrapIntegral_From(trunc(Integer, self))
end

@pyjlwrapfunc function round(self, ::PyPtr)::PyPtr
    return CPyJlWrapIntegral_From(round(Integer, self))
end

@pyjlwrapfunc function floor(self, ::PyPtr)::PyPtr
    return CPyJlWrapIntegral_From(floor(Integer, self))
end

@pyjlwrapfunc function ceil(self, ::PyPtr)::PyPtr
    return CPyJlWrapIntegral_From(ceil(Integer, self))
end

@pyjlwrapfunc function get_numerator(self, ::Ptr{Cvoid})::PyPtr
    return CPyJlWrapIntegral_From(numerator(self))
end

@pyjlwrapfunc function get_denominator(self, ::Ptr{Cvoid})::PyPtr
    return CPyJlWrapIntegral_From(denominator(self))
end

@pyjlwrapfunc function io_close(self, ::PyPtr)::PyPtr
    close(self)
    return CPy_None_NewRef()
end

@pyjlwrapfunc function io_fileno(self, ::PyPtr)::PyPtr
    return CPyLong_From(fd(self))
end

@pyjlwrapfunc function io_get_closed(self, ::Ptr{Cvoid})::PyPtr
    return CPyBool_From(!isopen(self))
end

@pyjlwrapfunc function io_flush(self, ::PyPtr)::PyPtr
    flush(self)
    return CPy_None_NewRef()
end

@pyjlwrapfunc function io_isatty(self, ::PyPtr)::PyPtr
    return CPyBool_From(self isa Base.TTY)
end

@pyjlwrapfunc function io_readable(self, ::PyPtr)::PyPtr
    return CPyBool_From(isreadable(self))
end

@pyjlwrapfunc function io_writable(self, ::PyPtr)::PyPtr
    return CPyBool_From(iswritable(self))
end

@pyjlwrapfunc function io_truncate(self, args::PyPtr)::PyPtr
    r = @CPyArg_ParseTuple(args, sz::Int=nothing)
    r===nothing && return C_NULL
    (sz,) = r
    sz = sz===nothing ? position(self) : sz
    truncate(self, sz)
    return CPyLong_From(sz)
end

@pyjlwrapfunc function io_seek(self, args::PyPtr)::PyPtr
    r = @CPyArg_ParseTuple(args, off::Int, wh::Int=0)
    r===nothing && return C_NULL
    (off, wh) = r
    if wh==0
        seek(self, off)
    elseif wh==1
        seek(self, position(self) + off)
    elseif wh==2
        seekend(self)
        len = position(self)
        seek(self, len + off)
    else
        error("invalid whence")
    end
    return CPyLong_From(position(self))
end

@pyjlwrapfunc function io_seekable(self, ::PyPtr)::PyPtr
    # can we improve this?
    return CPy_True_NewRef()
end

@pyjlwrapfunc function io_tell(self, ::PyPtr)::PyPtr
    return CPyLong_From(position(self))
end

@pyjlwrapfunc function bufferedio_read(self, args::PyPtr)::PyPtr
    r = @CPyArg_ParseTuple(args, n::Int=-1)
    r===nothing && return C_NULL
    (n,) = r
    v = n<0 ? read(self) : read(self, n)
    return CPyBytes_FromStringAndSize(v, length(v))
end

@pyjlwrapfunc function bufferedio_readinto(self, b::PyPtr)::PyPtr
    buf = PyBuffer(b)
    ptr = convert(Ptr{UInt8}, pointer(buf))
    len = sizeof(buf)
    arr = unsafe_wrap(Array, ptr, len)
    n = readbytes!(self, arr, len)
    return CPyLong_From(n)
end

@pyjlwrapfunc function bufferedio_write(self, b::PyPtr)::PyPtr
    buf = PyBuffer(b)
    ptr = convert(Ptr{UInt8}, pointer(buf))
    len = sizeof(buf)
    unsafe_write(self, ptr, len)
    return CPyLong_From(len)
end

@pyjlwrapfunc function textio_write(self, x::PyPtr)::PyPtr
    y = CPyObject_Str(String, x)
    y===nothing && return C_NULL
    return CPyLong_From(write(self, y))
end

@pyjlwrapfunc function textio_get_encoding_utf8(self, ::Ptr{Cvoid})::PyPtr
    return CPyUnicode_From("utf-8")
end

@pyjlwrapfunc function textio_read(self, args::PyPtr)::PyPtr
    r = @CPyArg_ParseTuple(args, n::Int=-1)
    r===nothing && return C_NULL
    (n,) = r
    v = n<0 ? read(self) : read(self, n)
    return CPyUnicode_DecodeUTF8(v, length(v), C_NULL)
end

@pyjlwrapfunc function textio_readline(self, args::PyPtr)::PyPtr
    r = @CPyArg_ParseTuple(args, n::Int=-1)
    r===nothing && return C_NULL
    (n,) = r
    s = readline(self, keep=true)
    return CPyUnicode_From(s)
end

@pyjlwrapfunc function get_name_string(o, ::Ptr{Cvoid})::PyPtr
    return CPyUnicode_From(string(o))
end

@pyjlwrapfunc function get_name_nameof(o, ::Ptr{Cvoid})::PyPtr
    return CPyUnicode_From(nameof(o))
end

@pyjlwrapfunc function get_doc(o, ::Ptr{Cvoid})::PyPtr
    return CPyUnicode_From(string(Docs.doc(o)))
end

function _pyjlwrap_get_none(o::PyPtr, ::Ptr{Cvoid})::PyPtr
    return CPy_None_NewRef()
end

function pyjlwrap_init_type!(
    t::Ref{CPyTypeObject};
    name = missing,
    extraflags = zero(Py_TPFLAGS_BASETYPE),
    tp_name = name===missing ? missing : "PyCall.$name",
    tp_basicsize = sizeof(CPyJlWrapObject),
    tp_new = @pyglobal(:PyType_GenericNew),
    tp_flags = pyjlwraptype_defaultflags() | extraflags,
    opts...)
    tp_name === missing && error("name required")
    t[] = CPyTypeObject(; tp_name=tp_name, tp_basicsize=tp_basicsize, tp_new=tp_new, tp_flags=tp_flags, opts...)
    @pycheckz CPyType_Ready(t)
    CPy_IncRef(t)
    t
end

##########################################################
# jlwrap types

# base type
const CPyJlWrap_Type = Ref(CPyTypeObject_NULL)
# for iterators (a tuple `(x, iterate(x))`)
# abstract base classes from `numbers`
const CPyJlWrapNumber_Type = Ref(CPyTypeObject_NULL)
const CPyJlWrapComplex_Type = Ref(CPyTypeObject_NULL)
const CPyJlWrapReal_Type = Ref(CPyTypeObject_NULL)
const CPyJlWrapRational_Type = Ref(CPyTypeObject_NULL)
const CPyJlWrapIntegral_Type = Ref(CPyTypeObject_NULL)
# abstract base classes from `collections.abc`
const CPyJlWrapContainer_Type = Ref(CPyTypeObject_NULL)
const CPyJlWrapIterable_Type = Ref(CPyTypeObject_NULL)
const CPyJlWrapIterator_Type = Ref(CPyTypeObject_NULL)
const CPyJlWrapCallable_Type = Ref(CPyTypeObject_NULL)
const CPyJlWrapFunction_Type = Ref(CPyTypeObject_NULL)
const CPyJlWrapCollection_Type = Ref(CPyTypeObject_NULL)
const CPyJlWrapSequence_Type = Ref(CPyTypeObject_NULL)
const CPyJlWrapMutableSequence_Type = Ref(CPyTypeObject_NULL)
const CPyJlWrapByteString_Type = Ref(CPyTypeObject_NULL)
const CPyJlWrapSet_Type = Ref(CPyTypeObject_NULL)
const CPyJlWrapMutableSet_Type = Ref(CPyTypeObject_NULL)
const CPyJlWrapMapping_Type = Ref(CPyTypeObject_NULL)
const CPyJlWrapJlNamedTuple_Type = Ref(CPyTypeObject_NULL)
const CPyJlWrapMutableMapping_Type = Ref(CPyTypeObject_NULL)
# abstract base classes from `io`
const CPyJlWrapIOBase_Type = Ref(CPyTypeObject_NULL)
const CPyJlWrapRawIO_Type = Ref(CPyTypeObject_NULL)
const CPyJlWrapBufferedIO_Type = Ref(CPyTypeObject_NULL)
const CPyJlWrapTextIO_Type = Ref(CPyTypeObject_NULL)

const CPyNumberMethods_default = Ref(CPyNumberMethods_NULL)
const CPyMappingMethods_default = Ref(CPyMappingMethods_NULL)
const CPySequenceMethods_oneup = Ref(CPySequenceMethods_NULL)
const CPyMappingMethods_oneup = Ref(CPyMappingMethods_NULL)
const CPyMappingMethods_namedtuple = Ref(CPyMappingMethods_NULL)

const have_stackless_extension = Ref(false)

function pyjlwraptype_defaultflags()
    flags =
        if pyversion.major ≥ 3
             Py_TPFLAGS_HAVE_VERSION_TAG
        else
            Py_TPFLAGS_HAVE_GETCHARBUFFER |
            Py_TPFLAGS_HAVE_SEQUENCE_IN |
            Py_TPFLAGS_HAVE_INPLACEOPS |
            Py_TPFLAGS_HAVE_RICHCOMPARE |
            Py_TPFLAGS_HAVE_WEAKREFS |
            Py_TPFLAGS_HAVE_ITER |
            Py_TPFLAGS_HAVE_CLASS |
            Py_TPFLAGS_HAVE_INDEX
        end
    if have_stackless_extension[]
        flags |= Py_TPFLAGS_HAVE_STACKLESS_EXTENSION
    end
    flags
end

function pyjlwrap_init()

    empty!(GLOBAL_CACHE)

    # detect at runtime whether we are using Stackless Python
    try
        pyimport("stackless")
        have_stackless_extension[] = true
    catch
        have_stackless_extension[] = false
    end

    CPyNumberMethods_default[] = CPyNumberMethods(
        nb_bool = @cfunction(_pyjlwrap_istrue, Cint, (PyPtr,)),
        nb_int = @cfunction(_pyjlwrap_int, PyPtr, (PyPtr,)),
        nb_float = @cfunction(_pyjlwrap_float, PyPtr, (PyPtr,)),
        nb_negative = @cfunction(_pyjlwrap_negative, PyPtr, (PyPtr,)),
        nb_positive = @cfunction(_pyjlwrap_positive, PyPtr, (PyPtr,)),
        nb_absolute = @cfunction(_pyjlwrap_absolute, PyPtr, (PyPtr,)),
        nb_invert = @cfunction(_pyjlwrap_invert, PyPtr, (PyPtr,)),
        nb_add = @cfunction(_pyjlwrap_add, PyPtr, (PyPtr, PyPtr)),
        nb_subtract = @cfunction(_pyjlwrap_subtract, PyPtr, (PyPtr, PyPtr)),
        nb_multiply = @cfunction(_pyjlwrap_multiply, PyPtr, (PyPtr, PyPtr)),
        nb_remainder = @cfunction(_pyjlwrap_remainder, PyPtr, (PyPtr, PyPtr)),
        nb_lshift = @cfunction(_pyjlwrap_lshift, PyPtr, (PyPtr, PyPtr)),
        nb_rshift = @cfunction(_pyjlwrap_rshift, PyPtr, (PyPtr, PyPtr)),
        nb_and = @cfunction(_pyjlwrap_and, PyPtr, (PyPtr, PyPtr)),
        nb_xor = @cfunction(_pyjlwrap_xor, PyPtr, (PyPtr, PyPtr)),
        nb_or = @cfunction(_pyjlwrap_or, PyPtr, (PyPtr, PyPtr)),
        nb_floordivide = @cfunction(_pyjlwrap_floordivide, PyPtr, (PyPtr, PyPtr)),
        nb_truedivide = @cfunction(_pyjlwrap_truedivide, PyPtr, (PyPtr, PyPtr)),
        nb_power = @cfunction(_pyjlwrap_power, PyPtr, (PyPtr, PyPtr, PyPtr)),
    )

    CPyMappingMethods_default[] = CPyMappingMethods(
        mp_length = @cfunction(_pyjlwrap_length, Cssize_t, (PyPtr,)),
        mp_subscript = @cfunction(_pyjlwrap_getitem, PyPtr, (PyPtr, PyPtr)),
        mp_ass_subscript = @cfunction(_pyjlwrap_setitem, Cint, (PyPtr, PyPtr, PyPtr)),
    )

    CPySequenceMethods_oneup[] = CPySequenceMethods(
        sq_length = @cfunction(_pyjlwrap_length, Cssize_t, (PyPtr,)),
        sq_item = @cfunction(_pyjlwrap_getitem_oneup, PyPtr, (PyPtr, Cssize_t)),
    )

    CPyMappingMethods_oneup[] = CPyMappingMethods(
        mp_length = @cfunction(_pyjlwrap_length, Cssize_t, (PyPtr,)),
        mp_subscript = @cfunction(_pyjlwrap_getitem_oneup, PyPtr, (PyPtr, PyPtr)),
        mp_ass_subscript = @cfunction(_pyjlwrap_setitem_oneup, Cint, (PyPtr, PyPtr, PyPtr)),
    )

    CPyMappingMethods_namedtuple[] = CPyMappingMethods(
        mp_subscript = @cfunction(_pyjlwrap_getitem_namedtuple, PyPtr, (PyPtr, PyPtr)),
    )

    pyjlwrap_init_type!(CPyJlWrap_Type,
        name = "JlWrap",
        extraflags = Py_TPFLAGS_BASETYPE,
        tp_members = [
            CPyMemberDef(name="__julia_value", typ=T_PYSSIZET, offset=fieldoffset(CPyJlWrapObject, 4), flags=READONLY, doc="Julia jl_value_t* (Any object)"),
        ],
        tp_dealloc = @cfunction(_pyjlwrap_dealloc, Cvoid, (PyPtr,)),
        tp_repr = @cfunction(_pyjlwrap_repr, PyPtr, (PyPtr,)),
        tp_str = @cfunction(_pyjlwrap_str, PyPtr, (PyPtr,)),
        tp_call = @cfunction(_pyjlwrap_call, PyPtr, (PyPtr,PyPtr,PyPtr)),
        tp_getattro = @cfunction(_pyjlwrap_getattr, PyPtr, (PyPtr,PyPtr)),
        tp_setattro = @cfunction(_pyjlwrap_setattr, Cint, (PyPtr,PyPtr,PyPtr)),
        tp_iter = @cfunction(_pyjlwrap_getiter, PyPtr, (PyPtr,)),
        tp_hash = sizeof(Py_hash_t) < sizeof(Int) ?
            @cfunction(_pyjlwrap_hash32, UInt32, (PyPtr,)) :
            @cfunction(_pyjlwrap_hash, UInt, (PyPtr,)),
        tp_weaklistoffset = fieldoffset(CPyJlWrapObject, 3),
        tp_richcompare = @cfunction(_pyjlwrap_richcompare, PyPtr, (PyPtr,PyPtr,Cint)),
        tp_as_number = CPyNumberMethods_default[],
        tp_as_mapping = CPyMappingMethods_default[],
        tp_methods = [
            CPyMethodDef(name="__dir__", meth=@cpymethod(_pyjlwrap_dir), flags=METH_NOARGS),
            CPyMethodDef(name="__complex__", meth=@cpymethod(_pyjlwrap_complex), flags=METH_NOARGS),
        ],
        tp_getset = [
            # CPyGetSetDef(name="__doc__", get=@cpygetfunc(_pyjlwrap_get_doc)),
            CPyGetSetDef(name="__name__", get=@cpygetfunc(_pyjlwrap_get_name_string)),
        ],
    )

    # ABSTRACT BASE CLASSES FROM `numbers`

    pyjlwrap_init_type!(CPyJlWrapNumber_Type,
        name = "JlWrapNumber",
        tp_base = CPyJlWrap_Type,
    )

    pyjlwrap_init_type!(CPyJlWrapComplex_Type,
        name = "JlWrapComplex",
        tp_base = CPyJlWrapNumber_Type,
        tp_getset = [
            CPyGetSetDef(name="real", get=@cpygetfunc(_pyjlwrap_get_real)),
            CPyGetSetDef(name="imag", get=@cpygetfunc(_pyjlwrap_get_imag)),
        ],
        tp_methods = [
            CPyMethodDef(name="conjugate", meth=@cpymethod(_pyjlwrap_conjugate), flags=METH_NOARGS),
        ]
    )

    pyjlwrap_init_type!(CPyJlWrapReal_Type,
        name = "JlWrapReal",
        tp_base = CPyJlWrapComplex_Type,
        tp_methods = [
            CPyMethodDef(name="trunc", meth=@cpymethod(_pyjlwrap_trunc), flags=METH_NOARGS),
            CPyMethodDef(name="round", meth=@cpymethod(_pyjlwrap_round), flags=METH_NOARGS),
            CPyMethodDef(name="floor", meth=@cpymethod(_pyjlwrap_floor), flags=METH_NOARGS),
            CPyMethodDef(name="ceil", meth=@cpymethod(_pyjlwrap_ceil), flags=METH_NOARGS),
        ]
    )

    pyjlwrap_init_type!(CPyJlWrapRational_Type,
        name = "JlWrapRational",
        tp_base = CPyJlWrapReal_Type,
        tp_getset = [
            CPyGetSetDef(name="numerator", get=@cpygetfunc(_pyjlwrap_get_numerator)),
            CPyGetSetDef(name="denominator", get=@cpygetfunc(_pyjlwrap_get_denominator)),
        ],
    )

    pyjlwrap_init_type!(CPyJlWrapIntegral_Type,
        name = "JlWrapIntegral",
        tp_base = CPyJlWrapRational_Type,
    )

    # ABSTRACT BASE CLASSES FROM `collections.abc`

    pyjlwrap_init_type!(CPyJlWrapIterable_Type,
        name = "JlWrapIterable",
        tp_base = CPyJlWrap_Type,
    )

    pyjlwrap_init_type!(CPyJlWrapIterator_Type,
        name = "JlWrapIterator",
        tp_base = CPyJlWrap_Type,
        tp_iter = @cfunction(pyincref_, PyPtr, (PyPtr,)),
        tp_iternext = @cfunction(_pyjlwrap_iternext, PyPtr, (PyPtr,)),
    )

    pyjlwrap_init_type!(CPyJlWrapCallable_Type,
        name = "JlWrapCallable",
        tp_base = CPyJlWrap_Type,
        tp_getset = [
            CPyGetSetDef(name="__qualname__", get=@cpygetfunc(_pyjlwrap_get_none)),
            CPyGetSetDef(name="__module__", get=@cpygetfunc(_pyjlwrap_get_none)),
            CPyGetSetDef(name="__closure__", get=@cpygetfunc(_pyjlwrap_get_none)),
        ],
    )

    pyjlwrap_init_type!(CPyJlWrapFunction_Type,
        name = "JlWrapFunction",
        tp_base = CPyJlWrapCallable_Type,
    )

    pyjlwrap_init_type!(CPyJlWrapContainer_Type,
        name = "JlWrapContainer",
        tp_base = CPyJlWrap_Type,
    )

    pyjlwrap_init_type!(CPyJlWrapCollection_Type,
        name = "JlWrapCollection",
        tp_base = CPyJlWrapIterable_Type,
    )

    pyjlwrap_init_type!(CPyJlWrapSequence_Type,
        name = "PyCall.JlWrapSequence",
        tp_base = CPyJlWrapCollection_Type,
        tp_as_sequence = CPySequenceMethods_oneup[],
        tp_as_mapping = CPyMappingMethods_oneup[],
    )

    pyjlwrap_init_type!(CPyJlWrapMutableSequence_Type,
        name = "JlWrapMutableSequence",
        tp_base = CPyJlWrapSequence_Type,
    )

    pyjlwrap_init_type!(CPyJlWrapByteString_Type,
        name = "JlWrapByteString",
        tp_base = CPyJlWrapSequence_Type,
    )

    pyjlwrap_init_type!(CPyJlWrapSet_Type,
        name = "JlWrapSet",
        tp_base = CPyJlWrapCollection_Type,
    )

    pyjlwrap_init_type!(CPyJlWrapMutableSet_Type,
        name = "JlWrapMutableSet",
        tp_base = CPyJlWrapSet_Type,
    )

    pyjlwrap_init_type!(CPyJlWrapMapping_Type,
        name = "JlWrapMapping",
        tp_base = CPyJlWrapCollection_Type,
        tp_as_mapping = CPyMappingMethods_default[],
        tp_iter = @cfunction(_pyjlwrap_getiter_keys, PyPtr, (PyPtr,)),
    )

    pyjlwrap_init_type!(CPyJlWrapJlNamedTuple_Type,
        name = "JlWrapJlNamedTuple",
        tp_base = CPyJlWrapMapping_Type,
        tp_as_sequence = CPySequenceMethods_oneup[],
        tp_as_mapping = CPyMappingMethods_namedtuple[],
    )

    pyjlwrap_init_type!(CPyJlWrapMutableMapping_Type,
        name = "JlWrapMutableMapping",
        tp_base = CPyJlWrapMapping_Type,
    )

    # ABSTRACT BASE CLASSES FROM `io`

    pyjlwrap_init_type!(CPyJlWrapIOBase_Type,
        name = "JlWrapIOBase",
        tp_base = CPyJlWrap_Type,
        tp_iter = @cfunction(_pyjlwrap_getiter_lines, PyPtr, (PyPtr,)),
        tp_methods = [
            CPyMethodDef(name="close", meth=@cpymethod(_pyjlwrap_io_close), flags=METH_NOARGS),
            CPyMethodDef(name="fileno", meth=@cpymethod(_pyjlwrap_io_fileno), flags=METH_NOARGS),
            CPyMethodDef(name="flush", meth=@cpymethod(_pyjlwrap_io_flush), flags=METH_NOARGS),
            CPyMethodDef(name="isatty", meth=@cpymethod(_pyjlwrap_io_isatty), flags=METH_NOARGS),
            CPyMethodDef(name="readable", meth=@cpymethod(_pyjlwrap_io_readable), flags=METH_NOARGS),
            CPyMethodDef(name="seek", meth=@cpymethod(_pyjlwrap_io_seek), flags=METH_VARARGS),
            CPyMethodDef(name="seekable", meth=@cpymethod(_pyjlwrap_io_seekable), flags=METH_NOARGS),
            CPyMethodDef(name="tell", meth=@cpymethod(_pyjlwrap_io_tell), flags=METH_NOARGS),
            CPyMethodDef(name="truncate", meth=@cpymethod(_pyjlwrap_io_truncate), flags=METH_VARARGS),
            CPyMethodDef(name="writable", meth=@cpymethod(_pyjlwrap_io_writable), flags=METH_NOARGS),
        ],
        tp_getset = [
            CPyGetSetDef(name="closed", get=@cpygetfunc(_pyjlwrap_io_get_closed)),
        ],
    )

    pyjlwrap_init_type!(CPyJlWrapRawIO_Type,
        name = "JlWrapRawIO",
        tp_base = CPyJlWrapIOBase_Type,
    )

    pyjlwrap_init_type!(CPyJlWrapBufferedIO_Type,
        name = "JlWrapBufferedIO",
        tp_base = CPyJlWrapIOBase_Type,
        tp_methods = [
            CPyMethodDef(name="read", meth=@cpymethod(_pyjlwrap_bufferedio_read), flags=METH_VARARGS),
            CPyMethodDef(name="readinto", meth=@cpymethod(_pyjlwrap_bufferedio_readinto), flags=METH_O),
            CPyMethodDef(name="write", meth=@cpymethod(_pyjlwrap_bufferedio_write), flags=METH_O),
        ],
    )

    pyjlwrap_init_type!(CPyJlWrapTextIO_Type,
        name = "JlWrapTextIO",
        tp_base = CPyJlWrapIOBase_Type,
        tp_methods = [
            CPyMethodDef(name="read", meth=@cpymethod(_pyjlwrap_textio_read), flags=METH_VARARGS),
            CPyMethodDef(name="readline", meth=@cpymethod(_pyjlwrap_textio_readline), flags=METH_VARARGS),
            CPyMethodDef(name="write", meth=@cpymethod(_pyjlwrap_textio_write), flags=METH_O),
        ],
        tp_getset = [
            CPyGetSetDef(name="encoding", get=@cpygetfunc(_pyjlwrap_textio_get_encoding_utf8))
        ],
    )

    # register with abstract base classes
    m = pyimport("numbers")
    m.Number.register(CPyJlWrapNumber_Type)
    m.Complex.register(CPyJlWrapComplex_Type)
    m.Real.register(CPyJlWrapReal_Type)
    m.Rational.register(CPyJlWrapRational_Type)
    m.Integral.register(CPyJlWrapIntegral_Type)

    m = pyimport("collections.abc")
    m.Iterable.register(CPyJlWrapIterable_Type)
    m.Iterator.register(CPyJlWrapIterator_Type)
    m.Callable.register(CPyJlWrapCallable_Type)
    m.Container.register(CPyJlWrapContainer_Type)
    m.Collection.register(CPyJlWrapCollection_Type)
    m.Sequence.register(CPyJlWrapSequence_Type)
    m.MutableSequence.register(CPyJlWrapMutableSequence_Type)
    m.ByteString.register(CPyJlWrapByteString_Type)
    m.Set.register(CPyJlWrapSet_Type)
    m.MutableSet.register(CPyJlWrapMutableSet_Type)
    m.Mapping.register(CPyJlWrapMapping_Type)
    m.MutableMapping.register(CPyJlWrapMutableMapping_Type)

    m = pyimport("io")
    m.IOBase.register(CPyJlWrapIOBase_Type)
    m.RawIOBase.register(CPyJlWrapRawIO_Type)
    m.TextIOBase.register(CPyJlWrapTextIO_Type)

end

##########################################################
# constructors for specific types

CPyJlWrap_From(x) = CPyJlWrap_New(CPyJlWrap_Type, x)
CPyJlWrap_From(x::Union{AbstractDict,AbstractArray,AbstractSet,NamedTuple,Tuple}) = CPyJlWrapIterable_From(x)
CPyJlWrap_From(x::Number) = CPyJlWrapNumber_From(x)
CPyJlWrap_From(x::IO) = CPyJlWrapIO_From(x)
CPyJlWrap_From(x::Union{Function,Type}) = CPyJlWrapCallable_From(x)

CPyJlWrapIterator_From(o) =
    let it = iterate(o)
        CPyJlWrap_New(CPyJlWrapIterator_Type, (o, Ref{Any}(it)))
    end

CPyJlWrapNumber_From(x) = CPyJlWrap_New(CPyJlWrapNumber_Type, x)
CPyJlWrapNumber_From(x::Complex) = CPyJlWrapComplex_From(x)
CPyJlWrapNumber_From(x::Real) = CPyJlWrapReal_From(x)

CPyJlWrapComplex_From(x) = CPyJlWrap_New(CPyJlWrapComplex_Type, x)
CPyJlWrapComplex_From(x::Real) = CPyJlWrapReal_From(x)

CPyJlWrapReal_From(x) = CPyJlWrap_New(CPyJlWrapReal_Type, x)
CPyJlWrapReal_From(x::Integer) = CPyJlWrapIntegral_From(x)
CPyJlWrapReal_From(x::Rational) = CPyJlWrapRational_From(x)

CPyJlWrapRational_From(x) = CPyJlWrap_New(CPyJlWrapRational_Type, x)
CPyJlWrapRational_From(x::Integer) = CPyJlWrapIntegral_From(x)

CPyJlWrapIntegral_From(x) = CPyJlWrap_New(CPyJlWrapIntegral_Type, x)

CPyJlWrapIterable_From(o) = CPyJlWrap_New(CPyJlWrapIterable_Type, o)
CPyJlWrapIterable_From(o::Union{AbstractDict,AbstractArray,AbstractSet,NamedTuple,Tuple}) = CPyJlWrapCollection_From(o)

CPyJlWrapCallable_From(o) = CPyJlWrap_New(CPyJlWrapCallable_Type, o)
CPyJlWrapCallable_From(o::Union{Function,Type}) = CPyJlWrapFunction_From(o)

CPyJlWrapFunction_From(o) = CPyJlWrap_New(CPyJlWrapFunction_Type, o)

CPyJlWrapCollection_From(o) = CPyJlWrap_New(CPyJlWrapCollection_Type, o)
CPyJlWrapCollection_From(o::Union{Tuple,AbstractArray}) = CPyJlWrapSequence_From(o)
CPyJlWrapCollection_From(o::AbstractSet) = CPyJlWrapSet_From(o)
CPyJlWrapCollection_From(o::Union{AbstractDict,NamedTuple}) = CPyJlWrapMapping_From(o)

CPyJlWrapSequence_From(o) = CPyJlWrap_New(CPyJlWrapSequence_Type, o)
CPyJlWrapSequence_From(o::AbstractArray) = CPyJlWrapMutableSequence_From(o)

CPyJlWrapMutableSequence_From(o) = CPyJlWrap_New(CPyJlWrapMutableSequence_Type, o)

CPyJlWrapMapping_From(o) = CPyJlWrap_New(CPyJlWrapMapping_Type, o)
CPyJlWrapMapping_From(o::Base.ImmutableDict) = CPyJlWrap_New(CPyJlWrapMapping_Type, o)
CPyJlWrapMapping_From(o::AbstractDict) = CPyJlWrapMutableMapping_From(o)
CPyJlWrapMapping_From(o::NamedTuple) = CPyJlWrapJlNamedTuple_From(o)

CPyJlWrapJlNamedTuple_From(o) = CPyJlWrap_New(CPyJlWrapJlNamedTuple_Type, o)

CPyJlWrapMutableMapping_From(o) = CPyJlWrap_New(CPyJlWrapMutableMapping_Type, o)

CPyJlWrapSet_From(o) = CPyJlWrap_New(CPyJlWrapSet_Type, o)
CPyJlWrapSet_From(o::AbstractSet) = CPyJlWrapMutableSet_From(o)

CPyJlWrapMutableSet_From(o) = CPyJlWrap_New(CPyJlWrapMutableSet_Type, o)

CPyJlWrapIO_From(o) = CPyJlWrap_New(CPyJlWrapIOBase_Type, o)
CPyJlWrapIO_From(o::IO) = CPyJlWrapBufferedIO_From(o)

CPyJlWrapBufferedIO_From(o) = CPyJlWrap_New(CPyJlWrapBufferedIO_Type, o)

CPyJlWrapRawIO_From(o) = CPyJlWrap_New(CPyJlWrapRawIO_Type, o)

CPyJlWrapTextIO_From(o) = CPyJlWrap_New(CPyJlWrapTextIO_Type, o)

export pyjlwrap, pyjlwrap_textio, pyjlwrap_rawio, pyjlwrap_bufferedio
pyjlwrap(x) = PyObject(CPyJlWrap_From(x))
pyjlwrap_rawio(x) = PyObject(CPyJlWrapRawIO_From(x))
pyjlwrap_textio(x) = PyObject(CPyJlWrapTextIO_From(x))
pyjlwrap_bufferedio(x) = PyObject(CPyJlWrapBufferedIO_From(x))

#########################################################################
# Precompilation: just an optimization to speed up initialization.
# Here, we precompile functions that are passed to cfunction by __init__,
# for the reasons described in JuliaLang/julia#12256.
precompile(_pyjlwrap_call, (PyPtr,PyPtr,PyPtr))
precompile(_pyjlwrap_dealloc, (PyPtr,))
precompile(_pyjlwrap_repr, (PyPtr,))
precompile(_pyjlwrap_hash, (PyPtr,))
precompile(_pyjlwrap_hash32, (PyPtr,))