text:
using Snopt
using SparseArrays

# ---- unconstrained, no dervatives provided. -----

function matyas(g, df, dg, x, deriv)
    f = 0.26 * (x[1]^2 + x[2]^2) - 0.48 * x[1] * x[2]
    fail = false
    return f, fail
end
x0 = [5.0; 7]
lx = [-10.0; -10]
ux = [10.0; 10]
lg = []
ug = []
rows = []
cols = []
options = Dict(
    "Derivative option" => 0
)
xopt, fopt, info, out = snopta(matyas, x0, lx, ux, lg, ug, rows, cols, options)

# ---- unconstrained, providing gradient -----

function rosenbrock(g, df, dg, x, deriv)
    f = (1 - x[1])^2 + 100*(x[2] - x[1]^2)^2
    fail = false

    if deriv
        df[1] = -2*(1 - x[1]) + 200*(x[2] - x[1]^2)*-2*x[1]
        df[2] = 200*(x[2] - x[1]^2)
    end

    return f, fail
end

x0 = [4.0; 4.0]
lx = [-5.0; -5.0]
ux = [5.0; 5.0]
lg = []
ug = []
rows = []
cols = []

xopt, fopt, info, out = snopta(rosenbrock, x0, lx, ux, lg, ug, rows, cols)

# --------- constrained, not providing derivatives ------

function barnes(g, df, dg, x, deriv)

    a1 = 75.196
    a3 = 0.12694
    a5 = 1.0345e-5
    a7 = 0.030234
    a9 = 3.5256e-5
    a11 = 0.25645
    a13 = 1.3514e-5
    a15 = -5.2375e-6
    a17 = 7.0e-10
    a19 = -1.6638e-6
    a21 = 0.0005
    a2 = -3.8112
    a4 = -2.0567e-3
    a6 = -6.8306
    a8 = -1.28134e-3
    a10 = -2.266e-7
    a12 = -3.4604e-3
    a14 = -28.106
    a16 = -6.3e-8
    a18 = 3.4054e-4
    a20 = -2.8673

    x1 = x[1]
    x2 = x[2]
    y1 = x1*x2
    y2 = y1*x1
    y3 = x2^2
    y4 = x1^2

    # --- function value ---

    f = a1 + a2*x1 + a3*y4 + a4*y4*x1 + a5*y4^2 +
        a6*x2 + a7*y1 + a8*x1*y1 + a9*y1*y4 + a10*y2*y4 +
        a11*y3 + a12*x2*y3 + a13*y3^2 + a14/(x2+1) +
        a15*y3*y4 + a16*y1*y4*x2 + a17*y1*y3*y4 + a18*x1*y3 +
        a19*y1*y3 + a20*exp(a21*y1)
    fail = false

    # --- constraints ---

    g[1] = 1 - y1/700.0
    g[2] = y4/25.0^2 - x2/5.0
    g[3] = (x1/500.0- 0.11) - (x2/50.0-1)^2

    return f, fail

end

x0 = [10.0; 10.0]
lx = [0.0; 0.0]
ux = [65.0; 70.0]
lg = -Inf*ones(3)
ug = zeros(3)
rows = [1, 2, 3, 1, 2, 3]
cols = [1, 1, 1, 2, 2, 2]
options = Dict(
    "Derivative option" => 0
)

xopt, fopt, info, out = snopta(barnes, x0, lx, ux, lg, ug, rows, cols, options)

# ----- constrained, providing derivatives (dense jacobian) ---------

function barnesgrad(g, df, dg, x, deriv)

    a1 = 75.196
    a3 = 0.12694
    a5 = 1.0345e-5
    a7 = 0.030234
    a9 = 3.5256e-5
    a11 = 0.25645
    a13 = 1.3514e-5
    a15 = -5.2375e-6
    a17 = 7.0e-10
    a19 = -1.6638e-6
    a21 = 0.0005
    a2 = -3.8112
    a4 = -2.0567e-3
    a6 = -6.8306
    a8 = -1.28134e-3
    a10 = -2.266e-7
    a12 = -3.4604e-3
    a14 = -28.106
    a16 = -6.3e-8
    a18 = 3.4054e-4
    a20 = -2.8673

    x1 = x[1]
    x2 = x[2]
    y1 = x1*x2
    y2 = y1*x1
    y3 = x2^2
    y4 = x1^2

    # --- function value ---

    f = a1 + a2*x1 + a3*y4 + a4*y4*x1 + a5*y4^2 +
        a6*x2 + a7*y1 + a8*x1*y1 + a9*y1*y4 + a10*y2*y4 +
        a11*y3 + a12*x2*y3 + a13*y3^2 + a14/(x2+1) +
        a15*y3*y4 + a16*y1*y4*x2 + a17*y1*y3*y4 + a18*x1*y3 +
        a19*y1*y3 + a20*exp(a21*y1)
    fail = false

    # --- constraints ---

    g[1] = 1 - y1/700.0
    g[2] = y4/25.0^2 - x2/5.0
    g[3] = (x1/500.0- 0.11) - (x2/50.0-1)^2

    if deriv
        # --- derivatives of f ---

        dy1 = x2
        dy2 = y1 + x1*dy1
        dy4 = 2*x1
        dfdx1 = a2 + a3*dy4 + a4*y4 + a4*x1*dy4 + a5*2*y4*dy4 +
            a7*dy1 + a8*y1 + a8*x1*dy1 + a9*y1*dy4 + a9*y4*dy1 + a10*y2*dy4 + a10*y4*dy2 +
            a15*y3*dy4 + a16*x2*y1*dy4 + a16*x2*y4*dy1 + a17*y3*y1*dy4 + a17*y3*y4*dy1 + a18*y3 +
            a19*y3*dy1 + a20*exp(a21*y1)*a21*dy1

        dy1 = x1
        dy2 = x1*dy1
        dy3 = 2*x2
        dfdx2 = a6 + a7*dy1 + a8*x1*dy1 + a9*y4*dy1 + a10*y4*dy2 +
            a11*dy3 + a12*x2*dy3 + a12*y3 + a13*2*y3*dy3 + a14*-1/(x2+1)^2 +
            a15*y4*dy3 + a16*y4*y1 + a16*y4*x2*dy1 + a17*y4*y1*dy3 + a17*y4*y3*dy1 + a18*x1*dy3 +
            a19*y3*dy1 + a19*y1*dy3 + a20*exp(a21*y1)*a21*dy1

        df[1] = dfdx1
        df[2] = dfdx2

        # --- derivatives of g ---

        dg[1] = -x2/700.0  # 1, 1
        dg[2] = 2*x1/25^2   # 2, 1
        dg[3] = 1.0/500  # 3, 1
        dg[4] = -x1/700.0  # 1, 2
        dg[5] = -1.0/5   # 2, 2
        dg[6] = -2*(x2/50.0-1)/50.0  # 3, 2
    end

    return f, fail
end

x0 = [10.0; 10.0]
lx = [0.0; 0.0]
ux = [65.0; 70.0]
lg = -Inf*ones(3)
ug = zeros(3)
rows = [1, 2, 3, 1, 2, 3]
cols = [1, 1, 1, 2, 2, 2]
options = Dict(
    "Verify level" => 1
)

xopt, fopt, info, out = snopta(barnesgrad, x0, lx, ux, lg, ug, rows, cols, options)

# ----- constrained, providing derivatives (sparse jacobian) ---------

function sparsegrad(g, df, dg, x, deriv)

    f = x[1]^2 - x[2]
    fail = false

    # --- constraints ---

    g[1] = x[2] - 2*x[1]
    g[2] = -x[2]

    # --- derivatives of f ---

    df[1] = 2*x[1]
    df[2] = -1.0

    # --- derivatives of g ---

    dg[1] = -2.0  # 1, 1
    dg[2] = 1.0  # 1, 2
    dg[3] = -1.0   # 2, 2

    return f, fail

end

x0 = [0.0; 0.0]
lx = [-10.0, -10.0]
ux = [10.0, 10.0]
lg = -Inf*ones(2)
ug = zeros(2)
rows = [1, 1, 2]
cols = [1, 2, 2]
options = Dict(
    "Verify level" => 1
)

xopt, fopt, info, out = snopta(sparsegrad, x0, lx, ux, lg, ug, rows, cols, options)

# ------- adding names in output file ------

x0 = [12.0; 10.0]
lx = [0.0; 0.0]
ux = [65.0; 70.0]
lg = -Inf*ones(3)
ug = zeros(3)

rows = [1, 2, 3, 1, 2, 3]
cols = [1, 1, 1, 2, 2, 2]
names = Snopt.Names("howdy", ["x1", "x2"], ["obj", "g1", "g2", "g3"])
xopt, fopt, info, out = snopta(barnesgrad, x0, lx, ux, lg, ug, rows, cols, names=names)

# -------- warm start --------------

x0 = [12.0; 10.0]
lx = [0.0; 0.0]
ux = [65.0; 70.0]
lg = -Inf*ones(3)
ug = zeros(3)
rows = [1, 2, 3, 1, 2, 3]
cols = [1, 1, 1, 2, 2, 2]

# artificially limiting the major iterations so we can restart
options = Dict(
    "Major iterations limit" => 15
)

xopt, fopt, info, out = snopta(barnesgrad, x0, lx, ux, lg, ug, rows, cols, options)
println("major iter = ", out.major_iter)

# warm start from where we stopped
warmstart = out.warm  # a WarmStart type
xopt, fopt, info, out = snopta(barnesgrad, warmstart, lx, ux, lg, ug, rows, cols)
println("major iter = ", out.major_iter)

# -------- linear functions --------------

function lp(g, df, dg, x, deriv)
    f = 0.0  # no nonlinear component
    g .= 0.0
    df .= 0.0
    fail = false

    return f, fail
end

x0 = ones(10)
lx = zeros(10)
ux = 4*ones(10)
lg = [5.0, 7, 1]
ug = [8.0, Inf, 10.0]
rows = []
cols = []

A = [
    0.46 0.54 0.40 0.39 0.49 0.03 0.66 0.26 0.05 0.60  # obj
    0.56 0.84 0.23 0.48 0.05 0.69 0.87 0.85 0.88 0.62  # constraint 1
    0.29 0.98 0.36 0.14 0.26 0.41 0.87 0.97 0.13 0.69  # constraint 2
    0.48 0.55 0.78 0.59 0.79 0.84 0.01 0.77 0.13 0.10  # constraint 3
]

xopt, fopt, info, out = snopta(lp, x0, lx, ux, lg, ug, rows, cols, A=A)


text:
NIFgen_lib_found = Libdl.find_library(libnifgen)
if ~isempty(NIFgen_lib_found)
#const lib = Libdl.dlopen(libnifgen) # Open the library explicitly. not working with dlsym in packages...
else
@warn "NI PXI libraries not found!\nPXI functionality broken\nPlease install NI drivers!"
end
#Libdl.dlclose(lib) # close dll
#sym = Libdl.dlsym(lib, :niScope_init)
#ccall(sym, ...) # Use the pointer `sym` instead of the (symbol, library) tuple
# ccall basic:
# ccall(:func, "dllname", return_type,(input types) ,input arguments)

include("wrapper/IviFgen_common.jl")
include("wrapper/IviFgen_api.jl")
include("wrapper/niFgen_common.jl")
include("wrapper/niFgen_api.jl")

ext: jl


text:
#=

    This is a port of the SPPL model hierarchical-markov-switching.pynb

    Note: I'm not sure I understand the real point with the SPPL model, since it
          includes has the p_transition and the two mu_x and mu_y in the model
          when predicting the values. The model has very nice images though. :-)
    The correspondance of the original zs and the generated Z is 92%.

    In this Turing model we get 81% using PG(5):

    As commented below, the SPPL model use the Bernoulli values (0 and 1)
    as indices in the matrices which - of course - don't work in Julia. I first tried
    to use Bernoulli(p) + 1 etc but finally settled with a more traditional if else
    approach.

    That said, the Z values that's created are mostly very close to 1 or 2 (the states), e.g
       [1.0 1.0 1.996 1.996 1.996 1.0 1.996 1.996 1.996 1.996 1.996 1.0 1.0 1.996 1.996
        1.0 1.0 1.0 1.996 1.996 1.996 1.996 1.996 1.996 1.996 1.996 1.996 1.996 1.996
        1.0 1.0 1.0 1.0 1.0 1.0 1.996 1.996 1.996 1.996 1.996 1.996 1.996 1.996 1.0 1.0
        1.0 1.0 1.0 1.0 1.0 1.996 1.996 1.996 1.996 1.996 1.996 1.996 1.0 1.134 1.0 1.0
        1.0 1.0 1.996 1.996 1.134 1.134 1.996 1.0 1.0 1.0 1.0 1.0 1.017 1.017 1.017 1.379
        1.948 1.996 1.996 1.996 1.778 1.76 1.228 1.0 1.271 1.98 1.985 1.943 1.971 1.995
        1.97 1.994 1.996 1.994 1.941 1.968 1.996 1.995 1.996])

    E.g. Z[2] = 1.0 means that all values are exactly 1.0. This - in turn - means that the
    chain summary seems very strange with ess and rhat values as Nan, since all the values are
    exactly the same (i.e. there is no variance at all).

    This is using Pg(5):
    parameters      mean       std   naive_se      mcse        ess      rhat   ess_per_sec
        Symbol   Float64   Float64    Float64   Float64    Float64   Float64       Float64

       separated    1.0040    0.0632     0.0020    0.0040   167.2813    1.0030        5.5092
            Z[1]    1.0000    0.0000     0.0000    0.0000        NaN       NaN           NaN
            Z[2]    1.0000    0.0000     0.0000    0.0000        NaN       NaN           NaN
            Z[3]    1.9960    0.0632     0.0020    0.0040   167.2813    1.0030        5.5092
            Z[4]    1.9960    0.0632     0.0020    0.0040   167.2813    1.0030        5.5092
            Z[5]    1.9960    0.0632     0.0020    0.0040   167.2813    1.0030        5.5092
            Z[6]    1.0000    0.0000     0.0000    0.0000        NaN       NaN           NaN
            Z[7]    1.9960    0.0632     0.0020    0.0040   167.2813    1.0030        5.5092
            Z[8]    1.9960    0.0632     0.0020    0.0040   167.2813    1.0030        5.5092
        ⋮           ⋮         ⋮         ⋮          ⋮         ⋮          ⋮           ⋮
                                                                            92 rows omitted

    Interestingly, using SMC() we got exactly 92% correspondance between zs and Z (as in the
    SPPL model) but to the price of some quite bad rhat values (and low ess values):

    parameters      mean       std   naive_se      mcse       ess      rhat   ess_per_sec
        Symbol   Float64   Float64    Float64   Float64   Float64   Float64       Float64

       separated    1.0000    0.0000     0.0000    0.0000       NaN       NaN           NaN
            Z[1]    1.0000    0.0000     0.0000    0.0000       NaN       NaN           NaN
            Z[2]    1.0000    0.0000     0.0000    0.0000       NaN       NaN           NaN
            Z[3]    1.6000    0.4901     0.0155    0.1165   12.3842    1.0975        1.5656
            Z[4]    2.0000    0.0000     0.0000    0.0000       NaN       NaN           NaN
            Z[5]    2.0000    0.0000     0.0000    0.0000       NaN       NaN           NaN
            Z[6]    1.6540    0.4759     0.0151    0.1210   15.2745    1.0035        1.9310
            Z[7]    1.9710    0.1679     0.0053    0.0290   34.6756    1.0293        4.3838
            Z[8]    2.0000    0.0000     0.0000    0.0000       NaN       NaN           NaN
        ⋮           ⋮         ⋮         ⋮          ⋮         ⋮         ⋮           ⋮
                                                                            92 rows omitted

    (Perhaps a more interesting application would be to recover the matrices from the data.
     I've tried that but it don't work well...
    )
