text:
"""sys_fr = freqresp(sys, w)

Evaluate the frequency response of a linear system

`w -> C*((iw*im -A)^-1)*B + D`

of system `sys` over the frequency vector `w`."""
function freqresp(sys::LTISystem, w_vec::AbstractVector{<:Real})
    # Create imaginary freq vector s
    if !iscontinuous(sys)
        Ts = sys.Ts == -1 ? 1.0 : sys.Ts
        s_vec = exp.(w_vec*(im*Ts))
    else
        s_vec = im*w_vec
    end
    if isa(sys, StateSpace)
        sys = _preprocess_for_freqresp(sys)
    end
    ny,nu = noutputs(sys), ninputs(sys)
    [evalfr(sys[i,j], s)[] for s in s_vec, i in 1:ny, j in 1:nu]
end

# Implements algorithm found in:
# Laub, A.J., "Efficient Multivariable Frequency Response Computations",
# IEEE Transactions on Automatic Control, AC-26 (1981), pp. 407-408.
function _preprocess_for_freqresp(sys::StateSpace)
    if isempty(sys.A) # hessfact does not work for empty matrices
        return sys
    end
    Tsys = numeric_type(sys)
    TT = promote_type(typeof(zero(Tsys)/norm(one(Tsys))), Float32)

    A, B, C, D = sys.A, sys.B, sys.C, sys.D
    F = hessenberg(A)
    T = F.Q
    P = C*T
    Q = T\B # TODO Type stability? # T is unitary, so mutliplication with T' should do the trick
    # FIXME; No performance improvement from Hessienberg structure, also weired renaming of matrices
    StateSpace(F.H, Q, P, D, sys.Ts)
end

#_preprocess_for_freqresp(sys::TransferFunction) = sys.matrix
#function _preprocess_for_freqresp(sys::TransferFunction)
#    map(sisotf -> _preprocess_for_freqresp(sisotf), sys.matrix)
#end

#_preprocess_for_freqresp(sys::SisoTf) = sys

"""
`evalfr(sys, x)` Evaluate the transfer function of the LTI system sys
at the complex number s=x (continuous-time) or z=x (discrete-time).

For many values of `x`, use `freqresp` instead.
"""
function evalfr(sys::StateSpace{T0}, s::Number) where {T0}
    T = promote_type(T0, typeof(one(T0)*one(typeof(s))/(one(T0)*one(typeof(s)))))
    try
        R = s*I - sys.A
        sys.D + sys.C*((R\sys.B)::Matrix{T})  # Weird type stability issue
    catch
        fill(convert(T, Inf), size(sys))
    end
end

function evalfr(G::TransferFunction{<:SisoTf}, s::Number)
    map(m -> evalfr(m,s), G.matrix)
end

"""
`F(s)`, `F(omega, true)`, `F(z, false)`

Notation for frequency response evaluation.
- F(s) evaluates the continuous-time transfer function F at s.
- F(omega,true) evaluates the discrete-time transfer function F at exp(i*Ts*omega)
- F(z,false) evaluates the discrete-time transfer function F at z
"""
function (sys::TransferFunction)(s)
    evalfr(sys,s)
end

function (sys::TransferFunction)(z_or_omega::Number, map_to_unit_circle::Bool)
    @assert !iscontinuous(sys) "It makes no sense to call this function with continuous systems"
    if map_to_unit_circle
        isreal(z_or_omega) ? evalfr(sys,exp(im*z_or_omega.*sys.Ts)) : error("To map to the unit circle, omega should be real")
    else
        evalfr(sys,z_or_omega)
    end
end

function (sys::TransferFunction)(z_or_omegas::AbstractVector, map_to_unit_circle::Bool)
    @assert !iscontinuous(sys) "It makes no sense to call this function with continuous systems"
    vals = sys.(z_or_omegas, map_to_unit_circle)# evalfr.(sys,exp.(evalpoints))
    # Reshape from vector of evalfr matrizes, to (in,out,freq) Array
    nu,ny = size(vals[1])
    [v[i,j]  for v in vals, i in 1:nu, j in 1:ny]
end

"""`mag, phase, w = bode(sys[, w])`

Compute the magnitude and phase parts of the frequency response of system `sys`
at frequencies `w`

`mag` and `phase` has size `(length(w), ny, nu)`"""
function bode(sys::LTISystem, w::AbstractVector; unit_circle_threshold = 1e-5, zero_threshold = 1e-8)
    resp = freqresp(sys, w)

    # use the number of integrators in the system to estimate initial phase at w=0
    count_unit_circle = x -> count(y -> abs(abs(y) - 1.0) < unit_circle_threshold, x)
    count_zeros = x -> count(y -> abs(y) < zero_threshold, x)

    phase0 = if isa(sys, DelayLtiSystem)
        0.0
    else
        zpk_sys = zpk(sys)
        if sys.Ts == 0.0
            pi/2 * (count_zeros(zpk_sys.matrix[1].z) - count_zeros(zpk_sys.matrix[1].p))
        else
            pi/2 * (count_unit_circle(zpk_sys.matrix[1].z) - count_unit_circle(zpk_sys.matrix[1].p))
        end
    end
    phase = rad2deg.(unwrap!(angle.(resp), 1, init = phase0))
    return abs.(resp), phase, w
end
bode(sys::LTISystem; kwargs...) = bode(sys, _default_freq_vector(sys, Val{:bode}()); kwargs...)

"""`re, im, w = nyquist(sys[, w])`

Compute the real and imaginary parts of the frequency response of system `sys`
at frequencies `w`

`re` and `im` has size `(length(w), ny, nu)`"""
function nyquist(sys::LTISystem, w::AbstractVector)
    resp = freqresp(sys, w)
    return real(resp), imag(resp), w
end
nyquist(sys::LTISystem) = nyquist(sys, _default_freq_vector(sys, Val{:nyquist}()))

"""`sv, w = sigma(sys[, w])`

Compute the singular values of the frequency response of system `sys` at
frequencies `w`

`sv` has size `(length(w), max(ny, nu))`"""
function sigma(sys::LTISystem, w::AbstractVector)
    resp = freqresp(sys, w)
    sv = dropdims(mapslices(svdvals, resp, dims=(2,3)),dims=3)
    return sv, w
end
sigma(sys::LTISystem) = sigma(sys, _default_freq_vector(sys, Val{:sigma}()))

function _default_freq_vector(systems::Vector{<:LTISystem}, plot)
    min_pt_per_dec = 60
    min_pt_total = 200
    bounds = map(sys -> _bounds_and_features(sys, plot)[1], systems)
    w1 = minimum(minimum.(bounds))
    w2 = maximum(maximum.(bounds))

    nw = round(Int, max(min_pt_total, min_pt_per_dec*(w2 - w1)))
    return exp10.(range(w1, stop=w2, length=nw))
end
_default_freq_vector(sys::LTISystem, plot) = _default_freq_vector(
        [sys], plot)

function _bounds_and_features(sys::LTISystem, plot::Val)
    # Get zeros and poles for each channel
    if !isa(plot, Val{:sigma})
        zs, ps = zpkdata(sys)
        # Compose vector of all zs, ps, positive conjugates only.
        zpType = promote_type(eltype(eltype(zs)), eltype(eltype(ps)))
        zp = vcat(zpType[], zs..., ps...) # Emty vector to avoid type unstable vcat()
        zp = zp[imag(zp) .>= 0.0]
    else
         # For sigma plots, use the MIMO poles and zeros
         zp = [tzero(sys); pole(sys)]
    end
    # Get the frequencies of the features, ignoring low frequency dynamics
    fzp = log10.(abs.(zp))
    fzp = fzp[fzp .> -4]
    fzp = sort!(fzp)
    # Determine the bounds on the frequency vector
    if !isempty(fzp)
        w1 = floor(fzp[1] - 0.2)
        w2 = ceil(fzp[end] + 0.2)
        # Expand the range for nyquist plots
        if plot isa Val{:nyquist}
            w1 -= 1.0
            w2 += 1.0
        end
    else
        w1 = 0.0
        w2 = 2.0
    end
    if !iscontinuous(sys) # Do not draw above Nyquist freq for disc. systems
        w2 = min(w2, log10(π/sys.Ts))
    end
    return [w1, w2], zp
end


text:
export Event, matrix, vector, isfilled, almost_zero, isalmostzero

mutable struct Event
    user::Int
    item::Int
    value::Float64 # e.g. rating, 0/1
end

accepted_types = Union{Nothing, Missing, Number}

function matrix(m::Int, n::Int;
                type::Type{<:accepted_types}=Union{Nothing, Float64}, initializer=nothing)
    Array{type}(initializer, m, n)
end

function vector(m::Int;
                type::Type{<:accepted_types}=Union{Nothing, Float64}, initializer=nothing)
    Array{type}(initializer, m)
end

function isfilled(a::AbstractArray; by_value=nothing)
    by_value ∉ Set(a)
end

almost_zero = 1e-256

function isalmostzero(x::Number)
    x <= almost_zero
end

ext: jl


text:
using DoublePendulum
using Test

@testset "DoublePendulum" begin
    include("data.jl")
end


text:
# This file is part of IntegerSequences.
# Copyright Peter Luschny. License is MIT.

(@__DIR__) ∉ LOAD_PATH && push!(LOAD_PATH, (@__DIR__))

module GeneralizedFibonacci
using Nemo, SeqUtils

export ModuleGeneralizedFibonacci
export Multinacci, V309896, V006053, V188021, V231181

"""

* The classical Fibonacci numbers are exported from the module Fibonacci. See I000045, F000045, L000045, V000045, R000045 and is000045.

* Fibonacci(n) is defined as the number of compositions of n with no part equal to 1. They are the special case Fibonacci(n) = Multinacci(2, n).

```
[m\n] 0  1  2   3   4   5    6    7     8     9    10    11
----------------------------------------------------------------
[0]   1, 0, 0,  0,  0,  0,   0,   0,    0,    0,    0,    0, ...
[1]   1, 1, 1,  1,  1,  1,   1,   1,    1,    1,    1,    1, ...
[2]   1, 1, 2,  3,  5,  8,  13,  21,   34,   55,   89,  144, ...
[3]   1, 1, 3,  4,  9, 14,  28,  47,   89,  155,  286,  507, ...
[4]   1, 1, 4,  5, 14, 20,  48,  75,  165,  274,  571,  988, ...
[5]   1, 1, 5,  6, 20, 27,  75, 110,  275,  429, 1001, 1637, ...
[6]   1, 1, 6,  7, 27, 35, 110, 154,  429,  637, 1638, 2548, ...
```

*  Multinacci, V309896, V006053, V188021, V231181
"""
const ModuleGeneralizedFibonacci = ""

"""

Return the ``n``-th term of the ``m``-th Multinacci sequence.
"""
function Multinacci(m::Int, n::Int)::fmpz
    haskey(CacheMultinacci, (m, n)) && return CacheMultinacci[(m, n)]
    n <  0 && return fmpz(0); n == 0 && return fmpz(1); m == 0 && return fmpz(m^n)
    a = sum((-1)^j*binomial(ZZ(m-1-j), ZZ(j ) )*Multinacci(m, n-1-2*j) for j ∈ 0:div(m-1,2))
    b = sum((-1)^j*binomial(ZZ(m-1-j), ZZ(j+1))*Multinacci(m, n-2-2*j) for j ∈ 0:div(m-2,2))
    r = a + b
    CacheMultinacci[(m, n)] = r
end

"""

Return the ``k``-th term of the ``n``-th Multinacci sequence.
"""
V309896(n, k) = Multinacci(n, k)

const CacheMultinacci = Dict{Tuple{Int, Int}, fmpz}()

"""

Return the ``n``-th term of the third Multinacci sequence.
"""
V006053(n) = Multinacci(3, n)

"""

Return the ``n``-th term of the fourth Multinacci sequence.
"""
V188021(n) = Multinacci(4, n)

"""

Return the ``n``-th term of the fiveth Multinacci sequence.
"""
V231181(n) = Multinacci(5, n)

#START-TEST-########################################################

using Test

function test()
    @testset "Multinacci" begin
        @test Multinacci(0, 0) == 1
        @test [Multinacci(2, n) for n ∈ 0:9] == [1, 1, 2, 3, 5,  8, 13, 21, 34,  55]
        @test [V006053(n) for n ∈ 0:9] == [1, 1, 3, 4, 9, 14, 28, 47, 89, 155]
        @test [V188021(n) for n ∈ 0:9] == [1, 1, 4, 5, 14, 20, 48, 75, 165, 274]
        @test [V231181(n) for n ∈ 0:9] == [1, 1, 5, 6, 20, 27, 75, 110, 275,429]
        @test [Multinacci(6, n) for n ∈ 0:9] == [1, 1, 6, 7, 27, 35, 110, 154,429,637]
        @test [Multinacci(7, n) for n ∈ 0:9] == [1, 1, 7, 8, 35, 44, 154, 208,637,910]
        @test [Multinacci(8, n) for n ∈ 0:9] == [1, 1, 8, 9, 44, 54,208,273,910, 1260]
    end
