text:
using OREnvironment
using ORInterface
using Test

#######################
# PRELIMINARIES FOR THE TESTS
#######################
mutable struct MyProblem <: Problem
    costs::Array{Float64,1};
    constraints::Array{<:Constraint,1};
    variablesConstraints::Array{Array{Int,1},1};
end

function constructorProblem()
    cost = collect(1.0:6.0);

    variables1 = [1, 3, 4, 6];
    variables2 = [1, 3, 5, 6];
    coefficients1 = [2.3, 3.2, 3.1, 12.34];
    coefficients2 = coefficients1 .+ 1.0;
    constraint1 = OREnvironment.constructConstraint(15.0, :lessOrEq, variables1, coefficients1);
    constraint2 = OREnvironment.constructConstraint(9.0, :lessOrEq, variables2, coefficients2);
    constraints = [constraint1, constraint2];

    variablesConstraints = [[1,2], Int[], [1,2], [1], [2], [1,2]];

    return MyProblem(cost, constraints, variablesConstraints);
end

function constructorSolution(numConstraints=2)
    status = OREnvironment.constructStatus(numConstraints);
    Tvariables = Float64; sizeArray = 6;
    sol  = OREnvironment.constructSolution(:FixedLengthArray, (Tvariables, sizeArray, status));
    return sol;
end

#######################
# TESTS
#######################
@testset "Construct constraint and getters and setters" begin
    variables = [1, 10, 23, 56];
    coefficients = [2.3, 3.2, 3.1, 12.34];
    constraint = OREnvironment.constructConstraint(10.0, :greaterOrEq, variables, coefficients);

    @test OREnvironment.get_rhs(constraint) == 10.0;
    @test OREnvironment.get_type(constraint) == :greaterOrEq;
    @test OREnvironment.get_coefficient(constraint, 1)  == 2.3;
    @test OREnvironment.get_coefficient(constraint, 56) == 12.34;
    @test OREnvironment.get_coefficient(constraint, 4)  == 0.0;
    @test OREnvironment.is_variable(constraint, 10) == true;
    @test OREnvironment.is_variable(constraint, 11) == false;

    OREnvironment.set_rhs!(constraint, 12.6);
    @test OREnvironment.get_rhs(constraint) == 12.6;
    OREnvironment.set_type!(constraint, :equal);
    @test OREnvironment.get_type(constraint) == :equal;
    OREnvironment.set_coefficient!(constraint, 23, 190.8);
    @test OREnvironment.get_coefficient(constraint, 23) == 190.8;
end

@testset "is_feasible" begin
    variables = [1, 10, 23, 56];
    coefficients = [2.3, 3.2, 3.1, 12.34];
    constraint = OREnvironment.constructConstraint(15.0, :lessOrEq, variables, coefficients);

    @testset " <= constraints" begin
       lhs = 13.2;
       @test OREnvironment.is_feasible(constraint, lhs) == true;
       lhs = 15.0;
       @test OREnvironment.is_feasible(constraint, lhs) == true;
       lhs = 20.0;
       @test OREnvironment.is_feasible(constraint, lhs) == false;
    end
    @testset " => constraints" begin
       lhs = 13.2;
       OREnvironment.set_type!(constraint, :greaterOrEq);
       @test OREnvironment.is_feasible(constraint, lhs) == false;
       lhs = 15.0;
       @test OREnvironment.is_feasible(constraint, lhs) == true;
       lhs = 20.0;
       @test OREnvironment.is_feasible(constraint, lhs) == true;
    end
    @testset " == constraints" begin
       lhs = 13.2;
       OREnvironment.set_type!(constraint, :equal);
       @test OREnvironment.is_feasible(constraint, lhs) == false;
       lhs = 15.0;
       @test OREnvironment.is_feasible(constraint, lhs) == true;
       lhs = 20.0;
       @test OREnvironment.is_feasible(constraint, lhs) == false;
    end
end

@testset "is_active" begin
    variables = [1, 10, 23, 56];
    coefficients = [2.3, 3.2, 3.1, 12.34];
    constraint = OREnvironment.constructConstraint(15.0, :lessOrEq, variables, coefficients);

    @testset " <= constraints" begin
       lhs = 13.2;
       @test OREnvironment.is_active(constraint, lhs) == false;
       lhs = 15.0;
       @test OREnvironment.is_active(constraint, lhs) == true;
    end
    @testset " => constraints" begin
       lhs = 13.2;
       OREnvironment.set_type!(constraint, :greaterOrEq);
       @test OREnvironment.is_active(constraint, lhs) == false;
       lhs = 15.0;
       @test OREnvironment.is_active(constraint, lhs) == true;
    end
    @testset " == constraints" begin
       lhs = 13.2;
       OREnvironment.set_type!(constraint, :equal);
       @test OREnvironment.is_active(constraint, lhs) == false;
       lhs = 15.0;
       @test OREnvironment.is_active(constraint, lhs) == true;
    end
end

@testset "add_constraint_index_to_variables" begin
    p = constructorProblem();
    variablesConstraints = [Int[] for i in 1:6];
    OREnvironment.add_constraint_index_to_variables!(p.constraints[1], 1, variablesConstraints);
    @test variablesConstraints[1] == [1];
    @test variablesConstraints[2] == Int[];
    @test variablesConstraints[3] == [1];
    @test variablesConstraints[4] == [1];
    @test variablesConstraints[5] == Int[];
    @test variablesConstraints[6] == [1];
    OREnvironment.add_constraint_index_to_variables!(p.constraints[2], 2, variablesConstraints);
    for i in 1:6
      @test variablesConstraints[i] == p.variablesConstraints[i];
    end
end

@testset "get_relationship_variables_constraints" begin
    p = constructorProblem();
    variablesConstraints = OREnvironment.get_relationship_variables_constraints(p.constraints, 6);

    for i in 1:6
        if length(variablesConstraints[i]) == 0
            @test variablesConstraints[i] == Int[];
        else
            for j in 1:length(variablesConstraints[i])
                @test variablesConstraints[i][j] == p.variablesConstraints[i][j];
            end
        end
    end

    # when passing no constraints
    dummy = Array{OREnvironment.DefaultConstraint,1}();
    variablesConstraints = OREnvironment.get_relationship_variables_constraints(dummy, 6);
    @test variablesConstraints == Array{Array{Int,1}, 1}();
end

@testset "parse_and_set_constraint_sign" begin
  l1 = "12x1 + 13x2 + 5x3 <= 10";
  l2 = "1x2 + 23.2x14 + 5.3x6 >= 12";
  l3 = "2x5 + 3x7 + 5x19 = 15";
  lwithError1 = "1x2 + 23x4 + 5x6 => 12";
  lwithError2 = "1x2 + 23x4 + 5x6 == 12";
  c = OREnvironment.constructConstraint(0.0, :equal, Int[], Float64[]);

  sign = OREnvironment.parse_and_set_constraint_sign!(c, l1, 1);
  @test sign == "<=";
  @test OREnvironment.get_type(c) == :lessOrEq;
  sign = OREnvironment.parse_and_set_constraint_sign!(c, l2, 2);
  @test sign == ">=";
  @test OREnvironment.get_type(c) == :greaterOrEq;
  sign = OREnvironment.parse_and_set_constraint_sign!(c, l3, 3);
  @test sign == "=";
  @test OREnvironment.get_type(c) == :equal;
  @test_throws ErrorException OREnvironment.parse_and_set_constraint_sign!(c, lwithError1, 4);
  @test_throws ErrorException OREnvironment.parse_and_set_constraint_sign!(c, lwithError2, 5);
end

@testset "parse_constraint" begin
  l1 = "12x1 + 13x2 + 5x3 <= 10";
  l2 = "1x2 + 23.2x14 + 5.3x6 >= 12";
  l3 = "2x5 + 3x7 + 5x19 = 15";
  lwithError1 = "1x2 + 23x4 + 5x6 => 12";
  lwithError2 = "1x2 + 23x4 + 5x6 == 12";
  c = OREnvironment.constructConstraint(0.0, :equal, Int[], Float64[]);

  OREnvironment.parse_constraint!(c, l1, 1);
  @test OREnvironment.get_rhs(c) == 10;
  @test OREnvironment.get_type(c) == :lessOrEq;
  @test OREnvironment.is_variable(c, 1) == true;
  @test OREnvironment.is_variable(c, 2) == true;
  @test OREnvironment.is_variable(c, 3) == true;
  @test OREnvironment.get_coefficient(c, 1) == 12.0;
  @test OREnvironment.get_coefficient(c, 2) == 13.0;
  @test OREnvironment.get_coefficient(c, 3) == 5.0;

  OREnvironment.parse_constraint!(c, l2, 2);
  @test OREnvironment.get_rhs(c) == 12;
  @test OREnvironment.get_type(c) == :greaterOrEq;
  @test OREnvironment.is_variable(c, 2) == true;
  @test OREnvironment.is_variable(c, 14) == true;
  @test OREnvironment.is_variable(c, 6) == true;
  @test OREnvironment.get_coefficient(c, 2) == 1.0;
  @test OREnvironment.get_coefficient(c, 14) == 23.2;
  @test OREnvironment.get_coefficient(c, 6) == 5.3;

  OREnvironment.parse_constraint!(c, l3, 3);
  @test OREnvironment.get_rhs(c) == 15;
  @test OREnvironment.get_type(c) == :equal;
  @test OREnvironment.is_variable(c, 5) == true;
  @test OREnvironment.is_variable(c, 7) == true;
  @test OREnvironment.is_variable(c, 19) == true;
  @test OREnvironment.get_coefficient(c, 5) == 2.0;
  @test OREnvironment.get_coefficient(c, 7) == 3.0;
  @test OREnvironment.get_coefficient(c, 19) == 5.0;

  @test_throws ErrorException OREnvironment.parse_constraint!(c, lwithError1, 4);
  @test_throws ErrorException OREnvironment.parse_constraint!(c, lwithError2, 5);
end

@testset "read_constraints" begin
  testdir = dirname(@__FILE__);
  file = joinpath(testdir, "sampleConstraints.txt");
  constraints = OREnvironment.read_constraints(file);

  c = constraints[1];
  @test OREnvironment.get_rhs(c) == 9;
  @test OREnvironment.get_type(c) == :lessOrEq;
  @test OREnvironment.is_variable(c, 1) == true;
  @test OREnvironment.is_variable(c, 5) == true;
  @test OREnvironment.is_variable(c, 17) == true;
  @test OREnvironment.get_coefficient(c, 1) == 2.0;
  @test OREnvironment.get_coefficient(c, 5) == 3.0;
  @test OREnvironment.get_coefficient(c, 17) == 2.5;

  c = constraints[2];
  @test OREnvironment.get_rhs(c) == 12;
  @test OREnvironment.get_type(c) == :greaterOrEq;
  @test OREnvironment.is_variable(c, 24) == true;
  @test OREnvironment.is_variable(c, 122) == true;
  @test OREnvironment.is_variable(c, 3) == true;
  @test OREnvironment.is_variable(c, 1) == true;
  @test OREnvironment.get_coefficient(c, 24) == -3.0;
  @test OREnvironment.get_coefficient(c, 122) == 5.2;
  @test OREnvironment.get_coefficient(c, 3) == -2.1;
  @test OREnvironment.get_coefficient(c, 1) == 1;

  c = constraints[3];
  @test OREnvironment.get_rhs(c) == 4;
  @test OREnvironment.get_type(c) == :equal;
  @test OREnvironment.is_variable(c, 9) == true;
  @test OREnvironment.is_variable(c, 15) == true;
  @test OREnvironment.is_variable(c, 32) == true;
  @test OREnvironment.get_coefficient(c, 9) == 3.0;
  @test OREnvironment.get_coefficient(c, 15) == 23.0;
  @test OREnvironment.get_coefficient(c, 32) == -10.0;
end

@testset "is_increment_feasible and compute_lhs_after_increment" begin
    s = constructorSolution();
    p = constructorProblem();

    @testset "providing a variable that is not in the constraint" begin
        idx = 1;
        variable = 2;
        Δvariable = 12.0;
        currentLHS = OREnvironment.get_constraint_consumption(s, idx);
        lhs = OREnvironment.compute_lhs_after_increment(variable, Δvariable, currentLHS, p.constraints[idx]);
        feasible = OREnvironment.is_increment_feasible(variable, Δvariable, currentLHS, p.constraints[idx]);
        @test feasible == true;
        @test lhs == 0.0;
        idx = 2;
        currentLHS = OREnvironment.get_constraint_consumption(s, idx);
        lhs = OREnvironment.compute_lhs_after_increment(variable, Δvariable, currentLHS, p.constraints[idx]);
        feasible = OREnvironment.is_increment_feasible(variable, Δvariable, currentLHS, p.constraints[idx]);
        @test feasible == true;
        @test lhs == 0.0;
    end

    @testset "variable that is in both constraints" begin
        idx = 1;
        variable = 1;
        Δvariable = 5.0;
        currentLHS = OREnvironment.get_constraint_consumption(s, idx);
        lhs = OREnvironment.compute_lhs_after_increment(variable, Δvariable, currentLHS, p.constraints[idx]);
        feasible = OREnvironment.is_increment_feasible(variable, Δvariable, currentLHS, p.constraints[idx]);
        @test feasible == true;
        @test lhs == 11.5;
        idx = 2;
        currentLHS = OREnvironment.get_constraint_consumption(s, idx);
        lhs = OREnvironment.compute_lhs_after_increment(variable, Δvariable, currentLHS, p.constraints[idx]);
        feasible = OREnvironment.is_increment_feasible(variable, Δvariable, currentLHS, p.constraints[idx]);
        @test feasible == false;
        @test lhs == 16.5;
    end
end

