text:
import HTTP:Response,   # For upgrade
            Request,    # For upgrade, target, origin, subprotocol
            Header,     # benchmarks, client_test, handshaketest
            header,     # For upgrade, subprotocol
            hasheader,  # For upgrade and check_upggrade
            setheader,  # For upgrade
            setstatus,  # For upgrade
            startwrite, # For upgrade
            startread,  # For _openstream
            handle,     # For _servercoroutine
            Connection, # For handshaketest
            Transaction,# For handshaketest
            URI,        # For open
            Handler,    # For WebSocketHandler, ServerWS, error_test
            RequestHandlerFunction,         # For ServerWS
            StatusError,                    # For open
            Servers,                        # For further imports
            Streams,                        # For further imports
            ConnectionPool                  # For further imports
import HTTP.ConnectionPool:
            getrawstream                    # For _openstream
import HTTP.Streams:
            Stream                          # For is_upgrade, handshaketest
import HTTP.Servers: MbedTLS                # For further imports
import HTTP.Servers.MbedTLS:
            SSLConfig                       # For ServerWS

"""
default_options() can be overruled by passing keyword arguments to the
ServerWS constructor.
# Options include
```
    sslconfig::Union{SSLConfig, Nothing}
    tcpisvalid::Function
    reuseaddr::Bool
    connection_count::Ref{Int}
    rate_limit::Union{Rational{Int}, Nothing}
    reuse_limit::Int
    readtimeout::Int
```
"""
default_options() = (in = Channel(1),
                    out = Channel(2),
                    sslconfig = nothing,
                    tcpisvalid = tcp->true,
                    reuseaddr = false,
                    connection_count = Ref(0),
                    rate_limit = nothing,
                    reuse_limit = typemax(Int),
                    readtimeout = 0)

"""
Initiate a websocket|client connection to server defined by url. If the server accepts
the connection and the upgrade to websocket, f is called with an open websocket|client

e.g. say hello, close and leave
```julia
using WebSockets
WebSockets.open("ws://127.0.0.1:8000") do ws
    write(ws, "Hello")
end;
```
If a server is listening and accepts, "Hello" is sent (as a CodeUnits{UInt8,String}).

On exit, a closing handshake is started. If the server is not currently reading
(which is a blocking function and it may be busy with other things), this side
will reset the underlying connection (ECONNRESET) after 10 seconds, and continue
execution without error.

"""
function open(f::Function, url; verbose=false, subprotocol = "", kw...)
    key = base64encode(rand(UInt8, 16))
    headers = [
        "Upgrade" => "websocket",
        "Connection" => "Upgrade",
        "Sec-WebSocket-Key" => key,
        "Sec-WebSocket-Version" => "13"
    ]
    if subprotocol != ""
        push!(headers, "Sec-WebSocket-Protocol" => subprotocol )
    end

    if in('#', url)
        throw(ArgumentError(" replace '#' with %23 in url: $url"))
    end
    uri = URI(url)
    if uri.scheme != "ws" && uri.scheme != "wss"
        throw(ArgumentError(" bad argument url: Scheme not ws or wss. Input scheme: $(uri.scheme)"))
    end
    openstream(stream) = _openstream(f, stream, key)
    try
        HTTP.open(
            openstream,
            "GET", uri, headers;
            reuse_limit=0, verbose=verbose ? 2 : 0, kw...)
    catch err
        if typeof(err) <: HTTP.IOExtras.IOError
            throw(WebSocketClosedError(" while open ws|client: $(string(err.e.msg))"))
        elseif typeof(err) <: HTTP.StatusError
            return err.response
        else
           rethrow(err)
        end
    end
end

"Called by open with a stream connected to a server, after handshake is initiated"
function _openstream(f::Function, stream, key::String)
    startread(stream)
    response = stream.message
    if response.status != 101
        return
    end
    check_upgrade(stream)
    if header(response, "Sec-WebSocket-Accept") != generate_websocket_key(key)
        throw(WebSocketError(0, "Invalid Sec-WebSocket-Accept\n" *
                                "$response"))
    end
    # unwrap the stream
    io = getrawstream(stream)
    ws = WebSocket(io, false)
    try
        f(ws)
    finally
        close(ws)
    end
end

"""
Not called from user code in most cases.

Called by serve > _servercoroutine, or in you own server coroutine definition if
    `serve(ServerWS,  host, port, verbose)`
does not suit your purpose. Call 'upgrade' if is_upgrade(stream.message) returns true.

Responds to a WebSocket handshake request.
If the connection is acceptable, sends status code 101
and headers according to RFC 6455, then calls
user's handler function f with the connection wrapped in
a WebSocket instance.

f(ws)           is called with the websocket and no client info
f(headers, ws)  also receives a dictionary of request headers for added security measures

On exit from f, a closing handshake is started. If the client is not currently reading
(which is a blocking function), this side will reset the underlying connection (ECONNRESET)
after 10 seconds and continue execution.

If the upgrade is not accepted, responds to client with '400'.

See /examples/server_listen_syntax.jl
```
"""
function upgrade(f::Function, stream)
    check_upgrade(stream)
    if !hasheader(stream, "Sec-WebSocket-Version", "13")
                    setheader(stream, "Sec-WebSocket-Version" => "13")
                    setstatus(stream, 400)
                    startwrite(stream)
        return
    end
    if hasheader(stream, "Sec-WebSocket-Protocol")
        requestedprotocol = header(stream, "Sec-WebSocket-Protocol")
        if !hasprotocol(requestedprotocol)
            setheader(stream, "Sec-WebSocket-Protocol" => requestedprotocol)
            setstatus(stream, 400)
            startwrite(stream)
            return
        else
            setheader(stream, "Sec-WebSocket-Protocol" => requestedprotocol)
        end
    end
    key = header(stream, "Sec-WebSocket-Key")
    decoded = UInt8[]
    try
        decoded = base64decode(key)
    catch
        setstatus(stream, 400)
        startwrite(stream)
        return
    end
    if length(decoded) != 16 # Key must be 16 bytes
        setstatus(stream, 400)
        startwrite(stream)
        return
    end
    # This upgrade is acceptable. Send the response.
    setheader(stream, "Sec-WebSocket-Accept" => generate_websocket_key(key))
    setheader(stream, "Upgrade" => "websocket")
    setheader(stream, "Connection" => "Upgrade")
    setstatus(stream, 101)
    startwrite(stream)
    # Pass the connection on as a WebSocket.
    io = getrawstream(stream)
    ws = WebSocket(io, true)
    # If the callback function f has two methods,
    # prefer the more secure one which takes (request, websocket)
    try
        if applicable(f, stream.message, ws)
            f(stream.message, ws)
        else
            f(ws)
        end
    catch err
        # Some errors will not reliably propagate when rethrown,
        # especially compile time errors.
        # On the server side, this function is running in a new task for every connection made
        # from outside. The rethrown errors might get lost or caught elsewhere, so we also
        # duplicate them to stderr here.
        # For working examples of error catching and reading them on the .out channel, see 'error_test.jl'.
        # If for some reason, the error messages from your 'f' cannot be read properly, here are
        # three alternative ways of finding them so you can correct:
        # 1) Include try..catch in your 'f', and print the errors to stderr.
        # 2) Turn the connection direction around, i.e. try to
        # provoke the error on the client side.
        # 3) Connect through a browser if that is not already what you are doing.
        # Some error messages may currently be shown there.
        # 4) use keyword argument loglevel = 3.
        # 5) modify the global logger to take control.
#        @warn("WebSockets.upgrade: Caught unhandled error while calling argument function f, the handler / gatekeeper:\n\t")
#        mt = typeof(f).name.mt
#        fnam = splitdir(string(mt.defs.func.file))[2]
#        printstyled(stderr, color= :yellow,"f = ", string(f) * " at " * fnam * ":" * string(mt.defs.func.line) * "\nERROR:\t")
#        showerror(stderr, err, stacktrace(catch_backtrace()))
         rethrow(err)
    finally
        close(ws)
    end
end

"""
Not normally called from user code.
Throws WebSocketError if the upgrade message is not basically valid.
Called from 'upgrade' for potential server side websockets,
and from `_openstream' for potential client side websockets.
"""
function check_upgrade(r)
    if !hasheader(r, "Upgrade", "websocket")
        throw(WebSocketError(0, "Check upgrade: Expected \"Upgrade => websocket\"!\n$(r)"))
    end
    if !(hasheader(r, "Connection", "upgrade") || hasheader(r, "Connection", "keep-alive, upgrade"))
        throw(WebSocketError(0, "Check upgrade: Expected \"Connection => upgrade or Connection => keep alive, upgrade\"!\n$(r)"))
    end
end

"""
Not normally called from user code.
Fast checking for websocket upgrade request vs content requests.
Called on all new connections in '_servercoroutine'.
"""
function is_upgrade(r::Request)
    if (r isa Request && r.method == "GET")  || (r isa Response && r.status == 101)
        if header(r, "Connection", "") != "keep-alive"
            # "Connection => upgrade" for most and "Connection => keep-alive, upgrade" for Firefox.
            if hasheader(r, "Connection", "upgrade") || hasheader(r, "Connection", "keep-alive, upgrade")
                if lowercase(header(r, "Upgrade", "")) == "websocket"
                    return true
                end
            end
        end
    end
    return false
end

is_upgrade(stream::Stream) = is_upgrade(stream.message)

# Inline docs in 'WebSockets.jl'
target(req::Request) = req.target
subprotocol(req::Request) = header(req, "Sec-WebSocket-Protocol")
origin(req::Request) = header(req, "Origin")

