text:
using Lazy: @>, takewhile
import Base: run
export Experiment, setup, run, addcolumn

const default_moment_resolution = 1.5ms
const default_input_resolution = (1/60)s
const exp_width = 1024
const exp_height = 768

const experiment_context = Array{Nullable{Experiment}}()
experiment_context[] = Nullable()

# internal function, usd to find the current experiment
function get_experiment()
  if isnull(experiment_context[])
    error("Unknown experiment context, call me inside `setup` or during an"*
          " experiment.")
  else
    get(experiment_context[])
  end
end

# internal function used to determine if there is an experiment running
function in_experiment()
  !isnull(experiment_context[])
end

function experiment_running()
  !isnull(experiment_context[]) && flags(get(experiment_context[])).processing
end

# internal functions used to update and retrieve the stack trace
# where the currently running moment was defined (improving error message
# readability)
const current_moment_trace = Array{Vector{StackFrame}}()
update_trace(m::AbstractMoment) =
  !isempty(moment_trace(m)) ? current_moment_trace[] = moment_trace(m) : nothing
update_trace(m::MomentSequence) =
  current_moment_trace[] = moment_trace(m.data[1])
moment_trace() = current_moment_trace[]
function moment_trace_string()
  if in_experiment()
    "\nOn trial $(Weber.trial()), offset $(Weber.offset())"*
    reduce(*,"",map(x -> string(x)*"\n",moment_trace()))
  else
    ""
  end
end

"""
    Weber.trial()

Returns the current trial of the experiment.
"""
trial(exp) = data(exp).trial
trial() = trial(get_experiment())

"""
    Weber.offset()

Returns the current offset. The offset represents a well defined time in the
experiment. The offset is typically incremented once for every call to
[`addpractice`](@ref) [`addtrial`](@ref) and [`addbreak`](@ref), unless you use
[`@addtrials`](@ref). You can use the offset to restart the experiment from a
well defined location.

!!! warning

    For offsets to be well defined, all calls to [`moment`](@ref) and
    [`@addtrials`](@ref) must follow the [guidlines](advanced.md) in the user
    guide. In particular, moments should not rely on state that changes during
    the experiment unless they are wrapped in an @addtrials macro.

"""
offset(exp) = data(exp).offset
offset() = offset(get_experiment())

"""
    Weber.tick()

With microsecond precision, this returns the number of elapsed seconds from the
start of the experiment to the start of the most recent moment.

If there is no experiment running, this returns the time since epoch with
microsecond precision.
"""
tick(exp) = data(exp).last_time
function tick()
  if isnull(experiment_context[])
    precise_time()
  else
    tick(get(experiment_context[]))
  end
end

"""
   Weber.metadata() = Dict{Symbol,Any}()

Returns metadata for this experiment. You can store
global state, specific to this experiment, in this dictionary.
"""
metadata(exp) = info(exp).meta
metadata() = metadata(get_experiment())

"""
    addcolumn(column::Symbol)

Adds a column to be recorded in the data file.

This function must be called during setup.  It cannot be called once the
experiment has begun. Repeatedly adding the same column only adds the column
once. After adding a column you can include that column as a keyword argument
to [`record`](@ref). You need not write to the column for every record.
If left out, the column will be empty in the resulting row of the data file.
"""
function addcolumn{T <: BaseExperiment}(exp::T,col::Symbol)
  if flags(exp).processing
    error("You cannot change the data file header once the experiment starts! ",
          "Make sure you call `addcolumn` during setup-time, not run-time.")
  end
  if col ∉ info(exp).header
    push!(info(exp).header,col)
  end
end
function addcolumn(exp::ExtendedExperiment,col::Symbol)
  addcolumn(next(exp),col)
end
addcolumn(col::Symbol) = addcolumn(get_experiment(),col)

"""
    Experiment([skip=0],[columns=[symbols...]],[debug=false],
               [moment_resolution=0.0015],[data_dir="data"],
               [width=1024],[height=768],[warn_on_trials_only=true],[extensions=[]])

Prepares a new experiment to be run.

# Keyword Arguments
* **skip** the number of offsets to skip. Allows restarting of an experiment
  somewhere in the middle. When an experiment is terminated, the most
  recent offset is reported. The offset is also recorded in each row
  of the resulting data file (also reported on exit).
* **columns** the names (as symbols) of columns that will be recorded during the
  experiment (using `record`). These can be set to fixed values (using :name =>
  value), or be filled in during a call to record (:name). The column `:value`
  is always included here, even if not specified, since there are number of
  events recorded automatically which make use of this column.
* **debug** if true, experiment will show in a windowed view
* **moment_resolution** the desired precision that moments
  should be presented at. Warnings will be printed for moments that
  lack this precision.
* **data_dir** the directory where data files should be stored (can be set to
  nothing to prevent a file from being created)
* **width** and **height** specified the screen resolution during the experiment
* **extensions** an array of Weber.Extension objects, which [extend](extend.md)
  the behavior of an experiment.
* **warn_on_trials_only** when true, latency warnings are only displayed
  when the trial count is greater than 0. Thus, practice and breaks
  that occur before the first trial do not raise latency warnings.

"""
function Experiment(;skip=0,columns=Symbol[],debug=false,
                    moment_resolution = default_moment_resolution,
                    data_dir = "data",
                    null_window = false,
                    hide_output = false,
                    input_resolution = default_input_resolution,
                    extensions = Extension[],
                    width=exp_width,height=exp_height,
                    warn_on_trials_only = true)
  if !sound_is_setup() && !null_window
    setup_sound()
    clear_sound_cache()
    empty!(_image_cache)
    empty!(convert_cache)
  end
  TimedSound.sound_setup_state.hooks = WeberSoundHooks()
  TimedSound.sound_setup_state.cache = true

  if !(data_dir == nothing || hide_output)
    mkpath(data_dir)
  elseif !hide_output
    warn("No directory specified for saving data. ALL DATA FROM THIS",
                  " EXPERIMENT WILL BE LOST!!! Refer to the documentation for",
                  "`Experiment`.")
  end

  moment_resolution_s = ustrip(inseconds(moment_resolution))
  input_resolution_s = ustrip(inseconds(input_resolution))
  if moment_resolution_s < approx_timer_resolution
    warn("The desired timing resolution of $moment_resolution "*
                  "seconds is probably not achievable on your system. The "*
                  "approximate minimum is $approx_timer_resolution seconds."*
                  " Try changing the moment_resolution to a higher value (see "*
                  "documentation for `Experiment`).")
  end

  meta = Dict{Symbol,Any}()
  start_time = precise_time()
  start_date = now()
  timestr = Dates.format(start_date,"yyyy-mm-dd__HH_MM_SS")
  info_values = filter(x -> x isa Pair,columns)
  reserved_columns = filter(x -> !(x isa Pair),columns)
  info_str = join(map(x -> x[2],info_values),"_")
  filename = (data_dir == nothing || hide_output ? Nullable() :
              Nullable(joinpath(data_dir,info_str*"_"*timestr*".csv")))
  einfo = ExperimentInfo(info_values,meta,input_resolution_s,
                         moment_resolution_s,start_date,reserved_columns,
                         filename,hide_output,warn_on_trials_only)

  offset = 0
  trial = 0
  trial_watcher = (e) -> nothing
  last_time = 0.0
  next_moment = 0.0
  pause_mode = Running
  moments = [MomentQueue()]
  streamers = Dict{Int,TimedSound.Streamer}()
  last_good_delta = -1.0
  last_bad_delta = -1.0
  data = ExperimentData(offset,trial,skip,last_time,next_moment,trial_watcher,
                        pause_mode,moments,streamers,last_good_delta,
                        last_bad_delta)

  running = processing = false
  flags = ExperimentFlags(running,processing)

  win = window(width,height,fullscreen=!debug,accel=!debug,null=null_window)

  final_exp = extend(UnextendedExperiment(einfo,data,flags,win),extensions)
  addcolumn(final_exp,:value)
  final_exp
end

"""
    setup(fn,experiment)

Setup the experiment, adding breaks, practice, and trials.

Setup creates the context necessary to generate elements of an experiment. All
calls to `addtrial`, `addbreak` and `addpractice` must be called inside of
`fn`. This function must be called before `run`.
"""
setup(fn::Function,exp::ExtendedExperiment) = setup(fn,next(exp))
function setup{T <: BaseExperiment}(fn::Function,exp::T)
  try
    # setup all trial moments for this experiment
    experiment_context[] = Nullable(top(exp))
    addmoment(top(exp),moment())
    fn()
    experiment_context[] = Nullable{Experiment}()
  catch e
    close(win(exp))
    # gc_enable(true)
    rethrow(e)
  end
  nothing
end

warmup_run(exp::Experiment{NullWindow}) = nothing
function warmup_run(exp::Experiment)
  # warm up JIT compilation
  warm_up = Experiment(null_window=true,hide_output=true)
  setup(warm_up) do
    addtrial(moment(record,"a"),
             moment(record,"b") >>
             moment(record,"d"),
             moment(record,"c"))
    @addtrials let i = 0
      @addtrials while i < 3
        addtrial(moment(() -> i+=1),
                 moment(record,"a"),
                 moment(record,"b"),
                 moment(record,"c"))
      end
    end
  end
  run(warm_up,await_input=false)
end

function pause(exp,message,time,firstpause=true)
  flags(exp).running = false
  record(top(exp),"paused")
  if firstpause
    save_display(win(exp))
    pause_sounds()
  end
  overlay = visual(colorant"gray",priority=Inf) + visual(message,priority=Inf)
  display(win(exp),overlay)
end

function unpause(exp,time)
  record(top(exp),"unpaused")
  data(exp).pause_mode = Running

  restore_display(win(exp))
  resume_sounds()

  data(exp).last_bad_delta = -1.0
  data(exp).last_good_delta = -1.0
  flags(exp).running = true
  process_event(exp,EndPauseEvent(time))
end

const Running = 0
const ToExit = 1
const Unfocused = 2
const Error = 3

function watch_pauses(exp,e)
  if data(exp).pause_mode == Running && iskeydown(e,key":escape:")
    pause(exp,"Exit? [Y for yes, or N for no]",time(e))
    data(exp).pause_mode = ToExit
  elseif data(exp).pause_mode == Running && isunfocused(e) && flags(exp).processing
    pause(exp,"Waiting for window focus...",time(e))
    data(exp).pause_mode = Unfocused
  elseif data(exp).pause_mode == ToExit && iskeydown(e,key"y")
    record(top(exp),"terminated")
    endexperiment(exp)
  elseif data(exp).pause_mode == ToExit && iskeydown(e,key"n")
    unpause(exp,time(e))
  elseif data(exp).pause_mode == Unfocused && isfocused(e)
    if flags(exp).processing
      pause(exp,"Paused. [To exit hit Y, to resume hit N]",time(e),false)
      data(exp).pause_mode = ToExit
    else
      data(exp).pause_mode = Running
      flags(exp).running = true
    end
  end
end

process_event(exp::Experiment,event::QuitEvent) = endexperiment(exp)
function process_event(exp::Experiment,event)
  if flags(exp).running
    data(exp).trial_watcher(event)
    process(exp,data(exp).moments,event)
  end
  if flags(exp).processing
    watch_pauses(exp,event)
  end
end

const sleep_resolution = 0.05
const sleep_amount = 0.002

const gc_time = 1

"""
    run(experiment;await_input=true)

Runs an experiment. You must call `setup` first.

By default, on windows, this function waits for user input before returning.
This prevents a console from closing at the end of an experiment, preventing the
user from viewing important messages. The exception is if run is called form
within Juno: await_input should never be set to true in this case.
"""
run(exp::ExtendedExperiment;keys...) = run(next(exp);keys...)
function run{T <: BaseExperiment}(
  exp::T;await_input=!Juno.isactive() && !is_apple())

  if Juno.isactive()
    if await_input
      error("`await_input` must be false when Juno is active.")
    else
      warn("""
        Running Weber experiment directly in Juno. Consider using @read_args.
        When runing directly in Juno, the experiment may not respond to input
        correctly.  To solve this problem on Windows, once the experiment
        begins, hit alt-tab to switch away from the expeirment, and then hit
        alt-tab again to switch back to the experiment. On Mac you can use
        command-tab to do the same thing.
      """)
    end
  end

  warmup_run(exp)
  # println("========================================")
  # println("Completed warmup run.")
  try
    record_header(exp)
    focus(win(exp))

    experiment_context[] = Nullable(top(exp))
    data(exp).pause_mode = Running
    flags(exp).processing = true
    flags(exp).running = true

    start = precise_time()
    tick = data(exp).last_time = last_input = last_delta = 0.0
    prepare!(data(exp).moments[1],Inf)
    while flags(exp).processing && !isempty(data(exp).moments)
      tick = data(exp).last_time = precise_time() - start

      # notify all moments about the new time
      if flags(exp).running
        process(exp,data(exp).moments,tick)
      end

      # handle all events (handles pauses, and notifys moments of events)
      if tick - last_input > info(exp).input_resolution
        poll_events(process_event,top(exp),tick)
        last_input = tick
      end

      # handle auditory streams
      next_stream = Inf
      for streamer in values(data(exp).streamers)
        if tick > streamer.next_stream
          process(streamer)
          next_stream = min(next_stream,streamer.next_stream)
        end
      end

      # # report on any irregularity in the timing of moments
      # if flags(exp).running && tick - last_delta > last_delta_resolution
      #   report_deltas(exp)
      #   last_delta = tick
      # end

      # refresh screen
      refresh_display(win(exp))

      # if after all this processing there's still plenty of time left
      # then sleep for a little while. (pausing also sleeps the loop)
      new_tick = precise_time() - start
      stream_len = ustrip(TimedSound.sound_setup_state.stream_unit/samplerate())
      if !flags(exp).running
        gc()
        sleep(sleep_amount)
      elseif ((new_tick - last_delta) > sleep_resolution &&
              (new_tick - last_input) > sleep_resolution &&
              new_tick + 0.2stream_len < next_stream &&
              (data(exp).next_moment - new_tick) > sleep_resolution)
        if (data(exp).next_moment - new_tick) > gc_time
          gc()
        end
        sleep(min(sleep_amount,0.05stream_len))
      end
    end
  catch e
    if !isempty(moment_trace())
      print(STDERR,"\nException during a moment defined")
      foreach(x -> println(STDERR,x),moment_trace())
    end
    if await_input
      showerror(STDERR,e,catch_backtrace())
    else
      rethrow(e)
    end
  finally
    record(top(exp),"closed")
    flags(exp).running = false
    flags(exp).processing = false
    experiment_context[] = Nullable()
    close(win(exp))
    # gc_enable(true)
    if !info(exp).hide_output
      info("Experiment terminated at offset $(data(exp).offset).")
      if !isnull(info(exp).file)
        info("Data recorded to: $(abspath(get(info(exp).file)))")
      end
    end
  end
  if await_input
    println("Hit enter to end experiment.")
    readline(STDIN)
  end
  nothing
end

function endexperiment(e::Experiment)
  flags(e).running = false
  flags(e).processing = false
end

function process(exp::Experiment,queues::Array{MomentQueue},x)
  filter!(queues) do queue
    !isempty(process(exp,queue,x))
  end
end

function skip_offsets(exp,queue)
  while !isempty(queue) && is_moment_skipped(exp,front(queue))
    dequeue!(queue)
  end
end

function prepare!(queue::MomentQueue,t::Float64)
  if required_delta_t(front(queue)) > 0.0 || very_first(queue)
    delta_t = required_delta_t(front(queue))
    prepare!(front(queue),t + delta_t)

    @_ queue begin
      drop(_,1)
      takewhile(m -> required_delta_t(m) == 0.0,_)
      foreach(m -> prepare!(m,t + delta_t),_)
    end
  end
end

function process(exp::Experiment,queue::MomentQueue,event::ExpEvent)
  skip_offsets(exp,queue)

  if !isempty(queue)
    moment = front(queue)
    update_trace(moment)
    handled = handle(exp,queue,moment,event)
    if handled
      prepare!(queue,time(event))
      data(exp).next_moment = minimum(map(next_moment_time,data(exp).moments))
    end
  end

  queue
end

