text:
silentSolver() = GurobiSolver(OutputFlag=0) # TODO: how to silence in a solver-agnostic way? The trick given in JuMP FAQ does not work...

facts("Solvers") do
  const year = 2000
  const a = 1.
  const a_discharge = a * 1.m^3/s
  const b = 2.
  const b_discharge = b * 1.m^3/s

  const riv_name = "Riv"
  const r_name = "Reservoir: tests for solvers"
  const r_min_capacity = 2.e8m^3
  const r_max_capacity = 8.e8m^3

  const time_horizon = 5
  const year = 2000
  const period = Day(1)

  context("Helpers: purpose (both purposes, dependence on period)") do
    purposes = DeterministicPurposes(drinkingWater=a_discharge, environmentalFlow=b_discharge)
    @fact hlp_JuMP_purposes_discharge(Model(), purposes, Day(1)) --> roughly((a + b) * 86400.e-6)
    @fact hlp_JuMP_purposes_discharge(Model(), purposes, Day(2)) --> roughly((a + b) * 86400.e-6 * 2)
    @fact hlp_JuMP_purposes_discharge(Model(), purposes, Week(1)) --> roughly((a + b) * 86400.e-6 * 7)
  end

  context("Helpers: purpose (both purposes, same role for both)") do
    # Both purposes have the same impact: switch the discharges, should get the same results as before.
    purposes = DeterministicPurposes(drinkingWater=b_discharge, environmentalFlow=a_discharge)
    @fact hlp_JuMP_purposes_discharge(Model(), purposes, Day(1)) --> roughly((a + b) * 86400.e-6)
    @fact hlp_JuMP_purposes_discharge(Model(), purposes, Day(2)) --> roughly((a + b) * 86400.e-6 * 2)
    @fact hlp_JuMP_purposes_discharge(Model(), purposes, Week(1)) --> roughly((a + b) * 86400.e-6 * 7)
  end

  context("Helpers: purpose (just drinking water)") do
    purposes = DeterministicPurposes(drinkingWater=a_discharge)
    @fact hlp_JuMP_purposes_discharge(Model(), purposes, Day(1)) --> roughly(a * 86400.e-6)
    @fact hlp_JuMP_purposes_discharge(Model(), purposes, Day(2)) --> roughly(a * 86400.e-6 * 2)
    @fact hlp_JuMP_purposes_discharge(Model(), purposes, Week(1)) --> roughly(a * 86400.e-6 * 7)
  end

  context("Helpers: purpose (just environmental flow)") do
    purposes = DeterministicPurposes(drinkingWater=b_discharge)
    @fact hlp_JuMP_purposes_discharge(Model(), purposes, Day(1)) --> roughly(b * 86400.e-6)
    @fact hlp_JuMP_purposes_discharge(Model(), purposes, Day(2)) --> roughly(b * 86400.e-6 * 2)
    @fact hlp_JuMP_purposes_discharge(Model(), purposes, Week(1)) --> roughly(b * 86400.e-6 * 7)
  end

  context("Helpers: output (nonlinear modes)") do
    # Recognised (non)linear modes.
    @fact _hlp_outputs_allowed_nonlinear_mode(:None) --> nothing
    @fact _hlp_outputs_allowed_nonlinear_mode(:LinearApprox) --> nothing
    @fact _hlp_outputs_allowed_nonlinear_mode(:PiecewiseLinear) --> nothing
    @fact _hlp_outputs_allowed_nonlinear_mode(:Cone) --> nothing
    @fact_throws _hlp_outputs_allowed_nonlinear_mode(:SomethingElseThatIsNotRecognised)
  end

  function outputHelpersTestLinear(period::Period, length::Int)
    riv = NaturalRiver(name=riv_name, scenarios=[TimeArray(gen_dates(length, year), gen_scenario(length), ["Discharge (m³/s)"], period)])
    outputs = ConstantDamOutputs(bottomOutlets=a_discharge, spillway=b_discharge)
    r = Reservoir(name=r_name, capacity=(r_min_capacity, r_max_capacity), rivers_in=[riv], outputs=outputs)

    # Smallest optimisation model so that this works.
    m = Model(solver=silentSolver())
    @variable(m, output[1:length] >= 0)
    @variable(m, ReservoirManagement._from_unitful(r_min_capacity) <= volume[1:length] <= ReservoirManagement._from_unitful(r_max_capacity))
    @objective(m, Max, sum(output))
    for t in 1:length
      hlp_JuMP_outputs_constraint(m, output[t], volume[t], r, nonlinear_mode=:None)
    end
    status = solve(m)

    # Check the output values: must be the maximum possible to release (the volume is free, there is zero constraint on it).
    @fact status --> :Optimal
    @fact getvalue(sum(output)) --> roughly((a + b) * 86400.e-6 * length * Dates.days(period))
  end

  context("Helpers: output (linear)") do
    # Test the constraints are the expected ones: maximise the output, that must be the input value. It depends on the
    # period and complete duration of optimisation.
    outputHelpersTestLinear(Day(1), 1)
    outputHelpersTestLinear(Day(2), 7)
    outputHelpersTestLinear(Week(1), 12)
  end

  function outputHelpersTestLinearConditional(period::Period, length::Int)
    # bathymetry: h = 280 [m] + V [10^6 m^3] * .1 [m / 10^6 m^3]
    # Bounds: V in [r_min_capacity, r_max_capacity], [m^3] --> 200 to 800 [10^6 m^3].
    # Hence for reservoir level: [300, 360] m. Condition in the middle
    geom = LinearReservoirBathymetry(280., .1)
    riv = NaturalRiver(name=riv_name, scenarios=[TimeArray(gen_dates(length, year), gen_scenario(length), ["Discharge (m³/s)"], period)])

    spillway_height = 330.
    r_cap_spillway = computeVolume(geom, spillway_height)
    outputs = ConstantDamOutputs(bottomOutlets=a_discharge, spillway=b_discharge, spillway_minLevel=spillway_height * 1.m)

    r = Reservoir(name=r_name, capacity=(r_min_capacity, r_max_capacity), rivers_in=[riv], outputs=outputs, bathymetry=geom)

    ## Disable the spillway.
    # Smallest optimisation model so that this works: completely disable the spillway by working on the maximum volume.
    m_1 = Model(solver=silentSolver())
    @variable(m_1, output_1[1:length] >= 0)
    @variable(m_1, ReservoirManagement._from_unitful(r_min_capacity) <= volume_1[1:length] <= r_cap_spillway - .1) # Need a small epsilon to be sure to be below the spillway!
    @objective(m_1, Max, sum(output_1))
    for t in 1:length
      hlp_JuMP_outputs_constraint(m_1, output_1[t], volume_1[t], r, nonlinear_mode=:None)
    end
    status_1 = solve(m_1)

    # Check the output values: must be the maximum possible to release (the volume is free,
    # there is zero constraint on it, except to disable the spillway).
    @fact status_1 --> :Optimal
    @fact getvalue(sum(output_1)) --> roughly(a * 86400.e-6 * length * Dates.days(period))

    ## Force enable the spillway.
    m_2 = Model(solver=silentSolver())
    @variable(m_2, output_2[1:length] >= 0)
    @variable(m_2, r_cap_spillway <= volume_2[1:length] <= ReservoirManagement._from_unitful(r_max_capacity)) # With the objective function, the spillway is always active.
    @objective(m_2, Max, sum(output_2))
    for t in 1:length
      hlp_JuMP_outputs_constraint(m_2, output_2[t], volume_2[t], r, nonlinear_mode=:None)
    end
    status_2 = solve(m_2)

    # Check the output values: must be the maximum possible to release (the volume is free,
    # there is zero constraint on it). As a consequence, the volume will be chosen to enable the spillway.
    @fact status_2 --> :Optimal
    @fact getvalue(sum(output_2)) --> roughly((a + b) * 86400.e-6 * length * Dates.days(period))

    ## Spillway free (but release maximised).
    m_3 = Model(solver=silentSolver())
    @variable(m_3, output_3[1:length] >= 0)
    @variable(m_3, ReservoirManagement._from_unitful(r_min_capacity) <= volume_3[1:length] <= ReservoirManagement._from_unitful(r_max_capacity))
    @objective(m_3, Max, sum(output_3))
    for t in 1:length
      hlp_JuMP_outputs_constraint(m_3, output_3[t], volume_3[t], r, nonlinear_mode=:None)
    end
    status_3 = solve(m_3)

    # Check the output values: must be the maximum possible to release (the volume is free,
    # there is zero constraint on it, except to force the spillway to be operational).
    @fact status_3 --> :Optimal
    @fact getvalue(sum(output_3)) --> roughly((a + b) * 86400.e-6 * length * Dates.days(period))
  end

  context("Helpers: output (linear, conditional [reservoir level])") do
    # Test the constraints are the expected ones: maximise the output, that must be the input value. It depends on the
    # period and complete duration of optimisation.
    outputHelpersTestLinearConditional(Day(1), 1)
    outputHelpersTestLinearConditional(Day(2), 7)
    outputHelpersTestLinearConditional(Week(1), 12)

    # If the maximum reservoir level is higher than the condition, show a warning.
    # TODO: no test framework allows this for now, it seems!
  end

  context("Helpers: inflow (one natural river)") do
    # Simplest case: one natural river.
    riv1 = NaturalRiver(name=riv_name, scenarios=[TimeArray(gen_dates(time_horizon, year), gen_scenario(time_horizon), ["Discharge (m³/s)"], period)])
    r = Reservoir(name=r_name, capacity=(r_min_capacity, r_max_capacity), rivers_in=[riv1])

    for t in 1:time_horizon
      @fact hlp_JuMP_inflow_discharge(Model(), r, 1, t) --> roughly(getScenario(riv1, 1, t) * Dates.days(period))
    end
  end

  context("Helpers: inflow (many natural rivers)") do
    riv1 = NaturalRiver(name=riv_name, scenarios=[TimeArray(gen_dates(time_horizon, year), gen_scenario(time_horizon), ["Discharge (m³/s)"], period)])
    riv2 = NaturalRiver(name=riv_name, scenarios=[TimeArray(gen_dates(time_horizon, year), gen_scenario(time_horizon), ["Discharge (m³/s)"], period)])
    riv3 = NaturalRiver(name=riv_name, scenarios=[TimeArray(gen_dates(time_horizon, year), gen_scenario(time_horizon), ["Discharge (m³/s)"], period)])
    r = Reservoir(name=r_name, capacity=(r_min_capacity, r_max_capacity), rivers_in=[riv1, riv2, riv3])

    for t in 1:time_horizon
      @fact hlp_JuMP_inflow_discharge(Model(), r, 1, t) --> roughly((getScenario(riv1, 1, t) + getScenario(riv2, 1, t) + getScenario(riv3, 1, t)) * Dates.days(period))
    end
  end

  context("Helpers: inflow (diverted rivers, imposed)") do
    riv1 = NaturalRiver(name=riv_name, scenarios=[TimeArray(gen_dates(time_horizon, year), gen_scenario(time_horizon), ["Discharge (m³/s)"], period)])
    riv2 = NaturalRiver(name=riv_name, scenarios=[TimeArray(gen_dates(time_horizon, year), gen_scenario(time_horizon), ["Discharge (m³/s)"], period)])
    riv3 = NaturalRiver(name=riv_name, scenarios=[TimeArray(gen_dates(time_horizon, year), gen_scenario(time_horizon), ["Discharge (m³/s)"], period)])
    driv1 = DivertedRiver(name=riv_name, scenarios=[TimeArray(gen_dates(time_horizon, year), gen_scenario(time_horizon), ["Discharge (m³/s)"], period)], environmental_flow=a_discharge, maximum_flow=b_discharge)
    driv2 = DivertedRiver(name=riv_name, scenarios=[TimeArray(gen_dates(time_horizon, year), gen_scenario(time_horizon), ["Discharge (m³/s)"], period)], environmental_flow=a_discharge, maximum_flow=b_discharge)
    r = Reservoir(name=r_name, capacity=(r_min_capacity, r_max_capacity), rivers_in=[riv1, riv2, riv3], rivers_diverted=[driv1, driv2])

    for t in 1:time_horizon
      @fact hlp_JuMP_inflow_discharge(Model(), r, 1, t, imposeDiverted=true) --> roughly((maxAllowableFlow(driv1, 1, t) + maxAllowableFlow(driv2, 1, t) + getScenario(riv1, 1, t) + getScenario(riv2, 1, t) + getScenario(riv3, 1, t)) * Dates.days(period))
    end
  end

  context("Helpers: inflow (diverted rivers, degree of freedom)") do
    # Must optimise, as the inflow becomes an optimisation variable!
    riv1 = NaturalRiver(name=riv_name, scenarios=[TimeArray(gen_dates(time_horizon, year), gen_scenario(time_horizon), ["Discharge (m³/s)"], period)])
    riv2 = NaturalRiver(name=riv_name, scenarios=[TimeArray(gen_dates(time_horizon, year), gen_scenario(time_horizon), ["Discharge (m³/s)"], period)])
    riv3 = NaturalRiver(name=riv_name, scenarios=[TimeArray(gen_dates(time_horizon, year), gen_scenario(time_horizon), ["Discharge (m³/s)"], period)])
    driv1 = DivertedRiver(name=riv_name, scenarios=[TimeArray(gen_dates(time_horizon, year), gen_scenario(time_horizon), ["Discharge (m³/s)"], period)], environmental_flow=a_discharge, maximum_flow=b_discharge)
    driv2 = DivertedRiver(name=riv_name, scenarios=[TimeArray(gen_dates(time_horizon, year), gen_scenario(time_horizon), ["Discharge (m³/s)"], period)], environmental_flow=a_discharge, maximum_flow=b_discharge)
    r = Reservoir(name=r_name, capacity=(r_min_capacity, r_max_capacity), rivers_in=[riv1, riv2, riv3], rivers_diverted=[driv1, driv2])

    m_m = Model(solver=silentSolver())
    @variable(m_m, v_divert[1:1, 1:time_horizon, getDivertedRiverNames(r)] >= 0)
    @variable(m_m, obj >= 0)
    @objective(m_m, Max, obj)
    @constraint(m_m, obj == sum{hlp_JuMP_inflow_discharge(m_m, r, driv1, 1, t, imposeDiverted=v_divert) + hlp_JuMP_inflow_discharge(m_m, r, driv2, 1, t, imposeDiverted=v_divert), t in 1:time_horizon})
    solve(m_m)

    @fact getvalue(obj) --> roughly(sum([(maxAllowableFlow(driv1, 1, t) + maxAllowableFlow(driv2, 1, t)) * Dates.days(period) for t in 1:time_horizon]))
  end

  context("Evaluation: purpose shortage (feasible solution)") do
    # 2 m^3/s available, withdraw 1 m^3/s, no output limit.
    ts = TimeArray(gen_dates(time_horizon, year), [2.e-6 * 86400 * Dates.days(period) for t in 1:time_horizon], ["Discharge (m³/s)"], period)
    riv = NaturalRiver(name=riv_name, scenarios=[ts])
    purpose = DeterministicPurposes(drinkingWater=1.0m^3/s)
    r = Reservoir(name=r_name, capacity=(r_min_capacity, r_max_capacity), rivers_in=[riv], purposes=purpose)

    demands = [sum([getNeed(p, period) for p in getPurposes(r)]) for t in 1:time_horizon]
    sol = minimumRuleCurve(r, verbosity=-1)
    eval = purposeShortage(sol, values(ts), verbosity=-1)

    @fact typeof(eval) --> ReservoirSolution
    @fact typeof(eval.solution) --> ShortageSolution
    @fact typeof(eval.solverOptions) --> ShortageOptions
    @fact typeof(eval.solverStatistics) --> ShortageSolverStatistics
    @fact eval.solver --> :ShortageSolver

    @fact eval.solution.infeasibilityCause --> :None
    @fact eval.solution.delta --> roughly(zeros(time_horizon)) # Structural property and actual expected value.
    @fact sum(eval.solution.delta .<= demands) --> length(eval.solution.delta) # Cannot lack more water than is needed, for each time step.

    @fact computePurposeShortageIndex(eval, purpose) --> roughly(0.) # Structural property and actual expected value.
    @fact computePurposeVulnerabilityIndex(eval, purpose) --> roughly(0.) # Structural property and actual expected value.
  end

  context("Evaluation: purpose shortage (shortage)") do
    # 0 m^3/s available, withdraw 2 m^3/s, no output limit, initial condition forced to the minimum reservoir level.
    ts = TimeArray(gen_dates(time_horizon, year), [0. for t in 1:time_horizon], ["Discharge (m³/s)"], period) # 1.e-6 * 86400 * Dates.days(period)
    riv = NaturalRiver(name=riv_name, scenarios=[ts])
    purpose = DeterministicPurposes(drinkingWater=2.m^3/s)
    r = Reservoir(name=r_name, capacity=(r_min_capacity, r_max_capacity), rivers_in=[riv], purposes=purpose)

    # Compute a rule curve and evaluate it.
    sol = minimumRuleCurve(r, verbosity=-1)
    eval = purposeShortage(sol, values(ts), forceInitialCondition=true, verbosity=-1)

    # Compute data required for the tests.
    minLevel = getSolution(sol)
    rule_curve_delta = [minLevel[t] - (t == 1 ? getMinimumCapacity(r) : minLevel[t - 1]) for t in 1:time_horizon]
    demands = [sum([getNeed(p, period) for p in getPurposes(r)]) + (rule_curve_delta[t] >= 0. ? rule_curve_delta[t] : 0.) for t in 1:time_horizon]

    # Actual tests.
    @fact typeof(eval) --> ReservoirSolution
    @fact typeof(eval.solution) --> ShortageSolution
    @fact typeof(eval.solverOptions) --> ShortageOptions
    @fact typeof(eval.solverStatistics) --> ShortageSolverStatistics
    @fact eval.solver --> :ShortageSolver

    @fact eval.solution.infeasibilityCause --> :Shortage
    @fact sum(eval.solution.delta) --> greater_than(0.1) # Structural property.
    @fact sum(eval.solution.delta) --> roughly(.864) # Actual expected value.
    @fact sum(eval.solution.delta .<= demands) --> length(eval.solution.delta) # Cannot lack more water than is needed, for each time step.

    @fact computePurposeShortageIndex(eval, purpose) --> greater_than(0.1) # Structural property.
    @fact computePurposeShortageIndex(eval, purpose) --> roughly(1.) # Actual expected value.
    @fact computePurposeVulnerabilityIndex(eval, purpose) --> greater_than(0.1) # Structural property.
    @fact computePurposeVulnerabilityIndex(eval, purpose) --> roughly(5.) # Actual expected value.
    # Why does those two indices give very different results? There is no water in, so the lack of water is very large.
    # The first one considers averages over years (for the only year in the problem, 100% of the water is missing).
    # The second one considers every time step (here, five). As the sum happens over the *squares* of ratios, this gives
    # an impressive difference.
  end

  context("Evaluation: purpose shortage (abundance)") do
    # Time series for deriving the rule curve (ts_lo), ensures there is a feasible solution.
    # Check for infeasibility against more extreme scenarios (ts_hi).
    ts_lo = TimeArray(gen_dates(time_horizon, year), [50.e-6 * 86400 * Dates.days(period) for t in 1:time_horizon], ["Discharge (m³/s)"], period) # 10^6 m^3 / period
    ts_hi = TimeArray(gen_dates(time_horizon, year), [50.e-2 * 86400 * Dates.days(period) for t in 1:time_horizon], ["Discharge (m³/s)"], period) # 10^6 m^3 / period
    riv = NaturalRiver(name=riv_name, scenarios=[ts_lo])
    purpose = DeterministicPurposes(drinkingWater=0.01m^3/s)
    outs = ConstantDamOutputs(bottomOutlets=1.m^3/s)
    r = Reservoir(name=r_name, capacity=(r_min_capacity, r_max_capacity), rivers_in=[riv], purposes=purpose, outputs=outs)

    demands = [sum([getNeed(p, period) for p in getPurposes(r)]) for t in 1:time_horizon]
    sol = minimumRuleCurve(r, useOutputs=getOutputIds(r), verbosity=-1)
    eval = purposeShortage(sol, values(ts_hi), useOutputs=getOutputIds(r), forceInitialCondition=true, verbosity=-1)

    @fact typeof(eval) --> ReservoirSolution
    @fact typeof(eval.solution) --> ShortageSolution
    @fact typeof(eval.solverOptions) --> ShortageOptions
    @fact typeof(eval.solverStatistics) --> ShortageSolverStatistics
    @fact eval.solver --> :ShortageSolver

    @fact eval.solution.infeasibilityCause --> :Surplus
    @fact sum(eval.solution.delta) --> less_than(0.0) # Structural property: surplus, hence negative.
    @fact abs(sum(eval.solution.delta)) --> greater_than(0.1) # Structural property: far from 0.
    @fact abs(sum(eval.solution.delta)) --> roughly(215999.56367) # Actual expected value.
    # No test between eval.solution.delta and demands, as purposes are fulfilled (surplus of water, not lack).

    @fact computePurposeShortageIndex(eval, purpose) --> roughly(0.) # Structural property and actual expected value.
    @fact computePurposeVulnerabilityIndex(eval, purpose) --> roughly(0.) # Structural property and actual expected value.
  end
