text:
const mjtNum = Cdouble
const mjMINVAL = 1.0e-15

const mjPI = 3.14159265358979323846
const mjMAXVAL = 1.0e10
const mjMINMU = 1.0e-5
const mjMINIMP = 0.0001
const mjMAXIMP = 0.9999
const mjMAXCONPAIR = 50
const mjMAXVFS = 200
const mjMAXVFSNAME = 100

const mjNEQDATA = 7
const mjNDYN = 10
const mjNGAIN = 10
const mjNBIAS = 10
const mjNREF = 2
const mjNIMP = 5
const mjNSOLVER = 1000

const mjtByte = Cuchar

@cenum mjtDisableBit::Cint begin      # disable default feature bitflags
    mjDSBL_CONSTRAINT = 1 << 0         # entire constraint solver
    mjDSBL_EQUALITY = 1 << 1         # equality constraints
    mjDSBL_FRICTIONLOSS = 1 << 2         # joint and tendon frictionloss constraints
    mjDSBL_LIMIT = 1 << 3         # joint and tendon limit constraints
    mjDSBL_CONTACT = 1 << 4         # contact constraints
    mjDSBL_PASSIVE = 1 << 5         # passive forces
    mjDSBL_GRAVITY = 1 << 6         # gravitational forces
    mjDSBL_CLAMPCTRL = 1 << 7         # clamp control to specified range
    mjDSBL_WARMSTART = 1 << 8         # warmstart constraint solver
    mjDSBL_FILTERPARENT = 1 << 9         # remove collisions with parent body
    mjDSBL_ACTUATION = 1 << 10        # apply actuation forces
    mjDSBL_REFSAFE = 1 << 11        # integrator safety: make ref[0]>=2*timestep

    mjNDISABLE = 12           # number of disable flags
end

@cenum mjtEnableBit::Cint begin       # enable optional feature bitflags
    mjENBL_OVERRIDE = 1 << 0         # override contact parameters
    mjENBL_ENERGY = 1 << 1         # energy computation
    mjENBL_FWDINV = 1 << 2         # record solver statistics
    mjENBL_SENSORNOISE = 1 << 3         # add noise to sensor data

    mjNENABLE = 4          # number of enable flags
end

@cenum mjtJoint::Cint begin           # type of degree of freedom
    mjJNT_FREE = 0            # global position and orientation (quat)       (7)
    mjJNT_BALL                         # orientation (quat) relative to parent        (4)
    mjJNT_SLIDE                        # sliding distance along body-fixed axis       (1)
    mjJNT_HINGE                        # rotation angle (rad) around body-fixed axis  (1)
end

@cenum mjtGeom::Cint begin            # type of geometric shape
    # regular geom types
    mjGEOM_PLANE = 0            # plane
    mjGEOM_HFIELD                      # height field
    mjGEOM_SPHERE                      # sphere
    mjGEOM_CAPSULE                     # capsule
    mjGEOM_ELLIPSOID                   # ellipsoid
    mjGEOM_CYLINDER                    # cylinder
    mjGEOM_BOX                         # box
    mjGEOM_MESH                        # mesh

    mjNGEOMTYPES                       # number of regular geom types

                                     # rendering-only geom types: not used in mjModel not counted in mjNGEOMTYPES
    mjGEOM_ARROW = 100          # arrow
    mjGEOM_ARROW1                      # arrow without wedges
    mjGEOM_ARROW2                      # arrow in both directions
    mjGEOM_LINE                        # line
    mjGEOM_SKIN                        # skin
    mjGEOM_LABEL                       # text label

    mGEOM_NONE = 1001         # missing geom type
end

@cenum mjtCamLight::Cint begin        # tracking mode for camera and light
    mjCAMLIGHT_FIXED = 0            # pos and rot fixed in body
    mjCAMLIGHT_TRACK                   # pos tracks body, rot fixed in global
    mjCAMLIGHT_TRACKCOM                # pos tracks subtree com, rot fixed in body
    mjCAMLIGHT_TARGETBODY              # pos fixed in body, rot tracks target body
    mjCAMLIGHT_TARGETBODYCOM           # pos fixed in body rot tracks target subtree com
end

@cenum mjtTexture::Cint begin         # type of texture
    mjTEXTURE_2D = 0            # 2d texture, suitable for planes and hfields
    mjTEXTURE_CUBE                     # cube texture, suitable for all other geom types
    mjTEXTURE_SKYBOX                   # cube texture used as skybox
end

@cenum mjtIntegrator::Cint begin      # integrator mode
    mjINT_EULER = 0            # semi-implicit Euler
    mjINT_RK4                          # 4th-order Runge Kutta
end

@cenum mjtCollision::Cint begin       # collision mode for selecting geom pairs
    mjCOL_ALL = 0            # test precomputed and dynamic pairs
    mjCOL_PAIR                         # test predefined pairs only
    mjCOL_DYNAMIC                      # test dynamic pairs only
end

@cenum mjtCone::Cint begin            # type of friction cone
    mjCONE_PYRAMIDAL = 0           # pyramidal
    mjCONE_ELLIPTIC                    # elliptic
end

@cenum mjtJacobian::Cint begin        # type of constraint Jacobian
    mjJAC_DENSE = 0           # dense
    mjJAC_SPARSE                       # sparse
    mjJAC_AUTO                         # dense if nv<60 sparse otherwise
end

@cenum mjtSolver::Cint begin          # constraint solver algorithm
    mjSOL_PGS = 0           # PGS    (dual)
    mjSOL_CG                           # CG     (primal)
    mjSOL_NEWTON                       # Newton (primal)
end

@cenum mjtEq::Cint begin              # type of equality constraint
    mjEQ_CONNECT = 0            # connect two bodies at a point (ball joint)
    mjEQ_WELD                          # fix relative position and orientation of two bodies
    mjEQ_JOINT                         # couple the values of two scalar joints with cubic
    mjEQ_TENDON                        # couple the lengths of two tendons with cubic
    mjEQ_DISTANCE                      # fix the contact distance betweent two geoms
end

@cenum mjtWrap::Cint begin            # type of tendon wrap object
    mjWRAP_NONE = 0            # null object
    mjWRAP_JOINT                       # constant moment arm
    mjWRAP_PULLEY                      # pulley used to split tendon
    mjWRAP_SITE                        # pass through site
    mjWRAP_SPHERE                      # wrap around sphere
    mjWRAP_CYLINDER                    # wrap around (infinite) cylinder
end

@cenum mjtTrn::Cint begin             # type of actuator transmission
    mjTRN_JOINT = 0            # force on joint
    mjTRN_JOINTINPARENT                # force on joint, expressed in parent frame
    mjTRN_SLIDERCRANK                  # force via slider-crank linkage
    mjTRN_TENDON                       # force on tendon
    mjTRN_SITE                         # force on site

    mjTRN_UNDEFINED = 1000         # undefined transmission type
end

@cenum mjtDyn::Cint begin             # type of actuator dynamics
    mjDYN_NONE = 0            # no internal dynamics; ctrl specifies force
    mjDYN_INTEGRATOR                   # integrator: da/dt = u
    mjDYN_FILTER                       # linear filter: da/dt = (u-a) / tau
    mjDYN_MUSCLE                       # piece-wise linear filter with two time constants
    mjDYN_USER                         # user-defined dynamics type
end

@cenum mjtGain::Cint begin            # type of actuator gain
    mjGAIN_FIXED = 0            # fixed gain
    mjGAIN_MUSCLE                      # muscle FLV curve computed by mju_muscleGain()
    mjGAIN_USER                        # user-defined gain type
end

@cenum mjtBias::Cint begin            # type of actuator bias
    mjBIAS_NONE = 0            # no bias
    mjBIAS_AFFINE                      # const + kp*length + kv*velocity
    mjBIAS_MUSCLE                      # muscle passive force computed by mju_muscleBias()
    mjBIAS_USER                        # user-defined bias type
end

@cenum mjtObj::Cint begin             # type of MujoCo object
    mjOBJ_UNKNOWN = 0            # unknown object type
    mjOBJ_BODY                         # body
    mjOBJ_XBODY                        # body, used to access regular frame instead of i-frame
    mjOBJ_JOINT                        # joint
    mjOBJ_DOF                          # dof
    mjOBJ_GEOM                         # geom
    mjOBJ_SITE                         # site
    mjOBJ_CAMERA                       # camera
    mjOBJ_LIGHT                        # light
    mjOBJ_MESH                         # mesh
    mjOBJ_SKIN                         # skin
    mjOBJ_HFIELD                       # heightfield
    mjOBJ_TEXTURE                      # texture
    mjOBJ_MATERIAL                     # material for rendering
    mjOBJ_PAIR                         # geom pair to include
    mjOBJ_EXCLUDE                      # body pair to exclude
    mjOBJ_EQUALITY                     # equality constraint
    mjOBJ_TENDON                       # tendon
    mjOBJ_ACTUATOR                     # actuator
    mjOBJ_SENSOR                       # sensor
    mjOBJ_NUMERIC                      # numeric
    mjOBJ_TEXT                         # text
    mjOBJ_TUPLE                        # tuple
    mjOBJ_KEY                          # keyframe
end

@cenum mjtConstraint::Cint begin      # type of constraint
    mjCNSTR_EQUALITY = 0            # equality constraint
    mjCNSTR_FRICTION_DOF               # dof friction
    mjCNSTR_FRICTION_TENDON            # tendon friction
    mjCNSTR_LIMIT_JOINT                # joint limit
    mjCNSTR_LIMIT_TENDON               # tendon limit
    mjCNSTR_CONTACT_FRICTIONLESS       # frictionless contact
    mjCNSTR_CONTACT_PYRAMIDAL          # frictional contact, pyramidal friction cone
    mjCNSTR_CONTACT_ELLIPTIC           # frictional contact elliptic friction cone
end

@cenum mjtConstraintState::Cint begin # constraint state
    mjCNSTRSTATE_SATISFIED = 0         # constraint satisfied, zero cost (limit, contact)
    mjCNSTRSTATE_QUADRATIC             # quadratic cost (equality, friction, limit, contact)
    mjCNSTRSTATE_LINEARNEG             # linear cost, negative side (friction)
    mjCNSTRSTATE_LINEARPOS             # linear cost, positive side (friction)
    mjCNSTRSTATE_CONE                  # squared distance to cone cost (elliptic contact)
end

@cenum mjtSensor::Cint begin          # type of sensor
                                     # common robotic sensors attached to a site
    mjSENS_TOUCH = 0            # scalar contact normal forces summed over sensor zone
    mjSENS_ACCELEROMETER               # 3D linear acceleration, in local frame
    mjSENS_VELOCIMETER                 # 3D linear velocity, in local frame
    mjSENS_GYRO                        # 3D angular velocity, in local frame
    mjSENS_FORCE                       # 3D force between site's body and its parent body
    mjSENS_TORQUE                      # 3D torque between site's body and its parent body
    mjSENS_MAGNETOMETER                # 3D magnetometer
    mjSENS_RANGEFINDER                 # scalar distance to nearest geom or site along z-axis

                                     # sensors related to scalar joints tendons, actuators
    mjSENS_JOINTPOS                    # scalar joint position (hinge and slide only)
    mjSENS_JOINTVEL                    # scalar joint velocity (hinge and slide only)
    mjSENS_TENDONPOS                   # scalar tendon position
    mjSENS_TENDONVEL                   # scalar tendon velocity
    mjSENS_ACTUATORPOS                 # scalar actuator position
    mjSENS_ACTUATORVEL                 # scalar actuator velocity
    mjSENS_ACTUATORFRC                 # scalar actuator force

                                     # sensors related to ball joints
    mjSENS_BALLQUAT                    # 4D ball joint quaterion
    mjSENS_BALLANGVEL                  # 3D ball joint angular velocity

                                     # joint and tendon limit sensors in constraint space
    mjSENS_JOINTLIMITPOS               # joint limit distance-margin
    mjSENS_JOINTLIMITVEL               # joint limit velocity
    mjSENS_JOINTLIMITFRC               # joint limit force
    mjSENS_TENDONLIMITPOS              # tendon limit distance-margin
    mjSENS_TENDONLIMITVEL              # tendon limit velocity
    mjSENS_TENDONLIMITFRC              # tendon limit force

                                     # sensors attached to an object with spatial frame: (x)body geom, site, camera
    mjSENS_FRAMEPOS                    # 3D position
    mjSENS_FRAMEQUAT                   # 4D unit quaternion orientation
    mjSENS_FRAMEXAXIS                  # 3D unit vector: x-axis of object's frame
    mjSENS_FRAMEYAXIS                  # 3D unit vector: y-axis of object's frame
    mjSENS_FRAMEZAXIS                  # 3D unit vector: z-axis of object's frame
    mjSENS_FRAMELINVEL                 # 3D linear velocity
    mjSENS_FRAMEANGVEL                 # 3D angular velocity
    mjSENS_FRAMELINACC                 # 3D linear acceleration
    mjSENS_FRAMEANGACC                 # 3D angular acceleration

                                     # sensors related to kinematic subtrees; attached to a body (which is the subtree root)
    mjSENS_SUBTREECOM                  # 3D center of mass of subtree
    mjSENS_SUBTREELINVEL               # 3D linear velocity of subtree
    mjSENS_SUBTREEANGMOM               # 3D angular momentum of subtree

                                     # user-defined sensor
    mjSENS_USER                        # sensor data provided by mjcb_sensor callback
end

@cenum mjtStage::Cint begin           # computation stage
    mjSTAGE_NONE = 0            # no computations
    mjSTAGE_POS                        # position-dependent computations
    mjSTAGE_VEL                        # velocity-dependent computations
    mjSTAGE_ACC                        # acceleration/force-dependent computations
end

@cenum mjtDataType::Cint begin        # data type for sensors
    mjDATATYPE_REAL = 0            # real values, no constraints
    mjDATATYPE_POSITIVE                # positive values; 0 or negative: inactive
    mjDATATYPE_AXIS                    # 3D unit vector
    mjDATATYPE_QUATERNION              # unit quaternion
end

@cenum mjtLRMode::Cint begin          # mode for actuator length range computation
    mjLRMODE_NONE = 0                # do not process any actuators
    mjLRMODE_MUSCLE                    # process muscle actuators
    mjLRMODE_MUSCLEUSER                # process muscle and user actuators
    mjLRMODE_ALL                       # process all actuators
