text:
# Copyright (c) 2018: Matthew Wilhelm & Matthew Stuber.
# This code is licensed under MIT license (see LICENSE.md for full details)
#############################################################################
# EAGO
# A development environment for robust and global optimization
# See https://github.com/PSORLab/EAGO.jl
#############################################################################
# src/eago_optimizer/optimizer.jl
# Defines optimizer structure used by EAGO. Namely, ObjectiveType, ProblemType
# EAGOParameters, InputProblem, ParsedProblem, and Optimizer.
#############################################################################

"""
$(TYPEDEF)

An abstract type the subtypes of which are associated with functions method
overloaded for for new extensions. An instance of the `DefaultExt <:ExtensionType`
structure to the `Optimizer` in the `ext_type` field.
"""
abstract type ExtensionType end
struct DefaultExt <: ExtensionType end
MOIU.map_indices(::Function, x::ExtensionType) = x
MOIU.map_indices(::Function, x::DefaultExt) = x

@enum(ObjectiveType, UNSET, SINGLE_VARIABLE, SCALAR_AFFINE, SCALAR_QUADRATIC, NONLINEAR)
@enum(ProblemType, UNCLASSIFIED, LP, MILP, SOCP, MISOCP, DIFF_CVX, MINCVX)

export EAGOParameters
"""
$(TYPEDEF)

Storage for parameters that do not change during a global solve.

$(TYPEDFIELDS)
"""
Base.@kwdef mutable struct EAGOParameters

    # Presolving options
    "Should EAGO attempt to remove type assert issues for user-defined functions (default = false)"
    presolve_scrubber_flag::Bool = false
    "Create and use DAG representations of user-defined function (default = false)."
    presolve_to_JuMP_flag::Bool = false
    "[FUTURE FEATURE, NOT CURRENTLY IMPLEMENTED] Apply the epigraph reformulation
    to the problem (default = false)."
    presolve_epigraph_flag::Bool = false
    "Rerranges the DAG using registered transformations (default = false)"
    presolve_flatten_flag::Bool = false

    # Conic reformulations
    "Attempt to bridge convex constraint to second order cone"
    conic_convert_quadratic::Bool = false

    # Iteration logging options
    "Turns logging on records global bounds, node count and run time. Additional
     options are available for recording information specific to subproblems (default = false)."
    log_on::Bool = false
    "Turns on logging of times and feasibility of subproblems (default = false)"
    log_subproblem_info::Bool = false
    "Log data every `log_interval` iterations (default = 1)."
    log_interval::Int64 = 1

    # Optimizer display options
    "The amount of information that should be printed to console while solving
    values range from 0 - 4: 0 is silent, 1 shows iteration summary statistics
    only, 2-4 show varying degrees of details about calculations within each
    iteration (default = 1)."
    verbosity::Int64 = 1
    "Display summary of iteration to console every `output_iterations` (default = 10)"
    output_iterations::Int64 = 1000
    "Display header for summary to console every `output_iterations` (default = 100)"
    header_iterations::Int64 = 10000

    # Node branching options
    "Convex coefficient used to select branch point. Branch point is given by
    `branch_cvx_factor*xmid + (1-branch_cvx_factor)*xsol` (default = 0.25)"
    branch_cvx_factor::Float64 = 0.25
    "Minimum distance from bound to have branch point normalized by width of
    dimension to branch on (default = 0.15)"
    branch_offset::Float64 = 0.15
    "Variables to branch on (default is all nonlinear)."
    branch_variable::Vector{Bool} = Bool[]
    "[FUTURE FEATURE, NOT CURRENTLY IMPLEMENTED] Number of times repeat node
    processing priorto branching (default = 4)."
    branch_max_repetitions::Int64 = 4
    "[FUTURE FEATURE, NOT CURRENTLY IMPLEMENTED] Volume ratio tolerance required
    to repeat processing the current node (default = 0.9)"
    branch_repetition_tol::Float64 = 0.9

    # Termination limits
    "Maximum number of nodes (default = 1E-7)"
    node_limit::Int64 = 1*10^7
    "Maximum CPU time in seconds (default = 1000)"
    time_limit::Float64 = 1000.0
    "Maximum number of iterations (default 3E6)"
    iteration_limit::Int64 = 3*10^6
    "Absolute tolerance for termination (default = 1E-3)"
    absolute_tolerance::Float64 = 1E-3
    "Relative tolerance for termination (default = 1E-3)"
    relative_tolerance::Float64 = 1E-3
    "Absolute constraint feasibility tolerance"
    absolute_constraint_feas_tolerance::Float64 = 1E-6
    "Perform only a local solve of the problem (default = false)."
    local_solve_only::Bool = false
    "[TO BE REMOVED] Flag stops B&B loop if feasible point found."
    feasible_local_continue::Bool = false

    # Options for constraint propagation
    "Depth in B&B tree above which constraint propagation should be disabled (default = 1000)"
    cp_depth::Int64 = 20
    "Number of times to repeat forward-reverse pass routine (default = 3)"
    cp_repetitions::Int64 = 3
    "Disable constraint propagation if the ratio of new node volume to beginning node volume exceeds
    this number (default = 0.99)"
    cp_tolerance::Float64 = 0.99
    "Use only valid interval bounds during constraint propagation (default = false)"
    cp_interval_only::Bool = false

    # obbt options
    "Depth in B&B tree above which OBBT should be disabled (default = 6)"
    obbt_depth::Int64 = 4
    "Number of repetitions of OBBT to perform in preprocessing (default = 3)"
    obbt_repetitions::Int64 = 20
    "Turn aggresive OBBT on (default = false)"
    obbt_aggressive_on::Bool = true
    "Maximum iteration to perform aggresive OBBT (default = 2)"
    obbt_aggressive_max_iteration::Int64 = 2
    "Minimum dimension to perform aggresive OBBT (default = 2)"
    obbt_aggressive_min_dimension::Int64 = 2
    "Tolerance to consider bounds equal (default = 1E-9)"
    obbt_tolerance::Float64 = 1E-9

    # Options for linear bound tightening
    "Depth in B&B tree above which linear FBBT should be disabled (default = 1000)"
    fbbt_lp_depth::Int64  = 1000
    "Number of repetitions of linear FBBT to perform in preprocessing (default = 3)"
    fbbt_lp_repetitions::Int64  = 3

    # Options for quadratic bound tightening
    "[FUTURE FEATURE, NOT CURRENTLY IMPLEMENTED] Depth in B&B tree above which univariate quadratic FBBT should be disabled (default = -1)"
    quad_uni_depth::Int64 = -1
    "[FUTURE FEATURE, NOT CURRENTLY IMPLEMENTED] Number of repetitions of univariate quadratic FBBT to perform in preprocessing (default = 2)"
    quad_uni_repetitions::Int64 = 2
    "[FUTURE FEATURE, NOT CURRENTLY IMPLEMENTED] Depth in B&B tree above which bivariate
    quadratic FBBT should be disabled (default = -1)"
    quad_bi_depth::Int64 = -1
    "[FUTURE FEATURE, NOT CURRENTLY IMPLEMENTED] Number of repetitions of bivariate quadratic FBBT to perform in preprocessing (default = 2)."
    quad_bi_repetitions::Int64 = 2

    # Duality-based bound tightening (DBBT) options
    "Depth in B&B tree above which duality-based bound tightening should be disabled (default = 1E10)"
    dbbt_depth::Int64 = 10^10
    "New bound is considered equal to the prior bound if within dbbt_tolerance (default = 1E-9)."
    dbbt_tolerance::Float64 = 1E-8

    # Subgradient tightening flag
    "Relax Tag used to specify type of McCormick operator"
    relax_tag::RelaxTag = NS()
    "Perform tightening of interval bounds using subgradients at each factor in
    each nonlinear tape during a forward pass (default = true)."
    subgrad_tighten::Bool = true
    "Perform tightening of interval bounds using subgradients at each factor in
    each nonlinear tape during a reverse pass (default = false)."
    reverse_subgrad_tighten::Bool = false
    "Outer round computed subgradient bounds by this amount"
    subgrad_tol::Float64 = 1E-10

    # Tolerance to add cuts and max number of cuts
    "Minimum number of cuts at each node to attempt (unsafe cuts not necessarily added)"
    cut_min_iterations::Int64 = 1
    "Maximum number of cuts at each node to attempt"
    cut_max_iterations::Int64 = 3
    "Convex coefficient used to select point for new added cuts. Branch point is
    given by `(1-cut_cvx)*xmid + cut_cvx*xsol` (default = 0.9)."
    cut_cvx::Float64 = 0.9
    "Add cut if the L1 distance from the prior cutting point to the new cutting
    point normalized by the box volume is greater than the tolerance (default = 0.05)."
    cut_tolerance::Float64 = 0.05
    "Adds an objective cut to the relaxed problem (default = true)."
    objective_cut_on::Bool = true

    "Use tolerances to determine safe cuts in a Khajavirad 2018 manner"
    cut_safe_on::Bool = true
    "Lower tolerance for safe-lp cut, Khajavirad 2018"
    cut_safe_l::Float64 = 1E-8
    "Upper tolerance for safe-lp cut, Khajavirad 2018"
    cut_safe_u::Float64 = 1E8
    "Constant tolerance for safe-lp cut, Khajavirad 2018"
    cut_safe_b::Float64 = 1E9

    "Solve upper problem for every node with depth less than `upper_bounding_depth`
    and with a probabilityof (1/2)^(depth-upper_bounding_depth) otherwise (default = 6)"
    upper_bounding_depth::Int64 = 6

    # handling for domain violations
    "Amount about a domain violation to ignore when propagating bounds."
    domain_violation_guard_on::Bool = false
    "Amount about a domain violation to ignore when propagating bounds."
    domain_violation_ϵ::Float64 = 1E-9
end

"""
$(TYPEDEF)

A structure used to hold objectives and constraints added to EAGO model.
The constraints generally aren't used for relaxations.
"""
Base.@kwdef mutable struct InputProblem

    # variables (set by MOI.add_variable in variables.jl)
    _variable_info::Vector{VariableInfo} = VariableInfo[]
    _variable_count::Int64 = 0

    # last constraint index added
    _last_constraint_index::Int = 0

    # linear constraint storage and count (set by MOI.add_constraint in moi_constraints.jl)
    _linear_leq_constraints::Vector{Tuple{SAF, LT}} = Tuple{SAF, LT}[]
    _linear_geq_constraints::Vector{Tuple{SAF, GT}} = Tuple{SAF, GT}[]
    _linear_eq_constraints::Vector{Tuple{SAF, ET}} = Tuple{SAF, ET}[]

    _linear_leq_count::Int = 0
    _linear_geq_count::Int = 0
    _linear_eq_count::Int = 0

    # quadratic constraint storage and count (set by MOI.add_constraint in moi_constraints.jl)
    _quadratic_leq_constraints::Vector{Tuple{SQF, LT}} = Tuple{SQF, LT}[]
    _quadratic_geq_constraints::Vector{Tuple{SQF, GT}} = Tuple{SQF, GT}[]
    _quadratic_eq_constraints::Vector{Tuple{SQF, ET}} = Tuple{SQF, ET}[]

    _quadratic_leq_count::Int = 0
    _quadratic_geq_count::Int = 0
    _quadratic_eq_count::Int = 0

    # conic constraint storage and count (set by MOI.add_constraint in moi_constraints.jl)
    _conic_second_order::Vector{Tuple{VECOFVAR, MOI.SecondOrderCone}} = Tuple{VECOFVAR, MOI.SecondOrderCone}[]

    _conic_second_order_count::Int = 0

    # nonlinear constraint storage
    _nonlinear_count::Int = 0

    # objective information (set by MOI.set(m, ::ObjectiveFunction...) in optimizer.jl)
    _objective_sv::SV = SV(VI(-1))
    _objective_saf::SAF = SAF(SAT[], 0.0)
    _objective_sqf::SQF = SQF(SAT[], SQT[], 0.0)
    _objective_type::ObjectiveType = UNSET

    # nlp constraints (set by MOI.set(m, ::NLPBlockData...) in optimizer.jl)
    _nlp_data::MOI.NLPBlockData = empty_nlp_data()

    # objective sense information (set by MOI.set(m, ::ObjectiveSense...) in optimizer.jl)
    _optimization_sense::MOI.OptimizationSense = MOI.MIN_SENSE
end

function Base.isempty(x::InputProblem)

    is_empty_flag = true
    new_input_problem = InputProblem()

    for field in fieldnames(InputProblem)

        field_value = getfield(x, field)

        if field_value isa Array
            if !isempty(field_value)
                is_empty_flag = false
                break
            end

        elseif field_value isa Number
            if getfield(new_input_problem, field) !== field_value
                is_empty_flag = false
                break
            end

        end
    end

    is_empty_flag &= x._nlp_data.evaluator isa EmptyNLPEvaluator
    is_empty_flag &= !x._nlp_data.has_objective
    is_empty_flag &= isempty(x._nlp_data.constraint_bounds)

    is_empty_flag &= isempty(x._objective_saf.terms)
    is_empty_flag &= x._objective_saf.constant === 0.0

    is_empty_flag &= isempty(x._objective_sqf.quadratic_terms)
    is_empty_flag &= isempty(x._objective_sqf.affine_terms)
    is_empty_flag &= x._objective_sqf.constant === 0.0

    return is_empty_flag
end