end


text:
# By default, Julia/LLVM does not use fused multiply-add operations (FMAs).
# Since these FMAs can increase the performance of many numerical algorithms,
# we need to opt-in explicitly.
# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details.
@muladd begin

# Save current mesh with some context information as an HDF5 file.
function save_mesh_file(mesh::TreeMesh, output_directory, timestep=0)
  save_mesh_file(mesh, output_directory, timestep, mpi_parallel(mesh))
end

function save_mesh_file(mesh::TreeMesh, output_directory, timestep,
                        mpi_parallel::Val{false})
  # Create output directory (if it does not exist)
  mkpath(output_directory)

  # Determine file name based on existence of meaningful time step
  if timestep > 0
    filename = joinpath(output_directory, @sprintf("mesh_%06d.h5", timestep))
  else
    filename = joinpath(output_directory, "mesh.h5")
  end

  # Open file (clobber existing content)
  h5open(filename, "w") do file
    # Add context information as attributes
    n_cells = length(mesh.tree)
    attributes(file)["mesh_type"] = get_name(mesh)
    attributes(file)["ndims"] = ndims(mesh)
    attributes(file)["n_cells"] = n_cells
    attributes(file)["n_leaf_cells"] = count_leaf_cells(mesh.tree)
    attributes(file)["minimum_level"] = minimum_level(mesh.tree)
    attributes(file)["maximum_level"] = maximum_level(mesh.tree)
    attributes(file)["center_level_0"] = mesh.tree.center_level_0
    attributes(file)["length_level_0"] = mesh.tree.length_level_0
    attributes(file)["periodicity"] = collect(mesh.tree.periodicity)

    # Add tree data
    file["parent_ids"] = @view mesh.tree.parent_ids[1:n_cells]
    file["child_ids"] = @view mesh.tree.child_ids[:, 1:n_cells]
    file["neighbor_ids"] = @view mesh.tree.neighbor_ids[:, 1:n_cells]
    file["levels"] = @view mesh.tree.levels[1:n_cells]
    file["coordinates"] = @view mesh.tree.coordinates[:, 1:n_cells]
  end

  return filename
end

# Save current mesh with some context information as an HDF5 file.
function save_mesh_file(mesh::TreeMesh, output_directory, timestep,
                        mpi_parallel::Val{true})
  # Create output directory (if it does not exist)
  mpi_isroot() && mkpath(output_directory)

  # Determine file name based on existence of meaningful time step
  if timestep >= 0
    filename = joinpath(output_directory, @sprintf("mesh_%06d.h5", timestep))
  else
    filename = joinpath(output_directory, "mesh.h5")
  end

  # Since the mesh is replicated on all ranks, only save from MPI root
  if !mpi_isroot()
    return filename
  end

  # Open file (clobber existing content)
  h5open(filename, "w") do file
    # Add context information as attributes
    n_cells = length(mesh.tree)
    attributes(file)["mesh_type"] = get_name(mesh)
    attributes(file)["ndims"] = ndims(mesh)
    attributes(file)["n_cells"] = n_cells
    attributes(file)["n_leaf_cells"] = count_leaf_cells(mesh.tree)
    attributes(file)["minimum_level"] = minimum_level(mesh.tree)
    attributes(file)["maximum_level"] = maximum_level(mesh.tree)
    attributes(file)["center_level_0"] = mesh.tree.center_level_0
    attributes(file)["length_level_0"] = mesh.tree.length_level_0
    attributes(file)["periodicity"] = collect(mesh.tree.periodicity)

    # Add tree data
    file["parent_ids"] = @view mesh.tree.parent_ids[1:n_cells]
    file["child_ids"] = @view mesh.tree.child_ids[:, 1:n_cells]
    file["neighbor_ids"] = @view mesh.tree.neighbor_ids[:, 1:n_cells]
    file["levels"] = @view mesh.tree.levels[1:n_cells]
    file["coordinates"] = @view mesh.tree.coordinates[:, 1:n_cells]
  end

  return filename
end

# Does not save the mesh itself to an HDF5 file. Instead saves important attributes
# of the mesh, like its size and the type of boundary mapping function.
# Then, within Trixi2Vtk, the StructuredMesh and its node coordinates are reconstructured from
# these attributes for plotting purposes
function save_mesh_file(mesh::StructuredMesh, output_directory)
  # Create output directory (if it does not exist)
  mkpath(output_directory)

  filename = joinpath(output_directory, "mesh.h5")

  # Open file (clobber existing content)
  h5open(filename, "w") do file
    # Add context information as attributes
    attributes(file)["mesh_type"] = get_name(mesh)
    attributes(file)["ndims"] = ndims(mesh)
    attributes(file)["size"] = collect(size(mesh))
    attributes(file)["mapping"] = mesh.mapping_as_string
  end

  return filename
end

# Does not save the mesh itself to an HDF5 file. Instead saves important attributes
# of the mesh, like its size and the corresponding `.mesh` file used to construct the mesh.
# Then, within Trixi2Vtk, the UnstructuredMesh2D and its node coordinates are reconstructured
# from these attributes for plotting purposes
function save_mesh_file(mesh::UnstructuredMesh2D, output_directory)
  # Create output directory (if it does not exist)
  mkpath(output_directory)

  filename = joinpath(output_directory, "mesh.h5")

  # Open file (clobber existing content)
  h5open(filename, "w") do file
    # Add context information as attributes
    attributes(file)["mesh_type"] = get_name(mesh)
    attributes(file)["ndims"] = ndims(mesh)
    attributes(file)["size"] = length(mesh)
    attributes(file)["mesh_filename"] = mesh.filename
    attributes(file)["periodicity"] = collect(mesh.periodicity)
  end

  return filename
end

# Does not save the mesh itself to an HDF5 file. Instead saves important attributes
# of the mesh, like its size and the type of boundary mapping function.
# Then, within Trixi2Vtk, the P4estMesh and its node coordinates are reconstructured from
# these attributes for plotting purposes
function save_mesh_file(mesh::P4estMesh, output_directory, timestep=0)
  # Create output directory (if it does not exist)
  mkpath(output_directory)

  # Determine file name based on existence of meaningful time step
  if timestep > 0
    filename = joinpath(output_directory, @sprintf("mesh_%06d.h5", timestep))
    p4est_filename = @sprintf("p4est_data_%06d", timestep)
  else
    filename = joinpath(output_directory, "mesh.h5")
    p4est_filename = "p4est_data"
  end

  p4est_file = joinpath(output_directory, p4est_filename)

  # Save the complete connectivity and `p4est` data to disk.
  save_p4est!(p4est_file, mesh.p4est)

  # Open file (clobber existing content)
  h5open(filename, "w") do file
    # Add context information as attributes
    attributes(file)["mesh_type"] = get_name(mesh)
    attributes(file)["ndims"] = ndims(mesh)
    attributes(file)["p4est_file"] = p4est_filename

    file["tree_node_coordinates"] = mesh.tree_node_coordinates
    file["nodes"] = Vector(mesh.nodes) # the mesh uses `SVector`s for the nodes
                                       # to increase the runtime performance
                                       # but HDF5 can only handle plain arrays
    file["boundary_names"] = mesh.boundary_names .|> String
  end

  return filename
end

"""
    load_mesh(restart_file::AbstractString; n_cells_max)

Load the mesh from the `restart_file`.
"""
function load_mesh(restart_file::AbstractString; n_cells_max=0, RealT=Float64)
  # Determine mesh filename
  mesh_file = get_restart_mesh_filename(restart_file, Val(mpi_isparallel()))

  if mpi_isparallel()
    return load_mesh_parallel(mesh_file; n_cells_max=n_cells_max, RealT=RealT)
  end

  load_mesh_serial(mesh_file; n_cells_max=n_cells_max, RealT=RealT)
end