roundstr(x,n=6) = x > 10.0^-n ? string(round(x,n)) : "≤1e-$n"
function process(exp::Experiment,queue::MomentQueue,t::Float64)
  skip_offsets(exp,queue)

  if !isempty(queue)
    start_time = precise_time()
    moment = front(queue)
    event_time = delta_t(moment) + queue.last
    if event_time - t <= info(exp).moment_resolution
      offset = t - start_time
      run_time = offset + precise_time()
      while event_time > run_time
        run_time = offset + precise_time()
      end
      data(exp).last_time = run_time
      update_trace(moment)
      if handle(exp,queue,moment,run_time)
        prepare!(queue,run_time)

        latency = run_time - event_time

        if (latency > info(exp).moment_resolution &&
            warn_delta_t(moment) &&
            show_sound_latency_warnings() &&
            !info(exp).hide_output)
          warn(
            "Delivered a moment with a high latency ($(roundstr(latency))
             seconds). This often happens at the start of an experiment, but
             should rarely, if ever, occur throughout the experiment. To reduce
             latency, reduce the amount of slow code in moments, close programs,
             or run on a faster machine. Or, if this amount of latency is
             acceptable, you should increase `moment_resolution` when you call
             `Experiment`.\nMoment: $moment \n\n"*moment_trace_string())
          record("high_latency",value=latency)
        end

        data(exp).next_moment = minimum(map(next_moment_time,data(exp).moments))
      end
    end
  end

  queue
end


text:
# this is for keeping the basic functionalities

using BitIntegers: @define_integers

import Base: unsafe_getindex, ==, show, promote_rule
using Base: @_inline_meta, @propagate_inbounds, @_propagate_inbounds_meta
import Base.GC: @preserve

struct ShortString{T} <: AbstractString where {T}
    size_content::T
end

"""Check if a string of size `sz` can be stored in ShortString{T}"""
function check_size(T, sz)
    max_len = sizeof(T) - size_bytes(T)  # the last few bytes are used to store the length
    if sz > max_len
        throw(ErrorException("sizeof(::$T) must be shorter than or equal to $(max_len) in length; you have supplied a string of size $sz"))
    end
end

"""Calculate the number of bytes required to store the size of the ShortString"""
size_bytes(::Type{T}) where {T} = (count_ones(sizeof(T)-1)+7)>>3

"""Calculate a mask to get the size stored in the ShortString"""
size_mask(T) = T((1<<(size_bytes(T)*8)) - 1)
size_mask(s::ShortString{T}) where {T} = size_mask(T)

"""The size of the chunk used to process String values"""
const CHUNKSZ   = sizeof(UInt)

"""The number of bits in the chunk type used to process String values"""
const CHUNKBITS = sizeof(UInt) == 4 ? 32 : 64

"""Internal function to pick up a byte at the given index in a ShortString"""
@inline _get_byte(s::ShortString, i::Int) = (s.size_content >>> (8*(sizeof(s) - i)))%UInt8

"""
Internal function to pick up a UInt32 (i.e. to contain 1 Char) at the given index
in a ShortString
"""
@inline function _get_word(s::ShortString{T}, i::Int) where {T}
    sz = sizeof(T)
    if sz <= 4
        # Shift up by 0-3 bytes
        (s.size_content%UInt32) << (8*(i + 3 - sz))
    else
        (s.size_content >>> (8*(sz - i - 3)))%UInt32
    end
end

"""Internal function to get the UInt32 representation of a Char from an index in a ShortString"""
@inline function _get_char(str::ShortString, pos::Int)
    chr = _get_word(str, pos)
    typ = chr >>> 28
    chr & ~ifelse(typ < 0x8, 0xffffff,
                  ifelse(typ < 0xe, 0x00ffff,
                         ifelse(typ < 0xf, 0x0000ff, 0x000000)))
end

"""Internal function, given a String and it's size in bytes, load it into a value of type T"""
@inline function _ss(::Type{T}, str::String, sz) where {T}
    if sizeof(T) <= sizeof(UInt)
        unsafe_load(reinterpret(Ptr{T}, pointer(str)))
    else
        pnt = reinterpret(Ptr{UInt}, pointer(str))
        fin = pnt + sz
        val = T(unsafe_load(pnt))
        off = CHUNKBITS
        while (pnt += CHUNKSZ) <= fin
            val |= T(unsafe_load(pnt)) << off
            off += CHUNKBITS
        end
        val
    end
end

# This can be optimized later, right now, it's just working 1 byte at a time
# to avoid accessing past the end
@inline function _ss(::Type{T}, str::SubString{String}, sz) where {T}
    pnt = pointer(str)
    fin = pnt + sz
    val = T(unsafe_load(pnt))
    off = 8
    while (pnt += 1) <= fin
        val |= T(unsafe_load(pnt)) << off
        off += 8
    end
    val
end

function ShortString{T}(s::Union{String,SubString{String}}) where {T}
    sz = sizeof(s)
    sz === 0 && return ShortString{T}(T(0))
    check_size(T, sz)
    bw = 8(sizeof(T) - sz)
    @preserve s ShortString{T}((ntoh(_ss(T, s, sz)) >>> bw) << bw | T(sz))
end

ShortString{T}(s::ShortString{T}) where {T} = s

function ShortString{T}(s::ShortString{S}) where {T, S}
    sz = sizeof(s)
    check_size(T, sz)
    # Flip it so empty bytes are at start, grow/shrink it, flip it back
    # size_mask(S) will return a mask for getting the size for Shorting Strings in (content size)
    # format, so something like 00001111 in binary.
    #  ~size_mask(S) will yield 11110000 which can be used as a mask to extract the content
    content = ntoh(T(ntoh(s.size_content & ~size_mask(S))))
    ShortString{T}(content | T(sz))
end

"""Amount to shift ShortString value by for each UInt sized chunk"""
const SHFT_INT = UInt === UInt32 ? 2 : 3

# Optimized conversion of a ShortString to a String
function String(s::ShortString{T}) where {T}
    len = sizeof(s)
    len === 0 && return ""
    val = ntoh(s.size_content)
    sv = Base.StringVector(len)
    # Loop over UInt64 sized chunks here
    pnt = reinterpret(Ptr{UInt}, pointer(sv))
    for i = 1:(len + sizeof(UInt) - 1) >>> SHFT_INT
        unsafe_store!(pnt, val % UInt)
        val >>>= 8*sizeof(UInt)
        pnt += sizeof(UInt)
    end
    String(sv)
end

Base.codeunit(s::ShortString) = UInt8
Base.codeunit(s::ShortString, i) = codeunits(String(s), i)
Base.codeunit(s::ShortString, i::Integer) = codeunit(String(s), i)
Base.codeunits(s::ShortString) = codeunits(String(s))

Base.convert(::ShortString{T}, s::String) where {T} = ShortString{T}(s)
Base.convert(::String, ss::ShortString) = String(ss)

Base.sizeof(s::ShortString) = Int(s.size_content & size_mask(s))

Base.firstindex(::ShortString) = 1
Base.lastindex(s::ShortString) = sizeof(s)
Base.ncodeunits(s::ShortString) = sizeof(s)

# Checks top two bits of first byte of character to see if valid position
isvalid(s::String, i::Integer) = (0 < i <= sizeof(s)) && ((_get_byte(s, i) & 0xc0) != 0x80)

@inline function Base.iterate(s::ShortString, i::Int=1)
    0 < i <= ncodeunits(s) || return nothing
    chr = _get_word(s, i)
    typ = chr >>> 24
    if typ < 0x80
        (reinterpret(Char, chr & 0xFF00_0000), i + 1)
    elseif typ < 0xe0
        (reinterpret(Char, chr & 0xFFFF_0000), i + 2)
    elseif typ < 0xf0
        (reinterpret(Char, chr & 0xFFFF_FF00), i + 3)
    else
        (reinterpret(Char, chr), i + 4)
    end
end

@inline function Base.isascii(s::ShortString{T}) where {T}
    val = s.size_content >>> (8*size_bytes(T))
    for i in 1:(sizeof(T)-size_bytes(T))
        iszero(val & 0x80) || return false
        val >>>= 8
    end
    return true
end

function Base.length(s::ShortString{T}) where T
    isascii(s) && return ncodeunits(s)

    # else have to do it the hard way:
    i = 0
    len = 0
    while i < ncodeunits(s)
        shifted = s.size_content >>> (8*(sizeof(T) - i))
        i += if shifted % UInt8 <= 0x7f  # 1 byte character
            1
        elseif shifted % UInt16 <= 0x7ff  # 2 byte character
            2
        elseif shifted % UInt32 <= 0xffff  # 3 byte character
            3
        else  # 4 byte character
            4
        end
        len += 1
    end
    return len
end

@propagate_inbounds function Base.getindex(str::ShortString, pos::Int=1)
    @_inline_meta()
    @boundscheck checkbounds(str, pos)
    reinterpret(Char, _get_char(str, pos))
end

function ==(s::ShortString{S}, b::Union{String, SubString{String}}) where {S}
    ncodeunits(b) == ncodeunits(s) || return false
    return s == ShortString{S}(b)
end
function ==(s::ShortString, b::AbstractString)
    # Could be a string type that might not use UTF8 encoding and that we don't have a
    # constructor for. Defer to equality that type probably has defined on `String`
    return String(s) == b
end

==(a::AbstractString, b::ShortString) = b == a
function ==(a::ShortString{S}, b::ShortString{S}) where {S}
    return a.size_content == b.size_content
end
function ==(a::ShortString{A}, b::ShortString{B}) where {A,B}
    ncodeunits(a) == ncodeunits(b) || return false
    # compare if equal after dropping size bits and
    # flipping so that the empty bytes are at the start
    ntoh(a.size_content & ~size_mask(A)) == ntoh(b.size_content & ~size_mask(B))
end

function Base.cmp(a::ShortString{S}, b::ShortString{S}) where {S}
    return cmp(a.size_content, b.size_content)
end

promote_rule(::Type{String}, ::Type{ShortString{S}}) where {S} = String

function promote_rule(::Type{ShortString{T}}, ::Type{ShortString{S}}) where {T,S}
    if sizeof(T) >= sizeof(S)
        return ShortString{promote_rule(T,S)}
    else
        return ShortString{promote_rule(S,T)}
    end
end

size_content(s::ShortString) = s.size_content

@define_integers 2048 Int2048 UInt2048

"""These are the default types used to for selecting the size of a ShortString"""
const def_types = (UInt32, UInt64, UInt128, UInt256, UInt512, UInt1024, UInt2048)

for T in def_types
    max_len = sizeof(T) - size_bytes(T)
    constructor_name = Symbol(:ShortString, max_len)
    macro_name = Symbol(:ss, max_len, :_str)

    @eval const $constructor_name = ShortString{$T}
    @eval macro $(macro_name)(s)
        Expr(:call, $constructor_name, s)
    end
end

"""
Return a ShortString type that can hold maxlen codeunits
The keyword parameter `types` can be used to pass a list of types
which can be used to store the string
If no type is large enough, then an `ArgumentError` is thrown
"""
function get_type(maxlen; types=def_types)
    for T in types
        maxlen <= sizeof(T) - size_bytes(T) && return ShortString{T}
    end
    throw(ArgumentError("$maxlen is too large to fit into any of the provided types: $types"))
end

ShortString(str::Union{String,SubString{String}}, maxlen = sizeof(str); types=def_types) =
    get_type(maxlen, types=types)(str)

macro ss_str(str, max="0")
    :( ShortString($str, $(parse(Int, max))) )
end

fsort(v::Vector{ShortString{T}}; rev = false) where {T} =
    sort(v, rev = rev, by = size_content, alg = RadixSort)
fsort!(v::Vector{ShortString{T}}; rev = false) where {T} =
    sort!(v, rev = rev, by = size_content, alg = RadixSort)

fsortperm(v::Vector{ShortString{T}}; rev = false) where {T} = sortperm(v, rev = rev)

ext: jl


text:
# This file is part of FlagSets.jl project, which is licensed under BDS-3-Clause.
# See file LICENSE.md for more information.

@enum MyFlagType Foo Bar Baz Quux

@testset "Basic operation" begin
    # Inline definition
    @flagset TFlag1 RoundDown RoundUp RoundNearest
    @test Int(TFlag1([RoundDown])) == 1
    @test Int(TFlag1([RoundUp])) == 2
    @test Int(TFlag1([RoundNearest])) == 4
    @test TFlag1(1) == TFlag1([RoundDown])
    @test TFlag1(2) == TFlag1([RoundUp])
    @test TFlag1(4) == TFlag1([RoundNearest])
    @test flags(TFlag1) == (RoundDown, RoundUp, RoundNearest)
    @test length(TFlag1([RoundDown, RoundUp])) == 2
    @test !iszero(TFlag1([RoundUp, RoundNearest]))
    @test iszero(TFlag1())
    @test eltype(TFlag1) == RoundingMode

    # Block definition
    @flagset TFlag2 begin
        down = RoundDown
        up = RoundUp
        near = RoundNearest
    end
    @test Int(TFlag2(down = true)) == 1
    @test Int(TFlag2(up = true)) == 2
    @test Int(TFlag2(near = true)) == 4
    @test TFlag2(down = true, up = false, near = true) == TFlag2([RoundDown, RoundNearest])
    @test length(TFlag2([RoundNearest])) == 1
    @test !iszero(TFlag2([RoundNearest]))
    @test iszero(TFlag2())
    @test eltype(TFlag2) == RoundingMode

    # Explicit numbering, inline
    @flagset TFlag3 {_,_} down = 2 --> RoundDown RoundUp near = 16 --> RoundNearest
    @test Int(TFlag3(down = true)) == 2
    @test Int(TFlag3([RoundUp])) == 4
    @test Int(TFlag3(near = true)) == 16
    @test isempty(typemin(TFlag3))
    @test Int(typemax(TFlag3)) == 22
    @test length(typemax(TFlag3)) == 3
    @test TFlag3(down = false, near = true) == TFlag3([RoundNearest])
    @test flags(TFlag3) == (RoundDown, RoundUp, RoundNearest)
    @test eltype(TFlag3) == RoundingMode

    # Construct from enums
    @flagset TFlag4 {} begin
        foo = 2 --> Foo
        bar = Bar
        16 --> Baz
        quux = 64 --> Quux
    end
    @test Int(TFlag4(foo = true)) == 2
    @test Int(TFlag4([Bar])) == 4
    @test Int(TFlag4([Baz])) == 16
    @test isempty(typemin(TFlag4))
    @test Int(typemax(TFlag4)) == 86
    @test length(typemax(TFlag4)) == 4
    @test TFlag4(foo = false, quux = true) == TFlag4([Quux])
    @test flags(TFlag4) == (Foo, Bar, Baz, Quux)
    @test eltype(TFlag4) == MyFlagType

    # UInt128
    @flagset TFlag4a a = UInt128(1) << 96 --> Foo b = UInt128(1) << 120 --> Bar
    @test length([TFlag4a([Foo, Bar])...]) == 2
    @test flags(TFlag4a([Foo, Bar])) === (Foo, Bar)
    @test all(s1 === s2 for (s1, s2) ∈ zip(TFlag4a([Foo, Bar]), [Foo, Bar]))

    # BigInt
    @flagset TFlag4b a = big(1) << 160 --> Int b = big(1) << 192 --> Inf
    @test length([TFlag4b([Int, Inf])...]) == 2
    @test flags(TFlag4b([Int, Inf])) === (Int, Inf)
    @test all(s1 === s2 for (s1, s2) ∈ zip(TFlag4b([Int, Inf]), [Int, Inf]))
end

@testset "Mask and set operations" begin
    # Mask operations
    @test TFlag1([RoundDown]) | TFlag1([RoundUp]) | TFlag1([RoundNearest]) ==
          TFlag1([RoundDown, RoundUp, RoundNearest])
    @test Int(TFlag1(7)) == 7
    @test TFlag1(7) == TFlag1([RoundDown]) | TFlag1([RoundUp]) | TFlag1([RoundNearest])
    @test TFlag1(7) & TFlag1([RoundDown]) == TFlag1([RoundDown])
    @test TFlag1(7) ⊻ TFlag1([RoundDown]) == TFlag1([RoundUp, RoundNearest])
    @test ~TFlag1([RoundDown]) == TFlag1([RoundUp, RoundNearest])
    @test Int(TFlag1([RoundDown])) < Int(TFlag1([RoundUp])) < Int(TFlag1([RoundNearest]))
    @test Int(TFlag1([RoundDown]) | TFlag1([RoundUp])) < Int(TFlag1([RoundNearest]))

    # Set operations
    @test TFlag1([RoundDown]) ∪ TFlag1([RoundUp]) ∪ TFlag1([RoundNearest]) ==
          TFlag1([RoundDown, RoundUp, RoundNearest])
    @test TFlag1(7) == TFlag1([RoundDown]) ∪ TFlag1([RoundUp]) ∪ TFlag1([RoundNearest])
    @test TFlag1(7) ∩ TFlag1([RoundDown]) == TFlag1([RoundDown])
    @test TFlag1(7) ⊻ TFlag1([RoundDown]) == TFlag1([RoundUp, RoundNearest])
    @test setdiff(typemax(TFlag1), TFlag1([RoundDown])) == TFlag1([RoundUp, RoundNearest])
    @test RoundDown ∈ TFlag1([RoundDown])
    @test RoundDown ∉ TFlag1([RoundUp, RoundNearest])
    @test TFlag1([RoundUp]) ⊊ TFlag1([RoundUp, RoundNearest])
    @test TFlag1([RoundNearest]) ⊊ TFlag1([RoundUp, RoundNearest])
    @test !(TFlag1([RoundUp, RoundNearest]) ⊊ TFlag1([RoundUp, RoundNearest]))
    @test !(TFlag1([RoundDown, RoundUp]) ⊊ TFlag1([RoundUp, RoundNearest]))
