text:
export
    init_env_1,
    init_env_2,
    init_env_3

init_env_1() = construct_regular_factory_world(;
    n_obstacles_x=2,
    n_obstacles_y=2,
    obs_width = [4;4], # obs_w = 8/n
    obs_offset = [4;4],
    env_pad = [1;1],
    # env_offset = [1,1],
    env_scale = 1,
    exclude_from_free = true,
)
init_env_2() = construct_regular_factory_world(;
    n_obstacles_x=4,
    n_obstacles_y=4,
    obs_width = [2;2],
    obs_offset = [2;2],
    env_pad = [1;1],
    # env_offset = [1,1],
    env_scale = 1,
    exclude_from_free = true,
)
init_env_3() = construct_regular_factory_world(;
    n_obstacles_x=8,
    n_obstacles_y=8,
    obs_width = [1;1],
    obs_offset = [1;1],
    env_pad = [1;1],
    # env_offset = [1,1],
    env_scale = 1,
    exclude_from_free = true,
)

export
    pctapf_problem,
    pctapf_problem_1,
    pctapf_problem_2,
    pctapf_problem_3,
    pctapf_problem_4,
    pctapf_problem_5,
    pctapf_problem_6,
    pctapf_problem_7,
    pctapf_problem_8

function get_zero_initial_conditions(G,N)
    to0_ = Dict{Int,Float64}(v=>0.0 for v in get_all_root_nodes(G))
    tr0_ = Dict{Int,Float64}(i=>0.0 for i in 1:N)
    return to0_, tr0_
end

function pctapf_problem(r0,config;Δt_op=0)
    tasks = config.tasks
    ops = config.ops
    s0 = map(t->t.first,tasks)
    sF = map(t->t.second,tasks)
    project_spec, _ = empty_pctapf_problem(r0,s0,sF)
    for op in ops
        add_operation!(project_spec,construct_operation(project_spec,-1,
            op.inputs,op.outputs,Δt_op))
    end
    def = SimpleProblemDef(project_spec,r0,s0,sF,config.shapes)
end

function empty_pctapf_problem(r0,s0,sF)
    object_ICs = [OBJECT_AT(o,x) for (o,x) in enumerate(s0)] # initial_conditions
    object_FCs = [OBJECT_AT(o,x) for (o,x) in enumerate(sF)] # final conditions
    robot_ICs = [ROBOT_AT(r,x) for (r,x) in enumerate(r0)]
    project_spec = ProjectSpec(object_ICs,object_FCs)
    project_spec, robot_ICs
end

function pctapf_problem(
        solver,
        sched::OperatingSchedule,
        problem_spec::ProblemSpec,
        env_graph,
        args...
    )
    env = construct_search_env(
        solver,
        sched,
        problem_spec,
        env_graph
        )
    PC_TAPF(env)
end
function pctapf_problem(solver,spec::ProjectSpec,env,robot_ics,prob_spec=ProblemSpec(D=env))
    sched = construct_partial_project_schedule(spec,prob_spec,robot_ics)
    return pctapf_problem(solver,sched,prob_spec,env)
end
function PC_TAPF(solver,def::SimpleProblemDef,env::GridFactoryEnvironment)
    sched, prob_spec = construct_task_graphs_problem(def,env)
    pctapf = pctapf_problem(solver,sched,prob_spec,env)
    return pctapf
end

function pcta_problem(
        # project_spec::ProjectSpec,
        sched::OperatingSchedule,
        problem_spec::ProblemSpec,
        # robot_ICs,
        env_graph,
        primary_objective = MakeSpan()
    )
    cache=initialize_planning_cache(sched)
    cost_model = typeof(primary_objective)(sched,cache)
    heuristic_model = NullHeuristic()
    layers = construct_environment_layer_dict(sched,env_graph,problem_spec)
    env = SearchEnv(
        schedule=sched,
        cache=cache,
        env_layers=layers,
        cost_model = cost_model,
        heuristic_model = heuristic_model
        )
    PC_TA(env)
end

function PC_TA(def::SimpleProblemDef,env::GridFactoryEnvironment,objective=MakeSpan())
    sched, prob_spec, env, _ = construct_task_graphs_problem(def,env)
    pcta = pcta_problem(sched,prob_spec,env,objective)
end

# Reusable problem initializers for testing
"""
    pctapf_problem_1

Optimal MakeSpan = 5
Optimal SumOfMakeSpans = 5
"""
function pctapf_problem_1(;cost_function=SumOfMakeSpans(),verbose=false,Δt_op=0,Δt_collect=[0,0,0],Δt_deposit=[0,0,0])
    N = 2                  # num robots
    M = 3                  # num delivery tasks
    vtx_grid = initialize_dense_vtx_grid(4,4)
    #  1   2   3   4
    #  5   6   7   8
    #  9  10  11  12
    # 13  14  15  16
    r0 = [1,4]
    s0 = [5,8,14]
    sF = [13,12,15]
    env_graph = construct_factory_env_from_vtx_grid(
        vtx_grid;
    )
    project_spec, robot_ICs = empty_pctapf_problem(r0,s0,sF)
    add_operation!(project_spec,construct_operation(project_spec,-1,[1,2],[3],Δt_op))
    add_operation!(project_spec,construct_operation(project_spec,-1,[3],  [], Δt_op))
    assignment_dict = Dict(1=>[1,3],2=>[2])

    def = SimpleProblemDef(project_spec,r0,s0,sF)
    sched, problem_spec = construct_task_graphs_problem(
        def,env_graph;cost_function=cost_function,Δt_collect=Δt_collect,Δt_deposit=Δt_deposit)
    return sched, problem_spec, env_graph, assignment_dict
end

"""
    pctapf_problem_2(;cost_function=SumOfMakeSpans(),verbose=false)

In this problem robot 1 will first do [1-9-17], then [17-21-35]
robot 2 will do [4-12-32]. The key thing is that robot 1 will need to wait
until robot 2 is finished before robot 1 can do its second task.

Optimal paths:
Optimal MakeSpan = 8
Optimal SumOfMakeSpans = 8
"""
function pctapf_problem_2(;cost_function=SumOfMakeSpans(),verbose=false)
    N = 2                  # num robots
    M = 3                  # num delivery tasks
    vtx_grid = initialize_dense_vtx_grid(8,4)
    #  1   2   3   4
    #  5   6   7   8
    #  9  10  11  12
    # 13  14  15  16
    # 17  18  19  20
    # 21  22  23  24
    # 25  26  27  28
    # 29  30  31  32
    r0 = [1,4]
    s0 = [9,12,21]
    sF = [17,32,25]
    env_graph = construct_factory_env_from_vtx_grid(
        vtx_grid;
    )
    project_spec, robot_ICs = empty_pctapf_problem(r0,s0,sF)

    add_operation!(project_spec,construct_operation(project_spec,-1,[1,2],[3],0))
    add_operation!(project_spec,construct_operation(project_spec,-1,[3],  [], 0))
    assignment_dict = Dict(1=>[1,3],2=>[2])

    def = SimpleProblemDef(project_spec,r0,s0,sF)
    sched, problem_spec = construct_task_graphs_problem(
        def,env_graph;cost_function=cost_function)
    return sched, problem_spec, env_graph, assignment_dict
end

"""
    pctapf_problem_3(;cost_function=SumOfMakeSpans(),verbose=false,Δt_op=0,Δt_collect=[0,0,0,0],Δt_deposit=[0,0,0,0])

In this problem robot 2 will need to yield to let robot 1 through.
First operation:
    robot 1 does [2-2-30]
Second operation:
    robot 1 does [30-30-32]
    robot 2 does [5-7-8]
Third operation:
    robot 2 does [8-12-16]
"""
function pctapf_problem_3(;cost_function=SumOfMakeSpans(),verbose=false,Δt_op=0,Δt_collect=[0,0,0,0],Δt_deposit=[0,0,0,0])
    N = 2                  # num robots
    M = 4                  # num delivery tasks
    vtx_grid = initialize_dense_vtx_grid(8,4)
    #  1   2   3   4
    #  5   6   7   8
    #  9  10  11  12
    # 13  14  15  16
    # 17  18  19  20
    # 21  22  23  24
    # 25  26  27  28
    # 29  30  31  32
    r0 = [2,5]
    s0 = [2,30,7,12]
    sF = [30,32,8,16]
    env_graph = construct_factory_env_from_vtx_grid(
        vtx_grid;
    )
    project_spec, robot_ICs = empty_pctapf_problem(r0,s0,sF)

    add_operation!(project_spec,construct_operation(project_spec,-1,[1],[2],0.0))
    add_operation!(project_spec,construct_operation(project_spec,-1,[2,3],[4],0.0))
    add_operation!(project_spec,construct_operation(project_spec,-1,[4],  [], 0.0))
    assignment_dict = Dict(1=>[1,2],2=>[3,4])

    def = SimpleProblemDef(project_spec,r0,s0,sF)
    sched, problem_spec = construct_task_graphs_problem(
        def,env_graph;cost_function=cost_function,Δt_collect=Δt_collect,Δt_deposit=Δt_deposit)

    return sched, problem_spec, env_graph, assignment_dict
end

"""
    pctapf_problem_4(;cost_function=SumOfMakeSpans(),verbose=false)

In this problem the cost of the task assignment problem is lower than the
true cost (which requires that one of the robots is delayed by a single time
step)
First operation:
    robot 1 does [2-2-8]
    robot 2 does [4-4-6]
"""
function pctapf_problem_4(;cost_function=SumOfMakeSpans(),verbose=false)
    N = 2                  # num robots
    M = 2                  # num delivery tasks
    vtx_grid = initialize_dense_vtx_grid(3,3)
    # 1  2  3
    # 4  5  6
    # 7  8  9
    r0 = [2,4]
    s0 = [2,4]
    sF = [8,6]
    env_graph = construct_factory_env_from_vtx_grid(
        vtx_grid;
    )

    project_spec, robot_ICs = empty_pctapf_problem(r0,s0,sF)
    add_operation!(project_spec,construct_operation(project_spec,-1,[1,2],[],0))
    assignment_dict = Dict(1=>[1],2=>[2])

    def = SimpleProblemDef(project_spec,r0,s0,sF)
    sched, problem_spec = construct_task_graphs_problem(
        def,env_graph;cost_function=cost_function)

    return sched, problem_spec, env_graph, assignment_dict
end

"""
    pctapf_problem_5(;cost_function=SumOfMakeSpans(),verbose=false)

In this problem the robots try to pass through each other in such a way that
an edge conflict is generated.

First operation:
    robot 1 does [3-11]
    robot 2 does [15-7]
"""
function pctapf_problem_5(;cost_function=SumOfMakeSpans(),verbose=false)
    N = 2                  # num robots
    M = 2                  # num delivery tasks
    vtx_grid = initialize_dense_vtx_grid(4,4)
    # 1   2   3   4
    # 5   6   7   8
    # 9  10  11  12
    # 13  14  15  16
    r0 = [3,15]
    s0 = [3,15]
    sF = [11,7]
    env_graph = construct_factory_env_from_vtx_grid(
        vtx_grid;
    )

    project_spec, robot_ICs = empty_pctapf_problem(r0,s0,sF)
    add_operation!(project_spec,construct_operation(project_spec,-1,[1,2],[],0))
    assignment_dict = Dict(1=>[1],2=>[2])

    def = SimpleProblemDef(project_spec,r0,s0,sF)
    sched, problem_spec = construct_task_graphs_problem(
        def,env_graph;cost_function=cost_function)

    return sched, problem_spec, env_graph, assignment_dict
end

"""
    pctapf_problem_6(;cost_function=SumOfMakeSpans(),verbose=false,Δt_op=1,Δt_collect=[0,0,0],Δt_deposit=[0,0,0])

Identical to `pctapf_problem_2`, but process time is non-zero.
In this problem robot 1 will first do [1-5-9], then [9-13-17]
robot 2 will do [4-8-32]. The key thing is that robot 1 will need to wait
until robot 2 is finished before robot 1 can do its second task
"""
function pctapf_problem_6(;cost_function=SumOfMakeSpans(),verbose=false,Δt_op=1,Δt_collect=[0,0,0],Δt_deposit=[0,0,0])
    N = 2                  # num robots
    M = 3                  # num delivery tasks
    vtx_grid = initialize_dense_vtx_grid(8,4)
    #  1   2   3   4
    #  5   6   7   8
    #  9  10  11  12
    # 13  14  15  16
    # 17  18  19  20
    # 21  22  23  24
    # 25  26  27  28
    # 29  30  31  32
    r0 = [1,4]
    s0 = [5,12,13]
    sF = [9,32,17]
    env_graph = construct_factory_env_from_vtx_grid(
        vtx_grid;
    )

    project_spec, robot_ICs = empty_pctapf_problem(r0,s0,sF)
    add_operation!(project_spec,construct_operation(project_spec,-1,[1,2],[3],Δt_op))
    add_operation!(project_spec,construct_operation(project_spec,-1,[3],  [], Δt_op))
    assignment_dict = Dict(1=>[1,3],2=>[2])

    def = SimpleProblemDef(project_spec,r0,s0,sF)
    sched, problem_spec = construct_task_graphs_problem(
        def,env_graph;cost_function=cost_function,Δt_collect=Δt_collect,Δt_deposit=Δt_deposit)

    return sched, problem_spec, env_graph, assignment_dict
end