end

@uninitable mutable struct mjLROpt # options for mj_setLengthRange()
    # flags
    mode::Cint         # which actuators to process (mjtLRMode)
    useexisting::Cint  # use existing length range if available
    uselimit::Cint     # use joint and tendon limits if available

    # algorithm parameters
    accel::mjtNum      # target acceleration used to compute force
    maxforce::mjtNum   # maximum force 0: no limit
    timeconst::mjtNum  # time constant for velocity reduction min 0.01
    timestep::mjtNum   # simulation timestep 0: use mjOption.timestep
    inttotal::mjtNum   # total simulation time interval
    inteval::mjtNum    # evaluation time interval (at the end)
    tolrange::mjtNum   # convergence tolerance (relative to range)
end

@uninitable mutable struct mjVFS
    nfile::Cint
    filename::SVector{mjMAXVFS,SVector{mjMAXVFSNAME,UInt8}}
    filesize::SVector{mjMAXVFS,Cint}
    filedata::SVector{mjMAXVFS,Ptr{Cvoid}}
end

@uninitable struct mjOption
    timestep::mjtNum
    apirate::mjtNum

    impratio::mjtNum
    tolerance::mjtNum
    noslip_tolerance::mjtNum
    mpr_tolerance::mjtNum

    gravity::SVector{3,mjtNum}
    wind::SVector{3,mjtNum}
    magnetic::SVector{3,mjtNum}
    density::mjtNum
    viscosity::mjtNum

    o_margin::mjtNum
    o_solref::SVector{mjNREF,mjtNum}
    o_solimp::SVector{mjNIMP,mjtNum}

    integrator::Cint
    collision::Cint
    cone::Cint
    jacobian::Cint
    solver::Cint
    iterations::Cint
    noslip_iterations::Cint
    mpr_iterations::Cint
    disableflags::Cint
    enableflags::Cint
end

# Julia doesn't allow struct definition's inside of structs like in C,
# so we define the inner structs of mjVisual outside, prefixed with `mjVisual_`
@uninitable struct mjVisual_global
    fovy::Cfloat
    ipd::Cfloat
    linewidth::Cfloat
    glow::Cfloat
    offwidth::Cint
    offheight::Cint
end

@uninitable struct mjVisual_quality
    shadowsize::Cint
    offsamples::Cint
    numslices::Cint
    numstacks::Cint
    numquads::Cint
end

@uninitable struct mjVisual_headlight
    ambient::SVector{3,Cfloat}
    diffuse::SVector{3,Cfloat}
    specular::SVector{3,Cfloat}
    active::Cint
end

@uninitable struct mjVisual_map
    stiffness::Cfloat
    stiffnessrot::Cfloat
    force::Cfloat
    torque::Cfloat
    alpha::Cfloat
    fogstart::Cfloat
    fogend::Cfloat
    znear::Cfloat
    zfar::Cfloat
    haze::Cfloat
    shadowclip::Cfloat
    shadowscale::Cfloat
    actuatortendon::Cfloat
end

@uninitable struct mjVisual_scale
    forcewidth::Cfloat
    contactwidth::Cfloat
    contactheight::Cfloat
    connect::Cfloat
    com::Cfloat
    camera::Cfloat
    light::Cfloat
    selectpoint::Cfloat
    jointlength::Cfloat
    jointwidth::Cfloat
    actuatorlength::Cfloat
    actuatorwidth::Cfloat
    framelength::Cfloat
    framewidth::Cfloat
    constraint::Cfloat
    slidercrank::Cfloat
end

@uninitable struct mjVisual_rgba
    fog::SVector{4,Cfloat}
    haze::SVector{4,Cfloat}
    force::SVector{4,Cfloat}
    inertia::SVector{4,Cfloat}
    joint::SVector{4,Cfloat}
    actuator::SVector{4,Cfloat}
    actuatornegative::SVector{4,Cfloat}
    actuatorpositive::SVector{4,Cfloat}
    com::SVector{4,Cfloat}
    camera::SVector{4,Cfloat}
    light::SVector{4,Cfloat}
    selectpoint::SVector{4,Cfloat}
    connect::SVector{4,Cfloat}
    contactpoint::SVector{4,Cfloat}
    contactforce::SVector{4,Cfloat}
    contactfriction::SVector{4,Cfloat}
    contacttorque::SVector{4,Cfloat}
    contactgap::SVector{4,Cfloat}
    rangefinder::SVector{4,Cfloat}
    constraint::SVector{4,Cfloat}
    slidercrank::SVector{4,Cfloat}
    crankbroken::SVector{4,Cfloat}
end

@uninitable struct mjVisual
    _global::mjVisual_global
    quality::mjVisual_quality
    headlight::mjVisual_headlight
    map::mjVisual_map
    scale::mjVisual_scale
    rgba::mjVisual_rgba
end

@uninitable struct mjStatistic
    meaninertia::mjtNum
    meanmass::mjtNum
    meansize::mjtNum
    extent::mjtNum
    center::SVector{3,mjtNum}
end

# TODO immutable
mutable struct mjModel
    nq::Cint
    nv::Cint
    nu::Cint
    na::Cint
    nbody::Cint
    njnt::Cint
    ngeom::Cint
    nsite::Cint
    ncam::Cint
    nlight::Cint
    nmesh::Cint
    nmeshvert::Cint
    nmeshtexvert::Cint
    nmeshface::Cint
    nmeshgraph::Cint
    nskin::Cint
    nskinvert::Cint
    nskintexvert::Cint
    nskinface::Cint
    nskinbone::Cint
    nskinbonevert::Cint
    nhfield::Cint
    nhfielddata::Cint
    ntex::Cint
    ntexdata::Cint
    nmat::Cint
    npair::Cint
    nexclude::Cint
    neq::Cint
    ntendon::Cint
    nwrap::Cint
    nsensor::Cint
    nnumeric::Cint
    nnumericdata::Cint
    ntext::Cint
    ntextdata::Cint
    ntuple::Cint
    ntupledata::Cint
    nkey::Cint
    nuser_body::Cint
    nuser_jnt::Cint
    nuser_geom::Cint
    nuser_site::Cint
    nuser_cam::Cint
    nuser_tendon::Cint
    nuser_actuator::Cint
    nuser_sensor::Cint
    nnames::Cint

    nM::Cint
    nemax::Cint
    njmax::Cint
    nconmax::Cint
    nstack::Cint
    nuserdata::Cint
    nmocap::Cint
    nsensordata::Cint

    nbuffer::Cint

    opt::mjOption
    vis::mjVisual
    stat::mjStatistic

    buffer::Ptr{Cvoid}

    qpos0::Ptr{mjtNum}
    qpos_spring::Ptr{mjtNum}

    body_parentid::Ptr{Cint}
    body_rootid::Ptr{Cint}
    body_weldid::Ptr{Cint}
    body_mocapid::Ptr{Cint}
    body_jntnum::Ptr{Cint}
    body_jntadr::Ptr{Cint}
    body_dofnum::Ptr{Cint}
    body_dofadr::Ptr{Cint}
    body_geomnum::Ptr{Cint}
    body_geomadr::Ptr{Cint}
    body_simple::Ptr{mjtByte}
    body_sameframe::Ptr{mjtByte}
    body_pos::Ptr{mjtNum}
    body_quat::Ptr{mjtNum}
    body_ipos::Ptr{mjtNum}
    body_iquat::Ptr{mjtNum}
    body_mass::Ptr{mjtNum}
    body_subtreemass::Ptr{mjtNum}
    body_inertia::Ptr{mjtNum}
    body_invweight0::Ptr{mjtNum}
    body_user::Ptr{mjtNum}

    jnt_type::Ptr{Cint}
    jnt_qposadr::Ptr{Cint}
    jnt_dofadr::Ptr{Cint}
    jnt_bodyid::Ptr{Cint}
    jnt_group::Ptr{Cint}
    jnt_limited::Ptr{mjtByte}
    jnt_solref::Ptr{mjtNum}
    jnt_solimp::Ptr{mjtNum}
    jnt_pos::Ptr{mjtNum}
    jnt_axis::Ptr{mjtNum}
    jnt_stiffness::Ptr{mjtNum}
    jnt_range::Ptr{mjtNum}
    jnt_margin::Ptr{mjtNum}
    jnt_user::Ptr{mjtNum}

    dof_bodyid::Ptr{Cint}
    dof_jntid::Ptr{Cint}
    dof_parentid::Ptr{Cint}
    dof_Madr::Ptr{Cint}
    dof_simplenum::Ptr{Cint}
    dof_solref::Ptr{mjtNum}
    dof_solimp::Ptr{mjtNum}
    dof_frictionloss::Ptr{mjtNum}
    dof_armature::Ptr{mjtNum}
    dof_damping::Ptr{mjtNum}
    dof_invweight0::Ptr{mjtNum}
    dof_M0::Ptr{mjtNum}

    geom_type::Ptr{Cint}
    geom_contype::Ptr{Cint}
    geom_conaffinity::Ptr{Cint}
    geom_condim::Ptr{Cint}
    geom_bodyid::Ptr{Cint}
    geom_dataid::Ptr{Cint}
    geom_matid::Ptr{Cint}
    geom_group::Ptr{Cint}
    geom_priority::Ptr{Cint}
    geom_sameframe::Ptr{mjtByte}
    geom_solmix::Ptr{mjtNum}
    geom_solref::Ptr{mjtNum}
    geom_solimp::Ptr{mjtNum}
    geom_size::Ptr{mjtNum}
    geom_rbound::Ptr{mjtNum}
    geom_pos::Ptr{mjtNum}
    geom_quat::Ptr{mjtNum}
    geom_friction::Ptr{mjtNum}
    geom_margin::Ptr{mjtNum}
    geom_gap::Ptr{mjtNum}
    geom_user::Ptr{mjtNum}
    geom_rgba::Ptr{Cfloat}

    site_type::Ptr{Cint}
    site_bodyid::Ptr{Cint}
    site_matid::Ptr{Cint}
    site_group::Ptr{Cint}
    site_sameframe::Ptr{mjtByte}
    site_size::Ptr{mjtNum}
    site_pos::Ptr{mjtNum}
    site_quat::Ptr{mjtNum}
    site_user::Ptr{mjtNum}
    site_rgba::Ptr{Cfloat}

    cam_mode::Ptr{Cint}
    cam_bodyid::Ptr{Cint}
    cam_targetbodyid::Ptr{Cint}
    cam_pos::Ptr{mjtNum}
    cam_quat::Ptr{mjtNum}
    cam_poscom0::Ptr{mjtNum}
    cam_pos0::Ptr{mjtNum}
    cam_mat0::Ptr{mjtNum}
    cam_fovy::Ptr{mjtNum}
    cam_ipd::Ptr{mjtNum}
    cam_user::Ptr{mjtNum}

    light_mode::Ptr{Cint}
    light_bodyid::Ptr{Cint}
    light_targetbodyid::Ptr{Cint}
    light_directional::Ptr{mjtByte}
    light_castshadow::Ptr{mjtByte}
    light_active::Ptr{mjtByte}
    light_pos::Ptr{mjtNum}
    light_dir::Ptr{mjtNum}
    light_poscom0::Ptr{mjtNum}
    light_pos0::Ptr{mjtNum}
    light_dir0::Ptr{mjtNum}
    light_attenuation::Ptr{Cfloat}
    light_cutoff::Ptr{Cfloat}
    light_exponent::Ptr{Cfloat}
    light_ambient::Ptr{Cfloat}
    light_diffuse::Ptr{Cfloat}
    light_specular::Ptr{Cfloat}

    mesh_vertadr::Ptr{Cint}
    mesh_vertnum::Ptr{Cint}
    mesh_texcoordadr::Ptr{Cint}
    mesh_faceadr::Ptr{Cint}
    mesh_facenum::Ptr{Cint}
    mesh_graphadr::Ptr{Cint}
    mesh_vert::Ptr{Cfloat}
    mesh_normal::Ptr{Cfloat}
    mesh_texcoord::Ptr{Cfloat}
    mesh_face::Ptr{Cint}
    mesh_graph::Ptr{Cint}

    skin_matid::Ptr{Cint}
    skin_rgba::Ptr{Cfloat}
    skin_inflate::Ptr{Cfloat}
    skin_vertadr::Ptr{Cint}
    skin_vertnum::Ptr{Cint}
    skin_texcoordadr::Ptr{Cint}
    skin_faceadr::Ptr{Cint}
    skin_facenum::Ptr{Cint}
    skin_boneadr::Ptr{Cint}
    skin_bonenum::Ptr{Cint}
    skin_vert::Ptr{Cfloat}
    skin_texcoord::Ptr{Cfloat}
    skin_face::Ptr{Cint}
    skin_bonevertadr::Ptr{Cint}
    skin_bonevertnum::Ptr{Cint}
    skin_bonebindpos::Ptr{Cfloat}
    skin_bonebindquat::Ptr{Cfloat}
    skin_bonebodyid::Ptr{Cint}
    skin_bonevertid::Ptr{Cint}
    skin_bonevertweight::Ptr{Cfloat}

    hfield_size::Ptr{mjtNum}
    hfield_nrow::Ptr{Cint}
    hfield_ncol::Ptr{Cint}
    hfield_adr::Ptr{Cint}
    hfield_data::Ptr{Cfloat}

    tex_type::Ptr{Cint}
    tex_height::Ptr{Cint}
    tex_width::Ptr{Cint}
    tex_adr::Ptr{Cint}
    tex_rgb::Ptr{mjtByte}

    mat_texid::Ptr{Cint}
    mat_texuniform::Ptr{mjtByte}
    mat_texrepeat::Ptr{Cfloat}
    mat_emission::Ptr{Cfloat}
    mat_specular::Ptr{Cfloat}
    mat_shininess::Ptr{Cfloat}
    mat_reflectance::Ptr{Cfloat}
    mat_rgba::Ptr{Cfloat}

    pair_dim::Ptr{Cint}
    pair_geom1::Ptr{Cint}
    pair_geom2::Ptr{Cint}
    pair_signature::Ptr{Cint}
    pair_solref::Ptr{mjtNum}
    pair_solimp::Ptr{mjtNum}
    pair_margin::Ptr{mjtNum}
    pair_gap::Ptr{mjtNum}
    pair_friction::Ptr{mjtNum}

    exclude_signature::Ptr{Cint}

    eq_type::Ptr{Cint}
    eq_obj1id::Ptr{Cint}
    eq_obj2id::Ptr{Cint}
    eq_active::Ptr{mjtByte}
    eq_solref::Ptr{mjtNum}
    eq_solimp::Ptr{mjtNum}
    eq_data::Ptr{mjtNum}

    tendon_adr::Ptr{Cint}
    tendon_num::Ptr{Cint}
    tendon_matid::Ptr{Cint}
    tendon_group::Ptr{Cint}
    tendon_limited::Ptr{mjtByte}
    tendon_width::Ptr{mjtNum}
    tendon_solref_lim::Ptr{mjtNum}
    tendon_solimp_lim::Ptr{mjtNum}
    tendon_solref_fri::Ptr{mjtNum}
    tendon_solimp_fri::Ptr{mjtNum}
    tendon_range::Ptr{mjtNum}
    tendon_margin::Ptr{mjtNum}
    tendon_stiffness::Ptr{mjtNum}
    tendon_damping::Ptr{mjtNum}
    tendon_frictionloss::Ptr{mjtNum}
    tendon_lengthspring::Ptr{mjtNum}
    tendon_length0::Ptr{mjtNum}
    tendon_invweight0::Ptr{mjtNum}
    tendon_user::Ptr{mjtNum}
    tendon_rgba::Ptr{Cfloat}

    wrap_type::Ptr{Cint}
    wrap_objid::Ptr{Cint}
    wrap_prm::Ptr{mjtNum}

    actuator_trntype::Ptr{Cint}
    actuator_dyntype::Ptr{Cint}
    actuator_gaintype::Ptr{Cint}
    actuator_biastype::Ptr{Cint}
    actuator_trnid::Ptr{Cint}
    actuator_group::Ptr{Cint}
    actuator_ctrllimited::Ptr{mjtByte}
    actuator_forcelimited::Ptr{mjtByte}
    actuator_dynprm::Ptr{mjtNum}
    actuator_gainprm::Ptr{mjtNum}
    actuator_biasprm::Ptr{mjtNum}
    actuator_ctrlrange::Ptr{mjtNum}
    actuator_forcerange::Ptr{mjtNum}
    actuator_gear::Ptr{mjtNum}
    actuator_cranklength::Ptr{mjtNum}
    actuator_acc0::Ptr{mjtNum}
    actuator_length0::Ptr{mjtNum}
    actuator_lengthrange::Ptr{mjtNum}
    actuator_user::Ptr{mjtNum}

    sensor_type::Ptr{Cint}
    sensor_datatype::Ptr{Cint}
    sensor_needstage::Ptr{Cint}
    sensor_objtype::Ptr{Cint}
    sensor_objid::Ptr{Cint}
    sensor_dim::Ptr{Cint}
    sensor_adr::Ptr{Cint}
    sensor_cutoff::Ptr{mjtNum}
    sensor_noise::Ptr{mjtNum}
    sensor_user::Ptr{mjtNum}

    numeric_adr::Ptr{Cint}
    numeric_size::Ptr{Cint}
    numeric_data::Ptr{mjtNum}

    text_adr::Ptr{Cint}
    text_size::Ptr{Cint}
    text_data::Ptr{UInt8}

    tuple_adr::Ptr{Cint}
    tuple_size::Ptr{Cint}
    tuple_objtype::Ptr{Cint}
    tuple_objid::Ptr{Cint}
    tuple_objprm::Ptr{mjtNum}

    key_time::Ptr{mjtNum}
    key_qpos::Ptr{mjtNum}
    key_qvel::Ptr{mjtNum}
    key_act::Ptr{mjtNum}

    name_bodyadr::Ptr{Cint}
    name_jntadr::Ptr{Cint}
    name_geomadr::Ptr{Cint}
    name_siteadr::Ptr{Cint}
    name_camadr::Ptr{Cint}
    name_lightadr::Ptr{Cint}
    name_meshadr::Ptr{Cint}
    name_skinadr::Ptr{Cint}
    name_hfieldadr::Ptr{Cint}
    name_texadr::Ptr{Cint}
    name_matadr::Ptr{Cint}
    name_pairadr::Ptr{Cint}
    name_excludeadr::Ptr{Cint}
    name_eqadr::Ptr{Cint}
    name_tendonadr::Ptr{Cint}
    name_actuatoradr::Ptr{Cint}
    name_sensoradr::Ptr{Cint}
    name_numericadr::Ptr{Cint}
    name_textadr::Ptr{Cint}
    name_tupleadr::Ptr{Cint}
    name_keyadr::Ptr{Cint}
    names::Ptr{UInt8}
