text:
# Author: Mathias Louboutin
# Date: June 2021
#

using JUDI, LinearAlgebra, Images, PyPlot, DSP, ImageGather

# Set up model structure
n = (601, 333)   # (x,y,z) or (x,z)
d = (15., 15.)
o = (0., 0.)

# Velocity [km/s]
v = ones(Float32,n) .+ 0.5f0
for i=1:12
    v[:,25*i+1:end] .= 1.5f0 + i*.25f0
end
v0 = imfilter(v, Kernel.gaussian(5))
v0_low = .95f0 .* imfilter(v, Kernel.gaussian(45))
v0_high = 1.05f0 .* imfilter(v, Kernel.gaussian(45))

# Slowness squared [s^2/km^2]
m = (1f0 ./ v).^2
m0 = (1f0 ./ v0).^2
m0_low = (1f0 ./ v0_low).^2
m0_high = (1f0 ./ v0_high).^2
dm = vec(m - m0)

# Setup info and model structure
nsrc = 51	# number of sources
model = Model(n, d, o, m; nb=40)
model0 = Model(n, d, o, m0; nb=40)
model0_low = Model(n, d, o, m0_low; nb=40)
model0_high = Model(n, d, o, m0_high; nb=40)

# Set up receiver geometry
nxrec = 401
xrec = range(0f0, stop=(n[1] -1)*d[1], length=nxrec)
yrec = 0f0
zrec = range(20f0, stop=20f0, length=nxrec)

# receiver sampling and recording time
timeR = 4000f0   # receiver recording time [ms]
dtR = 4f0    # receiver sampling interval [ms]

# Set up receiver structure
recGeometry = Geometry(xrec, yrec, zrec; dt=dtR, t=timeR, nsrc=nsrc)
# Set up source geometry (cell array with source locations for each shot)
xsrc = convertToCell(range(0f0, stop=(n[1] -1)*d[1], length=nsrc))
ysrc = convertToCell(range(0f0, stop=0f0, length=nsrc))
zsrc = convertToCell(range(20f0, stop=20f0, length=nsrc))

# source sampling and number of time steps
timeS = 4000f0  # ms
dtS = 4f0   # ms

# Set up source structure
srcGeometry = Geometry(xsrc, ysrc, zsrc; dt=dtS, t=timeS)
# setup wavelet
f0 = 0.015f0     # kHz
wavelet = ricker_wavelet(timeS, dtS, f0)
q = judiVector(srcGeometry, wavelet)

###################################################################################################
opt = Options(space_order=16)
# Setup operators
F = judiModeling(model, srcGeometry, recGeometry; options=opt)
# Nonlinear modeling
dD = F*q

# Common surface offset image gather
offsets = 0f0:150f0:8000f0
CIG = surface_gather(model0, q, dD; offsets=offsets, options=opt)
CIG_low = surface_gather(model0_low, q, dD; offsets=offsets, options=opt);
CIG_high = surface_gather(model0_high, q, dD; offsets=offsets, options=opt);

cc = 1e-1
figure(figsize=(20, 10))
subplot(231)
imshow(CIG[301, :, :], vmin=-cc, vmax=cc, cmap="Greys", aspect="auto", extent=[0, 8, 5, 0])
xlabel("Surface offset (km)")
ylabel("depth (km)")
title("Good velocity (CDP=4.5km)")
subplot(232)
imshow(CIG_low[301, :, :], vmin=-cc, vmax=cc, cmap="Greys", aspect="auto", extent=[0, 8, 5, 0])
title("Low velocity (CDP=4.5km)")
xlabel("Surface offset (km)")
ylabel("depth (km)")
subplot(233)
imshow(CIG_high[301, :, :], vmin=-cc, vmax=cc, cmap="Greys", aspect="auto", extent=[0, 8, 5, 0])
title("High velocity (CDP=4.5km)")
xlabel("Surface offset (km)")
ylabel("depth (km)")

subplot(234)
imshow(CIG[101, :, :], vmin=-cc, vmax=cc, cmap="Greys", aspect="auto", extent=[0, 8, 5, 0])
xlabel("Surface offset (km)")
ylabel("depth (km)")
title("Good velocity (CDP=1.5km)")
subplot(235)
imshow(CIG_low[101, :, :], vmin=-cc, vmax=cc, cmap="Greys", aspect="auto", extent=[0, 8, 5, 0])
title("Low velocity (CDP=1.5km)")
xlabel("Surface offset (km)")
ylabel("depth (km)")
subplot(236)
imshow(CIG_high[101, :, :], vmin=-cc, vmax=cc, cmap="Greys", aspect="auto", extent=[0, 8, 5, 0])
title("High velocity (CDP=1.5km)")
xlabel("Surface offset (km)")
ylabel("depth (km)")
savefig("./docs/img/cig_cdp.png", bbox_inches="tight")

# Plot gathers as a pseudo rtm image
cig_rtm_good = hcat([CIG[i, :, 1:25] for i=1:20:n[1]]...)
cig_rtm_low = hcat([CIG_low[i, :, 1:25] for i=1:20:n[1]]...)
cig_rtm_high = hcat([CIG_high[i, :, 1:25] for i=1:20:n[1]]...)

figure(figsize=(20, 10))
subplot(131)
imshow(cig_rtm_good, vmin=-cc, vmax=cc, cmap="Greys", aspect="auto", extent=[0, 8, 5, 0])
xticks([])
ylabel("depth (km)")
title("Good velocity")
subplot(132)
imshow(cig_rtm_low, vmin=-cc, vmax=cc, cmap="Greys", aspect="auto", extent=[0, 8, 5, 0])
xticks([])
ylabel("depth (km)")
title("Low velocity")
subplot(133)
imshow(cig_rtm_high, vmin=-cc, vmax=cc, cmap="Greys", aspect="auto", extent=[0, 8, 5, 0])
xticks([])
ylabel("depth (km)")
title("High velocity")
savefig("./docs/img/cig_line.png", bbox_inches="tight")


text:
export NoParams

struct NoParams{D,BT<:AbstractBlock{D}} <: TagBlock{BT,D}
    content::BT
end

YaoBlocks.parameters!(out, x::NoParams) = out
YaoBlocks.nparameters(x::NoParams) = 0
YaoBlocks.niparams(x::NoParams) = 0
YaoBlocks.mat(::Type{T}, x::NoParams) where {T} = mat(T, content(x))
YaoBlocks.PropertyTrait(x::NoParams) = PreserveAll()
YaoBlocks.dispatch!(f::Union{Function,Nothing}, x::NoParams, it::YaoBlocks.Dispatcher) = x
YaoBlocks.generic_dispatch!(f::Union{Function,Nothing}, x::NoParams, it::YaoBlocks.Dispatcher) = x
Base.adjoint(x::NoParams) = NoParams(adjoint(content(x)))
YaoBlocks.chsubblocks(x::NoParams, blk::AbstractBlock) = NoParams(blk)

ext: jl


text:
using DiffEqFlux, GalacticOptim, OrdinaryDiffEq

#A desired MWE for now, not a test yet.

function rober(du,u,p,t)
  y₁,y₂,y₃ = u
  k₁,k₂,k₃ = p
  du[1] = -k₁*y₁ + k₃*y₂*y₃
  du[2] =  k₁*y₁ - k₃*y₂*y₃ - k₂*y₂^2
  du[3] =  y₁ + y₂ + y₃ - 1
  nothing
end
M = [1. 0  0
     0  1. 0
     0  0  0]
prob_mm = ODEProblem(ODEFunction(rober,mass_matrix=M),[1.0,0.0,0.0],(0.0,10.0),(0.04,3e7,1e4))
sol = solve(prob_mm,Rodas5(),reltol=1e-8,abstol=1e-8)

dudt2 = Chain(x -> x.^3,Dense(6,50,tanh),Dense(50,2))

ndae = NeuralDAE(dudt2, (u,p,t) -> [u[1] + u[2] + u[3] - 1], tspan, M, DImplicitEuler(),
                        differential_vars = [true,true,false])
truedu0 = similar(u₀)
f(truedu0,u₀,p,0.0)

ndae(u₀,truedu0,Float64.(ndae.p))

function predict_n_dae(p)
    ndae(u₀,p)
end

function loss(p)
    pred = predict_n_dae(p)
    loss = sum(abs2,sol .- pred)
    loss,pred
end

p = p .+ rand(3) .* p

optfunc = GalacticOptim.OptimizationFunction((x, p) -> loss(x), GalacticOptim.AutoZygote())
optprob = GalacticOptim.OptimizationProblem(optfunc, p)
res = GalacticOptim.solve(optprob, BFGS(initial_stepnorm = 0.0001))


text:
using ModelTES, PyPlot

#following equation 4
# https://doi.org/10.1063/1.4984065 Morgan et al 2017
type TwoFluidRIT <: ModelTES.AbstractRIT
  Ic0::Float64
  cI::Float64
  cR::Float64

end
function ModelTES.R(I,T, RIT::TwoFluidRIT, Tc, Rn)
  a=(1-min(T,Tc)/Tc)^1.5
  b=RIT.cI*RIT.Ic0/I
  c=(1-a*b)
  Rn*RIT.cR*c
end
ModelTES.transitionwidth(RIT::TwoFluidRIT) = 0.001 # I dont think this is used much
tfrit = TwoFluidRIT(0.00742, 0.84, 0.45)
Tbath = 0.055; L = 45e-9; Rl = 0.0; C = 0.25e-12; G=42e-12; n=2.5;
Rn = 0.010; Tc = 0.0714; k=G/n/Tc^(n-1)
R0 = 0.15*Rn; Rpara = 0.0
tes_param = TESParams(n,Tc,Tbath,k,C,L,Rl,Rpara,Rn,tfrit)
bt = BiasedTES(tes_param, R0)

out = pulse(12000,1e-7, bt, 1000, 2000);

plot(times(out), out.I[1]-out.I)


text:
using Luxor

function sierpinski(txy, levelsyet)
    nxy = zeros(6)
    if levelsyet > 0
        for i in 1:6
            pos = i < 5 ? i + 2 : i - 4
            nxy[i] = (txy[i] + txy[pos]) / 2.0
        end
        sierpinski([txy[1],txy[2],nxy[1],nxy[2],nxy[5],nxy[6]], levelsyet-1)
        sierpinski([nxy[1],nxy[2],txy[3],txy[4],nxy[3],nxy[4]], levelsyet-1)
        sierpinski([nxy[5],nxy[6],nxy[3],nxy[4],txy[5],txy[6]], levelsyet-1)
    else
        poly([Point(txy[1],txy[2]),Point(txy[3],txy[4]),Point(txy[5],txy[6])], :fill ,close=true)
    end
end

Drawing(800, 800)
sierpinski([400., 100., 700., 500., 100., 500.], 7)
finish()
preview()

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
#-------------------------------------------------------------------
#* EMSO Model Library (EML) Copyright (C) 2004 - 2007 ALSOC.
#*
#* This LIBRARY is free software; you can distribute it and/or modify
#* it under the therms of the ALSOC FREE LICENSE as available at
#* http://www.enq.ufrgs.br/alsoc.
#*
#* EMSO Copyright (C) 2004 - 2007 ALSOC, original code
#* from http://www.rps.eng.br Copyright (C) 2002-2004.
#* All rights reserved.
#*
#* EMSO is distributed under the therms of the ALSOC LICENSE as
#* available at http://www.enq.ufrgs.br/alsoc.
#*--------------------------------------------------------------------
#* Author: Gerson Balbueno Bicca
#* $Id$
#*--------------------------------------------------------------------
type LMTD_Basic
	LMTD_Basic()=begin
		new(
			temp_delta(Dict{Symbol,Any}(
				:Brief=>"Temperature Difference at Inlet",
				:Lower=>1e-6,
				:Symbol=>"\\Delta T_0",
				:Protected=>true
			)),
			temp_delta(Dict{Symbol,Any}(
				:Brief=>"Temperature Difference at Outlet",
				:Lower=>1e-6,
				:Symbol=>"\\Delta T_L",
				:Protected=>true
			)),
			temp_delta(Dict{Symbol,Any}(
				:Brief=>"Logarithmic Mean Temperature Difference",
				:Lower=>1e-6,
				:Protected=>true
			)),
			positive(Dict{Symbol,Any}(
				:Brief=>"LMTD Correction Factor",
				:Lower=>0.1,
				:Protected=>true
			)),
			[
				:(LMTD = 0.5*(DT0+DTL)),
				:(LMTD= (DT0-DTL)/ln(DT0/DTL)),
				:(LMTD = 0.5*(DT0+DTL)*(1-(DT0-DTL)^2/(DT0*DTL)*(1+(DT0-DTL)^2/(DT0*DTL)/2)/12)),
			],
			[
				"Log Mean Temperature Difference","Log Mean Temperature Difference","Log Mean Temperature Difference",
			],
			[:DT0,:DTL,:LMTD,:Fc,]
		)
	end
	DT0::temp_delta
	DTL::temp_delta
	LMTD::temp_delta
	Fc::positive
	equations::Array{Expr,1}
	equationNames::Array{String,1}
	variables::Array{Symbol,1}
	attributes::Dict{Symbol,Any}
end
export LMTD_Basic
function setEquationFlow(in::LMTD_Basic)
	if DT0*DTL <= 0 * "K^2"
		addEquation(1)
	else
		if abs(DT0 - DTL) > 0.05*max(abs([DT0,DTL]))
			addEquation(2)
		else
			addEquation(3)
		end
	end
end
function atributes(in::LMTD_Basic,_::Dict{Symbol,Any})
	fields::Dict{Symbol,Any}=Dict{Symbol,Any}()
	fields[:Pallete]=false
	fields[:Brief]="Log Mean Temperature Difference Method."
	fields[:Info]="This model should be used as submodel when the LMTD needs to be calculating"
	drive!(fields,_)
	return fields
end
LMTD_Basic(_::Dict{Symbol,Any})=begin
	newModel=LMTD_Basic()
	newModel.attributes=atributes(newModel,_)
	newModel
end


text:
using FFTW
using MAT
using Random
using Printf
using StaticArrays

include("models.jl")

