text:
using Distributed
addprocs(2)
println("There are $(nprocs()) processes")

@everywhere begin
    using Pkg
    Pkg.activate("downstream")
    Pkg.develop(PackageSpec(path=joinpath(pwd(), "..")))
    Pkg.instantiate()
    using OrdinaryDiffEq
    prob = ODEProblem((u,p,t)->1.01u,0.5,(0.0,1.0))
    u0s = [rand()*prob.u0 for i in 1:2]
    function prob_func(prob,i,repeat)
        println("Running trajectory $i")
        ODEProblem(prob.f,u0s[i],prob.tspan)
    end
end

ensemble_prob = EnsembleProblem(prob, prob_func=prob_func)
sim = solve(ensemble_prob,Tsit5(),EnsembleSplitThreads(),trajectories=2)

@everywhere function lorenz!(du,u,p,t)
    du[1] = 10.0*(u[2]-u[1])
    du[2] = u[1]*(28.0-u[3]) - u[2]
    du[3] = u[1]*u[2] - (8/3)*u[3]
end

u0 = [1.0, 0.0, 0.0]
tspan = (0.0, 100.0)
p = [1, 2.0, 3]
prob = ODEProblem(lorenz!, u0, tspan, p)

@everywhere function prob_func(prob,i,repeat)
    prob = remake(prob, tspan=(rand(), 100.0), p=rand(3))
    return prob
end

ensemble_prob = EnsembleProblem(prob, prob_func = prob_func, safetycopy=true)

println("Running EnsembleSerial()")
@test length(solve(ensemble_prob, Tsit5(), EnsembleSerial(), trajectories=100)) == 100
println("Running EnsembleThreads()")
@test length(solve(ensemble_prob, Tsit5(), EnsembleThreads(), trajectories=100)) == 100
println("Running EnsembleDistributed()")
@test length(solve(ensemble_prob, Tsit5(), EnsembleDistributed(), trajectories=100)) == 100
println("Running EnsembleSplitThreads()")
@test length(solve(ensemble_prob, Tsit5(), EnsembleSplitThreads(), trajectories=100)) == 100

ext: jl


text:
mutable struct PIPSNLPData
    d::Union{JuMP.NLPEvaluator,Nothing}
    n::Int
    m::Int
    local_m::Int                    #rank row dimension
    jacnnz::Int                     #Jacobian num nonzeros
    hessnnz::Int                    #Hessian num nonzeros
    firstIeq::Vector{Int}           #first stage equality constraint jacobian nonzero rows
    firstJeq::Vector{Int}           #first stage equality constraint jacobian nonzero cols
    firstVeq::Vector{Float64}       #first stage equality constraint jacobian values
    secondIeq::Vector{Int}          #second stage equality constraint jacobian nonzero rows
    secondJeq::Vector{Int}          #second stage equality constraint jacobian nonzero cols
    secondVeq::Vector{Float64}      #second stage equality constraint jacobian values
    firstIineq::Vector{Int}         #first stage inequality constraint jacobian nonzero rows
    firstJineq::Vector{Int}
    firstVineq::Vector{Float64}
    secondIineq::Vector{Int}
    secondJineq::Vector{Int}
    secondVineq::Vector{Float64}
    num_eqconnect::Int              #num equalities connecting first to second stage
    num_ineqconnect::Int            #num inequalities connecting first to second stage
    eqconnect_lb::Vector{Float64}   #connection equality lower bounds
    eqconnect_ub::Vector{Float64}   #connection equality upper bounds
    ineqconnect_lb::Vector{Float64} #connection inequality lower bounds
    ineqconnect_ub::Vector{Float64} #connection inequality upper bounds
    eq_idx::Vector{Int}             #equality connect indices
    ineq_idx::Vector{Int}           #inequality connect indices
    firstJeqmat
    secondJeqmat
    firstJineqmat
    secondJineqmat
    linkIeq::Vector{Int}            #equality link jacobian rows
    linkJeq::Vector{Int}            #equality link jacobian cols
    linkVeq::Vector{Float64}        #equality link jacobian values
    linkIineq::Vector{Int}
    linkJineq::Vector{Int}
    linkVineq::Vector{Float64}
    x_sol::Vector{Float64}          #solution on rank
    coreid::Int                     #rank
    loaded::Bool
    local_unsym_hessnnz::Int        #local hessian num nonzeros
end
PIPSNLPData() = PIPSNLPData(nothing,0,0,0,0,0,Int[],Int[], Float64[], Int[], Int[], Float64[],Int[],Int[],Float64[], Int[], Int[],
                            Float64[], 0, 0, Float64[],Float64[],Float64[],Float64[], Int[], Int[],
                            nothing, nothing, nothing, nothing, Int[],Int[], Float64[], Int[], Int[], Float64[],Float64[], 0, false, 0)

#data that each worker needs
mutable struct PipsNLPWorkerData
	first_stage::OptiNode
    n_linkeq_cons::Int64
    n_linkineq_cons::Int64
    link_ineq_lower::Vector{Float64}
    link_ineq_upper::Vector{Float64}
    link_eq_lower::Vector{Float64}
    link_eq_upper::Vector{Float64}
end

function _get_pips_data(node::OptiNode)
    return node.ext[:pips_data]
end

function _setup_pips_nlp_data!(graph::OptiGraph)

	if !(has_subgraphs(graph))
		#the first stage is empty
		first_stage = OptiNode()
		submodels = all_nodes(graph)
	elseif length(graph.subgraphs) == 1
		num_nodes(graph) == 1 || error("Incompatible optigraph structure.  Optigraph contains a single subgraph, but has multiple optinodes in the highest level graph.")
		first_stage = getnode(graph,1)
		submodels = all_nodes(graph.subgraphs[1])
	end

    model_list = [first_stage;submodels]
	for (idx,node) in enumerate(model_list)
		node.ext[:pips_data] = PIPSNLPData()
		constraint_data = PipsNLP.PipsNlpInterface.get_constraint_data(node)
		node.ext[:constraint_data] = constraint_data
	end

	first_stage_data = _get_pips_data(first_stage)
	n_sub_models = length(submodels)
	model_list = [first_stage; submodels]

	#We assume every worker has the link constraint structure that connects across nodes
	#TODO: use MPI.ALLGATHER! to setup linkconstraints when they are different across workers?
	#NOTE: Each rank would need to know which links are equivalent e.g. using some sort of user-defined index.  For now,
	#we assume the user needs to model their problem correctly using ghost optinodes, or they use distribute_optigraph, which sets up
	#links correctly.  In the future, distributed optigraphs will make communicating the structure easier.

	#TODO: Check whether workers can use zeros for bounds on linkconstraints they don't need.

	#IDENTIFY LINK CONSTRAINTS
	link_connect_eq,link_connect_ineq,linkeqconstraints,linkineqconstraints  = _identify_linkconstraints(graph)

	#FIRST STAGE -> SECOND STAGE LINK CONSTRAINTS
	for link in link_connect_eq
		node = _get_subnode(first_stage,link)::OptiNode
		local_data = _get_pips_data(node)
		firstIeq = local_data.firstIeq
		firstJeq = local_data.firstJeq
		firstVeq = local_data.firstVeq
		secondIeq = local_data.secondIeq
		secondJeq = local_data.secondJeq
		secondVeq = local_data.secondVeq
		push!(local_data.eqconnect_lb, link.set.value)
		push!(local_data.eqconnect_ub, link.set.value)
		local_data.num_eqconnect += 1
		row = local_data.num_eqconnect
		for (var,coeff) in link.func.terms
			if getnode(var) == first_stage
				push!(firstIeq, row)
				push!(firstJeq, var.index.value)
				push!(firstVeq, coeff)
			else
				@assert getnode(var) == node
				push!(secondIeq, row)
				push!(secondJeq, var.index.value)
				push!(secondVeq, coeff)
			end
		end
	end

	for link in link_connect_ineq
		node = _get_subnode(first_stage,link)
		firstIineq = local_data.firstIineq
		firstJineq = local_data.firstJineq
		firstVineq = local_data.firstVineq
		secondIineq = local_data.secondIineq
		secondJineq = local_data.secondJineq
		secondVineq = local_data.secondVineq
		if isa(link.set,MOI.LessThan)
			push!(local_data.ineqconnect_lb, -Inf)
			push!(local_data.ineqconnect_ub, link.set.upper)
		elseif isa(link.set,MOI.GreaterThan)
			push!(local_data.ineqconnect_lb, link.set.lower)
			push!(local_data.ineqconnect_ub, Inf)
		elseif isa(link.set,MOI.Interval)
			push!(local_data.ineqconnect_lb, link.set.lower)
			push!(local_data.ineqconnect_ub, link.set.upper)
		end
		local_data.num_ineqconnect += 1
		row = local_data.num_ineqconnect
		for (var,coeff) in link.func.terms
			if getnode(var) == first_stage
				push!(firstIineq, row)
				push!(firstJineq, var.index.value)
				push!(firstVineq, ind)
			else
				@assert getnode(var) == node
				push!(secondIineq, row)
				push!(secondJineq, var.index.value)
				push!(secondVineq, ind)
			end
		end
	end

	#SECOND STAGE LINK CONSTRAINTS
	nlinkeq = length(linkeqconstraints)
	nlinkineq = length(linkineqconstraints)
	ineqlink_lb = zeros(nlinkineq)
	ineqlink_ub = zeros(nlinkineq)
	eqlink_lb = zeros(nlinkeq)
	eqlink_ub = zeros(nlinkeq)

	#INEQUALITY CONSTRAINTS
	#Inequality bounds
	#NOTE: linkconstraint bounds need to be on every rank.
	for (idx,link) in enumerate(linkineqconstraints)
		row = idx
		#Bounds
		if isa(link.set,MOI.LessThan)
			ineqlink_lb[row] = -Inf
			ineqlink_ub[row] = link.set.upper
		elseif isa(link.set,MOI.GreaterThan)
			ineqlink_lb[row] = link.set.lower
			ineqlink_ub[row] = Inf
		elseif isa(link.set,MOI.Interval)
			ineqlink_lb[row] = link.set.lower
			ineqlink_ub[row] = link.set.upper
		end
		#Inequality values
		for (var,coeff) in link.func.terms
			node = getnode(var)
			local_data = _get_pips_data(node)
			linkIineq = local_data.linkIineq
			linkJineq = local_data.linkJineq
			linkVineq = local_data.linkVineq
			push!(linkIineq, row)
			push!(linkJineq, var.index.value)
			push!(linkVineq, coeff)
		end
	end

	#EQUALITY CONSTRAINTS
	for (idx,link) in enumerate(linkeqconstraints)
		row = idx
		#Bounds
		eqlink_lb[row] = link.set.value
		eqlink_ub[row] = link.set.value
		#Equality values
		for (var,coeff) in link.func.terms
			node = getnode(var)
			local_data = _get_pips_data(node)
			linkIeq = local_data.linkIeq
			linkJeq = local_data.linkJeq
			linkVeq = local_data.linkVeq
			push!(linkIeq, row)              #the variable row
			push!(linkJeq, var.index.value)  #the variable column
			push!(linkVeq, coeff)            #the coefficient
		end
	end

    worker_data = PipsNLPWorkerData(
    first_stage,
    nlinkeq,
    nlinkineq,
    ineqlink_lb,
    ineqlink_ub,
    eqlink_lb,
    eqlink_ub)

	return worker_data
end

#categorize link constraints into hierarchical and distributed
function _identify_linkconstraints(graph::OptiGraph)
    links_connect_eq = LinkConstraint[]
    links_connect_ineq = LinkConstraint[]
    links_eq = LinkConstraint[]
    links_ineq = LinkConstraint[]

	if !(has_subgraphs(graph))
		for link in all_linkconstraints(graph)
			link_constraint = constraint_object(link)
			if isa(link_constraint.set,MOI.EqualTo)
				push!(links_eq,link_constraint)
			else
				push!(links_ineq,link_constraint)
			end
		end
	else
        first_stage_nodes = getnodes(graph)
	    sub_nodes = setdiff(all_nodes(graph),first_stage_nodes)
	    for link in all_linkconstraints(graph)
	        link_constraint = constraint_object(link)
	        link_nodes = getnodes(link_constraint)

	        #if first_stage_node in link_nodes
			if isempty(setdiff(first_stage_nodes,link_nodes))
	            @assert length(link_nodes) == 2
	            if isa(link_constraint.set,MOI.EqualTo)
	                push!(links_connect_eq,link_constraint)
	            else
	                push!(links_connect_ineq,link_constraint)
	            end
	        else
	            if isa(link_constraint.set,MOI.EqualTo)
	                push!(links_eq,link_constraint)
	            else
	                push!(links_ineq,link_constraint)
	            end
	        end
	    end
	end
    return links_connect_eq,links_connect_ineq,links_eq,links_ineq
end


text:
@testset "Physics" begin

@testset "Testing function definitions" begin
    @test isdefined(GRASS, :absorption_line)
    @test isdefined(GRASS, :gaussian_line)
    @test isdefined(GRASS, :lorentzian_line)
    @test isdefined(GRASS, :width_thermal)
    @test isdefined(GRASS, :quad_limb_darkening)
    @test isdefined(GRASS, :rotation_period)
    @test isdefined(GRASS, :patch_velocity_los)
end

@testset "Testing limb darkening" begin
    # disk center
    t1 = (0.0, 0.0)
    u1 = 0.4
    u2 = 0.26
    mu = GRASS.calc_mu(t1)
    @test GRASS.quad_limb_darkening(mu, u1, u2) == 1.0

    # off disk
    t2 = (1.25, 1.25)
    mu = GRASS.calc_mu(t2)
    @test GRASS.quad_limb_darkening(mu, u1, u2) == 0.0

    # no LD
    u1 = 0.0
    u2 = 0.0
    mu = range(0.1, 1.0, length=10)
    @test all(GRASS.quad_limb_darkening.(mu, u1, u2) .== 1.0)
end

@testset "Testing rotational velocity" begin
    # default pole
    @test GRASS.patch_velocity_los(0.0, 0.0) == 0.0
    @test GRASS.patch_velocity_los(1.0, 0.0) < 0.0
    @test GRASS.patch_velocity_los(-1.0, 0.0) > 0.0

    # pole at face
    pole = (0.0, 0.0, 1.0)
    @test GRASS.patch_velocity_los(0.0, 0.0, pole=pole) == 0.0
    @test GRASS.patch_velocity_los(1.0, 0.0, pole=pole) == 0.0
    @test GRASS.patch_velocity_los(-1.0, 0.0, pole=pole) == 0.0
end

@testset "Testing Gaussian line" begin
    xs = range(5434.0, 5435.0, step=0.01)
    ys = GRASS.gaussian_line.(xs, mid=5434.5, width=0.1, depth=1.0)
    @test minimum(ys) ≈ 0.0 atol=1e-5
    @test maximum(ys) ≈ 1.0 atol=1e-5

    ys = GRASS.gaussian_line.(xs, mid=5434.5, width=0.1, depth=0.5)
    @test minimum(ys) ≈ 0.5 atol=1e-5
    @test maximum(ys) ≈ 1.0 atol=1e-5
end

end

max_stars_repo_name: palumbom/GRASS

max_issues_repo_name: palumbom/GRASS

max_forks_repo_name: palumbom/GRASS


