text:
mutable struct NoCrashRewardModel <: AbstractMLRewardModel
    cost_dangerous_brake::Float64 # POSITIVE NUMBER
    reward_in_target_lane::Float64 # POSITIVE NUMBER

    brake_penalty_thresh::Float64 # (POSITIVE NUMBER) if the deceleration is greater than this cost_dangerous_brake will be accured
    target_lane::Int

#####################
    reward_in_first_lane::Float64 # POSITIVE NUMBER
    first_lane::Int
    target_point::Float64
    reward_for_control::Float64
    reward_for_collision::Float64
end

mutable struct EmergencyRewardModel <: AbstractMLRewardModel

    reward_for_control::Float64
    reward_for_collision::Float64
end

mutable struct RampRewardModel <: AbstractMLRewardModel
    cost_dangerous_brake::Float64 # POSITIVE NUMBER
    #reward_in_target_lane::Float64 # POSITIVE NUMBER
    reward_in_first_lane::Float64  # POSITIVE NUMBER
    reward_not_finish::Float64   # NEGATIVE NUMBER
    brake_penalty_thresh::Float64 # (POSITIVE NUMBER) if the deceleration is greater than this cost_dangerous_brake will be accured
    target_lane::Int
    first_lane::Int
end

#XXX temporary
#NoCrashRewardModel() = NoCrashRewardModel(100.,10.,2.5,4)
NoCrashRewardModel() = NoCrashRewardModel(100.,10.,2.5,4,1.,1,2500,10,100)
lambda(rm::NoCrashRewardModel) = rm.cost_dangerous_brake/rm.reward_in_target_lane

RampRewardModel() = RampRewardModel(100.,10.,-10000.,2.5,4,1)

mutable struct NoCrashIDMMOBILModel <: AbstractMLDynamicsModel
    nb_cars::Int
    phys_param::PhysicalParam

    # behaviors::Vector{BehaviorModel}
    # behavior_probabilities::WeightVec
    behaviors::BehaviorGenerator

    adjustment_acceleration::Float64
    lane_change_rate::Float64 # in LANES PER SECOND

    p_appear::Float64 # probability of a new car appearing if the maximum number are not on the road
    appear_clearance::Float64 # minimum clearance for a car to appear

    vel_sigma::Float64 # std of new car speed about v0

    lane_terminate::Bool # if true, terminate the simulation when the car has reached the desired lane
    brake_terminate_thresh::Float64 # terminate simulation if braking is above this (always positive)
    max_dist::Float64 # terminate simulation if distance becomes greater than this
end

#XXX temporary
function NoCrashIDMMOBILModel(nb_cars::Int,
                              pp::PhysicalParam;
                              vel_sigma = 0.5,
                              lane_terminate=false,
                              p_appear=0.5,
                              brake_terminate_thresh=Inf,
                              max_dist=Inf,
                              behaviors=DiscreteBehaviorSet(IDMMOBILBehavior[IDMMOBILBehavior(x[1],x[2],x[3],idx) for (idx,x) in
                                                 enumerate(Iterators.product(["cautious","normal","aggressive"],
                                                        [pp.v_slow+0.5;pp.v_med;pp.v_fast],
                                                        [pp.l_car]))], Weights(ones(9)))
                              )

    return NoCrashIDMMOBILModel(
        nb_cars,
        pp,
        behaviors,
        1., # adjustment accel
        1.0/(2.0*pp.dt), # lane change rate
        p_appear, # p_appear
        35.0, # appear_clearance
        vel_sigma, # vel_sigma
        lane_terminate,
        brake_terminate_thresh,
        max_dist
    )
end

const NoCrashMDP{R<:AbstractMLRewardModel} = MLMDP{MLState, MLAction, NoCrashIDMMOBILModel, R}
const NoCrashPOMDP{R<:AbstractMLRewardModel} = MLPOMDP{MLState, MLAction, MLObs, NoCrashIDMMOBILModel, R}
const SuccessMDP = NoCrashMDP{TargetLaneReward}
const SuccessPOMDP = NoCrashPOMDP{TargetLaneReward}
const RampMDP = NoCrashMDP{RampRewardModel}
const RampPOMDP = NoCrashPOMDP{RampRewardModel}

const NoCrashProblem{R<:AbstractMLRewardModel} = Union{NoCrashMDP{R}, NoCrashPOMDP{R}}
const SuccessProblem = Union{SuccessMDP, SuccessPOMDP}

# TODO issue here VVV need a different way to create observation
create_action(::NoCrashProblem) = MLAction()

# action space = {a in {accelerate,maintain,decelerate}x{left_lane_change,maintain,right_lane_change} | a is safe} U {brake}
immutable NoCrashActionSpace
    NORMAL_ACTIONS::Vector{MLAction} # all the actions except brake
    acceptable::IntSet
    brake::MLAction # this action will be EITHER braking at half the dangerous brake threshold OR the braking necessary to prevent a collision at all time in the future
    ca::MLAction #control_action
end
# TODO for performance, make this a macro?
const NB_NORMAL_ACTIONS = 9

function NoCrashActionSpace(mdp::NoCrashProblem)
    accels = (-mdp.dmodel.adjustment_acceleration, 0.0, mdp.dmodel.adjustment_acceleration)
    lane_changes = (-mdp.dmodel.lane_change_rate, 0.0, mdp.dmodel.lane_change_rate)
    NORMAL_ACTIONS = MLAction[MLAction(a,l) for (a,l) in Iterators.product(accels, lane_changes)] # this should be in the problem
    return NoCrashActionSpace(NORMAL_ACTIONS, IntSet(), MLAction(), MLAction()) # note: brake will be calculated later based on the state
end

function actions(mdp::NoCrashProblem)
    return NoCrashActionSpace(mdp)
end

actions(mdp::NoCrashProblem, s::Union{MLState, MLPhysicalState}) = actions(mdp, s, actions(mdp))

function actions(mdp::NoCrashProblem, s::Union{MLState, MLPhysicalState}, as::NoCrashActionSpace)
    acceptable = IntSet()
    for i in 1:NB_NORMAL_ACTIONS
        a = as.NORMAL_ACTIONS[i]
        ego_y = s.cars[1].y
        # prevent going off the road
        if ego_y == 1. && a.lane_change < 0. || ego_y == mdp.dmodel.phys_param.nb_lanes && a.lane_change > 0.0
            continue
        end
        if a.lane_change == 0.0
            continue
        end
        # prevent running into the person in front or to the side
        #if is_safe(mdp, s, as.NORMAL_ACTIONS[i])
        #    push!(acceptable, i)
        #end

        ###### for Emergency
        push!(acceptable, i)
    end
    #brake_acc = calc_brake_acc(mdp, s)
    #brake = MLAction(brake_acc, 0)
    brake = MLAction(-4, 0.0)
    ca = control_action_emergency(mdp,s,false)

    #ca = MLAction(ca, 0.0)  ####################
    return NoCrashActionSpace(as.NORMAL_ACTIONS, acceptable, brake, ca)
end

calc_brake_acc(mdp::NoCrashProblem{NoCrashRewardModel}, s::Union{MLState, MLPhysicalState}) = min(max_safe_acc(mdp,s), -mdp.rmodel.brake_penalty_thresh/2.0)
calc_brake_acc(mdp::NoCrashProblem{TargetLaneReward}, s::Union{MLState, MLPhysicalState}) = min(max_safe_acc(mdp, s), min(-mdp.dmodel.phys_param.brake_limit/2, -mdp.dmodel.brake_terminate_thresh/2))
calc_brake_acc(mdp::NoCrashProblem{RampRewardModel}, s::Union{MLState, MLPhysicalState}) = min(max_safe_acc(mdp,s), -mdp.rmodel.brake_penalty_thresh/2.0)

iterator(as::NoCrashActionSpace) = as
Base.start(as::NoCrashActionSpace) = 1
function Base.next(as::NoCrashActionSpace, state::Integer)
    while !(state in as.acceptable)
        if state == NB_NORMAL_ACTIONS+1
            return (as.brake, state+1)
        end
        if state > NB_NORMAL_ACTIONS+1
            return (as.ca, state+1)
        end
        state += 1
    end
    return (as.NORMAL_ACTIONS[state], state+1)
end
Base.done(as::NoCrashActionSpace, state::Integer) = state > NB_NORMAL_ACTIONS+2
Base.length(as::NoCrashActionSpace) = length(as.acceptable) + 2

function rand(rng::AbstractRNG, as::NoCrashActionSpace, a::MLAction=MLAction())
    nb_acts = length(as.acceptable)+2
    index = rand(rng,1:nb_acts)
    for (i,a) in enumerate(as) # is this efficient at all ?
        if i == index
            return a
        end
    end
end

"""
Calculate the maximum safe acceleration that will allow the car to avoid a collision if the car in front slams on its brakes
"""
function max_safe_acc(mdp::NoCrashProblem, s::Union{MLState,MLObs}, lane_change::Float64=0.0)
    dt = mdp.dmodel.phys_param.dt
    v_min = mdp.dmodel.phys_param.v_min
    l_car = mdp.dmodel.phys_param.l_car
    bp = mdp.dmodel.phys_param.brake_limit
    ego = s.cars[1]

    car_in_front = 0
    smallest_gap = Inf
    # find car immediately in front
    if length(s.cars) > 1
        for i in 2:length(s.cars)#nb_cars
            if occupation_overlap(s.cars[i].y, 0.0, ego.y, lane_change) # occupying same lane
                gap = s.cars[i].x - ego.x - l_car
                if gap >= -l_car && gap < smallest_gap
                    car_in_front = i
                    smallest_gap = gap
                end
            end
        end
        # calculate necessary acceleration
        if car_in_front == 0
            return Inf
        else
            n_brake_acc = nullable_max_safe_acc(smallest_gap, ego.vel, s.cars[car_in_front].vel, bp, dt)
            return get(n_brake_acc, -mdp.dmodel.phys_param.brake_limit)
        end
    end
    return Inf
end

"""
Return the maximum acceleration that the car behind can have on this step so that it won't hit the car in front if it slams on its brakes
"""
function max_safe_acc(gap, v_behind, v_ahead, braking_limit, dt)
    bp = braking_limit
    v = v_behind
    vo = v_ahead
    g = gap
    # VVV see mathematica notebook
    # return - (bp*dt + 2.*v - sqrt(8.*g*bp + bp^2*dt^2 - 4.*bp*dt*v + 4.*vo^2)) / (2.*dt)
    ret = 0.0
#########################

    t = vo/bp
    ret = (g+vo*t-bp*t*t/2+bp*(t-dt)*(t-dt)/2-v*dt-v*(t-dt)) / (dt*dt/2+(t-dt)*dt)
    return ret
###########################
#=
    try
        ret = - (bp*dt + 2.*v - sqrt(8.*g*bp + bp^2*dt^2 - 4.*bp*dt*v + 4.*vo^2)) / (2.*dt)
    catch ex
        @show gap
        @show v_behind
        @show v_ahead
        @show braking_limit
        @show dt
        rethrow(ex)
    end
    return ret
=#
end

"""
Return max_safe_acc or an empty Nullable if the discriminant is negative.
"""
function nullable_max_safe_acc(gap, v_behind, v_ahead, braking_limit, dt)
    bp = braking_limit
    v = v_behind
    vo = v_ahead
    g = gap
    # VVV see mathematica notebook
    #discriminant = 8.*g*bp + bp^2*dt^2 - 4.*bp*dt*v + 4.*vo^2
    discriminant = 8.*g*bp + bp^2*dt^2 - 2.*bp*dt*v + 4.*vo^2
    if discriminant >= 0.0
        return Nullable{Float64}(- (bp*dt + 2.*v - sqrt(discriminant)) / (2.*dt))
    else
        return Nullable{Float64}()
    end
end

"""
Return get a PID acc for lane change
"""
function get_PID_acc(mdp::NoCrashProblem, s::Union{MLState,MLObs})
    dt = mdp.dmodel.phys_param.dt
    v_min = mdp.dmodel.phys_param.v_min
    l_car = mdp.dmodel.phys_param.l_car
    bp = mdp.dmodel.phys_param.brake_limit
    ego = s.cars[1]

    pp = mdp.dmodel.phys_param
    neighborhood = get_neighborhood(pp, s, 1)

    car_in_front = neighborhood[2]
    if car_in_front == 0
        return 9999
    else
        smallest_gap = s.cars[car_in_front].x - ego.x - l_car
        de_acc = Gipps_acc(smallest_gap, ego.vel, s.cars[car_in_front].vel, bp, dt)
        RSS_acc = RSS_dis_acc(smallest_gap, ego.vel, s.cars[car_in_front].vel, bp, dt)
        if RSS_acc < de_acc
            de_acc = RSS_acc
        end
        #ttc = 4
        #desired_dis = ttc*(ego.vel-s.cars[car_in_front].vel)
        #desired_dis = 5+ego.vel*0.5
        return de_acc
        #=
        if smallest_gap > desired_dis
            return 9999
        else
            return -0.1*(desired_dis-smallest_gap)
        end
        =#
    end
end

function get_PID_lane_change_acc(mdp::NoCrashProblem, s::Union{MLState,MLObs})
    dt = mdp.dmodel.phys_param.dt
    v_min = mdp.dmodel.phys_param.v_min
    l_car = mdp.dmodel.phys_param.l_car
    bp = mdp.dmodel.phys_param.brake_limit
    ego = s.cars[1]

    pp = mdp.dmodel.phys_param
    neighborhood = get_neighborhood(pp, s, 1)

    car_in_front = neighborhood[1]
    if car_in_front == 0
        return 1.5
    else
        smallest_gap = s.cars[car_in_front].x - ego.x - l_car
        de_acc = desired_safe_dis(smallest_gap, ego.vel, s.cars[car_in_front].vel, bp, dt)
        return de_acc
    end
end

function desired_safe_dis(gap, v_behind, v_ahead, braking_limit, dt)
    bp = braking_limit
    v = v_behind
    vo = v_ahead
    g = gap
    # VVV see mathematica notebook
    p = 0.1
    g1 = ((2*v+bp*dt)^2-bp^2*dt+2*v*bp*dt-4*vo^2)/(8*bp)
    return p*(g-g1)