=#

using Turing, StatsPlots, Distributions, StatsBase
using CSV
include("jl_utils.jl")

"""
`flip2(p)`

As `flip(p)` (`Bernoulli(p)`) but returns 1 for false and 2 for two.
"""
function flip2(p)
    # flip(p) == true ? Dirac(2) : Dirac(1) # This only give 1
    # flip(p) == 1.0 ? Dirac(2) : Dirac(1)  # ibid.
    # Dirac((Bernoulli(p) == 1 ? 1 : 0)+1)  # ibid.
    Dirac(rand(Bernoulli(p)) + 1)
end

@model function hierarchical_markov_switching(X, Y, p_transition,mu_x,mu_y)
    n = length(X)
    Z = tzeros(n)
    # Z = Vector{Real}(undef, n)

    #=
    # Note: The model use separated as an index (in mu_x and mu_y) so we have
    #       to use 1..2 instead of 0..1
    separated ~ Categorical([0.6,0.4]) # false = 1, true = 2

    # Sample initial point.
    Z[1] ~ flip2(0.5)
    # And we have to convert Z[.] to Int64 since it has magically converted to a Float64...
    X[1] ~ Normal(mu_x[separated,Int64(Z[1])])
    Y[1] ~ Poisson(mu_y[separated,Int64(Z[1])])

    for i in 2:n
        Z[i] ~ flip2(p_transition[Int64(Z[i-1])])
        X[i] ~ Normal(mu_x[separated,Int64(Z[i])])
        Y[i] ~ Poisson(mu_y[separated,Int64(Z[i])])
    end
    =#

    #
    # But using separated and Z[.] as indices makes it slow and not good.
    # So I rewrote it like this instead.
    #
    separated ~ flip2(0.4)

    # Sample initial point.
    Z[1] ~ flip2(0.5)  # flip2(p) returns 1 for false and 2 for true
    if Z[1] == 1
        X[1] ~ Normal(mu_x[separated,1])
        Y[1] ~ Poisson(mu_y[separated,1])
    else
        X[1] ~ Normal(mu_x[separated,2])
        Y[1] ~ Poisson(mu_y[separated,2])
    end

    for i in 2:n
        if Z[i-1] == 1
            Z[i] ~ flip2(p_transition[1])
        else
            Z[i] ~ flip2(p_transition[2])
        end
        if Z[i] == 1
            X[i] ~ Normal(mu_x[separated,1])
            Y[i] ~ Poisson(mu_y[separated,1])
        else
            X[i] ~ Normal(mu_x[separated,2])
            Y[i] ~ Poisson(mu_y[separated,2])
        end
    end

end

p_transition = [0.2, 0.8]
mu_x = [5 7
        5 15]
mu_y = [5 8
        3 8]

# From a simulation of the SPPL model
xs = [3.880793246124543, 3.5058061097812354, 7.237774689951984, 6.905669457081878, 7.995126208729154, 4.6898044637528224, 6.263190493579016, 7.278231850683178, 7.852084821649775, 5.855167322724675, 5.720088619624132, 5.391958684807007, 7.169767164696265, 6.9726298890755345, 5.1058192848261035, 7.2162584421734435, 5.978129341306567, 4.632434607148869, 7.350778010868155, 7.163199329872764, 3.312995640232463, 6.984375601895885, 7.724455510712415, 5.347415852836971, 6.9430666068360685, 6.182619378597637, 4.474288262345756, 6.39310952727029, 5.016568585393731, 3.395680248255592, 6.485147290740601, 5.3963458371627215, 4.139882819360646, 4.449201069426469, 3.6512381850843667, 5.116290730592944, 6.937674534921111, 7.0500723954996225, 9.812754469939957, 6.196020900659568, 7.0346040847010425, 6.504549077277077, 7.854425984948853, 7.205493992191783, 6.249403045244166, 4.9176270039175245, 2.4866629773701137, 4.41806051644474, 4.578169073441404, 5.554708298980772, 7.41669328653965, 6.835240203522926, 6.309709638182985, 7.420028428918369, 7.731733741224992, 7.842460992547666, 5.988555563719187, 5.691336019473092, 6.2453880266204695, 6.122901245600784, 4.809476409451023, 4.155260934607323, 5.614429568237265, 6.195549626249333, 5.640085080509843, 6.436848039947654, 7.710450159939799, 7.635502275425415, 3.854963010979329, 3.6886646377960455, 4.584214524115992, 5.024753832591269, 3.2863538125483407, 6.787303600428101, 4.770602129374909, 5.706485809946122, 5.555718053917937, 8.060446637998837, 7.737160780287548, 8.301669108932963, 6.111263452419774, 6.36317050704917, 6.831994416446748, 5.245384044497105, 4.513669998612345, 5.847666590699994, 8.009842936472989, 7.511212171510765, 6.392442074959624, 7.417459874615795, 8.44176627215282, 4.99804279054474, 7.479966948085545, 7.666437579480905, 6.973682646784251, 6.613120207032197, 6.014530091460834, 7.8888624942168715, 6.469044373467464, 8.529156359409793]
ys = [5.0, 6.0, 4.0, 11.0, 9.0, 7.0, 11.0, 6.0, 6.0, 9.0, 3.0, 2.0, 8.0, 10.0, 3.0, 10.0, 10.0, 7.0, 6.0, 10.0, 3.0, 8.0, 12.0, 9.0, 9.0, 8.0, 6.0, 9.0, 7.0, 8.0, 5.0, 4.0, 7.0, 7.0, 3.0, 9.0, 6.0, 13.0, 7.0, 14.0, 7.0, 8.0, 11.0, 7.0, 4.0, 6.0, 7.0, 6.0, 3.0, 5.0, 10.0, 9.0, 4.0, 9.0, 5.0, 8.0, 9.0, 2.0, 7.0, 2.0, 4.0, 2.0, 6.0, 6.0, 6.0, 4.0, 5.0, 6.0, 3.0, 3.0, 4.0, 3.0, 3.0, 5.0, 4.0, 8.0, 4.0, 9.0, 5.0, 7.0, 8.0, 5.0, 7.0, 8.0, 2.0, 5.0, 7.0, 7.0, 5.0, 4.0, 11.0, 11.0, 7.0, 8.0, 13.0, 4.0, 9.0, 7.0, 12.0, 9.0]
# zs are the states. SPPL (Python) use 0..1 but we use 1..2 instead
zs = 1 .+ [0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 0, 0, 1, 1, 0, 1, 1, 0, 1, 1, 0, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 1, 1, 1, 1, 1, 0, 0, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1]
# n = length(xs)
n = length(xs)
println("p_transition: $p_transition")
println("mu_x: $mu_x")
println("mu_y: $mu_y")
model = hierarchical_markov_switching(xs,ys,p_transition,mu_x,mu_y)
# model = hierarchical_markov_switching_recover_matrices(xs,ys,zs,p_transition,mu_x,mu_y)

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

# chns = sample(model, HMC(0.1,6), 1_000)
# chns = sample(model, NUTS(), 1_000)

display(chns)

# The given zs
println("zs: $zs")

println("Z before rounding: $(mean(group(chns,:Z).value.data,dims=1)))")
zz = vcat(round.(Int64,mean(group(chns,:Z).value.data,dims=1))...)
println(" Z: $zz")
println("p(zs == Z): $(sum(zz .== zs )/n)")
# for i in 1:n
#     println("i:$i  zs:$(zs[i])  Z[i]:$(zz[i])  eq:$(zs[i] == zz[i])")
# end

# There are quite few that are in the "undecided zone"
println("Elements that differs with 0.25 from 1.5 ('the undecided zone')")
println(filter(i-> abs(i-1.5) < 0.25,   mean(group(chns,:Z).value.data,dims=1)))

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
"""Module for the modeling of dispersion of pollutants in bodies of water."""
module Lakes

using ArchGDAL
using ArgParse
using Dates

include(".\\Utils\\Functions.jl")

export run_lake

const agd = ArchGDAL

mutable struct Lake <: Functions.AbstractAnalysisObject
 # Parameters
    concentration::Float64    # Initial concentration of substance per depth of water.
    time::Int64               # Time from the begin of the diffusion.
    x::Float64                # Distance from secondary source in the X direction.
    y::Float64                # Distance from secondary source in the Y direction.
    fickian_x::Float64        # X direction Fickian transport coefficient.
    fickian_y::Float64        # Y direction Fickian transport coefficient.
    velocity_x::Float64       # Mean velocity of river in the X direction.
    velocity_y::Float64       # Mean velocity of river in the Y direction.
    direction::Int64          # Direction fo the flow.
    λk::Float64               # First order decadiment.
 # Computational results
    C::Float64                # Resulting concentration.

    Lake(concentration,time,x,y,fickian_x,fickian_y,velocity_x,velocity_y,direction,λk) = new(concentration,time,x,y,fickian_x,fickian_y,velocity_x,velocity_y,direction,λk,0.0)
end

function Functions.compute_concentration!( l::Lake )
    c1 =(l.x - (l.velocity_x * l.time))^2.0 / (4.0 * l.fickian_x * l.time)
    c2 =(l.y - (l.velocity_y * l.time))^2.0 / (4.0 * l.fickian_y * l.time)
    c3 = ℯ^(-(c1 + c2)) * ℯ^(-l.λk * l.time)
    c4 = l.concentration / ( 4.0 * π * l.time * √(l.fickian_x * l.fickian_y) )
    l.C = c3 * c4
    return l.C
end

"""
    run_lake(; dem_file::String, source_file::String, area_file::String="", wind_direction::Int64, contaminant_mass::Float64, tollerance::Int64=2,
               mean_flow_speed::Float64, resolution::Int64, hours::Float64, fickian_x::Float64=0.05, fickian_y::Float64=0.05, λk::Float64=0.0,
               output_path::String=".\\lake_otput_model.tiff" )

Create and save as `output_path` a raster containing the results of model of dispersion of pollutants in a lake.

#Arguments
- `dem_file::String`: path to the raster containing the height of the terrain in each cell.
- `source_file::String`: path to the shapefile containing the source point of the contaminants.
- `area_file::String=""`: path to the shapefile containing the poligon delimiting the area for the analysis.
- `wind_direction::Int64`: direction of the wind as an angle in degrees.
- `contaminant_mass::Float64`: initial mass of contaminants.
- `tollerance::Int64=2`: value used to determine wether the concentration of pollutant in a cell is relevant.
    Specifically, a concentration value is considered relevant if its value is within "tollerance" orders of magnitute from the concentration on other cells.
- `mean_flow_speed::Float64`:  mean flow speed of the water.
- `resolution::Float64`: size of the cell for the analysis.
- `hours::Int64`: time span of the analysis in hours.
- `fickian_x::Float64=0.05`: X direction Fickian transport coefficient.
- `fickian_y::Float64=0.05`: Y direction Fickian transport coefficient.
- `λk::::Float64=0.0`: First order decadiment.
- `output_path::AbstractString=".\\lake_otput_model.tiff"`: output file path.
"""#=
function run_lake(; dem_file::String, source_file::String, area_file::String="", wind_direction::Int64, contaminant_mass::Float64, tollerance::Int64=2,
                    mean_flow_speed::Float64, resolution::Float64, hours::Float64, fickian_x::Float64=0.05, fickian_y::Float64=0.05, λk::Float64=0.0,
                    output_path::String=".\\lake_otput_model.tiff" )

    geom = agd.getgeom(collect(agd.getlayer(agd.read(source_file), 0))[1])

    if agd.geomdim(geom) != 0
        throw(DomainError(source_file, "The shapefile must contain a point."))
    end

    dem = agd.read(dem_file)
    refsys = agd.getproj(dem)
    geotransform = agd.getgeotransform(dem)

    if refsys != agd.toWKT(agd.getspatialref(geom))
        throw(DomainError("The reference systems are not uniform. Aborting analysis."))
    end

    if !isempty(area_file)
        trg_geom = agd.getgeom(collect(agd.getlayer(agd.read(area_file), 0))[1])

        agd.geomdim(trg_geom) != 2 && throw(DomainError(area_file, "`area` must be a polygon."))
        agd.toWKT(agd.getspatialref(trg_geom)) != refsys && throw(DomainError("The reference systems are not uniform."))
        !agd.contains(trg_geom, src_geom) && throw(DomainError(area_file, "`area` polygon must contain the source."))
    end

    x_source = agd.getx(geom, 0)
    y_source = agd.gety(geom, 0)
    r_source, c_source = Functions.toIndexes(dem, x_source, y_source)

    velocity_x = √( round( mean_flow_speed * cos(deg2rad(wind_direction)), digits=3 )^2 )
    velocity_y = √( round( mean_flow_speed * sin(deg2rad(wind_direction)), digits=3 )^2 )

    lake = Lake(contaminant_mass, hours*3600.0, 0.0, 0.0, fickian_x, fickian_y, velocity_x, velocity_y, wind_direction, λk)
    points = !isempty(area_file) ? Functions.expand(r_source, c_source, contaminant_mass, tollerance, dem, trg_geom, lake) :
        Functions.expand(r_source, c_source, contaminant_mass, tollerance, dem, lake)

    maxR = maximum( point -> point[1], points )
    minR = minimum( point -> point[1], points )
    maxC = maximum( point -> point[2], points )
    minC = minimum( point -> point[2], points )

    geotransform = agd.getgeotransform(dem)
    geotransform[[1, 4]] = Functions.toCoords(dem, minR, minC)
    noData = Float32(agd.getnodatavalue(agd.getband(dem, 1)))
    data = fill(noData, maxR-minR+1, maxC-minC+1)
    for r in minR:maxR, c in minC:maxC
        match = findfirst( p -> p[1] == r && p[2] == c, points )
        if !isnothing(match)
            data[r-minR+1, c-minC+1] = points[match][3]
        end
    end
    Functions.writeRaster(data, agd.getdriver("GTiff"), geotransform, refsys, noData, output_path)
