text:
import CVRPSEP_jll
const LIB_CVRPSEP = CVRPSEP_jll.libcvrpsep

# /* Definition of constants for identification of constraint types. */
@enum ConstraintType begin
    CMGR_CT_MIN_ROUTES        = 101
    CMGR_CT_NODE_DEGREE       = 102  # /* Degree = 2 for each customer. */
    CMGR_CT_CAP               = 103  # /* Capacity constraint. */
    CMGR_CT_GENCAP            = 104  # /* Generalized capacity constraint. */
    # CMGR_CT_FCI               = 104  # /* For public version. */ # duplicate
    CMGR_CT_TWOMATCHING       = 105
    CMGR_CT_COMB              = 106
    CMGR_CT_STR_COMB          = 107  # /* Strengthened comb. */
    CMGR_CT_HYPOTOUR          = 108
    CMGR_CT_EXT_HYPOTOUR      = 109
    CMGR_CT_MSTAR             = 110  # /* Homogeneous multistar. */
    CMGR_CT_WMSTAR            = 111  # /* Weak multistar. */
    CMGR_CT_DJCUT             = 112  # /* Disjunctive cut. */
    CMGR_CT_GOMORY            = 113  # /* By variable numbers */
    CMGR_CT_TWOEDGES_HYPOTOUR = 114  # /* 2EH inequality */
    CMGR_BT_CLIQUE_DOWN       = 201  # /* x(S:S) <= RHS */
    CMGR_BT_CLIQUE_UP         = 202  # /* x(S:S) >= RHS */
    CMGR_BT_STAR_DOWN         = 301  # /* x(i:F) <=RHS */
    CMGR_BT_STAR_UP           = 302  # /* x(i:F) >=RHS */
    CMGR_CT_SLB               = 401  # /* x(F) >= RHS. Simple lower bound */
end

mutable struct CnstrRecord
    Ctype       :: Cint
    Key         :: Cint
    IntListSize :: Cint
    IntList     :: Ptr{Cint}
    ExtListSize :: Cint
    ExtList     :: Ptr{Cint}
    CListSize   :: Cint
    CList       :: Ptr{Cint}
    CoeffList   :: Ptr{Cdouble}
    A           :: Cint
    B           :: Cint
    L           :: Cint
    RHS         :: Cdouble
    BranchLevel :: Cint
    GlobalNr    :: Cint
end

CnstrPointer = Ptr{CnstrRecord}
CnstrPointerList = Ptr{CnstrPointer}

mutable struct CnstrMgrRecord
    CPL         :: CnstrPointerList
    Dim         :: Cint
    Size        :: Cint
end

CnstrMgrPointer = Ptr{CnstrMgrRecord}

# test code

Dim = 100
old_cmp_ref = Ref{CnstrMgrPointer}()
my_cmp_ref = Ref{CnstrMgrPointer}()

ccall(
    (:CMGR_CreateCMgr, LIB_CVRPSEP),
    Cvoid,
    (Ptr{CnstrMgrPointer}, Cint),
    my_cmp_ref, Dim
)
ccall(
    (:CMGR_CreateCMgr, LIB_CVRPSEP),
    Cvoid,
    (Ptr{CnstrMgrPointer}, Cint),
    old_cmp_ref, Dim
)

my_cmp = my_cmp_ref[]
old_cmp = old_cmp_ref[]

@show my_cmp_ref
@show my_cmp

################################
n_customers = 36
demand = [16, 18, 1, 13, 8, 23, 7, 27, 1, 3, 6, 24, 19, 2, 5, 16, 7, 4, 22, 7, 23, 16, 2, 2, 9, 2, 12, 1, 9, 23, 6, 19, 7, 7, 20, 20]
capacity = 100
edge_tail = [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 13, 13, 13, 13, 13, 13, 13, 13, 13, 13, 13, 13, 13, 13, 13, 13, 13, 13, 13, 13, 13, 13, 13, 13, 14, 14, 14, 14, 14, 14, 14, 14, 14, 14, 14, 14, 14, 14, 14, 14, 14, 14, 14, 14, 14, 14, 14, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 17, 17, 17, 17, 17, 17, 17, 17, 17, 17, 17, 17, 17, 17, 17, 17, 17, 17, 17, 17, 18, 18, 18, 18, 18, 18, 18, 18, 18, 18, 18, 18, 18, 18, 18, 18, 18, 18, 18, 19, 19, 19, 19, 19, 19, 19, 19, 19, 19, 19, 19, 19, 19, 19, 19, 19, 19, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 21, 21, 21, 21, 21, 21, 21, 21, 21, 21, 21, 21, 21, 21, 21, 21, 22, 22, 22, 22, 22, 22, 22, 22, 22, 22, 22, 22, 22, 22, 22, 23, 23, 23, 23, 23, 23, 23, 23, 23, 23, 23, 23, 23, 23, 24, 24, 24, 24, 24, 24, 24, 24, 24, 24, 24, 24, 24, 25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 26, 26, 26, 26, 26, 26, 26, 26, 26, 26, 26, 27, 27, 27, 27, 27, 27, 27, 27, 27, 27, 28, 28, 28, 28, 28, 28, 28, 28, 28, 29, 29, 29, 29, 29, 29, 29, 29, 30, 30, 30, 30, 30, 30, 30, 31, 31, 31, 31, 31, 31, 32, 32, 32, 32, 32, 33, 33, 33, 33, 34, 34, 34, 35, 35, 36]
edge_head = [2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 29, 30, 31, 32, 33, 34, 35, 36, 37, 30, 31, 32, 33, 34, 35, 36, 37, 31, 32, 33, 34, 35, 36, 37, 32, 33, 34, 35, 36, 37, 33, 34, 35, 36, 37, 34, 35, 36, 37, 35, 36, 37, 36, 37, 37]
edge_x = [1.0, 0.0, 2.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 2.0, 2.0, 0.0, 0.0, 0.0, 0.0, 2.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.5, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0]
n_edges = length(edge_head)
max_n_cuts = 10
integrality_tolerance = 1e-6
integer_and_feasible = Ref{Cchar}(-9)
max_violation = Ref{Cdouble}(-9.9)

ccall(
    (:CAPSEP_SeparateCapCuts, LIB_CVRPSEP),
    Cvoid,
    (Cint, Ptr{Cint}, Cint, Cint, Ptr{Cint}, Ptr{Cint}, Ptr{Cdouble}, CnstrMgrPointer, Cint, Cdouble, Ref{Cchar}, Ref{Cdouble}, CnstrMgrPointer),
    n_customers, demand, capacity, n_edges, edge_tail, edge_head, edge_x, old_cmp, max_n_cuts, integrality_tolerance, integer_and_feasible, max_violation, my_cmp
)
@show integer_and_feasible[]
@show max_violation[]

@show old_cmr = unsafe_load(old_cmp)
@show my_cmr = unsafe_load(my_cmp)

@show cpl = my_cmr.CPL

jl_cpl = unsafe_wrap(Array, cpl, my_cmr.Size)

for i in 1:length(jl_cpl)
    cr = unsafe_load(jl_cpl[i])
    # @show cr
    S = unsafe_wrap(Vector{Cint}, cr.IntList, cr.IntListSize)[2:end]
    @show S
    # @show unsafe_wrap(Vector{Cint}, cr.ExtList, cr.ExtListSize)
    # @show unsafe_wrap(Vector{Cint}, cr.CList, cr.CListSize)
    # @show unsafe_wrap(Vector{Cdouble}, cr.CoeffList, cr.CListSize)
    @show cr.RHS
end


text:
"""
    blockmodel(inblock, outblock[, backbone = blockbackbone(inblock)])

From a `backbone` model, construct a model suitable for learning
a mapping from `inblock` to `outblock`.
"""
function blockmodel end

blockmodel(inblock, outblock) = blockmodel(inblock, outblock, blockbackbone(inblock))

"""
    blockbackbone(inblock)

Create a default backbone that takes in block `inblock`.
"""
function blockbackbone end

"""
    blockmodel(inblock::ImageTensor{N}, outblock::OneHotTensor{0}, backbone)
    blockmodel(inblock::ImageTensor{N}, outblock::OneHotTensorMulti{0}, backbone)

Construct a model for N-dimensional image classification. `backbone` should
be a convolutional feature extractor taking in batches of image tensors with
`inblock.nch` color channels.
"""
function blockmodel(
        inblock::ImageTensor{N},
        outblock::Union{OneHotTensor{0}, OneHotTensorMulti{0}},
        backbone) where N
    outsz = Flux.outputsize(backbone, (ntuple(_ -> 256, N)..., inblock.nchannels, 1))
    outch = outsz[end-1]
    head = Models.visionhead(outch, length(outblock.classes), p = 0.)
    return Chain(backbone, head)
end

"""
    blockmodel(inblock::ImageTensor{N}, outblock::OneHotTensor{N}, backbone; kwargs...)

Construct a model for N-dimensional image segmentation. `backbone` should
be a convolutional feature extractor taking in batches of image tensors with
`inblock.nch` color channels. Keyword arguments are passed to [`UNetDynamic`](#).
"""
function blockmodel(inblock::ImageTensor{N}, outblock::OneHotTensor{N}, backbone; kwargs...) where N
    return UNetDynamic(
        backbone,
        (ntuple(_ -> 256, N)..., inblock.nchannels, 1),
        length(outblock.classes);
        kwargs...)
end

"""
    blockmodel(inblock::ImageTensor{N}, outblock::Keypoints{N}, backbone)

Construct a model for image to keypoint regression. `backbone` should
be a convolutional feature extractor taking in batches of image tensors with
`inblock.nch` color channels.
"""
function blockmodel(inblock::ImageTensor{N}, outblock::KeypointTensor{N}, backbone) where N
    outsz = Flux.outputsize(backbone, (ntuple(_ -> 256, N)..., inblock.nchannels, 1))
    outch = outsz[end-1]
    head = Models.visionhead(outch, prod(outblock.sz)*N, p = 0.)
    return Chain(backbone, head)
end

"""
    blockmodel(inblock::TableRow{M, N}, outblock::Union{Continuous, OneHotTensor{0}}, backbone=nothing) where {M, N}

Contruct a model for tabular classification or regression. `backbone` should be a
NamedTuple of categorical, continuous, and a finalclassifier layer, with
the first two taking in batches of corresponding row value matrices.
"""

"""
    blockmodel(::EncodedTableRow, ::OneHotTensor[, backbone])

Create a model for tabular classification. `backbone` should be named tuple
`(categorical = ..., continuous = ...)`. See [`TabularModel`](#) for more info.
"""
function blockmodel(inblock::EncodedTableRow, outblock::OneHotTensor{0}, backbone)
    TabularModel(
        backbone.categorical,
        backbone.continuous,
        Dense(100, length(outblock.classes))
    )
end

"""
    blockmodel(::EncodedTableRow, ::Continuous[, backbone])

Create a model for tabular regression. `backbone` should be named tuple
`(categorical = ..., continuous = ...)`. See [`TabularModel`](#) for more info.
"""
function blockmodel(inblock::EncodedTableRow, outblock::Continuous, backbone)
    TabularModel(
        backbone.categorical,
        backbone.continuous,
        Dense(100, outblock.size)
    )
end

"""
    blockbackbone(ImageTensor{2}(ch))

Construct a XResNet18 model that takes in an encoded image with `ch`
color channels.
"""
blockbackbone(inblock::ImageTensor{2}) = Models.xresnet18(c_in = inblock.nchannels)

function blockbackbone(inblock::EncodedTableRow{M, N}) where {M, N}
    embedszs = Models.get_emb_sz(collect(map(col->length(inblock.categorydict[col]), inblock.catcols)))
    catback = Models.tabular_embedding_backbone(embedszs)
    contback = Models.tabular_continuous_backbone(N)
    return (categorical = catback, continuous = contback)
end

ext: jl


text:
const LEAP_J2000 = round.(Int, (LeapSeconds.LS_EPOCHS .- value(J2000_TO_MJD)) * 86400)
const LEAP_TAI = LEAP_J2000 .+ round.(Int, LeapSeconds.LEAP_SECONDS) .- 1
const LEAP_TAI_SET = Set(LEAP_TAI)

"""
    from_utc(str::AbstractString, dateformat::Dates.DateFormat; scale=TAI)
    from_utc(dt::Dates.DateTime; scale=TAI)
    from_utc(year, month, day, hour=0, minute=0, second=0, fraction=0.0; scale=TAI)
    from_utc(year, month, day, hour, minute, seconds; scale=TAI)

Create an `Epoch` in `scale` based on a UTC timestamp, `Dates.DateTime` or date and
time components.

### Examples ###

```jldoctest; setup = :(using AstroTime; import Dates)
julia> from_utc(2016, 12, 31, 23, 59, 60, 0.0)
2017-01-01T00:00:36.000 TAI

julia> from_utc(2016, 12, 31, 23, 59, 60.0)
2017-01-01T00:00:36.000 TAI

julia> from_utc("2016-12-31T23:59:60.0")
2017-01-01T00:00:36.000 TAI

julia> from_utc("2016-12-31T23:59:60.0", scale=TDB)
2017-01-01T00:01:08.183 TDB
```
"""
from_utc