"""
Generate magnetic susceptibility maps and the corresponding fields due to
randomly placed spheres and randomly oriented infinite cylinders subject to a
static B-field inside a 3d rectangular volume.

Method
------
    1. Pick array size, voxel size for grid (N/n grid)
    2. for n iterations
        i.   pick number of tries, ntry, to check for valid config
        ii.  randomly place sources inside grid until no valid config can be
             found ntry times in a row
        iii. valid config means no overlapping structures. concentric allowed.
"""
function generate(
    sz::NTuple{N, Integer},
    vsz::NTuple{3, Real};
    outpath::Union{Nothing, AbstractString} = nothing,
    nsources::Integer = typemax(Int)-1,
    density::Real = 0.51,
    bdir::NTuple{3, Real} = (0, 0, 1),
    χμ::Float64 = 1e-6,
    χσ::Float64 = 0.07,
    seed::Int = 12345,
) where {N}
    @assert nsources > 0
    @assert density > 0 && density < 1
    Random.seed!(seed)

    if !(norm(bdir) ≈ 1)
        @warn "B direction vector is not normalized. Normalizing..."
        bdir = bdir ./ norm(bdir)
    end

    x = zeros(sz)
    f = zeros(sz)
    m = Array{Bool}(falses(sz))

    rmax = cld(minimum(sz), 5)
    rmax = 7
    Cs = ntuple(i -> rmax+1:sz[i]-(rmax+1)-1, Val(3))
    Cs = ntuple(i -> 10*rmax+1:sz[i]-(10*rmax+1)-1, Val(3))

    ns = 0
    while true
        model = rand((:sphere, :cylinder))
        model = :sphere

        r = rand(3:rmax)
        c = rand.(Cs)
        χ = χμ + χσ*randn()

        if model == :sphere
            x, f, m = sphere!(x, f, m, c, r, χ, bdir, vsz)
            ns += 1

        elseif model == :cylinder
            θ = π*rand()
            ϕ = 2π*rand()

            x, f, m = cylinder!(x, f, m, c, r, θ, ϕ, χ, bdir, vsz)
            ns += 1
        end

        if sum(m)/length(m) > density || ns >= nsources
            break
        end
    end

    if outpath !== nothing
        matwrite(outpath, Dict(:x => x, :f => f, :m => m), compress=true)
    end

    return x, f, m
end


text:
#############################################################
# Data Layers
#############################################################
include("layers/hdf5-data.jl")

#############################################################
# General Computation Layers
#############################################################
include("layers/inner-product.jl")
include("layers/convolution.jl")
include("layers/pooling.jl")
include("layers/softmax.jl")
include("layers/power.jl")
include("layers/split.jl")
include("layers/element-wise.jl")
include("layers/channel-pooling.jl")
include("layers/dropout.jl")
include("layers/argmax.jl")
include("layers/crop.jl")
include("layers/concat.jl")
include("layers/tied-inner-product.jl")
include("layers/index2onehot.jl")
include("layers/random-normal.jl")

#############################################################
# Loss Layers
#############################################################
include("layers/square-loss.jl")
include("layers/hinge-loss.jl")
include("layers/multinomial-logistic-loss.jl")
include("layers/softmax-loss.jl")
include("layers/softlabel-softmax-loss.jl")
include("layers/wasserstein-loss.jl")
include("layers/binary-cross-entropy-loss.jl")
include("layers/gaussian-kl-loss.jl")

#############################################################
# Statistics Layers
#############################################################
include("layers/accuracy.jl")
include("layers/binary-accuracy.jl")


text:
"""
    Create a mesh-like initial configuration for the simulation system.
It uses the information from the system such as the cut-off radius and the
density of the system to allocate the particles.

It modifies the original `positions` array, so it is an in-place operation.
"""
function init!(positions::AbstractArray, syst::System)
    dist = cbrt(1.0f0 / syst.ρ)
    dist_half = dist / 2.0f0

    # Define the first positions
    positions[1, 1] = -syst.rc + dist_half
    positions[2, 1] = -syst.rc + dist_half
    positions[3, 1] = -syst.rc + dist_half

    # Create a complete lattice
    @inbounds for i = 2:(syst.N - 1)
        positions[1, i] = positions[1, i - 1] + dist
        positions[2, i] = positions[2, i - 1]
        positions[3, i] = positions[3, i - 1]

        if positions[1, i] > syst.rc
            positions[1, i] = positions[1, 1]
            positions[2, i] = positions[2, i - 1] + dist

            if positions[2, i] > syst.rc
                positions[1, i] = positions[1, 1]
                positions[2, i] = positions[2, 1]
                positions[3, i] = positions[3, i - 1] + dist
            end
        end
    end

    return nothing
end

function snapshot(positions, filename, i)
    (N, M) = size(positions)
    snapshot_file = joinpath(filename, "system_$(i).csv")
    CSV.write(snapshot_file, DataFrame(reshape(positions, (M, N))))

    return nothing
end


text:
@testset "interpolations" begin

