text:
using Polynomials
using Plots
using PolynomialRoots
using ProgressMeter
using BenchmarkTools
using FileIO
using JLD2

gr()
#plotlyjs()
begin
	quot(v)=v[1]//v[2]
	vect(r)=[numerator(r),denominator(r)]
end

begin
	remainder(v) = v[1] % v[2]
	flr(v)=numerator(floor(v[1] // v[2]))
	newpair(v) = [v[2], remainder(v)]
	function continued_fraction(v)
		cf=[]
		while remainder(v) != 0
			push!(cf, flr(v))
			v=newpair(v)
		end
		push!(cf, flr(v))
		cf
	end
	function cf_value(cf)
		rev = reverse(cf)
		rat = popfirst!(rev)
		for j in rev
			rat = j+1//rat
		end
		[numerator(rat), denominator(rat)]
	end
end

function center_corner_index(v)
	v[1] == 1 && return [[0,1],[1,0],v[2]]
	center_cf = continued_fraction(v)[1:end-1]
	center = cf_value(center_cf)
	push!(center_cf, 1)
	k2 = cf_value(center_cf)
	corner = k2 - center
	ind = numerator((v[1] - corner[1]) // center[1])
	[center, corner, ind]
end

# testing updating q
function dQ(v)
		return Polynomial(push!(zeros(BigInt,v[2]), (-1)^v[1]))
end

function Q(v)
		v==[0,1] && return 1
		v==[1,0] && return 0
		v==[1,1] && return 1
		cci=center_corner_index(v)
		c,k,i=cci[1],cci[2],cci[3]
		return ([0 1; -dQ(c) Q(c)]^(i-1) * [Q(k); Q(k+c)])[2]
end

function bigroots(poly)
    PolynomialRoots.roots(convert(Vector{BigInt},coeffs(poly)))
end

function cleanpair(p)
    	r=p[1]//p[2]
    	return [numerator(r), denominator(r)]
end

function fractionline(n)  #fractions in (0,1/2) with denominator n.
    	top=ceil(n//2)
    	L=[cleanpair([j,n]) for j in 1:top-1 ]
    	return filter(rat-> rat[2]==n ,L)
end

function id(z)
	return z
end

function F64(z)
	return convert(Complex{Float64},z)
end

function generateDiscriminants(rational)
    q2 = Q(rational)^2
    dq  = dQ(rational)
    return (q2 - 4 * dq)
end

function generateQDict(max_denom)
	qDict = Dict{Vector{Int64}, Polynomial{BigInt}}([1,3] => Q([1,3]))

	for key in 4:max_denom
		rats = fractionline(key)
		for rat in rats
			qDict[rat] = Q(rat)
		end
	end
	return qDict
end

function generateDiscriminantsDict(qDict)
	discDict = Dict{Vector{Int64}, Polynomial{BigInt}}([1,3] => generateDiscriminants([1,3]))

	Threads.@threads for key in collect(keys(qDict))
		discDict[key] = generateDiscriminants(key)
	end
	return discDict
end

function generateRoots(dict)
	setprecision(BigFloat, 512) do
		n = length(dict)

		p = Progress(n, 1, "Generating roots...", 50)

		keysdict = collect(keys(dict))
		rootsDict = Dict(keysdict[1] => PolynomialRoots.roots(keysdict[1]))

		@Threads.threads for key in keysdict
			poly = coeffs(dict[key])
			rootsDict[key] = PolynomialRoots.roots(poly)
			next!(p)
		end
	return rootsDict
	end
end

function generateRootset(rootsdict)
	rootset = []

	denKey = collect(keys(rootsdict))
	n = length(denKey) +2

	p = Progress(n, 1, "Mapping unions...", 50)

	Threads.@threads for key in denKey
        	rootset=union(rootset, map(id,rootsdict[key]))
		next!(p)
	end
	return rootset
end

function runQPlot(rootset, discset, max_denom)

	c1=filter(z->abs2(z)<3, rootset);
	c2=filter(z->abs2(z)<3, discset);

	scatter([c1[j] for j in 1:length(c1)], markersize = 1,markerstrokewidth=0, c = :black, size = (5000,5000), label=false, aspect_ratio=1, framestyle= :none, background_color= :ivory)
	scatter!([c2[j] for j in 1:length(c1)], markersize = 1,markerstrokewidth=0, c = :red)

	savefig("plots/scatter_" * string(max_denom) * ".svg")
end

function main(max_denom)
	if isfile("data/roots_" * string(max_denom) * ".jld2")
		r = load("data/roots_" * string(max_denom) * ".jld2")
		rs = load("data/disc_roots_" * string(max_denom) * ".jld2")
		r = generateRootset(r)
		rs = generateRootset(rs)

		runQPlot(r,rs,max_denom)

	else
	q = generateQDict(max_denom)
	d = generateDiscriminantsDict(q)

	r = generateRoots(q)
	save(r, "data/roots_" * string(max_denom) * ".jld2")

	rs = generateRoots(d)
	save(rs, "data/disc_roots_" * string(max_denom) * ".jld2")

	r = generateRootset(r)
	rs = generateRootset(rs)

	runQPlot(r,rs,max_denom)
	end
end

main(200)

max_stars_repo_name: archmagethanos/turnip

max_issues_repo_name: archmagethanos/turnip

max_forks_repo_name: archmagethanos/turnip


text:
module SPADL

    """
        SPADL
    A Julia package for Event Stream Data (Wyscout) translation to the
    Soccer Player Action Description Language (SPADL).
    """

    using DataFrames
    using HTTP
    using ZipFile
    using JSON3
    using JSONTables
    using ProgressMeter
    using HDF5
    using LazyJSON

    export
    # Types
    WyscoutEvent,
    ActionsConfig,

    # Serializers
    get_events_data,
    create_match_index,
    get_matchs,
    get_df,
    competitions,
    games,
    teams,
    players,
    get_events,
    lineup,
    events,

    # SPADL
    get_tags,
    # make_new_positions,
    convert_to_actions,
    convert_duels,
    # fix_events,
    selector_duel_out_of_field,

    # utils
    to_hdf

    include("data/data.jl")
    # include("utils.jl")
    include("spadl/spadl.jl")

    function __init__()

        nothing
    end

end

ext: jl


text:
struct GraphicsPipelineCreateInfo
    mHandleRef::Ref{vk.VkGraphicsPipelineCreateInfo}
    mPreserve::Vector{Any}
    mLayout

    function GraphicsPipelineCreateInfo(;
        pNext::Ptr{Cvoid} = C_NULL,
        flags::vk.VkPipelineCreateFlags = vk.VkFlags(0),
        stages::Vector{vk.VkPipelineShaderStageCreateInfo} = Vector{vk.VkPipelineShaderStageCreateInfo}(),
        vertexInputState = nothing, #vk.VkPipelineVertexInputStateCreateInfo,
        inputAssemblyState = nothing, #vk.VkPipelineInputAssemblyStateCreateInfo,
        tessellationState = nothing, #vk.VkPipelineTessellationStateCreateInfo,
        viewportState = nothing, #vk.VkPipelineViewportStateCreateInfo,
        rasterizationState = nothing, #vk.VkPipelineRasterizationStateCreateInfo,
        multisampleState = nothing, #vk.VkPipelineMultisampleStateCreateInfo,
        depthStencilState = nothing, #vk.VkPipelineDepthStencilStateCreateInfo,
        colorBlendState = nothing, #vk.VkPipelineColorBlendStateCreateInfo,
        dynamicState = nothing, #vk.VkPipelineDynamicStateCreateInfo,
        layout, #required
        renderPass::vk.VkRenderPass, #required
        subpass::UInt32 = UInt32(0),
        basePipelineHandle::vk.VkPipeline = vk.VkPipeline(vk.VK_NULL_HANDLE),
        basePipelineIndex::Int32 = Int32(0),
    )
        this = new(Ref(vk.VkGraphicsPipelineCreateInfo(
            vk.VK_STRUCTURE_TYPE_GRAPHICS_PIPELINE_CREATE_INFO, #sType::VkStructureType
            pNext, #pNext::Ptr{Cvoid}
            flags, #flags::VkPipelineCreateFlags
            length(stages), #stageCount::UInt32
            pointer(stages), #pStages::Ptr{VkPipelineShaderStageCreateInfo}
            ref_to_pointer(vk.VkPipelineVertexInputStateCreateInfo, handleRef(vertexInputState)), #::Ptr{VkPipelineVertexInputStateCreateInfo}
            ref_to_pointer(vk.VkPipelineInputAssemblyStateCreateInfo, handleRef(inputAssemblyState)), #::Ptr{VkPipelineInputAssemblyStateCreateInfo}
            ref_to_pointer(vk.VkPipelineTessellationStateCreateInfo, handleRef(tessellationState)), #::Ptr{VkPipelineTessellationStateCreateInfo}
            ref_to_pointer(vk.VkPipelineViewportStateCreateInfo, handleRef(viewportState)), #::Ptr{VkPipelineViewportStateCreateInfo}
            ref_to_pointer(vk.VkPipelineRasterizationStateCreateInfo, handleRef(rasterizationState)), #::Ptr{VkPipelineRasterizationStateCreateInfo}
            ref_to_pointer(vk.VkPipelineMultisampleStateCreateInfo, handleRef(multisampleState)), #::Ptr{VkPipelineMultisampleStateCreateInfo}
            ref_to_pointer(vk.VkPipelineDepthStencilStateCreateInfo, handleRef(depthStencilState)), #::Ptr{VkPipelineDepthStencilStateCreateInfo}
            ref_to_pointer(vk.VkPipelineColorBlendStateCreateInfo, handleRef(colorBlendState)), #::Ptr{VkPipelineColorBlendStateCreateInfo}
            ref_to_pointer(vk.VkPipelineDynamicStateCreateInfo, handleRef(dynamicState)), #::Ptr{VkPipelineDynamicStateCreateInfo}
            handleRef(layout)[], #::VkPipelineLayout
            renderPass, #::VkRenderPass
            subpass, #::UInt32
            basePipelineHandle, #::VkPipeline
            basePipelineIndex #::Int32
            )),
            [stages, vertexInputState, inputAssemblyState, tessellationState,
                 viewportState, rasterizationState, multisampleState,
                 depthStencilState, colorBlendState, dynamicState], #mPreserve
            layout # mLayout
        )
        return this
    end
end

function handleRef(this::GraphicsPipelineCreateInfo)::Ref{vk.VkGraphicsPipelineCreateInfo}
    return this.mHandleRef
end

function defaults(::Type{GraphicsPipelineCreateInfo};
    stages::Vector{vk.VkPipelineShaderStageCreateInfo} = Vector{vk.VkPipelineShaderStageCreateInfo}(),
    layout, #required
    renderPass::vk.VkRenderPass, #required
    subpass::UInt32 = UInt32(0),
    depthTestEnable::vk.VkBool32 = VkExt.VK_FALSE,
    depthWriteEnable::vk.VkBool32 = VkExt.VK_FALSE,
    frontFace::vk.VkFrontFace = vk.VK_FRONT_FACE_COUNTER_CLOCKWISE,
    inputState::PipelineVertexInputStateCreateInfo = PipelineVertexInputStateCreateInfo()
)::GraphicsPipelineCreateInfo
    inputAssemblyState = PipelineInputAssemblyStateCreateInfo(
                                    topology = vk.VK_PRIMITIVE_TOPOLOGY_TRIANGLE_LIST, #::VkPrimitiveTopology
                                )
    tessellationState = PipelineTessellationStateCreateInfo()

    viewports = [vk.VkViewport(
                    0, #x::Cfloat
                    0, #y::Cfloat
                    float(typemax(Int32)), #width::Cfloat
                    float(typemax(Int32)), #height::Cfloat
                    0.0, #minDepth::Cfloat
                    1.0 #maxDepth::Cfloat
                )]
    scissors = [vk.VkRect2D(vk.VkOffset2D(0, 0), vk.VkExtent2D(typemax(Int32), typemax(Int32)))]
    viewportState = PipelineViewportStateCreateInfo(
                        viewports = viewports,
                        scissors = scissors
                    )

    rasterizationState = PipelineRasterizationStateCreateInfo(
                            polygonMode = vk.VK_POLYGON_MODE_FILL, #::VkPolygonMode
                            frontFace = frontFace, #::VkFrontFace
                            lineWidth = Cfloat(1.0) #::Cfloat
                        )

    depthStencilState = PipelineDepthStencilStateCreateInfo(
                            depthTestEnable = depthTestEnable,
                            depthWriteEnable = depthWriteEnable,
                            depthCompareOp = vk.VK_COMPARE_OP_GREATER, #::VkCompareOp
                        )

    multisampleState = PipelineMultisampleStateCreateInfo(
                            rasterizationSamples = vk.VK_SAMPLE_COUNT_1_BIT, #::VkSampleCountFlagBits
                            minSampleShading = Cfloat(1.0) #::Cfloat
                        )
    #addNoBlend
    attachments = [PipelineColorBlendAttachmentState(
                        colorBlendOp = vk.VK_BLEND_OP_ADD, #::VkBlendOp
                        alphaBlendOp = vk.VK_BLEND_OP_ADD, #::VkBlendOp
                        colorWriteMask = (vk.VK_COLOR_COMPONENT_R_BIT | vk.VK_COLOR_COMPONENT_G_BIT |
                                            vk.VK_COLOR_COMPONENT_B_BIT | vk.VK_COLOR_COMPONENT_A_BIT) #::VkColorComponentFlags
                    ).handle()]
    colorBlendState = PipelineColorBlendStateCreateInfo(
                            logicOp = vk.VK_LOGIC_OP_COPY, #::VkLogicOp
                            attachments = attachments #::Vector{vk.VkPipelineColorBlendAttachmentState}
                        )

    dynamicState = PipelineDynamicStateCreateInfo(states = [vk.VK_DYNAMIC_STATE_VIEWPORT])

    info = GraphicsPipelineCreateInfo(;
        stages = stages, #::Vector{vk.VkPipelineShaderStageCreateInfo
        vertexInputState = inputState, #::VkPipelineVertexInputStateCreateInfo
        inputAssemblyState = inputAssemblyState, #::VkPipelineInputAssemblyStateCreateInfo
        tessellationState = tessellationState, #::VkPipelineTessellationStateCreateInfo
        viewportState = viewportState, #::VkPipelineViewportStateCreateInfo
        rasterizationState = rasterizationState, #::VkPipelineRasterizationStateCreateInfo
        multisampleState = multisampleState, #::VkPipelineMultisampleStateCreateInfo
        depthStencilState = depthStencilState, #::VkPipelineDepthStencilStateCreateInfo
        colorBlendState = colorBlendState, #::VkPipelineColorBlendStateCreateInfo
        dynamicState = dynamicState, #
        layout = layout, #::PipelineLayout
        renderPass = renderPass, #::VkRenderPass
        subpass = subpass #::UInt32
    )
    return info
end

function deriveFrom(outFlagRef::Ref{vk.VkPipelineCreateFlags})
    outFlagRef[] |= vk.VK_PIPELINE_CREATE_DERIVATIVE_BIT
end

function allowDerivatives(outFlagRef::Ref{vk.VkPipelineCreateFlags}, val::Bool)
    if (val)
        outFlagRef[] |= vk.VK_PIPELINE_CREATE_ALLOW_DERIVATIVES_BIT
    else
        outFlagRef[] &= ~UInt32(vk.VK_PIPELINE_CREATE_ALLOW_DERIVATIVES_BIT)
    end
end

function disableOptimization(outFlagRef::Ref{vk.VkPipelineCreateFlags}, val::Bool)
    if (val)
        outFlagRef[] |= vk.VK_PIPELINE_CREATE_DISABLE_OPTIMIZATION_BIT
    else
        outFlagRef[] &= ~UInt32(vk.VK_PIPELINE_CREATE_DISABLE_OPTIMIZATION_BIT)
    end
end


text:
using SparseTransforms
using Test
include("test_signal.jl")
include("test_query.jl")
include("test_reconstruct.jl")
include("test_spright.jl")
include("test_random.jl")
include("test_codes.jl")

function tests()
    println("Testing SPRIGHT")
    test_signal()
    test_query()
    test_reconstruct()
    test_spright()
    test_random()
    test_codes()
end

tests()


text:
using BinaryBuilder, Pkg

name = "XGBoost"
version = v"1.6.1"

# Collection of sources required to build XGBoost
sources = [
    GitSource("https://github.com/dmlc/xgboost.git","5d92a7d936fc3fad4c7ecb6031c3c1c7da882a14"),
    DirectorySource("./bundled"),
    ArchiveSource("https://github.com/phracker/MacOSX-SDKs/releases/download/10.15/MacOSX10.15.sdk.tar.xz",
                  "2408d07df7f324d3beea818585a6d990ba99587c218a3969f924dfcc4de93b62")
]

# Bash recipe for building across all platforms
script = raw"""
cd ${WORKSPACE}/srcdir/xgboost
git submodule update --init

# Patch dmlc-core to use case-sensitive windows.h includes
(cd dmlc-core; atomic_patch -p1 "../../patches/dmlc_windows.patch")

if [[ "${target}" == x86_64-apple-darwin* ]]; then
    export MACOSX_DEPLOYMENT_TARGET=10.15
    #install a newer SDK which supports `___cxa_deleted_virtual`
    pushd $WORKSPACE/srcdir/MacOSX10.*.sdk
    rm -rf /opt/${target}/${target}/sys-root/System
    cp -ra usr/* "/opt/${target}/${target}/sys-root/usr/."
    cp -ra System "/opt/${target}/${target}/sys-root/."
    popd
fi

mkdir build && cd build
cmake .. -DCMAKE_INSTALL_PREFIX=${prefix} -DCMAKE_TOOLCHAIN_FILE="${CMAKE_TARGET_TOOLCHAIN}"
make -j${nproc}

# Manual installation, to avoid installing dmlc
cd ..
for header in include/xgboost/*.h; do
    install -Dv "${header}" "${includedir}/xgboost/$(basename ${header})"
done
install -Dvm 0755 xgboost${exeext} ${bindir}/xgboost${exeext}

if [[ ${target} == *mingw* ]]; then
    install -Dvm 0755 lib/xgboost.dll ${libdir}/xgboost.dll
else
    install -Dvm 0755 lib/libxgboost.${dlext} ${libdir}/libxgboost.${dlext}
fi

install_license LICENSE
"""

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

# The products that we will ensure are always built
products = [
    LibraryProduct(["libxgboost", "xgboost"], :libxgboost),
    ExecutableProduct("xgboost", :xgboost)
]

# Dependencies that must be installed before this package can be built
dependencies = [
    # For OpenMP we use libomp from `LLVMOpenMP_jll` where we use LLVM as compiler (BSD
    # systems), and libgomp from `CompilerSupportLibraries_jll` everywhere else.
    Dependency(PackageSpec(name="CompilerSupportLibraries_jll", uuid="e66e0078-7015-5450-92f7-15fbd957f2ae"); platforms=filter(!Sys.isbsd, platforms)),
    Dependency(PackageSpec(name="LLVMOpenMP_jll", uuid="1d63c593-3942-5779-bab2-d838dc0a180e"); platforms=filter(Sys.isbsd, platforms)),
]

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


text:
module  TestIterators

using PopGenCore
using Test

x = @nancycats

@testset "Iterators.jl" begin
    @testset "partition array" begin
        @test partitionarray(rand(20,5), [10,3,4,3]) .|> size == [(10, 5), (3, 5), (4, 5), (3, 5)]
        @test partitionarray(rand(20), [10,3,4,3]) .|> size == [(10, 1), (3, 1), (4, 1), (3, 1)]
    end

    @testset "pairiwse pairs" begin
        @test length(collect(pairwisepairs(1:3))) == 3
        @test collect(pairwisepairs(1:3)) == [(1, 2), (1, 3), (2, 3)]
        @test length(collect(pairwisepairs(1:10))) == 45
        @test length(collect(pairwisepairs(string.(1:3)))) == 3
        @test collect(pairwisepairs(string.(1:3))) == [("1", "2"), ("1", "3"), ("2", "3")]
        @test length(collect(pairwisepairs(string.(1:10)))) == 45
    end

    #TODO moved to PopGenSims.jl
    @testset "simulated pairs" begin
        @test simpairs(string.(1:4)) == [("1","2"), ("3","4")]
        @test_throws ArgumentError simpairs(string.(1:3))
    end

    @testset "skip iterators" begin
        @test collect(skipinf([1,2,3, Inf, -Inf, 4])) == [1,2,3,4]
        @test collect(skipnan([1,2,3, NaN,4])) == [1,2,3,4]
        @test collect(skipinfnan([1,2,3, Inf, -Inf, NaN, 4])) == [1,2,3,4]
    end
end

end


text:
using CSV, Plots
using DiseaseOutbreak

# Hubei and Shandong
hubeiFile = joinpath(@__DIR__, "..", "data",
     "covid19_china","time_series_covid19_confirmed_hubei.csv")
shandongFile = joinpath(@__DIR__, "..", "data",
          "covid19_china","time_series_covid19_confirmed_shandong.csv")

hubei = CSV.read(hubeiFile)
shandong = CSV.read(shandongFile)
hubeiPop = 57.0e6
shandongPop = 94.2e6

hubeiC = convert(Vector{Float64},hubei.ConfirmedCases[1:22])
hubeiFit = fitCaseModel(22,hubeiC,hubeiPop,
                        6.2,8.0,[0.1,0.1,2.0])
summary(hubeiFit)
plotFit(hubeiFit,40)
# savefig("hubei.pdf")
pyplotFit(hubeiFit, 40, log10="semilog", grid=true, fsize=[12, 8])
pyplotFit(hubeiFit, 40, log10="loglog", grid=true, fsize=[12, 8])
pyplotFit(hubeiFit, 40, log10="null", grid=true, fsize=[12, 8])

shandongC = convert(Vector{Float64},shandong.ConfirmedCases[1:30])
shandongFit = fitCaseModel(30,shandongC,shandongPop,
                           6.2,8.0,[0.5,0.045,15.0])
summary(shandongFit)

plotFit(shandongFit,40)
plot!(estimatedStates(30,shandongPop,shandongC[1],9.66,
                getParams(0.309,0.042,6.2,8.0,SIRX()))[!,:X],
                yaxis=:log,label="paper")

pyplotFit(shandongFit, 40, log10="semilog", grid=true, fsize=[12, 8])
paperRef = estimatedStates(30,shandongPop,shandongC[1],9.66,
                getParams(0.309,0.042,6.2,8.0,SIRX()))[!,:X]

pyplotFit(shandongFit, 30, ref=paperRef, log10="semilog", grid=true,
          fsize=[12, 8])

fitted(shandongFit,40)


text:
# Dummy test of sphere fitting

using RANSAC: FittedPlane, FittedSphere, fit
#using RANSAC: RANSACParameters

const EPSI = 0.1
const ALFI = deg2rad(10)
rp = ransacparameters()
#const defrp = RANSACParameters(rp, ϵ_sphere = EPSI, α_sphere = ALFI)
const defrp = ransacparameters(rp; sphere=(ϵ=EPSI, α=ALFI,))
rps(n) = [SVector{3,Float64}(rand(3)) for i in 1:n]
const unnecrc = RANSACCloud(rps(10), rps(10), 1)

@testset "true sphere 1" begin
    tv1 = [SVector(0,-1,0.0), SVector(0,0,-1.0), SVector(1,0,0.0), SVector(0,1,0.0)]
    tn1 = [SVector(0,-1,0.0), SVector(0,0,-1.0), SVector(1,0,0.0), SVector(0,1,0.0)]
    fs = fit(FittedSphere, tv1, tn1, unnecrc, defrp)
    #fp = fitplane(tv1, tn1, RANSACParameters(defrp, α_plane=π/2, collin_threshold=0.2))
    fp = fit(FittedPlane, tv1, tn1, unnecrc,  ransacparameters(defrp, plane=(α=π/2,), common=(collin_threshold=0.2,)))
    @test fs isa FittedSphere
    @test fp === nothing
end

@testset "true sphere 2" begin
    tv1 = [SVector(0,-0.99,0.0), SVector(0,0,-1.0), SVector(1.01,0,0.0), SVector(0,1,0.0)]
    tn1 = [SVector(0,-1,0.0), SVector(0,0,-1.0), SVector(1,0,0.0), SVector(0,1,0.0)]
    fs1 = fit(FittedSphere, tv1, tn1, unnecrc, defrp) # true
    #fs2 = fitsphere(tv1, tn1, RANSACParameters(defrp, ϵ_sphere=0.01)) # false, cause ϵ
    fs2 = fit(FittedSphere, tv1, tn1, unnecrc, ransacparameters(defrp, sphere=(ϵ=0.01,))) # false, cause ϵ
    #fp = fitplane(tv1, tn1, RANSACParameters(defrp, α_plane=π/2, collin_threshold=0.2))
    fp = fit(FittedPlane, tv1, tn1, unnecrc, ransacparameters(defrp, plane=(α=π/2,), common=(collin_threshold=0.2,)))
    @test fs1 isa FittedSphere
    @test fs2 === nothing
    @test fp === nothing
end

@testset "false sphere 1" begin
    tv1 = [SVector(0,1,0.0), SVector(0,0,-1.0), SVector(1,0,0.0), SVector(0,1,0.0)]
    tn1 = [SVector(0,-1,0.0), SVector(0,0,-1.0), SVector(1,0,0.0), SVector(0,1,0.0)]
    fs1 = fit(FittedSphere, tv1, tn1, unnecrc, defrp)
    # should be false even with large thresholds
    #fs2 = fitsphere(tv1, tn1, RANSACParameters(defrp, ϵ_sphere=10, α_sphere=π/2))
    fs2 = fit(FittedSphere, tv1, tn1, unnecrc, ransacparameters(defrp, sphere=(ϵ=10,α=π/2,)))
    #fp = fitplane(tv1, tn1, RANSACParameters(defrp, α_plane=π/2, collin_threshold=0.2))
    fp = fit(FittedPlane, tv1, tn1, unnecrc, ransacparameters(defrp, plane=(α=π/2,), common=(collin_threshold=0.2,)))
    @test fs1 === nothing
    @test fs2 === nothing
    @test fp === nothing
end


text:
const rot = [0 -1; 1 0]

export step!

function step!(sq, pc)

    dtmin, fant, i1, i2 = dt_min_position(pc)

    for i in 1:sq.n
        sq.q[i] = mod.(sq.q[i] .+ dtmin .* sq.v[i], 1)
    end

    if sq.v[i1]'sq.v[i2] == 0

        sq.v[[i1,i2]] = sq.v[[i2,i1]] # swap velocities

    elseif sq.v[i1]'sq.v[i2] == -1

        if (rot * sq.v[i1])'*(sq.q[i2] .+ offset[fant] .- sq.q[i1]) < 0
            sq.v[i1] = rot * sq.v[i1]
            sq.v[i2] = rot * sq.v[i2]
        elseif (rot * sq.v[i1])'*(sq.q[i2] .+ offset[fant] .- sq.q[i1]) > 0
            sq.v[i1] = - rot * sq.v[i1]
            sq.v[i2] = - rot * sq.v[i2]
        else # very rare case
            sq.v[[i1,i2]] = sq.v[[i2,i1]] # swap velocities
        end

    end

    pc.dt .-= dtmin
    reset!(pc, i1)
    reset!(pc, i2)

    compute_dt!(pc, i1, sq)
    compute_dt!(pc, i2, sq)

    return dtmin

end


text:
using PyCall, PyPlot, Random
Random.seed!(5)
sp = pyimport("scipy.interpolate")
x = LinRange(0,1,50)
y = (0.9 .+ 0.1rand(length(x))).*sin.(2*pi*(x.-0.5))
t = collect(x[2:2:end-1])
s1 = sp.LSQUnivariateSpline(x, y, t)
x2 = LinRange(0, 1, 201) # new x-grid
y2 = s1(x2) # evaluate spline on that new grid
figure()
plot(x,y,label="original")
plot(x2,y2,label="interp", color="k")
knots = s1.get_knots()
c = s1.get_coeffs()
newknots(knots, k) = vcat(fill(knots[1],k),knots,fill(knots[end],k))
forscipyknots = newknots(knots, 3)
s2 = sp.BSpline(forscipyknots, c, 3)
y3 = s2(x2)
plot(x2,y3,"--r", label="reconstructed \nfrom knots and coeff")
legend()
##
knots = [0,.4,.4,.4,.4,.7,1]
c = [2,-5,5,2,-3,-1,2,1.5]
forscipyknots = newknots(knots, 3)
s2 = sp.BSpline(forscipyknots, c, 3)
figure()
plot(x2, s2(x2))
##
using Images
img = convert(Array{Float64, 2}, Gray.(load("func.png")))


text:
struct Static{T} <: LineSearcher
    α::T
end

find_steplength(mstyle, ls::Static, φ, λ::T) where {T} = T(ls.α), φ(T(ls.α)), true

max_stars_repo_name: JuliaNLSolvers/NLSolvers

max_issues_repo_name: JuliaNLSolvers/NLSolvers

max_forks_repo_name: JuliaNLSolvers/NLSolvers


text:
function count_arrivals(limit::Int)
    destinations = destination_lookup(81 * Int(ceil(log10(limit))))

    count = 0
    for n in 1:limit
        count += destinations[digit_square(n)]
    end
    println(count)
end

function destination_lookup(size::Int)
    destinations = Dict()
    for n in 1:size
        destinations[n] = final_destination(n)
    end
    return destinations
end

function final_destination(n::Int)
    while n != 1 && n != 89
        n = digit_square(n)
    end

    if n == 1
        return 0
    end

    return 1
end

function digit_square(n::Int)
    return sum(d^2 for d in digits(n))
end

count_arrivals(10000000)

max_stars_repo_name: HarrisonGreen/Project-Euler-Solutions

max_issues_repo_name: HarrisonGreen/Project-Euler-Solutions

max_forks_repo_name: HarrisonGreen/Project-Euler-Solutions


text:
abstract type QuantumVariable{T} <: Convex.AbstractVariable{T} end

export ProbabilityVector, DensityMatrix, Choi

using Convex: ⪰

include("probability_vectors.jl")
include("choi.jl")
include("density_matrices.jl")


text:
"""
    artin_braid(f::FusionTree, i; inv::Bool = false) -> <:AbstractDict{typeof(t),<:Number}

Perform an elementary braid (Artin generator) of neighbouring uncoupled indices `i` and
`i+1` on a fusion tree `f`, and returns the result as a dictionary of output trees and
corresponding coefficients.

The keyword `inv` determines whether index `i` will braid above or below index `i+1`, i.e.
applying `artin_braid(f′, i; inv = true)` to all the outputs `f′` of
`artin_braid(f, i; inv = false)` and collecting the results should yield a single fusion
tree with non-zero coefficient, namely `f` with coefficient `1`. This keyword has no effect
if  `BraidingStyle(sectortype(f)) isa SymmetricBraiding`.
"""
function artin_braid(f::FusionTree{G,N}, i; inv::Bool = false) where {G<:Sector, N}
    1 <= i < N ||
        throw(ArgumentError("Cannot swap outputs i=$i and i+1 out of only $N outputs"))
    uncoupled = f.uncoupled
    coupled′ = f.coupled
    isdual′ = TupleTools.setindex(f.isdual, f.isdual[i], i+1)
    isdual′ = TupleTools.setindex(isdual′, f.isdual[i+1], i)
    inner = f.innerlines
    vertices = f.vertices
    if i == 1
        a, b = uncoupled[1], uncoupled[2]
        c = N > 2 ? inner[1] : coupled′
        uncoupled′ = TupleTools.setindex(uncoupled, b, 1)
        uncoupled′ = TupleTools.setindex(uncoupled′, a, 2)
        R = inv ? conj(Rsymbol(b,a,c)) : Rsymbol(a,b,c)
        f′ = FusionTree{G}(uncoupled′, coupled′, isdual′, inner, vertices)
        if FusionStyle(G) isa Abelian
            return SingletonDict(f′ => R)
        elseif FusionStyle(G) isa SimpleNonAbelian
            return FusionTreeDict(f′ => R)
        end
    end
    # case i > 1:
    b = uncoupled[i]
    d = uncoupled[i+1]
    a = i == 2 ? uncoupled[1] : inner[i-2]
    c = inner[i-1]
    e = i == N-1 ? coupled′ : inner[i]
    uncoupled′ = TupleTools.setindex(uncoupled, d, i)
    uncoupled′ = TupleTools.setindex(uncoupled′, b, i+1)
    if FusionStyle(G) isa Abelian
        inner′ = TupleTools.setindex(inner, first(a ⊗ d), i-1)
        R = inv ? conj(Rsymbol(d, b, first(b ⊗ d))) : Rsymbol(b, d, first(b ⊗ d))
        return SingletonDict(FusionTree{G}(uncoupled′, coupled′, isdual′, inner′) => R)
    elseif FusionStyle(G) isa SimpleNonAbelian
        local newtrees
        for c′ in a ⊗ d
            coeff = if inv
                    Rsymbol(a,b,c)*Fsymbol(b,a,d,e,c,c′)*conj(Rsymbol(c′,b,e))
                else
                    conj(Rsymbol(b,a,c))*Fsymbol(b,a,d,e,c,c′)*Rsymbol(b,c′,e)
                end
            inner′ = TupleTools.setindex(inner, c′, i-1)
            f′ = FusionTree{G}(uncoupled′, coupled′, isdual′, inner′)
            if coeff != zero(coeff)
                if @isdefined newtrees
                    newtrees[f′] = coeff
                else
                    newtrees = FusionTreeDict(f′ => coeff)
                end
            end
        end
        return newtrees
    else
        # TODO: implement DegenerateNonAbelian case
        throw(MethodError(artin_braid, (f, i)))
    end
end

# braid fusion tree
"""
    braid(f::FusionTree{<:Sector,N}, levels::NTuple{N,Int}, p::NTuple{N,Int})
    -> <:AbstractDict{typeof(t),<:Number}

Perform a braiding of the uncoupled indices of the fusion tree `f` and returns the result
as a `<:AbstractDict` of output trees and corresponding coefficients. The braiding is
specified by specifying that index `i` goes to position `perm[i]` and assinging to every
index a distinct level or depth `levels[i]`. This permutation is then decomposed into
elementary swaps between neighbouring indices, where the swaps are applied as braids such
that if `i` and `j` cross, ``τ_{i,j}`` is applied if `levels[i] < levels[j]` and
``τ_{j,i}^{-1}`` if `levels[i] > levels[j]`. This does not allow to encode the most general
braid, but a general braid can be obtained by combining such operations.
"""
function braid(f::FusionTree{G,N},
                levels::NTuple{N,Int},
                p::NTuple{N,Int}) where {G<:Sector, N}
    TupleTools.isperm(p) || throw(ArgumentError("not a valid permutation: $p"))
    if FusionStyle(G) isa Abelian && BraidingStyle(G) isa SymmetricBraiding
        coeff = Rsymbol(one(G), one(G), one(G))
        for i = 1:N
            for j = 1:i-1
                if p[j] > p[i]
                    a, b = f.uncoupled[j], f.uncoupled[i]
                    coeff *= Rsymbol(a, b, first(a ⊗ b))
                end
            end
        end
        uncoupled′ = TupleTools._permute(f.uncoupled, p)
        coupled′ = f.coupled
        isdual′ = TupleTools._permute(f.isdual, p)
        f′ = FusionTree{G}(uncoupled′, coupled′, isdual′)
        return SingletonDict(f′=>coeff)
    else
        coeff = Rsymbol(one(G), one(G), one(G))
        trees = FusionTreeDict(f=>coeff)
        newtrees = empty(trees)
        for s in permutation2swaps(p)
            inv = levels[s] > levels[s+1]
            for (f, c) in trees
                for (f′,c′) in artin_braid(f, s; inv = inv)
                    newtrees[f′] = get(newtrees, f′, zero(coeff)) + c*c′
                end
            end
            l = levels[s]
            levels = TupleTools.setindex(levels, levels[s+1], s)
            levels = TupleTools.setindex(levels, l, s+1)
            trees, newtrees = newtrees, trees
            empty!(newtrees)
        end
        return trees
    end
end

# permute fusion tree
"""
    permute(f::FusionTree, p::NTuple{N,Int}) -> <:AbstractDict{typeof(t),<:Number}

Perform a permutation of the uncoupled indices of the fusion tree `f` and returns the result
as a `<:AbstractDict` of output trees and corresponding coefficients; this requires that
`BraidingStyle(sectortype(f)) isa SymmetricBraiding`.
"""
function permute(f::FusionTree{G,N}, p::NTuple{N,Int}) where {G<:Sector, N}
    @assert BraidingStyle(G) isa SymmetricBraiding
    return braid(f, ntuple(identity, Val(N)), p)
end

"""
    split(f::FusionTree{G,N}, ::StaticLength(M))
    -> (::FusionTree{G,M}, ::FusionTree{G,N-M+1})

Split a fusion tree with the first M outgoing indices, and an incoming index corresponding
to the internal fusion tree index between outgoing indices N and N+1 of the original tree
`f`; and a second fusion tree whose first outgoing index is that same internal index. Its
remaining outgoing indices are the N-M outgoing indices of the original tree `f`, and also
the incoming index is the same. This is in the inverse of `insertat` in the sense that if
`f1, f2 = split(t, StaticLength(M)) ⇒ f == insertat(f2, 1, f1)`.
"""
split(f::FusionTree{G,N}, ::StaticLength{N}) where {G,N} =
    (f, FusionTree{G}((f.coupled,), f.coupled, (false,), (), ()))
split(f::FusionTree{G,N}, ::StaticLength{1}) where {G,N} =
    (FusionTree{G}((f.uncoupled[1],), f.uncoupled[1], (false,), (), ()), f)
function split(f::FusionTree{G,N}, ::StaticLength{0}) where {G,N}
    f1 = FusionTree{G}((), one(G), (), ())
    uncoupled2 = (one(G), f.uncoupled...)
    coupled2 = f.coupled
    isdual2 = (false, f.isdual...)
    innerlines2 = N >= 2 ? (f.uncoupled[1], f.innerlines...) : ()
    if FusionStyle(G) isa DegenerateNonAbelian
        vertices2 = (1, f.vertices...)
        return f1, FusionTree(uncoupled2, coupled2, isdual2, innerlines2, vertices2)
    else
        return f1, FusionTree(uncoupled2, coupled2, isdual2, innerlines2)
    end
end
function split(f::FusionTree{G,N}, ::StaticLength{M}) where {G,N,M}
    @assert 1 < M < N
    uncoupled1 = ntuple(n->f.uncoupled[n], Val(M))
    isdual1 = ntuple(n->f.isdual[n], Val(M))
    innerlines1 = M>2 ? ntuple(n->f.innerlines[n], Val(M-2)) : ()
    coupled1 = f.innerlines[M-1]
    vertices1 = ntuple(n->f.vertices[n], Val(M-1))
    f1 = FusionTree(uncoupled1, coupled1, isdual1, innerlines1, vertices1)

    uncoupled2 = (coupled1, ntuple(n->f.uncoupled[M+n], Val(N-M))...)
    isdual2 = (false, ntuple(n->f.isdual[M+n], Val(N-M))...)
    innerlines2 = ntuple(n->f.innerlines[M-1+n], Val(N-M-1))
    coupled2 = f.coupled
    vertices2 = ntuple(n->f.vertices[M-1+n], Val(N-M))
    f2 = FusionTree(uncoupled2, coupled2, isdual2, innerlines2, vertices2)
    return f1, f2
end

"""
    merge(f1::FusionTree{G,N₁}, f2::FusionTree{G,N₂}, c::G, μ = nothing)
    -> <:AbstractDict{<:FusionTree{G,N₁+N₂},<:Number}

Merge two fusion trees together to a linear combination of fusion trees whose uncoupled
sectors are those of `f1` followed by those of `f2`, and where the two coupled sectors of
`f1` and `f2` are further fused to `c`. In case of
`FusionStyle(G) == DegenerateNonAbelian()`, also a degeneracy label `μ` for the fusion of
the coupled sectors of `f1` and `f2` to `c` needs to be specified.
"""
function merge(f1::FusionTree{G,N₁}, f2::FusionTree{G,N₂},
                    c::G, μ = nothing) where {G,N₁,N₂}
    if !(c in f1.coupled ⊗ f2.coupled)
        throw(SectorMismatch("cannot fuse sectors $(f1.coupled) and $(f2.coupled) to $c"))
    end
    f0 = FusionTree((f1.coupled, f2.coupled), c, (false, false), (), (μ,))
    f, coeff = first(insertat(f0, 1, f1)) # takes fast path, single output
    @assert coeff == one(coeff)
    return insertat(f, N₁+1, f2)
end
function merge(f1::FusionTree{G,0}, f2::FusionTree{G,0}, c::G, μ =nothing) where {G}
    c == one(G) ||
        throw(SectorMismatch("cannot fuse sectors $(f1.coupled) and $(f2.coupled) to $c"))
    return SingletonDict(f1=>Fsymbol(one(G),one(G),one(G),one(G),one(G),one(G)))
end

"""
    insertat(f::FusionTree{G,N₁}, i, f2::FusionTree{G,N₂})
    -> <:AbstractDict{<:FusionTree{G,N₁+N₂-1},<:Number}

Attach a fusion tree `f2` to the uncoupled leg `i` of the fusion tree `f1` and bring it
into a linear combination of fusion trees in standard form. This requires that
`f2.coupled == f1.uncoupled[i]` and `f1.isdual[i] == false`.
"""
function insertat(f1::FusionTree{G}, i, f2::FusionTree{G,0}) where {G}
    # this actually removes uncoupled line i, which should be trivial
    (f1.uncoupled[i] == f2.coupled && !f1.isdual[i]) ||
        throw(SectorMismatch("cannot connect $(f2.uncoupled) to $(f1.uncoupled[i])"))
    coeff = Fsymbol(one(G), one(G), one(G), one(G), one(G), one(G))

    uncoupled = TupleTools.deleteat(f1.uncoupled, i)
    coupled = f1.coupled
    isdual = TupleTools.deleteat(f1.isdual, i)
    inner = TupleTools.deleteat(f1.innerlines, max(1,i-2))
    vertices = TupleTools.deleteat(t1.vertices, max(1, i-1))
    f = FusionTree(uncoupled, coupled, isdual, inner, vertices)
    if FusionStyle(G) isa Abelian
        return SingletonDict(f => coeff)
    elseif FusionStyle(G) isa SimpleNonAbelian
        return FusionTreeDict(f => coeff)
    end
end
function insertat(f1::FusionTree{G}, i, f2::FusionTree{G,1}) where {G}
    # identity operation
    (f1.uncoupled[i] == f2.coupled && !f1.isdual[i]) ||
        throw(SectorMismatch("cannot connect $(f2.uncoupled) to $(f1.uncoupled[i])"))
    coeff = Fsymbol(one(G), one(G), one(G), one(G), one(G), one(G))
    isdual′ = TupleTools.setindex(f1.isdual, f2.isdual[1], i)
    f = FusionTree{G}(f1.uncoupled, f1.coupled, isdual′, f1.innerlines, f1.vertices)
    if FusionStyle(G) isa Abelian
        return SingletonDict(f => coeff)
    elseif FusionStyle(G) isa SimpleNonAbelian
        return FusionTreeDict(f => coeff)
    end
end
function insertat(f1::FusionTree{G}, i, f2::FusionTree{G,2}) where {G}
    # elementary building block,
    (f1.uncoupled[i] == f2.coupled && !f1.isdual[i]) ||
        throw(SectorMismatch("cannot connect $(f2.uncoupled) to $(f1.uncoupled[i])"))
    uncoupled = f1.uncoupled
    coupled = f1.coupled
    inner = f1.innerlines
    b, c = f2.uncoupled
    isdual = f1.isdual
    isdualb, isdualc = f2.isdual
    if i == 1
        uncoupled′ = (b, c, tail(uncoupled)...)
        isdual′ = (isdualb, isdualc, tail(isdual)...)
        inner′ = (uncoupled[1], inner...)
        vertices′ = (f2.vertices..., f1.vertices...)
        coeff = Fsymbol(one(G), one(G), one(G), one(G), one(G), one(G))
        f′ = FusionTree(uncoupled′, coupled, isdual′, inner′, vertices′)
        if FusionStyle(G) isa Abelian
            return SingletonDict(f′ => coeff)
        elseif FusionStyle(G) isa SimpleNonAbelian
            return FusionTreeDict(f′ => coeff)
        end
    end
    uncoupled′ = TupleTools.insertafter(TupleTools.setindex(uncoupled, b, i), i, (c,))
    isdual′ = TupleTools.insertafter(TupleTools.setindex(isdual, isdualb, i), i, (isdualc,))
    a = i == 2 ? uncoupled[1] : inner[i-2]
    d = i == length(f1) ? coupled : inner[i-1]
    e′ = uncoupled[i]
    if FusionStyle(G) isa Abelian
        e = first(a ⊗ b)
        inner′ = TupleTools.insertafter(inner, i-2, (e,))
        f′ = FusionTree(uncoupled′, coupled, isdual′, inner′)
        coeff = conj(Fsymbol(a,b,c,d,e,e′))
        return SingletonDict(f′ => coeff)
    elseif FusionStyle(G) isa SimpleNonAbelian
        local newtrees
        for e in a ⊗ b
            inner′ = TupleTools.insertafter(inner, i-2, (e,))
            f′ = FusionTree(uncoupled′, coupled, isdual′, inner′)
            coeff = conj(Fsymbol(a,b,c,d,e,e′))
            if coeff != zero(coeff)
                if @isdefined newtrees
                    newtrees[f′] = coeff
                else
                    newtrees = FusionTreeDict(f′ => coeff)
                end
            end
        end
        return newtrees
    else
        # TODO: implement DegenerateNonAbelian case
        throw(MethodError(insertat, (f1, i, f2)))
    end
end
function insertat(f1::FusionTree{G}, i, f2::FusionTree{G}) where {G}
    (f1.uncoupled[i] == f2.coupled && !f1.isdual[i]) ||
        throw(SectorMismatch("cannot connect $(f2.uncoupled) to $(f1.uncoupled[i])"))
    if length(f1) == 1
        coeff = Fsymbol(one(G), one(G), one(G), one(G), one(G), one(G))
        if FusionStyle(G) isa Abelian
            return SingletonDict(f2 => coeff)
        elseif FusionStyle(G) isa SimpleNonAbelian
            return FusionTreeDict(f2 => coeff)
        end
    end
    if i == 1
        uncoupled = (f2.uncoupled..., tail(f1.uncoupled)...)
        isdual = (f2.isdual..., tail(f1.isdual)...)
        inner = (f2.innerlines..., f2.coupled, f1.innerlines...)
        vertices = (f2.vertices..., f1.vertices...)
        coupled = f1.coupled
        f′ = FusionTree(uncoupled, coupled, isdual, inner, vertices)
        coeff = Fsymbol(one(G), one(G), one(G), one(G), one(G), one(G))
        if FusionStyle(G) isa Abelian
            return SingletonDict(f′ => coeff)
        elseif FusionStyle(G) isa SimpleNonAbelian
            return FusionTreeDict(f′ => coeff)
        end
    else # recursive definition
        N2 = length(f2)
        f2′, f2′′ = split(f2, StaticLength(N2) - StaticLength(1))
        if FusionStyle(G) isa Abelian
            f, coeff = first(insertat(f1, i, f2′′))
            f′, coeff′ = first(insertat(f, i, f2′))
            return SingletonDict(f′=>coeff*coeff′)
        else
            local newtrees
            for (f, coeff) in insertat(f1, i, f2′′)
                if @isdefined newtrees
                    for (f′, coeff′) in insertat(f, i, f2′)
                        newtrees[f′] = get(newtrees, f′, zero(coeff′)) + coeff*coeff′
                    end
                else
                    newtrees = insertat(f, i, f2′)
                    for (f′, coeff′) in newtrees
                        newtrees[f′] = coeff*coeff′
                    end
                end
            end
            return newtrees
        end
    end
end

# repartition double fusion tree
"""
    repartition(f1::FusionTree{G,N₁}, f2::FusionTree{G,N₂},
                ::StaticLength{N}) where {G,N₁,N₂,N}
    -> <:AbstractDict{Tuple{FusionTree{G,N}, FusionTree{G,N₁+N₂-N}},<:Number}

Input is a double fusion tree that describes the fusion of a set of incoming uncoupled
sectors to a set of outgoing uncoupled sectors, represented using the individual trees of
outgoing (`f1`) and incoming sectors (`f2`) respectively (with identical coupled sector
`f1.coupled == f2.coupled`). Computes new trees and corresponding coefficients obtained from
repartitioning the tree by bending incoming to outgoing sectors (or vice versa) in order to
have `N` outgoing sectors.
"""
function repartition(f1::FusionTree{G,N₁},
                        f2::FusionTree{G,N₂},
                        V::StaticLength{N}) where {G<:Sector, N₁, N₂, N}
    f1.coupled == f2.coupled || throw(SectorMismatch())
    @assert 0 <= N <= N₁+N₂
    V1 = V
    V2 = StaticLength(N₁)+StaticLength(N₂)-V

    if FusionStyle(f1) isa Abelian || FusionStyle(f1) isa SimpleNonAbelian
        coeff = sqrt(dim(one(G)))*Bsymbol(one(G), one(G), one(G))
        uncoupled = (f1.uncoupled..., map(dual, reverse(f2.uncoupled))...)
        isdual = (f1.isdual..., map(!, reverse(f2.isdual))...)
        inner1ext = isa(StaticLength(N₁), StaticLength{0}) ? () :
                        (isa(StaticLength(N₁), StaticLength{1}) ? (one(G),) :
                            (one(G), first(uncoupled), f1.innerlines...))
        inner2ext = isa(StaticLength(N₂), StaticLength{0}) ? () :
                        (isa(StaticLength(N₂), StaticLength{1}) ? (one(G),) :
                            (one(G), dual(last(uncoupled)), f2.innerlines...))
        innerext = (inner1ext..., f1.coupled, reverse(inner2ext)...) # length N₁+N₂+1
        for n = N₁+1:N
             # map fusion vertex c<-(a,b) to splitting vertex (c,dual(b))<-a
            b = dual(uncoupled[n])
            a = innerext[n+1]
            c = innerext[n]
            coeff *= sqrt(dim(c)/dim(a))*conj(Bsymbol(a,b,c))
            if !isdual[n]
                coeff *= frobeniusschur(dual(b))
            end
        end
        for n = N₁:-1:N+1
            # map splitting vertex (a,b)<-c to fusion vertex a<-(c,dual(b))
            b = uncoupled[n]
            a = innerext[n]
            c = innerext[n+1]
            coeff *= sqrt(dim(c)/dim(a))*Bsymbol(a,b,c)
            if isdual[n]
                coeff *= conj(frobeniusschur(dual(b)))
            end
        end
        uncoupled1 = TupleTools.getindices(uncoupled, ntuple(n->n, V1))
        uncoupled2 = TupleTools.getindices(map(dual,uncoupled), ntuple(n->N₁+N₂+1-n, V2))
        isdual1 = TupleTools.getindices(isdual, ntuple(n->n, V1))
        isdual2 = TupleTools.getindices(map(!,isdual), ntuple(n->N₁+N₂+1-n, V2))
        innerlines1 = TupleTools.getindices(innerext, ntuple(n->n+2, V1 - StaticLength(2)))
        innerlines2 = TupleTools.getindices(innerext, ntuple(n->N₁+N₂-n, V2 - StaticLength(2)))
        c = innerext[N+1]
        f1′ = FusionTree{G}(uncoupled1, c, isdual1, innerlines1)
        f2′ = FusionTree{G}(uncoupled2, c, isdual2, innerlines2)
        return SingletonDict((f1′, f2′)=>coeff)
    else
        # TODO: implement DegenerateNonAbelian case
        throw(MethodError(repartition, (f1, f2, V)))
    end
end

# braid double fusion tree
const braidcache = LRU{Any,Any}(; maxsize = 10^5)
const usebraidcache_abelian = Ref{Bool}(false)
const usebraidcache_nonabelian = Ref{Bool}(true)

"""
    braid(f1::FusionTree{G}, f2::FusionTree{G},
            levels1::IndexTuple, levels2::IndexTuple,
            p1::IndexTuple{N₁}, p2::IndexTuple{N₂}) where {G<:Sector,N₁,N₂}
    -> <:AbstractDict{Tuple{FusionTree{G,N₁}, FusionTree{G,N₂}},<:Number}

Input is a fusion-splitting tree pair that describes the fusion of a set of incoming
uncoupled sectors to a set of outgoing uncoupled sectors, represented using the splitting
tree `f1` and fusion tree `f2`, such that the incoming sectors `f2.uncoupled` are fused to
`f1.coupled == f2.coupled` and then to the outgoing sectors `f1.uncoupled`. Compute new
trees and corresponding coefficients obtained from repartitioning and braiding the tree such
that sectors `p1` become outgoing and sectors `p2` become incoming. The uncoupled indices in
splitting tree `f1` and fusion tree `f2` have levels (or depths) `levels1` and `levels2`
respectively, which determines how indices braid. In particular, if `i` and `j` cross,
``τ_{i,j}`` is applied if `levels[i] < levels[j]` and ``τ_{j,i}^{-1}`` if `levels[i] >
levels[j]`. This does not allow to encode the most general braid, but a general braid can
be obtained by combining such operations.
"""
function braid(f1::FusionTree{G}, f2::FusionTree{G},
                levels1::IndexTuple, levels2::IndexTuple,
                p1::IndexTuple{N₁}, p2::IndexTuple{N₂}) where {G<:Sector,N₁,N₂}
    @assert length(f1) + length(f2) == N₁ + N₂
    @assert length(f1) == length(levels1) && length(f2) == length(levels2)
    @assert TupleTools.isperm((p1..., p2...))
    if FusionStyle(f1) isa Abelian &&
        BraidingStyle(f1) isa SymmetricBraiding
        if usebraidcache_abelian[]
            u = one(G)
            T = typeof(sqrt(dim(u))*Fsymbol(u,u,u,u,u,u)*Rsymbol(u,u,u))
            F₁ = fusiontreetype(G, StaticLength(N₁))
            F₂ = fusiontreetype(G, StaticLength(N₂))
            D = SingletonDict{Tuple{F₁,F₂}, T}
            return _get_braid(D, (f1, f2, levels1, levels2, p1, p2))
        else
            return _braid((f1, f2, levels1, levels2, p1, p2))
        end
    else
        if usebraidcache_nonabelian[]
            u = one(G)
            T = typeof(sqrt(dim(u))*Fsymbol(u,u,u,u,u,u)*Rsymbol(u,u,u))
            F₁ = fusiontreetype(G, StaticLength(N₁))
            F₂ = fusiontreetype(G, StaticLength(N₂))
            D = FusionTreeDict{Tuple{F₁,F₂}, T}
            return _get_braid(D, (f1, f2, levels1, levels2, p1, p2))
        else
            return _braid((f1, f2, levels1, levels2, p1, p2))
        end
    end
end

@noinline function _get_braid(::Type{D}, @nospecialize(key)) where D
    d::D = get!(braidcache, key) do
        _braid(key)
    end
    return d
end

const BraidKey{G<:Sector,N₁,N₂} = Tuple{<:FusionTree{G}, <:FusionTree{G},
                                        IndexTuple, IndexTuple,
                                        IndexTuple{N₁}, IndexTuple{N₂}}

function _braid((f1, f2, l1, l2, p1, p2)::BraidKey{G,N₁,N₂}) where {G<:Sector,N₁,N₂}
    p = linearizepermutation(p1, p2, length(f1), length(f2))
    levels = (l1..., reverse(l2)...)
    if FusionStyle(f1) isa Abelian
        (f,f0), coeff1 = first(repartition(f1, f2, StaticLength(N₁) + StaticLength(N₂)))
        f, coeff2 = first(braid(f, levels, p))
        (f1′,f2′), coeff3 = first(repartition(f, f0, StaticLength(N₁)))
        return SingletonDict((f1′,f2′)=>coeff1*coeff2*coeff3)
    elseif FusionStyle(f1) isa SimpleNonAbelian
        (f,f0), coeff1 = first(repartition(f1, f2, StaticLength(N₁) + StaticLength(N₂)))
        local newtrees
        for (f, coeff2) in braid(f, levels, p)
            (f1′, f2′), coeff3 = first(repartition(f, f0, StaticLength(N₁)))
            if @isdefined newtrees
                newtrees[(f1′,f2′)] = coeff1*coeff2*coeff3
            else
                newtrees = FusionTreeDict((f1′,f2′)=>coeff1*coeff2*coeff3)
            end
        end
        return newtrees
    else
        # TODO: implement DegenerateNonAbelian case
        throw(MethodError(braid, (f1, f2, l1, l2, p1, p2)))
    end
end

"""
    permute(f1::FusionTree{G}, f2::FusionTree{G},
            p1::NTuple{N₁,Int}, p2::NTuple{N₂,Int}) where {G,N₁,N₂}
    -> <:AbstractDict{Tuple{FusionTree{G,N₁}, FusionTree{G,N₂}},<:Number}

Input is a double fusion tree that describes the fusion of a set of incoming uncoupled
sectors to a set of outgoing uncoupled sectors, represented using the individual trees of
outgoing (`t1`) and incoming sectors (`t2`) respectively (with identical coupled sector
`t1.coupled == t2.coupled`). Computes new trees and corresponding coefficients obtained from
repartitioning and permuting the tree such that sectors `p1` become outgoing and sectors
`p2` become incoming.
"""
function permute(f1::FusionTree{G}, f2::FusionTree{G},
                    p1::IndexTuple{N₁}, p2::IndexTuple{N₂}) where {G<:Sector,N₁,N₂}
    @assert BraidingStyle(G) isa SymmetricBraiding
    levels1 = ntuple(identity, length(f1))
    levels2 = length(f1) .+ ntuple(identity, length(f2))
    return braid(f1, f2, levels1, levels2, p1, p2)
end


text:
include(joinpath("linReg", "research.jl"))
include(joinpath("linReg", "chosenVariables.jl"))

max_stars_repo_name: AlexLsn/CSO

max_issues_repo_name: AlexLsn/CSO

max_forks_repo_name: AlexLsn/CSO


text:
using Test
using LightGraphs
using SimpleWeightedGraphs
using Serialization
using ChemistryFeaturization
using Xtals

# NB: featurizing graphs is tested in ElementFeatureDescriptor and GraphNodeFeaturization tests

@testset "AtomGraph" begin

    wm_mp195 = [0.0 1.0 1.0 1.0; 1.0 0.0 1.0 1.0; 1.0 1.0 0.0 1.0; 1.0 1.0 1.0 0.0]
    lapl_mp195 = [
        1.0 -0.3333333 -0.3333333 -0.3333333
        -0.3333333 1.0 -0.3333333 -0.3333333
        -0.3333333 -0.3333333 1.0 -0.3333333
        -0.3333333 -0.3333333 -0.3333333 1.0
    ]
    els_mp195 = ["Ho", "Pt", "Pt", "Pt"]
    path_mp195 = abspath(@__DIR__, "..", "test_data", "strucs", "mp-195.cif")

    @testset "construct object" begin
        # build a cute little triangle graph
        adj = Float32.([0 1 1; 1 0 1; 1 1 0])
        g = SimpleWeightedGraph(adj)

        # add an element list that doesn't make sense
        @test_throws AssertionError AtomGraph(g, ["C"])

        # okay, now do it right, start with no features, checking both constructors
        ag = AtomGraph(g, ["C", "C", "C"])
        ag2 = AtomGraph(adj, ["C", "C", "C"])
        @test adjacency_matrix(ag.graph) == adjacency_matrix(ag2.graph)

        # and now from a file...
        ag = AtomGraph(path_mp195)
        # and from a Crystal object
        c = Crystal(path_mp195)
        ag2 = AtomGraph(c)

        @test weights(ag.graph) == weights(ag2.graph) == wm_mp195
        @test all(isapprox.(ag.laplacian, lapl_mp195, atol = 1e-7))
        @test all(isapprox.(ag2.laplacian, lapl_mp195, atol = 1e-7))
        @test ag.elements == ag2.elements == els_mp195

        # test that warning is thrown for NaNs in laplacian
        @test_throws ArgumentError AtomGraph(
            build_graph(
                abspath(@__DIR__, "..", "test_data", "strucs", "nanlaplstruc.cif"),
            )...,
        )
    end

    @testset "save/load" begin
        g = SimpleWeightedGraph(Float32.([0 1 1; 1 0 1; 1 1 0]))
        ag = AtomGraph(g, ["C", "C", "C"])
        serialize(abspath(@__DIR__, "..", "test_data", "strucs", "testgraph.jls"), ag)
        ag2 = deserialize(abspath(@__DIR__, "..", "test_data", "strucs", "testgraph.jls"))
        @test adjacency_matrix(ag.graph) == adjacency_matrix(ag2.graph)
        @test elements(ag) == ag2.elements
        @test ChemistryFeaturization.Atoms.normalized_laplacian(ag) == ag2.laplacian
    end

    @testset "batch processing" begin
        file_list = readdir(abspath(@__DIR__, "..", "test_data", "strucs"), join = true)
        ags = AtomGraph.(file_list)
        @test length(collect(skipmissing(ags))) == length(ags) - 1

        ids = map(f -> splitpath(f)[end], file_list)
        ags2 = collect(skipmissing(AtomGraph.(file_list, ids)))
        mp195s = [ag for ag in skipmissing(ags2) if contains(ag.id, "mp-195")]

        for ag in mp195s
            @test elements(ag) == els_mp195
            @test weights(ag.graph) == wm_mp195
            @test all(isapprox.(ag.laplacian, lapl_mp195, atol = 1e-7))
        end

    end

    @testset "visualize" begin
        # test lt_edge
        adj = Float32.([0 1 2; 1 0 1; 2 1 0])
        g = SimpleWeightedGraph{Int32}(adj)
        e_adj = [e for e in edges(g)]

        @test ChemistryFeaturization.Atoms.lt_edge(e_adj[2], e_adj[3]) == true # e1.src < e2.src
        @test ChemistryFeaturization.Atoms.lt_edge(e_adj[1], e_adj[2]) == true # e1.dest < e2.dest
        @test ChemistryFeaturization.Atoms.lt_edge(e_adj[3], e_adj[2]) == false
        adj = Float32.([0 1 1; 1 0 1; 1 1 0])
        g = SimpleWeightedGraph{Int32}(adj)
        ag = AtomGraph(g, ["C", "C", "C"])
        @test ChemistryFeaturization.Atoms.graph_edgewidths(ag) == [1.0, 1.0, 1.0]

        # test graph_colors
        ag = AtomGraph(g, ["O", "C", "O"])
        c1, c2, c3 = ChemistryFeaturization.Atoms.graph_colors(elements(ag))
        @test c1 == c3 != c2

        # TODO: figure out a way to test visualize_graph itself?
    end

end


text:
module KernelLibrary

export linear_kernel, polynomial_kernel, laplacian_kernel, gaussian_kernel, anova_kernel, rbf_kernel, fractrbf_kernel, rationalquad_kernel, multiquad_kernel, invmultiquad_kernel, circular_kernel, spherical_kernel, wave_kernel
export power_kernel, log_kernel, spline_kernel, cauchy_kernel, chisquare_kernel, histintersect_kernel, genhistintersect_kernel, gentstudent_kernel, fourier_kernel, sigmoid_kernel, wavelet_kernel, transinvwavelet_kernel
# Julia implementation of most of the kernel functions from
# crsouza.blogspot.com/2010/03/kernel-functions-for-machine-learning.html

function linear_kernel{T}(x::Array{T}, y::Vector{T}; c::Real=0.)
	c=convert(T,c)
	(y'*x .+ c)
end

function polynomial_kernel{T}(x::Array{T}, y::Vector{T}, d::Real; α::Real=1., c::Real=1.)
	d=convert(T,d); α=convert(T,α); c=convert(T,c);
	(α.*(y'*x) .+ c).^ d
end

function laplacian_kernel{T}(x::Array{T}, y::Vector{T}, σ::Real)
	σ=convert(T, σ)
	exp(-sum((x.-y).^2.,1).^.5 ./ σ)
end

function gaussian_kernel{T}(x::Array{T}, y::Vector{T}, σ::Real)
	σ=convert(T, σ)
	exp(-sum((x.-y).^2,1) ./ 2*σ.^ 2)
end

function anova_kernel{T}(x::Array{T}, y::Vector{T}, σ::Real, d::Real)
	d=convert(T,d); σ=convert(T, σ)
	sum(exp(-σ*(x.-y).^2).^d,1)
end

function rbf_kernel{T}(x::Array{T}, y::Vector{T}, σ::Real)
	σ=convert(T, σ)
	exp(-sum(abs(x.-y),1) ./ 2*σ.^ 2)
end

function fractrbf_kernel{T}(x::Array{T}, y::Vector{T}, σ::Real, fract::Real)
	σ=convert(T, σ); fract=convert(T, fract);
	exp(-sum(abs(x .- y).^fract,1).^(1/fract) ./ 2*σ.^ 2)
end

function rationalquad_kernel{T}(x::Array{T}, y::Vector{T}, c::Real)
	c = convert(T, c)
	n = sum((x.-y).^2,1)
	1 .- n./(n .+ c)
end

function multiquad_kernel{T}(x::Array{T}, y::Vector{T}, c::Real)
	c = convert(T, c)
	sqrt(sum((x.-y).^2,1) .+ c^2)
end

function invmultiquad_kernel{T}(x::Array{T}, y::Vector{T}, c::Real)
	c = convert(T, c)
	1 ./ sqrt(sum((x.-y).^2,1) .+ c^2)
end

function circular_kernel{T}(x::Array{T}, y::Vector{T}, σ::Real)
	σ=convert(T, σ)
	n = sqrt(sum((x.-y).^2,1))./σ
	out = (2/π).*(acos(-n) - n.*sqrt(1 .- n.^2))
	out[n.>=σ] = 0
	out
end

function spherical_kernel{T}(x::Array{T}, y::Vector{T}, σ::Real)
	σ=convert(T, σ)
	n = sqrt(sum((x.-y).^2,1))./σ
	out = 1 - 1.5.*n + .5.*n.^3
	out[n.>=σ] = 0
	out
end

function wave_kernel{T}(x::Array{T}, y::Vector{T}, θ::Real)
	θ=convert(T, θ)
	n = sqrt(sum((x.-y).^2,1))
	(θ./n).*sin(n./θ)
end

function power_kernel{T}(x::Array{T}, y::Vector{T}, d::Real)
	d=convert(T, d)
	-sqrt(sum((x.-y).^2,1)).^d
end

function log_kernel{T}(x::Array{T}, y::Vector{T}, d::Real)
	d=convert(T, d)
	-log(sqrt(sum((x.-y).^2.,1)).^d + 1)
end

function spline_kernel{T}(x::Array{T}, y::Vector{T})
	xy = x.*y
	mxy = broadcast(min, x,y)
    z = 1 .+ xy .+ xy.*mxy - ((x.+y)/2).*mxy.^2 + (mxy.^3)/3
    prod(z,1)
end

function cauchy_kernel{T}(x::Array{T}, y::Vector{T}, σ::Real)
	σ=convert(T, σ)
	1 ./ (1 .+ (sum((x.-y).^2,1) ./ σ.^2))
end

function chisquare_kernel{T}(x::Array{T}, y::Vector{T})
	1 .- sum((x.-y).^2 ./ ((x.+y)/2),1)
end

function histintersect_kernel{T}(x::Array{T}, y::Vector{T})
	sum(broadcast(min, x,y),1)
end

function genhistintersect_kernel{T}(x::Array{T}, y::Vector{T}, α::Real, β::Real)
	α=convert(T, α); β=convert(T, β);
	sum(broadcast(min, abs(x).^α,abs(y).^β),1)
end

function gentstudent_kernel{T}(x::Array{T}, y::Vector{T}, d::Real)
	d=convert(T, d)
	1 ./ (1.0 .+ sqrt(sum((x.-y).^2.,1)).^d)
end

function fourier_kernel{T}(x::Array{T}, y::Vector{T}, a::Real)
	a=convert(T, a)
    dist = x.-y;
    dist[dist.==0] = eps(eltype(x))
    z = sin((a + 1/2)*dist)./sin(dist/2);
    prod(z,1)
end

function sigmoid_kernel{T}(x::Array{T}, y::Vector{T}; α::Real=1, c::Real=0)
	α=convert(T, α); c=convert(T, c);
	tanh(α*y'*x + c)
end

motherwavelet(x::Array) = cos(1.75.*x).*exp(-(x.^2)./2)

function wavelet_kernel{T}(x::Array{T}, y::Vector{T}, a::Real, c::Real)
	a=convert(T, a); c=convert(T, c);
	prod(motherwavelet((x.-c)./a).*motherwavelet((y.-c)./a),1)
end

function transinvwavelet_kernel{T}(x::Array{T}, y::Vector{T}, a::Real)
	a=convert(T, a);
	prod(motherwavelet((x.-y)./a),1)
end

end # module


text:
#=
Conservative 7-color palette from Points of view: Color blindness, Bang Wong - Nature Methods
https://www.nature.com/articles/nmeth.1618?WT.ec_id=NMETH-201106
=#

const wong_colors = [
    RGB(230/255, 159/255, 0/255),
    RGB(86/255, 180/255, 233/255),
    RGB(0/255, 158/255, 115/255),
    RGB(240/255, 228/255, 66/255),
    RGB(0/255, 114/255, 178/255),
    RGB(213/255, 94/255, 0/255),
    RGB(204/255, 121/255, 167/255),
]

const default_palettes = Attributes(
    color = wong_colors,
    marker = [:circle, :xcross, :utriangle, :diamond, :dtriangle, :star8, :pentagon, :rect],
    linestyle = [nothing, :dash, :dot, :dashdot, :dashdotdot],
    side = [:left, :right]
)

const minimal_default = Attributes(
    palette = default_palettes,
    font = "Dejavu Sans",
    backgroundcolor = :white,
    color = :black,
    colormap = :viridis,
    marker = Circle,
    markersize = 0.1,
    linestyle = nothing,
    resolution = reasonable_resolution(),
    visible = true,
    clear = true,
    show_axis = true,
    show_legend = false,
    scale_plot = true,
    center = true,
    update_limits = true,
    axis = Attributes(),
    axis2d = Attributes(),
    axis3d = Attributes(),
    legend = Attributes(),
    axis_type = automatic,
    camera = automatic,
    limits = automatic,
    padding = Vec3f0(0.1),
    raw = false
)

const _current_default_theme = copy(minimal_default)

function current_default_theme(; kw_args...)
    return merge!(Attributes(kw_args), _current_default_theme)
end

function set_theme!(new_theme::Attributes)
    empty!(_current_default_theme)
    new_theme = merge!(new_theme, minimal_default)
    merge!(_current_default_theme, new_theme)
    return
end
function set_theme!(;kw_args...)
    set_theme!(Attributes(; kw_args...))
end


text:
using Test

using Pepper
import Pepper.CLI

let

    n = rand(3:8)
    cra = ["criterion-$i" for i=1:n]
    lvs = [["level-$i-$j" for j=1:rand(2:5)] for i=1:n]

    s1 = Pepper.System(cra, lvs)
    s2 = Pepper.System(cra, lvs)

    pts = [sort([rand(Int) for i in lv]) for lv in lvs]

    for i = 1:n, j = 1:length(pts[i])
        Pepper.setpts(s1, i, j, pts[i][j])
    end

    k = 2

    function oracle(_, base_idx, lhs, rhs)
        score1 = sum(Pepper.getpts(s1, i, j) for (i, j) in zip(base_idx, lhs))
        score2 = sum(Pepper.getpts(s1, i, j) for (i, j) in zip(base_idx, rhs))
        if score1 == score2
            return 1
        elseif score1 > score2
            return 2
        else
            return 3
        end
    end

    Pepper.paprika(oracle, s2, maxdegree=k)

end

let
    lines = split("""
    Cyanaurate:
        - mataeologue
        - neoclassicism
        - preponder
        - refractility
    Microwave:
        - polygrammatic
        - sugary
        - unpictorial
    Nonvariant:
        - telharmony
        - luteway
    Planeticose:
        - telharmony
        - luteway
    Archil:
        - facilitation
        - hairlock
        - anvilsmith
        - soursop
    """, "\n")
    s = Pepper.System()
    CLI.parse(lines, s)
    cases = []
    Pepper.paprika(s, maxdegree=3) do s, base_idx, lhs, rhs
        choice = rand(1:3)
        push!(cases, [base_idx, lhs, rhs, choice])
        choice
    end
    for case in cases
        base_idx, lhs, rhs, choice = case
        score1 = sum(Pepper.getpts(s, i, j) for (i, j) in zip(base_idx, lhs))
        score2 = sum(Pepper.getpts(s, i, j) for (i, j) in zip(base_idx, rhs))
        if choice == 1
            @test score1 == score2
        else
            if choice == 2
                score1, score2 = score2, score1
            end
            @test score1 < score2
        end
    end
end

let
    lines = split("""
    Cyanaurate:
        - mataeologue
        - neoclassicism
        - refractility
    Microwave:
        - polygrammatic
        - sugary
        - unpictorial: 5
    """, "\n")
    s = Pepper.System()
    CLI.parse(lines, s)

    # Pepper.setpts(s, "Microwave" => "unpictorial", 5)

    cases = []
    Pepper.paprika(s, ptsum=100) do s, base_idx, lhs, rhs
        choice = rand(1:2)
        push!(cases, [base_idx, lhs, rhs, choice])
        choice
    end

    Y = Pepper.shape(s)
    n = length(Y)

    @test sum(Pepper.getpts(s, i, Y[i]) for i = 1:n) == 100
    @test Pepper.getpts(s, "Microwave" => "unpictorial") == 5
end


text:
"""
    render(x::MyType)

Generic function that defines how a Julia object is rendered.
Should return a `Node` object.

# Examples
```julia
struct MyPlot
    s::Scope
end

WebIO.render(p::MyPlot) = WebIO.render(p.s)
```
"""
function render end
render(x::Union{Node, String}) = x
render(x::Text) = dom"pre"(x.content)
render(::Nothing) = ""
render(x::Any) =
    dom"div"(; setInnerHtml=richest_html(x))

const renderable_types = Type[]
"""
    `WebIO.register_renderable(MyType::Type)`

Registers that a WebIO.render method is available for instances of `MyType`.
Allows WebIO to hook into the display machinery of backends such as Atom and
IJulia to display the WebIO rendered version of the type as appropriate.

Also defines a `Base.show(io::IO, m::MIME"text/html", x::MyType)` as
`Base.show(io, m, WebIO.render(x))`
"""
function register_renderable(::Type{T}) where T
    @eval Base.show(io::IO, m::MIME"text/html", x::$T) = Base.show(io, m, WebIO.render(x))
    push!(renderable_types, T)
    return true
end

"""
Called after a provider is setup
"""
function re_register_renderables()
    for provider in providers_initialised
        for T in renderable_types
            register_renderable(T, Val(provider))
        end
    end
end

"""
Generic Conversion to Nodes
"""
const mime_order = map(MIME, [ "text/html", "text/latex", "image/svg+xml", "image/png", "image/jpeg", "text/markdown", "application/javascript", "text/plain" ])

function richest_mime(val)
    for mimetype in mime_order
        showable(mimetype, val) && return mimetype
    end
    error("value not writable for any mimetypes")
end

function richest_html(val)
    topmime = richest_mime(val)
    str_repr = stringmime(topmime, val)
    if topmime == MIME("text/html")
        return str_repr # |> WebIO.encode_scripts
    elseif topmime in map(MIME, ["image/png", "image/jpeg"])
        return ""
    elseif topmime == MIME("image/svg+xml")
        return str_repr
    else
        return "
$str_repr
"
    end
end

richest_html(::Nothing) = ""

function htmlstring(val)
    sprint() do io
        show(io, MIME"text/html"(), val)
    end
end

htmlstring(val::AbstractString) = val

"""
Wrap an observable in a scope to enable "live updating."

This method also contains distinct code paths for the cases where the observable
contains "non-simple" data types (in particular, observables that contain
Nodes, Scopes, or Widgets need specially handling).
"""
function observable_to_scope(obs::Observable)
    # Create scope that will contain our output observable for rendering.
    scope = Scope()

    # Create output observable (stored rich html representation of input obs).
    # We must create output via the Observable(scope, ...) constructor so that
    # the onjs call below works (the js is attached to the scope, not the
    # observable itself).

    # Do we need separate code paths for Scope/Widget and Node?
    if isa(obs[], Scope) || isa(obs[], AbstractWidget)
        output = Observable(scope, "obs-scope", render(obs[]))
        map!(render, output, obs)
        ensure_sync(scope, "obs-scope")
        scope.dom = node(ObservableNode(output.id, "obs-scope"))
        return scope
    end

    if isa(obs[], Node)
        output = Observable(
            scope,
            "obs-node",
            obs[],
        )
        map!(output, obs) do value
            if !isa(value, Node)
                @warn "A rendered observable (scope.id=$(scope.id), obs.id=$(obs.id)) changed from a WebIO node to $(typeof(value))."
                return nothing
            end
            return value
        end
        ensure_sync(scope, "obs-node")
        scope.dom = node(ObservableNode(output.id, "obs-node"))
        return scope
    end

    output = Observable{AbstractString}(
        scope,
        "obs-output",
        richest_html(obs[])
    )

    # Map new values of obs into the output observable and give it a label
    # within the scope.
    map!(output, obs) do data
        # data -> htmlstring(WebIO.render(data)),
        rich = richest_html(data)
        return rich
    end

    # Avoid not-executing


text:
 __precompile__()

module VoiceConversion

using StatsBase
using StatsFuns
using Distributions
using MelGeneralizedCepstrums
using HDF5, JLD
using Compat

export
    # Voice conversion
    FrameByFrameConverter,
    TrajectoryConverter,
    GMMMapParam,
    GMMMap,
    TrajectoryGMMMap,
    TrajectoryGVGMMMap,
    fvconvert,    # feature vector conversion
    vc,           # voice conversion routine
    ncomponents,  # number of mixture components
    dim,

    # Post filters
    VarianceScaling,
    fvpostf!,
    fvpostf,

    # Alignment
    align,
    align_mcep,

    # Datasets
    Dataset,
    ParallelDataset,
    GVDataset,
    push_delta

for fname in [
              "common",
              "dtw",
              "datasets",
              "align",
              "gmm",
              "gmmmap",
              "diffgmm",
              "trajectory_gmmmap",
              "gv"
    ]
    include(string(fname, ".jl"))
end

end # module


text:
@testset "multivar_predicate" begin
    let
        @syms w a b α::Real β::Real

        r1 = @rule sin(2(~x)) => 2sin(~x)*cos(~x)

        function define_success(rhs)
            function mypredicate(bindings)::Bool
                return bindings[:x] === a
            end

            function success(bindings, n)
                if n == 1
                    return mypredicate(bindings) ? rhs(bindings) : nothing
                else
                    return nothing
                end
            end
            return success
        end

        @eqtest r1(sin(2a)) == 2cos(a)*sin(a)
        @eqtest r1(sin(2a), define_success) == 2cos(a)*sin(a)
        @eqtest r1(sin(2b), define_success) == nothing
    end

    let
        @syms x a b z

        acr = @acrule((~a)^(~u) * (~a)^(~v) => (~a)^(~u + ~v))

        function define_success(rhs)
            function mypredicate(bindings)::Bool
                return bindings[:u] === a
            end

            function success(bindings, n)
                if n == 1
                    return mypredicate(bindings) ? rhs(bindings) : nothing
                else
                    return nothing
                end
            end
            return success
        end

        @eqtest acr(x^a * x^z) == x^(a + z)
        @eqtest acr(x^a * x^z, define_success) == x^(a + z)
        @eqtest acr(x^b * x^z, define_success) == nothing
    end

end


text:
mutable struct BinomialLeapMethod{F1,U} <: SafeLeapAlgorithm
  rates::Vector{Float64}
  total_rate::Float64
  stoichiometry::SparseMatrixCSC{Int,Int}
  leap_formula::F1
  available_reactants::U
end

# information we need:
# - binomial distributions
#     - maximum number of firings for reaction j
#     - sampling probability for reaction j
# - available reactants
# - leap size
#     - tau = max_j (max # of firings for reaction j) / rate_j


text:
@testset "Github Issues" begin
    @testset "Splitpath issue" begin
        # TODO(odow): seet the TODO in src/JuliaFormatter.jl. Remove once
        # JuliaFormatter.jl drops support for Julia 1.0.
        dirs = JuliaFormatter.splitpath(@__DIR__)
        @test length(dirs) > 2
        @test dirs[end] == "test"
        @test occursin("JuliaFormatter", dirs[end-1])
    end

    @testset "issue #137" begin
        str = """
        (
            let x = f() do
                    body
                end
                x
            end for x in xs
        )"""
        str_ = """
        (
               let x = f() do
                       body
                   end
                   x
               end for x in xs
         )"""
        @test fmt(str_) == str

        str = """
        (
            let
                x = f() do
                    body
                end
                x
            end for x in xs
        )"""
        str_ = """
        (
          let
              x = f() do
                  body
              end
              x
          end for x in xs)"""
        @test fmt(str_) == str

        str = """
        let n = try
                ..
            catch
                ..
            end
            ..
        end"""
        @test fmt(str) == str

        str = """
        let n = let
                ..
            end
            ..
        end"""
        @test fmt(str) == str

        str = """
        let n = begin
                ..
            end
            ..
        end"""
        @test fmt(str) == str
    end

    @testset "multiline / issue 139" begin
        str_ = """
        m = match(r\"""
                  (
                      pattern1 |
                      pattern2 |
                      pattern3
                  )
                  \"""x, aaa, str)"""
        str = """
        m = match(
            r\"""
            (
                pattern1 |
                pattern2 |
                pattern3
            )
            \"""x,
            aaa,
            str,
        )"""
        @test fmt(str_) == str

        str_ = """
        m = match(r```
                  (
                      pattern1 |
                      pattern2 |
                      pattern3
                  )
                  ```x, aaa, str)"""
        str = """
        m = match(
            r```
            (
                pattern1 |
                pattern2 |
                pattern3
            )
            ```x,
            aaa,
            str,
        )"""
        @test fmt(str_) == str

        str_ = """
        y = similar([
            1
            2
            3
        ], (4, 5))"""
        str = """
        y = similar(
            [
                1
                2
                3
            ],
            (4, 5),
        )"""
        @test fmt(str_) == str

        str_ = """
        y = similar(T[
            1
            2
            3
        ], (4, 5))"""
        str = """
        y = similar(
            T[
                1
                2
                3
            ],
            (4, 5),
        )"""
        @test fmt(str_) == str
    end

    @testset "issue #150" begin
        str_ = "const SymReg{B,MT} = ArrayReg{B,Basic,MT} where {MT <:AbstractMatrix{Basic}}"
        str = "const SymReg{B,MT} = ArrayReg{B,Basic,MT} where {MT<:AbstractMatrix{Basic}}"
        @test fmt(str_, whitespace_typedefs = false) == str

        str = "const SymReg{B, MT} = ArrayReg{B, Basic, MT} where {MT <: AbstractMatrix{Basic}}"
        @test fmt(str_, whitespace_typedefs = true) == str
    end

    @testset "issue #170 - block within comprehension" begin
        str_ = """
        ys = ( if p1(x)
                 f1(x)
        elseif p2(x)
            f2(x)
        else
            f3(x)
        end for    x in xs)
        """
        str = """
        ys = (
            if p1(x)
                f1(x)
            elseif p2(x)
                f2(x)
            else
                f3(x)
            end for x in xs
        )
        """
        @test fmt(str_) == str

        str = """
        ys = map(xs) do x
            if p1(x)
                f1(x)
            elseif p2(x)
                f2(x)
            else
                f3(x)
            end
        end
        """
        @test fmt(str) == str

        str_ = """
        y1 = Any[if true
            very_very_very_very_very_very_very_very_very_very_very_very_very_very_very_very_long_expr
        end for i in 1:1]"""
        str = """
        y1 = Any[
            if true
                very_very_very_very_very_very_very_very_very_very_very_very_very_very_very_very_long_expr
            end for i = 1:1
        ]"""
        @test fmt(str_) == str
        _, s = run_nest(str_, 100)
        @test s.line_offset == 1

        str_ = """
        y1 = [if true
            very_very_very_very_very_very_very_very_very_very_very_very_very_very_very_very_long_expr
        end for i in 1:1]"""
        str = """
        y1 = [
            if true
                very_very_very_very_very_very_very_very_very_very_very_very_very_very_very_very_long_expr
            end for i = 1:1
        ]"""
        @test fmt(str_) == str
        _, s = run_nest(str_, 100)
        @test s.line_offset == 1

        str_ = """
        y1 = [if true
            short_expr
        end for i in 1:1]"""
        str = """
        y1 = [
            if true
                short_expr
            end for i = 1:1
        ]"""
        @test fmt(str_) == str
        _, s = run_nest(str_, 100)
        @test s.line_offset == 1
    end

    @testset "issue #183" begin
        str_ = """
        function f(args...)

            next!(s.progress;
            # comment
            )
            nothing
        end"""
        str = """
        function f(args...)

            next!(
                s.progress,
                # comment
            )
            nothing
        end"""
        @test fmt(str_) == str
    end

    @testset "issue #189" begin
        str_ = """
    D2 = [
            (b_hat * y - delta_hat[i] * y) * gamma[i] + (b * y_hat - delta[i] * y_hat) *
                                                            gamma_hat[i] + (b_hat - y_hat) *
                                                                           delta[i] + (b - y) *
                                                                                delta_hat[i] - delta[i] * delta_hat[i]
            for i = 1:8
        ]"""
        str = """
        D2 = [
            (b_hat * y - delta_hat[i] * y) * gamma[i] +
            (b * y_hat - delta[i] * y_hat) * gamma_hat[i] +
            (b_hat - y_hat) * delta[i] +
            (b - y) * delta_hat[i] - delta[i] * delta_hat[i] for i = 1:8
        ]"""
        @test fmt(str_) == str
    end

    @testset "issue #193" begin
        str = """
        module Module
        # comment
        end"""
        @test fmt(str) == str

        str = """
        module Module
        # comment
        @test
        # comment
        end"""
        @test fmt(str) == str
    end

    @testset "issue #194" begin
        str_ = """
        function mystr( str::String )
        return SubString( str, 1:
        3 )
        end"""
        str = """
        function mystr(str::String)
            return SubString(str, 1:3)
        end"""
        @test fmt(str_) == str
    end

    @testset "issue #200" begin
        str_ = """
        begin
            f() do
                @info @sprintf \"\"\"
                Δmass   = %.16e\"\"\" abs(weightedsum(Q) - weightedsum(Qe)) / weightedsum(Qe)
            end
        end"""

        # NOTE: this looks slightly off because we're compensating for escaping quotes
        str = """
        begin
            f() do
                @info @sprintf \"\"\"
                Δmass   = %.16e\"\"\" abs(weightedsum(Q) - weightedsum(Qe)) /
                                   weightedsum(Qe)
            end
        end"""
        @test fmt(str_, m = 81) == str
        @test fmt(str, m = 82) == str_
    end

    @testset "issue #202" begin
        str_ = """
        @inline function _make_zop_getvalues(iterators)
            types = map(iterators) do itr
                t =     constructorof(typeof(itr))::Union{Iterators.ProductIterator,CartesianIndices}
                Val(t)
            end
            return function (xs) end
        end"""
        str = """
        @inline function _make_zop_getvalues(iterators)
            types = map(iterators) do itr
                t = constructorof(typeof(itr))::Union{Iterators.ProductIterator,CartesianIndices}
                Val(t)
            end
            return function (xs) end
        end"""
        @test fmt(str_, m = 92) == str

        str_ = """
        @vlplot(
            data = dataset("cars"),
            facet = {row = {field = :Origin, type = :nominal}},
            spec = {
                layer = [
                    {
                        mark = :point,
                        encoding =     {x = {field = :Horsepower}, y = {field = :Miles_per_Gallon}},
                    },
                    {
                        mark = {type = :rule, color = :red},
                        data = {values = [{ref = 10}]},
                        encoding = {y = {field = :ref, type = :quantitative}},
                    },
                ],
            }
        )"""
        str = """
        @vlplot(
            data = dataset("cars"),
            facet = {row = {field = :Origin, type = :nominal}},
            spec = {
                layer = [
                    {
                        mark = :point,
                        encoding = {x = {field = :Horsepower}, y = {field = :Miles_per_Gallon}},
                    },
                    {
                        mark = {type = :rule, color = :red},
                        data = {values = [{ref = 10}]},
                        encoding = {y = {field = :ref, type = :quantitative}},
                    },
                ],
            }
        )"""
        @test fmt(str_, m = 92) == str
    end

    @testset "issue #207" begin
        str_ = """
        @traitfn function predict_ar(m::TGP, p::Int = 3, n::Int = 1; y_past = get_y(m)) where {T,TGP<:AbstractGP{T};IsMultiOutput{TGP}}
        end"""

        str = """
        @traitfn function predict_ar(
            m::TGP,
            p::Int = 3,
            n::Int = 1;
            y_past = get_y(m),
        ) where {T,TGP<:AbstractGP{T};IsMultiOutput{TGP}} end"""
        @test fmt(str_, m = 92) == str

        str = """
        @traitfn function predict_ar(
            m::TGP,
            p::Int = 3,
            n::Int = 1;
            y_past = get_y(m),
        ) where {T, TGP <: AbstractGP{T}; IsMultiOutput{TGP}} end"""
        @test fmt(str_, m = 92, whitespace_typedefs = true) == str

        str_ = """
        @traitfn function predict_ar(m::TGP, p::Int = 3, n::Int = 1; y_past = get_y(m)) where C <: Union{T,TGP<:AbstractGP{T};IsMultiOutput{TGP}}
        end"""

        str = """
        @traitfn function predict_ar(
            m::TGP,
            p::Int = 3,
            n::Int = 1;
            y_past = get_y(m),
        ) where {C<:Union{T,TGP<:AbstractGP{T};IsMultiOutput{TGP}}} end"""
        @test fmt(str_, m = 92) == str

        str = """
        @traitfn function predict_ar(
            m::TGP,
            p::Int = 3,
            n::Int = 1;
            y_past = get_y(m),
        ) where {C <: Union{T, TGP <: AbstractGP{T}; IsMultiOutput{TGP}}} end"""
        @test fmt(str_, m = 92, whitespace_typedefs = true) == str
    end

    @testset "issue #218" begin
        str_ = raw"""
        for MT in GROUP_MANIFOLD_BASIS_DISAMBIGUATION
            eval(quote
                @invoke_maker 1 Manifold get_vector(M::$MT, e::Identity, X, B::VeeOrthogonalBasis)
            end)
        end"""
        str = raw"""
        for MT in GROUP_MANIFOLD_BASIS_DISAMBIGUATION
            eval(
                quote
                    @invoke_maker 1 Manifold get_vector(
                        M::$MT,
                        e::Identity,
                        X,
                        B::VeeOrthogonalBasis,
                    )
                end,
            )
        end"""
        @test fmt1(str_) == str
        @test fmt(str_) == str
    end

    @testset "issue 248" begin
        str_ = """
        var = call(a, @macrocall b)"""
        str = """
        var =
            call(
                a,
                @macrocall b
            )"""
        @test fmt(str_, 4, 1) == str
    end

    @testset "issue 260 - BracesCat" begin
        str = "{1; 2; 3}"
        @test fmt(str, 4, length(str)) == str

        str_ = "{1; 2; 3}"
        str = """
        {
          1;
          2;
          3;
        }"""
        @test fmt(str_, 2, length(str_) - 1) == str
        @test fmt(str, 2, length(str_)) == str_
    end

    @testset "issue 262 - removal of @ in nested macrocall" begin
        str = raw":($(@__MODULE__).@macro)"
        @test fmt(str) == str

        str = raw":($(@__MODULE__).property)"
        @test fmt(str) == str

        str = raw":($(@__MODULE__))"
        @test fmt(str) == str

        str = raw":($(@__MODULE__).@field.macro)"
        @test fmt(str) == str

        str_ = raw":($(@__MODULE__.macro).@field.macro)"
        str = raw":($(__MODULE__.@macro).@field.macro)"
        @test fmt(str_) == str
        @test fmt(str) == str

        str_ = raw"@a.b.c"
        str = raw"a.b.@c"
        @test fmt(str_) == str
        @test fmt(str) == str

        str_ = raw"@a.b.c"
        str = raw"a.b.@c"
        @test fmt(str_) == str
        @test fmt(str) == str

        str = raw"a.@b.c"
        @test fmt(str) == str

        str = raw"a.@b.c.d"
        @test fmt(str) == str

        str = raw"a.b.@c.d"
        @test fmt(str) == str

        str_ = raw"@a.b.c.d"
        str = raw"a.b.c.@d"
        @test fmt(str_) == str
        @test fmt(str) == str
    end

    @testset "issue 264 - `let` empty block body" begin
        str_ = "let; end"
        str = """
        let
        end"""
        @test fmt(str_) == str
    end

    @testset "issue 268 - whitespace around dot op if LHS is number literal" begin
        str = "xs[-5 .<= xs .& xs .<= 5]"
        @test fmt(str) == str
        str_ = "xs[(-5 .<= xs) .& (xs .<= 5)]"
        str = "xs[(-5 .<= xs).&(xs.<=5)]"
        @test fmt(str_) == str
    end

    @testset "issue 277 - flatten ops when no whitespace is allowed" begin
        # previously this would remove |>
        str_ = "get_actions(env)[env |> π.learner |> π.explorer]"
        str = "get_actions(env)[env|>π.learner|>π.explorer]"
        @test fmt(str_) == str
        @test fmt(str, whitespace_ops_in_indices = true) == str_
    end

    @testset "issue 286 - Float32 leading/trailing zeros" begin
        str_ = """
        a = 3.f0
        b = 3f0
        c = 30f0
        d = 30.0f0
        e = 30.123f0
        f = .123f0
        """
        str = """
        a = 3.0f0
        b = 3.0f0
        c = 30.0f0
        d = 30.0f0
        e = 30.123f0
        f = 0.123f0
        """
        @test fmt(str_) == str
    end

    @testset "issue 289 - no spaces/nesting for matrix elements" begin
        str_ = """
        A =  [0. 1 0 0
           -k/Jm -c/Jm k/Jm c/Jm
            0 0 0 1
            f(1,2) c/Ja -k/Ja -c/Ja]
        """

        str = """
        A =
            [
                0.0 1 0 0
                -k/Jm -c/Jm k/Jm c/Jm
                0 0 0 1
                f(1, 2) c/Ja -k/Ja -c/Ja
            ]
        """
        @test fmt(str_, 4, 1) == str
    end

    @testset "issue 317 - infinite recursion" begin
        str = raw"""
        SUITE["manifolds"][name]["tv = 2 * tv1 + 3 * tv2"] = @benchmarkable $tv =
            2 * $tv1 + 3 * $tv2
        """
        @test format_text(str, BlueStyle()) == str
    end

    @testset "issue 324 - bounds error when aligning binary op calls" begin
        # caused by the star operator
        str = """
        θ = eigvals(Matrix([0I(n^2) -I(n^2); P0 P1]), -Matrix([I(n^2) 0I(n^2); 0I(n^2) P2]))
        c = maximum(abs.(θ[(imag.(θ).==0).*(real.(θ).>0)]))
        """
        @test format_text(str, align_assignment = true) == str
    end

    @testset "issue 332" begin
        # this string has a nbsp after 'c'
        # so it should have an additional byte because
        # it's unicode
        str_ = """a = b || c ;
               f("A")"""
        str = """a = b || c;
               f("A")"""
        @test format_text(str_) == str
    end

    @testset "issue 336" begin
        str_ = """
        nzthis = _hessian_slice(d, ex, x, H, obj_factor, nzcount, recovery_tmp_storage, Val{1})::Int
        """
        str = """
        nzthis = _hessian_slice(
            d,
            ex,
            x,
            H,
            obj_factor,
            nzcount,
            recovery_tmp_storage,
            Val{1},
        )::Int
        """
        @test fmt(str_, 4, 80) == str
    end
end


text:
using RecipesBase

## Plug into RecipesBase so plots of expressions treated like plots of functions
@recipe f(::Type{T}, v::T) where {T<:Basic} = lambdify(v)
@recipe f(::Type{S}, ss::S) where {S<:AbstractVector{Basic}} = Function[lambdify(s) for s in ss]


text:
## Turing Complete Interpreter
### A Very Tiny Turing Complete Programming Language defined with denotational semantics

# semantica dalle dispense degano
using Metatheory

import Base.ImmutableDict
Mem = Dict{Symbol, Union{Bool, Int}}

read_mem = @theory begin
	(v::Symbol, σ::Mem) |> if v == :skip σ else σ[v] end
end

@testset "Reading Memory" begin
	ex = :((x), $(Mem(:x => 2)))
	@test true == areequal(read_mem, ex, 2)
end

arithm_rules = @theory begin
	(a + b, σ::Mem) 	=> (a, σ) + (b, σ)
	(a * b, σ::Mem) 	=> (a, σ) * (b, σ)
	(a - b, σ::Mem) 	=> (a, σ) - (b, σ)
	(n::Int, σ::Mem) 	=> n
	(a::Int + b::Int) 	|> a + b
	(a::Int * b::Int) 	|> a * b
	(a::Int - b::Int) 	|> a - b
end

@testset "Arithmetic" begin
	@test areequal(read_mem ∪ arithm_rules, :((2 + 3), $(Mem())), 5)
end

# don't need to access memory
bool_rules = @theory begin
	(a < b, σ::Mem) => (a, σ) < (b, σ)
	(a ∨ b, σ::Mem) => (a, σ) ∨ (b, σ)
	(a ∧ b, σ::Mem) => (a, σ) ∧ (b, σ)
	(¬(a), σ::Mem) => ¬((a, σ))

	(a::Bool, σ::Mem) 	|> a
	(¬a::Bool) 			|> !a
	(a::Bool ∨ b::Bool) 	|> (a || b)
	(a::Bool ∧ b::Bool) 	|> (a && b)
	(a::Int < b::Int) 	|> (a < b)
end

t = read_mem ∪ arithm_rules ∪ bool_rules

@testset "Booleans" begin
	@test areequal(t, :((false ∨ false), $(Mem())), false)

	exx = :((false ∨ false) ∨ ¬(false ∨ false), $(Mem(:x => 2)))
	g = EGraph(exx)
	saturate!(g, t)
	ex = extract!(g, astsize)
	@test ex == true
	params=SaturationParams(timeout=12)
	@test areequal(t, exx, true; params=params)

	@test areequal(t, :((2 < 3) ∧ (3 < 4), $(Mem(:x => 2))), true)
	@test areequal(t, :((2 < x) ∨ ¬(3 < 4), $(Mem(:x => 2))), false)
	@test  areequal(t, :((2 < x) ∨ ¬(3 < 4), $(Mem(:x => 4))), true)
end

if_rules = @theory begin
	(if guard; t end) => (if guard; t else :skip end)
	(if guard; t else f end, σ::Mem) => (if (guard, σ); t else f end, σ)
	(if true; t else f end, σ::Mem) => (t, σ)
	(if false; t else f end, σ::Mem) => (f, σ)
end

if_language = read_mem ∪ arithm_rules ∪ bool_rules ∪ if_rules

@testset "If Semantics" begin
	@test areequal(if_language, 2, :(if true x else 0 end, 	$(Mem(:x => 2))))
	@test areequal(if_language, 0, :(if false x else 0 end, 	$(Mem(:x => 2))))
	@test areequal(if_language, 2, :(if ¬(false) x else 0 end, $(Mem(:x => 2))))
	params=SaturationParams(timeout=10)
	@test areequal(if_language, 0, :(if ¬(2 < x) x else 0 end, $(Mem(:x => 3))); params=params)
end

while_rules = @theory begin
	(:skip, σ::Mem) => σ
	((c1; c2), σ::Mem) => ((c1, σ); c2)
	(c1::Int; c2) => c2
	(c1::Bool; c2) => c2
	(σ::Mem; c2) => (c2, σ)
	(while guard body end, σ::Mem) =>
		(if guard; (body; while guard body end) else :skip end, σ)
end

write_mem = @theory begin
	(sym::Symbol = val, σ::Mem) => (sym = (val, σ), σ)
	(sym::Symbol = val::Int, σ::Mem) |> merge(σ, Dict(sym => val))
end

while_language = if_language ∪ write_mem ∪ while_rules;

@testset "While Semantics" begin
	exx = :((x = 3), $(Mem(:x => 2)))
	g = EGraph(exx)
	saturate!(g, while_language)
	ex = extract!(g, astsize)

	@test areequal(while_language, Mem(:x => 3), exx)

	exx = :((x = 4; x = x + 1), $(Mem(:x => 3)))
	g = EGraph(exx)
	saturate!(g, while_language)
	ex = extract!(g, astsize)

	params=SaturationParams(timeout=10)
	@test areequal(while_language, Mem(:x => 5), exx; params=params)

	# FIXME bug!
	params=SaturationParams(timeout=14)
	exx = :((if x < 10 x = x + 1 else skip end), $(Mem(:x => 3)))
	# g = EGraph(exx)
	# params=SaturationParams(timeout=12)
	# saturate!(g, while_language, params)
	# println(extract!(g, astsize))
	@test areequal(while_language, Mem(:x => 4), exx; params=params)
	# exit(0)

	exx = :(
		(while x < 10
			x = x + 1
		end; x), $(Mem(:x => 3)))
	g = EGraph(exx)
	params=SaturationParams(timeout=100, scheduler=Schedulers.SimpleScheduler)
	saturate!(g, while_language, params)
	@test 10 == extract!(g, astsize)
end


text:
@testset "JLSOFile" begin

    withenv("JLSO_IMAGE" => "busybox") do
        jlso = JLSOFile(:data => "the image env variable is set")
        @test jlso.image == "busybox"
    end

    # Reset the cached image for future tests
    JLSO._CACHE[:IMAGE] = ""

    @testset "$(fmt): $k" for fmt in (:bson, :julia_serialize), (k, v) in datas
        jlso = JLSOFile(k => v; format=fmt, compression=:none)
        io = IOBuffer()
        bytes = fmt === :bson ? bson(io, Dict("object" => v)) : serialize(io, v)
        expected = take!(io)

        @test jlso.objects[k] == expected
    end

    @testset "kwarg constructor" begin
        jlso = JLSOFile(; a=collect(1:10), b="hello")
        @test jlso[:b] == "hello"
        if VERSION < v"1.7"
            @test haskey(jlso.manifest, "BSON")
        else
            @test haskey(jlso.manifest["deps"], "BSON")
        end
    end

    @testset "no-arg constructor" begin
        jlso = JLSOFile()
        @test jlso isa JLSOFile
        @test isempty(jlso.objects)
        if VERSION < v"1.7"
            @test haskey(jlso.manifest, "BSON")
        else
            @test haskey(jlso.manifest["deps"], "BSON")
        end
    end
end

@testset "unknown format" begin
    @test_throws(
        LOGGER,
        MethodError,
        JLSOFile("String" => "Hello World!", format=:unknown)
    )
end

@testset "show" begin
    jlso = JLSOFile(:string => datas[:String])
    expected = string(
        "JLSOFile([data]; version=v\"2.0.0\", julia=v\"$VERSION\", ",
        "format=:julia_serialize, image=\"\")"
    )
    @test sprint(show, jlso) == sprint(print, jlso)
end

@testset "activate" begin
    jlso = JLSOFile(:string => datas[:String])
    mktempdir() do d
        Pkg.activate(jlso, d) do
            @show Base.active_project()
        end

        @test ispath(joinpath(d, "Project.toml"))
        @test ispath(joinpath(d, "Manifest.toml"))

        @test Pkg.TOML.parsefile(joinpath(d, "Project.toml")) == jlso.project
        @test Pkg.TOML.parsefile(joinpath(d, "Manifest.toml")) == jlso.manifest
    end
end

@testset "keys/haskey" begin
    jlso = JLSOFile(:string => datas[:String])
    @test collect(keys(jlso)) == [:string]
    @test haskey(jlso, :string)
    @test !haskey(jlso, :other)
end

@testset "get/get!" begin
    v = datas[:String]
    jlso = JLSOFile(:str => v)
    @test get(jlso, :str, "fail") == v
    @test get!(jlso, :str, "fail") == v

    @test get(jlso, :other, v) == v
    @test !haskey(jlso, :other)

    @test get!(jlso, :other, v) == v
    @test jlso[:other] == v

    # key must be a Symbol
    @test get(jlso, "str", 999) == 999
    @test_throws MethodError get!(jlso, "str", 999)
end

max_stars_repo_name: invenia/JLSO

max_issues_repo_name: invenia/JLSO

max_forks_repo_name: invenia/JLSO


text:
using LinearAlgebra
using SparseArrays
using jInv.Mesh
using Multigrid
using Test

println("************************************************* AMG Setup For DivSigGrad ******************************************************");

domain = [0.0, 1.0, 0.0, 1.0];
n = [50,50];
Mr = getRegularMesh(domain,n)
m  = exp.(randn(prod(n)));
Ar = getNodalDivSigGradMatrix(Mr,vec(m));
Ar = Ar + 1e-8*norm(Ar,1)*sparse(1.0I, size(Ar,2),size(Ar,2));

levels      = 3;
numCores 	= 2;
maxIter     = 5;
relativeTol = 1e-4;
relaxType   = "SPAI";
relaxParam  = 1.0;
relaxPre 	= 1;
relaxPost   = 1;
cycleType   ='V';
coarseSolveType = "Julia";

MG = getMGparam(Float64,Int64,levels,numCores,maxIter,relativeTol,relaxType,relaxParam,relaxPre,relaxPost,cycleType,coarseSolveType);

N = size(Ar,2);

b = Ar*rand(N,3);
b = b/norm(b);
x = zeros(N,3);
SA_AMGsetup(Ar,MG,true,size(b,2),true);

println("****************************** Stand-alone SA-AMG: ******************************")
solveMG(MG,b,x,true);
@test norm(Ar*x - b) < 0.01;

println("****************************** CG preconditioned with SA-AMG: ******************************")
x[:] .= 0.0;
solveCG_MG(Ar,MG,b,x,true)
@test norm(Ar*x - b) < 0.005;

println("****************************** BiCGSTAB preconditioned with SA-AMG: ******************************")
x[:] .= 0.0;
solveBiCGSTAB_MG(Ar,MG,b,x,true);
@test norm(Ar*x - b) < 0.005;

println("****************************** transposeHierarchy *********************************");
transposeHierarchy(MG,true);

println("****************************** replaceMatrixHierarchy *********************************");
replaceMatrixInHierarchy(MG,Ar,true);

println("****************************** GMRES preconditioned with SA-AMG: (only one rhs...) ******************************")
x[:] .= 0.0;
b1 = vec(b[:,1]);
x1 = zeros(N);
MG = getMGparam(Float64,Int64,levels,numCores,maxIter,relativeTol,relaxType,relaxParam,relaxPre,relaxPost,cycleType,"GMRES");
SA_AMGsetup(Ar,MG,true,1,true);
solveGMRES_MG(Ar,MG,b1,x1,true,2,true)
@test norm(Ar*x1 - b1) < 0.001;
println("****************************** Classical AMG *********************************");
x1 = zeros(N);
println("Stand-alone Classical AMG one rhs:")
MG = getMGparam(Float64,Int64,levels,numCores,maxIter,relativeTol,relaxType,relaxParam,relaxPre,relaxPost,cycleType,coarseSolveType);
ClassicalAMGsetup(Ar,MG,true,size(b1,2),true);
solveMG(MG,b1,x1,true);
println("CG preconditioned with C-AMG in 3D:")
x1[:] .= 0.0;
solveCG_MG(Ar,MG,b1,x1,true)
@test norm(Ar*x1 - b1) < 0.005;

println("****************************** GMRES preconditioned with SA-AMG: (multiple rhs...) ******************************")

MG = getMGparam(Float64,Int64,levels,numCores,maxIter,relativeTol,"Jac-GMRES",relaxParam,relaxPre,relaxPost,'K',coarseSolveType);
SA_AMGsetup(Ar,MG,true,size(b,2),true);
x[:] .= 0.0;
solveGMRES_MG(Ar,MG,b,x,true,2,true)

Ar = 0;
b = 0;
x = 0;
Mr = 0;

println("************************************************* Example AMG 3D ******************************************************");

domain = [0.0, 1.0, 0.0, 1.0, 0.0, 1.0];
n = [32,32,16];
Mr = getRegularMesh(domain,n)
m  = exp.(randn(prod(n)));
Ar = getNodalDivSigGradMatrix(Mr,vec(m));
Ar = Ar + 1e-6*norm(Ar,1)*sparse(1.0I, size(Ar,2),size(Ar,2));

MG = getMGparam(Float64,Int64,levels,numCores,maxIter,relativeTol,relaxType,relaxParam,relaxPre,relaxPost,cycleType,"Julia");

N = size(Ar,2);
b = Ar*rand(N,3);
b = b/norm(b);

x = zeros(N,3);
println("Stand-alone SA-AMG:")
SA_AMGsetup(Ar,MG,true,size(b,2),true);
solveMG(MG,b,x,true);
@test norm(Ar*x - b) < 0.005;
println("CG preconditioned with SA-AMG in 3D:")
x[:] .= 0.0;
solveCG_MG(Ar,MG,b,x,true)
@test norm(Ar*x - b) < 0.005;

x = zeros(N,3);
println("Stand-alone Classical AMG:")
MG = getMGparam(Float64,Int64,levels,numCores,maxIter,relativeTol,relaxType,relaxParam,relaxPre,relaxPost,cycleType,"GMRES");
ClassicalAMGsetup(Ar,MG,true,size(b,2),true);
solveMG(MG,b,x,true);
@test norm(Ar*x - b) < 0.005;
println("CG preconditioned with C-AMG in 3D:")
x[:] .= 0.0;
solveCG_MG(Ar,MG,b,x,true)
@test norm(Ar*x - b) < 0.005;

Ar = 0;
b = 0;
x = 0;
Mr = 0;


text:
using Test
import DPFEHM

ns = [5, 10, 20]
mins = [0, 0, 0]
maxs = [1, 1, 1]

#test that the 2d grid first goes over the the y coordinate, then the x coordinate
coords, neighbors, areasoverlengths, volumes = DPFEHM.regulargrid2d(mins[1:2], maxs[1:2], ns[1:2], maxs[3] - mins[3])
@test coords[2, 1] != coords[2, 2]
@test coords[2, 1] == coords[2, 1 + ns[2]]
@test coords[1, 1] != coords[1, 1 + ns[2]]
@test all(coords[1, 1] .== coords[1, 1:ns[2]])
@test sum(volumes) ≈ prod(ns[1:2]) * prod((maxs[1:2] .- mins[1:2]) ./ (ns[1:2] .- 1)) * (maxs[3] - mins[3])#it goes a little higher than just prod(maxs - mins), because the center of the boundary cells are on the boundary, and half of those cells are outside the boundary

#test that the 3d grid first goes over the z coordinate, then the y coordinate, then the x coordinate
coords, neighbors, areasoverlengths, volumes = DPFEHM.regulargrid3d(mins, maxs, ns)
@test coords[3, 1] != coords[3, 2]
@test coords[3, 1] == coords[3, 1 + ns[3]]
@test coords[2, 1] != coords[2, 1 + ns[3]]
@test all(coords[2, 1] .== coords[2, 1:ns[3]])
@test coords[1, 1] != coords[1, 1 + ns[3] * ns[2]]
@test all(coords[1, 1] .== coords[1, 1:ns[2] * ns[3]])
@test sum(volumes) ≈ prod(ns) * prod((maxs .- mins) ./ (ns .- 1))#it goes a little higher than just prod(maxs - mins), because the center of the boundary cells are on the boundary, and half of those cells are outside the boundary

#test that it is computing the darcy velocity correctly for a constant pressure gradient along each of the coordinate axes
Ks = ones(reverse(ns)...)
h = [x for z = 1:ns[3], y = 1:ns[2], x = 1:ns[1]]
dv = DPFEHM.darcy_velocity(h, Ks, mins, maxs, ns)
@test dv(0.5, 0.5, 0.5) ≈ [ns[1] - 1, 0, 0]
h = [y for z = 1:ns[3], y = 1:ns[2], x = 1:ns[1]]
dv = DPFEHM.darcy_velocity(h, Ks, mins, maxs, ns)
@test dv(0.5, 0.5, 0.5) ≈ [0, ns[2] - 1, 0]
h = [z for z = 1:ns[3], y = 1:ns[2], x = 1:ns[1]]
dv = DPFEHM.darcy_velocity(h, Ks, mins, maxs, ns)
@test dv(0.5, 0.5, 0.5) ≈ [0, 0, ns[3] - 1]