end # testset

@testset "Type properties" begin
    # Default integer typing and compatibility with @symbol_flagset
    @flagset TFlag5 :flag5a :flag5b
    @test eltype(TFlag5) == Symbol
    @test TFlag5(flag5a = true) == TFlag5(:flag5a)
    @test typeof(Integer(TFlag5(:flag5a))) == UInt32
    @test typeof(TFlag5) == DataType
    @test typeof(TFlag5(:flag5a)) <: TFlag5 <: FlagSet <: AbstractSet
    @test isbitstype(Flag5)
    @test isbits(TFlag5(:flag5a))

    # Construct non-default and infer type
    @flagset TFlag6 {_, UInt16} begin
        foo = 2 --> Foo
        bar = Bar
        16 --> Baz
        quux = 64 --> Quux
    end
    @test typeof(Integer(TFlag6(foo = true))) == UInt16
    @test UInt8(TFlag6(foo = true)) === 0x02
    @test UInt16(TFlag6([Bar])) === 0x0004
    @test UInt128(TFlag6([Baz])) === 0x00000000000000000000000000000010
    @test typeof(TFlag6()) <: TFlag6 <: FlagSet <: AbstractSet

    # Explicit bits of non-default types
    @flagset TFlag7 {Int128, UInt8} sixty_four = big"1" --> 64 one_hundred = UInt8(2) --> 100
    @test Integer(TFlag7(sixty_four = true)) === UInt8(1)
    @test Integer(TFlag7(one_hundred = true)) === UInt8(2)
    @test eltype(TFlag7) == Int128
    @test TFlag7([64]) == TFlag7(sixty_four = true)
    @test TFlag7([100]) == TFlag7(one_hundred = true)
    @test 1 ∉ TFlag7([64])
    @test 64 ∈ TFlag7([64])

    # Correctly detect base type
    @flagset TFlag8 {} missing (Int64(1) << 60) --> nothing
    @test eltype(TFlag8) == Union{Missing,Nothing}
    @test typeof(Integer(TFlag8())) == UInt64
    @test typeof(Int(TFlag8())) == Int

    @flagset TFlag9 {Any,_} missing (Int128(1) << 120) --> nothing
    @test eltype(TFlag9) == Any
    @test typeof(Integer(TFlag9())) == UInt128
    @test typeof(Int(TFlag9())) == Int

    # UInt128 and BigInt
    @test typeof(Integer(TFlag4a())) == UInt128
    @test typeof(Integer(TFlag4b())) == BigInt
end # testset

@testset "Validate type" begin
    @test isvalid(TFlag1([RoundDown]))
    @test isvalid(TFlag1([RoundUp]))
    @test isvalid(TFlag1(0x0000_0003))
    @test !isvalid(TFlag1, UInt32(1 << 3))
    @test !isvalid(TFlag1, UInt32(1 << 4))
    @test !isvalid(TFlag1, UInt32(1 << 5 | 1))
    @test !isvalid(TFlag1, UInt32(1 << 29 | 1 << 3))

    @test isvalid(TFlag6(foo = true))
    @test isvalid(TFlag6([Baz]))
    @test isvalid(TFlag6(0x06))
    @test !isvalid(TFlag6, UInt8(1 << 3))
    @test !isvalid(TFlag6, UInt8(1 << 5))
    @test !isvalid(TFlag6, UInt8(1 << 4 | 1))
    @test !isvalid(TFlag6, UInt8(1 << 7 | 1 << 2))

    @test isvalid(TFlag1, [RoundDown])
    @test !isvalid(TFlag1, [RoundNearestTiesAway])
    @test isvalid(TFlag1, [RoundDown, RoundUp])
    @test isvalid(TFlag1, ([RoundDown, RoundUp]))
    @test isvalid(TFlag1, Set([RoundDown, RoundUp]))
    @test TFlag1(Set([RoundDown, RoundUp])) == TFlag1([RoundDown, RoundUp])
    @test !isvalid(TFlag1, [:foo, :bar])
end

@testset "Error conditions" begin
    @test_throws(
        ArgumentError("no arguments given for FlagSet Foo"),
        @macrocall(@flagset Foo),
    )
    @test_throws(
        invalid_type_error(:base, :Foo, :Float64),
        @macrocall(@flagset Foo {_, Float64} x = 1.0),
    )
    @test_throws(
        invalid_type_error(:flag, :Foo, 1234),
        @macrocall(@flagset Foo {1234} x = 1.0),
    )
    # Require uniqueness
    @test_throws(
        not_unique_error(:bit, :Foo, 1),
        @macrocall(@flagset Foo x = 1 --> :x y = 1 --> :y),
    )
    @test_throws(
        not_unique_error(:bit, :Foo, 2),
        @macrocall(@flagset Foo x = 2 --> :x y = 2 --> :y),
    )
    # Explicit bits must be powers of two
    @test_throws(
        invalid_bit_error(:Foo, :(_three = 3 --> nothing)),
        @macrocall(@flagset Foo _three = 3 --> nothing),
    )
    # Values must be integers
    @test_throws(
        invalid_bit_error(:Foo, :(_zero = "zero" --> nothing)),
        @macrocall(@flagset Foo _zero = "zero" --> nothing),
    )
    @test_throws(
        invalid_bit_error(:Foo, :('0' --> nothing)),
        @macrocall(@flagset Foo '0' --> nothing),
    )
    @test_throws(
        invalid_bit_error(:Foo, :(0.5 --> nothing)),
        @macrocall(@flagset Foo 0.5 --> nothing),
    )
    # Names must be valid identifiers
    @test_throws(
        invalid_flagspec_error(:Foo, :(1 = 2)),
        @macrocall(@flagset Foo 1 = 2),
    )
    # Disallow bit overflow
    @test_throws(
        overflow_error(:Foo, :(y = nothing)),
        @macrocall(@flagset Foo {nothing, UInt32} x = UInt32(1) << 31 --> missing y = nothing),
    )
end # testset

@testset "Input/Output" begin
    @flagset TFlag10 {nothing, UInt64} foo = 1 --> Foo bar = 256 --> Bar
    @flagset TFlag11 {nothing, UInt8} foo = 128 --> Foo bar = 2 --> Bar
    let io = IOBuffer()
        write(io, TFlag11([Foo, Bar]))
        write(io, TFlag10(bar = true))
        write(io, TFlag11([Foo]))
        seekstart(io)
        @test read(io, TFlag11) == TFlag11([Foo, Bar])
        @test read(io, TFlag10) == TFlag10(bar = true)
        @test read(io, TFlag11) == TFlag10([Foo])
    end
    let io = IOBuffer()
        serialize(io, TFlag10([Foo]))
        seekstart(io)
        @test deserialize(io) === TFlag10([Foo])
    end
end # testset

# In submodule and not imported to test prefix printing
module TSubModule
using ..FlagSets
@flagset TBits {Int, UInt8} 1 2 3 4
end

@testset "String representations" begin
    @flagset TFilePerms {String, UInt8} 4 --> "READ" 2 --> "WRITE" 1 --> "EXEC"

    @test string(TFilePerms) == "TFilePerms"
    @test string(TSubModule.TBits) == "Main.TSubModule.TBits"
    @test string(TFilePerms()) == "TFilePerms([])"
    @test string(TSubModule.TBits([1])) == "Main.TSubModule.TBits([1])"
    @test repr("text/plain", TFilePerms) ==
          "FlagSet $(string(TFilePerms)):\n" *
          " 0x01 --> \"EXEC\"\n" *
          " 0x02 --> \"WRITE\"\n" *
          " 0x04 --> \"READ\""
    @test repr("text/plain", TSubModule.TBits) ==
          "FlagSet Main.TSubModule.TBits:\n" *
          " 0x01 --> 1\n" *
          " 0x02 --> 2\n" *
          " 0x04 --> 3\n" *
          " 0x08 --> 4"
    @test repr(TFilePerms(["EXEC"])) == "TFilePerms([\"EXEC\"])"
    @test repr(TSubModule.TBits([1])) == "Main.TSubModule.TBits([1])"
    @test repr(TFilePerms(["EXEC", "READ"])) == "TFilePerms([\"EXEC\", \"READ\"])"
    @test repr(TSubModule.TBits([1, 4])) == "Main.TSubModule.TBits([1, 4])"
    @test repr(TFilePerms() | TFilePerms(["READ"])) == "TFilePerms([\"READ\"])"

    let io = IOBuffer(), ioc = IOContext(io, :compact => true, :module => Main)
        show(io, MIME"text/plain"(), FlagSet)
        @test String(take!(io)) == "FlagSet"
        show(ioc, MIME"text/plain"(), FlagSet)
        @test String(take!(io)) == "FlagSet"

        # Explicit :compact --> false required for consistency across Julia versions
        iof = IOContext(io, :compact => false, :module => Main)
        show(iof, MIME"text/plain"(), Union{TFilePerms,TSubModule.TBits})
        @test String(take!(io)) == "Union{TFilePerms, Main.TSubModule.TBits}"
        show(ioc, MIME"text/plain"(), Union{TFilePerms,TSubModule.TBits})
        @test String(take!(io)) == "Union{TFilePerms, TBits}"

        show(ioc, TFilePerms())
        @test String(take!(io)) == "TFilePerms(0x00)"
        show(ioc, TSubModule.TBits([1]))
        @test String(take!(io)) == "TBits(0x01)"
        show(ioc, TFilePerms(["EXEC", "READ"]))
        @test String(take!(io)) == "TFilePerms(0x05)"
        show(ioc, TSubModule.TBits([1, 4]))
        @test String(take!(io)) == "TBits(0x09)"
    end
end # testset


text:
abstract type NodeType end
struct InputNode <: NodeType end
struct HiddenNode <: NodeType end
struct OutputNode <: NodeType end
struct ExternalNode <: NodeType end

mutable struct Node
    string::String

    in_inhibitory::Vector{Node}
    out_inhibitory::Vector{Node}

    in_excitatory::Vector{Node}
    out_excitatory::Vector{Node}

    temp_in_inhibitory::Vector{Node}
    temp_out_inhibitory::Vector{Node}
    temp_in_excitatory::Vector{Node}
    temp_out_excitatory::Vector{Node}

    type::NodeType

    inhibitory_output::Float64
    excitatory_output::Float64

    temp_inhibitory_output::Float64
    temp_excitatory_output::Float64

    reachable_from_input::Bool

    io_value::Float64

    function Node(string, in_inhibitory, out_inhibitory, in_excitatory, out_excitatory)
        new(string, in_inhibitory, out_inhibitory, in_excitatory, out_excitatory,
            [], [], [], [], HiddenNode(), 0, 0, 0, 0, false, 0)
    end

    Node(string) = new(string, [], [], [], [], [], [], [], [], HiddenNode(), 0, 0, 0, 0, false, 0)
end

function fill_temp_links!(node)
    node.temp_in_inhibitory = []
    node.temp_out_inhibitory = []
    node.temp_in_excitatory = []
    node.temp_out_excitatory = []
    append!(node.temp_in_inhibitory, node.in_inhibitory)
    append!(node.temp_out_inhibitory, node.out_inhibitory)
    append!(node.temp_in_excitatory, node.in_excitatory)
    append!(node.temp_out_excitatory, node.out_excitatory)
end

function apply_temp_links!(node)
    node.in_inhibitory = []
    node.out_inhibitory = []
    node.in_excitatory = []
    node.out_excitatory = []

    append!(node.in_inhibitory, node.temp_in_inhibitory)
    append!(node.out_inhibitory, node.temp_out_inhibitory)
    append!(node.in_excitatory, node.temp_in_excitatory)
    append!(node.out_excitatory, node.temp_out_excitatory)
end

function apply_temp_activations!(node)
    node.inhibitory_output = node.temp_inhibitory_output
    node.excitatory_output = node.temp_excitatory_output
end

max_stars_repo_name: tochsner/ALifeBenchmark

max_issues_repo_name: tochsner/ALifeBenchmark

max_forks_repo_name: tochsner/ALifeBenchmark


text:
#!/usr/bin/env julia -- -p4
using Printf
using Distributed

function ζ(N::Int64,a::Int64)
  s = 0.0
  e = convert(Float64,a)

  #for i ∈ 1:N
  #  s += (1.0/i)^e
  #end

  @distributed (+) for i ∈ 1:N
     (1.0/i)^e
  end

end

#N = 16000000000
N = 1600000000
@time ζ_2 = ζ(N,2)
π_squared_over_6 = π^2 / 6.0;

@printf "      ζ(2) = %18.16f\n" ζ_2
@printf "π^2 / 6.0  = %18.16f\n" π_squared_over_6
@printf "     error = %18.16f\n" π_squared_over_6 - ζ_2

max_stars_repo_name: joelandman/rzf

max_issues_repo_name: joelandman/rzf

max_forks_repo_name: joelandman/rzf


text:
SubWorker{H,T} = RemoteChannel{Channel{Vector{SubProblem{H,T}}}}
ScenarioProblemChannel{S} = RemoteChannel{Channel{ScenarioProblems{S}}}
Work = RemoteChannel{Channel{Int}}

function load_subproblems!(subworkers::Vector{SubWorker{H,T}},
                           scenarioproblems::DistributedScenarioProblems,
                           decisions::Vector{DecisionChannel}) where {H <: AbstractFeasibilityHandler,
                                                                      T <: AbstractFloat}
    # Create subproblems on worker processes
    @sync begin
        for w in workers()
            subworkers[w-1] = RemoteChannel(() -> Channel{Vector{SubProblem{H,T}}}(1), w)
            prev = map(2:(w-1)) do p
                scenarioproblems.scenario_distribution[p-1]
            end
            start_id = isempty(prev) ? 0 : sum(prev)
            @async remotecall_fetch(initialize_subworker!,
                                    w,
                                    subworkers[w-1],
                                    scenarioproblems[w-1],
                                    decisions[w-1],
                                    start_id)
        end
    end
end

function initialize_subworker!(subworker::SubWorker{H,T},
                               scenarioproblems::ScenarioProblemChannel,
                               decisions::DecisionChannel,
                               start_id::Integer) where {H <: AbstractFeasibilityHandler, T <: AbstractFloat}
    sp = fetch(scenarioproblems)
    subproblems = Vector{SubProblem{H,T}}(undef, num_subproblems(sp))
    for i in 1:num_subproblems(sp)
        subproblems[i] = SubProblem(
            subproblem(sp, i),
            start_id + i,
            T(probability(scenario(sp, i))),
            all_known_decisions(fetch(decisions)),
            H)
    end
    put!(subworker, subproblems)
    return nothing
end

function restore_subproblems!(subworkers::Vector{<:SubWorker})
     @sync begin
        for w in workers()
            @async remotecall_fetch(w, subworkers[w-1]) do sw
                for subproblem in fetch(sw)
                    restore_subproblem!(subproblem)
                end
            end
        end
    end
    return nothing
end

function resolve_subproblems!(subworker::SubWorker{H,T},
                              decisions::DecisionChannel,
                              x::AbstractVector,
                              cutqueue::CutQueue{T},
                              aggregator::AbstractAggregator,
                              t::Integer,
                              metadata::MetaData) where {H <: AbstractFeasibilityHandler, T <: AbstractFloat}
    # Fetch all subproblems stored in worker
    subproblems::Vector{SubProblem{H,T}} = fetch(subworker)
    if isempty(subproblems)
        # Workers has nothing do to, return.
        return nothing
    end
    # Update subproblems
    update_known_decisions!(fetch(decisions), x)
    # Aggregation policy
    aggregation::AbstractAggregation = aggregator(length(subproblems), T)
    # Solve subproblems
    for subproblem in subproblems
        update_subproblem!(subproblem)
        cut::SparseHyperPlane{T} = subproblem(x)
        aggregate_cut!(cutqueue, aggregation, metadata, t, cut, x)
    end
    flush!(cutqueue, aggregation, metadata, t, x)
    return nothing
