text:
abstract MonteCarloScenarios
immutable FixedMonteCarloScenarios <: MonteCarloScenarios
  scenarios::Reservoir
end
immutable GeneratedMonteCarloScenarios <: MonteCarloScenarios
  method::Symbol
end

abstract ReportEvaluator
immutable CrossValidatedFeasibilityProbability <: ReportEvaluator end
immutable InverseModel <: ReportEvaluator end
immutable MonteCarloFeasibilityProbability <: ReportEvaluator
  scenarios::MonteCarloScenarios
end

"""
Prints a complete solution report for a reservoir, with multiple solutions and their evaluation.

Mandatory arguments:

  * `reservoir`: the reservoir to study
  * `ruleDuration`: the duration of the rule curves to produce
  * `guarantee`: the duration to guarantee for those rule curves

Optional arguments for input data:

  * `currentRuleCurve`: the existing rule curve, used as a benchmark for the newly computed ones.
    Its periodicity must be that of `reservoir`.
  * `climateChangeScenarios`: a set of inflow scenarios for some climate change, associating a description
    to the inflow. There is no distinction between the rivers in this argument!

Optional arguments to select the solvers and evaluators to use:

  * `scenarioGenerationAlgorithms`: algorithms used to generate scenarios (see Scenario generation module, the symbols
    are the argument to `scenarioGeneration`).
    Default: `[:Merge, :Mix]`
  * `ciAlgorithmsLevels`: algorithms used to compute the confidence intervals, associated with their confidence values
    (see Statistical module, the symbols are the argument to `generateScenarios`).
    Default: `Dict(:TStudent => [.95])`
  * `solvers`: solvers that are used to generate solutions (see Solver module, the symbols correspond to the function
    names of the solvers).
    Default: `[:minimumRuleCurve, :safeMinimumRuleCurve]`
  * `evaluators`: algorithms that are used to evaluate the new rule curves (see Solver module, the symbols correspond
    to the function names of the solvers).
    Default: `[]`
    TODO: not really flexible (think of parameters).

    The current list of possible evaluators is:

    * `:feasibilityProbabilitySyntheticScenarios`: computes the feasibility probability on synthetic scenarios
      (generated based on the statistical analysis of the data sets)
    * `:feasibilityProbabilityAllDatasets`: computes the feasibility probability on historical scenarios (all the given
       data sets), i.e. resubstitution
    * `:feasibilityProbabilityClimateChange`: computes the feasibility probability on provided climate change scenarios
      (where all rivers are merged)
    * `:hydropower`: evaluates the hydropower potential (with a generated price scenario)
      TODO: price scenario from outside.
    * `:worstCaseFunctions`:
    * `:maximumDeviation`:

  * `globalEvaluators`: algorithms that are used to evaluate algorithms that generate the rule curves, and not the
    rule curves themselves (see Solver module, the symbols correspond to the function names of the solvers).
    Default: `[]`
    TODO: not really flexible (think of parameters).
  * `disabledCombinations`: combinations of the above parameters that should *not* be run. This list contains pairs
    of a scenario generation algorithm and a solver. It is not possible to disable confidence intervals this way.
    Default: `[]`.
    Recommended value: `[(:Mix, :safeMinimumRuleCurve)]`, as this combination takes a lot of time to compute.
  * `verbosity`: passed as is to the solvers (`-1` to disable all output)
  * `saveFolder`: folder name to use with the following options
  * `savePlots`: whether the script should prepare some plots
  * `saveIntermediate`: whether the script should save the outputs for future reuse

TODO: how to add new solvers/evaluators? Work with a hash table and store lambdas with a fixed signature?
TODO: how to deal with solvers/evaluators parameters? Still via lambdas? Within the hash table?
TODO: provide more types to avoid typing Tuple{Tuple} too often. Use the defined data structures.
"""
function report(reservoir::Reservoir;
                ruleDuration::Period=Year(1), guarantee::Period=Year(2), ciPeriod::Period=Week(1),
                currentRuleCurve::Array{Float64, 1}=zeros(Float64, 0),
                climateChangeScenarios::Dict{Symbol, Array{TimeSeries.TimeArray{Float64,1,DateTime,Array{Float64,1}}, 1}}=Dict{Symbol, Array{TimeSeries.TimeArray{Float64,1,DateTime,Array{Float64,1}}, 1}}(),

                scenarioGenerationAlgorithms::Array{Symbol, 1}=Symbol[:Merge, :Mix],
                ciAlgorithmsLevels::Dict{Symbol, Array{Float64, 1}}=Dict{Symbol, Array{Float64, 1}}(:TStudent => [.95, .99]),
                solvers::Array{Symbol, 1}=[:minimumRuleCurve, :safeMinimumRuleCurve],
                evaluators::Array{Symbol, 1}=Symbol[],#:feasibilityProbabilityAllDatasets,:feasibilityProbabilityClimateChange], # :feasibilityProbability, # , :hydropower # , :worstCaseFunctions, :maximumDeviation
                globalEvaluators::Array{Symbol, 1}=Symbol[],#:crossValidatedFeasibilityProbability],
                disabledCombinations::Array{Tuple{Symbol, Symbol}, 1}=Tuple{Symbol, Symbol}[],#(:Mix, :safeMinimumRuleCurve)],
                saveFolder::AbstractString="", savePlots::Bool=false, saveIntermediate::Bool=false,
                verbosity::Int=1)
  # Some assumptions on the durations.
  @assert typeof(ruleDuration) <: Year
  @assert typeof(guarantee) <: Year

  # Ensure the output folder is ready.
  if length(saveFolder) > 0 && ! isdir(saveFolder)
    mkdir(saveFolder)
  end

  ### Generate parameters.
  hasCurrentRuleCurve = ! isempty(currentRuleCurve)

  generateDefaultLengthScenarios = in(:minimumRuleCurve, solvers) # Default length scenarios (i.e. guarantee).
  generateExtraLengthScenarios = in(:safeMinimumRuleCurve, solvers) # Longer scenarios (i.e. ruleDuration + guarantee).
  scenariosDefaultLength = Int(guarantee)
  scenariosExtraLength = Int(ruleDuration) + Int(guarantee)

  ruleLengthPeriods = Int(floor(Dates.days(ruleDuration) / Dates.days(getPeriod(reservoir))))
  guaranteePeriods = Int(floor(Dates.days(guarantee) / Dates.days(getPeriod(reservoir))))

  generateCIScenarios = length(ciAlgorithmsLevels) > 0
  ciAlgorithms = [algo[1] for algo in ciAlgorithmsLevels]

  # TODO: COMPLETE REFACTORING OF THE WHOLE FUNCTION! TODO TODO TODO TODO TODO TODO TODO TODO TODO TODO TODO TODO TODO TODO TODO TODO TODO TODO TODO TODO TODO TODO TODO TODO TODO TODO TODO TODO TODO TODO TODO TODO
  # A complete data structure of available solvers, presented with a uniform syntax:
  #     (reservoir, verbosity) -> f(...)::ReservoirSolution
  # TODO: Not yet useable throughout this script (just for cross-validation): for the solving process, needs information for
  # TODO:   the solution name (could just currify the functions for cross-validation and let that string as a parameter)
  # TODO:       (reservoir, verbosity, name) -> f(...)::ReservoirSolution
  # TODO:   however, this is still not enough, as minimumRuleCurve and safeMinimumRuleCurve do not take their scenarios from
  # TODO:   the same place (scgen_default, scgen_extra). Would need a more complex data structure to hold the solvers?
  # solversF = Dict{Symbol, Function}
  # for solver in solvers
  #   if solver == :minimumRuleCurve
  #     solversF[solver] = (sc, verbosity) -> minimumRuleCurve(sc, "", outputLength=ruleLengthPeriods, verbosity=verbosity)
  #   elseif solver == :safeMinimumRuleCurve
  #     solversF[solver] = (sc, verbosity) -> safeMinimumRuleCurve(sc, ruleLengthPeriods, guaranteePeriods, "", verbosity=verbosity)
  #   else
  #     error("Unrecognised solver: " * string(solver))
  #   end
  # end
  # TODO: COMPLETE REFACTORING OF THE WHOLE FUNCTION! TODO TODO TODO TODO TODO TODO TODO TODO TODO TODO TODO TODO TODO TODO TODO TODO TODO TODO TODO TODO TODO TODO TODO TODO TODO TODO TODO TODO TODO TODO TODO TODO

  ### Handle scenario generation.
  if verbosity >= 0; println("[=*=] Scenario generation: started... "); end
  cis = Dict{Tuple{Symbol, Float64}, Tuple{Reservoir, Reservoir, Reservoir}}()
  scgen_default = Dict{Symbol, Reservoir}()
  scgen_extra = Dict{Symbol, Reservoir}()

  # Confidence interval computations if needed.
  for algo in ciAlgorithms
    for level in ciAlgorithmsLevels[algo]
      cis[(algo, level)] = computeCI(reservoir, level, ci_period=ciPeriod, method=algo)
    end
  end

  # Actual scenario generation.
  if generateDefaultLengthScenarios
    for algo in scenarioGenerationAlgorithms
      scgen_default[algo] = scenarioGeneration(reservoir, algo, scenariosDefaultLength)
    end
  end

  if generateExtraLengthScenarios
    for algo in scenarioGenerationAlgorithms
      scgen_extra[algo] = scenarioGeneration(reservoir, algo, scenariosExtraLength)
    end
  end

  if generateCIScenarios
    for algo in ciAlgorithms
      for level in ciAlgorithmsLevels[algo]
        s = symbol("ci_", algo, "_", level)
        if generateDefaultLengthScenarios
          scgen_default[s] = scenarioGeneration(cis[(algo, level)][2], :Duplicate, scenariosDefaultLength)
        end
        if generateExtraLengthScenarios
          scgen_extra[s] = scenarioGeneration(cis[(algo, level)][2], :Duplicate, scenariosExtraLength)
        end
      end
    end
  end

  if verbosity >= 0; println("[=*=] Scenario generation: done! "); end

  ### Handle solving.
  if verbosity >= 0; println("[=*=] Solving: started... "); end
  solutions = Dict{Tuple{Symbol, Symbol}, ReservoirSolution}()
  if hasCurrentRuleCurve
    solutions[(:currentRuleCurve, :currentRuleCurve)] = wrapExistingSolution(reservoir, currentRuleCurve)
  end

  for solver in solvers
    if verbosity >= 0; println("[=*=]   Solving with the algorithm: " * string(solver)); end

    if solver == :minimumRuleCurve
      for (lbl, sc) in scgen_default
        if ! in((lbl, solver), disabledCombinations)
          if verbosity >= 0; println("[=*=]     Solving with the data set: " * string(lbl)); end
          sol_lbl = "stochastic, default solver, " * string(guarantee) * ", " * string(lbl)
          solutions[(solver, lbl)] = minimumRuleCurve(sc, sol_lbl, outputLength=ruleLengthPeriods, useOutputs=getOutputIds(sc), verbosity=verbosity)
        end
      end
    elseif solver == :safeMinimumRuleCurve
      for (lbl, sc) in scgen_extra
        if ! in((lbl, solver), disabledCombinations)
          if verbosity >= 0; println("[=*=]     Solving with the data set: " * string(lbl)); end
          sol_lbl = "stochastic, safety solver, " * string(guarantee) * " + " * string(ruleDuration) * ", " * string(lbl)
          solutions[(solver, lbl)] = safeMinimumRuleCurve(sc, ruleLengthPeriods, guaranteePeriods, sol_lbl, verbosity=verbosity)
        end
      end
    else
      error("Unrecognised solver: " * string(solver))
    end
  end

  if verbosity >= 0; println("[=*=] Solving: done! "); end

  ## Check the consistency of the solutions.
  # All must have the same number of timesteps.
  sol_n_ts = -1 # Will be initialised to a positive value when the first feasible solution is found.
  sol_ts_pb = Dict{Any, Int}()
  for (desc, solution) in solutions
    if solution.feasible
      l = countTimeSteps(solution)
      if sol_n_ts < 0
        sol_n_ts = l
      elseif sol_n_ts != l
        sol_ts_pb[desc] = l
      end
    end
  end
  if sol_n_ts > 0 && length(sol_ts_pb) > 0
    txt_sol = "A solution has"
    if length(sol_ts_pb) > 1
      txt_sol = string(length(sol_ts_pb)) * " solutions have"
    end
    error(txt_sol * " an inconsistent number of time steps: " * string(sol_ts_pb))
  end

  # Other tests to perform?

  if verbosity >= 1; println("[=*=] Solving: all solutions are consistent! "); end

  ### Handle evaluation.
  if verbosity >= 0; println("[=*=] Evaluating: started... "); end

  evaluations = Dict{Tuple{Tuple{Symbol, Symbol}, Symbol}, ReservoirSolution}() # ((solver, scenarioGenerator), evaluator) -> evaluation
  # Associates (solutionAlgorithm, evaluator) to the result of an evaluation.
  # A solutionAlgorithm is a pair (solver, scenarioGenerator).
  for evaluator in evaluators
    if verbosity >= 0; println("[=*=]   Evaluating with the algorithm: " * string(evaluator)); end

    for (desc, solution) in solutions
      if ! solution.feasible
        continue
      end

      solver = desc[1]
      sc_gen = desc[2]

      if verbosity >= 0; println("[=*=]     Evaluating the solution: solver " * string(solver) * ", scenario generation " * string(sc_gen)); end

      if evaluator == :feasibilityProbabilitySyntheticScenarios
        evaluations[(desc, evaluator)] = feasibilityProbability(solution, numberScenarios=1000, scenarioDuration=guarantee, useOutputs=getOutputIds(getReservoir(solution)), verbosity=verbosity)
      elseif evaluator == :feasibilityProbabilityAllDatasets
        if verbosity >= 0; println("[=*=]     Evaluating with default length data (" * string(guarantee) * ")."); end
        for (lbl, sc) in scgen_default
         evaluations[(desc, symbol(evaluator, "_default_", lbl))] = feasibilityProbability(solution, scenarios=sc, useOutputs=getOutputIds(getReservoir(solution)), verbosity=verbosity)
        end
        if verbosity >= 0; println("[=*=]     Evaluating with extra length data (" * string(ruleDuration + guarantee) * ")."); end
        for (lbl, sc) in scgen_extra
          evaluations[(desc, symbol(evaluator, "_extra_", lbl))] = feasibilityProbability(solution, scenarios=sc, useOutputs=getOutputIds(getReservoir(solution)), verbosity=verbosity)
        end
      elseif evaluator == :feasibilityProbabilityClimateChange
        if verbosity >= 0; println("[=*=]     Evaluating with climate change scenarios."); end
        if isempty(climateChangeScenarios); error("No climate change scenarios given, while the script was asked to test against them!"); end;
        for (lbl, sc) in climateChangeScenarios
          evaluations[(desc, symbol(evaluator, "_", lbl))] = feasibilityProbability(solution, scenarios=sc, useOutputs=getOutputIds(getReservoir(solution)), verbosity=verbosity)
        end
      elseif evaluator == :crossValidatedFeasibilityProbability
        error("Cross-validation is not a solution evaluator! ")
      elseif evaluator == :worstCaseFunctions
        evaluations[(desc, evaluator)] = worstCaseFunctions(solution, useOutputs=getOutputIds(getReservoir(solution)), verbosity=verbosity)
      elseif evaluator == :maximumDeviation
        inflow = vec(getSolution(worstCaseFunctions(solutions[(:currentRuleCurve, :currentRuleCurve)], useOutputs=getOutputIds(getReservoir(solution)), verbosity=verbosity)))
        evaluations[(desc, evaluator)] = maximumDeviation(solution, inflow, verbosity=verbosity)
      elseif evaluator == :hydropower
        if solver == :minimumRuleCurve
          price_scenario = [sin(2 * π / 20 * t) for t in 1:scenariosDefaultLength]
        elseif solver == :safeMinimumRuleCurve
          price_scenario = [sin(2 * π / 30 * t) for t in 1:scenariosExtraLength]
        else
          error("Unrecognised solver for generating price scenarios: " * string(solver))
        end
        evaluations[(desc, evaluator)] = hydropower(solution, useOutputs=getOutputIds(getReservoir(solution)), electricityPrices=price_scenario)
      else
        error("Unrecognised evaluator: " * string(evaluator))
      end
    end
  end

  globalEvaluations = Dict{Tuple{Symbol, Symbol, Symbol}, ReservoirSolution}() # (solver, sc_gen, globalEvaluator) -> evaluation result
  for globalEvaluator in globalEvaluators
    if verbosity >= 0; println("[=*=]   Globally evaluating with the algorithm: " * string(globalEvaluator)); end

    for solver in solvers
      if verbosity >= 0; println("[=*=]     Evaluating solver with " * string(globalEvaluator) * ": " * string(solver)); end

      if globalEvaluator == :crossValidatedFeasibilityProbability
        if solver == :minimumRuleCurve
          for sc_gen in scenarioGenerationAlgorithms
            if in((sc_gen, solver), disabledCombinations)
              continue
            end

            if verbosity >= 0; println("[=*=]     Evaluating solver with " * string(globalEvaluator) * ": " * string(solver) * ", scenario generation " * string(sc_gen)); end

            f = function (sc, verbosity)
              scgen = scenarioGeneration(sc, sc_gen, scenariosDefaultLength)
              return minimumRuleCurve(scgen, "", useOutputs=getOutputIds(sc), outputLength=ruleLengthPeriods, verbosity=verbosity)
            end
            globalEvaluations[(solver, sc_gen, globalEvaluator)] = crossValidatedFeasibilityProbability(reservoir, f, verbosity=verbosity, folds=10)
          end
          for algo in ciAlgorithms
            for level in ciAlgorithmsLevels[algo]
              sc_gen = symbol("ci_", algo, "_", level)

              if in((sc_gen, solver), disabledCombinations)
                continue
              end

              if verbosity >= 0; println("[=*=]     Evaluating solver with " * string(globalEvaluator) * ": " * string(solver) * ", scenario generation " * string(sc_gen)); end

              f = function (sc, verbosity)
                cis = computeCI(sc, level, ci_period=ciPeriod, method=algo)[2]
                scgen = scenarioGeneration(cis, :Duplicate, scenariosDefaultLength)
                return minimumRuleCurve(scgen, "", useOutputs=getOutputIds(sc), outputLength=ruleLengthPeriods, verbosity=verbosity)
              end
              globalEvaluations[(solver, sc_gen, globalEvaluator)] = crossValidatedFeasibilityProbability(reservoir, f, verbosity=verbosity, folds=10)
            end
          end
        elseif solver == :safeMinimumRuleCurve
          for sc_gen in scenarioGenerationAlgorithms
            if in((sc_gen, solver), disabledCombinations)
              continue
            end

            if verbosity >= 0; println("[=*=]     Evaluating solver with " * string(globalEvaluator) * ": " * string(solver) * ", scenario generation " * string(sc_gen)); end

            f = function (sc, verbosity)
              scgen = scenarioGeneration(sc, sc_gen, scenariosExtraLength)
              return safeMinimumRuleCurve(scgen, ruleLengthPeriods, guaranteePeriods, "", verbosity=verbosity)
            end
            globalEvaluations[(solver, sc_gen, globalEvaluator)] = crossValidatedFeasibilityProbability(reservoir, f, verbosity=verbosity, folds=10)
          end
          for algo in ciAlgorithms
            for level in ciAlgorithmsLevels[algo]
              sc_gen = symbol("ci_", algo, "_", level)

              if in((sc_gen, solver), disabledCombinations)
                continue
              end

              if verbosity >= 0; println("[=*=]     Evaluating solver with " * string(globalEvaluator) * ": " * string(solver) * ", scenario generation " * string(sc_gen)); end

              f = function (sc, verbosity)
                cis = computeCI(sc, level, ci_period=ciPeriod, method=algo)[2]
                scgen = scenarioGeneration(cis, :Duplicate, scenariosExtraLength)
                return safeMinimumRuleCurve(scgen, ruleLengthPeriods, guaranteePeriods, "", verbosity=verbosity)
              end
              globalEvaluations[(solver, sc_gen, globalEvaluator)] = crossValidatedFeasibilityProbability(reservoir, f, verbosity=verbosity, folds=10)
            end
          end
        else
          error("Unrecognised solver: " * string(solver))
        end
      else
        error("Unrecognised global evaluator: " * string(solver))
      end
    end
  end

  if verbosity >= 0; println("[=*=] Evaluating: done! "); end

  ### Finally, display report, generate figures, and save what should be saved.
  if verbosity >= 1
    for desc in keys(solutions)
      println("**************************************************************************************************************")
      println("** Solver: " * string(desc[1]))
      println("** Scenario generation: " * string(desc[2]))
      println("**")
      for evaluator in evaluators
        # The generated keys for those solvers are not detected by the following condition!
        if evaluator == :feasibilityProbabilityAllDatasets
          for (eval_desc, evaluation) in evaluations
            if desc == eval_desc[1] && startswith(string(eval_desc[2]), "feasibilityProbabilityAllDatasets")
              ev = evaluation.solution
              first_part_length = length("feasibilityProbabilityAllDatasets_")
              dataset = string(eval_desc[2])[first_part_length + 1:end]
              println("** Feasibility analysis on data set " * dataset * ": probability of being infeasible: " * string(ev.probabilityInfeasible * 100) * "%")
              println("** Feasibility analysis on data set " * dataset * ": probability of being infeasible due to water surplus: " * string(ev.probabilityInfeasibleSurplus * 100) * "%")
              println("** Feasibility analysis on data set " * dataset * ": probability of being infeasible due to water shortage: " * string(ev.probabilityInfeasibleShortage * 100) * "%")
              println("** Feasibility analysis on data set " * dataset * ": average water shortage: " * string(ev.averageShortage) * " 10^6 m^3")
              println("** Feasibility analysis on data set " * dataset * ": average shortage index: " * string(ev.averageShortageIndex))
              println("** Feasibility analysis on data set " * dataset * ": average shortage vulnerability index: " * string(ev.averageShortageVulnerability))
            end
          end
        end
        if evaluator == :feasibilityProbabilityClimateChange
          for (eval_desc, evaluation) in evaluations
            if desc == eval_desc[1] && startswith(string(eval_desc[2]), "feasibilityProbabilityClimateChange")
              ev = evaluation.solution
              first_part_length = length("feasibilityProbabilityClimateChange_")
              dataset = string(eval_desc[2])[first_part_length + 1:end]
              println("** Feasibility analysis on the climate change scenario " * dataset * ": probability of being infeasible: " * string(ev.probabilityInfeasible * 100) * "%")
              println("** Feasibility analysis on the climate change scenario " * dataset * ": probability of being infeasible due to water surplus: " * string(ev.probabilityInfeasibleSurplus * 100) * "%")
              println("** Feasibility analysis on the climate change scenario " * dataset * ": probability of being infeasible due to water shortage: " * string(ev.probabilityInfeasibleShortage * 100) * "%")
              println("** Feasibility analysis on the climate change scenario " * dataset * ": average water shortage: " * string(ev.averageShortage) * " 10^6 m^3")
              println("** Feasibility analysis on the climate change scenario " * dataset * ": average shortage index: " * string(ev.averageShortageIndex))
              println("** Feasibility analysis on the climate change scenario " * dataset * ": average shortage vulnerability index: " * string(ev.averageShortageVulnerability))
            end
          end
        end

        if haskey(evaluations, (desc, evaluator)) # Not all evaluators could work previously (e.g. cross-validation needs more scenarios than available with CIs, with more parameters than a sheer symbol).
          ev = evaluations[(desc, evaluator)].solution
          if evaluator == :feasibilityProbabilitySyntheticScenarios
            println("** Monte-Carlo feasibility analysis: probability of being infeasible: " * string(ev.probabilityInfeasible * 100) * "%")
            println("** Monte-Carlo feasibility analysis: probability of being infeasible due to water surplus: " * string(ev.probabilityInfeasibleSurplus * 100) * "%")
            println("** Monte-Carlo feasibility analysis: probability of being infeasible due to water shortage: " * string(ev.probabilityInfeasibleShortage * 100) * "%")
            println("** Monte-Carlo feasibility analysis: average water shortage: " * string(ev.averageShortage) * " 10^6 m^3")
            println("** Monte-Carlo feasibility analysis: average shortage index: " * string(ev.averageShortageIndex))
            println("** Monte-Carlo feasibility analysis: average shortage vulnerability index: " * string(ev.averageShortageVulnerability))
          elseif evaluator == :worstCaseFunctions
            println("** Inverse model analysis: annual volume below which infeasibility is guaranteed: " * string(ev.total) * " 10^6 m^3")
          elseif evaluator == :maximumDeviation
            println("** Maximum deviation: from the inverse solution for the current rule curve: " * string(ev.deviation) * " 10^6 m^3")
          else
            error("Unrecognised evaluator: " * string(evaluator))
          end
        end
      end
      println("**************************************************************************************************************")
    end

    for globalEvaluator in globalEvaluators
      if globalEvaluator == :crossValidatedFeasibilityProbability
        for (globalEvaluationDescription, globalEvaluation) in globalEvaluations
          solver = globalEvaluationDescription[1]
          sc_gen = globalEvaluationDescription[2]
          evaluator = globalEvaluationDescription[3]
          if evaluator == :crossValidatedFeasibilityProbability
            ev = globalEvaluation.solution

            println("**************************************************************************************************************")
            println("** Solver: ", solver)
            println("** Scenario generation: ", sc_gen)
            println("** ")
            println("** Cross-validation feasibility analysis: probability of being infeasible: " * string(ev.probabilityInfeasible * 100) * "%")
            println("** Cross-validation feasibility analysis: probability of being infeasible due to water surplus: " * string(ev.probabilityInfeasibleSurplus * 100) * "%")
            println("** Cross-validation feasibility analysis: probability of being infeasible due to water shortage: " * string(ev.probabilityInfeasibleShortage * 100) * "%")
            println("** Cross-validation feasibility analysis: probability of having no rule curve: " * string(ev.probabilityNoRuleCurve * 100) * "%")
            println("** Cross-validation feasibility analysis: average water shortage: " * string(ev.averageShortage) * " 10^6 m^3")
            println("** Cross-validation feasibility analysis: average shortage index: " * string(ev.averageShortageIndex))
            println("** Cross-validation feasibility analysis: average shortage vulnerability index: " * string(ev.averageShortageVulnerability))
            println("**************************************************************************************************************")
          end
        end
      end
    end
  end

  if saveIntermediate
    if verbosity >= 0; println("[=*=] Serialising: started... "); end

    folder_intermediate = saveFolder * "/intermediate/"
    if ! isdir(folder_intermediate)
      mkdir(folder_intermediate)
    end
    serialise_own(filename, o) = open(folder_intermediate * filename, "w") do stream; serialize(stream, o); end

    serialise_own("cis.jls", cis)
    serialise_own("scgen_default.jls", scgen_default)
    serialise_own("scgen_extra.jls", scgen_extra)
    serialise_own("solutions.jls", solutions)
    serialise_own("evaluations.jls", evaluations)
    serialise_own("global_evaluations.jls", globalEvaluations)

    if verbosity >= 0; println("[=*=] Serialising: done! "); end
  end

  if savePlots
    if verbosity >= 0; println("[=*=] Plotting: started... "); end

    folder_plots = saveFolder * "/plots/"
    if ! isdir(folder_plots)
      mkdir(folder_plots)
    end

    function exportPlot(plot::AbstractPlot, name::AbstractString)
      savefig(plot, folder_plots * name * ".png")
      savefig(plot, folder_plots * name * ".svg")
      nothing
    end

    function exportData(data::Array{Any, 2}, name::AbstractString)
      writedlm(folder_plots * name * ".tsv", data, '\t')
      nothing
    end

    period = getPeriod(reservoir)

    ## Data treatment.
    # Confidence intervals. First plots for each algorithm, then a global plot.
    # Indices in the code: once the average [1], then for each pair upper/lower bounds [2] and [3].
    if ! isempty(currentRuleCurve) && length(ciAlgorithms) > 0
      # Plot each algorithm independently (if at least two algorithms), all confidence levels at once.
      if length(ciAlgorithms) > 1
        for algo in ciAlgorithms
          # All confidence levels at once.
          data = vcat(cis[(algo, ciAlgorithmsLevels[algo][1])][1], [let lcis = cis[(algo, level)]; [lcis[2], lcis[3]]; end for level in ciAlgorithmsLevels[algo]]...)
          plt, dat = plotRiverAlongReservoirs(data, period, 1:ruleLengthPeriods, getName(reservoir), exportPlotData=true)
          exportPlot(plt, "confidence_intervals_" * string(algo))
          exportData(dat, "confidence_intervals_" * string(algo))
        end
      end

      # Each confidence level in its own plot (if more than one level requested), for each algorithm.
      for algo in ciAlgorithms
        if length(ciAlgorithmsLevels[algo]) > 1
          for level in ciAlgorithmsLevels[algo]
            data = vcat(cis[(algo, ciAlgorithmsLevels[algo][1])][1], cis[(algo, level)][2], cis[(algo, level)][3])
            plt, dat = plotRiverAlongReservoirs(data, period, 1:ruleLengthPeriods, getName(reservoir), exportPlotData=true)
            exportPlot(plt, "confidence_intervals_" * string(algo) * "_cl" * string(level))
            exportData(dat, "confidence_intervals_" * string(algo) * "_cl" * string(level))
          end
        end
      end

      # TODO: plot for a given confidence level (irrespectively of algorithms).

      # Global plot for all algorithms and confidence levels at once.
      data = vcat([[cis[(algo, ciAlgorithmsLevels[algo][1])][1]; vcat([let lcis = cis[(algo, level)]; [lcis[2]; lcis[3]]; end for level in ciAlgorithmsLevels[algo]]...)] for algo in ciAlgorithms]...)
      plt, dat = plotRiverAlongReservoirs(data, period, 1:ruleLengthPeriods, getName(reservoir), exportPlotData=true)
      exportPlot(plt, "confidence_intervals")
      exportData(dat, "confidence_intervals")
    end

    # Skip scenario generation output.

    ## Solutions.
    # Existing target, alone.
    if ! isempty(currentRuleCurve)
      plt, dat = plotRule(currentRuleCurve, period, title="Current rule curve", exportPlotData=true)
      exportPlot(plt, "current_rule")
      exportData(dat, "current_rule")
    end

    # Complete set of scenarios when merging them and compare them to the solution.
    if in(:Merge, scenarioGenerationAlgorithms) && in(:minimumRuleCurve, solvers) && ! isempty(currentRuleCurve)
      merge_sol = solutions[(:minimumRuleCurve, :Merge)]
      colourings = [:YearlyAverage, :WetSeasonAverage, :DrySeasonAverage, :DriestMonth, :DriestThreeMonths, :DriestSixMonths]

      for colouring in colourings
        plt, dat = plotSolutionVsScenarios(merge_sol, currentRuleCurve, colouring=colouring, exportPlotData=true)
        exportPlot(plt, "solution_safety_stochastic_merge_splitperscenario_" * string(colouring))
        exportData(dat, "solution_safety_stochastic_merge_splitperscenario_" * string(colouring))
      end
    end

    # Solutions for the various solvers: one plot for each solver first, then one plot with each and every solution.
    if ! isempty(currentRuleCurve)
      for solver in solvers
        solver_solution = filter((k, v) -> k[1] == solver, solutions)
        plt, dat = plotSolutions(collect(values(solver_solution)), currentRuleCurve, exportPlotData=true)
        exportPlot(plt, "solutions_" * string(solver))
        exportData(dat, "solutions_" * string(solver))
      end

      plt, dat = plotSolutions(collect(values(solutions)), currentRuleCurve, exportPlotData=true)
      exportPlot(plt, "solutions")
      exportData(dat, "solutions")
    end

    ## Evaluation.
    # Nothing for solution evaluators.
    # Nothing for global evaluators.

    if verbosity >= 0; println("[=*=] Plotting: done! "); end
  end

  return cis, scgen_default, scgen_extra, solutions, evaluations
