text:
function backtracking_line_search!(d::Union(DifferentiableFunction,
                                            TwiceDifferentiableFunction),
                                   x::Vector,
                                   p::Vector,
                                   new_x::Vector,
                                   new_gradient::Vector;
                                   c1::Real = 1e-6,
                                   c2::Real = 0.9,
                                   rho::Real = 0.9,
                                   iterations::Integer = 1_000)

    # Keep track of the number of iterations
    iteration = 0

    # Keep track of calls to f and g
    f_calls, g_calls = 0, 0

    # Calculate parameter space size
    n = length(x)

    # Evaluate the function and gradient at current position
    f_calls += 1
    g_calls += 1
    f_x = d.fg!(x, new_gradient)
    gxp = dot(new_gradient, p)

    # The default step-size is always 1.0
    alpha = 1.0

    # Propose a new x after moving length alpha in direction p
    for i in 1:n
        new_x[i] = x[i] + alpha * p[i]
    end

    # Expand step-size
    # NB: This is very expensive if the gradient costs much more than f(x)
    # In my tests, it also increases the number of steps, not decreases it
    # while (dot(new_gradient, p) < c2 * gxp) && (alpha < 65536.0)
    #     alpha *= 2.0
    #     for i in 1:n
    #         new_x[i] = x[i] + alpha * p[i]
    #     end
    #     g_calls += 1
    #     g!(new_x, new_gradient)
    # end

    # Keep coming closer to x until we find a point that is
    # as good as the gradient suggests we can achieve
    f_calls += 1
    f_new_x = d.f(new_x)
    while f_new_x > f_x + c1 * alpha * gxp
        iteration += 1
        if iteration > iterations
            error("Too many iterations in backtracking_line_search")
        end
        alpha *= rho
        for i in 1:n
            new_x[i] = x[i] + alpha * p[i]
        end
        f_calls += 1
        f_new_x = d.f(new_x)
    end

    return alpha, f_calls, g_calls
end


text:
"""
The `echo` integration test runs a fake WebSocket server that echoes all frames sent to it.

It runs in process, without any actual network sockets. It does this by faking the connection using
a channel instead of an actual socket. The vast majority of the code can still be executed, as is.
Only a small portion of code is specific to the actual network.
"""
using DandelionWebSockets: HTTPHandshake
using Random

"""
EchoServerHandler defines the behaviour of the fake in-process server that echoes anything sent to
it.
"""
struct EchoServerHandler <: WebSocketHandler
    connection::AbstractWSClient
end

on_text(e::EchoServerHandler, text::String) = send_text(e.connection, text)
on_binary(e::EchoServerHandler, data::AbstractVector{UInt8}) = send_binary(e.connection, data)
state_closed(e::EchoServerHandler) = nothing
state_closing(e::EchoServerHandler) = nothing
state_connecting(e::EchoServerHandler) = nothing
state_open(e::EchoServerHandler) = nothing

"""
EchoClientHandler defines the behaviour of the test, which is to send a bunch of messages to the
server and verify that we get an echo back for all of them.
"""
mutable struct EchoClientHandler
    client::WSClient
end

on_text(e::EchoClientHandler, text::String) = verifyecho(e, text)
on_binary(e::EchoClientHandler, data::AbstractVector{UInt8}) = verifyecho(e, data)
state_closed(e::EchoClientHandler) = nothing
state_closing(e::EchoClientHandler) = nothing
state_connecting(e::EchoClientHandler) = nothing
state_open(e::EchoClientHandler) = nothing

function verifyecho(e::EchoClientHandler, text::String)
    expected = "$(e.counter)"
    if expected == text
        if e.counter >= e.limit
            success(e)
        else
            e.counter += 1
            sendprobe(e, "$(e.counter)")
        end
    else
        println("ERROR: Expected '$expected', but got '$text'")
        failure(e)
    end
end

function sendprobe(e::EchoClientHandler, text::String)
    sendtext(e.client, text)
end

success(e::EchoClientHandler) = println("SUCCESS!")
failure(e::EchoClientHandler) = println("FAILURE!")

server = InProcessServer(EchoServerHandler)
# This replaces the `HTTPjlAdapter` type, which uses the `HTTP.jl` package to make an actual HTTP
# connection. This adapter instead asks the server for a channel, and uses that for returning a
# correct HTTP upgrade reply.
handshakeadapter = InProcessHandshakeAdapter(server)

client = WSClient(; handshake=HTTPHandshake(RandomDevice(), handshakeadapter))
echoclient = EchoClientHandler(client)

wsconnect(client, "ws://the/uri/does/not/matter/here", echoclient)

ext: jl


text:
abstract CorrectionMethod
immutable BonferroniMethod <: CorrectionMethod end
immutable BenjaminiHochbergMethod <: CorrectionMethod end
immutable BenjaminiYekutieliMethod <: CorrectionMethod end
immutable HochbergMethod <: CorrectionMethod end
immutable HolmMethod <: CorrectionMethod end

# some constants
Bonferroni = BonferroniMethod()
BenjaminiHochberg = BenjaminiHochbergMethod()
BenjaminiYekutieli = BenjaminiYekutieliMethod()
Hochberg = HochbergMethod()
Holm = HolmMethod()


text:
using CUDA
using LULESH
using MPI
using Printf

if length(ARGS) < 2
    printUsage()
    error("Wrong input arguments")
end

if !occursin("-u", ARGS[1]) && !occursin("-s", ARGS[1])
    printUsage()
    error("Wrong input arguments")
end

num_iters = 10 # AD_LENGTH of simulation usually set to "-1" but for debugging set to 10!
if length(ARGS) == 4
    num_iters = parse(Int, ARGS[4])
end

#   bool structured = ( strcmp(argv[1],"-s") == 0 );
structured = occursin("-s", ARGS[1])
structured = true
@show structured
# assume cube subdomain geometry for now (nx)
nx = parse(IndexT, ARGS[2])
# nx = 45
# TODO: change default nr to 11
nr = 1
balance = 1
cost = 1
devicetype = Vector
floattype = Float64
# To disable MPI, pass 'nothing' as the last argument
prob = LuleshProblem(num_iters, structured, nx, nr, balance, cost, Vector, Float64, nothing)

#   Int_t numRanks ;
#   Int_t myRank ;

# #if USE_MPI
#   Domain_member fieldData ;

#   MPI_Init(&argc, &argv) ;
#   MPI_Comm_size(MPI_COMM_WORLD, &numRanks) ;
#   MPI_Comm_rank(MPI_COMM_WORLD, &myRank) ;
# #else
#   numRanks = 1;
#   myRank = 0;
# #endif

#   cuda_init(myRank);

#   Domain *locDom ;

# Set up the mesh and decompose. Assumes regular cubes for now

# TODO: modify this constructor to account for new fields
# TODO: setup communication buffers

domain = Domain(prob)

# #if USE_MPI
#    // copy to the host for mpi transfer
#    locDom->h_nodalMass = locDom->nodalMass;

#    fieldData = &Domain::get_nodalMass;

#    // Initial domain boundary communication
#    CommRecv(*locDom, MSG_COMM_SBN, 1,
#             locDom->sizeX + 1, locDom->sizeY + 1, locDom->sizeZ + 1,
#             true, false) ;
#    CommSend(*locDom, MSG_COMM_SBN, 1, &fieldData,
#             locDom->sizeX + 1, locDom->sizeY + 1, locDom->sizeZ + 1,
#             true, false) ;
#    CommSBN(*locDom, 1, &fieldData) ;

#    // copy back to the device
#    locDom->nodalMass = locDom->h_nodalMass;

#    // End initialization
#    MPI_Barrier(MPI_COMM_WORLD);
# #endif

#   cudaDeviceSetCacheConfig(cudaFuncCachePreferL1);
#   cudaDeviceSetLimit(cudaLimitMallocHeapSize,1024*1024*1024);

@show domain.time_h
@show domain.stoptime
if getMyRank(prob.comm) == 0
    if (structured)
        @printf("Running until t=%f, Problem size=%dx%dx%d\n", domain.stoptime, nx, nx, nx)
    else
        @printf("Running until t=%f, Problem size=%d \n", domain.stoptime, domain.numElem)
        @warn "Unstructured setup not supported"
    end
end
# timestep to solution
its = 0
#   cudaProfilerStart();
start = getWtime(prob.comm)
while domain.time_h < domain.stoptime
    # this has been moved after computation of volume forces to hide launch latencies
    # TimeIncrement(locDom)
    lagrangeLeapFrog(domain)
    # checkErrors(domain, its, myRank)
    if getMyRank(prob.comm) == 0
        @printf("cycle = %d, time = %e, dt=%e\n", its+1, domain.time_h, domain.deltatime_h)
    end
    global its += 1
    # if its == num_iters
        break
    # end
end

#   // make sure GPU finished its work
#   cudaDeviceSynchronize();

# Use reduced max elapsed time
elapsed_time = getWtime(prob.comm) - start
elapsed_timeG = comm_max(elapsed_time, prob.comm)

#   cudaProfilerStop();

#   if (myRank == 0)
#     VerifyAndWriteFinalOutput(elapsed_timeG, *locDom, its, nx, numRanks, structured);

# #ifdef SAMI
#   DumpDomain(locDom) ;
# #endif
#   cudaDeviceReset();

MPI.Finalize()


text:
@inline isnearint(v, tol) = abs(v-round(v)) < tol

# TODO (cg 2021/10/22 17:04): The cutoffs chosen here are actually to some
# degree chosen as a balance between pointwise accuracy AND AD derivative
# accuracy. For example, the series gets worse faster for derivatives than it
# does for raw evals. Not entirely sure what to make of that, but here we are.
#
# TODO (cg 2021/10/27 10:28): the weak point is clearly with |z| between, say, 8
# and 15. We only get atols on the order of 1e-12 for that, which is not quite
# good enough to declare total victory. I can only manage to match the speed of
# AMOS in that part of the domain, too.
function _besselk(v, x, maxit=100, tol=1e-12, order=6)
  @assert x >= zero(x) DomainError("x needs to be non-negative.")
  iszero(x) && return Inf
  (real(v) < 0) && return _besselk(-v, x, maxit, tol, order)
  # TODO (cg 2021/11/16 16:44):  this is not the right way to test if you're
  # doing AD. But I don't want to hard-bake the ForwardDiff.Dual type in here.
  is_ad = !(v isa AbstractFloat)
  # Special cases, now just half-integers:
  #
  # TODO (cg 2021/12/16 12:41): for the moment, specifically at half
  # integer values the second derivatives are more accurate using the direct
  # series, even at values in which for direct evaluations v is large enough
  # that the series isn't so great. It isn't earth-shattering and the temme one
  # is much more expensive (300 ns vs 100 ns on my machine), but in the interest
  # of maximum safety I'm switching to this behavior.
  #
  # What would really be nice is some way if checking if (v isa Dual{T,V,N}
  # where T<: Dual). But again, I'm worried about getting too stuck with
  # ForwardDiff.
  if isinteger(v-1/2) && (x < 8.5)
    if is_ad
      return _besselk_ser(v, x, maxit, tol, false)
    else
      # Probably don't need the is_ad correction here.
      return _besselk_as(v, x, Int(ceil(v)), false)
    end
  end
  #
  # General cases:
  #
  # TODO (cg 2021/11/01 18:03): These branches are not perfect. If you go into
  # ./testing/accuracy.jl and track down the largest rtols between this code and
  # AMOS, you will be able to fiddle around with what version you use and get a
  # better rtol/atol. But I'm at the point where whenever I tweak something like
  # that, something else gets worse and makes it a wash. I think for the time
  # being I have to stop playing with things.
  if abs(x) < 8.5 # (x < 9)
    if (v > 2.85) || isnearint(v, 0.01) || ((x > 4) && is_ad)
      return _besselk_temme(v,x,maxit,tol,false) # direct series.
    else
      return _besselk_ser(v,x,maxit,tol,false) # direct series.
    end
  elseif abs(x) < 15.0
    return _besselk_asv(v,x,12) # uniform large order expn.
  elseif abs(x) < 30.0
    return _besselk_asv(v,x,8)
  else
    if abs(v) > 1.5 || is_ad
      return _besselk_asv(v,x,6)
    else
      return _besselk_as(v,x,order)
    end
  end
end

# A very simple wrapper that will use AMOS when possible, but fall back to this
# code. So now you can use AD on this but still get AMOS for direct evals.
#
# TODO (cg 2021/11/05 16:57): what's the most sensible naming thing here?
# Calling it besselk and not exporting it seems reasonable enough, but users
# will obviously want to import it. So not obvious what's best to do here.
function adbesselk(v, x, maxit=100, tol=1e-12, order=5)
  if v isa AbstractFloat
    SpecialFunctions.besselk(v, x)
  else
    _besselk(v, x, maxit, tol)
  end
end

# Not exactly a taylor series, but accurate enough.
#
# TODO (cg 2021/11/10 18:33): an enhancement here would be something that also
# worked for integers. But that seems hard. I did the whole Temme thing because
# integers are hard. But as it turns out, we really need it, because the Temme
# recursion depends on (x^v)*besselk(v,x) ->_{z->0} some finite number. But for
# v=0, which is needed to compute for _all_ integer v, that doesn't hold! An
# expansion like this that is valid for all v, including integer v, but is ALSO
# AD-compatible would take care of that entirely, because the K_{v+1}(x) term in
# the Temme series doesn't need f0.
@inline function besselkxv_t0(v, x)
  gv  = gamma(v)
  _2v = 2^(v-1)
  cof = (gv*_2v, zero(v), (_2v/4)*gv/(v-1), zero(v), (_2v/16)*gv/(v*v - 3*v + 2))
  evalpoly(x, cof)
end