text:
function regress(interpreter::Interpreter,
                 domain::Domain, state::State, action::Action, args;
                 check::Bool=true, fail_mode::Symbol=:error)
    # Check whether action is relevant
    if check && !relevant(interpreter, domain, state, action, args)
        if fail_mode == :no_op return state end
        error("Effect $(get_effect(action)) is not relevant.")
    end
    return regress!(interpreter, domain, copy(state), action, args; check=false)
end

function regress!(interpreter::Interpreter,
                  domain::Domain, state::State, action::Action, args;
                  check::Bool=true, fail_mode::Symbol=:error)
    # Check whether action is relevant
    if check && !relevant(interpreter, domain, state, action, args)
        if fail_mode == :no_op return state end
        error("Effect $(get_effect(action)) is not relevant.")
    end
    subst = Subst(var => val for (var, val) in zip(get_argvars(action), args))
    precond = substitute(get_precond(action), subst)
    effect = substitute(get_effect(action), subst)
    # Compute regression difference as Precond - Additions
    # TODO: Handle conditional effects, disjunctive preconditions, etc.
    pre_diff = precond_diff(domain, state, precond)
    eff_diff = effect_diff(domain, state, effect)
    append!(pre_diff.del, eff_diff.add)
    return update!(interpreter, domain, state, pre_diff)
end


text:
function cyclic(n)
    @var z[1:n]
    eqs = [sum(prod(z[(k-1)%n+1] for k = j:j+m) for j = 1:n) for m = 0:(n-2)]
    push!(eqs, prod(z) - 1)
    System(eqs, z)
end

function bacillus()
    @var w w2 w2v v w2v2 vP sigmaB w2sigmaB vPp phos
    poly = [
        (-1 * 0.7 * w + -2 * 3600.0 * (w^2 / 2) + 2 * 18.0 * w2) * (0.2 + sigmaB) +
        4.0 * 0.4 * (1 + 30.0sigmaB),
        -1 * 0.7 * w2 +
        3600.0 * (w^2 / 2) +
        -1 * 18.0 * w2 +
        -1 * 3600.0 * w2 * v +
        18.0w2v +
        36.0w2v +
        -1 * 3600.0 * w2 * sigmaB +
        18.0w2sigmaB,
        -1 * 0.7 * w2v +
        3600.0 * w2 * v +
        -1 * 18.0 * w2v +
        -1 * 3600.0 * w2v * v +
        18.0w2v2 +
        -1 * 36.0 * w2v +
        36.0w2v2 +
        1800.0 * w2sigmaB * v +
        -1 * 1800.0 * w2v * sigmaB,
        (
            -1 * 0.7 * v +
            -1 * 3600.0 * w2 * v +
            18.0w2v +
            -1 * 3600.0 * w2v * v +
            18.0w2v2 +
            -1 * 1800.0 * w2sigmaB * v +
            1800.0 * w2v * sigmaB +
            180.0vPp
        ) * (0.2 + sigmaB) + 4.5 * 0.4 * (1 + 30.0sigmaB),
        -1 * 0.7 * w2v2 + 3600.0 * w2v * v + -1 * 18.0 * w2v2 + -1 * 36.0 * w2v2,
        -1 * 0.7 * vP + 36.0w2v + 36.0w2v2 + -1 * 3600.0 * vP * phos + 18.0vPp,
        (
            -1 * 0.7 * sigmaB +
            -1 * 3600.0 * w2 * sigmaB +
            18.0w2sigmaB +
            1800.0 * w2sigmaB * v +
            -1 * 1800.0 * w2v * sigmaB
        ) * (0.2 + sigmaB) + 0.4 * (1 + 30.0sigmaB),
        -1 * 0.7 * w2sigmaB +
        3600.0 * w2 * sigmaB +
        -1 * 18.0 * w2sigmaB +
        -1 * 1800.0 * w2sigmaB * v +
        1800.0 * w2v * sigmaB,
        -1 * 0.7 * vPp + 3600.0 * vP * phos + -1 * 18.0 * vPp + -1 * 180.0 * vPp,
        (phos + vPp) - 2.0,
    ]
    System(poly, [w, w2, w2v, v, w2v2, vP, sigmaB, w2sigmaB, vPp, phos])
end