end


text:
function ralston3(f, t, y, dt, P)
    k1 = f(t, y, P)
    k2 = f(t + 1/2*dt, y + 1/2*dt*k1, P)
    k3 = f(t + 3/4*dt, y + 3/4*dt*k2, P)
    y + dt*(2/9*k1 + 1/3*k2 + 4/9*k3)
end

struct Ralston3 <: ODESolverType end

#=
        FOUR components
=#
function update(::Ralston3, fs, t, A, B, C, D, dt, P, tableau::Nothing=nothing)
    kA1 = update(fs[1], t, A, B, C, D, P)
    kB1 = update(fs[2], t, A, B, C, D, P)
    kC1 = update(fs[3], t, A, B, C, D, P)
    kD1 = update(fs[4], t, A, B, C, D, P)

    Ai = A + 1/2*dt*kA1
    Bi = B + 1/2*dt*kB1
    Ci = C + 1/2*dt*kC1
    Di = D + 1/2*dt*kD1

    kA2 = update(fs[1], t + 1/2*dt, Ai, Bi, Ci, Di, P)
    kB2 = update(fs[2], t + 1/2*dt, Ai, Bi, Ci, Di, P)
    kC2 = update(fs[3], t + 1/2*dt, Ai, Bi, Ci, Di, P)
    kD2 = update(fs[4], t + 1/2*dt, Ai, Bi, Ci, Di, P)

    Ai = A + 3/4*dt*kA2
    Bi = B + 3/4*dt*kB2
    Ci = C + 3/4*dt*kC2
    Di = D + 3/4*dt*kD2

    kA3 = update(fs[1], t + 3/4*dt, Ai, Bi, Ci, Di, P)
    kB3 = update(fs[2], t + 3/4*dt, Ai, Bi, Ci, Di, P)
    kC3 = update(fs[3], t + 3/4*dt, Ai, Bi, Ci, Di, P)
    kD3 = update(fs[4], t + 3/4*dt, Ai, Bi, Ci, Di, P)

    (A + dt*(2/9*kA1 + 1/3*kA2 + 4/9*kA3), B + dt*(2/9*kB1 + 1/3*kB2 + 4/9*kB3),
     C + dt*(2/9*kC1 + 1/3*kC2 + 4/9*kC3), D + dt*(2/9*kD1 + 1/3*kD2 + 4/9*kD3))