"""
    to_utc(ep)
    to_utc(::Type{DateTime}, ep)
    to_utc(::Type{Dates.DateTime}, ep)
    to_utc(::Type{String}, ep, dateformat=Dates.default_format(DateTime))

Create a UTC timestamp or `Dates.DateTime` from an `Epoch` `ep`.

### Examples ###

```jldoctest; setup = :(using AstroTime; import Dates)
julia> tai = from_utc(Dates.DateTime(2018, 2, 6, 20, 45, 0, 0))
2018-02-06T20:45:37.000 TAI

julia> to_utc(tai)
"2018-02-06T20:45:00.000"

julia> to_utc(String, tai, Dates.dateformat"yyyy-mm-dd")
"2018-02-06"

julia> to_utc(Dates.DateTime, tai)
2018-02-06T20:45:00
```
"""
to_utc

function from_utc(str::AbstractString,
        dateformat::Dates.DateFormat=Dates.default_format(AstroDates.DateTime);
        scale::TimeScale=TAI)
    dt = AstroDates.DateTime(str, dateformat)
    return from_utc(dt; scale=scale)
end

function from_utc(year::Integer, month::Integer, day::Integer,
    hour::Integer=0, minute::Integer=0, second::Integer=0, fraction=0.0;
    scale::TimeScale=TAI)
    dt = DateTime(year, month, day, hour, minute, second, fraction)
    return from_utc(dt; scale=scale)
end

function from_utc(year::Integer, month::Integer, day::Integer,
    hour::Integer, minute::Integer, seconds;
    scale::TimeScale=TAI)
    dt = DateTime(year, month, day, hour, minute, seconds)
    return from_utc(dt; scale=scale)
end

from_utc(dt::Dates.DateTime; scale::S=TAI) where {S} = from_utc(DateTime(dt); scale=scale)

function from_utc(dt::DateTime; scale::S=TAI) where S
    ep = TAIEpoch(dt)
    idx = searchsortedlast(LEAP_J2000, ep.second)
    if idx == 0
        jd1, jd2 = value.(julian_twopart(ep))
        offset = -offset_utc_tai(jd1, jd2)
    else
        offset = LeapSeconds.LEAP_SECONDS[idx]
        offset = ifelse(second(dt) >= 60, offset - 1, offset)
    end

    return Epoch{S}(TAIEpoch(offset, ep))
end

function to_utc(::Type{DateTime}, ep)
    tai = TAIEpoch(ep)
    idx = searchsortedlast(LEAP_TAI, tai.second)
    if idx == 0
        jd1, jd2 = value.(julian_twopart(tai))
        offset = offset_tai_utc(jd1, jd2)
    else
        offset = LeapSeconds.LEAP_SECONDS[idx]
    end
    dt = DateTime(TAIEpoch(-offset, tai))
    d = Date(dt)
    t = Time(dt)
    if tai.second in LEAP_TAI_SET
        t = Time(hour(t), minute(t), second(t) + 1, fractionofsecond(t))
    end

    return DateTime(d, t)
end

function to_utc(::Type{Dates.DateTime}, ep)
    dt = to_utc(DateTime, ep)
    return Dates.DateTime(dt)
end

function to_utc(::Type{String}, ep, dateformat=Dates.default_format(DateTime))
    dt = to_utc(DateTime, ep)
    return Dates.format(dt, dateformat)
end

to_utc(ep, args...) = to_utc(String, ep, args...)

"""
    now(::Type{Epoch})
    now(::Type{Epoch{S}}) where S<:TimeScale

Get the current date and time as an `Epoch`. The default time scale is TAI.

# Example

```julia
julia> now(Epoch)
2021-04-11T13:20:29.160 TAI

julia> now(TDBEpoch)
2021-04-11T13:21:21.518 TDB
```
"""
function Dates.now(::Type{Epoch{S}}) where S
    return from_utc(Dates.now(Dates.UTC); scale=S())
end

function Dates.now(::Type{Epoch})
    return from_utc(Dates.now(Dates.UTC); scale=TAI)
end

max_stars_repo_name: JuliaAstrodynamics/AstronomicalTime

max_issues_repo_name: JuliaAstrodynamics/AstronomicalTime

max_forks_repo_name: JuliaAstrodynamics/AstronomicalTime


text:
import MathOptInterface
const MOI = MathOptInterface

@testset "MOI wrapper" begin
    model = SatisfiabilityInterface.Optimizer{Bool}()

    x = MOI.add_constrained_variable(model, MOI.ZeroOne())
    y = MOI.add_constrained_variable(model, MOI.ZeroOne())
    z = MOI.add_constrained_variable(model, MOI.ZeroOne())

    MOI.add_constraint(model, MOI.VectorOfVariables([x, y, z]), SatisfiabilityInterface.CNF([true, true, true]))
    MOI.add_constraint(model, MOI.VectorOfVariables([x, y, z]), SatisfiabilityInterface.CNF([false, true, true]))
    MOI.add_constraint(model, MOI.VectorOfVariables([x, y]),    SatisfiabilityInterface.CNF([false, false]))
    MOI.add_constraint(model, MOI.VectorOfVariables([y, z]),    SatisfiabilityInterface.CNF([false, false]))
    MOI.add_constraint(model, MOI.VectorOfVariables([x, z]),    SatisfiabilityInterface.CNF([true, false]))

    MOI.optimize!(model)

    @test MOI.get(model, MOI.RawStatusString()) == "sat"
    @test MOI.get(model, MOI.TerminationStatus()) == MOI.OPTIMAL
    @test MOI.get(model, MOI.PrimalStatus()) == MOI.FEASIBLE_POINT
    @test !MOI.get(model, MOI.VariablePrimal(), x)
    @test MOI.get(model, MOI.VariablePrimal(), y)
    @test !MOI.get(model, MOI.VariablePrimal(), z)
end


text:
using HTTP, JSON, Sockets

include("routes.jl")

struct TimeoutError <: Exception end