function minors()
    @var x y z
    # This are the maximal minors of a 3 by 5 matrix
    exprs = [
        -(19 / 2) * x^6 - (306679 / 75600) * x^5 * y +
        (3566281 / 453600) * x^4 * y^2 +
        (2458541 / 113400) * x^3 * y^3 +
        (1057993 / 97200) * x^2 * y^4 - (4093 / 28350) * x * y^5 - (701 / 675) * y^6 +
        (156769 / 4704) * x^5 * z +
        (2136401 / 201600) * x^4 * y * z - (485549 / 317520) * x^3 * y^2 * z -
        (1228333069 / 38102400) * x^2 * y^3 * z - (444121 / 12960) * x * y^4 * z -
        (12491 / 2016) * y^5 * z + (5554901 / 1411200) * x^4 * z^2 -
        (1396343 / 453600) * x^3 * y * z^2 - (2815352089 / 50803200) * x^2 * y^2 * z^2 -
        (6727757957 / 127008000) * x * y^3 * z^2 - (4014469 / 378000) * y^4 * z^2 +
        (21830581 / 658560) * x^3 * z^3 +
        (208824821 / 14817600) * x^2 * y * z^3 - (105675343 / 9525600) * x * y^2 * z^3 -
        (744971 / 158760) * y^3 * z^3 +
        (4234567 / 352800) * x^2 * z^4 +
        (26384189 / 2116800) * x * y * z^4 - (153973 / 952560) * y^2 * z^4 +
        (1465123 / 176400) * x * z^5 +
        (537751 / 238140) * y * z^5 +
        (121 / 315) * z^6 +
        (160717 / 37800) * x^5 +
        (3742241 / 378000) * x^4 * y +
        (28611931 / 1134000) * x^3 * y^2 +
        (61601 / 36450) * x^2 * y^3 - (17656943 / 1701000) * x * y^4 -
        (20531 / 5400) * y^5 - (894547 / 84672) * x^4 * z -
        (2346079 / 470400) * x^3 * y * z - (16933223 / 396900) * x^2 * y^2 * z -
        (553452937 / 13608000) * x * y^3 * z - (206119 / 34020) * y^4 * z -
        (848783 / 705600) * x^3 * z^2 - (61574861 / 2540160) * x^2 * y * z^2 -
        (24620545003 / 457228800) * x * y^2 * z^2 - (92970251 / 9072000) * y^3 * z^2 -
        (1267237 / 423360) * x^2 * z^3 - (474156533 / 114307200) * x * y * z^3 -
        (30023869 / 6804000) * y^2 * z^3 +
        (186664 / 19845) * x * z^4 +
        (3533393 / 3175200) * y * z^4 +
        (2557 / 14175) * z^5 - (433 / 540) * x^4 +
        (3979387 / 1134000) * x^3 * y +
        (176181763 / 8164800) * x^2 * y^2 - (23612089 / 3402000) * x * y^3 -
        (1099939 / 97200) * y^4 - (3776473 / 423360) * x^3 * z -
        (57493058237 / 1143072000) * x^2 * y * z -
        (5521110767 / 54432000) * x * y^2 * z + (29454571 / 2041200) * y^3 * z -
        (313762747 / 6804000) * x^2 * z^2 - (995824699 / 38102400) * x * y * z^2 -
        (123893137 / 4665600) * y^2 * z^2 +
        (17925923 / 2381400) * x * z^3 +
        (5874637 / 1587600) * y * z^3 - (813 / 140) * z^4 - (5369 / 48600) * x^3 -
        (1596487 / 453600) * x^2 * y - (419513 / 32400) * x * y^2 -
        (113359 / 21600) * y^3 - (87364253 / 2721600) * x^2 * z -
        (246098857 / 10886400) * x * y * z + (104373149 / 3628800) * y^2 * z -
        (63677 / 3024) * x * z^2 - (583877 / 388800) * y * z^2 - (7757 / 3150) * z^3 +
        (14549 / 6480) * x^2 - (517663 / 75600) * x * y + (150341 / 8640) * y^2 -
        (10412771 / 302400) * x * z +
        (209327 / 181440) * y * z +
        (5506 / 945) * z^2 - (211 / 60) * x - (517 / 480) * y +
        (299 / 30) * z +
        47 / 16,
        (2 / 3) * x^6 + (220411 / 37800) * x^5 * y - (494933 / 113400) * x^4 * y^2 +
        (4298527 / 453600) * x^3 * y^3 +
        (794509 / 226800) * x^2 * y^4 +
        (4049 / 7200) * x * y^5 +
        (21 / 80) * y^6 - (4987 / 1575) * x^5 * z - (1746887 / 75600) * x^4 * y * z +
        (21906547 / 1360800) * x^3 * y^2 * z +
        (6301643 / 2041200) * x^2 * y^3 * z +
        (43001477 / 2721600) * x * y^4 * z +
        (47231 / 12960) * y^5 * z - (1141 / 2400) * x^4 * z^2 +
        (23865827 / 12700800) * x^3 * y * z^2 - (17041681 / 1190700) * x^2 * y^2 * z^2 +
        (55657337 / 23814000) * x * y^3 * z^2 +
        (306001 / 75600) * y^4 * z^2 +
        (4657139 / 470400) * x^3 * z^3 - (11497061 / 793800) * x^2 * y * z^3 +
        (3973489 / 6096384) * x * y^2 * z^3 +
        (4007189 / 1693440) * y^3 * z^3 +
        (389769 / 156800) * x^2 * z^4 +
        (98207447 / 19051200) * x * y * z^4 +
        (165523 / 79380) * y^2 * z^4 +
        (1689101 / 352800) * x * z^5 +
        (12493993 / 6350400) * y * z^5 +
        (2327 / 3600) * z^6 +
        (1129 / 450) * x^5 - (8123971 / 1134000) * x^4 * y +
        (5767 / 2268) * x^3 * y^2 +
        (30397 / 170100) * x^2 * y^3 +
        (745733 / 151200) * x * y^4 +
        (6509 / 2400) * y^5 - (63893 / 8400) * x^4 * z +
        (7910717 / 3175200) * x^3 * y * z - (312932947 / 9525600) * x^2 * y^2 * z +
        (126425483 / 9072000) * x * y^3 * z +
        (7200563 / 1360800) * y^4 * z - (33705709 / 2116800) * x^3 * z^2 -
        (729523649 / 12700800) * x^2 * y * z^2 -
        (1374377629 / 47628000) * x * y^2 * z^2 + (28807441 / 28576800) * y^3 * z^2 -
        (52739317 / 6350400) * x^2 * z^3 - (34461347 / 88905600) * x * y * z^3 +
        (64953487 / 91445760) * y^2 * z^3 - (10755347 / 6350400) * x * z^4 +
        (6332959 / 2540160) * y * z^4 +
        (161069 / 529200) * z^5 +
        (317 / 108) * x^4 +
        (36225143 / 1360800) * x^3 * y - (737491 / 129600) * x^2 * y^2 +
        (2872463 / 259200) * x * y^3 +
        (270253 / 129600) * y^4 - (1155239 / 64800) * x^3 * z -
        (91503641 / 3175200) * x^2 * y * z + (19841081 / 1270080) * x * y^2 * z -
        (18974807 / 5443200) * y^3 * z - (43379321 / 6350400) * x^2 * z^2 -
        (6952763359 / 381024000) * x * y * z^2 - (319039243 / 38102400) * y^2 * z^2 -
        (21814313 / 2646000) * x * z^3 - (497789 / 5443200) * y * z^3 -
        (25273 / 132300) * z^4 + (193469 / 10800) * x^3 -
        (12734041 / 2268000) * x^2 * y +
        (264791 / 64800) * x * y^2 +
        (525217 / 86400) * y^3 +
        (279821 / 31500) * x^2 * z +
        (108020863 / 10584000) * x * y * z - (19662193 / 907200) * y^2 * z +
        (19939433 / 3528000) * x * z^2 - (108970333 / 25401600) * y * z^2 -
        (32677 / 70560) * z^3 +
        (4453 / 10800) * x^2 +
        (2576759 / 100800) * x * y - (27121 / 2880) * y^2 +
        (59197 / 1400) * x * z +
        (1015583 / 302400) * y * z - (3773701 / 302400) * z^2 +
        (15217 / 720) * x +
        (1379 / 960) * y +
        (15359 / 5040) * z +
        223 / 16,
        (29 / 7) * x^6 +
        (23227 / 2520) * x^5 * y +
        (917771 / 63000) * x^4 * y^2 +
        (11548 / 945) * x^3 * y^3 - (515479 / 14400) * x^2 * y^4 -
        (25571 / 2592) * x * y^5 + (497 / 1080) * y^6 - (130651 / 8400) * x^5 * z -
        (21666527 / 504000) * x^4 * y * z - (6470071 / 1296000) * x^3 * y^2 * z -
        (4645133 / 403200) * x^2 * y^3 * z +
        (18227527 / 259200) * x * y^4 * z +
        (26387 / 1440) * y^5 * z +
        (456553 / 235200) * x^4 * z^2 - (178328503 / 10584000) * x^3 * y * z^2 +
        (5343791 / 846720) * x^2 * y^2 * z^2 +
        (114270983 / 1360800) * x * y^3 * z^2 +
        (12508007 / 453600) * y^4 * z^2 - (778381 / 78400) * x^3 * z^3 -
        (42594247 / 2646000) * x^2 * y * z^3 +
        (19687831 / 1036800) * x * y^2 * z^3 +
        (178310191 / 10886400) * y^3 * z^3 - (34949 / 7350) * x^2 * z^4 -
        (11336141 / 756000) * x * y * z^4 - (1164427 / 1701000) * y^2 * z^4 -
        (5001 / 800) * x * z^5 - (1089971 / 302400) * y * z^5 - (1039 / 1050) * z^6 +
        (3533 / 6300) * x^5 - (1921583 / 189000) * x^4 * y -
        (1492481 / 113400) * x^3 * y^2 - (180587 / 6048) * x^2 * y^3 +
        (966361 / 194400) * x * y^4 +
        (10747 / 2592) * y^5 - (11137493 / 352800) * x^4 * z -
        (641388757 / 15876000) * x^3 * y * z +
        (3712801 / 907200) * x^2 * y^2 * z +
        (601009153 / 4665600) * x * y^3 * z +
        (36446017 / 907200) * y^4 * z - (213277 / 14112) * x^3 * z^2 -
        (40330637 / 1134000) * x^2 * y * z^2 +
        (21871050209 / 571536000) * x * y^2 * z^2 +
        (29266459 / 907200) * y^3 * z^2 - (236459003 / 24696000) * x^2 * z^3 -
        (283084717 / 21168000) * x * y * z^3 + (681321199 / 254016000) * y^2 * z^3 -
        (325699 / 117600) * x * z^4 - (151170469 / 31752000) * y * z^4 -
        (49801 / 110250) * z^5 +
        (100679 / 3780) * x^4 +
        (10379 / 1296) * x^3 * y - (21673403 / 972000) * x^2 * y^2 +
        (4933739 / 194400) * x * y^3 - (684991 / 25920) * y^4 +
        (1054363 / 42336) * x^3 * z +
        (442296401 / 190512000) * x^2 * y * z +
        (727316399 / 27216000) * x * y^2 * z +
        (30058829 / 907200) * y^3 * z +
        (8978639 / 264600) * x^2 * z^2 +
        (1330320163 / 19051200) * x * y * z^2 +
        (6313 / 4032) * y^2 * z^2 +
        (210615067 / 10584000) * x * z^3 +
        (10294511 / 3024000) * y * z^3 +
        (1093333 / 441000) * z^4 +
        (77383 / 9450) * x^3 +
        (507847 / 40500) * x^2 * y +
        (52943 / 10800) * x * y^2 +
        (92641 / 12960) * y^3 +
        (60435509 / 5292000) * x^2 * z +
        (517865647 / 9072000) * x * y * z +
        (2229677 / 907200) * y^2 * z +
        (142591051 / 5292000) * x * z^2 +
        (4027 / 960) * y * z^2 +
        (26357 / 22050) * z^3 +
        (126251 / 9450) * x^2 +
        (17051 / 2700) * x * y - (134689 / 1440) * y^2 - (770923 / 151200) * x * z +
        (913747 / 75600) * y * z +
        (12853 / 2800) * z^2 - (1999 / 360) * x + (172 / 45) * y - (904 / 105) * z -
        209 / 8,
        -(10 / 3) * x^6 - (1153 / 504) * x^5 * y +
        (889 / 60) * x^4 * y^2 +
        (2360791 / 151200) * x^3 * y^3 +
        (10061 / 27000) * x^2 * y^4 - (509737 / 25200) * x * y^5 - (5917 / 450) * y^6 +
        (105163 / 8820) * x^5 * z - (9789209 / 2116800) * x^4 * y * z -
        (35670337 / 4233600) * x^3 * y^2 * z - (33335147 / 2646000) * x^2 * y^3 * z -
        (9131057 / 378000) * x * y^4 * z +
        (139639 / 25200) * y^5 * z +
        (192649 / 117600) * x^4 * z^2 +
        (8835811 / 846720) * x^3 * y * z^2 - (35411371 / 1587600) * x^2 * y^2 * z^2 -
        (48672359 / 21168000) * x * y^3 * z^2 +
        (80051 / 10800) * y^4 * z^2 +
        (356102 / 33075) * x^3 * z^3 - (8227481 / 3628800) * x^2 * y * z^3 +
        (63073973 / 15240960) * x * y^2 * z^3 - (197 / 75600) * y^3 * z^3 -
        (1484633 / 1058400) * x^2 * z^4 + (9985471 / 1016064) * x * y * z^4 -
        (12687223 / 1814400) * y^2 * z^4 + (150233 / 75600) * x * z^5 -
        (166951 / 37800) * y * z^5 - (20471 / 10800) * z^6 - (1852 / 945) * x^5 -
        (1189 / 11340) * x^4 * y +
        (839971 / 75600) * x^3 * y^2 +
        (2531383 / 324000) * x^2 * y^3 - (329003 / 25200) * x * y^4 - (227 / 20) * y^5 -
        (419001 / 39200) * x^4 * z - (274927927 / 13608000) * x^3 * y * z -
        (294227203 / 28576800) * x^2 * y^2 * z +
        (46111091 / 3528000) * x * y^3 * z +
        (21033239 / 453600) * y^4 * z +
        (26793971 / 11113200) * x^3 * z^2 - (1884015269 / 111132000) * x^2 * y * z^2 +
        (180244423 / 47628000) * x * y^2 * z^2 +
        (27223717 / 1270080) * y^3 * z^2 +
        (2575621 / 3175200) * x^2 * z^3 +
        (38388299 / 2540160) * x * y * z^3 +
        (41228297 / 4762800) * y^2 * z^3 +
        (413839 / 105840) * x * z^4 - (13710967 / 1411200) * y * z^4 -
        (154151 / 75600) * z^5 +
        (2465143 / 113400) * x^4 +
        (314857 / 162000) * x^3 * y +
        (43081193 / 6804000) * x^2 * y^2 - (33754031 / 907200) * x * y^3 -
        (507833 / 16200) * y^4 - (1463879 / 198450) * x^3 * z -
        (165791273 / 7938000) * x^2 * y * z - (369287183 / 6804000) * x * y^2 * z +
        (107246053 / 2721600) * y^3 * z +
        (55352327 / 4762800) * x^2 * z^2 +
        (718293743 / 63504000) * x * y * z^2 +
        (3794783 / 226800) * y^2 * z^2 +
        (26560811 / 3175200) * x * z^3 +
        (100644931 / 12700800) * y * z^3 +
        (266431 / 75600) * z^4 - (38197 / 17010) * x^3 + (3512851 / 1134000) * x^2 * y -
        (63289 / 1680) * x * y^2 + (1841 / 150) * y^3 - (18876703 / 1587600) * x^2 * z -
        (586122937 / 47628000) * x * y * z + (1279709 / 453600) * y^2 * z -
        (3089503 / 1270080) * x * z^2 +
        (15865867 / 1058400) * y * z^2 +
        (40233 / 5600) * z^3 - (256829 / 56700) * x^2 - (465601 / 21600) * x * y -
        (82267 / 1440) * y^2 - (1713977 / 50400) * x * z +
        (3229721 / 453600) * y * z +
        (7303 / 2240) * z^2 - (13469 / 720) * x + (383 / 240) * y - (11999 / 1008) * z -
        285 / 16,
        -(489 / 40) * x^6 - (77719 / 15120) * x^5 * y +
        (267929 / 31752) * x^4 * y^2 +
        (14208101 / 1587600) * x^3 * y^3 +
        (3097037 / 793800) * x^2 * y^4 - (2249 / 9450) * x * y^5 - (83 / 450) * y^6 +
        (11783 / 315) * x^5 * z - (380857 / 30240) * x^4 * y * z +
        (326371 / 21600) * x^3 * y^2 * z - (2740261 / 635040) * x^2 * y^3 * z -
        (731603 / 105840) * x * y^4 * z - (139 / 600) * y^5 * z +
        (156029 / 10080) * x^4 * z^2 +
        (19124683 / 635040) * x^3 * y * z^2 - (226728533 / 88905600) * x^2 * y^2 * z^2 +
        (39730727 / 9525600) * x * y^3 * z^2 +
        (5820139 / 2268000) * y^4 * z^2 +
        (3469769 / 66150) * x^3 * z^3 +
        (17835733 / 1190700) * x^2 * y * z^3 - (461024351 / 88905600) * x * y^2 * z^3 +
        (13947103 / 6350400) * y^3 * z^3 +
        (19490131 / 1058400) * x^2 * z^4 +
        (52893409 / 2381400) * x * y * z^4 +
        (357967 / 119070) * y^2 * z^4 +
        (3239909 / 264600) * x * z^5 +
        (22693901 / 4762800) * y * z^5 +
        (187 / 140) * z^6 +
        (11333 / 2700) * x^5 +
        (12163517 / 680400) * x^4 * y +
        (1459812847 / 71442000) * x^3 * y^2 +
        (373783 / 45360) * x^2 * y^3 - (39203 / 12960) * x * y^4 -
        (23771 / 18000) * y^5 - (114607 / 4725) * x^4 * z -
        (13299883 / 302400) * x^3 * y * z - (276735289 / 4762800) * x^2 * y^2 * z -
        (14432231 / 362880) * x * y^3 * z - (19137691 / 4536000) * y^4 * z -
        (380501 / 21600) * x^3 * z^2 - (260695451 / 3175200) * x^2 * y * z^2 -
        (630612179 / 9525600) * x * y^2 * z^2 - (216511 / 18144) * y^3 * z^2 -
        (1435591 / 635040) * x^2 * z^3 - (121932491 / 4762800) * x * y * z^3 -
        (58163557 / 12700800) * y^2 * z^3 +
        (2120249 / 529200) * x * z^4 +
        (267097 / 529200) * y * z^4 +
        (2167 / 3240) * z^5 +
        (11459 / 12960) * x^4 +
        (1289 / 189) * x^3 * y +
        (25911149 / 2381400) * x^2 * y^2 - (47656529 / 2721600) * x * y^3 -
        (17533 / 2700) * y^4 - (29615893 / 453600) * x^3 * z -
        (1438245331 / 13608000) * x^2 * y * z -
        (13010618809 / 142884000) * x * y^2 * z - (5065363 / 1814400) * y^3 * z -
        (112568063 / 2268000) * x^2 * z^2 - (14883632387 / 142884000) * x * y * z^2 -
        (15060991 / 362880) * y^2 * z^2 - (20820743 / 1270080) * x * z^3 -
        (14009717 / 2116800) * y * z^3 - (1723 / 1512) * z^4 +
        (93439 / 32400) * x^3 +
        (350801 / 100800) * x^2 * y - (15809887 / 604800) * x * y^2 - (535 / 48) * y^3 -
        (185643803 / 2721600) * x^2 * z - (45604103 / 403200) * x * y * z +
        (1052449 / 201600) * y^2 * z - (13916659 / 272160) * x * z^2 -
        (10998347 / 302400) * y * z^2 - (81359 / 10080) * z^3 + (1457 / 4320) * x^2 -
        (21497929 / 1209600) * x * y - (2233 / 240) * y^2 -
        (54671779 / 907200) * x * z + (1877693 / 80640) * y * z - (6997 / 504) * z^2 -
        (1701 / 320) * x +
        (9867 / 1280) * y +
        (7409 / 320) * z +
        837 / 64,
        -(561 / 140) * x^6 + (460223 / 58800) * x^5 * y - (7650443 / 352800) * x^4 * y^2 - (1992937 / 98000) * x^3 * y^3 -
        (39587 / 7875) * x^2 * y^4 - (28771 / 11340) * x * y^5 + (301 / 675) * y^6 -
        (37523 / 3360) * x^5 * z - (38302421 / 4233600) * x^4 * y * z +
        (55466501 / 705600) * x^3 * y^2 * z - (23967241 / 3175200) * x^2 * y^3 * z +
        (708251 / 113400) * x * y^4 * z +
        (297571 / 37800) * y^5 * z +
        (907777 / 117600) * x^4 * z^2 +
        (396988313 / 14817600) * x^3 * y * z^2 +
        (4007575717 / 63504000) * x^2 * y^2 * z^2 +
        (202935527 / 1587600) * x * y^3 * z^2 +
        (15421667 / 453600) * y^4 * z^2 - (1030669 / 70560) * x^3 * z^3 +
        (150426089 / 74088000) * x^2 * y * z^3 +
        (506748581 / 9072000) * x * y^2 * z^3 +
        (295785457 / 10584000) * y^3 * z^3 - (2894573 / 352800) * x^2 * z^4 -
        (1558981 / 117600) * x * y * z^4 + (373682 / 99225) * y^2 * z^4 -
        (52453 / 6300) * x * z^5 - (2875571 / 529200) * y * z^5 - (176 / 105) * z^6 -
        (55739 / 6300) * x^5 +
        (1501811 / 396900) * x^4 * y +
        (19103233 / 756000) * x^3 * y^2 +
        (105225241 / 4536000) * x^2 * y^3 +
        (973957 / 272160) * x * y^4 +
        (40121 / 16200) * y^5 - (5148797 / 151200) * x^4 * z -
        (139394203 / 1411200) * x^3 * y * z - (4800962771 / 381024000) * x^2 * y^2 * z +
        (35899531 / 4536000) * x * y^3 * z +
        (12601129 / 907200) * y^4 * z - (1605559 / 529200) * x^3 * z^2 -
        (28872007 / 992250) * x^2 * y * z^2 +
        (152989777 / 10584000) * x * y^2 * z^2 +
        (33709733 / 1360800) * y^3 * z^2 - (16073177 / 4233600) * x^2 * z^3 -
        (1299883783 / 190512000) * x * y * z^3 + (53429123 / 27216000) * y^2 * z^3 -
        (12347 / 31500) * x * z^4 - (828467 / 588000) * y * z^4 - (862 / 1575) * z^5 +
        (90577 / 15120) * x^4 +
        (15040213 / 6350400) * x^3 * y +
        (86912731 / 3265920) * x^2 * y^2 - (72777953 / 2721600) * x * y^3 +
        (115429 / 21600) * y^4 - (13141267 / 604800) * x^3 * z -
        (138795047659 / 2286144000) * x^2 * y * z -
        (1068020413 / 27216000) * x * y^2 * z +
        (43272007 / 1814400) * y^3 * z +
        (301622123 / 63504000) * x^2 * z^2 +
        (574755439 / 63504000) * x * y * z^2 +
        (3645611 / 259200) * y^2 * z^2 - (572167 / 181440) * x * z^3 -
        (38930821 / 3175200) * y * z^3 - (1913 / 900) * z^4 - (1726633 / 680400) * x^3 -
        (1261033 / 113400) * x^2 * y +
        (2136649 / 181440) * x * y^2 +
        (240227 / 7200) * y^3 - (22151083 / 5443200) * x^2 * z -
        (331428401 / 7257600) * x * y * z - (971543 / 201600) * y^2 * z +
        (5269 / 20160) * x * z^2 +
        (2209661 / 50400) * y * z^2 +
        (5588 / 1575) * z^3 +
        (82753 / 9072) * x^2 - (12559 / 151200) * x * y + (16127 / 1920) * y^2 -
        (13297943 / 604800) * x * z +
        (140491 / 10080) * y * z +
        (59687 / 1680) * z^2 +
        (11329 / 1440) * x - (71 / 80) * y - (23 / 80) * z - 333 / 32,
        -(43 / 8) * x^6 +
        (36529 / 5040) * x^5 * y +
        (5549 / 1890) * x^4 * y^2 +
        (432449 / 264600) * x^3 * y^3 - (357649 / 19845) * x^2 * y^4 -
        (53753 / 6750) * x * y^5 - (368 / 375) * y^6 + (1454093 / 141120) * x^5 * z -
        (2570719 / 529200) * x^4 * y * z - (1231037 / 231525) * x^3 * y^2 * z +
        (2643317 / 1587600) * x^2 * y^3 * z - (190277 / 11760) * x * y^4 * z -
        (51 / 1400) * y^5 * z - (534629 / 28224) * x^4 * z^2 +
        (2511667 / 176400) * x^3 * y * z^2 - (1103696501 / 44452800) * x^2 * y^2 * z^2 +
        (8494807 / 907200) * x * y^3 * z^2 +
        (1069823 / 31500) * y^4 * z^2 - (24217 / 33075) * x^3 * z^3 -
        (49530571 / 2116800) * x^2 * y * z^3 - (24642463 / 6350400) * x * y^2 * z^3 -
        (126629 / 75600) * y^3 * z^3 - (19461481 / 1058400) * x^2 * z^4 +
        (915533 / 3175200) * x * y * z^4 - (12689 / 1440) * y^2 * z^4 -
        (669097 / 226800) * x * z^5 - (601199 / 56700) * y * z^5 - (2057 / 540) * z^6 -
        (63557 / 7560) * x^5 - (6826709 / 793800) * x^4 * y +
        (17433953 / 1058400) * x^3 * y^2 +
        (5719663 / 1190700) * x^2 * y^3 +
        (17566163 / 1134000) * x * y^4 +
        (8767 / 1125) * y^5 - (5298761 / 1058400) * x^4 * z -
        (677579281 / 63504000) * x^3 * y * z +
        (399157901 / 95256000) * x^2 * y^2 * z +
        (505720987 / 31752000) * x * y^3 * z +
        (196961 / 15120) * y^4 * z +
        (221983 / 33600) * x^3 * z^2 +
        (167343703 / 7620480) * x^2 * y * z^2 +
        (1697257 / 47628) * x * y^2 * z^2 +
        (9103099 / 151200) * y^3 * z^2 - (2548711 / 1905120) * x^2 * z^3 +
        (160965253 / 5080320) * x * y * z^3 +
        (4294153 / 151200) * y^2 * z^3 - (7729 / 36288) * x * z^4 +
        (881 / 2160) * y * z^4 - (3581 / 1080) * z^5 +
        (1219501 / 272160) * x^4 +
        (385892197 / 57153600) * x^3 * y +
        (498801683 / 9525600) * x^2 * y^2 +
        (4714901 / 453600) * x * y^3 +
        (2053 / 180) * y^4 +
        (41170111 / 15240960) * x^3 * z - (329042429 / 10584000) * x^2 * y * z +
        (15214099 / 31752000) * x * y^2 * z +
        (2161913 / 45360) * y^3 * z +
        (196918811 / 8467200) * x^2 * z^2 +
        (711856841 / 15240960) * x * y * z^2 +
        (7240981 / 181440) * y^2 * z^2 +
        (2180833 / 113400) * x * z^3 +
        (47030719 / 1814400) * y * z^3 +
        (5569 / 1080) * z^4 +
        (368041 / 90720) * x^3 - (17094137 / 2721600) * x^2 * y +
        (20063299 / 907200) * x * y^2 +
        (7937 / 360) * y^3 - (4969597 / 907200) * x^2 * z -
        (16176199 / 1814400) * x * y * z +
        (53349 / 11200) * y^2 * z +
        (4098677 / 181440) * x * z^2 +
        (7060141 / 181440) * y * z^2 +
        (4529 / 288) * z^3 +
        (3439 / 1260) * x^2 - (1352627 / 201600) * x * y + (707 / 120) * y^2 -
        (3149401 / 362880) * x * z - (5095 / 2688) * y * z +
        (1523 / 72) * z^2 +
        (839 / 960) * x - (603 / 160) * y - (469 / 64) * z - 531 / 64,
        -(101 / 20) * x^6 - (126451 / 12600) * x^5 * y - (52649 / 2100) * x^4 * y^2 +
        (53153 / 31500) * x^3 * y^3 +
        (692599 / 75600) * x^2 * y^4 +
        (11197 / 5400) * x * y^5 - (7 / 36) * y^6 +
        (903 / 50) * x^5 * z +
        (3502991 / 75600) * x^4 * y * z - (50046163 / 2268000) * x^3 * y^2 * z -
        (680647 / 33600) * x^2 * y^3 * z - (17982049 / 680400) * x * y^4 * z -
        (37271 / 6480) * y^5 * z - (9647 / 50400) * x^4 * z^2 +
        (189563599 / 4233600) * x^3 * y * z^2 +
        (737187221 / 19051200) * x^2 * y^2 * z^2 -
        (174213217 / 6350400) * x * y^3 * z^2 - (4309649 / 453600) * y^4 * z^2 +
        (723553 / 50400) * x^3 * z^3 +
        (144053197 / 3528000) * x^2 * y * z^3 +
        (117944923 / 6804000) * x * y^2 * z^3 - (42815051 / 9525600) * y^3 * z^3 +
        (81071 / 10080) * x^2 * z^4 +
        (6978289 / 302400) * x * y * z^4 +
        (205894669 / 47628000) * y^2 * z^4 +
        (27813 / 5600) * x * z^5 +
        (3633277 / 1058400) * y * z^5 +
        (521 / 840) * z^6 +
        (19 / 36) * x^5 +
        (3344137 / 226800) * x^4 * y +
        (1452887 / 75600) * x^3 * y^2 +
        (35347231 / 1360800) * x^2 * y^3 +
        (5733463 / 907200) * x * y^4 - (10369 / 4320) * y^5 +
        (101089 / 2520) * x^4 * z +
        (11028331 / 317520) * x^3 * y * z - (54153679 / 9525600) * x^2 * y^2 * z -
        (438043297 / 4762800) * x * y^3 * z - (66679 / 2100) * y^4 * z +
        (118071 / 5600) * x^3 * z^2 +
        (1549714217 / 27216000) * x^2 * y * z^2 -
        (16287117607 / 381024000) * x * y^2 * z^2 - (613672261 / 21168000) * y^3 * z^2 +
        (3515423 / 220500) * x^2 * z^3 - (11743427 / 15876000) * x * y * z^3 -
        (89419703 / 10584000) * y^2 * z^3 +
        (2229499 / 470400) * x * z^4 +
        (781651 / 508032) * y * z^4 +
        (132877 / 176400) * z^5 - (64421 / 2160) * x^4 -
        (11103277 / 1360800) * x^3 * y + (3832951 / 453600) * x^2 * y^2 -
        (17641 / 189000) * x * y^3 - (18203 / 14400) * y^4 -
        (2501867 / 45360) * x^3 * z - (260838553 / 3175200) * x^2 * y * z -
        (5295087649 / 31752000) * x * y^2 * z - (111827287 / 1512000) * y^3 * z -
        (18942281 / 362880) * x^2 * z^2 - (6610196443 / 50803200) * x * y * z^2 -
        (271403081 / 5292000) * y^2 * z^2 - (34141027 / 846720) * x * z^3 -
        (31854803 / 1411200) * y * z^3 - (3474301 / 705600) * z^4 -
        (4121 / 2160) * x^3 + (2152697 / 129600) * x^2 * y -
        (9040487 / 907200) * x * y^2 + (235 / 72) * y^3 - (2816389 / 75600) * x^2 * z -
        (22377529 / 254016) * x * y * z - (6072571 / 100800) * y^2 * z -
        (20856991 / 604800) * x * z^2 - (184951531 / 5644800) * y * z^2 -
        (4480477 / 705600) * z^3 + (20401 / 4320) * x^2 - (62371 / 4200) * x * y +
        (7809 / 320) * y^2 - (1585813 / 151200) * x * z - (250301 / 4800) * y * z -
        (126821 / 9600) * z^2 +
        (3607 / 180) * x +
        (21371 / 320) * y +
        (557 / 112) * z +
        1071 / 16,
        (14 / 3) * x^6 + (2329 / 420) * x^5 * y - (699059 / 26460) * x^4 * y^2 +
        (3691 / 264600) * x^3 * y^3 - (899279 / 132300) * x^2 * y^4 +
        (7692 / 875) * x * y^5 +
        (31 / 12) * y^6 - (6137 / 420) * x^5 * z +
        (77009 / 8400) * x^4 * y * z +
        (39870697 / 3175200) * x^3 * y^2 * z - (104196787 / 4762800) * x^2 * y^3 * z +
        (22768199 / 3969000) * x * y^4 * z - (18019 / 1800) * y^5 * z -
        (546629 / 100800) * x^4 * z^2 - (740849 / 90720) * x^3 * y * z^2 -
        (59576507 / 9525600) * x^2 * y^2 * z^2 + (10123541 / 1058400) * x * y^3 * z^2 -
        (68791 / 10800) * y^4 * z^2 - (4633133 / 302400) * x^3 * z^3 -
        (74929 / 33600) * x^2 * y * z^3 + (79300723 / 9525600) * x * y^2 * z^3 -
        (69497 / 20160) * y^3 * z^3 - (231319 / 56700) * x^2 * z^4 -
        (531493 / 100800) * x * y * z^4 + (139271 / 90720) * y^2 * z^4 -
        (9563 / 3360) * x * z^5 +
        (5479 / 12096) * y * z^5 +
        (391 / 2160) * z^6 +
        (1853 / 360) * x^5 - (129821 / 113400) * x^4 * y -
        (5434711 / 1587600) * x^3 * y^2 +
        (206561 / 64800) * x^2 * y^3 +
        (101627 / 3780) * x * y^4 +
        (5923 / 600) * y^5 +
        (3780493 / 226800) * x^4 * z +
        (1125769457 / 23814000) * x^3 * y * z +
        (234841721 / 83349000) * x^2 * y^2 * z +
        (4883177 / 226800) * x * y^3 * z - (3355691 / 151200) * y^4 * z +
        (11462071 / 3175200) * x^3 * z^2 +
        (22573379 / 604800) * x^2 * y * z^2 +
        (567793459 / 13891500) * x * y^2 * z^2 - (22574443 / 1587600) * y^3 * z^2 +
        (31502839 / 6350400) * x^2 * z^3 +
        (709098209 / 38102400) * x * y * z^3 - (10823321 / 6350400) * y^2 * z^3 +
        (293119 / 635040) * x * z^4 +
        (46029653 / 5080320) * y * z^4 +
        (2951 / 1728) * z^5 - (767033 / 32400) * x^4 + (86393 / 75600) * x^3 * y -
        (52617421 / 1984500) * x^2 * y^2 +
        (1785493 / 45360) * x * y^3 +
        (11999 / 720) * y^4 +
        (284311 / 16200) * x^3 * z +
        (53483 / 3175200) * x^2 * y * z +
        (184785661 / 11907000) * x * y^2 * z - (10881973 / 151200) * y^3 * z -
        (362885699 / 25401600) * x^2 * z^2 + (46670083 / 2116800) * x * y * z^2 -
        (67828421 / 2116800) * y^2 * z^2 - (2004203 / 25401600) * x * z^3 -
        (103485733 / 5080320) * y * z^3 - (134527 / 30240) * z^4 +
        (148793 / 12960) * x^3 +
        (1782961 / 113400) * x^2 * y +
        (1557827 / 30240) * x * y^2 +
        (10981 / 480) * y^3 +
        (44419097 / 2721600) * x^2 * z +
        (246738319 / 6350400) * x * y * z - (137771 / 2240) * y^2 * z +
        (968323 / 396900) * x * z^2 - (2361599 / 84672) * y * z^2 -
        (224557 / 20160) * z^3 +
        (1199449 / 43200) * x^2 +
        (128017 / 2240) * x * y +
        (1739 / 80) * y^2 +
        (34483621 / 907200) * x * z - (1129 / 63) * y * z - (214547 / 40320) * z^2 +
        (2599 / 60) * x +
        (1035 / 32) * y +
        (1161 / 80) * z +
        639 / 16,
        (15 / 4) * x^6 - (7811 / 5880) * x^5 * y +
        (90403 / 29400) * x^4 * y^2 +
        (180401 / 11025) * x^3 * y^3 - (572491 / 12600) * x^2 * y^4 +
        (694333 / 37800) * x * y^5 - (469 / 90) * y^6 - (1399 / 1120) * x^5 * z -
        (310577 / 117600) * x^4 * y * z - (93897719 / 2646000) * x^3 * y^2 * z +
        (8777939 / 1764000) * x^2 * y^3 * z +
        (575129 / 21600) * x * y^4 * z - (205259 / 4200) * y^5 * z +
        (301643 / 50400) * x^4 * z^2 +
        (48892999 / 10584000) * x^3 * y * z^2 - (15110183 / 882000) * x^2 * y^2 * z^2 -
        (65666327 / 3969000) * x * y^3 * z^2 - (172567 / 6000) * y^4 * z^2 +
        (8803 / 2520) * x^3 * z^3 +
        (58366387 / 3528000) * x^2 * y * z^3 +
        (109429927 / 15876000) * x * y^2 * z^3 - (1029331 / 302400) * y^3 * z^3 +
        (96953 / 12600) * x^2 * z^4 +
        (91482443 / 6350400) * x * y * z^4 +
        (19795369 / 2268000) * y^2 * z^4 +
        (48443 / 12600) * x * z^5 +
        (2391707 / 302400) * y * z^5 +
        (37 / 24) * z^6 +
        (1343 / 140) * x^5 +
        (24557 / 5880) * x^4 * y +
        (485399 / 75600) * x^3 * y^2 +
        (416237 / 21000) * x^2 * y^3 +
        (43973123 / 2268000) * x * y^4 - (90113 / 5400) * y^5 +
        (383909 / 58800) * x^4 * z +
        (9587758951 / 222264000) * x^3 * y * z +
        (4095697 / 216000) * x^2 * y^2 * z +
        (555196927 / 10584000) * x * y^3 * z - (18948401 / 756000) * y^4 * z +
        (369821 / 25200) * x^3 * z^2 +
        (4701868229 / 333396000) * x^2 * y * z^2 +
        (39287599 / 3969000) * x * y^2 * z^2 - (150216781 / 3969000) * y^3 * z^2 +
        (36475 / 6048) * x^2 * z^3 +
        (3551689 / 181440) * x * y * z^3 - (34806349 / 1764000) * y^2 * z^3 +
        (60101 / 13230) * x * z^4 - (19839887 / 2822400) * y * z^4 +
        (5869 / 8400) * z^5 +
        (1846717 / 226800) * x^4 - (23974597 / 3175200) * x^3 * y -
        (690342211 / 6804000) * x^2 * y^2 - (13130807 / 1134000) * x * y^3 -
        (185219 / 3600) * y^4 - (116834413 / 12700800) * x^3 * z -
        (175142299 / 12700800) * x^2 * y * z - (28534591 / 756000) * x * y^2 * z -
        (94226747 / 756000) * y^3 * z - (41734883 / 4233600) * x^2 * z^2 -
        (309684919 / 8467200) * x * y * z^2 - (4786829 / 37800) * y^2 * z^2 -
        (2765183 / 264600) * x * z^3 - (22368259 / 564480) * y * z^3 -
        (100873 / 11200) * z^4 +
        (1267871 / 136080) * x^3 +
        (128353 / 2800) * x^2 * y - (386551 / 18900) * x * y^2 - (971 / 240) * y^3 -
        (40316993 / 3175200) * x^2 * z + (43529203 / 2116800) * x * y * z -
        (1349153 / 50400) * y^2 * z - (58790633 / 4233600) * x * z^2 -
        (54765671 / 1058400) * y * z^2 - (421891 / 33600) * z^3 +
        (4559393 / 226800) * x^2 +
        (3208603 / 151200) * x * y +
        (593 / 240) * y^2 +
        (104519 / 12600) * x * z - (26287 / 16800) * y * z - (88687 / 13440) * z^2 +
        (257 / 120) * x +
        (2153 / 80) * y +
        (5217 / 560) * z +
        171 / 16,
    ]

    System(horner.(exprs))