function load_mesh_serial(mesh_file::AbstractString; n_cells_max, RealT)
  ndims, mesh_type = h5open(mesh_file, "r") do file
    return read(attributes(file)["ndims"]),
           read(attributes(file)["mesh_type"])
  end

  if mesh_type == "TreeMesh"
    n_cells = h5open(mesh_file, "r") do file
      return read(attributes(file)["n_cells"])
    end
    mesh = TreeMesh(SerialTree{ndims}, max(n_cells, n_cells_max))
    load_mesh!(mesh, mesh_file)
  elseif mesh_type == "StructuredMesh"
    size_, mapping_as_string = h5open(mesh_file, "r") do file
      return read(attributes(file)["size"]),
             read(attributes(file)["mapping"])
    end

    size = Tuple(size_)

    # TODO: `@eval` is evil
    # A temporary workaround to evaluate the code that defines the domain mapping in a local scope.
    # This prevents errors when multiple restart elixirs are executed in one session, where one
    # defines `mapping` as a variable, while the other defines it as a function.
    #
    # This should be replaced with something more robust and secure,
    # see https://github.com/trixi-framework/Trixi.jl/issues/541).
    expr = Meta.parse(mapping_as_string)
    if expr.head == :toplevel
      expr.head = :block
    end

    if ndims == 1
      mapping = @eval function(xi)
        $expr
        mapping(xi)
      end
    elseif ndims == 2
      mapping = @eval function(xi, eta)
        $expr
        mapping(xi, eta)
      end
    else # ndims == 3
      mapping = @eval function(xi, eta, zeta)
        $expr
        mapping(xi, eta, zeta)
      end
    end

    mesh = StructuredMesh(size, mapping; RealT=RealT, unsaved_changes=false,
                      mapping_as_string=mapping_as_string)
  elseif mesh_type == "UnstructuredMesh2D"
    mesh_filename, periodicity_ = h5open(mesh_file, "r") do file
      return read(attributes(file)["mesh_filename"]),
             read(attributes(file)["periodicity"])
    end
    mesh = UnstructuredMesh2D(mesh_filename; RealT=RealT, periodicity=periodicity_,
                                unsaved_changes=false)
  elseif mesh_type == "P4estMesh"
    p4est_filename, tree_node_coordinates,
        nodes, boundary_names_ = h5open(mesh_file, "r") do file
      return read(attributes(file)["p4est_file"]),
             read(file["tree_node_coordinates"]),
             read(file["nodes"]),
             read(file["boundary_names"])
    end

    boundary_names = boundary_names_ .|> Symbol

    p4est_file = joinpath(dirname(mesh_file), p4est_filename)
    # Prevent Julia crashes when `p4est` can't find the file
    @assert isfile(p4est_file)

    p4est = load_p4est(p4est_file, Val(ndims))

    mesh = P4estMesh{ndims}(p4est, tree_node_coordinates,
                            nodes, boundary_names, "", false)
  else
    error("Unknown mesh type!")
  end

  return mesh
end

function load_mesh!(mesh::SerialTreeMesh, mesh_file::AbstractString)
  mesh.current_filename = mesh_file
  mesh.unsaved_changes = false

  # Read mesh file
  h5open(mesh_file, "r") do file
    # Set domain information
    mesh.tree.center_level_0 = read(attributes(file)["center_level_0"])
    mesh.tree.length_level_0 = read(attributes(file)["length_level_0"])
    mesh.tree.periodicity    = Tuple(read(attributes(file)["periodicity"]))

    # Set length
    n_cells = read(attributes(file)["n_cells"])
    resize!(mesh.tree, n_cells)

    # Read in data
    mesh.tree.parent_ids[1:n_cells] = read(file["parent_ids"])
    mesh.tree.child_ids[:, 1:n_cells] = read(file["child_ids"])
    mesh.tree.neighbor_ids[:, 1:n_cells] = read(file["neighbor_ids"])
    mesh.tree.levels[1:n_cells] = read(file["levels"])
    mesh.tree.coordinates[:, 1:n_cells] = read(file["coordinates"])
  end

  return mesh
end

function load_mesh_parallel(mesh_file::AbstractString; n_cells_max, RealT)
  if mpi_isroot()
    ndims_, n_cells = h5open(mesh_file, "r") do file
      read(attributes(file)["ndims"]),
      read(attributes(file)["n_cells"])
    end
    MPI.Bcast!(Ref(ndims_), mpi_root(), mpi_comm())
    MPI.Bcast!(Ref(n_cells), mpi_root(), mpi_comm())
  else
    ndims_ = MPI.Bcast!(Ref(0), mpi_root(), mpi_comm())[]
    n_cells = MPI.Bcast!(Ref(0), mpi_root(), mpi_comm())[]
  end

  mesh = TreeMesh(ParallelTree{ndims_}, max(n_cells, n_cells_max))
  load_mesh!(mesh, mesh_file)

  return mesh
end

function load_mesh!(mesh::ParallelTreeMesh, mesh_file::AbstractString)
  mesh.current_filename = mesh_file
  mesh.unsaved_changes = false

  if mpi_isroot()
    h5open(mesh_file, "r") do file
      # Set domain information
      mesh.tree.center_level_0 = read(attributes(file)["center_level_0"])
      mesh.tree.length_level_0 = read(attributes(file)["length_level_0"])
      mesh.tree.periodicity    = Tuple(read(attributes(file)["periodicity"]))
      MPI.Bcast!(collect(mesh.tree.center_level_0), mpi_root(), mpi_comm())
      MPI.Bcast!(collect(mesh.tree.length_level_0), mpi_root(), mpi_comm())
      MPI.Bcast!(collect(mesh.tree.periodicity),    mpi_root(), mpi_comm())

      # Set length
      n_cells = read(attributes(file)["n_cells"])
      MPI.Bcast!(Ref(n_cells), mpi_root(), mpi_comm())
      resize!(mesh.tree, n_cells)

      # Read in data
      mesh.tree.parent_ids[1:n_cells] = read(file["parent_ids"])
      mesh.tree.child_ids[:, 1:n_cells] = read(file["child_ids"])
      mesh.tree.neighbor_ids[:, 1:n_cells] = read(file["neighbor_ids"])
      mesh.tree.levels[1:n_cells] = read(file["levels"])
      mesh.tree.coordinates[:, 1:n_cells] = read(file["coordinates"])
      @views MPI.Bcast!(mesh.tree.parent_ids[1:n_cells],      mpi_root(), mpi_comm())
      @views MPI.Bcast!(mesh.tree.child_ids[:, 1:n_cells],    mpi_root(), mpi_comm())
      @views MPI.Bcast!(mesh.tree.neighbor_ids[:, 1:n_cells], mpi_root(), mpi_comm())
      @views MPI.Bcast!(mesh.tree.levels[1:n_cells],          mpi_root(), mpi_comm())
      @views MPI.Bcast!(mesh.tree.coordinates[:, 1:n_cells],  mpi_root(), mpi_comm())
    end
  else # non-root ranks
    # Set domain information
    mesh.tree.center_level_0 = MPI.Bcast!(collect(mesh.tree.center_level_0), mpi_root(), mpi_comm())
    mesh.tree.length_level_0 = MPI.Bcast!(collect(mesh.tree.length_level_0), mpi_root(), mpi_comm())[1]
    mesh.tree.periodicity    = Tuple(MPI.Bcast!(collect(mesh.tree.periodicity),    mpi_root(), mpi_comm()))

    # Set length
    n_cells = MPI.Bcast!(Ref(0), mpi_root(), mpi_comm())[]
    resize!(mesh.tree, n_cells)

    # Read in data
    @views MPI.Bcast!(mesh.tree.parent_ids[1:n_cells],      mpi_root(), mpi_comm())
    @views MPI.Bcast!(mesh.tree.child_ids[:, 1:n_cells],    mpi_root(), mpi_comm())
    @views MPI.Bcast!(mesh.tree.neighbor_ids[:, 1:n_cells], mpi_root(), mpi_comm())
    @views MPI.Bcast!(mesh.tree.levels[1:n_cells],          mpi_root(), mpi_comm())
    @views MPI.Bcast!(mesh.tree.coordinates[:, 1:n_cells],  mpi_root(), mpi_comm())
  end

  # Partition mesh
  partition!(mesh)

  return mesh
end

end # @muladd

ext: jl


text:
function show(io::IO, c::Colorant)
    show_colorant_string_with_eltype(io, typeof(c))
    _show_components(io, c)
end

# Special handling for Normed types: don't print the giant type name
function show(io::IO, c::ColorantNormed)
    show_colorant_string_with_eltype(io, typeof(c))
    if typeof(c) === base_colorant_type(typeof(c)) # non-parametric
        _show_components(io, c) # with trailing "N0f8" unless :compat=>true
    else
        _show_components(IOContext(io, :compact=>true), c)
    end
end

function _show_components(io::IO, c::Colorant{T,N}) where {T, N}
    print(io, '(')
    for i = 1:N
        i == 1 && show(io, comp1(c))
        i == 2 && show(io, comp2(c))
        i == 3 && show(io, comp3(c))
        i == 4 && show(io, comp4(c))
        i == 5 && show(io, comp5(c))
        print(io, i < N ? ',' : ')') # without spaces
    end
end

function Base.showarg(io::IO, a::Array{C}, toplevel) where C<:Colorant
    toplevel || print(io, "::")
    print(io, "Array{")
    show_colorant_string_with_eltype(io, C)
    print(io, ",$(ndims(a))}")
    toplevel && print(io, " with eltype ", C)
end

colorant_string(::Type{Union{}}) = "Union{}"
colorant_string(::Type{C}) where {C<:Colorant} = string(nameof(C))
function colorant_string_with_eltype(::Type{C}) where {C<:Colorant}
    io = IOBuffer()
    show_colorant_string_with_eltype(io, C)
    String(take!(io))
end

# for backward compatibility (Some other packages use this method.)
colorant_string_with_eltype(io::IO, ::Type{C}) where {C<:Colorant} =
    show_colorant_string_with_eltype(io, C)

show_colorant_string_with_eltype(io::IO, ::Type{Union{}}) = show(io, Union{})
function show_colorant_string_with_eltype(io::IO, ::Type{C}) where {C<:Colorant}
    if !isconcretetype(C)
        get(io, :compact, false) && return show(io, C)
        show(IOContext(io, :compact => true), C) # w/o module names
    elseif C === base_colorant_type(C) # non-parametric
        print(io, nameof(C))
    else
        print(io, nameof(C))
        print(io, '{')
        showcoloranttype(io, eltype(C))
        print(io, '}')
    end
end

showcoloranttype(io, ::Type{Union{}}) = show(io, Union{})
showcoloranttype(io, ::Type{T}) where {T<:FixedPoint} = FixedPointNumbers.showtype(io, T)
showcoloranttype(io, ::Type{T}) where {T} = show(io, T)


text:
using Test
using NeXLCore

@testset "MIP" begin
    @test all(map(z -> isapprox(J(Bloch1933, z) / z, 13.5, rtol = 1.0e-8), 1:99))
    @test all(
        map(
            z -> isapprox(
                J(Jensen1937, z) / z,
                9.0 * (1.0 + 0.5 * z^(-2.0 / 3.0)),
                rtol = 1.0e-8,
            ),
            1:99,
        ),
    )
    @test all(map(z -> isapprox(J(Wilson1941, z) / z, 11.5, rtol = 1.0e-8), 1:99))
    @test all(
        map(
            z -> isapprox(J(Sternheimer1964, z) / z, 9.76 + 58.82 * z^-1.19, rtol = 1.0e-8),
            12:99,
        ),
    )
    @test all(
        map(
            z -> isapprox(
                J(Springer1967, z) / z,
                9.0 * (1.0 + z^(-2.0 / 3.0)) + 0.03 * z,
                rtol = 1.0e-8,
            ),
            1:99,
        ),
    )
    @test all(
        map(
            z -> isapprox(
                J(Zeller1973, z) / z,
                10.04 + 8.25 * exp(-z / 11.22),
                rtol = 1.0e-8,
            ),
            1:99,
        ),
    )
    @test all(
        map(z -> isapprox(J(Brizuela1990, z) / z, 22.4 * z^-0.172, rtol = 1.0e-8), 1:99),
    )
    @test J(Berger1982, n"C") == 78.0
    @test J(Berger1982, n"Fe") == 286.0
    @test J(Berger1982, n"Pb") == 823.0
    @test J(Berger1982, n"Pu") == 921.0
end


text:
#################### Missing Values Sampler ####################

#################### Types and Constructors ####################

struct MISSTune
  dims::Tuple
  valueinds::Array
  distrinds::Array
end

function MISSTune(s::AbstractStochastic)
  MISSTune(s.distr, s.value)
end

function MISSTune(d::Distribution, v)
  MISSTune((), findall(isnan.(v)), Int[])
end

function MISSTune(D::Array{UnivariateDistribution}, v::Array)
  inds = findall(isnan.(v))
  MISSTune(dims(D), inds, inds)
end

function MISSTune(D::Array{MultivariateDistribution}, v::Array)
  isvalueinds = falses(size(v))
  isdistrinds = falses(size(D))
  for sub in CartesianIndices(size(D))
    n = length(D[sub])
    for i in 1:n
      if isnan(v[sub, i])
        isvalueinds[sub, i] = isdistrinds[sub] = true
      end
    end
  end
  MISSTune(dims(D), findall(isvalueinds), findall(isdistrinds))
end

#################### Sampler Constructor ####################
"""
    MISS(params::ElementOrVector{Symbol})

Construct a `Sampler` object to sampling missing output values. The constructor
should only be used to sample stochastic nodes upon which no other stochastic
node depends. So-called ‘output nodes’ can be identified with the `keys()`
function. Moreover, when the `MISS` constructor is included in a vector of
`Sampler` objects to define a sampling scheme, it should be positioned at the
beginning of the vector. This ensures that missing output values are updated
before any other samplers are executed.

Returns a `Sampler{Dict{Symbol, MISSTune}}` type object.

* `params`: stochastic node(s) that contain missing values (`NaN`) to be updated with the sampler.
"""
function MISS(params::ElementOrVector{Symbol})
  params = asvec(params)
  samplerfx = function(model::Model, block::Integer)
    tune = gettune(model, block)
    if model.iter == 1
      for key in params
        miss = MISSTune(model[key])
        if !isempty(miss.valueinds)
          tune[key] = miss
        end
      end
      params = intersect(keys(model, :dependent), keys(tune))
    end
    for key in params
      node = model[key]
      miss = tune[key]
      node[miss.valueinds] = rand(node, miss)
      update!(model, node.targets)
    end
    nothing
  end
  Sampler(params, samplerfx, Dict{Symbol, MISSTune}())
end

#################### Sampling Functions ####################

rand(s::AbstractStochastic, miss::MISSTune) = rand_sub(s.distr, miss)

function rand_sub(d::Distribution, miss::MISSTune)
  x = rand(d)
  Float64[x[i] for i in miss.valueinds]
end

function rand_sub(D::Array{UnivariateDistribution}, miss::MISSTune)
  Float64[rand(d) for d in D[miss.distrinds]]
end

function rand_sub(D::Array{MultivariateDistribution}, miss::MISSTune)
  X = Array{Float64}(undef, miss.dims)
  for i in miss.distrinds
    d = D[i]
    X[ind2sub(D, i)..., 1:length(d)] = rand(d)
  end
  X[miss.valueinds]
end


text:
struct UCI <: DataRepo
end

base_url(::UCI) = "https://archive.ics.uci.edu/ml/datasets/"

function to_cite(::UCI, mainpage)
    section = text_only(Gumbo.children(last(eachmatch(sel"p + p.normal", mainpage.root))))

    section = replace(section, "Available at: \n[Web Link]" => "")
    section = replace(section, r"\(?\s*\[Web Link\]\s*\)?" => "")
    section *= "\nLichman, M. (2013). UCI Machine Learning Repository [http://archive.ics.uci.edu/ml]. Irvine, CA: University of California, School of Information and Computer Science."

    replace(section, r"\n+" => "\n")