for interpolation in (Lagrange{1, RefCube, 1}(),
                      Lagrange{1, RefCube, 2}(),
                      Lagrange{2, RefCube, 1}(),
                      Lagrange{2, RefCube, 2}(),
                      Lagrange{2, RefTetrahedron, 1}(),
                      Lagrange{2, RefTetrahedron, 2}(),
                      Lagrange{3, RefCube, 1}(),
                      Serendipity{2, RefCube, 2}(),
                      Lagrange{3, RefTetrahedron, 1}(),
                      Lagrange{3, RefTetrahedron, 2}())

    # Test of utility functions
    ndim = JuAFEM.getdim(interpolation)
    r_shape = JuAFEM.getrefshape(interpolation)
    func_order = JuAFEM.getorder(interpolation)
    @test typeof(interpolation) <: Interpolation{ndim,r_shape,func_order}
    @test typeof(JuAFEM.getlowerdim(interpolation)) <: Interpolation{ndim-1}
    @test typeof(JuAFEM.getlowerorder(interpolation)) <: Interpolation{ndim,r_shape,func_order-1}

    n_basefuncs = getnbasefunctions(interpolation)
    x = rand(Tensor{1, ndim})
    f = (x) -> JuAFEM.value(interpolation, Tensor{1, ndim}(x))
    @test vec(ForwardDiff.jacobian(f, Array(x))') ≈
           reinterpret(Float64, JuAFEM.derivative(interpolation, x))
    @test sum(JuAFEM.value(interpolation, x)) ≈ 1.0

    coords = JuAFEM.reference_coordinates(interpolation)
    for node in 1:n_basefuncs
        N_node = JuAFEM.value(interpolation, coords[node])
        for k in 1:node
            if k == node
                @test N_node[k] ≈ 1.0
            else
                @test N_node[k] ≈ 0.0
            end
        end
    end
end

end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "MongoC"
version = v"1.19.1"

# Collection of sources required to complete build
sources = [
    ArchiveSource("https://github.com/mongodb/mongo-c-driver/releases/download/$(version)/mongo-c-driver-$(version).tar.gz",
                  "1732251e3f65bc02ce05c04ce34ef2819b154479108df669f0c045486952521d")
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir
cd mongo-c-driver-*
sed -i "s/Windows.h/windows.h/" src/libmongoc/src/mongoc/mongoc-client.c
sed -i "s/WinDNS.h/windns.h/" src/libmongoc/src/mongoc/mongoc-client.c
sed -i "s/Mstcpip.h/mstcpip.h/" src/libmongoc/src/mongoc/mongoc-client.c
sed -i "s/Mstcpip.h/mstcpip.h/" src/libmongoc/src/mongoc/mongoc-socket.c
sed -i "s/Dnsapi/dnsapi/" build/cmake/FindResSearch.cmake
mkdir cmake-build
cd cmake-build

if [[ "${nbits}" == 32 ]]; then
    export CFLAGS="-Wl,-rpath-link,/opt/${target}/${target}/lib"
elif [[ "${target}" != *-apple-* ]]; then
    export CFLAGS="-Wl,-rpath-link,/opt/${target}/${target}/lib64"
fi
if [[ "${target}" == *-mingw* ]]; then
    export CFLAGS="$CFLAGS -D__USE_MINGW_ANSI_STDIO=1"
fi

cmake \
    -DENABLE_AUTOMATIC_INIT_AND_CLEANUP=OFF \
    -DCMAKE_BUILD_TYPE=Release \
    -DENABLE_SSL=OPENSSL \
    -DENABLE_SASL=OFF \
    -DENABLE_EXAMPLES=OFF \
    -DENABLE_TESTS=OFF \
    -DENABLE_UNINSTALL=OFF \
    -DENABLE_STATIC=OFF \
    -DENABLE_SNAPPY=ON \
    -DENABLE_MONGOC=ON \
    -DENABLE_BSON=ON \
    -DOPENSSL_ROOT_DIR=$prefix \
    -DCMAKE_INSTALL_PREFIX=$prefix \
    -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} \
    ..
make -j${nproc}
make install

if [[ "${target}" == *-apple-darwin* ]]; then
    rm ${libdir}/libbson-1.0.0.dylib
    rm ${libdir}/libmongoc-1.0.0.dylib
    mv ${libdir}/libmongoc-1.0.0.0.0.dylib ${libdir}/libmongoc-1.0.0.dylib
    mv ${libdir}/libbson-1.0.0.0.0.dylib ${libdir}/libbson-1.0.0.dylib
fi
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms(; experimental=true)

# The products that we will ensure are always built
products = [
    LibraryProduct(["libbson-1", "libbson-1.0", "libbson"], :libbson),
    LibraryProduct(["libmongoc-1", "libmongoc-1.0", "libmongoc"], :libmongoc),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency(PackageSpec(name="OpenSSL_jll", uuid="458c3c95-2e84-50aa-8efc-19380b2a3a95"))
    Dependency(PackageSpec(name="Zlib_jll", uuid="83775a58-1f1d-513f-b197-d71354ab007a"))
    Dependency(PackageSpec(name="Zstd_jll", uuid="3161d3a3-bdf6-5164-811a-617609db77b4"))
    Dependency(PackageSpec(name="snappy_jll", uuid="fe1e1685-f7be-5f59-ac9f-4ca204017dfd"), v"1.1.9")
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; julia_compat = "1.6")

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
struct tSTE <: AbstractLoss
    α::Int
    constant::Float64

    function tSTE(;α::Int = 3)
        α ≥ 1 || throw(ArgumentError("α in tSTE loss must be ≥ 1"))
        new(α, (α + 1) / α)
    end

end

function kernel(loss::tSTE, X::Embedding)
    K = pairwise(SqEuclidean(), X.X, dims=2)
    Q = zeros(Float64, nitems(X), nitems(X))

    constant = ((loss.α + 1) / -2)

    for ij in eachindex(K)
        # We compute K[i,j] = (1 + ||x_i - x_j||_2^2/α) ^ (-(α+1)/2)
        base = (1 + K[ij] / loss.α)
        @inbounds Q[ij] = base ^ -1
        @inbounds K[ij] = base ^ constant
    end
    return K, Q
end

function gradient(loss::tSTE, triplets, X::AbstractMatrix)

    K, Q = kernel(loss, X) # Triplet kernel values (in the tSTE loss)

    # We need to create an array to prevent race conditions
    # This is the best average solution for small and big Embeddings
    nthreads = Threads.nthreads()
    triplets_range = partition(length(triplets), nthreads)

    C = zeros(Float64, nthreads)
    ∇C = [zeros(Float64, size(X)) for _ in 1:nthreads]

    Threads.@threads for tid in 1:nthreads
        C[tid] = tgradient!(∇C[tid], loss, triplets, X, K, Q, triplets_range[tid])
    end

    return sum(C), -sum(∇C)
end

function tgradient!(
    ∇C::Matrix{<:AbstractFloat},
    loss::tSTE,
    triplets::Triplets,
    X::AbstractMatrix,
    K::Matrix{<:AbstractFloat},
    Q::Matrix{<:AbstractFloat},
    triplets_range::UnitRange{Int64})

    C = 0.0

    for t in triplets_range
        @views @inbounds i, j, k = triplets[t][:i], triplets[t][:j], triplets[t][:k]

        @inbounds P = K[i,j] / (K[i,j] + K[i,k])
        C += -log(P)

        for d in 1:ndims(X)
            @inbounds ∂x_j = (1 - P) * Q[i,j] * (X[d,i] - X[d,j])
            @inbounds ∂x_k = (1 - P) * Q[i,k] * (X[d,i] - X[d,k])

            @inbounds ∇C[d,i] +=   loss.constant * (∂x_k - ∂x_j)
            @inbounds ∇C[d,j] +=   loss.constant *  ∂x_j
            @inbounds ∇C[d,k] += - loss.constant *  ∂x_k
        end
    end
    return C
end

function tgradient!(
    ∇C::Matrix{<:AbstractFloat},
    loss::tSTE,
    triplets::LabeledTriplets,
    X::AbstractMatrix,
    K::Matrix{<:AbstractFloat},
    Q::Matrix{<:AbstractFloat},
    triplets_range::UnitRange{Int64})

    C = 0.0

    for t in triplets_range
        @views @inbounds i, j, k = if triplets[t][:y] == -1
            triplets[t][:i], triplets[t][:j], triplets[t][:k]
        else
            triplets[t][:i], triplets[t][:k], triplets[t][:j]
        end

        @inbounds P = K[i,j] / (K[i,j] + K[i,k])
        C += -log(P)

        for d in 1:ndims(X)
            @inbounds ∂x_j = (1 - P) * Q[i,j] * (X[d,i] - X[d,j])
            @inbounds ∂x_k = (1 - P) * Q[i,k] * (X[d,i] - X[d,k])

            @inbounds ∇C[d,i] +=   loss.constant * (∂x_k - ∂x_j)
            @inbounds ∇C[d,j] +=   loss.constant *  ∂x_j
            @inbounds ∇C[d,k] += - loss.constant *  ∂x_k
        end
    end
    return C
end


text:
import Pkg, PackageCompiler
Pkg.activate("TowardsNotebook")
Pkg.instantiate()
PackageCompiler.create_app("TowardsNotebook", "TowardsNotebookCompiled")

max_stars_repo_name: Jblew/towards-julia-notebook

max_issues_repo_name: Jblew/towards-julia-notebook

max_forks_repo_name: Jblew/towards-julia-notebook


text:
module sky130_fd_pr
	# ERROR: TODO: Maybe have a julia-side SPICE model registry also,
	# to avoid requiring explicit .lib in the `.sp` files?
	__init__() = error("This is a SPICE model package, not a julia package. Don't attempt to load it like this.")
end


text:
# Formulas for representing and working with linear-model-type expressions
# Original by Harlan D. Harris.  Later modifications by John Myles White,
# Douglas M. Bates, and other contributors.

## Formulas are written as expressions and parsed by the Julia parser.
## For example :(y ~ a + b + log(c))
## In Julia the & operator is used for an interaction.  What would be written
## in R as y ~ a + b + a:b is written :(y ~ a + b + a&b) in Julia.
## The equivalent R expression, y ~ a*b, is the same in Julia

## The lhs of a one-sided formula is 'nothing'
## The rhs of a formula can be 1

is_call(ex::Expr) = Meta.isexpr(ex, :call)
is_call(ex::Expr, op::Symbol) = Meta.isexpr(ex, :call) && ex.args[1] == op
is_call(::Any) = false
is_call(::Any, ::Any) = false
check_call(ex) = is_call(ex) || throw(ArgumentError("non-call expression encountered: $ex"))

catch_dollar(ex::Expr) =
    Meta.isexpr(ex, :$) && throw(ArgumentError("interpolation with \$ not supported in @formula.  Use @eval @formula(...) instead."))

"""
    @formula(ex)

Capture and parse a formula expression as a `Formula` struct.

A formula is an abstract specification of a dependence between _left-hand_ and
_right-hand_ side variables as in, e.g., a regression model.  Each side
specifies at a high level how tabular data is to be converted to a numerical
matrix suitable for modeling.  This specification looks something like Julia
code, is represented as a Julia `Expr`, but uses special syntax.  The `@formula`
macro takes an expression like `y ~ 1 + a*b`, transforms it according to the
formula syntax rules into a lowered form (like `y ~ 1 + a + b + a&b`), and
constructs a `Formula` struct which captures the original expression, the
lowered expression, and the left- and right-hand-side.

Operators that have special interpretations in this syntax are

* `~` is the formula separator, where it is a binary operator (the first
  argument is the left-hand side, and the second is the right-hand side.
* `+` concatenates variables as columns when generating a model matrix.
* `&` representes an _interaction_ between two or more variables, which
  corresponds to a row-wise kronecker product of the individual terms
  (or element-wise product if all terms involved are continuous/scalar).
* `*` expands to all main effects and interactions: `a*b` is equivalent to
  `a+b+a&b`, `a*b*c` to `a+b+c+a&b+a&c+b&c+a&b&c`, etc.
* `1`, `0`, and `-1` indicate the presence (for `1`) or absence (for `0` and
  `-1`) of an intercept column.

The rules that are applied are

* The associative rule (un-nests nested calls to `+`, `&`, and `*`).
* The distributive rule (interactions `&` distribute over concatenation `+`).
* The `*` rule expands `a*b` to `a+b+a&b` (recursively).
* Subtraction is converted to addition and negation, so `x-1` becomes `x + -1`
  (applies only to subtraction of literal 1).
* Single-argument `&` calls are stripped, so `&(x)` becomes the main effect `x`.
"""
macro formula(ex)
    is_call(ex, :~) || throw(ArgumentError("expected formula separator ~, got $(ex.head)"))
    length(ex.args) == 3 ||  throw(ArgumentError("malformed expression in formula $ex"))
    terms!(sort_terms!(parse!(ex)))
end

"""
    abstract type FormulaRewrite end

Formula parsing is expressed as a bunch of expression re-writes, each of which
is a subtype of `FormulaRewrite`.  There are two methods that dispatch on these
types: `applies(ex, child_idx, rule::Type{<:FormulaRewrite})` checks whether the
re-write `rule` needs to be applied at argument `child_idx` of expression `ex`,
and `rewrite!(ex, child_idx, rule::Type{<:FormulaRewrite})` re-writes `ex`
according to `rule` at position `child_idx`, and returns the next `child_idx`
that needs to be checked.
"""
abstract type FormulaRewrite end

"""
    struct Star <: FormulaRewrite end

Expand `a*b` to `a + b + a&b` (`*(a,b)` to `+(a,b,&(a,b))`).  Applies
recursively to multiple `*` arguments, so needs a clean-up pass (from
distributive/associative).
"""
struct Star <: FormulaRewrite end
applies(ex::Expr, child_idx::Int, ::Type{Star}) = is_call(ex.args[child_idx], :*)
expand_star(a, b) = Expr(:call, :+, a, b, Expr(:call, :&, a, b))
function rewrite!(ex::Expr, child_idx::Int, ::Type{Star})
    child = ex.args[child_idx]
    @debug "  expand star: $ex -> "
    child.args = reduce(expand_star, child.args[2:end]).args
    @debug "               $ex"
    child_idx
end

"""
    struct AssociativeRule <: FormulaRewrite end

Apply associative rule: if in an expression headed by an associative operator
(`+,&,*`) and the sub-expression `child_idx` is headed by the same operator,
splice that child's children into it's location.
"""
struct AssociativeRule <: FormulaRewrite end
const ASSOCIATIVE = (:+, :&, :*)
applies(ex::Expr, child_idx::Int, ::Type{AssociativeRule}) =
    is_call(ex) &&
    ex.args[1] in ASSOCIATIVE &&
    is_call(ex.args[child_idx], ex.args[1])
function rewrite!(ex::Expr, child_idx::Int, ::Type{AssociativeRule})
    @debug "    associative: $ex -> "
    splice!(ex.args, child_idx, ex.args[child_idx].args[2:end])
    @debug "                 $ex"
    child_idx
end

"""
    struct Distributive <: FormulaRewrite end

Distributive propery: `&(a..., +(b...), c...)` to `+(&(a..., b_i, c...)_i...)`.
Replace outer call (:&) with inner call (:+), whose arguments are copies of the
outer call, one for each argument of the inner call.  For the ith new child, the
original inner call is replaced with the ith argument of the inner call.
"""
struct Distributive <: FormulaRewrite end
const DISTRIBUTIVE = Set([:& => :+])
applies(ex::Expr, child_idx::Int, ::Type{Distributive}) =
    is_call(ex) &&
    is_call(ex.args[child_idx]) &&
    (ex.args[1] => ex.args[child_idx].args[1]) in DISTRIBUTIVE
function rewrite!(ex::Expr, child_idx::Int, ::Type{Distributive})
    @debug "    distributive: $ex -> "
    new_args = deepcopy(ex.args[child_idx].args)
    for i in 2:length(new_args)
        new_child = deepcopy(ex)
        new_child.args[child_idx] = new_args[i]
        new_args[i] = new_child
    end
    ex.args = new_args
    @debug "                  $ex"
    # TODO: is it really necessary to re-check _every_argument after this?
    2
end

"""
    And1 <: FormulaRewrite

Remove numbers from interaction terms, so `1&x` becomes `&(x)` (which is later
cleaned up by `EmptyAnd`).
"""
struct And1 <: FormulaRewrite end
applies(ex::Expr, child_idx::Int, ::Type{And1}) =
    is_call(ex, :&) && ex.args[child_idx] isa Number
function rewrite!(ex::Expr, child_idx::Int, ::Type{And1})
    @debug "    &1: $ex ->"
    ex.args[child_idx] == 1 ||
        @warn "Number $(ex.args[child_idx]) removed from interaction term $ex"
    deleteat!(ex.args, child_idx)
    @debug "        $ex"
    child_idx
end

# default re-write is a no-op (go to next child)
rewrite!(ex::Expr, child_idx::Int, ::Nothing) = child_idx+1

# like `findfirst` but returns the first element where predicate is true, or
# nothing
function filterfirst(f::Function, a::AbstractArray)
    idx = findfirst(f, a)
    idx === nothing ? nothing : a[idx]
end

const SPECIALS = (:+, :&, :*, :~)

parse!(x) = parse!(x, [And1, Star, AssociativeRule, Distributive])
parse!(x, rewrites) = x
function parse!(ex::Expr, rewrites::Vector)
    @debug "parsing $ex"
    catch_dollar(ex)
    check_call(ex)

    # don't recurse into captured calls
    if is_call(ex, :capture_call) || is_call(ex, :(StatsModels.capture_call))
        @debug "  skipping capture_call"
        return ex
    end

    # parse a copy of non-special calls
    ex_parsed = ex.args[1] ∉ SPECIALS ? deepcopy(ex) : ex

    # iterate over children, checking for special rules
    child_idx = 2
    while child_idx <= length(ex_parsed.args)
        @debug "  ($(ex_parsed.args[1])) i=$child_idx: $(ex_parsed.args[child_idx])"
        # depth first: parse each child first
        parse!(ex_parsed.args[child_idx], rewrites)
        # find first rewrite rule that applies
        rule = filterfirst(r->applies(ex_parsed, child_idx, r), rewrites)
        # re-write according to that rule and update the child to position rewrite nex_parsedt
        child_idx = rewrite!(ex_parsed, child_idx, rule)
    end
    @debug "done: $ex_parsed"

    if ex.args[1] ∈ SPECIALS
        return ex_parsed
    else
        @debug "  capturing non-special call $ex"
        return capture_call_ex!(ex, ex_parsed)
    end
end

"""
    capture_call_ex!(ex::Expr, ex_parsed::Expr)

Capture a call to a function that is not part of the formula DSL.  This replaces
`ex` with a call to [`capture_call`](@ref).  `ex_parsed` is a copy of `ex` whose
arguments have been parsed according to the normal formula DSL rules and which
will be passed as the final argument to `capture_call`.
"""
function capture_call_ex!(ex::Expr, ex_parsed::Expr)
    symbols = extract_symbols(ex)
    symbols_ex = Expr(:tuple, symbols...)
    f_anon_ex = esc(Expr(:(->), symbols_ex, copy(ex)))
    f_orig = ex.args[1]
    ex.args = [:capture_call,
               esc(f_orig),
               f_anon_ex,
               tuple(symbols...),
               Meta.quot(deepcopy(ex)),
               :[$(ex_parsed.args[2:end]...)]]
    return ex
end

# generate Term expressions for symbols (including parsed args of non-special calls
terms!(::Nothing) = :(nothing)
terms!(s::Symbol) = :(Term($(Meta.quot(s))))
terms!(n::Number) = :(ConstantTerm($n))
function terms!(ex::Expr)
    if ex.args[1] ∈ SPECIALS
        ex.args[1] = esc(ex.args[1])
        ex.args[2:end] .= terms!.(ex.args[2:end])
    elseif is_call(ex, :capture_call)
        # final argument of capture_call holds parsed terms
        ex.args[end].args .= terms!.(ex.args[end].args)
    end
    return ex
end

function sort_terms!(ex::Expr)
    check_call(ex)
    if ex.args[1] ∈ ASSOCIATIVE
        sort!(view(ex.args, 2:length(ex.args)), by=degree)
    elseif is_call(ex, :capture_call)
        sort_terms!.(ex.args[end].args)
    else
        # recursively sort children
        sort_terms!.(ex.args)
    end
    ex
end
sort_terms!(x) = x

degree(i::Integer) = 0
degree(::Symbol) = 1
function degree(ex::Expr)
    check_call(ex)
    if ex.args[1] == :&
        length(ex.args) - 1
    elseif ex.args[1] == :|
        # put ranef terms at end
        typemax(Int)
    else
        # arbitrary functions are treated as main effect terms
        1
    end
end


text:
include("common.jl")

@testset "Optimize single objective problems in 5, 10, 30 and 100 dimensions with sNES" begin

@testset "$problem" for problem in ["Sphere", "Schwefel2.22", "Schwefel2.22"]
    p = BlackBoxOptim.Problems.examples[problem]

    @test fitness_for_opt(p, 5, 20,  1e2, separable_nes) < 0.01
    @test fitness_for_opt(p, 5, 20,  5e2, xnes) < 0.01

    @test fitness_for_opt(p, 10, 20, 3e2, separable_nes) < 0.01
    @test fitness_for_opt(p, 10, 20, 1e3, xnes) < 0.01

    @test fitness_for_opt(p, 30, 25, 1e3, separable_nes) < 0.01
    @test fitness_for_opt(p, 30, 25, 4e3, xnes) < 0.01

    @test fitness_for_opt(p, 100, 25, 5e3, separable_nes) < 0.01
    # Cannot run xnes in 100 dimensions since it scales badly.
end

@testset "Schwefel1.2" begin
    p = BlackBoxOptim.Problems.examples["Schwefel1.2"]

    @test fitness_for_opt(p, 30, 50, 4e3, separable_nes) < 10.0
    @test fitness_for_opt(p, 30, 50, 8e3, xnes) < 10.0
end

@testset "Rosenbrock" begin
    p = BlackBoxOptim.Problems.examples["Rosenbrock"]

    @test fitness_for_opt(p, 30, 40, 2e4, separable_nes) < 100.0
    @test fitness_for_opt(p, 30, 40, 5e4, xnes) < 100.0

    @test fitness_for_opt(p, 50, 40, 3e4, separable_nes) < 100.0
end

end


text:
# Extensive tests are in ../test/test_Entropy.jl
# TODO:  Write up a high-level description of the functions in this file
#using CSV
using DataFrames
export dist_check, pop_to_dist, pops_to_dist, pops_to_tbl, pop_counts_to_tbl, table_row_to_dist
export entropy, relative_entropy, conditional_entropy, joint_entropy, mutual_information
export sherwin_mutual_information, row_marginal, column_marginal
#include("../../information_theory/src/aliases.jl")

function dist_check( p::DIST_TYPE )
  sum = 0.0
  for x in keys(p)
    if p[x] < 0.0
      error(" value of distribution ", p, " is negative ")
    end
    sum += p[x]
  end
  if !isapprox(sum,1.0)
      error(" sum of distribution ", p, " is not 1.0")
  end
end

# Converts a population to a distribution (i. e., a Dict)
function pop_to_dist( pop::AbstractVector)
  N = length(pop)
  result = DIST_TYPE()
  for i in pop
    #pval = get(result,i,0.0)
    result[i] = get(result,i,0.0) + 1.0/N
    #println("i: ",i,"  pval: ",pval,"  result: ",result)
  end
  dist_check(result)
  result
end

# Assumes that both populations are indexed over the same set.
# Thus, both population must be the same length
# There is one entry in the table for each pair (P1[i],P2[i]).
function pops_to_dist( P1::Population, P2::Population )
  len1 = length(P1)
  len2 = length(P2)
  @assert( len1 == len2 )
  result = DIST_TYPE()
  for i = 1:len1
    entry = (convert(elt_type,P1[i]),convert(elt_type,P2[i]))
    result[entry] = get(result,entry,0.0) + 1.0/len1
  end
  result
end

# There is one column per "allele" of the combined population
# An "allele" is value of an element of the population.
# The sum of the entries in the result is 1.0.
# The sum of the entries in row i is length(P[i])/(sum(length(p) for p in P)
# Note that mutual_information(pops_to_tbl([P1,P2)) does not give the same result
# as mutual_informaiton(P1,P2) which is computed using pops_to_dist(P1,P2).
# The latter result should be the correct result.
function pops_to_tbl( P::PopVect )
  #P = map(x->convert(Vector{Any},x),P)
  m = length(P)   # the number of rows in result
  combined_pop = reduce(vcat, P[i] for i = 1:m )
  N = length(combined_pop)
  allele_list = unique(combined_pop)
  n = length(allele_list)   # The number columns in result
  result = zeros(Float64,m,n)
  for i = 1:length(P)
    dist = Dict( allele_list .=> zeros(Float64,n))
    for p in P[i]
      dist[p] += 1.0/N
    end
    result[i,:] = [ dist[a] for a in allele_list ]
  end
  result
end

# D is a DataFrame with 2 columns per population.
# Each odd numbered column is a list of Strings which are "allele names".
# The next even numbered column is an integer count for the corresponding allele name.
# Each allele name produces a column of the table.
# Example:
# julia> df = DataFrame(names1=["A","B","C"], counts1=[2,5,3], names2=["A","B","D"], counts2=[4,7,1])
# julia> pop_counts_to_tbl(df)
# 2×4 Array{Float64,2}:
#   0.0909091  0.227273  0.136364  0.0
#   0.181818   0.318182  0.0       0.0454545
function pop_counts_to_tbl( D::DataFrame, columns::Vector{Int64}=collect(1:(div(size(D)[2],2))))
  #m = div(length(names(D)),2)   # the number of rows in result
  m = length(columns)
  #println("m: ",m)
  # Convert columns of data frame D from CSV columns to columns of Strings or columns of Int64.
  for i = 1:m
    D[!,2*i-1] = convert(Vector{typeof(D[!,1][1])} ,D[!,2*i-1])
    D[!,2*i] = convert(Vector{Int64},D[!,2*i])
  end
  combined_pop = reduce(vcat, D[!,2*i-1] for i in columns)   # combine allele names
  combined_counts = reduce(vcat, D[!,2*i] for i in columns)  # combine counts
  N = length(combined_pop)
  combined_sum = sum(combined_counts)
  allele_list = unique(combined_pop)
  n = length(allele_list)   # The number columns in result
  #println("m: ",m,"  N: ",N,"  n: ",n,"  combined_sum: ",combined_sum)
  result = zeros(Float64,m,n)
  for i = 1:m
    ii = columns[i]
    #println("i: ",i,"  ii: ",ii)
    dist = Dict( allele_list .=> zeros(Float64,n))
    for j = 1:length(D[!,2*ii-1])
      dist[D[!,2*ii-1][j]] += D[!,2*ii][j]/combined_sum
    end
    result[i,:] = [ dist[a] for a in allele_list ]
  end
  result
end

# Converts a row of a joint probability distribution table to a DIST
# The DIST is normalized so that the probabilities sum to 1.0
function table_row_to_dist( tbl::Array{Float64,2}, row_index::Int64 )
  result = DIST_TYPE()
  (m,n) = size(tbl)
  ssum = sum(tbl[row_index,:])
  for j = 1:n
    result[j] = get(result,j,0.0) + tbl[row_index,j]/ssum
  end
  result
end

function row_marginal( tbl::Array{Float64,2} )
  [ sum(tbl[i,:]) for i = 1:size(tbl)[1]]
end

function column_marginal( tbl::Array{Float64,2} )
  [ sum(tbl[:,j]) for j = 1:size(tbl)[2]]
end

function entropy( p::DIST_TYPE; base::Float64=2.0 )
  result = 0.0
  for x in keys(p)
    result += p[x] > 0.0 ? -p[x]*log(base, p[x] ) : 0.0
    #println("x: ",x, "result_inc: ", p[x] * log(base, p[x]),"  result: ",result)
  end
  result
end

function entropy( p::Population; base::Float64=2.0 )
  entropy( pop_to_dist(p), base=base )
end

# Entropy of a sequence of probabilities.
# entropy(get_probs( v, numinputs))  computes the entropy of a vector v
#   where each v[i] is interpeted as a bit string of length 2^numinputs
function entropy( p::Vector{Float64}; base::Float64=2.0 )
  #println("bit prob entropy")
  entf(x) = x > 0.0 ? -x*log(base,x) : 0.0
  reduce(+,map(entf,p))
end

function entropy( tbl::Array{Float64,2}, row_index::Int64; base::Float64=2.0 )
  entropy( table_row_to_dist( tbl, row_index ), base = base )
end

# relative_entropy()
# Also Kullback Leibler divergence D( q || p )
# Note that if there are any keys x of p such that q[x] == 0, the result will be NaN
#   indicating that the result is not defined.
# Reversed p and q to agree with Cover & Thomas on 9/16/19
function relative_entropy( p::DIST_TYPE, q::DIST_TYPE; base::Float64=2.0 )
  result = 0.0
  for x in keys(p)  # Assume that any keys of q that aren't in p have value 0 for p
    try qval = q[x]
    catch
      return NaN
    end
    #pval = get(q,x,0.0)
    pval = get(p,x,0.0)
    result += pval > 0.0 ? pval*(log(base,pval)-log(base,q[x])) : 0.0
    #println("x: ",x,"  pval: ",pval,"  q[x]: ",q[x],"  result: ",result)
  end
  result
end

function relative_entropy( q::Population, p::Population; base::Float64=2.0 )
  relative_entropy( pop_to_dist(q), pop_to_dist(p), base=base )
end

function relative_entropy( p::Vector{Float64}, q::Vector{Float64}; base::Float64=2.0 )
  result = 0.0
  for i = 1:length(p)
    if q[i] == 0.0
      return NaN
    end
    result += p[i]*log(base,p[i]/q[i])
  end
  result
end

# Conditional entropy of a table.
# H(X|Y)  where  X  is the column-marginal of the table, and Y is the row marginal of the table
# The rows argument specifies the rows of the table that are used.  Order of rows is not important.
# The default is all rows.
# Verification in the file /home/evotech/information_theory/test/test_cover_mutint.jl
function conditional_entropy( tbl::Array{Float64,2}, rows::Vector{Int64}=collect(1:size(tbl)[1]); base::Float64=2.0 )
  #println("rows: ",rows)
  row_sums = map(x->sum(tbl[x,:]),rows)
  tbl_sum = sum(row_sums)   # sum of the part of the table that we are using.
  result = 0.0
  for x = 1:length(rows)
    #sum_px = sum(tbl[rows[x],:])
    for y = 1:size(tbl)[2]
      pyx = tbl[rows[x],y]/row_sums[x]
      result -= tbl[rows[x],y]/tbl_sum > 0.0 ? tbl[rows[x],y]/tbl_sum*log(base,pyx) : 0.0
      #println("(x,y): ",(x,y),"  pyx: ",pyx,"  inc: ",tbl[rows[x],y]>0.0 ? -tbl[rows[x],y]*log(2,pyx) : 0.0)
    end
  end
  result
end

# joint entropy by definition of Cover and Thomas equation 2.8
function joint_entropy( tbl::Array{Float64,2}; base::Float64=2.0 )
  if ! reduce( &, tbl .>= 0.0 )
    error("all entries of tbl must be non-netative in function joint entropy")
  end
  if !(sum(tbl) ≈ 1.0)
    tbl = tbl/sum(tbl)   # Normalize to sum 1
  end
  (m,n) = size(tbl)
  entf(x) = x != 0 ? -x*log(base,x) : 0.0
  reduce(+, map(entf,tbl))
end

# Assumes that both populations are indexed over the same set.
# Thus, both population must be the same length
# There is one entry in the table for each pair (P1[i],P2[j]).
function joint_entropy( P1::Population, P2::Population; base::Float64=2.0 )
  @assert length(P1) == length(P2)
  entropy(pops_to_dist(P1,P2),base=base)
end

# Mutual information by the defintion   I(P1;P2) = H(P1) + H(P2) - H(P1,P2) (eq. 2.45 of Cover)
# Uses the pops_to_dist() method in joint_entropy()
# This is the "standard" (non-Sherwin) definition of mutual information
function mutual_information( P1::Population, P2::Population; base::Float64=2.0 )
  @assert length(P1) == length(P2)
  entropy(P1,base=base) + entropy(P2,base=base) - joint_entropy(P1,P2,base=base)
end

# Mutual information by the "standard" defintion of eq. 2.28 of Cover and Thomas.
#   This is the mutual information of a joint distribution.
function mutual_information( tbl::Array{Float64,2}; base::Float64=2.0 )
  if ! reduce( &, tbl .>= 0.0 )
    error("all entries of tbl must be non-negative in function mutual information")
  end
  if !(sum(tbl) ≈ 1.0)
    tbl = tbl/sum(tbl)   # Normalize to sum 1
  end
  (m,n) = size(tbl)
  mif(i,j) = tbl[i,j] != 0.0 ? tbl[i,j]*log(base,tbl[i,j]/mr1[i]/mr2[j]) : 0.0
  mr1 = [ sum( tbl[i,:]) for i = 1:m ]
  mr2 = [ sum( tbl[:,j]) for j = 1:n ]
  reduce( +, [ mif(i,j) for i=1:m, j=1:n ] )
end

# tbl is the probability distribution for a total population which is made up of m subpopulations
#    all over the same allele set of size n.
# Row i of tbl (normalized to sum to 1) corresponds to a probability distribution over subpopulation i.
function sherwin_mutual_information( tbl::Array{Float64,2}; base::Float64=2.0 )
  @assert isapprox( sum(tbl), 1.0 )
  entf(x) = x != 0 ? -x*log(base,x) : 0.0
  (m,n) = size(tbl)
  row_weights = row_marginal(tbl)
  row_ents = [ sum( map(entf,tbl[i,:]/row_weights[i])) for i = 1:m ] # entropies of normalized rows
  #println("row_ents: ",row_ents)
  col_ents = map(entf,[sum(tbl[:,j]) for j = 1:n])   # entropies of the sums of columns
  #println("col_ents: ",col_ents)
  ent = sum(col_ents) - sum(row_ents[i]*row_weights[i] for i = 1:m)
  #println("ent: ",ent)
  ent
end

function mutual_information( D::DataFrame, columns::Vector{Int64}=collect(1:(div(size(D)[2],2))); base::Float64=2.0 )
  mutual_information(pop_counts_to_tbl( D, columns ), base=base )
end

function mutual_information( P::Population; base::Float64=2.0 )
  sherwin_mutual_information( pops_to_tbl( P ), base=base )
end

function mutual_information( P1::FPopulation, P2::FPopulation; base::Float64=2.0 )
  #println("MI pops_to_tbl")
  sherwin_mutual_information( pops_to_tbl( [P1, P2] ), base=base )
end

function mutual_information( P::PopVect; base::Float64=2.0 )
  #println("MI pops_to_tbl")
  sherwin_mutual_information( pops_to_tbl( P ), base=base )
end


text:
module TidyStanza

export combine, pivot_wider

include("Across.jl") # defines the Across type

include("relocate.jl")
include("any_of.jl")
include("last_col.jl")

using DataFrames

include("Pivot.jl") # defines pivoting functions

import DataFrames: combine

everything() = true
returntrue() = true

"""
    A function to get input col names from Across.innames
"""
function get_input_cols(df, across_innames::Where)
    names(df)[across_innames.fn.(eachcol(df))]
end

function get_input_cols(df::GroupedDataFrame, across_innames::Where)
    names(df)[across_innames.fn.(eachcol(parent(df)))]
end

function get_input_cols(df, across_innames)
    filter(across_innames, names(df))
end

function DataFrames.transform(df::AbstractDataFrame, da::Across)
    input_cols = get_input_cols(df, da.innames)

    if isnothing(da.outnames)
        pairs = [col => fn for col in input_cols, fn in da.fn]
        return transform(df, pairs...)
    else
        outcols = [
            replace(replace(da.outnames, "{col}" => col), "{fn}" => string(fn))
            for col in input_cols, fn in keys(da.fn)
        ]

        in_fn_out_iterator = zip(Iterators.product(input_cols, da.fn), outcols)

        pairs = [col => fn => outcol for ((col, fn), outcol) in in_fn_out_iterator]

        return transform(df, pairs...)
    end
end

function dplyr_across_in_combines(df, da)
    input_cols = get_input_cols(df, da.innames)

    if isnothing(da.outnames)
        pairs = [col => fn for col in input_cols, fn in da.fn]
        return combine(df, pairs...)
    else
        outcols = [
            replace(replace(da.outnames, "{col}" => col), "{fn}" => string(fn))
            for col in input_cols, fn in keys(da.fn)
        ]

        in_fn_out_iterator = zip(Iterators.product(input_cols, da.fn), outcols)

        pairs = [col => fn => outcol for ((col, fn), outcol) in in_fn_out_iterator]

        return combine(df, pairs...)
    end
end

DataFrames.combine(df, da::Across) = dplyr_across_in_combines(df, da)

DataFrames.combine(df::T, da::Across) where {T<:AbstractDataFrame} =
    dplyr_across_in_combines(df, da)

end


text:
#=

  Secret Santa problem II in Julia ConstraintSolver.jl

  From Maple Primes: "Secret Santa Graph Theory"
  http://www.mapleprimes.com/blog/jpmay/secretsantagraphtheory
  """
  Every year my extended family does a "secret santa" gift exchange.
  Each person draws another person at random and then gets a gift for
  them. At first, none of my siblings were married, and so the draw was
  completely random. Then, as people got married, we added the restriction
  that spouses should not draw each others names. This restriction meant
  that we moved from using slips of paper on a hat to using a simple
  computer program to choose names. Then people began to complain when
  they would get the same person two years in a row, so the program was
  modified to keep some history and avoid giving anyone a name in their
  recent history. This year, not everyone was participating, and so after
  removing names, and limiting the number of exclusions to four per person,
  I had data something like this:

  Name: Spouse, Recent Picks

  Noah: Ava. Ella, Evan, Ryan, John
  Ava: Noah, Evan, Mia, John, Ryan
  Ryan: Mia, Ella, Ava, Lily, Evan
  Mia: Ryan, Ava, Ella, Lily, Evan
  Ella: John, Lily, Evan, Mia, Ava
  John: Ella, Noah, Lily, Ryan, Ava
  Lily: Evan, John, Mia, Ava, Ella
  Evan: Lily, Mia, John, Ryan, Noah
  """

  Note: I interpret this as the following three constraints:
    1) One cannot be a Secret Santa of one's spouse
    2) One cannot be a Secret Santa for somebody two years in a row
    3) Optimization: maximize the time since the last time

  This model also handle single persons, something the original
  problem don't mention.

  Model created by Hakan Kjellerstrand, hakank@gmail.com
  See also my Julia page: http://www.hakank.org/julia/