end

#=
        THREE components
=#
function update(::Ralston3, fs, t, A, B, C, dt, P, tableau::Nothing=nothing)
    kA1 = update(fs[1], t, A, B, C, P)
    kB1 = update(fs[2], t, A, B, C, P)
    kC1 = update(fs[3], t, A, B, C, P)

    Ai = A + 1/2*dt*kA1
    Bi = B + 1/2*dt*kB1
    Ci = C + 1/2*dt*kC1

    kA2 = update(fs[1], t + 1/2*dt, Ai, Bi, Ci, P)
    kB2 = update(fs[2], t + 1/2*dt, Ai, Bi, Ci, P)
    kC2 = update(fs[3], t + 1/2*dt, Ai, Bi, Ci, P)

    Ai = A + 3/4*dt*kA2
    Bi = B + 3/4*dt*kB2
    Ci = C + 3/4*dt*kC2

    kA3 = update(fs[1], t + 3/4*dt, Ai, Bi, Ci, P)
    kB3 = update(fs[2], t + 3/4*dt, Ai, Bi, Ci, P)
    kC3 = update(fs[3], t + 3/4*dt, Ai, Bi, Ci, P)

    (A + dt*(2/9*kA1 + 1/3*kA2 + 4/9*kA3), B + dt*(2/9*kB1 + 1/3*kB2 + 4/9*kB3),
     C + dt*(2/9*kC1 + 1/3*kC2 + 4/9*kC3))