end

function six_revolute()
    @var z[1:6, 1:3]

    p = [1, 1, 0]
    α = [
        0.9606655201575953,
        -0.12634585123792536,
        -0.821994270849398,
        -0.4612529209681106,
        0.1914597581898848,
    ]
    a = [
        -0.6508165773155992,
        -0.6252070694344862,
        1.2753203087262615,
        -0.30371050526688304,
        0.45398822665915634,
        -0.4042256983795828,
        -1.1485430048515646,
        -0.22145594990414993,
        0.9590486265536351,
    ]
    ## define the system of polynomials
    f = [z[i, :] ⋅ z[i, :] for i = 2:5]
    g = [z[i, :] ⋅ z[i+1, :] for i = 1:5]
    h = sum(a[i] .* (z[i, :] × z[i+1, :]) for i = 1:3) + sum(a[i+4] .* z[i, :] for i = 2:5)
    F′ = [f .- 1; g .- cos.(α); h .- p]
    ## assign values to z₁ and z₆
    System(
        [subs(f, z[1, :] => [1, 0, 0], z[6, :] => [1, 0, 0]) for f in F′],
        variable_groups = [[z[2, :]; z[4, :]], [z[3, :]; z[5, :]]],
    )
end

function fano_quintic()
    q₀ = Complex{Float64}[
        0.3836468979754316+0.48294854343585913im,
        -0.3997845927222694-0.1872020963239582im,
        0.2877641371295415+0.13959986611131697im,
        1.0246085736811654-0.6706882552971819im,
        0.35467278956822795-0.6821462156768627im,
        -0.5186542729891018-0.6472352083267523im,
        -0.5803213495160932-1.1604897218476604im,
        0.2216615210358535+0.6447244163071337im,
        0.3650725307633961-1.423170740516543im,
        -0.07782388627393373+0.37375061147457606im,
        -0.267970179805653+0.5178907448772292im,
        -1.3039508977772771+0.9244564262790848im,
        0.05070983247723053-0.49222470085574327im,
        0.01585315985561874-1.5964641275845275im,
        -0.6702797151556827+0.899455810440721im,
        0.16123533462036213-0.14249720030175353im,
        -0.05801647381933075+0.02408769940754098im,
        -0.6428310495796878+0.3768196415291366im,
        0.4236920704992038-0.7662716901564756im,
        0.8308019425144062-0.33353922910187783im,
        -1.7206825953784448-0.023068433837692474im,
        -1.6674450156673142+0.1351816982748416im,
        1.016891843231412+1.5786636727636165im,
        0.3151651494431033+0.355704156731951im,
        0.7515693388748262+0.6299666874414362im,
        -0.36820201739484165-1.4285021667529485im,
        -1.2355050806442094-0.28609402812674734im,
        0.554735831573422-1.0149295902971325im,
        0.47248893311957996-0.23261176616023507im,
        -2.013652724617146-0.11466869282430801im,
        -1.5396892983622734-0.2949664781961386im,
        -0.4184209781159971+0.5449642760480936im,
        1.2144346526597238-0.16868272729235662im,
        0.9173137286532203+0.20477977881028425im,
        -0.4450853800437796+1.3593858094881657im,
        -0.5720658730688182+1.0763344735763516im,
        -1.9171254167925773-0.2747011694692946im,
        1.1921942553565317+0.44575106539716924im,
        1.0350508008515482+0.10009830817916106im,
        0.3901281635062643+0.08000123439730099im,
        -1.0714116835498375+1.157205798671797im,
        -0.6637087604243677+0.6461404584356966im,
        0.2214679807586311+0.22437921498443017im,
        0.8458782102538259-0.32826055277382443im,
        0.1025946064372987+1.037567452186008im,
        0.2377278791681354+0.6740608327834566im,
        0.49532004094565674+1.9952174010506525im,
        0.4568307118448063+0.9631852448527662im,
        0.4788652739248191-0.40455695253590435im,
        0.09221652636081985+0.7625458623610796im,
        -0.12425240651640111+1.6816982819349031im,
        0.7893904624781745+0.8629775422589433im,
        0.5721436751408631-0.4044177992453996im,
        -0.027875570104326758+0.012968741834617572im,
        -0.2858145461151917-0.04405199091539911im,
        -0.24069919475270554-0.060754487173562834im,
        0.7037314454480748-0.07741277292123394im,
        0.0372636983547525+0.6201252721270679im,
        0.7923164779358874-0.041045066455393846im,
        -0.7193308083445518+0.4030785967993166im,
        0.6878829645928711+0.7025043076112989im,
        0.8169863844877365-0.5911126929082672im,
        0.04782856084775355-0.8054508299944776im,
        -0.7027374989526586-1.0359616901446203im,
        0.6056428944077695+0.49930791024376064im,
        0.6412321091406482+1.314885022146572im,
        -0.331319266798614-0.10429521678269742im,
        -1.0913855443020988-0.5722558126816987im,
        -0.31090086608681744-0.12347310482357825im,
        0.0689555157196041-0.5398051922073547im,
        1.0748702515974689+1.5362064120254992im,
        -0.05290485395797417-0.3251779508261342im,
        -0.6463799980658426+0.44096989535814457im,
        -0.2743821036302145-0.7255061059609375im,
        0.23153545237088385-0.40167584002518003im,
        -0.418576962372166-0.3186349937115043im,
        -0.9477795397466203-0.8156397749171547im,
        0.047917328296541044+1.3174326463487218im,
        -1.2567854006147323-0.8193981780145082im,
        -0.3787559751573031-0.4669776616037851im,
        0.09351143095398255+0.4679778372533676im,
        0.02692492421670912+0.44818461911684im,
        -0.5027816524030012-0.3050313840292516im,
        0.8432030208496721+1.165070707855424im,
        0.3238182623736368+0.6538911888669904im,
        -0.40612081369460046+0.737432066207105im,
        -1.2268726064651367+0.4612772091237368im,
        -0.23840730064416424-0.6248741829761927im,
        0.5301718981960735+0.23444283568420737im,
        -0.07798620446200054+1.4654832951672472im,
        1.2064511977634496+0.5753299611026544im,
        -1.1107367739855873+0.4774660709904277im,
        1.6584330667420542+0.7443869745535754im,
        -0.05515230499908627+1.3317792685092573im,
        -0.06795807047004142-0.20267847847864062im,
        -0.4368275386474442-0.6988998583364666im,
        -0.17797910156645358-0.0032778322183904625im,
        2.2249253540741516-0.44105008069846247im,
        -0.8198149986755516-0.2481966019221811im,
        -0.011334319833761836+0.15159208016548192im,
        0.2940860308353497-0.22656517901536938im,
        -1.2121569288336493-0.4413596579347592im,
        1.5710495613404776+0.6487960514719261im,
        -0.6704069855852188-0.40288962747886303im,
        0.06692745796474238-0.0347684716960473im,
        -0.4729225318917847+0.4661833789708642im,
        0.8299215271862779-0.21461942234118758im,
        0.03397633025259013+0.012175348959508479im,
        -0.5977599028309415-1.1275243723226982im,
        0.048101857348091144-0.043643341785888815im,
        -0.4177427390228452+0.995034345527883im,
        -0.3188503572583698-0.16007925623999167im,
        -0.35482792595534385+0.06235995941860926im,
        0.14546523772036224+0.7779115806461815im,
        0.022246178997971426-0.03345409806041387im,
        -0.9103291782318836+0.28858060999716423im,
        0.7869434800704376-0.1974343252882416im,
        -0.2777362225119344-0.6552357376785277im,
        -0.9419716203488216-0.39726375447789025im,
        -0.6770033931124602+0.11295943019210838im,
        0.6894080970748474-1.9814830784431494im,
        0.7110844987437178-0.23660236950768132im,
        -0.41761833997405756-1.506846428650273im,
        -0.4513212474271783-1.0674194728486397im,
        -0.40897512850909873+0.054576738818524116im,
    ]
    gamma = 0.63 - 0.52im

    n = 4
    @var x[1:n]

    F, q = dense_poly(x, 5; coeff_name = :q)
    F = F(q[end] => 1)
    q = q[1:end-1]
    @var a[1:n-1] b[1:n-1] t
    L = [a; 1] .* t + [b; 0]
    FcapL = last(ModelKit.exponents_coefficients(subs(F, x => L), [t]))
    sys = System(ModelKit.horner.(FcapL), [a; b], q)
    sys, q₀, gamma