end

function Human_like_acc(gap, v_behind, v_ahead, braking_limit, dt)
    bp = -braking_limit
    br = bp/4
    bf = bp/2
    v = v_behind
    vo = v_ahead
    g = gap
    g1 = v_behind*1.8
    if g>g1
        g_acc = 9999
    else
        g_acc = 0.1*(g-g1)
    end

    if g_acc<-1.5
        g_acc = -1.5
    end
    return g_acc
end

function AEB_acc(gap, v_behind, v_ahead, braking_limit, dt)
    bp = braking_limit
    v = v_behind
    vo = v_ahead
    g = gap
    # VVV see mathematica notebook
    p = 0.1
    #g1 = ((2*v+bp*dt)^2-bp^2*dt+2*v*bp*dt-4*vo^2)/(8*bp)
    #g_RSS = v*dt+1/2*1*dt^2+(v+dt*1)^2/2/bp-vo^2/2/bp
    g_AEB = v_behind*0.9
    if g<g_AEB
        return -4
    else
        return 9999
    end
end

function Gipps_acc(gap, v_behind, v_ahead, braking_limit, dt)
    bp = -braking_limit
    br = bp/4
    bf = bp/2
    v = v_behind
    vo = v_ahead
    g = gap
    # VVV see mathematica notebook
    #p = 0.1
    #g1 = ((2*v+bp*dt)^2-bp^2*dt+2*v*bp*dt-4*vo^2)/(8*bp)
    m = br^2*dt^2-br*(2*gap-v*dt-vo^2/bf)
    if m > 0
        v_d = br*dt + sqrt(m)
        g_acc = (v_d-v)/dt
    else
        g_acc = -1.5
    end
    if g_acc<-1.5
        g_acc = -1.5
    end
    return g_acc
end

function RSS_dis_acc(gap, v_behind, v_ahead, braking_limit, dt)
    bp = braking_limit
    v = v_behind
    vo = v_ahead
    g = gap
    # VVV see mathematica notebook
    p = 0.1
    #g1 = ((2*v+bp*dt)^2-bp^2*dt+2*v*bp*dt-4*vo^2)/(8*bp)
    g_RSS = v*dt+1/2*1*dt^2+(v+dt*1)^2/2/bp-vo^2/2/bp
    if g<g_RSS
        return -4
    else
        return 9999
    end
end

function Detect_RSS_Braking(mdp::NoCrashProblem, s::Union{MLState,MLObs})
    dt = mdp.dmodel.phys_param.dt
    v_min = mdp.dmodel.phys_param.v_min
    l_car = mdp.dmodel.phys_param.l_car
    bp = mdp.dmodel.phys_param.brake_limit
    ego = s.cars[1]

    pp = mdp.dmodel.phys_param
    neighborhood = get_neighborhood(pp, s, 1)

    car_in_front = neighborhood[2]
    if car_in_front > 0
        smallest_gap = s.cars[car_in_front].x - ego.x - l_car
        RSS_acc = RSS_dis_acc(smallest_gap, ego.vel, s.cars[car_in_front].vel, bp, dt)
        if RSS_acc <-2.3
            return true
        end
    end
    return false
end

"""
Test whether, if the ego vehicle takes action a, it will always be able to slow down fast enough if the car in front slams on his brakes and won't pull in front of another car so close they can't stop
"""
function is_safe(mdp::NoCrashProblem, s::Union{MLState,MLObs}, a::MLAction)
    dt = mdp.dmodel.phys_param.dt
    if a.acc >= max_safe_acc(mdp, s, a.lane_change)
        return false
    end
    ############## check whether the lane change action is reasonable!
    if a.lane_change > 0.3
        return false
    end
    vel = s.cars[1].vel
    if vel>34 && a.acc>0
        return false
    end
    if s.x < 1500 && a.lane_change!=0.0 && isinteger(s.cars[1].y)
        return false
    end

    if s.x > 2500 && a.lane_change!=0.0 && isinteger(s.cars[1].y)
        return false
    end

    # check whether we will go into anyone else's lane so close that they might hit us or we might run into them
    if isinteger(s.cars[1].y) && a.lane_change != 0.0
        l_car = mdp.dmodel.phys_param.l_car
        for i in 2:length(s.cars)
            car = s.cars[i]
            ego = s.cars[1]
            if car.x < ego.x + l_car && occupation_overlap(ego.y, a.lane_change, car.y, 0.0)  # ego is in front of car
                # New definition of safe - the car behind can brake at max braking to avoid the ego if the ego
                # slams on his brakes
                # XXX IS THIS RIGHT??
                # I think I need a better definition of "safe" here
                gap = ego.x - car.x - l_car
                if gap <= 0.0
                    return false
                end
                n_braking_acc = nullable_max_safe_acc(gap, car.vel, ego.vel, mdp.dmodel.phys_param.brake_limit, dt)

                # if isnull(n_braking_acc) || get(n_braking_acc) < -mdp.dmodel.phys_param.brake_limit
                if isnull(n_braking_acc) || get(n_braking_acc) < max_accel(mdp.dmodel.behaviors)
                    return false
                end
            end
        end
    end
    return true
end

#XXX temp
#create_state(p::NoCrashProblem) = MLState(0.0, 0.0, Array(CarState, p.dmodel.nb_cars), nothing)
create_state(p::NoCrashProblem) = MLState(0.0, 0.0, Array{CarState}(p.dmodel.nb_cars), true, nothing)
create_observation(pomdp::Union{NoCrashPOMDP,SuccessPOMDP}) = MLObs(0.0, 0.0, Array(CarStateObs, pomdp.dmodel.nb_cars), nothing)

function generate_s(mdp::NoCrashProblem, s::MLState, a::MLAction, rng::AbstractRNG)

    @if_debug dbg_rng = copy(rng)

    sp::MLState=create_state(mdp)
    pp = mdp.dmodel.phys_param
    dt = pp.dt
    nb_cars = length(s.cars)
    resize!(sp.cars, nb_cars)
    #sp.x = 0
    sp.terminal = s.terminal

    ## Calculate deltas ##
    #====================#

    dxs = Array{Float64}(nb_cars)
    dvs = Array{Float64}(nb_cars)
    dys = Array{Float64}(nb_cars)
    lcs = Array{Float64}(nb_cars)

    # agent
    dvs[1] = a.acc*dt
    dxs[1] = s.cars[1].vel*dt + a.acc*dt^2/2.
    lcs[1] = a.lane_change
    dys[1] = a.lane_change*dt

    for i in 2:nb_cars
        neighborhood = get_neighborhood(pp, s, i)

        behavior = s.cars[i].behavior

        acc = gen_accel(behavior, mdp.dmodel, s, neighborhood, i, rng)
        dvs[i] = dt*acc
        dxs[i] = (s.cars[i].vel + dvs[i]/2.)*dt

        lcs[i] = gen_lane_change(behavior, mdp.dmodel, s, neighborhood, i, rng)
        dys[i] = lcs[i] * dt
    end

    ## Consistency checking ##
    #========================#

    # first prevent lane changes into each other
    changers = IntSet()
    for i in 1:nb_cars
        if lcs[i] != 0
            push!(changers, i)
        end
    end
    sorted_changers = sort!(collect(changers), by=i->s.cars[i].x, rev=true) # this might be slow because anonymous functions are slow
    # from front to back

    if length(sorted_changers) >= 2 #something to compare
        # iterate through pairs
        iter_state = start(sorted_changers)
        j, iter_state = next(sorted_changers, iter_state)
        while !done(sorted_changers, iter_state)
            i = j
            j, iter_state = next(sorted_changers, iter_state)
            car_i = s.cars[i] # front
            car_j = s.cars[j] # back

            # check if they are both starting to change lanes on this step
            if isinteger(car_i.y) && isinteger(car_j.y)

                # check if they are near each other lanewise
                if abs(car_i.y - car_j.y) <= 2.0

                    # make sure there is a conflict longitudinally
                    # if car_i.x - car_j.x <= pp.l_car || car_i.x + dxs[i] - car_j.x + dxs[j] <= pp.l_car
                    # if car_i.x - car_j.x <= mdp.dmodel.appear_clearance # made more conservative on 8/19
                    # if car_i.x - car_j.x <= get_idm_s_star(car_j.behavior.p_idm, car_j.vel, car_j.vel-car_i.vel) # upgraded to sstar on 8/19
                    ivp = car_i.vel + dt*dvs[i]
                    jvp = car_j.vel + dt*dvs[j]
                    ixp = car_i.x + dt*(car_i.vel + ivp)/2.0
                    jxp = car_j.x + dt*(car_j.vel + jvp)/2.0
                    n_max_acc_p = nullable_max_safe_acc(ixp-jxp-pp.l_car, jvp, ivp, pp.brake_limit,dt)
                    if ixp - jxp <= pp.l_car*0.2 || car_i.x - car_j.x <= pp.l_car*0.2 # || isnull(n_max_acc_p) || get(n_max_acc_p) < -pp.brake_limit

                        # check if they are moving towards each other
                        # if dys[i]*dys[j] < 0.0 && abs(car_i.y+dys[i] - car_j.y+dys[j]) < 2.0
                        if true # prevent lockstepping 8/19 (doesn't prevent it that well)

                            # make j stay in his lane
                            dys[j] = 0.0
                            lcs[j] = 0.0
                        end
                    end
                end
            end
        end
    end

    # second, prevent cars hitting each other due to noise
    sorted = sort!(collect(1:length(s.cars)), by=i->s.cars[i].x, rev=true)

    if length(sorted) >= 2 #something to compare
        # iterate through pairs
        iter_state = start(sorted)
        j, iter_state = next(sorted, iter_state)
        while !done(sorted, iter_state)
            i = j
            j, iter_state = next(sorted, iter_state)
            if j == 1
                continue # don't check for the ego since the ego does not have noise
            end
            car_i = s.cars[i]
            car_j = s.cars[j]

            # check if they overlap longitudinally
            if car_j.x + dxs[j] > car_i.x + dxs[i] - pp.l_car

                # check if they will be in the same lane
                if occupation_overlap(car_i.y + dys[i], car_j.y + dys[j])
                    # warn and nudge behind
                    #=
                    @if_debug begin
                        println("Conflict because of noise: front:$i, back:$j")
                        Gallium.@enter generate_s(mdp, s, a, dbg_rng)
                    end
                    if i == 1
                        # warn("Car nudged because noise would cause a crash (ego in front).")
                        #error("Car nudged because noise would cause a crash (ego in front).")
                    else
                        # warn("Car nudged because noise would cause a crash.")
                        #error("Car nudged because noise would cause a crash.")
                    end
                    =#
                    dxs[j] = car_i.x + dxs[i] - car_j.x - 1.01*pp.l_car
                    dvs[j] = 2.0*(dxs[j]/dt - car_j.vel)
                end
            end
        end
    end

    ## Dynamics and Exits ##
    #======================#

    exits = IntSet()
    for i in 1:nb_cars
        car = s.cars[i]
        xp = car.x + (dxs[i] - dxs[1])
        yp = car.y + dys[i]
        # velp = max(min(car.vel + dvs[i],pp.v_max), pp.v_min)
        velp = max(car.vel + dvs[i], 0.0) # removed speed limits on 8/13
        # note lane change is updated above

        if dvs[i]/dt < -mdp.dmodel.brake_terminate_thresh
            sp.terminal = Nullable{Symbol}(:brake)
        end

        # check if a lane was crossed and snap back to it
        if isinteger(car.y)
            # prevent a multi-lane change in a single timestep
            if abs(yp-car.y) > 1.
                yp = car.y + sign(dys[i])
            end
        else # car.y is not an integer
            if floor(yp) >= ceil(car.y)
                yp = ceil(car.y)
            end
            if ceil(yp) <= floor(car.y)
                yp = floor(car.y)
            end
        end

        # if yp < 1.0 || yp > pp.nb_lanes
        #     @show i
        #     @show yp
        #     println("mdp = $mdp")
        #     println("s = $s")
        #     println("a = $a")
        # end
        #print("yp")
        #println(yp)
        @assert yp >= 1.0 && yp <= pp.nb_lanes

        if xp < 0.0 || xp >= pp.lane_length
            push!(exits, i)
        else
            sp.cars[i] = CarState(xp, yp, velp, lcs[i], car.behavior, s.cars[i].id)
        end
    end

    next_id = maximum([c.id for c in s.cars]) + 1

    deleteat!(sp.cars, exits)
    nb_cars -= length(exits)

    ## Generate new cars ##
    #=====================#

    if nb_cars < mdp.dmodel.nb_cars && rand(rng) <= mdp.dmodel.p_appear

        behavior = rand(rng, mdp.dmodel.behaviors)
        vel = typical_velocity(behavior) + randn(rng)*mdp.dmodel.vel_sigma

        clearances = Array{Float64}(pp.nb_lanes)
        fill!(clearances, Inf)
        closest_cars = Array{Int}(pp.nb_lanes)
        fill!(closest_cars, 0)
        sstar_margins = Array{Float64}(pp.nb_lanes)
        if vel > s.cars[1].vel
            # put at back
            # sstar is the sstar of the new guy
            for i in 1:length(s.cars)
                lowlane, highlane = occupation_lanes(s.cars[i].y, lcs[i])
                back = s.cars[i].x - pp.l_car
                if back < clearances[lowlane]
                    clearances[lowlane] = back
                    closest_cars[lowlane] = i
                end
                if back < clearances[highlane]
                    clearances[highlane] = back
                    closest_cars[highlane] = i
                end
            end
            for j in 1:pp.nb_lanes
                other = closest_cars[j]
                if other == 0
                    sstar = 0
                else
                    sstar = get_idm_s_star(behavior.p_idm, vel, vel-s.cars[other].vel)
                end
                sstar_margins[j] = clearances[j] - sstar
            end
        else
            for i in 1:length(s.cars)
                lowlane, highlane = occupation_lanes(s.cars[i].y, lcs[i])
                front = pp.lane_length - (s.cars[i].x + pp.l_car) # l_car is half the length of the old car plus half the length of the new one
                if front < clearances[lowlane]
                    clearances[lowlane] = front
                    closest_cars[lowlane] = i
                end
                if front < clearances[highlane]
                    clearances[highlane] = front
                    closest_cars[highlane] = i
                end
            end
            for j in 1:pp.nb_lanes
                other = closest_cars[j]
                if other == 0
                    sstar = 0
                else
                    sstar = get_idm_s_star(s.cars[other].behavior.p_idm,
                                           s.cars[other].vel,
                                           s.cars[other].vel-vel)
                end
                sstar_margins[j] = clearances[j] - sstar
            end
        end

        margin, lane = findmax(sstar_margins)

        if margin > 0.0
            if vel > s.cars[1].vel
                # at back
                push!(sp.cars, CarState(0.0, lane, vel, 0.0, behavior, next_id))
            else
                push!(sp.cars, CarState(pp.lane_length, lane, vel, 0.0, behavior, next_id))
            end
        end
    end

    if mdp.dmodel.lane_terminate && sp.cars[1].y == mdp.rmodel.target_lane
        sp.terminal = Nullable{Any}(:lane)
    elseif sp.x > mdp.dmodel.max_dist
        sp.terminal = Nullable{Any}(:distance)
    end

    # sp.crashed = is_crash(mdp, s, sp, warning=false)
    # sp.crashed = false

    @assert sp.cars[1].x == s.cars[1].x # ego should not move
    sp.x = s.x + sp.cars[1].vel*0.75
    #ca = control_action_adapt(mdp,s)
    ca = control_action_emergency(mdp,s,false)
