text:
export simple_grid

function simple_grid(F::AbstractFiber, grid, generator::StartValueGenerator, avoid_check=_always_check, B=empty(grid); options=MembershipTestOptions())
    B = empty(grid)
    Fs = [F]
    for _=2:Threads.nthreads()
        push!(Fs, deepcopy(F))
    end

    partition = partition_work(length(grid))
    Threads.@threads for i=1:Threads.nthreads()
        simple_grid_kernel!(B, Fs[i], partition[i], grid, generator, avoid_check, options)
    end
    B
end

function simple_grid_kernel!(B, F::AbstractFiber, range, grid, generator::StartValueGenerator, avoid_check, options)
    for k in range
        if !B[k]
            w = grid[k]
            avoid, val = avoid_check(B, k)
            if !avoid
                B[k] = membershiptest(F, w, generator, options).successfull
            else
                B[k] = val
            end
        end
    end
end

_always_check(B, x) = (false, false)


text:
import Ghost: trace, V, Call, play!

inc_mul(a::Real, b::Real) = a * (b + 1.0)
inc_mul(A::AbstractArray, B::AbstractArray) = inc_mul.(A, B)
inc_mul2(A::AbstractArray, B::AbstractArray) = A .* (B .+ 1)

non_primitive(x) = 2x + 1
non_primitive_caller(x) = sin(non_primitive(x))

make_tuple(a, b) = (a, b)

function tuple_unpack(x)
    a, b = make_tuple(1, 2)
    x + a
end

function make_adder(x)
    return y -> x + y
end

function add(x, y)
    adder = make_adder(x)
    adder(y)
end

function constant_return_value(x)
    2x
    nothing
end

function no_input()
    print()
end

@testset "trace" begin
    # calls
    val, tape = trace(inc_mul, 2.0, 3.0)
    @test val == inc_mul(2.0, 3.0)
    @test length(tape) == 5
    @test tape[V(5)].args[1].id == 2

    # bcast
    A = rand(3)
    B = rand(3)
    val, tape = trace(inc_mul, A, B)
    @test val == inc_mul(A, B)
    # broadcasting may be lowered to different forms,
    # so making no assumptions regarding the tape

    val, tape = trace(inc_mul2, A, B)
    @test val == inc_mul2(A, B)

    # primitives
    x = 3.0
    val1, tape1 = trace(non_primitive_caller, x)
    val2, tape2 = trace(non_primitive_caller, x; primitives=Set([non_primitive, sin]))

    @test val1 == val2
    @test any(op isa Call && op.fn == (*) for op in tape1)
    @test tape2[V(3)].fn == non_primitive
    @test tape2[V(4)].fn == sin

    # tuple unpacking
    val, tape = trace(tuple_unpack, 4.0)
    @test val == tuple_unpack(4.0)

    # closures
    val, tape = trace(add, 2.0, 4.0)
    @test val == add(2.0, 4.0)
    @test play!(tape, add, 2.0, 5.0) == add(2.0, 5.0)
    @test play!(tape, add, 3.0, 5.0) == add(3.0, 5.0)

    # iterators
    xs = rand(3)
    f = xs -> [x + 1 for x in xs]
    val, tape = trace(f, xs)
    @test val == f(xs)
    xs2 = rand(5)
    @test play!(tape, nothing, xs2) == f(xs2)

    # constant return value
    _, tape = trace(constant_return_value, 1.0)
    @test play!(tape, 2.0) === nothing

    # no input
    _, tape = trace(no_input)
    @test tape[V(2)].fn == print
end

function return_in_branch(x)
    if x > 0
        return 2x
    end
    return 3x
end

@testset "trace: branches" begin
    # See https://github.com/dfdx/Ghost.jl/issues/8
    _, tape = trace(return_in_branch, 1.0)
    @test tape[V(length(tape))].args[1] == 2

    _, tape = trace(return_in_branch, -1.0)
    @test tape[V(length(tape))].args[1] == 3
end

function loop1(a, n)
    a = 2a
    for i in 1:n
        a = a * n
    end
    return a
end

function loop2(a, b)
    while b > 0
        a = a * b
        b = b - 1
    end
    return a
end

function loop3(a, b)
    while b > 1
        b = b - 1
        a = b
        while a < 100
            a = a * b + 1
        end
    end
    return a
end

function loop4(x, n, m)
    for i in 1:n
        for j in 1:m
            x = 2x
        end
    end
    return x
end

function loop5(a, n)
    for i=1:3
        a = loop1(a, n)
    end
    return a
end

function loop6(n)
    a = 0
    i = 1
    while true
        a += i
        if a > n
            break
        end
    end
    return a
end

function myabs(x)
    if x >= 0
        return x
    else
        return -x
    end
end

function myabs2(x)
    if x <= 0
        return -x
    else
        return 1*x
    end
end

@testset "trace: loops" begin
    should_trace_loops!(false)

    _, tape = trace(loop1, 1.0, 3)
    @test findfirst(op -> op isa Loop, tape.ops) === nothing
    # same number of iteration
    @test play!(tape, loop1, 2.0, 3) == loop1(2.0, 3)
    @test compile(tape)(loop1, 2.0, 3) == loop1(2.0, 3)
    # different number of iteration - with loop tracing off, should be incorrect
    @test play!(tape, loop1, 2.0, 4) != loop1(2.0, 4)
    @test compile(tape)(loop1, 2.0, 4) != loop1(2.0, 4)

    should_trace_loops!(true)

    _, tape = trace(loop1, 1.0, 3)
    @test play!(tape, loop1, 2.0, 4) == loop1(2.0, 4)
    @test compile(tape)(loop1, 2.0, 4) == loop1(2.0, 4)
    @test findfirst(op -> op isa Loop, tape.ops) !== nothing

    _, tape = trace(loop2, 1.0, 3)
    @test play!(tape, loop2, 2.0, 4) == loop2(2.0, 4)
    @test compile(tape)(loop2, 2.0, 4) == loop2(2.0, 4)
    @test findfirst(op -> op isa Loop, tape.ops) !== nothing

    _, tape = trace(loop3, 1.0, 3)
    @test play!(tape, loop3, 2.0, 4) == loop3(2.0, 4)
    @test compile(tape)(loop3, 2.0, 4) == loop3(2.0, 4)
    @test findfirst(op -> op isa Loop, tape.ops) !== nothing

    _, tape = trace(loop4, 1.0, 2, 3)
    @test play!(tape, loop4, 2.0, 3, 4) == loop4(2.0, 3, 4)
    @test compile(tape)(loop4, 2.0, 3, 4) == loop4(2.0, 3, 4)
    loop_idx = findfirst(op -> op isa Loop, tape.ops)
    @test loop_idx !== nothing
    subtape = tape[V(loop_idx)].subtape
    @test findfirst(op -> op isa Loop, subtape.ops) !== nothing

    _, tape = trace(loop5, 1.0, 3)
    @test play!(tape, loop5, 2.0, 4) == loop5(2.0, 4)
    @test compile(tape)(loop5, 2.0, 4) == loop5(2.0, 4)
    loop_idx = findfirst(op -> op isa Loop, tape.ops)
    @test loop_idx !== nothing
    subtape = tape[V(loop_idx)].subtape
    @test findfirst(op -> op isa Loop, subtape.ops) !== nothing

    should_trace_loops!()

    # Test with fixed boolean condition
    # Currently broken
    #=
    _, tape = trace(loop6, 3)
    @test play!(tape, loop6, 3) == loop6(3)
    @test compile(tape, loop6, 3) == loop6(3)
    =#
end

@testset "trace: branches" begin
    should_assert_branches!(true)

    _, tape = trace(myabs, 3)
    @test play!(tape, myabs, 3) == 3
    @test_throws AssertionError play!(tape, myabs, -3) # Throws assert exception due to bad branch

    _, tape = trace(myabs2, 3)
    @test play!(tape, myabs2, 3) == 3
    @test_throws AssertionError play!(tape, myabs2, -3) # Test with branching flipped

    should_assert_branches!(false)
end

ext: jl


text:
module MajoranaReps

export bx,by,bz,c,Sx,Sy,Sz,Ket, Down, Up
export OpSqrt, OpToStr, OpScaleType, OpProp, OpScale, OpAdd
export MajElem, MajScalar, MajoranaOp, OpReduce, ParityOp
export OpNonComutes, OpComutes, OpInnerProd, Majorana, MajoranaSp

if !@isdefined(TabLevel)
    TabLevel = ""
end
println(TabLevel*"Open MajoranaReps.jl")
TabLevel=TabLevel*"    "

include("MajoranaBase.jl")
include("OpInnerProd.jl")

####We will also overlaod the show
###We will fix the strings later
import Base.show
function Base.show(io::IO,object::MajoranaOp)
    print(io,OpToStr(object))
end
function Base.show(io::IO,object::MajoranaSp)
    print(io,OpToStr(object))
end

NumStr = Union{Number,String}

###The spin represenations
Majorana(A,i::NumStr) = MajoranaOp(1,[MajElem(A,i)])
bx(i::NumStr) = Majorana("bx",i); bx() = bx("")
by(i::NumStr) = Majorana("by",i); by() = by("")
bz(i::NumStr) = Majorana("bz",i); bz() = bz("")
c(i::NumStr) = Majorana("c",i); c() = c("")
Sx(i::NumStr) = im*bx(i)*c(i)
Sy(i::NumStr) = im*by(i)*c(i)
Sz(i::NumStr) = im*bz(i)*c(i)
Ket() = MajoranaOp(1,fill(MajElem("",0),0),true)
Down() = Ket()

function Up(iList...)
    OpCollect=Down()
    for i in iList
        OpCollect= Sx(i)*OpCollect
    end
    return OpCollect
end
###The parity operator
function ParityOp(i)
    bx(i)*by(i)*bz(i)*c(i)
end

function ValidateOp(Op)
    if !(Op=="bx" || Op=="by" || Op=="bz" ||Op=="c")
        thow(DomainError(Op,"Not a valid operator name"))
    end
end

function OpOrderedOld(ElemA,ElemB) ###Check if operators are ordered
    ###We implement the following ordering
    ## bx_j < c_j and (bx/c)_k < (bx/c)_j if k<j
    ## by_j < bz_j and (by/bz)_k < (by/bz)_j if k<j
    ## (bx/c)_j < (by/bz)_k for all j,k
    (OpA,iA)=ElemA
    (OpB,iB)=ElemB
    ValidateOp(OpA)
    ValidateOp(OpB)
    if OpA=="bx"
        if (OpB=="by" || OpB=="bz")
            return true ###correctly ordered
        elseif OpB=="bx" ###Check indexes
            return iA <= iB
        elseif OpB=="c" ###Check indexes
            return  iA <= iB
        end
    elseif OpA=="c"
        if (OpB=="by" || OpB=="bz")
            return true ###correctly ordered
        elseif OpB=="bx" ###Check indexes
            return iA < iB ##if indexes are the same the ordering is wrong
        elseif OpB=="c" ###Check indexes
            return  iA <= iB
        end
    elseif OpA=="by"
        if (OpB=="bx" || OpB=="c")
            return false ###incorrectly ordered
        elseif OpB=="by" ###Check indexes
            return iA <= iB
        elseif OpB=="bz" ###Check indexes
            return  iA <= iB
        end
    elseif OpA=="bz"
        if (OpB=="bx" || OpB=="c")
            return false ###incorrectly ordered
        elseif OpB=="by" ###Check indexes
            return iA < iB
        elseif OpB=="bz" ###Check indexes
            return  iA <= iB
        end
    end
    throw(DomainError("Case was never caught"))
end

##Individual terms cannor be reordered
function OpReOrder(a::MajoranaOp) return OpCopy(a) end

function OpReOrder(MajSp_in::MajoranaSp)
    verbose=false
    ### Loop though the function and apply the bubble sort algorithm to make sure the states are in order
    MajSp = OpCopy(MajSp_in)
    num_factors=length(MajSp.SpList)
    indx = 1
    while indx < num_factors
        if indx<=0 ; indx += 1 ; continue  end ##Sometimes we go back to long
        ElemA=MajSp.SpList[indx]
        ElemB=MajSp.SpList[indx+1]
        verbose && println("ElemA: ",OpToStr(ElemA))
        verbose && println("ElemB: ",OpToStr(ElemB))
        if OpProp(ElemA ,ElemB ) ###If poportional, simply merge the two
            verbose && println("Add the factors")
            MajSp.SpList[indx]=OpAdd(ElemA,ElemB)
            ##Then remove the B elems
            deleteat!(MajSp.SpList,indx+1)
            ##And step down the total count
            num_factors-=1
            if MajSp.SpList[indx].scale==0 ##IF the element is zero
                if num_factors!=1 ###And this is not the last element
                    verbose && println("Remoce sinze zero")
                    ###Remove this element also
                    deleteat!(MajSp.SpList,indx)
                    num_factors-=1
                    indx-=1 ##Take a step back
                else
                    verbose && println("Zero but keep")
                end
            end

        elseif ElemA.OpList < ElemB.OpList
            ###Correctly ordered Go to next pair
            indx += 1
        elseif ElemA.OpList > ElemB.OpList
            ###Reversly ordered ordered, fix and take one step back
            MajSp.SpList[indx]=OpCopy(ElemB)
            MajSp.SpList[indx+1]=OpCopy(ElemA)
            indx-=1
        else
            throw(DomainError("Should not happen"))
        end
    end

    return MajSp
end

function OpNonComutes(A::MajoranaOp,B::MajoranaOp;strict=true)
    commutes = OpComutes(A::MajoranaOp,B::MajoranaOp;strict=false)
    if strict && commutes
        throw(DomainError("A="*OpToStr(A)*" and B="*OpToStr(B)*" commutes!"))
    else return !commutes end
end

function OpComutes(A::MajoranaOp,B::MajoranaOp;strict=true)
    AB=A*B
    BA=B*A
    commutes= OpEqual(AB,BA,strict=false)
    if strict && !commutes
        throw(DomainError("A="*OpToStr(A)*" and B="*OpToStr(B)*" does not commute!\n"*
                          "AB="*OpToStr(AB)*"\nand\nBA="*OpToStr(BA)))
    else return commutes end
end

function BinToOcc(BinRepIn,NumOrbitals)
    BinRep=BinRepIn
    OccRep=BitArray(fill(false,NumOrbitals))
    for indx in 1:NumOrbitals
        OccRep[indx]=mod(BinRep,2)
        BinRep = BinRep >> 1
    end
    BinRep != 0 && error("BinRep = $BinRepIn encodes more than $NumOrbitals orbitals")
    return OccRep
end

TabLevel=TabLevel[1:end-4]
println(TabLevel*"Close MajoranaReps.jl")

end #Module


