text:
export
    newcontext, freecontext, # newcontext_puts,
    # checkoutlicense, releaselicense, get_release,
    load_param_file, save_param_file, load_tuner_file,
    set_param, get_param, # reset_params_to_defaults,
    # set_feastols, set_names, set_linearvars, set_honorbnds,
    # addcompcons, chgvarbnds,
    init_problem, solve_problem, restart_problem,
    mip_init_problem, mip_solve_problem,
    #set_branching_priorities, set_findiff_relstepsizes
    get_number_FC_evals, get_number_GA_evals,
    get_number_H_evals, get_number_HV_evals,
    get_number_iters, get_number_cg_iters,
    get_abs_feas_error, get_rel_feas_error,
    get_abs_opt_error, get_rel_opt_error,
    # get_solution, get_constraint_values,
    # get_objgrad_values, get_jacobian_values, get_hessian_values,
    get_mip_num_nodes, get_mip_num_solves,
    get_mip_abs_gap, get_mip_rel_gap,
    get_mip_incumbent_obj, get_mip_relaxation_bnd,
    get_mip_lastnode_obj, get_mip_incumbent_x #, check_first_ders,

"""
Returns a pointer to the solver object that is used in all other KNITRO API calls.

A new KNITRO license is acquired and held until `freecontext` has been called,
or until the calling program ends.
"""
function newcontext()
    ptr = @ktr_ccall(new,Ptr{Nothing},())
    if ptr == C_NULL
        error("KNITRO: Error creating solver")
    end
    ptr
end

"""
Returns a pointer to the solver object that is used in all other KNITRO API calls.

A new KNITRO license is acquired and held until `freecontext` has been called,
or until the calling program ends.

This function also takes an argument `f` that sets a 'put string' callback
function to handle output generated by the KNITRO solver, and an (optional) pointer
`userParams` for passing user-defined data.
"""
function newcontext_puts(f::Function, userParams=C_NULL)
    cb = cfunction(f, Cint, (Ptr{Cchar}, Ptr{Nothing}))
    ptr = @ktr_ccall(new,Ptr{Nothing}, (Ptr{Nothing}, Ptr{Nothing}),
                     cb, userParams)
    if ptr == C_NULL
        error("KNITRO: Error creating solver with put-string")
    end
    ptr
end

"Free all memory and release any KNITRO license acquired with `kp_env`"
function freecontext(kp_env::Ptr{Nothing})
    if kp_env != C_NULL
        println("KNITRO: calling freecontext on $(kp_env)")
        println("C_NULL (for reference): $(C_NULL)")
        return_code = @ktr_ccall(free, Int32, (Ptr{Nothing},), kp_env)
        if return_code != 0
            error("KNITRO: Error freeing memory")
        end
        println("KNITRO: freecontext successful")
    end
end

# /** Allocate memory for a license from the Artelys License Manager for high
#  *  volume KNITRO applications.  The license will be checked out the first
#  *  time KTR_new_zlm is called.  The license must be checked in later by
#  *  calling ZLM_release_license.
#  *  Returns NULL on error.
#  */
# ZLM_context_ptr  KNITRO_API ZLM_checkout_license (void);

# function checkoutlicense()
#   ptr = @zlm_ccall(checkout_license,Ptr{Nothing},())
#   if ptr == C_NULL
#     error("KNITRO: Error checking out license")
#   end
#   ZLMcontext(ptr)
# end

# /** Returns a pointer to the solver object that is used in all other KNITRO
#  *  API calls.  Pass the license acquired by calling ZLM_checkout_license.
#  *  This function also takes an argument that sets a "put string" callback
#  *  function to handle output generated by the KNITRO solver, and a pointer
#  *  for passing user-defined data.  See KTR_set_puts_callback for more
#  *  information.
#  *  Returns NULL on error.
#  */
# KTR_context_ptr  KNITRO_API KTR_new_zlm (KTR_puts    * const  fnPtr,
#                                          void        * const  userParams,
#                                          ZLM_context * const  pZLMcontext);

# /** Release the KNITRO license and free allocated memory.
#  *  KNITRO will set the address of the pointer to NULL after freeing
#  *  memory, to help avoid mistakes.
#  *  Returns 0 if OK, nonzero if error.
#  */
# int  KNITRO_API ZLM_release_license (ZLM_context *  pZLMcontext);

# function releaselicense(zc::ZLMcontext)
#   if zc.context != C_NULL
#     return_code = @ktr_ccall(free, Int32, (Ptr{Nothing},), zc.context)
#     if return_code != 0
#       error("KNITRO: Error releasing license and freeing memory")
#     end
#     zc.context = C_NULL
#   end
# end

"Reset all parameters to default values"
function reset_params_to_defaults(kp::KnitroProblem)
    return_code = @ktr_ccall(reset_params_to_defaults, Int32,
                             (Ptr{Nothing},), kp.env)
    if return_code != 0
        error("KNITRO: Error resetting parameters to default values")
    end
end

"Set all parameters specified in the given file"
function load_param_file(kp::KnitroProblem, filename::AbstractString)
    return_code = @ktr_ccall(load_param_file, Int32, (Ptr{Nothing}, Ptr{Cchar}),
                             kp.env, filename)
    if return_code != 0
        error("KNITRO: Error loading parameters from $(filename)")
    end
end

"Write all current parameter values to a file"
function save_param_file(kp::KnitroProblem, filename::AbstractString)
    return_code = @ktr_ccall(save_param_file, Int32, (Ptr{Nothing}, Ptr{Cchar}),
                             kp.env, filename)
    if return_code != 0
        error("KNITRO: Error writing parameters to $(filename)")
    end
end

"Set an integer parameter by its name"
function set_param(kp::KnitroProblem, name::AbstractString, value::Integer)
    return_code = @ktr_ccall(set_int_param_by_name, Int32, (Ptr{Nothing},
                             Ptr{Cchar}, Cint), kp.env, name, value)
    if return_code != 0
        error("KNITRO: Error setting int parameter by name")
    end
end

"Set a string parameter by its name"
function set_param(kp::KnitroProblem, name::AbstractString, value::AbstractString)
    return_code = @ktr_ccall(set_char_param_by_name, Int32, (Ptr{Nothing},
                             Ptr{Cchar}, Ptr{Cchar}), kp.env, name, value)
    if return_code != 0
        error("KNITRO: Error setting char parameter by name")
    end
end

function set_param(kp::KnitroProblem, name::AbstractString, value::Float64)
    return_code = @ktr_ccall(set_double_param_by_name, Int32, (Ptr{Nothing},
                             Ptr{Cchar}, Cdouble), kp.env, name, value)
    if return_code != 0
        error("KNITRO: Error setting float parameter by name")
    end
end

function set_param(kp::KnitroProblem, id::Integer, value::Integer)
    return_code = @ktr_ccall(set_int_param, Int32, (Ptr{Nothing}, Cint, Cint),
                             kp.env, id, value)
    if return_code != 0
        error("KNITRO: Error setting int parameter by id")
    end
end

function set_param(kp::KnitroProblem, id::Integer, value::AbstractString)
    return_code = @ktr_ccall(set_char_param, Int32, (Ptr{Nothing}, Cint,
                             Ptr{Cchar}), kp.env, id, value)
    if return_code != 0
        error("KNITRO: Error setting char parameter by id")
    end
end

function set_param(kp::KnitroProblem, id::Integer, value::Float64)
    return_code = @ktr_ccall(set_double_param, Int32, (Ptr{Nothing}, Cint,
                             Cdouble), kp.env, id, value)
    if return_code != 0
        error("KNITRO: Error setting float parameter by id")
    end
end

"Get an integer parameter by its name"
function get_param(kp::KnitroProblem, name::AbstractString, value::Vector{Int32})
    @ktr_ccall(get_int_param_by_name, Int32, (Ptr{Nothing}, Ptr{Cchar}, Ptr{Cint}),
               kp.env, name, value)
end

"Get a float parameter by its name"
function get_param(kp::KnitroProblem, name::AbstractString, value::Vector{Float64})
    @ktr_ccall(get_double_param_by_name, Int32, (Ptr{Nothing}, Ptr{Cchar},
               Ptr{Cdouble}), kp.env, name, value)
end

"Get a parameter by its name"
function get_param(kp::KnitroProblem, name::AbstractString)
    # try getting an integer parameter first
    int_result = Array(Cint, 1)
    float_result = Array(Cdouble, 1)
    if get_param(kp, name, int_result) == 0
        return int_result[1]
    elseif get_param(kp, name, float_result) != 0
        # otherwise, try getting a double parameter
        error("KNITRO: Error getting parameter by name")
    end
    float_result[1]
end

"Get an integer parameter by its id"
function get_param(kp::KnitroProblem, id::Integer, value::Vector{Int32})
    @ktr_ccall(get_int_param, Int32, (Ptr{Nothing}, Cint, Ptr{Cint}),
               kp.env, id, value)
end

"Get a float parameter by its id"
function get_param(kp::KnitroProblem, id::Integer, value::Vector{Float64})
    @ktr_ccall(get_double_param, Int32, (Ptr{Nothing}, Cint, Ptr{Cdouble}),
               kp.env, id, value)
end

"Get a parameter by its id"
function get_param(kp::KnitroProblem, id::Integer)
    # try getting an integer parameter first
    int_result = Array(Cint, 1)
    float_result = Array(Cdouble, 1)
    if get_param(kp, id, int_result) == 0
        return int_result[1]
    elseif get_param(kp, id, float_result) != 0
        # otherwise, try getting a double parameter
        error("KNITRO: Error getting parameter by id")
    end
    float_result[1]
end

### ===== NEW in KNITRO 10.0 =====
# /** Sets the string param_name with the name of parameter indexed by
#  *  param_id and returns 0. Returns an error if param_id does not
#  *  correspond to any parameter, or if the parameter output_size
#  *  (the size of char array param_name) is less than the size of the
#  *  parameter's description.
#  */
# int KNITRO_API KTR_get_param_name(      KTR_context_ptr kc,
#                                   const int             param_id,
#                                         char * const    param_name,
#                                   const size_t          output_size);

# /** Sets the string description with the description of the parameter
#  *  indexed by param_id and its possible values and returns 0. Returns an
#  *  error if param_id does not correspond to any parameter, or if the
#  *  parameter output_size (the size of char array description) is less than
#  *  the size of the parameter's description.
#  */
# int KNITRO_API KTR_get_param_doc(      KTR_context_ptr kc,
#                                  const int             param_id,
#                                        char * const    description,
#                                  const size_t          output_size);

# /** Sets the int * param_type to the parameter type of parameter indexed
#  *  by param_id. Possible values are KTR_PARAMTYPE_INTEGER, KTR_PARAMTYPE_FLOAT,
#  *  KTR_PARAMTYPE_STRING. Returns an error if param_id does not correspond to
#  *  any parameter.
#  */
# int KNITRO_API KTR_get_param_type(      KTR_context_ptr kc,
#                                   const int             param_id,
#                                         int * const     param_type);

# /** Set the int * num_param_values to the number of possible parameter
#  *  values for parameter indexed by param_id and returns 0. If there is
#  *  not a finite number of possible values, num_param_values will be zero.
#  *  Returns an error if param_id does not correspond to any parameter.
#  */
# int KNITRO_API KTR_get_num_param_values(      KTR_context_ptr kc,
#                                         const int             param_id,
#                                               int * const     num_param_values);

# /** Set string param_value_string to the description of parameter value
#  *  indexed by [param_id][value_id]. Returns an error if param_id does not
#  *  correspond to any parameter, or if value_id is greater than the number
#  *  of possible parameter values, or if there are not a finite number of
#  *  possible parameter values, or if the parameter output_size (the size
#  *  of char array param_value_string) is less than the size of the
#  *  parameter's description.
#  */
# int KNITRO_API KTR_get_param_value_doc(      KTR_context_ptr kc,
#                                        const int             param_id,
#                                        const int             value_id,
#                                              char * const    param_value_string,
#                                        const size_t          output_size);

# /** Gets the integer value corresponding to the parameter name input and
#  *  copies it into param_id input. Returns zero if successful and an error
#  *  code otherwise.
#  */
# int KNITRO_API KTR_get_param_id(      KTR_context_ptr kc,
#                                 const char * const    name,
#                                       int  * const    param_id);

### ===== END OF NEW in KNITRO 10.0 =====

"""
Similar to KTR_load_param_file but specifically allows user to
specify a file of options (and option values) to explore for
the KNITRO-Tuner.
"""
function load_tuner_file(kp::KnitroProblem, filename::AbstractString)
    return_code = @ktr_ccall(load_tuner_file, Int32, (Ptr{Nothing}, Ptr{Cchar}),
                             kp.env, filename)
    if return_code != 0
        error("KNITRO: Error loading tuner file $(filename)")
    end
end