#=
    if s.Control_Signal && (a == ca)
       sp.Control_Signal = true
    else
       sp.Control_Signal = false
    end
=#

    if a == ca && isinteger(sp.cars[1].y)
        sp.Control_Signal = true
    else
        sp.Control_Signal = false
    end

    if Detect_RSS_Braking(mdp,s)
        sp.t = 1
    else
        sp.t = 0
    end

#=
    if s.x > 2500 && isinteger(s.cars[1].y) && s.t == 0.0
        sp.Control_Signal = true
        sp.t = 1.0
    end
=#

    #if s.cars[1].y == 4.0 && s.t == 0.0
    #    sp.Control_Signal = true
    #    sp.t = 1.0
    #end
    return sp
end

function control_action(mdp::NoCrashProblem, s::Union{MLState, MLPhysicalState})
     vel = s.cars[1].vel
     acc = 0.0
     if vel < 26
         acc = 0.5
     end
     if vel > 27
         acc = -0.5
     end

     if acc >= max_safe_acc(mdp, s, 0.0)
         acc = max_safe_acc(mdp, s, 0.0)-0.1
     end
     lc = 0.0
     if is_safe(mdp,s,MLAction(acc,-mdp.dmodel.lane_change_rate)) && (s.cars[1].y != 1.0) && (s.x>500) && (s.x<2500)
         lc = -mdp.dmodel.lane_change_rate
     end
#=
     if is_safe(mdp,s,MLAction(acc,mdp.dmodel.lane_change_rate)) && (s.cars[1].y > 3.2) && (s.x>1200) && (s.x<2500)
         lc = -mdp.dmodel.lane_change_rate
     end
     if is_safe(mdp,s,MLAction(acc,mdp.dmodel.lane_change_rate)) && (s.cars[1].y > 2.2) && (s.x>2000) && (s.x<2500)
         lc = -mdp.dmodel.lane_change_rate
     end
     if is_safe(mdp,s,MLAction(acc,mdp.dmodel.lane_change_rate)) && (s.cars[1].y > 1.2) && (s.x>2250) && (s.x<2500)
         lc = -mdp.dmodel.lane_change_rate
     end
=#
     return MLAction(acc,lc)
end

function control_action_adapt2(mdp::NoCrashProblem, s::Union{MLState, MLPhysicalState})
     vel = s.cars[1].vel
     acc = 0.0
     if vel < 27
         acc = 0.5
     end
     if vel > 28
         acc = -0.5
     end
     if acc >= max_safe_acc(mdp, s, 0.0)
         acc = max_safe_acc(mdp, s, 0.0)-0.1
     end
     lc = 0.0
     if (s.cars[1].y != 1.0) && (s.x>1500) && (s.x<2500)
         if is_safe(mdp,s,MLAction(acc,-mdp.dmodel.lane_change_rate))
            lc = -mdp.dmodel.lane_change_rate
        elseif vel>22.3 && vel<31
           acc = -0.5
         end
     end
     return MLAction(acc,lc)
end

function control_action_adapt(mdp::NoCrashProblem, s::Union{MLState, MLPhysicalState})
     vel = s.cars[1].vel
     acc = 0.0
     if vel < 30
         acc = 0.5
     end
     if vel > 31
         acc = -0.5
     end
     if acc >= max_safe_acc(mdp, s, 0.0)
         acc = max_safe_acc(mdp, s, 0.0)-0.1
     end
     lc = 0.0
     pp = mdp.dmodel.phys_param
     neighborhood = get_neighborhood(pp, s, 1)
     PID_acc = get_PID_lane_change_acc(mdp,s)
     if (s.cars[1].y != 1.0) && (s.x>1500) && (s.x<2500)

         if is_safe(mdp,s,MLAction(acc,-mdp.dmodel.lane_change_rate))
            lc = -mdp.dmodel.lane_change_rate
        elseif  PID_acc < max_safe_acc(mdp, s, -mdp.dmodel.lane_change_rate) && vel<33
            acc = PID_acc
        #elseif  max_safe_acc(mdp, s, -mdp.dmodel.lane_change_rate) < 0.5 && max_safe_acc(mdp, s, 0.0)>0.5 && vel<33
    #        acc = -0.5
    #    elseif (neighborhood[1] !=0) && (neighborhood[4] != 0)
        elseif vel>22.3
            acc = -2
         end
     end
     return MLAction(acc,lc)
end

function control_action_emergency(mdp::NoCrashProblem, s::Union{MLState, MLPhysicalState}, RSS_signal::Bool)
     vel = s.cars[1].vel
     acc = 0.0
     des_vel = 19.5
     acc = (des_vel-vel)
     if acc>1.5
         acc = 1.5
     end
     if acc<-1.5
         acc = -1.5
     end
##########################
     dt = mdp.dmodel.phys_param.dt
     v_min = mdp.dmodel.phys_param.v_min
     l_car = mdp.dmodel.phys_param.l_car
     bp = mdp.dmodel.phys_param.brake_limit
     ego = s.cars[1]

     pp = mdp.dmodel.phys_param
     neighborhood = get_neighborhood(pp, s, 1)

     car_in_front = neighborhood[2]

     if car_in_front > 0
         smallest_gap = s.cars[car_in_front].x - ego.x - l_car
         #de_acc = Gipps_acc(smallest_gap, ego.vel, s.cars[car_in_front].vel, bp, dt)
         de_acc = Human_like_acc(smallest_gap, ego.vel, s.cars[car_in_front].vel, bp, dt)
         #RSS_acc = RSS_dis_acc(smallest_gap, ego.vel, s.cars[car_in_front].vel, bp, dt)
         g_AEB = AEB_acc(smallest_gap, ego.vel, s.cars[car_in_front].vel, bp, dt)
#=
         if RSS_signal && RSS_acc < de_acc
             de_acc = RSS_acc
         end
=#
         if RSS_signal && g_AEB < de_acc
             de_acc = g_AEB
         end

         if de_acc < acc
             acc = de_acc
         end
     end

     if acc<-4
         acc = -4
     end

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

     lc = 0.0
#=
     if is_safe(mdp,s,MLAction(acc,-mdp.dmodel.lane_change_rate)) && (s.cars[1].y >2.2)
         lc = -mdp.dmodel.lane_change_rate
     end

     if is_safe(mdp,s,MLAction(acc,mdp.dmodel.lane_change_rate)) && (s.cars[1].y <1.8)
         lc = mdp.dmodel.lane_change_rate
     end
=#
     return MLAction(acc,lc)
end

function control_action_emergency(mdp::NoCrashProblem, s::Union{MLState, MLPhysicalState})
    return control_action_emergency(mdp,s,false)
end

function reward(mdp::NoCrashProblem{NoCrashRewardModel}, s::MLState, a::MLAction, sp::MLState)
    r = 0.0
#=
    if sp.cars[1].y == mdp.rmodel.target_lane
        r += mdp.rmodel.reward_in_target_lane
    end

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

    if sp.cars[1].vel < 20
        r -= 10*(20-sp.cars[1].vel)
    end

    if s.Control_Signal
        r += mdp.rmodel.reward_for_control
    end

    if (s.x > mdp.rmodel.target_point) && (sp.cars[1].y != mdp.rmodel.target_lane)
           r -= mdp.rmodel.reward_in_first_lane*10
    end

    if (s.x > mdp.rmodel.target_point) && (sp.cars[1].y == mdp.rmodel.target_lane)
           r += mdp.rmodel.reward_in_first_lane
    end
#################################################

=#

###############################
    if s.Control_Signal
        r += 0
    else
        r -= mdp.rmodel.reward_for_control
    end

    if judge_collision(mdp,s)
        r -=mdp.rmodel.reward_for_collision
    end

    if a.acc < -mdp.rmodel.brake_penalty_thresh
        r -= mdp.rmodel.cost_dangerous_brake
    end
###############################

    return r
end

function reward(mdp::NoCrashProblem{RampRewardModel}, s::MLState, ::MLAction, sp::MLState)
    r = 0.0
    if sp.cars[1].y == mdp.rmodel.first_lane
        r += mdp.rmodel.reward_in_target_lane
    end
    nb_brakes = detect_braking(mdp, s, sp)
    r -= mdp.rmodel.cost_dangerous_brake*nb_brakes
    return r
end

"""
Return the number of braking actions that occured during this state transition
"""
function detect_braking(mdp::NoCrashProblem, s::MLState, sp::MLState, threshold::Float64)
    nb_brakes = 0
    nb_leaving = 0
    dt = mdp.dmodel.phys_param.dt
    for (i,c) in enumerate(s.cars)
        if length(sp.cars) >= i-nb_leaving
            cp = sp.cars[i-nb_leaving]
        else
            break
        end
        if cp.id != c.id
            nb_leaving += 1
            continue
        else
            acc = (cp.vel-c.vel)/dt
            if acc < -threshold
                nb_brakes += 1
            end
        end
    end
    # @assert nb_leaving <= 5 # sanity check - can remove this if it is violated as long as it doesn't happen all the time
    return nb_brakes
end

detect_braking(mdp::NoCrashProblem, s::MLState, sp::MLState) = detect_braking(mdp, s, sp, mdp.rmodel.brake_penalty_thresh)
detect_braking(mdp::SuccessProblem, s::MLState, sp::MLState) = detect_braking(mdp, s, sp, mdp.dmodel.brake_terminate_thresh)

"""
Return the ids of cars that brake during this state transition
"""
function braking_ids(mdp::NoCrashProblem, s::MLState, sp::MLState, threshold=mdp.rmodel.brake_penalty_thresh)
    braking = Int[]
    nb_leaving = 0
    dt = mdp.dmodel.phys_param.dt
    for (i,c) in enumerate(s.cars)
        if length(sp.cars) >= i-nb_leaving
            cp = sp.cars[i-nb_leaving]
        else
            break
        end
        if cp.id != c.id
            nb_leaving += 1
            continue
        else
            acc = (cp.vel-c.vel)/dt
            if acc < -threshold
                push!(braking, c.id)
            end
        end
    end
    # @assert nb_leaving <= 5 # sanity check - can remove this if it is violated as long as it doesn't happen all the time
    return braking
end

"""
Return the maximum braking (a positive number) by any car during the transition between s an sp
"""
function max_braking(mdp::NoCrashProblem, s::MLState, sp::MLState)
    nb_leaving = 0
    dt = mdp.dmodel.phys_param.dt
    min_acc = 0.0
    for (i,c) in enumerate(s.cars)
        if length(sp.cars) >= i-nb_leaving
            cp = sp.cars[i-nb_leaving]
        else
            break
        end
        if cp.id != c.id
            nb_leaving += 1
            continue
        else
            acc = (cp.vel-c.vel)/dt
            if acc < min_acc
                min_acc = acc
            end
        end
    end
    # @assert nb_leaving <= 5 # sanity check - can remove this if it is violated as long as it doesn't happen all the time
    return -min_acc
end

"""
Assign behaviors to a given physical state.
"""
function initial_state(mdp::NoCrashProblem, ps::MLPhysicalState, rng::AbstractRNG)
    s = MLState(ps, Array(CarState, length(ps.cars)))
    s.cars[1] = CarState(ps.cars[1], NORMAL)
    for i in 2:length(s.cars)
        behavior = rand(rng, mdp.dmodel.behaviors)
        s.cars[i] = CarState(ps.cars[i], behavior)
    end
    return s
end

function initial_state(mdp::NoCrashProblem, rng::AbstractRNG=Base.GLOBAL_RNG)
    return relaxed_initial_state(mdp, 200, rng)
end

function generate_o(mdp::NoCrashProblem, s::MLState, a::MLAction, sp::MLState)
    return MLObs(sp)
end

function generate_sor(pomdp::Union{NoCrashPOMDP, SuccessPOMDP}, s::MLState, a::MLAction, rng::AbstractRNG)
    sp, r = generate_sr(pomdp, s, a, rng)
    o = generate_o(pomdp, s, a, sp)
    return sp, o, r
end

discount(mdp::Union{MLMDP,MLPOMDP}) = mdp.discount
isterminal(mdp::Union{MLMDP,MLPOMDP},s::MLState) = !isnull(s.terminal)

max_stars_repo_name: zhcao/HighwayExit

max_issues_repo_name: zhcao/HighwayExit

max_forks_repo_name: zhcao/HighwayExit


