text:
"""
    BreakingNLI

A dataset of 8193 premise-hypothesis sentence-pairs for NLI.

Each example is annotated to entailment, contradiction, and neutral.
The premise and the hypothesis are identical except for one
word/phrase that was replaced. This dataset is meant for testing
methods trained to solve the natural language inference task, and it
requires some lexical and world knowledge to achieve reasonable
performance on it.

For details, see the [GitHub page](https://github.com/BIU-NLP/Breaking_NLI) or read the [2018 paper](https://aclweb.org/anthology/P18-2103/).

Available data:
```julia
BreakingNLI.test_jsonl()
```
"""
module BreakingNLI

import ...register_data
using DataDeps: @datadep_str, unpack

test_jsonl() = joinpath(datadep"Breaking_NLI", "data", "dataset.jsonl")

function __init__()
    register_data(
        "Breaking_NLI",
        """
        Each example is annotated to entailment, contradiction, and neutral.
        The premise and the hypothesis are identical except for one
        word/phrase that was replaced. This dataset is meant for testing
        methods trained to solve the natural language inference task, and it
        requires some lexical and world knowledge to achieve reasonable
        performance on it.

        For details, see the [GitHub page](https://github.com/BIU-NLP/Breaking_NLI) or read the [2018 paper](https://aclweb.org/anthology/P18-2103/).

        This dataset is distributed as a

        Licensed under a Creative Commons Attribution-ShareAlike 4.0
        International License:
        https://creativecommons.org/licenses/by-sa/4.0/
        """,
        "https://github.com/BIU-NLP/Breaking_NLI/archive/8a7658c1ce6b732f4e8af3b06560f1a13b8b18b0.zip",
        "0ccf9e308245036dec52e369316b94326ef06c96729b336d0269974b0814581e",
        postfetch = function (zip)
            unpack(zip)
            unpack(joinpath(datadep"Breaking_NLI",
                            "Breaking_NLI-8a7658c1ce6b732f4e8af3b06560f1a13b8b18b0",
                            "breaking_nli_dataset.zip"))
        end
    )
end

end # module


text:
# Create AtomicSolution for ODE.
function AtomicSolution(equation::AbstractEquationODE)
    AtomicSolutionODE(equation.t₀, equation.q₀[begin])
end

# Create AtomicSolution for ODE.
function AtomicSolution(solution::SolutionODE)
    AtomicSolutionODE(get_initial_conditions(solution, 1)...)
end

# Create AtomicSolution for partitioned ODE.
function AtomicSolution(equation::AbstractEquationPODE)
    AtomicSolutionPODE(equation.t₀, equation.q₀[begin], equation.p₀[begin])
end

# Create AtomicSolution for partitioned ODE.
function AtomicSolution(solution::SolutionPODE)
    AtomicSolutionPODE(get_initial_conditions(solution, 1)...)
end

# Create AtomicSolution for DAE.
function AtomicSolution(equation::AbstractEquationDAE)
    AtomicSolutionDAE(equation.t₀, equation.q₀[begin], equation.λ₀[begin])
end

# Create AtomicSolution for DAE.
function AtomicSolution(solution::SolutionDAE)
    AtomicSolutionDAE(get_initial_conditions(solution, 1)...)
end

# Create AtomicSolution for partitioned DAE.
function AtomicSolution(equation::AbstractEquationPDAE)
    AtomicSolutionPDAE(equation.t₀, equation.q₀[begin], equation.p₀[begin], equation.λ₀[begin])
end

# Create AtomicSolution for partitioned DAE.
function AtomicSolution(solution::SolutionPDAE)
    AtomicSolutionPDAE(get_initial_conditions(solution, 1)...)
end

# Print error for AtomicSolutions of equations not implemented, yet.
function AtomicSolution(equation::Equation)
    error("No AtomicSolution found for equation ", equation)
end

# Print error for AtomicSolutions of solution not implemented, yet.
function AtomicSolution(solution::Solution)
    error("No AtomicSolution found for solution ", solution)
end

ext: jl


text:
#SignalProcessing:
#-------------------------------------------------------------------------------

module SignalProcessing

using MDDatasets

#Functions to this module will be accessing/appending:
import MDDatasets: broadcastMD
import MDDatasets: value
import MDDatasets: DataF1
import MDDatasets: DS #To dispatch on a symbol & minimize namespace pollution
import MDDatasets: CastType
import FFTW: rfft, irfft

include("vectorop.jl")
include("base.jl")
include("timefreq.jl")
include("datasetop.jl")
include("dt_generators.jl")
include("ct_generators.jl")
include("show.jl")

#==Warnings
================================================================================
-To save time, many of the functions do not currently maintain numeric data
 types. For example, feeding in an array of Float32 might result in returing an
 array of Float64.  This should be corrected.

-Pole(:SYMBOL, 5) can construct objects with any invalid symbol (Neither :Hz or
 :rad).  This will cause akward error messages.
 TODO: Is there a better way?
==#

#==Exported symbols
===============================================================================#
#==Objects to help with multiple dispatch.
NOTE:
   In turn, the function names become more succinct, and fewer functions need to
   be exported (reduces namespace pollution).
==#
export Pole

#Time/Frequency domain signals:
export DataTime, DataFreq, DataZ

#Functions
#-------------------------------------------------------------------------------
export timedomain, freqdomain #Performs (i)fft, if necessary
export fspectrum #Returns sampled frequency spectrum (non-periodic signals)
export fcoeff #Returns Fourier-series coefficients (t-periodic signals)
export pulse #Gererate pulse response
export prbs #Generate PRBS sequence
export pattern #Generate pattern from a bit sequence
export datavec #getter function
export timespace #Generate a range representing time
#==Implements timespace(:PRIMARY, v1, :SECONDARY, v2; tstart = tstart)
	where PRIMARY/SECONDARY can be of the following combinations:
	(PRIMARY is targeted value - secondary might get rounded)
		(:ts, :tfund)
		(:tfund, :ts)
		(:ts, :n)        *TODO
		(:tfund, :n)     *TODO
		(:tsig, :n)      *TODO
	AND:
		ts: sampling period
		tfund: fundamental period of frequency-domain signal
		tspan: actual signal time span (instead of its *fundamental* period)
		n: number of sampling points (always exact)
==#

#==Other interface tools (symbols not exported to avoid collisions):
================================================================================
	MDDatasets.value... already part of MDDatasets
	Base.step... already part of base, so can't export it...
==#

end #SignalProcessing

#Last line


text:
using Test
using SwapStreams
using Documenter

@test_throws TypeError SwapStream{1}(IOBuffer())

strue = SwapStream(IOBuffer())
sfalse = SwapStream{false}(IOBuffer())

@test !isreadonly(strue)
@test !isreadonly(sfalse)
@test iswritable(strue)
@test isopen(strue)
@test isopen(sfalse)
@test isreadable(strue)
@test isreadable(sfalse)
@test position(strue) == 0
@test position(sfalse) == 0
write(strue, 1)
write(sfalse, 1)
@test position(strue) == 8
@test position(sfalse) == 8
mark(strue)
mark(sfalse)
@test ismarked(strue)
@test ismarked(sfalse)
write(strue, [2, 3])
write(sfalse, [2, 3])
@test position(strue) == 24
@test position(sfalse) == 24
reset(strue)
reset(sfalse)
@test position(strue) == 8
@test position(sfalse) == 8
@test !ismarked(strue)
@test !ismarked(sfalse)

mark(strue)
@test ismarked(strue)
unmark(strue)
@test !ismarked(strue)

seekend(strue)
seekend(sfalse)
@test eof(strue)
@test eof(sfalse)
seek(strue, 0)
seek(sfalse, 0)
@test read(strue, Int) == 1
@test read(sfalse, Int) == 1

close(strue)
@test !iswritable(strue)

@test SwapStreams.is_swapping(SwapStream(ifelse(ENDIAN_BOM == BigEndian, LittleEndian, BigEndian), IOBuffer()))
@test SwapStreams.is_swapping(typeof(SwapStream(ifelse(ENDIAN_BOM == BigEndian, LittleEndian, BigEndian), IOBuffer())))
@test !SwapStreams.is_swapping(SwapStream(ifelse(ENDIAN_BOM == BigEndian, BigEndian, LittleEndian), IOBuffer()))

x = UInt8[1,2,3]
s = SwapStream(IOBuffer())
write(s, x)
seek(s, 0)
@test read(s, UInt8) === UInt8(1)
@test read!(s, Vector{UInt8}(undef, 2)) == x[2:3]

x = Int8[1,2,3]
s = SwapStream(IOBuffer())
write(s, x)
seek(s, 0)
@test read(s, Int8) === Int8(1)
@test read!(s, Vector{Int8}(undef, 2)) == x[2:3]

seek(s, 0)
@test read(s, 1) == [Int8(1)]
skip(s, 1)
@test read(s, 1) == [Int8(3)]

###
### does it work with abstract arrays
###
@testset "AbstractArray I/O" begin

    struct IntVector <: AbstractVector{Int}
        data::Vector{Int}
    end

    Base.length(x::IntVector) = length(x.data)
    Base.size(x::IntVector) = size(x.data)
    Base.getindex(x::IntVector, i::Int) = x.data[i]

    Base.iterate(x::IntVector) = iterate(x.data)
    Base.iterate(x::IntVector, state) = iterate(x.data, state)

    s = SwapStream(IOBuffer())
    x = IntVector([1,2,3]);
    write(s, x)
    seek(s, 0)
    @test read!(s, Vector{Int}(undef, 3)) == [1, 2, 3]

    s = SwapStream{false}(IOBuffer())
    x = IntVector([1,2,3]);
    write(s, x)
    seek(s, 0)
    @test read!(s, Vector{Int}(undef, 3)) == [1, 2, 3]
end

@testset "NTuple I/O" begin
    s = SwapStream{false}(IOBuffer())
    write(s, [1, 2, 3])
    seek(s, 0)
    r = Ref{NTuple{3,Int}}()
    read!(s, r) === (1, 2, 3)

    s = SwapStream{true}(IOBuffer())
    write(s, [1, 2, 3])
    seek(s, 0)
    r = Ref{NTuple{3,Int}}()
    read!(s, r) === (1, 2, 3)
end

@test stat(SwapStream{false}(open("runtests.jl"))) isa Base.Filesystem.StatStruct

doctest(SwapStreams)


text:
# ------------------------------------------------------------
# Limiter functions
# ------------------------------------------------------------

"""
    linear(sL, sR)

linear average of slopes
"""
linear(sL, sR) = 0.5 * (sL + sR)

"""
    vanleer(sL, sR)
    vanleer(sL, s, sR, connect = 2)

van Leer limiter
"""
vanleer(sL, sR) =
    (fortsign(1.0, sL) + fortsign(1.0, sR)) * abs(sL) * abs(sR) /
    (abs(sL) + abs(sR) + 1.e-7)

#--- triangle ---#
function vanleer(sL, s, sR, connect = 2)
    δ = begin
        if connect == 2
            (
                (fortsign(1.0, sL) + fortsign(1.0, s)) * abs(sL) * abs(s) /
                (abs(sL) + abs(s) + 1.e-7),
                (fortsign(1.0, s) + fortsign(1.0, sR)) * abs(s) * abs(sR) /
                (abs(s) + abs(sR) + 1.e-7),
            )
        elseif connect == 3
            (
                (fortsign(1.0, sL) + fortsign(1.0, s)) * abs(sL) * abs(s) /
                (abs(sL) + abs(s) + 1.e-7),
                (fortsign(1.0, s) + fortsign(1.0, sR)) * abs(s) * abs(sR) /
                (abs(s) + abs(sR) + 1.e-7),
                (fortsign(1.0, sL) + fortsign(1.0, sR)) * abs(sL) * abs(sR) /
                (abs(sL) + abs(sR) + 1.e-7),
            )
        end
    end

    id = findmin(abs.(δ))[2]

    return δ[id]
end

"""
    minmod(sL, sR)
    minmod(sL, s, sR, connect = 2)

MinMod limiter
"""
minmod(sL, sR) = 0.5 * (fortsign(1.0, sL) + fortsign(1.0, sR)) * min(abs(sR), abs(sL))

#--- triangle ---#
function minmod(sL, s, sR, connect = 2)
    δ = begin
        if connect == 2
            (
                0.5 * (fortsign(1.0, sL) + fortsign(1.0, s)) * min(abs(s), abs(sL)),
                0.5 * (fortsign(1.0, s) + fortsign(1.0, sR)) * min(abs(sR), abs(s)),
            )
        elseif connect == 3
            (
                0.5 * (fortsign(1.0, sL) + fortsign(1.0, s)) * min(abs(s), abs(sL)),
                0.5 * (fortsign(1.0, s) + fortsign(1.0, sR)) * min(abs(sR), abs(s)),
                0.5 * (fortsign(1.0, sL) + fortsign(1.0, sR)) * min(abs(sR), abs(sL)),
            )
        end
    end

    id = findmin(abs.(δ))[2]

    return δ[id]
end

"""
    superbee(sL, sR)

SuperBee limiter
"""
function superbee(sL, sR)
    if sR >= 0.5 * sL && sR <= 2.0 * sL
        return 0.5 * (fortsign(1.0, sL) + fortsign(1.0, sR)) * max(abs(sL), abs(sR))
    elseif sR < 0.5 * sL && sR > 2.0 * sL
        return (fortsign(1.0, sL) + fortsign(1.0, sR)) * min(abs(sL), abs(sR))
    else
        return 0.0
    end
end

"""
    vanalbaba(sL, sR)

van Albaba limiter
"""
vanalbaba(sL, sR) = (sL^2 * sR + sL * sR^2) / (sL^2 + sR^2 + 1.e-7)