# Note that the special series cutofs are low compared to the above function. In
# general, the adbesselk* functions that are exported really should be pretty
# near machine precision here or should just fall back to AMOS when possible. It
# turns out that the *xv modifications in the code are really only helpful when
# the argument is pretty small.
function adbesselkxv(v, x, maxit=100, tol=1e-12, order=5)
  (iszero(v) && iszero(x)) && return Inf
  is_ad = !(v isa AbstractFloat)
  xcut  = is_ad ? 6.0 : 2.0
  if !isinteger(v) && (abs(x) <= 1e-8) # use Taylor at zero.
    return besselkxv_t0(v, x)
  elseif (x < xcut) && (v < 5.75) && !isnearint(v, 0.01)
    return _besselk_ser(v, x, maxit, tol, true)
  elseif (x < xcut)
    return _besselk_temme(v, x, maxit, tol, true)
  elseif is_ad && (x > xcut) && (x < 15.0)
    return _besselk_asv(v, x, 12, true)
  else
    return adbesselk(v, x, maxit, tol, order)*(x^v)
  end
end

# Unlike adbesselkxv, this function is pure BesselK.jl, including the cases in
# which special branches for (x^v)*besselk(v,x) don't come up. This is primarily
# used for testing.
function _besselkxv(v, x, maxit=100, tol=1e-12, order=5)
  (iszero(v) && iszero(x)) && return Inf
  is_ad = !(v isa AbstractFloat)
  xcut  = is_ad ? 6.0 : 2.0
  if !isinteger(v) && (abs(x) <= 1e-8) # use Taylor at zero.
    return besselkxv_t0(v, x)
  elseif (x < xcut) && (v < 5.75) && !isnearint(v, 0.01)
    return _besselk_ser(v, x, maxit, tol, true)
  elseif (x < xcut)
    return _besselk_temme(v, x, maxit, tol, true)
  elseif is_ad && (x > xcut) && (x < 15.0)
    return _besselk_asv(v, x, 12, true)
  else
    return _besselk(v, x, maxit, tol, order)*(x^v)
  end
end


text:
export Obstacle, StateSpaceObstacle, ConfigSpaceObstacle, WorkspaceObstacle
export RobotCollisionBody, RobotCollisionComponent, RobotWorkspaceComponent, WorkspaceRobotBody, PointRobot
export CollisionChecker, ContinuousCollisionChecker, DiscreteCollisionChecker
export is_free_state, is_free_motion, is_free_config, is_free_sweep, rand_free_state, intersecting, sweep_intersecting

# Obstacle Sets
struct Obstacle{T,S}
    set::S
end
const StateSpaceObstacle{S}  = Obstacle{:statespace,S}
const ConfigSpaceObstacle{S} = Obstacle{:configspace,S}
const WorkspaceObstacle{S}   = Obstacle{:workspace,S}
StateSpaceObstacle(set::S) where {S}  = Obstacle{:statespace,S}(set)
ConfigSpaceObstacle(set::S) where {S} = Obstacle{:configspace,S}(set)
WorkspaceObstacle(set::S) where {S}   = Obstacle{:workspace,S}(set)
intersecting(o::Obstacle, x) = intersecting(o.set, x)
sweep_intersecting(o::Obstacle, x0, xf) = sweep_intersecting(o.set, x0, xf)
@recipe function f(o::Obstacle; dims=(1, 2))#, obstacle_color=:red, obstacle_alpha=1)
    dims  --> dims
    # color :=  obstacle_color
    # alpha :=  obstacle_alpha
    # label --> ""
    o.set
end

# Abstract Robot Types (speculative for now)
abstract type RobotCollisionBody end
abstract type RobotCollisionComponent end
abstract type RobotWorkspaceComponent <: RobotCollisionComponent end

# Collision Checker
struct CollisionChecker{C,O,R<:RobotCollisionBody,S2C}
    obstacles::O
    robot::R
    state2config::S2C
    motion_count::Base.RefValue{Int}
    edge_count::Base.RefValue{Int}
end
const ContinuousCollisionChecker{O,R,S2C} = CollisionChecker{true,O,R,S2C}
const DiscreteCollisionChecker{O,R,S2C}   = CollisionChecker{false,O,R,S2C}
@recipe function f(CC::CollisionChecker; dims=(1, 2))#, obstacle_color=:red, obstacle_alpha=1)
    dims           --> dims
    # obstacle_color --> obstacle_color
    # obstacle_alpha --> obstacle_alpha
    foreach(o -> @series(begin o end), CC.obstacles)
end

# General Methods for Dispatch
## is_free_state
is_free_state(CC::CollisionChecker, x) = all(o -> is_free_state(o, CC.robot, CC.state2config, x), CC.obstacles)
is_free_state(o::Obstacle, robot, s2c, x) = is_free_config(o, robot, s2c(x))
is_free_state(o::StateSpaceObstacle,  robot, s2c, x) = !intersecting(o, x)

## is_free_motion
function is_free_motion(CC::CollisionChecker, x0, xf)
    CC.motion_count[] += 1
    all(o -> is_free_motion(o, CC.robot, CC.state2config, x0, xf), CC.obstacles)
end
is_free_motion(o::Obstacle, robot, s2c, x0, xf) = is_free_sweep(o, robot, s2c(x0), s2c(xf))
is_free_motion(o::StateSpaceObstacle,  robot, s2c, x0, xf) = !sweep_intersecting(o, x0, xf)

## is_free_config
is_free_config(o::ConfigSpaceObstacle, robot, q) = !intersecting(o, q)

## is_free_sweep
is_free_sweep(o::ConfigSpaceObstacle, robot, q0, qf) = !sweep_intersecting(o, q0, qf)

## rand_free_state
rand_free_state(CC::CollisionChecker, X) = (x = rand(X); is_free_state(CC, x) ? x : rand_free_state(CC, X))

# WorkspaceRobotBody
struct WorkspaceRobotBody{C2W} <: RobotCollisionBody
    config2workspace::C2W
end
is_free_config(o::WorkspaceObstacle, robot::WorkspaceRobotBody, q) = !intersecting(o, robot.config2workspace(q))
function is_free_sweep(o::WorkspaceObstacle, robot::WorkspaceRobotBody, q0, qf)
    !sweep_intersecting(o, robot.config2workspace(q0), robot.config2workspace(qf))
end

PointRobot() = WorkspaceRobotBody(identity)

# SingleLinkRobot
struct SingleLinkRobot{L,C2T} <: RobotCollisionBody
    link::L
    config2transform::C2T
end
function is_free_config(o::WorkspaceObstacle, robot::SingleLinkRobot, q)
    !intersecting(o, robot.link, robot.config2transform(q))
end
function is_free_sweep(o::WorkspaceObstacle, robot::SingleLinkRobot, q0, qf)
    !sweep_intersecting(o, robot.link, robot.config2transform(q0), robot.config2transform(qf))
end

## Convenience Constructors
_type_sort(obstacles) = obstacles
_type_sort(obstacles::Vector{T}) where {T} = isconcretetype(T) ? obstacles : TypeSortedCollection(obstacles)
function ContinuousCollisionChecker(_obstacles, robot::R=PointRobot(), state2config::S2C=identity) where {R,S2C}
    obstacles = _type_sort(_obstacles)
    CollisionChecker{true,typeof(obstacles),R,S2C}(obstacles, robot, state2config, Ref(0), Ref(0))
end
function DiscreteCollisionChecker(_obstacles, robot::R=PointRobot(), state2config::S2C=identity) where {R,S2C}
    obstacles = _type_sort(_obstacles)
    CollisionChecker{false,typeof(obstacles),R,S2C}(obstacles, robot, state2config, Ref(0), Ref(0))
end

# Super basic collision checking; `using SeparatingAxisTheorem2D` or ... to load more advanced capabilities
intersecting(X::Union{SimpleConvexSet,SetComplement}, x) = x in X

function sweep_intersecting(B::AxisAlignedBox{N}, x0, xf) where {N}
    all(B.lo[i] <= max(x0[i], xf[i]) && min(x0[i], xf[i]) <= B.hi[i] for i in 1:N) && _sweep_intersecting(B, x0, xf)
end
function _sweep_intersecting(B::AxisAlignedBox{N}, x0, xf) where {N}
    any((λ = (ifelse(x0[i] < B.lo[i], B.lo[i], B.hi[i]) - x0[i])/(xf[i] - x0[i]);
        0 <= λ <= 1 && all(i == j || B.lo[j] <= x0[j] + (xf[j] - x0[j])*λ <= B.hi[j] for j in 1:N)) for i in 1:N)
end

function sweep_intersecting(B::Ball, x0, xf)
    all(B.c[i] - B.r <= max(x0[i], xf[i]) &&
        min(x0[i], xf[i]) <= B.c[i] + B.r for i in 1:N) && _sweep_intersecting(B, x0, xf)
end
function _sweep_intersecting(B::Ball, x0, xf)
    x0 in B && return true
    xf in B && return true
    X = xf - x0
    Y = B.c - x0
    a = dot(X, X)
    b = dot(Y, Y)
    c = dot(X, Y)
    a*C.r*C.r < a*b - c*c && return false
    0 <= c <= a
end

sweep_intersecting(S::SetComplement{<:SimpleConvexSet}, x0, xf) = x0 in S || xf in S    # S.set is convex

max_stars_repo_name: pbpf/clean

max_issues_repo_name: pbpf/clean

max_forks_repo_name: pbpf/clean


text:
function  dbdt(r,x)

#   constructs sparse matrix Nsp for the time derivative equation
#
#   input
#   r        - radial grid (excludes boundaries)
#   x        - cos(theta) grid (excludes poles and equator)

   nr = length(r);
   nx = length(x);

   # initialize arrays
   ndim = nx * nr;
   row = Array{Int64}(undef,ndim);
   col = Array{Int64}(undef,ndim);
   val = Array{ComplexF64}(undef,ndim);

   #
   # time derivative of b_phi
   count = 0;
   for l = 1 : nx
       for k = 1 : nr

           count = count + 1;
           ipos = (l-1) * nr + k;
           row[count] = ipos;
           col[count] = ipos + ndim;
           val[count] = 1;
       end
   end

   nsp = sparse(row,col,val,2*ndim,2*ndim);

end

max_stars_repo_name: bbuffett/MacWaves

max_issues_repo_name: bbuffett/MacWaves

max_forks_repo_name: bbuffett/MacWaves


text:
if length(lowercase(strip(get(ENV, "PREDICTMD_TEST_PLOTS", "")))) == 0
    ENV["PREDICTMD_TEST_PLOTS"] = "false"
end

if isa(TEST_GROUP, TestGroupAll)
    ENV["PREDICTMD_TEST_PLOTS"] = "true"
end

if isa(TEST_GROUP, TestGroupTestPlots)
    ENV["PREDICTMD_TEST_PLOTS"] = "true"
end

if isa(TEST_GROUP, TestGroupTravis1)
    ENV["PREDICTMD_TEST_PLOTS"] = "true"
end
if isa(TEST_GROUP, TestGroupTravis2)
    ENV["PREDICTMD_TEST_PLOTS"] = "true"
end
if isa(TEST_GROUP, TestGroupTravis3)
    ENV["PREDICTMD_TEST_PLOTS"] = "true"
end
if isa(TEST_GROUP, TestGroupTravis4)
    ENV["PREDICTMD_TEST_PLOTS"] = "true"
end
if isa(TEST_GROUP, TestGroupTravis5)
    ENV["PREDICTMD_TEST_PLOTS"] = "true"
end
if isa(TEST_GROUP, TestGroupTravis6)
    ENV["PREDICTMD_TEST_PLOTS"] = "true"
end

if isa(TEST_GROUP, TestGroupDocker1)
    ENV["PREDICTMD_TEST_PLOTS"] = "true"
end
if isa(TEST_GROUP, TestGroupDocker2)
    ENV["PREDICTMD_TEST_PLOTS"] = "true"
end
if isa(TEST_GROUP, TestGroupDocker3)
    ENV["PREDICTMD_TEST_PLOTS"] = "true"
end
if isa(TEST_GROUP, TestGroupDocker4)
    ENV["PREDICTMD_TEST_PLOTS"] = "true"
end

@info(
    string(
        "PREDICTMD_TEST_PLOTS: \"",
        ENV["PREDICTMD_TEST_PLOTS"],
        "\"",
        )
    )


text:
# Inspired by https://mtk.sciml.ai/dev/tutorials/acausal_components/

# a point in the circuit with a given voltage and current
@connector function Pin(; name)
    sts = @variables v(t)=1.0 i(t)=1.0
    ODESystem(Equation[], t, sts, []; name)
end

# ground pin, where voltage is 0
function Ground(; name)
    @named g = Pin()
    eqs = [g.v ~ 0.]
    compose(ODESystem(eqs, t, [], []; name), g)
end

# abstraction for all 2-pin components
function OnePort(; name)
    # two pins
    @named p = Pin()
    @named n = Pin()
    # component has a voltage across it, and current through
    sts = @variables v(t)=1.0 i(t)=1.0
    eqs = [
        v ~ p.v - n.v    # KVL
        0. ~ p.i + n.i   # KCL
        i ~ p.i          # Current through component is current through +ve pin
    ]
    compose(ODESystem(eqs, t, sts, []; name), p, n)
end

function ModelingToolkit.connect(::Type{Pin}, ps...)
    eqs = [
        0. ~ sum(p->p.i, ps)    # KCL
    ]
    # KVL
    for i in 1:length(ps)-1
        push!(eqs, ps[i].v ~ ps[i+1].v)
    end
    return eqs
end

connect_star(comps...) = connect([comp.n for comp in comps]...)

function connect_series(comps...)
    eqs = Equation[]
    for i in 1:length(comps)-1
        eqs = vcat(eqs, connect(comps[i].p, comps[i+1].n))
    end
    return eqs
end

function connect_parallel(comps...)
    return [
        connect((comp.p for comp in comps)...)
        connect((comp.n for comp in comps)...)
    ]
