text:
###
# Main
###

using MAT
using Datetime
using Distributions
import NumericExtensions.sumabs
import NumericExtensions
const ne = NumericExtensions

include("config.jl")
include("drawcircle.jl")
include("gaussian_blur.jl")
include("utils.jl")

###
# set up initial configuration
###
function main(config=Dict(), disturbances=Dict(), cluster=false; enableVis :: Bool = false, fileName = "", loadTime = nothing, allow32Bit = false, forceJuliaImpl = false, resultFolder="results")
    useVis = enableVis && ((length(procs()) == 1) || (!(myid() in workers()))) && !cluster
    ###
    # Prepare GPU
    ###
    const P64BIT, USECL, ctx, queue = determineCapabilities(cluster, allow32Bit, forceJuliaImpl)

    ###
    # Prepare GUI
    ###

    guiproc = if useVis && (length(procs()) <= 1)
        first(addprocs(1))
    elseif useVis
        last(procs())
    else
        -1
    end

    rref = if useVis
       @at_proc guiproc using PyPlot
    end

    value = simulation(cluster, useVis, fileName, loadTime, USECL, P64BIT ? Float64 : Float32, config, disturbances, guiproc, ctx, queue, rref, resultFolder)

    ctx = nothing
    queue = nothing
    gc()
    return value
end

function simulation{T <: FloatingPoint}(cluster, enableVis, fileName, loadTime, USECL, :: Type{T}, config :: Dict, disturbances :: Dict, guiproc :: Int, ctx, queue, gui_rref, resultFolder)
worker = (length(procs()) > 1) && (myid() in workers())

# initialize membrane fields
Afield = nothing
Mfield = nothing
Ffield = nothing
Wfield = nothing
directionfield = nothing

# Prepare the configuration values and load the default settings
loadConfig(baseConfig)

# If we have a non empty fileName and we don't run on a cluster
# Load configuration values from data directory.
fileConfig = if fileName != "" && !cluster
    matread("data/$(fileName).mat")
else
    Dict()
end

# Load file config
loadConfig(fileConfig, dataVars)

# Overwrite configuration values from a
loadConfig(config, dataVars)
updateDependentValues() # Update the values that depend on the configuration

# Load initial state.
if loadTime != nothing
    history_W, history_A, history_M, history_F, history_dir =
        if cluster && checkVars(historyVars, collect(keys(config)))
            (config["history_W"], config["history_A"], config["history_M"], config["history_F"], config["history_dir"])
        elseif checkVars(historyVars, collect(keys(fileConfig)))
            (fileConfig["history_W"], fileConfig["history_A"], fileConfig["history_M"], fileConfig["history_F"], fileConfig["history_dir"])
        else
            error("Could not find historyVars despite being given a loadTime of $loadTime")
        end

    if loadTime == -1
        Wfield = history_W[:,:,end]
        Afield = history_A[:,:,end]
        Mfield = history_M[:,:,end]
        Ffield = history_F[:,:,end]
        directionfield = history_dir[:,:,end]
    else
        Wfield = history_W[:,:,loadTime]
        Afield = history_A[:,:,loadTime]
        Mfield = history_M[:,:,loadTime]
        Ffield = history_F[:,:,loadTime]
        directionfield = history_dir[:,:,loadTime]
    end
else
    Afield = zeros(T, fieldResY, fieldResX)
    Mfield = zeros(T, fieldResY, fieldResX)
    Ffield = avgF * ones(T, fieldResY, fieldResX)
    Wfield = ones(T, fieldResY, fieldResX)
    directionfield = pi/4 * ones(T, fieldResY, fieldResX)

    # draw membrane circle
    M_circ = zeros(T, fieldResY, fieldResX)

    for r in mR:0.01:(mR+mT-1)
        drawcircle!(M_circ, xc, yc, r)
    end
    Mfield += avgM * M_circ + 0.00001* rand(fieldResY, fieldResX)

    # fill in with autocatalyst
    A_circ = zeros(T, fieldResY, fieldResX)

    for r in 0:0.01:mR-1
        drawcircle!(A_circ, xc, yc, r)
    end

    Afield += avgA * A_circ

    Wfield -= (Mfield + Afield)

    ###
    # directionality initialization
    ###

    # make directionality in a ring
    for q1 in 1:fieldResX
        for q2 in 1:fieldResY
            directionfield[q2,q1] = atan((q2-yc)/(q1-xc))-pi/2

            # correct for 180 degrees
            while directionfield[q2, q1] > pi
                 directionfield[q2, q1] -= pi
            end
            while directionfield[q2, q1] <= 0 # ??? < 0
                 directionfield[q2, q1] += pi
            end
        end
    end

    directionfield[isnan(directionfield)] = 0
end

###
# Data storing
###

# create 3d matrices to store field activities
history_A = Array{T, 2}[]
history_F = Array{T, 2}[]
history_T = Array{T, 2}[]
history_M = Array{T, 2}[]
history_W = Array{T, 2}[]
history_dir = Array{T, 2}[]

#vectors to save global concentrations across time

Avec = T[]
Fvec = T[]
Mvec = T[]
Wvec = T[]
DAvec = T[]
DMvec = T[]

###
# Decide which compute context should be used
# Either OpenCL or Julia
###

if USECL
    include("cl/functions.jl")
else
    include("jl/functions.jl")
end

createComputeContext(T, ctx, queue)

###
# Prepare simulation
###

#define scaled diffusion
diffA = diffusionA*fieldRes/fieldSize
diffF = diffusionF*fieldRes/fieldSize
diffM = diffusionM*fieldRes/fieldSize
diffW = diffusionW*fieldRes/fieldSize

dF = zeros(T, fieldResY, fieldResX)

# create temp arrays
mpot = create(T)
mpot1 = create(T)
mpot2 = create(T)
wpot = create(T)
apot = create(T)
apot1 = create(T)

buff_Mfield = create(T)
buff_Wfield = create(T)
buff_Afield = create(T)
buff_Dfield = create(T)
buff_NewDfield = create(T)
buff_Ffield = create(T)

wlap = create(T)
alap = create(T)
mlap = create(T)
flap = create(T)

row1 = create(T)
row2 = create(T)
row3 = create(T)
row4 = create(T)
row5 = create(T)
row6 = create(T)

dW = create(T)
dA = create(T)
dM = create(T)
buff_dF = create(T)

temp = create(T)
area = create_n4(T)
const_area = create_const_n4(convert(T, 1.0/8.0))

# Make sure that the fields are in the correct DataFormat
Mfield = convert(Array{T}, Mfield)
Afield = convert(Array{T}, Afield)
Wfield = convert(Array{T}, Wfield)
Ffield = convert(Array{T}, Ffield)
directionfield = convert(Array{T}, directionfield)

# Stability and structure
stable = false
timeToStable = -1

old_Mfield = copy(Mfield)
old_Ffield = copy(Ffield)
old_Wfield = copy(Wfield)
old_Afield = copy(Afield)
old_directionfield = copy(directionfield)

###
# Simulation
###

t = 1

meanMField = mean(Mfield)
meanAField = mean(Afield)

###
# Disturbance steps
###
tT = timeTotal
outFileName = ""

###
# Set up dynamic precision
###

precision = 1.0
skipFlag = false
precisionCounter = 0
bufferStep = 2 # Minimum 1.0 increase to allow for higher precision over a second timeStep

previous_Mfield = copy(Mfield)
previous_Ffield = copy(Ffield)
previous_Wfield = copy(Wfield)
previous_Afield = copy(Afield)
previous_directionfield = copy(directionfield)

tx = Float64[]

while (t <= tT) && (meanAField < 1.0) !isnan(meanMField) && !isnan(meanAField)
    previous_Mfield = copy(Mfield)
    previous_Ffield = copy(Ffield)
    previous_Wfield = copy(Wfield)
    previous_Afield = copy(Afield)
    previous_directionfield = copy(directionfield)

	if t in keys(disturbances)
		val = disturbances[t]
		method = first(val)
		args = val[2:end]
		if method == :global && length(args) == 2
			(mu, sig) = args
            mask = !FrefillBinMask
			Mfield += mask .* rand(Normal(mu, sig), (fieldResY, fieldResX))
			Afield += mask .* rand(Normal(mu, sig), (fieldResY, fieldResX))
			Ffield += mask .* rand(Normal(mu, sig), (fieldResY, fieldResX))
			Wfield += mask .* rand(Normal(mu, sig), (fieldResY, fieldResX))
			directionfield += mask .* rand(Normal(mu, sig), (fieldResY, fieldResX))

			# Normalize
			Mfield[Mfield .< 0] = 0
			Afield[Afield .< 0] = 0
			Ffield[Ffield .< 0] = 0
			Wfield[Wfield .< 0] = 0

			directionfield = map(x -> mod(x, pi), directionfield)
		elseif method == :punch_local
			(x, y, alpha, beta) = args
			apply_punch_down!(Mfield, x, y, alpha, beta)
			apply_punch_down!(Afield, x, y, alpha, beta)
			apply_punch_down!(Ffield, x, y, alpha, beta)
			apply_punch_down!(Wfield, x, y, alpha, beta)
            # apply_punch_down!(directionfield, x, y, alpha, beta)
        elseif method == :punch_random
            (alpha, beta) = args
            lin_idx = sample(find((s) -> s > 0.5, Mfield))
            x, y = ind2sub(size(Mfield), lin_idx)
            apply_punch_down!(Mfield, x, y, alpha, beta)
            apply_punch_down!(Afield, x, y, alpha, beta)
            apply_punch_down!(Ffield, x, y, alpha, beta)
            apply_punch_down!(Wfield, x, y, alpha, beta)
		elseif method == :gaussian_blur && length(args) == 1
			sigma = args[1]
			imfilter_gaussian_no_nans!(Mfield, [sigma,sigma])
			imfilter_gaussian_no_nans!(Afield, [sigma,sigma])
			imfilter_gaussian_no_nans!(Ffield, [sigma,sigma])
			imfilter_gaussian_no_nans!(Wfield, [sigma,sigma])
			# directionfield = imfilter_gaussian(directionfield, [sigma,sigma])
		elseif method == :tear_membrane && length(args) == 1
			tearsize = args[1]

		    tS = iround(fieldRes/2+tearSize*fieldRes/fieldSize)
		    t0 = iround(fieldRes/2)
		    Mfield[fieldRes/2:fieldRes,t0:tS] = 0
		elseif method == :tear_dfield && length(args) == 1
			tearsize = args[1]

		    tS = iround(fieldRes/2+tearSize*fieldRes/fieldSize)
		    t0 = iround(fieldRes/2)
		    directionfield[fieldRes/2:fieldRes,t0:tS] = pi .* rand(iround(fieldRes-fieldRes/2+1),tS-t0+1)
		elseif method == :tear && length(args) == 1
			tearsize = args[1]

		    tS = iround(fieldRes/2+tearSize*fieldRes/fieldSize)
		    t0 = iround(fieldRes/2)

		    Mfield[fieldRes/2:fieldRes,t0:tS] = 0
		    directionfield[fieldRes/2:fieldRes,t0:tS] = pi .* rand(iround(fieldRes-fieldRes/2+1),tS-t0+1)
		else
			warn("Don't know how to handle $method with arguments $args, carry on.")
		end
	end

    ###
    # Pushing the fields on the GPU
    ###

    c_copy!(buff_Mfield, Mfield)
    c_copy!(buff_Wfield, Wfield)
    c_copy!(buff_Afield, Afield)
    c_copy!(buff_Ffield, Ffield)
    c_copy!(buff_Dfield, directionfield)

    ###
    # calculate potential based on repulsion
    ###
    area!(buff_Dfield, area)

    potential!(buff_Mfield, buff_Wfield, area, mpot, wpot, MW_repulsion)
    potential!(buff_Mfield, buff_Afield, area, mpot1, apot, MA_repulsion)
    potential!(buff_Mfield, buff_Mfield, area, mpot2, temp, MM_repulsion)

    potential!(buff_Afield, buff_Afield, const_area, apot1, temp, AA_repulsion)

    c_add!(apot, apot1)
    c_add!(mpot, mpot1)
    c_add!(mpot, mpot2)

    ###
    # move molecules and update directionality
    ###

    diffusion!(buff_Mfield, mpot, mlap)
    smul!(mlap, diffM)

    diffusion!(buff_Wfield, wpot, wlap)
    smul!(wlap, diffW)

    diffusion!(buff_Afield, apot, alap)
    smul!(alap, diffA)

    # Laplacian for diffusion
    # Todo calculate in on GPU to be coherent and to minimize data transfers.
    laplacian!(buff_Ffield, flap)
    smul!(flap, diffF)

    # update direction field based on alignment
    align!(buff_Mfield, buff_Dfield, buff_NewDfield)
    c_copy!(directionfield, buff_NewDfield)

    ###
    # reactions
    ##

    calcRow!(buff_Mfield, buff_Afield, buff_Ffield, buff_Wfield, row1, m11, a11, f11, w11)
    calcRow!(buff_Mfield, buff_Afield, buff_Ffield, buff_Wfield, row2, m12, a12, f12, w12)
    calcRow!(buff_Mfield, buff_Afield, buff_Ffield, buff_Wfield, row3, m21, a21, f21, w21)
    calcRow!(buff_Mfield, buff_Afield, buff_Ffield, buff_Wfield, row4, m22, a22, f22, w22)
    calcRow!(buff_Mfield, buff_Afield, buff_Ffield, buff_Wfield, row5, m31, a31, f31, w31)
    calcRow!(buff_Mfield, buff_Afield, buff_Ffield, buff_Wfield, row6, m32, a32, f32, w32)

    ##
    # Calculate dA
    ##

    delta!(row1, row2, row3, row4, row5, row6, dA,
            fsm(a12, a11, kf1), fsm(a11, a12, kb1), fsm(a22, a21, kf2), fsm(a21, a22, kb2), fsm(a32, a31, kf3), fsm(a31, a32, kb3))
    delta2!(dA, buff_Mfield, buff_Afield, alap, dA, decayA)

    ##
    # Calculate dM
    ##

    delta!(row1, row2, row3, row4, row5, row6, dM,
            fsm(m12, m11, kf1), fsm(m11, m12, kb1), fsm(m22, m21, kf2), fsm(m21, m22, kb2), fsm(m32, m31, kf3), fsm(m31, m32, kb3))
    delta2!(dM, buff_Afield, buff_Mfield, mlap, dM, decayM)

    ##
    # Calculate dW
    ##

    delta!(row1, row2, row3, row4, row5, row6, dW,
            fsm(w12, w11, kf1), fsm(w11, w12, kb1), fsm(w22, w21, kf2), fsm(w21, w22, kb2), fsm(w32, w31, kf3), fsm(w31, w32, kb3))

    c_add!(dW, wlap)

    ##
    # Calculate dF
    ##

    delta!(row1, row2, row3, row4, row5, row6, buff_dF,
            fsm(f12, f11, kf1), fsm(f11, f12, kb1), fsm(f22, f21, kf2), fsm(f21, f22, kb2), fsm(f32, f31, kf3), fsm(f31, f32, kb3))

    c_add!(buff_dF, flap)

    c_copy!(dF, buff_dF)

    dF[FrefillBinMask] += flowRateF * (saturationF .- Ffield[FrefillBinMask])

    dT = 0

    # update values
    smul!(dA, (stepIntegration / precision))
    c_add!(buff_Afield, dA)
    c_copy!(Afield, buff_Afield)

    Ffield += dF * (stepIntegration / precision)

    smul!(dM, (stepIntegration / precision))
    c_add!(buff_Mfield, dM)
    c_copy!(Mfield, buff_Mfield)

    smul!(dW, (stepIntegration / precision))
    c_add!(buff_Wfield, dW)
    c_copy!(Wfield, buff_Wfield)

    # calulate stability criteria

    sumabs_dA = sumabs(c_read(dA))
    sumabs_dM = sumabs(c_read(dM))
    sumabs_dF = sumabs(dF)
    sumabs_dW = sumabs(c_read(dW))

    negativeValues = any(Mfield .< 0) || any(Afield .< 0) || any(Ffield .< 0) || any(Wfield .< 0)

    if negativeValues
        warn("Integration error at t=$t increase precision from $precision to $(2 * precision)")
        precision = 2 * precision

        if precision > 1024
            error("Precision just exceeded 1024 clearly something is wrong here.")
        end

        Mfield = copy(previous_Mfield)
        Ffield = copy(previous_Ffield)
        Wfield = copy(previous_Wfield)
        Afield = copy(previous_Afield)
        directionfield = copy(previous_directionfield)
        skipFlag = true

        # If we increase it for the first time
        if precisionCounter == 0
            precisionCounter = bufferStep*precision
        else
            precisionCounter *= 2
        end
    end

    #save values for visualization
    meanAField = mean(Afield)
    meanMField = mean(Mfield)

    if !skipFlag
        stableCondition = (sumabs_dA < epsilon / precision) && (sumabs_dM < epsilon / precision) #&& (sumabs_dF < epsilon) && (sumabs_dW < epsilon)

        if !stable && stableCondition
            println("reached a stable config after $t")
            stable = true
            tT = t + stableTime
            timeToStable = t
        elseif stable && !stableCondition
            stable = false
            tT = t + timeTotal
            println("Lost stability")
        end

        push!(Avec, meanAField)
        push!(Mvec, meanMField)
        push!(Fvec, mean(Ffield))
        push!(Wvec, mean(Wfield))
        push!(DAvec, sumabs_dA)
        push!(DMvec, sumabs_dM)

        if (t % storeStep == 0)
          push!(history_A, Afield)
          push!(history_F, Ffield)
          push!(history_M, Mfield)
          push!(history_W, Wfield)
          push!(history_dir, directionfield)
        end

        push!(tx, t)
    end

    if enableVis && (t % visInterval == 0)
        if !isready(gui_rref)
            println()
            println("Waiting for GUI to be initialized.")
            @time wait(gui_rref)
        end

        @spawnat guiproc begin
          hold(false)
          # Timeseries plot

          subplot(241)
          plot(tx, Avec, "-", linewidth=2)
          title("Avec")

          subplot(242)
          plot(tx, Mvec, "-", linewidth=2)
          title("Mvec")

          subplot(243)
          plot(tx, DAvec, "-", linewidth=2)
          title("DAvec")
          axis([0, t, 0, 0.4])

          subplot(244)
          plot(tx, DMvec, "-", linewidth=2)
          title("DMvec")
          axis([0, t, 0, 0.4])

          subplot(245)
          pcolormesh(Mfield, vmin=0, vmax=0.6)
          title("Mfield")

          subplot(246)
          pcolormesh(Afield, vmin=0, vmax=0.6)
          title("Afield")

          subplot(247)
          pcolormesh(Ffield, vmin=0, vmax=1)
          title("Ffield")

            subplot(248)
            title("DirectionField")

            U = cos(directionfield)
            V = sin(directionfield)
            d1, d2 = size(directionfield)
            plt.quiver([1:d1], [1:d2], U, V, linewidth=1.5, headwidth = 0.5)
        end
    end

    yield() # to be able to answer question about the status
    if !skipFlag
        if precisionCounter > 0
            precisionCounter -= 1
            t += (stepIntegration / precision)
        else
            t += stepIntegration
        end
    else
        skipFlag = false
    end

    if precisionCounter == 0 && precision > 1.0
        if t % stepIntegration != 0
            warn("t=$t and integration=$stepIntegration t % integration should be zero but it is. $(t%stepIntegration)")
        end
        warn("Resetting precision to 1.0 at t=$t")
        precision = 1.0
    end