end

function demo()
    for n ∈ 0:9
        print("$n ")
        [Multinacci(n, k) for k ∈ 0:9] |> Println
    end
end

function perf()
    @time (for n ∈ 0:1000 Multinacci(3, n) end)
end

function main()
    test()
    demo()
    perf()
end

main()

end # module

max_stars_repo_name: OpenLibMathSeq/Sequences

max_issues_repo_name: OpenLibMathSeq/Sequences


text:
# Part 1
horizontal, depth = 0, 0

for line in eachline("input.txt")
    direction, distance = split(line, " ")
    distance = parse(Int, distance)
    if direction == "forward"
        horizontal += distance
    elseif direction == "down"
        depth += distance
    elseif direction == "up"
        depth -= distance
    end
end

@info "Final horizontal position times final depth = $(horizontal*depth)" horizontal depth

# Part 2
horizontal, depth, aim = 0, 0, 0

for line in eachline("input.txt")
    direction, distance = split(line, " ")
    distance = parse(Int, distance)
    if direction == "forward"
        horizontal += distance
        depth += aim*distance
    elseif direction == "down"
        aim += distance
    elseif direction == "up"
        aim -= distance
    end
end

@info "Final horizontal position times final depth = $(horizontal*depth)" horizontal depth aim


text:
using Test, SciMLBase, SciMLBase.EnsembleAnalysis

EA = SciMLBase.EnsembleAnalysis

# tests for https://github.com/SciML/DifferentialEquations.jl/issues/731
# make sure integer inputs work
@test all(EA.componentwise_mean([[1,1],[2,2],[3,3]]) .≈ [2.0,2.0])
m,v = EA.componentwise_meanvar([[1,1],[2,2],[3,3]])
@test all(m .≈ 2.0)
@test all(v .≈ 1.0)
mx,my,C = EA.componentwise_meancov([[1,1],[2,2],[3,3]], [[3,3],[2,2],[1,1]])
@test all(mx .≈ 2.0)
@test all(my .≈ 2.0)
@test all(C .≈ -1.0)
mx,my,C = EA.componentwise_weighted_meancov([[1,1],[2,2],[3,3]], [[3,3],[2,2],[1,1]], [[1,1],[2,2],[1,1]])
@test all(mx .≈ 2.0)
@test all(my .≈ 2.0)
@test all(C .≈ -.8)


text:
using ParallelKMeans, RDatasets, Plots
using Statistics, MAT, Clustering, JLD2
using FileIO
using Distances, Clustering, Statistics
k = 10

function clustering_yearwise(whichyear, k)
# whichyear = "2035"
        year = string("NAT_",whichyear, "_Generation")
        data_ip =  Array{Float64}(undef, 8760, 4)
        fname = string("C:\\Users\\djaykuma\\OneDrive - Energyville\\Freq_TNEP_paper\\MATLAB\\NAT_Results\\",year,"\\output.mat")
        vars = matread(fname)
        fname1 = string("C:\\Users\\djaykuma\\OneDrive - Energyville\\Freq_TNEP_paper\\MATLAB\\NAT_Results\\",year,"\\output_MC.mat")
        vars1 = matread(fname1)
        fname2 = string("C:\\Users\\djaykuma\\OneDrive - Energyville\\Freq_TNEP_paper\\MATLAB\\NAT_Results\\wind_sample.mat")
        vars2 = matread(fname2)
        for i = 1:8760
                     # if vars["Heq"][i] < 0.5
                     #         vars["Heq"][i] = 0.5
                     # end
                     data_ip[i,1] = vars["Heq"][i]/(maximum(vars["Heq"])) #inertia
                     data_ip[i,2] = vars1["MC1"][i]/(maximum(vars1["MC1"])) # gen cost
                     data_ip[i,3] = vars2["winddata"][i]/maximum(vars2["winddata"]) #wind gen
                     data_ip[i,4] = vars["Ptotal"][i]/(maximum(vars["Ptotal"])) #total gen and demand
        end
        D = pairwise(Euclidean(), data_ip', data_ip')
        result = Clustering.kmedoids(D, k)

        # kk = [mean(Clustering.silhouettes(Clustering.kmedoids(D, k), D)) for k = 2:50]

        # result.medoids[result.assignments]
        # result_kmeans = Clustering.kmeans(data_ip', k)
        # b1 = [Clustering.kmedoids(D, i).totalcost for i = 2:100]
        # plot(b1)
        # b2 = [ParallelKMeans.kmeans(data_ip', i, n_threads=1; tol=1e-6, max_iters=300, verbose=false).totalcost for i = 2:100]
        # plot(b2)
        # b3 = [ParallelKMeans.kmeans(data_ip', i; tol=1e-6, max_iters=300, verbose=false).totalcost for i = 2:10]
        # plot(b3)
        return result
end

function clustering_output(result, k)
        weights =[]
        scenarios = []
        for i = 1:k
                wgt = count(x->x == i, result.assignments)/8760
                push!(weights, wgt)
                scenarios = result.medoids
        end
        return weights, scenarios
end

year = ["2025","2030","2035","2040","2045","2050"]
for i = 1:length(year)
        whichyear = year[i]
        result = clustering_yearwise(whichyear, k)
        weights, scenarios = clustering_output(result, k)
        @save string("cluster_",whichyear,".jld2") scenarios weights
end

# scatter(data_ip[:,1], data_ip[:,2], color=:lightrainbow,marker_z=result.assignments, legend=false)
# plot(data_ip[:,2])

function collect_scenario(year)
        scenario = []
        weight = []
        for i = 1:length(year)
                file = string("cluster_", year[i], ".jld2")
                var = load(file)
                append!(scenario, var["scenarios"])
                append!(weight, var["weights"])
        end
        return scenario, weight
end
#
scenario, weight = collect_scenario(year)
#
@save "cluster.jld2" scenarios = scenario weights = weight
# example 1
# X1 = [rand(1, 10); rand(1, 10) ; rand(1, 10); rand(1, 10)]
# # X1 = round.(Int, X)
# result = Clustering.kmeans(X1, 5)
# D = pairwise(Euclidean(), X1, X1)
# kk = Clustering.silhouettes(result, D)
# mean(kk)
# # or
# b = [ParallelKMeans.kmeans(X1, i, n_threads=1; tol=1e-6, max_iters=300, verbose=false).totalcost for i = 2:10]
# plot(b)

# example 2
# iris = dataset("datasets", "iris"); # load the data
# features = collect(Matrix(iris[:, 1:4])'); # features to use for clustering
# result = kmeans(features, 5); # run K-means for the 3 clusters
# scatter(iris.PetalLength, iris.PetalWidth, marker_z=result.assignments,
#       color=:lightrainbow, legend=false)
# D = pairwise(Euclidean(), features, features)
# kk = silhouettes(result, D)
# mean(kk) # score reduces above certain value
# # or
# b = [ParallelKMeans.kmeans(features, i, n_threads=1; tol=1e-6, max_iters=300, verbose=false).totalcost for i = 2:20]

# using ParallelKMeans, RDatasets, Plots
# # load the data
# iris = dataset("datasets", "iris");
# # features to use for clustering
# features = collect(Matrix(iris[:, 1:4])');
# # various artifacts can be accessed from the result i.e. assigned labels, cost value etc
# result = kmeans(features, 3);
# # plot with the point color mapped to the assigned cluster index
# scatter(iris.PetalLength, iris.PetalWidth, marker_z=result.assignments,
#         color=:lightrainbow, legend=false)
# # Single Thread Implementation of Lloyd's Algorithm
# b = [ParallelKMeans.kmeans(features, i, n_threads=1; tol=1e-6, max_iters=300, verbose=false).totalcost for i = 2:10]
# # Multi-threaded Implementation of Lloyd's Algorithm by default
# c = [ParallelKMeans.kmeans(features, i; tol=1e-6, max_iters=300, verbose=false).totalcost for i = 2:10]

# scenarios = scenario_500[1:10]
# weights = [1,1,1,1,1,1,1,1,1,1]
# @save "cluster_2025.jld2" scenarios = scenarios weights = weights
# @save "cluster_2030.jld2" scenarios = scenarios weights = weights
# @save "cluster_2035.jld2" scenarios = scenarios weights = weights
# @save "cluster_2040.jld2" scenarios = scenarios weights = weights
# @save "cluster_2045.jld2" scenarios = scenarios weights = weights
# @save "cluster_2050.jld2" scenarios = scenarios weights = weights


text:
using CUSPARSE
using CUDArt
using Base.Test

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

@testset "gthr and gthrz" begin
    @testset for elty in [Float32,Float64,Complex64,Complex128]
        x = sparsevec(rand(1:m,k), rand(elty,k), m)
        y = rand(elty,m)
        @testset "gthr!" begin
            d_x = CudaSparseVector(x)
            d_y = CudaArray(y)
            d_y = CUSPARSE.gthr!(d_x,d_y,'O')
            h_x = to_host(d_x)
            @test h_x ≈ SparseVector(m,x.nzind,y[x.nzind])
        end

        @testset "gthr" begin
            d_x = CudaSparseVector(x)
            d_y = CudaArray(y)
            d_z = CUSPARSE.gthr(d_x,d_y,'O')
            h_z = to_host(d_z)
            @test h_z ≈ SparseVector(m,x.nzind,y[x.nzind])
        end

        @testset "gthrz!" begin
            d_x = CudaSparseVector(x)
            d_y = CudaArray(y)
            d_x,d_y = CUSPARSE.gthrz!(d_x,d_y,'O')
            h_x = to_host(d_x)
            h_y = to_host(d_y)
            @test h_x ≈ SparseVector(m,x.nzind,y[x.nzind])
            y[x.nzind] = zero(elty)
            @test h_y ≈ y
        end

        @testset "gthrz" begin
            d_x = CudaSparseVector(x)
            d_y = CudaArray(y)
            d_z,d_w = CUSPARSE.gthrz(d_x,d_y,'O')
            h_w = to_host(d_w)
            h_z = to_host(d_z)
            @test h_z ≈ SparseVector(m,x.nzind,y[x.nzind])
            y[x.nzind] = zero(elty)
            @test h_w ≈ y
        end
    end
end


text:
using simulate_quantum_memory

# Setup this simulation instance
timestamp  = convert(Int64, floor(time()*1000));
RATE       = float(ARGS[1]);
MAX_T      = convert(Int64, float(ARGS[2]));
NUM_GRID   = convert(Int64, float(ARGS[3]));
NUM_TRIALS = convert(Int64, float( ARGS[4]));

# Setup I/O for target and ffwd distributions
target_filename = string("target-", RATE, "-", MAX_T, "-", timestamp, ".txt")
target_stream   = open(target_filename, "w")
write(target_stream, "id, time, dt, weight\n")
ffwd_filename   = string("ffwd-", RATE, "-", MAX_T, "-", timestamp, ".txt")
ffwd_stream     = open(ffwd_filename, "w")
write(ffwd_stream, "time,dt,weight\n")

print("Simulate Fast Forward testing data\n");
print("==================================\n\n");
print("Starting up simulation...\n");

# Reseed the RNGs
srand(timestamp);
# Setup the [[15, 1, 3]] code
code = RMCode();
print("+ [[15, 1, 3]] code initialized\n")
# Setup a model with error rate and a top-out at MAX_T cycles
model = Model(RATE, NUM_GRID);
print("+ Single-angle dephase model initialized\n")
print("Setup complete.\n\n")

################
# Target model #
################
print("Starting up $NUM_TRIALS trials of `Target Model`...\n")

trials = 0;
runs   = 0;
while trials < NUM_TRIALS

	(time, dt, W_X, W_Z, DT) = single_angle(model, code, MAX_T;
	                                        adaptive=true,
																					simulate_decoder=false,
																					report=true);
	W_X_locs     = [1; nonzerolocs(W_X)];

	wick_lengths = [1; W_X_locs[2:end] - W_X_locs[1:end-1]];
	wick_dts     = DT[W_X_locs];
	wick_weights = sum_to_right_edge(W_Z, W_X_locs);

	for k=1:length(W_X_locs)
		write(target_stream, string(runs));            write(target_stream, ", ");
		write(target_stream, string(wick_lengths[k])); write(target_stream, ", ");
		write(target_stream, string(wick_dts[k]));     write(target_stream, ", ");
		write(target_stream, string(wick_weights[k])); write(target_stream, "\n");
	end
	write(target_stream, "\n");

	trials = trials + length(W_X_locs);
	runs = runs + 1;
end
print("Wrote $trials trials and $runs runs of `Target Model`.\n\n")

################
# Fast Forward #
################
print("Starting up $NUM_TRIALS trials of `Fast Forward`...\n")

for k=1:NUM_TRIALS
  dt = rand()*pi/8;
  (wick_length, spark, wick_weight) = fast_forward(code, model, dt;
																									 update_p_theta = false);
	write(ffwd_stream, string(wick_length)); write(ffwd_stream, ", ");
	write(ffwd_stream, string(dt));          write(ffwd_stream, ", ");
	write(ffwd_stream, string(wick_weight)); write(ffwd_stream, "\n");