text:
# of various kinds, and also makes plots that allows
# those plots to be examined in some useful way.
# The thing to plot is timing (possibly with variations)
# for each sha thingy, also registred with respect to
# calender date when the test was run.  Also, log what kind
# of machinery was used to run the tests (osx,linux,
# CPU stats (if known) etc.).
#
# Plan: Start out using straight (open) code for a couple of
#       benchmaks, then modify to use a macro that
#       instruments code (look at the @time macro for
#       inspiration), then make something that will be able
#       to run the tests.
#
# tic() toc() @timed @time @timev @timeit
# https://www.juliabloggers.com/timing-in-julia/# Time everything, make a file that logs benchmarks

ext: jl


text:
using Gadfly, BBOBFunctions

##

function plot_fun(f)
    r = 5

    c, n, nline = Geom.contour(levels=exp10.(range(-6, stop=10, length=40))), 800, 1000
    if any(f .== [BBOBFunctions.F16,BBOBFunctions.F17,BBOBFunctions.F18,BBOBFunctions.F19])
        c, n, nline = Geom.contour(levels=exp10.(range(-6, stop=8, length=20))), 500, 1000
    end

    p1 = plot(
        layer(x=[f.x_opt[1]],y=[f.x_opt[2]],Geom.point,Theme(default_color=colorant"red")),
        layer(
            z=(x,y)->f([x,y])-f.f_opt, x=range(-r, stop=r, length=n), y=range(-r, stop=r, length=n),
            c,
        ),
        Coord.cartesian(xmin=-r,xmax=r,ymin=-r,ymax=r),
        Guide.title(string(f))
    )

    # line plot
    x = collect(range(-r, stop=r, length=nline))

    p2 = plot(
        layer(x=x,y=[ f([xi,f.x_opt[2],f.x_opt[3]]) for xi in x]-f.f_opt, Geom.line),
        layer(x=x,y=[ f([f.x_opt[1],xi,f.x_opt[3]]) for xi in x]-f.f_opt, Geom.line,Theme(default_color=colorant"orange")),
        layer(x=x,y=[ f([f.x_opt[1],f.x_opt[2],xi]) for xi in x]-f.f_opt, Geom.line,Theme(default_color=colorant"red")),
#        layer(x=[f.x_opt[1]],y=[f([f.x_opt[1],0,0])]-f.f_opt,Geom.point),
#        layer(x=[f.x_opt[2]],y=[f([0,f.x_opt[2],0])]-f.f_opt,Geom.point,Theme(default_color=colorant"orange")),
#        layer(x=x,y=[f([xi,xi,xi]) for xi in x]-f.f_opt,Geom.line,Theme(default_color=colorant"red")),
        Coord.cartesian(xmin=-r,xmax=r),
        Guide.title(string(f)),
        Scale.y_log10
    )
    p1,p2
end

outdir = joinpath(Pkg.dir(),"BBOBFunctions","data","plots","functions")

for i in 20:length(enumerate(BBOBFunctions.BBOBFunction))

    f = getfield(BBOBFunctions,Symbol("F$i"))
    info(f)
    p1,p2 = plot_fun(f)

    draw(PDF(joinpath(outdir,"F$i.pdf"),16cm,16cm),p1)
    draw(SVG(joinpath(outdir,"F$i.svg"),16cm,16cm),p1)
    draw(PDF(joinpath(outdir,"1D","F$i.pdf"),16cm,12cm),p2)
    draw(SVG(joinpath(outdir,"1D","F$i.svg"),16cm,12cm),p2)
end

##


text:
JULIA_ENV = realpath(joinpath(@__DIR__, "..", ".."))
data_root = realpath(joinpath(@__DIR__, "..", "..", "data"))
data_raw_input_root = joinpath(data_root, "input", "raw")
data_input_root = joinpath(data_root, "input", "dami")
data_parameter_file = joinpath(data_input_root, "parameters.jser")
data_output_root = joinpath(data_root, "output")

worker_list = [("localhost", 1)]
exeflags = `--project="$JULIA_ENV"`
sshflags= `-i path/to/ssh/key/file`
SAMPLE_SIZE_LIMIT = 10_000

fmt_string = "[{name} | {date} | {level}]: {msg}"
loglevel = "debug"

### data directories ###
raw_data_dirs = Dict("semantic" => ["Arrhythmia", "Annthyroid", "Cardiotocography", "HeartDisease", "Hepatitis", "InternetAds", "PageBlocks", "Parkinson", "Pima", "SpamBase", "Stamps", "Wilt"],
                     "literature" => ["ALOI", "Glass", "Ionosphere", "KDDCup99", "Lymphography", "PenDigits", "Shuttle", "WBC", "WDBC", "WPBC", "Waveform"])
data_dirs = vcat(values(raw_data_dirs)...)

max_stars_repo_name: englhardt/ocs-evaluation

max_issues_repo_name: englhardt/ocs-evaluation

max_forks_repo_name: englhardt/ocs-evaluation


text:
export component, non_simple_extension

################################################################################
#
#  Constructor
#
################################################################################

@doc Markdown.doc"""
    NumberField(f::Vector{PolyElem{<:NumFieldElem}}, s::String="_\$", check = true)
                                              -> NumField, Vector{NumFieldElem}

Given a list $f_1, \ldots, f_n$ of univariate polynomials in $K[x]$ over
some number field $K$, constructs the extension $K[x_1, \ldots, x_n]/(f_1(x_1),
\ldots, f_n(x_n))$.

# Examples

```jldoctest
julia> Qx, x = QQ["x"];

julia> K, a = NumberField([x^2 - 2, x^2 - 3], "a")
(Non-simple number field with defining polynomials fmpq_mpoly[x1^2 - 2, x2^2 - 3], NfAbsNSElem[a1, a2])
```
"""
function _doc_stub_nf2 end

# To work around a bug in the built documentation.
#
abstract type DocuDummy2 end

@doc (@doc _doc_stub_nf2)
NumberField(::DocuDummy2)

@doc (@doc _doc_stub_nf2)
NumberField(::Vector{<:PolyElem{<:Union{NumFieldElem, fmpq}}}, ::String, check::Bool = true)

################################################################################
#
#  Basis
#
################################################################################

@doc Markdown.doc"""
    basis(L::NonSimpleNumField) -> Vector{NumFieldElem}

Returns the canonical basis of a non-simple extension $L/K$. If
$L = K(a_1,\dotsc,a_n)$ where each $a_i$ has degree $d_i$, then the basis will
be $a_1^{i_1}\dotsm a_d^{i_d}$ with $0 \leq i_j \leq d_j - 1$ for $1 \leq j \leq n$.

# Examples

```jldoctest
julia> Qx, x = QQ["x"];

julia> K, (a1, a2) = NumberField([x^2 - 2, x^2 - 3], "a");

julia> basis(K)
4-element Vector{NfAbsNSElem}:
 1
 a1
 a2
 a1*a2
```
"""
basis(::NonSimpleNumField)

################################################################################
#
#  Defining polynomials
#
################################################################################

export defining_polynomials

@doc Markdown.doc"""
    defining_polynomials(L::NonSimpleNumField) -> Vector{PolyElem}

Given a non-simple number field $L/K$, constructed as $L =
K[x]/(f_1,\dotsc,f_r)$, return the vector containing the $f_i$'s.
"""
defining_polynomials(::NonSimpleNumField)

defining_polynomials(K::NfRelNS) = K.abs_pol

defining_polynomials(K::NfAbsNS) = K.abs_pol

################################################################################
#
#  Conversion to simple extension
#
################################################################################

@doc Markdown.doc"""
    simple_extension(L::NonSimpleNumField) -> SimpleNumField, Map

Given a non-simple extension $L/K$, this function computes a simple extension
$M/K$ and a $K$-linear isomorphism $M \to L$.
"""
simple_extension(::NonSimpleNumField)

################################################################################
#
#  Check consistency
#
################################################################################

# This is used for the check in the constructor
function _check_consistency(K::NonSimpleNumField)
  QQz, z = PolynomialRing(base_field(K), "z")
  for i = 1:length(K.pol)
    v = [zero(QQz) for i in 1:length(K.pol)]
    v[i] = z
    p = evaluate(K.pol[i], v)
    if !is_irreducible(p)
      return false
    end
  end
  lg = gens(K)
  el = lg[1]
  f = minpoly(el)
  d = degree(f)
  for i = 2:length(lg)
    el += lg[i]
    f = minpoly(el)
    while !is_squarefree(f)
      el += lg[i]
      f = minpoly(el)
    end
    if degree(f) != prod(degree(K.pol[j], j) for j = 1:i)
      return false
    end
    if !is_irreducible(f)
      return false
    end
  end
  return true
end

################################################################################
#
#  Component
#
################################################################################

@doc Markdown.doc"""
    component(L::NonSimpleNumField, i::Int) -> SimpleNumField, Map

Given a non-simple extension $L/K$, this function returns the simple number field
corresponding to the $i$-th component of $L$ together with its embedding.
"""
function component(K::NonSimpleNumField, i::Int)
  fl = is_univariate(K.pol[i])
  @assert fl
  kx, _ = PolynomialRing(base_field(K), "x", cached = false)
  g = to_univariate(kx, K.pol[i])
  gK = gens(K)
  Ki, a = number_field(g, cached = false, check = false)
  mp = hom(Ki, K, gK[i])
  return Ki, mp
end

################################################################################
#
#  Non-simplify
#
################################################################################

function non_simple_extension(K::SimpleNumField)
  @assert base_field(K) isa FlintRationalField
  @assert is_normal(K)
  G, mG = automorphism_group(K)
  _subs = _subgroups_for_non_simple_extension(G)
  if length(_subs) == 0
    return [defining_polynomial(K)]
  end
  subf = Dict()
  for (subgrps, indice) in _subs
    for H in subgrps
      if !haskey(subf, H)
        subf[H] = defining_polynomial(fixed_field(K, [mG(H[2](h)) for h in H[1]])[1])
      end
    end
  end

  Qx = Globals.Qx

  l = Inf

  res = nothing

  for (subgrps, indice) in _subs
    v = [ subf[H] for H in subgrps ]
    if length(string(v)) < l
      res = v
    end
    L, = number_field(v)
    @assert is_isomorphic(simple_extension(L)[1], K)
  end

  return res
end

function _subgroups_for_non_simple_extension(G, maxorder = div(order(G), 2), maxnum = 5)
  sub = subgroups(G)
  n = order(G)
  res = []
  curmin = 0
  for i in 2:maxnum
    facs = _factorizations(n, i)
    for f in facs
      if curmin < sum(f) && !isempty(res)
        continue
      end

      newmin = sum(map(x -> divexact(n, x), f))

      possible_comb = [ [h for h in sub if order(h[1]) == f[i]] for i in 1:length(f) ]

      for it in Iterators.product(possible_comb...)
        maps = map(x -> x[2], collect(it))
        k = foldl((x, y) -> intersect(x, y)[2], maps)
        if order(domain(k)) != 1
          continue
        end

        if isempty(res)
          push!(res, (collect(it), collect(order.(first.(it)))))
          curmin = newmin
        else
          if newmin < curmin
            empty!(res)
            push!(res, (collect(it), collect(order.(first.(it)))))
            curmin = newmin
          elseif newmin == curmin
            push!(res, (collect(it), collect(order.(first.(it)))))
          end
        end
      end
    end
  end
  if length(res) == 0
    return res
  end
  maxlength = maximum(Int[length(x[2]) for x in res])
  res = [x for x in res if length(x[2]) == maxlength]
  return res
end

function _factorizations(n, parts)
  if parts == 1
    return Vector{Int}[Int[n]]
  end

  res = Vector{Vector{Int}}()
  for k in 2:isqrt(n)
    if mod(n, k) == 0
      nk = divexact(n, k)
      _f = _factorizations(nk, parts - 1)::Vector{Vector{Int}}
      for i in 1:length(_f)
        push!(_f[i], k)
      end
      append!(res, _f)
    end
  end
  res = map(sort!, res)
  return unique!(res)
end

################################################################################
#
#  Simplified simple extension
#
################################################################################

@doc Markdown.doc"""
    simplified_simple_extension(L::NonSimpleNumField) -> SimpleNumField, Map

Given a non-simple extension $L/K$, this function returns an isomorphic simple number field
with a "small" defining equation together with the isomorphism.
"""
function simplified_simple_extension(L::NonSimpleNumField; cached::Bool = true, is_abelian::Bool = false)
  OL = maximal_order(L)
  B = lll_basis(OL)
  B1 = _sieve_primitive_elements(B)
  a = B1[1]
  I = t2(a)
  for i = 2:min(50, length(B1))
    J = t2(B1[i])
    if J < I
      a = B1[i]
      I = J
    end
  end
  f = minpoly(a)
  Ls, gLs = number_field(f, cached = cached, check = false)
  mp = hom(Ls, L, a)
  return Ls, mp
end

function simplified_simple_extension(K::NfAbsNS; cached::Bool = true, is_abelian::Bool = false)
  OK = maximal_order(K)
  if is_abelian
    OS = _lll_CM(OK)
    OK.lllO = OS
  else
    OS = lll(OK)
  end
  el = _simplify(OS)
  f = minpoly(el)
  L, mL = number_field(f, cached = cached, check = false)
  mp = hom(L, K, el, check = false)
  return L, mp
end


text:
using PkgName, Base.Test

function equal(A::AbstractMatrix{T},V::AbstractMatrix{T},H::AbstractMatrix{T},m::Int) where{T}
    x = norm(A*V[:,1:m]-V*H)
    y = 10*eps(T)
    if x < y
        return true
    else
        return false
    end
end

A = sprand(100,100,0.3)
H, V = hessenberg(A,30)
@test equal(A,V,H,30) == true


text:
import ArgParse, TOML

function parse_commandline()
    s = ArgParse.ArgParseSettings()

    ArgParse.@add_arg_table! s begin
        "--parse", "-p"
            help = "parse a dotfile configuration"
            arg_type = AbstractString
    end

    return ArgParse.parse_args(s)
end

