text:
import JavaCall
include("java_types.jl")

struct JspritSolver
    params::Dict{String,String}
    JspritSolver() = new(Dict{String,String}())
end

# RADIAL_BEST("radial_best"),
# RADIAL_REGRET("radial_regret"),
# RANDOM_BEST("random_best"),
# RANDOM_REGRET("random_regret"),
# WORST_BEST("worst_best"),
# WORST_REGRET("worst_regret"),
# CLUSTER_BEST("cluster_best"),
# CLUSTER_REGRET("cluster_regret"),
# STRING_BEST("string_best"),
# STRING_REGRET("string_regret");
# FIXED_COST_PARAM("fixed_cost_param"),
# VEHICLE_SWITCH("vehicle_switch"),
# REGRET_TIME_WINDOW_SCORER("regret.tw_scorer"),
# REGRET_DISTANCE_SCORER("regret.distance_scorer"),
# INITIAL_THRESHOLD("initial_threshold"),
# ITERATIONS("iterations"),
# THREADS("threads"),
# RANDOM_REGRET_MIN_SHARE("random_regret.min_share"),
# RANDOM_REGRET_MAX_SHARE("random_regret.max_share"),
# RANDOM_BEST_MIN_SHARE("random_best.min_share"),
# RANDOM_BEST_MAX_SHARE("random_best.max_share"),
# RADIAL_MIN_SHARE("radial.min_share"),
# RADIAL_MAX_SHARE("radial.max_share"),
# CLUSTER_MIN_SHARE("cluster.min_share"),
# CLUSTER_MAX_SHARE("cluster.max_share"),
# WORST_MIN_SHARE("worst.min_share"),
# WORST_MAX_SHARE("worst.max_share"),
# THRESHOLD_ALPHA("threshold.alpha"),
# THRESHOLD_INI("threshold.ini"),
# THRESHOLD_INI_ABS("threshold.ini_abs"),
# INSERTION_NOISE_LEVEL("insertion.noise_level"),
# INSERTION_NOISE_PROB("insertion.noise_prob"),
# RUIN_WORST_NOISE_LEVEL("worst.noise_level"),
# RUIN_WORST_NOISE_PROB("worst.noise_prob"),
# FAST_REGRET("regret.fast"),
# MAX_TRANSPORT_COSTS("max_transport_costs"),
# CONSTRUCTION("construction"),
# BREAK_SCHEDULING("break_scheduling"),
# STRING_K_MIN("string_kmin"),
# STRING_K_MAX("string_kmax"),
# STRING_L_MIN("string_lmin"),
# STRING_L_MAX("string_lmax"),
# MIN_UNASSIGNED("min_unassigned"),
# PROPORTION_UNASSIGNED("proportion_unassigned");

function setup_jvm()
    cd(dirname(@__FILE__))
    JavaCall.addClassPath("./jsprit-core.jar")
    JavaCall.init()
end

################################################
############### Jsprit functions ###############
################################################

function build_routes(data::RvrpInstance, solver::JspritSolver,
                      computed_data::RvrpComputedData,
                      jsp_sol, jsp_types::JspritJavaTypes)
    routes = Route[]
    jsp_routes_collection = JavaCall.jcall(jsp_sol, "getRoutes", jsp_types.JCollection, ())
    jsp_routes = JavaCall.jcall(jsp_routes_collection, "toArray",
                     Array{JavaCall.JObject, 1}, ())

    nb_routes =  length(jsp_routes)
    for r_idx in 1:nb_routes

        r_id = string("route_", r_idx)
        jsp_route = jsp_types.JVehicleRoute(jsp_routes[r_idx].ptr)

        jsp_vehicle = JavaCall.jcall(jsp_route, "getVehicle",
                                     jsp_types.JVehicleI, ())
        jsp_vehicle_type = JavaCall.jcall(jsp_vehicle, "getType",
                                          jsp_types.JVTypeI, ())
        v_set_id = JavaCall.jcall(jsp_vehicle_type, "getTypeId",
                                  JavaCall.JString, ())

        sequence = Action[]
        # Start action
        jsp_start_act = JavaCall.jcall(jsp_route, "getStart",
                            jsp_types.JStartActivity, ())
        arr_time = JavaCall.jcall(jsp_start_act, "getArrTime", JavaCall.jdouble, ())
        jsp_loc = JavaCall.jcall(jsp_start_act, "getLocation",
                            jsp_types.JLocation, ())
        loc_id = JavaCall.jcall(jsp_loc, "getId", JavaCall.JString, ())
        name = JavaCall.jcall(jsp_start_act, "getName", JavaCall.JString, ())
        operation_type = 0
        if name in ["pickupShipment", "pickup"]
            operation_type = 1
        elseif name in ["deliverShipment", "deliver"]
            operation_type = 2
        end
        push!(sequence, Action("action_1", loc_id, 0, "", arr_time))
        jsp_action_list = JavaCall.jcall(jsp_route, "getActivities",
                                         jsp_types.JList, ())
        jsp_actions = JavaCall.jcall(jsp_action_list, "toArray",
                                     Array{JavaCall.JObject, 1}, ())
        nb_actions = length(jsp_actions)
        for action_idx in 1:nb_actions
            jsp_act = jsp_types.JJobActivity(jsp_actions[action_idx].ptr)
            jsp_loc = JavaCall.jcall(jsp_act, "getLocation",
                                     jsp_types.JLocation, ())
            loc_id = JavaCall.jcall(jsp_loc, "getId", JavaCall.JString, ())
            jsp_job = JavaCall.jcall(jsp_act, "getJob", jsp_types.JJob, ())
            req_id = JavaCall.jcall(jsp_job, "getId", JavaCall.JString, ())
            arr_time = JavaCall.jcall(jsp_act, "getArrTime", JavaCall.jdouble, ())
            req = data.requests[computed_data.request_id_2_index[req_id]]
            name = JavaCall.jcall(jsp_act, "getName", JavaCall.JString, ())
            if name in ["pickupShipment", "pickup"]
                tws = [tw.soft_range for tw in req.pickup_time_windows]
                tw = get_closest_tw(tws, arr_time)
                if arr_time < tw.lb
                    arr_time = tw.lb
                end
                operation_type = 1
            elseif name in ["deliverShipment", "deliver"]
                tws = [tw.soft_range for tw in req.delivery_time_windows]
                tw = get_closest_tw(tws, arr_time)
                if arr_time < tw.lb
                    arr_time = tw.lb
                end
                operation_type = 2
            end
            push!(sequence, Action(string("action_", action_idx+1), loc_id,
                                        operation_type, req_id, arr_time))
        end
        # End action -> Assumes for now that always returns to depot
        jsp_end_act = JavaCall.jcall(jsp_route, "getEnd",
                                     jsp_types.JEndActivity, ())
        jsp_loc = JavaCall.jcall(jsp_end_act, "getLocation",
                                 jsp_types.JLocation, ())
        loc_id = JavaCall.jcall(jsp_loc, "getId", JavaCall.JString, ())
        arr_time = JavaCall.jcall(jsp_end_act, "getArrTime", JavaCall.jdouble, ())
        name = JavaCall.jcall(jsp_end_act, "getName", JavaCall.JString, ())
        operation_type = 0
        if name in ["pickupShipment", "pickup"]
            operation_type = 1
        elseif name in ["deliverShipment", "deliver"]
            operation_type = 2
        end
        push!(sequence, Action(string("action_", nb_actions+2),
                                    loc_id, operation_type, "", arr_time))

        push!(routes, Route(
            r_id, v_set_id, sequence, 0
        ))
    end

    return routes
end

function jsprit_build_and_add_locations(locations::Vector{Location},
                                        vrp_builder, jsp_types::JspritJavaTypes)
    JLocationBuilder = jsp_types.JLocationBuilder
    jsp_locs = jsp_types.JLocation[]
    for l in locations
        jsp_loc_builder = JLocationBuilder(())
        JavaCall.jcall(jsp_loc_builder, "setId", JLocationBuilder,
                       (JavaCall.JString,), l.id)
        JavaCall.jcall(jsp_loc_builder, "setName", JLocationBuilder,
                       (JavaCall.JString,), l.id)
        JavaCall.jcall(jsp_loc_builder, "setIndex", JLocationBuilder,
                       (JavaCall.jint,), l.index-1)
        c = jsp_types.JCoord((JavaCall.jdouble, JavaCall.jdouble),
                             l.long_x, l.lat_y)
        JavaCall.jcall(jsp_loc_builder, "setCoordinate", JLocationBuilder,
                       (jsp_types.JCoord,), c)
        jsp_loc = jsp_types.JLocation((JLocationBuilder,), jsp_loc_builder)
        push!(jsp_locs, jsp_loc)
    end
    return jsp_locs
end

function jsprit_build_vehicle_type(v_cats::Vector{VehicleCategory},
                                   v_set::HomogeneousVehicleSet,
                                   computed_data::RvrpComputedData,
                                   jsp_types::JspritJavaTypes)
    # Assumptions:
    # 1. Properties are binary properties, thus can be transformed to skills
    # 1.1 Only the ids of the properties are useful

    vc = v_cats[computed_data.vehicle_category_id_2_index[v_set.vehicle_category_id]]

    JVTypeBuilder = jsp_types.JVTypeBuilder

    jsp_vt_builder = JVTypeBuilder((JavaCall.JString,), v_set.id)
    JavaCall.jcall(jsp_vt_builder, "setCostPerDistance", JVTypeBuilder,
                   (JavaCall.jdouble,), v_set.cost_periods[1].travel_distance_unit_cost)
    JavaCall.jcall(jsp_vt_builder, "setCostPerTransportTime", JVTypeBuilder,
                   (JavaCall.jdouble,), v_set.cost_periods[1].travel_time_unit_cost)
    JavaCall.jcall(jsp_vt_builder, "setCostPerServiceTime", JVTypeBuilder,
                   (JavaCall.jdouble,), v_set.cost_periods[1].service_time_unit_cost)
    JavaCall.jcall(jsp_vt_builder, "setCostPerWaitingTime", JVTypeBuilder,
                   (JavaCall.jdouble,), v_set.cost_periods[1].waiting_time_unit_cost)
    JavaCall.jcall(jsp_vt_builder, "setFixedCost", JVTypeBuilder,
                   (JavaCall.jdouble,), v_set.cost_periods[1].fixed_cost)
    for (k,v) in vc.capacities.of_vehicle
        cap_dim = computed_data.capacity_id_2_index[k]
        JavaCall.jcall(jsp_vt_builder, "addCapacityDimension", JVTypeBuilder,
                       (JavaCall.jint, JavaCall.jint), cap_dim-1, v)
    end
    for s in collect(keys(vc.properties.of_vehicle))
        JavaCall.jcall(jsp_vt_builder, "addSkill", JVTypeBuilder,
                       (JavaCall.JString), s)
    end

    jsp_vt = JavaCall.jcall(jsp_vt_builder, "build", jsp_types.JVType, ())
    return jsp_vt
end

function jsprit_add_vehicles(v_cats::Vector{VehicleCategory},
                             v_sets::Vector{HomogeneousVehicleSet},
                             locations::Vector{Location},
                             lgs::Vector{LocationGroup},
                             computed_data::RvrpComputedData,
                             vrp_builder, jsp_locs, jsp_types::JspritJavaTypes)

    loc_id_2_idx = computed_data.location_id_2_index
    JVehicleBuilder = jsp_types.JVehicleBuilder
    # Assumptions:
    # 1. A vehicle has a single departure point and a single arrival point
    # 2. JS working time window equals HomogeneousVehicleSet.work_periods[1].soft_range
    # 3. We use only soft ranges
    # 4. departure_depot_id is never empty

    for v_set in v_sets
        # JVType:
        jsp_vt = jsprit_build_vehicle_type(v_cats, v_set, computed_data,
                                           jsp_types)

        # For start location:
        lg_idx = computed_data.location_group_id_2_index[v_set.departure_location_group_id]
        start_depot_idx = computed_data.location_id_2_index[lgs[lg_idx].location_ids[1]]
        # For end location:
        lg_idx = computed_data.location_group_id_2_index[v_set.arrival_location_group_id]
        lg = lgs[lg_idx]
        if isempty(lg.location_ids)
            end_depot_idx = -1
            return_to_depot = false
        else
            end_depot_idx = computed_data.location_id_2_index[lg.location_ids[1]]
            return_to_depot = true
        end

        # Define start and end locations from jsp_locations
        start_loc = jsp_locs[start_depot_idx]
        end_loc = nothing
        if end_depot_idx != -1
            end_loc = jsp_locs[end_depot_idx]
        end

        early_start_time = v_set.work_periods[1].soft_range.lb
        late_arrival_time = v_set.work_periods[1].soft_range.ub

        # Add as many copies as there are in the class
        for idx in 1:v_set.nb_of_vehicles_range.soft_range.ub
            id = string(v_set.id, "_", idx)

            jsp_vehicle_builder = JavaCall.jcall(JVehicleBuilder, "newInstance",
                JVehicleBuilder, (JavaCall.JString,), id)

            JavaCall.jcall(jsp_vehicle_builder, "setType",
                JVehicleBuilder, (jsp_types.JVTypeI,), jsp_vt)
            JavaCall.jcall(jsp_vehicle_builder, "setEarliestStart",
                JVehicleBuilder, (JavaCall.jdouble,), early_start_time)
            JavaCall.jcall(jsp_vehicle_builder, "setLatestArrival",
                JVehicleBuilder, (JavaCall.jdouble,), late_arrival_time)
            JavaCall.jcall(jsp_vehicle_builder, "setReturnToDepot",
                JVehicleBuilder, (JavaCall.jboolean,), return_to_depot)
            JavaCall.jcall(jsp_vehicle_builder, "setStartLocation",
                JVehicleBuilder, (jsp_types.JLocation,), start_loc)
            if return_to_depot
                JavaCall.jcall(jsp_vehicle_builder, "setEndLocation",
                    JVehicleBuilder, (jsp_types.JLocation,), end_loc)
            end

            jsp_vehicle = JavaCall.jcall(jsp_vehicle_builder, "build",
                                         jsp_types.JVehicle, ())
            JavaCall.jcall(vrp_builder, "addVehicle", jsp_types.JVrpBuilder,
                           (jsp_types.JAbstractVehicle,), jsp_vehicle)
        end
    end
end