"""
Set an array of absolute feasibility tolerances (one for each
constraint and variable) to use for the termination tests.

The user options KTR_PARAM_FEASTOL/KTR_PARAM_FEASTOLABS define
a single tolerance that is applied equally to every constraint
and variable.  This API function allows the user to specify
separate feasibility termination tolerances for each constraint
and variable.  Values specified through this function will override
the value determined by KTR_PARAM_FEASTOL/KTR_PARAM_FEASTOLABS. The
tolerances should be positive values.  If a non-positive value is
specified, that constraint or variable will use the standard tolerances
based on  KTR_PARAM_FEASTOL/KTR_PARAM_FEASTOLABS.
    `c_tol` has length m,
    `x_tol` has length n, and
    `cc_tol` has length ncc, where
    `ncc` is the number of complementarity constraints added

The regular constraints are considered to be satisfied when
    c[i] - cUpBnds[i] <= cFeasTols[i]  for all i=1..m, and
    cLoBnds[i] - c[i] <= cFeasTols[i]  for all i=1..m
The variables are considered to be satisfied when
    x[i] - xUpBnds[i] <= xFeasTols[i]  for all i=1..n, and
    xLoBnds[i] - x[i] <= xFeasTols[i]  for all i=1..n
The complementarity constraints are considered to be satisfied when
    min(x1_i, x2_i) <= ccFeasTols[i]  for all i=1..ncc,
where x1 and x2 are the arrays of complementary pairs.

If there are no regular (or complementarity) constraints set
cFeasTols=NULL (or ccFeasTols=NULL).  If cFeasTols/xFeasTols/ccFeasTols=NULL,
then the standard tolerances will be used.
"""
function set_feastols(kp::KnitroProblem,
                      c_tol::Vector{Float64},
                      x_tol::Vector{Float64},
                      cc_tol::Vector{Float64})
    return_code = @ktr_ccall(set_feastols, Int32, (Ptr{Nothing}, Ptr{Cdouble},
                             Ptr{Cdouble}, Ptr{Cdouble}), kp.env, c_tol, x_tol,
                             cc_tol)
    if return_code != 0
        error("KNITRO: Error setting feasibility tolerances")
    end
end

"""
Set names for model components passed in by the user/modeling
language so that KNITRO can internally print out these names.

KNITRO makes a local copy of all inputs, so the application may
free memory after the call.

This routine must be called after calling KTR_init_problem /
KTR_mip_init_problem and before calling KTR_solve / KTR_mip_solve.
"""
function set_names(kp::KnitroProblem, objName::AbstractString,
                   varNames::Vector{AbstractString},
                   conNames::Vector{AbstractString})
    return_code = @ktr_ccall(set_names, Int32, (Ptr{Nothing}, Ptr{Cchar},
                             Ptr{Ptr{Cchar}}, Ptr{Ptr{Cchar}}),
                             kp.env, objName, varNames, conNames)
    if return_code != 0
        error("KNITRO: Error setting names for model components")
    end
end

### NEW in KNITRO 10.0
"""
This API function can be used to identify which variables only appear
linearly in the model (KTR_LINEARVAR_YES).

This information can be used by Knitro to perform more extensive preprocessing.
If a variable appears nonlinearly in any constraint or the objective (or if the
user does not know) then it should be marked as KTR_LINEARVAR_NO. Knitro makes a
local copy of all inputs, so the application may free memory after the call.

This routine must be called after calling KTR_init_problem / KTR_mip_init_problem,
and before calling KTR_solve / KTR_mip_solve.
"""
function set_linearvars(kp::KnitroProblem, linearVars::Vector{Cint})
    return_code = @ktr_ccall(set_linearvars, Int32, (Ptr{Nothing}, Ptr{Cint}),
                             kp.env, linearVars)
    # Returns 0 if OK, nonzero if error.
    if return_code != 0
        error("KNITRO: Error setting linear variables")
    end
end

"""
This API function can be used to identify which variables should satisfy its
variable bounds throughout the optimization process (KTR_HONORBNDS_ALWAYS).

The user option KTR_PARAM_HONORBNDS can be used to set ALL variables to honor
their bounds.  This routine takes precedence over the setting of
KTR_PARAM_HONORBNDS and is used to customize the settings for individual
variables. Knitro makes a local copy of all inputs, so the application may
free memory after the call.

This routine must be called after calling KTR_init_problem / KTR_mip_init_problem,
and before calling KTR_solve / KTR_mip_solve.
"""
function set_honorbnds(kp::KnitroProblem, honorBnds::Vector{Cint})
    return_code = @ktr_ccall(set_honorbnds, Int32, (Ptr{Nothing}, Ptr{Cint}),
                             kp.env, honorBnds)
    # Returns 0 if OK, nonzero if error.
    if return_code != 0
        error("KNITRO: Error setting honor bounds")
    end
end

### END OF NEW in KNITRO 10.0

# ----- Problem modification -----

"""
This function adds complementarity constraints to the problem.

It must be called after KTR_init_problem and before KTR_solve.
The two lists are of equal length, and contain matching pairs of
variable indices.  Each pair defines a complementarity constraint
between the two variables.  The function can be called more than once
to accumulate a long list of complementarity constraints in KNITRO's
internal problem definition.
"""
function addcompcons(kp::KnitroProblem,
                     ncons::Integer,
                     index1::Vector{Int32},
                     index2::Vector{Int32})
    return_code = @ktr_ccall(addcompcons, Int32, (Ptr{Nothing}, Cint, Ptr{Cint},
                             Ptr{Cint}), kp.env, ncons, index1, index2)
    # Returns 0 if OK, or a negative value on error.
    if return_code != 0
        error("KNITRO: Error adding complementary constraints")
    end
end

"""
Prepare KNITRO to re-optimize the current problem after
modifying the variable bounds from a previous solve.

It must be called after KTR_init_problem and precedes a call to KTR_solve.
"""
function chgvarbnds(kp::KnitroProblem,
                    x_L::Vector{Float64},
                    x_U::Vector{Float64})
    return_code = @ktr_ccall(chgvarbnds, Int32, (Ptr{Nothing},Ptr{Cdouble},
                             Ptr{Cdouble}), kp.env, x_L, x_U)
    # Returns 0 if OK, nonzero if error.
    if return_code != 0
        error("KNITRO: Error modifying variable bounds")
    end
end

# /* ----- Solving ----- */
"""
Initialize Knitro with a new problem.

Brief description of the arguments (consult the Knitro manual for details):
 *    objGoal        - objective goal (KTR_OBJGOAL_MINIMIZE or _MAXIMIZE)
 *    objType        - objective type (eg, KTR_OBJTYPE_GENERAL)
 *    x_L            - (length n) lower bounds on the variables
 *    x_U            - (length n) upper bounds on the variables
 *    c_Type         - (length m) constraint type (eg, KTR_CONTYPE_GENERAL)
 *    c_L            - (length m) lower bounds on the constraints
 *    c_U            - (length m) upper bounds on the constraints
 *    jac_var        - (length nnzJ) variable index (column) of each nonzero
 *    jac_cons       - (length nnzJ) constraint index (row) of each nonzero
 *    hess_rows      - (length nnzH) row index of each nonzero
 *    hess_cols      - (length nnzH) column index of each nonzero

[Optional] Keyword arguments:
 *    initial_x      - (length n) initial start point for the unknowns
 *    initial_lambda - (length m+n) initial start point for multipliers

Knitro makes a local copy of all inputs, so the application may free memory
after the call completes.

 - If `x_L` or `x_U` are not provided, then Knitro assumes all variables
 are unbounded in that direction.
 - If "initial_x" or "initial_lambda" are not provided, then Knitro determines
 the initial start point.
 - If a particular bound does not exist, then load the array element with
 `KTR_INFBOUND` for an upper bound, or `-KTR_INFBOUND` for a lower bound.
 Knitro interprets this value as infinity and ignores the bound.
 - To make the `j`-th constraint be an equality constraint, set `c_L[j]` equal
 to `c_U[j]`.
"""
function init_problem(kp::KnitroProblem,
                      objGoal::Integer,
                      objType::Integer,
                      x_L::Vector{Float64},     # length n
                      x_U::Vector{Float64},     # length n
                      c_Type::Vector{Int32},    # length m
                      c_L::Vector{Float64},     # length m
                      c_U::Vector{Float64},     # length m
                      jac_var::Vector{Int32},   # length nnzJ
                      jac_cons::Vector{Int32},  # length nnzJ
                      hess_rows::Vector{Int32}, # length nnzH
                      hess_cols::Vector{Int32}; # length nnzH
                      initial_x = C_NULL,         # length n
                      initial_lambda = C_NULL)    # length m+n
    n = length(x_L)
    m = length(c_Type)
    nnzJ = length(jac_var)
    nnzH = length(hess_rows)
    return_code = @ktr_ccall(init_problem, Int32, (Ptr{Nothing}, Cint, Cint,
                             Cint, Ptr{Cdouble}, Ptr{Cdouble}, Cint,
                             Ptr{Cint}, Ptr{Cdouble}, Ptr{Cdouble}, Cint,
                             Ptr{Cint}, Ptr{Cint}, Cint, Ptr{Cint},
                             Ptr{Cint}, Ptr{Nothing}, Ptr{Nothing}), kp.env, n,
                             objGoal, objType, x_L, x_U, m, c_Type, c_L, c_U,
                             nnzJ, jac_var, jac_cons, nnzH, hess_rows,
                             hess_cols, initial_x, initial_lambda)
    if return_code != 0
        error("KNITRO: Error initializing problem")
    end
end

function init_problem(kp::KnitroProblem,
                      objGoal::Integer,
                      objType::Integer,
                      x_L::Vector{Float64},     # length n
                      x_U::Vector{Float64},     # length n
                      c_Type::Vector{Int32},    # length m
                      c_L::Vector{Float64},     # length m
                      c_U::Vector{Float64},     # length m
                      jac_var::Vector{Int32},   # length nnzJ
                      jac_cons::Vector{Int32};  # length nnzJ
                      initial_x = C_NULL,       # length n
                      initial_lambda = C_NULL)  # length m+n
    n = length(x_L)
    m = length(c_Type)
    nnzJ = length(jac_var)
    return_code = @ktr_ccall(init_problem, Int32, (Ptr{Nothing}, Cint, Cint, Cint,
                             Ptr{Cdouble}, Ptr{Cdouble}, Cint, Ptr{Cint},
                             Ptr{Cdouble}, Ptr{Cdouble}, Cint, Ptr{Cint},
                             Ptr{Cint}, Cint, Ptr{Nothing}, Ptr{Nothing}, Ptr{Nothing},
                             Ptr{Nothing}), kp.env, n, objGoal, objType, x_L, x_U,
                             m, c_Type, c_L, c_U, nnzJ, jac_var, jac_cons, 0,
                             C_NULL, C_NULL, initial_x, initial_lambda)
    if return_code != 0
        error("KNITRO: Error initializing problem without exact hessians")
    end
end

"""
Call KNITRO to solve the problem.

If the application provides callback functions for evaluating the function,
constraints, and derivatives, then a single call to `solve_problem` returns
the solution.  Otherwise, KNITRO operates in reverse communications mode and
returns a status code that may request another call.

Returns one of the status codes KTR_RC_*. In particular:
<0 - Knitro terminated with a fatal error
0  - KNITRO is finished: x, lambda, and obj contain the optimal solution
1  - call `solve_problem` again (reverse comm) with obj and c containing
     the objective and constraints evaluated at x
2  - call `solve_problem` again (reverse comm) with objGrad and jac containing
     the objective and constraint first derivatives evaluated at x
3  - call `solve_problem` again (reverse comm) with hess containing
     H(x,lambda), the Hessian of the Lagrangian evaluated at x and lambda
7  - call `solve_problem` again (reverse comm) with hessVector containing
     the result of H(x,lambda) * hessVector

Brief description of the arguments (consult the KNITRO manual for details):
x          - output (length n) solution point estimate
lambda     - output (length m+n) Lagrange multiplier estimate
evalStatus - input  evaluation status (0=OK)
obj        - input  (length 1) objective at x
             output            optimal objective when finished
cons       - input  (length m) constraints at x
objGrad    - input  (length n) objective gradient at x
jac        - input  (length nnzJ) sparse constraint gradient at x
hess       - input  (length nnzH) sparse Hessian at x and lambda
hessVector - input  (length n) result of H(x,lambda) * hessVector
             output            vector to multiply Hessian by

If `gradopt` is set to compute finite differences for first derivatives,
then `solve_problem` will modify `objGrad` and `jac`; otherwise, these
arguments are not modified.
"""
function solve_problem(kp::KnitroProblem,
                       x::Vector{Float64},
                       lambda::Vector{Float64},
                       evalStatus::Integer,
                       obj::Vector{Float64},
                       cons::Vector{Float64},
                       objGrad::Vector{Float64},
                       jac::Vector{Float64},
                       hess::Vector{Float64},
                       hessVector::Vector{Float64})
    @ktr_ccall(solve, Int32, (Ptr{Nothing}, Ptr{Cdouble},
               Ptr{Cdouble}, Cint, Ptr{Cdouble}, Ptr{Cdouble},
               Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble},
               Ptr{Cdouble}, Any), kp.env, x, lambda, evalStatus,
               obj, cons, objGrad, jac, hess, hessVector, kp)
