text:
# 	This file is part of Maxima.jl. It is licensed under the MIT license
#   Copyright (c) 2016 Nathan Smith
#                 2021 Jingyue Wang

export string,
       show

import Base: string,
             show
#import Base.Multimedia: display

string(m::MExpr) = convert(String, m)

show(io::IO, m::MExpr) = print(io, convert(String, m))
function show(io::IO, ::MIME"text/plain", m::MExpr)
    check = "'(" * replace(convert(String, m), r";" => ");\n'(") * ")"
    write(ms.input, "$(replace(check, r";" => "\$"))\$\n print(ascii(4))\$")
    out = (readuntil(ms.output, EOT) |> String
                                     |> str -> rstrip(str, EOT))
    if occursin(synerr, out) || occursin(maxerr, out)
        @warn "Invalid Maxima expression"
        print(io, out)
    else
        mcall(mx"display2d: true")
        write(ms, replace(check, r";" => "; print(ascii(3))\$ "))
        str = rstrip(read(ms), '\n')
        mcall(mx"display2d: false")
        sp = split(str, ETX)
        for k in sp; print(io, k); end
    end
end

function show(io::IO, ::MIME"text/latex", m::MExpr)
    check = "'(" * replace(convert(String, m), r";" => ")\$\n'(") * ")"
    write(ms.input, "$check\$\n print(ascii(4))\$")
    out = (readuntil(ms.output, EOT) |> String
                                     |> str -> rstrip(str, EOT))
    if occursin(synerr, out) || occursin(maxerr, out)
        @warn "Invalid Maxima expression"
        print(io, out)
    else
        write(ms, "tex1('("*replace(convert(String, m), r";" => "))\$\ntex1('(") * "))")
        texstr = read(ms)
        print(io, replace(texstr, r"\nfalse\n" => ""))
    end
end

# A second approach for this function, which probably is better, is first
# converting `m` into a Julia expression, then use `MathMLRepr.jl`
function show(io::IO, ::MIME"text/mathml+xml", m::MExpr)
    check = "'(" * replace(convert(String, m), r";" => ")\$\n'(") * ")"
    write(ms.input, "$check\$\n print(ascii(4))\$")
    out = (readuntil(ms.output, EOT) |> String
                                     |> str -> rstrip(str, EOT))
    if occursin(synerr, out) || occursin(maxerr, out)
        @warn "Invalid Maxima expression"
        print(io, "Invalid Maxima Expression\n")
        print(io, out)
    else
        write(ms, "tex1('("*replace(convert(String, m), r";" => "))\$\ntex1('(") * "))")
        texstr = read(ms)
        s = replace(texstr, r"\nfalse\n" => "")
        # Warning! quick hacks inside this function
        s = clean_matrix_env(s)
        ## debug
        #write(io, s)
        ## Translate s into MathML.
        ts = TokenStream(s)
        nodes = exlist_to_mml(parse!(ts))
        if ts.top !== nothing && ts.top.val == "inline"
            o = mml_element(:math, nodes)
        else
            o = mml_element_attr(:math, nodes, "display=\"block\"")
        end
        write(io, o)
    end
end


text:
# GitHub
"""
    github_endpoint = "https://api.github.com/graphql"
Endpoint for the GitHub API v4 (GraphQL).
"""
const github_endpoint = "https://api.github.com/graphql";
"""
    client = GraphQLClient(github_endpoint,
                           auth = "bearer \$github_token",
                           headers = github_header)
Client for the GitHub API.
"""
const client = GraphQLClient(github_endpoint,
                             auth = "bearer $github_token",
                             headers = github_header)
"""
    github_api_query::String
Queries for finding open-sourced projects and their commit information from GitHub.
"""
const github_api_query = """
    query Search(\$license_created: String!) {
      ...RateLimit
      search(query: \$license_created,
             type: REPOSITORY,
             first: 100) {
        ...SearchLogic
      }
    }
    query SearchCursor(\$license_created: String!,
                 \$cursor: String!) {
      ...RateLimit
      search(query: \$license_created,
             type: REPOSITORY,
             first: 100,
             after: \$cursor) {
        ...SearchLogic
      }
    }
    query Repository(\$owner: String!,
                     \$name: String!,
                     \$until: GitTimestamp!) {
      ...RateLimit
      ...Repo
    }
    query RepositoryCursor(\$owner: String!,
                           \$name: String!,
                           \$until: GitTimestamp!,
                           \$cursor: String!) {
      ...RateLimit
      ...RepoCursor
    }
    fragment RateLimit on Query {
      rateLimit {
        remaining
        resetAt
      }
    }
    fragment SearchLogic on SearchResultItemConnection {
      repositoryCount
      pageInfo {
        endCursor
        hasNextPage
      }
      nodes {
        ... on Repository {
          isPrivate
          databaseId
          nameWithOwner
          createdAt
          isArchived
          isFork
          isMirror
        }
      }
    }
    fragment CommitHistory on CommitHistoryConnection {
      totalCount
      pageInfo {
        endCursor
        hasNextPage
      }
      nodes {
        author {
          user {
            login
          }
        }
        oid
        committedDate
        additions
        deletions
      }
    }
    fragment Repo on Query {
      repository(owner: \$owner, name: \$name) {
        defaultBranchRef {
          target {
            ... on Commit {
              history(first: 1, until: \$until) {
                ...CommitHistory
              }
            }
          }
        }
      }
    }
    fragment RepoCursor on Query {
      repository(owner: \$owner, name: \$name) {
        defaultBranchRef {
          target {
            ... on Commit {
              history(first: 100, until: \$until, after: \$cursor) {
                ...CommitHistory
              }
            }
          }
        }
      }
    }
    """;
"""
    binary_search_dt_interval(license::AbstractString,
                              interval::AbstractString)::data, as_of, created_at
Given a license and a datetime interval, it will use binary search to find
a datetime interval with no more than 1,000 results.
"""
@inline function binary_search_dt_interval(license::AbstractString,
                                           created_at::AbstractString)
    dt_start = match(r".*(?=\.{2})", created_at)
    if isnothing(dt_start)
        dt_start = replace(created_at, r"Z$" => "+00:00") |>
            (dt -> ZonedDateTime(dt, github_dtf))
    else
        dt_start = match(r".*(?=\.\.)", created_at).match |>
            (dt -> replace(dt, r"Z$" => "+00:00")) |>
            (dt -> ZonedDateTime(dt, github_dtf))
    end
    dt_end = match(r"(?<=\.{2}).*", created_at)
    if isnothing(dt_end)
        dt_end = until
    else
        dt_end = match(r"(?<=\.{2}).*", created_at).match |>
            (dt -> replace(dt, r"Z$" => "+00:00")) |>
            (dt -> ZonedDateTime(dt, github_dtf))
    end
    result = client.Query(github_api_query,
                          operationName = "Search",
                          vars = Dict("license_created" =>
                                      """
                                      license:$license
                                      created:$dt_start..$dt_end
                                      """))
    as_of = get_as_of(result.Info)
    json = JSON3.read(result.Data)
    @assert(haskey(json, :data))
    data = json.data
    github_wait_out(data.rateLimit)
    repositoryCount = data.search.repositoryCount
    while repositoryCount > 1_000
        dt_end = dt_start + (dt_end - dt_start) ÷ 2 |>
            (dt -> format(dt, github_dtf)) |>
            (dt -> ZonedDateTime(dt, github_dtf))
        result = client.Query(github_api_query,
                              operationName = "Search",
                              vars = Dict("license_created" =>
                                          """
                                          license:$license
                                          created:$dt_start..$dt_end
                                          """))
        as_of = get_as_of(result.Info)
        json = JSON3.read(result.Data)
        @assert(haskey(json, :data))
        data = json.data
        github_wait_out(data.rateLimit)
        repositoryCount = data.search.repositoryCount
    end
    data.search, as_of, "$dt_start..$dt_end"
end
"""
    get_as_of(response::Response)::String
Returns the zoned date time when the response was returned.
"""
get_as_of(response::Response) =
    response.headers[findfirst(x -> isequal("Date", x.first),
                               response.headers)].second[6:end] |>
        (dt -> ZonedDateTime(dt, response_dtf)) |>
        string
"""
    github_wait_out(rateLimit)
If the GitHub personal access token has exhausted the per hour limit,
it waits until it resets.
"""
github_wait_out(rateLimit) =
    iszero(rateLimit.remaining) &&
        ZonedDateTime(replace("2019-08-14T23:02:00-04:00", r"Z$" => "+00:00"),
                      "yyyy-mm-ddTHH:MM:SSzzzzz") - now(localzone()) |>
            sleep

ext: jl


text:
include("functions_log_test.jl")
include("handler_functions_events.jl")
include("handler_functions_websockets_general_test.jl")
include("handler_functions_websockets_subprotocol_test.jl")

"""
This function returns responses to http requests.
For serving the HTML javascript pages which then open WebSocket clients.
"""
function httphandle(request::Request, response::Response)
    global n_responders
    id = "server_functions.httphandle\t"
    clog(id, :cyan, request, "\n")
    if request.method=="GET"
        if request.resource =="/favicon.ico"
            response =  HttpServer.FileResponse(joinpath(@__DIR__, "favicon.ico"))
        elseif startswith(request.resource, "/browser")
            response =  HttpServer.FileResponse(joinpath(@__DIR__, splitdir(request.resource)[2]))
            if request.resource == "/browsertest.html"
                n_responders += 1
            end
        else
            response = Response(404, "$id, can't serve $request.resource.")
        end
    else
        response = Response(404, "$id, unexpected request.")
    end
    clog(id, response, "\n")
    push!(response.headers, "Connection" => "close")
    return response
end

"""
When a connection performs the handshake successfully, this function
is called as a separate task. Based on inspecting the request that was
accepted as an upgrade, it in turn calls one of three functions.
The function call does not return until it is time to close the websocket.
"""
function gatekeeper(wsrequest::Request, websocket::WebSocket)
    id = "server_functions.gatekeeper\t"
    clog(id, :cyan, wsrequest, "\t", :yellow, websocket, "\n")
    if subprotocol(wsrequest) == "websocket-testprotocol"
        ws_test_protocol(websocket)
    elseif subprotocol(wsrequest) == "websocket-test-binary"
        ws_test_binary(websocket)
    elseif subprotocol(wsrequest) != ""
        clog(id, :red, "Unknown sub protocol let through, not responding further. \n")
    else
        ws_general(websocket)
    end
    clog(id, "Exiting \n")
    nothing
end

function start_ws_server_async()
    id = "server_functions.start_ws_server\t"
    # Specify this subprotocol is to be let through to gatekeeper:
    WebSockets.addsubproto("websocket-testprotocol")
    WebSockets.addsubproto("websocket-test-binary")
    # Tell HttpHandler which functions to spawn when something happens.
    httpha = HttpHandler(httphandle)
    httpha.events["error"]  = ev_error
    httpha.events["listen"] = ev_listen
    httpha.events["connect"] = ev_connect
    httpha.events["close"] = ev_close
    httpha.events["write"] = ev_write
    httpha.events["reset"] = ev_reset
    # Pack the gatekeeper function in a recognizable function wrapper
    wsh = WebSocketHandler(gatekeeper)
    server = Server(httpha, wsh )
    clog(id, "Server to be started:\n", server )
    servertask = @async run( server, 8080)
    clog(id, "Server listening on 127.0.0.1:8080\n")
    server
end
"Closes websockets and server; other existing sockets, if any, remain open."
function closeall()
    if isdefined(:WEBSOCKETS)
        for va in values(WEBSOCKETS)
            if isopen(va)
                close(va)
            end
        end
        if isdefined(:server)
            close(server)
        end
    end
end
"Counts the number of open websocket in WEBSOCKETS dictionary"
function count_open_websockets()
    count = 0
    for w in values(WEBSOCKETS)
        if isopen(w)
            count += 1
        end
    end
    count
end


text:
const Tag = SmallString
const maxTagLength = smallLength
const MTagStorage = MSmallStringStorage # A mutable tag storage
const IntTag = IntSmallString  # An integer that can be cast to a Tag
const maxTags = 4
const TagSetStorage = SVector{maxTags,IntTag}
const MTagSetStorage = MVector{maxTags,IntTag}  # A mutable tag storage

struct TagSet
  data::TagSetStorage
  length::Int
  function TagSet()
    ts = TagSetStorage(ntuple(_ -> IntTag(0),Val(4)))
    new(ts,0)
  end
  TagSet(tags::TagSetStorage,len::Int) = new(tags,len)
end

TagSet(ts::TagSet) = ts

not(ts::Union{AbstractString,TagSet}) = Not(TagSet(ts))

function _hastag(ts::MTagSetStorage, ntags::Int, tag::IntTag)
  for n = 1:ntags
    @inbounds ts[n] == tag && return true
  end
  return false
end

function _addtag_ordered!(ts::MTagSetStorage, ntags::Int, tag::IntTag)
  if ntags == 0 || tag > @inbounds ts[ntags]
    @inbounds setindex!(ts,tag,ntags+1)
  else
    # check for repeated tags
    _hastag(ts,ntags,tag) && return ntags
    pos = ntags+1   # position new tag should go
    while pos > 1 && tag < @inbounds ts[pos-1]
      pos -= 1
      @inbounds setindex!(ts,ts[pos],pos+1)
    end
    @inbounds setindex!(ts,tag,pos)
  end
  return ntags+1
end

function _addtag!(ts::MTagSetStorage, ntags::Int, tag::IntTag)
  t = Tag(tag)
  if !isnull(t)
    if isint(t)
      error("Cannot use a bare integer as a tag.")
    else
      ntags = _addtag_ordered!(ts, ntags,tag)
    end
  end
  return ntags
end

function reset!(v::MTagStorage, nchar::Int)
  for i = 1:nchar
    @inbounds v[i] = IntChar(0)
  end
end