function jsprit_build_service(req::Request, jsp_loc, serv_time::Float64,
                              consumptions::Dict{String,Float64},
                              skills::Vector{String},
                              tws::Vector{FlexibleRange}, vrp_builder,
                              JSpecificService, JSpecificServiceBuilder,
                              computed_data::RvrpComputedData,
                              jsp_types::JspritJavaTypes)

    # Assumptions:
    # 2. Maximum time in vehicle is always infinite, so no need to set
    JService = jsp_types.JService
    JServiceBuilder = jsp_types.JServiceBuilder
    JLocation = jsp_types.JLocation
    JTW = jsp_types.JTW

    jsp_service_builder = JavaCall.jcall(JSpecificServiceBuilder, "newInstance",
        JSpecificServiceBuilder, (JavaCall.JString,), req.id)
    JavaCall.jcall(jsp_service_builder, "setServiceTime", JServiceBuilder,
                   (JavaCall.jdouble,), serv_time)
    JavaCall.jcall(jsp_service_builder, "setLocation", JServiceBuilder,
                   (JLocation,), jsp_loc)
    for tw in tws
        jsp_tw = JTW((JavaCall.jdouble, JavaCall.jdouble), tw.soft_range.lb,
                     tw.soft_range.ub)
        JavaCall.jcall(jsp_service_builder, "addTimeWindow", JServiceBuilder,
                       (JTW,), jsp_tw)
    end
    for (cap_id,cons) in consumptions
        cons_idx = computed_data.capacity_id_2_index[cap_id]
        JavaCall.jcall(jsp_service_builder, "addSizeDimension", JServiceBuilder,
                       (JavaCall.jint, JavaCall.jint), cons_idx-1, cons)
    end
    for s in skills
        JavaCall.jcall(jsp_service_builder, "addRequiredSkill",
                       JServiceBuilder, (JavaCall.JString), s)
    end
    return JavaCall.jcall(jsp_service_builder, "build", JSpecificService, ())
end

function jsprit_build_shipment(req::Request, jsp_pickup_loc,
                               jsp_delivery_loc,
                               consumptions::Dict{String,Float64},
                               skills::Vector{String},
                               vrp_builder, computed_data::RvrpComputedData,
                               jsp_types::JspritJavaTypes)
    # Assumptions:
    JShipment = jsp_types.JShipment
    JShipmentBuilder = jsp_types.JShipmentBuilder
    JLocation = jsp_types.JLocation
    JTW = jsp_types.JTW

    jsp_shipment_builder = JavaCall.jcall(JShipmentBuilder, "newInstance",
        JShipmentBuilder, (JavaCall.JString,), req.id)
    JavaCall.jcall(jsp_shipment_builder, "setPickupLocation", JShipmentBuilder,
                   (JLocation,), jsp_pickup_loc)
    JavaCall.jcall(jsp_shipment_builder, "setDeliveryLocation", JShipmentBuilder,
                   (JLocation,), jsp_delivery_loc)
    JavaCall.jcall(jsp_shipment_builder, "setPickupServiceTime",
        JShipmentBuilder, (JavaCall.jdouble,), req.pickup_service_time)
    JavaCall.jcall(jsp_shipment_builder, "setDeliveryServiceTime",
        JShipmentBuilder, (JavaCall.jdouble,), req.delivery_service_time)
    JavaCall.jcall(jsp_shipment_builder, "setMaxTimeInVehicle",
         JShipmentBuilder, (JavaCall.jdouble,), req.max_duration)
    for tw in req.pickup_time_windows
        jsp_tw = JTW((JavaCall.jdouble, JavaCall.jdouble), tw.soft_range.lb,
                     tw.soft_range.ub)
        JavaCall.jcall(jsp_shipment_builder, "addPickupTimeWindow",
                       JShipmentBuilder, (JTW,), jsp_tw)
    end
    for tw in req.delivery_time_windows
        jsp_tw = JTW((JavaCall.jdouble, JavaCall.jdouble), tw.soft_range.lb,
                     tw.soft_range.ub)
        JavaCall.jcall(jsp_shipment_builder, "addDeliveryTimeWindow",
                       JShipmentBuilder, (JTW,), jsp_tw)
    end
    for (cap_id,cons) in consumptions
        cons_idx = computed_data.capacity_id_2_index[cap_id]
        JavaCall.jcall(jsp_shipment_builder, "addSizeDimension",
            JShipmentBuilder, (JavaCall.jint, JavaCall.jint), cons_idx-1, cons)
    end
    for s in skills
        JavaCall.jcall(jsp_service_builder, "addRequiredSkill",
                       JServiceBuilder, (JavaCall.JString), s)
    end
    return JavaCall.jcall(jsp_shipment_builder, "build",
                          jsp_types.JShipment, ())
end

function jsprit_add_jobs(reqs::Vector{Request},
                         locations::Vector{Location},
                         lgs::Vector{LocationGroup},
                         p_spec_classes::Vector{ProductSpecificationClass},
                         computed_data::RvrpComputedData,
                         vrp_builder, jsp_locs, jsp_types::JspritJavaTypes)
    # Assumptions:
    # 1. Properties are binary properties, thus can be transformed to skills
    for req in reqs
        p_spec_class = p_spec_classes[computed_data.product_specification_class_id_2_index[req.product_specification_class_id]]
        skills = collect(keys(p_spec_class.property_requirements))
        c = get_capacity_consumptions(req, p_spec_classes, computed_data)
        # Decide if it is a shipment, pickup or delivery
        if req.request_type == 0 # Pickup and delivery
            pickup_loc_idx = computed_data.location_id_2_index[lgs[computed_data.location_group_id_2_index[req.pickup_location_group_id]].location_ids[1]]
            delivery_loc_idx = computed_data.location_id_2_index[lgs[computed_data.location_group_id_2_index[req.delivery_location_group_id]].location_ids[1]]
            jsp_job = jsprit_build_shipment(req, jsp_locs[pickup_loc_idx],
                 jsp_locs[delivery_loc_idx], c, skills, vrp_builder,
                 computed_data, jsp_types)
        else
            if req.request_type == 1 # Pickup only
                service_type = jsp_types.JPickup
                service_builder = jsp_types.JPickupBuilder
                loc_idx = computed_data.location_id_2_index[lgs[computed_data.location_group_id_2_index[req.pickup_location_group_id]].location_ids[1]]
                s = req.pickup_service_time
                tws = req.pickup_time_windows
            elseif req.request_type == 2 # Delivery only
                service_type = jsp_types.JDelivery
                service_builder = jsp_types.JDeliveryBuilder
                loc_idx = computed_data.location_id_2_index[lgs[computed_data.location_group_id_2_index[req.delivery_location_group_id]].location_ids[1]]
                s = req.delivery_service_time
                tws = req.delivery_time_windows
            end
            jsp_job = jsprit_build_service(req, jsp_locs[loc_idx], s, c,
                skills, tws, vrp_builder, service_type, service_builder,
                computed_data, jsp_types)
        end
        JavaCall.jcall(vrp_builder, "addJob", jsp_types.JVrpBuilder,
                       (jsp_types.JAbstractJob,), jsp_job)
    end

end

function jsprit_add_transport_costs(jsp_vrp_builder,
             locations::Vector{Location},
             travel_specifications::Vector{TravelSpecification},
             travel_periods::Vector{TravelPeriod},
             computed_data::RvrpComputedData,
             jsp_types::JspritJavaTypes)

    JCostsBuilder = jsp_types.JCostsBuilder

    nb_mats = length(travel_specifications)
    nb_periods = length(travel_periods)

    is_matrix_symetric = false
    # is_matrix_symetric = LinearAlgebra.issymetric(dist_mat)

    jsp_costs_builder = JavaCall.jcall(JCostsBuilder, "newInstance",
        JCostsBuilder, (JavaCall.jint, JavaCall.jint, JavaCall.jint,
        JavaCall.jboolean), length(locations), nb_mats, nb_periods,
        is_matrix_symetric)

    # Add the periods
    for p_idx in 1:nb_periods
        p = travel_periods[p_idx].period
        spec_id = travel_periods[p_idx].travel_specification_id
        mat_idx = computed_data.travel_specification_id_2_index[spec_id]
        JavaCall.jcall(jsp_costs_builder, "addTravelPeriod",
            JCostsBuilder, (JavaCall.jint, JavaCall.jdouble, JavaCall.jdouble,
            JavaCall.jint), p_idx-1, p.lb, p.ub, mat_idx-1)
    end

    for travel_spec in travel_specifications
        dist_mat = travel_spec.travel_distance_matrix
        time_mat = travel_spec.travel_time_matrix
        mat_idx = computed_data.travel_specification_id_2_index[travel_spec.id]
        for l1_idx in 1:length(locations)
            for l2_idx in 1:length(locations)
                if dist_mat != Array{Float64,2}(undef, 0, 0)
                    JavaCall.jcall(jsp_costs_builder, "addTransportDistance",
                        JCostsBuilder, (JavaCall.jint, JavaCall.jint,
                        JavaCall.jint, JavaCall.jdouble),
                        l1_idx-1, l2_idx-1, mat_idx-1, dist_mat[l1_idx, l2_idx])
                end
                if time_mat != Array{Float64,2}(undef, 0, 0)
                    JavaCall.jcall(jsp_costs_builder, "addTransportTime",
                        JCostsBuilder, (JavaCall.jint, JavaCall.jint,
                        JavaCall.jint, JavaCall.jdouble),
                        l1_idx-1, l2_idx-1, mat_idx-1, time_mat[l1_idx, l2_idx])
                end
            end
        end
    end

    jsp_costs = JavaCall.jcall(jsp_costs_builder, "build", jsp_types.JCosts, ())
    JavaCall.jcall(jsp_vrp_builder, "setRoutingCost", jsp_types.JVrpBuilder,
                   (jsp_types.JCostsI,), jsp_costs)
end

function jsprit_create_input(data::RvrpInstance,
                             computed_data::RvrpComputedData,
                             jsp_types::JspritJavaTypes)

    jsp_vrp_builder = jsp_types.JVrpBuilder(())
    FS = jsp_types.FleetSize
    # To set fleet size : get the possible values and set one of them
    # using function setFleetSize
    # println("fleetzise methods: ", JavaCall.listmethods(FS))
    vals = JavaCall.jcall(FS, "values", Vector{FS}, ())
    JavaCall.jcall(jsp_vrp_builder, "setFleetSize",  jsp_types.JVrpBuilder,
                   (FS,), vals[1]) # FINITE

    jsp_locs = jsprit_build_and_add_locations(data.locations, jsp_vrp_builder,
                                              jsp_types)
    jsprit_add_vehicles(
        data.vehicle_categories, data.vehicle_sets, data.locations,
        data.location_groups, computed_data, jsp_vrp_builder, jsp_locs,
        jsp_types
    )
    jsprit_add_jobs(
        data.requests, data.locations, data.location_groups,
        data.product_specification_classes, computed_data, jsp_vrp_builder,
        jsp_locs, jsp_types
    )
    jsprit_add_transport_costs(
        jsp_vrp_builder, data.locations, data.travel_specifications,
        data.travel_periods, computed_data, jsp_types
    )

    return JavaCall.jcall(jsp_vrp_builder, "build", jsp_types.JVrp, ())
end

function jsprit_transform_solution(data::RvrpInstance, solver::JspritSolver,
                                 computed_data::RvrpComputedData,
                                 jsp_sol, jsp_types::JspritJavaTypes)
    instance_id = string(data.id, "_jsprit_SOL_", rand(1:10000))
    problem_id = data.id
    cost = JavaCall.jcall(jsp_sol, "getCost", JavaCall.jdouble, ())
    routes = build_routes(data, solver, computed_data, jsp_sol, jsp_types)

    jsp_unassigned_collection = JavaCall.jcall(jsp_sol, "getUnassignedJobs",
                                               jsp_types.JCollection, ())
    jsp_unassigned = JavaCall.jcall(jsp_unassigned_collection, "toArray",
                                     Array{JavaCall.JObject, 1}, ())
    unassigned = String[]
    for idx in 1:length(jsp_unassigned)
        jsp_job = jsp_types.JJob(jsp_unassigned[idx].ptr)
        id = JavaCall.jcall(jsp_job, "getId", JavaCall.JString, ())
        push!(unassigned, id)
    end

    return RvrpSolution(instance_id, problem_id, cost, routes, unassigned)
end

function jsprit_create_algorithm(jsp_vrp, jsp_types::JspritJavaTypes,
                                 solver::JspritSolver)

    jsp_builder = JavaCall.jcall(jsp_types.JSPBuilder, "newInstance",
        jsp_types.JSPBuilder, (jsp_types.JVrp,), jsp_vrp)

    # Set params here
    for (k,v) in solver.params
        println("Setting param: ", k, " to value ", v)
        JavaCall.jcall(jsp_builder, "setProperty", jsp_types.JSPBuilder,
            (JavaCall.JString, JavaCall.JString), k, v)
    end

    return JavaCall.jcall(jsp_builder, "buildAlgorithm",
                          jsp_types.JAlgorithm, ())
end

function jsprit_search_best_solution(jsp_algo, jsp_types::JspritJavaTypes)
    jsp_sols = JavaCall.jcall(jsp_algo, "searchSolutions",
                              jsp_types.JCollection, ())
    return JavaCall.jcall(jsp_types.JSolutions, "bestOf", jsp_types.JSolution,
                          (jsp_types.JCollection,), jsp_sols)
end

############################ RVRP SOLVER functions #########################

function solve(data::RvrpInstance, solver::JspritSolver)
    if !check_initial_feasibility(data)
        id = string(data.id, "_jsprit_SOL_", rand(1:10000))
        return build_empty_sol(data, id)
    end
    if !JavaCall.isloaded()
        setup_jvm()
    end
    computed_data = build_computed_data(data)
    jsp_types = JspritJavaTypes()
    jsp_vrp = jsprit_create_input(data, computed_data, jsp_types)
    jsp_algo = jsprit_create_algorithm(jsp_vrp, jsp_types, solver)
    jsp_sol = jsprit_search_best_solution(jsp_algo, jsp_types)
    rvrp_sol = jsprit_transform_solution(data, solver, computed_data,
                                         jsp_sol, jsp_types)
    return rvrp_sol
end