end # While

###
# Find the reason we terminated
###

print("Finished because after $t: ")
if t > tT
  println("Time done")
elseif isnan(meanMField) || isnan(meanAField)
   println("Found NaNs")
elseif meanAField >= 1.0
    println("Mean of Afield is $meanAField >= 1.0")
end

###
# Export
###
try
    dt=now()
    stack(X) = reshape(hcat(X...), fieldResY, fieldResX, length(X))
    result = {
        "history_W" => stack(history_W),
        "history_A" => stack(history_A),
        "history_M" => stack(history_M),
        "history_F" => stack(history_F),
        "history_dir" => stack(history_dir),
        "Avec" => Avec,
        "Fvec" => Fvec,
        "Mvec" => Mvec,
        "Wvec" => Wvec,
        "DAvec" => DAvec,
        "DMvec" => DMvec,
        "tx" => tx
    }
    merge!(result, saveConfig(collect(keys(baseConfig))))
    prefix = if cluster
        "$resultFolder/$(myid())-"
    else
        "$resultFolder/"
    end
    outFileName = "$prefix$(year(dt))-$(month(dt))-$(day(dt))_$(hour(dt)):$(minute(dt)):$(second(dt)).mat"
    matwrite(outFileName, result)
catch e
    warn("Write failed because of $e")
end

if !(worker || cluster)
println("Press any key to exit program.")
readline(STDIN)
end

structM = sumabs(old_Mfield .- Mfield)
structA = sumabs(old_Afield .- Afield)
structF = sumabs(old_Ffield .- Ffield)
structW = sumabs(old_Wfield .- Wfield)
structd = sumabs(old_directionfield .- directionfield)

valueOfInterest =   if stable
                        if (meanMField >= deathEpsilon) || (meanAField >= deathEpsilon)
                            timeToStable # stable and living
                        else
                            -timeToStable # death through diffusion and decay
                        end
                    elseif (meanAField >= 1.0) || isnan(meanMField) || isnan(meanAField) # explosion
                        - t
                    end

destroyComputeContext()

return (t, timeToStable, valueOfInterest, stable, meanMField, meanAField, structM, structA, structF, structW, structd, outFileName)
end #Function

function fsm(x, y, k)
  if (k == zero(k)) || (x == y)
    return zero(k)
  else
    return ((x-y) * k )
  end
end


text:
# Eckstein, Bertsekas, "On the Douglas-Rachford Splitting Method and the
# Proximal Point Algorithm for Maximal Monotone Operators",
# Mathematical Programming, vol. 55, no. 1, pp. 293-318 (1989).

using Base.Iterators
using ProximalAlgorithms.IterationTools
using ProximalOperators: Zero
using LinearAlgebra
using Printf

struct DRS_iterable{R <: Real, C <: Union{R, Complex{R}}, T <: AbstractArray{C}, Tf, Tg}
    f::Tf
    g::Tg
    x::T
    gamma::R
end

Base.IteratorSize(::Type{<:DRS_iterable}) = Base.IsInfinite()

mutable struct DRS_state{T}
    x::T
    y::T
    r::T
    z::T
    res::T
end

DRS_state(iter::DRS_iterable) = DRS_state(copy(iter.x), zero(iter.x), zero(iter.x), zero(iter.x), zero(iter.x))

function Base.iterate(iter::DRS_iterable, state::DRS_state=DRS_state(iter))
    prox!(state.y, iter.f, state.x, iter.gamma)
    state.r .= 2 .*state.y .- state.x
    prox!(state.z, iter.g, state.r, iter.gamma)
    state.res .= state.y .- state.z
    state.x .-= state.res
    return state, state
end

# Solver

struct DouglasRachford{R}
    gamma::R
    maxit::Int
    tol::R
    verbose::Bool
    freq::Int

    function DouglasRachford{R}(; gamma::R, maxit::Int=1000, tol::R=R(1e-8),
        verbose::Bool=false, freq::Int=100
    ) where R
        @assert gamma > 0
        @assert maxit > 0
        @assert tol > 0
        @assert freq > 0
        new(gamma, maxit, tol, verbose, freq)
    end
end

function (solver::DouglasRachford{R})(
    x0::AbstractArray{C}; f=Zero(), g=Zero()
) where {R, C <: Union{R, Complex{R}}}

    stop(state::DRS_state) = norm(state.res, Inf) / solver.gamma <= solver.tol
    disp((it, state)) = @printf("%5d | %.3e\n", it, norm(state.res, Inf) / solver.gamma)

    iter = DRS_iterable(f, g, x0, solver.gamma)
    iter = take(halt(iter, stop), solver.maxit)
    iter = enumerate(iter)
    if solver.verbose iter = tee(sample(iter, solver.freq), disp) end

    num_iters, state_final = loop(iter)

    return state_final.y, state_final.z, num_iters

end

# Outer constructors

"""
    DouglasRachford([gamma, maxit, tol, verbose, freq])

Instantiate the Douglas-Rachford splitting algorithm (see [1]) for solving
convex optimization problems of the form

    minimize f(x) + g(x),

If `solver = DouglasRachford(args...)`, then the above problem is solved with

    solver(x0, [f, g])

Optional keyword arguments:

* `gamma::Real` (default: `1.0`), stepsize parameter.
* `maxit::Integer` (default: `1000`), maximum number of iterations to perform.
* `tol::Real` (default: `1e-8`), absolute tolerance on the fixed-point residual.
* `verbose::Bool` (default: `true`), whether or not to print information during the iterations.
* `freq::Integer` (default: `100`), frequency of verbosity.

References:

[1] Eckstein, Bertsekas, "On the Douglas-Rachford Splitting Method and the
Proximal Point Algorithm for Maximal Monotone Operators",
Mathematical Programming, vol. 55, no. 1, pp. 293-318 (1989).
"""
DouglasRachford(::Type{R}; kwargs...) where R = DouglasRachford{R}(; kwargs...)
DouglasRachford(; kwargs...) = DouglasRachford(Float64; kwargs...)

ext: jl


text:
function check_num_nodes(g::GNNGraph, x::AbstractArray)
    @assert g.num_nodes == size(x, ndims(x))
end
function check_num_edges(g::GNNGraph, e::AbstractArray)
    @assert g.num_edges == size(e, ndims(e))
end

sort_edge_index(eindex::Tuple) = sort_edge_index(eindex...)

function sort_edge_index(u, v)
    uv = collect(zip(u, v))
    p = sortperm(uv) # isless lexicographically defined for tuples
    return u[p], v[p]
end

cat_features(x1::Nothing, x2::Nothing) = nothing
cat_features(x1::AbstractArray, x2::AbstractArray) = cat(x1, x2, dims=ndims(x1))
cat_features(x1::Union{Number, AbstractVector}, x2::Union{Number, AbstractVector}) =
    cat(x1, x2, dims=1)

function cat_features(x1::NamedTuple, x2::NamedTuple)
    sort(collect(keys(x1))) == sort(collect(keys(x2))) || @error "cannot concatenate feature data with different keys"

    NamedTuple(k => cat_features(getfield(x1,k), getfield(x2,k)) for k in keys(x1))
end

# Turns generic type into named tuple
normalize_graphdata(data::Nothing; kws...) = NamedTuple()

normalize_graphdata(data; default_name::Symbol, kws...) =
normalize_graphdata(NamedTuple{(default_name,)}((data,)); default_name, kws...)