end

function work_on_subproblems!(subworker::SubWorker{H,T},
                              decisions::DecisionChannel,
                              work::Work,
                              finalize::Work,
                              cutqueue::CutQueue{T},
                              iterates::RemoteIterates{A},
                              metadata::MetaData,
                              aggregator::AbstractAggregator) where {H <: AbstractFeasibilityHandler,
                                                                     T <: AbstractFloat,
                                                                     A <: AbstractVector}
    subproblems::Vector{SubProblem{H,T}} = fetch(subworker)
    if isempty(subproblems)
       # Workers has nothing do to, return.
       return nothing
    end
    aggregation::AbstractAggregation = aggregator(length(subproblems), T)
    quit = false
    while true
        t::Int = try
            if isready(finalize)
                quit = true
                take!(finalize)
            else
                wait(work)
                take!(work)
            end
        catch err
            if err isa InvalidStateException
                # Master closed the work/finalize channel. Worker finished
                return nothing
            end
        end
        t == -1 && continue
        x::A = fetch(iterates, t)
        # Update subproblems
        update_known_decisions!(fetch(decisions), x)
        for subproblem in subproblems
            update_subproblem!(subproblem)
            cut::SparseHyperPlane{T} = subproblem(x)
            !quit && aggregate_cut!(cutqueue, aggregation, metadata, t, cut, x)
        end
        !quit && flush!(cutqueue, aggregation, metadata, t, x)
        if quit
            # Worker finished
            return nothing
        end
    end
end


text:
# Test MATLAB MAT file I/O

using MATLAB, Compat, Base.Test
import Compat: ASCIIString

a = Int32[1 2 3; 4 5 6]
b = [1.2, 3.4, 5.6, 7.8]
c = Any[[0., 1.], [1., 2.], [1., 2., 3.]]
d = @compat Dict{Any,Any}("name"=>"MATLAB", "score"=>100.)

immutable S
    x::Float64
    y::Bool
    z::Vector{Float64}
end

ss = S[S(1.0, true, [1., 2.]), S(2.0, false, [3., 4.])]

write_matfile("test.mat"; a=a, b=b, c=c, d=d, ss=mxstructarray(ss))

r = read_matfile("test.mat")
@test isa(r, Dict{ASCIIString, MxArray})
@test length(r) == 5

ra = jmatrix(r["a"])
rb = jvector(r["b"])
rc = jvariable(r["c"])
rd = jdict(r["d"])
rss = r["ss"]

gc()  # make sure that ra, rb, rc, rd remain valid

@test ra == a
@test rb == b
@test rc == c

@test rd["name"] == d["name"]
@test rd["score"] == d["score"]

@test is_struct(rss)
@test jscalar(get_field(rss, 1, "x")) == 1.0
@test jscalar(get_field(rss, 1, "y"))
@test jvector(get_field(rss, 1, "z")) == ss[1].z
@test jscalar(get_field(rss, 2, "x")) == 2.0
@test !jscalar(get_field(rss, 2, "y"))
@test jvector(get_field(rss, 2, "z")) == ss[2].z


text:
"""
--------------------------------------------------------------------------------------\n
Author: Huckleberry Febbo, Graduate Student, University of Michigan
Date Create: 2017, Last Modified: 2/06/2018 \n
--------------------------------------------------------------------------------------\n
"""
function checkCrash(n; kwargs...)
    kw = Dict(kwargs)  # NOTE kwargs are depreciated

    c = n.ocp.params[5]
    sm = c["misc"]["sm"]
    if !haskey(kw,:plant); plant=true;
    else; plant = get(kw,:plant,0);
    end
    @unpack_Vpara n.ocp.params[1]

    if isequal(c["misc"]["model"],:ThreeDOFv2) || isequal(c["misc"]["model"],:ThreeDOFv3)
        # check to see if the minimum vertical tire load was exceeded
        Fz_off = c["vehicle"][:Fz_off]  # should be checked with mode
        if isequal(c["misc"]["model"],:ThreeDOFv2)
            V = n.r.ip.dfsplant[end][:v]
            U = n.r.ip.dfsplant[end][:ux]
            Ax = n.r.ip.dfsplant[end][:ax]
            R = n.r.ip.dfsplant[end][:r]
            SA = n.r.ip.dfsplant[end][:sa]
        elseif isequal(c["misc"]["model"],:ThreeDOFv3)
            V = n.r.ip.dfsplant[end][:v]
            U = u0_*ones(length(V),1)
            Ax = zeros(length(V),1)
            R = n.r.ip.dfsplant[end][:r]
            SA = n.r.ip.dfsplant[end][:sa]
        end

        if any(@FZ_RL() .< Fz_off) || any(@FZ_RR() .< Fz_off) || any(@FZ_FR() .< Fz_off) || any(@FZ_RL() .< Fz_off)
            println("The vertical tire force went below Fz_off")
            return true, :tireOff
        end
    end
    # check to see if the vehicle crashed into and obstacle
    if plant
        t = n.r.ip.plant[:t]
        X = n.r.ip.plant[:x]
        Y = n.r.ip.plant[:y]
        crash_tmp = zeros(length(c["obstacle"]["radius"]),1)
        for obs in 1:length(c["obstacle"]["radius"])
            # obstacle postions after the initial postion
            X_obs= c["obstacle"]["x0"][obs] .+ c["obstacle"]["vx"][obs].*t
            Y_obs= c["obstacle"]["y0"][obs] .+ c["obstacle"]["vy"][obs].*t
            if minimum((X-X_obs).^2./(c["obstacle"]["radius"][obs]+sm).^2 + (Y-Y_obs).^2./(c["obstacle"]["radius"][obs]+sm).^2) < 1
                crash_tmp[obs] = 1
                println("the vehicle crashed! \n")
            end
        end
        if maximum(crash_tmp)>0
            return true, :crash
        else
            return false, NaN
        end
    else # in the NLOptControl.jl paper, there was no plant and the following code was used
        # also this will break, need to remove pts from function handle
       pts = 300
       # redo interpolation with desired number of points
      if n.s.integrationMethod == :ps
          interpolateLagrange!(n;numPts=Int64(pts/n.Ni))
      else
          interpolateLinear!(n;numPts=pts)
      end
      X = n.r.X_pts[:,1]; Y = n.r.X_pts[:,2];

      crash_tmp = zeros(pts,1);
      for i = 1:pts
        temp = (X[i]-c.o.X0[1])^2/(c.o.B[1]+sm)^2 + (Y[i]-c.o.Y0[1])^2/(c.o.A[1]+sm)^2;
       if temp < 1
          crash_tmp[i] = 1;
       else
          crash_tmp[i] = 0;
       end
      end

      if maximum(crash_tmp)>0
          crash = 1;
          print("the vehicle crashed! \n")
      else
          crash = 0;
      end
     return crash
  end
end

"""
# this funtion can be used to get the fieldnames in a parameter set
# example
FNames(Vpara())
--------------------------------------------------------------------------------------\n
Author: Huckleberry Febbo, Graduate Student, University of Michigan
Date Create: 4/08/2018, Last Modified: 4/08/2018 \n
--------------------------------------------------------------------------------------\n
"""
function FNames(pa)
# To find fieldnames
r = ()
for i in 1:length(fieldnames(p))
      r = (r...,fieldnames(p)[i])
end

return r
end


text:
# A minimalist's implementation of the ideas described in
# A Minimalist's Implementation of an Approximate Nearest Neighbor Algorithm in Fixed Dimensions.
# Paper: https://cs.uwaterloo.ca/~tmchan/sss.ps
# Additional references:
#  - http://en.wikipedia.org/wiki/Z-order_curve#Efficiently_building_quadtrees

lessmsb(x, y) = x < y && x < (x $ y)
function shuffdim(p, q)
	# Return the dimension with the largest most significant bit, which
	# decides which of p and q comes first in shuffle order.
	# Break ties in favor of the lower-index dimension (x, then y, then z)
	# since that's the order in which the bits are conceptually interleaved.
	k = 1
	kxor = p[k] $ q[k]
	for i in 2:min(size(p, 1), size(q, 1)) # Workaround for wonky bbox point stuff.
		ixor = p[i] $ q[i]
		if lessmsb(kxor, ixor)
			k, kxor = i, ixor
		end
	end
	k
end
shuffless(p, q) = (k = shuffdim(p, q); p[k] < q[k])
shuffmore(p, q) = (k = shuffdim(p, q); p[k] > q[k])
  shuffeq(p, q) = (k = shuffdim(p, q); p[k] == q[k])

sq(v) = v * v

function dist_sq_to_box(q, min, max)
	# Return the distance from q to the power-of-two-sized quadtree box containing min and max.
	k = shuffdim(min, max)

	# Handle the case when the bounding box is a single point separately.
	# Failing to do so will trigger a DomainError when calculating
	# i below, since the xor of two equal elements is zero.
	min[k] == max[k] && return dist_sq(q, min)::Uint

	# i is the 1-based index of the maximum differing bit
	# across min and max. It is also to the power of 2 of the
	# size of the quadtree box containing both min and max.
	i = 1 + exponent(float(min[k] $ max[k]))

	# We can compute the squared distance from q to the bounding box
	# by accumulating separately in each dimension
	d_sq::Uint = 0
	for j in 1:size(min, 1)
		# The left and right coordinates of the bounding box in this dimension
		lo = (min[j] >> i) << i # equivalent to floor(min[j] / 2^i) * 2^i
		hi = lo + (1 << i)      # equivalent to lo + 2^i

		# If q is outside the bounds, add to d_sq the amount by which it is so.
		if q[j] < lo
			d_sq += sq(q[j] - lo)
		elseif q[j] > hi
			d_sq += sq(q[j] - hi)
		end
	end
	d_sq
end

# Note: These can't stay Uint8s forever.
immutable BBoxPoint # {4}
	x::Uint8
	y::Uint8
	z::Uint8
	w::Uint8
end
Base.getindex(p::BBoxPoint, n::Int) = n == 1 ? p.x : n == 2 ? p.y : n == 3 ? p.z : n == 4 ? p.w : throw("BBoxPoint indexing error.")
Base.size(p::BBoxPoint) = (4,)
Base.size(p::BBoxPoint, n) = n == 1 ? 4 : throw("Invalid dimension.")

# Note: Due to the PointN types, this now assumes all coordinates are of the same type.
clampadd3{Q}(p::Q, r::Int) = BBoxPoint(min(p[1] + r, typemax(p[1])), min(p[2] + r, typemax(p[2])), min(p[3] + r, typemax(p[3])), 0)
clampsub3{Q}(p::Q, r::Int) = BBoxPoint(max(p[1] - r, typemin(p[1])), max(p[2] - r, typemin(p[2])), max(p[3] - r, typemin(p[3])), 0)
clampadd4{Q}(p::Q, r::Int) = BBoxPoint(min(p[1] + r, typemax(p[1])), min(p[2] + r, typemax(p[2])), min(p[3] + r, typemax(p[3])), min(p[4] + r, typemax(p[4])))
clampsub4{Q}(p::Q, r::Int) = BBoxPoint(max(p[1] - r, typemin(p[1])), max(p[2] - r, typemin(p[2])), max(p[3] - r, typemin(p[3])), max(p[4] - r, typemin(p[4])))

representative(p) = p

preprocess!(P) = sort!(P, lt=shuffless)

immutable Result{P, Q}
	r_sq::Uint
	point::P
	bbox_min::BBoxPoint # The bounding box enclosing the hypersphere with radius r
	bbox_max::BBoxPoint # centered on the query point (with side length 2r).
	function Result(r_sq::Uint, point::P, q::Q)
		r = iceil(sqrt(r_sq))
		new(r_sq, point, clampsub3(q, r), clampadd3(q, r))
	end
end

# P = point; Q = query
function nearest{P, Q}(t::TreapNode{P}, q::Q, R::Result{P, Q}, ε=0.0)
	isempty(t) && return R
	point = representative(key(t))
	r_sq = dist_sq(point, q)
	r_sq < R.r_sq && (R = Result{P, Q}(r_sq, key(t), q))

	if length(t) == 1 || dist_sq_to_box(q, representative(minimum(t)), representative(maximum(t))) * sq(1 + ε) > R.r_sq
		return R
	end

	if shuffless(q, point)
		R = nearest(left(t), q, R, ε)
		shuffmore(R.bbox_max, point) && (R = nearest(right(t), q, R, ε))
	else
		R = nearest(right(t), q, R, ε)
		shuffless(R.bbox_min, point) && (R = nearest(left(t), q, R, ε))
	end

	R
end

function nearest{P, Q}(t::Treap{P}, q::Q, ε=0.0)
	R = nearest(root(t), q, Result{P, Q}(uint(9999999), key(t), q), ε)
	@assert R.r_sq < uint(9999999)
	R.point
end

# BUG: This assumes you've already imported rgb.jl for no good reason.
# Define less-than and distances for pixels for testing purposes.
<(p::Pixel, q::Pixel) = shuffless(p, q)
dist_sq(p::Pixel, q::Pixel) = sq(p.r - q.r) + sq(p.g - q.g) + sq(p.b - q.b)

function test_sss(n=1000, ε=0.0)

	function linear_nearest(P, q)
		best = q
		best_dist::Float64 = Inf
		for p in P
			d = dist_sq(p, q)
			if d < best_dist
				best = p
				best_dist = d
			end
		end
		best
	end

	P = Treap{Pixel}()
	pixels = unique([Pixel(rand(Uint8), rand(Uint8), rand(Uint8)) for i in 1:1500])
	for p in pixels
		add!(P, p)
	end

	qs = [Pixel(rand(Uint8), rand(Uint8), rand(Uint8)) for i in 1:n]
	lps = [linear_nearest(P, q) for q in qs]
	min_distances = [sqrt(dist_sq(qs[i], lps[i])) for i in 1:n]

	num_matches = 0
	@time for i in 1:n
		q = qs[i]
		tp = nearest(P, q, ε)
		tree_distance = sqrt(dist_sq(q, tp))
		min_distance = min_distances[i]
		if min_distance != tree_distance
			lp = lps[i]
			println("Mismatch for $q: $tp, $lp. ($(int(min_distance)) vs. $(int(tree_distance))); $(int(100 * min_distance / tree_distance))%")
		else
			num_matches += 1
		end
	end
	println("SSS: $num_matches / $n matches.")
end

# test_sss()

max_stars_repo_name: yurivish/Julia-Playground

max_issues_repo_name: yurivish/Julia-Playground

max_forks_repo_name: yurivish/Julia-Playground


text:
using TimerOutputs, Printf
using Oceananigans

const timer = TimerOutput()

Ni = 2   # Number of iterations before benchmarking starts.
Nt = 10  # Number of iterations to use for benchmarking time stepping.

# Model resolutions to benchmarks. Focusing on 3D models for GPU.
Ns = [# (1, 1, 512),
      # (1, 128, 128), (128, 1, 128), (128, 128, 1),
      (32, 32, 32), (64, 64, 64), (128, 128, 128), (256, 256, 256)]

float_types = [Float32, Float64]  # Float types to benchmark.
archs = [CPU()]  # Architectures to benchmark on.

@hascuda archs = [CPU(), GPU()]  # Benchmark GPU on systems with CUDA-enabled GPUs.

arch_name(::CPU) = "CPU"
arch_name(::GPU) = "GPU"

benchmark_name(N, id)               = benchmark_name(N, id, nothing, nothing)
benchmark_name(N, id, arch::Symbol) = benchmark_name(N, id, arch, nothing)
benchmark_name(N, id, ft::DataType) = benchmark_name(N, id, nothing, ft)

function benchmark_name(N, id, arch, ft; npad=3)
    Nx, Ny, Nz = N
    print_arch = typeof(arch) <: Architecture ? true : false
    print_ft   = typeof(ft) == DataType && ft <: AbstractFloat ? true : false

    bn = ""
    bn *= lpad(Nx, npad, " ") * "x" * lpad(Ny, npad, " ") * "x" * lpad(Nz, npad, " ")
    bn *= " $id"

    if print_arch && print_ft
        arch = arch_name(arch)
        bn *= " ($arch, $ft)"
    elseif print_arch && !print_ft
        arch = arch_name(arch)
        bn *= " ($arch)"
    elseif !print_arch && print_ft
        bn *= " ($ft)"
    end

    return bn
end