end
=#
function run_lake(; dem_file::String, source_file::String, area_file::String="", wind_direction::Int64, contaminant_mass::Float64, tollerance::Int64=2,
                    mean_flow_speed::Float64, resolution::Float64, hours::Float64, fickian_x::Float64=0.05, fickian_y::Float64=0.05, λk::Float64=0.0,
                    output_path::String=".\\lake_otput_model.tiff" )

    geom = agd.getgeom(collect(agd.getlayer(agd.read(source_file), 0))[1])

    if agd.geomdim(geom) != 0
        throw(DomainError(source_file, "The shapefile must contain a point."))
    end

    dem = agd.read(dem_file)
    refsys = agd.getproj(dem)
    geotransform = agd.getgeotransform(dem)

    if refsys != agd.toWKT(agd.getspatialref(geom))
        throw(DomainError("The reference systems are not uniform. Aborting analysis."))
    end

    if !isempty(area_file)
        trg_geom = agd.getgeom(collect(agd.getlayer(agd.read(area_file), 0))[1])

        agd.geomdim(trg_geom) != 2 && throw(DomainError(area_file, "`area` must be a polygon."))
        agd.toWKT(agd.getspatialref(trg_geom)) != refsys && throw(DomainError("The reference systems are not uniform."))
        !agd.contains(trg_geom, src_geom) && throw(DomainError(area_file, "`area` polygon must contain the source."))
    end

    x_source = agd.getx(geom, 0)
    y_source = agd.gety(geom, 0)
    r_source, c_source = Functions.toIndexes(dem, x_source, y_source)

    velocity_x = √( round( mean_flow_speed * cos(deg2rad(wind_direction)), digits=3 )^2 )
    velocity_y = √( round( mean_flow_speed * sin(deg2rad(wind_direction)), digits=3 )^2 )

    lake = Lake(contaminant_mass, hours*3600.0, 0.0, 0.0, fickian_x, fickian_y, velocity_x, velocity_y, wind_direction, λk)
    points = !isempty(area_file) ? Functions.expand(r_source, c_source, contaminant_mass, tollerance, dem, trg_geom, lake) :
        Functions.expand(r_source, c_source, contaminant_mass, tollerance, dem, lake)

    maxR = maximum( point -> point[1], points )
    minR = minimum( point -> point[1], points )
    maxC = maximum( point -> point[2], points )
    minC = minimum( point -> point[2], points )

    geotransform[[1, 4]] = Functions.toCoords(dem, minR, minC)
    noData = Float32(agd.getnodatavalue(agd.getband(dem, 1)))
    data = fill(noData, maxR-minR+1, maxC-minC+1)
    for r in minR:maxR, c in minC:maxC
        match = findfirst( p -> p[1] == r && p[2] == c, points )
        if !isnothing(match)
            data[r-minR+1, c-minC+1] = points[match][3]
        end
    end
    Functions.writeRaster(data, agd.getdriver("GTiff"), geotransform, refsys, noData, output_path)
end

end # module

max_stars_repo_name: DavideFavaro/Tirocinio

max_issues_repo_name: DavideFavaro/Tirocinio

max_forks_repo_name: DavideFavaro/Tirocinio


text:
using Catalyst
using Test
using DelaySSAToolkit
rn = @reaction_network begin
    ρ, S+I --> E+I
    1/I, I --> R
    r, I --> R
end ρ r

jumpsys = convert(JumpSystem, rn, combinatoric_ratelaws=false)

u0 = [999,1,0,0] # S, I, E, R
tf = 400.
tspan = (0,tf)
ps = [1e-4, 1e-2] # parameters for ρ, r
τ = 20.
dprob = DiscreteProblem(jumpsys,u0,tspan,ps)

delay_trigger_affect! = function (integrator, rng)
    append!(integrator.de_chan[1], τ)
end
delay_trigger = Dict(1=>delay_trigger_affect!)
delay_interrupt = Dict()
delay_complete = Dict(1=>[2=>1, 3=>-1])
delayjumpset = DelayJumpSet(delay_trigger, delay_complete, delay_interrupt)

de_chan0 = [[]]

meanlist=[]
medianlist=[]
minlist=[]
maxlist=[]
stdlist=[]

for algo in algo_list
    djprob = DelayJumpProblem(jumpsys,dprob,algo,delayjumpset, de_chan0, save_positions=(false,false))
    a=@benchmark solve(djprob, SSAStepper())
    push!(meanlist,copy(mean(a).time/1e9))
    push!(medianlist,copy(median(a).time/1e9))
    push!(minlist,copy(minimum(a).time/1e9))
    push!(maxlist,copy(maximum(a).time/1e9))
    push!(stdlist,copy(std(a).time/1e9))
    print(djprob.aggregator," OK","\n")
end

df=DataFrame(mean=meanlist,median=medianlist,min=minlist,max=maxlist,std=stdlist)
CSV.write("test/zxytest/results/SEIR_save_positions=f.csv",df)

df=CSV.read("test/zxytest/results/SEIR_save_positions=f.csv",DataFrame)
medianlist=df.median

medianvalue=[string(round(mt,digits=6),"s") for mt in medianlist]

algo_name = ["DelayDirect","DelayRejection","DelayMNRM","DelayDirectCR"]
using Plots
p1=bar(algo_name,medianlist,legend=:false,title="SEIR",ylabel="mediantime")
scatter!(algo_name, 0.0002 .+ medianlist , markeralpha=0, series_annotations=medianvalue)
savefig(p1,"test/zxytest/results/SEIR_median.png")

#=
function runbenchmark(algo)
    djprob = DelayJumpProblem(dprob, algo, jumpset, delayjumpset, de_chan0, save_positions = (true,true),save_delay_channel = true)
    a=@benchmark solve(djprob, SSAStepper())
    #a=@bprofile solve(jprob, SSAStepper(), saveat = timestamps)
    return a
end

function testalgo(algo_list,times)
    for algo in algo_list
        onealgotext=zeros(1,times)
        for i=1:times
            bench=runbenchmark(algo)
            onealgotext[i]=median(bench).time/1e3
        end
        push!(test_all_algo,onealgotext)
    end
end

times=5
algo_list = [DelayMNRM(), DelayRejection(), DelayDirectCR()]

test_all_algo=[]
@time testalgo(algo_list,times)

test_all_algo
length(test_all_algo[1])
#mean and std
#DelayMNRM()
mean(test_all_algo[1])
std(test_all_algo[1])
#DelayRejection()
mean(test_all_algo[2])
std(test_all_algo[2])
#DelayDirectCR()
mean(test_all_algo[3])
std(test_all_algo[3])

algo_name = ["MNRM", "Rejection", "DirectCR"]
meanlist=[mean(test_all_algo[i]) for i=1:length(test_all_algo)]
stdlist=[std(test_all_algo[i]) for i=1:length(test_all_algo)]

using Printf
strmean = [@sprintf("%.3fus", yi) for yi in meanlist]
strstd = [@sprintf("%.3f", yi) for yi in stdlist]

using Plots
p1=bar(algo_name,meanlist,ylims=(16,22),text=strmean,bar_width=0.1,
       label="mean of time",text_position="outside")
p2=bar(algo_name,stdlist,ylims=(0.15,0.5),text=strstd,bar_width=0.1,
       label="std of time")
plot(p1, p2, layout = (2,1))
=#


text:
"""
    whoswork(blk, comm)

To which rank is block `blk` currently assigned to.

"""
function whoswork(blk, comm::MPI.Comm)
    blk % MPI.Comm_size(comm)
end

function whoswork(blk, comm::Nothing)
    error("Communicator Nothing should not ask this.")
    return nothing
end

"""
    is_my_work(blk, comm)

Returns a boolean whether the block `blk` is assigned to this rank.

"""
function is_my_work(blk, comm::MPI.Comm)
    whoswork(blk, comm) == MPI.Comm_rank(comm)
end

function is_my_work(blk, comm::Nothing)
    return true
end

abstract type AbstractCommPattern end
# Symmetric communication patterns
struct CommPatternTK <: AbstractCommPattern end
struct CommPatternT <: AbstractCommPattern end
struct CommPatternK <: AbstractCommPattern end

"""
    is_comm_pattern(t, tn, k, kn, pattern)

    Do period t, tn and contingencies k, kn match the communication pattern?

"""
function is_comm_pattern(t, tn, k, kn, ::CommPatternTK)
    return (
        # Neighboring periods and base case (k == 1)
        ((tn == t-1 || tn == t+1) && kn == 1 && k == 1) ||
        # From base case to contingencies
        (tn == t && k == 1 && kn != 1) ||
        # From contingencies to base case
        (tn == t && k != 1 && kn == 1)
    )
end

function is_comm_pattern(t, tn, k, kn, ::CommPatternT)
    return (
        # Neighboring periods and base case (k == 1)
        ((tn == t-1 || tn == t+1) && kn == 1 && k == 1)
    )
end

function is_comm_pattern(t, tn, k, kn, ::CommPatternK)
    return (
        # From base case to contingencies
        (tn == t && k == 1 && kn != 1) ||
        # From contingencies to base case
        (tn == t && k != 1 && kn == 1)
    )
end

"""
    comm_neighbors!(data, blocks, runinfo, pattern, comm)

Nonblocking communication with a given pattern. An array of requests is returned.

"""
function comm_neighbors!(data::AbstractArray{T,2}, blocks::AbstractBlocks, runinfo::ProxALProblem, pattern::AbstractCommPattern, comm::MPI.Comm) where {T}
	requests = MPI.Request[]
    for blk in runinfo.blkLinIndex
        block = blocks.blkIndex[blk]
        k = block[1]
        t = block[2]
        if is_my_work(blk, comm)
            for blkn in runinfo.blkLinIndex
                blockn = blocks.blkIndex[blkn]
                kn = blockn[1]
                tn = blockn[2]
                if is_comm_pattern(t, tn, k, kn, pattern) && !is_my_work(blkn, comm)
                    remote = whoswork(blkn, comm)
                    if isa(pattern, CommPatternTK)
                        sbuf = @view data[:,blk]
                        rbuf = @view data[:,blkn]
                    elseif isa(pattern, CommPatternT)
                        sbuf = @view data[:,t]
                        rbuf = @view data[:,tn]
                    else
                        error("Invalid communication pattern")
                    end
                    push!(requests, MPI.Isend(sbuf, remote, t, comm))
                    push!(requests, MPI.Irecv!(rbuf, remote, tn, comm))
                end
            end
        end
    end
    return requests
end

function comm_neighbors!(data::AbstractArray{T,3}, blocks::AbstractBlocks, runinfo::ProxALProblem, pattern::AbstractCommPattern, comm::MPI.Comm) where {T}
	requests = MPI.Request[]
    for blk in runinfo.blkLinIndex
        block = blocks.blkIndex[blk]
        k = block[1]
        t = block[2]
        if is_my_work(blk, comm)
            for blkn in runinfo.blkLinIndex
                blockn = blocks.blkIndex[blkn]
                kn = blockn[1]
                tn = blockn[2]
                if is_comm_pattern(t, tn, k, kn, pattern) && !is_my_work(blkn, comm)
                    remote = whoswork(blkn, comm)
                    if isa(pattern, CommPatternK)
                        sbuf = @view data[:,k,t]
                        rbuf = @view data[:,kn,tn]
                    else
                        error("Invalid communication pattern")
                    end
                    push!(requests, MPI.Isend(sbuf, remote, k, comm))
                    push!(requests, MPI.Irecv!(rbuf, remote, kn, comm))
                end
            end
        end
    end
    return requests
end

function comm_neighbors!(data::AbstractArray, blocks::AbstractBlocks, runinfo::ProxALProblem, pattern::AbstractCommPattern, comm::Nothing)
    return nothing
end

"""
    comm_wait!(requests)

Wait until the communciation requests `requests` have been fulfilled.

"""
function comm_wait!(requests::Vector{MPI.Request})
    return MPI.Waitall!(requests)
end

function comm_wait!(requests::Nothing)
    return nothing
end

"""
    comm_max(data, comm)

Collective to reduce and return the maximum over scalar `data`.

"""
function comm_max(data::Float64, comm::MPI.Comm)
    return MPI.Allreduce(data, MPI.MAX, comm)
end

function comm_max(data::Float64, comm::Nothing)
    return data
end

"""
    comm_sum(data::Float64, comm)

Collective to reduce and return the sum over scalar `data`.

"""
function comm_sum(data::Float64, comm::MPI.Comm)
    return MPI.Allreduce(data, MPI.SUM, comm)
end

function comm_sum(data::Float64, comm::Nothing)
    return data
end

"""
    comm_sum!(data, comm)

Collective to reduce the sum over array `data`.

"""
function comm_sum!(data::AbstractArray, comm::MPI.Comm)
    return MPI.Allreduce!(data, MPI.SUM, comm)
end

function comm_sum!(data::AbstractArray, comm::Nothing)
    return data
end

function comm_rank(comm::MPI.Comm)
    return MPI.Comm_rank(comm)
end

function comm_rank(comm::Nothing)
    return 0
end

function comm_barrier(comm::MPI.Comm)
    return MPI.Barrier(comm)
end

function comm_barrier(comm::Nothing)
    return nothing
end

max_stars_repo_name: exanauts/proxALM

max_issues_repo_name: exanauts/proxALM

max_forks_repo_name: exanauts/proxALM


text:
import Base:length, size, getindex, get, setindex!, push!, empty!
import StatsBase:sample
export capacity