end

function steiner()
    @var x[1:2] a[1:5] c[1:6] y[1:2, 1:5]

    #tangential conics
    f = a[1] * x[1]^2 + a[2] * x[1] * x[2] + a[3] * x[2]^2 + a[4] * x[1] + a[5] * x[2] + 1
    ∇ = differentiate(f, x)
    #5 conics
    g =
        c[1] * x[1]^2 +
        c[2] * x[1] * x[2] +
        c[3] * x[2]^2 +
        c[4] * x[1] +
        c[5] * x[2] +
        c[6]
    ∇_2 = differentiate(g, x)
    #the general system
    #f_a_0 is tangent to g_b₀ at x₀
    function Incidence(f, a₀, g, b₀, x₀)
        fᵢ = f(x => x₀, a => a₀)
        ∇ᵢ = ∇(x => x₀, a => a₀)
        Cᵢ = g(x => x₀, c => b₀)
        ∇_Cᵢ = ∇_2(x => x₀, c => b₀)
        [fᵢ; Cᵢ; det([∇ᵢ ∇_Cᵢ])]
    end
    @var v[1:6, 1:5]
    F = vcat(map(i -> Incidence(f, a, g, v[:, i], y[:, i]), 1:5)...)
    System(F, [a; vec(y)], vec(v))
end

function moments3()
    @var a[1:3] x[1:3] s[1:3] m[0:8]

    f0 = a[1] + a[2] + a[3]
    f1 = a[1] * x[1] + a[2] * x[2] + a[3] * x[3]
    f2 = a[1] * (x[1]^2 + s[1]) + a[2] * (x[2]^2 + s[2]) + a[3] * (x[3]^2 + s[3])
    f3 =
        a[1] * (x[1]^3 + 3 * s[1] * x[1]) +
        a[2] * (x[2]^3 + 3 * s[2] * x[2]) +
        a[3] * (x[3]^3 + 3 * s[3] * x[3])
    f4 =
        a[1] * (x[1]^4 + 6 * s[1] * x[1]^2 + 3 * s[1]^2) +
        a[2] * (x[2]^4 + 6 * s[2] * x[2]^2 + 3 * s[2]^2) +
        a[3] * (x[3]^4 + 6 * s[3] * x[3]^2 + 3 * s[3]^2)
    f5 =
        a[1] * (x[1]^5 + 10 * s[1] * x[1]^3 + 15 * x[1] * s[1]^2) +
        a[2] * (x[2]^5 + 10 * s[2] * x[2]^3 + 15 * x[2] * s[2]^2) +
        a[3] * (x[3]^5 + 10 * s[3] * x[3]^3 + 15 * x[3] * s[3]^2)
    f6 =
        a[1] * (x[1]^6 + 15 * s[1] * x[1]^4 + 45 * x[1]^2 * s[1]^2 + 15 * s[1]^3) +
        a[2] * (x[2]^6 + 15 * s[2] * x[2]^4 + 45 * x[2]^2 * s[2]^2 + 15 * s[2]^3) +
        a[3] * (x[3]^6 + 15 * s[3] * x[3]^4 + 45 * x[3]^2 * s[3]^2 + 15 * s[3]^3)
    f7 =
        a[1] * (x[1]^7 + 21 * s[1] * x[1]^5 + 105 * x[1]^3 * s[1]^2 + 105 * x[1] * s[1]^3) +
        a[2] * (x[2]^7 + 21 * s[2] * x[2]^5 + 105 * x[2]^3 * s[2]^2 + 105 * x[2] * s[2]^3) +
        a[3] * (x[3]^7 + 21 * s[3] * x[3]^5 + 105 * x[3]^3 * s[3]^2 + 105 * x[3] * s[3]^3)
    f8 =
        a[1] * (
            x[1]^8 +
            28 * s[1] * x[1]^6 +
            210 * x[1]^4 * s[1]^2 +
            420 * x[1]^2 * s[1]^3 +
            105 * s[1]^4
        ) +
        a[2] * (
            x[2]^8 +
            28 * s[2] * x[2]^6 +
            210 * x[2]^4 * s[2]^2 +
            420 * x[2]^2 * s[2]^3 +
            105 * s[2]^4
        ) +
        a[3] * (
            x[3]^8 +
            28 * s[3] * x[3]^6 +
            210 * x[3]^4 * s[3]^2 +
            420 * x[3]^2 * s[3]^3 +
            105 * s[3]^4
        )
    System([f0, f1, f2, f3, f4, f5, f6, f7, f8] - m, parameters = m)
end


text:
## Exercise 4-7
## Make a more general version of circle called arc that takes an additional parameter angle, which determines what fraction of a circle to draw. angle is in units of degrees, so when angle = 360, arc should draw a complete circle.
using ThinkJulia

println("Ans: ")

function arc(turtle::Turtle, radius::Int, angle)
    len = 2 * π * radius * angle / 360  # finding the length of the arc 2πr * a/360
    n = Int(trunc(len) / 3) + 1 # number of times the turtle will rotate and draw

    len /= n
    angle /= n

    for i in 1:n
        forward(turtle, len)
        turn(turtle, -angle)
    end
end

@svg begin
    arc(Turtle(), 200, 270)
end

println("End.")


text:
# FASTA Reader
# ============

"""
    FASTAReader(input::IO; index=nothing)
    FASTAReader{S}(input::IO; index=nothing)

Create a data reader of the FASTA file format.

When type parameter `S` is specified, the reader reads sequences in that type;
otherwise the reader tries to infer the sequence type based on the frequencies
of characters from the input.

# Arguments
* `input`: data source
* `index=nothing`: filepath to a random access index (currently *fai* is supported)
"""
type FASTAReader{S<:Sequence} <: Bio.IO.AbstractReader
    state::Ragel.State
    seqbuf::BufferedOutputStream{BufferedStreams.EmptyStream}
    index::Nullable{FASTAIndex}

    function FASTAReader(input::BufferedInputStream, index)
        return new(Ragel.State(fastaparser_start, input),
                   BufferedOutputStream(), index)
    end
end

function Bio.IO.stream(reader::FASTAReader)
    return reader.state.stream
end

function Base.eltype{S}(::Type{FASTAReader{S}})
    return FASTASeqRecord{S}
end

function FASTAReader(input::IO; index=nothing)
    # alphabet type will be inferred from data
    return FASTAReader{BioSequence}(input, index=index)
end

function (::Type{FASTAReader{S}}){S<:Sequence}(input::IO; index=nothing)
    if isa(index, AbstractString)
        index = FASTAIndex(index)
    else
        if index != nothing
            error("unrecognizable index argument")
        end
    end
    return FASTAReader{S}(BufferedInputStream(input), index)
end

function Base.getindex(reader::FASTAReader, name::AbstractString)
    if isnull(reader.index)
        error("no index")
    end
    seekrecord(reader.state.stream, get(reader.index), name)
    reader.state.cs = fastaparser_start
    return read(reader)
end

# Predict sequence type based on character frequencies in `seq[start:stop]`.
function predict(seq::Vector{UInt8}, start, stop)
    # count characters
    a = c = g = t = u = n = alpha = 0
    for i in start:stop
        @inbounds x = seq[i]
        if x == 0x41 || x == 0x61
            a += 1
        elseif x == 0x43 || x == 0x63
            c += 1
        elseif x == 0x47 || x == 0x67
            g += 1
        elseif x == 0x54 || x == 0x74
            t += 1
        elseif x == 0x55 || x == 0x75
            u += 1
        elseif x == 0x4e || x == 0x6e
            n += 1
        end
        if 0x41 ≤ x ≤ 0x5a || 0x61 ≤ x ≤ 0x7a
            alpha += 1
            if alpha ≥ 300 && t + u > 0 && a + c + g + t + u + n == alpha
                # pretty sure that the sequence is either DNA or RNA
                break
            end
        end
    end

    # the threshold (= 0.95) is somewhat arbitrary
    if (a + c + g + t + u + n) / alpha > 0.95
        if t ≥ u
            return DNASequence
        else
            return RNASequence
        end
    else
        return AminoAcidSequence
    end