end

#=
        TWO components
=#
function update(::Ralston3, fs, t, A, B, dt, P, tableau::Nothing=nothing)
    kA1 = update(fs[1], t, A, B, P)
    kB1 = update(fs[2], t, A, B, P)

    Ai = A + 1/2*dt*kA1
    Bi = B + 1/2*dt*kB1

    kA2 = update(fs[1], t + 1/2*dt, Ai, Bi, P)
    kB2 = update(fs[2], t + 1/2*dt, Ai, Bi, P)

    Ai = A + 3/4*dt*kA2
    Bi = B + 3/4*dt*kB2

    kA3 = update(fs[1], t + 3/4*dt, Ai, Bi, P)
    kB3 = update(fs[2], t + 3/4*dt, Ai, Bi, P)

    (A + dt*(2/9*kA1 + 1/3*kA2 + 4/9*kA3), B + dt*(2/9*kB1 + 1/3*kB2 + 4/9*kB3))
end

ext: jl


text:
using Statistics
using Photometry.Aperture
using FillArrays
using Rotations
using CoordinateTransformations
using ImageTransformations
using Interpolations
using StaticArrays

using PSFModels: PSFModel

export normalize_psf, normalize_psf!

const PSFType = Union{AbstractExtrapolation,PSFModel}

struct CubeGenerator{CT<:AbstractArray,PT<:PSFType,AT<:AbstractVector}
    cube::CT
    angles::AT
    psf::PT
end

"""
    CubeGenerator(cube, angles, psf; degree=Linear(), fill=0)

Creates a generator which can rapidly create synthetic data using a PSF. The `cube` and `angles` are used to define the geometry of the synthetic data cube. If the `psf` is a dense matrix, an extrapolator will be created with the given `degree` and `fill` value.

# Examples

```jldoctest
julia> using HCIDatasets: BetaPictoris

julia> cube, angles, psf = BetaPictoris[:cube, :pa, :psf];

julia> gen = CubeGenerator(cube, angles, psf); # using empirical PSF

julia> using PSFModels: Gaussian

julia> gen2 = CubeGenerator(cube, angles, Gaussian(4.7)); # using PSFModel
```
"""
function CubeGenerator(cube, angles, psf; kwargs...)
    if size(cube, 1) != length(angles)
        error("Number of frames in cube does not match number of parallactic angles")
    end
    psf′ = prep_psf(psf; kwargs...)
    return CubeGenerator(cube, angles, psf′)
end

function prep_psf(psf::AbstractMatrix{T}; degree=Linear(), fill=zero(T)) where T
    return ImageTransformations.box_extrapolation(psf, degree, fill)
end
prep_psf(psf::PSFModel; kwargs...) = psf

"""
    (::CubeGenerator)([T], pos; A=1)
    (::CubeGenerator)(base, pos; A=1)

Use the cube generator to generate a synthetic cube. The eltype can be manually specified with `T`, or otherwise a base can be given. If a base is given, the synthetic cube will be *injected*, which will add it to the `base`. `base` can be a matrix, in which case the generator will automatically flatten the output without extra allocations.

In the case that the generator is build on top of a geometry like [`AnnulusView`](@ref), using a matrix `base` will use the *filtered* matrix in order to avoid calculating out of bounds indices.

# Examples

```jldoctest
julia> using HCIDatasets: BetaPictoris

julia> cube, angles, psf = BetaPictoris[:cube, :pa, :psf];

julia> gen = CubeGenerator(cube, angles, psf); # using empirical PSF

julia> s = gen((51, 51)); # cube with PSF at center

julia> size(s) == size(cube)
true

julia> flat_s = gen(zero(flatten(cube)), (51, 51));

julia> flat_s ≈ flatten(s)
true

julia> gen(Polar(10, deg2rad(45)); A=10); # inject using Polar coords
```
"""
(gen::CubeGenerator)(pos; kwargs...) = gen(eltype(gen.cube), pos; kwargs...)

function (gen::CubeGenerator)(T::Type, pos; kwargs...)
    base = fill!(similar(gen.cube, T), zero(T))
    gen(base, pos; kwargs...)
end

function (gen::CubeGenerator)(base::AbstractArray{T,3}, pos; A=one(T)) where T
    ctr = reverse(center(base)[2:3])
    Threads.@threads for tidx in axes(base, 1)
        location = parse_position(pos, ctr, gen.angles[tidx])

        frame = @view base[tidx, :, :]
        tform = Translation(center(gen.psf) - location)
        tform = ImageTransformations.try_static(tform, frame)
        for idx in CartesianIndices(frame)
            pidx = tform(SVector(idx.I))
            if gen.psf isa PSFModel
                frame[idx] += A * gen.psf(Tuple(reverse(pidx))...)
            else
                frame[idx] += A * gen.psf(Tuple(pidx)...)
            end
        end
    end
    return base
end

function (gen::CubeGenerator)(base::AbstractMatrix{T}, pos; A=one(T)) where T
    ctr = reverse(center(gen.cube)[2:3])
    ny, nx = size(gen.cube)[2:3]
    Threads.@threads for tidx in axes(base, 1)
        location = parse_position(pos, ctr, gen.angles[tidx])
        tform = Translation(center(gen.psf) - location)
        for pidx in axes(base, 2)
            j = (pidx - 1) % ny + 1
            k = (pidx - 1) ÷ nx + 1
            I = tform(SVector(j, k))
            if gen.psf isa PSFModel
                base[tidx, pidx] += A * gen.psf(Tuple(reverse(I))...)
            else
                base[tidx, pidx] += A * gen.psf(Tuple(I)...)
            end
        end
    end
    return base
end

function (gen::CubeGenerator{<:AnnulusView})(base::AbstractArray{T,3}, pos; A=one(T)) where T
    mat = @view base[gen.cube.indices...]
    gen(mat, pos; A=A)
    return base
end

function (gen::CubeGenerator{<:AnnulusView})(base::AbstractMatrix{T}, pos; A=one(T)) where T
    ctr = reverse(center(gen.cube)[2:3])
    Threads.@threads for tidx in axes(base, 1)
        tidx′ = gen.cube.indices[1][tidx]
        location = parse_position(pos, ctr, gen.angles[tidx′])
        tform = Translation(center(gen.psf) - location)
        for (pidx, pidx′) in zip(axes(base, 2), gen.cube.indices[2])
            I = tform(SVector(pidx′.I))
            if gen.psf isa PSFModel
                base[tidx, pidx] += A * gen.psf(Tuple(reverse(I))...)
            else
                base[tidx, pidx] += A * gen.psf(Tuple(I)...)
            end
        end
    end
    return base
end

"""
    HCIToolbox.parse_position(pos::Tuple, ctr, pa=0)

Parse the position as the cartesian `(x, y)` pixel coordinates. If a parallactic angle is given the position will be rotated `pa` degrees clockwise.
"""
function parse_position(pos::Tuple, ctr, pa=0)
    xy = SVector(pos)
    ϕ = RotMatrix{2}(-deg2rad(pa))
    return recenter(ϕ, ctr)(xy) |> reverse
end

"""
    HCIToolbox.parse_position(pos::Polar, ctr, pa=0)

Parse the position as the polar `(r, θ)` pixel coordinates centered at `ctr`. If a parallactic angle is given the position will be rotated `pa` degrees clockwise.
"""
function parse_position(pos::Polar, ctr, pa=0)
    new_pos = Polar(pos.r, pos.θ - deg2rad(pa))
    Δxy = convert(SVector, new_pos)
    # recenter and switch to index layout
    return  Δxy |> Translation(ctr) |> reverse
end

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

"""
    inject(frame, psf, position; A=1, degree=Linear(), fill=0.0)

Injects the matrix `psf` into `frame` at the given position. If `position` is a Tuple, it will be parsed as pixel coordinates `(x, y)`. If `position` is a `Polar` (from CoordinateTransformations.jl) it will be parsed as polar coordinates centered on the frame center.

The amplitude of the injected frame can be changed with `A`, and in the case that `psf` is dense matrix, and extrapolator will be created with the given `degree` and `fill` value.

### Coordinate System
* Cartesian - Parsed as distance from the bottom-left corner of the image. Pixel convention is that `(1, 1)` is the center of the bottom-left pixel increasing right and up.
* Polar - Parsed as polar coordinates from the center of the image. `theta` is a position angle.

!!! note
    Due to the integral nature of array indices, frames or images with even-sized axes will have their center rounded to the nearest integer. This may cause unexpected results for small frames and images.

# Examples
```jldoctest
julia> inject(zeros(5, 5), ones(1, 1), (2, 1); A=2) # image coordinates (x=2, y=1)
5×5 Array{Float64,2}:
 0.0  2.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0

julia> inject(zeros(5, 5), ones(3, 3), Polar(1.5, deg2rad(90))) # polar coords
5×5 Array{Float64,2}:
 0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0
 0.0  1.0  1.0  1.0  0.0
 0.0  1.0  1.0  1.0  0.0
```
"""
inject(frame::AbstractMatrix, args...; kwargs...) = inject!(copy(frame), args...; kwargs...)

function inject!(frame::AbstractMatrix, kernel::AbstractMatrix{T}, pos; degree=Linear(), fill=zero(T), kwargs...) where T
    etp = ImageTransformations.box_extrapolation(kernel, degree, fill)
    return inject!(frame, etp, pos; kwargs...)
end

function inject!(frame::AbstractMatrix{T}, psf::PSFType, pos; A=one(T)) where T
    ctr = reverse(center(frame))
    location = parse_position(pos, ctr)
    tform = Translation(center(psf) - location)
    tform = ImageTransformations.try_static(tform, frame)
    for idx in CartesianIndices(frame)
        I = tform(SVector(idx.I))
        frame[idx] += A * psf(Tuple(I)...)
    end

    return frame
end

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

"""
    inject(cube, [angles], psf, position; A=1, degree=Linear(), fill=0.0)

Injects the matrix `psf` into each frame of `cube` at the given position. If `position` is a Tuple, it will be parsed as pixel coordinates `(x, y)`. If `position` is a `Polar` (from CoordinateTransformations.jl) it will be parsed as polar coordinates centered on the frame center.

If `angles` are provided, each position will be rotated by the give parallactic angles in degrees.

The amplitude of the injected frame can be changed with `A`, and in the case that `psf` is dense matrix, an extrapolator will be created with the given `degree` and `fill` value.
"""
inject(cube::AbstractArray{T,3}, args...; kwargs...) where T =
    inject!(copy(cube), args...; kwargs...)