for arch in archs, float_type in float_types, N in Ns
    Nx, Ny, Nz = N
    Lx, Ly, Lz = 100, 100, 100

    model = Model(N=(Nx, Ny, Nz), L=(Lx, Ly, Lz), arch=arch, float_type=float_type)
    time_step!(model, Ni, 1)  # First 1-2 iterations usually slower.

    bn =  benchmark_name(N, "static ocean", arch, float_type)
    @printf("Running benchmark: %s...\n", bn)
    for i in 1:Nt
        @timeit timer bn time_step!(model, 1, 1)
    end
end

print_timer(timer, title="Oceananigans.jl benchmarks")

bid = "static ocean"  # Benchmark ID. We only have one right now.

println("\n\nCPU Float64 -> Float32 speedup:")
for N in Ns
    bn32 = benchmark_name(N, bid, CPU(), Float32)
    bn64 = benchmark_name(N, bid, CPU(), Float64)
    t32  = TimerOutputs.time(timer[bn32])
    t64  = TimerOutputs.time(timer[bn64])
    @printf("%s: %.3f\n", benchmark_name(N, bid), t64/t32)
end

@hascuda begin
    println("\nGPU Float64 -> Float32 speedup:")
    for N in Ns
        bn32 = benchmark_name(N, bid, GPU(), Float32)
        bn64 = benchmark_name(N, bid, GPU(), Float64)
        t32  = TimerOutputs.time(timer[bn32])
        t64  = TimerOutputs.time(timer[bn64])
        @printf("%s: %.3f\n", benchmark_name(N, bid), t64/t32)
    end

    println("\nCPU -> GPU speedup:")
    for N in Ns, ft in float_types
        bn_cpu = benchmark_name(N, bid, CPU(), ft)
        bn_gpu = benchmark_name(N, bid, GPU(), ft)
        t_cpu  = TimerOutputs.time(timer[bn_cpu])
        t_gpu  = TimerOutputs.time(timer[bn_gpu])
        @printf("%s: %.3f\n", benchmark_name(N, bid, ft), t_cpu/t_gpu)
    end
end


text:
"""
A basis of Chebyshev polynomials of the first kind on the interval `[-1,1]`.
"""
struct ChebyshevT{T} <: OPS{T}
    n			::	Int
end

show(io::IO, d::ChebyshevT{Float64}) = print(io, "ChebyshevT($(length(d)))")
show(io::IO, d::ChebyshevT{T}) where T = print(io, "ChebyshevT{$(T)}($(length(d)))")

ChebyshevT(n::Int) = ChebyshevT{Float64}(n)

@deprecate ChebyshevT(n, a::Number, b::Number) ChebyshevT(n)→a..b

similar(b::ChebyshevT, ::Type{T}, n::Int) where {T} = ChebyshevT{T}(n)

convert(::Type{ChebyshevT{T}}, d::ChebyshevT) where {T} = ChebyshevT{T}(d.n)

"Mapped Chebyshev polynomials."
struct MappedChebyshev{T} <: MappedDict{T,T}
	superdict	::	ChebyshevT{T}
	map			::	ScalarAffineMap{T}
end

mapped_dict(dict::ChebyshevT{T}, map::ScalarAffineMap{S}) where {S,T} =
	MappedChebyshev{promote_type(S,T)}(dict, map)

hasinterpolationgrid(b::ChebyshevT) = true
hasderivative(b::ChebyshevT) = true
hasderivative(b::ChebyshevT, order::Int) = true
hasantiderivative(b::ChebyshevT) = true

hasgrid_transform(b::ChebyshevT, gb, ::ChebyshevNodes) = length(b) == length(gb)
hasgrid_transform(b::ChebyshevT, gb, ::ChebyshevExtremae) = length(b) == length(gb)
hasgrid_transform(b::ChebyshevT, gb, ::AbstractGrid) = false

first_moment(b::ChebyshevT{T}) where {T} = convert(T, pi)

interpolation_grid(b::ChebyshevT{T}) where {T} = ChebyshevNodes{T}(length(b))
iscompatible(dict::ChebyshevT, grid::ChebyshevNodes) = length(dict) == length(grid)
secondgrid(b::ChebyshevT{T}) where {T} = ChebyshevExtremae{T}(length(b))
transformgrid_extremae(b::ChebyshevT{T}) where {T} = ChebyshevExtremae{T}(length(b))

# extends the default definition at transform.jl
function transform_dict(s::ChebyshevT{T}; chebyshevpoints=:nodes, options...) where {T}
	if chebyshevpoints == :nodes
    	GridBasis(s)
	else
		@assert chebyshevpoints == :extremae
		GridBasis{T}(secondgrid(s))
	end
end

# Parameters alpha and beta of the corresponding Jacobi polynomial
jacobi_α(b::ChebyshevT{T}) where {T} = -one(T)/2
jacobi_β(b::ChebyshevT{T}) where {T} = -one(T)/2

# See DLMF, Table 18.9.1
# http://dlmf.nist.gov/18.9#i
rec_An(b::ChebyshevT{T}, n::Int) where {T} = n==0 ? one(T) : 2one(T)
rec_Bn(b::ChebyshevT{T}, n::Int) where {T} = zero(T)
rec_Cn(b::ChebyshevT{T}, n::Int) where {T} = one(T)

support(b::ChebyshevT{T}) where {T} = ChebyshevInterval{T}()

# We can define this O(1) evaluation method, but only for points that are
# real and lie in [-1,1]
# Note that if x is not Real, recurrence_eval will be called by the OPS supertype
function unsafe_eval_element(b::ChebyshevT, idx::PolynomialDegree, x::Real)
    abs(x) <= 1 ? cos(degree(idx)*acos(x)) : recurrence_eval(b, idx, x)
end

# The version below is safe for points outside [-1,1] too.
# If we don't define anything, evaluation will default to using the three-term
# recurence relation.
# unsafe_eval_element{T <: Real}(b::ChebyshevT, idx::Int, x::T) = real(cos((idx-1)*acos(x+0im)))

unsafe_moment(dict::ChebyshevT, idx; measure = lebesguemeasure(support(dict)), options...) =
	unsafe_moment(dict, idx, measure; options...)

unsafe_moment(dict::ChebyshevT, idx, measure; options...) =
	innerproduct(dict, idx, dict, PolynomialDegree(0), measure; options...)

function unsafe_moment(dict::ChebyshevT{T}, idx, ::LegendreWeight; options...) where {T}
    n = degree(idx)
    if n == 0
        T(2)
    else
        isodd(n) ? zero(T) : -T(2)/((n+1)*(n-1))
    end
end

## Inner products

hasmeasure(dict::ChebyshevT) = true
measure(dict::ChebyshevT{T}) where T = ChebyshevTWeight{T}()
issymmetric(::ChebyshevT) = true

innerproduct_native(b1::ChebyshevT, i::PolynomialDegree, b2::ChebyshevT, j::PolynomialDegree, m::ChebyshevTWeight;
			T = coefficienttype(b1), options...) =
	innerproduct_chebyshev_full(i, j, T)

function innerproduct_chebyshev_full(i, j, T)
	if i == j
		if i == PolynomialDegree(0)
			convert(T, pi)
		else
			convert(T, pi)/2
		end
	else
		zero(T)
	end
end

function innerproduct_native(b1::ChebyshevT, i::PolynomialDegree, b2::ChebyshevT, j::PolynomialDegree, measure::LegendreWeight; options...)
	n1 = degree(i)
	n2 = degree(j)
	(unsafe_moment(b1, PolynomialDegree(n1+n2), measure) + unsafe_moment(b1, PolynomialDegree(abs(n1-n2)), measure))/2
end

###################################################################################
# Methods to transform from ChebyshevT to ChebyshevNodes and ChebyshevExtremae
###################################################################################

evaluation(::Type{T}, dict::ChebyshevT, gb::GridBasis, grid::ChebyshevNodes; options...) where {T} =
	resize_and_transform(T, dict, gb, grid; chebyshevpoints = :nodes, options...)

evaluation(::Type{T}, dict::ChebyshevT, gb::GridBasis, grid::ChebyshevExtremae; options...) where {T} =
	resize_and_transform(T, dict, gb, grid; chebyshevpoints = :extremae, options...)

function chebyshev_transform_nodes(dict::ChebyshevT, T; options...)
	grid = interpolation_grid(dict)
	n = length(dict)
	# We compose the inverse DCT with a diagonal scaling
	F = inverse_chebyshev_operator(dict, GridBasis{T}(grid), T; options...)
	d = zeros(T, n)
	fill!(d, sqrt(T(n)/2))
	d[1] *= sqrt(T(2))
	for i in 1:length(d)
		d[i] *= (-1)^(i-1)
	end
	F * DiagonalOperator(dict, d)
end

function chebyshev_transform_extremae(dict::ChebyshevT, T; options...)
	grid = transformgrid_extremae(dict)
	F = inverse_chebyshevI_operator(dict, GridBasis{T}(grid), T; options...)
	d = zeros(T, length(dict))
	fill!(d, one(T)/2)
	d[1] *= 2
	d[end] *= 2
	F * DiagonalOperator(dict, d)
end

function transform_to_grid(T, src::ChebyshevT, dest::GridBasis, grid::ChebyshevNodes; options...)
	@assert length(src) == length(grid)
	chebyshev_transform_nodes(src, T; options...)
end

function transform_to_grid(T, src::ChebyshevT, dest::GridBasis, grid::ChebyshevExtremae; options...)
	@assert length(src) == length(grid)
	chebyshev_transform_extremae(src, T; options...)
end

transform_from_grid(T, src::GridBasis, dest::ChebyshevT, grid; options...) =
	inv(transform_to_grid(T, dest, src, grid; options...))

iscompatible(src1::ChebyshevT, src2::ChebyshevT) = true

function (*)(src1::ChebyshevT, src2::ChebyshevT, coef_src1, coef_src2)
    @assert domaintype(src1) == domaintype(src2)
    T = promote_type(eltype(coef_src1), eltype(coef_src2))
    dest = ChebyshevT{T}(length(src1)+length(src2))
    coef_dest = zeros(dest)
    for i = 1:length(src1)
        for j = 1:length(src2)
            coef_dest[i+j-1] += one(T)/2*coef_src1[i]*coef_src2[j]
            coef_dest[abs(i-j)+1] += one(T)/2*coef_src1[i]*coef_src2[j]
        end
    end
    (dest,coef_dest)
end

"A Chebyshev polynomial of the first kind."
struct ChebyshevTPolynomial{T} <: OrthogonalPolynomial{T}
    degree  ::  Int
end

ChebyshevTPolynomial{T}(p::ChebyshevTPolynomial) where {T} = ChebyshevTPolynomial{T}(p.degree)

show(io::IO, p::ChebyshevTPolynomial{Float64}) = print(io, "ChebyshevTPolynomial($(p.degree))")

convert(::Type{TypedFunction{T,T}}, p::ChebyshevTPolynomial) where {T} = ChebyshevTPolynomial{T}(p.degree)

support(::ChebyshevTPolynomial{T}) where {T} = ChebyshevInterval{T}()

(p::ChebyshevTPolynomial{T})(x) where {T} = eval_element(ChebyshevT{T}(degree(p)+1), degree(p)+1, x)

basisfunction(dict::ChebyshevT, idx) = basisfunction(dict, native_index(dict, idx))
basisfunction(dict::ChebyshevT{T}, idx::PolynomialDegree) where {T} = ChebyshevTPolynomial{T}(degree(idx))

dictionary(p::ChebyshevTPolynomial{T}) where {T} = ChebyshevT{T}(degree(p)+1)
index(p::ChebyshevTPolynomial) = degree(p)+1


text:
@testset "DirichletBC" begin
    var = Variable{Float64}(1, "u")

    dbc = DirichletBC(var, [2,3], 3.7)

    # Test the BoundaryCondition interface
    @test boundaryIDs(dbc) == [2,3]

    @test dbc.u.id == 1

    @test dbc.value == 3.7

    # Reinit the variable at a "Node"
    MOOSE.reinit!(var, 2, 2.2)

    # Test out the residual and Jacobian calculations
    res = MOOSE.computeResidual(dbc)

    @test abs(res - (2.2 - 3.7)) < 1e-9

    jac = MOOSE.computeJacobian(dbc, dbc.u)

    @test jac == 1.
end

@testset "DiriBCFD" begin
    var = Variable{Dual{4, Float64}}(1, "u")

    dbc = DirichletBC(var, [2,3], 3.7)

    # Test the BoundaryCondition interface
    @test boundaryIDs(dbc) == [2,3]

    @test dbc.u.id == 1

    @test dbc.value == 3.7

    # Reinit the variable at a "Node"
    MOOSE.reinit!(var, 2, MOOSE.dualVariable(2.2, var.id, 4))

    # Test out the residual and Jacobian calculations
    res = MOOSE.computeResidual(dbc)

    @test abs(res - (2.2 - 3.7)) < 1e-9

    jac = MOOSE.computeJacobian(dbc, dbc.u)

    @test jac == 1.
end


text:
"""
integrals on  tetrahedrons
"""

#identical tetrahedron
function k6p_ct(f, x1, x2, x3, y1, y2, y3)
    #sum = 0;

     #Domain 2/17
    #=1=#   ξ1 = x1
    #=2=#   ξ2 = x1 * x2
    #=3=#   ξ3 = x1 * x2 * x3
    #=4=#   η1 = x1 * x2 * x3 * y1
    #=5=#   η2 = x1 * x2 * x3 * y1 * y2
    #=6=#   η3 = x1 * x2 * x3 * y1 * y3

    det = x1^5*x2^4*x3^3*y1^2

    sum = f((1-ξ1,ξ1-ξ2,ξ3-η1+η2),(1-(ξ1-η3),ξ1-ξ2+η1-η3,ξ3-η1))*det

    sum += f((1-(ξ1-η3),ξ1-ξ2+η1-η3,ξ3-η1),(1-ξ1,ξ1-ξ2,ξ3-η1+η2))*det

    #=1=#   ξ1 = x1
    #=2=#   ξ2 = x1 * x2
    #=3=#   ξ3 = x1 * x2 * x3
    #=4=#   η1 = x1 * x2 * x3 * y1
    #=5=#   η2 = x1 * x2 * x3 * y1 * y2
    #=6=#   η3 = x1 * x2 * x3 * y1 * y2 * y3

    det = x1^5*x2^4*x3^3*y1^2*y2

    #Domain 1/18
    sum += f((1-ξ1,ξ1-ξ2,ξ3),(1-(ξ1-η2+η3),ξ1-ξ2+η3,ξ3-η1))*det

    sum += f((1-(ξ1-η2+η3),ξ1-ξ2+η3,ξ3-η1),(1-ξ1,ξ1-ξ2,ξ3))*det

    #Domain 3/16
    sum += f((1-ξ1,ξ1-ξ2,ξ3-η1),(1-(ξ1-η3),ξ1-ξ2+η2-η3,ξ3-η2))*det

    sum += f((1-(ξ1-η3),ξ1-ξ2+η2-η3,ξ3-η2),(1-ξ1,ξ1-ξ2,ξ3-η1),)*det

    #Domain 4/15
    sum += f((1-ξ1,ξ1-ξ2+η2-η3,ξ3-η2+η3),(1-(ξ1-η2),ξ1-ξ2,ξ3-η1))*det

    sum += f((1-(ξ1-η2),ξ1-ξ2,ξ3-η1),(1-ξ1,ξ1-ξ2+η2-η3,ξ3-η2+η3))*det

    #Domain 5/14
    sum += f((1-ξ1,ξ1-ξ2+η1-η2,ξ3-η1+η3),(1-(ξ1-η1),ξ1-ξ2,ξ3-η1))*det

    sum += f((1-(ξ1-η1),ξ1-ξ2,ξ3-η1),(1-ξ1,ξ1-ξ2+η1-η2,ξ3-η1+η3))*det

    #Domain 6/13
    sum += f((1-ξ1,ξ1-ξ2+η3,ξ3-η1),(1-(ξ1-η2),ξ1-ξ2,ξ3-η2))*det

    sum += f((1-(ξ1-η2),ξ1-ξ2,ξ3-η2),(1-ξ1,ξ1-ξ2+η3,ξ3-η1))*det

    #Domain 7/12
    sum += f((1-ξ1,ξ1-ξ2+η1-η3,ξ3-η1+η3),(1-(ξ1-η2+η3),ξ1-ξ2,ξ3-η1))*det

    sum += f((1-(ξ1-η2+η3),ξ1-ξ2,ξ3-η1),(1-ξ1,ξ1-ξ2+η1-η3,ξ3-η1+η3))*det

    #Domain 8/11
    sum += f((1-ξ1,ξ1-ξ2+η1,ξ3-η1),(1-(ξ1-η3),ξ1-ξ2,ξ3-η1+η2-η3))*det

    sum += f((1-(ξ1-η3),ξ1-ξ2,ξ3-η1+η2-η3),(1-ξ1,ξ1-ξ2+η1,ξ3-η1))*det

    #Domain 9/10
    sum += f((1-ξ1,ξ1-ξ2+η2,ξ3-η1),(1-(ξ1-η2+η3),ξ1-ξ2,ξ3-η2+η3))*det

    sum += f((1-(ξ1-η2+η3),ξ1-ξ2,ξ3-η2+η3),(1-ξ1,ξ1-ξ2+η2,ξ3-η1))*det

    return sum
end

function k6p_ct_A(f, x1, x2, x3, x4, y1, z1)
    #sum = 0;

    #Domain 2/17
    #=1=#   ξ1 = x1
    #=2=#   ξ2 = x2
    #=3=#   ξ3 = x3
    #=4=#   η1 = x4
    #=5=#   η2 = x4 * y1
    #=6=#   η3 = x4 * z1

    det = x4^2

    sum = f((1-ξ1,ξ1-ξ2,ξ3-η1+η2),(1-(ξ1-η3),ξ1-ξ2+η1-η3,ξ3-η1))*det

    sum += f((1-(ξ1-η3),ξ1-ξ2+η1-η3,ξ3-η1),(1-ξ1,ξ1-ξ2,ξ3-η1+η2))*det

    return sum
end

function k6p_ct_B(f, x1, x2, x3, x4, y1, y2)
    #sum = 0;
    #=1=#   ξ1 = x1
    #=2=#   ξ2 = x2
    #=3=#   ξ3 = x3
    #=4=#   η1 = x4
    #=5=#   η2 = x4 * y1
    #=6=#   η3 = x4 * y2

    det = x4^2

    #Domain 1/18
    sum = f((1-ξ1,ξ1-ξ2,ξ3),(1-(ξ1-η2+η3),ξ1-ξ2+η3,ξ3-η1))*det

    sum += f((1-(ξ1-η2+η3),ξ1-ξ2+η3,ξ3-η1),(1-ξ1,ξ1-ξ2,ξ3))*det

    #Domain 3/16
    sum += f((1-ξ1,ξ1-ξ2,ξ3-η1),(1-(ξ1-η3),ξ1-ξ2+η2-η3,ξ3-η2))*det

    sum += f((1-(ξ1-η3),ξ1-ξ2+η2-η3,ξ3-η2),(1-ξ1,ξ1-ξ2,ξ3-η1),)*det

    #Domain 4/15
    sum += f((1-ξ1,ξ1-ξ2+η2-η3,ξ3-η2+η3),(1-(ξ1-η2),ξ1-ξ2,ξ3-η1))*det

    sum += f((1-(ξ1-η2),ξ1-ξ2,ξ3-η1),(1-ξ1,ξ1-ξ2+η2-η3,ξ3-η2+η3))*det

    #Domain 5/14
    sum += f((1-ξ1,ξ1-ξ2+η1-η2,ξ3-η1+η3),(1-(ξ1-η1),ξ1-ξ2,ξ3-η1))*det

    sum += f((1-(ξ1-η1),ξ1-ξ2,ξ3-η1),(1-ξ1,ξ1-ξ2+η1-η2,ξ3-η1+η3))*det

    #Domain 6/13
    sum += f((1-ξ1,ξ1-ξ2+η3,ξ3-η1),(1-(ξ1-η2),ξ1-ξ2,ξ3-η2))*det

    sum += f((1-(ξ1-η2),ξ1-ξ2,ξ3-η2),(1-ξ1,ξ1-ξ2+η3,ξ3-η1))*det

    #Domain 7/12
    sum += f((1-ξ1,ξ1-ξ2+η1-η3,ξ3-η1+η3),(1-(ξ1-η2+η3),ξ1-ξ2,ξ3-η1))*det

    sum += f((1-(ξ1-η2+η3),ξ1-ξ2,ξ3-η1),(1-ξ1,ξ1-ξ2+η1-η3,ξ3-η1+η3))*det

    #Domain 8/11
    sum += f((1-ξ1,ξ1-ξ2+η1,ξ3-η1),(1-(ξ1-η3),ξ1-ξ2,ξ3-η1+η2-η3))*det

    sum += f((1-(ξ1-η3),ξ1-ξ2,ξ3-η1+η2-η3),(1-ξ1,ξ1-ξ2+η1,ξ3-η1))*det

    #Domain 9/10
    sum += f((1-ξ1,ξ1-ξ2+η2,ξ3-η1),(1-(ξ1-η2+η3),ξ1-ξ2,ξ3-η2+η3))*det

    sum += f((1-(ξ1-η2+η3),ξ1-ξ2,ξ3-η2+η3),(1-ξ1,ξ1-ξ2+η2,ξ3-η1))*det

    return sum
end

#overlapping face
function k6p_cf(f, x1, x2, x3, y1, y2, y3)
    #sum = 0;

    #Domain 1
    #=1=#   ξ1 = x1
    #=2=#   ξ2 = x1 * x2
    #=3=#   ξ3 = x1 * x2 * x3
    #=4=#   η1 = x1 * x2 * x3 * y1
    #=5=#   η2 = x1 * x2 * x3 * y1 * y2
    #=6=#   η3 = x1 * x2 * x3 * y3

    det = x1^5*x2^4*x3^3*y1

    #sum += f((ξ1,ξ2,η3),(ξ1-ξ3+η1,ξ2-ξ3+η2,η2))*det
    sum = f((1-ξ1,ξ1-ξ2,ξ3-η3),(1-(ξ1-ξ3+η1),ξ1-ξ2+η2,η1-η2))*det

    #Domain 3/5/13
    #=1=#   ξ1 = x1
    #=2=#   ξ2 = x1 * x2
    #=3=#   ξ3 = x1 * x2 * x3
    #=4=#   η1 = x1 * x2 * x3 * y1
    #=5=#   η2 = x1 * x2 * x3 * y1 * y2
    #=6=#   η3 = x1 * x2 * x3 * y1 * y3

    det = x1^5*x2^4*x3^3*y1^2

    sum += f((1-ξ1,ξ1-ξ2+ξ3-η1,η1-η2),(1-(ξ1-ξ3+η3),ξ1-ξ2,η3))*det

    sum += f((1-(ξ1-ξ3+η1),ξ1-ξ2,η1-η3),(1-ξ1,ξ1-ξ2+ξ3-η2,η2))*det

    sum += f((1-(ξ1-η2),ξ1-ξ2+η1-η2,ξ3-η1),(1-ξ1,ξ1-ξ2,η3))*det

    #Domain 2/4/6-12/14-15
    #=1=#   ξ1 = x1
    #=2=#   ξ2 = x1 * x2
    #=3=#   ξ3 = x1 * x2 * x3
    #=4=#   η1 = x1 * x2 * x3 * y1
    #=5=#   η2 = x1 * x2 * x3 * y1 * y2
    #=6=#   η3 = x1 * x2 * x3 * y1 * y2 * y3

    det = x1^5*x2^4*x3^3*y1^2*y2

    sum += f((1-ξ1,ξ1-ξ2,ξ3),(1-(ξ1-η1+η2),ξ1-ξ2+ξ3-η1,η3))*det

    sum += f((1-ξ1,ξ1-ξ2+ξ3-η1,η1),(1-(ξ1-ξ3+η2),ξ1-ξ2,η2-η3))*det

    sum += f((1-(ξ1-ξ3+η1),ξ1-ξ2,η1),(1-ξ1,ξ1-ξ2+ξ3-η2,η3))*det

    sum += f((1-(ξ1-η1),ξ1-ξ2,ξ3-η1),(1-ξ1,ξ1-ξ2+η1-η2,η2-η3))*det

    sum += f((1-(ξ1-η2),ξ1-ξ2,η1-η2),(1-ξ1,ξ1-ξ2+η2-η3,ξ3-η2+η3))*det

    sum += f((1-(ξ1-η2),ξ1-ξ2,ξ3-η2),(1-ξ1,ξ1-ξ2+η3,η1-η3))*det

    sum += f((1-ξ1,ξ1-ξ2+η1,ξ3-η1),(1-(ξ1-η3),ξ1-ξ2,η2-η3))*det

    sum += f((1-ξ1,ξ1-ξ2+η2,η1-η2),(1-(ξ1-η3),ξ1-ξ2,ξ3-η3))*det

    sum += f((1-ξ1,ξ1-ξ2+η2,ξ3-η2),(1-(ξ1-η2+η3),ξ1-ξ2,η1-η2+η3))*det

    sum += f((1-(ξ1-η2+η3),ξ1-ξ2+η3,η1-η2),(1-ξ1,ξ1-ξ2,ξ3))*det

    sum += f((1-(ξ1-η3),ξ1-ξ2+η2-η3,ξ3-η2),(1-ξ1,ξ1-ξ2,η1))*det

    return sum
end

function k6p_cf_A(f, x1, x2, x3, y1, y2, z1)
    #sum = 0;

    #Domain 1
    #=1=#   ξ1 = x1
    #=2=#   ξ2 = x2
    #=3=#   ξ3 = x3
    #=4=#   η1 = x3 * y1
    #=5=#   η2 = x3 * y2
    #=6=#   η3 = x3 * z1

    det = x3^3

    sum = f((1-ξ1,ξ1-ξ2,ξ3-η3),(1-(ξ1-ξ3+η1),ξ1-ξ2+η2,η1-η2))*det

    return sum
end

function k6p_cf_B(f, x1, x2, x3, y1, y2, z1)
    #sum = 0;

    #Domain 3/5/13
    #=1=#   ξ1 = x1
    #=2=#   ξ2 = x2
    #=3=#   ξ3 = x3
    #=4=#   η1 = x3 * y1
    #=5=#   η2 = x3 * y2
    #=6=#   η3 = x3 * y1 * z1

    det = x3^3*y1

    sum = f((1-ξ1,ξ1-ξ2+ξ3-η1,η1-η2),(1-(ξ1-ξ3+η3),ξ1-ξ2,η3))*det

    sum += f((1-(ξ1-ξ3+η1),ξ1-ξ2,η1-η3),(1-ξ1,ξ1-ξ2+ξ3-η2,η2))*det

    sum += f((1-(ξ1-η2),ξ1-ξ2+η1-η2,ξ3-η1),(1-ξ1,ξ1-ξ2,η3))*det

    return sum
end

function k6p_cf_C(f, x1, x2, x3, y1, y2, y3)
    #sum = 0;
    #Domain 2/4/6-12/14-15
    #=1=#   ξ1 = x1
    #=2=#   ξ2 = x2
    #=3=#   ξ3 = x3
    #=4=#   η1 = x3 * y1
    #=5=#   η2 = x3 * y2
    #=6=#   η3 = x3 * y3

    det = x3^3

    sum = f((1-ξ1,ξ1-ξ2,ξ3),(1-(ξ1-η1+η2),ξ1-ξ2+ξ3-η1,η3))*det

    sum += f((1-ξ1,ξ1-ξ2+ξ3-η1,η1),(1-(ξ1-ξ3+η2),ξ1-ξ2,η2-η3))*det

    sum += f((1-(ξ1-ξ3+η1),ξ1-ξ2,η1),(1-ξ1,ξ1-ξ2+ξ3-η2,η3))*det

    sum += f((1-(ξ1-η1),ξ1-ξ2,ξ3-η1),(1-ξ1,ξ1-ξ2+η1-η2,η2-η3))*det

    sum += f((1-(ξ1-η2),ξ1-ξ2,η1-η2),(1-ξ1,ξ1-ξ2+η2-η3,ξ3-η2+η3))*det

    sum += f((1-(ξ1-η2),ξ1-ξ2,ξ3-η2),(1-ξ1,ξ1-ξ2+η3,η1-η3))*det

    sum += f((1-ξ1,ξ1-ξ2+η1,ξ3-η1),(1-(ξ1-η3),ξ1-ξ2,η2-η3))*det

    sum += f((1-ξ1,ξ1-ξ2+η2,η1-η2),(1-(ξ1-η3),ξ1-ξ2,ξ3-η3))*det

    sum += f((1-ξ1,ξ1-ξ2+η2,ξ3-η2),(1-(ξ1-η2+η3),ξ1-ξ2,η1-η2+η3))*det

    sum += f((1-(ξ1-η2+η3),ξ1-ξ2+η3,η1-η2),(1-ξ1,ξ1-ξ2,ξ3))*det

    sum += f((1-(ξ1-η3),ξ1-ξ2+η2-η3,ξ3-η2),(1-ξ1,ξ1-ξ2,η1))*det

    return sum
end

#overlapping edge
function k6p_ce(f, x1, x2, x3, y1, y2, y3)

    #sum = 0;
    #=1=#   ξ1 = x1
    #=2=#   ξ2 = x1 * x2
    #=3=#   ξ3 = x1 * x2 * x3
    #=4=#   η1 = x1 * x2 * x3 *y1
    #=5=#   η2 = x1 * x2 * y2
    #=6=#   η3 = x1 * x2 * y2 * y3

    det = x1^5*x2^4*x3*y2
    #Domain 1
    sum = f((1-ξ1,ξ1-ξ2+η1,ξ3-η1),(1-(ξ1+η2-ξ2),ξ1-ξ2,η2-η3))*det

    #=1=#   ξ1 = x1
    #=2=#   ξ2 = x1 * x2
    #=3=#   ξ3 = x1 * x2 * x3
    #=4=#   η1 = x1 * x2 * x3 * y1
    #=5=#   η2 = x1 * x2 * x3 * y1 * y2
    #=6=#   η3 = x1 * x2 * y3

    det = x1^5*x2^4*x3^2*y1

    #Domain 2
    sum += f((1-ξ1,ξ1-ξ2,ξ2-η3),(1-(ξ1-ξ2+ξ3),ξ1-ξ2+ξ3-η1,η2))*det

    #=1=#   ξ1 = x1
    #=2=#   ξ2 = x1 * x2
    #=3=#   ξ3 = x1 * x2 * x3
    #=4=#   η1 = x1 * x2 * x3 * y1
    #=5=#   η2 = x1 * x2 * x3 * y1 * y2
    #=6=#   η3 = x1 * x2 * x3 * y1 * y2 * y3

    det = x1^5*x2^4*x3^3*y1^2*y2

    #Domain 3
    sum += f((1-(ξ1-η1),ξ1-ξ2,ξ3-η1),(1-ξ1,ξ1-η2,η3))*det

    #=1=#   ξ1 = x1
    #=2=#   ξ2 = x1 * x2
    #=3=#   ξ3 = x1 * x2 * x3
    #=4=#   η1 = x1 * x2 * y1
    #=5=#   η2 = x1 * x2 * y1 * y2
    #=6=#   η3 = x1 * x2 * y1 * y2 * y3

    det = x1^5*x2^4*y1^2*y2

    #Domain 4
    sum += f((1-(ξ1-η3),ξ1-η1,η2-η3),(1-ξ1,ξ1-ξ2,ξ3))*det

    #=1=#   ξ1 = x1
    #=2=#   ξ2 = x1 * x2
    #=3=#   ξ3 = x1 * x2 * x3
    #=4=#   η1 = x1 * x2 * y1
    #=5=#   η2 = x1 * x2 * x3 * y1 * y2
    #=6=#   η3 = x1 * x2 * y1 * y3

    det = x1^5*x2^4*x3*y1^2

    #Domain 5
    sum += f((1-(ξ1-η2),ξ1-ξ2,ξ2-ξ3),(1-ξ1,ξ1-η1,η3))*det

    return sum
end

function k6p_ce_A(f, x1, x2, y1, y2, z1, z2)

    #sum = 0;
    #=1=#   ξ1 = x1
    #=2=#   ξ2 = x2
    #=3=#   ξ3 = x2 * y1
    #=4=#   η1 = x2 * y2
    #=5=#   η2 = x2 * z1
    #=6=#   η3 = x2 * z2

    det = x2^4
    #Domain 1
    sum = f((1-ξ1,ξ1-ξ2+η1,ξ3-η1),(1-(ξ1+η2-ξ2),ξ1-ξ2,η2-η3))*det

    return sum
end

function k6p_ce_B(f, x1, x2, y1, y2, y3, z1)

    #sum = 0;
    #=1=#   ξ1 = x1
    #=2=#   ξ2 = x2
    #=3=#   ξ3 = x2 * y1
    #=4=#   η1 = x2 * y2
    #=5=#   η2 = x2 * y3
    #=6=#   η3 = x2 * z1

    det = x2^4

    #Domain 2
    sum = f((1-ξ1,ξ1-ξ2,ξ2-η3),(1-(ξ1-ξ2+ξ3),ξ1-ξ2+ξ3-η1,η2))*det

    return sum
end

function k6p_ce_C(f, x1, x2, y1, y2, y3, y4)

    #sum = 0;
    #=1=#   ξ1 = x1
    #=2=#   ξ2 = x2
    #=3=#   ξ3 = x2 * y1
    #=4=#   η1 = x2 * y2
    #=5=#   η2 = x2 * y3
    #=6=#   η3 = x2 * y4

    det = x2^4

    #Domain 3
    sum = f((1-(ξ1-η1),ξ1-ξ2,ξ3-η1),(1-ξ1,ξ1-η2,η3))*det
    return sum
end

function k6p_ce_D(f, x1, x2, y1, z1, z2, z3)

    #sum = 0;
    #=1=#   ξ1 = x1
    #=2=#   ξ2 = x2
    #=3=#   ξ3 = x2 * y1
    #=4=#   η1 = x2 * z1
    #=5=#   η2 = x2 * z2
    #=6=#   η3 = x2 * z3

    det = x2^4

    #Domain 4
    sum = f((1-(ξ1-η3),ξ1-η1,η2-η3),(1-ξ1,ξ1-ξ2,ξ3))*det
    return sum
end

function k6p_ce_E(f, x1, x2, y1, z1, z2, w1)

    #sum = 0;
    #=1=#   ξ1 = x1
    #=2=#   ξ2 = x2
    #=3=#   ξ3 = x2 * y1
    #=4=#   η1 = x2 * z1
    #=5=#   η2 = x2 * y1 * z2
    #=6=#   η3 = x2 * z1 * w1

    det = x2^4*y1*z1

    #Domain 5
    sum = f((1-(ξ1-η2),ξ1-ξ2,ξ2-ξ3),(1-ξ1,ξ1-η1,η3))*det

    return sum
end

#overlapping vertex
function k6p_cv(f, η1, η2, η3, η4, ξ1, ξ2)

    ξ1η1   = ξ1 * η1
    ξ1η1η2 = ξ1 * η1 * η2

    ξ1ξ2     = ξ1 * ξ2
    ξ1ξ2η3   = ξ1 * ξ2 * η3
    ξ1ξ2η3η4 = ξ1 * ξ2 * η3 * η4

    ((ξ1^5) * η1 * (ξ2^2) * η3) *(

    f((1-ξ1,ξ1η1-ξ1η1η2,ξ1η1η2),(1-ξ1ξ2,ξ1ξ2η3-ξ1ξ2η3η4,ξ1ξ2η3η4))+

    f((1-ξ1ξ2,ξ1ξ2η3-ξ1ξ2η3η4,ξ1ξ2η3η4),(1-ξ1,ξ1η1-ξ1η1η2,ξ1η1η2))

    )

end

function k6p_cv_A(f, x1, x2, x3, y1, y2, y3)
    #sum = 0;

    sum = f((1-x1,x1-x2,x3),(1-x1*y1,x1*y1-x1*y2,x1*y3))*(x1^3)

    sum += f((1-x1*y1,x1*y1-x1*y2,x1*y3),(1-x1,x1-x2,x3))*(x1^3)

    return sum
end

#positive distance
function k6p_pd(f, x1, x2, x3, y1, y2, y3)

    #sum = 0;
    #=1=#   ξ1 = x1
    #=2=#   ξ2 = x1 * x2
    #=3=#   ξ3 = x1 * x2 * x3
    #=4=#   η1 = y1
    #=5=#   η2 = y1 * y2
    #=6=#   η3 = y1 * y2 * y3

    det = x1^2*x2*y1^2*y2
    #Domain 1
    sum = f((1-ξ1,ξ1-ξ2,ξ3),(1-η1,η1-η2,η3))*det

    return sum
end

function k6p_pd_A(f, x1, x2, x3, y1, y2, y3)

    #sum = 0;
    #Domain 1
    sum = f((1-x1,x1-x2,x3),(1-y1,y1-y2,y3))

    return sum
end

"""
integrals on  tetrahedron and triagle
"""

### Common Face ###
function k5p_cf(f, x1, x2, x3, y1, y2)
    #sum = 0;

    #Domain 1
    #=1=#   ξ1 = x1
    #=2=#   ξ2 = x1 * x2
    #=3=#   ξ3 = x1 * x2 * x3
    #=4=#   ξ4 = x1 * x2 * x3 *y1
    #=5=#   ξ5 = x1 * x2 * x3 *y2

    det = x1^4*x2^3*x3^2

    sum = f((1-ξ1,ξ1-ξ2,ξ3-ξ4),(1-(ξ1-ξ3+ξ5),ξ2-ξ3))*det

    #Domain 2-9
    #=1=#   ξ1 = x1
    #=2=#   ξ2 = x1 * x2
    #=3=#   ξ3 = x1 * x2 * x3
    #=4=#   ξ4 = x1 * x2 * x3 *y1
    #=5=#   ξ5 = x1 * x2 * x3 *y1*y2

    det = x1^4*x2^3*x3^2*y1

    sum += f((1-ξ1,ξ1-ξ2,ξ3),(1-(ξ1-ξ5),ξ2-ξ3+ξ4-ξ5))*det

    sum += f((1-ξ1,ξ1-ξ2+ξ4-ξ5,ξ3-ξ4),(1-(ξ1-ξ3),ξ2-ξ3))*det

    sum += f((1-ξ1,ξ1-ξ2+ξ3-ξ4,ξ4),(1-(ξ1-ξ3+ξ5),ξ2-ξ3+ξ5))*det

    sum += f((1-(ξ1-ξ5),ξ1-ξ2,ξ3-ξ4),(1-ξ1,ξ2-ξ3))*det

    sum += f((1-(ξ1-ξ4+ξ5),ξ1-ξ2,ξ3-ξ4+ξ5),(1-ξ1,ξ2-ξ3+ξ5))*det

    sum += f((1-(ξ1-ξ4),ξ1-ξ2,ξ3-ξ4),(1-ξ1,ξ2-ξ4+ξ5))*det

    sum += f((1-ξ1,ξ1-ξ2+ξ3-ξ4+ξ5,ξ4-ξ5),(1-(ξ1-ξ3+ξ4),ξ2-ξ3+ξ4))*det

    sum += f((1-(ξ1-ξ4+ξ5),ξ1-ξ2+ξ5,ξ3-ξ4),(1-ξ1,ξ2))*det

    return sum
end

function k5p_cf_A(f, x1, x2, x3, y1, z1)
    #sum = 0;

    #Domain 1
    #=1=#   ξ1 = x1
    #=2=#   ξ2 = x2
    #=3=#   ξ3 = x3
    #=4=#   ξ4 = x3 * y1
    #=5=#   ξ5 = x3 * z1

    det = x3^2

    sum = f((1-ξ1,ξ1-ξ2,ξ3-ξ4),(1-(ξ1-ξ3+ξ5),ξ2-ξ3))*det

    return sum
end

function k5p_cf_B(f, x1, x2, x3, y1, y2)
    #sum = 0;

    #Domain 2-9
    #=1=#   ξ1 = x1
    #=2=#   ξ2 = x2
    #=3=#   ξ3 = x3
    #=4=#   ξ4 = x3 * y1
    #=5=#   ξ5 = x3 * y2

    det = x3^2

    sum = f((1-ξ1,ξ1-ξ2,ξ3),(1-(ξ1-ξ5),ξ2-ξ3+ξ4-ξ5))*det

    sum += f((1-ξ1,ξ1-ξ2+ξ4-ξ5,ξ3-ξ4),(1-(ξ1-ξ3),ξ2-ξ3))*det

    sum += f((1-ξ1,ξ1-ξ2+ξ3-ξ4,ξ4),(1-(ξ1-ξ3+ξ5),ξ2-ξ3+ξ5))*det

    sum += f((1-(ξ1-ξ5),ξ1-ξ2,ξ3-ξ4),(1-ξ1,ξ2-ξ3))*det

    sum += f((1-(ξ1-ξ4+ξ5),ξ1-ξ2,ξ3-ξ4+ξ5),(1-ξ1,ξ2-ξ3+ξ5))*det

    sum += f((1-(ξ1-ξ4),ξ1-ξ2,ξ3-ξ4),(1-ξ1,ξ2-ξ4+ξ5))*det

    sum += f((1-ξ1,ξ1-ξ2+ξ3-ξ4+ξ5,ξ4-ξ5),(1-(ξ1-ξ3+ξ4),ξ2-ξ3+ξ4))*det

    sum += f((1-(ξ1-ξ4+ξ5),ξ1-ξ2+ξ5,ξ3-ξ4),(1-ξ1,ξ2))*det

    return sum
end

### Common Edge ###
function k5p_ce(f, x1, x2, x3, y1, y2)
    #sum = 0;

    #Domain 1/2
    #=1=#   ξ1 = x1
    #=2=#   ξ2 = x1 * x2
    #=3=#   ξ3 = x1 * x2 * x3
    #=4=#   ξ4 = x1 * x2 * x3 *y1
    #=5=#   ξ5 = x1 * x2 * y2

    det = x1^4*x2^3*x3

    sum = f((1-ξ1,ξ1-ξ3,ξ4),(1-(ξ1-ξ5),ξ2-ξ5))*det

    sum += f((1-ξ1,ξ1-ξ2,ξ5),(1-(ξ1-ξ4),ξ3-ξ4))*det

    #Domain 3/4
    #=1=#   ξ1 = x1
    #=2=#   ξ2 = x1 * x2
    #=3=#   ξ3 = x1 * x2 * x3
    #=4=#   ξ4 = x1 * x2 * x3 *y1
    #=5=#   ξ5 = x1 * x2 * x3 *y1* y2

    det = x1^4*x2^3*x3^2*y1

    sum += f((1-(ξ1-ξ4),ξ1-ξ2,ξ3-ξ4),(1-ξ1,ξ5))*det

    sum += f((1-(ξ1-ξ5),ξ1-ξ3,ξ4-ξ5),(1-ξ1,ξ2))*det

    #Domain 5
    #=1=#   ξ1 = x1
    #=2=#   ξ2 = x1 * x2
    #=3=#   ξ3 = x1 * x2 * x3
    #=4=#   ξ4 = x1 * x2 * x3 *y1
    #=5=#   ξ5 = x1 * x2 * x3 *y2

    det = x1^4*x2^3*x3^2

    sum += f((1-(ξ1-ξ2+ξ3),ξ1-ξ2,ξ4),(1-ξ1,ξ2-ξ3+ξ5))*det

    return sum

end

function k5p_ce_A(f, x1, x2, y1, y2, z1)
    #sum = 0;

    #Domain 1/2
    #=1=#   ξ1 = x1
    #=2=#   ξ2 = x2
    #=3=#   ξ3 = x2 * y1
    #=4=#   ξ4 = x2 * y2
    #=5=#   ξ5 = x2 * z1

    det = x2^3

    sum = f((1-ξ1,ξ1-ξ3,ξ4),(1-(ξ1-ξ5),ξ2-ξ5))*det

    sum += f((1-ξ1,ξ1-ξ2,ξ5),(1-(ξ1-ξ4),ξ3-ξ4))*det

    return sum
end

function k5p_ce_B(f, x1, x2, y1, y2, y3)
    #sum = 0;

    #Domain 3/4
    #=1=#   ξ1 = x1
    #=2=#   ξ2 = x2
    #=3=#   ξ3 = x2 * y1
    #=4=#   ξ4 = x2 * y2
    #=5=#   ξ5 = x2 * y3

    det = x2^3

    sum = f((1-(ξ1-ξ4),ξ1-ξ2,ξ3-ξ4),(1-ξ1,ξ5))*det

    sum += f((1-(ξ1-ξ5),ξ1-ξ3,ξ4-ξ5),(1-ξ1,ξ2))*det

    return sum
end

function k5p_ce_C(f, x1, x2, y1, y2, z1)
    #sum = 0;

    #Domain 5
    #=1=#   ξ1 = x1
    #=2=#   ξ2 = x2
    #=3=#   ξ3 = x2 * y1
    #=4=#   ξ4 = x2 * y2
    #=5=#   ξ5 = x2 * y1 * z1

    det = x2^3*y1

    sum = f((1-(ξ1-ξ2+ξ3),ξ1-ξ2,ξ4),(1-ξ1,ξ2-ξ3+ξ5))*det

    return sum
end

### Common Vertex ###
function k5p_cv(f, x1, x2, x3, x4, x5)
    #sum = 0;

    sum = f((1-x1,x1-x1*x2,x1*x2*x3),(1-x1*x4,x1*x4*x5))*(x1^4*x2*x4)

    sum += f((1-x1*x3,x1*x3-x1*x3*x4,x1*x3*x4*x5),(1-x1,x1*x2))*(x1^4*x3^2*x4)

    return sum
end

function k5p_cv_A(f, x1, x2, x3, x4, x5)
    #sum = 0;

    sum = f((1-x1,x1-x2,x3),(1-x1*x4,x1*x5))*(x1^2)

    sum += f((1-x4*x1,x4*x1-x4*x2,x4*x3),(1-x4,x5))*(x4^3)

    return sum
end

### Positive Distance ###
function k5p_pd(f, x1, x2, x3, y1, y2)

    #u = [1-x1, x1-x1*x2, x1*x2*x3] #jacobian x1^2*x2
    #v = [1-y1, y1*y2] #jacobian y1

    return f((1-x1, x1-x1*x2, x1*x2*x3),(1-y1, y1*y2))*(x1^2*x2)*(y1)
end

function k5p_pd_A(f, x1, x2, x3, y1, y2)

    #u = [1-x1, x1-x2, x3]
    #v = [1-y1, y2]

    return f((1-x1, x1-x2, x3),(1-y1, y2))
end

"""
Integrals on triangles
"""
### Common Face ###
function k4p_cf(f, η1, η2, η3, ξ)

	(ξ^3)*((η1)^2)*(η2)*(

    f((1-ξ, ξ-ξ*η1+ξ*η1*η2), (1-(ξ-ξ*η1*η2*η3), ξ-ξ*η1)) +

	f((1-(ξ-ξ*η1*η2*η3), ξ-ξ*η1), (1-ξ, ξ-ξ*η1+ξ*η1*η2))  +

	f((1-ξ, ξ*η1*(1-η2+η2*η3)), (1-(ξ-ξ*η1*η2), ξ*η1*(1-η2))) +

	f((1-(ξ-ξ*η1*η2), ξ*η1*(1-η2)), (1-ξ, ξ*η1*(1-η2+η2*η3))) +

	f((1-(ξ-ξ*η1*η2*η3), ξ*η1*(1-η2*η3)), (1-ξ, ξ*η1*(1-η2))) +

	f((1-ξ, ξ*η1*(1-η2)), (1-(ξ-ξ*η1*η2*η3), ξ*η1*(1-η2*η3)))

    )
end

function k4p_cf_A(f, x1, x2, x3, y1)

	#sum = 0;
    #Domain 1-6
    #=1=#   ξ1 = x1
    #=2=#   ξ2 = x2
    #=3=#   ξ3 = x3
    #=4=#   ξ4 = x3 * y1

    det = x3

    sum=f((1-ξ1,ξ1-ξ2+ξ3),(1-(ξ1-ξ4),ξ1-ξ2))*det

    sum+=f((1-(ξ1-ξ4),ξ1-ξ2),(1-ξ1,ξ1-ξ2+ξ3))*det

    sum+=f((1-ξ1,ξ2-ξ3+ξ4),(1-(ξ1-ξ3),ξ2-ξ3))*det

    sum+=f((1-(ξ1-ξ3),ξ2-ξ3),(1-ξ1,ξ2-ξ3+ξ4))*det

    sum+=f((1-(ξ1-ξ4),ξ2-ξ4),(1-ξ1,ξ2-ξ3))*det

    sum+=f((1-ξ1,ξ2-ξ3),(1-(ξ1-ξ4),ξ2-ξ4))*det

    return sum
end

