text:
# ------------------------------------
# read noumber of nodes
# ------------------------------------
function read_nodenum(skipnum)
    """
    xmax: number of nodes in xi direction
    ymax: number of nodes in eta direction
    zmax: number of nodes in zeta direction
    """
    fff = []
    open("grid/nodesnum", "r") do f
        fff = read(f,String)
    end
    fff = split(fff,"\n",keepempty=false)
    num_nodes = length(fff) - skipnum

    for i in 1+skipnum:length(fff)
        fff[i] = replace(fff[i]," \r" => "")
    end

    temp = split(fff[2]," ")
    xmax = parse(Int64,temp[1])
    ymax = parse(Int64,temp[2])
    zmax = parse(Int64,temp[3])

    return xmax, ymax, zmax
end

# ------------------------------------
# read nodes
# ------------------------------------
function read_nodes(skipnum, xmax, ymax, zmax)
    fff = []
    open("grid/nodes", "r") do f
        fff = read(f,String)
    end
    fff = split(fff,"\n",keepempty=false)
    num_nodes = length(fff) - skipnum

    for i in 1+skipnum:length(fff)
        fff[i] = replace(fff[i]," \r" => "")
    end

    nodes = zeros(xmax, ymax, zmax, 3)
    for i in 1:num_nodes
        temp = split(fff[i+skipnum]," ")

        xnum = parse(Int64,temp[1])
        ynum = parse(Int64,temp[2])
        znum = parse(Int64,temp[3])
        nodes[xnum,ynum,znum,1] = parse(Float64,temp[4])
        nodes[xnum,ynum,znum,2] = parse(Float64,temp[5])
        nodes[xnum,ynum,znum,3] = parse(Float64,temp[6])
    end
    return nodes
end

# ------------------------------------
# read vecAx and vecAy
# ------------------------------------
function read_vecA(skipnum, xmax, ymax, zmax)

    #-------------------------------------
    # vecAx

    fff = []
    open("grid/vecAx", "r") do f
        fff = read(f,String)
    end
    fff = split(fff,"\n",keepempty=false)
    num_cell = length(fff) - skipnum

    for i in 1+skipnum:length(fff)
        fff[i] = replace(fff[i]," \r" => "")
    end

    vecAx = zeros(xmax, ymax-1, zmax-1, 3)
    for i in 1:num_cell
        temp = split(fff[i+skipnum]," ")

        x = parse(Int64,temp[1])
        y = parse(Int64,temp[2])
        z = parse(Int64,temp[3])
        vecAx[x,y,z,1] = parse(Float64,temp[4])
        vecAx[x,y,z,2] = parse(Float64,temp[5])
        vecAx[x,y,z,3] = parse(Float64,temp[6])
    end

    #-------------------------------------
    # vecAy
    fff = []
    open("grid/vecAy", "r") do f
        fff = read(f,String)
    end
    fff = split(fff,"\n",keepempty=false)
    num_cell = length(fff) - skipnum

    for i in 1+skipnum:length(fff)
        fff[i] = replace(fff[i]," \r" => "")
    end

    vecAy = zeros(xmax-1, ymax, zmax-1, 3)
    for i in 1:num_cell
        temp = split(fff[i+skipnum]," ")

        x = parse(Int64,temp[1])
        y = parse(Int64,temp[2])
        z = parse(Int64,temp[3])
        vecAy[x,y,z,1] = parse(Float64,temp[4])
        vecAy[x,y,z,2] = parse(Float64,temp[5])
        vecAy[x,y,z,3] = parse(Float64,temp[6])
    end

    #-------------------------------------
    # vecAz
    fff = []
    open("grid/vecAz", "r") do f
        fff = read(f,String)
    end
    fff = split(fff,"\n",keepempty=false)
    num_cell = length(fff) - skipnum

    for i in 1+skipnum:length(fff)
        fff[i] = replace(fff[i]," \r" => "")
    end

    vecAz = zeros(xmax-1, ymax-1, zmax, 3)
    for i in 1:num_cell
        temp = split(fff[i+skipnum]," ")

        x = parse(Int64,temp[1])
        y = parse(Int64,temp[2])
        z = parse(Int64,temp[3])
        vecAz[x,y,z,1] = parse(Float64,temp[4])
        vecAz[x,y,z,2] = parse(Float64,temp[5])
        vecAz[x,y,z,3] = parse(Float64,temp[6])
    end

    return vecAx, vecAy, vecAz
end

function read_allgrid()
    skip = 1
    xmax, ymax, zmax    = read_nodenum(skip)
    nodes               = read_nodes(skip, xmax, ymax, zmax)
    vecAx, vecAy, vecAz = read_vecA(skip, xmax, ymax, zmax)
    println("fin reading grid")
    return xmax, ymax, zmax, nodes, vecAx, vecAy, vecAz
end


text:
function mandelbrot(z::Complex{T}, numiters::Integer) where {T <: Real}
    z_i = copy(z)
    if abs(z) > 2
        return 0
    end
    for escapeiter in 1:numiters
        z_i = z_i * z_i + z
        if abs(z_i) > 2
            return escapeiter
        end
    end
    return numiters
end

function generateranges(numxrange)

    aspectratio = 16 / 9
    xrange = LinRange(-2.8, 1.6, numxrange)
    xrangeextent = xrange.stop - xrange.start
    yrangeextent = xrangeextent / aspectratio
    yrange = LinRange(-yrangeextent / 2, yrangeextent / 2, Integer(floor(numxrange / aspectratio)))

    return xrange, yrange
end

function generatezoomranges(centerx, centery, scalefactor, numyrange, aspectratio=16/9)
    expscale = 2 ^ -scalefactor
    yrangeextent = expscale
    xrangeextent = expscale * aspectratio
    numxrange = Int(floor(numyrange * aspectratio))
    xrange = LinRange(centerx - xrangeextent / 2, centerx + xrangeextent / 2, numxrange)
    yrange = LinRange(centery - yrangeextent / 2, centery + yrangeextent / 2, numyrange)
    return xrange, yrange
end

function plotmandelbulb(xrange, yrange, numiters::Integer, plot=true, showaxis=true)
    complexrange = [Complex{Float64}(x, y) for x in xrange, y in yrange]
    mandelbrotupto(z) = 1 - mandelbrot(z, numiters) / numiters
    mandelbrotimg = mandelbrotupto.(complexrange)

    if plot
        scene = Makie.image(xrange, yrange, mandelbrotimg, show_axis=showaxis, scale_plot=false)
        return mandelbrotimg, scene
    else
        return mandelbrotimg
    end

end

function rendermandelbulbpng(xrange, yrange, numiters::Integer, savename::AbstractString) #, plot=true, showaxis=true)
    mandelbrotimg = rendermandelbulbimg(xrange, yrange, numiters)
    save(savename, mandelbrotimg)
end

function rendermandelbulbimg(xrange, yrange, numiters::Integer) #, plot=true, showaxis=true)
    mandelbrotimg = zeros(UInt8, yrange.len, xrange.len)
    for (j, y) in enumerate(yrange), (i, x) in enumerate(xrange)
        mandelxy = 255 * (1 - mandelbrot(Complex{Float64}(x, y), numiters) / numiters)
        mandelbrotimg[j, i] = UInt8(floor(mandelxy))
    end
    return mandelbrotimg
end

function rendermandelbulbvideo(centerx::Number, centery::Number, scalefactorstart::Number,
     scalefactorend::Number, numiters::Integer, fps::Integer, seconds::Number, height::Integer, aspectratio::Number, savename::AbstractString)
    time = 0:1/fps:seconds
    numsteps = length(time)
    mandels = Array{Array{UInt8, 2}, 1}(undef, numsteps)
    for (i, t) in enumerate(time)
        scalefactor = (scalefactorend - scalefactorstart) * (t / seconds) + scalefactorstart
        xrange, yrange = generatezoomranges(centerx, centery, scalefactor, height, aspectratio)
        mandels[i] = rendermandelbulbimg(xrange, yrange, numiters)
    end
    props = [:priv_data => ("cruf"=>"22", "preset"=>"medium")]
    encodevideo(savename, mandels, framerate=fps, AVCodecContextProperties=props)

end

function rendermandelbrotimagecuda(centerx::Number, centery::Number, scalefactor::Number,
     numiters::Integer, height::Integer, aspectratio::Number)
    xrange, yrange = generatezoomranges(centerx, centery, scalefactor, height, aspectratio)
    complexrange = [Complex{Float64}(x, y) for x in xrange, y in yrange]
    mandelimg = map(z->mandelbrot(z, numiters) / numiters, complexrange)

    return xrange, yrange, mandelimg
end

function examplesettings()
    centerx = -0.5609882
    centery = 0.6409865
    centerx = 0.0
    centery = 0.0
    scalefactorstart = 0.0
    scalefactorend = 4.0
    height = 480
    #numiters = 200
    numiters = 250
    #numiters = 40
    fps = 10
    seconds = 10.0
    aspectratio = 16/9
    width = Int(floor(360 * aspectratio))
    #scene, layout = MakieLayout.layoutscene(1, 1, padding=100, show_axis=false)
    #ax = layout[1, 1] = MakieLayout.LAxis(scene, show_axis=false)

    #scene = Makie.Scene(raw = true, resolution = (width, height))
    scene = Makie.Scene()
    xrange, yrange, mandelimg = rendermandelbrotimagecuda(centerx, centery, scalefactorstart, numiters, height, aspectratio)
    #image = Makie.image!(scene, xrange, yrange, CUDA.CuArray(mandelimg), show_axis=false)
    #layout[1, 1] = Makie.image!(scene, xrange, yrange, mandelimg, show_axis=false, padding=0)
    #image = Makie.image!(ax, xrange, yrange, mandelimg, show_axis=false)
    #u = Makie.Node(0.0)
    #mandelimgnode = Makie.lift(x->mandelimg, u)
    mandelimgcolor = [ColorTypes.RGBA(Float32(color), Float32(color), Float32(color)) for color in mandelimg]
    image = Makie.image!(scene, mandelimgcolor, show_axis=false)
    #image = rendermandelbrotimageanimation(image, centerx, centery, scalefactorstart, scalefactorend, numiters, fps, seconds, height, aspectratio)
    return centerx, centery, scalefactorstart, scalefactorend, height, numiters, fps, seconds, aspectratio, scene#, u, mandelimgnode, mandelimg
    # centerx, centery, scalefactorstart, scalefactorend, height, numiters, fps, seconds, aspectratio, image = examplesettings(); image
end

function rendermandelbrotimageanimation(image, centerx::Number, centery::Number, scalefactorstart::Number,
     scalefactorend::Number, numiters::Integer, fps::Integer, seconds::Number, height::Integer, aspectratio::Number)
    time = 0:1/fps:seconds
    numsteps = length(time)
    for (i, t) in enumerate(time)
        scalefactor = (scalefactorend - scalefactorstart) * (t / seconds) + scalefactorstart
        xrange, yrange, mandelimg = rendermandelbrotimagecuda(centerx, centery, scalefactor, numiters, height, aspectratio)
        image.plots[1][:image][] = mandelimg
        #Makie.update!(image)
        sleep(0.0001)
    end
    return image
end