function TagSet(str::AbstractString)
  # Mutable fixed-size vector as temporary Tag storage
  current_tag = MTagStorage(ntuple(_ -> IntChar(0),Val(maxTagLength)))
  # Mutable fixed-size vector as temporary TagSet storage
  ts = MTagSetStorage(ntuple(_ -> IntTag(0),Val(maxTags)))
  nchar = 0
  ntags = 0
  for current_char in str
    if current_char == ','
      if nchar != 0
        ntags = _addtag!(ts,ntags,cast_to_uint(current_tag))
        # Reset the current tag
        reset!(current_tag,nchar)
        nchar = 0
      end
    elseif current_char != ' ' # TagSet constructor ignores whitespace
      nchar == maxTagLength && error("Currently, tags can only have up to $maxTagLength characters")
      nchar += 1
      @inbounds current_tag[nchar] = current_char
    end
  end
  # Store the final tag
  if nchar != 0
    ntags = _addtag!(ts,ntags,cast_to_uint(current_tag))
  end
  return TagSet(TagSetStorage(ts),ntags)
end

Base.convert(::Type{TagSet}, str::String) = TagSet(str)

"""
ITensors.data(T::TagSet)

Get the raw storage of the TagSet.

This is an internal function, please inform the
developers of ITensors.jl if there is functionality
you would like for TagSet that is not currently
available.
"""
data(T::TagSet) = T.data

Base.length(T::TagSet) = T.length
Base.getindex(T::TagSet,n::Int) = Tag(getindex(data(T),n))
Base.copy(ts::TagSet) = TagSet(data(ts),length(ts))

function Base.:(==)(ts1::TagSet,ts2::TagSet)
  l1 = length(ts1)
  l2 = length(ts2)
  l1 != l2 && return false
  for n in 1:l1
    @inbounds ts1[n] != ts2[n] && return false
  end
  return true
end

function hastag(ts::TagSet, tag)
  l = length(ts)
  l < 1 && return false
  for n = 1:l
    @inbounds Tag(tag) == ts[n] && return true
  end
  return false
end

function hastags(ts2::TagSet, tags1)
  ts1 = TagSet(tags1)
  l1 = length(ts1)
  l2 = length(ts2)
  l1 > l2 && return false
  for n1 = 1:l1
    @inbounds !hastag(ts2,ts1[n1]) && return false
  end
  return true
end

function addtags(ts::TagSet, tagsadd)
  if length(ts) == maxTags
    throw(ErrorException("Cannot add tag: TagSet already maximum size"))
  end
  tsadd = TagSet(tagsadd)
  res_ts = MVector(data(ts))
  ntags = length(ts)
  for n = 1:length(tsadd)
    @inbounds ntags = _addtag_ordered!(res_ts, ntags,IntSmallString(tsadd[n]))
  end
  return TagSet(TagSetStorage(res_ts),ntags)
end

function _removetag!(ts::MTagSetStorage, ntags::Int, t::Tag)
  for n = 1:ntags
    if @inbounds Tag(ts[n]) == t
      for j = n:ntags-1
        @inbounds ts[j] = ts[j+1]
      end
      @inbounds ts[ntags] = IntTag(0)
      return ntags -= 1
    end
  end
  return ntags
end

#TODO: optimize this function
function removetags(ts::TagSet, tagsremove)
  tsremove = TagSet(tagsremove)
  res_ts = MVector(data(ts))
  ntags = length(ts)
  for n=1:length(tsremove)
    @inbounds ntags = _removetag!(res_ts, ntags, tsremove[n])
  end
  return TagSet(TagSetStorage(res_ts),ntags)
end

#TODO: optimize this function
function replacetags(ts::TagSet, tagsremove, tagsadd)
  tsremove = TagSet(tagsremove)
  tsadd = TagSet(tagsadd)
  res_ts = MVector(data(ts))
  ntags = length(ts)
  # The TagSet must have the tags to be replaced
  !hastags(ts,tsremove) && return ts
  for n = 1:length(tsremove)
    @inbounds ntags = _removetag!(res_ts, ntags, tsremove[n])
  end
  for n = 1:length(tsadd)
    @inbounds ntags = _addtag_ordered!(res_ts, ntags,IntSmallString(tsadd[n]))
  end
  return TagSet(TagSetStorage(res_ts),ntags)
end

function tagstring(T::TagSet)
  res = ""
  N = length(T)
  N == 0 && return res
  for n=1:N-1
    res *= "$(Tag(T[n])),"
  end
  res *= "$(Tag(T[N]))"
  return res
end

function Base.show(io::IO, T::TagSet)
  print(io, "\"$(tagstring(T))\"")
end

function readcpp(io::IO,::Type{TagSet}; kwargs...)
  format = get(kwargs,:format,"v3")
  ts = TagSet()
  if format=="v3"
    mstore = MTagSetStorage(ntuple(_ -> IntTag(0),Val(maxTags)))
    ntags = 0
    for n=1:4
      t = readcpp(io,Tag;kwargs...)
      if t != Tag()
        ntags = _addtag_ordered!(mstore,ntags,IntSmallString(t))
      end
    end
    ts = TagSet(TagSetStorage(mstore),ntags)
  else
    throw(ArgumentError("read TagSet: format=$format not supported"))
  end
  return ts
end

function HDF5.write(parent::Union{HDF5File,HDF5Group},
                    name::AbstractString,
                    T::TagSet)
  g = g_create(parent,name)
  attrs(g)["type"] = "TagSet"
  attrs(g)["version"] = 1
  write(g,"tags", tagstring(T))
end

function HDF5.read(parent::Union{HDF5File,HDF5Group},
                   name::AbstractString,
                   ::Type{TagSet})
  g = g_open(parent,name)
  if read(attrs(g)["type"]) != "TagSet"
    error("HDF5 group '$name' does not contain TagSet data")
  end
  tstring = read(g,"tags")
  return TagSet(tstring)
end


text:
using MuJoCo
using LinearAlgebra

struct MPPIStrategy <: LearningStrategy
   costs::Vector{Float64}
   weights::Vector{Float64}

   Σ::Matrix{Float64}
   λ::Float64
   γ::Float64

   f::FunctionSet
   cfunc::Function
   rollout::Function
   c::Vector{Float64}
   o::Vector{Float64}
   function MPPIStrategy(model, sigma, lambda, gamma, funcset,
                         rollout=mjWrap.rollout2)
      dtype   = Float64
      K       = model.samples.K

      costs   = Vector{dtype}(undef, K)
      weights = Vector{dtype}(undef, K)

      new(costs, weights, sigma, lambda, gamma,
          funcset,
          (x...)->nothing,
          rollout,
          zeros(model.mjsys.nu), zeros(model.mjsys.ns))
   end
end

Base.show(io::IO, s::MPPIStrategy) = print(io, "MPPI")

function apply(f::Function, costs::Vector, states::Array{Float64,3})
   for k=1:size(states,3)
      costs[k] -= f(@view states[:,end,k])
   end
   #costs .-= f(states[:,end,:]) # block calculate
end

function setup!(s::MPPIStrategy, model::TrajOptModel)
   samples = model.samples

   # traj and samples both start from same place; send inital state to TrajOptModel
   if isnan(model.s0[1])
      @info("Setting initial state from experiment.")
      s.f.initstate!(model.mjsys, model.traj)
      model.s0 .= model.traj.state[:,1,1]
   else
      model.traj.state[:,1] .= model.s0
   end
   samples.state[:,1,:] .= model.s0

   push!(model.trace, :maxR => Vector{Float64}())
end

function update!(model::TrajOptModel, s::MPPIStrategy, iter, null)
   samples = model.samples
   T       = samples.T
   K       = samples.K
   mjsys   = model.mjsys

   #randn!(samples.ctrl)
   #rmul!(samples.ctrl, s.Σ[1]) # scale by one value
   #broadcast!(+, samples.ctrl, samples.ctrl, model.theta)
   #samples.ctrl[:,:,1] .= model.theta # have at least one thing to do a full theta rollout

   #samples.ctrl .= 0.0
   randn!( @view(samples.ctrl[:,(1+model.delay):end,:]) )
   rmul!(samples.ctrl, s.Σ[1]) # scale by one value
   broadcast!(+, samples.ctrl, samples.ctrl, model.theta)
   samples.ctrl[:,:,1] = model.theta # have at least one thing to do a full theta rollout

   mjWrap.limitcontrols(mjsys, samples)

   if iter == 1
      s.rollout(mjsys, samples,
                      s.cfunc, s.f.observe!, s.f.reward)
   else
      s.rollout(mjsys, samples,
                      s.cfunc, s.f.observe!, s.f.reward;
                      startT=1+model.delay)
   end
end

function applyweights!(theta::Matrix{Float64},
                       ctrl::Array{Float64,3}, weights::Vector{Float64})
   nu, T, K = size(ctrl)

   @inbounds for t=1:T
      #theta[:,t] = ctrl[:,t,:]*weights
      #@views A_mul_B!(theta[:,t], ctrl[:,t,:], weights)
      theta[:,t] .= 0.0
      for k=1:K
         for i=1:nu
            theta[i,t] += ctrl[i,t,k]*weights[k]
         end
      end
   end
end

function hook(s::MPPIStrategy, model::TrajOptModel, iter)

   LinearAlgebra.BLAS.set_num_threads(Threads.nthreads())
   BLAS.set_num_threads( 0 )

   mjsys   = model.mjsys
   samples = model.samples
   T       = samples.T
   K       = samples.K

   s.costs[:] = -1.0*Common.Rtau(samples.reward, s.γ)
   values = model.terminalvalue( samples.obs[:,end,:] )
   if values != nothing
      s.costs .-= values
   end

   baseline = minimum(s.costs)

   push!(model.trace[:stocR], mean(s.costs))
   push!(model.trace[:meanR], baseline)
   push!(model.trace[:maxR], maximum(s.costs))

   @. s.costs = exp((-(s.costs-baseline)/s.λ))
   η = sum(s.costs)
   @. s.weights = s.costs/η

   applyweights!(model.theta, samples.ctrl, s.weights)

   samples.ctrl[:,:,1] = model.theta # for visualization purposes

   if model.MPC && iter <= model.niter # get first controls, mj.step, store in model.traj
      traj    = model.traj

      # MPC: apply first controls to system; shift the rest around
      s.c[:] = model.theta[:,1]
      model.theta[:,1:end-1] .= @view model.theta[:,2:end]
      model.theta[:,end] .= 0.0

      # TODO
      # TODO
      # think about skip and how it should affect mpc mode
      # TODO
      # TODO
      mjWrap.reset(mjsys, model.s0, s.c, s.o)
      #for i=1:mjsys.skip
         mjWrap.step!(mjsys, s.c, model.s0, s.o) # get next state, next observation
      #end
      s.f.observe!(mjsys, model.s0, s.o) # arbitrary observation vector manipulation

      # advance our main trajectory
      traj.state[:,iter] = model.s0
      traj.ctrl[:,iter]  = s.c
      traj.obs[:,iter]   = s.o

      # update the mpc bundle of trajectories to start 'here'
      #for k=1:K
      #   samples.state[:,1,k] = model.s0
      #end
      samples.state[:,1+model.delay,:] .= model.s0
   end
end

function cleanup!(s::MPPIStrategy, model::TrajOptModel) # just rollout and store in traj
   if model.MPC == false
      #info("Running rollout for main trajectory.")
      mjsys   = model.mjsys
      traj    = model.traj

      model.traj.state .= model.s0
      traj.ctrl[:,:,1] = model.theta
      s.rollout(mjsys, traj, s.cfunc, s.f.observe!, s.f.reward)
   end
end


text:
import Base.:+
import Base.:*
import Base.:!

export Logic
export Proposition
export True
export False
export T
export F
export ↑
export ⟹
export ↔

## Data Structures ##

"""
    Represents a logical object.
"""
abstract type Logic end

"""
    Represents the value True.

field: value.
"""
struct True <: Logic
    value::True
    True() = new()
end

"""
    Represents the value False.

field: value.
"""
struct False <: Logic
    value::False
    False() = new()
end

"""
    Represents a proposition.

A variable that can take the value `True` or `False`. No others.

field: value.
"""
mutable struct Proposition <: Logic
    value::Logic
    Proposition() = new()
    function Proposition(p::Logic)
        instance = new()
        instance.value = p
        instance
    end
end

## Global Constants ##

const T = True()
const F = False()

## Generic Functions ##

"""
    value(p)

Return the logical value of `p`.
"""
value(p::True) = T
value(p::False) = F
value(p::Logic) = value(p.value)

"""
    show(p)

Print a string representation of the logical value `p`.
"""
Base.show(io::IO, p::True) = print(io, 'T')
Base.show(io::IO, p::False) = print(io, 'F')
Base.show(io::IO, p::Logic) = show(io, value(p))

## Binary Operators ##

"""
    p ↑ q

Not both. Alternative denial. Logical NAND.
"""
↑(p::True, q::True) = F
↑(p::True, q::False) = T
↑(p::False, q::True) = T
↑(p::False, q::False) = T
↑(p::Logic, q::Logic) = value(p) ↑ value(q)

"""
    p * q

Logical conjunction (AND).
"""
*(p::Logic, q::Logic) = (p ↑ q) ↑ (p ↑ q)

"""
    p + q

Logical disjunction (OR).
"""
+(p::Logic, q::Logic) = (p ↑ p) ↑ (q ↑ q)

"""
    p ⟹ q

Logical implication.
"""
⟹(p::Logic, q::Logic) = p ↑ (q ↑ q)

"""
    ↔(p, q)

Logical equality (=).
"""
↔(p::Logic, q::Logic) = (p ↑ q) ↑ ((p ↑ p) ↑ (q ↑ q))

## Unary Operators ##

"""
    !(p)

Negate the proposition `p`. Logical negation.
"""
!(p::Logic) = p ↑ p

"""
    +(p)

Returns a proposition identical to `p`. Logical identity.
"""
+(p::Logic) = p