"""
    SumTree(capacity::Int)

Efficiently sample and update weights.

For more detals, see the post at [here](https://jaromiru.com/2016/11/07/lets-make-a-dqn-double-learning-and-prioritized-experience-replay/).

Here we use a vector to represent the binary tree.
Suppose we will have `capacity` leaves at most.
Every time we `push!` new node into the tree, only the recent `capacity` node and their sum will be updated!

[------------Parent nodes------------][--------leaves--------]
[size: 2^ceil(Int, log2(capacity))-1 ][     size: capacity   ]

# Example
```julia
julia> t = SumTree(8)
0-element SumTree

julia> for i in 1:16
       push!(t, i)
       end

julia> t
8-element SumTree:
  9.0
 10.0
 11.0
 12.0
 13.0
 14.0
 15.0
 16.0

julia> sample(t)
(2, 10.0)

julia> sample(t)
(1, 9.0)

julia> inds, ps = sample(t,100000)
([8, 4, 8, 1, 5, 2, 2, 7, 6, 6  …  1, 1, 7, 1, 6, 1, 5, 7, 2, 7], [16.0, 12.0, 16.0, 9.0, 13.0, 10.0, 10.0, 15.0, 14.0, 14.0  …  9.0, 9.0, 15.0, 9.0, 14.0, 9.0, 13.0, 15.0, 10.0, 15.0])

julia> countmap(inds)
Dict{Int64,Int64} with 8 entries:
  7 => 14991
  4 => 12019
  2 => 10003
  3 => 11027
  5 => 12971
  8 => 16052
  6 => 13952
  1 => 8985

julia> countmap(ps)
Dict{Float64,Int64} with 8 entries:
  9.0  => 8985
  13.0 => 12971
  10.0 => 10003
  14.0 => 13952
  16.0 => 16052
  11.0 => 11027
  15.0 => 14991
  12.0 => 12019
```
"""
mutable struct SumTree <: AbstractArray{Int, 1}
    capacity::Int
    first::Int
    length::Int
    nparents::Int
    tree::Vector{Float64}
    function SumTree(capacity)
        nparents = 2^ceil(Int, log2(capacity)) - 1
        new(capacity, 1, 0, nparents, zeros(nparents + capacity))
    end
end

capacity(t::SumTree) = t.capacity
length(t::SumTree) = t.length
size(t::SumTree) = (length(t),)

function _index(t::SumTree, i::Int)
    ind = i + t.first - 1
    if ind > t.capacity
        ind -= t.capacity
    end
    ind
end

_tree_index(t::SumTree, i) = t.nparents + _index(t, i)

getindex(t::SumTree, i::Int) = t.tree[_tree_index(t, i)]

function setindex!(t::SumTree, p, i)
    tree_ind = _tree_index(t, i)
    change = p - t.tree[tree_ind]
    t.tree[tree_ind] = p
    while tree_ind != 1
        tree_ind = tree_ind ÷ 2
        t.tree[tree_ind] += change
    end
end

function push!(t::SumTree, p)
    if t.length == t.capacity
        t.first = (t.first == t.capacity ? 1 : t.first + 1)
    else
        t.length += 1
    end
    t[t.length] = p
end

function empty!(t::SumTree)
    t.length = 0.
    fill!(t.tree, 0.)
    # yes, no need to reset `t.first`
    # so, don't rely on that `t.first` is always 1 after `empty!`
    t
end

function get(t::SumTree, v)
    parent_ind = 1
    leaf_ind = parent_ind
    while true
        left_child_ind = parent_ind * 2
        right_child_ind = left_child_ind + 1
        if left_child_ind > length(t.tree)
            leaf_ind = parent_ind
            break
        else
            if v ≤ t.tree[left_child_ind]
                parent_ind = left_child_ind
            else
                v -= t.tree[left_child_ind]
                parent_ind = right_child_ind
            end
        end
    end
    p = t.tree[leaf_ind]
    ind = leaf_ind - t.nparents
    real_ind = ind >= t.first ? ind - t.first + 1 : ind + t.capacity - t.first + 1
    real_ind, p
end

sample(t::SumTree) = get(t, rand() * t.tree[1])

function sample(t::SumTree, n::Int)
    inds, priorities = Vector{Int}(undef, n), Vector{Float64}(undef, n)
    for i in 1:n
        ind, p = sample(t)
        inds[i] = ind
        priorities[i] = p
    end
    inds, priorities
end


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

mutable struct NetworkWatcherPropertiesFormat <: SwaggerModel
    provisioningState::Any # spec type: Union{ Nothing, ProvisioningState } # spec name: provisioningState

    function NetworkWatcherPropertiesFormat(;provisioningState=nothing)
        o = new()
        validate_property(NetworkWatcherPropertiesFormat, Symbol("provisioningState"), provisioningState)
        setfield!(o, Symbol("provisioningState"), provisioningState)
        o
    end
end # type NetworkWatcherPropertiesFormat

const _property_map_NetworkWatcherPropertiesFormat = Dict{Symbol,Symbol}(Symbol("provisioningState")=>Symbol("provisioningState"))
const _property_types_NetworkWatcherPropertiesFormat = Dict{Symbol,String}(Symbol("provisioningState")=>"ProvisioningState")
Base.propertynames(::Type{ NetworkWatcherPropertiesFormat }) = collect(keys(_property_map_NetworkWatcherPropertiesFormat))
Swagger.property_type(::Type{ NetworkWatcherPropertiesFormat }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_NetworkWatcherPropertiesFormat[name]))}
Swagger.field_name(::Type{ NetworkWatcherPropertiesFormat }, property_name::Symbol) =  _property_map_NetworkWatcherPropertiesFormat[property_name]

function check_required(o::NetworkWatcherPropertiesFormat)
    true
end

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


text:
"""
utils: a few miscellaneous helper functions
"""

module utils

#Get positions for proteins between i, j.
function define_region(i, j=0)
	if j == 0
		j = i
	end
	return (21 * (i - 1)) + 1 : 21 * j
end

#Lesser than for two inputs (used to sort through edges).
function two_lt(x::Tuple)
	return x[1] < x[2]
end

function big_vec(n)
	outp = Tuple[]
	for i in 1:n
		for j in 1:n
			push!(outp, (i, j))
		end
	end
	return outp
end

function ifilter(f::Function, itr)
	function _it()
		for i = itr
			if f(i)
				produce(i)
			end
		end
	end
	Task(_it)
end

function top_n(el_tuple, val_keys)
  a_max = findmax(el_tuple)
  max_key::Int = val_keys[a_max[2]]
  max_value = a_max[1]
  return (max_key, max_value)
end

function aa_vec_to_seq(vec)
	#Takes in a wide one-hot vector of protein sequence, converts it back to aa sequence.
	aa_seq = Char[]
	prot_len = div(length(vec), 21)
	for i in 1:prot_len
		sub_vec = vec[((i - 1) * 21) + 1 : i * 21]
		sub_aa = findmax(sub_vec)[2]
		push!(aa_seq, "ARNDCQEGHILKMFPSTWYV-"[sub_aa])
	end
	return join(aa_seq)
end

end


text:
abstract type Cartesian end
abstract type Polar end

mutable struct Point2D{T}
    x::Float64
    z::Float64
end

struct Grid{A,B,C,D}
    x::A
    z::A
    θ::A
    r::A
    e2n::B
    e2n_p1::B
    e2nP::B
    nθ::C
    nr::C
    nel::C
    nnod::C
    nVnod::C
    nvert::C
    neighbours::D
    neighbours_p1::D
end

struct ElementCoordinates{T}
    θ::T
    r::T
end

cartesian2polar(x,z) = (atan(x,z), sqrt(x^2+z^2))

polar2cartesian(x,z) = (z*sin(x), z*cos(x))

function init_grid(N; split=1)
    if split == 1
        nr = Int(1 + 2^N)
        nθ = Int(12 * 2^N)
        gr = Grid_split1(nθ, nr)
    else
        nr = Int(1+2^N)
        nθ = Int(12*2^N)
        if N == 5
            nr = Int(1 + 32)
            nθ = Int(256)
        end
        gr = Grid(nθ, nr)
    end
    IDs = point_ids(gr)

    return gr, IDs
end