# includet("mandelbrot.jl"); centerx, centery, scalefactorstart, scalefactorend, height, numiters, fps, seconds, aspectratio, image = examplesettings(); image
# includet("mandelbrot.jl"); @time rendermandelbrotimageanimation2(image, centerx, centery, scalefactorstart, scalefactorend, numiters, fps, seconds, height, aspectratio, "")
function rendermandelbrotimageanimation2(image, centerx::Float64, centery::Float64, scalefactorstart::Number,
     scalefactorend::Number, numiters::Integer, fps::Integer, seconds::Number, height::Integer, aspectratio::Number, savename::AbstractString)
    #time = 0:1/fps:seconds
    #numsteps = length(time)
    width = Int(floor(aspectratio * height))

    escape = CUDA.fill(0.0, width, height)
    escape_color = CUDA.fill(ColorTypes.RGBA{Float32}(0.0, 0.0, 0.0), width, height)
    escape_cpu = zeros(Float64, width, height)
    escape_cpu_color = [ColorTypes.RGBA(Float32(0.5), Float32(0.5), Float32(0.5)) for x in 1:width, y in 1:height]
    grayscalecolors = [ColorTypes.RGBA(i/255, i/255, i/255) for i in 0:255]
    for (l, z) in enumerate(escape_cpu)
        escape_cpu_color[l] = grayscalecolors[100]
    end
    #@show typeof(escape_cpu)

    numtype = Float64
    local scalefactor = numtype(copy(scalefactorstart))
    local centerxlocal = numtype(copy(centerx))
    local centerylocal = numtype(copy(centery))
    local rotation = numtype(0.0)

    local wdown = false
    local adown = false
    local sdown = false
    local ddown = false
    local jdown = false
    local kdown = false
    local udown = false
    local idown = false

    local endscene = false

    movementscale = numtype(0.2)
    zoomscale = numtype(1.0)
    rotscale = numtype(1.0)

    Makie.lift(image.events.keyboardbuttons) do but
        #@show but

        wdown = Makie.Keyboard.w in but
        adown = Makie.Keyboard.a in but
        sdown = Makie.Keyboard.s in but
        ddown = Makie.Keyboard.d in but

        jdown = Makie.Keyboard.j in but
        kdown = Makie.Keyboard.k in but

        udown = Makie.Keyboard.u in but
        idown = Makie.Keyboard.i in but

        endscene = Makie.Keyboard.left_control in but

        nothing
    end

    N = width * height
    #numblocks = ceil(Int, N/32)
    numthreads = (16, 16)
    numblocks = ceil(Int, width / numthreads[1]), ceil(Int, height / numthreads[2])
    #@show numblocks
    i = 0
    listener = image.plots[1][:image].listeners[1]
    makimg = image.plots[1][:image]
    a = Complex{numtype}(1., 0.)
    b = Complex{numtype}(0., 0.)
    c = Complex{numtype}(0., 0.)
    d = Complex{numtype}(1., 0.)

    mobiusparams = [a, b, c, d]

    z1 = Complex{numtype}(-1.4, 0.)
    z2 = Complex{numtype}(1., -1.)
    lr = 0.02

    lambdaangle = numtype(0.5)
    adisk = Complex{numtype}(numtype(0.2), numtype(0.3))
    lambdaangle, adisk = randomdiskmobius()

    lambdaangletarget, adisktarget = randomdiskmobius()

    #a, b, c, d = diskmobius(lambdaangle, adisk)

    starttime = time()
    local prevtime = starttime
    local mobiuswalktime = prevtime

    l = 0
    while !endscene
    #for k in 1:1
    #Makie.record(image, savename, enumerate(time); framerate=60) do (i, t)
        CUDA.fill!(escape, 0.0)

        curtime = time()
        deltatime = curtime - prevtime
        curtimeperframe = (curtime - starttime) / i
        #@show curtimeperframe
        curfps = 1 / curtimeperframe
        if i % 20 == 0
            @show curfps
        end

        expscale = numtype(2.0 ^ -scalefactor)
        yrangeextent = numtype(expscale)
        xrangeextent = numtype(expscale * aspectratio)
        xstart = numtype(-xrangeextent / 2)
        ystart = numtype(-yrangeextent / 2)

        cosrot = numtype(cos(rotation))
        sinrot = numtype(sin(rotation))

        if wdown
            centerxlocal += -sinrot * movementscale * expscale * deltatime
            centerylocal += cosrot * movementscale * expscale * deltatime
            #@show centerxlocal
            #@show centerylocal
        end
        if adown
            centerxlocal -= cosrot * movementscale * expscale * deltatime
            centerylocal -= sinrot * movementscale * expscale * deltatime
            #@show centerxlocal
            #@show centerylocal
        end
        if sdown
            centerxlocal -= -sinrot * movementscale * expscale * deltatime
            centerylocal -= cosrot * movementscale * expscale * deltatime
            #@show centerxlocal
            #@show centerylocal
        end
        if ddown
            centerxlocal += cosrot * movementscale * expscale * deltatime
            centerylocal += sinrot * movementscale * expscale * deltatime
            #@show centerxlocal
            #@show centerylocal
        end

        # note: scale factor for Float32 is approximately capped at 32
        if jdown
            scalefactor -= zoomscale * deltatime
            expscale = numtype(2.0 ^ -scalefactor)
            yrangeextent = numtype(expscale)
            xrangeextent = numtype(expscale * aspectratio)
            xstart = numtype(-xrangeextent / 2)
            ystart = numtype(-yrangeextent / 2)
            #@show scalefactor
        end
        if kdown
            scalefactor += zoomscale * deltatime
            expscale = numtype(2.0 ^ -scalefactor)
            yrangeextent = numtype(expscale)
            xrangeextent = numtype(expscale * aspectratio)
            xstart = numtype(-xrangeextent / 2)
            ystart = numtype(-yrangeextent / 2)
            #@show scalefactor
        end

        if udown
            rotation -= rotscale * deltatime
            cosrot = numtype(cos(rotation))
            sinrot = numtype(sin(rotation))
            #@show rotation
        end
        if idown
            rotation += rotscale * deltatime
            cosrot = numtype(cos(rotation))
            sinrot = numtype(sin(rotation))
            #@show rotation
        end

        if curtime - mobiuswalktime > 3
            lambdaangletarget, adisktarget = randomdiskmobius()
            mobiuswalktime = curtime
        end

        lambdaangle += (lambdaangletarget - lambdaangle) * deltatime
        adisk += (adisktarget - adisk) * deltatime
        #a, b, c, d = diskmobius(lambdaangle, adisk)

        #@show "optimization"
        begin
            distfunc = mobiusparams -> complexdist(mobiustransform(z1, mobiusparams...), z2)
            distparams = distfunc(mobiusparams)
            gradmobius = gradient(distfunc, mobiusparams)[1]
            mobiusparams .-= gradmobius .* lr .* distparams
        end
        if distparams < 1e-3 || l > 100
            z2 = Complex{numtype}(randn(), randn())
            l = 0
        end
        l += 1
        if i % 5 == 0
            #@show distfunc(mobiusparams)
            #@show z2
            #@show mobiusparams
            #@show gradmobius
        end

        #@show "cuda"
        #@time CUDA.@sync begin
        CUDA.@sync begin
            #@cuda threads=numthreads blocks=numblocks mandelbrotandregiongpu!(escape, centerxlocal, xstart, xrangeextent, centerylocal, ystart, yrangeextent, numiters, height, width, cosrot, sinrot, a, b, c, d)
            @cuda threads=numthreads blocks=numblocks mandelbrotandregiongpu!(escape_color, centerxlocal, xstart, xrangeextent, centerylocal, ystart, yrangeextent, numiters, height, width, cosrot, sinrot, mobiusparams...)
        end

        #@show "copy"
        #@time CUDA.copyto!(escape_cpu, escape)
        CUDA.copyto!(makimg[], escape_color)

        #@show "copytocolor"
        #for (l, z) in enumerate(escape_cpu)
            #escape_cpu_color[l] = ColorTypes.RGBA{Float32}(Float32(z), Float32(z), Float32(z))
        #end

        #@show "imageplotcopy"
        #makimg[] = escape_cpu_color
        Observables.notify!(makimg)
        #@time makimg.listeners[1](makimg.val)

        i += 1
        sleep(0.0001)
        prevtime = curtime
        #end
    end
    #@show length(time)
    #return escape
    return nothing
end

function mandelbrotandregiongpu!(escape_color, centerx::N, xstart::N, xrangeextent::N, centery::N, ystart::N, yrangeextent::N,
        numiters, height, width, cosrot::N, sinrot::N, a::Complex{N}, b::Complex{N}, c::Complex{N}, d::Complex{N}) where {N <: Real}
    indexx = (blockIdx().x - 1) * blockDim().x + threadIdx().x
    indexy = (blockIdx().y - 1) * blockDim().y + threadIdx().y
    stridex = blockDim().x * gridDim().x
    stridey = blockDim().y * gridDim().y
    #@cuprintln("thread $indexx, $indexy, stride $stridex, $stridey, weight $width $height")
    for i in indexx:stridex:width, j in indexy:stridey:height
        prex_ij = (i - 1) / (width - 1) * xrangeextent + xstart
        prey_ij = (j - 1) / (height - 1) * yrangeextent + ystart
        z_ij = Complex{N}(cosrot * prex_ij - sinrot * prey_ij + centerx, sinrot * prex_ij + cosrot * prey_ij + centery)
        z_ij_mobius = mobiustransform(z_ij, a, b, c, d)
        escape_ij = mandelbrot(z_ij_mobius, numiters) / numiters
        escape_ij_f32 = Float32(escape_ij)
        @inbounds escape_color[i, j] = ColorTypes.RGBA{Float32}(escape_ij_f32, escape_ij_f32, escape_ij_f32)
    end
    return nothing
end

function mandelbrotandregiongpurational!(escape_color, centerx::T, xstart::T, xrangeextent::T, centery::T, ystart::T, yrangeextent::T,
        numiters, height, width, cosrot::T, sinrot::T, rational::ComplexRational{T, N}) where {T <: Real, N}
    indexx = (blockIdx().x - 1) * blockDim().x + threadIdx().x
    indexy = (blockIdx().y - 1) * blockDim().y + threadIdx().y
    stridex = blockDim().x * gridDim().x
    stridey = blockDim().y * gridDim().y
    #@cuprintln("thread $indexx, $indexy, stride $stridex, $stridey, weight $width $height")
    for i in indexx:stridex:width, j in indexy:stridey:height
        prex_ij = (i - 1) / (width - 1) * xrangeextent + xstart
        prey_ij = (j - 1) / (height - 1) * yrangeextent + ystart
        z_ij = Complex{T}(cosrot * prex_ij - sinrot * prey_ij + centerx, sinrot * prex_ij + cosrot * prey_ij + centery)
        z_ij_rational = transform(rational, z_ij)
        escape_ij = mandelbrot(z_ij_rational, numiters) / numiters
        escape_ij_f32 = Float32(escape_ij)
        @inbounds escape_color[i, j] = ColorTypes.RGBA{Float32}(escape_ij_f32, escape_ij_f32, escape_ij_f32)
    end
    return nothing
end

function mobiustransform(z::Complex{N}, a::Complex{N}, b::Complex{N}, c::Complex{N}, d::Complex{N})::Complex{N} where {N <: Real}
    (z * a + b) / (z * c + d)
end

function complexdist(z1::Complex{N}, z2::Complex{N})::N where {N <: Real}
    abs(z1 - z2)
end

function complexdistsq(z1::Complex{N}, z2::Complex{N})::N where {N <: Real}
    real(((z1 - z2) * conj(z1 - z2)))
end

function diskmobius(lambdaangle::N, a::Complex{N}) where {N <: Real}
    # maps the unit disk to itself
    # example 4 in section 3.4 in Complex Variables 2nd Ed by Stephen Fisher book
    # f(z) = lambda * (a - z) / (1 - āz)
    lambda = Complex{N}(cos(lambdaangle), sin(lambdaangle))
    moba = -lambda
    mobb = lambda * a
    mobc = -conj(a)
    mobd = Complex{N}(1, 0)
    return moba, mobb, mobc, mobd
end

function randomdiskmobius()
    lambdaangle = rand(Float64) * τ
    aangle = rand(Float64) * τ
    aradius = rand(Float64)
    a = Complex{Float64}(aradius * cos(aangle), aradius * sin(aangle))
    return lambdaangle, a
end

function gpu_add3!(y, x)
    index = (blockIdx().x - 1) * blockDim().x + threadIdx().x
    stride = blockDim().x * gridDim().x
    for i in index:stride:length(y)
        @inbounds y[i] += x[i]
    end
    return nothing
end

function gpu_rational!(y, x, ct)
    index = (blockIdx().x - 1) * blockDim().x + threadIdx().x
    stride = blockDim().x * gridDim().x
    #numpol = @SVector ComplexF64[1 + 0.2im, 1, 0.2 - 0.25im]
    for i in index:stride:length(y)
        #@inbounds y[i] = evalpoly(x[i], ct.numerator) / evalpoly(x[i], ct.denominator)
        @inbounds y[i] = transform(ct, x[i])
    end
    return nothing
end

function cudastuff()
    N = 1024
    x = CUDA.fill(ComplexF64(1.3, 0), N)
    y = CUDA.fill(ComplexF64(0, 0), N)
    numblocks = ceil(Int, N/256)
    numpol = @SVector ComplexF64[1, 1 + 0.1im]
    denpol = @SVector ComplexF64[2, 0]
    rat = ComplexRational(numpol, denpol)
    @cuda threads=256 blocks=numblocks gpu_rational!(y, x, rat)
    @show y
end

function bench_gpu3!(y, x)
    numblocks = ceil(Int, length(y)/256)
    CUDA.@sync begin
        @cuda threads=256 blocks=numblocks gpu_add3!(y, x)
    end
end

function gpu_mandel!(escape, z, numiters)
    index = (blockIdx().x - 1) * blockDim().x + threadIdx().x
    stride = blockDim().x * gridDim().x
    for i in index:stride:length(z)
        @inbounds escape[i] += mandelbrot(z[i], numiters)
    end
    return
end

function cudastuffmandel()
    centerx = -0.5609882
    centery = 0.6409865
    scalefactor = 1.0
    height = 360
    numiters = 200
    aspectratio = 16/9
    width = Int(floor(aspectratio * height))

    xrange, yrange = generatezoomranges(centerx, centery, scalefactor, height, aspectratio)
    z = CUDA.CuArray([Complex{Float64}(x, y) for x in xrange, y in yrange])
    escape = CUDA.fill(0., width, height)

    N = width * height
    numblocks = ceil(Int, N/256)
    #@cuda threads=256 blocks=numblocks gpu_mandel!(escape, z, numiters)
    #@show escape
    function bench_mandelgpu!(escapecur, zcur)
        CUDA.@sync begin
            @cuda threads=256 blocks=numblocks gpu_mandel!(escapecur, zcur, numiters)
        end
    end
    return escape, z, bench_mandelgpu!

end

ext: jl

max_stars_repo_name: zoemcc/fractrixzeta

max_issues_repo_name: zoemcc/fractrixzeta

max_forks_repo_name: zoemcc/fractrixzeta


text:
using Pkg

Pkg.add("PyCall")
Pkg.add("POMDPs")
Pkg.add("QuickPOMDPs")
Pkg.add("POMDPSimulators")
Pkg.add("POMDPPolicies")

using POMDPs
POMDPs.add_registry()
Pkg.add("QMDP")

using QuickPOMDPs
using PyCall

max_stars_repo_name: logankilpatrick/quickpomdps

max_issues_repo_name: logankilpatrick/quickpomdps

max_forks_repo_name: logankilpatrick/quickpomdps


text:
using Compat
using Test
using UUIDs: UUID, uuid1, uuid_version

@test isempty(detect_ambiguities(Base, Core, Compat))

@testset "CartesianIndex" begin
    # https://github.com/JuliaLang/julia/pull/29440
    ci = CartesianIndex(1, 1)
    @test length(-ci:ci) == 9
    # https://github.com/JuliaLang/julia/pull/29442
    @test oneunit(ci) === ci
    # https://github.com/JuliaLang/julia/pull/30268
    A = randn(1,2,3)
    @test get(A, CartesianIndex(1,2,3), :some_default) === A[1,2,3]
    @test get(A, CartesianIndex(2,2,3), :some_default) === :some_default
    @test get(11:15, CartesianIndex(6), nothing) === nothing
    @test get(11:15, CartesianIndex(5), nothing) === 15