"""
    pctapf_problem_7(;cost_function=SumOfMakeSpans(),verbose=false,Δt_op=0,Δt_collect=[0,4,0],Δt_deposit=[0,0,0])

Robot 2 will have to sit and wait at the pickup station, meaning that robot 1
will have to go around if robot 2 is on the critical path
"""
function pctapf_problem_7(;cost_function=SumOfMakeSpans(),verbose=false,Δt_op=0,Δt_collect=[0,4,0],Δt_deposit=[0,0,0])
    N = 2                  # num robots
    M = 3                  # num delivery tasks
    vtx_grid = initialize_dense_vtx_grid(4,4)
    #  1   2   3   4
    #  5   6   7   8
    #  9  10  11  12
    # 13  14  15  16
    r0 = [2,9]
    s0 = [6,10,15]
    sF = [14,12,16]
    env_graph = construct_factory_env_from_vtx_grid(
        vtx_grid;
    )

    project_spec, robot_ICs = empty_pctapf_problem(r0,s0,sF)
    add_operation!(project_spec,construct_operation(project_spec,-1,[1,2],[3],Δt_op))
    add_operation!(project_spec,construct_operation(project_spec,-1,[3],  [], Δt_op))
    assignment_dict = Dict(1=>[1,3],2=>[2])

    def = SimpleProblemDef(project_spec,r0,s0,sF)
    sched, problem_spec = construct_task_graphs_problem(
        def,env_graph;cost_function=cost_function,Δt_collect=Δt_collect,Δt_deposit=Δt_deposit)
    return sched, problem_spec, env_graph, assignment_dict
end

"""
    pctapf_problem_8(;cost_function=SumOfMakeSpans(),verbose=false,Δt_op=0,Δt_collect=[0,0,0,0],Δt_deposit=[0,0,0,0])

Two-headed project. Robot 1 does the first half of the first head, and
robot 2 handles the first half of the second head, and then they swap.
Optimal MakeSpan = 8
Optimal SumOfMakeSpans = 16
"""
function pctapf_problem_8(;cost_function=SumOfMakeSpans(),verbose=false,Δt_op=0,Δt_collect=[0,0,0,0],Δt_deposit=[0,0,0,0])
    N = 2                  # num robots
    M = 4                  # num delivery tasks
    vtx_grid = initialize_dense_vtx_grid(8,4)
    #  1   2   3   4
    #  5   6   7   8
    #  9  10  11  12
    # 13  14  15  16
    # 17  18  19  20
    # 21  22  23  24
    # 25  26  27  28
    # 29  30  31  32
    r0 = [1,29]
    s0 = [5,25,12,24]
    sF = [8,28,9,21]
    env_graph = construct_factory_env_from_vtx_grid(
        vtx_grid;
    )

    project_spec, robot_ICs = empty_pctapf_problem(r0,s0,sF)
    add_operation!(project_spec,construct_operation(project_spec,-1,[1],[4],Δt_op))
    add_operation!(project_spec,construct_operation(project_spec,-1,[2],[3],Δt_op))
    add_operation!(project_spec,construct_operation(project_spec,-1,[4],[],Δt_op))
    add_operation!(project_spec,construct_operation(project_spec,-1,[3],[],Δt_op))
    assignment_dict = Dict(1=>[1,3],2=>[2,4])

    def = SimpleProblemDef(project_spec,r0,s0,sF)
    sched, problem_spec = construct_task_graphs_problem(
        def,env_graph;cost_function=cost_function,Δt_collect=Δt_collect,Δt_deposit=Δt_deposit)

    return sched, problem_spec, env_graph, assignment_dict
end

export pctapf_problem_9

"""
    pctapf_problem_9(;cost_function=SumOfMakeSpans(),verbose=false,Δt_op=0,Δt_collect=[0,0],Δt_deposit=[0,0])

Project with station-sharing. Station 5 needs to accessed by both robots for
picking up their objects.
"""
function pctapf_problem_9(;cost_function=SumOfMakeSpans(),verbose=false,Δt_op=0,Δt_collect=[0,0],Δt_deposit=[0,0])
    N = 2                  # num robots
    M = 2                  # num delivery tasks
    vtx_grid = initialize_dense_vtx_grid(4,4)
    #  1   2   3   4
    #  5   6   7   8
    #  9  10  11  12
    # 13  14  15  16
    r0 = [1,13]
    s0 = [5,5]
    sF = [8,6]
    env_graph = construct_factory_env_from_vtx_grid(
        vtx_grid;
    )

    project_spec, robot_ICs = empty_pctapf_problem(r0,s0,sF)
    add_operation!(project_spec,construct_operation(project_spec,-1,[1],[],Δt_op))
    add_operation!(project_spec,construct_operation(project_spec,-1,[2],[],Δt_op))
    assignment_dict = Dict(1=>[1],2=>[2])

    def = SimpleProblemDef(project_spec,r0,s0,sF)
    sched, problem_spec = construct_task_graphs_problem(
        def,env_graph;cost_function=cost_function,Δt_collect=Δt_collect,Δt_deposit=Δt_deposit)

    return sched, problem_spec, env_graph, assignment_dict
end

export pctapf_problem_10

"""
    pctapf_problem_10(;cost_function=MakeSpan(),verbose=false,Δt_op=0,Δt_collect=[0,0,0,0,0,0],Δt_deposit=[0,0,0,0,0,0])

Motivation for backtracking in ISPS
The makespan optimal solution is T = 8. However, the optimistic schedule
will always prioritize task route planning for tasks 1,2, and 3 before 4.
This leads to a double delay that will not be caught without backtracking
in ISPS. Hence, the solver will return a solution with T = 9.
"""
function pctapf_problem_10(;cost_function=MakeSpan(),verbose=false,Δt_op=0,Δt_collect=[0,0,0,0,0,0],Δt_deposit=[0,0,0,0,0,0])
    N = 4                  # num robots
    M = 4                  # num delivery tasks
    vtx_grid = initialize_dense_vtx_grid(13,11)

    #   1    2    3    4    5    6    7    8    9   10   11
    #  12   13   14   15   16   17   18   19   20   21   22
    #  23   24   25   26   27   28   29   30   31   32   33
    #  34   35   36   37   38   39   40   41   42   43   44
    #  45   46   47   48   49   50   51   52   53   54   55
    #  56   57   58   59   60   61   62   63   64   65   66
    #  67   68   69   70   71   72   73   74   75   76   77
    #  78   79   80   81   82   83   84   85   86   87   88
    #  89   90   91   92   93   94   95   96   97   98   99
    # 100  101  102  103  104  105  106  107  108  109  110
    # 111  112  113  114  115  116  117  118  119  120  121
    # 122  123  124  125  126  127  128  129  130  131  132
    # 133  134  135  136  137  138  139  140  141  142  143

    #            (2)
    #
    #
    #     (1)
    # (4)    [4]     [4] [3] [4]                         [3]
    #
    #
    #
    #            [2]
    #
    #     [1]

    r0 = [15,34,137, 5       ]
    s0 = [15,34,137, 5, 27,49 ]
    sF = [22,42,60,  27,49,71]
    env_graph = construct_factory_env_from_vtx_grid(
        vtx_grid;
    )

    project_spec, robot_ICs = empty_pctapf_problem(r0,s0,sF)
    add_operation!(project_spec,construct_operation(project_spec,-1,[1,2,3,6],[],Δt_op))
    # add_operation!(project_spec,construct_operation(project_spec,-1,[2],[],Δt_op))
    add_operation!(project_spec,construct_operation(project_spec,-1,[4],[5],Δt_op))
    add_operation!(project_spec,construct_operation(project_spec,-1,[5],[6],Δt_op))
    # add_operation!(project_spec,construct_operation(project_spec,-1,[4],[],Δt_op))
    assignment_dict = Dict(1=>[1],2=>[2],3=>[3],4=>[4,5,6])

    def = SimpleProblemDef(project_spec,r0,s0,sF)
    sched, problem_spec = construct_task_graphs_problem(
        def,env_graph;cost_function=cost_function,Δt_collect=Δt_collect,Δt_deposit=Δt_deposit)

    return sched, problem_spec, env_graph, assignment_dict
end

export pctapf_problem_11

"""
    pctapf_problem_11

Requires collaborative transport: Robots 1 and 2 transport object 1 while
robot 3 transports object 2. Robot 3 will need to move over to let the other
robots pass.
"""
function pctapf_problem_11(;
        cost_function = SumOfMakeSpans(),
        verbose = false,
        Δt_op=0,
        Δt_collect=[0,0,0],
        Δt_deposit=[0,0,0],
    )
    N = 3                  # num robots
    M = 3                  # num delivery tasks
    vtx_grid = initialize_dense_vtx_grid(4, 4)
    # 1   2   3   4
    # 5   6   7   8
    # 9  10  11  12
    # 13  14  15  16
    r0 = [1, 4, 15]
    s0 = [2, 11, 13]
    sF = [14, 1, 9]
    env_graph = construct_factory_env_from_vtx_grid(
        vtx_grid;
    )
    shapes = [(1, 2), (1, 1), (1, 1)]

    project_spec, robot_ICs = empty_pctapf_problem(r0, s0, sF)
    add_operation!(
        project_spec,
        construct_operation(project_spec, -1, [1, 2], [3], Δt_op),
    )
    add_operation!(
        project_spec,
        construct_operation(project_spec, -1, [3], [], Δt_op),
    )
    assignment_dict = Dict(1 => [1,3], 2 => [1], 3 => [2])

    def = SimpleProblemDef(project_spec,r0,s0,sF,shapes)
    sched, problem_spec = construct_task_graphs_problem(
        def,env_graph;cost_function=cost_function,Δt_collect=Δt_collect,Δt_deposit=Δt_deposit)

    return sched, problem_spec, env_graph, assignment_dict
end

export pctapf_problem_12

"""
    pctapf_problem_12(;

Robot 1 will plan a path first, but then that path will need to be extended by
one time step because robot 2 will get delayed by robot 3, which is on the
critical path.
"""
function pctapf_problem_12(;
        cost_function = MakeSpan(),
        verbose = false,
    )
    vtx_grid = initialize_dense_vtx_grid(4, 8)
    #  1   2   3   4   5   6   7   8
    #  9  10  11  12  13  14  15  16
    # 17  18  19  20  21  22  23  24
    # 25  26  27  28  29  30  31  32
    env_graph = construct_factory_env_from_vtx_grid(vtx_grid)
    r0 = [1, 4, 11]

    tasks = [17=>25,4=>20,15=>16]
    shapes = map(o->(1,1),tasks)
    ops = [
        (inputs=[2],outputs=[1]),
        (inputs=[1],outputs=[3]),
        (inputs=[3],outputs=[]),
    ]
    config = (tasks=tasks,shapes=shapes,ops=ops)
    assignment_dict = Dict(1=>[1],2=>[2],3=>[3])
    def = pctapf_problem(r0,config)
    sched, problem_spec = construct_task_graphs_problem(
        def,env_graph;cost_function=cost_function)

    return sched, problem_spec, env_graph, assignment_dict
end

export pctapf_problem_13

"""
    pctapf_problem_13

Same as `pctapf_problem_12`, except that there is a 4th robot who must collect
object 1 with robot 1.
"""
function pctapf_problem_13(;
        cost_function = MakeSpan(),
        verbose = false,
    )
    vtx_grid = initialize_dense_vtx_grid(4, 8)
    #  1   2   3   4   5   6   7   8
    #  9  10  11  12  13  14  15  16
    # 17  18  19  20  21  22  23  24
    # 25  26  27  28  29  30  31  32
    env_graph = construct_factory_env_from_vtx_grid(vtx_grid)
    r0 = [1, 4, 11, 2]
    tasks = [17=>25,4=>20,15=>16]
    shapes = [(1, 2), (1, 1), (1, 1)]
    ops = [
        (inputs=[2],outputs=[1]),
        (inputs=[1],outputs=[3]),
        (inputs=[3],outputs=[]),
    ]
    config = (tasks=tasks,shapes=shapes,ops=ops)
    assignment_dict = Dict(1=>[1],2=>[2],3=>[3],4=>[1])
    def = pctapf_problem(r0,config)
    sched, problem_spec = construct_task_graphs_problem(
        def,env_graph;cost_function=cost_function)

    return sched, problem_spec, env_graph, assignment_dict
end

export
    pctapf_test_problems,
    collaborative_pctapf_test_problems

pctapf_test_problems() = [
    pctapf_problem_1,
    pctapf_problem_2,
    pctapf_problem_3,
    pctapf_problem_4,
    pctapf_problem_5,
    pctapf_problem_6,
    pctapf_problem_7,
    pctapf_problem_8,
    pctapf_problem_9,
    pctapf_problem_10,
]
collaborative_pctapf_test_problems() = [
    pctapf_problem_11,
]

for op in [
        :pctapf_problem_1,
        :pctapf_problem_2,
        :pctapf_problem_3,
        :pctapf_problem_4,
        :pctapf_problem_5,
        :pctapf_problem_6,
        :pctapf_problem_7,
        :pctapf_problem_8,
        :pctapf_problem_9,
        :pctapf_problem_10,
        :pctapf_problem_11,
        :pctapf_problem_12,
        :pctapf_problem_13,
    ]
    @eval $op(solver,args...;kwargs...) = pctapf_problem(solver,$op(args...;kwargs...)...)
end

################################################################################
############################# Replanning Problems ##############################
################################################################################

export
    replanning_problem,
    replanning_problem_1,
    replanning_problem_2,
    replanning_problem_3,
    replanning_problem_4,
    replanning_test_problems