=#

using ConstraintSolver, JuMP
using Cbc, GLPK, Ipopt
const CS = ConstraintSolver
include("constraints_utils.jl")

function secret_santa2(n,M,people,spouses,rounds,print_solutions=true,all_solutions=true)

    cbc_optimizer = optimizer_with_attributes(Cbc.Optimizer, "logLevel" => 0)
    glpk_optimizer = optimizer_with_attributes(GLPK.Optimizer)
    ipopt_optimizer = optimizer_with_attributes(Ipopt.Optimizer)

    model = Model(optimizer_with_attributes(CS.Optimizer,   # "all_solutions"=> all_solutions,
                                                            "all_optimal_solutions"=>all_solutions,
                                                            "logging"=>[],

                                                            "traverse_strategy"=>:BFS,
                                                            # "traverse_strategy"=>:DFS,
                                                            # "traverse_strategy"=>:DBFS,

                                                            # "branch_split"=>:Smallest,
                                                            # "branch_split"=>:Biggest,
                                                            "branch_split"=>:InHalf,

                                                            # "simplify"=>false,
                                                            # "simplify"=>true, # default

                                                            "time_limit"=>6,

                                                            # "lp_optimizer" => cbc_optimizer,
                                                            "lp_optimizer" => glpk_optimizer,
                                                            # "lp_optimizer" => ipopt_optimizer,
                                        ))

    @variable(model, 1 <= santas[1:n] <= n, Int)
    @variable(model, 1 <= santas2[1:n] <= n, Int)

    @variable(model, 1 <= santa_distance[1:n] <= n+1, Int)
    @variable(model, 0 <= z <= 1000, Int)

    # Everyone gives and receives a Secret Santa
    @constraint(model, santas in CS.AllDifferentSet())

    # no Santa for a spouses (or him-/herself)
    for i in 1:n
        @constraint(model,santas[i] != i)
        if spouses[i] > 0
            @constraint(model,santas[i] != spouses[i])
        end
    end

    # optimize "distance" to earlier rounds:
    for i in 1:n
        # SantaDistance[I] == Rounds[I,Santas[I]]
        my_element(model,santas[i],rounds[i,:],santa_distance[i])
    end

    # Cannot be a Secret Santa for the same person two years in a row.

    for i in 1:n
       # Rounds[I,Santas2[I]] == 1,
       my_element(model,santas2[i],rounds[i,:],1)
       @constraint(model,santas[i] != santas2[i])
    end

    @constraint(model,z == sum(santa_distance))

    # Optimize the sum of santa_distance
    @objective(model,Max,z)

    # Solve the problem
    optimize!(model)

    status = JuMP.termination_status(model)
    if status == MOI.OPTIMAL
        num_sols = MOI.get(model, MOI.ResultCount())
        # println("num_sols:$num_sols\n")
        if print_solutions
            for sol in 1:num_sols
                println("solution #$sol")
                z_val = convert.(Integer,JuMP.value.(z; result=sol))
                santas_val = convert.(Integer,JuMP.value.(santas; result=sol))
                santas2_val = convert.(Integer,JuMP.value.(santas2; result=sol))
                santa_distance_val = convert.(Integer,JuMP.value.(santa_distance; result=sol))
                println("z:$z_val")
                println("santas:$santas_val")
                println("santas2:$santas2_val")
                println("santa_distance:$santa_distance_val")
                println()

            end
        end
        println("num_sols:$num_sols\n")
    else
        println("status:$status")
    end

    return status