"""
    inject!(cube, [angles], psf, position; A=1, degree=Linear(), fill=0.0)

In-place version of [`inject`](@ref) which modifies `cube`.
"""
function inject!(cube::AbstractArray{T,3}, kernel, pos; kwargs...) where T
    # All zeros position angles is actually 90° parallactic angle
    angles = Zeros(size(cube, 1))
    return inject!(cube, angles, kernel, pos; kwargs...)
end

function inject!(cube::AbstractArray{T,3}, angles::AbstractVector, kernel, pos; A = one(T), kwargs...) where T
    gen = CubeGenerator(cube, angles, kernel; kwargs...)
    return inject!(cube, gen, pos; A=A)
end

function inject!(cube::AbstractArray{T,3}, gen::CubeGenerator, pos; A=one(T)) where T
    return gen(cube, pos; A=A)
end

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

function normalize_psf!(psf::AbstractMatrix, fwhm, factor=1)
    ap = CircularAperture(center(psf), factor * fwhm / 2)
    area = photometry(ap, psf).aperture_sum
    return psf ./= area
end
normalize_psf(psf, fwhm, factor=1) = normalize_psf!(deepcopy(psf), fwhm, factor)


text:
module day25

using Base.Iterators

ex = """5764801
17807724"""

input = """10604480
4126658"""

parse_input(input) = parse.(BigInt, split(input, "\n"))

function tf(sn, ls)
    val = 1
    for i = 1:ls
        val *= sn
        val = rem(val, 20201227)
    end
    return val
end
@assert tf(7, 8) == 5764801
@assert tf(7, 11) == 17807724

function break_encryption(input)
    pub_keys = parse_input(input)

    # work out card's loop size
    lsc = find_ls(pub_keys[1])

    # @assert tf(pub_keys[1], lss[2]) == tf(pub_keys[2], lss[1])
    encryption_key = tf(pub_keys[2], lsc)
    return encryption_key
end

function find_ls(pub_key)
    val = 1
    for ls in countfrom(1)
        val *= 7
        val %= 20201227
        val == pub_key && return ls
    end
end

part1(input) = break_encryption(input)
@assert part1(ex) == 14897079
@info "Day 25, Part 1 answer: $(part1(input))"

end


text:
"""
TestDataSetsNumbering(DF::DataFrame;PERCENTAGE::Float64=0.6)
    テストデータ作成用行番号と検証用データ作成用行番号を作成するためのサンプリングを行なう

## 引数
* `DF`:サンプリング対象のデータセット
* `PERCENTAGE`:サンプリング割合（デフォルト：0.6）

## 戻り値
* 学習用データセット作成用行番号と検証用データセット作成用行番号
"""

#
function TestDataSetsNumbering(DF::DataFrame;PERCENTAGE::Float64=0.6)
  #元データの行数を求める
  FullNumbering=size(DF,1)

  #学習用データの行番号を求める
  #Num(100分率)抽出するため、抽出対象行を作成
  TrainingNumbering=sample(1:size(DF, 1), Int(round(PERCENTAGE * size(DF, 1),RoundNearest)) ,replace=false, ordered=true)

  #評価用データの行番号を求める
  TestNumbering=setdiff(1:FullNumbering,TrainingNumbering)

  return TrainingNumbering,TestNumbering
end

"""
CreateTestDataSets(DF::DataFrame;PERCENTAGE::Float64=0.6)
    テストデータと検証用データを作成する

## 引数
* `DF`:サンプリング対象のデータセット
* `PERCENTAGE`:サンプリング割合（デフォルト：0.6）

## 戻り値
* 学習用データセットと検証用データセット
"""

#
function CreateTestDataSets(DF::DataFrame;PERCENTAGE::Float64=0.6)
  #元データの行数を求める
  TrainingNumbering,TestNumbering=TestDataSetsNumbering(DF,PERCENTAGE=PERCENTAGE)

  #評価用データの行番号を求める
  Train=DF[TrainingNumbering,:]
  Test=DF[TestNumbering,:]

  return Train,Test
end


text:
#
println("Fd) Symbolic simplification of spherical tensor operators, matrix elements and spherical amplitudes.")
#
using SymEngine
j1  = Basic(:j1);    j2  = Basic(:j2);    j3  = Basic(:j3);    j4  = Basic(:j4);    j5  = Basic(:j5)
m1  = Basic(:m1);    m2  = Basic(:m2);    m3  = Basic(:m3);    m4  = Basic(:m4);    m5  = Basic(:m5)
ja  = Basic(:ja);    jb  = Basic(:jb);    jc  = Basic(:jc);    jd  = Basic(:jd);    je  = Basic(:je)
ma  = Basic(:ma);    mb  = Basic(:mb);    mc  = Basic(:mc);    md  = Basic(:md);    me  = Basic(:me)
k   = Basic(:k);     q   = Basic(:q );    K   = Basic(:K);     Q   = Basic(:Q);
kk  = Basic(:kk);    qq  = Basic(:qq);    KK  = Basic(:KK);    QQ  = Basic(:QQ);

if  true
# Define simple electronic and field operators
Ck      = SphericalTensor.CkOperator(k);                 println("$Ck");         @show Ck
Coulomb = SphericalTensor.CoulombOperator();             println("$Coulomb");    @show Coulomb
Dop     = SphericalTensor.DipoleOperator(false);         println("$Dop");        @show Dop
Tk      = SphericalTensor.TkOperator(k);                 println("$Tk");         @show Tk

gConst  = SphericalTensor.ScalarConstant(false, Basic(:g));                 println("$gConst");        @show gConst
uVector = SphericalTensor.UVector(false);                                   println("$uVector");       @show uVector

# Define spherical tensor product
prod1   = SphericalTensor.TensorProduct(false, j2, Ck, Tk);                 println("$prod1");         @show prod1
prod2   = SphericalTensor.TensorProduct(false, Basic(0), uVector, Dop);     println("$prod2");         @show prod2
prod3   = SphericalTensor.TensorProduct(false, K, prod1, prod2);            println("$prod3");         @show prod3
prod4   = SphericalTensor.TensorProduct(false, K, uVector, Tk);             println("$prod4");         @show prod4

# Define spherical and reduced states
jma     = SphericalTensor.SphericalState(false, ja, ma);                    println("$jma");           @show jma
jmb     = SphericalTensor.SphericalState(false, jb, mb);                    println("$jmb");           @show jmb
j2      = SphericalTensor.ReducedState(j2);                                 println("$j2");            @show j2
j3      = SphericalTensor.ReducedState(j3);                                 println("$j3");            @show j3

# Define spherical matrix elements
me1     = SphericalTensor.SphericalMatrixElement(jma, Tk, q, jmb);          println("$me1");           @show me1
me2     = SphericalTensor.SphericalMatrixElement(jma, prod2, q, jmb);       println("$me2");           @show me2
me3     = SphericalTensor.SphericalMatrixElement(jma, prod3, q, jmb);       println("$me3");           @show me3
me4     = SphericalTensor.SphericalMatrixElement(false, jma, gConst, prod3, q, jmb);
                                                                            println("$me4");           @show me4
me2s    = SphericalTensor.SphericalMatrixElement(true,  jma, gConst, prod2, q, jmb);
                                                                            println("$me2s");          @show me2s
me5     = SphericalTensor.SphericalMatrixElement(false, jma, gConst, prod4, q, jmb);
                                                                            println("$me5");           @show me5
me5s    = SphericalTensor.SphericalMatrixElement(true,  jma, gConst, prod4, q, jmb);
                                                                            println("$me5s");          @show me5s

# Expand spherical operators
we1     = SphericalTensor.expandSphericalTensorComponent(Coulomb, q);       println("$Coulomb");       @show we1
we2     = SphericalTensor.expandSphericalTensorComponent(uVector, q);       println("$uVector");       @show we2
we3     = SphericalTensor.expandSphericalTensorComponent(prod1, q);         println("$prod1");         @show we3
we4     = SphericalTensor.expandSphericalTensorComponent(prod2, q);         println("$prod2");         @show we4
we5     = SphericalTensor.expandSphericalTensorComponent(prod3, q);         println("$prod3");         @show we5

# Expand spherical operators
wm1     = SphericalTensor.expandSphericalMatrixElements([me1]);             println("\n $me1");        @show wm1
wm2     = SphericalTensor.expandSphericalMatrixElements([me2]);             println("\n $me2");        @show wm2
wm3     = SphericalTensor.expandSphericalMatrixElements([me3]);             println("\n $me3");        @show wm3
wm4     = SphericalTensor.expandSphericalMatrixElements([me4]);             println("\n $me4");        @show wm4
wm5     = SphericalTensor.expandSphericalMatrixElements([me2s, me2]);       println("\n $me2s $me2");  @show wm5
wm6     = SphericalTensor.expandSphericalMatrixElements([me5s, me5]);       println("\n $me5s $me5");  @show wm6

elseif  false

end


text:
#
# Multiplicative updating algorithm
#
# Reference:
#   Daniel D. Lee and H. Sebastian Seung. Algorithms for Non-negative
#   Matrix Factorization. Advances in NIPS, 2001.
#

mutable struct MultUpdate{T}
    obj::Symbol         # objective :mse or :div
    maxiter::Int        # maximum number of iterations
    verbose::Bool       # whether to show procedural information
    tol::T              # change tolerance upon convergence
    lambda::T           # regularization coefficient

    function MultUpdate{T}(;obj::Symbol=:mse,
                            maxiter::Integer=100,
                            verbose::Bool=false,
                            tol::Real=cbrt(eps(T)),
                            lambda::Real=sqrt(eps(T))) where T

        obj == :mse || obj == :div || throw(ArgumentError("Invalid value for obj."))
        maxiter > 1 || throw(ArgumentError("maxiter must be greater than 1."))
        tol > 0 || throw(ArgumentError("tol must be positive."))
        lambda >= 0 || throw(ArgumentError("lambda must be non-negative."))

        new{T}(obj, maxiter, verbose, tol, lambda)
    end
end

function solve!(alg::MultUpdate, X, W, H, nmfinit, nmfalg)

    if alg.obj == :mse
        nmf_skeleton!(MultUpdMSE(alg.lambda), X, W, H, alg.maxiter, alg.verbose, alg.tol, nmfinit, nmfalg)
    else # alg == :div
        nmf_skeleton!(MultUpdDiv(alg.lambda), X, W, H, alg.maxiter, alg.verbose, alg.tol, nmfinit, nmfalg)
    end
end

# the multiplicative updating algorithm for MSE objective

struct MultUpdMSE{T} <: NMFUpdater{T}
    lambda::T
end

struct MultUpdMSE_State{T}
    WH::Matrix{T}
    WtX::Matrix{T}
    WtWH::Matrix{T}
    XHt::Matrix{T}
    WHHt::Matrix{T}

    function MultUpdMSE_State{T}(X, W::Matrix{T}, H::Matrix{T}) where T
        p, n, k = nmf_checksize(X, W, H)
        new{T}(W * H,
               Matrix{T}(undef, k, n),
               Matrix{T}(undef, k, n),
               Matrix{T}(undef, p, k),
               Matrix{T}(undef, p, k))
    end
end

prepare_state(::MultUpdMSE{T}, X, W, H) where {T} = MultUpdMSE_State{T}(X, W, H)
evaluate_objv(::MultUpdMSE, s::MultUpdMSE_State, X, W, H) = sqL2dist(X, s.WH)