end

function solve_problem(kp::KnitroProblem,
                       x::Vector{Float64},
                       lambda::Vector{Float64},
                       evalStatus::Integer,
                       obj::Vector{Float64})
    # For callback mode
    @ktr_ccall(solve, Int32, (Ptr{Nothing}, Ptr{Cdouble},
               Ptr{Cdouble}, Cint, Ptr{Cdouble}, Ptr{Nothing},
               Ptr{Nothing}, Ptr{Nothing}, Ptr{Nothing}, Ptr{Nothing}, Any),
               kp.env, x, lambda, evalStatus, obj, C_NULL,
               C_NULL, C_NULL, C_NULL, C_NULL, kp)
end

"""
Prepare KNITRO to restart the current problem at start point [x_0, lambda_0].

If output to a file is enabled, this will erase the current file.
KNITRO parameter values are not changed by this call.
"""
function restart_problem(kp::KnitroProblem,
                         x_0::Vector{Cdouble},
                         lambda_0::Vector{Cdouble})
    return_code = @ktr_ccall(restart, Int32, (Ptr{Nothing}, Ptr{Cdouble},
                             Ptr{Cdouble}), kp.env, x_0, lambda_0)
    if return_code != 0
        error("KNITRO: Error restarting problem")
    end
end

"""
Prepare KNITRO to restart the current problem.

If output to a file is enabled, this will erase the current file.
KNITRO parameter values are not changed by this call.
"""
function restart_problem(kp::KnitroProblem)
    return_code = @ktr_ccall(restart, Int32, (Ptr{Nothing}, Ptr{Cdouble},
                             Ptr{Cdouble}), kp.env, C_NULL, C_NULL)
    if return_code != 0
        error("KNITRO: Error restarting problem")
    end
end

"""
Initialize KNITRO with a new MIP problem.  KNITRO makes a local copy of
all inputs, so the application may free memory after the call completes.
"""
function mip_init_problem(kp::KnitroProblem,
                          objGoal::Integer,
                          objType::Integer,
                          objFnType::Integer,
                          x_Type::Vector{Int32},    # length n
                          x_L::Vector{Float64},     # length n
                          x_U::Vector{Float64},     # length n
                          c_Type::Vector{Int32},    # length m
                          c_FnType::Vector{Int32},  # length m
                          c_L::Vector{Float64},     # length m
                          c_U::Vector{Float64},     # length m
                          jac_var::Vector{Int32},   # length nnzJ
                          jac_cons::Vector{Int32},  # length nnzJ
                          hess_rows::Vector{Int32}, # length nnzH
                          hess_cols::Vector{Int32}; # length nnzH
                          initial_x = C_NULL,       # length n
                          initial_lambda = C_NULL)  # length m+n
    n = length(x_L)
    m = length(c_Type)
    nnzJ = length(jac_var)
    nnzH = length(hess_rows)
    return_code = @ktr_ccall(mip_init_problem, Int32, (Ptr{Nothing}, Cint, Cint,
                             Cint, Cint, Ptr{Cint}, Ptr{Cdouble}, Ptr{Cdouble},
                             Cint, Ptr{Cint}, Ptr{Cint}, Ptr{Cdouble},
                             Ptr{Cdouble}, Cint, Ptr{Cint}, Ptr{Cint},Cint,
                             Ptr{Cint},Ptr{Cint}, Ptr{Nothing}, Ptr{Nothing}),
                             kp.env, n, objGoal, objType, objFnType, x_Type,
                             x_L, x_U, m, c_Type, c_FnType, c_L, c_U, nnzJ,
                             jac_var, jac_cons, nnzH, hess_rows, hess_cols,
                             initial_x, initial_lambda)
    if return_code != 0
        error("KNITRO: Error initializing MIP problem")
    end
end

function mip_init_problem(kp::KnitroProblem,
                          objGoal::Integer,
                          objType::Integer,
                          objFnType::Integer,
                          x_Type::Vector{Int32},    # length n
                          x_L::Vector{Float64},     # length n
                          x_U::Vector{Float64},     # length n
                          c_Type::Vector{Int32},    # length m
                          c_FnType::Vector{Int32},  # length m
                          c_L::Vector{Float64},     # length m
                          c_U::Vector{Float64},     # length m
                          jac_var::Vector{Int32},   # length nnzJ
                          jac_cons::Vector{Int32};  # length nnzJ
                          initial_x = C_NULL,       # length n
                          initial_lambda = C_NULL)  # length m+n
    n = length(x_L)
    m = length(c_Type)
    nnzJ = length(jac_var)
    return_code = @ktr_ccall(mip_init_problem, Int32, (Ptr{Nothing}, Cint, Cint,
                             Cint, Cint, Ptr{Cint}, Ptr{Cdouble}, Ptr{Cdouble},
                             Cint, Ptr{Cint}, Ptr{Cint}, Ptr{Cdouble},
                             Ptr{Cdouble},Cint,Ptr{Cint}, Ptr{Cint}, Cint,
                             Ptr{Nothing}, Ptr{Nothing}, Ptr{Nothing}, Ptr{Nothing}),
                             kp.env, n, objGoal, objType, objFnType, x_Type,
                             x_L, x_U, m, c_Type, c_FnType, c_L, c_U, nnzJ,
                             jac_var, jac_cons, 0, C_NULL, C_NULL,
                             initial_x, initial_lambda)
    if return_code != 0
        error("KNITRO: Error initializing MIP problem without exact hessians")
    end
end

"""
Set the branching priorities for integer variables.

Priorities must be positive numbers (variables with non-positive values
are ignored).  Variables with higher priority values will be considered
for branching before variables with lower priority values.  When
priorities for a subset of variables are equal, the branching rule is
applied as a tiebreaker.
"""
function set_branching_priorities(kp::KnitroProblem,
                                  xPriorities::Vector{Int})
    return_code = @ktr_ccall(mip_set_branching_priorities, Int32, (Ptr{Nothing},
                             Ptr{Cint}), kp.env, xPriorities)
    if return_code != 0
        error("KNITRO: Error setting MIP branching priorities")
    end
end

"""
Call KNITRO to solve the MIP problem.

If the application provides callback functions for evaluating the function,
constraints, and derivatives, then a single call to KTR_mip_solve returns
the solution.
Otherwise, KNITRO operates in reverse communications mode and returns a
status code that may request another call.

Returns one of the status codes KTR_RC_*. In particular:
0 - KNITRO is finished: x, lambda, and obj contain the optimal solution
1 - call KTR_solve again (reverse comm) with obj and c containing
    the objective and constraints evaluated at x
2 - call KTR_solve again (reverse comm) with objGrad and jac containing
    the objective and constraint first derivatives evaluated at x
3 - call KTR_solve again (reverse comm) with hess containing
    H(x,lambda), the Hessian of the Lagrangian evaluated at x and lambda
7 - call KTR_solve again (reverse comm) with hessVector containing
    the result of H(x,lambda) * hessVector

If `gradopt` is set to compute finite differences for first derivatives,
then KTR_mip_solve will modify objGrad and jac; otherwise, these arguments
are not modified.
"""
function mip_solve_problem(kp::KnitroProblem,
                           x::Vector{Float64},
                           lambda::Vector{Float64},
                           evalStatus::Integer,
                           obj::Vector{Float64},
                           cons::Vector{Float64},
                           objGrad::Vector{Float64},
                           jac::Vector{Float64},
                           hess::Vector{Float64},
                           hessVector::Vector{Float64})
    @ktr_ccall(mip_solve, Int32, (Ptr{Nothing}, Ptr{Cdouble},
               Ptr{Cdouble}, Cint, Ptr{Cdouble}, Ptr{Cdouble},
               Ptr{Cdouble},Ptr{Cdouble}, Ptr{Cdouble},
               Ptr{Cdouble}, Any), kp.env, x, lambda, evalStatus,
               obj, cons, objGrad, jac, hess, hessVector, kp)
end

function mip_solve_problem(kp::KnitroProblem,
                           x::Vector{Float64},
                           lambda::Vector{Float64},
                           evalStatus::Integer,
                           obj::Vector{Float64})
    # For callback mode
    @ktr_ccall(mip_solve, Int32, (Ptr{Nothing}, Ptr{Cdouble},
               Ptr{Cdouble}, Cint, Ptr{Cdouble}, Ptr{Nothing},
               Ptr{Nothing}, Ptr{Nothing}, Ptr{Nothing}, Ptr{Nothing}, Any),
               kp.env, x, lambda, evalStatus, obj, C_NULL,
               C_NULL, C_NULL, C_NULL, C_NULL, kp)
end