function main()
    parsed_args = parse_commandline()
    println("Parsed args:")
    for (arg,val) in parsed_args
        println("$arg => $val")
    end

    if !isnothing(parsed_args["parse"])
        dotfile_config = TOML.parsefile(parsed_args["parse"])
        println(dotfile_config)
        TOML.print(dotfile_config)
    else
        println("no config in sight!")
    end
end


text:
# compute kron(A, B)[Irow, Icol]
function kron_sub!(out, A, B, Irow, Icol)
  @assert size(out, 1) == length(Irow)
  @assert size(out, 2) == length(Icol)

  m, n = size(A)
  p, q = size(B)
  for col=1:length(Icol)
    j = Icol[col]
    ac = div(j-1, q) + 1
    bc = mod(j-1, q) + 1

    for row=1:length(Irow)
        i = Irow[row]
        ar = div(i-1, p) + 1
        br = mod(i-1, p) + 1

        out[row, col] = A[ar, ac] * B[br, bc]
    end
  end
  out
end
kron_sub(A, B, Irow, Icol) = kron_sub!(Matrix{Float64}(undef, length(Irow), length(Icol)), A, B, Irow, Icol)

# compute kron(A, B)[Irow, Icol]
function skron_sub!(out, A, B, Irow, Icol)
  @assert size(out, 1) == length(Irow)
  @assert size(out, 2) == length(Icol)

  m, n = size(A)
  p, q = size(B)
  for col=1:length(Icol)
    j = Icol[col]
    ac1 = div(j-1, q) + 1
    bc1 = mod(j-1, q) + 1
    ac2 = div(j-1, n) + 1
    bc2 = mod(j-1, n) + 1

    for row=1:length(Irow)
        i = Irow[row]
        ar1 = div(i-1, p) + 1
        br1 = mod(i-1, p) + 1
        ar2 = div(i-1, m) + 1
        br2 = mod(i-1, m) + 1

        out[row, col] = ( A[ar1, ac1] * B[br1, bc1] + B[ar2, ac2] * A[br2, bc2] ) / 2.
    end
  end
  out
end
skron_sub(A, B, Irow, Icol) = skron_sub!(Matrix{Float64}(undef, length(Irow), length(Icol)), A, B, Irow, Icol)

function _getElemKron(A::AbstractMatrix, B::AbstractMatrix, row, col)
    m, n = size(A)
    p, q = size(B)

    ac = div(col-1, p) + 1
    bc = mod(col-1, p) + 1
    ar = div(row-1, p) + 1
    br = mod(row-1, p) + 1

    @inbounds A[ar, ac] * B[br, bc]
end

function _getElemKron(S::AbstractMatrix, xi::AbstractVector, row, col)
    p = size(S, 1)

    ac = div(col-1, p) + 1
    bc = mod(col-1, p) + 1
    ar = div(row-1, p) + 1
    br = mod(row-1, p) + 1

    @inbounds S[ar, ac] * xi[br] * xi[bc]
end

function _getElemKron(xi::AbstractVector, S::AbstractMatrix, row, col)
    p = size(S, 1)

    ac = div(col-1, p) + 1
    bc = mod(col-1, p) + 1
    ar = div(row-1, p) + 1
    br = mod(row-1, p) + 1

    @inbounds S[br, bc] * xi[ar] * xi[ac]
end

function _getElemSKron(S::AbstractMatrix, xi::AbstractVector, row, col)
    p = size(S, 1)

    ac = div(col-1, p) + 1
    bc = mod(col-1, p) + 1
    ar = div(row-1, p) + 1
    br = mod(row-1, p) + 1

    @inbounds ( S[ar, ac] * xi[br] * xi[bc] + S[br, bc] * xi[ar] * xi[ac] ) / 2.
end

_getElemSKron(A::AbstractMatrix, B::AbstractMatrix, row, col) =
    ( _getElemKron(A, B, row, col) + _getElemKron(B, A, row, col) ) / 2.

# computes variance of Var(Sx - Sy)[indS, indS]
# where indS is a list of elements in [1:p, 1:p]
function _varSxSy(Sx, nx, Sy, ny, indS)

    varS = zeros(length(indS), length(indS))   # var(Sx) + var(Sy)

    I = CartesianIndices(Sx)
    for ia = 1:length(indS)
      a = indS[ia]
      for ib = 1:length(indS)
        b = indS[ib]

        i, j = Tuple( I[a] )
        k, l = Tuple( I[b] )

        varS[ia, ib] = (Sx[i, k] * Sx[j, l] + Sx[i, l] * Sx[j, k]) / (nx - 1) + (Sy[i, k] * Sy[j, l] + Sy[i, l] * Sy[j, k]) / (ny - 1)
      end
    end

    varS
end

function getSupport(x::SymmetricSparseIterate, τ::Float64=1e-3)
  p = size(x, 1)

  S = falses(p, p)
  for ci=1:p
      for ri=ci:p
          if abs(x[ri, ci] > τ)
              S[ri, ci] = true
              S[ci, ri] = true
          end
      end
  end
  S
end

function getSupport(x::SparseIterate, p::Int, τ::Float64=1e-3)
  S = falses(p, p)
  ind = 0
  for ci=1:p
      for ri=ci:p
          ind += 1
          if abs(x[ind] > τ)
              S[ri, ci] = true
              S[ci, ri] = true
          end
      end
  end
  S
end

getReducedSupport(x::SparseIterate, τ::Float64=1e-3) = abs.(x) .> τ

function getLinearSupport(row::Int, col::Int, S1::BitArray{2}, S2::BitArray{2})
  p = size(S1, 1)

  supp = Vector{Int64}()
  for i=1:p
    for j=1:p
      if S1[i, j] | S2[i,j]
        push!(supp, (j-1)*p + i)
      end
    end
  end
  makeFirst!(supp, p, row, col)
end

function getLinearSupport(row::Int, col::Int, S::BitArray{2})
  p = size(S, 1)

  supp = Vector{Int64}()
  for i=1:p
    for j=1:p
      if S[i, j]
        push!(supp, (j-1)*p + i)
      end
    end
  end
  makeFirst!(supp, p, row, col)
end

function makeFirst!(supp::Vector{Int64}, p, row, col)
  ind = (col - 1) * p + row
  pos = findfirst(isequal(ind), supp)
  if  pos === nothing
      pushfirst!(supp, ind)
  else
      supp[pos], supp[1] = supp[1], supp[pos]
  end

  if row != col
    ind = (row - 1) * p + col
    pos = findfirst(isequal(ind), supp)
    if  pos === nothing
        pushfirst!(supp, ind)
    else
        supp[pos], supp[2] = supp[2], supp[pos]
    end
  end
  supp
end

function makeFirst!(supp::Vector{Int64}, indElem)
  pos = findfirst(isequal(indElem), supp)
  if  pos === nothing
      pushfirst!(supp, indElem)
  else
      supp[pos], supp[1] = supp[1], supp[pos]
  end
  supp
end

function computeSimulationResult(res::Vector{DiffPrecResultNormal}, trueParam::Float64, α::Float64=0.05)
  Eestim = 0.
  bias = 0.
  coverage = 0
  lenCoverage = 0.

  num_rep = length(res)
  zα = quantile(Normal(), 1. - α / 2.)

  for rep=1:num_rep
    if res[rep].p - zα * res[rep].std < trueParam && res[rep].p + zα * res[rep].std > trueParam
      coverage += 1
    end
    lenCoverage += 2. * zα * res[rep].std
    Eestim += res[rep].p
  end

  bias = Eestim / num_rep - trueParam

  ep = [res[rep].p for rep=1:num_rep]
  mean_p = mean(ep)
  var_p = var(ep)

  SimulationResult(bias, coverage / num_rep, lenCoverage / num_rep, mean_p, var_p)
end

function computeSimulationResult(res::Vector{DiffPrecResultBoot}, trueParam::Float64, α::Float64=0.05)
  Eestim = 0.
  bias = 0.
  coverage = 0
  lenCoverage = 0.

  num_rep = length(res)
  z1 = α / 2.
  z2 = 1. - α / 2.

  for rep=1:num_rep
    lb, ub = quantile!(res[rep].boot_p, [z1, z2])
    if lb < trueParam && ub > trueParam
      coverage += 1
    end
    lenCoverage += (ub - lb)
    Eestim += res[rep].p
  end

  bias = Eestim / num_rep - trueParam

  ep = [res[rep].p for rep=1:num_rep]
  mean_p = mean(ep)
  var_p = var(ep)

  SimulationResult(bias, coverage / num_rep, lenCoverage / num_rep, mean_p, var_p)
end


text:
#===============================================================================

    DeconvolutionMatrices.jl

    types for K matrices and their inverses
    used in deconvolution.jl

    Author: Tom Price
    Date:   June 2019

===============================================================================#

import StaticArrays: SVector, SMatrix

#===============================================================================
    struct >>
===============================================================================#

mutable struct DeconvolutionMatrices
    k_s             ::Vector{SMatrix{C,C,Float_T} where {C}}
    k_inv_vec       ::Vector{SMatrix{C,C,Float_T} where {C}}
    inv_note        ::String
end

#===============================================================================
    null constructor >>
===============================================================================#

function DeconvolutionMatrices(
    calibration_data    ::CalibrationData{<: NumberOfChannels, <: Union{Int_T,Float_T}},
    calibration_args    ::CalibrationParameters
)
    s = size(calibration_data.array)
    w = s[1]
    c = s[2]
    v = calibration_args.k_method == well_proc_vec ? w : 1
    empty_matrix = SMatrix{c,c,Float_T}(fill(NaN_T,c,c))
    DeconvolutionMatrices(
        SVector{v}(fill(empty_matrix,v)),
        SVector{w}(fill(empty_matrix,w)),
        "")
end

#===============================================================================
    constant >>
===============================================================================#

## this constant was previously loaded using
## JLD.load("$LOAD_FROM_DIR/defines/k4dcv_ip84_calib79n80n81_vec.jld")["k4dcv"]
const K4DCV =
    DeconvolutionMatrices(
        Vector{SMatrix{2,2,Float_T}}([
            [0.645552  0.152178 ;  0.354448  0.847822],
            [0.673361  0.209691 ;  0.326639  0.790309],
            [0.66636   0.0408344;  0.33364   0.959166],
            [0.627177  0.144452 ;  0.372823  0.855548],
            [0.696393  0.05625  ;  0.303607  0.94375 ],
            [0.686193  0.0750364;  0.313807  0.924964],
            [0.661831  0.015839 ;  0.338169  0.984161],
            [0.678446  0.145038 ;  0.321554  0.854962],
            [0.664971  0.079949 ;  0.335029  0.920051],
            [0.702946 -0.0320236;  0.297054  1.03202 ],
            [0.676903  0.096115 ;  0.323097  0.903885],
            [0.721044  0.15577  ;  0.278956  0.84423 ],
            [0.702459  0.0888078;  0.297541  0.911192],
            [0.728335  0.167087 ;  0.271665  0.832913],
            [0.649742  0.0844896;  0.350258  0.91551 ],
            [0.664093  0.0787242;  0.335907  0.921276]]),
        Vector{SMatrix{2,2,Float_T}}([
            [1.71842  -0.308443 ; -0.718416  1.30844 ],
            [1.70447  -0.452244 ; -0.704466  1.45224 ],
            [1.53338  -0.0652802; -0.533375  1.06528 ],
            [1.77233  -0.299243 ; -0.772331  1.29924 ],
            [1.47428  -0.0878709; -0.474279  1.08787 ],
            [1.51346  -0.122778 ; -0.513464  1.12278 ],
            [1.52349  -0.0245189; -0.523488  1.02452 ],
            [1.60283  -0.271908 ; -0.602828  1.27191 ],
            [1.57268  -0.13666  ; -0.572677  1.13666 ],
            [1.40417   0.0435714; -0.404172  0.956429],
            [1.55631  -0.165491 ; -0.556308  1.16549 ],
            [1.49349  -0.275565 ; -0.493488  1.27556 ],
            [1.48487  -0.14472  ; -0.484871  1.14472 ],
            [1.48404  -0.297706 ; -0.484038  1.29771 ],
            [1.61965  -0.149472 ; -0.619648  1.14947 ],
            [1.57384  -0.134486 ; -0.573837  1.13449 ]]),
        "")

max_stars_repo_name: MakerButt/chaipcr

max_issues_repo_name: MakerButt/chaipcr

max_forks_repo_name: MakerButt/chaipcr


text:
__precompile__()

module JQuantmod

export get_price, get_irx, get_sp500, get_log_return, get_log_return_sp500, get_log_irx

using RCall
using Printf

function init()
    R"""
library(quantmod)

get_irx_ratio <- function(from ,to){
    IRX = getSymbols("^IRX", from = from, to=to,
            auto.assign=FALSE)
    IRX_Adj_Close = data.frame(IRX)$IRX.Adjust
    IRX_Annual_Return = (IRX_Adj_Close/100)[2:length(IRX_Adj_Close)]
    IRX_Annual_Log_Return = log(1+IRX_Annual_Return)
    IRX_Daily_Log_Return = IRX_Annual_Log_Return/252
    IRX_Daily_Log_Return
}

get_price <- function(s, to, from){
    INSTR = getSymbols(paste("^",s,sep=""), from = from, to=to,
            auto.assign=FALSE)
    cmd = sprintf("data.frame(INSTR)$%s.Adjusted", s)
    INSTR_Adj_Close=eval(parse(text=cmd))

    INSTR_Adj_Close
}

get_price2 <- function(s, to, from){
    INSTR = getSymbols(s, from = from, to=to,
            auto.assign=FALSE)
    cmd = sprintf("data.frame(INSTR)$%s.Adjusted", s)
    INSTR_Adj_Close=eval(parse(text=cmd))

    INSTR_Adj_Close
}
"""
end

"""
    get_log_return(s::String, from::String, to::String)

    Return the daily log price for symbol `s` from `from` to `to`
"""
function get_log_return(s::String, from::String, to::String)
    init()
    p = get_price(s, from, to)
    return log.(p[2:end])-log.(p[1:end-1])
end