function normalize_graphdata(data::NamedTuple; default_name, n, duplicate_if_needed=false)
    # This had to workaround two Zygote bugs with NamedTuples
    # https://github.com/FluxML/Zygote.jl/issues/1071
    # https://github.com/FluxML/Zygote.jl/issues/1072

    if n == 1
        # If last array dimension is not 1, add a new dimension.
        # This is mostly usefule to reshape globale feature vectors
        # of size D to Dx1 matrices.
        function unsqz(v)
            if v isa AbstractArray && size(v)[end] != 1
                v = reshape(v, size(v)..., 1)
            end
            v
        end

        data = NamedTuple{keys(data)}(unsqz.(values(data)))
    end

    sz = map(x -> x isa AbstractArray ? size(x)[end] : 0, data)
    if duplicate_if_needed
        # Used to copy edge features on reverse edges
        @assert all(s -> s == 0 ||  s == n || s == n÷2, sz)

        function duplicate(v)
            if v isa AbstractArray && size(v)[end] == n÷2
                v = cat(v, v, dims=ndims(v))
            end
            v
        end
        data = NamedTuple{keys(data)}(duplicate.(values(data)))
    else
        @assert all(s -> s == 0 ||  s == n, sz)
    end
    return data
end

ones_like(x::AbstractArray, T=eltype(x), sz=size(x)) = fill!(similar(x, T, sz), 1)
ones_like(x::SparseMatrixCSC, T=eltype(x), sz=size(x)) = ones(T, sz)
ones_like(x::CUMAT_T, T=eltype(x), sz=size(x)) = CUDA.ones(T, sz)

numnonzeros(a::AbstractSparseMatrix) = nnz(a)
numnonzeros(a::AbstractMatrix) = count(!=(0), a)

# each edge is represented by a number in
# 1:N^2
function edge_encoding(s, t, n; directed=true)
    if directed
        # directed edges and self-loops allowed
        idx = (s .- 1) .* n .+ t
        maxid = n^2
    else
        # Undirected edges and self-loops allowed
        maxid = n * (n + 1) ÷ 2

        mask = s .> t
        snew = copy(s)
        tnew = copy(t)
        snew[mask] .= t[mask]
        tnew[mask] .= s[mask]
        s, t = snew, tnew

        # idx = ∑_{i',i'<i} ∑_{j',j'>=i'}^n 1 + ∑_{j',i<=j'<=j} 1
        #     = ∑_{i',i'<i} ∑_{j',j'>=i'}^n 1 + (j - i + 1)
        #     = ∑_{i',i'<i} (n - i' + 1) + (j - i + 1)
        #     = (i - 1)*(2*(n+1)-i)÷2 + (j - i + 1)
        idx = @. (s-1)*(2*(n+1)-s)÷2 + (t-s+1)
    end
    return idx, maxid
end

# each edge is represented by a number in
# 1:N^2
function edge_decoding(idx, n; directed=true)
    if directed
        # g = remove_self_loops(g)
        s =  (idx .- 1) .÷ n .+ 1
        t =  (idx .- 1) .% n .+ 1
    else
        # We replace j=n in
        # idx = (i - 1)*(2*(n+1)-i)÷2 + (j - i + 1)
        # and obtain
        # idx = (i - 1)*(2*(n+1)-i)÷2 + (n - i + 1)

        # OR We replace j=i  and obtain??
        # idx = (i - 1)*(2*(n+1)-i)÷2 + 1

        # inverting we have
        s = @. ceil(Int, -sqrt((n + 1/2)^2 - 2*idx) + n + 1/2)
        t = @. idx - (s-1)*(2*(n+1)-s)÷2 - 1 + s
        # t =  (idx .- 1) .% n .+ 1
    end
    return s, t
end

@non_differentiable edge_encoding(x...)
@non_differentiable edge_decoding(x...)

####################################
# FROM MLBASE.jl
# https://github.com/JuliaML/MLBase.jl/pull/1/files
# remove when package is registered
##############################################

numobs(A::AbstractArray{<:Any, N}) where {N} = size(A, N)

# 0-dim arrays
numobs(A::AbstractArray{<:Any, 0}) = 1

function getobs(A::AbstractArray{<:Any, N}, idx) where N
    I = ntuple(_ -> :, N-1)
    return A[I..., idx]
end

getobs(A::AbstractArray{<:Any, 0}, idx) = A[idx]

function getobs!(buffer::AbstractArray, A::AbstractArray{<:Any, N}, idx) where N
    I = ntuple(_ -> :, N-1)
    buffer .= A[I..., idx]
    return buffer
end

# --------------------------------------------------------------------
# Tuples and NamedTuples

_check_numobs_error() =
    throw(DimensionMismatch("All data containers must have the same number of observations."))

function _check_numobs(tup::Union{Tuple, NamedTuple})
    length(tup) == 0 && return
    n1 = numobs(tup[1])
    for i=2:length(tup)
        numobs(tup[i]) != n1 && _check_numobs_error()
    end
end

function numobs(tup::Union{Tuple, NamedTuple})::Int
    _check_numobs(tup)
    return length(tup) == 0 ? 0 : numobs(tup[1])
end

function getobs(tup::Union{Tuple, NamedTuple}, indices)
    _check_numobs(tup)
    return map(x -> getobs(x, indices), tup)
end

function getobs!(buffers::Union{Tuple, NamedTuple},
                  tup::Union{Tuple, NamedTuple},
                  indices)
    _check_numobs(tup)

    return map(buffers, tup) do buffer, x
                getobs!(buffer, x, indices)
            end
end
#######################################################


text:
using Test
using FastFilter
using Random

utypes = [UInt8, UInt16, UInt32]
error_rates = [0.005, 1.5e-5, 2.33e-10]
rng = MersenneTwister(42)

for i in 1:length(utypes)
    @testset "binaryfusefilter.jl " begin
        n = Int(1e5)
        items = rand(rng, UInt64, n)
        filter = BinaryFuseFilter{utypes[i]}(items)

        @testset "true positives should always be inside" begin
            for i in items
                @test i in filter
            end
        end

        not_in_items = rand(rng, UInt64, Int(1e6))
        not_in_items = setdiff(not_in_items, items)

        @testset "false positives should be below threshold" begin
            fp = 0
            for i in not_in_items
                fp += i in filter
            end
            @test fp/length(not_in_items) < error_rates[i]
        end

        filter2 = BinaryFuseFilter{utypes[i]}(items; seed = UInt64(42))
        @testset "seed should influence construction" begin
            for i in items
                @test i in filter2
            end

            fp = 0
            for i in not_in_items
                fp += i in filter
            end
            @test fp/length(not_in_items) < error_rates[i]
            println("False positive rate for " * string(utypes[i]) *
                " filter: " * string(fp/length(not_in_items)))

            # and finally combined filters should have lower fp
            fp2 = 0
            for i in not_in_items
                fp2 += (i in filter) & (i in filter2)
            end

            @test fp2 <= fp
        end
    end
end


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule p3a_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("p3a")
JLLWrappers.@generate_main_file("p3a", UUID("31b8d6fc-6c78-56bd-b44f-ae5989ed963d"))
end  # module p3a_jll


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct LinuxConfiguration <: SwaggerModel
    disablePasswordAuthentication::Any # spec type: Union{ Nothing, Bool } # spec name: disablePasswordAuthentication
    ssh::Any # spec type: Union{ Nothing, SshConfiguration } # spec name: ssh
    provisionVMAgent::Any # spec type: Union{ Nothing, Bool } # spec name: provisionVMAgent

    function LinuxConfiguration(;disablePasswordAuthentication=nothing, ssh=nothing, provisionVMAgent=nothing)
        o = new()
        validate_property(LinuxConfiguration, Symbol("disablePasswordAuthentication"), disablePasswordAuthentication)
        setfield!(o, Symbol("disablePasswordAuthentication"), disablePasswordAuthentication)
        validate_property(LinuxConfiguration, Symbol("ssh"), ssh)
        setfield!(o, Symbol("ssh"), ssh)
        validate_property(LinuxConfiguration, Symbol("provisionVMAgent"), provisionVMAgent)
        setfield!(o, Symbol("provisionVMAgent"), provisionVMAgent)
        o
    end
end # type LinuxConfiguration

const _property_map_LinuxConfiguration = Dict{Symbol,Symbol}(Symbol("disablePasswordAuthentication")=>Symbol("disablePasswordAuthentication"), Symbol("ssh")=>Symbol("ssh"), Symbol("provisionVMAgent")=>Symbol("provisionVMAgent"))
const _property_types_LinuxConfiguration = Dict{Symbol,String}(Symbol("disablePasswordAuthentication")=>"Bool", Symbol("ssh")=>"SshConfiguration", Symbol("provisionVMAgent")=>"Bool")
Base.propertynames(::Type{ LinuxConfiguration }) = collect(keys(_property_map_LinuxConfiguration))
Swagger.property_type(::Type{ LinuxConfiguration }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_LinuxConfiguration[name]))}
Swagger.field_name(::Type{ LinuxConfiguration }, property_name::Symbol) =  _property_map_LinuxConfiguration[property_name]

function check_required(o::LinuxConfiguration)
    true
end

function validate_property(::Type{ LinuxConfiguration }, name::Symbol, val)
end


text:
type Treeplex{IntType<:Integer}
    dimension::IntType
    numinfosets::IntType
    simplexbegin::Array{IntType}
        simplexend::Array{IntType}
        simplexparentsequence::Array{IntType}
end

function prox(alpha::Real, g::Array{Real}, beta::Real, y::Array{Real})

end


text:
using Dash, DashHtmlComponents, DashCoreComponents
using PlotlyJS
using JSON
using Dates

using CSV
using DataFrames

df = CSV.read("/home/prabu/Research/mintsData/rawJld/001e06305a6c_df.csv", DataFrame)

function plotPM()

    OPCN2_pm1 = scatter(x = df.dateTime, y = df.OPCN2_pm1,
                        mode = "lines+markers", name = "OPCN2 PM1",
                        marker_color = :red, marker_size = 4)

    OPCN2_pm2_5 = scatter(x = df.dateTime, y = df.OPCN2_pm2_5,
                        mode = "lines+markers", name = "OPCN2 PM2.5",
                        marker_color = :blue, marker_size = 4)

    OPCN2_pm10 = scatter(x = df.dateTime, y = df.OPCN2_pm10,
                        mode = "lines+markers", name = "OPCN2 PM10",
                        marker_color = :black, marker_size = 4)

    layout = Layout(title = "PM1", xaxis_title = "DateTime", yaxis_title = "PM1",
                  plot_bgcolor = :transparent, paper_bgcolor = :transparent)

  plot([OPCN2_pm1, OPCN2_pm2_5, OPCN2_pm10], layout)
end

function plotCO2()

    SCD30_c02 = scatter(x = df.dateTime, y = df.SCD30_c02,
                        mode = "lines+markers", name = "SCD30",
                        marker_color = :red, marker_size = 4)

    layout = Layout(title = "CO2", xaxis_title = "DateTime", yaxis_title = "CO2",
                  plot_bgcolor = :transparent, paper_bgcolor = :transparent)

  plot([SCD30_c02], layout)
end

function plotTemp()

    AS7262_temp = scatter(x = df.dateTime, y = df.AS7262_temperature,
                  mode = "lines+markers",name = "AS7262",
                  marker_color = :red,marker_size = 4)

    BME280_temp = scatter(x = df.dateTime, y = df.BME280_temperature,
                  mode = "lines+markers", name = "BME280",
                  marker_color = :blue, marker_size = 4)

    SCD30_temp = scatter(x = df.dateTime, y = df.SCD30_temperature,
                 mode = "lines+markers", name = "SCD30",
                 marker_color = :black, marker_size = 4)

    layout = Layout(title = "Temperature", xaxis_title = "DateTime", yaxis_title = "Temperature",
                    plot_bgcolor = :transparent, paper_bgcolor = :transparent)

    plot([AS7262_temp, BME280_temp, SCD30_temp], layout)
end

function plotPressure()

    BME280_pressure = scatter(x = df.dateTime, y = df.BME280_pressure,
                        mode = "lines+markers", name = "BME280",
                        marker_color = :red, marker_size = 4)

    layout = Layout(title = "Pressure", xaxis_title = "DateTime", yaxis_title = "Pressure",
                  plot_bgcolor = :transparent, paper_bgcolor = :transparent)

  plot([BME280_pressure], layout)
end

function plotHumidity()

    BME280_humidity = scatter(x = df.dateTime, y = df.BME280_humidity,
                        mode = "lines+markers", name = "BME280",
                        marker_color = :red, marker_size = 4)

    SCD30_humidity = scatter(x = df.dateTime, y = df.SCD30_humidity,
                        mode = "lines+markers", name = "SCD30",
                        marker_color = :blue, marker_size = 4)

    layout = Layout(title = "Humidity", xaxis_title = "DateTime", yaxis_title = "Humidity",
                  plot_bgcolor = :transparent, paper_bgcolor = :transparent)

  plot([BME280_humidity, SCD30_humidity], layout)
end

app = dash()

app.layout = html_div() do
    html_h1(
        "MINTS Dashboards: UTD Nodes", style=Dict("color"=>:black, "textAlign"=>"center")),

    html_div(
        children = [
            dcc_graph(id="pm1-graph", figure = plotPM()),
            dcc_graph(id="CO2-graph", figure = plotCO2()),
            dcc_graph(id="temp-graph", figure = plotTemp()),
            dcc_graph(id="pressure-graph", figure = plotPressure()),
            dcc_graph(id="humidity-graph", figure = plotHumidity()),

            dcc_interval(id="interval-component", interval = 250, n_intervals=0),
        ],
    )
end

run_server(app, "0.0.0.0", debug=true)