"""
WSHandlerFunction(f::Function) <: Handler
 The provided argument should be one of the forms

    `f(WebSocket) => nothing`
    `f(Request, WebSocket) => nothing`

    The latter form is intended for gatekeeping, ref. RFC 6455 section 10.1
 f accepts a `WebSocket` and does interesting things with it, like reading, writing and exiting when finished
 and ready for starting the close handshake.
"""
struct WSHandlerFunction{F <: Function} <: Handler
    func::F # func(ws) or func(request, ws)
end

"""
    WebSockets.ServerWS(handler::Function, wshandler::Function, logger::IO)

WebSockets.ServerWS is an argument type for WebSockets.serve. Instances
include .in  and .out channels, see WebSockets.serve.

Special server options can be set using keyword arguments, and will overrule default_options().
"""
struct ServerWS
    handler::RequestHandlerFunction
    wshandler::WSHandlerFunction
    in::Channel{Any}
    out::Channel{Any}
    # the remaining parameters are server options,
    # documented in HTTP.jl
    sslconfig::Union{SSLConfig, Nothing}
    tcpisvalid::Function
    reuseaddr::Bool
    connection_count::Ref{Int}
    rate_limit::Union{Rational{Int}, Nothing}
    reuse_limit::Int
    readtimeout::Int
    # keywords only constructor
    ServerWS(;handler, wshandler, in, out, sslconfig, tcpisvalid, reuseaddr, connection_count, rate_limit, reuse_limit, readtimeout)=
          new(handler, wshandler, in, out, sslconfig, tcpisvalid, reuseaddr, connection_count, rate_limit, reuse_limit, readtimeout)
end
# User supplied keywords take precedence over default_options() here:
function ServerWS(handler::RequestHandlerFunction, wshandler::WSHandlerFunction; kwargs...)
    ServerWS(;handler = handler, wshandler = wshandler,  merge(default_options(), collect(kwargs))...)
end
# User supplied handler functions are packed into specialized types here.
# It also defines the most useful method ServerWS(h::Function,w::Function)
function ServerWS(h::Function, w::Function; kwargs...)
    ServerWS(RequestHandlerFunction(h), WSHandlerFunction(w); kwargs...)
end

to_ipaddr(ip::IPAddr) = ip
to_ipaddr(s::AbstractString) = parse(IPAddr, s)

"""
    WebSockets.serve(server::ServerWS, port)
    WebSockets.serve(server::ServerWS, host, port)
    WebSockets.serve(server::ServerWS, host, port, verbose)

A higher level interface for WebSockets.HTTP.listen.
If you are running serve asyncronously, the serve task will
exit when you close(ServerWS).

Puts any caught error and stacktrace on the server.out channel.
The server.in channel is used by close(ServerWS).

# Example
```julia
    @async WebSockets.serve(myserver_WS, "127.0.0.1", 8080)
```
After a suspected connection task failure:
```julia
    if isready(myserver_WS.out)
        stack_trace = take!(myserver_WS.out)
    end
```
"""
function serve(serverws::ServerWS, host, port, verbose)
    # An internal reference used for closing.
    tcpserver = Ref{Union{IOServer, Nothing}}(Sockets.listen(InetAddr(to_ipaddr(host), port)))
    # Start a couroutine that sleeps until something is put on the .in channel.
    @async begin
         # Next line will hold until something is put on the channel
         take!(serverws.in)
         if tcpserver[] != nothing
             close(tcpserver[])
             tcpserver[] = nothing
               GC.gc()
             yield()
        end
    end
    # We capture some variables in this inner function, which takes just one-argument.
    # The inner function will be called in a new task for every incoming connection.
    function _servercoroutine(stream::Stream)
        try
            if is_upgrade(stream.message)
                upgrade(serverws.wshandler.func, stream)
            else
                handle(serverws.handler, stream)
            end
        catch err
            put!(serverws.out, err)
            put!(serverws.out, stacktrace(catch_backtrace()))
        end
    end
    #
    # Call the constructor for the listen loop. When somebody tries to open a connection,
    # the listenloop
    # 1) Checks if we are willing to accept
    # 2) Spawns a new task, namely our coroutine _servercoroutine.
    HTTP.listen(_servercoroutine,
            host, port,
            verbose = verbose,
            sslconfig = serverws.sslconfig,
            tcpisvalid = serverws.tcpisvalid,
            server = tcpserver[],
            reuseaddr = serverws.reuseaddr,
            connection_count = serverws.connection_count,
            rate_limit = serverws.rate_limit,
            reuse_limit = serverws.reuse_limit,
            readtimeout = serverws.readtimeout)
    # We will only get to this point if the server is closed.
    # If this serve function is running as a coroutine, the server is closed
    # through the server.in channel, see above.
    return
end
serve(serverws::ServerWS; host= "127.0.0.1", port= "") =  serve(serverws, host, port, false)
serve(serverws::ServerWS, host, port) =  serve(serverws, host, port, false)
serve(serverws::ServerWS, port) =  serve(serverws, "127.0.0.1", port, false)

function Base.close(serverws::ServerWS)
    put!(serverws.in, "Close!")
    if isready(serverws.out)
        @warn "Server closed. Error messages exist on the server's .out channel. Check
            string(take!(yourserverWS.out))"
    end
    return
end


text:
# Before running this, please make sure to activate and instantiate the environment
# corresponding to [this `Project.toml`](https://raw.githubusercontent.com/alan-turing-institute/DataScienceTutorials.jl/master/Project.toml) and [this `Manifest.toml`](https://raw.githubusercontent.com/alan-turing-institute/DataScienceTutorials.jl/master/Manifest.toml)
# so that you get an environment which matches the one used to generate the tutorials:
#
# ```julia
# cd("DataScienceTutorials") # cd to folder with the *.toml
# using Pkg; Pkg.activate("."); Pkg.instantiate()
# ```

using CategoricalArrays

v = categorical(["AA", "BB", "CC", "AA", "BB", "CC"])

levels(v)

v = categorical([1, 2, 3, 1, 2, 3, 1, 2, 3], ordered=true)

levels(v)

v[1] < v[2]

v = categorical(["high", "med", "low", "high", "med", "low"], ordered=true)

levels(v)

v[1] < v[2]

levels!(v, ["low", "med", "high"])

v[1] < v[2]

v = categorical(["AA", "BB", missing, "AA", "BB", "CC"]);

levels(v)

# This file was generated using Literate.jl, https://github.com/fredrikekre/Literate.jl

ext: jl


text:
PoincareMap = ChaosTools.PoincareMap

mutable struct BasinInfo{I,F,D,Q}
    basin :: I
    xg :: F
    yg :: F
    iter_f! :: Function
    reinit_f! :: Function
    get_u :: Function
    current_color :: Int64
    next_avail_color :: Int64
    consecutive_match :: Int64
    consecutive_other_basins :: Int64
    prevConsecutives :: Int64
    prev_attr :: Int64
    prev_bas :: Int64
    prev_step :: Int64
    step :: Int64
    attractors :: D #Dict{Int16,Vector{Vector{Float64}}}
    Na :: Int64
    visited :: Q
    available :: Q
end

function Base.show(io::IO, bsn::BasinInfo)
    println(io, "Basin of attraction structure")
    println(io,  rpad(" size : ", 14),    size(bsn.basin))
    println(io,  rpad(" Number of attractors found: ", 14),   bsn.Na  )
end

"""
    basins_map2D(xg, yg, integ; kwargs...)
Compute an estimate of the basin of attraction on a two-dimensional plane using a map of the plane onto itself.
The dynamical system should be a discrete two dimensional system such as:
    * Discrete 2D map.
    * 2D poincaré map.
    * A 2D stroboscopic map.

[H. E. Nusse and J. A. Yorke, Dynamics: numerical explorations, Springer, New York, 1997]

## Arguments
* `xg`, `yg` : 1-dim range vector that defines the grid of the initial conditions to test.
* `integ` : A  integrator handle of the dynamical system. For a Poincaré map, the handle is a `pmap`
as defined in [`poincaremap`](@ref)

## Keyword Arguments
* `T` : Period of the stroboscopic map.
* `Ncheck` : A parameter that sets the number of consecutives hits of an attractor before deciding the basin
of the initial condition.

## Example:

```jl
using DynamicalSystems, ChaosTools
ds = Systems.rikitake(μ = 0.47, α = 1.0)
xg=range(-6.,6.,length=150); yg=range(-6.,6.,length=150)
pmap = poincaremap(ds, (3, 0.), Tmax=1e6; idxs = 1:2, rootkw = (xrtol = 1e-8, atol = 1e-8), reltol=1e-9)
bsn = basins_map2D(xg, yg, pmap)
```

"""
function basins_map2D(xg, yg, pmap::PoincareMap; Ncheck = 3)
    reinit_f! = (pmap,y) -> _init_map(pmap, y, pmap.i)
    get_u = (pmap) -> pmap.integ.u[pmap.i]
    basin = draw_basin!(xg, yg, pmap, step!, reinit_f!, get_u, Ncheck)
end

function _init_map(pmap::PoincareMap, y, idxs)
    u = zeros(1,length(pmap.integ.u))
    u[idxs] = y
    # all other coordinates are zero
    reinit!(pmap, u)
end
# nueva linea

function basins_map2D(xg, yg, integ; T=0., Ncheck = 2)
    if T>0
        iter_f! = (integ) -> step!(integ, T, true)
    else
        iter_f! = (integ) -> step!(integ)
    end
    reinit_f! =  (integ,y) -> reinit!(integ, y)
    get_u = (integ) -> integ.u

    return draw_basin!(xg, yg, integ, iter_f!, reinit_f!, get_u, Ncheck)
end