function supported_features(::Type{JspritSolver})
    features = BitSet()

    # Location based features  #
    # A location in jsprit has no time window
    # union!(features, HAS_OPENING_TIME_WINDOWS)
    # union!(features, HAS_MULTIPLE_OPENING_TIME_WINDOWS)
    # union!(features, HAS_FLEXIBLE_OPENING_TIME_WINDOWS)

    # Product based features #
    # union!(features, HAS_PRODUCT_CONFLICT_CLASSES)
    # union!(features, HAS_PRODUCT_PROHIBITED_PREDECESSOR_CLASSES)
    # union!(features, HAS_PRODUCT_PICKUPONLY_SHARING_CLASSES)
    # union!(features, HAS_PRODUCT_DELIVERYONLY_SHARING_CLASSES)
    # union!(features, HAS_PRODUCT_SHIPMENT_SHARING_CLASSES)
    union!(features, HAS_PRODUCT_PROPERTIES_REQUIREMENTS)
    union!(features, HAS_PRODUCT_CAPACITY_CONSUMPTIONS)
    union!(features, HAS_MULTIPLE_PRODUCT_CAPACITY_CONSUMPTIONS)

    # Request based features
    union!(features, HAS_SHIPMENT_REQUESTS)
    union!(features, HAS_PICKUPONLY_REQUESTS)
    union!(features, HAS_DELIVERYONLY_REQUESTS)
    union!(features, HAS_MAX_DURATION)
    union!(features, HAS_PICKUP_TIME_WINDOWS)
    union!(features, HAS_DELIVERY_TIME_WINDOWS)
    union!(features, HAS_MULTIPLE_PICKUP_TIME_WINDOWS)
    union!(features, HAS_MULTIPLE_DELIVERY_TIME_WINDOWS)

    # VehicleCategory based features
    union!(features, HAS_VEHICLE_CAPACITIES)
    union!(features, HAS_VEHICLE_PROPERTIES)
    union!(features, HAS_MULTIPLE_VEHICLE_CAPACITIES)
    union!(features, HAS_MULTIPLE_VEHICLE_PROPERTIES)

    # HomogeneousVehicleSet based features
    union!(features, HAS_TRAVEL_TIME_UNIT_COST)
    union!(features, HAS_SERVICE_TIME_UNIT_COST)
    union!(features, HAS_WAITING_TIME_UNIT_COST)
    union!(features, HAS_TRAVEL_DISTANCE_UNIT_COST)
    union!(features, HAS_MAX_NB_VEHICLES)
    union!(features, HAS_WORKING_TIME_WINDOW)
    union!(features, HAS_FIXED_COST_PER_VEHICLE)
    union!(features, HAS_OPEN_DEPARTURE)
    union!(features, HAS_OPEN_ARRIVAL)
    union!(features, HAS_ARRIVAL_DIFFERENT_FROM_DEPARTURE)

    # Instance based features
    union!(features, HAS_WORK_PERIODS)
    union!(features, HAS_MUTLIPLE_TRAVEL_TIME_PERIODS)
    union!(features, HAS_MULTIPLE_VEHICLE_CATEGORIES)
    union!(features, HAS_MULTIPLE_VEHICLE_SETS)

    return [features]
end


text:
@testset "TTR" begin
    @test ttr(["a","b","a","c","b"]) == 0.6
end

ext: jl


text:
using QuasiCopula, LinearAlgebra, Random, GLM
using DataFrames, Statistics
using BenchmarkTools, Test

BLAS.set_num_threads(1)
Threads.nthreads()

p_fixed = 3    # number of fixed effects, including intercept
m = 1    # number of variance components
# true parameter values
Random.seed!(12345)
# try next
βtrue = rand(Uniform(-2, 2), p_fixed)
θtrue = [0.5]
rtrue = 10.0

d = NegativeBinomial()
link = LogLink()
D = typeof(d)
Link = typeof(link)
T = Float64

samplesize = 10000
ni = 25

gcs = Vector{NBCopulaVCObs{T, D, Link}}(undef, samplesize)

Γ = θtrue[1] * ones(ni, ni)

# for reproducibility I will simulate all the design matrices here
Random.seed!(12345)
X_samplesize = [randn(ni, p_fixed - 1) for i in 1:samplesize]

for i in 1:samplesize
    X = [ones(ni) X_samplesize[i]]
    η = X * βtrue
    μ = exp.(η)
    p = rtrue ./ (μ .+ rtrue)
    vecd = Vector{DiscreteUnivariateDistribution}(undef, ni)
    vecd = [NegativeBinomial(rtrue, p[i]) for i in 1:ni]
    nonmixed_multivariate_dist = NonMixedMultivariateDistribution(vecd, Γ)
    # simuate single vector y
    y = Vector{Float64}(undef, ni)
    res = Vector{Float64}(undef, ni)
    rand(nonmixed_multivariate_dist, y, res)
    V = [ones(ni, ni)]
    gcs[i] = NBCopulaVCObs(y, X, V, d, link)
end

# form VarLmmModel
gcm = NBCopulaVCModel(gcs)
# precompile
println("precompiling NB VCM fit")
gcm2 = deepcopy(gcm);
QuasiCopula.fit!(gcm2, maxBlockIter = 1);

fittime = @elapsed QuasiCopula.fit!(gcm)
@show fittime
@show gcm.β
@show gcm.θ
@show gcm.r
@show gcm.∇β
@show gcm.∇θ
@show gcm.∇r

loglikelihood!(gcm, true, true)
vcov!(gcm)
@show QuasiCopula.confint(gcm)
# mse and time under our model
mseβ, mser, mseθ = MSE(gcm, βtrue, rtrue, θtrue)
@show mseβ
@show mser
@show mseθ

@test mseβ < 0.01
@test mseθ < 0.01
@test mser < 0.1


text:
function solve(prob::NEqProblem, x, method::Anderson, options::NEqOptions)
    if !(mstyle(prob) === InPlace())
        throw(ErrorException("solve() not defined for OutOfPlace() with Anderson"))
    end
    function fixedfromnleq(F, x)
        F .= value(prob, F, x) .+ x
    end
    fixedpoint!(fixedfromnleq, x, method;
                # kwargs
                Gx = similar(x),
                Fx = similar(x),
                f_abstol=sqrt(eps(eltype(x))),
                maxiter=options.maxiter)

end


text:
using Estimation
using Base.Test

include("frequency.jl")
include("APES.jl")


text:
using BinaryProvider

const JLTRILINOS_DIR = get(ENV, "JLTRILINOS_DIR", "")
const prefix = Prefix(JLTRILINOS_DIR == "" ? !isempty(ARGS) ? ARGS[1] : joinpath(@__DIR__,"usr") : JLTRILINOS_DIR)

products = Product[
    LibraryProduct(prefix, "libjltrilinos", :libjltrilinos)
]

if any(!satisfied(p; verbose=true) for p in products)
    # TODO: Add download code here
end

write_deps_file(joinpath(@__DIR__, "deps.jl"), products)


text:
using VectorSphericalWaves

n_max = 200

println("start testing `VectorSphericalWaves.γ_mn` ...")
for n = 1:n_max
    for m = -n:n
        try
            VectorSphericalWaves.γ_mn(m,n)
        catch e
            println("error for : n=$n, m=$m")
            """
            bt = backtrace()
            msg = sprint(showerror, e, bt)
            println(msg)
            """
            println(e)
            println()
            println()
        end
    end
end
println("end testing `VectorSphericalWaves.γ_mn` ...")
println()
println()
println("start testing `VectorSphericalWaves.γ_mn_dash` ...")
for n = 1:n_max
    for m = -n:n
        try
            VectorSphericalWaves.γ_mn_dash(m,n)
        catch e
            println("error for : n=$n, m=$m")
            """
            bt = backtrace()
            msg = sprint(showerror, e, bt)
            println(msg)
            """
            println(e)
            println()
            println()
        end
    end
end
println("end testing `VectorSphericalWaves.γ_mn_dash` ...")


text:
Marginalize = Command(
   "marginalize",
   "pangraph marginalize <options> [arguments]",
   "computes all pairwise marginalizations of a multiple sequence alignment graph",
   """multiple sequence alignment accepted in formats: [json]""",
   [
    Arg(
        String,
        "output path",
        (short="-o", long="--output-path"),
        "path to directory where all pairwise marginalizations will be stored\n\tif empty, will skip this computation",
        ""
    ),
    Arg(
        Bool,
        "reduce paralog paths",
        (short="-r", long="--reduce-paralog"),
        "collapse coparallel paths through duplications",
        false,
    ),
    Arg(
        String,
        "isolates to project onto",
        (short="-s", long="--strains"),
        "collapse the graph to only blocks contained by paths of the given isolates.\n\tcomma seperated list, no spaces",
        "",
    ),
   ],
   (args) -> let
       path = parse(Marginalize, args)
       path = if (path === nothing || length(path) == 0)
           nothing
       elseif length(path) == 1
           path
       else
           return 2
       end

       graph  = load(path, Marginalize)
       names  = collect(keys(graph.sequence))

       reduce = arg(Marginalize, "-r")
       output = arg(Marginalize, "-o")

       if length(output) > 0
           isdir(output) || mkpath(output)
           pairs  = [ (n₁,n₂) for n₁ in names for n₂ in names if n₁ < n₂ ]

           Threads.@threads for (name₁, name₂) in pairs
               G = deepcopy(graph)
               Graphs.keeponly!(G, name₁, name₂)
               if reduce
                   changed = true
                   while changed
                       l = length(G.block)
                       Graphs.deparalog!(G)
                       Graphs.detransitive!(G)
                       changed = l != length(G.block)
                   end
               else
                   Graphs.detransitive!(G)
               end

               Graphs.finalize!(G)
               open("$(output)/$(name₁)-$(name₂).json", "w") do io
                   marshal(io, G; fmt=:json)
               end
           end
       end

       isolates = arg(Marginalize, "-s")
       if length(isolates) > 0
           names = split(isolates,',')
           Graphs.keeponly!(graph, names...)

           if reduce
               changed = true
               while changed
                   l = length(graph.block)
                   Graphs.deparalog!(graph)
                   Graphs.detransitive!(graph)
                   changed = l != length(graph.block)
               end
           else
               Graphs.detransitive!(graph)
           end

           marshal(stdout, graph; fmt=:json)
       end
   end
)

max_stars_repo_name: neherlab/pangraph

max_issues_repo_name: neherlab/pangraph

max_forks_repo_name: neherlab/pangraph


text:
using Test
using LinearAlgebra
using SparseArrays
using Distributed

using jInv.Utils
using jInv.Mesh
using jInv.InverseSolve

# build domain and true image
domain = [0.0 1.0 0.0 1.0]
n      = [32,32]
Minv   = getRegularMesh(domain,n)
xc     = getCellCenteredGrid(Minv)
xtrue = sin.(2*pi*xc[:,1]).*sin.(pi*xc[:,2])

# get noisy data
A     = sparse(1.0I,Minv.nc,Minv.nc)
ids   = sort(rand(1:Minv.nc,round(Int64,Minv.nc*.8)))
A     = A[ids,:]
btrue = A*xtrue
bdata = btrue + .1*randn(size(btrue))/norm(btrue,Inf)
Wd    = ones(length(btrue))

# generate misfit param
i1     = (1:round(Int64,size(A,1)/4));
i2     = (round(Int64,size(A,1)/4)+1:round(Int64,2*size(A,1)/4));
i3 	   = (round(Int64,2*size(A,1)/4)+1 : round(Int64,3*size(A,1)/4));
i4 	   = (round(Int64,3*size(A,1)/4)+1 : size(A,1));
pFor1  = LSparam(A[i1,:],[])
pFor2  = LSparam(A[i2,:],[])

sigmaBack = zeros(length(xtrue))
Iact         = sparse(1.0I,Minv.nc,Minv.nc);
gl1          = getGlobalToLocal(1.0,sigmaBack)
gl2          = getGlobalToLocal(Iact,sigmaBack)

bd1          = bdata[i1]
bd2          = bdata[i2]
bd3          = bdata[i3]
bd4          = bdata[i4]
Wd1          = Wd[i1]
Wd2          = Wd[i2]
Wd3          = Wd[i3]
Wd4          = Wd[i4]

pMis =  getMisfitParam(pFor1,Wd1,bd1,SSDFun,fMod,gl1)

pMisRefs    = Array{RemoteChannel}(undef, 2)
pMisRefs[1] = initRemoteChannel(getMisfitParam,workers()[1],pFor1,Wd1,bd1,SSDFun,fMod,gl1)
pMisRefs[2] = initRemoteChannel(getMisfitParam,workers()[min(2,nworkers())],pFor2,Wd2,bd2,SSDFun,fMod,gl2)

pForRefs = Array{RemoteChannel}(undef, 4)
pForRefs[1] = initRemoteChannel(LSparam,workers()[1],A[i3,:],[])
pForRefs[2] = initRemoteChannel(LSparam,workers()[min(2,nworkers())],A[i4,:],[])
pForRefs[3] = initRemoteChannel(LSparam,workers()[min(3,nworkers())],A[i3,:],[])
pForRefs[4] = initRemoteChannel(LSparam,workers()[min(4,nworkers())],A[i4,:],[])

alpha        = 1.
x0           = (sum(xtrue)/length(xtrue))*ones(Minv.nc)

boundsLow    = minimum(xtrue)*ones(Minv.nc)
boundsHigh   = maximum(xtrue)*ones(Minv.nc)
pInv         = getInverseParam(Minv,fMod,diffusionReg,alpha,x0,boundsLow,boundsHigh)
display(pInv)

# Test with local misfit
# evaluate distance
sig,dsig = pInv.modelfun(x0)
Dc,F,dF,d2F,pMis,tMis = computeMisfit(sig,pMis,true)
# get Hessian
Hm = getHessian(sig,pMis,d2F)
# multiply with random direction and compare results
x = randn(size(Hm,2))
t1 = HessMatVec(x,pMis,sig,d2F)
t2 = Hm*x
@test norm(t1-t2,Inf)/norm(t1,Inf) < 1e-12

# Test with single remote ref misfit param
sig,dsig = pInv.modelfun(x0)
Dc,F,dF,d2F,pMisR,tMis = computeMisfit(sig,pMisRefs[1:1],true)
# get Hessian
Hm = getHessian(sig,pMisR[1],d2F[1])
# multiply with random direction and compare results
x = randn(size(Hm,2))
t1 = HessMatVec(x,pMisR,sig,d2F)
t2 = Hm*x
@test norm(t1-t2,Inf)/norm(t1,Inf) < 1e-12

# Test array of remote ref misfit param
# get Hessian
Hm = getHessian(sig,pMisR,d2F)
# multiply with random direction and compare results
x = randn(size(Hm,2))
t1 = HessMatVec(x,pMisR,sig,d2F)
t2 = Hm*x
@test norm(t1-t2,Inf)/norm(t1,Inf) < 1e-12

# Test array of misfit params
pMisA = [pMis;pMis;pMis]
sig,dsig = pInv.modelfun(x0)
Dc,F,dF,d2F,pMisA,tMis = computeMisfit(sig,pMisA,true)

# get Hessian
Hm = getHessian(sig,pMisA,d2F)
# multiply with random direction and compare results
x = randn(size(Hm,2))
t1 = HessMatVec(x,pMisA,sig,d2F)
t2 = Hm*x
@test norm(t1-t2,Inf)/norm(t1,Inf) < 1e-12