"""
    replanning_problem

Constructs a replanning problem, consisting of robot initial conditions, an
environment, and a sequence of project requests scheduled to arrive in the
factory at regular intervals.
Args:
- r0: list of integers specifying the start locations of the robots
- defs: a list of tuples, where each tuple is of the form

    `([start_1=>goal_1, ...], [([inputs],[outputs]),...])`

    where the `start=>goal` pairs define the start and end points for each
    object to be delivered, and the `([inputs],[outputs])` pairs define the
    objects that go in and out of each operation.
- env_graph: the environment (presumably a GridFactoryEnvironment)
Outputs:
- requests: a sequence of `ProjectRequest`s
- problem_spec: a `ProblemSpec`
- robot_ICs: Robot initial conditions `ROBOT_AT`
- env_graph: the environment
"""
function replanning_problem(solver,r0,defs,env_graph;
        cost_function=SumOfMakeSpans(),
        spacing=8,
        t0=0,
        Δt_op=0,
        )
    requests = Vector{ProjectRequest}()
    problem_spec = nothing
    robot_ICs = map(i->ROBOT_AT(i,r0[i]),1:length(r0))
    for (i,def) in enumerate(defs)
        s0 = map(i->i.first,def.tasks)
        sF = map(i->i.second,def.tasks)
        spec, _ = empty_pctapf_problem(r0,s0,sF)
        for op in def.ops
            add_operation!(spec,construct_operation(spec,-1,op.inputs,op.outputs,Δt_op))
        end
        if i == 1
            problem_def = SimpleProblemDef(spec,r0,s0,sF)
            _, problem_spec = construct_task_graphs_problem(
                problem_def,env_graph;
                cost_function=cost_function)
        end
        sched = construct_partial_project_schedule(
            spec,
            problem_spec,
            )
        t = t0+spacing*(i-1)
        push!(requests,ProjectRequest(sched,t,t))
    end

    base_schedule = construct_partial_project_schedule(
        robot_ICs,
        problem_spec,
        )
    base_env = construct_search_env(solver,base_schedule,problem_spec,env_graph)
    return RepeatedPC_TAPF(base_env,requests)
end

function replanning_problem_1(solver;kwargs...)
    vtx_grid = initialize_dense_vtx_grid(4,4)
    env_graph = construct_factory_env_from_vtx_grid(vtx_grid)
    #  1   2   3   4
    #  5   6   7   8
    #  9  10  11  12
    # 13  14  15  16
    r0 = [1,13]
    defs = [
        ( tasks=[1=>4,8=>6],            ops=[ (inputs=[1],outputs=[]), (inputs=[2],outputs=[]) ] ),
        ( tasks=[13=>1,16=>15],         ops=[ (inputs=[1],outputs=[2]), (inputs=[2],outputs=[]) ] ),
        ( tasks=[9=>11,8=>11,7=>16],    ops=[ (inputs=[1,2],outputs=[3]), (inputs=[3],outputs=[]) ] ),
        ( tasks=[2=>14,3=>15,11=>4],    ops=[ (inputs=[1,2],outputs=[3]), (inputs=[3],outputs=[]) ] ),
        ]
    return replanning_problem(solver,r0,defs,env_graph;kwargs...)
end

function replanning_problem_2(solver;kwargs...)
    vtx_grid = initialize_dense_vtx_grid(4,4)
    env_graph = construct_factory_env_from_vtx_grid(vtx_grid)
    #  1   2   3   4
    #  5   6   7   8
    #  9  10  11  12
    # 13  14  15  16
    r0 = [1,13]
    defs = [
        ( tasks=[13=>1,16=>15],         ops=[ (inputs=[1],outputs=[2]), (inputs=[2],outputs=[]) ] ),
        ( tasks=[1=>4,8=>6],            ops=[ (inputs=[1],outputs=[]), (inputs=[2],outputs=[]) ] ),
        ( tasks=[2=>14,3=>15,11=>4],    ops=[ (inputs=[1,2],outputs=[3]), (inputs=[3],outputs=[]) ] ),
        ( tasks=[9=>11,8=>11,7=>16],    ops=[ (inputs=[1,2],outputs=[3]), (inputs=[3],outputs=[]) ] ),
        ]
    return replanning_problem(solver,r0,defs,env_graph;kwargs...)
end

"""
    The robot should do better if it handles the single task in the second
    project prior to working on the third task of the first project.
"""
function replanning_problem_3(solver;kwargs...)
    vtx_grid = initialize_dense_vtx_grid(4,4)
    env_graph = construct_factory_env_from_vtx_grid(vtx_grid)
    #  1   2   3   4
    #  5   6   7   8
    #  9  10  11  12
    # 13  14  15  16
    r0 = [1]
    defs = [
        ( tasks=[2=>4,8=>5,12=>9],
            ops=[
                (inputs=[1],outputs=[2]),
                (inputs=[2],outputs=[3]),
                (inputs=[3],outputs=[])
                ] ),
        ( tasks=[13=>16],
            ops=[
                (inputs=[1],outputs=[]),
                ] ),
        ( tasks=[15=>1],
            ops=[
                (inputs=[1],outputs=[]),
                ] ),
        ]
    return replanning_problem(solver,r0,defs,env_graph;kwargs...)
end

"""
    Just intended to take longer so that the tests pass even if Julia hasn't
    finished warming up yet.
"""
function replanning_problem_4(solver;kwargs...)
    vtx_grid = initialize_dense_vtx_grid(4,16)
    env_graph = construct_factory_env_from_vtx_grid(vtx_grid)
    #  1   2   3   4   5   6   7   8   9  10  11  12  13  14  15  16
    # 17  18  19  20  21  22  23  24  25  26  27  28  29  30  31  32
    # 33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48
    # 49  50  51  52  53  54  55  56  57  58  59  60  61  62  63  64
    r0 = [1]
    defs = [
        ( tasks=[1=>16,32=>17,33=>48],
            ops=[
                (inputs=[1],outputs=[2]),
                (inputs=[2],outputs=[3]),
                (inputs=[3],outputs=[])
                ] ),
        ( tasks=[64=>49],
            ops=[
                (inputs=[1],outputs=[]),
                ] ),
        ]
    return replanning_problem(solver,r0,defs,env_graph;kwargs...)
end
replanning_test_problems() = [
    replanning_problem_1,
    replanning_problem_2,
    replanning_problem_3,
    replanning_problem_4,
]

################################################################################
##################### Random RepeatedPC_TAPF instantiation #####################
################################################################################

export
    random_pctapf_def,
    random_pctapf_def,
    random_repeated_pctapf_def

"""
    instantiate_random_pctapf_def(env,config)

Instantiate a random `PC_TAPF` problem based on the parameters of config.
"""
function random_pctapf_def(env::GridFactoryEnvironment,
        config;
        N                   = get(config,:N,30),
        M                   = get(config,:M,10),
        max_parents         = get(config,:max_parents,3),
        depth_bias          = get(config,:depth_bias,0.4),
        dt_min              = get(config,:dt_min,0),
        dt_max              = get(config,:dt_max,0),
        dt_collect          = get(config,:dt_collect,0),
        dt_deliver          = get(config,:dt_deliver,0),
        task_sizes          = get(config,:task_sizes,(1=>1.0,2=>0.0,4=>0.0)),
    )

    r0,s0,sF = get_random_problem_instantiation(N,M,get_pickup_zones(env),
        get_dropoff_zones(env),get_free_zones(env))
    project_spec = construct_random_project_spec(M,s0,sF;
        max_parents=max_parents,
        depth_bias=depth_bias,
        Δt_min=dt_min,
        Δt_max=dt_max)
    shapes = choose_random_object_sizes(M,Dict(task_sizes...))
    problem_def = SimpleProblemDef(project_spec,r0,s0,sF,shapes)
end

## Stochastic PCTAPF problems
# """
#     pctapf_problem_1
#
# Optimal MakeSpan = 5
# Optimal SumOfMakeSpans = 5
# """
# function spctapf_problem_1(;cost_function=SumOfMakeSpans(),verbose=false)
#     vtx_grid = initialize_dense_vtx_grid(4,4)
#     env_graph = construct_factory_env_from_vtx_grid(vtx_grid)
#     #  1   2   3   4
#     #  5   6   7   8
#     #  9  10  11  12
#     # 13  14  15  16
#     config = (
#         delivery_bots = [1],
#         cleanup_bots = [13],
#         tasks=[1=>4,8=>16],
#         ops=[ (inputs=[1],outputs=[]), (inputs=[2],outputs=[]) ] )
#     )
#     project_spec, robot_ICs = empty_pctapf_problem(r0,s0,sF)
#
#     add_operation!(project_spec,construct_operation(project_spec,-1,[1,2],[3],0))
#     add_operation!(project_spec,construct_operation(project_spec,-1,[3],  [], 0))
#     assignment_dict = Dict(1=>[1,3],2=>[2])
#
#     def = SimpleProblemDef(project_spec,r0,s0,sF)
#     sched, problem_spec = construct_task_graphs_problem(
#         def,env_graph;cost_function=cost_function)
#     return sched, problem_spec, env_graph, assignment_dict
# end


text:
module NAGD
export NAGDContext, NesterovGradientDescent, fit!
using ..Classification: OptimizerContext, Optimizer, σ, logcrossentropy, sigmoid!, logcrossentropy!, stopearly!, postprocess!
import ..Classification: fit!
using LinearAlgebra: mul!
#using Random: shuffle!
using Statistics: mean!

struct NAGDContext <: OptimizerContext
    learningrate::Float64
    decayrate::Float64
    tolerance::Float64
    maxiterations::Int
    batchsize::Int
end
NAGDContext(;learningrate=0.1, decayrate=0.8, tolerance=eps(Float64),
           maxiterations=1000, batchsize=8) = NAGDContext(learningrate, decayrate, tolerance, maxiterations, batchsize)

mutable struct NesterovGradientDescent <: Optimizer
    context::NAGDContext
    converged::Bool
    earlystopping::Bool
    hasvalidationset::Bool
    iterations::Int
    loss::Vector{Float64}
    validationloss::Vector{Float64}
    validationset::Tuple{Matrix{Float64}, Vector{Float64}}
end

function NesterovGradientDescent(context=NAGDContext())
    NesterovGradientDescent(context, false, false, false, 0, Float64[], Float64[],
                            (zeros(Float64, (1, 1)), Float64[]))
end

function fit!(optim::NesterovGradientDescent, β, X, y)
    context = optim.context
    previousloss = -Inf
    optim.loss = zeros(Float64, context.maxiterations)
    momentum = zeros(Float64, size(β))
    n, k = size(X)
    m = context.batchsize
    βnear     = zeros(Float64, length(β))
    ŷbatch    = zeros(Float64, m)
    ŷ         = zeros(Float64, n)
    residuals = zeros(Float64, m)
    ∇array    = zeros(Float64, (1, length(β)))
    ∇         = @view ∇array[1, :]
    tmp1batch = zeros(Float64, m)
    tmp1      = zeros(Float64, n)
    tmp2      = zeros(Float64, n)
    tmp       = zeros(Float64, (m, k))
    loss      = [0.0]
    Xbatch    = zeros(Float64, (m, k))
    ybatch    = zeros(Float64, m)
    batches   = partition(n, context.batchsize)
    batch     = 1:1

    @inbounds for i in 1:context.maxiterations
        optim.iterations += 1
        shuffle!(X, y)

        for batch in batches
            move!(X, y, Xbatch, ybatch, batch)

            @. βnear = β + context.decayrate * momentum
            mul!(tmp1batch, Xbatch, βnear)
            sigmoid!(ŷbatch, tmp1batch)

            # Gradient descent
            residuals .= ŷbatch .- ybatch
            tmp       .= Xbatch.*residuals
            mean!(∇array, tmp)

            @. momentum = context.decayrate*momentum - context.learningrate * ∇
            β         .+= momentum
        end

        mul!(tmp1, X, β)
        sigmoid!(ŷ, tmp1)

        logcrossentropy!(loss, tmp2, y, ŷ)
        optim.loss[i] = loss[1]
        if abs(previousloss - loss[1]) < context.tolerance
            optim.converged = true
            break
        else
            previousloss = loss[1]
        end

        if stopearly!(optim, β, i)
            break
        end
    end
    postprocess!(optim)
end

function partition(length::Number, batchsize::Number)::Vector{UnitRange{Int}}
    slices = Vector{UnitRange{Int}}(undef, length ÷ batchsize)
    lastaddon = length % batchsize
    start = 1
    i = 1
    while true
        stop = start + batchsize - 1
        if stop > length
            break
        end
        slices[i] = start:stop
        start = stop + 1
        i += 1
    end
    #slices[end] = first(slices[end]):(last(slices[end])+lastaddon)
    # I don't know how to do the last batch
    pop!(slices)
    slices
end

function shuffle!(X::Matrix{T} where T, y::Vector{V} where V)
    n = length(y)
    #@assert n == size(X, 1)
    #@assert n > 1
    j = 0
    @inbounds for i in n:-1:2
        @fastmath j = 1 + abs(rand(Int))%n
        swaprows!(X, i, j)
        y[i], y[j] = y[j], y[i]
    end
end

function swaprows!(X, i, j)
    for k in 1:size(X, 2)
        X[i, k], X[j, k] = X[j, k], X[i, k]
    end
end

function move!(X, y, X̃, ỹ, indices)
    j = 1
    for i in indices
        moverow!(X, X̃, i, j)
        ỹ[j] = y[i]
        j += 1
    end
end

function moverow!(X, X̃, i, j)
    for k in 1:size(X, 2)
        X̃[j, k] = X[i, k]
    end
end

end

ext: jl


text:
#
# Date created: 2022-03-14
# Author: aradclif
#
#
############################################################################################
function sizeblock(N::Int)
    block = Expr(:block)
    for d = 1:N
        ex = Expr(:(=), Symbol(:D_, d), Expr(:call, :size, :A, d))
        push!(block.args, ex)
    end
    block
end

function sizeproductsblock(N::Int)
    block = Expr(:block)
    for k = 3:N
        ex = Expr(:(=), Symbol(:D_, ntuple(identity, k - 1)...),
                  Expr(:call, :*, ntuple(d -> Symbol(:D_, d), k - 1)...))
        push!(block.args, ex)
    end
    block
end

function sumprodprecomputed2(N::Int)
    e = Expr(:call, :+, Symbol(:i_, 1))
    for k = 2:N
        if k == 2
            ex = Expr(:call, :*, Symbol(:D_, 1), Symbol(:i_, 2))
            push!(e.args, ex)
        else
            ex = Expr(:call, :*, Symbol(:D_, ntuple(identity, k - 1)...), Symbol(:i_, k))
            push!(e.args, ex)
        end
    end
    e