"""
    basins_general(xg, yg, integ; T=1., idxs=1:2)
Compute an estimate of the basin of attraction on a two-dimensional plane using a stroboscopic map.

## Arguments
* `xg`, `yg` : 1-dim range vector that defines the grid of the initial conditions to test.
* `integ` : integrator handle of the dynamical system.

## Keyword Arguments
* `dt` : Time step of the integrator. It recommended to use values above 1.
* `idxs = 1:D` : Optionally you can choose which variables to save. Defaults to the entire state.
* `Ncheck` : A parameter that sets the number of consecutives hits of an attractor before deciding the basin
of the initial condition.

## Example:
```jl
using DynamicalSystems, ChaosTools
ds = Systems.magnetic_pendulum(γ=1, d=0.2, α=0.2, ω=0.8, N=3)
integ = integrator(ds, u0=[0,0,0,0], reltol=1e-9)
xg=range(-4,4,length=150)
yg=range(-4,4,length=150)
@time bsn = basins_general(xg, yg, integ; dt=1., idxs=1:2)
"""
function basins_general(xg, yg, integ; dt=1., idxs=1:2, Ncheck = 10)
    i = typeof(idxs) <: Int ? i : SVector{length(idxs), Int}(idxs...)
    iter_f! = (integ) -> step!(integ, dt, true)
    reinit_f! =  (integ,y) -> _init_ds(integ, y, i)
    get_u = (integ) -> integ.u[i]
    return draw_basin!(xg, yg, integ, iter_f!, reinit_f!,get_u, Ncheck)
end

function _init_ds(integ, y, idxs)
    u = zeros(length(integ.u))
    u[idxs] = y
    # all other coordinates are zero
    reinit!(integ, u)
end

## Procedure described in  H. E. Nusse and J. A. Yorke, Dynamics: numerical explorations, Springer, New York, 1997
# The idea is to color the grid with the current color. When an attractor box is hit (even color), the initial condition is colored
# with the color of its basin (odd color). If the trajectory hits another basin 10 times in row the IC is colored with the same
# color as this basin.
function procedure!(bsn_nfo::BasinInfo, n, u, Ncheck::Int)
    max_check = 60
    #next_c = bsn_nfo.basin[n]
    next_c=get_data(bsn_nfo.basin, n)
    bsn_nfo.step += 1

    if iseven(next_c) && bsn_nfo.consecutive_match < max_check
        # check wether or not we hit an attractor (even color). Make sure we hit Ncheck consecutive times.
        if bsn_nfo.prev_attr == next_c
            bsn_nfo.prevConsecutives +=1
        else
            # reset prevConsecutives
            bsn_nfo.prev_attr = next_c
            bsn_nfo.prevConsecutives =1
            return 0;
        end

        if bsn_nfo.prevConsecutives >= Ncheck
            # Wait if we hit the attractor a Ncheck times in a row just to check if it is not a nearby trajectory
            c3 = next_c+1
            if Ncheck == 2
                # For maps we can color the previous steps as well. Every point of the trajectory lead
                # to the attractor
                find_and_replace!(bsn_nfo, bsn_nfo.current_color+1, c3)
            else
                # For higher dimensions we erase the past iterations.
                find_and_replace!(bsn_nfo, bsn_nfo.current_color+1, 1)
            end
            reset_bsn_nfo!(bsn_nfo)
            return c3
         end
    end

    if next_c == 1 && bsn_nfo.consecutive_match < max_check
        # uncolored box, color it with current odd color
        set_data!(bsn_nfo.basin, n, bsn_nfo.current_color+1)
        push!(bsn_nfo.visited,n)
        bsn_nfo.consecutive_match = 0
        return 0
    elseif next_c == 1 && bsn_nfo.consecutive_match >= max_check
        # Maybe chaotic attractor, perodic or long recursion.
        # Color this box as part of an attractor
        set_data!(bsn_nfo.basin, n, bsn_nfo.current_color)
        store_attractor!(bsn_nfo, u)
        bsn_nfo.consecutive_match = max_check
        return 0
    elseif next_c == bsn_nfo.current_color + 1
        # hit a previously visited box with the current color, possible attractor?
        if bsn_nfo.consecutive_match < max_check
            bsn_nfo.consecutive_match += 1
            return 0
        else
            set_data!(bsn_nfo.basin, n, bsn_nfo.current_color)
            store_attractor!(bsn_nfo, u)
            # We continue iterating until we hit again the same attractor. In which case we stop.
            return 0;
        end
    elseif isodd(next_c) && 0 < next_c < bsn_nfo.current_color &&  bsn_nfo.consecutive_match < max_check && Ncheck == 2
        # hit a colored basin point of the wrong basin, happens all the time, we check if it happens
        #10 times in a row or if it happens N times along the trajectory whether to decide if it is another basin.
        bsn_nfo.consecutive_other_basins += 1

        if bsn_nfo.prev_bas == next_c &&  bsn_nfo.prev_step == bsn_nfo.step-1
            bsn_nfo.prevConsecutives +=1
            bsn_nfo.prev_step += 1
        else
            bsn_nfo.prev_bas = next_c
            bsn_nfo.prev_step = bsn_nfo.step
            bsn_nfo.prevConsecutives =1
        end

        if bsn_nfo.consecutive_other_basins > 60 || bsn_nfo.prevConsecutives > 10
            find_and_replace!(bsn_nfo, bsn_nfo.current_color+1, next_c)
            reset_bsn_nfo!(bsn_nfo)
            return next_c
        end
        return 0
    elseif iseven(next_c) &&   (max_check <= bsn_nfo.consecutive_match < 2*max_check)
        # We make sure we hit the attractor 60 consecutive times
        bsn_nfo.consecutive_match+=1
        return 0
    elseif iseven(next_c) && bsn_nfo.consecutive_match >= max_check*2
        # We have checked the presence of an attractor: tidy up everything and get a new box.
        find_and_replace!(bsn_nfo, bsn_nfo.current_color+1, 1)
        # store color
        set_data!(bsn_nfo.basin, n, bsn_nfo.current_color)
        store_attractor!(bsn_nfo, u) # store attractor
        bsn_nfo.current_color = bsn_nfo.next_avail_color
        bsn_nfo.next_avail_color += 2
        println("even and > max check ", next_c)
        reset_bsn_nfo!(bsn_nfo)
        return next_c+1;
    else
        return 0
    end
end

function store_attractor!(bsn_nfo::BasinInfo, u)
    if haskey(bsn_nfo.attractors , bsn_nfo.current_color)
        push!(bsn_nfo.attractors[bsn_nfo.current_color],  u) # store attractor
    else
        bsn_nfo.attractors[bsn_nfo.current_color] = Dataset([SVector(u[1],u[2])])  # init dic
    end
end

"""
    draw_basin!(xg, yg, integ, iter_f!::Function, reinit_f!::Function)
Compute an estimate of the basin of attraction on a two-dimensional plane. This is a low level function,
for higher level functions see: `basin_poincare_map`, `basin_discrete_map`, `basin_stroboscopic_map`

## Arguments
* `xg`, `yg` : 1-dim range vector that defines the grid of the initial conditions to test.
* `integ` : integrator handle of the dynamical system.
* `iter_f!` : function that iterates the map or the system, see step! from DifferentialEquations.jl and
examples for a Poincaré map of a continuous system.
* `reinit_f!` : function that sets the initial condition to test on a two dimensional projection of the phase space.
"""
function draw_basin!(xg, yg, integ, iter_f!::Function, reinit_f!::Function, get_u::Function, Ncheck)

    complete = 0;

    bsn_nfo = BasinInfo(ones(Int16, length(xg), length(yg)), xg, yg, iter_f!, reinit_f!, get_u, 2,4,0,0,0,1,1,0,0,
                    Dict{Int16,Dataset{2,Float64}}(),0,Vector{CartesianIndex}(),Vector{CartesianIndex}())

    reset_bsn_nfo!(bsn_nfo)

    I = CartesianIndices(bsn_nfo.basin)
    j  = 1

    while complete == 0
        # pick the first empty box
        ind = 0
        for k in j:length(bsn_nfo.basin)
            if bsn_nfo.basin[I[k]] == 1
                j = k
                ind=I[k]
                break
            end
        end

        if ind == 0
            # We are done
            complete = 1
            break
        end

         ni = ind[1]; mi = ind[2];
         x0 = xg[ni]; y0 = yg[mi]

         # Tentatively assign a color: odd is for basins, even for attractors.
         # First color is one
         bsn_nfo.basin[ni,mi] = bsn_nfo.current_color + 1
         u0=[x0, y0]
         bsn_nfo.basin[ni,mi] = get_color_point!(bsn_nfo, integ, u0; Ncheck=Ncheck)
    end

    bsn_nfo.Na = Int((bsn_nfo.current_color-2)/2)

    return bsn_nfo
end

function get_color_point!(bsn_nfo::BasinInfo, integ, u0; Ncheck=2)
    # This routine identifies the attractor using the previously defined basin.
    # reinitialize integrator
    bsn_nfo.reinit_f!(integ, u0)
    reset_bsn_nfo!(bsn_nfo)

    done = 0;
    inlimbo = 0

    while done == 0
       old_u = bsn_nfo.get_u(integ)
       bsn_nfo.iter_f!(integ)
       new_u = bsn_nfo.get_u(integ)

       n = get_box(new_u, bsn_nfo)

       if !isnothing(n) # apply procedure only for boxes in the defined space
           done = procedure!(bsn_nfo, n, new_u, Ncheck)
           inlimbo = 0
       else
           # We are outside the defined grid
           inlimbo +=1
       end

       if inlimbo > 60
           done = check_outside_the_screen!(bsn_nfo, new_u, old_u, inlimbo)
       end
    end

    return done
end