"""
    weno5(wL2, wL1, wN, wR1, wR2)

5th-order WENO-JS interpolation
"""
function weno5(wL2, wL1, wN, wR1, wR2)
    ϵ = 1e-6

    β0 = 13.0 / 12.0 * (wN - 2.0 * wR1 + wR2)^2 + 1.0 / 4.0 * (3.0 * wN - 4.0 * wR1 + wR2)^2
    β1 = 13.0 / 12.0 * (wL1 - 2.0 * wN + wR1)^2 + 1.0 / 4.0 * (wL1 - wR1)^2
    β2 = 13.0 / 12.0 * (wL2 - 2.0 * wL1 + wN)^2 + 1.0 / 4.0 * (wL2 - 4.0 * wL1 + 3.0 * wN)^2

    # right interface
    dr0 = 0.3
    dr1 = 0.6
    dr2 = 0.1

    αr0 = dr0 / (ϵ + β0)^2
    αr1 = dr1 / (ϵ + β1)^2
    αr2 = dr2 / (ϵ + β2)^2

    ωr0 = αr0 / (αr0 + αr1 + αr2)
    ωr1 = αr1 / (αr0 + αr1 + αr2)
    ωr2 = αr2 / (αr0 + αr1 + αr2)

    qr0 = 1.0 / 3.0 * wN + 5.0 / 6.0 * wR1 - 1.0 / 6.0 * wR2
    qr1 = -1.0 / 6.0 * wL1 + 5.0 / 6.0 * wN + 1.0 / 3.0 * wR1
    qr2 = 1.0 / 3.0 * wL2 - 7.0 / 6.0 * wL1 + 11.0 / 6.0 * wN

    wR = ωr0 * qr0 + ωr1 * qr1 + ωr2 * qr2

    # left interface
    dl0 = 0.1
    dl1 = 0.6
    dl2 = 0.3

    αl0 = dl0 / (ϵ + β0)^2
    αl1 = dl1 / (ϵ + β1)^2
    αl2 = dl2 / (ϵ + β2)^2

    ωl0 = αl0 / (αl0 + αl1 + αl2)
    ωl1 = αl1 / (αl0 + αl1 + αl2)
    ωl2 = αl2 / (αl0 + αl1 + αl2)

    ql0 = 11.0 / 6.0 * wN - 7.0 / 6.0 * wR1 + 1.0 / 3.0 * wR2
    ql1 = 1.0 / 3.0 * wL1 + 5.0 / 6.0 * wN - 1.0 / 6.0 * wR1
    ql2 = -1.0 / 6.0 * wL2 + 5.0 / 6.0 * wL1 + 1.0 / 3.0 * wN

    wL = αl0 * ql0 + αl1 * ql1 + αl2 * ql2

    return wL, wR
end


text:
using DVIDArrays
using DVIDArrays.ImageTileArrays

a = ImageTileArray("130.211.134.112", 8000, "5c")

arr = a[2305:2604, 1473:1772,1281:1282]

using PyPlot
PyPlot.imshow( arr[:,:,1], cmap="gray" )


text:
"Tableau for explicit Euler method"
function getTableauExplicitEuler()
    a = zeros(Float64, 1, 1)
    b = [1.0]
    c = [0.0]
    o = 1

    TableauERK(:explicit_euler, o, a, b, c)
end

"Tableau for explicit midpoint method"
function getTableauExplicitMidpoint()
    a = [[0.0 0.0]
         [0.5 0.0]]
    b = [0.0, 1.0]
    c = [0.0, 0.5]
    o = 2

    TableauERK(:explicit_midpoint, o, a, b, c)
end

"Tableau for Heun's method"
function getTableauHeun()
    a = [[0.0 0.0]
         [1.0 0.0]]
    b = [0.5, 0.5]
    c = [0.0, 1.0]
    o = 2

    TableauERK(:heun, o, a, b, c)
end

"Tableau for Kutta's method of order three"
function getTableauKutta()
    a = [[ 0.0 0.0 0.0]
         [ 0.5 0.0 0.0]
         [-1.0 2.0 0.0]]
    b = [1/6, 4/6, 1/6]
    c = [0.0, 0.5, 1.0]
    o = 3

    TableauERK(:kutta, o, a, b, c)
end

"Tableau for explicit Runge-Kutta method of order four (1/6 rule)"
function getTableauERK4()
    a = [[0.0 0.0 0.0 0.0]
         [0.5 0.0 0.0 0.0]
         [0.0 0.5 0.0 0.0]
         [0.0 0.0 1.0 0.0]]
    b = [1/6, 1/3, 1/3, 1/6]
    c = [0.0, 0.5, 0.5, 1.0]
    o = 4

    TableauERK(:erk4, o, a, b, c)
end

"Tableau for explicit Runge-Kutta method of order four (3/8 rule)"
function getTableauERK438()
    a = [[ 0.0  0.0  0.0  0.0]
         [ 1/3  0.0  0.0  0.0]
         [-1/3  1.0  0.0  0.0]
         [ 1.0 -1.0  1.0  0.0]]
    b = [1/8, 3/8, 3/8, 1/8]
    c = [0.0, 1/3, 2/3, 1.0]
    o = 4

    TableauERK(:erk438, o, a, b, c)
end


text:
    formats = (
        "+##(#)##########",
        "+##(#)##########",
        "0##########",
        "0##########",
        "###-###-####",
        "(###)###-####",
        "1-###-###-####",
        "###.###.####",
        "###-###-####",
        "(###)###-####",
        "1-###-###-####",
        "###.###.####",
        "###-###-####x###",
        "(###)###-####x###",
        "1-###-###-####x###",
        "###.###.####x###",
        "###-###-####x####",
        "(###)###-####x####",
        "1-###-###-####x####",
        "###.###.####x####",
        "###-###-####x#####",
        "(###)###-####x#####",
        "1-###-###-####x#####",
        "###.###.####x#####"
    )

 function phone_number()

    t= formats[rand(1:length(formats))]

    for i=t
     t= replace(t,"#","$(rand(0:9))",1)
    end

    return t
end


text:
@testset "collatz" begin

  # A little puzzle inspired by the collatz conjecture, reach one by using the collatz conjecture moves (3n + 1 can be used even if the number is even)
  neighbours(n) = n % 2 == 0 ? [n ÷ 2, 3n + 1] : [3n + 1]

  # The best case is if the number is a power of 2, then the moves are just the log2, otherwise this is a very optimistic heuristic (it doesn't overestimate)
  heuristic(n, goal) = floor(Int64, log2(n))

  @testset "12" begin
    start = 12
    res = astar(neighbours, start, 1; heuristic)
    @test res.status == :success
    @test res.path == [12, 6, 3, 10, 5, 16, 8, 4, 2, 1]
    @test res.cost == 9
  end

  @testset "12 with maxcost" begin
    start = 12
    res = astar(neighbours, start, 1; heuristic, maxcost = 3)
    @test res.status == :nopath
    @test res.path == [12, 6, 3]
    @test res.cost == 2
  end

  @testset "12 with cost function" begin
    start = 12

    # dividing costs 1, multiplying costs 4
    cost(state, neighbour) = state / neighbour == 2 ? 1 : 4

    res = astar(neighbours, start, 1; heuristic, cost)
    @test res.status == :success
    @test res.path == [12, 6, 3, 10, 5, 16, 8, 4, 2, 1]
    @test res.cost == 15
  end

  @testset "maxdepth ignored" begin
    start = 12

    res = astar(neighbours, start, 1; heuristic, maxdepth = 0)
    @test res.status == :success
    @test res.path == [12, 6, 3, 10, 5, 16, 8, 4, 2, 1]
    @test res.cost == 9
  end

  @testset "all uninformed algorithms maxdepth" begin
    start = 12

    for algorithm in [depthfirst, breadthfirst, iterative_deepening]
      res = depthfirst(neighbours, start, 1; maxdepth = 0)
      @test res.status == :nopath
      @test res.path == [start]
    end
  end

  @testset "all algorithms timeout" begin
    start = 12

    for algo in [astar, depthfirst, breadthfirst, iterative_deepening]
      res = algo(neighbours, start, 1; timeout = 0.0)
      @test res.status == :timeout
      @test res.path == [start]
    end
  end

  @testset "all algorithms success" begin
    start = 12

    for algo in [astar, depthfirst, breadthfirst, iterative_deepening]
      res = algo(neighbours, start, 1)
      @test res.status == :success
      @test res.path == [12, 6, 3, 10, 5, 16, 8, 4, 2, 1]
    end
  end

  @testset "all algorithms disable closedset" begin
    start = 12

    for algo in [astar, depthfirst, breadthfirst, iterative_deepening]
      res = algo(neighbours, start, 1, enable_closedset = false)
      @test res.status == :success
      @test res.closedsetsize == 0
    end
  end
end