function update_wh!(upd::MultUpdMSE{T}, s::MultUpdMSE_State{T}, X, W::Matrix{T}, H::Matrix{T}) where T

    # fields
    lambda = upd.lambda
    WH = s.WH
    WtX = s.WtX
    WtWH = s.WtWH
    XHt = s.XHt
    WHHt = s.WHHt

    # update H
    Wt = transpose(W)
    mul!(WtX, Wt, X)
    mul!(WtWH, Wt, WH)

    @inbounds for i = 1:length(H)
        H[i] *= (WtX[i] / (WtWH[i] + lambda))
    end
    mul!(WH, W, H)

    # update W
    Ht = transpose(H)
    mul!(XHt, X, Ht)
    mul!(WHHt, WH, Ht)

    @inbounds for i = 1:length(W)
        W[i] *= (XHt[i] / (WHHt[i] + lambda))
    end
    mul!(WH, W, H)
end

# the multiplicative updating algorithm for divergence objective

struct MultUpdDiv{T} <: NMFUpdater{T}
    lambda::T
end

struct MultUpdDiv_State{T}
    WH::Matrix{T}
    sW::Matrix{T}     # sum(W, 1)
    sH::Matrix{T}     # sum(H, 2)
    Q::Matrix{T}      # X ./ (WH + lambda): size (p, n)
    WtQ::Matrix{T}    # W' * Q: size (k, n)
    QHt::Matrix{T}    # Q * H': size (p, k)

    function MultUpdDiv_State{T}(X, W::Matrix{T}, H::Matrix{T}) where T
        p, n, k = nmf_checksize(X, W, H)
        new{T}(W * H,
               Matrix{T}(undef, 1, k),
               Matrix{T}(undef, k, 1),
               Matrix{T}(undef, p, n),
               Matrix{T}(undef, k, n),
               Matrix{T}(undef, p, k))
    end
end

prepare_state(::MultUpdDiv{T}, X, W, H) where {T} = MultUpdDiv_State{T}(X, W, H)
evaluate_objv(::MultUpdDiv, s::MultUpdDiv_State, X, W, H) = gkldiv(X, s.WH)

function update_wh!(upd::MultUpdDiv{T}, s::MultUpdDiv_State{T}, X, W::Matrix{T}, H::Matrix{T}) where T

    p = size(X, 1)
    n = size(X, 2)
    k = size(W, 2)
    pn = p * n

    # fields
    lambda = upd.lambda
    sW = s.sW
    sH = s.sH
    WH = s.WH
    Q = s.Q
    WtQ = s.WtQ
    QHt = s.QHt

    @assert size(Q) == size(X)

    # update H
    @inbounds for i = 1:length(X)
        Q[i] = X[i] / (WH[i] + lambda)
    end
    mul!(WtQ, transpose(W), Q)
    sum!(fill!(sW, 0), W)
    @inbounds for j = 1:n, i = 1:k
        H[i,j] *= (WtQ[i,j] / sW[i])
    end
    mul!(WH, W, H)

    # update W
    @inbounds for i = 1:length(X)
        Q[i] = X[i] / (WH[i] + lambda)
    end
    mul!(QHt, Q, transpose(H))
    sum!(fill!(sH, 0), H)
    @inbounds for j = 1:k, i = 1:p
        W[i,j] *= (QHt[i,j] / sH[j])
    end
    mul!(WH, W, H)
end


text:
# Copyright (c) 2018: Matthew Wilhelm & Matthew Stuber.
# This code is licensed under MIT license (see LICENSE.md for full details)
#############################################################################
# McCormick.jl
# A McCormick operator library in Julia
# See https://github.com/PSORLab/McCormick.jl
#############################################################################
# src/McCormick.jl
# Main file for McCormick.jl. Contains constructors and helper functions.
#############################################################################

__precompile__()

module McCormick

using DocStringExtensions, LinearAlgebra
using DiffRules: diffrule
using StaticArrays: @SVector, SVector, zeros, ones
using ForwardDiff: Dual, Partials
using NaNMath, SpecialFunctions

import Base: +, -, *, /, convert, in, isempty, one, zero, real, eps, max, min,
             abs, inv, exp, exp2, exp10, expm1, log, log2, log10, log1p, acosh,
             sech, csch, coth, acsch, acoth, asech, rad2deg, deg2rad, sqrt, sin,
             cos, tan, min, max, sec, csc, cot, ^, step, sign, intersect,
             promote_rule, asinh, atanh, tanh, atan, asin, cosh, acos,
             sind, cosd, tand, asind, acosd, atand,
             secd, cscd, cotd, asecd, acscd, acotd, isone, isnan, empty,
             <, <=, ==, fma, cbrt, sinpi, cospi, union

using IntervalArithmetic
using IntervalArithmetic: @round
if isdefined(IntervalArithmetic, :big53)
    big_val(x) = IntervalArithmetic.big53(x)
else
    big_val(x) = IntervalArithmetic.bigequiv(x)
end

using IntervalRootFinding
import IntervalArithmetic: dist, mid, pow, +, -, *, /, convert, in, isempty,
                           one, zero, real, eps, max, min, abs, exp,
                           expm1, log, log2, log10, log1p, sqrt, ^,
                           sin, cos, tan, min, max, sec, csc, cot, step, sech,
                           csch, coth, acsch, acoth, asech,
                           sign, dist, mid, pow, Interval, interval, sinh, cosh,
                           ∩, IntervalBox, bisect, isdisjoint, length,
                           atan, asin, acos, AbstractInterval, atomic,
                           sind, cosd, tand, asind, acosd, atand,
                           secd, cscd, cotd, asecd, acscd, acotd, half_pi,
                           setrounding, diam, isthin, abs2

if ~(VERSION < v"1.1-")
    import SpecialFunctions: erf, erfc, erfinv, erfcinv
    export erf, erfinv, erfc, erfcinv, erf_kernel,
           erfinv_kernel, erfc_kernel, erfcinv_kernel
end

import Base.MathConstants.golden

# Export forward operators
export MC, cc, cv, Intv, lo, hi,  cc_grad, cv_grad, cnst, +, -, *, /, convert,
       one, zero, dist, real, eps, mid, exp, exp2, exp10, expm1, log, log2,
       log10, log1p, acosh, sqrt, sin, cos, tan, min, max, sec, csc, cot, ^,
       abs, step, sign, pow, in, isempty, intersect, length, mid3,
       acos, asin, atan, sinh, cosh, tanh, asinh, atanh, inv, sqr, sech,
       csch, coth, acsch, acoth, asech, rad2deg, deg2rad, diam,
       sind, cosd, tand, asind, acosd, atand, nan,
       sinhd, coshd, tanhd, asinhd, acoshd, atanhd,
       secd, cscd, cotd, asecd, acscd, acotd,
       secdh, cschd, cothd, asechd, acschd, acothd, isone, isnan, interval_MC,
       relu, param_relu, leaky_relu, maxsig, maxtanh, softplus, pentanh,
       sigmoid, bisigmoid, softsign, gelu, elu, selu, swish1,
       positive, negative, lower_bnd, upper_bnd, bnd, xlogx,
       <, <=, ==, fma, cbrt, abs2, sinpi, cospi, arh, xexpax, trilinear,
       xabsx, logcosh

# Export kernel operators
export plus_kernel, minus_kernel, mult_kernel, div_kernel, max_kernel,
       min_kernel, log_kernel, log2_kernel, log10_kernel, log1p_kernel,
       acosh_kernel, sqrt_kernel, exp_kernel, exp2_kernel, exp10_kernel,
       expm1_kernel, div_kernel, cos_kernel, sin_kernel, sinh_kernel,
       tanh_kernel, asinh_kernel, atanh_kernel, tan_kernel, acos_kernel,
       asin_kernel, atan_kernel, cosh_kernel, deg2rad_kernel, rad2deg_kernel,
       sec_kernel, csc_kernel, cot_kernel, asec_kernel, acsc_kernel,
       acot_kernel, sech_kernel, csch_kernel, coth_kernel, acsch_kernel,
       acoth_kernel, sind_kernel, cosd_kernel, tand_kernel, secd_kernel,
       cscd_kernel, cotd_kernel, asind_kernel, acosd_kernel,
       atand_kernel, asecd_kernel, acscd_kernel, acotd_kernel, relu_kernel,
       param_relu_kernel, leaky_relu_kernel, maxsig_kernel,
       maxtanh_kernel, softplus_kernel, pentanh_kernel, sigmoid_kernel,
       bisigmoid_kernel, softsign_kernel, gelu_kernel, elu_kernel, selu_kernel,
       swish1_kernel, positive_kernel, negative_kernel, lower_bnd_kernel,
       upper_bnd_kernel, bnd_kernel, xlogx_kernel, fma_kernel, cbrt_kernel,
       abs2_kernel, sinpi_kernel, cospi_kernel, arh_kernel, xexpax_kernel,
       trilinear_kernel, xabsx_kernel, logcosh_kernel

export seed_gradient, RelaxTag, NS, MV, Diff

export MCCallback, gen_expansion_params!, implicit_relax_h!, DenseMidInv,
       NewtonGS, KrawczykCW

"""
$(TYPEDEF)

An `abstract type` the subtypes of which define the manner of relaxation that will
be performed for each operator applied to the MC object. Currently, the `struct NS`
which specifies that standard (Mitsos 2009) are to be used is fully supported. Limited
support is provided for differentiable McCormick relaxations specified by `struct Diff`
(Khan 2017) and struct MV `struct MV` (Tsoukalas 2011.) A rounding-safe implementation of
the standard McCormick relaxations is specified by the struct `NSSafe` which is work in
progress.
"""
abstract type RelaxTag end

struct NS <: RelaxTag end
struct MV <: RelaxTag end
struct Diff <: RelaxTag end

const MC_ENV_MAX_INT = 100
const MC_ENV_TOL = 1E-10

const MC_INTERSECT_NOOP_FALLBACK = true
const MC_INTERSECT_TOL = 1E-13

const MC_DIFF_MU = 1
const MC_DIFF_MUT = convert(Float64, MC_DIFF_MU)
const MC_DIFF_MU1 = MC_DIFF_MU + 1
const MC_DIFF_MU1T = convert(Float64, MC_DIFF_MU1)
const MC_DIFF_MU1N = MC_DIFF_MU - 1
const MC_DIFF_DIV = MC_DIFF_MU1^(-1/MC_DIFF_MU)

const MC_MV_TOL = 1E-8
const MC_DEGEN_TOL = 1E-14
const MC_DOMAIN_TOL = 1E-10

const IntervalConstr = interval
const Half64 = Float64(0.5)
const Two64 = Float64(2.0)
const Three64 = Float64(3.0)
const EqualityTolerance = Float64(1E-12)
const DegToRadIntv = atomic(Interval{Float64}, pi)/Interval(180.0)
const one_intv = one(Interval{Float64})
const half_intv = Interval{Float64}(0.5)
const two_intv = Interval{Float64}(2.0)
const log2_intv = log(Interval{Float64}(2.0))
const log10_intv = log(Interval{Float64}(10.0))

const NumberNotRelax = Union{Bool, Float16, Float32, Signed, Unsigned, BigFloat,
                             Int8, Int64, Int32, Int16, Int128}

########### number list
int_list = [Int8,UInt8,Int16,UInt16,
            Int32,UInt32,Int64,UInt64,Int128,UInt128]
float_list = [Float16,Float32,Float64]