function get_box(u, bsn_nfo::BasinInfo)
    xg = bsn_nfo.xg; yg = bsn_nfo.yg; # aliases
    xstep = (xg[2]-xg[1])
    ystep = (yg[2]-yg[1])
    xu=u[1]
    yu=u[2]
    if typeof(bsn_nfo.basin) <: Cell
        if xu >= xg[1] && xu <= xg[end] && yu >= yg[1] && yu <= yg[end]
            return findleaf(bsn_nfo.basin,u)
        else
            return nothing
        end
    end
    n=0; m=0;
    # check boundary
    if xu >= xg[1] && xu <= xg[end] && yu >= yg[1] && yu <= yg[end]
        n = Int(round((xu-xg[1])/xstep))+1
        m = Int(round((yu-yg[1])/ystep))+1 # +1 for 1 based indexing
    else
        return nothing
    end
    return CartesianIndex(n,m)
end

function check_outside_the_screen!(bsn_nfo::BasinInfo, new_u, old_u, inlimbo)

    if norm(new_u-old_u) < 1e-5
        find_and_replace!(bsn_nfo,bsn_nfo.current_color+1, 1)
        reset_bsn_nfo!(bsn_nfo)
        # this CI goes to a attractor outside the screen, set to -1 (even color)
        return -1  # get next box
    elseif inlimbo > 60*20
        find_and_replace!(bsn_nfo,bsn_nfo.current_color+1, 1)
        reset_bsn_nfo!(bsn_nfo)
        # this CI is problematic or diverges, set to -1 (even color)
        return -1  # get next box
    end
    return 0
end

function find_and_replace!(bsn_nfo::BasinInfo, old_c, new_c)
    if typeof(bsn_nfo.basin) <: Matrix
        while !isempty(bsn_nfo.visited)
           ind = pop!(bsn_nfo.visited)
           if bsn_nfo.basin[ind] == old_c
               bsn_nfo.basin[ind] = new_c
           end
        end
        # ind = (bsn_nfo.basin .== old_c)
        # bsn_nfo.basin[ind] .= new_c
    else
        while !isempty(bsn_nfo.visited)
           lf = pop!(bsn_nfo.visited)
           if lf.data == old_c
               lf.data = new_c
           end
        end
        # for lf in allleaves(bsn_nfo.basin)
        #     if lf.data==old_c
        #         lf.data==new_c
        #     end
        # end

    end
end

function reset_bsn_nfo!(bsn_nfo::BasinInfo)
    #@show bsn_nfo.step
    bsn_nfo.consecutive_match = 0
    bsn_nfo.consecutive_other_basins = 0
    bsn_nfo.prevConsecutives = 0
    bsn_nfo.prev_attr = 1
    bsn_nfo.prev_bas = 1
    bsn_nfo.prev_step = 0
    bsn_nfo.step = 0
end

function basins_map2D_tree(xg, yg, integ; T=0., Ncheck = 2, r_init=0.2, r_max=0.1)
    if T>0
        iter_f! = (integ) -> step!(integ, T, true)
    else
        iter_f! = (integ) -> step!(integ)
    end
    reinit_f! =  (integ,y) -> reinit!(integ, y)
    get_u = (integ) -> integ.u

    return draw_basin_tree!(xg, yg, integ, iter_f!, reinit_f!, get_u, Ncheck, r_init, r_max)
end

import RegionTrees: AbstractRefinery, needs_refinement, refine_data

struct MyRefinery <: AbstractRefinery
    tolerance::Float64
end

# These two methods are all we need to implement
function needs_refinement(r::MyRefinery, cell)
        return maximum(cell.boundary.widths) > r.tolerance
end

function refine_data(r::MyRefinery, cell::Cell, indices)
        return 1
end

function is_in_cell(cell::Cell,  point::AbstractVector, indices)
    or = cell.boundary.origin
    wd = cell.boundary.widths
    dv = cell.divisions
    if point[1] > or[1] + wd[1] || point[1] < or[1]
        return false
    elseif point[2] > or[2] + wd[2] || point[2] < or[2]
        return false
    else
        ind = [1,1]
        if point[1] > dv[1]
            ind[1]=2
        end
        if point[2] > dv[2]
            ind[2]=2
        end
        if ind[1] == indices[1] && ind[2] == indices[2]
            return true
        else
            return false
        end
    end
end

function refine_data(cell::Cell, indices, bsn_nfo::BasinInfo)
    if iseven(cell.data)
        att = cell.data
        for p in bsn_nfo.attractors[att]
            if is_in_cell(cell,p, indices)
                #@show att,p,cell.boundary,indices
                return att
            end
        end
    end
    return 1
end

function get_neighbor_cells(root::Cell,cell::Cell)
    or = cell.boundary.origin
    wd = cell.boundary.widths
    dv = cell.divisions
    W=findleaf(root,dv + SVector(wd[1],0)).data
    E=findleaf(root,dv + SVector(-wd[1],0)).data
    N=findleaf(root,dv + SVector(0,wd[2])).data
    S=findleaf(root,dv + SVector(0,-wd[2])).data
    #@show [W E N S cell.data]
    return length(unique([W E N S cell.data])) > 1
end

function check_if_complete!(root::Cell, r_max, bsn_nfo::BasinInfo)
    complete = true
    refine_function = (cell, indices) -> refine_data(cell, indices, bsn_nfo)
    leaf_stack = Vector{Cell}()
    for lf in allleaves(root)
        if maximum(lf.boundary.widths) > r_max && get_neighbor_cells(root,lf)
            complete = false
            push!(leaf_stack,lf)
        end
    end

    while !isempty(leaf_stack)
        l = pop!(leaf_stack)
        split!(l, refine_function)
        for cl in children(l)
            if cl.data == 1
                push!(bsn_nfo.available,cl)
            end
        end
    end

    return complete
end

function draw_basin_tree!(xg, yg, integ, iter_f!::Function, reinit_f!::Function, get_u::Function, Ncheck, r_init, r_max)

    complete = 0;

    xi=xg[1]; yi=yg[1]; xf=xg[end]; yf=yg[end]

    bsn_nfo = BasinInfo(Cell(SVector(xi, yi), SVector(xf-xi, yf-yi),1), xg, yg, iter_f!, reinit_f!, get_u, 2,4,0,0,0,1,1,0,0,
                Dict{Int16,Dataset{2,Float64}}(),0,Vector{Cell}(),Vector{Cell}())

    reset_bsn_nfo!(bsn_nfo)

    # Initial refinement of the grid.
    r = MyRefinery(r_init)
    @show bsn_nfo.basin
    adaptivesampling!(bsn_nfo.basin, r)

    for lf in allleaves(bsn_nfo.basin)
        push!(bsn_nfo.available,lf)
    end

    while complete == 0

        # Get next available candidate
        lf = nothing
        while !isempty(bsn_nfo.available)
            l = pop!(bsn_nfo.available)
            if l.data == 1
                lf = l
                break
            end
        end

        if isempty(bsn_nfo.available)
            println("refinement!")
            if check_if_complete!(bsn_nfo.basin, r_max, bsn_nfo)
                complete=1
                break
            end
            # Get next available candidate
            while !isempty(bsn_nfo.available)
                l = pop!(bsn_nfo.available)
                if l.data == 1
                    lf = l
                    break
                end
            end

            println("end refinement! res=", lf.boundary.widths)

        end

        # Tentatively assign a color: odd is for basins, even for attractors.
        # First color is one
        lf.data = bsn_nfo.current_color + 1
        u0 = RegionTrees.center(lf)
        lf.data = get_color_point!(bsn_nfo, integ, u0; Ncheck=Ncheck)

    end
    bsn_nfo.Na = (bsn_nfo.current_color-2)÷2
    return bsn_nfo
end

function get_data(bsn::Matrix,n)
    return bsn[n]
end

function get_data(bsn::Cell,n)
    return n.data
end

function set_data!(bsn::Matrix,n,c)
    bsn[n]=c
end

function set_data!(bsn::Cell,n,c)
    n.data=c
end


text:
mutable struct SwitchAssessState{T}
    weight::Float64
    retval::T
    SwitchAssessState{T}(weight::Float64) where T = new{T}(weight)
end

function process!(gen_fn::Switch{C, N, K, T},
                  index::Int,
                  args::Tuple,
                  choices::ChoiceMap,
                  state::SwitchAssessState{T}) where {C, N, K, T}
    (weight, retval) = assess(getindex(gen_fn.branches, index), args, choices)
    state.weight = weight
    state.retval = retval
end

@inline process!(gen_fn::Switch{C, N, K, T}, index::C, args::Tuple, choices::ChoiceMap, state::SwitchAssessState{T}) where {C, N, K, T} = process!(gen_fn, getindex(gen_fn.cases, index), args, choices, state)

function assess(gen_fn::Switch{C, N, K, T},
                args::Tuple,
                choices::ChoiceMap) where {C, N, K, T}
    index = args[1]
    state = SwitchAssessState{T}(0.0)
    process!(gen_fn, index, args[2 : end], choices, state)
    return state.weight, state.retval
end

max_stars_repo_name: georgematheos/Gen

max_issues_repo_name: georgematheos/Gen

max_forks_repo_name: georgematheos/Gen


text:
mutable struct ResultChannel <: AbstractChannel{Result}
    current_result::Result
    ResultChannel(result::Result) = new(result)
end

function put!(channel::ResultChannel, result::Result)
    if result.cost_minimum < channel.current_result.cost_minimum
        channel.current_result = result
    end
    channel
end

function take!(channel::ResultChannel)
    fetch(channel)
end

function fetch(channel::ResultChannel)
    channel.current_result
end


text:
struct MyType
    a::Int
    b::Float64
end