end
print("Wrote $trials trials of `Fast Forward`.\n\n");

# Close file I/O
close(target_stream);
close(ffwd_stream);
print("Simulation complete. Results written to:\n");
print("+ Target model: $target_filename\n");
print("+ Fast forward: $ffwd_filename\n");


text:
function Mean()
    return Mean(())
end

function Mean(arg0::FirstMoment)
    return Mean((FirstMoment,), arg0)
end

function Mean(arg0::Mean)
    return Mean((Mean,), arg0)
end

function aggregate(obj::Mean, arg0::Mean)
    return jcall(obj, "aggregate", void, (Mean,), arg0)
end

function clear(obj::Mean)
    return jcall(obj, "clear", void, ())
end

function copy(obj::Mean)
    return jcall(obj, "copy", Mean, ())
end

function evaluate(obj::Mean, arg0::Vector{jdouble}, arg1::Vector{jdouble}, arg2::jint, arg3::jint)
    return jcall(obj, "evaluate", jdouble, (Vector{jdouble}, Vector{jdouble}, jint, jint), arg0, arg1, arg2, arg3)
end

function evaluate(obj::Mean, arg0::Vector{jdouble}, arg1::jint, arg2::jint)
    return jcall(obj, "evaluate", jdouble, (Vector{jdouble}, jint, jint), arg0, arg1, arg2)
end

function evaluate(obj::WeightedEvaluation, arg0::Vector{jdouble}, arg1::Vector{jdouble})
    return jcall(obj, "evaluate", jdouble, (Vector{jdouble}, Vector{jdouble}), arg0, arg1)
end

function get_n(obj::Mean)
    return jcall(obj, "getN", jlong, ())
end

function get_result(obj::Mean)
    return jcall(obj, "getResult", jdouble, ())
end

function increment(obj::Mean, arg0::jdouble)
    return jcall(obj, "increment", void, (jdouble,), arg0)
end


text:
"""
    b = int_F²(f)

Compute the time profile contribution to the b-value. To obtain the b-value,
multiply the result by `q^2 = (γg)^2`, where `γ` is the gyromagnetic ratio of the
water proton, `g` is the gradient amplitude, and `b = q^2 * int_F²(f)`.
"""
int_F²(f::TimeProfile{T}) where {T} = quadgk(s -> integral(f, s)^2, zero(T), echotime(f))
int_F²(f::PGSE) = f.δ^2 * (f.Δ - f.δ / 3)
int_F²(f::DoublePGSE) = 2f.δ^2 * (f.Δ - f.δ / 3)
int_F²(f::CosOGSE) = f.δ^3 / (2π * f.nperiod)^2
int_F²(f::SinOGSE) = 3f.δ^3 / (2π * f.nperiod)^2


text:
# fld-read-deprecated.jl

using MIRT: fld_header, fld_read
import MIRT: arg_get, string_to_array, fld_read_single
using Test: @test, @testset, @test_throws

dir = mktempdir()

@testset "fld eof" begin
	file = joinpath(dir, "fld-read-eof-test.fld")
	open(file, "w") do fid
	#	@show fid
	    println(fid, "# AVS field file eof test")
	end
	@test_throws String fld_read(file)

	open(file, "r") do fid
		@test_throws String fld_read_single(
			file, fid,
			(99,), Float32, fieldtype,
			false, # is_external_file,
			"", # extfile,
			Float32, "ieee-le", 0, 0,
		)
	end
end

@testset "fld misc" begin
	@test_throws String arg_get(["1", "2"], "3")
	@test string_to_array("0") isa AbstractArray{<:AbstractString}
end

@testset "fld ascii" begin
	tdir = mktempdir()
	file = joinpath(tdir, "fld-read-test.fld")
	extfilename = "fld-read-test.txt"
	chat = isinteractive()

	data = Float32.(1:5)
	head = [
		"# AVS field file"
		"ndim=1"
		"dim1=5"
		"nspace=1"
		"veclen=1"
		"data=float"
		"field=uniform"
		"variable 1 file=$extfilename filetype=ascii"
	]
	open(file, "w") do fid
		for line in head
			println(fid, line)
		end
	end

	@test fld_header(file) isa Array{<:AbstractString}

	# ensure that it throws if no extfile
	@test_throws String tmp = fld_read(file)

	# write ascii and test
	extfile = joinpath(tdir, extfilename)
	open(extfile, "w") do fid
		for d in data
			println(fid, d)
		end
	end

#	run(`op range $file`)
	tmp = fld_read(file ; chat=chat)
	@test tmp == data
	@test eltype(tmp) == Float32

	rm(file)
	rm(extfile)
end

@testset "read no ext" begin # missing extfile
	tdir = mktempdir()
	file = joinpath(tdir, "fld-read-test.fld")
	head = [
		"# AVS field file"
		"ndim=1"
		"dim1=4"
		"nspace=1"
		"veclen=1"
		"data=float"
		"field=uniform"
		"variable 1 file=bad.txt filetype=binary"
	]
	open(file, "w") do fid
		for line in head
			println(fid, line)
		end
	end

	@test_throws String fld_read(file)
	rm(file)
end

# todo: test read slices when that feature is added
# fld_write ...


text:
using DelimitedFiles: readdlm
llhlist = readdlm("../data/BrownianModels/llh_phylip.txt", header=true)[1];

@testset "Reference likelihoods" begin
	for i in 1:size(llhlist, 1)
		ft, fdm, ref_llh_fixedv, ref_llh_optimv = llhlist[i, :];

		dm = read_species_data("../data/BrownianModels/" * fdm, Float64);
		tree = read_newick("../data/BrownianModels/" * ft)[1]
		tree.dir = dm.dir;
		setspecies!(tree)

		init_model!(tree, dm, true)
		llh_fixedv = calc_llh!(tree, 1)
		@test ≈(ref_llh_fixedv, phylip_llh(tree, 1), atol = 0.01)
		delete_datablock!(tree, 1)
		pluck_nonsplitting!(tree)
		init_model!(tree, dm, true)
		calc_llh!(tree, 1)
		@test ≈(llh_fixedv, calc_llh!(tree, 1))

		init_model!(tree, dm, false)
		optimv!(tree, 2, niter = 5)
		llh_optimv = calc_llh!(tree, 2)
		@test llh_optimv ≥ llh_fixedv
		if ! isnan(ref_llh_optimv)
			@test ≈(phylip_llh(tree, 2), ref_llh_optimv, rtol=1)
		end
	end
end


text:
module HydrationOfCement

include("functionsForPackage.jl") #Include all the functionality
#import functionsForPackage

export coolfunc, coolfunc2 #Export the functions you want users to use
#coolFunc(2)

end


text:
# -------------------- FCS --------------------
abstract type FCS end

# Interface that must be implemented
function set_stick_lon! end
function set_stick_lat! end
function set_pedals! end
function set_thtl! end

"""
If allow_out_of_range, values outside the range will be allowed. This is useful in some
situations such as the trimming process. Otherwise, the values will be set to the limit
if throw_error is false (otherwise an error will be thrown)
"""
get_allow_out_of_range_inputs(fcs::FCS) = fcs.allow_out_of_range
"""
If allow_out_of_range, values outside the range will be allowed. This is useful in some
situations such as the trimming process. Otherwise, the values will be set to the limit
if throw_error is false (otherwise an error will be thrown)
"""
get_throw_error_on_out_of_range_inputs(fcs::FCS) = fcs.throw_error_on_out_of_range

"""
If allow_out_of_range, values outside the range will be allowed. This is useful in some
situations such as the trimming process. Otherwise, the values will be set to the limit
if throw_error is false (otherwise an error will be thrown)
"""
function set_allow_out_of_range_inputs!(fcs::FCS, val::Bool)
    fcs.allow_out_of_range = val
end

"""
If allow_out_of_range, values outside the range will be allowed. This is useful in some
situations such as the trimming process. Otherwise, the values will be set to the limit
if throw_error is false (otherwise an error will be thrown)
"""
function set_throw_error_on_out_of_range_inputs!(fcs::FCS, val::Bool)
    fcs.throw_error_on_out_of_range = val
end

# Interface to be trimmed
function get_controls_ranges_trimmer end

function copy(fcs::FCS)
    fcs_type = typeof(fcs)
    fcs_ = fcs_type(
        [copy(getfield(fcs, name)) for name in fieldnames(fcs_type)]...
        )
    return fcs_
end

# -------------------- Controls --------------------
"Controls are the input to a FCS."
abstract type Controls end

# Interface that must be implemented
"""
    set_controls!(fcs, controls)

Sets a FCS with the given controls.
"""
function set_controls! end

function get_controls_array end
function get_stick_lon end
function get_stick_lat end
function get_pedals end
function get_thrust_lever end

get_n_controls(c::Controls) = size(get_controls_array)

"""
    StickPedalsLeverControls(stick_lon, stick_lat, pedals, thrust_lever)

Classical controls stick longitudinal and lateral setting, pedals, and thrust lever.
"""
struct StickPedalsLeverControls <: Controls
    stick_lon :: Number
    stick_lat :: Number
    pedals :: Number
    thrust_lever :: Number
end

get_controls_array(c::StickPedalsLeverControls) = [
    get_stick_lon(c), get_stick_lat(c), get_pedals(c), get_thrust_lever(c)
    ]
get_stick_lon(c::StickPedalsLeverControls) = c.stick_lon
get_stick_lat(c::StickPedalsLeverControls) = c.stick_lat
get_pedals(c::StickPedalsLeverControls) = c.pedals
get_thrust_lever(c::StickPedalsLeverControls) = c.thrust_lever

function set_controls!(fcs::FCS, c::StickPedalsLeverControls)
    set_stick_lon!(fcs, get_stick_lon(c))
    set_stick_lat!(fcs, get_stick_lat(c))
    set_pedals!(fcs, get_pedals(c))
    set_thtl!(fcs, get_thrust_lever(c))
end

# -------------------- ControlsStream --------------------
"""
Controls are defined normally in a strem (function or callable) that returns the control
value a certain time.
"""
abstract type ControlsStream end

# Interface that must be implemented
"""
    get_controls(cs::ControlsStream, t)

Get the control values of a control stream at a certain instant t. Return Controls.
"""
function get_controls end

"""
    StickPedalsLeverStream(stick_lon, stick_lat, pedals, thrust_lever)

Classical controls stick longitudinal and lateral setting, pedals, and thrust lever.
Inputs are of type Input or Array{T, 1} where T<:Input
"""
struct StickPedalsLeverStream <: ControlsStream
    stick_lon_stream :: Union{Input, Array{T, 1}} where T<:Input
    stick_lat_stream :: Union{Input, Array{T, 1}} where T<:Input
    pedals_stream :: Union{Input, Array{T, 1}} where T<:Input
    thrust_lever_stream :: Union{Input, Array{T, 1}} where T<:Input
end

function get_controls(cs::StickPedalsLeverStream, t)
    slon = get_value(cs.stick_lon_stream, t)
    slat = get_value(cs.stick_lat_stream, t)
    ped = get_value(cs.pedals_stream, t)
    lev = get_value(cs.thrust_lever_stream, t)
    return StickPedalsLeverControls(slon, slat, ped, lev)
end


text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:light
#     text_representation:
#       extension: .jl
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.4.2
#   kernelspec:
#     display_name: Julia 1.4.0
#     language: julia
#     name: julia-1.4
# ---

# # Strong Landau Damping
#
# Electrostatic example of strong Landau damping
# $$
# f(x,v) =\frac{1}{2\pi\sigma^2}  \exp
# \Big( - \frac{v_1^2 + v_2^2}{2\sigma^2} \Big)
# ( 1+\alpha \cos(k x)􏰁),
# $$
#

using ProgressMeter, Plots
using CSV, Dates, FFTW
using GEMPIC
using Sobol
using Random
using Distributions
using JLD2
using FileIO

function main_sample!( pg   :: ParticleGroup{1,1},
                  ps   :: ParticleSampler,
                  df   :: GEMPIC.CosSumGaussianSpin,
                  mesh :: Mesh )

    s = zeros( pg.n_spin )

    theta = 0.0
    phi = 0.0
    n_rnds = 0
    if df.params.n_gaussians > 1
       n_rnds = 1
    end

    δ = zeros(df.params.n_gaussians)
    for i_v=1:df.params.n_gaussians
       δ[i_v] = sum(df.params.δ[1:i_v])
    end

    n_rnds += 2
    rdn = zeros(3)

    if ps.sampling_type == :sobol
       rng_sobol  = SobolSeq(1)
    end

    rng_random = MersenneTwister(ps.seed)

    d = Normal()

    for i_part = 1:(pg.n_particles)

       if ps.sampling_type == :sobol
           x = mesh.xmin[1] + Sobol.next!(rng_sobol)[1] * mesh.Lx[1]
       else
           x = mesh.xmin[1] + rand(rng_random) * mesh.Lx[1]
       end

       v = rand(rng_random, d)

       # For multiple Gaussian, draw which one to take
       rnd_no = rdn[3]

       i_gauss = 1
       while( rnd_no > δ[i_gauss] )
          i_gauss += 1
       end
       v = v * df.params.σ[i_gauss] + df.params.μ[i_gauss]