"""
    True(p)

Returns a Proposition that is True for all `p`. Logical true.
"""
True(p::Logic) = T

"""
    False(p)

Returns a Proposition that is False for all `p`. Logical false.
"""
False(p::Logic) = F

max_issues_repo_name: iamazadi/Porta

max_forks_repo_name: iamazadi/Porta


text:
using OrdinaryDiffEq, DiffEqDevTools, DiffEqProblemLibrary, DiffEqBase

## Setup Tests

prob = prob_ode_linear

probs = Vector{ODETestProblem}(2)
probs[1] = prob
probs[2] = prob_ode_2Dlinear

tspan = [0,1]
tspans = Vector{Vector{Int64}}(2)
tspans[1] = tspan; tspans[2] = tspan
abstols = 1./10.^(3:10)
reltols = 1./10.^(3:10)

setups = [Dict(:alg=>RK4());Dict(:alg=>Euler());Dict(:alg=>BS3());
          Dict(:alg=>Midpoint());Dict(:alg=>BS5());Dict(:alg=>DP5())]

t1 = @elapsed sol = solve(prob,RK4(),dt=1/2^(4))
t2 = @elapsed sol2 = solve(prob,setups[1][:alg],dt=1/2^(4))

@test (sol2[end] == sol[end])

## Shootout Tests

println("Shootout Tests")

shoot = ode_shootout(prob,setups,dt=1/2^(4))

#show(shoot)
#println(shoot)
shoot[end]

set = ode_shootoutset(probs,setups;dt=1/2^(4))

#println(set[1])
#println(set[:])
set[end]
set[1][:]

## WorkPrecision Tests
println("WorkPrecision Tests")
wp = ode_workprecision(prob,DP5(),abstols,reltols;name="Dormand-Prince 4/5")

wp[1]
wp[:]
wp[end]
#show(wp)

wp_set = ode_workprecision_set(prob,abstols,reltols,setups;dt=1/2^4,numruns=2)

wp_set[1]
wp_set[:]
wp_set[end]
#println(wp_set)
#show(wp_set)
@test (minimum(diff(wp_set[1].errors).==0)) # The errors for a fixed timestep method should be constant

max_stars_repo_name: davidmuaddib/DiffEqDevTools

max_issues_repo_name: davidmuaddib/DiffEqDevTools

max_forks_repo_name: davidmuaddib/DiffEqDevTools


text:
#############
## NormIID ##
#############

type NormIID <: AbstrUnivarModel
    μ::Float64
    σ::Float64

    function NormIID(μ::Float64, σ::Float64)
        if σ <= 0
            error("σ parameter must be larger than 0")
        end
        return new(μ, σ)
    end
end

function NormIID(params::Array{Float64, 1})
    return NormIID(params...)
end

import Base.Multimedia.display
function display(mod::NormIID)
    nDigits = 3
    println("\nIID, constant univariate normal distribution:")
    println("    μ: $(round(mod.μ, nDigits))")
    println("    σ: $(round(mod.σ, nDigits))")
end

function description(mod::NormIID)
    return "IID, normally distributed"
end

## getParams
##----------

function getParams(mod::NormIID)
    return (mod.μ, mod.σ)
end

## simulate
##---------

function simulate(mod::NormIID, nObs::Int)
    ## simulate from NormIID model
    return randn(nObs)*mod.σ + mod.μ
end

## estimate
##---------

function estimate(dt::Type{NormIID}, data::Array{Float64, 1})
    normFit = fit_mle(Normal, data)
    muHat, sigmaHat = params(normFit)
    return NormIID(muHat, sigmaHat)
end

## resimulate
##-----------

@doc doc"""
Fit model and resimulate data from it.
"""->
function resimulate(data::AbstractTimenum,
                    dt::Type{Econometrics.NormIID},
                    nPaths::Int = 1)
    ## get returns without NAs
    pureRets = dropna(data.vals[1])
    nRets = size(pureRets, 1)

    ## fit model to pure returns
    modFit = estimate(dt, pureRets)

    return resimulate(data, modFit, nPaths)
end


text:
using Statistics
using NCDatasets

using ArgParse
using JSON
using Formatting

include("LinearRegression.jl")
include("nanop.jl")
include("CESMReader.jl")

using .CESMReader

correlation = (x1, x2) -> x1' * x2 / (sum(x1.^2)*sum(x2.^2)).^0.5

function parse_commandline()

    s = ArgParseSettings()
    @add_arg_table s begin

        "--data-file-prefix"
            help = "Data filename prefix including folder and path until the timestamp. File extension `nc` is assumed."
            arg_type = String
            required = true

        "--data-file-timestamp-form"
            help = "Data filename timestamp form. Either `YEAR` or `YEAR_MONTH`."
            arg_type = String
            required = true

        "--domain-file"
            help = "Domain file."
            arg_type = String
            required = true

        "--output-file"
            help = "Output file."
            arg_type = String
            required = true

        "--beg-year"
            help = "Year of begin."
            arg_type = Int64
            required = true

        "--end-year"
            help = "Year of end."
            arg_type = Int64
            required = true

        "--varname"
            help = "Variable name. If it is 2D then it should be (lon, lat, time), if it is 3D then it should be (lon, lat, lev/depth, time)"
            arg_type = String
            required = true

        "--no-detrend"
            help = "Whether to detrend the data or not."
            action = :store_true

        "--output-monthly-anomalies"
            help = "Whether to output the monthly anomalies for all months."
            arg_type = Bool
            default = false

        "--dims"
            help = "The form of dimensions. Now can be `XYZT`, `XYT`, `YZT`"
            arg_type = String
            required = true

    end

    return parse_args(ARGS, s)
end

parsed = parse_commandline()
print(json(parsed, 4))

output_file = parsed["output-file"]

if parsed["data-file-timestamp-form"] == "YEAR"
    filename_format = format("{:s}{{:04d}}.nc", joinpath(parsed["data-file-prefix"]))
    form = :YEAR
elseif parsed["data-file-timestamp-form"] == "YEAR_MONTH"
    filename_format = format("{:s}{{:04d}}-{{:02d}}.nc", joinpath(parsed["data-file-prefix"]))
    form = :YEAR_MONTH
end

Dataset(parsed["domain-file"], "r") do ds
    global mask = replace(ds["mask"], missing=>NaN)
end

if ! ( parsed["dims"] in ["XYZT", "XYT", "YZT", "YZ"] )
    ErrorException("Unknown dimension: " * parsed["dims"]) |> throw
end

# Creating preknowledge of dimensions
Dataset(format(filename_format, parsed["beg-year"], 1), "r") do ds

    var = ds[parsed["varname"]]
    dims = size(var)

    global beg_t = (parsed["beg-year"] - 1) * 12 + 1
    global end_t = (parsed["end-year"] - 1) * 12 + 12

    if parsed["dims"] == "XYZT"
        global (Nx, Ny, Nz, _ ) = size(var)
    elseif parsed["dims"] == "XYT"
        global (Nx, Ny, _ ) = size(var)
        Nz = 1
    elseif parsed["dims"] == "YZT"
        global (Ny, Nz, _ ) = size(var)
        Nx = 1
    elseif parsed["dims"] == "YZ"
        global (Ny, Nz ) = size(var)
        Nx = 1
    end

    global Nt = end_t - beg_t + 1

    if mod(Nt, 12) != 0
        ErrorException("Time record is not multiple of 12") |> throw
    end

    global nyears = Int64(Nt / 12)
    global Ns = Int64(Nt / 3)
    global Na = Int64(Nt / 12)
end

# M = monthly
# S = seasonal
# A = annual

data_MM    = zeros(Float64, Nx, Ny, Nz, 12)
data_MA    = zeros(Float64, Nx, Ny, Nz, Nt)
data_MAVAR = zeros(Float64, Nx, Ny, Nz, 12)
data_MASTD = zeros(Float64, Nx, Ny, Nz, 12)

data_SM    = zeros(Float64, Nx, Ny, Nz, 4)
data_SA    = zeros(Float64, Nx, Ny, Nz, Ns)
data_SAVAR = zeros(Float64, Nx, Ny, Nz, 4)
data_SASTD = zeros(Float64, Nx, Ny, Nz, 4)

data_AM    = zeros(Float64, Nx, Ny, Nz, 1)
data_AA    = zeros(Float64, Nx, Ny, Nz, Na)
data_AAVAR = zeros(Float64, Nx, Ny, Nz, 1)
data_AASTD = zeros(Float64, Nx, Ny, Nz, 1)

data_ZONAL_MM    = zeros(Float64, Ny, Nz, 12)
data_ZONAL_MA    = zeros(Float64, Ny, Nz, Nt)
data_ZONAL_MAVAR = zeros(Float64, Ny, Nz, 12)
data_ZONAL_MASTD = zeros(Float64, Ny, Nz, 12)

data_ZONAL_SM    = zeros(Float64, Ny, Nz, 4)
data_ZONAL_SA    = zeros(Float64, Ny, Nz, Ns)
data_ZONAL_SAVAR = zeros(Float64, Ny, Nz, 4)
data_ZONAL_SASTD = zeros(Float64, Ny, Nz, 4)

data_ZONAL_AM    = zeros(Float64, Ny, Nz, 1)
data_ZONAL_AA    = zeros(Float64, Ny, Nz, Na)
data_ZONAL_AAVAR = zeros(Float64, Ny, Nz, 1)
data_ZONAL_AASTD = zeros(Float64, Ny, Nz, 1)

months  = collect(Float64, 1:Nt)
seasons = collect(Float64, 1:Ns)
years   = collect(Float64, 1:Na)

fh = FileHandler(filename_format=filename_format, form=form)

function doit!(
    t    :: AbstractArray{Float64, 1},
    d     :: AbstractArray{Float64, 1},
    M    :: AbstractArray{Float64},
    A    :: AbstractArray{Float64},
    AVAR :: AbstractArray{Float64},
    ASTD :: AbstractArray{Float64},
    nyears :: Integer,
    period :: Integer,
    idx  :: Tuple,
)

    if ! parsed["no-detrend"]
        d = detrend(t, d, order=2)
    end

    M[idx..., :] = mean( reshape(d, period, :), dims=2 )[:, 1]
    A[idx..., :] = d - repeat( M[idx..., :], outer=nyears)

    for m = 1:period
        d_yy = view(A, idx..., m:period:(m+nyears*period-1))
        AVAR[idx..., m] = var(d_yy)
        ASTD[idx..., m] = std(d_yy)
    end

end

if parsed["dims"] == "XYZT"

    for k=1:Nz
        spatial_rng = (:, :, k)
        global data = reshape( getData(fh, parsed["varname"], (parsed["beg-year"], parsed["end-year"]), spatial_rng), Nx, Ny, Nt)

        for i=1:Nx, j=1:Ny
            d = view(data, i, j, :)
            doit!(months,           d, data_MM, data_MA, data_MAVAR, data_MASTD, nyears, 12, i, j, k)

            seasonal_d = mean(reshape( circshift(d, -2), 3, :), dims=1)[1, :] # spring, summer, autumn, winter
            doit!(seasons, seasonal_d, data_SM, data_SA, data_SAVAR, data_SASTD, nyears,  4, i, j, k)

            annual_d = mean(reshape( d, 12, :), dims=1)[1, :]
            doit!(years, annual_d, data_AM, data_AA, data_AAVAR, data_AASTD, nyears,  1, i, j, k)
        end
    end

elseif  parsed["dims"] == "XYT"

    spatial_rng = (:, :)
    global data = reshape( getData(fh, parsed["varname"], (parsed["beg-year"], parsed["end-year"]), spatial_rng), Nx, Ny, Nt)
    for i=1:Nx, j=1:Ny

        d = view(data, i, j, :)
        doit!(months,           d, data_MM, data_MA, data_MAVAR, data_MASTD, nyears, 12, (i, j, 1) )

        seasonal_d = mean(reshape( circshift(d, -2), 3, :), dims=1)[1, :]
        doit!(seasons, seasonal_d, data_SM, data_SA, data_SAVAR, data_SASTD, nyears,  4, (i, j, 1) )

        annual_d = mean(reshape( d, 12, :), dims=1)[1, :]
        doit!(years, annual_d, data_AM, data_AA, data_AAVAR, data_AASTD,     nyears,  1, (i, j, 1) )

    end

    # === zonal mean ====
    data_zonal = nanmean( data, dims=(1,) )
    for j=1:Ny

        monthly_d = view(data_zonal, 1, j, :)
        doit!(months, monthly_d, data_ZONAL_MM, data_ZONAL_MA, data_ZONAL_MAVAR, data_ZONAL_MASTD, nyears, 12, (j, 1))

        seasonal_d = mean(reshape( circshift(monthly_d, -2), 3, :), dims=1)[1, :]
        doit!(seasons, seasonal_d, data_ZONAL_SM, data_ZONAL_SA, data_ZONAL_SAVAR, data_ZONAL_SASTD, nyears, 4, (j, 1))

        annual_d = mean(reshape( monthly_d, 12, : ), dims=1)[1, :]
        doit!(years,   annual_d, data_ZONAL_AM, data_ZONAL_AA, data_ZONAL_AAVAR, data_ZONAL_AASTD, nyears, 1, (j, 1))

    end

