text:
using CUDAdrv, CUDAnative, CuArrays
using RandomNumbers.MersenneTwisters

#Small function required by twister
mt_magic(y) = ((y % Int32) << 31 >> 31) & 0x9908b0df

#Function that runs in parallel
function kernel(backing_array, seedArr, outputArray, N_rand_to_make, N)

    #set some constants
    M = 397
    UPPER_MASK = 0x80000000
    LOWER_MASK = 0x7fffffff

    #Set tid and totalThreads
    tid = ((threadIdx().x) + (blockIdx().x - 1) * blockDim().x)
    totalThreads = blockDim().x * gridDim().x

    #Set MT19937 array (mt) and integer (mti) from seedArr[tid]
    backing_array[threadIdx().x, blockIdx().x, 1] = seedArr[tid]
    for i in 2:N
        backing_array[threadIdx().x, blockIdx().x, i] = 0x6c078965 * (backing_array[threadIdx().x, blockIdx().x, i-1] ⊻ (backing_array[threadIdx().x, blockIdx().x, i-1] >> 30)) + (i - 1) % UInt32
    end
    mti = N + 1

    #This while loop allows more random numbers to be generated than there are threads initialised on the GPU
    while tid <= N_rand_to_make
        ############## I do this inline so I can mutate MT19937 directly
        if mti > N
            for i in 1:N-M
                y = (backing_array[threadIdx().x, blockIdx().x, i] & UPPER_MASK) | (backing_array[threadIdx().x, blockIdx().x, i+1] & LOWER_MASK)
                backing_array[threadIdx().x, blockIdx().x, i] = backing_array[threadIdx().x, blockIdx().x, i + M] ⊻ (y >> 1) ⊻ mt_magic(y)
            end
            for i in N-M+1:N-1
                y = (backing_array[threadIdx().x, blockIdx().x, i] & UPPER_MASK) | (backing_array[threadIdx().x, blockIdx().x, i+1] & LOWER_MASK)
                backing_array[threadIdx().x, blockIdx().x, i] = backing_array[threadIdx().x, blockIdx().x, i + M - N] ⊻ (y >> 1) ⊻ mt_magic(y)
            end
            begin
                y = (backing_array[threadIdx().x, blockIdx().x, N] & UPPER_MASK) | (backing_array[threadIdx().x, blockIdx().x, 1] & LOWER_MASK)
                backing_array[threadIdx().x, blockIdx().x, N] = backing_array[threadIdx().x, blockIdx().x, M] ⊻ (y >> 1) ⊻ mt_magic(y)
            end
            mti = 1
        end
        k = backing_array[threadIdx().x, blockIdx().x, mti]
        k ⊻= (k >> 11)
        k ⊻= (k << 7) & 0x9d2c5680
        k ⊻= (k << 15) & 0xefc60000
        k ⊻= (k >> 18)

        mti += 1
        ############
        outputArray[tid] = k
        tid += totalThreads
    end

    return nothing
end

function gpuRand(threadNum, blockNum, randNum)

    #Make an array of random states to seed with
    N = 624
    seed_arr = []
    for i in 1:((randNum ÷ N) + 1)
        append!(seed_arr, MT19937().mt)

    end
    seed_array = CuArrays.CuArray(UInt32.(seed_arr))

    #Make other args
    outputArray = CuArrays.CuArray(fill(0, randNum))
    backing_array = CuArrays.CuArray(fill(UInt32(0), threadNum, blockNum, N))

    #Make random numbers on GPU
    @cuda threads = threadNum blocks = blockNum kernel(backing_array, seed_array, outputArray, randNum, N)

    return Array(outputArray)
end


text:
using OptionalUnits
using Test
using Unitful