end
function sumprodprecomputed2(N::Int, D)
    params = D.parameters
    e = params[1] !== Static.One ? Expr(:call, :+, Symbol(:i_, 1)) : Expr(:call, :+, Symbol(:j_, 1))
    for k = 2:N
        if k == 2
            sym = params[2] !== Static.One ? :i_ : :j_
            ex = Expr(:call, :*, Symbol(:D_, 1), Symbol(sym, 2))
            push!(e.args, ex)
        else
            sym = params[k] !== Static.One ? :i_ : :j_
            ex = Expr(:call, :*, Symbol(:D_, ntuple(identity, k - 1)...), Symbol(sym, k))
            push!(e.args, ex)
        end
    end
    e
end

function sumprodconstant(N::Int)
    Expr(:(=), :D_sp,
         Expr(:call, :+, ntuple(d -> Expr(:call, :*, ntuple(i -> Symbol(:D_, i), d)...), N - 1)...))
end

function outerloopgen(N::Int, D)
    # loops = Expr(:for)
    block = Expr(:block)
    params = D.parameters
    for d = N:-1:1
        if params[d] !== Static.One
            ex = Expr(:(=), Symbol(:i_, d), Expr(:call, :axes, :A, d))
            push!(block.args, ex)
        end
    end
    # push!(loops.args, block)
    Expr(:for, block)
end
function innerloopgen(N::Int, D)
    # loops = Expr(:for)
    block = Expr(:block)
    params = D.parameters
    for d = N:-1:1
        if params[d] === Static.One
            ex = Expr(:(=), Symbol(:i_, d), Expr(:call, :axes, :A, d))
            push!(block.args, ex)
        end
    end
    # push!(loops.args, block)
    # loops
    Expr(:for, block)
end

function innerpost(N::Int, D)
    params = D.parameters
    # block = Expr(:block)
    b = Expr(:ref, :B, ntuple(d -> params[d] === Static.One ? 1 : Symbol(:i_, d), N)...)
    c = Expr(:ref, :C, ntuple(d -> params[d] === Static.One ? 1 : Symbol(:i_, d), N)...)
    # e1 = Expr(:(=), b, :m)
    # e2 = Expr(:(=), c, :j)
    # push!(block.args, e1)
    # push!(block.args, e2)
    # block
    Expr(:block, Expr(:(=), b, :m), Expr(:(=), c, :j))
end
function innerpost2(N::Int, D)
    params = D.parameters
    block = Expr(:block)
    b = Expr(:ref, :B, (Symbol(:i_, d) for d = 1:N if params[d] !== Static.One)...)
    c = Expr(:ref, :C, (Symbol(:i_, d) for d = 1:N if params[d] !== Static.One)...)
    e1 = Expr(:(=), b, :m)
    e2 = Expr(:(=), c, :j)
    push!(block.args, e1)
    push!(block.args, e2)
    block
end

function innerpost2_leading1(N::Int, D)
    params = D.parameters
    block = Expr(:block)
    b = Expr(:ref, :B, ntuple(d -> params[d] === Static.One ? 1 : Symbol(:i_, d), N)[2:N]...)
    c = Expr(:ref, :C, ntuple(d -> params[d] === Static.One ? 1 : Symbol(:i_, d), N)[2:N]...)
    e1 = Expr(:(=), b, :m)
    e2 = Expr(:(=), c, :j)
    push!(block.args, e1)
    push!(block.args, e2)
    block
end

function maxcompareblock(N::Int)
    a = Expr(:ref, :A, ntuple(d -> Symbol(:i_, d), N)...)
    d = sumprodprecomputed2(N)
    push!(d.args, :D_sp)
    yₑ = Expr(:(=), :y, Expr(:call, :(>), a, :m)) # f should only be > or <
    mₑ = Expr(:(=), :m, Expr(:if, :y, a, :m))
    jₑ = Expr(:(=), :j, Expr(:if, :y, d, :j))
    Expr(:block, yₑ, mₑ, jₑ)
end

function maxcompareblock2(N::Int)
    a = Expr(:ref, :A, ntuple(d -> Symbol(:i_, d), N)...)
    d = sumprodprecomputed2(N)
    push!(d.args, :D_sp)
    yₑ = Expr(:(=), :y, Expr(:call, :(>), a, :m)) # f should only be > or <
    mₑ = Expr(:(=), :m, Expr(:call, :ifelse, :y, a, :m))
    jₑ = Expr(:(=), :j, Expr(:call, :ifelse, :y, d, :j))
    # jₑ = Expr(:(=), :j, Expr(:call, :+, Expr(:call, :zero, :Tₒ), Expr(:call, :ifelse, :y, d, :j)))
    Expr(:block, yₑ, mₑ, jₑ)
end

# function maxcompareblock3(N::Int, D)
#     block = Expr(:block)
#     params = D.parameters
#     a = Expr(:ref, :A, ntuple(d -> Symbol(:i_, d), N)...)
#     # b = Expr(:ref, :B, ntuple(d -> params[d] === Static.One ? 1 : Symbol(:i_, d), N)...)
#     d = sumprodprecomputed2(N)
#     push!(d.args, :D_sp)
#     jₑ = Expr(:(=), :j, Expr(:call, :ifelse, Expr(:call, :(>=), a, :m), d, :j))
#     push!(block.args, jₑ)
#     block
# end

function maxcompareblock4(N::Int)
    a = Expr(:ref, :A, ntuple(d -> Symbol(:i_, d), N)...)
    d = sumprodprecomputed2(N)
    push!(d.args, :D_sp)
    jₑ = Expr(:(+=), :j, Expr(:call, :ifelse, Expr(:call, :(>), a, :m), d, 0))
    mₑ = Expr(:(=), :m, Expr(:call, :max, a, :m))
    Expr(:block, jₑ, mₑ)
end
function lvfindmax_quote4(N::Int, D)
    b1 = sizeblock(N)
    b2 = sizeproductsblock(N)
    b3 = Expr(:block, sumprodconstant(N), Expr(:(=), :D_sp, Expr(:call, :-, :D_sp)))
    outer = outerloopgen(N, D)
    b4 = Expr(:block, Expr(:(=), :j, Expr(:call, :zero, :Tₒ)),
              Expr(:(=), :m, Expr(:call, :typemin, :T)))
    inner = innerloopgen(N, D)
    b5 = maxcompareblock4(N)
    push!(inner.args, b5)
    push!(b4.args, inner)
    b6 = innerpost(N, D)
    push!(b4.args, b6.args...)
    push!(outer.args, b4)
    return quote
        $b1
        $b2
        $b3
        @turbo $outer
    end
end

# function maxcompareblock5(N::Int, D)
#     a = Expr(:ref, :A, ntuple(d -> Symbol(:i_, d), N)...)
#     yₑ = Expr(:(=), :y, Expr(:call, :(>), a, :m)) # f should only be > or <
#     mₑ = Expr(:(=), :m, Expr(:call, :ifelse, :y, a, :m))
#     # jₑ = Expr(:(=), Symbol(:j_, d), Expr(:call, :ifelse, Expr(:y), Symbol(:i_, d), Symbol(:j_, d)))
#     js = Any[]
#     for d = 1:N
#         if params[d] === Static.One
#             ex = Expr(:(=), Symbol(:j_, d),
#                       Expr(:call, :ifelse, Expr(:y), Symbol(:i_, d), Symbol(:j_, d)))
#             push!(js, ex)
#         end
#     end
#     Expr(:block, yₑ, mₑ, js...)
# end

# function innerpost5(N::Int, D)
#     params = D.parameters
#     d = sumprodprecomputed2(N, D)
#     push!(d.args, :D_sp)
#     b = Expr(:ref, :B, ntuple(d -> params[d] === Static.One ? 1 : Symbol(:i_, d), N)...)
#     c = Expr(:ref, :C, ntuple(d -> params[d] === Static.One ? 1 : Symbol(:i_, d), N)...)
#     Expr(:block, Expr(:(=), b, :m), Expr(:(=), c, d))
# end

function lvfindmax_quote(N::Int, D)
    b1 = sizeblock(N)
    b2 = sizeproductsblock(N)
    b3 = Expr(:block, sumprodconstant(N), Expr(:(=), :D_sp, Expr(:call, :-, :D_sp)))
    outer = outerloopgen(N, D)
    b4 = Expr(:block, Expr(:(=), :j, 1), Expr(:(=), :m, Expr(:call, :typemin, :T)))
    inner = innerloopgen(N, D)
    b5 = maxcompareblock(N)
    push!(inner.args, b5)
    push!(b4.args, inner)
    b6 = innerpost(N, D)
    push!(b4.args, b6.args...)
    push!(outer.args, b4)
    return quote
        $b1
        $b2
        $b3
        @turbo $outer
    end
end
@generated function _lvfindmax!(C::AbstractArray{Tₒ, N}, A::AbstractArray{T, N},
                             B::AbstractArray{T, N}, dims::D) where {Tₒ, T, N, D}
    lvfindmax_quote(N, D)
end
function lvfindmax(A::AbstractArray{T, N}, dims::NTuple{M, Int}) where {T, N, M}
    Dᴮ′ = ntuple(d -> d ∈ dims ? StaticInt(1) : size(A, d), N)
    B = similar(A, Dᴮ′)
    C = similar(B, Int)
    _lvfindmax!(C, A, B, Dᴮ′)
    B, CartesianIndices(A)[C]
end
function lvfindmax_quote2(N::Int, D)
    b1 = sizeblock(N)
    b2 = sizeproductsblock(N)
    b3 = Expr(:block, sumprodconstant(N), Expr(:(=), :D_sp, Expr(:call, :-, :D_sp)))
    outer = outerloopgen(N, D)
    # b4 = Expr(:block, Expr(:(=), :j, Expr(:call, :zero, :Tₒ)),
    #           Expr(:(=), :m, Expr(:call, :typemin, :T)))
    b4 = Expr(:block, Expr(:(=), :j, Expr(:call, :zero, :Tₒ)),
              Expr(:(=), :m, Expr(:call, :typemin, :T)))
    inner = innerloopgen(N, D)
    b5 = maxcompareblock2(N)
    push!(inner.args, b5)
    push!(b4.args, inner)
    b6 = innerpost(N, D)
    push!(b4.args, b6.args...)
    push!(outer.args, b4)
    return quote
        $b1
        $b2
        $b3
        @turbo $outer
    end
end
@generated function _lvfindmax2!(C::AbstractArray{Tₒ, N}, A::AbstractArray{T, N},
                             B::AbstractArray{T, N}, dims::D) where {Tₒ, T, N, D}
    lvfindmax_quote2(N, D)
end
function lvfindmax2(A::AbstractArray{T, N}, dims::NTuple{M, Int}) where {T, N, M}
    Dᴮ′ = ntuple(d -> d ∈ dims ? StaticInt(1) : size(A, d), N)
    B = similar(A, Dᴮ′)
    C = similar(B, Int)
    _lvfindmax2!(C, A, B, Dᴮ′)
    B, CartesianIndices(A)[C]
end

# function lvfindmax_quote3(N::Int, D)
#     b1 = sizeblock(N)
#     b2 = sizeproductsblock(N)
#     b3 = Expr(:block, sumprodconstant(N), Expr(:(=), :D_sp, Expr(:call, :-, :D_sp)))
#     outer = outerloopgen(N, D)
#     b4 = Expr(:block, Expr(:(=), :j, 1), Expr(:(=), :m, Expr(:call, :typemin, :T)))
#     inner = innerloopgen(N, D)
#     b5 = maxcompareblock2(N)
#     push!(inner.args, b5)
#     push!(b4.args, inner)
#     b6 = D.parameters[1] === StaticInt{1} ? innerpost2_leading1(N, D) : innerpost(N, D)
#     push!(b4.args, b6.args...)
#     push!(outer.args, b4)
#     return quote
#         $b1
#         $b2
#         $b3
#         @turbo $outer
#     end
# end
# @generated function _lvfindmax3!(C::AbstractArray{Tₒ, M}, A::AbstractArray{T, N},
#                                  B::AbstractArray{T, M}, dims::D) where {Tₒ, M, T, N, D}
#     lvfindmax_quote3(N, D)
# end
# function lvfindmax3(A::AbstractArray{T, N}, dims::NTuple{M, Int}) where {T, N, M}
#     if 1 ∈ dims
#         Dᴮ′ = ntuple(d -> d ∈ dims ? StaticInt(1) : size(A, d), N)
#         Dᴮ = ntuple(d -> d + 1 ∈ dims ? StaticInt(1) : size(A, d), N - 1)
#         B = similar(A, Dᴮ)
#         C = similar(B, Int)
#         _lvfindmax3!(C, A, B, Dᴮ′)
#         return reshape(B, Dᴮ′), CartesianIndices(A)[reshape(C, Dᴮ′)]
#     else
#         Dᴮ′ = ntuple(d -> d ∈ dims ? StaticInt(1) : size(A, d), N)
#         B = similar(A, Dᴮ′)
#         C = similar(B, Int)
#         _lvfindmax2!(C, A, B, Dᴮ′)
#         return B, CartesianIndices(A)[C]
#     end
# end

function lvfindmax_quote5(N::Int, D)
    b1 = sizeblock(N)
    b2 = sizeproductsblock(N)
    b3 = Expr(:block, sumprodconstant(N), Expr(:(=), :D_sp, Expr(:call, :-, :D_sp)))
    outer = outerloopgen(N, D)
    params = D.parameters
    js = Any[]
    for d = 1:N
        if params[d] === Static.One
            ex = Expr(:(=), Symbol(:j_, d), Expr(:call, :zero, :Tₒ))
            push!(js, ex)
        end
    end
    b4 = Expr(:block, js..., Expr(:(=), :m, Expr(:call, :typemin, :T)))
    inner = innerloopgen(N, D)
    b5 = maxcompareblock5(N, D)
    push!(inner.args, b5)
    push!(b4.args, inner)
    b6 = innerpost5(N, D)
    push!(b4.args, b6.args...)
    push!(outer.args, b4)
    return quote
        $b1
        $b2
        $b3
        @turbo $outer
    end
end