text:
#-------------------------------------------------------------------------------
# File: instances.jl
# Description: This files contains the instance sructure.
# Date: November 17, 2019
# Authors: Killian Fretaud, Rémi Garcia,
#-------------------------------------------------------------------------------

mutable struct Instance
    nb_jobs::Int # Nb of jobs
    nb_machines::Int # Nb of machines
    nb_robot::Int # Nb of position on the robot
    travel_times::Array{Int, 1} # t_j is the time taken to travel from j-1 to j
    processing_times::Array{Int, 2} # p_ij is the processing time for the job i on machine j

    Instance(N, M, T, P) = new(N, M, N*(M+1), T, P)
end

max_stars_repo_name: k-fretaud/Flowshop-robotic-cell

max_issues_repo_name: k-fretaud/Flowshop-robotic-cell

max_forks_repo_name: k-fretaud/Flowshop-robotic-cell


text:
using Pkg
include("JockeyEditor.jl")
function parse_bash_command(s)
    dims = split(s, " ")
    command = dims[1]
    if length(dims) != 1
        param = dims[2]
    end
    if command == "cd"
        cd(param)
    elseif command == "ls"
        readdir();
    elseif command == "pwd"
        pwd()
    elseif command == "rmdir"
        rmdir(param)
    elseif command == "mkdir"
        mkdir(param)
    elseif command == "source"
        Pkg.activate(param)
    else
        try
            run(Cmd(s))
        catch
            println("Could not determine command")
        end
    end
end

function install()
    Pkg.add("ReplMaker")
    home = homedir()
    startup = string(home, "/.julia/config/startup.jl")
    try
        touch(startup)
    catch
        mkdir(string(home,"/.julia/config"))
        touch(startup)
    end
    open(startup, "w") do io
        write(io, "# ========================\n")
        write(io, "#   Jockey Installation\n")
        write(io, "# ========================\n")
        write(io, "\n")
           write(io, "using Jockey\n")
           write(io, "using ReplMaker\n")
           write(io, "atreplinit() do repl\n")
               write(io, "try\n")
               write(io, "@eval using ReplMaker\n")
               write(io, "@async initrepl(\n")
               write(io, "parse_bash_command;\n")
               write(io, "prompt_text=\"jockey> \",\n")
               write(io, "prompt_color = :cyan,\n")
               write(io, "start_key = '[',\n")
               write(io, "mode_name = \"Jockey\");")
               write(io, "catch\n end\n end\n")

    end;

end


text:
@testset "Hubbard Model" begin

    # Type hierarchy
    @test HubbardModelAttractive <: HubbardModel
    @test HubbardModelRepulsive <: HubbardModel

    # constructors
    @test HubbardModel(1, 1, U =  1.0) isa HubbardModelRepulsive
    @test HubbardModel(1, 1, U = -1.0) isa HubbardModelAttractive

    for T in (HubbardModelAttractive, HubbardModelRepulsive)
        m = T(8, 1)
        @test length(m.l) == 8
        @test typeof(m) == T{Chain}

        l = SquareLattice(5)
        m = T(l = l);
        @test length(m.l) == 25
        @test typeof(m) == T{SquareLattice}

        m = T(l, t = 0.5);
        @test length(m.l) == 25
        @test m.t == 0.5
        @test typeof(m) == T{SquareLattice}

        d = Dict(:l => l, :U => 0.5)
        m = T(d)
        @test typeof(m) == T{SquareLattice}
        @test m.U == 0.5
    end
end


text:
"""
    polyinterp(t,y)

Construct a callable polynomial interpolant through the points in
vectors `t`,`y` using the barycentric interpolation formula.
"""
function polyinterp(t,y)
    n = length(t)-1
    C = (t[n+1]-t[1]) / 4           # scaling factor to ensure stability
    tc = t/C

    # Adding one node at a time, compute inverses of the weights.
    ω = ones(n+1)
    for m in 0:n-1
        d = tc[1:m+1] .- tc[m+2]    # vector of node differences
        @. ω[1:m+1] *= d            # update previous
        ω[m+2] = prod( -d )         # compute the new one
    end
    w = 1 ./ ω                      # go from inverses to weights

    # This function evaluates the interpolant at given x.
    p = function (x)
        terms = @. w / (x - t)
        if any(isinf.(terms))     # there was division by zero
            # return the node's data value
            idx = findfirst(x.==t)
            f = y[idx]
        else
            f = sum(y.*terms) / sum(terms)
        end
    end
    return p
end

"""
    triginterp(t,y)

Construct the trigonometric interpolant for the points defined by
vectors `t` and `y`.
"""
function triginterp(t,y)
    N = length(t)

    function τ(x)
        if x==0
            return 1.0
        else
            denom = isodd(N) ? N*sin(π*x/2) : N*tan(π*x/2)
            return sin(N*π*x/2)/denom
        end
    end

    return function (x)
        sum( y[k]*τ(x-t[k]) for k in eachindex(y) )
    end
end

"""
    ccint(f,n)

Perform Clenshaw-Curtis integration for the function `f` on `n`+1
nodes in [-1,1]. Returns the integral estimate and a vector of the
nodes used. Note: `n` must be even.
"""
function ccint(f,n)
    @assert iseven(n) "Value of `n` must be an even integer."
    # Find Chebyshev extreme nodes.
    θ = [ i*π/n for i in 0:n ]
    x = -cos.(θ)

    # Compute the C-C weights.
    c = similar(θ)
    c[[1,n+1]] .= 1/(n^2-1)
    s = sum( cos.(2k*θ[2:n])/(4k^2-1) for k in 1:n/2-1 )
    v = @. 1 - 2s - cos(n*θ[2:n])/(n^2-1)
    c[2:n] = 2v/n

    # Evaluate integrand and integral.
    I = dot(c,f.(x))   # vector inner product
    return I,x
end

"""
    glint(f,n)

Perform Gauss-Legendre integration for the function `f` on `n` nodes
in (-1,1). Returns the integral estimate and a vector of the nodes used.
"""
function glint(f,n)
    # Nodes and weights are found via a tridiagonal eigenvalue problem.
    β = @. 0.5/sqrt(1-(2*(1:n-1))^(-2))
    T = diagm(-1=>β,1=>β)
    λ,V = eigen(T)
    p = sortperm(λ)
    x = λ[p]               # nodes
    c = @. 2V[1,p]^2       # weights

    # Evaluate the integrand and compute the integral.
    I = dot(c,f.(x))      # vector inner product
    return I,x
end

"""
    intinf(f,tol)

Perform adaptive doubly-exponential integration of function `f`
over (-Inf,Inf), with error tolerance `tol`. Returns the integral
estimate and a vector of the nodes used.
"""
function intinf(f,tol)
    x = t -> sinh(sinh(t))
    dx_dt = t -> cosh(t)*cosh(sinh(t))
    g = t -> f(x(t))*dx_dt(t)

    # Find where to truncate the integration interval.
    M = 3
    while (abs(g(-M)) > tol/100) || (abs(g(M)) > tol/100)
        M += 0.5
        if isinf(x(M))
            @warn "Function may not decay fast enough."
            M -= 0.5
            break
        end
    end

    I,t = intadapt(g,-M,M,tol)
	return I,x.(t)
end

"""
    intsing(f,tol)

Adaptively integrate function `f` over (0,1), where `f` may be
singular at zero, with error tolerance `tol`. Returns the
integral estimate and a vector of the nodes used.
"""
function intsing(f,tol)
    x = t -> 2/(1+exp(2sinh(t)))
	dx_dt = t -> cosh(t)/cosh(sinh(t))^2
	g = t -> f(x(t))*dx_dt(t)

    # Find where to truncate the integration interval.
    M = 3
    while abs(g(M)) > tol/100
        M += 0.5
        if iszero(x(M))
            @warn "Function may grow too rapidly."
            M -= 0.5
            break
        end
    end

    I,t = intadapt(g,0,M,tol)
	return I,x.(t)
end


text:
function arterialpressure!(system::LegSystem,n::Int64)
    for i = 1:length(system.branches.ID)
        # for j = 1:system.solverparams.JL
        #     if system.branches.A[i][j,n+2] < 0
        #         error("Area negative in pressure calculation for artery $i, node $j,
        #             time step $n. Area = $(system.branches.A[i][j,n+2]) m2.")
        #     end
        # end
        system.branches.P[i][:,n+2] .= system.branches.beta[i][end].*
            (system.branches.A[i][:,n+2].^0.5.-
            sqrt.(system.branches.A0[i][end]))./mmHgToPa;
    end
end

max_stars_repo_name: dcanuto/leg-hemodynamics

max_issues_repo_name: dcanuto/leg-hemodynamics

max_forks_repo_name: dcanuto/leg-hemodynamics


text:
using TensorTrains

L=10
h,V = hubbard_1D(L)
H_tto = hV_to_mpo(h,V,tol=1e-5)
ψ_0 = half_filling(L)
E_als,ψ_tt = als_eigsolv(H_tto,ψ_0,sweep_schedule=[2,5],rmax_schedule=[16,32])
E_mals,ψ_tt,r = mals_eigsolv(H_tto,ψ_tt,tol=1e-6,sweep_schedule=[2,5],rmax_schedule=[32,48])

max_stars_repo_name: msdupuy/Tensor-Train-Julia

max_issues_repo_name: msdupuy/Tensor-Train-Julia

max_forks_repo_name: msdupuy/Tensor-Train-Julia


text:
using JuliaInterpreter
using JuliaInterpreter: enter_call_expr
using Test

function summer(A)
    s = zero(eltype(A))
    for a in A
        s += a
    end
    return s
end

A = [0.12, -.99]
frame = JuliaInterpreter.enter_call(summer, A)
frame2 = JuliaInterpreter.enter_call(summer, A)
@test summer(A) == something(runframe(frame)) == something(runstack(frame2))

A = rand(1000)
@test @interpret(sum(A)) ≈ sum(A)  # note: the compiler can leave things in registers to increase accuracy, doesn't happen with interpreted
fapply() = (Core.apply_type)(Base.NamedTuple, (), Tuple{})
@test @interpret(fapply()) == fapply()
function fbc()
    bc = Broadcast.broadcasted(CartesianIndex, 6, [1, 2, 3])
    copy(bc)
end
@test @interpret(fbc()) == fbc()
@test @interpret(repr("hi")) == repr("hi")  # this tests kwargs and @generated functions

fkw(x::Int8; y=0, z="hello") = y
@test @interpret(fkw(Int8(1); y=22, z="world")) == fkw(Int8(1); y=22, z="world")

# Throwing exceptions across frames
function f_exc_inner()
    error("inner")
end

f_exc_inner2() = f_exc_inner()

const caught = Ref(false)
function f_exc_outer1()
    try
        f_exc_inner()
    catch err    # with an explicit err capture
        caught[] = true
        rethrow(err)
    end
end

function f_exc_outer2()
    try
        f_exc_inner()
    catch        # implicit err capture
        caught[] = true
        rethrow()
    end
end

function f_exc_outer3(f)
    try
        f()
    catch err
        return err
    end
end

@test !caught[]
ret = @interpret f_exc_outer3(f_exc_outer1)
@test ret == ErrorException("inner")
@test caught[]

caught[] = false
ret = @interpret f_exc_outer3(f_exc_outer2)
@test ret == ErrorException("inner")
@test caught[]

caught[] = false
ret = @interpret f_exc_outer3(f_exc_inner2)
@test ret == ErrorException("inner")
@test !caught[]

stc = try f_exc_outer1() catch
    stacktrace(catch_backtrace())
end
sti = try @interpret(f_exc_outer1()) catch
    stacktrace(catch_backtrace())
end
@test_broken stc == sti

# issue #3
@test @interpret(joinpath("/home/julia/base", "sysimg.jl")) == "/home/julia/base/sysimg.jl"
@test @interpret(10.0^4) == 10.0^4
# issue #6
@test @interpret(Array.body.body.name) === Array.body.body.name
@test @interpret(Vararg.body.body.name) === Vararg.body.body.name
frame = JuliaInterpreter.prepare_toplevel(Main, :(Vararg.body.body.name))
@test JuliaInterpreter.finish_and_return!(JuliaStackFrame[], frame, true) === Vararg.body.body.name

# issue #8
ex = quote
    if sizeof(JLOptions) === ccall(:jl_sizeof_jl_options, Int, ())
    else
        ccall(:jl_throw, Cvoid, (Any,), "Option structure mismatch")
    end
end
frame = JuliaInterpreter.prepare_toplevel(Base, ex)
JuliaInterpreter.finish_and_return!(JuliaStackFrame[], frame, true)

# ccall with two Symbols
ex = quote
    @testset "Some tests" begin
       @test 2 > 1
    end
end
frame = JuliaInterpreter.prepare_toplevel(Main, ex)
JuliaInterpreter.finish_and_return!(JuliaStackFrame[], frame, true)

@test @interpret Base.Math.DoubleFloat64(-0.5707963267948967, 4.9789962508669555e-17).hi ≈ -0.5707963267948967

# ccall with cfunction
fcfun(x::Int, y::Int) = 1
ex = quote   # in lowered code, cf is a Symbol
    cf = @eval @cfunction(fcfun, Int, (Int, Int))
    ccall(cf, Int, (Int, Int), 1, 2)
end
frame = JuliaInterpreter.prepare_toplevel(Main, ex)
@test JuliaInterpreter.finish_and_return!(JuliaStackFrame[], frame, true) == 1
ex = quote
    let   # in lowered code, cf is a SlotNumber
        cf = @eval @cfunction(fcfun, Int, (Int, Int))
        ccall(cf, Int, (Int, Int), 1, 2)
    end