end


text:
filename() = "day15.txt"

"""
Read data and return a pair of data structures:
1. A dictionary of tuples where K=spoken number, V=(lastpos,recentpos)
2. An array of spoken numbers
"""
function read_data()
    seq = parse.(Int, split(readline(filename()), ","))
    return Dict(v => (0, k) for (k,v) in enumerate(seq)), seq
end

"""
Solve the puzzle with a simple loop. Since we always spoken the last
spoken number, the Dict is guaranteed to have a value.

If the first value of the tuple is 0, then it is a starting value and say 0;
otherwise, say the difference of values in that tuple. Then, update the
Dict for the spoken number with a new tuple.
"""
function solve(year)
    dct, seq = read_data()
    start = length(seq) + 1
    spoken = seq[end]
    for i in start:year
        tup = dct[spoken]
        val = tup[1] == 0 ? 0 : tup[2] - tup[1]
        if haskey(dct, val)
            dct[val] = (dct[val][2], i)
        else
            dct[val] = (0, i)
        end
        spoken = val
    end
    return spoken
end

part1() = solve(2020)
part2() = solve(30000000)

#=
julia> @benchmark part2()
BenchmarkTools.Trial:
  memory estimate:  395.84 MiB
  allocs estimate:  84
  --------------
  minimum time:     2.576 s (0.32% GC)
  median time:      2.583 s (1.69% GC)
  mean time:        2.583 s (1.69% GC)
  maximum time:     2.591 s (3.05% GC)
  --------------
  samples:          2
  evals/sample:     1
=#

using Plots

function animation(year)
    anim = Animation()
    dct, seq = read_data()
    start = length(seq) + 1
    spoken = seq[end]
    plot(seq, background_color = :black, legend = :none,
        xlabel = "Year", ylabel = "Number Spoken",
        title = "Advent of Code (Day 15)",
        size = (500,300),
        xlims = (1, 2100),
        ylims = (1, 1600))
    for i in start:year
        tup = dct[spoken]
        val = tup[1] == 0 ? 0 : tup[2] - tup[1]
        if haskey(dct, val)
            dct[val] = (dct[val][2], i)
        else
            dct[val] = (0, i)
        end
        spoken = val
        push!(seq, spoken)
        if i % 20 == 0
            plot!(seq[1:i])
            frame(anim)
        end
    end
    return gif(anim, "day15_anim.gif"; fps = 25)
end

ext: jl


text:
struct Ship
    x::Int64
    y::Int64
    waypoint::NamedTuple{(:rx, :ry), Tuple{Int64, Int64}}
end

function parse_instructions(lines)
    instructions = NamedTuple{(:action, :value), Tuple{Char, Int64}}[]
    for line in lines
        a = line[1]
        v = parse(Int64, line[2:end])
        push!(instructions, (action=a, value=v))
    end
    instructions
end

function rotate_waypoint(ship, value)
    deg = value < 0 ? value + 360 : value
    if deg == 0
        ship
    elseif deg == 90
        Ship(ship.x, ship.y, (rx=-ship.waypoint.ry, ry=ship.waypoint.rx))
    elseif deg == 180
        Ship(ship.x, ship.y, (rx=-ship.waypoint.rx, ry=-ship.waypoint.ry))
    elseif deg == 270
        Ship(ship.x, ship.y, (rx=ship.waypoint.ry, ry=-ship.waypoint.rx))
    end
end

function move(ship, instruction)
    if instruction.action == 'N'
        Ship(ship.x, ship.y, (rx=ship.waypoint.rx, ry=ship.waypoint.ry + instruction.value))
    elseif instruction.action == 'S'
        Ship(ship.x, ship.y, (rx=ship.waypoint.rx, ry=ship.waypoint.ry - instruction.value))
    elseif instruction.action == 'E'
        Ship(ship.x, ship.y, (rx=ship.waypoint.rx + instruction.value, ry=ship.waypoint.ry))
    elseif instruction.action == 'W'
        Ship(ship.x, ship.y, (rx=ship.waypoint.rx - instruction.value, ry=ship.waypoint.ry))
    elseif instruction.action == 'L'
        rotate_waypoint(ship, instruction.value)
    elseif instruction.action == 'R'
        rotate_waypoint(ship, -instruction.value)
    elseif instruction.action == 'F'
        x = ship.x + instruction.value * ship.waypoint.rx
        y = ship.y + instruction.value * ship.waypoint.ry
        Ship(x, y, (rx=ship.waypoint.rx, ry=ship.waypoint.ry))
    end
end

function navigate(instructions)
    initial = Ship(0, 0, (rx=10, ry=1))
    ship = foldl(move, instructions; init=initial)
    abs(ship.x - initial.x) + abs(ship.y - initial.y)
end

function main()
    lines = readlines("../resources/input_12.txt")
    instructions = parse_instructions(lines)
    result = navigate(instructions)
    println(result)
end

main()


text:
mutable struct Explicit <: TuningStrategy end

struct ExplicitState{R, N}
    range::R # a model-generating iterator
    next::N # to hold output of `iterate(range)`
end

ExplictState(r::R, n::N) where {R,N} = ExplicitState{R, Union{Nothing, N}}(r, n)

function MLJTuning.setup(tuning::Explicit, model, range, n, verbosity)
    next = iterate(range)
    return ExplicitState(range, next)
end

# models! returns all available models in the range at once:
function MLJTuning.models(tuning::Explicit,
                           model,
                           history,
                           state,
                           n_remaining,
                           verbosity)

    range, next  = state.range, state.next

    next === nothing && return nothing, state

    m, s = next
    models = [m, ]

    next = iterate(range, s)

    i = 1 # current length of `models`
    while i < n_remaining
        next === nothing && break
        m, s = next
        push!(models, m)
        i += 1
        next = iterate(range, s)
    end

    new_state = ExplicitState(range, next)

    return models, new_state

end

function default_n(tuning::Explicit, range)
    try
        length(range)
    catch MethodError
        DEFAULT_N
    end
end


text:
# functions from GasChromatographyTools.jl:
"""
	common(s_1, s_1)

Compare two arrays and return the common elements.
"""
function common(s_1, s_2)
	common_s = String[]
	if length(s_1) >= length(s_2)
		for i=1:length(s_1)
			if s_1[i] in s_2
				push!(common_s, s_1[i])
			end
		end
	else
		for i=1:length(s_2)
			if s_2[i] in s_1
				push!(common_s, s_2[i])
			end
		end
	end
	return common_s
end

"""
	compare_peaklist(pl_1, pl_2)

Compare two peaklists (results of GasChromatographySimulator.jl) and calculate
absolute and relative differences of retention times and peak widths.
"""
function compare_peaklist(pl_1, pl_2)
	name = pl_1.Name
	tR1 = pl_1.tR
	τR1 = pl_1.τR
	tR2 = Array{Float64}(undef, size(pl_1)[1])
	τR2 = Array{Float64}(undef, size(pl_1)[1])
	for i=1:size(pl_1)[1]
		i2 = findfirst(name[i].==pl_2.Name)
		tR2[i] = pl_2.tR[i2]
		τR2[i] = pl_2.τR[i2]
	end
	ΔtR = tR1 .- tR2
	ΔτR = τR1 .- τR2
	rel_tR = ΔtR.*100.0./tR1
	rel_τR = ΔτR.*100.0./τR1
	compare_pl = DataFrame(Name=name, tR1=tR1, tR2=tR2, ΔtR=ΔtR, rel_tR=rel_tR, τR1=τR1, τR2=τR2, ΔτR=ΔτR, rel_τR=rel_τR)
	return compare_pl
end

"""
	compare_measurement_simulation(meas, peaklist)

Compare the retention times of measured and simulated substances.

# Arguments
* `meas`: DataFrame consisting at least of the columns `:Name` and `:RT`
  (measured retention time in s)
* `peaklist`: DataFrame as result from GasChromatographySimulator.jl with the
  columns `:Name` and `:tR` (simulated retention time in s).

The comparison is done by searching the same `Name` of the substance in both
DataFrames and calculating the absolute difference of the retention times (in s)
and the relative difference (in %).
"""
function compare_measurement_simulation(meas, peaklist)
	name = meas.Name
	tRm = meas.RT
	tRs = Array{Float64}(undef, size(meas)[1])
	for i=1:size(meas)[1]
		i2 = findfirst(name[i].==peaklist.Name)
		if typeof(i2) == Nothing
			tRs[i] = NaN
		else
			tRs[i] = peaklist.tR[i2]
		end
	end
	compare_df = DataFrame(Name=name, measured_tR=tRm, simulated_tR=tRs, ΔtR=tRm.-tRs, rel_tR=(tRm.-tRs)./tRm.*100.0)
	return compare_df
end


text:
function run_naive(num_sims::Integer) :: Float64

    # `x` and `y` are vectors
    x = rand(Float64, num_sims)
    y = rand(Float64, num_sims)

    # dot-syntax: broadcast
    r = sqrt.( x.^2 + y.^2 )

    area = mapreduce(
                     el -> el <= 1.0 ? 1.0 : 0.0, # map
                     +, # reduce
                     r) / num_sims
    return 4*area
end


text:
function mountaincar_run_online_model_search(rock_c::Float64,
                                             num_eval_samples::Int64;
                                             max_likelihood = false,
                                             optimistic = false,
                                             model_evaluation = false,
                                             debug = false,
                                             seed = 0)
    rng = MersenneTwister(seed)
    mountaincar = MountainCar(rock_c)
    model = MountainCar(0.0)
    agent = MountainCarOnlineModelSearchAgent(mountaincar, model,
                                              num_eval_samples)
    n_steps = run(agent, rng; debug = debug, max_likelihood = max_likelihood,
                  optimistic = optimistic, model_evaluation = model_evaluation,)
    n_steps