text:
# The first two consecutive numbers to have two distinct prime factors are:
#
# 14 = 2 × 7
# 15 = 3 × 5
#
# The first three consecutive numbers to have three distinct prime factors are:
#
# 644 = 2² × 7 × 23
# 645 = 3 × 5 × 43
# 646 = 2 × 17 × 19.
#
# Find the first four consecutive integers to have four distinct prime factors
# each. What is the first of these numbers?

using ProjectEulerSolutions

# Calculate prime factors iterating through each number, maintaining an array
# of four counts of prime factors.
function p047solution_upfront(n::Integer=100)::Integer
    q = zeros(Integer, 4)
    index = 1
    start = 645

    # Initialize it with first numbers
    for i in start:start+3
        q[index] = length(Set(primefactors(i)))
        index = mod(index, 4) + 1
    end

    # Start cycling through checking for every set for length 4
    for i in start+4:n
        q[index] = length(Set(primefactors(i)))
        index = mod(index, 4) + 1
        # This is faster than all(q .== 4)
        if q[1] == 4 && q[2] == 4 && q[3] == 4 && q[4] == 4
            return i - 3
        end
    end

    return -1
end

# Build up prime factors by iterating through.
function p047solution_iterative(n::Integer=100)::Integer
    factors = zeros(Integer, n)
    for i in 2:n-4
        if factors[i] == 0
            factors[i:i:n] .+= 1
        end
        # This is faster than all(factors[i:i+3] .== 4)
        if factors[i] == 4 && factors[i+1] == 4 && factors[i+2] == 4 && factors[i+3] == 4
            return i
        end
    end
    return -1
end

p047 = Problems.Problem(Dict("Upfront" => p047solution_upfront,
                             "Iterative" => p047solution_iterative))

Problems.benchmark(p047, 150_000)

max_stars_repo_name: gnujosh/julia-euler-project

max_issues_repo_name: gnujosh/julia-euler-project

max_forks_repo_name: gnujosh/julia-euler-project


text:
using BenchmarkTools
using OffsetArrays

const dim = 1000

x = Array{Float64}(undef, 2*dim)
y = OffsetArray{Float64}(undef, -dim + 1 : dim)
x2d = Array{Float64}(undef, 2*dim, 2*dim)
y2d = OffsetArray{Float64}(undef, -dim + 1 : dim, -dim + 1 : dim)

fill1d(x) = for i in axes(x,1); x[i] = i; end
fill2d(x) = for j in axes(x,2); for i in axes(x,1); x[i,j] = i + j; end; end
update(x) = for i in axes(x,1); x[i] = x[i] + i; end
update2d(x) = for j in axes(x,2); for i in axes(x,1); x[i,j] = x[i,j] + i + j; end; end
update_eachindex(x) = for i in eachindex(x); x[i] = x[i] + i; end

unsafe_fill(x) = @inbounds(for i in axes(x,1); x[i] = i; end)
unsafe_fill2d(x) = @inbounds(for j in axes(x,2); for i in axes(x,1); x[i,j] = i + j; end; end)
unsafe_update(x) = @inbounds(for i in axes(x,1); x[i] = x[i] + i; end)
unsafe_update2d(x) = @inbounds(for j in axes(x,2); for i in axes(x,1); x[i,j] = x[i,j] + i + j; end; end)
unsafe_update_eachindex(x) = @inbounds(for i in eachindex(x); x[i] = x[i] + i; end)

@show @benchmark fill1d(x)
@show @benchmark fill1d(y)
@show @benchmark fill2d(x2d)
@show @benchmark fill2d(y2d)
@show @benchmark update(x)
@show @benchmark update(y)
@show @benchmark update2d(x2d)
@show @benchmark update2d(y2d)
@show @benchmark update_eachindex(x)
@show @benchmark update_eachindex(y)

@show @benchmark unsafe_fill(x)
@show @benchmark unsafe_fill(y)
@show @benchmark unsafe_fill2d(x2d)
@show @benchmark unsafe_fill2d(y2d)
@show @benchmark unsafe_update(x)
@show @benchmark unsafe_update(y)
@show @benchmark unsafe_update2d(x2d)
@show @benchmark unsafe_update2d(y2d)
@show @benchmark unsafe_update_eachindex(x)
@show @benchmark unsafe_update_eachindex(y)


text:
countprimes(to; segment_length = 1024 * 32) = countprimes(1, to, segment_length)

function generate_siever_primes(small_sieve::SmallSieve, segment_lo)
    xs = small_sieve.xs
    sievers = Vector{Siever}(undef, vec_count_ones(xs))
    j = 0
    @inbounds for i = eachindex(xs)
        x = xs[i]
        while x != 0x00
            sievers[j += 1] = Siever(compute_prime(x, i), segment_lo)
            x &= x - 0x01
        end
    end
    return sievers
end

function countprimes(from, to, segment_length = 1024 * 32)
    # 1 is not a prime number anyway
    from = max(from, 2)

    first_byte = cld(from, 30)
    last_byte = cld(to, 30)

    sievers = generate_siever_primes(SmallSieve(isqrt(to)), 30 * (first_byte - 1) + 1)
    small_buffer = create_presieve_buffer()

    # Now create a chunk.
    xs = Vector{UInt8}(undef, segment_length)

    segment_index_start = first_byte

    count = if from <= 2
        3
    elseif from <= 5
        2
    elseif from <= 7
        1
    else
        0
    end

    @inbounds while segment_index_start <= last_byte
        segment_index_next = min(segment_index_start + segment_length, last_byte + 1)
        segment_curr_len = segment_index_next - segment_index_start

        # Presieve
        # fill!(xs, 0xFF)
        apply_presieve_buffer!(xs, small_buffer, segment_index_start, segment_index_next - 1)

        # @show vec_count_ones(xs, segment_curr_len)

        # Set the preceding so many bits before `from` to 0
        if segment_index_start == first_byte
            if first_byte === 1
                xs[1] = 0b11111110 # just make 1 not a prime.
            end
            for i = 1 : 8
                30 * (segment_index_start - 1) + ps[i] >= from && break
                xs[1] &= wheel_mask(ps[i])
            end
        end

        # Set the remaining so many bits after `to` to 0
        if segment_index_next == last_byte + 1
            for i = 8 : -1 : 1
                to >= 30 * (segment_index_next - 2) + ps[i] && break
                xs[segment_curr_len] &= wheel_mask(ps[i])
            end
        end

        # Sieve the interval, but skip the pre-sieved primes
        for p_idx in 5:length(sievers)
            p            = sievers[p_idx]
            last_idx     = 0
            n_bytes      = segment_index_next - segment_index_start
            byte_idx     = p.byte_index - segment_index_start + 1
            wheel_idx    = p.wheel_index
            increment    = p.prime_div_30
            @sieve_loop :unroll :save_on_exit
            sievers[p_idx] = Siever(increment, segment_index_start + byte_idx - 1, last_idx)
        end

        count += vec_count_ones(xs, segment_curr_len)

        segment_index_start += segment_length
    end

    return count
end


text:
using ..CommonFunctions

@doc raw"""
# Hermite's Interpolating Polynomials

Here, we implement a two point Hermite interpolation function which passes
through the function and its first derivative for the interval ``[0,1]``.
The polynomial is determined by four constraint equations, matching the
function and its derivative at the points ``0`` and ``1``.

Start by defining the 3rd degree polynomial and its derivative by
```math
\begin{align*}
g(x) &= a_0 + a_1 x + a_2 x^2 + a_3 x^3 , \\
g'(x) &= a_1 + 2 a_2 x + 3 a_3 x^2 ,
\end{align*}
```
and apply the constraints
```math
\begin{align*}
g(0) &= f_0 & & \Rightarrow & a_0 &= f_0 , \\
g(1) &= f_1 & & \Rightarrow & a_0 + a_1 + a_2 + a_3 &= f_1 , \\
g'(0) &= f'_0 & & \Rightarrow & a_1 &= f'_0 , \\
g'(1) &= f'_1 & & \Rightarrow & a_1 + 2 a_2 + 3 a_3 &= f'_1 . \\
\end{align*}
```
Solving for ``a_0, a_1, a_2, a_3`` leads to
```math
\begin{align*}
a_0 &= f_0 , &
a_1 &= f'_0 , &
a_2 &= - 3 f_0 + 3 f_1 - 2 f'_0 - f'_1 , &
a_3 &= 2 f_0 - 2 f_1 + f'_0 + f'_1 ,
\end{align*}
```
so that the polynomial ``g(x)`` reads
```math
g(x) = f_0 + f'_0 x + (- 3 f_0 + 3 f_1 - 2 f'_0 - f'_1) x^2 + (2 f_0 - 2 f_1 + f'_0 + f'_1) x^3 .
```
The function and derivative values can be factored out, so that ``g(x)`` can be rewritten as
```math
g(x) = f_0 (1 - 3 x^2 + 2 x^3) + f_1 (3 x^2 - 2 x^3) + f'_0 (x - 2 x^2 + x^3) + f'_1 (- x^2 + x^3) ,
```
or in generic form as
```math
g(x) = f_0 a_0(x) + f_1 a_1(x) + f'_0 b_0(x) + f'_1 b_1(x) ,
```
with basis functions
```math
\begin{align*}
a_0 (x) &= 1 - 3 x^2 + 2 x^3 , &
b_0 (x) &= x - 2 x^2 + x^3 , \\
a_1 (x) &= 3 x^2 - 2 x^3 , &
b_1 (x) &= - x^2 + x^3 .
\end{align*}
```
The derivative ``g'(x)`` accordingly reads
```math
g'(x) = f_0 a'_0(x) + f_1 a'_1(x) + f'_0 b'_0(x) + f'_1 b'_1(x) ,
```
with
```math
\begin{align*}
a'_0 (x) &= - 6 x + 6 x^2 , &
b'_0 (x) &= 1 - 4 x + 3 x^2 , \\
a'_1 (x) &= 6 x - 6 x^2 , &
b'_1 (x) &= - 2 x + 3 x^2 .
\end{align*}
```
The basis functions ``a_0``and ``a_1`` are associated with the function
values at ``x_0`` and ``x_1``, respectively, while the basis functions
``b_0`` and ``b_1`` are associated with the derivative values at
``x_0`` and ``x_1``.
The basis functions satisfy the following relations,
```math
\begin{align*}
a_i (x_j) &= \delta_{ij} , &
b_i (x_j) &= 0 , &
a'_i (x_j) &= 0 , &
b'_i (x_j) &= \delta_{ij} , &
i,j &= 0, 1 ,
\end{align*}
```
where ``\delta_{ij}`` denotes the Kronecker-delta, so that
```math
\begin{align*}
g(0) &= f_0 , &
g(1) &= f_1 , &
g'(0) &= f'_0 , &
g'(1) &= f'_1 .
\end{align*}
```
"""
struct HermiteInterpolation{T} <: Interpolator{T}
    x₀::T
    x₁::T
    Δx::T

    function HermiteInterpolation{T}(x₀, x₁, Δx, d) where {T}
        new(x₀, x₁, Δx)
    end
end

function HermiteInterpolation(x₀::T, x₁::T, Δx::T, d::Int) where {T}
    HermiteInterpolation{T}(x₀, x₁, Δx, d)
end

function CommonFunctions.evaluate!(int::HermiteInterpolation{T}, y₀::Vector{T}, y₁::Vector{T}, f₀::Vector{T}, f₁::Vector{T}, x::T, y::Vector{T}) where {T}
    local a₀::T
    local a₁::T
    local b₀::T
    local b₁::T

    # Interpolate y values at required locations
    if x == int.x₀
        y .= y₀
    elseif x == int.x₁
        y .= y₁
    else
        a₁ = 3x^2 - 2x^3
        a₀ = 1 - a₁
        b₁ = x^2*(x-1)
        b₀ = x*(1-x)+b₁
        y .= a₀ .* y₀ .+ a₁ .* y₁ .+ b₀ .* int.Δx .* f₀ .+ b₁ .* int.Δx .* f₁
    end
end

function CommonFunctions.evaluate!(int::HermiteInterpolation{T}, y₀::Vector{T}, y₁::Vector{T}, f₀::Vector{T}, f₁::Vector{T}, x::T, y::Vector{T}, f::Vector{T}) where {T}
    local a₀::T
    local a₁::T
    local b₀::T
    local b₁::T

    evaluate!(int, y₀, y₁, f₀, f₁, x, y)

    # Interpolate f values at required locations
    if x == int.x₀
        f .= f₀
    elseif x == int.x₁
        f .= f₁
    else
        a₁ = (6x - 6x^2) / int.Δx
        a₀ = - a₁
        b₁ = x*(3x-2)
        b₀ = 1-2x+b₁
        f .= b₀ .* f₀ .+ b₁ .* f₁ .+ a₀ .* y₀ .+ a₁ .* y₁
    end
end


text:
module LinkedLists

export LinkedList, Node

mutable struct Node{T}
  data::T
  next::Union{Node{T},Nothing}
end

function Node(data::T) where {T}
    Node(data, nothing)
end

mutable struct LinkedList{T}
    head::Union{Node{T},Nothing}
    tail::Union{Node{T},Nothing}
end

LinkedList{T}() where {T} = LinkedList{T}(nothing, nothing)
LinkedList() = LinkedList{Any}(nothing, nothing)

function get_testlist(start::Int, finish::Int)
    head = current = Node(start)

    for i in (start+1):finish
        node = Node(i)
        current.next = node
        current = node
    end

    tail = current

    LinkedList(head, tail)
end

function Base.iterate(list::LinkedList, state::Union{Node,Nothing} = list.head)
    (state == nothing || list.head == nothing) && return nothing
    (state, state.next)
end

Base.eltype(::Type{LinkedList{T}}) where T = Node{T}

function Base.length(list::LinkedList)
    count = 0

    for _ in list
        count += 1
    end

    count
end

function Base.getindex(list::LinkedList, index::Int)
    count = 1

    for node in list
        count == index && return node
        count += 1
    end

    throw(BoundsError("$(count-1) element LinkedList", index))
end

function Base.firstindex(list::LinkedList)
   1
end

function Base.lastindex(list::LinkedList)
   length(list)
end

function Base.setindex!(list::LinkedList{T}, node::Node{T}, index::Int) where {T}
    (index > lastindex(list) || index < firstindex(list)) &&
        throw(BoundsError("$(length(list)) element LinkedList", index))

    if index < lastindex(list)
        node.next = list[index+1]
    else
        node.next = nothing
        list.tail = node
    end

    if index > firstindex(list)
        list[index-1].next = node
    else
        list.head = node
    end
end

Base.show(io::IO, node::Node{T}) where {T} = print(io, "[$(node.data)::$(T)]→")