"""
$(TYPEDEF)

A structure used to expressions and problem descriptions EAGO uses to formulate
relaxed problems.
"""
Base.@kwdef mutable struct ParsedProblem

    # Problem classification (set in parse_classify_problem!)
    _problem_type::ProblemType = UNCLASSIFIED

    # objectives (set in initial_parse)
    _objective_sv::SV = SV(VI(-1))
    "_objective_saf stores the objective and is used for constructing linear affine cuts
     of any ObjectiveType"
    _objective_saf::SAF = SAF(SAT[], 0.0)
    _objective_saf_parsed::AffineFunctionIneq = AffineFunctionIneq()
    _objective_sqf::BufferedQuadraticIneq = BufferedQuadraticIneq()
    _objective_nl::BufferedNonlinearFunction = BufferedNonlinearFunction()
    _objective_type::ObjectiveType = UNSET

    # objective sense information (set by convert_to_min in parse.jl)
    _optimization_sense::MOI.OptimizationSense = MOI.MIN_SENSE

    # non-single variable constraints (set in initial_parse)
    _saf_leq::Vector{AffineFunctionIneq} = AffineFunctionIneq[]
    _saf_eq::Vector{AffineFunctionEq} = AffineFunctionEq[]
    _sqf_leq::Vector{BufferedQuadraticIneq} = BufferedQuadraticIneq[]
    _sqf_eq::Vector{BufferedQuadraticEq} = BufferedQuadraticEq[]
    _conic_second_order::Vector{BufferedSOC} = BufferedSOC[]

    # count of non-single variable constraints (set in initial_parse)
    _saf_leq_count::Int = 0
    _saf_eq_count::Int = 0
    _sqf_leq_count::Int = 0
    _sqf_eq_count::Int = 0
    _conic_second_order_count::Int = 0

    # nlp constraints (set in initial_parse)
    _nlp_data::MOI.NLPBlockData = empty_nlp_data()

    # storage for nonlinear functions
    _nonlinear_constr::Vector{BufferedNonlinearFunction} = BufferedNonlinearFunction[]

    # nonlinear constraint storage
    _nonlinear_count::Int = 0

    # nonlinear evaluator
    _relaxed_evaluator = Evaluator()

    # variables (set in initial_parse)
    _variable_info::Vector{VariableInfo} = VariableInfo[]
    _variable_count::Int = 0

    # count of single variable constraint types (set in load_relaxed_problem!)
    _var_leq_count::Int = 0
    _var_geq_count::Int = 0
    _var_eq_count::Int = 0
end

function Base.isempty(x::ParsedProblem)

    is_empty_flag = true

    new_input_problem = ParsedProblem()
    for field in fieldnames(ParsedProblem)

        field_value = getfield(x, field)
        if field_value isa Array
            if !isempty(field_value)
                is_empty_flag = false
                break
            end

        elseif field_value isa Number
            if getfield(new_input_problem, field) !== field_value
                is_empty_flag = false
                break
            end
        end
    end

    is_empty_flag &= x._nlp_data.evaluator isa EmptyNLPEvaluator
    is_empty_flag &= !x._nlp_data.has_objective
    is_empty_flag &= isempty(x._nlp_data.constraint_bounds)

    is_empty_flag &= isempty(x._objective_saf.terms)
    is_empty_flag &= x._objective_saf.constant === 0.0

    is_empty_flag &= isempty(x._objective_saf.terms)
    is_empty_flag &= x._objective_saf.constant === 0.0
    is_empty_flag &= isempty(x._objective_saf_parsed.terms)
    is_empty_flag &= x._objective_saf_parsed.constant === 0.0
    is_empty_flag &= x._objective_saf_parsed.len === 0

    is_empty_flag &= isempty(x._objective_sqf.func.quadratic_terms)
    is_empty_flag &= isempty(x._objective_sqf.func.affine_terms)
    is_empty_flag &= x._objective_sqf.func.constant === 0.0

    is_empty_flag &= isempty(x._objective_sqf.buffer)
    is_empty_flag &= isempty(x._objective_sqf.saf.terms)
    is_empty_flag &= x._objective_sqf.saf.constant === 0.0
    is_empty_flag &= x._objective_sqf.len === 0

    return is_empty_flag
end

function default_nlp_solver()

    upper_optimizer = Ipopt.Optimizer()

    MOI.set(upper_optimizer, MOI.RawParameter("max_iter"),3000)
    MOI.set(upper_optimizer, MOI.RawParameter("acceptable_tol"), 1E30)
    MOI.set(upper_optimizer, MOI.RawParameter("acceptable_iter"), 300)
    MOI.set(upper_optimizer, MOI.RawParameter("constr_viol_tol"), 0.000001)
    MOI.set(upper_optimizer, MOI.RawParameter("acceptable_compl_inf_tol"), 0.000001)
    MOI.set(upper_optimizer, MOI.RawParameter("acceptable_dual_inf_tol"), 1.0)
    MOI.set(upper_optimizer, MOI.RawParameter("acceptable_constr_viol_tol"), 0.000001)
    MOI.set(upper_optimizer, MOI.RawParameter("print_level"), 0)

    return upper_optimizer
end

export Optimizer
"""
$(TYPEDEF)

The main optimizer object used by EAGO to solve problems during the optimization
routine. The following commonly used options are described below and can be set
via keyword arguments in the JuMP/MOI model. The raw parameter interface however
is likely preferable. The Optimizer is organized in the following manner. Parameters
which are expected to be constant over the entire solve are stored in
`_parameters::EAGOParameters` field. User-facing keywords not in EAGOParameters field:
- `relaxed_optimizer::MOI.AbstractOptimizer`: An instance of the optimizer used to solve the relaxed subproblems (default = GLPK.Optimizer())
- `obbt_variable_values::Vector{Bool}`: Variables to perform OBBT on (default: all variables in nonlinear expressions).
- `upper_optimizer::MOI.AbstractOptimizer`: Optimizer used to solve upper bounding problems. (default = Ipopt.Optimizer)
- `enable_optimize_hook::Bool`: Specifies that the optimize_hook! function should be called rather than throw the problem to the standard B&B routine (default = false).
- `ext::Dict{Symbol, Any}`: Holds additional storage needed for constructing extensions to EAGO (default = Dict{Symbol,Any}).
- `ext_type::ExtensionType`: Holds an instance of a subtype of `EAGO.ExtensionType` used to define new custom subroutines (default = DefaultExt()).
"""
Base.@kwdef mutable struct Optimizer <: MOI.AbstractOptimizer

    # Options for optimality-based bound tightening
    # set as a user-specified option
    relaxed_optimizer::MOI.AbstractOptimizer = GLPK.Optimizer()

    # set as a user-specified option (if empty set to all nonlinear by TODO in TODO)
    obbt_variable_values::Vector{Bool} = Bool[]

    # Upper bounding options (set as a user-specified option)
    upper_optimizer::MOI.AbstractOptimizer = default_nlp_solver()

    # Extensions (set as user-specified option)
    enable_optimize_hook::Bool = false
    ext::Dict{Symbol, Any} = Dict{Symbol,Any}()
    ext_type::ExtensionType = DefaultExt()

    # set as user-specified option
    _parameters::EAGOParameters = EAGOParameters()

    # set by MOI manipulations (see Input problem structure)
    _input_problem::InputProblem = InputProblem()

    # loaded from _input_problem by TODO
    _working_problem::ParsedProblem = ParsedProblem()

    _termination_status_code::MOI.TerminationStatusCode = MOI.OPTIMIZE_NOT_CALLED
    _result_status_code::MOI.ResultStatusCode = MOI.OTHER_RESULT_STATUS

    _stack::BinaryMinMaxHeap{NodeBB} = BinaryMinMaxHeap{NodeBB}()

    # set in node_selection!
    _current_node::NodeBB = NodeBB()

    _first_relax_point_set::Bool = false
    _current_xref::Vector{Float64} = Float64[]
    _candidate_xref::Vector{Float64} = Float64[]

    _use_prior_objective_xref::Bool = false
    _current_objective_xref::Vector{Float64} = Float64[]
    _prior_objective_xref::Vector{Float64} = Float64[]

    # set in label_branch_variables! and label_fixed_variables! respectively in parse.jl
    _user_branch_variables::Bool = false
    _fixed_variable::Vector{Bool} = Bool[]
    _branch_variable_count::Int = 0
    _branch_to_sol_map::Vector{Int} = Int[]
    _sol_to_branch_map::Vector{Int} = Int[]

    _continuous_solution::Vector{Float64} = Float64[]

    # all subproblem immutable subproblem status are set in global_solve in corresponding routines
    # in optimize_nonconvex.jl
    _preprocess_feasibility::Bool = true
    _preprocess_result_status::MOI.ResultStatusCode = MOI.OTHER_RESULT_STATUS
    _preprocess_termination_status::MOI.TerminationStatusCode = MOI.OPTIMIZE_NOT_CALLED

    _lower_result_status::MOI.ResultStatusCode = MOI.OTHER_RESULT_STATUS
    _lower_termination_status::MOI.TerminationStatusCode = MOI.OPTIMIZE_NOT_CALLED
    _lower_feasibility::Bool = true
    _lower_objective_value::Float64 = -Inf

    # set in TODO
    _lower_solution::Vector{Float64} = Float64[]
    _lower_lvd::Vector{Float64} = Float64[]
    _lower_uvd::Vector{Float64} = Float64[]

    _cut_result_status::MOI.ResultStatusCode = MOI.OTHER_RESULT_STATUS
    _cut_termination_status::MOI.TerminationStatusCode = MOI.OPTIMIZE_NOT_CALLED
    _cut_objective_value::Float64 = -Inf
    _cut_feasibility::Bool = true

    # set in TODO
    _cut_solution::Vector{Float64} = Float64[]

    _upper_result_status::MOI.ResultStatusCode = MOI.OTHER_RESULT_STATUS
    _upper_termination_status::MOI.TerminationStatusCode = MOI.OPTIMIZE_NOT_CALLED
    _upper_feasibility::Bool = true
    _upper_objective_value::Float64 = Inf

    # array is initialized to correct size in TODO, reset in single_nlp_solve! in optimize_convex.jl
    _upper_variables::Vector{VI} =  VI[]

    # set in TODO
    _upper_solution::Vector{Float64} = Float64[]

    _postprocess_feasibility::Bool = true

    # set to time limit in initial_parse! in parse.jl, decremented throughout global_solve in optimize_nonconvex.jl
    _time_left::Float64 = 1000.0

    # set constructor reset on empty! and  to zero in initial parse! in parse.jl
    _start_time::Float64 = 0.0
    _run_time::Float64 = 0.0
    _parse_time::Float64 = 0.0
    _presolve_time::Float64 = 0.0
    _last_preprocess_time::Float64 = 0.0
    _last_lower_problem_time::Float64 = 0.0
    _last_upper_problem_time::Float64 = 0.0
    _last_postprocessing_time::Float64 = 0.0

    # reset in initial_parse! in parse.jl
    _min_converged_value::Float64 = Inf
    _global_lower_bound::Float64 = -Inf
    _global_upper_bound::Float64 = Inf
    _maximum_node_id::Int64 = 0
    _iteration_count::Int64 = 0
    _node_count::Int64 = 0

    # Storage for output, reset in initial_parse! in parse.jl
    _solution_value::Float64 = 0.0
    _feasible_solution_found::Bool = false
    _first_solution_node::Int64 = -1
    _objective_value::Float64 = -Inf
    _best_upper_value::Float64 = Inf

    # Optimality-Based Bound Tightening (OBBT) Options
    # set in TODO
    _obbt_working_lower_index::Vector{Bool} = Bool[]
    _obbt_working_upper_index::Vector{Bool} = Bool[]
    _lower_indx_diff::Vector{Bool} = Bool[]
    _upper_indx_diff::Vector{Bool} = Bool[]
    _old_low_index::Vector{Bool} = Bool[]
    _old_upp_index::Vector{Bool} = Bool[]
    _new_low_index::Vector{Bool} = Bool[]
    _new_upp_index::Vector{Bool} = Bool[]
    _obbt_variables::Vector{VI} = VI[]
    _obbt_variable_count::Int = 0
    _obbt_performed_flag::Bool = false

    # Buffers for fbbt, set in presolve, used in preprocess
    _lower_fbbt_buffer::Vector{Float64} = Float64[]
    _upper_fbbt_buffer::Vector{Float64} = Float64[]

    # Feasibility-Based Bound Tightening Options
    # set in set_constraint_propagation_fbbt in domain_reduction.jl
    _cp_improvement::Float64 = 0.0
    _cp_evaluation_reverse::Bool = false

    _cut_iterations::Int64 = 0
    _cut_add_flag::Bool = false

    # Options for Repetition (If DBBT Performed Well)
    # set in within preprocess in optimize_nonconvex.jl
    _node_repetitions::Int64 = 0
    _initial_volume::Float64 = 0.0
    _final_volume::Float64 = 0.0

    # Log
    _log::Log = Log()

    # set in TODO
    _buffered_quadratic_ineq_ci::Vector{CI{SAF,LT}} = CI{SAF,LT}[]
    _buffered_quadratic_eq_ci::Vector{CI{SAF,LT}} = CI{SAF,LT}[]

    _buffered_nonlinear_ci::Vector{CI{SAF,LT}} = CI{SAF,LT}[]

    # set initially in TODO, reset in objective_cut in relax.jl
    _objective_cut_ci_sv::CI{SV,LT} = CI{SV,LT}(-1)

    # initialized to empty in constructor (or via MOI.empty), filled in objective_cut in relax.jl
    # called by obbt in domain_reduction.jl, lower_problem, and add_cut in optimize_nonconvex.jl,
    # emptied in delete_objective_cuts! in relax.jl
    _objective_cut_ci_saf::Vector{CI{SAF,LT}} = CI{SAF,LT}[]

    # need to retreive primal _relaxed_variable_index
    # set in TODO
    #"Number of variables actively branched on in B&B routine (excludes linear and fixed)"
    _relaxed_variable_number::Int = 0
    _relaxed_variable_index::Vector{VI} = VI[]
    _relaxed_variable_eq::Vector{Tuple{CI{SV, ET}, Int}} = Tuple{CI{SV, ET}, Int}[]
    _relaxed_variable_lt::Vector{Tuple{CI{SV, LT}, Int}} = Tuple{CI{SV, LT}, Int}[]
    _relaxed_variable_gt::Vector{Tuple{CI{SV, GT}, Int}} = Tuple{CI{SV, GT}, Int}[]

    # set as user-input
    _branch_variables::Vector{Bool} = Bool[]

    _new_eval_constraint::Bool = false
    _new_eval_objective::Bool = false

    _node_to_sv_leq_ci::Vector{CI{SV,LT}} = CI{SV,LT}[]
    _node_to_sv_geq_ci::Vector{CI{SV,GT}} = CI{SV,GT}[]

    #"Set to true if a nonlinear evaluator was created (NLconstraint or NLobjective specified)"
    _nonlinear_evaluator_created::Bool = false

    #_relaxed_evaluator::Evaluator = Evaluator{1,NS}()
    #_relaxed_constraint_bounds::Vector{MOI.NLPBoundsPair} = Vector{MOI.NLPBoundsPair}[]