end


text:
import Roots.fzero

# test robustness of derivative free algorithms

## tests:http://people.sc.fsu.edu/~jburkardt/cpp_src/test_zero/test_zero.html
function newton_baffler(x)
    if ( x - 0.0 ) < -0.25
        0.75 * ( x - 0 ) - 0.3125
    elseif  ( x - 0 ) < 0.25
        2.0 * ( x - 0 )
    else
        0.75 * ( x - 0 ) + 0.3125
    end
end

mutable struct CallableFunction
f
n
end
(F::CallableFunction)(x) = (F.n+=1; F.f(x))

pathological = [
                (x ->sin( x ) - x / 2, .1),
                (x-> 2 * x - exp( - x ), 1),
                (x -> x * exp( - x ), 1),
                (x -> exp( x ) - 1 / ( 10 * x )^2, 1),
                (x -> ( x + 3 ) * ( x - 1 )^2, 1),

                (x -> exp( x ) - 2 - 1 / ( 10 * x )^2 + 2 / ( 100 * x )^3, 1),
                (x -> cbrt(x), 1),
                (x -> cos( x ) - x, 1),
                (x-> newton_baffler(x), 8),
                (x -> 20.0 * x / ( 100.0 * x^2 + 1.0), 0.095),

                (x ->  ( 4.0 + x^2) * ( 2.0 + x ) * ( 2.0 - x )  / ( 16.0 * x * x * x * x + 0.00001 ), 1),
                (x -> (x == 1.0) ? float(0) : sign(x-1.0) * exp(log(1e4) + log(abs(x - 1.0)) - 1.0/(x-1.0)^2), 1),
                (x -> 0.00000000001 * (x - 100.0), 1),
                (x -> 1.0 / ( ( x - 0.3 ) * ( x - 0.3 ) + 0.01 ) + 1.0 / ( ( x - 0.9 ) * ( x - 0.9 ) + 0.04 ) + 2.0 * x - 5.2, -1),
                (x -> ( 1 - 6x^2) * cbrt(x) * exp(-x^2) / (3*x), -0.25),

                (x -> ( pi * ( x - 5.0 ) / 180.0 ) - 0.8 * sin( pi * x / 180.0 ), 1),
                (x -> x^3 - 2*x - 5, 2),
                (x -> 1e6 * (x^7 -7x^6 +21x^5 -35x^4 +35x^3-21x^2+7x-1),  0.990),
                (x -> cos(100*x)-4*erf(30*x-10), 0.0)
                ]

function howgood(f1, xstar)
    if isnan(xstar)
        "x"
    elseif iszero(f1(xstar))       # perfect
        "."
    elseif f1(prevfloat(xstar)) * f1(nextfloat(xstar)) < 0 # perfect
        "."
    elseif abs(f1(xstar)) < 8 * eps(xstar) # nearly perfect
        "."
    elseif abs(f1(xstar)) <= 10eps(xstar)^(1/3)
        "~"
    elseif abs(xstar) > 5   # ran away yet converged
        "!"
    else                    # not close?
        "*"
    end
end

## make a graphic comparing values
function run_robustness_test(bigfloats=false)
    orders = [Order0(),Order1(),Order2(),Order5(),Order8(),Order16()]
    N = length(orders) + 1
    k,n = length(pathological), 50
    m = ["" for i=1:(N + 2)*k, j=1:n]
    ## how to mark values. If error we also use "x"

    marks = ["HP", "0 ", "1 ", "2 ", "5 ", "8 ", "16"]
    io = IOBuffer()

    for i in 1:k

        println("")
        println("-- $i ---")
        println("")

        f1,x = pathological[i]
        if bigfloats
            xs = -big(5):10/50:5
        else
            xs = linspace(-5, 5, n)
        end

        calls = Int[0]
        hits = 0
        j = 1
        print(io, marks[j], " ")
        for k in 1:n
            F = CallableFunction(f1, 0)
            xstar = Roots.dfree(F, xs[k])
            push!(calls, F.n)
            val = howgood(f1, xstar)
            val in [".", "~"] && (hits += 1)
            print(io, val)
        end
        print(String(take!(io)))
        print(" h=$hits n = $(sum(calls)) max=$(maximum(calls))")
        println("")

        for j in 2:N
            calls = Int[0]
            hits = 0
            print(io, marks[j], " ")
            for k in 1:n
                M = orders[j-1]
                F = CallableFunction(f1, 0)
                state = Roots.init_state(M, F, xs[k])
                options = Roots.init_options(M, state)
                tracks = Roots.NullTracks()
                xstar = find_zero(M, F, options, state, tracks)
                push!(calls, F.n)
                val = howgood(f1, xstar)
                val in [".", "~"] && (hits += 1)
                print(io, val)
            end
            print(String(take!(io)))
            print(" h=$hits n = $(sum(calls)) max=$(maximum(calls))")
            println("")
        end
    end
end


text:
using NPZ
using PyCall
using Healpix

np = pyimport("numpy")

spath = "/group/cmb/litebird/usr/ytakase/beam_convolution/devided_maps/outputs/mapdata"

nside = 512
npix = nside2npix(nside)

RunRange = 1536
outmap = zeros(2, npix)
for i in 1:RunRange
    outmap .+= npzread(spath*"/idx=$i.npz")["convolved_map"]
end

np.savez_compressed("/group/cmb/litebird/usr/ytakase/beam_convolution/devided_maps/result/hiroaki_conv", outmap = outmap)


text:
"""
    RandomWalk{S,T,U<:Sampleable} <: AbstractMarkovChain

##### Fields
- `statespace::S` : the state-space of the random walk (e.g. an integer lattice or a graph).
- `initialstate::T` : the state at which the random walk begins (at time-step 1).
- `stepsampler::U` : sampler that allows i.i.d. steps to be drawn.
"""
struct RandomWalk{S,T,U<:Sampleable} <: AbstractMarkovChain
    statespace::S
    initialstate::T
    stepsampler::U
end

"""
    Base.show(io::IO, randomwalk::RandomWalk) -> nothing

Show a `RandomWalk` in a human-friendly manner.
"""
function Base.show(io::IO, randomwalk::RandomWalk)
    println(io, typeof(randomwalk))
    println(io, " State Space:")
    println(io, "  ", randomwalk.statespace)
    println(io, " Initial state:")
    println(io, "  ", randomwalk.initialstate)
    println(io, " Sampler that allows i.i.d. steps to be drawn:")
    println(io, "  ", randomwalk.stepsampler)
    return nothing
end

"""
    Base.eltype(randomwalk::RandomWalk) -> DataType

Return the type associated with `randomwalk`'s domain.
"""
function Base.eltype(randomwalk::RandomWalk)
    return typeof(randomwalk.initialstate)
end

"""
    RandomWalk(latticedim::Integer, initialstate=zeros(Int, latticedim),
               stepsampler::Sampleable = ElementVectorSampler(latticedim) -> RandomWalk

Construct and return a `RandomWalk` on the integer lattice with dimension `latticedim`.

##### Complexity
- Time complexity: `O(1)`.
"""
function RandomWalk(latticedim::Integer=1, initialstate=zeros(Int,latticedim),
                    stepsampler::Sampleable = ElementVectorSampler(latticedim))
    latticedim < 1 && throw(ArgumentError("latticedim must be positive."))
    length(initialstate) != latticedim &&
        throw(DimensionMismatch("initial state's length should equal latticedim"))
    if latticedim == 1
        return RandomWalk(Int, initialstate[1], stepsampler)
    else
        return RandomWalk(repeat([Int], latticedim), initialstate, stepsampler)
    end
end

"""
    RandomWalk(graph::AbstractGraph, initialvertex=rand(1:nv(graph));
               checkargs::Bool = true) -> RandomWalk

Construct and return a `RandomWalk` on a graph.

##### Complexity
- Time complexity: `O(n)` when `checkargs == true`; otherwise, `O(1)`.
"""
function RandomWalk(graph::AbstractGraph, initialvertex=rand(1:nv(graph));
                    checkargs::Bool = true)
    if checkargs && !(initialvertex in vertices(graph))
        throw(ArgumentError("initialvertex must be a vertex in graph)."))
    else
        return RandomWalk(graph, initialvertex, NeighborSampler(graph, initialvertex))
    end
end

"""
    randomwalkmc(graph::AbstractGraph, initialvertex=rand(1:nv(graph));
                 checkargs::Bool = true) -> MarkovChain

Construct a `MarkovChain` representation of a random walk by reducing to the `MarkovChain`
API.

Observe that this reduction only works for random walks on finite graphs, for `MarkovChain`
employs a finite transition matrix. If one needs a random walk on an infinite graph, then
use the `RandomWalk` constructor designed to dispatch on `AbstractGraph`s.

##### Complexity
- Time complexity: `O(n)` when `checkargs == true`; otherwise, `O(1)`.
"""
function randomwalkmc(graph::AbstractGraph, initialvertex=rand(1:nv(graph));
                      checkargs::Bool = true)
    if checkargs && !(initialvertex in vertices(graph))
        throw(ArgumentError("initialvertex must be a vertex in graph)."))
    else
        initialdist = zeros(nv(graph))
        initialdist[initialvertex] = 1
        # broadcast division
        return MarkovChain(Diagonal(1 ./ degree(graph))*adjacency_matrix(graph), initialdist)
    end
end

max_stars_repo_name: LoganDGraham/StochasticProcesses

max_issues_repo_name: LoganDGraham/StochasticProcesses

max_forks_repo_name: LoganDGraham/StochasticProcesses


text:
# write your own tests here
include("singlethreaded.jl")


text:
"""
    makeAxisArray(I::AbstractArray{T,5}, acqData::AcquisitionData) where T

creates an axies array with properly characterized axis from the image `I`.
For this, it uses the information in `acqData`.
The axies respectively describe the following coordinates:
x, y, z, acqData.numEchoes, acqData.numCoils
"""
function ImageUtils.makeAxisArray(I::AbstractArray{T,5}, acqData::AcquisitionData) where T

  offset = [0.0, 0.0, 0.0]*Unitful.mm
  spacing = pixelspacing(acqData)

  im = AxisArray(I,
		   Axis{:x}(range(offset[1], step=spacing[1], length=size(I,1))),
		   Axis{:y}(range(offset[2], step=spacing[2], length=size(I,2))),
		   Axis{:z}(range(offset[3], step=spacing[3], length=size(I,3))),
		   Axis{:echos}(1:size(I,4)),
       Axis{:coils}(1:size(I,5)))

  #imMeta = ImageMeta(im, Dict{String,Any}())
  #return imMeta
  return im
end

function ImageUtils.makeAxisArray(I::AbstractArray{T,5}, spacing::Vector{Float64}) where T

  offset = [0.0, 0.0, 0.0]*Unitful.mm

  sp = uconvert.(Unitful.mm, spacing*Unitful.m)

  im = AxisArray(I,
		   Axis{:x}(range(offset[1], step=sp[1], length=size(I,1))),
		   Axis{:y}(range(offset[2], step=sp[2], length=size(I,2))),
		   Axis{:z}(range(offset[3], step=sp[3], length=size(I,3))),
		   Axis{:echos}(1:size(I,4)),
       Axis{:coils}(1:size(I,5)))

  return im
end


text:
# # Code derived from metabolicEP (https://github.com/anna-pa-m/Metabolic-EP)

# function epconverge!(epfields::EPFields, epmat::M, epalg::EPAlg, eponesweep!::T) where {T<:Function,M<:AbstractEPMat}

#     @extract epalg : maxiter verbose epsconv alpha beta_vec

#     returnstatus = :unconverged
#     iter = 0

#     # sweep ep till maxiter is reached or max(errav, errvar) < epsconv
#     # alpha = epalg.alpha epsconv maxiter
#     prog = ProgressThresh(epsconv)# ; desc = "MaxEnt EP... ")
#     max_beta = findmax(beta_vec)
#     while iter < maxiter
#         iter += 1
#         # eponesweep! will be eponesweepT0! or eponesweep depending on alpha
#         (errav, errvar, errμ, errs) = eponesweep!(epfields, epalg, epmat)

#         max_err = max(errav, errvar)
#         if max_err < epsconv
#             returnstatus = :converged
#             break
#         end

#         verbose && update!(prog, max_err; showvalues = [
#             (:iter, iter),
#             (:maxiter, maxiter),
#             (:alpha, alpha),
#             (:max_beta, max_beta)
#         ])
#     end

#     verbose && (finish!(prog); flush(stderr))

#     return returnstatus, iter
# end


text:
using Peaks
using Test, OffsetArrays

include("minmax.jl")
include("peakprom.jl")
include("peakwidth.jl")

max_stars_repo_name: halleysfifthinc/Peaks

max_issues_repo_name: halleysfifthinc/Peaks

max_forks_repo_name: halleysfifthinc/Peaks


text:
module UltraDark

using Base.Threads: @threads
using LinearAlgebra
using AbstractFFTs: fftfreq, rfftfreq
using PencilFFTs, MPI

using FFTW

export simulate
export Grids, PencilGrids
export Config, SimulationConfig, constant_scale_factor, TimeStepOptions
export OutputConfig

include("grids.jl")
include("pencil_grids.jl")
include("phase_diff.jl")
include("time_step.jl")
include("output.jl")
include("config.jl")

import .Output: OutputConfig
import .Output: output_summary_row, output_summary_header, output_grids
import .Config: SimulationConfig, constant_scale_factor, TimeStepOptions

function outer_step!(Δt, grids, constants; a=1.0)
    psi_whole_step!(Δt, grids, constants)
end

"""
    psi_half_step

Deprecated.
Use `outer_step` with Δt/2 instead.
"""
function psi_half_step!(Δt, grids, constants)
    @inbounds @threads for i in eachindex(grids.ψx)
        grids.ψx[i] *= exp(- im * Δt / 2 * grids.Φx[i])
    end
end

"""
    psi_half_step

Deprecated.
Use `outer_step` instead.
"""
function psi_whole_step!(Δt, grids, constants)
    @inbounds @threads for i in eachindex(grids.ψx)
        grids.ψx[i] *= exp(- im * Δt / 1 * grids.Φx[i])
    end
end

function inner_step!(Δt, grids, constants; a=1.0)
    phi_whole_step!(Δt, grids, constants; a=1.0)
end

"""
    phi_whole_step

Deprecated.
Use `inner_step` instead.
"""
function phi_whole_step!(Δt, grids, constants; a=1.0)
    # TODO: not all part of Φ update

    mul!(grids.ψk, grids.fft_plan, grids.ψx)
    @inbounds @threads for i in eachindex(grids.ψk)
        grids.ψk[i] *= exp(-im * Δt/2 * grids.k[i]^2 / a^2)
    end
    ldiv!(grids.ψx, grids.fft_plan, grids.ψk)

end

function update_gravitational_potential!(grids, constants; a=1.0)
    @inbounds @threads for i in eachindex(grids.ρx)
        grids.ρx[i] = abs2(grids.ψx[i])
    end

    mul!(grids.Φk, grids.rfft_plan, grids.ρx)
    @inbounds @threads for i in eachindex(grids.Φk)
        grids.Φk[i] *= -4 * π / (a * grids.rk[i]^2)
    end
    grids.Φk[grids.rk_vanish_indices] .= 0
    ldiv!(grids.Φx, grids.rfft_plan, grids.Φk)