function Base.show(io::IO, list::LinkedList{T}) where {T}
    print(io, "LinkedList::$(T)")
    isempty(list) || print(io, " ")
    for node in list
        print(io, "[$(node.data)]→")
    end
end

function Base.push!(list::LinkedList{T}, nodes::Vararg{Node{T}}) where {T}
    head_node = current = nodes[1]

    for (index, node) in enumerate(nodes)
        index == 1 && continue
        current.next = node
        current = node
    end
    list[end].next = head_node
    list.tail = head_node

    list
end

function withnodes(nodes::Vararg{Node{T}}) where {T}
    head_node = current = nodes[1]

    for (index, node) in enumerate(nodes)
        index == 1 && continue
        current.next = node
        current = node
    end

    LinkedList(head_node, current)
end

import Base: +

function (+)(l1::LinkedList{T}, l2::LinkedList{T}) where {T}
    _l1 = deepcopy(l1)
    _l2 = deepcopy(l2)

    isempty(_l1) && return _l2
    isempty(_l2) && return _l1

    _l1[end].next = _l2[1]
    _l1.tail = _l2.tail

    _l1
end

function Base.hash(list::LinkedList)
    h = ""

    for node in list
        h *= string(node.data)
    end

    hash(h)
end

import Base.==

function ==(l1::LinkedList{T}, l2::LinkedList{T}) where {T}
  hash(l1) == hash(l2)
end

function push(list::LinkedList{T}, nodes::Vararg{Node{T}}) where {T}
    list + withnodes(nodes...)
end

function Base.push!(list::LinkedList{T}, nodes::Vararg{Node{T}}) where {T}
    list[end].next = withnodes(nodes...)[1]
    list
end

function pushfirst(list::LinkedList{T}, nodes::Vararg{Node{T}}) where {T}
    withnodes(nodes...) + list
end

function Base.pushfirst!(list::LinkedList{T}, nodes::Vararg{Node{T}}) where {T}
    list.head = (withnodes(nodes...) + list)[1]
    list
end

function Base.insert!(list::LinkedList{T}, pos::Int, nodes::Vararg{Node{T}}) where {T}
    (pos < 1 || pos > length(list)) &&
        throw(BoundsError("$(length(list)) element LinkedList", pos))

    nodes_list = withnodes(nodes...)

    if pos < length(list)
        nodes_list[end].next = list[pos]
    end
    if pos > 1
        list[pos-1].next = nodes_list[1]
    else
        list.head = nodes_list[1]
    end
    if pos == length(list)
      list[pos].next = nodes_list[1]
      list.tail = nodes_list.tail
    end

    list
end

function Base.pop!(list::LinkedList)
    last_node = list[end]

    if length(list) == 1
      list.head = nothing
      list.tail = nothing
    else
      list[end-1].next = nothing
      list.tail = list[end]
    end

    last_node
end

function Base.popfirst!(list::LinkedList)
    first_node = list[1]
    list.head = first_node.next
    first_node.next = nothing
    first_node
end

function Base.delete!(list::LinkedList{T}, remove_node::Node{T}) where {T}
    for (index, node) in enumerate(list)
        node == remove_node || continue

        if index > 1 && index < length(list)
            list[index-1].next = list[index+1]
        elseif index == 1 && index < length(list)
            list.head = list[index+1]
        elseif index > 1 && index == length(list)
            list[index-1].next = nothing
            list.tail = list[index-1]
        end

        return node
    end
end

function Base.deleteat!(list::LinkedList, index::Int) where {T}
    delete!(list, list[index])
end

function Base.empty!(list::LinkedList)
    list.head = nothing

    list
end

function Base.reverse(list::LinkedList)
    _list = deepcopy(list)
    last_node = _list[end]
    first_node = _list.head

    for index in length(_list):-1:1
        if index == 1
            _list[1].next = nothing
            break
        end

        _list[index].next = _list[index-1]
    end
    _list.head = last_node
    _list.tail = first_node

    _list
end

end

max_stars_repo_path: Chapter 6/LinkedLists.jl

max_issues_repo_path: Chapter 6/LinkedLists.jl

max_forks_repo_path: Chapter 6/LinkedLists.jl


text:
using LinearAlgebra

struct ValueOne; end
ValueOne()

# Compute X <- a X + b I.
function matfun_axpby!(X,a,b,Y::UniformScaling)
    m,n=size(X)
    if ~(a isa ValueOne)
        rmul!(X,a)
    end
    @inbounds for i=1:n
        X[i,i]+=(b isa ValueOne) ? 1 : b
    end
end

@inline function exp_ps_m6_opt_rho2_7(A)
    T=promote_type(eltype(A),Float64)
    A_copy=similar(A,T); A_copy .= A;
    return exp_ps_m6_opt_rho2_7!(A_copy)
end

@inline function exp_ps_m6_opt_rho2_7!(A)
    T=promote_type(eltype(A),Float64) # Make it work for many 'bigger' types (matrices and scalars)
    # max_memslots=9
    n=size(A,1)
    # The first slots are precomputed nodes [:A]
    memslots2 = similar(A,T)
    memslots3 = similar(A,T)
    memslots4 = similar(A,T)
    memslots5 = similar(A,T)
    memslots6 = similar(A,T)
    memslots7 = similar(A,T)
    memslots8 = similar(A,T)
    memslots9 = similar(A,T)
    # Assign precomputed nodes memslots
    memslots1=A # overwrite A
    # Uniform scaling is exploited.
    # No matrix I explicitly allocated.
    # Computation order: Bb2 Ba2 B2 Ba3 Bb3 B3 Ba4 Bb4 B4 Ba5 Bb5 B5 Ba6 Bb6 B6 Ba7 Bb7 B7 T2k9
    # Computing Bb2 = x*I+x*A
    coeff1=0.729254790966769
    coeff2=0.36467163874605407
    memslots2 .= coeff2.*memslots1
    mul!(memslots2,true,I*coeff1,true,true)
    # Computing Ba2 = x*I+x*A
    coeff1=0.7280078546464023
    coeff2=0.3646378676083862
    memslots3 .= coeff2.*memslots1
    mul!(memslots3,true,I*coeff1,true,true)
    # Computing B2 with operation: mult
    mul!(memslots4,memslots3,memslots2)
    # Deallocating Ba2 in slot 3
    # Deallocating Bb2 in slot 2
    # Computing Ba3 = x*I+x*A+x*B2
    coeff1=0.19178054915267911
    coeff2=0.7596053740559234
    coeff3=0.722565841223338
    memslots2 .= coeff2.*memslots1 .+ coeff3.*memslots4
    mul!(memslots2,true,I*coeff1,true,true)
    # Computing Bb3 = x*I+x*A+x*B2
    coeff1=1.0908348608940441
    coeff2=0.26673427641075104
    coeff3=0.1599439900856757
    memslots3 .= coeff2.*memslots1 .+ coeff3.*memslots4
    mul!(memslots3,true,I*coeff1,true,true)
    # Computing B3 with operation: mult
    mul!(memslots5,memslots2,memslots3)
    # Deallocating Ba3 in slot 2
    # Deallocating Bb3 in slot 3
    # Computing Ba4 = x*I+x*A+x*B2+x*B3
    coeff1=-0.0050940656098555845
    coeff2=0.3344727118113129
    coeff3=0.9346441756079398
    coeff4=0.789504699569824
    memslots2 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5
    mul!(memslots2,true,I*coeff1,true,true)
    # Computing Bb4 = x*I+x*A+x*B2+x*B3
    coeff1=1.2113351965247428
    coeff2=0.08511915910470753
    coeff3=-0.0039775493520180115
    coeff4=0.006605177298124698
    memslots3 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5
    mul!(memslots3,true,I*coeff1,true,true)
    # Computing B4 with operation: mult
    mul!(memslots6,memslots2,memslots3)
    # Deallocating Ba4 in slot 2
    # Deallocating Bb4 in slot 3
    # Computing Ba5 = x*I+x*A+x*B2+x*B3+x*B4
    coeff1=1.3379066745135616e-7
    coeff2=-0.002462022732959875
    coeff3=0.03688690867443445
    coeff4=-0.20340771723284812
    coeff5=0.986721033920639
    memslots2 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5 .+ coeff5.*memslots6
    mul!(memslots2,true,I*coeff1,true,true)
    # Computing Bb5 = x*I+x*A+x*B2+x*B3+x*B4
    coeff1=2.3798863906382593e-5
    coeff2=1.0837308004804432e-6
    coeff3=2.2489586532561667e-6
    coeff4=5.0091806253298255e-9
    coeff5=4.380604588051657e-8
    memslots3 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5 .+ coeff5.*memslots6
    mul!(memslots3,true,I*coeff1,true,true)
    # Computing B5 with operation: mult
    mul!(memslots7,memslots2,memslots3)
    # Deallocating Ba5 in slot 2
    # Deallocating Bb5 in slot 3
    # Computing Ba6 = x*I+x*A+x*B2+x*B3+x*B4+x*B5
    coeff1=0.0006283647624375745
    coeff2=-0.0856504866246125
    coeff3=-0.06273050414017863
    coeff4=0.8231044911539176
    coeff5=0.954079109272231
    coeff6=3.689858409164275e-5
    memslots2 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5 .+ coeff5.*memslots6 .+ coeff6.*memslots7
    mul!(memslots2,true,I*coeff1,true,true)
    # Computing Bb6 = x*I+x*A+x*B2+x*B3+x*B4+x*B5
    coeff1=0.009818461198297107
    coeff2=0.0010390323767811615
    coeff3=0.000350687159799907
    coeff4=8.23863768936227e-5
    coeff5=4.618679223329589e-5
    coeff6=0.9999998069262623
    memslots3 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5 .+ coeff5.*memslots6 .+ coeff6.*memslots7
    mul!(memslots3,true,I*coeff1,true,true)
    # Computing B6 with operation: mult
    mul!(memslots8,memslots2,memslots3)
    # Deallocating Ba6 in slot 2
    # Deallocating Bb6 in slot 3
    # Computing Ba7 = x*I+x*A+x*B2+x*B3+x*B4+x*B5+x*B6
    coeff1=0.9919381605205452
    coeff2=0.26220671903598886
    coeff3=0.26972239418076976
    coeff4=0.5440682102221873
    coeff5=1.01874072744531
    coeff6=3.3499625112901884e-5
    coeff7=0.05804569963149164
    memslots2 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5 .+ coeff5.*memslots6 .+ coeff6.*memslots7 .+ coeff7.*memslots8
    mul!(memslots2,true,I*coeff1,true,true)
    # Computing Bb7 = x*I+x*A+x*B2+x*B3+x*B4+x*B5+x*B6
    coeff1=-0.05660472123179401
    coeff2=0.021982855815076786
    coeff3=0.032263911106579925
    coeff4=0.021016177521655365
    coeff5=-0.021400262174087047
    coeff6=0.0006732366416760824
    coeff7=1.0044248526362387
    memslots3 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5 .+ coeff5.*memslots6 .+ coeff6.*memslots7 .+ coeff7.*memslots8
    mul!(memslots3,true,I*coeff1,true,true)
    # Computing B7 with operation: mult
    mul!(memslots9,memslots2,memslots3)
    # Deallocating Ba7 in slot 2
    # Deallocating Bb7 in slot 3
    # Computing T2k9 = x*I+x*A+x*B2+x*B3+x*B4+x*B5+x*B6+x*B7
    coeff1=0.6845260319725605
    coeff2=0.1283513925026888
    coeff3=0.021508466841232733
    coeff4=0.055068936428514716
    coeff5=0.2914936642490679
    coeff6=1.0432353938528486e-5
    coeff7=0.00037551432005694267
    coeff8=1.0045369223085259
    # Smart lincomb recycle A
    memslots1 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5 .+ coeff5.*memslots6 .+ coeff6.*memslots7 .+ coeff7.*memslots8 .+ coeff8.*memslots9
    mul!(memslots1,true,I*coeff1,true,true)
    # Deallocating B2 in slot 4
    # Deallocating B3 in slot 5
    # Deallocating B4 in slot 6
    # Deallocating B5 in slot 7
    # Deallocating B6 in slot 8
    # Deallocating B7 in slot 9
    return memslots1 # Returning T2k9
end

max_stars_repo_name: matrixfunctions/GraphMatFunData

max_issues_repo_name: matrixfunctions/GraphMatFunData

max_forks_repo_name: matrixfunctions/GraphMatFunData


text:
"""
abstract type Geometry end
"""

struct NoGeometry{C} <: Geometry
    capacitance::C
end
capacitance(g::NoGeometry) = g.capacitance

######### Geometries ###########
# is_geometric(::Compartment{T}) = T
# is_geometric(::Component) = false
"""
Every geometry needs methods for:
    capacitance()
    maybe reversals as dictionary?
    initial_conditions(f::Species)
"""


text:
@testset "Unit tests" begin
    @testset "#659 Volume errors if data is not a cube" begin
        fig, ax, vplot = volume(1:8, 1:8, 1:10, rand(8, 8, 10))
        lims = AbstractPlotting.data_limits(vplot)
        lo, hi = extrema(lims)
        @test all(lo .<= 1)
        @test all(hi .>= (8,8,10))
    end

    include("conversions.jl")
    include("quaternions.jl")
    include("projection_math.jl")
    include("liftmacro.jl")
    include("makielayout.jl")
    include("figures.jl")
end


text:
module Sherlogs

    export Sherlog64,Sherlog32,Sherlog16,
            DrWatson64,DrWatson32,DrWatson16,
            LogBook,reset_logbook,reset_logbooks,
            get_logbook,get_logbooks,
            get_stacktrace,get_stacktraces, get_st_index,
            entropy

    import Base: (+), (-), (*), (/), (^)

    import StatsBase: entropy

    include("logbook.jl")
    include("Sherlog64.jl")
    include("Sherlog32.jl")
    include("Sherlog16.jl")
    include("DrWatson64.jl")
    include("DrWatson32.jl")
    include("DrWatson16.jl")
    include("stacktraces.jl")

end


text:
# [function](https://docs.julialang.org/en/v1/manual/functions/#man-functions-1) to calculate the volume of a sphere
function sphere_vol(r)
    # julia allows [Unicode names](https://docs.julialang.org/en/v1/manual/unicode-input/#Unicode-Input-1) (in UTF-8 encoding)
    # so either "pi" or the symbol π can be used
    return 4/3*pi*r^3
end

# functions can also be defined more succinctly
quadratic(a, sqr_term, b) = (-b + sqr_term) / 2a