########### differentiable functions unitary functions
CVList = Symbol[:cosh,:exp,:exp2,:exp10]        ### function is convex
CCList = Symbol[:acosh,:log,:log2,:log10,:sqrt] ### function is concave
CCtoCVList = Symbol[:asin,:sinh,:atanh,:tan]    ### function is concave then convex
CVtoCCList = Symbol[:atan,:acos,:tanh,:asinh]   ### function is convex then concave
Template_List = union(CVList,CCList,CCtoCVList,CVtoCCList)

########### non differentiable and non-unitary functions
OtherList = Symbol[:sin,:cos,:min,:max,:abs,:step, :sign, :inv, :*, :+, :-, :/,
:promote_rule, :convert, :one, :zero, :real, :dist, :eps, :fma, :^]

NaNsqrt(x) = x < 0.0 ? NaN : Base.sqrt(x)

"""
$(TYPEDSIGNATURES)

Creates a `x::SVector{N,Float64}` object that is one at `x[j]` and zero everywhere else.
"""
function seed_gradient(j::Int64, x::Val{N}) where N
  return SVector{N,Float64}(ntuple(i -> i == j ? 1.0 : 0.0, Val{N}()))
end

"""
$(TYPEDSIGNATURES)

Calculates the middle of three numbers returning the value and the index where `x >= y`.
"""
function mid3(x::Float64, y::Float64, z::Float64)
  (((x >= y) && (y >= z)) || ((z >= y) && (y >= x))) && (return y, 2)
  (((y >= x) && (x >= z)) || ((z >= x) && (x >= y))) && (return x, 1)
  return z, 3
end

"""
$(TYPEDSIGNATURES)

Calculates the middle of three numbers (x,y,z) returning the value where `x <= y`.
"""
function mid3v(x::Float64, y::Float64, z::Float64)
    z <= x && (return x)
    y <= z && (return y)
    return z
end

"""
$(TYPEDSIGNATURES)

Takes the concave relaxation gradient 'cc_grad', the convex relaxation gradient
'cv_grad', and the index of the midpoint returned 'id' and outputs the appropriate
gradient according to McCormick relaxation rules.
"""
function mid_grad(cc_grad::SVector{N,Float64}, cv_grad::SVector{N,Float64}, id::Int64) where N
    if (id == 1)
        return cc_grad
    elseif (id == 2)
        return cv_grad
    end
    return zero(SVector{N,Float64})
end

"""
$(TYPEDSIGNATURES)

Calculates the value of the slope line segment between `(xL, f(xL))` and `(xU, f(xU))`
defaults to evaluating the derivative of the function if the interval is tight.
"""
@inline function dline_seg(f::Function, df::Function, x::Float64, xL::Float64, xU::Float64)
    delta = xU - xL
    if delta == 0.0
        return f(x), df(x)
    else
        yL = f(xL)
        yU = f(xU)
        return (yL*(xU - x) + yU*(x - xL))/delta, (yU - yL)/delta
    end
end
@inline function dline_seg(f::Function, df::Function, x::Float64, xL::Float64, xU::Float64, n::Int64)
    delta = xU - xL
    if delta == 0.0
        return f(x, n), df(x, n)
    else
        yL = f(xL, n)
        yU = f(xU, n)
        return (yL*(xU - x) + yU*(x - xL))/delta, (yU - yL)/delta
    end
end
@inline function dline_seg(f::Function, df::Function, x::Float64, xL::Float64, xU::Float64, c::Float64)
    delta = xU - xL
    if delta == 0.0
        return f(x, c), df(x, c)
    else
        yL = f(xL, c)
        yU = f(xU, c)
        return (yL*(xU - x) + yU*(x - xL))/delta, (yU - yL)/delta
    end
end

"""
$(TYPEDSIGNATURES)

Refines convex/concave relaxations `cv` and `cc` with associated subgradients
`cv_grad` and `cc_grad` by intersecting them with the interval boudns `xL`
and `xU`.
"""
function cut(xL::Float64, xU::Float64, cv::Float64, cc::Float64,
             cv_grad::SVector{N,Float64}, cc_grad::SVector{N,Float64}) where N

    if cc > xU
        cco = xU
        cc_grado = zero(SVector{N,Float64})
    else
        cco = cc
        cc_grado = cc_grad
    end
    if cv < xL
        cvo = xL
        cv_grado = zero(SVector{N,Float64})
    else
        cvo = cv
        cv_grado = cv_grad
    end
    return cvo, cco, cv_grado, cc_grado
end

lo(x::Interval{Float64}) = x.lo
hi(x::Interval{Float64}) = x.hi

function step(x::Interval{Float64})
     isempty(x) && return emptyinterval(x)
     xmin::Float64 = ((x.lo) < 0.0) ? 0.0 : 1.0
     xmax::Float64 = ((x.hi) >= 0.0) ? 1.0 : 0.0
     return Interval{Float64}(xmin,xmax)
end

#abstract type AbstractMC <: Real end
"""
$(TYPEDEF)

`MC{N, T <: RelaxTag} <: Real` is the McCormick (w/ (sub)gradient) structure which is used to overload
standard calculations. The fields are:
$(TYPEDFIELDS)
"""
struct MC{N, T <: RelaxTag} <: Real
    "Convex relaxation"
    cv::Float64
    "Concave relaxation"
    cc::Float64
    "Interval bounds"
    Intv::Interval{Float64}
    "(Sub)gradient of convex relaxation"
    cv_grad::SVector{N,Float64}
    "(Sub)gradient of concave relaxation"
    cc_grad::SVector{N,Float64}
    "Boolean indicating whether the relaxations are constant over the domain. True if bounding an interval/constant.
     False, otherwise. This may change over the course of a calculation `cnst` for `zero(x)` is `true` even if `x.cnst`
     is `false`."
    cnst::Bool
    function MC{N,T}(cv1::Float64, cc1::Float64, Intv1::Interval{Float64},
                     cv_grad1::SVector{N,Float64}, cc_grad1::SVector{N,Float64},
                     cnst1::Bool) where {N, T <: RelaxTag}
        new(cv1, cc1, Intv1, cv_grad1, cc_grad1, cnst1)
    end
end

"""
MC{N,T}(y::Interval{Float64})

Constructs McCormick relaxation with convex relaxation equal to `y.lo` and
concave relaxation equal to `y.hi`.
"""
function MC{N,T}(y::Interval{Float64}) where {N, T <: RelaxTag}
    MC{N,T}(y.lo, y.hi, y, zero(SVector{N,Float64}),
                           zero(SVector{N,Float64}), true)
end

"""
MC{N,T}(y::Float64)

Constructs McCormick relaxation with convex relaxation equal to `y` and
concave relaxation equal to `y`.
"""
MC{N,T}(y::Float64) where {N, T <: RelaxTag} = MC{N,T}(Interval{Float64}(y))
function MC{N,T}(y::Y) where {N, T <: RelaxTag, Y <: AbstractIrrational}
    MC{N,T}(Interval{Float64}(y))
end
MC{N,T}(y::Q) where {N, T <: RelaxTag, Q <: NumberNotRelax} = MC{N,T}(Interval{Float64}(y))

"""
MC{N,T}(cv::Float64, cc::Float64)

Constructs McCormick relaxation with convex relaxation equal to `cv` and
concave relaxation equal to `cc`.
"""
function MC{N,T}(cv::Float64, cc::Float64) where {N, T <: RelaxTag}
    MC{N,T}(cv, cc, Interval{Float64}(cv, cc), zero(SVector{N,Float64}),
            zero(SVector{N,Float64}), true)
end

"""
MC{N,T}(val::Float64, Intv::Interval{Float64}, i::Int64)

Constructs McCormick relaxation with convex relaxation equal to `val`,
concave relaxation equal to `val`, interval bounds of `Intv`, and a unit subgradient
with nonzero's ith dimension of length N.
"""
function MC{N,T}(val::Float64, Intv::Interval{Float64}, i::Int64) where {N, T <: RelaxTag}
    MC{N,T}(val, val, Intv, seed_gradient(i, Val{N}()), seed_gradient(i, Val{N}()), false)
end
function MC{N,T}(x::MC{N,T}) where {N, T <: RelaxTag}
    MC{N,T}(x.cv, x.cc, x.Intv, x.cv_grad, x.cc_grad, x.cnst)
end

Intv(x::MC) = x.Intv
lo(x::MC) = x.Intv.lo
hi(x::MC) = x.Intv.hi
cc(x::MC) = x.cc
cv(x::MC) = x.cv
cc_grad(x::MC) = x.cc_grad
cv_grad(x::MC) = x.cv_grad
cnst(x::MC) = x.cnst
length(x::MC) = length(x.cc_grad)

diam(x::MC) = diam(x.Intv)
isthin(x::MC) = isthin(x.Intv)

function isone(x::MC)
    flag = true
    flag &= (x.Intv.lo == 1.0)
    flag &= (x.Intv.hi == 1.0)
    flag &= x.cnst
    return flag
end

"""
$(TYPEDSIGNATURES)

Defines a local 1D newton method to solve for the root of `f` between the bounds
`xL` and `xU` using `x0` as a starting point. The derivative of `f` is `df`. The
inputs `envp1` and `envp2` are the envelope calculation parameters.
"""
function newton(x0::Float64, xL::Float64, xU::Float64, f::Function, df::Function,
               envp1::Float64, envp2::Float64)

    dfk = 0.0
    xk = max(xL, min(x0, xU))
    fk::Float64 = f(xk, envp1, envp2)

    for i = 1:MC_ENV_MAX_INT
        dfk = df(xk, envp1, envp2)
        if (abs(fk) < MC_ENV_TOL)
            return (xk, false)
        end
        (dfk == 0.0) && return (0.0, true)
        if (xk == xL && fk/dfk > 0.0)
            return (xk, false)
        elseif (xk == xU && fk/dfk < 0.0)
            return (xk, false)
        end
        xk = max(xL, min(xU, xk - fk/dfk))
        fk = f(xk, envp1, envp2)
    end

    (0.0, true)
end

"""
$(TYPEDSIGNATURES)

Defines a local 1D secant method to solve for the root of `f` between
the bounds `xL` and `xU` using `x0` and `x1` as a starting points. The inputs
`envp1` and `envp2` are the envelope calculation parameters.
"""
function secant(x0::Float64, x1::Float64, xL::Float64, xU::Float64, f::Function,
                envp1::Float64, envp2::Float64)

    xkm = max(xL, min(xU, x0))
    xk = max(xL, min(xU, x1))
    fkm::Float64 = f(xkm, envp1, envp2)

    for i = 1:MC_ENV_MAX_INT
        fk::Float64  = f(xk, envp1, envp2)
        Bk = (fk - fkm)/(xk - xkm)
        if (abs(fk) < MC_ENV_TOL)
            return (xk, false)
        end
        (Bk == 0.0) && return (0.0, true)
        if (xk == xL) && (fk/Bk > 0.0)
            return (xk, false)
        elseif (xk == xU) && (fk/Bk < 0.0)
            return (xk, false)
        end
        xkm = xk
        fkm = fk
        xk = max(xL, min(xU, xk - fk/Bk))
    end

    (0.0,true)
end

"""
$(TYPEDSIGNATURES)

Defines a local 1D golden section method to solve for the root of `f` between
the bounds `xL` and `xU` using `x0` as a starting point. Define iteration used
in golden section method. The inputs `envp1` and `envp2` are the envelope
calculation parameters.
"""
function golden_section(xL::Float64, xU::Float64, f::Function, envp1::Float64,
                        envp2::Float64)
  fL::Float64 = f(xL, envp1, envp2)
  fU::Float64 = f(xU, envp1, envp2)

  fL*fU > 0.0 && (return NaN)
  xm = xU - (2.0 - golden)*(xU - xL)
  fm::Float64 = f(xm,envp1,envp2)
  return golden_section_it(1, xL, fL, xm, fm, xU, fU, f, envp1, envp2)