@testset "OptionalUnits.jl" begin
    @testset "BasicFunctions" begin
        @optionalunits function f1(a::Unitful.Length → u"m", b::Unitful.Length → u"mm")
            a + b
        end
        @test f1(1.0, 2.0) == 1.002u"m"
        @test f1(1, 2u"m") == 3u"m"
        @test f1(1, 2.0u"m") == 3.0u"m"
        @test f1(1u"mm", 2) == 3u"mm"
        @test f1(1u"km", 2u"km") == 3u"km"
        @test_throws MethodError f1(1u"m", 2u"m/s")
    end
    @testset "BasicFunctionsVectors" begin
        @optionalunits function f2(a::Vector{Unitful.Length → u"m"}, b::Unitful.Length → u"mm")
            a[1] += b
            return a
        end
        @test f2([1.0,2.0], 2.0) == [1.002u"m",2.0u"m"]
        @test f2([1,2], 2u"m") == [3u"m",2u"m"]
        @test f2([1,2], 2.0u"m") == [3.0u"m",2.0u"m"]
        @test f2([1u"mm",2u"mm"], 2) == [3u"mm",2u"mm"]
        @test f2([1,2]u"km", 2u"km") == [3u"km",2u"km"]
        @test_throws MethodError f2([1u"m",2], 3)
    end
    @testset "FunctionsWithDefaults" begin
        @optionalunits function f3(a::Unitful.Length → u"m", b::Unitful.Length → u"mm"=1u"mm")
            a + b
        end
        @test f3(1.0) == 1.001u"m"
        @test f3(1, 2u"m") == 3u"m"
        @test f3(1u"mm") == 2u"mm"
        @test f3(1u"km", 2u"km") == 3u"km"
        @test_throws MethodError f3(1u"m", 2u"m/s")
    end
    @testset "FunctionsWithDefaultsVector" begin
        @optionalunits function f4(a::Unitful.Length → u"m", b::Vector{Unitful.Length → u"mm"}=[1u"mm",2u"mm"])
            b[1] += a
            return b
        end
        @test f4(1) == [1001u"mm",2u"mm"]
        @test f4(1, [2u"m",3u"m"]) == [3u"m", 3u"m"]
        @test f4(1u"mm") == [2u"mm",2u"mm"]
        @test f4(1u"km", [2u"km", 4u"mm"]) == [3u"km",4u"mm"]
        @test_throws MethodError f4(1u"m", [2u"m",2u"m/s"])
    end
    @testset "Struct" begin
        @optionalunits struct Point
            x::typeof(1u"m")
            y::typeof(1.0u"mm")
            size::Int
        end

        @test Point(1,1,1)==Point(1u"m",1.0u"mm",1)
        @test Point(2u"m",5u"m",2)==Point(2u"m",5000.0u"mm",2)

        @test_throws InexactError Point(1u"mm",1,1)
        @test_throws InexactError Point(1.5,1,1)
    end
    @testset "StructVector" begin
        @optionalunits struct MyVector
            values::Vector{typeof(1.0u"mm")}
        end

        @test MyVector([1.0,1.0,1.0]).values==MyVector([1.0u"mm",1.0u"mm",1.0u"mm"]).values
        @test MyVector([2u"m",5u"m",2u"mm"]).values==MyVector([2000.0u"mm",5000.0u"mm",2.0u"mm"]).values

        @test_throws AssertionError MyVector([1u"mm",1,1])
    end
end

ext: jl


text:
# utility functions

tofloat(array::AbstractArray{T}) where T = map(Float32, array)
tofloat(array::AbstractArray{Float32}) = array

fft_frequencies(samplerate::Real, nfft::Int) = collect(linspace(0f0, samplerate / 2f0, (nfft >> 1) + 1))

"""returns the number of frames when the signal is partitioned into overlapping frames"""
nframes(length::Int, framesize::Int, hopsize::Int) = div(length - framesize, hopsize) + 1

"""Provides a view of a signal partitioned into overlapping frames, similar to MATLAB's buffer"""
struct FrameView{T<:AbstractVector} <: AbstractVector{Vector}
    original::T
    framesize::Int
    hopsize::Int
    nframes::Int

    function FrameView(original, framesize, hopsize)
        # n = noverlap is a problem - the algorithm will not terminate.
        #@boundscheck((0 < hopsize <= framesize), error("hopsize must be between zero and framesize"))
        new{T}(original, framesize, hopsize, nframes(length(original), framesize, hopsize))
    end
end