end

#
# Test for n=8 (without single) or n= 9 with a single person
#
function test_secret_santa2(n=9,print_solutions=true,all_solutions=false)
    println("n:$n")
    # n = 8 # Without Single person
    # n = 9 # With a Single person
    noah,ava,ryan,mia,ella,john,lily,evan,single = 1:9
    people = ["noah","ava","ryan","mia","ella","john","lily","evan","single"]

    spouses = [
             ava,  # noa
             noah, # ava
             mia,  # rya
             ryan, # mia
             john, # ella
             ella, # john
             evan, # lily
             lily  # evan
             , 0   # single has no spouse
       ]

     M = n+1 # "large M" to indicate no earlier history

    #
    # The matrix version of earlier rounds.
    # M means that no earlier Santa.
    # Note: Ryan and Mia has the same recipient for years 3 and 4,
    #       and Ella and John has for year 4.
    #       This seems to be caused by modification of
    #       original data.
    #
    # rounds with a single person (fake data)
    #
    rounds = [
    #N  A  R  M  El J  L  Ev S
    [0, M, 3, M, 1, 4, M, 2, 2], # Noah
    [M, 0, 4, 2, M, 3, M, 1, 1], # Ava
    [M, 2, 0, M, 1, M, 3, 4, 4], # Ryan
    [M, 1, M, 0, 2, M, 3, 4, 3], # Mia
    [M, 4, M, 3, 0, M, 1, 2, M], # Ella
    [1, 4, 3, M, M, 0, 2, M, M], # John
    [M, 3, M, 2, 4, 1, 0, M, M], # Lily
    [4, M, 3, 1, M, 2, M, 0, M], # Evan
    [1, 2, 3, 4, M, 2, M, M, 0]  # Single
    ]

    secret_santa2(n,M,people,spouses,resize_matrix(rounds),print_solutions,all_solutions)

end

# First we solve for first optimal solutions
#=
z:67
santas:[4, 7, 1, 6, 2, 8, 3, 5]
santas2:[5, 8, 5, 2, 7, 1, 6, 4]
santa_distance:[9, 9, 9, 9, 4, 9, 9, 9]
  1.516136 seconds (24.50 M allocations: 552.591 MiB, 10.32% gc time)

z:90
santas:[7, 5, 6, 1, 9, 4, 3, 2, 8]
santas2:[5, 9, 5, 2, 7, 1, 6, 4, 1]
santa_distance:[10, 10, 10, 10, 10, 10, 10, 10, 10]
  1.831563 seconds (29.13 M allocations: 652.884 MiB, 11.81% gc time)

3.794155 seconds (54.49 M allocations: 1.221 GiB, 10.56% gc time)
=#
@time test_secret_santa2(8,true,false)
@time test_secret_santa2(9,true,false)

# Then the total number of solutions
#=
n:8
num_sols:8
  2.428113 seconds (39.62 M allocations: 904.140 MiB, 9.26% gc time)

n:9
num_sols:24
  3.991528 seconds (58.10 M allocations: 1.314 GiB, 17.37% gc time)

10.229509 seconds (152.62 M allocations: 3.439 GiB, 12.86% gc time)
=#
@time test_secret_santa2(8,false,true)
@time test_secret_santa2(9,false,true)

max_stars_repo_name: Wikunia/hakank

max_issues_repo_name: Wikunia/hakank

max_forks_repo_name: Wikunia/hakank


text:
# testset for the Fridge.jl package
@testset "Fridge" begin
    using Fridge

    testList = ["cheese","potato","tomato","cabbage"]
    testDict = Dict("fries" => ["potato","salt"],
        "macaroni" => ["cheese", "salt", "jambon"],
        "salad" => ["cabbage","tomato","oil","horseradish"],
        "chocolatemilk" => ["chocolate","milk"])

    testDict2 = Dict("fries" => ["potato","salt"],
        "macaroni" => ["cheese", "salt", "jambon"])

    testDict3 = Dict("fries" => ["potato","salt"],
        "macaroni" => ["cheese", "salt", "jambon"],
        "salad" => ["cabbage","tomato","oil","horseradish"])

    testDict4 = Dict("fries" => ["potato","salt"],
    "macaroni" => ["cheese", "salt", "jambon"],
    "salad" => ["cabbage","tomato","oil","horseradish"],
    "chocolatemilk" => ["chocolate","milk"],
    "waffle" => ["sugar", "eggs", "milk"],
    "thea" => ["herbs", "water"],
    "coffee" => ["coffeebeans", "water"])

    testSolution = Dict("salad" => [0,0,1,1,2], "fries" => [0,1,0,0,1])

    @testset "checkIngredients" begin
        # check if all ingredients are passed if no replacement is needed
        @test checkIngredients(testList,["cheese","salt","tomato","chocolate"]) == testList

        # test for type
        @test checkIngredients(testList,["cheese","salt","tomato","chocolate"]) isa Vector{String}
    end

    @testset "greedyFindCombo" begin
        # test max number of recipes end
        @test greedyFindCombo(testList, testDict, 2) == Dict("salad" => [0,0,1,1,2], "fries" => [0,1,0,0,1])

        # test all ingredients used end
        @test greedyFindCombo(testList, testDict, 3) == Dict("macaroni" => [1,0,0,0,2], "salad" => [0,0,1,1,2], "fries" => [0,1,0,0,1])

        # test no recipes left end
        @test greedyFindCombo(testList, testDict2, 3) == Dict("macaroni" => [1,0,0,0,2], "fries" => [0,1,0,0,1])

        # test no recipes left that use ingredients from the fridge
        @test greedyFindCombo(testList, testDict, 4) == Dict("macaroni" => [1,0,0,0,2], "salad" => [0,0,1,1,2], "fries" => [0,1,0,0,1])
    end

    @testset "randomCombo" begin
        # test all ingredients used end
        @test randomCombo(testList, testDict3, 3) == Dict("macaroni" => [1,0,0,0,2], "salad" => [0,0,1,1,2], "fries" => [0,1,0,0,1])

        # test no recipes left end
        @test randomCombo(testList, testDict2, 3) == Dict("macaroni" => [1,0,0,0,2], "fries" => [0,1,0,0,1])

        # test length of max number recipes end
        @test length(values(randomCombo(testList, testDict, 2))) == 2
    end

    @testset "Neighbour" begin
        # test not random with empty tabuList
        @test "macaroni" in keys(Neighbour(testSolution, testList, testDict, 2, [], false))

        # test not random with tabuList
        @test !in("macaroni",keys(Neighbour(testSolution, testList, testDict, 2, ["macaroni"], false)))

        # test length of random with empty tabuList
        @test length(values(Neighbour(testSolution, testList, testDict, 2, [], true))) == 2
    end

    @testset "SAFindCombo" begin
        # test output type
        @test typeof(SAFindCombo(testSolution,  testList, testDict4, 2, true, tabuLength=1)) == Dict{Any, Any}

        # test output length
        @test length(values(SAFindCombo(testSolution,  testList, testDict4, 2, true, tabuLength=1))) == 2
    end

    @testset "scrapeRecipe" begin
        # check if .csv file is produced
        scrapeRecipe(106541,106542,"./testDB.csv")
        @test isfile("./testDB.csv")
    end

    @testset "loadRecipeDBCSV" begin
        # check if downloaded recipes are correct
        testDictDB = loadRecipeDBCSV("./testDB.csv")
        @test keys(testDictDB) == keys(Dict("Belgian Chocolate Mousse" => ["test"], "Belgian Buns" => ["test"]))
    end

    @testset "findBestRecipe" begin
        scrapeRecipe(106543,106555,"./testDB.csv")
        # check if the combination of all recipes works
        @test length(values(findBestRecipe(testList, "./testDB.csv", numRecipes=2, randRecipe=true, testmode=true))) == 0
    end
end


text:
# # Splitting a string into words and counting them in parallel
#
# We start from the parallel algorithm presented in Guy Steele's 2009
# ICFP talk [(video)](https://vimeo.com/6624203).  It splits a
# space-separated string into list of strings (words).  The repeating
# theme in the talk was to build "singleton solutions" and then merge
# them together using an associative function.  We will follow this
# guideline and slightly extend the algorithm.
#
# It is highly recommended to just watch the talk for understanding
# the algorithm.  However, we briefly describe how it works.
#
# When a certain contiguous region of a string is processed, we either
# already have seen at least one space or not.  These two states are
# tracked using following two types.  If there is no space so far, we
# only have a chunk of a possibly larger word (see example below):

struct Chunk
    s::String
end

# If there are one or more spaces, (possibly zero) words that are
# already determined and left/right "chunks" have to be tracked
# separately:

struct Segment
    l::String
    A::Vector{String}
    r::String
end

# Here is an example taken from the talk:
#
# ```
# Segment("Here", ["is", "a"], "")
#    |
#    |       Segment("lian", [], "string")
#  __|_____            _|______
# |        |          |        |
# Here is a sesquipedalian string of words
#           |________|          |________|
#    Chunk("sesquipeda")        Segment("g", ["of"], "words")
# ```

# We then need a way to merge two results which can independently in
# one of the above two states.

⊕(x::Chunk, y::Chunk) = Chunk(x.s * y.s)
⊕(x::Chunk, y::Segment) = Segment(x.s * y.l, y.A, y.r)
⊕(x::Segment, y::Chunk) = Segment(x.l, x.A, x.r * y.s)
⊕(x::Segment, y::Segment) =
    Segment(x.l,
            append!(append!(x.A, maybewordv(x.r * y.l)), y.A),
            y.r)