#for a peaked initial condition in the s direction
       s = [0, 0, 1]
#for a uniformly distributed initial condition on the sphere
#       for tt = 1:10
#            s[1] = randn()
#            s[2] = randn()
#            s[3] = randn()
#            if norm(s) > 10^(-4)
#                break
#            end
#        end
#        s .= s./norm(s)

       # Set weight according to value of perturbation
       w  = GEMPIC.eval_x_density(df, x) * prod(mesh.Lx)

       # Copy the generated numbers to the particle
       GEMPIC.set_x(pg, i_part, x[1])
       GEMPIC.set_v(pg, i_part, v)
       GEMPIC.set_spin(pg, i_part, 1, s[1])
       GEMPIC.set_spin(pg, i_part, 2, s[2])
       GEMPIC.set_spin(pg, i_part, 3, s[3])
       # Set weights.
       GEMPIC.set_weights(pg, i_part, w)

    end

end

# +
function run( steps :: Int64)

   σ, μ = 0.02, 0.0
   kx, α = 1.004355, 0.001
   xmin, xmax = 0, 2π/kx
   nx = 32 #NC 1024
   n_particles = 1000 #NC 100000
   mesh = Mesh( xmin, xmax, nx)
   spline_degree = 3

   df = GEMPIC.CosSumGaussianSpin([[kx]], [α], [[σ]], [[μ]] )

   particle_group2 = ParticleGroup{1,1}( n_particles, n_spin=3)
   sampler = ParticleSampler{1,1}( :sobol, n_particles)

   main_sample!(particle_group2, sampler, df, mesh)

   particle_group = ParticleGroup{1,1}( n_particles, n_spin=3)

   for  i_part = 1:n_particles

       x = zeros( 1 )
       v = zeros( 1 )
       s = zeros( 3 )
       w = zeros( 1 )
       x = GEMPIC.get_x(particle_group2, i_part)
       v = GEMPIC.get_v(particle_group2, i_part)

       s1 = GEMPIC.get_spin(particle_group2, i_part, 1)
       s2 = GEMPIC.get_spin(particle_group2, i_part, 2)
       s3 = GEMPIC.get_spin(particle_group2, i_part, 3)

       w = GEMPIC.get_weights(particle_group2, i_part)

       GEMPIC.set_x(particle_group, i_part, x[1])
       GEMPIC.set_v(particle_group, i_part, v[1])
       GEMPIC.set_spin(particle_group, i_part, 1, s1)
       GEMPIC.set_spin(particle_group, i_part, 2, s2)
       GEMPIC.set_spin(particle_group, i_part, 3, s3)
       GEMPIC.set_weights(particle_group, i_part, w[1])

   end

   kernel_smoother2 = ParticleMeshCoupling( mesh, n_particles, spline_degree-2, :galerkin)
   kernel_smoother1 = ParticleMeshCoupling( mesh, n_particles, spline_degree-1, :galerkin)
   kernel_smoother0 = ParticleMeshCoupling( mesh, n_particles, spline_degree, :galerkin)

   rho = zeros(Float64, nx)
   efield_poisson = zeros(Float64, nx)

   # Init!ialize the field solver
   maxwell_solver = Maxwell1DFEM(mesh, spline_degree)
   # efield by Poisson
   solve_poisson!( efield_poisson, particle_group, kernel_smoother0, maxwell_solver, rho )

   # Initialize the arrays for the spline coefficients of the fields
   k0 = 12.0523
   E0 = 10
   ww = 12.104827940833333

   #Test 1 with only Ey and Ay (and with HH=0)
   Ey(x) = E0*cos(k0*x)
   Ez(x) = 0.0*E0*sin(k0*x)
   Ay(x) = -E0/ww*sin(k0*x)
   Az(x) = 0.0*E0/ww*cos(k0*x)

   efield_dofs = [efield_poisson, zeros(Float64, nx), zeros(Float64, nx)]
   afield_dofs = [zeros(Float64, nx), zeros(Float64, nx)]

   l2projection!( efield_dofs[2], maxwell_solver, Ey, spline_degree)
   l2projection!( efield_dofs[3], maxwell_solver, Ez, spline_degree)
   l2projection!( afield_dofs[1], maxwell_solver, Ay, spline_degree)
   l2projection!( afield_dofs[2], maxwell_solver, Az, spline_degree)

   propagator = HamiltonianSplittingSpin( maxwell_solver,
                                      kernel_smoother0,
                                      kernel_smoother1,
                                      kernel_smoother2,
                                      particle_group,
                                      efield_dofs,
                                      afield_dofs)

   efield_dofs_n = propagator.e_dofs

   thdiag = TimeHistoryDiagnosticsSpin( particle_group, maxwell_solver,
                           kernel_smoother0, kernel_smoother1 );

#    mode = zeros(ComplexF64,steps,nx)
     th_modes = Vector{ComplexF64}[]
    elec_tmp = zeros(Float64,nx)

   Δt = 0.002

   @showprogress 1 for jstep = 1:steps # loop over time

       # Strang splitting
       strang_splitting!(propagator, Δt, 1)

       solve_poisson!( efield_poisson, particle_group,
                       kernel_smoother0, maxwell_solver, rho)

       write_step!(thdiag, jstep * Δt, spline_degree,
                       efield_dofs,  afield_dofs,
                       efield_dofs_n, efield_poisson, propagator)

       #diagnostics

       #store particles at some specific times
       if (jstep % 1000 == 0)
           GEMPIC.save( "save_particles", jstep, particle_group)
       end

       #Fourier modes of the longitudinal electric field
       for i = 1:nx
           elec_tmp[i] = GEMPIC.evaluate(thdiag.kernel_smoother_1, (i-1)*propagator.delta_x[1], efield_dofs[1])
       end
       push!(th_modes,fft(elec_tmp))

   end

   thdiag.data, th_modes

end
# +
thdiag, th_modes = run(100) # choose number of steps

CSV.write("thdiag-$(now()).csv", thdiag)

file="th_modes"
#datafile = @sprintf("%s.jld2", file)

FileIO.save("th_modes.jld2", Dict("modes" => th_modes))

# -


text:
###############################################################################
#
#   ArbTypes.jl : Parent and object types for Arb
#
#   Copyright (C) 2015 Tommy Hofmann
#   Copyright (C) 2015 Fredrik Johansson
#
###############################################################################

export ArbField, arb
export AcbField, acb
export ArbPolyRing, arb_poly
export AcbPolyRing, acb_poly
export ArbMatSpace, arb_mat
export AcbMatSpace, acb_mat

arb_check_prec(p::Int) = (p >= 2 && p < (typemax(Int) >> 4)) || throw(ArgumentError("invalid precision"))

# Rounding modes

const ARB_RND_DOWN = Cint(0)   # towards zero
const ARB_RND_UP = Cint(1)     # away from zero
const ARB_RND_FLOOR = Cint(2)  # towards -infinity
const ARB_RND_CEIL = Cint(3)   # towards +infinity
const ARB_RND_NEAR = Cint(4)   # to nearest

################################################################################
#
#  Types and memory management for ArbField
#
################################################################################

mutable struct ArbField <: Field
  prec::Int

  function ArbField(p::Int = 256, cached::Bool = true)
    arb_check_prec(p)
    if haskey(ArbFieldID, p)
      return ArbFieldID[p]
    else
      z = new(p)
      if cached
        ArbFieldID[p] = z
      end
      return z
    end
  end
end

const ArbFieldID = Dict{Int, ArbField}()

prec(x::ArbField) = x.prec

# these may be used for shallow operations
mutable struct arf_struct
  exp::Int # fmpz
  size::UInt # mp_size_t
  d1::UInt # mantissa_struct
  d2::UInt
end

mutable struct mag_struct
  exp::Int # fmpz
  man::UInt # mp_limb_t
end

mutable struct arb_struct
  mid_exp::Int # fmpz
  mid_size::UInt # mp_size_t
  mid_d1::UInt # mantissa_struct
  mid_d2::UInt
  rad_exp::Int # fmpz
  rad_man::UInt
end

mutable struct acb_struct
  real_mid_exp::Int # fmpz
  real_mid_size::UInt # mp_size_t
  real_mid_d1::Int # mantissa_struct
  real_mid_d2::Int
  real_rad_exp::Int # fmpz
  real_rad_man::UInt
  imag_mid_exp::Int # fmpz
  imag_mid_size::UInt # mp_size_t
  imag_mid_d1::Int # mantissa_struct
  imag_mid_d2::Int
  imag_rad_exp::Int # fmpz
  imag_rad_man::UInt
end

mutable struct arb <: FieldElem
  mid_exp::Int # fmpz
  mid_size::UInt # mp_size_t
  mid_d1::UInt # mantissa_struct
  mid_d2::UInt
  rad_exp::Int # fmpz
  rad_man::UInt
  parent::ArbField

  function arb()
    z = new()
    ccall((:arb_init, :libarb), Void, (Ref{arb}, ), z)
    finalizer(z, _arb_clear_fn)
    return z
  end

  function arb(x::Union{Int, UInt, Float64, fmpz, fmpq,
                        BigFloat, AbstractString, arb}, p::Int)
    z = new()
    ccall((:arb_init, :libarb), Void, (Ref{arb}, ), z)
    _arb_set(z, x, p)
    finalizer(z, _arb_clear_fn)
    return z
  end

  function arb(x::Union{Int, UInt, Float64, fmpz, BigFloat, arb})
    z = new()
    ccall((:arb_init, :libarb), Void, (Ref{arb}, ), z)
    _arb_set(z, x)
    finalizer(z, _arb_clear_fn)
    return z
  end

  function arb(mid::arb, rad::arb)
    z = new()
    ccall((:arb_init, :libarb), Void, (Ref{arb}, ), z)
    ccall((:arb_set, :libarb), Void, (Ref{arb}, Ref{arb}), z, mid)
    ccall((:arb_add_error, :libarb), Void, (Ref{arb}, Ref{arb}), z, rad)
    finalizer(z, _arb_clear_fn)
    return z
  end

  #function arb(x::arf)
  #  z = new()
  #  ccall((:arb_init, :libarb), Void, (Ref{arb}, ), z)
  #  ccall((:arb_set_arf, :libarb), Void, (Ref{arb}, Ref{arf}), z, x)
  #  finalizer(z, _arb_clear_fn)
  #  return z
  #end
end

function _arb_clear_fn(x::arb)
  ccall((:arb_clear, :libarb), Void, (Ref{arb}, ), x)
end

################################################################################
#
#  Types and memory management for AcbField
#
################################################################################

mutable struct AcbField <: Field
  prec::Int

  function AcbField(p::Int = 256, cached::Bool = true)
    arb_check_prec(p)
    if haskey(AcbFieldID, p)
      return AcbFieldID[p]
    else
      z = new(p)
      if cached
        AcbFieldID[p] = z
      end
      return z
    end
  end
end

const AcbFieldID = Dict{Int, AcbField}()

prec(x::AcbField) = x.prec

mutable struct acb <: FieldElem
  real_mid_exp::Int     # fmpz
  real_mid_size::UInt # mp_size_t
  real_mid_d1::UInt    # mantissa_struct
  real_mid_d2::UInt
  real_rad_exp::Int     # fmpz
  real_rad_man::UInt
  imag_mid_exp::Int     # fmpz
  imag_mid_size::UInt # mp_size_t
  imag_mid_d1::UInt    # mantissa_struct
  imag_mid_d2::UInt
  imag_rad_exp::Int     # fmpz
  imag_rad_man::UInt
  parent::AcbField

  function acb()
    z = new()
    ccall((:acb_init, :libarb), Void, (Ref{acb}, ), z)
    finalizer(z, _acb_clear_fn)
    return z
  end

  function acb(x::Union{Int, UInt, Float64, fmpz, BigFloat, arb, acb})
    z = new()
    ccall((:acb_init, :libarb), Void, (Ref{acb}, ), z)
    _acb_set(z, x)
    finalizer(z, _acb_clear_fn)
    return z
  end

  function acb(x::Union{Int, UInt, Float64, fmpz, fmpq,
                        BigFloat, arb, acb, AbstractString}, p::Int)
    z = new()
    ccall((:acb_init, :libarb), Void, (Ref{acb}, ), z)
    _acb_set(z, x, p)
    finalizer(z, _acb_clear_fn)
    return z
  end

  #function acb{T <: Union{Int, UInt, Float64, fmpz, BigFloat, arb}}(x::T, y::T)
  #  z = new()
  #  ccall((:acb_init, :libarb), Void, (Ref{acb}, ), z)
  #  _acb_set(z, x, y)
  #  finalizer(z, _acb_clear_fn)
  #  return z
  #end

  function acb(x::T, y::T, p::Int) where {T <: Union{Int, UInt, Float64, fmpz, fmpq, BigFloat, AbstractString, arb}}
    z = new()
    ccall((:acb_init, :libarb), Void, (Ref{acb}, ), z)
    _acb_set(z, x, y, p)
    finalizer(z, _acb_clear_fn)
    return z
  end