end

"""
    auxiliary_step!(Δt, grids, t, constants)

Do an auxiliary inner step.
By default this does nothing, but can be overridden in multiple dispatch.
"""
function auxiliary_step!(Δt, grids, t, constants)
end

"""
    add_external_potential!(t, grids, constants)

Add a gravitational potential to the grid.
By default this does nothing, but can be overridden in multiple dispatch.
"""
function add_external_potential!(t, grids, constants)
end

"""
Take `n` steps with time step `Δt`

# Examples

```jldoctest
julia> using UltraDark: take_steps!, Grids, OutputConfig, Config

julia> take_steps!(Grids(1.0, 16), 0, 0.5, 10, OutputConfig(mktempdir(), []), Config.constant_scale_factor, nothing)
5.0

```
"""
function take_steps!(grids, t_start, Δt, n, output_config, a, constants)

    t::Float64 = t_start

    half_step = true

    for step in 1:n
        if half_step
            outer_step!(Δt/2, grids, constants; a=a(t))
            t += Δt / 2
            half_step = false
        else
            outer_step!(Δt, grids, constants; a=a(t))
            t += Δt
        end

        inner_step!(Δt, grids, constants; a=a(t))
        update_gravitational_potential!(grids, constants; a=a(t))
        add_external_potential!(t, grids, constants)
        auxiliary_step!(Δt, grids, t, constants)

        output_summary_row(grids, output_config, t, a(t), Δt)
    end

    outer_step!(Δt/2, grids, constants)
    t += Δt / 2

    output_summary_row(grids, output_config, t, a(t), Δt)

    t
end

"""
Evolve `grids` forward from `t_start` to `t_end`
"""
function evolve_to!(t_start, t_end, grids, output_config, sim_config; constants=nothing)

    @assert t_start < t_end

    t = t_start

    while (t < t_end) && ~(t ≈ t_end)

        if ~good_phase_diff(grids, sim_config)
            output_grids(grids, output_config, -1)
            throw("Phase gradient is too large to continue")
        end

        Δt, n_steps = actual_time_step(
            max_time_step(grids, sim_config.a(t), sim_config.time_step_options),
            t_end - t,
            sim_config.time_step_options,
        )

        t = take_steps!(grids, t, Δt, n_steps, output_config, sim_config.a, constants)
    end

    t
end

function simulate(grids, sim_config, output_config::OutputConfig; constants=nothing)

    # Setup output
    mkpath(output_config.directory)
    output_summary_header(output_config)

    t_begin = output_config.output_times[1]

    # Initialise vars other than ψx
    # This is required so the initial time step can be calculated
    t_initial = output_config.output_times[1]
    update_gravitational_potential!(grids, constants; a=sim_config.a(t_initial))
    add_external_potential!(t_initial, grids, constants)

    # Output initial conditions
    output_grids(grids, output_config, 1)
    if ~good_phase_diff(grids, sim_config)
        throw("Initial phase gradient is already in untrusted regime")
    end

    for (index, t_end) in enumerate(output_config.output_times[2:end])
        t_begin = evolve_to!(
            t_begin,
            t_end,
            grids,
            output_config,
            sim_config;
            constants=constants,
        )
        @info "Reached time $t_begin"
        output_grids(grids, output_config, index + 1)
    end

    if ~good_phase_diff(grids, sim_config)
        throw("Phase gradient is too large to end")
    end

end

end # module


text:
export lineargradient, back!, GradAlg, value, grad

"Gradient Algorithm"
abstract type GradAlg end

"`lineargradient(rv), ω::AbstractΩ, ::GradAlg)` Returns as vector gradient of ω components"
function lineargradient end

"`back!(rv, ω::AbstractΩ, ::GradAlg)` update ω s.t ω[id].grad is gradient"
function back! end

"`grad(rv, ω, ::GradAlg)` returns ω where gradients in components"
function grad end

# "Value associated with dual/tracked/etc number/array"
# function value end

max_stars_repo_name: zenna/expect

max_issues_repo_name: zenna/expect

max_forks_repo_name: zenna/expect


text:
using LinearAlgebra, Statistics
# commandline arguments.  In reality, use ArgParse equivalents
num = parse(Int64,ARGS[2])
ϵ = randn(num)
println("mean = $(mean(ϵ)), std = $(std(ϵ)) of $num values")


text:
# we store the actual julia values here
# the `value` field of `PyJuliaValueObject` indexes into here
const PYJLVALUES = []
# unused indices in PYJLVALUES
const PYJLFREEVALUES = Int[]

function _pyjl_new(t::PyPtr, ::PyPtr, ::PyPtr)
    o = ccall(UnsafePtr{PyTypeObject}(t).alloc[!], PyPtr, (PyPtr, Py_ssize_t), t, 0)
    o == PyNULL && return PyNULL
    UnsafePtr{PyJuliaValueObject}(o).weaklist[] = PyNULL
    UnsafePtr{PyJuliaValueObject}(o).value[] = 0
    return o
end

function _pyjl_dealloc(o::PyPtr)
    idx = UnsafePtr{PyJuliaValueObject}(o).value[]
    if idx != 0
        PYJLVALUES[idx] = nothing
        push!(PYJLFREEVALUES, idx)
    end
    UnsafePtr{PyJuliaValueObject}(o).weaklist[!] == PyNULL || PyObject_ClearWeakRefs(o)
    ccall(UnsafePtr{PyTypeObject}(Py_Type(o)).free[!], Cvoid, (PyPtr,), o)
    nothing
end

const PYJLMETHODS = Vector{Any}()

function PyJulia_MethodNum(f)
    @nospecialize f
    push!(PYJLMETHODS, f)
    return length(PYJLMETHODS)
end

function _pyjl_isnull(o::PyPtr, ::PyPtr)
    ans = PyJuliaValue_IsNull(o) ? POINTERS._Py_TrueStruct : POINTERS._Py_FalseStruct
    Py_IncRef(ans)
    ans
end

function _pyjl_callmethod(o::PyPtr, args::PyPtr)
    nargs = PyTuple_Size(args)
    @assert nargs > 0
    num = PyLong_AsLongLong(PyTuple_GetItem(args, 0))
    num == -1 && return PyNULL
    f = PYJLMETHODS[num]
    # this form gets defined in jlwrap/base.jl
    return _pyjl_callmethod(f, o, args, nargs)::PyPtr
end

const PYJLBUFCACHE = Dict{Ptr{Cvoid},Any}()

@kwdef struct PyBufferInfo{N}
    # data
    ptr::Ptr{Cvoid}
    readonly::Bool
    # items
    itemsize::Int
    format::String
    # layout
    shape::NTuple{N,Int}
    strides::NTuple{N,Int}
    suboffsets::NTuple{N,Int} = ntuple(i -> -1, N)
end

_pyjl_get_buffer_impl(obj::PyPtr, buf::Ptr{Py_buffer}, flags::Cint, x, f) = _pyjl_get_buffer_impl(obj, buf, flags, f(x)::PyBufferInfo)

function _pyjl_get_buffer_impl(obj::PyPtr, buf::Ptr{Py_buffer}, flags::Cint, info::PyBufferInfo{N}) where {N}
    b = UnsafePtr(buf)
    c = []

    # not influenced by flags: obj, buf, len, itemsize, ndim
    b.obj[] = C_NULL
    b.buf[] = info.ptr
    b.itemsize[] = info.itemsize
    b.len[] = info.itemsize * prod(info.shape)
    b.ndim[] = N

    # readonly
    if !info.readonly
        b.readonly[] = 0
    elseif Utils.isflagset(flags, PyBUF_WRITABLE)
        PyErr_SetString(POINTERS.PyExc_BufferError, "not writable")
        return Cint(-1)
    else
        b.readonly[] = 1
    end

    # format
    if Utils.isflagset(flags, PyBUF_FORMAT)
        push!(c, info.format)
        b.format[] = pointer(info.format)
    else
        b.format[] = C_NULL
    end

    # shape
    if Utils.isflagset(flags, PyBUF_ND)
        shape = Py_ssize_t[info.shape...]
        push!(c, shape)
        b.shape[] = pointer(shape)
    else
        b.shape[] = C_NULL
    end

    # strides
    if Utils.isflagset(flags, PyBUF_STRIDES)
        strides = Py_ssize_t[info.strides...]
        push!(c, strides)
        b.strides[] = pointer(strides)
    elseif Utils.size_to_cstrides(info.itemsize, info.shape) == info.strides
        b.strides[] = C_NULL
    else
        PyErr_SetString(POINTERS.PyExc_BufferError, "not C contiguous and strides not requested")
        return Cint(-1)
    end

    # suboffsets
    if all(==(-1), info.suboffsets)
        b.suboffsets[] = C_NULL
    elseif Utils.isflagset(flags, PyBUF_INDIRECT)
        suboffsets = Py_ssize_t[info.suboffsets...]
        push!(c, suboffsets)
        b.suboffsets[] = pointer(suboffsets)
    else
        PyErr_SetString(POINTERS.PyExc_BufferError, "indirect array and suboffsets not requested")
        return Cint(-1)
    end

    # check contiguity
    if Utils.isflagset(flags, PyBUF_C_CONTIGUOUS)
        if Utils.size_to_cstrides(info.itemsize, info.shape) != info.strides
            PyErr_SetString(POINTERS.PyExc_BufferError, "not C contiguous")
            return Cint(-1)
        end
    end
    if Utils.isflagset(flags, PyBUF_F_CONTIGUOUS)
        if Utils.size_to_fstrides(info.itemsize, info.shape) != info.strides
            PyErr_SetString(POINTERS.PyExc_BufferError, "not Fortran contiguous")
            return Cint(-1)
        end
    end
    if Utils.isflagset(flags, PyBUF_ANY_CONTIGUOUS)
        if Utils.size_to_cstrides(info.itemsize, info.shape) != info.strides &&
           Utils.size_to_fstrides(info.itemsize, info.shape) != info.strides
            PyErr_SetString(POINTERS.PyExc_BufferError, "not contiguous")
            return Cint(-1)
        end
    end

    # internal
    cptr = Base.pointer_from_objref(c)
    PYJLBUFCACHE[cptr] = c
    b.internal[] = cptr

    # obj
    Py_IncRef(obj)
    b.obj[] = obj
    Cint(0)
end

function _pyjl_get_buffer(o::PyPtr, buf::Ptr{Py_buffer}, flags::Cint)
    num_ = PyObject_GetAttrString(o, "_jl_buffer_info")
    num_ == C_NULL && (PyErr_Clear(); PyErr_SetString(POINTERS.PyExc_BufferError, "not a buffer"); return Cint(-1))
    num = PyLong_AsLongLong(num_)
    Py_DecRef(num_)
    num == -1 && return Cint(-1)
    try
        f = PYJLMETHODS[num]
        x = PyJuliaValue_GetValue(o)
        return _pyjl_get_buffer_impl(o, buf, flags, x, f)::Cint
    catch exc
        @debug "error getting the buffer" exc
        PyErr_SetString(POINTERS.PyExc_BufferError, "some error occurred getting the buffer")
        return Cint(-1)
    end
end

function _pyjl_release_buffer(xo::PyPtr, buf::Ptr{Py_buffer})
    delete!(PYJLBUFCACHE, UnsafePtr(buf).internal[!])
    nothing
end

function _pyjl_reduce(self::PyPtr, ::PyPtr)
    v = _pyjl_serialize(self, PyNULL)
    v == PyNULL && return PyNULL
    args = PyTuple_New(1)
    args == PyNULL && (Py_DecRef(v); return PyNULL)
    err = PyTuple_SetItem(args, 0, v)
    err == -1 && (Py_DecRef(args); return PyNULL)
    red = PyTuple_New(2)
    red == PyNULL && (Py_DecRef(args); return PyNULL)
    err = PyTuple_SetItem(red, 1, args)
    err == -1 && (Py_DecRef(red); return PyNULL)
    f = PyObject_GetAttrString(self, "_jl_deserialize")
    f == PyNULL && (Py_DecRef(red); return PyNULL)
    err = PyTuple_SetItem(red, 0, f)
    err == -1 && (Py_DecRef(red); return PyNULL)
    return red
end

function _pyjl_serialize(self::PyPtr, ::PyPtr)
    try
        io = IOBuffer()
        serialize(io, PyJuliaValue_GetValue(self))
        b = take!(io)
        return PyBytes_FromStringAndSize(pointer(b), sizeof(b))
    catch e
        PyErr_SetString(POINTERS.PyExc_Exception, "error serializing this value")
        return PyNULL
    end
end

function _pyjl_deserialize(t::PyPtr, v::PyPtr)
    try
        ptr = Ref{Ptr{Cchar}}()
        len = Ref{Py_ssize_t}()
        err = PyBytes_AsStringAndSize(v, ptr, len)
        err == -1 && return PyNULL
        io = IOBuffer(unsafe_wrap(Array, Ptr{UInt8}(ptr[]), Int(len[])))
        x = deserialize(io)
        return PyJuliaValue_New(t, x)
    catch e
        PyErr_SetString(POINTERS.PyExc_Exception, "error deserializing this value")
        return PyNULL
    end
end

const _pyjlbase_name = "juliacall.ValueBase"
const _pyjlbase_type = fill(C.PyTypeObject())
const _pyjlbase_isnull_name = "_jl_isnull"
const _pyjlbase_callmethod_name = "_jl_callmethod"
const _pyjlbase_reduce_name = "__reduce__"
const _pyjlbase_serialize_name = "_jl_serialize"
const _pyjlbase_deserialize_name = "_jl_deserialize"
const _pyjlbase_methods = Vector{PyMethodDef}()
const _pyjlbase_as_buffer = fill(PyBufferProcs())

function init_jlwrap()
    empty!(_pyjlbase_methods)
    push!(_pyjlbase_methods,
        PyMethodDef(
            name = pointer(_pyjlbase_callmethod_name),
            meth = @cfunction(_pyjl_callmethod, PyPtr, (PyPtr, PyPtr)),
            flags = Py_METH_VARARGS,
        ),
        PyMethodDef(
            name = pointer(_pyjlbase_isnull_name),
            meth = @cfunction(_pyjl_isnull, PyPtr, (PyPtr, PyPtr)),
            flags = Py_METH_NOARGS,
        ),
        PyMethodDef(
            name = pointer(_pyjlbase_reduce_name),
            meth = @cfunction(_pyjl_reduce, PyPtr, (PyPtr, PyPtr)),
            flags = Py_METH_NOARGS,
        ),
        PyMethodDef(
            name = pointer(_pyjlbase_serialize_name),
            meth = @cfunction(_pyjl_serialize, PyPtr, (PyPtr, PyPtr)),
            flags = Py_METH_NOARGS,
        ),
        PyMethodDef(
            name = pointer(_pyjlbase_deserialize_name),
            meth = @cfunction(_pyjl_deserialize, PyPtr, (PyPtr, PyPtr)),
            flags = Py_METH_O | Py_METH_CLASS,
        ),
        PyMethodDef(),
    )
    _pyjlbase_as_buffer[] = PyBufferProcs(
        get = @cfunction(_pyjl_get_buffer, Cint, (PyPtr, Ptr{Py_buffer}, Cint)),
        release = @cfunction(_pyjl_release_buffer, Cvoid, (PyPtr, Ptr{Py_buffer})),
    )
    _pyjlbase_type[] = PyTypeObject(
        name = pointer(_pyjlbase_name),
        basicsize = sizeof(PyJuliaValueObject),
        # new = POINTERS.PyType_GenericNew,
        new = @cfunction(_pyjl_new, PyPtr, (PyPtr, PyPtr, PyPtr)),
        dealloc = @cfunction(_pyjl_dealloc, Cvoid, (PyPtr,)),
        flags = Py_TPFLAGS_BASETYPE | Py_TPFLAGS_HAVE_VERSION_TAG,
        weaklistoffset = fieldoffset(PyJuliaValueObject, 3),
        # getattro = POINTERS.PyObject_GenericGetAttr,
        # setattro = POINTERS.PyObject_GenericSetAttr,
        methods = pointer(_pyjlbase_methods),
        as_buffer = pointer(_pyjlbase_as_buffer),
    )
    o = POINTERS.PyJuliaBase_Type = PyPtr(pointer(_pyjlbase_type))
    if PyType_Ready(o) == -1
        PyErr_Print()
        error("Error initializing 'juliacall.ValueBase'")
    end