"""
    get_price(s::String, from::String, to::String)

    Return the daily price for symbol `s` from `from` to `to`
"""
function get_price(s::String, from::String, to::String)
    init()
    if s[1]=='^'
        s = (s[2:end])
        reval(@sprintf("r = get_price('%s',from='%s',to='%s')", s, from, to))
    else
        reval(@sprintf("r = get_price2('%s',from='%s',to='%s')", s, from, to))
    end
    @rget r
    r
end

"""
    get_log_irx(from::String, to::String)

    Return the daily risk-free log return from `from` to `to`
"""
function get_log_irx(from::String, to::String)
    init()
    reval(@sprintf("r = get_irx_ratio('%s','%s')", from, to))
    @rget r
    r
end

"""
    get_irx(from::String, to::String)

    Return the daily risk-free return from `from` to `to`
"""
function get_irx(from::String, to::String)
    init()
    r = get_log_irx(from, to)
    return exp.(r) .- 1
end

"""
    get_sp500(from::String, to::String)

    Return the daily return for S&P500 from `from` to `to`
"""
function get_sp500(from::String, to::String)
    init()
    return get_price("^GSPC", from, to)
end

"""
    get_log_return_sp500(from::String, to::String)

    Return the daily log return for S&P500 from `from` to `to`
"""
function get_log_return_sp500(from::String, to::String)
    init()
    return get_log_return("^GSPC", from, to)
end
end


text:
#!/usr/bin/env julia

unshift!(LOAD_PATH, joinpath(@__DIR__, "..", "ext"))

include("loadmeta.jl")
include("generate.jl")


text:
@forwardplotfun plot_autocorr
@forwardplotfun plot_bpv
@forwardplotfun plot_compare
@forwardplotfun plot_density
@forwardplotfun plot_dist
@forwardplotfun plot_dist_comparison
@forwardplotfun plot_elpd
@forwardplotfun plot_energy
@forwardplotfun plot_ess
@forwardplotfun plot_forest
@forwardplotfun plot_hdi
@forwardplotfun plot_kde
@forwardplotfun plot_khat
@forwardplotfun plot_loo_pit
@forwardplotfun plot_mcse
@forwardplotfun plot_pair
@forwardplotfun plot_parallel
@forwardplotfun plot_posterior
@forwardplotfun plot_ppc
@forwardplotfun plot_rank
@forwardplotfun plot_separation
@forwardplotfun plot_trace
@forwardplotfun plot_violin

@deprecate plot_joint(args...; kwargs...) plot_pair(args...; kwargs...)

function convert_arguments(::typeof(plot_compare), df, args...; kwargs...)
    pdf = topandas(Val(:DataFrame), df; index_name=:name)
    return tuple(pdf, args...), kwargs
end

function convert_arguments(::typeof(plot_elpd), data, args...; kwargs...)
    dict = Dict(k => try
        topandas(Val(:ELPDData), v)
    catch
        convert_to_inference_data(v)
    end for (k, v) in pairs(data))
    return tuple(dict, args...), kwargs
end
function convert_arguments(::typeof(plot_khat), df, args...; kwargs...)
    return tuple(topandas(Val(:ELPDData), df), args...), kwargs
end

for f in (
    :plot_autocorr,
    :plot_ess,
    :plot_mcse,
    :plot_pair,
    :plot_parallel,
    :plot_posterior,
    :plot_trace,
    :plot_violin,
)
    @eval begin
        function convert_arguments(::typeof($(f)), data, args...; kwargs...)
            idata = convert_to_inference_data(data; group=:posterior)
            return tuple(idata, args...), kwargs
        end
    end
end

for f in (:plot_autocorr, :plot_ess, :plot_mcse, :plot_posterior, :plot_violin)
    @eval begin
        function convert_arguments(::typeof($(f)), data::AbstractArray, args...; kwargs...)
            return tuple(data, args...), kwargs
        end
    end
end

function convert_arguments(::typeof(plot_energy), data, args...; kwargs...)
    dataset = convert_to_dataset(data; group=:sample_stats)
    return tuple(dataset, args...), kwargs
end

for f in (:plot_density, :plot_forest, :plot_rank)
    @eval begin
        function convert_arguments(
            ::typeof($(f)), data, args...; transform=identity, group=:posterior, kwargs...
        )
            tdata = transform(data)
            dataset = convert_to_dataset(tdata; group=group)
            return tuple(dataset, args...), kwargs
        end
    end
end

for f in (:plot_density, :plot_forest)
    @eval begin
        function convert_arguments(
            ::typeof($(f)),
            data::Union{AbstractVector,Tuple},
            transform=identity,
            group=:posterior,
            args...;
            kwargs...,
        )
            tdata = transform(data)
            datasets = map(tdata) do datum
                return convert_to_dataset(datum; group=group)
            end
            return tuple(datasets, args...), kwargs
        end
        function convert_arguments(
            ::typeof($(f)),
            data::AbstractVector{<:Real},
            transform=identity,
            group=:posterior,
            args...;
            kwargs...,
        )
            tdata = transform(data)
            dataset = convert_to_dataset(tdata; group=group)
            return tuple(dataset, args...), kwargs
        end
    end
end


text:
"""
    load_function(path) --> compiled_f

load a `StaticCompiledFunction` from a given path. This object is callable.
"""
function load_function(path; filename="obj")
    instantiate(deserialize(joinpath(path, "$filename.cjl")))
end

struct LazyStaticCompiledFunction{rt, tt}
    f::Symbol
    path::String
    name::String
    filename::String
    reloc::IdDict{Any,String}
end

"""
    unsafe_pointer_from_objref(x)

Sometimes Julia embeds immutables like `Base.string` into code, and julia
will error if you call `pointer_from_objref(string)`, claiming that it
doesn't have a pointer even though that's a lie.
"""
unsafe_pointer_from_objref(x) = ccall(:jl_value_ptr, Ptr{Cvoid}, (Any,), x)

function instantiate(p::LazyStaticCompiledFunction{rt, tt}) where {rt, tt}
    # LLVM.load_library_permantly(dirname(Libdl.dlpath(Libdl.dlopen("libjulia"))))
    lljit = LLVM.LLJIT(;tm=LLVM.JITTargetMachine())
    jd = LLVM.JITDylib(lljit)
    flags = LLVM.API.LLVMJITSymbolFlags(LLVM.API.LLVMJITSymbolGenericFlagsExported, 0)
    ofile = LLVM.MemoryBufferFile(joinpath(p.path, "$(p.filename).o")) #$(Libdl.dlext)

    # Set all the uninitialized global variables to point to julia values from the relocation table
    for (val, name) ∈ p.reloc
        address = LLVM.API.LLVMOrcJITTargetAddress(reinterpret(UInt, unsafe_pointer_from_objref(val)))

        symbol = LLVM.API.LLVMJITEvaluatedSymbol(address, flags)
        gv = LLVM.API.LLVMJITCSymbolMapPair(LLVM.mangle(lljit, name), symbol)
        mu = absolute_symbols(Ref(gv))
        LLVM.define(jd, mu)
    end
    # consider switching to one mu for all gvs instead of one per gv.
    # I tried that already, but I got an error saying
    # JIT session error: Symbols not found: [ __Type_Vector_Float64___274 ]

    # Link to libjulia
    prefix = LLVM.get_prefix(lljit)
    dg = LLVM.CreateDynamicLibrarySearchGeneratorForProcess(prefix)
    LLVM.add!(jd, dg)

    LLVM.add!(lljit, jd, ofile)
    fptr = pointer(LLVM.lookup(lljit, "julia_" * p.name))

    StaticCompiledFunction{rt, tt}(p.f, fptr, lljit, p.reloc)
end

function absolute_symbols(symbols)
    ref = LLVM.API.LLVMOrcAbsoluteSymbols(symbols, length(symbols))
    LLVM.MaterializationUnit(ref)
end

struct StaticCompiledFunction{rt, tt}
    f::Symbol
    ptr::Ptr{Nothing}
    jit::LLVM.LLJIT
    reloc::IdDict{Any, String}
end

function Base.show(io::IO, f::StaticCompiledFunction{rt, tt}) where {rt, tt}
    types = [tt.parameters...]
    print(io, String(f.f), "(", join(("::$T" for T ∈ tt.parameters), ',')  ,") :: $rt")
end

function (f::StaticCompiledFunction{rt, tt})(args...) where {rt, tt}
    Tuple{typeof.(args)...} == tt || error("Input types don't match compiled target $((tt.parameters...,)). Got arguments of type $(typeof.(args))")
    out = RefValue{rt}()
    refargs = Ref(args)
    ccall(f.ptr, Nothing, (Ptr{rt}, Ref{tt}), pointer_from_objref(out), refargs)
    out[]
end

instantiate(f::StaticCompiledFunction) = f


text:
module DiscreteTimeSEIRD
using DifferentialEquations

function seird(du, u, parms, t)
    S,E,Ir,Id,R,D = u
    β,δ,γ,Γ,μ,ϵ,ω = parms
    I = Ir + Id
    N = S + E + I + R + D
    dS = S - β*S*I/N + ω*R
    dE = E + β*S*I/N - γ*E + ϵ
    dIr = Ir + γ*(1-μ)E - γ*Ir + ϵ
    dId = Id + γ*μ*E - Γ*Id + ϵ
    dR = R + γ*Ir - ω*R
    dD = D + Γ*Id
    du = [dS, dE, dIr, dId, dR, dD]
end

end


text:
"""
A macro `@replace` that replaces one function with another in an expression
"""
module Replace

using Cassette

export @replace

function _replace!(ctx::T, mapping, EvalModule) where T
  for (key, value) in mapping
    @eval EvalModule function Cassette.overdub(ctx::$T, fn::typeof($key), args...)
      return $value(args...)
    end
  end
  return ctx
end

"""
`@replace foo bar` replaces occurences in expression `bar` of the first item of
each iteration of `foo` with the second item in each iteration of `foo`.

```julia
1 == @replace Dict(sin=>cos) sin(0.0)
exp(0.0) + 1 == @replace Dict(sin=>cos, tan=>exp) sin(0.0) + tan(0.0)
```
"""
macro replace(mapping, ex, EvalModule::Module=__module__)
  @eval EvalModule using Cassette
  @eval EvalModule Cassette.@context Ctx
  quote
    ctx = $EvalModule.Ctx()
    _replace!(ctx, $(esc(mapping)), $EvalModule)
    Cassette.@overdub(ctx, $(esc(ex)))
  end
end

"""
`@replace foo bar qux` replaces the function `foo` with the function `bar` where it occurs in expression `qux`

```julia
1 == @replace sin cos sin(0.0)
```
"""
macro replace(this, forthat, inhere, EvalModule::Module=__module__)
  quote
    @replace Dict($(esc(this))=>$(esc(forthat))) $(esc(inhere)) $EvalModule
  end
end

end


text:
using MarketData

c = convert(TimestampArray, cl)

facts("Timestamp construction") do

    context("fields are constructed") do
        @fact Timestamp(Date(2000,1,1),12).timestamp => Date(2000,1,1)
        @fact Timestamp(Date(2000,1,1),12).value     => 12
    end
end

facts("Getindex works on TimeType") do

    context("getindex on single Date, DateTime") do
        @fact length(c[(Date(2000,1,1))])    => 0
        @fact c[(Date(2000,1,3))][1].value   => 111.94
    end

    context("getindex on Vector Date, DateTime") do
        @fact length(c[[(Date(2000,1,1)), Date(2000,1,3)]])  => 1
        @fact c[[(Date(2000,1,1)), Date(2000,1,3)]][1].value => 111.94
    end

    context("getindex on StepRange Date, DateTime") do
        @fact length(c[[(Date(2000,1,1)):Date(2000,1,3)]])  => 1
        @fact c[[(Date(2000,1,1)):Date(2000,1,3)]][1].value => 111.94
    end
end


text:
using VehicleCalculate,GLMakie
include("zlsb250.jl")
zs650zhudong = PrimTransmission(44,72,"曲轴")
zs650zhuzhou = DrivingTransmission(15,16,20,21,23,27,"主轴")
zs650fuzhou = DrivenTransmission(39,33,34,31,31,34,"副轴")
zs650flywheel = Flywheel(14,41,"大小飞")

zs650chuan = Transmission(zs650zhudong,zs650zhuzhou,zs650fuzhou,zs650flywheel)

zs650 = vehiclespeed("ZS650",zs650chuan)
zs650.second_ratio
zlsb250.second_ratio

fg = Figure(resolution = (600, 400))
ax1 = Axis(fg[1,2])
ax2 = Axis(fg[1,1])
x_zs650 = 1:6
y_zs650 = zs650.second_ratio
lines!(ax1,x_zs650,y_zs650,label = "zs650")
lines!(ax1,x_sl250,y_sl250,label = "sl250")

scatter!(ax1,x_zs650,y_zs650)
scatter!(ax1,x_sl250,y_sl250)
xb = [1,2,3,4]
yb = [4,4,4,4]
hist!(ax2,yb,direction=:x,bins = 4,offset = 0.1)
fg

h = ["$i" * "档" for i in ["一","二","三","四","五","六"]]
v = ["zs650","SL250"]
c = hcat(zs650.second_ratio,zlsb250.second_ratio)

tableplot(h,v,c,reso = (600,400))

max_stars_repo_name: powerzf/VehicleCalculate

max_issues_repo_name: powerzf/VehicleCalculate

max_forks_repo_name: powerzf/VehicleCalculate


text:
export cplx2real, pb_cplx2real, real2cplx

varname_cplx2real(varname::String) = (varname*"_Re",varname*"_Im")