# calculates x for 0 = a*x^2+b*x+c, [arguments types](https://docs.julialang.org/en/v1/manual/functions/#Further-Reading-1) can be defined in function definitions
function quadratic2(a::Float64, b::Float64, c::Float64)
    # unlike other languages 2a is equivalent to 2*a
    # a^2 is used instead of a**2 or pow(a,2)
    sqr_term = sqrt(b^2-4a*c)
    r1 = quadratic(a, sqr_term, b)
    r2 = quadratic(a, -sqr_term, b)
    # multiple values can be returned from a function using tuples
    # if the [return](https://docs.julialang.org/en/v1/manual/functions/#The-return-Keyword-1) keyword is omitted, the last term is returned
    r1, r2
end

vol = sphere_vol(3)
# @printf allows number formatting but does not automatically append the \n to statements, see below
using Printf
@printf "volume = %0.3f\n" vol
#> volume = 113.097

quad1, quad2 = quadratic2(2.0, -2.0, -12.0)
println("result 1: ", quad1)
#> result 1: 3.0
println("result 2: ", quad2)
#> result 2: -2.0

max_stars_repo_name: talgatomarov/JuliaByExample

max_issues_repo_name: talgatomarov/JuliaByExample

max_forks_repo_name: talgatomarov/JuliaByExample


text:
using Test
using GeoParams
@testset "Density.jl" begin

# This tests the MaterialParameters structure
CharUnits_GEO   =   GEO_units(viscosity=1e19, length=1000km);

# Define a linear viscous creep law
x1      =   ConstantDensity(ρ=2900kg/m^3)
@test x1.ρ.val == 2900kg/m^3

Nondimensionalize!(x1,CharUnits_GEO)
@test x1.ρ.val ≈ 2.9e-16

x2      =   PT_Density()
@test x2.α.val==3e-5/K
@test x2.ρ0.val==2900kg/m^3

Nondimensionalize!(x2,CharUnits_GEO)
@test x2.T0.val≈0.21454659702313156

# Compute with density
@test ComputeDensity(1.0,1.0, x2) ≈ 2.8419999999999996e-16
@test ComputeDensity(1.0,1.0, x1) ≈ 2.9e-16

# Read Phase diagram interpolation object
fname   =   "./test_data/Peridotite.in"
PD_data =   PerpleX_LaMEM_Diagram(fname);
@test PD_data.meltFrac(1500,1e7) ≈ 0.24368492372485706
@test PD_data.Rho(1500,1e7) ≈ 3042.836820256982
@test PD_data.meltRho(1500,1e7) ≈ 2662.227167592414
@test PD_data.rockRho(1500,1e7) ≈ 3165.467673917775

@test ComputeDensity(1e7, 1500, PD_data) ≈ 3042.836820256982

# Do the same but non-dimensionalize the result
CharDim  =  GEO_units();
PD_data1 =  PerpleX_LaMEM_Diagram(fname, CharDim=CharDim);

rho_ND   =  PD_data1.Rho(Nondimensionalize(1500K,CharDim),Nondimensionalize(1e8*Pa,CharDim))
Vp_ND    =  PD_data1.Vp(Nondimensionalize(1500K,CharDim),Nondimensionalize(1e8*Pa,CharDim))
Vs_ND    =  PD_data1.Vs(Nondimensionalize(1500K,CharDim),Nondimensionalize(1e8*Pa,CharDim))

# redimensionalize and check with value from original structure that did not use non-dimensionalization
@test   ustrip(Dimensionalize(rho_ND,kg/m^3,CharDim)) ≈ PD_data.Rho(1500,1e8)
@test   ustrip(Dimensionalize(Vp_ND, km/s,  CharDim)) ≈ PD_data.Vp(1500,1e8)
@test   ustrip(Dimensionalize(Vs_ND, km/s,  CharDim)) ≈ PD_data.Vs(1500,1e8)

# Test computation of density for the whole computational domain, using arrays
MatParam    =   Array{MaterialParams, 1}(undef, 3);
MatParam[1] =   SetMaterialParams(Name="Mantle", Phase=0,
                        CreepLaws= (PowerlawViscous(), LinearViscous(η=1e23Pa*s)),
                        Density   = PerpleX_LaMEM_Diagram("test_data/Peridotite.in"));

MatParam[2] =   SetMaterialParams(Name="Crust", Phase=1,
                        CreepLaws= (PowerlawViscous(), LinearViscous(η=1e23Pa*s)),
                        Density   = ConstantDensity(ρ=2900kg/m^3));

MatParam[3] =   SetMaterialParams(Name="UpperCrust", Phase=2,
                        CreepLaws= (PowerlawViscous(), LinearViscous(η=1e23Pa*s)),
                        Density   = PT_Density());

# test computing material properties
Phases              = ones(Int64,400,400)*0;
Phases[:,20:end] .= 1
Phases[:,300:end] .= 2

rho     = zeros(size(Phases))
T       =  ones(size(Phases))
P       =  ones(size(Phases))*10

ComputeDensity!(rho, Phases, P,T, MatParam)
@test sum(rho)/400^2 ≈ 2920.6148898225

# test computing material properties when we have PhaseRatios, instead of Phase numbers
PhaseRatio  = zeros(size(Phases)...,length(MatParam));
for i in CartesianIndices(Phases)
    iz = Phases[i]
    I = CartesianIndex(i,iz+1)
    PhaseRatio[I] = 1.0
end

ComputeDensity!(rho, PhaseRatio, P,T, MatParam)
@test sum(rho)/400^2 ≈ 2920.6148898225

end


text:
include("learners/learners.jl")
include("explorers/explorers.jl")
include("q_based_policy.jl")


text:
######## Checking to see if the values in Figure 2, of rates varying with temperature, align with our calcs

#### all species 0.01g
Z = 1
M = BodyMasses(Z=Z)

#### fertilisation level 3
I_K = 3

# T range
temps = [0:1:40;] .+ 273.15

######### 2a - K

K2 = map(x -> (carryingcapacity_BA(I_K, param, M, x)[1]), temps) #get b0 for interaction 2 -> 1
pltK2 = plot(temps .- 273.15, K2, c = :orange , lw = 3
    , ylabel = "K", xlabel = "temperature (C)"
    , label = "2a", leg = :top
)

# carryingcapacity_BA(I_K, param, M, Tmin)  # 6.469796324965876
# carryingcapacity_BA(I_K, param, M, T10)   # 2.229581330642422
# carryingcapacity_BA(I_K, param, M, T20)   # 0.8262686110014499
# carryingcapacity_BA(I_K, param, M, T30)   # 0.32693445475490657
# carryingcapacity_BA(I_K, param, M, Tmax)  # 0.1372514415844402
# # CORRECT #

######### 2b - r
r2 = map(x -> (growth_BA(param, M, x)[1])/1e-6, temps) #get b0 for interaction 2 -> 1
pltr2 = plot(temps .- 273.15, r2, c = :orange , lw = 3
    , ylabel = "r(x 10^-6)", xlabel = "temperature (C)"
    , label = "2b", leg = :topleft
)

# growth_BA(param, M, Tmin) # 4.293834089901734e-8
# growth_BA(param, M, T20)  # 4.900749719809166e-7
# growth_BA(param, M, T30)  # 1.4677418185488954e-6
# growth_BA(param, M, Tmax) # 4.098371477948643e-6
# # CORRECT #

######### 2c - x relative to r

xrel2 = map(x -> (metabolism_BA(param, M, x)[2]/growth_BA(param, M, x)[1]), temps) #get b0 for interaction 2 -> 1
pltxrel2 = plot(temps .- 273.15, xrel2, c = :orange , lw = 3
    , ylabel = "x_rel", xlabel = "temperature (C)"
    , label = "2c", leg = :top
)

# metabolism_BA(param, M, Tmin) # 3.6997852501357905e-8
# 3.6997852501357905e-8/4.293834089901734e-8            # 0.8616507234960403

# metabolism_BA(param, M, T10) # 1.0418637435157166e-7

# metabolism_BA(param, M, T20)  # 2.7338159268322533e-7
# 2.7338159268322533e-7/4.900749719809166e-7            # 0.5578362665169335

# metabolism_BA(param, M, T30)  # 6.731119190535962e-7

# metabolism_BA(param, M, Tmax) # 1.5646345778326476e-6
# 1.5646345778326476e-6/4.098371477948643e-6            # 0.3817698288823232
# # CORRECT #

######### 2d - y relative to x

yrel2 = map(x -> (max_ingestion_BA(param, M, Z, x)[2,1]/metabolism_BA(param, M, x)[2]), temps) #get b0 for interaction 2 -> 1
pltyrel2 = plot(temps .- 273.15, yrel2, c = :orange , lw = 3
    , ylabel = "y_rel", xlabel = "temperature (C)"
    , label = "2d", leg = :top
)

# max_ingestion_BA(param, M, Z, Tmin) # 2.92703e-6
# 2.92703e-6/3.6997852501357905e-8                      # 79.11351070694093

# max_ingestion_BA(param, M, Z, T10) # 6.58049e-6
# 6.58049e-6/1.0418637435157166e-7                      # 63.16075437844174

# max_ingestion_BA(param, M, Z, T20)  #  1.11074e-5
# 1.11074e-5/2.7338159268322533e-7                      # 40.62965575326955

# max_ingestion_BA(param, M, Z, T30)  #  1.41136e-5
# 1.41136e-5/6.731119190535962e-7                      # 20.96768694847047

# max_ingestion_BA(param, M, Z, Tmax) #  1.35309e-5
# 1.35309e-5/1.5646345778326476e-6                      # 8.647961761616685
# # CORRECT #

######### 2e - B0 relative to K

B0rel2 = map(x -> (half_saturation_BA(param, M, Z, x)[2,1])/(carryingcapacity_BA(I_K, param, M, x)[1]), temps) #get b0 for interaction 2 -> 1
pltB0rel2 = plot(temps .- 273.15, B0rel2, c = :orange , lw = 3
    , ylabel = "B0_rel", xlabel = "temperature (C)"
    , label = "2e", leg = :topleft
)

# half_saturation_BA(param, M, Z, Tmin) # 2.08984
# 2.08984 / 6.469796324965876                            # 0.3230148052629806

# half_saturation_BA(param, M, Z , T10) #  2.65656
# 2.65656 /2.229581330642422                             # 1.1915062094794946

# half_saturation_BA(param, M, Z , T20) # 2.63599
# 2.63599/0.8262686110014499                            # 3.1902337386447983

# half_saturation_BA(param, M, Z , T30) # 2.0392
# 2.0392 /0.32693445475490657                          # 6.237335864550373

# half_saturation_BA(param, M, Z , Tmax) # 1.22858
# 1.22858/ 0.1372514415844402                           # 8.951308531387262
# # CORRECT #

l = @layout [a ; b c; d e]
plot(pltK2, pltr2, pltxrel2, pltyrel2, pltB0rel2,layout = l)

######### 2f - biomass extremes
# time span for each simulation
tspan = (0.0, 315360000000.0)
# define extinction threshold
threshold = 10e-12
##### data frame to store outputs
df = DataFrame(
    Z = [],
    I_K = [],
    T = [],
    K = [],
    r = [],
    x_2 = [],
    y_21 = [],
    B0_21 = [],
    maxB1 = [],
    maxB2 = [],
    maxB3 = [],
    minB1 = [],
    minB2 = [],
    minB3 = [],
    sp_survival = []
)

### Loop over temperature gradient

for T in temps
    # generate parameters specific to this combo of Z, I_K and T
    p = ModelParameters(param, T, I_K, Z)

    # Calculate starting biomasses
    u0 = [p.K[1]/2, p.K[1]/2, p.K[1]/2]

    # Define the problem
    prob = ODEProblem(BEFW, u0, tspan, p)

    # Solve the problem
    sol = solve(prob, maxiters = 1e7)

    ###### Output metrics ######
    # turn vector output of biomasses into matrix
    matrix_u = hcat(sol.u...)' # 56×3 adjoint(::Matrix{Float64}) with eltype Float64:

    ### Biomass extremes
    maxB1 = maximum(matrix_u[2 * Int64(round(length(matrix_u[:,1])/3)):Int64(length(matrix_u[:,1])),1])
    maxB2 = maximum(matrix_u[2 * Int64(round(length(matrix_u[:,1])/3)):Int64(length(matrix_u[:,1])),2])
    maxB3 = maximum(matrix_u[2 * Int64(round(length(matrix_u[:,1])/3)):Int64(length(matrix_u[:,1])),3])

    minB1 = minimum(matrix_u[2 * Int64(round(length(matrix_u[:,1])/3)):Int64(length(matrix_u[:,1])),1])
    minB2 = minimum(matrix_u[2 * Int64(round(length(matrix_u[:,1])/3)):Int64(length(matrix_u[:,1])),2])
    minB3 = minimum(matrix_u[2 * Int64(round(length(matrix_u[:,1])/3)):Int64(length(matrix_u[:,1])),3])

    ### Number of surviving species
    # final biomasses
    final_u = matrix_u[length(sol.t),:]

    # count number of species remaining that are above extinction threshold
    survived = count(final_u .> threshold)

    ### Push outputs to DataFrame
    push!(df, [
        Z,
        I_K,
        T,
        p.K[1],
        p.r[1],
        p.x[2],
        p.y[2,1],
        p.B0[2,1],
        maxB1,
        maxB2,
        maxB3,
        minB1,
        minB2,
        minB3,
        survived
    ])

    # print some stuff - see how the simulation is progressing
    println(("Z = $Z", "I_K = $I_K", "T = $T"))
end

df
# plot
# initialise an empty plot
pltextr2 = Plots.plot([NaN], [NaN],
                label = "",
                xlims = (0,40),
                ylims = (0,6.5),
                foreground_colour_legend = nothing,
                xlabel = "temperature",
                ylabel = "biomass extremes",
                leg = :topleft)

tempsC = [0:1:40;]
plot!(pltextr2, tempsC, df.maxB1, ls = :dot, c = :green, lw = 3, label = "basal")
plot!(pltextr2, tempsC, df.minB1, ls = :dot, c = :green, lw = 3, label = "")
plot!(pltextr2, tempsC, df.maxB2, ls = :dot, c = :blue, lw = 3, label = "intermediate")
plot!(pltextr2, tempsC, df.minB2, ls = :dot, c = :blue, lw = 3, label = "")
plot!(pltextr2, tempsC, df.maxB3, ls = :dot, c = :red, lw = 3, label = "top")
plot!(pltextr2, tempsC, df.minB3, ls = :dot, c = :red, lw = 3, label = "")


text:
#################################################################
#
#    Common defs for samplers
#
#################################################################

export EmpMCTuner