end

function description(::UCI, mainpage)
    desc =  text_only(first(eachmatch(sel"p.normal", mainpage.root)))
    desc * "\n\n" *  to_cite(UCI(), mainpage)
end

function data_fullname(::UCI, mainpage)
    data_fullname = text_only(eachmatch(sel".heading b", mainpage.root))
    data_fullname *= " (UCI ML Repository)"
end

function get_urls(::UCI, mainpage)
    datapage_link = first(eachmatch(sel"tr tr a", mainpage.root))
    datapage_url = joinpath(base_url(UCI()), attrs(datapage_link)["href"])

    data_urls = get_dataurls_from_webserver_index(datapage_url)
end


text:
# plotting.jl
using Gadfly
using Cairo
using Fontconfig
using ArrayViews

include("./b-splines.jl")

function evaluateBasis(U, u, order, nctl, concerned_index, N)

  active_bases = zeros(Float64, order)

  for i = 1:length(u)
    span = findSpan(u[i], nctl, U, order)
    basisFunctions(U, order, u[i], span, active_bases)
    for j = 1:length(active_bases)
      index = span - order + j
      if index == concerned_index
        N[i] = active_bases[j]
      end # end if index == concerned_index
    end # end for j = 1:length(active_bases)
  end  # end for i = 1:length(u)

  return nothing
end

U = [0,0,0,0, 0.1,0.2,0.5,1.0,1.0,1.0,1.0]
u = 0:0.01:1
order = 4
degree = order - 1
nctl = length(U) - order
nbases = length(U) - order # Total number of basis functions for a knot vector
# concerned_index = 3 # Index of the basis function to be plotted
N = zeros(Float64, length(u), nbases)
# println("sixe(N,2) = ", size(N,2))
for i = 1:size(N,2)
  concerned_index = i
  evaluateBasis(U, u, order, nctl, concerned_index,sview(N,:,i))
end

# println("N = \n$N")

myplot = plot(
  layer(x=u, y=N[:,1], Geom.line, Theme(default_color=colorant"green")),
  layer(x=u, y=N[:,2], Geom.line, Theme(default_color=colorant"blue")),
  layer(x=u, y=N[:,3], Geom.line, Theme(default_color=colorant"red")),
  layer(x=u, y=N[:,4], Geom.line, Theme(default_color=colorant"black")),
  layer(x=u, y=N[:,5], Geom.line, Theme(default_color=colorant"brown")),
  layer(x=u, y=N[:,6], Geom.line, Theme(default_color=colorant"pink")),
  layer(x=u, y=N[:,7], Geom.line, Theme(default_color=colorant"yellow")),
  Guide.XLabel("u"),
  Guide.YLabel("Basis Function Values"),
  Guide.yticks(ticks=collect(0:0.2:1))
)

draw(PDF("Bases.pdf", 10inch, 10inch), myplot)


text:
"""
## mocouq

This subroutine forms the derivatives of a Mohr-Coulomb potential
function with respect to the three invariants (psi in degrees).

### Function
```julia
(dq1,dq2,dq3) = mocouq(psi,dsbar,theta)
```

### Arguments
```julia
* psi::Float64               : Local force and momemts (Updated)
* dsbar::Float64             : Globale forces and moments
* theta::Float64             : Element orientation angle (3D)
```

### Return values
```julia
* dq1::Float64               : Local force and momemts (Updated)
* dq2::Float64               : Globale forces and moments
* dq3::Float64               : Element orientation angle (3D)
```
"""
function mocouq(psi,dsbar,theta)
 snth=sin(theta)
 snps=sind(psi)
 sq3=sqrt(3.0)
 dq1=snps
 if abs(snth) > 0.49
   c1 = 1.0
   snth < 0.0 && (c1 = -1.0)
   dq2 = (sq3*0.5-c1*snps*0.5/sq3)*sq3*0.5/dsbar
   dq3 = 0.0
 else
   csth = cos(theta)
   cs3th = cos(3theta)
   tn3th = tan(3theta)
   tnth = snth/csth
   dq2 = sq3*csth/dsbar*((1.0+tnth*tn3th)+snps*(tn3th-tnth)/sq3)*0.5
   dq3 = 0.5*3.0*(sq3*snth+snps*csth)/(cs3th*dsbar*dsbar)
 end
 (dq1, dq2, dq3)
end


text:
"""
    KalmanProcessUtils

Utils for the construction and handling of Kalman Process structs.
"""
module KalmanProcessUtils

using LinearAlgebra
using Statistics
using JLD2
using DocStringExtensions
using EnsembleKalmanProcesses
import EnsembleKalmanProcesses: Process, Unscented, Inversion, Sampler, SparseInversion
import EnsembleKalmanProcesses: SampleSuccGauss, IgnoreFailures
using EnsembleKalmanProcesses.ParameterDistributions

using ..ReferenceModels
using ..ReferenceStats
using ..ModelTypes
using ..LESUtils
using ..HelperFuncs
using ..DistributionUtils

export generate_ekp, generate_tekp
export get_sparse_indices, get_regularized_indices
export get_Δt, PiecewiseConstantDecay, PiecewiseConstantGrowth

abstract type LearningRateScheduler end

"""
    PiecewiseConstantDecay{FT <: Real, IT <: Int} <: LearningRateScheduler

Piecewise constant decay learning rate scheduler.

Halves the time step periodically with period `τ`.

# Fields

$(TYPEDFIELDS)

# Constructors

$(METHODLIST)
"""
struct PiecewiseConstantDecay{FT <: Real, IT <: Int} <: LearningRateScheduler
    "Initial learning rate"
    Δt_init::FT
    "Halving time"
    τ::IT
end

"""
    PiecewiseConstantGrowth{FT <: Real, IT <: Int} <: LearningRateScheduler

Piecewise constant growth learning rate scheduler.

Doubles the time step periodically with period `τ`.

# Fields

$(TYPEDFIELDS)

# Constructors

$(METHODLIST)
"""
struct PiecewiseConstantGrowth{FT <: Real, IT <: Int} <: LearningRateScheduler
    "Initial learning rate"
    Δt_init::FT
    "Doubling time"
    τ::IT
end

"Retrieve learning rate `Δt` from a LearningRateScheduler"
get_Δt(Δt::FT, iteration::IT) where {FT <: Real, IT <: Int} = Δt
get_Δt(lrs::PiecewiseConstantDecay, iteration::IT) where {IT <: Int} = lrs.Δt_init * 2^(-floor(iteration / lrs.τ))
get_Δt(lrs::PiecewiseConstantGrowth, iteration::IT) where {IT <: Int} = lrs.Δt_init * 2^(floor(iteration / lrs.τ))

"""
    generate_ekp(
        ref_stats::ReferenceStatistics,
        process::Process,
        u::Union{Matrix{T}, T} = nothing;
        failure_handler::String = "ignore_failures",
        outdir_path::String = pwd(),
        to_file::Bool = true,
    ) where {T}

Generates, and possible writes to file, an EnsembleKalmanProcess
from a parameter ensemble and reference statistics.

Inputs:
 - ref_stats :: ReferenceStatistics defining the inverse problem.
 - process :: Type of EnsembleKalmanProcess used to evolve the ensemble.
 - u :: An ensemble of parameter vectors, used if !isa(process, Unscented).
 - failure_handler :: String describing what failure handler to use.
 - outdir_path :: Output path.
 - to_file :: Whether to write the serialized prior to a JLD2 file.

Output:
 - The generated EnsembleKalmanProcess.
"""
function generate_ekp(
    ref_stats::ReferenceStatistics,
    process::Process,
    u::Union{Matrix{T}, T} = nothing;
    failure_handler::String = "ignore_failures",
    outdir_path::String = pwd(),
    to_file::Bool = true,
) where {T}

    @assert isa(process, Unscented) || !isnothing(u) "Incorrect EKP constructor."
    @assert failure_handler in ["ignore_failures", "high_loss", "sample_succ_gauss"]
    if failure_handler == "sample_succ_gauss"
        fh = SampleSuccGauss()
    else
        fh = IgnoreFailures()
    end
    ekp =
        isnothing(u) ? EnsembleKalmanProcess(ref_stats.y, ref_stats.Γ, process, failure_handler_method = fh) :
        EnsembleKalmanProcess(u, ref_stats.y, ref_stats.Γ, process, failure_handler_method = fh)
    if to_file
        jldsave(ekobj_path(outdir_path, 1); ekp)
    end
    return ekp
end

"""
    generate_tekp(
        ref_stats::ReferenceStatistics,
        priors::ParameterDistribution,
        process::Process,
        u::Union{Matrix{T}, T} = nothing;
        l2_reg::Union{Dict{String, Vector{R}}, R} = nothing,
        failure_handler::String = "ignore_failures",
        outdir_path::String = pwd(),
        to_file::Bool = true,
    ) where {T, R}

Generates, and possible writes to file, a Tikhonov EnsembleKalmanProcess
from a parameter ensemble and reference statistics.

Tikhonov regularization is implemented through output state augmentation
with the input parameter vector. The input L2 regularization hyperparameter
should be interpreted as the inverse of the variance of our prior belief in
the magnitude of the parameters.

Inputs:
 - ref_stats :: ReferenceStatistics defining the inverse problem.
 - priors :: Parameter priors used for L2 (i.e., Tikhonov) regularization
 - process :: Type of EnsembleKalmanProcess used to evolve the ensemble.
 - u :: An ensemble of parameter vectors, used if !isa(process, Unscented).
 - l2_reg :: L2 regularization hyperparameter driving parameter values toward prior.
        May be a float (isotropic regularization) or a dictionary of regularizations
        per parameter.
 - failure_handler :: String describing what failure handler to use.
 - outdir_path :: Output path.
 - to_file :: Whether to write the serialized prior to a JLD2 file.

Output:
 - The generated augmented EnsembleKalmanProcess.
"""
function generate_tekp(
    ref_stats::ReferenceStatistics,
    priors::ParameterDistribution,
    process::Process,
    u::Union{Matrix{T}, T} = nothing;
    l2_reg::Union{Dict{String, Vector{R}}, R} = nothing,
    failure_handler::String = "ignore_failures",
    outdir_path::String = pwd(),
    to_file::Bool = true,
) where {T, R}

    @assert isa(process, Unscented) || !isnothing(u) "Incorrect TEKP constructor."
    @assert failure_handler in ["ignore_failures", "high_loss", "sample_succ_gauss"]
    if failure_handler == "sample_succ_gauss"
        fh = SampleSuccGauss()
    else
        fh = IgnoreFailures()
    end

    μ = vcat(mean(priors)...)
    if isa(l2_reg, Dict)
        # flatten l2_reg dict
        if any(1 .< [length(val) for val in collect(values(l2_reg))])
            _, l2_reg_values = flatten_config_dict(l2_reg)
        else
            l2_reg_values = collect(values(l2_reg))
        end
        l2_reg_values = vcat(l2_reg_values...)

        # dict must be complete to preserve ordering
        @assert length(μ) == length(l2_reg_values) "Dictionary of regularizations l2_reg must include all parameters."

        # Augment exclusively with nonzero weights
        aug_indices = get_regularized_indices(l2_reg)
        l2_reg_values = l2_reg_values[aug_indices]
        μ = μ[aug_indices]
        Γ_θ = inv(Diagonal(l2_reg_values))

    elseif !isnothing(l2_reg)
        @assert l2_reg > eps(R) "If system is augmented, provide nonzero l2_reg."
        Γ_θ = Diagonal(repeat([inv(l2_reg)], length(μ)))
    else
        Γ_θ = cov(priors)
    end

    # Augment system with regularization towards prior mean
    y_aug = vcat([ref_stats.y, μ]...)
    Γ_aug_list = [ref_stats.Γ, Array(Γ_θ)]
    Γ_aug = cat(Γ_aug_list..., dims = (1, 2))

    ekp =
        isnothing(u) ? EnsembleKalmanProcess(y_aug, Γ_aug, process, failure_handler_method = fh) :
        EnsembleKalmanProcess(u, y_aug, Γ_aug, process, failure_handler_method = fh)
    if to_file
        jldsave(ekobj_path(outdir_path, 1); ekp)
    end
    return ekp
end

"Returns the sparse parameter indices given the sparsity configuration and the number of parameters."
function get_sparse_indices(sparsity_config, p)
    if isa(sparsity_config, Dict{String, Vector{Bool}})
        @assert sum(map(length, collect(values(sparsity_config)))) == p
        return flat_dict_keys_where(sparsity_config)
    elseif isa(sparsity_config, Bool) && sparsity_config
        return Colon()
    else
        throw(ArgumentError("Sparsity config entry not recognized, pass a Bool or Dict{String, Vector{Bool}}."))
    end
end

"Returns the indices of parameters to be regularized, given the l2 regularization configuration dictionary."
get_regularized_indices(l2_config::Dict) = flat_dict_keys_where(l2_config, above_eps)

end # module


text:
module SCIP

# assorted utility functions
include("util.jl")

# load deps, version check
include("init.jl")

# wrapper of SCIP library
include("wrapper.jl")

# memory management
include("managed_scip.jl")

# constraints from nonlinear expressions
include("nonlinear.jl")

# constraint handlers
include("conshdlr.jl")

# separators
include("sepa.jl")

# implementation of MOI
include("MOI_wrapper.jl")

# convenience functions
include("convenience.jl")

# warn about rewrite
include("compat.jl")

end


text:
# store directory with trailing slash to ensure correct order
# https://github.com/mirage/irmin/issues/352

struct GitMutableTree
    entries::SortedDict{String,Tuple{GitMode,Any},Base.Order.ForwardOrdering}
    GitMutableTree() = new(SortedDict{String,Tuple{GitMode,Any}}(Base.Order.Forward))
end

function Base.setindex!(mtree::GitMutableTree, value, path::String)
    n = findfirst(isequal('/'), path)
    if n === nothing
        # file
        mtree.entries[path] = value
    else
        base = path[1:n]
        rpath = path[n+1:end]
        # waiting on:
        #   https://github.com/JuliaCollections/DataStructures.jl/pull/336
        # rmode, rtree = get!(() -> (mode_dir, GitMutableTree()), mtree.entries, base)
        if !haskey(mtree.entries, base)
            mtree.entries[base] = (mode_dir, GitMutableTree())
        end
        rmode, rtree = mtree.entries[base]
        rtree[rpath] = value
    end
end

function Base.convert(::Type{GitTree}, mtree::GitMutableTree)
    entries = GitTreeEntry[]
    for (name, (mode, obj)) in mtree.entries
        if endswith(name, '/')
            name = chop(name)
        end
        if obj isa SHA1
            hash = obj
        else
            hash = oid(convert(GitTree, obj))
        end
        push!(entries,GitTreeEntry(mode, name, hash))
    end
    GitTree(entries)
end

oid(mtree::GitMutableTree) = oid(convert(GitTree, mtree))