#######################
# TESTS WHEN DEALING WITH SOLUTIONS
#######################
@testset "is_increment_feasible" begin
    s = constructorSolution();
    p = constructorProblem();

    variable = 1;
    Δvariable = 1.0;
    @test OREnvironment.is_increment_feasible(s, p.constraints, variable, Δvariable, p.variablesConstraints[variable]) == true;
    Δvariable = 3.0;
    @test OREnvironment.is_increment_feasible(s, p.constraints, variable, Δvariable, p.variablesConstraints[variable]) == false;

    variable = 5;
    Δvariable = 1.0;
    @test OREnvironment.is_increment_feasible(s, p.constraints, variable, Δvariable, p.variablesConstraints[variable]) == true;
    Δvariable = 3.0;
    @test OREnvironment.is_increment_feasible(s, p.constraints, variable, Δvariable, p.variablesConstraints[variable]) == false;

    # case with no constraints
    variable = 2;
    Δvariable = 10.0;
    @test OREnvironment.is_increment_feasible(s, p.constraints, variable, Δvariable, p.variablesConstraints[variable]) == true;
    Δvariable = 30.0;
    @test OREnvironment.is_increment_feasible(s, p.constraints, variable, Δvariable, p.variablesConstraints[variable]) == true;
end

@testset "is_current_consumption_feasible" begin
    s = constructorSolution();
    p = constructorProblem();

    OREnvironment.set_constraint_consumption!(s, 1, 7.0);
    OREnvironment.set_constraint_consumption!(s, 2, 7.0);
    feasible = OREnvironment.is_current_consumption_feasible(s, p.constraints);
    @test feasible == true;

    OREnvironment.set_constraint_consumption!(s, 1, 17.0);
    feasible = OREnvironment.is_current_consumption_feasible(s, p.constraints);
    @test feasible == false;

    OREnvironment.set_constraint_consumption!(s, 1, 7.0);
    OREnvironment.set_constraint_consumption!(s, 2, 17.0);
    feasible = OREnvironment.is_current_consumption_feasible(s, p.constraints);
    @test feasible == false;

    OREnvironment.set_constraint_consumption!(s, 1, 37.0);
    feasible = OREnvironment.is_current_consumption_feasible(s, p.constraints);
    @test feasible == false;

    OREnvironment.set_constraint_consumption!(s, 1, 6.0);
    OREnvironment.set_constraint_consumption!(s, 2, 5.0);
    feasible = OREnvironment.is_current_consumption_feasible(s, p.constraints);
    @test feasible == true;
end

@testset "compute_lhs and is_feasible" begin
    s = constructorSolution();
    p = constructorProblem();

    OREnvironment.add_solution!(s, 1, 2.0);
    lhs = OREnvironment.compute_lhs(p.constraints[1], s);
    @test lhs == 4.6;
    lhs = OREnvironment.compute_lhs(p.constraints[2], s);
    @test lhs == 6.6;
    @test OREnvironment.is_feasible(s, p.constraints) == true;

    OREnvironment.add_solution!(s, 2, 4.5);
    lhs = OREnvironment.compute_lhs(p.constraints[1], s);
    @test lhs == 4.6;
    lhs = OREnvironment.compute_lhs(p.constraints[2], s);
    @test lhs == 6.6;
    @test OREnvironment.is_feasible(s, p.constraints) == true;

    OREnvironment.add_solution!(s, 3, 1.5);
    lhs = OREnvironment.compute_lhs(p.constraints[1], s);
    @test lhs == 9.4;
    lhs = OREnvironment.compute_lhs(p.constraints[2], s);
    @test lhs == 12.9;
    @test OREnvironment.is_feasible(s, p.constraints) == false;

    OREnvironment.remove_all_solutions!(s);
    lhs = OREnvironment.compute_lhs(p.constraints[1], s);
    @test lhs == 0.0;
    lhs = OREnvironment.compute_lhs(p.constraints[2], s);
    @test lhs == 0.0;
    @test OREnvironment.is_feasible(s, p.constraints) == true;
end

@testset "compute_lhs and is_active" begin
    s = constructorSolution();
    p = constructorProblem();

    OREnvironment.add_solution!(s, 1, 2.0);
    @test OREnvironment.is_active(p.constraints[1], s) == false; # lhs = 4.6
    @test OREnvironment.is_active(p.constraints[2], s) == false; # lhs = 6.6

    OREnvironment.add_solution!(s, 3, 1.5);
    @test OREnvironment.is_active(p.constraints[1], s) == false; # lhs = 9.4
    @test OREnvironment.is_active(p.constraints[2], s) == false; # lhs = 12.9

    OREnvironment.add_solution!(s, 4, 5.6/3.1);
    @test OREnvironment.is_active(p.constraints[1], s) == true; # lhs = 15.0
    @test OREnvironment.is_active(p.constraints[2], s) == false; # lhs = 12.9
end

@testset "is_active_under_current_consumption" begin
    s = constructorSolution();
    p = constructorProblem();

    OREnvironment.set_constraint_consumption!(s, 1, 12.3);
    OREnvironment.set_constraint_consumption!(s, 2, 9.0);
    @test OREnvironment.is_active_under_current_consumption(p.constraints[1], 1, s) == false;
    @test OREnvironment.is_active_under_current_consumption(p.constraints[2], 2, s) == true;
end


text:
# Based on https://www.tandfonline.com/doi/abs/10.1080/01966324.1987.10737221
# TODO: make a package out of this function. Kombinator?

function _ossb_exact_convex_combination(instance::CombinatorialInstance{T}, t::Dict{T, Float64}) where T
  # Input: a combinatorial problem, a (continuous) solution to the OSSB problem,
  # *scaled so that v=1*. This ensures that the point is within the polytope.

  # The solution to return: a set of solutions (each of them being a set of arms),
  # their weights in the convex combination.
  # This representation assumes that all vertices of the polytope are actual solutions.
  points = Vector{Int}[]
  weights = Float64[]

  # Convert the current point t as a vector (this is required to perform linear algebra).
  arms = all_arm_indices(instance) # Assimilate an arm and its index.
  idx_arm = Dict(i => arms[i] for i in 1:length(arms))
  d = instance.n_arms
  x = zeros(Float64, d)
  for (arm, weight) in t
    x[idx_arm[arm]] = weight
  end

  # Counters for the algorithm.
  point = t # What remains to decompose.
  δ = 1.0 # Weight still to assign.

  # Start the algorithm.

  # Find the currently tight constraints to create the matrix G.
  _, tight_matrix = _ossb_find_tight_constraints(instance, x, idx_arm)

  while rank(tight_matrix) < d
    # Find a new vertex by decomposing (what's left of) x.
    y_bar = x
    while rank(tight_matrix) < d
      # Find a solution dir such that G dir = 0 with dir ≠ 0.
      if length(tight_matrix) == 0 # Easier to compute than rank(G) == 0
        # No constraint to satisfy. Take an arbitrary "solution". This is a pathologically
        # underdetermined system.
        # Normalise by dividing by the norm.
        dir = ones(d) ./ sqrt(d)
      else
        # This system is underdetermined, so that the usual syntax A\b is not guaranteed to work
        # (Julia will return a minimum-norm solution in the underdetermined case, which is very likely zero).
        dir = vec(nullspace(tight_matrix)[:, 1])

        # Also normalise this vector, in case the previous computation gave strange results.
        n = norm(dir, 2)
        @assert n > 1.0e-6 # Just in case...
        dir ./= n
      end

      # Compute a point y = y_bar - α d that touches a face of the polyhedron (at least one more
      # constraint is tight there).
      α_val, y_bar = _ossb_touch_facet(copy(instance), y_bar, dir, idx_arm)

      # Update tight_matrix for this vertex-to-be (or already a vertex).
      _, tight_matrix = _ossb_find_tight_constraints(instance, y_bar, idx_arm)
    end

    # Retrieve that point from x to continue the decomposition.
    dir = x - y_bar
    # @show dir
    γ_val, x = _ossb_touch_facet(copy(instance), x, dir, idx_arm)
    δ_scale = γ_val / (1 + γ_val)

    # This point y is a vertex that plays a role in the convex decomposition.
    push!(points, round.(Int, y_bar))
    push!(weights, δ_scale * δ)

    δ *= 1 - δ_scale

    # Update tight_matrix for the rest of the point to decompose.
    _, tight_matrix = _ossb_find_tight_constraints(instance, x, idx_arm)
  end

  # The last iterate x must be an extreme point, as the rank of the tight constraints is d.
  push!(points, round.(Int, x))
  push!(weights, δ)

  return points, weights
end