end

function mountaincar_run_cmax(rock_c::Float64;
                              debug = false,
                              seed = 0)
    rng = MersenneTwister(seed)
    mountaincar = MountainCar(rock_c)
    model = MountainCar(0.0)
    agent = MountainCarOnlineModelSearchAgent(mountaincar, model,
                                              10)
    n_steps = run_cmax(agent, rng; debug = debug)
    n_steps
end

function mountaincar_run_true(rock_c::Float64;
                              debug = false,
                              seed = 0)
    rng = MersenneTwister(seed)
    mountaincar = MountainCar(rock_c)
    model = MountainCar(0.0)
    agent = MountainCarOnlineModelSearchAgent(mountaincar, model, 10)
    n_steps = run_true(agent, rng; debug = debug)
    n_steps
end

function mountaincar_run_online_experiments()
    num_eval_samples = 50
    seeds = collect(1:10)
    # ml_steps = []
    # oms_steps = []
    # cmax_steps = []
    # true_steps = []
    # rbms_steps = []
    me_steps = []
    rock_c_values = collect(0.02:0.002:0.03)
    for rock_c in rock_c_values
        # ml_sub_steps = []
        # oms_sub_steps = []
        # cmax_sub_steps = []
        # true_sub_steps = []
        # rbms_sub_steps = []
        me_sub_steps = []
        for seed in seeds
            # println("MLE with rock_c ", rock_c, " seed ", seed)
            # push!(ml_sub_steps, mountaincar_run_online_model_search(
            #     rock_c, num_eval_samples; max_likelihood = true, seed = seed,
            #     debug = false,
            # ))
            # println("OMS with rock_c ", rock_c, " seed ", seed)
            # push!(oms_sub_steps, mountaincar_run_online_model_search(
            #     rock_c, num_eval_samples; optimistic = true, seed = seed,
            #     debug = false,
            # ))
            # println("CMAX with rock_c ", rock_c, " seed ", seed)
            # push!(cmax_sub_steps, mountaincar_run_cmax(
            #     rock_c; seed = seed, debug = false,
            # ))
            # println("True with rock_c ", rock_c, " seed ", seed)
            # push!(true_sub_steps, mountaincar_run_true(
            #     rock_c; seed = seed, debug = false,
            # ))
            # println("RBMS with rock_c ", rock_c, " seed ", seed)
            # push!(rbms_sub_steps, mountaincar_run_online_model_search(
            #     rock_c, num_eval_samples; optimistic = false, seed = seed,
            #     debug = false,
            # ))
            println("ME with rock_c ", rock_c, " seed ", seed)
            push!(me_sub_steps, mountaincar_run_online_model_search(
                rock_c, num_eval_samples; model_evaluation = true, seed = seed,
                debug = false,
            ))
        end
        # push!(ml_steps, ml_sub_steps)
        # push!(oms_steps, oms_sub_steps)
        # push!(cmax_steps, cmax_sub_steps)
        # push!(true_steps, true_sub_steps)
        # push!(rbms_steps, rbms_sub_steps)
        push!(me_steps, me_sub_steps)
    end
    # ml_mean_steps = [mean(x) for x in ml_steps]
    # ml_std_steps = [std(x) for x in ml_steps] ./ sqrt(length(seeds))
    # oms_mean_steps = [mean(x) for x in oms_steps]
    # oms_std_steps = [std(x) for x in oms_steps] ./ sqrt(length(seeds))
    # cmax_mean_steps = [mean(x) for x in cmax_steps]
    # cmax_std_steps = [std(x) for x in cmax_steps] ./ sqrt(length(seeds))
    # true_mean_steps = [mean(x) for x in true_steps]
    # true_std_steps = [std(x) for x in true_steps] ./ sqrt(length(seeds))
    # rbms_mean_steps = [mean(x) for x in rbms_steps]
    # rbms_std_steps = [std(x) for x in rbms_steps] ./ sqrt(length(seeds))
    me_mean_steps = [mean(x) for x in me_steps]
    me_std_steps = [std(x) for x in me_steps] ./ sqrt(length(seeds))

    # println("MLE")
    # println(ml_steps)
    # println("OMS")
    # println(oms_steps)
    # println("CMAX")
    # println(cmax_steps)
    # println("TRUE")
    # println(true_steps)
    # println("RBMS")
    # println(rbms_steps)
    println("ME")
    println(me_steps)

    # println("MLE")
    # println(ml_mean_steps)
    # println("OMS")
    # println(oms_mean_steps)
    # println("CMAX")
    # println(cmax_mean_steps)
    # println("TRUE")
    # println(true_mean_steps)
    # println("RBMS")
    # println(rbms_mean_steps)
    println("ME")
    println(me_mean_steps)

    # println("MLE")
    # println(ml_std_steps)
    # println("OMS")
    # println(oms_std_steps)
    # println("CMAX")
    # println(cmax_std_steps)
    # println("TRUE")
    # println(true_std_steps)
    # println("RBMS")
    # println(rbms_std_steps)
    println("ME")
    println(me_std_steps)

    # plot(rock_c_values, ml_mean_steps, grid=true, yerror=ml_std_steps, lw=3,
    #      label="MLE")
    # plot!(rock_c_values, oms_mean_steps, grid=true, yerror=oms_std_steps, lw=3,
    #      label="TOMS")
    # plot!(rock_c_values, cmax_mean_steps, grid=true, yerror=cmax_std_steps, lw=3,
    #       label="CMAX")
    # plot!(rock_c_values, true_mean_steps, grid=true, yerror=true_std_steps, lw=3,
    #      label="TRUE")
    # plot!(rock_c_values, rbms_mean_steps, grid=true, yerror=rbms_std_steps, lw=3,
    #      label="RBMS")
    # xlabel!("Misspecification")
    # ylabel!("Number of steps to reach goal")
end


text:
"""
	LatticeModule
Store information about the Lattice

# Type list

* LatticeInfo
* Latttice
* ClassicalLattice,
* Classical2dLattice,
* Classical2dSquareLattice,
* Classical2dFractalLattice,
* QuantumLattice,
* Quantum2dLattice,
* Quantum2dSquareLattice,
* Classical3dLattice,
* Classical3dSquareLattice,
* Quantum2dFractalLattice,
* FractalLattice,

# Method list

* getDimension
* getGeometry
* getCoarserate
* getLegextension
* getQuantumOrClassical
* isClassical,
* isTwoDimension,
* isSquareLattice,
* isFractalLattice,
* buildLattice,
* initialize!,
* isZeroTemperature,
* isClassical,
* getModelname,
* getStates,
* isSymmetricFactorization,
* isPottsModel,
* getHamiltonian,
* getTemperature,
* getExternalfield,
* getEnvParameters,
* setTemperature!,
* setExternalfield!,
* setEnvParameters!,
* initialize!,
* getTensorT,
* getFactorW,
* getFactorWp,
* setTensorT!,
* makeTensorT,
* makeFactorWp,
* makeTensorQ,
* buildClassicalLattice,
* getSpaceDimension,
* buildQuantumLattice,
* getTrotterparameter,
* getTrotteriteration,
* getTrotterlayers,
* getBeta,
* setTrotterparameter!,
* setTrotterIteration!,
* getTensorP,
* getTensorQ,
* getLegextension,
* getHausdorffDim,
* getFractalDim,
* setTensorP!,
* setTensorQ!

# Test list

* testlatticeinfo,
* testClassical2dSquareLattice,
* testClassical2dFractalLattice,
* testQuantum2dSquareLattice,
* testClassical3dSquareLattice,
* testQuantum2dFractalLattice,

"""
module LatticeModule
	using TensorOperations
	using TensorMatrices_lemon
	using HoTRG_lemon.SpinModule

	export
		# lattice_info.jl
		LatticeInfo,
		getDimension,
		getGeometry,
		getCoarserate,
		getLegextension,
		getQuantumOrClassical,
		isClassical,
		isTwoDimension,
		isSquareLattice,
		isFractalLattice,
		testlatticeinfo,
		# lattice.jl
		Lattice,
		buildLattice,
			# <inherit from spinmodel>
		initialize!,
		isZeroTemperature,
		isClassical,
		getModelname,
		getStates,
		isSymmetricFactorization,
		isPottsModel,
		getHamiltonian,
		getTemperature,
		getExternalfield,
		getEnvParameters,
		setTemperature!,
		setExternalfield!,
		setEnvParameters!,
		initialize!,
		getTensorT,
		getFactorW,
		getFactorWp,
		setTensorT!,
		# lattice_classical_2d_localtensor.jl
		makeTensorT,
		makeFactorWp,
		makeTensorQ,
		# lattice_classical.jl
		ClassicalLattice,
		buildClassicalLattice,
		# lattice_classical_2d.jl
		Classical2dLattice,
		getSpaceDimension,
		# lattice_classical_2d_square.jl
		Classical2dSquareLattice,
		testClassical2dSquareLattice,
		# lattice_classical_2d_fractal.jl
		Classical2dFractalLattice,
		testClassical2dFractalLattice,
		# lattice_quantum.jl
		QuantumLattice,
		buildQuantumLattice,
			# <inherit from QuantumSpinModel>
		getTrotterparameter,
		getTrotteriteration,
		getTrotterlayers,
		getBeta,
		setTrotterparameter!,
		setTrotterIteration!,
		# lattice_quantum_2d.jl
		Quantum2dLattice,
		# lattice_quantum_2d_square.jl
		Quantum2dSquareLattice,
		testQuantum2dSquareLattice,
		# lattice_classical_3d.jl
		Classical3dLattice,
		# lattice_classical_3d_square.jl
		Classical3dSquareLattice,
		testClassical3dSquareLattice,
		# lattice_quantum_2d_fractal.jl
		Quantum2dFractalLattice,
		testQuantum2dFractalLattice,
		# lattice_fractal.jl
		FractalLattice,
		getTensorP,
		getTensorQ,
		getLegextension,
		getHausdorffDim,
		getFractalDim,
		setTensorP!,
		setTensorQ!

	include("lattice_info.jl")
	include("lattice.jl")
	include("lattice_classical_2d_localtensor.jl")
	include("lattice_classical.jl")
	include("lattice_classical_2d.jl")
	include("lattice_classical_2d_square.jl")
	include("lattice_classical_2d_fractal.jl")
	include("lattice_quantum.jl")
	include("lattice_quantum_2d.jl")
	include("lattice_quantum_2d_square.jl")
	include("lattice_classical_3d.jl")
	include("lattice_classical_3d_square.jl")
	include("lattice_quantum_2d_fractal.jl")
	include("lattice_fractal.jl")
end

# inherit getSpaceDimension to the Simulator
#		getTensorP,
#		getTensorQ,
#		getLegextension,
#		getHausdorffDim,
#		getFractalDim,
#		setTensorP!,
#		setTensorQ!


text:
export ft_setbitmode
export FT_BITMODE_RESET, FT_BITMODE_ASYNC_BITBANG
export FT_BITMODE_MPSSE, FT_BITMODE_SYNC_BITBANG
export FT_BITMODE_MCU_HOST, FT_BITMODE_FAST_SERIAL
export FT_BITMODE_CBUS_BITBANG, FT_BITMODE_SYNC_FIFO

#Bit Modes (see ft_setbitmode)
const FT_BITMODE_RESET = 0x00
const FT_BITMODE_ASYNC_BITBANG = 0x01
const FT_BITMODE_MPSSE = 0x02
const FT_BITMODE_SYNC_BITBANG = 0x04
const FT_BITMODE_MCU_HOST = 0x08
const FT_BITMODE_FAST_SERIAL = 0x10
const FT_BITMODE_CBUS_BITBANG = 0x20
const FT_BITMODE_SYNC_FIFO = 0x40

function ft_setbitmode(ft_handle::Culong, mask::UInt8, mode::UInt8)
  @assert mode in keys(ftbitmodedict)
  ft_status = ccall((:FT_SetBitMode, d2xx),
                     Cuint,
                     (Culong, Cuchar, Cuchar),
                     ft_handle, mask, mode)
  checkstatus(ft_status)
  return nothing
end


text:
export Dirac

struct Dirac{X} <: AbstractMeasure{X}
    supp :: X
end

Dirac(x::X) = Dirac{X}(x)

eltype(μ::Dirac{X}) = X

function (μ::Dirac{X})(s) where {X}
    μ.supp ∈ s && return 1
    return 0
end

function logdensity(μ::Dirac{X}, x::X)
    μ.supp ∈ s && return 0.0
    return -Inf
end


text:
#--------------------------------------------------------------------
# Spacetime Discretization methods in Julia
# Soham 08-2018
# Define the Abstract Types
#--------------------------------------------------------------------
abstract type Manifold{Tag} end
abstract type Space{Tag, N} <: Manifold{Tag} end

struct U end
struct V end
struct Singularity end

"""
    Struct to hold a field
    in space S, dimension D and eltype T.
"""
struct Field{S, D, T}
    space::S
    value::AbstractArray{T, D}
end

"""
    Struct to hold a field
    in space S, dimension D and eltype T.
"""
struct Operator{S, D, T}
    space::S
    value::AbstractArray{T, D}
end

struct ProductSpace{T1, T2}
    S1::T1
    S2::T2
end

mutable struct Parameters{T}
    npoints::NTuple{2, Int}
    npatchs::NTuple{2, Int}
    ubounds::NTuple{2, T}
    vbounds::NTuple{2, T}
    nfields::Int
end

struct ChebyshevGL{Tag, N, T} <: Space{Tag, N}
    min::T
    max::T
    ChebyshevGL{Tag, N, T}(min, max) where {Tag, N, T} = max > min ? new{Tag, N, T}(min, max) : error("bounds are out of order")
end

struct Chebyshev{Tag, N, T} <: Space{Tag, N}
    min::T
    max::T
    Chebyshev{Tag, N, T}(min, max) where {Tag, N, T} = max > min ? new{Tag, N, T}(min, max) : error("bounds are out of order")
end

Cardinal{Tag, N, T} = Union{ChebyshevGL{Tag, N, T},
                              Chebyshev{Tag, N, T}}