end

#####
#####
##### General MOI utilities required
#####
#####

const EAGO_OPTIMIZER_ATTRIBUTES = Symbol[:relaxed_optimizer, :relaxed_optimizer_kwargs, :upper_optimizer,
                                         :enable_optimize_hook, :ext, :ext_type, :_parameters]
const EAGO_MODEL_STRUCT_ATTRIBUTES = Symbol[:_stack, :_log, :_current_node, :_working_problem, :_input_problem]
const EAGO_MODEL_NOT_STRUCT_ATTRIBUTES = setdiff(fieldnames(Optimizer), union(EAGO_OPTIMIZER_ATTRIBUTES,
                                                                              EAGO_MODEL_STRUCT_ATTRIBUTES))

function MOI.empty!(m::Optimizer)

    # create a new empty optimizer and copy fields to m
    new_optimizer = Optimizer()
    for field in union(EAGO_MODEL_STRUCT_ATTRIBUTES, EAGO_MODEL_NOT_STRUCT_ATTRIBUTES)
        setfield!(m, field, getfield(new_optimizer, field))
    end

    return nothing
end

function MOI.is_empty(m::Optimizer)

    is_empty_flag = uninitialized(m._current_node)
    is_empty_flag &= isempty(m._stack)
    is_empty_flag &= isempty(m._log)
    is_empty_flag &= isempty(m._input_problem)
    is_empty_flag &= isempty(m._working_problem)

    new_optimizer = Optimizer()
    for field in EAGO_MODEL_NOT_STRUCT_ATTRIBUTES
        if getfield(m, field) != getfield(new_optimizer, field)
            is_empty_flag = false
            break
        end
    end

    return is_empty_flag
end

function MOI.copy_to(model::Optimizer, src::MOI.ModelLike; copy_names = false)
    return MOI.Utilities.default_copy_to(model, src, copy_names)
end

#####
#####
##### Utilities for checking that JuMP model contains variables used in expression
#####
#####

function check_inbounds!(m::Optimizer, vi::VI)
    if !(1 <= vi.value <= m._input_problem._variable_count)
        error("Invalid variable index $vi. ($(m._input_problem._variable_count) variables in the model.)")
    end
    return nothing
end

check_inbounds!(m::Optimizer, var::SV) = check_inbounds!(m, var.variable)

function check_inbounds!(m::Optimizer, aff::SAF)
    for term in aff.terms
        check_inbounds!(m, term.variable_index)
    end
    return nothing
end

function check_inbounds!(m::Optimizer, quad::SQF)
    for term in quad.affine_terms
        check_inbounds!(m, term.variable_index)
    end
    for term in quad.quadratic_terms
        check_inbounds!(m, term.variable_index_1)
        check_inbounds!(m, term.variable_index_2)
    end
    return nothing
end

function check_inbounds!(m::Optimizer, vov::VECOFVAR)
    for vi in vov.variables
        check_inbounds!(m, vi)
    end
    return
end

#####
#####
##### Set & get attributes of model
#####
#####

function MOI.set(m::Optimizer, ::MOI.Silent, value)

     m._parameters.verbosity = 0
     m._parameters.log_on = false
     return nothing

end

function MOI.set(m::Optimizer, p::MOI.RawParameter, value)

    if p.name isa String
        psym = Symbol(p.name)
    elseif p.name isa Symbol
        psym = p.name
    else
        error("EAGO only supports raw parameters with Symbol or String names.")
    end

    if psym in fieldnames(EAGOParameters)
        setfield!(m._parameters, psym, value)
    else
        setfield!(m, psym, value)
    end

    return nothing
end

function MOI.set(m::Optimizer, ::MOI.TimeLimitSec, value::Nothing)
    m._parameters.time_limit = Inf
    return nothing
end

function MOI.set(m::Optimizer, ::MOI.TimeLimitSec, value::Float64)
    m._parameters.time_limit = value
    return nothing
end

function MOI.get(m::Optimizer, ::MOI.ListOfVariableIndices)
    return [MOI.VariableIndex(i) for i = 1:length(m._input_problem._variable_info)]
end

function MOI.get(m::Optimizer, ::MOI.ObjectiveValue)
    mult = 1.0
    if m._input_problem._optimization_sense === MOI.MAX_SENSE
        mult *= -1.0
    end
    return mult*m._objective_value
end

MOI.get(m::Optimizer, ::MOI.NumberOfVariables) = m._input_problem._variable_count

function MOI.get(m::Optimizer, ::MOI.ObjectiveBound)
    if m._input_problem._optimization_sense === MOI.MAX_SENSE
        bound = -m._global_lower_bound
    else
        bound = m._global_upper_bound
    end
    return bound
end

function MOI.get(m::Optimizer, ::MOI.RelativeGap)
    LBD = m._global_lower_bound
    UBD = m._global_upper_bound
    if m._input_problem._optimization_sense === MOI.MAX_SENSE
        gap = abs(UBD - LBD)/abs(LBD)
    else
        gap = abs(UBD - LBD)/abs(UBD)
    end
    return gap
end

MOI.get(m::Optimizer, ::MOI.SolverName) = "EAGO: Easy Advanced Global Optimization"
MOI.get(m::Optimizer, ::MOI.TerminationStatus) = m._termination_status_code
MOI.get(m::Optimizer, ::MOI.PrimalStatus) = m._result_status_code
MOI.get(m::Optimizer, ::MOI.SolveTime) = m._run_time
MOI.get(m::Optimizer, ::MOI.NodeCount) = m._maximum_node_id
MOI.get(m::Optimizer, ::MOI.ResultCount) = (m._result_status_code === MOI.FEASIBLE_POINT) ? 1 : 0
MOI.get(m::Optimizer, ::MOI.TimeLimitSec) = m.time_limit

function MOI.get(model::Optimizer, ::MOI.VariablePrimal, vi::MOI.VariableIndex)
    check_inbounds!(model, vi)
    return model._continuous_solution[vi.value]
end

function MOI.get(m::Optimizer, p::MOI.RawParameter)
    if p.name isa String
        psym = Symbol(p.name)
    elseif p.name isa Symbol
        psym = p.name
    else
        error("EAGO only supports raw parameters with Symbol or String names.")
    end

    if psym in fieldnames(EAGOParameters)
        return getfield(m._parameters, psym)
    else
        return getfield(m, psym)
    end
end

#####
#####
##### Support, set, and evaluate objective functions
#####
#####
MOI.supports(::Optimizer, ::MOI.TimeLimitSec) = true
MOI.supports(::Optimizer, ::MOI.ObjectiveSense) = true
MOI.supports(::Optimizer, ::MOI.ObjectiveFunction{F}) where {F <: Union{SV, SAF, SQF}} = true

function MOI.set(m::Optimizer, ::MOI.NLPBlock, nlp_data::MOI.NLPBlockData)
    if nlp_data.has_objective
        m._input_problem._objective_type = NONLINEAR
    end
    m._input_problem._nlp_data = nlp_data
    return nothing
end

function MOI.set(m::Optimizer, ::MOI.ObjectiveFunction{SV}, func::SV)
    check_inbounds!(m, func)
    m._input_problem._objective_sv = func
    m._input_problem._objective_type = SINGLE_VARIABLE
    return nothing
end

function MOI.set(m::Optimizer, ::MOI.ObjectiveFunction{SAF}, func::SAF)
    check_inbounds!(m, func)
    m._input_problem._objective_saf = func
    m._input_problem._objective_type = SCALAR_AFFINE
    return nothing
end

function MOI.set(m::Optimizer, ::MOI.ObjectiveFunction{SQF}, func::SQF)
    check_inbounds!(m, func)
    m._input_problem._objective_sqf = func
    m._input_problem._objective_type = SCALAR_QUADRATIC
    return nothing
end

function MOI.set(m::Optimizer, ::MOI.ObjectiveSense, sense::MOI.OptimizationSense)
    m._input_problem._optimization_sense = sense
    return nothing
end


text:
###############################################################################
#
#   fmpz_rel_series.jl : Power series over flint fmpz integers
#
###############################################################################

export fmpz_rel_series, FmpzRelSeriesRing

###############################################################################
#
#   Data type and parent object methods
#
###############################################################################

function O(a::fmpz_rel_series)
   val = pol_length(a) + valuation(a) - 1
   val < 0 && throw(DomainError())
   z = fmpz_rel_series(Array{fmpz}(0), 0, val, val)
   z.parent = parent(a)
   return z
end

elem_type(::Type{FmpzRelSeriesRing}) = fmpz_rel_series

parent_type(::Type{fmpz_rel_series}) = FmpzRelSeriesRing

base_ring(R::FmpzRelSeriesRing) = R.base_ring

var(a::FmpzRelSeriesRing) = a.S

###############################################################################
#
#   Basic manipulation
#
###############################################################################

max_precision(R::FmpzRelSeriesRing) = R.prec_max

function normalise(a::fmpz_rel_series, len::Int)
   if len > 0
      c = fmpz()
      ccall((:fmpz_poly_get_coeff_fmpz, :libflint), Void,
         (Ref{fmpz}, Ref{fmpz_rel_series}, Int), c, a, len - 1)
   end
   while len > 0 && iszero(c)
      len -= 1
      if len > 0
         ccall((:fmpz_poly_get_coeff_fmpz, :libflint), Void,
            (Ref{fmpz}, Ref{fmpz_rel_series}, Int), c, a, len - 1)
      end
   end
   return len
end

function pol_length(x::fmpz_rel_series)
   return ccall((:fmpz_poly_length, :libflint), Int, (Ref{fmpz_rel_series},), x)
end

precision(x::fmpz_rel_series) = x.prec

function polcoeff(x::fmpz_rel_series, n::Int)
   if n < 0
      return fmpz(0)
   end
   z = fmpz()
   ccall((:fmpz_poly_get_coeff_fmpz, :libflint), Void,
         (Ref{fmpz}, Ref{fmpz_rel_series}, Int), z, x, n)
   return z
end

zero(R::FmpzRelSeriesRing) = R(0)

one(R::FmpzRelSeriesRing) = R(1)

function gen(R::FmpzRelSeriesRing)
   z = fmpz_rel_series([fmpz(1)], 1, max_precision(R) + 1, 1)
   z.parent = R
   return z
end

function deepcopy_internal(a::fmpz_rel_series, dict::ObjectIdDict)
   z = fmpz_rel_series(a)
   z.prec = a.prec
   z.val = a.val
   z.parent = parent(a)
   return z
end

function renormalize!(z::fmpz_rel_series)
   i = 0
   zlen = pol_length(z)
   zval = valuation(z)
   zprec = precision(z)
   while i < zlen && polcoeff(z, i) == 0
      i += 1
   end
   z.prec = zprec
   if i == zlen
      z.val = zprec
   else
      z.val = zval + i
      ccall((:fmpz_poly_shift_right, :libflint), Void,
            (Ref{fmpz_rel_series}, Ref{fmpz_rel_series}, Int), z, z, i)
   end
   return nothing
end

###############################################################################
#
#   String I/O
#
###############################################################################

function show(io::IO, a::FmpzRelSeriesRing)
   print(io, "Univariate power series ring in ", var(a), " over ")
   show(io, base_ring(a))
end

show_minus_one(::Type{fmpz_rel_series}) = show_minus_one(fmpz)

###############################################################################
#
#   Unary operators
#
###############################################################################

function -(x::fmpz_rel_series)
   z = parent(x)()
   ccall((:fmpz_poly_neg, :libflint), Void,
                (Ref{fmpz_rel_series}, Ref{fmpz_rel_series}),
               z, x)
   z.prec = x.prec
   z.val = x.val
   return z
end

###############################################################################
#
#   Binary operators
#
###############################################################################

function +(a::fmpz_rel_series, b::fmpz_rel_series)
   check_parent(a, b)
   lena = pol_length(a)
   lenb = pol_length(b)
   prec = min(a.prec, b.prec)
   val = min(a.val, b.val)
   lena = min(lena, prec - a.val)
   lenb = min(lenb, prec - b.val)
   z = parent(a)()
   if a.val < b.val
      lenz = max(lena, lenb + b.val - a.val)
      ccall((:fmpz_poly_set_trunc, :libflint), Void,
            (Ref{fmpz_rel_series}, Ref{fmpz_rel_series}, Int),
            z, b, max(0, lenz - b.val + a.val))
      ccall((:fmpz_poly_shift_left, :libflint), Void,
            (Ref{fmpz_rel_series}, Ref{fmpz_rel_series}, Int),
            z, z, b.val - a.val)
      ccall((:fmpz_poly_add_series, :libflint), Void,
                (Ref{fmpz_rel_series}, Ref{fmpz_rel_series}, Ref{fmpz_rel_series}, Int),
               z, z, a, lenz)
   elseif b.val < a.val
      lenz = max(lena + a.val - b.val, lenb)
      ccall((:fmpz_poly_set_trunc, :libflint), Void,
            (Ref{fmpz_rel_series}, Ref{fmpz_rel_series}, Int),
            z, a, max(0, lenz - a.val + b.val))
      ccall((:fmpz_poly_shift_left, :libflint), Void,
            (Ref{fmpz_rel_series}, Ref{fmpz_rel_series}, Int),
            z, z, a.val - b.val)
      ccall((:fmpz_poly_add_series, :libflint), Void,
                (Ref{fmpz_rel_series}, Ref{fmpz_rel_series}, Ref{fmpz_rel_series}, Int),
               z, z, b, lenz)
   else
      lenz = max(lena, lenb)
      ccall((:fmpz_poly_add_series, :libflint), Void,
                (Ref{fmpz_rel_series}, Ref{fmpz_rel_series}, Ref{fmpz_rel_series}, Int),
               z, a, b, lenz)
   end
   z.prec = prec
   z.val = val
   renormalize!(z)
   return z