maybewordv(s::String) = isempty(s) ? String[] : [s]
nothing  # hide

# Input is a sequence of `Char`s.  Each of them has to be converted
# into a "singleton solution" which can be merged with already
# aggregated (or another singleton) solution with `⊕`:

segmentorchunk(c::Char) = c == ' ' ? Segment("", [], "") : Chunk(string(c))
nothing  # hide

# Putting them together, we get:

function collectwords(s::String)
    g = mapfoldl(segmentorchunk, ⊕, s; init=Segment("", [], ""))
    if g isa Char
        return maybewordv(g.s)
    else
        return append!(append!(maybewordv(g.l), g.A), maybewordv(g.r))
    end
end
nothing  # hide

# Let's run a few tests covering some edge cases:

using Test
@testset begin
    @test collectwords("This is a sample") == ["This", "is", "a", "sample"]
    @test collectwords(" Here is another sample ") == ["Here", "is", "another", "sample"]
    @test collectwords("JustOneWord") == ["JustOneWord"]
    @test collectwords(" ") == []
    @test collectwords("") == []
end
nothing  # hide

# ## String-splitting transducer
#
# Let's try to make it re-usable by packaging it into transducers.

using Transducers

# Rather than accumulating words into a vector, we are going to write
# a transducer that "emits" a word as soon as it is ready.  The
# downstream transducer may choose to record everything or only
# aggregate, e.g., reduced statistics.  To this end, we replace
# `Segment` in the original algorithm to

struct Vacant
    l::String
    r::String
end

# and output the words in the "middle" without accumulating it.  We
# use [`ScanEmit`](@ref) which requires an operator/function like `⊕`
# above but returning a pair of output and next state.  This function
# (`extract` below) must have the signature `(S, S) -> (O, S)` where
# `S` is the type for accumulated state and input and `O` is the
# output type.

extract(x::Chunk, y::Chunk) = (), Chunk(x.s * y.s)
extract(x::Chunk, y::Vacant) = (), Vacant(x.s * y.l, y.r)
extract(x::Vacant, y::Chunk) = (), Vacant(x.l, x.r * y.s)
extract(x::Vacant, y::Vacant) = maybewordt(x.r * y.l), Vacant(x.l, y.r)

maybewordt(s) = isempty(s) ? () : (s,)
nothing  # hide

# `maybewordt(x.r * y.l)` in `extract(x::Vacant, y::Vacant)` is the
# "emission".
#
# The words at the beginning and/or the end are not handled by
# `extract`.  This must be handled separately:

lastword(x::Chunk) = maybewordt(x.s)
lastword(x::Vacant) = (maybewordt(x.r)..., maybewordt(x.l)...)

vacantorchunk(c::Char) = c == ' ' ? Vacant("", "") : Chunk(string(c))

wordsxf = Map(vacantorchunk) |> ScanEmit(extract, Chunk(""), lastword) |> Cat()

# Test:

@testset begin
    @test collect(wordsxf, "This is a sample") == ["is", "a", "sample", "This"]
    @test collect(wordsxf, " Here is another sample ") == ["Here", "is", "another", "sample"]
    @test collect(wordsxf, "JustOneWord") == ["JustOneWord"]
    @test collect(wordsxf, " ") == []
    @test collect(wordsxf, "") == []
end
nothing  # hide

# Side note: In the first example, the first word `This` comes last.
# This is actually expected since both `.l` and `.r` are flushed in
# `lastword` which is called at the very end.  Here, `This` is stored
# in `.l` field.  If the order of the words is important, there are
# many possible fixes.  For example, `extract` and `lastword` can
# bundle information about the origin of the word (left vs
# middle-or-right).  Alternatively, perhaps the easiest solution is to
# insert a space in the beginning of input data.
#
# ## Word-counting transducer
#
# We can pipe the resulting words into various transducers.

processcount(word) = Base.ImmutableDict(word => 1)
countxf = wordsxf |> Map(processcount)

# Transducer `countxf` constructs a "singleton solution" as a
# dictionary which then accumulated with the associative reducing step
# function `mergecont!`:

mergecont!(a, b) = merge!(+, a, b)
nothing  # hide

# Putting the transducer and reducing function together, we get:

countwords(s; kwargs...) =
    reduce(mergecont!,
           countxf,
           collect(s);
           init = CopyInit(Dict{String,Int}()),
           kwargs...)
nothing  # hide

# Side note: Since [`reduce`](@ref) does not support string, the input
# string is converted to a `Vector{Char}` first by `collect`.
#
# Side note 2: We use [`CopyInit`](@ref) to create a fresh initial
# state for each sub-reduce to avoid overwriting mutable data between
# threads.
#
# Side note 3: [`reduce`](@ref) wraps `mergecont!` automatically with
# [`Completing`](@ref).  This is why `mergecont!` does not have to
# have the unary method.

# Let's run some tests with different `basesize` (`length(s) /
# basesize` corresponds to number of tasks to be used):

@testset for basesize in [1, 2, 4]
    @test countwords("This is a sample", basesize=basesize) ==
        Dict("This" => 1, "is" => 1, "a" => 1, "sample" => 1)
    @test countwords(" Here is another sample ", basesize=basesize) ==
        Dict("Here" => 1, "is" => 1, "another" => 1, "sample" => 1)
    @test countwords("JustOneWord", basesize=basesize) ==
        Dict("JustOneWord" => 1)
    @test countwords(" ", basesize=basesize) == Dict()
    @test countwords("", basesize=basesize) == Dict()
    @test countwords("aaa bb aaa aaa bb bb aaa", basesize=basesize) ==
        Dict("aaa" => 4, "bb" => 3)
    @test countwords("あああ いい あああ あああ いい いい あああ", basesize=basesize) ==
        Dict("あああ" => 4, "いい" => 3)
end
nothing  # hide


text:
xml = readxml("data/math.xml").root

io = IOBuffer()
print_tree(io, xml)
ser = String(take!(io))
str = """
math
└─ apply
   ├─ times
   ├─ ci
   ├─ ci
   └─ ci
"""
@test strip(ser) == strip(str)


text:
"""
    g1_DA_semiinf_CW(τ::AbstractVector, ρ, μa, μsp; BFi = 2e-8, β = 1.0, n_ext = 1.0, n_med = 1.0, z = 0.0, λ = 700.0)

Compute the electric field autocorrelation function function g1 in a semi-infinite geometry (CW).

# Arguments
- `τ`: time lags (s)
- `ρ`: the source detector separation (cm⁻¹)
- `μa`: absorption coefficient (cm⁻¹)
- `μsp`: reduced scattering coefficient (cm⁻¹)

# Keyword arguments
- `BFi`: Blood flow index ~αDb (cm²/s)
- `β`: constant in Siegert relation dependent on collection optics
- `n_ext`: the boundary's index of refraction (air or detector)
- `n_med`: the sample medium's index of refraction
- `z`: the z-depth orthogonal from the boundary (cm)
- `λ`: wavelength (nm)
"""
function g1_DA_semiinf_CW(τ::AbstractVector, ρ, μa, μsp; BFi = 2e-8, n_ext = 1.0, n_med = 1.0, z = 0.0, λ = 700.0)
    g1 = similar(τ)
    k0 = 2 * n_med * π / (λ * 1e-7) # this should be in cm
    tmp = 2 * μsp * BFi * k0^2

    G0 = fluence_DA_semiinf_CW(ρ, μa, μsp, n_ext = n_ext, n_med = n_med, z = z)

    for ind in eachindex(τ)
        μa_dynamic = μa + tmp * τ[ind]
        G1 = fluence_DA_semiinf_CW(ρ, μa_dynamic, μsp, n_ext = n_ext, n_med = n_med, z = z)
        g1[ind] = G1 / G0
    end

    return g1
end

"""
    g1_DA_semiinf_TD(τ::AbstractVector, t::AbstractFloat, ρ, μa, μsp; BFi = 2e-8, β = 1.0, n_ext = 1.0, n_med = 1.0, z = 0.0, λ = 700.0)
    g1_DA_semiinf_TD(τ::AbstractVector, t::Abstractector, ρ, μa, μsp; BFi = 2e-8, β = 1.0, n_ext = 1.0, n_med = 1.0, z = 0.0, λ = 700.0)

Compute the time-domain electric field autocorrelation function function g1 in a semi-infinite geometry.

If t is a AbstractVector then g1(tau) will be calculated by ∫P(t)g1(tau, t)dt with integration bounds given by t[1] and t[end].
This is needed for time-gated applications.

# Arguments
- `τ`: time lags (s)
- `t`: photon time of flight (related to photon path length)
- `ρ`: the source detector separation (cm⁻¹)
- `μa`: absorption coefficient (cm⁻¹)
- `μsp`: reduced scattering coefficient (cm⁻¹)

# Keyword arguments
- `BFi`: Blood flow index ~αDb (cm²/s)
- `β`: constant in Siegert relation dependent on collection optics
- `n_ext`: the boundary's index of refraction (air or detector)
- `n_med`: the sample medium's index of refraction
- `z`: the z-depth orthogonal from the boundary (cm)
- `λ`: wavelength (nm)
- `N_quad`: number of nodes in gauss-legendre integration
"""
function g1_DA_semiinf_TD(τ::AbstractVector, t::AbstractFloat, ρ, μa, μsp; BFi = 2e-8, n_ext = 1.0, n_med = 1.0, z = 0.0, λ = 700.0)
    g1 = similar(τ)
    k0 = 2 * n_med * π / (λ * 1e-7) # this should be in cm
    tmp = 2 * μsp * BFi * k0^2

    G0 = fluence_DA_semiinf_TD(t, ρ, μa, μsp, n_ext = n_ext, n_med = n_med, z = z)

    for ind in eachindex(τ)
        μa_dynamic = μa + tmp * τ[ind]
        G1 = fluence_DA_semiinf_TD(t, ρ, μa_dynamic, μsp, n_ext = n_ext, n_med = n_med, z = z)
        g1[ind] = G1 / G0
    end

    return g1
end
function g1_DA_semiinf_TD(τ::AbstractVector, t::AbstractVector, ρ, μa, μsp; BFi = 2e-8, n_ext = 1.0, n_med = 1.0, z = 0.0, λ = 700.0, N_quad = 50)
    g1 = similar(τ)
    k0 = 2 * n_med * π / (λ * 1e-7) # this should be in cm
    tmp = 2 * μsp * BFi * k0^2
    x, w = gausslegendre(N_quad)

    tpsf_norm = integrate(t -> fluence_DA_semiinf_TD(t, ρ, μa, μsp, n_ext = n_ext, n_med = n_med, z = z), x, w, t[1], t[2])

    Threads.@threads for ind in eachindex(τ)
        μa_dynamic = μa + tmp * τ[ind]
        g1[ind] = integrate(t -> (fluence_DA_semiinf_TD(t, ρ, μa_dynamic, μsp, n_ext = n_ext, n_med = n_med, z = z) / tpsf_norm), x, w, t[1], t[2])
    end

    return g1
end

function g2_DA_semiinf_CW(τ::AbstractVector, ρ, μa, μsp; BFi = 2e-8, n_ext = 1.0, n_med = 1.0, z = 0.0, λ = 700.0, β = 1.0)
    g1 = g1_DA_semiinf_CW(τ, ρ, μa, μsp, BFi = BFi, n_ext = n_ext, n_med = n_med, z = z, λ = λ)
    for ind in eachindex(τ)
        g1[ind] = 1 + β * abs(g1[ind])^2
    end
    return g1
end

function g2_DA_semiinf_TD(τ::AbstractVector, t::Union{AbstractVector, AbstractFloat}, ρ, μa, μsp; BFi = 2e-8, n_ext = 1.0, n_med = 1.0, z = 0.0, λ = 700.0, β = 1.0, N_quad = 50)
    if t isa AbstractVector
        g1 = g1_DA_semiinf_TD(τ, t, ρ, μa, μsp, BFi = BFi, n_ext = n_ext, n_med = n_med, z = z, λ = λ, N_quad = N_quad)
    elseif t isa AbstractFloat
        g1 = g1_DA_semiinf_TD(τ, t, ρ, μa, μsp, BFi = BFi, n_ext = n_ext, n_med = n_med, z = z, λ = λ)
    end

    for ind in eachindex(τ)
        g1[ind] = 1 + β * abs(g1[ind])^2
    end
    return g1
end

function integrate(f::Function, x, w, lb, ub)
    out = zero(eltype(x))
    for m in eachindex(x)
        out += f(x[m] * (ub - lb) / 2 + (lb + ub) / 2) * w[m]
    end
    return out * (ub - lb) / 2
end

# use this function when it is more efficient to compute y points all together for layered media (Laplace contour is fixed for each time point)
function integrate_compute_y_first(f::Function, x, w, lb, ub)
    x = @. x * (ub - lb) / 2 + (lb + ub) / 2
    y = f(x) # the function should be able to take a vector and return a vector
    out = zero(eltype(x))
    for m in eachindex(x)
        out += y[m] * w[m]
    end
    return out * (ub - lb) / 2
end


text:
@testset "Covering" begin
  @testset "Rectangle" begin
    r = cover(RegularGrid{Float64}(100, 200), RectangleCoverer())
    @test r == RectangleRegion((0.,0.), (99.,199.))

    r = cover(PointSet([0. 1. 2.; 0. 2. 1.]), RectangleCoverer())
    @test r == RectangleRegion((0.,0.), (2.,2.))
  end
end


text:
export TensileBar, TensileNotchedBar

function TensileNotchedBar(args...; kwargs...)
    TensileBar(args...; notched=true, kwargs...)
end