end

function _acb_clear_fn(x::acb)
  ccall((:acb_clear, :libarb), Void, (Ref{acb}, ), x)
end

################################################################################
#
#  Types and memory management for ArbPolyRing
#
################################################################################

mutable struct ArbPolyRing <: PolyRing{arb}
  base_ring::ArbField
  S::Symbol

  function ArbPolyRing(R::ArbField, S::Symbol, cached::Bool = true)
    if haskey(ArbPolyRingID, (R, S))
      return ArbPolyRingID[R, S]
    else
      z = new(R, S)
      if cached
        ArbPolyRingID[R, S] = z
      end
      return z
    end
  end
end

const ArbPolyRingID = Dict{Tuple{ArbField, Symbol}, ArbPolyRing}()

mutable struct arb_poly <: PolyElem{arb}
  coeffs::Ref{Void}
  length::Int
  alloc::Int
  parent::ArbPolyRing

  function arb_poly()
    z = new()
    ccall((:arb_poly_init, :libarb), Void, (Ref{arb_poly}, ), z)
    finalizer(z, _arb_poly_clear_fn)
    return z
  end

  function arb_poly(x::arb, p::Int)
    z = new()
    ccall((:arb_poly_init, :libarb), Void, (Ref{arb_poly}, ), z)
    ccall((:arb_poly_set_coeff_arb, :libarb), Void,
                (Ref{arb_poly}, Int, Ref{arb}), z, 0, x)
    finalizer(z, _arb_poly_clear_fn)
    return z
  end

  function arb_poly(x::Array{arb, 1}, p::Int)
    z = new()
    ccall((:arb_poly_init, :libarb), Void, (Ref{arb_poly}, ), z)
    for i = 1:length(x)
        ccall((:arb_poly_set_coeff_arb, :libarb), Void,
                (Ref{arb_poly}, Int, Ref{arb}), z, i - 1, x[i])
    end
    finalizer(z, _arb_poly_clear_fn)
    return z
  end

  function arb_poly(x::arb_poly)
    z = new()
    ccall((:arb_poly_init, :libarb), Void, (Ref{arb_poly}, ), z)
    ccall((:arb_poly_set, :libarb), Void, (Ref{arb_poly}, Ref{arb_poly}), z, x)
    finalizer(z, _arb_poly_clear_fn)
    return z
  end

  function arb_poly(x::arb_poly, p::Int)
    z = new()
    ccall((:arb_poly_init, :libarb), Void, (Ref{arb_poly}, ), z)
    ccall((:arb_poly_set_round, :libarb), Void,
                (Ref{arb_poly}, Ref{arb_poly}, Int), z, x, p)
    finalizer(z, _arb_poly_clear_fn)
    return z
  end

  function arb_poly(x::fmpz_poly, p::Int)
    z = new()
    ccall((:arb_poly_init, :libarb), Void, (Ref{arb_poly}, ), z)
    ccall((:arb_poly_set_fmpz_poly, :libarb), Void,
                (Ref{arb_poly}, Ref{fmpz_poly}, Int), z, x, p)
    finalizer(z, _arb_poly_clear_fn)
    return z
  end

  function arb_poly(x::fmpq_poly, p::Int)
    z = new()
    ccall((:arb_poly_init, :libarb), Void, (Ref{arb_poly}, ), z)
    ccall((:arb_poly_set_fmpq_poly, :libarb), Void,
                (Ref{arb_poly}, Ref{fmpq_poly}, Int), z, x, p)
    finalizer(z, _arb_poly_clear_fn)
    return z
  end
end

function _arb_poly_clear_fn(x::arb_poly)
  ccall((:arb_poly_clear, :libarb), Void, (Ref{arb_poly}, ), x)
end

parent(x::arb_poly) = x.parent

elem_type(x::ArbPolyRing) = arb_poly

var(x::ArbPolyRing) = x.S

prec(x::ArbPolyRing) = prec(x.base_ring)

base_ring(a::ArbPolyRing) = a.base_ring

################################################################################
#
#  Types and memory management for AcbPolyRing
#
################################################################################

mutable struct AcbPolyRing <: PolyRing{acb}
  base_ring::AcbField
  S::Symbol

  function AcbPolyRing(R::AcbField, S::Symbol, cached::Bool = true)
    if haskey(AcbPolyRingID, (R, S))
      return AcbPolyRingID[R, S]
    else
      z = new(R, S)
      if cached
        AcbPolyRingID[R, S] = z
      end
      return z
    end
  end
end

const AcbPolyRingID = Dict{Tuple{AcbField, Symbol}, AcbPolyRing}()

mutable struct acb_poly <: PolyElem{acb}
  coeffs::Ref{Void}
  length::Int
  alloc::Int
  parent::AcbPolyRing

  function acb_poly()
    z = new()
    ccall((:acb_poly_init, :libarb), Void, (Ref{acb_poly}, ), z)
    finalizer(z, _acb_poly_clear_fn)
    return z
  end

  function acb_poly(x::acb, p::Int)
    z = new()
    ccall((:acb_poly_init, :libarb), Void, (Ref{acb_poly}, ), z)
    ccall((:acb_poly_set_coeff_acb, :libarb), Void,
                (Ref{acb_poly}, Int, Ref{acb}), z, 0, x)
    finalizer(z, _acb_poly_clear_fn)
    return z
  end

  function acb_poly(x::Array{acb, 1}, p::Int)
    z = new()
    ccall((:acb_poly_init, :libarb), Void, (Ref{acb_poly}, ), z)
    for i = 1:length(x)
        ccall((:acb_poly_set_coeff_acb, :libarb), Void,
                (Ref{acb_poly}, Int, Ref{acb}), z, i - 1, x[i])
    end
    finalizer(z, _acb_poly_clear_fn)
    return z
  end

  function acb_poly(x::acb_poly)
    z = new()
    ccall((:acb_poly_init, :libarb), Void, (Ref{acb_poly}, ), z)
    ccall((:acb_poly_set, :libarb), Void, (Ref{acb_poly}, Ref{acb_poly}), z, x)
    finalizer(z, _acb_poly_clear_fn)
    return z
  end

  function acb_poly(x::arb_poly, p::Int)
    z = new()
    ccall((:acb_poly_init, :libarb), Void, (Ref{acb_poly}, ), z)
    ccall((:acb_poly_set_arb_poly, :libarb), Void,
                (Ref{acb_poly}, Ref{arb_poly}, Int), z, x, p)
    ccall((:acb_poly_set_round, :libarb), Void,
                (Ref{acb_poly}, Ref{acb_poly}, Int), z, z, p)
    finalizer(z, _acb_poly_clear_fn)
    return z
  end

  function acb_poly(x::acb_poly, p::Int)
    z = new()
    ccall((:acb_poly_init, :libarb), Void, (Ref{acb_poly}, ), z)
    ccall((:acb_poly_set_round, :libarb), Void,
                (Ref{acb_poly}, Ref{acb_poly}, Int), z, x, p)
    finalizer(z, _acb_poly_clear_fn)
    return z
  end

  function acb_poly(x::fmpz_poly, p::Int)
    z = new()
    ccall((:acb_poly_init, :libarb), Void, (Ref{acb_poly}, ), z)
    ccall((:acb_poly_set_fmpz_poly, :libarb), Void,
                (Ref{acb_poly}, Ref{fmpz_poly}, Int), z, x, p)
    finalizer(z, _acb_poly_clear_fn)
    return z
  end

  function acb_poly(x::fmpq_poly, p::Int)
    z = new()
    ccall((:acb_poly_init, :libarb), Void, (Ref{acb_poly}, ), z)
    ccall((:acb_poly_set_fmpq_poly, :libarb), Void,
                (Ref{acb_poly}, Ref{fmpq_poly}, Int), z, x, p)
    finalizer(z, _acb_poly_clear_fn)
    return z
  end
end

function _acb_poly_clear_fn(x::acb_poly)
  ccall((:acb_poly_clear, :libarb), Void, (Ref{acb_poly}, ), x)
end

parent(x::acb_poly) = x.parent

elem_type(x::AcbPolyRing) = acb_poly

var(x::AcbPolyRing) = x.S

prec(x::AcbPolyRing) = prec(x.base_ring)

base_ring(a::AcbPolyRing) = a.base_ring

################################################################################
#
#  Types and memory management for ArbMatSpace
#
################################################################################

mutable struct ArbMatSpace <: MatSpace{arb}
  rows::Int
  cols::Int
  base_ring::ArbField

  function ArbMatSpace(R::ArbField, r::Int, c::Int, cached::Bool = true)
    if haskey(ArbMatSpaceID, (R, r, c))
      return ArbMatSpaceID[(R, r, c)]
    else
      z = new(r, c, R)
      if cached
        ArbMatSpaceID[(R, r, c)] = z
      end
      return z::ArbMatSpace
    end
  end
end

const ArbMatSpaceID = Dict{Tuple{ArbField, Int, Int}, ArbMatSpace}()

mutable struct arb_mat <: MatElem{arb}
  entries::Ref{Void}
  r::Int
  c::Int
  rows::Ref{Void}
  base_ring::ArbField

  function arb_mat(r::Int, c::Int)
    z = new()
    ccall((:arb_mat_init, :libarb), Void, (Ref{arb_mat}, Int, Int), z, r, c)
    finalizer(z, _arb_mat_clear_fn)
    return z
  end

  function arb_mat(a::fmpz_mat)
    z = new()
    ccall((:arb_mat_init, :libarb), Void,
                (Ref{arb_mat}, Int, Int), z, a.r, a.c)
    ccall((:arb_mat_set_fmpz_mat, :libarb), Void,
                (Ref{arb_mat}, Ref{fmpz_mat}), z, a)
    finalizer(z, _arb_mat_clear_fn)
    return z
  end

  function arb_mat(a::fmpz_mat, prec::Int)
    z = new()
    ccall((:arb_mat_init, :libarb), Void,
                (Ref{arb_mat}, Int, Int), z, a.r, a.c)
    ccall((:arb_mat_set_round_fmpz_mat, :libarb), Void,
                (Ref{arb_mat}, Ref{fmpz_mat}, Int), z, a, prec)
    finalizer(z, _arb_mat_clear_fn)
    return z
  end

  function arb_mat(r::Int, c::Int, arr::Array{T, 2}) where {T <: Union{Int, UInt, fmpz, Float64, BigFloat, arb}}
    z = new()
    ccall((:arb_mat_init, :libarb), Void,
                (Ref{arb_mat}, Int, Int), z, r, c)
    finalizer(z, _arb_mat_clear_fn)
    for i = 1:r
      for j = 1:c
        el = ccall((:arb_mat_entry_ptr, :libarb), Ref{arb},
                    (Ref{arb_mat}, Int, Int), z, i - 1, j - 1)
        Nemo._arb_set(el, arr[i, j])
      end
    end
    return z
  end

  function arb_mat(r::Int, c::Int, arr::Array{T, 1}) where {T <: Union{Int, UInt, fmpz, Float64, BigFloat, arb}}
    z = new()
    ccall((:arb_mat_init, :libarb), Void,
                (Ref{arb_mat}, Int, Int), z, r, c)
    finalizer(z, _arb_mat_clear_fn)
    for i = 1:r
      for j = 1:c
        el = ccall((:arb_mat_entry_ptr, :libarb), Ref{arb},
                    (Ref{arb_mat}, Int, Int), z, i - 1, j - 1)
        Nemo._arb_set(el, arr[(i-1)*c+j])
      end
    end
    return z
  end

  function arb_mat(r::Int, c::Int, arr::Array{T, 2}, prec::Int) where {T <: Union{Int, UInt, fmpz, fmpq, Float64, BigFloat, arb, AbstractString}}
    z = new()
    ccall((:arb_mat_init, :libarb), Void,
                (Ref{arb_mat}, Int, Int), z, r, c)
    finalizer(z, _arb_mat_clear_fn)
    for i = 1:r
      for j = 1:c
        el = ccall((:arb_mat_entry_ptr, :libarb), Ref{arb},
                    (Ref{arb_mat}, Int, Int), z, i - 1, j - 1)
        _arb_set(el, arr[i, j], prec)
      end
    end
    return z
  end

  function arb_mat(r::Int, c::Int, arr::Array{T, 1}, prec::Int) where {T <: Union{Int, UInt, fmpz, fmpq, Float64, BigFloat, arb, AbstractString}}
    z = new()
    ccall((:arb_mat_init, :libarb), Void,
                (Ref{arb_mat}, Int, Int), z, r, c)
    finalizer(z, _arb_mat_clear_fn)
    for i = 1:r
      for j = 1:c
        el = ccall((:arb_mat_entry_ptr, :libarb), Ref{arb},
                    (Ref{arb_mat}, Int, Int), z, i - 1, j - 1)
        _arb_set(el, arr[(i-1)*c+j], prec)
      end
    end
    return z
  end

  function arb_mat(a::fmpq_mat, prec::Int)
    z = new()
    ccall((:arb_mat_init, :libarb), Void,
                (Ref{arb_mat}, Int, Int), z, a.r, a.c)
    ccall((:arb_mat_set_fmpq_mat, :libarb), Void,
                (Ref{arb_mat}, Ref{fmpq_mat}, Int), z, a, prec)
    finalizer(z, _arb_mat_clear_fn)
    return z
  end