function Base.getindex(x::FrameView, i::Int)
    #@boundscheck((i >= 1 && i <= x.k), Base.throw_boundserror(x, i))
    slice(x.original, (i-1) * x.hopsize + (1:x.framesize))
end

Base.size(x::FrameView) = (x.nframes,)

function Base.iterate(x::FrameView, i::Int)
   if i > x.nframes
       nothing
   else
       (x[i], i+1)
   end
end

"""returns the DCT filters"""
function dct(nfilters::Int, ninput::Int)
    basis = Array(Float32, nfilters, ninput)
    samples = (1f0:2f0:2ninput) * π / 2ninput
    for i = 1:nfilters
        basis[i, :] = cos(i * samples)
    end

    basis *= sqrt(2f0/ninput)
    basis
end

"""returns the number of times that the value changes its sign"""
function zero_crossings(array::AbstractVector{T}, size::Int = length(array), offset::Int = 0) where {T<:Real}
    result = 0
    previous = 0
    for i in offset .+ (1:size)
        number = array[i]
        sgn = number == 0 ? 0 : number > 0 ? 1 : -1
        if sgn != previous
            result += 1
        end
        previous = sgn
    end
    result
end

# units from SampledSignals

const frames = Hz*s
const FrameQuant = DimensionlessQuantity

"""
    inframes([Type,]quantity[, rate])

Translate the given quantity to a (unitless) number of time or frequency frames,
given a particular samplerate. Note that this isn't quantized to integer numbers
of frames. If given a `Type`, the result will first be coerced to the given type.

If the given quantity is Unitful, we use the given units. If it is not we assume
it is already a value in frames.

# Example

julia> inframes(0.5s, 44100Hz)
22050.0

julia> inframes(1000Hz, 2048/44100Hz)
46.439909297052154

"""
inframes(::Type{T}, frame::FrameQuant, rate=nothing) where T <: Integer =
    round(T, ustrip(uconvert(frames, frame)))
inframes(frame::FrameQuant, rate=nothing) = ustrip(uconvert(frames, frame))
inframes(::Type{T}, time::Unitful.Time, rate) where T <: Integer =
    round(T, inseconds(time)*inHz(rate))
inframes(time::Unitful.Time, rate) = inseconds(time)*inHz(rate)
inframes(::Type{T}, freq::Unitful.Frequency, rate) where T <: Integer =
    round(T, inHz(freq)*inseconds(rate))
inframes(freq::Unitful.Frequency, rate) = inHz(freq)*inseconds(rate)
inframes(::Type, frame::Unitful.AbstractQuantity) = error("Unknown sample rate")
inframes(frame::Unitful.AbstractQuantity) = error("Unknown sample rate")
inframes(::Type{T}, frame::Real) where T = T(frame)
inframes(frame::Real) = frame

"""
    inHz(quantity[, rate])

Translate a particular quantity (usually a frequency) to a (unitless) value in
Hz.

If the given quantity is Unitful, we use the given units. If it is not we assume
it is already a value in Hz.

For some units (e.g. frames) you will need to specify a sample rate:

# Example

julia> inHz(1.0kHz)
1000.0

"""
inHz(x::Unitful.Frequency, rate=nothing) = ustrip(uconvert(Hz, x))
inHz(x::FrameQuant) = error("Unknown sample rate")
# assume we have a spectrum buffer with a sample rate in seconds
inHz(x::FrameQuant, rate) = inHz(inframes(x) / rate)
inHz(x::Real, rate) = x
inHz(x::Real) = x

"""
   inseconds(quantity[, rate])

Translate a particular quantity (usually a time) to a (unitless) value in
seconds.

If the given quantity is Unitful, we use the given units. If it is not we assume
it is already a value in seconds.

For some units (e.g. frames) you will need to specify a sample rate:

# Examples
julia> inseconds(50.0ms)
0.05

julia> inseconds(441frames, 44100Hz)
0.01

"""
inseconds(x::Unitful.Time, rate=nothing) = ustrip(uconvert(s,x))
inseconds(x::FrameQuant) = error("Unknown sample rate")
# assume we have a time buffer with sample rate in hz
inseconds(x::FrameQuant, rate) = inseconds(inframes(x) / rate)
inseconds(x::Real, rate) = x
inseconds(x::Real) = x