function run(host, port)
  println("Start to server: listen to $host:$port")
  HTTP.serve(host, port, reuseaddr=true, rate_limit=100000//1) do request::HTTP.Request
    path = request.target
    println("Request $path")
    payload = Nothing
    try
      payload = JSON.parse(IOBuffer(HTTP.payload(request)))
    catch e
      println("500 Internal Server Error: $e")
      Base.showerror(Base.stdout, e, catch_backtrace())
      return HTTP.Response(500, "Error: $e")
    end

    channel = Channel(spawn=false) do ch
      try
        answer = routes(path, payload)
        isopen(ch) && put!(ch, HTTP.Response(answer))
      catch e
        println("500 $e")
        Base.showerror(Base.stdout, e, catch_backtrace())
        isopen(ch) && put!(ch, HTTP.Response(500, "Error: $e"))
      end
    end

    timedwait(() -> isready(channel), 5.0)

    if !isready(channel)
      isopen(channel) && close(channel)
      println("408 Timeout")
      return HTTP.Response(408, "Error: Request Timeout")
    end
    result = take!(channel)
    isopen(channel) && close(channel)
    return result
  end
end

function main()
  host = "0.0.0.0"
  port = 80
  run(host, port)
end; main()


text:
## RCAM using dedicated owners to store L and R factors with min-norm
# Function Definitions

export dclrMN

function dclrMN(b::AbstractArray{ET}, eta::Number, r::Int, K::Int) where ET<:Number
# Performs RCAM on a single matricised monochromatic slice
# Calls: minFrobQR

    # Get worker IDs
    wids = sort(workers())
    nworkers = length(wids)

    if nworkers < 3
        error(string("DCLR method requires atleast three workers, got ", nworkers))
    end

    # Dedicate two workers to distribute/store L and R
    ownerL = wids[1]
    ownerR = wids[2]

    # Dedicate remaining workers to as solvers
    solvers = wids[3:nworkers]
    nsolvers = length(solvers)

    # Give this info to every solver
    [sendto(pid, ownerL = ownerL, ownerR = ownerR) for pid in procs()] #change to procs()[solvers]

    # Get size of data inorder to generate L0
    m::Int,n::Int = size(b)

    #Create L and R on their owners
    sendto(ownerL, L=collect(ET, randn(m,r)))
    sendto(ownerR, R=zeros(ET,n,r))

    # Find non-zero inds for each col, reduce data into collections
    indn = (1:n)::UnitRange{Int64}

    # This global declaration is a workaround to open issue #15276
    indnzd = distribute([findall(!iszero,b[:,i]) for i::Int in indn],
                        procs = procs()[solvers])
    bsubd = distribute([b[indnzd[i],i] for i in indn],
                        procs = procs()[solvers])

    # Find non zero inds for each row, reduce data into collections
    indm = (1:m)::UnitRange{Int64}
    indnz_trand = distribute([findall(!iszero,b[i,:]) for i::Int in indm],
                        procs = procs()[solvers])
    bsub_trand = distribute([b[i,indnz_trand[i]] for i in indm],
                        procs = procs()[solvers])

    # Clear data
    b = []

    # Begin Alterations
    for k = 1:K

        @time @sync for pid in procs()[solvers]
            @async fetch(@spawnat pid altupdateR!(ownerL, ownerR, bsubd, indnzd, eta))
        end

        @time @sync for pid in procs()[solvers]
            @async fetch(@spawnat pid altupdateL!(ownerL, ownerR, bsub_trand,indnz_trand, eta))
        end

    end

   return ownerL, ownerR

end

function altupdateL!(ownerL::Int, ownerR::Int, bsubd::AbstractArray{AT}, indnzd, eta::Number) where AT

    ET = eltype(AT)

    #Get info about local part of data
    m = localindices(bsubd)[1]

    # Get updated version of R
    R_tmp = fetch(@spawnat(ownerR, Main.R))

    # Update the ith row of Rd
    for i in m

        local_i = i - m[1] + 1

        # Ask for required part of L
        req_rows = DistributedArrays.localpart(indnzd)[local_i]

        # Solve for what this row or R will be updated to
        L_update = collect(ET,conj(pminnorm(R_tmp[req_rows,:],
                            DistributedArrays.localpart(bsubd)[local_i], eta)))

        # Create a future for this result
        L_future = Future(ownerL)
        put!(L_future, L_update)

        # Update the i-th row of R with the future
        @spawnat(ownerL, setindex!(Main.L, fetch(L_future), i, :))

    end

end

function altupdateR!(ownerL::Int, ownerR::Int, bsubd::AbstractArray{AT}, indnzd, eta::Number) where AT

    ET = eltype(AT)

    #Get info about local part of data
    m = localindices(bsubd)[1]

    # Get updated version of L
    L_tmp = fetch(@spawnat(ownerL, Main.L))

    # Update the ith row of Rd
    for i in m

        local_i = i - m[1] + 1

        # Ask for required part of L
        req_rows = DistributedArrays.localpart(indnzd)[local_i]

        # Solve for what this row or R will be updated to
        R_update = collect(ET,conj(pminnorm(L_tmp[req_rows,:],
                            DistributedArrays.localpart(bsubd)[local_i], eta)))

        # Create a future for this result
        R_future = Future(ownerR)
        put!(R_future, R_update)

        # Update the i-th row of R with the future
        @spawnat(ownerR, setindex!(Main.R, fetch(R_future), i, :))

    end

end

function pminnorm( A::TA, b::Tb, eta::Number) where {ETA<:Number, ETb<:Number, TA<:AbstractArray{ETA}, Tb<:AbstractArray{ETb}}
# Solves a single row, meant to be mapped over a DArray

    Q::TA = svd(A).U
    d = eta^2 - norm(b - Q*Q'*b)^2

    # Force type of v
    v::Array{ETA,1} = Array{ETA,1}()
    if d < 0
        v = A\b
    else

    G = Q'*A
    H = G*G'
    b2 = Q'*b
    v1 = norm(H\b2)
    c = sqrt(d)/v1
    tmp = G'*G
    M = tmp + c*Matrix{eltype(tmp)}(I, size(tmp)[1], size(tmp)[2])

    tmp2 = G'*b2
    v = M\tmp2

    end
    return v
end
export pminnorm


text:
################################################################################
#
#   KAGOME LATTICE
#
################################################################################

# Implementation
# - implementation 1
# - labels <: Any
function getUnitcellKagome(
            unitcell_type  :: Type{U},
            implementation :: Val{1}
        ) :: U where {LS,LB,S<:AbstractSite{LS,2},B<:AbstractBond{LB,2}, U<:AbstractUnitcell{S,B}}

    # return a new Unitcell
    return newUnitcell(
        # Type of the unitcell
        U,
        # Bravais lattice vectors
        Vector{Float64}[
            Float64[sqrt(3.0)/2, -0.5],
            Float64[sqrt(3.0)/2, +0.5]
        ],
        # Sites
        S[
            newSite(S, Float64[     0.0   ,  0.0 ], getDefaultLabelN(LS,1)),
            newSite(S, Float64[sqrt(3.0)/4, -0.25], getDefaultLabelN(LS,2)),
            newSite(S, Float64[sqrt(3.0)/4, +0.25], getDefaultLabelN(LS,3))
        ],
        # Bonds
        B[
            newBond(B, 1, 2, getDefaultLabel(LB), ( 0, 0)),
            newBond(B, 1, 2, getDefaultLabel(LB), (-1, 0)),
            newBond(B, 1, 3, getDefaultLabel(LB), ( 0, 0)),
            newBond(B, 1, 3, getDefaultLabel(LB), ( 0,-1)),
            newBond(B, 2, 1, getDefaultLabel(LB), ( 0, 0)),
            newBond(B, 2, 1, getDefaultLabel(LB), (+1, 0)),
            newBond(B, 2, 3, getDefaultLabel(LB), ( 0, 0)),
            newBond(B, 2, 3, getDefaultLabel(LB), (+1,-1)),
            newBond(B, 3, 1, getDefaultLabel(LB), ( 0, 0)),
            newBond(B, 3, 1, getDefaultLabel(LB), ( 0,+1)),
            newBond(B, 3, 2, getDefaultLabel(LB), ( 0, 0)),
            newBond(B, 3, 2, getDefaultLabel(LB), (-1,+1))
        ]
    )
end


text:
###############################################################################
# other stuff, trivia and non-trivia
###############################################################################

function Nemo.PolynomialRing(R::Nemo.Ring, n::Int, s::String="x";
                             cached::Bool = false, ordering::Symbol = :lex)
  return Nemo.PolynomialRing(R, ["$s$i" for i=1:n], cached = cached,
                                                    ordering = ordering)
end

#TODO: makes only sense if f is univ (uses only one var)
function (Rx::FmpzPolyRing)(f::fmpq_mpoly)
  fp = Rx()
  R = base_ring(Rx)
  d = denominator(f)
  @assert d == 1
  for t = terms(f)
    e = total_degree(t)
    @assert length(t) == 1
    c = coeff(t, 1)
    setcoeff!(fp, e, numerator(c*d))
  end
  return fp
end

function (Rx::FmpqPolyRing)(f::fmpq_mpoly)
  fp = Rx()
  R = base_ring(Rx)
  for t = terms(f)
    e = total_degree(t)
    @assert length(t) == 1
    c = coeff(t, 1)
    setcoeff!(fp, e, c)
  end
  return fp
end

function (Rx::GFPPolyRing)(f::fmpq_mpoly)
  fp = Rx()
  R = base_ring(Rx)
  d = denominator(f)
  for t = terms(f)
    e = total_degree(t)
    @assert length(t) == 1
    c = coeff(t, 1)
    setcoeff!(fp, e, R(numerator(c*d)))
  end
  return fp * inv(R(d))
end

function derivative(f::Generic.Frac{T}, x::T) where {T <: MPolyElem}
  return derivative(f, var_index(x))
end

function derivative(f::Generic.Frac{T}, i::Int) where {T <: MPolyElem}
  n = numerator(f)
  d = denominator(f)
  return (derivative(n, i)*d - n*derivative(d, i))//d^2
end

function evaluate(f::Generic.Frac{T}, V::Vector{U}) where {T <: RingElem, U <: RingElem}
  return evaluate(numerator(f), V)//evaluate(denominator(f), V)
end

function evaluate(f::Generic.Frac{T}, v::U) where {T <: RingElem, U <: RingElem}
  return evaluate(numerator(f), v)//evaluate(denominator(f), v)
end

function Hecke.lead(f::AbstractAlgebra.MPolyElem)
  iszero(f) && error("zero poly")
  return coeff(f, 1)
end

function Hecke.coefficients(f::AbstractAlgebra.MPolyElem)
  return Hecke.PolyCoeffs(f)
end

function Base.iterate(PC::Hecke.PolyCoeffs{<:AbstractAlgebra.MPolyElem}, st::Int = 0)
  st += 1
  if st > length(PC.f)
    return nothing
  else
    return coeff(PC.f, st), st
  end
end

Base.length(M::Hecke.PolyCoeffs{<:AbstractAlgebra.MPolyElem}) = length(M.f)


text:
module Transformers

using Flux
using Requires
using Requires: @init

export Transformer, TransformerDecoder
export Stack, @nntopo_str, @nntopo

export dataset, datafile, get_batch, get_vocab

export todevice
export Gpt
export Bert

const Abstract3DTensor{T} = AbstractArray{T, 3}
const Container{T} = Union{NTuple{N, T}, Vector{T}} where N

"move data to device, when CuArrays is loaded, basically = `CuArrays.cu` except `AbstractArray{Int}` become `CuArray{Int}`"
todevice(x) = x
todevice(x, xs...) = (x, xs...)

@init @require CuArrays="3a865a2d-5b23-5a0f-bc46-62713ec82fae" begin
    import .CuArrays

    todevice(x) = CuArrays.cu(x)
    todevice(x, xs...) = (todevice(x), todevice.(xs)...)
    todevice(x::Union{Tuple, NamedTuple}) = map(todevice, x)
    todevice(x::AbstractArray{Int}) = CuArrays.CuArray(x)
end

#implement batchmul for flux
include("./fix/batchedmul.jl")

#dropout noise shape impl
include("./fix/dropout.jl")

#scatter/gather with atomic ops
include("./fix/atomic.jl")
include("./fix/scatter.jl")

include("./basic/Basic.jl")
include("./stacks/Stacks.jl")
include("./datasets/Datasets.jl")

include("./pretrain/Pretrain.jl")

include("./gpt/GenerativePreTrain.jl")
include("./bert/BidirectionalEncoder.jl")

using .Basic
using .Stacks
using .Datasets
using .Pretrain
using .GenerativePreTrain
using .BidirectionalEncoder

end # module


text:
#****************************************************************************
#                     Molecular Dynamics Potentials (MDP)
#                            CESMIX-MIT Project
#
#  Contributing authors: Ngoc-Cuong Nguyen (cuongng@mit.edu, exapde@gmail.com)
#****************************************************************************

function writeconfig(config, filename)

# disp('Writing configurations into a binary file...');

nconfigs = config.nconfigs; # number of configurations
dim = config.dim; # physical dimension
ncx = config.ncx; # number of compoments of x
ncv = config.ncv; # number of compoments of v
nce = config.nce; # number of compoments of e
ncf = config.ncf; # number of compoments of f
ncq = config.ncq; # number of compoments of q

tmp = [nconfigs dim ncq ncv nce ncf];

nsize = zeros(20,1);
nsize[1] = length(tmp[:]);
nsize[2] = length(config.natom[:]);
nsize[3] = length(config.a[:]);
nsize[4] = length(config.b[:]);
nsize[5] = length(config.c[:]);
nsize[6] = length(config.e[:]);
nsize[7] = length(config.t[:]);
nsize[8] = length(config.x[:]);
nsize[9] = length(config.q[:]);
nsize[10] = length(config.v[:]);
nsize[11] = length(config.f[:]);
nsize[12] = length(config.we[:]);
nsize[13] = length(config.wf[:]);

fileID = open(filename,"w");

write(fileID,Float64(length(nsize[:])));
write(fileID,Float64.(nsize[:]));
write(fileID,Float64.(tmp[:]));
write(fileID,Float64.(config.natom));
write(fileID,Float64.(config.a));
write(fileID,Float64.(config.b));
write(fileID,Float64.(config.c));
write(fileID,Float64.(config.e));
write(fileID,Float64.(config.t));
write(fileID,Float64.(config.x));
write(fileID,Float64.(config.q));
write(fileID,Float64.(config.v));
write(fileID,Float64.(config.f));
write(fileID,Float64.(config.we));
write(fileID,Float64.(config.wf));
close(fileID);

end

max_stars_repo_name: cesmix-mit/MDP

max_issues_repo_name: cesmix-mit/MDP

max_forks_repo_name: cesmix-mit/MDP


text:
using DifferentialEquations
using Plots

function macrop_ode(dY,Y,p,t)
 dY[1] = (p[1]-p[2])*Y[1] - p[3]*Y[2]
 dY[2] = p[4]*Y[1]*Y[3] - (p[5]+p[3]+p[2])*Y[2] - (p[3]*((Y[2]^2)/Y[1])*((p[6]+1)/p[6]))
 dY[3] = p[7]*Y[2] - (p[8]*Y[3]) - (p[4]*Y[1]*Y[3])
end

par=[1.4,1.05,0.0003,0.01,0.5,0.1,10.0,10.0]
init=[100.0,10.0,10.0]
tspan=(0.0,100.0)

macro_par = ODEProblem(macrop_ode,init,tspan,par)

sol=solve(macro_par);

# Plot
plot(sol,xlabel="Time",yscale=:log10)

max_stars_repo_name: epimodels/epicookbook

max_issues_repo_name: epimodels/epicookbook

max_forks_repo_name: epimodels/epicookbook


text:
"""
"""
function print_vector(x::Vector{Float64}, msg::String = "")
    @printf("%s\n", msg)
    for (i,v) in enumerate(x)
        @printf("    %+.6f", v)
        if i % 5 == 0
            @printf("\n")
        end
    end
    @printf("\n")
end

"""
"""
function dropzeros!(x::Vector{Float64}, eps::Float64 = 1.0e-10)
    for (i,v) in enumerate(x)
        if abs(v) < eps
            x[i] = 0.0
        end
    end
end

function print_matrix(A::SparseMatrixCSC{Float64, Int64})
    for i = 1:A.m
        a = A[i,:]
        SparseArrays.droptol!(a, 1.0e-10)
        if length(a.nzind) > 0
            @printf("row%6d", i)
            for (k,j) in enumerate(a.nzind)
                @printf("\tcol%6d\t%+.6e", j, a.nzval[k])
            end
            @printf("\n")
        end
    end
end


text:
"""
    bounded_by(f, a::Arf, b::Arf, C::Arf; degree, abs_value, log_bisection, depth_start, maxevals, depth, threaded, verbose)

Return `true` if the function `f` can be shown to be bounded by `C` on
the interval `[a, b]`, i.e. `f(x) <= C` for all `x ∈ [a, b]`,
otherwise return `false`.

This function is similar to first computing the maximum with
[`maximum_enclosure`](@ref) and then check if the computed maximum
satisfies the bound. However if the only thing needed is to check if
the bound holds this method has a number of benefits
- it aborts early if the bound is shown to not hold
- it doesn't try to compute an accurate enclosure of the maximum, it
  only bisects as much as is needed for getting the bound.
- the implementation is simpler and easier to check for correctness.

The maximum of `f` is enclosed by the use of Taylor series through
[`maximum_series`](@ref). The degree of the expansion is kept constant
and the interval is bisected until we can either conclude that the
bound holds on the whole interval or there is subinterval where it
doesn't hold.

The function `f` should support evaluation on both `Arb` and
`ArbSeries` and should return an enclosure of the result in both
cases. The degree used for `ArbSeries` can be set with the `degree`
argument (defaults to 8). If `degree` is negative then it will fall
back to direct evaluation with `Arb` and not make use of
[`maximum_series`](@ref), this is usually much slower but does not
require the function to be implemented for `ArbSeries`.

If `abs_value = true` the instead consider the function `abs(f(x))` on
the interval `[a, b]`.

If `log_bisection = true` then the intervals are bisected in a
logarithmic scale, see [`bisect_interval`](@ref) for details.

The argument `depth_start` bisect the interval using
[`bisect_interval_recursive`](@ref) before starting to compute the
extrema. This can be useful if it is known beforehand that a certain
number of bisections will be necessary before the enclosures get good
enough. It defaults to `0` which corresponds to not bisecting the
interval at all before starting.

The arguments `maxevals` and `depth` can be used to limit the number
of function evaluations and the number of bisections of the interval
respectively. Notice that `depth` takes `depth_start` into account, so
the maximum number of iterations is `depth - depth_start`.

If `threaded = true` then evaluate `f` in parallel on the intervals
using [`Threads.@threads`](@ref).

If `verbose = true` then output information about the process.
"""
function bounded_by(
    f,
    a::Arf,
    b::Arf,
    C::Arf;
    degree::Integer = 8,
    abs_value = false,
    log_bisection = false,
    depth_start::Integer = 0,
    maxevals::Integer = 1000,
    depth::Integer = 20,
    threaded = false,
    verbose = false,
)
    check_interval(a, b)

    maybe_abs = abs_value ? abs : identity

    if a == b
        return maybe_abs(f(Arb(a))) <= C
    end

    # List of intervals
    intervals = bisect_interval_recursive(a, b, depth_start, log_midpoint = log_bisection)

    iterations = 0
    evals = 0
    while true
        iterations += 1
        evals += length(intervals)

        # Evaluate f on the remaining intervals
        values = similar(intervals, Arb)
        if degree < 0
            if threaded
                Threads.@threads for i in eachindex(intervals)
                    values[i] = maybe_abs(f(Arb(intervals[i])))
                end
            else
                for i in eachindex(intervals)
                    values[i] = maybe_abs(f(Arb(intervals[i])))
                end
            end
        else
            if threaded
                Threads.@threads for i in eachindex(intervals)
                    values[i], _ = maximum_series(f, intervals[i]...; degree, abs_value)
                end
            else
                for (i, (a, b)) in enumerate(intervals)
                    values[i], _ = maximum_series(f, a, b; degree, abs_value)
                end
            end
        end

        # Check for each interval if the bound is satisfied
        next_intervals = sizehint!(empty(intervals), 2length(intervals))
        for i in eachindex(intervals)
            if values[i] > C
                verbose && @info "bound doesn't hold on the interval x = $(intervals[i])"
                verbose && @info "got the maximum of f(x) to be $(values[i])"
                return false
            elseif !(values[i] <= C)
                push!(
                    next_intervals,
                    bisect_interval(intervals[i]..., log_midpoint = log_bisection)...,
                )
            end
        end
        intervals = next_intervals

        verbose && @info "iteration: $(lpad(iterations, 2)), " *
              "remaining intervals: $(lpad(length(intervals) ÷ 2, 3)), "

        non_finite = count(!isfinite, values)
        verbose && non_finite > 0 && @info "non-finite intervals: $non_finite"

        isempty(intervals) && break

        # Check if we have done the maximum number of function
        # evaluations or reached the maximum depth
        if evals >= maxevals || iterations >= depth - depth_start
            if verbose
                evals >= maxevals &&
                    @warn "reached maximum number of evaluations $evals >= $maxevals"
                iterations >= depth - depth_start && @warn "reached maximum depth $depth"
            end
            return false
        end
    end

    return true
end


text:
# NOTE: We have to circumvent JuMP.NLexpressions to be able to support nonlinear expressions
#       without breaking everything...
#       Thus why all operations on JuMP-incompatible functions is on Exprs.

"""
    fetch_variable(model::JuMP.Model, varkey::Union{Symbol, String, VariableRef, Array})

Returns JuMP.VariableRefs that match a given Symbol, String, VariableRef,
or array of these.
"""
function fetch_variable(model::JuMP.Model, varkey::Union{Symbol, String, VariableRef})
    if varkey isa Symbol
        try
            return model[varkey]
        catch UndefVarError
            throw(KeyError(string("Varkey ", varkey,
                              " is invalid for Model.")))
        end
    elseif varkey isa VariableRef
        is_valid(model, varkey) && return varkey
        throw(KeyError(string("Varkey ", varkey,
                              " is invalid for Model.")))
    elseif varkey isa String
        ret = JuMP.variable_by_name(model, varkey)
        !isnothing(ret) && return ret
        throw(KeyError(string("Varkey ", varkey,
                              " is invalid for Model.")))
    end
end

function fetch_variable(model::JuMP.Model, varkey::Array)
    return [fetch_variable(model, key) for key in varkey]
end

""" Helper function for finding bounds. """
function get_bound(var::JuMP.VariableRef)
    if JuMP.has_lower_bound(var)
        if JuMP.has_upper_bound(var)
            return [JuMP.lower_bound(var), JuMP.upper_bound(var)]
        else
            return [JuMP.lower_bound(var), Inf]
        end
    else
        if JuMP.has_upper_bound(var)
            return [-Inf, JuMP.upper_bound(var)]
        else
            return [-Inf, Inf]
        end
    end
end

"""
    get_bounds(var::JuMP.VariableRef})
    get_bounds(vars::Array{JuMP.VariableRef})

Returns bounds of JuMP variables.
"""
function get_bounds(var::JuMP.VariableRef)
    return var => get_bound(var)
end

get_bounds(vars::Array{JuMP.VariableRef}) = Dict(get_bounds(var) for var in vars)

"""
    get_unbounds(var::JuMP.VariableRef)
    get_unbounds(vars::Array{VariableRef})
    get_unbounds(gm::Union{JuMP.Model, GlobalModel})

Returns variables with no lower and/or upper bounds.
"""
function get_unbounds(var::JuMP.VariableRef)
   if JuMP.has_lower_bound(var)
        if !JuMP.has_upper_bound(var)
            return var => [JuMP.lower_bound(var), Inf]
        else
            return
        end
    else
        if JuMP.has_upper_bound(var)
            return var => [-Inf, JuMP.upper_bound(var)]
        else
            return var => [-Inf, Inf]
        end
    end
end

function get_unbounds(vars::Array{JuMP.VariableRef})
     unbounds = Dict()
     for var in vars
        res = get_unbounds(var)
        if !isa(res, Nothing)
            unbounds[res.first] = res.second
        end
    end
    if length(unbounds) == 0
        return
    else
        return unbounds
    end
end

"""
    data_to_DataFrame(data::Union{Dict, DataFrame, DataFrameRow})

Gets data with different keys, and returns a DataFrame with string headers.
"""
function data_to_DataFrame(data::Union{Dict, DataFrame, DataFrameRow})
    if data isa DataFrame
        return data
    elseif data isa DataFrameRow
        return DataFrame(data)
    else
        newdata = Dict(string(key) => val for (key, val) in data)
        return DataFrame(newdata)
    end
end

"""
    distance_to_set(val::Union{Array{<:Real},<:Real}, set::MOI.AbstractSet)

Wrapper around MathOptSetDistances.distance_to_set.
Distance 0 if val ∈ set. Otherwise, returns Float64.
"""
function distance_to_set(val::Union{Array{<:Real},<:Real}, set::MOI.AbstractSet)

    return MathOptSetDistances.distance_to_set(MathOptSetDistances.DefaultDistance(), val, set)
end

"""
    $(TYPEDSIGNATURES)

Returns constant of MOI.AbstractSet, if it exists.
"""
function get_constant(set::MOI.AbstractSet)
    try
        MOI.constant(set)
    catch
        nothing
    end
end

"""
    $(TYPEDSIGNATURES)

Returns an "evaluate-able" function from an Expr, or nothing for a JuMP.ConstraintRef.
Can extend in the future to other elements.
"""
function functionify(constraint)
    if constraint isa Expr
        f = eval(constraint)
        f isa Function && return f
        throw(OCTHaGOnException(string("functionify(", f, ") is not a valid function.")))
    else
        return nothing
    end
end

"""
    vars_from_expr(expr::Expression, model::JuMP.Model)

Returns the JuMP Variables that are associated with a given expression.
Note: Expr's must be defined with a single input or a tuple of inputs, eg:

    ex = :(x -> 5*x)
    ex = :((x, y, z) -> sum(x[i] for i=1:4) - y[1] * y[2] + z)
"""
function vars_from_expr(expr::Expr, model::JuMP.Model)
    if expr.args[1] isa Symbol
        return [fetch_variable(model, expr.args[1])]
    else
        return [fetch_variable(model, arg) for arg in expr.args[1].args]
    end
end

"""
    vars_from_constraint(con::JuMP.ConstraintRef)

Returns the JuMP Variables that are associated with a JuMP.ConstraintRef.
Note: Currently only works for affine and quadratic constraints.
"""
function vars_from_constraint(con::JuMP.ConstraintRef)
    confunc = constraint_object(con).func
    if confunc isa JuMP.GenericAffExpr
        return collect(keys(confunc.terms))
    elseif confunc isa JuMP.GenericQuadExpr
        vars = [term for (term,val) in confunc.aff.terms]
        for (term, val) in confunc.terms
            append!(vars, [term.a, term.b])
        end
        return unique(vars)
    elseif confunc isa VariableRef
        return [confunc]
    else
        throw(OCTHaGOnException("Only affine or quadratic constraints are current supported."))
    end
end

"""
    gradientify(expr::Expr, expr_vars::Array)
    gradientify(expr::JuMP.ConstraintRef, expr_vars::Array)

Turns an expression into a gradient-able (via ForwardDiff), flattened function.
TODO: add option to not use gradients for certain functions or data driven constraints.
"""
function gradientify(expr::Expr, expr_vars::Array)
    var_ranges = get_var_ranges(expr_vars)
    f = functionify(expr)
    gradable_fn = x -> Base.invokelatest(f, [x[i] for i in var_ranges]...)
    return x -> ForwardDiff.gradient(gradable_fn, x)
end

function gradientify(con::JuMP.ConstraintRef, expr_vars::Array)
    confunc = constraint_object(con).func
    if confunc isa JuMP.GenericAffExpr
        gradvals = [confunc.terms[var] for var in flat(expr_vars)]
        gradfn = let gvals = gradvals
            function (x)
                return gradvals
            end
        end
        return gradfn
    elseif confunc isa JuMP.GenericQuadExpr
        flatvars = flat(expr_vars)
        garr = zeros(length(flatvars))
        gdict = Dict()
        for (term, val) in confunc.aff.terms # Affine component
            idx = findall(x -> x == term, flatvars)
            if !isempty(idx)
                garr[idx[1]] += val
            end
        end
        for (term, val) in confunc.terms # Quadratic component
            idx1 = findall(x -> x == term.a, flatvars)[1]
            idx2 = findall(x -> x == term.b, flatvars)[1]
            gdict[(idx1, idx2)] = val
        end
        gradfn = let gdict = gdict, garr = garr
            function (x)
                q_terms = zeros(length(garr))
                for (term, val) in gdict
                    q_terms[term[1]] += val * x[term[2]]
                    q_terms[term[2]] += val * x[term[1]]
                end
                return garr + q_terms
            end
        end
        return gradfn
    else
        throw(
        OCTHaGOnException("Currently, only supporting gradients of affine or quadratic JuMP constraints. Please " *
                     "submit your constraint as a Expr instead and try again. "))
    end
end

gradientify(con::Nothing, expr_vars::Array) = nothing

"""
    linearize_objective!(model::JuMP.Model)
Makes sure that the objective function is affine.
"""
function linearize_objective!(model::JuMP.Model)
    objtype = JuMP.objective_function(model)
    objsense = string(JuMP.objective_sense(model))
    if objtype isa Union{VariableRef, GenericAffExpr} || objsense == "FEASIBILITY_SENSE"
        return
    else
        aux = @variable(model)
        @objective(model, Min, aux)
        # Default optimization problem is always a minimization
        coeff = 1;
        objsense == "MAX_SENSE" && (coeff = -1)
        try
            @constraint(model, aux >= coeff*JuMP.objective_function(model))
        catch
            @NLconstraint(model, aux >= coeff*JuMP.objective_function(model))
        end
        return
    end
end

""" Checks whether any defined bounds are infeasible by given Model. """
function check_infeasible_bound(bound::Pair)
    key = bound.first
    val = bound.second
    @assert key isa JuMP.VariableRef
    @assert val isa Array
    model_bounds = get_bound(key)
    if minimum(model_bounds) >= maximum(val) || maximum(model_bounds) <= minimum(val)
        throw(OCTHaGOnException("Infeasible bounds."))
    else
        return [maximum([minimum(model_bounds), minimum(val)]),
                    minimum([maximum(model_bounds), maximum(val)])]
    end
end

"""
    bound!(model::JuMP.Model, bound::Pair)
    bound!(model::JuMP.Model, bounds::Dict)
    bound!(model::GlobalModel, bounds::Union{Pair,Dict})

Adds outer bounds to JuMP Model from Dict or Pair of data.
"""
function bound!(model::JuMP.Model, bound::Pair)
    if isa(bound.first, JuMP.VariableRef)
        tightest_bound = check_infeasible_bound(bound)
        b_min = minimum(tightest_bound)
        b_max = maximum(tightest_bound)
        !isinf(b_min) && JuMP.set_lower_bound(bound.first, b_min)
        !isinf(b_max) && JuMP.set_upper_bound(bound.first, b_max)
    else
        vars = fetch_variable(model, bound.first)
        if vars isa JuMP.VariableRef
            bound!(model, vars => bound.second)
        elseif vars isa Array
            for var in vars
                bound!(model, var => bound.second)
            end
        else
            throw(OCTHaGOnException("Bound with fetch_variable has failed. Try bounding using JuMP.VariableRefs!"))
        end
    end
    return
end

function bound!(model::JuMP.Model, bounds::Dict)
    for bd in collect(bounds)
        bound!(model, bd)
    end
end

"""
    restrict_to_set(var::JuMP.VariableRef, s::Union{Set, Array})

Restricts variable to a set s. Useful for non-integer sets or when taking log of integer variables.
"""
function restrict_to_set(var::JuMP.VariableRef, s::Union{Set, Array})
    int = @variable(var.model, [1:length(s)], Bin)
    @constraint(var.model, sum(int) == 1)
    @constraint(var.model, var == sum(s .* int))
    return
end

""" Counts the number of constraints in an optimization model."""
function count_constraints(model::JuMP.Model)
    types = JuMP.list_of_constraint_types(model)
    return sum(length(all_constraints(model, type[1], type[2])) for type in types)
end

""" Counts the number of constraint types in an optimization model."""
function count_types(model::JuMP.Model)
    return length(JuMP.list_of_constraint_types(model))
end

""" Counts the number of variables in an optimization model. """
function count_variables(model::JuMP.Model)
    return length(JuMP.all_variables(model))
end


text:
using VideoIO, Images, DifferentialEquations, DataStructures, ImageTransformations
include("control.jl")
include("cal_contr.jl")

"""
    _start()

Main function which starts the video feed
    and calls subsequent functions
"""
function _start(display, arsol)
    cam = VideoIO.opencamera()

    # set the previous brightness to 1
    prevBrightness = 1.0
    _setBrightness(display, prevBrightness)

    n = 0
    while true
        n = (n+1)%10
        rgbImg = read(cam)
        if n != 0
            continue
        end

        # decrease the image size by a factor of 8 to
        # decrease the load of computation
        rgbImg = imresize(rgbImg, ratio=1/8)

        # convet rgb to hsl color space
        singleChannel = _processImage(rgbImg, "HSL")

        # Find what brightness should have been at this stage
        toBeBrightness = max(0.05, min(1.0, sum(singleChannel)/length(singleChannel)))

        # find the start and end from solution
        endIndex = searchsortedfirst(arsol, toBeBrightness*10)
        startIndex = searchsortedfirst(arsol, prevBrightness*10)

        # Change the brightness gradually
        if startIndex < endIndex
            for i in startIndex:endIndex
                _setBrightness(display, min(arsol[i]/10 + 0.3, 1.0))
            end
        else
            for i in endIndex:startIndex
                _setBrightness(display, min(arsol[i]/10 + 0.3, 1.0))
            end
        end
        # _setBrightness(display, brightness)
        prevBrightness = toBeBrightness
    end
end

"""
    _processImage()

Handles the image processing part which will
be taken by contrast finder
pType: Eithe gray or HSL with default value gray
returns: A single channel matrix
"""
function _processImage(img, pType="gray")
    if pType == "gray"
        return Gray.(img)
    elseif pType == "HSL"
        hslImg = HSL.(img)
        channels = channelview(float.(hslImg))
        return channels[3,:,:]
    end
end

# display = _getDisplay()
# _start(display)


text:
# Autogenerated wrapper script for CImGui_jll for i686-w64-mingw32
export compile_commands, libcimgui, libimgui

JLLWrappers.@generate_wrapper_header("CImGui")
JLLWrappers.@declare_file_product(compile_commands)
JLLWrappers.@declare_library_product(libcimgui, "libcimgui.dll")
JLLWrappers.@declare_library_product(libimgui, "libimgui-cpp.dll")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_file_product(
        compile_commands,
        "share\\compile_commands.json",
    )

    JLLWrappers.@init_library_product(
        libcimgui,
        "bin\\libcimgui.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libimgui,
        "bin\\libimgui-cpp.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
"""

`SymPy` package to interface with Python's [SymPy library](http://www.sympy.org) through `PyCall`.

The basic idea is that a new type -- `Sym` -- is made to hold symbolic
objects.  For this type, the basic functions from SymPy and appropriate functions
of `Julia` are overloaded for `Sym` objects so that the expressions
are treated symbolically and not evaluated immediately. Instances of
this type are created by the constructor `Sym`, the function `symbols` or the macro
`@vars`.

On loading, a priviledged set of the functions from the `sympy` module
are defined as generic functions with their first argument narrowed to
symbolic types. Others may be accessed by qualification, as in
`sympy.trigsimp`. Calling `import_from(sympy)` will import the
rest. SymPy methods are called through Python's dot-call syntax.  To
find documentation on SymPy functions and methods, one should refer to
SymPy's [website](http://docs.sympy.org/latest/index.html).

Plotting is provided through  `Plots` recipes. For details, see
the help page for [`sympy_plotting`](@ref).

The package documentation provides many examples.

"""
module SymPy

if isdefined(Base, :Experimental) && isdefined(Base.Experimental, Symbol("@optlevel"))
    @eval Base.Experimental.@optlevel 1
end

## minimal version of conveniently using SymPy in Julia
## uses getfield overloading to access sympy methods of Sym object

using PyCall

using SpecialFunctions
using LinearAlgebra
using Markdown

import Base: show
import Base: convert, promote_rule
import Base: getproperty
import Base: hash, ==
import Base: length, size
import Base.iterate
import Base: +, -, *, /, //, \, ^

export @vars, Sym, sympify, symbols, @symfuns, @syms
export SymMatrix, SymFunction
export PI, IM, oo, zoo, True, False
export N, subs

export sympy, sympy_core, sympy_matrices, import_from#, import_sympy
export free_symbols

include("types.jl")
include("decl.jl")
include("constructors.jl")
include("utils.jl")
include("numbers.jl")
include("mathops.jl")
include("mathfuns.jl")
include("generic.jl")
include("matrix.jl")
include("sets.jl")
include("symfunction.jl")
include("assumptions.jl")
include("lambdify.jl")
include("patternmatch.jl")
include("permutations.jl")
include("plot_recipes.jl")

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

pynull() = PyCall.PyNULL()
"""
    sympy

Variable from `pyimport("sympy")`. Numerous methods are available through Python's dot-call syntax.
"""
const sympy  = PyCall.PyNULL()
"""
    sympy_core

Variable from `pyimport("sympy.core")`.
"""
const sympy_core =  PyCall.PyNULL()
"""
    sympy_matrices

Variable from `pyimport("sympy.matrices")`.
"""
const sympy_matrices =  PyCall.PyNULL()
const mpmath = PyCall.PyNULL()
const combinatorics  = PyCall.PyNULL()

# core.sympy.numbers.*
"PI is symbolic `pi`"
global PI = Sym(pynull())
"IM is a symbolic `im`"
global IM = Sym(pynull())
"oo is a symbolic infinity. Example: `integrate(exp(-x), x, 0, oo)`."
global oo = Sym(pynull())
"zoo complex inifinity"
global zoo = Sym(pynull())
"True from SymPy"
global True = Sym(pynull())
"False from SymPy"
global False = Sym(pynull())

# Can not actually initiate many things until `sympy` is defined, so not until runtime
function __init__()

    ## Define sympy, mpmath, ...
    copy!(sympy, PyCall.pyimport_conda("sympy", "sympy"))
    copy!(sympy_core, PyCall.pyimport("sympy.core"))
    copy!(sympy_matrices, PyCall.pyimport("sympy.matrices"))
    copy!(PI.__pyobject__,  sympy.pi)
    copy!(IM.__pyobject__, sympy.I)
    copy!(oo.__pyobject__, sympy.oo)
    copy!(zoo.__pyobject__, sympy.zoo)
    copy!(True.__pyobject__, PyCall.PyObject(true))
    copy!(False.__pyobject__, PyCall.PyObject(false))

    # mpmath
    try
        PyCall.mpmath_init()
        copy!(mpmath, PyCall.pyimport_conda("mpmath", "mpmath"))
    catch err
       # can't load
    end

    # pull in alibrary
    copy!(combinatorics, PyCall.pyimport_conda("sympy.combinatorics", "sympy"))

    ## mappings from PyObjects to types.
    ## order here is fussy, as we needed ImmutableMatrix before Basic
    pytype_mapping(sympy_matrices.ImmutableMatrix, SymMatrix)
    pytype_mapping(sympy_matrices.ImmutableDenseMatrix, SymMatrix)

    basictype = sympy_core.basic.Basic
    pytype_mapping(basictype, Sym)

    pytype_mapping(sympy_matrices.Matrix, Array{Sym})
    pytype_mapping(sympy_matrices.MatrixBase, Array{Sym})

    pytype_mapping(sympy.FiniteSet, Set)

    pytype_mapping(sympy.combinatorics.permutations.Permutation, SymPermutation)
    pytype_mapping(combinatorics.perm_groups.PermutationGroup, SymPermutationGroup)

    if mpmath != PyCall.PyNULL()
        ## ignore warnings for now...
        mpftype = mpmath."mpf"
        pytype_mapping(mpftype, BigFloat) ## Raises warning!!!
        mpctype = mpmath."mpc"
        pytype_mapping(mpctype, Complex{BigFloat})
    end

    ## We *would* like to import some methods from the `sympy` module, but we can't as calling
    ## @eval to import from a different package causes precompilation issues.
    ## Rather, we only call `import_sympy` to print a list of commands to add to this file (as
    ## are added below). This is done by uncommenting the `println` statements in `import_from`.
    #import_sympy()

end

## On load, generic functions in the base modules that match a function in sympy.* are
## defined on `Sym` objects, as are those listed here:
##
priviledged = (:And, :Or, :Not, :Xor,
               #
               :apart, :cancel, :cse, :expand, :factor, :flatten, :nsimplify,
               :isolate, :simplify, :together, :unflatten,
               #
               :srepr,:doit,
               #
               :integrate, :line_integrate, :interpolate, :limit, :series, :summation,
               :hessian,
               #
               :prime, :multiplicity, :degree, :coeffs,
               #
               :DiracDelta, :Heaviside,
               #
               :linsolve, :nonlinsolve, :nroots, :nsolve, :pdsolve, :real_root,
               :real_roots, :root, :rootof, :roots, :rsolve, :solve, :solveset,
               :ode_order,
               #
               :Min, :Max, :Abs,:numer, :denom, :conjugate, :ln,
               #
               :intersection, :intervals, :isprime

               )

function in_base(uv)
    u = Symbol(uv[1])
    for M in base_Ms
        isdefined(M, u) && return true
    end
    false
end
is_function(uv) = SymPy.pycall_hasproperty(uv[2], :__class__) &&  occursin("unction", uv[2].__class__.__name__)

"""
    import_sympy

This method imports all functions from `mpmath` and a priviledged set
of functions from `sympy`, as well as the relational operators.

These functions are narrowed on their first argument being of type `SymbolicObject`.

A few modules are checked for namespace collisions.

If a function naturally takes an non-Symbolic argument as a first argument, then it can be qualified: e.g.
`sympy.sin(2)` (as opposed to `sin(Sym(2))`).

If a constructor is needed (which is not a function) then it must be
qualified. (E.g. `sympy.Function("F")`, though for this particular
case, there is `SymFunction` defined for convenience.)

"""
function import_sympy()
    if mpmath != PyCall.PyNULL()
        mps = _get_member_functions(mpmath, base_exclude)
        sps = Introspection.getmembers(sympy)
        imp = Tuple(Symbol.(setdiff(keys(mps), keys(sps))))
        import_from(mpmath, imp)
    end
    ## import from
    ## import_from(sympy)
    d = Introspection.getmembers(sympy)
    d1 = filter(uv -> in_base(uv) && is_function(uv), d)
    import_from(sympy, setdiff(Symbol.(collect(keys(d1))),  Symbol.(base_exclude)))
    import_from(sympy, priviledged)
    import_from(sympy, (:Ne,  :Le, :Eq, :Ge, :Gt,
                        :GreaterThan, :LessThan,
                        :StrictGreaterThan, :StrictLessThan,
                        :Equality, :Unequality
                        ), typ=:Number)
end

## :Lt is in Base.Order
import Base.Order: Lt
Lt(x::Number, args...;kwargs...) = sympy.Lt(x, args...; kwargs...)
export(Lt)

### Add generic methods and new methods
include("importexport.jl")

end # module


text:
"""
This function checks if value is in input dictionary.
"""
function check_input_for_variable(dict, name, default = 0)
  var = default
  if haskey(dict, name)
    var = dict[name]
  end
  return var
end

"""
This function calculates the adjusted capital cost of a project.
"""
function calculate_capital_costs(inputs)
  #capital costs are the cost of getting the plant to completion, minus incentive benefits, minus tax benefits
  plant_cost = calculate_plant_cost(inputs["cost_inputs"])
  incentives = calculate_incentive_benefits(inputs["incentive_inputs"], plant_cost)
  tax_benefits = calculate_tax_benefits(inputs["tax_inputs"], plant_cost, inputs["incentive_inputs"])
  return(plant_cost - incentives - tax_benefits)
end

"""
This function adds the total construction and development costs of a project.
"""
function calculate_plant_cost(cost_inputs)
  return sum(cost_inputs)
end

"""
This function returns the MACRS depreciation schedule of a project.
"""
function get_macrs_schedule(num_depreciation_years, macrs_df)
  return macrs_df[:, Symbol("Y"*string(num_depreciation_years))]
end

"""
This function calculates the incentive benefits which can come from grants and investment tax credits
   grants are in dollars/KW of capacity
   itc is in terms of percent of plant cost
   itc_recovery_factor scales the percent of capital costs which can be covered by the itc. This can be due to transaction frictions.
"""
function calculate_incentive_benefits(incentive_inputs, plant_cost)
  grants = check_input_for_variable(incentive_inputs, "grants", 0)
  itc = check_input_for_variable(incentive_inputs, "itc", 0)
  itc_recovery_factor = check_input_for_variable(incentive_inputs, "itc_recovery_factor", 1)
  return grants + itc*itc_recovery_factor*plant_cost
end
"""

This function calculates the tax benefits which are from depreciating the capital asset. tax inputs include:
    macrs:                            the depreciation schedule in percent for the asset
    depreciable_percent:              fraction of plant costs which can be depreciated
    discount_rate:                    used to transform future tax benefits to present value
    tax_rate:                         tax percent which is deducted.
"""
function calculate_tax_benefits(tax_inputs, plant_cost, incentive_inputs)
  #
  macrs = check_input_for_variable(tax_inputs, "macrs", [0])
  itc = check_input_for_variable(incentive_inputs, "itc", 0)
  depreciable_percent = check_input_for_variable(tax_inputs, "depreciable_percent", 1)
  tax_rate = check_input_for_variable(tax_inputs, "tax_rate", .25)
  discount_rate = check_input_for_variable(tax_inputs, "discount_rate", 0.08)
  #
  return calculate_present_value(plant_cost*depreciable_percent*(1-itc/2)*tax_rate*macrs, discount_rate)
end

"""
This function calculates the present value which takes yearly values and a discount rate and returns the present value.
"""
function calculate_present_value(values, discount_rate)
  return sum([values[i]*(1+discount_rate)^(-i) for i in (1:length(values))])
end

"""
This function calculates calculates the weighted average cost of capital
 Takes a list with: Tax rate, debt interest rate, debt fraction, and either the equity discount rate or the required values to calculate a CAPM model
"""
function calculate_wacc(wacc_inputs)
  tax_rate =  check_input_for_variable(wacc_inputs, "tax_rate", 0.25)
  debt_fraction =  check_input_for_variable(wacc_inputs, "debt_fraction", 0.6)
  debt_interest_rate =  check_input_for_variable(wacc_inputs, "debt_interest_rate", 0.05)
  equity_discount_rate =  check_input_for_variable(wacc_inputs, "equity_discount_rate", 0.12)
  #
  debt_contribution = (1- tax_rate)*debt_interest_rate*debt_fraction
  equity_contribution = equity_discount_rate*(1-debt_fraction)
  return (debt_contribution + equity_contribution)
end
#_________________________________________________________________________________________________________________________________________________
"""
This function calculates a CAPM model equity rate = risk-free rate + beta*(market return - risk_free rate) + project-specific adder
"""
function calculate_equity_discount_rate(inputs)
  project_specific_adder = check_input_for_variable(inputs, "project_specific_adder", 0)
  risk_free_rate = check_input_for_variable(inputs, "risk_free_rate", 0.04)
  beta = check_input_for_variable(inputs, "beta", 1)
  market_return = check_input_for_variable(inputs, "market_return", .12)
  return (risk_free_rate + beta*(market_return - risk_free_rate) + project_specific_adder)
end

"""
This function utilizes all the finance functions to return the WACC and adjusted capital cost of a project.
"""
function calculate_capcost_and_wacc(project_capex::Vector{Float64},
                                    category::String,
                                    investor_dir::String,
                                    online_year::Int64)
  finance_data =  read_data(joinpath(investor_dir, "finance_params.csv"))
  investor_characteristics = read_data(joinpath(investor_dir, "characteristics.csv"))

  project_row = findfirst(x-> x == category, finance_data[:, "Category"])
  macrs_df = read_data(joinpath(investor_dir, "MACRS Schedule.csv"))

  equity_discount_rate = finance_data[project_row, "EQUITY_RATE"]

  if online_year <= 2020
      itc = finance_data[project_row, "ITC 2020"]
  elseif online_year == 2021
      itc = finance_data[project_row, "ITC 2021"]
  else
      itc = finance_data[project_row, "ITC post 2021"]
  end

  incentive_inputs = Dict("itc" => itc,
                          "grants" => finance_data[project_row, "GRANTS"],
                          "itc_recovery_factor" => finance_data[project_row, "ITC_RECOVERY_FACTOR"])

  macrs_schedule = get_macrs_schedule(finance_data[project_row, "MACRS_DEPRECIATION_YEARS"], macrs_df)
  wacc_inputs = Dict("tax_rate" => finance_data[project_row, "TAX_RATE"],
                      "debt_fration" => finance_data[project_row, "DEBT_FRACTION"],
                      "debt_interest_rate" => finance_data[project_row, "DEBT_INTEREST_RATE"],
                      "equity_discount_rate" => finance_data[project_row, "EQUITY_RATE"])

  discount_rate = calculate_wacc(wacc_inputs) + investor_characteristics[1, "discount_rate_adder"]

  tax_inputs = Dict("macrs" => macrs_schedule,
                    "tax_rate" => finance_data[project_row, "TAX_RATE"],
                    "discount_rate" => discount_rate,
                    "depreciable_percent"=> finance_data[project_row, "DEPRECIABLE_PERCENT"])

  cap_cost = zeros(length(project_capex))
  for i in 1:length(project_capex)
      COSTS = [project_capex[i]]
      capital_cost_inputs = Dict("cost_inputs" => [project_capex[i]], "incentive_inputs" => incentive_inputs, "tax_inputs" => tax_inputs)
      cap_cost[i] = calculate_capital_costs(capital_cost_inputs)
  end

  return cap_cost, discount_rate
end


text:
__precompile__()

module OperatingSystems

export
    CURRENT_OS, # Could be Sys.OS if this goes into Base

    # Supertypes
    AbstractOS,
    Unix,
    BSD,

    # Specific systems
    DragonFly,
    FreeBSD,
    Fuchsia,
    Linux,
    MacOS,
    NetBSD,
    OpenBSD,
    Windows,

    # Unknown system
    UnknownOS

abstract type AbstractOS end
abstract type Unix <: AbstractOS end
abstract type BSD <: Unix end

struct DragonFly <: BSD end
struct FreeBSD <: BSD end
struct Fuchsia <: AbstractOS end # lol
struct Linux <: Unix end
struct MacOS <: BSD end
struct NetBSD <: BSD end
struct OpenBSD <: BSD end
struct UnknownOS <: AbstractOS end
struct Windows <: AbstractOS end

const CURRENT_OS = let kern = Sys.KERNEL
    if kern === :Linux
        Linux
    elseif kern === :Windows || kern === :NT
        Windows
    elseif kern === :Darwin || kern === :Apple
        MacOS
    elseif kern === :FreeBSD || kern === :OpenBSD || kern === :NetBSD || kern === :DragonFly
        getfield(current_module(), kern)::BSD
    elseif kern === :Magenta
        Fuchsia
    else
        UnknownOS
    end
end

end # module


text:
@rule typeof(dot)(:in1, Marginalisation) (m_out::UnivariateNormalDistributionsFamily, m_in2::PointMass, meta::AbstractCorrection) = begin
    return @call_rule typeof(dot)(:in2, Marginalisation) (m_out = m_out, m_in1 = m_in2, meta = meta)
end


text:
using PsychoJulia
using GeometryTypes
using Random
using Test

@testset "States" begin
    fstate = PsychoJulia.FixationState(Rectangle(-10,-10,20,20),
                                       PsychoJulia.ClockTimer(time(), 5.0, false),Point2f0(0.0, 0.0))
    b =  PsychoJulia.is_inside(fstate, Point2f0(0.0, 0.0))
    @test b

    rstate = PsychoJulia.ResponseState(Rectangle(50.0,50.0, 20.0, 20.0),
                                       PsychoJulia.ClockTimer(time(), 5.0, false), Point2f0(0.0, 0.0))
    b = PsychoJulia.is_inside(rstate, Point2f0(51.0, 51.0))
    @test b

    estate = PsychoJulia.TrialEndState(PsychoJulia.ClockTimer(time(), 5.0, false))
    trans = PsychoJulia.StateTransitions([fstate, rstate, estate], [2 3 1;3 3 3], 1)
    done = false
    record = PsychoJulia.ExperimentRecord()
    ntrials = 0
    PsychoJulia.reset!(fstate.clock)
    while !done
        state = trans.states[trans.current]
        aa = PsychoJulia.get_area(state)
        PsychoJulia.set_position!(state,Point2f0(aa.x + 1.0, aa.y+1.0))
        PsychoJulia.next!(trans, record)
        ntrials = 0
        for ee in record.event
            if occursin("-> 3", ee)
                ntrials += 1
            end
        end
        done = ntrials >= 2
    end
    @test record.event == ["1 -> 2", "2 -> 3", "3 -> 1", "1 -> 2", "2 -> 3"]

    dd = diff(record.timestamp)
    #check that we get the correct timestamps to whitin 1ms precision
    @test dd ≈ [5.0, 5.0, 5.0, 5.0] atol=0.001
end


text:
using Test

#
# These tests mostly try to replicate examples from
# World Bank 'Handbook on poverty and inequality'
# by Haughton and Khandker
# http://documents.worldbank.org/curated/en/488081468157174849/Handbook-on-poverty-and-inequality
#

# ...
primitive type Amount <: AbstractFloat 64 end

const  TOL = 0.00001; # for nearly equal

"
This just creates an array which is `times` vcat copies of `a`
"
function vcn( a :: Array{Float64}, times :: Int64 )
    nrows = size( a )[1]
    ncols = size( a )[2]
    newrows = nrows*times
    out = zeros( Float64, newrows, ncols )
    p = 0
    for row in 1:nrows
        for i in 1:times
            p += 1
            out[p,:] .= a[row,:]
        end
    end
    out
end

NONE=Array{Symbol,1}();

function removeIgnored( d :: Dict{ Symbol, <:Any}, ignore::Array{Symbol,1} ):: Dict{ Symbol, <:Any}
    for i in ignore
        if haskey( d, i )
            delete!( d, i )
        end
    end
    return d;
end

"
element - by element compare of the type of dicts we use for poverty and inequality output
"
function comparedics( left :: Dict{ Symbol, <:Any}, right :: Dict{ Symbol, <:Any}, ignore::Array{Symbol,1} = NONE ) :: Bool
    left  = removeIgnored( left, ignore )
    right = removeIgnored( right, ignore )
    lk = keys( left )
    if lk != keys( right )
        return false
    end
    for k in lk
        # try catch here in case types are way off
        try
            if !( left[k] ≈ right[k] )
                l = left[k]
                r = right[k]
                print( "comparison failed '$k' : left = $l right = $r")
                return false
            end
       catch e
            return false
       end
    end
    return true
end

@testset "WB Chapter 4 - Poverty " begin
    country_a = [ 1.0 100; 1.0 100; 1 150; 1 150 ]
    country_b = copy( country_a )
    country_b[1:2,2] .= 124
    country_c = [ 1.0 100; 1.0 110; 1 150; 1 160 ]
    #
    # Ch 4 doesn't discuss weighting issues, so
    # we'll add some simple checks for that.
    # a_2 and b_2 should be the same as country_a and _b,
    # but with 2 obs of weight 2 rather than 4 of weight 1
    #
    country_a_2 = [2.0 100; 2.0 150 ]
    country_b_2 = [2.0 124; 2.0 150 ]
    # d should be a big version of a and also produce same result
    country_d = vcn( country_a, 50 )
    # attempt to blow things up with huge a clone
    country_d = vcn( country_c, 100_000 )

    # very unbalanced copy of dataset 1 with 10,000 weight1 1:2 and 2 weight 10,000 7:10
    country_e = vcn( country_c[1:2,:], 10_000 )
    cx = copy(country_c[3:4,:])
    cx[:,1] .= 10_000
    country_e = vcat( country_e, cx )

    line = 125.0
    growth = 0.05

    country_a_pov = makepoverty( country_a, line, growth )
    print("country A " );println( country_a_pov )
    country_a_2_pov = makepoverty( country_a_2, line,growth )
    country_b_pov = makepoverty( country_b, line, growth )
    country_c_pov = makepoverty( country_c, line, growth )
    print("country C " );println( country_c_pov )
    country_d_pov = makepoverty( country_d, line, growth )
    print("country D " );println( country_d_pov )
    country_e_pov = makepoverty( country_e, line, growth )
    print("country E " );println( country_e_pov )

    @test comparedics( country_a_pov, country_a_2_pov )
    @test comparedics( country_c_pov, country_e_pov )
    # @test comparedics( country_c_pov, country_d_pov )

    # numbers from WP ch. 4
    @test country_a_pov[:headcount] ≈ 0.5
    @test country_b_pov[:headcount] ≈ 0.5
    @test country_b_pov[:gap] ≈ 1.0/250.0
    @test country_c_pov[:watts] ≈ 0.0877442307
    # some of these are hand-calculations, for from Ada version
    @test isapprox( country_c_pov[:gap], 0.080000, atol = TOL )
    @test isapprox( country_c_pov[:foster_greer_thorndyke][ 1 ], 0.5000000, atol = TOL ) # pov level
    @test isapprox( country_c_pov[:foster_greer_thorndyke][ 2 ], 0.1984059, atol = TOL )
    @test isapprox( country_c_pov[:foster_greer_thorndyke][ 3 ], 0.0800000, atol = TOL )
    @test isapprox( country_c_pov[:foster_greer_thorndyke][ 4 ], 0.0327530, atol = TOL )
    @test isapprox( country_c_pov[:foster_greer_thorndyke][ 5 ], 0.0136000, atol = TOL )
    @test isapprox( country_c_pov[:foster_greer_thorndyke][ 6 ], 0.0057192, atol = TOL )
    @test isapprox( country_c_pov[:sen ], 0.0900000, atol = TOL )
    @test isapprox( country_c_pov[:shorrocks ], 0.0625000, atol = TOL )
    @test isapprox( country_c_pov[:watts ], 0.0877442, atol = TOL )
    @test isapprox( country_c_pov[:time_to_exit ], 1.7548846, atol = TOL )
    @test isapprox( country_c_pov[:gini_amongst_poor ], 0.0238095, atol = TOL )
    @test isapprox( country_c_pov[:poverty_gap_gini ], 0.5625000, atol = TOL )
end # poverty testset

#
# reproduce WB Table 6.3 with various cominations ofweights & data
# note table has errors:
#
#    1. uses log10 not ln for  theil##
#    2. has N in wrong place for ge(2) - outside bracket
#
@testset "WB Chapter 6 - Inequality " begin
    c1 = [1.0 10; 1 15; 1 20; 1 25; 1 40; 1 20; 1 30; 1 35; 1 45; 1 90 ]
    # these next are copies of c1 intended
    # to check we haven't screwed up the weighting
    c2 = vcn( c1, 2 )
    c3 = copy( c1 )
    c3[:,1] .= 10_000.0
    c4 = copy( c1 )
    c4[:,1] .= 2.0
    # very unbalanced copy of dataset 1 with 100,000 weight1 1:6 and 4 weight 100,000 7:10
    c64k = vcn( c1[1:6,:], 100_000 )
    cx = copy(c1[7:10,:])
    cx[:,1] .= 100_000
    c64k = vcat( c64k, cx )
    iq1 = makeinequality( c1 )
    iq2 = makeinequality( c2 )
    iq3 = makeinequality( c3 )
    iq4 = makeinequality( c4 )
    iq64k = makeinequality( c64k )
    # weighting and multiplying should make no difference
    @test comparedics( iq1, iq2, [:total_income, :total_population, :deciles] )
    println( "iq1");println( iq1 )
    println( "iq2");println( iq2 )
    println( "iq3");println( iq3 )
    println( "iq64k");println( iq64k )

    @test comparedics( iq1, iq3, [:total_income, :total_population, :deciles] )
    @test comparedics( iq1, iq4, [:total_income, :total_population, :deciles] )
    @test comparedics( iq1, iq64k, [:total_income, :total_population, :deciles] )
    @test isapprox( iq1[:gini ], 0.3272727, atol = TOL )
    @test isapprox( iq1[:theil][1], 0.1792203, atol = TOL )
    @test isapprox( iq1[:theil][2],  0.1830644, atol = TOL )
    @test isapprox( iq1[:generalised_entropy][ 1 ], 0.1883288, atol = TOL )
    @test isapprox( iq1[:generalised_entropy][ 2 ], 0.1954897, atol = TOL )
    @test isapprox( iq1[:generalised_entropy][ 3 ], 0.2047211, atol = TOL )
    @test isapprox( iq1[:generalised_entropy][ 4 ], 0.2162534, atol = TOL )
    @test isapprox( iq1[:generalised_entropy][ 5 ], 0.2303812, atol = TOL )
    @test isapprox( iq1[:generalised_entropy][ 6 ], 0.2474728, atol = TOL )
    @test isapprox( iq1[:atkinson][ 1 ], 0.0446396, atol = TOL )
    @test isapprox( iq1[:atkinson][ 2 ], 0.0869155, atol = TOL )
    @test isapprox( iq1[:atkinson][ 3 ], 0.1267328, atol = TOL )
    @test isapprox( iq1[:atkinson][ 4 ], 0.1640783, atol = TOL )
    @test isapprox( iq1[:atkinson][ 5 ], 0.1989991, atol = TOL )
    @test isapprox( iq1[:atkinson][ 6 ], 0.2315817, atol = TOL )
    @test isapprox( iq1[:atkinson][ 7 ], 0.2619332, atol = TOL )
    @test isapprox( iq1[:atkinson][ 8 ], 0.2901688, atol = TOL )
    @test isapprox( iq1[:atkinson][ 9 ], 0.3164032, atol = TOL )
    @test isapprox( iq1[:hoover], 0.2363636, atol = TOL )
    print( iq1 )
end # inequality testset


text:
# function decon(method::DeconvolutionMethod, X::Vector, σ_X::Vector)
    # ϵ_distr = Normal.(0.0, σ_X)
    # deconvolved = decon(method, X, ϵ_distr)
    # return deconvolved
# end

# """ Sample from the CDF """
# function cdf_sample(xx, F)
    # u = rand() # random uniform
    # i = searchsortedfirst(F, u)
    # return xx[i]
# end


text:
using Revise

using DataFrames, Statistics, Distributions, StatsBase,
  LinearAlgebra, BenchmarkTools, Random, Optim, Zygote, UnPack

MFloat64 = Union{Missing, Float64}

Random.seed!(11)

# a simple AR(1) model
struct AR1{T<:Real, Tσ²<:Union{Nothing,Real}}
  ϕ::T
  b::T
  μ::T
  σ²::Tσ² #note- this is σ²e, not the overall σ² of the series
end

AR1(;ϕ, μ=missing, b=missing, σ²=nothing) = AR1(
  ϕ, coalesce(b, μ*(1-ϕ)), coalesce(μ, b/(1-ϕ)),σ²)

centerar1(ar1::AR1) = AR1(ar1.ϕ, 0.0, 0.0, ar1.σ²)

function estimatear1ols(x::Vector{T}) where T
  x̃ = @view x[2:end]
  Lx̃ = @view x[1:(end-1)]

  ϕ = cov(x̃, Lx̃)/(var(Lx̃))
  b = mean(x̃) - ϕ * mean(Lx̃)
  μ = b/(1-ϕ)

  ε = similar(x)
  ε[1] = x[1] - μ
  ε[2:end] = x̃ .- ϕ .* Lx̃
  σ² = var(ε)

  @assert b ≈ mean(x̃) - ϕ * mean(Lx̃)
  @assert ϕ ≈ cov(x̃ .- mean(x̃), Lx̃ .- mean(Lx̃))/var(
    Lx̃ .- mean(Lx̃))

  return AR1(;ϕ,b,σ²)
end

function estimatear1mle(x; tol=10^-10)

  #start by centering
  N=length(x)
  μ = mean(x)
  x̃ = @view(x[2:end]) .- μ
  Lx̃ = @view(x[1:(end-1)]) .- μ

  x₁ = x[1]
  #computes the next iteration of ϕ given a previous ϕ
  cσ²(ϕ) = (x₁^2*(1-ϕ^2) + sum((x̃ .- ϕ*Lx̃).^2))/N
  function objϕ(ϕ)
    σ² = cσ²(ϕ)
    return (sum((x̃ .- ϕ * Lx̃) .* Lx̃) - (σ²/(1-ϕ^2) - x₁^2)*ϕ)^2
  end

  ϕ = optimize(objϕ, -1.01, 1.01, Brent()) |> Optim.minimizer

  return AR1(;ϕ, μ, σ² = cσ²(ϕ))
end

function estimatear1mlecheck(x; tol=10^-10)

  #start by centering
  N=length(x)
  μ = mean(x)
  x̃ = @view(x[2:end]) .- μ
  Lx̃ = @view(x[1:(end-1)]) .- μ

  x₁ = x[1]

  l(ϕ, σ²) = (-N/2*log(2*π*σ²)
    + 1/2*log(1-ϕ^2)
    - x₁^2*(1-ϕ^2)/(2*σ²)
    - sum((x̃ .- ϕ .* Lx̃).^2)/(2*σ²)) * -1

  l(Φ) = l(Φ[1], Φ[2])

  ϕ₀ = cov(x̃, Lx̃)/var(Lx̃) #start by guessing the least squares solution
  ε₀= x̃ - ϕ₀*Lx̃
  σ²₀ = var(ε₀)

  ϕ, σ² = optimize(l,  [-1.0,0.0], [1.0,Inf], [ϕ₀, σ²₀], Fminbox(LBFGS()),
    autodiff = :forward) |> Optim.minimizer

  return AR1(;ϕ, μ, σ²)
end

AR1(x::AbstractVector; estimationmethod=estimatear1mle) = estimationmethod(x)

function checkmath2(x)
  if length(x) > 250
    x = x[1:250] |> deepcopy
    @warn "length(x) > 250 so truncating length to 250. This is only for the purposes
    of computing the full (not log) likelihood"
  end

  N=length(x)
  μ = mean(x)
  x̃ = @view(x[2:end]) .- μ
  Lx̃ = @view(x[1:(end-1)]) .- μ

  ϕ₀ = 0.75
  σ²₀ = 2.2
  x₁ = x̃[1]

  L(ϕ, σ²) =  (1-ϕ^2)^(1/2)/σ²^(1/2)*pdf(Normal(), x₁/(σ²/(1-ϕ^2))^(1/2))*prod(
    ((x̃ᵢ,Lx̃ᵢ)->1/σ²^(1/2)*pdf(Normal(), (x̃ᵢ-ϕ*Lx̃ᵢ)/σ²^(1/2))).(x̃,Lx̃))

  l(ϕ, σ²) = (-N/2*log(2*π*σ²)
    + 1/2*log(1-ϕ^2)
    - x₁^2*(1-ϕ^2)/(2*σ²)
    - sum((x̃ .- ϕ .* Lx̃).^2)/(2*σ²))

  @assert log(L(ϕ₀, σ²₀)) ≈ l(ϕ₀, σ²₀,) "log(L(ϕ, σ²)): $(log(L(ϕ₀, σ²₀)))
    l(ϕ, σ²): $(l(ϕ₀, σ²₀))"

  #now verify the derivitive and conditional wrt sigma^2
  lσ²(ϕ=ϕ₀, σ²=σ²₀,) = (-N/(2σ²)
    + x₁^2*(1-ϕ^2)/(2σ²^2)
    + sum((x̃ .- ϕ .* Lx̃).^2)/(2σ²^2))
  lσ²z(ϕ=ϕ₀, σ²=σ²₀,) = gradient(σ²->l(ϕ, σ²,), σ²)[1]
  lσ²n() =  (l(ϕ₀, σ²₀+1e-6,)-l(ϕ₀, σ²₀,))/1e-6
  @assert lσ²() ≈ lσ²z() "lσ²(...): $(lσ²())
    lσ²z(...): $(lσ²z())"

  cσ²(ϕ) = (x₁^2*(1-ϕ^2) + sum((x̃ .- ϕ*Lx̃).^2))/N
  @assert 1 + lσ²z(ϕ₀, cσ²(ϕ₀)) ≈ 1.0 "cσ²(ϕ₀): $(cσ²(ϕ₀))
    lσ²z(ϕ₀, cσ²(ϕ₀)): $(lσ²z(ϕ₀, cσ²(ϕ₀)))"

  #finally, verify the derivitive wrt ϕ
  lϕ(ϕ=ϕ₀, σ²=σ²₀) = (-ϕ/(1-ϕ^2) +
    x₁^2*ϕ/σ² +
    sum((x̃ .- ϕ * Lx̃) .* Lx̃)/σ²)
  lϕz(ϕ=ϕ₀, σ²=σ²₀,) = gradient(ϕ->l(ϕ, σ²), ϕ)[1]
  @assert lϕ() ≈ lϕz() "lϕ(...): $(lϕ())
    lϕz(...): $(lϕz())"

  ϕcerr(ϕ, σ²=σ²₀)= sum((x̃ .- ϕ * Lx̃) .* Lx̃) - (σ²/(1-ϕ^2) - x₁^2)*ϕ
  cϕ(σ²=σ²₀,) = optimize(ϕ->ϕcerr(ϕ)^2, -1.01, 1.01, Brent(), abs_tol=10^-14) |> Optim.minimizer
  @assert isapprox(lϕz(cϕ()), 0.0,atol=1e-5) "cϕ(ϕ₀): 2$(cϕ())
    lϕz(ϕcerr()): $(lϕz(ϕcerr()))"
  #lμ()
end

function simulatear1!(Xt::AbstractMatrix;
    ϕ::TV, σ::TV, μ::TV) where TV<:AbstractVector
  #@unpack ϕ, b, σ², μ = ar1

  #first simulate the centered distribution
  #we will do this all in place

  rand!(Normal(), Xt)
  Xt .*= σ

  Xt[:,1] .*= 1 ./ (1 .- ϕ.^2) .^ 0.5 # the first entry is simulated from the unconditional distribution

  #=Threads.@threads for c ∈ 1:size(X,2)
    for t ∈ 2:size(X,1)
      X[t,c] += ϕ[c] * X[t-1,c]
    end
  end=#

  for t ∈ 2:size(Xt,2)
    for c ∈ 1:size(Xt,1)
      Xt[c,t] += ϕ[c] * Xt[c,t-1]
    end
  end

  Xt .+= μ

  return nothing
end

#simulates an AR1 series
function simulatear1!(x::AbstractVector; ϕ, σ, μ)

  rand!(Normal(0,σ), x)
  x[1] *= 1/(1-ϕ^2)^0.5 # the first entry is simulated from the unconditional distribution
  for t ∈ 2:length(x)
    x[t] += ϕ * x[t-1]
  end
  x .+= μ

  return nothing
end

#non-inplace version for testing
function simulatear1(ar1::AR1{T,Tσ²}, N::Int) where {T,Tσ²}
  res = Vector{Float64}(undef, N)
  simulatear1!(res; σ=ar1.σ²^0.5, μ=ar1.μ, ϕ=ar1.ϕ)
  return res
end

function compareoncear1xy!(Xt::AbstractMatrix,
  yar1::AR1,  y::AbstractVector,
  comparefunc::Tcomparefunc;
  ϕx::Vector{Float64},σx::Vector{Float64}, μx::Vector{Float64}) where Tcomparefunc

  simulatear1!(Xt; ϕ=ϕx, σ=σx, μ=μx)
  simulatear1!(y, ϕ=yar1.ϕ, σ=yar1.σ²^0.5, μ=yar1.μ)
  comparefunc(Xt',y)
end

function compareoncear1xy!(Xt::AbstractVector,
  yar1::AR1,  y::AbstractVector,
  comparefunc::Tcomparefunc;
  ϕx::Vector{Float64},σx::Vector{Float64}, μx::Vector{Float64}) where Tcomparefunc

  simulatear1!(Xt; ϕ=ϕx[1], σ=σx[1], μ=μx[1])
  simulatear1!(y, ϕ=yar1.ϕ, σ=yar1.σ²^0.5, μ=yar1.μ)
  comparefunc(Xt,y)
end

#this runs the simulations- basically a parametric bootstrap under the null
#statfunc is computed for each pass, while aggfunc is computed on the statistics
#at the end of the pass
function comparear1undernull(X::Tx,y::AbstractVector,
    statfunc::Tstatfunc, aggfunc::Taggfunc; Nsimulations::Int,
    checkstationarity::Bool = true, stationaritytol::Float64=10^-4,
    Xts = (size(X,2) > 1 ?
      [similar(X' |> Matrix{Float64}) for t ∈ 1:Threads.nthreads()] :
      [similar(X |> vec |> Vector{Float64}) for t ∈ 1:Threads.nthreads()])
    ) where {Tx<:AbstractMatrix, Tstatfunc, Taggfunc}

  x = collect(eachcol(X))
  K = length(x)
  xar1 = (AR1.(x)) .|> centerar1
  yar1 = AR1(y) |> centerar1

  results = Dict(s => missings(Float64, Nsimulations, K) for s ∈ keys(statfunc(X,y)))
  @assert (k->ismissing.(results[k]) |> all).(keys(results)) |> all

  #this assumes stationary results are meaningless
  if checkstationarity && max((xar1 .|> xar1i -> abs(xar1i.ϕ))..., abs(yar1.ϕ)) ≥ 1.0-stationaritytol
    xar1 .= ((xar1i)->AR1(NaN, NaN, xar1i.μ, xar1i.σ²)).(xar1)
    yar1 = AR1(NaN, NaN, yar1.μ, yar1.σ²)
  end

  ys = [similar(y) for t ∈ 1:Threads.nthreads()]

  ϕx::Vector{Float64} = xar1 .|> xar1i->xar1i.ϕ
  σx::Vector{Float64} = xar1 .|> xar1i->xar1i.σ²^0.5
  μx::Vector{Float64} = xar1 .|> xar1i->xar1i.μ

  #draw a sample and compute the aggregate function, sotring the reuslts in an array
  Threads.@threads for n ∈ 1:Nsimulations
  #for n ∈ 1:Nsimulations

    t::Int = Threads.threadid()
    #res = compareoncear1xy!(xar1, xs[1], yar1, ys[1], statfunc)
    res = compareoncear1xy!(Xts[t], yar1, ys[t], statfunc; ϕx, σx, μx)

    #this structure allows us to compute multiple statistics
    #println(res)
    for s ∈ keys(res)

      @assert ismissing.(results[s][n, :]) |> all
      results[s][n,:] .= res[s]
    end
  end

  return (; Dict(s=>aggfunc(results[s], s) for s ∈ keys(results))...)
end

#this is just a wrapper for the 2-vector case
function compareundernull(x::Tx,y::AbstractVector,
    statfunc::Tstatfunc, aggfunc::Taggfunc, comparefunc::Tcomparefunc; Nsimulations::Int,
    checkstationarity::Bool = true, stationaritytol::Float64=10^-4
    ) where {Tx<:AbstractVector, Tstatfunc, Taggfunc, Tcomparefunc}

  #=function vstatfunc(X,y)
    #println("X $(size(X))")
    rawstats = statfunc(X |> vec, y)
    return Dict(k=>rawstats[k] for (k,v) ∈ zip(keys(rawstats), values(rawstats)))
  end
  vaggfunc(vr,s) = aggfunc(vr |> vec, s)=#

  return comparefunc(reshape(x,length(x),1) |> Matrix{Float64}, y, statfunc, aggfunc;
    Nsimulations, checkstationarity, stationaritytol)
end

#wrapper for the wrapper
comparear1undernull(x::Tx,y::AbstractVector,
    statfunc::Tstatfunc, aggfunc::Taggfunc; Nsimulations::Int,
    checkstationarity::Bool = true, stationaritytol::Float64=10^-4
    ) where {Tx<:AbstractVector, Tstatfunc, Taggfunc} = compareundernull(
      x,y,statfunc,aggfunc, comparear1undernull;
      Nsimulations, checkstationarity, stationaritytol)

#=function comparearima11undernull_old(x::Tx,y::AbstractVector,
    statfunc::Tstatfunc, aggfunc::Taggfunc; Nsimulations::Int,
    checkstationarity::Bool = true, stationaritytol::Float64=10^-4) where {Tx, Tstatfunc, Taggfunc}

    #difference everything
    dx = x[2:end] .- x[1:(end-1)]
    dy = x[2:end] .- y[1:(end-1)]

    dxar1 = AR1(dx) |> centerar1
    dyar1 = AR1(dy) |> centerar1

    results = Dict(s => Vector{MFloat64}(undef, Nsimulations) for s ∈ keys(statfunc(dx,dy)))

    #this assumes stationary results are meaningless
    if checkstationarity && max(abs(dxar1.ϕ), abs(dyar1.ϕ)) ≥ 1.0-stationaritytol
      dxar1 = AR1(NaN, NaN, dxar1.μ, dxar1.σ²)
      dyar1 = AR1(NaN, NaN, dxar1.μ, dxar1.σ²)
    end

    dxs = [similar(dx) for t ∈ 1:Threads.nthreads()]
    dys = [similar(dy) for t ∈ 1:Threads.nthreads()]

    #draw a sample and compute the aggregate function, sotring the reuslts in an array
    Threads.@threads for n ∈ 1:Nsimulations
      t::Int = Threads.threadid()
      res = compareoncear1xy!(dxar1, dxs[t], dyar1, dys[t], statfunc)

      #this structure allows us to compute multiple statistics
      for s ∈ keys(res)
        @assert results[s][n] ≡ missing
        results[s][n] = res[s]
      end
    end

    return (; Dict(s=>aggfunc(skipmissing(results[s]), s) for s ∈ keys(results))...)
end=#

function comparearima11undernull(X::Tx,y::AbstractVector,
    statfunc::Tstatfunc, aggfunc::Taggfunc; Nsimulations::Int,
    checkstationarity::Bool = true, stationaritytol::Float64=10^-4
    ) where {Tx<:AbstractMatrix, Tstatfunc, Taggfunc}

  #difference everything
  dX = reduce(hcat, (xi->xi[2:end] .- xi[1:(end-1)]).(eachcol(X)))
  dx = dX |> eachcol |> collect
  dy = y[2:end] .- y[1:(end-1)]

  K = length(dx)

  dxar1 = (AR1.(dx)) .|> centerar1
  dyar1 = AR1(dy) |> centerar1

  results = Dict(s => missings(Float64, Nsimulations, K) for s ∈ keys(statfunc(dX,dy)))
  @assert (k->ismissing.(results[k]) |> all).(keys(results)) |> all

  #this assumes stationary results are meaningless
  if checkstationarity && max((dxar1i -> abs(dxar1i.ϕ)).(dxar1)..., abs(dyar1.ϕ)) ≥ 1.0-stationaritytol
    dxar1 .= ((dxar1i)->AR1(NaN, NaN, dxar1i.μ, dxar1i.σ²)).(dxar1)
    dyar1 = AR1(NaN, NaN, dyar1.μ, dyar1.σ²)
  end

  dXs = [similar(dX) for t ∈ 1:Threads.nthreads()]
  dxs = dXs .|> eachcol .|> collect
  dys = [similar(dy) for t ∈ 1:Threads.nthreads()]

  #draw a sample and compute the aggregate function, sotring the reuslts in an array
  Threads.@threads for n ∈ 1:Nsimulations
  #for n ∈ 1:Nsimulations

    t::Int = Threads.threadid()
    #res = compareoncear1xy!(xar1, xs[1], yar1, ys[1], statfunc)
    res = compareoncear1xy!(dxar1, dXs[t], dxs[t], dyar1, dys[t], statfunc)

    #this structure allows us to compute multiple statistics
    for s ∈ keys(res)

      @assert ismissing.(results[s][n, :]) |> all
      results[s][n,:] .= res[s]
    end
  end

  return (; Dict(s=>aggfunc(results[s], s) for s ∈ keys(results))...)
end

comparearima11undernull(x::Tx,y::AbstractVector,
    statfunc::Tstatfunc, aggfunc::Taggfunc; Nsimulations::Int,
    checkstationarity::Bool = true, stationaritytol::Float64=10^-4
    ) where {Tx<:AbstractVector, Tstatfunc, Taggfunc} = compareundernull(
      x,y,statfunc,aggfunc, comparearima11undernull;
      Nsimulations, checkstationarity, stationaritytol)

nonparametricnullp(v::AbstractVector{T}, x) where T<:Real = (
  x ≥ 0 ? sum(v .> x)/length(v) : sum(v .< x)/length(v))
nonparametricnullp(v::AbstractVector, x) = nonparametricnullp(v |> skipmissing |> collect, x)
nonparametricnullp(M::AbstractMatrix, x) = (
  v->nonparametricnullp(v |> skipmissing |> collect, x)).(eachcol(M))

function testar1(;N)
  ϕ = 0.8
  ε = rand(Normal(), N) + rand(Normal(1,1), N)
  x = missings(Float64, N)
  x[1] = ε[1]

  global x = x
  for t ∈ 2:N
    x[t] = ϕ*x[t-1] + ε[t]
  end
  @info "ϕ=$ϕ, b=$(mean(ε)), μ=$(mean(x))"

  checkmath2(x)

  #=estolsar1 = @btime AR1($x, estimationmethod=estimatear1ols)
  @info "ols conditional: $estolsar1"

  estmlear1 = @btime AR1($x, estimationmethod=estimatear1mle)
  @info "mle: $estmlear1"

  estmle2ar1 = @btime AR1($x, estimationmethod=estimatear1mlecheck)
  @info "mle2: $estmle2ar1"=#

  actar1=AR1(;ϕ, b=0.5, σ²=2)
  xsim = simulatear1(actar1, 10^6)
  @assert isapprox(mean(xsim), 2.5, atol=0.01)  "mean(xsim): $(mean(xsim))"
  @assert isapprox(var(xsim), 2/(1-ϕ^2), atol=0.01) "var(xsim): $(var(xsim))"
  ar1sim = AR1(xsim)
  @info "ar1sim: $(ar1sim)"
  @assert isapprox(ar1sim.ϕ, ϕ, atol=0.01)
  @assert isapprox(ar1sim.b, 0.5, atol=0.01)
  @assert isapprox(ar1sim.μ, 2.5, atol=0.01)
  @assert isapprox(ar1sim.σ², 2, atol=0.01)

  #=xsim! = simulatear1(actar1, Vector{Float64}(undef,10^6))
  @assert isapprox(mean(xsim), 2.5, atol=0.01)  "mean(xsim): $(mean(xsim))"
  @assert isapprox(var(xsim), 2/(1-ϕ^2), atol=0.01) "var(xsim): $(var(xsim))"
  ar1sim = AR1(xsim!)
  @info "ar1sim: $(ar1sim)"
  @assert isapprox(ar1sim.ϕ, ϕ, atol=0.01)
  @assert isapprox(ar1sim.b, 0.5, atol=0.01)
  @assert isapprox(ar1sim.μ, 2.5, atol=0.01)
  @assert isapprox(ar1sim.σ², 2, atol=0.01)=#

  Nsimulations=10^4
  dmeantest = 0.2
  #recall the compare func will pass the name of the relevant stat func when it is called
  aggfunc(v, s) = nonparametricnullp(v, Dict(:cor=>0.1, :dmean=>dmeantest)[s])
  statfunc(x,y) = Dict(:cor=>cor(x,y), :dmean=>mean(x)-mean(y))

  @info "std(x): $(std(x))"
  @info "expected std(x): $((2/(1-ϕ^2))^0.5)"
  #derive the distribution for the difference in means of two AR(1) dists
  R = [0.8^abs(i-j) for i ∈ 1:N, j ∈ 1:N]
  Sdmean = (2/N^2*(ones(N)'* ( R .* 2/(1-ϕ^2)) * ones(N)))^0.5
  Epdmean = cdf(Normal(), - dmeantest/Sdmean)
  res = comparear1undernull(x, deepcopy(x), statfunc, aggfunc; Nsimulations)
  @info "compareundernull result: $res (Epdmean: $Epdmean)"
  y = deepcopy(x)

  @btime comparear1undernull($x, $y, $statfunc, $aggfunc; Nsimulations=$Nsimulations)
  throw("simulation complete")

  X = hcat(x,x .* 2.0)
  mstatfunc(X,y) = Dict(:cor=>cor(X,y) |> vec, :cov=>cov(X,y) |> vec)
  maggfunc(M,s) = (v->nonparametricnullp(v, Dict(:cor=>0.1, :cov=>0.1)[s])).(eachcol(M))

  res = comparear1undernull(X, y, mstatfunc, maggfunc; Nsimulations)
  @info "compareundernull matrix result: $res"

  ix = cumsum(x .+ 0.2)
  iy = cumsum(x .+ 0.2)
  res = comparear1undernull(ix, deepcopy(ix), statfunc, aggfunc; Nsimulations)
  @info "integrated compareundernull (mispecified) dmean: $res"
  res = comparearima11undernull(ix, deepcopy(ix), statfunc, aggfunc; Nsimulations)
  @info "integrated compareundernull dmean: $res"
  res = comparear1undernull(X, deepcopy(x), mstatfunc, maggfunc; Nsimulations)
  @info "integrated compareundernull matrix result: $res"

end

@time testar1(N=1000)

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
"""
eyediag(wf, Tp, sps::Integer=10 ; traces = 100, n = 3, tini = 0.0)

Return the eye diagram of the waveform `wf`.

This function returns a matrix with `n*sps` rows and `traces` columns. Each column
corresponds to one trace of the waveform `wf`; each trace includes `n` symbols,
each of duration `Tp`. Each symbol is sampled `sps` times. The eye diagram is
evaluated starting at time `tini`.

To obtain a plot of the eye diagram, plot the columns of the returned matrix.
"""
function eyediag(wf, Tp, sps::Integer=10 ; traces = 100, n = 3, tini = 0.0)
    ed = zeros(ComplexF64, sps*n, traces)
    tend = tini + n*Tp
    for tr in 1:traces
        t = range(tini, tend, length=sps*n)
        ed[:,tr] = wf.(t)
        tini = tend
        tend = tini + n*Tp
    end
    return ed
end