end

# julia#29679
@test !isnothing(1)
@test isnothing(nothing)

# https://github.com/JuliaLang/julia/pull/29749
@testset "row/column/slice iterators" begin
    # Simple ones
    M = [1 2 3; 4 5 6; 7 8 9]
    @test collect(eachrow(M)) == collect(eachslice(M, dims = 1)) == [[1, 2, 3], [4, 5, 6], [7, 8, 9]]
    @test collect(eachcol(M)) == collect(eachslice(M, dims = 2)) == [[1, 4, 7], [2, 5, 8], [3, 6, 9]]
    @test_throws DimensionMismatch eachslice(M, dims = 4)

    # Higher-dimensional case
    M = reshape([(1:16)...], 2, 2, 2, 2)
    @test_throws MethodError collect(eachrow(M))
    @test_throws MethodError collect(eachcol(M))
    @test collect(eachslice(M, dims = 1))[1][:, :, 1] == [1 5; 3 7]
end

# Support for positional `stop`
@test range(0, 5, length = 6) == 0.0:1.0:5.0
@test range(0, 10, step = 2) == 0:2:10
@test_throws ArgumentError range(0, 10)

mutable struct TLayout
    x::Int8
    y::Int16
    z::Int32
end
tlayout = TLayout(5,7,11)
@test hasfield(TLayout, :y)
@test !hasfield(TLayout, :a)
@test hasproperty(tlayout, :x)
@test !hasproperty(tlayout, :p)

@test merge((a=1,b=1)) == (a=1,b=1)
@test merge((a=1,), (b=2,), (c=3,)) == (a=1,b=2,c=3)

@testset "only" begin
    @test only([3]) === 3
    @test_throws ArgumentError only([])
    @test_throws ArgumentError only([3, 2])

    @test @inferred(only((3,))) === 3
    @test_throws ArgumentError only(())
    @test_throws ArgumentError only((3, 2))

    @test only(Dict(1=>3)) === (1=>3)
    @test_throws ArgumentError only(Dict{Int,Int}())
    @test_throws ArgumentError only(Dict(1=>3, 2=>2))

    @test only(Set([3])) === 3
    @test_throws ArgumentError only(Set(Int[]))
    @test_throws ArgumentError only(Set([3,2]))

    @test @inferred(only((;a=1))) === 1
    @test_throws ArgumentError only(NamedTuple())
    @test_throws ArgumentError only((a=3, b=2.0))

    @test @inferred(only(1)) === 1
    @test @inferred(only('a')) === 'a'
    @test @inferred(only(Ref([1, 2]))) == [1, 2]
    @test_throws ArgumentError only(Pair(10, 20))

    @test only(1 for ii in 1:1) === 1
    @test only(1 for ii in 1:10 if ii < 2) === 1
    @test_throws ArgumentError only(1 for ii in 1:10)
    @test_throws ArgumentError only(1 for ii in 1:10 if ii > 2)
    @test_throws ArgumentError only(1 for ii in 1:10 if ii > 200)
end

# https://github.com/JuliaLang/julia/pull/32628
@testset "mod with ranges" begin
    for n in -10:10
        @test mod(n, 0:4) == mod(n, 5)
        @test mod(n, 1:5) == mod1(n, 5)
        @test mod(n, 2:6) == 2 + mod(n-2, 5)
        @test mod(n, Base.OneTo(5)) == mod1(n, 5)
    end
    @test mod(Int32(3), 1:5) == 3
    @test mod(big(typemax(Int))+99, 0:4) == mod(big(typemax(Int))+99, 5)
    @test_throws MethodError mod(3.141, 1:5)
    @test_throws MethodError mod(3, UnitRange(1.0,5.0))
    @test_throws MethodError mod(3, 1:2:7)
    @test_throws DivideError mod(3, 1:0)
end

# https://github.com/JuliaLang/julia/pull/31664
@testset "search character in strings" begin
    let astr = "Hello, world.\n",
        u8str = "∀ ε > 0, ∃ δ > 0: |x-y| < δ ⇒ |f(x)-f(y)| < ε"
        @test_throws BoundsError findnext('z', astr, 0)
        @test_throws BoundsError findnext('∀', astr, 0)
        @test findfirst('x', astr) == nothing
        @test findfirst('\0', astr) == nothing
        @test findfirst('\u80', astr) == nothing
        @test findfirst('∀', astr) == nothing
        @test findfirst('∀', u8str) == 1
        @test findfirst('ε', u8str) == 5
        @test findfirst('H', astr) == 1
        @test findfirst('l', astr) == 3
        @test findfirst('e', astr) == 2
        @test findfirst('u', astr) == nothing
        @test findnext('l', astr, 4) == 4
        @test findnext('l', astr, 5) == 11
        @test findnext('l', astr, 12) == nothing
        @test findfirst(',', astr) == 6
        @test findnext(',', astr, 7) == nothing
        @test findfirst('\n', astr) == 14
        @test findnext('\n', astr, 15) == nothing
        @test_throws BoundsError findnext('ε', astr, nextind(astr,lastindex(astr))+1)
        @test_throws BoundsError findnext('a', astr, nextind(astr,lastindex(astr))+1)
        @test findlast('x', astr) == nothing
        @test findlast('\0', astr) == nothing
        @test findlast('\u80', astr) == nothing
        @test findlast('∀', astr) == nothing
        @test findlast('∀', u8str) == 1
        @test findlast('ε', u8str) == 54
        @test findlast('H', astr) == 1
        @test findprev('H', astr, 0) == nothing
        @test findlast('l', astr) == 11
        @test findprev('l', astr, 5) == 4
        @test findprev('l', astr, 4) == 4
        @test findprev('l', astr, 3) == 3
        @test findprev('l', astr, 2) == nothing
        @test findlast(',', astr) == 6
        @test findprev(',', astr, 5) == nothing
        @test findlast('\n', astr) == 14
    end
end

using LinearAlgebra, Random