function treehash(filelist)
    mtree = GitMutableTree()
    for filename in filelist
        # TODO: deal with executable & symlinks
        hash = oid(filename)
        mtree[filename] = (mode_normal, hash)
    end
    return oid(mtree)
end

max_stars_repo_name: simonbyrne/GitX

max_issues_repo_name: simonbyrne/GitX

max_forks_repo_name: simonbyrne/GitX


text:
# Adapted from libtiff v4.0.9

#   Copyright (c) 1988-1997 Sam Leffler
#   Copyright (c) 1991-1997 Silicon Graphics, Inc.

#   Permission to use, copy, modify, distribute, and sell this software and
#   its documentation for any purpose is hereby granted without fee, provided
#   that (i) the above copyright notices and this permission notice appear in
#   all copies of the software and related documentation, and (ii) the names of
#   Sam Leffler and Silicon Graphics may not be used in any advertising or
#   publicity relating to the software without the specific, prior written
#   permission of Sam Leffler and Silicon Graphics.

#   THE SOFTWARE IS PROVIDED "AS-IS" AND WITHOUT WARRANTY OF ANY KIND,
#   EXPRESS, IMPLIED OR OTHERWISE, INCLUDING WITHOUT LIMITATION, ANY
#   WARRANTY OF MERCHANTABILITY OR FITNESS FOR A PARTICULAR PURPOSE.

#   IN NO EVENT SHALL SAM LEFFLER OR SILICON GRAPHICS BE LIABLE FOR
#   ANY SPECIAL, INCIDENTAL, INDIRECT OR CONSEQUENTIAL DAMAGES OF ANY KIND,
#   OR ANY DAMAGES WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS,
#   WHETHER OR NOT ADVISED OF THE POSSIBILITY OF DAMAGE, AND ON ANY THEORY OF
#   LIABILITY, ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE
#   OF THIS SOFTWARE.

#   NB: In the comments below,
#    - items marked with a + are obsoleted by revision 5.0,
#    - items marked with a ! are introduced in revision 6.0.
#    - items marked with a % are introduced post revision 6.0.
#    - items marked with a $ are obsoleted by revision 6.0.
#    - items marked with a & are introduced by Adobe DNG specification.

"""
	$(TYPEDEF)

List of many common named TIFF Tags. This is not an exhaustive list but should
cover most cases.

$(FIELDS)

"""
@enum TiffTag begin
	SUBFILETYPE = 254 			 #  subfile data descriptor
#define	    FILETYPE_REDUCEDIMAGE	0x1	/* reduced resolution version */
#define	    FILETYPE_PAGE		0x2	/* one page of many */
#define	    FILETYPE_MASK		0x4	/* transparency mask */
	OSUBFILETYPE = 255 			 #  +kind of data in subfile
#define	    OFILETYPE_IMAGE		1	/* full resolution image data */
#define	    OFILETYPE_REDUCEDIMAGE	2	/* reduced size image data */
#define	    OFILETYPE_PAGE		3	/* one page of many */
	IMAGEWIDTH = 256 			 #  image width in pixels
	IMAGELENGTH = 257 			 #  image height in pixels
	BITSPERSAMPLE = 258 			 #  bits per channel (sample)
	COMPRESSION = 259 			 #  data compression technique
	PHOTOMETRIC = 262 			 #  photometric interpretation
	THRESHHOLDING = 263 			 #  +thresholding used on data
#define	    THRESHHOLD_BILEVEL		1	/* b&w art scan */
#define	    THRESHHOLD_HALFTONE		2	/* or dithered scan */
#define	    THRESHHOLD_ERRORDIFFUSE	3	/* usually floyd-steinberg */
	CELLWIDTH = 264 			 #  +dithering matrix width
	CELLLENGTH = 265 			 #  +dithering matrix height
	FILLORDER = 266 			 #  data order within a byte
#define	    FILLORDER_MSB2LSB		1	/* most significant -> least */
#define	    FILLORDER_LSB2MSB		2	/* least significant -> most */
	DOCUMENTNAME = 269 			 #  name of doc. image is from
	IMAGEDESCRIPTION = 270 			 #  info about image
	MAKE = 271 			 #  scanner manufacturer name
	MODEL = 272 			 #  scanner model name/number
	STRIPOFFSETS = 273 			 #  offsets to data strips
	ORIENTATION = 274 			 #  +image orientation
#define	    ORIENTATION_TOPLEFT		1	/* row 0 top, col 0 lhs */
#define	    ORIENTATION_TOPRIGHT	2	/* row 0 top, col 0 rhs */
#define	    ORIENTATION_BOTRIGHT	3	/* row 0 bottom, col 0 rhs */
#define	    ORIENTATION_BOTLEFT		4	/* row 0 bottom, col 0 lhs */
#define	    ORIENTATION_LEFTTOP		5	/* row 0 lhs, col 0 top */
#define	    ORIENTATION_RIGHTTOP	6	/* row 0 rhs, col 0 top */
#define	    ORIENTATION_RIGHTBOT	7	/* row 0 rhs, col 0 bottom */
#define	    ORIENTATION_LEFTBOT		8	/* row 0 lhs, col 0 bottom */
	SAMPLESPERPIXEL = 277 			 #  samples per pixel
	ROWSPERSTRIP = 278 			 #  rows per strip of data
	STRIPBYTECOUNTS = 279 			 #  bytes counts for strips
	MINSAMPLEVALUE = 280 			 #  +minimum sample value
	MAXSAMPLEVALUE = 281 			 #  +maximum sample value
	XRESOLUTION = 282 			 #  pixels/resolution in x
	YRESOLUTION = 283 			 #  pixels/resolution in y
	PLANARCONFIG = 284 			 #  storage organization
#define	    PLANARCONFIG_CONTIG		1	/* single image plane */
#define	    PLANARCONFIG_SEPARATE	2	/* separate planes of data */
	PAGENAME = 285 			 #  page name image is from
	XPOSITION = 286 			 #  x page offset of image lhs
	YPOSITION = 287 			 #  y page offset of image lhs
	FREEOFFSETS = 288 			 #  +byte offset to free block
	FREEBYTECOUNTS = 289 			 #  +sizes of free blocks
	GRAYRESPONSEUNIT = 290 			 #  $gray scale curve accuracy
#define	    GRAYRESPONSEUNIT_10S	1	/* tenths of a unit */
#define	    GRAYRESPONSEUNIT_100S	2	/* hundredths of a unit */
#define	    GRAYRESPONSEUNIT_1000S	3	/* thousandths of a unit */
#define	    GRAYRESPONSEUNIT_10000S	4	/* ten-thousandths of a unit */
#define	    GRAYRESPONSEUNIT_100000S	5	/* hundred-thousandths */
	GRAYRESPONSECURVE = 291 			 #  $gray scale response curve
	T4OPTIONS = 292 			 #  TIFF 6.0 proper name alias
#define	    GROUP3OPT_2DENCODING	0x1	/* 2-dimensional coding */
#define	    GROUP3OPT_UNCOMPRESSED	0x2	/* data not compressed */
#define	    GROUP3OPT_FILLBITS		0x4	/* fill to byte boundary */
	T6OPTIONS = 293 			 #  TIFF 6.0 proper name
#define	    GROUP4OPT_UNCOMPRESSED	0x2	/* data not compressed */
	RESOLUTIONUNIT = 296 			 #  units of resolutions
#define	    RESUNIT_NONE		1	/* no meaningful units */
#define	    RESUNIT_INCH		2	/* english */
#define	    RESUNIT_CENTIMETER		3	/* metric */
	PAGENUMBER = 297 			 #  page numbers of multi-page
	COLORRESPONSEUNIT = 300 			 #  $color curve accuracy
#define	    COLORRESPONSEUNIT_10S	1	/* tenths of a unit */
#define	    COLORRESPONSEUNIT_100S	2	/* hundredths of a unit */
#define	    COLORRESPONSEUNIT_1000S	3	/* thousandths of a unit */
#define	    COLORRESPONSEUNIT_10000S	4	/* ten-thousandths of a unit */
#define	    COLORRESPONSEUNIT_100000S	5	/* hundred-thousandths */
	TRANSFERFUNCTION = 301 			 #  !colorimetry info
	SOFTWARE = 305 			 #  name & release
	DATETIME = 306 			 #  creation date and time
	ARTIST = 315 			 #  creator of image
	HOSTCOMPUTER = 316 			 #  machine where created
	PREDICTOR = 317 			 #  prediction scheme w/ LZW
#define     PREDICTOR_NONE		1	/* no prediction scheme used */
#define     PREDICTOR_HORIZONTAL	2	/* horizontal differencing */
#define     PREDICTOR_FLOATINGPOINT	3	/* floating point predictor */
	WHITEPOINT = 318 			 #  image white point
	PRIMARYCHROMATICITIES = 319 			 #  !primary chromaticities
	COLORMAP = 320 			 #  RGB map for palette image
	HALFTONEHINTS = 321 			 #  !highlight+shadow info
	TILEWIDTH = 322 			 #  !tile width in pixels
	TILELENGTH = 323 			 #  !tile height in pixels
	TILEOFFSETS = 324 			 #  !offsets to data tiles
	TILEBYTECOUNTS = 325 			 #  !byte counts for tiles
	BADFAXLINES = 326 			 #  lines w/ wrong pixel count
	CLEANFAXDATA = 327 			 #  regenerated line info
#define	    CLEANFAXDATA_CLEAN		0	/* no errors detected */
#define	    CLEANFAXDATA_REGENERATED	1	/* receiver regenerated lines */
#define	    CLEANFAXDATA_UNCLEAN	2	/* uncorrected errors exist */
	CONSECUTIVEBADFAXLINES = 328 			 #  max consecutive bad lines
	SUBIFD = 330 			 #  subimage descriptors
	INKSET = 332 			 #  !inks in separated image
#define	    INKSET_CMYK			1	/* !cyan-magenta-yellow-black color */
#define	    INKSET_MULTIINK		2	/* !multi-ink or hi-fi color */
	INKNAMES = 333 			 #  !ascii names of inks
	NUMBEROFINKS = 334 			 #  !number of inks
	DOTRANGE = 336 			 #  !0% and 100% dot codes
	TARGETPRINTER = 337 			 #  !separation target
	EXTRASAMPLES = 338 			 #  !info about extra samples
#define	    EXTRASAMPLE_UNSPECIFIED	0	/* !unspecified data */
#define	    EXTRASAMPLE_ASSOCALPHA	1	/* !associated alpha data */
#define	    EXTRASAMPLE_UNASSALPHA	2	/* !unassociated alpha data */
	SAMPLEFORMAT = 339 			 #  !data sample format
#define	    SAMPLEFORMAT_UINT		1	/* !unsigned integer data */
#define	    SAMPLEFORMAT_INT		2	/* !signed integer data */
#define	    SAMPLEFORMAT_IEEEFP		3	/* !IEEE floating point data */
#define	    SAMPLEFORMAT_VOID		4	/* !untyped data */
#define	    SAMPLEFORMAT_COMPLEXINT	5	/* !complex signed int */
#define	    SAMPLEFORMAT_COMPLEXIEEEFP	6	/* !complex ieee floating */
	SMINSAMPLEVALUE = 340 			 #  !variable MinSampleValue
	SMAXSAMPLEVALUE = 341 			 #  !variable MaxSampleValue
    CLIPPATH = 343           # %ClipPath [Adobe TIFF technote 2]
    XCLIPPATHUNITS = 344     # %XClipPathUnits [Adobe TIFF technote 2]
    YCLIPPATHUNITS = 345	 # %YClipPathUnits [Adobe TIFF technote 2]
    INDEXED	= 346        # %Indexed [Adobe TIFF Technote 3]
	JPEGTABLES = 347 			 #  %JPEG table stream
	OPIPROXY = 351 			 #  %OPI Proxy [Adobe TIFF technote]
# Tags 400-435 are from the TIFF/FX spec
	GLOBALPARAMETERSIFD = 400 			 #  !
	PROFILETYPE = 401 			 #  !
#define     PROFILETYPE_UNSPECIFIED	0	/* ! */
#define     PROFILETYPE_G3_FAX		1	/* ! */
	FAXPROFILE = 402 			 #  !
#define     FAXPROFILE_S			1	/* !TIFF/FX FAX profile S */
#define     FAXPROFILE_F			2	/* !TIFF/FX FAX profile F */
#define     FAXPROFILE_J			3	/* !TIFF/FX FAX profile J */
#define     FAXPROFILE_C			4	/* !TIFF/FX FAX profile C */
#define     FAXPROFILE_L			5	/* !TIFF/FX FAX profile L */
#define     FAXPROFILE_M			6	/* !TIFF/FX FAX profile LM */
	CODINGMETHODS = 403 			 #  !TIFF/FX coding methods
#define     CODINGMETHODS_T4_1D		(1 << 1)	/* !T.4 1D */
#define     CODINGMETHODS_T4_2D		(1 << 2)	/* !T.4 2D */
#define     CODINGMETHODS_T6		(1 << 3)	/* !T.6 */
#define     CODINGMETHODS_T85 		(1 << 4)	/* !T.85 JBIG */
#define     CODINGMETHODS_T42 		(1 << 5)	/* !T.42 JPEG */
#define     CODINGMETHODS_T43		(1 << 6)	/* !T.43 colour by layered JBIG */
	VERSIONYEAR = 404 			 #  !TIFF/FX version year
	MODENUMBER = 405 			 #  !TIFF/FX mode number
	DECODE = 433 			 #  !TIFF/FX decode
	IMAGEBASECOLOR = 434 			 #  !TIFF/FX image base colour
	T82OPTIONS = 435 			 #  !TIFF/FX T.82 options

# Tags 512-521 are obsoleted by Technical Note #2 which specifies a
# revised JPEG-in-TIFF scheme.

	JPEGPROC = 512 			 #  !JPEG processing algorithm
#define	    JPEGPROC_BASELINE		1	/* !baseline sequential */
#define	    JPEGPROC_LOSSLESS		14	/* !Huffman coded lossless */
	JPEGIFOFFSET = 513 			 #  !pointer to SOI marker
	JPEGIFBYTECOUNT = 514 			 #  !JFIF stream length
	JPEGRESTARTINTERVAL = 515 			 #  !restart interval length
	JPEGLOSSLESSPREDICTORS = 517 			 #  !lossless proc predictor
	JPEGPOINTTRANSFORM = 518 			 #  !lossless point transform
	JPEGQTABLES = 519 			 #  !Q matrix offsets
	JPEGDCTABLES = 520 			 #  !DCT table offsets
	JPEGACTABLES = 521 			 #  !AC coefficient offsets
	YCBCRCOEFFICIENTS = 529 			 #  !RGB -> YCbCr transform
	YCBCRSUBSAMPLING = 530 			 #  !YCbCr subsampling factors
	YCBCRPOSITIONING = 531 			 #  !subsample positioning
