text:
"""
Returns a handle for the SPI device on channel.  Data will be
transferred at baud bits per second.  The flags may be used to
modify the default behaviour of 4-wire operation, mode 0,
active low chip select.

An auxiliary SPI device is available on all models but the
A and B and may be selected by setting the A bit in the
flags. The auxiliary device has 3 chip selects and a
selectable word size in bits.

spi_channel:= 0-1 (0-2 for the auxiliary SPI device).
 baud:= 32K-125M (values above 30M are unlikely to work).
spi_flags:= see below.

Normally you would only use the [*spi_**] functions if
you are or will be connecting to the Pi over a network.  If
you will always run on the local Pi use the standard SPI
module instead.

spi_flags consists of the least significant 22 bits.

. .
21 20 19 18 17 16 15 14 13 12 11 10  9  8  7  6  5  4  3  2  1  0
b  b  b  b  b  b  R  T  n  n  n  n  W  A u2 u1 u0 p2 p1 p0  m  m
. .

mm defines the SPI mode.

WARNING: modes 1 and 3 do not appear to work on
the auxiliary device.

. .
Mode POL PHA
0    0   0
1    0   1
2    1   0
3    1   1
. .

px is 0 if CEx is active low (default) and 1 for active high.

ux is 0 if the CEx GPIO is reserved for SPI (default)
and 1 otherwise.

A is 0 for the standard SPI device, 1 for the auxiliary SPI.

W is 0 if the device is not 3-wire, 1 if the device is 3-wire.
Standard SPI device only.

nnnn defines the number of bytes (0-15) to write before
switching the MOSI line to MISO to read data.  This field
is ignored if W is not set.  Standard SPI device only.

T is 1 if the least significant bit is transmitted on MOSI
first, the default (0) shifts the most significant bit out
first.  Auxiliary SPI device only.

R is 1 if the least significant bit is received on MISO
first, the default (0) receives the most significant bit
first.  Auxiliary SPI device only.

bbbbbb defines the word size in bits (0-32).  The default (0)
sets 8 bits per word.  Auxiliary SPI device only.

The [*spi_read*], [*spi_write*], and [*spi_xfer*] functions
transfer data packed into 1, 2, or 4 bytes according to
the word size in bits.

For bits 1-8 there will be one byte per character.
For bits 9-16 there will be two bytes per character.
For bits 17-32 there will be four bytes per character.

E.g. 32 12-bit words will be transferred in 64 bytes.

The other bits in flags should be set to zero.

...
# open SPI device on channel 1 in mode 3 at 50000 bits per second

h = spi_open(pi, 1, 50000, 3)
...
"""
function spi_open(self::Pi, spi_channel, baud, spi_flags=0)
    # I p1 spi_channel
    # I p2 baud
    # I p3 4
    ## extension ##
    # I spi_flags
    extents=IOBuffer()
    write(extents, spi_flags::Cint)
    return _u2i(_pigpio_command_ext(
        self.sl, _PI_CMD_SPIO, spi_channel, baud, 4, extents))
end

"""
Closes the SPI device associated with handle.

handle:= >=0 (as returned by a prior call to [*spi_open*]).

...
spi_close(pi, h)
...
"""
function spi_close(self, handle)
    return _u2i(_pigpio_command(self.sl, _PI_CMD_SPIC, handle, 0))
end

"""
Reads count bytes from the SPI device associated with handle.

handle:= >=0 (as returned by a prior call to [*spi_open*]).
count:= >0, the number of bytes to read.

The returned value is a tuple of the number of bytes read and a
bytearray containing the bytes.  If there was an error the
number of bytes read will be less than zero (and will contain
the error code).

...
(b, d) = spi_read(pi, h, 60) # read 60 bytes from device h
if b == 60
# process read data
else
# error path
...
"""
function spi_read(self::Pi, handle, count)
    # Don't raise exception.  Must release lock.
    bytes = u2i(_pigpio_command(
    self.sl, _PI_CMD_SPIR, handle, count, false))
    if bytes > 0
        data = rxbuf(bytes)
    else
        data = ""
    end
    unlock(self.sl.l)
    return bytes, data
end

"""
Writes the data bytes to the SPI device associated with handle.

handle:= >=0 (as returned by a prior call to [*spi_open*]).
data:= the bytes to write.

...
spi_write(pi, 0, b'\\x02\\xc0\\x80') # write 3 bytes to device 0

spi_write(pi, 0, b'defgh')        # write 5 bytes to device 0

spi_write(pi, 0, "def")           # write 3 bytes to device 0

spi_write(pi, 1, [2, 192, 128])   # write 3 bytes to device 1
...
"""
function spi_write(self::Pi, handle, data)
    # I p1 handle
    # I p2 0
    # I p3 len
    ## extension ##
    # s len data bytes
    return _u2i(_pigpio_command_ext(
        self.sl, _PI_CMD_SPIW, handle, 0, length(data), data))
end

"""
Writes the data bytes to the SPI device associated with handle,
returning the data bytes read from the device.

handle:= >=0 (as returned by a prior call to [*spi_open*]).
data:= the bytes to write.

The returned value is a tuple of the number of bytes read and a
bytearray containing the bytes.  If there was an error the
number of bytes read will be less than zero (and will contain
the error code).

...
(count, rx_data) = spi_xfer(pi, h, b'\\x01\\x80\\x00')

(count, rx_data) = spi_xfer(pi, h, [1, 128, 0])

(count, rx_data) = spi_xfer(pi, h, b"hello")

(count, rx_data) = spi_xfer(pi, h, "hello")
...
"""
function spi_xfer(self::Pi, handle, data)
    # I p1 handle
    # I p2 0
    # I p3 len
    ## extension ##
    # s len data bytes

    # Don't raise exception.  Must release lock.
    bytes = u2i(_pigpio_command_ext(
    self.sl, _PI_CMD_SPIX, handle, 0, length(data), data, false))
    if bytes > 0
        data = rxbuf(bytes)
    else
        data = ""
    end
    unlock(self.sl.l)
    return bytes, data
end

"""
Returns a handle for the serial tty device opened
at baud bits per second.

tty:= the serial device to open.
baud:= baud rate in bits per second, see below.
ser_flags:= 0, no flags are currently defined.

Normally you would only use the [*serial_**] functions if
you are or will be connecting to the Pi over a network.  If
you will always run on the local Pi use the standard serial
module instead.

The baud rate must be one of 50, 75, 110, 134, 150,
200, 300, 600, 1200, 1800, 2400, 4800, 9600, 19200,
38400, 57600, 115200, or 230400.

...
h1 = serial_open(pi, "/dev/ttyAMA0", 300)

h2 = serial_open(pi, "/dev/ttyUSB1", 19200, 0)
...
"""
function serial_open(self::Pi, tty, baud, ser_flags=0)
    # I p1 baud
    # I p2 ser_flags
    # I p3 len
    ## extension ##
    # s len data bytes
    return _u2i(_pigpio_command_ext(
        self.sl, _PI_CMD_SERO, baud, ser_flags, length(tty), [tty]))
end

"""
Closes the serial device associated with handle.

handle:= >=0 (as returned by a prior call to [*serial_open*]).

...
serial_close(pi, h1)
...
"""
function serial_close(self::Pi, handle)
    return _u2i(_pigpio_command(self.sl, _PI_CMD_SERC, handle, 0))
end

"""
Returns a single byte from the device associated with handle.

handle:= >=0 (as returned by a prior call to [*serial_open*]).

...
b = serial_read_byte(pi, h1)
...
"""
function serial_read_byte(self::Pi, handle)
    return _u2i(_pigpio_command(self.sl, _PI_CMD_SERRB, handle, 0))
end

"""
Writes a single byte to the device associated with handle.

handle:= >=0 (as returned by a prior call to [*serial_open*]).
byte_val:= 0-255, the value to write.

...
serial_write_byte(pi, h1, 23)

serial_write_byte(h1, ord(pi, 'Z'))
...
"""
function serial_write_byte(self::Pi, handle, byte_val)
    return _u2i(
        _pigpio_command(self.sl, _PI_CMD_SERWB, handle, byte_val))
end

"""
Reads up to count bytes from the device associated with handle.

handle:= >=0 (as returned by a prior call to [*serial_open*]).
count:= >0, the number of bytes to read.

The returned value is a tuple of the number of bytes read and a
bytearray containing the bytes.  If there was an error the
number of bytes read will be less than zero (and will contain
the error code).

...
(b, d) = serial_read(pi, h2, 100)
if b > 0
# process read data
...
"""
function serial_read(self::Pi, handle, count)
    # Don't raise exception.  Must release lock.
    bytes = u2i(
    _pigpio_command(self.sl, _PI_CMD_SERR, handle, count, false))
    if bytes > 0
        data = rxbuf(bytes)
    else
        data = ""
    end
    unlock(self.sl.l)
    return bytes, data
end

"""
Writes the data bytes to the device associated with handle.

handle:= >=0 (as returned by a prior call to [*serial_open*]).
data:= the bytes to write.

...
serial_write(pi, h1, b'\\x02\\x03\\x04')

serial_write(pi, h2, b'help')

serial_write(pi, h2, "hello")

serial_write(pi, h1, [2, 3, 4])
...
"""
function serial_write(self::Pi, handle, data)
    # I p1 handle
    # I p2 0
    # I p3 len
    ## extension ##
    # s len data bytes

    return _u2i(_pigpio_command_ext(
        self.sl, _PI_CMD_SERW, handle, 0, length(data), [data]))
end

"""
Returns the number of bytes available to be read from the
device associated with handle.

handle:= >=0 (as returned by a prior call to [*serial_open*]).

...
rdy = serial_data_available(pi, h1)

if rdy > 0
(b, d) = serial_read(pi, h1, rdy)
...
"""
function serial_data_available(self::Pi, handle)
    return _u2i(_pigpio_command(self.sl, _PI_CMD_SERDA, handle, 0))
end

"""
Opens a GPIO for bit bang reading of serial data.

user_gpio:= 0-31, the GPIO to use.
baud:= 50-250000, the baud rate.
bb_bits:= 1-32, the number of bits per word, default 8.

The serial data is held in a cyclic buffer and is read using
[*bb_serial_read*].

It is the caller's responsibility to read data from the cyclic
buffer in a timely fashion.

...
status = bb_serial_read_open(pi, 4, 19200)
status = bb_serial_read_open(pi, 17, 9600)
...
"""
function bb_serial_read_open(self, user_gpio, baud, bb_bits=8)
    # pigpio message format

    # I p1 user_gpio
    # I p2 baud
    # I p3 4
    ## extension ##
    # I bb_bits
    extents = IOBuffer()
    write(extents, bb_bits::Cuint)
    return _u2i(_pigpio_command_ext(
        self.sl, _PI_CMD_SLRO, user_gpio, baud, 4, extents))
end

"""
Returns data from the bit bang serial cyclic buffer.

user_gpio:= 0-31 (opened in a prior call to [*bb_serial_read_open*])

The returned value is a tuple of the number of bytes read and a
bytearray containing the bytes.  If there was an error the
number of bytes read will be less than zero (and will contain
the error code).

The bytes returned for each character depend upon the number of
data bits [*bb_bits*] specified in the [*bb_serial_read_open*]
command.

For [*bb_bits*] 1-8 there will be one byte per character.
For [*bb_bits*] 9-16 there will be two bytes per character.
For [*bb_bits*] 17-32 there will be four bytes per character.

...
(count, data) = bb_serial_read(pi, 4)
...
"""
function bb_serial_read(self, user_gpio)
    # Don't raise exception.  Must release lock.
    bytes = u2i(
        _pigpio_command(self.sl, _PI_CMD_SLR, user_gpio, 10000, false))
    if bytes > 0
        data = self._rxbuf(bytes)
    else
        data = ""
    end
    unlock(self.sl.l)
    return bytes, data
end

"""
Closes a GPIO for bit bang reading of serial data.

user_gpio:= 0-31 (opened in a prior call to [*bb_serial_read_open*])

...
status = bb_serial_read_close(pi, 17)
...
"""
function bb_serial_read_close(self, user_gpio)
    return _u2i(_pigpio_command(self.sl, _PI_CMD_SLRC, user_gpio, 0))
end

"""
Invert serial logic.

user_gpio:= 0-31 (opened in a prior call to [*bb_serial_read_open*])
invert:= 0-1 (1 invert, 0 normal)

...
status = bb_serial_invert(pi, 17, 1)
...
"""
function bb_serial_invert(self, user_gpio, invert)
    return _u2i(_pigpio_command(self.sl, _PI_CMD_SLRI, user_gpio, invert))
end


text: include("selfrecursiveinclude.jl")

ext: jl


text:
# 2520 is the smallest number that can be divided by each of the numbers from 1 to 10 without any remainder.

# What is the smallest positive number that is evenly divisible by all of the numbers from 1 to 20?

primes = [19, 17, 13, 11, 7, 5, 3, 2]
added_factors = [3, 2, 2, 2]

number = prod(primes)*prod(added_factors)

for i in 1:20
    if number % i != 0
        println("Not divisible by: ", i)
    end
end

max_stars_repo_name: pedvide/project-euler

max_issues_repo_name: pedvide/project-euler

max_forks_repo_name: pedvide/project-euler


text:
#### Woodbury matrices ####
# This type provides support for the Woodbury matrix identity
type Woodbury{T} <: AbstractMatrix{T}
    A
    U::Matrix{T}
    C
    Cp
    V::Matrix{T}
    tmpN1::Vector{T}
    tmpN2::Vector{T}
    tmpk1::Vector{T}
    tmpk2::Vector{T}

    function Woodbury(A::AbstractMatrix{T}, U::Matrix{T}, C, V::Matrix{T})
        N = size(A, 1)
        k = size(U, 2)
        if size(A, 2) != N || size(U, 1) != N || size(V, 1) != k || size(V, 2) != N
            throw(DimensionMismatch(""))
        end
        if k > 1
            if size(C, 1) != k || size(C, 2) != k
                throw(DimensionMismatch(""))
            end
        end
        Cp = inv(inv(C) + V*(A\U))
        # temporary space for allocation-free solver
        tmpN1 = Array(T, N)
        tmpN2 = Array(T, N)
        tmpk1 = Array(T, k)
        tmpk2 = Array(T, k)
        # don't copy A, it could be huge
        new(A, copy(U), copy(C), Cp, copy(V), tmpN1, tmpN2, tmpk1, tmpk2)
    end
