text:
abstract type TuningStrategy <: MLJ.MLJType end
const ParameterName=Union{Symbol,Expr}

"""
    Grid(resolution=10, parallel=true)

Define a grid-based hyperparameter tuning strategy, using the
specified `resolution` for numeric hyperparameters. For use with a
`TunedModel` object.

Individual hyperparameter resolutions can also be specified, as in

    Grid(resolution=[:n => r1, :(atom.max_depth) => r2])

where `r1` and `r2` are `NumericRange` objects.

See also [TunedModel](@ref), [range](@ref).

"""
mutable struct Grid <: TuningStrategy
    resolution::Union{Int,Vector{<:Pair{<:ParameterName,Int}}}
    acceleration::AbstractResource
end

# Constructor with keywords
Grid(; resolution=10, acceleration::AbstractResource=DEFAULT_RESOURCE[]) =
    Grid(resolution, acceleration)

MLJBase.show_as_constructed(::Type{<:Grid}) = true

"""
$TYPEDEF

Container for a deterministic tuning strategy.
"""
mutable struct DeterministicTunedModel{T,M<:Deterministic} <: MLJ.Deterministic
    model::M
    tuning::T  # tuning strategy
    resampling # resampling strategy
    measure
    weights::Union{Nothing,Vector{<:Real}}
    operation
    ranges::Union{Vector{<:ParamRange},ParamRange}
    full_report::Bool
    train_best::Bool
end

"""
$TYPEDEF

Container for a probabilistic tuning strategy.
"""
mutable struct ProbabilisticTunedModel{T,M<:Probabilistic} <: MLJ.Probabilistic
    model::M
    tuning::T  # tuning strategy
    resampling # resampling strategy
    measure
    weights::Union{Nothing,AbstractVector{<:Real}}
    operation
    ranges::Union{Vector{<:ParamRange},ParamRange}
    full_report::Bool
    train_best::Bool
end

const EitherTunedModel{T,M} =
    Union{DeterministicTunedModel{T,M},ProbabilisticTunedModel{T,M}}

MLJBase.is_wrapper(::Type{<:EitherTunedModel}) = true

"""
    tuned_model = TunedModel(; model=nothing,
                             tuning=Grid(),
                             resampling=Holdout(),
                             measure=nothing,
                             weights=nothing,
                             operation=predict,
                             ranges=ParamRange[],
                             full_report=true)

Construct a model wrapper for hyperparameter optimization of a
supervised learner.

Calling `fit!(mach)` on a machine `mach=machine(tuned_model, X, y)` or
`mach=machine(tuned_model, task)` will:

- Instigate a search, over clones of `model`, with the hyperparameter
  mutations specified by `ranges`, for a model optimizing the specified
  `measure`, using performance evaluations carried out using the specified
  `tuning` strategy and `resampling` strategy.

- Fit an internal machine, based on the optimal model `fitted_params(mach).best_model`,
  wrapping the optimal `model` object in *all* the provided data `X, y`
  (or in `task`). Calling `predict(mach, Xnew)` then returns predictions
  on `Xnew` of this internal machine.

*Important.* If a custom measure `measure` is used, and the measure is
a score, rather than a loss, be sure to check that
`MLJ.orientation(measure) == :score` to ensure maximization of the
measure, rather than minimization. Overide an incorrect value with
`MLJ.orientation(::typeof(measure)) = :score`.

If `measure` supports sample weights (`MLJ.supports_weights(measure)
== true`) then these can be passed to the measure as `weights`.

In the case of two-parameter tuning, a Plots.jl plot of performance
estimates is returned by `plot(mach)` or `heatmap(mach)`.

"""
function TunedModel(;model=nothing,
                    tuning=Grid(),
                    resampling=Holdout(),
                    measure=nothing,
                    weights=nothing,
                    operation=predict,
                    ranges=ParamRange[],
                    minimize=true,
                    full_report=true,
                    train_best=true)

    !isempty(ranges) || error("You need to specify ranges=... ")
    model !== nothing || error("You need to specify model=... ")
    model isa Supervised || error("model must be a SupervisedModel. ")

    message = clean!(model)
    isempty(message) || @info message

    if model isa Deterministic
        return DeterministicTunedModel(model, tuning, resampling,
           measure, weights, operation, ranges, full_report, train_best)
    elseif model isa Probabilistic
        return ProbabilisticTunedModel(model, tuning, resampling,
           measure, weights, operation, ranges, full_report, train_best)
    end
    error("$model does not appear to be a Supervised model.")
end

function MLJBase.clean!(model::EitherTunedModel)
    message = ""
    if model.measure === nothing
        model.measure = default_measure(model)
        message *= "No measure specified. Using measure=$(model.measure). "
    end
    return message
end

## GRID SEARCH

function MLJBase.fit(tuned_model::EitherTunedModel{Grid,M}, verbosity::Int, X, y) where M

    if tuned_model.ranges isa AbstractVector
        ranges = tuned_model.ranges
    else
        ranges = [tuned_model.ranges,]
    end

    ranges isa AbstractVector{<:ParamRange} ||
        error("ranges must be a ParamRange object or a vector of " *
              "ParamRange objects. ")

    # Build a vector of resolutions, one element per range. In case of
    # OrdinalRange provide a dummy value of 5. In case of a dictionary
    # with missing keys for the NumericRange`s, use fallback of 5.
    resolution = tuned_model.tuning.resolution
    if resolution isa Vector
        val_given_field = Dict(resolution...)
        fields = keys(val_given_field)
        resolutions = map(ranges) do range
            if range.field in fields
                return val_given_field[range.field]
            else
                if range isa MLJ.NumericRange && verbosity > 0
                    @warn "No resolution specified for "*
                    "$(range.field). Will use a value of 5. "
                end
                return 5
            end
        end
    else
        resolutions = fill(resolution, length(ranges))
    end

    if tuned_model.measure isa AbstractVector
        measure = tuned_model.measure[1]
        verbosity >=0 &&
            @warn "Provided `meausure` is a vector. Using first element only. "
    else
        measure = tuned_model.measure
    end

    minimize = ifelse(orientation(measure) == :loss, true, false)

    if verbosity > 0 && tuned_model.train_best
        if minimize
            @info "Mimimizing $measure. "
        else
            @info "Maximizing $measure. "
        end
    end

    parameter_names = [string(r.field) for r in ranges]
    scales = [scale(r) for r in ranges]

    # the mutating model:
    clone = deepcopy(tuned_model.model)

    resampler = Resampler(model=clone,
                          resampling=tuned_model.resampling,
                          measure=measure,
                          weights=tuned_model.weights,
                          operation=tuned_model.operation)

    resampling_machine = machine(resampler, X, y)

    # tuple of iterators over hyper-parameter values:
    iterators = map(eachindex(ranges)) do j
        range = ranges[j]
        if range isa MLJ.NominalRange
            MLJ.iterator(range)
        elseif range isa MLJ.NumericRange
            MLJ.iterator(range, resolutions[j])
        else
            throw(TypeError(:iterator, "", MLJ.ParamRange, range))
        end
    end

#    nested_iterators = copy(tuned_model.ranges, iterators)

    n_iterators = length(iterators) # same as number of ranges
    A = MLJ.unwind(iterators...)
    N = size(A, 1)

    if tuned_model.full_report
        measurements = Vector{Float64}(undef, N)
    end

    # initialize search for best model:
    best_model = deepcopy(tuned_model.model)
    best_measurement = ifelse(minimize, Inf, -Inf)
    s = ifelse(minimize, 1, -1)

    # evaluate all the models using specified resampling:
    # TODO: parallelize!

    meter = Progress(N+1, dt=0, desc="Iterating over a $N-point grid: ",
                     barglyphs=BarGlyphs("[=> ]"), barlen=25, color=:yellow)
    verbosity != 1 || next!(meter)

    for i in 1:N
        verbosity != 1 || next!(meter)

        A_row = Tuple(A[i,:])

 #       new_params = copy(nested_iterators, A_row)

        # mutate `clone` (the model to which `resampler` points):
        for k in 1:n_iterators
            field = ranges[k].field
            recursive_setproperty!(clone, field, A_row[k])
        end

        if verbosity == 2
            fit!(resampling_machine, verbosity=0)
        else
            fit!(resampling_machine, verbosity=verbosity-1)
        end
        e = evaluate(resampling_machine).measurement[1]

        if verbosity > 1
            text = prod("$(parameter_names[j])=$(A_row[j]) \t" for j in 1:length(A_row))
            text *= "measurement=$e"
            println(text)
        end

        if s*(best_measurement - e) > 0
            best_model = deepcopy(clone)
            best_measurement = e
        end

        if tuned_model.full_report
 #           models[i] = deepcopy(clone)
            measurements[i] = e
        end

    end

    if tuned_model.train_best
        verbosity < 1 || @info "Training best model on all supplied data."

        # train best model on all the data:
        # TODO: maybe avoid using machines here and use model fit/predict?
        fitresult = machine(best_model, X, y)
        fit!(fitresult, verbosity=verbosity-1)
    else
        fitresult = tuned_model.model
    end

    if tuned_model.full_report
        report = (# models=models,
                  # best_model=best_model,
                  parameter_names= permutedims(parameter_names), # row vector
                  parameter_scales=permutedims(scales),  # row vector
                  parameter_values=A,
                  measurements=measurements,
                  best_measurement=best_measurement)
    else
        report = (# models=[deepcopy(clone),][1:0],         # empty vector
                  # best_model=best_model,
                  parameter_names= permutedims(parameter_names), # row vector
                  parameter_scales=permutedims(scales),   # row vector
                  parameter_values=A[1:0,1:0],            # empty matrix
                  measurements=[best_measurement, ][1:0], # empty vector
                  best_measurement=best_measurement)
    end

    cache = nothing

    return fitresult, cache, report

end

MLJBase.fitted_params(::EitherTunedModel, fitresult) = (best_model=fitresult.model,)

MLJBase.predict(tuned_model::EitherTunedModel, fitresult, Xnew) = predict(fitresult, Xnew)
MLJBase.best(model::EitherTunedModel, fitresult) = fitresult.model

## METADATA

MLJBase.load_path(::Type{<:DeterministicTunedModel}) = "MLJ.DeterministicTunedModel"
MLJBase.package_name(::Type{<:DeterministicTunedModel}) = "MLJ"
MLJBase.package_uuid(::Type{<:DeterministicTunedModel}) = ""
MLJBase.package_url(::Type{<:DeterministicTunedModel}) = "https://github.com/alan-turing-institute/MLJ.jl"
MLJBase.is_pure_julia(::Type{<:DeterministicTunedModel{T,M}}) where {T,M} = MLJBase.is_pure_julia(M)
MLJBase.input_scitype(::Type{<:DeterministicTunedModel{T,M}}) where {T,M} = MLJBase.input_scitype(M)
MLJBase.target_scitype(::Type{<:DeterministicTunedModel{T,M}}) where {T,M} = MLJBase.target_scitype(M)

MLJBase.load_path(::Type{<:ProbabilisticTunedModel}) = "MLJ.ProbabilisticTunedModel"
MLJBase.package_name(::Type{<:ProbabilisticTunedModel}) = "MLJ"
MLJBase.package_uuid(::Type{<:ProbabilisticTunedModel}) = ""
MLJBase.package_url(::Type{<:ProbabilisticTunedModel}) = "https://github.com/alan-turing-institute/MLJ.jl"
MLJBase.is_pure_julia(::Type{<:ProbabilisticTunedModel{T,M}}) where {T,M} = MLJBase.is_pure_julia(M)
MLJBase.input_scitype(::Type{<:ProbabilisticTunedModel{T,M}}) where {T,M} = MLJBase.input_scitype(M)
MLJBase.target_scitype(::Type{<:ProbabilisticTunedModel{T,M}}) where {T,M} = MLJBase.target_scitype(M)

## LEARNING CURVES

