text:
function multisolve(m::Model, mdata::MultiData, ::WeightedSum, ::LinearProblem)
    objectives = mdata.objectives
    sensemap = Dict(MOI.MIN_SENSE => 1.0, MOI.MAX_SENSE => -1.0)

    vararr = all_variables(m)

    numobj = length(objectives)
    Phi = zeros(numobj, numobj)

    # Individual minimisations
    for (i, objective) in enumerate(objectives)
        @objective(m, objective.sense, objective.f)
        for (key, value) in objective.initialvalue
            set_start_value(variable_by_name(m, key), value)
        end
        optimize!(m, ignore_optimize_hook=true);
        if !(termination_status(m) in [MOI.OPTIMAL, MOI.LOCALLY_SOLVED])
            return termination_status(m)
        end

        push!(mdata.utopiavarvalues, [value(var) for var in vararr])
        push!(mdata.paretovarvalues, [value(var) for var in vararr])

        Phi[:,i] = sensevalue.(objectives)

        push!(mdata.paretofront, value.(objectives))
    end
    Fmax = maximum(Phi, dims=2)
    Fmin = minimum(Phi, dims=2) # == diag(Phi)?
    if Fmax == Fmin
        error("The Nadir and Utopia points are equal") # I think that's what this means?
    end

    mdata.Phi = Phi

    beta = zeros(numobj); beta[end] = 1.0

    betatree = betas(numobj, mdata.pointsperdim - 1)

    for betaval in betatree
        if count(t -> t != 0, betaval) == 1
            # Skip individual optimisations as
            # they are already performed
            continue
        end
        @objective(m, MOI.MIN_SENSE,
            sum(betaval[i] * (sensemap[objectives[i].sense] * objectives[i].f -
            Fmin[i]) / (Fmax[i] - Fmin[i]) for i in Base.OneTo(numobj))
        )

        optimize!(m, ignore_optimize_hook=true);
        if !(termination_status(m) in [MOI.OPTIMAL, MOI.LOCALLY_SOLVED])
            return termination_status(m)
        end

        push!(mdata.paretofront, value.(objectives))
        push!(mdata.paretovarvalues, [value(var) for var in vararr])
    end
    return MOI.OPTIMAL
end

function multisolve(m::Model, mdata::MultiData, ::EpsilonCons, ::LinearProblem)
    objectives = mdata.objectives
    sensemap = Dict(MOI.MIN_SENSE => 1.0, MOI.MAX_SENSE => -1.0)

    vararr = all_variables(m)

    numobj = length(objectives)
    Phi = zeros(numobj, numobj)

    if numobj > 2
        # TODO:
        # The logic  here becomes difficult, as the feasible region will require
        # a dependency between the constraints
        Base.error(":EPS is thought through for > 2 objectives yet")
    end

    # Individual minimisations
    for (i, objective) in enumerate(objectives)
        @objective(m, objective.sense, objective.f)
        for (key, value) in objective.initialvalue
            set_start_value(variable_by_name(m, key), value)
        end
        optimize!(m, ignore_optimize_hook=true);
        if !(termination_status(m) in [MOI.OPTIMAL, MOI.LOCALLY_SOLVED])
            return termination_status(m)
        end

        push!(mdata.utopiavarvalues, [value(var) for var in vararr])
        push!(mdata.paretovarvalues, [value(var) for var in vararr])

        Phi[:,i] = sensevalue.(objectives)

        push!(mdata.paretofront, value.(objectives))
    end
    Fmax = maximum(Phi, dims=2)
    Fmin = minimum(Phi, dims=2) # == diag(Phi)?

    mdata.Phi = Phi
    # maintain last objective, other as constraints
    @objective(m, objectives[end].sense, objectives[end].f)

    @constraint(m, objconstr,
                  sensemap[objectives[1].sense] * objectives[1].f
                  <= Fmax[1])

    for betaval in 0.0:0.001:1.0

        set_normalized_rhs(objconstr, betaval * Fmin[1] + (1 - betaval) * Fmax[1])

        optimize!(m, ignore_optimize_hook=true);
        if !(termination_status(m) in [MOI.OPTIMAL, MOI.LOCALLY_SOLVED])
            return termination_status(m)
        end

        push!(mdata.paretofront, value.(objectives))
        push!(mdata.paretovarvalues, [value(var) for var in vararr])
    end

    return MOI.OPTIMAL
end


text:
# symlink scripts so we can call via ~/.pope/debugscript.jl
println("running Pope/build.jl")
scriptsdir = normpath(joinpath(@__DIR__, "..", "scripts"))
targetdir = expanduser("~/.pope")
if !isdir(targetdir)
    println("creating $targetdir")
    mkdir(targetdir)
end
@show scriptsdir
@show targetdir
for n in readdir(scriptsdir)
    p = joinpath(scriptsdir, n)
    t = joinpath(targetdir, n)
    if isfile(p)
        s = "adding $p $t"
        if isfile(t)
            s*=" (overwriting)"
            @show t
            rm(t)
        end
        println(s)
        open(t,"w") do f
            write(f,"julia $p \$@")
        end
        chmod(t,0o755) # all can rx, user can rwx
    end
end

ext: jl


text:
a = 0

for i in 1:100000000
	a += 1
end

println(a)


text:
# Autogenerated wrapper script for Gingko_jll for armv7l-linux-gnueabihf-cxx11
export libgingko, libgingko_cuda, libgingko_device, libgingko_dpcpp, libgingko_hip, libgingko_omp, libgingko_reference