end
frame = JuliaInterpreter.prepare_toplevel(Main, ex)
@test JuliaInterpreter.finish_and_return!(JuliaStackFrame[], frame, true) == 1

# From Julia's test/ambiguous.jl. This tests whether we renumber :enter statements correctly.
ambig(x, y) = 1
ambig(x::Integer, y) = 2
ambig(x, y::Integer) = 3
ambig(x::Int, y::Int) = 4
ambig(x::Number, y) = 5
ex = quote
    let
        cf = @eval @cfunction(ambig, Int, (UInt8, Int))
        @test_throws(MethodError, ccall(cf, Int, (UInt8, Int), 1, 2))
    end
end
frame = JuliaInterpreter.prepare_toplevel(Main, ex)
JuliaInterpreter.finish_and_return!(JuliaStackFrame[], frame, true)

# Core.Compiler
ex = quote
    length(code_typed(fcfun, (Int, Int)))
end
frame = JuliaInterpreter.prepare_toplevel(Main, ex)
@test JuliaInterpreter.finish_and_return!(JuliaStackFrame[], frame, true) == 1


text:
using ProblemParser
using Utils

using Chain

const file = find_input(@__FILE__)
const input = parse(Rectangular(Convert(Int8)), slurp(file))

const NEIGHBORS = CartesianIndices((-1:1, -1:1))

function flash(input, locations)
    length(locations) == 0 && return 0
    cascade = Set{CartesianIndex}()
    for loc in locations
        for neighbor in NEIGHBORS .+ [loc]
            input[neighbor] += 1
            if input[neighbor] == 10
                push!(cascade, neighbor)
            end
        end
    end
    input[begin,:] .= 0
    input[end,:] .= 0
    input[:,begin] .= 0
    input[:,end] .= 0
    return length(locations) + flash(input, cascade)
end

struct Aoc2111
    matrix::Matrix{Int8}
    inner::CartesianIndices
end
Aoc2111(input) = Aoc2111(input,
    CartesianIndices((firstindex(input, 1) + 1:lastindex(input, 2) + 1, firstindex(input, 1) + 1:lastindex(input, 2) + 1))
)

function Base.iterate(aoc::Aoc2111, state=boundaryconditions(aoc.matrix, 0))
    state[aoc.inner] .+= 1
    flashers = findall(>(9), state)
    flashes = flash(state, flashers)
    state[state .> 9] .= 0
    ((grid=view(state, aoc.inner), flashes=flashes), state)
end

Base.IteratorSize(::Aoc2111) = Base.IsInfinite()
Base.eltype(::Aoc2111) = @NamedTuple{grid::AbstractMatrix{Int8}, flashes::Int}

function part1(input, steps)
    @chain input begin
        Aoc2111
        Iterators.take(steps)
        map(tup -> tup.flashes, _)
        sum
    end
end

assertequal(part1(input, 100), 1608)

function part2(input)
    @chain input begin
        Aoc2111
        enumerate
        Iterators.filter(_) do (i,state)
            all(state.grid .== 0)
        end
        first
        first
    end
end

assertequal(part2(input), 214)

max_stars_repo_name: CmdQ/AoC


text:
function dice_result(dice)
    if dice[1] == dice[2] == dice[3]
        if dice[1] == 1
            return 7
        end
        return dice[1]
    end
    sorted_dice = sort([dice...])
    if sorted_dice[1] == 1 && sorted_dice[2] == 2 && sorted_dice[3] == 4
        if dice[1] == 2
            return 10
        end
        return 8
    end
    if (sorted_dice[2] == sorted_dice[1] + 1 &&
        sorted_dice[3] == sorted_dice[2] + 1)
        return 3
    end
    1
end

type Reroll
    score::Float64
    dice::(Int64, Int64)
end

DICE_SCORE = Dict{Array{Float64}, Reroll}()

for i = 1:6
    for j = 1:6
        for k = 1:6
            DICE_SCORE[[i, j, k]] = Reroll(dice_result((i, j, k)), (0, 0))
        end
    end
end

function score_cdf(score)
    num_under = count(x -> x <= score, values(DICE_SCORE))
    return num_under / length(DICE_SCORE)
end

function new_roll(dice, i, val)
    roll = [dice...]
    roll[i] = val
    roll
end

function max_expected_value(dice, dice_score, bias = 0)
    current_score = DICE_SCORE[dice]
    max_score = current_score.score + bias
    best_i = 0
    best_j = 0
    for i = 1:3
        for j = 1:3
            if i == j
                # rolls = [new_roll(dice, i, k).score for k = 1:6]
                # println(rolls)
                results = [dice_score[new_roll(dice, i, k)].score for k = 1:6]
                # println(results)
                score = sum(results) / 6
            else
                score = sum([dice_score[new_roll(new_roll(dice, i, x), j, y)].score
                             for x = 1:6, y = 1:6]) / 36
            end
            if score > max_score
                max_score = score
                best_i = i
                best_j = j
            end
        end
    end
    return Reroll(max_score, (best_i, best_j))
end

BEST_ONE_STEP = Dict{Array{Int64}, Reroll}()
BEST_TWO_STEP = Dict{Array{Int64}, Reroll}()

for i = 1:6
    for j = 1:6
        for k = 1:6
            BEST_ONE_STEP[[i, j, k]] = max_expected_value([i, j, k], DICE_SCORE)
        end
    end
end

for i = 1:6
    for j = 1:6
        for k = 1:6
            BEST_TWO_STEP[[i, j, k]] = max_expected_value([i, j, k], BEST_ONE_STEP, 1)
        end
    end
end

println(BEST_TWO_STEP)

# println(max_expected_value([4,2,5], DICE_SCORE))
# println(score_cdf(10))

println(BEST_TWO_STEP[[2,2,2]])


text:
# This file is auto-generated by AWSMetadata.jl
using AWS
using AWS.AWSServices: ec2_instance_connect
using AWS.Compat
using AWS.UUIDs