function TensileBar(; gen_mat=nothing, spc_mat=nothing, resolution=0.1, solver_=:qs, notched=false, file_prefix="TensileBar")
    if solver_ in [:qs]
        Steps = 10
        sargs = (Steps, 0.0001)
        skwargs = Dict(:max_iter=>500)
        fwf = 1
    else
        Steps = 1000
        sargs = (Steps, )
        skwargs = Dict()
        fwf = 10
    end
    solver = PDBenchmark.NameParam(solver_, sargs, Dict(:filewrite_freq=>fwf, :neigh_update_freq=>100, :file_prefix=>file_prefix, :start_at=>0, skwargs...))

    if notched
        obj = () -> NotchedBar(bounds=[0 12; 0 1; 0 1], notch=[5.75 6.25; 0.0 0.25; 0.0 1.0])
    else
        obj = () -> Bar(bounds=[0 12; 0 1; 0 1])
    end

    reso = resolution
    geom = PDBenchmark.NameParam(obj, (), Dict(:resolution=>reso))

    if solver_ in [:qs]
        bc_f1 = env -> (vec(env.y[1,:] .<= 01.0),  [-0.01, 0, 0])
        bc_f2 = env -> (vec(env.y[1,:] .>= 11.0), [0.01, 0, 0])
        bc = [  PDBenchmark.NameParam(:MoveBC, (bc_f1,)),
                PDBenchmark.NameParam(:MoveBC, (bc_f2,)),
                ]
    else
        bc_f1 = env -> (vec(env.y[1,:] .<= 1.0),  [-0.0001, 0, 0])
        bc_f2 = env -> (vec(env.y[1,:] .>= 11.0), [0.0001, 0, 0])
        bc = [PDBenchmark.NameParam(:MoveBC, (bc_f1,)), PDBenchmark.NameParam(:MoveBC, (bc_f2,))]
    end
    RM = [PDBenchmark.NameParam(:LinearRepulsionModel, (100.0,), Dict(:blocks=>(1,), :distanceX=>3, :max_neighs=>200))]

    f = env -> begin
        left_ = (env.y[1,:] .<= 1.0)
        right_ = (env.y[1,:] .>= 11.0)
        right_max = maximum(env.y[1, right_])
        env.Params = Dict("left" => left_, "right"=>right_, "right_max"=>right_max)
        env.Out = Dict("Force" => zeros(3, Steps), "Displacement" => zeros(1, Steps))
        env.Collect! = (env, step) -> begin
                env.Out["Force"][:, step] = sum(env.f[:,env.Params["left"]], dims=2)
                env.Out["Displacement"][1, step] = 2*(maximum(env.y[1,env.Params["right"]]) - env.Params["right_max"])
                if step%fwf==0
                    println("E = $(10.0*env.Out["Force"][1, step]/env.Out["Displacement"][1, step])")
                end
            end
    end

    test = PDBenchmark.Test(solver, geom, gen_mat, spc_mat, bc, RM, f)

    return test
end

max_stars_repo_name: ravinderbhattoo/PDBenchmark

max_issues_repo_name: ravinderbhattoo/PDBenchmark

max_forks_repo_name: ravinderbhattoo/PDBenchmark


text:
"""
    Schedule(ts, values, interpolations; unit = :epoch; kwargs...)

Describes how the values of a hyperparameter change over the
training.

## Examples

```julia
# one-cycle schedule over 10 epochs from 0.01 over 0.1 to 0.001
Schedule([0, 2, 10], [0.01, 0.1, 0.001], [sineio(), sineio()])
```

`Schedule` is an alias for `Animations.Animation`, see the
[documentation](https://jkrumbiegel.github.io/Animations.jl/dev)
for more detailed information.
"""
const Schedule = Animation

"""
    Scheduler(schedules...)

Callback for hyperparameter scheduling. Takes pairs of [`HyperParameter`](#)
types and [`Schedule`](#)s.

See [the tutorial](../docs/tutorials/hyperparameters.md).

## Example

```julia
lrschedule = Schedule([0, 10], [0.1, 0.001], Animations.sineio())
scheduler = Scheduler(
    LearningRate => lrschedule
)
```

See also [`Schedule`](#).
"""
struct Scheduler <: Callback
    schedules::Dict{Type{<:HyperParameter}, Animations.FiniteLengthAnimation}
    Scheduler(args...; unit = :epoch, kwargs...) = new(Dict(args...; kwargs...))
end

Base.show(io::IO, scheduler::Scheduler) = print(
    io, "Scheduler(", join(keys(scheduler.schedules), ", "), ")")

function stateaccess(scheduler::Scheduler)
    # TODO: implement proper merging of permissions
    if length(keys(scheduler.schedules)) == 0
        hpstateaccess = (;)
    else
        hpstateaccess = merge(stateaccess.(keys(scheduler.schedules))...)
    end
    return (
        data = Read(), cbstate = (history = Read(), hyperparams = Write()),
        hpstateaccess...
    )
end

function on(::Init, phase, scheduler::Scheduler, learner)
    learner.cbstate.hyperparams = MVHistory()
end

function on(::BatchBegin, ::AbstractTrainingPhase, scheduler::Scheduler, learner)
    step = learner.cbstate.history.steps + 1
    for (H, animation) in scheduler.schedules
        value = Animations.at(animation, step)
        sethyperparameter!(learner, H, value)
        push!(learner.cbstate.hyperparams, Symbol(H), step, value)
    end
end

"""
    onecycle(nepochs, max_val, [start_val, end_val; start_pctg])

Creates a one-cycle [`Schedule`](#) over `nepochs` epochs from `start_val`
over `max_val` to `end_val`.
"""
function onecycle(
        nepochs, max_val,
        start_val = max_val / 10,
        end_val = max_val / 30;
        start_pctg = 0.1)
    return Animation(
        [0, nepochs * start_pctg, nepochs],
        [start_val, max_val, end_val],
        [Animations.sineio(), Animations.sineio()]
    )

end


text:
__precompile__(true)

module Anamnesis

using MacroTools

using Core.Compiler: return_type

include("utils.jl")
include("scribe.jl")
include("macros.jl")

export Scribe
export call
export @scribeof, @mem, @localmem, @scribe, @localscribe, @rawcall, @rawfunc, @scribeofrawfunc, @anamnesis

end


text:
##
## Gamma Function
##

SpecialFunctions.gamma(z::Union{ArbOrRef,AcbOrRef}) = gamma!(zero(z), z)
SpecialFunctions.gamma(z::Series) = gamma_series!(zero(z), z, length(z))

SpecialFunctions.loggamma(x::Union{ArbOrRef,AcbOrRef}) = lgamma!(zero(x), x)
SpecialFunctions.loggamma(x::Series) = lgamma_series!(zero(x), x, length(x))

#SpecialFunctions.logabsgamma(x)
# Not implemented by Arb

#SpecialFunctions.logfactorial(x)
# Only relevant for Arblib

SpecialFunctions.digamma(x::Union{ArbOrRef,AcbOrRef}) = digamma!(zero(x), x)
SpecialFunctions.digamma(x::Series) = digamma_series!(zero(x), x, length(x))

#SpecialFunctions.invdigamma(x)
# Not implemented by Arb

function SpecialFunctions.trigamma(x::AcbOrRef)
    res = one(x)
    return polygamma!(res, res, x)
end

SpecialFunctions.polygamma(s::AcbOrRef, x::AcbOrRef) = polygamma!(zero(x), s, x)

SpecialFunctions.gamma(a::ArbOrRef, z::ArbOrRef) =
    hypgeom_gamma_upper!(Arb(0, prec = _precision((a, z))), a, z, 0)
SpecialFunctions.gamma(a::AcbOrRef, z::AcbOrRef) =
    hypgeom_gamma_upper!(Acb(0, prec = _precision((a, z))), a, z, 0)
SpecialFunctions.gamma(a::ArbOrRef, z::ArbSeries) =
    hypgeom_gamma_upper_series!(zero(z), a, z, 0, length(z))
SpecialFunctions.gamma(a::AcbOrRef, z::AcbSeries) =
    hypgeom_gamma_upper_series!(zero(z), a, z, 0, length(z))

#loggamma(a,z)
# Not implemented by Arb

function SpecialFunctions.gamma_inc(a::ArbOrRef, x::ArbOrRef)
    Γ = hypgeom_gamma_upper!(Arb(0, prec = _precision((a, x))), a, x, 1)
    # γ = 1 - Γ
    γ = neg!(Arb(0, prec = _precision((a, x))), Γ)
    add!(γ, γ, 1)
    return (γ, Γ)
end
function SpecialFunctions.gamma_inc(a::AcbOrRef, x::AcbOrRef)
    Γ = hypgeom_gamma_upper!(Acb(0, prec = _precision((a, x))), a, x, 1)
    # γ = 1 - Γ
    γ = neg!(Acb(0, prec = _precision((a, x))), Γ)
    add!(γ, γ, 1)
    return (γ, Γ)
end
function SpecialFunctions.gamma_inc(a::ArbOrRef, x::ArbSeries)
    Γ = hypgeom_gamma_upper_series!(zero(x), a, x, 1, length(x))
    γ = 1 - Γ
    return (γ, Γ)
end
function SpecialFunctions.gamma_inc(a::AcbOrRef, x::AcbSeries)
    Γ = hypgeom_gamma_upper_series!(zero(x), a, x, 1, length(x))
    γ = 1 - Γ
    return (γ, Γ)
end

function SpecialFunctions.beta_inc(a::ArbOrRef, b::ArbOrRef, x::ArbOrRef)
    β = hypgeom_beta_lower!(Arb(prec = _precision((a, x))), a, b, x, 1)
    # Β = 1 - β
    Β = neg!(Arb(0, prec = _precision((a, x))), β)
    add!(Β, Β, 1)
    return (β, Β)
end
function SpecialFunctions.beta_inc(a::AcbOrRef, b::AcbOrRef, x::AcbOrRef)
    β = hypgeom_beta_lower!(Acb(prec = _precision((a, x))), a, b, x, 1)
    # Β = 1 - β
    Β = neg!(Acb(0, prec = _precision((a, x))), β)
    add!(Β, Β, 1)
    return (β, Β)
end
function SpecialFunctions.beta_inc(a::ArbOrRef, b::ArbOrRef, x::ArbSeries)
    β = hypgeom_beta_lower_series!(zero(x), a, b, x, 1, length(x))
    Β = 1 - β
    return (β, Β)
end
function SpecialFunctions.beta_inc(a::AcbOrRef, b::AcbOrRef, x::AcbSeries)
    β = hypgeom_beta_lower_series!(zero(x), a, b, x, 1, length(x))
    Β = 1 - β
    return (β, Β)
end

#gamma_inc_inv(a, p, q)
# Not implemented by Arb

#beta(x,y)
# Not implemented directly by Arb, could use beta_inc or gamma to implement it?

#logbeta(x,y)
#Not implemented by Arb

#logabsbeta(x,y)
#Not implemented by Arb

#logabsbinomial(x,y)
#Not implemented by Arb

##
## Exponential and Trigonometric Integrals
##

# The version when ν is not specified could be defined directly, but
# it would likely be better if it was defined in SpecialFunctions
# directly.
SpecialFunctions.expint(ν::ArbOrRef, x::ArbOrRef) =
    hypgeom_expint!(Arb(prec = _precision((ν, x))), ν, x)
SpecialFunctions.expint(ν::AcbOrRef, x::AcbOrRef) =
    hypgeom_expint!(Acb(prec = _precision((ν, x))), ν, x)
SpecialFunctions.expint(ν::ArbOrRef, x::ArbSeries) =
    hypgeom_gamma_upper_series!(zero(x), 1 - ν, x, 2, length(x))
SpecialFunctions.expint(ν::AcbOrRef, x::AcbSeries) =
    hypgeom_gamma_upper_series!(zero(x), 1 - ν, x, 2, length(x))

function SpecialFunctions.expint(x::Union{ArbOrRef,AcbOrRef})
    ν = one(x)
    return hypgeom_expint!(ν, ν, x)
end
SpecialFunctions.expint(x::Series) =
    SpecialFunctions.expint(eltype(x)(1, prec = precision(x)), x)

SpecialFunctions.expinti(x::Union{ArbOrRef,AcbOrRef}) = hypgeom_ei!(zero(x), x)
SpecialFunctions.expinti(x::Series) = hypgeom_ei_series!(zero(x), x, length(x))

#expintx(x)
# Not implemented by arb

SpecialFunctions.sinint(x::Union{ArbOrRef,AcbOrRef}) = hypgeom!(zero(x), x)
SpecialFunctions.sinint(x::Series) = hypgeom_si_series!(zero(x), x, length(x))

SpecialFunctions.cosint(x::Union{ArbOrRef,AcbOrRef}) = hypgeom_ci!(zero(x), x)
SpecialFunctions.cosint(x::Series) = hypgeom_ci_series!(zero(x), x, length(x))

##
## Error Functions, Dawson’s and Fresnel Integrals
##

SpecialFunctions.erf(x::Union{ArbOrRef,AcbOrRef}) = hypgeom_erf!(zero(x), x)
SpecialFunctions.erf(x::Series) = hypgeom_erf_series!(zero(x), x, length(x))

#SpecialFunctions.erf(x,y)
# Not implemented by Arb

SpecialFunctions.erfc(x::Union{ArbOrRef,AcbOrRef}) = hypgeom_erfc!(zero(x), x)
SpecialFunctions.erfc(x::Series) = hypgeom_erfc_series!(zero(x), x, length(x))

#SpecialFunctions.erfcinv(x)
# Not implemented by Arb

#SpecialFunctions.erfcx(x)
# Not implemented by Arb

#SpecialFunctions.logerfc(x)
# Not implemented by Arb

#SpecialFunctions.logerfcx(x)
# Not implemented by Arb

#SpecialFunctions.erfi(x)
SpecialFunctions.erfi(x::Union{ArbOrRef,AcbOrRef}) = hypgeom_erfi!(zero(x), x)
SpecialFunctions.erfi(x::Series) = hypgeom_erfi_series!(zero(x), x, length(x))

#SpecialFunctions.erfinv(x)
# Not implemented by Arb

#SpecialFunctions.dawson(x)
# Not implemented by Arb

##
## Airy and Related Functions
##

function SpecialFunctions.airyai(z::ArbOrRef)
    ai = zero(z)
    ccall(
        @libarb(arb_hypgeom_airy),
        Cvoid,
        (Ref{arb_struct}, Ref{Cvoid}, Ref{Cvoid}, Ref{Cvoid}, Ref{arb_struct}, Int64),
        ai,
        C_NULL,
        C_NULL,
        C_NULL,
        z,
        precision(z),
    )
    return ai