text:
@block KimFung ["attributes"] begin

    @cell "align" [scatter, text, align, "2d"] begin
        scene = scatter(rand(10), color=:red)
        text!(scene,"adding text",textsize = 0.6, align = (:center, :center))
    end

    @cell "fillrange" [contour, fillrange, "2d"] begin
        x = LinRange(-1, 1, 20)
        y = LinRange(-1, 1, 20)
        z = x .* y'
        contour(x, y, z, levels = 0, linewidth = 0, fillrange = true)
    end

    @cell "font" [text, scatter, font, "2d"] begin
        scene = Scene()
        scatter!(scene, rand(10), color=:red)
        text!(scene,"adding text",textsize = 0.6, align = (:center, :center), font = "Blackchancery")
    end

    @cell "glowcolor, glowwidth" [scatter, glowcolor, glowwidth, "2d"] begin
        scatter(randn(10),color=:blue, glowcolor = :orange, glowwidth = 10)
    end

    @cell "isorange, isovalue" [volume, algorithm, isorange, isovalue] begin
        r = range(-1, stop = 1, length = 100)
        matr = [(x.^2 + y.^2 + z.^2) for x = r, y = r, z = r]
        volume(matr .* (matr .> 1.4), algorithm = :iso, isorange = 0.05, isovalue = 1.7)
    end

    @cell "levels" [contour, colormap, levels, "2d"] begin
        x = LinRange(-1, 1, 20)
        y = LinRange(-1, 1, 20)
        z = x .* y'
        contour(x, y, z, linewidth = 3, colormap = :colorwheel, levels = 50)
    end

    @cell "position" [scatter, text, position, "2d"] begin
        scene = Scene()
        scatter!(scene, rand(10), color=:red)
        text!(scene,"adding text",textsize = 0.6, position = (5.0, 1.1))
    end

    @cell "rotation" [text, rotation, "2d"] begin
        text("Hello World", rotation = 1.1)
    end

    @cell "shading" [mesh, sphere, shading, "2d"] begin
        mesh(Sphere(Point3f0(0), 1f0), color = :orange, shading = false)
    end

    @cell "visible" [scatter, visible, "2d"] begin
        vbox(
            scatter(randn(20), color = to_colormap(:deep, 20), markersize = 1, visible = true),
            scatter(randn(20), color = to_colormap(:deep, 20), markersize = 1, visible = false)
        )
    end
end


text:
function ifthenelse_singlevariable_test(
    model::MOI.ModelLike,
    config::MOIT.Config,
)
    MOI.empty!(model)

    @test MOI.supports_constraint(model, MOI.VariableIndex, MOI.Integer)
    @test MOI.supports_constraint(model, MOI.VariableIndex, MOI.EqualTo{Int})
    @test MOI.supports_constraint(
        model,
        MOI.VectorOfVariables,
        CP.IfThenElse{MOI.EqualTo{Int}, MOI.EqualTo{Int}, MOI.EqualTo{Int}},
    )

    x1, _ = MOI.add_constrained_variable(model, MOI.Integer())
    x2, _ = MOI.add_constrained_variable(model, MOI.Integer())
    x3, _ = MOI.add_constrained_variable(model, MOI.Integer())

    c1 = MOI.add_constraint(model, x1, MOI.EqualTo(1))
    c2 = MOI.add_constraint(
        model,
        MOI.VectorOfVariables([x1, x2, x3]),
        CP.IfThenElse(MOI.EqualTo(1), MOI.EqualTo(1), MOI.EqualTo(0)),
    )

    @test MOI.is_valid(model, x1)
    @test MOI.is_valid(model, x2)
    @test MOI.is_valid(model, x3)
    @test MOI.is_valid(model, c1)
    @test MOI.is_valid(model, c2)

    if config.solve
        MOI.optimize!(model)
        @test MOI.get(model, MOI.TerminationStatus()) == MOI.OPTIMAL
        @test MOI.get(model, MOI.PrimalStatus()) == MOI.FEASIBLE_POINT

        @test MOI.get(model, MOI.ResultCount()) >= 1
        @test MOI.get(model, MOI.VariablePrimal(), x1) == 1
        @test MOI.get(model, MOI.VariablePrimal(), x2) == 1
    end
end

function ifthenelse_scalaraffinefunction_test(
    model::MOI.ModelLike,
    config::MOIT.Config,
)
    MOI.empty!(model)

    @test MOI.supports_constraint(model, MOI.VariableIndex, MOI.Integer)
    @test MOI.supports_constraint(
        model,
        MOI.ScalarAffineFunction{Int},
        MOI.EqualTo{Int},
    )
    @test MOI.supports_constraint(
        model,
        MOI.VectorAffineFunction{Int},
        CP.IfThenElse{MOI.EqualTo{Int}, MOI.EqualTo{Int}, MOI.EqualTo{Int}},
    )

    x1, _ = MOI.add_constrained_variable(model, MOI.Integer())
    x2, _ = MOI.add_constrained_variable(model, MOI.Integer())
    x3, _ = MOI.add_constrained_variable(model, MOI.Integer())

    c1 = MOI.add_constraint(model, 1 * x1, MOI.EqualTo(1))
    c2 = MOI.add_constraint(
        model,
        MOIU.vectorize([x1, x2, x3]),
        CP.IfThenElse(MOI.EqualTo(1), MOI.EqualTo(1), MOI.EqualTo(0)),
    )

    @test MOI.is_valid(model, x1)
    @test MOI.is_valid(model, x2)
    @test MOI.is_valid(model, x3)
    @test MOI.is_valid(model, c1)
    @test MOI.is_valid(model, c2)

    if config.solve
        MOI.optimize!(model)
        @test MOI.get(model, MOI.TerminationStatus()) == MOI.OPTIMAL
        @test MOI.get(model, MOI.PrimalStatus()) == MOI.FEASIBLE_POINT

        @test MOI.get(model, MOI.ResultCount()) >= 1
        @test MOI.get(model, MOI.VariablePrimal(), x1) == 1
        @test MOI.get(model, MOI.VariablePrimal(), x2) == 1
    end
end

const ifthenelsetests = Dict(
    "ifthenelse_singlevariable" => ifthenelse_singlevariable_test,
    "ifthenelse_scalaraffinefunction" =>
        ifthenelse_scalaraffinefunction_test,
)

MOIT.@moitestset ifthenelse


text:
using FieldDicts
using Documenter: doctest
using Test

@test iterate(FieldDict(nothing)) === nothing

@test first(FieldDict((x=1,))) == Pair(:x, 1)

mutable struct Foo
    x::Int
    y::Float64
    z::Foo

    Foo() = new()
end

x = Foo()
d = FieldDict(x)

d[:x] = 1

@test get!(d, :x, 2) == 1
@test get!(d, :z, x).x == 1

x = Foo()
d = FieldDict(x)

d[:x] = 1

@test get!(() -> 2, d, :x) == 1
@test get(() -> 2, d, :x) == 1
@test get(() -> 2, d, :z) == 2
get!(() -> x, d, :z)
@test x.x == 1

doctest(FieldDicts)


text:
Base.round(x::Double64) = round(x, RoundNearest)
Base.round(::Type{T}, x::Double64) where {T<:Integer} = T(round(x, RoundNearest))

function Base.round(x::Double64, ::RoundingMode{:Up})
    return ceil(x)
end
Base.round(::Type{T}, x::Double64, ::RoundingMode{:Up}) where {T<:Integer} = T(round(x, RoundUp))

function Base.round(x::Double64, ::RoundingMode{:Down})
    return floor(x)
end
Base.round(::Type{T}, x::Double64, ::RoundingMode{:Down}) where {T<:Integer} = T(round(x, RoundDown))

function Base.round(x::Double64, ::RoundingMode{:ToZero})
    return signbit(x) ? ceil(x) : floor(x)
end
Base.round(::Type{T}, x::Double64, ::RoundingMode{:ToZero}) where {T<:Integer} = T(round(x, RoundToZero))

function Base.round(x::Double64, ::RoundingMode{:RoundFromZero})
    return signbit(x) ? floor(x) : ceil(x)
end
Base.round(::Type{T}, x::Double64, ::RoundingMode{:FromZero}) where {T<:Integer} = T(round(x, RoundFromZero))

function Base.round(x::Double64, ::RoundingMode{:Nearest})
    signbit(x) && return -round(-x, RoundNearest)
    a = trunc(x + 0.5)
    return iseven(a) ? a : trunc(x - 0.5)
end
Base.round(::Type{T}, x::Double64, ::RoundingMode{:Nearest}) where {T<:Integer} = T(round(x, RoundNearest))

function Base.round(x::Double64, ::RoundingMode{:NearestTiesAway})
    signbit(x) && return -round(-x, RoundNearestTiesAway)
    !isinteger(x - 0.5) && return round(x, RoundNearest)
    return round(x + 0.5, RoundNearest)
end
Base.round(::Type{T}, x::Double64, ::RoundingMode{:NearestTiesAway}) where {T<:Integer} = T(round(x, RoundNearestTiesAway))

function Base.round(x::Double64, ::RoundingMode{:NearestTiesUp})
    signbit(x) && return -round(-x, RoundNearestTiesUp)
    !isinteger(x - 0.5) && return round(x, RoundUp)
    return round(x + 0.5, RoundNearest)
end
Base.round(::Type{T}, x::Double64, ::RoundingMode{:NearestTiesUp}) where {T<:Integer} = T(round(x, RoundNearestTiesUp))

#=
Base.round(x::Double64; digits::Integer=0, base = 10) = Double64(round(Float128(x); digits=digits, base=base))
Base.round(x::Double64; sigdigits::Integer=0, base = 10) = Double64(round(Float128(x); sigdigits=sigdigits, base=base))
Base.round(x::Double64, r::RoundingMode; digits=Integer=0, base = 10) = Double64(round(Float128(x), r; digits=digits, base=base))
Base.round(x::Double64, r::RoundingMode; sigdigits=Integer=0, base = 10) = Double64(round(Float128(x), r; sigdigits=sigdigits, base=base))
=#

#=
# NOTE: this relies on the current keyword dispatch behaviour (#9498).
function Base.round(x::Double64, r::RoundingMode=RoundNearest;
               digits::Union{Nothing,Integer}=nothing, sigdigits::Union{Nothing,Integer}=nothing, base::Union{Nothing,Integer}=nothing)
    if digits === nothing
        if sigdigits === nothing
            if base === nothing
                # avoid recursive calls
                throw(MethodError(round, (x,r)))
            else
                round(x,r)
                # or throw(ArgumentError("`round` cannot use `base` argument without `digits` or `sigdigits` arguments."))
            end
        else
            isfinite(x) || return float(x)
            _round_sigdigits(x, r, sigdigits, base === nothing ? 10 : base)
        end
    else
        if sigdigits === nothing
            isfinite(x) || return float(x)
            _round_digits(x, r, digits, base === nothing ? 10 : base)
        else
            throw(ArgumentError("`round` cannot use both `digits` and `sigdigits` arguments."))
        end
    end
end
=#

function Base._round_sigdigits(x::Double64, r::RoundingMode, sigdigits::Integer, base)
    h = Base.hidigit(x, base)
    Base._round_digits(x, r, sigdigits-h, base)
end

Base.hidigit(x::Double64, base) = Base.hidigit(Hi(x), base)

#=
  round([T,] x, [r::RoundingMode])
  round(x, [r::RoundingMode]; digits::Integer=0, base = 10)
  round(x, [r::RoundingMode]; sigdigits::Integer, base = 10)
  Rounds the number x.
  Without keyword arguments, x is rounded to an integer value, returning a value of type T, or of the same type of x
  if no T is provided. An InexactError will be thrown if the value is not representable by T, similar to convert.
  If the digits keyword argument is provided, it rounds to the specified number of digits after the decimal place (or
  before if negative), in base base.
  If the sigdigits keyword argument is provided, it rounds to the specified number of significant digits, in base
  base.
  The RoundingMode r controls the direction of the rounding; the default is RoundNearest, which rounds to the nearest
  integer, with ties (fractional values of 0.5) being rounded to the nearest even integer. Note that round may give
  incorrect results if the global rounding mode is changed (see rounding).
=#


text:
pushfirst!(LOAD_PATH, joinpath(@__DIR__, "..", "packages"))
import VSCodeLiveUnitTesting
popfirst!(LOAD_PATH)

VSCodeLiveUnitTesting.live_unit_test(ARGS[1], ARGS[2])

max_stars_repo_name: novitk/julia-vscode

max_issues_repo_name: novitk/julia-vscode

max_forks_repo_name: novitk/julia-vscode


text:
function comp_data(
    conn::Union{LibPQ.Connection, DBInterface.Connection},
    gvkeys::AbstractArray{String},
    dateStart::Union{Date,Int}=1950,
    dateEnd::Union{Date,Int}=Dates.today();
    annual::Bool=true,
    filters::Dict{String,<:Any}=Dict(
        "datafmt" => "STD",
        "indfmt" => "INDL",
        "consol" => "C",
        "popsrc" => "D"
    ),
    cols::Array{String}=["gvkey", "datadate", "fyear", "sale", "revt", "xopr"],
)
    if 0 < length(gvkeys) < 1000 # I limit this to 1000 since larger strings
        # slow down the query more than necessary
        temp_filters = Dict{String, Any}()
        for (key, val) in filters
            temp_filters[key] = val
        end
        filters = temp_filters
        filters["gvkey"] = gvkeys
    end

    return comp_data(
            conn,
            dateStart,
            dateEnd;
            annual,
            filters,
            cols
        )
end