"""
    send_serial_console_sshpublic_key(instance_id, sshpublic_key)
    send_serial_console_sshpublic_key(instance_id, sshpublic_key, params::Dict{String,<:Any})

Pushes an SSH public key to the specified EC2 instance. The key remains for 60 seconds,
which gives you 60 seconds to establish a serial console connection to the instance using
SSH. For more information, see EC2 Serial Console in the Amazon EC2 User Guide.

# Arguments
- `instance_id`: The ID of the EC2 instance.
- `sshpublic_key`: The public key material. To use the public key, you must have the
  matching private key. For information about the supported key formats and lengths, see
  Requirements for key pairs in the Amazon EC2 User Guide.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"SerialPort"`: The serial port of the EC2 instance. Currently only port 0 is supported.
  Default: 0
"""
function send_serial_console_sshpublic_key(
    InstanceId, SSHPublicKey; aws_config::AbstractAWSConfig=global_aws_config()
)
    return ec2_instance_connect(
        "SendSerialConsoleSSHPublicKey",
        Dict{String,Any}("InstanceId" => InstanceId, "SSHPublicKey" => SSHPublicKey);
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function send_serial_console_sshpublic_key(
    InstanceId,
    SSHPublicKey,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return ec2_instance_connect(
        "SendSerialConsoleSSHPublicKey",
        Dict{String,Any}(
            mergewith(
                _merge,
                Dict{String,Any}(
                    "InstanceId" => InstanceId, "SSHPublicKey" => SSHPublicKey
                ),
                params,
            ),
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    send_sshpublic_key(availability_zone, instance_id, instance_osuser, sshpublic_key)
    send_sshpublic_key(availability_zone, instance_id, instance_osuser, sshpublic_key, params::Dict{String,<:Any})

Pushes an SSH public key to the specified EC2 instance for use by the specified user. The
key remains for 60 seconds. For more information, see Connect to your Linux instance using
EC2 Instance Connect in the Amazon EC2 User Guide.

# Arguments
- `availability_zone`: The Availability Zone in which the EC2 instance was launched.
- `instance_id`: The ID of the EC2 instance.
- `instance_osuser`: The OS user on the EC2 instance for whom the key can be used to
  authenticate.
- `sshpublic_key`: The public key material. To use the public key, you must have the
  matching private key.

"""
function send_sshpublic_key(
    AvailabilityZone,
    InstanceId,
    InstanceOSUser,
    SSHPublicKey;
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return ec2_instance_connect(
        "SendSSHPublicKey",
        Dict{String,Any}(
            "AvailabilityZone" => AvailabilityZone,
            "InstanceId" => InstanceId,
            "InstanceOSUser" => InstanceOSUser,
            "SSHPublicKey" => SSHPublicKey,
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function send_sshpublic_key(
    AvailabilityZone,
    InstanceId,
    InstanceOSUser,
    SSHPublicKey,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return ec2_instance_connect(
        "SendSSHPublicKey",
        Dict{String,Any}(
            mergewith(
                _merge,
                Dict{String,Any}(
                    "AvailabilityZone" => AvailabilityZone,
                    "InstanceId" => InstanceId,
                    "InstanceOSUser" => InstanceOSUser,
                    "SSHPublicKey" => SSHPublicKey,
                ),
                params,
            ),
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end


text:
using CorpusLoaders
using Test
using Base.Iterators
using MultiResolutionIterators

@testset "Basic use" begin
    wk_gen = load(Senseval3())
    docs = collect(take(wk_gen,10));

    words = consolidate(flatten_levels(docs, (!lvls)(Senseval3, :word)))
    @test length(words) > length(docs)
    @test length(words) > sum(length.(docs))
    @test typeof(words) == Vector{CorpusLoaders.TaggedWord}

    plain_words = word.(words)
    @test typeof(plain_words) <: Vector{String}
end

@testset "Index" begin
    corp = load(Senseval3())
    index = Dict{String, Vector{Vector{String}}}()
    for sentence in flatten_levels(corp, (!lvls)(Senseval3, :sent, :word))
        multisense_words = filter(x->x isa SenseAnnotatedWord,  sentence)
        context_words = word.(sentence)
        for target in multisense_words
            word_sense = sensekey(target)
            uses = get!(Vector{Vector{String}}, index, word_sense)
            push!(uses, context_words)
        end
    end

    index

    index["answer%1:04:00::"]
    index["answer%1:10:00::"]
    index["answer%2:32:00::"]
end


text:
# Getting Started -
## Developing a New ACE Model with ACE1.jl
## without IPFitting

# First import ACE1.jl
using ACE1, JuLIP, IPFitting
using LinearAlgebra: norm

# [1] create an ACE1 basis for a Si potential
# rpi = rotation and permutation invariant
r0 = rnn(:Si)
basis = rpi_basis(;
      species = :Si,
      N = 3,                        # correlation order = body-order - 1
      maxdeg = 10,                  # polynomial degree
      r0 = r0,                      # estimate for NN distance
      D = SparsePSHDegree(; wL=1.3, csp=1.0),
      rin = 0.65*r0, rcut = 5.0,    # domain for radial basis (cf documentation)
      pin = 0)
@show length(basis)

# Remarks:
#  SparsePSHDegree(; wL=1.3, csp=1.0) selects a specific sparse subset of the
#  infinite lattice of possible ACE basis functions. See the documentation
#  for more details:
#     ?SparsePSHDegree

# [2] create a random training set (usually should load one e.g. from xyz)

function gen_dat()
   sw = StillingerWeber()
   at = rattle!(bulk(:Si, cubic=true) * 3, 0.3)
   return Dat(at, "diax3"; E = energy(sw, at), F = forces(sw, at) )
end

train = [ gen_dat() for _=1:50 ]

# [3] estimate the parameters

# assemble the least squares system; using "" for the filename means that
# it will not be saved to disk; this evaluates the basis on the training
# but doesn't really assemble the lsq system yet. This allows one to experimate
# with different fitting parameters
dB = LsqDB("", basis, train)

# specify weights for lsq system, and a reference potential.
weights = Dict("default" => Dict("E" => 15.0, "F" => 1.0 , "V" => 1.0 ))

# now assemble and solve the LSQ system
IP, lsqinfo = lsqfit(dB; weights = weights, asmerrs = true)

# [4] Perform some tests, here we do it manually just to demonstrate usage
#     of the potential
# note IP is now an ACE model whose energy, force, virial etc can now be
# evaluated using JuLIP.

test =  [ gen_dat() for _=1:20 ]
Etest = [ dat.D["E"][1] for dat in test ]
Emodel = [ energy(IP, dat.at) for dat in test ]
rmse_E = norm(Etest - Emodel) / sqrt(length(test))


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct NetworkSecurityGroupPropertiesFormat <: SwaggerModel
    securityRules::Any # spec type: Union{ Nothing, Vector{SecurityRule} } # spec name: securityRules
    defaultSecurityRules::Any # spec type: Union{ Nothing, Vector{SecurityRule} } # spec name: defaultSecurityRules
    networkInterfaces::Any # spec type: Union{ Nothing, Vector{NetworkInterface} } # spec name: networkInterfaces
    subnets::Any # spec type: Union{ Nothing, Vector{Subnet} } # spec name: subnets
    flowLogs::Any # spec type: Union{ Nothing, Vector{FlowLog} } # spec name: flowLogs
    resourceGuid::Any # spec type: Union{ Nothing, String } # spec name: resourceGuid
    provisioningState::Any # spec type: Union{ Nothing, ProvisioningState } # spec name: provisioningState

    function NetworkSecurityGroupPropertiesFormat(;securityRules=nothing, defaultSecurityRules=nothing, networkInterfaces=nothing, subnets=nothing, flowLogs=nothing, resourceGuid=nothing, provisioningState=nothing)
        o = new()
        validate_property(NetworkSecurityGroupPropertiesFormat, Symbol("securityRules"), securityRules)
        setfield!(o, Symbol("securityRules"), securityRules)
        validate_property(NetworkSecurityGroupPropertiesFormat, Symbol("defaultSecurityRules"), defaultSecurityRules)
        setfield!(o, Symbol("defaultSecurityRules"), defaultSecurityRules)
        validate_property(NetworkSecurityGroupPropertiesFormat, Symbol("networkInterfaces"), networkInterfaces)
        setfield!(o, Symbol("networkInterfaces"), networkInterfaces)
        validate_property(NetworkSecurityGroupPropertiesFormat, Symbol("subnets"), subnets)
        setfield!(o, Symbol("subnets"), subnets)
        validate_property(NetworkSecurityGroupPropertiesFormat, Symbol("flowLogs"), flowLogs)
        setfield!(o, Symbol("flowLogs"), flowLogs)
        validate_property(NetworkSecurityGroupPropertiesFormat, Symbol("resourceGuid"), resourceGuid)
        setfield!(o, Symbol("resourceGuid"), resourceGuid)
        validate_property(NetworkSecurityGroupPropertiesFormat, Symbol("provisioningState"), provisioningState)
        setfield!(o, Symbol("provisioningState"), provisioningState)
        o
    end
end # type NetworkSecurityGroupPropertiesFormat

const _property_map_NetworkSecurityGroupPropertiesFormat = Dict{Symbol,Symbol}(Symbol("securityRules")=>Symbol("securityRules"), Symbol("defaultSecurityRules")=>Symbol("defaultSecurityRules"), Symbol("networkInterfaces")=>Symbol("networkInterfaces"), Symbol("subnets")=>Symbol("subnets"), Symbol("flowLogs")=>Symbol("flowLogs"), Symbol("resourceGuid")=>Symbol("resourceGuid"), Symbol("provisioningState")=>Symbol("provisioningState"))
const _property_types_NetworkSecurityGroupPropertiesFormat = Dict{Symbol,String}(Symbol("securityRules")=>"Vector{SecurityRule}", Symbol("defaultSecurityRules")=>"Vector{SecurityRule}", Symbol("networkInterfaces")=>"Vector{NetworkInterface}", Symbol("subnets")=>"Vector{Subnet}", Symbol("flowLogs")=>"Vector{FlowLog}", Symbol("resourceGuid")=>"String", Symbol("provisioningState")=>"ProvisioningState")
Base.propertynames(::Type{ NetworkSecurityGroupPropertiesFormat }) = collect(keys(_property_map_NetworkSecurityGroupPropertiesFormat))
Swagger.property_type(::Type{ NetworkSecurityGroupPropertiesFormat }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_NetworkSecurityGroupPropertiesFormat[name]))}
Swagger.field_name(::Type{ NetworkSecurityGroupPropertiesFormat }, property_name::Symbol) =  _property_map_NetworkSecurityGroupPropertiesFormat[property_name]

function check_required(o::NetworkSecurityGroupPropertiesFormat)
    true
end

function validate_property(::Type{ NetworkSecurityGroupPropertiesFormat }, name::Symbol, val)
end


text:
maxval = 1.0  # maximum absolute value of perturbations

# Large eddy size (m).
# Above this scale, the spectrum is
# white.  If you don't want to use
# this, make it larger than
# max(ni,nj)*dx.
large_eddy_scale = 1000000000.0

function ring_filter(Nx, Ny, Δx, Δy, L₁)
    # center point indices
    i_c = Nx/2 + 1
    j_c = Ny/2 + 1

    # Compute wavenumbers
    κ = zeros(Nx, Ny)
    for j in 1:Ny, i in 1:Nx
        κ[i, j] = √( (2π*(i-i_c) / (Nx*Δx))^2 + (2π*(j-j_c) / (Ny*Δy))^2 )
    end

    # Generate spectrum
    A = 1
    c = zeros(ComplexF64, Nx, Ny)
    for j in 1:Ny, i in 1:Nx
        r1 = rand()
        r2 = rand()
        c[i, j] = A * exp(-2π*im*Nx*r1) * exp(-2π*im*Ny*r2)
    end

    i_min = j_min = 0
    v_min = Inf
    for j in 1:Ny, i in 1:Nx
        val = abs(κ[i, j] - 2π/L₁)
        if val < v_min
            v_min = val
            i_min = i
            j_min = j
        end
    end

    threshold = max(
        abs(κ[i_min, j_min+1] - 2π/L₁),
        abs(κ[i_min, j_min-1] - 2π/L₁),
        abs(κ[i_min+1, j_min] - 2π/L₁),
        abs(κ[i_min-1, j_min] - 2π/L₁)
    )

    for j in 1:Ny, i in 1:Nx
        if abs(κ[i, j] - 2π/L₁) > threshold
            c[i, j] = 0
        end
    end

    spectrum_2d = copy(c)

    c = circshift(c, (Nx/2, Ny/2))
    KE = ifft(c)
    c = circshift(fft(KE), (Nx/2, Ny/2))

    Nk = min(Nx, Ny)
    sp = zeros(Nk)
    np = zeros(Nk)

    for j in 1:Ny, i in 1:Nx
        if Nx > Ny
            ik = √( ((i-i_c)*Ny/Nx)^2 + (j-j_c)^2 )
            ik = round(Int, ik)
        elseif Ny > Nx
            ik = √( (i-i_c)^2 + ((j-j_c)*Nx/Ny)^2 )
            ik = round(Int, ik)
        else
            ik = √( (i-i_c)^2 + (j-j_c)^2 )
            ik = round(Int, ik)
        end

        if ik > 0 && ik < Nk/2
            sp[ik] += abs(ccc[i, j])^2
            np[ik] += 1
        end
    end

    spectrum_1d = copy(sp)

    return abs2.(KE), spectrum_2d, spectrum_1d
end

Nx = Ny = 192
Δx = Δy = 125
L₁ = 1000

x = range(0, Nx*Δx, length=Nx)
y = range(0, Ny*Δy, length=Ny)

KE, s2, s1 = ring_filter(Nx, Ny, Δx, Δy, L₁)

l = @layout [
    a{0.7w} [b{0.5h}
             c{0.5h}]
]

pert_plot = contourf(x/1000, y/1000, KE, color=:dense, linewidth=0, xticks=0:4:24, xlabel="x (km)", yticks=0:4:24, ylabel="y (km)", title="Perturbation field", colorbar=nothing)
s2_plot = contourf(abs2.(s2), color=:dense, linewidth=0, xticks=[], yticks=[], colorbar=nothing, framestyle=:box, title="2D spectrum", xlabel="kx", ylabel="ky")
s1_plot = plot(log10.(s1), linewidth=2, label="", xlims=(0, 96), xticks=[], yticks=[], xlabel="k", ylabel="log E(k)", title="1D spectrum")
plot(pert_plot, s2_plot, s1_plot, layout=l)

Nx = Ny = 192
Δx = Δy = 1/192
L₁ = 0.3

x = range(0, Nx*Δx, length=Nx)
y = range(0, Ny*Δy, length=Ny)

KE, s2, s1 = ring_filter(Nx, Ny, Δx, Δy, L₁)

l = @layout [
    a{0.7w} [b{0.5h}
             c{0.5h}]
]

pert_plot = contourf(x, y, KE, color=:dense, linewidth=0, xlabel="x", ylabel="y", title="Perturbation field", colorbar=nothing)
s2_plot = contourf(abs2.(s2), color=:dense, linewidth=0, xticks=[], yticks=[], colorbar=nothing, framestyle=:box, title="2D spectrum", xlabel="kx", ylabel="ky")
s1_plot = plot(log10.(s1), linewidth=2, label="", xlims=(0, 96), xticks=[], yticks=[], xlabel="k", ylabel="log E(k)", title="1D spectrum")
plot(pert_plot, s2_plot, s1_plot, layout=l)


text:
# LLVM IR optimization

function addTargetPasses!(pm, tm, triple)
    add_library_info!(pm, triple)
    add_transform_info!(pm, tm)
end

# Based on Julia's optimization pipeline, minus the SLP and loop vectorizers.
function addOptimizationPasses!(pm, opt_level=2)
    # compare with the using Julia's optimization pipeline directly:
    #ccall(:jl_add_optimization_passes, Cvoid,
    #      (LLVM.API.LLVMPassManagerRef, Cint, Cint),
    #      pm, opt_level, #=lower_intrinsics=# 0)
    #return

    # NOTE: LLVM 12 disabled the hoisting of common instruction
    #       before loop vectorization (https://reviews.llvm.org/D84108).
    #
    #       This is re-enabled with calls to cfg_simplify here,
    #       to merge allocations and sometimes eliminate them,
    #       since AllocOpt does not handle PhiNodes.
    #       Enable this instruction hoisting because of this and Union benchmarks.

    constant_merge!(pm)

    if opt_level < 2
        cpu_features!(pm)
        if opt_level == 1
            instruction_simplify!(pm)
        end
        if LLVM.version() >= v"12"
            cfgsimplification!(pm; hoist_common_insts=true)
        else
            cfgsimplification!(pm)
        end
        if opt_level == 1
            scalar_repl_aggregates!(pm)
            instruction_combining!(pm)
            early_cse!(pm)
            # maybe add GVN?
            # also try GVNHoist and GVNSink
        end
        mem_cpy_opt!(pm)
        always_inliner!(pm) # Respect always_inline
        lower_simdloop!(pm) # Annotate loop marked with "loopinfo" as LLVM parallel loop
        return
    end

    propagate_julia_addrsp!(pm)
    scoped_no_alias_aa!(pm)
    type_based_alias_analysis!(pm)
    if opt_level >= 3
        basic_alias_analysis!(pm)
    end
    if LLVM.version() >= v"12"
        cfgsimplification!(pm; hoist_common_insts=true)
    else
        cfgsimplification!(pm)
    end
    dce!(pm)
    scalar_repl_aggregates!(pm)

    #mem_cpy_opt!(pm)

    always_inliner!(pm) # Respect always_inline

    # Running `memcpyopt` between this and `sroa` seems to give `sroa` a hard
    # time merging the `alloca` for the unboxed data and the `alloca` created by
    # the `alloc_opt` pass.

    alloc_opt!(pm)
    # consider AggressiveInstCombinePass at optlevel > 2
    instruction_combining!(pm)
    if LLVM.version() >= v"12"
        cfgsimplification!(pm; hoist_common_insts=true)
    else
        cfgsimplification!(pm)
    end
    cpu_features!(pm)
    scalar_repl_aggregates!(pm)
    instruction_simplify!(pm)
    jump_threading!(pm)
    correlated_value_propagation!(pm)

    reassociate!(pm)

    early_cse!(pm)

    # Load forwarding above can expose allocations that aren't actually used
    # remove those before optimizing loops.
    alloc_opt!(pm)
    loop_rotate!(pm)
    # moving IndVarSimplify here prevented removing the loop in perf_sumcartesian(10:-1:1)
    loop_idiom!(pm)

    # LoopRotate strips metadata from terminator, so run LowerSIMD afterwards
    lower_simdloop!(pm) # Annotate loop marked with "loopinfo" as LLVM parallel loop
    licm!(pm)
    julia_licm!(pm)
    loop_unswitch!(pm)
    licm!(pm)
    julia_licm!(pm)
    inductive_range_check_elimination!(pm)
    # Subsequent passes not stripping metadata from terminator
    instruction_simplify!(pm)
    ind_var_simplify!(pm)
    loop_deletion!(pm)
    loop_unroll!(pm) # TODO: in Julia createSimpleLoopUnroll

    # Run our own SROA on heap objects before LLVM's
    alloc_opt!(pm)
    # Re-run SROA after loop-unrolling (useful for small loops that operate,
    # over the structure of an aggregate)
    scalar_repl_aggregates!(pm)
    # might not be necessary:
    instruction_simplify!(pm)

    gvn!(pm)
    mem_cpy_opt!(pm)
    sccp!(pm)

    # Run instcombine after redundancy elimination to exploit opportunities
    # opened up by them.
    # This needs to be InstCombine instead of InstSimplify to allow
    # loops over Union-typed arrays to vectorize.
    instruction_combining!(pm)
    jump_threading!(pm)
    correlated_value_propagation!(pm)
    dead_store_elimination!(pm)

    # More dead allocation (store) deletion before loop optimization
    # consider removing this:
    alloc_opt!(pm)
    # see if all of the constant folding has exposed more loops
    # to simplification and deletion
    # this helps significantly with cleaning up iteration
    cfgsimplification!(pm)  # See note above, don't hoist instructions before LV
    loop_deletion!(pm)
    instruction_combining!(pm)
    loop_vectorize!(pm)
    loop_load_elimination!(pm)
    # Cleanup after LV pass
    if LLVM.version() >= v"12"
        cfgsimplification!(pm; # Aggressive CFG simplification
            forward_switch_cond_to_phi=true,
            convert_switch_to_lookup_table=true,
            need_canonical_loop=true,
            hoist_common_insts=true,
            sink_common_insts=true) # FIXME: Causes assertion in llvm-late-lowering
    else
        cfgsimplification!(pm)
    end

    aggressive_dce!(pm)
end

function optimize!(@nospecialize(job::CompilerJob), mod::LLVM.Module)
    triple = llvm_triple(job.target)
    tm = llvm_machine(job.target)

    global current_job
    current_job = job

    ModulePassManager() do pm
        addTargetPasses!(pm, tm, triple)
        addOptimizationPasses!(pm)
        run!(pm, mod)
    end

    # NOTE: we need to use multiple distinct pass managers to force pass ordering;
    #       intrinsics should never get lowered before Julia has optimized them.
    # XXX: why doesn't the barrier noop pass work here?

    # lower intrinsics
    ModulePassManager() do pm
        addTargetPasses!(pm, tm, triple)

        add!(pm, FunctionPass("LowerGCFrame", lower_gc_frame!))

        if job.source.kernel
            # GC lowering is the last pass that may introduce calls to the runtime library,
            # and thus additional uses of the kernel state intrinsic.
            add!(pm, FunctionPass("LowerKernelState", lower_kernel_state!))
            add!(pm, ModulePass("CleanupKernelState", cleanup_kernel_state!))
        end

        # remove dead uses of ptls
        aggressive_dce!(pm)
        add!(pm, ModulePass("LowerPTLS", lower_ptls!))

        # the Julia GC lowering pass also has some clean-up that is required
        late_lower_gc_frame!(pm)

        remove_ni!(pm)
        remove_julia_addrspaces!(pm)

        # Julia's operand bundles confuse the inliner, so repeat here now they are gone.
        # FIXME: we should fix the inliner so that inlined code gets optimized early-on
        always_inliner!(pm)

        # some of Julia's optimization passes happen _after_ lowering intrinsics
        combine_mul_add!(pm)
        div_rem_pairs!(pm)

        run!(pm, mod)
    end

    # target-specific optimizations
    optimize_module!(job, mod)

    # we compile a module containing the entire call graph,
    # so perform some interprocedural optimizations.
    #
    # for some reason, these passes need to be distinct from the regular optimization chain,
    # or certain values (such as the constant arrays used to populare llvm.compiler.user ad
    # part of the LateLowerGCFrame pass) aren't collected properly.
    #
    # these might not always be safe, as Julia's IR metadata isn't designed for IPO.
    ModulePassManager() do pm
        addTargetPasses!(pm, tm, triple)

        # - remove unused kernel state arguments
        # - simplify function calls that don't use the returned value
        dead_arg_elimination!(pm)

        run!(pm, mod)
    end

    # compare to Clang by using the pass manager builder APIs:
    #LLVM.clopts("-print-after-all", "-filter-print-funcs=$(LLVM.name(entry))")
    #ModulePassManager() do pm
    #    addTargetPasses!(pm, tm, triple)
    #    PassManagerBuilder() do pmb
    #        optlevel!(pmb, 2)
    #        populate!(pm, pmb)
    #    end
    #    run!(pm, mod)
    #end

    return
end

## lowering intrinsics
cpu_features!(pm::PassManager) = add!(pm, ModulePass("LowerCPUFeatures", cpu_features!))
function cpu_features!(mod::LLVM.Module)
    job = current_job::CompilerJob
    ctx = context(mod)
    changed = false

    argtyps = Dict(
        "f32" => Float32,
        "f64" => Float64,
    )

    # have_fma
    for f in functions(mod)
        ft = eltype(llvmtype(f))
        fn = LLVM.name(f)
        startswith(fn, "julia.cpu.have_fma.") || continue
        typnam = fn[20:end]

        # determine whether this back-end supports FMA on this type
        has_fma = if haskey(argtyps, typnam)
            typ = argtyps[typnam]
            have_fma(job.target, typ)
        else
            # warn?
            false
        end
        has_fma = ConstantInt(return_type(ft), has_fma)

        # substitute all uses of the intrinsic with a constant
        materialized = LLVM.Value[]
        for use in uses(f)
            val = user(use)
            replace_uses!(val, has_fma)
            push!(materialized, val)
        end

        # remove the intrinsic and its uses
        for val in materialized
            @assert isempty(uses(val))
            unsafe_delete!(LLVM.parent(val), val)
        end
        @assert isempty(uses(f))
        unsafe_delete!(mod, f)
    end

    return changed
end

# lower object allocations to to PTX malloc
#
# this is a PoC implementation that is very simple: allocate, and never free. it also runs
# _before_ Julia's GC lowering passes, so we don't get to use the results of its analyses.
# when we ever implement a more potent GC, we will need those results, but the relevant pass
# is currently very architecture/CPU specific: hard-coded pool sizes, TLS references, etc.
# such IR is hard to clean-up, so we probably will need to have the GC lowering pass emit
# lower-level intrinsics which then can be lowered to architecture-specific code.
function lower_gc_frame!(fun::LLVM.Function)
    job = current_job::CompilerJob
    mod = LLVM.parent(fun)
    ctx = context(fun)
    changed = false

    # plain alloc
    if haskey(functions(mod), "julia.gc_alloc_obj")
        alloc_obj = functions(mod)["julia.gc_alloc_obj"]
        alloc_obj_ft = eltype(llvmtype(alloc_obj))
        T_prjlvalue = return_type(alloc_obj_ft)
        T_pjlvalue = convert(LLVMType, Any; ctx, allow_boxed=true)

        for use in uses(alloc_obj)
            call = user(use)::LLVM.CallInst

            # decode the call
            ops = arguments(call)
            sz = ops[2]

            # replace with PTX alloc_obj
            Builder(ctx) do builder
                # NOTE: this happens late during the pipeline, where we may have to
                #       pass a kernel state arguments to the runtime function.
                state = kernel_state_type(job)

                position!(builder, call)
                ptr = call!(builder, Runtime.get(:gc_pool_alloc), [sz]; state)
                replace_uses!(call, ptr)
            end

            unsafe_delete!(LLVM.parent(call), call)

            changed = true
        end

        @compiler_assert isempty(uses(alloc_obj)) job
    end

    # we don't care about write barriers
    if haskey(functions(mod), "julia.write_barrier")
        barrier = functions(mod)["julia.write_barrier"]

        for use in uses(barrier)
            call = user(use)::LLVM.CallInst
            unsafe_delete!(LLVM.parent(call), call)
            changed = true
        end

        @compiler_assert isempty(uses(barrier)) job
    end

    return changed
end

# lower the `julia.ptls_states` intrinsic by removing it, since it is GPU incompatible.
#
# this assumes and checks that the TLS is unused, which should be the case for most GPU code
# after lowering the GC intrinsics to TLS-less code and having run DCE.
#
# TODO: maybe don't have Julia emit actual uses of the TLS, but use intrinsics instead,
#       making it easier to remove or reimplement that functionality here.
function lower_ptls!(mod::LLVM.Module)
    job = current_job::CompilerJob
    changed = false

    intrinsic = VERSION >= v"1.7.0-DEV.1205" ? "julia.get_pgcstack" : "julia.ptls_states"

    if haskey(functions(mod), intrinsic)
        ptls_getter = functions(mod)[intrinsic]

        for use in uses(ptls_getter)
            val = user(use)
            if isempty(uses(val))
                unsafe_delete!(LLVM.parent(val), val)
                changed = true
            else
                # the validator will detect this
            end
        end
     end

    return changed
end


text:
using Test
using HMMBase
using Distributions
using Random

Random.seed!(2019)

hmms = [
    HMM([0.9 0.1; 0.1 0.9], [Normal(10,1), Gamma(1,1)]),
    HMM([0.9 0.1; 0.1 0.9], [Categorical([0.1, 0.2, 0.7]), Categorical([0.5, 0.5])]),
    HMM([0.9 0.1; 0.1 0.9], [MvNormal([0.0,0.0], [1.0,1.0]), MvNormal([10.0,10.0], [1.0,1.0])])
]

@testset "Integration $(typeof(hmm))" for hmm in hmms
    # HMM API
    @test hmm !== copy(hmm)

    z, y = rand(hmm, 1000, seq = true)
    @test size(z, 1) == size(y, 1)
    @test size(y, 2) == size(hmm, 2)

    yp = rand(hmm, z)
    @test size(z, 1) == size(y, 1)
    @test size(y, 2) == size(hmm, 2)

    L = likelihoods(hmm, y)
    LL = likelihoods(hmm, y, logl = true)
    @test size(L) == size(LL)

    # Forward/Backward
    α1, logtot1 = forward(hmm, y)
    α2, logtot2 = forward(hmm, y, logl = true)

    @test logtot1 ≈ logtot2
    @test α1 ≈ α2

    β1, logtot3 = backward(hmm, y)
    β2, logtot4 = backward(hmm, y, logl = true)

    @test logtot3 ≈ logtot4
    @test β1 ≈ β2

    logtot5 = loglikelihood(hmm, y)

    @test size(α1) == size(α2) == size(β1) == size(β2)
    @test logtot1 ≈ logtot2 ≈ logtot3 ≈ logtot4 ≈ logtot5

    γ1 = posteriors(hmm, y)
    γ2 = posteriors(hmm, y, logl = true)

    @test γ1 ≈ γ2

    # Viterbi
    zv1 = viterbi(hmm, y)
    zv2 = viterbi(hmm, y; logl = true)
    @test size(zv1) == size(zv2) == size(z)

    # MLE
    hmm2, _ = fit_mle(hmm, y, maxiter = 1)
    @test size(hmm2) == size(hmm)
    @test typeof(hmm2) == typeof(hmm)

    hmm2, _ = fit_mle(hmm, y, init = nothing)
    @test size(hmm2) == size(hmm)
    @test typeof(hmm2) == typeof(hmm)

    hmm2, _ = fit_mle(hmm, y, init = :kmeans, robust = true)
    @test size(hmm2) == size(hmm)
    @test typeof(hmm2) == typeof(hmm)
end

# Test edge cases (no observations, one observation)
@testset "Integration T=$T" for T in [0, 1]
    hmm = HMM([0.9 0.1; 0.1 0.9], [Normal(10,1), Normal(1,1)])

    z, y = rand(hmm, 0, seq = true)
    @test size(z, 1) == size(y, 1)
    @test size(y, 2) == size(hmm, 2)

    yp = rand(hmm, z)
    @test size(z, 1) == size(y, 1)
    @test size(y, 2) == size(hmm, 2)

    L = likelihoods(hmm, y)
    LL = likelihoods(hmm, y, logl = true)
    @test size(L) == size(LL)

    # Forward/Backward
    α1, logtot1 = forward(hmm, y)
    α2, logtot2 = forward(hmm, y, logl = true)

    @test logtot1 ≈ logtot2
    @test α1 ≈ α2

    β1, logtot3 = backward(hmm, y)
    β2, logtot4 = backward(hmm, y, logl = true)

    @test logtot3 ≈ logtot4
    @test β1 ≈ β2

    logtot5 = loglikelihood(hmm, y)

    @test size(α1) == size(α2) == size(β1) == size(β2)
    @test logtot1 ≈ logtot2 ≈ logtot3 ≈ logtot4 ≈ logtot5

    γ1 = posteriors(hmm, y)
    γ2 = posteriors(hmm, y, logl = true)

    @test γ1 ≈ γ2

    # Viterbi
    zv1 = viterbi(hmm, y)
    zv2 = viterbi(hmm, y; logl = true)
    @test size(zv1) == size(zv2) == size(z)

    # Utilities
    @test_nowarn gettransmat(z)
    @test_nowarn gettransmat(z, relabel = true)
end


text:
@testset "2.add-two-numbers.jl" begin
    @test merge_k_lists([
        ListNode{Int}([1, 4, 5]), ListNode{Int}([1, 3, 4]), ListNode{Int}([2, 6])
    ]) == ListNode{Int}([1, 1, 2, 3, 4, 4, 5, 6])
    @test merge_k_lists(ListNode{Int}[]) === nothing
    @test merge_k_lists([nothing]) === nothing
end


text:
using Random
using DelimitedFiles

function benchmarkDFTFunction(f::Function, maxTime, repititions::Int, resultFile)
    n = 2
    times = Array{Float64, 1}(undef, 0)
    sizes = Array{Int, 1}(undef, 0)

    Random.seed!(30011995) #my birthday...
    duration = 0.0
    while duration <= maxTime
        n *= 2
        inputArr = rand(n)

        startTime = time_ns()
        for i = 1:repititions
            f(inputArr)
        end
        endTime = time_ns()

        duration = (endTime - startTime) / repititions * 1.0e-9
        push!(sizes, n)
        push!(times, duration)
    end

    open(resultFile, "w") do io
        writedlm(io, [sizes, times], ',')
    end
end

max_stars_repo_name: JulienKluge/ComputationalPhotonics

max_issues_repo_name: JulienKluge/ComputationalPhotonics

max_forks_repo_name: JulienKluge/ComputationalPhotonics


text:
module MathToolsMod
using LeafPhotosynthesisMod:fluxes
using Leaf:leaf_params
#-----------------------------------------------------------------------
# !DESCRIPTION:
# Math tools
#
# !USES:
#use shr_kind_mod, only : r8 => shr_kind_r8
#use clm_varctl, only : iulog
#use abortutils, only : endrun
#use CanopyFluxesMultilayerType, only : mlcanopy_type
#!
#! !PUBLIC TYPES:
#implicit none
#!
#! !PUBLIC MEMBER FUNCTIONS:
export quadratic,beta_function,log_gamma_function, hybrid, zbrent
# hybrid             ! Solve for the root of a function using secant and Brent's methods
# zbrent             ! Use Brent's method to find the root of a function
# quadratic          ! Solve a quadratic equation for its two roots
# tridiag            ! Solve a tridiagonal system of equations
# beta_function      ! Evaluate the beta function at p and q: B(p,q)
# log_gamma_function ! Evaluate the log natural of the gamma function at x: ln(G(x))

#-----------------------------------------------------------------------
function quadratic(a, b, c)
  discr = b^2 - 4*a*c
  discr >= 0 ?   ( (-b + sqrt(discr))/(2a), (-b - sqrt(discr))/(2a) ) : 0.0 #error("Only complex roots")
end # function

#-----------------------------------------------------------------------
function beta_function(p, q)
#
# !DESCRIPTION:
# Return the value of the beta function evaluated at p and q: B(p,q)
#
# !USES:
#
exp(log_gamma_function(p) + log_gamma_function(q) - log_gamma_function(p+q))

end #function beta_function

#-----------------------------------------------------------------------
function log_gamma_function(x)
#
# !DESCRIPTION:
# Return the value of the log natural of the gamma function evaluated at x: ln(G(x))
#
  coef = [76.18009172947146, -86.50532032941677,24.01409824083091, -1.231739572450155, 0.1208650973866179e-02, -0.5395239384953e-05];
  stp = 2.5066282746310005;
  y = x
  tmp = x + 5.5
  tmp = (x + 0.5) * log(tmp) - tmp
  ser = 1.000000000190015;
  for co in coef
    y = y + 1.0
    ser = ser + co / y
  end
  return tmp + log(stp*ser/x)
end # function log_gamma_function

function hybrid(flux::fluxes,leaf::leaf_params, func::Function, xa::Number, xb::Number, tol::Number)
    #
    # !DESCRIPTION:
    # Solve for the root of a function, given initial estimates xa and xb.
    # The root is updated until its accuracy is tol.
    #
    # !USES:
    #
    # !ARGUMENTS:
    #procedure (xfunc) :: func         ! Function to solve
    #real(r8), intent(in) :: xa, xb    ! Initial estimates of root
    #real(r8), intent(in) :: tol       ! Error tolerance
    #type(mlcanopy_type), intent(inout) :: mlcanopy_inst
    itmax = 40
    x0 = x1 = f0 = f1 = 0.0
    x0 = xa
    f0 = func(x0,flux, leaf)
    if (f0 == 0.0) return  x0 end

    x1 = xb
    f1 = func(x1,flux, leaf)
    if (f1 == 0.0) return x1 end

    if (f1 < f0)
       minx = x1
       minf = f1
    else
       minx = x0
       minf = f0
    end

    # First use the secant method, and then use the brent method as a backup
    for iter = 0:itmax
       iter = iter + 1
       dx = -f1 * (x1 - x0) / (f1 - f0)
       x = x1 + dx
       if (abs(dx) < tol)
          x0 = x
          break
       end
       x0 = x1
       f0 = f1
       x1 = x
       f1 = func(x1,flux, leaf)
       if (f1 < minf)
          minx = x1
          minf = f1
       end

       # If a root zone is found, use the brent method for a robust backup strategy
       if (f1 * f0 < 0.0)
          x = zbrent(flux,  leaf,func,x0, x1, xtol=tol)
          x0 = x
          break
       end

       # In case of failing to converge within itmax iterations stop at the minimum function
       if (iter > itmax)
          f1 = func(minx,flux, leaf)
          x0 = minx
          break
       end

    end

    return x0
end #function hybrid
function zbrent(flux::fluxes,leaf::leaf_params,f::Function, x0::Number, x1::Number, args::Tuple=();
               xtol::AbstractFloat=1e-7, ytol=2eps(Float64),
               maxiter::Integer=50)
    EPS = eps(Float64)
    y0 = f(x0,flux,leaf)
    y1 = f(x1,flux,leaf)
    if abs(y0) < abs(y1)
        # Swap lower and upper bounds.
        x0, x1 = x1, x0
        y0, y1 = y1, y0
    end
    x2 = x0
    y2 = y0
    x3 = x2
    bisection = true
    for _ in 1:maxiter
        # x-tolerance.
        if abs(x1-x0) < xtol
            return x1
        end

        # Use inverse quadratic interpolation if f(x0)!=f(x1)!=f(x2)
        # and linear interpolation (secant method) otherwise.
        if abs(y0-y2) > ytol && abs(y1-y2) > ytol
            x = x0*y1*y2/((y0-y1)*(y0-y2)) +
                x1*y0*y2/((y1-y0)*(y1-y2)) +
                x2*y0*y1/((y2-y0)*(y2-y1))
        else
            x = x1 - y1 * (x1-x0)/(y1-y0)
        end

        # Use bisection method if satisfies the conditions.
        delta = abs(2EPS*abs(x1))
        min1 = abs(x-x1)
        min2 = abs(x1-x2)
        min3 = abs(x2-x3)
        if (x < (3x0+x1)/4 && x > x1) ||
           (bisection && min1 >= min2/2) ||
           (!bisection && min1 >= min3/2) ||
           (bisection && min2 < delta) ||
           (!bisection && min3 < delta)
            x = (x0+x1)/2
            bisection = true
        else
            bisection = false
        end

        y = f(x,flux,leaf)
        # y-tolerance.
        if abs(y) < ytol
            return x
        end
        x3 = x2
        x2 = x1
        if sign(y0) != sign(y)
            x1 = x
            y1 = y
        else
            x0 = x
            y0 = y
        end
        if abs(y0) < abs(y1)
            # Swap lower and upper bounds.
            x0, x1 = x1, x0
            y0, y1 = y1, y0
        end
    end
    error("Max iteration exceeded")
end
end # module MathToolsMod

max_stars_repo_name: ruillercoaster/LSM-SPAM

max_issues_repo_name: ruillercoaster/LSM-SPAM

max_forks_repo_name: ruillercoaster/LSM-SPAM


text:
module TestBackup

using LibCloud.Backup
using Base.Test
using PyCall

function dummy()
    cp = BackupDriver(BackupProvider.DUMMY, "key", "secret")
    println("Regions:\n", list_regions(cp))
end

end # module TestBackup


text:
using StatsBase
import StatsBase.entropy

function entropy{T<:Integer}(x::Array{T,1})
    N = sum(x)
    K = length(x)
    P = x/N #convert to probabilitiy
    return -dot(P,log2(P + (P.==0)))
end

function entropy(NN::Hist2d)
    Nc = NN.wordcnt
    N = sum(Nc[:])
    H1 = log2(Nc + (Nc.==0))
    H2 = log2(sum(Nc,1))
    H3 = log2(sum(Nc,2))
    H4 = log2(N)
    H = Nc.*(broadcast(-,broadcast(-, H1, H2),H3) + H4)
    H[isnan(H)] = 0
    H = 1/N*sum(sum(H))
    return H
end

#functions operating on tress
function entropy(X::Dict{Int64, Dict{Int64, Dict{Int64,Int64}}})
    E = 0
    N = 0
    for (k,v) in X
        ee,n = entropy(v)
        E += ee*n
        N += n
    end
    return E/N,N
end

function entropy(X::Dict{Int64, Dict{Int64,Int64}})
    E = 0
    N = 0
    for (k,v) in X
        ee,n = entropy(v)
        E += ee*n
        N += n
    end
    return E/N,N
end

function entropy(X::Dict{Int64,Int64})
    N = sum(values(X))
    P = float(map(v -> v/N, values(X)))
    return (-sum(P.*log2(P + float(P.==0))),N)
end


text:
"""
    preproc(A, B) is the preprocessing step for computing the
    Generalized Singular Value Decomposition (GSVD), it computes:

```math
    A = U * RA * Q' ,    B = V * RB * Q'
```
where
- U is m-by-m orthogonal matrix,
- V is p-by-p orthogonal matrix,
- Q is n-by-n orthogonal matrix,
- RA and RB have the following structures:

(1) if m-k-l >= 0:

                   n-k-l   k    l                   n-k-l  k    l
       RA =     k ( 0     A12  A13 )     RB =   l ( 0      0   B13 )
                l ( 0      0   A23 )          p-l ( 0      0    0  )
            m-k-l ( 0      0    0  ),

(2) if m-k-l < 0:

                   n-k-l   k    l                   n-k-l  k    l
       RA =     k ( 0     A12  A13 )     RB =   l ( 0      0   B13 )
              m-k ( 0      0   A23 ),          p-l ( 0      0    0  )

where
- k-by-k matrix A12 and l-by-l matrix B13 are nonsingular upper triangular,
- A23 is l-by-l upper triangular if m-k-l >= 0,
  otherwise A23 is (m-k)-by-l upper trapezoidal,
- k+l is the effective numerical rank of the (m+p)-by-n matrix [A;B],
  l is the effective numerical rank of B.

# arguments:
A: m-by-n Float 32 or Float64 matrix,
B: p-by-n Float 32 or Float64 matrix.

# returns:
- U: m-by-m matrix,
- V: p-by-p matrix,
- Q: n-by-n matrix,
- k: interger, equals rank([A;B]) - rank(B),
- l: interger, equals rank(B),
- RA: m-by-n matrix, overwrites A in preproc!(),
- RB: p-by-n matrix, overwrites B in preproc!().
"""
function preproc!(A::StridedMatrix{T}, B::StridedMatrix{T}) where T
    m, n = size(A)
    p = size(B)[1]
    tola = max(m,n)*opnorm(A, 1)*eps(T)
    tolb = max(p,n)*opnorm(B, 1)*eps(T)

    # Step 1: QR decomposition with
    # col pivoting of B
    #               l   n-l
    # B * P = V * (B11  B12) l
    #             ( 0    0 ) p-l
    F1 = qr!(B, Val(true))
    if p > n
        V = F1.Q*Matrix{T}(I, p, p)
    else
        V = Matrix{T}(F1.Q)
    end

    # Step 2: Update A
    # A = A*P
    A = A*F1.P

    # Step 3: Set Q
    Q = F1.P

    # Determine the numerical rank of B
    l = 0
    for i = 1:min(p,n)
        if abs(B[i,i]) > tolb
            l += 1
        end
    end

    # Clean up B to make it upper triangular
    @views triu!(B[1:l,1:l])
    if p > l
        @views B[l+1:p,1:n] .= T(0.0)
    end

    # Step 2: RQ decomposition of (B11  B12)
    # (B11  B12) = (0  B22) * Z
    if p >= l && n != l
        B[1:l,:], tau1 = @views LAPACK.gerqf!(B[1:l,:])
        @views LAPACK.ormrq!('R', 'T', B[1:l,:], tau1, A)
        @views LAPACK.ormrq!('R', 'T', B[1:l,:], tau1, Q)
        # Clean up B
        @views B[1:l,1:n-l] .= T(0.0)
        @views triu!(B[1:l,n-l+1:n])
    end

    # Step 3: QR decomposition with
    # col pivoting of A11, if A11 is not empty
    #            n-l    l
    #       A = (A11  A12) m
    #
    #       A11 * P = U * (T11  T12)
    #                     ( 0    0 )
    F2 = @views qr!(A[:,1:n-l], Val(true))

    # Determine the numerical rank of A11
    k = 0
    for i = 1:min(m, n-l)
        if abs(A[i,i]) > tola
            k += 1
        end
    end

    if m > n-l
        U = F2.Q*Matrix{T}(I, m, m)
    else
        U = Matrix{T}(F2.Q)
    end

    A[1:m, n-l+1:n] = @views U' * A[1:m, n-l+1:n]
    Q[:,1:n-l] = @views Q[:,1:n-l]*F2.P

    # Clean up A to make it upper triangular
    @views triu!(A[1:k,1:k])
    if m > k
        @views A[k+1:m,1:n-l] .= T(0.0)
    end

    # Step 4: RQ decomposition of (T11  T12)
    # (T11  T12) = (0  T12) * Z
    if n-l > k
        A[1:k,1:n-l], tau2 = @views LAPACK.gerqf!(A[1:k,1:n-l])
        @views LAPACK.ormrq!('R', 'T', A[1:k,1:n-l], tau2, Q[1:n,1:n-l])
        # Clean up A
        @views A[1:k,1:n-l-k] .= T(0.0)
        @views triu!(A[1:k,n-l-k+1:n-l])
    end

    # Step 5: QR decomposition of A[k+1:m,n-l+1:n]
    if m > k
        # Weird, you have to add @views for submatrix input, otherwise,
        # it will not be overwritten.
        # A[k+1:m,n-l+1:n], tau3 = @views LAPACK.geqrf!(A[k+1:m,n-l+1:n])
        # @views LAPACK.ormqr!('R', 'N', A[k+1:m,n-l+1:n], tau3, U[:,k+1:m])
        @views F3 = qr!(A[k+1:m,n-l+1:n])
        if m - k > l
            U[:,k+1:m] = @views U[:,k+1:m]*(F3.Q*Matrix{T}(I,m-k,m-k))
        else
            U[:,k+1:m] = @views U[:,k+1:m]*Matrix{T}(F3.Q)
        end
        @views triu!(A[k+1:m,n-l+1:n])
    end

    return U, V, Q, k, l, A, B
end

# preproc() calls preproc!() before making a copy of input matrices A and B,
# thus A, B will not be overwritten
function preproc(A::StridedMatrix{T}, B::StridedMatrix{T}) where T
    return preproc!(copy(A),copy(B))
end


text:
function sum_diff(x)
n = length(x); d = 1/(n-1)
s = zero(eltype(x))
s = s + (x[2] - x[1]) / d
for i = 2:length(x)-1
s = s + (x[i+1] - x[i+1]) / (2*d)
end
s = s + (x[n] - x[n-1])/d
end
function sum_diff_fast(x)
n=length(x); d = 1/(n-1)
s = zero(eltype(x))
@fastmath s = s + (x[2] - x[1]) / d
@fastmath for i = 2:n-1
s = s + (x[i+1] - x[i+1]) / (2*d)
end
@fastmath s = s + (x[n] - x[n-1])/d
end

max_stars_repo_path: Module 2/Chapter05/fastmath_mac.jl

max_issues_repo_path: Module 2/Chapter05/fastmath_mac.jl

max_forks_repo_path: Module 2/Chapter05/fastmath_mac.jl