elseif parsed["dims"] in ["YZT", "YZ"]

    spatial_rng = (:, :)
    global data = reshape( getData(fh, parsed["varname"], (parsed["beg-year"], parsed["end-year"]), spatial_rng), Ny, Nz, Nt)
    for j=1:Ny, k=1:Nz

        d = view(data, j, k, :)
        doit!(months,           d, data_MM, data_MA, data_MAVAR, data_MASTD, nyears, 12, (1, j, k)    )

        seasonal_d = mean(reshape( circshift(d, -2 ), 3, :), dims=1)[1, :]
        doit!(seasons, seasonal_d, data_SM, data_SA, data_SAVAR, data_SASTD, nyears,  4, (1, j, k)    )

        annual_d = mean(reshape( d, 12, :), dims=1)[1, :]
        doit!(years, annual_d, data_AM, data_AA, data_AAVAR, data_AASTD,     nyears,  1, (1, j, k)    )

    end

    data_ZONAL_MM .= data_MM[1, :, :, :]
    data_ZONAL_MA .= data_MA[1, :, :, :]
    data_ZONAL_MAVAR .= data_MAVAR[1, :, :, :]
    data_ZONAL_MASTD .= data_MASTD[1, :, :, :]

    data_ZONAL_SM .= data_SM[1, :, :, :]
    data_ZONAL_SA .= data_SA[1, :, :, :]
    data_ZONAL_SAVAR .= data_SAVAR[1, :, :, :]
    data_ZONAL_SASTD .= data_SASTD[1, :, :, :]

    data_ZONAL_AM .= data_AM[1, :, :, :]
    data_ZONAL_AA .= data_AA[1, :, :, :]
    data_ZONAL_AAVAR .= data_AAVAR[1, :, :, :]
    data_ZONAL_AASTD .= data_AASTD[1, :, :, :]

end

Dataset(output_file, "c") do ds

    defDim(ds, "months", 12)
    defDim(ds, "seasons", 4)
    defDim(ds, "years", 1)
    defDim(ds, "time", Nt)
    defDim(ds, "Nx", Nx)
    defDim(ds, "Ny", Ny)
    defDim(ds, "Nz", Nz)

    datas =  convert(Array{Any}, [

        (format("{:s}_MM",    parsed["varname"]),       data_MM,    ("Nx", "Ny", "Nz", "months"), Dict()),
        (format("{:s}_MAVAR", parsed["varname"]),       data_MAVAR, ("Nx", "Ny", "Nz", "months"), Dict()),
        (format("{:s}_MASTD", parsed["varname"]),       data_MASTD, ("Nx", "Ny", "Nz", "months"), Dict()),
        (format("{:s}_ZONAL_MM", parsed["varname"]),    data_ZONAL_MM, ("Ny", "Nz", "months"), Dict()),
        (format("{:s}_ZONAL_MAVAR", parsed["varname"]), data_ZONAL_MAVAR, ("Ny", "Nz", "months"), Dict()),
        (format("{:s}_ZONAL_MASTD", parsed["varname"]), data_ZONAL_MASTD, ("Ny", "Nz", "months"), Dict()),

        (format("{:s}_SM",    parsed["varname"]),       data_SM,    ("Nx", "Ny", "Nz", "seasons"), Dict()),
        (format("{:s}_SAVAR", parsed["varname"]),       data_SAVAR, ("Nx", "Ny", "Nz", "seasons"), Dict()),
        (format("{:s}_SASTD", parsed["varname"]),       data_SASTD, ("Nx", "Ny", "Nz", "seasons"), Dict()),
        (format("{:s}_ZONAL_SM", parsed["varname"]),    data_ZONAL_SM, ("Ny", "Nz", "seasons"), Dict()),
        (format("{:s}_ZONAL_SAVAR", parsed["varname"]), data_ZONAL_SAVAR, ("Ny", "Nz", "seasons"), Dict()),
        (format("{:s}_ZONAL_SASTD", parsed["varname"]), data_ZONAL_SASTD, ("Ny", "Nz", "seasons"), Dict()),

        (format("{:s}_AM",    parsed["varname"]),       data_AM,    ("Nx", "Ny", "Nz", "years"), Dict()),
        (format("{:s}_AAVAR", parsed["varname"]),       data_AAVAR, ("Nx", "Ny", "Nz", "years"), Dict()),
        (format("{:s}_AASTD", parsed["varname"]),       data_AASTD, ("Nx", "Ny", "Nz", "years"), Dict()),
        (format("{:s}_ZONAL_AM", parsed["varname"]),    data_ZONAL_AM, ("Ny", "Nz", "years"), Dict()),
        (format("{:s}_ZONAL_AAVAR", parsed["varname"]), data_ZONAL_AAVAR, ("Ny", "Nz", "years"), Dict()),
        (format("{:s}_ZONAL_AASTD", parsed["varname"]), data_ZONAL_AASTD, ("Ny", "Nz", "years"), Dict()),

    ])

    if parsed["output-monthly-anomalies"]
        push!(datas, (format("{:s}_MA",    parsed["varname"]),       data_MA,    ("Nx", "Ny", "Nz", "time"),   Dict()) )
    end

    if  parsed["dims"] in [ "XYT", "XYZT" ]
        push!(datas, ("mask",  mask,    ("Nx", "Ny"),   Dict()))
    end

    for (varname, vardata, vardim, attrib) in datas
        if ! haskey(ds, varname)
            var = defVar(ds, varname, Float64, vardim)
            var.attrib["_FillValue"] = 1e20
        end

        println("Writing variable:  ", varname)
        var = ds[varname]

        for (k, v) in attrib
            var.attrib[k] = v
        end

        rng = []
        for i in 1:length(vardim)-1
            push!(rng, Colon())
        end
        push!(rng, 1:size(vardata)[end])
        var[rng...] = vardata

    end

end

max_stars_repo_name: meteorologytoday/CESM-diagnostic

max_issues_repo_name: meteorologytoday/CESM-diagnostic

max_forks_repo_name: meteorologytoday/CESM-diagnostic


text:
# Deprecations v0.1
@deprecate stack(x, dims) stack(x; dims=dims)
@deprecate unstack(x, dims) unstack(x; dims=dims)
@deprecate unsqueeze(x::AbstractArray, dims::Int) unsqueeze(x; dims=dims)
@deprecate unsqueeze(dims::Int) unsqueeze(dims=dims)
@deprecate labelmap(x) group_indices(x)
@deprecate frequencies(x) group_counts(x)
@deprecate eachbatch(data, batchsize; kws...) eachobs(data; batchsize, kws...)
@deprecate eachbatch(data; size=1, kws...) eachobs(data; batchsize=size, kws...)


text:
module LakeModelParametersModule

struct LakeModelParameters
  instant_throughflow::Bool
  lake_retention_constant::Float64
  function LakeModelParameters()
    instant_throughflow::Bool = true
    lake_retention_constant::Float64 = 0.1
    return new(instant_throughflow,lake_retention_constant)
  end
end

end

max_stars_repo_name: ThomasRiddick/DynamicHD

max_issues_repo_name: ThomasRiddick/DynamicHD

max_forks_repo_name: ThomasRiddick/DynamicHD


text:
__precompile__() # this module is safe to precompile
module ODINN

###############################################
############       PACKAGES     ##############
###############################################

using Statistics, LinearAlgebra, Random, Polynomials
using JLD2
using OrdinaryDiffEq, DiffEqFlux
using Zygote: @ignore
using Flux
using Tullio, RecursiveArrayTools
using Infiltrator
using Plots
using Makie, CairoMakie
import Pkg
using Distributed
using ProgressMeter
using PyCall

###############################################
#############    PARAMETERS     ###############
###############################################

include("helpers/parameters.jl")

###############################################
#############  ODINN LIBRARIES  ###############
###############################################

cd(@__DIR__)
global root_dir = dirname(Base.current_project())
global root_plots = joinpath(root_dir, "plots")

include("helpers/utils.jl")

#### Plotting functions  ###
include("helpers/plotting.jl")
### Iceflow modelling functions  ###
# (includes utils.jl as well)
include("helpers/iceflow.jl")

###############################################
#############  PYTHON LIBRARIES  ##############
###############################################

const netCDF4 = PyNULL()
const cfg = PyNULL()
const utils = PyNULL()
const workflow = PyNULL()
const tasks = PyNULL()
const graphics = PyNULL()
const bedtopo = PyNULL()
const MBsandbox = PyNULL()

# Essential Python libraries
const np = PyNULL()
const xr = PyNULL()

###############################################
######### PYTHON JULIA INTERACTIONS  ##########
###############################################

include(joinpath(ODINN.root_dir, "src/helpers/config.jl"))
### Climate data processing  ###
include(joinpath(ODINN.root_dir, "src/helpers/climate.jl"))
### OGGM configuration settings  ###
include(joinpath(ODINN.root_dir, "src/helpers/oggm.jl"))

end # module


text:
using Test
using Random
include("../../../../examples/findminmax.jl")

@testset for seed in 1:100
    CUDA.seed!(seed)
    rng = MersenneTwister(seed)
    # ns = [1:10; [255, 256, 257, 511, 512, 513]; rand(11:1000, 100); ]
    ns = [1:1100; [10^4, 10^5, 10^6]]
    @testset for n in ns
        xs = CUDA.rand(n)
        xs_cpu = Array(xs)
        @test findminmax(xs) == findminmax(xs_cpu) == findminmax_base(xs_cpu)
    end
end

@testset "large" begin
    seed = 123456789
    CUDA.seed!(seed)
    rng = MersenneTwister(seed)
    @testset for n in [10^7, 10^8]
        xs = CUDA.rand(n)
        xs_cpu = Array(xs)
        @test findminmax(xs) == findminmax(xs_cpu) == findminmax_base(xs_cpu)
    end
end


text:
# Remove when package is registered
push!(LOAD_PATH, "../src")

using Optini
using Documenter

makedocs(;
    modules=[Optini],
    authors="Long Nguyen",
    repo="https://github.com/lhnguyen-vn/Optini.jl/blob/{commit}{path}#L{line}",
    sitename="Optini.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://lhnguyen-vn.github.io/Optini.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
        "Manual" => [
            "Univariate" => "man/univariate.md",
            "Multivariate" => [
                "man/multivariate/linesearch.md",
                "man/multivariate/trustregion.md"
            ]
        ],
        "tutorials.md"
    ],
)

deploydocs(;
    repo="github.com/lhnguyen-vn/Optini.jl.git",
    devbranch = "main"
)


text:
# tests for find_zeros

using Roots
using Test

mutable struct CallableFunction
    f
    n
end
(F::CallableFunction)(x) = (F.n += 1; F.f(x))

@testset "find_zeros" begin
    function azero(f, x)
        fx = f(x)
        iszero(fx) && return true
        sign(fx) * sign(f(nextfloat(x))) < 0 && return true
        sign(fx) * sign(f(prevfloat(x))) < 0 && return true
        abs(fx) <= 8eps(x) && return true
        false
    end

    ## easier ones, counting steps
    F = CallableFunction(x -> exp(x) - x^4, 0)
    xrts = find_zeros(F, -5, 20)
    @test length(xrts) == 3
    @test all(azero.((F,), xrts))
    @test F.n <= 3000

    F = CallableFunction(x -> cos(x) + cos(2x), 0)
    xrts = find_zeros(F, 0, 4pi)
    @test length(xrts) == 6
    @test all(azero.((F,), xrts))
    @test F.n <= 5000

    T11(x) = 1024x^11 - 2816x^9 + 2816x^7 - 1232x^5 + 220x^3 - 11x
    U9(x) = 512x^9 - 1024x^7 + 672x^5 - 160x^3 + 10x

    F = CallableFunction(T11, 0)
    xrts = find_zeros(F, -1, 1)
    @test length(xrts) == 11
    @test all(azero.((F,), xrts))
    @test F.n <= 10_000

    F = CallableFunction(U9, 0)
    xrts = find_zeros(F, -1, 1)
    @test length(xrts) == 9
    @test all(azero.((F,), xrts))
    @test F.n <= 10_000

    W(n) = x -> prod(x - i for i in 1:n)
    Wi(n) = x -> prod((x - i)^i for i in 1:n)

    F = CallableFunction(W(20), 0)
    xrts = find_zeros(F, -1, 21)
    @test length(xrts) == 20
    @test all(azero.((F,), xrts))
    @test F.n <= 20_000

    F = CallableFunction(Wi(6), 0)
    xrts = find_zeros(F, -1, 7)
    @test length(xrts) == 6
    @test all(azero.((F,), xrts))
    @test F.n <= 10_000

    ## Harder ones
    f1(x) = 2 * exp(0.5 * x) * (sin(5 * x) + sin(101 * x))
    tiger_tail(x) = f1(x) - round(f1(x)) # (-2,1) then filter

    f2(x) = (x - 0.5) * (x - 0.5001) * (x - 4) * (x - 4.05) * (x - 9.3) # (0,10)
    f3(x) = (x - 3)^2 * (x - 4)^2 # (0,5)
    delta = 0.001
    f4(x) = (x - 0.5)^3 * (x - (0.5 + delta)) * (x - 1)
    f5(x) = (x - 0.5)^3 * (x - (0.5 + delta))^3 * (x - 4) * (x - (4 + delta)) * (x - 4.2)^2
    M(n) = x -> prod((x - (0.5 - (1 / 10)^i)) for i in 1:n)
    f6 = M(4)
    f7 = M(5) # too much

    xrts = find_zeros(tiger_tail, -2.0, 1.0)
    xrts = filter(u -> -1 / 4 < tiger_tail(u) < 1 / 4, xrts)
    @test 345 <= length(xrts) <= 345
    @test all(azero.((tiger_tail,), xrts))

    xrts = find_zeros(f2, 0.0, 10.0)
    @test length(xrts) == 5
    @test all(azero.((f2,), xrts))

    xrts = find_zeros(f3, 0.0, 5.0)
    @test length(xrts) == 2
    @test all(azero.((f3,), xrts))

    xrts = find_zeros(f4, 0.0, 1.5)
    @test length(xrts) == 3
    @test all(azero.((f4,), xrts))

    xrts = find_zeros(f5, 0.0, 5.0)
    @test length(xrts) >= 3          # too hard to get 5 w/o luck, as with no_pts=21/k=4
    @test all(azero.((f5,), xrts))

    xrts = find_zeros(f6, 0.0, 10.0)
    @test length(xrts) == 4
    @test all(azero.((f6,), xrts))

    xrts = find_zeros(f7, 0.0, 10.0)  # too sensitive to interval
    @test length(xrts) >= 3           # should be 5
    @test all(azero.((f7,), xrts))

    # Issue #141 solve over [a,b], not (a,b)
    @test length(find_zeros(p -> p * (1 - p), 0, 1)) == 2
    @test length(find_zeros(sin, 0, 5pi)) == 5 + 1

    # test different ways to specify an interval (just need `extrema` defined)
    @test find_zeros(sin, (3, 4)) ≈ [float(pi)]
    @test find_zeros(sin, [3, 4]) ≈ [float(pi)]
    @test find_zeros(sin, 3:4) ≈ [float(pi)]
    @test find_zeros(sin, SomeInterval(3, 4)) ≈ [float(pi)]
    @test find_zeros(sin, range(3, stop=4, length=20)) ≈ [float(pi)]

    # test with constant function
    @test isempty(find_zeros(x -> 4, -10, 10))