end

function -(a::fmpz_rel_series, b::fmpz_rel_series)
   check_parent(a, b)
   lena = pol_length(a)
   lenb = pol_length(b)
   prec = min(a.prec, b.prec)
   val = min(a.val, b.val)
   lena = min(lena, prec - a.val)
   lenb = min(lenb, prec - b.val)
   lenz = max(lena, lenb)
   z = parent(a)()
   if a.val < b.val
      lenz = max(lena, lenb + b.val - a.val)
      ccall((:fmpz_poly_set_trunc, :libflint), Void,
            (Ref{fmpz_rel_series}, Ref{fmpz_rel_series}, Int),
            z, b, max(0, lenz - b.val + a.val))
      ccall((:fmpz_poly_shift_left, :libflint), Void,
            (Ref{fmpz_rel_series}, Ref{fmpz_rel_series}, Int),
            z, z, b.val - a.val)
      ccall((:fmpz_poly_neg, :libflint), Void,
            (Ref{fmpz_rel_series}, Ref{fmpz_rel_series}), z, z)
      ccall((:fmpz_poly_add_series, :libflint), Void,
                (Ref{fmpz_rel_series}, Ref{fmpz_rel_series}, Ref{fmpz_rel_series}, Int),
               z, z, a, lenz)
   elseif b.val < a.val
      lenz = max(lena + a.val - b.val, lenb)
      ccall((:fmpz_poly_set_trunc, :libflint), Void,
            (Ref{fmpz_rel_series}, Ref{fmpz_rel_series}, Int),
            z, a, max(0, lenz - a.val + b.val))
      ccall((:fmpz_poly_shift_left, :libflint), Void,
            (Ref{fmpz_rel_series}, Ref{fmpz_rel_series}, Int),
            z, z, a.val - b.val)
      ccall((:fmpz_poly_sub_series, :libflint), Void,
                (Ref{fmpz_rel_series}, Ref{fmpz_rel_series}, Ref{fmpz_rel_series}, Int),
               z, z, b, lenz)
   else
      lenz = max(lena, lenb)
      ccall((:fmpz_poly_sub_series, :libflint), Void,
                (Ref{fmpz_rel_series}, Ref{fmpz_rel_series}, Ref{fmpz_rel_series}, Int),
               z, a, b, lenz)
   end
   z.prec = prec
   z.val = val
   renormalize!(z)
   return z
end

function *(a::fmpz_rel_series, b::fmpz_rel_series)
   check_parent(a, b)
   lena = pol_length(a)
   lenb = pol_length(b)
   aval = valuation(a)
   bval = valuation(b)
   prec = min(a.prec - aval, b.prec - bval)
   lena = min(lena, prec)
   lenb = min(lenb, prec)
   z = parent(a)()
   z.val = a.val + b.val
   z.prec = prec + z.val
   if lena == 0 || lenb == 0
      return z
   end
   lenz = min(lena + lenb - 1, prec)
   ccall((:fmpz_poly_mullow, :libflint), Void,
                (Ref{fmpz_rel_series}, Ref{fmpz_rel_series}, Ref{fmpz_rel_series}, Int),
               z, a, b, lenz)

   return z
end

###############################################################################
#
#   Ad hoc binary operators
#
###############################################################################

function *(x::Int, y::fmpz_rel_series)
   z = parent(y)()
   z.prec = y.prec
   z.val = y.val
   ccall((:fmpz_poly_scalar_mul_si, :libflint), Void,
                (Ref{fmpz_rel_series}, Ref{fmpz_rel_series}, Int),
               z, y, x)
   return z
end

*(x::fmpz_rel_series, y::Int) = y * x

function *(x::fmpz, y::fmpz_rel_series)
   z = parent(y)()
   z.prec = y.prec
   z.val = y.val
   ccall((:fmpz_poly_scalar_mul_fmpz, :libflint), Void,
                (Ref{fmpz_rel_series}, Ref{fmpz_rel_series}, Ref{fmpz}),
               z, y, x)
   return z
end

*(x::fmpz_rel_series, y::fmpz) = y * x

###############################################################################
#
#   Shifting
#
###############################################################################

function shift_left(x::fmpz_rel_series, len::Int)
   len < 0 && throw(DomainError())
   xlen = pol_length(x)
   z = fmpz_rel_series(x)
   z.prec = x.prec + len
   z.val = x.val + len
   z.parent = parent(x)
   return z
end

function shift_right(x::fmpz_rel_series, len::Int)
   len < 0 && throw(DomainError())
   xlen = pol_length(x)
   xval = valuation(x)
   z = parent(x)()
   if len >= xlen + xval
      z.prec = max(0, x.prec - len)
      z.val = max(0, x.prec - len)
   else
      z.prec = max(0, x.prec - len)
      z.val = max(0, xval - len)
      zlen = min(xlen + xval - len, xlen)
      ccall((:fmpz_poly_shift_right, :libflint), Void,
                (Ref{fmpz_rel_series}, Ref{fmpz_rel_series}, Int),
               z, x, xlen - zlen)
      renormalize!(z)
   end
   return z
end

###############################################################################
#
#   Truncation
#
###############################################################################

function truncate(x::fmpz_rel_series, prec::Int)
   prec < 0 && throw(DomainError())
   xlen = pol_length(x)
   xprec = precision(x)
   xval = valuation(x)
   if xprec + xval <= prec
      return x
   end
   z = parent(x)()
   z.prec = prec
   if prec <= xval
      z = parent(x)()
      z.val = prec
      z.prec = prec
   else
      z.val = xval
      ccall((:fmpz_poly_set_trunc, :libflint), Void,
                (Ref{fmpz_rel_series}, Ref{fmpz_rel_series}, Int),
               z, x, min(prec - xval, xlen))
   end
   return z
end

###############################################################################
#
#   Powering
#
###############################################################################

function ^(a::fmpz_rel_series, b::Int)
   b < 0 && throw(DomainError())
   if isgen(a)
      z = parent(a)()
      z = setcoeff!(z, 0, fmpz(1))
      z.prec = a.prec + b - 1
      z.val = b
   elseif pol_length(a) == 0
      z = parent(a)()
      z.prec = b*valuation(a)
      z.val = b*valuation(a)
   elseif pol_length(a) == 1
      return parent(a)([polcoeff(a, 0)^b], 1,
                           (b - 1)*valuation(a) + precision(a), b*valuation(a))
   elseif b == 0
      return one(parent(a))
   else
      z = parent(a)()
      z.prec = a.prec + (b - 1)*valuation(a)
      z.val = b*valuation(a)
      ccall((:fmpz_poly_pow_trunc, :libflint), Void,
                (Ref{fmpz_rel_series}, Ref{fmpz_rel_series}, Int, Int),
               z, a, b, z.prec - z.val)
   end
   return z
end

###############################################################################
#
#   Comparison
#
###############################################################################

function ==(x::fmpz_rel_series, y::fmpz_rel_series)
   check_parent(x, y)
   prec = min(x.prec, y.prec)
   if prec <= x.val && prec <= y.val
      return true
   end
   if x.val != y.val
      return false
   end
   xlen = normalise(x, min(pol_length(x), prec - x.val))
   ylen = normalise(y, min(pol_length(y), prec - y.val))
   if xlen != ylen
      return false
   end
   return Bool(ccall((:fmpz_poly_equal_trunc, :libflint), Cint,
                (Ref{fmpz_rel_series}, Ref{fmpz_rel_series}, Int),
               x, y, xlen))
end

function isequal(x::fmpz_rel_series, y::fmpz_rel_series)
   if parent(x) != parent(y)
      return false
   end
   if x.prec != y.prec || x.val != y.val || pol_length(x) != pol_length(y)
      return false
   end
   return Bool(ccall((:fmpz_poly_equal, :libflint), Cint,
                (Ref{fmpz_rel_series}, Ref{fmpz_rel_series}, Int),
               x, y, pol_length(x)))
end

###############################################################################
#
#   Ad hoc comparisons
#
###############################################################################

function ==(x::fmpz_rel_series, y::fmpz)
   if precision(x) == 0
      return true
   elseif pol_length(x) > 1
      return false
   elseif pol_length(x) == 1
      if x.val == 0
         z = fmpz()
         ccall((:fmpz_poly_get_coeff_fmpz, :libflint), Void,
                       (Ref{fmpz}, Ref{fmpz_rel_series}, Int), z, x, 0)
         return ccall((:fmpz_equal, :libflint), Bool,
               (Ref{fmpz}, Ref{fmpz}, Int), z, y, 0)
      else
         return false
      end
   else
      return iszero(y)
   end
end

==(x::fmpz, y::fmpz_rel_series) = y == x

==(x::fmpz_rel_series, y::Integer) = x == fmpz(y)

==(x::Integer, y::fmpz_rel_series) = y == x

###############################################################################
#
#   Exact division
#
###############################################################################

function divexact(x::fmpz_rel_series, y::fmpz_rel_series)
   check_parent(x, y)
   iszero(y) && throw(DivideError())
   yval = valuation(y)
   xval = valuation(x)
   if yval != 0
      if xval >= yval
         x = shift_right(x, yval)
         y = shift_right(y, yval)
      end
   end
   !isunit(y) && error("Unable to invert power series")
   prec = min(x.prec - x.val, y.prec - y.val)
   z = parent(x)()
   z.val = xval - yval
   z.prec = prec + z.val
   if prec != 0
      ccall((:fmpz_poly_div_series, :libflint), Void,
                (Ref{fmpz_rel_series}, Ref{fmpz_rel_series}, Ref{fmpz_rel_series}, Int),
               z, x, y, prec)
   end
   return z
end

###############################################################################
#
#   Ad hoc exact division
#
###############################################################################

function divexact(x::fmpz_rel_series, y::Int)
   y == 0 && throw(DivideError())
   z = parent(x)()
   z.prec = x.prec
   z.val = x.val
   ccall((:fmpz_poly_scalar_divexact_si, :libflint), Void,
                (Ref{fmpz_rel_series}, Ref{fmpz_rel_series}, Int),
               z, x, y)
   return z
end

function divexact(x::fmpz_rel_series, y::fmpz)
   iszero(y) && throw(DivideError())
   z = parent(x)()
   z.prec = x.prec
   z.prec = x.prec
   z.val = x.val
   ccall((:fmpz_poly_scalar_divexact_fmpz, :libflint), Void,
                (Ref{fmpz_rel_series}, Ref{fmpz_rel_series}, Ref{fmpz}),
               z, x, y)
   return z
end

divexact(x::fmpz_rel_series, y::Integer) = divexact(x, fmpz(y))

###############################################################################
#
#   Inversion
#
###############################################################################

function inv(a::fmpz_rel_series)
   a == 0 && throw(DivideError())
   !isunit(a) && error("Unable to invert power series")
   ainv = parent(a)()
   ainv.prec = a.prec
   ainv.val = 0
   ccall((:fmpz_poly_inv_series, :libflint), Void,
                (Ref{fmpz_rel_series}, Ref{fmpz_rel_series}, Int),
               ainv, a, a.prec)
   return ainv
end

###############################################################################
#
#   Unsafe functions
#
###############################################################################

function zero!(x::fmpz_rel_series)
  ccall((:fmpz_poly_zero, :libflint), Void,
                   (Ref{fmpz_rel_series},), x)
  x.prec = parent(x).prec_max
  return x
end

function setcoeff!(z::fmpz_rel_series, n::Int, x::fmpz)
   ccall((:fmpz_poly_set_coeff_fmpz, :libflint), Void,
                (Ref{fmpz_rel_series}, Int, Ref{fmpz}),
               z, n, x)
   return z
end

function mul!(z::fmpz_rel_series, a::fmpz_rel_series, b::fmpz_rel_series)
   lena = pol_length(a)
   lenb = pol_length(b)
   aval = valuation(a)
   bval = valuation(b)
   prec = min(a.prec - aval, b.prec - bval)
   lena = min(lena, prec)
   lenb = min(lenb, prec)
   z.val = a.val + b.val
   z.prec = prec + c.val
   lenz = min(lena + lenb - 1, prec)
   if lena <= 0 || lenb <= 0
      lenz = 0
   end
   ccall((:fmpz_poly_mullow, :libflint), Void,
                (Ref{fmpz_rel_series}, Ref{fmpz_rel_series}, Ref{fmpz_rel_series}, Int),
               z, a, b, lenz)
   return z
end

function addeq!(a::fmpz_rel_series, b::fmpz_rel_series)
   lena = pol_length(a)
   lenb = pol_length(b)
   prec = min(a.prec, b.prec)
   val = min(a.val, b.val)
   lena = min(lena, prec - a.val)
   lenb = min(lenb, prec - b.val)
   if a.val < b.val
      z = fmpz_rel_series()
      lenz = max(lena, lenb + b.val - a.val)
      ccall((:fmpz_poly_set_trunc, :libflint), Void,
            (Ref{fmpz_rel_series}, Ref{fmpz_rel_series}, Int),
            z, b, max(0, lenz - b.val + a.val))
      ccall((:fmpz_poly_shift_left, :libflint), Void,
            (Ref{fmpz_rel_series}, Ref{fmpz_rel_series}, Int),
            z, z, b.val - a.val)
      ccall((:fmpz_poly_add_series, :libflint), Void,
                (Ref{fmpz_rel_series}, Ref{fmpz_rel_series}, Ref{fmpz_rel_series}, Int),
               a, a, z, lenz)
   elseif b.val < a.val
      lenz = max(lena + a.val - b.val, lenb)
      ccall((:fmpz_poly_truncate, :libflint), Void,
            (Ref{fmpz_rel_series}, Int),
            a, max(0, lenz - a.val + b.val))
      ccall((:fmpz_poly_shift_left, :libflint), Void,
            (Ref{fmpz_rel_series}, Ref{fmpz_rel_series}, Int),
            a, a, a.val - b.val)
      ccall((:fmpz_poly_add_series, :libflint), Void,
                (Ref{fmpz_rel_series}, Ref{fmpz_rel_series}, Ref{fmpz_rel_series}, Int),
               a, a, b, lenz)
   else
      lenz = max(lena, lenb)
      ccall((:fmpz_poly_add_series, :libflint), Void,
                (Ref{fmpz_rel_series}, Ref{fmpz_rel_series}, Ref{fmpz_rel_series}, Int),
               a, a, b, lenz)
   end
   a.prec = prec
   a.val = val
   renormalize!(a)
   return a