#define	    YCBCRPOSITION_CENTERED	1	/* !as in PostScript Level 2 */
#define	    YCBCRPOSITION_COSITED	2	/* !as in CCIR 601-1 */
	REFERENCEBLACKWHITE = 532 			 #  !colorimetry info
	STRIPROWCOUNTS = 559 			 #  !TIFF/FX strip row counts
	XMLPACKET = 700	         # %XML packet [Adobe XMP Specification, January 2004]
    OPIIMAGEID = 32781	     # %OPI ImageID [Adobe TIFF technote]
# tags 32952-32956 are private tags registered to Island Graphics
	REFPTS = 32953 			 #  image reference points
	REGIONTACKPOINT = 32954 			 #  region-xform tack point
	REGIONWARPCORNERS = 32955 			 #  warp quadrilateral
	REGIONAFFINE = 32956 			 #  affine transformation mat
# tags 32995-32999 are private tags registered to SGI
	MATTEING = 32995 			 #  $use ExtraSamples
	DATATYPE = 32996 			 #  $use SampleFormat
	IMAGEDEPTH = 32997 			 #  z depth of image
	TILEDEPTH = 32998 			 #  z depth/data tile
# tags 33300-33309 are private tags registered to Pixar

# TIFFTAG_PIXAR_IMAGEFULLWIDTH and TIFFTAG_PIXAR_IMAGEFULLLENGTH
# are set when an image has been cropped out of a larger image.
# They reflect the size of the original uncropped image.
# The TIFFTAG_XPOSITION and TIFFTAG_YPOSITION can be used
# to determine the position of the smaller image in the larger one.

    PIXAR_IMAGEFULLWIDTH = 33300     # full image size in x
    PIXAR_IMAGEFULLLENGTH = 33301    # full image size in y

# tag 33405 is a private tag registered to Eastman Kodak
	WRITERSERIALNUMBER = 33405 			 #  device serial number
	CFAREPEATPATTERNDIM = 33421 			 #  dimensions of CFA pattern
	CFAPATTERN = 33422 			 #  color filter array pattern
# tag 33432 is listed in the 6.0 spec w/ unknown ownership
	COPYRIGHT = 33432 			 #  copyright string
# IPTC TAG from RichTIFF specifications
    RICHTIFFIPTC = 33723
# 34016-34029 are reserved for ANSI IT8 TIFF/IT <dkelly@apago.com)
	IT8SITE = 34016 			 #  site name
	IT8COLORSEQUENCE = 34017 			 #  color seq. [RGB,CMYK,etc]
	IT8HEADER = 34018 			 #  DDES Header
	IT8RASTERPADDING = 34019 			 #  raster scanline padding
	IT8BITSPERRUNLENGTH = 34020 			 #  # of bits in short run
	IT8BITSPEREXTENDEDRUNLENGTH = 34021 			 #  # of bits in long run
	IT8COLORTABLE = 34022 			 #  LW colortable
	IT8IMAGECOLORINDICATOR = 34023 			 #  BP/BL image color switch
	IT8BKGCOLORINDICATOR = 34024 			 #  BP/BL bg color switch
	IT8IMAGECOLORVALUE = 34025 			 #  BP/BL image color value
	IT8BKGCOLORVALUE = 34026 			 #  BP/BL bg color value
	IT8PIXELINTENSITYRANGE = 34027 			 #  MP pixel intensity value
	IT8TRANSPARENCYINDICATOR = 34028 			 #  HC transparency switch
	IT8COLORCHARACTERIZATION = 34029 			 #  color character. table
	IT8HCUSAGE = 34030 			 #  HC usage indicator
    IT8TRAPINDICATOR = 34031	 # Trapping indicator (untrapped=0, trapped=1)
	IT8CMYKEQUIVALENT = 34032 			 #  CMYK color equivalents
# tags 34232-34236 are private tags registered to Texas Instruments
	FRAMECOUNT = 34232 			 #  Sequence Frame Count
# tag 34377 is private tag registered to Adobe for PhotoShop
    PHOTOSHOP = 34377
# tags 34665, 34853 and 40965 are documented in EXIF specification
	EXIFIFD = 34665 			 #  Pointer to EXIF private directory
# tag 34750 is a private tag registered to Adobe?
	ICCPROFILE = 34675 			 #  ICC profile data
	IMAGELAYER = 34732 			 #  !TIFF/FX image layer information
# tag 34750 is a private tag registered to Pixel Magic
	JBIGOPTIONS = 34750 			 #  JBIG options
	GPSIFD = 34853 			 #  Pointer to GPS private directory
# tags 34908-34914 are private tags registered to SGI
	FAXRECVPARAMS = 34908 			 #  encoded Class 2 ses. parms
	FAXSUBADDRESS = 34909 			 #  received SubAddr string
	FAXRECVTIME = 34910 			 #  receive time (secs)
	FAXDCS = 34911 			 #  encoded fax ses. params, Table 2/T.30
# tags 37439-37443 are registered to SGI <gregl@sgi.com>
	STONITS = 37439 			 #  Sample value to Nits
# tag 34929 is a private tag registered to FedEx
    FEDEX_EDR = 34929	# unknown use
	INTEROPERABILITYIFD = 40965 			 #  Pointer to Interoperability private directory
# Adobe Digital Negative (DNG) format tags */
	DNGVERSION = 50706 			 #  &DNG version number
	DNGBACKWARDVERSION = 50707 			 #  &DNG compatibility version
	UNIQUECAMERAMODEL = 50708 			 #  &name for the camera model
	CFALAYOUT = 50711 			 #  &spatial layout of the CFA
	LINEARIZATIONTABLE = 50712 			 #  &lookup table description
	BLACKLEVEL = 50714 			 #  &zero light encoding level
	DEFAULTSCALE = 50718 			 #  &default scale factors
	CAMERACALIBRATION1 = 50723 			 #  &calibration matrix 1
	CAMERACALIBRATION2 = 50724 			 #  &calibration matrix 2
	BASELINENOISE = 50731 			 #  &relative noise level
	LINEARRESPONSELIMIT = 50734 			 #  &non-linear encoding range
	CAMERASERIALNUMBER = 50735 			 #  &camera's serial number
	LENSINFO = 50736 			 #  info about the lens
	CHROMABLURRADIUS = 50737 			 #  &chroma blur radius
#define TIFFTAG_ANTIALIASSTRENGTH	50738	&relative strength of the camera's anti-alias filter
	SHADOWSCALE = 50739 			 #  &used by Adobe Camera Raw
	DNGPRIVATEDATA = 50740 			 #  &manufacturer's private data
#define TIFFTAG_MAKERNOTESAFETY		50741	whether the EXIF MakerNote tag is safe to preserve along with the rest of the EXIF data
	CALIBRATIONILLUMINANT1 = 50778 			 #  &illuminant 1
	CALIBRATIONILLUMINANT2 = 50779 			 #  &illuminant 2
	BESTQUALITYSCALE = 50780 			 #  &best quality multiplier
#define TIFFTAG_RAWDATAUNIQUEID		50781	/* &unique identifier for
#define TIFFTAG_ORIGINALRAWFILENAME	50827	/* &file name of the original
#define TIFFTAG_ORIGINALRAWFILEDATA	50828	/* &contents of the original
#define TIFFTAG_ACTIVEAREA		50829	/* &active (non-masked) pixels
#define TIFFTAG_MASKEDAREAS		50830	/* &list of coordinates
	ASSHOTICCPROFILE = 50831 			 #  &these two tags used to
#define TIFFTAG_ASSHOTPREPROFILEMATRIX	50832	/* map cameras's color space
	CURRENTICCPROFILE = 50833 			 #  &
	CURRENTPREPROFILEMATRIX = 50834 			 #  &
# tag 65535 is an undefined tag used by Eastman Kodak
	DCSHUESHIFTVALUES = 65535 			 #  hue shift correction data

# The following are ``pseudo tags'' that can be used to control
# codec-specific functionality.  These tags are not written to file.
# Note that these values start at 0xffff+1 so that they'll never
# collide with Aldus-assigned tags.

# If you want your private pseudo tags ``registered'' (i.e. added to
# this file), please post a bug report via the tracking system at
# http://www.remotesensing.org/libtiff/bugs.html with the appropriate
# C definitions to add.

	FAXMODE = 65536 			 #  Group 3/4 format control
#define	    FAXMODE_CLASSIC	0x0000		/* default, include RTC */
#define	    FAXMODE_NORTC	0x0001		/* no RTC at end of data */
#define	    FAXMODE_NOEOL	0x0002		/* no EOL code at end of row */
#define	    FAXMODE_BYTEALIGN	0x0004		/* byte align row */
#define	    FAXMODE_WORDALIGN	0x0008		/* word align row */
#define	    FAXMODE_CLASSF	FAXMODE_NORTC	/* TIFF Class F */
	JPEGQUALITY = 65537 			 #  Compression quality level
# Note: quality level is on the IJG 0-100 scale.  Default value is 75
	JPEGCOLORMODE = 65538 			 #  Auto RGB<=>YCbCr convert?
#define	    JPEGCOLORMODE_RAW	0x0000		/* no conversion (default) */
#define	    JPEGCOLORMODE_RGB	0x0001		/* do auto conversion */
	JPEGTABLESMODE = 65539 			 #  What to put in JPEGTables
#define	    JPEGTABLESMODE_QUANT 0x0001		/* include quantization tbls */
#define	    JPEGTABLESMODE_HUFF	0x0002		/* include Huffman tbls */
# Note: default is JPEGTABLESMODE_QUANT | JPEGTABLESMODE_HUFF */
	FAXFILLFUNC = 65540 			 #  G3/G4 fill function
	PIXARLOGDATAFMT = 65549 			 #  PixarLogCodec I/O data sz
#define	    PIXARLOGDATAFMT_8BIT	0	/* regular u_char samples */
#define	    PIXARLOGDATAFMT_8BITABGR	1	/* ABGR-order u_chars */
#define	    PIXARLOGDATAFMT_11BITLOG	2	/* 11-bit log-encoded (raw) */
#define	    PIXARLOGDATAFMT_12BITPICIO	3	/* as per PICIO (1.0==2048) */
#define	    PIXARLOGDATAFMT_16BIT	4	/* signed short samples */
#define	    PIXARLOGDATAFMT_FLOAT	5	/* IEEE float samples */
# 65550-65556 are allocated to Oceana Matrix <dev@oceana.com> */
	DCSIMAGERTYPE = 65550 			 #  imager model & filter
#define     DCSIMAGERMODEL_M3           0       /* M3 chip (1280 x 1024) */
#define     DCSIMAGERMODEL_M5           1       /* M5 chip (1536 x 1024) */
#define     DCSIMAGERMODEL_M6           2       /* M6 chip (3072 x 2048) */
#define     DCSIMAGERFILTER_IR          0       /* infrared filter */
#define     DCSIMAGERFILTER_MONO        1       /* monochrome filter */
#define     DCSIMAGERFILTER_CFA         2       /* color filter array */
#define     DCSIMAGERFILTER_OTHER       3       /* other filter */
	DCSINTERPMODE = 65551 			 #  interpolation mode
#define     DCSINTERPMODE_NORMAL        0x0     /* whole image, default */
#define     DCSINTERPMODE_PREVIEW       0x1     /* preview of image (384x256) */
	DCSBALANCEARRAY = 65552 			 #  color balance values
	DCSCORRECTMATRIX = 65553 			 #  color correction values
	DCSGAMMA = 65554 			 #  gamma value
	DCSTOESHOULDERPTS = 65555 			 #  toe & shoulder points
	DCSCALIBRATIONFD = 65556 			 #  calibration file desc
# Note: quality level is on the ZLIB 1-9 scale. Default value is -1
	ZIPQUALITY = 65557 			 #  compression quality level
	PIXARLOGQUALITY = 65558 			 #  PixarLog uses same scale
# 65559 is allocated to Oceana Matrix <dev@oceana.com>
	DCSCLIPRECTANGLE = 65559 			 #  area of image to acquire
	SGILOGDATAFMT = 65560 			 #  SGILog user data format
#define     SGILOGDATAFMT_FLOAT		0	/* IEEE float samples */
#define     SGILOGDATAFMT_16BIT		1	/* 16-bit samples */
#define     SGILOGDATAFMT_RAW		2	/* uninterpreted data */
#define     SGILOGDATAFMT_8BIT		3	/* 8-bit RGB monitor values */
	SGILOGENCODE = 65561 			 #  SGILog data encoding control
#define     SGILOGENCODE_NODITHER	0     /* do not dither encoded values*/
#define     SGILOGENCODE_RANDITHER	1     /* randomly dither encd values */
	LZMAPRESET = 65562 			 #  LZMA2 preset (compression level)
	PERSAMPLE = 65563 			 #  interface for per sample tags
#define     PERSAMPLE_MERGED        0	/* present as a single value */
#define     PERSAMPLE_MULTI         1	/* present as multiple values */