text:
using Gtk
using Sockets
params = LakeShoreGaussMeterParams(ip=ip"192.168.2.2")
gauss = LakeShoreGaussMeter(deviceID="lakeshoref", params=params, dependencies=Dict{String, Union{Device, Missing}}())

if ask_dialog("You are trying to start the LakeShore gaussmeter hardware test. Please ensure that the device is turned on and working.","Cancel","Start test")
    @error "Not yet implemented!"

    @test ask_dialog("Did everything seem normal?")
else
    @test_broken 0
end


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

JLLWrappers.@generate_main_file_header("XZ")
JLLWrappers.@generate_main_file("XZ", UUID("ffd25f8a-64ca-5728-b0f7-c24cf3aae800"))
end  # module XZ_jll


text:
include("../../util/src/util_module.jl")

module dev

module matrix_vector_prod
using LinearAlgebra

function matrix_vector_prod_default(A::Matrix{T}, x::Vector{T}) where {T<:Float64}
    sz_a = size(A)
    sz_b = size(x)

    if length(sz_a) != 2
        error("A is not a matrix")
    end
    if length(sz_b) != 1
        error("B is not a vector")
    end
    if sz_a[2] != sz_b[1]
        error("A and B are not compatible (size mismatch) ")
    end
    c = A * x
    return c
end

function matrix_vector_prod_inner(A, x)
    # A in R^m-by-p
    # x in R^p-by-1
    # c = A * x in R^m-by-1
    sz_a = size(A)
    sz_b = size(x)

    if length(sz_a) != 2
        error("A is not a matrix")
    end
    if length(sz_b) != 1
        error("B is not a vector")
    end
    if sz_a[2] != sz_b[1]
        error("A and B are not compatible (size mismatch) ")
    end

    c = zeros(eltype(A), sz_a[1])
    @inbounds @simd for ii = 1:sz_a[1]
        # c[ii] = A[ii, :]' * x
        c[ii] =  dot( @view(A[ii, :]), x)
    end
    return c
end

function matrix_vector_prod_inner_loop(A, x)
    # A in R^m-by-p
    # x in R^p-by-1
    # c = A * x in R^m-by-1
    sz_a = size(A)
    sz_b = size(x)

    if length(sz_a) != 2
        error("A is not a matrix")
    end
    if length(sz_b) != 1
        error("B is not a vector")
    end
    if sz_a[2] != sz_b[1]
        error("A and B are not compatible (size mismatch) ")
    end

    c = zeros(eltype(A), sz_a[1])
    @inbounds @simd for ii = 1:sz_a[1]
        inner_sum = zero(eltype(x))
        @inbounds @simd for jj = 1:sz_a[2]
            inner_sum += A[ii, jj] * x[jj]
        end
        c[ii] = inner_sum
    end
    return c
end

function matrix_vector_prod_outer(A, x)
    # A in R^m-by-p
    # x in R^p-by-1
    # c = A * x in R^m-by-1
    sz_a = size(A)
    sz_b = size(x)

    if length(sz_a) != 2
        error("A is not a matrix")
    end
    if length(sz_b) != 1
        error("B is not a vector")
    end
    if sz_a[2] != sz_b[1]
        error("A and B are not compatible (size mismatch) ")
    end

    c = zeros(eltype(A), sz_a[1])
    @inbounds @simd for jj = 1:sz_a[2]
        c .+= @view( A[:, jj] ) .* x[jj]
    end
    return c
end

function matrix_vector_prod_outer_loop(A, x)
    # A in R^m-by-p
    # x in R^p-by-1
    # c = A * x in R^m-by-1
    sz_a = size(A)
    sz_b = size(x)

    if length(sz_a) != 2
        error("A is not a matrix")
    end
    if length(sz_b) != 1
        error("B is not a vector")
    end
    if sz_a[2] != sz_b[1]
        error("A and B are not compatible (size mismatch) ")
    end

    c = zeros(eltype(A), sz_a[1])
    @inbounds @simd for jj = 1:sz_a[2]
        @inbounds @simd for ii = 1:sz_a[1]
            c[ii] += A[ii, jj] * x[jj]
        end
    end
    return c
end

function matrix_mult_inner(A, B)
    # A in R^m-by-p
    # B in R^p-by-n
    # C = A * B in R^m-by-n
    # [ma, p] = size(A)
    # [, p] = size(A)
    sz_a = size(A)
    sz_b = size(B)

    if length(sz_a) != 2
        error("A is not a matrix")
    end
    if length(sz_b) != 2
        error("B is not a matrix")
    end
    if sz_a[2] != sz_b[1]
        error("A and B are not compatible (size mismatch) ")
    end

    C = zeros(eltype(A), sz_a[1], sz_b[2])
    for ii = 1:sz_a[1]
        for jj = 1:sz_b[2]
            C[ii, jj] = dot( A[ii, :],  B[:, jj] )
            # C[ii, jj] = A[ii, :]' *  B[:, jj]
        end
    end
    return C
end

function matrix_mult_outer(A, B)
    # A in R^m-by-p
    # B in R^p-by-n
    # C = A * B in R^m-by-n
    # [ma, p] = size(A)
    # [, p] = size(A)
    sz_a = size(A)
    sz_b = size(B)

    if length(sz_a) != 2
        error("A is not a matrix")
    end
    if length(sz_b) != 2
        error("B is not a matrix")
    end
    if sz_a[2] != sz_b[1]
        error("A and B are not compatible (size mismatch) ")
    end

    C = zeros(eltype(A), sz_a[1], sz_b[2])
    for kk = 1:sz_a[2]
        C += A[:, kk] * B[kk, :]'
    end
    return C
end

function matrix_mult_default(A, B)
    # A in R^m-by-p
    # B in R^p-by-n
    # C = A * B in R^m-by-n
    # [ma, p] = size(A)
    # [, p] = size(A)
    sz_a = size(A)
    sz_b = size(B)

    if length(sz_a) != 2
        error("A is not a matrix")
    end
    if length(sz_b) != 2
        error("B is not a matrix")
    end
    if sz_a[2] != sz_b[1]
        error("A and B are not compatible (size mismatch) ")
    end
    C = A * B
    return C
end

function matrix_mult_inner_loop(A, B)
    # A in R^m-by-p
    # B in R^p-by-n
    # C = A * B in R^m-by-n
    # [ma, p] = size(A)
    # [, p] = size(A)
    sz_a = size(A)
    sz_b = size(B)

    if length(sz_a) != 2
        error("A is not a matrix")
    end
    if length(sz_b) != 2
        error("B is not a matrix")
    end
    if sz_a[2] != sz_b[1]
        error("A and B are not compatible (size mismatch) ")
    end
    C = zeros(eltype(A), sz_a[1], sz_b[2])

    for ii = 1:sz_a[1]
        for jj = 1:sz_b[2]
            inner_sum = zero(eltype(A))
            for kk = 1:sz_a[2]
                inner_sum += A[ii, kk] *  B[kk, jj]
            end
            C[ii, jj] = inner_sum
        end
    end
    return C
end

function matrix_mult_outer_loop(A, B)
    #
    # Algo probably wrong
    #
    # A in R^m-by-p
    # B in R^p-by-n
    # C = A * B in R^m-by-n
    # [ma, p] = size(A)
    # [, p] = size(A)
    sz_a = size(A)
    sz_b = size(B)

    if length(sz_a) != 2
        error("A is not a matrix")
    end
    if length(sz_b) != 2
        error("B is not a matrix")
    end
    if sz_a[2] != sz_b[1]
        error("A and B are not compatible (size mismatch) ")
    end
    C = zeros(eltype(A), sz_a[1], sz_b[2])

    #=
    for kk = 1:sz_a[2]
        for ii = 1:sz_a[1]
            C[ii, kk] += A[ii, kk] * B[kk, ii]
        end
    end
    =#
    for ii = 1:sz_b[2]
        for kk = 1:sz_a[2]
            C[:, ii] += A[:, kk] * B[kk, ii]
        end
    end

    return C
end

end # dev

module devtest
using ..dev
using BenchmarkTools

function benchmark_matrix_vector_prod_default(A, x)
    b = @benchmark dev.matrix_vector_prod_default($A, $x)

    println("=== matrix_vector_prod_default ===")
    show(
            stdout,
            MIME("text/plain"),
            b )
    println("\n\n")
end

function benchmark_matrix_vector_prod_inner(A, x)
    b = @benchmark dev.matrix_vector_prod_inner($A, $x)

    println("=== matrix_vector_prod_inner ===")
    show(
            stdout,
            MIME("text/plain"),
            b )
    println("\n\n")
end

function benchmark_matrix_vector_prod_outer(A, x)
    b = @benchmark dev.matrix_vector_prod_outer($A, $x)

    println("=== matrix_vector_prod_outer ===")
    show(
            stdout,
            MIME("text/plain"),
            b )
    println("\n\n")
end

function benchmark_matrix_vector_prod_inner_loop(A, x)
    b = @benchmark dev.matrix_vector_prod_inner_loop($A, $x)

    println("=== matrix_vector_prod_inner_loop ===")
    show(
            stdout,
            MIME("text/plain"),
            b )
    println("\n\n")
end

function benchmark_matrix_vector_prod_outer_loop(A, x)
    b = @benchmark dev.matrix_vector_prod_outer_loop($A, $x)

    println("=== matrix_vector_prod_outer_loop ===")
    show(
            stdout,
            MIME("text/plain"),
            b )
    println("\n\n")
end

function benchmark_all()

    # size definitions
    # m = 1000
    # p = 5000

    # m = 10
    # p = 500000

    m = 100000
    p = 50

    # allocate the matrices
    A = rand(m, p)
    x = rand(p)

    benchmark_matrix_vector_prod_default(A, x)
    benchmark_matrix_vector_prod_inner(A, x)
    benchmark_matrix_vector_prod_outer(A, x)
    benchmark_matrix_vector_prod_inner_loop(A, x)
    benchmark_matrix_vector_prod_outer_loop(A, x)

end

end

module devtest_matrix
    using ..dev
    using BenchmarkTools

function benchmark_matrix_mult_inner(A, B)
    b = @benchmark dev.matrix_mult_inner($A, $B)

    println("=== matrix_mult_inner ===")
    show(
            stdout,
            MIME("text/plain"),
            b )
    println("\n\n")
   end

function benchmark_matrix_mult_outer(A, B)
    b = @benchmark dev.matrix_mult_outer($A, $B)

    println("=== matrix_mult_outer ===")
    show(
            stdout,
            MIME("text/plain"),
            b )
    println("\n\n")
end

function benchmark_matrix_mult_default(A, B)
    b = @benchmark dev.matrix_mult_default($A, $B)

    println("=== matrix_mult_default ===")
    show(
            stdout,
            MIME("text/plain"),
            b )
    println("\n\n")
end

function benchmark_matrix_mult_inner_loop(A, B)
    b = @benchmark dev.matrix_mult_inner_loop($A, $B)

    println("=== matrix_mult_inner_loop ===")
    show(
            stdout,
            MIME("text/plain"),
            b )
    println("\n\n")
end

function benchmark_matrix_mult_outer_loop(A, B)
    b = @benchmark dev.matrix_mult_outer_loop($A, $B)

    println("=== matrix_mult_outer_loop ===")
    show(
            stdout,
            MIME("text/plain"),
            b )
    println("\n\n")
end

function benchmark_all()
    # size definitions
    m = 100
    p = 500
    n = 100

    # allocate the matrices
    A = rand(m, p)
    B = rand(p, n)

    benchmark_matrix_mult_inner(A, B)
    benchmark_matrix_mult_outer(A, B)
    benchmark_matrix_mult_default(A, B)
    benchmark_matrix_mult_inner_loop(A, B)
    # benchmark_matrix_mult_outer_loop(A, B)
end

    # multiply
    #=
    Ci = dev.matrix_mult_inner(A, B)
    Co = dev.matrix_mult_outer(A, B)
    Cd = A * B;

    println( all( abs.( Ci .- Co ) .<= 1e-10 ) ) =#

end

#=
module plottest

using PyPlot
PyPlot.pygui(true)

x = LinRange(-2, 2, 30) # [1, 2, 3]
y = LinRange(-4, 1, 30) # [10, 20, 30]
z = [i.^2 + j.^2 for i in x, j in y]
z2 = [50, 51, 52]

PyPlot.figure()
# PyPlot.plot3D(x, y, z2)
PyPlot.plot_surface(x, y, z, alpha = 0.5)

end =#

max_stars_repo_name: HomoModelicus/julia

max_issues_repo_name: HomoModelicus/julia

max_forks_repo_name: HomoModelicus/julia


text:
const REGISTRY_LIST = MirrorUpdater.Types.Registry[
    MirrorUpdater.Types.Registry(
        ;
        owner = "JuliaRegistries",
        name = "General",
        uuid = "23338594-aafe-5451-b93e-139f81909106",
        url = "https://github.com/JuliaRegistries/General",
        ),
    MirrorUpdater.Types.Registry(
        ;
        owner = "fredrikekre",
        name = "Registry",
        uuid = "ae0cb698-197b-42ec-a0a0-4f871aea6013",
        url = "https://github.com/fredrikekre/Registry",
        ),
    ]


text:
@testset "Overorders" begin
   Qx,  x = FlintQQ["x"]

  f = x^3 - 1000*x^2 + 1000*x - 1000

   K,  a = NumberField(f, "a");

  E = EquationOrder(K)

  #@test length(overorders(E)) == 16

  f = x^3 + 31 * x^2 + 43* x + 77

   K,  a = NumberField(f, "a");

  E = EquationOrder(K)

  #@test length(overorders(E)) == 15
  #@test length(Hecke.overorders_naive(E)) == 15
end


text:
using Outliers
using Base.Test

@testset "Tukey Fences" begin
    tukey_fences([1,100,101,102,200]) == (97.0,105.0)
end


text:
# Autogenerated wrapper script for NEO_jll for x86_64-linux-gnu-cxx11
export libigdrcl, libze_intel_gpu, ocloc