end

function add!(c::fmpz_rel_series, a::fmpz_rel_series, b::fmpz_rel_series)
   lena = length(a)
   lenb = length(b)

   prec = min(a.prec, b.prec)

   lena = min(lena, prec)
   lenb = min(lenb, prec)

   lenc = max(lena, lenb)
   c.prec = prec
   ccall((:fmpz_poly_add_series, :libflint), Void,
                (Ref{fmpz_rel_series}, Ref{fmpz_rel_series}, Ref{fmpz_rel_series}, Int),
               c, a, b, lenc)
   return c
end

###############################################################################
#
#   Promotion rules
#
###############################################################################

promote_rule(::Type{fmpz_rel_series}, ::Type{T}) where {T <: Integer} = fmpz_rel_series

promote_rule(::Type{fmpz_rel_series}, ::Type{fmpz}) = fmpz_rel_series

###############################################################################
#
#   Parent object call overload
#
###############################################################################

function (a::FmpzRelSeriesRing)()
   z = fmpz_rel_series()
   z.prec = a.prec_max
   z.val = a.prec_max
   z.parent = a
   return z
end

function (a::FmpzRelSeriesRing)(b::Integer)
   if b == 0
      z = fmpz_rel_series()
      z.prec = a.prec_max
   else
      z = fmpz_rel_series([fmpz(b)], 1, a.prec_max, 0)
   end
   z.parent = a
   return z
end

function (a::FmpzRelSeriesRing)(b::fmpz)
   if b == 0
      z = fmpz_rel_series()
      z.prec = a.prec_max
   else
      z = fmpz_rel_series([b], 1, a.prec_max, 0)
   end
   z.parent = a
   return z
end

function (a::FmpzRelSeriesRing)(b::fmpz_rel_series)
   parent(b) != a && error("Unable to coerce power series")
   return b
end

function (a::FmpzRelSeriesRing)(b::Array{fmpz, 1}, len::Int, prec::Int, val::Int)
   z = fmpz_rel_series(b, len, prec, val)
   z.parent = a
   return z
end

ext: jl


text:
function setupβLSsolverw(optim_algorithm,
    Es,
    Bs,
    As,
    LS_inds,
    U_rad_cost,
    y_cost::Vector{Complex{T}},
    κs_β_DOFs, κs_β_orderings;
    w_lb_default = 0.2,
    w_ub_default = 5.0,
    β_max_iters = 50,
    β_xtol_rel = 1e-9,
    β_ftol_rel = 1e-9,
    β_maxtime = Inf) where T

    N_β = sum( getNβ(κs_β_DOFs[n], Bs[n]) for n = 1:length(Bs) )

    β_lb = ones(T, N_β) .* (-π)
    β_ub = ones(T, N_β) .* (π)

    p_lb = β_lb
    p_ub = β_ub

    q, updateβfunc, updatewfunc, E_BLS, w_BLS, b_BLS,
    getβfunc = setupcostβLSw(Es, Bs, As, LS_inds, U_rad_cost, y_cost,
        κs_β_DOFs, κs_β_orderings;
        w_lb_default = w_lb_default,
        w_ub_default = w_ub_default)

    # q is simulated spectra. f is cost function.
    f = pp->costβLSw(U_rad_cost, y_cost, updateβfunc, updatewfunc, pp, Es,
    E_BLS, w_BLS, b_BLS, q)

    df = xx->FiniteDiff.finite_difference_gradient(f, xx)

    opt = NLopt.Opt(optim_algorithm, N_β)

    run_optim = pp->runNLopt!(opt,
        pp,
        f,
        df,
        p_lb,
        p_ub;
        max_iters = β_max_iters,
        xtol_rel = β_xtol_rel,
        ftol_rel = β_ftol_rel,
        maxtime = β_maxtime)

    return run_optim, f, E_BLS, w_BLS, b_BLS, updateβfunc, updatewfunc, q
end

"""
w version.
"""
function setupcostβLSw(Es::Vector{NMRSpectraSimulator.καFIDModelType{T, SST}},
    Bs,
    As,
    LS_inds,
    U0_rad,
    y0::Vector{Complex{T}},
    κs_β_DOFs, κs_β_orderings;
    w_lb_default = 0.2,
    w_ub_default = 5.0) where {T <: Real, SST}

    #N_β = sum( getNβ(κs_β_DOFs[n], Bs[n]) for n = 1:length(Bs) )
    N_β = sum( getNβ(Bs[n]) for n = 1:length(Bs) )

    st_ind_β = 1
    fin_ind_β = st_ind_β + N_β - 1
    #updateβfunc = pp->updateβ!(Bs, κs_β_orderings, κs_β_DOFs, pp, st_ind_β)
    updateβfunc = pp->updateβ!(Bs, pp, st_ind_β)

    ### LS w.
    U_rad_LS = U0_rad[LS_inds]

    N_w_vars = length(Es)
    E_BLS, w_BLS, b_BLS = setupupdateLS(length(U_rad_LS), N_w_vars, y0[LS_inds])

    w_lb = ones(N_w_vars) .* w_lb_default
    w_ub = ones(N_w_vars) .* w_ub_default

    updatewfunc = xx->updatew!(E_BLS, b_BLS, w_BLS,
    U_rad_LS, Es, As, w_lb, w_ub)

    #### extract parameters from p.
    getβfunc = pp->pp[st_ind_β:fin_ind_β]

    f = uu->NMRSpectraSimulator.evalitpproxymixture(uu, As, Es; w = w_BLS)

    return f, updateβfunc, updatewfunc, E_BLS, w_BLS, b_BLS, getβfunc
end

function costβLSw(U_rad,
    S_U::Vector{Complex{T}},
    updateβfunc, updatewfunc,
    p::Vector{T}, Es,
    E_BLS::Matrix{Complex{T}}, w_BLS::Vector{T}, b_BLS,
    f)::T where T <: Real

    updateβfunc(p)
    updatewfunc(1.0)

    # ## l-2 costfunc.
    # cost = zero(T)
    # for m = 1:length(S_U)
    #
    #     f_u = f(U_rad[m])
    #
    #     cost += abs2( f_u - S_U[m] )
    # end
    #
    # # faster l-2 cost compute.
    # B = reinterpret(T, E_BLS)
    # tmp = B*w_BLS - b_BLS
    # cost = dot(tmp, tmp)
    cost = norm(reinterpret(T, E_BLS)*w_BLS - b_BLS)^2 # compact version.

    return cost
end


text:
@testset "NGramizer" begin

    sample_text = "this is some sample text"
    tkns = TextAnalysis.tokenize(Languages.English(), sample_text)
    ngs = TextAnalysis.ngramize(Languages.English(), tkns, 1)
    @test isequal(ngs, Dict{String,Int}("some" => 1,
    	                                     "this" => 1,
    	                                     "is" => 1,
    	                                     "sample" => 1,
    	                                     "text" => 1))
    ngs = TextAnalysis.ngramize(Languages.English(), tkns, 2)
    @test isequal(ngs, Dict{String,Int}("this is" => 1,
                                             "sample text" => 1,
                                             "is some" => 1,
                                              "some sample" => 1))
    ngs = TextAnalysis.ngramize(Languages.English(), tkns, 1, 2)
    @test isequal(ngs, Dict{String,Int}("this is" => 1,
                                             "is some" => 1,
                                             "some sample" => 1,
                                             "sample text" => 1,
                                             "this" => 1,
                                             "is" => 1,
                                             "some" => 1,
                                             "sample" => 1,
                                             "text" => 1))
end


text:
# Solvers corresponding to solving a closed-form expression

"""
$SIGNATURES

Fit a least square regression either with no penalty (OLS) or with a L2 penalty
(Ridge).

## Complexity

Assuming `n` dominates `p`,

* non-iterative (full solve):     O(np²) - dominated by the construction of the
                                  Hessian X'X.
* iterative (conjugate gradient): O(κnp) - with κ the number of CG steps
                                  (κ ≤ p).
"""
function _fit(glr::GLR{L2Loss,<:L2R}, solver::Analytical, X, y, scratch)
    # full solve
    if !solver.iterative
        λ  = getscale(glr.penalty)
        if iszero(λ)
            # standard LS solution
            return augment_X(X, glr.fit_intercept) \ y
        else
            # Ridge case -- form the Hat Matrix then solve
            H = form_XtX(X, glr.fit_intercept, λ)
            b = X'y
            glr.fit_intercept && (b = vcat(b, sum(y)))
            return cholesky!(H) \ b
        end
    end
    # Iterative case, note that there is no augmentation here
    # it is done implicitly in the application of the Hessian to
    # avoid copying X
    # The number of CG steps to convergence is at most `p`
    _,p,_ = npc(scratch)
    max_cg_steps = min(solver.max_inner, p)
    # Form the Hessian map, cost of application H*v is O(np)
    Hm = LinearMap(Hv!(glr, X, y, scratch), p;
                   ismutating=true, isposdef=true, issymmetric=true)
    b  = X'y
    glr.fit_intercept && (b = vcat(b, sum(y)))
    return cg(Hm, b; maxiter=max_cg_steps)
end