@testset "TreeArray" begin
    for NodeType in (Node{Node{LeafNode{Float64, 2, 2}, 2, 2}, 2, 2},
                     HashNode{Node{LeafNode{Float64, 2, 2}, 2, 2}, 2, 2},
                     DynamicNode{Node{LeafNode{Float64, 2, 2}, 2, 2}, 2},
                     DynamicHashNode{Node{LeafNode{Float64, 2, 2}, 2, 2}, 2})
        A = @inferred TreeArray(NodeType, 16, 16)

        @test size(A) == (16, 16)
        if NodeType <: Union{Node, HashNode}
            @test size(rootnode(A.tree)) == (4, 4)
        else
            @test size(rootnode(A.tree)) == (1, 1)
        end

        A[1:3, 1:3] .= reshape(1:9, 3, 3)
        @test A[1:3, 1:3] == reshape(1:9, 3, 3)

        A[11:13, 11:13] .= reshape(11:19, 3, 3)
        @test A[11:13, 11:13] == reshape(11:19, 3, 3)

        # check activation
        mask = map(eachindex(A)) do i
            isactive(A, Tuple(i)...)
        end
        @test all(mask[1:3, 1:3])
        @test all(mask[11:13, 11:13])
        @test !any(mask[1:3, 4:10])
        @test !any(mask[11:13, 1:3])

        # allocate!
        A.tree .= nothing
        mask = similar(A, Bool)
        TreeArrays.allocate!(A, mask)
        @test map(i -> isactive(A, i), eachindex(A)) == mask
    end
end

@testset "StructTreeArray" begin
    for NodeType in (Node{Node{@StructLeafNode{MyType, 2, 2}, 2, 2}, 2, 2},
                     HashNode{Node{@StructLeafNode{MyType, 2, 2}, 2, 2}, 2, 2},
                     DynamicNode{Node{@StructLeafNode{MyType, 2, 2}, 2, 2}, 2},
                     DynamicHashNode{Node{@StructLeafNode{MyType, 2, 2}, 2, 2}, 2})
        A = @inferred StructTreeArray(NodeType, 16, 16)

        # allocate!
        mask = rand(Bool, size(A))
        mask[1] = true # make it true in at least one element
        TreeArrays.allocate!(A, mask)
        @test map(i -> isactive(A, i), eachindex(A)) == mask

        # isallocated
        # @test map(i -> TreeArrays.isallocated(A, i), eachindex(A)) == mask
        A .= nothing
        # @test map(i -> TreeArrays.isallocated(A, i), eachindex(A)) == mask
        @test any(map(i -> TreeArrays.isallocated(A, i), eachindex(A)))
        TreeArrays.cleanup!(A)
        @test map(i -> TreeArrays.isallocated(A, i), eachindex(A)) == falses(16, 16)

        for i in 1:length(A)
            A.a[i] = i
            A.b[i] = 2i
        end
        @test A.a == reshape(1:16*16, 16, 16)
        @test A.b == 2*reshape(1:16*16, 16, 16)

        @test continuousview(A, 3:11, 9:14) == OffsetArray(A[3:11, 9:14], 3:11, 9:14)
        @test continuousview(A, 3:11, 9:14).a == OffsetArray(A.a[3:11, 9:14], 3:11, 9:14)
    end
end


text:
module ThreadedSparseArrays

export ThreadedSparseMatrixCSC

using LinearAlgebra
import LinearAlgebra: mul!
using SparseArrays
import SparseArrays: getcolptr, AbstractSparseMatrixCSC
const AdjOrTransDenseMatrix = if VERSION < v"1.6.0-rc2"
    SparseArrays.AdjOrTransStridedOrTriangularMatrix
else
    SparseArrays.AdjOrTransDenseMatrix
end

# * Threading utilities
struct RangeIterator
    k::Int
    d::Int
    r::Int
end

"""
    RangeIterator(n::Int,k::Int)

Returns an iterator splitting the range `1:n` into `min(k,n)` parts of (almost) equal size.
"""
RangeIterator(n::Int, k::Int) = RangeIterator(min(n,k),divrem(n,k)...)
Base.length(it::RangeIterator) = it.k
endpos(it::RangeIterator, i::Int) = i*it.d+min(i,it.r)
Base.iterate(it::RangeIterator, i::Int=1) = i>it.k ? nothing : (endpos(it,i-1)+1:endpos(it,i), i+1)

# * ThreadedSparseMatrixCSC

"""
    ThreadedSparseMatrixCSC(A)

Thin container around `A::SparseMatrixCSC` that will enable certain
threaded multiplications of `A` with dense matrices.
"""
struct ThreadedSparseMatrixCSC{Tv,Ti,At} <: AbstractSparseMatrixCSC{Tv,Ti}
    A::At
    ThreadedSparseMatrixCSC(A::At) where {Tv,Ti,At<:AbstractSparseMatrixCSC{Tv,Ti}} =
        new{Tv,Ti,At}(A)
end

Base.size(A::ThreadedSparseMatrixCSC, args...) = size(A.A, args...)

for f in [:rowvals, :nonzeros, :getcolptr]
    @eval SparseArrays.$(f)(A::ThreadedSparseMatrixCSC) = SparseArrays.$(f)(A.A)
end

# sparse * sparse multiplications are not (currently) threaded, but we want to keep the return type
for (T1,t1) in ((ThreadedSparseMatrixCSC,identity), (Adjoint{<:Any,<:ThreadedSparseMatrixCSC},adjoint), (Transpose{<:Any,<:ThreadedSparseMatrixCSC},transpose))
    for (T2,t2) in ((ThreadedSparseMatrixCSC,identity), (Adjoint{<:Any,<:ThreadedSparseMatrixCSC},adjoint), (Transpose{<:Any,<:ThreadedSparseMatrixCSC},transpose))
        @eval Base.:(*)(A::$T1, B::$T2) = ThreadedSparseMatrixCSC($t1($t1(A).A)*$t2($t2(B).A))
    end
end

function mul!(C::StridedVecOrMat, A::ThreadedSparseMatrixCSC, B::Union{StridedVector,AdjOrTransDenseMatrix}, α::Number, β::Number)
    size(A, 2) == size(B, 1) || throw(DimensionMismatch())
    size(A, 1) == size(C, 1) || throw(DimensionMismatch())
    size(B, 2) == size(C, 2) || throw(DimensionMismatch())
    nzv = nonzeros(A)
    rv = rowvals(A)
    if β != 1
        β != 0 ? rmul!(C, β) : fill!(C, zero(eltype(C)))
    end
    @sync for r in RangeIterator(size(C,2), Threads.nthreads())
        Threads.@spawn for k in r
            @inbounds for col in 1:size(A, 2)
                αxj = B[col,k] * α
                for j in nzrange(A, col)
                    C[rv[j], k] += nzv[j]*αxj
                end
            end
        end
    end
    C
end

for (T, t) in ((Adjoint, adjoint), (Transpose, transpose))
    @eval function mul!(C::StridedVecOrMat, xA::$T{<:Any,<:ThreadedSparseMatrixCSC}, B::AdjOrTransDenseMatrix, α::Number, β::Number)
        A = xA.parent
        size(A, 2) == size(C, 1) || throw(DimensionMismatch())
        size(A, 1) == size(B, 1) || throw(DimensionMismatch())
        size(B, 2) == size(C, 2) || throw(DimensionMismatch())
        nzv = nonzeros(A)
        rv = rowvals(A)
        if β != 1
            β != 0 ? rmul!(C, β) : fill!(C, zero(eltype(C)))
        end
        @sync for r in RangeIterator(size(C,2), Threads.nthreads())
            Threads.@spawn for k in r
                @inbounds for col in 1:size(A, 2)
                    tmp = zero(eltype(C))
                    for j in nzrange(A, col)
                        tmp += $t(nzv[j])*B[rv[j],k]
                    end
                    C[col,k] += tmp * α
                end
            end
        end
        C
    end

    @eval function mul!(C::StridedVecOrMat, xA::$T{<:Any,<:ThreadedSparseMatrixCSC}, B::StridedVector, α::Number, β::Number)
        A = xA.parent
        size(A, 2) == size(C, 1) || throw(DimensionMismatch())
        size(A, 1) == size(B, 1) || throw(DimensionMismatch())
        size(B, 2) == size(C, 2) || throw(DimensionMismatch())
        @assert size(B,2)==1
        nzv = nonzeros(A)
        rv = rowvals(A)
        if β != 1
            β != 0 ? rmul!(C, β) : fill!(C, zero(eltype(C)))
        end
        @sync for r in RangeIterator(size(A,2), Threads.nthreads())
            Threads.@spawn @inbounds for col in r
                tmp = zero(eltype(C))
                for j in nzrange(A, col)
                    tmp += $t(nzv[j])*B[rv[j]]
                end
                C[col] += tmp * α
            end
        end
        C
    end
end

function mul!(C::StridedVecOrMat, X::AdjOrTransDenseMatrix, A::ThreadedSparseMatrixCSC, α::Number, β::Number)
    mX, nX = size(X)
    nX == size(A, 1) || throw(DimensionMismatch())
    mX == size(C, 1) || throw(DimensionMismatch())
    size(A, 2) == size(C, 2) || throw(DimensionMismatch())
    rv = rowvals(A)
    nzv = nonzeros(A)
    if β != 1
        β != 0 ? rmul!(C, β) : fill!(C, zero(eltype(C)))
    end
    # TODO: split in X isa DenseMatrixUnion and X isa Adjoint/Transpose so we can use @simd in the first case (see original code in SparseArrays)
    @sync for r in RangeIterator(size(A,2), Threads.nthreads())
        Threads.@spawn for col in r
            @inbounds for k in nzrange(A, col)
                Aiα = nzv[k] * α
                rvk = rv[k]
                for multivec_row in 1:mX
                    C[multivec_row, col] += X[multivec_row, rvk] * Aiα
                end
            end
        end
    end
    C
end