"""
    function comp_data(
        conn::Union{LibPQ.Connection, DBInterface.Connection}[,
        gvkeys::AbstractArray{String},]
        dateStart::Union{Date,Int}=1950,
        dateEnd::Union{Date,Int}=Dates.today();
        annual::Bool=true,
        filters::Dict{String,<:Any}=Dict(
            "datafmt" => "STD",
            "indfmt" => "INDL",
            "consol" => "C",
            "popsrc" => "D"
        ),
        cols::Array{String}=["gvkey", "datadate", "fyear", "sale", "revt", "xopr"]
    )

Downloads data from Compustat for firms (if list of gvkeys is provided,
filters to those firms) available firms over a period. Data can be annual
(set `annual=true`) or quarterly (set `annual=false`). For
quarterly data, you also likely need to change the columns
that are downloaded (ie, sales is "saleq" in quarterly data).

Filters is a dictionary of String => String (or array of String) pairings that
will be applied to the SQL query.
"""
function comp_data(
    conn::Union{LibPQ.Connection, DBInterface.Connection},
    dateStart::Union{Date,Int}=1950,
    dateEnd::Union{Date,Int}=Dates.today();
    annual::Bool=true,
    filters::Dict{String,<:Any}=Dict(
        "datafmt" => "STD",
        "indfmt" => "INDL",
        "consol" => "C",
        "popsrc" => "D"
    ),
    cols::AbstractArray{String}=["gvkey", "datadate", "fyear", "sale", "revt", "xopr"]
)

    if typeof(dateStart) == Int
        dateStart = Dates.Date(dateStart, 1, 1)
    end
    if typeof(dateEnd) == Int
        dateEnd = Dates.Date(dateEnd, 12, 31)
    end

    tab = annual ? default_tables["comp_funda"] : default_tables["comp_fundq"]

    colString = join(cols, ", ")
    filterString = create_filter(filters, "WHERE datadate between '$dateStart' AND '$dateEnd'")
    query = """
        select $colString
        from $tab
        $filterString
    """

    return run_sql_query(conn, query) |> DataFrame
end


text:
# Calculate various quantities for the Gaussian task (2.3.1)

calcχ(mu, ν, sigma) = mu * ν / sigma

calcν(sigma, sigma_n) = sigma^2 / sigma_n^2

calch_n(learnerG::GLearner, y) = exp(-y^2/(2 * learnerG.sigma^2)) / sqrt(2. * π * learnerG.sigma^2)

calcf(χ, ν) = exp(-χ^2/(2. * ν)) / sqrt(2. * π / ν)

calcφ(learnerG::Union{GLearner, GBayesFilter, GMPN, GSOR}, y) = y / learnerG.sigma

calcfratio(χnom, νnom, χdenom, νdenom) = exp(- (χnom^2)/(2. * νnom) + (χdenom^2)/(2. * νdenom)) * sqrt(νnom / νdenom)

calcsqrtterm(νnom1, νnom2, νdenom1, νdenom2) = sqrt((νnom1 * νnom2) / (νdenom1 * νdenom2))

function calcexpterm(χnom1, χnom2, χdenom1, χdenom2,
                    νnom1, νnom2, νdenom1, νdenom2)

    expterm = - (χnom1^2) / (2. * νnom1)
    expterm -= (χnom2^2) / (2. * νnom2)
    expterm += (χdenom1^2) / (2. * νdenom1)
    expterm += (χdenom2^2) / (2. * νdenom2)
end
function calcdoublefratio(χnom1, χnom2, χdenom1, χdenom2,
                    νnom1, νnom2, νdenom1, νdenom2)

    expterm = calcexpterm(χnom1, χnom2, χdenom1, χdenom2,
                        νnom1, νnom2, νdenom1, νdenom2)
    sqrtterm = calcsqrtterm(νnom1, νnom2, νdenom1, νdenom2)
    doublefratio = exp(expterm) * sqrtterm
end
function checkInf(variable)
    if isinf(variable)
        #println("I'm Inf")
        variable = 1e16
    end
    variable
end
function computemu_n!(learnerG::GLearner)
    learnerG.mu_n[1] = learnerG.sigma * learnerG.χ_n[1] / learnerG.ν_n[1]
end
# Calculate Bayes Factor Surprise
function calcSgm!(learnerG::GLearner, y)
    learnerG.Sgm_n[1] = calcdoublefratio(learnerG.χ_0,
                        learnerG.χ_n[1] + calcφ(learnerG, y),
                        learnerG.χ_0 + calcφ(learnerG, y),
                        learnerG.χ_n[1],
                        learnerG.ν_0, learnerG.ν_n[1] + 1.,
                        learnerG.ν_0 + 1., learnerG.ν_n[1])
    learnerG.Sgm_n[1] = checkInf(learnerG.Sgm_n[1])
end
function computevar_n!(learnerG::GLearner)
    learnerG.var_n[1] = learnerG.sigma / learnerG.ν_n[1]
end
function calcγ!(learnerG::Union{GLearner, GBayesFilter, GMPN, GSOR})
    learnerG.γ_n[1] = learnerG.m * learnerG.Sgm_n[1]/(1. + learnerG.m * learnerG.Sgm_n[1])
end
# Calculate Shannon Surprise
function calcShannon!(learnerG::Union{GLearner, GBayesFilter, GMPN, GSOR}, y)
    # Py_0 = ~N(y; m_0, σ0^2 + σ^2)
    expterm = exp(- (y^2)/(2. * learnerG.sigma^2) -
                (learnerG.χ_0^2)/(2. * learnerG.ν_0) +
                (learnerG.χ_0 + calcφ(learnerG, y))^2/(2. * (learnerG.ν_0 + 1.))
                )
    sqrtterm = sqrt(learnerG.ν_0 /
                    (2. * pi * (learnerG.ν_0 + 1.))
                    )
    Py_0 = expterm * sqrtterm / learnerG.sigma
    Py_0 = checkInf(Py_0)
    learnerG.Ssh_n[1] = - log10(Py_0) + log10(1. /
                                        (learnerG.changeprobability +
                                        ((1. - learnerG.changeprobability) / learnerG.Sgm_n[1])
                                        )
                                    )
    learnerG.Ssh_n[1] = checkInf(learnerG.Ssh_n[1])
end


text:
################################################################################
#
#     Hecke.jl : Hecke main file
#
# This file is part of Hecke.
#
# Copyright (c) 2015: Claus Fieker, Tommy Hofmann
# All rights reserved.
#
# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided that the following conditions are met:
# * Redistributions of source code must retain the above copyright notice, this
#   list of conditions and the following disclaimer.
#
# * Redistributions in binary form must reproduce the above copyright notice,
#   this list of conditions and the following disclaimer in the documentation
#   and/or other materials provided with the distribution.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
# DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
# FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
# DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
# SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
# CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
# OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
# OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
#
# (C) 2015, 2016 Claus Fieker, Tommy Hofmann
#
################################################################################

__precompile__()

module Hecke

################################################################################
#
#  Load FPlll if available
#
################################################################################

if isdir(joinpath(Pkg.dir(),"FPlll"))
  using FPlll
end

################################################################################
#
#  Import/export
#
################################################################################

import Base: show, minimum, rand, prod, copy, rand!, rand, ceil, round,
             size, dot, in, powermod, ^, getindex, ==, <, >, +, *, /, \, -, !=,
             getindex, setindex!, transpose, getindex, //, colon, div,
             floor, max, BigFloat, precision, dot,
             first, StepRange, show, one, zero, inv, iseven, isodd, convert,
             angle, abs2, isless, exponent, base, isfinite, zeros, rem,
             min, numerator, denominator, exp, maximum

# To make all exported Nemo functions visible to someone using "using Hecke"
# we have to export everything again
# dong it the "import" route, we can pick & choose...

import AbstractAlgebra

import Nemo

exclude = [:Nemo, :AbstractAlgebra, :RealField,
           :factor,
           :call, :factors, :parseint, :strongequal, :window, :xgcd]

for i in names(Nemo)
  i in exclude && continue
  eval(parse("import Nemo." * string(i)))
  eval(Expr(:export, i))
end

import Nemo: acb_struct, Ring, Group, Field, NmodRing, nmod, arf_struct,
       elem_to_mat_row!, elem_from_mat_row,
       acb_vec, array, acb_vec_clear

export @vprint, @hassert, @vtime, add_verbose_scope, get_verbose_level,
       set_verbose_level, add_assert_scope, get_assert_level, set_assert_level,
       update, show, StepRange, domain, codomain, image, preimage,
       modord, resultant, @test_and_infer, next_prime, ispower

###############################################################################
#
#   Library initialisation
#
###############################################################################

#const pkgdir = realpath(joinpath(dirname(@__FILE__), ".."))
const pkgdir = const pkgdir = Pkg.dir("Hecke")
const libhecke = joinpath(pkgdir, "local", "lib", "libhecke")
const libdir = joinpath(pkgdir, "local", "lib")