# EXIF tags

	EXIFTAG_EXPOSURETIME = 33434 			 #  Exposure time
	EXIFTAG_FNUMBER = 33437 			 #  F number
	EXIFTAG_EXPOSUREPROGRAM = 34850 			 #  Exposure program
	EXIFTAG_SPECTRALSENSITIVITY = 34852 			 #  Spectral sensitivity
	EXIFTAG_ISOSPEEDRATINGS = 34855 			 #  ISO speed rating
	EXIFTAG_OECF = 34856 			 #  Optoelectric conversion factor
	EXIFTAG_EXIFVERSION = 36864 			 #  Exif version
	EXIFTAG_DATETIMEORIGINAL = 36867 			 #  Date and time of original data generation
	EXIFTAG_DATETIMEDIGITIZED = 36868 			 #  Date and time of digital data generation
	EXIFTAG_COMPONENTSCONFIGURATION = 37121 			 #  Meaning of each component
	EXIFTAG_COMPRESSEDBITSPERPIXEL = 37122 			 #  Image compression mode
	EXIFTAG_SHUTTERSPEEDVALUE = 37377 			 #  Shutter speed
	EXIFTAG_APERTUREVALUE = 37378 			 #  Aperture
	EXIFTAG_BRIGHTNESSVALUE = 37379 			 #  Brightness
	EXIFTAG_EXPOSUREBIASVALUE = 37380 			 #  Exposure bias
	EXIFTAG_MAXAPERTUREVALUE = 37381 			 #  Maximum lens aperture
	EXIFTAG_SUBJECTDISTANCE = 37382 			 #  Subject distance
	EXIFTAG_METERINGMODE = 37383 			 #  Metering mode
	EXIFTAG_LIGHTSOURCE = 37384 			 #  Light source
	EXIFTAG_FLASH = 37385 			 #  Flash
	EXIFTAG_FOCALLENGTH = 37386 			 #  Lens focal length
	EXIFTAG_SUBJECTAREA = 37396 			 #  Subject area
	EXIFTAG_MAKERNOTE = 37500 			 #  Manufacturer notes
	EXIFTAG_USERCOMMENT = 37510 			 #  User comments
	EXIFTAG_SUBSECTIME = 37520 			 #  DateTime subseconds
	EXIFTAG_SUBSECTIMEORIGINAL = 37521 			 #  DateTimeOriginal subseconds
	EXIFTAG_SUBSECTIMEDIGITIZED = 37522 			 #  DateTimeDigitized subseconds
	EXIFTAG_FLASHPIXVERSION = 40960 			 #  Supported Flashpix version
	EXIFTAG_COLORSPACE = 40961 			 #  Color space information
	EXIFTAG_PIXELXDIMENSION = 40962 			 #  Valid image width
	EXIFTAG_PIXELYDIMENSION = 40963 			 #  Valid image height
	EXIFTAG_RELATEDSOUNDFILE = 40964 			 #  Related audio file
	EXIFTAG_FLASHENERGY = 41483 			 #  Flash energy
	EXIFTAG_SPATIALFREQUENCYRESPONSE = 41484 			 #  Spatial frequency response
	EXIFTAG_FOCALPLANEXRESOLUTION = 41486 			 #  Focal plane X resolution
	EXIFTAG_FOCALPLANEYRESOLUTION = 41487 			 #  Focal plane Y resolution
	EXIFTAG_FOCALPLANERESOLUTIONUNIT = 41488 			 #  Focal plane resolution unit
	EXIFTAG_SUBJECTLOCATION = 41492 			 #  Subject location
	EXIFTAG_EXPOSUREINDEX = 41493 			 #  Exposure index
	EXIFTAG_SENSINGMETHOD = 41495 			 #  Sensing method
	EXIFTAG_FILESOURCE = 41728 			 #  File source
	EXIFTAG_SCENETYPE = 41729 			 #  Scene type
	EXIFTAG_CFAPATTERN = 41730 			 #  CFA pattern
	EXIFTAG_CUSTOMRENDERED = 41985 			 #  Custom image processing
	EXIFTAG_EXPOSUREMODE = 41986 			 #  Exposure mode
	EXIFTAG_WHITEBALANCE = 41987 			 #  White balance
	EXIFTAG_DIGITALZOOMRATIO = 41988 			 #  Digital zoom ratio
	EXIFTAG_FOCALLENGTHIN35MMFILM = 41989 			 #  Focal length in 35 mm film
	EXIFTAG_SCENECAPTURETYPE = 41990 			 #  Scene capture type
	EXIFTAG_GAINCONTROL = 41991 			 #  Gain control
	EXIFTAG_CONTRAST = 41992 			 #  Contrast
	EXIFTAG_SATURATION = 41993 			 #  Saturation
	EXIFTAG_SHARPNESS = 41994 			 #  Sharpness
	EXIFTAG_DEVICESETTINGDESCRIPTION = 41995 			 #  Device settings description
	EXIFTAG_SUBJECTDISTANCERANGE = 41996 			 #  Subject distance range
	EXIFTAG_IMAGEUNIQUEID = 42016 			 #  Unique image ID
end

@enum SampleFormats begin
	SAMPLEFORMAT_UINT = 1				# !unsigned integer data
	SAMPLEFORMAT_INT = 2				# !signed integer data
	SAMPLEFORMAT_IEEEFP = 3				# !IEEE floating point data
	SAMPLEFORMAT_VOID = 4				# !untyped data
	SAMPLEFORMAT_COMPLEXINT = 5			# !complex signed int
	SAMPLEFORMAT_COMPLEXIEEEFP = 6		# !complex ieee floating
end

@enum PhotometricInterpretations begin
	PHOTOMETRIC_MINISWHITE = 0		# min value is white
	PHOTOMETRIC_MINISBLACK = 1		# min value is black
	PHOTOMETRIC_RGB = 2				# RGB color model
	PHOTOMETRIC_PALETTE = 3			# color map indexed
	PHOTOMETRIC_MASK = 4			# $holdout mask
	PHOTOMETRIC_SEPARATED = 5		# !color separations
	PHOTOMETRIC_YCBCR = 6			# !CCIR 601
	PHOTOMETRIC_CIELAB = 8			# !1976 CIE L*a*b*
	PHOTOMETRIC_ICCLAB = 9			# ICC L*a*b* [Adobe TIFF Technote 4]
	PHOTOMETRIC_ITULAB = 10			# ITU L*a*b*
	PHOTOMETRIC_CFA = 32803			# color filter array
	PHOTOMETRIC_LOGL = 32844		# CIE Log2(L)
	PHOTOMETRIC_LOGLUV = 32845		# CIE Log2(L) (u',v')
end

@enum CompressionType begin
	COMPRESSION_NONE = 1			# No compression
	COMPRESSION_CCITTRLE = 2		# CCITT modified Huffman RLE
	COMPRESSION_CCITT_T4 = 3        # CCITT T.4 (TIFF 6 name)
	COMPRESSION_CCITT_T6 = 4 		# CCITT T.6 (TIFF 6 name)
	COMPRESSION_LZW = 5				# Lempel-Ziv  & Welch
	COMPRESSION_OJPEG = 6			# !6.0 JPEG
	COMPRESSION_JPEG = 7			# %JPEG DCT compression
	COMPRESSION_T85	= 9				# !TIFF/FX T.85 JBIG compression
	COMPRESSION_T43 = 10			# !TIFF/FX T.43 colour by layered JBIG compression
	COMPRESSION_NEXT = 32766		# NeXT 2-bit RLE
	COMPRESSION_CCITTRLEW = 32771	#1 w/ word alignment
	COMPRESSION_PACKBITS = 32773	# Macintosh RLE
	COMPRESSION_THUNDERSCAN	= 32809	# ThunderScan RLE
# codes 32895-32898 are reserved for ANSI IT8 TIFF/IT <dkelly@apago.com)
	COMPRESSION_IT8CTPAD = 32895 	# IT8 CT w/padding
	COMPRESSION_IT8LW = 32896   	# IT8 Linework RLE
	COMPRESSION_IT8MP = 32897   	# IT8 Monochrome picture
	COMPRESSION_IT8BL = 32898   	# IT8 Binary line art
# compression codes 32908-32911 are reserved for Pixar
	COMPRESSION_PIXARFILM = 32908   # Pixar companded 10bit LZW
	COMPRESSION_PIXARLOG = 32909   	# Pixar companded 11bit ZIP
	COMPRESSION_DEFLATE = 32946		# Deflate compression
	COMPRESSION_ADOBE_DEFLATE = 8	# Deflate compression, as recognized by Adobe
# compression code 32947 is reserved for Oceana Matrix <dev@oceana.com>
	COMPRESSION_JBIG = 34661		# ISO JBIG
	COMPRESSION_SGILOG = 34676		# SGI Log Luminance RLE
	COMPRESSION_SGILOG24 = 34677	# SGI Log 24-bit packed
	COMPRESSION_JP2000 = 34712   	# Leadtools JPEG2000
	COMPRESSION_LZMA = 34925		# LZMA2
end

@enum ExtraSamples begin
	EXTRASAMPLE_UNSPECIFIED = 0		# unspecified data
	EXTRASAMPLE_ASSOCALPHA = 1 		# associated alpha data
	EXTRASAMPLE_UNASSALPHA = 2 		# unassociated alpha data
end


text:
using Rotations
using StaticArrays
using Printf

function readBundleFile(bundle_file)
    fin = open(bundle_file)

    # read the header
    header = readline(fin)
    nCams, nPts = parse.(Int, split(readline(fin)))
    println("found $nCams cameras and $nPts points")

    # read the camera blocks
    camera_params = []
    for i=1:nCams
        f,   k1,  k2 = parse.(Float64, split(readline(fin)))
        r11, r12, r13 = parse.(Float64, split(readline(fin)))
        r21, r22, r23 = parse.(Float64, split(readline(fin)))
        r31, r32, r33 = parse.(Float64, split(readline(fin)))
        t1,  t2,  t3 = parse.(Float64, split(readline(fin)))
        R = @SMatrix [r11 r12 r13; r21 r22 r23; r31 r32 r33]
        aa = AngleAxis(RotMatrix{3}(R))
        v = rotation_angle(aa) .* rotation_axis(aa)
        push!(camera_params, [v[1], v[2], v[3], t1, t2, t3, f, k1, k2])
    end

    # read the point blocks
    observations = []
    point_params = []
    for i=1:nPts
        X1, X2, X3 = parse.(Float64, split(readline(fin))) # position
        readline(fin) # color (ignore)
        tokens = split(readline(fin))

        num_observations = parse(Int, tokens[1])
        for k=1:num_observations
            cam_id = parse(Int, tokens[(k-1)*4 + 2])
            key_id = parse(Int, tokens[(k-1)*4 + 3])
            xcoord = parse(Float64, tokens[(k-1)*4 + 4])
            ycoord = parse(Float64, tokens[(k-1)*4 + 5])

            push!(observations, [cam_id, i-1, xcoord, ycoord])
        end
        # Note to self: 0,0 is the center of the image.
        # (-w/2, -h/2) is the lower-left corner.
        push!(point_params, [X1, X2, X3])
    end

    close(fin)
    return camera_params, point_params, observations
end

function writeBALFile(bal_file, camera_params, point_params, observations)
    fout = open(bal_file, "w")
    @printf(fout, "%d %d %d\n", length(camera_params), length(point_params), length(observations))
    for obs in observations
        @printf(fout, "%d %d %lf %lf\n", obs...)
    end
    for cam in camera_params
        @printf(fout, "%lf %lf %lf %lf %lf %lf %lf %lf %lf\n", cam...)
    end
    for pt in point_params
        @printf(fout, "%lf %lf %lf\n", pt...)
    end
    close(fout)
end

function reindexProblem(camera_params, point_params, observations)

    changed_something = true
    while changed_something
        println("Start of pruning round.")
        changed_something = false

        # look for cameras that see fewer than 5 points
        # count how many obs per camera
        camera_counts = zeros(length(camera_params))
        for obs in observations
            cam_index = convert(Int, obs[1])
            camera_counts[cam_index+1] += 1
        end
        mask = camera_counts .< 5
        if any(mask)
            changed_something = true

            # compute the reindexing scheme
            old2new = Dict()
            count = 0
            for i=1:length(camera_params)
                if camera_counts[i] >= 5
                    old2new[i-1] = count
                    count += 1
                end
            end

            # apply the reindexing:
            new_observations = []
            for obs in observations
                cam_index = convert(Int, obs[1])
                if camera_counts[cam_index+1] >= 5
                    push!(new_observations, [old2new[cam_index], obs[2], obs[3], obs[4]])
                end
            end
            println("got rid of $(length(observations) - length(new_observations)) cameras.")
            observations = new_observations
            camera_params = camera_params[.!mask]
        end

        # look for points seen by fewer than two cameras
        # count how many obs per point
        point_counts = zeros(length(point_params))
        for obs in observations
            pt_index = convert(Int, obs[2])
            point_counts[pt_index+1] += 1
        end
        mask = point_counts .< 2
        if any(mask)
            changed_something = true

            # compute the reindexing scheme
            old2new = Dict()
            count = 0
            for i=1:length(point_params)
                if point_counts[i] >= 2
                    old2new[i-1] = count
                    count += 1
                end
            end

            # apply the reindexing:
            new_observations = []
            for obs in observations
                pt_index = convert(Int, obs[2])
                if point_counts[pt_index+1] >= 2
                    push!(new_observations, [obs[1], old2new[pt_index], obs[3], obs[4]])
                end
            end
            println("got rid of $(length(observations) - length(new_observations)) points.")
            observations = new_observations
            point_params = point_params[.!mask]
        end

    end
    println("Done pruning.")
    return camera_params, point_params, observations
end

bundle_file = ARGS[1]
println("reading in bundle file: $bundle_file")
bal_file = ARGS[2]
println("writing to bal file: $bal_file")
camera_params, point_params, observations = readBundleFile(bundle_file)
camera_params, point_params, observations = reindexProblem(camera_params, point_params, observations)
writeBALFile(bal_file, camera_params, point_params, observations)

max_stars_repo_name: wilsonkl/sfmflex-release

max_issues_repo_name: wilsonkl/sfmflex-release

max_forks_repo_name: wilsonkl/sfmflex-release


text:
using Test
using MDCDL
using FFTW
using ImageFiltering
using Random
using LinearAlgebra

@testset "Multidimensionals" begin
    Random.seed!(39273529)

    @testset "representation of linear operator" begin
        @testset "FFT" begin
            for d = 1:5
                for n = 1:10
                    szx = rand(1:8,d)
                    x = rand(szx...)

                    A = MDCDL.representationmatrix(fft, szx...)
                    @test A*vec(x) ≈ vec(fft(x))
                end
            end
        end

        @testset "DCT" begin
            for d = 1:5
                for n = 1:10
                    szx = rand(1:8,d)
                    x = rand(szx...)

                    A = MDCDL.representationmatrix(dct, szx...)
                    @test A*vec(x) ≈ vec(dct(x))
                end
            end
        end

        @testset "imfilter" begin
            for d = 1:3
                for n = 1:10
                    szx = rand(1:8,d)
                    x = rand(szx...)
                    szker = rand.(UnitRange.(1, szx))
                    ker = centered(rand(szker...))

                    myfcn(t) = imfilter(t, ker, Pad(:circular))

                    A = MDCDL.representationmatrix(myfcn, szx...)
                    @test A*vec(x) ≈ vec(myfcn(x))
                end
            end
        end
    end

    @testset "permutated DCT" begin
        @testset "unitarity" begin
            for d = 1:5
                for n = 1:10
                    szx = rand(1:8, d)
                    x = rand(szx...)

                    A = MDCDL.permdctmtx(szx...)

                    @test A' * A ≈ I
                end
            end
        end
        @testset "symmetry" begin
            for d = 1:5
                for n = 1:10
                    szx = rand(1:8, d)
                    x = rand(szx...)

                    A = MDCDL.permdctmtx(szx...)

                    nch = size(A, 1)
                    Γ = diagm(0 => [ ones(cld(nch, 2)); -ones(fld(nch, 2)) ])
                    @test Γ * reverse(A, dims=2) ≈ A
                end
            end
        end
        # for d = 1:5
        #     for n = 1:10
        #         szx = rand(1:8,d)
        #         x = rand(szx...)
        #
        #         A = MDCDL.permdctmtx(szx...)
        #         @test sort(A*vec(x)) ≈ sort(vec(dct(x)/sqrt(prod(szx))))
        #     end
        # end
    end

    @testset "centered DFT" begin
        @testset "unitarity" begin
            for d = 1:5
                for n = 1:10
                    szx = rand(1:8, d)
                    x = rand(szx...)

                    A = MDCDL.cdftmtx(szx...)

                    @test A' * A ≈ I
                end
            end
        end
        @testset "symmetry" begin
            for d = 1:5
                for n = 1:10
                    szx = rand(1:8, d)
                    x = rand(szx...)

                    A = MDCDL.cdftmtx(szx...)

                    @test conj(reverse(A, dims=2)) ≈ A
                end
            end
        end
        # for d = 1:5
        #     for n = 1:10
        #         szx = rand(1:8,d)
        #         x = rand(szx...)
        #
        #         A = MDCDL.permdctmtx(szx...)
        #         @test sort(A*vec(x)) ≈ sort(vec(dct(x)/sqrt(prod(szx))))
        #     end
        # end
    end

    @testset "polyphase matrix" begin
        for d = 1:5
            for n = 1:10
                df = (rand(1:4, d)...,)
                szx = ((rand(1:3, d) .* df)...,)
                x = rand(szx...)

                pvx = MDCDL.mdarray2polyphase(x, df)

                @test size(pvx.data, 1) == prod(df)
                @test size(pvx.data, 2) == fld(length(x), prod(df))

                rx = MDCDL.polyphase2mdarray(pvx, df)
                @test rx == x
            end
        end
    end