"""
Set an array of relative stepsizes to use for the finite-difference
gradient/Jacobian computations when using finite-difference
first derivatives.

Finite-difference step sizes `delta` in KNITRO are
computed as:

          delta[i] = relStepSizes[i]*max(abs(x[i]),1)

The default relative step sizes for each component of `x` are `sqrt(eps)`
for forward finite differences, and `eps^(1/3)` for central finite
differences.  Use this function to overwrite the default values.
Array relStepSizes has length n and all values should be non-zero.
"""
function set_findiff_relstepsizes(kp::KnitroProblem,
                                  relStepSizes::Vector{Float64})
    return_code = @ktr_ccall(set_findiff_relstepsizes, Int32, (Ptr{Nothing},
                             Ptr{Cdouble}), kp.env, relStepSizes)
    if return_code != 0
        error("KNITRO: Error setting relative stepsizes for the
              finite-difference gradient/Jacobian computations")
    end
end

# ----- Reading solution properties -----

"""
Return the number of function evaluations requested by KTR_solve.

A single request evaluates the objective and all constraint functions.
"""
function get_number_FC_evals(kp::KnitroProblem)
    n = @ktr_ccall(get_number_FC_evals, Int32, (Ptr{Nothing},), kp.env)
    if n < 0
        error("KNITRO: Error returning number of function evaluations")
    end
    n
end

"""
Return the number of gradient evaluations requested by KTR_solve.

A single request evaluates first derivatives of the objective and
all constraint functions.
"""
function get_number_GA_evals(kp::KnitroProblem)
    n = @ktr_ccall(get_number_GA_evals, Int32, (Ptr{Nothing},), kp.env)
    if n < 0
        error("KNITRO: Error returning number of gradient evaluations")
    end
    n
end

"""
Return the number of Hessian evaluations requested by KTR_solve.

A single request evaluates second derivatives of the objective and
all constraint functions.
"""
function get_number_H_evals(kp::KnitroProblem)
    n = @ktr_ccall(get_number_H_evals, Int32, (Ptr{Nothing},), kp.env)
    if n < 0
        error("KNITRO: Error returning number of hessian evaluations")
    end
    n
end

"""
Return the number of Hessian-vector products requested by KTR_solve.

A single request evaluates the product of the Hessian of the
Lagrangian with a vector submitted by KNITRO.
"""
function get_number_HV_evals(kp::KnitroProblem)
    n = @ktr_ccall(get_number_HV_evals, Int32, (Ptr{Nothing},), kp.env)
    if n < 0
        error("KNITRO: Error returning number of hessian-vector evaluations")
    end
    n
end

# /* ----- Solution properties for continuous problems only ----- */

"Return the number of iterations made by KTR_solve."
function get_number_iters(kp::KnitroProblem)
    n = @ktr_ccall(get_number_iters, Int32, (Ptr{Nothing},), kp.env)
    if n < 0
        error("KNITRO: Error returning number of iterations")
    end
    n
end

"Return the number of conjugate gradient (CG) iterations by KTR_solve"
function get_number_cg_iters(kp::KnitroProblem)
    n = @ktr_ccall(get_number_cg_iters, Int32, (Ptr{Nothing},), kp.env)
    if n < 0
        error("KNITRO: Error returning no. of conjugate gradient iterations")
    end
    n
end

"""
Return the absolute feasibility error at the solution

Refer to the KNITRO manual section on Termination Tests for a
detailed definition of this quantity.
"""
function get_abs_feas_error(kp::KnitroProblem)
    n = @ktr_ccall(get_abs_feas_error, Float64, (Ptr{Nothing},), kp.env)
    if n < 0
        error("KNITRO: Error getting absolute feasibility error at solution")
    end
    n
end

"""
Return the relative feasibility error at the solution

Refer to the KNITRO manual section on Termination Tests for a
detailed definition of this quantity.
"""
function get_rel_feas_error(kp::KnitroProblem)
    n = @ktr_ccall(get_rel_feas_error, Float64, (Ptr{Nothing},), kp.env)
    if n < 0
        error("KNITRO: Error getting relative feasibility error at solution")
    end
    n
end

"""
Return the absolute optimality error at the solution.

Refer to the KNITRO manual section on Termination Tests for a
detailed definition of this quantity.
"""
function get_abs_opt_error(kp::KnitroProblem)
    n = @ktr_ccall(get_abs_opt_error, Float64, (Ptr{Nothing},), kp.env)
    if n < 0
        error("KNITRO: Error getting absolute optimality error at solution")
    end
    n
end

"""
Return the relative optimality error at the solution.

Refer to the KNITRO manual section on Termination Tests for a
detailed definition of this quantity.
"""
function get_rel_opt_error(kp::KnitroProblem)
    n = @ktr_ccall(get_rel_opt_error, Float64, (Ptr{Nothing},), kp.env)
    if n < 0
        error("KNITRO: Error getting relative optimality error at solution")
    end
    n
end

"""
Return the solution status, objective, primal and dual variables.

The `status`, `obj`, `x`, and `lambda` vectors will be modified in-place
with the values returned by the routine.
"""
function get_solution(kp::KnitroProblem,
                      status::Vector{Int32},
                      obj::Vector{Float64},
                      x::Vector{Float64},
                      lambda::Vector{Float64})
    return_code = @ktr_ccall(get_solution, Int32, (Ptr{Nothing}, Ptr{Cint},
                             Ptr{Cdouble},Ptr{Cdouble},Ptr{Cdouble}), kp.env,
                             status, obj, x, lambda)
    if return_code < 0
        error("KNITRO: Error getting the solution status and values")
    end
end

"Return the values of the constraint vector c(x) in `c` through `cons`"
function get_constraint_values(kp::KnitroProblem, cons::Vector{Float64})
    return_code = @ktr_ccall(get_constraint_values, Int32, (Ptr{Nothing},
                             Ptr{Cdouble}), kp.env, cons)
    if return_code < 0
        error("KNITRO: Error getting the values of the constraint vector c(x)")
    end
end

"Return the values of the objective gradient vector through `objGrad`"
function get_objgrad_values(kp::KnitroProblem, objGrad::Vector{Float64})
    return_code = @ktr_ccall(get_objgrad_values, Int32, (Ptr{Nothing},
                             Ptr{Cdouble}), kp.env, objGrad)
    if return_code < 0
        error("KNITRO: Error getting the objective gradient vector")
    end
end

"Return the values of the (non-zero) sparse Jacobian through `jac`"
function get_jacobian_values(kp::KnitroProblem, jac::Vector{Float64})
    return_code = @ktr_ccall(get_jacobian_values, Int32, (Ptr{Nothing},
                             Ptr{Cdouble}), kp.env, jac)
    if return_code < 0
        error("KNITRO: Error getting the values of the sparse Jacobian")
    end
end

"""
Return the values of the (non-zero sparse) Hessian (or possibly Hessian
approximation) through `hess`.  This routine is currently only valid
if 1 of the 2 following cases holds:
  1) KTR_HESSOPT_EXACT (presolver on or off), or;
  2) KTR_HESSOPT_BFGS or KTR_HESSOPT_SR1, but only with the
     KNITRO presolver off (i.e. KTR_PRESOLVE_NONE).

In all other cases, either KNITRO does not have an internal
representation of the Hessian (or Hessian approximation),
or the internal Hessian approximation corresponds only to
the presolved problem form and may not be valid for the
original problem form.  In these cases `hess` is left
unmodified, and the routine has return code 1.

Note that in case 2 above (KTR_HESSOPT_BFGS or KTR_HESSOPT_SR1)
the values returned in `hess` are the upper triangular values
of the dense quasi-Newton Hessian approximation stored row-wise.
There are ((n*n - n)/2 + n) such values (where `n` is the number
of variables in the problem. These values may be quite different
from the values of the exact Hessian.

When KTR_HESSOPT_EXACT (case 1 above) the Hessian values
returned correspond to the non-zero sparse Hessian indices
provided by the user in `init_problem()`.

Returns 0 if call is successful;
        1 if `hess` was not set because KNITRO does not
        have a valid Hessian for the model stored.
"""
function get_hessian_values(kp::KnitroProblem,
                            hess::Vector{Float64})
    return_code = @ktr_ccall(get_hessian_values, Int32, (Ptr{Nothing},
                             Ptr{Cdouble}), kp.env, hess)
    if return_code < 0
        error("KNITRO: Error getting the values of the Hessian")
    end
    return_code
end

# /* ----- Solution properties for MIP problems only ----- */

"Return the number of nodes processed in MIP solve"
function get_mip_num_nodes(kp::KnitroProblem)
    n = @ktr_ccall(get_mip_num_nodes, Int32, (Ptr{Nothing},), kp.env)
    if n < 0
        error("KNITRO: Error getting the number of nodes processed in
              MIP solve")
    end
    n
end

"Return the number of continuous subproblems processed in MIP solve."
function get_mip_num_solves(kp::KnitroProblem)
    n = @ktr_ccall(get_mip_num_solves, Int32, (Ptr{Nothing},), kp.env)
    if n < 0
        error("KNITRO: Error getting the number of subproblems processed in
              MIP solve")
    end
    n
end

"""
Return the final absolute integrality gap in the MIP solve.

Returns KTR_INFBOUND if no incumbent (i.e., integer feasible) point found.

Refer to the KNITRO manual section on Termination Tests for
a detailed definition of this quantity.
"""
function get_mip_abs_gap(kp::KnitroProblem)
    return_code = @ktr_ccall(get_mip_abs_gap, Float64, (Ptr{Nothing},), kp.env)
    if return_code == KTR_RC_BAD_KCPTR
        error("KNITRO: Error getting the final absolute integrality gap in
              MIP solve.")
    end
    return_code
end

"""
Return the final relative integrality gap in the MIP solve.

Refer to the KNITRO manual section on Termination Tests for
a detailed definition of this quantity.
"""
function get_mip_rel_gap(kp::KnitroProblem)
    return_code = @ktr_ccall(get_mip_rel_gap, Float64, (Ptr{Nothing},), kp.env)
    if return_code == KTR_RC_BAD_KCPTR
        error("KNITRO: Error getting the final relative integrality gap in
              MIP solve.")
    end
    return_code
end

"Return the objective value of the MIP incumbent solution"
function get_mip_incumbent_obj(kp::KnitroProblem)
    return_code = @ktr_ccall(get_mip_incumbent_obj, Float64, (Ptr{Nothing},),
                             kp.env)
    if return_code == KTR_RC_BAD_KCPTR
        error("KNITRO: Error getting the objective value of the MIP incumbent
              solution.")
    end
    return_code
end

"Return the value of the current MIP relaxation bound"
function get_mip_relaxation_bnd(kp::KnitroProblem)
    return_code = @ktr_ccall(get_mip_relaxation_bnd, Float64, (Ptr{Nothing},),
                             kp.env)
    if return_code == KTR_RC_BAD_KCPTR
        error("KNITRO: Error getting the value of the current MIP relaxation
              bound.")
    end
    return_code
end

"Return the objective value of the most recent MIP node subproblem"
function get_mip_lastnode_obj(kp::KnitroProblem)
    return_code = @ktr_ccall(get_mip_lastnode_obj, Float64, (Ptr{Nothing},),
                             kp.env)
    if return_code == KTR_RC_BAD_KCPTR
        error("KNITRO: Error getting the objective value of the most recent
              MIP subproblem.")
    end
    return_code
end

"""
Return the MIP incumbent solution in 'x' if it exists

Returns 1 if incumbent solution exists and call is successful;
        0 if no incumbent (i.e., integer feasible) exists
"""
function get_mip_incumbent_x(kp::KnitroProblem, x::Vector{Float64})
    return_code = @ktr_ccall(get_mip_incumbent_x, Int32, (Ptr{Nothing},
                             Ptr{Cdouble}), kp.env, x)
    if return_code < 0
        error("KNITRO: Error getting the MIP incumbent solution 'x'.")
    end
    return_code
end

"""
Compare the application's analytic first derivatives to a finite
difference approximation at x.  The objective and all constraint
functions are checked.

Returns one of the status codes KTR_RC_*. In particular:
0 - routine is finished
1 - call routine again (reverse comm) with obj and c containing
    the objective and constraints evaluated at x
2 - call routine again (reverse comm) with objGrad and jac containing
    the objective and constraint first derivatives evaluated at x

Description of the arguments:
x                - input  (length n) point at which to check derivatives
                   output            point at which to evaluate obj and c
finiteDiffMethod - 1 = forward differences, 2 = central differences
absThreshold     - print when |estimate - analytic| > threshold
relThreshold     - print when |estimate - analytic| > threshold * scale
                   where scale = max{1, |analytic|}
evalStatus       - input            evaluation status (0=OK)
obj              - input            objective at x
c                - input (length m) constraints at x
objGrad          - input (length n) analytic gradient at x
jac              - input (length nnzJ) analytic constraint Jacobian at x
"""
function check_first_ders(kp::KnitroProblem,
                          x::Vector{Float64},
                          finiteDiffMethod::Integer,
                          absThreshold::Float64,
                          relThreshold::Float64,
                          evalStatus::Integer,
                          obj::Float64,
                          cons::Vector{Float64},
                          objGrad::Vector{Float64},
                          jac::Vector{Float64})
    return_code = @ktr_ccall(check_first_ders, Int32, (Ptr{Nothing}, Ptr{Cdouble},
                             Cint, Cdouble, Cdouble, Cint, Cdouble,
                             Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble},
                             Ptr{Nothing}), kp.env, x, finiteDiffMethod,
                             absThreshold, relThreshold, evalStatus,
                             obj, cons, objGrad, jac, C_NULL)
    if return_code < 0
        error("KNITRO: Error comparing the application's analytic first
              derivatives to a finite difference approximation at x")
    end
end


text:
__precompile__(true)
module XLSX

<<<<<<< HEAD
import ZipFile
import EzXML
using Dates, Printf
=======
import Dates
using Printf

import ZipFile, EzXML
>>>>>>> upstream/master

# https://github.com/fhs/ZipFile.jl/issues/39
if !hasmethod(Base.bytesavailable, Tuple{ZipFile.ReadableFile})
    Base.bytesavailable(f::ZipFile.ReadableFile) = f.uncompressedsize - f._pos
end

const SPREADSHEET_NAMESPACE_XPATH_ARG = [ "xpath" => "http://schemas.openxmlformats.org/spreadsheetml/2006/main" ]

include("structs.jl")
include("cellref.jl")
include("sst.jl")
include("stream.jl")
include("table.jl")
include("relationship.jl")
include("read.jl")
include("workbook.jl")
include("worksheet.jl")
include("cell.jl")
include("styles.jl")
include("write.jl")

end # module XLSX

ext: jl


text:
export summarize_psd, get_freqs_and_periods

using Statistics: mean, median

function decibel(value::Real)
    return 10 * log10(value)
end

function calculate_freq_range(sampling_rate::Float64, N::Int64)
    delta = 1. / sampling_rate
    total_duration = delta * N
    range_arr = Array{Float64}(undef, N)
    for i in 1:N
        range_arr[i] = (i - 1) / total_duration
    end

    return range_arr
end

function get_freqs_and_periods(sampling_rate::Float64, window_length::Int64, smoothing_width_factor::Float64)
    nyquiest_freq = sampling_rate / 2.0
    octave_step = 2 ^ 0.125
    count = 0

    # Short period (i.e. high freqnency)
    fh = nyquiest_freq
    ts = 1. / fh
    # Long period (i.e. low freqnency)
    fl = 0.0
    tl = smoothing_width_factor * ts
    # Longest resolvable period
    tr = window_length / 5.0

    # Count the lenght of left/right frequency array length
    while ts <= tr && tl <= tr
        count += 1
        ts *= octave_step
        tl *= octave_step
    end

    left_freqs = Array{Float64}(undef, count)
    right_freqs = Array{Float64}(undef, count)
    center_periods = Array{Float64}(undef, count)

    # Set the left and right frequency of each octave
    ts = 1. / fh
    tl = smoothing_width_factor * ts
    for i in 1:count
        fh = 1. / ts
        fl = 1. / tl
        left_freqs[i] = fh;
        right_freqs[i] = fl;
        center_periods[i] = sqrt(ts * tl)
        ts *= octave_step
        tl *= octave_step
    end

    return left_freqs, right_freqs, center_periods
end

function summarize_psd(psd_bin::AbstractArray{<:Real, 2}, sampling_rate::Float64, smoothing_width_factor::Float64)
    num_segments, N = size(psd_bin)
    """psd_result = Array{eltype(psd_bin)}(undef, 4, N)
    for i in 1:N
        psd_result[1, i] = mean(psd_bin[:, i]) # mean

        min_and_max = extrema(psd_bin[:, i])
        psd_result[2, i] = min_and_max[1] # min
        psd_result[3, i] = min_and_max[2] # max

        psd_result[4, i] = median(psd_bin[:, i]) # median
    end"""

    # Dimension reduction technique escribed in McMarana 2004
    left_freqs, right_freqs, center_periods = get_freqs_and_periods(sampling_rate, (Int64)(N / sampling_rate), smoothing_width_factor)
    len_freqs = size(left_freqs, 1)
    estimated_freqs = calculate_freq_range(sampling_rate, N)
    #psd_bin_reduced = zeros(eltype(psd_bin), num_segments, len_freqs)
    psd_bin_reduced = Array{eltype(psd_bin)}(undef, num_segments, len_freqs)
    for i in 1:num_segments
        for j in 1:len_freqs
            count = 0
            psd_bin_reduced[i, j] = 0
            for k in 1:N
                if left_freqs[j] >= estimated_freqs[k] && estimated_freqs[k] >= right_freqs[j]
                    psd_bin_reduced[i, j] += psd_bin[i, k]
                    count += 1
                end
            end
            psd_bin_reduced[i, j] /= count
        end
    end

    # Set unit to decibel (dB)
    psd_result = Array{eltype(psd_bin_reduced)}(undef, 4, len_freqs)
    for i in 1:len_freqs
        psd_result[1, i] = decibel(mean(psd_bin_reduced[:, i])) # mean

        min_and_max = extrema(psd_bin_reduced[:, i])
        psd_result[2, i] = decibel(min_and_max[1]) # min
        psd_result[3, i] = decibel(min_and_max[2]) # max

        psd_result[4, i] = decibel(median(psd_bin_reduced[:, i])) # median
    end

    return psd_result, center_periods
end


text:
module Vec

using StaticArrays
using LinearAlgebra
using Printf
import LinearAlgebra: ⋅, ×

export
    AbstractVec,
    VecE,  # an abstract euclidean-group vector
    VecSE, # an abstract special euclidean-group element

    VecE2, # two-element Float64 vector, {x, y}
    VecE3, # three-element Float64 vector, {x, y, z}
    VecSE2, # point in special euclidean group of order 2, {x, y, θ}

    polar, # construct a polar vector with (r,θ)

    proj,  # vector projection
           # proj(a::vec, b::vec, ::Type{Float64}) will do scalar projection of a onto b
           # proj(a::vec, b::vec, ::Type{Vec}) will do vector projection of a onto b

    lerp,  # linear interpolation between two vec's
    invlerp,
    lerp_angle,

    normsquared,         # see docstrings below
    scale_euclidean,
    clamp_euclidean,
    normalize_euclidian,

    dist,  # scalar distance between two vec's
    dist2, # squared scalar distance between two vec's

    rot,   # rotate the vector, always using the Right Hand Rule
    rot_normalized, # like rot, but assumes axis is normalized

    deltaangle, # signed delta angle
    angledist,  # distance between two angles

    inertial2body,
    body2inertial,

    orientation,
    are_collinear,
    get_intersection,

    Circ,
    AABB,
    OBB

abstract type AbstractVec{N, R} <: FieldVector{N, R} end
abstract type VecE{N, R} <: AbstractVec{N, R} end
abstract type VecSE{N, R} <: AbstractVec{N, R} end

lerp(a::Real, b::Real, t::Real) = a + (b-a)*t
invlerp(a::Real, b::Real, c::Real) = (c - a)/(b-a)

"The L2 norm squared."
normsquared(a::AbstractVec) = sum(x^2 for x in a)

"Scale the euclidean part of the vector by a factor b while leaving the orientation part unchanged."
scale_euclidean(a::VecE, b::Real) = b.*a

"Clamp each element of the euclidean part of the vector while leaving the orientation part unchanged."
clamp_euclidean(a::VecE, lo::Real, hi::Real) = clamp.(a, lo, hi)

"Normalize the euclidean part of the vector while leaving the orientation part unchanged."
normalize_euclidian(a::VecE, p::Real=2) = normalize(a, p)

include("common.jl")
include("vecE2.jl")
include("vecE3.jl")
include("vecSE2.jl")

include("geom/geom.jl")
include("coordinate_transforms.jl")
include("quat.jl")

Base.vec(v::Union{AbstractVec, Quat, RPY}) = convert(Vector{Float64}, v)

end # module


text:
using BenchmarkTools,Knet

const N=1000

const sizes = (1,10,100,1000)

function f2(x,r)
    nx = length(x)
    if r==0
        for i=1:N; y=ccall(("sum_32_20",Knet.libknet8),Float32,(Cint,Ptr{Float32}),nx,x); end
    elseif r==1
        y = similar(x, 1, size(x,2)); ny = length(y); sy = stride(x,2)
        for i=1:N; ccall(("sum_32_21",Knet.libknet8),Nothing,(Cint,Ptr{Float32},Cint,Cint,Ptr{Float32}),nx,x,sy,ny,y); end
    elseif r==2
        y = similar(x, size(x,1), 1); ny = length(y); sy = stride(x,1)
        for i=1:N; ccall(("sum_32_21",Knet.libknet8),Nothing,(Cint,Ptr{Float32},Cint,Cint,Ptr{Float32}),nx,x,sy,ny,y); end
    else
        error("not supported")
    end
    Knet.cudaDeviceSynchronize()
    return y
end

for region in (0,1,2)
    println("sum(a,$region)")
    for s in sizes; print("\t$s"); end; println()
    for nrows in sizes
        print(nrows)
        for ncols in sizes
            a = KnetArray(rand(Float32,nrows,ncols))
            bm = @benchmark f2($a,$region) seconds=1
            mm = round(Int, minimum(bm.times)/N)
            print("\t$mm")
            Knet.gc()
        end
        println()
    end
end

# commit 359d3646 2018-09-22, julia 1.0.0 vs commit 4aa5f92f 2018-08-14, julia 0.6.4
# GPU:V100, CPU:Intel(R) Xeon(R) CPU E5-2680 v4 @ 2.40GHz
#
# sum(a,0)                      		sum(a,0)
# 	1	10	100	1000            	1	10	100	1000
# 1	17834	17643	18922	19007           1	17134	17165	18364	18395
# 10	18977	18978	18977	19006           10	18448	18362	18411	18404
# 100	18974	18965	18989	19064           100	18406	17920	17924	17543
# 1000	18954	18945	19013	27356           1000	17560	17205	17169	25933
# sum(a,1)                                      sum(a,1)
# 	1	10	100	1000            	1	10	100	1000
# 1	3335	3326	3322	5727            1	3188	3190	3191	5721
# 10	3327	3339	3326	6611            10	3189	3192	3191	6596
# 100	3319	3595	3588	6796            100	3492	3482	3486	6777
# 1000	3615	3637	3630	11536           1000	3501	3503	3380	11583
# sum(a,2)                                      sum(a,2)
# 	1	10	100	1000            	1	10	100	1000
# 1	3641	3615	3623	3620            1	3347	3203	3208	3300
# 10	3533	3441	3458	3842            10	3200	3197	3206	3804
# 100	3535	3447	3341	4300            100	3029	3033	3026	4295
# 1000	5782	5974	6923	24663           1000	5767	5967	7001	24780

# COMMIT 3d32e16 2017-05-17 merge
#
# sum(a,0)
# 	1	10	100	1000
# 1	31804	31817	31814	31822
# 10	31790	31852	31827	32158
# 100	31848	31876	32085	35014
# 1000	31831	32072	35023	65961
# sum(a,1)
# 	1	10	100	1000
# 1	4689	4689	4931	13799
# 10	4936	4938	5436	16528
# 100	4936	4938	5470	17198
# 1000	6457	6608	7302	48299
# sum(a,2)
# 	1	10	100	1000
# 1	4690	4689	4691	6779
# 10	4690	4748	4885	8096
# 100	4935	5208	6864	26480
# 1000	13809	15307	32166	216076

# COMMIT: ccad0cb 2017-04-15 (added general reductions, same as last)
# sum(a,0)
# 	1	10	100	1000
# 1	31741	31853	31864	31853
# 10	31779	31804	31899	32097
# 100	31818	31895	32155	35172
# 1000	31871	32094	35185	66605
# sum(a,1)
# 	1	10	100	1000
# 1	4619	4689	4892	13813
# 10	4885	4905	5404	16503
# 100	4848	4914	5442	17207
# 1000	6390	6579	7300	48356
# sum(a,2)
# 	1	10	100	1000
# 1	4611	4621	4606	6692
# 10	4689	4731	4830	8054
# 100	4887	5104	6839	26471
# 1000	13841	15296	32159	217215

# COMMIT: c78d436 2017-04-04 (fixed large array bug, but slowed last col sum(a,1) and last row sum(a,2))
# sum(a,0)
# 	1	10	100	1000
# 1	32333	32391	32305	32404
# 10	32388	32403	32440	32703
# 100	32411	32407	32668	35571
# 1000	32362	32618	35586	66601
# sum(a,1)
# 	1	10	100	1000
# 1	4690	4690	4924	13812
# 10	4939	4941	5427	16532
# 100	4932	4938	5470	17179
# 1000	6506	6583	7299	48242
# sum(a,2)
# 	1	10	100	1000
# 1	4690	4691	4691	6773
# 10	4690	4750	4884	8085
# 100	4929	5197	6866	26475
# 1000	13841	15303	32141	216153

# COMMIT: 89c2026 2017-04-02 (buggy with large arrays)
# sum(a,0)
# 	1	10	100	1000
# 1	32326	32417	32342	32467
# 10	32390	32346	32232	32597
# 100	32212	32345	32571	35466
# 1000	32363	32547	35519	66645
# sum(a,1)
# 	1	10	100	1000
# 1	4470	4528	4777	10954
# 10	4713	4777	5142	12531
# 100	4699	4778	5200	12915
# 1000	6245	6435	7126	42014
# sum(a,2)
# 	1	10	100	1000
# 1	4480	4481	4483	6546
# 10	4541	4601	4695	7929
# 100	4777	4975	6781	26433
# 1000	10970	12016	28019	208529


text:
using KernelAbstractions, Test, CUDA

if has_cuda_gpu()
    CUDA.allowscalar(false)
end

# Simple kernel for matrix multiplication
@kernel function vadd_kernel!(c, a, b)
    i, j = @index(Global, NTuple)

    c[i,j] = a[i,j] + b[i,j]
end

# Creating a wrapper kernel for launching with error checks
function vadd!(a, b, c)
    if isa(a, Array)
        kernel! = vadd_kernel!(CPU(),4)
    else
        kernel! = vadd_kernel!(CUDADevice(),256)
    end
    kernel!(c, a, b, ndrange=size(c))
end

a = ones(1024, 1024)
b = ones(1024, 1024)
c = zeros(1024, 1024)

# beginning CPU tests, returns event
ev = vadd!(a,b,c)
wait(ev)

@test isapprox(c, a.+b)

# beginning GPU tests
d_a = CuArray(a)
d_b = CuArray(b)
d_c = CuArray(c)

ev = vadd!(d_a, d_b, d_c)
wait(ev)

@test isapprox(Array(d_c), a.+b)


text:
using UnitCellPlotter
using Test

@testset "UnitCellPlotter.jl" begin
    # Write your tests here.
end


text:
using StaticArrays
using LinearAlgebra
using LoopVectorization
using Distributed
include("constants.jl")
include("smear.jl")

# get rid of all the wasted space with the α and damping_fac and make SVector
mutable struct Electrostatics
    natom::Int
    q::Vector{Float64} # 4 * nW
    dipoles::Vector{MVector{3, Float64}} # 3 * natom
    damping_fac::SVector{4, Float64} # 4
    α::SVector{4, Float64} # 4
    ϕ::Vector{Float64} # 4 * nW
    E_field_q::Vector{MVector{3, Float64}} # 3 * natom elements
    E_field_dip::Vector{MVector{3, Float64}} # 3 * natom elements
    dip_dip_tensor::Matrix{SMatrix{3, 3, Float64, 9}} # natom * natom * 3 * 3 tensor
    previous_dipoles::Vector{MVector{3, Float64}} # 3 * natom
    smear::Smear # smearing method
    dipoles_converged::Bool # convergence option
end

Electrostatics(charges::Vector{Float64}, C::TTM_Constants) =
Electrostatics(length(charges), charges, [@MVector zeros(3) for _ in  1:length(charges)], [C.damping_factor_O, C.damping_factor_H, C.damping_factor_H, C.damping_factor_M], [C.α_O, C.α_H, C.α_H, C.α_M], zeros(length(charges)), [@MVector zeros(3) for _ in 1:length(charges)],
[@MVector zeros(3) for _ in 1:length(charges)], vecvec_to_matrix2([[@MMatrix zeros(3,3) for _ in 1:length(charges)] for _ in 1:length(charges)]), [@MVector zeros(3) for _ in  1:length(charges)], get_smear_type(C.name), false)

function vecvec_to_matrix2(vecvec::AbstractVector{T}) where T <: AbstractVector
        dim1 = length(vecvec)
        dim2 = length(vecvec[1])
        my_array = Array{eltype(vecvec[1]), 2}(undef, dim1, dim2)
        @inbounds for i in 1:dim1, j in 1:dim2
        my_array[i,j] = vecvec[i][j]
    end
    return my_array
end

@inline function index_mod_4(i::Int)
    return ((i-1) % 4) + 1
end

function reset_electrostatics!(elec_data::Electrostatics)
    """
    Resets the electrostatics struct so as to not have problems repeating the calculations
    with the same container.
    """
    elec_data.ϕ -= elec_data.ϕ
    elec_data.E_field_q -= elec_data.E_field_q
    elec_data.E_field_dip -= elec_data.E_field_dip
    elec_data.dip_dip_tensor -= elec_data.dip_dip_tensor
end

function get_α_sqrt(α::Vector{Float64})
    α_sqrt = zeros(4 * length(α))
    for i in 1:length(α)
        for j in 1:3
            @inbounds α_sqrt[3*(i-1)+j] = sqrt(α[i])
        end
    end
    return α_sqrt
end

@inline function skip_ij(i::Int, j::Int, only_M_sites::Bool=false)
    """
    Determine if i and j are a part of the same water.
    Since everything has to be in OHHM order, we can just infer from the indices.
    Note that this won't check if your input is valid (all indices should be >= 1).
    We also provide an option for filtering based on if the M-site is the only polarizable considered. This is the case in the TTM3-F potential.
    """
    if only_M_sites
        return (i % 4 != 0) || (j % 4 != 0) || ((i-1) ÷ 4 == (j-1) ÷ 4)
    else
        return (i-1) ÷ 4 == (j-1) ÷ 4
    end
end

function form_dipole_dipole_tensor!(elec_data::Electrostatics, coords::Vector{SVector{3, Float64}}, only_M_sites::Bool)
    atom_stride::Int = 1
    if only_M_sites
        atom_stride = 4
    end

    for i in atom_stride:atom_stride:(elec_data.natom-1)
        i3::Int = 3 * (i - 1) + 1
        for j in (i+atom_stride):atom_stride:elec_data.natom
            j3::Int = 3 * (j - 1) + 1
            @views r_ij = coords[i] - coords[j]

            aDD::Float64 = skip_ij(i, j) ? elec_data.smear.aDD_intramolecular : elec_data.smear.aDD_intermolecular

            ts1, ts2 = elec_data.smear.smear2(r_ij ⋅ r_ij, elec_data.damping_fac[index_mod_4(i)] * elec_data.damping_fac[index_mod_4(j)], aDD)
            elec_data.dip_dip_tensor[i, j]::SMatrix{3, 3, Float64, 9} = 3.0 * ts2 * r_ij * r_ij' - diagm(ones(3)) * ts1
            elec_data.dip_dip_tensor[j, i] = elec_data.dip_dip_tensor[i, j]
        end
    end
end

function charge_charge_interactions!(elec_data::Electrostatics, coords::Vector{SVector{3, Float64}})
    for i in 1:(elec_data.natom-1)
        for j in (i+1):elec_data.natom
            # charge-charge
            if !skip_ij(i, j)
                @views r_ij = coords[i] - coords[j]

                ts0, ts1 = elec_data.smear.smear1(r_ij ⋅ r_ij, elec_data.damping_fac[index_mod_4(i)] * elec_data.damping_fac[index_mod_4(j)], elec_data.smear.aCC)

                @inbounds elec_data.ϕ[i] += ts0 * elec_data.q[j]
                @inbounds elec_data.ϕ[j] += ts0 * elec_data.q[i]

                @inbounds elec_data.E_field_q[i] += ts1 * elec_data.q[j] * r_ij
                @inbounds elec_data.E_field_q[j] -= ts1 * elec_data.q[i] * r_ij
            end
        end
    end
end

function charge_charge_gradients!(elec_data::Electrostatics, grads_E::Vector{MVector{3, Float64}})
    for i in 1:elec_data.natom
        @inbounds grads_E[i] = -elec_data.q[i] * elec_data.E_field_q[i]
    end
end

function charge_dipole_gradients!(elec_data::Electrostatics, coords::Vector{SVector{3, Float64}}, grads_E::Vector{MVector{3, Float64}}, only_M_sites::Bool)
    for i in 1:elec_data.natom
        if only_M_sites
            if (i % 4) != 0
                continue
            end
        elseif !(!only_M_sites && (i % 4) != 0)
            continue
        end
        for j in 1:elec_data.natom
            if (skip_ij(i, j))
                continue # skip this (i, j) pair
            end
            Rij = coords[i] - coords[j]
            diR = elec_data.dipoles[i] ⋅ Rij

            @inbounds ts1, ts2 = elec_data.smear.smear2(Rij ⋅ Rij, elec_data.damping_fac[index_mod_4(i)] * elec_data.damping_fac[index_mod_4(j)], elec_data.smear.aCD)

            if (j % 4) != 1
                @inbounds derij = elec_data.q[j] * (3 * ts2 * diR * Rij - ts1 * elec_data.dipoles[i])

                @inbounds grads_E[i] += derij
                @inbounds grads_E[j] -= derij
            end

            elec_data.ϕ[j] -= ts1 * diR
        end
    end
end

function dipole_dipole_gradients!(elec_data::Electrostatics, coords::Vector{SVector{3, Float64}}, grads_E::Vector{MVector{3, Float64}}, only_M_sites::Bool, atom_stride::Int)
    for i in atom_stride:atom_stride:elec_data.natom-1
        if (i-1) % 4 == 3 && only_M_sites
            continue
        end
        for j in i+atom_stride:atom_stride:elec_data.natom
            if (j-1) % 4 == 3
                continue
            end
            Rij = coords[i] - coords[j]
            @inbounds diR  = elec_data.dipoles[i] ⋅ Rij
            @inbounds djR  = elec_data.dipoles[j] ⋅ Rij
            @inbounds didj = elec_data.dipoles[i] ⋅ elec_data.dipoles[j]

            aDD = skip_ij(i, j) ? elec_data.smear.aDD_intramolecular : elec_data.smear.aDD_intermolecular

            _, ts2, ts3 = elec_data.smear.smear3(Rij ⋅ Rij, elec_data.damping_fac[index_mod_4(i)] * elec_data.damping_fac[index_mod_4(j)], aDD)

            @inbounds derij = - 3 * ts2 * (didj * Rij + djR * elec_data.dipoles[i] + diR * elec_data.dipoles[j]) + 15 * ts3 * diR * djR * Rij
            @inbounds grads_E[i] += derij
            @inbounds grads_E[j] -= derij
        end
    end
end

# OPTIMIZE EVERYTHING HERE
# get the properly sized matrices
# use static vectors wherever possible
# where is all the allocation happening??
# get rid of all the intermediates and loops over x,y,z
#
# Some of the interactions are completely independent. For instance, try separating the charge-charge interactions and formation of the ddt into separate tasks handled by their own threads.
function electrostatics(elec_data::Electrostatics, coords::Vector{SVector{3, Float64}}, grads_E::Union{Vector{MVector{3, Float64}}, Nothing}=nothing, only_M_sites::Bool=false, dipole_tolerance::Float64=1.0e-15, dipole_maxiter::Int=5000)
    reset_electrostatics!(elec_data)

    atom_stride::Int = 1
    if only_M_sites
        atom_stride = 4
    end

    # @SPEED make these separate tasks that run independently
    # and then join afterwards
    charge_charge_interactions!(elec_data, coords)
    form_dipole_dipole_tensor!(elec_data, coords, only_M_sites)

    # This seems to be a not much better guess than just
    # starting from zeros.
    # Using the dipoles from last call to the potential is the best
    # when one can be sure the dipoles are similar from call to call.
    # Perhaps there should be "MD/MC Mode" to hint that the
    # dipoles will be similar from call to call.

    # Calculate induced dipoles iteratively
    #for i in 1:elec_data.natom
    #    @inbounds elec_data.dipoles[i] = elec_data.α[1:3] .* elec_data.E_field_q[i]
    #    @inbounds elec_data.previous_dipoles[i] = elec_data.α[1:3] .* elec_data.E_field_q[i]
    #end

    dmix::Float64 = 0.75
    stath::Float64 = debye() / chargecon() / sqrt(elec_data.natom)

    elec_data.dipoles_converged = false

    for iter in 1:dipole_maxiter
        for k in atom_stride:atom_stride:elec_data.natom
            elec_data.E_field_dip[k] = @distributed (+) for l in atom_stride:atom_stride:elec_data.natom
                @inbounds @views elec_data.dip_dip_tensor[k, l] * elec_data.dipoles[l]
            end
        end

        deltadip::Float64 = 0.0

        for i in atom_stride:atom_stride:elec_data.natom
            @inbounds elec_data.dipoles[i] = elec_data.α[index_mod_4(i)] * (elec_data.E_field_q[i] + elec_data.E_field_dip[i])
        end

        @inbounds elec_data.dipoles = dmix * elec_data.dipoles + (1.0 - dmix) * elec_data.previous_dipoles
        delta = sum(sum.(elec_data.dipoles - elec_data.previous_dipoles))

        deltadip = sqrt(delta^2) * stath

        if deltadip < dipole_tolerance
            elec_data.dipoles_converged = true
            break # converged!
        else
            copy!.(elec_data.previous_dipoles, elec_data.dipoles)
            elec_data.E_field_dip -= elec_data.E_field_dip
        end
    end

    E_elec = 0.5 * (elec_data.ϕ ⋅ elec_data.q)

    E_ind::Float64 = 0.0
    for i in atom_stride:atom_stride:elec_data.natom
        @inbounds E_ind -= elec_data.dipoles[i] ⋅ elec_data.E_field_q[i]
    end
    E_ind *= 0.5

    # if no gradients requested, then return the energy here
    if grads_E === nothing
        return E_elec + E_ind
    end

    ### CALCULATE THE GRADIENTS ###
    charge_charge_gradients!(elec_data, grads_E)
    charge_dipole_gradients!(elec_data, coords, grads_E, only_M_sites)
    dipole_dipole_gradients!(elec_data, coords, grads_E, only_M_sites, atom_stride)

    return E_elec + E_ind
end

max_stars_repo_name: heindelj/MolecularPotentials

max_issues_repo_name: heindelj/MolecularPotentials

max_forks_repo_name: heindelj/MolecularPotentials


text:
"""
    Channel

Struct containing asymptotic kernels for a channel.
* `q1   :: Matrix{Float64}`   : kernel with both fermionic frequencies -> Inf
* `q2_1 :: Array{Float64, 3}` : kernel with second fermionic frequency -> Inf
* `q2_2 :: Array{Float64, 3}` : kernel with first fermionic frequency -> Inf
* `q3   :: Array{Float64, 4}` : full channel
"""
struct Channel
    q1   :: Matrix{Float64}
    q2_1 :: Array{Float64, 3}
    q2_2 :: Array{Float64, 3}
    q3   :: Array{Float64, 4}
end

# generate channel dummy
function get_channel_empty(
    r :: Reduced_lattice,
    m :: Mesh
    ) :: Channel

    num_sites = length(r.sites)

    # init kernels
    q1   = zeros(Float64, num_sites, m.num_Ω)
    q2_1 = zeros(Float64, num_sites, m.num_Ω, m.num_ν)
    q2_2 = zeros(Float64, num_sites, m.num_Ω, m.num_ν)
    q3   = zeros(Float64, num_sites, m.num_Ω, m.num_ν, m.num_ν)

    # build channel
    ch = Channel(q1, q2_1, q2_2, q3)

    return ch
end

# get interpolated value of q1
function get_q1(
    site :: Int64,
    p    :: Param,
    ch   :: Channel
    )    :: Float64

    val  = p.lower_weight * ch.q1[site, p.lower_index]
    val += p.upper_weight * ch.q1[site, p.upper_index]

    return val
end

# get interpolated value of q1 on all lattice sites
function get_q1_avx!(
    r        :: Reduced_lattice,
    p        :: Param,
    ch       :: Channel,
    temp     :: SubArray{Float64, 1, Array{Float64, 3}},
    exchange :: Bool,
    sgn      :: Float64
    )        :: Nothing

    # deref channel
    q1 = ch.q1

    # deref param
    lower_index  = p.lower_index
    upper_index  = p.upper_index
    lower_weight = p.lower_weight
    upper_weight = p.upper_weight

    # check for site exchange
    if exchange
        indices = r.exchange

        @turbo unroll = 1 for i in eachindex(temp)
            val      = lower_weight * q1[indices[i], lower_index]
            val     += upper_weight * q1[indices[i], upper_index]
            temp[i] += sgn * val
        end
    else
        @turbo unroll = 1 for i in eachindex(temp)
            val      = lower_weight * q1[i, lower_index]
            val     += upper_weight * q1[i, upper_index]
            temp[i] += sgn * val
        end
    end

    return nothing
end

# get interpolated value of q2_1
function get_q2_1(
    site :: Int64,
    p1   :: Param,
    p2   :: Param,
    ch   :: Channel
    )    :: Float64

    val  = p1.lower_weight * p2.lower_weight * ch.q2_1[site, p1.lower_index, p2.lower_index]
    val += p1.upper_weight * p2.lower_weight * ch.q2_1[site, p1.upper_index, p2.lower_index]
    val += p1.lower_weight * p2.upper_weight * ch.q2_1[site, p1.lower_index, p2.upper_index]
    val += p1.upper_weight * p2.upper_weight * ch.q2_1[site, p1.upper_index, p2.upper_index]

    return val
end

# get interpolated value of q2_1 on all lattice sites
function get_q2_1_avx!(
    r        :: Reduced_lattice,
    p1       :: Param,
    p2       :: Param,
    ch       :: Channel,
    temp     :: SubArray{Float64, 1, Array{Float64, 3}},
    exchange :: Bool,
    sgn      :: Float64
    )        :: Nothing

    # deref channel
    q2_1 = ch.q2_1

    # deref param
    lower_index1  = p1.lower_index
    upper_index1  = p1.upper_index
    lower_weight1 = p1.lower_weight
    upper_weight1 = p1.upper_weight

    lower_index2  = p2.lower_index
    upper_index2  = p2.upper_index
    lower_weight2 = p2.lower_weight
    upper_weight2 = p2.upper_weight

    # check for site exchange
    if exchange
        indices = r.exchange

        @turbo unroll = 1 for i in eachindex(temp)
            val      = lower_weight1 * lower_weight2 * q2_1[indices[i], lower_index1, lower_index2]
            val     += upper_weight1 * lower_weight2 * q2_1[indices[i], upper_index1, lower_index2]
            val     += lower_weight1 * upper_weight2 * q2_1[indices[i], lower_index1, upper_index2]
            val     += upper_weight1 * upper_weight2 * q2_1[indices[i], upper_index1, upper_index2]
            temp[i] += sgn * val
        end
    else
        @turbo unroll = 1 for i in eachindex(temp)
            val      = lower_weight1 * lower_weight2 * q2_1[i, lower_index1, lower_index2]
            val     += upper_weight1 * lower_weight2 * q2_1[i, upper_index1, lower_index2]
            val     += lower_weight1 * upper_weight2 * q2_1[i, lower_index1, upper_index2]
            val     += upper_weight1 * upper_weight2 * q2_1[i, upper_index1, upper_index2]
            temp[i] += sgn * val
        end
    end

    return nothing
end

# get interpolated value of q2_2
function get_q2_2(
    site :: Int64,
    p1   :: Param,
    p2   :: Param,
    ch   :: Channel
    )    :: Float64

    val  = p1.lower_weight * p2.lower_weight * ch.q2_2[site, p1.lower_index, p2.lower_index]
    val += p1.upper_weight * p2.lower_weight * ch.q2_2[site, p1.upper_index, p2.lower_index]
    val += p1.lower_weight * p2.upper_weight * ch.q2_2[site, p1.lower_index, p2.upper_index]
    val += p1.upper_weight * p2.upper_weight * ch.q2_2[site, p1.upper_index, p2.upper_index]

    return val
end

# get interpolated value of q2_2 on all lattice sites
function get_q2_2_avx!(
    r        :: Reduced_lattice,
    p1       :: Param,
    p2       :: Param,
    ch       :: Channel,
    temp     :: SubArray{Float64, 1, Array{Float64, 3}},
    exchange :: Bool,
    sgn      :: Float64
    )        :: Nothing

    # deref channel
    q2_2 = ch.q2_2

    # deref param
    lower_index1  = p1.lower_index
    upper_index1  = p1.upper_index
    lower_weight1 = p1.lower_weight
    upper_weight1 = p1.upper_weight

    lower_index2  = p2.lower_index
    upper_index2  = p2.upper_index
    lower_weight2 = p2.lower_weight
    upper_weight2 = p2.upper_weight

    # check for site exchange
    if exchange
        indices = r.exchange

        @turbo unroll = 1 for i in eachindex(temp)
            val      = lower_weight1 * lower_weight2 * q2_2[indices[i], lower_index1, lower_index2]
            val     += upper_weight1 * lower_weight2 * q2_2[indices[i], upper_index1, lower_index2]
            val     += lower_weight1 * upper_weight2 * q2_2[indices[i], lower_index1, upper_index2]
            val     += upper_weight1 * upper_weight2 * q2_2[indices[i], upper_index1, upper_index2]
            temp[i] += sgn * val
        end
    else
        @turbo unroll = 1 for i in eachindex(temp)
            val      = lower_weight1 * lower_weight2 * q2_2[i, lower_index1, lower_index2]
            val     += upper_weight1 * lower_weight2 * q2_2[i, upper_index1, lower_index2]
            val     += lower_weight1 * upper_weight2 * q2_2[i, lower_index1, upper_index2]
            val     += upper_weight1 * upper_weight2 * q2_2[i, upper_index1, upper_index2]
            temp[i] += sgn * val
        end
    end

    return nothing
end

# get interpolated value of q3
function get_q3(
    site :: Int64,
    p1   :: Param,
    p2   :: Param,
    p3   :: Param,
    ch   :: Channel
    )    :: Float64

    val  = p1.lower_weight * p2.lower_weight * p3.lower_weight * ch.q3[site, p1.lower_index, p2.lower_index, p3.lower_index]
    val += p1.upper_weight * p2.lower_weight * p3.lower_weight * ch.q3[site, p1.upper_index, p2.lower_index, p3.lower_index]
    val += p1.lower_weight * p2.upper_weight * p3.lower_weight * ch.q3[site, p1.lower_index, p2.upper_index, p3.lower_index]
    val += p1.upper_weight * p2.upper_weight * p3.lower_weight * ch.q3[site, p1.upper_index, p2.upper_index, p3.lower_index]
    val += p1.lower_weight * p2.lower_weight * p3.upper_weight * ch.q3[site, p1.lower_index, p2.lower_index, p3.upper_index]
    val += p1.upper_weight * p2.lower_weight * p3.upper_weight * ch.q3[site, p1.upper_index, p2.lower_index, p3.upper_index]
    val += p1.lower_weight * p2.upper_weight * p3.upper_weight * ch.q3[site, p1.lower_index, p2.upper_index, p3.upper_index]
    val += p1.upper_weight * p2.upper_weight * p3.upper_weight * ch.q3[site, p1.upper_index, p2.upper_index, p3.upper_index]

    return val
end

# get interpolated value of q3 on all lattice sites
function get_q3_avx!(
    r        :: Reduced_lattice,
    p1       :: Param,
    p2       :: Param,
    p3       :: Param,
    ch       :: Channel,
    temp     :: SubArray{Float64, 1, Array{Float64, 3}},
    exchange :: Bool,
    sgn      :: Float64
    )        :: Nothing

    # deref channel
    q3 = ch.q3

    # deref param
    lower_index1  = p1.lower_index
    upper_index1  = p1.upper_index
    lower_weight1 = p1.lower_weight
    upper_weight1 = p1.upper_weight

    lower_index2  = p2.lower_index
    upper_index2  = p2.upper_index
    lower_weight2 = p2.lower_weight
    upper_weight2 = p2.upper_weight

    lower_index3  = p3.lower_index
    upper_index3  = p3.upper_index
    lower_weight3 = p3.lower_weight
    upper_weight3 = p3.upper_weight

    # check for site exchange
    if exchange
        indices = r.exchange

        @turbo unroll = 1 for i in eachindex(temp)
            val      = lower_weight1 * lower_weight2 * lower_weight3 * q3[indices[i], lower_index1, lower_index2, lower_index3]
            val     += upper_weight1 * lower_weight2 * lower_weight3 * q3[indices[i], upper_index1, lower_index2, lower_index3]
            val     += lower_weight1 * upper_weight2 * lower_weight3 * q3[indices[i], lower_index1, upper_index2, lower_index3]
            val     += upper_weight1 * upper_weight2 * lower_weight3 * q3[indices[i], upper_index1, upper_index2, lower_index3]
            val     += lower_weight1 * lower_weight2 * upper_weight3 * q3[indices[i], lower_index1, lower_index2, upper_index3]
            val     += upper_weight1 * lower_weight2 * upper_weight3 * q3[indices[i], upper_index1, lower_index2, upper_index3]
            val     += lower_weight1 * upper_weight2 * upper_weight3 * q3[indices[i], lower_index1, upper_index2, upper_index3]
            val     += upper_weight1 * upper_weight2 * upper_weight3 * q3[indices[i], upper_index1, upper_index2, upper_index3]
            temp[i] += sgn * val
        end
    else
        @turbo unroll = 1 for i in eachindex(temp)
            val      = lower_weight1 * lower_weight2 * lower_weight3 * q3[i, lower_index1, lower_index2, lower_index3]
            val     += upper_weight1 * lower_weight2 * lower_weight3 * q3[i, upper_index1, lower_index2, lower_index3]
            val     += lower_weight1 * upper_weight2 * lower_weight3 * q3[i, lower_index1, upper_index2, lower_index3]
            val     += upper_weight1 * upper_weight2 * lower_weight3 * q3[i, upper_index1, upper_index2, lower_index3]
            val     += lower_weight1 * lower_weight2 * upper_weight3 * q3[i, lower_index1, lower_index2, upper_index3]
            val     += upper_weight1 * lower_weight2 * upper_weight3 * q3[i, upper_index1, lower_index2, upper_index3]
            val     += lower_weight1 * upper_weight2 * upper_weight3 * q3[i, lower_index1, upper_index2, upper_index3]
            val     += upper_weight1 * upper_weight2 * upper_weight3 * q3[i, upper_index1, upper_index2, upper_index3]
            temp[i] += sgn * val
        end
    end

    return nothing
end

# get interpolated value of channel for a given frequency buffer
function get_channel(
    site :: Int64,
    b    :: Buffer,
    ch   :: Channel
    )    :: Float64

    val = 0.0

    if b.kernel == 1
        val = get_q1(site, b.p1, ch)
    elseif b.kernel == 2
        val = get_q2_1(site, b.p1, b.p2, ch)
    elseif b.kernel == 3
        val = get_q2_2(site, b.p1, b.p3, ch)
    elseif b.kernel == 4
        val = get_q3(site, b.p1, b.p2, b.p3, ch)
    end

    return val
end

# get interpolated value of channel for a given frequency buffer on all lattice sites
function get_channel_avx!(
    r        :: Reduced_lattice,
    b        :: Buffer,
    ch       :: Channel,
    temp     :: SubArray{Float64, 1, Array{Float64, 3}},
    exchange :: Bool,
    sgn      :: Float64
    )        :: Nothing

    if b.kernel == 1
        get_q1_avx!(r, b.p1, ch, temp, exchange, sgn)
    elseif b.kernel == 2
        get_q2_1_avx!(r, b.p1, b.p2, ch, temp, exchange, sgn)
    elseif b.kernel == 3
        get_q2_2_avx!(r, b.p1, b.p3, ch, temp, exchange, sgn)
    elseif b.kernel == 4
        get_q3_avx!(r, b.p1, b.p2, b.p3, ch, temp, exchange, sgn)
    end

    return nothing
end

# replace channel with another channel
function replace_with!(
    ch1 :: Channel,
    ch2 :: Channel
    )   :: Nothing

    @turbo ch1.q1   .= ch2.q1
    @turbo ch1.q2_1 .= ch2.q2_1
    @turbo ch1.q2_2 .= ch2.q2_2
    @turbo ch1.q3   .= ch2.q3

    return nothing
end

# multiply channel with factor
function mult_with!(
    ch  :: Channel,
    fac :: Float64
    )   :: Nothing

    @turbo ch.q1   .*= fac
    @turbo ch.q2_1 .*= fac
    @turbo ch.q2_2 .*= fac
    @turbo ch.q3   .*= fac

    return nothing
end

# multiply channel with factor and add to other channel
function mult_with_add_to!(
    ch2 :: Channel,
    fac :: Float64,
    ch1 :: Channel
    )   :: Nothing

    @turbo ch1.q1   .+= fac .* ch2.q1
    @turbo ch1.q2_1 .+= fac .* ch2.q2_1
    @turbo ch1.q2_2 .+= fac .* ch2.q2_2
    @turbo ch1.q3   .+= fac .* ch2.q3

    return nothing
end

"""
    get_abs_max(
        ch :: Channel
        )  :: Float64

Returns maximum absolute value of a channel.
"""
function get_abs_max(
    ch :: Channel
    )  :: Float64

    max_ch = maximum(abs.(ch.q3))

    return max_ch
end

# resample a channel to new meshes via trilinear interpolation
function resample_from_to!(
    Ω_old  :: Vector{Float64},
    ν_old  :: Vector{Float64},
    ch_old :: Channel,
    Ω_new  :: Vector{Float64},
    ν_new  :: Vector{Float64},
    ch_new :: Channel
    )      :: Nothing

    # get dimensions
    num_sites = size(ch_old.q2_1, 1)
    num_Ω     = size(ch_old.q2_1, 2)
    num_ν     = size(ch_old.q2_1, 3)

    # resample q1
    for w in 1 : num_Ω
        # get interpolation parameters
        p1 = get_param(Ω_new[w], Ω_old)

        for site in 1 : num_sites
            ch_new.q1[site, w] = get_q1(site, p1, ch_old)
        end
    end

    # resample q2_1 and q2_2
    for v in 1 : num_ν
        for w in 1 : num_Ω
            # get interpolation parameters
            p1 = get_param(Ω_new[w], Ω_old)
            p2 = get_param(ν_new[v], ν_old)

            for site in 1 : num_sites
                ch_new.q2_1[site, w, v] = get_q2_1(site, p1, p2, ch_old)
                ch_new.q2_2[site, w, v] = get_q2_2(site, p1, p2, ch_old)
            end
        end
    end

    # resample q3
    for vp in 1 : num_ν
        for v in 1 : num_ν
            for w in 1 : num_Ω
                # get interpolation parameters
                p1 = get_param(Ω_new[w],  Ω_old)
                p2 = get_param(ν_new[v],  ν_old)
                p3 = get_param(ν_new[vp], ν_old)

                for site in 1 : num_sites
                    ch_new.q3[site, w, v, vp] = get_q3(site, p1, p2, p3, ch_old)
                end
            end
        end
    end

    return nothing
end


text:
function fields(M::T, r, z) where T<:AbstractEquilibrium
    B = Bfield(M, r, z)
    E = Efield(M, r, z)
    return B, E
end

function fields(M::AbstractEquilibrium, x, y, z)
    r = hypot(x,y)
    phi = atan(y,x)
    B, E = fields(M,r,z)
    sp, cp = sincos(phi)
    B = SVector{3}(B[1]*cp - B[2]*sp, B[1]*sp + B[2]*cp, B[3])
    E = SVector{3}(E[1]*cp - E[2]*sp, E[1]*sp + E[2]*cp, E[3])
    return B, E
end

function Bfield(M::T, r, z) where T<:AbstractEquilibrium
    psi = M(r,z)
    gval = poloidal_current(M,psi)
    grad_psi = psi_gradient(M,r,z)

    # Calculate B-Field
    cc = cocos(M)
    invr = inv(r)
    #invr not part of cocos factor but added here for performance
    cocos_factor = invr*cc.sigma_RpZ*cc.sigma_Bp/((2pi)^cc.exp_Bp)

    BR =  cocos_factor*grad_psi[2]
    Bz = -cocos_factor*grad_psi[1]
    Bt = gval*invr

    B = SVector{3}(cylindrical_cocos(cc, BR, Bt, Bz))
    return B
end

function Bfield(M::T, x, y, z) where T<:AbstractEquilibrium
    r = hypot(x,y)
    phi = atan(y,x)
    B = Bfield(M,r,z)
    sp, cp = sincos(phi)
    B_xyz = SVector{3}(B[1]*cp - B[2]*sp, B[1]*sp + B[2]*cp, B[3])
    return B_xyz
end

function poloidal_Bfield(M::T, r, z) where T<:AbstractEquilibrium
    rmaxis, zmaxis = magnetic_axis(M)
    cc = cocos(M)
    B = Bfield(M, r, z)

    ir, iphi, iz = cylindrical_cocos_indices(cc)
    sign_theta = cc.sigma_RpZ*cc.sigma_rhotp # + CW, - CCW
    sign_Bp    = sign_theta*sign((z-zmaxis)*B[ir] - (r-rmaxis)*B[iz]) # sign(theta)*sign(r x B)
    Bpol = sign_Bp*sqrt(B[ir]^2 + B[iz]^2)
    return Bpol
end

function Jfield(M::T, r, z; curl=false) where T<:AbstractEquilibrium
    if curl
        return curlB(M, r, z)/mu0
    end

    psi = M(r,z)
    gval = poloidal_current(M,psi)
    grad_psi = psi_gradient(M,r,z)

    gp = poloidal_current_gradient(M,psi)
    pp = pressure_gradient(M,psi)

    cc = cocos(M)

    invrmu0 = inv(r*mu0)
    cocos_factor = cc.sigma_RpZ*invrmu0
    jr = -cocos_factor*gp*grad_psi[2]
    jz =  cocos_factor*gp*grad_psi[1]

    cocos_factor = -cc.sigma_Bp*((2pi)^cc.exp_Bp)
    jt = cocos_factor*(r*pp + gval*gp*invrmu0)

    return SVector{3}(cylindrical_cocos(cc,jr,jt,jz))
end

function Jfield(M::T, x, y, z; kwargs...) where T<:AbstractEquilibrium
    r = hypot(x,y)
    phi = atan(y,x)
    J = Jfield(M,r,z; kwargs...)
    sp, cp = sincos(phi)
    J_xyz = SVector{3}(J[1]*cp - J[2]*sp, J[1]*sp + J[2]*cp, J[3])
    return J_xyz
end

function poloidal_Jfield(M::T, r, z) where T<:AbstractEquilibrium
    rmaxis, zmaxis = magnetic_axis(M)
    cc = cocos(M)
    B = Jfield(M, r, z)

    ir, iphi, iz = cylindrical_cocos_indices(cc)
    sign_theta = cc.sigma_RpZ*cc.sigma_rhotp # + CW, - CCW
    sign_Jp    = sign_theta*sign((z-zmaxis)*J[ir] - (r-rmaxis)*J[iz]) # sign(theta)*sign(r x B)
    Jpol = sign_Jp*sqrt(J[ir]^2 + J[iz]^2)
    return Jpol
end

function Efield(M::T, r, z) where T<:AbstractEquilibrium

    psi = M(r,z)
    phi_grad = electric_potential_gradient(M, psi)
    if phi_grad == 0
        return SVector{3}(zero(r),zero(r),zero(r))
    end

    Bpol = poloidal_Bfield(M, r, z) #signed

    grad_psi = psi_gradient(M,r,z)
    grad_psi = grad_psi/norm(grad_psi)

    # is Bpol signed in this equation?
    Er = -r*Bpol*phi_grad

    ER = Er*grad_psi[1] # Er*dpsi/dR = (-dphi/dpsi)*(dpsi/dR)
    Ez = Er*grad_psi[2] # Er*dpsi/dz = (-dphi/dpsi)*(dpsi/dz)
    Et = zero(Ez)

    return SVector{3}(cylindrical_cocos(cocos(M), ER, Et, Ez))
end

function Efield(M::T, x, y, z) where T<:AbstractEquilibrium
    r = hypot(x,y)
    phi = atan(y,x)
    E = Efield(M,r,z)
    sp, cp = sincos(phi)
    E_xyz = SVector{3}(E[1]*cp - E[2]*sp, E[1]*sp + E[2]*cp, E[3])
    return E_xyz
end

function Efield(M::T, r, z, vrot::AbstractVector) where T<:AbstractEquilibrium
    cc = cocos(M)
    psi = M(r,z)
    gval = poloidal_current(M,psi)
    B = Bfield(M, r, z)

    grad_psi = psi_gradient(M,r,z)
    grad_psi = grad_psi/norm(grad_psi)

    qval = safety_factor(M,psi)
    dp = pressure_gradient(M, psi)

    gradp = SVector{3}(cylindrical_cocos(cc, dp*grad_psi[1], zero(grad_psi[1]), dp*grad_psi[2]))

    E = cross(vrot,B) .+ gradp/qval

    return E
end

function Efield(M::T, x, y, z, vrot::AbstractVector) where T<:AbstractEquilibrium
    r = hypot(x,y)
    phi = atan(y,x)
    E = Efield(M,r,z,vrot)
    sp, cp = sincos(phi)
    E_xyz = SVector{3}(E[1]*cp - E[2]*sp, E[1]*sp + E[2]*cp, E[3])
    return E_xyz
end

function rho_p(M::T, psi) where T<:AbstractEquilibrium
    psimag,psibry = psi_limits(M)
    flux = abs(psibry - psimag)
    return sqrt((psi - psimag)/flux)
end

function rho_p(M::T, r, z) where T<:AbstractEquilibrium
    rho_p(M,M(r,z))
end

function rho_p(M::T, x, y, z) where T<:AbstractEquilibrium
    r = hypot(x,y)
    rho_p(M,r,z)
end

function toroidal_flux(M::T,psi) where T<:AbstractEquilibrium
    cc = cocos(M)
    cocos_factor = cc.sigma_Bp*cc.sigma_rhotp*(2pi)^(1 - cc.exp_Bp)

    psimag, psibry = psi_limits(M)
    phi = cocos_factor*hquadrature(x->safety_factor(M,x), psimag, psi)[1]
    return phi
end

function safety_factor(M::AbstractEquilibrium, psi)
    g = poloidal_current(M,psi)
    fs = boundary(M, psi)
    q = (g/(2pi))*average(fs,(x,y)->inv(poloidal_Bfield(M,x,y)*x^2))*circumference(fs)
    return q
end

function plasma_current(M::AbstractEquilibrium)
    cc = cocos(M)

    bdry = boundary(M)
    dr = diff(bdry.r)
    dz = diff(bdry.z)

    ll = vcat(0,cumsum(hypot.(dr, dz)))

    Bp = poloidal_Bfield.(M, bdry.r, bdry.z)

    Ip = cc.sigma_rhotp * trapz(ll, Bp) / mu0

    return Ip
end

function beta_n(M::AbstractEquilibrium)
    Ip = plasma_current(M)*1e-6
    bdry = boundary(M)
    rmin, rmax = extrema(bdry.r)
    a = (rmax - rmin)/2
    beta_n = M.beta_t / abs(Ip/a/M.B0) * 100
    return beta_n
end

function gradB_autodiff(M::AbstractEquilibrium, r, z)
    gB_rz = ForwardDiff.gradient(x->norm(Bfield(M,x[1],x[2])), SVector{2}(r,z))
    return SVector{3}(cylindrical_cocos(cocos(M), gB_rz[1], 0.0, gB_rz[2]))
end

function gradB_autodiff(M::AbstractEquilibrium, x, y, z)
    gB_xyz = ForwardDiff.gradient(x->norm(Bfield(M,x[1],x[2],x[3])), SVector{3}(x,y,z))
    return SVector{3}(gB_xyz[1], gB_xyz[2], gB_xyz[3])
end

gradB(M::AbstractEquilibrium, r, z) = gradB_autodiff(M, r, z)

function gradB(M::AbstractEquilibrium, x, y, z)
    r = hypot(x,y)
    phi = atan(y,x)
    gB = gradB(M, r, z)
    sp, cp = sincos(phi)
    gB_xyz = SVector{3}(gB[1]*cp - gB[2]*sp, gB[1]*sp + gB[2]*cp, gB[3])
end

function curlB_autodiff(M::AbstractEquilibrium, r, z)
    cc = cocos(M)
    B = Bfield(M,r,z)

    # Complete version
    #J = ForwardDiff.jacobian(x->Bfield(M,x[1],x[3]),SVector{3}(r,zero(r),z))
    #return cc.sigma_RpZ*SVector{3}(cylindrical_cocos(cc, J[3,2]/r - J[2,3], J[1,3] - J[3,1], (B[2]/r + J[2,1]) - J[1,2]/r))

    # Exploit axisymmetry to simplify
    J = ForwardDiff.jacobian(x->Bfield(M,x[1],x[2]),SVector{2}(r,z))
    return cc.sigma_RpZ*SVector{3}(cylindrical_cocos(cc, -J[2,2], J[1,2] - J[3,1], (B[2]/r + J[2,1])))
end

function curlB_autodiff(M::AbstractEquilibrium, x, y, z)
    J = ForwardDiff.jacobian(x->Bfield(M,x[1],x[2],x[3]),SVector{3}(x,y,z))
    return SVector{3}(J[3,2] - J[2,3], J[1,3]-J[3,1], J[2,1] - J[1,2])
end

curlB(M::AbstractEquilibrium, r, z) = curlB_autodiff(M, r, z)

function curlB(M::AbstractEquilibrium, x, y, z)
    r = hypot(x,y)
    phi = atan(y,x)
    cB = curlB(M,r,z)
    sp, cp = sincos(phi)
    cB_xyz = SVector{3}(cB[1]*cp - cB[2]*sp, cB[1]*sp + cB[2]*cp, cB[3])
    return cB_xyz
end


text:
using Contractions
using Documenter

DocMeta.setdocmeta!(Contractions, :DocTestSetup, :(using Contractions); recursive=true)

makedocs(;
    modules=[Contractions],
    authors="JCharis <jcharistech.com> and contributors",
    repo="https://github.com/jcharis/Contractions.jl/blob/{commit}{path}#{line}",
    sitename="Contractions.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://jcharis.github.io/Contractions.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/jcharis/Contractions.jl",
)