######### sample type returned by samplers  ############
immutable MCMCSample
	ppars::Vector{Float64}                  # proposed parameter vector
	plogtarget::Float64                     # log target	at proposed parameter vector
  pgrads::Union(Vector{Float64}, Nothing) # grad of log target at proposed parameter vector
	pars::Vector{Float64}                   # parameter vector before proposal
	logtarget::Float64                      # log target before proposal
  grads::Union(Vector{Float64}, Nothing)  # grad of log target before proposal
	diagnostics::Dict{Any,Any}              # sampler-dependant diagnostic variables
end

MCMCSample(ppars::Vector{Float64}, plogtarget::Float64,
  pars::Vector{Float64}, logtarget::Float64, diagnostics::Dict{Any,Any}) =
  MCMCSample(ppars, plogtarget, nothing, pars, logtarget, nothing, diagnostics)

MCMCSample(ppars::Vector{Float64}, plogtarget::Float64, pars::Vector{Float64}, logtarget::Float64) =
  MCMCSample(ppars, plogtarget, pars, logtarget, Dict())

MCMCSample(ppars::Vector{Float64}, plogtarget::Float64, pgrads::Union(Vector{Float64}, Nothing),
  pars::Vector{Float64}, logtarget::Float64, grads::Union(Vector{Float64}, Nothing)) =
  MCMCSample(ppars, plogtarget, pgrads, pars, logtarget, grads, Dict())

######### Empirical tuner common for various samplers ##
immutable EmpiricalMCMCTuner <: MCMCTuner
  adaptStep::Int
  maxStep::Int
  targetPath::Float64
  targetRate::Float64
  verbose::Bool

  function EmpiricalMCMCTuner(adaptStep::Int, maxStep::Int, targetPath::Float64, targetRate::Float64, verbose::Bool)
    @assert adaptStep > 0 "Adaptation step size ($adaptStep) should be > 0"
    @assert maxStep > 0 "Adaptation step size ($maxStep) should be > 0"
    @assert 0 < targetRate < 1 "Target acceptance rate ($targetRate) should be between 0 and 1"
    new(adaptStep, maxStep, targetPath, targetRate, verbose)
  end
end

typealias EmpMCTuner EmpiricalMCMCTuner

EmpMCTuner(targetRate::Float64; adaptStep::Int=100, maxStep::Int=200, targetPath::Float64=1., verbose::Bool=false) =
  EmpMCTuner(adaptStep, maxStep, targetPath, targetRate, verbose)

# Sampling task launcher
spinTask(m::MCMCModel, s::MCMCSampler, r::MCMCRunner) = MCMCTask(Task(() -> SamplerTask(m, s, r)), m, s, r)


text:
export scorePkg

const REQUIRE_EXISTS = 20.
const REQUIRE_VERSION = 20.
const LICENSE_EXISTS = 10.
const LICENSE = 0.
const LICENSE_FILE = 0.
const TEST_EXISTS = 20.
const TEST_RUNTESTS = 10.
const TEST_PASSES = 20.
const TEST_NOWARNING = 10.
const TEST_TRAVIS = 10.
const MAX_PKG_SCORE = REQUIRE_EXISTS + REQUIRE_VERSION + LICENSE_EXISTS +
                      TEST_EXISTS + TEST_RUNTESTS + TEST_PASSES + TEST_NOWARNING + TEST_TRAVIS

const URL_EXISTS = 20.
const DESC_EXISTS = 20.
const REQUIRES_OK = 20.
const REQUIRES_FAILS = 0.
const REQUIRES_PASSES = 0.
const MAX_METADATA_SCORE = URL_EXISTS + DESC_EXISTS + REQUIRES_OK

macro scoreMsg(key, msg, fatal)
  esc(quote
    write(o, $msg)
    if features[$key]
      total_score += eval($key)
      write(o, "    - ✓ Passed (+$(eval($key)))\n")
    else
      write(o, "    - ✗ Failed!\n")
      fatal_error = fatal_error || (true && $fatal)
    end
  end)
end

function scorePkg(features, pkg_name, pkg_path, metadata_path, o = STDOUT)

  total_score = 0.
  max_score = 0.
  fatal_error = false

  write(o, "# $pkg_name\n\n")
  #write(o, "- Local path: $pkg_path\n")
  #write(o, "- METADATA path: $metadata_path\n")

  if pkg_path != ""
    max_score += MAX_PKG_SCORE
    write(o, "\n## Package Itself\n")
    write(o, "\n### REQUIRE file\n")
    @scoreMsg(:REQUIRE_EXISTS,  "- Requirement: packages must have a REQUIRE file\n", true)
    @scoreMsg(:REQUIRE_VERSION, "- Requirement: REQUIRE file specifies a Julia version\n", true)

    write(o, "\n### Licensing\n")
    @scoreMsg(:LICENSE_EXISTS,  "- Recommendation: Packages should have a license\n", false)
    if features[:LICENSE_EXISTS]
      write(o, " - License detected in $(features[:LICENSE_FILE]): $(features[:LICENSE])\n")
    end

    write(o, "\n### Testing\n")
    @scoreMsg(:TEST_EXISTS,   "- Recommendation: Packages should have tests\n", false)
    if features[:TEST_EXISTS]
      @scoreMsg(:TEST_PASSES,    "- Requirement: If tests exist, they should pass.\n", true)
      #@scoreMsg(:TEST_NOWARNING, "- Recommendation: If tests exist, they should not have any warnings.\n", false)
    end
    @scoreMsg(:TEST_RUNTESTS, "- Recommendation: Packages should have a test/runtests.jl file\n", false)
    @scoreMsg(:TEST_TRAVIS,   "- Recommendation: Packages should have TravisCI support\n", false)
  end

  if metadata_path != ""
    max_score += MAX_METADATA_SCORE
    write(o, "\n## Package METADATA Entry\n")
    write(o, "\n### url file\n")
    @scoreMsg(:URL_EXISTS, "- Requirement: Packages must have a url file\n", true)

    write(o, "\n### DESCRIPTION.md file\n")
    @scoreMsg(:DESC_EXISTS, "- Recommendation: Packages should have a DESCRIPTION.md file\n", false)

    write(o, "\n### requires files\n")
    @scoreMsg(:REQUIRES_OK, "- Requirement: Each package version requires file must specify a Julia version\n", true)
    if !features[:REQUIRES_OK]
      write(o, "- Failed versions:\n")
      for version in features[:REQUIRES_FAILS]
        write(o, "    - $version\n")
      end
      write(o, "- Passed versions:\n")
      for version in features[:REQUIRES_PASSES]
        write(o, "    - $version\n")
      end
    end
  end

  write(o, "\n---\n")
  write(o, "\n## Summary\n")

  write(o, " - Total score: $total_score out of $max_score\n")
  if fatal_error
    write(o, " - One or more requirements failed - please fix and try again.\n\n")
  end

  write(o, "\n---\n\n## Raw Features\n")
  for k in keys(features)
    if k == :REQUIRES_PASSES || k == :REQUIRES_FAILS
      write(o, "    $k == [$(join(features[k]," "))]\n")
    else
      write(o, "    $k == $(features[k])\n")
    end
  end

  return total_score / max_score
end


text:
function hand_tuned_step!(out::AbstractVector,
                          bc::BorderArray,
                          bu::BorderArray,
                          alg)
    c = bc.inner
    u = bu.inner
    i1 = first(eachindex(c))
    iN = last( eachindex(c))
    pm = alg.variant
    innerinds = (i1+1):(iN-1)

    @inbounds @simd for i in innerinds
        ∇u2 = (u[i+1] - u[i])^2
        K = get_K(pm, i)
        c[i] = c_from_grad2(pm, K, ∇u2)
    end
    cstencil = CStencil(pm)
    mapstencil_border!(cstencil, c, bu)

    @inbounds @simd for i in innerinds
        u₀ = u[i]
        u₊ = u[i+1]
        u₋ = u[i-1]
        c₀ = c[i]
        c₊ = c[i+1]

        ∇u = u₊ - u₀
        ∇c = c₊ - c₀
        Δu = u₊ - 2u₀ + u₋
        dudt = c₀ * Δu + ∇c * ∇u
        out[i] = u[i] + alg.lambda * dudt
    end
    pmstencil = PMStencil(alg.lambda)
    mapstencil_border!(pmstencil, out, bu, bc)
    out
end

function hand_tuned_step!(out::AbstractMatrix,
                          bc::BorderArray,
                          bu::BorderArray,
                          alg)
    c = bc.inner
    u = bu.inner
    iax, jax = axes(c)
    i1 = first(iax)
    iN = last(iax)
    j1 = first(jax)
    jN = last(jax)
    pm = alg.variant
    iinnerinds = (i1+1):(iN-1)
    jinnerinds = (j1+1):(jN-1)

    @inbounds @simd for j in jinnerinds
        for i in iinnerinds
            ∇u2 = (u[i+1,j] - u[i,j])^2 + (u[i,j+1] - u[i,j])^2
            K = get_K(pm, i,j)
            c[i,j] = c_from_grad2(pm, K, ∇u2)
        end
    end
    cstencil = CStencil(pm)
    mapstencil_border!(cstencil, c, bu)

    @inbounds @simd for j in jinnerinds
        for i in iinnerinds
            u00 = u[i,  j  ]
            up0 = u[i+1,j  ]
            um0 = u[i-1,j  ]
            u0p = u[i  ,j+1]
            u0m = u[i  ,j-1]
            c00 = c[i  ,j  ]
            cp0 = c[i+1,j  ]
            c0p = c[i  ,j+1]

            ∇c∇u = (cp0 - c00)*(up0 - u00) + (c0p - c00)*(u0p - u00)
            Δu = u0p + up0 -4u00 + u0m + um0
            dudt = c00 * Δu + ∇c∇u
            out[i,j] = u00 + alg.lambda * dudt
        end
    end
    pmstencil = PMStencil(alg.lambda)
    mapstencil_border!(pmstencil, out, bu, bc)
    out
end


text:
include("../src/AlphaShapes.jl")
using LinearAlgebraicRepresentation
Lar = LinearAlgebraicRepresentation

if VERSION < VersionNumber("1.0.0")
	using Base.Test
else
	using Test
end

include("deWall.jl")
include("delaunay.jl")
include("alpha_complex.jl")
include("geometry.jl")


text:
function setup_q_profile_2D_serial!(w_inv_q::Array{T,2}, freesurface::Bool, nsponge::Int64,
    dt::T, q_freq::T, q_min::T, q_interior::T) where T<:Float32
    # check for unphysical q pairing
    if (q_min < eps(T)) && (q_interior < eps(T))
        w_inv_q .= zero(T)
        return
    end
    nz,nx = size(w_inv_q)
    lqmin = log(q_min)
    lqmax = log(q_interior)
    qprof = zeros(Float32, nsponge)

    for ksponge = 1:nsponge
        dk = (ksponge - 1) / (nsponge - 1)
        lq = lqmin + dk * (lqmax - lqmin)
        qprof[ksponge] = T(exp(lq))
    end

    for kz = 1:nz
        ksz = (freesurface) ? (nz - 1 - (kz - 1)) : min((kz - 1), (nz - 1 - (kz - 1)))
        for kx = 1:nx
            ksx = min((kx - 1),   (nx - 1 - (kx - 1)))
            ksponge = min(ksx, ksz) + 1 # add one to get back to julia indexing
            @inbounds w_inv_q[kz,kx] = (dt * 2.0 * π * q_freq / q_interior)
            if ksponge ∈ 1:nsponge
                @inbounds w_inv_q[kz,kx] = dt * 2.0 * π * q_freq / qprof[ksponge]
            end
        end
    end
end

function make_qprof(nsponge::Int64, q_min::T, q_interior::T) where T<:Float32
    lqmin = log(q_min)
    lqmax = log(q_interior)
    qprof = zeros(Float32, nsponge)

    @threads for ksponge = 1:nsponge
        dk = (ksponge - 1) / (nsponge - 1)
        lq = lqmin + dk * (lqmax - lqmin)
        @inbounds qprof[ksponge] = exp(lq)
    end

    return qprof
end

function setup_q_profile_2D_threaded!(w_inv_q::Array{T,2}, freesurface, nsponge::Int64,
        dt::T, q_freq::T, q_min::T, q_interior::T) where T<:Float32
    # check for unphysical q pairing
    if (q_min < eps(T)) && (q_interior < eps(T))
        w_inv_q .= zero(T)
        return
    end

    function compute_w_inv_q_2D!(w_inv_q::Array{T,2}, qprof::Array{T,1}, freesurface::Bool, nsponge::Int64, q_freq::T, q_interior::T) where T<:Float32
        nz,nx = size(w_inv_q)

        @threads for kz = 1:nz
            ksz = (freesurface) ? (nz - 1 - (kz - 1)) : min((kz - 1), (nz - 1 - (kz - 1)))
            for kx = 1:nx
                ksx = min((kx - 1),   (nx - 1 - (kx - 1)))
                ksponge = min(ksx, ksz) + 1 # add one to get back to julia indexing
                @inbounds w_inv_q[kz,kx] = dt * 2.0 * π * q_freq / q_interior
                if ksponge ∈ 1:nsponge
                    @inbounds w_inv_q[kz,kx] = dt * 2.0 * π * q_freq / qprof[ksponge]
                end
            end
        end
    end

    qprof = make_qprof(nsponge, q_min, q_interior)
    compute_w_inv_q_2D!(w_inv_q, qprof, freesurface, nsponge, q_freq, q_interior)
end

function setup_q_profile_3D_serial!(w_inv_q::Array{T,3}, freesurface::Bool, nsponge::Int64,
    dt::T, q_freq::T, q_min::T, q_interior::T) where T<:Float32
    # check for unphysical q pairing
    if (q_min < eps(T)) && (q_interior < eps(T))
        w_inv_q .= zero(T)
        return
    end

    nz,ny,nx = size(w_inv_q)
    lqmin = log(q_min)
    lqmax = log(q_interior)
    qprof = zeros(Float32, nsponge)

    for ksponge = 1:nsponge
        dk = (ksponge - 1) / (nsponge - 1)
        lq = lqmin + dk * (lqmax - lqmin)
        qprof[ksponge] = T(exp(lq))
    end

    for kz = 1:nz
        ksz = (freesurface) ? (nz - 1 - (kz - 1)) : min((kz - 1), (nz - 1 - (kz - 1)))

        for ky = 1:ny
            ksy = min((ky - 1),   (ny - 1 - (ky - 1)))
            for kx = 1:nx
                ksx = min((kx - 1),   (nx - 1 - (kx - 1)))
                ksponge = min(ksx,min(ksy, ksz)) + 1 # add one to get back to julia indexing
                @inbounds w_inv_q[kz,ky,kx] = (dt * 2.0 * π * q_freq / q_interior)
                if ksponge ∈ 1:nsponge
                    @inbounds w_inv_q[kz,ky,kx] = dt * 2.0 * π * q_freq / qprof[ksponge]
                end
            end
        end
    end