# inseconds(x, rate::Quantity) = inseconds(x,inHz(rate))
# inseconds(x::FrameQuant, rate::Real) = (ustrip(x) / rate)
# inseconds(x::Quantity, rate::Real) = ustrip(uconvert(s,x))


text:
# https://juliadocs.github.io/Documenter.jl/stable/man/guide/#Package-Guide
# push!(LOAD_PATH,"../src/")

# Run these locally to build docs/build folder:
# julia --color=yes --project=docs/ -e 'using Pkg; pkg\"add https://github.com/bankofcanada/TimeSeriesEcon.jl.git; add https://github.com/bankofcanada/ModelBaseEcon.jl.git\"'
# julia --color=yes --project=docs/ -e 'using Pkg; Pkg.develop(PackageSpec(path=pwd())); Pkg.instantiate()'
# julia --project=docs/ docs/make.jl

using Documenter, StateSpaceEcon

# Workaround for JuliaLang/julia/pull/28625
if Base.HOME_PROJECT[] !== nothing
    Base.HOME_PROJECT[] = abspath(Base.HOME_PROJECT[])
end

makedocs(sitename = "StateSpaceEcon.jl",
         format = Documenter.HTML(prettyurls = get(ENV, "CI", nothing) == "true"),
         modules = [StateSpaceEcon],
         doctest = false,
         pages = [
        "Reference" => "index.md",

        # "Examples" => "examples.md"
    ]
)

deploydocs(
    repo = "github.com/bankofcanada/StateSpaceEcon.jl.git",
)


text:
export perturbative_transmon_Ueigen

# The following matrices were taken from the source of arXiv:1706.06566.
# The regexes below were used to construct a matrix in Julia syntax
# from the LaTeX source syntax:

# \\ => newline
# \^\{([0-9]+)\}([^\}]) => ^$1*$2
# \\sqrt\{([0-9]+)\} => (√$1)
# \\frac\{([√0-9\^\*\{\}\(\)]+)\}\{([√0-9\^\*\{\}\(\)]+)\} => ($1)/($2)
# \^\{([0-9]+)\} => ^$1
# \)\( => )*(

# Additionally, Julia didn't recognize the minus signs (incorrect glyph?)
# I replaced those.

const U0 = collect(adjoint([
 1 0 0 0 0
 0 1 0 0 0
 0 0 1 0 0
 0 0 0 1 0
 0 0 0 0 1
 0 0 0 0 0
 0 0 0 0 0
 0 0 0 0 0
 0 0 0 0 0
 0 0 0 0 0
 0 0 0 0 0
 0 0 0 0 0
 0 0 0 0 0
 0 0 0 0 0
 0 0 0 0 0
 0 0 0 0 0
 0 0 0 0 0
 0 0 0 0 0
 0 0 0 0 0
 0 0 0 0 0
 0 0 0 0 0
 0 0 0 0 0
 0 0 0 0 0
 0 0 0 0 0
 0 0 0 0 0]))

const U1 = collect(adjoint([
 0 0 -(1)/((√2)) 0 -(1)/((√2)^3*(√3))
 0 0 0 -(5)/((√2)*(√3)) 0
 (1)/((√2)) 0 0 0 -(7)/((√3))
 0 (5)/((√2)*(√3)) 0 0 0
 (1)/((√2)^3*(√3)) 0 (7)/((√3)) 0 0
 0 ((√5))/((√2)^3*(√3)) 0 (3(√5))/(1) 0
 0 0 ((√5))/((√2)^3) 0 (11(√5))/((√2)*(√3))
 0 0 0 ((√35))/((√2)^3*(√3)) 0
 0 0 0 0 ((√35))/(2(√3))
 0 0 0 0 0
 0 0 0 0 0
 0 0 0 0 0
 0 0 0 0 0
 0 0 0 0 0
 0 0 0 0 0
 0 0 0 0 0
 0 0 0 0 0
 0 0 0 0 0
 0 0 0 0 0
 0 0 0 0 0
 0 0 0 0 0
 0 0 0 0 0
 0 0 0 0 0
 0 0 0 0 0
 0 0 0 0 0])) / 2^3