for (T, t) in ((Adjoint, adjoint), (Transpose, transpose))
    @eval function mul!(C::StridedVecOrMat, X::AdjOrTransDenseMatrix, xA::$T{<:Any,<:ThreadedSparseMatrixCSC}, α::Number, β::Number)
        A = xA.parent
        mX, nX = size(X)
        nX == size(A, 2) || throw(DimensionMismatch())
        mX == size(C, 1) || throw(DimensionMismatch())
        size(A, 1) == size(C, 2) || throw(DimensionMismatch())
        rv = rowvals(A)
        nzv = nonzeros(A)
        if β != 1
            β != 0 ? rmul!(C, β) : fill!(C, zero(eltype(C)))
        end

        # transpose of Threaded * Dense algorithm above
        @sync for r in RangeIterator(size(C,1), Threads.nthreads())
            Threads.@spawn for k in r
                @inbounds for col in 1:size(A, 2)
                    αxj = X[k,col] * α
                    for j in nzrange(A, col)
                        C[k, rv[j]] += $t(nzv[j])*αxj
                    end
                end
            end
        end
        C
    end
end

end # module


text:
#=move_user_move.jl > moves the piece at the source into the target area. If promote
is equal to ”T” then promote the piece after it has completed the move.
Accepts 6 command line args : <filename> => database, <xsource> => xSource
<ysource> => ySource <xtarget> =>xTarget <ytarget> =>yTarget <promote> => promote
=#
#include("square.jl")
include("dParse.jl")
module move_user_move

  using SQLite
  function moveUserMove(database,sourcex,sourcey,targetx,targety,promote::Bool)

    board = Parse(database)
  #=  database = ARGS[1] #/path/to/database/file {string}
    sourcey = parse(Int,chomp(ARGS[2])) #x coordinate of piece {Int}
    sourcex = parse(Int,chomp(ARGS[3])) #y coordinate of piece {Int}
    targetx = parse(Int,chomp(ARGS[4])) #x coordinate to place piece {Int}
    targety = parse(Int,chomp(ARGS[5])) #y coordinate to place piece {Int}
    promoteDict = Dict("T" => true, "F" => false)
    shouldPromote = promoteDict[ARGS[6]] #Promote the piece after moving {bool}=#
    db = SQLite.DB(database) #Opens the database gamefile

      #= ---- Determines the move_number ---- =#

      res = SQLite.query(db,"""SELECT MAX("move_number") FROM moves;""") #Finds the last played move (maximum move_number)
      lastMoveIDNullable = res[1,1] #SQL query with max move_number (POSSIBLY "NULL" if no moves have been made)

      if (!isnull(lastMoveIDNullable)) #Checks that lastMoveID was not NULL
        lastMoveID = get(res[1,1]) #Parses the last move move_number as an Int

      else #lastMoveID is NULL
        lastMoveID = 0 #move_number "0" is unsused. This implies no moves have been made

      end
      move_number = lastMoveID+1 #Current move number

    #=-----UPDATE DATABASE & BOARD W/MOVE-----=#
    if shouldPromote #Option will be '!'
      SQLite.query(db,"""INSERT INTO "moves" (move_number, move_type, sourcex, sourcey, targetx, targety, option)
      VALUES ($(move_number),'move',$sourcex, $sourcey,$targetx, $targety, !);""")
      board[(targetx,targety)] = board[(sourcex,sourcey)] #Updates the board before next move
      board[(targetx,targety)][1]="promoted$(board[(targetx,targety)][1])"
      delete!(board,(sourcex,sourcey))
    else #Option will be NULL
      SQLite.query(db,"""INSERT INTO "moves" (move_number, move_type, sourcex, sourcey, targetx, targety)
      VALUES ($(move_number),'move',$sourcex, $sourcey,$targetx, $targety);""")
      board[(targetx,targety)] = board[(sourcex,sourcey)] #Updates the board before next move
      delete!(board,(sourcex,sourcey))
    end
  #ST.saveBoard(board)
end
end

max_stars_repo_name: kinoroy/Shogi

max_issues_repo_name: kinoroy/Shulia