end

"""
$(TYPEDSIGNATURES)

Define iteration used in golden section method. The inputs `fa`,`fb`, and `fc`,
are the function `f` evaluated at `a`,`b`, and `c` respectively. The inputs
`envp1` and `envp2` are the envelope calculation parameters. The value `init` is
the iteration number of the golden section method.
"""
function golden_section_it(init::Int, a::Float64, fa::Float64, b::Float64,
                           fb::Float64, c::Float64, fc::Float64, f::Function,
                           envp1::Float64, envp2::Float64)
    flag = (c - b > b - a)
    x = flag ? b + (2.0 - golden)*(c - b) : b - (2.0 - golden)*(b - a)
    itr = init
    if abs(c-a) < MC_ENV_TOL*(abs(b) + abs(x)) || (itr > MC_ENV_MAX_INT)
        return (c + a)/2.0
    end
    itr += 1
    fx::Float64 = f(x, envp1, envp2)
    if flag
        if fa*fx < 0.0
            out = golden_section_it(itr, a, fa, b, fb, x, fx, f, envp1, envp2)
        else
            out = golden_section_it(itr, b, fb, x, fx, c, fc, f, envp1, envp2)
        end
    else
        if fa*fb < 0.0
            out = golden_section_it(itr, a, fa, x, fx, b, fb, f, envp1, envp2)
        else
            out = golden_section_it(itr, x, fx, b, fb, c, fc, f, envp1, envp2)
        end
    end
    return out
end

function check_relaxation_error!(x::MC{N,T}) where {N, T<:RelaxTag}
    if x.cc < x.cv
        @error "After performing calculations, a convex relaxation greater
                than the associated concave relaxation was encountered. Please
                open an issue at https://github.com/PSORLab/McCormick.jl/issues/new
                so we may address this error in the future."
    end
    return nothing
end

include("forward_operators/forward.jl")
include("implicit_routines/implicit.jl")

using Reexport
@reexport using IntervalArithmetic
@reexport using StaticArrays
@reexport using SpecialFunctions

end


text:
include("../../module/CustomStencil.jl")
include("../../misc/misc.jl")
include("../../misc/cpu_stencils.jl")

## Star Stencil definition with radius = 4
coefs = [0.75;0.5]
star_stencil = @def_stencil_expression c[0]D[x,y,z] + @sum(i, 1,1, c[i]*(
            D[x+i,y,z] + D[x,y+i,z] + D[x,y,z+i] +
            D[x-i,y,z] + D[x,y-i,z] + D[x,y,z-i]))
a = 2
star_stencil = def_stencil_epxpression(:(c[0]D[x,y,z] + @sum(i, 1,$(a), c[i]*(
            D[x+i,y,z] + D[x,y+i,z] + D[x,y,z+i] +
            D[x-i,y,z] + D[x,y-i,z] + D[x,y,z-i]))))
st_def = CreateStencilDefinition(star_stencil, coefs)
st_inst1 = NewStencilInstance(st_def, m_step=false)

## Input Data size Definition
radius = 4
nx = 8
ny = 8
nz = 8
dx = 1<<(nx)
dy = 1<<(ny)
dz = 1<<(nz)
data = CreateData(dx,dy,dz)

t_steps = 16
## Compile Functions

gpu_out = ApplyFFTstencil(st_inst1, data, t_steps)
CUDA.cuProfilerStart()
NVTX.@range "r1" begin

     t_steps = 32
    gpu_out = ApplyFFTstencil(st_inst1, data, t_steps,32) end
NVTX.@range "r2" begin
	t_steps = 48
    gpu_out = ApplyFFTstencil(st_inst1, data, t_steps,48) end
NVTX.@range "r3" begin
t_steps = 64
    gpu_out = ApplyFFTstencil(st_inst1, data, t_steps,64) end
NVTX.@range "r4" begin

t_steps = 80
    gpu_out = ApplyFFTstencil(st_inst1, data, t_steps,80) end
NVTX.@range "r5" begin
t_steps = 96
    gpu_out = ApplyFFTstencil(st_inst1, data, t_steps,96) end
exit()


text:
# Julia wrapper for header: /usr/local/include/libavutil/twofish.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

export
    av_twofish_alloc,
    av_twofish_init,
    av_twofish_crypt

function av_twofish_alloc()
    ccall((:av_twofish_alloc,libavutil),Ptr{AVTWOFISH},())
end

function av_twofish_init(ctx,key,key_bits::Integer)
    ccall((:av_twofish_init,libavutil),Cint,(Ptr{AVTWOFISH},Ptr{UInt8},Cint),ctx,key,key_bits)
end

function av_twofish_crypt(ctx,dst,src,count::Integer,iv,decrypt::Integer)
    ccall((:av_twofish_crypt,libavutil),Void,(Ptr{AVTWOFISH},Ptr{UInt8},Ptr{UInt8},Cint,Ptr{UInt8},Cint),ctx,dst,src,count,iv,decrypt)
end

max_stars_repo_name: asafasban/AAVideoIO

max_issues_repo_name: asafasban/AAVideoIO

max_forks_repo_name: asafasban/AAVideoIO


text:
"""
    initsim!(sim;flags=FFTW.MEASURE)

Initialize all arrays, measures and transform libraries for a
particular simulation.
"""
function initsim!(sim;flags=FFTW.MEASURE)
    @unpack L,N = sim
    X,K,dX,dK,DX,DK,T = makearraystransforms(L,N)
    espec = 0.5*k2(L,N)
    @pack! sim = T,X,K,espec
    return nothing
end

"""
    V(x,t) = ...

Define the system potential, with default zero. Should be defined as a scalar
function (without `.`), suitable for broadcasting on spatial arrays via `V.(...)`.
"""
V(x,t) = 0.0
V(x,y,t) = 0.0
V(x,y,z,t) = 0.0

"""
    χ = nlin(ϕ,sim,t)

Evalutes nonlinear terms in `x`-space, returning to `k`-space.
"""
function nlin(ϕ,sim,t)
    @unpack g,x,y = sim
    ψ = xspace(ϕ,sim)
    ψ .*= @. g*abs2(ψ) + V(x,y',t)
    return kspace(ψ,sim)
end

"""
    nlin!(ϕ,sim::Sim{D},t)

Mutating evaluation of position space nonlinear terms.
Dispatches on dimension `D`, using potential `V(x...,t)`.
"""
function nlin!(dϕ,ϕ,sim::Sim{1},t)
    @unpack g,X,V0 = sim; x = X[1]
    dϕ .= ϕ
    xspace!(dϕ,sim)
    @. dϕ *= V0 + V(x,t) + g*abs2(dϕ)
    kspace!(dϕ,sim)
    return nothing
end

function nlin!(dϕ,ϕ,sim::Sim{2},t)
    @unpack g,X,V0 = sim; x,y = X
    y = y'
    dϕ .= ϕ
    xspace!(dϕ,sim)
    @. dϕ *= V0 + V(x,y,t) + g*abs2(dϕ)
    kspace!(dϕ,sim)
    return nothing
end

function nlin!(dϕ,ϕ,sim::Sim{3},t)
    @unpack g,X,V0 = sim; x,y,z = X
    y = y'; z = reshape(z,(1,1,length(z)))
    dϕ .= ϕ
    xspace!(dϕ,sim)
    @. dϕ *= V0 + V(x,y,z,t) + g*abs2(dϕ)
    kspace!(dϕ,sim)
    return nothing
end

"""
    χ = Lgp!(dϕ,ϕ,sim,t)

In-place evaluation of Gross-Pitaevskii equation defined by `nlin!` and parameters in `sim`.
Allows imaginary time `γ`, and evolves in rotating frame defined by chemical potential `μ`.
"""
function Lgp!(dϕ,ϕ,sim,t)
    @unpack γ,μ,espec = sim
    nlin!(dϕ,ϕ,sim,t)
    @. dϕ = -im*(1.0 - im*γ)*(dϕ + (espec - μ)*ϕ)
    return nothing
end

"""
    χ = Lgp(ϕ,sim,t)

Evaluate Gross-Pitaevskii equation defined by `nlin` and parameters in `sim`.
Allows imaginary time `γ`, and evolves in rotating frame defined by chemical potential `μ`.
"""
function Lgp(ϕ,sim,t)
    @unpack μ,γ,espec = sim
    chi = nlin(ϕ,sim,t)
    return @. -im*(1.0 - im *γ)*((espec - μ)*ϕ + chi )
end

function internalnorm(u,t)
    return sum((abs2.(u) .> 1e-6*maximum(abs2.(u))).*abs2.(u))
end

"""
    runsim(sim,ϕ;info,tplot,nfiles)

Call DifferentialEquations to solve Gross-Pitaevskii equation.
"""
function runsim(sim,ϕ=sim.ϕi;info=true,tplot=false,nfiles=false)
    @unpack nfiles,path,filename = sim

    function savefunction(ψ...)
        isdir(path) || mkpath(path)
        i = findfirst(x->x== ψ[2],sim.t)
        padi = lpad(string(i),ndigits(length(sim.t)),"0")
        info && println("⭆ Save $i at t = $(trunc(ψ[2];digits=3))")
        # tofile = path*"/"*filename*padi*".jld2"
        tofile = joinpath(path,filename*padi*".jld2")
        save(tofile,"ψ",ψ[1],"t",ψ[2])
    end

    savecb = FunctionCallingCallback(savefunction;
                     funcat = sim.t, # times to save at
                     func_everystep=false,
                     func_start = true,
                     tdir=1)

    prob = ODEProblem(Lgp!,ϕ,(sim.ti,sim.tf),sim)
    info && @info "⭆ 𝒅𝜳 Evolving in kspace"
    info && @info "⭆ Damping γ = $(sim.γ)"
    (info && nfiles) && @info "⭆ Saving to "*path
    nfiles ?
    (sol = solve(prob,alg=sim.alg,saveat=sim.t[end],reltol=sim.reltol,callback=savecb,dense=false,maxiters=1e10,progress=true)) :
    (sol = solve(prob,alg=sim.alg,saveat=sim.t,reltol=sim.reltol,dense=false,maxiters=1e10,progress=true))
    info && @info "⭆ Finished."
return sol
end

# TODO vortex lattice in 2D, persistent current in 3D examples

function testsim(sim)
    err = false
    sol = try
            runsim(sim;info=false)
        catch e
            err = true
        end
return sol,err
end


text:
#ДАНО: Робот - в юго-западном углу поля, на котором расставлено некоторое количество маркеров.
#РЕЗУЛЬТАТ: Функция вернула значение средней температуры всех замаркированных клеток.
tempsum = 0
numofcells = 0
function count_avg_temp(r)
    dir = Ost
    move_while_can!(r,dir)
    dir = inverse(dir)
    while isborder(r,Nord)==0
        move!(r,Nord)
        move_while_can!(r,dir)
        dir = inverse(dir)
    end
    if numofcells > 0
        return tempsum/numofcells
    else
        return 0
    end
end
function move_while_can!(r,side)
    checkcell(r)
    while isborder(r,side)==0
        move!(r,side)
        checkcell(r)
    end
end
function checkcell(r::Robot)
   global tempsum, numofcells
   if ismarker(r)==1
        tempsum += temperature(r)
        numofcells += 1
   end
end
function inverse(side)
    return HorizonSide(mod(Int(side) + 2, 4))
end
count_avg_temp(r)