function maxcompareblock5(N::Int, D)
    a = Expr(:ref, :A, ntuple(d -> Symbol(:i_, d), N)...)
    yₑ = Expr(:(=), :y, Expr(:call, :(>), a, :m)) # f should only be > or <
    mₑ = Expr(:(=), :m, Expr(:call, :ifelse, :y, a, :m))
    # jₑ = Expr(:(=), Symbol(:j_, d), Expr(:call, :ifelse, Expr(:y), Symbol(:i_, d), Symbol(:j_, d)))
    params = D.parameters
    js = Any[]
    for d = 1:N
        if params[d] === Static.One
            ex = Expr(:(=), Symbol(:j_, d),
                      Expr(:call, :ifelse, :y, Symbol(:i_, d), Symbol(:j_, d)))
            push!(js, ex)
        end
    end
    Expr(:block, yₑ, mₑ, js...)
end

function innerpost5(N::Int, D)
    params = D.parameters
    d = sumprodprecomputed2(N, D)
    push!(d.args, :D_sp)
    b = Expr(:ref, :B, ntuple(d -> params[d] === Static.One ? 1 : Symbol(:i_, d), N)...)
    c = Expr(:ref, :C, ntuple(d -> params[d] === Static.One ? 1 : Symbol(:i_, d), N)...)
    Expr(:block, Expr(:(=), b, :m), Expr(:(=), c, d))
end

@generated function _lvfindmax5!(C::AbstractArray{Tₒ, N}, A::AbstractArray{T, N},
                                 B::AbstractArray{T, N}, dims::D) where {Tₒ, T, N, D}
    lvfindmax_quote5(N, D)
end
function lvfindmax5(A::AbstractArray{T, N}, dims::NTuple{M, Int}) where {T, N, M}
    Dᴮ′ = ntuple(d -> d ∈ dims ? StaticInt(1) : size(A, d), N)
    B = similar(A, Dᴮ′)
    C = similar(B, Int)
    _lvfindmax5!(C, A, B, Dᴮ′)
    B, CartesianIndices(A)[C]
end

A = rand(4, 3, 5);

@benchmark findmax(A, dims=dims)

@benchmark bfindmax(A, dims)

@benchmark lvfindmax(A, dims)
@benchmark lvfindmax2(A, dims)
@benchmark lvfindmax3(A, dims)
@benchmark lvfindmax5(A, dims)

N = ndims(A)
DD = typeof(ntuple(d -> d ∈ dims ? StaticInt(1) : size(A, d), N))
lvfindmax_quote(N, DD)
lvfindmax_quote2(N, DD)
lvfindmax_quote3(N, DD)
lvfindmax_quote4(N, DD)
lvfindmax_quote5(N, DD)
lvfindmax2(A, dims) == findmax(A, dims=dims)
lvfindmax5(A, dims) == findmax(A, dims=dims)
lvfindmax(A, (1,))
findmax(A, dims=dims)

dims=(2,)
dims=(1,)
dims=(2,3)
dims=(2,3,4,5)

for d₂ = 2:ndims(A), d₁ = 2:ndims(A)
    dims = (d₁, d₂)
    @assert lvfindmax2(A, dims) == findmax(A, dims=dims)
    @assert lvfindmax3(A, dims) == findmax(A, dims=dims)
end

@benchmark bfindmax(A, dims)
@benchmark findmax(A, dims=dims)
A = rand(10, 100, 1000);
dims = (1,)


text:
module HashDicts

export HashDict

include("bucket.jl")
include("hashdict.jl")

end # module


text:
# module MySQL_Octo.SQL

function iterate(query::QueryBuilders.Queryable)
    MySQLRowIterator(query.adapter.handle, repr(query))
end


text:
import LinearAlgebra
import FileIO
import GeometryBasics
import Observables
import Makie
import GLMakie

using Porta

startframe = 1
FPS = 60
resolution = (3840, 2160)
segments1 = 48
segments2 = 9
speed = 1
factor = 0.0025
scale = 1.0
radius = 0.002
number = 720
frames = 3600
basemapradius = 0.325
basemapsegments = 2segments1
basepointradius = 0.004
basepointtransparency = false
exportmode = ["gif", "frames", "video"]
exportmode = exportmode[2]
modelname = "twospinorsalpha3"

# The scene object that contains other visual objects
scene = Makie.Scene(backgroundcolor = :black,
                               show_axis = false,
                               resolution = resolution)

quaternionicrgb(p::Quaternion, q::Quaternion, s::Real, v::Real) = begin
    hue = acos(dot(p, q)) / π
    hsvtorgb([hue * 360; s; v])
end

"""
    getflower(;N=4, A=.5, B=-pi/7, P=pi/2, Q=0, number=300)

Calculates the x, y and z points of a flower in the base space.
with the given number of petals N, the fattness of the petals A,
the height of the petals B, the latitude of the flower P,
the rotation of the flower Q, and the total number of points in the grid.
"""
function getflower(;N=4, A=.5, B=-pi/7, P=pi/2, Q=0, number=300)
    N = 6
    A = .5
    B = -pi/6
    P = π / 2 - π / 6
    Q = -π / 3
    t = range(0, stop = 2pi, length = number)
    az = 2pi .* t + A .* cos.(N .* 2pi .* t) .+ Q
    po = B .* sin.(N .* 2pi .* t) .+ P
    x = cos.(az).*sin.(po)
    y = sin.(az).*sin.(po)
    z = cos.(po)
    points = Array{Geographic}(undef, number)
    for i in 1:number
        points[i] = Geographic(Cartesian(x[i], y[i], z[i]))
    end
    points
end

A(s, p::S²) = begin
    z = s(p)
    z₀, z₁ = vec(ComplexPlane(z))
    r = vec(z) # radial
    f = [-r[2]; r[1]; -r[4]; r[3]] # fiber tangent
    X₀, X₁ = vec(ComplexPlane(Quaternion(f)))
    1/2 * (conj(z₀) * X₀ - z₀ * conj(X₀) + conj(z₁) * X₁ - z₁ * conj(X₁))
end

A(p::ComplexPlane) = begin
    z₀, z₁ = vec(p)
    r = vec(Quaternion(p)) # radial
    f = [-r[2]; r[1]; -r[4]; r[3]] # fiber tangent
    X₀, X₁ = vec(ComplexPlane(Quaternion(f)))
    1/2 * (conj(z₀) * X₀ - z₀ * conj(X₀) + conj(z₁) * X₁ - z₁ * conj(X₁))
end

basemapconfig = Biquaternion(Quaternion(1, 0, 0, 0), ℝ³(-0.9 + basemapradius / 2, -√2 - basemapradius / 2, 0))
basemapcolor = Makie.RGBAf(0.25, 0.25, 0.25, 0.25)
basemaptransparency = true
basemap = Sphere(basemapconfig,
                 scene,
                 radius = basemapradius,
                 segments = basemapsegments,
                 color = basemapcolor,
                 transparency = basemaptransparency)

curvepoints = getflower(number = number)
solidcolors = Array{Makie.RGBAf,1}(undef, number)
ghostcolors = similar(solidcolors)
basepoints = Array{Sphere,1}(undef, number)
for i in 1:number
    config = Biquaternion(ℝ³(Cartesian(curvepoints[i])) * basemapradius)
    q₁ = λ⁻¹map(curvepoints[i])
    q₂ = Quaternion(1, 0, 0, 0)
    rgb = quaternionicrgb(q₁, q₂, 1, 1)
    solidcolors[i] = Makie.RGBAf(rgb..., 0.7)
    ghostcolors[i] = Makie.RGBAf(rgb..., 0.1)
    sphere = Sphere(basemapconfig * config,
                    scene,
                    radius = basepointradius,
                    segments = segments2,
                    color = solidcolors[i],
                    transparency = basepointtransparency)
    basepoints[i] = sphere
end

bundleconfig = Biquaternion(ℝ³(0, 0, 0))
solidspinors = Array{Twospinor,1}(undef, number)
ghostspinors = Array{Twospinor,1}(undef, number)
solidgauge1 = U1(0)
solidgauge2 = U1(2π)
ghostgauge1 = U1(0)
ghostgauge2 = U1(2π)
spinors = Array{ComplexPlane,1}(undef, number)
for i in 1:number
    transparency = false
    point = ComplexPlane(λ⁻¹map(curvepoints[i]))
    spinors[i] = point
    solidspinor = Twospinor(scene,
                            point,
                            gauge1 = solidgauge1,
                            gauge2 = solidgauge2,
                            configuration = bundleconfig,
                            radius = radius,
                            segments1 = segments1,
                            segments2 = segments2,
                            color = solidcolors[i],
                            transparency = transparency)
    solidspinors[i] = solidspinor
    transparency = true
    ghostspinor = Twospinor(scene,
                            point,
                            gauge1 = ghostgauge1,
                            gauge2 = ghostgauge2,
                            configuration = bundleconfig,
                            radius = radius,
                            segments1 = segments1,
                            segments2 = segments2,
                            color = ghostcolors[i],
                            transparency = transparency)
    ghostspinors[i] = ghostspinor
end

basespacepoints = Array{ComplexLine,1}(undef, number)

width = 2basepointradius
tail, head = ℝ³(rand(3)), ℝ³(rand(3))
arrows = []
for i in 1:number
    push!(arrows, Arrow(tail, head, scene, width = width, color = solidcolors[i], transparency = false))
end

tori = []
for i in 1:number
    push!(tori, Torus(Biquaternion(tail), scene, r = basepointradius, R = 0.15,
                      segments = 30, color = solidcolors[i], transparency = true))
end

group2algebra(A::SU2, r::Float64) = begin
    τ₃ = -0.5 .* [im 0; 0 -im]
    (-2 * r) .* (A.a * τ₃ * LinearAlgebra.inv(A.a))
end

exponentiate(M::Array{<:Complex,2}, t::Float64) = begin
    eigenvalues = LinearAlgebra.eigvals(M)
    eigenvectors = LinearAlgebra.eigvecs(M)
    S = eigenvectors
    S⁻¹ = LinearAlgebra.inv(S)
    e = LinearAlgebra.Diagonal(eigenvalues) .* t
    S = convert.(Complex{BigFloat}, S)
    S⁻¹ = convert.(Complex{BigFloat}, S⁻¹)
    e = convert.(Complex{BigFloat}, e)
    G = S * ℯ.^e * S⁻¹
    z₁, z₂ = G[1,1], G[1,2]
    v = LinearAlgebra.normalize([real(z₁); imag(z₁); real(z₂); imag(z₂)])
    z₁, z₂ = Complex{Float64}(v[1] + im * v[2]), Complex{Float64}(v[3] + im * v[4])
    G = [z₁ -conj(z₂); z₂ conj(z₁)]
    SU2(G)
end

L₁ = [Complex(0) Complex(1); Complex(1) Complex(0)]
L₂ = [Complex(0) Complex(-im); Complex(im) Complex(0)]
L₃ = [Complex(1) Complex(0); Complex(0) Complex(-1)]
rotation = 4π
# operator = exponentiate(L₃, 0.0)
ϵ = float(π / 2)

"""
    animate(i)

Update the state of observables with the given frame number `i`.
"""
function animate(i::Int)
    step = i / frames
    τ = step * rotation
    ψ = Int(τ ÷ ϵ)
    θ = τ % ϵ

    for index in 1:number
        solidspinor = solidspinors[index]
        ghostspinor = ghostspinors[index]
        p = spinors[index]
        operator = exponentiate(group2algebra(SU2(p), 0.0), 0.0)
        q₁ = [exponentiate(group2algebra(operator, 1.0), ϵ) for x in 1:ψ]
        q = exponentiate(group2algebra(operator, 1.0), θ)
        for item in q₁
            q = SU2(normalize(Quaternion(item * q)))
        end
        p = ComplexPlane(normalize(Quaternion(SU2(p) * q)))
        α = A(p)
        g₁ = U1(0)
        g₂ = U1(α)
        g₃ = U1(2π)
        update(solidspinor, p, g₁, g₂)
        update(ghostspinor, p, g₂, g₃)
        rgb = quaternionicrgb(Quaternion(p), Quaternion(1, 0, 0, 0), 0.75 + 0.25 * rand(), 0.75 + 0.25 * rand())
        rgb = 0.8 .* rgb + 0.2 .* rand(3)
        solidcolor = Makie.RGBAf(rgb..., 0.5 + 0.2 * rand())
        ghostcolor = Makie.RGBAf(rgb..., 0.2 + 0.2 * rand())
        update(solidspinor, solidcolor)
        update(ghostspinor, ghostcolor)

        basepointmarker = basepoints[index]
        p = πmap(Quaternion(p))
        p = ℝ³(Cartesian(p)) * basemapradius
        configuration = Biquaternion(p)
        update(basepointmarker, basemapconfig * configuration)
        update(basepointmarker, solidcolor)

        arrow = arrows[index]
        initial = ℝ³(0, 0, basemapradius)
        q₁ = Biquaternion(getrotation(normalize(initial), normalize(p)))
        q₂ = Biquaternion(Quaternion(angle(α), p))
        q₃ = Biquaternion(initial + ℝ³(0.15, 0, 0))
        tail = p
        head = gettranslation(q₁ * q₂ * q₃)
        tail, head = applyconfig([tail; head], basemapconfig)
        update(arrow, tail, head - tail)
        update(arrow, solidcolor)

        torus = tori[index]
        q₁ = Biquaternion(p)
        q₂ = Biquaternion(getrotation(normalize(p), normalize(initial)))
        update(torus, basemapconfig * q₁ * q₂)
        update(torus, ghostcolor)
    end

    # n = ℝ³(1, 0, 0)
    # v = ℝ³(0, 0, 2) * √2
    # q = Quaternion(τ / 4, ℝ³(0, 1, 0))
    # v = rotate(v, q)
    # # update eye position
    # # scene.camera.eyeposition.val
    # upvector = GeometryBasics.Vec3f0(vec(n)...)
    # eyeposition = GeometryBasics.Vec3f0(vec(v)...)
    # lookat = GeometryBasics.Vec3f0(0, 0, 0)
    # Makie.update_cam!(scene, eyeposition, lookat, upvector)
    # scene.center = false # prevent scene from recentering on display