"""
	pb = pb_cplx2real(pb_C::Problem)

Convert a complex Polynomial Optimization Problem into a real polynomial
problem in real variables.
"""
function pb_cplx2real(pb_C::Problem)
	pb = Problem()
	for (varName, varType) in get_variables(pb_C)
		if varType <: Complex
			varName_real, varName_imag = varname_cplx2real(varName)
			add_variable!(pb, Variable(varName_real, Real))
			add_variable!(pb, Variable(varName_imag, Real))
		else
			add_variable!(pb, Variable(varName, varType))
		end
	end

	realPart, imagPart = cplx2real(pb_C.objective)
	set_objective!(pb, realPart)

	for (cstrName, cstr) in get_constraints(pb_C)
		realPart, imagPart = cplx2real(cstr.p)
		cstrName_real, cstrName_imag = varname_cplx2real(cstrName)

		if length(realPart) != 0
			# Set bounds to proper infty, easier to detect... TODO : Missing attribute ctr_kind ?
			lb = real(cstr.lb)==-Inf ? -Inf-im*Inf : real(cstr.lb)
			ub = real(cstr.ub)== Inf ? +Inf+im*Inf : real(cstr.ub)
			cstrreal = lb << realPart << ub
			cstr.precond != :none && (cstrreal.precond = cstr.precond)
			add_constraint!(pb, cstrName_real, cstrreal)
		end
		if length(imagPart) != 0
			lb = imag(cstr.lb)==-Inf ? -Inf-im*Inf : imag(cstr.lb)
			ub = imag(cstr.ub)== Inf ? +Inf+im*Inf : imag(cstr.ub)
			cstrimag = lb << imagPart << ub
			cstr.precond != :none && (cstrimag.precond = cstr.precond)
			add_constraint!(pb, cstrName_imag, cstrimag)
		end
	end
	return pb
end

# Conversion of all complex variables to real ones in the Poly structures

"""
	realPart, imagPart = cplx2real(expo::Exponent)

	Convert a complex Exponent in complex variables into `realPart` and
	`imagPart` polynomials of twice as many variables, real and imag parts of
	`expo` variables. Done recursively with the `cplx2real_rec` function.
"""
function cplx2real(expo::Exponent, λ::Number)
	a=real(λ)
	b=imag(λ)
	if expo.degree.explvar  == 1 && expo.degree.conjvar  == 1
		if length(expo.expo)==1
			realPart = Polynomial(); add!(realPart, 0)
			imagPart = Polynomial(); add!(imagPart, 0)

			z1 = first(expo.expo)
			# println("Square : ", z1[1])
			x1 = Variable(z1[1].name*"_Re", Real)
			y1 = Variable(z1[1].name*"_Im", Real)
			x1x1 = SortedDict{Variable, Degree}(x1=>Degree(2,0))
			y1y1 = SortedDict{Variable, Degree}(y1=>Degree(2,0))
			# (a+ib) * conj(z1) * z1 = ax1^2+ay1^2+ibx1^2+iby1^2
			if a != 0
				realPart = Polynomial(SortedDict{Exponent, Number}(Exponent(x1x1)=>a, Exponent(y1y1)=>a))
			end
			if b != 0
				imagPart = Polynomial(SortedDict{Exponent, Number}(Exponent(x1x1)=>b, Exponent(y1y1)=>b))
			end
			return realPart, imagPart
		else
			if first(expo.expo)[2].conjvar==1
				z1 = first(expo.expo)
				z2 = last(expo.expo)
			else
				z2 = first(expo.expo)
				z1 = last(expo.expo)
			end
			# println("Bilinear : ", z1, ", ", z2)
			x1 = Variable(z1[1].name*"_Re", Real)
			y1 = Variable(z1[1].name*"_Im", Real)
			x2 = Variable(z2[1].name*"_Re", Real)
			y2 = Variable(z2[1].name*"_Im", Real)
			# (a+ib)*conj(z1) * z2  = (a+ib)*(x1-iy1)*(x2+iy2)
			# = (a+ib)*(x1.x2 + y1.y2 + ix1.y2-ix2.y1)
			# = a.x1.x2 + a.y1.y2 + ia.x1.y2-ia.x2.y1+ib.x1.x2 + ib.y1.y2 - b.x1.y2+b.x2.y1
			# =
			#    a.x1.x2+a.y1.y2-b.x1.y2+b.x2.y1
			#    +
			#    i(a.x1.y2-a.x2.y1+b.x1.x2+b.y1.y2)
			x1x2 = SortedDict{Variable, Degree}(x1=>Degree(1,0), x2=>Degree(1,0))
			y1y2 = SortedDict{Variable, Degree}(y1=>Degree(1,0), y2=>Degree(1,0))
			x1y2 = SortedDict{Variable, Degree}(x1=>Degree(1,0), y2=>Degree(1,0))
			x2y1 = SortedDict{Variable, Degree}(x2=>Degree(1,0), y1=>Degree(1,0))

			realPart = Polynomial(SortedDict{Exponent, Number}(Exponent(x1x2)=>a, Exponent(y1y2)=>a, Exponent(x1y2)=>-b, Exponent(x2y1)=>b))
			imagPart = Polynomial(SortedDict{Exponent, Number}(Exponent(x1x2)=>b, Exponent(y1y2)=>b, Exponent(x1y2)=>a, Exponent(x2y1)=>-a))
			return realPart, imagPart
		end
	elseif expo.degree.explvar  == 0 && expo.degree.conjvar  == 0
		realPart = Polynomial(); add!(realPart, a)
		imagPart = Polynomial(); add!(imagPart, b)
		return realPart, imagPart
	else
		realPart = Polynomial(); add!(realPart, 1)
		imagPart = Polynomial(); add!(imagPart, 0)

		vars, inds = collect(keys(expo.expo)), collect(values(expo.expo))
		return cplx2real_rec(vars, inds, realPart, imagPart, length(expo)+1, Degree(0,0))
	end
end

"""
realPart, imagPart = cplx2real_rec(vars::Array{Variable}, degs::Array{Degree}, realPart::Polynomial, imagPart::Polynomial, cur_ind::Int, cur_deg::Degree)

Transform recursively the complex exponent represented by the `vars` and `degs`
arrays into its real and imag parts, functions of its imag and real part
variables.
`cur_ind` decreases to 0, `cur_deg` decreases to Degree(0,0) for each step of
cur_ind. Terminaison case is reached at 0, Degree(0,0).
Initial arrays `vars` and `degs` are read only.

### Arguments
- vars::Array{Variable}
- degs::Array{Degree}
- realPart::Polynomial
- imagPart::Polynomial
- cur_ind::Int
- cur_deg::Degree
"""
function cplx2real_rec(vars::Array{Variable}, degs::Array{Degree}, realPart::Polynomial, imagPart::Polynomial, cur_ind::Int, cur_deg::Degree)
	## Final case:
	if cur_ind == 1 && cur_deg == Degree(0,0)
		return (realPart, imagPart)
	## One less variable to deal with:
	elseif cur_deg == Degree(0,0)
		return cplx2real_rec(vars, degs, realPart, imagPart, cur_ind-1, degs[cur_ind-1])
	## Recursion rule, decrease the current variable exponent until it reaches Degree(0,0):
	else
		var = vars[cur_ind]
		if iscomplex(var)
			var_real, var_imag = varname_cplx2real(var.name)
			var_R, var_I = Variable(var_real, Real), Variable(var_imag, Real)
			if cur_deg.explvar > 0
				realPart_new = product(var_R, realPart)
				add!(realPart_new, product(-1, product(var_I, imagPart)))

				imagPart_new = product(var_R, imagPart)
				add!(imagPart_new,  product(var_I, realPart))

				return cplx2real_rec(vars, degs, realPart_new, imagPart_new, cur_ind, Degree(cur_deg.explvar-1, cur_deg.conjvar))
			elseif cur_deg.conjvar > 0
				cur_deg.explvar == 0 || warn(LOGGER, "cur_deg.explvar should be 0 (and not $(cur_deg.explvar)), set to this value")
				realPart_new = product(var_R, realPart)
				add!(realPart_new,  product(var_I, imagPart))

				imagPart_new = product(var_R, imagPart)
				add!(imagPart_new, product(-1, product(var_I, realPart)))

				return cplx2real_rec(vars, degs, realPart_new, imagPart_new, cur_ind, Degree(0, cur_deg.conjvar-1))
			end
		elseif isbool(var)
			return cplx2real_rec(vars, degs, var*realPart, var*imagPart, cur_ind, Degree(0,0))
		else
			return cplx2real_rec(vars, degs, product(Exponent(SortedDict(var=>cur_deg)), realPart),
											 product(Exponent(SortedDict(var=>cur_deg)), imagPart),
											 cur_ind,
											 Degree(0,0))
		end
	end
end

"""
	realPart, imagPart = cplx2real(pol::Polynomial)

	Convert a complex polynomial in complex variables into `realPart` and
	`imagPart` polynomials of twice as many variables, real and imag parts of
	`pol` variables.
"""
function cplx2real(pol::Polynomial)
  realPart = Polynomial()
  imagPart = Polynomial()
  for (expo, λ) in pol
	if (expo.degree.explvar  == 1 && expo.degree.conjvar  == 1 )||( expo.degree.explvar  == 0 && expo.degree.conjvar  == 0)
		realexpo, imagexpo = cplx2real(expo, λ)

		add!(realPart, realexpo)
		add!(imagPart, imagexpo)
	else
	    realexpo, imagexpo = cplx2real(expo, λ)

		add!(realPart, product(realexpo,   real(λ)))
		add!(realPart, product(imagexpo, - imag(λ)))
		add!(imagPart, product(imagexpo,   real(λ)))
		add!(imagPart, product(realexpo,   imag(λ)))
	end
  end
  return (realPart, imagPart)
end

function cplx2real(pt_C::Point)
	pt = Point()
	for (var, val) in pt_C
		if var.kind <: Complex
			var_real, var_imag = varname_cplx2real(var.name)
			if real(val) != 0
				pt[Variable(var_real, Real)] = real(val)
			end
			if imag(val) != 0
				pt[Variable(var_imag, Real)] = imag(val)
			end
		else
			if val != 0
				pt[var] = real(val)
			end
		end
	end
	return pt
end

function real2cplx(pt::Point)
	ptC = Point()
	for (var, val) in pt
		if ismatch(r"_Re$", var.name)
			var_c = Variable(var.name[1:end-3], Complex)
			if !haskey(ptC, var_c)
				ptC[var_c] = val
			else
				ptC[var_c] += val
			end

		elseif ismatch(r"_Im$", var.name)
			var_c = Variable(var.name[1:end-3], Complex)
			if !haskey(ptC, var_c)
				ptC[var_c] = val*im
			else
				ptC[var_c] += val*im
			end

		else
			@assert !iscomplex(var)
			ptC[var] = val
		end
	end
	return ptC
end


text:
abstract type CoordinateSystem end

(tr::Transformed{<:CoordinateSystem})(p) = tr.obj(tr.transf(p))

struct Cartesian <: CoordinateSystem end

(c::Cartesian)(p) = p


text:
module ArarForecast

using Statistics
using TimeSeries
using Dates

include("arar.jl")

export arar

end


text:
using Random
using DABenchmarks

Random.seed!(0)

using Distributed
machine = [("MACHINENAME", 8)]
addprocs(machine, dir = "/home/USER", exename="/path/to/julia", exeflags="--project")
# addprocs(8, exeflags="--project") # For local workers

@everywhere using StochasticPrograms
@everywhere using LShapedSolvers
@everywhere using HydroModels
@everywhere using Gurobi

function prepare_dbenchmark(nscenarios::Integer, nsamples::Integer; timeout::Integer = 1000)
    dtr_1 = LShapedSolver(GurobiSolver(OutputFlag=0), regularization = :tr, log=false, distributed = true, subsolver = ()->GurobiSolver(OutputFlag=0), κ=1.0, bundle = 1)
    dtr_10 = LShapedSolver(GurobiSolver(OutputFlag=0), regularization = :tr, log=false, distributed = true, subsolver = ()->GurobiSolver(OutputFlag=0), κ=1.0, bundle = 10)
    dtr_50 = LShapedSolver(GurobiSolver(OutputFlag=0), regularization = :tr, log=false, distributed = true, subsolver = ()->GurobiSolver(OutputFlag=0), κ = 1.0, bundle = 50)
    dtr_100 = LShapedSolver(GurobiSolver(OutputFlag=0), regularization = :tr, log=false, distributed = true, subsolver = ()->GurobiSolver(OutputFlag=0), κ = 1.0, bundle = 100)
    dtr_125 = LShapedSolver(GurobiSolver(OutputFlag=0), regularization = :tr, log=false, distributed = true, subsolver = ()->GurobiSolver(OutputFlag=0), κ=1.0, bundle = 125)
    solvers = [dtr_1, dtr_10, dtr_50, dtr_100, dtr_125]
    solvernames = ["1", "10", "50", "100", "125"]
    # Create Day-ahead benchmark
    return Scaling(solvernames, solvers, nscenarios, nsamples; timeout = timeout)
end
@info "Running bundled L-shaped benchmarks on 16 workers..."
dab = prepare_dbenchmark(1000, 10); warmup_benchmarks!(dab); run_benchmarks!(dab); save_results!(dab, "tr_bundle_8.json")
exit()

max_stars_repo_name: martinbiel/mpc-reproducibility

max_issues_repo_name: martinbiel/mpc-reproducibility

max_forks_repo_name: martinbiel/mpc-reproducibility