const U2 = collect(adjoint([
-(13)/(2^8*3) 0 -(5)/((√2)^7*3) 0 ((√3))/((√2)^9)
 0 -(35)/(2^8) 0 -(13)/((√2)^7*(√3)) 0
 (13)/((√2)^11*3) 0 -(419)/(2^8*3) 0 -(67(√3))/(2^6)
 0 (37)/((√2)^11*(√3)) 0 -(405)/(2^8) 0
 (1)/((√2)^5*(√3)) 0 (145)/(2^6*(√3)) 0 -(961)/(2^8)
 0 (41)/((√2)^9*(√15)) 0 (79(√5))/(2^6) 0
 (23)/(2^6*3(√5)) 0 (103)/((√2)^7*3(√5)) 0 (29(√15))/((√2)^9)
 0 (11(√7))/(2^6*(√5)) 0 (103(√7))/((√2)^9*(√15)) 0
 ((√35))/((√2)^15*3) 0 (43(√7))/(2^5*3(√5)) 0 (3(√21))/((√5))
 0 ((√35))/((√2)^15) 0 (53(√7))/(2^5*(√15)) 0
 0 0 (5(√7))/((√2)^15) 0 (21(√21))/((√2)^11)
 0 0 0 (5(√77))/((√2)^15*(√3)) 0
 0 0 0 0 (5(√77))/((√2)^15)
 0 0 0 0 0
 0 0 0 0 0
 0 0 0 0 0
 0 0 0 0 0
 0 0 0 0 0
 0 0 0 0 0
 0 0 0 0 0
 0 0 0 0 0
 0 0 0 0 0
 0 0 0 0 0
 0 0 0 0 0
 0 0 0 0 0])) / 2^2

 const U3 = collect(adjoint([
 -(17)/(2^3*3) 0 -(725)/((√2)^11*3) 0 (1747)/((√2)^13*(√3))
  0 -(113)/(2^2*3) 0 -(869(√3))/((√2)^11) 0
  (347)/((√2)^11*3) 0 -(1241)/(2^3*3) 0 -(6217)/(2^5*(√3))
  0 (515(√3))/((√2)^11) 0 -(189)/(1) 0
  (605)/((√2)^13*(√3)) 0 (3793)/(2^5*(√3)) 0 -(2165)/(2^2)
  0 (2009(√5))/((√2)^13*(√3)) 0 (2225(√5))/(2^5) 0
  (409)/(2^4*3(√5)) 0 (68981)/((√2)^13*3(√5)) 0 (44689)/((√2)^11*(√15))
  0 (461(√5))/(2^4*(√7)) 0 (8663(√35))/((√2)^13*(√3)) 0
  (239)/((√2)^5*3(√35)) 0 (52681)/(2^5*3(√35)) 0 (36119(√7))/(2^5*(√15))
  0 (127(√5))/((√2)^3*3(√7)) 0 (52763(√5))/(2^5*3(√21)) 0
  (31(√7))/(2^5*3) 0 (1963)/((√2)^5*3(√7)) 0 (202283)/((√2)^9*3(√21))
  0 (41(√77))/(2^5*3) 0 (233(√11))/((√2)*(√21)) 0
  (5(√77))/(2^6*3(√3)) 0 (17(√231))/((√2)^9) 0 (3769(√11))/((√2)^5*3(√7))
  0 (5(√1001))/(2^6*3(√3)) 0 (61(√1001))/((√2)^9*3) 0
  0 0 (35(√143))/(2^6*3(√3)) 0 (497(√143))/(2^5*3)
  0 0 0 (35(√715))/(2^6*3(√3)) 0
  0 0 0 0 (35(√715))/(2^5*3(√3))
  0 0 0 0 0
  0 0 0 0 0
  0 0 0 0 0
  0 0 0 0 0
  0 0 0 0 0
  0 0 0 0 0
  0 0 0 0 0
  0 0 0 0 0])) / 2^8