"""
    curve = learning_curve!(mach; resolution=30, resampling=Holdout(), measure=rms, operation=predict, range=nothing, n=1)

Given a supervised machine `mach`, returns a named tuple of objects
needed to generate a plot of performance measurements, as a function
of the single hyperparameter specified in `range`. The tuple `curve`
has the following keys: `:parameter_name`, `:parameter_scale`,
`:parameter_values`, `:measurements`.

For `n` not equal to 1, multiple curves are computed, and the value of
`curve.measurements` is an array, one column for each run. This is
useful in the case of models with indeterminate fit-results, such as a
random forest.

````julia
using CSV
X, y = datanow()
atom = RidgeRegressor()
ensemble = EnsembleModel(atom=atom)
mach = machine(ensemble, X, y)
r_lambda = range(ensemble, :(atom.lambda), lower=0.1, upper=100, scale=:log10)
curve = MLJ.learning_curve!(mach; range=r_lambda)
using Plots
plot(curve.parameter_values, curve.measurements, xlab=curve.parameter_name, xscale=curve.parameter_scale)
````

If the specified hyperparameter is the number of iterations in some
iterative model (and that model has an appropriately overloaded
`MLJBase.update` method) then training is not restarted from scratch
for each increment of the parameter, ie the model is trained
progressively.

````julia
atom.lambda=1.0
r_n = range(ensemble, :n, lower=2, upper=150)
curves = MLJ.learning_curve!(mach; range=r_n, verbosity=3, n=5)
plot(curves.parameter_values, curves.measurements, xlab=curves.parameter_name)
````

"""
function learning_curve!(mach::Machine{<:Supervised};
                         resolution=30, resampling=Holdout(),
                         measure=rms, operation=predict,
                         range=nothing, verbosity=1, n=1)

    range !== nothing || error("No param range specified. Use range=... ")

    tuned_model = TunedModel(model=mach.model, ranges=range,
                             tuning=Grid(resolution=resolution),
                             resampling=resampling, measure=measure,
                             full_report=true, train_best=false)
    tuned = machine(tuned_model, mach.args...)

    measurements = reduce(hcat, [(fit!(tuned, verbosity=verbosity, force=true);
                                  tuned.report.measurements) for c in 1:n])
    report = tuned.report
    parameter_name=report.parameter_names[1]
    parameter_scale=report.parameter_scales[1]
    parameter_values=[report.parameter_values[:, 1]...]
    measurements_ = (n == 1) ? [measurements...] : measurements

    return (parameter_name=parameter_name,
            parameter_scale=parameter_scale,
            parameter_values=parameter_values,
            measurements = measurements_)
end


text:
"""
Generic implementation of an unconstrained single-field FE space
Private fields and type parameters
"""
struct UnconstrainedFESpace <: SingleFieldFESpace
  nfree::Int
  ndirichlet::Int
  cell_dofs::AbstractArray
  cell_basis::CellField
  cell_dof_basis::CellDofBasis
  dirichlet_dof_tag::Vector{Int8}
  dirichlet_cells::Vector{Int}
  ntags::Int
end

# FESpace interface

constraint_style(::Type{<:UnconstrainedFESpace}) = Val{false}()

function num_free_dofs(f::UnconstrainedFESpace)
  f.nfree
end

function get_cell_basis(f::UnconstrainedFESpace)
  f.cell_basis
end

function zero_free_values(f::UnconstrainedFESpace)
  zeros(num_free_dofs(f))
end

# SingleFieldFESpace interface

function get_cell_dofs(f::UnconstrainedFESpace)
  f.cell_dofs
end

function get_cell_dof_basis(f::UnconstrainedFESpace)
  f.cell_dof_basis
end

function num_dirichlet_dofs(f::UnconstrainedFESpace)
  f.ndirichlet
end

function num_dirichlet_tags(f::UnconstrainedFESpace)
  f.ntags
end

function zero_dirichlet_values(f::UnconstrainedFESpace)
  zeros(num_dirichlet_dofs(f))
end

function get_dirichlet_dof_tag(f::UnconstrainedFESpace)
  f.dirichlet_dof_tag
end

function scatter_free_and_dirichlet_values(f::UnconstrainedFESpace,free_values,dirichlet_values)
  cell_dofs = get_cell_dofs(f)
  LocalToGlobalPosNegArray(cell_dofs,free_values,dirichlet_values)
end

function gather_free_and_dirichlet_values!(free_vals,dirichlet_vals,f::UnconstrainedFESpace,cell_vals)

  cell_dofs = get_cell_dofs(f)
  cache_vals = array_cache(cell_vals)
  cache_dofs = array_cache(cell_dofs)
  cells = 1:length(cell_vals)

  _free_and_dirichlet_values_fill!(
    free_vals,
    dirichlet_vals,
    cache_vals,
    cache_dofs,
    cell_vals,
    cell_dofs,
    cells)

  (free_vals,dirichlet_vals)
end

function gather_dirichlet_values!(dirichlet_vals,f::UnconstrainedFESpace,cell_vals)

  cell_dofs = get_cell_dofs(f)
  cache_vals = array_cache(cell_vals)
  cache_dofs = array_cache(cell_dofs)
  free_vals = zero_free_values(f)
  cells = f.dirichlet_cells

  _free_and_dirichlet_values_fill!(
    free_vals,
    dirichlet_vals,
    cache_vals,
    cache_dofs,
    cell_vals,
    cell_dofs,
    cells)

  dirichlet_vals
end

function  _free_and_dirichlet_values_fill!(
  free_vals,
  dirichlet_vals,
  cache_vals,
  cache_dofs,
  cell_vals,
  cell_dofs,
  cells)

  for cell in cells
    vals = getindex!(cache_vals,cell_vals,cell)
    dofs = getindex!(cache_dofs,cell_dofs,cell)
    for (i,dof) in enumerate(dofs)
      val = vals[i]
      if dof > 0
        free_vals[dof] = val
      elseif dof < 0
        dirichlet_vals[-dof] = val
      else
        @unreachable "dof ids either positibe or negative, not zero"
      end
    end
  end

end

ext: jl


text:
using Distributions, Random

function simulateSDT(;d=2.,c=0.,Nd,kwargs...)
    θhit=cdf(Normal(0,1),d/2-c)
    θfa=cdf(Normal(0,1),-d/2-c)
    hits = rand(Binomial(Nd,θhit))
    fas = rand(Binomial(Nd,θfa))
      return (hits=hits,fas=fas,Nd=Nd)
 end


text:
using HDF5, Interpolations

fname = joinpath(dirname(pathof(F16Model)), "data/F16AeroData.h5");
fid = h5open(fname, "r");

function createAeroFunction(data,indepVars...)
    len = length.(indepVars);
    y = reshape(data,len);
    itp = interpolate(indepVars,y,Gridded(Linear()));
    return itp;
end

# Read independent variables
alpha1 = read(fid["alpha1"])[:];
alpha2 = read(fid["alpha2"])[:];
beta1 = read(fid["beta1"])[:];
dh1 = read(fid["dh1"])[:];
dh2 = read(fid["dh2"])[:];

# Create aero functions
_Cx = createAeroFunction(read(fid["_Cx"])[:],alpha1,beta1,dh1);
_Cy = createAeroFunction(read(fid["_Cy"])[:],alpha1,beta1);
_Cz = createAeroFunction(read(fid["_Cz"])[:],alpha1,beta1,dh1);

_Cl = createAeroFunction(read(fid["_Cl"])[:],alpha1,beta1,dh2);
_Cm = createAeroFunction(read(fid["_Cm"])[:],alpha1,beta1,dh1);
_Cn = createAeroFunction(read(fid["_Cn"])[:],alpha1,beta1,dh2);

# Leading Edge Influence
# ======================
_Cx_lef = createAeroFunction(read(fid["_Cx_lef"])[:],alpha2,beta1);
_Cy_lef = createAeroFunction(read(fid["_Cy_lef"])[:],alpha2,beta1);
_Cz_lef = createAeroFunction(read(fid["_Cz_lef"])[:],alpha2,beta1);

_Cl_lef = createAeroFunction(read(fid["_Cl_lef"])[:],alpha2,beta1);
_Cm_lef = createAeroFunction(read(fid["_Cm_lef"])[:],alpha2,beta1);
_Cn_lef = createAeroFunction(read(fid["_Cn_lef"])[:],alpha2,beta1);

# Stability Derivatives
# =====================
_Cxq = createAeroFunction(read(fid["_Cxq"])[:],alpha1);
_Cyp = createAeroFunction(read(fid["_Cyp"])[:],alpha1);
_Czq = createAeroFunction(read(fid["_Czq"])[:],alpha1);
_Cmq = createAeroFunction(read(fid["_Cmq"])[:],alpha1);

_Cyr = createAeroFunction(read(fid["_Cyr"])[:],alpha1);
_Cnr = createAeroFunction(read(fid["_Cnr"])[:],alpha1);

_Cnp = createAeroFunction(read(fid["_Cnp"])[:],alpha1);
_Clp = createAeroFunction(read(fid["_Clp"])[:],alpha1);
_Clr = createAeroFunction(read(fid["_Clr"])[:],alpha1);

_deltaCxq_lef = createAeroFunction(read(fid["_deltaCxq_lef"])[:],alpha2);
_deltaCyr_lef = createAeroFunction(read(fid["_deltaCyr_lef"])[:],alpha2);
_deltaCyp_lef = createAeroFunction(read(fid["_deltaCyp_lef"])[:],alpha2);

_deltaCzq_lef = createAeroFunction(read(fid["_deltaCzq_lef"])[:],alpha2);
_deltaClr_lef = createAeroFunction(read(fid["_deltaClr_lef"])[:],alpha2);
_deltaClp_lef = createAeroFunction(read(fid["_deltaClp_lef"])[:],alpha2);

_deltaCmq_lef = createAeroFunction(read(fid["_deltaCmq_lef"])[:],alpha2);
_deltaCnr_lef = createAeroFunction(read(fid["_deltaCnr_lef"])[:],alpha2);
_deltaCnp_lef = createAeroFunction(read(fid["_deltaCnp_lef"])[:],alpha2);

_Cy_r30 = createAeroFunction(read(fid["_Cy_r30"])[:],alpha1,beta1);
_Cn_r30 = createAeroFunction(read(fid["_Cn_r30"])[:],alpha1,beta1);
_Cl_r30 = createAeroFunction(read(fid["_Cl_r30"])[:],alpha1,beta1);

_Cy_a20 = createAeroFunction(read(fid["_Cy_a20"])[:],alpha1,beta1);
_Cy_a20_lef = createAeroFunction(read(fid["_Cy_a20_lef"])[:],alpha2,beta1);

_Cn_a20 = createAeroFunction(read(fid["_Cn_a20"])[:],alpha1,beta1);
_Cn_a20_lef = createAeroFunction(read(fid["_Cn_a20_lef"])[:],alpha2,beta1);

_Cl_a20 = createAeroFunction(read(fid["_Cl_a20"])[:],alpha1,beta1);
_Cl_a20_lef = createAeroFunction(read(fid["_Cl_a20_lef"])[:],alpha2,beta1);

_deltaCnbeta = createAeroFunction(read(fid["_deltaCnbeta"])[:],alpha1);
_deltaClbeta = createAeroFunction(read(fid["_deltaClbeta"])[:],alpha1);
_deltaCm = createAeroFunction(read(fid["_deltaCm"])[:],alpha1);

_eta_el = createAeroFunction(read(fid["_eta_el"])[:],dh1);

close(fid);

# ------------- Other derived aero functions ----------------

function Delta_lef(alpha::Real,beta::Real)::Tuple{Real,Real,Real,Real,Real,Real}
    delta_Cx_lef = _Cx_lef(alpha,beta) - _Cx(alpha,beta,0.0);
    delta_Cy_lef = _Cy_lef(alpha,beta) - _Cy(alpha,beta);
    delta_Cz_lef = _Cz_lef(alpha,beta) - _Cz(alpha,beta,0.0);

    delta_Cl_lef = _Cl_lef(alpha,beta) - _Cl(alpha,beta,0.0);
    delta_Cm_lef = _Cm_lef(alpha,beta) - _Cm(alpha,beta,0.0);
    delta_Cn_lef = _Cn_lef(alpha,beta) - _Cn(alpha,beta,0.0);

    return (delta_Cx_lef, delta_Cy_lef, delta_Cz_lef,delta_Cm_lef, delta_Cn_lef, delta_Cl_lef);
end

# Damping Terms
function _Damping(alpha::Real)::Tuple{Real,Real,Real,Real,Real,Real,Real,Real,Real}
    Cxq = _Cxq(alpha);
    Cyr = _Cyr(alpha);
    Cyp = _Cyp(alpha);
    Czq = _Czq(alpha);
    Clr = _Clr(alpha);
    Clp = _Clp(alpha);
    Cmq = _Cmq(alpha);
    Cnr = _Cnr(alpha);
    Cnp = _Cnp(alpha);

    return (Cxq, Cyr, Cyp, Czq, Clr, Clp, Cmq, Cnr, Cnp);
end

function _Damping_lef(alpha::Real)::Tuple{Real,Real,Real,Real,Real,Real,Real,Real,Real}
    delta_Cxq_lef = _deltaCxq_lef(alpha);
    delta_Cyr_lef = _deltaCyr_lef(alpha);
    delta_Cyp_lef = _deltaCyp_lef(alpha);
    delta_Czq_lef = _deltaCzq_lef(alpha);
    delta_Clr_lef = _deltaClr_lef(alpha);
    delta_Clp_lef = _deltaClp_lef(alpha);
    delta_Cmq_lef = _deltaCmq_lef(alpha);
    delta_Cnr_lef = _deltaCnr_lef(alpha);
    delta_Cnp_lef = _deltaCnp_lef(alpha);

    return (delta_Cxq_lef, delta_Cyr_lef, delta_Cyp_lef, delta_Czq_lef, delta_Clr_lef, delta_Clp_lef, delta_Cmq_lef, delta_Cnr_lef, delta_Cnp_lef);
end