function Grid(nθ::Int64, nr::Int64; r_out=2.22, r_in=1.22)

    r_out = r_out
    r_in = r_in
    nn = nr*(nθ)# total number of nodes
    nels = (nr-1)*nθ
    r = fill(0.0, nn)
    θ = fill(0.0, nn)
    angle = deg2rad(0)
    # angle   = 0
    dθ  = 2*π/nθ

    # -- Nodal positions
    Ir = fill(1.0, nr)
    @inbounds for ii in 1:nθ
        idx = @. (1:nr) + nr*(ii-1)
        r[idx] = LinRange(r_in, r_out, nr)
        @. (θ[idx] = angle + dθ*(ii-1)*Ir)
    end

    # -- Quadrilateral elements
    id_el = Array{Int64,2}(undef,nels,4)
    @inbounds for ii in 1:nθ
        if ii < nθ
            idx  = @. (1:nr-1) + (nr-1)*(ii-1)
            idx1 = @. (1:nr-1) + nr*(ii-1)
            idx2 = @. (1:nr-1) + nr*(ii)
            id_el[idx,:] .= [idx1 idx2 idx2.+1 idx1.+1]
        else
            idx  = @. (1:nr-1) + (nr-1)*(ii-1)
            idx1 = @. (1:nr-1) + nr*(ii-1)
            idx2 = @. 1:nr-1
            id_el[idx,:] .= [idx1 idx2 idx2.+1 idx1.+1]
        end
    end

    # Split quadrilateral in 2 triangles
    ntri = nels * 4
    id_els_tri = Array{Int64,2}(undef,ntri,3)
    dθ /= 2
    angle = dθ
    dr = r[2]-r[1]
    r_c, θ_c = fill(0.0, nels), fill(0.0, nels)
    id_centers = nn .+ (1 : nels)
    #         id_el       = [id_el id_centers']
    i1 = [1, 2]
    i2 = [2, 3]
    i3 = [3, 4]
    i4 = [4, 1]
    @inbounds for ii in 1 : nels
        # connectivity matrix
        idx = @. (1:4) + 4*(ii-1)
        id_els_tri[idx[1],:] = [id_el[ii,i1]' id_centers[ii]]
        id_els_tri[idx[2],:] = [id_el[ii,i2]' id_centers[ii]]
        id_els_tri[idx[3],:] = [id_el[ii,i3]' id_centers[ii]]
        id_els_tri[idx[4],:] = [id_el[ii,i4]' id_centers[ii]]
        # coords of central node
        θ_c[ii] = θ[id_el[ii,1]] + dθ
        r_c[ii] = r[id_el[ii,1]] + dr/2
    end

    r = vcat(r, r_c)
    θ = vcat(θ, θ_c)
    EL2NOD = Array(id_els_tri')
    neighbours = element_neighbours(EL2NOD)
    θ, r, EL2NOD = add_midpoints(EL2NOD, neighbours, θ, r)
    EL2NOD_P1 = p2top1(EL2NOD)
    neighbours_p1 = element_neighbours(EL2NOD_P1)
    nel = size(EL2NOD,2)
    nVnod = maximum(view(EL2NOD,1:3,:))
    # EL2NOD, θ, r =  sixth_node(EL2NOD, θ, r)
    x, z = polar2cartesian(θ,r)
    nnod = length(x)
    # EL2NODP = Matrix(reshape(1:3*nel,3,nel))
    EL2NODP = EL2NOD[1:3,:]

    Grid(x,
        z,
        θ,
        r,
        EL2NOD,
        EL2NOD_P1,
        EL2NODP,
        nθ,
        nr,
        nel,
        nnod,
        nVnod,
        3,
        neighbours,
        neighbours_p1)
end

function Grid_split1(nθ::Int64, nr::Int64)
        # 1 = split rectangle into 2 triangles
        # 2 = split rectangle into 4 triangles
    R  = 2.22
    r_out = R
    r_in = 1.22
    nn = nr*(nθ)# total number of nodes
    nels = (nr-1)*nθ
    r = fill(0.0, nn)
    θ = fill(0.0, nn)
    angle = deg2rad(0)
    # angle   = 0
    dθ  = 2*π/nθ

    # -- Nodal positions
    Ir = fill(1.0, nr)
    @inbounds for ii in 1:nθ
        idx = @. (1:nr) + nr*(ii-1)
        r[idx] = LinRange(r_in, r_out, nr)
        @. (θ[idx] = angle + dθ*(ii-1)*Ir)
    end

    # -- Quadrilateral elements
    id_el = Array{Int64,2}(undef,nels,4)
    @inbounds for ii in 1:nθ
        if ii < nθ
            idx  = @. (1:nr-1) + (nr-1)*(ii-1)
            idx1 = @. (1:nr-1) + nr*(ii-1)
            idx2 = @. (1:nr-1) + nr*(ii)
            id_el[idx,:] .= [idx1 idx2 idx2.+1 idx1.+1]
        else
            idx  = @. (1:nr-1) + (nr-1)*(ii-1)
            idx1 = @. (1:nr-1) + nr*(ii-1)
            idx2 = @. 1:nr-1
            id_el[idx,:] .= [idx1 idx2 idx2.+1 idx1.+1]
        end
    end

    # -- Split quadrilateral in 2 triangles
    ntri = nels * 2
    id_els_tri = Array{Int64,2}(undef,ntri,3)
    i1 = [1, 2, 4]
    i2 = [2, 3, 4]
    @inbounds for ii in 1 : nels
        # -- connectivity matrix
        idx = @. (1:2) + 2*(ii-1)
        id_els_tri[idx[1],:] .= id_el[ii,i1]
        id_els_tri[idx[2],:] .= id_el[ii,i2]
    end

    EL2NOD = Array(id_els_tri')
    neighbours = element_neighbours(EL2NOD)
    θ, r, EL2NOD = add_midpoints(EL2NOD, neighbours, θ, r)
    EL2NOD_P1 = p2top1(EL2NOD)
    neighbours_p1 = element_neighbours(EL2NOD_P1)
    nel = size(EL2NOD,2)
    nVnod = maximum(view(EL2NOD,1:3,:))
    # EL2NOD, θ, r =  sixth_node(EL2NOD, θ, r)
    x, z = polar2cartesian(θ,r)
    nnod = length(x)
    # EL2NODP = Matrix(reshape(1:3*nel,3,nel))
    EL2NODP = EL2NOD[1:3,:]

    nvert = 3

    Grid(x,
        z,
        θ,
        r,
        EL2NOD,
        EL2NOD_P1,
        EL2NODP,
        nθ,
        nr,
        nel,
        nnod,
        nVnod,
        nvert,
        neighbours_p1)
end

function element_area(gr::Grid)
    xel = view(gr.x, gr.e2n)
    zel = view(gr.z, gr.e2n)
    @. @views 0.5*(xel[1,:]*(zel[2,:] - zel[3,:]) +
                   xel[2,:]*(zel[3,:] - zel[1,:]) +
                   xel[3,:]*(zel[1,:] - zel[2,:]))
end

function element_area(x, z, e2n)
    xel = view(x, e2n)
    zel = view(z, e2n)
    @. @views 0.5*(xel[1,:]*(zel[2,:] - zel[3,:]) +
                   xel[2,:]*(zel[3,:] - zel[1,:]) +
                   xel[3,:]*(zel[1,:] - zel[2,:]))
end

function inradius(gr::Grid)
    # https://en.wikibooks.org/wiki/Trigonometry/Circles_and_Triangles/The_Incircle
    x, z, nel, e2n = gr.x, gr.z, gr.nel, gr.e2n_p1
    r = Vector{Float64}(undef,nel)
    @inbounds @fastmath for iel in 1:nel
        xv = ntuple(i->x[e2n[i,iel]], 3) # verteices x-coords
        zv = ntuple(i->z[e2n[i,iel]], 3) # verteices z-coords
        a = distance(xv[1], zv[1], xv[2], zv[2]) # distance of side AB
        b = distance(xv[1], zv[1], xv[3], zv[3]) # distance of side AC
        c = distance(xv[3], zv[3], xv[2], zv[2]) # distance of side BC
        s = (a+b+c)*0.5 # semiperimeter
        r[iel] = √((s-a)*(s-b)*(s-c)/s) # inradius
    end
    minimum(r)
end

function inrectangle(gr::Grid)
    # define a percentage of the diagonal
    x, z, nel, e2n = gr.x, gr.z, gr.nel*4, gr.e2n_p1
    # x, z, nel, e2n = gr.x, gr.z, gr.nel, gr.e2n
    r = Vector{Float64}(undef,nel)
    @inbounds @fastmath for iel in 1:nel
        xv = ntuple(i->x[e2n[i,iel]], 3) # verteices x-coords
        zv = ntuple(i->z[e2n[i,iel]], 3) # verteices z-coords
        xmin, xmax = extrema(xv)
        zmin, zmax = extrema(zv)
        r[iel] = √( (xmin-xmax)^2 + (zmin-zmax)^2) # inradius
    end

    minimum(r)
end

function sixth_node(EL2NOD, θ, r)
    nel = size(EL2NOD, 2)
    nnod = maximum(EL2NOD)
    θel = θ[EL2NOD[1:3,:]]
    rel = r[EL2NOD[1:3,:]]
    fixangles!(θel)
    θ_new = mean(θel, dims=1)
    r_new = mean(rel, dims=1)
    θ = vcat(θ, vec(θ_new))
    r = vcat(r, vec(r_new))
    return vcat(EL2NOD, transpose(nnod.+(1:nel))), θ, r
end

function p2top1(EL2NOD)
    nel = size(EL2NOD, 2)
    EL2NOD_P1 = Array{Int64,2}(undef,3,nel*4)
    local_map = [1 4 6; 4 2 5; 6 5 3; 4 5 6]'
    @inbounds for iparent = 1:nel
        local_element = view(EL2NOD,:, iparent)
        child_elements = view(local_element,local_map)
        iglobal = (1:4) .+ 4*(iparent-1)
        EL2NOD_P1[:, iglobal] .= child_elements
    end
    return EL2NOD_P1
end

function element_neighbours(e2n)
    els = size(e2n,1) == 3 ? e2n : view(e2n, 1:3,:)
    nel = size(els,2)

    # Incidence matrix
    I, J, V = Int[], Int[], Bool[]
    @inbounds for i in axes(els,1), j in axes(els,2)
        node = els[i,j]
        push!(I, node)
        push!(J, j)
        push!(V, true)
    end
    incidence_matrix = sparse(J, I, V)

    # Find neighbouring elements
    neighbour = [Int64[] for _ in 1:nel]
    nnod = I[end]
    @inbounds for node in 1:nnod
        # Get elements sharing node
        r = nzrange(incidence_matrix, node)
        el_neighbour = incidence_matrix.rowval[r]
        # Add neighbouring elements to neighbours map
        for iel1 in el_neighbour
            current_neighbours = neighbour[iel1]
            for iel2 in el_neighbour
                # check for non-self neighbour and avoid repetitions
                if (iel1!=iel2) && (iel2 ∉ current_neighbours)
                    push!(neighbour[iel1],iel2)
                end
            end
        end
    end

    return neighbour
end

function edge_connectivity(EL2NOD, neighbours)
    nel = size(EL2NOD,2)
    local_map = [ 1 2 3; 2 3 1]
    global_idx = 0
    el2edge = fill(0,3,nel)
    edge2node = Vector{Int64}[]
    edge_neighbour = Array{Int64,2}(undef, 2, 3)

    @inbounds for iel in 1:nel
        # local edge
        edge = EL2NOD[local_map, iel]
        sort!(edge, dims=1)

        # neighbours of local element
        el_neighbours = neighbours[iel]
        edge_neighbours = [EL2NOD[local_map, i] for i in el_neighbours]

        # check edge by edge
        for iedge in 1:3

            if el2edge[iedge, iel] == 0
                global_idx += 1
                el2edge[iedge, iel] = global_idx
                push!(edge2node, view(edge,:, iedge))

                # edges the neighbours
                for (c, ieln) in enumerate(el_neighbours)
                    edge_neighbour .= edge_neighbours[c]
                    sort!(edge_neighbour, dims=1)

                    # check wether local edge is in neighbouring element
                    for i in 1:3
                        if (edge[1, iedge] == edge_neighbour[1,i]) && (edge[2, iedge] ∈  edge_neighbour[2,i])
                            el2edge[i, ieln] = global_idx
                            break
                        end
                    end

                end

            end
        end
    end

    return el2edge, edge2node

end

function add_midpoints(EL2NOD, neighbours, θ, r)
    el2edge, edges2node = edge_connectivity(EL2NOD, neighbours)
    # make mid points
    nnods = size(edges2node, 1)
    θmid = Vector{Float64}(undef,nnods)
    rmid = similar(θmid)
    @inbounds for i in 1:nnods
        θbar = @MVector [θ[edges2node[i][1]], θ[edges2node[i][2]]]
        if abs.(θbar[1] - θbar[2]) > π
            for j in 1:2
                if θbar[j] < π
                    θbar[j] += 2π
                end
            end
        end
        θmid[i] = 0.5*(θbar[1]+θbar[2])
        rmid[i] = 0.5*(r[edges2node[i][1]] + r[edges2node[i][2]])
    end

    old_nnods = maximum(EL2NOD)
    EL2NOD_new = el2edge .+ old_nnods

    return vcat(θ, θmid), vcat(r, rmid), vcat(EL2NOD, EL2NOD_new)
end

function point_ids(M::Grid)

    top = "outter"
    bot = "inner"
    inner = "inside"

    nnod = length(M.r)
    IDs = Vector{String}(undef,nnod)

    rmin, rmax = extrema(M.r)

    @inbounds for (i, ri) in enumerate(M.r)
        if ri == rmax
            IDs[i] = top
        elseif ri == rmin
            IDs[i] = bot
        else
            IDs[i] = inner
        end
    end

    return IDs
end

max_stars_repo_name: albert-de-montserrat/Persephone

max_issues_repo_name: albert-de-montserrat/Persephone

max_forks_repo_name: albert-de-montserrat/Persephone


text:
@testset "Quadratic" begin
    p1 = SVector(1.0, 1.0, 1.0); p2 = SVector(2.0, 2.0, 0.0); p3 = SVector(3.0, 3.5, 1.0);
    q = Quadratic(p1, p2, p3)

    @test q[0] ≈ p1
    @test q[1.405563856998] ≈ p2
    @test q[end] ≈ p3

    @test all(map(isapprox, q[[0, 1.405563856998, length(q)]], [p1, p2, p3]))

    # TODO test transformations
end


text:
using Images
using FileIO

# load our julia function
include("julia.jl")

h = 2500
w = 5000
iter = 10000
m = Array{Float64, 3}(undef, 3, h, w)
m = fill!(m, 0.0)

# time measurements
print("starting...\n")
tStart=time()

# for every pixel
for y=1:h, x=1:w
    # translate numbers [1:w, 1:h] -> -2:2 + -1:1 im
    c = complex((x-w/2)/(h/2), (y-h/2)/(h/2))
    # call our julia function
    color = julia(c, iter)/convert(Float64, iter)

    # color adjust weights
    red_color, blue_color, green_color = 1.0, 1.0, 1.0
    if color < 0.01
	    red_color, blue_color, green_color = 1.0, 1.0, 1.0
    else
	    red_color = color * 0.80
	    blue_color = color * 0.80
	    green_color = color * 0.15
    end

     m[1, y,x]  = red_color
     m[2, y,x]  = blue_color
     m[3, y,x]  = green_color
end

tStop = time()
print("done. took ", tStop-tStart, " seconds\n");

img = colorview(RGB, m)
save("julia.jpg", img)

max_stars_repo_name: rabarar/juliaFractals

max_issues_repo_name: rabarar/juliaFractals

max_forks_repo_name: rabarar/juliaFractals


text:
function tcm_predictForest(N, K, lambda, forest_left, forest_right, forest_product, forest_isLeaf, assortments)
	M = size(assortments,1)
	v = zeros(size(assortments))
	v_by_class = Array{ Array{Float64,2}, 1}(K)
	for t in 1:K
		v_one_class = zeros(Float64,size(assortments))
	    for m = 1 : M
	        cn = 1;
	        while (!forest_isLeaf[t][cn])
	        	if (assortments[m, forest_product[t][cn]] > 0)
	        		cn = forest_left[t][cn]
	        	else
	        		cn = forest_right[t][cn]
	        	end
	        end
	        v_one_class[m, forest_product[t][cn]] = 1.0
	    end
		v += lambda[t] * v_one_class

		v_by_class[t] = v_one_class
	end

	return v, v_by_class
end

max_stars_repo_licenses: MIT Unlicense

max_issues_repo_licenses: MIT Unlicense

max_forks_repo_licenses: MIT Unlicense


text:
function LeastSquaresConverter(arg0::MultivariateVectorFunction, arg1::Vector{jdouble})
    return LeastSquaresConverter((MultivariateVectorFunction, Vector{jdouble}), arg0, arg1)
end

function LeastSquaresConverter(arg0::MultivariateVectorFunction, arg1::Vector{jdouble}, arg2::RealMatrix)
    return LeastSquaresConverter((MultivariateVectorFunction, Vector{jdouble}, RealMatrix), arg0, arg1, arg2)
end

function LeastSquaresConverter(arg0::MultivariateVectorFunction, arg1::Vector{jdouble}, arg2::Vector{jdouble})
    return LeastSquaresConverter((MultivariateVectorFunction, Vector{jdouble}, Vector{jdouble}), arg0, arg1, arg2)
end

function value(obj::LeastSquaresConverter, arg0::Vector{jdouble})
    return jcall(obj, "value", jdouble, (Vector{jdouble},), arg0)
end


text:
function preprocess_datafiles()
  # read unprocessed datafiles
  data = readtable("../data/raw/atussum.csv")

  # define activity-to-columns matching
  # format:
  # symbol for building dataframe, array of symbols corresponding to activity in original data
  headers=string.(names(data))
  activities = ((:sleep,         Symbol.(["t010101"])),
                (:work,          Symbol.(["t050101", "t050102"])),
                (:tv,            Symbol.(["t120305", "t120303", "t120304", "t120308"])),
                (:edu,           Symbol.(filter(x->ismatch(r"t06",x),headers))),
                (:household,     Symbol.(filter(x->ismatch(r"t02",x),headers))),
                (:relax,         Symbol.(filter(x->ismatch(r"t12",x),headers))),
                (:carehh,        Symbol.(filter(x->ismatch(r"t03",x),headers))),
                (:eat,           Symbol.(filter(x->ismatch(r"t11",x),headers))),
                (:travel,        Symbol.(filter(x->ismatch(r"t18",x),headers))),
                (:personalcare,  Symbol.(filter(x->ismatch(r"t01.[^1]",x),headers))),
                (:sport,         Symbol.(filter(x->ismatch(r"t13",x),headers))),
                (:buy,           Symbol.(filter(x->ismatch(r"t07",x),headers))),
                (:religion,      Symbol.(filter(x->ismatch(r"t14",x),headers))),
                (:telephone,     Symbol.(filter(x->ismatch(r"t16",x),headers))),
                (:carenhh,       Symbol.(filter(x->ismatch(r"t04",x),headers))),
                (:volunteer,     Symbol.(filter(x->ismatch(r"t15",x),headers)))
            );
  # build new dataframe
  data2=DataFrame()
  data2[:weight] = data[:tufnwgtp]
  data2[:age] = data[:teage]
  data2[:employment] = data[:trdpftpt]
  for act in activities
    sym = act[1]
    data2[sym] = 0
    for origsym in act[2]
      data2[sym] += data[origsym]
    end
  end

  #write processed data
  mkpath("../data/processed")
  writetable("../data/processed/processed.csv",data2)
end


text:
# WE routines

""" `run_we`: Run a serial WE simulation, optionally returning the ensemble at
each, step with

### Arguments
* `E₀` - initial particle ensemble
* `B₀` - initial bin data structure
* `sampler` - WE sampler functions data structure
* `n_we_steps` - number of steps in the WE run
### Optional Arguments
* `n_save_iters = 1` - save the ensemble and bins every `n_save_iters` iterations.  Set `n_save_iters=n_we_steps` to only save the final values.
"""
function run_we(E₀::TE, B₀::TB, sampler::TWES, n_we_steps::Int; n_save_iters = 1) where
   {TE<:Ensemble, TB<:Bins, TWES<:WEsampler}

   E = deepcopy(E₀);
   B = deepcopy(B₀);
   E_trajectory = TE[];
   B_trajectory = TB[];

   for t in 0:n_we_steps-1
      # first selection is at t = 0
      sampler.selection!(E, B, t);
      copy!(E.ω, E.ω̂);
      @. E.ξ = deepcopy(E.ξ̂);
      @. E.d = deepcopy(E.d̂);
      sampler.mutation!.(E.ξ);
      # after mutation, time is t ↦ t+1
      sampler.rebin!(E, B, t+1);
      # analysis
      sampler.analysis!(E, B, t+1);

      if(mod(t+1, n_save_iters)==0)
         push!(E_trajectory, deepcopy(E))
         push!(B_trajectory, deepcopy(B))
      end
   end

   return E_trajectory, B_trajectory

end

"""
`run_we_observables`: Run a serial WE simulation, returning the values a
specified fucntion, `f`, along the trajecotry.

### Arguments
* `E₀` - initial particle ensemble
* `B₀` - initial bin data structure
* `sampler` - WE sampler functions data structure
* `n_we_steps` - number of steps in the WE run
* `observables` - Tuple of scalar observable functions for the ergodic average
"""
function run_we_observables(E₀::TE, B₀::TB, sampler::TWES, n_we_steps::Int, observables::Tuple{Vararg{<:Function,NO}}) where
   {TE<:Ensemble, TB<:Bins, TWES<:WEsampler, NO}

   E = deepcopy(E₀);
   B = deepcopy(B₀);
   observables_trajectory = zeros(NO, n_we_steps);

   for t in 0:n_we_steps-1
      # first selection is at t = 0
      sampler.selection!(E, B, t);
      copy!(E.ω, E.ω̂);
      @. E.ξ = deepcopy(E.ξ̂);
      @. E.d = deepcopy(E.d̂);
      sampler.mutation!.(E.ξ);
      # after mutation, time is t ↦ t+1
      sampler.rebin!(E, B, t+1);
      # analysis
      sampler.analysis!(E, B, t+1);

      ntuple(k-> observables_trajectory[k,t+1] =(observables[k]).(E.ξ) ⋅ E.ω, NO)
   end

   return observables_trajectory

end

"""
`run_we!`: Run an in place serial WE simulation with

### Arguments
* `E` - particle ensemble
* `B` - bin data structure
* `sampler` - WE sampler functions data structure
* `n_we_steps` - number of steps in the WE run
"""
function run_we!(E::TE, B::TB, sampler::TWES, n_we_steps::Int) where
   {TE<:Ensemble, TB<:Bins, TWES<:WEsampler}

   for t in 0:n_we_steps-1
      # first selection is at t = 0
      sampler.selection!(E, B, t);
      copy!(E.ω, E.ω̂);
      @. E.ξ = deepcopy(E.ξ̂);
      @. E.d = deepcopy(E.d̂);
      sampler.mutation!.(E.ξ);
      # after mutation, time is t ↦ t+1
      sampler.rebin!(E, B, t+1);
      # analysis
      sampler.analysis!(E, B, t+1);

   end
   E, B

end


text:
"""
GenX: An Configurable Capacity Expansion Model
Copyright (C) 2021,  Massachusetts Institute of Technology
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.
A complete copy of the GNU General Public License v2 (GPLv2) is available
in LICENSE.txt.  Users uncompressing this from an archive may not have
received this license file.  If not, see <http://www.gnu.org/licenses/>.
"""

@doc raw"""
	morris(EP::Model, path::AbstractString, setup::Dict, inputs::Dict, outpath::AbstractString, OPTIMIZER)

We apply the Method of Morris developed by [Morris, M., 1991](https://www.jstor.org/stable/1269043) in order to identify the input parameters that produce the largest change on total system cost. Method of Morris falls under the simplest class of one-factor-at-a-time (OAT) screening techniques. It assumes l levels per input factor and generates a set of trajectories through the input space. As such, the Method of Morris generates a grid of uncertain model input parameters, $x_i, i=1, ..., k$,, where the range $[x_i^{-}, x_i^{+}$ of each uncertain input parameter i is split into l intervals of equal length. Each trajectory starts at different realizations of input parameters chosen at random and are built by successively selecting one of the inputs randomly and moving it to an adjacent level. These trajectories are used to estimate the mean and the standard deviation of each input parameter on total system cost. A high estimated mean indicates that the input parameter is important; a high estimated standard deviation indicates important interactions between that input parameter and other inputs.
"""
struct MatSpread{T1,T2}
    mat::T1
    spread::T2
end

struct MorrisResult{T1,T2}
    means::T1
    means_star::T1
    variances::T1
    elementary_effects::T2
end
function generate_design_matrix(p_range, p_steps, rng;len_design_mat,groups)
    ps = [range(p_range[i][1], stop=p_range[i][2], length=p_steps[i]) for i in 1:length(p_range)]
    indices = [rand(rng, 1:i) for i in p_steps]
    all_idxs_original = Vector{typeof(indices)}(undef,len_design_mat)

    for i in 1:len_design_mat
        j = rand(rng, 1:length(p_range))
        indices[j] += (rand(rng) < 0.5 ? -1 : 1)
        if indices[j] > p_steps[j]
            indices[j] -= 2
        elseif indices[j] < 1.0
            indices[j] += 2
        end
        all_idxs_original[i] = copy(indices)
    end

    df_all_idx_original = DataFrame(all_idxs_original,:auto)
    println(df_all_idx_original)
    all_idxs = similar(df_all_idx_original)
    for g in unique(groups)
        temp = findall(x->x==g, groups)
        for k in temp
            all_idxs[k,:] = df_all_idx_original[temp[1],:]
        end
    end
    println(all_idxs)

    B = Array{Array{Float64}}(undef,len_design_mat)
    for j in 1:len_design_mat
        cur_p = [ps[u][(all_idxs[:,j][u])] for u in 1:length(p_range)]
        B[j] = cur_p
    end
    reduce(hcat, B)
end

function calculate_spread(matrix)
    spread = 0.0
    for i in 2:size(matrix,2)
        spread += sqrt(sum(abs2.(matrix[:,i] - matrix[:,i-1])))
    end
    spread
end

function sample_matrices(p_range,p_steps, rng;num_trajectory,total_num_trajectory,len_design_mat,groups)
    matrix_array = []
    println(num_trajectory)
    println(total_num_trajectory)
    if total_num_trajectory<num_trajectory
        error("total_num_trajectory should be greater than num_trajectory preferably atleast 3-4 times higher")
    end
    for i in 1:total_num_trajectory
        mat = generate_design_matrix(p_range, p_steps, rng;len_design_mat,groups)
        spread = calculate_spread(mat)
        push!(matrix_array,MatSpread(mat,spread))
    end
    sort!(matrix_array,by = x -> x.spread,rev=true)
    matrices = [i.mat for i in matrix_array[1:num_trajectory]]
    reduce(hcat,matrices)
end

function my_gsa(f, p_steps, num_trajectory, total_num_trajectory, p_range::AbstractVector,len_design_mat,groups)
    rng = Random.default_rng()
    design_matrices_original = sample_matrices(p_range, p_steps, rng;num_trajectory,
                                        total_num_trajectory,len_design_mat,groups)
    println(design_matrices_original)
    L = DataFrame(design_matrices_original,:auto)
    println(L)

    distinct_trajectories = Array{Int64}(undef,num_trajectory)
    design_matrices = Matrix(DataFrame(unique(last, pairs(eachcol(L[!,1:len_design_mat])))))
    distinct_trajectories[1] = length(design_matrices[1,:])
    if num_trajectory > 1
        for i in 2:num_trajectory
            design_matrices = hcat(design_matrices, Matrix(DataFrame(unique(last, pairs(eachcol(L[!,(i-1)*len_design_mat+1:i*len_design_mat]))))))
            distinct_trajectories[i] = length(Matrix(DataFrame(unique(last, pairs(eachcol(L[!,(i-1)*len_design_mat+1:i*len_design_mat])))))[1,:])
        end
    end
    println(distinct_trajectories)
    println(design_matrices)

    multioutput = false
    desol = false
    local y_size

    _y = [f(design_matrices[:,i]) for i in 1:size(design_matrices,2)]
    multioutput = !(eltype(_y) <: Number)
    if eltype(_y) <: RecursiveArrayTools.AbstractVectorOfArray
        y_size = size(_y[1])
        _y = vec.(_y)
        desol = true
    end
    all_y = multioutput ? reduce(hcat,_y) : _y
    println(all_y)
    effects = []
    while(length(effects) < length(groups))
        push!(effects,Vector{Float64}[])
    end
    for i in 1:num_trajectory
        len_design_mat = distinct_trajectories[i]
        y1 = multioutput ? all_y[:,(i-1)*len_design_mat+1] : all_y[(i-1)*len_design_mat+1]
        for j in (i-1)*len_design_mat+1:(i*len_design_mat)-1
            y2 = y1
            del = design_matrices[:,j+1] - design_matrices[:,j]
            change_index = 0
            for k in 1:length(del)
                if abs(del[k]) > 0
                    change_index = k
                    break
                end
            end
            del = sum(del)
            y1 = multioutput ? all_y[:,j+1] : all_y[j+1]
            effect = @. (y1-y2)/(del)
            elem_effect = typeof(y1) <: Number ? effect : mean(effect, dims = 2)
            temp_g_index = findall(x->x==groups[change_index], groups)
            for g in temp_g_index
                println(effects)
                println(elem_effect)
                push!(effects[g],elem_effect)
            end
        end
    end
    means = eltype(effects[1])[]
    means_star = eltype(effects[1])[]
    variances = eltype(effects[1])[]
    for k in effects
        if !isempty(k)
            push!(means, mean(k))
            push!(means_star, mean(x -> abs.(x), k))
            push!(variances, var(k))
        else
            push!(means, zero(effects[1][1]))
            push!(means_star, zero(effects[1][1]))
            push!(variances, zero(effects[1][1]))
        end
    end
    if desol
        f_shape = x -> [reshape(x[:,i],y_size) for i in 1:size(x,2)]
        means = map(f_shape,means)
        means_star = map(f_shape,means_star)
        variances = map(f_shape,variances)
    end
    MorrisResult(reduce(hcat, means),reduce(hcat, means_star),reduce(hcat, variances),effects)
end
function morris(EP::Model, path::AbstractString, setup::Dict, inputs::Dict, outpath::AbstractString, OPTIMIZER)

    # Reading the input parameters
    Morris_range = DataFrame(CSV.File(joinpath(path, "Method_of_morris_range.csv"), header=true), copycols=true)
    groups = Morris_range[!,:Group]
    p_steps = Morris_range[!,:p_steps]
    total_num_trajectory = Morris_range[!,:total_num_trajectory][1]
    num_trajectory = Morris_range[!,:num_trajectory][1]
    len_design_mat = Morris_range[!,:len_design_mat][1]
    uncertain_columns = unique(Morris_range[!,:Parameter])
    #save_parameters = zeros(length(Morris_range[!,:Parameter]))

    # Creating the range of uncertain parameters in terms of absolute values
    sigma = zeros((1, 2))
    for column in uncertain_columns
        sigma = [sigma; [inputs["dfGen"][!,Symbol(column)] .* (1 .+ Morris_range[Morris_range[!,:Parameter] .== column, :Lower_bound] ./100) inputs["dfGen"][!,Symbol(column)] .* (1 .+ Morris_range[Morris_range[!,:Parameter] .== column, :Upper_bound] ./100)]]
    end
    sigma = sigma[2:end,:]

    p_range = mapslices(x->[x], sigma, dims=2)[:]

    # Creating a function for iteratively solving the model with different sets of input parameters
    f1 = function(sigma)
        #print(sigma)
        print("\n")
        #save_parameters = hcat(save_parameters, sigma)

        for column in uncertain_columns
            index = findall(s -> s == column, Morris_range[!,:Parameter])
            inputs["dfGen"][!,Symbol(column)] = sigma[first(index):last(index)]
        end

        EP = generate_model(setup, inputs, OPTIMIZER)
        #EP, solve_time = solve_model(EP, setup)
        redirect_stdout((()->optimize!(EP)),open("/dev/null", "w"))
        [objective_value(EP)]
    end

    # Perform the method of morris analysis
    m = my_gsa(f1,p_steps,num_trajectory,total_num_trajectory,p_range,len_design_mat,groups)
    println(m.means)
    println(DataFrame(m.means', :auto))
    #save the mean effect of each uncertain variable on the objective fucntion
    Morris_range[!,:mean] = DataFrame(m.means', :auto)[!,:x1]
    println(DataFrame(m.variances', :auto))
    #save the variance of effect of each uncertain variable on the objective function
    Morris_range[!,:variance] = DataFrame(m.variances', :auto)[!,:x1]

    CSV.write(joinpath(outpath, "morris.csv"), Morris_range)

end

max_stars_repo_name: Betristor/LTESOM

max_issues_repo_name: Betristor/LTESOM

max_forks_repo_name: Betristor/LTESOM


text:
export create_scs_matrix, create_scs_settings, create_scs_data, create_scs_cone

# Takes a vector or matrix or sparse matrix A and creates an SCSMatrix
function create_scs_matrix(m::Int, n::Int, A::SCSVecOrMatOrSparse)
    A_sparse = sparse(A)

    values = A_sparse.nzval * 1.0
    rowval = convert(Array{Int, 1}, A_sparse.rowval .- 1)
    colptr = convert(Array{Int, 1}, A_sparse.colptr .- 1)

    return SCSMatrix(pointer(values), pointer(rowval), pointer(colptr), m, n)
end

function create_scs_settings(normalize=1::Int, scale=converter(Cdouble, 5.0)::Cdouble, rho_x=convert(Cdouble,1e-3)::Cdouble,
                        max_iters=2500::Int, eps=converter(Cdouble, 1e-3)::Cdouble, alpha=convert(Cdouble, 1.8)::Cdouble,
                        cg_rate=convert(Cdouble,2)::Cdouble, verbose=1::Int, warm_start=0::Int)
    return SCSSettings(normalize, scale, rho_x, max_iters, eps, alpha, cg_rate, verbose, warm_start)
end

# Create an SCSData type
# We assume we are solving a problem of the form
# minimize        c' * x
# subject to      A * x + s = b
#                 s in K
# A is the matrix with m rows and n cols
# b is of length m x 1
# c is of length n x 1
# refer to create_scs_cone for K
function create_scs_data(;m::Int=nothing, n::Int=nothing, A::Ptr{SCSMatrix}=nothing,
        b::Ptr{Cdouble}=nothing,  c::Ptr{Cdouble}=nothing, max_iters=20000::Int,
        eps=convert(Cdouble, 1e-4)::Cdouble, alpha=convert(Cdouble, 1.8)::Cdouble,
        rho_x=convert(Cdouble, 1e-3)::Cdouble, scale=convert(Cdouble, 5.0)::Cdouble,
        cg_rate=convert(Cdouble, 2)::Cdouble, verbose=1::Int,
        normalize=1::Int, warm_start=0::Int, options...)

    for (k, v) in options
        @eval(($k) = ($v))
    end
    stgs = create_scs_settings(normalize, scale, rho_x, max_iters, eps, alpha, cg_rate, verbose, warm_start)
    return SCSData(m, n, A, b, c, pointer([stgs]))
end

# Refer to comment above
function create_scs_data(m::Int, n::Int, A::Ptr{SCSMatrix}, b::Ptr{Cdouble}, c::Ptr{Cdouble}; options...)
    return create_scs_data(m=m, n=n, A=A, b=b, c=c; options...)
end

# Refer to comment above
function create_scs_data(m::Int, n::Int, A::SCSVecOrMatOrSparse, b::Array{Float64,},
        c::Array{Float64,}; options...)
    if size(b, 1) != m || size(b, 2) != 1 || size(c, 1) != n || size(c, 2) != 1
        error("Size of b must be m x 1 and size of c must be n x 1")
    end
    A = [create_scs_matrix(m, n, A)]
    return create_scs_data(m=m, n=n, A=pointer(A), b=pointer(b), c=pointer(c); options...)
end

# Create an SCSCone type
# We assume we are solving a problem of the form
# minimize        c' * x
# subject to      A * x + s = b
#                 s in K
# where K is a product cone of
# zero cones,
# linear cones { x | x >= 0 },
# second-order cones { (t,x) | ||x||_2 <= t },
# semi-definite cones { X | X psd }, and
# exponential cones {(x,y,z) | y e^(x/y) <= z, y>0 }.
#
# f (num primal zero / dual free cones, i.e. primal equality constraints)
# l (num linear cones)
# q (array of SOCs sizes)
# s (array of SDCs sizes)
# ep (num primal exponential cones)
# ed (num dual exponential cones).
function create_scs_cone(f::Int, l::Int, q::Ptr{Int}, qsize::Int, s::Ptr{Int},
        ssize::Int, ep::Int, ed::Int, p::Ptr{Float64}, psize::Int)
    return SCSCone(f, l, q, qsize, s, ssize, ep, ed, p, psize)
end

# Refer to comment above
function create_scs_cone(f::Int, l::Int, q::Array{Int,}, qsize::Int, s::Array{Int,},
        ssize::Int, ep::Int, ed::Int, p::Array{Float64,}, psize::Int)
    return SCSCone(f, l, pointer(q), qsize, pointer(s), ssize, ep, ed, pointer(p), psize)
end

# We assume we are solving a problem of the form
# minimize        c' * x
# subject to      A * x + s = b
#                 s in K
# where K is a product cone of
# zero cones,
# linear cones { x | x >= 0 },
# second-order cones { (t,x) | ||x||_2 <= t },
# semi-definite cones { X | X psd }, and
# exponential cones {(x,y,z) | y e^(x/y) <= z, y>0 }.
#
#
# A is the matrix with m rows and n cols
# b is of length m x 1
# c is of length n x 1
#
# f (num primal zero / dual free cones, i.e. primal equality constraints)
# l (num linear cones)
# q (array of SOCs sizes)
# s (array of SDCs sizes)
# ep (num primal exponential cones)
# ed (num dual exponential cones).
#
# Returns object of type Solution
# type Solution with
# x, y, s, status (ASCII string), ret_val (numerical status)
#
function SCS_solve(m::Int, n::Int, A::SCSVecOrMatOrSparse, b::Array{Float64,},
        c::Array{Float64,}, f::Int, l::Int, q::Array{Int,}, qsize::Int, s::Array{Int,},
        ssize::Int, ep::Int, ed::Int, p::Array{Float64,}, psize::Int,
        primal_sol::Vector{Float64}=Float64[],
        dual_sol::Vector{Float64}=Float64[],
        slack::Vector{Float64}=Float64[];
        options...)

    options_dict = Dict(options)
    if :linearsolver in keys(options_dict)
        T = options_dict[:linearsolver]
        options = [(a,b) for (a,b) in options if a≠:linearsolver]
    else
        T = SCS.Direct
    end

    data = create_scs_data(m, n, A, b, c; options...)
    cone = create_scs_cone(f, l, q, qsize, s, ssize, ep, ed, p, psize)

    if (:warm_start, true) in options && length(primal_sol) == n && length(dual_sol) == m && length(slack) == m
        x = primal_sol
        y = dual_sol
        s = slack
    else
        x = zeros(n)
        y = zeros(m)
        s = zeros(m)
    end
    solution = SCSSolution(pointer(x), pointer(y), pointer(s))
    status, solution, info, p_work = SCS_solve(T, data, cone, solution)
    SCS_finish(T, p_work)
    return Solution(x, y, s, status)

end

# for legacy API maintenance, sets power cones to zero
function SCS_solve(m::Int, n::Int, A::SCSVecOrMatOrSparse, b::Array{Float64,},
        c::Array{Float64,}, f::Int, l::Int, q::Array{Int,}, qsize::Int, s::Array{Int,},
        ssize::Int, ep::Int, ed::Int,
        primal_sol::Vector{Float64}=Float64[],
        dual_sol::Vector{Float64}=Float64[],
        slack::Vector{Float64}=Float64[];
        options...)
        return SCS_solve(m, n, A, b, c, f, l, q, qsize, s, ssize,
        ep, ed, Float64[], 0, primal_sol, dual_sol, slack; options...)
end


text:
file = joinpath(dir,"testfile.txt")
file2 = joinpath(dir,"sal.txt")

@testset "Line Count Testing" begin
    # Malformed line
    ml = "aaaa\nbbb\ncccc\n"
    #No Ending NL
    nonl = "aaaa\nbbbb\ncccc"
    #Extra junk at end
    extra = "aaaa\nbbbb\ncccc\ndd"
    #Carrige Retruns
    cr = "aaaa\r\nbbbb\r\ncccc\r\n"
    #Normal
    b = "aaaa\nbbbb\ncccc\n"
    # no data
    nodata = ""
    #test base.countlines for fix to counting so we can adjust code
    @test countlines(IOBuffer(nonl)) == 2
    @test countlines(IOBuffer(b)) == 3
    @test FWF.mod_countlines(IOBuffer(nonl)) == 3
    @test FWF.mod_countlines(IOBuffer(b)) == 3
    @test FWF.mod_countlines(IOBuffer(" ")) == 1
    @test FWF.mod_countlines(IOBuffer("")) == 0
    # Not error until malformed is back
    #@test_throws FWF.ParsingException FWF.row_countlines(IOBuffer(ml))
    @test FWF.row_countlines(IOBuffer(ml)) == (3, 0)
    @test FWF.row_countlines(IOBuffer(nonl)) == (3, 0)
    # Used to be an error condition not now.
    #@test_throws FWF.ParsingException FWF.row_countlines(IOBuffer(extra))
    @test FWF.row_countlines(IOBuffer(extra)) == (4, 0)
    @test FWF.row_countlines(IOBuffer(b)) == (3, 0)
    @test FWF.row_countlines(IOBuffer(cr)) == (3, 1)
    @test FWF.row_countlines(IOBuffer(nodata)) == (0, 0)
end

@testset "readsplitline! Testing" begin
    s = Vector{Union{Missing,String}}()
    tmp = FWF.Source(file, [4,4,8])
    FWF.readsplitline!(s, tmp)
    @test s[1] == "abcd"
    @test s[2] == "1234"
    @test s[3] == "10102017"
    @test_throws ArgumentError FWF.readsplitline!(s, tmp.io, Vector{UnitRange{Int}}())
    FWF.readsplitline!(s, tmp)
    @test s[1] == "efgh"
    @test s[2] == "5678"
    @test s[3] == "10112017"
    @test_throws ArgumentError FWF.readsplitline!(s, tmp)

    tmp = FWF.Source(file, [4,4,4])
    FWF.readsplitline!(s, tmp)
    @test s[1] == "abcd"
    @test s[2] == "1234"
    @test s[3] == "1010"
    @test_throws ArgumentError FWF.readsplitline!(s, tmp.io, Vector{UnitRange{Int}}())
    FWF.readsplitline!(s, tmp)
    @test s[1] == "efgh"
    @test s[2] == "5678"
    @test s[3] == "1011"
    @test_throws ArgumentError FWF.readsplitline!(s, tmp)

    tmp = FWF.Source(file, [1:4,5:8,9:16])
    FWF.readsplitline!(s, tmp)
    @test s[1] == "abcd"
    @test s[2] == "1234"
    @test s[3] == "10102017"
    @test_throws ArgumentError FWF.readsplitline!(s, tmp.io, Vector{UnitRange{Int}}())
    FWF.readsplitline!(s, tmp)
    @test s[1] == "efgh"
    @test s[2] == "5678"
    @test s[3] == "10112017"
    @test_throws ArgumentError FWF.readsplitline!(s, tmp)

    tmp = FWF.Source(file, [1:2,5:6,9:10])
    FWF.readsplitline!(s, tmp)
    @test s[1] == "ab"
    @test s[2] == "12"
    @test s[3] == "10"
    @test_throws ArgumentError FWF.readsplitline!(s, tmp.io, Vector{UnitRange{Int}}())
    FWF.readsplitline!(s, tmp)
    @test s[1] == "ef"
    @test s[2] == "56"
    @test s[3] == "10"
    @test_throws ArgumentError FWF.readsplitline!(s, tmp)
    Data.reset!(tmp)
    @test Data.accesspattern(tmp) == Data.Sequential() # not sure why this isn't a type
    FWF.readsplitline!(s, tmp)
    @test s[1] == "ab"
    @test s[2] == "12"
    @test s[3] == "10"
    @test_throws ArgumentError FWF.readsplitline!(s, tmp.io, Vector{UnitRange{Int}}())

    b="""
    aaaa
    bbb
    cccc"""
    tmp = FWF.Source(IOBuffer(b),[4], errorlevel=:skip)
    @test FWF.readsplitline!(s, tmp) == 4
    @test s[1] == "aaaa"
    @test FWF.readsplitline!(s, tmp) == 4
    @test s[1] == "cccc"
    # Not error until malformed is back
    #@test_throws FWF.ParsingException FWF.Source(IOBuffer(b),[4],skiponerror=false)
    tmp = FWF.Source(IOBuffer(b),[4], errorlevel=:skip)
    @test FWF.readsplitline!(s, tmp) == 4
    @test s[1] == "aaaa"
    # @test tmp.schema.rows == 3
    #@test_throws FWF.ParsingException FWF.readsplitline!(s, tmp)
    @test FWF.readsplitline!(s, tmp) == 4
    @test s[1] == "cccc"

    # test overlong line
    b="""
    aaaa
    bbb
    ccccc"""
    tmp = FWF.Source(IOBuffer(b),[4], errorlevel=:skip)
    @test FWF.readsplitline!(s, tmp) == 4
    @test s[1] == "aaaa"
    @test FWF.readsplitline!(s, tmp) == 5
    @test s[1] == "cccc"

    # test overlong line
    b="""
    aaaaa
    bbb
    cccc"""
    tmp = FWF.Source(IOBuffer(b),[4], errorlevel=:skip)
    @test FWF.readsplitline!(s, tmp) == 5
    @test s[1] == "aaaa"
    @test FWF.readsplitline!(s, tmp) == 4
    @test s[1] == "cccc"

    #ensure utf 8 doesn't mess us up
    b = """
    abcd
    \u263ae
    fghi
    """
    tmp = FWF.Source(IOBuffer(b),[4])
    @test FWF.readsplitline!(s, tmp) == 4
    @test s[1] == "abcd"
    @test FWF.readsplitline!(s, tmp) == 4
    @test s[1] == "\u263ae"
    @test FWF.readsplitline!(s, tmp) == 4
    @test s[1] == "fghi"
    @test_throws ArgumentError FWF.readsplitline!(s, tmp)
end

@testset "FWF.read Testing" begin
    tmp = FWF.read(IOBuffer("abc12310102017\ndef45610112017\n"), [3,3,8], types=[String,Int,DateFormat("mmddyyyy")])
    @test tmp[1,1] == "abc"
    @test tmp[1,2] == 123
    @test tmp[1,3] == Date(2017,10,10)
    @test tmp[2,1] == "def"
    @test tmp[2,2] == 456
    @test tmp[2,3] == Date(2017,10,11)
    strrep(s::String, r::UnitRange) = [repeat(s, n) for n in r]
    naValues = vcat(strrep("*", 1:23), strrep("#", 1:23), "NAME WITHHELD BY AGENCY", "NAME WITHHELD BY OPM", "NAME UNKNOWN", "UNSP", "<NA>", "000000", "999999", "")

    #    ~FieldName,~Length,~Type,
    format = DataFrame(["PSEUDO_ID" 9 Int;
    "EMPLOYEE_NAME" 23 String;
    "FILE_DATE" 8  DateFormat("yyyymmdd");
    "AGENCY" 2 String;
    "SUB_AGENCY" 2 String;
    "DUTY_STATION" 9 String;
    "AGE" 6 String;
    "EDUCATION_LEVEL" 2 String;
    "PAY_PLAN" 2 String;
    "GRADE" 2 String;
    "LOS_LEVEL" 6 String;
    "OCCUPATION" 4 String;
    "OCCUPATIONAL_CATEGORY" 1 String;
    "ADJUSTED_BASIC_PAY" 6 Int;
    "SUPERVISORY_STATUS" 1 String;
    "TYPE_OF_APPOINTMENT" 2 String;
    "WORK_SCHEDULE" 1 String;
    "NSFTP_IND"  1 String])

    tmp = FWF.read(file2, convert(Array{Int},format[:x2]),
        header=convert(Array{String}, format[:x1]), types=convert(Array{Union{Type, DateFormat}},format[:x3]),
        missings=naValues)

    @test ismissing(tmp[2,8])
    @test tmp[150651,4] == "CM"
    @test ismissing(tmp[15321,13])
    @test tmp[15333,13] == "P"
    @test tmp[153242,16] == "15"
    @test tmp[9999,4] == "AG"
    @test tmp[36, 2] == "NAME WITHHELD BY ÄĜË"

    # Simple UTF-8 test
    tmp = FWF.read(IOBuffer("α1x\na2y\n∀∅z"), [1,1,1], unitbytes=false)
    @test tmp[1,1] == "α"
    @test tmp[1,2] == "1"
    @test tmp[2,2] == "2"
    @test tmp[3,2] == "∅"
    @test tmp[3,3] == "z"
    tmp = FWF.read(IOBuffer("α1x\na2y\n∀∅z"), [2,1], unitbytes=false)
    @test tmp[1,1] == "α1"
    @test tmp[2,1] == "a2"
    @test tmp[3,1] == "∀∅"
    @test tmp[3,2] == "z"
    tmp = FWF.read(IOBuffer("α1x\na2y\n∀∅z"), [3], unitbytes=false)
    @test tmp[1,1] == "α1x"
    @test tmp[2,1] == "a2y"
    @test tmp[3,1] == "∀∅z"

    tmp = FWF.read(IOBuffer("α1\na\n∀∅z"), [1,1,1], unitbytes=false)
    @test tmp[1,1] == "α"
    @test tmp[1,2] == "1"
    @test ismissing(tmp[2,2])
    @test tmp[3,2] == "∅"
    @test tmp[3,3] == "z"
    @test ismissing(tmp[1,3])
    @test ismissing(tmp[2,3])

    @test_throws FWF.ParsingException FWF.read(IOBuffer("α1\na\n∀∅z"), [1,1,1], unitbytes=false, errorlevel=:error)

    tmp = FWF.read(IOBuffer("α1\na\n∀∅z"), [1,1,1], unitbytes=false, errorlevel=:skip)
    @test nrow(tmp) == 1
    @test tmp[1,1] == "∀"
    @test tmp[1,2] == "∅"
    @test tmp[1,3] == "z"

    tmp = FWF.read(IOBuffer("αąx\na2y\n∀∅z"), [1,1,1], unitbytes=false, header=true)
    @test names(tmp) == [:α, :ą, :x]
    @test nrow(tmp) == 2

    tmp = FWF.read(IOBuffer("αąx\naby\n∀∅z"), [1,1,1], unitbytes=false, header=true, skip = 1)
    @test names(tmp) == [:a, :b, :y]
    @test nrow(tmp) == 1

    @test_throws ArgumentError FWF.read(IOBuffer("αąx\na2y\n∀∅z"), [1,1,1], unitbytes=false, header=true, skip = 2)

    @test_throws ArgumentError FWF.read(IOBuffer("αąx\na2y\n∀∅z"), [1,1,1], unitbytes=false, header=false, skip = 3)
    @test nrow(FWF.read(IOBuffer("αąx\na2y\n∀∅z"), [1,1,1], unitbytes=false, header=false, skip = 2)) == 1
    @test nrow(FWF.read(IOBuffer("αąx\na2y\n∀∅z"), [1,1,1], unitbytes=false, header=false, skip = 1)) == 2
end


text:
# AUTO GENERATED FILE - DO NOT EDIT

export image

"""
    image(;kwargs...)

An Image component.
Image with optional placeholder for loading and error state. For more information, see: https://mantine.dev/core/image/
Keyword arguments:
- `id` (String; optional): The ID of this component, used to identify dash components in callbacks
- `alt` (String; optional): Image alt text, used as title for placeholder if image was not loaded
- `caption` (String; optional): Image figcaption, displayed bellow image
- `className` (String; optional): Often used with CSS to style elements with common properties
- `fit` (a value equal to: "cover", "contain"; optional): Image object-fit property
- `height` (String | Real; optional): Image height, defaults to original image height adjusted to given width
- `radius` (a value equal to: "xs", "sm", "md", "lg", "xl" | Real; optional): Predefined border-radius value from theme.radius or number for border-radius in px
- `src` (String; optional): Image src
- `style` (Dict; optional): Inline style override
- `width` (String | Real; optional): Image width, defaults to 100%, cannot exceed 100%
- `withPlaceholder` (Bool; optional): Enable placeholder when image is loading and when image fails to load
"""
function image(; kwargs...)
        available_props = Symbol[:id, :alt, :caption, :className, :fit, :height, :radius, :src, :style, :width, :withPlaceholder]
        wild_props = Symbol[]
        return Component("image", "Image", "dash_mantine_components", available_props, wild_props; kwargs...)
end


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

function assure_ordering!(ring, exps, coeffs, metainfo)
    if ring.ord != metainfo.computeord
        sort_input_to_change_ordering(exps, coeffs, metainfo.computeord)
    end
    ring.ord = metainfo.computeord
end

# TODO: change project structure

function input_statistics(exps)
    sz = length(exps)
    deg = maximum(sum(e[1]) for e in exps)
    sz, deg
end

function select_tablesize(ring, exps)
    nvars = ring.nvars
    sz = length(exps)

    tablesize = 2^10
    if nvars > 4
        tablesize = 2^14
    end
    if nvars > 7
        tablesize = 2^16
    end

    if sz < 3
        tablesize = div(tablesize, 2)
    end
    if sz < 2
        tablesize = div(tablesize, 2)
    end

    tablesize
end

function cleanup_gens!(ring, gens_ff, prime)
    gens_ff.ch = prime
    ring.ch = prime
    normalize_basis!(gens_ff)
end

# Mutate everything!

#------------------------------------------------------------------------------
# Finite field groebner

function groebner_ff(
            ring::PolyRing,
            exps::Vector{Vector{Vector{UInt16}}},
            coeffs::Vector{Vector{UInt64}},
            reduced::Bool,
            meta::GroebnerMetainfo) where {Rng<:Random.AbstractRNG}

    # select hashtable size
    tablesize = select_tablesize(ring, exps)

    basis, ht = initialize_structures(
                        ring, exps, coeffs, meta.rng, tablesize)

    f4!(ring, basis, ht, reduced, meta.linalg, meta.rng)

    gbexps = hash_to_exponents(basis, ht)
    gbexps, basis.coeffs
end

#------------------------------------------------------------------------------
# Rational field groebner

function groebner_qq(
            ring::PolyRing,
            exps::Vector{Vector{ExponentVector}},
            coeffs::Vector{Vector{CoeffQQ}},
            reduced::Bool,
            meta::GroebnerMetainfo) where {Rng<:Random.AbstractRNG}

    # we can mutate coeffs and exps here

    # select hashtable size
    tablesize = select_tablesize(ring, exps)
    @info "Selected tablesize $tablesize"

    # initialize hashtable and finite field generators structs
    gens_temp_ff, ht = initialize_structures_ff(ring, exps,
                                        coeffs, meta.rng, tablesize)
    gens_ff = copy_basis_thorough(gens_temp_ff)

    # now hashtable is filled correctly,
    # and gens_temp_ff exponents are correct and in correct order.
    # gens_temp_ff coefficients are filled with random stuff and
    # gens_temp_ff.ch is 0

    # to store integer and rational coefficients of groebner basis
    coeffaccum  = CoeffAccum()
    # to store BigInt buffers and reduce overall memory usage
    coeffbuffer = CoeffBuffer()

    # scale coefficients of input to integers
    coeffs_zz = scale_denominators(coeffbuffer, coeffs)

    # keeps track of used prime numbers
    primetracker = PrimeTracker(coeffs_zz)

    #=
        exps, coeffs and coeffs_zz must be *not changed* during whole computation
    =#

    i = 1

    # copy basis so that we initial exponents dont get lost
    gens_ff = copy_basis_thorough(gens_temp_ff)

    prime = nextluckyprime!(primetracker)
    @info "$i: selected lucky prime $prime"

    # perform reduction and store result in gens_ff
    reduce_modulo!(coeffbuffer, coeffs_zz, gens_ff.coeffs, prime)

    # do some things to ensure generators are correct
    cleanup_gens!(ring, gens_ff, prime)

    # compute groebner basis in finite field
    #=
    Need to make sure input invariants in f4! are satisfied, f4.jl for details
    =#

    tracer = Tracer()

    pairset = initialize_pairset()

    # ADDED
    global F4TIME
    # ADDED
    F4TIME += @elapsed f4trace!(ring, gens_ff, tracer, pairset, ht, reduced, meta.linalg, meta.rng)
    # f4!(ring, gens_ff, ht, reduced, meta.linalg)

    # reconstruct into integers
    @info "CRT modulo ($(primetracker.modulo), $(prime))"

    # ADDED
    global RECTIME
    # ADDED
    RECTIME += @elapsed reconstruct_crt!(coeffbuffer, coeffaccum, primetracker, gens_ff.coeffs, prime)
    # reconstruct_crt!(coeffbuffer, coeffaccum, primetracker, gens_ff.coeffs, prime)

    # reconstruct into rationals
    @info "Reconstructing modulo $(primetracker.modulo)"
    # ADDED
    RECTIME += @elapsed reconstruct_modulo!(coeffbuffer, coeffaccum, primetracker)
    # reconstruct_modulo!(coeffbuffer, coeffaccum, primetracker)

    correct = false
    # ADDED
    global CORRTIME
    t = time()
    if correctness_check!(coeffaccum, coeffbuffer, primetracker, meta,
                            ring, exps, coeffs, coeffs_zz, gens_temp_ff, gens_ff, ht)
        @info "Success!"
        CORRTIME += time() - t
        correct = true
    end
    CORRTIME += time() - t

    gap = 1
    primegaps = (1, 1, 1, 1, 1)
    primemult = 2
    # if first prime was not successfull
    while !correct

        if i <= length(primegaps)
            gap = primegaps[i]
        else
            gap = gap*primemult
        end

        for j in 1:gap
            i += 1

            # copy basis so that we initial exponents dont get lost
            gens_ff = copy_basis_thorough(gens_temp_ff)
            prime = nextluckyprime!(primetracker)
            @info "$i: selected lucky prime $prime"
            # perform reduction and store result in gens_ff
            reduce_modulo!(coeffbuffer, coeffs_zz, gens_ff.coeffs, prime)
            # do some things to ensure generators are correct
            cleanup_gens!(ring, gens_ff, prime)
            # compute groebner basis in finite field
            #=
            Need to make sure input invariants in f4! are satisfied, f4.jl for details
            =#
            # ADDED
            global F4TIME
            # ADDED

            # @error "tracer" tracer

            F4TIME += @elapsed f4trace!(ring, gens_ff, tracer, pairset, ht, reduced, meta.linalg, meta.rng)
            # f4!(ring, gens_ff, ht, reduced, meta.linalg)
            # reconstruct into integers
            @info "CRT modulo ($(primetracker.modulo), $(prime))"
            # ADDED
            global RECTIME
            # ADDED

            if meta.linalg === :prob
                # check if probabilistic linear algebra failed
                shape_control(coeffaccum.gb_coeffs_qq, gens_ff.coeffs)
            end

            # ADDED
            # prevcoeffs = deepcopy(coeffaccum.gb_coeffs_zz)

            RECTIME += @elapsed reconstruct_crt!(coeffbuffer, coeffaccum,
                                        primetracker, gens_ff.coeffs, prime)
            # reconstruct_crt!(coeffbuffer, coeffaccum, primetracker, gens_ff.coeffs, prime)

            #=
            -----------------------------------
            Given r1, r2, and m1, m2, so that
                r1 ≡ 1/b1 (mod m1)
                r2 ≡ 1/b2 (mod m2)
            Determine if b1 = b2
            -----------------------------------

            Easy:
            b1 is inverse of r1 mod m1;
            b2 is inverse of r2 mod m2;

            then check equality
            =#

            #=
            -----------------------------------
            We go further.
            Given r1, r2, and m1, m2, so that
                r1 ≡ a1/b1 (mod m1)
                r2 ≡ a2/b2 (mod m2)
            Determine if a1 = a2 and b1 = b2
            -----------------------------------
            r1*b1 ≡ a1 (mod m1),
            r2*b2 ≡ a2 (mod m2),

            r1*b1 * r2*b2 = a1*a2 (mod m1*m2)
            =#

            #=
            if j != gap
                if crt_issame_check!(coeffbuffer, coeffaccum,
                                gens_ff.coeffs, primetracker.modulo, prime)
                    @info "CRT same check passed!"
                    break
                end
            end
            =#
        end

        # reconstruct into rationals
        @info "Reconstructing modulo $(primetracker.modulo)"
        # ADDED
        stats = @timed reconstruct_modulo!(coeffbuffer, coeffaccum, primetracker)
        RECTIME += stats.time

        # ADDED TODO
        if !stats.value
            continue
        end

        # reconstruct_modulo!(coeffbuffer, coeffaccum, primetracker)
        # ADDED
        global CORRTIME
        t = time()
        if correctness_check!(coeffaccum, coeffbuffer, primetracker, meta,
                                ring, exps, coeffs, coeffs_zz, gens_temp_ff, gens_ff, ht)
            @info "Success!"
            CORRTIME += time() - t
            correct = true
        end
        CORRTIME += time() - t
        #=
        if correctness_check!(coeffaccum, coeffbuffer, primetracker, meta,
                                ring, exps, coeffs, coeffs_zz, gens_temp_ff, gens_ff, ht)
            @info "Success!"
            break
        end
        =#

        if meta.linalg === :prob
            if i > 1024
                coeff_control(primetracker, coeffaccum)
            end
        end
    end

    #=
    if meta.usefglm
        fglm_f4!(ring, gens_ff, ht)
    end
    =#

    # ADDED
    global PRIMES
    PRIMES = i

    # normalize_coeffs!(gbcoeffs_qq)
    gb_exps = hash_to_exponents(gens_ff, ht)
    gb_exps, coeffaccum.gb_coeffs_qq
end

#######################################
# Finite field isgroebner

# UWU!
function isgroebner_ff(
            ring::PolyRing,
            exps::Vector{Vector{ExponentVector}},
            coeffs::Vector{Vector{CoeffFF}},
            meta)

    isgroebner_f4(ring, exps, coeffs, meta.rng)
end

#######################################
# Rational field groebner

# TODO
function isgroebner_qq(
            ring::PolyRing,
            exps::Vector{Vector{ExponentVector}},
            coeffs::Vector{Vector{CoeffQQ}},
            meta)

    # if randomized result is ok
    if !meta.guaranteedcheck
        coeffs_zz = scale_denominators(coeffs)
        goodprime = nextgoodprime(coeffs_zz, Int[], 2^30 + 3)
        ring_ff, coeffs_ff = reduce_modulo(coeffs_zz, ring, goodprime)
        isgroebner_f4(ring_ff, exps, coeffs_ff, meta.rng)
    else # if proved result is needed, compute in rationals
        isgroebner_f4(ring, exps, coeffs, meta.rng)
    end
end