function _ossb_find_tight_constraints(instance::CombinatorialInstance{T}, x::Vector{Float64}, idx_arm::Dict{T, Int}) where T
  # Assumption: point within the polytope; the combinatorial instance has an LP formulation.

  # Get the LP formulation, without objective function (_).
  m, _, vars = get_lp_formulation(instance, Dict(arm => 0.0 for arm in all_arm_indices(instance)))

  # vars is a dictionary of arms to variables, but we will need the reverse operation.
  var_to_arm = Dict(values(vars) .=> keys(vars))

  # Iterate over the constraints.
  tight_single = JuMP.ConstraintRef[]
  tight_linear = JuMP.ConstraintRef[]
  for (F, S) in list_of_constraint_types(m)
    # Ignore the case where this is an integrality constraint (it does not add bounds). Simple bounds are required!
    if S == MOI.Integer
      continue
    end

    if S == MOI.ZeroOne
      @assert F == JuMP.VariableRef # TODO: is MOI.SingleVariable possible here? Other functions should not be allowed.
      for c in all_constraints(m, F, S)
        lhs = value(jump_function(constraint_object(c)), vr -> x[var_to_arm[vr]])
        if abs(lhs) <= 1.0e-6 || abs(lhs - 1.0) <= 1.0e-6
          push!(tight_single, c)
        end
      end
    else
      for c in all_constraints(m, F, S)
        # If the constraint is tight at t, keep it.
        lhs = value(c, vr -> x[var_to_arm[vr]])
        rhs = MOI.constant(moi_set(constraint_object(c)))
        if abs(lhs - rhs) <= 1.0e-6
          push!(tight_linear, c)
        end
      end
    end
  end

  # Retrieve the coefficients for each constraint.
  var_to_idx = Dict(vars[arm] => idx for (arm, idx) in idx_arm) # Must be rebuild for each new JuMP model.

  tight_to_coefficients = Dict{JuMP.ConstraintRef, Vector{Float64}}()
  sizehint!(tight_to_coefficients, length(tight_single) + length(tight_linear))

  nvars = MOI.get(m, MOI.NumberOfVariables())

  for c in tight_single
    tight_to_coefficients[c] = spzeros(nvars)
    var = jump_function(constraint_object(c))
    tight_to_coefficients[c][var_to_idx[var]] = 1.0
  end
  for c in tight_linear
    tight_to_coefficients[c] = spzeros(nvars)
    for (var, coeff) in jump_function(constraint_object(c)).terms
      tight_to_coefficients[c][var_to_idx[var]] = coeff
    end
  end

  # The returned set may very well be empty, if the solution hit the bound on v.
  tight_matrix = (length(tight_to_coefficients) > 0) ? vcat(collect(values(tight_to_coefficients))'...) : zeros(0, 0)
  return tight_to_coefficients, tight_matrix
end

function _ossb_touch_facet(instance::CombinatorialInstance{T}, x::Vector{Float64}, dir::Vector{Float64}, idx_arm::Dict{T, Int}) where T
  # Maximise α such that ȳ = y - 2 α dir ∈ polytope(instance), where y = x - dir.
  # Return this optimum α and ȳ = y - 2 α⋆ dir.
  d = length(idx_arm)

  # First, get the standard formulation, without objective function (_).
  msub, _, varssub = get_lp_formulation(instance, Dict(arm => 0.0 for arm in all_arm_indices(instance)))
  var_to_idx = Dict(varssub[arm] => idx for (arm, idx) in idx_arm)

  # Then, build a similar model for this specific problem.
  # TODO: could also do it by inspection of all constraints.
  m = Model(instance.solver.solver) # TODO: only for LP-based solvers.
  # TODO: write a dedicated solver for this (highly simplified) problem? It could work by simply shortening an interval, iterating through the constraints.
  set_silent(m) # TODO: maybe allow to configure this?
  @variable(m, α)
  @objective(m, Max, α)

  # Iterate through the constraints and build the new ones.
  y_bar = x .+ α .* dir

  for (F, S) in list_of_constraint_types(msub)
    # Ignore the case where this is an integrality constraint (it does not add bounds).
    if S == MOI.Integer
      continue
    end

    if S == MOI.ZeroOne
      @assert F == JuMP.VariableRef # TODO: is MOI.SingleVariable possible here? Other functions should not be allowed.
      for c in all_constraints(msub, F, S)
        var = jump_function(constraint_object(c)) # A simple variable reference.
        lhs = y_bar[var_to_idx[var]]
        @constraint(m, 0 <= lhs <= 1)
      end
    else
      for c in all_constraints(msub, F, S)
        lhs = AffExpr() + jump_function(constraint_object(c)).constant # The constant should be zero (only present in MOI set).
        for (var, coeff) in jump_function(constraint_object(c)).terms
          lhs += coeff * y_bar[var_to_idx[var]]
        end
        @constraint(m, lhs ∈ moi_set(constraint_object(c)))
      end
    end
  end

  # Solve the new problem. (It should be really easy, as there is only one variable.)
  optimize!(m)
  α_val = value(α)
  y_bar_val = value.(y_bar)

  return α_val, y_bar_val
end

ext: jl


text:
function petk06()
#
# M. Petkovic testing polynomials, page 146
#
    y = [-1.0*[1;1;1;1];3*[1;1;1];-im;-im];

    p1 = reverse(poly(y).a)
    p2 = [1.0;-2;5];
    p2 = conv(p2,p2);
    p = conv(p1,p2);
    z = [-1.0 4; 3 3; -im 2; 1+2*im 2; 1-2*im 2];

    p, PolyZeros(z)
end


text:
"""
    Link

An abstract type whose subtypes determine methods for [`linkfun`](@ref), [`linkinv`](@ref),
[`mueta`](@ref), and [`inverselink`](@ref).
"""
abstract type Link end

# Make links broadcast like a scalar
Base.Broadcast.broadcastable(l::Link) = Ref(l)

"""
    Link01

An abstract subtype of [`Link`](@ref) which are links defined on (0, 1)
"""
abstract type Link01 <: Link end

"""
    CauchitLink

A [`Link01`](@ref) corresponding to the standard Cauchy distribution,
[`Distributions.Cauchy`](https://juliastats.github.io/Distributions.jl/stable/univariate.html#Distributions.Cauchy).
"""
mutable struct CauchitLink <: Link01 end

"""
    CloglogLink

A [`Link01`](@ref) corresponding to the extreme value (or log-Wiebull) distribution.  The
link is the complementary log-log transformation, `log(1 - log(-μ))`.
"""
mutable struct CloglogLink  <: Link01 end

"""
    IdentityLink

The canonical [`Link`](@ref) for the `Normal` distribution, defined as `η = μ`.
"""
mutable struct IdentityLink <: Link end

"""
    InverseLink

The canonical [`Link`](@ref) for [`Distributions.Gamma`](https://juliastats.github.io/Distributions.jl/stable/univariate.html#Distributions.Gamma) distribution, defined as `η = inv(μ)`.
"""
mutable struct InverseLink  <: Link end

"""
    InverseSquareLink

The canonical [`Link`](@ref) for [`Distributions.InverseGaussian`](https://juliastats.github.io/Distributions.jl/stable/univariate.html#Distributions.InverseGaussian) distribution, defined as `η = inv(abs2(μ))`.
"""
mutable struct InverseSquareLink  <: Link end

"""
    LogitLink

The canonical [`Link01`](@ref) for [`Distributions.Bernoulli`](https://juliastats.github.io/Distributions.jl/stable/univariate.html#Distributions.Bernoulli) and [`Distributions.Binomial`](https://juliastats.github.io/Distributions.jl/stable/univariate.html#Distributions.Binomial).
The inverse link, [`linkinv`](@ref), is the c.d.f. of the standard logistic distribution,
[`Distributions.Logistic`](https://juliastats.github.io/Distributions.jl/stable/univariate.html#Distributions.Logistic).
"""
mutable struct LogitLink <: Link01 end

"""
    LogLink

The canonical [`Link`](@ref) for [`Distributions.Poisson`](https://juliastats.github.io/Distributions.jl/stable/univariate.html#Distributions.Poisson), defined as `η = log(μ)`.
"""
mutable struct LogLink <: Link end

"""
    NegativeBinomialLink

The canonical [`Link`](@ref) for [`Distributions.NegativeBinomial`](https://juliastats.github.io/Distributions.jl/stable/univariate.html#Distributions.NegativeBinomial) distribution, defined as `η = log(μ/(μ+θ))`.
The shape parameter θ has to be fixed for the distribution to belong to the exponential family.
"""
mutable struct NegativeBinomialLink  <: Link
    θ::Float64
end

"""
    ProbitLink

A [`Link01`](@ref) whose [`linkinv`](@ref) is the c.d.f. of the standard normal
distribution, [`Distributions.Normal()`](https://juliastats.github.io/Distributions.jl/stable/univariate.html#Distributions.Normal).
"""
mutable struct ProbitLink <: Link01 end

"""
    SqrtLink

A [`Link`](@ref) defined as `η = √μ`
"""
mutable struct SqrtLink <: Link end

"""
    linkfun(L::Link, μ)

Return `η`, the value of the linear predictor for link `L` at mean `μ`.

# Examples
```jldoctest
julia> μ = inv(10):inv(5):1
0.1:0.2:0.9

julia> show(linkfun.(LogitLink(), μ))
[-2.19722, -0.847298, 0.0, 0.847298, 2.19722]

```
"""
function linkfun end

"""
    linkinv(L::Link, η)

Return `μ`, the mean value, for link `L` at linear predictor value `η`.

# Examples
```jldoctest
julia> μ = 0.1:0.2:1
0.1:0.2:0.9

julia> η = logit.(μ);

julia> linkinv.(LogitLink(), η) ≈ μ
true
```
"""
function linkinv end

"""
    mueta(L::Link, η)

Return the derivative of [`linkinv`](@ref), `dμ/dη`, for link `L` at linear predictor value `η`.

# Examples
```jldoctest
julia> mueta(LogitLink(), 0.0)
0.25

julia> mueta(CloglogLink(), 0.0) ≈ 0.36787944117144233
true

julia> mueta(LogLink(), 2.0) ≈ 7.38905609893065
true
```
"""
function mueta end

"""
    inverselink(L::Link, η)

Return a 3-tuple of the inverse link, the derivative of the inverse link, and when appropriate, the variance function `μ*(1 - μ)`.

The variance function is returned as NaN unless the range of μ is (0, 1)

# Examples
```jldoctest
julia> inverselink(LogitLink(), 0.0)
(0.5, 0.25, 0.25)

julia> μ, oneminusμ, variance = inverselink(CloglogLink(), 0.0);

julia> μ + oneminusμ ≈ 1
true

julia> μ*(1 - μ) ≈ variance
true

julia> isnan(last(inverselink(LogLink(), 2.0)))
true
```
"""
function inverselink end

"""
    canonicallink(D::Distribution)

Return the canonical link for distribution `D`, which must be in the exponential family.

# Examples
```jldoctest
julia> canonicallink(Bernoulli())
LogitLink()
```
"""
function canonicallink end

linkfun(::CauchitLink, μ) = tan(pi * (μ - oftype(μ, 1/2)))
linkinv(::CauchitLink, η) = oftype(η, 1/2) + atan(η) / pi
mueta(::CauchitLink, η) = one(η) / (pi * (one(η) + abs2(η)))
function inverselink(::CauchitLink, η)
    μlower = atan(-abs(η)) / π
    μlower += oftype(μlower, 1/2)
    η > 0 ? 1 - μlower : μlower, inv(π * (1 + abs2(η))), μlower * (1 - μlower)
end

linkfun(::CloglogLink, μ) = log(-log1p(-μ))
function linkinv(::CloglogLink, η::T) where T<:Real
    clamp(-expm1(-exp(η)), eps(T), one(T) - eps(T))
end
function mueta(::CloglogLink, η::T) where T<:Real
    max(eps(T), exp(η) * exp(-exp(η)))
end
function inverselink(::CloglogLink, η)
    expη = exp(η)
    μ = -expm1(-expη)
    omμ = exp(-expη)   # the complement, 1 - μ
    μ, max(floatmin(μ), expη * omμ), max(floatmin(μ), μ * omμ)
end

linkfun(::IdentityLink, μ) = μ
linkinv(::IdentityLink, η) = η
mueta(::IdentityLink, η) = one(η)
inverselink(::IdentityLink, η) = η, one(η), oftype(η, NaN)

linkfun(::InverseLink, μ) = inv(μ)
linkinv(::InverseLink, η) = inv(η)
mueta(::InverseLink, η) = -inv(abs2(η))
function inverselink(::InverseLink, η)
    μ = inv(η)
    μ, -abs2(μ), oftype(μ, NaN)
end

linkfun(::InverseSquareLink, μ) = inv(abs2(μ))
linkinv(::InverseSquareLink, η) = inv(sqrt(η))
mueta(::InverseSquareLink, η) = -inv(2η*sqrt(η))
function inverselink(::InverseSquareLink, η)
    μ = inv(sqrt(η))
    μ, -μ / (2η), oftype(μ, NaN)
end

linkfun(::LogitLink, μ) = logit(μ)
linkinv(::LogitLink, η) = logistic(η)
function mueta(::LogitLink, η)
    expabs = exp(-abs(η))
    denom = 1 + expabs
    (expabs / denom) / denom
end
function inverselink(::LogitLink, η)
    expabs = exp(-abs(η))
    opexpabs = 1 + expabs
    deriv = (expabs / opexpabs) / opexpabs
    η ≤ 0 ? expabs / opexpabs : inv(opexpabs), deriv, deriv
end

linkfun(::LogLink, μ) = log(μ)
linkinv(::LogLink, η) = exp(η)
mueta(::LogLink, η) = exp(η)
function inverselink(::LogLink, η)
    μ = exp(η)
    μ, μ, oftype(μ, NaN)
end

linkfun(nbl::NegativeBinomialLink, μ) = log(μ / (μ + nbl.θ))
linkinv(nbl::NegativeBinomialLink, η) = ℯ^η * nbl.θ / (1-ℯ^η)
mueta(nbl::NegativeBinomialLink, η) = ℯ^η * nbl.θ / (1-ℯ^η)
function inverselink(nbl::NegativeBinomialLink, η)
    μ = ℯ^η * nbl.θ / (1-ℯ^η)
    deriv = μ * (1 + μ / nbl.θ)
    μ, deriv, oftype(μ, NaN)
end

linkfun(::ProbitLink, μ) = -sqrt2 * erfcinv(2μ)
linkinv(::ProbitLink, η) = erfc(-η / sqrt2) / 2
mueta(::ProbitLink, η) = exp(-abs2(η) / 2) / sqrt2π
function inverselink(::ProbitLink, η)
    μlower = erfc(abs(η) / sqrt2) / 2
    μupper = 1 - μlower
    η < 0 ? μlower : μupper, exp(-abs2(η) / 2 ) / sqrt2π, μlower * μupper
end

linkfun(::SqrtLink, μ) = sqrt(μ)
linkinv(::SqrtLink, η) = abs2(η)
mueta(::SqrtLink, η) = 2η
inverselink(::SqrtLink, η) = abs2(η), 2η, oftype(η, NaN)

canonicallink(::Bernoulli) = LogitLink()
canonicallink(::Binomial) = LogitLink()
canonicallink(::Gamma) = InverseLink()
canonicallink(::InverseGaussian) = InverseSquareLink()
canonicallink(d::NegativeBinomial) = NegativeBinomialLink(d.r)
canonicallink(::Normal) = IdentityLink()
canonicallink(::Poisson) = LogLink()

"""
    glmvar(D::Distribution, μ)

Return the value of the variance function for `D` at `μ`

The variance of `D` at `μ` is the product of the dispersion parameter, ϕ, which does not
depend on `μ` and the value of `glmvar`.  In other words `glmvar` returns the factor of the
variance that depends on `μ`.

# Examples
```jldoctest
julia> μ = 1/6:1/3:1;

julia> glmvar.(Normal(), μ)    # constant for Normal()
3-element Array{Float64,1}:
 1.0
 1.0
 1.0

julia> glmvar.(Bernoulli(), μ) ≈ μ .* (1 .- μ)
true

julia> glmvar.(Poisson(), μ) == μ
true
```
"""
function glmvar end

glmvar(::Union{Bernoulli,Binomial}, μ) = μ * (1 - μ)
glmvar(::Gamma, μ) = abs2(μ)
glmvar(::InverseGaussian, μ) = μ^3
glmvar(d::NegativeBinomial, μ) = μ * (1 + μ/d.r)
glmvar(::Normal, μ) = one(μ)
glmvar(::Poisson, μ) = μ

"""
    mustart(D::Distribution, y, wt)

Return a starting value for μ.

For some distributions it is appropriate to set `μ = y` to initialize the IRLS algorithm but
for others, notably the Bernoulli, the values of `y` are not allowed as values of `μ` and
must be modified.

# Examples
```jldoctest
julia> mustart(Bernoulli(), 0.0, 1) ≈ 1/4
true

julia> mustart(Bernoulli(), 1.0, 1) ≈ 3/4
true

julia> mustart(Binomial(), 0.0, 10) ≈ 1/22
true

julia> mustart(Normal(), 0.0, 1) ≈ 0
true
```
"""
function mustart end

mustart(::Bernoulli, y, wt) = (y + oftype(y, 1/2)) / 2
mustart(::Binomial, y, wt) = (wt * y + oftype(y, 1/2)) / (wt + one(y))
function mustart(::Union{Gamma, InverseGaussian}, y, wt)
    fy = float(y)
    iszero(y) ? oftype(y, 1/10) : fy
end
function mustart(::NegativeBinomial, y, wt)
    fy = float(y)
    iszero(y) ? fy + oftype(fy, 1/6) : fy
end
mustart(::Normal, y, wt) = y
function mustart(::Poisson, y, wt)
    fy = float(y)
    fy + oftype(fy, 1/10)
end

"""
    devresid(D, y, μ)

Return the squared deviance residual of `μ` from `y` for distribution `D`

The deviance of a GLM can be evaluated as the sum of the squared deviance residuals.  This
is the principal use for these values.  The actual deviance residual, say for plotting, is
the signed square root of this value
```julia
sign(y - μ) * sqrt(devresid(D, y, μ))
```

# Examples
```jldoctest
julia> devresid(Normal(), 0, 0.25) ≈ abs2(0.25)
true

julia> devresid(Bernoulli(), 1, 0.75) ≈ -2*log(0.75)
true

julia> devresid(Bernoulli(), 0, 0.25) ≈ -2*log1p(-0.25)
true
```
"""
function devresid end

function devresid(::Bernoulli, y, μ)
    if y == 1
        return -2 * log(μ)
    elseif y == 0
        return -2 * log1p(-μ)
    end
    throw(ArgumentError("y should be 0 or 1 (got $y)"))
end
function devresid(::Binomial, y, μ)
    if y == 1
        return -2 * log(μ)
    elseif y == 0
        return -2 * log1p(-μ)
    else
        return 2 * (y * (log(y) - log(μ)) + (1 - y)*(log1p(-y) - log1p(-μ)))
    end
end
devresid(::Gamma, y, μ) = -2 * (log(y / μ) - (y - μ) / μ)
devresid(::InverseGaussian, y, μ) = abs2(y - μ) / (y * abs2(μ))
function devresid(d::NegativeBinomial, y, μ)
    θ = d.r
    v = 2 * (xlogy(y, y / μ) + xlogy(y + θ, (μ + θ)/(y + θ)))
    return μ == 0 ? oftype(v, NaN) : v
end
devresid(::Normal, y, μ) = abs2(y - μ)
devresid(::Poisson, y, μ) = 2 * (xlogy(y, y / μ) - (y - μ))

"""
    dispersion_parameter(D)  # not exported

Does distribution `D` have a separate dispersion parameter, ϕ?

Returns `false` for the `Bernoulli`, `Binomial` and `Poisson` distributions, `true` otherwise.

# Examples
```jldoctest
julia> show(GLM.dispersion_parameter(Normal()))
true
julia> show(GLM.dispersion_parameter(Bernoulli()))
false
```
"""
dispersion_parameter(D) = true
dispersion_parameter(::Union{Bernoulli, Binomial, Poisson}) = false

"""
    loglik_obs(D, y, μ, wt, ϕ)  # not exported

Returns `wt * logpdf(D(μ, ϕ), y)` where the parameters of `D` are derived from `μ` and `ϕ`.

The `wt` argument is a multiplier of the result except in the case of the `Binomial` where
`wt` is the number of trials and `μ` is the proportion of successes.

The loglikelihood of a fitted model is the sum of these values over all the observations.
"""
function loglik_obs end

loglik_obs(::Bernoulli, y, μ, wt, ϕ) = wt*logpdf(Bernoulli(μ), y)
loglik_obs(::Binomial, y, μ, wt, ϕ) = logpdf(Binomial(Int(wt), μ), Int(y*wt))
loglik_obs(::Gamma, y, μ, wt, ϕ) = wt*logpdf(Gamma(inv(ϕ), μ*ϕ), y)
loglik_obs(::InverseGaussian, y, μ, wt, ϕ) = wt*logpdf(InverseGaussian(μ, inv(ϕ)), y)
loglik_obs(::Normal, y, μ, wt, ϕ) = wt*logpdf(Normal(μ, sqrt(ϕ)), y)
loglik_obs(::Poisson, y, μ, wt, ϕ) = wt*logpdf(Poisson(μ), y)
# We use the following parameterization for the Negative Binomial distribution:
#    (Γ(θ+y) / (Γ(θ) * y!)) * μ^y * θ^θ / (μ+θ)^{θ+y}
# The parameterization of NegativeBinomial(r=θ, p) in Distributions.jl is
#    Γ(θ+y) / (y! * Γ(θ)) * p^θ(1-p)^y
# Hence, p = θ/(μ+θ)
loglik_obs(d::NegativeBinomial, y, μ, wt, ϕ) = wt*logpdf(NegativeBinomial(d.r, d.r/(μ+d.r)), y)


text:
# Autogenerated wrapper script for aiger_jll for aarch64-linux-gnu
export aiginfo, aigsplit, aigtoblif, aigmiter, bliftoaig, aigand, aigfuzz, aigtocnf, aigtodot, smvtoaig, wrapstim, aigreset, aigstrip, aigtosmv, aigmove, aigsim, andtoaig, soltostim, aigunroll, aigtoaig, aigunconstraint, aigdd, aigflip, aigor, aignm, aigjoin

## Global variables
PATH = ""
LIBPATH = ""
LIBPATH_env = "LD_LIBRARY_PATH"

# Relative path to `aiginfo`
const aiginfo_splitpath = ["bin", "aiginfo"]

# This will be filled out by __init__() for all products, as it must be done at runtime
aiginfo_path = ""

# aiginfo-specific global declaration
function aiginfo(f::Function; adjust_PATH::Bool = true, adjust_LIBPATH::Bool = true)
    global PATH, LIBPATH
    env_mapping = Dict{String,String}()
    if adjust_PATH
        if !isempty(get(ENV, "PATH", ""))
            env_mapping["PATH"] = string(PATH, ':', ENV["PATH"])
        else
            env_mapping["PATH"] = PATH
        end
    end
    if adjust_LIBPATH
        if !isempty(get(ENV, LIBPATH_env, ""))
            env_mapping[LIBPATH_env] = string(LIBPATH, ':', ENV[LIBPATH_env])
        else
            env_mapping[LIBPATH_env] = LIBPATH
        end
    end
    withenv(env_mapping...) do
        f(aiginfo_path)
    end
end

# Relative path to `aigsplit`
const aigsplit_splitpath = ["bin", "aigsplit"]

# This will be filled out by __init__() for all products, as it must be done at runtime
aigsplit_path = ""

# aigsplit-specific global declaration
function aigsplit(f::Function; adjust_PATH::Bool = true, adjust_LIBPATH::Bool = true)
    global PATH, LIBPATH
    env_mapping = Dict{String,String}()
    if adjust_PATH
        if !isempty(get(ENV, "PATH", ""))
            env_mapping["PATH"] = string(PATH, ':', ENV["PATH"])
        else
            env_mapping["PATH"] = PATH
        end
    end
    if adjust_LIBPATH
        if !isempty(get(ENV, LIBPATH_env, ""))
            env_mapping[LIBPATH_env] = string(LIBPATH, ':', ENV[LIBPATH_env])
        else
            env_mapping[LIBPATH_env] = LIBPATH
        end
    end
    withenv(env_mapping...) do
        f(aigsplit_path)
    end
end

# Relative path to `aigtoblif`
const aigtoblif_splitpath = ["bin", "aigtoblif"]

# This will be filled out by __init__() for all products, as it must be done at runtime
aigtoblif_path = ""

# aigtoblif-specific global declaration
function aigtoblif(f::Function; adjust_PATH::Bool = true, adjust_LIBPATH::Bool = true)
    global PATH, LIBPATH
    env_mapping = Dict{String,String}()
    if adjust_PATH
        if !isempty(get(ENV, "PATH", ""))
            env_mapping["PATH"] = string(PATH, ':', ENV["PATH"])
        else
            env_mapping["PATH"] = PATH
        end
    end
    if adjust_LIBPATH
        if !isempty(get(ENV, LIBPATH_env, ""))
            env_mapping[LIBPATH_env] = string(LIBPATH, ':', ENV[LIBPATH_env])
        else
            env_mapping[LIBPATH_env] = LIBPATH
        end
    end
    withenv(env_mapping...) do
        f(aigtoblif_path)
    end
end

# Relative path to `aigmiter`
const aigmiter_splitpath = ["bin", "aigmiter"]

# This will be filled out by __init__() for all products, as it must be done at runtime
aigmiter_path = ""

# aigmiter-specific global declaration
function aigmiter(f::Function; adjust_PATH::Bool = true, adjust_LIBPATH::Bool = true)
    global PATH, LIBPATH
    env_mapping = Dict{String,String}()
    if adjust_PATH
        if !isempty(get(ENV, "PATH", ""))
            env_mapping["PATH"] = string(PATH, ':', ENV["PATH"])
        else
            env_mapping["PATH"] = PATH
        end
    end
    if adjust_LIBPATH
        if !isempty(get(ENV, LIBPATH_env, ""))
            env_mapping[LIBPATH_env] = string(LIBPATH, ':', ENV[LIBPATH_env])
        else
            env_mapping[LIBPATH_env] = LIBPATH
        end
    end
    withenv(env_mapping...) do
        f(aigmiter_path)
    end
end

# Relative path to `bliftoaig`
const bliftoaig_splitpath = ["bin", "bliftoaig"]

# This will be filled out by __init__() for all products, as it must be done at runtime
bliftoaig_path = ""

# bliftoaig-specific global declaration
function bliftoaig(f::Function; adjust_PATH::Bool = true, adjust_LIBPATH::Bool = true)
    global PATH, LIBPATH
    env_mapping = Dict{String,String}()
    if adjust_PATH
        if !isempty(get(ENV, "PATH", ""))
            env_mapping["PATH"] = string(PATH, ':', ENV["PATH"])
        else
            env_mapping["PATH"] = PATH
        end
    end
    if adjust_LIBPATH
        if !isempty(get(ENV, LIBPATH_env, ""))
            env_mapping[LIBPATH_env] = string(LIBPATH, ':', ENV[LIBPATH_env])
        else
            env_mapping[LIBPATH_env] = LIBPATH
        end
    end
    withenv(env_mapping...) do
        f(bliftoaig_path)
    end
end

# Relative path to `aigand`
const aigand_splitpath = ["bin", "aigand"]

# This will be filled out by __init__() for all products, as it must be done at runtime
aigand_path = ""

# aigand-specific global declaration
function aigand(f::Function; adjust_PATH::Bool = true, adjust_LIBPATH::Bool = true)
    global PATH, LIBPATH
    env_mapping = Dict{String,String}()
    if adjust_PATH
        if !isempty(get(ENV, "PATH", ""))
            env_mapping["PATH"] = string(PATH, ':', ENV["PATH"])
        else
            env_mapping["PATH"] = PATH
        end
    end
    if adjust_LIBPATH
        if !isempty(get(ENV, LIBPATH_env, ""))
            env_mapping[LIBPATH_env] = string(LIBPATH, ':', ENV[LIBPATH_env])
        else
            env_mapping[LIBPATH_env] = LIBPATH
        end
    end
    withenv(env_mapping...) do
        f(aigand_path)
    end
end

# Relative path to `aigfuzz`
const aigfuzz_splitpath = ["bin", "aigfuzz"]

# This will be filled out by __init__() for all products, as it must be done at runtime
aigfuzz_path = ""

# aigfuzz-specific global declaration
function aigfuzz(f::Function; adjust_PATH::Bool = true, adjust_LIBPATH::Bool = true)
    global PATH, LIBPATH
    env_mapping = Dict{String,String}()
    if adjust_PATH
        if !isempty(get(ENV, "PATH", ""))
            env_mapping["PATH"] = string(PATH, ':', ENV["PATH"])
        else
            env_mapping["PATH"] = PATH
        end
    end
    if adjust_LIBPATH
        if !isempty(get(ENV, LIBPATH_env, ""))
            env_mapping[LIBPATH_env] = string(LIBPATH, ':', ENV[LIBPATH_env])
        else
            env_mapping[LIBPATH_env] = LIBPATH
        end
    end
    withenv(env_mapping...) do
        f(aigfuzz_path)
    end
end

# Relative path to `aigtocnf`
const aigtocnf_splitpath = ["bin", "aigtocnf"]

# This will be filled out by __init__() for all products, as it must be done at runtime
aigtocnf_path = ""

# aigtocnf-specific global declaration
function aigtocnf(f::Function; adjust_PATH::Bool = true, adjust_LIBPATH::Bool = true)
    global PATH, LIBPATH
    env_mapping = Dict{String,String}()
    if adjust_PATH
        if !isempty(get(ENV, "PATH", ""))
            env_mapping["PATH"] = string(PATH, ':', ENV["PATH"])
        else
            env_mapping["PATH"] = PATH
        end
    end
    if adjust_LIBPATH
        if !isempty(get(ENV, LIBPATH_env, ""))
            env_mapping[LIBPATH_env] = string(LIBPATH, ':', ENV[LIBPATH_env])
        else
            env_mapping[LIBPATH_env] = LIBPATH
        end
    end
    withenv(env_mapping...) do
        f(aigtocnf_path)
    end
end

# Relative path to `aigtodot`
const aigtodot_splitpath = ["bin", "aigtodot"]

# This will be filled out by __init__() for all products, as it must be done at runtime
aigtodot_path = ""

# aigtodot-specific global declaration
function aigtodot(f::Function; adjust_PATH::Bool = true, adjust_LIBPATH::Bool = true)
    global PATH, LIBPATH
    env_mapping = Dict{String,String}()
    if adjust_PATH
        if !isempty(get(ENV, "PATH", ""))
            env_mapping["PATH"] = string(PATH, ':', ENV["PATH"])
        else
            env_mapping["PATH"] = PATH
        end
    end
    if adjust_LIBPATH
        if !isempty(get(ENV, LIBPATH_env, ""))
            env_mapping[LIBPATH_env] = string(LIBPATH, ':', ENV[LIBPATH_env])
        else
            env_mapping[LIBPATH_env] = LIBPATH
        end
    end
    withenv(env_mapping...) do
        f(aigtodot_path)
    end
end

# Relative path to `smvtoaig`
const smvtoaig_splitpath = ["bin", "smvtoaig"]

# This will be filled out by __init__() for all products, as it must be done at runtime
smvtoaig_path = ""

# smvtoaig-specific global declaration
function smvtoaig(f::Function; adjust_PATH::Bool = true, adjust_LIBPATH::Bool = true)
    global PATH, LIBPATH
    env_mapping = Dict{String,String}()
    if adjust_PATH
        if !isempty(get(ENV, "PATH", ""))
            env_mapping["PATH"] = string(PATH, ':', ENV["PATH"])
        else
            env_mapping["PATH"] = PATH
        end
    end
    if adjust_LIBPATH
        if !isempty(get(ENV, LIBPATH_env, ""))
            env_mapping[LIBPATH_env] = string(LIBPATH, ':', ENV[LIBPATH_env])
        else
            env_mapping[LIBPATH_env] = LIBPATH
        end
    end
    withenv(env_mapping...) do
        f(smvtoaig_path)
    end
end

# Relative path to `wrapstim`
const wrapstim_splitpath = ["bin", "wrapstim"]

# This will be filled out by __init__() for all products, as it must be done at runtime
wrapstim_path = ""

# wrapstim-specific global declaration
function wrapstim(f::Function; adjust_PATH::Bool = true, adjust_LIBPATH::Bool = true)
    global PATH, LIBPATH
    env_mapping = Dict{String,String}()
    if adjust_PATH
        if !isempty(get(ENV, "PATH", ""))
            env_mapping["PATH"] = string(PATH, ':', ENV["PATH"])
        else
            env_mapping["PATH"] = PATH
        end
    end
    if adjust_LIBPATH
        if !isempty(get(ENV, LIBPATH_env, ""))
            env_mapping[LIBPATH_env] = string(LIBPATH, ':', ENV[LIBPATH_env])
        else
            env_mapping[LIBPATH_env] = LIBPATH
        end
    end
    withenv(env_mapping...) do
        f(wrapstim_path)
    end
end

# Relative path to `aigreset`
const aigreset_splitpath = ["bin", "aigreset"]

# This will be filled out by __init__() for all products, as it must be done at runtime
aigreset_path = ""

# aigreset-specific global declaration
function aigreset(f::Function; adjust_PATH::Bool = true, adjust_LIBPATH::Bool = true)
    global PATH, LIBPATH
    env_mapping = Dict{String,String}()
    if adjust_PATH
        if !isempty(get(ENV, "PATH", ""))
            env_mapping["PATH"] = string(PATH, ':', ENV["PATH"])
        else
            env_mapping["PATH"] = PATH
        end
    end
    if adjust_LIBPATH
        if !isempty(get(ENV, LIBPATH_env, ""))
            env_mapping[LIBPATH_env] = string(LIBPATH, ':', ENV[LIBPATH_env])
        else
            env_mapping[LIBPATH_env] = LIBPATH
        end
    end
    withenv(env_mapping...) do
        f(aigreset_path)
    end
end

# Relative path to `aigstrip`
const aigstrip_splitpath = ["bin", "aigstrip"]

# This will be filled out by __init__() for all products, as it must be done at runtime
aigstrip_path = ""

# aigstrip-specific global declaration
function aigstrip(f::Function; adjust_PATH::Bool = true, adjust_LIBPATH::Bool = true)
    global PATH, LIBPATH
    env_mapping = Dict{String,String}()
    if adjust_PATH
        if !isempty(get(ENV, "PATH", ""))
            env_mapping["PATH"] = string(PATH, ':', ENV["PATH"])
        else
            env_mapping["PATH"] = PATH
        end
    end
    if adjust_LIBPATH
        if !isempty(get(ENV, LIBPATH_env, ""))
            env_mapping[LIBPATH_env] = string(LIBPATH, ':', ENV[LIBPATH_env])
        else
            env_mapping[LIBPATH_env] = LIBPATH
        end
    end
    withenv(env_mapping...) do
        f(aigstrip_path)
    end
end

# Relative path to `aigtosmv`
const aigtosmv_splitpath = ["bin", "aigtosmv"]

# This will be filled out by __init__() for all products, as it must be done at runtime
aigtosmv_path = ""

# aigtosmv-specific global declaration
function aigtosmv(f::Function; adjust_PATH::Bool = true, adjust_LIBPATH::Bool = true)
    global PATH, LIBPATH
    env_mapping = Dict{String,String}()
    if adjust_PATH
        if !isempty(get(ENV, "PATH", ""))
            env_mapping["PATH"] = string(PATH, ':', ENV["PATH"])
        else
            env_mapping["PATH"] = PATH
        end
    end
    if adjust_LIBPATH
        if !isempty(get(ENV, LIBPATH_env, ""))
            env_mapping[LIBPATH_env] = string(LIBPATH, ':', ENV[LIBPATH_env])
        else
            env_mapping[LIBPATH_env] = LIBPATH
        end
    end
    withenv(env_mapping...) do
        f(aigtosmv_path)
    end
end

# Relative path to `aigmove`
const aigmove_splitpath = ["bin", "aigmove"]

# This will be filled out by __init__() for all products, as it must be done at runtime
aigmove_path = ""

# aigmove-specific global declaration
function aigmove(f::Function; adjust_PATH::Bool = true, adjust_LIBPATH::Bool = true)
    global PATH, LIBPATH
    env_mapping = Dict{String,String}()
    if adjust_PATH
        if !isempty(get(ENV, "PATH", ""))
            env_mapping["PATH"] = string(PATH, ':', ENV["PATH"])
        else
            env_mapping["PATH"] = PATH
        end
    end
    if adjust_LIBPATH
        if !isempty(get(ENV, LIBPATH_env, ""))
            env_mapping[LIBPATH_env] = string(LIBPATH, ':', ENV[LIBPATH_env])
        else
            env_mapping[LIBPATH_env] = LIBPATH
        end
    end
    withenv(env_mapping...) do
        f(aigmove_path)
    end
end

# Relative path to `aigsim`
const aigsim_splitpath = ["bin", "aigsim"]

# This will be filled out by __init__() for all products, as it must be done at runtime
aigsim_path = ""

# aigsim-specific global declaration
function aigsim(f::Function; adjust_PATH::Bool = true, adjust_LIBPATH::Bool = true)
    global PATH, LIBPATH
    env_mapping = Dict{String,String}()
    if adjust_PATH
        if !isempty(get(ENV, "PATH", ""))
            env_mapping["PATH"] = string(PATH, ':', ENV["PATH"])
        else
            env_mapping["PATH"] = PATH
        end
    end
    if adjust_LIBPATH
        if !isempty(get(ENV, LIBPATH_env, ""))
            env_mapping[LIBPATH_env] = string(LIBPATH, ':', ENV[LIBPATH_env])
        else
            env_mapping[LIBPATH_env] = LIBPATH
        end
    end
    withenv(env_mapping...) do
        f(aigsim_path)
    end
end

# Relative path to `andtoaig`
const andtoaig_splitpath = ["bin", "andtoaig"]

# This will be filled out by __init__() for all products, as it must be done at runtime
andtoaig_path = ""

# andtoaig-specific global declaration
function andtoaig(f::Function; adjust_PATH::Bool = true, adjust_LIBPATH::Bool = true)
    global PATH, LIBPATH
    env_mapping = Dict{String,String}()
    if adjust_PATH
        if !isempty(get(ENV, "PATH", ""))
            env_mapping["PATH"] = string(PATH, ':', ENV["PATH"])
        else
            env_mapping["PATH"] = PATH
        end
    end
    if adjust_LIBPATH
        if !isempty(get(ENV, LIBPATH_env, ""))
            env_mapping[LIBPATH_env] = string(LIBPATH, ':', ENV[LIBPATH_env])
        else
            env_mapping[LIBPATH_env] = LIBPATH
        end
    end
    withenv(env_mapping...) do
        f(andtoaig_path)
    end
end

# Relative path to `soltostim`
const soltostim_splitpath = ["bin", "soltostim"]

# This will be filled out by __init__() for all products, as it must be done at runtime
soltostim_path = ""

# soltostim-specific global declaration
function soltostim(f::Function; adjust_PATH::Bool = true, adjust_LIBPATH::Bool = true)
    global PATH, LIBPATH
    env_mapping = Dict{String,String}()
    if adjust_PATH
        if !isempty(get(ENV, "PATH", ""))
            env_mapping["PATH"] = string(PATH, ':', ENV["PATH"])
        else
            env_mapping["PATH"] = PATH
        end
    end
    if adjust_LIBPATH
        if !isempty(get(ENV, LIBPATH_env, ""))
            env_mapping[LIBPATH_env] = string(LIBPATH, ':', ENV[LIBPATH_env])
        else
            env_mapping[LIBPATH_env] = LIBPATH
        end
    end
    withenv(env_mapping...) do
        f(soltostim_path)
    end
end

# Relative path to `aigunroll`
const aigunroll_splitpath = ["bin", "aigunroll"]

# This will be filled out by __init__() for all products, as it must be done at runtime
aigunroll_path = ""

# aigunroll-specific global declaration
function aigunroll(f::Function; adjust_PATH::Bool = true, adjust_LIBPATH::Bool = true)
    global PATH, LIBPATH
    env_mapping = Dict{String,String}()
    if adjust_PATH
        if !isempty(get(ENV, "PATH", ""))
            env_mapping["PATH"] = string(PATH, ':', ENV["PATH"])
        else
            env_mapping["PATH"] = PATH
        end
    end
    if adjust_LIBPATH
        if !isempty(get(ENV, LIBPATH_env, ""))
            env_mapping[LIBPATH_env] = string(LIBPATH, ':', ENV[LIBPATH_env])
        else
            env_mapping[LIBPATH_env] = LIBPATH
        end
    end
    withenv(env_mapping...) do
        f(aigunroll_path)
    end
end

# Relative path to `aigtoaig`
const aigtoaig_splitpath = ["bin", "aigtoaig"]

# This will be filled out by __init__() for all products, as it must be done at runtime
aigtoaig_path = ""

# aigtoaig-specific global declaration
function aigtoaig(f::Function; adjust_PATH::Bool = true, adjust_LIBPATH::Bool = true)
    global PATH, LIBPATH
    env_mapping = Dict{String,String}()
    if adjust_PATH
        if !isempty(get(ENV, "PATH", ""))
            env_mapping["PATH"] = string(PATH, ':', ENV["PATH"])
        else
            env_mapping["PATH"] = PATH
        end
    end
    if adjust_LIBPATH
        if !isempty(get(ENV, LIBPATH_env, ""))
            env_mapping[LIBPATH_env] = string(LIBPATH, ':', ENV[LIBPATH_env])
        else
            env_mapping[LIBPATH_env] = LIBPATH
        end
    end
    withenv(env_mapping...) do
        f(aigtoaig_path)
    end
end

# Relative path to `aigunconstraint`
const aigunconstraint_splitpath = ["bin", "aigunconstraint"]

# This will be filled out by __init__() for all products, as it must be done at runtime
aigunconstraint_path = ""

# aigunconstraint-specific global declaration
function aigunconstraint(f::Function; adjust_PATH::Bool = true, adjust_LIBPATH::Bool = true)
    global PATH, LIBPATH
    env_mapping = Dict{String,String}()
    if adjust_PATH
        if !isempty(get(ENV, "PATH", ""))
            env_mapping["PATH"] = string(PATH, ':', ENV["PATH"])
        else
            env_mapping["PATH"] = PATH
        end
    end
    if adjust_LIBPATH
        if !isempty(get(ENV, LIBPATH_env, ""))
            env_mapping[LIBPATH_env] = string(LIBPATH, ':', ENV[LIBPATH_env])
        else
            env_mapping[LIBPATH_env] = LIBPATH
        end
    end
    withenv(env_mapping...) do
        f(aigunconstraint_path)
    end
end

# Relative path to `aigdd`
const aigdd_splitpath = ["bin", "aigdd"]

# This will be filled out by __init__() for all products, as it must be done at runtime
aigdd_path = ""

# aigdd-specific global declaration
function aigdd(f::Function; adjust_PATH::Bool = true, adjust_LIBPATH::Bool = true)
    global PATH, LIBPATH
    env_mapping = Dict{String,String}()
    if adjust_PATH
        if !isempty(get(ENV, "PATH", ""))
            env_mapping["PATH"] = string(PATH, ':', ENV["PATH"])
        else
            env_mapping["PATH"] = PATH
        end
    end
    if adjust_LIBPATH
        if !isempty(get(ENV, LIBPATH_env, ""))
            env_mapping[LIBPATH_env] = string(LIBPATH, ':', ENV[LIBPATH_env])
        else
            env_mapping[LIBPATH_env] = LIBPATH
        end
    end
    withenv(env_mapping...) do
        f(aigdd_path)
    end
end

# Relative path to `aigflip`
const aigflip_splitpath = ["bin", "aigflip"]

# This will be filled out by __init__() for all products, as it must be done at runtime
aigflip_path = ""

# aigflip-specific global declaration
function aigflip(f::Function; adjust_PATH::Bool = true, adjust_LIBPATH::Bool = true)
    global PATH, LIBPATH
    env_mapping = Dict{String,String}()
    if adjust_PATH
        if !isempty(get(ENV, "PATH", ""))
            env_mapping["PATH"] = string(PATH, ':', ENV["PATH"])
        else
            env_mapping["PATH"] = PATH
        end
    end
    if adjust_LIBPATH
        if !isempty(get(ENV, LIBPATH_env, ""))
            env_mapping[LIBPATH_env] = string(LIBPATH, ':', ENV[LIBPATH_env])
        else
            env_mapping[LIBPATH_env] = LIBPATH
        end
    end
    withenv(env_mapping...) do
        f(aigflip_path)
    end
end

# Relative path to `aigor`
const aigor_splitpath = ["bin", "aigor"]

# This will be filled out by __init__() for all products, as it must be done at runtime
aigor_path = ""

# aigor-specific global declaration
function aigor(f::Function; adjust_PATH::Bool = true, adjust_LIBPATH::Bool = true)
    global PATH, LIBPATH
    env_mapping = Dict{String,String}()
    if adjust_PATH
        if !isempty(get(ENV, "PATH", ""))
            env_mapping["PATH"] = string(PATH, ':', ENV["PATH"])
        else
            env_mapping["PATH"] = PATH
        end
    end
    if adjust_LIBPATH
        if !isempty(get(ENV, LIBPATH_env, ""))
            env_mapping[LIBPATH_env] = string(LIBPATH, ':', ENV[LIBPATH_env])
        else
            env_mapping[LIBPATH_env] = LIBPATH
        end
    end
    withenv(env_mapping...) do
        f(aigor_path)
    end
end

# Relative path to `aignm`
const aignm_splitpath = ["bin", "aignm"]

# This will be filled out by __init__() for all products, as it must be done at runtime
aignm_path = ""

# aignm-specific global declaration
function aignm(f::Function; adjust_PATH::Bool = true, adjust_LIBPATH::Bool = true)
    global PATH, LIBPATH
    env_mapping = Dict{String,String}()
    if adjust_PATH
        if !isempty(get(ENV, "PATH", ""))
            env_mapping["PATH"] = string(PATH, ':', ENV["PATH"])
        else
            env_mapping["PATH"] = PATH
        end
    end
    if adjust_LIBPATH
        if !isempty(get(ENV, LIBPATH_env, ""))
            env_mapping[LIBPATH_env] = string(LIBPATH, ':', ENV[LIBPATH_env])
        else
            env_mapping[LIBPATH_env] = LIBPATH
        end
    end
    withenv(env_mapping...) do
        f(aignm_path)
    end
end

# Relative path to `aigjoin`
const aigjoin_splitpath = ["bin", "aigjoin"]

# This will be filled out by __init__() for all products, as it must be done at runtime
aigjoin_path = ""

# aigjoin-specific global declaration
function aigjoin(f::Function; adjust_PATH::Bool = true, adjust_LIBPATH::Bool = true)
    global PATH, LIBPATH
    env_mapping = Dict{String,String}()
    if adjust_PATH
        if !isempty(get(ENV, "PATH", ""))
            env_mapping["PATH"] = string(PATH, ':', ENV["PATH"])
        else
            env_mapping["PATH"] = PATH
        end
    end
    if adjust_LIBPATH
        if !isempty(get(ENV, LIBPATH_env, ""))
            env_mapping[LIBPATH_env] = string(LIBPATH, ':', ENV[LIBPATH_env])
        else
            env_mapping[LIBPATH_env] = LIBPATH
        end
    end
    withenv(env_mapping...) do
        f(aigjoin_path)
    end
end

"""
Open all libraries
"""
function __init__()
    global artifact_dir = abspath(artifact"aiger")

    # Initialize PATH and LIBPATH environment variable listings
    global PATH_list, LIBPATH_list
    # We first need to add to LIBPATH_list the libraries provided by Julia
    append!(LIBPATH_list, [joinpath(Sys.BINDIR, Base.LIBDIR, "julia"), joinpath(Sys.BINDIR, Base.LIBDIR)])
    global aiginfo_path = normpath(joinpath(artifact_dir, aiginfo_splitpath...))

    push!(PATH_list, dirname(aiginfo_path))
    global aigsplit_path = normpath(joinpath(artifact_dir, aigsplit_splitpath...))

    push!(PATH_list, dirname(aigsplit_path))
    global aigtoblif_path = normpath(joinpath(artifact_dir, aigtoblif_splitpath...))

    push!(PATH_list, dirname(aigtoblif_path))
    global aigmiter_path = normpath(joinpath(artifact_dir, aigmiter_splitpath...))

    push!(PATH_list, dirname(aigmiter_path))
    global bliftoaig_path = normpath(joinpath(artifact_dir, bliftoaig_splitpath...))

    push!(PATH_list, dirname(bliftoaig_path))
    global aigand_path = normpath(joinpath(artifact_dir, aigand_splitpath...))

    push!(PATH_list, dirname(aigand_path))
    global aigfuzz_path = normpath(joinpath(artifact_dir, aigfuzz_splitpath...))

    push!(PATH_list, dirname(aigfuzz_path))
    global aigtocnf_path = normpath(joinpath(artifact_dir, aigtocnf_splitpath...))

    push!(PATH_list, dirname(aigtocnf_path))
    global aigtodot_path = normpath(joinpath(artifact_dir, aigtodot_splitpath...))

    push!(PATH_list, dirname(aigtodot_path))
    global smvtoaig_path = normpath(joinpath(artifact_dir, smvtoaig_splitpath...))

    push!(PATH_list, dirname(smvtoaig_path))
    global wrapstim_path = normpath(joinpath(artifact_dir, wrapstim_splitpath...))

    push!(PATH_list, dirname(wrapstim_path))
    global aigreset_path = normpath(joinpath(artifact_dir, aigreset_splitpath...))

    push!(PATH_list, dirname(aigreset_path))
    global aigstrip_path = normpath(joinpath(artifact_dir, aigstrip_splitpath...))

    push!(PATH_list, dirname(aigstrip_path))
    global aigtosmv_path = normpath(joinpath(artifact_dir, aigtosmv_splitpath...))

    push!(PATH_list, dirname(aigtosmv_path))
    global aigmove_path = normpath(joinpath(artifact_dir, aigmove_splitpath...))

    push!(PATH_list, dirname(aigmove_path))
    global aigsim_path = normpath(joinpath(artifact_dir, aigsim_splitpath...))

    push!(PATH_list, dirname(aigsim_path))
    global andtoaig_path = normpath(joinpath(artifact_dir, andtoaig_splitpath...))

    push!(PATH_list, dirname(andtoaig_path))
    global soltostim_path = normpath(joinpath(artifact_dir, soltostim_splitpath...))

    push!(PATH_list, dirname(soltostim_path))
    global aigunroll_path = normpath(joinpath(artifact_dir, aigunroll_splitpath...))

    push!(PATH_list, dirname(aigunroll_path))
    global aigtoaig_path = normpath(joinpath(artifact_dir, aigtoaig_splitpath...))

    push!(PATH_list, dirname(aigtoaig_path))
    global aigunconstraint_path = normpath(joinpath(artifact_dir, aigunconstraint_splitpath...))

    push!(PATH_list, dirname(aigunconstraint_path))
    global aigdd_path = normpath(joinpath(artifact_dir, aigdd_splitpath...))

    push!(PATH_list, dirname(aigdd_path))
    global aigflip_path = normpath(joinpath(artifact_dir, aigflip_splitpath...))

    push!(PATH_list, dirname(aigflip_path))
    global aigor_path = normpath(joinpath(artifact_dir, aigor_splitpath...))

    push!(PATH_list, dirname(aigor_path))
    global aignm_path = normpath(joinpath(artifact_dir, aignm_splitpath...))

    push!(PATH_list, dirname(aignm_path))
    global aigjoin_path = normpath(joinpath(artifact_dir, aigjoin_splitpath...))

    push!(PATH_list, dirname(aigjoin_path))
    # Filter out duplicate and empty entries in our PATH and LIBPATH entries
    filter!(!isempty, unique!(PATH_list))
    filter!(!isempty, unique!(LIBPATH_list))
    global PATH = join(PATH_list, ':')
    global LIBPATH = join(LIBPATH_list, ':')

    # Add each element of LIBPATH to our DL_LOAD_PATH (necessary on platforms
    # that don't honor our "already opened" trick)
    #for lp in LIBPATH_list
    #    push!(DL_LOAD_PATH, lp)
    #end
end  # __init__()


text:
module TestAstTools

using FLoops: unbound_rhs, is_dotcall
using Test

function test_unbound_rhs()
    @test unbound_rhs(:a) == [:a]
    @test unbound_rhs(:(
        let x = 1
            (x, y)
        end
    )) == [:y]
    @test unbound_rhs(quote end) == []
    @test unbound_rhs(quote
        ()
    end) == []
    @test unbound_rhs(quote
        []
    end) == []
    @test unbound_rhs(nothing) == []
end

function test_is_dotcall()
    @test is_dotcall(:(a .+ b))
    @test is_dotcall(:(.√a))
    @test is_dotcall(:(f.(x)))
    @test is_dotcall(:((f.a).(b)))
    @test !is_dotcall(:(a + b))
    @test !is_dotcall(:(√a))
    @test !is_dotcall(:(f(x)))
    @test !is_dotcall(:(f.a))
    @test !is_dotcall(:((f.a)(b)))
    @test !is_dotcall(:(f.a(b)))
end

function test_is_dotcall_nargs()
    @test is_dotcall(:(a .+ b), 2)
    @test is_dotcall(:(.√a), 1)
    @test is_dotcall(:(f.(x)), 1)
    @test is_dotcall(:((f.a).(b)), 1)
    @test !is_dotcall(:(a .+ b), 1)
    @test !is_dotcall(:(.√a), 2)
    @test !is_dotcall(:(f.(x)), 2)
    @test !is_dotcall(:((f.a).(b)), 2)
    @test !is_dotcall(:(a + b), 2)
    @test !is_dotcall(:(√a), 1)
    @test !is_dotcall(:(f(x)), 1)
    @test !is_dotcall(:(f.a), 1)
    @test !is_dotcall(:((f.a)(b)), 1)
    @test !is_dotcall(:(f.a(b)), 1)
end

end  # module


text:
using Test
using Random

include("yahtzee.jl")

@testset "test roll_dice" begin
    # fix random seed
    Random.seed!(1234);
    # known sequence for given seed
    exp_first_throw = [1, 3, 5, 1, 6]
    exp_second_throw = [1, 6, 4, 5, 4]
    exp_third_throw = [1, 2, 5, 3, 2]

    n_dice = 5
    sides=(1,2,3,4,5,6)
    throw = roll_dice(n_dice, sides)
    @test throw == exp_first_throw
    throw = roll_dice(n_dice, sides)
    @test throw == exp_second_throw
    throw = roll_dice(n_dice, sides)
    @test throw == exp_third_throw
end

@testset "testing yahtzee" begin
    # try a million throws and see if we get close
    # to the statistical average of 4.6%
    # (https://en.wikipedia.org/wiki/Yahtzee#Probabilities)
    n_yahtzees = 0
    n_tests = 1_000_000
    for i in 1:n_tests
        collected_side, n_collected = yahtzee()
        if n_collected == 5
            n_yahtzees += 1
        end
    end
    @test n_yahtzees/n_tests ≈ 0.046 atol=0.001
end


text:
"""
$(TYPEDEF)

$(TYPEDFIELDS)

A simple wrapper for a `DiscreteUniform` distribution over non-numerical arrays.

"""
struct CategoricalFactor <: DiscreteUnivariateDistribution
    values::Array{N, 1} where N
    distribution::DiscreteUniform
end

"""
$(TYPEDSIGNATURES)

```jldoctest
julia> a = CategoricalFactor([:a, :b, 2, 1.0])
CategoricalFactor(
values: Any[:a, :b, 2, 1.0]
distribution: Distributions.DiscreteUniform(a=1, b=4)
)

```
"""
function CategoricalFactor(values::Array{N, 1}) where N
    CategoricalFactor(values, DiscreteUniform(1, length(values)))
end

"""
$(TYPEDSIGNATURES)

```jldoctest
julia> rand(CategoricalFactor([:a, :b, 2, 1.0]), 6)
6-element Array{Any,1}:
  :b
 2
 1.0
 2
  :a
 2

```
"""
function rand(distribution::CategoricalFactor, n::Int = 1)
    [distribution.values[s] for s in rand(distribution.distribution, n)]
end


text:
# Functions for loading images from the Sloan Digital Sky Survey.
# Everything specific to SDSS should be confined to this file.

module SDSSIO

import Base.convert, Base.getindex

import FITSIO
import WCS
using CodecZlib
using CodecBzip2

import ..Celeste: SurveyDataSet, BoundingBox, load_images
import ..Model: Image, CatalogEntry, AbstractPSFMap
import ..PSF

export SDSSDataSet

# types of things to mask in read_mask().
const DEFAULT_MASK_PLANES = ["S_MASK_INTERP",  # bad pixel (was interpolated)
                             "S_MASK_SATUR",  # saturated
                             "S_MASK_CR",  # cosmic ray
                             "S_MASK_GHOST"]  # electronics artifacts

const BAND_CHAR_TO_NUM = Dict('u'=>1, 'g'=>2, 'r'=>3, 'i'=>4, 'z'=>5)

# -----------------------------------------------------------------------------
# Background / sky

"""
    SDSSBackground(sky_small, sky_x, sky_y, calibration) <: AbstractArray{Float32,2}

The 2-d background ("sky") in an SDSS image, supporting `size` and
`getindex`.  `sky_small` is a small version of the sky image in counts
(ADU).  `sky_x` and `sky_y` map indicies in the full image to
coordinates in `sky_small`, where bilinear interpolation is then used
to get values at non-integer coordinates. The interpolated value is
then multiplied by `calibration`.

For example, suppose `bkg` is a `SDSSBackground` instance. `bkg[i, j]`
is determined by first finding the non-integer coordinates `(sky_x[i],
sky_y[j])`.  Suppose `sky_x[i] = 3.3` and `sky_y[j] = 4.7`. `bkg[i,
j]` will be the result of linear interpolation between the values
`sky_small[3:4, 4:5]`, multiplied by `calibration[i]`.

For coordinates that are out-of-bounds (e.g., `sky_x[i] < 1.0` or
`sky_x[i] > size(sky_small,1)` where the interpolation would index
sky_small values outside the array, the nearest values in the
sky_small array are used (in other words, constant extrapolation),
consistent with the SDSS data model documentation.

For more information, see
https://github.com/jeff-regier/Celeste.jl/wiki/About-SDSS-and-Stripe-82
"""
struct SDSSBackground <: AbstractArray{Float32, 2}
    sky_small::Matrix{Float32} # background flux per pixel, in DNs
    sky_x::Vector{Float32} # interpolation coordinates
    sky_y::Vector{Float32}
    calibration::Vector{Float32} # nMgy per DN for each row

    function SDSSBackground(sky_small, sky_x, sky_y, calibration)
        @assert length(calibration) == length(sky_x)
        new(sky_small, sky_x, sky_y, calibration)
    end
end

Base.size(bkg::SDSSBackground) = (length(bkg.sky_x), length(bkg.sky_y))

function Base.getindex(sky::SDSSBackground, i::Int, j::Int)
    nx, ny = size(sky.sky_small)

    y0 = floor(Int, sky.sky_y[j])
    y1 = y0 + 1
    yw0 = sky.sky_y[j] - y0
    yw1 = 1.0f0 - yw0

    # modify out-of-bounds indicies to 1 or ny
    y0 = min(max(y0, 1), ny)
    y1 = min(max(y1, 1), ny)

    x0 = floor(Int, sky.sky_x[i])
    x1 = x0 + 1
    xw0 = sky.sky_x[i] - x0
    xw1 = 1.0f0 - xw0

    # modify out-of-bounds indicies to 1 or nx
    x0 = min(max(x0, 1), nx)
    x1 = min(max(x1, 1), nx)

    # bi-linear interpolation
    sky_dns = (xw0 * yw0 * sky.sky_small[x0, y0]
             + xw1 * yw0 * sky.sky_small[x1, y0]
             + xw0 * yw1 * sky.sky_small[x0, y1]
             + xw1 * yw1 * sky.sky_small[x1, y1])

    # return sky intensity in nMgy
    sky_dns * sky.calibration[i]
end

"""
read_sky(hdu)

Construct an image of the sky from the sky HDU of a SDSS \"frame\" file.
This is typically the 3rd HDU of the file.
"""
function read_sky(hdu::FITSIO.TableHDU)
    sky_small = squeeze(read(hdu, "ALLSKY"), 3)::Array{Float32, 2}
    sky_x = vec(read(hdu, "XINTERP"))::Vector{Float32}
    sky_y = vec(read(hdu, "YINTERP"))::Vector{Float32}

    # convert sky interpolation coordinates from 0-indexed to 1-indexed
    for i=1:length(sky_x)
        sky_x[i] += 1.0f0
    end
    for i=1:length(sky_y)
        sky_y[i] += 1.0f0
    end

    # sky intensity has to be strictly greater than 0 for the elbo
    # to be defined
    @assert all((x)-> x > 1e-12, sky_small)

    return sky_small, sky_x, sky_y
end

"""
read_frame(f::FITSIO.FITS)

Read an SDSS \"frame\" FITS file and return a 4-tuple:

- `image`: sky-subtracted and calibrated 2-d image.
- `calibration`: 1-d \"calibration\".
- `sky`: 2-d sky that was subtracted.
- `wcs`: WCSTransform constructed from image header.
"""
function read_frame(f::FITSIO.FITS)
    hdr = FITSIO.read_header(f[1], String)::String
    image = read(f[1])::Array{Float32, 2}  # sky-subtracted & calibrated data, in nMgy
    calibration = read(f[2])::Vector{Float32}
    sky_small, sky_x, sky_y = read_sky(f[3])

    sky = SDSSBackground(sky_small, sky_x, sky_y, calibration)

    wcs = WCS.from_header(hdr)[1]

    return image, calibration, sky, wcs
end

"""
read_field_gains(fname, fieldnum)

Return the image gains for field number `fieldnum` in an SDSS
\"photoField\" file `f`.
"""
function read_field_gains(f::FITSIO.FITS, fieldnum)
    fieldnums = read(f[2], "FIELD")::Vector{Int32}
    gains = read(f[2], "GAIN")::Array{Float32, 2}

    # Find first occurance of `fieldnum` and return the corresponding gain.
    for i=1:length(fieldnums)
        if fieldnums[i] == fieldnum
            bandgains = gains[:, i]
            return Dict(zip("ugriz", bandgains))
        end
    end

    error("field number $fieldnum not found in file: $fname")
end

"""
read_mask(f::FITSIO.FITS[, mask_planes])

Read a \"fpM\"-format SDSS file and return masked image ranges,
based on `mask_planes`. Returns two `Vector{UnitRange{Int}}`,
giving the range of x and y indicies to be masked.
"""
function read_mask(f::FITSIO.FITS, mask_planes=DEFAULT_MASK_PLANES)
    # The last (12th) HDU contains a key describing what each of the
    # other HDUs are. Use this to find the indicies of all the relevant
    # HDUs (those with attributeName matching a value in `mask_planes`).
    value = read(f[12], "Value")::Vector{Int32}
    def = read(f[12], "defName")::Vector{String}
    attribute = read(f[12], "attributeName")::Vector{String}

    # initialize return values
    xranges = UnitRange{Int}[]
    yranges = UnitRange{Int}[]

    # Loop over keys and check if each is a mask plane we're interested in
    # (those with defName == "S_MASKTYPE" and "attributeName" in mask_planes).
    # If so, read from the corresponding HDU and construct ranges to mask.
    for i=1:length(value)
        if (def[i] == "S_MASKTYPE" && attribute[i] in mask_planes)

            # `value` starts from 0, but first table hdu is hdu number 2
            hdunum = value[i] + 2

            cmin = read(f[hdunum], "cmin")::Vector{Int32}
            cmax = read(f[hdunum], "cmax")::Vector{Int32}
            rmin = read(f[hdunum], "rmin")::Vector{Int32}
            rmax = read(f[hdunum], "rmax")::Vector{Int32}

            # "c" ("column") refers to mask's NAXIS1 (x axis);
            # "r" ("row") refers to mask's NAXIS2 (y axis).
            # cmin/cmax and rmin/rmax are 0-based and inclusive, so we add 1
            # to both.
            for j=1:length(cmin)
                push!(xranges, (cmin[j] + 1):(cmax[j] + 1))
                push!(yranges, (rmin[j] + 1):(rmax[j] + 1))
            end
        end
    end

    return xranges, yranges
end

# -----------------------------------------------------------------------------
# PSF-related functions

"""
    SDSSPSFMap <: AbstractPSFMap

SDSS representation of a spatially variable PSF. The PSF is represented as
a weighted combination of eigenimages (stored in `rrows`), where the weights
vary smoothly across the image as a polynomial of the form

```
weight[k](x, y) = sum_{i,j} cmat[i, j, k] * (rcs * x)^i (rcs * y)^j
```

where `rcs` is a coordinate transformation and `x` and `y` are zero-indexed.
"""
struct SDSSPSFMap <: AbstractPSFMap
    rrows::Array{Float64,2}  # A matrix of flattened eigenimages.
    rnrow::Int  # The number of rows in an eigenimage.
    rncol::Int  # The number of columns in an eigenimage.
    cmat::Array{Float64,3}  # The coefficients of the weight polynomial

    function SDSSPSFMap(rrows::Array{Float64, 2},
                        rnrow::Integer, rncol::Integer,
                        cmat::Array{Float64, 3})
        # rrows contains eigen images. Each eigen image is along the first
        # dimension in a flattened form. Check that dimensions match up.
        @assert size(rrows, 1) == rnrow * rncol

        # The second dimension is the number of eigen images, which should
        # match the number of coefficient arrays.
        @assert size(rrows, 2) == size(cmat, 3)

        return new(rrows, Int(rnrow), Int(rncol), cmat)
    end
end

"""
    (psfmap::SDSSPDFMap)(x, y)

Evaluate the PSF at the given image coordinates. The size of the result is
will be `(psf.rnrow, psf.rncol)`, with the PSF (presumably) centered in the
stamp.

This function was originally based on the function sdss_psf_at_points
in astrometry.net:
https://github.com/dstndstn/astrometry.net/blob/master/util/sdss_psf.py
"""
function (psfmap::SDSSPSFMap)(x::Real, y::Real)
    const RCS = 0.001  # A coordinate transform to keep polynomial
                       # coefficients to a reasonable size.
    nk = size(psfmap.rrows, 2)  # number of eigen images.

    # initialize output stamp
    stamp = zeros(psfmap.rnrow, psfmap.rncol)

    # Loop over eigen images
    for k=1:nk
        # calculate the weight for the k-th eigen image from psfmap.cmat.
        # Note that the image coordinates and coefficients are intended
        # to be zero-indexed.
        w = 0.0
        for j=1:size(psfmap.cmat, 2), i=1:size(psfmap.cmat, 1)
            w += (psfmap.cmat[i, j, k] *
                  (RCS * (x - 1.0))^(i-1) * (RCS * (y - 1.0))^(j-1))
        end

        # add the weighted k-th eigen image to the output stamp
        for i=1:length(stamp)
            stamp[i] += w * psfmap.rrows[i, k]
        end
    end

    return stamp
end

"""
read_psfmap(fitsfile, band)

Read spatially variable PSF map for the given band ('u', 'g', 'r', 'i', or 'z')
from the open FITSIO.FITS instance `fitsfile`, which should be a SDSS
\"psField\" file, and return a SDSSPSFMap instance representing the spatially
variable PSF.
"""
function read_psfmap(fitsfile::FITSIO.FITS, band::Char)

    # get the HDU, assuming PSF hdus are in order after primary header
    extnum = 1 + BAND_CHAR_TO_NUM[band]
    hdu = fitsfile[extnum]::FITSIO.TableHDU

    nrows = FITSIO.read_key(hdu, "NAXIS2")[1]::Int
    nrow_b = Int((read(hdu, "nrow_b")::Vector{Int32})[1])
    ncol_b = Int((read(hdu, "ncol_b")::Vector{Int32})[1])
    rnrow = (read(hdu, "rnrow")::Vector{Int32})[1]
    rncol = (read(hdu, "rncol")::Vector{Int32})[1]
    cmat_raw = read(hdu, "c")::Array{Float32, 3}
    rrows_raw = read(hdu, "rrows")::Array{Array{Float32,1},1}

    # Only the first (nrow_b, ncol_b) submatrix of cmat is used for
    # reasons obscure to the author.
    cmat = Array{Float64,3}(nrow_b, ncol_b, size(cmat_raw, 3))
    for k=1:size(cmat_raw, 3), j=1:nrow_b, i=1:ncol_b
        cmat[i, j, k] = cmat_raw[i, j, k]
    end

    # convert rrows to Array{Float64, 2}, assuming each row is the same length.
    rrows = Matrix{Float64}(length(rrows_raw[1]), length(rrows_raw))
    for i=1:length(rrows_raw)
        rrows[:, i] = rrows_raw[i]
    end

    return SDSSPSFMap(rrows, rnrow, rncol, cmat)
end

"""
Read a source catalog from an SDSS \"photoObj\" FITS file for a given
run/camcol/field combination, returning a Vector of columns and a Vector
of column names.

This is currently pretty specific to Celeste because it returns only columns
that we're interested in, and it does some transformations of some of the
columns.

The photoObj file format is documented here:
https://data.sdss.org/datamodel/files/BOSS_PHOTOOBJ/RERUN/RUN/CAMCOL/photoObj.html
"""
function read_photoobj(f::FITSIO.FITS, band::Char='r')
    b = BAND_CHAR_TO_NUM[band]

    # sometimes the expected table extension is only an empty generic FITS
    # header, indicating no objects. We check explicitly for this case here
    # and return an empty table
    if (length(f) < 2) || (!isa(f[2], FITSIO.TableHDU))
        catalog = Dict("objid"=>String[],
                       "ra"=>Float64[],
                       "dec"=>Float64[],
                       "thing_id"=>Int32[],
                       "mode"=>UInt8[],
                       "is_star"=>BitVector(),
                       "is_gal"=>BitVector(),
                       "frac_dev"=>Float32[],
                       "ab_exp"=>Float32[],
                       "theta_exp"=>Float32[],
                       "phi_exp"=>Float32[],
                       "ab_dev"=>Float32[],
                       "theta_dev"=>Float32[],
                       "phi_dev"=>Float32[],
                       "phi_offset"=>Float32[])
        for (b, n) in BAND_CHAR_TO_NUM
            catalog["psfflux_$b"] = Float32[]
            catalog["compflux_$b"] = Float32[]
            catalog["expflux_$b"] = Float32[]
            catalog["devflux_$b"] = Float32[]
        end
        return catalog
    end

    hdu = f[2]::FITSIO.TableHDU

    objid = read(hdu, "objid")::Vector{String}
    ra = read(hdu, "ra")::Vector{Float64}
    dec = read(hdu, "dec")::Vector{Float64}
    mode = read(hdu, "mode")::Vector{UInt8}
    thing_id = read(hdu, "thing_id")::Vector{Int32}

    # Get "bright" objects.
    # (In objc_flags, the bit pattern Int32(2) corresponds to bright objects.)
    objc_flags = read(hdu, "objc_flags")::Vector{Int32}
    objc_flags2 = read(hdu, "objc_flags2")::Vector{Int32}

    # bright, saturated, or large
    bad_flags1 = objc_flags .& UInt32(1^2 + 2^18 + 2^24) .!= 0

    # nopeak, DEBLEND_DEGENERATE, or saturated center
    bad_flags2 = objc_flags2 .& UInt32(2^14 + 2^18 + 2^11) .!= 0

    has_child = read(hdu, "nchild")::Vector{Int16} .> 0

    # 6 = star, 3 = galaxy, others can be ignored.
    objc_type = read(hdu, "objc_type")::Vector{Int32}
    is_star = objc_type .== 6
    is_gal = objc_type .== 3
    bad_type = ((!).(is_star .| is_gal))

    fracdev = vec((read(hdu, "fracdev")::Matrix{Float32})[b, :])

    # determine mask for rows
    has_dev = fracdev .> 0.
    has_exp = fracdev .< 1.
    is_comp = (has_dev .& has_exp)
    is_bad_fracdev = ((fracdev .< 0.) .| (fracdev .> 1))

    # TODO: We don't really want to exclude objects entirely just for being
    # bright: we just don't want to use for scoring (since
    # they're very saturated, presumably).
    mask = (!).(is_bad_fracdev .| bad_type .| bad_flags1 .| bad_flags2 .| has_child)

    # Read the fluxes.
    # Record the cmodelflux if the galaxy is composite, otherwise use
    # the flux for the appropriate type.
    psfflux = read(hdu, "psfflux")::Matrix{Float32}
    cmodelflux = read(hdu, "cmodelflux")::Matrix{Float32}
    devflux = read(hdu, "devflux")::Matrix{Float32}
    expflux = read(hdu, "expflux")::Matrix{Float32}

    # We actually only store the following properties for the given band
    # (see below in Dict construction) NB: the phi quantites in tractor are
    # multiplied by -1.
    phi_dev_deg = read(hdu, "phi_dev_deg")::Matrix{Float32}
    phi_exp_deg = read(hdu, "phi_exp_deg")::Matrix{Float32}
    phi_offset = read(hdu, "phi_offset")::Matrix{Float32}

    theta_dev = read(hdu, "theta_dev")::Matrix{Float32}
    theta_exp = read(hdu, "theta_exp")::Matrix{Float32}

    ab_exp = read(hdu, "ab_exp")::Matrix{Float32}
    ab_dev = read(hdu, "ab_dev")::Matrix{Float32}

    # construct result catalog
    catalog = Dict("objid"=>objid[mask],
                   "ra"=>ra[mask],
                   "dec"=>dec[mask],
                   "thing_id"=>thing_id[mask],
                   "mode"=>mode[mask],
                   "is_star"=>is_star[mask],
                   "is_gal"=>is_gal[mask],
                   "frac_dev"=>fracdev[mask],
                   "ab_exp"=>vec(ab_exp[b, mask]),
                   "theta_exp"=>vec(theta_exp[b, mask]),
                   "phi_exp"=>vec(phi_exp_deg[b, mask]),
                   "ab_dev"=>vec(ab_dev[b, mask]),
                   "theta_dev"=>vec(theta_dev[b, mask]),
                   "phi_dev"=>vec(phi_dev_deg[b, mask]),
                   "phi_offset"=>vec(phi_offset[b, mask]))
    for (b, n) in BAND_CHAR_TO_NUM
        catalog["psfflux_$b"] = vec(psfflux[n, mask])
        catalog["compflux_$b"] = vec(cmodelflux[n, mask])
        catalog["expflux_$b"] = vec(expflux[n, mask])
        catalog["devflux_$b"] = vec(devflux[n, mask])
    end

    # left over note about fluxes above: Use the cmodelflux if the
    # galaxy is composite, otherwise use the flux for the appropriate type.

    return catalog
end

"""
Convert from a catalog in dictionary-of-arrays, as returned by
read_photoobj to Vector{CatalogEntry}.
"""
function convert(::Type{Vector{CatalogEntry}}, catalog::Dict)
    out = CatalogEntry[]

    for i=1:length(catalog["objid"])
        worldcoords = [catalog["ra"][i], catalog["dec"][i]]
        frac_dev = catalog["frac_dev"][i]

        # Fill star and galaxy flux
        star_fluxes = zeros(5)
        gal_fluxes = zeros(5)
        for (j, band) in enumerate(['u', 'g', 'r', 'i', 'z'])
            # Make negative fluxes positive.
            psfflux = max(catalog["psfflux_$band"][i], 1e-6)
            devflux = max(catalog["devflux_$band"][i], 1e-6)
            expflux = max(catalog["expflux_$band"][i], 1e-6)

            # galaxy flux is a weighted sum of the two components.
            star_fluxes[j] = psfflux
            gal_fluxes[j] = frac_dev * devflux + (1.0 - frac_dev) * expflux
        end

        # For shape parameters, we use the dominant component (dev or exp)
        usedev = frac_dev > 0.5
        fits_ab = usedev? catalog["ab_dev"][i] : catalog["ab_exp"][i]
        fits_phi = usedev? catalog["phi_dev"][i] : catalog["phi_exp"][i]
        fits_theta = usedev? catalog["theta_dev"][i] : catalog["theta_exp"][i]

        # effective radius
        re_arcsec = max(fits_theta, 1. / 30)
        re_pixel = re_arcsec / 0.396

        fits_phi -= catalog["phi_offset"][i]

        # fits_phi is now degrees counter-clockwise from vertical
        # in pixel coordinates.

        # Celeste's phi measures radians counter-clockwise from vertical.
        celeste_phi_rad = fits_phi * (pi / 180)

        entry = CatalogEntry(worldcoords, catalog["is_star"][i], star_fluxes,
                             gal_fluxes, frac_dev, fits_ab, celeste_phi_rad,
                             re_pixel)
        push!(out, entry)
    end

    return out
end

# -----------------------------------------------------------------------------
# SDSSDataSet

"""
    SDSSDataSet(basedir; dirlayout = :celeste, compressed = false,
                slurp = false, iostrategy = :plain)

Construct an SDSSDataSet instance with options for where the data is located
(`stagedir`), how the directory structure is laid out (`dirlayout`), whether
the images are compressed (`compressed`), and how to read the images (`slurp`
and `iostrategy`).

`dirlayout` is either `:celeste`, or `:sdss`. The latter is the directory
structure used on the SDSS file server.

`slurp=true` results in reading FITS files into memory all at once, bypassing
the buffering used in the FITSIO library.

`iostrategy` can be `:plain` or `:masterrpc`.
"""
struct SDSSDataSet <: SurveyDataSet
    basedir::String
    dirlayout::Symbol
    compressed::Bool
    slurp::Bool
    iostrategy::Symbol

    function SDSSDataSet(basedir::String; dirlayout::Symbol = :celeste,
                         compressed::Bool = false, slurp::Bool = false,
                         iostrategy::Symbol = :plain)
        if !(dirlayout == :celeste || dirlayout == :sdss)
            error("unknown dirlayout")
        end
        if !(iostrategy == :plain || iostrategy == :masterrpc)
            error("unknown iostrategy")
        end
        return new(basedir, dirlayout, compressed, slurp, iostrategy)
    end
end

Base.:(==)(x::SDSSDataSet, y::SDSSDataSet) =
    (x.basedir == y.basedir &&
     x.dirlayout == y.dirlayout &&
     x.compressed == y.compressed &&
     x.slurp == y.slurp &&
     x.iostrategy == y.iostrategy)

struct RunCamcolField
    run::Int16
    camcol::UInt8
    field::Int16
end

function RunCamcolField(run::String, camcol::String, field::String)
    RunCamcolField(
        parse(Int16, run),
        parse(UInt8, camcol),
        parse(Int16, field))
end

# Helpers for SDSSDataSet file paths: used internally only.
abstract type SDSSFileDesc; end

struct PhotoObj <: SDSSFileDesc
    rcf::RunCamcolField
end
struct PhotoField <: SDSSFileDesc
    run::Int16
    camcol::UInt8
end
struct PsField <: SDSSFileDesc
    rcf::RunCamcolField
end
struct Frame <: SDSSFileDesc
    rcf::RunCamcolField
    band::Char
end
struct Mask <: SDSSFileDesc
    rcf::RunCamcolField
    band::Char
end
struct FieldExtents <: SDSSFileDesc; end
rcf(img::SDSSFileDesc) = img.rcf
rcf(img::PhotoField) = RunCamcolField(img.run, img.camcol, -1)

filename(p::PhotoObj) = "photoObj-$(dec(p.rcf.run,6))-$(dec(p.rcf.camcol))-$(dec(p.rcf.field,4)).fits"
filename(p::PhotoField) = "photoField-$(dec(p.run,6))-$(p.camcol).fits"
filename(p::PsField) = "psField-$(dec(p.rcf.run,6))-$(p.rcf.camcol)-$(dec(p.rcf.field,4)).fit"
filename(p::Frame) = "frame-$(p.band)-$(dec(p.rcf.run,6))-$(p.rcf.camcol)-$(dec(p.rcf.field,4)).fits"
filename(p::Mask) = "fpM-$(dec(p.rcf.run,6))-$(p.band)$(p.rcf.camcol)-$(dec(p.rcf.field,4)).fit"
filename(p::FieldExtents) = "field_extents.fits"

function compute_fname(dataset::SDSSDataSet, img::SDSSFileDesc)
    isa(img, FieldExtents) && return (joinpath(dataset.basedir, filename(img)), nothing)
    imgrcf = rcf(img)
    basedir = dataset.basedir
    if dataset.dirlayout == :celeste
        subdir = "$basedir/$(imgrcf.run)/$(imgrcf.camcol)"
        !isa(img, PhotoField) && (subdir = joinpath(subdir, string(imgrcf.field)))
    elseif dataset.dirlayout == :sdss
        if typeof(img) <: Union{Mask, PsField}
            # Photometric reductions
            subdir = joinpath(basedir, "boss/photo/redux/301", "$(imgrcf.run)/objcs/$(imgrcf.camcol)")
        elseif typeof(img) <: PhotoField
            subdir = joinpath(basedir, "boss/photoObj/301", "$(imgrcf.run)")
        elseif typeof(img) <: PhotoObj
            subdir = joinpath(basedir, "boss/photoObj/301", "$(imgrcf.run)/$(imgrcf.camcol)")
        elseif typeof(img) <: Frame
            subdir = joinpath(basedir, "boss/photoObj/frames/301", "$(imgrcf.run)/$(imgrcf.camcol)")
        end
    end
    fname = joinpath(subdir, filename(img))
    compression = nothing
    (dataset.compressed && isa(img, Frame)) && (fname *= ".bz2"; compression = Bzip2Decompression())
    (dataset.compressed && isa(img, Mask)) && (fname *= ".gz"; compression = GzipDecompression())
    fname, compression
end

function slurp_fits(fname::String)
    is_linux() || is_bsd() || error("Slurping not implemented on this OS")
    # Do this with explicit POSIX calls, to make sure to avoid
    # any well-intended, but ultimately unhelpful intermediate
    # buffering.
    fd = ccall(:open, Cint, (Ptr{UInt8}, Cint, Cint), fname, Base.Filesystem.JL_O_RDONLY, 0)
    systemerror("open($fname)", fd == -1)
    stat_struct = zeros(UInt8, ccall(:jl_sizeof_stat, Int32, ()))
    # Can't use Base.stat because threading
    ret = ccall(:jl_fstat, Cint, (Cint, Ptr{UInt8}), fd, stat_struct)
    ret == 0 || throw(Base.UVError("stat",r))
    size = Base.Filesystem.StatStruct(stat_struct).size
    data = Array{UInt8}(size)
    rsize = ccall(:read, Cint, (Cint, Ptr{UInt8}, Csize_t), fd, data, size)
    systemerror("read", rsize != size)
    r = ccall(:close, Cint, (Cint,), fd)
    systemerror("close", r == -1)
    data
end

function _read_raw_data(dataset, img)
    fname, compression = compute_fname(dataset, img)
    data = slurp_fits(fname)
    data, compression
end

function open_fits(dataset::SDSSDataSet, img::SDSSFileDesc)
    if dataset.iostrategy == :plain
        fname, compression = compute_fname(dataset, img)
        if dataset.slurp || compression != nothing
            data = slurp_fits(fname)
            (compression != nothing) && (data = transcode(compression, data))
            return FITSIO.FITS(data)
        else
            return FITSIO.FITS(fname)
        end
    elseif dataset.iostrategy == :masterrpc
        data, compression = remotecall_fetch(_read_raw_data, 1, dataset, img)
        (compression != nothing) && (data = transcode(compression, data))
        return FITSIO.FITS(data)
    end
end

"""
Query the SDSS database for all fields that overlap the given RA, Dec range.
"""
function _get_overlapping_field_extents(dataset::SDSSDataSet, box::BoundingBox)
    f = open_fits(dataset, FieldExtents())

    hdu = f[2]::FITSIO.TableHDU

    # read in the entire table.
    all_run = read(hdu, "run")::Vector{Int16}
    all_camcol = read(hdu, "camcol")::Vector{UInt8}
    all_field = read(hdu, "field")::Vector{Int16}
    all_ramin = read(hdu, "ramin")::Vector{Float64}
    all_ramax = read(hdu, "ramax")::Vector{Float64}
    all_decmin = read(hdu, "decmin")::Vector{Float64}
    all_decmax = read(hdu, "decmax")::Vector{Float64}

    close(f)

    ret = Tuple{RunCamcolField, BoundingBox}[]

    # The ramin, ramax, etc is a bit unintuitive because we're looking
    # for any overlap.
    for i in eachindex(all_ramin)
        if (all_ramax[i] > box.ramin && all_ramin[i] < box.ramax &&
                all_decmax[i] > box.decmin && all_decmin[i] < box.decmax)
            cur_box = BoundingBox(all_ramin[i], all_ramax[i],
                                  all_decmin[i], all_decmax[i])
            cur_fe = RunCamcolField(all_run[i], all_camcol[i], all_field[i])
            push!(ret, (cur_fe, cur_box))
        end
    end

    return ret
end

"""
    load_field_images(strategy::IOStrategy, rcf::RunCamcolField)

Read a SDSS run/camcol/field into a vector of 5 Celeste Images (one in each
SDSS bandpass).
"""
function load_field_images(dataset::SDSSDataSet, rcf::RunCamcolField)
    # read gain for each band
    f = open_fits(dataset, PhotoField(rcf.run, rcf.camcol))
    gains = read_field_gains(f, rcf.field)
    close(f)

    # open FITS file containing PSF for each band
    psffile = open_fits(dataset, PsField(rcf))

    images = Vector{Image}(5)

    for (b, band) in enumerate(['u', 'g', 'r', 'i', 'z'])
        # load image data
        f = open_fits(dataset, Frame(rcf, band))
        pixels, calibration, sky, wcs = read_frame(f)
        close(f)

        # read mask and apply it to pixel data
        f = open_fits(dataset, Mask(rcf, band))
        mask_xranges, mask_yranges = read_mask(f)
        close(f)
        for i=1:length(mask_xranges)
            pixels[mask_xranges[i], mask_yranges[i]] = NaN
        end

        psfmap = read_psfmap(psffile, band)

        # fit Celeste PSF at center of image
        nx, ny = size(pixels)
        psfstamp = psfmap(nx / 2., ny / 2.)
        celeste_psf = PSF.fit_raw_psf_for_celeste(psfstamp, 2)[1]

        nelec_per_nmgy = gains[band] ./ calibration

        # scale pixels to raw electron counts
        @assert length(nelec_per_nmgy) == nx
        @inbounds for j=1:ny, i=1:nx
            pixels[i, j] = nelec_per_nmgy[i] * (pixels[i, j] + sky[i, j])
        end

        images[b] = Image(pixels, b, wcs, celeste_psf, sky, nelec_per_nmgy,
                          psfmap)

    end

    close(psffile)

    return images
end

function load_field_images(dataset::SDSSDataSet,
                           rcfs::Vector{RunCamcolField})
    images = Image[]
    for rcf in rcfs
        append!(images, load_field_images(dataset, rcf))
    end
    return images
end

# TODO: We'd like to remove this method in favor of load_field_catalog, but it
# is still being used in the benchmark code. First the benchmark code should
# be updated to use a Celeste Catalog rather than a raw SDSS catalog, if that
# makes sense there.
function read_photoobj(dataset::SDSSDataSet, rcf::RunCamcolField)
    f = open_fits(dataset, PhotoObj(rcf))
    rawcatalog = read_photoobj(f, 'r')
    close(f)
    return rawcatalog
end

function load_field_catalog(dataset::SDSSDataSet, rcf::RunCamcolField)
    rawcatalog = read_photoobj(dataset, rcf)
    return convert(Vector{CatalogEntry}, rawcatalog)
end

# load all images overlapping a box
function load_images(dataset::SDSSDataSet, box::BoundingBox)
    fes = _get_overlapping_field_extents(dataset, box)
    rcfs = [fe[1] for fe in fes]
    return load_field_images(dataset, rcfs)
end

end  # module