end

function setup_q_profile_3D_threaded!(w_inv_q::Array{T,3}, freesurface, nsponge::Int64,
    dt::T, q_freq::T, q_min::T, q_interior::T) where T<:Float32

    # check for unphysical q pairing
    if (q_min < eps(T)) && (q_interior < eps(T))
        w_inv_q .= zero(T)
        return
    end

    function compute_w_inv_q_3D!(w_inv_q::Array{T,3}, qprof::Array{T,1}, freesurface::Bool, nsponge::Int64, q_freq::T, q_interior::T) where T<:Float32
        nz,ny,nx = size(w_inv_q)

        @threads for kz = 1:nz
            ksz = (freesurface) ? (nz - 1 - (kz - 1)) : min((kz - 1), (nz - 1 - (kz - 1)))

            for ky = 1:ny
                ksy = min((ky - 1),   (ny - 1 - (ky - 1)))
                for kx = 1:nx
                    ksx = min((kx - 1),   (nx - 1 - (kx - 1)))
                    ksponge = min(ksx,min(ksy, ksz)) + 1 # add one to get back to julia indexing
                    @inbounds w_inv_q[kz,ky,kx] = dt * 2.0 * π * q_freq / q_interior
                    if ksponge ∈ 1:nsponge
                        @inbounds w_inv_q[kz,ky,kx] = dt * 2.0 * π * q_freq / qprof[ksponge]
                    end
                end
            end
        end
    end

    qprof = make_qprof(nsponge, q_min, q_interior)
    compute_w_inv_q_3D!(w_inv_q, qprof, freesurface, nsponge, q_freq, q_interior)
end

max_stars_repo_name: ChevronETC/WaveFD

max_issues_repo_name: ChevronETC/WaveFD

max_forks_repo_name: ChevronETC/WaveFD


text:
struct DisplayCell{T}
    label::Symbol
    elem::T
end
Base.show(io::IO, c::DisplayCell) = print(io, "$(c.label)=$(c.elem)")
function Base.typeinfo_prefix(io::IO, x::AbstractArray{DisplayCell{T}}) where {T}
    # print [a=1, b=2] instead of DisplayCell[a=1, b=2]
    Base.typeinfo_prefix(io, similar(x, T))
end

function Base.show(io::IO, x::LArray{T,N,D,Syms}) where {T,N,D,Syms}
    cells = Array{DisplayCell{T}, N}(undef, size(x)...)
    for (i, label) in enumerate(Syms)
        cells[i] = DisplayCell(label, x[i])
    end
    show(io, cells)
end

function Base.show(io::IO, x::SLArray{S,T,N,L,Syms}) where {S,T,N,L,Syms}
    cells = Array{DisplayCell{T}, N}(undef, size(x)...)
    for (i, label) in enumerate(Syms)
        cells[i] = DisplayCell(label, x[i])
    end
    show(io, cells)
end


text:
test_corpus = joinpath(Base.Filesystem.pwd(), "test_corpus.txt")
@testset "Document read tests" begin
    test_tokens = String["Hello", "there", "This", "is", "a", "sentence"]
    @test tokens(SimpleDocument(test_corpus)) == test_tokens
    @test tokens(ComplexityDocument(test_corpus, "", "", "", "", "", 0, "", "", "")) == test_tokens
end


text:
"""
Likelihood-free inference for the blow-fly model was introduced by [Simon N.
Wood](http://dx.doi.org/10.1038/nature09319). We model here the discrete time
stochastic dynamics of the size ``N`` of an adult blowfly population as given in [section 1.2.3 of the supplementary
information](https://static-content.springer.com/esm/art%3A10.1038%2Fnature09319/MediaObjects/41586_2010_BFnature09319_MOESM302_ESM.pdf).
```math
N_{t+1} = P N_{t-τ} exp(-N_{t-τ}/N₀)eₜ + Nₜ exp(-δ ϵₜ)
```
where ``eₜ`` and ``ϵₜ`` are independent Gamma random deviates with
mean 1 and variance ``σ_p²`` and ``σ_d²``, respectively.
"""
module BlowFlyModel
using Random, Distributions, StatsBase

Base.@kwdef struct BlowFly{R,S}
    burnin::Int = 50
    T::Int = 1000
    rng::R = MersenneTwister()
    statistics::S = histogram_summary_statistics
end
function (m::BlowFly)(P, N₀, σd, σp, τ, δ)
    p1 = Gamma(1 / σp^2, σp^2)
    p2 = Gamma(1 / σd^2, σd^2)
    T = m.T + m.burnin + τ
    N = fill(180.0, T)
    for t = τ+1:T-1
        N[t+1] =
            P * N[t-τ] * exp(-N[t-τ] / N₀) * rand(m.rng, p1) +
            N[t] * exp(-δ * rand(m.rng, p2))
    end
    m.statistics(N[end-m.T+1:end])
end

histogram_summary_statistics(N) = fit(Histogram, N, 140:16:16140).weights

# We will use a normal prior on log-transformed parameters.
function parameter(logparams)
    lP, lN₀, lσd, lσp, lτ, lδ = logparams
    (
        P = round(exp(2 + 2lP)),
        N₀ = round(exp(4 + 0.5lN₀)),
        σd = exp(-.5 + lσd),
        σp = exp(-.5 + lσp),
        τ = round(Int, max(1, min(500, exp(2 + lτ)))),
        δ = exp(-1 + 0.4lδ),
    )
end
(m::BlowFly)(logparams) = m(parameter(logparams)...)

end # module

import .BlowFlyModel: BlowFly
default_parameters(::Type{BlowFly}) = [
    (log(29) - 2) / 2,
    (log(260) - 4) * 2,
    log(0.6) + 0.5,
    log(0.3) + 0.5,
    log(7) - 2,
    (log(0.2) + 1) / 0.4,
]
prior(::Type{BlowFly}) = MultivariateNormal(zeros(6), ones(6))
function ApproxInferenceProblem(
    ::Type{BlowFly};
    prior = prior(BlowFly),
    target = default_parameters(BlowFly),
    kwargs...,
)
    model = BlowFly(; kwargs...)
    ApproxInferenceProblem(model, model(target), prior, target)
end


text:
# This file is part of IntegerSequences.
# Copyright Peter Luschny. License is MIT.

(@__DIR__) ∉ LOAD_PATH && push!(LOAD_PATH, (@__DIR__))

module LandauConstants

using Nemo, NumberTheory, GeneralizedBinomial, Products
export V000120, V102376, V011371, V120738, V060818, V124399, V056982
export V001511, V046161, V001316, V163590, V327492, V005187, V001790
export V277233, V327495, V327493, V327491, V007814, V001803, V327494
export V006519, V000265, Bin, BitCount, EvenPart, OddPart

"""

* V277233(n) / 4^V005187(n) are the Landau constants. These constants are defined as G(n)`` = ∑_{j=0..n} g(j)^2`` with the normalized central binomial
```
    g(n) = (2n)! / (2^n n!)^2 = V001790(n) / V046161(n).
```
* V327495(n) / 4^V327492(n) are the rationals considered here. These numbers are defined as H(n)`` = ∑_{j=0..n} h(j)^2`` with the normalized swinging factorial
```
    h(n) = n! / (2^n ⌊n/2⌋!)^2 = V163590(n) / V327493(n).
```
* In particular, this means that we have the pure integer representations
```
    V277233(n) = ∑{k=0..n}(V001790(k) 2^{V005187(n) - V005187(k)})^2;
    V327495(n) = ∑{k=0..n}(V163590(k) 2^{V327492(n) - V327492(k)})^2.
```

* V000120, V102376, V011371, V120738, V060818, V124399, V056982, V001511, V046161, V001316, V163590, V327492, V005187, V001790, V277233, V327495, V327493, V327491, V007814, V001803, V327494, V006519, V000265, Bin, BitCount, EvenPart, OddPart
"""
const ModuleLandauConstants = ""

"""

Return the 2-adic valuation of ``n``.
"""
Bin(n::Int) = n == 0 ? 0 : valuation(n, 2)

"""

Return the even part of ``n``.
"""
EvenPart(n::Int) = 1 << Bin(n)

"""

Return the even part of ``n``.
"""
V006519(n::Int) = EvenPart(n)

"""

Return the odd part of ``n``.
"""
OddPart(n::Int) = n >> Bin(n)

"""

Return the odd part of ``n``.
"""
V000265(n::Int) = OddPart(n)

"""

Return the odd part of the swinging factorial ``n! / ⌊n/2⌋!^2``.
"""
function V163590(n::Int)
    (r = get(CacheV163590, n, 0)) > 0 && return r
    h = V163590(n - 1) << Bin(n)
    CacheV163590[n] = isEven(n) ? div(h, n) : h * n
end

const CacheV163590 = Dict{Int,Int}(0 => 1)

"""

Return the odd part of the swinging factorial of ``2n+1``.
"""
V001803(n) = V163590(2n+1)
#V001803(n) = sum(<<(V001790(k), V005187(n) - V005187(k)) for k ∈ 0:n)

"""

Return the odd part of the swinging factorial of ``2n``.
"""
V001790(n::Int) = V163590(2n)
# V001790(n) = Binomial(2n, n) >> BitCount(2n)

"""

Return valuation ``(n, 2)`` if ``4`` divides ``n`` else ``1 + (n mod 2)``.
"""
V327491(n::Int) = Divides(4, n) ? Bin(n) : mod(n, 2) + 1

"""

Return ``2n - bc(n) + mod(n, 2)`` where bc(n) is BitCount(n).
"""
V327492(n::Int) = 2n - BitCount(n) + mod(n, 2)

"""

Return ``2^{2n - bc(n) + mod(n, 2)}`` where bc(n) is BitCount(n).
"""
V327493(n::Int) = 2^(2n - BitCount(n) + mod(n, 2))
#V327493(n::Int) = 2^V327492(n)

"""

Return the 2-adic valuation of ``n``.
"""
V007814(n::Int) = Bin(n)

"""

Return the 2-adic valuation of ``2n``.
"""
V001511(n::Int) = Bin(2n)

"""

Return the number of ``1``'s in binary expansion of ``n``.
"""
BitCount(n::Int) = count_ones(n)
#sum(digits(n, base = 2))

"""

Return the number of ``1``'s in binary expansion of ``n``.
"""
V000120(n::Int) = BitCount(n)

"""

Return 2^BitCount(n). Should be called Glaisher's sequence since James Glaisher showed that odd binomial coefficients are counted by this sequence.
"""
V001316(n::Int) = 2^BitCount(n)

"""

Return 4^BitCount(n).
"""
V102376(n::Int) = 4^BitCount(n)

"""

Return ``n`` - BitCount(n).
"""
V011371(n::Int) = n - BitCount(n)

"""

Return ``2n`` - BitCount(n).
"""
V005187(n::Int) = 2n - BitCount(n)

"""

Return ``4n`` - BitCount(n).
"""
V120738(n::Int) = 4n - BitCount(n)

"""

Return ``2^{n - bc(n)}`` where bc(n) is BitCount(n).
"""
V060818(n::Int) = 2^(n - BitCount(n))

"""

Return ``2^{2n - bc(n)}`` where bc(n) is BitCount(n).
"""
V046161(n::Int) = 2^(2n - BitCount(n))
# V046161(n) = 2^V005187(n)

"""

Return ``4^{n - bc(n)}`` where bc(n) is BitCount(n).
"""
V124399(n::Int) = 4^(n - BitCount(n))

"""

Return ``4^{2n - bc(n)}`` where bc(n) is BitCount(n).
"""
V056982(n::Int) = 4^(2n - BitCount(n))
# V056982(n) = 4^V005187(n)

"""

Return the numerator of ``∑_{j=0..n} j!/(2^j ⌊j/2⌋!)^2``.
"""
V327494(n) = sum(<<(V163590(k), V327492(n) - V327492(k)) for k ∈ 0:n)

"""

Return the numerators of the Landau constants.
"""
V277233(n) = sum(<<(V001790(k), V005187(n) - V005187(k))^2 for k ∈ 0:n)

"""

Return the numerators of the swinging Landau constants.
"""
V327495(n) = sum(<<(V163590(k), V327492(n) - V327492(k))^2 for k ∈ 0:n)

#START-TEST-########################################################

using Test

function test()
    @testset "Landau" begin
        @test V006519(0) == 1
        @test V000265(0) == 0
        @test V327494(0) == 1
        @test V001803(0) == 1
        @test V007814(0) == 0
        @test V327491(0) == 0
        @test V327493(0) == 1
        @test V000120(0) == 0
        @test V102376(0) == 1
        @test V011371(0) == 0
        @test V120738(0) == 0
        @test V060818(0) == 1
        @test V124399(0) == 1
        @test V056982(0) == 1
        @test V001511(0) == 0
        @test V046161(0) == 1
        @test V001316(0) == 1
        @test V163590(0) == 1
        @test V327492(0) == 0
        @test V005187(0) == 0
        @test V001790(0) == 1
        @test V277233(0) == 1
        @test V327495(0) == 1

        @test V006519(10) == 2
        @test V000265(10) == 5
        @test V327494(10) == 403341
        @test V001803(10) == 969969
        @test V007814(10) == 1
        @test V327491(10) == 1
        @test V327493(10) == 262144
        @test V000120(10) == 2
        @test V102376(10) == 16
        @test V011371(10) == 8
        @test V120738(10) == 38
        @test V060818(10) == 256
        @test V124399(10) == 65536
        @test V056982(10) == 68719476736
        @test V001511(10) == 2
        @test V046161(10) == 262144
        @test V001316(10) == 4
        @test V001511(10) == 2
        @test V163590(10) == 63
        @test V327492(10) == 18
        @test V005187(10) == 18
        @test V001790(10) == 46189
        @test V277233(10) == 125233642041
        @test V327495(10) == 74796032037
    end
end

function demo() end

function perf()
    "V001511" |> println
    @time for n ∈ 0:100 V001511(n) end
    "V163590" |> println
    @time for n ∈ 0:100 V163590(n) end
    "V327492" |> println
    @time for n ∈ 0:100 V327492(n) end
    "V005187" |> println
    @time for n ∈ 0:100 V005187(n) end
    "V001790" |> println
    @time for n ∈ 0:100 V001790(n) end
    "V277233" |> println
    @time for n ∈ 0:100 V277233(n) end
    "V327495" |> println
    @time for n ∈ 0:100 V327495(n) end
end

function main()
    test()
    demo()
    perf()
end

main()

end # module

max_stars_repo_name: OpenLibMathSeq/Sequences

max_issues_repo_name: OpenLibMathSeq/Sequences