end


text:
# This Source Code Form is subject to the terms of the Mozilla Public
# License, v. 2.0. If a copy of the MPL was not distributed with this
# file, You can obtain one at http://mozilla.org/MPL/2.0/.

module SubmodularUtilities

using LinearAlgebra

export lazy_greedy
using Base.Order
using DataStructures
"""
    lazy_greedy(f, ground_set, k)

`ground_set` should be an array of unique integers, which denotes the ground
set of monotone submodular function `f`. `f` represents a monotone submodular
function. If `T` is an array whose elements constitute a subset of `ground_set`,
`f(T)` is the function value of `T`. `k` is the cardinality constraint.
`lazy_greedy` outputs an array of length `k` using the lazy greedy algorithm.
"""
function lazy_greedy(f, ground_set, k)
    uppers = [(Inf, i) for i in ground_set]
    myorder = ReverseOrdering(By(x->x[1]))
    heapify!(uppers, myorder)
    s = []
    current = 0.
    for i in 1:k
        while true
            v, e = heappop!(uppers, myorder)
            new_upper_bound = f(vcat(s, e)) - current
            if new_upper_bound >= uppers[1][1]
                push!(s, e)
                current += new_upper_bound
                break
            else
                heappush!(uppers, (new_upper_bound, e), myorder)
            end
        end
    end
    return s
end

mutable struct FPL
    x
    eta
    accumulated_gradient
    n0
end

export Constraint
"""
    Constraint(A, sense, b, lb, ub)

The constraints are A*x sense b and lb <= x <= ub.
`sense` is a vector of constraint sense characters `<`, `=`, and `>`.
`A` is a 2-dimensional array, which denotes the constraint matrix.
`b` is the right-hand side vector.
"""
struct Constraint
    A
    sense
    b
    lb
    ub
end

using MathProgBase, Clp
function update(fpl, gradient, constraint)
    fpl.accumulated_gradient .-= gradient
    sol = linprog(fpl.eta * fpl.accumulated_gradient + fpl.n0,
        constraint.A, constraint.sense, constraint.b, constraint.lb, constraint.ub, ClpSolver())
    if sol.status == :Optimal
        fpl.x = sol.sol
    else
        error("No solution was found.")
    end
end

function get_vector(fpl)
    return fpl.x
end

export meta_frank_wolfe
"""
    meta_frank_wolfe(dim_x, f_list, g_list, eta, K, constraint; variance_reduction = false, n0 = rand(dim_x))

Args:
- `dim_x`: The dimension of the variables that we optimize.
- `f_list`: List of continuous DR-submodular functions `[f_1, f_2, ..., f_T]`.
- `g_list`: List of (unbiased estimates of) the gradients of functions in `f_list`.
- `eta`: This implementation uses Follow-Perturbed-Leader (FPL) for linear losses (which can be found in Hazan's book Introduction to Online Convex Optimization) as the online linear maximization algorithm. `eta` is the parameter of FPL.
- `K`: Number of instances of the chosen online linear maximization algorithm.
- `constraint`: Constraint set. It is an object of type `Constraint`.
- `variance_reduction`: Boolean variable which is `true` if variance reduction is enabled
- `n0`: the value of n0 in FPL for linear losses.

Returns:
- Array of rewards `[f_1(x_1), f_2(x_2), ..., f_T(f_T)]`.
"""
function meta_frank_wolfe(dim_x, f_list, g_list, eta, K, constraint;
    variance_reduction = false, n0 = rand(dim_x))
    sol = linprog(-n0,
        constraint.A, constraint.sense, constraint.b, constraint.lb, constraint.ub, ClpSolver())
    if sol.status != :Optimal
        error("No solution was found.")
    end
    x0 = sol.sol
    fpl = [FPL(x0, eta, zeros(dim_x), n0) for k in 1:K]
    @assert length(f_list) == length(g_list)
    T = length(f_list)
    reward = zeros(T)

    for t in 1:T
        v = [get_vector(fpl[k]) for k in 1:K]
        x = zeros(K + 1, dim_x)
        for k in 1:K
            x[k + 1, :] = x[k, :] + v[k] / K
        end
        reward[t] = f_list[t](x[K + 1, :])
        d = zeros(dim_x)
        for k in 1:K
            if variance_reduction
                rho = 0.5 * k^(-2/3)
                d = (1 - rho) * d + rho * g_list[t](x[k, :])
            else
                d = g_list[t](x[k, :])
            end
            update(fpl[k], d, constraint)
        end
    end
    return reward
end

using Ipopt
export get_projection_operator
"""
    get_projection_operator(constraint; dim_x = size(constraint.A)[2])

Args:
- `constraint`: Constraint set. It is an object of type `Constraint`.
- `dim_x`: the dimension of the point. The default value is deduced from the size of `constraint.A`.

Returns:
- A function `projection` such that `projection(x)` outputs the projection of `x`
onto the set `constraint`.
"""
function get_projection_operator(constraint; dim_x = size(constraint.A)[2])
    function projection(x0)
        sol = quadprog(-x0, eye(dim_x), constraint.A, constraint.sense, constraint.b, constraint.lb, constraint.ub, IpoptSolver(print_level=0))
        if sol.status == :Optimal
            return sol.sol
        end
        error("No projection of $x0 is found.")
    end
    return projection
end

export online_gradient_ascent
"""
    online_gradient_ascent(x0, f_list, g_list, projection; step_size = t->1/t)

Args:
- `x0`: the initial value
- `f_list`: List of continuous DR-submodular functions `[f_1, f_2, ..., f_T]`.
- `g_list`: List of (unbiased estimates of) the gradients of functions in `f_list`.
- `projection`: Projection operator of the constraint set. You may want to get it by calling `get_projection_operator`.
- `step_size`: A function such that `step_size(t)` is the step size at the t-th iteration.

Returns:
- Array of rewards `[f_1(x_1), f_2(x_2), ..., f_T(f_T)]`.
"""
function online_gradient_ascent(x0, f_list, g_list, projection; step_size = t->1/t)
    @assert length(f_list) == length(g_list)
    T = length(f_list)
    reward = zeros(T)
    x = copy(x0)
    for t in 1:T
        reward[t] = f_list[t](x)
        x = projection(x + g_list[t](x) * step_size(t))
    end
    return reward
end

export pipage_round
"""
    pipage_round(x)

Given an array `x` whose every entry is between 0 and 1, round x into a binary vector
using the pipage rounding algorithm presented in [^CCPV].

[^CCPV]: Calinescu, Gruia, et al. "Maximizing a monotone submodular function subject to a matroid constraint." SIAM Journal on Computing 40.6 (2011): 1740-1766.
"""
function pipage_round(x)
    function pipage_round_raw(x)
        y = copy(x)
        T = []
        for i in 1:length(y)
            if 0 < y[i] < 1
                push!(T, i)
            end
        end

        try
            while length(T) > 0
                i, j = T[1], T[2]
                if y[i] + y[j] < 1
                    p = y[j] / (y[i] + y[j])
                    if rand() < p
                        y[j] += y[i]
                        y[i] = 0.
                        deleteat!(T, 1)
                    else
                        y[i] += y[j]
                        y[j] = 0.
                        deleteat!(T, 2)
                    end
                else
                    p = (1. - y[i]) / (2. - y[i] - y[j])
                    if rand() < p
                        y[i] += y[j] - 1.
                        y[j] = 1.
                        deleteat!(T, 2)
                        if y[i] == 0.
                            deleteat!(T, 1)
                        end
                    else
                        y[j] += y[i] - 1.
                        y[i] = 1.
                        deleteat!(T, 1)
                        if y[j] == 0.
                            deleteat!(T, 1)
                        end
                    end
                end
            end
        catch
            return y
        end
        return y
    end
    function regularize_answer(ans)
        pp = pipage_round_raw(ans)
        pp[pp .< 1] = 0.
        return pp
    end
    px = regularize_answer(x)
    v = [i for i in 1:length(x) if px[i] != 0]
    return v
end

export random_round
"""
    random_round(x)

Given an array `x` whose every entry is between 0 and 1, round `x` into a binary vector.
The i-th entry of the output vector is 1 with probability `x[i]` and is 0 otherwise.
"""
function random_round(x)
    p = rand(size(x));
    return find(p .< x)
end

export get_random_evaluation_of_multilinear_extension
"""
    get_random_evaluation_of_multilinear_extension(f_discrete)

The input `f_discrete` is assumed to be a monotone submodular function defined
on the ground set `1:n`, where n is some positive integer. Given an array `T`
whose elements are chosen from `1:n`, `f_discrete(T)` is the function value of
the subset `T`. `get_random_evaluation_of_multilinear_extension` returns a
function that takes a size-n array `x` as input and outputs an unbiased estimate of
the function value of the multilinear extension of `f_discrete` at `x`. In other
 words, the output is a function that maps `x` to `f_discrete(random_round(x))`.
"""
function get_random_evaluation_of_multilinear_extension(f_discrete)
    return x->f_discrete(random_round(x))
end

export get_random_gradient_of_multilinear_extension
"""
    get_random_gradient_of_multilinear_extension(f_discrete)

The input `f_discrete` is assumed to be a monotone submodular function defined
on the ground set `1:n`, where n is some positive integer. Given an array `T`
whose elements are chosen from `1:n`, `f_discrete(T)` is the function value of
the subset `T`. `get_random_gradient_of_multilinear_extension` returns a
function that takes a size-n array `x` as input and outputs an unbiased estimate of
the gradient of the multilinear extension of `f_discrete` at `x`.

In other words, the output function first computes `S = random_round(x)`. Then
the i-th partial derivative is
`f_discrete(vcat(S, i)) - f_discrete(setdiff(S, i))`.
"""
function get_random_gradient_of_multilinear_extension(f_discrete)
    function stochastic_gradient(x)
        S = random_round(x);
        grad = zeros(x);
        for i in 1:length(x)
            grad[i] = f_discrete(vcat(S, i)) - f_discrete(setdiff(S, i));
        end
        return grad
    end
    return stochastic_gradient
end

export get_function_exemplar_based_clustering
"""
    get_function_exemplar_based_clustering(data)

Returns the objective function of exemplar-based clustering defined on the
ground set `1:size(data)[1]`. The argument `data` is a 2-dimensional array.
Each row represents a data point and each column represents an attribute.
"""
function get_function_exemplar_based_clustering(data)
    function f_exemplar(S)
        n_V = size(data)[1];
        function L(S)
            return sum([minimum([norm(data[e, :] - (v > 0 ? data[v, :] : 0))^2 for v in S])
                for e in 1:n_V]) / n_V
        end
        return L([0]) - L(vcat(S, 0))
    end
    return f_exemplar
end

export get_function_active_set_selection
"""
    get_function_active_set_selection(data; sigma = 1., h = 0.75)

Returns the objective function of active set selection defined on the ground set
`1:size(data)[2]`. The argument `data` is a 2-dimensional array.
Each row represents a data point and each column represents an attribute.
The covariance matrix uses a Gaussian kernel
`exp(-norm(data[:, i] - data[:, j])^2 / h^2)`. The function value of the
objective function evaluated at the subset `S` is
`0.5 * log(det(eye(length(S)) + cov_matrix[S, S] / sigma^2))`.
"""
function get_function_active_set_selection(data; sigma = 1., h = 0.75)
    n_attr = size(data)[2];
    cov_matrix = zeros(n_attr, n_attr);
    for i in 1:n_attr
        for j in 1:i
            cov_matrix[i, j] = exp(-norm(data[:, i] - data[:, j])^2 / h^2)
            if j != i
                cov_matrix[j, i] = cov_matrix[i, j]
            end
        end
    end
    function f_active_set_discrete(S)
        if isempty(S)
            return 0.;
        end
        return 0.5 * log(det(eye(length(S)) + cov_matrix[S, S] / sigma^2))
    end
    return f_active_set_discrete
end
end # end of module SubmodularUtilities


text:
using REPL

export
    # wrapping
    CSI,
    displaysize,
    raw!,
    cmove_up,
    cmove_down,
    cmove_left,
    cmove_right,
    cmove_line_up,
    cmove_line_down,
    cmove_col,
    clear,
    clear_line,
    beep,
    enable_bracketed_paste,
    disable_bracketed_paste,
    end_keypad_transmit_mode,

    # extensions
    displaysize,
    cmove,
    clear_line,
    cmove_line_last,
    cshow,
    csave,
    crestore,
    alt_screen,

    # utils
    read_next_char,
    init_term,
    read_stream_bytes,
    read_stream,
    flush,
    buffered

# +----------------------------+
# | wrapping of REPL.Terminals |
# +----------------------------+

const CSI = REPL.Terminals.CSI

displaysize(t=term) = REPL.Terminals.displaysize(t)
raw!(t::REPL.Terminals.UnixTerminal, enable::Bool) = REPL.Terminals.raw!(t, enable)
raw!(enable::Bool) = raw!(term, enable)

cmove_up(stream::IO, n::Int) = Base.write(stream, "$(CSI)$(n)A")
cmove_up(n::Int) = cmove_up(out_stream, n)
cmove_up(stream::IO) = cmove_up(stream, 1)
cmove_up()  = cmove_up(1)

cmove_down(stream::IO, n::Int) = Base.write(stream, "$(CSI)$(n)B")
cmove_down(n::Int) = cmove_down(out_stream, n)
cmove_down(stream::IO) = cmove_down(stream, 1)
cmove_down() = cmove_down(1)

cmove_right(stream::IO, n::Int) = Base.write(stream, "$(CSI)$(n)C")
cmove_right(n::Int) = cmove_right(out_stream, n)
cmove_right(stream::IO) = cmove_right(stream, 1)
cmove_right() = cmove_right(1)