end

Woodbury{T}(A::AbstractMatrix{T}, U::Matrix{T}, C, V::Matrix{T}) = Woodbury{T}(A, U, C, V)
Woodbury{T}(A::AbstractMatrix{T}, U::Vector{T}, C, V::Matrix{T}) = Woodbury{T}(A, reshape(U, length(U), 1), C, V)

size(W::Woodbury) = size(W.A)

function show(io::IO, W::Woodbury)
    println(io, summary(W), ":")
    print(io, "A:\n", W.A)
    print(io, "\nU:\n")
    print_matrix(io, W.U)
    if isa(W.C, Matrix)
        print(io, "\nC:\n")
        print_matrix(io, W.C)
    else
        print(io, "\nC: ", W.C)
    end
    print(io, "\nV:\n")
    print_matrix(io, W.V)
end

full{T}(W::Woodbury{T}) = convert(Matrix{T}, W)
convert{T}(::Type{Matrix{T}}, W::Woodbury{T}) = full(W.A) + W.U*W.C*W.V

function similar(W::Woodbury, T, dims::Dims)
    (length(dims) == 2 && dims[1] == dims[2] ) || throw(DimensionMismatch(""))
    n = size(W, 1)
    k = size(W.U, 2)
    Woodbury{T}(similar(W.A), Array(T, n, k), Array(T, k, k), Array(T, k, n))
end

copy(W::Woodbury) = Woodbury(W.A, W.U, W.C, W.V)

## Woodbury matrix routines ##

*(W::Woodbury, B::StridedVecOrMat)=W.A*B + W.U*(W.C*(W.V*B))

function \(W::Woodbury, R::StridedVecOrMat)
    AinvR = W.A\R
    AinvR - W.A\(W.U*(W.Cp*(W.V*AinvR)))
end

det(W::Woodbury)=det(W.A)*det(W.C)/det(W.Cp)

# Allocation-free solver for arbitrary strides (requires that W.A has a
# non-aliasing "solve" routine, e.g., is Tridiagonal)
function solve!(x::AbstractArray, xrng::Range{Int}, W::Woodbury, rhs::AbstractArray, rhsrng::Range{Int})
    solve!(W.tmpN1, 1:length(W.tmpN1), W.A, rhs, rhsrng)
    A_mul_B!(W.tmpk1, W.V, W.tmpN1)
    A_mul_B!(W.tmpk2, W.Cp, W.tmpk1)
    A_mul_B!(W.tmpN2, W.U, W.tmpk2)
    solve!(W.tmpN2, W.A, W.tmpN2)
    indx = first(xrng)
    xinc = step(xrng)
    for i = 1:length(W.tmpN2)
        x[indx] = W.tmpN1[i] - W.tmpN2[i]
        indx += xinc
    end
    nothing
end

function solve!(x::AbstractVector, W::Woodbury, rhs::AbstractVector)
    solve!(x, 1:length(x), W, rhs, 1:length(rhs))
    x
end
solve(W::Woodbury, rhs::AbstractVector) = solve!(similar(rhs), W, rhs)
solve(W::Woodbury, B::StridedMatrix)=solve!(similar(B), W, B)
function solve!(X::StridedMatrix, W::Woodbury, B::StridedMatrix)
    size(B, 1) == size(W, 1) || throw(DimensionMismatch(""))
    size(X) == size(B) || throw(DimensionMismatch(""))
    m, n = size(B)
    r = 1:m
    for j = 1:n
        r.start = (j-1)*m+1
        solve!(X, r, W, B, r)
    end
    X
end

max_stars_repo_name: alyst/julia

max_issues_repo_name: alyst/julia

max_forks_repo_name: alyst/julia


text:
function get_auto_traverse_strategy(com::CS.CoM)
    return com.sense == MOI.FEASIBILITY_SENSE ? :DFS : :BFS
end

function get_traverse_strategy(com; options = SolverOptions())
    if options.traverse_strategy == :DBFS
        return isempty(com.solutions) ? Val(:DFS) : Val(:BFS)
    end
    return Val(options.traverse_strategy)
end

function get_auto_branch_strategy(com::CS.CoM)
    return :IMPS # Infeasible and Minimum Possiblity Search
end

function get_branch_strategy(; options = SolverOptions())
    strategy = options.branch_strategy
    return Val(strategy)
end

function get_branch_split(; options = SolverOptions())
    strategy = options.branch_split
    return Val(strategy)
end

const POSSIBLE_OPTIONS = Dict(
    :traverse_strategy => [:Auto, :BFS, :DFS, :DBFS],
    :branch_strategy => [:Auto, :ABS, :IMPS],
    :branch_split => [:Auto, :Smallest, :Biggest, :InHalf],
)

function SolverOptions()
    logging = [:Info, :Table]
    table = TableSetup(
        [
            CS.TableCol(:open_nodes, "#Open", Int, 10, :center),
            CS.TableCol(:closed_nodes, "#Closed", Int, 10, :center),
            CS.TableCol(:incumbent, "Incumbent", Float64, 20, :center),
            CS.TableCol(:best_bound, "Best Bound", Float64, 20, :center),
            CS.TableCol(:duration, "Time [s]", Float64, 10, :center),
        ],
        Dict(:min_diff_duration => 5.0),
    )
    seed = 1
    traverse_strategy = :Auto
    branch_strategy = :Auto
    branch_split = :Auto
    backtrack = true
    max_bt_steps = typemax(Int)
    backtrack_sorting = true
    keep_logs = false
    rtol = 1e-6
    atol = 1e-6
    solution_type = Float64
    all_solutions = false
    all_optimal_solutions = false
    lp_optimizer = nothing
    time_limit = Inf
    no_prune = false
    decay = 0.999
    max_probes = 10
    max_confidence_deviation = 20
    simplify = true

    return SolverOptions(
        logging,
        table,
        time_limit,
        seed,
        traverse_strategy,
        branch_strategy,
        branch_split,
        backtrack,
        max_bt_steps,
        backtrack_sorting,
        keep_logs,
        rtol,
        atol,
        solution_type,
        all_solutions,
        all_optimal_solutions,
        lp_optimizer,
        no_prune,
        ActivityOptions(decay, max_probes, max_confidence_deviation),
        simplify,
    )
end

function combine_options(options)
    defaults = SolverOptions()
    options_dict = Dict{Symbol,Any}()
    for kv in options
        if !in(kv[1], fieldnames(SolverOptions))
            @error "The option $(kv[1]) doesn't exist."
        else
            moi_key = MOI.RawParameter(string(kv[1]))
            if is_possible_option_value(moi_key, kv[2])
                options_dict[kv[1]] = kv[2]
            else
                @error "The option $(kv[1]) doesn't have $(kv[2]) as a possible value. Possible values are: $(POSSIBLE_OPTIONS[moi_key])"
            end
        end
    end

    for fname in fieldnames(SolverOptions)
        if haskey(options_dict, fname)
            setfield!(
                defaults,
                fname,
                convert(fieldtype(SolverOptions, fname), options_dict[fname]),
            )
        end
    end
    return defaults
end

function is_possible_option_value(option_param::MOI.RawParameter, value)
    option = Symbol(option_param.name)
    if haskey(POSSIBLE_OPTIONS, option)
        return value in POSSIBLE_OPTIONS[option]
    end
    return true
end


text:
import Base: put!, wait, isready, take!

isready(D::MyChannel) = all(D.d)

function wait(D::MyChannel)
    while !isready(D)
        wait(D.cond_take)
    end
end

function put!(D::MyChannel, k::Int64)
    @inbounds D.d[k] = true
    if isready(D)
        notify(D.cond_take)
    else
        wait(D)
    end
    D
end

function take!(D::MyChannel)
    if isready(D)
        for i=1:length(D.d)
            @inbounds D.d[i]=false
        end
    end
end

max_stars_repo_name: paulmthompson/JNeuron

max_issues_repo_name: paulmthompson/JNeuron

max_forks_repo_name: paulmthompson/JNeuron


text:
using DrWatson
@quickactivate "ChaosThroughBilliards"
include(srcdir("style.jl"))

x = 1.0
y = 0.54

re = billiard_rectangle(x, y)
si = billiard_sinai(y/4, x, y)
st = billiard_stadium(x-y, y)

ps = [Particle(4x/5, y/5, π/9)]

# interactive_billiard(re, ps;
# backgroundcolor = :black, res = MAXRES, particle_size = 2.0, add_controls = false)

for (bd, name) in zip((re, si), ("rect", "sinai"))
# for (bd, name) in zip((si,), ("sinai",))

billiard_video(
    videodir("mirror_$(name).mp4"), bd, deepcopy(ps);
    frames = 1200, backgroundcolor = :black, colors = COLORS,
    tailwidth = 3.5, particle_size = 2, res = (1200*x, 1200*y),
    plot_particles = true,
    dt = 0.0001, speed = 150, tail = 10000, # this makes ultra fine temporal resolution
)

end

max_stars_repo_name: JuliaDynamics/ChaosThroughBilliards

max_issues_repo_name: JuliaDynamics/ChaosThroughBilliards

max_forks_repo_name: JuliaDynamics/ChaosThroughBilliards


text:
#Written By Johnathan Bizzano
module Optimizer
    using ..Utils
    import ..ConstantPool
    import ..Error
    using ..Computation

    @Fun(Problem, output::Float64, new_constants::AbstractArray)

    include("DeRand1Bin.jl")

    function copywraperror(cons, err_func::Error.WrappedFunc)
        Problem(
            function (constants)
                copy!(cons, constants)
                return err_func()
            end)
    end

    function de_rand_1_bin(; population_size, iterations, scalefactor = .8, crossoverrate = .7, rnginterval = 1E-8, IsVerbose=true)
        return MemoryWrapper(sett ->
                    function (cons_pool, error_func, lower_bound, upper_bound)
                        return diffevorandomsinglebin(sett.device, copywraperror(cons_pool, error_func), copy(cons_pool),
                                lower_bound, upper_bound, population_size,
                                iterations, scalefactor=scalefactor, crossoverrate=crossoverrate,
                                rnginterval=rnginterval, IsVerbose=IsVerbose)
                    end)
    end

end


text:
get_vars = Symbolics.get_variables

function isin(var::Num, vars)
    for varible in vars
        if isequal(varible, var)
            return true
        end
    end
    return false
end

function get_Num(vars::Dict)
    return collect(keys(vars))
end

function check_eqs(eqs)
    if typeof(eqs) == Equation
        return toexpr([eqs])
    elseif typeof(eqs) == Vector{Equation}
        return toexpr(eqs)
    else
        error("Error: type of eqs should not be $(typeof(eqs)).It must be Equation or Vector{Equation}")
    end
end

function toexpr(eqs)
    expr = Vector{Num}([])
    for eq in eqs
        push!(expr,eq.lhs - eq.rhs)
    end
    return expr
end

function get_all_vars(eqs_expr)
    vars = Vector{Num}([])
    for expr in eqs_expr
        vars = union(vars,get_vars(expr))
    end
    return vars
end