text:
# OMG WHAT IS THIS LANGUAGE?!?!?!
name = "The ManBearPig!"
code = "Something Else"
println("Hello World! My name is $name")
println("Hi, $code)
println("Just a new line to create conflict")

max_stars_repo_name: ricardoprins/manbearpig

max_issues_repo_name: ricardoprins/manbearpig

max_forks_repo_name: ricardoprins/manbearpig


text:
include("libavfilter_h.jl")

#include("avfilter.jl")
#include("avfiltergraph.jl")
#include("buffersrc.jl")
#include("vsrc_buffer.jl")


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule XLA_Tools_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("XLA_Tools")
JLLWrappers.@generate_main_file("XLA_Tools", UUID("5debfaa8-cf0e-5577-9f8f-ffadd6421daa"))
end  # module XLA_Tools_jll


text:
@add_problem constant function constant_Issue_166(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test}
    # Issue #166
    α = Variable(5)
    fix!(α, ones(5,1))

    # has const vexity, but not at the head
    c = (rand(5,5) * α) * ones(1,5)

    β = Variable(5)
    β.value = ones(5)

    problem = minimize(sum(c * β), [β >= 0])
    handle_problem!(problem)
    if test
        @test problem.optval ≈ evaluate(sum(c * β)) atol=atol rtol=rtol
        @test problem.optval ≈ 0.0 atol=atol rtol=rtol
        @test β.value ≈ zeros(5) atol=atol rtol=rtol
    end
end

@add_problem constant function constant_Issue_228(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test}
    x = Variable(2)
    y = Variable(2)
    fix!(x, [1 1]')
    prob = minimize(y'*(x+[2 2]'), [y>=0])
    handle_problem!(prob)
    if test
        @test prob.optval ≈ 0.0 atol=atol rtol=rtol
    end

    prob = minimize(x'*y, [y>=0])
    handle_problem!(prob)
    if test
        @test prob.optval ≈ 0.0 atol=atol rtol=rtol
    end
end

@add_problem constant function constant_Test_double_fix!(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test}
    x = Variable()
    y = Variable()
    fix!(x, 1.0)
    prob = minimize(y*x, [y >= x, x >= 0.5])
    handle_problem!(prob)
    if test
        @test prob.optval ≈ 1.0 atol=atol rtol=rtol
    end

    fix!(x, 2.0)
    handle_problem!(prob)
    if test
        @test prob.optval ≈ 4.0 atol=atol rtol=rtol
    end

    free!(x)
    fix!(y, 1.0)
    handle_problem!(prob)
    if test
        @test prob.optval ≈ 0.5 atol=atol rtol=rtol
    end
end

@add_problem constant function constant_fix!_and_multiply(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test}
    p = Variable()
    fix!(p, 1.0)
    x = Variable(2,2)
    prob = minimize( tr(p*x), [ x >= 1 ])
    handle_problem!(prob)
    if test
        @test prob.optval ≈ 2.0 atol=atol rtol=rtol
        @test evaluate( tr(p*x) ) ≈ 2.0 atol=atol rtol=rtol
    end
end

@add_problem constant function constant_fix!_with_complex_numbers(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test}
    x = ComplexVariable()
    fix!(x, 1.0 + im*1.0)
    y = Variable()
    prob = minimize( real(x*y), [ y >= .5, real(x) >= .5, imag(x) >= 0])
    handle_problem!(prob)
    if test
        @test prob.optval ≈ .5 atol=atol rtol=rtol
        @test evaluate(real(x*y)) ≈ .5 atol=atol rtol=rtol
        @test evaluate(y) ≈ 0.5 atol=atol rtol=rtol
    end

    free!(x)
    if test
        fix!(y)
    else # if we haven't solved the problem, we need to fix! it to a number
        fix!(y, 0.5)
    end
    handle_problem!(prob)
    if test
        @test prob.optval ≈ 0.25 atol=atol rtol=rtol
        @test evaluate(real(x*y)) ≈ 0.25 atol=atol rtol=rtol
        @test real(evaluate(x)) ≈ 0.5 atol=atol rtol=rtol
        @test evaluate(y) ≈ 0.5 atol=atol rtol=rtol
    end

    if test
        @test_throws DimensionMismatch fix!(x, rand(2))
        @test_throws DimensionMismatch fix!(x, rand(2,2))
    end
end

@add_problem constant function constant_fix!_with_vectors(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test}
    x = ComplexVariable(5)
    fix!(x, ones(5) + im*ones(5))
    y = Variable()
    prob = minimize( real(y*sum(x)), [ y >= .5, real(x) >= .5, imag(x) >= 0])
    handle_problem!(prob)
    if test
        @test prob.optval ≈ 2.5 atol=atol rtol=rtol
        @test evaluate(real(y*sum(x))) ≈ 2.5 atol=atol rtol=rtol
        @test evaluate(y) ≈ 0.5 atol=atol rtol=rtol
    end

    free!(x)
    if test
        fix!(y)
    else # if we haven't solved the problem, we need to fix! it to a number
        fix!(y, 0.5)
    end

    handle_problem!(prob)
    if test
        @test prob.optval ≈ 1.25 atol=atol rtol=rtol
        @test evaluate(real(y*sum(x))) ≈ 1.25 atol=atol rtol=rtol
        @test real(evaluate(x)) ≈ 0.5*ones(5) atol=atol rtol=rtol
        @test evaluate(y) ≈ 0.5 atol=atol rtol=rtol
    end

    if test
        @test_throws DimensionMismatch fix!(x, rand(5,5))
        @test_throws DimensionMismatch fix!(x, rand(4))
    end

end


text:
t = Template(;
    user="tomouellette",
    authors="Tom W Ouellette",
    julia=v"1.5.4",
    plugins=[
        License(; name="MIT"),
        Git(; manifest=true, ssh=true),
        TravisCI(; linux = true, osx = true, x64 = true, x86 = true),
    ],
)

generate("CanEvolve.jl",t)


text:
module Imagine

using Images
using SIUnits, SIUnits.ShortUnits, Compat

import Images.imread
import Base.convert

export imagine2nrrd, Micron

Micron = SIUnits.NonSIUnit{typeof(Meter),:µm}()
convert(::Type{SIUnits.SIQuantity},::typeof(Micron)) = Micro*Meter

function imread{S<:IO}(s::S, ::Type{Images.ImagineFile})
    h = parse_header(s, Images.ImagineFile)
    filename = s.name[7:end-1]
    basename, ext = splitext(filename)
    camfilename = basename*".cam"
    T = h["pixel data type"]
    sz = [h["image width"], h["image height"], h["frames per stack"], h["nStacks"]]
    if sz[4] == 1
        sz = sz[1:3]
        if sz[3] == 1
            sz = sz[1:2]
        end
    end
    havez = h["frames per stack"] > 1
    havet = h["nStacks"] > 1
    # Check that the file size is consistent with the expected size
    if !isfile(camfilename)
        warn("Cannot open ", camfilename)
        data = Array(T, sz[1], sz[2], sz[3], 0)
    else
        fsz = filesize(camfilename)
        n_stacks = sz[end]
        if fsz != sizeof(T)*prod(int64(sz))  # guard against overflow on 32bit
            warn("Size of image file is different from expected value")
            n_stacks = ifloor(fsz / sizeof(T) / prod(sz[1:end-1]))
        end
        if sizeof(T)*prod(int64(sz[1:end-1]))*n_stacks > typemax(Uint)
            warn("File size is too big to mmap on 32bit")
            n_stacks = ifloor(fsz / sizeof(T) / typemax(Uint))
        end
        if n_stacks < sz[end]
            println("Truncating to ", n_stacks, length(sz) == 4 ? " stacks" : " frames")
            sz[end] = n_stacks
        end
        sc = open(camfilename, "r")
        data = mmap_array(T, ntuple(length(sz), i->sz[i]), sc)
    end
    um_per_pixel = h["um per pixel"]*µm
    pstart = h["piezo"]["stop position"]
    pstop = h["piezo"]["start position"]
    if length(sz)>2
        dz = abs(pstart - pstop)/sz[3]
    else dz = 0.0 end

    props = @compat Dict(
        "spatialorder" => havez ? ["x", "l", "z"] : ["x", "l"],
        "colorspace" => "Gray",
        "pixelspacing" => havez ? [um_per_pixel, um_per_pixel, dz] : [um_per_pixel, um_per_pixel],
        "limits" => (uint16(0), uint16(2^h["original image depth"]-1)),
        "imagineheader" => h,
        "suppress" => Set(Any["imagineheader"])
    )
    if havet
        props["timedim"] = havez ? 4 : 3
    end
    Image(data, props)
end

abstract Endian
type LittleEndian <: Endian; end
type BigEndian <: Endian; end
const endian_dict = @compat Dict("l"=>LittleEndian, "b"=>BigEndian)
const nrrd_endian_dict = @compat Dict(LittleEndian=>"little",BigEndian=>"big")
parse_endian(s::ASCIIString) = endian_dict[lowercase(s)]

function parse_vector_int(s::String)
    ss = split(s, r"[ ,;]", false)
    v = Array(Int, length(ss))
    for i = 1:length(ss)
        v[i] = int(ss[i])
    end
    return v
end

const bitname_dict = @compat Dict(
  "int8"      => Int8,
  "uint8"     => Uint8,
  "int16"     => Int16,
  "uint16"    => Uint16,
  "int32"     => Int32,
  "uint32"    => Uint32,
  "int64"     => Int64,
  "uint64"    => Uint64,
  "float16"   => Float16,
  "float32"   => Float32,
  "single"    => Float32,
  "float64"   => Float64,
  "double"    => Float64)

parse_bittypename(s::ASCIIString) = bitname_dict[lowercase(s)]

function float64_or_empty(s::ASCIIString)
    if isempty(s)
        return NaN
    else
        return float64(s)
    end
end

function parse_quantity_or_empty(s::ASCIIString)
    if isempty(s)
        return NaN
    else
        return parse_quantity(s)
    end
end

_unit_string_dict = @compat Dict("um" => Micro*Meter, "s" => Second, "us" => Micro*Second, "MHz" => Mega*Hertz)
function parse_quantity(s::String, strict::Bool = true)
    # Find the last character of the numeric component
    m = match(r"[0-9\.\+-](?![0-9\.\+-])", s)
    if m == nothing
        error("String does not have a 'value unit' structure")
    end
    val = float64(s[1:m.offset])
    ustr = strip(s[m.offset+1:end])
    if isempty(ustr)
        if strict
            error("String does not have a 'value unit' structure")
        else
            return val
        end
    end
    val * _unit_string_dict[ustr]
end

# Read and parse a *.imagine file (an Imagine header file)
const compound_fields = Any["piezo", "binning"]
const field_key_dict = @compat Dict{String,Function}(
    "header version"               => float64,
    "app version"                  => identity,
    "date and time"                => identity,
    "rig"                          => identity,
    "byte order"                   => parse_endian,
    "stimulus file content"        => identity,  # stimulus info parsed separately
    "comment"                      => identity,
    "ai data file"                 => identity,
    "image data file"              => identity,
    "start position"               => parse_quantity,
    "stop position"                => parse_quantity,
    "bidirection"                  => x->bool(int(x)),
    "output scan rate"             => x->parse_quantity(x, false),
    "nscans"                       => int,
    "channel list"                 => parse_vector_int,
    "label list"                   => identity,
    "scan rate"                    => x->parse_quantity(x, false),
    "min sample"                   => int,
    "max sample"                   => int,
    "min input"                    => float64,
    "max input"                    => float64,
    "original image depth"         => int,
    "saved image depth"            => int,
    "image width"                  => int,
    "image height"                 => int,
    "number of frames requested"   => int,
    "nStacks"                      => int,
    "idle time between stacks"     => parse_quantity,
    "pre amp gain"                 => float64_or_empty,
    "EM gain"                      => float64_or_empty,
    "gain"                         => float64_or_empty,
    "exposure time"                => parse_quantity,
    "vertical shift speed"         => parse_quantity_or_empty,
    "vertical clock vol amp"       => float64,
    "readout rate"                 => parse_quantity_or_empty,
    "pixel order"                  => identity,
    "frame index offset"           => int,
    "frames per stack"             => int,
    "pixel data type"              => parse_bittypename,
    "camera"                       => identity,
    "um per pixel"                 => float64,
    "hbin"                         => int,
    "vbin"                         => int,
    "hstart"                       => int,
    "hend"                         => int,
    "vstart"                       => int,
    "vend"                         => int,
    "angle from horizontal (deg)"  => float64_or_empty)

function parse_header(s::IOStream, ::Type{Images.ImagineFile})
    headerdict = Dict{ASCIIString, Any}()
    for this_line = eachline(s)
        this_line = strip(this_line)
        if !isempty(this_line) && !ismatch(r"\[.*\]", this_line)
            # Split on =
            m = match(r"=", this_line)
            if m.offset < 2
                error("Line does not contain =")
            end
            k = this_line[1:m.offset-1]
            v = this_line[m.offset+1:end]
            if in(k, compound_fields)
                thisdict = Dict{ASCIIString, Any}()
                # Split on semicolon
                strs = split(v, r";")
                for i = 1:length(strs)
                    substrs = split(strs[i], r":")
                    @assert length(substrs) == 2
                    k2 = strip(substrs[1])
                    func = field_key_dict[k2]
                    v2 = strip(substrs[2])
                    try
                        thisdict[k2] = func(v2)
                    catch err
                        println("Error processing key ", k2, " with value ", v2)
                        rethrow(err)
                    end
                end
                headerdict[k] = thisdict
            else
                func = field_key_dict[k]
                try
                    headerdict[k] = func(v)
                catch err
                    println("Error processing key ", k, " with value ", v)
                    rethrow(err)
                end
            end
        end
    end
    return headerdict
end

function imagine2nrrd(sheader::IO, h::Dict{ASCIIString, Any}, datafilename = nothing)
    println(sheader, "NRRD0001")
    T = h["pixel data type"]
    if T<:FloatingPoint
        println(sheader, "type: ", (T == Float32) ? "float" : "double")
    else
        println(sheader, "type: ", lowercase(string(T)))
    end
    sz = [h["image width"], h["image height"], h["frames per stack"], h["nStacks"]]
    kinds = ["space", "space", "space", "time"]
    if sz[end] == 1
        sz = sz[1:3]
        kinds = kinds[[1,2,4]]
    end
    println(sheader, "dimension: ", length(sz))
    print(sheader, "sizes:")
    for z in sz
        print(sheader, " ", z)
    end
    print(sheader, "\nkinds:")
    for k in kinds
        print(sheader, " ", k)
    end
    print(sheader, "\n")
    println(sheader, "encoding: raw")
    println(sheader, "endian: ", nrrd_endian_dict[h["byte order"]])
    if isa(datafilename, String)
        println(sheader, "data file: ", datafilename)
    end
    sheader
end

function imagine2nrrd(nrrdname::String, h::Dict{ASCIIString, Any}, datafilename = nothing)
    sheader = open(nrrdname, "w")
    imagine2nrrd(sheader, h, datafilename)
    close(sheader)
end

end