end

n = ℝ³(1, 0, 0)
v = ℝ³(0, 0, 2) * √2
# update eye position
# scene.camera.eyeposition.val
upvector = GeometryBasics.Vec3f0(vec(n)...)
eyeposition = GeometryBasics.Vec3f0(vec(v)...)
lookat = GeometryBasics.Vec3f0(0, 0, 0)
Makie.update_cam!(scene, eyeposition, lookat, upvector)
scene.center = false # prevent scene from recentering on display
if exportmode ∈ ["gif", "video"]
    outputextension = exportmode == "gif" ? "gif" : "mkv"
    Makie.record(scene, "gallery/$modelname.$outputextension",
                            framerate = FPS) do io
        for i in startframe:frames
            animate(i) # animate the scene
            Makie.recordframe!(io) # record a new frame
            step = (i - 1) / frames
            println("Completed step $(100step).\n")
        end
    end
elseif exportmode == "frames"
    directory = joinpath("gallery", modelname)
    !isdir(directory) && mkdir(directory)
    directory = joinpath("gallery", modelname)
    !isdir(directory) && mkdir(directory)
    for i in startframe:frames
        start = time()
        animate(i) # animate the scene
        elapsed = time() - start
        sleep(0.1)
        println("Generating frame $i took $elapsed (s).")

        start = time()

        paddingnumber = length(digits(frames))
        imageid = lpad(i, paddingnumber, "0")
        imagename = "$(modelname)_$(resolution[2])p_$(FPS)fps_$imageid.jpeg"
        filepath = joinpath(directory, imagename)
        FileIO.save(filepath,
                    scene;
                    resolution = resolution,
                    pt_per_unit = 400.0,
                    px_per_unit = 400.0)
        elapsed = time() - start
        sleep(0.1)
        println("Saving file $filepath took $elapsed (s).")
        step = i / frames
        println("Completed step $(100step).\n")

        stitch() = begin
            part = i ÷ (frames / 3)
            exportdir = joinpath(directory, "export")
            !isdir(exportdir) && mkdir(exportdir)
            WxH = "$(resolution[1])x$(resolution[2])"
            #WxH = "1920x1080"
            commonpart = "$(modelname)_$(resolution[2])p_$(FPS)fps"
            inputname = "$(commonpart)_%0$(paddingnumber)d.jpeg"
            inputpath = joinpath(directory, inputname)
            outputname = "$(commonpart)_$(part).mp4"
            outputpath = joinpath(exportdir, outputname)
            command1 = `ffmpeg -y -f image2 -framerate $FPS -i $inputpath -s $WxH -pix_fmt yuvj420p $outputpath`
            run(command1)
        end

        if isapprox(i % (frames / 3), 0)
            stitch()
            sleep(1)
        end
        # break
    end
end

max_stars_repo_name: iamazadi/geometrizer

max_issues_repo_name: iamazadi/Porta

max_forks_repo_name: iamazadi/Porta


text:
module HCAS

using POMDPs
using POMDPModelTools
using GridInterpolations
using LocalFunctionApproximation
using StaticArrays

include("constants.jl")
include("hCAS.jl")
include("transitions.jl")
include("rewards.jl")

end # module

max_stars_repo_name: cherrywoods/HorizontalCAS

max_issues_repo_name: ahmed-irfan/HorizontalCAS

max_forks_repo_name: ahmed-irfan/HorizontalCAS


text:
# test 3: number of trial paths tests
using Distributed;
addprocs(20);
@everywhere include("loadMod.jl");
@everywhere const GUROBI_ENV = Gurobi.Env();
#pmap(i -> importIpopt(),1:30);

NList = [1,5,10,15,20];
dataList = Dict();
iterMax = 20;
NN = 5;

caseList = [13,33,123];
T = 96;
τ = Int64(1/6*T);
Δt = 0.25;
pathTrain = load("pathHist_600.jld");

for ci in 1:length(caseList)
    pathDict = pathTrain["pathDict"][ci][T];
    for j in procs()
        remotecall_fetch(readInData,j,ci,caseList,T,τ);
    end

    cutDict,LBHist,UBHist,UBuHist,UBlHist,timeHist = solveMain(T, Δt, 30, false,false, 2, 2);

    startPGT = time();
    cutDictPG = preGen(T, Δt, NN, iterMax);
    preGenT = time() - startPGT;
    cutDictPGOri = deepcopy(cutDictPG);

    for N in NList
        startT = time();
        # cutDict,LBHist,UBHist,UBuHist,UBlHist,timeHist = solveMain(τ, T, Δt, N, true, false,
        #     Int64(round(100/N)),Int64(round(100/N)), cutDictPG, false, 200, [], 0, pathDict);
        cutDict,LBHist,UBHist,UBuHist,UBlHist,timeHist = solveMain(T, Δt, N, false, false,
            max(Int64(round(500/N)),20), max(Int64(round(500/N)),20), cutDictPG, false, 200, [], 0, pathDict);
        elapsedT = time() - startT;
        dataList[N] = [LBHist,UBHist,UBuHist,UBlHist,timeHist,elapsedT,preGenT];

        save("NResults_$(ci).jld","NOut",dataList);
        cutDictPG = deepcopy(cutDictPGOri);
    end
end


text:
fname = :ARCqKdense
c = Combi(hessian_dense,PDataK,solve_modelKARC,preprocessKARC,decreaseKARC,TparamsKARC())
include("Template.jl")


text:
using Dates

function initialize_data(ideas)
    data = Dict{Symbol, Any}()
    data[:ideas] = sort!(ideas)
    data[:choices] = Vector{Vector{String}}()
    data[:dates] = Vector{DateTime}()
    data[:counts] = zeros(Int, length(ideas))
    return data
end

function update_data(result_idxs, ideas)
    if isfile(joinpath(@__DIR__, "data.bson"))
        data = BSON.load("data.bson")
    else
        data = initialize_data(ideas)
    end
    if sort(data[:ideas]) ≠ sort(ideas)
        error("The test was done while the `ideas.txt` file has changed."*
        "We cannot update the results! Start a new project!")
    end
    push!(data[:choices], [ideas[i] for i in result_idxs])
    push!(data[:dates], Dates.now())
    for i in result_idxs
        data[:counts][i] += 1
    end
    BSON.bson(joinpath(@__DIR__, "data.bson"), copy(data))
    return data
end

# Contributed by Lyndon White, Github tag @oxinabox
# on the official Julia Slack.
module Cursor
    # note that  \u1b  is in decimal 033, which is an escape code
    const ERASE_TO_EOL = "\u1b[K"
    const UP1 = "\u1b[A"
    up(n) = "\u1b[$(n)A"
    const MOVE_TO_SOL = "\r" #Carriage Return
end

function move_cursor_up_while_clearing_lines(io, numlinesup)
    for _ in 1:numlinesup
        print(io, "\r" * "\u1b[K" * "\u1b[A")
    end
    print(io, "\r" * "\u1b[K" * "\r")
end

max_stars_repo_name: JuliaTagBot/WhatDoYouLike

max_issues_repo_name: JuliaTagBot/WhatDoYouLike

max_forks_repo_name: JuliaTagBot/WhatDoYouLike


text:
using Revise
using UArrayInts
using BenchmarkTools

@btime UArrayInts.add_carry($(0xff), $(0xff))
@btime 0x00ff + 0x00ff

@btime UArrayInts.add!($([0xff, 0xff]), $([0xff, 0xff]))
@btime $(0xffff) + $(0xffff)
@btime UArrayInts.add!($([0xff, 0xff]), $([0xff, 0xff]), 1)
@btime UArrayInts.add!($([0xffff, 0xffff]), $([0xffff, 0xffff]))
@btime $(0xffffffff) + $(0xffffffff)
@btime UArrayInts.add!($([0xffff]), $([0xffff, 0xffff]))
@btime (0x00000000000ffff + 0x0000000ffffffff) % UInt16

UArrayInts.mul!([0xffff, 0xffff], 0xffff)

@btime UArrayInts.mul!($([0xffff, 0xffff]), $(0xffff))
@code_warntype UArrayInts.mul!(Array{UInt8}(undef, 2), [0xff, 0xff], [0xff, 0xff])
@code_llvm UArrayInts.mul!(Array{UInt8}(undef, 2), [0xff, 0xff], [0xff, 0xff])
0xffffffff * 0xffff
@btime UArrayInts.mul!($(Array{UInt8}(undef, 2)), $([0xff, 0xff]), $([0xff, 0xff]))
0xffff * 0xffff
@btime $(0xffff) * $(0xffff)
@btime UArrayInts.mul!($([0xffff, 0xffff]), $([0xffff, 0xffff]))
@btime $(0xffffffff) * $(0xffffffff)
@btime UArrayInts.mul!($([0xffff]), $([0xffff, 0xffff]))
@btime (0x00000000000ffff * 0x0000000ffffffff) % UInt16


text:
using DelimitedFiles

export juliaify_address, process_intcode!, elf_intcode_input!, determine_elf_intcode_inputs

juliaify_address(address::Int) = address + 1

function process_intcode!(intcode::Vector{Int})
    instruction_pointer = juliaify_address(0)
    while true
        opcode = intcode[instruction_pointer]
        if opcode == 99
            return
        elseif opcode == 1
            inputaddr1 = juliaify_address(intcode[instruction_pointer + 1])
            inputaddr2 = juliaify_address(intcode[instruction_pointer + 2])
            outputaddr = juliaify_address(intcode[instruction_pointer + 3])
            intcode[outputaddr] = intcode[inputaddr1] + intcode[inputaddr2]
            instruction_pointer += 4
        elseif opcode == 2
            inputaddr1 = juliaify_address(intcode[instruction_pointer + 1])
            inputaddr2 = juliaify_address(intcode[instruction_pointer + 2])
            outputaddr = juliaify_address(intcode[instruction_pointer + 3])
            intcode[outputaddr] = intcode[inputaddr1] * intcode[inputaddr2]
            instruction_pointer += 4
        else
            throw(ArgumentError("Unexpected optcode: $opcode"))
        end
    end
end

elf_intcode_input!(intcode, index, fix) = intcode[index+1] = fix

function determine_elf_intcode_inputs(intcode, result)
    for position1 = 0:99
        for position2 = 0:99
            intcode_ = copy(intcode)
            elf_intcode_input!(intcode_, 1, position1)
            elf_intcode_input!(intcode_, 2, position2)
            process_intcode!(intcode_)
            if intcode_[1] == result
                return position1, position2
            end
        end
    end
end

d02_inputpath = joinpath(@__DIR__, "input.txt")
d02_inputread(path) = vec(readdlm(path, ',', Int))

function d02_answer1(input)
    input = copy(input)
    elf_intcode_input!(input, 1, 12)
    elf_intcode_input!(input, 2, 2)
    process_intcode!(input)
    return input[1]
end

function d02_answer2(input)
    noun, verb = determine_elf_intcode_inputs(input, 19690720)
    return 100 * noun + verb
end


text:
using NewtonsMethod
using Test
using LinearAlgebra, Statistics, Compat, ForwardDiff

@testset "NewtonsMethod.jl" begin

    f(x) = (x-1.0)^3
    f′(x) = 3.0*(x-1.0)^2
    x_0 = 0.1
    @test newtonroot(f, f′, x_0 = x_0).root ≈ 0.9999998168636032
    @test newtonroot(f, x_0 = x_0).root ≈ 0.9999998168636032
    @test newtonroot(f, f′, x_0 = BigFloat(x_0), tolerance = 1E-80).root ≈ BigFloat(1.0)
    @test newtonroot(f, x_0 = BigFloat(x_0), tolerance = 1E-80).root ≈ BigFloat(1.0)

    #maxiter
    @test newtonroot(f, f′, x_0 = x_0, maxiter = 5).root == nothing
    @test newtonroot(f, x_0 = x_0, maxiter = 5).root == nothing
    @test newtonroot(f, f′, x_0 = BigFloat(x_0), tolerance = 1E-80, maxiter = 10).root == nothing
    @test newtonroot(f, x_0 = BigFloat(x_0), tolerance = 1E-80, maxiter = 10).root == nothing

    g(x) = log(x)
    g′(x) = 1.0/x
    x_0 = 2.0
    @test newtonroot(g, g′, x_0 = x_0).root ≈ 1.0
    @test newtonroot(g, x_0 = x_0).root ≈ 1.0
    @test newtonroot(g, g′, x_0 = BigFloat(x_0), tolerance = 1E-50).root ≈ BigFloat(1.0)
    @test newtonroot(g, x_0 = BigFloat(x_0), tolerance = 1E-50).root ≈ BigFloat(1.0)

    @test !(newtonroot(g, g′, x_0 = x_0, tolerance = 1E-2).root ≈ 1.0)
    @test !(newtonroot(g, x_0 = x_0, tolerance = 1E-2).root ≈ 1.0)

    #non-convergence

    h(x) = x^2 + 1.0
    h′(x) = 2x
    x_0 = 10.0
    @test newtonroot(h, h′, x_0 = x_0).root == nothing
    @test newtonroot(h, x_0 = x_0).root == nothing
    @test newtonroot(h, h′, x_0 = BigFloat(x_0), tolerance = 1E-100).root == nothing
    @test newtonroot(h, x_0 = BigFloat(x_0), tolerance = 1E-100).root == nothing

end


text:
#=
This file is part of Pseudospectra.jl.

Julia translation
copyright 2017 Ralph Smith

Portions from original MATLAB code (via EigTool)
Copyright © 2002-2014 Michael Overton and Emre Mengi

SPDX-License-Identifier: BSD-3-Clause
License-Filename: LICENSES/BSD-3-Clause_Eigtool
=#