end

@testset "find_zeros: not Float64 types" begin
    for T in [Float16, Float32, BigFloat]
        rts = find_zeros(x -> cos(x) - x / 10, T(0.0), T(10.0))
        @test eltype(rts) == T
    end
end


text:
using BenchmarkTools
using SPH

function M6_hardcoded(q)
    if q < 0.0
        return 0.0
    elseif q < 1.0
        return (3 - q)^5 - 6*(2 - q)^5 + 15*(1 - q)^5
    elseif q < 2.0
        return (3 - q)^5 - 6*(2 - q)^5
    elseif q < 3.0
        return (3 - q)^5
    else
        return 0.0
    end
end

println("Kernel function benchmarks")
println("Hardcoded M6")
@btime $M6_hardcoded(x) setup=(x=3*rand())

println("Programmatic M6 (and 1st and second derivatives")
@btime $M6[1](x) setup=(x=3*rand())
@btime $(differentiate(M6[1]))(x) setup=(x=3*rand())
@btime $(differentiate(M6[1], 2))(x) setup=(x=3*rand())


text:
export rotate, rotate_corrs, rotate_corrs!, get_corr

"""
  rotate(corrpaths, ROTOUT)

Rotate all correlations from ENZ to RTZ in directory `CORROUT`.
# Arguments
- `corrpaths::Array{String,1}`: Array of filepaths for 9 component correlations.
- `ROTOUT::String`: Output directory for rotated correlations.
- `Overwrite::Bool`: Overwrite non-rotated correlations with rotated correlations.
"""
function rotate(corrpaths::Array{String,1}, ROTOUT::String; overwrite::Bool=false)

    # read correlations
    N = length(corrpaths)
    comps = Array{String}(undef,N)
    for ii = 1:N
        path, file = splitdir(corrpaths[ii])
        file = replace(file,".jld2"=>"")
        net1,sta1,chan1,net2,sta2,chan2 = split(file,'.')
        comps[ii] = chan1[end] * chan2[end]
    end

    # check if CORROUT exists and make new directory
    ROTOUT = expanduser(ROTOUT)
    isdir(ROTOUT) || mkpath(ROTOUT)

    # load data
    corrs = [load_corr(corrpaths[ii],comps[ii]) for ii = 1:N]
    azi,baz = corrs[1].azi, corrs[1].baz

    # do not rotate if no azi and baz in corrs
    if (azi == 0.) & (baz == 0.)
        println("No location data for $(corrs[1].name) $(corrs[1].id)")
    else
        # rotate and save correlations
        rotate_corrs!(corrs,azi,baz)
        for corr in corrs
            save_corr(corr,ROTOUT)
        end

        # delete non-rotated correlations
        if overwrite
            for ii = 1:N
                rm(corrpaths[ii])
            end
        end
    end

end

"""
    rotate_corrs(corrs,azi,baz)

Rotate cross-correlations from E,N,Z to R,T,Z coordinate system.
# Arguments
- `corrs::Array{CorrData,1}`: Array of EE, EN, EZ, NE, NN, NZ, ZE, ZN, ZZ CorrData
- `azi::Float64`: Azimuth (in radians) from station 1 to station 2
- `baz::Float64`: Back azimuth (in radians) from station 1 to station 2
"""
function rotate_corrs!(corrs::Array{CorrData,1},azi::Float64,baz::Float64)
    N = length(corrs)
    # get time intersect of all components
    ts = [c.t for c in corrs]
    inter = intersect(ts)[1]

    # subset so all corrs match in time
    for ii = 1:N
        ind = findall(x -> x ∈ inter, corrs[ii].t)
        corrs[ii].t = corrs[ii].t[ind]
        corrs[ii].corr = corrs[ii].corr[:,ind]
    end

    comp_dict = Dict()
    for ii = 1:N
        comp_dict[corrs[ii].comp] = ii
    end

    # get EE EN NN NE comps for rotation to TT RR TR RT
    # see Lin et al., 2008: Surface wave tomography of the western
    # United States from ambient seismic noise, GJI
    big_rot = [[-cos(azi)*cos(baz) cos(azi)*sin(baz) -sin(azi)*sin(baz) sin(azi)*cos(baz)];
               [-sin(azi)*sin(baz) -sin(azi)*cos(baz) -cos(azi)*cos(baz) -cos(azi)*sin(baz)];
               [-cos(azi)*sin(baz) -cos(azi)*cos(baz) sin(azi)*cos(baz) sin(azi)*sin(baz)];
               [-sin(azi)*cos(baz) sin(azi)*sin(baz) cos(azi)*sin(baz) -cos(azi)*cos(baz)]]
    to_rot = Array{eltype(corrs[1].corr)}(undef,4,size(corrs[1].corr)[1],length(inter))
    to_rot[1,:,:] = corrs[comp_dict["EE"]].corr
    to_rot[2,:,:] = corrs[comp_dict["EN"]].corr
    to_rot[3,:,:] = corrs[comp_dict["NN"]].corr
    to_rot[4,:,:] = corrs[comp_dict["NE"]].corr
    rotated = rotation_kernel(big_rot,to_rot)

    # update correlations from EE EN NN NE to TT RR TR RT
    update_corr!(corrs[comp_dict["EE"]],rotated[1,:,:],"TT")
    update_corr!(corrs[comp_dict["EN"]],rotated[2,:,:],"RR")
    update_corr!(corrs[comp_dict["NN"]],rotated[3,:,:],"TR")
    update_corr!(corrs[comp_dict["NE"]],rotated[3,:,:],"RT")

    # rotate from EZ NZ to RZ TZ
    azi_rot = [[cos(azi) -sin(azi)];
               [sin(azi) cos(azi)]]
    to_rot = Array{eltype(corrs[1].corr)}(undef,2,size(corrs[1].corr)[1],length(inter))
    to_rot[1,:,:] = corrs[comp_dict["NZ"]].corr
    to_rot[2,:,:] = corrs[comp_dict["EZ"]].corr
    rotated = rotation_kernel(azi_rot,to_rot)
    update_corr!(corrs[comp_dict["NZ"]],rotated[1,:,:],"RZ")
    update_corr!(corrs[comp_dict["EZ"]],rotated[2,:,:],"TZ")

    # rotate from ZE ZN to ZR ZT
    baz_rot = [[-cos(baz) -sin(baz)];
               [sin(baz) -cos(baz)]]
    to_rot = Array{eltype(corrs[1].corr)}(undef,2,size(corrs[1].corr)[1],length(inter))
    to_rot[1,:,:] = corrs[comp_dict["ZN"]].corr
    to_rot[2,:,:] = corrs[comp_dict["ZE"]].corr
    rotated = rotation_kernel(baz_rot,to_rot)
    update_corr!(corrs[comp_dict["ZN"]],rotated[1,:,:],"ZR")
    update_corr!(corrs[comp_dict["ZE"]],rotated[2,:,:],"ZT")

    return corrs
end
rotate_corrs(corrs::Array{CorrData,1},azi::Float64,baz::Float64) =
                      (U = deepcopy(corrs);rotate_corrs!(U,azi,baz);return U)

"""
  rotation_kernel(rot_matrix, to_rotate)

Rotation function for 2D correlation functions.

# Arguments
- `rot_matrix::Array`: Rotation matrix (either 2x2 or 4x4).
- `to_rotate::Array`: Array to be rotated column-wise.
"""
function rotation_kernel(rot_matrix::Array, to_rotate::Array)
    rotated = similar(to_rotate)
    for ii = 1:size(to_rotate,2)
        for jj = 1:size(to_rotate,3)
            rotated[:,ii,jj] .= rot_matrix * to_rotate[:,ii,jj]
        end
    end
    return rotated
end

"""
  update_corr!(C,comp)

Change component, name and corr of rotated correlation `C` to component `comp`.

# Arguments
- `C::CorrData`: CorrData object.
- `new_corr::Array`: Rotated correlation matrix.
- `comp::String`: New component after rotation, e.g. "RR"
"""
function update_corr!(C::CorrData,new_corr::Array,comp::String)
    # change comp of correlation
    C.comp = comp

    # change name of correlation
    newname = split(C.name,'.')
    newname[4] = newname[4][1:end-1] * comp[1]
    newname[8] = newname[8][1:end-1] * comp[2]
    C.name = join(newname,'.')

    # change rotated flag
    C.rotated = true

    # update correlation
    C.corr = new_corr
    return C
end

"""
  get_corr(CORROUT)

Get name of all correlations in directory `CORROUT` for rotation.

# Arguments
- `CORROUT::String`: Path to directory of correaltion files, e.g. /home/CORR/
"""
function get_corr(CORROUT::String; autocorr::Bool=false)
    files = glob("*",CORROUT)
    corrs = [splitdir(f)[2] for f in files]
    corrs = [replace(c,".jld2"=>"") for c in corrs]

    # get net1.sta1.net2.sta2, comp and if autocorrelation from name of
    # correlation jld2 file
    N = length(corrs)
    NETSTA = Array{String,1}(undef,N)
    COMP = Array{String,1}(undef,N)
    AUTO = Array{Bool,1}(undef,N)

    for ii = 1:N
        net1,sta1,chan1,net2,sta2,chan2 = split(corrs[ii],'.')
        NETSTA[ii] = join([net1,sta1,net2,sta2],'.')
        COMP[ii] = chan1[end] * chan2[end]
        AUTO[ii] = net1*sta1 == net2*sta2
    end

    # create dataframe
    df = DataFrame(NETSTA = NETSTA, COMP = COMP, AUTO=AUTO, FILE = files)

    # remove autocorrelations
    df = df[df[:AUTO] .== autocorr,:]

    return df2files(df)
end

"""
  df2files(df)

Helper function to group correlation pairs.

Returns an array of file paths for each correlation pair in dataframe.

# Arguments
- `df::DataFrame`: DataFrame with correlation pair name, filepath and component.
"""
function df2files(df::DataFrame)
    grouped = groupby(df,:NETSTA)
    N = length(grouped)
    out = Array{Array{String,1},1}(undef,N)
    for ii = 1:N
        out[ii] = collect(grouped[ii][:FILE])
    end
    return out
end


text:
#!/usr/bin/env julia

using Shadowsocks
using Sockets
using ArgParse
using FileWatching

function main(args)
    s = ArgParseSettings(description = "This is a julia implementation of shadowsocks",
        commands_are_required = false,
        version = "0.0.2",
        add_version = true)

    @add_arg_table s begin
        "--server", "-s"
            help = "shadowsocks server"
        "--port", "-p"
            help = "shadowsocks server port"
            arg_type = Int
        "--listen", "-l"
            help = "shadowsocks listen port"
            arg_type = Int
        "--method", "-m"
            help = "encryption method\nchacha20-poly1305-ietf, xchacha20-poly1305-ietf"
        "--password", "-k"
            help = "access password, like Pass1234"
        "--uri", "-c"
            help = "ss://chacha20-poly1305-ietf:imgk0000@127.0.0.1:8388\nor ss://chacha20-poly1305-ietf:imgk0000@:8388"
        "--config-file", "-f"
            help = "use config file"
    end

    config = nothing

    if args == String[]
        main(String["--help"])
    else
        parsed_args = parse_args(args, s)
        if parsed_args["config-file"] != nothing
            config = Base.Filesystem.abspath(parsed_args["config-file"])
        elseif parsed_args["uri"] != nothing
            config, err = Shadowsocks.parseURI(parsed_args["uri"])
            if err != nothing
                main(String["--help"])
            end
        else
            if haskey(parsed_args, "listen")
                config = Shadowsocks.SSClient(
                    parse(Sockets.IPAddr, parsed_args["server"]),
                    parsed_args["port"],
                    parsed_args["listen"],
                    uppercase(parsed_args["method"]),
                    parsed_args["password"]
                )
            else
                config = Shadowsocks.SSServer(
                    parse(Sockets.IPAddr, parsed_args["server"]),
                    parsed_args["port"],
                    uppercase(parsed_args["method"]),
                    parsed_args["password"]
                )
            end
        end
    end

    if config isa String
        if (local configs = Shadowsocks.readConfigFile(config); configs isa Array{Shadowsocks.SSConfig})
            Shadowsocks.@terminal("running shadowsocks client")
            Shadowsocks.run(configs)
        elseif configs isa Dict{String, Any}
            ch = Channel{Dict{String, Any}}(1)
            put!(ch, configs)

            @async while true # watch config file
                watch_file(config, -1)
                put!(ch, Shadowsocks.readConfigFile(config))
            end

            Shadowsocks.@terminal("running shadowsocks server")
            Shadowsocks.@terminal("watching config file - $config")
            Shadowsocks.run(ch)
        end
    elseif config isa Shadowsocks.SSConfig && config.lisPort == nothing
        Shadowsocks.@terminal("running shadowsocks server")
        Shadowsocks.runServer(config)
    elseif config isa Shadowsocks.SSConfig && config.lisPort != nothing
        Shadowsocks.@terminal("running shadowsocks client")
        Shadowsocks.runClient(config)
    end

    main(String["--help"])
end

main(ARGS)
# main(["-c", "ss://chacha20-poly1305-ietf:imgk0000@:8388"])
# main(["-c", "ss://chacha20-poly1305-ietf:imgk0000@0.0.0.0:8388"])