text:
@testset "Field arithmetic" begin
    @field(A) do
        I₀ = 1.0
        T = 2.0
        σ = 3.0
        Tmax = 3.0
    end

    @field(A2) do
        I₀ = 1.0
        T = 2.0
        σ = 3.0
        Tmax = 3.0
        ξ = 0.6
    end

    @testset "Summed fields" begin
        @field(B) do
            I₀ = 0.5
            T = 1.0
            σ = 3.0
            Tmax = 3.0
        end

        C = A + B
        @test span(C) == span(A)
        @test span(B) ⊆ span(C)
        @test steps(C, 100) == steps(A, 100/austrip(period(B)))

        @test polarization(C) == LinearPolarization()
        @test dimensions(C) == 1

        @test vector_potential(C, 0.4) == vector_potential(A, 0.4) + vector_potential(B, 0.4)
        @test field_amplitude(C, 0.4) == field_amplitude(A, 0.4) + field_amplitude(B, 0.4)

        # These functions only make sense when adding two field of the same wavelength
        for fun in [wavelength, period, frequency, wavenumber, fundamental, photon_energy]
            @test fun(A+A) == fun(A)
        end

        @test max_frequency(C) == frequency(B)

        @test continuity(C) == Inf

        @test field_amplitude(phase_shift(C, π/3), 0.4) ≈ field_amplitude(phase_shift(A, π/3), 0.4) + field_amplitude(phase_shift(B, π/3), 0.4)

        @test field_amplitude(-A, 0.4) == -field_amplitude(A, 0.4)
        @test field_amplitude(A-B, 0.4) == field_amplitude(A, 0.4) - field_amplitude(B, 0.4)

        withenv("UNITFUL_FANCY_EXPONENTS" => true) do
            @test string(C) == """
                               ┌ Linearly polarized field with
                               │   - I₀ = 1.0000e+00 au = 3.5094452e16 W cm⁻² =>
                               │     - E₀ = 1.0000e+00 au = 514.2207 GV m⁻¹
                               │     - A₀ = 0.3183 au
                               │   – a Fixed carrier @ λ = 14.5033 nm (T = 48.3777 as, ω = 3.1416 Ha = 85.4871 eV)
                               │   – and a Gaussian envelope of duration 170.8811 as (intensity FWHM; ±2.00σ)
                               │   – Uₚ = 0.0253 Ha = 689.2724 meV => α = 0.1013 Bohr = 5.3617 pm
                               ⊕
                               │ Linearly polarized field with
                               │   - I₀ = 5.0000e-01 au = 1.7547226e16 W cm⁻² =>
                               │     - E₀ = 7.0711e-01 au = 363.6089 GV m⁻¹
                               │     - A₀ = 0.1125 au
                               │   – a Fixed carrier @ λ = 7.2516 nm (T = 24.1888 as, ω = 6.2832 Ha = 170.9742 eV)
                               │   – and a Gaussian envelope of duration 170.8811 as (intensity FWHM; ±1.00σ)
                               └   – Uₚ = 0.0032 Ha = 86.1591 meV => α = 0.0179 Bohr = 947.8211 fm
                               """
        end
        withenv("UNITFUL_FANCY_EXPONENTS" => false) do
            @test string(C) == """
                               ┌ Linearly polarized field with
                               │   - I₀ = 1.0000e+00 au = 3.5094452e16 W cm^-2 =>
                               │     - E₀ = 1.0000e+00 au = 514.2207 GV m^-1
                               │     - A₀ = 0.3183 au
                               │   – a Fixed carrier @ λ = 14.5033 nm (T = 48.3777 as, ω = 3.1416 Ha = 85.4871 eV)
                               │   – and a Gaussian envelope of duration 170.8811 as (intensity FWHM; ±2.00σ)
                               │   – Uₚ = 0.0253 Ha = 689.2724 meV => α = 0.1013 Bohr = 5.3617 pm
                               ⊕
                               │ Linearly polarized field with
                               │   - I₀ = 5.0000e-01 au = 1.7547226e16 W cm^-2 =>
                               │     - E₀ = 7.0711e-01 au = 363.6089 GV m^-1
                               │     - A₀ = 0.1125 au
                               │   – a Fixed carrier @ λ = 7.2516 nm (T = 24.1888 as, ω = 6.2832 Ha = 170.9742 eV)
                               │   – and a Gaussian envelope of duration 170.8811 as (intensity FWHM; ±1.00σ)
                               └   – Uₚ = 0.0032 Ha = 86.1591 meV => α = 0.0179 Bohr = 947.8211 fm
                               """
        end
    end

    @testset "Delayed fields" begin
        B = delay(A, 0.4)
        @test field_amplitude(delay(A,-0.4), 0.0) == field_amplitude(B, 0.8) == field_amplitude(A, 0.4)
        @test field_amplitude(delay(A, π*u"rad"), 0.0) == field_amplitude(A, -1.0)

        @test delay(B) == 0.4
        @test iszero(delay(A))

        @test span(B) == -5.6..6.4

        withenv("UNITFUL_FANCY_EXPONENTS" => true) do
            @test string(B) == """
                               Linearly polarized field with
                                 - I₀ = 1.0000e+00 au = 3.5094452e16 W cm⁻² =>
                                   - E₀ = 1.0000e+00 au = 514.2207 GV m⁻¹
                                   - A₀ = 0.3183 au
                                 – a Fixed carrier @ λ = 14.5033 nm (T = 48.3777 as, ω = 3.1416 Ha = 85.4871 eV)
                                 – and a Gaussian envelope of duration 170.8811 as (intensity FWHM; ±2.00σ)
                                 – Uₚ = 0.0253 Ha = 689.2724 meV => α = 0.1013 Bohr = 5.3617 pm
                                 – delayed by 0.4000 jiffies = 9.6755 as"""
        end
        withenv("UNITFUL_FANCY_EXPONENTS" => false) do
            @test string(B) == """
                               Linearly polarized field with
                                 - I₀ = 1.0000e+00 au = 3.5094452e16 W cm^-2 =>
                                   - E₀ = 1.0000e+00 au = 514.2207 GV m^-1
                                   - A₀ = 0.3183 au
                                 – a Fixed carrier @ λ = 14.5033 nm (T = 48.3777 as, ω = 3.1416 Ha = 85.4871 eV)
                                 – and a Gaussian envelope of duration 170.8811 as (intensity FWHM; ±2.00σ)
                                 – Uₚ = 0.0253 Ha = 689.2724 meV => α = 0.1013 Bohr = 5.3617 pm
                                 – delayed by 0.4000 jiffies = 9.6755 as"""
        end
    end

    @testset "Combined" begin
        @field(X) do
            I₀ = 1.0
            ω = 1.0
            cycles = 6.0
            env = :cos²
            ξ = 1.0
        end

        @field(Y) do
            I₀ = 1.0
            ω = 2.0
            cycles = 6.0
            env = :cos²
            ξ = -1.0
        end

        Z = X + delay(Y, 3/2π)

        @test field_amplitude(Z, 4.0) ≈ field_amplitude(X, 4.0) + field_amplitude(Y, 4.0 - 3/2π)

        W = X + delay(Y, 35.0)
        @test span(W) == span(X).left..(span(Y).right+35.0)
    end

    @testset "Padded fields" begin
        @test_throws ArgumentError PaddedField(A, -3, 1)
        @test_throws ArgumentError PaddedField(A, 3, -1)

        B = PaddedField(A, 3.0u"fs", 8.0u"fs")
        B2 = PaddedField(A2, 3.0u"fs", 8.0u"fs")

        @test parent(B) == A
        @test parent(B2) == A2

        for (F1,F2) in [(A,B), (A2,B2)]
            @test field_amplitude(F2, 0.4) == field_amplitude(F1, 0.4)
            @test field_amplitude(F2, 7.0) == zero(field_amplitude(F1, 0.4))
        end

        @test all(endpoints(span(B)) .≈ (-6-austrip(3.0u"fs"), 6+austrip(8.0u"fs")))

        @test time_integral(B) == time_integral(A)

        withenv("UNITFUL_FANCY_EXPONENTS" => true) do
            @test string(B) == """
                               Padding before 124.0241 jiffies = 3.0000 fs and after 330.7310 jiffies = 8.0000 fs of
                               Linearly polarized field with
                                 - I₀ = 1.0000e+00 au = 3.5094452e16 W cm⁻² =>
                                   - E₀ = 1.0000e+00 au = 514.2207 GV m⁻¹
                                   - A₀ = 0.3183 au
                                 – a Fixed carrier @ λ = 14.5033 nm (T = 48.3777 as, ω = 3.1416 Ha = 85.4871 eV)
                                 – and a Gaussian envelope of duration 170.8811 as (intensity FWHM; ±2.00σ)
                                 – Uₚ = 0.0253 Ha = 689.2724 meV => α = 0.1013 Bohr = 5.3617 pm"""
        end
        withenv("UNITFUL_FANCY_EXPONENTS" => false) do
            @test string(B) == """
                               Padding before 124.0241 jiffies = 3.0000 fs and after 330.7310 jiffies = 8.0000 fs of
                               Linearly polarized field with
                                 - I₀ = 1.0000e+00 au = 3.5094452e16 W cm^-2 =>
                                   - E₀ = 1.0000e+00 au = 514.2207 GV m^-1
                                   - A₀ = 0.3183 au
                                 – a Fixed carrier @ λ = 14.5033 nm (T = 48.3777 as, ω = 3.1416 Ha = 85.4871 eV)
                                 – and a Gaussian envelope of duration 170.8811 as (intensity FWHM; ±2.00σ)
                                 – Uₚ = 0.0253 Ha = 689.2724 meV => α = 0.1013 Bohr = 5.3617 pm"""
        end
    end

    @testset "Windowed fields" begin
        B = WindowedField(A, -0.1u"fs", 0.05u"fs")
        B2 = WindowedField(A2, -0.1u"fs", 0.05u"fs")

        @test parent(B) == A
        @test parent(B2) == A2

        for fun in [vector_potential, field_amplitude, intensity]
            for (F1,F2) in [(A,B), (A2,B2)]
                @test field_amplitude(F2, 0.4) == field_amplitude(F1, 0.4)
                @test field_amplitude(F2, -5.0) == zero(field_amplitude(F1, 0.4))
                @test field_amplitude(F2, 3.0) == zero(field_amplitude(F1, 0.4))
            end
        end

        @test all(endpoints(span(B)) .≈ (-austrip(0.1u"fs"), austrip(0.05u"fs")))

        @test field_amplitude(phase_shift(B, 2.0), 0.4) == field_amplitude(phase_shift(A, 2.0), 0.4)
        @test field_amplitude(phase_shift(B, 2.0), 3.0) == zero(field_amplitude(phase_shift(A, 2.0), 0.4))

        withenv("UNITFUL_FANCY_EXPONENTS" => true) do
            @test string(B) == """
                               Window from -4.1341 jiffies = -100.0000 as to 2.0671 jiffies = 50.0000 as of
                               Linearly polarized field with
                                 - I₀ = 1.0000e+00 au = 3.5094452e16 W cm⁻² =>
                                   - E₀ = 1.0000e+00 au = 514.2207 GV m⁻¹
                                   - A₀ = 0.3183 au
                                 – a Fixed carrier @ λ = 14.5033 nm (T = 48.3777 as, ω = 3.1416 Ha = 85.4871 eV)
                                 – and a Gaussian envelope of duration 170.8811 as (intensity FWHM; ±2.00σ)
                                 – Uₚ = 0.0253 Ha = 689.2724 meV => α = 0.1013 Bohr = 5.3617 pm"""
        end
        withenv("UNITFUL_FANCY_EXPONENTS" => false) do
            @test string(B) == """
                               Window from -4.1341 jiffies = -100.0000 as to 2.0671 jiffies = 50.0000 as of
                               Linearly polarized field with
                                 - I₀ = 1.0000e+00 au = 3.5094452e16 W cm^-2 =>
                                   - E₀ = 1.0000e+00 au = 514.2207 GV m^-1
                                   - A₀ = 0.3183 au
                                 – a Fixed carrier @ λ = 14.5033 nm (T = 48.3777 as, ω = 3.1416 Ha = 85.4871 eV)
                                 – and a Gaussian envelope of duration 170.8811 as (intensity FWHM; ±2.00σ)
                                 – Uₚ = 0.0253 Ha = 689.2724 meV => α = 0.1013 Bohr = 5.3617 pm"""
        end
    end
end


text:
@inline function _AudzeEglaisDist(LHC)
    n,d = size(LHC)
    dist = 0.0
    dist_tmp = 0.0

    # Squared l-2 norm of distances between all (unique) points
    for i = 2:n
        for j = 1:i-1
            for k = 1:d
                @inbounds dist_tmp += (LHC[i,k]-LHC[j,k])^2
            end
            @inbounds dist += 1/dist_tmp
            dist_tmp = 0.0
        end
    end
    output = 1/dist
    return output
end

function _AudzeEglaisObjective(dim::Continuous,LHC)
    output = _AudzeEglaisDist(LHC)
    return output
end

function _AudzeEglaisObjective(dim::Categorical,LHC)
    output = _AudzeEglaisDist(LHC)
    output == Inf ? 0 : output
end

"""
    function AudzeEglaisObjective!(LHC::T) where T <: AbstractArray
Return the scalar which should be maximized when using the Audze-Eglais
distance as the objective function. Note this is the inverse of the typical
Audze-Eglais distance which normally is minimized.
"""
function AudzeEglaisObjective(LHC::T; dims::Array{V,1} =[Continuous() for i in 1:size(LHC,2)],
                                            interSampleWeight::Float64=1.0,
                                            ) where T <: AbstractArray where V <: LHCDimension

    out = 0.0

    #Compute the objective function among all points
    out += _AudzeEglaisObjective(Continuous(),LHC)*interSampleWeight

    #Compute the objective function within each categorical dimension
    categoricalDimInds = findall(x->typeof(x)==Categorical,dims)
    for i in categoricalDimInds
        for j = 1:dims[i].levels
            subLHC = @view LHC[LHC[:,i] .== j,:]
            out += _AudzeEglaisObjective(dims[i],subLHC)*dims[i].weight
        end
    end

    return out
end

# Remove depwarning in release 2.x.x
function AudzeEglaisObjective!(dist,LHC::T; dims::Array{V,1} =[Continuous() for i in 1:size(LHC,2)],
                                            interSampleWeight::Float64=1.0,
                                            ) where T <: AbstractArray where V <: LHCDimension
    @warn "AudzeEglaisObjective!(dist,LHC) is deprecated and does not differ from AudzeEglaisObjective(LHC)"
    AudzeEglaisObjective(LHC; dims = dims, interSampleWeight = interSampleWeight)
end


text:
using CSV, DataFrames, Gadfly, Compose

# import files
blocks = CSV.read("./results/blocks.csv", DataFrame)
bumps = CSV.read("./results/bumps.csv", DataFrame)
heavisine = CSV.read("./results/heavisine.csv", DataFrame)
doppler = CSV.read("./results/doppler.csv", DataFrame)
quadchirp = CSV.read("./results/quadchirp.csv", DataFrame)
mishmash = CSV.read("./results/mishmash.csv", DataFrame)

# build plots
function makeplot(df::DataFrame, main::AbstractString)
    # start theme
    fonts = Theme(
        point_size=4pt,
        background_color="white",
        key_title_font_size=12pt, key_title_color="black",
        key_label_font_size=10pt, key_label_color="black",
        major_label_font_size=12pt, major_label_color="black",
        minor_label_font_size=10pt, minor_label_color="black"
    )
    Gadfly.push_theme(fonts)

    # plots
    set_default_plot_size(10inch, 7inch)
    psnrplot = plot(
        df[df[:,:transform] .!= "None", :],
        y=:transform, x=:PSNR, color=:bestTH, ygroup=:method,
        xintercept=df[df[:, :transform] .== "None", :PSNR],
        Geom.subplot_grid(
            Geom.point, Geom.vline(color=["gray"], size=[2mm]), free_y_axis=true
        ),
        Guide.title("Peak Signal-to-Noise Ratio (PSNR)")
    )
    ssimplot = plot(
        df[df[:,:transform] .!= "None", :],
        y=:transform, x=:SSIM, color=:bestTH, ygroup=:method,
        xintercept=df[df[:, :transform] .== "None", :SSIM],
        Geom.subplot_grid(
            Geom.point, Geom.vline(color="gray", size=[2mm]), free_y_axis=true
        ),
        Guide.title("Structural Similarity Index Measure (SSIM)"))
    # stack plots
    result = gridstack([psnrplot ssimplot])
    # add title
    result = vstack(
        compose(
            context(0, 0, 1, 0.04),
            (context(0, 0, 1, 1), text(0.5, 1.0, main, hcenter, vbottom), Compose.fontsize(20pt), fill("black")),
            (context(), rectangle(), fill(colorant"white"))
        ),                                                          # title
        compose(context(0, 0, 1, 0.96), result))                    # plots
    return result
end

blocksplot = makeplot(blocks, "Blocks")
bumpsplot = makeplot(bumps, "Bumps")
heavisineplot = makeplot(heavisine, "Heavisine")
dopplerplot = makeplot(doppler, "Doppler")
quadchirpplot = makeplot(quadchirp, "Quadchirp")
mishmashplot = makeplot(mishmash, "Mishmash")

draw(SVG("./figures/blocks.svg", 10inch, 7inch), blocksplot)
draw(SVG("./figures/bumps.svg", 10inch, 7inch), bumpsplot)
draw(SVG("./figures/heavisine.svg", 10inch, 7inch), heavisineplot)
draw(SVG("./figures/doppler.svg", 10inch, 7inch), dopplerplot)
draw(SVG("./figures/quadchirp.svg", 10inch, 7inch), quadchirpplot)
draw(SVG("./figures/mishmash.svg", 10inch, 7inch), mishmashplot)

max_stars_repo_name: zengfung/WaveletsDenoisingExperiment

max_issues_repo_name: zengfung/WaveletsDenoisingExperiment

max_forks_repo_name: zengfung/WaveletsDenoisingExperiment


text:
function tss(y::AbstractVector, hasintercept::Bool, weights::AbstractWeights)
    if hasintercept
        m = mean(y, weights)
        sum(@inbounds (y[i] - m)^2 * weights[i] for i in eachindex(y))
    else
        sum(@inbounds y[i]^2 * weights[i] for i in eachindex(y))
    end
end