end

function _arb_mat_clear_fn(x::arb_mat)
  ccall((:arb_mat_clear, :libarb), Void, (Ref{arb_mat}, ), x)
end

################################################################################
#
#  Types and memory management for AcbMatSpace
#
################################################################################

mutable struct AcbMatSpace <: MatSpace{acb}
  rows::Int
  cols::Int
  base_ring::AcbField

  function AcbMatSpace(R::AcbField, r::Int, c::Int, cached::Bool = true)
    if haskey(AcbMatSpaceID, (R, r, c))
      return AcbMatSpaceID[(R, r, c)]
    else
      z = new(r, c, R)
      if cached
        AcbMatSpaceID[(R, r, c)] = z
      end
      return z::AcbMatSpace
    end
  end
end

const AcbMatSpaceID = Dict{Tuple{AcbField, Int, Int}, AcbMatSpace}()

mutable struct acb_mat <: MatElem{acb}
  entries::Ref{Void}
  r::Int
  c::Int
  rows::Ref{Void}
  base_ring::AcbField

  function acb_mat(r::Int, c::Int)
    z = new()
    ccall((:acb_mat_init, :libarb), Void, (Ref{acb_mat}, Int, Int), z, r, c)
    finalizer(z, _acb_mat_clear_fn)
    return z
  end

  function acb_mat(a::fmpz_mat)
    z = new()
    ccall((:acb_mat_init, :libarb), Void,
                (Ref{acb_mat}, Int, Int), z, a.r, a.c)
    ccall((:acb_mat_set_fmpz_mat, :libarb), Void,
                (Ref{acb_mat}, Ref{fmpz_mat}), z, a)
    finalizer(z, _acb_mat_clear_fn)
    return z
  end

  function acb_mat(a::fmpz_mat, prec::Int)
    z = new()
    ccall((:acb_mat_init, :libarb), Void,
                (Ref{acb_mat}, Int, Int), z, a.r, a.c)
    ccall((:acb_mat_set_round_fmpz_mat, :libarb), Void,
                (Ref{acb_mat}, Ref{fmpz_mat}, Int), z, a, prec)
    finalizer(z, _acb_mat_clear_fn)
    return z
  end

  function acb_mat(a::arb_mat)
    z = new()
    ccall((:acb_mat_init, :libarb), Void,
                (Ref{acb_mat}, Int, Int), z, a.r, a.c)
    ccall((:acb_mat_set_arb_mat, :libarb), Void,
                (Ref{acb_mat}, Ref{arb_mat}), z, a)
    finalizer(z, _acb_mat_clear_fn)
    return z
  end

  function acb_mat(a::arb_mat, prec::Int)
    z = new()
    ccall((:acb_mat_init, :libarb), Void,
                (Ref{acb_mat}, Int, Int), z, a.r, a.c)
    ccall((:acb_mat_set_round_arb_mat, :libarb), Void,
                (Ref{acb_mat}, Ref{arb_mat}, Int), z, a, prec)
    finalizer(z, _acb_mat_clear_fn)
    return z
  end

  function acb_mat(r::Int, c::Int, arr::Array{T, 2}) where {T <: Union{Int, UInt, Float64, fmpz}}
    z = new()
    ccall((:acb_mat_init, :libarb), Void,
                (Ref{acb_mat}, Int, Int), z, r, c)
    finalizer(z, _acb_mat_clear_fn)
    for i = 1:r
      for j = 1:c
        el = ccall((:acb_mat_entry_ptr, :libarb), Ref{acb},
                    (Ref{acb_mat}, Int, Int), z, i - 1, j - 1)
        _acb_set(el, arr[i, j])
      end
    end
    return z
  end

  function acb_mat(r::Int, c::Int, arr::Array{T, 2}) where {T <: Union{BigFloat, acb, arb}}
    z = new()
    ccall((:acb_mat_init, :libarb), Void,
                (Ref{acb_mat}, Int, Int), z, r, c)
    finalizer(z, _acb_mat_clear_fn)
    for i = 1:r
      for j = 1:c
        el = ccall((:acb_mat_entry_ptr, :libarb), Ref{acb},
                    (Ref{acb_mat}, Int, Int), z, i - 1, j - 1)
        _acb_set(el, arr[i, j])
      end
    end
    return z
  end

  function acb_mat(r::Int, c::Int, arr::Array{T, 1}) where {T <: Union{Int, UInt, Float64, fmpz}}
    z = new()
    ccall((:acb_mat_init, :libarb), Void,
                (Ref{acb_mat}, Int, Int), z, r, c)
    finalizer(z, _acb_mat_clear_fn)
    for i = 1:r
      for j = 1:c
        el = ccall((:acb_mat_entry_ptr, :libarb), Ref{acb},
                    (Ref{acb_mat}, Int, Int), z, i - 1, j - 1)
        _acb_set(el, arr[(i-1)*c+j])
      end
    end
    return z
  end

  function acb_mat(r::Int, c::Int, arr::Array{T, 1}) where {T <: Union{BigFloat, acb, arb}}
    z = new()
    ccall((:acb_mat_init, :libarb), Void,
                (Ref{acb_mat}, Int, Int), z, r, c)
    finalizer(z, _acb_mat_clear_fn)
    for i = 1:r
      for j = 1:c
        el = ccall((:acb_mat_entry_ptr, :libarb), Ref{acb},
                    (Ref{acb_mat}, Int, Int), z, i - 1, j - 1)
        _acb_set(el, arr[(i-1)*c+j])
      end
    end
    return z
  end

  function acb_mat(r::Int, c::Int, arr::Array{T, 2}, prec::Int) where {T <: Union{Int, UInt, fmpz, fmpq, Float64}}
    z = new()
    ccall((:acb_mat_init, :libarb), Void,
                (Ref{acb_mat}, Int, Int), z, r, c)
    finalizer(z, _acb_mat_clear_fn)
    for i = 1:r
      for j = 1:c
        el = ccall((:acb_mat_entry_ptr, :libarb), Ref{acb},
                    (Ref{acb_mat}, Int, Int), z, i - 1, j - 1)
        _acb_set(el, arr[i, j], prec)
      end
    end
    return z
  end

  function acb_mat(r::Int, c::Int, arr::Array{T, 2}, prec::Int) where {T <: Union{BigFloat, arb, AbstractString, acb}}
    z = new()
    ccall((:acb_mat_init, :libarb), Void,
                (Ref{acb_mat}, Int, Int), z, r, c)
    finalizer(z, _acb_mat_clear_fn)
    for i = 1:r
      for j = 1:c
        el = ccall((:acb_mat_entry_ptr, :libarb), Ref{acb},
                    (Ref{acb_mat}, Int, Int), z, i - 1, j - 1)
        _acb_set(el, arr[i, j], prec)
      end
    end
    return z
  end

  function acb_mat(r::Int, c::Int, arr::Array{T, 1}, prec::Int) where {T <: Union{Int, UInt, fmpz, fmpq, Float64}}
    z = new()
    ccall((:acb_mat_init, :libarb), Void,
                (Ref{acb_mat}, Int, Int), z, r, c)
    finalizer(z, _acb_mat_clear_fn)
    for i = 1:r
      for j = 1:c
        el = ccall((:acb_mat_entry_ptr, :libarb), Ref{acb},
                    (Ref{acb_mat}, Int, Int), z, i - 1, j - 1)
        _acb_set(el, arr[(i-1)*c+j], prec)
      end
    end
    return z
  end

  function acb_mat(r::Int, c::Int, arr::Array{T, 1}, prec::Int) where {T <: Union{BigFloat, arb, AbstractString, acb}}
    z = new()
    ccall((:acb_mat_init, :libarb), Void,
                (Ref{acb_mat}, Int, Int), z, r, c)
    finalizer(z, _acb_mat_clear_fn)
    for i = 1:r
      for j = 1:c
        el = ccall((:acb_mat_entry_ptr, :libarb), Ref{acb},
                    (Ref{acb_mat}, Int, Int), z, i - 1, j - 1)
        _acb_set(el, arr[(i-1)*c+j], prec)
      end
    end
    return z
  end

  function acb_mat(r::Int, c::Int, arr::Array{Tuple{T, T}, 2}, prec::Int) where {T <: Union{Int, UInt, Float64, fmpz}}

    z = new()
    ccall((:acb_mat_init, :libarb), Void,
                (Ref{acb_mat}, Int, Int), z, r, c)
    finalizer(z, _acb_mat_clear_fn)
    for i = 1:r
      for j = 1:c
        el = ccall((:acb_mat_entry_ptr, :libarb), Ref{acb},
                    (Ref{acb_mat}, Int, Int), z, i - 1, j - 1)
        _acb_set(el, arr[i, j][1], arr[i,j][2], prec)
      end
    end
    return z
  end

  function acb_mat(r::Int, c::Int, arr::Array{Tuple{T, T}, 2}, prec::Int) where {T <: Union{fmpq, BigFloat, arb, AbstractString}}

    z = new()
    ccall((:acb_mat_init, :libarb), Void,
                (Ref{acb_mat}, Int, Int), z, r, c)
    finalizer(z, _acb_mat_clear_fn)
    for i = 1:r
      for j = 1:c
        el = ccall((:acb_mat_entry_ptr, :libarb), Ref{acb},
                    (Ref{acb_mat}, Int, Int), z, i - 1, j - 1)
        _acb_set(el, arr[i, j][1], arr[i,j][2], prec)
      end
    end
    return z
  end

  function acb_mat(r::Int, c::Int, arr::Array{Tuple{T, T}, 1}, prec::Int) where {T <: Union{Int, UInt, Float64, fmpz}}

    z = new()
    ccall((:acb_mat_init, :libarb), Void,
                (Ref{acb_mat}, Int, Int), z, r, c)
    finalizer(z, _acb_mat_clear_fn)
    for i = 1:r
      for j = 1:c
        el = ccall((:acb_mat_entry_ptr, :libarb), Ref{acb},
                    (Ref{acb_mat}, Int, Int), z, i - 1, j - 1)
        _acb_set(el, arr[(i-1)*c+j][1], arr[(i-1)*c+j][2], prec)
      end
    end
    return z
  end

  function acb_mat(r::Int, c::Int, arr::Array{Tuple{T, T}, 1}, prec::Int) where {T <: Union{fmpq, BigFloat, arb, AbstractString}}

    z = new()
    ccall((:acb_mat_init, :libarb), Void,
                (Ref{acb_mat}, Int, Int), z, r, c)
    finalizer(z, _acb_mat_clear_fn)
    for i = 1:r
      for j = 1:c
        el = ccall((:acb_mat_entry_ptr, :libarb), Ref{acb},
                    (Ref{acb_mat}, Int, Int), z, i - 1, j - 1)
        _acb_set(el, arr[(i-1)*c+j][1], arr[(i-1)*c+j][2], prec)
      end
    end
    return z
  end

  function acb_mat(a::fmpq_mat, prec::Int)
    z = new()
    ccall((:acb_mat_init, :libarb), Void,
                (Ref{acb_mat}, Int, Int), z, a.r, a.c)
    ccall((:acb_mat_set_fmpq_mat, :libarb), Void,
                (Ref{acb_mat}, Ref{fmpq_mat}, Int), z, a, prec)
    finalizer(z, _acb_mat_clear_fn)
    return z
  end
end

function _acb_mat_clear_fn(x::acb_mat)
  ccall((:acb_mat_clear, :libarb), Void, (Ref{acb_mat}, ), x)
end


text:
using Revise
using FastGroupBy
const M=1000; const K=100
srand(1)
# @time svec1 = rand(["i"*dec(k,rand(1:7)) for k in 1:M÷K], M)
@time svec1 = rand([string(rand(Char.(32:126), rand(1:8))...) for k in 1:M÷K], M)
@time radixsort!(svec1)
issorted(svec1)

# include("src/sorttwo_lsd16.jl")
const M=10_000_000; const K=100
srand(1)
# @time svec1 = rand(["i"*dec(k,7) for k in 1:M÷K], M)
@time svec1 = rand([string(rand(Char.(32:126), rand(1:8))...) for k in 1:M÷K], M)
# @time fast(svec1)
@time radixsort!(svec1) # 3 seconds on 10m; 41 seconds on 100m;
# @code_warntype radixsort!(svec1)
issorted(svec1)