text:
function multivariate_normal_lkj_quote(M, L, T, (track_y, track_μ, track_Σ), partial, lkjinv)
    if partial
        q = quote
            $(Expr(:meta, :inline))
            δ = y - μ
        end
        if !lkjinv
            if partial && track_Σ
                push!(q.args, :(R, ∂R∂L = StructuredMatrices.∂inv′(L)) )
            else
                push!(q.args, :(R = StructuredMatrices.inv′(L)) )
            end
        end
        push!(q.args, :(Rδ = R*δ))
        if track_Σ
            if (Σ_type <: AbstractFloat) && ((y_type <: AbstractVector) || (μ_type <: AbstractVector))
                if y_type <: AbstractVector
                    logdet_expr = :($(M) * LinearAlgebra.logdet(R))
                else # μ_type <: AbstractVector
                    logdet_expr = :($(M) * LinearAlgebra.logdet(R))
                end
            else
                logdet_expr = :(dR = det(R); log(dR))
            end
            out_expr = Expr(:tuple, :($logdet_expr - $(T(0.5)) * (Rδ' * Rδ)))
        else
            out_expr = Expr(:tuple, :($(T(-0.5)) * (Rδ' * Rδ)))
        end
        if track_y | track_μ
            push!(q.args, :(∂Rδ_∂δ = R' * Rδ))
            push!(q.args, :(∂Rδ_∂δ = ∂Rδ_∂δ + ∂Rδ_∂δ))
            track_y && push!(out_expr.args, :(-0.5∂Rδ_∂δ))
            track_μ && push!(out_expr.args, :(0.5∂Rδ_∂δ))
        end
        if track_Σ
            push!(q.args, :(∂Rδ_∂R = Rδ * δ'))
            push!(q.args, :(∂Rδ_∂R = ∂Rδ_∂R + ∂Rδ_∂R))
            push!(q.args, quote
                ∂out∂R = MutableFixedSizeVector{$L,$T}(undef)
                @vectorize $T for m ∈ 1:$M
                    ∂out∂R[m] = 1/dR[m] - 0.5 * ∂Rδ_∂R[m]
                end
                @vectorize $T for l ∈ 1:$(L-M)
                    ∂out∂R[l+$M] = - 0.5 * ∂Rδ_∂R[l + $M]
                end
            end)
            if lkjinv
                push!(out_expr.args, :(∂out∂R))
            else
                push!(out_expr.args, :(∂out∂R*∂R∂L))
            end
        end
        push!(q.args, out_expr)
        return q
    else
        if track_Σ
            if (Σ_type <: AbstractFloat) && ((y_type <: AbstractVector) || (μ_type <: AbstractVector))
                if y_type <: AbstractVector
                    logdet_expr = :($(M) * LinearAlgebra.logdet(R))
                else # μ_type <: AbstractVector
                    logdet_expr = :($(M) * LinearAlgebra.logdet(R))
                end
            else
                logdet_expr = :(logdet(R))
            end
            out_expr = :($logdet_expr - $(T(0.5)) * (Rδ' * Rδ))
        else
            out_expr = :($(T(-0.5)) * (Rδ' * Rδ))
        end
        quote
            $(Expr(:meta, :inline))
            δ = y - μ
            R = PaddedMatrices.invchol(Σ)
            Rδ = R*δ
            $out_expr
        end
    end
end
@generated function Normal(y::AbstractFixedSizeVector{M,T},
                                    μ, L::StructuredMatrices.AbstractLowerTriangularMatrix{M,T,LL},
                                    ::Val{track}) where {M, T <: Real, LL, track}
    multivariate_normal_lkj_quote(M, LL, T, track, false, false)
end
@generated function Normal(y::AbstractFixedSizeVector{M,T},
                                    μ, R::StructuredMatrices.AbstractUpperTriangularMatrix{M,T,LL},
                                    ::Val{track}) where {M, T <: Real, LL, track}

    multivariate_normal_lkj_quote(M, LL, T, track, false, true)
end

@inline function Normalc(
    sp::StackPointer,
    Y, μ,
    Σ::AbstractFixedSizeCovarianceMatrix{R,T,R},
    ::Val{track}
) where {R, K, T, track}
    Σcopy = DistributionParameters.MutableFixedSizeCovarianceMatrix{R,T,R}(undef)
    copyto!(Σcopy, Σ)
    Normal(sp, Y, μ, Σcopy, Val{track}())
end
@inline function Normalc(
    sp::StackPointer,
    Y, μ,
    Σ::AbstractFixedSizeCovarianceMatrix{R,T,R},
    Σcopy::AbstractFixedSizeCovarianceMatrix{R,T,R},
    ::Val{track}
) where {R, K, T, track}
    copyto!(Σcopy, Σ)
    Normal(sp, Y, μ, Σcopy, Val{track}())
end
@inline function ∂Normalc(
    sp::StackPointer,
    Y, μ,
    Σ::AbstractFixedSizeCovarianceMatrix{R,T,R},
    ::Val{track}
) where {R, K, T, track}
    Σcopy = DistributionParameters.MutableFixedSizeCovarianceMatrix{R,T,R}(undef)
    copyto!(Σcopy, Σ)
    ∂Normal(sp, Y, μ, Σcopy, Val{track}())
end
@inline function ∂Normalc(
    sp::StackPointer,
    Y, μ,
    Σ::AbstractFixedSizeCovarianceMatrix{R,T,R},
    Σcopy::AbstractFixedSizeCovarianceMatrix{R,T,R},
    ::Val{track}
) where {R, K, T, track}
    copyto!(Σcopy, Σ)
    ∂Normal(sp, Y, μ, Σcopy, Val{track}())
end

@generated function ∂Normal(y::AbstractFixedSizeVector{M,T},
                                    μ, L::StructuredMatrices.AbstractLowerTriangularMatrix{M,T,LL},
                                    ::Val{track}) where {M, T <: Real, LL, track}
    multivariate_normal_lkj_quote(M, LL, T, track, true, false)
end
@generated function ∂Normal(y::AbstractFixedSizeVector{M,T},
                                    μ, R::StructuredMatrices.AbstractUpperTriangularMatrix{M,T,LL},
                                    ::Val{track}) where {M, T <: Real, LL, track}

    multivariate_normal_lkj_quote(M, LL, T, track, true, true)
end

function matrix_normal_ar_lkj_quote(M, N, T, (track_y, track_μ, track_Λ, track_L), partial)
    W, Wshift = VectorizationBase.pick_vector_width_shift(T)
    V = Vec{W,T}
    @assert track_y == false
    initialize_block = quote
        Λᵥ = StructuredMatrices.cache(Λ)
        Ny = length(Y)
        qf = vbroadcast(Vec{$W,$T}, zero($T))
        δ = MutableFixedSizeMatrix{$M,$N,$V}(undef)
        δU = MutableFixedSizeMatrix{$M,$N,$V}(undef)
        Yᵥ = vectorizable(Y)
        i = 0
    end
    loop_block = quote
        Yblock = vload($V, Yᵥ + i)
        PaddedMatrices.diff!(δ, μ, Yblock)
        mul!(δU, δ, U)
    end
    closing_block = quote end
    return_expr = Expr(:tuple, :qfscalar)
    if partial
        # First, we look at initializations
        if track_μ
            push!(initialize_block.args, :(∂qf∂δ = zero(PaddedMatrices.MutableFixedSizeMatrix{$M,$N,Vec{$W,$T}}) ))
            push!(return_expr.args, :( vsum(∂qf∂δ) ) )
        end
        if ( track_μ || track_L )
            push!(initialize_block.args, :(Λᵥ′ΛᵥδU = PaddedMatrices.MutableFixedSizeMatrix{$M,$N,Vec{$W,$T}}(undef) ))
        end
        if track_Λ
            push!(initialize_block.args, :(∂qf∂Λ = vbroadcast(Vec{$W,$T}, zero($T))))
            push!(closing_block.args, :((logdetΛ, ∂logdetΛ) = StructuredMatrices.∂logdet(Λ)))
            push!(return_expr.args, :(vsum(∂qf∂Λ) + $N * Ny * ∂logdetΛ))
        end
        if track_L
            push!(initialize_block.args, :((U, ∂U∂L) = StructuredMatrices.∂inv(L)))
            push!(closing_block.args, :((logdetL, ∂logdetL) = StructuredMatrices.∂logdet(L)))
            push!(initialize_block.args, :(∂qf∂U = zero(StructuredMatrices.MutableUpperTriangularMatrix{N,Vec{$W,$T},$(StructuredMatrices.binomial2(N+1))})))
            push!(return_expr.args, :∂qf∂L)
        else
            push!(initialize_block.args, :(U = inv(L)) )
        end

        # Now, at the loop block
        if (track_L || track_μ)
            if track_Λ
                push!(loop_block.args, quote
                    qfᵢ, ∂qf∂Λᵢ = StructuredMatrices.∂selfcrossmul_and_quadform!(Λᵥ′ΛᵥδU, Λᵥ, δU)
                    qf = vadd(qfᵢ, qf)
                    ∂qf∂Λ = vadd(∂qf∂Λᵢ, ∂qf∂Λ)
                end)
            else
                push!(loop_block.args, quote
                    qf = vadd(StructuredMatrices.selfcrossmul_and_quadform!(Λᵥ′ΛᵥδU, Λᵥ, δU), qf)
                end)
            end
        elseif track_Λ
            push!(loop_block.args, quote
                qfᵢ, ∂qf∂Λᵢ = StructuredMatrices.∂quadform(Λᵥ, δU)
                qf = vadd(qfᵢ, qf)
                ∂qf∂Λ = vadd(∂qf∂Λᵢ, ∂qf∂Λ)
            end)
        else
            push!(loop_block.args, :(qf = vadd(StructuredMatrices.quadform(Λᵥ, δU), qf)))
        end
        track_L && push!(loop_block.args, :(StructuredMatrices.submul!(∂qf∂U, δ', Λᵥ′ΛᵥδU)))
        track_μ && push!(loop_block.args, :(StructuredMatrices.submul!(∂qf∂δ, Λᵥ′ΛᵥδU, U')))
    else # We are not taking partials
        track_L && push!(closing_block.args, :(logdetL = logdet(L)))
        track_Λ && push!(closing_block.args, :(logdetΛ = logdet(Λ)))

        push!(loop_block.args, :(qf = vadd(StructuredMatrices.quadform(Λᵥ, δU), qf)))
    end
    push!(loop_block.args, :(i += $W))
    # Here we handle the log determinants
    if track_L
        if track_Λ # track_L and track_Λ
            push!(closing_block.args, :(@fastmath qfscalar = Ny * ( $N * logdetΛ - $M * logdetL) - 0.5 * vsum(qf) ))
        else # track_L but not Λ
            push!(closing_block.args, :(@fastmath qfscalar = -Ny * $M * logdetL - 0.5 * vsum(qf) ))
        end
        if partial
            push!(closing_block.args, quote
                ∂qf∂L_part = StructuredMatrices.vsumvec(∂qf∂U)' * ∂U∂L
                ∂qf∂L = vmuladd($T(-Ny * $M), ∂logdetL', ∂qf∂L_part)
            end)
        end
    elseif track_Λ
        push!(closing_block.args, :(@fastmath qfscalar = Ny * $N * logdetΛ - 0.5 * vsum(qf) ))
    else
        push!(closing_block.args, :(@fastmath qfscalar = - 0.5 * vsum(qf) ))
    end
    quote
        $initialize_block
        Ysize = size(Y.data, 3)
        remmask = Y.mask
        @inbounds for ifrac ∈ 1:Ysize
            Yblock = vload($V, Yᵥ + i)
            PaddedMatrices.diff!(δ, μ, Yblock)
            if ifrac == Ysize
                PaddedMatrices.mask!(δ, remmask)
            end
        end
        $closing_block
        $(return_expression(return_expr))
    end
end

@generated function Normal(
            Y::AbstractScatteredArray{T,2,<: Union{SMatrix{M,N,T},AbstractFixedSizeMatrix{M,N,T}},1,2},
            μ::Union{<:SMatrix{M,N,T},<:AbstractFixedSizeMatrix{M,N,T}},
            Λ::AbstractAutoregressiveMatrix{T,V},
            L::LKJCorrCholesky{N,T}, ::Val{track}
        ) where {M,N,T,V,track}
    matrix_normal_ar_lkj_quote(M,N,T,track,false)
end
@generated function ∂Normal(
            Y::AbstractScatteredArray{T,2,<: Union{SMatrix{M,N,T},AbstractFixedSizeMatrix{M,N,T}},1,2},
            μ::Union{<:SMatrix{M,N,T},<:AbstractFixedSizeMatrix{M,N,T}},
            Λ::AbstractAutoregressiveMatrix{T,V},
            L::LKJCorrCholesky{N,T}, ::Val{track}
        ) where {M,N,T,V,track}
    matrix_normal_ar_lkj_quote(M,N,T,track,true)
end
function matrix_normal_ar_lkjinv_quote(M, N, T, (track_y, track_μ, track_Λ, track_U), partial)
    W, Wshift = VectorizationBase.pick_vector_width_shift(T)
    V = Vec{W,T}
    WT = VectorizationBase.REGISTER_SIZE
    @assert track_y == false

    δ = MutableFixedSizeMatrix{M,N,V,M,(M*N)}(undef)
    δU = MutableFixedSizeMatrix{M,N,V,M,(M*N)}(undef)
    # buffer_size = 2M*N*WT
    initialize_block = quote
        Λᵥ = StructuredMatrices.cache(Λ)
        Ny = length(Y)
        qf = vbroadcast(Vec{$W,$T}, zero($T))
        Yᵥ = vectorizable(Y)
        # δ  = MutableFixedSizeMatrix{$M,$N,$V,$M,$(M*N)}(undef)
        δ = $δ
        # δU = MutableFixedSizeMatrix{$M,$N,$V,$M,$(M*N)}(undef)
        δU = $δU
        i = 0
    end
    # partial ||
    # push!(initialize_block.args, quote
    #     println("Mean:")
    #     println(μ)
    #     println("AR matrix:")
    #     println(Λ.spacing.rinvOmρ²ᵗ)
    #     println(Λ.spacing.ρᵗ .* Λ.spacing.rinvOmρ²ᵗ)
    #     println("U:")
    #     println(U)
    # end)
    loop_block = quote
        mul!(δU, δ, U)
    end
    closing_block = quote end
    return_expr = Expr(:tuple, :qfscalar)
    if partial
        # First, we look at initializations
        if track_μ
            # push!(initialize_block.args, :(∂qf∂δ = zero(PaddedMatrices.MutableFixedSizeMatrix{$M,$N,Vec{$W,$T}}) ))
            ∂qf∂δ = PaddedMatrices.MutableFixedSizeMatrix{M,N,Vec{W,T}}(undef)
            push!(initialize_block.args, quote
                ∂qf∂δ = $∂qf∂δ
                PaddedMatrices.zero!(∂qf∂δ)
            end)
            push!(return_expr.args, :( vsum(∂qf∂δ) ) )
        end
        if ( track_μ || track_U )
            # push!(initialize_block.args, :(Λᵥ′ΛᵥδU = PaddedMatrices.MutableFixedSizeMatrix{$M,$N,Vec{$W,$T}}(undef) ))
            Λᵥ′ΛᵥδU = PaddedMatrices.MutableFixedSizeMatrix{M,N,Vec{W,T}}(undef)
            push!(initialize_block.args, :(Λᵥ′ΛᵥδU = $Λᵥ′ΛᵥδU ))
        end
        if track_Λ
            push!(initialize_block.args, :(∂qf∂Λ = vbroadcast(Vec{$W,$T}, zero($T))))
            push!(closing_block.args, :((logdetΛ, ∂logdetΛ) = StructuredMatrices.∂logdet(Λ)))
            push!(return_expr.args, :(vsum(∂qf∂Λ) + $N * Ny * ∂logdetΛ))
        end
        if track_U
            push!(closing_block.args, :((logdetU, ∂logdetU) = StructuredMatrices.∂logdet(U)))
            # push!(initialize_block.args, :(∂qf∂U = zero(StructuredMatrices.MutableUpperTriangularMatrix{$N,Vec{$W,$T},$(StructuredMatrices.binomial2(N+1))})))
            ∂qf∂U = StructuredMatrices.MutableUpperTriangularMatrix{N,Vec{W,T},StructuredMatrices.binomial2(N+1)}(undef)
            push!(initialize_block.args, quote
                ∂qf∂U = $∂qf∂U
                PaddedMatrices.zero!(∂qf∂U)
            end)
            push!(return_expr.args, :∂out∂U)
        end

        # Now, at the loop block
        if (track_U || track_μ)
            if track_Λ
                push!(loop_block.args, quote
                    qfᵢ, ∂qf∂Λᵢ = StructuredMatrices.∂selfcrossmul_and_quadform!(Λᵥ′ΛᵥδU, Λᵥ, δU)
                    qf = vadd(qfᵢ, qf)
                    ∂qf∂Λ = vadd(∂qf∂Λᵢ, ∂qf∂Λ)
                end)
            else
                push!(loop_block.args, quote
                    qf = vadd(StructuredMatrices.selfcrossmul_and_quadform!(Λᵥ′ΛᵥδU, Λᵥ, δU), qf)
                end)
            end
        elseif track_Λ
            push!(loop_block.args, quote
                qfᵢ, ∂qf∂Λᵢ = StructuredMatrices.∂quadform(Λᵥ, δU)
                qf = vadd(qfᵢ, qf)
                ∂qf∂Λ = vadd(∂qf∂Λᵢ, ∂qf∂Λ)
            end)
        else
            push!(loop_block.args, :(qf = vadd(StructuredMatrices.quadform(Λᵥ, δU), qf)))
        end
        track_U && push!(loop_block.args, :(StructuredMatrices.submul!(∂qf∂U, δ', Λᵥ′ΛᵥδU)))
        # if track_U
        #     push!(loop_block.args, quote
        #         StructuredMatrices.submul!(∂qf∂U, δ', Λᵥ′ΛᵥδU)
        #         if !StructuredMatrices.all_finite(∂qf∂U)
        #             println("Not all elements of ∂qf∂U were finite!\n")
        #             println(∂qf∂U)
        #             println("\nArguments, δ':\n")
        #             println(δ')
        #             println("\nArguments, Λᵥ′ΛᵥδU:\n")
        #             println(Λᵥ′ΛᵥδU)
        #             throw("\nAll elements of ∂qf∂U should be finite!")
        #         end
        #     end)
        # end
        track_μ && push!(loop_block.args, :(StructuredMatrices.submul!(∂qf∂δ, Λᵥ′ΛᵥδU, U')))
    else # We are not taking partials
        track_U && push!(closing_block.args, :(logdetU = logdet(U)))
        track_Λ && push!(closing_block.args, :(logdetΛ = logdet(Λ)))

        push!(loop_block.args, :(qf = vadd(StructuredMatrices.quadform(Λᵥ, δU), qf)))
    end
    push!(loop_block.args, :(i += $W))
    # Here we handle the log determinants
    if track_U
        if track_Λ # track_U and track_Λ
            # if !partial
            #     push!(closing_block.args, :(@show Ny, $M, logdetU, $N, logdetΛ))
            #     push!(closing_block.args, :(@show -0.5qfvsum))
            #     push!(closing_block.args, :(@show qf))
            #     push!(closing_block.args, :(@show  Ny * ( $M * logdetU + $N * logdetΛ)))
            # end
            push!(closing_block.args, :(@fastmath qfscalar = Ny * ( $M * logdetU + $N * logdetΛ) - 0.5qfvsum))
        else # track_U but not Λ
            push!(closing_block.args, :(@fastmath qfscalar = Ny * $M * logdetU - 0.5qfvsum ))
        end
        if partial
            push!(closing_block.args, quote
                ∂qf∂Usummed = StructuredMatrices.vsum(∂qf∂U)
                coef = $T(Ny * $M)
                @vectorize $T for n ∈ 1:$N
                    ∂qf∂Usummed[n] = coef * ∂logdetU[n] + ∂qf∂Usummed[n]
                end
                ∂out∂U = UpperTriangularMatrix(∂qf∂Usummed)
            end)
        end
    elseif track_Λ
        push!(closing_block.args, :(@fastmath qfscalar = Ny * $N * logdetΛ - 0.5qfvsum ))
    else
        push!(closing_block.args, :(qfscalar = - 0.5qfvsum ))
    end
    quote
        # @show $((track_y, track_μ, track_Λ, track_U))
        $initialize_block
        Ysize = size(Y.data, 3)
        remmask = Y.mask
        @inbounds for ifrac ∈ 1:Ysize
            # Yblock = vload($V, Yᵥ + i)
            # PaddedMatrices.diff!(δ, μ, Yblock)
            # PaddedMatrices.vload!(δ, Yᵥ + i)
            Yᵢ = PaddedMatrices.vload($V, Yᵥ + i)
            PaddedMatrices.diff!(δ, μ, Yᵢ)
            ifrac == Ysize && PaddedMatrices.mask!(δ, remmask)
            $loop_block
        end
        qfvsum = vsum(qf)
        $closing_block
        $(return_expression(return_expr))
    end
end
@generated function Normal(
            Y::AbstractScatteredArray{T,2,<: Union{SMatrix{M,N,T},AbstractFixedSizeMatrix{M,N,T}},1},
            μ::Union{<:SMatrix{M,N,T},<:AbstractFixedSizeMatrix{M,N,T}},
            Λ::AbstractAutoregressiveMatrix{T,V},
            U::StructuredMatrices.AbstractUpperTriangularMatrix{N,T}, ::Val{track}
        # ) where {M,N,V,T,track}
        ) where {M,N,T,V,track}
    matrix_normal_ar_lkjinv_quote(M,N,T,track,false)
end
@generated function ∂Normal(
            Y::AbstractScatteredArray{T,2,<: Union{SMatrix{M,N,T},AbstractFixedSizeMatrix{M,N,T}},1},
            μ::Union{<:SMatrix{M,N,T},<:AbstractFixedSizeMatrix{M,N,T}},
            Λ::AbstractAutoregressiveMatrix{T,V},
            U::StructuredMatrices.AbstractUpperTriangularMatrix{N,T}, ::Val{track}
        # ) where {M,N,V,T,track}
        ) where {M,N,T,V,track}
    matrix_normal_ar_lkjinv_quote(M,N,T,track,true)
end


text:
module Epochs

using LeapSeconds: offset_tai_utc, offset_utc_tai
using EarthOrientation: getΔUT1_TAI, EOP_DATA, get
using LeapSeconds: offset_tai_utc, offset_utc_tai
using MuladdMacro

import Dates
import Dates: format
import Dates: year, day, month, hour, minute, second
import Dates: millisecond, microsecond, nanosecond
import Dates: yearmonthday, dayofyear, now

import LeapSeconds

import ..EPOCH_ISO_FORMAT

import ..AstroDates: Date, DateTime, Time
import ..AstroDates: calendar, fractionofday, fractionofsecond, secondinday, subsecond
import ..AstroDates: j2000, julian, julian_twopart
import ..AccurateArithmetic: apply_offset, two_sum

using ..TimeScales: find_path

export CCSDS_EPOCH, FIFTIES_EPOCH, FUTURE_INFINITY, GALILEO_EPOCH, GPS_EPOCH, J2000_EPOCH
export J2000_TO_JULIAN, J2000_TO_MJD, JULIAN_EPOCH, MODIFIED_JULIAN_EPOCH, PAST_INFINITY
export UNIX_EPOCH
export Epoch
export NoOffsetError, NoPathError
export year, month, day, hour, minute, second
export millisecond, microsecond, nanosecond, subsecond
export yearmonthday, dayofyear, fractionofday, fractionofsecond, secondinday
export timescale, now
export getoffset, insideleap
export j2000, julian, julian_period, julian_twopart, modified_julian
export from_utc, to_utc
export -

using ..TimeScales
using ..AstroDates
using ..Periods

include("types.jl")
include("julian.jl")
include("tdb_constants.jl")
include("offsets.jl")
include("accessors.jl")
include("operations.jl")
include("dates.jl")
include("aliases.jl")
include("ranges.jl")
include("epoch_constants.jl")
include("io.jl")
include("utc.jl")

end

max_stars_repo_name: JuliaAstrodynamics/AstronomicalTime

max_issues_repo_name: JuliaAstrodynamics/AstronomicalTime

max_forks_repo_name: JuliaAstrodynamics/AstronomicalTime


text:
using Documenter, EspyInsideFunction

push!(LOAD_PATH,"../src")
makedocs(sitename ="EspyInsideFunction.jl",
        modules   = [EspyInsideFunction],
        format    = Documenter.HTML(prettyurls = false),
        pages     = ["index.md"],
        source  = "src",
        build   = "build"
        )

here = @__DIR__
mv(here*"\\build\\index.html",here*"\\index.html", force=true)
mv(here*"\\build\\search.html",here*"\\search.html", force=true)
mv(here*"\\build\\search_index.js",here*"\\search_index.js", force=true)
mv(here*"\\build\\assets",here*"\\assets", force=true)
#deploydocs(repo = "github.com/PhilippeMaincon/EspyInsideFunction.jl.git")

max_stars_repo_name: PhilippeMaincon/Espy


text:
function test_zygote(f, xs...; kws...)
    config = ZygoteRuleConfig()
    test_rrule(config, f, xs...; kws..., rrule_f = rrule_via_ad)
end

function rr(seed)
    rng = Random.seed!(seed)
    return rng
end

max_stars_repo_name: JuliaML/MLBase

max_issues_repo_name: JuliaML/MLBase

max_forks_repo_name: JuliaML/MLBase


text:
include("taxi_models_and_data.jl")

ld = 2
batch_size = [128, 256, 1024]
epochs = [500, 750, 1000]
learning_rates = [1f-4, 7f-4, 2f-3]
losses = [DCGANLoss(), LSLoss(), WLossGP()]
losses_name = ["BCE", "LS", "WLoss"]
dirs = []

for (loss, loss_name) in zip(losses, losses_name)
	for (bs, nepoch) in zip(batch_size, epochs)
		for lr in learning_rates
			output_dir = "$(loss_name)_BS$(bs)_LR$(lr)"
			push!(dirs, output_dir)
			println("starting $(output_dir)")
			G, D, Ghist, Dhist = train(Settings(G=TaxiGConv,
								  D=TaxiDConvSpectral,
								  epochs=nepoch,
								  batch_size=bs,
								  rand_input=Taxi_input,
								  loss=loss,
								  img_fun=gen_Taxi_images,
								  nclasses=2,
								  latent_dim=ld,
								  verbose_freq=10,
								  optD = ADAM(lr, (0.5, 0.99)),
								  optG = ADAM(lr, (0.5, 0.99)),
								  output_dir = "$(output_dir)",))

			G = TaxiGConv(G.g_labels |> cpu, G.g_latent |> cpu, G.g_common |> cpu)
			BSON.@save "$(output_dir)/conv_generator_ld$(ld).bson" G

			D = TaxiDConvSpectral(D.d_labels |> cpu, D.d_common |> cpu)
			BSON.@save "$(output_dir)/conv_discriminator_ld$(ld).bson" D

			BSON.@save "$(output_dir)/d_hist_ld$(ld).bson" Dhist
			BSON.@save "$(output_dir)/g_hist_ld$(ld).bson" Ghist
			BSON.@save "dirs.bson" dirs
		end
	end
end

max_stars_repo_name: sisl/VerifyGAN

max_issues_repo_name: sisl/VerifyGAN

max_forks_repo_name: sisl/VerifyGAN


text:
using DrWatson
quickactivate(@__DIR__)
@show projectdir()
##
using StatsBase
using StatsPlots
using DataFrames
using Dates
using CSV
plotlyjs()
##
fname = projectdir("data/smoothed_contact_rates.csv")
surv = CSV.File(fname) |> DataFrame
##
fname = projectdir("data/google/mobility_denmark.csv")
mobil = CSV.File(fname) |> DataFrame
##
fname = "/home/and/data/covidsurvey/contacts.csv"
rt = CSV.File(fname) |> DataFrame
fname = projectdir("reports/2021-07-21/Rt_cases=_chains=1_hospit=false_model=v1_sero=false_steps=1000_warmup=1000.csv")
fname = projectdir("reports/2021-07-21/Rt_cases=2020-06-01_chains=5_hospit=false_model=hospit_sero=false_steps=1000_warmup=1000.csv")
rt = CSV.File(fname) |> DataFrame
fname = projectdir("reports/2021-07-21/Rt_cases=_chains=1_hospit=false_model=hospit_sero=false_steps=1000_warmup=1000.csv")
rt_d = CSV.File(fname) |> DataFrame
fname = projectdir("reports/2021-07-21/Rt_cases=2020-06-01_chains=5_hospit=false_model=cases_sero=false_steps=1000_warmup=1000.csv")
fname = projectdir("reports/national_prediction_end_of_peak/Rt_chains=5_model=hospit_observ=2021-01-13_predstart=2020-11-10_regional=false_semipar=true_steps=200_warmup=1000.csv")
rt = CSV.File(fname) |> DataFrame
rename!(rt, :dates=>:date)
##
function StatsBase.zscore(df::DataFrame)
    y = similar(df)
    for n in names(df)
        y[!,n] = zscore(df[:,n])
    end
    return y
end

function conform(args...; normalize = false)
    s = maximum( first.( select.(args, Ref(:date) ) ) ).date
    e = minimum( last.( select.(args, Ref(:date) ) ) ).date
    x = s:Day(1):e
    y = filter.(Ref(:date => >=(s)), args)
    y = filter.(Ref(:date => <=(e)), y)
    normalize && ( y = zscore.( select.( y, Ref(Not(:date)) ) ) )
    (x, y...)
end

let
    x,yr,ys,ym = conform(rt[:,[:mq,:date]], surv, mobil[:,[:nonresidential, :date]], normalize=true)#rt_ssi[:,[:mq,:date]])

    p = plot(xticks=:native, size=(1300,800), lw=2)
    plot!(x, yr.mq, label="Rt (hospit)", lw=2)
    # plot!(x, yssi.mq, label="Rt (SSI)", lw=2)
    plot!(x, ys.strangers, lab="contacts to strangers", lw=2)
    plot!(x, ym.nonresidential, lab="nonresidential mobility", lw=2)
    # plot!(x, y2.friends, lab="friends", lw=2)
    # plot!(x, y2.family, lab="family", lw=2)
    # plot!(x, y2.colleagues, lab="colleague", lw=2)

    fname = "/home/and/.tmp/tmp.html"
    savefig(p, fname)
    run(`firefox $(fname)`, wait=false)
end

## ============================================================================
plotlyjs()
let
    obs = Dict()
    National.observables!(obs)

    p = plot(xticks=:native, size=(1000,800))
    x, y = eachcol(obs["cases"])
    plot!(x,y|>zscore,lab="cases")
    x, y = eachcol(obs["hospit"])
    plot!(x,y|>zscore,lab="hospit")
    x, y = eachcol(obs["deaths"])
    plot!(x,y|>zscore,lab="deaths")

    fname = "/home/and/.tmp/tmp.html"
    savefig(p, fname)
    run(`firefox $(fname)`, wait=false)
    p
end

ys = []
for n in ["cases", "hospit", "deaths"]
    x, y = eachcol(obs[n])
    i = findall(x-> Date("2020-11-30")<x<Date("2021-02-01"), x)
    push!(ys, y[i]|>zscore)
end

lags = 0:50
c = crosscor(ys[2], ys[3], lags)
plot(lags, c)
argmax(c)

## =============================================================================
plot(NegativeBinomial2(1., 5.41))
plot!(NegativeBinomial2(4.5, 5.41))
plot!(GammaMeanCv(5.1, 0.86))
mean(GammaMeanCv(5.1, 0.86))

## =============================================================================
# compare with SSI
p = let
    fname = projectdir("reports/2021-07-21/Rt_cases=2020-06-01_chains=5_hospit=false_model=cases_sero=false_steps=1000_warmup=1000.csv")
    rt = CSV.File(fname) |> DataFrame
    rt_ssi = projectdir("data/Rt_SSI.csv") |> load |> DataFrame
    date, yrt, yrt_ssi = conform(rt[:,[:mq,:date]], rt_ssi[:,[:mq,:date]]; normalize = false)
    p = plot(legend=:topleft, ylabel="R_t", size=(800, 500), xticks=:native)
    plot!(date, yrt_ssi.mq, lw=2, label="SSI", c=:sienna)
    plot!(date, yrt.mq, lw=2, label="model without predictors", c=:royalblue)
end

## =============================================================================
# sero
data = National.load_data()
plot!(data.cases.date  .- Day(16), cumsum(data.cases.country) / data.turing_data.population*2)
scatter!(data.seroprev.date , data.seroprev.country)

max_stars_repo_name: andreaskoher/CovidSurvey

max_issues_repo_name: andreaskoher/CovidSurvey

max_forks_repo_name: andreaskoher/CovidSurvey


text:
import NodeModule.Node

function create_distance_matrix(nodes::Array{Node})
    """Creates and returns the matrix of distances between given nodes.
    Uses euclidean_distance function to calculate distance.
    """
    node_num = length(nodes)
    result_matrix = zeros(node_num, node_num)

    for n1 in nodes
        for n2 in nodes
            result_matrix[n1.id,n2.id] = euclidean_distance(n1, n2)
        end
    end
    return result_matrix
end

max_stars_repo_name: zahidkizmaz/firefly-algorithm-tsp

max_issues_repo_name: zahidkizmaz/firefly-algorithm-tsp

max_forks_repo_name: zahidkizmaz/firefly-algorithm-tsp


text:
my_f(2,1)
my_f(2,3)

derivate_of_my_f(2,1)

@testset "Sonic.jl" begin
    @test my_f(2, 1) == 7
    @test my_f(2, 3) == 31
end

@testset "Derivative Tests" begin
    @test derivate_of_my_f(2, 1) == 2
end


text:
# svrg implementation

export SVRG

"""
# struct SVRG

The structure describing iterations for stochastic variance reduced gradient.

## Method

We have the following sequence:
   - a batch_gradient as the full gradient at current position
   - storing gradient and position before the mini batch sequence
   - then for `nb_mini_batch` :
       - uniform sampling of **one** term of the summation
       - computation of the gradient of the term at current position and the gradient
                        at position before mini batch
       - computation of direction of propagation as the batch gradient + gradient of term at current
                        position -  gradient of term at position before mini batch sequence
       - update of position with adequate step size

   The step size used in the algorithm is constant and according to the ref paper it should be of the order of
   L/4 where L is the lipschitz constant of the function to minimize

## FIELDS

- nb\\_minibatch : The number of minibatch to run after one full batch run
- stepsize: the step to use in
"""
mutable struct SVRG
    # number of minibatch of one term.
    nb_minibatch :: Int64
    # step_size
    stepsize :: Float64
end

"""
# function minimize(svrgpb::SVRG, evaluation::Evaluator{F,G}, max_iter::Int64, initialposition::Array{Float64,N})

Generic function used in minimisations for all algorithms.

The function has 3 type parameters F, G and N:

- F and G corresponds to the types of function used in instantiating the structure TermFunction{F} and TermmGradient{G}
(see in test and examples directories the cases TestLinearRegression or TestLogisticRegression).

- N is the dimension of array for searched parameters/gradients which depend upon the problem and the way we modelize it.
(See Logistic regression where we used N=2)

## Args
- svrgpb : the structure describing main parameters of the batch strategy
- Evaluator{F,G}
- max_iter : maximum number of iterations
- initialposition : initial position at iteration beginning

"""
function minimize(svrgpb::SVRG, evaluation::Evaluator{F,G}, max_iterations, initial_position::Array{Float64,N}) where {F,G,N}
    direction = zeros(Float64, size(initial_position))
    full_batch_gradient = zeros(Float64, size(initial_position))
    mini_batch_gradient_current = zeros(Float64, size(initial_position))
    mini_batch_gradient_origin = zeros(Float64, size(initial_position))
    nbterms = get_nbterms(evaluation)
    #
    mt = MersenneTwister(117);
    position = Array{Float64,N}(initial_position)
    iteration = 0
    initial_value = 0.
    more = true
    while more
        iteration += 1
        # get iteration parameters
        # compute gradient on large batch index set and store initial position
        compute_gradient!(evaluation, position ,full_batch_gradient)
        # get number of mini batch
        nb_mini_batch = svrgpb.nb_minibatch
        position_before_mini_batch = Array{Float64,N}(position)
        # loop on small batch iterations
        for i in 1:nb_mini_batch
            # sample mini batch terms
            term = 1 + floor(Int64, nbterms * rand(mt))
            compute_gradient!(evaluation, position , term, mini_batch_gradient_current)
            compute_gradient!(evaluation, position_before_mini_batch , term, mini_batch_gradient_origin)
            direction = mini_batch_gradient_current - mini_batch_gradient_origin + full_batch_gradient;
            position = position - svrgpb.stepsize * direction;
        end
        value = compute_value(evaluation, position)
        if iteration == 1
            initial_value = value
        end
        if iteration >= max_iterations
            @info "Reached maximal number of iterations required, initial_value final value:" initial_value value
            return position, value
        end
    end

end


text:
using DSGE, JLD2, DataStructures
import Test: @test, @testset
import DSGE: @test_matrix_approx_eq

# Load true values
true_inverse_basis = load("../../test_outputs/solve/solve.jld2", "inverse_basis")
G1 = load("../../test_outputs/solve/solve.jld2", "G1")
impact = load("../../test_outputs/solve/solve.jld2", "impact")

# Test
m = KrusellSmith()
### Steady state

# Load test values
params = load("../../test_outputs/krusell_smith/saved_params.jld2", "params")
grids = load("../../test_outputs/krusell_smith/saved_params.jld2", "grids")
init_params = load("../../test_outputs/krusell_smith/saved_params.jld2", "init_params")
approx_params = load("../../test_outputs/krusell_smith/saved_params.jld2", "approx_params")
vars_ss_mat = load("../../test_outputs/krusell_smith/vars_ss.jld2")
vars_ss_mat = vars_ss_mat["varsSS"]

# Update parameters
m[:γ] = params[:γ]
m[:ρ] = params[:ρ]
m[:δ] = params[:δ]
m[:α] = params[:α]
m[:σ_tfp] = params[:σ_tfp]
m[:ρ_tfp] = params[:ρ_tfp]
m[:λ1] = params[:λ1]
m[:λ2] = params[:λ2]
m[:μ] = params[:μ]
m[:τ] = params[:τ]

# Update settings
update!(m.settings[:I], Setting(:I, 100))
update!(m.settings[:amin], Setting(:amin, 0.))
update!(m.settings[:amax], Setting(:amax, 100.))
update!(m.settings[:a], Setting(:a, collect(linspace(0., 100., 100))))
update!(m.settings[:da], Setting(:da, 100./(100 - 1)))
update!(m.settings[:aaa], Setting(:aaa, reshape(get_setting(m, :aa), 2 * 100, 1)))
update!(m.settings[:J], Setting(:J, 2))
update!(m.settings[:z], Setting(:z, [0., 1.]))
update!(m.settings[:dz], Setting(:dz, 1.))
update!(m.settings[:zz], Setting(:zz, ones(100, 1) * get_setting(m, :z)'))
update!(m.settings[:zzz], Setting(:zzz, reshape(get_setting(m, :zz), 2 * 100, 1)))
update!(m.settings[:n_jump_vars], Setting(:n_jump_vars, length(m.endogenous_states[:value_function])))
update!(m.settings[:n_state_vars], Setting(:n_state_vars, length(m.endogenous_states[:distribution]) + n_shocks_exogenous(m)))
update!(m.settings[:n_state_vars_unreduce], Setting(:n_state_vars_unreduce, 0))
update!(m.settings[:rmax], Setting(:rmin, .0001))
update!(m.settings[:rmax], Setting(:rmax, m[:ρ].value))
update!(m.settings[:r0], Setting(:r0, .005))
update!(m.settings[:maxit_HJB], Setting(:maxit_HJB, 100))
update!(m.settings[:crit_HJB], Setting(:crit_HJB, 1e-6))
update!(m.settings[:Δ_HJB], Setting(:Δ_HJB, 1e4))
update!(m.settings[:Ir], Setting(:Ir, 100))
update!(m.settings[:crit_S], Setting(:crit_S, 1e-5))
update!(m.settings[:reduce_state_vars], Setting(:reduce_state_vars, true))
update!(m.settings[:reduce_v], Setting(:reduce_v, true))
update!(m.settings[:krylov_dim], Setting(:krylov_dim, 5))
update!(m.settings[:n_knots], Setting(:n_knots, 12))
update!(m.settings[:c_power], Setting(:c_power, 7))
update!(m.settings[:knots_dict], Setting(:knots_dict, Dict(1 => collect(linspace(get_setting(m, :amin), get_setting(m, :amax), 12 - 1)))))
update!(m.settings[:spline_grid], Setting(:spline_grid, get_setting(m, :a)))
update!(m.settings[:n_prior], Setting(:n_prior, 1))
update!(m.settings[:n_post], Setting(:n_post, 2))

# Perform solve
steadystate!(m)
TTT, RRR, ~, inverse_basis, ~ = solve(m; sparse_mat = false)
@testset "Simulating Impulse Responses (non-sparse)" begin
    @test @test_matrix_approx_eq full(TTT) full(G1)
    @test @test_matrix_approx_eq full(RRR) full(impact)
    @test @test_matrix_approx_eq full(true_inverse_basis) full(inverse_basis)
end

steadystate!(m)
TTT, RRR, ~, inverse_basis, ~ = solve(m; sparse_mat = true)
@testset "Simulating Impulse Responses (sparse)" begin
    @test @test_matrix_approx_eq full(TTT) full(G1)
    @test @test_matrix_approx_eq full(RRR) full(impact)
    @test @test_matrix_approx_eq full(true_inverse_basis) full(inverse_basis)
end


text:
#function FindNeighbourObjs(universe::Universe, obj::DefectObj)
#    universe.objs
#end

function FindNeighbourObjs(universe::Universe, obj::DefectObj)
    objs = Vector{DefectObj}()
    for dCellIndex in universe.cellContainer.dCellIndexs
        cellIndex = obj.cellIndex + dCellIndex
        try
            objs =vcat(objs, universe.cellContainer.cells[CartesianIndex(Tuple(cellIndex))].objs)
        catch
            continue
        end
        #objs = vcat(objs, theObjs) #need modify
    end
    objs
end

function DCellIndexs(nDimension)
    dCellIndexs = Vector{Vector{Int64}}()
    dCellIndex = Vector{Int64}()
    IterDCellIndexs(dCellIndexs, dCellIndex, 1, nDimension)
    dCellIndexs
end

function IterDCellIndexs(dCellIndexs::Vector{Vector{Int64}}, dCellIndex::Vector{Int64},
                             d::Int64, nDimension::Int64)
    for dI in [-1,0,1]
        if d < nDimension
            IterDCellIndexs(dCellIndexs,vcat(dCellIndex,dI), d+1, nDimension)
        else
            push!(dCellIndexs, copy(vcat(dCellIndex,dI)))
        end
    end
end

max_stars_repo_name: permissionx/mc-ke

max_issues_repo_name: permissionx/mc-ke

max_forks_repo_name: permissionx/mc-ke


text:
# transmcp.jl
# translated by Changhyun Kwon from transmcp.gms
# http://www.gams.com/modlib/libhtml/transmcp.htm
#
# Transportation model as equilibrium problem (TRANSMCP,SEQ=126)
#
#    Dantzig's original transportation model (TRNSPORT) is
#    reformulated as a linear complementarity problem.  We first
#    solve the model with fixed demand and supply quantities, and
#    then we incorporate price-responsiveness on both sides of the
#    market.
#
#
# Dantzig, G B, Chapter 3.3. In Linear Programming and Extensions.
# Princeton University Press, Princeton, New Jersey, 1963.

using Complementarity, JuMP

plants = ["seattle", "san-diego"]
markets = ["new-york", "chicago", "topeka"]

capacity = [350, 600]
a = Dict(zip(plants, capacity))

demand = [325, 300, 275]
b = Dict(zip(markets, demand))

elasticity = [1.5, 1.2, 2.0]
esub = Dict(zip(markets, elasticity))

distance = [ 2.5 1.7 1.8 ;
             2.5 1.8 1.4  ]
d = Dict()
for i in 1:length(plants), j in 1:length(markets)
    d[plants[i], markets[j]] = distance[i,j]
end

f = 90

m = MCPModel()
@variable(m, w[i in plants] >= 0)
@variable(m, p[j in markets] >= 0)
@variable(m, x[i in plants, j in markets] >= 0)

@NLexpression(m, c[i in plants, j in markets], f * d[i,j] / 1000)

@NLexpression(m, profit[i in plants, j in markets],    w[i] + c[i,j] - p[j])
@NLexpression(m, supply[i in plants],                  a[i] - sum{x[i,j], j in markets})
@NLexpression(m, fxdemand[j in markets],               sum{x[i,j], i in plants} - b[j])

complements(m, profit, x)
complements(m, supply, w)
complements(m, fxdemand, p)

PATHSolver.path_options(
                "convergence_tolerance 1e-8",
                "output yes",
                "time_limit 3600"
                )

status = solveMCP(m)

@show getvalue(x)
@show getvalue(w)
@show getvalue(p)

@show status
@assert status == :Solved
@assert getvalue(x["seattle", "chicago"]) == 300.0
@assert getvalue(p["topeka"]) == 0.126