function __init__()

  if myid() == 1
    println("")
    print("Welcome to \n")
    print_with_color(:red, "
    _    _           _
   | |  | |         | |
   | |__| | ___  ___| | _____
   |  __  |/ _ \\/ __| |/ / _ \\
   | |  | |  __/ (__|   <  __/
   |_|  |_|\\___|\\___|_|\\_\\___|
    ")

    println()
    print("Version")
    print_with_color(:green, " $VERSION_NUMBER ")
    print("... \n ... which comes with absolutely no warranty whatsoever")
    println()
    println("(c) 2015-2018 by Claus Fieker, Tommy Hofmann and Carlo Sircana")
    println()
  else
    println("Hecke $VERSION_NUMBER ...")
  end

  if "HOSTNAME" in keys(ENV) && ENV["HOSTNAME"] == "juliabox"
    push!(Libdl.DL_LOAD_PATH, "/usr/local/lib")
  elseif is_linux()
    push!(Libdl.DL_LOAD_PATH, libdir)
    Libdl.dlopen(libhecke)
  else
    push!(Libdl.DL_LOAD_PATH, libdir)
    ENV["PATH"] = ENV["PATH"] * ";" * joinpath(Pkg.dir("Nemo"), "local", "lib")
  end

  t = create_accessors(AnticNumberField, acb_root_ctx, get_handle())
  global _get_nf_conjugate_data_arb = t[1]
  global _set_nf_conjugate_data_arb = t[2]

  t = create_accessors(AnticNumberField, Dict{Int, acb_roots}, get_handle())
  global _get_nf_conjugate_data_arb_roots = t[1]
  global _set_nf_conjugate_data_arb_roots = t[2]

  t = create_accessors(AnticNumberField,
                       Tuple{Array{nf_elem, 1}, nf_elem},
                       get_handle())
  global _get_nf_torsion_units = t[1]
  global _set_nf_torsion_units = t[2]

  t = create_accessors(AnticNumberField, NfOrd, get_handle())

  global _get_maximal_order_of_nf = t[1]
  global _set_maximal_order_of_nf = t[2]

  t = create_accessors(NfOrd, ClassGrpCtx, get_handle())

  global _get_ClassGrpCtx_of_order = t[1]
  global _set_ClassGrpCtx_of_order = t[2]

  t = create_accessors(NfOrd, UnitGrpCtx, get_handle())

  global _get_UnitGrpCtx_of_order = t[1]
  global _set_UnitGrpCtx_of_order = t[2]

  t = create_accessors(NfOrd, Array, get_handle())

  global _get_carlos_units_of_order = t[1]
  global _set_carlos_units_of_order = t[2]

  t = create_accessors(AnticNumberField, roots_ctx, get_handle())

  global _get_roots_ctx_of_nf = t[1]
  global _set_roots_ctx_of_nf = t[2]

  t = create_accessors(AnticNumberField, Array, get_handle())

  global _get_cyclotomic_ext_nf = t[1]
  global _set_cyclotomic_ext_nf = t[2]

  t = create_accessors(RelativeExtension, NfRelOrd, get_handle())

  global _get_maximal_order_of_nf_rel = t[1]
  global _set_maximal_order_of_nf_rel = t[2]

  global R = _RealRing()

  # Stuff for elliptic curves
  # polynomial rings Zx = ZZ[x] and _Zxy = ZZ[x,y]
  # will be removed eventually
  global const _Zx = PolynomialRing(FlintZZ, "_x")[1]
  global const _Zxy = PolynomialRing(_Zx, "_y")[1]
  global const _x = gen(_Zx)
  global const _y = gen(_Zxy)

  global const flint_rand_ctx = flint_rand_state()
end

function conjugate_data_arb(K::AnticNumberField)
  try
    c = _get_nf_conjugate_data_arb(K)::acb_root_ctx
    return c
  catch
    c = acb_root_ctx(K.pol)
    _set_nf_conjugate_data_arb(K, c)
    return c
  end
end

function conjugate_data_arb_roots(K::AnticNumberField, p::Int)
  already_set = false
  local c
  try
    c = _get_nf_conjugate_data_arb_roots(K)::Dict{Int, acb_roots}
    already_set = true
  catch
    c = Dict{Int, acb_roots}()
  end

  if already_set && haskey(c, p)
    return c[p]
  end

  #if p > 2^18
  #  Base.show_backtrace(STDOUT, backtrace())
  #end
  rootc = conjugate_data_arb(K)
  q = rootc.prec
  while q < p
    refine(rootc)
    q = rootc.prec
  end
  @assert p <= q
  rall = deepcopy(rootc.roots)
  rreal = deepcopy(rootc.real_roots)
  rcomplex = deepcopy(rootc.complex_roots)
  for z in rall
    expand!(z, -p)
  end
  for z in rreal
    expand!(z, -p)
  end
  for z in rcomplex
    expand!(z, -p)
  end
  c[p] = acb_roots(p, rall, rreal, rcomplex)
  if !already_set
    _set_nf_conjugate_data_arb_roots(K, c)
  end
  return c[p]::acb_roots
end

function _torsion_units(K::AnticNumberField)
  try
    c = _get_nf_torsion_units(K)::Tuple{Array{nf_elem, 1}, nf_elem}
    return c
  catch
    O = maximal_order(K)
    tor, gen = _torsion_units(O)
    tornf = [ elem_in_nf(x) for x in tor]
    gennf = elem_in_nf(gen)
    _set_nf_torsion_units(K, (tornf, gennf))
    return tornf, gennf
  end
end

function torsion_units(K::AnticNumberField)
  ar, g = _torsion_units(K)
  return ar
end

function torsion_units_gen(K::AnticNumberField)
  ar, g = _torsion_units(K)
  return g
end

################################################################################
#
#  Version number
#
################################################################################

global VERSION_NUMBER = v"0.4.5"

################################################################################
#
#  Verbose printing
#
################################################################################

global hecke = Hecke

global VERBOSE_SCOPE = Symbol[]

global VERBOSE_LOOKUP = Dict{Symbol, Int}()

global VERBOSE_PRINT_INDENT = [ 0 ]

function add_verbose_scope(s::Symbol)
  !(s in VERBOSE_SCOPE) && push!(VERBOSE_SCOPE, s)
  nothing
end

function pushindent()
  a = VERBOSE_PRINT_INDENT[1]
  VERBOSE_PRINT_INDENT[1] = a + 1
  nothing
end

function clearindent()
  VERBOSE_PRINT_INDENT[1] = 0
  nothing
end

function popindent()
  a = VERBOSE_PRINT_INDENT[1]
  VERBOSE_PRINT_INDENT[1] = a - 1
  nothing
end

function _global_indent()
  s = "  "^VERBOSE_PRINT_INDENT[1]
  return s
end

macro vprint(args...)
  if length(args) == 2
    quote
      if get_verbose_level($(args[1])) >= 1
        print(_global_indent())
        print($(esc((args[2]))))
      end
    end
  elseif length(args) == 3
    quote
      if get_verbose_level($(args[1])) >= $(args[2])
        print(_global_indent())
        print($(esc((args[3]))))
      end
    end
  end
end

macro v_do(args...)
  if length(args) == 2
    quote
      if get_verbose_level($(args[1])) >= 1
       $(esc(args[2]))
      end
    end
  elseif length(args) == 3
    quote
      if get_verbose_level($(args[1])) >= $(args[2])
        $(esc(args[3]))
      end
    end
  end
end

function set_verbose_level(s::Symbol, l::Int)
  !(s in VERBOSE_SCOPE) && error("Not a valid symbol")
  VERBOSE_LOOKUP[s] = l
end

function get_verbose_level(s::Symbol)
  !(s in VERBOSE_SCOPE) && error("Not a valid symbol")
  return get(VERBOSE_LOOKUP, s, 0)::Int
end

################################################################################
#
#  Assertions
#
################################################################################

global ASSERT_SCOPE = Symbol[]

global ASSERT_LOOKUP = Dict{Symbol, Int}()

function add_assert_scope(s::Symbol)
  !(s in ASSERT_SCOPE) && push!(ASSERT_SCOPE, s)
  nothing
end

function set_assert_level(s::Symbol, l::Int)
  !(s in ASSERT_SCOPE) && error("Not a valid symbol")
  ASSERT_LOOKUP[s] = l
end

function get_assert_level(s::Symbol)
  !(s in ASSERT_SCOPE) && error("Not a valid symbol")
  return get(ASSERT_LOOKUP, s, 0)::Int
end

macro hassert(args...)
  if length(args) == 2
    quote
      if get_assert_level($(args[1])) >= 1
        @assert $(esc(args[2]))
      end
    end
  elseif length(args) == 3
    quote
      if get_assert_level($(args[1])) >= $(args[2])
        @assert $(esc(args[3]))
      end
    end
  end
end

################################################################################
#
#  Custom test function
#
################################################################################

function test_module(x, new::Bool = true)
   julia_exe = Base.julia_cmd()
   if x == "all"
     test_file = joinpath(pkgdir, "test/runtests.jl")
   else
     test_file = joinpath(pkgdir, "test/$x.jl")
   end

   if new
     cmd = "using Base.Test; using Hecke; include(\"$test_file\");"
     info("spawning ", `$julia_exe -e \"$cmd\"`)
     run(`$julia_exe -e $cmd`)
   else
     info("Running tests for $x in same session")
     try
       include(test_file)
     catch e
       if isa(e, LoadError)
         println("You need to do \"using Base.Test\"")
       else
         rethrow(e)
       end
     end
   end
end

################################################################################
#
#   Do @infert and @test simultanously
#
################################################################################

macro test_and_infer(f,args,res)
  quote
    if isa($(esc(args)), Tuple)
      Base.Test.@inferred $(esc(f))($(esc(args))...)
      Base.Test.@test $(esc(f))($(esc(args))...) == $(esc(res))
    else
      Base.Test.@inferred $(esc(f))($(esc(args)))
      local t = $(esc(res))
      Base.Test.@test $(esc(f))($(esc(args))) == t
    end
  end
end

################################################################################
#
#  Verbose time printing
#
################################################################################

macro vtime(args...)
  if length(args) == 2
    msg = string(args[2])
    quote
      if get_verbose_level($(args[1])) >= 1
        local t0 = time_ns()
        local val = $(esc(args[2]))
        println((time_ns()-t0)/1e9, " @ ", $msg)
        val
      else
        local val2 = $(esc(args[2]))
        val2
      end
    end
  elseif length(args) == 3
    msg = string(args[3])
    quote
      if get_verbose_level($(args[1])) >= $(args[2])
        local t0 = time_ns()
        local val = $(esc(args[3]))
        println((time_ns()-t0)/1e9, " @ ", $msg)
        val
      else
        local val2 = $(esc(args[3]))
        val2
      end
    end
  end
end

#usage
# @vtime_add_ellapsed :ClassGroup 2 clg :saturate  s= hnf(a)
# @vtime_add :ClassGroup 2 clg :saturate  0.5
# -> clg.time[:saturate] +=
function _vtime_add(D::Dict, k::Any, v::Any)
  if haskey(D, k)
    D[k] += v
  else
    D[k] = v
  end
end

macro vtime_add(flag, level, var, key, value)
  quote
    if get_verbose_level($flag) >= $level
      _vtime_add($(esc(var)).time, $key, $(esc(value)))
    end
  end
end

macro vtime_add_elapsed(flag, level, var, key, stmt)
  quote
    tm = @elapsed $(esc(stmt))
    if get_verbose_level($flag) >= $level
      _vtime_add($(esc(var)).time, $key, tm)
    end
  end
end

################################################################################
#
#  Exception types
#
################################################################################

mutable struct LowPrecisionCholesky <: Exception end

Base.showerror(io::IO, e::LowPrecisionCholesky) =
    print(io, """
    Negative diagonal in Cholesky decomposition, probably a precision issue""")

mutable struct LowPrecisionLLL <: Exception end

Base.showerror(io::IO, e::LowPrecisionLLL) =
    print(io, """
    Transformation matrix has too large entries relative to precision in LLL""")

# what is this function doing here?
function checkbounds(a::Int, b::Int) nothing; end;

################################################################################
add_assert_scope(:PID_Test)
set_assert_level(:PID_Test, 0)

################################################################################
#
#  "Submodules"
#
################################################################################

include("HeckeTypes.jl")
include("NfRel/NfRelTypes.jl")
include("AlgAss/AlgAssTypes.jl")
include("Map.jl")
include("Misc.jl")
include("GrpAb.jl")
include("LinearAlgebra.jl")
include("NfAbs.jl")
include("NfOrd.jl")
include("Sparse.jl")
include("BigComplex.jl")
include("conjugates.jl")
include("NfRel.jl")
include("analytic.jl")
include("helper.jl")
include("EllCrv.jl")
include("LargeField.jl")
include("RCF.jl")
include("Grp.jl")
include("AlgAss.jl")

for T in subtypes(Map(HeckeMap))
  (M::T)(a) = image(M, a)
end

################################################################################
#
#  Verbose printing
#
################################################################################

doc"""
***
    vshow(A) -> Void

> Prints all fields of $A$.
"""
function vshow(A)
  for i in fieldnames(typeof(A))
    if isdefined(A, i)
      println("$i: ")
      println(getfield(A, i), "\n")
    else
      println("$i: Not definied")
    end
  end
end
################################################################################
#
#  Element types for parent types
#
################################################################################

# Nemo only provides element_types for parent objects

elem_type{T}(::Type{FacElemMon{T}}) = FacElem{elem_type(T), T}

elem_type{T}(::Type{Generic.ResRing{T}}) = Generic.Res{T}

################################################################################
#
#  Aliases
#
################################################################################

hasroot = ispower

################################################################################
#
#  An update function
#
################################################################################

function update()

  olddir = pwd()

  println("Updating Hecke ... ")
  cd(Pkg.dir("Hecke"))
  run(`git pull`)

  pkgdir = Pkg.dir("Nemo")

  println("Updating Nemo ... ")
  cd("$pkgdir")
  run(`git pull`)

  println("Updating antic ... ")
  cd("$pkgdir/deps/antic")
  run(`git pull`)

  println("Updating flint ... ")
  cd("$pkgdir/deps/flint2")
  run(`git pull`)
  run(`make -j`)
  run(`make install`)

  println("Updating arb ... ")
  cd("$pkgdir/deps/arb")
  run(`git pull`)
  run(`make -j`)
  run(`make install`)

  cd(olddir)
end

function whos(io::IO=STDOUT, m::Module=current_module(), pattern::Regex=r"")
    maxline = Base.tty_size()[2]
    line = zeros(UInt8, maxline)
    head = PipeBuffer(maxline + 1)
    for v in sort!(names(m, true)) # show also NON exported stuff!
        s = string(v)
        if isdefined(m, v) && ismatch(pattern, s)
            value = getfield(m, v)
            @printf head "%30s " s
            try
                bytes = Base.summarysize(value)
                if bytes < 10_000
                    @printf(head, "%6d bytes  ", bytes)
                else
                    @printf(head, "%6d KB     ", bytes ÷ (1024))
                end
                print(head, Base.summary(value))
            catch e
                print(head, "#=ERROR: unable to show value=#")
                println(e)
            end

            newline = search(head, UInt8('\n')) - 1
            if newline < 0
                newline = nb_available(head)
            end
            if newline > maxline
                newline = maxline - 1 # make space for ...
            end
            line = resize!(line, newline)
            line = read!(head, line)

            Base.write(io, line)
            if nb_available(head) > 0 # more to read? replace with ...
                print(io, '\u2026') # hdots
            end
            println(io)
            seekend(head) # skip the rest of the text
        end
    end
end

whos(m::Module, pat::Regex=r"") = whos(STDOUT, m, pat)
whos(pat::Regex) = whos(STDOUT, current_module(), pat)

function print_cache()
  sym = [];
  for a in collect(names(Nemo, true));
    d = parse("Nemo." * string(a));
      try z = eval(d); push!(sym, (d, z));
    catch e;
    end;
  end
  for f in sym;
    #if f[2] isa Array || f[2] isa Dict || f[2] isa ObjectIdDict;
    try
      print(f[1], " ", length(f[2]), "\n");
    end;
  end

  sym = [];
  for a in collect(names(Nemo.Generic, true));
    d = parse("Nemo.Generic." * string(a));
      try z = eval(d); push!(sym, (d, z));
    catch e;
    end;
  end
  for f in sym;
    #if f[2] isa Array || f[2] isa Dict || f[2] isa ObjectIdDict;
    try
      print(f[1], " ", length(f[2]), "\n");
    catch e
    end;
  end

  sym = [];
  for a in collect(names(Hecke, true));
    d = parse("Hecke." * string(a));
      try z = eval(d); push!(sym, (d, z));
    catch e;
    end;
  end
  for f in sym;
    #if f[2] isa Array || f[2] isa Dict || f[2] isa ObjectIdDict;
    try
      print(f[1], " ", length(f[2]), "\n");
    end;
  end

end

################################################################################
#
#  Testing only "submodules"
#
################################################################################

#
# stuff for 0.5
#

@inline __get_rounding_mode() = Base.MPFR.rounding_raw(BigFloat)

#precompile(maximal_order, (AnticNumberField, ))
#precompile(class_group, (NfOrd, ))

end


text:
using Peaks
using Glob
using FFTW
using DSP
using LsqFit

"""
	SPeaks

A struct to store the peaks found by peak finding algorithm

# Fields

- `peaks::Vector{Int64}`          : indexes of the peaks
- `proms::Vector{Float64}`        : prominence of the peaks
- `widths::Vector{Float64}`       : width of the peaks
- `leftedge::Vector{Float64}`     : position of the left edge
- `rightedge::Vector{Float64}`    : position of the right edge
- `xs::Vector{Float64}`           : vector of `x` positions (times)
- `ys::Vector{Float64}`           : vector of `y` positions (voltages)
- `xbase::Vector{Float64}`        : pedestal
- `promsel::Float64`              : selection value of prmominence

"""
struct SPeaks
	peaks::Vector{Int64}
	proms::Vector{Float64}
	widths::Vector{Float64}
	leftedge::Vector{Float64}
	rightedge::Vector{Float64}
	xs::Vector{Float64}
	ys::Vector{Float64}
	xbase::Vector{Float64}
	promsel::Float64
end

"""
	read_waveform(wdir::String, wname::String; wmus::Float64=5, symmetric::Bool=false)

reads a weveform specifying directory and file name

# Fields

- `wdir`   : name of directory
- `wname`  : name of file
- `wmus`   : width of time window in μs.
- `wsym`   : If true the window is symmetric (centered at zero)
- `sign`   : reverses (the default) or not the waveform sign
- `vf`     : voltage factor, (the default is volts to mlVolts, 1000 )
- `tf`     : time factor, (the default is seconds to microseconds, 10^6 )

# Returns

Waveform as a data frame. Defaults: amplitude in mV (signal positive), time in μs.

"""
function read_waveform(wdir::String, wname::String; wmus::Float64, wsym::Bool=true,
	                   sign=-1.0, vf= 1e+3, tf=1e+6)

	wvfm = load_df_from_csv(wdir, wname, enG,5)
	shift_waveform!(wvfm, wmus, wsym, sign, vf, tf)

end

"""
	read_waveform(csvf::Vector{String}, fnumber::Integer; wmus::Float64, wsym::Bool)

  Like previous method but takes a vector of strings with the name of the files
  and an integer with the index of the file to be read

"""
function read_waveform(csvf::Vector{String}, fnumber::Integer; wmus::Float64,
	                   wsym::Bool)

	wvfm = DataFrame(CSV.File(csvf[fnumber], header=5, delim=enG.delim, decimal=enG.decimal))
	shift_waveform!(wvfm, wmus, wsym)
end

"""
	shift_waveform!(wvfm:DataFrame, wmus::Float64, wsym::Bool)

 Shifts waveform to mV in Y (reverses sign) and to μs in x.
 If the window is centered in zero, shifts it, so that first time is zero.

 # Fields

- `wvfm`   : Waveform (a Data Frame with time/amplitude columns)
- `wmus`   : width of time window in μs.
- `wsym`   : If true the window is symmetric (centered at zero)
- `sign`   : reverses (the default) or not the waveform sign
- `vf`     : voltage factor, (the default is volts to mlVolts, 1000 )
- `tf`     : time factor, (the default is seconds to microseconds, 10^6 )

# Returns

Shifted waveform

"""
function shift_waveform!(wvfm::DataFrame, wmus::Float64, wsym::Bool, sign=-1.0, vf= 1e+3, tf=1e+6)
	wvfm[!, :Ampl] = sign * vf * wvfm[!, :Ampl]
	wvfm[!, :Time] = tf  .* wvfm[!, :Time]
	if wsym
		wvfm[!, :Time] = wmus/2.0  .+ wvfm[!, :Time]
	end
	wvfm
end

"""
		sampling_rate(wvfm::DataFrame, wmus::Float64)

# Fields

- `wvfm`  : Waveform (a Data Frame with time/amplitude columns)
- `wmus`  : width of time window in μs.

# Returns

The sampling rate

"""
function sampling_rate(wvfm::DataFrame, wmus::Float64)
	uconvert(GHz, (length(wvfm.Time) / (wmus*1.0e-6)) * Hz)
end

"""
sampling_period(wvfm::DataFrame, wmus::Float64)

# Fields

- `wvfm`  : Waveform (a Data Frame with time/amplitude columns)
- `wmus`  : width of time window in μs.

# Returns

The sampling period

"""
function sampling_period(wvfm::DataFrame, wmus::Float64)
	uconvert(ns, 1.0 / sampling_rate(wvfm, wmus))
end

"""
	wstats(wvfm::DataFrame; nsigma=1, waveform='filtered')

Return the mean, std, and (+-) threshold at nsigma, for raw or filtered waveform

"""
function wstats(wvfm::DataFrame; nsigma=1, waveform="filtered")
	if waveform=="filtered"
		ampl = "fAmpl"
	else
		ampl = "Ampl"
	end

	meanfs = mean(wvfm[!, ampl])
	stdfs = std(wvfm[!, ampl])
	thrp = meanfs + stdfs * nsigma
	thrn = meanfs - stdfs * nsigma
	meanfs, stdfs, thrp, thrn
end

"""
	fourier_transform(wvfm1::DataFrame, ws::Float64)

	Computes the fourier transform of wvfm1

	# Fields

- `wvfm1`   : Waveform
- `ws`      : width of the data window in μs

"""
function fourier_transform(wvfm1::DataFrame, ws::Float64)
	ts = sampling_period(wvfm1, ws)
	tss = uconvert(s, ts)/s           # tss is the sampling period in seconds
	F = fft(wvfm1.Ampl) |> fftshift   # fourier transform
	freqs = fftfreq(length(wvfm1.Time), 1.0/tss) |> fftshift  # frequencies in Hz
	ts, F, freqs
end

"""
	filter_signal_lowpass(wvfm::DataFrame,
	                      flh::Float64, fs::Float64, filtertype::String,
	                      n=4, ϵ=1.0)

`flh`         : higer frequency accepted by the filter in Hz
`fs`          : sampling frequency
`filtertype`  : two possible filters: Butterworth or Chebyshev1
`n`           : filter order
`ϵ`           : ripple factor for Chebyshev1 filter

Returns

frequency response, filter (Butterworth, or Chebyshev) and filtered signal

"""
function filter_signal_lowpass(wvfm::DataFrame,
	                           flh::Float64, fs::Float64, filtertype::String,
	                           n=4, ϵ=1.0)

	responsetype = Lowpass(flh; fs=fs)

	if filtertype == "Butterworth"
		designmethod = Butterworth(n)
	elseif filtertype == "Chebyshev1"
		ripple = 10.0*log10(ϵ^2 + 1.0)
		designmethod = Chebyshev1(n, ripple)
	end
	filter = digitalfilter(responsetype, designmethod)
	H, _ = freqresp(filter)
	H, filter, filt(filter, wvfm.Ampl)
end

"""
	filter_signal_lp!(wvfm::DataFrame, wmus::Float64;
                           filtertype::String, flhz::Float64,
                           n=4, ϵ=1.0)
Fields

`wvfm`        : Waveform (a Data Frame with time/amplitude columns)
`wmus`        : width of time window in μs.
`filtertype`  : two possible filters: Butterworth or Chebyshev1
`flhz`        : higer frequency accepted by the filter in Hz
`n`           : filter order
`ϵ`           : ripple factor for Chebyshev1 filter

Returns

A column with the filtered amplitude is added to the waveform
"""
function filter_signal_lp!(wvfm::DataFrame, wmus::Float64;
                           filtertype::String, flhz::Float64,
                           n=4, ϵ=1.0)
	fs = uconvert(Hz, sampling_rate(wvfm, wmus))/Hz
	_, _, fsgncb = filter_signal_lowpass(wvfm, flhz, fs, filtertype, n, ϵ)
	wvfm[!, "fAmpl"] = fsgncb
end

"""
	select_filtered_peaks(wvfm::DataFrame, thrp::Float64;
                          promsel::Float64, wsel=0.0)

	Fields

`wvfm`        : Waveform (a Data Frame with time/amplitude columns)
`thrp`        : threshold for peak search.
`promsel`     : cut on peak prominence
`wsel`        : cut on peak width

Returns

A filtered waveform and a SPeaks struct

"""
function select_filtered_peaks(wvfm::DataFrame, thrp::Float64;
                               promsel::Number, wsel=0.0)

	wvfmFlt = filter(row -> row[:fAmpl] > thrp, wvfm)
	pks, _ = findmaxima(wvfmFlt.fAmpl)
	peaks2, proms = peakproms(pks, wvfmFlt.fAmpl; minprom=promsel)

	if length(peaks2) == 0
		return Nothing
	end

	peaks2, widths, leftedge, rightedge = peakwidths(peaks2, wvfmFlt.fAmpl, proms; minwidth=wsel)

	ys = [wvfmFlt.fAmpl[i] for i in peaks2]
	xs = [wvfmFlt.Time[i] for i in peaks2]
	npeaks = length(peaks2)
	xbase = ones(npeaks,1)*[promsel]

	return wvfmFlt, SPeaks(peaks2, proms, widths, leftedge, rightedge, xs,ys, xbase, promsel)
end

"""
	select_peaks(csvf::Vector{String}, i0::Integer, il::Integer;
	             wmus::Float64, wsym::Bool, flhz::Float64, promsel::Number,
	             nsigma=2.0, filtertype="Butterworth")

Fields

`csvf`        : A vector of strings with the names of the files to be readout
`i0`          : index of first file to read
`il`          : index of last file to read
`wmus`        : width of time window in μs.
`wsym`        : If true the window is symmetric (centered at zero)
`flhz`        : higer frequency accepted by the filter in Hz
`promsel`     : cut on peak prominence
`nsigma`      : number of sigmas to compute stats
`filtertype`  : filter type

Returns

A vector of SPeaks

"""
function select_peaks(csvf::Vector{String}, i0::Integer, il::Integer;
	                  wmus::Float64, wsym::Bool, flhz::Float64, promsel::Number,
	                  nsigma=2.0, filtertype="Butterworth")

	SPK = []
	for fileNumber in i0:il
		wvfm = read_waveform(csvf,fileNumber; wmus=wmus, wsym=wsym)
		filter_signal_lp!(wvfm, wmus; filtertype=filtertype, flhz=flhz, n=4)
		_, _, thrp, _ = wstats(wvfm; nsigma=nsigma, waveform="filtered")
		result = select_filtered_peaks(wvfm, thrp; promsel = promsel, wsel=0.0)

		if result != Nothing
			_, speaks = result
			push!(SPK,speaks)
		end
	end
	SPK
end

"""
	proms(csvf::Vector{String}, fi::Integer, fe::Integer;
	      wmus::Float64, wsym::Bool, flhz::Float64, promsel::Number,
	      nsigma=2.0, filtertype="Butterworth")

Fields

`csvf`        : A vector of strings with the names of the files to be readout
`i0`          : index of first file to read
`il`          : index of last file to read
`wmus`        : width of time window in μs.
`wsym`        : If true the window is symmetric (centered at zero)
`flhz`        : higer frequency accepted by the filter in Hz
`promsel`     : cut on peak prominence
`nsigma`      : number of sigmas to compute stats
`filtertype`  : filter type

Returns

A vector of proms and a vector of widths

"""
function proms(csvf::Vector{String}, fi::Integer, fe::Integer;
	wmus::Float64, wsym::Bool, flhz::Float64, promsel::Number,
	nsigma=2.0, filtertype="Butterworth")

	spks0 = select_peaks(csvf, fi, fe; wmus=wmus, wsym=wsym, flhz=flhz,
				promsel=promsel, nsigma=nsigma, filtertype=filtertype)

	sP = reduce(vcat,[spks0[i].proms for i in 1:length(spks0)])
	sW = reduce(vcat,[spks0[i].widths for i in 1:length(spks0)])
	return sP, sW
end

"""
	runproms(csvf::Vector{String}, fi::Integer, fe::Integer;
             wmus::Float64, wsym::Bool=true, flhz::Float64, promsel::Number)

Fields

`csvf`        : A vector of strings with the names of the files to be readout
`fi`          : index of first file to read
`fe`          : index of last file to read
`wmus`        : width of time window in μs.
`wsym`        : If true the window is symmetric (centered at zero)
`flhz`        : higer frequency accepted by the filter in Hz
`promsel`     : cut on peak prominence

Returns

A vector of proms and a vector of widths

"""
function runsel(csvf::Vector{String}, fi::Integer, fe::Integer;
	            wmus::Float64, wsym::Bool, flhz::Float64, promsel::Number, tlmus::Float64)

	spks = select_peaks(csvf, fi, fe; wmus=wmus, wsym=wsym, flhz=flhz, promsel=promsel)
	df = DataFrame()
	df.pr = reduce(vcat,[spks[i].proms  for i in 1:length(spks)])
	df.ws    = reduce(vcat,[spks[i].widths  for i in 1:length(spks)])
	df.ts   = reduce(vcat,[subtract_laser(spks[i].xs, tlmus) for i in 1:length(spks)])
	df.vs   = reduce(vcat,[spks[i].ys  for i in 1:length(spks)])
	df
end

"""
	function subtract_laser(tv::Vector{Float64}, ltus::Float64)

Fields

`tv`   : A vector measured times
`ltus` : time (period) of laser, in μs

Returns

A vector of times where the time of the laser has been subtracted.

"""
function subtract_laser(tv::Vector{Float64}, ltus::Float64)
	tv - ltus * floor.(tv/ltus)
end

"""
	fit_peaks(htime; pa0=[100.0, 0.5], i0=1)
Fields

	`htime` : An histogram of times
	`pa0`   : initial value of fit parameters
	`i0`    : first bin to fit

Returns

Fits an exponential to the histogramm, returns the coefficients
of the fit, the errors of the fit and a vector of predictions from the fit
"""
function fit_peaks(htime; pa0=[100.0, 0.5], i0=1)
	expo(t, N, λ) = N*exp(-t/λ)
	mexp(t, p) = p[1] * exp.(-t/p[2])
	tdata = centers(htime)
	vdata = htime.weights
	il = length(tdata)
	fit = curve_fit(mexp, tdata[i0:il], vdata[i0:il], pa0)
	coef(fit), stderror(fit), expo.(tdata, coef(fit)...)
end

"""
	plot_waveform(wvfm; window=false, waveform='filtered', wstart=-1, wend=-1, thrp=0, thrn=0)

"""
function plot_waveform(wvfm::DataFrame;
	                   window=false, waveform="filtered", wstart=-1, wend=-1, thrp=0, thrn=0, scatter=true)

	if waveform=="filtered"
		ampl = "fAmpl"
	else
		ampl = "Ampl"
	end

	if window
		p1 = plot(wvfm.Time[wstart:wend], wvfm[wstart:wend, ampl], lw=2, label=waveform, fmt = :png)
		if scatter
			scatter!(wvfm.Time[wstart:wend], wvfm[wstart:wend, ampl],  markersize=2,
					  color = :black, label=false, fmt = :png)
		end
	else
		p1 = plot(wvfm.Time, wvfm[!, ampl], lw=2, label=waveform, fmt = :png)
		if scatter
			scatter!(wvfm.Time, wvfm[!, ampl], markersize=2,
					  color = :black, label=false, fmt = :png)
		end
	end

	hline!([thrp], label=false, fmt = :png)
	hline!([thrn], label=false, fmt = :png)
	xlabel!("Time (μs or index)")
	ylabel!("Amplitude")
	p1
end

"""
	plot_fft(wvfm1, ws; window=false, is=-1, il=-1)

"""
function plot_fft(wvfm1::DataFrame, ws::Float64; window=false, is=-1, il=-1)
	ts, F, freqs = fourier_transform(wvfm1, ws)
	xl = Integer(length(freqs))
	xi = Integer(2 + xl/2)
	fmhz = freqs[xi:xl] * 1e-6      # frequencies in MHz
	ampl = abs.(F)[xi:xl]
	if window
		scatter(fmhz[is:il], ampl[is:il],  markersize=2,
				  color = :black,
	    		  label=false,
				  fmt = :png)
	else
		plot(fmhz, ampl)
	end

	xlabel!("frequency (MHz)")
	ylabel!("Amplitude")
end

"""
	plot_filtered_peaks(wvfmFlt, spks)

"""
function plot_filtered_peaks(wvfmFlt::DataFrame, spks::SPeaks)
	plot(wvfmFlt.fAmpl, lw=2, label=false)
	scatter!(spks.peaks, spks.ys, label="peak")
	scatter!(spks.leftedge, spks.xbase, label="leftedge")
	scatter!(spks.rightedge, spks.xbase, label="rightedge")
	hline!([spks.promsel], label=false)
	xlabel!("peak index")
	ylabel!("peak amplitude")
end

"""
	plot_fit(htime, coeff, tft)

"""
function plot_fit(htime, coeff, tft; savefig=false, fn="")
	tdata = centers(htime)
	vdata = htime.weights
	ps1 = scatter(tdata, vdata, yerr=sqrt.(vdata), markersize=2,
				  color = :black,
	    		  label="data",
				  fmt = :png)
	pp = plot(ps1, tdata, tft, lw=2, label="μ = $(round(coeff[2]*1000, sigdigits=2)) ns", fmt = :png)
	xlabel!("t (μs)")
	ylabel!("frequency")
	if savefig
		png(pp, fn)
	end
	return pp
end

max_stars_repo_name: jjgomezcadenas/laserlab

max_issues_repo_name: jjgomezcadenas/laserlab

max_forks_repo_name: jjgomezcadenas/laserlab