"""
    psa_radius(A,ϵ [,d]) -> r,z

Compute ϵ-pseudospectral radius for a dense matrix.

Quadratically convergent two-way method to compute the
ϵ-pseudospectral radius `r` of a dense matrix `A`. Also returns a vector
`z` of points where the pseudospectrum intersects the circle of radius `r`.
Uses the "criss-cross" algorithm of Overton and Mengi.

The ϵ-pseudospectral radius is
```
   maximum(abs(z)) for z s.t. minimum(σ(A-zI)) == ϵ
```

Optional arg:

* `d`: eigenvalues of A, if known in advance

Keyword args:

* `verbosity: 0 for quiet, 1 for noise
"""
function psa_radius end

#=
Undocumented keyword args:
* `plotfig`: if a +ve integer, a figure number to use for diagnostic plots (**uses PyPlot**)
* `iterprompt`: whether to pause between steps; for plot visibility

Plotting only works if PyPlot bindings are imported into this module, which we
won't do by default.
=#

function psa_radius(A,epsln,eA=zeros(complex(eltype(A)),0);
                   verbosity=0,plotfig=0,iterprompt=false)

# This is the version from EigTool.
#  It is suited to epsilon not too small (otherwise see Overton's web
#  pages, or perhaps Mengi's) and moderate size dense A
#  For large (esp. sparse) A, see Overton's psapsr code

    if (plotfig > 0)
        if (!isdefined(Main, :PyPlot) || !isdefined(:plot)
            || !(plot === Main.PyPlot.plot))
            @warn("psa_radius: plotting is only implemented for PyPlot")
            plotfig = 0
        end
    end

    n,m = size(A)
    (n==m) || throw(ArgumentError("matrix must be square"))
    (isa(epsln,Real) && (epsln >= 0)) || throw(ArgumentError("ϵ must be >= 0"))

    isempty(eA) && (eA = eigvals(A))

    if plotfig > 0
        figure(plotfig)
        clf()
        plot(real(eA),imag(eA),"k.")
    end

    smalltol = 1e-10 * max(norm(A),epsln)
    purespr = maximum(abs.(eA)) # spectral radius

    if epsln == 0 # degenerate case
        idx = sortperm(-abs.(eA))
        eA = eA[idx]
        f = purespr
        z = [ew for ew in eA if (abs(ew) >= (f-smalltol))]
        return f,z
    else
        if (epsln < 1e-9)
            # TODO: divert to specialized Hamiltonian method
            @warn("epsilon too small; expect poor accuracy")
        end
        rold = -1e-6
        r,ind = findmax(abs.(eA)) # initial iterate
        theta = angle.(eA[ind])
        (length(theta) == 1) && (theta = [theta])
        thetaold = theta
        iter = 0
        no_imageig = false
        thetabest = NaN
        mE = Matrix(epsln * I, size(A)...)
        Areal = (eltype(A) <: Real)
        realtol = smalltol # used to detect zero real parts
        radtol = smalltol # determines how far ew magnitudes can be from unity

        while !no_imageig && (r > rold)
            if verbosity > 0
                @printf("\npsa_radius: r=%22.15f   ",r)
                @printf("\npsa_radius: theta=")
                print(theta)
            end
            iter += 1
            (iter > 20) && error("psa_radius: too many steps")
            thetabestt = thetabest
            rold = r
            # given the resulting directions in theta (computed in previous
            # iteration except when iter=1), look for the circle with the
            # greatest radius intersecting pseudo-spectrum boundary.
            # Note: input theta is a vector, but r is scalar
            r,thetabest = pspr_2way_rad(A,mE,theta,realtol,iter,rold,plotfig,verbosity)
            ptout = "Computing Pseudospectral Radius..(iteration $iter)"
            if r > rold
                # given current r, look for all relevant intersections of the
                # circle with radius r and the pseudospectrum, process and
                # return pair midpoints
                thetaold = theta
                theta = pspr_2way_theta(A,mE,epsln,r,thetabest,iter,radtol,
                                        smalltol, plotfig,verbosity)
                if iterprompt
                    print(" [RET]")
                    readline()
                end
                no_thetaeig = isempty(theta)
                no_thetaeig && (rold = r)
            else
                thetabest = thetabestt
                r = rold
            end
        end # while

        if isempty(thetabest)
            # TODO: fall back to Hamiltonian scheme
            error("Failed in first iteration")
        end

        f = r
        z = [r * (cos(thetabest) + 1im*sin(thetabest))]
        if Areal && !isreal(z[1])
            z = [z[1],conj(z[1])]
        end
        return f,z
    end
end