# Rudder Influence
function _RudderInfluence(alpha::Real,beta::Real)::Tuple{Real,Real,Real}
    delta_Cy_r30 = _Cy_r30(alpha,beta) - _Cy(alpha,beta);
    delta_Cn_r30 = _Cn_r30(alpha,beta) - _Cn(alpha,beta,0.0);
    delta_Cl_r30 = _Cl_r30(alpha,beta) - _Cl(alpha,beta,0.0);

    return (delta_Cy_r30, delta_Cn_r30, delta_Cl_r30);
end

 function _AileronInfluence(alpha::Real,beta::Real)::Tuple{Real,Real,Real,Real,Real,Real}
    delta_Cy_a20 = _Cy_a20(alpha,beta) - _Cy(alpha,beta);
    delta_Cy_a20_lef = _Cy_a20_lef(alpha,beta) - _Cy_lef(alpha,beta) - delta_Cy_a20;

    delta_Cn_a20 = _Cn_a20(alpha,beta) - _Cn(alpha,beta,0.0);
    delta_Cn_a20_lef = _Cn_a20_lef(alpha,beta) - _Cn_lef(alpha,beta) - delta_Cn_a20;

    delta_Cl_a20 = _Cl_a20(alpha,beta) - _Cl(alpha,beta,0.0);
    delta_Cl_a20_lef = _Cl_a20_lef(alpha,beta) - _Cl_lef(alpha,beta) - delta_Cl_a20;

    return (delta_Cy_a20, delta_Cy_a20_lef, delta_Cn_a20, delta_Cn_a20_lef, delta_Cl_a20, delta_Cl_a20_lef);
 end

 function _OtherCoefficients(alpha,el)::Tuple{Real,Real,Real,Real}
    delta_Cnbeta = _deltaCnbeta(alpha);
    delta_Clbeta = _deltaClbeta(alpha);
    delta_Cm = _deltaCm(alpha);
    eta_el = _eta_el(el);

    return (delta_Cnbeta, delta_Clbeta, delta_Cm, eta_el);
 end


text:
__precompile__()
# Set of helper routines for WOS fractional Laplacian
#
#
module KOSHELPER
#
using PyPlot
using Interpolations
using Dierckx
#
export init, MyTriMesh,get_tri_mesh,Final_Arnoldi, Add_Arnoldi!,Construct_Arnoldi,
Arnoldi,Matern,mat_cov,prior,new_data!,L2norm,L2L2norm,TheData,MyNormal,norm_gradient, dyda_compare, do_plot, Compare, smooth_interpolate!, smooth_interpolate2!
function __init__()
end
###################### begin: TriMesh
type MyTriMesh
  x::Array{Float64,1}
  y::Array{Float64,1}
  mg_x::Array{Float64,2}
  mg_y::Array{Float64,2}
  hx::Float64
  hy::Float64
  nx::Int64
  ny::Int64
  n_int::Int64
  n_ext::Int64
  get_pt::Function
  get_pt_ext::Function
  get_pt_all::Function
  get_ext::Function
  reshape::Function
  evalfn::Function # evaluate (f,g) on TriMesh, f interior, g exterior.
  int::Array{Int64,1}
  ext::Array{Int64,1}
  get_levels::Function
  parents::Array{Int64,2}
  get_parents::Function
  finest_parents::Array{Int64,1}
  get_level_ind::Function
  move_to_TriMesh::Function
  sq_int::Function
  sq_int2::Function
  scale::Function
end
###
function form_triangular_TriMesh2(pts_in,line0,tri0,no_levels)
  # pts_in, initial pts, e.g., 5 corners for square (inc centre)
  # tri0, initial triang as [n,3] matrix of point indices
  # assume triangles are well-ordered so [1,2,3]
  # means lines 1,2,3 are order [a,b][b,c][c,a],
  # where a,b,c are point numbers
  #
  level=zeros(Int64,no_levels);
  level[1]=size(pts_in)[1];
  #
  ndim=8*4^(no_levels)
  pts=zeros(ndim,2);   pts[1:level[1],:]=pts_in
  parents=zeros(Int64,ndim,2)
  #
  cur_tri_no=size(tri0)[1]
  cur_lin_no=size(line0)[1]
  #
  for l=1:(no_levels-2)
    k=0;    line_no=0
    line=zeros(Int64,4*cur_lin_no,2)
    tri =zeros(Int64,4*cur_tri_no,3)
    tri_no=0
    midpoint=zeros(Int64,cur_lin_no)
    halflines=zeros(Int64,4*cur_lin_no,2)
    # divide lines to get new lines and points
    for i=1:cur_lin_no # split lines
      pt_no=level[l]+i
      n1=line0[i,1]; n2=line0[i,2] # point numbers
      pts[pt_no,:]=(pts[n1,:]+pts[n2,:])/2
      # add 2 lines
      line[line_no+1,:]=[n1,   pt_no]; #  order is important
      line[line_no+2,:]=[pt_no,n2] # preserve clockwise ordering
      line_no+=2
      # record parents
      parents[pt_no,:]=[n1,n2]
      # need map from lines to midpts
      midpoint[i]=pt_no # line0 num to pt_num
      halflines[i,:]=[line_no-1, line_no] # map from cur_lin_no to two half line_no
    end
    new_pts=(1+level[l]):(level[l]+cur_lin_no)
    # loop over current triangles
    for k=1:cur_tri_no#
      tri_line_no=tri0[k,:] # line0 numbers for this tri
      hl=halflines[tri_line_no,:]
      # sprintln("\nstart", midpoint[tri_line_no]," ",hl, " ", line[hl[:],:] )
      # for each successive pair of lines
      for i=1:3
        # common point from
        common_ind=line0[tri_line_no[i],2] # need triangle is well-ordered
        j=mod(i, 3)+1
        assert(line0[tri_line_no[i],2]==line0[tri_line_no[j],1])
        # midpoints
        mid1=midpoint[tri_line_no[i]] # =pt_no for line i
        mid2=midpoint[tri_line_no[j]] # =pt_no for line j
        # add line
        line[line_no+1,:]=[mid2,mid1]
        line_no+=1
        # add trinagle
        # [mid2,mid1] [pt_no i,n2 i] [ n1 j, pt_no j]
        assert(line[line_no,2]==line[hl[i,2],1])
        assert(line[line_no,1]==line[hl[j,1],2])
        tri[tri_no+1,:]=[line_no, hl[i,2], hl[j,1]]
        tri_no+=1
      end
      # add the interior triangle
      tri[tri_no+1,:]=[line_no-2,line_no,line_no-1]
      tri_no+=1
    end
    # store number of lines
    level[l+1]=level[l]+cur_lin_no
    cur_lin_no=line_no
    line0=line
    tri0=tri
    cur_tri_no=tri_no
  end
  return pts,level,parents,line0,cur_lin_no