text:
## Program to run PPSeq on a section of data

# INPUTS:

# PPSeq_file = the PPSeq.jl file to use
# Data_directory = where to find the data to use
# Time = how long the segment is MUST BE FLOAT
# Neurons = how many neurons being used MUST BE INTEGER

if size(ARGS) != (4,)
    error("Not correct number of command line arguments")
end

PPSeq_file = ARGS[1]
Data_directory = ARGS[2]
max_time = parse(Float64, ARGS[3])
num_neurons = parse(Int64, ARGS[4])

# Import PPSeq
include(PPSeq_file)
const seq = PPSeq

# Other Imports
import DelimitedFiles: readdlm
import Random
import StatsBase: quantile
using CSV, DataFrames, Dates

function save_results(results, run_id)

    # ---------- anneal_latent_event_hist -----------
    latent_event_hist = results[:anneal_latent_event_hist]
    init_array = zeros(Float64, 0, 5)
    for (i, s) in enumerate(latent_event_hist)
        aux_array = zeros(Float64, length(s)+1, 5)
        for (j, z) in enumerate(s)
            aux_array[j, 1] = z.assignment_id
            aux_array[j, 2] = z.timestamp
            aux_array[j, 3] = z.seq_type
            aux_array[j, 4] = z.seq_warp
            aux_array[j, 5] = z.amplitude
        end
        aux_array[length(s)+1, :] = -1*ones(5)
        init_array = vcat(init_array, aux_array)
    end
    anneal_latent_event_hist_array = init_array
    ANNEAL_LATENT_EVENT_SHAPE = size(init_array)

    # ----------- latent_event_hist ------------
    latent_event_hist = results[:latent_event_hist]
    init_array = zeros(Float64, 0, 5)
    for (i, s) in enumerate(latent_event_hist)
        aux_array = zeros(Float64, length(s)+1, 5)
        for (j, z) in enumerate(s)
            aux_array[j, 1] = z.assignment_id
            aux_array[j, 2] = z.timestamp
            aux_array[j, 3] = z.seq_type
            aux_array[j, 4] = z.seq_warp
            aux_array[j, 5] = z.amplitude
        end
        aux_array[length(s)+1, :] = -1*ones(5)
        init_array = vcat(init_array, aux_array)
    end
    latent_event_hist_array = init_array
    LATENT_EVENT_SHAPE = size(init_array)

    latent_event_hist_frame = DataFrame(vcat(anneal_latent_event_hist_array, latent_event_hist_array),
        Symbol.(["assignment_id", "timestamp", "seq_type", "seq_warp", "amplitude"]))

    # --------- global_hist ----------
    global_vars_list = results[:globals_hist]
    N = size(global_vars_list[1].neuron_response_log_proportions)[1]
    R = size(global_vars_list[1].neuron_response_log_proportions)[2]
    init_array = zeros(Float64, 0, 3*R)
    seq_type_log_proportions = zeros(Float64, length(global_vars_list), R)
    bkgd_log_proportions = zeros(Float64, length(global_vars_list), N)
    for (i, s) in enumerate(global_vars_list)
        neuron_response = hcat(s.neuron_response_log_proportions, s.neuron_response_offsets,
        s.neuron_response_widths)
        init_array = vcat(init_array, neuron_response)
        seq_type_log_proportions[i, :] = s.seq_type_log_proportions
        bkgd_log_proportions[i, :] = s.bkgd_log_proportions
    end
    neuron_response_array = init_array
    seq_type_log_proportions_array = seq_type_log_proportions
    bkgd_log_proportions_array = bkgd_log_proportions

    # --------- anneal_global_hist ----------
        global_vars_list = results[:anneal_globals_hist]
    N = size(global_vars_list[1].neuron_response_log_proportions)[1]
    R = size(global_vars_list[1].neuron_response_log_proportions)[2]
    init_array = zeros(Float64, 0, 3*R)
    seq_type_log_proportions = zeros(Float64, length(global_vars_list), R)
    bkgd_log_proportions = zeros(Float64, length(global_vars_list), N)
    for (i, s) in enumerate(global_vars_list)
        neuron_response = hcat(s.neuron_response_log_proportions, s.neuron_response_offsets,
        s.neuron_response_widths)
        init_array = vcat(init_array, neuron_response)
        seq_type_log_proportions[i, :] = s.seq_type_log_proportions
        bkgd_log_proportions[i, :] = s.bkgd_log_proportions
    end
    anneal_neuron_response_array = init_array
    anneal_seq_type_log_proportions_array = seq_type_log_proportions
    anneal_bkgd_log_proportions_array = bkgd_log_proportions

    neuron_response_frame = DataFrame(vcat(anneal_neuron_response_array, neuron_response_array))
    seq_type_log_proportions_frame = DataFrame(vcat(anneal_seq_type_log_proportions_array, seq_type_log_proportions_array))
    bkgd_log_proportions_frame = DataFrame(vcat(anneal_bkgd_log_proportions_array, bkgd_log_proportions_array))

    # --------- assignment_hist ----------
    assigment_hist_frame = DataFrame(hcat(results[:anneal_assignment_hist], results[:assignment_hist]))

    aux_array = zeros(Float64, length(results[:log_p_hist]), 1)
    aux_array[:, 1] = results[:log_p_hist]
    log_p_hist_array = aux_array

    aux_array = zeros(Float64, length(results[:anneal_log_p_hist]), 1)
    aux_array[:, 1] = results[:anneal_log_p_hist]
    anneal_log_p_hist_array = aux_array
    log_p_hist_frame = DataFrame(vcat(anneal_log_p_hist_array, log_p_hist_array))
    # return log_p_hist_frame

    aux_array = zeros(Float64, length(results[:initial_assignments]), 1)
    aux_array[:, 1] = results[:initial_assignments]
    initial_assignments_frame = DataFrame(aux_array)

    results_directory = run_id*"_results/"
    if !isdir(results_directory)
        mkdir(results_directory)
    end

    CSV.write(results_directory*"latent_event_hist.csv", latent_event_hist_frame, writeheader=true)
    CSV.write(results_directory*"neuron_response.csv", neuron_response_frame, writeheader=true)
    CSV.write(results_directory*"seq_type_log_proportions.csv", seq_type_log_proportions_frame, writeheader=true)
    CSV.write(results_directory*"bkgd_log_proportions_array.csv", bkgd_log_proportions_frame, writeheader=true)
    CSV.write(results_directory*"log_p_hist.csv", log_p_hist_frame, writeheader=true)
    CSV.write(results_directory*"initial_assignments.csv", latent_event_hist_frame, writeheader=true)

    open(results_directory*"data_readme.txt", "w") do io
        write(io, "\n")
        write(io, "Latent_events: [N_spikes, iter] \n")
        write(io, "anneal_latent_event: "*string(ANNEAL_LATENT_EVENT_SHAPE)*"\n")
        write(io, "after_anneal_latent_event: "*string(LATENT_EVENT_SHAPE)*"\n")
        write(io, "full_latent_event: "*string(size(latent_event_hist_frame))*"\n")
        write(io, "Iterations separated by a row of -1 \n")
        write(io, "\n")
        write(io, "Neuron response: [log_proportions, offsets, widths] (neurons*iters) × 3R \n")
        write(io, "anneal_neuron_response: "*string(size(anneal_neuron_response_array))*"\n")
        write(io, "after_neuron_response: "*string(size(neuron_response_array))*"\n")
        write(io, "full_neuron_response: "*string(size(neuron_response_frame))*"\n")
        write(io, "\n")
        write(io, "assigment_hist: n_spikes x (100 iters anneal + 200 iters) \n")
        write(io, "assigment_hist: "*string(size(assigment_hist_frame))*"\n")
        write(io, "log_p_hist: "*string(size(log_p_hist_frame))*", (iters x 1)\n")
        write(io, "bkgd_log_proportions: "*string(size(bkgd_log_proportions_frame))*"\n")
        write(io, "seq_type_log_proportions: "*string(size(seq_type_log_proportions_frame))*"\n")
        write(io, "(anneal iters) + (without anneal iters) \n")
        write(io, "\n")
        write(io, "initial_assignments, "*string(size(initial_assignments_frame))*"\n")
        write(io, "For more details on globals and events, see \n https://github.com/rodrigcd/PPSeq.jl/blob/84d8c1da3f7fe55f93555ae94242d1309c7d44fb/src/model/structs.jl")
    end

    return latent_event_hist_frame

end

# Load spikes.
spikes = seq.Spike[]
file_name = Data_directory
for (n, t) in eachrow(readdlm(file_name, '\t', Float64, '\n'))
    push!(spikes, seq.Spike(Int(n), t))
end

# Setup the config details
config = Dict(

    # Model hyperparameters
    :num_sequence_types =>  5,
    :seq_type_conc_param => 1.0,
    :seq_event_rate => 1.0,

    :mean_event_amplitude => 100.0,
    :var_event_amplitude => 1000.0,

    :neuron_response_conc_param => 0.1,
    :neuron_offset_pseudo_obs => 1.0,
    :neuron_width_pseudo_obs => 1.0,
    :neuron_width_prior => 0.5,

    :num_warp_values => 1,
    :max_warp => 1.0,
    :warp_variance => 1.0,

    :mean_bkgd_spike_rate => 30.0,
    :var_bkgd_spike_rate => 30.0,
    :bkgd_spikes_conc_param => 0.3,
    :max_sequence_length => Inf,

    # MCMC Sampling parameters.
    :num_threads => 10,
    :num_anneals => 10,
    :samples_per_anneal => 100,
    :max_temperature => 40.0,
    :save_every_during_anneal => 10,
    :samples_after_anneal => 2000,
    :save_every_after_anneal => 10,
    :split_merge_moves_during_anneal => 10,
    :split_merge_moves_after_anneal => 10,
    :split_merge_window => 1.0,

);

# Then train the PPSeq Model
# Initialize all spikes to background process.
init_assignments = fill(-1, length(spikes))

# Construct model struct (PPSeq instance).
model = seq.construct_model(config, max_time, num_neurons)

# Run Gibbs sampling with an initial annealing period.
results = seq.easy_sample!(model, spikes, init_assignments, config);

# Save the results
save_results(results, "../Simple"*Dates.format(Dates.now(),"HH:MM"))


text:
"""
Unit tests for XYZ package.

------------------------------------------------------------------------------
COPYRIGHT/LICENSE. This file is part of the XYZ package. It is subject to
the license terms in the LICENSE file found in the top-level directory of
this distribution. No part of the XYZ package, including this file, may be
copied, modified, propagated, or distributed except according to the terms
contained in the LICENSE file.
------------------------------------------------------------------------------
"""
# --- Imports

# External packages
using Documenter
using Test
using TestTools: jltest

# Boid.jl
using Boid

# --- Run tests

jltest.run_tests(@__DIR__)


text:
@testset " Hydro Tests" begin
    PSI.activepower_variables(ps_model, generators_hg, 1:24)
    PSI.activepower_constraints(ps_model, generators_hg, PSI.HydroDispatchRunOfRiver, PM.DCPlosslessForm, 1:24)
    PSI.activepower_constraints(ps_model, generators_hg, PSI.HydroDispatchRunOfRiver, PM.StandardACPForm, 1:24)
    PSI.reactivepower_variables(ps_model, generators_hg, 1:24)
    PSI.reactivepower_constraints(ps_model, generators_hg, PSI.HydroDispatchRunOfRiver, PM.StandardACPForm, 1:24)
end

@testset " Hydro Tests" begin
    ps_model = PSI._canonical_model_init(length(sys5b.buses), nothing, PM.AbstractPowerFormulation, sys5b.time_periods)
    PSI.activepower_variables(ps_model, generators_hg, 1:24)
    PSI.commitment_variables(ps_model, generators_hg, 1:24);
    PSI.activepower_constraints(ps_model, generators_hg, PSI.HydroCommitmentRunOfRiver, PM.DCPlosslessForm, 1:24)
    PSI.activepower_constraints(ps_model, generators_hg, PSI.HydroCommitmentRunOfRiver, PM.StandardACPForm, 1:24)
    PSI.reactivepower_variables(ps_model, generators_hg, 1:24)
    PSI.reactivepower_constraints(ps_model, generators_hg, PSI.HydroCommitmentRunOfRiver, PM.StandardACPForm, 1:24)
end


text:
struct BinfileMeltRate <: AbstractMeltRate
    input_filename::String    #specify melt filename
end

function BinfileMeltRate(;
                        input_filename = nothing)

    #input file exceptions
    ~(input_filename === nothing) || throw(ArgumentError("You must pass an input filename"))
    isfile(input_filename) || throw(ArgumentError("Did not find the specified binary file"))

    return BinfileMeltRate(input_filename)
end

function update_melt_rate!(melt_model::BinfileMeltRate, fields, grid)
    @unpack input_filename = melt_model
    @unpack basal_melt = fields.gh

    file_size = stat(input_filename).size #bytes in input file must match matrix returned
    melt_rate = zeros(grid.nx,grid.ny)
    (file_size == sizeof(melt_rate)) || throw(DimensionMismatch("Size of input file incompatible with specified nx, ny"))
    try
        read!(input_filename, melt_rate)
    catch
        Error("Input file read error")
    end

   # melt_rate .= ntoh.(melt_rate)
    basal_melt[:] = melt_rate[:]
    return nothing
end


text:
"""
    load_OPFproblems(input_type, instance_path::String)

Read provided files in `instance_path` function of `input_type` and create a structure ds of type DataSource for each scenario. Also create empty structures gs and mp for each scenario
For each scenario, treat data in ds to complete gs : from bus and links in ds, generate node_linksin and node_linksout
For each scenario, complete mp structure
Return a dictonary Scenario => Scenario structure

# Arguments
- `input_type::T` : type of entry format for example MatpowerInput, GOCInput
- `instance_path::String` : path to the instance to treat, provided the instance path must be associated to an instance of format input_type

# Output
- `OPFpbs`

# Examples
```
julia > OPFpbs = load_OPFproblems(MatpowerInput, "instances\\matpower\\WB2.m")
julia > OPFpbs
SortedDict{String,Scenario} with 1 entry:
  "BaseCase" => Scenario
```
"""