end

nothing

max_stars_repo_name: nixir/MDCDL

max_issues_repo_name: nixir/MDCDL

max_forks_repo_name: nixir/MDCDL


text:
#= cd(@__DIR__)
using Pkg
Pkg.activate(".")

using Roots
using Plots
using Latexify
using Unitful
using UnitfulAtomic
using UnitfulEquivalences
using UnitfulRecipes
pgfplotsx()

kC =  -14.0*ustrip(auconvert(1.60217662e-19u"C"/8.85418782e-12u"F/m"/4.0/π/11.7))

m = ustrip(auconvert(9.109e-31u"kg"))
V₀ = ustrip(auconvert(-1e6u"eV"))
a = ustrip(auconvert(5.4e-10u"m"))
b0 = ustrip(auconvert(3.9e-15u"m"))
ħ = ustrip(auconvert(6.626070715e-34/2.0/π * u"J*s"))
E = ustrip(auconvert.(collect(LinRange(-200.0,200.0,1000000)).*u"eV")).*(1+0.0im)
q(E) = sqrt(2.0*m*E/(ħ^2))

N = 5
Ve(x,p) = kC * 1/abs(x-p)
r = -a*2.0:b0*50.0:a*2.0
Pot = similar(r).*0.0
for i in 1:N
    Pot += Ve.(r,((i-1)-(N-1)/2)*a)
end
plot(r,Pot,lw=2.0,dpi=300,legend=false)
ylims!(-25,1.0)
xlims!(-1*a,1*a)
rect(x,x0,b,V0) = V0*(((x-x0) > -b/2) & ((x-x0) < b/2))
Ry = .-rect.(r,-a,10000*b0,25.0).-rect.(r,0.0,10000*b0,25.0).-rect.(r,a,10000*b0,25.0)
plot!(r,Ry ,color=:black,legend=false)
k₂(E) = sqrt(2.0*m*(V₀-E)/(ħ^2))

F(E)=cos(q(E)*a)*cosh(k₂(E)*b0) + (k₂(E)^2-q(E)^2)/(2*k₂(E)*q(E))*sin(q(E)*a)*sinh(k₂(E)*b0)
y = F.(E)
plot(real.(y))
gui()
yvalid = copy(y)
yvalid[ab0s.(real.(y)) .> 1.0] .= NaN64

Cf = ustrip(auconvert(1.0u"eV"))
x = real.(acos.(yvalid))./pi
y = real.(E)./Cf

x = vcat(-x,x,-x,x)
y = vcat(y,y,-y,-y)

plot(x,y)
ylims!(-5,5) =#


text:
using Distributions
using StatsBase
using StatsFuns
using FastGaussQuadrature
using KernelFunctions
using UnPack
using LinearAlgebra
using PDMats
using Plots, StatsPlots
using FiniteDiff
using Random123
using RandomNumbers
using Roots
using Optim

struct PBODataPoint{F<:Real}
    α::F
    βs::AbstractVector{F}
    ξ::AbstractVector{F}
    x::AbstractVector{F}
end

struct PBOHyperparameters{F<:Real}
    ℓσ::F
    ℓl::AbstractVector{F}
    ℓσn::F
end

function pgp_ℓpl(Δ, hypers::PBOHyperparameters)
    n_itgr = 32
    itgr_x, itgr_w = gausshermite(n_itgr)
    σ = exp(hypers.ℓσn)
    m = size(Δ, 2)

    # convolution of a gaussian pdf and cdf
    ϕ★ϕ(x) = dot(itgr_w, normcdf.(x .- itgr_x*√2)) / √π

    -sum(ϕ★ϕ.(Δ)) / m
end

function pgp_ℓpθ(hypers::PBOHyperparameters)
    @unpack ℓσ, ℓl, ℓσn = hypers
    pℓσ  = logpdf(Normal(0.0, 2.0),          ℓσ)
    pℓl  = logpdf(MvNormal(length(ℓl), 2.0), ℓl)
    pℓσn = logpdf(Normal(0.0, 2.0),         ℓσn)
    pℓσ + pℓl + pℓσn
end

function pgp_kernel(hypers::PBOHyperparameters)
    @unpack ℓσ, ℓl, ℓσn = hypers
    k     = Matern52Kernel()
    l⁻¹   = exp.(-ℓl)
    t_ard = ARDTransform(l⁻¹)
    σ²    = exp(2*ℓσ)
    σ² * transform(k, t_ard)
end

function pgp_gram_matrix(data::Array{<:PBODataPoint},
                         hypers::PBOHyperparameters)
    n_dims = length(data[1].x)
    n_data = length(data)
    n_β    = length(data[1].βs)
    n_f    = n_data * (1 + n_β)
    X      = Array{Float64}(undef, n_dims, n_f)

    for i = 1:n_data
        @unpack α, βs, ξ, x = data[i]
        f_idx      = (i - 1)*(1 + n_β) + 1
        X[:,f_idx] = x + ξ*α
        for j = 1:n_β
            X[:,f_idx+j] = x + ξ*βs[j]
        end
    end

    k      = pgp_kernel(hypers)
    Σ      = kernelmatrix(k, X, obsdim=2)
    display(heatmap(Σ))
    PDMats.PDMat(Σ), X
end

# function pgp_pmℓl(latents::AbstractVector{<:PBODataPoint},
#                   data::AbstractVector{<:PBODataPoint},
#                   hypers::PBOHyperparameters)
#     Σ   = PDMats.PDMat(Σ)
#     ℓpf = -1/2*invquad(Σ, f)

#     ℓpθ = pgp_ℓpθ(hypers)
#     ℓq  = ℓpf + ℓpy + ℓpθ

#     Σ, ℓq
# end

function pgp_is_α(idx, Δ_shape)
    n_α = Δ_shape[1]
    n_β = Δ_shape[2]
    mod(idx, n_β+1) == 1
end

function pgp_α_idx(idx, Δ_shape)
    n_β = Δ_shape[2]
    div(idx - 1, n_β+1) + 1
end

function pgp_β_idx(idx, Δ_shape)
    n_β = Δ_shape[2]
    mod(idx - 1, n_β+1)
end

function pgp_∇T_Λ(Σ, f, Δ, hypers::PBOHyperparameters)
    σ    = exp(hypers.ℓσ)
    ϕ₂   = Normal(0, √2)
    m    = size(Δ, 2)
    n_α  = size(Δ, 1)
    n_β  = m
    n_f  = prod(size(Δ)) + size(Δ, 1)
    ϕ₂Δ  = pdf(ϕ₂, Δ)

    # Gradient
    ∇q  = zeros(n_f)
    idx = 1
    for idx = 1:n_f
        if(pgp_is_α(idx, size(Δ)))
            i = pgp_α_idx(idx, size(Δ))
            ∇q[idx] = sum(view(ϕ₂Δ,i,:)) / σ / m
        else
            i = pgp_α_idx(idx, size(Δ))
            j = pgp_β_idx(idx, size(Δ))
            ∇q[idx] = -ϕ₂Δ[i,j] / σ / m
        end
        idx += 1
    end
    ∇T = -(Σ \ f) + ∇q

    # Hessian
    Λ = zeros(n_f, n_f)
    for idx = 1:(1+n_β):n_f
        # fᵢ ∈ α, fⱼ ∈ α
        i = pgp_α_idx(idx, size(Δ))
        Λ[idx,idx] = dot(view(ϕ₂Δ,i,:), view(Δ,i,:)) / 2 / σ^2 / m

        # fᵢ ∈ α, fⱼ ∈ β
        for j = 1:n_β
            β_idx = idx + j
            Λ[idx, β_idx] = -Δ[i,j] * ϕ₂Δ[i,j] / 2 / σ^2 / m
        end

        # fᵢ ∈ β, fⱼ ∈ β
        for j = 1:n_β
            β_idx = idx + j
            Λ[β_idx, β_idx] = Δ[i,j] * ϕ₂Δ[i,j] / 2 / σ^2 / m
        end
    end
    Λ = Symmetric(Λ)
    ∇T, Λ
end

function pgp_Δ(f, data, hypers)
    σ   = exp(hypers.ℓσ)
    n_α = length(data)
    n_β = length(data[1].βs)
    Δ   = zeros(n_α, n_β)
    for i = 1:n_α
        for j = 1:n_β
            α_idx = (i-1) * n_β + 1
            β_idx = (i-1) * n_β + j + 1
            Δ[i,j] = (f[β_idx] - f[α_idx]) / σ
        end
    end
    Δ
end

function derivative_test()
    ℓθ  = PBOHyperparameters(0.0, zeros(2), 0.0)
    k   = pgp_kernel(ℓθ)
    n_f = (1+5)*4
    X   = rand(2, n_f)
    Σ   = PDMats.PDMat(kernelmatrix(k, X, obsdim=2))

    σ = exp(ℓθ.ℓσ)
    f = randn(n_f)*10
    Δ = zeros(4, 5)
    for i = 1:4
        for j = 1:5
            α_idx = (i-1) * 6 + 1
            β_idx = (i-1) * 6 + j + 1
            Δ[i,j] = (f[β_idx] - f[α_idx]) / σ
        end
    end

    ∇T, Λ = pgp_∇T_Λ(Σ, f, Δ, ℓθ)
    ∇²T   = -inv(Σ) + Λ

    T(f_in) = begin
        σ    = exp(ℓθ.ℓσ)
        Δ_in = zeros(4, 5)
        for i = 1:4
            for j = 1:5
                α_idx = (i-1) * 6 + 1
                β_idx = (i-1) * 6 + j + 1
                Δ_in[i,j] = (f_in[β_idx] - f_in[α_idx]) / σ
            end
        end
        -0.5*invquad(Σ, f_in) + pgp_ℓpl(Δ_in, ℓθ)
    end
    ∇T̂  = FiniteDiff.finite_difference_gradient(T, f)
    ∇²T̂ = FiniteDiff.finite_difference_hessian(T, f)

    println("∇²T norm: ", norm(∇²T - ∇²T̂))
    println("∇T norm:  ", norm(∇T  - ∇T̂))
end

function pgp_laplace(Σ, data, f_init, hypers)
#=
    Variant of the Newton's method based mode-locating algorithm (GPML, Algorithm 3.1)
    Utilizes the Woodburry identity for avoiding two cholesky factorizations
    per Newton iteration.
    Reduces the stepsize whenever the marginal likelhood gets stuck
    Algortihm 3.1 utilizes the fact that W is diagonal which is not for our case.
    Note: ( K^{-1} + W )^{-1} = K ( I - ( I + W K )^{-1} ) W K
=##
    α        = 0.5
    max_iter = 50
    f        = f_init
    Δ        = pgp_Δ(f, data, hypers)
    ∇T, Λ    = pgp_∇T_Λ(Σ, f, Δ, hypers)
    for i = 1:max_iter
        W     = -Λ
        WK    = W*Σ
        b     = W*f + ∇T
        B     = I + WK
        Blu   = lu(B)
        a     = (b - Blu \ (WK*b))
        p     = (Σ*a) - f

        fₖ     = deepcopy(f)
        α      = 1.0
        c₂     = 0.9
        c₂pᵀ∇T = c₂*dot(p, ∇T)
        while(true)
            f     = fₖ + α*p
            Δ     = pgp_Δ(f, data, hypers)
            ∇T, Λ = pgp_∇T_Λ(Σ, f, Δ, hypers)

            if(-dot(p, ∇T) ≤ -c₂pᵀ∇T || α < 1e-4)
                break
            end
            α /= 2.0
        end
        if(α < 1e-4)
            break
        end

        T = -0.5*invquad(Σ, f) + pgp_ℓpl(Δ, hypers)
        @info "Newton step" iter=i norm∇T=norm(∇T) T=T α=α
    end
    f
end

function ackley(x::Vector{<:Real}, a = 20, b = 0.2, c = 2*π)
    d = length(x);
    term1 = -a * exp(-b * √(sum(x.^2) / d))
    term2 = -exp(sum(cos.(c*x)) / d)
    return term1 + term2 + a + exp(1)
end

function pgp_find_bounds(x, ξ)
    # Find the upper/lower bound of α
    # where x + αξ is within the hyper-cube [0, 1]ᴰ
    # by solving the system of equations:
    # xᵢ + αξᵢ = 1, if α ≥ 0, α ∈ A+, else then α ∈ A-
    # xᵢ + αξᵢ = 0, if α ≥ 0, α ∈ A+, else then α ∈ A-
    # α upper bound = max{α | α ∈ A+}
    # α lower bound = min{α | α ∈ A-}

    lb   = -Inf
    ub   = Inf
    dims = length(x)

    for i = 1:dims
        α = (1 - x[i]) / ξ[i]
        if(α > 0)
            ub = min(α, ub)
        else
            lb = max(α, lb)
        end

        α = -x[i] / ξ[i]
        if(α > 0)
            ub = min(α, ub)
        else
            lb = max(α, lb)
        end
    end
    lb, ub
end

function laplace_test()
    dims     = 2
    n_points = 10
    f        = x->ackley((x .- 0.5) * 20)
    prng     = Philox4x(UInt64, (0x22cf80fe510ae4bd, 0x1b23a6a16922848f,), 10)
    Random123.set_counter!(prng, 1)

    data = PBODataPoint{Real}[]
    for i = 1:n_points
        ξ  = randn(prng, dims)
        ξ /= maximum(ξ)
        x  = rand(prng, dims)

        lb, ub = pgp_find_bounds(x, ξ)

        res = optimize(α -> f(ξ*α + x), lb, ub)
        α   = Optim.minimizer(res)[1]
        β   = rand(prng, Uniform(lb, ub), 5)

        datum = PBODataPoint{Real}(α, β, ξ, x)
        push!(data, datum)
    end

    hypers = PBOHyperparameters(-1.0, fill(-1.5, 2), -1.0)
    Σ, X   = pgp_gram_matrix(data, hypers)
    f      = pgp_laplace(Σ, data, zeros(size(Σ, 1)), hypers)
    display(scatter(X[1,:], X[2,:], f))
end


text:
"""
`f` from page `a.jl`.

Links:

- [`ccall`](@ref)
- [`while`](@ref)
- [`@time(x)`](@ref)
- [`T(x)`](@ref)
- [`T(x, y)`](@ref)
- [`f(::Integer)`](@ref)
- [`f(::Any)`](@ref)
- [`f(::Any, ::Any)`](@ref)
- [`f(x, y, z)`](@ref)

[^footnote]:

    Footnote contents. [^footnote]

"""
f(x) = x