end
function SpecialFunctions.airyai(z::AcbOrRef)
    ai = zero(z)
    ccall(
        @libarb(acb_hypgeom_airy),
        Cvoid,
        (Ref{acb_struct}, Ref{Cvoid}, Ref{Cvoid}, Ref{Cvoid}, Ref{acb_struct}, Int64),
        ai,
        C_NULL,
        C_NULL,
        C_NULL,
        z,
        precision(z),
    )
    return ai
end
function SpecialFunctions.airyai(z::ArbSeries)
    ai = zero(z)
    ccall(
        @libarb(arb_hypgeom_airy_series),
        Cvoid,
        (
            Ref{arb_poly_struct},
            Ref{Cvoid},
            Ref{Cvoid},
            Ref{Cvoid},
            Ref{arb_poly_struct},
            Int64,
            Int64,
        ),
        ai,
        C_NULL,
        C_NULL,
        C_NULL,
        z,
        length(z),
        precision(z),
    )
    return ai
end
function SpecialFunctions.airyai(z::AcbSeries)
    ai = zero(z)
    ccall(
        @libarb(acb_hypgeom_airy_series),
        Cvoid,
        (
            Ref{acb_poly_struct},
            Ref{Cvoid},
            Ref{Cvoid},
            Ref{Cvoid},
            Ref{acb_poly_struct},
            Int64,
            Int64,
        ),
        ai,
        C_NULL,
        C_NULL,
        C_NULL,
        z,
        length(z),
        precision(z),
    )
    return ai
end

function SpecialFunctions.airyaiprime(z::ArbOrRef)
    ai_prime = zero(z)
    ccall(
        @libarb(arb_hypgeom_airy),
        Cvoid,
        (Ref{Cvoid}, Ref{arb_struct}, Ref{Cvoid}, Ref{Cvoid}, Ref{arb_struct}, Int64),
        C_NULL,
        ai_prime,
        C_NULL,
        C_NULL,
        z,
        precision(z),
    )
    return ai_prime
end
function SpecialFunctions.airyaiprime(z::AcbOrRef)
    ai_prime = zero(z)
    ccall(
        @libarb(acb_hypgeom_airy),
        Cvoid,
        (Ref{Cvoid}, Ref{acb_struct}, Ref{Cvoid}, Ref{Cvoid}, Ref{acb_struct}, Int64),
        C_NULL,
        ai_prime,
        C_NULL,
        C_NULL,
        z,
        precision(z),
    )
    return ai_prime
end
function SpecialFunctions.airyaiprime(z::ArbSeries)
    ai_prime = zero(z)
    ccall(
        @libarb(arb_hypgeom_airy_series),
        Cvoid,
        (
            Ref{Cvoid},
            Ref{arb_poly_struct},
            Ref{Cvoid},
            Ref{Cvoid},
            Ref{arb_poly_struct},
            Int64,
            Int64,
        ),
        C_NULL,
        ai_prime,
        C_NULL,
        C_NULL,
        z,
        length(z),
        precision(z),
    )
    return ai_prime
end
function SpecialFunctions.airyaiprime(z::AcbSeries)
    ai_prime = zero(z)
    ccall(
        @libarb(acb_hypgeom_airy_series),
        Cvoid,
        (
            Ref{Cvoid},
            Ref{acb_poly_struct},
            Ref{Cvoid},
            Ref{Cvoid},
            Ref{acb_poly_struct},
            Int64,
            Int64,
        ),
        C_NULL,
        ai_prime,
        C_NULL,
        C_NULL,
        z,
        length(z),
        precision(z),
    )
    return ai_prime
end

function SpecialFunctions.airybi(z::ArbOrRef)
    bi = zero(z)
    ccall(
        @libarb(arb_hypgeom_airy),
        Cvoid,
        (Ref{Cvoid}, Ref{Cvoid}, Ref{arb_struct}, Ref{Cvoid}, Ref{arb_struct}, Int64),
        C_NULL,
        C_NULL,
        bi,
        C_NULL,
        z,
        precision(z),
    )
    return bi
end
function SpecialFunctions.airybi(z::AcbOrRef)
    bi = zero(z)
    ccall(
        @libarb(acb_hypgeom_airy),
        Cvoid,
        (Ref{Cvoid}, Ref{Cvoid}, Ref{acb_struct}, Ref{Cvoid}, Ref{acb_struct}, Int64),
        C_NULL,
        C_NULL,
        bi,
        C_NULL,
        z,
        precision(z),
    )
    return bi
end
function SpecialFunctions.airybi(z::ArbSeries)
    bi = zero(z)
    ccall(
        @libarb(arb_hypgeom_airy_series),
        Cvoid,
        (
            Ref{Cvoid},
            Ref{Cvoid},
            Ref{arb_poly_struct},
            Ref{Cvoid},
            Ref{arb_poly_struct},
            Int64,
            Int64,
        ),
        C_NULL,
        C_NULL,
        bi,
        C_NULL,
        z,
        length(z),
        precision(z),
    )
    return bi
end
function SpecialFunctions.airybi(z::AcbSeries)
    bi = zero(z)
    ccall(
        @libarb(acb_hypgeom_airy_series),
        Cvoid,
        (
            Ref{Cvoid},
            Ref{Cvoid},
            Ref{acb_poly_struct},
            Ref{Cvoid},
            Ref{acb_poly_struct},
            Int64,
            Int64,
        ),
        C_NULL,
        C_NULL,
        bi,
        C_NULL,
        z,
        length(z),
        precision(z),
    )
    return bi
end

function SpecialFunctions.airybiprime(z::ArbOrRef)
    bi_prime = zero(z)
    ccall(
        @libarb(arb_hypgeom_airy),
        Cvoid,
        (Ref{Cvoid}, Ref{Cvoid}, Ref{Cvoid}, Ref{arb_struct}, Ref{arb_struct}, Int64),
        C_NULL,
        C_NULL,
        C_NULL,
        bi_prime,
        z,
        precision(z),
    )
    return bi_prime
end
function SpecialFunctions.airybiprime(z::AcbOrRef)
    bi_prime = zero(z)
    ccall(
        @libarb(acb_hypgeom_airy),
        Cvoid,
        (Ref{Cvoid}, Ref{Cvoid}, Ref{Cvoid}, Ref{acb_struct}, Ref{acb_struct}, Int64),
        C_NULL,
        C_NULL,
        C_NULL,
        bi_prime,
        z,
        precision(z),
    )
    return bi_prime
end
function SpecialFunctions.airybiprime(z::ArbSeries)
    bi_prime = zero(z)
    ccall(
        @libarb(arb_hypgeom_airy_series),
        Cvoid,
        (
            Ref{Cvoid},
            Ref{Cvoid},
            Ref{Cvoid},
            Ref{arb_poly_struct},
            Ref{arb_poly_struct},
            Int64,
            Int64,
        ),
        C_NULL,
        C_NULL,
        C_NULL,
        bi_prime,
        z,
        length(z),
        precision(z),
    )
    return bi_prime
end
function SpecialFunctions.airybiprime(z::AcbSeries)
    bi_prime = zero(z)
    ccall(
        @libarb(acb_hypgeom_airy_series),
        Cvoid,
        (
            Ref{Cvoid},
            Ref{Cvoid},
            Ref{Cvoid},
            Ref{acb_poly_struct},
            Ref{acb_poly_struct},
            Int64,
            Int64,
        ),
        C_NULL,
        C_NULL,
        C_NULL,
        bi_prime,
        z,
        length(z),
        precision(z),
    )
    return bi_prime
end

#SpecialFunctions.airyaix(z)
# Not implemented by Arb

#SpecialFunctions.airyaiprimex(z)
# Not implemented by Arb

#SpecialFunctions.airybix(z)
# Not implemented by Arb

#SpecialFunctions.airybiprimex(z)
# Not implemented by Arb

##
## Bessel Functions
##

SpecialFunctions.besselj(ν::ArbOrRef, z::ArbOrRef) =
    hypgeom_bessel_j!(Arb(0, prec = _precision((ν, z))), ν, z)
SpecialFunctions.besselj(ν::AcbOrRef, z::AcbOrRef) =
    hypgeom_bessel_j!(Acb(0, prec = _precision((ν, z))), ν, z)

function SpecialFunctions.besselj0(z::Union{ArbOrRef,AcbOrRef})
    res = zero(z)
    return hypgeom_bessel_j!(res, res, z)
end

function SpecialFunctions.besselj1(z::Union{ArbOrRef,AcbOrRef})
    res = one(z)
    return hypgeom_bessel_j!(res, res, z)
end

#SpecialFunctions.besseljx(nu,z)
# Not implemented by Arb

#SpecialFunctions.sphericalbesselj(nu,z)
# The general method implemented by SpecialFunctions is not completely
# rigorous since it makes a cutoff for small values.
# TODO: We could check for the special case x = 0
function SpecialFunctions.sphericalbesselj(ν::ArbOrRef, x::ArbOrRef)
    # res = besselj(ν + 1 // 2, x)
    res = Arb(1 // 2, prec = _precision((ν, x)))
    add!(res, res, ν)
    hypgeom_bessel_j!(res, res, x)

    # factor = sqrt(π / 2x)
    factor = Arb(π, prec = _precision((ν, x)))
    div!(factor, factor, x)
    mul_2exp!(factor, factor, -1)
    sqrt!(factor, factor)

    return mul!(res, factor, res)
end

SpecialFunctions.bessely(ν::ArbOrRef, z::ArbOrRef) =
    hypgeom_bessel_y!(Arb(0, prec = _precision((ν, z))), ν, z)
SpecialFunctions.bessely(ν::AcbOrRef, z::AcbOrRef) =
    hypgeom_bessel_y!(Acb(0, prec = _precision((ν, z))), ν, z)

function SpecialFunctions.bessely0(z::Union{ArbOrRef,AcbOrRef})
    res = zero(z)
    return hypgeom_bessel_y!(res, res, z)
end

function SpecialFunctions.bessely1(z::Union{ArbOrRef,AcbOrRef})
    res = one(z)
    return hypgeom_bessel_y!(res, res, z)
end

#SpecialFunctions.besselyx(nu,z)
# Not implemented by Arb

#SpecialFunctions.sphericalbessely(nu,z)
# Aliased to √((float(T))(π)/2x) * bessely(nu + one(nu)/2, x) which works fine

SpecialFunctions.besselh(ν::ArbOrRef, k::Integer, z::ArbOrRef) =
    SpecialFunctions.besselh(Acb(ν), k, Acb(z))
function SpecialFunctions.besselh(ν::AcbOrRef, k::Integer, z::AcbOrRef)
    J, Y = Acb(prec = _precision((ν, z))), Acb(prec = _precision((ν, z)))
    hypgeom_bessel_jy!(J, Y, ν, z)
    mul_onei!(Y, Y)
    if k == 1
        return add!(J, J, Y)
    elseif k == 2
        return sub!(J, J, Y)
    else
        throw(SpecialFunctions.AmosException(1)) # This is what SpecialFunctions throw
    end
end

#SpecialFunctions.besselhx(nu,k,z)
# Not implemented by Arb

#SpecialFunctions.hankelh1(nu,z)
# Aliased to besselh(nu, 1, z)

#SpecialFunctions.hankelh1x(nu,z)
# Aliased to besselhx(nu, 1, z)

#SpecialFunctions.hankelh2(nu,z)
# Aliased to besselh(nu, 2, z)

#SpecialFunctions.hankelh2x(nu,z)
# Aliased to besselhx(nu, 2, z)

SpecialFunctions.besseli(ν::ArbOrRef, z::ArbOrRef) =
    hypgeom_bessel_i!(Arb(prec = _precision((ν, z))), ν, z)
SpecialFunctions.besseli(ν::AcbOrRef, z::AcbOrRef) =
    hypgeom_bessel_i!(Acb(prec = _precision((ν, z))), ν, z)

#SpecialFunctions.besselix(nu,z)
# The scaling used by Arb seems to be different from that used by
# SpecialFunctions.

SpecialFunctions.besselk(ν::ArbOrRef, z::ArbOrRef) =
    hypgeom_bessel_k!(Arb(prec = _precision((ν, z))), ν, z)
SpecialFunctions.besselk(ν::AcbOrRef, z::AcbOrRef) =
    hypgeom_bessel_k!(Acb(prec = _precision((ν, z))), ν, z)

SpecialFunctions.besselkx(ν::ArbOrRef, z::ArbOrRef) =
    hypgeom_bessel_k_scaled!(Arb(prec = _precision((ν, z))), ν, z)
SpecialFunctions.besselkx(ν::AcbOrRef, z::AcbOrRef) =
    hypgeom_bessel_k_scaled!(Acb(prec = _precision((ν, z))), ν, z)

#jinc(x)
# Aliased to 2 * besselj1(π*x) / (π*x) which works fine

##
## Elliptic Integrals
##

SpecialFunctions.ellipk(m::AcbOrRef) = elliptic_k!(zero(m), m)
SpecialFunctions.ellipk(m::AcbSeries) = elliptic_k_series!(zero(m), m, length(m))

SpecialFunctions.ellipe(m::AcbOrRef) = elliptic_e!(zero(m), m)

##
## Zeta and Related Functions
##

SpecialFunctions.eta(x::AcbOrRef) = dirichlet_eta!(zero(x), x)

SpecialFunctions.zeta(s::Union{ArbOrRef,AcbOrRef}) = zeta!(zero(s), s)
SpecialFunctions.zeta(s::Series) =
    SpecialFunctions.zeta(s, eltype(s)(1, prec = precision(s)))
SpecialFunctions.zeta(s::ArbOrRef, z::ArbOrRef) =
    hurwitz_zeta!(Arb(prec = _precision((s, z))), s, z)
SpecialFunctions.zeta(s::AcbOrRef, z::AcbOrRef) =
    hurwitz_zeta!(Acb(prec = _precision((s, z))), s, z)
SpecialFunctions.zeta(s::ArbSeries, z::ArbOrRef) = zeta_series!(zero(s), s, z, 0, length(s))
SpecialFunctions.zeta(s::AcbSeries, z::AcbOrRef) = zeta_series!(zero(s), s, z, 0, length(s))