function load_OPFproblems(input_type, instance_path::String)
  # raw = "powersystem.raw"
  # gen = "generator.csv"
  # con = "contingency.csv"
  #
  # rawfile = joinpath(instance_path,raw)
  # genfile = joinpath(instance_path, gen)
  # confile = joinpath(instance_path, con)
  #
  # OPFpbs = read_GOCfiles(rawfile, genfile, confile)

  OPFpbs = read_input(input_type, instance_path)
  # baseMVA, bus, link, node_elems, link_elems, elem_to_typekey = read_input(input_type, filenames)
  # ds = DataSource(baseMVA, bus, link)

  for scenario in keys(OPFpbs)
    node_linksin, node_linksout = SortedDict{String, SortedSet{Link}}(), SortedDict{String, SortedSet{Link}}()

    for (linkname, _) in OPFpbs[scenario].ds.link
      if !haskey(node_linksin, linkname.dest)
        node_linksin[linkname.dest] = SortedSet{Link}()
      end
      union!(node_linksin[linkname.dest], [linkname])

      if !haskey(node_linksout, linkname.orig)
        node_linksout[linkname.orig] = SortedSet{Link}()
      end
      union!(node_linksout[linkname.orig], [linkname])
    end

    OPFpbs[scenario].gs.node_linksin = node_linksin
    OPFpbs[scenario].gs.node_linksout = node_linksout

    node_vars=SortedDict{String, SortedDict{String, Variable}}(busname => SortedDict{String, Variable}() for busname in keys(OPFpbs[scenario].ds.bus))
    link_vars=SortedDict{Link, SortedDict{String, Variable}}(link => SortedDict{String, Variable}() for link in keys(OPFpbs[scenario].ds.link))

    node_formulations = SortedDict{String, SortedDict{String, Symbol}}()
    link_formulations = SortedDict{Link, SortedDict{String, Symbol}}()
    for bus in keys(OPFpbs[scenario].ds.bus)
      node_formulations[bus] = SortedDict{String, Symbol}(elem => :NbMinVar for elem in keys(OPFpbs[scenario].ds.bus[bus]))
    end
    for link in keys(OPFpbs[scenario].ds.link)
      link_formulations[link] = SortedDict{String, Symbol}(elem => :NbMinVar for elem in keys(OPFpbs[scenario].ds.link[link]))
    end

    OPFpbs[scenario].mp = MathematicalProgramming(node_formulations, link_formulations, node_vars, link_vars)
  end

  return OPFpbs
end

function load_OPFproblems(rawfile::String, genfile::String, confile::String)
  OPFpbs = read_GOCfiles(rawfile, genfile, confile)

  for scenario in keys(OPFpbs)
    node_linksin, node_linksout = SortedDict{String, SortedSet{Link}}(), SortedDict{String, SortedSet{Link}}()

    for (linkname, _) in OPFpbs[scenario].ds.link
      if !haskey(node_linksin, linkname.dest)
        node_linksin[linkname.dest] = SortedSet{Link}()
      end
      union!(node_linksin[linkname.dest], [linkname])

      if !haskey(node_linksout, linkname.orig)
        node_linksout[linkname.orig] = SortedSet{Link}()
      end
      union!(node_linksout[linkname.orig], [linkname])
    end

    OPFpbs[scenario].gs.node_linksin = node_linksin
    OPFpbs[scenario].gs.node_linksout = node_linksout

    node_vars=SortedDict{String, SortedDict{String, Variable}}(busname => SortedDict{String, Variable}() for busname in keys(OPFpbs[scenario].ds.bus))
    link_vars=SortedDict{Link, SortedDict{String, Variable}}(link => SortedDict{String, Variable}() for link in keys(OPFpbs[scenario].ds.link))

    node_formulations = SortedDict{String, SortedDict{String, Symbol}}()
    link_formulations = SortedDict{Link, SortedDict{String, Symbol}}()
    for bus in keys(OPFpbs[scenario].ds.bus)
      node_formulations[bus] = SortedDict{String, Symbol}(elem => :NbMinVar for elem in keys(OPFpbs[scenario].ds.bus[bus]))
    end
    for link in keys(OPFpbs[scenario].ds.link)
      link_formulations[link] = SortedDict{String, Symbol}(elem => :NbMinVar for elem in keys(OPFpbs[scenario].ds.link[link]))
    end

    OPFpbs[scenario].mp = MathematicalProgramming(node_formulations, link_formulations, node_vars, link_vars)
  end

  return OPFpbs
end


text:
#=
There is an N by M matrix of zeroes. Given N and M, write a function to count the number of ways of starting at the top-left corner and getting to the bottom-right corner. You can only move right or down.

For example, given a 2 by 2 matrix, you should return 2, since there are two ways to get to the bottom-right:

    Right, then down
    Down, then right

Given a 5 by 5 matrix, there are 70 ways to get to the bottom-right.
=#

function topleft_to_bottomright(matrix)
    if size(matrix,1) == 1 && size(matrix,2) == 1
        return 1
    elseif size(matrix,1) > 1 && size(matrix,2) > 1
        return topleft_to_bottomright(matrix[2:end,:]) + topleft_to_bottomright(matrix[:,2:end])
    else
        return 1
    end
end

max_stars_repo_name: DominiqueCaron/daily-coding-problem

max_issues_repo_name: DominiqueCaron/daily-coding-problem

max_forks_repo_name: DominiqueCaron/daily-coding-problem


text:
#!/usr/bin/env julia

# Create unstructured grid VTK file.

using WriteVTK

const FloatType = Float32
const vtk_filename_noext = "unstructured"

# 3D mesh
function mesh_data(::Val{3})
    # This is basically a structured grid, but defined as an unstructured one.
    # Based on the structured.jl example.
    Ni, Nj, Nk = 40, 50, 20
    indices = LinearIndices((1:Ni, 1:Nj, 1:Nk))
    Npts = length(indices)

    # Create points and point data.
    pts_ijk = Array{FloatType}(undef, 3, Ni, Nj, Nk)
    pdata_ijk = Array{FloatType}(undef, Ni, Nj, Nk)

    for k = 1:Nk, j = 1:Nj, i = 1:Ni
        r = 1 + (i - 1)/(Ni - 1)
        th = 3*pi/4 * (j - 1)/(Nj - 1)
        pts_ijk[1, i, j, k] = r * cos(th)
        pts_ijk[2, i, j, k] = r * sin(th)
        pts_ijk[3, i, j, k] = (k - 1)/(Nk - 1)
        pdata_ijk[i, j, k] = i*i + k*sqrt(j)
    end

    # Create cells (all hexahedrons in this case) and cell data.
    celltype = VTKCellTypes.VTK_HEXAHEDRON
    cells = MeshCell[]
    cdata = FloatType[]

    for k = 2:Nk, j = 2:Nj, i = 2:Ni
        # Define connectivity of cell.
        inds = Array{Int32}(undef, 8)
        inds[1] = indices[i-1, j-1, k-1]
        inds[2] = indices[i  , j-1, k-1]
        inds[3] = indices[i  , j  , k-1]
        inds[4] = indices[i-1, j  , k-1]
        inds[5] = indices[i-1, j-1, k  ]
        inds[6] = indices[i  , j-1, k  ]
        inds[7] = indices[i  , j  , k  ]
        inds[8] = indices[i-1, j  , k  ]

        # Define cell.
        c = MeshCell(celltype, inds)

        push!(cells, c)
        push!(cdata, i*j*k)
    end

    pts = reshape(pts_ijk, 3, Npts)
    pdata = reshape(pdata_ijk, Npts)

    return pts, cells, pdata, cdata
end

# 2D mesh
function mesh_data(::Val{2})
    Ni, Nj = 40, 50
    indices = LinearIndices((1:Ni, 1:Nj))
    Npts = length(indices)

    # Create points and point data.
    pts_ijk = Array{FloatType}(undef, 2, Ni, Nj)
    pdata_ijk = Array{FloatType}(undef, Ni, Nj)

    for j = 1:Nj, i = 1:Ni
        r = 1 + (i - 1)/(Ni - 1)
        th = 3*pi/4 * (j - 1)/(Nj - 1)
        pts_ijk[1, i, j] = r * cos(th)
        pts_ijk[2, i, j] = r * sin(th)
        pdata_ijk[i, j] = i*i + sqrt(j)
    end

    # Create cells (all quads in this case) and cell data.
    celltype = VTKCellTypes.VTK_QUAD
    cells = MeshCell[]
    cdata = FloatType[]

    for j = 2:Nj, i = 2:Ni
        # Define connectivity of cell.
        inds = Array{Int32}(undef, 4)
        inds[1] = indices[i-1, j-1]
        inds[2] = indices[i  , j-1]
        inds[3] = indices[i  , j  ]
        inds[4] = indices[i-1, j  ]

        # Define cell.
        c = MeshCell(celltype, inds)

        push!(cells, c)
        push!(cdata, i*j)
    end

    pts = reshape(pts_ijk, 2, Npts)
    pdata = reshape(pdata_ijk, Npts)

    return pts, cells, pdata, cdata
end

# 1D mesh
function mesh_data(::Val{1})
    Ni = 40
    Npts = Ni

    # Create points and point data.
    pts_ijk = Array{FloatType}(undef, 1, Ni)
    pdata_ijk = Array{FloatType}(undef, Ni)

    for i = 1:Ni
        pts_ijk[1, i] = (i-1)^2 / (Ni-1)^2
        pdata_ijk[i] = i*i
    end

    # Create cells (all lines in this case) and cell data.
    celltype = VTKCellTypes.VTK_LINE
    cells = MeshCell[]
    cdata = FloatType[]

    for i = 2:Ni
        # Define connectivity of cell.
        inds = Array{Int32}(undef, 2)
        inds[1] = i-1
        inds[2] = i

        # Define cell.
        c = MeshCell(celltype, inds)

        push!(cells, c)
        push!(cdata, i)
    end

    pts = reshape(pts_ijk, 1, Npts)
    pdata = reshape(pdata_ijk, Npts)

    return pts, cells, pdata, cdata
end

function main()
    outfiles = String[]
    for dim in 1:3
        pts, cells, pdata, cdata = mesh_data(Val{dim}())

        @time begin
            # Initialise new vtu file (unstructured grid).
            vtk = vtk_grid(vtk_filename_noext*"_$(dim)D", pts, cells,
                           compress=3)

            # This is also accepted in 3D:
            # vtk = vtk_grid(vtk_filename_noext, pts[1, :], pts[2, :], pts[3, :],
            #                cells)

            # Add some point and cell data.
            vtk_point_data(vtk, pdata, "my_point_data")
            vtk_cell_data(vtk, cdata, "my_cell_data")

            # Save and close vtk file.
            append!(outfiles, vtk_save(vtk))
        end
    end

    println("Saved:  ", join(outfiles, "  "))
    return outfiles::Vector{String}
end

main()


text:
#=

# Comparison with DifferentialEquations.jl
([source code](SOURCE_URL))

Solve the nonlinear diffusion equation

```math
\partial_t u -\Delta u^m = 0
```
in $\Omega=(-1,1)$ with homogeneous Neumann boundary conditons using the implicit Euler method.

This equation is also called  "porous medium equation".
The Barenblatt solution is an exact solution of this problem which for m>1 has a finite support.
We initialize this problem with the exact solution for $t=t_0=0.001$.

(see Barenblatt, G. I. "On nonsteady motions of gas and fluid in porous medium." Appl. Math. and Mech.(PMM) 16.1 (1952): 67-78.)