using CompilerSupportLibraries_jll
JLLWrappers.@generate_wrapper_header("Gingko")
JLLWrappers.@declare_library_product(libgingko, "libginkgo.so.1.4.0")
JLLWrappers.@declare_library_product(libgingko_cuda, "libginkgo_cuda.so.1.4.0")
JLLWrappers.@declare_library_product(libgingko_device, "libginkgo_device.so.1.4.0")
JLLWrappers.@declare_library_product(libgingko_dpcpp, "libginkgo_dpcpp.so.1.4.0")
JLLWrappers.@declare_library_product(libgingko_hip, "libginkgo_hip.so.1.4.0")
JLLWrappers.@declare_library_product(libgingko_omp, "libginkgo_omp.so.1.4.0")
JLLWrappers.@declare_library_product(libgingko_reference, "libginkgo_reference.so.1.4.0")
function __init__()
    JLLWrappers.@generate_init_header(CompilerSupportLibraries_jll)
    JLLWrappers.@init_library_product(
        libgingko,
        "lib/libginkgo.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libgingko_cuda,
        "lib/libginkgo_cuda.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libgingko_device,
        "lib/libginkgo_device.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libgingko_dpcpp,
        "lib/libginkgo_dpcpp.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libgingko_hip,
        "lib/libginkgo_hip.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libgingko_omp,
        "lib/libginkgo_omp.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libgingko_reference,
        "lib/libginkgo_reference.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
# TODO deprecate testshuffle
_checkErrorCCWNumerics(ccwl::Union{CommonConvWrapper{F},CommonConvWrapper{Mixture{N_,F,S,T}}}, testshuffle::Bool=false)  where {N_,F<:AbstractRelativeMinimize,S,T} = nothing
_checkErrorCCWNumerics(ccwl::Union{CommonConvWrapper{F},CommonConvWrapper{Mixture{N_,F,S,T}}}, testshuffle::Bool=false)  where {N_,F<:AbstractManifoldMinimize,S,T} = nothing

function _checkErrorCCWNumerics(ccwl::Union{CommonConvWrapper{F},CommonConvWrapper{Mixture{N_,F,S,T}}},
                                testshuffle::Bool=false)  where {N_,F<:AbstractRelativeRoots,S,T}
  #
  # @info "ccwl zDim and xDim" ccwl.zDim ccwl.xDim
  if testshuffle || ccwl.partial || (!ccwl.partial && ccwl.zDim < ccwl.xDim)
    error("<:AbstractRelativeRoots factors with less measurement dimensions than variable dimensions have been discontinued, easy conversion to <:AbstractRelativeMinimize is the better option.")
  elseif !( ccwl.zDim >= ccwl.xDim && !ccwl.partial )
    error("Unresolved numeric <:AbstractRelativeRoots solve case")
  end
  nothing
end

_perturbIfNecessary(fcttype::Union{F,<:Mixture{N_,F,S,T}},
                    len::Int=1,
                    perturbation::Real=1e-10 ) where {N_,F<:AbstractRelativeMinimize,S,T} = 0

_perturbIfNecessary(fcttype::Union{F,<:Mixture{N_,F,S,T}},
                    len::Int=1,
                    perturbation::Real=1e-10 ) where {N_,F<:AbstractManifoldMinimize,S,T} = 0
#

_perturbIfNecessary(fcttype::Union{F,<:Mixture{N_,F,S,T}},
                    len::Int=1,
                    perturbation::Real=1e-10 ) where {N_,F<:AbstractRelativeRoots,S,T} = perturbation*randn(len)
#

# internal use only, and selected out from approxDeconv functions
_solveLambdaNumeric(fcttype::AbstractPrior,
                    objResX::Function,
                    residual::AbstractVector{<:Real},
                    u0::AbstractVector{<:Real},
                    islen1::Bool=false;
                    perturb::Real=1e-10 ) = u0
#

function _solveLambdaNumeric( fcttype::Union{F,<:Mixture{N_,F,S,T}},
                              objResX::Function,
                              residual::AbstractVector{<:Real},
                              u0::AbstractVector{<:Real},
                              islen1::Bool=false )  where {N_,F<:AbstractRelativeRoots,S,T}
  #

  #
  r = NLsolve.nlsolve( (res, x) -> res .= objResX(x), u0, inplace=true) #, ftol=1e-14)

  #
  return r.zero
end

function _solveLambdaNumeric( fcttype::Union{F,<:Mixture{N_,F,S,T}},
                              objResX::Function,
                              residual::AbstractVector{<:Real},
                              u0::AbstractVector{<:Real},
                              islen1::Bool=false  )  where {N_,F<:AbstractRelativeMinimize,S,T}
                              # retries::Int=3 )
  #
  # wrt #467 allow residual to be standardize for Roots and Minimize and Parametric cases.
  r = if islen1
    Optim.optimize((x) -> (residual .= objResX(x); sum(residual.^2)), u0, Optim.BFGS() )
  else
    Optim.optimize((x) -> (residual .= objResX(x); sum(residual.^2)), u0)
  end

  #
  return r.minimizer
end

function _solveLambdaNumeric( fcttype::Union{F,<:Mixture{N_,F,S,T}},
                              objResX::Function,
                              residual::AbstractVector{<:Real},
                              u0,#::AbstractVector{<:Real},
                              variableType::InferenceVariable,
                              islen1::Bool=false)  where {N_,F<:AbstractManifoldMinimize,S,T}
  #
  M = getManifold(variableType)#fcttype.M
  # the variable is a manifold point, we are working on the tangent plane in optim for now.
  #
  #TODO this is not general to all manifolds, should work for lie groups.
  # ϵ = identity_element(M, u0)
  ϵ = getPointIdentity(variableType)
  # X0c = get_coordinates(M, u0, log(M, ϵ, u0), DefaultOrthogonalBasis())
  X0c = vee(M, u0, log(M, ϵ, u0))

  # objResX(p) returns tangent vector at p, X=log(M, p, ...)
  # norm(M, p, X) == distance(M, p, X)
  #TODO fix closure for performance
  fM = getManifold(fcttype)
  function cost(p, X, Xc)
    hat!(M, X, ϵ, Xc)
    exp!(M, p, ϵ, X)
    # X = objResX(p)
    # return norm(fM, p, X)^2 #TODO the manifold of p and X are not always the same
    #options getPointIdentity or leave it to factor
    residual = objResX(p)
    return sum(residual.^2)
  end

  # # separate statements to try improve type-stability
  # r = if islen1
  #   Optim.optimize(cost, X0c, Optim.BFGS())
  # else
  #   Optim.optimize(cost, X0c, Optim.NelderMead())
  # end
  alg = islen1 ? Optim.BFGS() : Optim.NelderMead()
  X0 = hat(M, ϵ, X0c)
  p0 = exp(M, ϵ, X0)
  r = Optim.optimize(Xc->cost(p0, X0, Xc), X0c, alg)
  if !Optim.converged(r)
    @debug "Optim did not converge:" r
  end
  return exp(M, ϵ, hat(M, ϵ, r.minimizer))

end

#TODO Consolidate with _solveLambdaNumeric, see #1374
function _solveLambdaNumericMeas( fcttype::Union{F,<:Mixture{N_,F,S,T}},
                                  objResX::Function,
                                  residual::AbstractVector{<:Real},
                                  u0,#::AbstractVector{<:Real},
                                  variableType::InferenceVariable,
                                  islen1::Bool=false)  where {N_,F<:AbstractManifoldMinimize,S,T}
  #
  # Assume measurement is on the tangent
  M = getManifold(variableType)#fcttype.M
  # the variable is a manifold point, we are working on the tangent plane in optim for now.
  ϵ = getPointIdentity(variableType)
  X0c = vee(M, ϵ, u0)

  function cost(X, Xc)
    hat!(M, X, ϵ, Xc)
    residual = objResX(X)
    return sum(residual.^2)
  end

  alg = islen1 ? Optim.BFGS() : Optim.NelderMead()
  X0 = hat(M, ϵ, X0c)
  r = Optim.optimize(Xc->cost(X0, Xc), X0c, alg)
  if !Optim.converged(r)
    @debug "Optim did not converge:" r
  end

  return hat(M, ϵ, r.minimizer)

end

## ================================================================================================
## Heavy dispatch for all AbstractFactor / Mixture cases below
## ================================================================================================

# internal function to dispatch view on either vector or matrix, rows are dims and samples are columns
_getindextuple(tup::Tuple, ind1::Int) = [getindex(t, ind1) for t in tup]

# TODO, likely a shortlived function, and should be replaced with ccw.hypoParams::Tuple(hypo1, hypo2,...), made at construction and allows direct hypo lookup
# DevNotes, also see new `hyporecipe` approach (towards consolidation CCW CPT FMd CF...)
function _view(nt::NamedTuple, idxs::AbstractVector{<:Integer})
  varParams = tuple([nt[i] for i in idxs]...)
  tup = tuple(varParams...)
  nms = keys(nt)[idxs]
  return NamedTuple{nms,typeof(tup)}(tup)
  # tuple([nt[i] for i in idxs]...)
end

function _buildCalcFactorMixture( ccwl::CommonConvWrapper,
                                  _fmd_,
                                  smpid,
                                  measurement_,
                                  varParams )
  #
  # FIXME, make thread safe (cache)
  CalcFactor( ccwl.usrfnc!, _fmd_, smpid,
              length(measurement_), measurement_, varParams, true, ccwl.dummyCache)
end

function _buildCalcFactorMixture( ccwl::CommonConvWrapper{Mixture{N_,F,S,T}},
                                  _fmd_,
                                  smpid,
                                  measurement_,
                                  varParams ) where {N_,F <: AbstractFactor,S,T}
  #
  # just a passthrough similar to pre-v0.20
  # FIXME, make thread safe (cache)
  CalcFactor( ccwl.usrfnc!.mechanics, _fmd_, smpid,
              length(measurement_), measurement_, varParams, true, ccwl.dummyCache)
end

"""
    $SIGNATURES
Internal function to build lambda pre-objective function for finding factor residuals.

Notes
- Unless passed in as separate arguments, this assumes already valid in `cpt_`:
  - `cpt_.p`
  - `cpt_.activehypo`
  - `cpt_.factormetadata`
  - `ccwl.params`
  - `ccwl.measurement`

DevNotes
- TODO refactor relationship and common fields between (CCW, FMd, CPT, CalcFactor)
"""
function _buildCalcFactorLambdaSample(ccwl::CommonConvWrapper,
                                      smpid::Integer,
                                      cpt_::ConvPerThread = ccwl.cpt[Threads.threadid()],
                                      target = view(ccwl.params[ccwl.varidx][smpid], ccwl.partialDims),
                                      measurement_ = ccwl.measurement,
                                      fmd_::FactorMetadata = cpt_.factormetadata;
                                      _slack = nothing  )
  #

  # build a view to the decision variable memory
  varValsHypo = _view(ccwl.params, cpt_.activehypo)
  # tup = tuple(varParams...)
  # nms = keys(ccwl.params)[cpt_.activehypo]
  # varValsHypo = NamedTuple{nms,typeof(tup)}(tup)

  # prepare fmd according to hypo selection
  # FIXME must refactor (memory waste) and consolidate with CCW CPT FMd CF
  # FIXME move up out of smpid loop and only update bare minimal fields
  _fmd_ = FactorMetadata( view(fmd_.fullvariables, cpt_.activehypo),
                          view(fmd_.variablelist, cpt_.activehypo),
                          varValsHypo, #varParams, # view(fmd_.arrRef, cpt_.activehypo),
                          fmd_.solvefor,
                          fmd_.cachedata  )
  #
  # get the operational CalcFactor object
  cf = _buildCalcFactorMixture(ccwl, _fmd_, smpid, measurement_, varValsHypo)
  # new dev work on CalcFactor
  # cf = CalcFactor(ccwl.usrfnc!, _fmd_, smpid,
  #                 length(measurement_), measurement_, varValsHypo)
  #

  # reset the residual vector
  fill!(cpt_.res, 0.0) # Roots->xDim | Minimize->zDim

  # build static lambda
  unrollHypo! = if _slack === nothing
    () -> cf( measurement_[smpid], map(vvh->getindex(vvh, smpid), varValsHypo)... )
  else
    # slack is used to shift the residual away from the natural "zero" tension position of a factor,
    # this is useful when calculating factor gradients at a variety of param locations resulting in "non-zero slack" of the residual.
    # see `IIF.calcFactorResidualTemporary`
    # NOTE this minus operation assumes _slack is either coordinate or tangent vector element (not a manifold or group element)
    () -> cf( measurement_[smpid], map(vvh->getindex(vvh, smpid), varValsHypo)... ) .- _slack
  end

  return unrollHypo!, target
end

"""
    $(SIGNATURES)

Solve free variable x by root finding residual function `fgr.usrfnc(res, x)`.  This is the
penultimate step before calling numerical operations to move actual estimates, which is
done by an internally created lambda function.

ccw.X must be set to memory ref the param[varidx] being solved, at creation of ccw

Notes
- Assumes `cpt_.p` is already set to desired X decision variable dimensions and size.
- Assumes only `ccw.particleidx` will be solved for
- small random (off-manifold) perturbation used to prevent trivial solver cases, div by 0 etc.
  - perturb is necessary for NLsolve cases, and smaller than 1e-10 will result in test failure
- Also incorporates the active hypo lookup

DevNotes
- TODO testshuffle is now obsolete, should be removed
- TODO perhaps consolidate perturbation with inflation or nullhypo
"""
function _solveCCWNumeric!( ccwl::Union{CommonConvWrapper{F},
                                        CommonConvWrapper{Mixture{N_,F,S,T}}};
                            perturb::Real=1e-10,
                            # testshuffle::Bool=false,
                            _slack=nothing  ) where {N_,F<:AbstractRelative,S,T}
  #

  #
  thrid = Threads.threadid()
  cpt_ = ccwl.cpt[thrid]

  smpid = cpt_.particleidx
  # cannot Nelder-Mead on 1dim, partial can be 1dim or more but being conservative.
  islen1 = length(ccwl.partialDims) == 1 || ccwl.partial
  # islen1 = length(cpt_.X[:, smpid]) == 1 || ccwl.partial

  # build the pre-objective function for this sample's hypothesis selection
  unrollHypo!, target = _buildCalcFactorLambdaSample(ccwl, smpid, cpt_, _slack=_slack )

  # broadcast updates original view memory location
  ## using CalcFactor legacy path inside (::CalcFactor)
  _hypoObj = (x) -> (target.=x; unrollHypo!())

  # TODO small off-manifold perturbation is a numerical workaround only, make on-manifold requires RoME.jl #244
  # use all element dimensions : ==> 1:ccwl.xDim
  target .+= _perturbIfNecessary(getFactorType(ccwl), length(target), perturb)

  # do the parameter search over defined decision variables using Minimization
  retval = _solveLambdaNumeric(getFactorType(ccwl), _hypoObj, cpt_.res, cpt_.X[smpid][ccwl.partialDims], islen1 )

  # Check for NaNs
  if sum(isnan.(retval)) != 0
    @error "$(ccwl.usrfnc!), ccw.thrid_=$(thrid), got NaN, smpid = $(smpid), r=$(retval)\n"
    return nothing
  end

  # insert result back at the correct variable element location
  cpt_.X[smpid][ccwl.partialDims] .= retval

  nothing
end
# brainstorming
# should only be calling a new arg list according to activehypo at start of particle
# Try calling an existing lambda
# sensitive to which hypo of course , see #1024
# need to shuffle content inside .cpt.fmd as well as .params accordingly
#

function _solveCCWNumeric!( ccwl::Union{CommonConvWrapper{F},
                                        CommonConvWrapper{Mixture{N_,F,S,T}}};
                            perturb::Real=1e-10,
                            # testshuffle::Bool=false,
                            _slack=nothing  ) where {N_,F<:AbstractManifoldMinimize,S,T}
  #
  #   # FIXME, move this check higher and out of smpid loop
  # _checkErrorCCWNumerics(ccwl, testshuffle)

  #
  thrid = Threads.threadid()
  cpt_ = ccwl.cpt[thrid]

  smpid = cpt_.particleidx
  # cannot Nelder-Mead on 1dim, partial can be 1dim or more but being conservative.
  islen1 = length(ccwl.partialDims) == 1 || ccwl.partial

  # build the pre-objective function for this sample's hypothesis selection
  unrollHypo!, target = _buildCalcFactorLambdaSample(ccwl, smpid, cpt_, view(ccwl.params[ccwl.varidx], smpid), _slack=_slack)

  # broadcast updates original view memory location
  ## using CalcFactor legacy path inside (::CalcFactor)
  # _hypoObj = (x) -> (target.=x; unrollHypo!())
  function _hypoObj(x)
    target[] .= x
    return unrollHypo!()
  end

  # TODO small off-manifold perturbation is a numerical workaround only, make on-manifold requires RoME.jl #244
  # use all element dimensions : ==> 1:ccwl.xDim
  # F <: AbstractRelativeRoots && (target .+= _perturbIfNecessary(getFactorType(ccwl), length(target), perturb))

  # do the parameter search over defined decision variables using Minimization
  # retval = _solveLambdaNumeric(getFactorType(ccwl), _hypoObj, cpt_.res, cpt_.X[smpid][cpt_.p], islen1 )
  sfidx = ccwl.varidx
  retval = _solveLambdaNumeric(getFactorType(ccwl), _hypoObj, cpt_.res, cpt_.X[smpid], ccwl.vartypes[sfidx](), islen1)

  # Check for NaNs
  # FIXME isnan for manifold ProductRepr
  # if sum(isnan.(retval)) != 0
    # @error "$(ccwl.usrfnc!), ccw.thrid_=$(thrid), got NaN, smpid = $(smpid), r=$(retval)\n"
    # return nothing
  # end

  # FIXME insert result back at the correct variable element location
  cpt_.X[smpid] .= retval

  nothing
end

#


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md

abstract AbstractField

abstract Discrete <: AbstractField
abstract Continuous <: AbstractField
abstract Constant <: AbstractField
abstract Variable <: AbstractField
abstract TimeVariant <: AbstractField
abstract TimeInvariant <: AbstractField

type Field{A<:Union{Discrete,Continuous}, B<:Union{Constant,Variable}, C<:Union{TimeVariant,TimeInvariant}}
    data
end

typealias FieldSet Dict{AbstractString, Field}

### Basic data structure for discrete field

type Increment{T}
    time :: Float64
    data :: T
end

function Base.convert{T}(::Type{Increment{T}}, data::Pair{Float64,T})
    return Increment{T}(data[1], data[2])
end

function Base.convert{T}(::Type{Increment{Vector{Vector{T}}}}, data::Pair{Float64, Matrix{T}})
    time = data[1]
    content = data[2]
    return Increment(time, Vector{T}[content[:,i] for i=1:size(content,2)])
end

function Base.getindex{T}(increment::Increment{Vector{T}}, i::Int64)
    return increment.data[i]
end

function Base.(:*)(d, increment::Increment)
    return d*increment.data
end

### Basic data structure for continuous field

type Basis
    basis :: Function
    dbasis :: Function
end

function call(basis::Basis, xi::Vector)
    basis.basis(xi)
end

function call(basis::Basis, xi::Vector, ::Type{Val{:grad}})
    basis.dbasis(xi)
end

### Different field combinations and other typealiases

typealias DCTI Field{Discrete,   Constant, TimeInvariant}
typealias DVTI Field{Discrete,   Variable, TimeInvariant}
typealias DCTV Field{Discrete,   Constant, TimeVariant}
typealias DVTV Field{Discrete,   Variable, TimeVariant}
typealias CCTI Field{Continuous, Constant, TimeInvariant}
typealias CVTI Field{Continuous, Variable, TimeInvariant} # can be used to interpolate in spatial dimension
typealias CCTV Field{Continuous, Constant, TimeVariant} # can be used to interpolate in time
typealias CVTV Field{Continuous, Variable, TimeVariant}

typealias ScalarIncrement{T} Increment{T}
typealias VectorIncrement{T} Increment{Vector{T}}
typealias TensorIncrement{T} Increment{Matrix{T}}

typealias DiscreteField      Union{DCTI, DVTI, DCTV, DVTV}
typealias ContinuousField    Union{CCTI, CVTI, CCTV, CVTV}
typealias ConstantField      Union{DCTI, DCTV, CCTI, CCTV}
typealias VariableField      Union{DVTI, DVTV, CVTI, CVTV}
typealias TimeInvariantField Union{DCTI, DVTI, CCTI, CVTI}
typealias TimeVariantField   Union{DCTV, DVTV, CCTV, CVTV}

### Convenient functions to create fields

#function Base.convert(::Type{Field}, data)
#    return Field(data)
#end

function Field(data)
    return DCTI(data)
end

function Field(data::Vector)
    return DVTI(data)
end

function Field{T}(data::Pair{Float64, T}...)
    return DCTV([Increment{T}(d[1], d[2]) for d in data])
end
#=
function Field{T}(data::Pair{Float64, Vector{T}}...)
    return DVTV([Increment{Vector{T}}(d[1], d[2]) for d in data])
end

function Field{T}(data::Pair{Float64, Dict{Int64, T}}...)
    return DVTV([Increment{Dict{Int64, T}}(d[1], d[2]) for d in data])
end
=#

function Field{T<:Union{Vector, Dict}}(data::Pair{Float64, T}...)
    return DVTV([Increment{T}(d[1], d[2]) for d in data])
end

function Field(data::Dict)
    return DVTI(data)
end

function convert{T}(::Type{DCTV}, data::Pair{Real, Vector{T}}...)
    return DCTV([Increment{Vector{T}}(d[1], d[2]) for d in data])
end

""" Create new discrete, constant, time variant field.

Examples
--------
julia> t0 = 0.0; t1=1.0; y0 = 0.0; y1 = 1.0
julia> f = DCTV(t0 => y0, t1 => y1)

"""
function convert{T,v<:Real}(::Type{DCTV}, data::Pair{v, T}...)
    return DCTV([Increment(d[1],d[2]) for d in data])
end
#function Base.convert(::Type{DCTV}, data::Pair{Real, Any}...)
#    return DCTV([Increment{Vector}(d[1], d[2]) for d in data])
#end

function Field(func::Function)
    if method_exists(func, Tuple{})
        return CCTI(func)
    elseif method_exists(func, Tuple{Float64})
        return CCTV(func)
    elseif method_exists(func, Tuple{Vector})
        return CVTI(func)
    elseif method_exists(func, Tuple{Vector, Number})
        return CVTV(func)
    else
        error("no proper definition found for function: check methods.")
    end
end

function CVTI(basis::Function, dbasis::Function)
    return CVTI(Basis(basis, dbasis))
end

function Field(basis::Function, dbasis::Function)
    return CVTI(basis, dbasis)
end

### Accessing and manipulating discrete fields

function getindex(field::DVTV, i::Int64)
    return field.data[i]
end

function push!(field::DCTV, data::Pair)
    push!(field.data, data)
end

function push!(field::DVTV, data::Pair)
    push!(field.data, data)
end

function getindex(field::DVTV, i::Int64)
    return field.data[i]
end

function getindex(field::DVTI, i::Int64)
    return field.data[i]
end

function getindex(field::DCTV, i::Int64)
    return field.data[i]
end

function getindex(field::Field, i::Int64)
    return field.data[i]
end

function length(field::DVTI)
    return length(field.data)
end

function length(field::DCTI)
    return 1
end

function length(field::DVTV)
    return length(field.data)
end

function length(field::DCTV)
    return length(field.data)
end

function first(field::Union{DCTV, DVTV})
    return field[1]
end

function isapprox(f1::DCTI, f2::DCTI)
    isapprox(f1.data, f2.data)
end

for op = (:+, :*, :/, :-)
    @eval ($op)(increment::Increment, field::DCTI) = ($op)(increment.data, field.data)
    @eval ($op)(field::DCTI, increment::Increment) = ($op)(increment.data, field.data)
    @eval ($op)(field1::DCTI, field2::DCTI) = ($op)(field1.data, field2.data)
    @eval ($op)(field::DCTI, k::Number) = ($op)(field.data, k)
    @eval ($op)(k::Number, field::DCTI) = ($op)(field.data, k)
end

function Base.(:+)(f1::DVTI, f2::DVTI)
    return DVTI(f1.data + f2.data)
end

function Base.(:-)(f1::DVTI, f2::DVTI)
    return DVTI(f1.data - f2.data)
end

function Base.(:*){T<:Real}(c::T, field::DVTI)
    return DVTI(c*field.data)
end

function Base.(:*)(N::Matrix, f::DCTI)
    return f.data*N'
end

#
#   Multiply DVTI field with another vector T. Vector length
#   must match to the field length and this can be used mainly
#   for interpolation purposes, i.e., u = ∑ Nᵢuᵢ
#
function Base.(:*)(T::Vector, f::DVTI)
    @assert length(T) == length(f)
    return sum([T[i]*f[i] for i=1:length(f)])
end

function vec(field::DVTI)
    return [field.data...;]
end

function vec(field::DCTV)
    error("trying to vectorize $field does not make sense")
end

function endof(field::Field)
    return endof(field.data)
end

#function Base.similar{T}(field::DVTI, data::Vector{T})
#    return Increment(reshape(data, round(Int, length(data)/length(increment)), length(increment)))
#end

function similar{T}(field::DVTI, data::Vector{T})
    n = length(field.data)
    data = reshape(data, round(Int, length(data)/n), n)
    newdata = Vector[data[:,i] for i=1:n]
    return typeof(field)(newdata)
end

function start(::DVTI)
    return 1
end

function next(f::DVTI, state)
    return f.data[state], state+1
end

function done(f::DVTI, s)
    return s > length(f.data)
end

""" Update time-dependent fields with new values.

Examples
--------

julia> f = Field(0.0 => 1.0)
julia> update!(f, 1.0 => 2.0)

Now field has two (time, value) pairs: (0.0, 1.0) and (1.0, 2.0)

Notes
-----
Time vector is assumed to be ordered t_i-1 < t_i < t_i+1. If updating
field with already existing time the old value is replaced with new one.

"""
function update!{T}(field::Union{DCTV, DVTV}, val::Pair{Float64, T})
    time, data = val
    if isapprox(last(field).time, time)
        last(field).data = data
    else
        push!(field.data, Increment(val...))
    end
end

function update!{T}(field::Union{DCTI, DVTI}, val::T)
    field.data = val
end

### Accessing continuous fields

function call(field::CVTI, xi::Vector)
    return field.data(xi)
end

function call(field::CVTV, xi, time::Float64)
    return field.data(xi, time)
end

function call(field::CVTI, xi::Vector, ::Type{Val{:Grad}})
    return field.data(xi, Val{:Grad})
end

function call(field::CCTV, time::Float64)
    return field.data(time)
end

function convert(::Type{Basis}, field::CVTI)
    return field.data
end

### Interpolation

""" Interpolate time-invariant field in time direction. """
function call(field::DVTI, time::Float64)
    return field
end
function call(field::DCTI, time::Float64)
    return field
end
function call(field::CVTI, time::Float64)
    return field.data()
end
function call(field::CCTI, time::Float64)
    return field.data()
end

""" Interpolate constant time-variant field in time direction. """
function call(field::DCTV, time::Real)
    time < first(field).time && return DCTI(first(field).data)
    time > last(field).time && return DCTI(last(field).data)
    for i=reverse(1:length(field))
        isapprox(field[i].time, time) && return DCTI(field[i].data)
    end
    for i=reverse(2:length(field))
        t0 = field[i-1].time
        t1 = field[i].time
        if t0 < time < t1
            y0 = field[i-1].data
            y1 = field[i].data
            dt = t1-t0
            new_data = y0*(1-(time-t0)/dt) + y1*(1-(t1-time)/dt)
            return DCTI(new_data)
        end
    end
    error("interpolate DCTV: unknown failure when interpolating $(field.data) for time $time")
end

function call(field::DVTV, time::Float64)
    time < first(field).time && return DVTI(first(field).data)
    time > last(field).time && return DVTI(last(field).data)
    for i=reverse(1:length(field))
        isapprox(field[i].time, time) && return DVTI(field[i].data)
    end
    for i=reverse(2:length(field))
        t0 = field[i-1].time
        t1 = field[i].time
        if t0 < time < t1
            y0 = field[i-1].data
            y1 = field[i].data
            dt = t1-t0
            new_data = y0*(1-(time-t0)/dt) + y1*(1-(t1-time)/dt)
            return DVTI(new_data)
        end
    end
    error("interpolate DVTV: unknown failure when interpolating $(field.data) for time $time")
end

""" Interpolate constant field in spatial dimension. """
function call(basis::CVTI, field::DCTI, xi::Vector)
    return field.data
end

""" Interpolate variable field in spatial dimension. """
function call(basis::CVTI, values::DVTI, xi::Vector)
    N = basis(xi)
    return sum([N[i]*values[i] for i=1:length(N)])
end

function call(basis::CVTI, geometry::DVTI, xi::Vector, ::Type{Val{:grad}})
    dbasis = basis(xi, Val{:grad})
#    J = sum([dbasis[:,i]*geometry[i]' for i=1:length(geometry)])
    J = sum([kron(dbasis[:,i], geometry[i]') for i=1:length(geometry)])
    invJ = isa(J, Vector) ? inv(J[1]) : inv(J)
    grad = invJ * dbasis
    return grad
end

function call(basis::CVTI, geometry::DVTI, values::DVTI, xi::Vector, ::Type{Val{:grad}})
    grad = basis(geometry, xi, Val{:grad})
#    gradf = sum([grad[:,i]*values[i]' for i=1:length(geometry)])'
    gradf = sum([kron(grad[:,i], values[i]') for i=1:length(values)])'
    return length(gradf) == 1 ? gradf[1] : gradf
end

function call(basis::CVTI, xi::Vector, time::Number)
    basis(xi)
end

function Base.(:*)(grad::Matrix, field::DVTI)
    return sum([kron(grad[:,i], field[i]') for i=1:length(field)])'
end

function DVTV(data::Pair{Float64, Vector}...)
    return DVTV([Increment(d[1], d[2]) for d in data])
end

function start(f::DVTV)
    return start(f.data)
end

function next(f::DVTV, state)
    return next(f.data, state)
end

function done(f::DVTV, state)
    return done(f.data, state)
end

""" Return time vector from time variable field. """
function keys(field::DVTV)
    return Float64[increment.time for increment in field]
end

function setindex!(field::Field, val, idx::Int64)
    field.data[idx] = val
end


text:
### Struct variables
struct Path
    packagepath::String
    fullpath::String
    dataname::String
    extension::String
end

struct PowerSystem
    bus::Array{Float64,2}
    branch::Array{Float64,2}
    generator::Array{Float64,2}
    generatorcost::Array{Float64,2}
    basePower::Float64
end

struct PowerSystemNum
    Nbus::Int64
    Nbranch::Int64
    Ngen::Int64
    Ncost::Int64
end

struct FlowSettings
    algorithm::String
    solve::String
    main::Bool
    flow::Bool
    generation::Bool
    save::String
    maxIter::Int64
    stopping::Float64
    reactive::Array{Bool,1}
end

struct Measurements
    pmuVoltage::Array{Float64,2}
    pmuCurrent::Array{Float64,2}
    legacyFlow::Array{Float64,2}
    legacyCurrent::Array{Float64,2}
    legacyInjection::Array{Float64,2}
    legacyVoltage::Array{Float64,2}
end

struct MeasurementsNum
    pmuNv::Int64
    pmuNc::Int64
    legacyNf::Int64
    legacyNc::Int64
    legacyNi::Int64
    legacyNv::Int64
end

struct GeneratorSettings
    runflow::Int64
    set::Dict{Symbol,Array{Float64,1}}
    variance::Dict{Symbol,Array{Float64,1}}
    save::String
end

struct EstimationSettings
    algorithm::String
    main::Bool
    flow::Bool
    estimate::Bool
    error::Bool
    maxIter::Int64
    stopping::Float64
    start::Dict{Symbol,Float64}
    bad::Dict{Symbol,Float64}
    lav::Dict{Symbol,Float64}
    observe::Dict{Symbol,Float64}
    covariance::Bool
    solve::String
    save::String
    saveextension::String
end

struct OptimalFlowSettings
    algorithm::String
    main::Bool
    flow::Bool
    generation::Bool
    save::String
end

### Load path
function loadpath(args)
    pathtoJuliaGrid = Base.find_package("JuliaGrid")
    if pathtoJuliaGrid == nothing
        throw(ErrorException("JuliaGrid not found in install packages"))
    end
    packagepath = abspath(joinpath(dirname(pathtoJuliaGrid), ".."))
    extension = ".h5"; path = ""; dataname = ""; fullpath = ""
    for i = 1:length(args)
        try
            extension = string(match(r"\.[A-Za-z0-9]+$", args[i]).match)
        catch
            extension = ""
        end
        if extension == ".h5" || extension == ".xlsx" || extension == ".m"
            fullpath = args[i]
            path = dirname(args[i])
            dataname = basename(args[i])
            break
        end
    end

    if isempty(extension)
        throw(ErrorException("the input DATA extension is not found"))
    elseif extension != ".h5" && extension != ".xlsx" && extension != ".m"
        throw(DomainError(extension, "the input DATA extension is not supported"))
    end

    if path == ""
        path = joinpath(packagepath, "src/data/")
        fullpath = joinpath(packagepath, "src/data/", dataname)
    end

    if dataname in cd(readdir, path)
        println("The input power system: $dataname")
    else
        throw(DomainError(dataname, "the input DATA is not found"))
    end

    return Path(packagepath, fullpath, dataname, extension)
end

### Load power system data
function loadsystem(path)
    if path.extension == ".h5"
        fid = h5open(path.fullpath, "r")
            bus::Array{Float64,2} = h5read(path.fullpath, "/bus")
            branch::Array{Float64,2} = h5read(path.fullpath, "/branch")
            Nbus = datastruct(bus, 13; var = "bus")
            Nbranch = datastruct(branch, 14; var = "branch")
        if exists(fid, "generator")
            generator::Array{Float64,2} = h5read(path.fullpath, "/generator")
            Ngen = datastruct(generator, 21; var = "generator")
        else
            generator = zeros(1, 21)
            Ngen = 0
        end
        if exists(fid, "generatorcost")
            generatorcost::Array{Float64,2} = h5read(path.fullpath, "/generatorcost")
            Ncost = datastruct(generatorcost, 5; var = "generatorcost")
        else
            generatorcost = zeros(1, 5)
            Ncost = 0
        end
        if exists(fid, "basePower")
            basePower::Float64 = h5read(path.fullpath, "/basePower")[1]
        else
            basePower = 100.0
            println("The variable basePower not found. The algorithm proceeds with default value: 100 MVA.")
        end
        if exists(fid, "info")
            info::Array{String,2} = h5read(path.fullpath, "/info")
        else
            info = Array{String}(undef, 0, 1)
        end
        close(fid)
    end

    if path.extension == ".xlsx"
        xf = XLSX.openxlsx(path.fullpath, mode = "r")
        if "bus" in XLSX.sheetnames(xf)
            start = startxlsx(xf["bus"])
            bus = xf["bus"][:][start:end, :]
            Nbus = datastruct(bus, 13; var = "bus")
        else
            throw(ErrorException("error opening sheet bus"))
        end
        if "branch" in XLSX.sheetnames(xf)
            start = startxlsx(xf["branch"])
            branch = xf["branch"][:][start:end, :]
            Nbranch = datastruct(branch, 14; var = "branch")
        else
            throw(ErrorException("error opening sheet branch"))
        end
        if "generator" in XLSX.sheetnames(xf)
            start = startxlsx(xf["generator"])
            generator = xf["generator"][:][start:end, :]
            Ngen = datastruct(generator, 21; var = "generator")
        else
            generator = zeros(1, 21)
            Ngen = 0
        end
        if "generatorcost" in XLSX.sheetnames(xf)
            start = startxlsx(xf["generatorcost"])
            generatorcost = xf["generatorcost"][:][start:end, :]
            Ncost = datastruct(generatorcost, 5; var = "generatorcost")
        else
            generatorcost = zeros(1, 5)
            Ncost = 0
        end
        if "basePower" in XLSX.sheetnames(xf)
            start = startxlsx(xf["basePower"])
            basePower = Float64(xf["basePower"][:][start])
        else
            basePower = 100.0
            println("The variable basePower not found. The algorithm proceeds with default value: 100 MVA.")
        end
        if "info" in XLSX.sheetnames(xf)
            info = convert(Array{String,2}, coalesce.(xf["info"][:], ""))
        else
            info = Array{String}(undef, 0, 1)
        end
        # close(xf)
    end

    if path.extension == ".m"
        busline = String[]; busflag = false
        branchline = String[]; branchflag = false
        genline = String[]; genflag = false
        gencostline = String[]; gencostflag = false

        datafile = open(path.fullpath, "r")
        lines = readlines(datafile)
        basePower = 0.0

        for (i, line) in enumerate(lines)
            if occursin("mpc.branch", line) && occursin("[", line)
                branchflag = true
            elseif occursin("mpc.bus", line) && occursin("[", line)
                busflag = true
            elseif occursin("mpc.gen", line) && occursin("[", line) && !occursin("mpc.gencost", line)
                genflag = true
            elseif occursin("mpc.gencost", line) && occursin("[", line)
                gencostflag = true
            elseif occursin("mpc.baseMVA", line)
                line = split(line, "=")[end]
                line = split(line, ";")[1]
                basePower = parse(Float64, line)
            end

            if branchflag
                branchflag, branchline = parse_line(i, line, branchflag, branchline)
            elseif busflag
                busflag, busline = parse_line(i, line, busflag, busline)
            elseif genflag
                genflag, genline = parse_line(i, line, genflag, genline)
            elseif gencostflag
                gencostflag, gencostline = parse_line(i, line, gencostflag, gencostline)
            end
        end

        if basePower == 0
            basePower = 100.0
            println("The variable basePower not found. The algorithm proceeds with default value: 100 MVA.")
        end
        if isempty(busline)
            throw(ErrorException("error opening variable mpc.bus"))
        end
        if isempty(branchline)
            throw(ErrorException("error opening variable mpc.branch"))
        end

        Ncol = length(parse.(Float64, split(busline[1])))
        Nbus = length(busline)
        bus = zeros(Nbus, Ncol)
        for (k, line) in enumerate(busline)
            bus[k, :] = parse.(Float64, split(line))
        end

        Ncol = length(parse.(Float64, split(branchline[1]))) + 1
        Nbranch = length(branchline)
        branch = zeros(Nbranch, Ncol)
        for (k, line) in enumerate(branchline)
            branch[k, :] = [k; parse.(Float64, split(line))]
        end

        if !isempty(genline)
            Ncol = length(parse.(Float64, split(genline[1])))
            Ngen = length(genline)
            generator = zeros(Ngen, Ncol)
            for (k, line) in enumerate(genline)
                generator[k, :] = parse.(Float64, split(line))
            end
        else
            generator = zeros(1, 21)
            Ngen = 0
        end

        if !isempty(gencostline)
            Ncol = length(parse.(Float64, split(gencostline[1])))
            Ncost = length(gencostline)
            generatorcost = zeros(Ncost, Ncol)
            for (k, line) in enumerate(gencostline)
                generatorcost[k, :] = parse.(Float64, split(line))
            end
        else
            generatorcost = zeros(1, 5)
            Ncost = 0
        end

        info = Array{String}(undef, 0, 1)
    end

    info = infogrid(bus, branch, generator, info, path.dataname, Nbranch, Nbus, Ngen)

    return PowerSystem(bus, branch, generator, generatorcost, basePower), PowerSystemNum(Nbus, Nbranch, Ngen, Ncost), info
end

### Power flow settings
@inbounds function pfsettings(args, max, stop, react, solve, save, system, num)
    algorithm = "false"
    main = false
    flow = false
    generation = false
    reactive = [false; true; false]

    for i = 1:length(args)
        if args[i] in ["nr", "gs", "fnrxb", "fnrbx", "dc"]
            algorithm = args[i]
        end
        if args[i] == "main"
            main = true
        end
        if args[i] == "flow"
            flow = true
        end
        if args[i] == "generation" && num.Ngen != 0
            generation = true
        end
    end

    if react == 1 && num.Ngen != 0
        reactive = [true; true; false]
    end

    if algorithm == "false"
        algorithm = "nr"
        println("Invalid power flow METHOD key. The algorithm proceeds with the AC power flow.")
    end

    if !isempty(save)
        path = dirname(save)
        data = basename(save)
        if isempty(data)
            dataname = join(replace(split(system.dataname, ""), "."=>""))
            data = string(dataname, "_results", system.extension)
        end
        save = joinpath(path, data)
    end

    return FlowSettings(algorithm, solve, main, flow, generation, save, max, stop, reactive)
end

### Load measurement data
@inbounds function loadmeasurement(path, system, numsys; pmuvar = [], legvar = [], runflow = 0)
    if runflow == 1
        pmuNv = numsys.Nbus
        pmuNc = 2 * numsys.Nbranch
        legacyNf = 2 * numsys.Nbranch
        legacyNc = 2 * numsys.Nbranch
        legacyNi = numsys.Nbus
        legacyNv = numsys.Nbus

        pmuVoltage = fill(0.0, pmuNv, 9)
        pmuCurrent = fill(0.0, pmuNc, 11)
        legacyFlow = fill(0.0, legacyNf, 11)
        legacyCurrent = fill(0.0, legacyNc, 7)
        legacyInjection = fill(0.0, legacyNi, 9)
        legacyVoltage = fill(0.0, legacyNv, 5)
    else
        flag = isempty(pmuvar) || isempty(legvar)

        if path.extension == ".h5"
            fid = h5open(path.fullpath, "r")
            if exists(fid, "pmuVoltage")
                pmuVoltage::Array{Float64,2} = h5read(path.fullpath, "/pmuVoltage")
                if flag
                    pmuNv = datastruct(pmuVoltage, 7; var = "pmuVoltage")
                else
                    pmuNv = datastruct(pmuVoltage, 9; var = "pmuVoltage")
                end
            else
                pmuVoltage = zeros(1, 9); pmuNv = 0
            end
            if exists(fid, "pmuCurrent")
                pmuCurrent::Array{Float64,2} = h5read(path.fullpath, "/pmuCurrent")
                if flag
                    pmuNc = datastruct(pmuCurrent, 9; var = "pmuCurrent")
                else
                    pmuNc = datastruct(pmuCurrent, 11; var = "pmuCurrent")
                end
            else
                pmuCurrent = zeros(1, 11); pmuNc = 0
            end
            if exists(fid, "legacyFlow")
                legacyFlow::Array{Float64,2} = h5read(path.fullpath, "/legacyFlow")
                if flag
                    legacyNf = datastruct(legacyFlow, 9; var = "legacyFlow")
                else
                    legacyNf = datastruct(legacyFlow, 11; var = "legacyFlow")
                end
            else
                legacyFlow = zeros(1, 11); legacyNf = 0
            end
            if exists(fid, "legacyCurrent")
                legacyCurrent::Array{Float64,2} = h5read(path.fullpath, "/legacyCurrent")
                if flag
                    legacyNc = datastruct(legacyCurrent, 6; var = "legacyCurrent")
                else
                    legacyNc = datastruct(legacyCurrent, 7; var = "legacyCurrent")
                end
            else
                legacyCurrent = zeros(1, 7); legacyNc = 0
            end
            if exists(fid, "legacyInjection")
                legacyInjection::Array{Float64,2} = h5read(path.fullpath, "/legacyInjection")
                if flag
                    legacyNi = datastruct(legacyInjection, 7; var = "legacyInjection")
                else
                    legacyNi = datastruct(legacyInjection, 9; var = "legacyInjection")
                end
            else
                legacyInjection = zeros(1, 9); legacyNi = 0
            end
            if exists(fid, "legacyVoltage")
                legacyVoltage::Array{Float64,2} = h5read(path.fullpath, "/legacyVoltage")
                if flag
                    legacyNv = datastruct(legacyVoltage, 4; var = "legacyVoltage")
                else
                    legacyNv = datastruct(legacyVoltage, 5; var = "legacyVoltage")
                end
            else
                legacyVoltage = zeros(1, 9); legacyNv = 0
            end
            close(fid)
        end

        if path.extension == ".xlsx"
            xf = XLSX.openxlsx(path.fullpath, mode = "r")
            if "pmuVoltage" in XLSX.sheetnames(xf)
                start = startxlsx(xf["pmuVoltage"])
                pmuVoltage = xf["pmuVoltage"][:][start:end, :]
                if flag
                    pmuNv = datastruct(pmuVoltage, 7; var = "pmuVoltage")
                else
                    pmuNv = datastruct(pmuVoltage, 9; var = "pmuVoltage")
                end
            else
                pmuVoltage = zeros(1, 9); pmuNv = 0
            end
            if "pmuCurrent" in XLSX.sheetnames(xf)
                start = startxlsx(xf["pmuCurrent"])
                pmuCurrent = xf["pmuCurrent"][:][start:end, :]
                if flag
                    pmuNc = datastruct(pmuCurrent, 9; var = "pmuCurrent")
                else
                    pmuNc = datastruct(pmuCurrent, 11; var = "pmuCurrent")
                end
            else
                pmuCurrent = zeros(1, 11); pmuNc = 0
            end
            if "legacyFlow" in XLSX.sheetnames(xf)
                start = startxlsx(xf["legacyFlow"])
                legacyFlow = xf["legacyFlow"][:][start:end, :]
                if flag
                    legacyNf = datastruct(legacyFlow, 9; var = "legacyFlow")
                else
                    legacyNf = datastruct(legacyFlow, 11; var = "legacyFlow")
                end
            else
                legacyFlow = zeros(1, 11); legacyNf = 0
            end
            if "legacyCurrent" in XLSX.sheetnames(xf)
                start = startxlsx(xf["legacyCurrent"])
                legacyCurrent = xf["legacyCurrent"][:][start:end, :]
                if flag
                    legacyNc = datastruct(legacyCurrent, 6; var = "legacyCurrent")
                else
                    legacyNc = datastruct(legacyCurrent, 7; var = "legacyCurrent")
                end
            else
                legacyCurrent = zeros(1, 7); legacyNc = 0
            end
            if "legacyInjection" in XLSX.sheetnames(xf)
                start = startxlsx(xf["legacyInjection"])
                legacyInjection = xf["legacyInjection"][:][start:end, :]
                if flag
                    legacyNi = datastruct(legacyInjection, 7; var = "legacyInjection")
                else
                    legacyNi = datastruct(legacyInjection, 9; var = "legacyInjection")
                end
            else
                legacyInjection = zeros(1, 9); legacyNi = 0
            end
            if "legacyVoltage" in XLSX.sheetnames(xf)
                start = startxlsx(xf["legacyVoltage"])
                legacyVoltage = xf["legacyVoltage"][:][start:end, :]
                if flag
                    legacyNv = datastruct(legacyVoltage, 4; var = "legacyVoltage")
                else
                    legacyNv = datastruct(legacyVoltage, 5; var = "legacyVoltage")
                end
            else
                legacyVoltage = zeros(1, 9); legacyNv = 0
            end
            # close(xf)
        end
    end

    return Measurements(pmuVoltage, pmuCurrent, legacyFlow, legacyCurrent, legacyInjection, legacyVoltage),
        MeasurementsNum(pmuNv, pmuNc, legacyNf, legacyNc, legacyNi, legacyNv)
end

### Measurement generator power flow settings
@inbounds function gepfsettings(max, stop, react, solve)
    save = ""
    algorithm = "nr"
    main = false
    flow = false
    generator = false

    reactive = [false; true; false]
    if react == 1
        reactive = [true; true; false]
    end

    return FlowSettings(algorithm, solve, main, flow, generator, save, max, stop, reactive)
end

#### Measurement generator settings
@inbounds function gesettings(num, pmuset, pmuvariance, legacyset, legacyvariance, runflow, save)
    set = Dict{Symbol, Array{Float64,1}}()
    variance = Dict{Symbol, Array{Float64,1}}()

    ########## PMU set ##########
    if !isa(pmuset, Array)
        pmuset = [pmuset]
    end
    pmuset = vec(pmuset)

    onebyone = false
    complete = false
    for i in pmuset
        if i in ["Vi", "Ti", "Iij", "Dij"]
            onebyone = true
        end
        if i == "complete"
            complete = true
        end
    end
    if complete && onebyone
        deleteat!(pmuset, findall(x->x=="complete", pmuset))
        miss = setdiff(["Vi", "Ti", "Iij", "Dij"], pmuset)
        for i in miss
            pmuset = [pmuset; i; -200.0]
        end
    end

    if onebyone
        for i in pmuset
            if i in ["Vi", "Ti"] && num.pmuNv != 0
                push!(set, :pmuVoltage => [-100.0, -100.0])
            end
            if i in ["Iij", "Dij"] && num.pmuNc != 0
                push!(set, :pmuCurrent => [-100.0, -100.0])
            end
        end
        for (k, i) in enumerate(pmuset)
            if i == "Vi" && num.pmuNv != 0
                set[:pmuVoltage][1] = round(nextelement(pmuset, k))
            end
            if i == "Ti" && num.pmuNv != 0
                set[:pmuVoltage][2] = round(nextelement(pmuset, k))
            end
            if i == "Iij" && num.pmuNc != 0
                set[:pmuCurrent][1] = round(nextelement(pmuset, k))
            end
            if i == "Dij" && num.pmuNc != 0
                set[:pmuCurrent][2] = round(nextelement(pmuset, k))
            end
        end
    end
    if !onebyone && (num.pmuNv != 0 || num.pmuNc != 0)
        for (k, i) in enumerate(pmuset)
            if i == "redundancy"
                push!(set, Symbol(:pmu, i) => [nextelement(pmuset, k)])
                break
            end
            if i == "device"
                push!(set, Symbol(:pmu, i) => [nextelement(pmuset, k)])
                break
            end
            if i in ["complete" "optimal"]
                push!(set, Symbol(:pmu, i) => [0.0])
                break
            end
        end
    end

    ########## Legacy set ##########
    if !isa(legacyset, Array)
        legacyset = [legacyset]
    end
    legacyset = vec(legacyset)

    onebyone = false
    complete = false
    for i in legacyset
        if i in ["Pij", "Qij", "Iij", "Pi", "Qi", "Vi"]
            onebyone = true
        end
        if i == "complete"
            complete = true
        end
    end
    if complete && onebyone
        deleteat!(legacyset, findall(x->x=="complete", legacyset))
        miss = setdiff(["Pij", "Qij", "Iij", "Pi", "Qi", "Vi"], legacyset)
        for i in miss
            legacyset = [legacyset; i; -200.0]
        end
    end

    if onebyone
        for i in legacyset
            if i in ["Pij", "Qij"] && num.legacyNf != 0
                push!(set, :legacyFlow => [-100.0, -100.0])
            end
            if i == "Iij" && num.legacyNc != 0
                push!(set, :legacyCurrent => [-100.0])
            end
            if i in ["Pi", "Qi"] && num.legacyNi != 0
                push!(set, :legacyInjection => [-100.0, -100.0])
            end
            if i == "Vi" && num.legacyNv != 0
                push!(set, :legacyVoltage => [-100.0])
            end
        end
        for (k, i) in enumerate(legacyset)
            if i == "Pij" && num.legacyNf != 0
                set[:legacyFlow][1] = round(nextelement(legacyset, k))
            end
            if i == "Qij" && num.legacyNf != 0
                set[:legacyFlow][2] = round(nextelement(legacyset, k))
            end
            if i == "Iij" && num.legacyNc != 0
                set[:legacyCurrent] = [round(nextelement(legacyset, k))]
            end
            if i == "Pi" && num.legacyNi != 0
                set[:legacyInjection][1] = round(nextelement(legacyset, k))
            end
            if i == "Qi" && num.legacyNi != 0
                set[:legacyInjection][2] = round(nextelement(legacyset, k))
            end
            if i == "Vi" && num.legacyNv != 0
                set[:legacyVoltage] = [round(nextelement(legacyset, k))]
            end
        end
    end
    if !onebyone && any([num.legacyNf, num.legacyNc, num.legacyNi, num.legacyNv] .!= 0)
        for (k, i) in enumerate(legacyset)
            if i == "redundancy"
                push!(set, Symbol(:legacy, i) => [nextelement(legacyset, k)])
                break
            end
            if i == "complete"
                push!(set, Symbol(:legacy, i) => [0.0])
                break
            end
        end
    end

    ########## PMU variance ##########
    pmuvariance = vec(pmuvariance)
    onebyone = false
    complete = false
    for i in pmuvariance
        if i in ["Vi", "Ti", "Iij", "Dij"]
            onebyone = true
        end
        if i == "complete"
            complete = true
        end
    end
    if complete && onebyone
        valall = 0.0
        for (k, i) in enumerate(pmuvariance)
            if i == "complete"
                valall = nextelement(pmuvariance, k)
            end
        end
        miss = setdiff(["Vi", "Ti", "Iij", "Dij"], pmuvariance)
        for i in miss
            pmuvariance = [pmuvariance; i; valall]
        end
    end
    if onebyone
        for i in pmuvariance
            if i in ["Vi", "Ti"] && num.pmuNv != 0
                push!(variance, :pmuVoltage => [-100.0, -100.0])
            end
            if i in ["Iij", "Dij"] && num.pmuNc != 0
                push!(variance, :pmuCurrent => [-100.0, -100.0])
            end
        end
        for (k, i) in enumerate(pmuvariance)
            if i == "Vi" && num.pmuNv != 0
                variance[:pmuVoltage][1] = nextelement(pmuvariance, k)
            end
            if i == "Ti" && num.pmuNv != 0
                variance[:pmuVoltage][2] = nextelement(pmuvariance, k)
            end
            if i == "Iij" && num.pmuNc != 0
                variance[:pmuCurrent][1] = nextelement(pmuvariance, k)
            end
            if i == "Dij" && num.pmuNc != 0
                variance[:pmuCurrent][2] = nextelement(pmuvariance, k)
            end
        end
    end
    if !onebyone && (num.pmuNv != 0 || num.pmuNc != 0)
        for (k, i) in enumerate(pmuvariance)
            if i == "random"
                push!(variance, Symbol(:pmu, i) => [nextelement(pmuvariance, k), nextelement(pmuvariance, k + 1)])
                break
            end
            if i == "complete"
                push!(variance, Symbol(:pmu, i) => [nextelement(pmuvariance, k)])
                break
            end
        end
    end

    ########## Legacy variance ##########
    legacyvariance = vec(legacyvariance)
    onebyone = false
    complete = false
    for i in legacyvariance
        if i in ["Pij", "Qij", "Iij", "Pi", "Qi", "Vi"]
            onebyone = true
        end
        if i == "complete"
            complete = true
        end
    end
    if complete && onebyone
        valall = 0.0
        for (k, i) in enumerate(legacyvariance)
            if i == "complete"
                valall = nextelement(legacyvariance, k)
            end
        end
        miss = setdiff(["Pij", "Qij", "Iij", "Pi", "Qi", "Vi"], legacyvariance)
        for i in miss
            legacyvariance = [legacyvariance; i; valall]
        end
    end
    if onebyone
        for i in legacyvariance
            if i in ["Pij", "Qij"] && num.legacyNf != 0
                push!(variance, :legacyFlow => [-100.0, -100.0])
            end
            if i == "Iij" && num.legacyNc != 0
                push!(variance, :legacyCurrent => [-100.0])
            end
            if i in ["Pi", "Qi"] && num.legacyNi != 0
                push!(variance, :legacyInjection => [-100.0, -100.0])
            end
            if i == "Vi" && num.legacyNi != 0
                push!(variance, :legacyVoltage => [-100.0])
            end
        end
        for (k, i) in enumerate(legacyvariance)
            if i == "Pij" && num.legacyNf != 0
                variance[:legacyFlow][1] = nextelement(legacyvariance, k)
            end
            if i == "Qij" && num.legacyNf != 0
                variance[:legacyFlow][2] = nextelement(legacyvariance, k)
            end
            if i == "Iij" && num.legacyNc != 0
                variance[:legacyCurrent] = [nextelement(legacyvariance, k)]
            end
            if i == "Pi" && num.legacyNi != 0
                variance[:legacyInjection][1] = nextelement(legacyvariance, k)
            end
            if i == "Qi" && num.legacyNi != 0
                variance[:legacyInjection][2] = nextelement(legacyvariance, k)
            end
            if i == "Vi" && num.legacyNv != 0
                variance[:legacyVoltage] = [nextelement(legacyvariance, k)]
            end
        end
    end
    if !onebyone  && any([num.legacyNf, num.legacyNc, num.legacyNi, num.legacyNv] .!= 0)
        for (k, i) in enumerate(legacyvariance)
            if i == "random"
                push!(variance, Symbol(:legacy, i) =>  [nextelement(legacyvariance, k), nextelement(legacyvariance, k + 1)])
                break
            end
            if i == "complete"
                push!(variance, Symbol(:legacy, i) => [nextelement(legacyvariance, k)])
                break
            end
        end
    end

    return GeneratorSettings(runflow, set, variance, save)
end

### Load state estimation data from generator or input data
function loadse(args)
    data = true
    for i in args
        if typeof(i) == Tuple{Measurements, PowerSystem, Array{String,2}}
            data = false
            break
        end
    end

    return data
end

### State estimation settings
@inbounds function sesettings(args, system, max, stop, start, badset, lavset, observeset, covariance, solve, save)
    algorithm = "false"
    main = false; flow = false; estimate = false; error = false
    lavkey = 0.0; badkey = 0.0; observekey = 0.0

    for i in args
        if i in ["dc", "nonlinear", "pmu"]
            algorithm = i
        end
        if i == "lav"
            lavkey = 1.0
        end
        if i == "bad"
            badkey = 1.0
        end
        if i == "observe"
            observekey = 1.0
        end
        if i == "main"
            main = true
        end
        if i == "flow"
            flow = true
        end
        if i == "estimate"
            estimate = true
        end
        if i == "error"
            error = true
        end
    end

    if algorithm == "false"
        algorithm = "nonlinear"
        println("Invalid state estimation METHOD key. The algorithm proceeds with the nonlinear state estimation.")
    end

    startkey = 1.0; Vmin = 0.98; Vmax = 1.02; Tmin = -0.5; Tmax = 0.5
    if !isempty(start)
        if isa(start, String)
            if start == "warm"
                startkey = 1.0
            elseif start == "flat"
                startkey = 2.0
            elseif start == "random"
                startkey = 3.0
            else
                println("Invalid state estimation START key. The algorithm proceeds with the warm initial point.")
            end
        else
            for (k, i) in enumerate(start)
                if i == "Vi"
                    startV = minmax(nextelement(start, k), nextelement(start, k + 1))
                    Vmin = startV[1]; Vmax = startV[2]
                    startkey = 3.0
                end
                if i == "Ti"
                    startT = minmax(nextelement(start, k), nextelement(start, k + 1))
                    Tmin = startT[1]; Tmax = startT[2]
                    startkey = 3.0
                end
            end
        end
    end
    start = Dict(:start => startkey, :Vmin => Vmin, :Vmax => Vmax, :Tmin => Tmin, :Tmax => Tmax)

    threshold = 3.0; pass = 1.0; critical = 1e-10
    if !isempty(badset)
        badkey = 1.0
        for (k, i) in enumerate(badset)
            if i == "threshold"
                threshold = nextelement(badset, k)
            end
            if i == "pass"
                pass = nextelement(badset, k)
            end
            if i == "critical"
                critical = nextelement(badset, k)
            end
        end
    end
    bad = Dict(:bad => badkey, :treshold => threshold, :pass => pass, :critical => critical)

    optimize = 1.0
    if !isempty(lavset)
        lavkey = 1.0
        for i in lavset
            if i == "Ipopt"
                optimize = 2.0
            end
        end
    end
    lav = Dict(:lav => lavkey, :optimize => optimize)

    islands = 1; islandflow = 0; islandMax = 2000; islandBreak = 10; islandStop = 1.0; islandTreshold = 1e5
    restoreMax = 100; pivot = 1e-10; restore = 1; Pij = 0.0; Pi = 0.0; Ti = 0.0
    if !isempty(observeset)
        observekey = 1.0
        for (k, i) in enumerate(observeset)
            if i == "islandBP"
                islands = 2
            end
            if i == "flow"
                islandflow = 1
            end
            if i == "islandMax"
                islandMax = trunc(Int64, nextelement(observeset, k))
            end
            if i == "islandBreak"
                islandBreak = trunc(Int64, nextelement(observeset, k))
            end
            if i == "islandStop"
                islandStop = nextelement(observeset, k)
            end
            if i == "islandTreshold"
                islandTreshold = nextelement(observeset, k)
            end
            if i == "restoreBP"
                restore = 2
            end
            if i == "restoreMax"
                restoreMax = trunc(Int64, nextelement(observeset, k))
            end
            if i == "pivot"
                pivot = nextelement(observeset, k)
            end
            if i == "Pij"
                Pij = nextelement(observeset, k)
            end
            if i == "Pi"
                Pi = nextelement(observeset, k)
            end
            if i == "Ti"
                Ti = nextelement(observeset, k)
            end
        end
    end
    if Pij == 0.0 && Pi == 0.0 && Ti == 0.0
        Pi = 1e5
    end
    if islandflow == 1
        islands = 3
    end

    observe = Dict(:observe => observekey,
                   :islands => islands, :islandMax => islandMax, :islandBreak => islandBreak, :islandStop => islandStop,  :islandTreshold => islandTreshold,
                   :restore => restore, :restoreMax => restoreMax, :pivot => pivot, :Pij => Pij, :Pi => Pi, :Ti => Ti)

    covarinace = false
    if algorithm == "pmu" && covariance == 1
        covarinace = true
    end

    saveextension = ""
    if !isempty(save)
        path = dirname(save)
        data = basename(save)
        if isempty(data)
            dataname = (join(replace(split(system.dataname, ""), "."=>"")))
            data = string(dataname, "_results", system.extension)
        end
        saveextension = string(match(r"\.[A-Za-z0-9]+$", data).match)
        save = joinpath(path, data)
    end

    return EstimationSettings(algorithm, main, flow, estimate, error, max, stop,
        start, bad, lav, observe, covarinace, solve, save, saveextension)
end

### Load state estimation data from measurement generator
function loadsedirect(args)
    idx = 0
    for (k, i) in enumerate(args)
        if typeof(i) == Tuple{Measurements, PowerSystem, Array{String,2}}
            idx = k
            break
        end
    end
    measurements::Measurements, system::PowerSystem, info::Array{String,2} = args[idx]

    pmuNv = size(measurements.pmuVoltage, 1)
    if pmuNv == 1 && measurements.pmuVoltage[1, 1] == 0
        pmuNv = 0
    end
    pmuNc = size(measurements.pmuCurrent, 1)
    if pmuNc == 1 && measurements.pmuCurrent[1, 1] == 0
        pmuNc = 0
    end
    legacyNf = size(measurements.legacyFlow, 1)
    if legacyNf == 1 && measurements.legacyFlow[1, 1] == 0
        legacyNf = 0
    end
    legacyNc = size(measurements.legacyCurrent, 1)
    if legacyNc == 1 && measurements.legacyCurrent[1, 1] == 0
        legacyNc = 0
    end
    legacyNi = size(measurements.legacyInjection, 1)
    if legacyNi == 1 && measurements.legacyInjection[1, 1] == 0
        legacyNi = 0
    end
    legacyNv = size(measurements.legacyVoltage, 1)
    if legacyNv == 1 && measurements.legacyVoltage[1, 1] == 0
        legacyNv = 0
    end
    Nbus = size(system.bus, 1)
    Nbranch = size(system.branch, 1)
    Ngen = size(system.generator, 1)
    if Ngen == 1 && system.generator[1, 1] == 0
        Ngen = 0
    end
    Ncost = size(system.generatorcost, 1)
    if Ncost == 1 && system.generatorcost[1, 1] == 0
        Ncost = 0
    end

    return system, PowerSystemNum(Nbus, Nbranch, Ngen, Ncost),
        measurements, MeasurementsNum(pmuNv, pmuNc, legacyNf, legacyNc, legacyNi, legacyNv), info
end

### Optimal power flow settings
@inbounds function opfsettings(args, save, system, num)
    if num.Ncost == 0
        throw(ErrorException("error opening generatorcost"))
    end
    if num.Ncost != num.Ngen
        throw(ErrorException("generator and generatorcost dimensions mismatch"))
    end

    algorithm = "false"
    main = false
    flow = false
    generation = false

    for i = 1:length(args)
        if args[i] in ["dc", "ac"]
            algorithm = args[i]
        end
        if args[i] == "main"
            main = true
        end
        if args[i] == "flow"
            flow = true
        end
        if args[i] == "generation" && num.Ngen != 0
            generation = true
        end
    end

    if algorithm == "false"
        algorithm = "ac"
        println("Invalid optimal power flow METHOD key. The algorithm proceeds with the AC optimal power flow.")
    end

    if !isempty(save)
        path = dirname(save)
        data = basename(save)
        if isempty(data)
            dataname = join(replace(split(system.dataname, ""), "."=>""))
            data = string(dataname, "_results", system.extension)
        end
        save = joinpath(path, data)
    end

    return OptimalFlowSettings(algorithm, main, flow, generation, save)
end

max_stars_repo_name: mcosovic/JuliaGrid

max_issues_repo_name: mcosovic/JuliaGrid

max_forks_repo_name: mcosovic/JuliaGrid


text:
import Base.LinAlg.BLAS: axpy!
export idrs, idrs!

####

@inline function omega(t, s)
    angle = sqrt(2.)/2
    ns = vecnorm(s)
    nt = vecnorm(t)
    ts = vecdot(t,s)
    rho = abs(ts/(nt*ns))
    om = ts/(nt*nt)
    if rho < angle
        om = om*convert(typeof(om),angle)/rho
    end
    om
end

@inline linsys_op(x, A) = A*x

"""
IDRsSolver.jl
============================

The Induced Dimension Reduction method is a family of simple and fast Krylov
subspace algorithms for solving large nonsymmetric linear systems. The idea
behind the IDR(s) variant is to generate residuals that are in the nested
subspaces of shrinking dimension s.

Syntax
------

        X, h::ConvergenceHistory = idrs(A, b; s = 8, tol=sqrt(eps(typeof(real(b[1])))), maxiter = length(b)^2))

        The operator A must support multiplication with b,
        the right hand side b must support vecnorm, vecdot, copy, rand! and axpy!.
        .

Arguments
---------

       s -- dimension reduction number. Normally, a higher s gives faster convergence,
            but also  makes the method more expensive.
       tol -- tolerance of the method.
       maxiter -- maximum number of iterations
       smoothing -- (Logical) specifies if residual smoothing must be applied. Default is false.

Output
------

        X -- Approximated solution by IDR(s)
        h -- Convergence history

        The [`ConvergenceHistory`](https://github.com/JuliaLang/IterativeSolvers.jl/issues/6) type provides information about the iteration history.
            - `isconverged::Bool`, a flag for whether or not the algorithm is converged.
            - `threshold`, the convergence threshold
            - `residuals::Vector`, the value of the convergence criteria at each iteration

References
----------

    [1] IDR(s): a family of simple and fast algorithms for solving large
        nonsymmetric linear systems. P. Sonneveld and M. B. van Gijzen
        SIAM J. Sci. Comput. Vol. 31, No. 2, pp. 1035--1062, 2008
    [2] Algorithm 913: An Elegant IDR(s) Variant that Efficiently Exploits
        Bi-orthogonality Properties. M. B. van Gijzen and P. Sonneveld
        ACM Trans. Math. Software,, Vol. 38, No. 1, pp. 5:1-5:19, 2011
    [3] This file is a translation of the following MATLAB implementation:
        http://ta.twi.tudelft.nl/nw/users/gijzen/idrs.m
    [4] IDR(s)' webpage http://ta.twi.tudelft.nl/nw/users/gijzen/IDR.html
"""
idrs(A, b; s = 8, tol=sqrt(eps(typeof(real(b[1])))), maxiter = length(b)^2, smoothing=false) =
    idrs!(zerox(A, b), A, b; s=s, tol=tol, maxiter=maxiter, smoothing=smoothing)

function idrs!(x, A, b; s = 8, tol=sqrt(eps(typeof(real(b[1])))), maxiter=length(x)^2, smoothing=false)
    history = ConvergenceHistory()
    history[:tol] = tol
    reserve!(history,:resnorm, maxiter)
    idrs_method!(history, x, linsys_op, (A,), b, s, tol, maxiter; smoothing=smoothing)
    x, history
end

function idrs_method!{T}(log::ConvergenceHistory, X, op, args, C::T,
    s::Number, tol::Number, maxiter::Number; smoothing::Bool=false)

    R = C - op(X, args...)::T
    normR = vecnorm(R)
	iter = 0

    if smoothing
        X_s = copy(X)
        R_s = copy(R)
        T_s = zeros(R)
    end

    if normR <= tol           # Initial guess is a good enough solution
        return X, ConvergenceHistory(0<= res[end] < tol, tol, length(res), res)
    end

    Z = zero(C)

    P = T[rand!(copy(C)) for k in 1:s]
    U = T[copy(Z) for k in 1:s]
    G = T[copy(Z) for k in 1:s]
    Q = copy(Z)
    V = copy(Z)

    M = eye(eltype(C),s,s)
    f = zeros(eltype(C),s)
    c = zeros(eltype(C),s)

    om::eltype(C) = 1
    while normR > tol && iter < maxiter
        for i in 1:s,
            f[i] = vecdot(P[i], R)
        end
        for k in 1:s
            nextiter!(log,mvps=1)

            # Solve small system and make v orthogonal to P

            c = LowerTriangular(M[k:s,k:s])\f[k:s]
            copy!(V, G[k])
            scale!(c[1], V)

            copy!(Q, U[k])
            scale!(c[1], Q)
            for i = k+1:s
                axpy!(c[i-k+1], G[i], V)
                axpy!(c[i-k+1], U[i], Q)
            end

            # Compute new U[:,k] and G[:,k], G[:,k] is in space G_j

            #V = R - V
            scale!(-1., V)
            axpy!(1., R, V)

            copy!(U[k], Q)
            axpy!(om, V, U[k])
            G[k] = op(U[k], args...)

            # Bi-orthogonalise the new basis vectors

            for i in 1:k-1
                alpha = vecdot(P[i],G[k])/M[i,i]
                axpy!(-alpha, G[i], G[k])
                axpy!(-alpha, U[i], U[k])
            end

            # New column of M = P'*G  (first k-1 entries are zero)

            for i in k:s
                M[i,k] = vecdot(P[i],G[k])
            end

            #  Make r orthogonal to q_i, i = 1..k

            beta = f[k]/M[k,k]
            axpy!(-beta, G[k], R)
            axpy!(beta, U[k], X)

            normR = vecnorm(R)
            if smoothing
                # T_s = R_s - R
                copy!(T_s, R_s); axpy!(-1., R, T_s)

                gamma = vecdot(R_s, T_s)/vecdot(T_s, T_s)

                # R_s = R_s - gamma*T_s
                axpy!(-gamma, T_s, R_s)
                # X_s = X_s - gamma*(X_s - X)
                copy!(T_s, X_s); axpy!(-1., X, T_s); axpy!(-gamma, T_s, X_s)

                normR = vecnorm(R_s)
            end
            push!(log, :resnorm, normR)
            iter += 1
            if normR < tol || iter > maxiter
                shrink!(log)
                setconv(log, 0<=normR<tol)
                return X
            end
            if k < s
                f[k+1:s] = f[k+1:s] - beta*M[k+1:s,k]
            end

        end

        # Now we have sufficient vectors in G_j to compute residual in G_j+1
        # Note: r is already perpendicular to P so v = r
        copy!(V, R)
        Q = op(V, args...)::T
        om = omega(Q, R)
        axpy!(-om, Q, R)
        axpy!(om, V, X)

        normR = vecnorm(R)
        if smoothing
            # T_s = R_s - R
            copy!(T_s, R_s); axpy!(-1., R, T_s)

            gamma = vecdot(R_s, T_s)/vecdot(T_s, T_s)

            # R_s = R_s - gamma*T_s
            axpy!(-gamma, T_s, R_s)
            # X_s = X_s - gamma*(X_s - X)
            copy!(T_s, X_s); axpy!(-1., X, T_s); axpy!(-gamma, T_s, X_s)

            normR = vecnorm(R_s)
        end
        iter += 1
        nextiter!(log,mvps=1)
        push!(log, :resnorm, normR)
    end
    if smoothing
        copy!(X, X_s)
    end
    shrink!(log)
    setconv(log, 0<=normR<tol)
    X
end


text:
module MCMCBenchmarks

using Reexport

@reexport using Turing, CmdStan, StatsPlots, Pkg, Tracker
@reexport using Statistics, DataFrames, Random, Parameters, DynamicHMC
@reexport using LogDensityProblems, TransformVariables, Dates, FillArrays
@reexport using AdvancedHMC, ForwardDiff, Distributed, CSV, SpecialFunctions
import Base: get
import MCMCChains: Chains

include("plotting.jl")
include("Utilities.jl")
include("TOML/print.jl")
include("TOML/TOML.jl")
include("TOML/parser.jl")
include("core_benchmarking.jl")
include("helper_functions.jl")

export
  modifyConfig!,
  addChainSummary!,
  removeBurnin,
  toDict,
  addKW!,
  setreps,
  compileStanModel,
  addPerformance!,
  updateResults!,
  runSampler,
  benchmark!,
  benchmark,
  pbenchmark,
  MCMCSampler,
  AHMCNUTS,
  CmdStanNUTS,
  DHMCNUTS,
  DNNUTS,
  plotdensity,
  plotsummary,
  plotscatter,
  plotrecovery,
  nptochain,
  setprocs,
  save,
  Permutation,
  initChains,
  initStan,
  getMetadata,
  run_all_benchmarks
end # module


text:
module TestReport

using Test
using TMLE
using MLJ

@testset "Test influencecurve" begin
    @test TMLE.influencecurve([1, 1, 1], [1, 0, 1], [0.8, 0.1, 0.8], [0.8, 0.2, 0.8], 1) ==
        [0.0
        -0.9
        0.0]
end

@testset "Test summary" begin
    r1 = TMLE.QueryReport((t=[0, 1],), [1, 2, 3, 4], 1, 0.8)
    s1 = briefreport(r1)
    @test s1.query == (t = [0, 1],)
    @test s1.pvalue ≈ 5.88e-8 atol=1e-2
    @test s1.confint[1] ≈ 2.234 atol=1e-2
    @test s1.confint[2] ≈ 4.765 atol=1e-2
    @test s1.estimate == 1.0
    @test s1.initial_estimate == 0.8
    @test s1.stderror ≈ 0.645 atol=1e-2
    @test s1.mean_inf_curve == 2.5

end
end

true

max_stars_repo_name: olivierlabayle/TMLE

max_issues_repo_name: olivierlabayle/TMLE

max_forks_repo_name: olivierlabayle/TMLE


text:
d = split.(readlines("input.txt"))
println(sum(allunique.(d)))
println(sum(allunique.(map.(x -> reduce(*, sort(split(x, ""))), d))))


text:
function test_most_popular()
    println("-- Testing MostPopular recommender")

    data = DataAccessor(sparse([1 2 3; 4 5 0]))
    recommender = MostPopular(data)
    build!(recommender)
    @test ranking(recommender, 1, 1) == 2.0
    @test ranking(recommender, 1, 3) == 1.0
end

test_most_popular()


text:
module ContinuedFractions

import Base: iterate, IteratorSize, IteratorEltype, HasLength, HasEltype,
length, eltype, collect

export ContinuedFraction, quotients, convergents, ConvergentIterator

struct ContinuedFraction{T}
	quotients::Vector{T}
	ContinuedFraction{T}(qs::Vector{T}) where {T<:Integer} = new{T}(qs)
end

eltype(it::ContinuedFraction{T}) where {T} = T

quotients(cf::ContinuedFraction) = cf.quotients

struct ConvergentIterator{T}
    qs::Vector{T}
	ConvergentIterator{T}(qs::Vector{T}) where {T<:Integer} = new{T}(qs)
end
ConvergentIterator(qs::Vector{T}) where {T<:Integer} = ConvergentIterator{T}(qs)

_done(it::ConvergentIterator, state::Int) = state > length(it.qs)
function Base.iterate(it::ConvergentIterator, state::Int = 1)
	_done(it, state) && return nothing
	convergent = Rational(ContinuedFraction(it.qs[1:state]))
    convergent, state + 1
end

IteratorSize(::Type{ConvergentIterator}) = HasLength()
IteratorEltype(::Type{ConvergentIterator}) = HasEltype()

length(it::ConvergentIterator) = length(it.qs)
eltype(::Type{ConvergentIterator{T}}) where {T} = Rational{T}
eltype(it::ConvergentIterator{T})     where {T} = Rational{T}
collect(it::ConvergentIterator{T})    where {T} = collect(Rational{T}, it)

convergents(cf::ContinuedFraction) = convergents(quotients(cf))
convergents(qs::Vector{<:Integer}) = ConvergentIterator(qs)

function Base.Rational(cf::ContinuedFraction)
    qs = quotients(cf)
    isempty(qs) && return 0 // 1
    length(qs) == 1 && return qs[1] // 1

    remainder = qs[2:end]
    rat = Rational(ContinuedFraction(remainder))
    (qs[1] * rat.num + rat.den) // rat.num
end

function ContinuedFraction(rat::Rational{T}) where {T<:Integer}
    a = div(rat.num, rat.den)
    a * rat.den == rat.num && return ContinuedFraction(T[a])  # Exact!

    cf = ContinuedFraction(rat.den//(rat.num - a*rat.den))
    pushfirst!(quotients(cf), a) # insert at index 1
    cf
end

# let rationalize handle conversion from floating point
ContinuedFraction(x::BigFloat) = ContinuedFraction(rationalize(BigInt, x))
ContinuedFraction(x::AbstractFloat) = ContinuedFraction(rationalize(x))

ContinuedFraction(x::Integer) = ContinuedFraction([x])
ContinuedFraction(qs::Vector{T}) where {T<:Integer} = ContinuedFraction{T}(qs)

end # module ContinuedFractions


text:
info("Loading required packages")
using advanced_timer, MATLAB, MUMPS

info("Loading julia code")

include("linear_system_solvers/linear_system_solvers.jl")
include("core/core.jl")
include("helpers/eigenvalues.jl")
include("helpers/woodbury.jl")
include("newton_direction/newton_direction.jl")
include("line_search/line_search.jl")
include("strategies/strategies.jl")
include("homogeneous.jl")
include("benchmarking/benchmarking.jl")

max_forks_repo_name: ohinder/Everest


text:
struct ConTab{T1<:AbstractMatrix, T2<:AbstractVector}
    X::T1
    mx::T2
    my::T2
end

function contab(x, y; normalize = true, correction = 0.5, verbose = true)
    t = [count((i .=== x) .& (j .=== y)) for i in sort(unique(skipmissing(x))), j in sort(unique(skipmissing(y)))]
    t = Matrix{Float64}(t)
    t[t .== 0] .= correction / sum(t)
    if normalize
        t = t ./ sum(t)
        if any(t .< eps())
            @warn "Some normalized frequency are smaller than `eps()`. Use larger `correction` value."
        end
        if any(t .< 0.01) && verbose
            @warn "Probabilities of some elements in the contigency tabel, $(findall(x -> x < 0.01, t)) are less than 0.01. The estimation of polychoric correlation may be unstable. For `snpolyc`, use a smaller number of `step` to solve it."
        end
    end
    return ConTab(t, sum(t; dims = 2)[:], sum(t; dims = 1)[:])
end

function contab(t; normalize = true, correction = 0.5, verbose = true)
    t = Matrix(t)
    t[t .== 0] .= correction / sum(t)
    if normalize
        t = t ./ sum(t)
        if any(t .< eps())
            @warn "Some normalized frequency are smaller than `eps()`. Use more large `correction` value."
        end
        if any(t .< 0.01) && verbose
            @warn "Some elements in the contigency tabel, $(findall(x -> x < 0.01, t)). The estimation of polychoric correlation may be unstable. For `snpolyc`, use a larger number of `step` to solve it."
        end
    end
    return contab(t, sum(t; dims = 2)[:], sum(t; dims = 1)[:])
end

# Calculate all mass, not density.

"""
    ppdf(f, max, min, step)

Compute the volume of the rectangular area ∈ [min, max] by using trapezoidal rule.
"""
function ppdf(f, max, min, step)
    x = range(min[1], stop = max[1], step = step)
    y = range(min[2], stop = max[2], step = step)
    return trapz((x, y), [f([i, j]) for i in x, j in y])
    # return hcubature(f, max, min)[1]
end

"""
Compute, more efficiently than using the trapezoidal rule, the volume of the bivariate normal distributions.
"""
function ppdf(f, m)
    g(h, k, θ) = exp(-0.5 * (h^2 + k^2 - 2*h*k*θ)/(1-θ^2)) / sqrt(1-θ^2 )
    Φ(x) = cdf(Normal(0, 1), x)
    return quadgk(t->g(m[1], m[2], t), 0, f.Σ[1,2])[1] / 2π + Φ(m[1])*Φ(m[2])
    # v = range(0, f.Σ[1,2], length = 20)
    # V = @trapz v x g(m[1], m[2], x)
    # return V / 2π + Φ(m[1])*Φ(m[2])
end

# Numerical integration with sophisticated 1dim version
function H(a, b, c, d, mvd::MvNormal)
    A1 = a == -Inf || c == -Inf ?  0.0 : ppdf(mvd, [a, c])
    A2 = a == -Inf ? 0.0 : ppdf(mvd, [a, d])
    A3 = c == -Inf ? 0.0 : ppdf(mvd, [b, c])
    A4 = d == Inf && b == Inf ? 1.0 : ppdf(mvd, [b, d])
    return A1 - A2 - A3 + A4
end

# Loss function
function loss(t::ConTab, ξx, ξy, d::MultivariateDistribution)
    I, J = size(t.X)
    h = zeros(Float64, I, J)
    @fastmath for i in 1:I, j in 1:J
        if i == I && j == J && sum(h) < 1.0
            h[i, j] = 1.0 - sum(h)
        elseif i == I
            h[i, j] = t.my[j] - sum(h[1:end-1, j])[1]
        elseif j == J
            h[i, j] = t.mx[i] - sum(h[i, 1:end-1])[1]
        else
            h[i, j] = H(ξx[i], ξx[i+1], ξy[j], ξy[j+1], d)
        end
    end
    if any(h .< 0.0)
        h = abs.(h)
    end
    h = h ./ sum(h)
    return - sum(t.X .* log.(h))
end

function polyc(t::ConTab)
    # Fix marginal freq
    cumx = [0.0; cumsum(t.mx)[1:end-1]; 1.0]
    cumy = [0.0; cumsum(t.my)[1:end-1]; 1.0]
    ξx = quantile.(Normal(0, 1), cumx)
    ξy = quantile.(Normal(0, 1), cumy)
    opt = optimize(ρ -> loss(t, ξx, ξy, MvNormal([0, 0], [1.0 ρ; ρ 1.0])), -1.0, 1.0, Brent())
    return (ρ = opt.minimizer, τ₁ = ξx, τ₂ = ξy)
end

function polyc(x, y; verbose = true)
    tab = contab(x, y; verbose = verbose)
    # Fix marginal freq
    polyc(tab)
end

"""
    polycor(X::Union{AbstractDataFrame, AbstractMatrix}; lower_tri = false)
Estimate polychoric correlation matrix from X.

# Arguments
- `X` Matrix contains categorical vectors.
- `lower_tri` logical.

# Examples
```julia

using Distributions
# Support function
function discretize(x, th)
    sum(th .< x)
end
# Gen sim data.
raw = rand(MvNormal([0, 0], [1 .5; .5 1]), 10000)'
τ1 = sort!(rand(Uniform(-2, 2), 3))
τ2 = sort!(rand(Uniform(-2, 2), 4))
obs = [discretize.(raw[:, 1], Ref(τ1)) discretize.(raw[:, 2], Ref(τ2)) ]

polycor(obs)
```
"""
function polycor(X)
    J = size(X, 2)
    r = Matrix{AbstractFloat}(undef, J, J)
    r[diagind(r)] .= 1.0
    for i in 1:J
        x = @view X[:, i]
        for j in i:J
            y = @view X[:, j]
            r[j, i] = polyc(x, y; verbose = false).ρ
            r[i, j] = r[j, i]
        end
    end
    return r
end

function replace_diagonal!(M)
    offdiagind(A) = (ι for ι in CartesianIndices(A) if ι[1] ≠ ι[2])
    L = maximum(M[collect(offdiagind(M))])
    M[diagind(M)] .= L
end


text:
"""
```
dualarray(a::AbstractArray, b::AbstractArray)

```

returns an `Array` that has the size of input `a` and an element type consistent with the element types
of the inputs `a` and `b`. For example, suppose you want to write a function of the form

```
julia> function f(a, b)
           F = similar(a)
           F[1] = a[1] * b[1]
           F[2] = a[2] + b[2]
       end
```

If you were to automatically differentiate `f` with respect to `b`, then `F` would not have the correct element type since its type
will be the same as `a`. Rather than require the user to write

```
julia> F = if eltype(b) <: ForwardDiff.Dual
           similar(a, eltype(b))
       else
           similar(a)
       end
```

the user can use `dualvector` to write

```jldoctest
julia> a = rand(3)
julia> b = ones(ForwardDiff.Dual, 5)
julia> F = dualvector(a, b)
3-element Array{ForwardDiff.Dual,1}:
 #undef
 #undef
 #undef
```

Note that the element types of `a` and `b` must be subtypes of `Real` (or else `ForwardDiff` will not work).
"""
@inline dualarray(a::AbstractArray{<: ForwardDiff.Dual}, b::AbstractArray{<: ForwardDiff.Dual}) = similar(a)
@inline dualarray(a::AbstractArray{<: ForwardDiff.Dual}, b::AbstractArray{<: Real})             = similar(a)
@inline dualarray(a::AbstractArray{<: Real}, b::AbstractArray{<: Real})                         = similar(a)
@inline dualarray(a::AbstractArray{<: Real}, b::AbstractArray{<: ForwardDiff.Dual})             = similar(a, eltype(b))

"""
```
dualvector(a::AbstractVector, b::AbstractVector)
```

has the same behavior of `dualarray` but acts specifically on vectors. This function
is primarily for user convenience.
"""
@inline dualvector(a::AbstractVector{<: ForwardDiff.Dual}, b::AbstractVector{<: ForwardDiff.Dual}) = similar(a)
@inline dualvector(a::AbstractVector{<: ForwardDiff.Dual}, b::AbstractVector{<: Real})             = similar(a)
@inline dualvector(a::AbstractVector{<: Real}, b::AbstractVector{<: Real})                         = similar(a)
@inline dualvector(a::AbstractVector{<: Real}, b::AbstractVector{<: ForwardDiff.Dual})             = similar(a, eltype(b))


text:
# Julia wrapper for header: /usr/include/scip/pub_sol.h
# Automatically generated using Clang.jl wrap_c

function SCIPsolGetOrigin(sol)
    ccall((:SCIPsolGetOrigin, libscip), SCIP_SOLORIGIN, (Ptr{SCIP_SOL},), sol)
end

function SCIPsolIsOriginal(sol)
    ccall((:SCIPsolIsOriginal, libscip), UInt32, (Ptr{SCIP_SOL},), sol)
end

function SCIPsolIsPartial(sol)
    ccall((:SCIPsolIsPartial, libscip), UInt32, (Ptr{SCIP_SOL},), sol)
end

function SCIPsolGetOrigObj(sol)
    ccall((:SCIPsolGetOrigObj, libscip), Cdouble, (Ptr{SCIP_SOL},), sol)
end

function SCIPsolGetTime(sol)
    ccall((:SCIPsolGetTime, libscip), Cdouble, (Ptr{SCIP_SOL},), sol)
end

function SCIPsolGetRunnum(sol)
    ccall((:SCIPsolGetRunnum, libscip), Cint, (Ptr{SCIP_SOL},), sol)
end

function SCIPsolGetNodenum(sol)
    ccall((:SCIPsolGetNodenum, libscip), Clonglong, (Ptr{SCIP_SOL},), sol)
end

function SCIPsolGetDepth(sol)
    ccall((:SCIPsolGetDepth, libscip), Cint, (Ptr{SCIP_SOL},), sol)
end

function SCIPsolGetHeur(sol)
    ccall((:SCIPsolGetHeur, libscip), Ptr{SCIP_HEUR}, (Ptr{SCIP_SOL},), sol)
end

function SCIPsolSetHeur(sol, heur)
    ccall((:SCIPsolSetHeur, libscip), Cvoid, (Ptr{SCIP_SOL}, Ptr{SCIP_HEUR}), sol, heur)
end

function SCIPsolGetIndex(sol)
    ccall((:SCIPsolGetIndex, libscip), Cint, (Ptr{SCIP_SOL},), sol)
end

function SCIPsolGetAbsBoundViolation(sol)
    ccall((:SCIPsolGetAbsBoundViolation, libscip), Cdouble, (Ptr{SCIP_SOL},), sol)
end

function SCIPsolGetRelBoundViolation(sol)
    ccall((:SCIPsolGetRelBoundViolation, libscip), Cdouble, (Ptr{SCIP_SOL},), sol)
end

function SCIPsolGetAbsIntegralityViolation(sol)
    ccall((:SCIPsolGetAbsIntegralityViolation, libscip), Cdouble, (Ptr{SCIP_SOL},), sol)
end

function SCIPsolGetAbsLPRowViolation(sol)
    ccall((:SCIPsolGetAbsLPRowViolation, libscip), Cdouble, (Ptr{SCIP_SOL},), sol)
end

function SCIPsolGetRelLPRowViolation(sol)
    ccall((:SCIPsolGetRelLPRowViolation, libscip), Cdouble, (Ptr{SCIP_SOL},), sol)
end

function SCIPsolGetAbsConsViolation(sol)
    ccall((:SCIPsolGetAbsConsViolation, libscip), Cdouble, (Ptr{SCIP_SOL},), sol)
end

function SCIPsolGetRelConsViolation(sol)
    ccall((:SCIPsolGetRelConsViolation, libscip), Cdouble, (Ptr{SCIP_SOL},), sol)
end


text:
testOnAllBackends("Create Array") do af
	println("\tconstant")
	afArr = constant(af, 1.1f0, 2, 2)
	@test host(afArr) == [[1.1f0, 1.1f0] [1.1f0, 1.1f0]]

	afArr = constant(af, one(Int64), 4)
	@test host(afArr) == [one(Int64), one(Int64), one(Int64), one(Int64)]

	afArr = constant(af, zero(UInt64), 1, 2)
	@test host(afArr) == [zero(UInt64) zero(UInt64)]

	println("\trandu")
	afArr = randu(af, Float32, 1, 2)
	@test host(afArr .>= 0) == [true true]
	@test host(afArr .<= 1) == [true true]

	println("\trandn")
	afArr = randn(af, Float32, 1, 2)
	@test host(afArr .>= -4.0f0) == [true true]
	@test host(afArr .<= 4.0f0) == [true true]

	println("\tlookup")
	afArr = array(af, [[1, 2] [3, 4]])
	result = lookup(afArr, array(af, [1]), 0)
	@test host(result) == [2 4]
	result = lookup(afArr, array(af, [0]), 1)
	@test host(result) == [1, 2]

	println("\ttranspose")
	afArr = array(af, [[1, 2] [3, 4]])
	result = host(transpose(afArr))

	@test result == [[1, 3] [2, 4]]

	afArr = array(af, [1, 2, 3, 4])
	result = host(transpose(afArr))

	@test result == [1 2 3 4]

	afArr = array(af, [[1, 2] [3, 4]])
	transpose!(afArr)

	@test host(afArr) == [[1, 3] [2, 4]]

	println("\trange")
	afArr = range(af, Float32, [5])
	@test host(afArr) == [0.0f0,1.0f0,2.0f0,3.0f0,4.0f0]

	afArr = range(af, Float32, [3, 3])
	@test host(afArr) == [[0.0f0, 1.0f0, 2.0f0] [0.0f0, 1.0f0, 2.0f0] [0.0f0, 1.0f0, 2.0f0]]

	afArr = range(af, Float32, [3, 3], 1)
	@test host(afArr) == [[0.0f0, 0.0f0, 0.0f0] [1.0f0, 1.0f0, 1.0f0] [2.0f0, 2.0f0, 2.0f0]]

end

max_stars_repo_name: unbornchikken/julia-ml-proto

max_issues_repo_name: unbornchikken/julia-ml-proto

max_forks_repo_name: unbornchikken/julia-ml-proto


text:
# Example code to print list of built HVDC branches and converters
function display_results_tnep(result)
    built_cv = []
    built_br = []
    for (c, conv) in result["solution"]["convdc_ne"]
        if isapprox(conv["isbuilt"] , 1; atol = 0.01)
            print("Conv: $c \n")
            push!(built_cv,c)
        end
    end
    for (b, branch) in result["solution"]["branchdc_ne"]
        if isapprox(branch["isbuilt"] , 1; atol = 0.01)
            print("Branch: $b \n")
            push!(built_br,b)
        end
    end
    return built_cv, built_br
end


text:
module SubModuleRosenBrock

using Test
using F1Method
using LinearAlgebra
using SciMLBase
using ForwardDiff

# Set up:
# - overload `SteadyStateProblem` constructor
# - overload `solve` function
# - define solver algorithm (basic Newton here)
# - define type for that algorithm (here `MyAlg`)

function newton_solve(F, ∇ₓF, x; Ftol=1e-10)
    while norm(F(x)) ≥ Ftol
        x .-= ∇ₓF(x) \ F(x)
    end
    return x
end

# Create a type for the solver's algorithm
struct MyAlg <: SciMLBase.AbstractSteadyStateAlgorithm end

# Overload SciMLBase's solve function
function SciMLBase.solve(prob::SciMLBase.AbstractSteadyStateProblem,
                          alg::MyAlg;
                          Ftol=1e-10)
    # Define the functions according to SciMLBase.SteadyStateProblem type
    p = prob.p
    x0 = copy(prob.u0)
    dx, df = copy(x0), copy(x0)
    F(x) = prob.f(x, p)
    ∇ₓF(x) = prob.f.jac(x, p)
    # Compute `u_steady` and `resid` as per SciMLBase using my algorithm
    x_steady = newton_solve(F, ∇ₓF, x0, Ftol=Ftol)
    resid = F(x_steady)
    # Return the common SciMLBase solution type
    SciMLBase.build_solution(prob, alg, x_steady, resid; retcode=:Success)
end

# This is quasi derived from the Rosenbrock function
# whereby F(x,p) is almost the gradient of the Rosenbrock function
# except the Hessian matroix is made non-symmetric by slightly altering
# F(x,p)

# State function ≈ gradient of Rosenbrock.
# With p = [a, b] the parameters defining the constants
# in Rosenbrocks function, and the slight alteration to make the
# Hessian matrix non-symmetric.
# Specifically, if r(x,y) denotes the rosenbrock function,
#     F₁([x, y], p) = ∂r/∂x
#     F₂([x, y], p) = 0.5 * ∂r/∂y
statefun(x,p) = [
    -2 * (p[1] - x[1]) - 4 * p[2] * (x[2] - x[1]^2) * x[1]
    p[2] * (x[2] - x[1]^2)
]
F = ODEFunction(statefun, jac = (x,p) -> ForwardDiff.jacobian(x -> statefun(x, p), x))

# Define mismatch function f(x,p) and its derivative ∇ₓf(x,p)
# (Note ∇ₓF and ∇ₓf are required by the F1 method)
function state_mismatch(x)
    δ(x) = x .- 1
    return 0.5δ(x)'δ(x)
end
function parameter_mismatch(p)
    δ(p) = log.(p)
    return 0.5δ(p)'δ(p)
end
f(x,p) = state_mismatch(x) + parameter_mismatch(p)
∇ₓf(x,p) = ForwardDiff.jacobian(x -> [f(x,p)], x)

# Initialize the starting state and parameters
x₀ = rand(2)
p₀ = rand(2)
# Initialize the memory cache for storing reusable objects
mem = F1Method.initialize_mem(F, ∇ₓf, x₀, p₀, MyAlg())

# Define the functions via the F1 method
F1_objective(p) = F1Method.objective(f, F, mem, p, MyAlg())
F1_gradient(p) = F1Method.gradient(f, F, ∇ₓf, mem, p, MyAlg())
F1_Hessian(p) = F1Method.hessian(f, F, ∇ₓf, mem, p, MyAlg())

# Define the exact solution and objective for comparison
exact_solution(p) = [p[1], p[1]^2]
exact_objective(p) = f(exact_solution(p), p)

# Run the tests
@testset "objective" begin
    @test exact_objective(p₀) ≈ F1_objective(p₀)
    @test exact_objective(2p₀) ≈ F1_objective(2p₀)
end
@testset "gradient" begin
    @test ForwardDiff.jacobian(p -> [exact_objective(p)], p₀) ≈ F1_gradient(p₀)
    @test ForwardDiff.jacobian(p -> [exact_objective(p)], 2p₀) ≈ F1_gradient(2p₀)
end
@testset "Hessian" begin
    @test ForwardDiff.hessian(p -> exact_objective(p), p₀) ≈ F1_Hessian(p₀)
    @test ForwardDiff.jacobian(p -> [exact_objective(p)], 2p₀) ≈ F1_gradient(2p₀)
end

end # submodule


text:
module BenchOverhead

using BenchmarkTools
using Transducers

const SUITE = BenchmarkGroup()

# https://github.com/JuliaFolds/Transducers.jl/issues/307
SUITE["foldl"] = @benchmarkable foldl(+, Map(identity), $([1, 2, 3, 4]))
SUITE["reduce_basecase"] =
    @benchmarkable reduce(+, Map(identity), $([1, 2, 3, 4]), basesize = 1000)

end  # module
BenchOverhead.SUITE


text:
function fib(n)
    n1 = BigInt(1)
    n2 = BigInt(1)
    for i in 1:(n-2)
        n1, n2 = n2, n1 + n2
    end
    n2
end

n = parse(Int64, ARGS[1])
println("res = ", fib(n))

max_stars_repo_name: juliancoffee/fib

max_issues_repo_name: juliancoffee/fib

max_forks_repo_name: juliancoffee/fib