end
###
function get_tri_mesh(xy,no_levels::Int64,d)
  print("Creating triangular mesh...")
  # for plotting TriMesh
  nx=1+2^(no_levels-1); ny=nx
  #
  xTriMesh=linspace(xy[1],xy[2],nx)
  yTriMesh=linspace(xy[3],xy[4],ny)
  #
  hx=(xy[2]-xy[1])/(nx-1)
  hy=(xy[4]-xy[3])/(ny-1)
  #
  mg_x = repmat(xTriMesh[:]',ny,1) # used for plotting
  mg_y = repmat(yTriMesh[:] ,1,nx)
  #
  pts0=zeros(5,2)
  pts0[1,:]=[xy[1],xy[3]]; pts0[2,:]=[xy[1],xy[4]];
  pts0[3,:]=[xy[2],xy[3]]; pts0[4,:]=[xy[2],xy[4]];
  pts0[5,:]=(pts0[1,:]+pts0[4,:])/2
  line0=[1 2; 1 3;  4 2; 4 3;        5 1; 2 5; 3 5;5 4 ]
  tri0=[1 6 5; 2 7 5; 7 8 4; 3 6 8] # 4x3
  #
  pts,level,parents,line,cur_lin_no=form_triangular_TriMesh2(pts0,
                                      line0, tri0,
                                      no_levels)
  qpts,qlevel,qparents,=form_triangular_TriMesh2(pts0,
                                      line0, tri0,
                                      no_levels+2)
  # add points on finest level to make a rectangular TriMesh
  pt_no=level[no_levels-1]
  for i=1:cur_lin_no
    n1=line[i,1]; n2=line[i,2]
    if (pts[n1,1]==pts[n2,1] || pts[n1,2]==pts[n2,2])
      pt_no+=1
      pts[pt_no,:]=(pts[n1,:]+pts[n2,:])/2
      parents[pt_no,:]=[n1,n2]
    end
  end
  level[no_levels]=pt_no
  ndim=level[end]
  #
  ds=zeros(ndim); xTriMesh=deepcopy(ds); yTriMesh=deepcopy(ds)
  for i=1:level[no_levels]
    ds[i]=d(pts[i,:])[1]
    xTriMesh[i]=pts[i,1]; yTriMesh[i]=pts[i,2]
  end
  #
  int=find(x-> x>0, ds)[:];  n_int=length(int)
  ext=find(x-> x<=0,ds)[:];  n_ext=length(ext)
  #
  function get_pt_all(i) # return ith pt
    return [xTriMesh[i] yTriMesh[i]]
  end
  #
  function get_pt_ext(i) # return ith exterior pt
    return get_pt_all(ext[i])[:]
  end
  #
  function get_pt_int(i) # return ith interior point
    return get_pt_all(int[i])[:]
  end
  #
  function get_parents(i)
    if i>level[1] && i<=level[no_levels]
      return parents[i,:]
    else
      return [0,0]
    end
  end
  ############################
  # get_finest_parents
    r=(1+level[no_levels-1]):level[no_levels]
    pa=zeros(Int64,length(r),2)
    for i in r # update on finest interior
      j=i-level[no_levels-1]
      pa[j,:]     = parents[i,:]
    end
    fpa=union(unique(pa[:]),r)
  ##############################
  #######################
  function get_ext(out,fn)# interior vector, ext fn
    # for i=1:n_int
    #      ri=int[i]
    #      #out[ri]=inm[i]
    # end
    for i=1:n_ext
       ri=ext[i]
       out[ri]=fn(get_pt_ext(i))
    end
    return out
  end
  #######################
    function get_reshape(inm,fn)# interior vector, ext fn
    out=zeros(level[no_levels])
    for i=1:n_int
         ri=int[i]
         out[ri]=inm[i]
    end
    for i=1:n_ext
       ri=ext[i]
       out[ri]=fn(get_pt_ext(i))
    end
    return out
  end
  ##############
  function get_fn(fin,fout)
    out=zeros(level[no_levels])
    for i=1:n_int
      ri=int[i]
      out[ri]=fin(get_pt_int(i))
    end
    for i=1:n_ext
      ri=ext[i]
      out[ri]= fout(get_pt_ext(i))
    end
    return out
  end
  ############################
  function get_levels(l)
    if (l<1 || l>no_levels)
      error("level number out of range")
    end
    r=1:level[l]
    out=zeros(2,level[l])
    out[1,:]=xTriMesh[r]
    out[2,:]=yTriMesh[r]
    return out# colums are pints
  end
  #
  function only_level_ind(l)
    if (l>1)
      r=(1+level[l-1]):level[l]
    else
      r=1:level[1]
    end
    return r
  end
  # ###########################
  # function move_to_TriMesh(Z)
  #   s=Spline2D(xTriMesh,yTriMesh,Z; kx=1, ky=1, s=0.2)
  #   # on a TriMesh evaluate spline
  #   return evalTriMesh(s,mg_x[1,:],mg_y[:,1])
  # end
  ###########################
  function to_TriMesh(in)
    Z=zeros(nx,ny)
    for i=1:level[no_levels]
      #println( (pts[i,1]-xy[1])/hx ," ",(pts[i,2]-xy[3])/hy)
      j=1+Int64(round((pts[i,1]-xy[1])/hx ))
      k=1+Int64(round((pts[i,2]-xy[3])/hy ))
      #passed test!
      #println(pts[i,:], " ", mg_x[1,j], " ",mg_y[k,1])
      Z[j,k]=in[i]
    end
    return Z
  end
  ################################
  function sq_int(Z,ell)
    if ell<no_levels
      sum=quad_sq_tri(Z,ell)
    else
      H1=(xy[2]-xy[1])*(xy[4]-xy[3])/4
      H=H1*4.0^(2-ell)#area of square quarters at each level
      sum=quad_sq_TriMesh(to_TriMesh(Z),H)
    end
    return sum
  end
  ######################################
  function quad_sq_tri(Z, ell)
    # Let f be the pw linear interpolant
    # on triangles on level ell
    # compute int f^2 dx exact
    # note level ell=1 integrates the 5-point initial 4-triangle TriMesh
    # by looping over all level-2 points and working out averages of parents
    # (which are level-1 points)
    H1=(xy[2]-xy[1])*(xy[4]-xy[3])/4
    H=H1*(4.0^(1-ell))#triangle's areas quarters at each level
    sum=0.
    t = x-> (x[1]>xy[1] && x[1]<xy[2] && x[2]>xy[3] && x[2]<xy[4]) # test interior
    for i=(1+qlevel[ell]):qlevel[ell+1]# all points on level ell
      # average all parents
      pa=qparents[i,:]
      tv1=t(pts[pa[1],:]); tv2=t(pts[pa[2],:])
      if tv1 || tv2
        sum+=mean(Z[pa]).^2*2
      else
        sum+=Z[pa[1]]^2
      end
    end
    return  sum*H*(1./3)
  end
  ##################
  function quad_sq_TriMesh(Z, H)
    # Let f be the  pw linear interpolant of
    # a TriMesh function Z (defined by  pyramid function on each square
    # of the TriMesh using midpoint and average of each corner)
    # compute int f^2 \,dx exactly
    # H=size of square
    mysum=0.
    nx,ny=size(Z); v=zeros(4)
    for i=1:nx-1
      for j=1:ny-1
        x=i*hx+xy[1]; y=j+hy*xy[3];
        v[1]=Z[i,j];   v[2]=Z[i+1,j]
        v[3]=Z[i+1,j+1]; v[4]=Z[i,j+1]
        tmp=(v+circshift(v,[1]))/2
        c=sum(v.^2)+4*sum(tmp.^2)+16*mean(v)^2
        mysum+=c
      end
    end
    return mysum*H*(1/36)
  end
  ################
  function quad_sq_TriMesh2(Z_, f,ell)
    # Let f be the  pw linear interpolant of
    # a TriMesh function Z (defined by  pyramid function on each square
    # of the TriMesh using midpoint and average of each corner)
    # compute int f^2 \,dx exactly
    # H=size of square
    H1=(xy[2]-xy[1])*(xy[4]-xy[3])/4
    H=H1*4.0^(2-ell)#area of square quarters at each level
    Z=to_TriMesh(Z_)
    mysum=0.
    nx,ny=size(Z); v=zeros(4)
    for i=1:nx-1
      for j=1:ny-1
        x=(i-1)*hx+xy[1]; y=(j-1)*hy+xy[3];
        apts=[[x, y], [x+hx, y], [x+hx, y+hy], [x, y+hy]]; vf=f.(apts)
        hpts=(apts+circshift(apts,[1]))/2; hf=f.(hpts)
        mpts=mean(apts); mf=f(mpts)
        v[1]=Z[i,j];   v[2]=Z[i+1,j]
        v[3]=Z[i+1,j+1]; v[4]=Z[i,j+1]
        tmp=(v+circshift(v,[1]))/2
        c=sum((v-vf).^2)+4*sum((tmp-hf).^2)+16*(mean(v)-mf)^2
        mysum+=c
      end
    end
    return mysum*H*(1/36)
  end
  #################################
  function scale1(ell) #scaling for norm
    return sqrt((xy[2]-xy[1])*(xy[4]-xy[3]))*2.0^(-ell)/2
  end
  println("......finished.")
  ##################################
  return MyTriMesh(xTriMesh, yTriMesh, mg_x, mg_y, hx, hy, nx,ny, n_int,n_ext,
   get_pt_int,get_pt_ext,get_pt_all,get_ext,get_reshape,get_fn,int,ext,
   get_levels,parents,get_parents, fpa,only_level_ind,
    to_TriMesh, sq_int,quad_sq_TriMesh2,scale1)
end
###################################################
###################################################
###################################################
function norm_gradient(G,itp)
  g=zeros(2)
  out=zeros(G.n_int)
  for i=1:G.n_int
    pt=G.get_pt(i)
    gradient!(g,itp,pt[1],pt[2])
    out[i]=norm(g)*sqrt(G.hx*G.hy)
  end
  return out
end
################# end: TriMesh
#########################################################
#########################################################
#########################################################

# function L2L2norm(dist::MyNormal, TriMesh::MyTriMesh)
#   mysum=trace(dist.C)*TriMesh.hx*TriMesh.hy
#   return sqrt(mysum)
# end
# #
# function L2norm(dist::MyNormal, TriMesh::MyTriMesh)
#   mysum=sum(dist.mu.^2)*TriMesh.hx*TriMesh.hy
#   return sqrt(mysum)
# end
###################### begin: Arnoldi
###
type Arnoldi # method
  KSS::Set # Krylov sub-space
  H::Array # H matrix
  k::Int32 # dimension of H matrix (k+1 -dim Krylov space)
  N::Int32 # dimension of state space
  rhs::Array{Float64,1} # for next solve
  ascale::Float64
end
#
function Construct_Arnoldi(x::Array{Float64,1})
  ascale=sqrt(length(x)*1.0)
  xl=norm(x)
  tmp=x/xl
  S=Set()
  push!(S,[1,tmp])
  return Arnoldi(S,Array{Float64}(1,1),0,length(x),tmp*ascale,ascale)
end

######################
function SmartRestart!(arnoldi_data::Arnoldi,p)
  # Reduce the Arnoldi matrix using implicitly shifted Arnoldi
  # n denotes the dimension of Arnoldi H-matrix (n+1 dim Krylov basis)
  # reduced to dimension  p=n-length(mu) dim H-matrix
  #
  k=arnoldi_data.k # dimension of H-matrix
  if (k>p)
    println("Implicit restart: reducing dimension ", k," to ", p,".")
    rhs=arnoldi_data.rhs*arnoldi_data.H[k+1,k]/arnoldi_data.ascale
    # get Ritz values of H (for shifts)
    H=arnoldi_data.H[1:k,1:k]
    d=eigvals(H) #
    d=sort(d,by=abs)#ordered with increasing magnitue
    # number of shifts required is k-p
    # so throw  away p smallest evalues
    # mu is empty if p \geq k.%
    mu=d[1:(k-p)] #
    println("Ritz values=", d, " and using shifts=", mu)
    # initialise and do implicit shifts
    # load matrix V with first-k Krylov vectors
    V=zeros(arnoldi_data.N, k)
    for vs in arnoldi_data.KSS
      if (vs[1]<=k)
        V[:,vs[1]]=vs[2]
     end
    end #
    # initialise values of v
    v=zeros(k); v[end]=1
    # loop over shifts
    for i=1:length(mu)
      Q,R=qr(H-mu[i]*eye(k))
      H=Q'*H*Q # ' is already Hermitian transpose
      V=V*Q;      v=Q'*v
    end
    # define p-dimensional H matrix
    arnoldi_data.H=zeros(p+1,p+1) # yes it has dim p+1
    arnoldi_data.H[1:p,1:p]=H[1:p,1:p] # copy from above
    arnoldi_data.k=p # new dimension
    empty!(arnoldi_data.KSS) # Krylov space
    for i=1:p
      push!(arnoldi_data.KSS,[i,V[:,i]])
    end #
    # last Krylov vector and new rhs
    z=V[:,p+1]*H[p+1,p]+rhs*v[p]; norm_z=norm(z); z=z/norm_z
    push!(arnoldi_data.KSS,[p+1,z])
    arnoldi_data.H[p+1,p]=norm_z # next rhs
    arnoldi_data.rhs=z*arnoldi_data.ascale
  end
end
#######################
function Add_Arnoldi!(arnoldi_data::Arnoldi,zin::Array{Float64,1})
  # arnoldi_data (Krylov Subspace as array of [k,matrix], each matrix of unit F-norm, size of psace,
  # current H matrix)
  # initialise Arnoldi  with (v,[0],0) where v is unit rhs and data is solve of Ax=rhs
  k=arnoldi_data.k+1
  H=deepcopy(arnoldi_data.H)
  z=deepcopy(zin)/arnoldi_data.ascale
  # H should be k x k matrix already
  for vs in arnoldi_data.KSS,
    H[vs[1],k]=vecdot(z,vs[2])
    z=z-H[vs[1],k]*vs[2]
  end # end Gram-Schmidt
  d,v=eig(arnoldi_data.H[1:k,1:k])
  perm=sortperm(d,by=abs);  d=d[perm] #ordered with increasing magnitue
  ev_max=d[end]  # largest eigenvalue
  if k>1
    sep=abs(ev_max-d[end-1]) # spectral gap
  else
    sep=1.0
  end
  #
  value=0.
  # increase size of H
  arnoldi_data.H=zeros(Float64,k+1,k+1)
  arnoldi_data.H[1:k,1:k]=H
  arnoldi_data.H[k+1,k]=norm(z)
  #
  z=z/arnoldi_data.H[k+1,k]
  push!(arnoldi_data.KSS,[k+1,z])
  arnoldi_data.k=k
  #
  ev_max, evec,all_lam,sep=Final_Arnoldi(arnoldi_data)
  value=abs(arnoldi_data.H[k+1,k])*abs(evec[k])
  arnoldi_data.rhs=z*arnoldi_data.ascale
  # eigenvalue approx, new RHS
  return ev_max,  value*arnoldi_data.ascale, sep
end
#
function Final_Arnoldi(arnoldi_data::Arnoldi)
  # get eigenfunction aswell and return
  k=arnoldi_data.k
  d,v=eig(arnoldi_data.H[1:k,1:k])
  perm=sortperm(d,by=abs)#ordered with increasing magnitue
  d=d[perm]; lam=d[end]; ev=v[:,perm[end]]; #largest eigenvalue
  if k>1
    sep=abs(lam-d[end-1]) # spectral gap
  else
    sep=1.
  end
  evec=zeros(arnoldi_data.N)
  for vs in arnoldi_data.KSS
    if (vs[1]<=length(ev))
      evec=evec+ev[vs[1]]*vs[2]
    end
  end # end Gram-Schmidt
  return lam,evec,d,sep
end
# end: Arnoldi
########################
function dyda_lam(alp,d=2)
  # Dyda lower bound for smallest eigenvalue
  lam=2^alp*gamma(alp/2+1)*gamma((alp+d)/2)*(alp+2)*(alp+d)*(6-alp)/gamma(d/2)/(12*d+(16-2*d)*alp)
  return lam
end
########################
function dyda_lam_up(alp)
  data=[1.05096, 1.10993, 1.34374, 2.00612, 3.27594, 4.56719, 5.13213]
  valp=([0.1, 0.2, 0.5, 1, 1.5, 1.8, 1.9],)
  itp = interpolate(valp,data, Gridded(Linear()))
  return itp[alp]
end
##
function dyda_compare(alpha, lam)
  lam=real(lam)
  if (lam>dyda_lam(alpha) && lam<dyda_lam_up(alpha))
     println("   >>>> eigenvalue in Dyda's interval <<<<")
  else
    println("Outside Dyda's interval=(",dyda_lam(alpha),",",dyda_lam_up(alpha),").")
  end
  println("Relative error of lam=",lam," is ",
      min(abs(dyda_lam(alpha)-lam), abs(dyda_lam_up(alpha)-lam))/lam)
end

###################### begin: plotting
function do_plot(TriMesh::MyTriMesh,Z,mytitle,filename)

  PyPlot.matplotlib[:rc]("text", usetex=true) # allow tex rendering
  PyPlot.matplotlib[:rc]("xtick", labelsize=7)
  PyPlot.matplotlib[:rc]("axes", labelsize=8)
  PyPlot.matplotlib[:rc]("axes", labelpad=1)
 ticklabel_format(style="sci")
  fig = figure("pyplot_surfaceplot",figsize=(4.5,3.5))
  clf()
  ax = fig[:add_subplot](2,1,1, projection = "3d")
  ax[:plot_surface](TriMesh.mg_x, TriMesh.mg_y, Z, edgecolors="k",
   cmap=ColorMap("gray"), alpha=0.8, linewidth=0.25)
  xlabel(L"x");  ylabel(L"y"); title(mytitle)
  fig[:canvas][:draw]
  tmp=gcf()
  locator_params(nbins=5)
  ticklabel_format(style="sci",axis="z",scilimits=(0,0))
  tmp[:set_size_inches](4.7,3.5)
  savefig(filename)#,bbox_inches="tight")
end
#
function do_plot(TriMesh::MyTriMesh,Z1,Z2,mytitle1,mytitle2)
  fig = figure("pyplot_surfaceplot",figsize=(10,10))
  #
  ax = fig[:add_subplot](2,1,1, projection = "3d")
  ax[:plot_surface](TriMesh.mg_x, TriMesh.mg_y, Z1, edgecolors="k",
   cmap=ColorMap("gray"), alpha=0.8, linewidth=0.25)
  xlabel("X")
  ylabel("Y")
  title(mytitle1)
  #
  ax = fig[:add_subplot](2,1,2, projection = "3d")
  ax[:plot_surface](TriMesh.mg_x, TriMesh.mg_y, Z2, edgecolors="k",
   cmap=ColorMap("gray"), alpha=0.8, linewidth=0.25)
  xlabel("X")
  ylabel("Y")
  title(mytitle2)
  #
  tmp=fig[:canvas][:draw]
  tmp=gcf()
end
###
function Compare(G, f1, Z, N,mytitle1,mytitle2)
  gx=linspace(G.x[1],G.x[end],N)
  gy=linspace(G.y[1],G.y[end],N)
  mg_x = repmat(gx[:]',N,1)
  mg_y = repmat(gy[:],1,N)
  itp = interpolate((G.x,G.y),Z,Gridded(Linear()))
  f1v=zeros(N,N); f2v=zeros(N,N)
  for i=1:N
    for j=1:N
      if norm([gx[i],gy[j]])<1
        f1v[i,j]=f1([gx[i],gy[j]])
      end
      f2v[i,j]=itp[gx[i],gy[j]]
    end
  end
  fig = figure("pyplot_surfaceplot",figsize=(10,10))
  #
  ax = fig[:add_subplot](2,1,1, projection = "3d")
  ax[:plot_surface](mg_x, mg_y, f2v, edgecolors="k",
   cmap=ColorMap("gray"), alpha=0.8, linewidth=0.25)
  xlabel("X")
  ylabel("Y")
  title(mytitle1)
  #
  ax = fig[:add_subplot](2,1,2, projection = "3d")
  ax[:plot_surface](mg_x, mg_y, abs(f1v-f2v), edgecolors="k",
   cmap=ColorMap("gray"), alpha=0.8, linewidth=0.25)
  xlabel("X")
  ylabel("Y")
  title(mytitle2)
  #
  fig[:canvas][:draw]
  gcf()
end
# end: plotting

# ######################
function smooth_interpolate!(G,dict,p,data_in,sig_sq)
  N=G.n_int
  new_rhs_full=G.reshape(data_in,x->0.);  Z=G.move_to_TriMesh(new_rhs_full)
  if (true) # linear interpolant
    # Interpolations module
      itp=interpolate((G.mg_x[1,:],G.mg_y[:,1] ),
                      Z,
                      Gridded(Linear()))
      p.fn_rhs=x->[itp[x[i,1],x[i,2]] for i=1:size(x)[1]]#*max.(0,sign.(p.d(x)))
  else # spline interpolant
      itp0=interpolate(Z,BSpline(Cubic(Line())), OnTriMesh())
      itp = scale(itp0, G.x, G.y)
      p.fn_rhs=x->itp[x[1],x[2]]
  end
  if haskey(dict,"extra_plot")
    do_plot(G,Z,"rhs")
  end
end
######################
function smooth_interpolate2!(G,dict,p,data_in,sig_sq)
  prior=deepcopy(dict["prior"])
  N=G.n_int;   data=TheData(zeros(N),zeros(N),zeros(N))
  data.mu=data_in;
  #println(data_in) good
  #println(sig_sq)good
  data.var=ones(data_in)*sig_sq
  data.no=ones(Int32,N)
  if (false)
    new_data!(prior, data)
    #println(prior.mu)#bad
    new_rhs_full=G.reshape(prior.mu,KOS.fn_zero) # full object
  else
    for i=1:G.n_int
      data_in[i]/=p.d(G.get_pt(i))
    end
    new_rhs_full=G.reshape(data_in,x->0.)
  end
  #println(new_rhs_full)bad

  Z=G.move_to_TriMesh(new_rhs_full)
  itp=interpolate((G.mg_x[1,:],G.mg_y[:,1] ),
                      Z,
                      Gridded(Linear()))
  #p.fn_rhs=x->itp[x[1],x[2]]*max(0,p.d(x))
  p.fn_rhs=x->max(0,p.d(x))
  do_plot(G,Z,"rhs")
end
######################
#
end # module KOSHELPER


text:
"""
scipy.stats module

- [Statistical functions (scipy.stats) Reference Guide](https://docs.scipy.org/doc/scipy/reference/stats.html)

# Examples

You can calculate Pearson correlation coefficient and p-value:

```julia-repl
julia> a = [0, 0, 0, 1, 1, 1, 1];

julia> b = collect(0:6);

julia> SciPy.stats.pearsonr(a, b)
(0.8660254037844386, 0.011724811003954649)
```
"""
module stats

using PyCall
import PyCall: hasproperty # Base.hasproperty in Julia 1.2

import ..pystats
import .._generate_docstring
import ..LazyHelp

const _ignore_funcs = ["stats", "describe"]

@doc LazyHelp(pystats, "describe") function describe(args...; kws...)
   desc = pycall(pystats.describe, PyAny, args...; kws...)
   return (nobs=desc[1], minmax=desc[2], mean=desc[3], variance=desc[4], skewness=desc[5], kurtosis=desc[6])
end

# This list can be generated by print_scipy_api_list.py script
all_properties = ["F_onewayBadInputSizesWarning", "F_onewayConstantInputWarning", "PearsonRConstantInputWarning", "PearsonRNearConstantInputWarning", "SpearmanRConstantInputWarning", "alpha", "anderson", "anderson_ksamp", "anglit", "ansari", "arcsine", "argus", "bartlett", "bayes_mvs", "bernoulli", "beta", "betabinom", "betaprime", "binned_statistic", "binned_statistic_2d", "binned_statistic_dd", "binom", "binom_test", "boltzmann", "boxcox", "boxcox_llf", "boxcox_normmax", "boxcox_normplot", "bradford", "brunnermunzel", "burr", "burr12", "cauchy", "chi", "chi2", "chi2_contingency", "chisquare", "circmean", "circstd", "circvar", "combine_pvalues", "contingency", "cosine", "cramervonmises", "crystalball", "cumfreq", "describe", "dgamma", "dirichlet", "distributions", "dlaplace", "dweibull", "energy_distance", "entropy", "epps_singleton_2samp", "erlang", "expon", "exponnorm", "exponpow", "exponweib", "f", "f_oneway", "fatiguelife", "find_repeats", "fisher_exact", "fisk", "fligner", "foldcauchy", "foldnorm", "friedmanchisquare", "gamma", "gausshyper", "gaussian_kde", "genexpon", "genextreme", "gengamma", "genhalflogistic", "geninvgauss", "genlogistic", "gennorm", "genpareto", "geom", "gilbrat", "gmean", "gompertz", "gstd", "gumbel_l", "gumbel_r", "halfcauchy", "halfgennorm", "halflogistic", "halfnorm", "hmean", "hypergeom", "hypsecant", "invgamma", "invgauss", "invweibull", "invwishart", "iqr", "itemfreq", "jarque_bera", "johnsonsb", "johnsonsu", "kappa3", "kappa4", "kde", "kendalltau", "kruskal", "ks_1samp", "ks_2samp", "ksone", "kstat", "kstatvar", "kstest", "kstwo", "kstwobign", "kurtosis", "kurtosistest", "laplace", "laplace_asymmetric", "levene", "levy", "levy_l", "levy_stable", "linregress", "loggamma", "logistic", "loglaplace", "lognorm", "logser", "loguniform", "lomax", "mannwhitneyu", "matrix_normal", "maxwell", "median_abs_deviation", "median_absolute_deviation", "median_test", "mielke", "mode", "moment", "mood", "morestats", "moyal", "mstats", "mstats_basic", "mstats_extras", "multinomial", "multiscale_graphcorr", "multivariate_hypergeom", "multivariate_normal", "multivariate_t", "mvn", "mvsdist", "nakagami", "nbinom", "ncf", "nct", "ncx2", "nhypergeom", "norm", "normaltest", "norminvgauss", "obrientransform", "ortho_group", "pareto", "pearson3", "pearsonr", "percentileofscore", "planck", "pointbiserialr", "poisson", "power_divergence", "powerlaw", "powerlognorm", "powernorm", "ppcc_max", "ppcc_plot", "probplot", "randint", "random_correlation", "rankdata", "ranksums", "rayleigh", "rdist", "recipinvgauss", "reciprocal", "relfreq", "rice", "rv_continuous", "rv_discrete", "rv_histogram", "rvs_ratio_uniforms", "scoreatpercentile", "sem", "semicircular", "shapiro", "siegelslopes", "sigmaclip", "skellam", "skew", "skewnorm", "skewtest", "spearmanr", "special_ortho_group", "statlib", "stats", "t", "theilslopes", "tiecorrect", "tmax", "tmean", "tmin", "trapezoid", "trapz", "triang", "trim1", "trim_mean", "trimboth", "truncexpon", "truncnorm", "tsem", "tstd", "ttest_1samp", "ttest_ind", "ttest_ind_from_stats", "ttest_rel", "tukeylambda", "tvar", "uniform", "unitary_group", "variation", "vonmises", "vonmises_line", "wald", "wasserstein_distance", "weibull_max", "weibull_min", "weightedtau", "wilcoxon", "wishart", "wrapcauchy", "yeojohnson", "yeojohnson_llf", "yeojohnson_normmax", "yeojohnson_normplot", "yulesimon", "zipf", "zmap", "zscore"]

for f in all_properties
    f in _ignore_funcs && continue

    sf = Symbol(f)
    @eval @doc LazyHelp(pystats, $f) $sf(args...; kws...) = pycall(pystats.$f, PyAny, args...; kws...)
end

function __init__()
    copy!(pystats, pyimport_conda("scipy.stats", "scipy"))
end

end # module


text:
module SurfaceReactions
using LightXML
using DifferentialEquations
include("Constants.jl")
include("Reactions.jl")
include("Utils.jl")
include("IdealGas.jl")

"""
Function to create composite type of SurfaceMechanism and
SpeciesRxnMap. In general the application programs shall only
call the compile_mech function and not directly the create_mech
and species_rxn_map function
#   Usage
    compile_mech(file_path,gas_species)
-   file_path::String : path to the mechanism xml file
-   gas_species::Array{String,1} : list of gasphase species
"""
function compile_mech(file_path::T, thermo_obj::IdealGas.SpeciesThermoObj, gas_species::Array{T,1}) where T <: AbstractString
    sm = create_mech(file_path,gas_species,thermo_obj.molwt)
    rxn_map = species_rxn_map(gas_species,sm)
    return MechanismDefinition(sm,rxn_map)
end

"""
Function to read a surface reaction mechanism file
#   Usage:
    create_mech(file_path,gas_species)
-   file_path::String : file name including the relative path
-   gas_species::Array{String,1} : gasphase species list
"""
function create_mech(file_path::String, gas_species::Array{String,1},molwt::Array{Float64,1})
    xmldoc = parse_file(file_path)
    mech_root = root(xmldoc)
    mech_unit = attribute(mech_root,"unit")
    energy_factor = convert2si(mech_unit)

    #get the surface species
    surface_species = map(x->uppercase(x) ,split(content(get_elements_by_tagname(mech_root,"species")[1])))

    #default initial values
    site_coordination_vector = ones(length(surface_species))
    covg_vec = zeros(length(surface_species))
    site_density = 1.0
    site_name = ""
    #read the site node
    site_node = get_elements_by_tagname(mech_root,"site")
    for elements in site_node
        site_name = attribute(elements,"name")
        #get the site coordination
        coordination = get_elements_by_tagname(elements,"coordination")[1]
        if is_elementnode(coordination)
            site_coord_species = split(content(coordination),",")
            for items in site_coord_species
                sp,val = split(items,"=")
                sp_id = get_index(String(sp),surface_species)
                site_coordination_vector[sp_id] = parse(Float64,String(val))
            end
        end
        #get the site density
        density = get_elements_by_tagname(elements,"density")[1]
        if is_elementnode(density)
            unit = strip(attribute(density,"unit"))
            site_density = parse(Float64, content(density))
            if uppercase(unit) == "MOL/M2"
                site_density *= 1e4
            end
        else
            throw(error("Surface site density not specified\n"))
        end

        #get the intitial coverages
        covg_node = get_elements_by_tagname(elements,"initial")[1]
        if is_elementnode(covg_node)
            ini_covg = split(content(covg_node),",")
            for items in ini_covg
                sp,val = split(items,"=")
                sp_id = get_index(String(sp),surface_species)
                covg_vec[sp_id] = parse(Float64,String(val))
            end
        else
            throw(error("Initial coverage not specified\n"))
        end

    end

    site_info = SiteInfo(site_name,site_density,site_coordination_vector,covg_vec)

    #Read the coverage dependencies if any. This must be read before reading the reactions
    covg_nodes = get_elements_by_tagname(mech_root,"coverage")
    cov_dep_rxns = Dict{Int64,Dict{Int,Float64}}() #rxn_id=>{species_id=>cov_act_energy}
    for covg in covg_nodes
        rxn_ids = split(attribute(covg,"id"))
        sp,val = split(content(covg),"=")
        sp_id = get_index(String(sp),surface_species)
        cov_dict = Dict(sp_id=>parse(Float64,val))

        for rxn in rxn_ids
           cov_dep_rxns[parse(Int64,rxn)] = cov_dict
        end
    end
    #Read the order dependencies if any. This must be read before reading the reactions
    order_nodes = get_elements_by_tagname(mech_root,"order")
    order_dep_rxns = Dict{Int64,Dict{Int,Float64}}() #rxn_id=>{species_id=>order}
    for order in order_nodes
        rxn_ids = split(attribute(order,"id"))
        sp,val = split(content(order),"=")
        sp_id = get_index(String(sp),surface_species)
        order_dict = Dict(sp_id=>parse(Float64,val))
        for rxn in rxn_ids
            order_dep_rxns[parse(Int64,rxn)] = order_dict
        end
    end

    #Read Motz-Wise Correction
    mwc_node = get_elements_by_tagname(mech_root,"mwc")
    mwc_rxns = Array{Int64,1}()
    if length(mwc_node) > 0
        if is_elementnode(mwc_node[1])
            mwc_rxns = [ parse(Int64,rxn) for rxn in split(content(mwc_node[1]))]
        end
    end

    #Define the reaction array to be returned
    rxn_array= Array{ElementaryReactions,1}()

    rct_ids = Array{Int64,1}()
    prdt_ids = Array{Int64,1}()
    all_species = Array{String,1}()
    all_species = copy(gas_species)
    append!(all_species,surface_species)

    #Read sticking reactions
    site_index = get_index(site_name,all_species)

    stick_node = get_elements_by_tagname(mech_root,"stick")[1]
    if is_elementnode(stick_node)
        stick_rxns = get_elements_by_tagname(stick_node,"rxn")
        for rxn in stick_rxns
            rxn_id = parse(Int64,strip(attribute(rxn,"id")))
            stoic, params = split(content(rxn),"@")
            reversible, r_species, p_species = parase_rxn_string(String(stoic))
            rct_ids = [get_index(sp,all_species) for sp in r_species]
            prdt_ids = [get_index(sp,all_species) for sp in p_species]
            s,β,E = parse_rxn_params(String(params))
            #Check for coverage, order and mwc dependency
            cov,ord,mwc = get_dependencies(cov_dep_rxns,order_dep_rxns,mwc_rxns,rxn_id)
            #convert sticking coefficient to rate constant
            n_site = count(x->x==site_index, rct_ids)
            gas_species_id = filter(x->x != site_index, rct_ids)[1]
            sqrt_term = 100*sqrt(R/(2*π*molwt[gas_species_id]))
            s_by_gamma = s/site_density^n_site
            k = sqrt_term*s_by_gamma #rate constant
            #Apply the MWC
            if mwc
                k /= (1-0.5*s)
            end
            p = Parameters(s,k,β,E)

            stoic = Stoichiometry(reversible,rct_ids,prdt_ids,cov,ord)
            elm_rxn = ElementaryReactions(rxn_id,stick,p,mwc,stoic)
            push!(rxn_array,elm_rxn)

        end
    end

    #Read Arrhenius reactions
    arrhenius_node = get_elements_by_tagname(mech_root,"arrhenius")[1]
    if is_elementnode(arrhenius_node)
        arrhenius_rxns = get_elements_by_tagname(arrhenius_node,"rxn")
        for rxn in arrhenius_rxns
            rxn_id = parse(Int64,strip(attribute(rxn,"id")))
            stoic,params = split(content(rxn),"@")
            reversible, r_species, p_species = parase_rxn_string(String(stoic))
            rct_ids = [get_index(sp,all_species) for sp in r_species]
            prdt_ids = [get_index(sp,all_species) for sp in p_species]
            k,β,E = parse_rxn_params(String(params))
            p = Parameters(0.0,k,β,E)
            #Check for coverage, order and mwc dependency. MWC not applicable to Arrhenius type
            cov,ord,mwc = get_dependencies(cov_dep_rxns,order_dep_rxns,mwc_rxns,rxn_id)
            stoic = Stoichiometry(reversible,rct_ids,prdt_ids,cov,ord)
            elm_rxn = ElementaryReactions(rxn_id,arrhenius,p,mwc,stoic)
            push!(rxn_array,elm_rxn)
        end
    end

    return SurfaceMechanism(energy_factor,site_info,surface_species,rxn_array)
end

"""
Create array of structures of SpeciesRxnMap. The SpeciesRxnMap is a
composit structure of species_id in participating reactions and the
stoichiometric coefficients
#   Usage
    species_rxn_map(gas_species,sm)
-   gas_species::Array{String,1} :  list of gas phase species
-   sm::SurfaceMechanism : SurfaceMechanism composite type
"""
function species_rxn_map(gas_species::Array{String,1},sm::SurfaceMechanism)
    #Array of SpeciesRxnMap
    srm_array = Array{SpeciesRxnMap,1}()
    #create list of all species
    all_species = Array{String,1}()
    all_species = copy(gas_species)
    append!(all_species,sm.species)
    #loop over all species
    for sp in all_species
        sp_id = get_index(sp,all_species)
        #check if the species is present in the reactant_ids of a reaction
        rxns = Array{Int64,1}()
        stc = Array{Int64,1}()
        for rxn in sm.reactions
            #count the number of occurances of species in reactants list (using ids)
            n = count(x->x==sp_id,rxn.rxn_stoic.reactant_ids)
            if n > 0
                push!(rxns,rxn.id)
                push!(stc,-n)
            end
            #count the number of occurances of species in the products list (using ids)
            n = count(x->x==sp_id,rxn.rxn_stoic.product_ids)
            if n > 0
                push!(rxns,rxn.id)
                push!(stc,n)
            end
        end
        push!(srm_array,SpeciesRxnMap(sp_id,rxns,stc))
    end
    return srm_array
end

"""
Function to calculate the molar production rate of gas-phase and
surface species based on the surface reaction mechanism. The function
will alter the field s_rate after the calculation.
#   Usage:
    calculate_molar_production_rate!(state,md)
-   state::State : state of the reaction, which is of the composite type State
-   thermo_obj::SpeciesThermoObj : Composite type SpeciesThermoObj. Required in case of reversible reactions
-   md::MechanismDefinition : the mechanism, which is of the composite type MechanismDefinition
"""
function calculate_molar_production_rates!(state::ReactionState, thermo_obj::IdealGas.SpeciesThermoObj,md::SurfaceMechDefinition)
    #get the mole fractions from state
    mf = state.mole_frac
    #convert mole fractions to concentrations
    all_conc = (state.p/R/state.T) .* mf
    all_conc *= 1e-6 # convert to mol/cm3
    #get the surface coverages from state
    coverage = state.covg
    #convert coverages to surface concentrations mol/cm2
    surf_conc =  (md.sm.si.density .* coverage) ./ md.sm.si.site_coordination
    #merge the concentrations
    append!(all_conc,surf_conc)
    rxn_rate = Array{Float64,1}()

    #calcuate the rate of individual reactions
    for rxn in md.sm.reactions
        local prdt_conc = prod(all_conc[rxn.rxn_stoic.reactant_ids])
        #check for coverage dependency of the reaction
        local cov_act_E = 0
        #cov_act_E = sum(coverage[collect(keys(rxn.rxn_stoic.cov_dep))] .* collect(values(rxn.rxn_stoic.cov_dep)))
        if !isempty(rxn.rxn_stoic.cov_dep)
            #calculate coverage dependent activation energy
            for (sp_id,cov_e) in rxn.rxn_stoic.cov_dep
                cov_act_E += coverage[sp_id]*cov_e
            end
        end
        #check for order dependency
        local m_order = 1.0
        if !isempty(rxn.rxn_stoic.order_dep)
            #calculate order change
            for (sp_id,order) in rxn.rxn_stoic.order_dep
                m_order *= coverage[sp_id]^order
            end
        end
        #calculate the forward reaction rate
        E = (rxn.params.E + cov_act_E)*md.sm.energy_factor
        k = rxn.params.k0 * state.T^rxn.params.β * exp(-E/R/state.T)
        if rxn.type == stick
            k *= sqrt(state.T)
        end
        push!(rxn_rate,k*prdt_conc*m_order*1e4)#e4 is to convert the units to mol/m2

        if rxn.rxn_stoic.reversible
            throw(error("Reverse reaction not implemented\n"))
        end

    end
    #calculate the species molar production rates
    for srm in md.srm_array
        state.s_rate[srm.species_id] = sum(rxn_rate[srm.rxns] .* srm.stoic_coeff)
    end
end

"""
Function to calculate the molar production rate of gas-phase species.
This function integrates the rate equations so that \frac{d\theta}{dt}=0
#   Usage:
    calculate_molar_production_rate!(state,md)
-   state::State : state of the reaction, which is of the composite type State
-   thermo_obj::SpeciesThermoObj : Composite type SpeciesThermoObj. Required in case of reversible reactions
-   md::MechanismDefinition : the mechanism, which is of the composite type MechanismDefinition
"""
function calculate_ss_molar_production_rates!(state::ReactionState, thermo_obj::IdealGas.SpeciesThermoObj,md::SurfaceMechDefinition, time =1.0)
    t_span = (0,time)
    coverage = state.covg
    p = (state,thermo_obj,md)
    prob = ODEProblem(covg_integration!,coverage,t_span,p)
    sol = solve(prob, alg_hints=[:stiff] , reltol=1e-6, abstol=1e-10, save_everystep=false,save_start=false)
    state.covg = sol.u[lastindex(sol.u)]
    return sol.t, state
end

function covg_integration!(du, u, p, t)
    state, thermo_obj,md = p
    state.covg = u
    ng = length(state.mole_frac)
    calculate_molar_production_rates!(state,thermo_obj,md)
    @inbounds for i in eachindex(u)
        du[i] = (md.sm.si.site_coordination[i]*state.s_rate[ng+i])/(md.sm.si.density*1e4)
    end
end

"""
Function to update the rate parameters. This function is only used for
sensitivity analysis
#   Usage
    update_params(rate_params,md)
-   rate_params::Array  : One dimentional array of floating point variables
-   md:SurfaceMechDefinition : SurfaceMechDefinition struct
"""
function update_params(rate_params::Array{Float64,1},md::SurfaceMechDefinition)
    #assign the new parameters to SurfaceMechDefinition object
    for i in eachindex(rate_params)
        if Int(md.sm.reactions[i].type) == Int(stick)
            k = md.sm.reactions[i].params.k0 #This is the current value of converted rate constant including MWC
            if md.sm.reactions[i].mwc == true
                #take out the MWC factor
                k *= (1-0.5*md.sm.reactions[i].s)
            end
            #new k value for the purturbed sticking coefficient
            k = (k/md.sm.reactions[i].params.s)*rate_params[i]
            md.sm.reactions[i].params.s = rate_params[i]
            if md.sm.reactions[i].mwc ==  true
                #Apply the MWC back
                k /= (1-0.5*rate_params[i])
            end
            md.sm.reactions[i].params.k0 = k
        else
            md.sm.reactions[i].params.k0 = rate_params[i]
        end
    end

end

"""
A function for use in global sensitivity analysis. The function updates the parameters with the
    supplied rate_params
#   Usage
    update_params!(md, rate_params,rxn_ids, constraint_ids)
-   md : SurfaceMechDefinition structure
-   rate_params : parameters to that needs to be set to the mechanism object
-   rxn_ids : id of the reactions considered in the sensitivity analysis
-   constraint_ids : id of the reverse reaction ids considered
"""
function update_params!(md, rate_params::Array{Float64,1},rxn_ids::Array{Int64,1}, constraint_ids::Array{Int64,1}, pratio::Array{Float64,1})
    for id in eachindex(rxn_ids)
        if Int(md.sm.reactions[rxn_ids[id]].type) == Int(stick)
            # The current value of sticking coefficient converted to rate constant including MWC
            k = md.sm.reactions[rxn_ids[id]].params.k0
            if md.sm.reactions[rxn_ids[id]].mwc == true
                #take out the MWC factor
                k *= (1-0.5*md.sm.reactions[rxn_ids[id]].params.s)
            end
            # multiplication factor used to convert sticking coefficient to rate constant
            m_factor = k/md.sm.reactions[rxn_ids[id]].params.s
            # new sticking coefficient
            md.sm.reactions[rxn_ids[id]].params.s = rate_params[id]
            #new k value for the purturbed sticking coefficient
            k  = md.sm.reactions[rxn_ids[id]].params.s * m_factor
            #Apply the MWC back
            if md.sm.reactions[id].mwc == true
                k /= (1-0.5*rate_params[id])
            end
            md.sm.reactions[rxn_ids[id]].params.k0 = k

            # Now if there is a constraint, i.e. if a reverse reaction present, then the
            # parameter for that reaction needs to be changed accordingly.
            if constraint_ids[id] != 0
                md.sm.reactions[constraint_ids[id]].params.k0 = k/pratio[id]
            end
        else
            md.sm.reactions[rxn_ids[id]].params.k0 = rate_params[id]
            if constraint_ids[id] != 0
                md.sm.reactions[constraint_ids[id]].params.k0 = rate_params[id]/pratio[id]
            end
        end
    end
end

#end of module
end


text:
"""
Given a remote repo URL and a git tree spec, get a `Project` object
for the project file in that tree and a hash string for the tree.
"""
# function get_project(remote_url::AbstractString, tree_spec::AbstractString)
#     # TODO?: use raw file downloads for GitHub/GitLab
#     mktempdir(mkpath("packages")) do tmp
#         # bare clone the package repo
#         @debug("bare clone the package repo")
#         repo = LibGit2.clone(remote_url, joinpath(tmp, "repo"), isbare=true)
#         tree = try
#             LibGit2.GitObject(repo, tree_spec)
#         catch err
#             err isa LibGit2.GitError && err.code == LibGit2.Error.ENOTFOUND || rethrow(err)
#             error("$remote_url: git object $(repr(tree_spec)) could not be found")
#         end
#         tree isa LibGit2.GitTree || (tree = LibGit2.peel(LibGit2.GitTree, tree))
#
#         # check out the requested tree
#         @debug("check out the requested tree")
#         tree_path = abspath(tmp, "tree")
#         GC.@preserve tree_path begin
#             opts = LibGit2.CheckoutOptions(
#                 checkout_strategy = LibGit2.Consts.CHECKOUT_FORCE,
#                 target_directory = Base.unsafe_convert(Cstring, tree_path)
#             )
#             LibGit2.checkout_tree(repo, tree, options=opts)
#         end
#
#         # look for a project file in the tree
#         @debug("look for a project file in the tree")
#         project_file = Pkg.Types.projectfile_path(tree_path)
#         project_file !== nothing && isfile(project_file) ||
#             error("$remote_url: git tree $(repr(tree_spec)) has no project file")
#
#         # parse the project file
#         @debug("parse the project file")
#         project = Pkg.Types.read_project(project_file)
#         project.name === nothing &&
#             error("$remote_url $(repr(tree_spec)): package has no name")
#         project.uuid === nothing &&
#             error("$remote_url $(repr(tree_spec)): package has no UUID")
#         project.version === nothing &&
#             error("$remote_url $(repr(tree_spec)): package has no version")
#
#         return project, string(LibGit2.GitHash(tree))
#     end
# end

# These can compromise the integrity of the registry and cannot be
# opted out of.
const mandatory_errors = [:version_exists,
                          :change_package_name,  # not implemented
                          :change_package_uuid,  # not allowed
                          :package_self_dep,
                          :name_mismatch,
                          :wrong_stdlib_uuid,
                          :package_url_missing
                          ]

# These are considered errors by Registrator and are default for the
# `register` function. The caller can override this selection.
const registrator_errors = [:version_zero,
                            :version_less_than_all_existing,
                            :change_package_url,
                            :dependency_not_found,
                            :julia_before_07_in_compat,
                            :invalid_compat,
                            :unexpected_registration_error
                            ]

mutable struct ReturnStatus
    triggered_checks::Vector{NamedTuple}
    errors::Set{Symbol}
    error_found::Bool
end

function ReturnStatus(errors::Vector{Symbol} = Symbol[])
    ReturnStatus(NamedTuple[], union!(Set(mandatory_errors), errors), false)
end

function add!(status::ReturnStatus, check::Symbol,
              data::Union{Nothing, NamedTuple}=nothing)
    if isnothing(data)
        push!(status.triggered_checks, (id = check,))
    else
        push!(status.triggered_checks, merge((id = check,), data))
    end
    status.error_found |= check in status.errors
end

haserror(status::ReturnStatus) = status.error_found

struct RegBranch
    name::String
    version::VersionNumber
    branch::String

    metadata::Dict{String, Any} # "error", "warning", kind etc.

    function RegBranch(pkg::Pkg.Types.Project, branch::AbstractString)
        new(pkg.name, pkg.version, branch, Dict{String,Any}())
    end
end

function add_label!(regbr::RegBranch, label)
    if !haskey(regbr.metadata, "labels")
        regbr.metadata["labels"] = String[]
    end
    pushfirst!(regbr.metadata["labels"], label)
end

# error in regbr.metadata["errors"]
# warning in regbr.metadata["warning"]
# version labels for the PR in in regbr.metadata["labels"]
function set_metadata!(regbr::RegBranch, status::ReturnStatus)
    empty!(regbr.metadata)
    for triggered_check in reverse(status.triggered_checks)
        check = triggered_check.id
        data = triggered_check
        if check == :version_zero
            regbr.metadata["error"] = "Package version must be greater than 0.0.0"
        elseif check == :new_package_label
            add_label!(regbr, "new package")
        elseif check == :not_standard_first_version
            regbr.metadata["warning"] =
                """This looks like a new registration that registers version $(data.version).
                Ideally, you should register an initial release with 0.0.1, 0.1.0 or 1.0.0 version numbers"""
        elseif check == :version_less_than_all_existing
            regbr.metadata["error"] = "Version $(data.version) less than least existing version $(data.least)"
        elseif check == :version_exists
            regbr.metadata["error"] = "Version $(data.version) already exists"
        elseif check == :major_release
            add_label!(regbr, "major_release")
        elseif check == :minor_release
            add_label!(regbr, "minor_release")
        elseif check == :patch_release
            add_label!(regbr, "patch_release")
        elseif check == :breaking
            add_label!(regbr, "BREAKING")
        elseif check == :version_skip
            regbr.metadata["warning"] = "Version $(data.version) skips over $(data.next)"
        elseif check == :change_package_name
            regbr.metadata["error"] = "Changing package names not supported yet"
        elseif check == :change_package_url
            regbr.metadata["error"] = "Changing package repo URL not allowed, please submit a pull request with the URL change to the target registry and retry."
        elseif check == :new_version
            regbr.metadata["kind"] = "New version"
        elseif check == :new_package
            regbr.metadata["kind"] = "New package"
        elseif check == :change_package_uuid
            regbr.metadata["error"] = "Changing UUIDs is not allowed"
        elseif check == :package_self_dep
            regbr.metadata["error"] = "Package $(data.name) mentions itself in `[deps]`"
        elseif check == :name_mismatch
            regbr.metadata["error"] = "Error in (Julia)Project.toml: UUID $(data.uuid) refers to package '$(data.reg_name)' in registry but Project.toml has '$(data.project_name)'"
        elseif check == :wrong_stdlib_uuid
            regbr.metadata["error"] = "Error in (Julia)Project.toml: UUID $(data.project_uuid) for package $(data.name) should be $(data.stdlib_uuid)"
        elseif check == :dependency_not_found
            regbr.metadata["error"] = "Error in (Julia)Project.toml: Package '$(data.name)' with UUID: $(data.uuid) not found in registry or stdlib"
        elseif check == :julia_before_07_in_compat
            regbr.metadata["error"] = "Julia version < 0.7 not allowed in `[compat]`"
        elseif check == :invalid_compat
            regbr.metadata["error"] = "Following packages are mentioned in `[compat]` but not found in `[deps]` or `[extras]`:\n" * join(data.invalid_compats, "\n")
        elseif check == :package_url_missing
            regbr.metadata["error"] = "No repo URL provided for a new package"
        elseif check == :unexpected_registration_error
            regbr.metadata["error"] = "Unexpected error in registration"
        end
    end
    return regbr
end

get_backtrace(ex) = sprint(Base.showerror, ex, catch_backtrace())

function write_registry(registry_path::AbstractString, reg::RegistryData)
    open(registry_path, "w") do io
        TOML.print(io, reg)
    end
end

function check_version!(version::VersionNumber, existing::Vector{VersionNumber},
                        status::ReturnStatus)
    if version == v"0"
        add!(status, :version_zero)
        haserror(status) && return
    end

    @assert issorted(existing)
    if isempty(existing)
        add!(status, :new_package_label)
        if !(version in [v"0.0.1", v"0.1", v"1"])
            add!(status, :not_standard_first_version, (version = version,))
        end
    else
        idx = searchsortedlast(existing, version)
        if idx <= 0
            add!(status, :version_less_than_all_existing,
                 (version = version, least = first(existing)))
            return
        end

        previous = existing[idx]
        if version == previous
            add!(status, :version_exists, (version = version,))
            haserror(status) && return
        end
        if version.major != previous.major
            add!(status, :major_release)
            add!(status, :breaking)
            next = VersionNumber(previous.major + 1, 0, 0)
        elseif version.minor != previous.minor
            add!(status, :minor_release)
            if version.major == 0
                add!(status, :breaking)
            end
            next = VersionNumber(previous.major, previous.minor + 1, 0)
        else
            add!(status, :patch_release)
            if version.major == version.minor == 0
                add!(status, :breaking)
            end
            next = VersionNumber(previous.major, previous.minor, previous.patch + 1)
        end
        if version > next
            add!(status, :version_skip, (version = version, next = next))
        end
    end

    return
end

findpackageerror!(name::AbstractString, uuid::Base.UUID,
                  regdata::Array{RegistryData}, status::ReturnStatus) =
    findpackageerror!(name, string(uuid), regdata, status)

# Check that `uuid` is found in `regdata` OR `name` is found in
# `BUILTIN_PKGS` (i.e. stdlibs). In the former case, check that `name`
# matches the found uuid and in the latter case that `uuid` matches
# the found name.
function findpackageerror!(name::AbstractString, uuid::AbstractString,
                           regdata::Array{RegistryData}, status::ReturnStatus)
    for registry_data in regdata
        if haskey(registry_data.packages, uuid)
            name_in_reg = registry_data.packages[uuid]["name"]
            if name_in_reg != name
                @debug(:name_mismatch)
                add!(status, :name_mismatch,
                     (uuid = uuid, reg_name = name_in_reg, project_name = name))
            end
            return
        end
    end

    if haskey(BUILTIN_PKGS, name)
        if BUILTIN_PKGS[name] != uuid
            @debug(:wrong_stdlib_uuid)
            add!(status, :wrong_stdlib_uuid,
                 (project_uuid = uuid, name = name,
                  stdlib_uuid = BUILTIN_PKGS[name]))
        end
    else
        @debug(:dependency_not_found)
        add!(status, :dependency_not_found, (name = name, uuid = uuid))
    end

    return
end

import Pkg.Types: VersionRange, VersionBound, VersionSpec

function versionrange(lo::VersionBound, hi::VersionBound)
    lo.t == hi.t && (lo = hi)
    return VersionRange(lo, hi)
end

function find_package_in_registry(pkg::Pkg.Types.Project,
                                  registry_file::AbstractString,
                                  registry_path::AbstractString,
                                  registry_data::RegistryData,
                                  status::ReturnStatus)
    uuid = string(pkg.uuid)
    if haskey(registry_data.packages, uuid)
        package_data = registry_data.packages[uuid]
        if package_data["name"] != pkg.name
            err = :change_package_name
            @debug(err)
            add!(status, err)
            haserror(status) && return nothing
        end
        package_path = joinpath(registry_path, package_data["path"])
        add!(status, :new_version)
    else
        @debug("Package with UUID: $uuid not found in registry, checking if UUID was changed")
        for (k, v) in registry_data.packages
            if v["name"] == pkg.name
                err = :change_package_uuid
                @debug(err)
                add!(status, err)
                haserror(status) && return nothing
            end
        end

        @debug("Creating directory for new package $(pkg.name)")
        package_path = joinpath(registry_path, package_relpath(pkg))
        mkpath(package_path)

        @debug("Adding package UUID to registry")
        push!(registry_data, pkg)
        write_registry(registry_file, registry_data)
        add!(status, :new_package)
    end

    return package_path
end

# If `package_repo` is an empty string, replace it with what is
# already stored in Package.toml.
function check_package!(package_repo::AbstractString,
                        package_path::AbstractString,
                        status::ReturnStatus)
    package_file = joinpath(package_path, "Package.toml")
    # If this is a registration of a new package, the package file has
    # not been created yet.
    if isfile(package_file)
        repo = TOML.parsefile(package_file)["repo"]
        if isempty(package_repo)
            package_repo = repo
        elseif repo != package_repo
            err = :change_package_url
            @debug(err)
            add!(status, err)
        end
    elseif isempty(package_repo)
        add!(status, :package_url_missing)
    end

    return package_repo
end

function update_package_file(pkg::Pkg.Types.Project,
                             package_repo::AbstractString,
                             package_path::AbstractString)
    package_info = Dict("name" => pkg.name,
                        "uuid" => string(pkg.uuid),
                        "repo" => package_repo)
    package_file = joinpath(package_path, "Package.toml")
    open(package_file, "w") do io
        TOML.print(io, package_info; sorted=true,
            by = x -> x == "name" ? 1 : x == "uuid" ? 2 : 3)
    end
    nothing
end

function get_versions_file(package_path::AbstractString)
    filename = joinpath(package_path, "Versions.toml")
    data = isfile(filename) ? TOML.parsefile(filename) : Dict{String, Any}()
    return filename, data
end

function check_versions!(pkg::Pkg.Types.Project,
                         versions_data::Dict{String, Any},
                         status::ReturnStatus)
    versions = sort!([VersionNumber(v) for v in keys(versions_data)])
    check_version!(pkg.version, versions, status)
    return versions
end

function update_versions_file(pkg::Pkg.Types.Project,
                              versions_file::AbstractString,
                              versions_data::Dict{String, Any},
                              tree_hash::AbstractString)
    version_info = Dict{String, Any}("git-tree-sha1" => string(tree_hash))
    versions_data[string(pkg.version)] = version_info

    open(versions_file, "w") do io
        # TOML.print with sorted=true sorts recursively
        # so this by function needs to handle the outer dict
        # with version number keys, and the inner dict with
        # git-tree-sha1, yanked, etc as keys.
        function by(x)
            if occursin(Base.VERSION_REGEX, x)
                return VersionNumber(x)
            else
                if x == "git-tree-sha1"
                    return 1
                elseif x == "yanked"
                    return 2
                else
                    return 3
                end
            end
        end
        TOML.print(io, versions_data; sorted=true, by=by)
    end
end

function check_deps!(pkg::Pkg.Types.Project,
                     regdata::Vector{RegistryData},
                     status::ReturnStatus)
    if pkg.name in keys(pkg.deps)
        err = :package_self_dep
        @debug(err)
        add!(status, err, (name = pkg.name,))
        haserror(status) && return
    end

    @debug("Verifying package name and uuid in deps")
    for (name, uuid) in pkg.deps
        findpackageerror!(name, uuid, regdata, status)
    end
end

# Note that Compress.load should load with respect to Versions.toml
# before update and Compress.save should save with respect to
# Versions.toml after update. This is handled with the `old_versions'
# argument and the assumption that Versions.toml has been updated with
# the new version before calling this function.
function update_deps_file(pkg::Pkg.Types.Project,
                          package_path::AbstractString,
                          old_versions::Vector{VersionNumber})
    deps_file = joinpath(package_path, "Deps.toml")
    if isfile(deps_file)
        deps_data = Compress.load(deps_file, old_versions)
    else
        deps_data = Dict()
    end

    deps_data[pkg.version] = pkg.deps
    Compress.save(deps_file, deps_data)
end

function check_compat!(pkg::Pkg.Types.Project,
                       regdata::Vector{RegistryData},
                       regpaths::Vector{String},
                       status::ReturnStatus)
    if haskey(pkg.compat, "julia")
        ver = Pkg.Types.semver_spec(pkg.compat["julia"])
        if any(map(x -> !isempty(intersect(Pkg.Types.VersionRange("0-0.6"), x)), ver.ranges))
            err = :julia_before_07_in_compat
            @debug(err)
            add!(status, err)
            haserror(status) && return
        end
    end

    # Note: These checks are meaningless for Julia >= 1.2 since
    # Pkg.Types.load_project will give an error if there are compat
    # entries not mentioned in deps, nor in extras.
    invalid_compats = []
    for name in keys(pkg.compat)
        indeps = haskey(pkg.deps, name)
        inextras = haskey(pkg.extras, name)
        if !indeps && !inextras && name != "julia"
            push!(invalid_compats, name)
        end
    end
    if !isempty(invalid_compats)
        err = :invalid_compat
        @debug(err)
        add!(status, err, (invalid_compats = invalid_compats,))
        haserror(status) && return
    end

    # Note: `findpackageerror` has already been run for all entries in
    # deps so doing it again for the intersection of compat and deps
    # is redundant. Doing it for the intersection of compat and extras
    # is meaningful but it is unclear why not just do it for all
    # entries in extras at the same time it's done for all entries in
    # deps. Alternatively it can be skipped entirely since nothing
    # that is in extras only will be stored in the registry files
    # anyway.
    for name in keys(pkg.compat)
        if name != "julia"
            indeps = haskey(pkg.deps, name)
            inextras = haskey(pkg.extras, name)

            if indeps
                uuidofdep = string(pkg.deps[name])
                findpackageerror!(name, uuidofdep, regdata, status)
            elseif inextras
                uuidofdep = string(pkg.extras[name])
                findpackageerror!(name, uuidofdep, regdata, status)
            end

            haserror(status) && return
        end
    end

    return
end

# See the comments for `update_deps_file` for the rationale for the
# `old_versions` argument.
function update_compat_file(pkg::Pkg.Types.Project,
                            package_path::AbstractString,
                            old_versions::Vector{VersionNumber})
    @debug("update package data: compat file")
    compat_file = joinpath(package_path, "Compat.toml")
    if isfile(compat_file)
        compat_data = Compress.load(compat_file, old_versions)
    else
        compat_data = Dict()
    end

    d = Dict()
    for (name, version) in pkg.compat
        if !haskey(pkg.deps, name) && name != "julia"
            @debug("$name is a test-only dependency; omitting from Compat.toml")
            continue
        end

        spec = Pkg.Types.semver_spec(version)
        # The call to `map(versionrange, )` can be removed
        # once Pkg is updated to a version including
        # https://github.com/JuliaLang/Pkg.jl/pull/1181
        # and support for older versions is dropped.
        ranges = map(r->versionrange(r.lower, r.upper), spec.ranges)
        ranges = VersionSpec(ranges).ranges # this combines joinable ranges
        d[name] = length(ranges) == 1 ? string(ranges[1]) : map(string, ranges)
    end

    compat_data[pkg.version] = d
    Compress.save(compat_file, compat_data)
end

function get_registrator_tree_sha()
    regtreesha = nothing
    reg_pkgs = Pkg.Display.status(Pkg.Types.Context(),
                   [Pkg.PackageSpec("Registrator", Base.UUID("4418983a-e44d-11e8-3aec-9789530b3b3e"))])
    if !isempty(reg_pkgs)
        regtreesha = reg_pkgs[1].new.hash
    end
    if regtreesha === nothing
        regpath = abspath(joinpath(@__DIR__, "..", ".."))
        if isdir(joinpath(regpath, ".git"))
            regtreesha = LibGit2.head(regpath)
        else
            regtreesha = "unknown"
        end
    end

    return regtreesha
end

function check_and_update_registry_files(pkg, package_repo, tree_hash,
                                         registry_path, registry_deps_paths,
                                         status)
    # find package in registry
    @debug("find package in registry")
    registry_file = joinpath(registry_path, "Registry.toml")
    registry_data = parse_registry(registry_file)
    package_path = find_package_in_registry(pkg, registry_file, registry_path,
                                            registry_data, status)
    haserror(status) && return

    # update package data: package file
    @debug("update package data: package file")
    package_repo = check_package!(package_repo, package_path, status)
    haserror(status) && return
    update_package_file(pkg, package_repo, package_path)

    # update package data: versions file
    @debug("update package data: versions file")
    versions_file, versions_data = get_versions_file(package_path)
    old_versions = check_versions!(pkg, versions_data, status)
    haserror(status) && return
    update_versions_file(pkg, versions_file, versions_data, tree_hash)

    # update package data: deps file
    @debug("update package data: deps file")
    registry_deps_data = map(registry_deps_paths) do registry_path
        parse_registry(joinpath(registry_path, "Registry.toml"))
    end
    regdata = [registry_data; registry_deps_data]
    check_deps!(pkg, regdata, status)
    haserror(status) && return
    update_deps_file(pkg, package_path, old_versions)

    # update package data: compat file
    @debug("check compat section")
    regpaths = [registry_path; registry_deps_paths]
    check_compat!(pkg, regdata, regpaths, status)
    haserror(status) && return
    update_compat_file(pkg, package_path, old_versions)
end

"""
    register(package_repo, pkg, tree_hash; registry, registry_deps, push, gitconfig)

Register the package at `package_repo` / `tree_hash` in `registry`.
Returns a `RegEdit.RegBranch` which contains information about the registration and/or any
errors or warnings that occurred.

# Arguments

* `package_repo::AbstractString`: The git repository URL for the package to be registered. If empty, keep the stored repository URL.
* `pkg::Pkg.Types.Project`: the parsed (Julia)Project.toml file for the package to be registered
* `tree_hash::AbstractString`: the tree hash (not commit hash) of the package revision to be registered

# Keyword Arguments

* `registry::AbstractString="$DEFAULT_REGISTRY_URL"`: the git repository URL for the registry
* `registry_deps::Vector{String}=[]`: the git repository URLs for any registries containing
    packages depended on by `pkg`
* `push::Bool=false`: whether to push a registration branch to `registry` for consideration
* `gitconfig::Dict=Dict()`: dictionary of configuration options for the `git` command
"""
function register(
    package_repo::AbstractString, pkg::Pkg.Types.Project, tree_hash::AbstractString;
    registry::AbstractString = DEFAULT_REGISTRY_URL,
    registry_deps::Vector{<:AbstractString} = AbstractString[],
    checks_triggering_error = registrator_errors,
    push::Bool = false,
    force_reset::Bool = true,
    branch::String = registration_branch(pkg),
    cache::RegistryCache=REGISTRY_CACHE,
    gitconfig::Dict = Dict()
)
    # get info from package registry
    @debug("get info from package registry")
    if !isempty(package_repo)
        package_repo = GitTools.normalize_url(package_repo)
    end

    # return object
    regbr = RegBranch(pkg, branch)

    # status object
    status = ReturnStatus(checks_triggering_error)

    # get up-to-date clone of registry
    @debug("get up-to-date clone of registry")
    registry = GitTools.normalize_url(registry)
    registry_repo = get_registry(registry; gitconfig=gitconfig, force_reset=force_reset, cache=cache)
    registry_path = LibGit2.path(registry_repo)

    isempty(registry_deps) || @debug("get up-to-date clones of registry dependencies")
    registry_deps_paths = map(registry_deps) do registry
        LibGit2.path(get_registry(GitTools.normalize_url(registry); gitconfig=gitconfig, force_reset=force_reset, cache=cache))
    end

    clean_registry = true
    err = nothing
    try
        # branch registry repo
        @debug("branch registry repo")
        git = gitcmd(registry_path, gitconfig)
        run(pipeline(`$git checkout -f master`; stdout=devnull))
        if branch != "master"
            run(pipeline(`$git branch -f $branch`; stdout=devnull))
            run(pipeline(`$git checkout -f $branch`; stdout=devnull))
        end

        check_and_update_registry_files(pkg, package_repo, tree_hash,
                                        registry_path, registry_deps_paths,
                                        status)
        haserror(status) && return set_metadata!(regbr, status)

        regtreesha = get_registrator_tree_sha()
        # To get "kind" information. Otherwise redundant to do it here.
        set_metadata!(regbr, status)

        # commit changes
        @debug("commit changes")
        message = """
        $(regbr.metadata["kind"]): $(pkg.name) v$(pkg.version)

        UUID: $(pkg.uuid)
        Repo: $(package_repo)
        Tree: $(string(tree_hash))

        Registrator tree SHA: $(regtreesha)
        """
        registry_file = joinpath(registry_path, "Registry.toml")
        package_path = joinpath(registry_path, package_relpath(pkg))
        run(pipeline(`$git add -- $package_path`; stdout=devnull))
        run(pipeline(`$git add -- $registry_file`; stdout=devnull))
        run(pipeline(`$git commit -m $message`; stdout=devnull))

        # push -f branch to remote
        if push
            @debug("push -f branch to remote")
            run(pipeline(`$git push -f -u origin $branch`; stdout=devnull))
        else
            @debug("skipping git push")
        end

        clean_registry = false
    catch ex
        println(get_backtrace(ex))
        add!(status, :unexpected_registration_error)
    finally
        if clean_registry
            @debug("cleaning up possibly inconsistent registry", registry_path=showsafe(registry_path), err=showsafe(err))
            rm(registry_path; recursive=true, force=true)
        end
    end
    return set_metadata!(regbr, status)
end

struct RegisterParams
    package_repo::AbstractString
    pkg::Pkg.Types.Project
    tree_sha::AbstractString
    registry::AbstractString
    registry_deps::Vector{<:AbstractString}
    push::Bool
    gitconfig::Dict

    function RegisterParams(package_repo::AbstractString,
                            pkg::Pkg.Types.Project,
                            tree_sha::AbstractString;
                            registry::AbstractString=DEFAULT_REGISTRY_URL,
                            registry_deps::Vector{<:AbstractString}=[],
                            push::Bool=false,
                            gitconfig::Dict=Dict())
        new(package_repo, pkg, tree_sha,
            registry, registry_deps,
            push, gitconfig)
    end
end

register(regp::RegisterParams) = register(regp.package_repo, regp.pkg, regp.tree_sha;
                                          registry=regp.registry, registry_deps=regp.registry_deps,
                                          push=regp.push, gitconfig=regp.gitconfig)

"""
    find_registered_version(pkg, registry_path)

If the package and version specified by `pkg` exists in the registry
at `registry_path`, return its tree hash. Otherwise return the empty
string.
"""
function find_registered_version(pkg::Pkg.Types.Project,
                                 registry_path::AbstractString)
    registry_file = joinpath(registry_path, "Registry.toml")
    registry_data = parse_registry(registry_file)
    # Cannot use find_package_in_registry since it may add paths in
    # the registry.
    if !haskey(registry_data.packages, string(pkg.uuid))
        return ""
    end
    package_data = registry_data.packages[string(pkg.uuid)]
    package_path = joinpath(registry_path, package_data["path"])
    _, versions_data = get_versions_file(package_path)
    if !haskey(versions_data, string(pkg.version))
        return ""
    end
    return versions_data[string(pkg.version)]["git-tree-sha1"]
end