max_forks_repo_name: kinoroy/Shulia


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:4] <= 1)

                     @variable(m, -11.08168115564476 <= q <= 12.629895565298366)

                     add_NL_constraint(m, :(swish(-0.8013086439641905 + -0.577380120181807*swish(-0.15161094843442724 + 0.3901524114721284*$(x[1]) + -0.7620628362671806*$(x[2]) + 0.5531251925343885*$(x[3]) + 0.8877751701397409*$(x[4])) + -0.9073468716494122*swish(-0.9517329762307463 + 0.7324043290359472*$(x[1]) + -0.6011940599882695*$(x[2]) + 0.8132214951552479*$(x[3]) + 0.2016959253369457*$(x[4])) + -0.34890000179172853*swish(0.5700564024677015 + -0.6508434167160906*$(x[1]) + 0.3531203702734036*$(x[2]) + 0.9592602519851474*$(x[3]) + -0.2457804299770685*$(x[4]))) + swish(-0.6400914085711236 + 0.5999574270201746*swish(-0.15161094843442724 + 0.3901524114721284*$(x[1]) + -0.7620628362671806*$(x[2]) + 0.5531251925343885*$(x[3]) + 0.8877751701397409*$(x[4])) + -0.3302091884595031*swish(-0.9517329762307463 + 0.7324043290359472*$(x[1]) + -0.6011940599882695*$(x[2]) + 0.8132214951552479*$(x[3]) + 0.2016959253369457*$(x[4])) + 0.9262651374928534*swish(0.5700564024677015 + -0.6508434167160906*$(x[1]) + 0.3531203702734036*$(x[2]) + 0.9592602519851474*$(x[3]) + -0.2457804299770685*$(x[4]))) + swish(0.7171327203675486 + 0.9873737246227612*swish(-0.15161094843442724 + 0.3901524114721284*$(x[1]) + -0.7620628362671806*$(x[2]) + 0.5531251925343885*$(x[3]) + 0.8877751701397409*$(x[4])) + -0.03055409136192777*swish(-0.9517329762307463 + 0.7324043290359472*$(x[1]) + -0.6011940599882695*$(x[2]) + 0.8132214951552479*$(x[3]) + 0.2016959253369457*$(x[4])) + 0.20254351925716962*swish(0.5700564024677015 + -0.6508434167160906*$(x[1]) + 0.3531203702734036*$(x[2]) + 0.9592602519851474*$(x[3]) + -0.2457804299770685*$(x[4]))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
using ApproxFun, LinearAlgebra, Test
    import ApproxFun: testbandedblockbandedoperator, testblockbandedoperator, testraggedbelowoperator, Block

@testset "PDE" begin
    @testset "Zero Dirichlet" begin
        S = JacobiWeight(1.,1.,Jacobi(1.,1.))^2
        Δ = Laplacian(S)

        testbandedblockbandedoperator(Δ)

        u = Fun((x,y)->sin(π*x)*sin(π*y),S)
        f = -2π^2*u

        F = qr(Δ)
        ApproxFun.resizedata!(F,:,1000)
        @time v=F\f
        @test norm((u-v).coefficients)<100eps()

        F=qr(Δ)
        ApproxFun.resizedata!(F.R_cache,:,100)
        ApproxFun.resizedata!(F.R_cache,:,1000)
        @time v=F \ f
        @test norm((u-v).coefficients)<100eps()

        F=qr(Δ)
        @time v=F\f
        @test norm((u-v).coefficients)<100eps()
    end

    @testset "Rectangle Laplace/Poisson" begin
        dx = dy = ChebyshevInterval()
        d = dx × dy
        g = Fun((x,y)->exp(x)*cos(y),∂(d))

        B = Dirichlet(d)

        testblockbandedoperator(B)
        testbandedblockbandedoperator(Laplacian(d))
        testbandedblockbandedoperator(Laplacian(d)+0.0I)

        A=[Dirichlet(d);Laplacian(d)]

        @time u=A\[g,0.]
        @test u(.1,.2) ≈ real(exp(0.1+0.2im))

        A=[Dirichlet(d);Laplacian(d)+0.0I]
        @time u=A\[g,0.]

        @test u(.1,.2) ≈ real(exp(0.1+0.2im))

        ## Poisson
        f=Fun((x,y)->exp(-10(x+.2)^2-20(y-.1)^2),ChebyshevInterval()^2,500)  #default is [-1,1]^2
        d=domain(f)
        A=[Dirichlet(d);Laplacian(d)]
        @time  u=\(A,[zeros(∂(d));f];tolerance=1E-7)
        @test ≈(u(.1,.2),-0.04251891975068446;atol=1E-5)
    end

    @testset "Bilaplacian" begin
        dx = dy = ChebyshevInterval()
        d = dx × dy
        Dx = Derivative(dx); Dy = Derivative(dy)
        L = Dx^4⊗I + 2*Dx^2⊗Dy^2 + I⊗Dy^4

        testbandedblockbandedoperator(L)

        B = Dirichlet(dx) ⊗ Operator(I,dy)
        testraggedbelowoperator(B)

        A=[Dirichlet(dx) ⊗ Operator(I,dy);
                Operator(I,dx)  ⊗ Dirichlet(dy);
                Neumann(dx) ⊗ Operator(I,dy);
                Operator(I,dx) ⊗ Neumann(dy);
                 L]

        testraggedbelowoperator(A)

        @time u=\(A,[[1,1],[1,1],[0,0],[0,0],0];tolerance=1E-5)
        @test u(0.1,0.2) ≈ 1.0
    end

    @testset "Periodic x Interval" begin
        d=PeriodicSegment() × ChebyshevInterval()

        u_ex=Fun((x,y)->real(cos(x+im*y)),d)
        @test u_ex(1.0,0.1) ≈ real(cos(1.0+im*0.1)) atol=10eps()

        B=Dirichlet(Space(d))

        @test B.order == 0  # tests stupid bug
        g=Fun((x,y)->real(cos(x+im*y)),rangespace(B))  # boundary data

        @test norm((B*u_ex-g).coefficients) < 100eps()

        testbandedblockbandedoperator(Laplacian(d))

        @time u=[B;Laplacian(d)]\[g;0.]

        @test u(.1,.2) ≈ real(cos(.1+.2im))
    end

    @testset "Schrodinger" begin
        dx=0..1; dt=0.0..0.001
        C=Conversion(Chebyshev(dx)*Ultraspherical(1,dt),Ultraspherical(2,dx)*Ultraspherical(1,dt))
        testbandedblockbandedoperator(C)
        testbandedblockbandedoperator(Operator{ComplexF64}(C))

        d = dx × dt

        x,y = Fun(d)
        @test x(0.1,0.0001) ≈ 0.1
        @test y(0.1,0.0001) ≈ 0.0001

        V = x^2

        Dt=Derivative(d,[0,1]);Dx=Derivative(d,[1,0])

        ϵ = 1.
        u0 = Fun(x->exp(-100*(x-.5)^2)*exp(-1/(5*ϵ)*log(2cosh(5*(x-.5)))),dx)
        L = ϵ*Dt+(.5im*ϵ^2*Dx^2)
        testbandedblockbandedoperator(L)

        @time u = \([timedirichlet(d);L],[u0,[0.,0.],0.];tolerance=1E-5)
        @test u(0.5,0.001) ≈ 0.857215539785593+0.08694948835021317im  # empircal from ≈ schurfact
    end

    @testset "check we dispatch correctly to get fast build" begin
        S = JacobiWeight(1.,1.,Jacobi(1.,1.))^2
        Δ = Laplacian(S)
        @time S = view(Δ.op.ops[1].ops[1].op,Block.(1:40), Block.(1:40))
        @test typeof(S.parent.domaintensorizer) == ApproxFun.Trivial2DTensorizer
    end
end


text:
module io

function guess_type(x)
    try
        return eval(parse(x))
    catch
        return x
    end
end

function split(word::AbstractString, delimiter::AbstractString)
    words = []
    lastsearchpos = -1:0
    while true
        searchpos = search(word, delimiter, lastsearchpos[end]+1)
        if searchpos==0:-1
            push!(words, word[lastsearchpos[end]+1:end])
            break
        end
        push!(words, word[lastsearchpos[end]+1:searchpos[1]-1])
        lastsearchpos = searchpos
    end
    return words
end

function join(delimiter::AbstractString, words::Vector)
    result = words[1]
    for i=2:length(words)
        result = result*delimiter*words[i]
    end
    return result
end

function dict2filename(d::Dict; extension=".dat")
    pairs = ["$(key)=$(item)" for (key, item) in d]
    return join(";", sort(pairs))*extension
end

function filename2dict(name::AbstractString; extension=".dat")
    @assert endswith(name, extension)
    items = split(name[1:end-length(extension)], ";")
    D = Dict()
    for item in items
        array = split(item, "=")
        key = array[1]
        value = array[2]
        D[key] = guess_type(value)
    end
    return D
end

end # module


text:
@testset "basic" begin
   @test get_E10(1.0) ≈ 0.3760301626167376u"eV"
   @test get_E10(1u"nm") ≈ 0.3760301626167376u"eV"
   @test get_E10(1u"nm", me=1.0) ≈ 0.3760301626167376u"eV"
end


text:
using ReverseGeocode
using Test
using StaticArrays
using NearestNeighbors

test_locs = [
(tag = (country="Czechia", country_code="CZ", city="Olomouc"), latlon = SA[49.5863897, 17.2627342]),
(tag = (country="Germany", country_code="DE", city="Bad Elster"), latlon = SA[50.3005700, 12.2091950]),
(tag = (country="Czechia", country_code="CZ", city="Hranice"), latlon = SA[50.3050656, 12.1865356]),
(tag = (country="Norway", country_code="NO", city="Meråker"), latlon = SA[63.3342550, 12.0280064]),
(tag = (country="Norway", country_code="NO", city="Meråker"), latlon = SA[63.2887794, 12.1626800]), # cordinates are in fact 5 km from the Norwegian border in Sweden, by the town of Storlien
]

@testset "ReverseGeocode.jl" begin
    gc = Geocoder(;data_dir="./data", geo_file="test_cities")

    # Test individually
    for loc in test_locs
        @test loc.tag == ReverseGeocode.decode(gc, loc.latlon)
    end

    # Test with array of SAs
    @test [loc.tag for loc in test_locs] == ReverseGeocode.decode(gc, [loc.latlon for loc in test_locs])

    # Test with Matrix
    test_loc_matrix = Matrix(hcat([l.latlon for l in test_locs]...))
    @test [loc.tag for loc in test_locs] == ReverseGeocode.decode(gc, test_loc_matrix)
end


text:
using Documenter, MessyTimeSeries;
makedocs(
    sitename="MessyTimeSeries.jl",
    pages = [
        "index.md",
        "man/getting_started.md",
        "Manual" => ["man/kalman.md", "man/subsampling.md"],
        "Appendix" => ["man/methods.md"],
    ]
);

deploydocs(
    repo = "github.com/fipelle/MessyTimeSeries.jl.git",
)


text:
using CUDA
using ExaPF
using KernelAbstractions
using Test
using Printf
using TimerOutputs

import ExaPF: PowerSystem, LinearSolvers, AutoDiff

# Newton-Raphson tolerance
ntol = 1e-6
# For debugging in REPL use the following lines
# empty!(ARGS)
# push!(ARGS, "KrylovBICGSTAB")
# push!(ARGS, "CPU")
# push!(ARGS, "case300.m")
# push!(ARGS, "caseGO30R-025.raw")

# We do need the time in ms, and not with time units all over the place
function TimerOutputs.prettytime(t)
    value = t / 1e6 # "ms"

    if round(value) >= 100
        str = string(@sprintf("%.0f", value))
    elseif round(value * 10) >= 100
        str = string(@sprintf("%.1f", value))
    elseif round(value * 100) >= 100
        str = string(@sprintf("%.2f", value))
    else
        str = string(@sprintf("%.3f", value))
    end
    return lpad(str, 6, " ")
end

function printtimer(timers, key::String)
   prettytime(timers[key].accumulated_data.time)
end

linsolver = eval(Meta.parse("LinearSolvers.$(ARGS[1])"))
device = eval(Meta.parse("$(ARGS[2])()"))
datafile = joinpath(dirname(@__FILE__), ARGS[3])
pf = PowerSystem.PowerNetwork(datafile)
polar = PolarForm(pf, device)
cache = ExaPF.get(polar, ExaPF.PhysicalState())
jx = AutoDiff.Jacobian(polar, ExaPF.power_balance, State())
npartitions = ceil(Int64,(size(jx.J,1)/64))
if npartitions < 2
    npartitions = 2
end
precond = ExaPF.LinearSolvers.BlockJacobiPreconditioner(jx.J, npartitions, device)
# Retrieve initial state of network
x0 = ExaPF.initial(polar, State())
u0 = ExaPF.initial(polar, Control())

algo = linsolver(precond)
powerflow_solver = NewtonRaphson(tol=ntol)
nlp = ExaPF.ReducedSpaceEvaluator(polar;
                                    linear_solver=algo, powerflow_solver=powerflow_solver)
convergence = ExaPF.update!(nlp, u0)
ExaPF.reset!(nlp)
convergence = ExaPF.update!(nlp, u0)
ExaPF.reset!(nlp)
TimerOutputs.reset_timer!(ExaPF.TIMER)
convergence = ExaPF.update!(nlp, u0)

# Make sure we are converged
@assert(convergence.has_converged)

# Output
prettytime = TimerOutputs.prettytime
timers = ExaPF.TIMER.inner_timers
inner_timer = timers["Newton"]
if ARGS[1] == "DirectSolver"
    println("$(ARGS[1]), $(ARGS[2]), $(ARGS[3]),",
            printtimer(timers, "Newton"),
            ", $(convergence.has_converged)")
else
    println("$(ARGS[1]), $(ARGS[2]), $(ARGS[3]),",
            printtimer(timers, "Newton"),",",
            printtimer(inner_timer, "Jacobian"),",",
            printtimer(inner_timer, "Linear Solver"),
            ", $(convergence.has_converged)")
end


text:
"""
APIInvoker holds the transport and format used for a remote api call.
"""
struct APIInvoker{T<:AbstractTransport,F<:AbstractMsgFormat}
    transport::T
    format::F
end

"""
APIInvoker holds the transport and format used for a remote api call.
This method creates an APIInvoker over ZMQ transport using JSON message format.
(provided for backward compatibility)
"""
APIInvoker(addr::String, ctx::Context=Context()) = APIInvoker(ZMQTransport(addr, REQ, false, ctx), JSONMsgFormat())
APIInvoker(ip::IPv4, port::Int, ctx::Context=Context()) = APIInvoker(ZMQTransport(ip, port, REQ, false, ctx), JSONMsgFormat())

"""
Calls a remote api `cmd` with `args...` and `data...`.
The response is formatted as specified by the formatter specified in `conn`.
"""
function apicall(conn::APIInvoker{T,F}, cmd, args...; data...) where {T<:AbstractTransport, F<:AbstractMsgFormat}
    req = wireformat(conn.format, String(cmd), args...; data...)
    resp = sendrecv(conn.transport, req)
    juliaformat(conn.format, resp)
end


text:
mutable struct LIP
    dt::Float64
    t::Float64
    T_sup::Float64 #support time
    zc::Float64

    #desired foot location
    p_x::Float64
    p_y::Float64
    p_x_star::Float64 #modified foot location x
    p_y_star::Float64

    #Gait parameters
    s_x::Float64
    s_y::Float64

    #COM init conditions
    x_0::Float64
    vx_0::Float64
    y_0::Float64
    vy_0::Float64

    #COM real-time state
    x_t::Float64
    vx_t::Float64
    y_t::Float64
    vy_t::Float64

    #COM desired state
    x_d::Float64
    vx_d::Float64
    y_d::Float64
    vy_d::Float64

    #final state for one gait unit
    bar_x::Float64
    bar_vx::Float64
    bar_y::Float64
    bar_vy::Float64

    support_leg::String
    left_foot_pos::Array{Float64}
    right_foot_pos::Array{Float64}
    COM_pos::Array{Float64}
    T_c::Float64
    C::Float64
    S::Float64

    #physical params
    m::Float64
    h::Float64

    function LIP(;dtt=0.001, Ts=1.0, sl="left_leg", mass=1.0, height=1.0)
        t = 0.0
        T_sup = Ts

        p_x=0
        p_y = 0.0
        p_x_star = 0.0 #modified foot location x
        p_y_star = 0.0
        support_leg = sl

        #Gait parameters
        s_y = 0.0
        s_x = 0.0

        #COM init conditions
        x_0 = 0.0
        vx_0 = 0.0
        y_0 = 0.0
        vy_0 = 0.0

        #COM real-time state
        x_t = 0.0
        vx_t = 0.0
        y_t = 0.0
        vy_t = 0.0
        zc = 0.0

        #COM desired state
        x_d = 0.0
        vx_d = 0.0
        y_d = 0.0
        vy_d = 0.0

        #final state for one gait unit
        bar_x = 0.0
        bar_vx = 0.0
        bar_y = 0.0
        bar_vy = 0.0

        left_foot_pos = SA[0.0,0.0,0.0]
        right_foot_pos = SA[0.0,0.0,0.0]
        COM_pos = SA[0.0,0.0,0.0]
        T_c = 0.0
        C = 0.0
        S = 0.0
        dt = dtt

        m = mass
        h = height

        new(dt, t, T_sup, zc, p_x, p_y, p_x_star, p_y_star, s_x, s_y, x_0, vx_0,
        y_0, vy_0, x_t, vx_t, y_t, vy_t, x_d, vx_d, y_d, vy_d, bar_x, bar_vx,
        bar_y, bar_vy, support_leg, left_foot_pos, right_foot_pos, COM_pos, T_c,
        C, S, m, h)

    end
end

function initializeModel!(model, COM_pos, left_foot_pos, right_foot_pos)
        model.COM_pos = COM_pos
        model.left_foot_pos = left_foot_pos
        model.right_foot_pos = right_foot_pos

        model.zc = model.h
        model.T_c = sqrt(model.zc / 9.81) # set gravity parameter as 9.81
        model.C = cosh(model.T_sup/model.T_c)
        model.S = sinh(model.T_sup/model.T_c)
end

function updateParameters!(model, T_sup)
        model.T_sup = T_sup
        model.C = cosh(model.T_sup/model.T_c)
        model.S = sinh(model.T_sup/model.T_c)
end

function step!(model)
        model.t += model.dt
        t = model.t
        T_c = model.T_c

        model.x_t = model.x_0*cosh(t/T_c) + T_c*model.vx_0*sinh(t/T_c)
        model.vx_t = model.x_0/T_c*sinh(t/T_c) + model.vx_0*cosh(t/T_c)

        model.y_t = model.y_0*cosh(t/T_c) + T_c*model.vy_0*sinh(t/T_c)
        model.vy_t = model.y_0/T_c*sinh(t/T_c) + model.vy_0*cosh(t/T_c)
end

function calculateXtVt!(model, t)
        T_c = model.T_c
        x_t = model.x_0*cosh(t/T_c) + T_c*model.vx_0*sinh(t/T_c)
        vx_t = model.x_0/T_c*sinh(t/T_c) + model.vx_0*cosh(t/T_c)

        y_t = model.y_0*cosh(t/T_c) + T_c*model.vy_0*sinh(t/T_c)
        vy_t = model.y_0/T_c*sinh(t/T_c) + model.vy_0*cosh(t/T_c)

        return x_t, vx_t, y_t, vy_t
end

function nextReferenceFootLocation!(model, s_x, s_y, theta=0)
        if model.support_leg == "left_leg" # then the next support leg is the right leg
            p_x_new = model.p_x + cos(theta)*s_x - sin(theta)*s_y
            p_y_new = model.p_y + sin(theta)*s_x + cos(theta)*s_y
        elseif model.support_leg == "right_leg" # then the next support leg is the left leg
            p_x_new = model.p_x + cos(theta)*s_x + sin(theta)*s_y
            p_y_new = model.p_y + sin(theta)*s_x - cos(theta)*s_y
        end

        return p_x_new, p_y_new
end

function nextState!(model, s_x, s_y, theta=0.0)
        if model.support_leg == "left_leg"
            bar_x_new = cos(theta)*s_x/2.0 - sin(theta)*s_y/2.0
            bar_y_new = sin(theta)*s_x/2.0 + cos(theta)*s_y/2.0
        elseif model.support_leg == "right_leg"
            bar_x_new = cos(theta)*s_x/2.0 + sin(theta)*s_y/2.0
            bar_y_new = sin(theta)*s_x/2.0 - cos(theta)*s_y/2.0
        end
        return bar_x_new, bar_y_new
end

function nextVel!(model, bar_x=0, bar_y=0, theta=0)
        C = model.C
        S = model.S
        T_c = model.T_c

        bar_vx_new = cos(theta)*(1+C)/(T_c*S)*bar_x -
                                sin(theta)*(C-1)/(T_c*S)*bar_y
        bar_vy_new = sin(theta)*(1+C)/(T_c*S)*bar_x +
                                cos(theta)*(C-1)/(T_c*S)*bar_y

        return bar_vx_new, bar_vy_new
end

function targetState!(model, p_x, bar_x, bar_vx)
    x_d = p_x + bar_x
    vx_d = bar_vx

    return x_d, vx_d
end

function modifiedFootLocation!(model, a=1.0, b=1.0, x_d=0, vx_d=0,
                                x_0=0, vx_0=0)
    C = model.C
    S = model.S
    T_c = model.T_c
    D = a*(C - 1)^2 + b*(S/T_c)^2
    p_x_star = -a*(C-1)*(x_d - C*x_0 - T_c*S*vx_0)/D - b*S*(vx_d -
                                        S*x_0/T_c - C*vx_0)/(T_c*D)
    return p_x_star
end

function calculateFootLocationForNextStep!(model, s_x=0.0, s_y=0.0, a=1.0,
                     b=1.0, theta=0.0, x_0=0.0, vx_0=0.0, y_0=0.0, vy_0=0.0)
    model.s_x = s_x
    model.s_y = s_y
    p_x_new, p_y_new = nextReferenceFootLocation!(model,s_x, s_y, theta)
    bar_x, bar_y = nextState!(model,s_x, s_y, theta)
    bar_vx, bar_vy = nextVel!(model,bar_x, bar_y, theta)

    model.x_d, model.vx_d = targetState!(model,p_x_new, bar_x, bar_vx)
    model.y_d, model.vy_d = targetState!(model,p_y_new, bar_y, bar_vy)
    model.p_x_star = modifiedFootLocation!(model,a, b, model.x_d,
                                            model.vx_d, x_0, vx_0)
    model.p_y_star = modifiedFootLocation!(model,a, b, model.y_d,
                                            model.vy_d, y_0, vy_0)
end

function switchSupportLeg!(model)
    if model.support_leg == "left_leg"
        model.support_leg = "right_leg"
        COM_pos_x = model.x_t + model.left_foot_pos[1]
        COM_pos_y = model.y_t + model.left_foot_pos[2]
        model.x_0 = COM_pos_x - model.right_foot_pos[1]
        model.y_0 = COM_pos_y - model.right_foot_pos[2]
    elseif model.support_leg == "right_leg"
        model.support_leg = "left_leg"
        COM_pos_x = model.x_t + model.right_foot_pos[1]
        COM_pos_y = model.y_t + model.right_foot_pos[2]
        model.x_0 = COM_pos_x - model.left_foot_pos[1]
        model.y_0 = COM_pos_y - model.left_foot_pos[2]
    end

    model.t = 0
    model.vx_0 = model.vx_t
    model.vy_0 = model.vy_t
end

function at_goal(current_position, goal_position)
    if current_position[1] >= goal_position[1]
        return true
    else
        return false
    end
end


text:
using Flux
using BSON: @load

# Play game in different environment only by
# importing Flux models
representation = Flux.@load "ttt_representation.bson"
dynamics = Flux.@load "ttt_dynamics.bson"
policyvalue = Flux.@load "ttt_policyvalue.bson"

tic_tac_toe_image = zeros(Float32, 3, 3, 3)

hidden_state = representation(tic_tac_toe_image)
policy_logits, value = policyvalue(hidden_state)

first_action = 5
hidden_state_action = cat(hidden_state, ones(3,3); dims=3)
next_hidden_state, reward = dynamics(hidden_state_action)

second_action = 1
hidden_state_action = cat(next_hidden_state, ones(3,3); dims=3)
next_hidden_state, reward = dynamics(hidden_state_action)


text:
function make_vertical_periodic!(mesh)
    bottomcellids = bottom_cellids(mesh)
    topcellids = top_cellids(mesh)

    @assert length(bottomcellids) == length(topcellids)
    for idx in 1:length(bottomcellids)
        update_cell_connectivity!(mesh,1,bottomcellids[idx],topcellids[idx])
        update_cell_connectivity!(mesh,3,topcellids[idx],bottomcellids[idx])
    end
end

function make_horizontal_periodic!(mesh)
    leftcellids = left_cellids(mesh)
    rightcellids = right_cellids(mesh)

    @assert length(leftcellids) == length(rightcellids)
    for idx in 1:length(leftcellids)
        update_cell_connectivity!(mesh,4,leftcellids[idx],rightcellids[idx])
        update_cell_connectivity!(mesh,2,rightcellids[idx],leftcellids[idx])
    end
end


text: include("linearoperator.jl")