"""
Given a set of angles in theta, in each direction finds the point
on the pseudospectrum boundary with the largest radius, i.e.

  `rnew[j] = max |z| s.t. `σ_min(A - z I) = epsln and angle(z) = theta[j]`

`rbest` is the maximum of the largest radius in any direction

          `rbest = max rnew[j], 1<=j<=length(theta)`

`thetabest` is the angle of a point on the boudary with radius `rbest` in one of
the directions in `theta`, i.e. let

 `rbest = rnew[k]`, for some k, 1<=k<=length(theta),`

then `thetabest = theta[k]`.
"""
function pspr_2way_rad(A,mE,theta,realtol,iter,rold,plotfig,verbosity)

    n = size(A,1)
    Aisreal = (eltype(A) <: Real)
    rnew = zeros(size(theta)...)
    for j in 1:length(theta)
        K = vcat(hcat(im*exp(im*theta[j])*A',mE),hcat(-mE,im*exp(-im*theta[j])*A))
        eK = eigvals(K)

        if minimum(abs.(real.(eK))) <= realtol # check for real eigenvalue
            rnew[j] = maximum([imag(ew) for ew in eK if (abs(real(ew)) < realtol)])
            if plotfig > 0
                figure(plotfig)
                # plot radial line between old and new points
                hold(true)
                plot(rold*cos(theta[j]), rold*sin(theta[j]),"bx")
                plot([rold*cos(theta[j]), rnew[j]*cos(theta[j])],
                     [rold*sin(theta[j]), rnew[j]*sin(theta[j])],"m-")
                plot([rnew*cos(theta[j])], [rnew[j]*sin(theta[j])],"b+")
                if Aisreal
                    plot(rold*cos(theta[j]), -rold*sin(theta[j]),"bx")
                    plot([rold*cos(theta[j]), rnew[j]*cos(theta[j])],
                         [-rold*sin(theta[j]), -rnew[j]*sin(theta[j])],"m-")
                    plot([rnew*cos(theta[j])], [-rnew[j]*sin(theta[j])],"b+")
                end
            end
        else
            # there may be no point on the boundary in this direction
            rnew[j] = -Inf

        end
    end
    if isempty(rnew)  # CHECKME: should this be !isfinite(rnew)
        error("horizontal search failed: no intersection points found for one of the midpoints")
        # TODO: fall back to alternate version
    end
    rbest, ind = findmax(rnew) # outermost
    thetabest = theta[ind]
    if plotfig > 0
        plot(rbest*cos(thetabest),rbest*sin(thetabest),"b*")
    end
    return rbest, thetabest
end

"""
Given a radius `r`, it first computes the intersection points
of eps-pseudospectrum boundary with the circle with radius `r`.
This is achieved by finding the generalized eigenvalues of the
matrix pencil `F - lambda*G` where

       `F=[-eps*I A;r*I 0], G=[0 r*I; A' -eps*I]`

and then performing a singular value test. The singular value test
is neccessary to eliminate the points for which eps is a singular
value but not the smallest one. Finally the midpoints of two
succesive intersection points on the circle with radius r is
calculated, i.e. let `re^(i*theta_i)` and `re^(i*theta_(i+1))` be the
ith and (i+1)th intersection points, then ith midpoint is
re^(i*(theta_i+theta_(i+1))/2). Special attention is paid to
keep the angles in the interval [-pi,pi). Furthermore, as it
was the case for pseudo-absicca code, we specifically considered
the case when the angle from the previous iteration is contained
in one of the intervals. At the exit thetanew contains the
angles of the midpoints.

"""
function pspr_2way_theta(A,mE,epsln,r,thetawant,iter,radtol,smalltol,
                         plotfig,verbosity)
    svd_check = true
    Areal = (eltype(A) <: Real)
    n = size(A,1)

    # compute generalized eigenvalues of matrix pencil F - λ G
    R = Matrix(r * I, n, n)
    O = zeros(n,n)
    F = vcat(hcat(mE,A),hcat(R,O))
    G = vcat(hcat(O,R),hcat(A',mE))
    eM = eigvals(F,G)

    # extract ews with magnitude (nearly) 1
    eM = [ew for ew in eM if (abs(abs(ew)-1) < radtol)]
    if isempty(eM)
        thetanew = zeros(0)
    else
        idx = sortperm(angle.(eM))
        theta = angle.(eM[idx])
        # Perform singular value test on the points probably on the
        # ϵ-pseudospectrum boundary. The ones actually there have
        # smallest singular value equal to ϵ
        idx2 = Int[]
        for j in 1:length(theta)
            (theta[j] < 0) && (theta[j] += 2π)
            Ashift = A - (r*(cos(theta[j])+im*sin(theta[j]))) * I
            s = svdvals(Ashift)
            minσ,minind = findmin(abs.(s .- epsln))
            (minind == n) && push!(idx2,j)
        end
        nremoved = length(theta) - length(idx2)
        if (nremoved > 0)
            if (verbosity > 0)
                @printf("\npspr_2way_theta: singular value test removed %d ew",nremoved)
            end
            theta = theta[idx2]
        end
        if isempty(theta)
            error("singular value test removed all intersection points: please try a smaller ϵ")
        end

        theta = sort(theta)
        # organize in pairs and take midpoints
        thetanew = zeros(0)
        ind = 0
        # shift thetawant (angle from previous iteration) into canon. interval
        (thetawant < 0) && (thetawant += 2π)

        for j in 1:length(theta)
            thetalow = theta[j]
            if j < length(theta)
                thetahigh = theta[j+1]
            else
                thetahigh = theta[1]+2π
            end
            # before taking the midpoint, if this interval is not very short,
            # check and see if thetawant is in this interval, well away from the
            # end points.  If so, break this pair into two pairs, one above
            # and one below thetawant, and take midpoints of both.
            inttol = 0.01 * (thetahigh - thetalow)

            # needed for last interval
	    if ((thetawant+2π > thetalow + inttol)
                && (thetawant+2*pi < thetahigh - inttol))
	        thetawantt = thetawant + 2π
	    else
	        thetawantt = thetawant
	    end

            if (thetawantt > thetalow + inttol) && (thetawantt < thetahigh - inttol)
                # lower midpoint
	        thetamid = (thetalow + thetawantt)/2

	        # shift thetamid into the interval [-π,π] again
	        (thetamid >= 2π) && (thetamid -= 2π)
	        (thetamid >= π) && (thetamid -= 2π)

	        # remove the midpoint if the minimum singular value is greater
	        # than ϵ, since in this case the midpoint should lie outside the
	        # ϵ-pseudospectrum.
	        if (minimum(svdvals(A-r*exp(im*thetamid)*I)) <= epsln)
		    push!(thetanew, thetamid)
	        end

                # upper midpoint
                thetamid = (thetawantt + thetahigh)/2

	        # shift thetamid into the interval [-π,π] again
	        (thetamid >= 2π) && (thetamid -= 2π)
	        (thetamid >= π) && (thetamid -= 2π)

	        # remove the midpoint if the minimum singular value is greater
	        # than ϵ, since in this case the midpoint should lie outside the
	        # ϵ-pseudospectrum.
	        if (minimum(svdvals(A-r*exp(im*thetamid)*I)) <= epsln)
		    push!(thetanew, thetamid)
	        end

            else
                # otherwise, if thetawant is not in the interval
                # take the midpoint of thetalow and thetahigh
                thetamid = (thetalow + thetahigh)/2

	        # shift thetamid into the interval [-π,π] again
	        (thetamid >= 2π) && (thetamid -= 2π)
	        (thetamid >= π) && (thetamid -= 2π)

	        # remove the midpoint if the minimum singular value is greater than
	        # ϵ, since in this case the midpoint should lie outside the
	        # ϵ-pseudospectrum.
	        if (minimum(svdvals(A-r*exp(im*thetamid)*I)) <= epsln)
		    push!(thetanew, thetamid)
	        end
            end
        end

        if plotfig > 0
            figure(plotfig)
            plot(r*cos(0:0.04:6.4),r*sin(0:0.04:6.4))
            # plot intersection points of circle (radius r) and pseudospectrum
            pointsoncircle = r * (cos.(theta) + im*sin.(theta))
            plot(real(pointsoncircle),imag(pointsoncircle),"ro")
        end

        # DEVNOTE: upstream has this logic but doesn't use ynew
        # Did they mean to replace thetanew?
        #=
        if Areal
            # discard midpoints in lower half plane
            ynew = [yy for yy in thetanew if ((yy >= 0) || (yy == -π))]
        end
        =#
    end
    return thetanew
end


text:
const BESSEL_ORDER_INFO = """
derivatives of Bessel functions with respect to the order are not implemented currently:
https://github.com/JuliaMath/SpecialFunctions.jl/issues/160
"""

const INCOMPLETE_GAMMA_INFO = """
derivatives of the incomplete Gamma functions with respect to parameter `a` are not
implemented currently:
https://github.com/JuliaMath/SpecialFunctions.jl/issues/317
"""

const INCOMPLETE_EXPINT_INFO = """
derivatives of the exponential integral with respect to parameter `ν` are not implemented
currently:
https://github.com/JuliaMath/SpecialFunctions.jl/issues/321
"""

ChainRulesCore.@scalar_rule(airyai(x), airyaiprime(x))
ChainRulesCore.@scalar_rule(airyaiprime(x), x * airyai(x))
ChainRulesCore.@scalar_rule(airybi(x), airybiprime(x))
ChainRulesCore.@scalar_rule(airybiprime(x), x * airybi(x))
ChainRulesCore.@scalar_rule(besselj0(x), -besselj1(x))
ChainRulesCore.@scalar_rule(
    besselj1(x),
    (besselj0(x) - besselj(2, x)) / 2,
)
ChainRulesCore.@scalar_rule(bessely0(x), -bessely1(x))
ChainRulesCore.@scalar_rule(
    bessely1(x),
    (bessely0(x) - bessely(2, x)) / 2,
)
ChainRulesCore.@scalar_rule(dawson(x), 1 - (2 * x * Ω))
ChainRulesCore.@scalar_rule(digamma(x), trigamma(x))
ChainRulesCore.@scalar_rule(erf(x), (2 / sqrt(π)) * exp(-x * x))
ChainRulesCore.@scalar_rule(erfc(x), -(2 / sqrt(π)) * exp(-x * x))
ChainRulesCore.@scalar_rule(erfcinv(x), -(sqrt(π) / 2) * exp(Ω^2))
ChainRulesCore.@scalar_rule(erfcx(x), (2 * x * Ω) - (2 / sqrt(π)))
ChainRulesCore.@scalar_rule(erfi(x), (2 / sqrt(π)) * exp(x * x))
ChainRulesCore.@scalar_rule(erfinv(x), (sqrt(π) / 2) * exp(Ω^2))
ChainRulesCore.@scalar_rule(gamma(x), Ω * digamma(x))
ChainRulesCore.@scalar_rule(
    gamma(a, x),
    (
        ChainRulesCore.@not_implemented(INCOMPLETE_GAMMA_INFO),
        - exp(-x) * x^(a - 1),
    ),
)
ChainRulesCore.@scalar_rule(
    gamma_inc(a, x, IND),
    @setup(z = exp(-x) * x^(a - 1) / gamma(a)),
    (
        ChainRulesCore.@not_implemented(INCOMPLETE_GAMMA_INFO),
        z,
        ChainRulesCore.NoTangent(),
    ),
    (
        ChainRulesCore.@not_implemented(INCOMPLETE_GAMMA_INFO),
        -z,
        ChainRulesCore.NoTangent(),
    ),
)
ChainRulesCore.@scalar_rule(
    invdigamma(x),
    inv(trigamma(invdigamma(x))),
)
ChainRulesCore.@scalar_rule(trigamma(x), polygamma(2, x))

# binary
ChainRulesCore.@scalar_rule(
    besselj(ν, x),
    (
        ChainRulesCore.@not_implemented(BESSEL_ORDER_INFO),
        (besselj(ν - 1, x) - besselj(ν + 1, x)) / 2
    ),
)
ChainRulesCore.@scalar_rule(
    besseli(ν, x),
    (
        ChainRulesCore.@not_implemented(BESSEL_ORDER_INFO),
        (besseli(ν - 1, x) + besseli(ν + 1, x)) / 2,
    ),
)
ChainRulesCore.@scalar_rule(
    bessely(ν, x),
    (
        ChainRulesCore.@not_implemented(BESSEL_ORDER_INFO),
        (bessely(ν - 1, x) - bessely(ν + 1, x)) / 2,
    ),
)
ChainRulesCore.@scalar_rule(
    besselk(ν, x),
    (
        ChainRulesCore.@not_implemented(BESSEL_ORDER_INFO),
        -(besselk(ν - 1, x) + besselk(ν + 1, x)) / 2,
    ),
)
ChainRulesCore.@scalar_rule(
    hankelh1(ν, x),
    (
        ChainRulesCore.@not_implemented(BESSEL_ORDER_INFO),
        (hankelh1(ν - 1, x) - hankelh1(ν + 1, x)) / 2,
    ),
)
ChainRulesCore.@scalar_rule(
    hankelh2(ν, x),
    (
        ChainRulesCore.@not_implemented(BESSEL_ORDER_INFO),
        (hankelh2(ν - 1, x) - hankelh2(ν + 1, x)) / 2,
    ),
)
ChainRulesCore.@scalar_rule(
    polygamma(m, x),
    (
        ChainRulesCore.NoTangent(),
        polygamma(m + 1, x),
    ),
)
# todo: setup for common expr
ChainRulesCore.@scalar_rule(
    beta(a, b),
    (Ω*(digamma(a) - digamma(a + b)),
     Ω*(digamma(b) - digamma(a + b)),)
)
ChainRulesCore.@scalar_rule(
    logbeta(a, b),
    (digamma(a) - digamma(a + b),
     digamma(b) - digamma(a + b),)
)

# actually is the absolute value of the logorithm of gamma paired with sign gamma
ChainRulesCore.@scalar_rule(
    logabsgamma(x),
    digamma(x),
    ChainRulesCore.ZeroTangent()
)

ChainRulesCore.@scalar_rule(loggamma(x), digamma(x))
ChainRulesCore.@scalar_rule(
    loggamma(a, x),
    (
        ChainRulesCore.@not_implemented(INCOMPLETE_GAMMA_INFO),
        -exp(- (x + Ω)) * x^(a - 1),
    )
)

# exponential integrals
ChainRulesCore.@scalar_rule(expint(z), - exp(-z) / z)
ChainRulesCore.@scalar_rule(
    expint(ν, z),
    (
        ChainRulesCore.@not_implemented(INCOMPLETE_EXPINT_INFO),
        - expint(ν - 1, z),
    )
)
ChainRulesCore.@scalar_rule(expintx(z), Ω - inv(z))
ChainRulesCore.@scalar_rule(
    expintx(ν, z),
    (
        ChainRulesCore.@not_implemented(INCOMPLETE_EXPINT_INFO),
        Ω - expintx(ν - 1, z),
    )
)
ChainRulesCore.@scalar_rule(expinti(x), exp(x) / x)
ChainRulesCore.@scalar_rule(sinint(x), sinc(x / π))
ChainRulesCore.@scalar_rule(cosint(x), cos(x) / x)


text:
type Nil end

type List{T}
    ar :: Union{T, List{T}, Nil}
    dr :: Union{T, List{T}, Nil}
end

cons = List
car(l) = l.ar
cdr(l) = l.dr

max_stars_repo_name: sclereid/collection

max_issues_repo_name: sclereid/collection

max_forks_repo_name: sclereid/collection


text:
# misc constraint tests
@testset "node bc" begin
    grid = generate_grid(Triangle, (1, 1))
    addnodeset!(grid, "nodeset", x-> x[2] == -1 || x[1] == -1)
    dh = DofHandler(grid)
    push!(dh, :u, 2)
    push!(dh, :p, 1)
    close!(dh)
    ch = ConstraintHandler(dh)
    dbc1 = Dirichlet(:u, getnodeset(grid, "nodeset"), (x,t) -> x, [1, 2])
    dbc2 = Dirichlet(:p, getnodeset(grid, "nodeset"), (x,t) -> 0, 1)
    add!(ch, dbc1)
    add!(ch, dbc2)
    close!(ch)
    update!(ch)

    @test ch.prescribed_dofs == collect(1:9)
    @test ch.inhomogeneities == [-1, -1, 1, -1, -1, 1, 0, 0, 0]

    ## test node bc with mixed dof handler
   dim = 2
   mesh = generate_grid(QuadraticQuadrilateral, (2,1))
   addcellset!(mesh, "set1", Set(1))
   addcellset!(mesh, "set2", Set(2))
   addnodeset!(mesh, "bottom", Set(1:5))

   dh  = MixedDofHandler(mesh)

   ip_quadratic = Lagrange{dim, RefCube, 2}()
   ip_linear = Lagrange{dim, RefCube, 1}()
   field_u = Field(:u, ip_quadratic, dim)
   field_c = Field(:c, ip_linear, 1)
   push!(dh, FieldHandler([field_u, field_c], getcellset(mesh, "set1")))
   push!(dh, FieldHandler([field_u], getcellset(mesh, "set2")))

   close!(dh)

   ch = ConstraintHandler(dh)
   add!(ch, dh.fieldhandlers[1], Dirichlet(:u, getnodeset(mesh, "bottom"), (x,t)->1.0, 1))
   add!(ch, dh.fieldhandlers[2], Dirichlet(:u, getnodeset(mesh, "bottom"), (x,t)->1.0, 1))
   add!(ch, dh.fieldhandlers[1], Dirichlet(:c, getnodeset(mesh, "bottom"), (x,t)->2.0, 1))
   close!(ch)
   update!(ch)

   @test ch.prescribed_dofs == [1,3,9,19,20,23,27]
   @test ch.inhomogeneities == [1.0, 1.0, 1.0, 2.0, 2.0, 1.0, 1.0]

   ## MixedDofHandler: let first FieldHandler not have all fields
   dim = 2
   mesh = generate_grid(Quadrilateral, (2,1))
   addcellset!(mesh, "set1", Set(1))
   addcellset!(mesh, "set2", Set(2))

   ip = Lagrange{dim, RefCube, 1}()
   field_u = Field(:u, ip, dim)
   field_c = Field(:c, ip, 1)

   dh = MixedDofHandler(mesh)
   push!(dh, FieldHandler([field_u], getcellset(mesh, "set1")))
   push!(dh, FieldHandler([field_u, field_c], getcellset(mesh, "set2")))
   close!(dh)

   ch = ConstraintHandler(dh)
   add!(ch, dh.fieldhandlers[1], Dirichlet(:u, getfaceset(mesh, "bottom"), (x,t)->1.0, 1))
   add!(ch, dh.fieldhandlers[2], Dirichlet(:u, getfaceset(mesh, "bottom"), (x,t)->1.0, 1))
   add!(ch, dh.fieldhandlers[2], Dirichlet(:c, getfaceset(mesh, "bottom"), (x,t)->2.0, 1))
   close!(ch)
   update!(ch)

   @test ch.prescribed_dofs == [1,3,9,13,14]
   @test ch.inhomogeneities == [1.0, 1.0, 1.0, 2.0, 2.0]
end

@testset "edge bc" begin
    grid = generate_grid(Hexahedron, (1, 1, 1))
    addedgeset!(grid, "edge", x-> x[1] ≈ -1.0 && x[3] ≈ -1.0)

    dh = DofHandler(grid)
    push!(dh, :u, 3)
    push!(dh, :p, 1)
    close!(dh)

    ch = ConstraintHandler(dh)
    dbc1 = Dirichlet(:u, getedgeset(grid, "edge"), (x,t) -> x, [1, 2, 3])
    add!(ch, dbc1)
    close!(ch)
    update!(ch)

    @test ch.prescribed_dofs == [1,2,3,10,11,12]
    @test ch.inhomogeneities == [-1.0, -1.0, -1.0, -1.0, 1.0, -1.0]

    #Shell mesh edge bcs
    nodes = [Node{3,Float64}(Vec(0.0,0.0,0.0)), Node{3,Float64}(Vec(1.0,0.0,0.0)),
             Node{3,Float64}(Vec(1.0,1.0,0.0)), Node{3,Float64}(Vec(0.0,1.0,0.0)),
             Node{3,Float64}(Vec(2.0,0.0,0.0)), Node{3,Float64}(Vec(2.0,2.0,0.0))]

    cells = [Quadrilateral3D((1,2,3,4)), Quadrilateral3D((2,5,6,3))]
    grid = Grid(cells,nodes)

    #3d quad with 1st order 2d interpolation
    dh = DofHandler(grid)
    push!(dh, :u, 1, Lagrange{2,RefCube,2}())
    push!(dh, :θ, 1, Lagrange{2,RefCube,2}())
    close!(dh)

    addedgeset!(grid, "edge", x -> x[2] ≈ 0.0) #bottom edge
    ch = ConstraintHandler(dh)
    dbc1 = Dirichlet(:θ, getedgeset(grid, "edge"), (x,t) -> (0.0,), [1])
    add!(ch, dbc1)
    close!(ch)
    update!(ch)

    @test ch.prescribed_dofs == [10, 11, 14, 25, 27]
end

@testset "affine constraints" begin

    grid = generate_grid(Line, (10,))
    dh = DofHandler(grid)
    push!(dh, :u, 1)
    close!(dh)

    test_acs = [
        # Simple homogeneous constraint
        [AffineConstraint(4, [(7 => 1.0)], 0.0)],
        # Two dofs and inhomogeneity
        [AffineConstraint(2, [(5 => 1.0), (6 =>2.0)], 1.0)],
        # Two linear constraints
        [AffineConstraint(2, [9=>1.0], 0.0),
         AffineConstraint(3, [9=>1.0], 0.0)],
        #
        [AffineConstraint(2, [7=>3.0, 8=>1.0], -1.0),
         AffineConstraint(4, [9=>-1.0], 2.0)]
    ]

    for acs in test_acs
        ch = ConstraintHandler(dh)
        add!(ch, Dirichlet(:u, getfaceset(grid, "left"), (x,t)->0.0))
        for lc in acs
            add!(ch, lc)
        end
        close!(ch)
        update!(ch, 0.0)
        C, g = Ferrite.create_constraint_matrix(ch)

        # Assemble
        K = create_sparsity_pattern(dh, ch)
        f = zeros(ndofs(dh)); f[end] = 1.0
        for cell in CellIterator(dh)
            K[celldofs(cell), celldofs(cell)] += 2.0 * [1 -1; -1 1]
        end
        copies = (K = copy(K), f1 = copy(f), f2 = copy(f))

        # Solve by actually condensing the matrix
        ff  = C' * (f - K * g)
        KK = C' * K * C
        _aa = KK \ ff
        aa = C * _aa + g

        # Solving by modifying K inplace
        apply!(K, f, ch)
        a = K \ f
        apply!(a, ch)

        # Solve with extracted RHS data
        rhs = get_rhs_data(ch, copies.K)
        apply!(copies.K, ch)
        apply_rhs!(rhs, copies.f1, ch)
        a_rhs1 = apply!(copies.K \ copies.f1, ch)
        apply_rhs!(rhs, copies.f2, ch)
        a_rhs2 = apply!(copies.K \ copies.f2, ch)

        @test a ≈ aa ≈ a_rhs1 ≈ a_rhs2
    end

end