const M=10_000_000; const K=100
srand(1)
# @time svec1 = rand(["i"*dec(k,7) for k in 1:M÷K], M)
@time svec1 = rand([string(rand(Char.(32:126), rand(1:16))...) for k in 1:M÷K], M)
# @time fast(svec1)
@time radixsort!(svec1) # 4 seconds on 10m; 41 seconds on 100m;
# @code_warntype radixsort!(svec1)
issorted(svec1)

const M=100_000_000; const K=100;
srand(1);
# @time svec1 = rand(["i"*dec(k,7) for k in 1:M÷K], M)
svec1 = rand([string(rand(Char.(32:126), rand(1:8))...) for k in 1:M÷K], M);
# @time fast(svec1)
@time radixsort!(svec1) # 19
# @code_warntype radixsort!(svec1)
issorted(svec1)

const M=100_000_000; const K=100;
srand(1);
# @time svec1 = rand(["i"*dec(k,7) for k in 1:M÷K], M)
svec1 = rand([string(rand(Char.(32:126), rand(1:8))...) for k in 1:M÷K], M);
# @time fast(svec1)
@time FastGroupBy.radixsort8!(svec1) # 17 seconds
# @code_warntype radixsort!(svec1)
issorted(svec1)

const M=200_000_000; const K=100;
srand(1);
# @time svec1 = rand(["i"*dec(k,7) for k in 1:M÷K], M)
svec1 = rand([string(rand(Char.(32:126), rand(1:8))...) for k in 1:M÷K], M);
# @time fast(svec1)
@time radixsort!(svec1) # 41 seconds
# @code_warntype radixsort!(svec1)
issorted(svec1)

const M=200_000_000; const K=100;
srand(1);
# @time svec1 = rand(["i"*dec(k,7) for k in 1:M÷K], M)
svec1 = rand([string(rand(Char.(32:126), rand(1:8))...) for k in 1:M÷K], M);
# @time fast(svec1)
@time FastGroupBy.radixsort8!(svec1) # 25-29-26 seconds slightly slower
# @code_warntype radixsort!(svec1)
issorted(svec1)

# 7seconds on 10million; 2mins on 100m and only 12 seconds in
# srand(1)
# @time svec1 = rand([string(rand(Char.(32:126), rand(1:8))...) for k in 1:M÷K], M)
# @time sort(svec1)

# srand(1)
# @time svec1 = rand(["i"*dec(k,rand(1:7)) for k in 1:M÷K], M)
# @time radixsort!(svec1)
# issorted(length.(svec1))
# issorted(svec1)

# srand(1)
# @time svec1 = rand(["i"*dec(k,rand(1:7)) for k in 1:M÷K], M)
# @time sort(svec1) # takes about 2 mins for 100m

srand(1)
# @time svec1 = rand(["id"*dec(k,rand(1:14)) for k in 1:M÷K], M)
@time svec1 = rand([string(rand(Char.(32:126), rand(1:16))...) for k in 1:M÷K], M)
@time radixsort!(svec1)
issorted(svec1)
#
# srand(1)
# @time svec1 = rand(["id"*dec(k,22) for k in 1:M÷K], M)
# @time radixsort!(svec1)
# issorted(svec1)
#
# x = svec1[1:end-1] .> svec1[2:end]
# x = vcat(false, x) .| vcat(x, false)
# (1:length(x))[x]
#
# svec1[9999248:9999249]
#
#
#
# @time svec1 = rand(["id"*dec(k,rand(1:24)) for k in 1:M÷K], M)
# #@time svec1 = ["id"*dec(k,rand(1:22)) for k in 1:M÷K]
# @time sort!(svec1)
# issorted(svec1)

# @time cc = load_bits.(svec1)
# @time sort(cc)
# include("sorttwo2.jl")
# @time sorttwo2!(cc, svec1)

# include("sorttwo2.jl")
const M=100_000_000; const K=100
srand(1)
# @time svec1 = rand(["i"*dec(k,rand(1:7)) for k in 1:M÷K], M)
@time svec1 = rand([string(rand(Char.(32:126), rand(1:8))...) for k in 1:M÷K], M)
@time lens = maximum(sizeof.(svec1)) # 3-5 econds
iters = ceil(lens/8)
@time cc = load_bits.(svec1) # 9 seconds
ii = [1:length(cc)...]
# @time sorttwo2!(cc, svec1) # 11 seconds
@time sorttwo2!(cc, ii) ## 5 seconds

@time svec1 = svec1[ii]
issorted(svec1)

# return svec
# sort the longest vector using LSD

# skipbytes = 0
# bitsrep = load_bits.(svec)
# sorttwo2!(bitsrep, svec)
# i = 0
# while !issorted(svec)
#     i = i + 1
#     skipbytes += 8
#     if i == 4
#         throw(ErrorException("wassup"))
#     end
#     sorttwo2!(load_bits.(svec, skipbytes), svec)
# end
# return svec


text:
using Test, LinearAlgebra
using DiffEqSensitivity, StochasticDiffEq
using Random

@info "Test SDE Transformation"

seed = 100
tspan = (0.0, 0.1)
p = [1.01,0.87]

# scalar
f(u,p,t) = p[1]*u
σ(u,p,t) = p[2]*u

Random.seed!(seed)
u0 = rand(1)
linear_analytic(u0,p,t,W) = @.(u0*exp((p[1]-p[2]^2/2)*t+p[2]*W))

prob = SDEProblem(SDEFunction(f,σ,analytic=linear_analytic),σ,u0,tspan,p)
sol = solve(prob,SOSRI(),adaptive=false, dt=0.001, save_noise=true)

@test isapprox(sol.u_analytic,sol.u, atol=1e-4)

du = zeros(size(u0))
u = sol.u[end]
transformed_function = StochasticTransformedFunction(sol,sol.prob.f,sol.prob.g)
#transformed_function(du,u,p,tspan[2])
du2 = transformed_function(u,p,tspan[2])

#@test du[1] == (p[1]*u[1]-p[2]^2*u[1])
@test du2[1] == (p[1]*u[1]-p[2]^2*u[1])
#@test du2 == du
transformed_function = StochasticTransformedFunction(sol,sol.prob.f,sol.prob.g,(u,p,t)->p[2]^2*u)
du2 = transformed_function(u,p,tspan[2])
@test du2[1] == (p[1]*u[1]-p[2]^2*u[1])

linear_analytic_strat(u0,p,t,W) = @.(u0*exp((p[1])*t+p[2]*W))