cmove_left(stream::IO, n::Int) = Base.write(stream, "$(CSI)$(n)D")
cmove_left(n::Int) = cmove_left(out_stream, n)
cmove_left(stream::IO) = cmove_left(stream, 1)
cmove_left() = cmove_left(1)

cmove_col(stream::IO, n::Int) = (Base.write(stream, '\r'); n > 1 && cmove_right(stream, n-1)) # SCI n G
cmove_col(n::Int) = cmove_col(out_stream, n)
cmove_col(stream::IO) = cmove_col(stream, 1)
cmove_col() = cmove_col(1)

cmove_line_up(stream::IO, n::Int) = (cmove_up(stream, n); cmove_col(stream)) # CSI n F
cmove_line_up(n::Int) = cmove_line_up(out_stream, n)
cmove_line_up(stream::IO) = cmove_line_up(stream, 1)
cmove_line_up() = cmove_line_up(1)

cmove_line_down(stream::IO, n::Int) = (cmove_down(stream, n); cmove_col(stream)) # SCI n E
cmove_line_down(n::Int) = cmove_line_down(out_stream, n)
cmove_line_down(stream::IO) = cmove_line_down(stream, 1)
cmove_line_down() = cmove_line_down(1)

@eval clear(stream::IO) = Base.write(stream, $"$(CSI)H$(CSI)2J")
clear() = clear(out_stream)

@eval clear_line(stream::IO) = Base.write(stream, $"\r$(CSI)0K")
clear_line() = clear_line(out_stream)

beep(stream::IO) = Base.write(stream,"\x7")
beep() = beep(err_stream)

@eval enable_bracketed_paste(stream::IO) = Base.write(stream, $"$(CSI)?2004h")
enable_bracketed_paste() = enable_bracketed_paste(out_stream)

@eval disable_bracketed_paste(stream::IO) = Base.write(stream, $"$(CSI)?2004l")
disable_bracketed_paste() = disable_bracketed_paste(out_stream)

@eval end_keypad_transmit_mode(stream::IO) = Base.write(stream, $"$(CSI)?1l\x1b>")
end_keypad_transmit_mode() = end_keypad_transmit_mode(out_stream)

# +------------+
# | extensions |
# +------------+

displaysize(stream::IO, height::Int, width::Int) = Base.write(stream, "$(CSI)8;$(height);$(width)t")
displaysize(height::Int, width::Int) = displaysize(out_stream, height, width)

cmove(stream::IO, y::Int, x::Int) = Base.write(stream, "$(CSI)$(y);$(x)H")
cmove(y::Int, x::Int) = cmove(out_stream, y, x)

cmove_line_last(stream::IO) = Base.write(stream, "$(CSI)$(displaysize()[1]);1H")
cmove_line_last() = cmove_line_last(out_stream)

clear_line(stream::IO, row::Int) = (Base.write(stream, "$(CSI)$(row);1H"); clear_line())
clear_line(row::Int) = clear_line(out_stream, row)

cshow(stream::IO, enable=true) = (enable ? Base.write(stream, "$(CSI)?25h") : Base.write(stream, "$(CSI)?25l"))
cshow(enable=true) = cshow(out_stream, enable)

csave(stream::IO) = Base.write(stream, "$(CSI)s")
csave() = csave(out_stream)

crestore(stream::IO) = Base.write(stream, "$(CSI)u")
crestore() = crestore(out_stream)

alt_screen(stream::IO, enable::Bool) = (enable ? Base.write(stream, "$(CSI)?1049h") : Base.write(stream, "$(CSI)?1049l"))
alt_screen(enable::Bool) = alt_screen(out_stream, enable)

# +----+
# | IO |
# +----+

write(stream::IO, args...) = Base.write(stream, args...)
write(args...) = Base.write(out_stream, args...)

print(stream::IO, args...) = Base.print(stream, args...)
print(args...) = Base.print(out_stream, args...)

println(stream::IO, args...) = Base.println(stream, args...)
println(args...) = Base.println(out_stream, args...)

join(stream::IO, args...) = Base.join(stream, args...)
join(args...) = Base.join(out_stream, args...)

# +-------+
# | utils |
# +-------+

function init_term(in_stream=stdin, out_stream=stdout, err_stream=stderr)
    return REPL.Terminals.TTYTerminal(
        get(ENV, "TERM", Sys.iswindows() ? "" : "dumb"),
        in_stream, out_stream, err_stream
    )
end

read_next_byte(io::IO) = read(io, 1)[1]

function read_stream_bytes(stream::IO)
    queue = UInt8[]

    push!(queue, read_next_byte(stream))

    stream_size = stream.buffer.size
    if stream_size > 1
        for i=1:stream_size-1
            push!(queue, read_next_byte(stream))
        end
    end

    return queue
end

read_next_char(io::IO) = Char(read_next_byte(io))

read_stream(stream::IO) = String(read_stream_bytes(stream))
read_stream() = read_stream(in_stream)

flush(stream::IO, buffer::IOBuffer) = Base.write(stream, take!(buffer))

function buffered(f, stream::IO, argv...)
    buffer=IOBuffer()
    f(buffer, argv...)
    flush(stream, buffer)

    return
end

buffered(f, argv...) = buffered(f, out_stream, argv...)


text:
MissingFloat=Union{Missing,Float64}

function mean1(x)
    out=zeros(size(x,2))
    for i=1:size(x,2)
        out[i]=mean(skipmissing(x[:,i]))
    end
    return out
end

function std1(x)
    out=zeros(size(x,2))
    for i=1:size(x,2)
        out[i]=std(skipmissing(x[:,i]))
    end
    return out
end

function mean2(x)
    out=zeros(size(x,1))
    for i=1:size(x,1)
        out[i]=mean(skipmissing(x[i,:]))
    end
    return out
end

small_value = 1e-12

function F1(N)
    # N[i,j] is the number of times model predicts class i and true is class j
    # I use class 1 for the "positive" class and class 2 for the "negative" class
    true_positive = N[1,1]
    false_positive = N[1,2]
    false_negative = N[2,1]
    return true_positive/(true_positive + 0.5*(false_positive+false_negative))
end

function Precision(N)
    # N[i,j] is the number of times model predicts class i and true is class j
    # I use class 1 for the "positive" class and class 2 for the "negative" class
    true_positive = N[1,1]
    false_positive = N[1,2]
    false_negative = N[2,1]
    return true_positive/(true_positive + false_positive)
end

function Recall(N)
    # N[i,j] is the number of times model predicts class i and true is class j
    # I use class 1 for the "positive" class and class 2 for the "negative" class
    true_positive = N[1,1]
    false_positive = N[1,2]
    false_negative = N[2,1]
    return true_positive/(true_positive + false_negative)
end

function fitBeta(mean_val,var_val)
    mn=mean_val
    s2=min(var_val,0.5^2)
    a=(1-mn)*mn^2/s2 - mn
    b=a*(1/mn-1)
    return a,b
end

function confidenceBeta(mean_fhat,var_fhat,confidence)
#    println([mean_fhat,var_fhat,confidence])
    if ismissing(mean_fhat)
        return missing,missing,missing
    end
    if var_fhat<1e-20
        return mean_fhat,mean_fhat,mean_fhat
    end
    var_fhat = var_fhat + small_value # try to avoid numerical issues

    tail=0.5*(1-confidence) # tail probabilities (for each tail)

    a,b = fitBeta(mean_fhat,var_fhat)

    a = max(a,0.001) # ugly hack to avoid problems
    b = max(b,0.001) # ugly hack to avoid problems

    dist=Beta(a,b)
    mn = a/(a+b) # mean posterior error

    upper_limit = quantile(dist,confidence+tail)
    lower_limit = quantile(dist,tail)

    return mean_fhat, lower_limit, upper_limit
end

function f(c_pred,c_true;metric,metric_par=0)
    @match metric begin
        "F" =>  return  (c_pred .==1) .* (c_true .==1)
        "Accuracy" => return  (c_pred .==1) .* (c_true .==1) + (c_pred .==0) .* (c_true .==0)
        "Specificity" => return  (c_pred .==0) .* (c_true .==0)
    end
end

function g(c_pred,c_true;metric,metric_par=0)
    alpha=metric_par
    @match metric begin
        "F" =>  return alpha.*(c_pred .==1) + (1 - alpha) .* (c_true .== 1)
        "Accuracy" => return ones(size(c_pred))
        "Specificity" => return (c_true .==0)
    end
end

micro_add=1e-10

function f_micro(c_pred,c_true)
    #return  sum((c_pred .==1) .* (c_true .==1),dims=2)
    return  (c_pred .==1) .* (c_true .==1) .+ micro_add
end

function g_micro(c_pred,c_true)
    alpha=0.5
    #return sum(alpha.*(c_pred .==1) + (1 - alpha) .* (c_true .== 1),dims=2)
    return alpha.*(c_pred .==1) + (1 - alpha) .* (c_true .== 1) .+ micro_add
end

function fgtilde(c_pred::Vector{Int};metric,metric_par=0)
    ftilde1 = f(c_pred,ones(size(c_pred));metric,metric_par)
    gtilde1 = g(c_pred,ones(size(c_pred));metric,metric_par)
    ftilde0 = f(c_pred,zeros(size(c_pred));metric,metric_par)
    gtilde0 = g(c_pred,zeros(size(c_pred));metric,metric_par)
    return (ftilde1,ftilde0,gtilde1,gtilde0)
end

function FhatIS(c_pred::Vector{Int},c_true::Vector{Int},q;metric,metric_par=0)
    ftilde = f(c_pred,c_true;metric,metric_par)
    gtilde = g(c_pred,c_true;metric,metric_par)
    Fapprox = (small_value + sum(ftilde./q))./(small_value + sum(gtilde./q))
    mn = CheckMissing(Fapprox)

    M = length(c_pred) # number of samples
    (ftilde1,ftilde0,gtilde1,gtilde0) = fgtilde(c_pred;metric,metric_par)

    lambda=0.999 # NEED TO CHECK THIS!!!!
    c_true_s=lambda*c_true .+ (1-lambda)*0.5
    gtilde = c_true_s.*gtilde1 + (1 .- c_true_s).*gtilde0
    htilde2 = c_true_s.*(ftilde1 -Fapprox.*gtilde1).^2 + (1 .- c_true_s).*(ftilde0 .- Fapprox.*gtilde0).^2

    yhat = mean( gtilde./q)/Ntest
    numhat = mean( htilde2./(q.*q))/Ntest

    var = (small_value.^2 .+ numhat) ./ (small_value.^2 .+  yhat^2)/ (M*Ntest)

    if ismissing(mn) | ismissing(var)
        return missing,missing
    else
        return mn,var
    end

end

function FhatBS(c_pred::Vector{Int},c_true::Vector{Int},b;metric,metric_par=0)
    ftilde = f(c_pred,c_true;metric,metric_par)
    gtilde = g(c_pred,c_true;metric,metric_par)
    Fapprox = (small_value + sum(ftilde./b))./(small_value + sum(gtilde./b))
    mn  = CheckMissing(Fapprox)

    (ftilde1,ftilde0,gtilde1,gtilde0) = fgtilde(c_pred;metric,metric_par)

    lambda=0.999 # slightly smooth to avoid numerical issues
    c_true_s=lambda*c_true .+ (1-lambda)*0.5
    gtilde = c_true_s.*gtilde1 + (1 .- c_true_s).*gtilde0
    htilde2 = c_true_s.*(ftilde1 -Fapprox.*gtilde1).^2 + (1 .- c_true_s).*(ftilde0 .- Fapprox.*gtilde0).^2
    htilde = c_true_s.*(ftilde1 -Fapprox.*gtilde1) + (1 .- c_true_s).*(ftilde0 .- Fapprox.*gtilde0)

    var = sum( (htilde2./b .- htilde.^2)./b) / (sum(gtilde./b)).^2

    if ismissing(mn) | ismissing(var)
        return missing,missing
    else
        return mn,var
    end
end

function FhatBS_micro(c_pred,c_true,b)
    ftilde = sum(f_micro(c_pred,c_true),dims=2)
    gtilde = sum(g_micro(c_pred,c_true),dims=2)
    Fapprox = (small_value + sum(ftilde./b))./(small_value + sum(gtilde./b))
    mn  = CheckMissing(Fapprox)

    lambda=0.999
    c_true_s=lambda*c_true .+ (1-lambda)*0.5

    ftilde1 = sum(f_micro(c_pred,ones(size(c_true))),dims=2)
    ftilde0 = sum(f_micro(c_pred,zeros(size(c_true))),dims=2)
    gtilde1 = sum(g_micro(c_pred,ones(size(c_true))),dims=2)
    gtilde0 = sum(g_micro(c_pred,zeros(size(c_true))),dims=2)

    htilde2 = c_true_s.*(ftilde1 -Fapprox.*gtilde1).^2 + (1 .- c_true_s).*(ftilde0 .- Fapprox.*gtilde0).^2
    htilde = c_true_s.*(ftilde1 -Fapprox.*gtilde1) + (1 .- c_true_s).*(ftilde0 .- Fapprox.*gtilde0)

    gtilde = c_true.*gtilde1 + (1 .- c_true).*gtilde0

    var = sum( (htilde2./b .- htilde.*htilde)./b) / (sum(gtilde./b)).^2

    if ismissing(mn) | ismissing(var)
        return missing,missing
    else
        return mn,var
    end
end