using gmmlib_jll
using libigc_jll
using oneAPI_Level_Zero_Headers_jll
JLLWrappers.@generate_wrapper_header("NEO")
JLLWrappers.@declare_library_product(libigdrcl, "libigdrcl.so")
JLLWrappers.@declare_library_product(libze_intel_gpu, "libze_intel_gpu.so.1")
JLLWrappers.@declare_executable_product(ocloc)
function __init__()
    JLLWrappers.@generate_init_header(gmmlib_jll, libigc_jll, oneAPI_Level_Zero_Headers_jll)
    JLLWrappers.@init_library_product(
        libigdrcl,
        "lib64/intel-opencl/libigdrcl.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libze_intel_gpu,
        "lib64/libze_intel_gpu.so.1",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_executable_product(
        ocloc,
        "bin/ocloc",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using HuffmanCode
using Test

@testset "HuffmanCode.jl" begin
    # Write your tests here.
end

max_stars_repo_name: Colcactus/Huffman

max_issues_repo_name: Colcactus/Huffman

max_forks_repo_name: Colcactus/Huffman


text:
"""
    default_symbolics_mapping :: Dict{String,Any}

Default mapping of SBML function names to Julia functions, represented as a
dictionary from Strings (SBML names) to anything `eval`uable as Julia&Symbolics
functions, such as symbols and expressions.

The default mapping only contains the basic SBML functions that are
unambiguously represented in Julia; it is supposed to be extended by the user
if more functions need to be supported.
"""
const default_symbolics_mapping = Dict{String,Any}(
    "*" => :*,
    "+" => :+,
    "-" => :-,
    "/" => :/,
    "abs" => :abs,
    "arccos" => :acos,
    "arccosh" => :acosh,
    "arccot" => :acot,
    "arccoth" => :acoth,
    "arccsc" => :acsc,
    "arccsch" => :acsch,
    "arcsec" => :asec,
    "arcsech" => :asech,
    "arcsin" => :asin,
    "arcsinh" => :asinh,
    "arctan" => :atan,
    "arctanh" => :atanh,
    "ceiling" => :ceil,
    "cos" => :cos,
    "cosh" => :cosh,
    "cot" => :cot,
    "coth" => :coth,
    "csc" => :csc,
    "csch" => :csch,
    "eq" => :isequal,
    "exp" => :exp,
    "factorial" => :factorial,
    "floor" => :floor,
    "geq" => :>=,
    "gt" => :>,
    "leq" => :<=,
    "ln" => :log,
    "log" => :sbmlLog,
    "lt" => :<,
    "piecewise" => :(Core.ifelse),
    "power" => :^,
    "root" => :sbmlRoot,
    "sech" => :sech,
    "sec" => :sec,
    "sinh" => :sinh,
    "sin" => :sin,
    "tanh" => :tanh,
    "tan" => :tan,
)

sbmlLog(x) = log(x, 10)
sbmlLog(base, x) = log(base, x)

sbmlRoot(x) = sqrt(x)
sbmlRoot(power, x) = x^(1 / power)

allowed_sym(x, allowed_funs) =
    haskey(allowed_funs, x) ? allowed_funs[x] :
    throw(DomainError(x, "Unknown SBML function"))

"""
    const default_symbolics_constants::Dict{String, Any}

A dictionary of default constants filled in place of SBML Math constants in the
symbolics conversion.
"""
const default_symbolics_constants =
    Dict{String,Any}("true" => true, "false" => false, "pi" => pi, "e" => exp(1))

"""
    Base.convert(
        ::Type{Num},
        x::SBML.Math;
        mapping = default_symbolics_mapping,
        convert_time = (x::SBML.MathTime) -> Num(Variable(Symbol(x.id))).val,
        convert_const = (x::SBML.MathConst) -> Num(default_symbolics_constants[x.id]),
    )

Convert SBML.[`Math`](@ref) to `Num` type from Symbolics package. The
conversion of functions can be customized by supplying a custom mapping; if
nothing is supplied, [`default_symbolics_mapping`](@ref) that translates basic
functions to their Julia equivalents is assumed.

Translation of [`MathLambda`](@ref) is not supported by Symbolics.

[`MathTime`](@ref) is handled specially, the function from the argument
`convert_time` is called to possibly specify any desired behavior. By default,
it just creates a variable with the same name as the time variable name stored
in SBML.
"""
function Base.convert(
    ::Type{Num},
    x::SBML.Math;
    mapping = default_symbolics_mapping,
    convert_time = (x::SBML.MathTime) -> Num(Variable(Symbol(x.id))).val,
    convert_const = (x::SBML.MathConst) -> Num(default_symbolics_constants[x.id]),
)
    conv(x::SBML.MathApply) = eval(allowed_sym(x.fn, mapping))(conv.(x.args)...)
    conv(x::SBML.MathTime) = convert_time(x)
    conv(x::SBML.MathConst) = convert_const(x)
    conv(x::SBML.MathIdent) = Num(Variable(Symbol(x.id))).val
    conv(x::SBML.MathVal) = Num(x.val)
    conv(x::SBML.MathLambda) = throw(DomainError(x, "can't translate lambdas to symbolics"))
    conv(x)
end


text:
using Test;

using AdalmPluto;

@testset "libIIO/channel.jl" begin
    toggleNoAssertions(true);

    try
        global uri          = AdalmPluto.scan("usb", 1, false);
        global ctx          = AdalmPluto.createContext(uri);
        global tx, rx       = AdalmPluto.findTRXDevices(ctx);
        global txc, rxc     = AdalmPluto.findTRXChannels(ctx);
        global txloc, rxloc = AdalmPluto.findLoChannels(ctx);
        global txi, txq     = AdalmPluto.findIQChannels(tx, "voltage0", "voltage1", true);
        global rxi, rxq     = AdalmPluto.findIQChannels(rx, "voltage0", "voltage1", false);
        global channels     = [rxc, rxloc, rxi, rxq, txc, txloc, txi, txq];
    catch
        @error "Could not get a valid uri or context or channels, skipping the remaining tests.";
        C_iio_context_destroy(ctx);
        return;
    end

    # uncomment the following section to show all the available attributes and properties
    #  for channel in channels
        #  @show channel;
        #  nb_attrs = C_iio_channel_get_attrs_count(channel);
        #  for idx = 0:nb_attrs-1
            #  @show attr = C_iio_channel_get_attr(channel, UInt32(idx));
            #  @show C_iio_channel_attr_read(channel, attr);
            #  @show C_iio_channel_attr_read_bool(channel, attr);
            #  @show C_iio_channel_attr_read_double(channel, attr);
            #  @show C_iio_channel_attr_read_longlong(channel, attr);
            #  @show C_iio_channel_attr_get_filename(channel, attr);
            #  @show C_iio_channel_get_device(channel);
            #  @show C_iio_channel_get_id(channel);
            #  @show C_iio_channel_get_type(channel);
            #  @show C_iio_channel_get_name(channel);
            #  @show C_iio_channel_get_modifier(channel);
            #  @show C_iio_channel_is_enabled(channel);
            #  @show C_iio_channel_is_output(channel);
            #  @show C_iio_channel_is_scan_element(channel);
        #  end
    #  end

    @testset "Channel properties" begin
        #  C_iio_channel_get_device (channel::Ptr, …)
        @test C_iio_channel_get_device(rxc) != C_NULL;
        #  C_iio_channel_get_id (channel::Ptr, …)
        @test C_iio_channel_get_id(rxc) == "voltage0";
        #  C_iio_channel_get_modifier (channel::Ptr, …)
        @test C_iio_channel_get_modifier(rxc) == AdalmPluto.libIIO_jl.IIO_NO_MOD;
        #  C_iio_channel_get_name (channel::Ptr, …)
        @test C_iio_channel_get_name(rxloc) == "RX_LO";
        #  C_iio_channel_get_type (channel::Ptr, …)
        @test C_iio_channel_get_type(rxc) == AdalmPluto.libIIO_jl.IIO_VOLTAGE;
        #  C_iio_channel_is_output (channel::Ptr, …)
        @test C_iio_channel_is_output(rxc) == false;
        @test C_iio_channel_is_output(txc) == true;
        #  C_iio_channel_is_scan_element (channel::Ptr, …)
        @test C_iio_channel_is_scan_element(rxc) == false;
        @test C_iio_channel_is_scan_element(rxi) == true;

        #  C_iio_channel_disable (channel::Ptr, …)
        #  C_iio_channel_enable (channel::Ptr, …)
        #  C_iio_channel_is_enabled (channel::Ptr, …)
        enabled = C_iio_channel_is_enabled(rxi);
        if enabled
            @test C_iio_channel_disable(rxi) === nothing;
            @test C_iio_channel_is_enabled(rxi) == false;
            @test C_iio_channel_enable(rxi) === nothing;
            @test C_iio_channel_is_enabled(rxi) == true;
        else
            @test C_iio_channel_enable(rxi) === nothing;
            @test C_iio_channel_is_enabled(rxi) == true;
            @test C_iio_channel_disable(rxi) === nothing;
            @test C_iio_channel_is_enabled(rxi) == false;
        end

        saved_data_ptr = C_iio_channel_get_data(rxc);
        test_data = b"Test data\0";
        #  C_iio_channel_set_data (channel::Ptr, …)
        @test C_iio_channel_set_data(rxc, Ptr{Cvoid}(pointer(test_data))) === nothing;
        #  C_iio_channel_get_data (channel::Ptr, …)
        @test C_iio_channel_get_data(rxc) == Ptr{Cvoid}(pointer(test_data));
        C_iio_channel_set_data(rxc, saved_data_ptr);
    end

    @testset "Channel attributes" begin
        #  C_iio_channel_get_attrs_count (channel::Ptr, …)
        @test (global nb_attrs = C_iio_channel_get_attrs_count(rxc)) > 0;

        if nb_attrs < 1
            @warn "No channel attributes, skipping the related tests";
        else
            #  C_iio_channel_get_attr (channel::Ptr, …)
            @test (global first_attr = C_iio_channel_get_attr(rxc, UInt32(0))) != "";
            #  C_iio_channel_find_attr (channel::Ptr, …)
            @test C_iio_channel_find_attr(rxc, first_attr) == first_attr;
            #  C_iio_channel_attr_get_filename (channel::Ptr, …)
            @test C_iio_channel_attr_get_filename(rxc, "hardwaregain") == "in_voltage0_hardwaregain";

            attr = "gain_control_mode";
            #  C_iio_channel_attr_read (channel::Ptr, …)
            @test C_iio_channel_attr_read(rxc, attr)[1] > 0;
            @test C_iio_channel_attr_read(rxc, "dummy")[1] < 0;
            saved_value = C_iio_channel_attr_read(rxc, attr)[2];
            #  C_iio_channel_attr_write (channel::Ptr, …)
            @test C_iio_channel_attr_write(rxc, attr, "manual") > 0 && C_iio_channel_attr_read(rxc, attr)[2] == "manual";
            @test C_iio_channel_attr_write(rxc, attr, "slow_attack") > 0 && C_iio_channel_attr_read(rxc, attr)[2] == "slow_attack";
            new_gainmode = b"manual\0";
            #  C_iio_channel_attr_write_raw (channel::Ptr, …)
            @test (C_iio_channel_attr_write_raw(rxc, attr, Ptr{Cvoid}(pointer(new_gainmode)), UInt64(sizeof(new_gainmode))) > 0
                && C_iio_channel_attr_read(rxc, attr)[2] == "manual");
            # put back the original value
            C_iio_channel_attr_write(rxc, attr, saved_value);

            #  C_iio_channel_attr_read_bool (channel::Ptr, …)
            @test C_iio_channel_attr_read_bool(rxc, "filter_fir_en")[1] == 0;
            ## Doesn't work, don't know why
            #  saved_fir_en = C_iio_channel_attr_read_bool(rxc, "filter_fir_en")[2];
            #  @test (C_iio_channel_attr_write_bool(rxc, "filter_fir_en", !saved_fir_en) == 0
            #      && C_iio_channel_attr_read_bool(rxc, "filter_fir_en")[2] == !saved_fir_en);
            ##
            #  C_iio_channel_attr_write_bool (channel::Ptr, …)
            phy = C_iio_context_find_device(ctx, AdalmPluto.PHY_DEVICE_NAME);
            r, c, a = C_iio_device_identify_filename(phy, "in_out_voltage_filter_fir_en");
            saved_value = C_iio_channel_attr_read_bool(c, a)[2];
            @test_skip (C_iio_channel_attr_write_bool(c, a, !saved_value) == 0
                   & C_iio_channel_attr_read_bool(c, a)[2] == !saved_value);
            # put back the original value
            C_iio_channel_attr_write_bool(c, a, saved_value);

            #  C_iio_channel_attr_read_double (channel::Ptr, …)
            @test C_iio_channel_attr_read_double(txc, "hardwaregain")[1] == 0;
            saved_value = C_iio_channel_attr_read_double(txc, "hardwaregain")[2];
            #  C_iio_channel_attr_write_double (channel::Ptr, …)
            # assuming saved_gain > minimum available gain
            @test (C_iio_channel_attr_write_double(txc, "hardwaregain", saved_value - 0.25) == 0
                && C_iio_channel_attr_read_double(txc, "hardwaregain")[2] == saved_value - 0.25);
            # put back the original value
            C_iio_channel_attr_write_double(txc, "hardwaregain", saved_value);

            #  C_iio_channel_attr_read_longlong (channel::Ptr, …)
            @test C_iio_channel_attr_read_longlong(rxc, "rf_bandwidth")[1] == 0;
            saved_value = C_iio_channel_attr_read_longlong(rxc, "rf_bandwidth")[2];
            #  C_iio_channel_attr_write_longlong (channel::Ptr, …)
            @test (C_iio_channel_attr_write_longlong(rxc, "rf_bandwidth", 200000)[1] == 0
                && C_iio_channel_attr_read_longlong(rxc, "rf_bandwidth")[2] == 200000);
            sleep(1)
            @test (C_iio_channel_attr_write_longlong(rxc, "rf_bandwidth", 4000000)[1] == 0)

            # put back the original value
            C_iio_channel_attr_write_longlong(rxc, "rf_bandwidth", saved_value);
        end
    end

    # Those functions are not yet implemented as they need either data using the correct format or a function pointer.
    # C_iio_channel_read_raw, C_iio_channel_write, and C_iio_channel_write_raw also need a buffer initialized.
    @testset "Placeholders" begin
        #  C_iio_channel_attr_read_all
        @test_skip C_iio_channel_attr_read_all() === nothing;
        #  C_iio_channel_attr_write_all
        @test_skip C_iio_channel_attr_write_all() === nothing;
        #  C_iio_channel_read_raw
        @test_skip C_iio_channel_read_raw() === nothing;
        #  C_iio_channel_write
        @test_skip C_iio_channel_write() === nothing;
        #  C_iio_channel_write_raw
        @test_skip C_iio_channel_write_raw() === nothing;
    end

    #  C_iio_channel_read! (chn::Ptr, …)
    # This function is not tested here as a working buffer is needed

    C_iio_context_destroy(ctx);
end


text:
include("../src/pressurecorrelations.jl")
include("../src/tempcorrelations.jl")

@testset "Temperature correlations" begin

q_o = 2219
q_w = 11
APIoil = 24
sg_water = 1
GLR = 1.762 * 1000^2 / (2219 + 11)
sg_gas = 0.7
id = 2.992
whp = 150 #psig
A = Shiu_Beggs_relaxationfactor(q_o, q_w, GLR, APIoil, sg_water, sg_gas, id, whp)

@test A ≈ 2089 atol = 1

g_g = 0.0106 * 100
bht = 173
z = 5985 - 0

@test Ramey_temp(z, bht, A, g_g) ≈ 130 atol = 2
@test Ramey_temp(1, bht, A, g_g) ≈ 173 atol = 0.1

#= visual test
depths = range(1, stop = z, length = 100) ;

test_temps = Ramey_wellboretemp.(depths, 0, bht, A, g_g)
depths_plot = z .- depths |> collect

using Gadfly

plot(y = depths_plot, x = test_temps, Geom.line, Coord.cartesian(yflip = true))
=#

end #Temperature correlations


text:
export Hamming
include("Hamming.jl")

max_stars_repo_name: fercarozzi/myseismicjulia

max_issues_repo_name: fercarozzi/myseismicjulia

max_forks_repo_name: fercarozzi/myseismicjulia


text:
@testset "Levenberg-Marquardt: NIST Dataset (Eckerle4)" begin
    print_head("Levenberg-Marquardt: NIST Dataset (Eckerle4)")

    obj = optimizer(nd=3, ny=35, method="lm")
    fun = Eckerle4()
    sol = Vector{Float64}(undef, 3)

    parμ = [1.0, 10.0, 460.0]
    βvec = zeros(Float64, 35)
    βmat = zeros(Float64, 35, 35)

    print_body("\033[1m\033[34mLikelihood Precision Matrix: Identity\033[0m")

    @simd for i in eachindex(βvec)
        @inbounds βvec[i] = 1.0
    end

    @simd for i in eachindex(βvec)
        @inbounds βmat[i,i] = 1.0
    end

    ans = [1.55438286, 4.08883291, 451.541218]

    print_body(string("Precision Const.: ", optimize!(sol, fun, parμ, 1.0, obj)))
    @inbounds for i in eachindex(ans)
        @test sol[i] ≈ ans[i] rtol=1e-7
    end

    print_body(string("Precision Vector: ", optimize!(sol, fun, parμ, βvec, obj)))
    @inbounds for i in eachindex(ans)
        @test sol[i] ≈ ans[i] rtol=1e-7
    end

    print_body(string("Precision Matrix: ", optimize!(sol, fun, parμ, βmat, obj)))
    @inbounds for i in eachindex(ans)
        @test sol[i] ≈ ans[i] rtol=1e-7
    end

    print_body("\033[1m\033[34mLikelihood Precision Matrix: Custom\033[0m")

    @simd for i in eachindex(βvec)
        @inbounds βvec[i] = 1.0 + 0.1 * (i-1)
    end

    @simd for i in eachindex(βvec)
        @inbounds βmat[i,i] = 1.0 + 0.1 * (i-1)
    end

    @inbounds ans[1], ans[2], ans[3] = 1.55104316, 4.07147586, 451.543588

    print_body(string("Precision Vector: ", optimize!(sol, fun, parμ, βvec, obj)))
    @inbounds for i in eachindex(ans)
        @test sol[i] ≈ ans[i] rtol=1e-7
    end

    print_body(string("Precision Matrix: ", optimize!(sol, fun, parμ, βmat, obj)))
    @inbounds for i in eachindex(ans)
        @test sol[i] ≈ ans[i] rtol=1e-7
    end
end


text:
using Distributed, Pkg
Pkg.activate(joinpath(@__DIR__, ".."))
if haskey(ENV, "FOOD_CACHING_PARALLEL")
    nprocs() < Sys.CPU_THREADS && addprocs(Sys.CPU_THREADS - nprocs(),
                                           exeflags = "--project",
                                           dir = joinpath(@__DIR__, ".."))
end
using FoodCachingModels
import FoodCachingFitting: fit, parse_args

fit(; parse_args()...)


text:
using Documenter, Pope, Pope.LJH

makedocs(
    sitename = "Pope.jl",
    modules = [Pope, LJH],
    format = Documenter.HTML(),
    pages = ["Home" => "index.md",
             "Scripts" => "scripts.md",
             "LJH" => "ljh.md",
             "Pope Internals" => "pope.md",
             "News" => "news.md",
             "Noise Analysis" => "noise.md"
             ]
)

deploydocs(
    repo   = "github.com/usnistgov/Pope.jl.git",
    deps   = nothing,
    target = "build",
    make= nothing
)

println("finished make.jl")


text:
using  jInv.Mesh
using  jInv.Utils
using  jInv.LinearSolvers
using  jInv.InverseSolve
using  jInvVis
using  EikonalInv
using  MAT
using  FWI
using  ForwardHelmholtz
using  PyPlot
#############################################################################################################

####################################################################################################################################
####################################################################################################################################
#######################   Problem definition and data creation #####################################################################
####################################################################################################################################
####################################################################################################################################

include("ex2DFWI/readModelAndGenerateMeshMref.jl");
include("ex2DFWI/prepareFWIDataFiles.jl");

m = readdlm("model/SEGmodel2Dsalt.dat");
m = m*1e-3;
m = m';

newSize          = [256,128];
pad     	     = 10;
ABLPad 		     = pad + 8;
jumpSrc 	 	 = 5
maxBatchSize     = 256;
omega   	     = [0.5,0.75,1.25,1.75]*2*pi;

offset  = ceil(Int64,(newSize[1]*(10.0/13.5)));
println("Offset is: ",offset)
domain = [0.0,13.5,0.0,4.2];

(m,Minv,mref,boundsHigh,boundsLow) = readModelAndGenerateMeshMref(m,pad,newSize,domain);

	limits = [1.5,4.5];
	figure(1,figsize = (22,10))
	plotModel(m,true,Minv,pad,limits);
	figure(2,figsize = (22,10))# ,figsize = (22,10)
	plotModel(mref,true,Minv,pad,limits);

useFilesForFields = false;

# ###################################################################################################################
dataFilenamePrefix = "ex2DFWI/DATA_SEG2D";

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

######################## DIRECT SOLVER #################################################
numCores 	= 4;
BLAS.set_num_threads(numCores);
Ainv = getMUMPSsolver([],0,0,2); # Alternatives: Ainv = getJuliaSolver();

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

println("omega*maximum(h): ",omega*maximum(Minv.h)*sqrt(maximum(1./(boundsLow.^2))));

# This is a list of workers for FWI. Ideally they should be on different machines.
workersFWI = [workers()[1]];
prepareFWIDataFiles(m,Minv,mref,boundsHigh,boundsLow,dataFilenamePrefix,omega,one(Complex128)*ones(size(omega)), pad,ABLPad,jumpSrc,
					offset,workersFWI,maxBatchSize,Ainv,useFilesForFields);

####################################################################################################################################
####################################################################################################################################
####################### Inversion setup for a single machine  ######################################################################
####################################################################################################################################
####################################################################################################################################
dataFilenamePrefix = "ex2DFWI/DATA_SEG2D";
file = matread(string(dataFilenamePrefix,"_PARAM.mat"));
n_cells 		= file["n"];
domain 	        = file["domain"];
gamma 			= file["gamma"];
omega 			= file["omega"];
waveCoef 		= file["waveCoef"];
boundsLow 		= file["boundsLow"];
boundsHigh 		= file["boundsHigh"];
mref 			= file["mref"];
file = 0;

Minv = getRegularMesh(domain,n_cells);
if length(omega)==1 # matlab saves a 1 variable array as scalar.
	omega = [omega];
	waveCoef = [waveCoef];
end

### Read receivers and sources files
RCVfile = string(dataFilenamePrefix,"_rcvMap.dat");
SRCfile = string(dataFilenamePrefix,"_srcMap.dat");
srcNodeMap = readSrcRcvLocationFile(SRCfile,Minv);
rcvNodeMap = readSrcRcvLocationFile(RCVfile,Minv);
Q = generateSrcRcvProjOperators(Minv.n+1,srcNodeMap);
Q = Q.*1/(norm(Minv.h)^2);
println("We have ",size(Q,2)," sources");
P = generateSrcRcvProjOperators(Minv.n+1,rcvNodeMap);

### Read the data files to an array of array pointers
Wd   = Array{Array{Complex128,2}}(length(omega))
dobs = Array{Array{Complex128,2}}(length(omega))
for k = 1:length(omega)
	omRound = string(round((omega[k]/(2*pi))*100.0)/100.0);
	(DobsFWIwk,WdFWIwk) =  readDataFileToDataMat(string(dataFilenamePrefix,"_freq",omRound,".dat"),srcNodeMap,rcvNodeMap);
	Wd[k] 	= WdFWIwk;
	dobs[k] = DobsFWIwk;
end

# setup active cells and background for the inversion.
N = prod(Minv.n+1);
Iact = sparse(I,N,N);
mback   = zeros(Float64,N);

########################################################################################################
##### Set up remote workers ############################################################################
########################################################################################################

batch = size(Q,2);         # We solve all the sources at once.
useFilesForFields = false; # a flag whether to store the fields on the disk.
## Set up a MUMPS direct solver #Ainv = getMUMPSsolver([],0,0,2); #
Ainv = getJuliaSolver();
## Choose the workers for FWI (here, its a single worker)
workersFWI = [1];
## Set up workers and division to tasks per frequencies ######################
(pFor,contDiv,SourcesSubInd) = getFWIparam(omega,waveCoef,vec(gamma),Q,P,Minv,Ainv,workersFWI,batch,useFilesForFields);
misfun = SSDFun;
pMis = getMisfitParam(pFor, Wd, dobs, misfun, Iact,mback);

########################################################################################################
# Set up Inversion #################################################################################
########################################################################################################

## models are usually given in velocity (km/sec). We invert here for the slowness, which is 1/v.
mref 				= velocityToSlow(mref)[1];
t    				= copy(boundsLow);
boundsLow 			= velocityToSlow(boundsHigh)[1];
boundsHigh 			= velocityToSlow(t)[1]; t = 0;
modfun 				= slowToSlowSquared;

maxStep				=0.05*maximum(boundsHigh);

regparams 			= [1.0,1.0,1.0,1e-5];
cgit 				= 8;
alpha 				= 1e-10;
pcgTol 				= 1e-1;
maxit 				= 8;
HesPrec 			= getExactSolveRegularizationPreconditioner();
regfun(m,mref,M) 	= wdiffusionRegNodal(m,mref,M,Iact=Iact,C=regparams);

pInv = getInverseParam(Minv,modfun,regfun,alpha,mref[:],boundsLow,boundsHigh,
                         maxStep=maxStep,pcgMaxIter=cgit,pcgTol=pcgTol,
						 minUpdate=1e-3, maxIter = maxit,HesPrec=HesPrec);

plotting = true;
function plotIntermediateResults(mc,Dc,iter,pInv,PMis,resultsFilename="")
	# Models are usually shown in velocity.
	fullMc = slowSquaredToVelocity(reshape(Iact*pInv.modelfun(mc)[1] + mback,tuple((pInv.MInv.n+1)...)))[1];
	if plotting
		close(888);
		figure(888);
		plotModel(fullMc,false,[],0,[1.5,4.8]);
	end
end

# Run one sweep of a frequency continuation procedure.
mc,Dc,flag,His = freqCont(copy(mref[:]), pInv, pMis,contDiv, 4, "",plotIntermediateResults,"Joint",1,1,"projGN");

subplot(1,3,1)
semilogy(His[end].F,"-o")
xlabel("PGNCG iterations")
title("misfit")

subplot(1,3,2)
semilogy(His[end].Rc,"-o")
xlabel("PGNCG iterations")
title("regularizer");

subplot(1,3,3)
semilogy(His[end].dJ,"-o")
xlabel("PGNCG iterations")
title("norm of proj grad");


text:
if VERSION < v"1.1-DEV"
    isnothing(::Nothing) = true
    isnothing(::Any) = false
end

function error_plot(fun::Function,
                    Rf, rₘₐₓ, Ns::AbstractVector{<:Integer},
                    ρcol::Int,
                    errors::Vector{<:Tuple{Int,String,<:Real,<:Real}},
                    elapsed_col::Int,
                    headers::Vector{String},
                    other_vectors::AbstractVector...;
                    sa = subplot(311),
                    sb = subplot(312),
                    sc = subplot(313),
                    kwargs...)
    plot_fun = (i,data) -> begin
        ρ = data[:,ρcol]

        line = nothing
        for ((ecol,elabel,_,_),style) in zip(errors,Iterators.cycle(["-", ":", "--"]))
            sca(sa)
            line = loglog(ρ,abs.(data[:,ecol]), "s$style",
                          color = !isnothing(line) ? line.get_color() : nothing)[1]
            sca(sc)
            loglog(abs.(data[:,ecol]),data[:,elapsed_col], "s$style",
                   color = !isnothing(line) ? line.get_color() : nothing)
        end

        sca(sb)
        loglog(ρ,data[:,elapsed_col],"s-")
    end

    result = test_convergence_rates(fun, Rf, rₘₐₓ, Ns, ρcol, errors, headers, other_vectors...;
                                    plot_fun=plot_fun, kwargs...)

    sca(sa)
    axes_labels_opposite(:x)
    ylabel("Error")
    xlabel(L"\rho")

    sca(sb)
    no_tick_labels()
    ylabel("Elapsed time [s]")

    sca(sc)
    xlabel("Error")
    ylabel("Elapsed time [s]")
    tight_layout()

    result
end


text:
# ## ---------------------------------------------------------
# function plot_frees(net::MetNet, frees;
#         δ = 40,
#         verbose = true,
#         m = (8, :square),
#         pltkwargs...
#     )
#     freeis = rxnindex(net, frees)
#     frees = rxns(net, freeis)

#     box = BoxGrid(net, freeis, δ)
#     V = length(box)

#     # Start collecting, TODO: make threaded
#     tolf = 0.01
#     tol = Δv(net, freeis) .* tolf

#     boxc = 0
#     polc = 0

#     cid = ("plot_frees", hash(net), δ, hash(freeis))
#     polV = lcache(cid; verbose) do

#         polV_ = [Float64[] for id in frees]

#         prog = Progress(V; dt = 0.1, desc = "Collecting... ")
#         for freev in box
#             fullv = fixxing(net, freeis, freev; tol) do
#                 fba(net)
#             end

#             # info
#             boxc += 1
#             verbose && next!(prog; showvalues = [(:polfrac, polc/boxc)])

#             # feasible
#             isfea = any(fullv .!= 0.0)
#             !isfea && continue

#             # collect
#             for (v, id) in zip(freev, eachindex(frees))
#                 push!(polV_[id], v)
#             end

#             # info
#             polc += 1

#         end
#         verbose && finish!(prog)

#         return polV_
#     end

#     # --------------------------------------------
#     println()
#     verbose && @info("Ploting")
#     ps = Plots.Plot[]
#     defaults = (;label = "", m, alpha = 0.5, stroke = false)

#     for i in 1:length(frees), j in (i + 1):length(frees)
#         freex, freey = frees[i], frees[j]
#         polVx, polVy = polV[i], polV[j]

#         verbose && @show freex, freey

#         p = plot(;xlabel = freex, ylabel = freey)

#         # reduce redundancy
#         pool = Set{Tuple{Float64, Float64}}()
#         curr_vx = nothing
#         for si in sortperm(polVx)
#             vx, vy = polVx[si], polVy[si]

#             if vx != curr_vx # new vx
#                 # plot if not empty
#                 !isempty(pool) && scatter!(p, first.(pool), last.(pool); color = :red,
#                     defaults..., pltkwargs...)

#                 empty!(pool)
#                 curr_vx = vx
#             end
#             push!(pool, (vx, vy))
#         end

#         push!(ps, p)
#     end
#     return ps

# end

## ------------------------------------------------------
function marginal_plots(f::Function, V; normalize = true, bins = length(V))
    ps = Plots.Plot[]
    # marginals
    freeids = free_ids(V)
    ridxs = view(shuffle(eachindex(V)), 1:bins)
    for id in freeids
        Vi_ = Vi(V, id)
        Vi_ = view(Vi_, ridxs)
        P = f.(Vi_)
        normalize && normalizeP!(P)
        p = bar(hist(Vi_, P); label = "", xlabel = string(id), color = :black)
        push!(ps, p)
    end
    ps
end
marginal_plots(V; normalize = true, bins = length(V)) =
    marginal_plots((v) -> 1.0, V; normalize, bins)

## ------------------------------------------------------
function pol_plots(V; bins = 10000)
    frees = free_ids(V)
    ps = Plots.Plot[]
    ridxs = view(shuffle(eachindex(V)), 1:bins)
    for i1 in 1:length(frees)
        for i2 in (i1 + 1):length(frees)
            id1, id2 = frees[i1], frees[i2]
            V1, V2 = Vi(V, id1), Vi(V, id2)
            p = scatter( view(V1, ridxs), view(V2, ridxs);
                label = "",
                xlabel = string(id1),
                ylabel = string(id2)
            )
            push!(ps, p)
        end
    end
    ps
end


text:
using StaticArrays
# -----------------------------------------------------
#
# Datasets from Time Series Analysis by Box and Jenkins
#
# -----------------------------------------------------

const seriesA = SA[17.0, 16.6, 16.3, 16.1, 17.1, 16.9, 16.8, 17.4, 17.1, 17.0, 16.7, 17.4, 17.2, 17.4,
                17.4, 17.0, 17.3, 17.2, 17.4, 16.8, 17.1, 17.4, 17.4, 17.5, 17.4, 17.6, 17.4, 17.3,
                17.0, 17.8, 17.5, 18.1, 17.5, 17.4, 17.4, 17.1, 17.6, 17.7, 17.4, 17.8, 17.6, 17.5,
                16.5, 17.8, 17.3, 17.3, 17.1, 17.4, 16.9, 17.3, 17.6, 16.9, 16.7, 16.8, 16.8, 17.2,
                16.8, 17.6, 17.2, 16.6, 17.1, 16.9, 16.6, 18.0, 17.2, 17.3, 17.0, 16.9, 17.3, 16.8,
                17.3, 17.4, 17.7, 16.8, 16.9, 17.0, 16.9, 17.0, 16.6, 16.7, 16.8, 16.7, 16.4, 16.5,
                16.4, 16.6, 16.5, 16.7, 16.4, 16.4, 16.2, 16.4, 16.3, 16.4, 17.0, 16.9, 17.1, 17.1,
                16.7, 16.9, 16.5, 17.2, 16.4, 17.0, 17.0, 16.7, 16.2, 16.6, 16.9, 16.5, 16.6, 16.6,
                17.0, 17.1, 17.1, 16.7, 16.8, 16.3, 16.6, 16.8, 16.9, 17.1, 16.8, 17.0, 17.2, 17.3,
                17.2, 17.3, 17.2, 17.2, 17.5, 16.9, 16.9, 16.9, 17.0, 16.5, 16.7, 16.8, 16.7, 16.7,
                16.6, 16.5, 17.0, 16.7, 16.7, 16.9, 17.4, 17.1, 17.0, 16.8, 17.2, 17.2, 17.4, 17.2,
                16.9, 16.8, 17.0, 17.4, 17.2, 17.2, 17.1, 17.1, 17.1, 17.4, 17.2, 16.9, 16.9, 17.0,
                16.7, 16.9, 17.3, 17.8, 17.8, 17.6, 17.5, 17.0, 16.9, 17.1, 17.2, 17.4, 17.5, 17.9,
                17.0, 17.0, 17.0, 17.2, 17.3, 17.4, 17.4, 17.0, 18.0, 18.2, 17.6, 17.8, 17.7, 17.2,
                17.4]

const seriesB = SA[460, 457, 452, 459, 462, 459, 463, 479, 493, 490, 492, 498, 499, 497, 496,
                490, 489, 478, 487, 491, 487, 482, 479, 478, 479, 477, 479, 475, 479, 476,
                476, 478, 479, 477, 476, 475, 475, 473, 474, 474, 474, 465, 466, 467, 471,
                471, 467, 473, 481, 488, 490, 489, 489, 485, 491, 492, 494, 499, 498, 500,
                497, 494, 495, 500, 504, 513, 511, 514, 510, 509, 515, 519, 523, 519, 523,
                531, 547, 551, 547, 541, 545, 549, 545, 549, 547, 543, 540, 539, 532, 517,
                527, 540, 542, 538, 541, 541, 547, 553, 559, 557, 557, 560, 571, 571, 569,
                575, 580, 584, 585, 590, 599, 603, 599, 596, 585, 587, 585, 581, 583, 592,
                592, 596, 596, 595, 598, 598, 595, 595, 592, 588, 582, 576, 578, 589, 585,
                580, 579, 584, 581, 581, 577, 577, 578, 580, 586, 583, 581, 576, 571, 575,
                575, 573, 577, 582, 584, 579, 572, 577, 571, 560, 549, 556, 557, 563, 564,
                567, 561, 559, 553, 553, 553, 547, 550, 544, 541, 532, 525, 542, 555, 558,
                551, 551, 552, 553, 557, 557, 548, 547, 545, 545, 539, 539, 535, 537, 535,
                536, 537, 543, 548, 546, 547, 548, 549, 553, 553, 552, 551, 550, 553, 554,
                551, 551, 545, 547, 547, 537, 539, 538, 533, 525, 513, 510, 521, 521, 521,
                523, 516, 511, 518, 517, 520, 519, 519, 519, 518, 513, 499, 485, 454, 462,
                473, 482, 486, 475, 459, 451, 453, 446, 455, 452, 457, 449, 450, 435, 415,
                398, 399, 361, 383, 393, 385, 360, 364, 365, 370, 374, 359, 335, 323, 306,
                333, 330, 336, 328, 316, 320, 332, 320, 333, 344, 339, 350, 351, 350, 345,
                350, 359, 375, 379, 376, 382, 370, 365, 367, 372, 373, 363, 371, 369, 376,
                387, 387, 376, 385, 385, 380, 373, 382, 377, 376, 379, 386, 387, 386, 389,
                394, 393, 409, 411, 409, 408, 393, 391, 388, 396, 387, 383, 388, 382, 384,
                382, 383, 383, 388, 395, 392, 386, 383, 377, 364, 369, 355, 350, 353, 340,
                350, 349, 358, 360, 360, 366, 359, 356, 355, 367, 357, 361, 355, 348, 343,
                330, 340, 339, 331, 345, 352, 346, 352, 357]

const seriesB2 = SA[445, 448, 450, 447, 451, 453, 454, 454, 459, 440, 446, 443, 443, 440, 439, 435,
                435, 436, 435, 435, 435, 433, 429, 428, 425, 427, 425, 422, 409, 407, 423, 422,
                417, 421, 424, 414, 419, 429, 426, 425, 424, 425, 425, 424, 425, 421, 414, 410,
                411, 406, 406, 413, 411, 410, 405, 409, 410, 405, 401, 401, 401, 414, 419, 425,
                423, 411, 414, 420, 412, 415, 412, 412, 411, 412, 409, 407, 408, 415, 413, 413,
                410, 405, 410, 412, 413, 411, 411, 409, 406, 407, 410, 408, 408, 409, 410, 409,
                405, 406, 405, 407, 409, 407, 409, 425, 425, 428, 436, 442, 442, 433, 435, 433,
                435, 429, 439, 437, 439, 438, 435, 433, 437, 437, 444, 441, 440, 441, 439, 439,
                438, 437, 441, 442, 441, 437, 427, 423, 424, 428, 428, 431, 425, 423, 420, 426,
                418, 416, 419, 418, 416, 419, 425, 421, 422, 422, 417, 420, 417, 418, 419, 419,
                417, 419, 422, 423, 422, 421, 421, 419, 418, 421, 420, 413, 413, 408, 409, 415,
                415, 420, 420, 424, 426, 423, 423, 425, 431, 436, 436, 440, 436, 443, 445, 439,
                443, 445, 450, 461, 471, 467, 462, 456, 464, 463, 465, 464, 456, 460, 458, 453,
                453, 449, 447, 453, 450, 459, 457, 453, 455, 453, 450, 456, 461, 463, 463, 461,
                465, 473, 473, 475, 499, 485, 491, 496, 504, 504, 509, 511, 524, 525, 541, 531,
                529, 530, 531, 527, 525, 519, 514, 509, 505, 513, 525, 519, 519, 522, 522]

const seriesC = SA[26.6, 27.0, 27.1, 27.1, 27.1, 27.1, 26.9, 26.8, 26.7, 26.4, 26.0, 25.8, 25.6, 25.2,
                25.0, 24.6, 24.2, 24.0, 23.7, 23.4, 23.1, 22.9, 22.8, 22.7, 22.6, 22.4, 22.2, 22.0,
                21.8, 21.4, 20.9, 20.3, 19.7, 19.4, 19.3, 19.2, 19.1, 19.0, 18.9, 18.9, 19.2, 19.3,
                19.3, 19.4, 19.5, 19.6, 19.6, 19.6, 19.6, 19.6, 19.7, 19.9, 20.0, 20.1, 20.2, 20.3,
                20.6, 21.6, 21.9, 21.7, 21.3, 21.2, 21.4, 21.7, 22.2, 23.0, 23.8, 24.6, 25.1, 25.6,
                25.8, 26.1, 26.3, 26.3, 26.2, 26.0, 25.8, 25.6, 25.4, 25.2, 24.9, 24.7, 24.5, 24.4,
                24.4, 24.4, 24.4, 24.4, 24.3, 24.4, 24.4, 24.4, 24.4, 24.4, 24.5, 24.5, 24.4, 24.3,
                24.2, 24.2, 24.0, 23.9, 23.7, 23.6, 23.5, 23.5, 23.5, 23.5, 23.5, 23.7, 23.8, 23.8,
                23.9, 23.9, 23.8, 23.7, 23.6, 23.4, 23.2, 23.0, 22.8, 22.6, 22.4, 22.0, 21.6, 21.3,
                21.2, 21.2, 21.1, 21.0, 20.9, 21.0, 21.0, 21.1, 21.2, 21.1, 20.9, 20.8, 20.8, 20.8,
                20.8, 20.9, 20.8, 20.8, 20.7, 20.7, 20.8, 20.9, 21.2, 21.4, 21.7, 21.8, 21.9, 22.2,
                22.5, 22.8, 23.1, 23.4, 23.8, 24.1, 24.6, 24.9, 24.9, 25.1, 25.0, 25.0, 25.0, 25.0,
                24.9, 24.8, 24.7, 24.6, 24.5, 24.5, 24.5, 24.5, 24.5, 24.5, 24.5, 24.4, 24.4, 24.2,
                24.2, 24.1, 24.1, 24.0, 24.0, 24.0, 23.9, 23.8, 23.8, 23.7, 23.7, 23.6, 23.7, 23.6,
                23.6, 23.6, 23.5, 23.5, 23.4, 23.3, 23.3, 23.3, 23.4, 23.4, 23.3, 23.2, 23.3, 23.3,
                23.2, 23.1, 22.9, 22.8, 22.6, 22.4, 22.2, 21.8, 21.3, 20.8, 20.2, 19.7, 19.3, 19.1,
                19.0, 18.8]

const seriesD = SA[8.0, 8.0, 7.4, 8.0, 8.0, 8.0, 8.0, 8.8, 8.4, 8.4, 8.0, 8.2, 8.2, 8.2,
                8.4, 8.4, 8.4, 8.6, 8.8, 8.6, 8.6, 8.6, 8.6, 8.6, 8.8, 8.9, 9.1, 9.5,
                8.5, 8.4, 8.3, 8.2, 8.1, 8.3, 8.4, 8.7, 8.8, 8.8, 9.2, 9.6, 9.0, 8.8,
                8.6, 8.6, 8.8, 8.8, 8.6, 8.6, 8.4, 8.3, 8.4, 8.3, 8.3, 8.1, 8.2, 8.3,
                8.5, 8.1, 8.1, 7.9, 8.3, 8.1, 8.1, 8.1, 8.4, 8.7, 9.0, 9.3, 9.3, 9.5,
                9.3, 9.5, 9.5, 9.5, 9.5, 9.5, 9.5, 9.9, 9.5, 9.7, 9.1, 9.1, 8.9, 9.3,
                9.1, 9.1, 9.3, 9.5, 9.3, 9.3, 9.3, 9.9, 9.7, 9.1, 9.3, 9.5, 9.4, 9.0,
                9.0, 8.8, 9.0, 8.8, 8.6, 8.6, 8.0, 8.0, 8.0, 8.0, 8.6, 8.0, 8.0, 8.0,
                7.6, 8.6, 9.6, 9.6, 10.0, 9.4, 9.3, 9.2, 9.5, 9.5, 9.5, 9.9, 9.9, 9.5,
                9.3, 9.5, 9.5, 9.1, 9.3, 9.5, 9.3, 9.1, 9.3, 9.1, 9.5, 9.4, 9.5, 9.6,
                10.2, 9.8, 9.6, 9.6, 9.4, 9.4, 9.4, 9.4, 9.6, 9.6, 9.4, 9.4, 9.0, 9.4,
                9.4, 9.6, 9.4, 9.2, 8.8, 8.8, 9.2, 9.2, 9.6, 9.6, 9.8, 9.8, 10.0, 10.0,
                9.4, 9.8, 8.8, 8.8, 8.8, 8.8, 9.6, 9.6, 9.6, 9.2, 9.2, 9.0, 9.0, 9.0,
                9.4, 9.0, 9.0, 9.4, 9.4, 9.6, 9.4, 9.6, 9.6, 9.6, 10.0, 10.0, 9.6, 9.2,
                9.2, 9.2, 9.0, 9.0, 9.6, 9.8, 10.2, 10.0, 10.0, 10.0, 9.4, 9.2, 9.6, 9.7,
                9.7, 9.8, 9.8, 9.8, 10.0, 10.0, 8.6, 9.0, 9.4, 9.4, 9.4, 9.4, 9.4, 9.6,
                10.0, 10.0, 9.8, 9.8, 9.7, 9.6, 9.4, 9.2, 9.0, 9.4, 9.6, 9.6, 9.6, 9.6,
                9.6, 9.6, 9.0, 9.4, 9.4, 9.4, 9.6, 9.4, 9.6, 9.6, 9.8, 9.8, 9.8, 9.6,
                9.2, 9.6, 9.2, 9.2, 9.6, 9.6, 9.6, 9.6, 9.6, 9.6, 10.0, 10.0, 10.4, 10.4,
                9.8, 9.0, 9.6, 9.8, 9.6, 8.6, 8.0, 8.0, 8.0, 8.0, 8.4, 8.8, 8.4, 8.4,
                9.0, 9.0, 9.4, 10.0, 10.0, 10.0, 10.2, 10.0, 10.0, 9.6, 9.0, 9.0, 8.6,
                9.0, 9.6, 9.6, 9.0, 9.0, 8.9, 8.8, 8.7, 8.6, 8.3, 7.9, 8.5, 8.7, 8.9, 9.1,
                9.1, 9.1]

const seriesE = SA[101, 82, 66, 35, 31, 7, 20, 92, 154, 125,
                85, 68, 38, 23, 10, 24, 83, 132, 131, 118,
                90, 67, 60, 47, 41, 21, 16, 6, 4, 7,
                14, 34, 45, 43, 48, 42, 28, 10, 8, 2,
                0, 1, 5, 12, 14, 35, 46, 41, 30, 24,
                16, 7, 4, 2, 8, 17, 36, 50, 62, 67,
                71, 48, 28, 8, 13, 57, 122, 138, 103, 86,
                63, 37, 24, 11, 15, 40, 62, 98, 124, 96,
                66, 64, 54, 39, 21, 7, 4, 23, 55, 94,
                96, 77, 59, 44, 47, 30, 16, 7, 37, 74]

const seriesF = SA[47, 64, 23, 71, 38, 64, 55, 41, 59, 48, 71, 35, 57, 40, 58,
                44, 80, 55, 37, 74, 51, 57, 50, 60, 45, 57, 50, 45, 25, 59,
                50, 71, 56, 74, 50, 58, 45, 54, 36, 54, 48, 55, 45, 57, 50,
                62, 44, 64, 43, 52, 38, 59, 55, 41, 53, 49, 34, 35, 54, 45,
                68, 38, 50, 60, 39, 59, 40, 57, 54, 23]

const seriesG = SA[112, 118, 132, 129, 121, 135, 148, 148, 136, 119, 104, 118,
                115, 126, 141, 135, 125, 149, 170, 170, 158, 133, 114, 140,
                145, 150, 178, 163, 172, 178, 199, 199, 184, 162, 146, 166,
                171, 180, 193, 181, 183, 218, 230, 242, 209, 191, 172, 194,
                196, 196, 236, 235, 229, 243, 264, 272, 237, 211, 180, 201,
                204, 188, 235, 227, 234, 264, 302, 293, 259, 229, 203, 229,
                242, 233, 267, 269, 270, 315, 364, 347, 312, 274, 237, 278,
                284, 277, 317, 313, 318, 374, 413, 405, 355, 306, 271, 306,
                315, 301, 356, 348, 355, 422, 465, 467, 404, 347, 305, 336,
                340, 318, 362, 348, 363, 435, 491, 505, 404, 359, 310, 337,
                360, 342, 406, 396, 420, 472, 548, 559, 463, 407, 362, 405,
                417, 391, 419, 461, 472, 535, 622, 606, 508, 461, 390, 432]

# ------------------------------------------------------
#
# Datasets from Introduction Time Series and Forecasting
# by Brockwell and Davis
#
# ------------------------------------------------------

const dowj = SA[110.94, 110.69, 110.43, 110.56, 110.75, 110.84, 110.46, 110.56, 110.46, 110.05,
            109.6, 109.31, 109.31, 109.25, 109.02, 108.54, 108.77, 109.02, 109.44, 109.38,
            109.53, 109.89, 110.56, 110.56, 110.72, 111.23, 111.48, 111.58, 111.9, 112.19,
            112.06, 111.96, 111.68, 111.36, 111.42, 112, 112.22, 112.7, 113.15, 114.36,
            114.65, 115.06, 115.86, 116.4, 116.44, 116.88, 118.07, 118.51, 119.28, 119.79,
            119.7, 119.28, 119.66, 120.14, 120.97, 121.13, 121.55, 121.96, 122.26, 123.79,
            124.11, 124.14, 123.37, 123.02, 122.86, 123.02, 123.11, 123.05, 123.05, 122.83,
            123.18, 122.67, 122.73, 122.86, 122.67, 122.09, 122, 121.23]

const wine = SA[464, 675, 703, 887, 1139, 1077, 1318, 1260, 1120, 963, 996, 960, 530, 883,
            894, 1045, 1199, 1287, 1565, 1577, 1076, 918, 1008, 1063, 544, 635, 804, 980,
            1018, 1064, 1404, 1286, 1104, 999, 996, 1015, 615, 722, 832, 977, 1270, 1437,
            1520, 1708, 1151, 934, 1159, 1209, 699, 830, 996, 1124, 1458, 1270, 1753, 2258,
            1208, 1241, 1265, 1828, 809, 997, 1164, 1205, 1538, 1513, 1378, 2083, 1357, 1536,
            1526, 1376, 779, 1005, 1193, 1522, 1539, 1546, 2116, 2326, 1596, 1356, 1553, 1613,
            814, 1150, 1225, 1691, 1759, 1754, 2100, 2062, 2012, 1897, 1964, 2186, 966, 1549,
            1538, 1612, 2078, 2137, 2907, 2249, 1883, 1739, 1828, 1868, 1138, 1430, 1809, 1763,
            2200, 2067, 2503, 2141, 2103, 1972, 2181, 2344, 970, 1199, 1718, 1683, 2025, 2051,
            2439, 2353, 2230, 1852, 2147, 2286, 1007, 1665, 1642, 1525, 1838, 1892, 2920, 2572,
            2617, 2047]

const lake = SA[10.38, 11.86, 10.97, 10.8, 9.79, 10.39, 10.42, 10.82, 11.4, 11.32, 11.44, 11.68, 11.17,
            10.53, 10.01, 9.91, 9.14, 9.16, 9.55, 9.67, 8.44, 8.24, 9.1, 9.09, 9.35, 8.82, 9.32, 9.01,
            9, 9.8, 9.83, 9.72, 9.89, 10.01, 9.37, 8.69, 8.19, 8.67, 9.55, 8.92, 8.09, 9.37, 10.13,
            10.14, 9.51, 9.24, 8.66, 8.86, 8.05, 7.79, 6.75, 6.75, 7.82, 8.64, 10.58, 9.48, 7.38,
            6.9, 6.94, 6.24, 6.84, 6.85, 6.9, 7.79, 8.18, 7.51, 7.23, 8.42, 9.61, 9.05, 9.26, 9.22,
            9.38, 9.1, 7.95, 8.12, 9.75, 10.85, 10.41, 9.96, 9.61, 8.76, 8.18, 7.21, 7.13, 9.1, 8.25,
            7.91, 6.89, 5.96, 6.8, 7.68, 8.38, 8.52, 9.74, 9.31, 9.89, 9.96, ]


text:
module LowStorageRungeKuttaMethod
export LowStorageRungeKutta2N
export LSRK54CarpenterKennedy, LSRK144NiegemannDiehlBusch, LSRKEulerMethod

using GPUifyLoops
include("LowStorageRungeKuttaMethod_kernels.jl")

using ..ODESolvers
ODEs = ODESolvers
using ..SpaceMethods
using ..MPIStateArrays: device, realview

"""
    LowStorageRungeKutta2N(f, RKA, RKB, RKC, Q; dt, t0 = 0)

This is a time stepping object for explicitly time stepping the differential
equation given by the right-hand-side function `f` with the state `Q`, i.e.,

```math
  \\dot{Q} = f(Q, t)
```

with the required time step size `dt` and optional initial time `t0`.  This
time stepping object is intended to be passed to the `solve!` command.

The constructor builds a low-storage Runge-Kutta scheme using 2N
storage based on the provided `RKA`, `RKB` and `RKC` coefficient arrays.

The available concrete implementations are:

  - [`LSRK54CarpenterKennedy`](@ref)
  - [`LSRK144NiegemannDiehlBusch`](@ref)
"""
mutable struct LowStorageRungeKutta2N{T, RT, AT, Nstages} <: ODEs.AbstractODESolver
  "time step"
  dt::RT
  "time"
  t::RT
  "rhs function"
  rhs!
  "Storage for RHS during the LowStorageRungeKutta update"
  dQ::AT
  "low storage RK coefficient vector A (rhs scaling)"
  RKA::NTuple{Nstages, RT}
  "low storage RK coefficient vector B (rhs add in scaling)"
  RKB::NTuple{Nstages, RT}
  "low storage RK coefficient vector C (time scaling)"
  RKC::NTuple{Nstages, RT}

  function LowStorageRungeKutta2N(rhs!, RKA, RKB, RKC,
                                  Q::AT; dt=0, t0=0) where {AT<:AbstractArray}

    T = eltype(Q)
    RT = real(T)

    dQ = similar(Q)
    fill!(dQ, 0)

    new{T, RT, AT, length(RKA)}(RT(dt), RT(t0), rhs!, dQ, RKA, RKB, RKC)
  end
end

function LowStorageRungeKutta2N(spacedisc::AbstractSpaceMethod, RKA, RKB, RKC,
                                Q::AT; dt=0, t0=0) where {AT<:AbstractArray}
  rhs! = (x...; increment) -> SpaceMethods.odefun!(spacedisc, x..., increment = increment)
  LowStorageRungeKutta2N(rhs!, RKA, RKB, RKC, Q; dt=dt, t0=t0)
end

ODEs.updatedt!(lsrk::LowStorageRungeKutta2N, dt) = (lsrk.dt = dt)
ODEs.updatetime!(lsrk::LowStorageRungeKutta2N, time) = (lsrk.t = time)

"""
    ODESolvers.dostep!(Q, lsrk::LowStorageRungeKutta2N, p, timeend::Real,
                       adjustfinalstep::Bool)

Use the 2N low storage Runge--Kutta method `lsrk` to step `Q` forward in time
from the current time, to the time `timeend`. If `adjustfinalstep == true` then
`dt` is adjusted so that the step does not take the solution beyond the
`timeend`.
"""
function ODEs.dostep!(Q, lsrk::LowStorageRungeKutta2N, p, timeend::Real,
                      adjustfinalstep::Bool)
  time, dt = lsrk.t, lsrk.dt
  if adjustfinalstep && time + dt > timeend
    dt = timeend - time
  end
  @assert dt > 0

  ODEs.dostep!(Q, lsrk, p, time, dt)

  if dt == lsrk.dt
    lsrk.t += dt
  else
    lsrk.t = timeend
  end

end

"""
    ODESolvers.dostep!(Q, lsrk::LowStorageRungeKutta2N, p, time::Real,
                       dt::Real, [slow_δ, slow_rv_dQ, slow_scaling])

Use the 2N low storage Runge--Kutta method `lsrk` to step `Q` forward in time
from the current time `time` to final time `time + dt`.

If the optional parameter `slow_δ !== nothing` then `slow_rv_dQ * slow_δ` is
added as an additionall ODE right-hand side source. If the optional parameter
`slow_scaling !== nothing` then after the final stage update the scaling
`slow_rv_dQ *= slow_scaling` is performed.
"""
function ODEs.dostep!(Q, lsrk::LowStorageRungeKutta2N, p, time::Real,
                      dt::Real, slow_δ = nothing, slow_rv_dQ = nothing,
                      in_slow_scaling = nothing)
  RKA, RKB, RKC = lsrk.RKA, lsrk.RKB, lsrk.RKC
  rhs!, dQ = lsrk.rhs!, lsrk.dQ

  rv_Q = realview(Q)
  rv_dQ = realview(dQ)

  threads = 256
  blocks = div(length(rv_Q) + threads - 1, threads)

  for s = 1:length(RKA)
    rhs!(dQ, Q, p, time + RKC[s] * dt, increment = true)

    slow_scaling = nothing
    if s == length(RKA)
      slow_scaling = in_slow_scaling
    end
    # update solution and scale RHS
    @launch(device(Q), threads=threads, blocks=blocks,
            update!(rv_dQ, rv_Q, RKA[s%length(RKA)+1], RKB[s], dt,
                    slow_δ, slow_rv_dQ, slow_scaling))
  end
end

"""
    LSRKEulerMethod(f, Q; dt, t0 = 0)

This function returns a [`LowStorageRungeKutta2N`](@ref) time stepping object
for explicitly time stepping the differential
equation given by the right-hand-side function `f` with the state `Q`, i.e.,

```math
  \\dot{Q} = f(Q, t)
```

with the required time step size `dt` and optional initial time `t0`.  This
time stepping object is intended to be passed to the `solve!` command.

This method uses the LSRK2N framework to implement a simple Eulerian forward time stepping scheme for the use of debugging.

### References

"""
function LSRKEulerMethod(F, Q::AT; dt=nothing, t0=0) where {AT <: AbstractArray}
  T = eltype(Q)
  RT = real(T)

  RKA = (RT(0),)

  RKB = (RT(1),)

  RKC = (RT(0),)

  LowStorageRungeKutta2N(F, RKA, RKB, RKC, Q; dt=dt, t0=t0)
end

"""
    LSRK54CarpenterKennedy(f, Q; dt, t0 = 0)

This function returns a [`LowStorageRungeKutta2N`](@ref) time stepping object
for explicitly time stepping the differential
equation given by the right-hand-side function `f` with the state `Q`, i.e.,

```math
  \\dot{Q} = f(Q, t)
```

with the required time step size `dt` and optional initial time `t0`.  This
time stepping object is intended to be passed to the `solve!` command.

This uses the fourth-order, low-storage, Runge--Kutta scheme of Carpenter
and Kennedy (1994) (in their notation (5,4) 2N-Storage RK scheme).

### References

    @TECHREPORT{CarpenterKennedy1994,
      author = {M.~H. Carpenter and C.~A. Kennedy},
      title = {Fourth-order {2N-storage} {Runge-Kutta} schemes},
      institution = {National Aeronautics and Space Administration},
      year = {1994},
      number = {NASA TM-109112},
      address = {Langley Research Center, Hampton, VA},
    }
"""
function LSRK54CarpenterKennedy(F, Q::AT; dt=0, t0=0) where {AT <: AbstractArray}
  T = eltype(Q)
  RT = real(T)

  RKA = (RT(0),
         RT(-567301805773  // 1357537059087),
         RT(-2404267990393 // 2016746695238),
         RT(-3550918686646 // 2091501179385),
         RT(-1275806237668 // 842570457699 ))

  RKB = (RT(1432997174477 // 9575080441755 ),
         RT(5161836677717 // 13612068292357),
         RT(1720146321549 // 2090206949498 ),
         RT(3134564353537 // 4481467310338 ),
         RT(2277821191437 // 14882151754819))

  RKC = (RT(0),
         RT(1432997174477 // 9575080441755),
         RT(2526269341429 // 6820363962896),
         RT(2006345519317 // 3224310063776),
         RT(2802321613138 // 2924317926251))

  LowStorageRungeKutta2N(F, RKA, RKB, RKC, Q; dt=dt, t0=t0)
end

"""
    LSRK144NiegemannDiehlBusch((f, Q; dt, t0 = 0)

This function returns a [`LowStorageRungeKutta2N`](@ref) time stepping object
for explicitly time stepping the differential
equation given by the right-hand-side function `f` with the state `Q`, i.e.,

```math
  \\dot{Q} = f(Q, t)
```

with the required time step size `dt` and optional initial time `t0`.  This
time stepping object is intended to be passed to the `solve!` command.

This uses the fourth-order, 14-stage, low-storage, Runge--Kutta scheme of
Niegemann, Diehl, and Busch (2012) with optimized stability region

### References

    @article{niegemann2012efficient,
      title={Efficient low-storage Runge--Kutta schemes with optimized stability regions},
      author={Niegemann, Jens and Diehl, Richard and Busch, Kurt},
      journal={Journal of Computational Physics},
      volume={231},
      number={2},
      pages={364--372},
      year={2012},
      publisher={Elsevier}
    }
"""
function LSRK144NiegemannDiehlBusch(F, Q::AT; dt=0,
                                    t0=0) where {AT <: AbstractArray}
  T = eltype(Q)
  RT = real(T)

  RKA = (RT(0),
         RT(-0.7188012108672410),
         RT(-0.7785331173421570),
         RT(-0.0053282796654044),
         RT(-0.8552979934029281),
         RT(-3.9564138245774565),
         RT(-1.5780575380587385),
         RT(-2.0837094552574054),
         RT(-0.7483334182761610),
         RT(-0.7032861106563359),
         RT( 0.0013917096117681),
         RT(-0.0932075369637460),
         RT(-0.9514200470875948),
         RT(-7.1151571693922548))

  RKB = (RT(0.0367762454319673),
         RT(0.3136296607553959),
         RT(0.1531848691869027),
         RT(0.0030097086818182),
         RT(0.3326293790646110),
         RT(0.2440251405350864),
         RT(0.3718879239592277),
         RT(0.6204126221582444),
         RT(0.1524043173028741),
         RT(0.0760894927419266),
         RT(0.0077604214040978),
         RT(0.0024647284755382),
         RT(0.0780348340049386),
         RT(5.5059777270269628))

  RKC = (RT(0),
         RT(0.0367762454319673),
         RT(0.1249685262725025),
         RT(0.2446177702277698),
         RT(0.2476149531070420),
         RT(0.2969311120382472),
         RT(0.3978149645802642),
         RT(0.5270854589440328),
         RT(0.6981269994175695),
         RT(0.8190890835352128),
         RT(0.8527059887098624),
         RT(0.8604711817462826),
         RT(0.8627060376969976),
         RT(0.8734213127600976))

  LowStorageRungeKutta2N(F, RKA, RKB, RKC, Q; dt=dt, t0=t0)
end

end

max_stars_repo_name: akshaysridhar/CLIMA

max_issues_repo_name: akshaysridhar/CLIMA

max_forks_repo_name: akshaysridhar/CLIMA


text:
using DiffEqFlux, GalacticOptim, OrdinaryDiffEq

function newtons_cooling(du, u, p, t)
    temp = u[1]
    k, temp_m = p
    du[1] = dT = -k*(temp-temp_m)
  end

function true_sol(du, u, p, t)
    true_p = [log(2)/8.0, 100.0]
    newtons_cooling(du, u, true_p, t)
end

function dudt_(u,p,t)
    ann(u, p).* u
end

cb = function (p,l,pred;doplot=false) #callback function to observe training
    display(l)
    # plot current prediction against data
    if doplot
      pl = scatter(t,ode_data[1,:],label="data")
      scatter!(pl,t,pred[1,:],label="prediction")
      display(plot(pl))
    end
    return false
end

u0 = Float32[200.0]
datasize = 30
tspan = (0.0f0, 1.5f0)

t = range(tspan[1], tspan[2], length=datasize)
true_prob = ODEProblem(true_sol, u0, tspan)
ode_data = Array(solve(true_prob, Tsit5(), saveat=t))

ann = FastChain(FastDense(1,8,tanh), FastDense(8,1,tanh))
pp = initial_params(ann)
prob = ODEProblem{false}(dudt_, u0, tspan, pp)

function predict_adjoint(fullp, time_batch)
    Array(solve(prob, Tsit5(), p = fullp, saveat = time_batch))
end

function loss_adjoint(fullp, batch, time_batch)
    pred = predict_adjoint(fullp,time_batch)
    sum(abs2, batch .- pred), pred
end

k = 10
train_loader = Flux.Data.DataLoader((ode_data, t), batchsize = k)

numEpochs = 300
l1 = loss_adjoint(pp, train_loader.data[1], train_loader.data[2])[1]

optfun = OptimizationFunction((θ, p, batch, time_batch) -> loss_adjoint(θ, batch, time_batch), GalacticOptim.AutoZygote())
optprob = OptimizationProblem(optfun, pp)
using IterTools: ncycle
res1 = GalacticOptim.solve(optprob, ADAM(0.05), ncycle(train_loader, numEpochs), cb = cb, maxiters = numEpochs)
@test 10res1.minimum < l1


text:
# 2bit Reader
# ===========

"""
    TwoBitReader(input::IO)

Create a data reader of the 2bit file format.

# Arguments
* `input`: data source
"""
type TwoBitReader{T<:IO} <: Bio.IO.AbstractReader
    # input stream
    input::T

    # sequence names
    names::Vector{String}

    # file offsets
    offsets::Vector{UInt32}

    # byte-swapped or not
    swapped::Bool
end

function TwoBitReader(input::IO)
    swapped, seqcount = read_2bit_header(input)
    names, offsets = read_2bit_index(input, swapped, seqcount)
    @assert seqcount == length(names) == length(offsets)
    return TwoBitReader(input, names, offsets, swapped)
end

function Bio.IO.stream(reader::TwoBitReader)
    return reader.input
end

function Base.eltype{S}(::Type{TwoBitReader{S}})
    return SeqRecord{ReferenceSequence,Vector{UnitRange{Int}}}
end

function Base.length(p::TwoBitReader)
    return length(p.names)
end

function Base.start(p::TwoBitReader)
    return 1
end

function Base.done(p::TwoBitReader, k)
    return k > endof(p)
end

function Base.next(p::TwoBitReader, k)
    return p[k], k + 1
end

# Random access
# -------------

function Base.checkbounds(p::TwoBitReader, k::Integer)
    if 1 ≤ k ≤ endof(p)
        return true
    end
    throw(BoundsError(p, k))
end

function Base.checkbounds(p::TwoBitReader, r::Range)
    if isempty(r) || (1 ≤ first(r) && last(r) ≤ endof(p))
        return true
    end
    throw(BoundsError(p, k))
end

function Base.endof(p::TwoBitReader)
    return length(p)
end

function Base.getindex(p::TwoBitReader, k::Integer)
    checkbounds(p, k)
    seek(p.input, p.offsets[k])
    seq, mblocks = read_2bit_seq(p.input, p.swapped)
    return SeqRecord{ReferenceSequence,Vector{UnitRange{Int}}}(p.names[k], seq, mblocks)
end

function Base.getindex(p::TwoBitReader, name::AbstractString)
    k = findfirst(p.names, name)
    if k == 0
        throw(KeyError(name))
    end
    return p[k]
end

function Base.getindex(p::TwoBitReader, r::Range)
    checkbounds(p, r)
    return [p[k] for k in r]
end

function Base.getindex{S<:AbstractString}(p::TwoBitReader, names::AbstractVector{S})
    return [p[name] for name in names]
end

function read_2bit_header(input)
    signature = read(input, UInt32)
    if signature == 0x1A412743
        swapped = false
    elseif signature == 0x4327411A
        swapped = true
    else
        error("invalid 2bit signature")
    end
    version  = read32(input, swapped)
    seqcount = read32(input, swapped)
    reserved = read32(input, swapped)
    @assert version == reserved == 0
    return swapped, seqcount
end

function read_2bit_index(input, swapped, seqcount)
    names = String[]
    offsets = UInt32[]
    for i in 1:seqcount
        namesize = read(input, UInt8)
        name = read(input, UInt8, namesize)
        offset = read32(input, swapped)
        push!(names, String(name))
        push!(offsets, offset)
    end
    return names, offsets
end

function read_2bit_seq(input, swap)
    # read record header
    dnasize = read32(input, swap)
    blockcount = read32(input, swap)  # N blocks
    blockstarts = read32(input, blockcount, swap)
    blocksizes = read32(input, blockcount, swap)
    mblockcount = read32(input, swap)  # masked blocks
    mblockstarts = read32(input, mblockcount, swap)
    mblocksizes = read32(input, mblockcount, swap)
    reserved = read32(input, swap)
    @assert reserved == 0

    # read packed DNAs
    data = zeros(UInt64, cld(dnasize, 32))
    i_stop = BitIndex(dnasize, 2)
    i = BitIndex(1, 2)
    while i ≤ i_stop
        x = read(input, UInt8)
        data[index(i)] |= twobit2refseq(x) << offset(i)
        i += 8
    end

    # make an N vector
    nmask = falses(Int(dnasize))
    for j in 1:blockcount
        s = blockstarts[j]
        for k in s+1:s+blocksizes[j]
            nmask[k] = true
        end
    end

    # make masked blocks
    mblocks = UnitRange{Int}[]
    for j in 1:mblockcount
        s = mblockstarts[j]
        push!(mblocks, s+1:s+mblocksizes[j])
    end

    return ReferenceSequence(data, nmask, 1:dnasize), mblocks
end

# mapping table from .2bit DNA encoding to Bio.jl DNA encoding
const twobit2refseq_table = let
    # T: 00, C: 01, A: 10, G: 11
    f(x) = x == 0b00 ? UInt64(3) :
           x == 0b01 ? UInt64(1) :
           x == 0b10 ? UInt64(0) :
           x == 0b11 ? UInt64(2) : error()
    tcag = 0b00:0b11
    tbl = UInt64[]
    for a in tcag, b in tcag, c in tcag, d in tcag
        x::UInt64 = 0
        x |= f(a) << 0
        x |= f(b) << 2
        x |= f(c) << 4
        x |= f(d) << 6
        push!(tbl, x)
    end
    tbl
end

twobit2refseq(x::UInt8) = twobit2refseq_table[Int(x)+1]

# read 32 bits and swap bytes if necessary
function read32(input, swap)
    x = read(input, UInt32)
    if swap
        x = bswap(x)
    end
    return x
end

function read32(input, n, swap)
    xs = read(input, UInt32, n)
    if swap
        map!(bswap, xs)
    end
    return xs
end


text:
using BinaryBuilder, Pkg

name = "Gotree"
version = v"0.4.2"
hash = "24bf544c9bde1c125968dd35a5bf63df9ad6b4ce"

# Collection of sources required to build pprof
sources = [
    GitSource("https://github.com/evolbioinfo/gotree.git", hash),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/gotree/

VERSION_PACKAGE=github.com/evolbioinfo/gotree/cmd.Version
NAME=${bindir}/gotree${exeext}

mkdir -p ${bindir}
go build -o ${NAME} -ldflags "-X ${VERSION_PACKAGE}=${version}"
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms()

# The products that we will ensure are always built
products = [
    ExecutableProduct("gotree", :gotree),
]

# Dependencies that must be installed before this package can be built
dependencies = Dependency[
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; compilers=[:c, :go], julia_compat="1.6")

max_stars_repo_name: MichelJuillard/Yggdrasil

max_issues_repo_name: MichelJuillard/Yggdrasil

max_forks_repo_name: MichelJuillard/Yggdrasil


text:
import PlutoSliderServer
import PlutoSliderServer.Pluto
import PlutoSliderServer.HTTP

using Test
using UUIDs
using Base64

function poll(query::Function, timeout::Real=Inf64, interval::Real=1 / 20)
    start = time()
    while time() < start + timeout
        if query()
            return true
        end
        sleep(interval)
    end
    return false
end

select(f::Function, xs) =
    for x in xs
        if f(x)
            return x
        end
    end

@testset "Folder watching" begin
    test_dir = tempname(cleanup=false)
    mkdir(test_dir)

    try
        # open the folder on macos:
        run(`open $(test_dir)`)
    catch
    end

    notebook_paths_to_copy = ["basic2.jl"]

    for p in notebook_paths_to_copy
        cp(joinpath(@__DIR__, p), joinpath(test_dir, p))
    end

    port = rand(12345:65000)

    still_booting = Ref(true)
    ready_result = Ref{Any}(nothing)
    function on_ready(result)
        ready_result[] = result
        still_booting[] = false
    end

    t = Pluto.@asynclog begin
        try
            PlutoSliderServer.run_directory(
                test_dir;
                Export_enabled=false,
                Export_output_dir=test_dir,
                SliderServer_port=port,
                SliderServer_watch_dir=true,
                on_ready,
            )
        catch e
            if !(e isa TaskFailedException)
                showerror(stderr, e, stacktrace(catch_backtrace()))
            end
        end
    end

    while still_booting[]
        sleep(0.1)
    end

    notebook_sessions = ready_result[].notebook_sessions

    @testset "Adding a file" begin

        cp(joinpath(test_dir, "basic2.jl"), joinpath(test_dir, "basic2 copy.jl"))

        @test poll(10, 1 / 20) do
            length(notebook_sessions) == length(notebook_paths_to_copy) + 1
        end

        newsesh = () -> select(s -> s.path == "basic2 copy.jl", notebook_sessions)

        @test !isnothing(newsesh())
        @test newsesh().current_hash != newsesh().desired_hash

        @test poll(60, 1 / 20) do
            newsesh().current_hash == newsesh().desired_hash
        end

        @test isfile(joinpath(test_dir, "basic2 copy.html"))

        @test !occursin(
            "slider_server_url = undefined",
            read(joinpath(test_dir, "basic2 copy.html"), String),
        )
        @test occursin(
            "slider_server_url = \".\"",
            read(joinpath(test_dir, "basic2 copy.html"), String),
        )
    end

    @testset "Removing a file" begin

        rm(joinpath(test_dir, "basic2 copy.jl"))

        @test !isfile(joinpath(test_dir, "basic2 copy.jl"))

        @test poll(30, 1 / 20) do
            length(notebook_sessions) == length(notebook_paths_to_copy)
        end
        @test !isfile(joinpath(test_dir, "basic2 copy.html"))

    end

    coolsesh = () -> select(s -> s.path == "subdir/cool.jl", notebook_sessions)
    coolcontents() = read(joinpath(test_dir, "subdir", "cool.html"), String)

    @testset "Adding a file (again)" begin

        mkdir(joinpath(test_dir, "subdir"))

        cp(joinpath(test_dir, "basic2.jl"), joinpath(test_dir, "subdir", "cool.jl"))

        @test poll(60, 1 / 20) do
            isfile(joinpath(test_dir, "subdir", "cool.html"))
        end
        @test poll(5, 1 / 20) do
            coolsesh().current_hash == coolsesh().desired_hash
        end

        @test !occursin("slider_server_url = undefined", coolcontents())
        @test occursin("slider_server_url = \"..\"", coolcontents())
    end

    @testset "Update an existing file" begin

        function coolconnectionkeys()
            response = HTTP.get(
                "http://localhost:$(port)/bondconnections/$(HTTP.URIs.escapeuri(coolsesh().current_hash))/",
            )
            result = Pluto.unpack(response.body)
            keys(result) |> collect |> sort
        end

        @test coolconnectionkeys() == sort(["x", "y", "s", "s2"])

        old_html_contents = coolcontents()

        Pluto.readwrite(
            joinpath(@__DIR__, "parallelpaths4.jl"),
            joinpath(test_dir, "subdir", "cool.jl"),
        )

        @test poll(5, 1 / 20) do
            coolsesh().current_hash != coolsesh().desired_hash
        end
        @test isfile(joinpath(test_dir, "subdir", "cool.html"))
        @test poll(60, 1 / 20) do
            coolsesh().current_hash == coolsesh().desired_hash
        end
        @test isfile(joinpath(test_dir, "subdir", "cool.html"))
        @test coolcontents() != old_html_contents

        @test coolconnectionkeys() == sort([
            "x",
            "y",
            "show_dogs",
            "b",
            "c",
            "five1",
            "five2",
            "six1",
            "six2",
            "six3",
            "cool1",
            "cool2",
            "world",
            "boring",
        ])
    end

    sleep(2)
    close(ready_result[].serversocket)

    try
        wait(t)
    catch e
        if !(e isa TaskFailedException)
            rethrow(e)
        end
    end
    # schedule(t, InterruptException(), error=true)
    @info "DONEZO"
end


text:
###################################################
## networktypes.jl
## Types representing network timings problem (NetworkData) and network timings estimates (Network Timings)
###################################################

"""
    Represents the data about the trip orig=>dest in the network
    Origin and destination are represented as (edge, weight) pairs for generality.
    If origin and destination are nodes, then the edge is just a loop from that node to that
    node, with weight 0. So a path from node 1 to node 5 will be indicated as
     (1,1,0.0) to (5,5,0.0)
"""
struct NetworkTrip
    "Origin edge and weight of that edge"
    orig::Tuple{Int,Int,Float64}
    "Destination edge and weight of that edge"
    dest::Tuple{Int,Int,Float64}
    "Travel time"
    time::Float64
    "Number of rides"
    weight::Float64
    "Boolean that identifies if trip results from road projection"
    roadProj::Bool
end
function Base.show(io::IO, t::NetworkTrip)
    if t.roadProj
        @printf(io,"(%d,%d)=>(%d,%d), w=%.1f, t=%.1fs",t.orig[1],t.orig[2],t.dest[1],t.dest[2],t.weight,t.time)
    else
        @printf(io,"%d=>%d, w=%.2f, t=%.1fs",t.orig[1],t.dest[1],t.weight,t.time)
    end
end

"""
    `NetworkData`
    Represent all data to use for timings inference in Network
"""
mutable struct NetworkData
    "The city routing network"
    network::Network
    "The trips data"
    trips::Vector{NetworkTrip}
    "lower bounds on roads timing (maximum speed)"
    minTimes::AbstractArray{Float64,2}
end

NetworkData(n::Network, trips::Vector{NetworkTrip}) = NetworkData(n, trips, maxSpeedTimes(n))

function Base.show(io::IO, n::NetworkData)
    density = 100*length(n.trips)/((nv(n.network.graph)-1)*nv(n.network.graph))
    println(io,"NetworkData: trip information in network")
    @printf(io,"%d trips (%.2f%% density)",length(n.trips), density)
end

"""
    Contains timings from anywhere to anywhere in a network
    - represents a solution to a NetworkData problem
    - for now, just corresponds to RoutingPaths object
    useful attributes:
    - `times` : link times
    - `pathTimes`: all timings
"""
const NetworkTimings = RoutingNetworks.RoutingPaths


text:
using Test
using Cerberus
using SparseArrays
import DisjunctiveConstraints, Gurobi, MathOptInterface
const MOI = MathOptInterface
const MOIU = MOI.Utilities

include("util.jl")
include("algorithm/util.jl")

const CONFIG = Cerberus.AlgorithmConfig(silent = true)
const GRB_ENV = isdefined(Main, :GRB_ENV) ? Main.GRB_ENV : Gurobi.Env()

for (root, dirs, files) in walkdir(@__DIR__)
    for _file in filter(f -> endswith(f, ".jl"), files)
        file = relpath(joinpath(root, _file), @__DIR__)
        if file in ["runtests.jl", "util.jl"]
            continue
        end

        @testset "$(file)" begin
            include(file)
        end
    end
end


text:
@doc raw"""
    JONSWAP{K}(α,ωₚ,γ,r)
    JONSWAP{K}(x)

4 parameter JONSWAP model for vertical displacement.

# Type parameter
The type parameter `K` is a non-negative integer representing the ammount of aliasing to be done.
This is useful as some devices have high-pass filters which limit the ammount of aliasing that will be present in a recorded series.
In particular, a `JONSWAP{K}` model for a series recorded every `Δ` seconds would be appropriate for a series sampled every `Δ` seconds and filtered with a high-pass filter at `(2K-1)π/Δ`.

# Arguments
- `α`: The scale parameter.
- `ωₚ`: The peak angular frequency.
- `γ`: The peak enhancement factor.
- `r`: The tail decay index.

# Vector constructor arguments
- `x`: Vector of parameters `[α,ωₚ,γ,r]`.

# Background
The JONSWAP spectral density function is
```math
f(ω) = αω^{-r}\exp\left \{-\frac{r}{4}\left(\frac{|ω|}{ωₚ}\right)^{-4}\right \}γ^{δ(|ω|)}
```
where
```math
δ(ω) = \exp\left\{-\tfrac{1}{2 (0.07+0.02\cdot\mathbb{1}_{ωₚ>|ω|})^2}\left (\tfrac{|ω|}{ωₚ}-1\right )^2\right\}.
```

"""
struct JONSWAP{K} <: WhittleLikelihoodInference.UnknownAcvTimeSeriesModel{1,Float64}
    α::Float64
    ωₚ::Float64
    γ::Float64
    r::Float64

    r_over4::Float64
    ωₚ²::Float64
    ωₚ³::Float64
    ωₚ⁴::Float64
    ωₚ⁶::Float64
    logγ::Float64
    function JONSWAP{K}(α,ωₚ,γ,r) where {K}
        α > 0 || throw(ArgumentError("JONSWAP requires α > 0"))
        ωₚ > 0 || throw(ArgumentError("JONSWAP requires ωₚ > 0"))
        γ >= 1 || throw(ArgumentError("JONSWAP requires γ > 1"))
        r > 1 || throw(ArgumentError("JONSWAP requires r > 1"))
        new{K}(α,ωₚ,γ,r,r/4,ωₚ^2,ωₚ^3,ωₚ^4,ωₚ^6,log(γ))
    end
    function JONSWAP{K}(x::AbstractVector{Float64}) where {K}
        @boundscheck checkparameterlength(x,JONSWAP{K})
        @inbounds JONSWAP{K}(x[1], x[2], x[3], x[4])
    end
end

# functions to throw informative error if type parameter not provided
JONSWAP(x::AbstractVector{Float64}) = JONSWAP(1,1,1,1)
JONSWAP(α,ωₚ,γ,r) = error("JONSWAP process requires the ammount of aliasing specified as a type parameter. Use JONSWAP{K}() where K ∈ N.")

WhittleLikelihoodInference.npars(::Type{JONSWAP{K}}) where {K} = 4
WhittleLikelihoodInference.nalias(::JONSWAP{K}) where {K} = K

WhittleLikelihoodInference.lowerbounds(::Type{JONSWAP{K}}) where {K} = [0,0,1,1]
WhittleLikelihoodInference.upperbounds(::Type{JONSWAP{K}}) where {K} = [Inf,Inf,Inf,Inf]

@inline @fastmath function WhittleLikelihoodInference.sdf(model::JONSWAP{K}, ω::Real) where {K}
    α,ωₚ,γ,r = model.α,model.ωₚ,model.γ,model.r
    ω = abs(ω)
    if ω > 1e-10
        σ1² = 0.0049 + 0.0032 * (ω > ωₚ)
        δ = exp(-1 / (2σ1²) * (ω / ωₚ - 1)^2)
        ωₚ⁴_over_ω⁴ = model.ωₚ⁴ / (ω^4)
        return (α * ω^(-r) * exp(-(model.r_over4) * ωₚ⁴_over_ω⁴) * γ^δ) / 2
    else
        return 0.0
    end
end

@propagate_inbounds @fastmath function WhittleLikelihoodInference.grad_add_sdf!(out, model::JONSWAP{K}, ω::Real) where {K}
    @boundscheck checkbounds(out,1:4)
    @inbounds begin
        ω = abs(ω)
        if ω > 1e-10
            α,ωₚ,γ,r = model.α,model.ωₚ,model.γ,model.r
            σ1² = 0.0049 + 0.0032 * (ω > ωₚ)
            δ = exp(-1 / (2σ1²) * (ω / ωₚ - 1)^2)
            ω⁻⁴ = ω^(-4)
            ωₚ⁴_over_ω⁴ = model.ωₚ⁴ * ω⁻⁴
            sdf = (α * ω^(-r) * exp(-(model.r_over4) * ωₚ⁴_over_ω⁴) * γ^δ) / 2

            ∂S∂α = sdf / α
            ∂S∂ωₚ = sdf * (δ*model.logγ * ω / σ1² *(ω-ωₚ) / model.ωₚ³ - r*model.ωₚ³*ω⁻⁴)
            ∂S∂γ = sdf * δ / γ
            ∂S∂r = sdf * (-log(ω)-ωₚ⁴_over_ω⁴/4)

            out[1] += ∂S∂α
            out[2] += ∂S∂ωₚ
            out[3] += ∂S∂γ
            out[4] += ∂S∂r
        end # 0 otherwise
    end
    return nothing
end

@propagate_inbounds @fastmath function WhittleLikelihoodInference.hess_add_sdf!(out, model::JONSWAP{K}, ω::Real) where {K}
    @boundscheck checkbounds(out,1:10)
    @inbounds begin
        ω = abs(ω)
        if ω > 1e-10
            α,ωₚ,γ,r = model.α,model.ωₚ,model.γ,model.r
            σ1² = 0.0049 + 0.0032 * (ω > ωₚ)
            δ = exp(-1 / (2σ1²) * (ω / ωₚ - 1)^2)
            ω⁻⁴ = ω^(-4)
            ωₚ⁴_over_ω⁴ = model.ωₚ⁴ * ω⁻⁴
            ωₚ³_over_ω⁴ = model.ωₚ³ * ω⁻⁴
            sdf = (α * ω^(-r) * exp(-(model.r_over4) * ωₚ⁴_over_ω⁴) * γ^δ) / 2

            δωlogγ_over_σ1² = δ*model.logγ * ω / σ1²
            ∂S∂ωₚUsepart = (δ*model.logγ * ω / σ1² *(ω-ωₚ) / model.ωₚ³ - r*ωₚ³_over_ω⁴)
            ∂S∂ωₚ = sdf * ∂S∂ωₚUsepart
            δ_over_γ = δ / γ
            ∂S∂γ = sdf * δ_over_γ
            ∂r_part = (-log(ω)-ωₚ⁴_over_ω⁴/4)
            ∂S∂r = sdf * ∂r_part

            ∂S∂αωₚ = ∂S∂ωₚ / α
            ∂S∂αγ = ∂S∂γ / α
            ∂S∂αr = ∂S∂r / α

            ∂S∂ωₚ2 = ∂S∂ωₚ * ∂S∂ωₚUsepart + sdf * (δωlogγ_over_σ1² *((-3ω + 2ωₚ)/model.ωₚ⁴ + ω * (ω-ωₚ)^2/model.ωₚ⁶/σ1²) - 3r*model.ωₚ²*ω⁻⁴)
            ∂S∂ωₚγ = ∂S∂γ * ∂S∂ωₚUsepart + sdf * δ_over_γ * ω / σ1²*(ω-ωₚ)/model.ωₚ³
            ∂S∂ωₚr = ∂S∂r * ∂S∂ωₚUsepart - sdf * model.ωₚ³*ω⁻⁴
            ∂S∂γ2 = δ_over_γ * (∂S∂γ - sdf/γ)
            ∂S∂γr = ∂S∂r * δ_over_γ
            ∂S∂r2 = ∂S∂r * ∂r_part

            # out[1] += 0
            out[2] += ∂S∂αωₚ
            out[3] += ∂S∂αγ
            out[4] += ∂S∂αr
            out[5] += ∂S∂ωₚ2
            out[6] += ∂S∂ωₚγ
            out[7] += ∂S∂ωₚr
            out[8] += ∂S∂γ2
            out[9] += ∂S∂γr
            out[10] += ∂S∂r2
        end # 0 otherwise
    end

    return nothing
end