end

PyJuliaValue_IsNull(o::PyPtr) = UnsafePtr{PyJuliaValueObject}(o).value[] == 0

PyJuliaValue_GetValue(o::PyPtr) = PYJLVALUES[UnsafePtr{PyJuliaValueObject}(o).value[]]

PyJuliaValue_SetValue(o::PyPtr, @nospecialize(v)) = begin
    idx = UnsafePtr{PyJuliaValueObject}(o).value[]
    if idx == 0
        if isempty(PYJLFREEVALUES)
            push!(PYJLVALUES, v)
            idx = length(PYJLVALUES)
        else
            idx = pop!(PYJLFREEVALUES)
            PYJLVALUES[idx] = v
        end
        UnsafePtr{PyJuliaValueObject}(o).value[] = idx
    else
        PYJLVALUES[idx] = v
    end
    nothing
end

PyJuliaValue_New(t::PyPtr, @nospecialize(v)) = begin
    if PyType_IsSubtype(t, POINTERS.PyJuliaBase_Type) != 1
        PyErr_SetString(POINTERS.PyExc_TypeError, "Expecting a subtype of 'juliacall.ValueBase'")
        return PyNULL
    end
    o = PyObject_CallObject(t, PyNULL)
    o == PyNULL && return PyNULL
    PyJuliaValue_SetValue(o, v)
    return o
end


text:
# -------------------------------------------------------------------------------------------------------------------------------------------------------------
# Script creating monthly data files
# -------------------------------------------------------------------------------------------------------------------------------------------------------------

using DataFrames, CSV, XLSX, LinearAlgebra, Statistics

dir_io = "C:/Users/u0148308/data/comext/" # location of input/output (io)

# -------------------------------------------------------------------------------------------------------------------------------------------------------------

# timespan
years = string.(2015:2021)
months = lpad.(1:12, 2, '0')

# initial cleaning of the data
function initial_cleaning(year::String, month::String)

    path = dir_io * "raw/" * "full" * year * month * ".dat"
    df = CSV.read(path, DataFrame)

    # correct column types
    # Notes:
    #   - treat PERIOD as numeric to allow sorting
    transform!(df, [:DECLARANT, :PARTNER, :PERIOD] .=> ByRow(Int64), renamecols=false)

    # different additional PROD_BEC column from 2017 onwards
    col_strings = [:DECLARANT_ISO, :PARTNER_ISO, :TRADE_TYPE, :PRODUCT_NC, :PRODUCT_SITC, :PRODUCT_CPA2002, :PRODUCT_CPA2008, :PRODUCT_CPA2_1,
                    :PRODUCT_BEC, :PRODUCT_SECTION, :FLOW, :STAT_REGIME]
    col_strings = ifelse(year in ["2015", "2016"], col_strings, [col_strings; :PRODUCT_BEC5])
    transform!(df,  col_strings .=> ByRow(string), renamecols=false)

    transform!(df, [:VALUE_IN_EUROS, :QUANTITY_IN_KG, :SUP_QUANTITY] .=> ByRow(Float64), renamecols=false)

    # rename some indicators
    df[:, :TRADE_TYPE] .= ifelse.(df[:, :TRADE_TYPE] .== "I", "intra", "extra")
    df[:, :FLOW] .= ifelse.(df[:, :FLOW] .== "1", "imports", "exports")
    df[:, :SUPP_UNIT] .= ifelse.(ismissing.(df[:, :SUPP_UNIT]), missing, string.(df[:, :SUPP_UNIT])) # to make sure type is of Union{Missing, String}

    # missing data
    # Notes:
    #   - if value or quanity is reported as 0 use missing instead
    transform!(df, [:VALUE_IN_EUROS, :QUANTITY_IN_KG, :SUP_QUANTITY] .=> ByRow(x->ifelse(x == 0.0, missing, x)), renamecols=false)

    # aggregate over STAT_REGIME, SUPP_UNIT
    # Notes:
    #   - only lose few rows ~1%
    #   - disregard SUPP_UNIT and SUP_QUANTITY for the time being
    cols_grouping = ["DECLARANT_ISO", "PARTNER_ISO", "TRADE_TYPE", "PRODUCT_NC", "FLOW", "PERIOD"]
    gdf = groupby(df, cols_grouping)
    df = combine(gdf, [:VALUE_IN_EUROS, :QUANTITY_IN_KG] .=> sum, renamecols=false)

    return df
end

function pct_change(input::AbstractVector)
    [i == 1 ? missing : (input[i]-input[i-1])/input[i-1]*100 for i in eachindex(input)]
end

function mom_change(period::AbstractVector, input::AbstractVector)
    [i == 1 ? missing : period[i]-period[i-1] != 1 ? missing : (input[i]-input[i-1])/input[i-1]*100 for i in eachindex(input)]
end

g(x) = x[end-1:end]

periods = repeat(years, inner=12) .* repeat(months, outer=length(years))
periods = parse.(Int64, periods) # convert to Int64 (actually needed? think we can do all with String)

# create monthly files
# Notes:
#   - import data for 3 months, compute prices and shares
#       + could use 2 months only but then we can only compute MOM. In this setting we obtain also UNIT_PRICE_CHANGE between Jan and Mar,
#         this might be useful later if we realize MOM is too restrictive
#   - subset to middle month (i.e. Feb if Jan, Feb, Mar) and export
#       + means we lose first and last month (i.e. cannot compute changes since no previous/following data)

for (i, time) in enumerate(periods[1:end-2])

    year = string.(time)[1:4]
    month = g.(string.(periods[[i, i+1, i+2]]))
    sub_month = month[2]

    df = DataFrame(DECLARANT_ISO=String[], PARTNER_ISO=String[], TRADE_TYPE=String[], PRODUCT_NC=String[], FLOW=String[], PERIOD=Int64[],
                VALUE_IN_EUROS=Union{Missing, Float64}[], QUANTITY_IN_KG=Union{Missing, Float64}[])

    append!(df, initial_cleaning(year, month[1]))
    append!(df, initial_cleaning(year, sub_month))
    append!(df, initial_cleaning(year, month[3]))

    sort!(df)

    # -------------------------------------------------------------------------------------------------------------------------------------------------------------
    # Unit prices
    # -------------------------------------------------------------------------------------------------------------------------------------------------------------

    # compute UNIT_PRICE
    transform!(df, [:VALUE_IN_EUROS, :QUANTITY_IN_KG] => ByRow((v,q) -> v/q) => :UNIT_PRICE)

    # compute UNIT_PRICE_CHANGE
    cols_grouping = ["DECLARANT_ISO", "PARTNER_ISO", "PRODUCT_NC", "FLOW"]
    gdf = groupby(df, cols_grouping)
    df = transform(gdf, :UNIT_PRICE => pct_change => :UNIT_PRICE_CHANGE)

    # compute MOM_PRICE_CHANGE
    gdf = groupby(df, cols_grouping)
    df = transform(gdf, [:PERIOD, :UNIT_PRICE] => mom_change => :MOM)

    # export data
    df_export = subset(df, :PERIOD => ByRow(x-> string(x) == year * sub_month))
    CSV.write(dir_io * "clean/" * "df_prices_" * year * sub_month * ".csv", df_export)

    println(" ✓ Price data for $sub_month/$year were successfully created.")

    # -------------------------------------------------------------------------------------------------------------------------------------------------------------
    # Share of Belgium
    # -------------------------------------------------------------------------------------------------------------------------------------------------------------

    # compute EU_WORLD (EU imports/exports to the entire world)
    cols_grouping = ["PRODUCT_NC", "FLOW", "PERIOD"]
    gdf = groupby(df, cols_grouping)
    df_sum_EU_WORLD = combine(gdf, [:VALUE_IN_EUROS, :QUANTITY_IN_KG] .=> (x -> sum(skipmissing(x))), renamecols=false)
    df_sum_EU_WORLD.TRADE_TYPE .= "total"

    # compute EU_EXTRA and EU_INTRA (EU imports/exports to extra/intra EU countries)
    cols_grouping = ["TRADE_TYPE", "PRODUCT_NC", "FLOW", "PERIOD"]
    gdf = groupby(df, cols_grouping)
    df_sum_EU_EXTRA_INTRA = combine(gdf, [:VALUE_IN_EUROS, :QUANTITY_IN_KG] .=> (x -> sum(skipmissing(x))), renamecols=false)

    # compute BE_WORLD (BE imports/exports to the entire world)
    df_BE = subset(df, :DECLARANT_ISO => ByRow(x -> x == "BE"))
    cols_grouping = ["PRODUCT_NC", "FLOW", "PERIOD"]
    gdf = groupby(df_BE, cols_grouping)
    df_sum_BE_WORLD = combine(gdf, [:VALUE_IN_EUROS, :QUANTITY_IN_KG] .=> (x -> sum(skipmissing(x))), renamecols=false)
    df_sum_BE_WORLD.TRADE_TYPE .= "total"

    # compute BE_EXTRA and BE_INTRA (BE imports/exports to extra/intra EU countries)
    cols_grouping = ["TRADE_TYPE", "PRODUCT_NC", "FLOW", "PERIOD"]
    gdf = groupby(df_BE, cols_grouping)
    df_sum_BE_EXTRA_INTRA = combine(gdf, [:VALUE_IN_EUROS, :QUANTITY_IN_KG] .=> (x -> sum(skipmissing(x))), renamecols=false)

    # merge tables
    df_sum_EU = vcat(df_sum_EU_EXTRA_INTRA, df_sum_EU_WORLD)
    df_sum_BE = vcat(df_sum_BE_EXTRA_INTRA, df_sum_BE_WORLD)

    # join, rename and compute shares in percentages (?)
    df_join = leftjoin(df_sum_EU, df_sum_BE, on=[:TRADE_TYPE, :PRODUCT_NC, :FLOW, :PERIOD], makeunique=true)
    rename!(df_join, ["TRADE_TYPE", "PRODUCT_NC", "FLOW", "PERIOD", "VALUE_EU", "QUANITY_EU", "VALUE_BE", "QUANITY_BE"])
    df_join.VALUE_SHARE = df_join.VALUE_BE ./ df_join.VALUE_EU .* 100
    df_join.QUANTITY_SHARE = df_join.QUANITY_BE ./ df_join.QUANITY_EU .* 100

    # possibly need to sort before grouping, seems that grouping sorts incorrectly sometimes?
    sort!(df_join, [:PRODUCT_NC, :FLOW, :TRADE_TYPE, :PERIOD])

    # compute SHARE_CHANGE
    cols_grouping = ["TRADE_TYPE", "PRODUCT_NC", "FLOW"]
    gdf = groupby(df_join, cols_grouping)
    df_share_BE = transform(gdf, [:VALUE_SHARE, :QUANTITY_SHARE] .=> pct_change .=> [:VALUE_SHARE_CHANGE, :QUANTITY_SHARE_CHANGE])

    # compute MOM_CHANGE
    gdf = groupby(df_share_BE, cols_grouping)
    df_share_BE = transform(gdf, [:PERIOD, :VALUE_SHARE] => mom_change => :MOM_VALUE,
                                 [:PERIOD, :QUANTITY_SHARE] => mom_change => :MOM_QUANTITY)

    # export data
    df_export = subset(df_share_BE, :PERIOD => ByRow(x-> string(x) == year * sub_month))
    CSV.write(dir_io * "clean/" * "df_share_BE_" * year * sub_month * ".csv", df_export)

    println(" ✓ Belgian shares data for $sub_month/$year were successfully created. \n")

end


text:
function setindex!(cur_matrix::FourBlockMatrix,cur_value::Any,i::Int,j::Int)
  n = cur_matrix.n
  m = cur_matrix.m

  @assert i <= n + m
  @assert j <= n + m

  if i <= n
    if j <= n
      cur_matrix.A[i,j] = cur_value
    else
      cur_matrix.B[i,j-n] = cur_value
    end
  else
    if j <= n
      cur_matrix.C[i-n,j] = cur_value
    else
      cur_matrix.D[i-n,j-n] = cur_value
    end
  end
end


text:
using Test

@testset "Nesterov accel. ($R)" for R in [Float32, Float64]

using LinearAlgebra
using ProximalAlgorithms: NesterovAcceleration, update!

n, k = 20, 10

Q = randn(R, n, k)
H = Q*Q' + I
l = randn(R, n)

f(x) = dot(x, H*x)/2 + dot(x, l)

x_star = -H\l
f_star = f(x_star)
norm_x_star = norm(x_star)

Lip = opnorm(H)
gamma = 1/Lip
x = zeros(n)
res = zeros(n)
x_prev = zeros(n)
res_prev = zeros(n)

acc = NesterovAcceleration(x)

mul!(res, H, x)
res .+= l
res .*= gamma

d = acc * res

for it in 1:500
    # store iterate and residual for the operator update later
    res_prev .= res
    x_prev .= x

    # compute accelerated direction
    mul!(d, acc, res)

    # update iterate
    x .-= d

    # compute new residual
    mul!(res, H, x)
    res .+= l
    res .*= gamma

    # update operator
    update!(acc, x - x_prev, res - res_prev)

    # test that iterates satisfy Thm 4.4 from Beck, Teboulle (2009)
    @test f(x) - f_star <= 2/(gamma * (it+1)^2) * norm_x_star^2
end

end


text:
using TimeZones
using Test

include("../src/Cn2An.jl")

function run_an2cn_core_ten_thousand_times()
    t_start = now()
    for _ in 1:1000000
        @test Cn2An.an2cn(9876543298765432) == "九千八百七十六万五千四百三十二亿九千八百七十六万五千四百三十二"
    end
    t_end = now()
    println(t_end-t_start)
end

run_an2cn_core_ten_thousand_times()


text:
module BSR

try
    include("paths.jl") # This file is generated by the package build script
catch
    error("Run Pkg.build(\"BSR\") first")
end

include("common.jl")
include("knot_set.jl")
include("target.jl")
include("prep_state.jl")
include("calc.jl")
include("load_matrices.jl")
include("load_states.jl")
include("load_transitions.jl")

end # module