function check_vars(eqs, vars)
    copy_vars = copy(vars)
    eqs_expr = eqs
    eqs_var = get_all_vars(eqs_expr)
    if length(eqs_var) != length(collect(keys(vars)))
        error("Error: number of varialbes in vars don't equal to
        number of variables in eqs")
    end
    for var in eqs_var
        try
            pop!(copy_vars, var)
        catch
            error("Error: $(var) not in vars")
        end
    end
    rest_var = collect(keys(copy_vars))
    if length(rest_var) > 0
        error("$(rest_var) not in eqs")
    end
    return true
end

function islinear(eqs, vars)
    for eq in eqs
        if !Symbolics.islinear(eq.lhs, get_Num(vars))
            return false
        end
    end
    return true
end


text:
const Blend = CairoPattern

"""
    blend(from::Point, to::Point)

Create an empty linear blend.

A blend is a specification of how one color changes into another. Linear blends are
defined by two points: parallel lines through these points define the start and stop
locations of the blend. The blend is defined relative to the current axes origin. This means
that you should be aware of the current axes when you define blends, and when you use them.

To add colors, use `addstop()`.
"""
function blend(from::Point, to::Point)
    return Cairo.pattern_create_linear(from.x, from.y, to.x, to.y)
end

"""
    blend(centerpos1, rad1, centerpos2, rad2, color1, color2)

Create a radial blend.

Example:

    redblue = blend(
        pos, 0,                   # first circle center and radius
        pos, tiles.tilewidth/2,   # second circle center and radius
        "red",
        "blue"
        )
"""
function blend(centerpos1::Point, rad1, centerpos2::Point, rad2, color1, color2)
    newblend = blend(centerpos1, rad1, centerpos2, rad2)
    addstop(newblend, 0, color1)
    addstop(newblend, 1, color2)
    return newblend
end

"""
    blend(pt1::Point, pt2::Point, color1, color2)

Create a linear blend.

Example:

    redblue = blend(pos, pos, "red", "blue")
"""
function blend(pt1::Point, pt2::Point, color1, color2)
    newblend = blend(pt1, pt2)
    addstop(newblend, 0, color1)
    addstop(newblend, 1, color2)
    return newblend
end

"""
    blend(from::Point, startradius, to::Point, endradius)

Create an empty radial blend.

Radial blends are defined by two circles that define the start and stop locations. The
first point is the center of the start circle, the first radius is the radius of the first
circle.

A new blend is empty. To add colors, use `addstop()`.
"""
function blend(from::Point, startradius, to::Point, endradius)
    return Cairo.pattern_create_radial(from.x, from.y, startradius, to.x, to.y, endradius)
end

"""
    addstop(b::Blend, offset, col)
    addstop(b::Blend, offset, (r, g, b, a))
    addstop(b::Blend, offset, string)

Add a color stop to a blend. The offset specifies the location along the blend's 'control
vector', which varies between 0 (beginning of the blend) and 1 (end of the blend). For
linear blends, the control vector is from the start point to the end point. For radial
blends, the control vector is from any point on the start circle, to the corresponding point
on the end circle.

Examples:

    blendredblue = blend(Point(0, 0), 0, Point(0, 0), 1)
    addstop(blendredblue, 0, setcolor(sethue("red")..., .2))
    addstop(blendredblue, 1, setcolor(sethue("blue")..., .2))
    addstop(blendredblue, 0.5, sethue(randomhue()...))
    addstop(blendredblue, 0.5, setcolor(randomcolor()...))

"""
function addstop(b::Blend, offset, col::ColorTypes.Colorant)
    temp = convert(RGBA,  col)
    Cairo.pattern_add_color_stop_rgba(b, offset, temp.r, temp.g, temp.b, temp.alpha)
end

function addstop(b::Blend, offset, col::AbstractString)
    temp = parse(RGBA, col)
    Cairo.pattern_add_color_stop_rgba(b, offset, temp.r, temp.g, temp.b, temp.alpha)
end

function addstop(b::Blend, offset, col::NTuple{4, Number})
    Cairo.pattern_add_color_stop_rgba(b, offset, col[1], col[2], col[3], col[4])
end

function addstop(b::Blend, offset, col::NTuple{3, Number})
    currentopacity = get_current_alpha()
    Cairo.pattern_add_color_stop_rgba(b, offset, col[1], col[2], col[3], currentopacity)
end

"""
    setblend(blend)

Start using the named blend for filling graphics.

This aligns the original coordinates of the blend definition with the current axes.
"""
function setblend(b::Blend)
    Cairo.set_source(get_current_cr(), b)
end

"""
    blendadjust(ablend, center::Point, xscale, yscale, rot=0)

Modify an existing blend by scaling, translating, and rotating it so that it will fill a
shape properly even if the position of the shape is nowhere near the original location of
the blend's definition.

For example, if your blend definition was this (notice the `1`)

    blendgoldmagenta = blend(
            Point(0, 0), 0,                   # first circle center and radius
            Point(0, 0), 1,                   # second circle center and radius
            "gold",
            "magenta"
            )

you can use it in a shape that's 100 units across and centered at `pos`, by calling this:

    blendadjust(blendgoldmagenta, Point(pos.x, pos.y), 100, 100)

then use `setblend()`:

    setblend(blendgoldmagenta)
"""
function blendadjust(ablend, center::Point, xscale, yscale, rot=0)
    blendmatrix(ablend,
        juliatocairomatrix(
            rotationmatrix(-rot) *
            scalingmatrix(1/xscale, 1/yscale) *
            translationmatrix(-center.x, -center.y) *
        cairotojuliamatrix([1 0 0 1 0 0])))
end

"""
    blendmatrix(b::Blend, m)

Set the matrix of a blend.

To apply a sequence of matrix transforms to a blend:

```
A = [1 0 0 1 0 0]
Aj = cairotojuliamatrix(A)
Sj = scalingmatrix(2, .2) * Aj
Tj = translationmatrix(10, 0) * Sj
A1 = juliatocairomatrix(Tj)
blendmatrix(b, As)
```
"""
function blendmatrix(b::Blend, m)
    cm = Cairo.CairoMatrix(m[1], m[2], m[3], m[4], m[5], m[6])
    Cairo.set_matrix(b, cm)
end


text:
using PorousMaterialResampling
using CorrelationFunctions
using XUnit

@testset "Test PorousMaterialResampling" begin include("reconstruct.jl") end


text:
# Key methods to add for a new dimensional data type

"""
    rebuild(x, args...)
    rebuild(x; kw...)

Rebuild an object struct with updated field values.

This is an abstraction that alows inbuilt and custom types to be rebuilt
functionally to update them, as most objects in DimensionalData are immutable.

`x` can be a `AbstractDimArray`, a `Dimension`, `IndexMode` or other custom types.

The arguments reuired are defined for the abstract type that has a `rebuild` method.
"""
function rebuild end
rebuild(x; kw...) = ConstructionBase.setproperties(x, (; kw...))

"""
    dims(x) => Tuple{Vararg{<:Dimension}}
    dims(x, dims::Tuple) => Tuple{Vararg{<:Dimension}}
    dims(x, dim) => Dimension

Return a tuple of `Dimension`s for an object, in the order that matches
the axes or columns etc. of the underlying data.

`dims` can be `Dimension`, `Dimension` types, or `Symbols` for `Dim{Symbol}`.

The default is to return `nothing`.
"""
function dims end

"""
    refdims(x) => Tuple{Vararg{<:Dimension}}
    refdims(x, dims::Tuple) => Tuple{Vararg{<:Dimension}}
    refdims(x, dim) => Dimension

Reference dimensions for an array that is a slice or view of another
array with more dimensions.

`slicedims(a, dims)` returns a tuple containing the current new dimensions
and the new reference dimensions. Refdims can be stored in a field or disgarded,
as it is mostly to give context to plots. Ignoring refdims will simply leave some
captions empty.

The default is to return an empty `Tuple` `()`.
"""
function refdims end
refdims(x, lookup) = dims(refdims(x), lookup)

"""
    val(x)
    val(dims::Tuple) => Tuple

Return the contained value of a wrapper object.

`dims` can be `Dimension`, `Dimension` types, or `Symbols` for `Dim{Symbol}`.

Objects that don't define a `val` method are returned unaltered.
"""
function val end

"""
    index(dim::Dimension{<:Val}) => Tuple
    index(dim::Dimension{<:AbstractArray}) => AbstractArray
    index(dims::NTuple{N}) => Tuple{Vararg{Union{AbstractArray,Tuple},N}}

Return the contained index of a `Dimension`.

Only valid when a `Dimension` contains an `AbstractArray`
or a Val tuple like `Val{(:a, :b)}()`. The `Val` is unwrapped
to return just the `Tuple`

`dims` can be a `Dimension`, or a tuple of `Dimension`.
"""
function index end

"""
    mode(dim:Dimension) => IndexMode
    mode(dims::Tuple) => Tuple{Vararg{<:IndexMode,N}}
    mode(A::AbstractDimArray, [dims::Tuple]) => Tuple

Returns the [`IndexMode`](@ref) of a dimension. This dictates
properties of the dimension such as array axis and index order,
and sampling properties.

`dims` can be a `Dimension`, a dimension type, or a tuple of either.
"""
function mode end

"""
    metadata(dim::Dimension)
    metadata(dims::Tuple{<:Dimension,Vararg})
    metadata(A::AbstractDimArray, dims::Tuple)  => (Dim metadata)
    metadata(A::AbstractDimArray) => (Array metadata)

Returns the metadata for an array or the specified dimension(s).
`dims` can be a `Symbol` (with `Dim{X}`, a `Dimension`, a `Dimension` type,
or a mixed tuple.

`dims` can be a `Dimension`, a dimension type, or a tuple of either.
"""
function metadata end

"""
    bounds(dim::Dimension) => Tuple{T,T}}
    bounds(dims::Tuple{<:Dimension,Vararg}) => Tuple{Vararg{<:Tuple{T,T}}}
    bounds(A::AbstractArray, [dims]) => Tuple{Vararg{Tuple{T,T},N}}
    bounds(A::AbstractArray, dim) => Tuple{T,T}

Return the bounds of all dimensions of an object, of a specific dimension,
or of a tuple of dimensions.

`dims` can be a `Dimension`, a dimension type, or a tuple of either.
"""
function bounds end

"""
    name(x) => String
    name(xs::NTuple{N,<:Dimension}) => NTuple{N,String}
    name(A::AbstractDimArray, dims::NTuple{N}) => NTuple{N,String}

Get the name of an array or Dimension, or a tuple of of either.

`dims` can be `Dimension`s, `Dimension` types, or `Symbols` for `Dim{Symbol}`.
"""
function name end

"""
    units(x) => Union{Nothing,Any}
    units(::NTuple{N}) => NTuple{N}
    unit(A::AbstractDimArray, dims::NTuple{N}) => NTuple{N,String}

Get the units of an array or `Dimension`, or a tuple of of either.

Units do not have a set field, and may or may not be included in `metadata`.
This method is to facilitate use in labels and plots when units are available,
not a guarantee that they will be. If not available, `nothing` is returned.

`dims` can be `Dimension`s, `Dimension` types, or `Symbols` for `Dim{Symbol}`.
"""
function units end

"""
    label(x) => String
    label(dims::NTuple{N,<:Dimension}) => NTuple{N,String}
    label(A::AbstractDimArray, dims::NTuple{N,<:Dimension}) => NTuple{N,String}

Get a plot label for data or a dimension. This will include the name and units
if they exist, and anything else that should be shown on a plot.

`dims` can be `Dimension`s, `Dimension` types, or `Symbols` for `Dim{Symbol}`.
"""
function label end

"""
    order(dim:Dimension) => Order
    order(dims::Tuple) => Tuple{Vararg{<:Order,N}}
    order(A::AbstractDimArray, [dims::Tuple]) => Tuple

Return the [`Order`](@ref) for each dimension.

`dims` can be `Dimension`s, `Dimension` types, or `Symbols` for `Dim{Symbol}`.
"""
function order end

"""
    sampling(dim:Dimension) => Sampling
    sampling(dims::Tuple) => Tuple{Vararg{<:Sampling,N}}
    sampling(A::AbstractDimArray, [dims::Tuple]) => Tuple

Return the [`Sampling`](@ref) for each dimension.

`dims` can be `Dimension`s, `Dimension` types, or `Symbols` for `Dim{Symbol}`.
"""
function sampling end

"""
    span(dim:Dimension) => Span
    span(dims::Tuple) => Tuple{Vararg{<:Span,N}}
    span(A::AbstractDimArray, [dims::Tuple]) => Tuple

Return the [`Span`](@ref) for each dimension.

`dims` can be `Dimension`s, `Dimension` types, or `Symbols` for `Dim{Symbol}`.
"""
function span end

"""
    locus(dim:Dimension) => Locus
    locus(dims::Tuple) => Tuple{Vararg{<:Locus,N}}
    locus(A::AbstractDimArray, [dims::Tuple]) => Tuple

Return the [`Locus`](@ref) for each dimension.

`dims` can be `Dimension`s, `Dimension` types, or `Symbols` for `Dim{Symbol}`.
"""
function locus end

"""
    arrayorder(dim:Dimension) => Union{Forward,Reverse}
    arrayorder(dims::Tuple) => Tuple{Vararg{<:Union{Forward,Reverse},N}}
    arrayorder(A::AbstractDimArray, [dims::Tuple]) => Tuple

Return the [`Order`](@ref) (`Forward` or `Reverse`) of the array,
for each dimension.

`dims` can be `Dimension`s, `Dimension` types, or `Symbols` for `Dim{Symbol}`.
"""
function arrayorder end

"""
    indexorder(dim:Dimension) => Union{Forward,Reverse}
    indexorder(dims::Tuple) => Tuple{Vararg{<:Union{Forward,Reverse},N}}
    indexorder(A::AbstractDimArray, [dims::Tuple]) => Tuple

Return the [`Order`](@ref) (`Forward` or `Reverse`) of the dimension index,
for each dimension.

`dims` can be `Dimension`s, `Dimension` types, or `Symbols` for `Dim{Symbol}`.
"""
function indexorder end

"""
    relation(dim:Dimension) => Union{Forward,Reverse}
    relation(dims::Tuple) => Tuple{Vararg{<:Union{Forward,Reverse},N}}
    relation(A::AbstractDimArray, [dims::Tuple]) => Tuple

Return the relation (`Forward` or `Reverse`) between the dimension index
and the array axis, for each dimension.

`dims` can be `Dimension`s, `Dimension` types, or `Symbols` for `Dim{Symbol}`.
"""
function relation end


text:
using Iris
using Test

@testset "BoundaryConditions" begin
    bc = noBC{1}()
    @test BCHermiticity(bc)<:HermitianBC
    @test Iris.Common.BCLocality(bc)<:Iris.Common.LocalBC
    @test Iris.Common.getside(bc) == 1
    bc = DirichletBC{2}()
    @test BCHermiticity(bc)<:HermitianBC
    @test Iris.Common.BCLocality(bc)<:Iris.Common.LocalBC
    @test Iris.Common.getside(bc) == 2
    bc = NeumannBC{3}()
    @test BCHermiticity(bc)<:HermitianBC
    @test Iris.Common.BCLocality(bc)<:Iris.Common.LocalBC
    @test Iris.Common.getside(bc) == 3
    bc = FloquetBC{4}()
    @test BCHermiticity(bc)<:HermitianBC
    @test Iris.Common.BCLocality(bc)<:Iris.Common.NonLocalBC
    @test Iris.Common.getside(bc) == 4
    bc = MatchedBC{5}(in=[1])
    @test BCHermiticity(bc)<:NonHermitianBC
    @test Iris.Common.BCLocality(bc)<:Iris.Common.NonLocalBC
    @test Iris.Common.getside(bc) == 5
    conj(bc).out == bc.in
end


text:
"""
Abstract evolutionary optimization result type
"""
abstract type OptimizationResults end

"""
    summary(result)

Shows the optimization algorithm that produced this `result`.
"""
summary(or::OptimizationResults) = summary(or.method)

"""
    minimizer(result)

A minimizer object from the optimization `result`.
"""
minimizer(r::OptimizationResults) = r.minimizer

"""
    minimum(result)

A minimum value from the optimization `result`.
"""
minimum(r::OptimizationResults) = r.minimum

"""
    iterations(result)

A number of iterations to reach the minimum.
"""
iterations(r::OptimizationResults) = r.iterations

"""
    iteration_limit_reached(result)

Returns `true` if the iteration limit was reached.
"""
iteration_limit_reached(r::OptimizationResults) = r.iteration_converged

"""
    trace(result)

Returns a trace of optimization states from the optimization `result`.
"""
trace(r::OptimizationResults) = length(r.trace) > 0 ? r.trace : error("No trace in optimization results. To get a trace, run optimize() with store_trace = true.")

"""
    f_calls(result)

Returns a number of an objective function calls.
"""
f_calls(r::OptimizationResults) = r.f_calls

"""
    tol(result)

Returns an absolute tollerance value of the optimization `result`.
"""
tol(r::OptimizationResults) = error("`tol` is not implemented for $(summary(r)).")

time_limit(r::OptimizationResults) = r.time_limit
time_run(  r::OptimizationResults) = r.time_run

"""
Evolutionary optimization result type
"""
mutable struct EvolutionaryOptimizationResults{O<:AbstractOptimizer, T, Tx, Tf} <: OptimizationResults
    method::O
    minimizer::Tx
    minimum::Tf
    iterations::Int
    iteration_converged::Bool
    converged::Bool
    abstol::T
    trace::OptimizationTrace
    f_calls::Int
    time_limit::Float64
    time_run::Float64
end

"""
    converged(result)

Returns `true` if the optimization sucesfully coverged to a minimum value.
"""
converged(r::EvolutionaryOptimizationResults) = r.converged

tol(r::EvolutionaryOptimizationResults) = r.abstol

function show(io::IO, r::EvolutionaryOptimizationResults)
    failure_string = "failure"
    if iteration_limit_reached(r)
        failure_string *= " (reached maximum number of iterations)"
    end
    if time_run(r) > time_limit(r)
        failure_string *= " (exceeded time limit of $(time_limit(r)))"
    end
    print(io, "\n")
    print(io, " * Status: ", converged(r) ? "success" : failure_string, "\n\n")
    print(io, " * Candidate solution\n")
    mzr = minimizer(r)
    if isa(mzr, Array)
        nx = length(mzr)
        str_x_elements = ["$_x" for _x in Iterators.take(mzr, min(nx, 3))]
        if nx >= 4
            push!(str_x_elements, " ...")
        end
        print(io, "    Minimizer:  [", join(str_x_elements, ", "),  "]\n")
    else
        print(io, "    Minimizer:  $mzr\n")
    end
    print(io, "    Minimum:    $(minimum(r))\n")
    print(io, "    Iterations: $(iterations(r))\n")
    print(io, "\n")
    print(io, " * Found with\n")
    print(io, "    Algorithm: $(summary(r))\n")
    print(io, "\n")
    print(io, " * Work counters\n")
    tr = round(time_run(r); digits=4)
    tl = isnan(time_limit(r)) ? Inf : round(time_limit(r); digits=4)
    print(io, "    Seconds run:   $tr (vs limit $tl)\n")
    print(io, "    Iterations:    $(iterations(r))\n")
    print(io, "    f(x) calls:    $(f_calls(r))\n")
    return
end


text:
using Parameters

struct DSGE_Model{data <: AbstractArray, alfa <: Float64, bet <: Float64, phi <:Float64, var_pars <: Array{Float64}}
    "Data"
    data::data
    "Alpha"
    alfa::alfa
    "Beta"
    bet::bet
    "Phi"
    phi::phi
    "VAR EUA"
    var_pars::var_pars
end

function (problem::DSGE_Model)(pars)
    @unpack alfa,bet,phi,var_pars, data = problem
    @unpack h,sig,eta,delta_D,theta_D,delta_I,theta_I, chi, rho_a, rho_gamma, rho_cp, rho_phi, sigma_a,sigma_gamma, sigma_cp,sigma_phi, rho_1, rho_2, lambda_pi, lambda_y, lambda_s, sigma_i = pars
    llh = log_like_dsge([alfa;bet;phi;var_pars;h;sig;eta;delta_D;theta_D;delta_I;theta_I; chi; rho_a; rho_gamma; rho_cp; rho_phi; sigma_a;sigma_gamma; sigma_cp;sigma_phi; rho_1; rho_2; lambda_pi; lambda_y; lambda_s; sigma_i],data)
        return llh+logpdf(prior_h,h)+logpdf(prior_sig,sig)+logpdf(prior_eta,eta)+logpdf(prior_delta_D,delta_D)+logpdf(prior_theta_D,theta_D)+logpdf(prior_delta_I,delta_I)+logpdf(prior_theta_I,theta_I)+logpdf(prior_chi,chi)+logpdf(prior_rho_a,rho_a)+logpdf(prior_rho_gamma,rho_gamma)+logpdf(prior_rho_cp,rho_cp)+logpdf(prior_rho_phi,rho_phi)+logpdf(prior_sigma_a,sigma_a)+logpdf(prior_sigma_gamma,sigma_gamma)+logpdf(prior_sigma_cp,sigma_cp)+logpdf(prior_sigma_phi,sigma_phi)+logpdf(prior_rho_1,rho_1)+logpdf(prior_rho_2,rho_2)+logpdf(prior_lambda_pi,lambda_pi)+logpdf(prior_lambda_y,lambda_y)+logpdf(prior_lambda_s,lambda_s)+logpdf(prior_sigma_i,sigma_i)
end

function problem_transform(p::DSGE_Model)
    as((h = as_unit_interval, sig = as_positive_real, eta = as_positive_real, delta_D = as_unit_interval, theta_D = as_unit_interval, delta_I = as_unit_interval, theta_I = as_unit_interval, chi = as_positive_real, rho_a = as_unit_interval, rho_gamma = as_unit_interval, rho_cp = as_unit_interval, rho_phi = as_unit_interval, sigma_a = as_positive_real, sigma_gamma = as_positive_real, sigma_cp = as_positive_real, sigma_phi = as_positive_real, rho_1 = as_unit_interval,rho_2 = as_unit_interval,lambda_pi = as_positive_real,lambda_y = as_positive_real,lambda_s = as_positive_real,sigma_i = as_positive_real,))
end

max_stars_repo_name: danmrc/DSGE-HMC

max_issues_repo_name: danmrc/DSGE-HMC

max_forks_repo_name: danmrc/DSGE-HMC


text:
module Precompile_TypeDB_Protocol

using ProtoBuf
using ..TypeDBClient

g = TypeDBClient.typedb.protocol

function _precompile_()
    ccall(:jl_generating_output, Cint, ()) == 1 || return nothing
    @assert Base.precompile(Tuple{typeof(meta),Type{g.Relation_RemovePlayer_Req}})   # time: 0.11327494
    @assert Base.precompile(Tuple{typeof(meta),Type{g.Thing_SetHas_Req}})   # time: 0.04059317
    @assert Base.precompile(Tuple{typeof(meta),Type{g.Type_ResPart}})   # time: 0.03636061
    @assert Base.precompile(Tuple{typeof(meta),Type{g.QueryManager_ResPart}})   # time: 0.030095141
    @assert Base.precompile(Tuple{typeof(meta),Type{g.ConceptManager_Res}})   # time: 0.02858994
    @assert Base.precompile(Tuple{typeof(meta),Type{g.QueryManager_Res}})   # time: 0.027888825
    @assert Base.precompile(Tuple{typeof(meta),Type{g.Transaction_Client}})   # time: 0.024184799
    @assert Base.precompile(Tuple{typeof(meta),Type{g.ThingType_GetOwns_Req}})   # time: 0.023652244
    @assert Base.precompile(Tuple{typeof(meta),Type{g.ConceptMapGroup}})   # time: 0.023117358
    @assert Base.precompile(Tuple{typeof(meta),Type{g.Thing_GetHas_ResPart}})   # time: 0.023006914
    @assert Base.precompile(Tuple{typeof(meta),Type{g.Type_SetSupertype_Res}})   # time: 0.022904683
    @assert Base.precompile(Tuple{typeof(meta),Type{g.Thing_Req}})   # time: 0.02200334
    @assert Base.precompile(Tuple{typeof(meta),Type{g.QueryManager_Undefine_Req}})   # time: 0.021988356
    @assert Base.precompile(Tuple{typeof(meta),Type{g.ConceptManager_GetThing_Req}})   # time: 0.021817593
    @assert Base.precompile(Tuple{typeof(meta),Type{g.AttributeType_Put_Req}})   # time: 0.021352673
    @assert Base.precompile(Tuple{typeof(meta),Type{g.LogicManager_PutRule_Req}})   # time: 0.021125427
    @assert Base.precompile(Tuple{typeof(meta),Type{g.AttributeType_Get_Res}})   # time: 0.021080498
    @assert Base.precompile(Tuple{typeof(meta),Type{g.QueryManager_Delete_Req}})   # time: 0.020728838
    @assert Base.precompile(Tuple{typeof(meta),Type{g.QueryManager_Req}})   # time: 0.02057298
    @assert Base.precompile(Tuple{typeof(meta),Type{g.AttributeType_GetOwners_Req}})   # time: 0.020349976
    @assert Base.precompile(Tuple{typeof(meta),Type{g.EntityType_Create_Req}})   # time: 0.01956494
    @assert Base.precompile(Tuple{typeof(meta),Type{g.Thing_Res}})   # time: 0.017860284
    @assert Base.precompile(Tuple{Core.kwftype(typeof(Type)),NamedTuple{(:reqs,), Tuple{Vector{g.Transaction_Req}}},Type{g.Transaction_Client}})   # time: 0.017820457
    @assert Base.precompile(Tuple{typeof(meta),Type{g.Thing_ResPart}})   # time: 0.017225703
    @assert Base.precompile(Tuple{typeof(meta),Type{g.ConceptManager_Req}})   # time: 0.015448299
    @assert Base.precompile(Tuple{typeof(meta),Type{g.Explanation}})   # time: 0.015126186
    @assert Base.precompile(Tuple{typeof(meta),Type{g._Type}})   # time: 0.014537138
    @assert Base.precompile(Tuple{typeof(meta),Type{g.Thing}})   # time: 0.013407678
    @assert Base.precompile(Tuple{typeof(meta),Type{g.Attribute_Value}})   # time: 0.012723284
    @assert Base.precompile(Tuple{typeof(meta),Type{g.RelationType_SetRelates_Req}})   # time: 0.012556498
    @assert Base.precompile(Tuple{typeof(meta),Type{g.ThingType_SetOwns_Req}})   # time: 0.01229427
    @assert Base.precompile(Tuple{typeof(meta),Type{g.Type_SetSupertype_Req}})   # time: 0.012243458
    @assert Base.precompile(Tuple{typeof(meta),Type{g.ConceptMap}})   # time: 0.012237117
    @assert Base.precompile(Tuple{typeof(meta),Type{g.Rule_Delete_Req}})   # time: 0.011915595
    @assert Base.precompile(Tuple{typeof(meta),Type{g.Type_SetLabel_Req}})   # time: 0.011890562
    @assert Base.precompile(Tuple{typeof(meta),Type{g.LogicManager_Req}})   # time: 0.011888661
    @assert Base.precompile(Tuple{typeof(meta),Type{g.ConceptManager_GetThingType_Res}})   # time: 0.011882486
    @assert Base.precompile(Tuple{typeof(meta),Type{g.Explanation_VarList}})   # time: 0.011782851
    @assert Base.precompile(Tuple{typeof(meta),Type{g.ThingType_GetPlays_ResPart}})   # time: 0.011676633
    @assert Base.precompile(Tuple{typeof(meta),Type{g.Relation_AddPlayer_Req}})   # time: 0.011581464
    @assert Base.precompile(Tuple{typeof(meta),Type{g.LogicManager_GetRules_ResPart}})   # time: 0.011531885
    @assert Base.precompile(Tuple{Type{g.Transaction_Server}})   # time: 0.011474489
    @assert Base.precompile(Tuple{typeof(meta),Type{g.QueryManager_MatchGroup_Req}})   # time: 0.011434047
    @assert Base.precompile(Tuple{typeof(meta),Type{g.QueryManager_MatchAggregate_Res}})   # time: 0.011292823
    @assert Base.precompile(Tuple{typeof(meta),Type{g.ThingType_GetInstances_ResPart}})   # time: 0.011253795
    @assert Base.precompile(Tuple{typeof(meta),Type{g.Explainables}})   # time: 0.011244047
    @assert Base.precompile(Tuple{typeof(meta),Type{g.QueryManager_Update_ResPart}})   # time: 0.011241352
    @assert Base.precompile(Tuple{typeof(meta),Type{g.Attribute_GetOwners_Req}})   # time: 0.011203552
    @assert Base.precompile(Tuple{typeof(meta),Type{g.Type_GetSubtypes_ResPart}})   # time: 0.011081024
    @assert Base.precompile(Tuple{Type{g.Session_Pulse_Res}})   # time: 0.011076559
    @assert Base.precompile(Tuple{typeof(meta),Type{g.QueryManager_Match_ResPart}})   # time: 0.011035354
    @assert Base.precompile(Tuple{typeof(meta),Type{g.RelationType_Create_Res}})   # time: 0.011021406
    @assert Base.precompile(Tuple{typeof(meta),Type{g.ThingType_SetPlays_Req}})   # time: 0.010988181
    @assert Base.precompile(Tuple{typeof(meta),Type{g.Transaction_Stream_ResPart}})   # time: 0.010976053
    @assert Base.precompile(Tuple{typeof(meta),Type{g.RelationType_GetRelatesForRoleLabel_Req}})   # time: 0.010974543
    @assert Base.precompile(Tuple{typeof(meta),Type{g.Transaction_Rollback_Res}})   # time: 0.010892521
    @assert Base.precompile(Tuple{typeof(meta),Type{g.NumericGroup}})   # time: 0.010849708
    @assert Base.precompile(Tuple{typeof(meta),Type{g.Relation_GetPlayers_ResPart}})   # time: 0.010761713
    @assert Base.precompile(Tuple{typeof(meta),Type{g.ConceptManager_GetThingType_Req}})   # time: 0.010752866
    @assert Base.precompile(Tuple{typeof(meta),Type{g.Rule_Req}})   # time: 0.010743545
    @assert Base.precompile(Tuple{typeof(meta),Type{g.Explainable}})   # time: 0.010725368
    @assert Base.precompile(Tuple{typeof(meta),Type{g.LogicManager_Res}})   # time: 0.010704533
    @assert Base.precompile(Tuple{typeof(meta),Type{g.QueryManager_MatchGroupAggregate_Req}})   # time: 0.010661521
    @assert Base.precompile(Tuple{typeof(meta),Type{g.EntityType_Create_Res}})   # time: 0.010562908
    @assert Base.precompile(Tuple{typeof(meta),Type{g.Type_IsAbstract_Req}})   # time: 0.010539572
    @assert Base.precompile(Tuple{typeof(meta),Type{g.Rule_Res}})   # time: 0.010471803
    @assert Base.precompile(Tuple{typeof(meta),Type{g.RelationType_UnsetRelates_Res}})   # time: 0.010455552
    @assert Base.precompile(Tuple{typeof(meta),Type{g.QueryManager_Update_Req}})   # time: 0.010436815
    @assert Base.precompile(Tuple{typeof(meta),Type{g.RelationType_GetRelates_ResPart}})   # time: 0.010421858
    @assert Base.precompile(Tuple{typeof(meta),Type{g.LogicManager_GetRule_Res}})   # time: 0.010417575
    @assert Base.precompile(Tuple{Type{g.Session_Open_Res}})   # time: 0.010406066
    @assert Base.precompile(Tuple{typeof(meta),Type{g.LogicManager_GetRule_Req}})   # time: 0.010353709
    @assert Base.precompile(Tuple{typeof(meta),Type{g.ConceptManager_PutAttributeType_Req}})   # time: 0.010277795
    @assert Base.precompile(Tuple{typeof(meta),Type{g.ThingType_UnsetPlays_Req}})   # time: 0.010258366
    @assert Base.precompile(Tuple{typeof(meta),Type{g.Rule}})   # time: 0.010229583
    @assert Base.precompile(Tuple{typeof(meta),Type{g.Thing_GetPlaying_ResPart}})   # time: 0.010186644
    @assert Base.precompile(Tuple{typeof(meta),Type{g.LogicManager_ResPart}})   # time: 0.01015848
    @assert Base.precompile(Tuple{typeof(meta),Type{g.QueryManager_MatchGroupAggregate_ResPart}})   # time: 0.010028499
    @assert Base.precompile(Tuple{typeof(meta),Type{g.QueryManager_Insert_ResPart}})   # time: 0.010002855
    @assert Base.precompile(Tuple{typeof(meta),Type{g.ConceptManager_PutAttributeType_Res}})   # time: 0.009968781
    @assert Base.precompile(Tuple{typeof(meta),Type{g.QueryManager_Undefine_Res}})   # time: 0.009943785
    @assert Base.precompile(Tuple{typeof(meta),Type{g.QueryManager_Explain_Req}})   # time: 0.0099278
    @assert Base.precompile(Tuple{typeof(meta),Type{g.ThingType_UnsetAbstract_Req}})   # time: 0.009927308
    @assert Base.precompile(Tuple{typeof(meta),Type{g.Thing_GetHas_Req}})   # time: 0.009924366
    @assert Base.precompile(Tuple{typeof(meta),Type{g.RoleType_GetRelationTypes_ResPart}})   # time: 0.009915132
    @assert Base.precompile(Tuple{typeof(meta),Type{g.Thing_GetType_Req}})   # time: 0.009914403
    @assert Base.precompile(Tuple{typeof(meta),Type{g.Attribute_GetOwners_ResPart}})   # time: 0.009897542
    @assert Base.precompile(Tuple{typeof(meta),Type{g.Type_GetSupertypes_Req}})   # time: 0.009879283
    @assert Base.precompile(Tuple{typeof(meta),Type{g.AttributeType_SetRegex_Req}})   # time: 0.009803332
    @assert Base.precompile(Tuple{typeof(meta),Type{g.Type_IsAbstract_Res}})   # time: 0.009788234
    @assert Base.precompile(Tuple{typeof(meta),Type{g.Numeric}})   # time: 0.009760826
    @assert Base.precompile(Tuple{typeof(meta),Type{g.ConceptManager_PutEntityType_Req}})   # time: 0.009753091
    @assert Base.precompile(Tuple{typeof(meta),Type{g.QueryManager_MatchGroup_ResPart}})   # time: 0.009750567
    @assert Base.precompile(Tuple{typeof(meta),Type{g.RelationType_GetRelatesForRoleLabel_Res}})   # time: 0.009745792
    @assert Base.precompile(Tuple{typeof(meta),Type{g.AttributeType_GetRegex_Res}})   # time: 0.009720303
    @assert Base.precompile(Tuple{typeof(meta),Type{g.Thing_GetRelations_Req}})   # time: 0.009717738
    @assert Base.precompile(Tuple{typeof(meta),Type{g.AttributeType_Get_Req}})   # time: 0.009703181
    @assert Base.precompile(Tuple{typeof(meta),Type{g.Relation_GetPlayersByRoleType_ResPart}})   # time: 0.009663366
    @assert Base.precompile(Tuple{typeof(meta),Type{g.Transaction_Commit_Req}})   # time: 0.009637717
    @assert Base.precompile(Tuple{typeof(meta),Type{g.RoleType_GetPlayers_Req}})   # time: 0.009627365
    @assert Base.precompile(Tuple{typeof(meta),Type{g.ThingType_GetOwns_ResPart}})   # time: 0.009579233
    @assert Base.precompile(Tuple{typeof(meta),Type{g.Thing_GetRelations_ResPart}})   # time: 0.009575339
    @assert Base.precompile(Tuple{typeof(meta),Type{g.Type_Delete_Req}})   # time: 0.00953032
    @assert Base.precompile(Tuple{typeof(meta),Type{g.RelationType_UnsetRelates_Req}})   # time: 0.009518828
    @assert Base.precompile(Tuple{typeof(meta),Type{g.Relation_GetPlayersByRoleType_RoleTypeWithPlayer}})   # time: 0.009517489
    @assert Base.precompile(Tuple{typeof(meta),Type{g.Relation_GetRelating_ResPart}})   # time: 0.009454011
    @assert Base.precompile(Tuple{typeof(meta),Type{g.Type_GetSupertype_Res}})   # time: 0.009450192
    @assert Base.precompile(Tuple{typeof(meta),Type{g.ConceptManager_PutRelationType_Req}})   # time: 0.009445112
    @assert Base.precompile(Tuple{typeof(meta),Type{g.QueryManager_Explain_ResPart}})   # time: 0.009409025
    @assert Base.precompile(Tuple{typeof(meta),Type{g.Transaction_Commit_Res}})   # time: 0.009379997
    @assert Base.precompile(Tuple{typeof(meta),Type{g.Type_GetSupertype_Req}})   # time: 0.009365733
    @assert Base.precompile(Tuple{typeof(meta),Type{g.RoleType_GetPlayers_ResPart}})   # time: 0.009362556
    @assert Base.precompile(Tuple{typeof(meta),Type{g.ThingType_UnsetOwns_Req}})   # time: 0.009352733
    @assert Base.precompile(Tuple{typeof(meta),Type{g.QueryManager_MatchAggregate_Req}})   # time: 0.009337294
    @assert Base.precompile(Tuple{typeof(meta),Type{g.AttributeType_GetOwners_ResPart}})   # time: 0.009325083
    @assert Base.precompile(Tuple{typeof(meta),Type{g.Relation_GetPlayersByRoleType_Req}})   # time: 0.009319417
    @assert Base.precompile(Tuple{typeof(meta),Type{g.QueryManager_Match_Req}})   # time: 0.009304448
    @assert Base.precompile(Tuple{typeof(meta),Type{g.ConceptManager_GetThing_Res}})   # time: 0.009292919
    @assert Base.precompile(Tuple{typeof(meta),Type{g.Type_GetSupertypes_ResPart}})   # time: 0.009276615
    @assert Base.precompile(Tuple{typeof(meta),Type{g.Explainables_Owned}})   # time: 0.009266806
    @assert Base.precompile(Tuple{typeof(meta),Type{g.ConceptManager_PutRelationType_Res}})   # time: 0.009264035
    @assert Base.precompile(Tuple{typeof(meta),Type{g.Relation_GetPlayers_Req}})   # time: 0.009262358
    @assert Base.precompile(Tuple{typeof(meta),Type{g.ThingType_SetOwns_Res}})   # time: 0.009207241
    @assert Base.precompile(Tuple{typeof(meta),Type{g.RelationType_Create_Req}})   # time: 0.009176729
    @assert Base.precompile(Tuple{typeof(meta),Type{g.ConceptManager_PutEntityType_Res}})   # time: 0.009174211
    @assert Base.precompile(Tuple{typeof(meta),Type{g.ThingType_GetPlays_Req}})   # time: 0.00915192
    @assert Base.precompile(Tuple{typeof(meta),Type{g.LogicManager_GetRules_Req}})   # time: 0.009127014
    @assert Base.precompile(Tuple{typeof(meta),Type{g.Transaction_Rollback_Req}})   # time: 0.009115244
    @assert Base.precompile(Tuple{typeof(meta),Type{g.Thing_GetType_Res}})   # time: 0.009109693
    @assert Base.precompile(Tuple{typeof(meta),Type{g.Concept}})   # time: 0.009107701
    @assert Base.precompile(Tuple{typeof(meta),Type{g.QueryManager_Define_Req}})   # time: 0.00910414
    @assert Base.precompile(Tuple{typeof(meta),Type{g.Rule_SetLabel_Req}})   # time: 0.009058628
    @assert Base.precompile(Tuple{typeof(meta),Type{g.ThingType_SetPlays_Res}})   # time: 0.009030789
    @assert Base.precompile(Tuple{typeof(meta),Type{g.QueryManager_Insert_Req}})   # time: 0.008967989
    @assert Base.precompile(Tuple{typeof(meta),Type{g.LogicManager_PutRule_Res}})   # time: 0.008959993
    @assert Base.precompile(Tuple{typeof(meta),Type{g.QueryManager_Define_Res}})   # time: 0.008932782
    @assert Base.precompile(Tuple{typeof(meta),Type{g.Thing_UnsetHas_Req}})   # time: 0.008894681
    @assert Base.precompile(Tuple{typeof(meta),Type{g.AttributeType_Put_Res}})   # time: 0.008858936
    @assert Base.precompile(Tuple{typeof(meta),Type{g.ThingType_GetInstances_Req}})   # time: 0.008833896
    @assert Base.precompile(Tuple{typeof(meta),Type{g.Transaction_Open_Res}})   # time: 0.008832352
    @assert Base.precompile(Tuple{typeof(meta),Type{g.Thing_Delete_Res}})   # time: 0.008705167
    @assert Base.precompile(Tuple{typeof(meta),Type{g.ThingType_UnsetPlays_Res}})   # time: 0.008704196
    @assert Base.precompile(Tuple{typeof(meta),Type{g.ThingType_SetAbstract_Res}})   # time: 0.008673189
    @assert Base.precompile(Tuple{typeof(meta),Type{g.Relation_RemovePlayer_Res}})   # time: 0.008617327
    @assert Base.precompile(Tuple{typeof(meta),Type{g.AttributeType_SetRegex_Res}})   # time: 0.008589685
    @assert Base.precompile(Tuple{typeof(meta),Type{g.RoleType_GetRelationTypes_Req}})   # time: 0.008556931
    @assert Base.precompile(Tuple{typeof(meta),Type{g.ThingType_UnsetAbstract_Res}})   # time: 0.008529544
    @assert Base.precompile(Tuple{typeof(meta),Type{g.Thing_SetHas_Res}})   # time: 0.008495999
    @assert Base.precompile(Tuple{typeof(meta),Type{g.Type_SetLabel_Res}})   # time: 0.008488814
    @assert Base.precompile(Tuple{typeof(meta),Type{g.Type_Delete_Res}})   # time: 0.008479692
    @assert Base.precompile(Tuple{typeof(meta),Type{g.Type_GetSubtypes_Req}})   # time: 0.008458761
    @assert Base.precompile(Tuple{typeof(meta),Type{g.Relation_AddPlayer_Res}})   # time: 0.008457315
    @assert Base.precompile(Tuple{typeof(meta),Type{g.ThingType_UnsetOwns_Res}})   # time: 0.008415489
    @assert Base.precompile(Tuple{typeof(meta),Type{g.Rule_Delete_Res}})   # time: 0.00839782
    @assert Base.precompile(Tuple{typeof(meta),Type{g.Thing_UnsetHas_Res}})   # time: 0.008370067
    @assert Base.precompile(Tuple{typeof(meta),Type{g.RelationType_GetRelates_Req}})   # time: 0.008363739
    @assert Base.precompile(Tuple{typeof(meta),Type{g.Thing_Delete_Req}})   # time: 0.00835378
    @assert Base.precompile(Tuple{typeof(meta),Type{g.Thing_GetPlaying_Req}})   # time: 0.008339361
    @assert Base.precompile(Tuple{typeof(meta),Type{g.Rule_SetLabel_Res}})   # time: 0.008334654
    @assert Base.precompile(Tuple{typeof(meta),Type{g.Relation_GetRelating_Req}})   # time: 0.00832416
    @assert Base.precompile(Tuple{typeof(meta),Type{g.ThingType_SetAbstract_Req}})   # time: 0.008243179
    @assert Base.precompile(Tuple{typeof(meta),Type{g.AttributeType_GetRegex_Req}})   # time: 0.008231799
    @assert Base.precompile(Tuple{typeof(meta),Type{g.RelationType_SetRelates_Res}})   # time: 0.008213344
    @assert Base.precompile(Tuple{typeof(meta),Type{g.QueryManager_Delete_Res}})   # time: 0.008168623
    @assert Base.precompile(Tuple{Type{g.Transaction_Res}})   # time: 0.001181481
    @assert Base.precompile(Tuple{Type{g.Transaction_Open_Res}}) # time: 0.001030904

    @info "Precompile protocol done!"
end

end


text:
### A Pluto.jl notebook ###
# v0.14.4

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ 42b4e6fe-a897-11eb-2676-f7fd96f35a22
begin
	using NeidSolarScripts
	using RvSpectMLBase
	using EchelleInstruments
	using PlutoUI
	using Plots
end

# ╔═╡ 5313b674-a897-11eb-3820-218f26a14d4d
begin
	fits_path = "/mnt/data_simons/NEID/DRPv0.7-fixedflatfielding2/"
	fits_fn = "neidL1_20210104T182500.fits"
	continuum_afs_path = joinpath(fits_path,"output","continuum")
	continuum_afs_fn = "neidL1_20210104T182500_continuum=afs.jld2"

	spec = NEID.read_data(joinpath(fits_path,fits_fn))
end;

# ╔═╡ 28aced5a-3100-4523-9f1a-f7b14016791b
md"Order index  $(@bind order_idx NumberField(1:100; default=60))
"

# ╔═╡ 8352fbc6-483c-41fa-8649-1d3d40a33ef6
md"   Stretch factor $(@bind stretch_fac NumberField(0.05:20; default=1.0))
   Merging threshold $(@bind merge_thresh NumberField(0:1; default=0.05))
   Ready $(@bind calc_continuum CheckBox(default=false))
"

# ╔═╡ f36dee02-9a27-4a58-938f-5bad6a9ba286
md"Stretch Factor = $stretch_fac   Merging Factor = $merge_thresh"

# ╔═╡ 37b20ff3-e93a-4ce1-82eb-29c23581add3
if 1<=order_idx<=100
	 lambda = spec.λ[:,order_idx]
	 f_obs = spec.flux[:,order_idx]
	 lambda_plt = lambda
	 λmin_default = minimum(lambda)
	 λmax_default = maximum(lambda)
end;

# ╔═╡ 851d39ba-a897-11eb-28e5-1191c9c9e172
if calc_continuum && 0<stretch_fac<=40 && 0 <= merge_thresh <10
	 (anchors, continuum, f_filtered_internal) = Continuum.calc_continuum(lambda,f_obs,
	      stretch_factor=stretch_fac, merging_threshold=merge_thresh, verbose=true)
	 println("# num_anchors = ", length(anchors))
	 #f_filtered = Continuum.calc_rolling_median(lambda,f_obs, width=4)
end;

# ╔═╡ f99f1957-0799-4296-a681-874b7407bb58
if calc_continuum
	f_filtered = f_filtered_internal
	continuum_plt = continuum
	anchors_plt = anchors
else
	f_filtered = Continuum.calc_rolling_median(lambda,f_obs, width=4)
	continuum_plt = nothing
	anchors_plt = nothing
end;

# ╔═╡ abbc4880-fedd-4bdd-a1db-b94e8a4d5654
@bind λmin NumberField(4000:10000; default=λmin_default)

# ╔═╡ c12d7cda-68be-4b78-b3e0-26a93cb6ac39
@bind λmax NumberField(4000:10000; default=λmax_default)

# ╔═╡ ce053df6-a899-11eb-2024-a9d147acc990
#= md"λ_min  $(@bind λmin NumberField(4000:10000; default=λmin_default))
   λ_max  $(@bind λmax NumberField(4000:10000; default=λmax_default))
"
=#

# ╔═╡ ab0fc462-a897-11eb-1727-e71788522c02
if 10000>λmax>λmin>4000
  plt = scatter(lambda,f_obs, ms=1.2, color=:grey,legend=:none)
	if f_filtered != nothing && anchors_plt != nothing
  		plot!(lambda_plt,f_filtered,color=:cyan)
  		plot!(lambda,continuum_plt,color=:magenta)
  		scatter!(lambda[anchors_plt], f_filtered[anchors_plt], color=:blue, ms=4)
	end
  #ylims!(0.95,1.15)
  #xmin = 5365
  #xmax = xmin+35
  xlims!(plt,λmin,λmax)
  ylims!(plt,extrema(f_filtered[λmin.<=lambda.<=λmax]))
  plt
end

# ╔═╡ Cell order:
# ╠═42b4e6fe-a897-11eb-2676-f7fd96f35a22
# ╠═5313b674-a897-11eb-3820-218f26a14d4d
# ╟─28aced5a-3100-4523-9f1a-f7b14016791b
# ╠═8352fbc6-483c-41fa-8649-1d3d40a33ef6
# ╟─f36dee02-9a27-4a58-938f-5bad6a9ba286
# ╟─37b20ff3-e93a-4ce1-82eb-29c23581add3
# ╟─851d39ba-a897-11eb-28e5-1191c9c9e172
# ╟─f99f1957-0799-4296-a681-874b7407bb58
# ╠═abbc4880-fedd-4bdd-a1db-b94e8a4d5654
# ╠═c12d7cda-68be-4b78-b3e0-26a93cb6ac39
# ╟─ce053df6-a899-11eb-2024-a9d147acc990
# ╠═ab0fc462-a897-11eb-1727-e71788522c02


text:
# Problem evaluation #
# ========================== #
function _eval_first_stage(stochasticprogram::StochasticProgram, x::AbstractVector)
    first_stage = get_stage_one(stochasticprogram)
    return eval_objective(first_stage.obj, x)
end
function _eval_second_stage(stochasticprogram::TwoStageStochasticProgram, x::AbstractVector, scenario::AbstractScenario, solver::MPB.AbstractMathProgSolver)
    outcome = outcome_model(stochasticprogram, x, scenario, solver)
    solve(outcome)
    return probability(scenario)*getobjectivevalue(outcome)
end
function _eval_second_stages(stochasticprogram::TwoStageStochasticProgram{S,SP},
                             x::AbstractVector,
                             solver::MPB.AbstractMathProgSolver) where {S, SP <: ScenarioProblems}
    outcome_generator = scenario -> outcome_model(stochasticprogram, x, scenario; solver = solver)
    return outcome_mean(outcome_generator, scenarios(stochasticprogram))
end
function _eval_second_stages(stochasticprogram::TwoStageStochasticProgram{S,SP},
                             x::AbstractVector,
                             solver::MPB.AbstractMathProgSolver) where {S, SP <: DScenarioProblems}
    Qs = Vector{Float64}(undef, nworkers())
    outcome_generator = scenario -> outcome_model(stochasticprogram, x, scenario; solver = solver)
    @sync begin
        for (i,w) in enumerate(workers())
            @async Qs[i] = remotecall_fetch((sp,outcome_generator)->begin
                scenarioproblems = fetch(sp)
                isempty(scenarioproblems.scenarios) && return 0.0
                return outcome_mean(outcome_generator, scenarioproblems.scenarios)
            end,
            w,
            stochasticprogram.scenarioproblems[w-1],
            outcome_generator)
        end
    end
    return sum(Qs)
end
function _stat_eval_second_stages(stochasticprogram::TwoStageStochasticProgram{S,SP},
                                  x::AbstractVector,
                                  solver::MPB.AbstractMathProgSolver) where {S, SP <: ScenarioProblems}
    outcome_generator = scenario -> outcome_model(stochasticprogram, x, scenario; solver = solver)
    𝔼Q, σ² = welford(outcome_generator, scenarios(stochasticprogram))
    return 𝔼Q, sqrt(σ²)
end
function _stat_eval_second_stages(stochasticprogram::TwoStageStochasticProgram{S,SP},
                                  x::AbstractVector,
                                  solver::MPB.AbstractMathProgSolver) where {S, SP <: DScenarioProblems}
    partial_welfords = Vector{Tuple{Float64,Float64,Int}}(undef, nworkers())
    @sync begin
        for (i,w) in enumerate(workers())
            @async partial_welfords[i] = remotecall_fetch((sp,stage_one_generator,stage_two_generator,stage_one_params,stage_two_params,x,solver)->begin
                scenarioproblems = fetch(sp)
                isempty(scenarioproblems.scenarios) && return zero(eltype(x)), zero(eltype(x))
                    outcome_generator = scenario -> begin
                        outcome_model = Model(solver = solver)
                        _outcome_model!(outcome_model,
                                        stage_one_generator,
                                        stage_two_generator,
                                        stage_one_params,
                                        stage_two_params,
                                        x,
                                        scenario)
                        return outcome_model
                    end
                return (welford(outcome_generator, scenarioproblems.scenarios)..., length(scenarioproblems.scenarios))
            end,
            w,
            stochasticprogram.scenarioproblems[w-1],
            stochasticprogram.generator[:stage_1_vars],
            stochasticprogram.generator[:stage_2],
            stage_parameters(stochasticprogram, 1),
            stage_parameters(stochasticprogram, 2),
            x,
            solver)
        end
    end
    𝔼Q, σ², _ = reduce(aggregate_welford, partial_welfords)
    return 𝔼Q, sqrt(σ²)
end
function _eval(stochasticprogram::StochasticProgram{2}, x::AbstractVector, solver::MPB.AbstractMathProgSolver)
    xlength = decision_length(stochasticprogram)
    length(x) == xlength || error("Incorrect length of given decision vector, has ", length(x), " should be ", xlength)
    all(.!(isnan.(x))) || error("Given decision vector has NaN elements")
    cᵀx = _eval_first_stage(stochasticprogram, x)
    𝔼Q = _eval_second_stages(stochasticprogram, x, solver)
    return cᵀx+𝔼Q
end
# Mean/variance calculations #
# ========================== #
function outcome_mean(outcome_generator::Function, scenarios::Vector{<:AbstractScenario})
    Qs = zeros(length(scenarios))
    for (i,scenario) in enumerate(scenarios)
        let outcome = outcome_generator(scenario)
            status = solve(outcome, suppress_warnings = true)
            if status != :Optimal
                if status == :Infeasible
                    Qs[i] = outcome.objSense == :Max ? -Inf : Inf
                elseif status == :Unbounded
                    Qs[i] = outcome.objSense == :Max ? Inf : -Inf
                else
                    error("Outcome model could not be solved, returned status: $status")
                end
            else
                Qs[i] = probability(scenario)*getobjectivevalue(outcome)
            end
        end
    end
    return sum(Qs)
end
function welford(generator::Function, scenarios::Vector{<:AbstractScenario})
    Q̄ₖ = 0
    Sₖ = 0
    N = length(scenarios)
    for k = 1:N
        Q̄ₖ₋₁ = Q̄ₖ
        let problem = generator(scenarios[k])
            status = solve(problem, suppress_warnings = true)
            Q = if status != :Optimal
                Q = if status == :Infeasible
                    problem.objSense == :Max ? -Inf : Inf
                elseif status == :Unbounded
                    problem.objSense == :Max ? Inf : -Inf
                else
                    error("Outcome model could not be solved, returned status: $status")
                end
            else
                Q = getobjectivevalue(problem)
            end
            Q̄ₖ = Q̄ₖ + (Q-Q̄ₖ)/k
            Sₖ = Sₖ + (Q-Q̄ₖ)*(Q-Q̄ₖ₋₁)
        end
    end
    return Q̄ₖ, Sₖ/(N-1)
end
function aggregate_welford(left::Tuple, right::Tuple)
    x̄ₗ, σₗ², nₗ = left
    x̄ᵣ, σᵣ², nᵣ = right
    δ = x̄ᵣ-x̄ₗ
    N = nₗ+nᵣ
    x̄ = (nₗ*x̄ₗ+nᵣ*x̄ᵣ)/N
    Sₗ = σₗ²*(nₗ-1)
    Sᵣ = σᵣ²*(nᵣ-1)
    S = Sₗ+Sᵣ+nₗ*nᵣ/N*δ^2
    return (x̄, S/(N-1), N)
end
# Evaluation API #
# ========================== #
"""
    evaluate_decision(stochasticprogram::TwoStageStochasticProgram,
                      decision::AbstractVector;
                      solver = JuMP.UnsetSolver())

Evaluate the first-stage `decision` in `stochasticprogram`.

In other words, evaluate the first-stage objective at `decision` and solve outcome models of `decision` for every available scenario. Optionally, supply a capable `solver` to solve the outcome models. Otherwise, any previously set solver will be used.
"""
function evaluate_decision(stochasticprogram::StochasticProgram{2}, decision::AbstractVector; solver::SPSolverType = JuMP.UnsetSolver())
    # Use cached solver if available
    supplied_solver = pick_solver(stochasticprogram, solver)
    # Abort if no solver was given
    if isa(supplied_solver, JuMP.UnsetSolver)
        error("Cannot evaluate decision without a solver.")
    end
    return _eval(stochasticprogram, decision, internal_solver(supplied_solver))
end
"""
    evaluate_decision(stochasticprogram::TwoStageStochasticProgram,
                      decision::AbstractVector,
                      scenario::AbstractScenario;
                      solver = JuMP.UnsetSolver())

Evaluate the result of taking the first-stage `decision` if `scenario` is the actual outcome in `stochasticprogram`.
"""
function evaluate_decision(stochasticprogram::StochasticProgram{2},
                           decision::AbstractVector,
                           scenario::AbstractScenario;
                           solver::SPSolverType = JuMP.UnsetSolver())
    # Use cached solver if available
    supplied_solver = pick_solver(stochasticprogram, solver)
    # Abort if no solver was given
    if isa(supplied_solver, JuMP.UnsetSolver)
        error("Cannot evaluate decision without a solver.")
    end
    outcome = outcome_model(stochasticprogram, decision, scenario; solver = solver)
    status = solve(outcome)
    if status == :Optimal
        return _eval_first_stage(stochasticprogram, decision) + getobjectivevalue(outcome)
    end
    error("Outcome model could not be solved, returned status: $status")
end
"""
    evaluate_decision(stochasticmodel::StochasticModel{2},
                      decision::AbstractVector,
                      sampler::AbstractSampler;
                      solver = JuMP.UnsetSolver(),
                      confidence = 0.95,
                      N = 1000)

Return a statistical estimate of the objective of the two-stage `stochasticmodel` at `decision` in the form of a confidence interval at level `confidence`, over the scenario distribution induced by `sampler`.

In other words, evaluate `decision` on a sampled model of size `N`. Generate an confidence interval using the sample variance of the evaluation.

See also: [`confidence_interval`](@ref)
"""
function evaluate_decision(stochasticmodel::StochasticModel{2},
                           decision::AbstractVector,
                           sampler::AbstractSampler;
                           solver::SPSolverType = JuMP.UnsetSolver(),
                           confidence::AbstractFloat = 0.95,
                           Ñ::Integer = 1000,
                           kw...)
    CI = let eval_model = sample(stochasticmodel, sampler, Ñ, defer = true; kw...)
        # Condidence level
        α = 1-confidence
        cᵀx = _eval_first_stage(eval_model, decision)
        𝔼Q, σ = _stat_eval_second_stages(eval_model, decision, internal_solver(solver))
        z = quantile(Normal(0,1), 1-α)
        L = cᵀx + 𝔼Q - z*σ/sqrt(Ñ)
        U = cᵀx + 𝔼Q + z*σ/sqrt(Ñ)
        remove_scenarios!(eval_model)
        return ConfidenceInterval(L, U, confidence)
    end
    return CI
end
"""
    lower_bound(stochasticmodel::StochasticModel{2},
                sampler::AbstractSampler;
                solver = JuMP.UnsetSolver(),
                confidence = 0.95,
                N = 100,
                M = 10)

Generate a confidence interval around a lower bound on the true optimum of the two-stage `stochasticmodel` at level `confidence`, over the scenario distribution induced by `sampler`.

`N` is the size of the sampled models used to generate the interval and generally governs how tight it is. `M` is the number of sampled models.
"""
function lower_bound(stochasticmodel::StochasticModel{2}, sampler::AbstractSampler; solver::SPSolverType = JuMP.UnsetSolver(), confidence::AbstractFloat = 0.95, N::Integer = 100, M::Integer = 10, log = true, keep = true, offset = 0, indent::Int = 0, kw...)
    # Condidence level
    α = 1-confidence
    # Lower bound
    Qs = Vector{Float64}(undef, M)
    progress = Progress(M, 0.0, "$(repeat(" ", indent))Lower CI    ")
    log && sleep(0.1)
    log && ProgressMeter.update!(progress, 0, keep = false, offset = offset)
    for i = 1:M
        let sampled_model = sample(stochasticmodel, sampler, N; kw...)
            Qs[i] = VRP(sampled_model, solver = solver)
            remove_scenarios!(sampled_model)
            remove_subproblems!(sampled_model)
        end
        log && ProgressMeter.update!(progress, i, keep = keep, offset = offset)
    end
    Q̂ = mean(Qs)
    σ = std(Qs)
    t = quantile(TDist(M-1), 1-α)
    L = Q̂ - t*σ/sqrt(M)
    U = Q̂ + t*σ/sqrt(M)
    return ConfidenceInterval(L, U, 1-α)
end
"""
    upper_bound(stochasticmodel::StochasticModel{2},
                sampler::AbstractSampler;
                solver = JuMP.UnsetSolver(),
                confidence = 0.95,
                N = 100,
                T = 10,
                n = 1000)

Generate a confidence interval around an upper of the true optimum of the two-stage `stochasticmodel` at level `confidence`, over the scenario distribution induced by `sampler`.

`N` is the size of the sampled model used to generate a candidate decision. `Ñ` is the size of each sampled model and `T` is the number of sampled models.
"""
function upper_bound(stochasticmodel::StochasticModel{2}, sampler::AbstractSampler; solver::SPSolverType = JuMP.UnsetSolver(), confidence::AbstractFloat = 0.95, N::Integer = 100, T::Integer = 10, Ñ::Integer = 1000, log = true, keep = true, offset = 0, indent::Int = 0, kw...)
    # Condidence level
    α = 1-confidence
    # decision generation
    sampled_model = sample(stochasticmodel, sampler, N; kw...)
    optimize!(sampled_model, solver = solver)
    x̂ = optimal_decision(sampled_model)
    return upper_bound(stochasticmodel, x̂, sampler; solver = solver, confidence = confidence, T = T, Ñ = Ñ, log = log, keep = keep, offset = offset, indent = indent, kw...)
end
"""
    upper_bound(stochasticmodel::StochasticModel{2},
                x::AbstractVector,
                sampler::AbstractSampler;
                solver = JuMP.UnsetSolver(),
                confidence = 0.95,
                T = 10,
                Ñ = 1000)

Generate a confidence interval around an upper bound of the expected value of the decision `x` in the two-stage `stochasticmodel` at level `confidence`, over the scenario distribution induced by `sampler`.

`Ñ` is the size of each sampled model and `T` is the number of sampled models.
"""
function upper_bound(stochasticmodel::StochasticModel{2}, x::AbstractVector, sampler::AbstractSampler; solver::SPSolverType = JuMP.UnsetSolver(), confidence::AbstractFloat = 0.95, T::Integer = 10, Ñ::Integer = 1000, log = true, keep = true, offset = 0, indent::Int = 0, kw...)
    # Condidence level
    α = 1-confidence
    Qs = Vector{Float64}(undef, T)
    progress = Progress(T, 0.0, "$(repeat(" ", indent))Upper CI    ")
    log && sleep(0.1)
    log && ProgressMeter.update!(progress, 0, keep = false, offset = offset)
    for i = 1:T
        let eval_model = sample(stochasticmodel, sampler, Ñ, defer = true; kw...)
            Qs[i] = evaluate_decision(eval_model, x; solver = solver)
            remove_scenarios!(eval_model)
        end
        log && ProgressMeter.update!(progress, i, keep = keep, offset = offset)
    end
    Q̂ = mean(Qs)
    σ = std(Qs)
    t = quantile(TDist(T-1), 1-α)
    L = Q̂ - t*σ/sqrt(T)
    U = Q̂ + t*σ/sqrt(T)
    return ConfidenceInterval(L, U, 1-α)
end
"""
    confidence_interval(stochasticmodel::StochasticModel{2},
                        sampler::AbstractSampler;
                        solver = JuMP.UnsetSolver(),
                        confidence = 0.9,
                        N = 100,
                        M = 10,
                        T = 10)

Generate a confidence interval around the true optimum of the two-stage `stochasticmodel` at level `confidence` using SAA, over the scenario distribution induced by `sampler`.

`N` is the size of the sampled models used to generate the interval and generally governs how tight it is. `M` is the number of sampled models used in the lower bound calculation, and `T` is the number of sampled models used in the upper bound calculation.
"""
function confidence_interval(stochasticmodel::StochasticModel{2}, sampler::AbstractSampler; solver::SPSolverType = JuMP.UnsetSolver(), confidence::AbstractFloat = 0.9, N::Integer = 100, M::Integer = 10, T::Integer = 10, Ñ::Integer = 1000, log = true, keep = true, offset = 0, indent::Int = 0, kw...)
    # Condidence level
    α = (1-confidence)/2
    # Lower bound
    lower_CI = lower_bound(stochasticmodel, sampler; solver = solver, confidence = 1-α, N = N, M = M, log = log, keep = keep, offset = offset, indent = indent, kw...)
    L = lower(lower_CI)
    # Upper bound
    upper_CI = upper_bound(stochasticmodel, sampler; solver = solver, confidence = 1-α, N = N, T = T, Ñ = Ñ, log = log, keep = keep, offset = offset, indent = indent, kw...)
    U = upper(upper_CI)
    return ConfidenceInterval(L, U, confidence)
end
"""
    gap(stochasticmodel::StochasticModel{2},
        x::AbstractVector,
        sampler::AbstractSampler;
        solver = JuMP.UnsetSolver(),
        confidence = 0.9,
        N = 100,
        M = 10,
        T = 10)

Generate a confidence interval around the gap between the result of using decison `x` and true optimum of the two-stage `stochasticmodel` at level `confidence` using SAA, over the scenario distribution induced by `sampler`.

`N` is the size of the SAA models used to generate the interval and generally governs how tight it is. `M` is the number of sampled models used in the lower bound calculation, and `T` is the number of sampled models used in the upper bound calculation.
"""
function gap(stochasticmodel::StochasticModel{2}, x::AbstractVector, sampler::AbstractSampler; solver::SPSolverType = JuMP.UnsetSolver(), confidence::AbstractFloat = 0.9, N::Integer = 100, M::Integer = 10, T::Integer = 10, Ñ::Integer = 1000, log = true, keep = true, offset = 0, indent::Int = 0, kw...)
    # Condidence level
    α = (1-confidence)/2
    # Lower bound
    lower_CI = lower_bound(stochasticmodel, sampler; solver = solver, confidence = 1-α, N = N, M = M, log = log, keep = keep, offset = offset, indent = indent, kw...)
    L = lower(lower_CI)
    # Upper bound
    upper_CI = upper_bound(stochasticmodel, x, sampler; solver = solver, confidence = 1-α, N = N, T = T, Ñ = Ñ, log = log, keep = keep, offset = offset, indent = indent, kw...)
    U = upper(upper_CI)
    return ConfidenceInterval(0., U-L, confidence)
end
# ========================== #


text:
using BenchmarkTools
using PkgBenchmark

using TextStylometry

# Define a parent BenchmarkGroup to contain our suite
const SUITE = BenchmarkGroup()

SUITE["SimpleDocument"] = BenchmarkGroup(["english", "japanese"])
SUITE["ComplexityDocument"] = BenchmarkGroup(["english", "japanese"])

test_file = "test_corpus.txt"
SUITE["SimpleDocument"]["file"] = @benchmarkable SimpleDocument($test_file)
SUITE["SimpleDocument"]["tokens"] = @benchmarkable tokens(SimpleDocument($test_file))


text:
using PyPlot

# Global parameters
K_CONST =  9.869232667160130e-16 * 86400 * 1e3
K = 20.0 .* ones(n,m) # millidarcy
K[8:10,:] .= 120.0

K *= K_CONST
n = 15
m = 30
h = 30.0 # meter

x = (1:m)*h .- 0.5h
y = (1:n)*h .- 0.5h
pcolormesh(x, y, K, rasterized=true)
plot([3*h], [9*h], ">r")
plot([28*h], [9*h], "<g")
plot(x, 0.45h*ones(m), "v")
gca().invert_yaxis()
xlabel("x")
ylabel("y")
colorbar()
axis("scaled")
gca().set(frame_on=false)
savefig("data/setting.png")
savefig("data/setting.pdf")

max_stars_repo_name: ADCMEMarket/ADCMEImages

max_issues_repo_name: ADCMEMarket/ADCMEImages

max_forks_repo_name: ADCMEMarket/ADCMEImages


text:
function read_cnf(filename)

    lines = filter(line -> !startswith(line, "c"), readlines(filename))

    if !startswith(lines[1], "p cnf")
        error("Incorrect format: First line must be `p cnf`")
    end

    num_vars, num_clauses = parse.(Int, split(lines[1])[3:4])

    clauses = Vector{Int}[]

    line = 2

    for i in 1:num_clauses
        clause = Int[]
        done = false

        while true
            current_literals = parse.(Int, split(lines[line]))
            if current_literals[end] == 0
                done = true
                pop!(current_literals)
            end

            line += 1

            append!(clause, current_literals)

            done && break
        end

        push!(clauses, clause)

    end

    return SATProblem(clauses)
end


text:
using ADCME
using PyCall
using LinearAlgebra
using PyPlot
using Random
using AdFem
Random.seed!(233)

function fem_stiffness(hmat,m,n,h)
    fem_stiffness_ = load_op_and_grad("./build/libFemStiffness","fem_stiffness", multiple=true)
    hmat,m_,n_,h = convert_to_tensor([hmat,m,n,h], [Float64,Int32,Int32,Float64])
    ii, jj, vv = fem_stiffness_(hmat,m_,n_,h)
    SparseTensor(ii, jj, vv, 2(m+1)*(n+1), 2(m+1)*(n+1))
    # ii, jj, vv
end

# TODO: specify your input parameters
m = n = 20
h = 0.1
D = diagm(0=>[1,1,0.5])
A = fem_stiffness(D,m,n,h)
Av = compute_fem_stiffness_matrix(D, m, n, h)
sess = Session(); init(sess)
@show maximum(abs.(run(sess, A) - Av))
# run(sess, A[1])

# uncomment it for testing gradients
# error()

# TODO: change your test parameter to `m`
#       in the case of `multiple=true`, you also need to specify which component you are testings
# gradient check -- v
fix_rhs = constant(rand(2(m+1)*(n+1)))
function scalar_function(hmat)
    res = find(fem_stiffness(hmat,m,n,h))[3]
    return sum(res^2)
end

# TODO: change `m_` and `v_` to appropriate values
m_ = constant(diagm(0=>[1,1,0.5]))
v_ = rand(3,3)
y_ = scalar_function(m_)
dy_ = gradients(y_, m_)
ms_ = Array{Any}(undef, 5)
ys_ = Array{Any}(undef, 5)
s_ = Array{Any}(undef, 5)
w_ = Array{Any}(undef, 5)
gs_ =  @. 1 / 10^(1:5)

for i = 1:5
    g_ = gs_[i]
    ms_[i] = m_ + g_*v_
    ys_[i] = scalar_function(ms_[i])
    s_[i] = ys_[i] - y_
    w_[i] = s_[i] - g_*sum(v_.*dy_)
end

sess = Session(); init(sess)
sval_ = run(sess, s_)
wval_ = run(sess, w_)
close("all")
loglog(gs_, abs.(sval_), "*-", label="finite difference")
loglog(gs_, abs.(wval_), "+-", label="automatic differentiation")
loglog(gs_, gs_.^2 * 0.5*abs(wval_[1])/gs_[1]^2, "--",label="\$\\mathcal{O}(\\gamma^2)\$")
loglog(gs_, gs_ * 0.5*abs(sval_[1])/gs_[1], "--",label="\$\\mathcal{O}(\\gamma)\$")

plt.gca().invert_xaxis()
legend()
xlabel("\$\\gamma\$")
ylabel("Error")


text:
# Generate notebooks using literate
using Pkg
Pkg.instantiate()

using SIIPExamples
using PowerSystemCaseBuilder

PowerSystemCaseBuilder.clear_all_serialized_system()
SIIPExamples.literate_all(execute = false) # ensures that all notebooks exist
SIIPExamples.literate_all(force = true)


text:
using IntervalArithmetic
using Test

@test IntervalArithmetic.round_expr(:(a + b), RoundDown) == :(($(Expr(:escape, :+)))($(Expr(:escape, :a)), $(Expr(:escape, :b)), RoundingMode{:Down}()))

@test IntervalArithmetic.round_expr(:(sin(a)), RoundUp) ==
    :(sin($(Expr(:escape, :a)), RoundingMode{:Up}()))


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

JLLWrappers.@generate_main_file_header("MPSolve")
JLLWrappers.@generate_main_file("MPSolve", UUID("184ad7ad-f978-5c9a-833c-fb4ab31fd3fb"))
end  # module MPSolve_jll


text:
using FixedPointNumbers
using Base.Test
using Compat

@test reinterpret(0xa2uf8)  == 0xa2
@test reinterpret(0xa2uf10) == 0xa2
@test reinterpret(0xa2uf12) == 0xa2
@test reinterpret(0xa2uf14) == 0xa2
@test reinterpret(0xa2uf16) == 0xa2

@test reinterpret(UFixed8, 0xa2) == 0xa2uf8
@test reinterpret(UFixed10, 0x1fa2) == 0x1fa2uf10
@test reinterpret(UFixed12, 0x1fa2) == 0x1fa2uf12
@test reinterpret(UFixed14, 0x1fa2) == 0x1fa2uf14
@test reinterpret(UFixed16, 0x1fa2) == 0x1fa2uf16

@test UFixed8(1.0) == 0xffuf8
@test UFixed8(0.5) == 0x80uf8
@test UFixed14(1.0) == 0x3fffuf14
v = @compat UFixed12.([2])
@test v == UFixed12[0x1ffeuf12]
@test isa(v, Vector{UFixed12})

UF2 = (UFixed{UInt32,16}, UFixed{UInt64,3}, UFixed{UInt64,51}, UFixed{UInt128,7}, UFixed{UInt128,51})

for T in (FixedPointNumbers.UF..., UF2...)
    @test zero(T) == 0
    @test one(T) == 1
    @test one(T) * one(T) == one(T)
    @test typemin(T) == 0
    @test realmin(T) == 0
    @test eps(zero(T)) == eps(typemax(T))
    @test sizeof(T) == sizeof(FixedPointNumbers.rawtype(T))
end
@test typemax(UFixed8) == 1
@test typemax(UFixed10) == typemax(UInt16)//(2^10-1)
@test typemax(UFixed12) == typemax(UInt16)//(2^12-1)
@test typemax(UFixed14) == typemax(UInt16)//(2^14-1)
@test typemax(UFixed16) == 1
@test typemax(UFixed10) == typemax(UInt16) // (2^10-1)
@test typemax(UFixed12) == typemax(UInt16) // (2^12-1)
@test typemax(UFixed14) == typemax(UInt16) // (2^14-1)
@test typemax(UFixed{UInt32,16}) == typemax(UInt32) // (2^16-1)
@test typemax(UFixed{UInt64,3}) == typemax(UInt64) // (2^3-1)
@test typemax(UFixed{UInt128,7}) == typemax(UInt128) // (2^7-1)
@test typemax(UFixed{UInt128,100}) == typemax(UInt128) // (UInt128(2)^100-1)

# TODO: change back to InexactError when it allows message strings
@test_throws ArgumentError UFixed8(2)
@test_throws ArgumentError UFixed8(255)
@test_throws ArgumentError UFixed8(0xff)
@test_throws ArgumentError UFixed16(2)
@test_throws ArgumentError UFixed16(0xff)
@test_throws ArgumentError UFixed16(0xffff)
@test_throws ArgumentError convert(UFixed8,  typemax(UFixed10))
@test_throws ArgumentError convert(UFixed16, typemax(UFixed10))
@test_throws ArgumentError convert(UFixed{UInt128,100}, 10^9)
@test_throws ArgumentError convert(UFixed{UInt128,100}, 10.0^9)

x = UFixed8(0.5)
@test isfinite(x) == true
@test isnan(x) == false
@test isinf(x) == false

@test convert(UFixed8,  1.1/typemax(UInt8)) == eps(UFixed8)
@test convert(UFixed10, 1.1/typemax(UInt16)*64) == eps(UFixed10)
@test convert(UFixed12, 1.1/typemax(UInt16)*16) == eps(UFixed12)
@test convert(UFixed14, 1.1/typemax(UInt16)*4)  == eps(UFixed14)
@test convert(UFixed16, 1.1/typemax(UInt16))    == eps(UFixed16)
@test convert(UFixed{UInt32,16}, 1.1/typemax(UInt32)*2^16) == eps(UFixed{UInt32,16})
@test convert(UFixed{UInt64,3},  1.1/typemax(UInt64)*UInt64(2)^61)  == eps(UFixed{UInt64,3})
@test convert(UFixed{UInt128,7}, 1.1/typemax(UInt128)*UInt128(2)^121) == eps(UFixed{UInt128,7})

@test convert(UFixed8,  1.1f0/typemax(UInt8)) == eps(UFixed8)

@test convert(Float64, eps(UFixed8)) == 1/typemax(UInt8)
@test convert(Float32, eps(UFixed8)) == 1.0f0/typemax(UInt8)
@test convert(BigFloat, eps(UFixed8)) == BigFloat(1)/typemax(UInt8)
for T in (FixedPointNumbers.UF..., UF2...)
    @test convert(Bool, zero(T)) == false
    @test convert(Bool, one(T))  == true
    @test convert(Bool, convert(T, 0.2)) == true
    @test convert(Int, one(T)) == 1
    @test convert(Rational, one(T)) == 1
end
@test convert(Rational, convert(UFixed8, 0.5)) == 0x80//0xff
@test convert(UFixed16, one(UFixed8)) === one(UFixed16)
@test convert(UFixed16, UFixed8(0.5)).i === 0x8080
@test convert(UFixed{UInt16,7}, UFixed{UInt8,7}(0.504)) === UFixed{UInt16,7}(0.504)

@test  UFixed8(0.2) % UFixed8  === UFixed8(0.2)
@test UFixed14(1.2) % UFixed16 === UFixed16(0.20002)
@test UFixed14(1.2) % UFixed8  === UFixed8(0.196)

for i = 0.0:0.1:1.0
    @test i % UFixed8 === UFixed8(i)
end
@test ( 1.5 % UFixed8).i == round(Int,  1.5*255) % UInt8
@test (-0.3 % UFixed8).i == round(Int, -0.3*255) % UInt8

for i = 0.0:0.1:64.0
    @test i % UFixed10 === UFixed10(i)
end
@test (65.2 % UFixed10).i == round(Int, 65.2*1023) % UInt16
@test (-0.3 % UFixed10).i == round(Int, -0.3*1023) % UInt16

@test 1 % UFixed8 == 1
@test 2 % UFixed8 == UFixed8(0.996)

x = UFixed8(0b01010001, 0)
@test ~x == UFixed8(0b10101110, 0)
@test -x == 0xafuf8

@test isa(float(one(UFixed{UInt8,7})),   Float32)
@test isa(float(one(UFixed{UInt32,18})), Float64)
@test isa(float(one(UFixed{UInt32,25})), Float64)

for T in (FixedPointNumbers.UF..., UF2...)
    x = T(0x10,0)
    y = T(0x25,0)
    fx = float(x)
    fy = float(y)
    @test y > x
    @test y != x
    @test typeof(x+y) == T
    @test typeof((x+y)-y) == T
    @test typeof(x*y) == T
    @test typeof(x/y) == T
    @test_approx_eq(x+y, T(0x35,0))
    @test_approx_eq((x+y)-x, fy)
    @test_approx_eq((x-y)+y, fx)
    @test_approx_eq(x*y, convert(T, fx*fy))
    @test_approx_eq(x/y, convert(T, fx/fy))
    @test_approx_eq(x^2, convert(T, fx^2))
    @test_approx_eq(x^2.1f0, fx^2.1f0)
    @test_approx_eq(x^2.1, convert(Float64, x)^2.1)
end

function testtrunc{T}(inc::T)
    incf = convert(Float64, inc)
    tm = reinterpret(typemax(T))/reinterpret(one(T))
    x = zero(T)
    for i = 0 : min(1e6, reinterpret(typemax(T))-1)
        xf = incf*i
        try
            @test typeof(trunc(x)) == T
            @test trunc(x) == trunc(xf)
            @test typeof(round(x)) == T
            @test round(x) == round(xf)
            cxf = ceil(xf)
            if cxf < tm
                @test typeof(ceil(x)) == T
                @test ceil(x) == ceil(xf)
            end
            @test typeof(floor(x)) == T
            @test floor(x) == floor(xf)
            @test trunc(Int,x) == trunc(Int,xf)
            @test round(Int,x) == round(Int,xf)
            @test floor(Int,x) == floor(Int,xf)
            if cxf < tm
                @test ceil(Int,x) == ceil(Int,xf)
            end
        catch err
            println("Failed on x = ", x, ", xf = ", xf)
            rethrow(err)
        end
        x = convert(T, x+inc)
    end
end

for T in (FixedPointNumbers.UF..., UF2...)
    testtrunc(eps(T))
end

function testapprox{T}(::Type{T})
    for x = typemin(T):eps(T):typemax(T)-eps(T)
        y = x+eps(T)
        @test x ≈ y
        @test y ≈ x
        @test !(x ≈ y+eps(T))
    end
end
for T in FixedPointNumbers.UF
    testapprox(T)
end

@test !(UFixed8(0.5) < UFixed8(0.5))
@test UFixed8(0.5) <= UFixed8(0.5)

@test div(0x10uf8, 0x02uf8) == fld(0x10uf8, 0x02uf8) == 8
@test div(0x0fuf8, 0x02uf8) == fld(0x0fuf8, 0x02uf8) == 7
@test Base.fld1(0x10uf8, 0x02uf8) == 8
@test Base.fld1(0x0fuf8, 0x02uf8) == 8
@test mod(0x10uf8, 0x02uf8) == rem(0x10uf8, 0x02uf8) == 0
@test mod(0x0fuf8, 0x02uf8) == rem(0x0fuf8, 0x02uf8) == 0x01uf8
@test mod1(0x10uf8, 0x02uf8) == 0x02uf8
@test mod1(0x0fuf8, 0x02uf8) == 0x01uf8

r = 1uf8:1uf8:48uf8
@test length(r) == 48

counter = 0
for x in UFixed8(0):eps(UFixed8):UFixed8(1)
    counter += 1
end
@test counter == 256

# Promotion within UFixed
@test @inferred(promote(UFixed8(0.2), UFixed8(0.8))) ===
    (UFixed8(0.2), UFixed8(0.8))
@test @inferred(promote(UFixed{UInt16,3}(0.2), UFixed{UInt8,3}(0.86))) ===
    (UFixed{UInt16,3}(0.2), UFixed{UInt16,3}(0.86))
@test @inferred(promote(UFixed{UInt8,7}(0.197), UFixed{UInt8,4}(0.8))) ===
    (UFixed{UInt16,7}(0.197), UFixed{UInt16,7}(0.8))

@test UFixed{UInt16,16}(1)   == UFixed{UInt8,8}(1)
@test UFixed{UInt16,16}(0.2) == UFixed{UInt8,8}(0.2)
@test UFixed{UInt16,8}(1)   == UFixed{UInt8,8}(1)
@test UFixed{UInt16,8}(0.2) == UFixed{UInt8,8}(0.2)
@test UFixed{UInt16,16}(1)   == UFixed{UInt8,6}(1)
@test UFixed{UInt16,16}(0.20635) == UFixed{UInt8,6}(0.20635)
@test UFixed{UInt16,4}(1)   == UFixed{UInt8,6}(1)
@test UFixed{UInt16,4}(0.2) == UFixed{UInt8,6}(0.2)

@test promote_type(UFixed8,Float32,Int) == Float32
@test promote_type(UFixed8,Int,Float32) == Float32
@test promote_type(Int,UFixed8,Float32) == Float32
@test promote_type(Int,Float32,UFixed8) == Float32
@test promote_type(Float32,Int,UFixed8) == Float32
@test promote_type(Float32,UFixed8,Int) == Float32

# Show
x = 0xaauf8
iob = IOBuffer()
show(iob, x)
str = takebuf_string(iob)
@test str == "UFixed{UInt8,8}(0.667)"

# scaledual
function generic_scale!(C::AbstractArray, X::AbstractArray, s::Number)
    length(C) == length(X) || error("C must be the same length as X")
    for i = 1:length(X)
        @inbounds C[i] = X[i]*s
    end
    C
end

a = rand(UInt8, 10)
rfloat = similar(a, Float32)
rfixed = similar(rfloat)
af8 = reinterpret(UFixed8, a)

b = 0.5
bd, eld = scaledual(b, af8[1])
@assert b*a[1] == bd*eld

b, ad = scaledual(0.5, a)
@test b == 0.5
@test ad == a
b, ad = scaledual(0.5, ad)
generic_scale!(rfloat, a, 0.5)
generic_scale!(rfixed, ad, b)
@test rfloat == rfixed

# reductions
a = UFixed8[0xffuf8, 0xffuf8]
@test sum(a) == 2.0
@test sum(a, 1) == [2.0]

a = UFixed14[3.2, 2.4]
acmp = Float64(a[1])*Float64(a[2])
@test prod(a) == acmp
@test prod(a, 1) == [acmp]

x = UFixed8(0.3)
for T in (Float16, Float32, Float64, BigFloat)
    y = convert(T, x)
    @test isa(y, T)
end

for T in (UFixed{UInt8,8}, UFixed{UInt8,6},
          UFixed{UInt16,16}, UFixed{UInt16,14},
          UFixed{UInt32,32}, UFixed{UInt32,30},
          UFixed{UInt64,64}, UFixed{UInt64,62})
    a = rand(T)
    @test isa(a, T)
    a = rand(T, (3, 5))
    @test isa(a, Array{T,2})
    @test size(a) == (3,5)
end


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

name = "igraphBuilder"
version = v"0.0.1"

# Collection of sources required to build igraphBuilder
sources = [
    "https://github.com/igraph/igraph/releases/download/0.7.1/igraph-0.7.1.tar.gz" =>
    "d978030e27369bf698f3816ab70aa9141e9baf81c56cc4f55efbe5489b46b0df",

]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir
cd igraph-0.7.1/
./configure --prefix=$prefix --host=$target
make
make install

"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    Linux(:i686, libc=:glibc),
    Linux(:x86_64, libc=:glibc),
#    Linux(:aarch64, libc=:glibc),
    Linux(:armv7l, libc=:glibc, call_abi=:eabihf),
    Linux(:powerpc64le, libc=:glibc),
    MacOS(:x86_64)
]

# The products that we will ensure are always built
products(prefix) = [
    LibraryProduct(prefix, "libigraph", :libigraph)
]

# Dependencies that must be installed before this package can be built
dependencies = [
    "https://github.com/bicycle1885/ZlibBuilder/releases/download/v1.0.4/build_Zlib.v1.2.11.jl",
    "https://github.com/bicycle1885/XML2Builder/releases/download/v1.0.2/build_XML2Builder.v2.9.9.jl"
]

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

max_stars_repo_name: jstrube/igraphBuilder

max_issues_repo_name: jstrube/igraphBuilder

max_forks_repo_name: jstrube/igraphBuilder