function FhatBS_macro(c_pred,c_true,b)

    den = length(b)

    a_bar = vec(sum(  ((c_pred.==1) .& (c_true.==1))./b, dims=1)./den)
    b_bar = vec(sum(  ((c_pred.==1) .& (c_true.==0))./b, dims=1)./den)
    c_bar = vec(sum(  ((c_pred.==0) .& (c_true.==1))./b, dims=1)./den)

    Fapprox = Fabc(a_bar,b_bar,c_bar)
    mn  = CheckMissing(Fapprox)

    Ga = gradient(a -> Fabc(a,b_bar,c_bar), a_bar)[1]
    Gb = gradient(b -> Fabc(a_bar,b,c_bar), b_bar)[1]
    Gc = gradient(c -> Fabc(a_bar,b_bar,c), c_bar)[1]

    Ga2= Ga.*Ga
    Gb2= Gb.*Gb
    Gc2= Gc.*Gc

    Ntest=size(c_pred,1)
    vA=zeros(size(c_pred,2))
    vB=zeros(size(c_pred,2))
    vC=zeros(size(c_pred,2))
    lambda=0.9
    c_pred_s =lambda*c_pred .+ (1-lambda)*0.5
    c_true_s =lambda*c_true .+ (1-lambda)*0.5 # this isn't really correct since we need to include additional terms
    c_true_s =1.0*c_true
    for k=1:size(c_pred,2)
        for n=1:size(c_pred,1)
            vA[k] += c_pred_s[n,k]*c_true_s[n,k]*Ga2[k]*(1 ./ b[n] -1)./b[n]
            vB[k] += c_pred_s[n,k]*(1-c_true_s[n,k])*Gb2[k]*(1 ./ b[n] -1)./b[n]
            vC[k] += (1-c_pred_s[n,k])*c_true_s[n,k]*Gc2[k]*(1 ./ b[n] -1)./b[n]
        end
        # Note that there is no term E[DeltaA*DeltaB] since for a given c_true this term is zero
    end
    vA=vA/den^2
    vB=vB/den^2
    vC=vC/den^2

    var = CheckMissing(sum(vA+vB+vC))

    if ismissing(mn) | ismissing(var)
        return missing,missing
    else
        return mn,var
    end

end

function Fabc(a::Vector,b::Vector,c::Vector)
    C=length(a)
    aa=a .+ small_value
    bb=b .+ small_value
    cc=c .+ small_value

    p= sum(aa./(aa + bb))/C
    r= sum(aa./(aa + cc))/C

    return 2*p*r/(p+r)
end

function FapproxFn(c_pred::Vector{Int},p1_true_guess;metric,metric_par=0)
    (ftilde1,ftilde0,gtilde1,gtilde0) = fgtilde(c_pred;metric,metric_par)
    ftilde = p1_true_guess.*ftilde1 + (1 .- p1_true_guess).*ftilde0
    gtilde = p1_true_guess.*gtilde1 + (1 .- p1_true_guess).*gtilde0
    return sum(ftilde)/sum(gtilde)
end

function happroxFn(c_pred::Vector{Int},p1_true_guess;metric,metric_par=0)
    Fapprox = FapproxFn(c_pred,p1_true_guess;metric,metric_par)
    (ftilde1,ftilde0,gtilde1,gtilde0) = fgtilde(c_pred;metric,metric_par)
    htilde2 = p1_true_guess.*(ftilde1 -Fapprox.*gtilde1).^2 + (1 .- p1_true_guess).*(ftilde0 .- Fapprox.*gtilde0).^2
    return sqrt.(htilde2)
end

function MultiMetrics(c_pred,c_true)
    true_positive=sum((c_pred.==1) .& (c_true.==1))
    false_positive=sum((c_pred.==1) .& (c_true.==0))
    false_negative=sum((c_pred.==0) .& (c_true.==1))
    true_negative=sum((c_pred.==0) .& (c_true.==0))
    F1_micro=true_positive/(true_positive + 0.5*(false_positive+false_negative)) # micro F1

    accuracy_individual=(true_positive+true_negative)/(true_positive+true_negative+false_positive+false_negative) # fraction of individually correct entries
    accuracy_joint=mean(all(c_true .== c_pred,dims=2)) # fraction of correct vector predictions

    true_positive=mean((c_pred.==1) .& (c_true.==1),dims=1)
    false_positive=mean((c_pred.==1) .& (c_true.==0),dims=1)
    false_negative=mean((c_pred.==0) .& (c_true.==1),dims=1)
    true_negative=mean((c_pred.==0) .& (c_true.==0),dims=1)
    Prec = mean(true_positive./(true_positive .+ false_positive))
    Rec = mean(true_positive./(true_positive .+ false_negative))

    F1_macro=2*Prec*Rec/(Prec+Rec)

    return F1_micro,F1_macro,accuracy_individual,accuracy_joint
end

function FapproxFn_multi_micro(c_pred,p1_true_guess)
    alpha=0.5
    C=size(c_pred,2)
    num=0.0
    den=0.0
    for n=1:size(c_pred,1)
        for k=1:C
            num=num+ c_pred[n,k]*p1_true_guess[n,k]
            den=den+ alpha*c_pred[n,k] + (1-alpha)*p1_true_guess[n,k]
        end
    end
    return num/den
end

function h_multi_micro(cp,ct,F)
    return f_micro(cp,ct)-F.*g_micro(cp,ct)
end

function happroxFn_multi_micro(c_pred,p1_true_guess)
    F =FapproxFn_multi_micro(c_pred,p1_true_guess)
    p1=1.0*p1_true_guess
    p0=1.0 .- p1
    h2=zeros(size(c_pred,1))
    avh=zeros(size(c_pred,1))
    for n=1:size(c_pred,1)
        tmp1=0
        for k=1:size(c_pred,2)
            tmp1 = tmp1 + p0[n,k]*h_multi_micro(c_pred[n,k],0,F)^2 +
                p1[n,k]*h_multi_micro(c_pred[n,k],1,F)^2
        end
        tmp2=0
        for k=1:size(c_pred,2)
            tmp2 = tmp2 + (p0[n,k]*h_multi_micro(c_pred[n,k],0,F) +
                p1[n,k]*h_multi_micro(c_pred[n,k],1,F))^2
        end
        tmp3=0
        for k=1:size(c_pred,2)
            tmp3 = tmp3 + p0[n,k]*h_multi_micro(c_pred[n,k],0,F) +
                    p1[n,k]*h_multi_micro(c_pred[n,k],1,F)
        end
        h2[n]=tmp1-tmp2+tmp3^2

        for k=1:size(c_pred,2)
            avh[n] += p0[n,k]*h_multi_micro(c_pred[n,k],0,F) +
                p1[n,k]*h_multi_micro(c_pred[n,k],1,F)
        end

    end
    return sqrt.(h2), avh
end

function happroxFn_multi_macro(c_pred,p1_true_guess)

    a_bar = vec(mean( p1_true_guess .* ((c_pred.==1) .& (c_true.==1)), dims=1)) # check!
    b_bar = vec(mean( (1 .- p1_true_guess) .* ((c_pred.==1) .& (c_true.==0)), dims=1)) # check!
    c_bar = vec(mean( p1_true_guess .* ((c_pred.==0) .& (c_true.==1)), dims=1)) # check!

    Ga = gradient(a -> Fabc(a,b_bar,c_bar), a_bar)[1]
    Gb = gradient(b -> Fabc(a_bar,b,c_bar), b_bar)[1]
    Gc = gradient(c -> Fabc(a_bar,b_bar,c), c_bar)[1]

    Ga2= Ga.*Ga
    Gb2= Gb.*Gb
    Gc2= Gc.*Gc

    p1=1.0*p1_true_guess
    p0=1.0 .- p1
    h2=zeros(size(c_pred,1))
    for n=1:size(c_pred,1)
        tmpa=0
        for k=1:size(c_pred,2)
            tmpa = tmpa + c_pred[n,k]*p1[n,k]*Ga2[k]
        end
        tmpb=0
        for k=1:size(c_pred,2)
            tmpb = tmpb + c_pred[n,k]*p0[n,k]*Gb2[k]
        end
        tmpc=0
        for k=1:size(c_pred,2)
            tmpc = tmpc + (1-c_pred[n,k])*p1[n,k]*Gc2[k] #  check!
        end

        h2[n]=tmpa+tmpb+tmpc
    end
    return sqrt.(h2)
end

function CheckMissing(x)
    if isnan(x) | ismissing(x)
        return missing
    else
        return x
    end
end

function optimalb(h,M)
    totalh=sum(h)
    sumind=0
    order=sortperm(-h)
    b=zeros(length(h))
    if M*h[order[1]]./totalh <=1
        return M*h/totalh
    end

    count=0
    sumh=0
    for ind in order
        count += 1
        sumh += h[ind]
        D=totalh - sumh
        N=M-count
        G=h[ind]*N/D
        if G <=1
            b[order[1:count-1]]=ones(count-1)
            b[order[count:end]]=h[order[count:end]]*N/D
            break
        end
    end
   return b
end

function optimalb!(h,M)
    totalh=sum(h)
    sumind=0
    order=sortperm(-h)

    if M*h[order[1]]./totalh <=1
        copy!(b,M*h/totalh)
    end

    count=0
    sumh=0
    for ind in order
        count += 1
        sumh += h[ind]
        D=totalh - sumh
        N=M-count
        G=h[ind]*N/D
        if G <=1
            b[order[1:count-1]]=ones(count-1)
            b[order[count:end]]=h[order[count:end]]*N/D
            break
        end
    end

end

max_stars_repo_name: davidbarber/TestEstimationPublic

max_issues_repo_name: davidbarber/TestEstimationPublic

max_forks_repo_name: davidbarber/TestEstimationPublic


text:
function sbatch(file; kwargs...)
    for (k,v) in kwargs
        k = replace(string(k), "_"=>"-")
        v = string(v)
        kw = if length(k) == 1
            "-$(k) $(v)"
        elseif length(v) == 0
            "--$(k)"
        else
            "--$(k)=$(v)"
        end
        write(file, "#SBATCH $(kw)\n")
    end
end

function setup_sweep(tasks::AbstractVector{String}, name::String;
                     port=12385, log_dir="logs",
                     extra_cmds = [],
                     kwargs...)
    open("tasks", "w") do file
        write(file, join(tasks, "\n"))
    end
    mkpath(log_dir)

    open("slurm.sh", "w") do file
        write(file, "#!/bin/sh\n")
        sbatch(file;
               J=name,
               o="$(log_dir)/$(name)_%j.out",
               e="$(log_dir)/$(name)_%j.out",
               kwargs...)

        write(file, "set -e\n")
        for cmd in extra_cmds
            write(file, "$(cmd)\n")
        end

        write(file, "srun -Q -N \${SLURM_NNODES} -n \${SLURM_NPROCS} julia -e 'using Pkg; Pkg.activate(\".\"); import SlurmSweeps; SlurmSweeps.instance(port=$(port), log_dir=\"$(log_dir)\")'")
    end
end

export setup_sweep


text:
using AbstractGPs
using ConjugateComputationVI
using Distributions
using Optim
using ParameterHandling
using Plots
using Random
using RDatasets
using StatsFuns
using Zygote

using ConjugateComputationVI:
    batch_quadrature,
    GaussHermiteQuadrature,
    UnivariateFactorisedLikelihood

# Prepare data.

# Get the (infamous) Iris dataset, and only keep two labels -- we'll just do binary
# classification for now.
iris = filter(
    row -> row.Species == "setosa" || row.Species == "virginica",
    dataset("datasets", "iris"),
)

setosa = filter(row -> row.Species == "setosa", iris)
virginica = filter(row -> row.Species == "virginica", iris)

scatter(setosa.SepalLength, setosa.SepalWidth; label="Setosa")
scatter!(virginica.SepalLength, virginica.SepalWidth; label="Virginica")

# Check that there are only two features.
@assert length(unique(iris.Species)) == 2

# Encode the specifies using Ints. Could also have used Booleans.
y = map(x -> x == "setosa" ? 1 : 0, iris.Species)