@testset "generalized dot #32739" begin
    Random.seed!(42) # https://github.com/JuliaLang/Compat.jl/issues/712

    # stdlib/LinearAlgebra/test/generic.jl
    for elty in (Int, Float32, Float64, BigFloat, Complex{Float32}, Complex{Float64}, Complex{BigFloat})
        n = 10
        if elty <: Int
            A = rand(-n:n, n, n)
            x = rand(-n:n, n)
            y = rand(-n:n, n)
        elseif elty <: Real
            A = convert(Matrix{elty}, randn(n,n))
            x = rand(elty, n)
            y = rand(elty, n)
        else
            A = convert(Matrix{elty}, complex.(randn(n,n), randn(n,n)))
            x = rand(elty, n)
            y = rand(elty, n)
        end
        @test dot(x, A, y) ≈ dot(A'x, y) ≈ *(x', A, y) ≈ (x'A)*y
        @test dot(x, A', y) ≈ dot(A*x, y) ≈ *(x', A', y) ≈ (x'A')*y
        elty <: Real && @test dot(x, transpose(A), y) ≈ dot(x, transpose(A)*y) ≈ *(x', transpose(A), y) ≈ (x'*transpose(A))*y
        B = reshape([A], 1, 1)
        x = [x]
        y = [y]
        @test dot(x, B, y) ≈ dot(B'x, y)
        @test dot(x, B', y) ≈ dot(B*x, y)
        elty <: Real && @test dot(x, transpose(B), y) ≈ dot(x, transpose(B)*y)
    end

    # stdlib/LinearAlgebra/test/symmetric.jl
    n = 10
    areal = randn(n,n)/2
    aimg  = randn(n,n)/2
    @testset for eltya in (Float32, Float64, ComplexF32, ComplexF64, BigFloat, Int)
        a = eltya == Int ? rand(1:7, n, n) : convert(Matrix{eltya}, eltya <: Complex ? complex.(areal, aimg) : areal)
        asym = transpose(a) + a                 # symmetric indefinite
        aherm = a' + a                 # Hermitian indefinite
        apos  = a' * a                 # Hermitian positive definite
        aposs = apos + transpose(apos)        # Symmetric positive definite
        ε = εa = eps(abs(float(one(eltya))))
        x = randn(n)
        y = randn(n)
        b = randn(n,n)/2
        x = eltya == Int ? rand(1:7, n) : convert(Vector{eltya}, eltya <: Complex ? complex.(x, zeros(n)) : x)
        y = eltya == Int ? rand(1:7, n) : convert(Vector{eltya}, eltya <: Complex ? complex.(y, zeros(n)) : y)
        b = eltya == Int ? rand(1:7, n, n) : convert(Matrix{eltya}, eltya <: Complex ? complex.(b, zeros(n,n)) : b)

        @testset "generalized dot product" begin
            for uplo in (:U, :L)
                @test dot(x, Hermitian(aherm, uplo), y) ≈ dot(x, Hermitian(aherm, uplo)*y) ≈ dot(x, Matrix(Hermitian(aherm, uplo)), y)
                @test dot(x, Hermitian(aherm, uplo), x) ≈ dot(x, Hermitian(aherm, uplo)*x) ≈ dot(x, Matrix(Hermitian(aherm, uplo)), x)
            end
            if eltya <: Real
                for uplo in (:U, :L)
                    @test dot(x, Symmetric(aherm, uplo), y) ≈ dot(x, Symmetric(aherm, uplo)*y) ≈ dot(x, Matrix(Symmetric(aherm, uplo)), y)
                    @test dot(x, Symmetric(aherm, uplo), x) ≈ dot(x, Symmetric(aherm, uplo)*x) ≈ dot(x, Matrix(Symmetric(aherm, uplo)), x)
                end
            end
        end
    end

    # stdlib/LinearAlgebra/test/uniformscaling.jl
    @testset "generalized dot" begin
        x = rand(-10:10, 3)
        y = rand(-10:10, 3)
        λ = rand(-10:10)
        J = UniformScaling(λ)
        @test dot(x, J, y) == λ*dot(x, y)
    end

    # stdlib/LinearAlgebra/test/bidiag.jl
    # The special method for this is not in Compat #683, so this tests the generic fallback
    @testset "generalized dot" begin
        for elty in (Float64, ComplexF64)
            dv = randn(elty, 5)
            ev = randn(elty, 4)
            x = randn(elty, 5)
            y = randn(elty, 5)
            for uplo in (:U, :L)
                B = Bidiagonal(dv, ev, uplo)
                @test dot(x, B, y) ≈ dot(B'x, y) ≈ dot(x, Matrix(B), y)
            end
        end
    end

    # Diagonal -- no such test in Base.
    @testset "diagonal" begin
        x = rand(-10:10, 3) .+ im
        y = rand(-10:10, 3) .+ im
        d = Diagonal(rand(-10:10, 3) .+ im)
        @test dot(x,d,y) == dot(x,collect(d),y) == dot(x, d*y)
    end
end

# https://github.com/JuliaLang/julia/pull/33568
@testset "function composition" begin
    @test ∘(x -> x-2, x -> x-3, x -> x+5)(7) == 7
    fs = [x -> x[1:2], uppercase, lowercase]
    @test ∘(fs...)("ABC") == "AB"

    # https://github.com/JuliaLang/julia/pull/34251
    @testset "unary" begin
        @test ∘(identity) === identity
        @test ∘(inv) === inv
    end
end

# https://github.com/JuliaLang/julia/pull/33128
@testset "pkgdir" begin
    @test pkgdir(Main) === nothing
    @test joinpath(pkgdir(Compat), "") == abspath(joinpath(@__DIR__, ".."))
end

# https://github.com/JuliaLang/julia/pull/33736/
@testset "ReverseOrdering constructor" begin
    @test Base.Order.ReverseOrdering() == Base.Order.Reverse
end

# https://github.com/JuliaLang/julia/pull/32968
@testset "filter on Tuples" begin
    @test filter(isodd, (1,2,3)) == (1, 3)
    @test filter(isequal(2), (true, 2.0, 3)) === (2.0,)
    @test filter(i -> true, ()) == ()
    @test filter(identity, (true,)) === (true,)
    longtuple = ntuple(identity, 20)
    @test filter(iseven, longtuple) == ntuple(i->2i, 10)
    @test filter(x -> x<2, (longtuple..., 1.5)) === (1, 1.5)
end

# https://github.com/JuliaLang/julia/pull/34652
@testset "ismutable" begin
    @test ismutable(1) == false
    @test ismutable([]) == true
end

# https://github.com/JuliaLang/julia/pull/28761
@testset "uuid5" begin
    u1 = uuid1()
    u5 = uuid5(u1, "julia")
    @test uuid_version(u5) == 5
    @test u5 == UUID(string(u5)) == UUID(GenericString(string(u5)))
    @test u5 == UUID(UInt128(u5))

    following_uuids = [
        UUID("22b4a8a1-e548-4eeb-9270-60426d66a48e"),
        UUID("30ea6cfd-c270-569f-b4cb-795dead63686"),
        UUID("31099374-e3a0-5fde-9482-791c639bf29b"),
        UUID("6b34b357-a348-53aa-8c71-fb9b06c3a51e"),
        UUID("fdbd7d4d-c462-59cc-ae6a-0c3b010240e2"),
        UUID("d8cc6298-75d5-57e0-996c-279259ab365c"),
    ]

    for (idx, init_uuid) in enumerate(following_uuids[1:end-1])
        next_id = uuid5(init_uuid, "julia")
        @test next_id == following_uuids[idx+1]
    end

    # Some UUID namespaces provided in the appendix of RFC 4122
    # https://tools.ietf.org/html/rfc4122.html#appendix-C
    namespace_dns  = UUID(0x6ba7b8109dad11d180b400c04fd430c8) # 6ba7b810-9dad-11d1-80b4-00c04fd430c8
    namespace_url  = UUID(0x6ba7b8119dad11d180b400c04fd430c8) # 6ba7b811-9dad-11d1-80b4-00c04fd430c8
    namespace_oid  = UUID(0x6ba7b8129dad11d180b400c04fd430c8) # 6ba7b812-9dad-11d1-80b4-00c04fd430c8
    namespace_x500 = UUID(0x6ba7b8149dad11d180b400c04fd430c8) # 6ba7b814-9dad-11d1-80b4-00c04fd430c8

    # Python-generated UUID following each of the standard namespaces
    standard_namespace_uuids = [
        (namespace_dns,  UUID("00ca23ad-40ef-500c-a910-157de3950d07")),
        (namespace_oid,  UUID("b7bf72b0-fb4e-538b-952a-3be296f07f6d")),
        (namespace_url,  UUID("997cd5be-4705-5439-9fe6-d77b18d612e5")),
        (namespace_x500, UUID("993c6684-82e7-5cdb-bd46-9bff0362e6a9")),
    ]

    for (init_uuid, next_uuid) in standard_namespace_uuids
        result = uuid5(init_uuid, "julia")
        @test next_uuid == result
    end
end

@testset "Irrational zero and one" begin
    @test one(pi) === true
    @test zero(pi) === false
    @test one(typeof(pi)) === true
    @test zero(typeof(pi)) === false
end

# https://github.com/JuliaLang/julia/pull/32753
@testset "evalpoly real" begin
    for x in -1.0:2.0, p1 in -3.0:3.0, p2 in -3.0:3.0, p3 in -3.0:3.0
        evpm = @evalpoly(x, p1, p2, p3)
        @test evalpoly(x, (p1, p2, p3)) == evpm
        @test evalpoly(x, [p1, p2, p3]) == evpm
    end
end
@testset "evalpoly complex" begin
    for x in -1.0:2.0, y in -1.0:2.0, p1 in -3.0:3.0, p2 in -3.0:3.0, p3 in -3.0:3.0
        z = x + im * y
        evpm = @evalpoly(z, p1, p2, p3)
        @test evalpoly(z, (p1, p2, p3)) == evpm
        @test evalpoly(z, [p1, p2, p3]) == evpm
    end
    @test evalpoly(1+im, (2,)) == 2
    @test evalpoly(1+im, [2,]) == 2
end

# https://github.com/JuliaLang/julia/pull/35298
begin
    # A custom linear slow sparse-like array that relies upon Dict for its storage
    struct TSlow{T,N} <: AbstractArray{T,N}
        data::Dict{NTuple{N,Int}, T}
        dims::NTuple{N,Int}
    end
    TSlow(::Type{T}, dims::Int...) where {T} = TSlow(T, dims)
    TSlow(::Type{T}, dims::NTuple{N,Int}) where {T,N} = TSlow{T,N}(Dict{NTuple{N,Int}, T}(), dims)

    TSlow{T,N}(X::TSlow{T,N})         where {T,N  } = X
    TSlow(     X::AbstractArray{T,N}) where {T,N  } = TSlow{T,N}(X)
    TSlow{T  }(X::AbstractArray{_,N}) where {T,N,_} = TSlow{T,N}(X)
    TSlow{T,N}(X::AbstractArray     ) where {T,N  } = begin
        A = TSlow(T, size(X))
        for I in CartesianIndices(X)
            A[Tuple(I)...] = X[Tuple(I)...]
        end
        A
    end
    Base.size(A::TSlow) = A.dims
    Base.similar(A::TSlow, ::Type{T}, dims::Dims) where {T} = TSlow(T, dims)
    Base.IndexStyle(::Type{A}) where {A<:TSlow} = IndexCartesian()
    Base.getindex(A::TSlow{T,N}, i::Vararg{Int,N}) where {T,N} = get(A.data, i, zero(T))
    Base.setindex!(A::TSlow{T,N}, v, i::Vararg{Int,N}) where {T,N} = (A.data[i] = v)
end

# https://github.com/JuliaLang/julia/pull/35304
@testset "similar(PermutedDimsArray)" begin
    x = PermutedDimsArray([1 2; 3 4], (2, 1))
    @test similar(x, 3,3) isa Array
    z = TSlow([1 2; 3 4])
    x_slow = PermutedDimsArray(z, (2, 1))
    @test similar(x_slow, 3,3) isa TSlow
end

# https://github.com/JuliaLang/julia/pull/34548
@testset "@NamedTuple" begin
    @test (@NamedTuple {a::Int, b::String}) === NamedTuple{(:a, :b),Tuple{Int,String}} ===
        @NamedTuple begin
            a::Int
            b::String
        end
    @test (@NamedTuple {a::Int, b}) === NamedTuple{(:a, :b),Tuple{Int,Any}}
end

struct NonFunctionCallable end
(::NonFunctionCallable)(args...) = +(args...)

@testset "mergewith" begin
    d1 = Dict("A" => 1, "B" => 2)
    d2 = Dict("B" => 3.0, "C" => 4.0)
    @test mergewith(+, d1, d2) == Dict("A" => 1, "B" => 5, "C" => 4)
    @test mergewith(*, d1, d2) == Dict("A" => 1, "B" => 6, "C" => 4)
    @test mergewith(-, d1, d2) == Dict("A" => 1, "B" => -1, "C" => 4)
    @test mergewith(NonFunctionCallable(), d1, d2) == Dict("A" => 1, "B" => 5, "C" => 4)
    @test foldl(mergewith(+), [d1, d2]; init=Dict{Union{},Union{}}()) ==
        Dict("A" => 1, "B" => 5, "C" => 4)
end

@testset "mergewith!" begin
    d1 = Dict("A" => 1, "B" => 3, "C" => 4)
    d2 = Dict("B" => 3, "C" => 4)
    mergewith!(+, d1, d2)
    @test d1 == Dict("A" => 1, "B" => 6, "C" => 8)
    mergewith!(*, d1, d2)
    @test d1 == Dict("A" => 1, "B" => 18, "C" => 32)
    mergewith!(-, d1, d2)
    @test d1 == Dict("A" => 1, "B" => 15, "C" => 28)
    mergewith!(NonFunctionCallable(), d1, d2)
    @test d1 == Dict("A" => 1, "B" => 18, "C" => 32)
    @test foldl(mergewith!(+), [d1, d2]; init=empty(d1)) ==
        Dict("A" => 1, "B" => 21, "C" => 36)
end

# https://github.com/JuliaLang/julia/pull/34427
@testset "isdisjoint" begin
    for S in (Set, BitSet, Vector)
        for (l,r) in ((S([1,2]),     S([3,4])),
                      (S([5,6,7,8]), S([7,8,9])),
                      (S([1,2]),     S([3,4])),
                      (S([5,6,7,8]), S([7,8,9])),
                      (S([1,2,3]),   S()),
                      (S(),          S()),
                      (S(),          S([1,2,3])),
                      (S([1,2,3]),   S([1])),
                      (S([1,2,3]),   S([1,2])),
                      (S([1,2,3]),   S([1,2,3])),
                      (S([1,2,3]),   S([4])),
                      (S([1,2,3]),   S([4,1])))
            @test isdisjoint(l,l) == isempty(l)
            @test isdisjoint(l,r) == isempty(intersect(l,r))
        end
    end
end

# https://github.com/JuliaLang/julia/pull/35577
@testset "union on OneTo" begin
    @test union(Base.OneTo(3), Base.OneTo(4)) === Base.OneTo(4)
end

# https://github.com/JuliaLang/julia/pull/35929
# https://github.com/JuliaLang/julia/pull/29135
@testset "strided transposes" begin
    for t in (Adjoint, Transpose)
        @test strides(t(rand(3))) == (3, 1)
        @test strides(t(rand(3,2))) == (3, 1)
        @test strides(t(view(rand(3, 2), :))) == (6, 1)
        @test strides(t(view(rand(3, 2), :, 1:2))) == (3, 1)

         A = rand(3)
        @test pointer(t(A)) === pointer(A)
        B = rand(3,1)
        @test pointer(t(B)) === pointer(B)
    end
    @test_throws MethodError strides(Adjoint(rand(3) .+ rand(3).*im))
    @test_throws MethodError strides(Adjoint(rand(3, 2) .+ rand(3, 2).*im))
    @test strides(Transpose(rand(3) .+ rand(3).*im)) == (3, 1)
    @test strides(Transpose(rand(3, 2) .+ rand(3, 2).*im)) == (3, 1)

     C = rand(3) .+ rand(3).*im
    @test_throws ErrorException pointer(Adjoint(C))
    @test pointer(Transpose(C)) === pointer(C)
    D = rand(3,2) .+ rand(3,2).*im
    @test_throws ErrorException pointer(Adjoint(D))
    @test pointer(Transpose(D)) === pointer(D)
end

# https://github.com/JuliaLang/julia/pull/27516
@testset "two arg @inferred" begin
    g(a) = a < 10 ? missing : 1
    @test ismissing(g(9))
    @test g(10) == 1
    @inferred Missing g(9)
    @inferred Missing g(10)
end

# https://github.com/JuliaLang/julia/pull/36360
@testset "get_set_num_threads" begin
    default = Compat.get_num_threads()
    @test default isa Int # seems dodgy, could be nothing!
    @test default > 0
    Compat.set_num_threads(1)
    @test Compat.get_num_threads() === 1
    Compat.set_num_threads(default)
    @test Compat.get_num_threads() === default

    # Run the ::Nothing method, to check no error:
    Compat.set_num_threads(nothing)
    Compat.set_num_threads(default)

    if VERSION < v"1.6.0-DEV.322"
        # These tests from PR rely on internal functions which would be BLAS. not Compat.
        @test_logs (:warn,) match_mode=:any Compat._set_num_threads(1, _blas=:unknown)
        if Compat.guess_vendor() !== :osxblas
            # test osxblas which is not covered by CI
            withenv("VECLIB_MAXIMUM_THREADS" => nothing) do
                @test @test_logs(
                    (:warn,),
                    (:warn,),
                    match_mode=:any,
                    Compat._get_num_threads(_blas=:osxblas),
                ) === nothing
                @test_logs Compat._set_num_threads(1, _blas=:osxblas)
                @test @test_logs(Compat._get_num_threads(_blas=:osxblas)) === 1
                @test_logs Compat._set_num_threads(2, _blas=:osxblas)
                @test @test_logs(Compat._get_num_threads(_blas=:osxblas)) === 2
            end
        end
    end
end

# https://github.com/JuliaLang/julia/pull/30915
@testset "curried comparisons" begin
    eql5 = (==)(5)
    neq5 = (!=)(5)
    gte5 = (>=)(5)
    lte5 = (<=)(5)
    gt5  = (>)(5)
    lt5  = (<)(5)

    @test eql5(5) && !eql5(0)
    @test neq5(6) && !neq5(5)
    @test gte5(5) && gte5(6)
    @test lte5(5) && lte5(4)
    @test gt5(6) && !gt5(5)
    @test lt5(4) && !lt5(5)
end

@testset "contains" begin
    @test contains("foo", "o")
    @test contains("o")("foo")
end

# https://github.com/JuliaLang/julia/pull/35052
@testset "curried startswith/endswith" begin
    @test startswith("a")("abcd")
    @test endswith("d")("abcd")
end

# https://github.com/JuliaLang/julia/pull/37517
@testset "ComposedFunction" begin
    @test sin ∘ cos isa Compat.ComposedFunction
    @test sin ∘ cos === Compat.ComposedFunction(sin, cos)
    c = sin ∘ cos
    @test c.outer === sin
    @test c.inner === cos
    if VERSION < v"1.6.0-DEV.1037"
        @test c.f === sin
        @test c.g === cos
        @test propertynames(c) == (:f, :g, :outer, :inner)
    else
        @test propertynames(c) == (:outer, :inner)
    end
end

# From spawn.jl
shcmd = `sh`
havebb = false
if Sys.iswindows()
    busybox = download("https://cache.julialang.org/https://frippery.org/files/busybox/busybox.exe", joinpath(tempdir(), "busybox.exe"))
    havebb = try # use busybox-w32 on windows, if available
        success(`$busybox`)
        true
    catch
        false
    end
    if havebb
        shcmd = `$busybox sh`
    end
end

# https://github.com/JuliaLang/julia/pull/37244
@testset "addenv()" begin
    cmd = Cmd(`$shcmd -c "echo \$FOO \$BAR"`, env=Dict("FOO" => "foo"))
    @test strip(String(read(cmd))) == "foo"
    cmd = addenv(cmd, "BAR" => "bar")
    @test strip(String(read(cmd))) == "foo bar"
    cmd = addenv(cmd, Dict("FOO" => "bar"))
    @test strip(String(read(cmd))) == "bar bar"
    cmd = addenv(cmd, ["FOO=baz"])
    @test strip(String(read(cmd))) == "baz bar"
end

# https://github.com/JuliaLang/julia/pull/37559
@testset "reinterpred(reshape, ...)" begin
    # simplified from PR
    Ar = Int64[1 3; 2 4]
    @test @inferred(ndims(reinterpret(reshape, Complex{Int64}, Ar))) == 1
    @test @inferred(axes(reinterpret(reshape, Complex{Int64}, Ar))) === (Base.OneTo(2),)
    @test @inferred(size(reinterpret(reshape, Complex{Int64}, Ar))) == (2,)

    _B = Complex{Int64}[5+6im, 7+8im, 9+10im]
    @test @inferred(ndims(reinterpret(reshape, Int64, _B))) == 2
    @test @inferred(axes(reinterpret(reshape, Int64, _B))) === (Base.OneTo(2), Base.OneTo(3))
    @test @inferred(size(reinterpret(reshape, Int64, _B))) == (2, 3)
    @test @inferred(ndims(reinterpret(reshape, Int128, _B))) == 1
    @test @inferred(axes(reinterpret(reshape, Int128, _B))) === (Base.OneTo(3),)
    @test @inferred(size(reinterpret(reshape, Int128, _B))) == (3,)

    A = Int64[1, 2, 3, 4]
    Av = [Int32[1,2], Int32[3,4]]

    @test_throws ArgumentError reinterpret(Vector{Int64}, A) # ("cannot reinterpret `Int64` as `Vector{Int64}`, type `Vector{Int64}` is not a bits type")
    @test_throws ArgumentError reinterpret(Int32, Av) # ("cannot reinterpret `Vector{Int32}` as `Int32`, type `Vector{Int32}` is not a bits type")
    @test_throws ArgumentError("cannot reinterpret a zero-dimensional `Int64` array to `Int32` which is of a different size") reinterpret(Int32, reshape([Int64(0)]))
    @test_throws ArgumentError("cannot reinterpret a zero-dimensional `Int32` array to `Int64` which is of a different size") reinterpret(Int64, reshape([Int32(0)]))
    @test_throws ArgumentError reinterpret(Tuple{Int,Int}, [1,2,3,4,5]) # ("""cannot reinterpret an `$Int` array to `Tuple{$Int, $Int}` whose first dimension has size `5`.
                               # The resulting array would have non-integral first dimension.
                               # """)
    @test_throws ArgumentError("`reinterpret(reshape, Complex{Int64}, a)` where `eltype(a)` is Int64 requires that `axes(a, 1)` (got Base.OneTo(4)) be equal to 1:2 (from the ratio of element sizes)") reinterpret(reshape, Complex{Int64}, A)
    @test_throws ArgumentError("`reinterpret(reshape, T, a)` requires that one of `sizeof(T)` (got 24) and `sizeof(eltype(a))` (got 16) be an integer multiple of the other") reinterpret(reshape, NTuple{3, Int64}, _B)
    @test_throws ArgumentError reinterpret(reshape, Vector{Int64}, Ar) # ("cannot reinterpret `Int64` as `Vector{Int64}`, type `Vector{Int64}` is not a bits type")
    @test_throws ArgumentError("cannot reinterpret a zero-dimensional `UInt8` array to `UInt16` which is of a larger size") reinterpret(reshape, UInt16, reshape([0x01]))

    # getindex
    _A = A
    @test reinterpret(Complex{Int64}, _A) == [1 + 2im, 3 + 4im]
    @test reinterpret(Float64, _A) == reinterpret.(Float64, A)
    @test reinterpret(reshape, Float64, _A) == reinterpret.(Float64, A)

    Ars = Ar
    @test reinterpret(reshape, Complex{Int64}, Ar) == [1 + 2im, 3 + 4im]
    @test reinterpret(reshape, Float64, Ar) == reinterpret.(Float64, Ars)

    # setindex
    A3 = collect(reshape(1:18, 2, 3, 3))
    A3r = reinterpret(reshape, Complex{Int}, A3)
    @test A3r[4] === A3r[1,2] === A3r[CartesianIndex(1, 2)] === 7+8im
    A3r[2,3] = -8-15im
    @test A3[1,2,3] == -8
    @test A3[2,2,3] == -15
    A3r[4] = 100+200im
    @test A3[1,1,2] == 100
    @test A3[2,1,2] == 200
    A3r[CartesianIndex(1,2)] = 300+400im
    @test A3[1,1,2] == 300
    @test A3[2,1,2] == 400

    # Test 0-dimensional Arrays
    A = zeros(UInt32)
    B = reinterpret(Int32,A)
    Brs = reinterpret(reshape,Int32,A)
    @test size(B) == size(Brs) == ()
    @test axes(B) == axes(Brs) == ()
    B[] = Int32(5)
    @test B[] === Int32(5)
    @test Brs[] === Int32(5)
    @test A[] === UInt32(5)

    # reductions
    a = [(1,2,3), (4,5,6)]
    ars = reinterpret(reshape, Int, a)
    @test sum(ars) == 21
    @test sum(ars; dims=1) == [6 15]
    @test sum(ars; dims=2) == reshape([5,7,9], (3, 1))
    @test sum(ars; dims=(1,2)) == reshape([21], (1, 1))
    # also test large sizes for the pairwise algorithm
    a = [(k,k+1,k+2) for k = 1:3:4000]
    ars = reinterpret(reshape, Int, a)
    @test sum(ars) == 8010003
end

# https://github.com/JuliaLang/julia/pull/29634
@testset "5-arg mul!" begin
    A = [1.0 2.0; 3.0 4.0]
    B = [10.0, 2.0]
    C = [100.0, 1000.0]
    x = -20.0
    alpha = 0.1
    beta = 0.1

    Cmut = copy(C)
    @test  Cmut == mul!(Cmut, A, B, alpha, beta) ≈ ((A * B * alpha) + (C * beta))

    Cmut = copy(C)
    @test  Cmut == mul!(Cmut, x, B, alpha, beta) ≈ ((x * B * alpha) + (C * beta))

    Cmut = copy(C)
    @test  Cmut == mul!(Cmut, B, x, alpha, beta) ≈ ((B * x * alpha) + (C * beta))
end

# https://github.com/JuliaLang/julia/pull/35243
@testset "parseatom and parseall" begin
    @test Compat.parseatom(raw"foo$(@bar)baz", 5; filename="foo") ==
        (Expr(:macrocall, Symbol("@bar"), LineNumberNode(1, :foo)), 11)

    ex = Compat.parseall(
        raw"""
        begin
            @a b
            @c() + 1
        end
        """;
        filename="foo",
    )
    @test ex == Expr(:toplevel,
        LineNumberNode(1, :foo),
        Expr(:block,
            LineNumberNode(2, :foo),
            Expr(:macrocall, Symbol("@a"), LineNumberNode(2, :foo), :b),
            LineNumberNode(3, :foo),
            Expr(:call,
                 :+,
                 Expr(:macrocall, Symbol("@c"), LineNumberNode(3, :foo)),
                 1,
            ),
        ),
    )

    ex = Compat.parseall(
        raw"""

        begin a = 1 end

        begin
            b = 2
        end
        """;
        filename="foo",
    )
    @test ex == Expr(:toplevel,
        LineNumberNode(2, :foo),
        Expr(:block,
            LineNumberNode(2, :foo),
            :(a = 1),
        ),
        LineNumberNode(4, :foo),
        Expr(:block,
            LineNumberNode(5, :foo),
            :(b = 2),
        ),
    )
end

include("iterators.jl")

nothing


text:
# indexing.jl - low-level utils for working with indexed expressions
# TODO: rename since indexing isn't in focus of this package anymore

# const IDX_NAMES = [:i, :j, :k, :m, :n, :p, :q, :r, :s, :l]

## utils

isref(ex::Expr) = ex.head == :ref
isref(x) = false

isindexed(ex::Expr) = isref(ex) || any(isindexed, ex.args)
isindexed(x) = false

"""
Split possibly indexed variable into a name and indices. Examples:

    split_indexed(:x)         ==> (:x, [])
    split_indexed(:(x[i,j]))  ==> (:x, [:i,:j])

See also: make_indexed
"""
split_indexed(var) = (var, [])

function split_indexed(ex::Expr)
    @assert ex.head == :ref
    return convert(Symbol, ex.args[1]), ex.args[2:end]
end

"""
Make indexed variable. Examples:

    make_indexed(:x, [])       ==> :x
    make_indexed(:x, [:i,:j])  ==> :(x[i,j])

See also: split_indexed
"""
function make_indexed(vname::Symbol, vidxs::Vector)
    return isempty(vidxs) ? vname : Expr(:ref, vname, vidxs...)
end

"""Generate index names and make an indexed variable using them"""
function with_indices(x::Symbol, start_idx::Int, num_idxs::Int)
    return make_indexed(x, IDX_NAMES[start_idx:start_idx+num_idxs-1])
end

with_indices(x::Symbol, num_idxs::Int) = with_indices(x, 1, num_idxs)

## get_vars, get_var_names, get_indices

function get_vars!(x, rec::Bool, result::Vector{Union{Symbol, Expr}})
    # do nothing
end

function get_vars!(x::Symbol, rec::Bool, result::Vector{Union{Symbol, Expr}})
    push!(result, x)
end

function get_vars!(ex::Expr, rec::Bool, result::Vector{Union{Symbol, Expr}})
    get_vars!(ExH(ex), rec, result)
end

function get_vars!(ex::ExH{:ref}, rec::Bool, result::Vector{Union{Symbol, Expr}})
    push!(result, Expr(ex))
end

function get_vars!(ex::ExH{:call}, rec::Bool, result::Vector{Union{Symbol, Expr}})
    for arg in ex.args[2:end]
        if isa(arg, Symbol) || isref(arg)
            push!(result, arg)
        elseif rec
            get_vars!(arg, rec, result)
        end
    end
end

function get_vars!(ex::ExH{:.}, rec::Bool, result::Vector{Union{Symbol, Expr}})
    @assert ex.args[2].head == :tuple "Dot is only allowed in broadcasting: $ex"
    for arg in ex.args[2].args
        if isa(arg, Symbol) || isref(arg)
            push!(result, arg)
        elseif rec
            get_vars!(arg, rec, result)
        end
    end
end

function get_vars!(ex::ExH{Symbol("'")}, rec::Bool, result::Vector{Union{Symbol, Expr}})
    arg = ex.args[1]
    if isa(arg, Symbol) || isref(arg)
        push!(result, arg)
    elseif rec
        get_vars!(arg, rec, result)
    end
end

function get_vars!(ex::ExH{:(=)}, rec::Bool, result::Vector{Union{Symbol, Expr}})
    get_vars!(ex.args[1], rec, result)
    get_vars!(ex.args[2], rec, result)
end

function get_vars!(ex::ExH{:tuple}, rec::Bool, result::Vector{Union{Symbol, Expr}})
    for arg in ex.args
        get_vars!(arg, rec, result)
    end
end

function get_vars!(ex::ExH{:block}, rec::Bool, result::Vector{Union{Symbol, Expr}})
    for subex in ex.args
        get_vars!(subex, rec, result)
    end
end

"""Get variables (`Symbol` or `Expr(:ref)`) involved in exprssion"""
function get_vars(ex::Expr; rec::Bool=false)
    result = Vector{Union{Symbol, Expr}}()
    get_vars!(ex, rec, result)
    return result
end

function get_vars(x::Symbol; rec::Bool=false)
    return Union{Symbol, Expr}[x]
end

function get_vars(x; rec::Bool=false)
    return Union{Symbol, Expr}[]
end

function get_var_names(ex; rec::Bool=false)
    vars = get_vars(ex; rec=rec)
    return [split_indexed(var)[1] for var in vars]
end

function get_indices(ex; rec::Bool=false)
    vars = get_vars(ex; rec=rec)
    return [split_indexed(var)[2] for var in vars]
end

"""
find_vars(ex; rec=true)

Same as `get_vars()`, but recursive by default
"""
find_vars(ex; rec::Bool=true) = get_vars(ex; rec=rec)
find_var_names(ex; rec::Bool=true) = get_var_names(ex; rec=rec)
find_indices(ex; rec::Bool=true) = get_indices(ex; rec=rec)


text:
# Definition of the composite type Alm

"An array of a_ℓm numbers."
mutable struct Alm{T <: Number}
    alm::Array{T}
    lmax::Int
    mmax::Int
    tval::Int

    Alm{T}(lmax, mmax) where {T <: Number} = new(zeros(T, numberOfAlms(lmax, mmax)),
                                                 lmax,
                                                 mmax,
                                                 2lmax + 1)

    function Alm{T}(lmax, mmax, arr::Array{T}) where {T <: Number}
        (numberOfAlms(lmax, mmax) == length(arr)) || throw(DomainError())

        new{T}(arr, lmax, mmax, 2lmax + 1)
    end
end

################################################################################

"""
    numberOfAlms(lmax::Integer, mmax::Integer) -> Integer
    numberOfAlms(lmax::Integer) -> Integer

Return the size of the array of complex numbers needed to store the
a_lm coefficients in the range of ℓ and m specified by `lmax` and
`mmax`. If `mmax` is not specified, it is assumed to be equal to
`lmax`. If `lmax` and `mmax` are inconsistent or negative, a
`DomainError` exception is thrown.
"""
function numberOfAlms(lmax, mmax)
    (lmax >= 0) || throw(DomainError(lmax, "`lmax` is not positive or zero"))
    (mmax >= 0) || throw(DomainError(mmax, "`mmax` is not positive or zero"))
    (0 ≤ mmax ≤ lmax) || throw(DomainError((lmax, mmax), "`lmax` and `mmax` are inconsistent"))

    div((mmax + 1) * (mmax + 2), 2) + (mmax + 1) * (lmax - mmax)
end

numberOfAlms(lmax) = numberOfAlms(lmax, lmax)

shr(x, y) = x >> y
shr(x::Array{T}, y) where {T} = [a >> y for a in x]

almIndexL0(alm::Alm{T}, m) where {T} = shr((m .* (alm.tval .- m)), 1) .+ 1
almIndex(alm::Alm{T}, l, m) where {T} = almIndexL0(alm, m) .+ l

################################################################################

"""
    readAlmFromFITS{T <: Complex}(f::FITSIO.FITSFile, t::Type{T}) -> Alm{T}
    readAlmFromFITS{T <: Complex}(fileName::String, t::Type{T}) -> Alm{T}

Read a set of a_ℓm coefficients from a FITS file. If the code fails,
FITSIO will raise an exception. (Refer to the FITSIO library for more
information.)
"""
function readAlmFromFITS(f::FITSIO.FITSFile,
                         t::Type{T}) where {T <: Complex}
    numOfRows = FITSIO.fits_get_num_rows(f)

    idx = Array{Int64}(undef, numOfRows)
    almReal = Array{Float64}(undef, numOfRows)
    almImag = Array{Float64}(undef, numOfRows)

    FITSIO.fits_read_col(f, 1, 1, 1, idx)
    FITSIO.fits_read_col(f, 2, 1, 1, almReal)
    FITSIO.fits_read_col(f, 3, 1, 1, almImag)

    l = floor.(Int64, sqrt.(idx .- 1))
    m = idx .- l.^2 .- l .- 1
    if count(x -> x < 0, m) > 0
        throw(DomainError())
    end

    result = Alm{T}(maximum(l), maximum(m))
    i = almIndex(result, l, m)
    result.alm = complex.(almReal[i], almImag[i])
end

function readAlmFromFITS(fileName,
                         t::Type{T}) where {T <: Complex}
    f = FITSIO.fits_open_table(fileName)
    try
        result = readAlmFromFITS(f, t)
        return result
    finally
        FITSIO.fits_close_file(f)
    end
end


text:
using Plots, Dates

export qwckplot

##  qwckplot function

"""
    qwckplot(fld::MeshArray)

Call qwckplot(fld::MeshArray) with date as title. Example:

```
!isdir("GRID_LLC90") ? error("missing files") : nothing
GridVariables=GridLoad(GridSpec("LLC90"))
qwckplot(GridVariables["Depth"])
```
"""
function qwckplot(fld::MeshArray)
    tmp1=Dates.now()
    tmp1=Dates.format(tmp1, "yyyy-mm-dd HH:MM:SS")
    qwckplot(fld,"Plotted at time "*tmp1)
end

"""
    qwckplot(fld::MeshArray,ttl::String)

Plot input using convert2array and heatmap + add title
"""
function qwckplot(fld::MeshArray,ttl::String)
    arr=convert2array(fld)
    arr=permutedims(arr,[2 1])
    #This uses Plots.jl:
    p=heatmap(arr,title=ttl)
end

"""
    qwckplot(fld::MeshArray,clims::NTuple{2, Number})

Plot input using convert2array and heatmap w. chosen clims
"""
function qwckplot(fld::MeshArray,clims::NTuple{2, Number})
    arr=convert2array(fld)
    arr=permutedims(arr,[2 1])
    #This uses Plots.jl:
    p=heatmap(arr,clims=clims)
end


text:
print("Enter ΔH in kJ/mol: ")
ΔH = tryparse(Float64, readline())
print("Enter ΔS in kJ/(mol*K): ")
ΔS = tryparse(Float64, readline())

ΔSPositive = ΔS > 0 # BOOLEAN VALUE (TRUE OR FALSE)
ΔHPositive = ΔH > 0

if ΔHPositive && !ΔSPositive
    println("Reaction is not spontaneous at any temperature")
elseif !ΔHPositive  && ΔSPositive
    println("Reaction is spontaneous at any temperature")
else
    temperatureDecription = ""
    if !ΔHPositive  && !ΔSPositive
        temperatureDecription = "low"
    else
        temperatureDecription = "high"
    end
    println("Spontaneous at $temperatureDecription Temperatures")
end

max_stars_repo_name: j-hayes/intro-to-computing-sci-eng

max_issues_repo_name: j-hayes/intro-to-computing-sci-eng

max_forks_repo_name: j-hayes/intro-to-computing-sci-eng


text:
JSONRPC.parse_params(::Type{Val{Symbol("\$/cancelRequest")}}, params) = CancelParams(params)
function process(r::JSONRPC.Request{Val{Symbol("\$/cancelRequest")},CancelParams}, server) end

JSONRPC.parse_params(::Type{Val{Symbol("\$/setTraceNotification")}}, params) = params
function process(r::JSONRPC.Request{Val{Symbol("\$/setTraceNotification")},Dict{String,Any}}, server) end

function JSONRPC.parse_params(::Type{Val{Symbol("julia/lint-package")}}, params) end
function process(r::JSONRPC.Request{Val{Symbol("julia/lint-package")},Nothing}, server) end

JSONRPC.parse_params(::Type{Val{Symbol("julia/toggle-lint")}}, params) = TextDocumentIdentifier(params["textDocument"])
function process(r::JSONRPC.Request{Val{Symbol("julia/toggle-lint")},TextDocumentIdentifier}, server)
    doc = server.documents[URI2(r.uri)]
    doc._runlinter = !doc._runlinter
end

function JSONRPC.parse_params(::Type{Val{Symbol("julia/reload-modules")}}, params) end
function process(r::JSONRPC.Request{Val{Symbol("julia/reload-modules")},Nothing}, server) end

JSONRPC.parse_params(::Type{Val{Symbol("julia/toggleFileLint")}}, params) = params
function process(r::JSONRPC.Request{Val{Symbol("julia/toggleFileLint")}}, server) end

function JSONRPC.parse_params(::Type{Val{Symbol("julia/toggle-log")}}, params) end
function process(r::JSONRPC.Request{Val{Symbol("julia/toggle-log")},Nothing}, server)
    server.debug_mode = !server.debug_mode
end

JSONRPC.parse_params(::Type{Val{Symbol("julia/getCurrentBlockOffsetRange")}}, params) = TextDocumentPositionParams(params)
function process(r::JSONRPC.Request{Val{Symbol("julia/getCurrentBlockOffsetRange")}}, server)
    if !haskey(server.documents, URI2(r.params.textDocument.uri))
        send(JSONRPC.Response(r.id, CancelParams(r.id)), server)
        return
    end
    tdpp = r.params
    doc = server.documents[URI2(tdpp.textDocument.uri)]
    offset = get_offset(doc, tdpp.position)
    x = getcst(doc)
    loc = 0
    p1, p2, p3 = 1, x.span, x.fullspan
    if typof(x) === CSTParser.FileH
        (offset > x.fullspan || x.args === nothing) && return 1, x.span, x.fullspan
        for a in x.args
            if loc <= offset < loc + a.fullspan
                if typof(a) === CSTParser.ModuleH
                    if loc + a.args[1].fullspan + a.args[2].fullspan < offset < loc + a.args[1].fullspan + a.args[2].fullspan + a.args[3].fullspan
                        loc0 = loc +  a.args[1].fullspan + a.args[2].fullspan
                        loc += a.args[1].fullspan + a.args[2].fullspan
                        for b in a.args[3].args
                            if loc <= offset < loc + b.fullspan
                                p1, p2, p3 = loc + 1, loc + b.span, loc + b.fullspan
                                break
                            end
                            loc += b.fullspan
                        end
                    else
                        p1, p2, p3 = loc + 1, loc + a.span, loc + a.fullspan
                    end
                elseif typof(a) === CSTParser.TopLevel
                    p1, p2, p3 = loc + 1, loc + a.span, loc + a.fullspan
                    for b in a.args
                        if loc <= offset < loc + b.fullspan
                            p1, p2, p3 = loc + 1, loc + b.span, loc + b.fullspan
                        end
                        loc += b.fullspan
                    end
                else
                    p1, p2, p3 = loc + 1, loc + a.span, loc + a.fullspan
                end
            end
            loc += a.fullspan
        end
    end

    response = JSONRPC.Response(r.id, (isempty(get_text(doc)) ? 0 : length(get_text(doc), 1, max(1, p1)), length(get_text(doc), 1, p2), length(get_text(doc), 1, p3)))

    send(response, server)
end

JSONRPC.parse_params(::Type{Val{Symbol("julia/activateenvironment")}}, params) = params
function process(r::JSONRPC.Request{Val{Symbol("julia/activateenvironment")}}, server)
    server.env_path = r.params
    server.symbol_server = StaticLint.SymbolServer.SymbolServerProcess(depot = server.depot_path, environment = server.env_path)
    @info "Restarted symbol server"
    StaticLint.SymbolServer.getstore(server.symbol_server)
    @info "StaticLint store set"
    kill(server.symbol_server)

    for (uri, doc) in server.documents
        parse_all(doc, server)
    end
    @info "Finished reparsing everything"
end

# SymbolServer:parallel branch ################################################

# function process(r::JSONRPC.Request{Val{Symbol("julia/activateenvironment")}}, server::LanguageServerInstance)
#     _set_worker_env(r.params, server)
# end

# function _set_worker_env(env_path, server::LanguageServerInstance)
#     wid = last(procs())
#     server.debug_mode && @info "Activating: ", server.env_path
#     @fetchfrom wid SymbolServer.Pkg.activate(env_path)

#     server.symbol_server.context = @fetchfrom wid SymbolServer.Pkg.Types.Context()
#     server.debug_mode && @info "New project file: ",  server.symbol_server.context.env.project_file

#     missing_pkgs = SymbolServer.disc_load_project(server.symbol_server)
#     if isempty(missing_pkgs)
#         server.ss_task = nothing
#     else
#         pkg_uuids = collect(keys(missing_pkgs))
#         server.debug_mode && @info "Missing or outdated package caches: ", collect(missing_pkgs)
#         server.ss_task = @spawnat wid SymbolServer.cache_packages_and_save(server.symbol_server.context, pkg_uuids)
#     end
# end
###############################################################################


text:
function MOI.get(regularizer::AbstractPenalizer, param::RawPenalizationParameter)
    name = Symbol(param.name)
    if !(name in fieldnames(typeof(penalizer)))
        error("Unrecognized parameter name: $(name) for penalizer $(typeof(penalizer)).")
    end
    return getfield(penalizer, name)
end

function MOI.set(penalizer::AbstractPenalizer, param::RawPenalizationParameter, value)
    name = Symbol(param.name)
    if !(name in fieldnames(typeof(penalizer)))
        error("Unrecognized parameter name: $(name) for penalizer $(typeof(penalizer)).")
    end
    setfield!(penalizer, name, value)
    return nothing
end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

__precompile__(true)
module Pkg

import Random
import REPL
using REPL.TerminalMenus

depots() = Base.DEPOT_PATH
logdir() = joinpath(depots()[1], "logs")
devdir() = get(ENV, "JULIA_PKG_DEVDIR", joinpath(homedir(), ".julia", "dev"))

# load snapshotted dependencies
include("../ext/TOML/src/TOML.jl")

include("GitTools.jl")
include("PlatformEngines.jl")
include("Types.jl")
include("Display.jl")
include("Pkg2/Pkg2.jl")
include("GraphType.jl")
include("Resolve.jl")
include("Operations.jl")
include("API.jl")
include("REPLMode.jl")

# Define new variables so tab comleting Pkg. works.
const add         = API.add
const rm          = API.rm
const up          = API.up
const test        = API.test
const gc          = API.gc
const init        = API.init
const build       = API.build
const installed   = API.installed
const pin         = API.pin
const free        = API.free
const checkout    = API.checkout
const develop     = API.develop
const generate    = API.generate
const instantiate = API.instantiate
const resolve     = API.resolve
const status      = Display.status
const update      = up
const activate    = API.activate

# legacy CI script support
import .API: clone, dir

import .REPLMode: @pkg_str
export @pkg_str

function __init__()
    if isdefined(Base, :active_repl)
        REPLMode.repl_init(Base.active_repl)
    else
        atreplinit() do repl
            if isinteractive() && repl isa REPL.LineEditREPL
                isdefined(repl, :interface) || (repl.interface = REPL.setup_interface(repl))
                REPLMode.repl_init(repl)
            end
        end
    end
end

using .Types
using UUIDs
import LibGit2
import Dates

METADATA_compatible_uuid(pkg::String) = Types.uuid5(Types.uuid_package, pkg)

include("precompile.jl")

end # module

max_stars_repo_name: helgee/julia

max_issues_repo_name: helgee/julia

max_forks_repo_name: helgee/julia


text:
function EpochDerivativesEquations(arg0::JString, arg1::NumericalPropagator)
    return EpochDerivativesEquations((JString, NumericalPropagator), arg0, arg1)
end

function compute_derivatives(obj::EpochDerivativesEquations, arg0::SpacecraftState, arg1::Vector{jdouble})
    return jcall(obj, "computeDerivatives", Vector{jdouble}, (SpacecraftState, Vector{jdouble}), arg0, arg1)
end

function get_mapper(obj::EpochDerivativesEquations)
    return jcall(obj, "getMapper", AbsoluteJacobiansMapper, ())
end

function get_name(obj::EpochDerivativesEquations)
    return jcall(obj, "getName", JString, ())
end

function get_selected_parameters(obj::EpochDerivativesEquations)
    return jcall(obj, "getSelectedParameters", ParameterDriversList, ())
end

function set_initial_jacobians(obj::EpochDerivativesEquations, arg0::SpacecraftState)
    return jcall(obj, "setInitialJacobians", SpacecraftState, (SpacecraftState,), arg0)
end

function set_initial_jacobians(obj::EpochDerivativesEquations, arg0::SpacecraftState, arg1::Vector{Vector{jdouble}}, arg2::Vector{Vector{jdouble}})
    return jcall(obj, "setInitialJacobians", SpacecraftState, (SpacecraftState, Vector{Vector{jdouble}}, Vector{Vector{jdouble}}), arg0, arg1, arg2)
end


text:
module PolynomialGTMTests

using PolynomialGTM, Test

@testset "Constructors" begin

    try
        GTM()
        @test true
    catch e
        @test throw(e)
    end

    try
        GTM(; stm=true)
        @test true
    catch e
        @test throw(e)
    end

end

@testset "Functions" begin

    try
        f = GTMFunction()
        f(randn(4), randn(2), NaN)
        @test true
    catch e
        @test throw(e)
    end

    try
        f = GTMFunction(; stm=true)
        f(randn(20), randn(2), NaN)
        @test true
    catch e
        @test throw(e)
    end

end

end

max_stars_repo_name: cadojo/PolynomialGTMModel

max_issues_repo_name: cadojo/PolynomialGTMModel

max_forks_repo_name: cadojo/PolynomialGTMModel


text:
@defcomp grosseconomy begin
	GA		= Variable(index=[time])	#Growth rate of productivity
	AL		= Variable(index=[time])	#Level of total factor productivity
    K       = Variable(index=[time])    #Capital stock (trillions 2010 US dollars)
    YGROSS  = Variable(index=[time])    #Gross world product GROSS of abatement and damages (trillions 2010 USD per year)

	a0		= Parameter()				#Initial level of total factor productivity
	ga0		= Parameter()				#Initial growth rate for TFP per 5 years
	dela	= Parameter()				#Decline rate of TFP per 5 years
    I       = Parameter(index=[time])   #Investment (trillions 2010 USD per year)
    l       = Parameter(index=[time])   #Level of population and labor
    dk      = Parameter()               #Depreciation rate on capital (per year)
    gama    = Parameter()               #Capital elasticity in production function
    k0      = Parameter()               #Initial capital value (trill 2010 USD)

    function run_timestep(p, v, d, t)
		#Define function for GA
		for tm = 1:100
			v.GA[tm] = p.ga0 * exp(-p.dela * 5 * (tm-1))
		end

		#Define function for AL
		if is_first(t)
			v.AL[t] = p.a0
		else
			v.AL[t] = v.AL[t-1]/(1 - v.GA[t-1])
		end

        #Define function for K
        if is_first(t)
            v.K[t] = p.k0
        else
            v.K[t] = (1 - p.dk)^5 * v.K[t-1] + 5 * p.I[t-1]
        end

        #Define function for YGROSS
        v.YGROSS[t] = (v.AL[t] * (p.l[t]/1000)^(1-p.gama)) * (v.K[t]^p.gama)
    end
end


text:
@adjoint StepRange{Int64, Int64}(a::Int64, b::Int64) = a:b, Δ -> (nothing, nothing)
@adjoint StepRange{Int64, Int64}(a::Int64, b::Int64, c::Int64) = a:b:c, Δ -> (nothing, nothing, nothing)

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
using OrdinaryDiffEq, DiffEqProblemLibrary, Test

using DiffEqProblemLibrary.ODEProblemLibrary: importodeproblems; importodeproblems()
import DiffEqProblemLibrary.ODEProblemLibrary: prob_ode_2Dlinear,
                                               prob_ode_large2Dlinear,
                                               prob_ode_linear,
                                               prob_ode_2Dlinear_notinplace

u0=rand(300,20).*ones(300,20)/2
prob = prob_ode_2Dlinear_notinplace
prob2 = prob_ode_2Dlinear

sol =solve(prob,Euler(),dt=1//2^(4),save_everystep=false)
sol =solve(prob2,Euler(),dt=1//2^(4),save_everystep=false)

alloc1 = @allocated sol =solve(prob,Euler(),dt=1//2^(6),save_everystep=false)
alloc2 = @allocated sol2 =solve(prob2,Euler(),dt=1//2^(6),save_everystep=false)

alloc1 = @allocated sol =solve(prob,Euler(),dt=1//2^(6),save_everystep=false)
alloc2 = @allocated sol2 =solve(prob2,Euler(),dt=1//2^(6),save_everystep=false)

@test alloc2 <= alloc1

sol = solve(prob_ode_linear,Euler(),dt=1//2^(6),save_everystep=true)
sol2 = solve(prob_ode_linear,Euler(),sol[:],sol.t,sol.k;dt=1//2^(8),save_everystep=true)

sol = solve(prob_ode_large2Dlinear,Euler(),dt=1//2^(6),save_everystep=true)
sol2 = solve(prob_ode_large2Dlinear,Euler(),sol[:],sol.t,sol.k;dt=1//2^(8),save_everystep=true)

sol = solve(prob_ode_large2Dlinear,Euler(),dt=1//2^(6),save_everystep=true)
alloc1 = @allocated sol = solve(prob_ode_large2Dlinear,Euler(),dt=1//2^(8),save_everystep=true)
alloc2 = @allocated sol2 = solve(prob_ode_large2Dlinear,Euler(),sol[:],sol.t,sol.k;dt=1//2^(8),save_everystep=true)

@test alloc2 <= alloc1

sol =solve(prob,Tsit5())


text:
module Leapfrog

import Base.length
export init_cost, init, init_all, solve

include("Types.jl")
include("Utilities.jl")
include("Last.jl")
include("Upper.jl")

function solve(cost::CostParameters, mp::RLSModel)
    ess = ESS(cost.nC) # Start with feasible ESS; we're solving for all equilibria
    number_of_equilibria=0
    τ = Array(Stage, cost.nC)
    for i in reverse(1:cost.nC)
        τ[i] = Stage(Array(EQ, i, i, 5))
    end

    solve_last_corner!(cost, mp, τ)
    solve_last_edge!(cost, mp, τ)
    solve_last_interior!(cost, mp, τ, ess)
    for ic = reverse(1:cost.nC-1)
        solve_corner!(cost, mp, ic, τ, ess)
        if ic != cost.nC
            solve_edges!(cost, mp, ic, τ, ess)
            solve_interior!(cost, mp, ic, τ, ess)
        end
    end

    number_of_equilibria += 1
    jump!(ess)
    if ess[1,1,1] != 0
        for void=1:1000000000
            if ess.ess_to_ic[ess.active][3] == ess.nC
                solve_last_interior!(cost, mp, ess.ess_to_ic[ess.active][1], ess.ess_to_ic[ess.active][2], τ, ess)
            elseif ess.ess_to_ic[ess.active][3] < ess.nC
                solve_interior!(cost, mp, ess.ess_to_ic[ess.active][3], ess.ess_to_ic[ess.active][1], ess.ess_to_ic[ess.active][2], τ, ess)
            end

            for ic = reverse(1:ess.ess_to_ic[ess.active][3]-1) # already solved stage!
                solve_corner!(cost, mp, ic, τ, ess)
                if ic != cost.nC
                    solve_edges!(cost, mp, ic, τ, ess)
                    solve_interior!(cost, mp, ic, τ, ess)
                end
            end
            number_of_equilibria+=1
            jump!(ess)
            mod(void, 100000) == 0 ? println(void) :
            if ess[1,1,1] == 0
                break
            end
        end
    end
    return τ, ess, number_of_equilibria
end

function init_cost(switch::Int64)
    if switch == 1
        k1 = 8.3
        k2 = 1.
        C  = linspace(5, 0, 4)
        Cmin = minimum(C)
        nC = length(C)
        p = [1.; 1.; 1.; 0]
        π = 1-p
        return CostParameters(k1, k2, C, Cmin, nC, p, π)
    elseif switch == 2
        k1 = 8.3
        k2 = 1.
        C  = linspace(5, 0, 5)
        Cmin = minimum(C)
        nC = length(C)
        p = [1.; 1.; 1.; 1.; 0.]
        π = 1-p
        return CostParameters(k1, k2, C, Cmin, nC, p, π)
    elseif switch == 3
        k1 = 8.3
        k2 = 1.
        C  = linspace(5, 0, 3)
        Cmin = minimum(C)
        nC = length(C)
        p = [1.; 1.; 0.]
        π = 1-p
        return CostParameters(k1, k2, C, Cmin, nC, p, π)
    end
end

function init_mp(cost)
    R = 0.05
    Δt = 1
    β = exp(- R*Δt) # R is continuous time iterest rate
    T = (cost.nC-1)*3 + 1
    return RLSModel(R, β, Δt, T, number_of_strings(cost.nC))
end

function init_all(switch::Int64)
    cost = init_cost(switch)
    return init_mp(cost), cost
end
end


text:
################################################################################
#
#	CONCRETE TYPE
#
#   ReciprocalPath{P} <: AbstractReciprocalPath{P}
#   --> P is the reciprocal point type (<: AbstractReciprocalPoint{D})
#       --> D is the dimension of embedding space
#
#   FILE CONTAINS
#       - concrete struct definition
#       - interface implementation
#
################################################################################

################################################################################
#
#   STRUCT DEFINITION
#
################################################################################
mutable struct ReciprocalPath{P} <: AbstractReciprocalPath{P}

    # the only field is the array of points
    points :: Vector{P}

end

# export the type
export ReciprocalPath

################################################################################
#
#	INTERFACING / ACCESSING RECIPROCAL PATHS
#
################################################################################

# default constructor interface
# used for creation of new reciprocal paths
function newReciprocalPath(
            ::Type{ReciprocalPath{RPO}},
            points :: Vector{<:RPO}
        ) :: ReciprocalPath{RPO} where {D,RPO<:AbstractReciprocalPoint{D}}

    # return the newly created object
    return ReciprocalPath{RPO}(points)
end

# export the newUnitcell interface
export newReciprocalPath

# access the list of points
function points(
            p :: ReciprocalPath{RPO}
        ) :: Vector{RPO} where {D,RPO<:AbstractReciprocalPoint{D}}

    # return the respective field
    return p.points
end
# set the list of points
function points!(
            p :: ReciprocalPath{RPO},
            points :: Vector{<:RPO}
        ) where {D,RPO<:AbstractReciprocalPoint{D}}

    # set the respective field
    p.points = points
end

# export access to points
export points, points!


text:
# Implementation of the SKAN neuron model from:
#   "Racing to Learn: Statistical Inference and Learning in a Single Spiking Neuron with Adaptive Kernels" Afshar et al 2014
#   "Turn Down That Noise: Synaptic Encoding of Afferent SNR in a Single Spiking Neuron" Afshar et al 2014

# "A component of a SKAN architecture"
# abstract AbstractSkan

export
    SkanParams,
    InputNeuron,
    SkanNeuron,
    SkanGroup,
    SkanSynapse,
    membrane_potential,
    step!

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

using Parameters

@with_kw immutable SkanParams{T<:Real}
    # Δt::Int = 1
    # ramp_step_rise::T =
    ramp_step_adjustment::T = 1  # ddr
    ramp_step_min::T = 0         #
    ramp_step_max::T = 400
    ramp_step_initfunc::Function = () -> T(rand(100:200))
    # threshold_rise_initfunc::Function = (n) -> 40n
    # threshold_fall_initfunc::Function = (n) -> 100n
    # inhibitor_decay::T
    weight_init::T = 10000
    # weight_initfunc::Function = () -> 10000
    # weight_rise::T
    # weight_fall::T
end

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

"Represents an input spike train"
type InputNeuron <: SpikingNeuron
    spike::Bool
end
InputNeuron() = InputNeuron(false)
step!(neuron::InputNeuron) = nothing

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

"""
A SkanNeuron is modeled on Afshar et al's Synapto-Dendritic Kernel Adaptation (SKAN) framework,
which is a straightforward and efficient approximation of delayed spiking characteristics, and a linearly
updating alternative to the Spike Response Model (SRM), Hodgen-Huxley, or other similarly complex neuronal
models.  The model avoids multiplication and division, and has a smart normalization step using only left/right
shifts, thus compute power can be used for additional model expressivity instead of biological correctness.
"""
type SkanNeuron{T<:Real, S} <: SpikingNeuron
    incoming_synapses::S
    # outgoing_synapses::S
    pulse::Bool         # s(t): the soma output (pulse)... binary: on or off
    spike::Bool         # u(t): the initial spike. true when: !s(t-1) && s(t)
    potential::T        # ∑rᵢ(t): membrane potential
    threshold::T        # θ(t):   threshold voltage
    threshold_rise::T   # threshold change when increasing
    threshold_fall::T   # threshold change when decreasing
end

function membrane_potential{T<:Real}(neuron::SkanNeuron{T})
    sumrt = zero(T)
    for synapse in neuron.incoming_synapses
        sumrt += synapse.ramp
    end
    sumrt
end

# update the neurons' membrane potentials and update pulse/spike/threshold/potential
#   pulse(t) = sum(rᵢ(t)) > θ(t-1)  (for i ∈ incoming_synapses)
#   spike(t) = pulse(t) && !pulse(t-1)
function step!(neuron::SkanNeuron)
    potential = membrane_potential(neuron)

    # pulse while potential is above threshold
    # spike only in the first timestep of a pulse
    pulse = potential >= neuron.threshold
    neuron.spike = pulse && !neuron.pulse
    neuron.pulse = pulse

    # during a pulse, increase the threshold
    if pulse
        neuron.threshold += neuron.threshold_rise

    # when potential returns to rest, decrease the threshold
    elseif potential <= 0 && neuron.potential > 0
        neuron.threshold -= neuron.threshold_fall
    end

    # save the potential
    neuron.potential = potential
end

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

# group of SKAN neurons which have a shared inhibitory signal
type SkanGroup{T<:Real} <: AbstractVector{SkanNeuron{T}}
    neurons::Vector{SkanNeuron{T}}
    inhibitor::T
end
Base.size(group::SkanGroup) = size(group.neurons)
Base.getindex(group::SkanGroup, i::Integer) = group.neurons[i]

# connects a presynaptic SkanNeuron to a postsynaptic SkanNeuron
type SkanSynapse{T<:Real, S1<:SpikingNeuron, S2<:SpikingNeuron}
    presynaptic::S1
    postsynaptic::S2
    weight::T           # wᵢ(t): synaptic weight
    weight_adj_flag::T  # dᵢ(t): 1 on uᵢ(t)
    ramp::T             # rᵢ(t): ramp height
    ramp_step::T        # Δrᵢ(t): ramp height change per time period Δt
    ramp_flag::Int8     # pᵢ(t): state of the ramp-up/ramp-down flag {-1, 0, 1}
end

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

"One time step in a simulation"
function step!{N<:SpikingNeuron, S<:SkanSynapse}(params::SkanParams,
                                              neurons::AbstractArray{N},
                                              synapses::AbstractArray{S})

    # update the ramp vars for each synapse
    #   r(t) += p(t-1) * Δr(t-1)
    #   Δr(t) += ddr * s(t-1)
    for synapse in synapses
        synapse.ramp += synapse.ramp_flag * synapse.ramp_step
        synapse.ramp = max(0, synapse.ramp)
        if synapse.postsynaptic.pulse
            synapse.ramp_step += synapse.ramp_flag * params.ramp_step_adjustment
        end
    end

    # update the neurons
    for neuron in neurons
        step!(neuron)
    end

    # update the synaptic ramp flags
    for synapse in synapses
        if synapse.ramp_flag > 0 && synapse.ramp >= synapse.weight
            synapse.ramp_flag = -1
        elseif synapse.ramp_flag < 0 && synapse.ramp <= 0
            synapse.ramp_flag = 0
        elseif synapse.ramp_flag == 0 && synapse.presynaptic.spike
            synapse.ramp_flag = 1
        end
    end

end

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


text:
juser@juliabox:~$ /opt/julia-0.5/bin/julia eight_queen_puzzle.jl
n = 1
elapsed:  0.000001 seconds
solutions = 1

n = 2
elapsed:  0.000001 seconds
solutions = 0

n = 3
elapsed:  0.000001 seconds
solutions = 0

n = 4
elapsed:  0.000001 seconds
solutions = 2

n = 5
elapsed:  0.000003 seconds
solutions = 10

n = 6
elapsed:  0.000008 seconds
solutions = 4

n = 7
elapsed:  0.000028 seconds
solutions = 40

n = 8
elapsed:  0.000108 seconds
solutions = 92

n = 9
elapsed:  0.000463 seconds
solutions = 352

n = 10
elapsed:  0.002146 seconds
solutions = 724

n = 11
elapsed:  0.010646 seconds
solutions = 2680

n = 12
elapsed:  0.057603 seconds
solutions = 14200

n = 13
elapsed:  0.334600 seconds
solutions = 73712

n = 14
elapsed:  2.055078 seconds
solutions = 365596

n = 15
elapsed: 13.480449 seconds
solutions = 2279184

n = 16
elapsed: 97.192552 seconds
solutions = 14772512

n = 17
elapsed:720.314676 seconds
solutions = 95815104


text:
#julia data structure
#---------array----------#
#julia index of array is start with 1 ,not 0
d = ["fsds","fdsfs",1]
print(d)
print(d[1])
print(d[2])
print(d[3])

a = Array{String,1}

#------------priority queue-------------#
#pq = Collections.priorityQueue()
#pq["a"] = 15
#pq["b"] = 332
#pq["a"] = 0x32

#print(pq["a"])

#pd

max_stars_repo_name: zqsheng/snippet

max_issues_repo_name: zqsheng/snippet

max_forks_repo_name: zqsheng/snippet


text:
@doc raw"""
    device_range(canonical_model::CanonicalModel,
                        range_data::Vector{NamedMinMax},
                        cons_name::Symbol,
                        var_name::Symbol)

Constructs min/max range constraint from device variable.

# Constraints
If min and max within an epsilon width:

``` variable[r[1], t] == r[2].max ```

Otherwise:

``` r[2].min <= variable[r[1], t] <= r[2].max ```

where r in range_data.

# LaTeX

`` x = r^{max}, \text{ for } |r^{max} - r^{min}| < \varepsilon ``

`` r^{min} \leq x \leq r^{max}, \text{ otherwise } ``

# Arguments
* canonical_model::CanonicalModel : the canonical model built in PowerSimulations
* range_data::Vector{NamedMinMax} : contains name of device (1) and its min/max (2)
* cons_name::Symbol : name of the constraint
* var_name::Symbol : the name of the continuous variable
"""
function device_range(canonical_model::CanonicalModel,
                    range_data::Vector{NamedMinMax},
                    cons_name::Symbol,
                    var_name::Symbol)

    time_steps = model_time_steps(canonical_model)
    variable = var(canonical_model, var_name)
    set_name = (r[1] for r in range_data)
    _add_cons_container!(canonical_model, cons_name, set_name, time_steps)
    constraint = con(canonical_model, cons_name)

    for r in range_data
          if abs(r[2].min - r[2].max) <= eps()
            @warn("The min - max values in range constraint with eps() distance to each other. Range Constraint will be modified for Equality Constraint")
                for t in time_steps
                    constraint[r[1], t] = JuMP.@constraint(canonical_model.JuMPmodel, variable[r[1], t] == r[2].max)
                end
          else
                for t in time_steps
                    constraint[r[1], t] = JuMP.@constraint(canonical_model.JuMPmodel, r[2].min <= variable[r[1], t] <= r[2].max)
                end
            end
    end

    return

end

@doc raw"""
    device_semicontinuousrange(canonical_model::CanonicalModel,
                                    scrange_data::Vector{NamedMinMax},
                                    cons_name::Symbol,
                                    var_name::Symbol,
                                    binvar_name::Symbol)

Constructs min/max range constraint from device variable and on/off decision variable.

# Constraints
If device min = 0:

``` varcts[r[1], t] <= r[2].max*varbin[r[1], t]) ```

``` varcts[r[1], t] >= 0.0 ```

Otherwise:

``` varcts[r[1], t] <= r[2].max*varbin[r[1], t] ```

``` varcts[r[1], t] >= r[2].min*varbin[r[1], t] ```

where r in range_data.

# LaTeX

`` 0 \leq x^{cts} \leq r^{max} x^{bin}, \text{ for } r^{min} = 0 ``

`` r^{min} x^{bin} \leq x^{cts} \leq r^{max} x^{bin}, \text{ otherwise } ``

# Arguments
* canonical_model::CanonicalModel : the canonical model built in PowerSimulations
* scrange_data::Vector{NamedMinMax} : contains name of device (1) and its min/max (2)
* cons_name::Symbol : name of the constraint
* var_name::Symbol : the name of the continuous variable
* binvar_name::Symbol : the name of the binary variable
"""
function device_semicontinuousrange(canonical_model::CanonicalModel,
                                    scrange_data::Vector{NamedMinMax},
                                    cons_name::Symbol,
                                    var_name::Symbol,
                                    binvar_name::Symbol)

    time_steps = model_time_steps(canonical_model)
    ub_name = _middle_rename(cons_name, "_", "ub")
    lb_name = _middle_rename(cons_name, "_", "lb")

    varcts = var(canonical_model, var_name)
    varbin = var(canonical_model, binvar_name)

    #MOI has a semicontinous set, but after some tests is not clear most MILP solvers support it.
    #In the future this can be updated

    set_name = (r[1] for r in scrange_data)
    _add_cons_container!(canonical_model, ub_name, set_name, time_steps)
    _add_cons_container!(canonical_model, lb_name, set_name, time_steps)
    con_ub = con(canonical_model, ub_name)
    con_lb = con(canonical_model, lb_name)

    for t in time_steps, r in scrange_data

        # If the variable was a lower bound != 0, not removing the LB can cause infeasibilities
        if JuMP.has_lower_bound(varcts[r[1], t])
            JuMP.set_lower_bound(varcts[r[1], t], 0.0)
        end

        if r[2].min == 0.0

            con_ub[r[1], t] = JuMP.@constraint(canonical_model.JuMPmodel, varcts[r[1], t] <= r[2].max*varbin[r[1], t])
            con_lb[r[1], t] = JuMP.@constraint(canonical_model.JuMPmodel, varcts[r[1], t] >= 0.0)

        else

            con_ub[r[1], t] = JuMP.@constraint(canonical_model.JuMPmodel, varcts[r[1], t] <= r[2].max*varbin[r[1], t])
            con_lb[r[1], t] = JuMP.@constraint(canonical_model.JuMPmodel, varcts[r[1], t] >= r[2].min*varbin[r[1], t])

        end

    end

    return

end

@doc raw"""
    reserve_device_semicontinuousrange(canonical_model::CanonicalModel,
                                    scrange_data::Vector{NamedMinMax},
                                    cons_name::Symbol,
                                    var_name::Symbol,
                                    binvar_name::Symbol)

Constructs min/max range constraint from device variable and on/off decision variable.

# Constraints
If device min = 0:

``` varcts[r[1], t] <= r[2].max*(1-varbin[r[1], t]) ```

``` varcts[r[1], t] >= 0.0 ```

Otherwise:

``` varcts[r[1], t] <= r[2].max*(1-varbin[r[1], t]) ```

``` varcts[r[1], t] >= r[2].min*(1-varbin[r[1], t]) ```

where r in range_data.

# LaTeX

`` 0 \leq x^{cts} \leq r^{max} (1 - x^{bin} ), \text{ for } r^{min} = 0 ``

`` r^{min} (1 - x^{bin} ) \leq x^{cts} \leq r^{max} (1 - x^{bin} ), \text{ otherwise } ``

# Arguments
* canonical_model::CanonicalModel : the canonical model built in PowerSimulations
* scrange_data::Vector{NamedMinMax} : contains name of device (1) and its min/max (2)
* cons_name::Symbol : name of the constraint
* var_name::Symbol : the name of the continuous variable
* binvar_name::Symbol : the name of the binary variable
"""
function reserve_device_semicontinuousrange(canonical_model::CanonicalModel,
                                            scrange_data::Vector{NamedMinMax},
                                            cons_name::Symbol,
                                            var_name::Symbol,
                                            binvar_name::Symbol)

    time_steps = model_time_steps(canonical_model)
    ub_name = _middle_rename(cons_name, "_", "ub")
    lb_name = _middle_rename(cons_name, "_", "lb")

    varcts = var(canonical_model, var_name)
    varbin = var(canonical_model, binvar_name)

    # MOI has a semicontinous set, but after some tests is not clear most MILP solvers support it.
    # In the future this can be updated

    set_name = (r[1] for r in scrange_data)
    _add_cons_container!(canonical_model, ub_name, set_name, time_steps)
    _add_cons_container!(canonical_model, lb_name, set_name, time_steps)
    con_ub = con(canonical_model, ub_name)
    con_lb = con(canonical_model, lb_name)

    for t in time_steps, r in scrange_data

            if r[2].min == 0.0

                con_ub[r[1], t] = JuMP.@constraint(canonical_model.JuMPmodel, varcts[r[1], t] <= r[2].max*(1-varbin[r[1], t]))
                con_lb[r[1], t] = JuMP.@constraint(canonical_model.JuMPmodel, varcts[r[1], t] >= 0.0)

            else

                con_ub[r[1], t] = JuMP.@constraint(canonical_model.JuMPmodel, varcts[r[1], t] <= r[2].max*(1-varbin[r[1], t]))
                con_lb[r[1], t] = JuMP.@constraint(canonical_model.JuMPmodel, varcts[r[1], t] >= r[2].min*(1-varbin[r[1], t]))

            end

    end

    return

 end


text:
#=
  What is the color of the taxi?

  From https://www.bayesia.com/2018-03-02-probabilistic-reasoning-under-uncertainty

  Originally from Kahnemann, Slovic, Tversky "Judgement under uncertainty"

  There has been an accicent involving a taxi.
  There are two kinds of taxis:
    - yellow taxi: 85%
    - white taxi: 15%

  A witness say: It was a white taxi involved in the accident.

  Researcher:
    - 80% of all witness statements are true
    - 20% of all witness statements are false.

  What is the probability that it was a white taxi involved in
  the accident?

  Answer: The probability that it was a white taxi involved in the accident
          is about 41-42%. And it's the same as in the talk.

  This model:
  Distributions of variable involved
  yellow     =>    5897  (0.589700)
  white      =>    4103  (0.410300)

  Distributions of variable witness_white (num:0)
  1.00000 =>    9999  (0.999900)
  0.00000 =>       1  (0.000100)

  Distributions of variable witness_yellow (num:0)
  1.00000 =>    5614  (0.561400)
  0.00000 =>    4386  (0.438600)

  Cf ~/blog/color_of_the_taxi.blog
     ~/psi/color_of_the_taxi.psi
     ~/webppl/color_of_the_taxi.wppl

=#

using Turing, StatsPlots, DataFrames
include("jl_utils.jl")

@model function color_of_the_taxi()
    white = 1
    yellow = 2
    #  Prior distributions of the different taxis.
    involved ~ Categorical([0.15,0.85]) # [white,yellow]

    #  Witness says color but is is only x percent reliable.
    #  Witness experts states that a witness can only be 80% reliable
    #  (given the same circumstances as the accidents).
    witness = tzeros(2)
    for c in 1:2
        if c == involved
            witness[c] ~ flip(0.80)
        else
            witness[c] ~ flip(0.20)
        end
    end

    true ~ Dirac(witness[white] == true)
    # true ~ Dirac(witness[yellow] == true)

    witness_white ~ Dirac(witness[white])
    witness_yellow ~ Dirac(witness[yellow])

end

model = color_of_the_taxi()
num_chains = 4

# chns = sample(model, Prior(), 10_000)
# chns = sample(model, MH(), 10_000)
chns = sample(model, PG(5), 10_000)
# chns = sample(model, SMC(), 10_000)
# chns = sample(model, IS(), 10_000)

display(chns)
# display(plot(chns))

println("white: 1 yellow: 2")
show_var_dist_pct(chns,:involved,["white","yellow"])
show_var_dist_pct(chns,:witness_white)
show_var_dist_pct(chns,:witness_yellow)

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank