text:
module Touchstone

import Base: ==, ≈

export Options, DataPoint, NoiseDataPoint, TouchstoneData
export freqs, params, param, mags, dBmags, angs, reals, imags
export noiseFreqs, minNoiseFigures, reflCoeffs, effNoiseRess
export version, ports, refs
export parse_touchstone_stream, parse_touchstone_string, parse_touchstone_file
export write_touchstone_stream, write_touchstone_string, write_touchstone_file
export cmp

"Holds values from option line."
struct Options
  unit::Float64
  parameter::Symbol
  format::Symbol
  resistance::Float64
end
Options() = Options( units[ default_units ], parameters[ default_parameter ], formats[ default_format ], default_resistance )
"Compares Options."
function ==( o1, o2 :: Options )
  o1.unit == o2.unit &&
  o1.parameter == o2.parameter &&
  o1.format == o2.format &&
  o1.resistance == o2.resistance
end
function ≈( o1, o2 :: Options )
  o1.unit ≈ o2.unit &&
  o1.parameter ≈ o2.parameter &&
  o1.format ≈ o2.format &&
  o1.resistance ≈ o2.resistance
end
function cmp( o1, o2::Options )
  if o1 != o2
    println( "Options differs:" )
    println( "o1: " )
    println( o1 )
    println( "o2: " )
    println( o2 )
  end
end

"Holds scalar frequency and parameter matrix for one data point."
struct DataPoint
  frequency::Float64
  parameter::Matrix{Complex{Float64}}
end
"Compares data points."
function ==( dp1, dp2 :: DataPoint )
  dp1.frequency == dp2.frequency &&
  dp1.parameter == dp2.parameter
end
"Compares data points approximately."
function ≈( dp1, dp2 :: DataPoint )
  dp1.frequency ≈ dp2.frequency &&
  dp1.parameter ≈ dp2.parameter
end
function cmp( dp1, dp2 :: DataPoint )
  if dp1 != dp2
    println( "DataPoint differs:" )
    println( "dp1: " )
    println( dp1 )
    println( "dp2: " )
    println( dp2 )
    cmp
  end
end
function cmp( dp1, dp2 :: Vector{ T } ) where T
  for pair in zip( dp1, dp2 )
    cmp( pair... )
  end
end
struct NoiseDataPoint
  frequency::Float64
  minNoiseFigure::Float64
  reflCoeff::Complex{Float64}
  effNoiseRes::Float64
end
function ==( dp1, dp2 :: NoiseDataPoint )
  dp1.frequency == dp2.frequency &&
  dp1.minNoiseFigure == dp2.minNoiseFigure &&
  dp1.reflCoeff == dp2.reflCoeff &&
  dp1.effNoiseRes == dp2.effNoiseRes
end
function ≈( dp1, dp2 :: NoiseDataPoint )
  dp1.frequency ≈ dp2.frequency &&
  dp1.minNoiseFigure ≈ dp2.minNoiseFigure &&
  dp1.reflCoeff ≈ dp2.reflCoeff &&
  dp1.effNoiseRes ≈ dp2.effNoiseRes
end
function cmp( dp1, dp2 :: NoiseDataPoint )
  if dp1 != dp2
    println( "NoiseDataPoint differs:" )
    println( "dp1: " )
    println( dp1 )
    println( "dp2: " )
    println( dp2 )
  end
end
"Holds data for Touchstone file."
struct TouchstoneData
  data::Vector{ DataPoint }
  noiseData::Vector{ NoiseDataPoint }
  options::Options
  comments::Vector{ String }
  keywordparams::Dict{ Symbol, Any }
  function TouchstoneData(
    data,
    noiseData = Vector{ NoiseDataPoint }(),
    options = Options(),
    comments = Vector{ String }(),
    keywordparams  = Dict{ Symbol, Any }()
  )
    if !issorted( data, by=x -> x.frequency )
      error( "Frequencies of data vector not in ascending order." )
    end
    if !issorted( noiseData, by=x -> x.frequency )
      error( "Frequencies of noise data vector not in ascending order." )
    end
    if length( data ) > 0 && length( noiseData ) > 0 && data[ end ].frequency < noiseData[ 1 ].frequency
      error( "Last frequency of data vector is smaller than first frequency of noise data vector." )
    end
    v = version( keywordparams )
    if v == "1.0"
      p = ports( data )
    elseif v == "2.0"
      p = ports( keywordparams )
      if haskey( keywordparams, :NumberOfFrequencies )
        n = keywordparams[ :NumberOfFrequencies ]
        if length( data ) != n
          error( "V2.0: $( n ) data rows expected, $( length( data ) ) found." )
        end
      else
        error( "V2.0: [Number of Frequencies] keyword expected." )
      end
      if haskey( keywordparams, :Reference ) && length( keywordparams[ :Reference ] ) != p
        error( "V2.0: $( p ) parameters for [Reference] expected." )
      end
      if haskey( keywordparams, :MixedModeOrder ) && options.parameter in [ :HybridHParameters, :HybridGParameters ]
        error( "[Mixed-Mode Order] keyword for $( parameterstrings[ options.parameter ] ) parameter format not allowed." )
      end
      if p == 2
        if haskey( keywordparams, :TwoPortDataOrder )
          dataOrder = keywordparams[ :TwoPortDataOrder ]
          if dataOrder != "12_21"  && dataOrder != "21_12"
            error( "V2.0: Unknown parameter '$dataOrder' for [Two-Port Data Order] keyword." )
          end
        else
          error( "V2.0: [Two-Port Data Order] expected for two port data." )
        end
      else
        if haskey( keywordparams, :NoiseData )
          error( "V2.0: [NoiseData] keyword not allowed for $( p ) port data." )
        end
        if haskey( keywordparams, :TwoPortDataOrder )
          error( "V2.0: [Two-Port Data Order] keyword not allowed for $( p ) port data." )
        end
      end
    end
    if p != 2
      if options.parameter in [ :HybridHParameters, :HybridGParameters ]
        error( "$( parameterstrings[ options.parameter ] ) parameter format for $( p )-port files not allowed." )
      end
    end
    new( data, noiseData, options, comments, keywordparams )
  end
end
function cmp( ts1, ts2 :: TouchstoneData )
  if ts1 != ts2
    println( "TouchstoneData differs:" )
    println( "ts1: " )
    println( ts1 )
    println( "ts2: " )
    println( ts2 )
    cmp( ts1.data, ts2.data )
    cmp( ts1.noiseData, ts2.noiseData )
    cmp( ts1.options, ts2.options )

  end
end
"Compares Touchstone data."
function ==( ts1, ts2 :: TouchstoneData )
    ts1.data == ts2.data &&
    ts1.noiseData == ts2.noiseData &&
    ts1.options == ts2.options &&
    ts1.comments == ts2.comments &&
    ts1.keywordparams == ts2.keywordparams
end
function ≈( ts1, ts2 :: TouchstoneData )
    ts1.data == ts2.data || ts1.data ≈ ts2.data &&
    ts1.noiseData == ts2.noiseData || ts1.noiseData ≈ ts2.noiseData &&
    ts1.options == ts2.options || ts1.options ≈ ts2.options &&
    ts1.comments == ts2.comments &&
    ts1.keywordparams == ts2.keywordparams
end
"""
    freqs( ts )

Gets the frequency vector from Touchstone data.
"""
freqs( dp::DataPoint ) = dp.frequency
freqs( ts::TouchstoneData ) = map( freqs, ts.data )

noiseFreqs( ndp::NoiseDataPoint ) = ndp.frequency
noiseFreqs( ts::TouchstoneData ) = map( noiseFreqs, ts.noiseData )
minNoiseFigures( ndp::NoiseDataPoint ) = ndp.minNoiseFigure
minNoiseFigures( ts::TouchstoneData ) = map( minNoiseFigures, ts.noiseData )
reflCoeffs( ndp::NoiseDataPoint ) = ndp.reflCoeff
reflCoeffs( ts::TouchstoneData ) = map( reflCoeffs, ts.noiseData )
effNoiseRess( ndp::NoiseDataPoint ) = ndp.effNoiseRes
effNoiseRess( ts::TouchstoneData ) = map( effNoiseRess, ts.noiseData )

"""
    params( ts )

Gets the vector of parameter matrices from Touchstone data.
"""
params( dp::DataPoint ) =  dp.parameter
params( ts::TouchstoneData ) = map( params, ts.data )

"""
    param( ts, p1, p2 )

Gets the vector of parameters with indices p1, p2 from Touchstone data.
"""
param( dp::DataPoint, p1 = 1, p2 = 1 ) =  dp.parameter[ p1, p2 ]
param( ts::TouchstoneData, p1 = 1, p2 = 1 ) = map( dp -> param( dp, p1, p2 ), ts.data )

"""
    mags( ts, p1, p2 )

Gets the vector of the magnitudes of parameters with indices p1, p2 from Touchstone data.
"""
mags( dp::DataPoint, p1 = 1, p2 = 1 ) = abs( param( dp, p1, p2 ) )
mags( ts::TouchstoneData, p1 = 1, p2 = 1 ) = map( dp -> mags( dp, p1, p2 ), ts.data )

"""
    dBmags( ts, p1, p2 )

Gets the vector of the magnitudes in dB of parameters with indices p1, p2 from Touchstone data.
"""
dBmags( dp::DataPoint, p1 = 1, p2 = 1 ) = 20log10( mags( dp, p1, p2 ) )
dBmags( ts::TouchstoneData, p1 = 1, p2 = 1 ) = map( dp -> dBmags( dp, p1, p2 ), ts.data )

"""
    angs( ts, p1, p2 )

Gets the vector of the angles in radians of parameters with indices p1, p2 from Touchstone data.
"""
angs( dp::DataPoint, p1 = 1, p2 = 1 ) = rad2deg( angle( param( dp, p1, p2 ) ) )
angs( ts::TouchstoneData, p1 = 1, p2 = 1 ) = map( dp -> angs( dp, p1, p2 ), ts.data )

"""
    reals( ts, p1, p2 )

Gets the vector of the real parts of parameters with indices p1, p2 from Touchstone data.
"""
reals( dp::DataPoint, p1 = 1, p2 = 1 ) = real( param( dp, p1, p2 ) )
reals( ts::TouchstoneData, p1 = 1, p2 = 1 ) = map( dp -> reals( dp, p1, p2 ), ts.data )

"""
    imags( ts, p1, p2 )

Gets the vector of the imaginary parts of parameters with indices p1, p2 from Touchstone data.
"""
imags( dp::DataPoint, p1 = 1, p2 = 1 ) = imag( param( dp, p1, p2 ) )
imags( ts::TouchstoneData, p1 = 1, p2 = 1 ) = map( dp -> imags( dp, p1, p2 ), ts.data )

function version( keywordparams::Dict{ Symbol, Any } )
  if haskey( keywordparams, :Version )
    if keywordparams[ :Version ] == 2
      return "2.0"
    else
      error( "Unknown version." )
    end
  else
    return "1.0"
  end
end
version( ts::TouchstoneData ) = version( ts.keywordparams )

#V2.0
ports( keywordparams::Dict{ Symbol, Any } ) = keywordparams[ :NumberOfPorts ]
#V1.0
ports( data::Vector{ DataPoint } ) = length( data ) > 0 ? size( data[ 1 ].parameter, 1 ) : 0
#V1.0, V2.0
ports( ts::TouchstoneData ) = version( ts ) == "2.0" ? ports( ts.keywordparams ) : ports( ts.data )

function refs( ts::TouchstoneData )
  if version( ts ) == "2.0" && haskey( ts.keywordparams, :Reference )
    refs = ts.keywordparams[ :Reference ]
  else
    refs = fill( ts.options.resistance, ports( ts ) )
  end
  return refs
end

include( "ParseTouchstone.jl" )
include( "WriteTouchstone.jl" )

end


text:
import OrdinaryDiffEq, DiffEqCallbacks

"""
integrate()
"""
function integrate(T::Vector, f:: Function, state0::S, fout::Function;
                    alg::OrdinaryDiffEq.OrdinaryDiffEqAlgorithm = OrdinaryDiffEq.DP5(),
                    callback = nothing, kwargs...) where S

    if isa(state0, Vector{<:Real})
        x0 = state0
    else
        x0 = state0.data
    end

    N = state0.N
    fout_diff(u, t, integrator) = fout(t, S(N, deepcopy(u)))

    out_type = pure_inference(fout, Tuple{eltype(T),typeof(state0)})
    out = DiffEqCallbacks.SavedValues(eltype(T),out_type)
    scb = DiffEqCallbacks.SavingCallback(fout_diff,out,saveat=T,
                                        save_everystep=false,
                                        save_start = false)

    prob = OrdinaryDiffEq.ODEProblem(f, x0, (T[1], T[end]))

    full_cb = OrdinaryDiffEq.CallbackSet(callback, scb)

    sol = OrdinaryDiffEq.solve(prob, alg;
            reltol=1.0e-6,
            abstol=1.0e-8,
            save_everystep = false,
            save_start = false,
            save_end = false,
            callback=full_cb,
            kwargs...)

    out.t, out.saveval

end

Base.@pure pure_inference(fout,T) = Core.Compiler.return_type(fout, T)

ext: jl


text:
using Logging: with_logger
using Test: @test, @testset, @test_throws

using GitHubActions
using SimpleMock: called_once_with, mock
using Suppressor: @capture_out

const GHA = GitHubActions

@testset "GitHubActions.jl" begin
    @test (@capture_out GHA.command("a", (), "")) == "::a::\n"
    @test (@capture_out GHA.command("a", (), nothing)) == "::a::\n"
    @test (@capture_out GHA.command("a", (), "bar")) == "::a::bar\n"
    @test (@capture_out GHA.command("a", (), ())) == "::a::[]\n"
    @test (@capture_out GHA.command("a", (b=1,), "")) == "::a b=1::\n"
    @test (@capture_out GHA.command("a", (b=1, c=2), "")) == "::a b=1,c=2::\n"
    @test (@capture_out GHA.command("a", (b="c%d\re\nf",), "")) == "::a b=c%25d%0De%0Af::\n"
    @test (@capture_out GHA.command("a", (), "a%b\rc\nd:e,f")) == "::a::a%25b%0Dc%0Ad:e,f\n"

    @test (@capture_out end_group()) == "::endgroup::\n"

    withenv(() -> (@test get_state("a") == ""), "STATE_a" => nothing)
    withenv(() -> (@test get_state("a") == ""), "STATE_a" => "")
    withenv(() -> (@test get_state("a") == "b"), "STATE_a" => "b")

    @test (@capture_out log_debug("a")) == "::debug::a\n"
    @test (@capture_out log_error("a")) == "::error::a\n"
    @test (@capture_out log_warning("a")) == "::warning::a\n"

    @test (@capture_out save_state("a", "b")) == "::save-state name=a::b\n"

    @test (@capture_out set_command_echo(true)) == "::echo::on\n"
    @test (@capture_out set_command_echo(false)) == "::echo::off\n"

    @test (@capture_out set_output("a", "b")) == "::set-output name=a::b\n"

    @test (@capture_out set_secret("a")) == "::add-mask::a\n"

    @test (@capture_out start_group("a")) == "::group::a\n"

    mktemp() do file, io
        withenv("GITHUB_PATH" => file, "PATH" => "/bin") do
            sep = Sys.iswindows() ? ';' : ':'
            add_path("a")
            @test ENV["PATH"] == string("a", sep, "/bin")
            @test read(file, String) == "a\n"
            add_path("b")
            @test ENV["PATH"] == string("b", sep, "a", sep, "/bin")
            @test read(file, String) == "a\nb\n"
        end
    end

    withenv(() -> (@test get_input("A") == ""), "INPUT_A" => "")
    withenv(() -> (@test get_input("A") == ""), "INPUT_A" => nothing)
    withenv(() -> (@test get_input("A") == "b"), "INPUT_A" => "b")
    withenv(() -> (@test get_input("a") == "b"), "INPUT_A" => "b")
    withenv(() -> (@test get_input("a b") == "c"), "INPUT_A_B" => "c")
    withenv(() -> (@test_throws GHA.MissingInputError get_input("a"; required=true)), "INPUT_A" => "")
    withenv(() -> (@test_throws GHA.MissingInputError get_input("a"; required=true)), "INPUT_A" => nothing)

    @test (@capture_out group(() -> println("!"), "a")) == "::group::a\n!\n::endgroup::\n"

    mktemp() do file, io
        withenv("GITHUB_ENV" => file, map(c -> string(c) => nothing, 'a':'z')...) do
            set_env("a", "b")
            @test ENV["a"] == "b"
            @test read(file, String) == "a<<EOF\nb\nEOF\n"
            set_env("b", "fooEOFbar")
            @test read(file, String) == "a<<EOF\nb\nEOF\nb<<EOFEOF\nfooEOFbar\nEOFEOF\n"
            rm(file)
            set_env("c", [])
            @test ENV["c"] == "[]"
            @test read(file, String) == "c<<EOF\n[]\nEOF\n"
            set_env("d", nothing)
            @test ENV["d"] == ""
            @test read(file, String) == "c<<EOF\n[]\nEOF\nd<<EOF\n\nEOF\n"
        end
    end

    if VERSION.minor < 6
        mock(atexit) do ae
            @test (@capture_out set_failed("a")) == "::error::a\n"
            @test called_once_with(ae, GHA.fail)
        end
    end

    rx(level) = Regex("^::$level file=$(@__FILE__),line=\\d+::a")
    with_logger(GitHubActionsLogger()) do
        @test match(rx("debug"), (@capture_out @debug "a")) !== nothing
        @test match(rx("warning"), (@capture_out @warn "a")) !== nothing
        @test match(rx("error"), (@capture_out @error "a")) !== nothing
        @test (@capture_out @info "a") == "a\n"

        @test (@capture_out @info "a" b=1 c=2 d=Text("e\nf")) == "a\n  b = 1\n  c = 2\n  d = \n    e\n    f\n"
        @test endswith((@capture_out @warn "a" b=1 c=2), "::a%0A  b = 1%0A  c = 2\n")

        expected = "::warning file=bar,line=1::foo\n"
        @test (@capture_out @warn "foo" location=("bar", 1)) == expected
    end
end


text:
##############################################################################
##
## @eachrow
##
##############################################################################

# Recursive function that traverses the syntax tree of e, replaces instances of
# ":(:(x))" with ":x[row]".
eachrow_replace(x) = x
eachrow_replace(e::QuoteNode) = Expr(:ref, e, :row)

function eachrow_replace(e::Expr)
    if onearg(e, :^)
        return e.args[2]
    end

    # Traverse the syntax tree of e
    col = get_column_expr(e)
    if col !== nothing
        return :($e[row])
    # equivalent to protect_replace_syms
    elseif e.head == :.
        x_new = eachrow_replace(e.args[1])
        y = e.args[2]
        y_new = y isa Expr ? eachrow_replace(y) : y

        return Expr(:., x_new, y_new)
    else
        mapexpr(eachrow_replace, e)
    end
end

protect_eachrow_replace(e) = e
protect_eachrow_replace(e::Expr) = eachrow_replace(e)

function eachrow_replace_dotted(e, membernames)
    x_new = eachrow_repalce(e.args[1])
    y_new = protect_eachrow_repalce(e.args[2])
    Expr(:., x_new, y_new)
end

function eachrow_find_newcols(e::Expr, newcol_decl)
    if e.head == :macrocall && e.args[1] == Symbol("@newcol")
        ea = e.args[3]
        # expression to assign a new column to df
        return (nothing, Any[Expr(:(=), ea.args[1], Expr(:call, ea.args[2], :undef, :_N))])
    else
        if isempty(e.args)
            return (e.args, Any[])
        end
        newargs = Any[]
        for ea in e.args
            (nea, newcol) = eachrow_find_newcols(ea, newcol_decl)
            nea != nothing && push!(newargs, nea)
            nea == nothing && length(newcol) > 0 && append!(newcol_decl, newcol)
        end
        return (Expr(e.head, newargs...), newcol_decl)
    end
end

eachrow_find_newcols(x, newcol_decl) = (x, Any[])

function eachrow_helper(df, body, deprecation_warning)
    # @deprecate cannot be used because eachrow is a macro, and the @warn should not be in
    # eachrow itself because then it will be displayed when the macro is evaluated.
    deprecation_warning && @warn "`@byrow!` and `@byrow` are deprecated, use `@eachrow` instead."
    e_body, e_newcols = eachrow_find_newcols(body, Any[])
    _df = gensym()
    quote
        let $_df = $df
            local _N = nrow($_df)
            local _DF = @transform($_df, $(e_newcols...))
            $(with_helper(:_DF, :(for row = 1:_N
                $(eachrow_replace(e_body))
            end)))
            _DF
        end
    end
end

"""
    @eachrow(df, body)

Act on each row of a data frame, producing a new dataframe.
Similar to

```
for row in eachrow(copy(df))
    ...
end
```

Includes support for control flow and `begin end` blocks. Since the
"environment" induced by `@eachrow df` is implicitly a single row of `df`,
use regular operators and comparisons instead of their elementwise counterparts
as in `@with`. Note that the scope within `@eachrow` is a hard scope.

`eachrow` also supports special syntax for allocating new columns. The syntax
`@newcol x::Vector{Int}` allocates a new uninitialized column `:x` with an `Vector` container
with eltype `Int`.This feature makes it easier to use `eachrow` for data
transformations. `_N` is introduced to represent the number of rows in the data frame,
`_DF` represents the `DataFrame` including added columns, and `row` represents
the index of the current row.

Changes to the rows do not affect `df` but instead a freshly allocated data frame is returned
by `@eachrow`. Also note that the returned data frame does not share columns
with `df`. See [`@eachrow!`](@ref) which employs the same syntax but modifies
the data frame in-place.

Like with `@transform`, `@eachrow` supports the use of `$DOLLAR` to work with column names
stored as variables. Using `$DOLLAR` with a multi-column selector, such as a `Vector` of
`Symbol`s, is currently unsupported.

### Arguments

* `df` : an `AbstractDataFrame`
* `expr` : expression operated on row by row

### Returns

The modified `AbstractDataFrame`.

### Examples

```julia
julia> using DataFramesMeta

julia> df = DataFrame(A = 1:3, B = [2, 1, 2]);

julia> let x = 0
            @eachrow df begin
                if :A + :B == 3
                    x += 1
                end
            end  #  This doesn't work without the let
            x
        end
2

julia> @eachrow df begin
            if :A > :B
                :A = 0
            end
        end
3×2 DataFrame
 Row │ A      B
     │ Int64  Int64
─────┼──────────────
   1 │     1      2
   2 │     0      1
   3 │     0      2

julia> df2 = @eachrow df begin
           @newcol :colX::Vector{Float64}
           :colX = :B == 2 ? pi * :A : :B
       end
3×3 DataFrame
 Row │ A      B      colX
     │ Int64  Int64  Float64
─────┼───────────────────────
   1 │     1      2  3.14159
   2 │     2      1  1.0
   3 │     3      2  9.42478

julia> varA = :A; varB = :B;

julia> df2 = @eachrow df begin
           @newcol :colX::Vector{Float64}
           :colX = $(DOLLAR)varB == 2 ? pi * $(DOLLAR)varA : $(DOLLAR)varB
       end
3×3 DataFrame
 Row │ A      B      colX
     │ Int64  Int64  Float64
─────┼───────────────────────
   1 │     1      2  3.14159
   2 │     2      1  1.0
   3 │     3      2  9.42478

julia> x = [1, 1, 1];

julia> @eachrow df begin
           x[row] = :A
       end;

julia> x
3-element Vector{Int64}:
 1
 2
 3

julia> @eachrow df begin
           @newcol :m::Vector{Float64}
           :m = mean(_DF[:, row])
       end
3×3 DataFrame
 Row │ A      B      m
     │ Int64  Int64  Float64
─────┼───────────────────────
   1 │     1      2  2.0
   2 │     2      1  1.66667
   3 │     3      2  1.22222

```
"""
macro eachrow(df, body)
    esc(eachrow_helper(df, body, false))
end

function eachrow!_helper(df, body)
    e_body, e_newcols = eachrow_find_newcols(body, Any[])
    _df = gensym()
    quote
        let $_df = $df
            local _N = nrow($_df)
            local _DF = @transform!($_df, $(e_newcols...))
            $(with_helper(:_DF, :(for row = 1:_N
                $(eachrow_replace(e_body))
            end)))
            _DF
        end
    end
end

"""
    @eachrow!(df, body)

Act on each row of a data frame in-place, similar to

```
for row in eachrow(df)
    ... # Actions that modify `df`.
end
```

Includes support for control flow and `begin end` blocks. Since the
"environment" induced by `@eachrow! df` is implicitly a single row of `df`,
use regular operators and comparisons instead of their elementwise counterparts
as in `@with`. Note that the scope within `@eachrow!` is a hard scope.

`eachrow!` also supports special syntax for allocating new columns. The syntax
`@newcol x::Vector{Int}` allocates a new uninitialized column `:x` with an `Vector` container
with eltype `Int`.This feature makes it easier to use `eachrow` for data
transformations. `_N` is introduced to represent the number of rows in the data frame,
`_DF` represents the `dataframe` including added columns, and `row` represents
the index of the current row.

Changes to the rows directly affect `df`. The operation will modify the
data frame in place. See [`@eachrow`](@ref) which employs the same syntax but allocates
a fresh data frame.

Like with `@transform!`, `@eachrow!` supports the use of `$DOLLAR` to work with column names
stored as variables. Using `$DOLLAR` with a multi-column selector, such as a `Vector` of
`Symbol`s, is currently unsupported.

### Arguments

* `df` : an `AbstractDataFrame`
* `expr` : expression operated on row by row

### Returns

The modified `AbstractDataFrame`.

### Examples

```julia
julia> using DataFramesMeta

julia> df = DataFrame(A = 1:3, B = [2, 1, 2]);

julia> let x = 0
            @eachrow! df begin
                if :A + :B == 3
                    x += 1
                end
            end  #  This doesn't work without the let
            x
        end
2

julia> df2 = copy(df);

julia> @eachrow! df2 begin
           if :A > :B
               :A = 0
           end
       end;

julia> df2
3×2 DataFrame
 Row │ A      B
     │ Int64  Int64
─────┼──────────────
   1 │     1      2
   2 │     0      1
   3 │     0      2

julia> df2 = copy(df);

julia> @eachrow! df2 begin
           @newcol :colX::Vector{Float64}
           :colX = :B == 2 ? pi * :A : :B
       end
3×3 DataFrame
 Row │ A      B      colX
     │ Int64  Int64  Float64
─────┼───────────────────────
   1 │     1      2  3.14159
   2 │     2      1  1.0
   3 │     3      2  9.42478

julia> varA = :A; varB = :B;

julia> df2 = copy(df);

julia> @eachrow! df2 begin
           @newcol :colX::Vector{Float64}
           :colX = $(DOLLAR)varB == 2 ? pi * $(DOLLAR)varA : $(DOLLAR)varB
       end
3×3 DataFrame
 Row │ A      B      colX
     │ Int64  Int64  Float64
─────┼───────────────────────
   1 │     1      2  3.14159
   2 │     2      1  1.0
   3 │     3      2  9.42478

julia> x = [1, 1, 1];

julia> @eachrow! df begin
           x[row] = :A
       end;

julia> x
3-element Vector{Int64}:
 1
 2
 3

julia> @eachrow! df begin
           @newcol :m::Vector{Float64}
           :m = mean(_DF[:, row])
       end
3×3 DataFrame
 Row │ A      B      m
     │ Int64  Int64  Float64
─────┼───────────────────────
   1 │     1      2  2.0
   2 │     2      1  1.66667
   3 │     3      2  1.22222
```
"""
macro eachrow!(df, body)
    esc(eachrow!_helper(df, body))
end


text:
function get_edges(obj::EdgesWithNodeInfoBuilder)
    return jcall(obj, "getEdges", List, ())
end

function visit_internal_node(obj::EdgesWithNodeInfoBuilder, arg0::BSPTree)
    return jcall(obj, "visitInternalNode", void, (BSPTree,), arg0)
end

function visit_leaf_node(obj::EdgesWithNodeInfoBuilder, arg0::BSPTree)
    return jcall(obj, "visitLeafNode", void, (BSPTree,), arg0)
end

function visit_order(obj::EdgesWithNodeInfoBuilder, arg0::BSPTree)
    return jcall(obj, "visitOrder", BSPTreeVisitor_Order, (BSPTree,), arg0)
end


text:
using Documenter, DynamicHMC

makedocs(modules = [DynamicHMC],
         format = Documenter.HTML(; prettyurls = get(ENV, "CI", nothing) == "true"),
         clean = true,
         sitename = "DynamicHMC.jl",
         authors = "Tamás K. Papp",
         checkdocs = :exports,
         # strict = true,
         pages = ["Introduction" => "index.md",
                  "A worked example" => "worked_example.md",
                  "Documentation" => "interface.md"])

deploydocs(repo = "github.com/tpapp/DynamicHMC.jl.git",
           push_preview = true)


text:
# Conway's Game Of Life in a GtkCanvas

using Gtk, Images
using Profile, BenchmarkTools

gtkCanvas = @GtkCanvas()
win = GtkWindow(gtkCanvas, "CA", 1000, 1000, false, true)
show(gtkCanvas)

# call back for mouse events
gtkCanvas.mouse.button1press = @guarded (widget, event) -> begin
    ctx = getgc(widget)
    set_source_rgb(ctx, 0, 1, 0)
    x = event.x
    y = event.y

    arc(ctx, x, y, 5, 0, 2pi)
    rectangle(ctx, x - 15, y - 15, 30, 30)
    stroke(ctx)
    reveal(widget)
end

@guarded draw(gtkCanvas) do widget
    # Called by the main loop at specific frequency

    #   memory estimate:  16 bytes
    #   allocs estimate:  1
    #   --------------
    #   minimum time:     11.210 μs (0.00% GC)
    #   median time:      16.250 μs (0.00% GC)
    #   mean time:        374.226 μs (0.00% GC)
    #   maximum time:     64.180 ms (0.00% GC)
    #   --------------
    #   samples:          10000
    #   evals/sample:     1

    ctx = getgc(gtkCanvas)
    h = height(gtkCanvas)
    w = width(gtkCanvas)

    clearColor = RGB(0.14, 0.18, 0.22)
    set_source_rgb(ctx, clearColor.r, clearColor.g, clearColor.b) # 56/255, 60/255, 74/255) - - - ---
    rectangle(ctx, 0, 0, w, h)
    fill(ctx)

   # draw_grid(gtkCanvas, 5, 5)
   # stroke(ctx)

    draw_ca(ctx, h, w)
    fill(ctx)
end

function draw_grid(canvas, grid_width, grid_height)
    ctx = getgc(canvas)
    h_canvas = height(canvas)
    w_canvas = width(canvas)
    w_cell = round(w_canvas / grid_width)
    h_cell = round(h_canvas / grid_height)
    cell_size = min(w_cell, h_cell)
    w_grid = cell_size * grid_width
    h_grid = cell_size * grid_height

    set_source_rgb(ctx, 0, 0, 0);
    set_line_width(ctx, 0.5);

    for x in 0:cell_size:w_grid
        move_to(ctx, x, 0)
        line_to(ctx, x, h_grid)
    end

    set_source_rgb(ctx, 0.23, 0.25, 0.29)
    for y  in 0:cell_size:h_grid
        move_to(ctx, 0, y)
        line_to(ctx, h_grid, y)
    end
    stroke(ctx)
end

function draw_ca(ctx, h_canvas, w_canvas)
    # called within @guarded draw(gtkCanvas) above
    matrix_h, matrix_w = size(ca.matrix)
    w_cell = round(w_canvas / matrix_w)
    h_cell = round(h_canvas / matrix_h)
    color = RGB(0.6, 0.32, 0.32)
    set_source_rgb(ctx, color.r, color.g, color.b)
    for i in 0:matrix_h-1, j in 0:matrix_w-1
        value = ca.matrix[i + 1, j + 1]
        if value == 1
            rectangle(ctx, i * w_cell,  j * h_cell, w_cell - 1, h_cell - 1)
        end
    end
end

println("""=========================================================================================

                                ███▄ ▄███▓ ▄▄▄       ██▓ ███▄     █
                                ▓██▒▀█▀ ██▒▒████▄    ▓██▒ ██ ▀█   █
                                ▓██    ▓██░▒██  ▀█▄  ▒██▒▓██  ▀█ ██▒
                                ▒██    ▒██ ░██▄▄▄▄██ ░██░▓██▒  ▐▌██▒
                                ▒██▒   ░██▒ ▓█   ▓██▒░██░▒██░   ▓██░
                                ░ ▒░   ░  ░ ▒▒   ▓▒█░░▓  ░ ▒░   ▒ ▒
                                ░  ░      ░  ▒   ▒▒ ░ ▒ ░░ ░░   ░ ▒░
                                ░      ░     ░   ▒    ▒ ░   ░   ░ ░
                                    ░         ░  ░ ░           ░

===============================================================================================""")

# cd = Condition()
# endit(ww) = notify(cd)
# signal_connect(endit, win, :destroy)
# showall(win)
# wait(cd)

mutable struct CellularAutomata
    matrix
end

function gol_number_of_adjacent_cells_slow(i, j, previous_state_matrix)
    # elegant implementation but very slow :

    #   memory estimate:  42.81 KiB
    #   allocs estimate:  9
    #   --------------
    #   minimum time:     4.889 μs (0.00% GC)
    #   median time:      5.860 μs (0.00% GC)
    #   mean time:        7.435 μs (14.42% GC)
    #   maximum time:     374.442 μs (93.94% GC)
    #   --------------
    #   samples:          10000
    #   evals/sample:     7

    h,w = size(previous_state_matrix)

    empty_row = reshape( [ 0 for i = 1:w ] , (1,w) )
    empty_col = [ 0 for i = 1:h+2]

    tmp = [ empty_row ; previous_state_matrix ; empty_row ]
    tmp = [ empty_col  tmp  empty_col]

    sub = tmp[i:i+2,j:j+2]
    sub[2,2] = 0 # don't count the center cell
    sum(sub)
end

function gol_number_of_adjacent_cells(i, j, matrix)
    # faster implementation :

    #   memory estimate:  0 bytes
    #   allocs estimate:  0
    #   --------------
    #   minimum time:     34.048 ns (0.00% GC)
    #   median time:      35.196 ns (0.00% GC)
    #   mean time:        35.230 ns (0.00% GC)
    #   maximum time:     101.239 ns (0.00% GC)
    #   --------------
    #   samples:          10000
    #   evals/sample:     993

    h, w = size(matrix)

    sum = 0
    if i == 1
        if j == 1
            #  __..
            # |0x
            # |xx

            sum = matrix[ 2, 1 ] + matrix[ 1, 2 ] + matrix[ 2, 2 ]
        elseif j == w
            #         ..__
            #           xo|
            #           xx|

            sum = matrix[ i+1, j ] + matrix[ i,j-1 ] + matrix[ i+1, j-1 ]
        else
            #   ..___..
            #     xox
            #     xxx

            sum = matrix[ i+1, j-1 ] + matrix[ i+1, j ] + matrix[ i+1, j+1 ] +
                        matrix[ i, j-1 ] + matrix[ i, j+1 ]
        end
    elseif i == h
        if j == 1
            # :
            # |xx
            # |ox
            #  --

            sum = matrix[ i-1, j ] + matrix[ i-1, j+1 ] + matrix[ i, j+1 ]
        elseif j == w
            #               :
            #             xx|
            #             xo|
            #             --

            sum = matrix[ i, j-1 ] + matrix[ i-1, j ] + matrix[ i-1, j-1 ]
        else
            #
            #      xxx
            #      xox
            #    ..---..

            sum = matrix[ i-1, j-1 ] + matrix[ i-1, j ] + matrix[ i-1, j+1 ] +
                matrix[ i, j-1 ] + matrix[ i, j+1 ]
        end
    else # 1 < i < h
        if j == 1
            #  :
            #  |xx
            #  |ox
            #  |xx
            #  :

            sum = matrix[ i-1, j+1 ] + matrix[ i, j+1 ] + matrix[ i+1, j+1] +
                 matrix[ i-1, j ] + matrix[ i+1, j ]
        elseif j == w
            #          :
            #        xx|
            #        xo|
            #        xx|
            #          :

            sum = matrix[ i-1, j-1 ] + matrix[ i, j-1 ] + matrix[ i+1, j-1 ] +
                matrix[ i-1, j ] + matrix[ i+1, j ]
        else
            #     xxx
            #     xox
            #     xxx

            sum = matrix[ i-1, j-1 ] + matrix[ i-1, j ] + matrix[ i-1, j+1 ] +
                  matrix[ i,   j-1 ]                    + matrix[ i,   j+1 ] +
                  matrix[ i+1, j-1 ] + matrix[ i+1, j ] + matrix[ i+1, j+1 ]
        end
    end
    sum
end

# GOL rules
function gol_compute_next_state(i, j, previous_state_matrix)
	adjacent_live_cell_count = gol_number_of_adjacent_cells(i, j, previous_state_matrix)
    previous_state = previous_state_matrix[i,j]
    next_state = 0
    if previous_state == 1
        if adjacent_live_cell_count < 2
            next_state = 0
        elseif adjacent_live_cell_count <= 3
            next_state = 1
        elseif adjacent_live_cell_count > 3
            next_state = 0
        end
    elseif previous_state == 0
        if adjacent_live_cell_count == 3
            next_state = 1
        else
            next_state = 0
        end
    end
    next_state
end

function tick_ca!(ca::CellularAutomata)
    h, w = size(ca.matrix)
    next_matrix = similar(ca.matrix)
    for i=1:w, j=1:h
        next_matrix[i,j] = gol_compute_next_state(i, j, ca.matrix)
    end
    ca.matrix = next_matrix
end

function tick_ca!(ca::CellularAutomata, steps)
    for _ = 1:steps
        tick_ca!(ca)
    end
end

# init ca
target_frequency =10
s = 600
w, h = s, s
matrix = [ rand([0,1]) for i = 1:h, j = 1:w ]
ca = CellularAutomata(matrix)

function main_loop()
    t=time()

    tick_ca!(ca)

    tick_time = time()-t
    t = time()

    draw(gtkCanvas)

    draw_time = time() - t

    time_to_sleep = 1 / target_frequency - draw_time - tick_time
    if time_to_sleep > 0
        sleep(1 / target_frequency - draw_time)
    end
    println("* ca tick  time : $(round(tick_time, digits=9)) ")
    println("- grk draw time : $(round(draw_time, digits=9)) ")

    actual_frequency = 1/(time()-t)
    println("                             actual tick frequency : $actual_frequency")
end

while true
    main_loop()
end

max_stars_repo_name: leroy-walton/playgrounds

max_issues_repo_name: leroy-walton/playgrounds

max_forks_repo_name: leroy-walton/playgrounds


text:
"""
    RMSProp(::Vector{XVariable}; lr=1e-2, inertia=0.99, eps=1e-8, L1decay=0.001, L2decay=0.01)

Implements RMSProp algorithm.
"""
mutable struct RMSProp <: Optimizer
    xparams::Vector{XVariable}
    w::Vector
    lr::AbstractFloat
    ϵ::AbstractFloat
    inertia::AbstractFloat
    L1decay::AbstractFloat
    L2decay::AbstractFloat
    name::String
    function RMSProp(xparams::Vector{XVariable}; lr=1e-2, inertia=0.99, eps=1e-8, L1decay=0.001, L2decay=0.01)
        n = length(xparams)
        w = Vector(undef,n)
        for i = 1:n
            c , θ = xparams[i]
            w[i] = Zeros(typeof(ᵛ(θ)), θ.shape)
        end
        new(xparams, w, lr, eps, inertia, L1decay, L2decay, "RMSProp")
    end
end

function Base.show(io::IO, O::RMSProp)
    print("RMSProp(lr=$(O.lr), ϵ=$(O.ϵ), inertia=$(O.inertia), L1decay=$(O.L1decay), L2decay=$(O.L2decay))")
end

function update!(O::RMSProp; clipfn::Function=LPInfNormClip, clipvalue=10.0)
    w = O.w
    ϵ = O.ϵ
    μ = - O.lr
    ρ = O.inertia
    λ₁ = O.L1decay
    λ₂ = O.L2decay

    for i = 1:length(O.xparams)
        c , θ = O.xparams[i]
        ∇ = clipfn(setNanInfZero(δ(θ)), clipvalue)
        𝒗 = ᵛ(θ)
        @. w[i] = ρ * w[i] + (1-ρ) * ∇ * ∇
        if c == 'w'
            if λ₁==0 && λ₂==0
                @. 𝒗 += μ / (sqrt(w[i])+ϵ) * ∇
            elseif λ₁==0 && λ₂!=0
                @. 𝒗 += μ / (sqrt(w[i])+ϵ) * (∇ + λ₂ * 𝒗)
            elseif λ₁!=0 && λ₂==0
                @. 𝒗 += μ / (sqrt(w[i])+ϵ) * (∇ + λ₁ * sign(𝒗))
            else  # λ₁!=0 && λ₂!=0
                @. 𝒗 += μ / (sqrt(w[i])+ϵ) * (∇ + λ₁ * sign(𝒗) + λ₂ * 𝒗)
            end
        else
            @. 𝒗 += μ / (sqrt(w[i])+ϵ) * ∇
        end
    end
end

max_stars_repo_name: sonosole/Delta

max_issues_repo_name: sonosole/Delta

max_forks_repo_name: sonosole/Delta


text:
## EXPORTING LEARNING NETWORKS AS MODELS WITH @from_network

# closure to generate the fit methods for exported composite. Here `mach`
function fit_method(mach, models...)

    signature = mach.fitresult
    mach_args = mach.args

    function _fit(model::M, verbosity::Integer, args...) where M
        length(args) > length(mach_args) &&
            throw(ArgumentError("$M does not support more than "*
                                "$(length(mach_args)) training arguments"))
        replacement_models = [getproperty(model, fld)
                              for fld in fieldnames(M)]
        model_replacements = [models[j] => replacement_models[j]
                              for j in eachindex(models)]
        source_replacements = [mach_args[i] => source(args[i])
                               for i in eachindex(args)]
        replacements = vcat(model_replacements, source_replacements)

        new_mach =
            replace(mach, replacements...; empty_unspecified_sources=true)

        return!(new_mach, model, verbosity)
    end

    return _fit
end

net_error(message) = throw(ArgumentError("Learning network export error.\n"*
                                     string(message)))
net_error(k::Int) = throw(ArgumentError("Learning network export error $k. "))

_insert_subtyping(ex, subtype_ex) =
    Expr(:(<:), ex, subtype_ex)

_exported_type(::Probabilistic) = :ProbabilisticComposite
_exported_type(::Deterministic) = :DeterministicComposite
_exported_type(::Unsupervised) = :UnsupervisedComposite
_exported_type(::Static) = :StaticComposite

function eval_and_reassign(modl, ex)
    s = gensym()
    evaluated = modl.eval(ex)
    if evaluated isa Symbol
        hack = String(evaluated)
        modl.eval(:($s = Symbol($hack)))
    else
        modl.eval(:($s = $evaluated))
    end
    return s, evaluated
end

function without_line_numbers(block_ex)
    block_ex.head == :block || throw(ArgumentError)
    args = filter(block_ex.args) do arg
        !(arg isa LineNumberNode)
    end
    return Expr(:block, args...)
end

function from_network_preprocess(modl, mach_ex, block_ex)

    mach_ex, mach  = eval_and_reassign(modl, mach_ex)
    mach isa Machine{<:Surrogate} ||
        net_error("$mach is not a learning network machine. ")
    if block_ex.head == :block
        block_ex = without_line_numbers(block_ex)
        struct_ex = block_ex.args[1]
        trait_declaration_exs = block_ex.args[2:end]
    elseif block_ex.head == :struct
        struct_ex = block_ex
        trait_declaration_exs = []
    else
        net_error("Expected `struct`, `mutable struct` or "*
                  "`begin ... end` block, but got `$block_ex` ")
    end

    # if necessary add or modify struct subtyping:
    if struct_ex.args[2] isa Symbol
        struct_ex.args[2] = _insert_subtyping(struct_ex.args[2],
                                              _exported_type(mach.model))
        modeltype_ex = struct_ex.args[2].args[1]
    elseif struct_ex.args[2] isa Expr
        struct_ex.args[2].head == :(<:) ||
                    net_error("Badly formed `struct` subtying. ")
        modeltype_ex = struct_ex.args[2].args[1]
        super = eval(struct_ex.args[2].args[2])
        inferred_super_ex = _exported_type(mach.model)
        if !(super <: Composite)
            @warn "New composite type must subtype `Composite` but "*
            "`$super` does not. Instead declaring "*
            "`$modeltype_ex <: $inferred_super_ex`. "
            struct_ex.args[2].args[2] = inferred_super_ex
        end
    else
        net_error(41)
    end

    # test if there are no fields:
    field_exs = without_line_numbers(struct_ex.args[3]).args
    no_fields = isempty(field_exs)

    # extract trait definitions:
    trait_ex_given_name_ex = Dict{Symbol,Any}()

    ne() = net_error("Bad trait declaration. ")
    for ex in trait_declaration_exs
        ex isa Expr           || ne()
        ex.head == :(=)       || ne()
        ex.args[1] isa Symbol || ne()
        ex.args[1] in MLJModelInterface.MODEL_TRAITS ||
            net_error("Expected a model trait as keywork but "*
                      "got $(ex.args[2]). Options are:\n"*
                      "$MLJModelInterface.MODEL_TRAIES. ")
        length(ex.args) == 2  || ne()
        trait_ex_given_name_ex[ex.args[1]] = ex.args[2]
    end

    return mach_ex, modeltype_ex, struct_ex, no_fields, trait_ex_given_name_ex

end

function from_network_(modl,
                       mach_ex,
                       modeltype_ex,
                       struct_ex,
                       no_fields,
                       trait_ex_given_name_ex)

    args = gensym(:args)
    models = gensym(:models)
    instance = gensym(:instance)

    # Define the new model type with keyword constructor:
    if no_fields
        modl.eval(struct_ex)
    else
        modl.eval(MLJBase.Parameters.with_kw(struct_ex, modl, false))
    end

    # Test that an instance can be created:
    try
        modl.eval(:($modeltype_ex()))
    catch e
        @error "Problem instantiating a default instance of the "*
        "new composite type. Each field name in the struct expression "*
        "must have a corresponding model instance (that also appears "*
        "somewhere in the network). "*
        "Perhaps you forgot to specify one of these?"
        throw(e)
    end

    # code defining fit method:
    program1 = quote

        $(isdefined(modl, :MLJ) ? :(import MLJ.MLJBase) : :(import MLJBase))
        $(isdefined(modl, :MLJ) ? :(import MLJ.MLJBase.MLJModelInterface) :
          :(import MLJBase.MLJModelInterface))

        $instance = $modeltype_ex()
        $models = [getproperty($instance, name)
                   for name in fieldnames($modeltype_ex)]

        MLJModelInterface.fit(model::$modeltype_ex, verb::Integer, $args...) =
            MLJBase.fit_method($mach_ex, $models...)(model, verb, $args...)

    end

    modl.eval(program1)

    # define composite model traits:
    for (name_ex, value_ex) in trait_ex_given_name_ex
        program = quote
            MLJBase.$name_ex(::Type{<:$modeltype_ex}) = $value_ex
        end
        modl.eval(program)
    end

    return nothing

end

"""

    @from_network mach [mutable] struct NewCompositeModel
           ...
    end

or

    @from_network mach begin
        [mutable] struct NewCompositeModel
           ...
        end
        <optional trait declarations>
    end

Create a new stand-alone model type called `NewCompositeModel`, using
the specified learning network machine `mach` as a blueprint.

For more on learning network machines, see [`machine`](@ref).

### Example

Consider the following simple learning network for training a decision
tree after one-hot encoding the inputs, and forcing the predictions to
be point-predictions (rather than probabilistic):

```julia
Xs = source()
ys = source()

hot = OneHotEncoder()
tree = DecisionTreeClassifier()

W = transform(machine(hot, Xs), Xs)
yhat = predict_mode(machine(tree, W, ys), W)
```

A learning network machine is defined by

```julia
mach = machine(Deterministic(), Xs, ys; predict=yhat)
```

To specify a new `Deterministic` composite model type `WrappedTree` we
specify the model instances appearing in the network as "default"
values in the following decorated struct definition:

```julia
@from_network mach struct WrappedTree
    encoder=hot
    decision_tree=tree
end
```
and create a new instance with `WrappedTree()`.

To allow the second model component to be replaced by any other
probabilistic model we instead make a mutable struct declaration and,
if desired, annotate types appropriately.  In the following code
illustration some model trait declarations have also been added:

```julia
@from_network mach begin
    mutable struct WrappedTree
        encoder::OneHotEncoder=hot
        classifier::Probabilistic=tree
    end
    input_scitype = Table(Continuous, Finite)
    is_pure_julia = true
end
```

"""
macro from_network(exs...)
    args = from_network_preprocess(__module__, exs...)
    modeltype_ex = args[2]
    from_network_(__module__, args...)
end


text:
import Pkg
import Random

@inline function _precompile_execution_content(pkgs::AbstractVector{<:Symbol};
                                               notest::AbstractVector{<:Symbol} = Symbol[],
                                               test_stdlibs::Bool = false)
    stdlib_names = Symbol.(collect(values(Pkg.Types.stdlibs())))
    result = "\n"
    result *= "import Pkg\n"
    result *= "import Test\n"
    for pkg in pkgs
        result *= "import $(pkg)\n"
        pkg_is_stdlib = pkg in stdlib_names
        pkg_is_not_stdlib = !pkg_is_stdlib
        pkg_is_in_notest = pkg in notest
        pkg_is_not_in_notest = !pkg_is_in_notest
        if pkg_is_not_in_notest && (test_stdlibs || pkg_is_not_stdlib)
            module_name = "Test_$(pkg)_$(Random.randstring(16))"
            result *= "module $(module_name)\n"
            result *= "import $(pkg)\n"
            result *= "if Base.pkgdir($(pkg)) !== nothing && isfile(joinpath(Base.pkgdir($(pkg)), \"test\", \"runtests.jl\"))\n"
            result *= "include(joinpath(Base.pkgdir($(pkg)), \"test\", \"runtests.jl\"))\n"
            result *= "end # end if\n"
            result *= "end # end module $(module_name)\n"
            result *= "import .$(module_name)\n"
        end
    end
    return result
end

@inline function _precompile_execution(filename::AbstractString,
                                       pkgs::AbstractVector{<:Symbol};
                                       notest::AbstractVector{<:Symbol} = Symbol[],
                                       test_stdlibs::Bool = false)
    precompile_execution_content =_precompile_execution_content(pkgs;
                                                                notest = notest,
                                                                test_stdlibs = test_stdlibs)
    rm(filename; force = true, recursive = true)
    mkpath(dirname(filename))
    open(filename, "w") do io
        print(io, precompile_execution_content)
    end
    return filename
end

@inline function _precompile_execution(pkgs::AbstractVector{<:Symbol};
                                       notest::AbstractVector{<:Symbol} = Symbol[],
                                       test_stdlibs::Bool = false)
    my_temp_dir = mktempdir(; cleanup=true)
    my_temp_name = joinpath(my_temp_dir, "precompile_execution_file.jl")
    return _precompile_execution(my_temp_name,
                                 pkgs;
                                 notest = notest,
                                 test_stdlibs = test_stdlibs)
end


text:
# RolloutSimulator
# maintained by @zsunberg

"""
A fast simulator that just returns the reward

The simulation will be terminated when either
1) a terminal state is reached (as determined by `isterminal()` or
2) the discount factor is as small as `eps` or
3) max_steps have been executed

Keyword Arguments:
    - eps
    - max_steps

<<<<<<< HEAD
    initial_state

    eps

    max_steps
=======
Usage (optional arguments in brackets):
    ro = RolloutSimulator()
    history = simulate(ro, pomdp, policy, [updater [, init_belief [, init_state]]])
>>>>>>> 4dac41f2cbdbeeb4e91d93f1b69bd000bce163b0
"""
struct RolloutSimulator{RNG<:AbstractRNG} <: Simulator
    rng::RNG

    # optional: if these are null, they will be ignored
    max_steps::Nullable{Int}
    eps::Nullable{Float64}

    # DEPRECATED: remove in v0.2.8 or higher
    initial_state::Nullable{Any}
end

# These are the only safe constructors to use
RolloutSimulator(rng::AbstractRNG, d::Int=typemax(Int)) = RolloutSimulator(rng, Nullable{Int}(d), Nullable{Float64}(), Nullable{Any}())
function RolloutSimulator(;rng=MersenneTwister(rand(UInt32)),
                           initial_state=Nullable{Any}(),
                           eps=Nullable{Float64}(),
                           max_steps=Nullable{Int}())
    if !isnull(initial_state)
        warn("The initial_state argument for RolloutSimulator is deprecated. The initial state should be specified as the last argument to simulate(...).")
    end
    return RolloutSimulator{typeof(rng)}(rng, max_steps, eps, initial_state)
end

# COMPATIBILITY CONSTRUCTOR: DO NOT USE!
# Once version 0.2.7 is registered, start having this throw a warning in master
function RolloutSimulator(rng::AbstractRNG, is::Nullable{Any}, eps::Nullable{Float64}, ms::Nullable{Int})
    return RolloutSimulator(rng, ms, eps, is)
end

@POMDP_require simulate(sim::RolloutSimulator, pomdp::POMDP, policy::Policy) begin
    @req updater(::typeof(policy))
    bu = updater(policy)
    @subreq simulate(sim, pomdp, policy, bu)
end

@POMDP_require simulate(sim::RolloutSimulator, pomdp::POMDP, policy::Policy, bu::Updater) begin
    @req initial_state_distribution(::typeof(pomdp))
    dist = initial_state_distribution(pomdp)
    @subreq simulate(sim, pomdp, policy, bu, dist)
end

function simulate(sim::RolloutSimulator, pomdp::POMDP, policy::Policy, bu::Updater=updater(policy))
    dist = initial_state_distribution(pomdp)
    return simulate(sim, pomdp, policy, bu, dist)
end

@POMDP_require simulate(sim::RolloutSimulator, pomdp::POMDP, policy::Policy, updater::Updater, initial_belief) begin
    @req rand(::typeof(sim.rng), ::typeof(initial_belief))
    @subreq simulate(sim, pomdp, policy, updater, initial_belief, s)
end

function simulate{S}(sim::RolloutSimulator, pomdp::POMDP{S}, policy::Policy, updater::Updater, initial_belief)

    if !isnull(sim.initial_state)
        s = convert(S, get(sim.initial_state))::S
    else
        s = rand(sim.rng, initial_belief)::S
    end

    return simulate(sim, pomdp, policy, updater, initial_belief, s)
end

@POMDP_require simulate(sim::RolloutSimulator, pomdp::POMDP, policy::Policy, updater::Updater, initial_belief, s) begin
    P = typeof(pomdp)
    S = state_type(P)
    A = action_type(P)
    O = obs_type(P)
    @req initialize_belief(::typeof(updater), ::typeof(initial_belief))
    @req isterminal(::P, ::S)
    @req discount(::P)
    @req generate_sor(::P, ::S, ::A, ::typeof(sim.rng))
    b = initialize_belief(updater, initial_belief)
    @req action(::typeof(policy), ::typeof(b))
    @req update(::typeof(updater), ::typeof(b), ::A, ::O)
end

function simulate(sim::RolloutSimulator, pomdp::POMDP, policy::Policy, updater::Updater, initial_belief, s)
    if isterminal(pomdp, s)
        return 0.0
    end

    eps = get(sim.eps, 0.0)
    max_steps = get(sim.max_steps, typemax(Int))

    disc = 1.0
    r_total = 0.0

    b = initialize_belief(updater, initial_belief)
    a = nothing
    o = nothing

    step = 1

    while disc > eps && !isterminal(pomdp, s, a, o) && step <= max_steps # TODO also check for terminal observation
        a = action(policy, b)

        sp, o, r = generate_sor(pomdp, s, a, sim.rng)

        r_total += disc*r
        # if isterminal(pomdp, sp, a, o)
        #     break
        # end

        s = sp

        bp = update(updater, b, a, o)
        b = bp

        disc *= discount(pomdp)
        step += 1
    end

    return r_total
end

@POMDP_require simulate(sim::RolloutSimulator, mdp::MDP, policy::Policy) begin
    if isnull(sim.initial_state)
        @req initial_state(::typeof(mdp), ::typeof(sim.rng))
    end
    istate = initial_state(mdp, sim.rng)
    @subreq simulate(sim, mdp, policy, istate)
end

@POMDP_require simulate(sim::RolloutSimulator, mdp::MDP, policy::Policy, initial_state) begin
    P = typeof(mdp)
    S = typeof(initial_state)
    A = action_type(mdp)
    @req isterminal(::P, ::S)
    @req action(::typeof(policy), ::S)
    @req generate_sr(::P, ::S, ::A, ::typeof(sim.rng))
    @req discount(::P)
end

function simulate(sim::RolloutSimulator, mdp::MDP, policy::Policy)
    istate=get(sim.initial_state, initial_state(mdp, sim.rng))
    simulate(sim, mdp, policy, istate)
end

function simulate{S}(sim::RolloutSimulator, mdp::Union{MDP{S}, POMDP{S}}, policy::Policy, initial_state::S)
    if isterminal(mdp, initial_state)
        return 0.0
    end

    eps = get(sim.eps, 0.0)
    max_steps = get(sim.max_steps, typemax(Int))

    s = initial_state
    a = nothing

    disc = 1.0
    r_total = 0.0
    step = 1

    while disc > eps && !isterminal(mdp, s, a) && step <= max_steps
        a = action(policy, s)

        sp, r = generate_sr(mdp, s, a, sim.rng)

        r_total += disc*r
        # if isterminal(mdp, sp, a)
        #     break
        # end

        s = sp

        disc *= discount(mdp)
        step += 1
    end

    return r_total
end


text:
module SymbolicAD

import Base.Meta: isexpr
import JuMP
import LinearAlgebra
import MathOptInterface
import SparseArrays
import Symbolics
import Test

const MOI = MathOptInterface

include("nonlinear_oracle.jl")
include("MOI_wrapper.jl")
include("jump_test_utilities.jl")
end


text:
using BigNums
using BigNums: ArbDigit, DoubleArbDigit
using Test
using PropCheck
using Random

### Function Defs

include("genfuncs.jl")
include("properties.jl")

### Tests ###

@testset "Tests" begin
@testset "Construction" begin
    @test check(isnormalized, genArbUInt())
    @test checkOne()
    @test checkZero()
    @test check(setOne, genArbUInt())
    @test check(setZero, genArbUInt())
    @testset "constructNormalized($T)" for (gen,T) in ((PropCheck.vector(5, igen(ArbDigit)),Vector{ArbDigit}),
                                                        (igen(ArbDigit),ArbDigit),
                                                        (igen(DoubleArbDigit),DoubleArbDigit))
        @test check(constructNormalized, gen)
    end
end
@testset "Addition" begin
    @testset "$f" for (n,f) in ((2,commutative_add),
                                (3,associative_add),
                                (1,identity_add))
        gen = PropCheck.tuple(n, genArbUInt())
        @test check(Base.splat(f), gen)
    end
end
@testset "Multiplication" begin
    @testset "long_kara_eq" begin
        gen = PropCheck.tuple(2, PropCheck.vector(Integrated(Generator{UInt}(rng -> rand(rng,33:65))), igen(ArbDigit)))
        @test check(Base.splat(long_karatsuba_eq), gen)
    end
    @testset "$f" for (n,f) in ((2,commutative_mul),
                            (3,associative_mul),
                            (1,identity_mul),
                            (1,zero_mul))
        @testset for i in (1:32, 33:64) #, 65:1024)
            sizeGen = Integrated(Generator{UInt}(
                rng -> unsigned(rand(rng, i))
            ))
            gen = PropCheck.tuple(n, genArbUInt(sizeGen))
            @test check(Base.splat(f), gen)
        end
    end
end
@testset "Bitwise" begin
    @testset "Shifts" begin
        @testset "$leftShiftNoDataLoss" begin
            gen = interleave(genArbUInt(), igen(UInt8))
            @test check(Base.splat(leftShiftNoDataLoss), gen)
        end
        @testset "$leftShiftTrailLeadCorrect" begin
            gen = interleave(genArbUInt(), igen(UInt8)) # shifting more would be irresponsible
            @test check(Base.splat(leftShiftTrailLeadCorrect), gen)
        end
        @testset "$additionOfShifts" begin
            gen = interleave(genArbUInt(), igen(Tuple{UInt8, UInt8}))
            @test check(Base.splat(additionOfShifts), gen)
        end
        @testset "$shiftAntiInverse" begin
            gen = interleave(genArbUInt(), igen(UInt8)) # shifting more would be irresponsible
            @test check(Base.splat(shiftAntiInverse), gen)
        end
        @testset "$shiftValueSubtraction" begin
            shiftGen = filter(issorted, igen(Tuple{UInt8, UInt8}))
            gen = interleave(genArbUInt(), shiftGen)
            @test check(Base.splat(shiftValueSubtraction), gen)
        end
        @testset "$rightShiftSame" begin
            gen = interleave(genArbUInt(), igen(UInt8)) # shifting more would be irresponsible
            @test check(Base.splat(rightShiftSame), gen) broken=true
        end
        @testset "$shiftIdentity" begin
            gen = genArbUInt()
            @test check(shiftIdentity, gen)
        end
        @testset "$lossyShift" begin
            gen = interleave(genArbUInt(), igen(UInt8))
            @test check(Base.splat(lossyShift), gen)
        end
    end
    @testset "negation" begin
        @testset "$inverse" begin
            gen = genArbUInt()
            @test check(inverse, gen)
        end
    end
    @testset "and" begin
        @testset "$f" for (n,f) in ((2,commutative(&)),
                                    (1,zeroing_and),
                                    (1,identity_and))
            gen = PropCheck.tuple(n, genArbUInt())
            @test check(Base.splat(f), gen)
        end
    end
    @testset "or" begin
        @testset "$f" for (n,f) in ((2,identity_or),
                                    (2,commutative(|)))
            gen = PropCheck.tuple(n, genArbUInt())
            @test check(Base.splat(f), gen)
        end
    end
    @testset "xor" begin
        @testset "$zeroing_xor" begin
            gen = genArbUInt()
            @test check(zeroing_xor, gen)
        end
        @testset "$identity_xor" begin
            gen = genArbUInt()
            @test check(identity_xor, gen)
        end
        @testset "$commutative" begin
            gen = PropCheck.tuple(2, genArbUInt())
            @test check(Base.splat(commutative(⊻)), gen)
        end
    end
    @testset "nor" begin
        @testset "$f" for (n,f) in ((2,chaining_nor),
                                    (2,commutative(⊽)))
            gen = PropCheck.tuple(n, genArbUInt())
            @test check(Base.splat(f), gen)
        end
    end
    @testset "nand" begin
        @testset "$f" for (n,f) in ((2,chaining_nand),
                                    (2,commutative(⊼)))
            gen = PropCheck.tuple(n, genArbUInt())
            @test check(Base.splat(f), gen)
        end
    end
end
@testset "Ordering" begin
    @testset "$orderTransitive" begin
        gen = interleave(genArbUInt(), genArbUInt(), genArbUInt())
        @test check(Base.splat(orderTransitive), gen)
    end
    @testset "$orderReversal" begin
        gen = interleave(genArbUInt(), genArbUInt())
        @test check(Base.splat(orderReversal), gen)
    end
    @testset "$trichotomy" begin
        gen = interleave(genArbUInt(), genArbUInt())
        @test check(Base.splat(trichotomy), gen)
    end
    @testset "$orderPreservedAddition" begin
        genPair = Integrated(Generator{Tuple{ArbUInt, ArbUInt}}() do rng
            upper = generate(rng, ArbUInt)
            lower = rand(rng, zero(ArbUInt):upper)
            return (lower, upper + 0x1) # guarantee that the second one is larger
        end)
        gen = interleave(genPair, genArbUInt())
        @test check(Base.splat(orderPreservedAddition), gen)
    end
end
@testset "Strings & parsing" begin
    @testset "parse(string(..))" begin
        gen = genArbUInt(igen(0x100))
        @test check(x -> stringparse(x, 16), gen)
    end
    hexchar(c) = map(x -> Char(x < 0xa ? x+0x30 : x-0xa+0x61), c)
    elgen = Integrated(Generator{UInt8}(rng -> rand(rng, 0x0:0xf)))
    len = igen(unsigned(1000))
    elVectorGen = filter(x -> !iszero(length(x)), PropCheck.vector(len, elgen))
    @testset "string(parse(..))" begin
        validHexString = map(String ∘ hexchar, elVectorGen)

        @test check(x -> parsestring(x, 16), validHexString)
    end
    @testset "throwing inputs" begin
        @test_throws ArgumentError parse(ArbUInt, "")

        nonHex = filter(!isxdigit, '\u00':'\uff')
        invalidHexString = map(elVectorGen, String) do charVec
            retVec = map(hexchar, charVec)
            targetIdx = rand(eachindex(retVec))
            retVec[targetIdx] = rand(nonHex)
            join(retVec)
        end
        @test check(throwNonHex, invalidHexString)
    end
end
end


text:
fname = :ST_TRdense
c = Combi(hessian_dense,PDataST,solve_modelST_TR,preprocessST,decreaseGen,TparamsST())
include("Template.jl")


text:
module UnROOT

using LazyArrays
import Mmap: mmap
export ROOTFile, LazyBranch, LazyTree

import Base: close, keys, get, getindex, getproperty, show, length, iterate, position, ntoh, reinterpret
ntoh(b::Bool) = b

import AbstractTrees: children, printnode, print_tree

using CodecLz4, CodecXz, CodecZstd, StaticArrays, LorentzVectors, ArraysOfArrays
using Mixers, Parameters, Memoization, LRUCache
import IterTools: groupby

import LibDeflate: zlib_decompress!, Decompressor

import Tables, TypedTables, PrettyTables

@static if VERSION < v"1.6"
    Base.first(itr, n::Integer) = collect(Iterators.take(itr, n))
    Base.first(a::AbstractVector{S}, n::Integer) where S<: AbstractString = a[1:(length(a) > n ? n : end)]
    Base.first(a::S, n::Integer) where S<: AbstractString = a[1:(length(a) > n ? n : end)]
end

"""
    OffsetBuffer

Works with seek, position of the original file. Think of it as a view of IOStream that can be
indexed with original positions.
"""
struct OffsetBuffer{T}
    io::T
    offset::Int
end
Base.read(io::OffsetBuffer, nb) = Base.read(io.io, nb)
Base.seek(io::OffsetBuffer, i) = Base.seek(io.io, i - io.offset)
Base.skip(io::OffsetBuffer, i) = Base.skip(io.io, i)
Base.position(io::OffsetBuffer) = position(io.io) + io.offset

include("constants.jl")
include("streamsource.jl")
include("io.jl")
include("types.jl")
include("utils.jl")
include("streamers.jl")
include("bootstrap.jl")
include("root.jl")
include("iteration.jl")
include("custom.jl")
include("displays.jl")

end # module


text:
using Nemo, GibbsTypePriors, JLD

grid_k = [collect(1:25:10000); 10000]

Pkn_numerical_accuracy = GibbsTypePriors.Pkn_NGG_arb.(grid_k, 10000, 1.2, 0.6) |> x -> accuracy_bits.(x)
save("test/graphical_tests/saves_for_graphical_tests/accuracy_Pkn_10000.jld", "Pkn_numerical_accuracy", Pkn_numerical_accuracy)

max_stars_repo_name: konkam/GibbsTypePriors

max_issues_repo_name: konkam/GibbsTypePriors

max_forks_repo_name: konkam/GibbsTypePriors


text:
function test_mocha_kernels(backend::Backend, data_type)
  println("-- Testing Mocha CUDA kernels on $data_type")

  println("    > logistic loss forward")
  eps = 1e-5

  h, w, c, n = (5, 6, 7, 128)
  prob = abs(rand(data_type, (w,h,c,n))) + 0.1
  label = abs(rand(Int, (w, h, 1, n))) % c
  label = convert(Array{data_type}, label)

  prob_blob = make_blob(backend, data_type, size(prob))
  copy!(prob_blob, prob)
  label_blob = make_blob(backend, data_type, size(label))
  copy!(label_blob, label)

  # This should always be float32
  loss_blob = make_blob(backend, Float32, 1, 1, 1, 1)
  copy!(loss_blob, Float32[0])

  spatial_dim = h*w
  prob_dim = c

  x_block = round(Int, ceil(Float64(n)/CUDA.THREADS_PER_BLOCK_X))
  y_block = spatial_dim

  if data_type == Float32
    kernel = backend.mocha.logistic_loss_forward_float
  else
    kernel = backend.mocha.logistic_loss_forward_double
  end
  weights = convert(Ptr{data_type}, 0)

  CUDA.launch(kernel, (x_block, y_block), (CUDA.THREADS_PER_BLOCK_X, 1),
      (prob_blob.ptr.p, label_blob.ptr.p, weights, n, spatial_dim, prob_dim, loss_blob.ptr.p))

  loss = Float32[0]
  copy!(loss, loss_blob)
  loss = loss[1] / (n*h*w)

  # simulate CUDA blocking, because float is very inaccurate
  # adding 512 float32 of -log(0.1) will produce noticeable different
  # results (1e-2) from adding 256 of them and then add the two sums
  expected_loss = 0f0
  for j = 1:w
    for k = 1:h
      local_loss = 0f0
      for i = 1:n
        local_loss += -log(prob[j, k, round(Int, label[j,k,1,i])+1, i])
      end
      expected_loss += local_loss
    end
  end
  expected_loss /= (n*h*w)

  #println("loss = $loss")
  #println("expected_loss = $expected_loss")
  @test (-eps < expected_loss - loss < eps)
end

function test_mocha_kernels(backend::Backend)
  test_mocha_kernels(backend, Float32)
  test_mocha_kernels(backend, Float64)
end

test_mocha_kernels(backend_gpu)


text:
# Make some strings
name1 = "Viki"
name2 = "Noah"
phrase = " is friends with "

# Concatenate strings
first_sentence = name1 * phrase * name2
println(first_sentence)

# Concatenate a different string
second_sentence = name2 * phrase * name1
println(second_sentence)

# Now let's repeat strings
rep1 = "my " * "my " * "my "
println(rep1)
rep2 = "my "^3
println(rep2)

# Use both the * and ^ operator
ops = "my "^3 * "what big teeth you have"
println(ops)

# Concatenate variables
x, y, z = "what's this? ", "what's this? ", "there's magic in the air!"
println(x*y*z)

# We can also repeat multiple strings
println((x*y*z)^3)

# If using command line, use ; to suppress outputs like so
# x, y, z = "string 1", "this is a string", "but this? what's this?";
# x*y*z

# We can also use characters outside of those available on international keyboards
# On the command line, type
# \alpha
# Then click tab
# You can directly use the unicode value by typing the following (single quotes must be included)
# '\u03b1'
# The second method allows for use in scripts
println('\u03b1')

# In Julia, we can add numbers to characters to give us the corresponding letter in the alphabet
# This is based on the unicode, where 'a' is \u0061 and 'b' is \u0062
println('a' + 1)

max_stars_repo_name: KorfLab/julie

max_issues_repo_name: KorfLab/julie

max_forks_repo_name: KorfLab/julie


text:
using Pluto
Pluto.run(;host="0.0.0.0", port=1234, launch_browser=false, require_secret_for_access=false)


text:
@testset "Rastrigin" begin

    function test_result(result::Evolutionary.EvolutionaryOptimizationResults, N::Int, tol::Float64)
        fitness = minimum(result)
        extremum = Evolutionary.minimizer(result)
        if round(fitness) == 0
            @test extremum ≈ zeros(N) atol=tol
            @test fitness ≈ 0.0 atol=tol
        else
            # @warn("Found local minimum!!!")
            @test sum(abs, extremum) < N
        end
    end

    # Objective function
    function rastrigin(x::AbstractVector{T}) where {T <: AbstractFloat}
        n = length(x)
        return 10n + sum([ x[i]^2 - 10cos(convert(T,2π*x[i])) for i in 1:n ])
    end

    # Parameters
    N = 3
    P = 100
    initState = ()->rand(N)

    # Testing: (μ/μ_I,λ)-σ-Self-Adaptation-ES
    # with non-isotropic mutation operator y' := y + (σ_1 N_1(0, 1), ..., σ_N N_N(0, 1))
    m = ES(mu = 15, lambda = P)
    @test m.μ == 15
    @test m.λ == P
    result = Evolutionary.optimize( rastrigin,
        initState,
        ES(
            initStrategy = AnisotropicStrategy(N),
            recombination = average, srecombination = average,
            mutation = gaussian, smutation = gaussian,
            selection=:comma,
            μ = 15, λ = P
        ),Evolutionary.Options(iterations=1000, show_trace=false)
    )
    println("(15,$(P))-σ-SA-ES => F: $(minimum(result)), C: $(Evolutionary.iterations(result))")
    test_result(result, N, 1e-1)

    # Testing: CMA-ES
    result = Evolutionary.optimize(rastrigin, initState, CMAES(mu = 15, lambda = P))
    println("(15/15,$(P))-CMA-ES => F: $(minimum(result)), C: $(Evolutionary.iterations(result))")
    test_result(result, N, 1e-1)

    # Testing: GA
    m = GA(epsilon = 10.0)
    @test m.ɛ == 10

    selections = [:roulette=>rouletteinv, :sus=>susinv, :rank=>ranklinear(1.5)]
    crossovers = [:discrete=>discrete, :intermediate0=>intermediate(0.), :intermediate0_25=>intermediate(0.5), :line=>line(0.2)]
    mutations = [:domrng0_5=>domainrange(fill(0.5,N)), :uniform=>uniform(3.0), :gaussian=>gaussian(0.6)]

    @testset "GA settings" for (sn,ss) in selections, (xn,xovr) in crossovers, (mn,ms) in mutations
        xn == :discrete && (mn == :uniform || mn == :gaussian) && continue # bad combination
        xn == :line && mn == :gaussian && continue # bad combination
        result = Evolutionary.optimize( rastrigin, initState,
            GA(
                populationSize = P,
                ɛ = 0.1,
                selection = ss,
                crossover = xovr,
                mutation = ms
            ), Evolutionary.Options(iterations=1000, successive_f_tol=25)
        )
        println("GA:$(sn):$(xn):$(mn)(N=$(N),P=$(P),x=.8,μ=.1,ɛ=0.1) => F: $(minimum(result)), C: $(Evolutionary.iterations(result))")
        test_result(result, N, 1e-1)
    end

    # Testing: DE
    selections = [:rand=>random, :perm=>permutation, :rndoff=>randomoffset, :best=>best]
    mutations = [:exp=>exponential(0.5), :bin=>uniformbin(0.5)]

    @testset "DE settings" for (sn,ss) in selections, (mn,ms) in mutations, n in 1:2
        result = Evolutionary.optimize( rastrigin, initState,
            DE(
                populationSize = P,
                n=n,
                selection = ss,
                recombination = ms,
                F = 0.9
            )
        )
        println("DE/$sn/$n/$mn(F=0.9,Cr=0.5) => F: $(minimum(result)), C: $(Evolutionary.iterations(result))")
        test_result(result, N, 1e-2)
    end

end


text:
function ludic_filter{T<:Integer}(n::T)
    0 < n || throw(DomainError())
    slud = trues(n)
    for i in 2:(n-1)
        slud[i] || continue
        x = 0
        for j in (i+1):n
            slud[j] || continue
            x += 1
            x %= i
            x == 0 || continue
            slud[j] = false
        end
    end
    return slud
end

ludlen = 10^5
slud = ludic_filter(ludlen)
ludics = collect(1:ludlen)[slud]

n = 25
println("Generate and show here the first ", n, " ludic numbers.")
print("    ")
crwid = 76
wid = 0
for i in 1:(n-1)
    s = @sprintf "%d, " ludics[i]
    wid += length(s)
    if crwid < wid
        print("\n    ")
        wid = 0
    end
    print(s)
end
println(ludics[n])

n = 10^3
println()
println("How many ludic numbers are there less than or equal to ", n, "?")
println("    ", sum(slud[1:n]))

lo = 2000
hi = lo+5
println()
println("Show the ", lo, "..", hi, "'th ludic numbers.")
for i in lo:hi
    println("    Ludic(", i, ") = ", ludics[i])
end

n = 250
println()
println("Show all triplets of ludic numbers < ", n)
for i = 1:n-7
    slud[i] || continue
    j = i+2
    slud[j] || continue
    k = i+6
    slud[k] || continue
    println("    ", i, ", ", j, ", ", k)
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
"A `HyperRectangle` is the cartesian product of intervals."
abstract type HyperRectangle{T} <: ProductDomain{T} end

boundingbox(d::HyperRectangle) = d

"Compute all corners of the hyperrectangle."
function corners(d::HyperRectangle)
	left = infimum(d)
	right = supremum(d)
    N = length(left)
    corners = [zeros(numtype(d), N) for i in 1:2^N]
    # All possible permutations of the corners
    for i=1:2^length(left)
        for j=1:N
            corners[i][j] = ((i>>(j-1))%2==0) ? left[j] : right[j]
        end
    end
    corners
end

"A `Cube` is a hyperrectangle with equal side lengths in each dimension."
abstract type Cube{T} <: HyperRectangle{T} end

"A cube in a fixed N-dimensional Euclidean space."
const EuclideanCube{N,T} = Cube{SVector{N,T}}

"A cube with vector elements of variable length."
const VectorCube{T} = Cube{Vector{T}}

# This will cause a warning if vectors of different length are used with cubes
iscompatiblepair(x::AbstractVector, d::VectorCube) = length(x) == dimension(d)

"The unit cube is the domain `[0,1]^d`."
abstract type UnitCube{T} <: Cube{T} end

component(d::UnitCube, i::Int) =
    (1 <= i <= dimension(d) || throw(BoundsError); UnitInterval{numtype(d)}())

volume(d::UnitCube) = 1

"A unit cube that is specified by the element type `T`."
struct StaticUnitCube{T} <: UnitCube{T}
end

StaticUnitCube() = StaticUnitCube{SVector{3,Float64}}()
StaticUnitCube(::Val{N}) where {N} = StaticUnitCube{SVector{N,Float64}}()

StaticUnitCube{T}(n::Int) where {T} =
    (@assert n == euclideandimension(T); StaticUnitCube{T}())
StaticUnitCube{T}(::Val{N}) where {N,T} =
    (@assert N == euclideandimension(T); StaticUnitCube{T}())

similardomain(d::StaticUnitCube, ::Type{T}) where {T<:StaticTypes} =
    StaticUnitCube{T}()
similardomain(d::StaticUnitCube, ::Type{T}) where {T} =
    DynamicUnitCube{T}(dimension(d))

components(d::StaticUnitCube{SVector{N,T}}) where {N,T} =
    ntuple(x->UnitInterval{T}(), Val(N))

"The unit cube in a fixed N-dimensional space."
const EuclideanUnitCube{N,T} = StaticUnitCube{SVector{N,T}}

EuclideanUnitCube{N}() where {N} = EuclideanUnitCube{N,Float64}()

const UnitSquare{T} = EuclideanUnitCube{2,T}

"A unit cube whose dimension is specified by a field."
struct DynamicUnitCube{T} <: UnitCube{T}
    dimension   ::  Int

    DynamicUnitCube{T}(n::Int) where {T} = new(n)
    DynamicUnitCube{T}(n::Int) where {T<:StaticTypes} =
        (@assert n == euclideandimension(T); new(n))
end

DynamicUnitCube(n::Int) = DynamicUnitCube{Vector{Float64}}(n)

dimension(d::DynamicUnitCube) = d.dimension

components(d::DynamicUnitCube) = map(x->UnitInterval{numtype(d)}(), 1:dimension(d))

similardomain(d::DynamicUnitCube, ::Type{T}) where {T} =
    DynamicUnitCube{T}(d.dimension)
similardomain(d::DynamicUnitCube, ::Type{T}) where {T <: StaticTypes} =
    StaticUnitCube{T}()

"The unit cube with vector elements of a given dimension."
const VectorUnitCube{T} = DynamicUnitCube{Vector{T}}

VectorUnitCube(n::Int = 3) = VectorUnitCube{Float64}(n)
VectorUnitSquare() = VectorUnitCube(2)

UnitCube(n::Int) = DynamicUnitCube(n)
UnitCube(::Val{N} = Val(3)) where {N} = EuclideanUnitCube{N}()

UnitCube{T}(n::Int) where {T <: StaticTypes} = StaticUnitCube{T}(n)
UnitCube{T}(::Val{N}) where {N,T} = StaticUnitCube{T}(Val(N))
UnitCube{T}() where {T <: StaticTypes} = StaticUnitCube{T}()
UnitCube{T}(n::Int) where {T} = DynamicUnitCube{T}(n)

UnitCube(domains::UnitInterval...) = UnitCube(domains)
UnitCube(domains::NTuple{N,UnitInterval{T}}) where {N,T} =
    UnitCube{SVector{N,T}}(domains)
UnitCube(domains::SVector{N,UnitInterval{T}}) where {N,T} =
    UnitCube{SVector{N,T}}(domains)

UnitCube{T}(domains::UnitInterval...) where {T} = UnitCube{T}(domains)
UnitCube{T}(domain::NTuple{N,<:UnitInterval}) where {N,T<:SVector{N}} =
    StaticUnitCube{T}()
UnitCube{T}(domain::SVector{N,<:UnitInterval}) where {N,T<:SVector{N}} =
    StaticUnitCube{T}()

# Constructor: careful about ambiguities with FixedInterval arguments below
ProductDomain(domains::UnitInterval...) = UnitCube(domains...)
ProductDomain(domains::NTuple{N,<:UnitInterval}) where {N} = UnitCube(domains)
ProductDomain(domains::SVector{N,<:UnitInterval}) where {N} = UnitCube(domains)
ProductDomain(domains::AbstractVector{<:UnitInterval{T}}) where {T} =
    VectorUnitCube{T}(length(domains))
ProductDomain{T}(domains::UnitInterval...) where {N,S,T<:SVector{N,S}} =
	UnitCube{T}(domains...)
ProductDomain{T}(domains::NTuple{N,<:UnitInterval}) where {N,S,T<:SVector{N,S}} =
	UnitCube{T}(domains)
ProductDomain{T}(domains::SVector{N,<:UnitInterval}) where {N,S,T<:SVector{N,S}} =
	UnitCube{T}(domains)
ProductDomain{T}(domains::AbstractVector{<:UnitInterval}) where {S,T<:Vector{S}} =
	VectorUnitCube{S}(length(domains))

## Display:
show(io::IO, d::EuclideanUnitCube{3,Float64}) = print(io, "UnitCube()")
show(io::IO, d::EuclideanUnitCube{N,Float64}) where {N} = print(io, "UnitCube(Val($(N)))")
show(io::IO, d::UnitSquare{Float64}) = print(io, "UnitSquare()")
show(io::IO, d::UnitSquare{T}) where {T} = print(io, "UnitSquare{$(T)}()")
show(io::IO, d::VectorUnitCube{Float64}) = print(io, "UnitCube($(dimension(d)))")

# set the display stencils to [] to opt-out of composite display for the types above
Display.displaystencil(d::EuclideanUnitCube{N,Float64}) where {N} = []
Display.displaystencil(d::UnitSquare) = []
Display.displaystencil(d::VectorUnitCube{Float64}) = []
Display.object_parentheses(d::EuclideanUnitCube{N,Float64}) where {N} = false
Display.object_parentheses(d::UnitSquare) = false
Display.object_parentheses(d::VectorUnitCube{Float64}) = false

"An N-dimensional rectangle is the cartesian product of closed intervals."
struct Rectangle{T} <: HyperRectangle{T}
    a   ::  T
    b   ::  T

    function Rectangle{T}(a::S,b::S) where {S,T}
        @assert length(a)==length(b)
        new(a,b)
    end
end

dimension(d::Rectangle) = length(d.a)
component(d::Rectangle, i::Int) = ClosedInterval(d.a[i],d.b[i])
components(d::Rectangle) = map(ClosedInterval, d.a, d.b)

Rectangle(a, b) = Rectangle(promote(a,b)...)
Rectangle(a::T, b::T) where {T} = Rectangle{T}(a, b)
Rectangle(a::NTuple{N,T}, b::NTuple{N,T}) where {N,T} =
    Rectangle(SVector{N,T}(a), SVector{N,T}(b))

Rectangle(domains::Tuple) = Rectangle(domains...)
Rectangle(domains::ClosedInterval...) = Rectangle(promote_domains(domains)...)
Rectangle(domains::ClosedInterval{T}...) where {T} =
    Rectangle(map(infimum, domains), map(supremum, domains))
Rectangle(domains::AbstractVector{<:ClosedInterval}) =
    Rectangle(map(infimum, domains), map(supremum, domains))
Rectangle(domains::Domain...) =
    error("The Rectangle constructor expects two points or a list of intervals (closed).")

Rectangle{T}(domains::Tuple) where {T} = Rectangle{T}(domains...)
Rectangle{T}(domains::ClosedInterval...) where {T} =
	Rectangle{T}(map(infimum, domains), map(supremum, domains))
Rectangle{T}(domains::AbstractVector{<:ClosedInterval}) where {T} =
    Rectangle{T}(map(infimum, domains), map(supremum, domains))
Rectangle{T}(domains::Domain...) where {T} =
    error("The Rectangle constructor expects two points or a list of intervals (closed).")

ProductDomain(domains::ClosedInterval...) = Rectangle(domains...)
ProductDomain(domains::AbstractVector{<:ClosedInterval}) =
    Rectangle(map(infimum, domains), map(supremum, domains))
ProductDomain{T}(domains::ClosedInterval...) where {T} = Rectangle{T}(domains...)
ProductDomain{T}(domains::AbstractVector{<:ClosedInterval}) where {T} =
    Rectangle{T}(map(infimum, domains), map(supremum, domains))

"The N-fold cartesian product of a fixed interval."
struct FixedIntervalProduct{N,T,D} <: Cube{SVector{N,T}}
end

component(d::FixedIntervalProduct{N,T,D}, i::Int) where {N,T,D} =
    (1 <= i <= dimension(d) || throw(BoundsError); D())
components(d::FixedIntervalProduct{N,T,D}) where {N,T,D} =
    ntuple(x->D(), Val(N))

volume(d::FixedIntervalProduct{N,T,D}) where {N,T,D} = volume(D())^N

FixedIntervalProduct(domains::NTuple{N,D}) where {N,D <: FixedInterval} =
	FixedIntervalProduct{N,eltype(D),D}()
FixedIntervalProduct(domains::SVector{N,D}) where {N,D <: FixedInterval} =
	FixedIntervalProduct{N,eltype(D),D}()

const ChebyshevProductDomain{N,T} = FixedIntervalProduct{N,T,ChebyshevInterval{T}}
ChebyshevProductDomain(::Val{N}) where {N} = ChebyshevProductDomain{N}()
ChebyshevProductDomain{N}() where {N} = ChebyshevProductDomain{N,Float64}()

ProductDomain(domains::D...) where {D <: FixedInterval} =
	FixedIntervalProduct(domains)
ProductDomain(domains::NTuple{N,D}) where {N,D <: FixedInterval} =
	FixedIntervalProduct(domains)
ProductDomain(domains::SVector{N,<:FixedInterval}) where {N} =
	FixedIntervalProduct(domains)
ProductDomain{T}(domains::D...) where {N,S,T<:SVector{N,S},D <: FixedInterval} =
	FixedIntervalProduct(convert.(Ref(Domain{S}), domains))
ProductDomain{T}(domains::NTuple{N,D}) where {N,S,T<:SVector{N,S},D <: FixedInterval} =
	FixedIntervalProduct(convert.(Ref(Domain{S}), domains))
ProductDomain{T}(domains::SVector{N,<:FixedInterval}) where {N,S,T<:SVector{N}} =
	FixedIntervalProduct(domains)


text:
# vim: set et ts=2 sw=2:
r0() = Random.seed!(it.seed)

function ok()
  @testset "ALL" begin
    _lib(); _some(); _sym() end end

function _lib()
  r0()
  @testset "lib" begin
    @test few("abcdefgh",2) == ['b','c']
    @test thing("string") == "string"
    @test thing(11.5) == 11.5
    all = [row for row in csv("data/weather.csv")]
    @test 5       == length(all[1])
    @test 15      == length(all)
    @test Float64 == typeof(all[2][2]) end end

function _sym()
  @testset "some" begin
    s=  col(txt="a")
    @test typeof(s)==Sym
    for i in "aaaabbc" inc!(s,i) end
    @test s.mode == 'a'
    @test s.most == 4 end end

function _some()
  r0()
  @testset "some" begin
    s=  col(txt="<a")
    @test typeof(s)==Some
    s.max = 32
    for i in 1:1000 inc!(s,int(100*rand())) end
    lst = all(s)
    @test lst[1] < lst[end]
    @test mid(s,lo=16) == 69
    @test 32.37 < sd(s) < 32.38
    @test s.w == -1
    @test .78 < norm!(s,75) < .79 end end

max_stars_repo_name: timm/jut

max_issues_repo_name: timm/jut

max_forks_repo_name: timm/jut


text:
abstract type PassType end

struct OpaqueShadingPass <: PassType end
struct LightingPass <: PassType end
struct TransparentShadingPass <: PassType end

function Vk.PipelineColorBlendStateCreateInfo(pass::PassType)
    Vk.PipelineColorBlendStateCreateInfo(
        false,
        Vk.LOGIC_OP_AND,
        [Vk.PipelineColorBlendAttachmentState(pass)],
    )
end

function Vk.PipelineColorBlendAttachmentState(pass::OpaqueShadingPass)
    Vk.PipelineColorBlendAttachmentState(
        false,
        Vk.BLEND_FACTOR_ONE,
        Vk.BLEND_FACTOR_ZERO,
        Vk.BLEND_OP_ADD,
        Vk.BLEND_FACTOR_ONE,
        Vk.BLEND_FACTOR_ZERO,
        Vk.BLEND_OP_ADD;
        color_write_mask = Vk.COLOR_COMPONENT_R_BIT | Vk.COLOR_COMPONENT_G_BIT | Vk.COLOR_COMPONENT_B_BIT
    )
end

function Vk.PipelineColorBlendAttachmentState(pass::TransparentShadingPass)
    Vk.PipelineColorBlendAttachmentState(
        true,
        Vk.BLEND_FACTOR_SRC_ALPHA,
        Vk.BLEND_FACTOR_ONE_MINUS_SRC_ALPHA,
        Vk.BLEND_OP_ADD,
        Vk.BLEND_FACTOR_SRC_ALPHA,
        Vk.BLEND_FACTOR_ONE_MINUS_SRC_ALPHA,
        Vk.BLEND_OP_ADD;
        color_write_mask = Vk.COLOR_COMPONENT_R_BIT | Vk.COLOR_COMPONENT_G_BIT | Vk.COLOR_COMPONENT_B_BIT | Vk.COLOR_COMPONENT_A_BIT
    )
end


text:
using DynamicGrids, Test
using DynamicGrids: inbounds, isinbounds, _inbounds, _isinbounds, _cyclic_index,
    SimData, _unwrap, ismasked

@testset "boundary boundary checks are working" begin
    @testset "inbounds with Remove() returns index and false for an boundaryed index" begin
        @test _inbounds(Remove(), (4, 5), 1, 1) == ((1,1), true)
        @test _inbounds(Remove(), (4, 5), 2, 3) == ((2,3), true)
        @test _inbounds(Remove(), (4, 5), 4, 5) == ((4,5), true)
        @test _inbounds(Remove(), (4, 5), 0, 0) == ((0,0), false)
        @test _inbounds(Remove(), (3, 2), 2, 3) == ((2,3), false)
        @test _inbounds(Remove(), (1, 4), 2, 3) == ((2,3), false)
        @test _inbounds(Remove(), (2, 3), 200, 300) == ((200,300), false)
        @test _inbounds(Remove(), (4, 5), -3, -100) == ((-3,-100), false)
    end
    @testset "inbounds with Wrap() returns new index and true for an boundaryed index" begin
        @test _inbounds(Wrap(), (10, 10),  -2, 3) == ((8,  3), true)
        @test _inbounds(Wrap(), (10, 10),   2, 0) == ((2, 10), true)
        @test _inbounds(Wrap(), (10, 10),  22, 0) == ((2, 10), true)
        @test _inbounds(Wrap(), (10, 10), -22, 0) == ((8, 10), true)
    end
    @testset "isinbounds" begin
        @test _isinbounds((4, 5), 4, 5) == true
        @test _isinbounds((2, 3), 200, 300) == false
        @test _isinbounds((10, 10), -22, 0) == false
    end
    @testset "boundscheck objects" begin
        output = ArrayOutput(zeros(Int, 10, 10); tspan=1:10)
        sd = SimData(output.extent, Ruleset())
        @test inbounds(sd, 5, 5) == ((5, 5), true)
        @test inbounds(first(sd),  5, 5) == ((5, 5), true)
        @test inbounds(sd, 12, 5) == ((12, 5), false)
        sd_wrap = SimData(output.extent, Ruleset(; boundary=Wrap()))
        @test inbounds(sd_wrap, 5, 5) == ((5, 5), true)
        @test inbounds(sd_wrap, 12, 5) == ((2, 5), true)
        @test inbounds(first(sd_wrap), 12, 5) == ((2, 5), true)
    end
end

@testset "isinferred" begin
    @testset "unstable conditional" begin
        rule = let threshold = 20
            Cell() do x
                x > 1 ? 2 : 0.0
            end
        end
        output = ArrayOutput(rand(Int, 10, 10); tspan=1:10)
        @test_throws ErrorException isinferred(output, rule)
    end

    @testset "return type" begin
        rule = Neighbors{:a,:a}(Moore{1}(zeros(Bool, 3, 3))) do hood, x
            round(Int, x + sum(hood))
        end
        output = ArrayOutput((a=rand(Int, 10, 10),); tspan=1:10)
        @test isinferred(output, rule)
        output = ArrayOutput((a=rand(Bool, 10, 10),); tspan=1:10)
        @test_throws ErrorException isinferred(output, rule)
    end

    @testset "let blocks" begin
        a = 0.7
        rule = SetCell() do data, x, I
            add!(first(data), round(Int, a + x), I...)
        end
        output = ArrayOutput(zeros(Int, 10, 10); tspan=1:10)
        @test_throws ErrorException isinferred(output, Ruleset(rule))
        a = 0.7
        rule = let a = a
            SetCell() do data, x, I
                add!(first(data), round(Int, a), I...)
            end
        end
        output = ArrayOutput(zeros(Int, 10, 10); tspan=1:10)
        @test isinferred(output, Ruleset(rule))
    end

end

@testset "ismasked" begin
    output = ArrayOutput(zeros(2, 2); mask=Bool[1 0; 0 1], tspan=1:10)
    sd = SimData(output.extent, Ruleset())
    @test ismasked(sd, 1, 2) == true
    @test ismasked(sd, 2, 2) == false
    output_nomask = ArrayOutput(zeros(2, 2); tspan=1:10)
    sd = SimData(output_nomask.extent, Ruleset())
    @test ismasked(sd, 1, 2) == false
end

@testset "unwrap" begin
    @test _unwrap(1) == 1
    @test _unwrap(Val(:a)) == :a
    @test _unwrap(Aux(:a)) == :a
    @test _unwrap(Grid(:a)) == :a
    @test _unwrap(Aux{:x}) == :x
    @test _unwrap(Grid{:x}) == :x
    @test _unwrap(Val{:x}) == :x
end


text:
#DEPRECATED FILE

module CONCEPT

export CompositeType, Composite, AbstractComposite, TimeStamp, PreciseTime

end
push!(LOAD_PATH,joinpath(pwd(),"src","CONCEPT","common.jl"))


text:
"""
Truncation `IndividualsSelector`, with size.

"""
struct TruncationSelector <: IndividualsSelector
    k::Integer
end

#function select(::TruncationSelector, population::FitPopulation, n::Integer)
#    #TODO implement
#end


text:
module AbstractWrappedDataFrames

export AbstractWrappedDataFrame

using DataFrames
using Lazy: @forward

abstract type AbstractWrappedDataFrame <: AbstractDataFrame end

# Functions for AbstractWrappedDataFrame

# TODO Not ideal to hard code `df` field
dataframe(sf::AbstractWrappedDataFrame) = getfield(sf, :df)

# Implement the "unofficial" AbstractDataFrame interface
# See https://github.com/invenia/KeyedFrames.jl/issues/19#issuecomment-674753267

# Using Lazy.jl, we can forward these functions to the underlying data frame.

@forward AbstractWrappedDataFrame.df Base.getindex
@forward AbstractWrappedDataFrame.df Base.setindex!
@forward AbstractWrappedDataFrame.df Base.propertynames
@forward AbstractWrappedDataFrame.df Base.push!
@forward AbstractWrappedDataFrame.df Base.copy
@forward AbstractWrappedDataFrame.df Base.empty!
@forward AbstractWrappedDataFrame.df Base.hcat
@forward AbstractWrappedDataFrame.df Base.vcat
@forward AbstractWrappedDataFrame.df Base.sort!
@forward AbstractWrappedDataFrame.df Base.append!
@forward AbstractWrappedDataFrame.df Base.delete!
@forward AbstractWrappedDataFrame.df Base.parent
@forward AbstractWrappedDataFrame.df Base.parentindices

@forward AbstractWrappedDataFrame.df DataFrames.insertcols!
@forward AbstractWrappedDataFrame.df DataFrames.ncol
@forward AbstractWrappedDataFrame.df DataFrames.nrow
@forward AbstractWrappedDataFrame.df DataFrames.transform!
@forward AbstractWrappedDataFrame.df DataFrames.select!
@forward AbstractWrappedDataFrame.df DataFrames.index

@forward AbstractWrappedDataFrame.df DataFrames.DataFrameRow
# @forward AbstractWrappedDataFrame.df DataFrames.DataFrameRows
@forward AbstractWrappedDataFrame.df DataFrames.SubDataFrame
# @forward AbstractWrappedDataFrame.df DataFrames.GroupedDataFrame
# @forward AbstractWrappedDataFrame.df DataFrames.DataFrameColumns

# TODO This is not ideal because I would rather return only the properties
# of the underlying data frame. But if I don't support getting the object's
# own properties then I run into other problems with the `parent` function.
function Base.getproperty(sf::T, s::Symbol) where {T <: AbstractWrappedDataFrame}
    if s ∈ fieldnames(T)
        return getfield(sf, s)
    else
        return getproperty(dataframe(sf), s)
    end
end

# Custom forwarders to avoid ambiguity

# Avoid ambiguity since
# 1) Base defines setproperty!(::Any, ::Symbol, ::Any)
# 2) @forward defines setproperty!(::AbstractWrappedDataFrame, args...; kwargs...)
function Base.setproperty!(sf::AbstractWrappedDataFrame, s::Symbol, v::Any)
    return setproperty!(dataframe(sf), s, v)
end

# required by REPL completion.
# or get this error "propertynames(::SF, ::Bool) is ambiguous."
function Base.propertynames(sf::AbstractWrappedDataFrame, private::Bool)
    return propertynames(dataframe(sf), private)
end

# required by all joins
function Base.convert(::Type{S}, df::D) where {S <: AbstractWrappedDataFrame, D <: AbstractDataFrame}
    return S(df)
end

# required by eachcol
function Base.convert(::Type{D}, sf::S) where {S <: AbstractWrappedDataFrame, D <: AbstractDataFrame}
    return dataframe(sf)
end

# required by eachcol
function Base.convert(::Type{S}, sf::S) where {S <: AbstractWrappedDataFrame}
    return sf
end

# required by row indexing to return DataFrameRow
#=
julia> wdf[1, :]
ERROR: MethodError: getindex(::WDF, ::Int64, ::Colon) is ambiguous. Candidates:
  getindex(df::AbstractDataFrame, rowind::Integer, ::Colon) in DataFrames at /Users/tomkwong/.julia/packages/DataFrames/cdZCk/src/dataframerow/dataframerow.jl:128
  getindex(x::AbstractWrappedDataFrame, args...; kwargs...) in AbstractWrappedDataFrames at /Users/tomkwong/.julia/packages/Lazy/9Xnd3/src/macros.jl:297
  getindex(df::AbstractDataFrame, rowind::Integer, colinds::Union{Colon, Regex, AbstractArray{T,1} where T, All, Between, InvertedIndex}) in DataFrames at /Users/tomkwong/.julia/packages/DataFrames/cdZCk/src/dataframerow/dataframerow.jl:126
Possible fix, define
  getindex(::AbstractWrappedDataFrame, ::Integer, ::Colon)
=#
for v in [Colon, DataFrames.MultiColumnIndex]
    @eval function Base.getindex(s::AbstractWrappedDataFrame, i::Integer, arg::$v)
        return dataframe(s)[i, arg]
    end
end

# Custom extensions as these functions do not take AbstractDataFrame

using DataFrames: RowGroupDict, MultiColumnIndex

# required by semi/anti-joins
function DataFrames.findrow(gd::RowGroupDict,
                 sf::AbstractWrappedDataFrame,
                 args...)
    return DataFrames.findrow(gd, dataframe(sf), args...)
end

# required by inner/outer/left/right-joins
function DataFrames.findrows(gd::RowGroupDict,
                  sf::AbstractWrappedDataFrame,
                  args...)
    return DataFrames.findrows(gd, dataframe(sf), args...)
end

# required by select and many others
function DataFrames.manipulate(sf::AbstractWrappedDataFrame, args...; kwargs...)
    return DataFrames.manipulate(dataframe(sf), args...; kwargs...)
end

# required for groupby
# function DataFrames.SubDataFrame(sf::AbstractWrappedDataFrame, args...; kwargs...)
#     return DataFrames.SubDataFrame(dataframe(sf), args...; kwargs...)
# end

# required for `wdf[1, Between(:x, :y)] .= 4`
# function DataFrames.DataFrameRow(sf::AbstractWrappedDataFrame, args...; kwargs...)
#     return DataFrames.DataFrameRow(dataframe(sf), args...; kwargs...)
# end

end


text:
"""
    mutable struct DomainUpdater <: AbstractFieldUpdater

    DomainUpdater(; nper = 5, nskip = 100, lam = 0.1, pmax = 5.0, rlxn = 0.9)

"""
mutable struct DomainUpdater <: AbstractFieldUpdater
    nper      :: Int
    nskip     :: Int
    lam       :: Float64
    pmax      :: Float64
    rlxn      :: Float64

    oldparams :: Vector{Float64}
    step      :: Vector{Float64} # Newton step
    stress    :: Vector{Float64} # First derivative of free energy
    dstress   :: Matrix{Float64} # Jacobian matrix for cell parameters
    system    :: Option{FieldSystem}

    function DomainUpdater(; nper::Integer = 10, nskip::Integer = 200,
        lam::Real = 0.1, pmax::Real = 1.0, rlxn::Real = 0.5)
        return new(nper, nskip, lam, pmax, rlxn, [], [], [], zeros(0,0), nothing)
    end
end

#==============================================================================#

Base.show(io::IO, domain::DomainUpdater) = @printf(io, "DomainUpdater(nper = %d, nskip = %d)", domain.nper, domain.nskip)

stress_error(domain::DomainUpdater) = norm(domain.stress)

function setup!(domain::DomainUpdater, sys::FieldSystem)
    domain.system = sys

    npr = nparams(sys.cell)
    domain.stress = scfstress(sys)
    domain.oldparams = zeros(npr)
    domain.step = zeros(npr)
    domain.dstress = zeros(npr, npr)

    return nothing
end

function scfstress!(domain::DomainUpdater)
    @assert !isnothing(domain.system)
    domain.stress .= scfstress(domain.system)
    return nothing
end

function step!(domain::DomainUpdater; nsteps::Integer = 1, tol::Real = 1e-5)
    @assert !isnothing(domain.system)
    sys = domain.system
    cell = sys.cell

    # Calculate current stress
    scfstress!(domain)
    err = stress_error(domain)

    converged = false
    iter = 0

    while iter < nsteps && !converged
        dstress!(domain)
        if cell.dim == 1
            domain.step .= -1.0 * domain.stress[1] / domain.dstress[1]
        else
            mul!(domain.step, inv(domain.dstress), domain.stress)
            domain.step .*= -1.0
        end

        # Rescale step-size if needed
        step_size = norm(domain.step)
        if step_size > domain.pmax
            domain.step .*= domain.pmax / step_size
        end

        # Step cell parameters, update cell and stress
        @. cell.params += domain.lam * domain.step
        @. cell.params = domain.rlxn*cell.params + (1 - domain.rlxn)*domain.oldparams

        # Update to new stress
        update!(cell)
        density!(sys)
        residuals!(sys)
        scfstress!(domain)

        err = stress_error(domain)
        iter += 1
    end

    return nothing
end

function dstress!(domain::DomainUpdater)
    @assert !isnothing(domain.system)
    sys = domain.system
    cell = sys.cell
    domain.oldparams .= cell.params

    # Finite difference step size
    dx = 1e-8
    x = 1.0 / dx

    # For each cell param, calculate response to small perturbation
    for k = 1:nparams(cell)
        cell.params .= domain.oldparams
        cell.params[k] += dx

        update!(cell)
        density!(sys)
        pstress = scfstress(sys)

        @. domain.dstress[:,k] = x * (pstress - domain.stress)
    end

    # Return to original state
    cell.params .= domain.oldparams
    density!(sys)

    return nothing
end


text:
using ElectronDisplay

const tmp_dir = join([expanduser("~"),".DeepShiba"], "/")
const dot_file_path = join([tmp_dir,"tmp_graph.dot"], "/")

struct PNGContainer
    content
end

function Base.show(io::IO, ::MIME"image/png", c::PNGContainer)
    write(io, c.content)
end

function _dot_var(var::Variable, show_value)
    name = (var.name == "") ? "" : var.name * ":"
    value = (show_value == 0 ? var.data : var.grad.data)
    if var.data !== nothing
        var_size = size(value)
        if isempty(var_size)
            try var.data !== nothing
                name *= "$(value)"
            catch
                name *= "nothing"
            end
        else
            name *= "shape: $(var_size) \n type: $(get_value_type(value))"
        end
    end
    dot_var = "$(objectid(var)) [label=\"$name\", color=orange, style=filled]\n"
end

function _dot_func(f::Func)
    f_type = typeof(f)
    if f_type == DeepShiba.Pow
        f_type = "$(f_type)($(f.c))"
    end
    txt = "$(objectid(f)) [label=\"$(f_type)\", color=lightblue, style=filled, shape=box]\n"
    for x in f.inputs
        txt *= "$(objectid(x)) -> $(objectid(f))\n"
    end
    for y in f.outputs
        txt *= "$(objectid(f)) -> $(objectid(y))\n"
    end
    return txt
end

function get_dot_graph(var, show_value, title)
    txt = ""
    funcs = []
    seen_set = Set{Func}()
    push!(funcs, var.creator)
    txt = _dot_var(var, show_value)
    while !(isempty(funcs))
        f = pop!(funcs)
        txt *= _dot_func(f)
        for x in f.inputs
            txt *= _dot_var(x, show_value)
            if x.creator !== nothing && (!(x.creator in seen_set))
                push!(seen_set, x.creator)
                push!(funcs, x.creator)
            end
        end
    end
    return "digraph g {
            graph [
                labelloc=\"t\";
                label= \"$(title)\"
            ];
                $txt
            }"
end

function plot_tmp_dir(extension)
    extension = "png"
    to_file = join([tmp_dir, "graph.png"], "/")
    cmd = `dot $(dot_file_path) -T $(extension) -o $(to_file)`
    run(cmd)
    return to_file
end

function plot(var::Variable; to_file="", show_value="data", title="")
    if show_value == "grad"
        show_value = 1
    elseif show_value == "data"
        show_value = 0
    else
        error("The argument to show_value is invalid.")
    end

    dot_graph = get_dot_graph(var, show_value, title)
    (!(ispath(tmp_dir))) && (mkdir(tmp_dir))
    open(dot_file_path, "w") do io
        write(io, dot_graph)
    end

    if to_file == ""
        if split(PROGRAM_FILE, "/")[end] == "kernel.jl"
            png_file_path = plot_tmp_dir(".png")
            c = open(png_file_path) do io
                PNGContainer(read(io))
            end
            return c
        elseif PROGRAM_FILE == ""
            png_file_path = plot_tmp_dir(".png")
            c = open(png_file_path) do io
                PNGContainer(read(io))
            end
            display(c)
        else
            return c
        end
    else
        extension = split(to_file, ".")[2]
        cmd = `dot $(dot_file_path) -T $(extension) -o $(to_file)`
        run(cmd)
    end
end