Here, we compare the implicit Euler approach in VoronoiFVM with the ODE solvers in [DifferentialEquations.jl](https://diffeq.sciml.ai/stable/) and demonstrate the possibility to use VoronoiFVM to define differential
operators compatible with its [ODEFunction](https://diffeq.sciml.ai/stable/features/performance_overloads/#ODEFunction)
interface.

Note that this example requires  PyPlot and DifferentialEquations to be installed.
=#

module TestDiffEq

using VoronoiFVM
using DifferentialEquations
using LinearAlgebra
using Printf
using PyPlot

import Base:push!

function barenblatt(x,t,m)
    tx=t^(-1.0/(m+1.0))
    xx=x*tx
    xx=xx*xx
    xx=1- xx*(m-1)/(2.0*m*(m+1));
    if xx<0.0
        xx=0.0
    end
    return tx*xx^(1.0/(m-1.0))
end

function create_porous_medium_problem(n,m,unknown_storage)
    h=1.0/convert(Float64,n/2)
    X=collect(-1:h:1)
    grid=VoronoiFVM.Grid(X)

    function flux!(f,u0,edge)
        u=unknowns(edge,u0)
        f[1]=u[1,1]^m-u[1,2]^m
    end

    storage!(f,u,node)= f[1]=u[1]

    physics=VoronoiFVM.Physics(flux=flux!,storage=storage!)

    sys=VoronoiFVM.System(grid,physics,unknown_storage=unknown_storage)
    enable_species!(sys,1,[1])
    sys,X
end

function run_vfvm(;n=20,m=2,t0=0.001, tend=0.01,tstep=1.0e-6,unknown_storage=:dense)

    sys,X=create_porous_medium_problem(n,m,unknown_storage)

    inival=unknowns(sys)
    inival[1,:].=map(x->barenblatt(x,t0,m),X)

    solution=unknowns(sys)
    control=VoronoiFVM.NewtonControl()
    control.verbose=false
    control.Δt=tstep
    control.Δu_opt=0.05
    control.Δt_min=tstep

    times=collect(t0:t0:tend)
    times=[t0,tend]
    sol=VoronoiFVM.solve(inival,sys,times,control=control,store_all=true)
    err=norm(sol[1,:,end]-map(x->barenblatt(x,tend,m),X))
    sol,X,err
end

function run_diffeq(;n=20,m=2, t0=0.001,tend=0.01, unknown_storage=:dense,solver=nothing)
    sys,X=create_porous_medium_problem(n,m,unknown_storage)
    inival=unknowns(sys)
    inival[1,:].=map(x->barenblatt(x,t0,m),X)
    tspan = (t0,tend)
    sol=VoronoiFVM.solve(DifferentialEquations,inival,sys,tspan,solver=solver)
    err=norm(sol[1,:,end]-map(x->barenblatt(x,tend,m),X))
    sol, X,err
end

function main(;m=2,n=20, solver=nothing, unknown_storage=:dense)
    function plotsol(sol,X)
        f=sol[1,:,:]'
        contourf(X,sol.t,f,0:0.1:10,cmap=:summer)
        contour(X,sol.t,f,0:1:10,colors=:black)
    end

    clf()
    subplot(121)

    t=@elapsed begin
        sol,X,err=run_vfvm(m=m,n=n, unknown_storage=unknown_storage)
    end
    title(@sprintf("VoronoiFVM: %.0f ms e=%.2e",t*1000,err))
    plotsol(sol,X)

    subplot(122)
    t=@elapsed begin
        sol,X,err=run_diffeq(m=m,n=n,solver=solver, unknown_storage=unknown_storage)
    end
    plotsol(sol,X)
    title(@sprintf("DifferentialEq: %.0f ms, e=%.2e",t*1000,err))

    gcf().set_size_inches(8,4)
    gcf()
end
end


text:
@testset "utils tests" begin

end


text:
using Plots
using NPZ
using ProgressMeter

# load the simulated data
data = npzread("../lorenz_sim_data.npy")
n = 2000

# initialize 3d plot
plt = path3d(1, xlim=(-25,25), ylim=(-25,25), zlim=(0,50),
                xlab = "x", ylab = "y", zlab = "z",
                markercolor=RGBA(254/255,190/255,31/255,250/255),
                title = "Lorenz Attractor", marker = 1)

# plot simulated data
p = Progress(n, 0.01)
anim = @animate for i=1:n
    push!(plt, data[1,i], data[2,i], data[3,i])
    update!(p, i)
end every 3

gif(anim, "./simulated_lorenz_attractor.gif", fps = 15)


text:
# examples:
#   search_wdi("countries","name",r"united"i)
#   search_wdi("indicators","description",r"gross national"i)
function search_wdi(data::ASCIIString, entry::Symbol, regx::Regex)
    data_opts = ["countries","indicators"]
    if !(data in data_opts)
        error("unsupported data source:", data, ". supported are:\n",data_opts)
    end
    if data == "countries"
        return search_countries(entry,regx)
    end
    if data == "indicators"
        return search_indicators(entry,regx)
    end
end

function search_countries(entry::Symbol,regx::Regex)
    # test whether chosen column name is valid
    entries = [:iso3c,
               :iso2c,
               :name,
               :region,
               :regionId,
               :capital,
               :longitude,
               :latitude,
               :income,
               :incomeId,
               :lending,
               :lendingId
               ]

    if !(entry in entries)
        error("unsupported country entry: \"",entry,"\". supported are:\n",entries)
    end
    country_match(entry,regx)
end

function search_indicators(entry::Symbol,regx::Regex)
    # test whether chosen column name is valid
    entries = [:indicator,
               :name,
               :description,
               :source_database,
               :source_databaseId,
               :source_organization
               ]

    if !(entry in entries)
        error("unsupported indicator entry: \"",entry,"\". supported are\n",entries)
    end
    indicator_match(entry,regx)
end

function country_match(entry::Symbol,regex::Regex)
    df = loadWBMeta("countries")
    df_match(df,entry,regex)
end

function indicator_match(entry::Symbol,regex::Regex)
    df = loadWBMeta("indicators")
    df_match(df,entry,regex)
end

function regex_match(da::DataArray{UTF8String,1},
                     regex::Regex)
    toCheck = !isna(da)
    nObs = size(da, 1)
    res = DataArray(Bool, nObs)
    res[isna(da)] = false
    checkRes = map(x -> ismatch(regex,x), da[toCheck])
    res[toCheck] = checkRes
    return res
end

function regex_match(da::DataArray{ASCIIString,1},
                     regex::Regex)
    toCheck = !isna(da)
    nObs = size(da, 1)
    res = DataArray(Bool, nObs)
    res[isna(da)] = false
    checkRes = map(x -> ismatch(regex,x), da[toCheck])
    res[toCheck] = checkRes
    return res
end

function df_match(df::AbstractDataFrame,
                  entry::Symbol,
                  regex::Regex)
    return df[regex_match(df[entry],regex),:] # use symbol to get
                                        # column
end


text:
module WaveGPU

# using CairoMakie
# using GLMakie
using ArrayFire: AFArray, allowslow, swap!
using DDF
using DiffEqCallbacks
using DifferentialForms
using LinearAlgebra
using Logging: global_logger
using OrdinaryDiffEq
using Printf
using ProgressLogging
using SparseArrays
using StaticArrays
using TerminalLoggers: TerminalLogger

################################################################################

walltime() = time_ns() / 1.0e+9

function Base.copyto!(y::AFArray, x::AFArray)
    @assert length(y) == length(x)
    z = copy(x)
    swap!(y, z)
    finalize(z)
    return y
end
function LinearAlgebra.mul!(y::AFArray, A::AFArray, x::AFArray)
    z = A * x
    swap!(y, z)
    finalize(z)
    return y
end

global_logger(TerminalLogger())

################################################################################

function wave(::Val{D}, levels::Int; gpu::Bool=false) where {D}
    D::Int
    S = Float32 #Float64
    T = Float32 #Float64

    println("Create manifold...")

    n = 2^levels
    nelts = Tuple(n for d in 1:D)
    mfd = large_hypercube_manifold(Val(D), S; nelts=nelts, optimize_mesh=false)

    ############################################################################

    println("Sample initial conditions...")

    # TODO: use first-order formulation

    k = SVector{D,S}(1 for d in 1:D)
    ω = norm(k)
    function u⁼(t, x)
        local u = cospi(ω * t) * prod(sinpi(k[d] * x[d]) for d in 1:D)
        return Form{D,0}((u,))
    end
    function ∂ₜu⁼(t, x)
        local ∂ₜu = -S(π) *
                    ω *
                    sinpi(ω * t) *
                    prod(sinpi(k[d] * x[d]) for d in 1:D)
        return Form{D,0}((∂ₜu,))
    end

    t₀ = S(0)
    u₀ = sample(Fun{D,Pr,0,D,S,T}, x -> u⁼(t₀, x), mfd)
    ∂ₜu₀ = sample(Fun{D,Pr,0,D,S,T}, x -> ∂ₜu⁼(t₀, x), mfd)

    ############################################################################

    println("Set up equations of motion...")

    # ∂ₜₜu = Δu
    # u̇ := ∂ₜu

    # ∂ₜu = u̇
    # ∂ₜu̇ = Δu

    combine(u, u̇) = [u.values.vec; u̇.values.vec]
    function split(U)
        local n = length(u₀)
        @assert length(U) == 2n
        local uv = U[1:n]
        local u̇v = U[(n + 1):(2n)]
        return (Fun{D,Pr,0,D,S,T}(mfd, IDVector{0}(uv)),
                Fun{D,Pr,0,D,S,T}(mfd, IDVector{0}(u̇v)))
    end

    U₀ = combine(u₀, ∂ₜu₀)

    N = zero(Op{D,Pr,0,Pr,0}, mfd)
    E = one(Op{D,Pr,0,Pr,0}, mfd)
    B = isboundary(Val(Pr), Val(0), mfd)
    L = laplace(Val(Pr), Val(0), mfd)

    function f(U, p, t)
        local u, u̇ = split(U)
        local ∂ₜu = u̇
        local ∂ₜu̇ = L * u
        # local ∂ₜuₜ = sample(Fun{D,Pr,0,D,S,T}, x -> ∂ₜu⁼(t, x), mfd)
        # local ∂ₜₜuₜ = sample(Fun{D,Pr,0,D,S,T}, x -> ∂ₜₜu⁼(t, x), mfd)
        ∂ₜu = (E - B) * ∂ₜu # + B * ∂ₜuₜ
        ∂ₜu̇ = (E - B) * ∂ₜu̇ # + B * ∂ₜₜuₜ
        return combine(∂ₜu, ∂ₜu̇)
    end

    F = [N.values sparse(E.values); L.values N.values]
    F = [(E - B).values N.values; N.values (E - B).values] * F
    f!(∂ₜU, U, p, t) = mul!(∂ₜU, F, U)

    ############################################################################

    println("Set up GPU...")

    if gpu
        allowslow(AFArray, false)
        # copy data to GPU
        U₀ = AFArray(U₀)
        F = AFArray(F)
    end

    ############################################################################

    println("Solve...")

    Δx = minimum(get_volumes(mfd, 1))
    # @show Δx

    t₁ = S(1)

    # TODO: Experiment with linear ODE, Hamiltonian systems

    # TODO: Don't use `saveat`, use a `SavingCallback` instead

    # prob = ODEProblem(f, U₀, (t₀, t₁); saveat=(t₁ - t₀) / 8, save_end=false)
    prob = ODEProblem(f!, U₀, (t₀, t₁); saveat=(t₁ - t₀) / 8, save_end=false)

    dt = Δx / 2

    sol = nothing
    nextoutput = walltime()
    @withprogress name = "ODE" begin
        function progress(U, t, integrator)
            if t == t₁ || walltime() ≥ nextoutput
                iter = integrator.iter
                time = integrator.t
                nsol = norm(integrator.u, Inf)
                @printf "    iter=%d   time=%#.6g   ‖sol‖∞=%#.6g\n" iter time nsol
                @logprogress (t - t₀) / (t₁ - t₀)
                nextoutput = walltime() + 1
            end
        end

        sol = solve(prob, Tsit5(); adaptive=false, dt=dt, alias_u0=true,
                    internalnorm=(u, t) -> S(0),
                    # unstable_check=(dt,u,p,t)->false,
                    callback=FunctionCallingCallback(progress))
    end

    ############################################################################

    println("Analyse solution...")

    # @progress name = "Plotting" for i in axes(sol.t, 1)
    #     # t = (n - i) / S(n) * t₀ + i / S(n) * t₁
    #     # U = sol(t)
    #     t = sol.t[i]
    #     U = sol[i]
    #     u, ∂ₜu = split(U)
    #     uₜ = sample(Fun{D,Pr,0,D,S,T}, x -> u⁼(t, x), mfd)
    #     e = u - uₜ
    #     outdir = "wave$(D)d-$levels"
    #     mkpath(outdir)
    #     plot_function(u, "$outdir/wave$(D)d-u.$(@sprintf "%04d" i - 1).png")
    #     plot_function(e, "$outdir/wave$(D)d-e.$(@sprintf "%04d" i - 1).png")
    # end

    i = last(axes(sol.t, 1))
    t = sol.t[i]
    U = sol[i]
    if gpu
        # copy data from GPU
        U = Array(U)
    end
    u, ∂ₜu = split(U)
    uₜ = sample(Fun{D,Pr,0,D,S,T}, x -> u⁼(t, x), mfd)
    e = u - uₜ
    ne = sqrt(norm(e)^2 / length(e))
    println("    ‖u-u⁼‖₂=$ne")

    ############################################################################

    println("Done.")

    return nothing
end

# D = 2:
#     levels    error ‖u-u⁼‖₂
#     1         0.02805463824876254
#     2         0.0020292411087007196
#     3         0.00013532326532406398
#     4         8.727862238311897e-6
#     5         5.540717968813277e-7
#     6         3.49003314447833e-8
#     7         2.189781209413498e-9
#     8         1.3712800517004272e-10
#     9         8.580900735370499e-12
#     10        5.369060952547811e-13
#     11        3.1564728009539904e-14
#     12        6.75775317464105e-15
#     13        1.5417520244384893e-14

end


text:
using MathProgBase

function knapsack()
    f = Float64[3, 2, 3, 1, 2]
    A = [7. 2. 9. 3. 1.]
    capacity = 10.
    solution = mixintprog(-f, A, '<', capacity, 'I', 0, 1, MIPSolver(:Cbc,LogLevel=0))
    println("Solution status: ", solution.status)
    println("Optimal value: ", solution.objval) # should be -7
    print("Solution vector: ")
    show(solution.sol)
    println()

end

knapsack()


text:
module Reindex

export ShiftReindex,OddReindex,Reflect

struct ShiftReindex{T,N,O}
  array :: Vector{T}
end

function ShiftReindex(v::Vector{T},o::Int) where T
  N = length(v)+o
  ShiftReindex{T,N,o}(v)
end
(vec::ShiftReindex{T,N,O})(i::Int) where {T,N,O} =
                    (i > N || i < 1+O) ? T(0) : vec.array[i-O]
(vec::ShiftReindex{T,N,O})(ir::AbstractArray{M}) where {T,N,O,M} =
                    (vec::ShiftReindex{T,N,O}).(collect(ir))

struct OddReindex{T,N}
   array :: Vector{T}
end
function OddReindex(v::Vector{T}) where T
    N = length(v)-1
    OddReindex{T,N}(v)
end
(vec::OddReindex{T,N})(i::Int) where {T,N} =
        abs(i) > N ? T(0) : i < 0 ? conj(vec.array[1-i]) : vec.array[1+i]
(vec::OddReindex{T,N})(ir::AbstractArray{M}) where {T,N,M} =
                    (vec::OddReindex{T,N}).(collect(ir))
ReindexedArray = Union{ShiftReindex,OddReindex}

struct Reflect{T}
   array :: T
end
(vec::Reflect{<:ReindexedArray})(i) = vec.array(-i)
(vec::Reflect{<:AbstractArray})(i)  = vec.array[-i]

end