### Common edge ###
function k4p_ce(f, η1, η2, η3, ξ)

	ξη1 = ξ * η1
	η1η2 = η1 * η2
	η2η3 = η2 * η3
	η1η2η3 = η1η2 * η3

	(ξ^3)*((η1)^2) *f((1-ξ, ξη1*η3), (1-ξ*(1-η1η2), ξη1*(1-η2))) +

	(ξ^3)*((η1)^2)*(η2)*(
	f((1-ξ, ξη1), (1-ξ*(1-η1η2η3), ξη1*η2*(1-η3))) +

	f((1-ξ*(1-η1η2), ξη1*(1-η2)), (1-ξ, ξη1*η2η3)) +

	f((1-ξ*(1-η1η2η3), ξη1*η2*(1-η3)), (1-ξ, ξη1)) +

	f((1-ξ*(1-η1η2η3), ξη1*(1-η2η3)), (1-ξ, ξη1*η2))
	)
end

function k4p_ce_A(f, x1, x2, y1, z1)
   # sum = 0;

    #Domain 1
    #=1=#   ξ1 = x1
    #=2=#   ξ2 = x2
    #=3=#   ξ3 = x2*y1
    #=4=#   ξ4 = x2 *z1

    det = x2^2

    sum=f((1-ξ1,ξ4),(1-(ξ1-ξ3),ξ2-ξ3))*det

    return sum
end

function k4p_ce_B(f, x1, x2, y1, y2)
    #sum = 0;
    #Domain 2-5
    #=1=#   ξ1 = x1
    #=2=#   ξ2 = x2
    #=3=#   ξ3 = x2 * y1
    #=4=#   ξ4 = x2 * y2

    det = x2^2

    sum=f((1-ξ1,ξ2),(1-(ξ1-ξ4),ξ3-ξ4))*det

    sum+=f((1-(ξ1-ξ3),(ξ2-ξ3)),(1-ξ1,ξ4))*det

    sum+=f((1-(ξ1-ξ4),(ξ3-ξ4)),(1-ξ1,ξ2))*det

    sum+=f((1-(ξ1-ξ4),(ξ2-ξ4)),(1-ξ1,ξ3))*det

    return sum
end

### Common vertex ###
function k4p_cv(f,η1, η2, η3, ξ)

	ξη1 = ξ * η1
	ξη2 = ξ * η2

	(ξ^3)*η2*(

	f((1-ξ, ξη1), (1-ξη2, ξη2*η3)) +  f((1-ξη2, ξη2*η3), (1-ξ, ξη1))

	)
end

function k4p_cv_A(f,x1, x2, y1, y2)

	(x1^2)*(

	f((1-x1, x2), (1-x1*y1, x1*y2)) +  f((1-x1*y1, x1*y2), (1-x1, x2))

	)
end

### Positive Distance ###
function k4p_pd(f,x1, x2, y1, y2)
	(x1*y1)*f((1-x1, x1*x2), (1-y1, y1*y2))
end

function k4p_pd_A(f,x1, x2, y1, y2)
	f((1-x1, x2), (1-y1, y2))
end

#=
function k3p_cf_orig2(f, η1, η2, η3, ξ)

	(ξ^3)*((η1)^2)*(η2)*(

    f((1-ξ, ξ*η1), (1-(ξ-ξ*η1*η2*η3), ξ*η1-ξ*η1*η2))  +

	f((1-(ξ-ξ*η1*η2*η3), ξ*η1-ξ*η1*η2), (1-ξ, ξ*η1)) +

	f((1-ξ, ξ*η1*(1-η2+η2*η3)), (1-(ξ-ξ*η1*η2), ξ*η1*(1-η2))) +

	f((1-(ξ-ξ*η1*η2), ξ*η1*(1-η2)), (1-ξ, ξ*η1*(1-η2+η2*η3))) +

	f((1-(ξ-ξ*η1*η2*η3), ξ*η1*(1-η2*η3)), (1-ξ, ξ*η1*(1-η2))) +

	f((1-ξ, ξ*η1*(1-η2)), (1-(ξ-ξ*η1*η2*η3), ξ*η1*(1-η2*η3)))

	)
end

function k3p_cf_new(f, ξ1, ξ2, η1, η2)

	(ξ1^3)*(ξ2^2)*(η1)*(

    #OK
    f((1-ξ1, ξ1*ξ2), (1-(ξ1-ξ1*ξ2*η1+ξ1*ξ2*η1*η2), ξ1*ξ2-ξ1*ξ2*η1))  +
    #OK
    f((1-ξ1, ξ1*ξ2-ξ1*ξ2*η1+ξ1*ξ2*η1*η2), (1-(ξ1-ξ1*ξ2*η1), ξ1*ξ2-ξ1*ξ2*η1)) +
    #OK
    f((1-ξ1, ξ1*ξ2-ξ1*ξ2*η1), (1-(ξ1-ξ1*ξ2*η1*η2), ξ1*ξ2-ξ1*ξ2*η1*η2)) +
    #OK
    f((1-(ξ1-ξ1*ξ2*η1*η2), ξ1*ξ2-ξ1*ξ2*η1*η2),(1-ξ1, ξ1*ξ2-ξ1*ξ2*η1)) +
    #OK
    f((1-(ξ1-ξ1*ξ2*η1), ξ1*ξ2-ξ1*ξ2*η1),(1-ξ1, ξ1*ξ2-ξ1*ξ2*η1+ξ1*ξ2*η1*η2)) +
	#OK
    f( (1-(ξ1-ξ1*ξ2*η1+ξ1*ξ2*η1*η2), ξ1*ξ2-ξ1*ξ2*η1),(1-ξ1, ξ1*ξ2))
	)
end

function k3p_cf_gm(f, ξ1, ξ2, η1, η2)

	(η1)*(

    #OK
    f((1-ξ1, ξ2), (1-(ξ1-η1+η1*η2), ξ2-η1))  +
    #OK
    f((1-ξ1, ξ2-η1+η1*η2), (1-(ξ1-η1), ξ2-η1)) +
    #OK
    f((1-ξ1, ξ2-η1), (1-(ξ1-η1*η2), ξ2-η1*η2)) +
    #OK
    f((1-(ξ1-η1*η2), ξ2-η1*η2),(1-ξ1, ξ2-η1)) +
    #OK
    f((1-(ξ1-η1), ξ2-η1),(1-ξ1, ξ2-η1+η1*η2)) +
	#OK
    f( (1-(ξ1-η1+η1*η2), ξ2-η1),(1-ξ1, ξ2))
	)
end
=#
#=
function k3p_ce_new(f, ξ1, ξ2, η1, η2)

    #OK
    (ξ1^3)*(ξ2^2) * f((1-ξ1, ξ1*ξ2-ξ1*ξ2*η2), (1-(ξ1*ξ2*η1-ξ1*ξ2+ξ1), ξ1*ξ2*η1))  +

	(ξ1^3)*(ξ2^2)*(η1)*(
    #OK
	f((1-ξ1, ξ1*ξ2*η1-ξ1*ξ2+ξ1), (1-ξ1*ξ2, ξ1*ξ2*η1-ξ1*ξ2*η1*η2)) +
    #OK
    f((1-(ξ1-ξ1*ξ2*η1), ξ1*ξ2-ξ1*ξ2*η1), (1-ξ1, ξ1*ξ2*η1*η2)) +
    #OK
	f((1-(ξ1-ξ1*ξ2*η1*η2), ξ1*ξ2*η1-ξ1*ξ2*η1*η2), (1-ξ1, ξ1*ξ2)) +
    #OK
	f((1-(ξ1-ξ1*ξ2*η1+ξ1*ξ2*η1*η2), ξ1*ξ2*η1*η2-ξ1*ξ2*η1+ξ1*ξ2), (1-ξ1, ξ1*ξ2*η1))

	)
end

function k3p_ce_gm(f, ξ1, ξ2, η1, η2)

    #OK
    (ξ2^2) * f((1-ξ1, ξ2-ξ2*η2), (1-(ξ2*η1-ξ2+ξ1), ξ2*η1))  +

	(ξ2^2)*(η1)*(
    #OK
	f((1-ξ1, ξ2*η1-ξ2+ξ1), (1-ξ2, ξ2*η1-ξ2*η1*η2)) +
    #OK
    f((1-(ξ1-ξ2*η1), ξ2-ξ2*η1), (1-ξ1, ξ2*η1*η2)) +
    #OK
	f((1-(ξ1-ξ2*η1*η2), ξ2*η1-ξ2*η1*η2), (1-ξ1, ξ2)) +
    #OK
	f((1-(ξ1-ξ2*η1+ξ2*η1*η2), ξ2*η1*η2-ξ2*η1+ξ2), (1-ξ1, ξ2*η1))

	)
end
=#

#=
function k3p_cv_gm(f,ξ1, ξ2, η1, η2)

    (ξ1^2)*(

	f((1-ξ1, ξ2), (1-ξ1*η1, ξ1*η2)) +  f((1-ξ1*η1, ξ1*η2), (1-ξ1, ξ2))

    )

    #=
	(ξ1^2)*η1*(

	f((1-ξ1, ξ2), (1-ξ1*η1, ξ1*η1*η2)) +  f((1-ξ1*η1, ξ1*η1*η2), (1-ξ1, ξ2))

    )
    =#
end
=#


text:
using LinearAlgebra
using ProgressMeter
using StatsFuns

"""
Binomial Likelihood Bayesian Graph Fused Elastic Net Model
"""
mutable struct BayesianBinomialGFEN
    # trail info
    edges::Vector{Tuple{Int, Int}}

    num_nodes::Int
    num_edges::Int
    num_trails::Int

    # model parameters and computed quantities
    nbrs::Vector{Vector{Int}}
    tv1::Vector{Vector{Float64}}
    tv2::Vector{Vector{Float64}}
    ridge::Vector{Float64}
    lasso::Vector{Float64}

    # convergence and steps
    trained::Bool
    steps::Int

    # main constructor
    function BayesianBinomialGFEN(
        edges::Vector{Tuple{Int, Int}};
        tv1::Union{Vector{Float64},Float64} = 0.0,
        tv2::Union{Vector{Float64},Float64} = 0.0,
        lasso::Union{Vector{Float64},Float64} = 0.0,
        ridge::Union{Vector{Float64},Float64} = 0.0
    )
        # new object
        x = new()

        # assign trails
        x.edges = edges
        x.num_edges = length(edges)
        x.num_nodes = maximum(maximum(s) for s in edges)

        # assign penalties
        tv1 = fill_if_scalar(tv1, x.num_edges)
        tv2 = fill_if_scalar(tv2, x.num_edges)
        lasso = fill_if_scalar(lasso, x.num_nodes)
        ridge = fill_if_scalar(ridge, x.num_nodes)

        @assert length(tv1) == x.num_edges
        @assert length(tv2) == x.num_edges
        @assert length(lasso) == x.num_nodes
        @assert length(ridge) == x.num_nodes

        x.lasso = lasso
        x.ridge = ridge

        # make neighbors map
        x.nbrs = [Int[] for _ in 1:x.num_nodes]
        x.tv1 = [Float64[] for _ in 1:x.num_nodes]
        x.tv2 = [Float64[] for _ in 1:x.num_nodes]

        for (i, (s, t)) in enumerate(edges)
            push!(x.nbrs[s], t)
            push!(x.nbrs[t], s)

            push!(x.tv1[s], tv1[i])
            push!(x.tv1[t], tv1[i])

            push!(x.tv2[s], tv2[i])
            push!(x.tv2[t], tv2[i])
        end

        for i in 1:x.num_nodes
            @assert length(x.nbrs[i]) == length(x.tv1[i])
            @assert length(x.nbrs[i]) == length(x.tv2[i])
        end

        # convergence and steps
        x.trained = false
        x.steps = 0

        # end constructor
        return x
    end
end
Base.show(io::IO, x::BayesianBinomialGFEN) = print(io, "BayesianBinomialGFEN")
Base.print(io::IO, x::BayesianBinomialGFEN) = print(io, "BayesianBinomialGFEN")

function binomial_gibbs_step(
    s::Float64,  # successes
    a::Float64,  # attempts
    nbr_values::Vector{Float64},  # values of neighbors
    tv1::Vector{Float64},  # edge penalties l1
    tv2::Vector{Float64},  # edge penalties l2
    lasso::Float64,  # l1 reg
    ridge::Float64;  # l2 reg
    clamp_value::Float64 = 10.0
)::Float64
    # add a small constant for numerical stability
    n_nbrs =  length(nbr_values)
    # heuristic_scale = clamp(2^(1.0 / a), 1.0, 1e6)

    # define target loglikelihood, any target is ok as long as it is concave
    # and it is easy to provide one point with positive and negative slope
    function target_dens(θ)::Tuple{Float64, Float64}
        if θ >= 0.0
            z = exp(-θ)
            ω = z / (1.0 + z)
            logll = - a * log(1.0 + z) - (a - s) * θ
            ∇logll = a * ω - (a - s)
        else
            z = exp(θ)
            ω = z / (1.0 + z)
            logll = s * θ - a * log(1.0 + z)
            ∇logll = s - a * ω
        end

        tv1_reg, tv2_reg, ∇tv1_reg, ∇tv2_reg = 0., 0., 0., 0.
        if n_nbrs > 0
            δ_nbrs = θ .- nbr_values
            tv1_reg = dot(abs.(δ_nbrs), tv1)
            tv2_reg = 0.5 * dot(δ_nbrs.^2, tv2)
            ∇tv1_reg = dot(sign.(δ_nbrs), tv1)
            ∇tv2_reg = dot(δ_nbrs, tv2)
        end

        lasso_reg = lasso * abs(θ)
        ∇lasso_reg = lasso * sign(θ)
        ridge_reg = 0.5ridge * θ^2
        ∇ridge_reg = ridge * θ
        f = (logll - tv1_reg - tv2_reg - lasso_reg - ridge_reg) # / heuristic_scale
        g = (∇logll - ∇tv1_reg - ∇tv2_reg - ∇lasso_reg - ∇ridge_reg)  #/ heuristic_scale
        f, g
    end

    # it's easy to define points with positive and negative slopes
    # based on the minimum and max posible values
    # efficiency of the sample increase with good envelope points
    lower = max(min(logit(s / a),  minimum(nbr_values) - 1e-6), -clamp_value)
    upper = min(max(logit(s / a), maximum(nbr_values) + 1e-6), clamp_value)
    # lower = 0.0
    # upper = 0.0
    # build sampler
    support = (-clamp_value, clamp_value)
    init = (lower, upper)
    sampler = RejectionSampler(
        target_dens, support, init, autograd=false, use_secants=true, apply_log=false
    )

    # return one sample
    out = run_sampler!(sampler, 1)[1]
    num_attempts = length(sampler.envelop)
    # println("Num attempts: $num_attempts")
    # error("hi")
    return out
end

function binomial_gibbs_sweep(
    θ::Vector{Float64},  # current
    s::Vector{Float64},  # successes
    a::Vector{Float64},  # attempts
    nbrs::Vector{Vector{Int}}, # neighbors
    tv1::Vector{Vector{Float64}},
    tv2::Vector{Vector{Float64}},
    lasso::Vector{Float64},
    ridge::Vector{Float64};
    async::Bool=true
)
    T = length(θ)
    θnew = zeros(Float64, T)

    # async updates using multiple threads
    if async
        @inbounds @threads for t in 1:T
            nbr_values = [θ[i] for i in nbrs[t]]
            θnew[t] = binomial_gibbs_step(
                s[t], a[t], nbr_values, tv1[t], tv2[t], lasso[t], ridge[t]
            )
        end
    else
        @inbounds for t in 1:T
            nbr_values = [θ[i] for i in nbrs[t]]
            θnew[t] = binomial_gibbs_step(
                s[t], a[t], nbr_values, tv1[t], tv2[t], lasso[t], ridge[t]
            )
        end
    end
    return θnew
end

function sample_chain(
    m::BayesianBinomialGFEN,  # model
    s::Vector{Float64},  # successes
    a::Vector{Float64},  # attempts
    n::Int;  # chain iterations (full sweeps),
    init::Union{Vector{Float64},Nothing} = nothing,  # initial values for chain,
    init_eps::Float64=0.0,
    async::Bool,
    thinning::Int = 1,
    verbose::Bool = false,
    burnin::Int = 0
)
    # store each gibbs sweep in matrix
    T = length(s)
    θ = zeros(Float64, T, (n - burnin) ÷ thinning)

    num_saved = 0
    if isnothing(init)
        init = logistic.((s .+ init_eps) ./ (a .+ 2init_eps))
    end
    θ_curr = init

    verbose && (pbar = Progress(n))
    for i in 1:n
        θ_curr = binomial_gibbs_sweep(
            θ_curr, s, a, m.nbrs, m.tv1, m.tv2, m.lasso, m.ridge, async=async
        )
        if i > burnin && i % thinning == 0
            θ[:, num_saved + 1] = θ_curr
            num_saved += 1
        end
        verbose && (next!(pbar))
    end

    return θ
end