function Fstat(coef::Vector{Float64}, matrix_vcov::AbstractMatrix{Float64}, has_intercept::Bool)
    coefF = copy(coef)
    # TODO: check I can't do better
    length(coef) == has_intercept && return NaN
    if has_intercept
        coefF = coefF[2:end]
        matrix_vcov = matrix_vcov[2:end, 2:end]
    end
    return (coefF' * (matrix_vcov \ coefF)) / length(coefF)
end


text:
#######################################################################
# setting up
#######################################################################
using Revise
using FastGroupBy, BenchmarkTools, SortingLab, CategoricalArrays, Test
using Random

# import Base: getindex, similar, setindex!, size
N = 100_000; K = 100
Random.seed!(1);
val = rand(1:5, N);
pool = "id".*lpad.(1:100,3, '0');
fn = sum;

#######################################################################
# convert to CategoricalVector
#######################################################################
y = CategoricalArray{String, 1}(rand(UInt32(1):UInt32(length(pool)), N), CategoricalPool(pool, true));
y = compress(y);
# @benchmark sort($y)

z = CategoricalArray{String, 1}(rand(UInt32(1):UInt32(length(pool)), N), CategoricalPool(pool, true));
z = compress(z);
byveccv = (y, z);
 # 2mins for 2b length vectors

@benchmark fgroupreduce($+, $byveccv, $val) # 7.5 seconds for 2billion
res1 = fgroupreduce(+, byveccv, val)
res1max = fgroupreduce(max, byveccv, val)

@benchmark fby($sum, $byveccv, $val)
res = fby(sum, byveccv, val)
# @code_warntype fby(sum, byveccv, val)

@test all(res[1] .== res1[1])
@test all(res[2] .== res1[2])
@test all(res[3] .== res1[3])

sort!(dfa, cols=[:y, :z])

@test all(dfa[:y] .== res1[1])
@test all(dfa[:z] .== res1[2])
@test all(dfa[:val_sum] .== res1[3])

import FastGroupBy: fgroupreduce
fgroupreduce(fn, df, bysyms::Tuple{Symbol, Symbol}, val::Symbol) = DataFrame([fgroupreduce(fn, (df[bysyms[1]], df[bysyms[2]]), df[val])...], [bysyms..., val])
@time a = fgroupreduce(+, df, (:y, :z), :val)

DataFrame([a...])

if false
    using DataFrames

    df = DataFrame(y = y, z = z, val = val)

    @benchmark aggregate($df, $[:y,:z], $sum)
    # BenchmarkTools.Trial:
    #   memory estimate:  6.28 GiB
    #   allocs estimate:  101028306
    #   --------------
    #   minimum time:     20.904 s (4.59% GC)
    #   median time:      20.904 s (4.59% GC)
    #   mean time:        20.904 s (4.59% GC)
    #   maximum time:     20.904 s (4.59% GC)
    #   --------------
    #   samples:          1
    #   evals/sample:     1

    @benchmark by($df, $[:y,:z], df1->sum(df1[:val]))
    # BenchmarkTools.Trial:
    # memory estimate:  6.25 GiB
    # allocs estimate:  100628265
    # --------------
    # minimum time:     29.894 s (3.82% GC)
    # median time:      29.894 s (3.82% GC)
    # mean time:        29.894 s (3.82% GC)
    # maximum time:     29.894 s (3.82% GC)
    # --------------
    # samples:          1
    # evals/sample:     1

    @benchmark fgroupreduce($+, $byveccv, $val) samples = 5 seconds = 20 # 7.5 seconds for 2billion
    # BenchmarkTools.Trial:
    # memory estimate:  311.27 KiB
    # allocs estimate:  1338
    # --------------
    # minimum time:     366.854 ms (0.00% GC)
    # median time:      387.955 ms (0.00% GC)
    # mean time:        401.369 ms (0.00% GC)
    # maximum time:     443.034 ms (0.00% GC)
    # --------------
    # samples:          13
    # evals/sample:     1

    using IndexedTables, IterableTables
    t = table(df)
    ti = reindex(t,(:y,:z))

    @benchmark IndexedTables.groupreduce(+, ti, (:y,:z), select = :val) samples = 5 seconds = 20
    # BenchmarkTools.Trial:
    # memory estimate:  1.07 MiB
    # allocs estimate:  10243
    # --------------
    # minimum time:     2.465 s (0.00% GC)
    # median time:      2.561 s (0.00% GC)
    # mean time:        2.555 s (0.00% GC)
    # maximum time:     2.660 s (0.00% GC)
    # --------------
    # samples:          5
    # evals/sample:     1

    using RCall

    r_timing = R"""
    memory.limit(2^31-1)
    N=1e8; K=100
    set.seed(1)
    library(data.table)
    DT <- data.table(
    id1 = sample(sprintf("id%03d",1:K), N, TRUE),      # large groups (char)
    id2 = sample(sprintf("id%03d",1:K), N, TRUE),      # large groups (char)
    v1 =  sample(5, N, TRUE)                          # int in range [1,5]
    )

    system.time( DT[, sum(v1), keyby="id1,id2"] )
    """
    # RCall.RObject{RCall.RealSxp}
    # user  system elapsed
    # 4.75    0.29    5.12

    r_timing2 = R"""
    setkey(DT, id1, id2)
    system.time( DT[, sum(v1), keyby="id1,id2"])
    """
    # user  system elapsed
    # 1.24    0.07    1.33

    using Plots
    bar(
        ["FastGroupBy.jl", "DataFrames.jl","IndexedTables.jl (indexed)","R data.table", "R data.table (indexed)"],
        [0.401, 20.9, 2.555, r_timing[1], r_timing2[3]],
        title = "Group-by 2 factors with uniques: 100; len: 100m",
        label = "seconds",
        ylabel = "mean run time"
    )
    savefig("Group_by2_cate_perf.png")

end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder

name = "AzStorage"
version = v"0.2.1"

# Collection of sources required to build AzStorage
sources = [
    GitSource(
        "https://github.com/ChevronETC/AzStorage.jl.git",
        "db5e73ef1f6a2b5a58d258bab76202df3e29ff56"
    )
]

# Bash recipe for building across all platforms
script = raw"""
cd ${WORKSPACE}/srcdir/AzStorage.jl/src

# We need to tell the makefile where to find libssh2 on windows
if [[ ${target} == *mingw* ]]; then
    export LDFLAGS="-L${libdir}"
fi

make

cp libAzStorage.so ${libdir}/libAzStorage.${dlext}
"""

# 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
# TODO - add libgomp dependency
products = [
    LibraryProduct("libAzStorage", :libAzStorage)
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency("CompilerSupportLibraries_jll"),
    # libcurl changed compatibility version for macOS from v7.71 to v7.73 (v11
    # to v12)
    Dependency("LibCURL_jll", v"7.71.1"),
    # The following libraries are dependencies of LibCURL_jll which is now a
    # stdlib, but the stdlib doesn't explicitly list its dependencies
    Dependency("LibSSH2_jll"),
    Dependency("MbedTLS_jll", v"2.16.0"),
    Dependency("nghttp2_jll"),
    Dependency("Zlib_jll"),
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: blenessy/Yggdrasil

max_issues_repo_name: blenessy/Yggdrasil

max_forks_repo_name: blenessy/Yggdrasil


text:
using BioSequences, Distributed, DataFrames, nucpos, BioBackgroundModels, CSV, Serialization
import StatsBase:sample

#JOB PATHS AND CONSTANTS

sample_rows = 7092 #2 million bases

sib_pos = "/bench/PhD/danpos_results/pooled/sib.Fnor.smooth.positions.xls"
rys_pos = "/bench/PhD/danpos_results/pooled/rys.Fnor.smooth.positions.xls"

danio_genome_path = "/bench/PhD/seq/GRCz11/GCA_000002035.4_GRCz11_genomic.fna"
danio_gen_index_path = "/bench/PhD/seq/GRCz11/GCA_000002035.4_GRCz11_genomic.fna.fai"
danio_gff_path = "/bench/PhD/seq/GRCz11/Danio_rerio.GRCz11.94.gff3"

refined_folders_path = "/bench/PhD/NGS_binaries/BBM/refined_folders"
selected_hmms = "/bench/PhD/NGS_binaries/BBM/selected_hmms"

sib_df_binary = "/bench/PhD/NGS_binaries/BMI/sib_diff_positions"
rys_df_binary = "/bench/PhD/NGS_binaries/BMI/rys_diff_positions"
combined_df_binary = "/bench/PhD/NGS_binaries/BMI/combined_diff_positions"

sib_code_binary = "/bench/PhD/NGS_binaries/BMI/sib_diff_codes"
rys_code_binary = "/bench/PhD/NGS_binaries/BMI/rys_diff_codes"
combined_code_binary = "/bench/PhD/NGS_binaries/BMI/combined_diff_codes"

sib_diff_bg = "/bench/PhD/NGS_binaries/BMI/sib_diff_bg"
rys_diff_bg = "/bench/PhD/NGS_binaries/BMI/rys_diff_bg"
combined_diff_bg = "/bench/PhD/NGS_binaries/BMI/combined_diff_bg"

sib_fa = "/bench/PhD/thicweed_results/sib_nuc_position_sequences.fa"
sib_arch = "/bench/PhD/thicweed_results/sib_nuc_position_sequences.fa_archs.txt"

rys_fa = "/bench/PhD/thicweed_results/rys_nuc_position_sequences.fa"
rys_arch = "/bench/PhD/thicweed_results/rys_nuc_position_sequences.fa_archs.txt"

@info "Reading from position.xls..."
sib_df=CSV.read(sib_pos)
rys_df=CSV.read(rys_pos)

@info "Mapping..."
nucpos.map_positions!(sib_df, rys_df)
nucpos.map_positions!(rys_df, sib_df)

@info "Making differential position dataframes..."
sib_diff_df = deepcopy(sib_df[findall(iszero,sib_df.mapped_pos),:])
sib_diff_df = sib_diff_df[findall(!isequal("MT"), sib_diff_df.chr),:]
rys_diff_df = deepcopy(rys_df[findall(iszero,rys_df.mapped_pos),:])
rys_diff_df = rys_diff_df[findall(!isequal("MT"), rys_diff_df.chr),:]

nucpos.add_position_sequences!(sib_diff_df, danio_genome_path, danio_gen_index_path)
nucpos.add_position_sequences!(rys_diff_df, danio_genome_path, danio_gen_index_path)

@info "Filtering ambiguous sequences..."
deleterows!(sib_diff_df, [hasambiguity(seq) for seq in sib_diff_df.seq])
deleterows!(rys_diff_df, [hasambiguity(seq) for seq in rys_diff_df.seq])

@info "Sampling subsets of sequences..."

sib_diff_df=sib_diff_df[sample(axes(sib_diff_df, 1), sample_rows; replace=false, ordered=true), :]
rys_diff_df=rys_diff_df[sample(axes(rys_diff_df, 1), sample_rows; replace=false, ordered=true), :]

@info "Masking positions by genome partition and strand..."
BioBackgroundModels.add_partition_masks!(sib_diff_df, danio_gff_path, 500, (:chr,:seq,:start))
BioBackgroundModels.add_partition_masks!(rys_diff_df, danio_gff_path, 500, (:chr,:seq,:start))

@info "Obtaining cluster data..."
nucpos.get_cluster!(sib_diff_df, sib_fa, sib_arch)
nucpos.get_cluster!(rys_diff_df, rys_fa, rys_arch)

@info "Serializing dataframes..."
combined_diff_df = vcat(sib_diff_df, rys_diff_df)

serialize(sib_df_binary, sib_diff_df)
serialize(rys_df_binary, rys_diff_df)
serialize(combined_df_binary, combined_diff_df)

@info "Creating coded observation sets..."
sib_codes = nucpos.observation_setup(sib_diff_df, order=0, symbol=:seq)
rys_codes = nucpos.observation_setup(rys_diff_df, order=0, symbol=:seq)
combined_codes = nucpos.observation_setup(combined_diff_df, order=0, symbol=:seq)

@info "Serializing coded observation sets..."
serialize(sib_code_binary, Matrix(transpose(sib_codes)))
serialize(rys_code_binary, Matrix(transpose(rys_codes)))
serialize(combined_code_binary, Matrix(transpose(combined_codes)))

@info "Setting up for BBM likelihood calculations..."
BHMM_dict = Dict{String,BHMM}()
refined_folders=deserialize(refined_folders_path)
for (part, folder) in refined_folders
    BHMM_dict[part]=folder.partition_report.best_model[2]
end

serialize(selected_hmms,BHMM_dict)

@info "Performing calculations..."
sib_lh_matrix = BGHMM_likelihood_calc(sib_diff_df, BHMM_dict, symbol=:seq)
rys_lh_matrix = BGHMM_likelihood_calc(rys_diff_df, BHMM_dict, symbol=:seq)
combined_lh_matrix = hcat(sib_lh_matrix,rys_lh_matrix)

@info "Serializing background matrices..."
serialize(sib_diff_bg, sib_lh_matrix)
serialize(rys_diff_bg, rys_lh_matrix)
serialize(combined_diff_bg, combined_lh_matrix)

@info "Job done!"


text:
# ---
# title: Geometries
# cover: assets/geometries.png
# description: Visualizing geometries
# author: "[Pietro Vertechi](https://github.com/piever)"
# ---

using AlgebraOfGraphics, CairoMakie
using GeometryBasics
set_aog_theme!() #src

geometry = [Rect(Vec(i, j), Vec(1, 1)) for i in 0:7 for j in 0:7]
group = [isodd(i + j) ? "light square" : "dark square" for i in 0:7 for j in 0:7]
df = (; geometry, group)

plt = data(df) * visual(Poly) * mapping(:geometry, color = :group)
fg = draw(plt; axis=(aspect=1,))

# save cover image #src
mkpath("assets") #src
save("assets/geometries.png", fg) #src


text:
using ClearMethods
using Test

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


text:
# This file is a part of SimilaritySearch.jl
# License is Apache 2.0: https://www.apache.org/licenses/LICENSE-2.0.txt

export IHCSearch

# Iterated Hill Climbing Search
struct IHCSearch <: LocalSearchAlgorithm
    restarts::Int
    use_local_improvement::Bool
    IHCSearch() = new(1, false)
    IHCSearch(r; use_local_improvement=false) = new(r, use_local_improvement)
end

"""
    hill_climbing(index::SearchGraph, dist::Function, q, res::KnnResult, navigation_state, nodeID::Int64, use_local_improvement::Bool)

Runs a single hill climbing search process starting in vertex `nodeID`
"""
function hill_climbing(index::SearchGraph, dist::Function, q, res::KnnResult, navigation_state, nodeID::Int64, use_local_improvement::Bool)
    omin::Int = -1
    dmin::Float32 = typemax(Float32)

    while true
        dmin = typemax(Float32)
        omin = -1
        navigation_state[nodeID] = EXPLORED
        @inbounds for childID in index.links[nodeID]
            S = get(navigation_state, childID, UNKNOWN)
            S != UNKNOWN && continue
            navigation_state[childID] = VISITED
            d = convert(Float32, dist(index.db[childID], q))
            if use_local_improvement  ## this yields to better quality but can't be tuned for early stopping
                push!(res, childID, d)
                if d < dmin
                    dmin = d
                    omin = childID
                end
            else
                if push!(res, childID, d) && d < dmin
                    dmin = d
                    omin = childID
                end
            end
        end

        if omin < 0
            break
        else
            nodeID = omin
        end
    end
end

"""
    search(isearch::IHCSearch, index::SearchGraph, dist::Function, q, res::KnnResult, navigation_state, hints=EMPTY_INT_VECTOR)

Performs an iterated hill climbing search for `q`. The given `hints` are used as starting points of the search; a random
selection is performed otherwise.
"""
function search(isearch::IHCSearch, index::SearchGraph, dist::Function, q, res::KnnResult, navigation_state, hints=EMPTY_INT_VECTOR)
    n = length(index.db)
    restarts = min(isearch.restarts, n)
    if length(hints) == 0
        hints = rand(1:n, isearch.restarts)
    end

    @inbounds for start_point in hints
        # start_point = rand(range)
        S = get(navigation_state, start_point, UNKNOWN)
        if S == UNKNOWN
            navigation_state[start_point] = VISITED
            d = convert(Float32, dist(q, index.db[start_point]))
            push!(res, start_point, d)
            hill_climbing(index, dist, q, res, navigation_state, start_point, isearch.use_local_improvement)
        end
    end

    res
end

"""
    opt_expand_neighborhood(fun, algo::IHCSearch, n::Integer, iter::Integer, probes::Integer)

Generates configurations of the IHCSearch that feed the `optimize!` function (internal function)
"""
function opt_expand_neighborhood(fun, algo::IHCSearch, n::Integer, iter::Integer, probes::Integer)
    logn = ceil(Int, log(2, n+1))
    probes = probes == 0 ? logn : probes
    f(x) = max(1, x + rand(-logn:logn))

    for i in 1:probes
        IHCSearch(f(algo.restarts)) |> fun
    end
end


text:
@testset "Internal Sorting" begin
    data = [21, 12, 6, 26, 4, 15, 2, 8, 24, 22, 13, 30, 7, 9, 27, 1, 23, 19, 11, 14, 18, 28, 20, 25, 16, 17, 5, 3, 10, 29]
    @test insertionsort!(copy(data)) == 1:length(data)
    @test selectionsort!(copy(data)) == 1:length(data)
    @test shellsort!(copy(data)) == 1:length(data)

    @test mergesort!(copy(data)) == 1:length(data)
    @test quicksort!(copy(data)) == 1:length(data)

    @test heapsort!(copy(data)) == 1:length(data)
end


text:
module Dispersions

export AbstractDispersion,
TwoLevelSystem

abstract type AbstractDispersion end

"""
    tls = TwoLevelSystem(tls; :key1 => value1, :key2 => value2, ...)

new tls object from old, with modified fields
"""
struct TwoLevelSystem <: AbstractDispersion
    D₀::Float64
    k₀::Float64
    γp::Float64
    ω₀::Float64

    TwoLevelSystem(D₀=0, k₀=Inf, γp=1e8, ω₀=k₀) = new(D₀, k₀, γp, ω₀)
    TwoLevelSystem(tls::TwoLevelSystem; D₀=tls.D₀, k₀=tls.k₀, γp=tls.γp, ω₀=tls.ω₀) = TwoLevelSystem(tls.D₀, tls.k₀, tls.γp, tls.ω₀)

    function Base.show(io::IO, tls::TwoLevelSystem)
        if !get(io, :sub, false)
            print(io, "Two Level System:\n")
        end
        print(io, "\tD₀: ", tls.D₀, "\n",
        "\tω₀: ", tls.ω₀, "\n",
        "\tγ⟂: ", tls.γp)
    end
end

(tls::TwoLevelSystem)(ω::Number) = tls.γp/(ω-tls.ω₀+1im*tls.γp)
function (tls::TwoLevelSystem)(ω::Number,ωs::Array,ψ::Array)
    γω = tls(ω)
    h = zeros(Float64,size(ψ,1))
    for i ∈ eachindex(ωs)
        h += abs2(ωs[i])*abs2.(ψ[:,i])
    end
    return γω./(1 .+ h)
end

end # module


text:
using NeuralPDE, Flux, ModelingToolkit, GalacticOptim, Optim, DiffEqFlux
using Quadrature, Cuba, CUDA, QuasiMonteCarlo
import ModelingToolkit: Interval, infimum, supremum
using DaggerFlux

@parameters t x y
@variables u(..)
Dxx = Differential(x)^2
Dyy = Differential(y)^2
Dt = Differential(t)
t_min = 0.0
t_max = 2.0
x_min = 0.0
x_max = 2.0
y_min = 0.0
y_max = 2.0

# 2D PDE
eq  = Dt(u(t,x,y)) ~ Dxx(u(t,x,y)) + Dyy(u(t,x,y))

analytic_sol_func(t,x,y) = exp(x+y)*cos(x+y+4t)
# Initial and boundary conditions
bcs = [u(t_min,x,y) ~ analytic_sol_func(t_min,x,y),
       u(t,x_min,y) ~ analytic_sol_func(t,x_min,y),
       u(t,x_max,y) ~ analytic_sol_func(t,x_max,y),
       u(t,x,y_min) ~ analytic_sol_func(t,x,y_min),
       u(t,x,y_max) ~ analytic_sol_func(t,x,y_max)]

# Space and time domains
domains = [t ∈ Interval(t_min,t_max),
           x ∈ Interval(x_min,x_max),
           y ∈ Interval(y_min,y_max)]

# Neural network

# Automatically calls `collect` on a `DaggerChain`
struct AutoChain
    inner
end
(ac::AutoChain)(x...) = collect(ac.inner(x...))

inner = 25
chain = Chain(Dense(3,inner,Flux.σ),
              Dense(inner,inner,Flux.σ),
              Dense(inner,inner,Flux.σ),
              Dense(inner,inner,Flux.σ),
              Dense(inner,1))
chain = AutoChain(DaggerChain(chain))
#= TODO: Use FastChain
chain = FastChain(FastDense(3,inner,Flux.σ),
                  FastDense(inner,inner,Flux.σ),
                  FastDense(inner,inner,Flux.σ),
                  FastDense(inner,inner,Flux.σ),
                  FastDense(inner,1))
chain = AutoChain(DaggerChain(Chain(chain)))

#initθ = CuArray(Float64.(DiffEqFlux.initial_params(chain)))
initθ = Float64.(DiffEqFlux.initial_params(chain))
=#

strategy = GridTraining(0.05)
discretization = PhysicsInformedNN(chain,
                                   strategy)#;
                                   #init_params = initθ)

@named pde_system = PDESystem(eq,bcs,domains,[t,x,y],[u(t, x, y)])
prob = discretize(pde_system,discretization)
symprob = symbolic_discretize(pde_system,discretization)

cb = function (p,l)
    println("Current loss is: $l")
    return false
end

res = GalacticOptim.solve(prob, ADAM(0.01); cb=cb, maxiters=2500)

max_stars_repo_name: jpsamaroo/DaggerFluxPINNs

max_issues_repo_name: jpsamaroo/DaggerFluxPINNs

max_forks_repo_name: jpsamaroo/DaggerFluxPINNs


text:
#  ------------------------------------------------------------------------------------------
#  Copyright (c) Microsoft Corporation. All rights reserved.
#  Licensed under the MIT License (MIT). See LICENSE in the repo root for license information.
#  ------------------------------------------------------------------------------------------

# Exports
export create_batch_client, create_blob_client, create_queue_client, create_clients
export create_blob_containers, create_pool_and_resource_file, resize_pool, create_pool
export create_batch_resource_from_file, create_batch_resource_from_bytes, create_batch_resource_from_blob
export create_batch_job, create_batch_task, submit_batch_job, create_batch_env
export wait_for_tasks_to_complete, wait_for_task_to_complete#, wait_for_one_task_to_complete
export create_batch_output_file, create_task_constraint, enable_auto_scale, create_batch_envs
export wait_for_one_task_from_multi_jobs, wait_for_one_task_from_multi_pool
export upload_bytes_to_container, create_blob_url, create_batch_resource_from_blob_url

###################################################################################################
# Wrappers around Python credential functions

function create_batch_client(credentials::Dict{String, Any})
    return azureclusterlesshpc.create_batch_client(credentials)
end

function create_blob_client(credentials::Dict{String, Any})
    return azureclusterlesshpc.create_blob_client(credentials)
end

function create_queue_client(credentials::Dict{String, Any})
    return azureclusterlesshpc.queue_blob_client(credentials)
end

function create_clients(credentials::Dict{String, Any}; batch::Bool=false, blob::Bool=false, queue::Bool=false)

    batch ? (batch_client = create_batch_client(credentials)) : batch_client = nothing
    blob ? (blob_client = create_blob_client(credentials)) : blob_client = nothing
    queue ? (queue_client = create_queue_client(credentials)) : queue_client = nothing

    return [Dict(
        "batch_client" => batch_client,
        "blob_client" => blob_client,
        "queue_client" => queue_client
        )]
end

function create_clients(credential_list::Array; batch::Bool=false, blob::Bool=false, queue::Bool=false)

    # List of all clients
    clients = Array{Any}(undef, 0)

    for credential in credential_list

        batch ? (batch_client = create_batch_client(credential)) : batch_client = nothing
        blob ? (blob_client = create_blob_client(credential)) : blob_client = nothing
        queue ? (queue_client = create_queue_client(credential)) : queue_client = nothing

        cred = Dict(
            "batch_client" => batch_client,
            "blob_client" => blob_client,
            "queue_client" => queue_client
            )

        push!(clients, cred)
    end
    return clients
end

###################################################################################################
# Blob stuff

create_blob_url(blob_client::PyObject, container_name, blob_list) =
    azureclusterlesshpc.create_blob_url(blob_client, container_name, blob_list)

upload_bytes_to_container(blob_client::PyObject, container_name, blob_name, blob; verbose=true) =
    azureclusterlesshpc.upload_bytes_to_container(blob_client, container_name, blob_name, blob; verbose=verbose)

# Create containers given a list of container names
create_blob_containers(blob_client::PyObject, container_name_list::Array{String, 1}) =
    azureclusterlesshpc.create_blob_containers(blob_client, container_name_list)

create_batch_output_file(blob_client, storage_account_name, container_name, filename) =
    azureclusterlesshpc.create_batch_output_file(blob_client, storage_account_name, container_name, filename)

###################################################################################################
# Batch stuff

# Create pool
function create_pool(batch_service_client, pool_id, pool_vm_size, pool_node_count, node_os_publisher,
    node_os_offer, node_os_sku; image_resource_id=nothing, enable_inter_node=false, resource_files=nothing,
    enable_auto_scale=false, auto_scale_formula=nothing, auto_scale_evaluation_interval_minutes=nothing,
    container=nothing, container_registry=nothing)

    azureclusterlesshpc.create_pool(batch_service_client, pool_id, pool_vm_size, pool_node_count, node_os_publisher,
        node_os_offer, node_os_sku, image_resource_id=image_resource_id, enable_inter_node=enable_inter_node,
        resource_files=resource_files, enable_auto_scale=enable_auto_scale, auto_scale_formula=auto_scale_formula,
        auto_scale_evaluation_interval_minutes=auto_scale_evaluation_interval_minutes, container=container,
        container_registry=container_registry)
end

# Create pool and resource file
function create_pool_and_resource_file(clients, pool_id, pool_vm_size, pool_node_count, node_os_publisher,
    node_os_offer, node_os_sku, file_name, container_name; image_resource_id=nothing, enable_inter_node=false,
    enable_auto_scale=false, auto_scale_formula=nothing, auto_scale_evaluation_interval_minutes=15, container=false,
    container_registry=nothing)

    azureclusterlesshpc.create_pool_and_resource_file(clients, pool_id, pool_vm_size, pool_node_count,
        node_os_publisher, node_os_offer, node_os_sku, file_name, container_name, image_resource_id=image_resource_id,
        enable_inter_node=enable_inter_node, enable_auto_scale=enable_auto_scale, auto_scale_formula=auto_scale_formula,
        auto_scale_evaluation_interval_minutes=auto_scale_evaluation_interval_minutes, container=container,
        container_registry=container_registry)
end

# Enable auto scaling for batch pool
enable_auto_scale(batch_client, pool_id, auto_scale_formula; auto_scale_evaluation_interval_minutes=5) =
    azureclusterlesshpc.enable_auto_scale(batch_client, pool_id, auto_scale_formula,
        auto_scale_evaluation_interval_minutes=auto_scale_evaluation_interval_minutes)

# Resize pool
resize_pool(batch_client, pool_id, target_dedicated_nodes, target_low_priority_nodes;
    resize_timeout_minutes=nothing, node_deallocation_option=nothing, pool_resize_options=nothing) =
    azureclusterlesshpc.resize_pool(batch_client, pool_id, target_dedicated_nodes, target_low_priority_nodes,
        resize_timeout_minutes=resize_timeout_minutes, node_deallocation_option=node_deallocation_option,
        pool_resize_options=pool_resize_options)

# Create resource from url
create_batch_resource_from_blob_url(shared_url, shared_blob) =
    azureclusterlesshpc.create_batch_resource_from_blob_url(shared_url, shared_blob)

# Create resource file from file
create_batch_resource_from_file(blob_client, container, file; verbose=true) =
    azureclusterlesshpc.create_batch_resource_from_file(blob_client, container, file, verbose=verbose)

# Create resource file from existing blob
create_batch_resource_from_blob(blob_client, container, blob) =
    azureclusterlesshpc.create_batch_resource_from_blob(blob_client, container, blob)

# Create resource file from bytes
create_batch_resource_from_bytes(blob_client, container, blob_name, blob; verbose=true) =
    azureclusterlesshpc.create_batch_resource_from_bytes(blob_client, container, blob_name, blob, verbose=verbose)

# Create batch job
create_batch_job(batch_client, job_id, pool_id; uses_task_dependencies=false, priority=0, verbose=true) =
    azureclusterlesshpc.create_batch_job(batch_client, job_id, pool_id, uses_task_dependencies=uses_task_dependencies,
        priority=priority, verbose=verbose)

# Create tasks for batch job
create_batch_task(; resource_files=nothing, environment_variables=nothing, application_cmd=nothing,
    output_files=nothing, taskname="task", task_constraints=nothing, num_nodes_per_task=1, docker_container=nothing) =
    azureclusterlesshpc.create_batch_task(resource_files=resource_files, environment_variables=environment_variables,
        application_cmd=application_cmd, output_files=output_files, taskname=taskname, task_constraints=task_constraints,
        num_nodes_per_task=num_nodes_per_task, docker_container=docker_container)

# Create task constraints
create_task_constraint(; max_wall_clock_time=nothing, retention_time=nothing, max_task_retry_count=0) =
    azureclusterlesshpc.create_task_constraint(max_wall_clock_time=max_wall_clock_time, retention_time=retention_time,
        max_task_retry_count=max_task_retry_count)

# Wait for all tasks to complete
wait_for_tasks_to_complete(batch_service_client, job_id;  task_timeout=60, fetch_timeout=60, verbose=true, num_restart=0) =
    azureclusterlesshpc.wait_for_tasks_to_complete(batch_service_client, job_id, task_timeout=task_timeout,
        fetch_timeout=fetch_timeout, verbose=verbose, num_restart=num_restart)

# Wait for specified task to complete
wait_for_task_to_complete(batch_service_client, job_id, task_id, timeout; verbose=true, num_restart=0) =
    azureclusterlesshpc.wait_for_task_to_complete(batch_service_client, job_id, task_id, timeout, verbose=verbose,
    num_restart=num_restart)

# Wait for one task from a list of tasks to complete
wait_for_one_task_from_multi_pool(batch_service_client, job_id, task_id_list;
    task_timeout=60, fetch_timeout=60, verbose=true, num_restart=0) =
    azureclusterlesshpc.wait_for_one_task_from_multi_pool(batch_service_client, job_id, task_id_list,
    task_timeout=task_timeout, fetch_timeout=fetch_timeout, verbose=verbose, num_restart=num_restart)

wait_for_one_task_from_multi_jobs(batch_service_client, job_id_list, task_id_list;
    task_timeout=60, fetch_timeout=60, verbose=true, num_restart=0) =
    azureclusterlesshpc.wait_for_one_task_from_multi_jobs(batch_service_client, job_id_list, task_id_list,
    task_timeout=task_timeout, fetch_timeout=fetch_timeout, verbose=verbose, num_restart=num_restart)

# Create batch environment variable
create_batch_env(name, value) = azureclusterlesshpc.create_batch_env(name, value)
create_batch_envs(names, values) = azureclusterlesshpc.create_batch_envs(names, values)


text:
export submesh

function submesh(mesh::VolumeMesh{Dim,T,U}, name::String) where {Dim,T,U}
    # Submesh elements prior to point id remap
    element_ids = mesh.groups[name]
    nelements = length(element_ids)
    types = Vector{U}(undef, nelements)
    offsets = Vector{U}(undef, nelements)
    materials = Vector{UInt8}(undef, nelements)
    connectivity_len = 1
    for (i, id) in enumerate(element_ids)
        types[i] = mesh.types[id]
        offsets[i] = connectivity_len
        materials[i] = mesh.materials[id]
        connectivity_len += points_in_vtk_type(types[i])
    end
    connectivity = Vector{U}(undef, connectivity_len - 1)
    for (i, id) in enumerate(element_ids)
        npts = points_in_vtk_type(types[i])
        sm_offset = offsets[i]
        m_offset = mesh.offsets[id]
        connectivity[sm_offset:sm_offset+npts-1] = mesh.connectivity[m_offset:m_offset+npts-1]
    end

    # Submesh points
    point_ids = BitSet()
    for vid in connectivity
        push!(point_ids, vid)
    end
    npoints = length(point_ids)
    points = Vector{Point{Dim,T}}(undef, npoints)
    for (i, vid) in enumerate(point_ids)
        points[i] = mesh.points[vid]
    end

    # Remap point ids in connectivity
    point_ids_vec = collect(point_ids)
    for (i, vid) in enumerate(connectivity)
        connectivity[i] = searchsortedfirst(point_ids_vec, vid)
    end

    # Submesh groups prior to element id remap
    groups = Dict{String,BitSet}()
    for group_name in keys(mesh.groups)
        set_intersection = mesh.groups[group_name] ∩ element_ids
        if length(set_intersection) !== 0 && name !== group_name
            groups[group_name] = set_intersection
        end
    end

    # Remap element ids in groups
    element_ids_vec = collect(element_ids)
    for group_name in keys(groups)
        groups[group_name] = BitSet(
                                searchsortedfirst.(Ref(element_ids_vec),
                                                   groups[group_name]
                                )
                             )
    end
    return VolumeMesh{Dim,T,U}(points, offsets, connectivity, types,
                               materials, mesh.material_names, name, groups)
end

function submesh(mesh::PolytopeVertexMesh, name::String)
    # Submesh polytopes prior to vertex id remap
    polytope_ids = mesh.groups[name]
    polytopes = [mesh.polytopes[pid] for pid in polytope_ids]

    # Submesh points
    vertex_ids = BitSet()
    for polytope in polytopes
        for vid in vertices(polytope)
            push!(vertex_ids, vid)
        end
    end
    npoints = length(vertex_ids)
    point_type = typeof(vertices(mesh)[1])
    points = Vector{point_type}(undef, npoints)
    for (i, vid) in enumerate(vertex_ids)
        points[i] = mesh.vertices[vid]
    end

    # Remap vertex ids in polytopes
    vertex_ids_vec = collect(vertex_ids)
    for i in eachindex(polytopes)
        polytopes[i] = searchsortedfirst.(Ref(vertex_ids_vec), vertices(polytopes[i]))
    end

    # Submesh groups prior to polytope id remap
    groups = Dict{String,BitSet}()
    for group_name in keys(mesh.groups)
        set_intersection = mesh.groups[group_name] ∩ polytope_ids
        if length(set_intersection) !== 0 && name !== group_name
            groups[group_name] = set_intersection
        end
    end

    # Remap polytope ids in groups
    polytope_ids_vec = collect(polytope_ids)
    for group_name in keys(groups)
        groups[group_name] = BitSet(
                                searchsortedfirst.(Ref(polytope_ids_vec),
                                                   groups[group_name]
                                )
                             )
    end

    materials = Vector{UInt8}(undef, length(polytope_ids_vec))
    for (i, id) in enumerate(polytope_ids_vec)
        materials[i] = mesh.materials[id]
    end
    return typeof(mesh)(points, polytopes, materials,
                        mesh.material_names, mesh.name, groups)
end

max_stars_repo_name: khurrumsaleem/MOCNeutronTransport

max_issues_repo_name: khurrumsaleem/MOCNeutronTransport

max_forks_repo_name: khurrumsaleem/MOCNeutronTransport


text:
include("Utilities.jl")
include("VEM.jl")
include("SuSiEGLM.jl")

function init(yt::Vector{Float64},Xt₀::Matrix{Float64},S::Vector{Float64},β::Vector{Float64},ξ::Vector{Float64},τ²::Float64;tol=1e-4)

     res= emGLMM(yt,Xt₀,S,τ²,β,ξ;tol=tol)

    return res

end

"""

    initialization(y::Vector{Float64},X₀::Union{Matrix{Float64},Vector{Float64}},T::Matrix{Float64},
        S::Vector{Float64};tol=1e-4)

Returns initial values for parameters `τ2, β, ξ` to run fine-mapping for SuSiEGLMM.

# Arguments

- `y` : a n x 1 vector of binary trait
- `X₀` : a n x c matrix of covariates.  The intercept is default if no covariates is added.
- `T` : a matrix of eigen-vectors by eigen decomposition to K (kinship)
- `S` : a vecor of eigen-values by eigen decomposition to K (kinship)

## Keyword Arguments

- `tol` : tolerance. Default value is `1e-4`.

# Output

- `Xt₀`: a n x c transformed covariate matrix
- `yt`: a transformed trait vector
-  Type of struct, `Null_est` :  initial values of parameters obtained by EM algorithm, or a null model of the SuSiE-GLMM excluding SuSiE implementation.  It includes following estimates:

    - `τ2` : a constant variance in the variance component of genetic relatedness, `τ²K`, where K is a kinship matrix
    -  `β` : a c x 1 vector of fixed effects for covariates
    -  `ξ` : a n x 1 vector of variational parameters to fit a mixed logitstic function

"""
function initialization(y::Vector{Float64},X::Matrix{Float64},X₀::Union{Matrix{Float64},Vector{Float64}},T::Matrix{Float64},
        S::Vector{Float64};tol=1e-4)

        n=length(y)
    # check if covariates are added as input and include the intercept.
    if(X₀!= ones(n,1)) #&&(size(X₀,2)>1)
        X₀ = hcat(ones(n),X₀)
    end

#                    Xt₀ = rotateX(X₀,T)
#                    yt = rotateY(y,T)
     Xt, Xt₀, yt = rotate(y,X,X₀,T)
    #  y0= getXy('N',T,y) # rotate w/o centering for β0
    #initialization
     τ0 = 0.0001 #rand(1)[1]; #arbitray
    # τ0=1.2
    β0 = glm(X₀,y,Binomial()) |> coef
    ξ0 =sqrt.(getXy('N',Xt₀,β0).^2+ τ0*S)

    init_est= init(yt,Xt₀,S,β0,ξ0,τ0;tol=tol)

    return Xt, Xt₀, yt, init_est
end

# version 1
function susieGLMM(L::Int64,Π::Vector{Float64},yt::Vector{Float64},Xt::Matrix{Float64},Xt₀::Matrix{Float64},
        S::Vector{Float64},est0::Null_est;tol=1e-4)

    # n, p = size(Xt)
    #initialization :
     σ0 = 0.1*ones(L);

     result = emGLMM(L,yt,Xt,Xt₀,S,est0.τ2,est0.β,est0.ξ,σ0,Π;tol=tol)

    return result

end

# version 2
function susieGLMM(L::Int64,Π::Vector{Float64},y::Vector{Float64},X::Matrix{Float64},X₀::Union{Matrix{Float64},Vector{Float64}},T::Matrix{Float64},
    S::Vector{Float64};tol=1e-4)

# check if covariates are added as input and include the intercept.
    #  n=length(y)
    #     if(X₀!= ones(n,1)) #&&(size(X₀,2)>1)
    #         X₀ = hcat(ones(n),X₀)
    #     end

    Xt, Xt₀, yt = rotate(y,X,X₀,T)
    # #initialization :
    σ0 = 0.1*ones(L);
     τ0 = 0.001  #rand(1)[1]; #arbitray
    β0 = glm(X₀,y,Binomial()) |> coef
    ν0 =sum(repeat(Π,outer=(1,L)).*σ0',dims=2)[:,1] ; #ν²0
    ξ0 =sqrt.(getXy('N',Xt.^2.0,ν0)+ getXy('N',Xt₀,β0).^2+ τ0*S )

    # res0 = susieGLM(L,Π,y,X,X₀;tol=1e-2)
    # ξ0 = sqrt.(res0.ξ.^2+τ0*S)

    result = emGLMM(L,yt,Xt,Xt₀,S,τ0,β0,ξ0,σ0,Π;tol=tol)
    # result = emGLMM(L,yt,Xt,Xt₀,S,τ0,res0.β,ξ0,res0.σ0,Π;tol=tol)

return result

end

# compute score test statistic : need to check again
function computeT(init0::Null_est,yt::Vector{Float64},Xt₀::Matrix{Float64},Xt::Matrix{Float64})

        m=axes(Xt,2)
        r₀ =  2*yt.*(getXy('N',Xt₀,init0.β)+init0.μ)
        p̂ = logistic.(r₀)
        Γ  = p̂.*(1.0.-p̂)
        # XX=Xt₀'Diagonal(Γ)
        proj= I - Xt₀*(symXX('T',sqrt.(Γ).*Xt₀)\(Xt₀'Diagonal(Γ)))
        # proj = I - Xt₀*(getXX('N',XX,'N',Xt₀))\XX
        G̃ = getXX('N',proj,'N',Xt)

        Tstat= zeros(m)

         ĝ = getXy('T',G̃,yt-p̂).^2

           for j = m

              Tstat[j] = ĝ[j]/(G̃[:,j]'*(Γ.*G̃[:,j]))
           end

    return Tstat
end

"""

    scoreTest(K::Union{Array{Float64,3},Matrix{Float64}},G::GenoInfo,y::Vector{Float64},X::Matrix{Float64},
    X₀::Union{Matrix{Float64},Vector{Float64}}=ones(length(y),1);LOCO::Bool=true,tol=1e-4)

Returns 1-df score test statistic for case-control association tests that follows Chi-square distribution, `T²∼Χ²₍₁₎` under `H₀: β=0`.

# Arguments

- `K` : a 3d-array of a symmetric positive definite kinship, `size(K) = (n,n, # of Chromosomes)` if `LOCO = true` (default).  Or a n x n matrix if false.
- `G` : a Type of struct, `GenoInfo`. See [`GenoInfo`](@ref).
- `y` : a n x 1 vector of  binary trait
- `X` : a n x p matrix of genetic markers (SNPs)
- `X₀`: a n x c matrix of covariates.  The intercept is default if no covariates is added.

## Keyword Arguments

- `LOCO`: Boolean. Default is `true` performs score test according to the Leave One Chromosome Out (LOCO) scheme.
- `tol` : tolerance. Default is `1e-4`.

# Output

- `T_stat` : a p x 1 vector of test statistics for case-control association tests.  The  test statistic is computed based on [reginie](https://www.nature.com/articles/s41588-021-00870-7).
- `p_value`: a p x 1 vector of the corresponding p-values (by Chi^2 with 1 df) for case-control association tests.

"""
function scoreTest(K::Union{Array{Float64,3},Matrix{Float64}},G::GenoInfo,y::Vector{Float64},X::Matrix{Float64},X₀::Union{Matrix{Float64},Vector{Float64}}=ones(length(y),1)
        ;LOCO::Bool=true,tol=1e-4)

        Chr=sort(unique(G.chr));
        T, S = svdK(K;LOCO=LOCO)
        # println("Eigen-decomposition is completed.")

    if (LOCO)

        # T_stat=zeros(size(X,2))

         T_stat = @distributed (vcat) for j= eachindex(Chr)
                       midx= findall(G.chr.== Chr[j])
                       Xt, Xt₀, yt, init0 = initialization(y,X[:,midx],X₀,T[:,:,j],S[:,j];tol=tol)
                       tstat = computeT(init0,yt,Xt₀,Xt)
                        tstat
                     end

    else #no loco

        Xt, Xt₀, yt, init0 = initialization(y,X,X₀,T,S;tol=tol)

        T_stat = @distributed (vcat) for j= eachindex(Chr)
            midx= findall(G.chr.== Chr[j])
            tstat = computeT(init0,yt,Xt₀,Xt[:,midx])
            tstat
         end

    end

    p_value = ccdf.(Chisq(1),T_stat) #check later for parallelization

    return T_stat, p_value
end

"""

     fineQTL_glmm(G::GenoInfo,y::Vector{Float64},X::Matrix{Float64},
        X₀::Union{Matrix{Float64},Vector{Float64}},
        T::Union{Array{Float64,3},Matrix{Float64}},S::Union{Matrix{Float64},Vector{Float64}};
        LOCO::Bool=true,L::Int64=10,Π::Vector{Float64}=[1/size(X,2)],tol=1e-4)

Performs SuSiE (Sum of Single Effects model) GLMM fine-mapping analysis for a binary trait (logistic mixed model).

# Arguments

- `G` : a Type of struct, `GenoInfo`. See [`GenoInfo`](@ref).
- `y` : a n x 1 vector of  binary trait
- `X` : a n x p matrix of genetic markers selected from QTL analysis (per Chromosome for LOCO)
- `X₀`: a n x 1 vector or n x c matrix of covariates.  The intercept is default if no covariates is added.
- `T` : a matrix of eigen-vectors by eigen decomposition to K (kinship)
- `S` : a vecor of eigen-values by eigen decomposition to K (kinship)

## Keyword Arguments

- `LOCO`: Boolean. Default is `true` indicates fine-mapping performs according to the Leave One Chromosome Out (LOCO) scheme.
- `L` : an integer. The number of single effects for SuSiE implementation. Default is `10`.
- `Π` : a p x 1 vector of prior inclusion probabilities for SuSiE.  Default is `1/p`, where `p = size(X,2)`. If different probabilities are added to SNPs, the length of Π should be `p`.
- `tol`: tolerance. Default is `1e-4`.

# Output

 Returns a Type of struct, `Result`, which includes

-  `ξ` : a n x 1 vector of variational parameters to fit a mixed logitstic function
-  `β` : a c x 1 vector of fixed effects for covariates
-  `σ0` : a L x 1 vector of hyper-parameters for prior variances for SuSiE
-  `τ2` : a constant variance in the variance component of genetic relatedness, `τ²K`, where K is a kinship matrix
-  `α` : p x L matrix of posterior inclusion probabilities for SuSiE
-  `ν` : p x L matrix of posterior mean of SuSiE
-  `Σ` : p x L matrix of posterior variances of SuSiE

"""
function fineQTL_glmm(G::GenoInfo,y::Vector{Float64},X::Matrix{Float64},
        X₀::Union{Matrix{Float64},Vector{Float64}},
        T::Union{Array{Float64,3},Matrix{Float64}},S::Union{Matrix{Float64},Vector{Float64}};
        LOCO::Bool=true,L::Int64=10,Π::Vector{Float64}=[1/size(X,2)],tol=1e-4)

        Chr=sort(unique(G.chr));

     if(LOCO)

   est= @distributed (vcat) for j= eachindex(Chr)
                midx= findall(G.chr.== Chr[j])
                Xt, Xt₀, yt, init0 = initialization(y,X[:,midx],X₀,T[:,:,j],S[:,j];tol=tol)
                           #check size of Π

                          if (Π==[1/size(X,2)]) #default value
                              m=length(midx)
                              Π1 =ones(m)/m #adjusting πⱼ
                             elseif (length(Π)!= size(X,2))
                                println("Error. The length of Π should match $(size(X,2)) SNPs!")
                             else
                              Π1 = Π[midx]
                           end

                est0= susieGLMM(L,Π1,yt,Xt,Xt₀,S[:,j],init0;tol=tol)
                        est0
                           end

    else #no loco for one genomic region

            est = @distributed (vcat) for j= eachindex(Chr)
               midx= findall(G.chr.== Chr[j])

               # check prior probabilities
                 if (Π==[1/size(X,2)])
                    m=length(midx)
                    Π1 =ones(m)/m
                  elseif (length(Π)!= size(X,2))
                    println("Error. The length of Π should match $(size(X,2)) SNPs!")
                   else
                    Π1 = Π[midx]
                 end

                 Xt, Xt₀, yt, init0 = initialization(y,X[:,midx],X₀,T,S;tol=tol)
                 est0 = susieGLMM(L,Π1,yt,Xt,Xt₀,S,init0;tol=tol)
                 est0
            end

    end # loco

    return est

end

#version 2
function fineQTL_glmm(K::Union{Array{Float64,3},Matrix{Float64}},G::GenoInfo,y::Vector{Float64},X::Matrix{Float64},
    X₀::Union{Matrix{Float64},Vector{Float64}}=ones(length(y),1);
    LOCO::Bool=true,L::Int64=10,Π::Vector{Float64}=[1/size(X,2)],tol=1e-4)

    Chr=sort(unique(G.chr));
    T, S = svdK(K;LOCO=LOCO)

 if(LOCO)

     est= @distributed (vcat) for j= eachindex(Chr)
            midx= findall(G.chr.== Chr[j])

                       #check size of Π

                      if (Π==[1/size(X,2)]) #default value
                          m=length(midx)
                          Π1 =ones(m)/m #adjusting πⱼ
                         elseif (length(Π)!= size(X,2))
                            println("Error. The length of Π should match $(size(X,2)) SNPs!")
                         else
                          Π1 = Π[midx]
                       end

            est0= susieGLMM(L,Π1,y,X[:,midx],X₀,T[:,:,j],S[:,j];tol=tol)
            est0
     end

else #no loco

        est = @distributed (vcat) for j= eachindex(Chr)
           midx= findall(G.chr.== Chr[j])

           # check prior probabilities
             if (Π==[1/size(X,2)])
                m=length(midx)
                Π1 =ones(m)/m
              elseif (length(Π)!= size(X,2))
                println("Error. The length of Π should match $(size(X,2)) SNPs!")
               else
                Π1 = Π[midx]
             end

             est0 = susieGLMM(L,Π1,y,X[:,midx],X₀,T,S;tol=tol)
             est0
        end

end # loco

return est

end

function randomizeData!(y::Vector{Float64},X::Matrix{Float64},X₀::Matrix{Float64},
    T::Array{Float64,3},S::Matrix{Float64})

    #   Random.seed!(seed)
    ivec = Vector(1:length(yt))
    #randomize minibatches
    shuffle!(MersenneTwister(123), ivec) # shuffle index
    y[:] = y[ivec]
    X[:,:] = X[ivec,:]
    X₀[:,:] = X₀[ivec,:]
    T[:,:,:] = T[ivec,:,:]
    S[:,:] = S[ivec,:]

end

function randomizeData!(y::Vector{Float64},X::Matrix{Float64},X₀::Matrix{Float64},
    T::Array{Float64,2},S::Vector{Float64})

    #   Random.seed!(seed)
    ivec = Vector(1:length(yt))
    #randomize minibatches
    shuffle!(MersenneTwister(123), ivec) # shuffle index
    y[:] = y[ivec]
    X[:,:] = X[ivec,:]
    X₀[:,:] = X₀[ivec,:]
    T[:,:] = T[ivec,:]
    S[:] = S[ivec]

end

# function fineMapping_miniBatch(y,X,X₀,K;bsize::Int64=256,LOCO::Bool=true,L::Int64=10,Π::Vector{Float64}=[1/size(X,2)],tol=1e-5)

#     M= length(y)
#     num_epoch = floor(M/bsize)

#     T,S = svdK(K;LOCO=LOCO)
#     println("Eigen-decomposition is completed.")
#     randomizeData!(y,X,X₀,T,S)

# end

function fineMapping(G::GenoInfo,y::Vector{Float64},X::Matrix{Float64},X₀::Union{Matrix{Float64},Vector{Float64}}=ones(length(y),1);
        K::Union{Array{Float64,3},Matrix{Float64}}=Matrix(1.0I,1,1),L::Int64=10,Π::Vector{Float64}=[1/size(X,2)],LOCO::Bool=true,
        model=["susieglmm","susieglm","susie","mvsusie"],tol=1e-4)

    #need to work more

    if(model=="susieglmm")

        T, S = svdK(K;LOCO=LOCO)

        println("Eigen-decomposition is completed.")

        est = fineQTL_glmm(G,y,X,X₀,T,S;L=L,Π=Π,LOCO=LOCO,tol=tol)
            println("SuSiEGLMM is completed.")

    elseif(model=="susieglm")
        est = fineQTL_glm(G,y,X,X₀,L=L,Π=Π,tol=tol)

    else #model=mvsusie

        println("it is not ready yet.")

    end #model

    return est
end

function fineMapping1(f::Function,args...;kwargs...)

  #need to work more : using splats to give more freedom gets slow in performance.
        res = f(args...;kwargs...)

    return res

end

export init, initialization, fineQTL_glmm, susieGLMM, computeT, scoreTest, GenoInfo


text:
# Give all the arguments the names they have in the signature
function insert_function_preamble(body, def, call_args_typetuple)
    body = Expr(:block, body)
    if !isempty(def[:args])  # give args their local names
        argnames = first.(MacroTools.splitarg.(def[:args]))
        pushfirst!(body.args,
            Expr(:(=),
                Expr(:tuple, argnames...),
                :call_args
            )
        )
    end

    if !isempty(def[:whereparams])  # give typeparams local names
        type_params = determine_type_params(meth, call_args_typetuple)
        pushfirst!(body.args,
            Expr(:(=),
                Expr(:tuple, wherevar.(def[:whereparams])...),
                :type_params
            )
        )
    end
    return body
end

function recurse!(expr, reflection)
    @assert isexpr(expr, :call)
    expr.args[1] = qualify_name(
        reflection.meth.module, expr.args[1]
    )
    return Expr(:block,
        Expr(:call,
            :graft, reflection.grafter,
            expr.args...
        )
    )
end

# Can't graft intrinsics
graft(grafter, f::Core.IntrinsicFunction, args...) = f(args...)
# for prototyping purpose ban grafting on constructors
graft(grafter, f::Type, args...) = f(args...)

struct Reflection
    grafter
    meth::Method
    ast
    def::Dict{Symbol}
end

"""
As generated functions can't call things defined after them
we call this function to redefine `graft`.
Do this after you have defined your grafters and everything
they call.
"""
function redeclare_graft()
    @eval @inline @generated function graft(grafter, f, call_args...)
#        Core.println(string("*** ", f))
        Core.println(call_args)
        call_args_typetuple = Tuple{call_args...}
        meth = get_method(f, call_args_typetuple)
        Core.println(string("*** ", meth))

        ast = get_ast(meth)
        if ast===nothing
            return quote
                # for dev purposes
                #@warn "No AST got. Directly invoking."
                f(call_args...)
            end
        end

        def = MacroTools.splitdef(ast)
        body = def[:body]

        #TODO: Replace the anon functions

        # Do the actual replacement (grafting)
        reflection = Reflection(grafter.instance, meth, body, def)
        body = grafter.instance(reflection)

        body = qualify_calls!(body, meth.module)
        body = insert_function_preamble(body, def, call_args_typetuple)
        body = MacroTools.flatten(body)

        #Core.println(sprint((io,x)->dump(io, x;maxdepth=100), body)); Core.println()
        Core.println(repr(body))
        Core.println()

        return body
    end
end

redeclare_graft() #We are actually defining it the first time here


text:
function goToDirect(r, side)
    while !isborder(r, side)
        move!(r, side)
        putmarker!(r)
    end
end

function goToStart(r)
    while !isborder(r, Sud)
        move!(r, Sud)
    end
    while !isborder(r, West)
        move!(r, West)
    end
end
function drawAll(r)
    top = 0
    right = 0

    while !isborder(r, Sud) || !isborder(r, West)
        if !isborder(r, Sud)
            top += 1
            move!(r, Sud)
        end
        if !isborder(r, West)
            right += 1
            move!(r, West)
        end
    end
    putmarker!(r)

    k = 0
    directions = (Ost, West)

    while !isborder(r, Nord)
        goToDirect(r, directions[k%2+1])
        move!(r, Nord)
        putmarker!(r)
        k += 1
    end
    goToDirect(r, directions[k%2+1])

    goToStart(r)

    for _  = 1:top
        move!(r, Nord)
    end
    for _ = 1:right
        move!(r, Ost)
    end
end

max_stars_repo_name: Trilonka/MireaAl

max_issues_repo_name: Trilonka/MireaAl

max_forks_repo_name: Trilonka/MireaAl


text:
const ValidPathInput{TF<:AbstractFloat} = Union{Pair{Symbol,<:AbstractVecOrMat{TF}},
    Vector{<:Pair{Symbol,<:AbstractVecOrMat{TF}}},
    Dict{Symbol,<:AbstractVecOrMat{TF}}}
const PathType{TF<:AbstractFloat} = Union{TF,Pair{Int,Vector{TF}},Pair{Int,Matrix{TF}}}

struct Transition{TF<:AbstractFloat}
    jacs::TotalJacobian{TF}
    exovars::Vector{Symbol}
    unknowns::Vector{Symbol}
    varpaths::Dict{Symbol,PathType{TF}}
    H_U::LU{TF,Matrix{TF}}
    upaths::Vector{TF}
    resids::Vector{TF}
end

function Transition(jacs::TotalJacobian{TF}, exopaths::ValidPathInput,
        initials::Union{ValidPathInput,Nothing}=nothing;
        H_U::Union{LU,Nothing}=nothing) where TF
    nT = jacs.nT
    varvals = jacs.varvals
    exopaths isa Pair && (exopaths = (exopaths,))
    initials !== nothing && initials isa Pair && (initials = (initials,))
    varpaths = Dict{Symbol,PathType{TF}}()
    if initials !== nothing
        for (k, v) in initials
            varpaths[k] = v isa Pair ? v : 0=>v
        end
    end
    exovars = Symbol[]
    for (k, v) in exopaths
        k in jacs.srcs || throw(ArgumentError("$k is not a source variable"))
        push!(exovars, k)
        varpaths[k] = v isa Pair ? v : 0=>v
    end
    upaths = TF[]
    unknowns = collect(setdiff(jacs.srcs, exovars))
    for u in unknowns
        upath = get(varpaths, u, nothing)
        if upath === nothing
            upath = fill(varvals[u], nT)
            varpaths[u] = 0=>upath
            append!(upaths, upath)
        else
            upath isa Pair && (upath = upath[2])
            if upath isa Matrix
                throw(ArgumentError("only scalar variable is supported for the unknowns"))
            elseif upath isa Vector
                l = length(upath)
                if l < nT
                    resize!(upath, nT)
                    upath[l+1:end] .= upath[l]
                end
                append!(upaths, upath)
            end
        end
    end
    nU = length(unknowns)
    ntar = length(jacs.tars)
    if H_U === nothing
        H_U = Matrix(hvcat(((nU for _ in 1:ntar)...,),
            (jacs.totals[v][t] for t in jacs.tars for v in unknowns)...))
        H_U = lu!(H_U)
    end
    for blk in jacs.blks
        # Any output must have not been met before unless provided by initials
        for ov in outputs(blk)
            outpath = get(varpaths, ov, nothing)
            if outpath === nothing
                varpaths[ov] = 0=>fill(varvals[ov], nT)
            end
        end
        for inv in invars(blk)
            n = name(inv)
            inpath = get(varpaths, n, nothing)
            if inpath === nothing
                varpaths[n] = varvals[n]
            else
                s = shift(inv)
                if !iszero(s)
                    if inpath isa Pair
                        i0 = inpath[1]
                        inpath = inpath[2]
                    else
                        i0 = 0
                    end
                    if i0 == s
                        continue
                    elseif s < 0
                        if i0 + s >= 0
                            continue
                        else
                            pushfirst!(inpath, (inpath[1] for i in 1:-i0-s)...)
                            varpaths[n] = -s=>inpath
                        end
                    else
                        ilast = i0 + s + nT
                        if ilast <= length(inpath)
                            continue
                        else
                            push!(inpath, (inpath[end] for i in length(inpath)+1:ilast)...)
                        end
                    end
                end
            end
        end
    end
    resids = Vector{TF}(undef, nT*ntar)
    return Transition(jacs, exovars, unknowns, varpaths, H_U, upaths, resids)
end

function _inputs!(tr::Transition, inputs::AbstractVector)
    nT = tr.jacs.nT
    i0 = 0
    @inbounds for u in tr.unknowns
        path = tr.varpaths[u]
        if path isa Pair
            s, path = (path[1], path[2])
        else
            s = 0
        end
        for i in 1:nT
            path[s+i] = inputs[i0+i]
        end
        i0 += nT
    end
end

function _resids!(resids::AbstractVector, tr::Transition)
    nT = tr.jacs.nT
    varvals = tr.jacs.varvals
    i0 = 0
    @inbounds for tar in tr.jacs.tars
        path = tr.varpaths[tar]
        if path isa Pair
            s, path = (path[1], path[2])
        else
            s = 0
        end
        # It is important to subtract the steady-state residuals of the targets
        val = varvals[tar]
        for i in 1:nT
            resids[i0+i] = path[s+i] - val
        end
        i0 += nT
    end
end

function residuals!(resids::AbstractVector, tr::Transition)
    for b in tr.jacs.blks
        transition!(tr.varpaths, b, tr.jacs.nT)
    end
    _resids!(resids, tr)
    return resids
end

residuals!(tr::Transition) = residuals!(tr.resids, tr)

function residuals!(resids::AbstractVector, tr::Transition, inputs::AbstractVector)
    _inputs!(tr, inputs)
    return residuals!(resids, tr)
end

residuals!(tr::Transition, inputs::AbstractVector) =
    residuals!(tr.resids, tr, inputs)

function update!(tr::Transition)
    residuals!(tr, tr.upaths)
    tr.upaths .-= ldiv!(tr.H_U, tr.resids)
    return tr.resids
end

function solve!(tr::Transition; tol::Real=1e-8, maxit::Int=50, pnorm::Real=Inf,
        verbose::Bool=true)
    count = 0
    diffnorm = Inf
    for i in 1:maxit
        count += 1
        diffnorm = norm(update!(tr), pnorm)
        verbose && println("  Iteration: $i   $pnorm-norm: $diffnorm")
        diffnorm < tol && break
    end
    if diffnorm < tol
        verbose && println("Converged in $count iterations!")
    else
        @warn "Convergence failed with $pnorm-norm being $diffnorm"
    end
end


text:
if Sys.isapple()
    function _font_paths()
        [
            "/Library/Fonts", # Additional fonts that can be used by all users. This is generally where fonts go if they are to be used by other applications.
            joinpath(homedir(), "Library/Fonts"), # Fonts specific to each user.
            "/Network/Library/Fonts", # Fonts shared for users on a network
        ]
    end
elseif Sys.iswindows()
    _font_paths() = [joinpath(get(ENV, "SYSTEMROOT", "C:\\Windows"), "Fonts")]
else
    function add_recursive(result, path)
        for p in readdir(path)
            pabs = joinpath(path, p)
            if isdir(pabs)
                push!(result, pabs)
                add_recursive(result, pabs)
            end
        end
    end
    function _font_paths()
        result = String[]
        for p in ("/usr/share/fonts", joinpath(homedir(), "/.fonts"), "/usr/local/share/fonts",)
            if isdir(p)
                add_recursive(result, p)
            end
        end
        result
    end
end

freetype_extensions() = (".FON", ".OTC", ".FNT", ".BDF", ".PFR", ".OTF", ".TTF", ".TTC", ".CFF", ".WOFF")
function freetype_can_read(font::String)
    fontname, ext = splitext(font)
    uppercase(ext) in freetype_extensions()
end

function loaded_faces()
    if isempty(loaded_fonts)
        for path in fontpaths()
            for font in readdir(path)
                # There doesn't really seem to be a reliable pattern here.
                # there are fonts that should be supported and dont load
                # and fonts with an extension not on the FreeType website, which
                # load just fine. So we just try catch it!
                #freetype_can_read(font) || continue
                fpath = joinpath(path, font)
                try
                    push!(loaded_fonts, newface(fpath)[1])
                catch
                end
            end
        end
    end
    return loaded_fonts
end

family_name(x::String) = replace(lowercase(x), ' ' => "") # normalize

function family_name(x)
    fname = x.family_name
    fname == C_NULL && return ""
    family_name(unsafe_string(fname))
end

function style_name(x)
    sname = x.style_name
    sname == C_NULL && return ""
    lowercase(unsafe_string(sname))
end

function match_font(face, name, italic, bold)
    ft_rect = unsafe_load(face)
    fname = family_name(ft_rect)
    sname = style_name(ft_rect)
    italic = italic == (sname == "italic")
    bold = bold == (sname == "bold")
    perfect_match = (fname == name) && italic && bold
    fuzzy_match = occursin(name, fname)
    score = fuzzy_match + bold + italic
    return perfect_match, fuzzy_match, score
end

function try_load(fpath)
    try
        newface(fpath)[]
    catch e
        return nothing
    end
end

function findfont(
        name::String;
        italic = false, bold = false, additional_fonts::String = ""
    )
    font_folders = copy(fontpaths())
    normalized_name = family_name(name)
    isempty(additional_fonts) || pushfirst!(font_folders, additional_fonts)
    candidates = Pair{Ptr{FreeType.FT_FaceRec}, Int}[]
    for folder in font_folders
        for font in readdir(folder)
            fpath = joinpath(folder, font)
            face = try_load(fpath)
            face === nothing && continue
            perfect_match, fuzzy_match, score = match_font(
                face, normalized_name, italic, bold
            )
            perfect_match && return face
            if fuzzy_match
                push!(candidates, face => score)
            else
                FT_Done_Face(face)
            end
        end
    end
    if !isempty(candidates)
        sort!(candidates, by = last)
        final_candidate = pop!(candidates)
        foreach(FT_Done_Face, candidates)
        return final_candidate
    end
    return nothing
end


text:
function build_model(name::String, A)

  num_fifos = 2

  f = open(name * ".json")
  s = JSON.parse(f)

  arch = TopLevel{A,2}(s["name"])
  row = s["dimensions"]["row"]
  col = s["dimensions"]["column"]
  lev = s["dimensions"]["level"]
  proc_add = s["processor_tile"]["add"]
  proc_remove = s["processor_tile"]["remove"]
  mem_proc_add = s["memory_processor_tile"]["add"]
  mem_proc_remove = s["memory_processor_tile"]["remove"]
  one_port = s["memory"]["1_port"]
  two_port = s["memory"]["2_port"]
  num_links = s["configurations"]["num_links"]
  in_handler = s["input_handler"]
  out_handler = s["output_handler"]

  close(f)

  processor = build_processor_tile(num_links)
  if lev > 0
    for r in 0:row-1, c in 0:col-1, l in 0:lev-1
      (in("($r,$c,$l)",proc_remove) || in("($r,$c,$l)",one_port)
       || in("($r,$c,$l)",two_port) || in("($r,$c,$l)",mem_proc_add)
       || in("($r,$c,$l)",mem_proc_remove)) && continue
      add_child(arch, processor, CartesianIndex(r,c,l))
    end
    for proc in proc_add
      r, c, l = stringtoindex3(proc)
      add_child(arch, processor, CartesianIndex(r,c,l))
    end
  elseif lev == 0
    for r in 0:row-1, c in 0:col-1
      (in("($r,$c)",proc_remove) || in("($r,$c)",one_port)
       || in("($r,$c)",two_port) || in("($r,$c)",mem_proc_add)
       || in("($r,$c)",mem_proc_remove)) && continue
      add_child(arch, processor, CartesianIndex(r,c))
    end
    for proc in proc_add
      r, c = stringtoindex2(proc)
      add_child(arch, processor, CartesianIndex(r,c))
    end
  end

  memory_processor = build_processor_tile(num_links, include_memory = true)
  if lev > 0
    for proc in mem_proc_add
      r, c, l = stringtoindex3(proc)
      add_child(arch, memory_processor, CartesianIndex(r,c,l))
    end
  elseif lev == 0
    for proc in mem_proc_add
      r, c = stringtoindex2(proc)
      add_child(arch, memory_processor, CartesianIndex(r,c))
    end
  end

  memory_2port = build_memory(2)
  if lev > 0
    for mem in two_port
      r, c, l = stringtoindex3(mem)
      add_child(arch, memory_2port, CartesianIndex(r,c,l))
    end
  elseif lev == 0
    for mem in two_port
      r, c = stringtoindex2(mem)
      add_child(arch, memory_2port, CartesianIndex(r,c))
    end
  end

  memory_1port = build_memory(2)
  if lev > 0
    for mem in one_port
      r, c, l = stringtoindex3(mem)
      add_child(arch, memory_1port, CartesianIndex(r,c,l))
    end
  elseif lev == 0
    for mem in one_port
      r, c = stringtoindex2(mem)
      add_child(arch, memory_1port, CartesianIndex(r,c))
    end
  end

  input_handler = build_input_handler(1)
  if lev > 0
    for handler in in_handler
      r,c,l = stringtoindex3(handler)
      add_child(arch, input_handler, CartesianIndex(r,c,l))
    end
  elseif lev == 0
    for handler in in_handler
      r,c = stringtoindex2(handler)
      add_child(arch, input_handler, CartesianIndex(r,c))
    end
  end

  output_handler = build_output_handler(1)
  if lev > 0
    for handler in out_handler
      r,c,l = stringtoindex3(handler)
      add_child(arch, output_handler, CartesianIndex(r,c,l))
    end
  elseif lev == 0
    for handler in out_handler
      r,c = stringtoindex2(handler)
      add_child(arch, output_handler, CartesianIndex(r,c))
    end
  end

  connect_processors(arch, num_links)
  connect_io(arch, num_links)
  connect_memories(arch)
  return arch
end

function stringtoindex2(string)

  row,column = split(string,",")
  row = parse(Int64,split(row,"(")[2])
  column = parse(Int64,split(column,")")[1])

  return row, column
end

function stringtoindex3(string)

  row,column,level = split(string,",")
  row = parse(Int64,split(row,"(")[2])
  column = parse(Int64,column)
  level = parse(Int64,split(column,")")[1])

  return row, column, level
end