prob_strat = SDEProblem{false}(SDEFunction((u,p,t)->p[1]*u-1//2*p[2]^2*u,σ,analytic=linear_analytic_strat),σ,u0,tspan,p)
Random.seed!(seed)
sol_strat = solve(prob_strat,RKMil(interpretation=:Stratonovich),adaptive=false, dt=0.0001, save_noise=true)
prob_strat1 = SDEProblem{false}(SDEFunction((u,p,t)->transformed_function(u,p,t).+1//2*p[2]^2*u[1],σ,analytic=linear_analytic),σ,u0,tspan,p)
Random.seed!(seed)
sol_strat1 = solve(prob_strat1,RKMil(interpretation=:Stratonovich),adaptive=false, dt=0.0001, save_noise=true)

# Test if we recover Ito solution in Stratonovich sense
@test isapprox(sol_strat.u, sol_strat1.u, atol=1e-4) # own transformation and custom function agree
@test !isapprox(sol_strat.u_analytic,sol_strat.u, atol=1e-4) # we don't get the stratonovich solution for the linear SDE
@test isapprox(sol_strat1.u_analytic,sol_strat.u, atol=1e-4) # we do recover the analytic solution from the Ito sense

# inplace

f!(du,u,p,t) = @.(du = p[1]*u)
σ!(du,u,p,t) = @.(du = p[2]*u)

prob = SDEProblem(SDEFunction(f!,σ!,analytic=linear_analytic),σ!,u0,tspan,p)
sol = solve(prob,SOSRI(),adaptive=false, dt=0.001, save_noise=true)

@test isapprox(sol.u_analytic,sol.u, atol=1e-4)

du = zeros(size(u0))
u = sol.u[end]
transformed_function = StochasticTransformedFunction(sol,sol.prob.f,sol.prob.g)
transformed_function(du,u,p,tspan[2])

@test isapprox(du[1], (p[1]*u[1]-p[2]^2*u[1]), atol=1e-15)
# @test isapprox(du2[1], (p[1]*u[1]-p[2]^2*u[1]), atol=1e-15)
# @test isapprox(du2, du,  atol=1e-15)

transformed_function = StochasticTransformedFunction(sol,sol.prob.f,sol.prob.g,(du,u,p,t)-> (du.=p[2]^2*u))
transformed_function(du,u,p,tspan[2])
@test du[1] == (p[1]*u[1]-p[2]^2*u[1])

# diagonal noise

u0 = rand(3)

prob = SDEProblem(SDEFunction(f,σ,analytic=linear_analytic),σ,u0,tspan,p)
sol = solve(prob,SOSRI(),adaptive=false, dt=0.001, save_noise=true)
u = sol.u[end]

transformed_function = StochasticTransformedFunction(sol,sol.prob.f,sol.prob.g)
du2 = transformed_function(u,p,tspan[2])
@test isapprox(du2,(p[1]*u-p[2]^2*u), atol=1e-15)

transformed_function = StochasticTransformedFunction(sol,sol.prob.f,sol.prob.g,(u,p,t)->p[2]^2*u)
du2 = transformed_function(u,p,tspan[2])
@test du2[1] == (p[1]*u[1]-p[2]^2*u[1])

prob = SDEProblem(SDEFunction(f!,σ!,analytic=linear_analytic),σ!,u0,tspan,p)
sol = solve(prob,SOSRI(),adaptive=false, dt=0.001, save_noise=true)

du = zeros(size(u0))
u = sol.u[end]
transformed_function = StochasticTransformedFunction(sol,sol.prob.f,sol.prob.g)
transformed_function(du,u,p,tspan[2])
@test isapprox(du,(p[1]*u-p[2]^2*u), atol=1e-15)

transformed_function = StochasticTransformedFunction(sol,sol.prob.f,sol.prob.g,(du,u,p,t)-> (du.=p[2]^2*u))
transformed_function(du,u,p,tspan[2])
@test isapprox(du,(p[1]*u-p[2]^2*u), atol=1e-15)

#  non-diagonal noise torus
u0 = rand(2)
p =  rand(1)

fnd(u,p,t) = 0*u
function σnd(u,p,t)
  du = [cos(p[1])*sin(u[1])  cos(p[1])*cos(u[1])   -sin(p[1])*sin(u[2])   -sin(p[1])*cos(u[2])
   sin(p[1])*sin(u[1])   sin(p[1])*cos(u[1])    cos(p[1])*sin(u[2])   cos(p[1])*cos(u[2]) ]
  return du
end

prob = SDEProblem(fnd,σnd,u0,tspan,p,noise_rate_prototype=zeros(2,4))
sol = solve(prob,EM(),adaptive=false, dt=0.001, save_noise=true)

transformed_function = StochasticTransformedFunction(sol,sol.prob.f,sol.prob.g)
du2 = transformed_function(u0,p,tspan[2])
@test isapprox(du2,zeros(2), atol=1e-15)

transformed_function = StochasticTransformedFunction(sol,sol.prob.f,sol.prob.g,(u,p,t)->false*u)
du2 = transformed_function(u0,p,tspan[2])
@test isapprox(du2,zeros(2), atol=1e-15)

fnd!(du,u,p,t) = du .= false
function σnd!(du,u,p,t)
  du[1,1] = cos(p[1])*sin(u[1])
  du[1,2] = cos(p[1])*cos(u[1])
  du[1,3] = -sin(p[1])*sin(u[2])
  du[1,4] = -sin(p[1])*cos(u[2])
  du[2,1] = sin(p[1])*sin(u[1])
  du[2,2] = sin(p[1])*cos(u[1])
  du[2,3] = cos(p[1])*sin(u[2])
  du[2,4] = cos(p[1])*cos(u[2])
  return nothing
end

prob = SDEProblem(fnd!,σnd!,u0,tspan,p,noise_rate_prototype=zeros(2,4))
sol = solve(prob,EM(),adaptive=false, dt=0.001, save_noise=true)

du = zeros(size(u0))
u = sol.u[end]
transformed_function = StochasticTransformedFunction(sol,sol.prob.f,sol.prob.g)
transformed_function(du,u,p,tspan[2])
@test isapprox(du,zeros(2), atol=1e-15)

transformed_function = StochasticTransformedFunction(sol,sol.prob.f,sol.prob.g,(du,u,p,t)-> (du.=false*u))
transformed_function(du,u,p,tspan[2])
@test isapprox(du,zeros(2), atol=1e-15)

t = sol.t[end]

"""
Check compatibility of StochasticTransformedFunction with vjp for adjoints
"""

###
# Check general compatibility of StochasticTransformedFunction() with Zygote
###

using Zygote

# scalar case
Random.seed!(seed)
u0 = rand(1)
p = rand(2)
λ = rand(1)

_dy, back = Zygote.pullback(u0, p) do u, p
  vec(f(u, p, t)-p[2]^2*u)
end
∇1,∇2 = back(λ)

@test isapprox(∇1, (p[1]-p[2]^2)*λ, atol=1e-15)
@test isapprox(∇2, (@. [1,-2*p[2]]*u0*λ[1]), atol=1e-15)

prob = SDEProblem(f,σ,u0,tspan,p)
sol = solve(prob,SOSRI(),adaptive=false, dt=0.001, save_noise=true)
transformed_function = StochasticTransformedFunction(sol,sol.prob.f,sol.prob.g)

# Zygote doesn't allow nesting
_dy, back = Zygote.pullback(u0, p) do u, p
  vec(transformed_function(u, p, t))
end
@test_broken back(λ)
# @test isapprox(∇1, (p[1]-p[2]^2)*λ, atol=1e-15)
# @test isapprox(∇2, (@. [1,-2*p[2]]*u0*λ[1]), atol=1e-15)

transformed_function = StochasticTransformedFunction(sol,sol.prob.f,sol.prob.g, (u,p,t)->p[2]^2*u)
_dy, back = Zygote.pullback(u0, p) do u, p
  vec(transformed_function(u, p, t))
end
∇1,∇2 = back(λ)
@test isapprox(∇1, (p[1]-p[2]^2)*λ, atol=1e-15)
@test isapprox(∇2, (@. [1,-2*p[2]]*u0*λ[1]), atol=1e-15)

###
# Check general compatibility of StochasticTransformedFunction() with ReverseDiff
###

using ReverseDiff

# scalar

tape = ReverseDiff.GradientTape((u0, p, [t])) do u,p,t
  vec(f(u, p, t)-p[2]^2*u)
end

tu, tp, tt = ReverseDiff.input_hook(tape)

output = ReverseDiff.output_hook(tape)
ReverseDiff.unseed!(tu) # clear any "leftover" derivatives from previous calls
ReverseDiff.unseed!(tp)
ReverseDiff.unseed!(tt)

ReverseDiff.value!(tu, u0)
ReverseDiff.value!(tp, prob.p)
ReverseDiff.value!(tt, [t])
ReverseDiff.forward_pass!(tape)
ReverseDiff.increment_deriv!(output, λ)
ReverseDiff.reverse_pass!(tape)

@test isapprox(ReverseDiff.deriv(tu), (p[1]-p[2]^2)*λ, atol=1e-15)
@test isapprox(ReverseDiff.deriv(tp), (@. [1,-2*p[2]]*u0*λ[1]), atol=1e-15)

tape = ReverseDiff.GradientTape((u0, p, [t])) do u,p,t
  _dy, back = Zygote.pullback(u, p) do u, p
    vec(σ(u, p, t))
  end
  tmp1,tmp2 = back(_dy)
  return f(u, p, t) - vec(tmp1)
end

tu, tp, tt = ReverseDiff.input_hook(tape)

output = ReverseDiff.output_hook(tape)
ReverseDiff.unseed!(tu) # clear any "leftover" derivatives from previous calls
ReverseDiff.unseed!(tp)
ReverseDiff.unseed!(tt)

ReverseDiff.value!(tu, u0)
ReverseDiff.value!(tp, prob.p)
ReverseDiff.value!(tt, [t])
ReverseDiff.forward_pass!(tape)
ReverseDiff.increment_deriv!(output, λ)
ReverseDiff.reverse_pass!(tape)

@test isapprox(ReverseDiff.deriv(tu), (p[1]-p[2]^2)*λ, atol=1e-15)
@test isapprox(ReverseDiff.deriv(tp), (@. [1,-2*p[2]]*u0*λ[1]), atol=1e-15)

transformed_function = StochasticTransformedFunction(sol,sol.prob.f,sol.prob.g)

tape = ReverseDiff.GradientTape((u0, p, [t])) do u,p,t
  vec(transformed_function(u, p, first(t)))
end

tu, tp, tt = ReverseDiff.input_hook(tape)

output = ReverseDiff.output_hook(tape)
ReverseDiff.unseed!(tu) # clear any "leftover" derivatives from previous calls
ReverseDiff.unseed!(tp)
ReverseDiff.unseed!(tt)

ReverseDiff.value!(tu, u0)
ReverseDiff.value!(tp, prob.p)
ReverseDiff.value!(tt, [t])
ReverseDiff.forward_pass!(tape)
ReverseDiff.increment_deriv!(output, λ)
ReverseDiff.reverse_pass!(tape)

@test isapprox(ReverseDiff.deriv(tu), (p[1]-p[2]^2)*λ, atol=1e-15)
@test isapprox(ReverseDiff.deriv(tp), (@. [1,-2*p[2]]*u0*λ[1]), atol=1e-15)

# diagonal
Random.seed!(seed)
u0 = rand(3)
λ = rand(3)

_dy, back = Zygote.pullback(u0, p) do u, p
  vec(f(u, p, t)-p[2]^2*u)
end
∇1,∇2 = back(λ)

@test isapprox(∇1, (p[1]-p[2]^2)*λ, atol=1e-15)
@test isapprox(∇2[1], dot(u0,λ), atol=1e-15)
@test isapprox(∇2[2], -2*p[2]*dot(u0,λ), atol=1e-15)

tape = ReverseDiff.GradientTape((u0, p, [t])) do u,p,t
  vec(transformed_function(u, p, first(t)))
end

tu, tp, tt = ReverseDiff.input_hook(tape)

output = ReverseDiff.output_hook(tape)
ReverseDiff.unseed!(tu) # clear any "leftover" derivatives from previous calls
ReverseDiff.unseed!(tp)
ReverseDiff.unseed!(tt)

ReverseDiff.value!(tu, u0)
ReverseDiff.value!(tp, p)
ReverseDiff.value!(tt, [t])
ReverseDiff.forward_pass!(tape)
ReverseDiff.increment_deriv!(output, λ)
ReverseDiff.reverse_pass!(tape)

tmptp = ReverseDiff.deriv(tp)

@test isapprox(ReverseDiff.deriv(tu), (p[1]-p[2]^2)*λ, atol=1e-15)
@test isapprox(tmptp[1], dot(u0,λ), atol=1e-15)
@test isapprox(tmptp[2], -2*p[2]*dot(u0,λ), atol=1e-15)

# non-diagonal
Random.seed!(seed)
u0 = rand(2)
p =  rand(1)
λ =  rand(2)

_dy, back = Zygote.pullback(u0, p) do u, p
  vec(fnd(u, p, t))
end
∇1,∇2 = back(λ)

@test isapprox(∇1, zero(∇1), atol=1e-15)

prob = SDEProblem(fnd,σnd,u0,tspan,p,noise_rate_prototype=zeros(2,4))
sol = solve(prob,EM(),adaptive=false, dt=0.001, save_noise=true)
transformed_function = StochasticTransformedFunction(sol,sol.prob.f,sol.prob.g)

tape = ReverseDiff.GradientTape((u0, p, [t])) do u,p,t
  vec(transformed_function(u, p, first(t)))
end

tu, tp, tt = ReverseDiff.input_hook(tape)

output = ReverseDiff.output_hook(tape)
ReverseDiff.unseed!(tu) # clear any "leftover" derivatives from previous calls
ReverseDiff.unseed!(tp)
ReverseDiff.unseed!(tt)

ReverseDiff.value!(tu, u0)
ReverseDiff.value!(tp, p)
ReverseDiff.value!(tt, [t])
ReverseDiff.forward_pass!(tape)
ReverseDiff.increment_deriv!(output, λ)
ReverseDiff.reverse_pass!(tape)

@test isapprox(ReverseDiff.deriv(tu), zero(u0), atol=1e-15)
@test isapprox(ReverseDiff.deriv(tp), zero(p), atol=1e-15)

###
# Check Mutating functions
###
# scalar
Random.seed!(seed)
u0 = rand(1)
p = rand(2)
λ = rand(1)

prob = SDEProblem(SDEFunction(f!,σ!,analytic=linear_analytic),σ!,u0,tspan,p)
sol = solve(prob,SOSRI(),adaptive=false, dt=0.001, save_noise=true)

du = zeros(size(u0))
u = sol.u[end]
transformed_function = StochasticTransformedFunction(sol,sol.prob.f,sol.prob.g)

function inplacefunc!(du,u,p,t)
  du .= p[2]^2*u
  return nothing
end

tape = ReverseDiff.GradientTape((u0, p, [t])) do u,p,t
  du1 = similar(u, size(u))
  du2 = similar(u, size(u))
  f!(du1, u, p, first(t))
  inplacefunc!(du2, u, p, first(t))
  return vec(du1-du2)
end

tu, tp, tt = ReverseDiff.input_hook(tape) # u0

output = ReverseDiff.output_hook(tape) # p[1]*u0 -p[2]^2*u0
ReverseDiff.unseed!(tu) # clear any "leftover" derivatives from previous calls
ReverseDiff.unseed!(tp)
ReverseDiff.unseed!(tt)

ReverseDiff.value!(tu, u0)
ReverseDiff.value!(tp, p)
ReverseDiff.value!(tt, [t])
ReverseDiff.forward_pass!(tape)
ReverseDiff.increment_deriv!(output, λ)
ReverseDiff.reverse_pass!(tape)

@test isapprox(ReverseDiff.deriv(tu), (p[1]-p[2]^2)*λ, atol=1e-15) # -0.016562475307537294
@test isapprox(ReverseDiff.deriv(tp), (@. [1,-2*p[2]]*u0*λ[1]), atol=1e-15)  #[0.017478629739736098, -0.023103635221731166]

tape = ReverseDiff.GradientTape((u0, p, [t])) do u,p,t
  _dy, back = Zygote.pullback(u, p) do u, p
    out_ = Zygote.Buffer(similar(u))
    σ!(out_, u, p, t)
    vec(copy(out_))
  end
  tmp1,tmp2 = back(λ)
  du1 = similar(u, size(u))
  f!(du1, u, p, first(t))
  return  vec(du1-tmp1)
end

tu, tp, tt = ReverseDiff.input_hook(tape)

output = ReverseDiff.output_hook(tape)
ReverseDiff.unseed!(tu) # clear any "leftover" derivatives from previous calls
ReverseDiff.unseed!(tp)
ReverseDiff.unseed!(tt)

ReverseDiff.value!(tu, u0)
ReverseDiff.value!(tp, p)
ReverseDiff.value!(tt, [t])
ReverseDiff.forward_pass!(tape)
ReverseDiff.increment_deriv!(output, λ)
ReverseDiff.reverse_pass!(tape)

@test_broken isapprox(ReverseDiff.deriv(tu), (p[1]-p[2]^2)*λ, atol=1e-15)
@test_broken isapprox(ReverseDiff.deriv(tp), (@. [1,-2*p[2]]*u0*λ[1]), atol=1e-15)

tape = ReverseDiff.GradientTape((u0, p, [t])) do u1,p1,t1
  du1 = similar(u1, size(u1))
  transformed_function(du1, u1, p1, first(t1))
  return  vec(du1)
end

tu, tp, tt = ReverseDiff.input_hook(tape) # p[1]*u0

output = ReverseDiff.output_hook(tape) # p[1]*u0 -p[2]^2*u0
ReverseDiff.unseed!(tu) # clear any "leftover" derivatives from previous calls
ReverseDiff.unseed!(tp)
ReverseDiff.unseed!(tt)

ReverseDiff.value!(tu, u0)
ReverseDiff.value!(tp, p)
ReverseDiff.value!(tt, [t])
ReverseDiff.forward_pass!(tape)
ReverseDiff.increment_deriv!(output, λ)
ReverseDiff.reverse_pass!(tape)

@test_broken isapprox(ReverseDiff.deriv(tu), (p[1]-p[2]^2)*λ, atol=1e-15)
@test_broken isapprox(ReverseDiff.deriv(tp), (@. [1,-2*p[2]]*u0*λ[1]), atol=1e-15)

transformed_function = StochasticTransformedFunction(sol,sol.prob.f,sol.prob.g,(du,u,p,t)-> (du.=p[2]^2*u))

tape = ReverseDiff.GradientTape((u0, p, [t])) do u1,p1,t1
  du1 = similar(u1, size(u1))
  transformed_function(du1, u1, p1, first(t1))
  return  vec(du1)
end

tu, tp, tt = ReverseDiff.input_hook(tape) # p[1]*u0

output = ReverseDiff.output_hook(tape) # p[1]*u0 -p[2]^2*u0
ReverseDiff.unseed!(tu) # clear any "leftover" derivatives from previous calls
ReverseDiff.unseed!(tp)
ReverseDiff.unseed!(tt)

ReverseDiff.value!(tu, u0)
ReverseDiff.value!(tp, p)
ReverseDiff.value!(tt, [t])
ReverseDiff.forward_pass!(tape)
ReverseDiff.increment_deriv!(output, λ)
ReverseDiff.reverse_pass!(tape)

@test isapprox(ReverseDiff.deriv(tu), (p[1]-p[2]^2)*λ, atol=1e-15)
@test isapprox(ReverseDiff.deriv(tp), (@. [1,-2*p[2]]*u0*λ[1]), atol=1e-15)