Julia^2: Julia sets in Julia, using Interact. I'd like it faster.

JuliaSquared.jl

using Plots, Interact, Base.Threads, BenchmarkTools

# only called by resize.
function batchify(maxIndex, numBatches)
    if mod(maxIndex, numBatches) == 0
        n = numBatches
    else
        n = numBatches - 1
    end
    batchSize = div(maxIndex, n)
    lastBatch = batchSize*n+1:maxIndex
    batches = [(i-1)*batchSize+1:i*batchSize for i in 1:n]
    if length(lastBatch) > 0
        push!(batches, lastBatch)
    end
    batches
end

numThreads = Threads.nthreads() - 1

# scale should change with picture size
scale = 350.0 # the greater the more details
shift = 0.0+0.0im

function resize(w, h)
    global width  = w
    global height = h
    # Complex number grid, to compute Julia(., c) for.
    # axis have odd numbers by construction. good?
#    global xAxis = [1.0y/scale for y in -div(width,  2):1:div(width,  2)]
#    global yAxis = [1.0x/scale for x in -div(height, 2):1:div(height, 2)]
#    global valueGrid = [x+y*im for y in yAxis, x in xAxis]
#    # Here we will put our function values.
#    global dataGrid = [UInt8(0) for x = 1:height, y = 1:width]
#    global batches = batchify(height, numThreads)
    :ok
end

# start values for picture size
resize(800, 700)

const MAXITER = UInt32(255)

function julia(z0, c, maxiter) :: UInt32
    z = z0
    for i = UInt32(1):maxiter
        z = z*z + c
        if abs2(z) >= 4.0
            return i
        end
    end
    return UInt32(0)
end

# I'm not happy with the parallelization yet
function juliaSet(batches, valueGrid, dataGrid, c)
    @threads for batch in batches
        for x in batch
            for y = 1:width
                @inbounds dataGrid[x,y] = julia(valueGrid[x,y], c, MAXITER)
            end
        end
    end
end

# juliaP and juliaSetP by Simon Danisch. May contain errors added by zenon ...
function juliaP(z0, c, maxiter)
    z = z0
    for i in 1:maxiter
        abs2(z) > 4f0 && return (i - 1) % UInt32
        z = z * z + c
    end
    return 0 % UInt32 # % is used to convert without overflow check
end

function juliaSetP(valueGrid, dataGrid, c)
    Threads.@threads for i in eachindex(dataGrid)
        @inbounds dataGrid[i] = juliaP(valueGrid[i], c, MAXITER)
    end
    return dataGrid
end

function j(c)
    xAxis = [1.0y/scale for y in -div(width,  2):1:div(width,  2)][1:end-1]
    yAxis = [1.0x/scale for x in -div(height, 2):1:div(height, 2)][1:end-1]
    valueGrid = [x+y*im for y in yAxis, x in xAxis]
    # Here we will put our function values.
    dataGrid = [UInt8(0) for x = 1:height, y = 1:width]
    batches = batchify(height, numThreads)
    println("value : $(size(valueGrid))")
    println("data : $(size(dataGrid))")
    println("simon:")
    @btime juliaSetP($valueGrid, $dataGrid, $c)
    println("zenon:")
    @btime juliaSet($batches, $valueGrid, $dataGrid, $c)
    # color = :ice ?
    #heatmap(xAxis, yAxis, dataGrid, size=(width,height), nbins=256, aspect_ratio=1, title="c = $c")
    :ok
end

j() = j(-0.62+0.42im)
# -0.62+0.42im
# 0.39+0.6im
# 0.0+0.64im
# 0.0+0.75im
# -0.512511498387847167+ 0.521295573094847167im https://en.wikipedia.org/wiki/File:Julia_set,_plotted_with_Matplotlib.svg

function p()
    xAxis = [1.0y/scale for y in -div(width,  2):1:div(width,  2)][1:end-1]
    yAxis = [1.0x/scale for x in -div(height, 2):1:div(height, 2)][1:end-1]
    valueGrid = [x+y*im for y in yAxis, x in xAxis]
    # Here we will put our function values.
    dataGrid = [UInt8(0) for x = 1:height, y = 1:width]
    batches = batchify(height, numThreads)
    @manipulate throttle = 0.1 for x in -1.0:0.01:1.0, y in -1.0:0.01:1.0
        c = x+y*im
        # @time juliaSet(batches, valueGrid, dataGrid, c)
        @time juliaSetP(valueGrid, dataGrid, c)
        heatmap(xAxis, yAxis, dataGrid, size=(width,height), nbins=256, aspect_ratio=1, title="c = $c")
    end
end

#display(p())
#display(j())