const U4 = collect(adjoint([
-(62333)/((√2)^21*15) 0 -(9911)/(2^5*15) 0 (39007)/(2^6*3(√3))
 0 -(1313189)/((√2)^21*15) 0 -(137483)/(2^5*15(√3)) 0
 (16181)/(2^7*15) 0 -(651579)/((√2)^21) 0 -(3086339)/((√2)^15*15(√3))
 0 (110401)/(2^7*5(√3)) 0 -(44292461)/((√2)^21*15) 0
 (269(√3))/(2^5) 0 (227829(√3))/((√2)^15*5) 0 -(149454457)/((√2)^21*15)
 0 (8281(√5))/(2^6*(√3)) 0 (254291(√5))/((√2)^15*3) 0
 (7601)/((√2)^15*(√5)) 0 (5833)/(3(√5)) 0 (55231)/(2^2*3(√15))
 0 (85841(√7))/((√2)^15*3(√5)) 0 (5911741)/(2^6*3(√105)) 0
 (164443)/(2^9*3(√35)) 0 (1543667)/((√2)^13*3(√35)) 0 (11114269)/((√2)^11*3(√105))
 0 (4110871)/(2^9*9(√35)) 0 (1935007(√5))/((√2)^13*3(√21)) 0
 (135697)/((√2)^13*45(√7)) 0 (6278647(√7))/(2^9*45) 0 (1265399(√7))/(2^6*3(√3))
 0 (283201(√11))/((√2)^13*45(√7)) 0 (26982119(√11))/(2^9*15(√21)) 0
 (2759(√11))/((√2)^9*15(√21)) 0 (199855(√11))/(2^7*3(√21)) 0 (135190589(√11))/(2^9*45(√7))
 0 (8843(√143))/((√2)^11*15(√21)) 0 (1592063(√143))/(2^7*45(√7)) 0
 (91(√143))/(2^7*3(√3)) 0 (3217(√143))/((√2)^7*15(√3)) 0 (73963(√143))/((√2)^5*45)
 0 (343(√715))/(2^7*9) 0 (17593(√143))/((√2)^11*3(√15)) 0
 (35(√715))/(2^9*9) 0 (413(√715))/((√2)^11*9) 0 (11509(√143))/((√2)^7*3(√15))
 0 (35(√12155))/(2^9*9) 0 (161(√12155))/((√2)^11*3(√3)) 0
 0 0 (35(√12155))/(2^9*3) 0 (553(√12155))/(2^6*3(√3))
 0 0 0 (35(√230945))/(2^9*3(√3)) 0
 0 0 0 0 (175(√46189))/(2^9*3(√3))
 0 0 0 0 0
 0 0 0 0 0
 0 0 0 0 0
 0 0 0 0 0])) / (√2)^21

const U5 = collect(adjoint([
-(153917)/((√2)^9*45) 0 -(10765201)/(2^8*45) 0 (202767029)/(2^9*15(√3))
 0 -(41527)/((√2)^3*3) 0 -(10916291)/(2^8*3(√3)) 0
 (4600889)/(2^8*45) 0 -(50029111)/((√2)^9*45) 0 -(51091643)/((√2)^15*5(√3))
 0 (9434411)/(2^8*3(√3)) 0 -(28185541)/((√2)^7*9) 0
 (9851179)/(2^9*15(√3)) 0 (188993981)/((√2)^15*15(√3)) 0 -(187170527)/((√2)^7*15)
 0 (25184897)/(2^9*(√15)) 0 (394451371)/((√2)^15*9(√5)) 0
 (994363)/((√2)^11*9(√5)) 0 (887375369)/(2^9*9(√5)) 0 (204498013)/(2^8*3(√15))
 0 (1873273(√7))/((√2)^11*3(√5)) 0 (5612971213)/(2^9*3(√105)) 0
 (558181)/(2^3*9(√35)) 0 (530587109)/((√2)^13*9(√35)) 0 (3942561341)/((√2)^17*(√105))
 0 (350365(√5))/(29(√7)) 0 (451113919)/((√2)^13*(√105)) 0
 (4552733)/((√2)^13*15(√7)) 0 (99431659)/(2^3*45(√7)) 0 (7205463341)/(2^6*15(√21))
 0 (94247731)/((√2)^13*9(√77)) 0 (7581187(√11))/(2^2*3(√21)) 0
 (117428699)/((√2)^15*45(√231)) 0 (3349510421)/(2^6*45(√231)) 0 (301458499(√11))/(2^4*15(√7))
 0 (6249499(√13))/((√2)^15*(√231)) 0 (449418449(√13))/(2^6*9(√77)) 0
 (108659(√143))/(2^5*45(√3)) 0 (555486781(√13))/((√2)^15*45(√33)) 0 (1352574373(√13))/((√2)^13*15(√11))
 0 (203167(√143))/(2^5*9(√5)) 0 (108310393(√13))/((√2)^15*(√165)) 0
 (2633(√143))/(2^3*9(√5)) 0 (59661(√715))/((√2)^13) 0 (1051619747(√13))/((√2)^15*3(√165))
 0 (491(√2431))/(9(√5)) 0 (454661(√2431))/((√2)^13*(√15)) 0
 (329(√12155))/((√2)^13*9) 0 (16369(√2431))/(2^3*9(√5)) 0 (2945699(√2431))/(2^6*3(√15))
 0 (133(√230945))/((√2)^13*3) 0 (401(√138567))/(2^2*(√5)) 0
 (35(√46189))/((√2)^15*9) 0 (2345(√46189))/(2^6*9) 0 (3071(√46189))/(2^3*(√3))
 0 (35(√323323))/((√2)^15*3(√3)) 0 (2695(√323323))/(2^6*9) 0
 0 0 (385(√29393))/((√2)^15*3(√3)) 0 (11165(√29393))/((√2)^13*3)
 0 0 0 (385(√676039))/((√2)^15*9) 0
 0 0 0 0 (385(√676039))/(2^7*3(√3))])) / (√2)^33

"""
   perturbative_transmon_Ueigen(ξ)
Returns a matrix representing a unitary transformation taking Fock states to
perturbative transmon eigenstates. The transformation is fifth-order in ξ and
uses the first 5 transmon eigenstates.
"""
perturbative_transmon_Ueigen(ξ) = @evalpoly(ξ, U0, U1, U2, U3, U4, U5)


text:
# Everying Euler trails/tours

export euler

# This function finds an Euler trail/tour given initial
# vertices. Returns the trail if it exists of [] if it does not.

"""
`euler(G,u,v)` finds an Euler trail in `G` starting at `u` and ending
at `v` returned as a list of vertices (in the order they are visited
on the trail).

`euler(G,u)` finds an Euler tour beginning and ending at
`u`. Alternatively, call `euler(G)` and the initial/final vertex will
be selected for you.

Note: The algorithm will succeed even if there are isolated vertices
in the graph (just don't choose an isolated vertex as the first/last).

If no Euler trail/tour exists, an empty list is returned.
"""
function euler(G::SimpleGraph{T}, u::T, v::T) where {T}
    notrail = T[]

    # perform basic checks:
    if !(has(G, u) && has(G, v))
        error("One or both of these vertices is not in this graph")
    end

    # special case: if all vertices have degree zero
    if all([deg(G, x) == 0 for x in G.V])
        if u == v
            return [u]
        else
            return notrail
        end
    end

    # if either vertex has degree zero, we're doomed
    if deg(G, u) == 0 || deg(G, v) == 0
        return notrail
    end

    # vertex degree checks
    if u == v
        for x in G.V
            if deg(G, x) % 2 == 1
                return notrail
            end
        end
    else        # u != v
        for x in G.V
            if x == u || x == v
                if deg(G, x) % 2 == 0
                    return notrail
                end
            else
                if deg(G, x) % 2 == 1
                    return notrail
                end
            end
        end
    end

    # Remove isolates and check for connectivity
    GG = trim(G)
    if !is_connected(GG)
        return notrail
    end

    # all tests have been satisfied. Now find the trail in GG using a
    # helper function.
    return euler_work!(GG, u)
end

# special case: find an Euler tour from a specified vertex
function euler(G::SimpleGraph{T}, u::T) where {T}
    return euler(G, u, u)
end

# special case: find any Euler tour. If the graph is connected, any
# vertex will do but if there are isolated vertices, we don't want to
# pick one of those!
function euler(G::SimpleGraph{T}) where {T}
    if cache_check(G, :euler)
        return cache_recall(G, :euler)
    end
    if NV(G) == 0
        return T[]
    end

    verts = vlist(G)

    # search for a vertex that isn't isolated
    for u in verts
        if deg(G, u) > 0
            return euler(G, u, u)
        end
    end

    # if we reach here, the graph only has isolated vertics. Let it
    # work from any isolated vertex (and likely fail).
    u = verts[1]
    tour = euler(G, u, u)
    cache_save(G, :euler, tour)
    return tour
end

# private helper function for euler()
function euler_work!(G::SimpleGraph{T}, u::T) where {T}
    trail = T[]

    while true
        # if last vertex
        if NV(G) == 1
            append!(trail, [u])
            return trail
        end

        # get the neighbors of u
        Nu = G[u]

        # if only one neighbor delete and move on
        if length(Nu) == 1
            w = Nu[1]
            SimpleGraphs.delete!(G, u)
            append!(trail, [u])
            u = w
        else
            for w in Nu
                if !is_cut_edge(G, u, w)
                    SimpleGraphs.delete!(G, u, w)
                    append!(trail, [u])
                    u = w
                    break
                end
            end
        end
    end
    error("This can't happen")
end


text:
module BISC195Bioinformatics

export normalizeDNA,
        composition,
        gc_content,
        complement,
        reverse_complement,
        parse_fasta,
        LongDNASeq,
        uniquekmers,
        distance

using BioSequences
import BioSequences: composition, gc_content, complement, reverse_complement, LongDNASeq

"""
    normalizeDNA(::AbstractString)

Ensures that a sequence only contains valid bases
(or `'N'` for unknown bases).
Returns a String.
"""
function normalizeDNA(seq)
    seq = uppercase(string(seq))
    for base in seq
        # note: `N` indicates an unknown base
        occursin(base, "AGCTN") || error("invalid base $base")
    end
    return LongDNASeq(seq) # change to `return LongDNASeq(seq)` if you want to try to use BioSequences types
end

# Your code here.
# Don't forget to export your functions!

function composition(sequence::AbstractString)
    BioSequences.composition(LongDNASeq(sequence))
end

function gc_content(sequence::AbstractString)
    BioSequences.gc_content(LongDNASeq(sequence))
end

function complement(sequence::AbstractString)
    BioSequences.complement(LongDNASeq(sequence))
end

function complement(base::Char)
    nt = convert(DNA, base)
    Char(BioSequences.complement(nt))
end

function reverse_complement(sequence::AbstractString)
    BioSequences.reverse_complement(LongDNASeq(sequence))
end

function parse_fasta(path)
    headers = []
    sequences = []
    str = ""
    for line in eachline(path)
        if '>' ∈ line
            push!(headers, string(lstrip(line, '>')))
            if str != ""
                push!(sequences, LongDNASeq(str))
            end
            str = ""
        elseif line != ""
            # for base in line
            #     occursin(base, "AGCTNagctn") || error("invalid base $base")
            # end
            str = str * line
        end
    end
    push!(sequences, LongDNASeq(str))
    parsedfile = (headers, sequences)
    return parsedfile
end

function uniquekmers(seq, k)
    sequence = String(seq)
    1 <= k <= length(sequence) || error("k must be a positive integer less than the length of the sequence")
    kmers = []
    stopindex = length(sequence)-(k-1)
    for i in 1:stopindex
        kmer = sequence[i:i+k-1]
        if match(r"[^ACGT]", kmer) != nothing # if kmer contains any ambiguous base
            continue
        end
        push!(kmers, LongDNASeq(kmer))
    end
    return Set(kmers)
end

function distance(set1, set2)
    1 - (length(intersect(set1, set2)) / length(union(set1, set2)))
end

end # module BISC195Bioinformaticsp