# Construct input matrix from feature columns.
X = hcat(iris.SepalLength, iris.SepalWidth);
x = ColVecs(collect(X'));

# Split into train and test sets.
train_indices = randperm(MersenneTwister(123456), length(y))[1:50];
test_indices = setdiff(eachindex(y), train_indices);
x_tr = x[train_indices];
x_te = x[test_indices];
y_tr = y[train_indices];
y_te = y[test_indices];

# Specify and perform inference in model.

θ_init = (scale=fixed(1.0), stretch=positive.(ones(2)));

θ_init_flat, unflatten = ParameterHandling.flatten(θ_init);

function build_latent_gp(θ::AbstractVector{<:Real})
    return build_latent_gp(ParameterHandling.value(unflatten(θ)))
end
function build_latent_gp(θ::NamedTuple)
    gp = GP(θ.scale * SEKernel() ∘ AbstractGPs.ScaleTransform(θ.stretch))
    lik = UnivariateFactorisedLikelihood(f -> Bernoulli(logistic(f)))
    return LatentGP(gp, lik, 1e-9)
end

f_approx_post, results_summary = ConjugateComputationVI.optimize_elbo(
    build_latent_gp,
    GaussHermiteQuadrature(10),
    x_tr,
    y_tr,
    θ_init_flat,
    BFGS(
        alphaguess = Optim.LineSearches.InitialStatic(scaled=true),
        linesearch = Optim.LineSearches.BackTracking(),
    ),
    Optim.Options(
        show_trace = true,
        iterations=25,
        f_calls_limit=50,
    ),
);

# Visualise results of model.

function post_mean_prediction(x::AbstractVector)
    m, v = mean_and_var(f_approx_post(x))
    conditionals = fill(f -> pdf(Bernoulli(logistic(f)), 1), eachindex(m))
    return batch_quadrature(conditionals, m, sqrt.(v), 25)
end

# Compute predictions using the optimal hyperparameters.
function MAP_predict(x::AbstractVector)
    posterior_mean_predictions = post_mean_prediction(x)
    MAP_predictions = map(x -> x > 0.5 ? 1 : 0, posterior_mean_predictions)
    return MAP_predictions
end

# Compute training set accuracy as a %age.
100 * mean(MAP_predict(x_tr) .== y_tr)

# Compute test set accuracy as a %age.
100 * mean(MAP_predict(x_te) .== y_te)

# Compute prediction surface.
x1_range = range(minimum(X[:, 1]) - 3, maximum(X[:, 1]) + 3; length=100);
x2_range = range(minimum(X[:, 2]) - 3, maximum(X[:, 2]) + 3; length=100);
x_pr_vecs = [[x1, x2] for x1 in x1_range, x2 in x2_range];
X_pr = reduce(hcat, x_pr_vecs);
x_pr = ColVecs(X_pr);

# Plot classification.
MAP_pr = MAP_predict(x_pr);
p1 = heatmap(
    x1_range, x2_range, reshape(MAP_pr, length(x2_range), length(x1_range))';
    label="classifier",
);
scatter!(p1, setosa.SepalLength, setosa.SepalWidth; label="");
scatter!(p1, virginica.SepalLength, virginica.SepalWidth; label="");

# Plot posterior mean prob.
post_mean_pr = post_mean_prediction(x_pr);
p2 = heatmap(
    x1_range, x2_range, reshape(post_mean_pr, length(x2_range), length(x1_range))';
    label="prob",
);
scatter!(p2, setosa.SepalLength, setosa.SepalWidth; label="");
scatter!(p2, virginica.SepalLength, virginica.SepalWidth; label="");

# Plot posterior marginals.
ms = marginals(f_approx_post(x_pr));
p3 = heatmap(
    x1_range, x2_range, reshape(mean.(ms), length(x2_range), length(x1_range))';
    label="latent mean",
);
scatter!(p3, setosa.SepalLength, setosa.SepalWidth; label="");
scatter!(p3, virginica.SepalLength, virginica.SepalWidth; label="");

p4 = heatmap(
    x1_range, x2_range, reshape(std.(ms), length(x2_range), length(x1_range))';
    label="latent std",
);
scatter!(p4, setosa.SepalLength, setosa.SepalWidth; label="");
scatter!(p4, virginica.SepalLength, virginica.SepalWidth; label="");

plot(p1, p2, p3, p4; layout=(2, 2))


text:
"""
```
transform_data(m::AbstractModel, levels::DataFrame; verbose::Symbol = :low)
```

Transform data loaded in levels and order columns appropriately for the DSGE model. Returns
DataFrame of transformed data.

The DataFrame `levels` is output from `load_data_levels`. The series in levels are
transformed as specified in `m.data_transforms`.
- To prepare for per-capita transformations, population data are filtered using
    `hpfilter`. The series in `levels` to use as the population series is given by the
    `population_mnemonic` setting. If `use_population_forecast` is `true`, a population
    forecast is appended to the recorded population levels before the filtering. Both
    filtered and unfiltered population levels and growth rates are added to the `levels`
    data frame.
- The transformations are applied for each series using the `levels` DataFrame as input.
"""
function transform_data(m::AbstractModel, levels::DataFrame; verbose::Symbol = :low)

    population_mnemonic = get_setting(m, :population_mnemonic)
    n_obs, _ = size(levels)

    # Step 1: HP filter population forecasts, if they're being used

    # population_recorded: historical population, unfiltered
    # population_all: full unfiltered series (including forecast)
    # dlpopulation_forecast: growth rates of population forecasts pre-filtering

    population_recorded = levels[:,[:date, population_mnemonic]]
    population_all, dlpopulation_forecast, n_population_forecast_obs = if use_population_forecast(m)
        if VERBOSITY[verbose] >= VERBOSITY[:high]
            println("Loading population forecast...")
        end

        # load population forecast
        population_forecast_file = inpath(m, "data", "population_forecast_$(data_vintage(m)).csv")
        pop_forecast = readtable(population_forecast_file)

        rename!(pop_forecast, :POPULATION,  population_mnemonic)
        DSGE.na2nan!(pop_forecast)
        DSGE.format_dates!(:date, pop_forecast)

        # use our "real" series as current value
        pop_all = [population_recorded; pop_forecast[2:end,:]]

        # return values
        pop_all[population_mnemonic],
        difflog(pop_forecast[population_mnemonic]),
        length(pop_forecast[population_mnemonic])
    else
        population_recorded[:,population_mnemonic], [NaN], 1
    end

    # hp filter
    population_all = convert(Array, population_all)
    filtered_population, _ = hpfilter(population_all, 1600)

    # filtered series (levels)
    filtered_population_recorded = filtered_population[1:end-n_population_forecast_obs+1]
    filtered_population_forecast = filtered_population[end-n_population_forecast_obs+1:end]

    # filtered growth rates
    dlpopulation_recorded          = difflog(population_recorded[population_mnemonic])
    dlfiltered_population_recorded = difflog(filtered_population_recorded)

    levels[:filtered_population]          = filtered_population_recorded
    levels[:filtered_population_growth]   = dlfiltered_population_recorded
    levels[:unfiltered_population_growth] = dlpopulation_recorded

    # Step 2: apply transformations to each series
    transformed = DataFrame()
    transformed[:date] = levels[:date]

    for series in keys(m.data_transforms)
        if VERBOSITY[verbose] >= VERBOSITY[:high]
            println("Transforming series " * string(series) * "...")
        end
        f = m.data_transforms[series]
        transformed[series] = f(levels)
    end

    sort!(transformed, cols = :date)
end


text:
module FieldModule

using Printf: @printf
using UserExceptionModule: UserError
using CoordsModule: Coords,LatLonCoords,DirectionIndicator,Generic1DCoords
using GridModule: Grid,LatLonGrid, for_all, wrap_coords, coords_in_grid,for_all_with_line_breaks
using GridModule: UnstructuredGrid
using InteractiveUtils: subtypes
using SpecialDirectionCodesModule
using MergeTypesModule
import Base.maximum
import Base.*
import Base./
import Base.+
import Base.==
import Base.isapprox
import Base.fill!
import Base.show
import Base.round
import Base.convert
import Base.sum
import Base.Broadcast
using InteractiveUtils

abstract type Field{T} end

function get(field::Field{T},coords::Coords) where {T}
  throw(UserError())
end

function set!(field::Field,coords::Coords)
  throw(UserError())
end

function fill!(field::Field)
  throw(UserError())
end

function maximum(field::Field)
  throw(UserError())
end

function +(lfield::Field,rfield::Field)
  throw(UserError())
end

function *(lfield::Field,value::T) where {T<:Number}
  throw(UserError())
end

function /(lfield::Field,value::T) where {T<:Number}
  throw(UserError())
end

function ==(lfield::Field,rfield::Field)
  throw(UserError())
end

function sum(field::Field)
  throw(UserError())
end

# This could be done by overloading broadcast but this is overcomplicated
# Assumes result is fractional and hence returns a field of floating point numbers
function elementwise_divide(lfield::Field{T},rfield::Field{T}) where {T}
  return Field{Float64}(lfield.grid,lfield.data./rfield.data)::Field{Float64}
end

# This could be done by overloading broadcast but this is overcomplicated
function equals(lfield::Field,value::T) where {T}
  return Field{Bool}(lfield.grid,Bool.(lfield.data .== value))::Field{Bool}
end

function isapprox(lfield::Field,rfield::Field;
                  rtol::Real=atol>0 ? 0 : sqrt(eps),
                  atol::Real,nans::Bool)
  throw(UserError())
end

function invert(field::Field)
    throw(UserError())
end

function add_offset(field::Field,offset::T,skip_offset_on_values::Vector{T}) where {T}
    throw(UserError())
end

function round!(type::DataType,field::Field{T}) where {T}
  throw(UserError())
end

function add_offset(field::T2,offset::T,skip_offset_on_values::Vector{T}) where {T,T2<:Field{T}}
  for_all(get_grid(field)) do coords::Coords
    if ! (field(coords) in skip_offset_on_values)
      set!(field,coords,field(coords)+offset)
    end
  end
end

function divide(field::T2,divisor::T) where {T,T2<:Field{T}}
  for_all(get_grid(field)) do coords::Coords
    set!(field,coords,field(coords)/divisor)
  end
end

get_grid(obj::T) where {T <: Field} =
  obj.grid::Grid

get_data(obj::T) where {T <: Field} =
  obj.data::Array

function repeat(field::Field{T}, count::Integer) where {T}
  array::Array{Field{T},1} = Array{Field{T},1}(undef,count)
  for i in 1:count
    array[i] = deepcopy(field)
  end
  return array
end

struct LatLonField{T} <: Field{T}
  data::Array{T,2}
  grid::LatLonGrid
  function LatLonField{T}(grid::LatLonGrid,value::T) where {T}
    data::Array{T,2} = value == zero(T) ? zeros(T,grid.nlat,grid.nlon) :
                                          ones(T,grid.nlat,grid.nlon)*value
    return new(data,grid)
  end
  function LatLonField{MergeTypes}(grid::LatLonGrid,value::MergeTypes)
    data::Array{MergeTypes,2} =  zeros(MergeTypes,grid.nlat,grid.nlon)
    fill!(data,value)
    return new(data,grid)
  end
  function LatLonField{T}(grid::LatLonGrid,values::AbstractArray{T,2}) where {T}
    if size(values,1) != grid.nlat ||
       size(values,2) != grid.nlon
       error("Values provided don't match selected grid")
    end
    return new(values,grid)
  end
  function LatLonField{MergeTypes}(grid::LatLonGrid,
                                   integer_field::LatLonField{T}) where {T<:Signed}
    grid = get_grid(integer_field)
    merge_type_field = LatLonField{MergeTypes}(grid,null_mtype)
    for_all(grid) do coords::Coords
      element::MergeTypes = MergeTypes(get(integer_field,coords))
      set!(merge_type_field,coords,element)
    end
    return new(merge_type_field.data,grid)
  end
end

struct UnstructuredField{T} <: Field{T}
  data::Array{T,1}
  grid::UnstructuredGrid
  function UnstructuredField{T}(grid::UnstructuredGrid,value::T) where {T}
    data::Array{T,1} = value == zero(T) ? zeros(T,grid.ncells) :
                                          ones(T,grid.ncells)*value
    return new(data,grid)
  end
  function UnstructuredField{MergeTypes}(grid::UnstructuredGrid,value::MergeTypes)
    data::Array{MergeTypes,1} =  zeros(MergeTypes,grid.ncells)
    fill!(data,value)
    return new(data,grid)
  end
  function UnstructuredField{T}(grid::UnstructuredGrid,values::AbstractArray{T,1}) where {T}
    if size(values,1) != grid.ncells
       error("Values provided don't match selected grid")
    end
    return new(values,grid)
  end
  function UnstructuredField{MergeTypes}(grid::UnstructuredGrid,
                                   integer_field::UnstructuredField{T}) where {T<:Signed}
    grid = get_grid(integer_field)
    merge_type_field = UnstructuredField{MergeTypes}(grid,null_mtype)
    for_all(grid) do coords::Coords
      element::MergeTypes = MergeTypes(get(integer_field,coords))
      set!(merge_type_field,coords,element)
    end
    return new(merge_type_field.data,grid)
  end
end

LatLonField{T}(grid::LatLonGrid) where {T} = LatLonField{T}(grid,zeros(T,grid.nlat,grid.nlon))

UnstructuredField{T}(grid::UnstructuredGrid) where {T} = UnstructuredField{T}(grid,zeros(T,grid.ncells))

Field{T}(grid::LatLonGrid,value::T) where {T} = LatLonField{T}(grid::LatLonGrid,value::T)
Field{T}(grid::LatLonGrid) where {T} = LatLonField{T}(grid::LatLonGrid)

Field{T}(grid::UnstructuredGrid,value::T) where {T} = UnstructuredField{T}(grid::UnstructuredGrid,
                                                                           value::T)
Field{T}(grid::UnstructuredGrid) where {T} = UnstructuredField{T}(grid::UnstructuredGrid)

Field{T}(grid::LatLonGrid,values::AbstractArray{T,2}) where {T} =
  LatLonField{T}(grid::LatLonGrid,values::AbstractArray{T,2})
Field{T}(grid::UnstructuredGrid,values::AbstractArray{T,1}) where {T} =
  UnstructuredField{T}(grid::UnstructuredGrid,
                       values::AbstractArray{T,1})

Field{MergeTypes}(grid::LatLonGrid,integer_field::LatLonField{T}) where {T<:Signed} =
  LatLonField{MergeTypes}(grid::LatLonGrid,integer_field::LatLonField{T})
Field{MergeTypes}(grid::UnstructuredGrid,integer_field::UnstructuredField{T}) where {T<:Signed} =
  UnstructuredField{MergeTypes}(grid::UnstructuredGrid,integer_field::UnstructuredField{T})

convert(::Type{Field{MergeTypes}},x::T2) where {T<:Signed,T2<:Field{T}} = Field{MergeTypes}(get_grid(x),x)

function (latlon_field::LatLonField{T})(coords::LatLonCoords) where {T}
    return latlon_field.data[coords.lat,coords.lon]::T
end

function get(latlon_field::LatLonField{T},coords::LatLonCoords) where {T}
    return latlon_field(coords)::T
end

function set!(latlon_field::LatLonField{T},coords::LatLonCoords,value::T) where {T}
    latlon_field.data[coords.lat,coords.lon] = value
end

function +(lfield::LatLonField{T},rfield::LatLonField{T}) where {T}
  return LatLonField{T}(lfield.grid,lfield.data + rfield.data)
end

function *(lfield::Field,value::T) where {T<:Number}
  return LatLonField{T}(lfield.grid,lfield.data * value)
end

function /(lfield::Field,value::T) where {T<:Number}
  return LatLonField{T}(lfield.grid,lfield.data / value)
end

function ==(lfield::LatLonField{T},rfield::LatLonField{T}) where {T}
  return (lfield.data == rfield.data)::Bool
end

LatLonFieldOrUnstructuredField = Union{LatLonField{T},UnstructuredField{T}} where {T}

function sum(field::T) where {T<:LatLonFieldOrUnstructuredField}
  return sum(field.data)
end

function isapprox(lfield::LatLonField{T},rfield::LatLonField{T};
                  rtol::Real=atol>0 ? 0 : sqrt(eps),
                  atol::Real=0.0,nans::Bool=false) where {T}
  return isapprox(lfield.data,rfield.data,
                 rtol=rtol,atol=atol,nans=nans)::Bool
end

function invert(field::LatLonField{T}) where {T}
  return LatLonField{T}(field.grid,.!field.data)
end

function round(type::DataType,field::LatLonField{T}) where {T}
  return LatLonField{type}(field.grid,round.(type,field.data))
end

function maximum(field::LatLonField)
  return maximum(field.data)
end

function show(io::IO,field::LatLonField{T}) where {T <: Number}
  for_all_with_line_breaks(field.grid) do coords::Coords
    if T <: AbstractFloat
      @printf(io,"%6.2f",field(coords))
    elseif T == Bool
      print(io,field(coords) ? " X " : " - ")
    else
      @printf(io,"%4.0f",field(coords))
    end
    print(io," ")
  end
end

function fill!(field::LatLonField{T},value::T) where {T}
  return fill!(field.data,value)
end

function (unstructured_field::UnstructuredField{T})(coords::Generic1DCoords) where {T}
    return unstructured_field.data[coords.index]::T
end

function get(unstructured_field::UnstructuredField{T},coords::Generic1DCoords) where {T}
    return unstructured_field(coords)::T
end

function set!(unstructured_field::UnstructuredField{T},coords::Generic1DCoords,value::T) where {T}
    unstructured_field.data[coords.index] = value
end

function +(lfield::UnstructuredField{T},rfield::UnstructuredField{T}) where {T}
  return UnstructuredField{T}(lfield.grid,lfield.data + rfield.data)
end

function *(lfield::UnstructuredField,value::T) where {T<:Number}
  return UnstructuredField{T}(lfield.grid,lfield.data * value)
end

function /(lfield::UnstructuredField,value::T) where {T<:Number}
  return UnstructuredField{T}(lfield.grid,lfield.data / value)
end

function ==(lfield::UnstructuredField{T},rfield::UnstructuredField{T}) where {T}
  return (lfield.data == rfield.data)::Bool
end

function isapprox(lfield::UnstructuredField{T},rfield::UnstructuredField{T};
                  rtol::Real=atol>0 ? 0 : sqrt(eps),
                  atol::Real=0.0,nans::Bool=false) where {T}
  return isapprox(lfield.data,rfield.data,
                 rtol=rtol,atol=atol,nans=nans)::Bool
end

function invert(field::UnstructuredField{T}) where {T}
  return UnstructuredField{T}(field.grid,.!field.data)
end

function round(type::DataType,field::UnstructuredField{T}) where {T}
  return UnstructuredField{type}(field.grid,round.(type,field.data))
end

function maximum(field::UnstructuredField)
  return maximum(field.data)
end

function fill!(field::UnstructuredField{T},value::T) where {T}
  return fill!(field.data,value)
end

function show(io::IO,field::UnstructuredField{T}) where {T <: Number}
  for_all_with_line_breaks(field.grid) do coords::Coords
    if T <: AbstractFloat
      @printf(io,"%6.2f",field(coords))
    elseif T == Bool
      print(io,field(coords) ? " X " : " - ")
    else
      @printf(io,"%4.0f",field(coords))
    end
    print(io," ")
  end
end

abstract type DirectionIndicators end

function get(field::DirectionIndicators,coords::Coords)
  throw(UserError())
end

function set!(field::DirectionIndicators,coords::Coords)
  throw(UserError())
end

struct LatLonDirectionIndicators <: DirectionIndicators
  next_cell_coords::Array{LatLonField{Int64},1}
  function LatLonDirectionIndicators(dir_based_rdirs::LatLonField{Int64})
    grid = dir_based_rdirs.grid
    nlat = dir_based_rdirs.grid.nlat
    nlon = dir_based_rdirs.grid.nlon
    lat_indices::LatLonField{Int64} =
      LatLonField{Int64}(grid,zeros(Int64,nlat,nlon))
    lon_indices::LatLonField{Int64} =
      LatLonField{Int64}(grid,zeros(Int64,nlat,nlon))
    for_all(grid) do coords::LatLonCoords
      dir_based_rdir::Int64 = dir_based_rdirs(coords)
      if dir_based_rdir == dir_based_indicator_ocean
        set!(lat_indices,coords,coord_base_indicator_ocean)
        set!(lon_indices,coords,coord_base_indicator_ocean)
      elseif dir_based_rdir == dir_based_indicator_outflow
        set!(lat_indices,coords,coord_base_indicator_outflow)
        set!(lon_indices,coords,coord_base_indicator_outflow)
      elseif dir_based_rdir == dir_based_indicator_truesink
        set!(lat_indices,coords,coord_base_indicator_truesink)
        set!(lon_indices,coords,coord_base_indicator_truesink)
      elseif dir_based_rdir == dir_based_indicator_lake
        set!(lat_indices,coords,coord_base_indicator_lake)
        set!(lon_indices,coords,coord_base_indicator_lake)
      elseif 1 <= dir_based_rdir <= 9
        lat_offset::Int64 = dir_based_rdir in [7,8,9] ? -1 :
                           (dir_based_rdir in [1,2,3] ?  1 : 0)
        lon_offset::Int64 = dir_based_rdir in [7,4,1] ? -1 :
                           (dir_based_rdir in [9,6,3] ?  1 : 0)
        destination_coords::LatLonCoords = LatLonCoords(coords.lat + lat_offset,
                                                        coords.lon + lon_offset)
        destination_coords = wrap_coords(grid,destination_coords)
        if ! coords_in_grid(grid,destination_coords)
          error("Direction based direction indicator points out of grid")
        end
        set!(lat_indices,coords,destination_coords.lat)
        set!(lon_indices,coords,destination_coords.lon)
      else
        error("Invalid direction based direction indicator")
      end
    end
    new(LatLonField{Int64}[lat_indices,lon_indices])
  end
end

struct UnstructuredFieldDirectionIndicators <: DirectionIndicators
  next_cell_coords::UnstructuredField{Int64}
end

function (latlon_direction_indicators::LatLonDirectionIndicators)(coords::LatLonCoords)
  return DirectionIndicator(LatLonCoords(latlon_direction_indicators.next_cell_coords[1](coords),
                                         latlon_direction_indicators.next_cell_coords[2](coords)))
end

function get(latlon_direction_indicators::LatLonDirectionIndicators,coords::LatLonCoords)
  return latlon_direction_indicators(coords)
end

function set!(latlon_direction_indicators::DirectionIndicators,
               coords::LatLonCoords,direction::DirectionIndicator)
  next_cell_coords =  get_next_cell_coords(direction,coords)
  set!(latlon_direction_indicators.next_cell_coords[1],coords,next_cell_coords.lat)
  set!(latlon_direction_indicators.next_cell_coords[2],coords,next_cell_coords.lon)
end

function fill!(field::LatLonField{T},value::T) where {T}
  return fill!(field.data,value)
end

struct UnstructuredDirectionIndicators <: DirectionIndicators
  next_cell_coords::UnstructuredField{Int64}
end

function (unstructured_direction_indicators::UnstructuredDirectionIndicators)(coords::Generic1DCoords)
  return DirectionIndicator(Generic1DCoords(unstructured_direction_indicators.next_cell_coords(coords)))
end

function get(unstructured_direction_indicators::UnstructuredDirectionIndicators,
             coords::Generic1DCoords)
  return unstructured_direction_indicators(coords)
end

function set!(unstructured_direction_indicators::UnstructuredDirectionIndicators,
              coords::Generic1DCoords,direction::DirectionIndicator)
  next_cell_coords =  get_next_cell_coords(direction,coords)
  set!(unstructured_direction_indicators.next_cell_coords,coords,next_cell_coords.index)
end

#Cannot define functors on abstract types in Julia yet. Use this
#workaround - needs to be at end of file as it only acts on subtypes of
#Field/DirectionIndicators that are already defined

for T in subtypes(Field)
 @eval function (field::$T)(coords::Coords) throw(UserError()) end
end

for T in subtypes(DirectionIndicators)
  @eval function (direction_indicators::$T)(coords::Coords) throw(UserError()) end
end

end

max_stars_repo_name: ThomasRiddick/DynamicHD

max_issues_repo_name: ThomasRiddick/DynamicHD

max_forks_repo_name: ThomasRiddick/DynamicHD


text:
import GPLinearODEMaker.powers_of_negative_one

"""
cos_kernel function created by kernel_coder(). Requires 1 hyperparameters. Likely created using cos_kernel_base() as an input.
Use with include("src/kernels/cos_kernel.jl").
hyperparameters == ["λ"]
"""
function cos_kernel(
    hyperparameters::Vector{<:Real},
    δ::Real,
    dorder::Vector{<:Integer};
    shift_ind::Integer=0)

    include_shift = shift_ind!=0
    @assert length(hyperparameters)==1+Int(include_shift) "hyperparameters is the wrong length"
    dorder_len = 3 + Int(include_shift)
    @assert length(dorder)==dorder_len "dorder is the wrong length"
    @assert maximum(dorder) < 3 "No more than two derivatives for each time or hyperparameter can be calculated"

    dorder2 = dorder[2]
    dorder[2] += dorder[1]

    if include_shift
        dorder[2] += dorder[shift_ind+2]
        δ += hyperparameters[shift_ind]
        hyperparameters_view = view(hyperparameters, 1:2 .!= shift_ind)
        dorder_view = view(dorder, 2:dorder_len)
        dorder_view = view(dorder_view, 1:(dorder_len-1) .!= (shift_ind+1))
    else
        hyperparameters_view = hyperparameters
        dorder_view = view(dorder, 2:dorder_len)
    end

    λ = hyperparameters_view[1]

    if dorder_view==[6, 2]
        func = -2584214.15233042*cos(6.28318530717959*δ/λ)/λ^8 + 5412365.46417601*δ*sin(6.28318530717959*δ/λ)/λ^9 + 2429063.94011407*δ^2*cos(6.28318530717959*δ/λ)/λ^10
    end

    if dorder_view==[5, 2]
        func = -293778.89739387*sin(6.28318530717959*δ/λ)/λ^7 - 738346.900665833*δ*cos(6.28318530717959*δ/λ)/λ^8 + 386597.533155429*δ^2*sin(6.28318530717959*δ/λ)/λ^9
    end

    if dorder_view==[4, 2]
        func = 31170.9091308808*cos(6.28318530717959*δ/λ)/λ^6 - 97926.29913129*δ*sin(6.28318530717959*δ/λ)/λ^7 - 61528.9083888195*δ^2*cos(6.28318530717959*δ/λ)/λ^8
    end

    if dorder_view==[3, 2]
        func = 2976.60256130878*sin(6.28318530717959*δ/λ)/λ^5 + 12468.3636523523*δ*cos(6.28318530717959*δ/λ)/λ^6 - 9792.629913129*δ^2*sin(6.28318530717959*δ/λ)/λ^7
    end

    if dorder_view==[2, 2]
        func = -236.870505626145*cos(6.28318530717959*δ/λ)/λ^4 + 1488.30128065439*δ*sin(6.28318530717959*δ/λ)/λ^5 + 1558.54545654404*δ^2*cos(6.28318530717959*δ/λ)/λ^6
    end

    if dorder_view==[1, 2]
        func = -12.5663706143592*sin(6.28318530717959*δ/λ)/λ^3 - 157.91367041743*δ*cos(6.28318530717959*δ/λ)/λ^4 + 248.050213442399*δ^2*sin(6.28318530717959*δ/λ)/λ^5
    end

    if dorder_view==[0, 2]
        func = -12.5663706143592*δ*sin(6.28318530717959*δ/λ)/λ^3 - 39.4784176043574*δ^2*cos(6.28318530717959*δ/λ)/λ^4
    end

    if dorder_view==[6, 1]
        func = 369173.450332917*cos(6.28318530717959*δ/λ)/λ^7 - 386597.533155429*δ*sin(6.28318530717959*δ/λ)/λ^8
    end

    if dorder_view==[5, 1]
        func = 48963.149565645*sin(6.28318530717959*δ/λ)/λ^6 + 61528.9083888195*δ*cos(6.28318530717959*δ/λ)/λ^7
    end

    if dorder_view==[4, 1]
        func = -6234.18182617615*cos(6.28318530717959*δ/λ)/λ^5 + 9792.629913129*δ*sin(6.28318530717959*δ/λ)/λ^6
    end

    if dorder_view==[3, 1]
        func = -744.150640327196*sin(6.28318530717959*δ/λ)/λ^4 - 1558.54545654404*δ*cos(6.28318530717959*δ/λ)/λ^5
    end

    if dorder_view==[2, 1]
        func = 78.9568352087149*cos(6.28318530717959*δ/λ)/λ^3 - 248.050213442399*δ*sin(6.28318530717959*δ/λ)/λ^4
    end

    if dorder_view==[1, 1]
        func = 6.28318530717959*sin(6.28318530717959*δ/λ)/λ^2 + 39.4784176043574*δ*cos(6.28318530717959*δ/λ)/λ^3
    end

    if dorder_view==[0, 1]
        func = 6.28318530717959*δ*sin(6.28318530717959*δ/λ)/λ^2
    end

    if dorder_view==[6, 0]
        func = -61528.9083888195*cos(6.28318530717959*δ/λ)/λ^6
    end

    if dorder_view==[5, 0]
        func = -9792.629913129*sin(6.28318530717959*δ/λ)/λ^5
    end

    if dorder_view==[4, 0]
        func = 1558.54545654404*cos(6.28318530717959*δ/λ)/λ^4
    end

    if dorder_view==[3, 0]
        func = 248.050213442399*sin(6.28318530717959*δ/λ)/λ^3
    end

    if dorder_view==[2, 0]
        func = -39.4784176043574*cos(6.28318530717959*δ/λ)/λ^2
    end

    if dorder_view==[1, 0]
        func = -6.28318530717959*sin(6.28318530717959*δ/λ)/λ
    end

    if dorder_view==[0, 0]
        func = cos(6.28318530717959*δ/λ)
    end

    dorder[2] = dorder2  # resetting dorder[2]
    return powers_of_negative_one(dorder2) * float(func)  # correcting for amount of t2 derivatives

end

return cos_kernel, 1  # the function handle and the number of kernel hyperparameters


text:
using Test
using LogicCircuits
using ProbabilisticCircuits
using DataStructures

@testset "Sampling Test" begin
    EPS = 1e-2;
    prob_circuit = zoo_psdd("little_4var.psdd");
    flow_circuit = FlowΔ(prob_circuit, 16, Bool);

    N = 4;
    data_all = XData(generate_data_all(N));

    calc_prob_all = log_likelihood_per_instance(flow_circuit, data_all);
    calc_prob_all = exp.(calc_prob_all);

    using DataStructures
    hist = DefaultDict{AbstractString,Float64}(0.0)

    Nsamples = 1000 * 1000
    for i = 1:Nsamples
        cur = join(Int.(sample(prob_circuit)))
        hist[cur] += 1
    end

    for k in keys(hist)
        hist[k] /= Nsamples
    end

    for k in keys(hist)
        cur = parse(Int32, k, base=2) + 1 # cause Julia arrays start at 1 :(
        @test calc_prob_all[cur] ≈ hist[k] atol= EPS;
    end

end

@testset "Sampling With Evidence" begin
    # TODO (pashak) this test should be improved by adding few more cases
    EPS = 1e-3;
    prob_circuit = zoo_psdd("little_4var.psdd");

    opts= (compact⋀=false, compact⋁=false)
    flow_circuit = UpFlowΔ(prob_circuit, 1, Float64, opts);

    N = 4;
    data = XData(Int8.([0 -1 0 -1]));
    calc_prob = marginal_log_likelihood_per_instance(flow_circuit, data);
    calc_prob = exp.(calc_prob);

    flow_circuit_all = UpFlowΔ(prob_circuit, 4, Float64, opts);
    data_all = XData(Int8.([
                            0 0 0 0;
                            0 0 0 1;
                            0 1 0 0;
                            0 1 0 1;
                        ]));
    calc_prob_all = marginal_log_likelihood_per_instance(flow_circuit_all, data_all);
    calc_prob_all = exp.(calc_prob_all);

    calc_prob_all ./= calc_prob[1]

    hist = DefaultDict{AbstractString,Float64}(0.0)

    Nsamples = 1000 * 1000
    for i = 1:Nsamples
        cur = join(Int.(sample(flow_circuit)))
        hist[cur] += 1
    end

    for k in keys(hist)
        hist[k] /= Nsamples
    end

    for ind = 1:4
        cur = join(data_all.x[ind, :])
        @test calc_prob_all[ind] ≈ hist[cur] atol= EPS;
    end
end


text:
# Include the `Benchmarker` code.
#
# Push through the `test_vgg` to make sure that whole pipeline works.
include("../../experiments/Benchmarker/src/Benchmarker.jl")
using .Benchmarker: Benchmarker

@testset "Testing Benchmarker" begin
    Benchmarker.kernel_profile(Benchmarker.test_vgg())
    Benchmarker.run_conventional(
        Benchmarker.test_vgg(),
        [AutoTM.Optimizer.Static, AutoTM.Optimizer.Synchronous, AutoTM.Optimizer.Numa],
        Benchmarker.common_ratios(),
    )

    # Generate Plots
    Benchmarker.plot_speedup(
        models = [Benchmarker.test_vgg()],
    )

    Benchmarker.plot_conventional_error(
        models = [Benchmarker.test_vgg()],
    )

    Benchmarker.plot_costs(
        pairs = [Benchmarker.test_vgg() => "synchronous"],
    )
end

max_stars_repo_name: darchr/AutoTM

max_issues_repo_name: darchr/AutoTM

max_forks_repo_name: darchr/AutoTM


text:
abstract type AbstractNoise end

struct StandardGaussian{T} <: AbstractNoise
    size::T
end

# Must include RNG support here
noise(s::StandardGaussian{<:Integer}, t) = randn(s.size)
noise(s::StandardGaussian{<:Tuple}, t) = randn(s.size...)

struct DefinedNoise{T} <: AbstractNoise
    values::T
end

noise(s::DefinedNoise{<:Vector}, t) = s.values[t]


text:
@testset "laplace.jl" begin
    # affine transformations
    test_affine_transformations(Laplace, randn(), randn()^2)
end


text:
"""
    StatisticalModel

Abstract supertype for all statistical models.
"""
abstract type StatisticalModel end

"""
    coef(model::StatisticalModel)

Return the coefficients of the model.
"""
function coef end

"""
    coefnames(model::StatisticalModel)

Return the names of the coefficients.
"""
function coefnames end

"""
    coeftable(model::StatisticalModel; level::Real=0.95)

Return a table with coefficients and related statistics of the model.
`level` determines the level for confidence intervals (by default, 95%).

The returned `CoefTable` object implements the
[Tables.jl](https://github.com/JuliaData/Tables.jl/) interface, and can be
converted e.g. to a `DataFrame` via `using DataFrames; DataFrame(coeftable(model))`.
"""
function coeftable end

"""
    confint(model::StatisticalModel; level::Real=0.95)

Compute confidence intervals for coefficients, with confidence level `level` (by default 95%).
"""
function confint end

"""
    deviance(model::StatisticalModel)

Return the deviance of the model relative to a reference, which is usually when applicable
the saturated model. It is equal, *up to a constant*, to ``-2 \\log L``, with ``L``
the likelihood of the model.
"""
function deviance end

"""
    islinear(model::StatisticalModel)

Indicate whether the model is linear.
"""
function islinear end

"""
    nulldeviance(model::StatisticalModel)

Return the deviance of the null model, that is the one including only the intercept.
"""
function nulldeviance end

"""
    loglikelihood(model::StatisticalModel)
    loglikelihood(model::StatisticalModel, observation)

Return the log-likelihood of the model.

With an `observation` argument, return the contribution of `observation` to the
log-likelihood of `model`.

If `observation` is a `Colon`, return a vector of each observation's contribution
to the log-likelihood of the model. In other words, this is the vector of the
pointwise log-likelihood contributions.

In general, `sum(loglikehood(model, :)) == loglikelihood(model)`.
"""
function loglikelihood end

"""
    nullloglikelihood(model::StatisticalModel)

Return the log-likelihood of the null model corresponding to `model`.
This is usually the model containing only the intercept.
"""
function nullloglikelihood end

"""
    score(model::StatisticalModel)

Return the score of the model, that is the gradient of the
log-likelihood with respect to the coefficients.
"""
function score end

"""
    nobs(model::StatisticalModel)

Return the number of independent observations on which the model was fitted. Be careful
when using this information, as the definition of an independent observation may vary
depending on the model, on the format used to pass the data, on the sampling plan
(if specified), etc.
"""
function nobs end

"""
    dof(model::StatisticalModel)

Return the number of degrees of freedom consumed in the model, including
when applicable the intercept and the distribution's dispersion parameter.
"""
function dof end

"""
    mss(model::StatisticalModel)

Return the model sum of squares.
"""
function mss end

"""
    rss(model::StatisticalModel)

Return the residual sum of squares of the model.
"""
function rss end

"""
    informationmatrix(model::StatisticalModel; expected::Bool = true)

Return the information matrix of the model. By default the Fisher information matrix
is returned, while the observed information matrix can be requested with `expected = false`.
"""
function informationmatrix end

"""
    stderror(model::StatisticalModel)

Return the standard errors for the coefficients of the model.
"""
function stderror end

"""
    vcov(model::StatisticalModel)

Return the variance-covariance matrix for the coefficients of the model.
"""
function vcov end

"""
    weights(model::StatisticalModel)

Return the weights used in the model.
"""
function weights end

"""
    isfitted(model::StatisticalModel)

Indicate whether the model has been fitted.
"""
function isfitted end

"""
Fit a statistical model.
"""
function fit end

"""
Fit a statistical model in-place.
"""
function fit! end

"""
    aic(model::StatisticalModel)

Akaike's Information Criterion, defined as ``-2 \\log L + 2k``, with ``L`` the likelihood
of the model, and `k` its number of consumed degrees of freedom
(as returned by [`dof`](@ref)).
"""
aic(model::StatisticalModel) = -2loglikelihood(model) + 2dof(model)

"""
    aicc(model::StatisticalModel)

Corrected Akaike's Information Criterion for small sample sizes (Hurvich and Tsai 1989),
defined as ``-2 \\log L + 2k + 2k(k-1)/(n-k-1)``, with ``L`` the likelihood of the model,
``k`` its number of consumed degrees of freedom (as returned by [`dof`](@ref)),
and ``n`` the number of observations (as returned by [`nobs`](@ref)).
"""
function aicc(model::StatisticalModel)
    k = dof(model)
    n = nobs(model)
    -2loglikelihood(model) + 2k + 2k*(k+1)/(n-k-1)
end

"""
    bic(model::StatisticalModel)

Bayesian Information Criterion, defined as ``-2 \\log L + k \\log n``, with ``L``
the likelihood of the model,  ``k`` its number of consumed degrees of freedom
(as returned by [`dof`](@ref)), and ``n`` the number of observations
(as returned by [`nobs`](@ref)).
"""
bic(model::StatisticalModel) = -2loglikelihood(model) + dof(model)*log(nobs(model))

"""
    r2(model::StatisticalModel)
    r²(model::StatisticalModel)

Coefficient of determination (R-squared).

For a linear model, the R² is defined as ``ESS/TSS``, with ``ESS`` the explained sum of squares
and ``TSS`` the total sum of squares.
"""
function r2 end

"""
    r2(model::StatisticalModel, variant::Symbol)
    r²(model::StatisticalModel, variant::Symbol)

Pseudo-coefficient of determination (pseudo R-squared).

For nonlinear models, one of several pseudo R² definitions must be chosen via `variant`.
Supported variants are:
- `:MacFadden` (a.k.a. likelihood ratio index), defined as ``1 - \\log (L)/\\log (L_0)``;
- `:CoxSnell`, defined as ``1 - (L_0/L)^{2/n}``;
- `:Nagelkerke`, defined as ``(1 - (L_0/L)^{2/n})/(1 - L_0^{2/n})``.
- `:devianceratio`, defined as ``1 - D/D_0``.

In the above formulas, ``L`` is the likelihood of the model,
``L_0`` is the likelihood of the null model (the model with only an intercept),
``D`` is the deviance of the model (from the saturated model),
``D_0`` is the deviance of the null model,
``n`` is the number of observations (given by [`nobs`](@ref)).

The Cox-Snell and the deviance ratio variants both match the classical definition of R²
for linear models.
"""
function r2(model::StatisticalModel, variant::Symbol)
    loglikbased = (:McFadden, :CoxSnell, :Nagelkerke)
    if variant in loglikbased
        ll = loglikelihood(model)
        ll0 = nullloglikelihood(model)
        if variant == :McFadden
            1 - ll/ll0
        elseif variant == :CoxSnell
            1 - exp(2 * (ll0 - ll) / nobs(model))
        elseif variant == :Nagelkerke
            (1 - exp(2 * (ll0 - ll) / nobs(model))) / (1 - exp(2 * ll0 / nobs(model)))
        end
    elseif variant == :devianceratio
        dev  = deviance(model)
        dev0 = nulldeviance(model)
        1 - dev/dev0
    else
        error("variant must be one of $(join(loglikbased, ", ")) or :devianceratio")
    end
end

const r² = r2

"""
    adjr2(model::StatisticalModel)
    adjr²(model::StatisticalModel)

Adjusted coefficient of determination (adjusted R-squared).

For linear models, the adjusted R² is defined as ``1 - (1 - (1-R^2)(n-1)/(n-p))``, with ``R^2``
the coefficient of determination, ``n`` the number of observations, and ``p`` the number of
coefficients (including the intercept). This definition is generally known as the Wherry Formula I.
"""
function adjr2 end

"""
    adjr2(model::StatisticalModel, variant::Symbol)
    adjr²(model::StatisticalModel, variant::Symbol)

Adjusted pseudo-coefficient of determination (adjusted pseudo R-squared).
For nonlinear models, one of the several pseudo R² definitions must be chosen via `variant`.
The only currently supported variants are `:MacFadden`, defined as ``1 - (\\log (L) - k)/\\log (L0)`` and
`:devianceratio`, defined as ``1 - (D/(n-k))/(D_0/(n-1))``.
In these formulas, ``L`` is the likelihood of the model, ``L0`` that of the null model
(the model including only the intercept), ``D`` is the deviance of the model,
``D_0`` is the deviance of the null model, ``n`` is the number of observations (given by [`nobs`](@ref)) and
``k`` is the number of consumed degrees of freedom of the model (as returned by [`dof`](@ref)).
"""
function adjr2(model::StatisticalModel, variant::Symbol)
    k = dof(model)
    if variant == :McFadden
        ll = loglikelihood(model)
        ll0 = nullloglikelihood(model)
        1 - (ll - k)/ll0
    elseif variant == :devianceratio
        n = nobs(model)
        dev  = deviance(model)
        dev0 = nulldeviance(model)
        1 - (dev*(n-1))/(dev0*(n-k))
    else
        error("variant must be one of :McFadden or :devianceratio")
    end
end

const adjr² = adjr2

max_stars_repo_name: JuliaStats/StatsAPI

max_issues_repo_name: JuliaStats/StatsAPI

max_forks_repo_name: JuliaStats/StatsAPI