text:
function maketablea1(panel::DataFrame; trialrunonly::Bool=trialrunonly,
    runpanela::Bool=true,
    runpanelb::Bool=true,
    runpanelc::Bool=true,
    runpaneld::Bool=true
    )::Nothing
  #load the data we will work with
  #panel::DataFrame = constructportfolios(refreshportfolios=false, trialrunonly=trialrunonly)
  sort!(panel, [:permno, :date])
  apanel::DataFrame = panel[panel.annflag,:]

  ################common data
  cqports::PortfolioConstructions = leverageconstructions()
  TA1_FNAMES::Vector{Symbol} = CQUART_FNAMES_D
  TA1_FVALS::Vector{Symbol} = CQUART_FVALS_D
  TA1_FGROUPS::Vector{Symbol} = CQUART_FGROUPS_D
  TA1_FNAMES_ALT::Vector{Symbol} = CQUART_FNAMES_ALT_D
  TA1_FGROUPS_ALT::Vector{Symbol} = CQUART_FGROUPS_ALT_D

  ###################PANEL C############
  TA1_PANEL_A_ROWKEYS::Vector{Symbol} = [
    :constantbrsd, :plusxmktbrsd, :plusumdbrsd,
    :constantbrse, :plusxmktbrse, :plusumdbrse,
    :constantbrsx, :plusxmktbrsx, :plusumdbrsx,
  ]

  TA1_PANEL_A_ROWLABEL::Dict = Dict(
    :constantbrsd=>"BRS \$\\Delta\$Debt",
    :plusxmktbrsd=>"BRS+XMKT",
    :plusumdbrsd=>"BRS+FF5+UMD",
    :constantbrse=>"\\midrule BRS \$\\Delta\$Equity",
    :plusxmktbrse=>"BRS+XMKT",
    :plusumdbrse=>"BRS+FF5+UMD",
    :constantbrsx=>"\\midrule BRS \$\\Delta\$External",
    :plusxmktbrsx=>"BRS+XMKT",
    :plusumdbrsx=>"BRS+FF5+UMD",
    )

  TA1_PANEL_A_CONTROLS::Dict = Dict(
    :constantbrsd=>[:P_brsd],
    :plusxmktbrsd=>[:P_brsd, :mkt],
    :plusumdbrsd=>[:P_brsd, :mkt, :P_ff5_smb,	:P_ff5_hml,	:P_ff5_rmw,	:P_ff5_cma, :P_ff3m_umd],
    :constantbrse=>[:P_brse],
    :plusxmktbrse=>[:P_brse, :mkt],
    :plusumdbrse=>[:P_brse, :mkt, :P_ff5_smb,	:P_ff5_hml,	:P_ff5_rmw,	:P_ff5_cma, :P_ff3m_umd],
    :constantbrsx=>[:P_brsx],
    :plusxmktbrsx=>[:P_brsx, :mkt],
    :plusumdbrsx=>[:P_brsx, :mkt, :P_ff5_smb,	:P_ff5_hml,	:P_ff5_rmw,	:P_ff5_cma, :P_ff3m_umd],
    )

  TA1_PANEL_A_FOCALS::Dict = Dict(
  :constantbrsd=>TA1_FNAMES,
  :plusxmktbrsd=>TA1_FNAMES,
  :plusumdbrsd=>TA1_FNAMES,
  :constantbrse=>TA1_FNAMES,
  :plusxmktbrse=>TA1_FNAMES,
  :plusumdbrse=>TA1_FNAMES,
  :constantbrsx=>TA1_FNAMES,
  :plusxmktbrsx=>TA1_FNAMES,
  :plusumdbrsx=>TA1_FNAMES)

  TA1_PANEL_A_NAME::String = "ta1panela-$(REPLICATION_TYPE[])-$(OUT_SUFFIX[])"

  TA1_PANEL_A_MIN_ROWS::Int = 24

  runpanela && tablexpanelC(panel,
    fvals=TA1_FVALS,
    fgroups=TA1_FGROUPS,
    Fnames = TA1_FNAMES,
    rowkeys=TA1_PANEL_A_ROWKEYS,
    rowlabel=TA1_PANEL_A_ROWLABEL,
    controls=TA1_PANEL_A_CONTROLS,
    focals=TA1_PANEL_A_FOCALS,
    tablename=TA1_PANEL_A_NAME,
    minrowspergroup = TA1_PANEL_A_MIN_ROWS
  )

  #########################PANEL B#########################
  TA1_PANEL_B_CONTROLS::Dict = Dict(
    :brsdquartiles=>[:brsd],
    :brsdquartiles10ctrl=>[T1_CONTROL10; :brsd],
    :brsequartiles=>[:brse],
    :brsequartiles10ctrl=>[T1_CONTROL10; :brse],
    :brsxquartiles=>[:brsx],
    :brsxquartiles10ctrl=>[T1_CONTROL10; :brsx],
    :brsdvalues=>[:brsd],
    :brsdvalues10ctrl=>[T1_CONTROL10; :brsd],
    :brsevalues=>[:brse],
    :brsevalues10ctrl=>[T1_CONTROL10; :brse],
    :brsxvalues=>[:brsx],
    :brsxvalues10ctrl=>[T1_CONTROL10; :brsx],)

  TA1_PANEL_B_RET::Dict = Dict(
    :brsdquartiles=>:ret,
    :brsdquartiles10ctrl=>:ret,
    :brsequartiles=>:ret,
    :brsequartiles10ctrl=>:ret,
    :brsxquartiles=>:ret,
    :brsxquartiles10ctrl=>:ret,
    :brsdvalues=>:ret,
    :brsdvalues10ctrl=>:ret,
    :brsevalues=>:ret,
    :brsevalues10ctrl=>:ret,
    :brsxvalues=>:ret,
    :brsxvalues10ctrl=>:ret)

  TA1_PANEL_B_ROWKEYS::Vector{Symbol} = [
    :brsdquartiles,
    :brsdquartiles10ctrl,
    :brsequartiles,
    :brsequartiles10ctrl,
    :brsxquartiles,
    :brsxquartiles10ctrl,
    :brsdvalues,
    :brsdvalues10ctrl,
    :brsevalues,
    :brsevalues10ctrl,
    :brsxvalues,
    :brsxvalues10ctrl,]

  TA1_PANEL_B_ROWLABEL::Dict = Dict(
    :brsdquartiles=>"Leverage Quartiles+\$\\Delta\$Debt",
    :brsdquartiles10ctrl=>"... +10 vars",
    :brsequartiles=>"Leverage Quartiles+\$\\Delta\$Equity",
    :brsequartiles10ctrl=>"... +10 vars",
    :brsxquartiles=>"Leverage Quartiles+\$\\Delta\$External",
    :brsxquartiles10ctrl=>"... +10 vars",
    :brsdvalues=>"\\midrule Leverage Values+\$\\Delta\$Debt",
    :brsdvalues10ctrl=>"... +10 vars",
    :brsevalues=>"Leverage Values+\$\\Delta\$Equity",
    :brsevalues10ctrl=>"... +10 vars",
    :brsxvalues=>"Leverage Values+\$\\Delta\$External",
    :brsxvalues10ctrl=>"... +10 vars")

  TA1_PANEL_B_FOCALS::Dict = Dict(
    :brsdquartiles=>TA1_FGROUPS,
    :brsdquartiles10ctrl=>TA1_FGROUPS,
    :brsequartiles=>TA1_FGROUPS,
    :brsequartiles10ctrl=>TA1_FGROUPS,
    :brsxquartiles=>TA1_FGROUPS,
    :brsxquartiles10ctrl=>TA1_FGROUPS,
    :brsdvalues=>TA1_FVALS,
    :brsdvalues10ctrl=>TA1_FVALS,
    :brsevalues=>TA1_FVALS,
    :brsevalues10ctrl=>TA1_FVALS,
    :brsxvalues=>TA1_FVALS,
    :brsxvalues10ctrl=>TA1_FVALS,  )

  TA1_PANEL_B_NAME::String = "ta1panelb-$(REPLICATION_TYPE[])-$(OUT_SUFFIX[])"

  TA1_PANEL_B_MIN_ROWS::Int = 24

  runpanelb && tablexpanelD(panel,
      fvals=TA1_FVALS,
      fgroups=TA1_FGROUPS,
      ret = TA1_PANEL_B_RET,
      rowkeys=TA1_PANEL_B_ROWKEYS,
      rowlabel=TA1_PANEL_B_ROWLABEL,
      controls=TA1_PANEL_B_CONTROLS,
      focals=TA1_PANEL_B_FOCALS,
      tablename=TA1_PANEL_B_NAME,
      minrowspergroup = TA1_PANEL_B_MIN_ROWS
    )

  return nothing
end

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
# License for this file: MIT (expat)
# Copyright 2017-2018, DLR Institute of System Dynamics and Control

"""
    module Simulate_StateSelection

DAE-model to test manual state selection of simple multibody system.
"""
module Simulate_StateSelection

include(joinpath("models", "StateSelection.jl"))
import .StateSelection
import ModiaMath

model  = StateSelection.Model()
result = ModiaMath.simulate!(model, stopTime=1.0, log=true)

ModiaMath.plot(result, [:s, :sd, ("f[1]", "f[2]", "f[3]")], heading="simulationWithoutMacro/Simulate_StateSelection.jl")

end

ext: jl


text:
using AccretionFormulae
using Plots
using Measures


text:
"""
    AdditiveWhiteGaussianNoise <: AbstractImageNoise
    AdditiveWhiteGaussianNoise([μ=0.0], σ)

apply white gaussian noise to image

For gray images, it uses the following formula:

    out = clamp01.(in .+ σ .* randn(size(in)) .+ μ)

RGB images are treated as 3D-gray images, generic Color3 images will be
converted to RGB images first.

# Examples
```julia
img = testimage("lena_gray_256")
n = AdditiveWhiteGaussianNoise(0.1)
out = apply_noise(img, n)
```

See also: [`apply_noise`](@ref), [`apply_noise!`](@ref)

# References
[1] Wikipedia contributors. (2019, March 8). Additive white Gaussian noise. In _Wikipedia, The Free Encyclopedia_. Retrieved 14:32, June 9, 2019, from https://en.wikipedia.org/w/index.php?title=Additive_white_Gaussian_noise&oldid=886818982
"""
struct AdditiveWhiteGaussianNoise <: AbstractImageNoise
    """mean"""
    μ::Float64
    """standard deviation"""
    σ::Float64

    function AdditiveWhiteGaussianNoise(μ, σ)
        σ >= zero(σ) || throw(ArgumentError("std σ should be non-negative"))
        new(μ, σ)
    end
end
const AWGN = AdditiveWhiteGaussianNoise
AWGN(σ::Real) = AWGN(zero(σ), σ)

show(io::IO, n::AWGN) = println(io, "AdditiveWhiteGaussianNoise(μ=", n.μ, ", σ=", n.σ, ")")

function (n::AWGN)(out::AbstractArray{T},
                   in::GenericGrayImage;
                   rng::Union{AbstractRNG, Nothing} = nothing
                   ) where T<:Number
    if rng === nothing
        _rand = x->randn(floattype(T), size(x))
    else
        _rand = x->randn(rng, floattype(T), size(x))
    end
    out .= clamp01.(in .+ n.σ .* _rand(in) .+ n.μ)
end

for T in (AbstractGray, AbstractRGB)
    @eval function (n::AWGN)(out::AbstractArray{<:$T},
                             in::GenericImage;
                             rng::Union{AbstractRNG, Nothing} = nothing)
        n(channelview(out), channelview(in); rng=rng)
    end
end

# Since generic Color3 aren't vector space, they are converted to RGB images
(n::AWGN)(out::AbstractArray{<:Color3},
         in::GenericImage;
         rng::Union{AbstractRNG, Nothing} = nothing) =
    n(channelview(of_eltype(RGB, out)),
      channelview(of_eltype(RGB, in));
      rng=rng)


text:
module Reachable2

struct B end

end


text:
using HttpServer
using WebSockets
using JSON
using JLD
using Jukai
using Jukai.Tokenization
using Jukai.POSTagging
using Merlin

include("conf.jl")

const filepath = dirname(@__FILE__)
const clients = Dict()
const tokenizer = load(joinpath(filepath,"data/tokenizer_50.jld"), "tokenizer")
const conf = begin
  e = readconf(joinpath(filepath,"conf/visual.conf"))
  d = Dict("entity_types" => e)
  JSON.json(d)
end

const postagger = begin
  println("training postagger...")
  p = POSTagger()
  POSTagging.train(p, joinpath(filepath,"data/wsj_00-18.conll"), joinpath(filepath,"data/wsj_00-18.conll"))
  println("finish.")
  p
end

wsh = WebSocketHandler() do req, client
  println("Client: $(client.id) is connected.")
  write(client, conf)
  while true
    #println("Request from $(client.id) recieved.")
    msg = bytestring(read(client))
    length(msg) > 1000 && continue
    chars = convert(Vector{Char}, msg)
    doc = Tokenization.decode(tokenizer, chars)
    POSTagging.decode!(postagger, doc)
    doc = map(x -> map(to_dict, x), doc)
    res = JSON.json(doc)
    #println(res)
    write(client, res)
  end
end

onepage = readall(joinpath(dirname(@__FILE__), "index.html"))
httph = HttpHandler() do req::Request, res::Response
  Response(onepage)
end
server = Server(httph, wsh)
run(server, 3000)

#=
const sample_tokens = [
  Token(1, 6, 1, "Pierre", 1, "NNP", -2, 2),
  Token(2, "Vinken", 2, "NNP", -1, 8),
  Token(3, ",", 3, ",", -2, 2),
  Token(4, "61", 4, "CD", -4, 5),
  Token(5, "years", 5, "NNS", -3, 6),
  Token(6, "old", 6, "JJ", -2, 2),
  Token(7, ",", 7, ",", -2, 2),
  Token(8, "will", 8, "MD", 0, 0),
  Token(9, "join", 9, "VB", -1, 8),
  Token(10, "the", 10, "DT", -3, 11),
  Token(11, "board", 11, "NN", -2, 9),
  Token(12, "Nov.", 12, "NNP", -2, 9),
  Token(13, "29", 13, "CD", -3, 12),
  Token(14, ".", 13, ".", -1, 8)]
=#

max_stars_repo_name: hshindo/JukaiWeb

max_issues_repo_name: hshindo/JukaiWeb

max_forks_repo_name: hshindo/JukaiWeb


text:
for (func, typ) in [(:clblasChemv, CL_float2),
                    (:clblasZhemv, CL_double2)]

    @eval @blasfun $func(order::clblasOrder, uplo::clblasUplo,
                              N::Csize_t,
                              alpha::$typ,
                              A::cl.CL_mem, offa::Csize_t, lda::Csize_t,
                              X::cl.CL_mem, offx::Csize_t, incx::Cint,
                              beta::$typ,
                              Y::cl.CL_mem, offy::Csize_t, incy::Cint,
                              n_queues::cl.CL_uint,
                              queues::Ptr{cl.CL_command_queue},
                              n_events_in_wait_list::cl.CL_uint,
                              event_wait_list::Ptr{cl.CL_event},
                              events::Ptr{cl.CL_event})

    @eval @blasfun2 $func(order::clblasOrder, uplo::clblasUplo,
                               N::Csize_t,
                               alpha::$typ,
                               A::cl.CL_mem, offa::Csize_t, lda::Csize_t,
                               X::cl.CL_mem, offx::Csize_t, incx::Cint,
                               beta::$typ,
                               Y::cl.CL_mem, offy::Csize_t, incy::Cint)

end


text:
# Automatically generated using Clang.jl wrap_c, version 0.0.0

using Compat

const SUNDIALS_PACKAGE_VERSION = "2.5.0"
const SUNDIALS_DOUBLE_PRECISION = Cint(1)
const SUNDIALS_BLAS_LAPACK = Cint(0)
const SUNDIALS_DENSE = Cint(1)
const SUNDIALS_BAND = Cint(2)

# Skipping MacroDefinition: DENSE_COL ( A , j ) ( ( A -> cols ) [ j ] )
# Skipping MacroDefinition: DENSE_ELEM ( A , i , j ) ( ( A -> cols ) [ j ] [ i ] )
# Skipping MacroDefinition: BAND_COL ( A , j ) ( ( ( A -> cols ) [ j ] ) + ( A -> s_mu ) )
# Skipping MacroDefinition: BAND_COL_ELEM ( col_j , i , j ) ( col_j [ ( i ) - ( j ) ] )
# Skipping MacroDefinition: BAND_ELEM ( A , i , j ) ( ( A -> cols ) [ j ] [ ( i ) - ( j ) + ( A -> s_mu ) ] )

typealias realtype Cdouble

type _DlsMat
    _type::Cint
    M::Clong
    N::Clong
    ldim::Clong
    mu::Clong
    ml::Clong
    s_mu::Clong
    data::Ptr{realtype}
    ldata::Clong
    cols::Ptr{Ptr{realtype}}
end

typealias DlsMat Ptr{_DlsMat}

# begin enum ANONYMOUS_1
typealias ANONYMOUS_1 UInt32
const PREC_NONE = (UInt32)(0)
const PREC_LEFT = (UInt32)(1)
const PREC_RIGHT = (UInt32)(2)
const PREC_BOTH = (UInt32)(3)
# end enum ANONYMOUS_1

# begin enum ANONYMOUS_2
typealias ANONYMOUS_2 UInt32
const MODIFIED_GS = (UInt32)(1)
const CLASSICAL_GS = (UInt32)(2)
# end enum ANONYMOUS_2

typealias ATimesFn Ptr{Void}
typealias PSolveFn Ptr{Void}

type _generic_N_Vector_Ops
    nvclone::Ptr{Void}
    nvcloneempty::Ptr{Void}
    nvdestroy::Ptr{Void}
    nvspace::Ptr{Void}
    nvgetarraypointer::Ptr{Void}
    nvsetarraypointer::Ptr{Void}
    nvlinearsum::Ptr{Void}
    nvconst::Ptr{Void}
    nvprod::Ptr{Void}
    nvdiv::Ptr{Void}
    nvscale::Ptr{Void}
    nvabs::Ptr{Void}
    nvinv::Ptr{Void}
    nvaddconst::Ptr{Void}
    nvdotprod::Ptr{Void}
    nvmaxnorm::Ptr{Void}
    nvwrmsnorm::Ptr{Void}
    nvwrmsnormmask::Ptr{Void}
    nvmin::Ptr{Void}
    nvwl2norm::Ptr{Void}
    nvl1norm::Ptr{Void}
    nvcompare::Ptr{Void}
    nvinvtest::Ptr{Void}
    nvconstrmask::Ptr{Void}
    nvminquotient::Ptr{Void}
end

typealias N_Vector_Ops Ptr{_generic_N_Vector_Ops}

type _generic_N_Vector
    content::Ptr{Void}
    ops::Ptr{_generic_N_Vector_Ops}
end

typealias N_Vector Ptr{_generic_N_Vector}
typealias N_Vector_S Ptr{N_Vector}

const SPBCG_SUCCESS = Cint(0)
const SPBCG_RES_REDUCED = Cint(1)
const SPBCG_CONV_FAIL = Cint(2)
const SPBCG_PSOLVE_FAIL_REC = Cint(3)
const SPBCG_ATIMES_FAIL_REC = Cint(4)
const SPBCG_PSET_FAIL_REC = Cint(5)
const SPBCG_MEM_NULL = Cint(-1)
const SPBCG_ATIMES_FAIL_UNREC = Cint(-2)
const SPBCG_PSOLVE_FAIL_UNREC = Cint(-3)
const SPBCG_PSET_FAIL_UNREC = Cint(-4)

# Skipping MacroDefinition: SPBCG_VTEMP ( mem ) ( mem -> r )

type SpbcgMemRec
    l_max::Cint
    r_star::N_Vector
    r::N_Vector
    p::N_Vector
    q::N_Vector
    u::N_Vector
    Ap::N_Vector
    vtemp::N_Vector
end

typealias SpbcgMem Ptr{Void}

const SPGMR_SUCCESS = Cint(0)
const SPGMR_RES_REDUCED = Cint(1)
const SPGMR_CONV_FAIL = Cint(2)
const SPGMR_QRFACT_FAIL = Cint(3)
const SPGMR_PSOLVE_FAIL_REC = Cint(4)
const SPGMR_ATIMES_FAIL_REC = Cint(5)
const SPGMR_PSET_FAIL_REC = Cint(6)
const SPGMR_MEM_NULL = Cint(-1)
const SPGMR_ATIMES_FAIL_UNREC = Cint(-2)
const SPGMR_PSOLVE_FAIL_UNREC = Cint(-3)
const SPGMR_GS_FAIL = Cint(-4)
const SPGMR_QRSOL_FAIL = Cint(-5)
const SPGMR_PSET_FAIL_UNREC = Cint(-6)

# Skipping MacroDefinition: SPGMR_VTEMP ( mem ) ( mem -> vtemp )

type _SpgmrMemRec
    l_max::Cint
    V::Ptr{N_Vector}
    Hes::Ptr{Ptr{realtype}}
    givens::Ptr{realtype}
    xcor::N_Vector
    yg::Ptr{realtype}
    vtemp::N_Vector
end

type SpgmrMemRec
    l_max::Cint
    V::Ptr{N_Vector}
    Hes::Ptr{Ptr{realtype}}
    givens::Ptr{realtype}
    xcor::N_Vector
    yg::Ptr{realtype}
    vtemp::N_Vector
end

typealias SpgmrMem Ptr{_SpgmrMemRec}

const SPTFQMR_SUCCESS = Cint(0)
const SPTFQMR_RES_REDUCED = Cint(1)
const SPTFQMR_CONV_FAIL = Cint(2)
const SPTFQMR_PSOLVE_FAIL_REC = Cint(3)
const SPTFQMR_ATIMES_FAIL_REC = Cint(4)
const SPTFQMR_PSET_FAIL_REC = Cint(5)
const SPTFQMR_MEM_NULL = Cint(-1)
const SPTFQMR_ATIMES_FAIL_UNREC = Cint(-2)
const SPTFQMR_PSOLVE_FAIL_UNREC = Cint(-3)
const SPTFQMR_PSET_FAIL_UNREC = Cint(-4)

# Skipping MacroDefinition: SPTFQMR_VTEMP ( mem ) ( mem -> vtemp1 )

type SptfqmrMemRec
    l_max::Cint
    r_star::N_Vector
    q::N_Vector
    d::N_Vector
    v::N_Vector
    p::N_Vector
    r::Ptr{N_Vector}
    u::N_Vector
    vtemp1::N_Vector
    vtemp2::N_Vector
    vtemp3::N_Vector
end

typealias SptfqmrMem Ptr{Void}

# Skipping MacroDefinition: RCONST ( x ) x

const BIG_REAL = DBL_MAX
const SMALL_REAL = DBL_MIN
const UNIT_ROUNDOFF = DBL_EPSILON
const FALSE = Cint(0)
const TRUE = Cint(1)
const CV_ADAMS = Cint(1)
const CV_BDF = Cint(2)
const CV_FUNCTIONAL = Cint(1)
const CV_NEWTON = Cint(2)
const CV_NORMAL = Cint(1)
const CV_ONE_STEP = Cint(2)
const CV_SUCCESS = Cint(0)
const CV_TSTOP_RETURN = Cint(1)
const CV_ROOT_RETURN = Cint(2)
const CV_WARNING = Cint(99)
const CV_TOO_MUCH_WORK = Cint(-1)
const CV_TOO_MUCH_ACC = Cint(-2)
const CV_ERR_FAILURE = Cint(-3)
const CV_CONV_FAILURE = Cint(-4)
const CV_LINIT_FAIL = Cint(-5)
const CV_LSETUP_FAIL = Cint(-6)
const CV_LSOLVE_FAIL = Cint(-7)
const CV_RHSFUNC_FAIL = Cint(-8)
const CV_FIRST_RHSFUNC_ERR = Cint(-9)
const CV_REPTD_RHSFUNC_ERR = Cint(-10)
const CV_UNREC_RHSFUNC_ERR = Cint(-11)
const CV_RTFUNC_FAIL = Cint(-12)
const CV_MEM_FAIL = Cint(-20)
const CV_MEM_NULL = Cint(-21)
const CV_ILL_INPUT = Cint(-22)
const CV_NO_MALLOC = Cint(-23)
const CV_BAD_K = Cint(-24)
const CV_BAD_T = Cint(-25)
const CV_BAD_DKY = Cint(-26)
const CV_TOO_CLOSE = Cint(-27)

typealias CVRhsFn Ptr{Void}

CVRhsFn_wrapper(fp::CVRhsFn) = fp
CVRhsFn_wrapper(f) = cfunction(f,Cint,(realtype,N_Vector,N_Vector,Ptr{Void}))

typealias CVRootFn Ptr{Void}

CVRootFn_wrapper(fp::CVRootFn) = fp
CVRootFn_wrapper(f) = cfunction(f,Cint,(realtype,N_Vector,Ptr{realtype},Ptr{Void}))

typealias CVEwtFn Ptr{Void}
typealias CVErrHandlerFn Ptr{Void}

const CVDLS_SUCCESS = Cint(0)
const CVDLS_MEM_NULL = Cint(-1)
const CVDLS_LMEM_NULL = Cint(-2)
const CVDLS_ILL_INPUT = Cint(-3)
const CVDLS_MEM_FAIL = Cint(-4)
const CVDLS_JACFUNC_UNRECVR = Cint(-5)
const CVDLS_JACFUNC_RECVR = Cint(-6)

typealias CVDlsDenseJacFn Ptr{Void}
typealias CVDlsBandJacFn Ptr{Void}

const CVSPILS_SUCCESS = Cint(0)
const CVSPILS_MEM_NULL = Cint(-1)
const CVSPILS_LMEM_NULL = Cint(-2)
const CVSPILS_ILL_INPUT = Cint(-3)
const CVSPILS_MEM_FAIL = Cint(-4)
const CVSPILS_PMEM_NULL = Cint(-5)
const CVSPILS_MAXL = Cint(5)
const CVSPILS_MSBPRE = Cint(50)

# Skipping MacroDefinition: CVSPILS_DGMAX RCONST ( 0.2 )
# Skipping MacroDefinition: CVSPILS_EPLIN RCONST ( 0.05 )

typealias CVSpilsPrecSetupFn Ptr{Void}
typealias CVSpilsPrecSolveFn Ptr{Void}
typealias CVSpilsJacTimesVecFn Ptr{Void}

const CV_SIMULTANEOUS = Cint(1)
const CV_STAGGERED = Cint(2)
const CV_STAGGERED1 = Cint(3)
const CV_CENTERED = Cint(1)
const CV_FORWARD = Cint(2)
const CV_HERMITE = Cint(1)
const CV_POLYNOMIAL = Cint(2)
const CV_NO_QUAD = Cint(-30)
const CV_QRHSFUNC_FAIL = Cint(-31)
const CV_FIRST_QRHSFUNC_ERR = Cint(-32)
const CV_REPTD_QRHSFUNC_ERR = Cint(-33)
const CV_UNREC_QRHSFUNC_ERR = Cint(-34)
const CV_NO_SENS = Cint(-40)
const CV_SRHSFUNC_FAIL = Cint(-41)
const CV_FIRST_SRHSFUNC_ERR = Cint(-42)
const CV_REPTD_SRHSFUNC_ERR = Cint(-43)
const CV_UNREC_SRHSFUNC_ERR = Cint(-44)
const CV_BAD_IS = Cint(-45)
const CV_NO_QUADSENS = Cint(-50)
const CV_QSRHSFUNC_FAIL = Cint(-51)
const CV_FIRST_QSRHSFUNC_ERR = Cint(-52)
const CV_REPTD_QSRHSFUNC_ERR = Cint(-53)
const CV_UNREC_QSRHSFUNC_ERR = Cint(-54)
const CV_NO_ADJ = Cint(-101)
const CV_NO_FWD = Cint(-102)
const CV_NO_BCK = Cint(-103)
const CV_BAD_TB0 = Cint(-104)
const CV_REIFWD_FAIL = Cint(-105)
const CV_FWD_FAIL = Cint(-106)
const CV_GETY_BADT = Cint(-107)

typealias CVQuadRhsFn Ptr{Void}
typealias CVSensRhsFn Ptr{Void}
typealias CVSensRhs1Fn Ptr{Void}
typealias CVQuadSensRhsFn Ptr{Void}
typealias CVRhsFnB Ptr{Void}
typealias CVRhsFnBS Ptr{Void}
typealias CVQuadRhsFnB Ptr{Void}
typealias CVQuadRhsFnBS Ptr{Void}

type CVadjCheckPointRec
    my_addr::Ptr{Void}
    next_addr::Ptr{Void}
    t0::realtype
    t1::realtype
    nstep::Clong
    order::Cint
    step::realtype
end

const CVDLS_NO_ADJ = Cint(-101)
const CVDLS_LMEMB_NULL = Cint(-102)

typealias CVDlsDenseJacFnB Ptr{Void}
typealias CVDlsBandJacFnB Ptr{Void}

const CVSPILS_NO_ADJ = Cint(-101)
const CVSPILS_LMEMB_NULL = Cint(-102)

typealias CVSpilsPrecSetupFnB Ptr{Void}
typealias CVSpilsPrecSolveFnB Ptr{Void}
typealias CVSpilsJacTimesVecFnB Ptr{Void}

const IDA_NORMAL = Cint(1)
const IDA_ONE_STEP = Cint(2)
const IDA_YA_YDP_INIT = Cint(1)
const IDA_Y_INIT = Cint(2)
const IDA_SUCCESS = Cint(0)
const IDA_TSTOP_RETURN = Cint(1)
const IDA_ROOT_RETURN = Cint(2)
const IDA_WARNING = Cint(99)
const IDA_TOO_MUCH_WORK = Cint(-1)
const IDA_TOO_MUCH_ACC = Cint(-2)
const IDA_ERR_FAIL = Cint(-3)
const IDA_CONV_FAIL = Cint(-4)
const IDA_LINIT_FAIL = Cint(-5)
const IDA_LSETUP_FAIL = Cint(-6)
const IDA_LSOLVE_FAIL = Cint(-7)
const IDA_RES_FAIL = Cint(-8)
const IDA_REP_RES_ERR = Cint(-9)
const IDA_RTFUNC_FAIL = Cint(-10)
const IDA_CONSTR_FAIL = Cint(-11)
const IDA_FIRST_RES_FAIL = Cint(-12)
const IDA_LINESEARCH_FAIL = Cint(-13)
const IDA_NO_RECOVERY = Cint(-14)
const IDA_MEM_NULL = Cint(-20)
const IDA_MEM_FAIL = Cint(-21)
const IDA_ILL_INPUT = Cint(-22)
const IDA_NO_MALLOC = Cint(-23)
const IDA_BAD_EWT = Cint(-24)
const IDA_BAD_K = Cint(-25)
const IDA_BAD_T = Cint(-26)
const IDA_BAD_DKY = Cint(-27)

typealias IDAResFn Ptr{Void}

IDAResFn_wrapper(fp::IDAResFn) = fp
IDAResFn_wrapper(f) = cfunction(f,Cint,(realtype,N_Vector,N_Vector,N_Vector,Ptr{Void}))

typealias IDARootFn Ptr{Void}

IDARootFn_wrapper(fp::IDARootFn) = fp
IDARootFn_wrapper(f) = cfunction(f,Cint,(realtype,N_Vector,N_Vector,Ptr{realtype},Ptr{Void}))

typealias IDAEwtFn Ptr{Void}
typealias IDAErrHandlerFn Ptr{Void}

const IDADLS_SUCCESS = Cint(0)
const IDADLS_MEM_NULL = Cint(-1)
const IDADLS_LMEM_NULL = Cint(-2)
const IDADLS_ILL_INPUT = Cint(-3)
const IDADLS_MEM_FAIL = Cint(-4)
const IDADLS_JACFUNC_UNRECVR = Cint(-5)
const IDADLS_JACFUNC_RECVR = Cint(-6)

typealias IDADlsDenseJacFn Ptr{Void}
typealias IDADlsBandJacFn Ptr{Void}

const IDASPILS_SUCCESS = Cint(0)
const IDASPILS_MEM_NULL = Cint(-1)
const IDASPILS_LMEM_NULL = Cint(-2)
const IDASPILS_ILL_INPUT = Cint(-3)
const IDASPILS_MEM_FAIL = Cint(-4)
const IDASPILS_PMEM_NULL = Cint(-5)

typealias IDASpilsPrecSetupFn Ptr{Void}
typealias IDASpilsPrecSolveFn Ptr{Void}
typealias IDASpilsJacTimesVecFn Ptr{Void}

const IDA_SIMULTANEOUS = Cint(1)
const IDA_STAGGERED = Cint(2)
const IDA_CENTERED = Cint(1)
const IDA_FORWARD = Cint(2)
const IDA_HERMITE = Cint(1)
const IDA_POLYNOMIAL = Cint(2)
const IDA_NO_QUAD = Cint(-30)
const IDA_QRHS_FAIL = Cint(-31)
const IDA_FIRST_QRHS_ERR = Cint(-32)
const IDA_REP_QRHS_ERR = Cint(-33)
const IDA_NO_SENS = Cint(-40)
const IDA_SRES_FAIL = Cint(-41)
const IDA_REP_SRES_ERR = Cint(-42)
const IDA_BAD_IS = Cint(-43)
const IDA_NO_QUADSENS = Cint(-50)
const IDA_QSRHS_FAIL = Cint(-51)
const IDA_FIRST_QSRHS_ERR = Cint(-52)
const IDA_REP_QSRHS_ERR = Cint(-53)
const IDA_NO_ADJ = Cint(-101)
const IDA_NO_FWD = Cint(-102)
const IDA_NO_BCK = Cint(-103)
const IDA_BAD_TB0 = Cint(-104)
const IDA_REIFWD_FAIL = Cint(-105)
const IDA_FWD_FAIL = Cint(-106)
const IDA_GETY_BADT = Cint(-107)

typealias IDAQuadRhsFn Ptr{Void}
typealias IDASensResFn Ptr{Void}
typealias IDAQuadSensRhsFn Ptr{Void}
typealias IDAResFnB Ptr{Void}
typealias IDAResFnBS Ptr{Void}
typealias IDAQuadRhsFnB Ptr{Void}
typealias IDAQuadRhsFnBS Ptr{Void}

type IDAadjCheckPointRec
    my_addr::Ptr{Void}
    next_addr::Ptr{Void}
    t0::realtype
    t1::realtype
    nstep::Clong
    order::Cint
    step::realtype
end

const IDADLS_NO_ADJ = Cint(-101)
const IDADLS_LMEMB_NULL = Cint(-102)

typealias IDADlsDenseJacFnB Ptr{Void}
typealias IDADlsBandJacFnB Ptr{Void}

const IDASPILS_NO_ADJ = Cint(-101)
const IDASPILS_LMEMB_NULL = Cint(-102)

typealias IDASpilsPrecSetupFnB Ptr{Void}
typealias IDASpilsPrecSolveFnB Ptr{Void}
typealias IDASpilsJacTimesVecFnB Ptr{Void}

const KIN_SUCCESS = Cint(0)
const KIN_INITIAL_GUESS_OK = Cint(1)
const KIN_STEP_LT_STPTOL = Cint(2)
const KIN_WARNING = Cint(99)
const KIN_MEM_NULL = Cint(-1)
const KIN_ILL_INPUT = Cint(-2)
const KIN_NO_MALLOC = Cint(-3)
const KIN_MEM_FAIL = Cint(-4)
const KIN_LINESEARCH_NONCONV = Cint(-5)
const KIN_MAXITER_REACHED = Cint(-6)
const KIN_MXNEWT_5X_EXCEEDED = Cint(-7)
const KIN_LINESEARCH_BCFAIL = Cint(-8)
const KIN_LINSOLV_NO_RECOVERY = Cint(-9)
const KIN_LINIT_FAIL = Cint(-10)
const KIN_LSETUP_FAIL = Cint(-11)
const KIN_LSOLVE_FAIL = Cint(-12)
const KIN_SYSFUNC_FAIL = Cint(-13)
const KIN_FIRST_SYSFUNC_ERR = Cint(-14)
const KIN_REPTD_SYSFUNC_ERR = Cint(-15)
const KIN_ETACHOICE1 = Cint(1)
const KIN_ETACHOICE2 = Cint(2)
const KIN_ETACONSTANT = Cint(3)
const KIN_NONE = Cint(0)
const KIN_LINESEARCH = Cint(1)

typealias KINSysFn Ptr{Void}

KINSysFn_wrapper(fp::KINSysFn) = fp
KINSysFn_wrapper(f) = cfunction(f,Cint,(N_Vector,N_Vector,Ptr{Void}))

typealias KINErrHandlerFn Ptr{Void}
typealias KINInfoHandlerFn Ptr{Void}

const KINDLS_SUCCESS = Cint(0)
const KINDLS_MEM_NULL = Cint(-1)
const KINDLS_LMEM_NULL = Cint(-2)
const KINDLS_ILL_INPUT = Cint(-3)
const KINDLS_MEM_FAIL = Cint(-4)
const KINDLS_JACFUNC_UNRECVR = Cint(-5)
const KINDLS_JACFUNC_RECVR = Cint(-6)

typealias KINDlsDenseJacFn Ptr{Void}
typealias KINDlsBandJacFn Ptr{Void}

const KINSPILS_SUCCESS = Cint(0)
const KINSPILS_MEM_NULL = Cint(-1)
const KINSPILS_LMEM_NULL = Cint(-2)
const KINSPILS_ILL_INPUT = Cint(-3)
const KINSPILS_MEM_FAIL = Cint(-4)
const KINSPILS_PMEM_NULL = Cint(-5)
const KINSPILS_MAXL = Cint(10)

typealias KINSpilsPrecSetupFn Ptr{Void}
typealias KINSpilsPrecSolveFn Ptr{Void}
typealias KINSpilsJacTimesVecFn Ptr{Void}

# Skipping MacroDefinition: NV_CONTENT_S ( v ) ( ( N_VectorContent_Serial ) ( v -> content ) )
# Skipping MacroDefinition: NV_LENGTH_S ( v ) ( NV_CONTENT_S ( v ) -> length )
# Skipping MacroDefinition: NV_OWN_DATA_S ( v ) ( NV_CONTENT_S ( v ) -> own_data )
# Skipping MacroDefinition: NV_DATA_S ( v ) ( NV_CONTENT_S ( v ) -> data )
# Skipping MacroDefinition: NV_Ith_S ( v , i ) ( NV_DATA_S ( v ) [ i ] )

type _N_VectorContent_Serial
    length::Clong
    own_data::Cint
    data::Ptr{realtype}
end

typealias N_VectorContent_Serial Ptr{_N_VectorContent_Serial}

const FCMIX_CVODE = Cint(1)
const FCMIX_IDA = Cint(2)
const FCMIX_KINSOL = Cint(3)

# Skipping MacroDefinition: MIN ( A , B ) ( ( A ) < ( B ) ? ( A ) : ( B ) )
# Skipping MacroDefinition: MAX ( A , B ) ( ( A ) > ( B ) ? ( A ) : ( B ) )
# Skipping MacroDefinition: SQR ( A ) ( ( A ) * ( A ) )

typealias CVLocalFn Ptr{Void}
typealias CVCommFn Ptr{Void}

const CVDIAG_SUCCESS = Cint(0)
const CVDIAG_MEM_NULL = Cint(-1)
const CVDIAG_LMEM_NULL = Cint(-2)
const CVDIAG_ILL_INPUT = Cint(-3)
const CVDIAG_MEM_FAIL = Cint(-4)
const CVDIAG_INV_FAIL = Cint(-5)
const CVDIAG_RHSFUNC_UNRECVR = Cint(-6)
const CVDIAG_RHSFUNC_RECVR = Cint(-7)

typealias CVLocalFnB Ptr{Void}
typealias CVCommFnB Ptr{Void}

const CVDIAG_NO_ADJ = Cint(-101)

typealias IDABBDLocalFn Ptr{Void}
typealias IDABBDCommFn Ptr{Void}
typealias IDABBDLocalFnB Ptr{Void}
typealias IDABBDCommFnB Ptr{Void}

const KINBBDPRE_SUCCESS = Cint(0)
const KINBBDPRE_PDATA_NULL = Cint(-11)
const KINBBDPRE_FUNC_UNRECVR = Cint(-12)

typealias KINCommFn Ptr{Void}
typealias KINLocalFn Ptr{Void}


text:
inpt = """NOT dq -> dr
kg OR kf -> kh
ep OR eo -> eq
3176 -> b
NOT gs -> gt
dd OR do -> dp
eg AND ei -> ej
y AND ae -> ag
jx AND jz -> ka
lf RSHIFT 2 -> lg
z AND aa -> ac
dy AND ej -> el
bj OR bi -> bk
kk RSHIFT 3 -> km
NOT cn -> co
gn AND gp -> gq
cq AND cs -> ct
eo LSHIFT 15 -> es
lg OR lm -> ln
dy OR ej -> ek
NOT di -> dj
1 AND fi -> fj
kf LSHIFT 15 -> kj
NOT jy -> jz
NOT ft -> fu
fs AND fu -> fv
NOT hr -> hs
ck OR cl -> cm
jp RSHIFT 5 -> js
iv OR jb -> jc
is OR it -> iu
ld OR le -> lf
NOT fc -> fd
NOT dm -> dn
bn OR by -> bz
aj AND al -> am
cd LSHIFT 15 -> ch
jp AND ka -> kc
ci OR ct -> cu
gv AND gx -> gy
de AND dk -> dm
x RSHIFT 5 -> aa
et RSHIFT 2 -> eu
x RSHIFT 1 -> aq
ia OR ig -> ih
bk LSHIFT 1 -> ce
y OR ae -> af
NOT ca -> cb
e AND f -> h
ia AND ig -> ii
ck AND cl -> cn
NOT jh -> ji
z OR aa -> ab
1 AND en -> eo
ib AND ic -> ie
NOT eh -> ei
iy AND ja -> jb
NOT bb -> bc
ha OR gz -> hb
1 AND cx -> cy
NOT ax -> ay
ev OR ew -> ex
bn RSHIFT 2 -> bo
er OR es -> et
eu OR fa -> fb
jp OR ka -> kb
ea AND eb -> ed
k AND m -> n
et RSHIFT 3 -> ev
et RSHIFT 5 -> ew
hz RSHIFT 1 -> is
ki OR kj -> kk
NOT h -> i
lv LSHIFT 15 -> lz
as RSHIFT 1 -> bl
hu LSHIFT 15 -> hy
iw AND ix -> iz
lf RSHIFT 1 -> ly
fp OR fv -> fw
1 AND am -> an
ap LSHIFT 1 -> bj
u LSHIFT 1 -> ao
b RSHIFT 5 -> f
jq AND jw -> jy
iu RSHIFT 3 -> iw
ih AND ij -> ik
NOT iz -> ja
de OR dk -> dl
iu OR jf -> jg
as AND bd -> bf
b RSHIFT 3 -> e
jq OR jw -> jx
iv AND jb -> jd
cg OR ch -> ci
iu AND jf -> jh
lx -> a
1 AND cc -> cd
ly OR lz -> ma
NOT el -> em
1 AND bh -> bi
fb AND fd -> fe
lf OR lq -> lr
bn RSHIFT 3 -> bp
bn AND by -> ca
af AND ah -> ai
cf LSHIFT 1 -> cz
dw OR dx -> dy
gj AND gu -> gw
jg AND ji -> jj
jr OR js -> jt
bl OR bm -> bn
gj RSHIFT 2 -> gk
cj OR cp -> cq
gj OR gu -> gv
b OR n -> o
o AND q -> r
bi LSHIFT 15 -> bm
dy RSHIFT 1 -> er
cu AND cw -> cx
iw OR ix -> iy
hc OR hd -> he
0 -> c
db OR dc -> dd
kk RSHIFT 2 -> kl
eq LSHIFT 1 -> fk
dz OR ef -> eg
NOT ed -> ee
lw OR lv -> lx
fw AND fy -> fz
dz AND ef -> eh
jp RSHIFT 3 -> jr
lg AND lm -> lo
ci RSHIFT 2 -> cj
be AND bg -> bh
lc LSHIFT 1 -> lw
hm AND ho -> hp
jr AND js -> ju
1 AND io -> ip
cm AND co -> cp
ib OR ic -> id
NOT bf -> bg
fo RSHIFT 5 -> fr
ip LSHIFT 15 -> it
jt AND jv -> jw
jc AND je -> jf
du OR dt -> dv
NOT fx -> fy
aw AND ay -> az
ge LSHIFT 15 -> gi
NOT ak -> al
fm OR fn -> fo
ff AND fh -> fi
ci RSHIFT 5 -> cl
cz OR cy -> da
NOT ey -> ez
NOT ju -> jv
NOT ls -> lt
kk AND kv -> kx
NOT ii -> ij
kl AND kr -> kt
jk LSHIFT 15 -> jo
e OR f -> g
NOT bs -> bt
hi AND hk -> hl
hz OR ik -> il
ek AND em -> en
ao OR an -> ap
dv LSHIFT 1 -> ep
an LSHIFT 15 -> ar
fo RSHIFT 1 -> gh
NOT im -> in
kk RSHIFT 1 -> ld
hw LSHIFT 1 -> iq
ec AND ee -> ef
hb LSHIFT 1 -> hv
kb AND kd -> ke
x AND ai -> ak
dd AND do -> dq
aq OR ar -> as
iq OR ip -> ir
dl AND dn -> do
iu RSHIFT 5 -> ix
as OR bd -> be
NOT go -> gp
fk OR fj -> fl
jm LSHIFT 1 -> kg
NOT cv -> cw
dp AND dr -> ds
dt LSHIFT 15 -> dx
et RSHIFT 1 -> fm
dy RSHIFT 3 -> ea
fp AND fv -> fx
NOT p -> q
dd RSHIFT 2 -> de
eu AND fa -> fc
ba AND bc -> bd
dh AND dj -> dk
lr AND lt -> lu
he RSHIFT 1 -> hx
ex AND ez -> fa
df OR dg -> dh
fj LSHIFT 15 -> fn
NOT kx -> ky
gk OR gq -> gr
dy RSHIFT 2 -> dz
gh OR gi -> gj
lj AND ll -> lm
x OR ai -> aj
bz AND cb -> cc
1 AND lu -> lv
as RSHIFT 3 -> au
ce OR cd -> cf
il AND in -> io
dd RSHIFT 1 -> dw
NOT lo -> lp
c LSHIFT 1 -> t
dd RSHIFT 3 -> df
dd RSHIFT 5 -> dg
lh AND li -> lk
lf RSHIFT 5 -> li
dy RSHIFT 5 -> eb
NOT kt -> ku
at OR az -> ba
x RSHIFT 3 -> z
NOT lk -> ll
lb OR la -> lc
1 AND r -> s
lh OR li -> lj
ln AND lp -> lq
kk RSHIFT 5 -> kn
ea OR eb -> ec
ci AND ct -> cv
b RSHIFT 2 -> d
jp RSHIFT 1 -> ki
NOT cr -> cs
NOT jd -> je
jp RSHIFT 2 -> jq
jn OR jo -> jp
lf RSHIFT 3 -> lh
1 AND ds -> dt
lf AND lq -> ls
la LSHIFT 15 -> le
NOT fg -> fh
at AND az -> bb
au AND av -> ax
kw AND ky -> kz
v OR w -> x
kk OR kv -> kw
ks AND ku -> kv
kh LSHIFT 1 -> lb
1 AND kz -> la
NOT kc -> kd
x RSHIFT 2 -> y
et OR fe -> ff
et AND fe -> fg
NOT ac -> ad
jl OR jk -> jm
1 AND jj -> jk
bn RSHIFT 1 -> cg
NOT kp -> kq
ci RSHIFT 3 -> ck
ev AND ew -> ey
1 AND ke -> kf
cj AND cp -> cr
ir LSHIFT 1 -> jl
NOT gw -> gx
as RSHIFT 2 -> at
iu RSHIFT 1 -> jn
cy LSHIFT 15 -> dc
hg OR hh -> hi
ci RSHIFT 1 -> db
au OR av -> aw
km AND kn -> kp
gj RSHIFT 1 -> hc
iu RSHIFT 2 -> iv
ab AND ad -> ae
da LSHIFT 1 -> du
NOT bw -> bx
km OR kn -> ko
ko AND kq -> kr
bv AND bx -> by
kl OR kr -> ks
1 AND ht -> hu
df AND dg -> di
NOT ag -> ah
d OR j -> k
d AND j -> l
b AND n -> p
gf OR ge -> gg
gg LSHIFT 1 -> ha
bn RSHIFT 5 -> bq
bo OR bu -> bv
1 AND gy -> gz
s LSHIFT 15 -> w
NOT ie -> if
as RSHIFT 5 -> av
bo AND bu -> bw
hz AND ik -> im
bp AND bq -> bs
b RSHIFT 1 -> v
NOT l -> m
bp OR bq -> br
g AND i -> j
br AND bt -> bu
t OR s -> u
hz RSHIFT 5 -> ic
gk AND gq -> gs
fl LSHIFT 1 -> gf
he RSHIFT 3 -> hg
gz LSHIFT 15 -> hd
hf OR hl -> hm
1 AND gd -> ge
fo OR fz -> ga
id AND if -> ig
fo AND fz -> gb
gr AND gt -> gu
he OR hp -> hq
fq AND fr -> ft
ga AND gc -> gd
fo RSHIFT 2 -> fp
gl OR gm -> gn
hg AND hh -> hj
NOT hn -> ho
gl AND gm -> go
he RSHIFT 5 -> hh
NOT gb -> gc
hq AND hs -> ht
hz RSHIFT 3 -> ib
hz RSHIFT 2 -> ia
fq OR fr -> fs
hx OR hy -> hz
he AND hp -> hr
gj RSHIFT 5 -> gm
hf AND hl -> hn
hv OR hu -> hw
NOT hj -> hk
gj RSHIFT 3 -> gl
fo RSHIFT 3 -> fq
he RSHIFT 2 -> hf"""

function buildCode(sentence)
	tokens = split(sentence)
	expression = uppercase(tokens[end])
	expression *= "xxx = "
	for x in 1:length(tokens)-2
		tokens[x] = uppercase(tokens[x])
		if tokens[x] == "AND"
			tokens[x] = "&"
		elseif tokens[x] == "OR"
			tokens[x] = "|"
		elseif tokens[x] == "NOT"
			tokens[x] = "~"
		elseif tokens[x] == "LSHIFT"
			tokens[x] = "<<"
		elseif tokens[x] == "RSHIFT"
			tokens[x] = ">>"
		elseif !(typeof(parse(tokens[x])) <: Number)
			tokens[x] = tokens[x] * "xxx" # to protect variable names from clashes with builtins
		end
	end
	expression *= join(tokens[1:end-2], " ")
	println(expression)
	return parse(expression)
end

grid = Expr[]
for sen in split(inpt, "\n")
	sent = buildCode(sen)
	push!(grid, sent)
end

function isDefined(expr)
	if typeof(expr) == Symbol
		# println(expr)
		if !(expr in names(Main) || expr in names(Base))
			return false
		end
	elseif typeof(expr) == Expr
		for a in expr.args
			if !isDefined(a)
				return false
			end
		end
	# else
		# println(typeof(expr), expr)
	end
	return true
end

while length(grid) > 0
	for (i,e) in enumerate(grid)
		if isDefined(grid[i].args[2])
			p = splice!(grid, i)
			println(p)
			eval(p)
			break
		end
	end
end
println(Axxx)

max_stars_repo_name: jstrube/adventofcode

max_issues_repo_name: jstrube/adventofcode

max_forks_repo_name: jstrube/adventofcode


text:
"""Example showing how to use SALib (via PyCall) to generate, run, and analyze results"""

using Pkg

ENV["PYTHON"] = raw"C:\programs\miniconda3\envs\salib"  # Path to environment python.exe
ENV["CONDA_JL_HOME"] = raw"C:\programs\miniconda3\envs\salib"  # Path to conda environment folder

# Pkg.build("Conda")
Pkg.build("PyCall")

using PyCall
using Agtor
using Dates
using DataFrames

# Start month/day of growing season
# Represents when water allocations start getting announced.
const gs_start = (5, 1)

function setup_basin(data_dir, climate_data)
    basin_spec_dir::String = "$(data_dir)basins"
    basin_spec::Dict = load_spec(basin_spec_dir)[:TestBasin]

    basin_name = basin_spec[:name]
    zone_specs = basin_spec[:zone_spec]

    collated_specs = []
    agparams = collect_agparams!(zone_specs, collated_specs)

    Manager_A = BaseManager("optimizing")
    Manager_B = RigidManager("rigid", 0.05)
    manage_zones = ((Manager_A, ("Zone_1", )), (Manager_B, ("Zone_2", )))

    b = Basin(name=basin_name, zone_spec=zone_specs, climate_data=climate_data, managers=manage_zones)

    return b, agparams
end

function run_model(scenario_id, basin)

    # Example annual water allocations
    allocs = (surface_water=150.0, groundwater=40.0)

    for dt_i in basin.climate.time_steps
        # Resetting allocations for example run
        if monthday(dt_i) == gs_start
            for z in basin.zones
                update_available_water!(z, allocs)
            end
        end

        run_timestep!(basin)
    end

    # Save results as they complete
    all_res = []
    @sync for z in basin.zones
        results = collect_results(z)

        # Uncomment the below if intermediate results are desired
        # run_id = "$(scenario_id)/$(z.name)"
        # @async save_results!("basin_async_save.jld2", run_id, results; mode="a+")

        push!(all_res, results[1])
    end

    return calculate_metrics(all_res)
end

function scenario_run(samples, basin)::Array{Array{Float64}}
    basin_wide_res = []
    for (row_id, r) in enumerate(eachrow(samples))
        tmp_b = deepcopy(basin)
        update_model!(tmp_b, r)
        push!(basin_wide_res, run_model(row_id, tmp_b))
    end

    return basin_wide_res
end

function pd_to_df(df_pd)
    colnames = map(Symbol, df_pd.columns)
    df = DataFrame(Any[Array(df_pd[c].values) for c in colnames], colnames)

    return df
end

function calculate_metrics(collated_results::Union{Tuple, Array})::Array{Float64}

    z_res = DataFrame(fill(Float64, 4), Symbol.((i for i in range(1,4))), length(collated_results))

    for (idx, zone_res) in enumerate(collated_results)
        z_res[idx, :] = calculate_metrics(zone_res)
    end

    return [mean(z_res[:, i]) for i in names(z_res)]
end

data_dir = "test/data/"
example_basin, agparams = setup_basin(data_dir, "test/data/climate/basin_climate_data.csv")
# scen_data = joinpath(data_dir, "scenarios", "sampled_params.csv")
# samples = DataFrame!(CSV.File(scen_data))

param_names = agparams[:, :name]
param_bounds = agparams[:, [:min_val, :max_val]]

py"""
from SALib import ProblemSpec
import numpy as np
import pandas as pd

sp = ProblemSpec({
  'names': $(param_names),
  'bounds': list($(Array(param_bounds))),
  'output': ['CV', 'Mean ML per Yield [ML/t]', 'Mean Income [$]',
             'Mean Water Use [ML]']
})

sp.sample_latin(4, seed=101)

orig_samples = sp.samples

tmp = pd.DataFrame(sp.samples, columns=sp['names'])
# tmp.to_csv("salib_example.csv", index=False)

# Attach Julia DataFrame, bypassing any checks...
sp._samples = $(pd_to_df)(tmp)

# Run model
sp.evaluate($(scenario_run), $(example_basin))

# Revert back to original samples as numpy array
sp.samples = orig_samples

# Have to coerce results to numpy array as these
# are returned as a Python list by Julia
sp.results = np.array(sp.results)

# Perform analysis
sp.analyze_rbd_fast()
"""

@info "Analysis:" py"sp.to_df()"

using PyCall

py"""
from SALib import ProblemSpec
import numpy as np
import pandas as pd

sp = ProblemSpec({
  'names': $(param_names),
  'bounds': list($(Array(param_bounds))),
  'output': ['CV', 'Mean ML per Yield [ML/t]', 'Mean Income [$]',
             'Mean Water Use [ML]']
})

sp.sample_latin(4, seed=101)

orig_samples = sp.samples

# Convert Pandas DF to Julia DataFrame, bypassing any checks...
tmp = pd.DataFrame(sp.samples, columns=sp['names'])
sp._samples = $(pd_to_df)(tmp)

# Run model
sp.evaluate($(scenario_run), $(example_basin))

# Revert back to original samples as numpy array
sp.samples = orig_samples

# Have to coerce results to numpy array as these
# are returned as a Python list by Julia
sp.results = np.array(sp.results)

# Perform analysis
sp.analyze_rbd_fast()
"""

@info "Analysis:" py"sp.to_df()"

max_stars_repo_name: ConnectedSystems/agtor-jl

max_issues_repo_name: ConnectedSystems/agtor-jl

max_forks_repo_name: ConnectedSystems/agtor-jl


text:
#Time scale analysis
using Sundials
using FastGaussQuadrature

function getfractionbelow(histo,n,N)
    s = 0
    for w in histo.weights
        if w < n
            s += w
        else
            s += n
        end
    end
    return s/N
end

function getfractionleft(histo,tol,N)
    s = 0
    for i in 2:length(histo.edges[1])
        if histo.edges[1][i] < tol
            s += histo.weights[i-1]
        end
    end
    return s/N
end

function splitnum(histo,n)
    c = 0
    boo = histo.weights[1] < n
    for (i,x) in enumerate(histo.weights)
        if boo && x > n
            c += 1
            boo = false
        elseif !boo && x < n
            c += 1
            boo = true
        end
    end
    return c
end

function splitnuminds(histo,n)
    c = 0
    inds = Array{Int64,1}()
    boo = histo.weights[1] < n
    for (i,x) in enumerate(histo.weights)
        if boo && x > n
            c += 1
            boo = false
            push!(inds,i)
        elseif !boo && x < n
            c += 1
            boo = true
            push!(inds,i)
        end
    end
    return c,inds
end

"""
Identify the timescale that best divides the fastest relaxing species (usually radicals)
from the slower relaxing species (usually stable species)
This is done by analysis of the timescale histogram
If there are more than two timescale peaks this algorithm will pick
the timescale that separates the two fastest peaks

timescales are grained multiplicatively with taumax,taumin, with taures being the
multiplicative factor
usediag=true tells the algorithm to determine timescales from the diagonal of the jacobian
otherwise it will compute and use the eigenvalues, this hasn't been observed to have
a significant impact on dividing time scale prediction

Because this relies on constructing the timescale distribution of the species
It doesn't work for very small mechanisms with ~12 or less species and
may not be able to find a multimodal slice in these cases
"""
function getdividingtimescale(Jy;taumax=1e4,taumin=1e-18,taures=10.0^0.5,usediag=true)
    if usediag
        taus = 1.0./abs.(diag(Jy))
    else
        taus = 1.0./abs.(eigvals(Jy))
    end
    N = length(taus)
    histo = fit(StatsBase.Histogram, taus, 10.0.^(log10(taumin):log10(taures):log10(taumax)), closed=:left)
    Nmax = maximum(histo.weights)
    splits = zeros(Nmax)
    fracbelow = zeros(Nmax)
    for i = 1:Nmax
        splits[i] = splitnum(histo,i)
        fracbelow[i] = getfractionbelow(histo,i,N)
    end
    fracleft = cumsum(histo.weights)./N
    fourinds = Array{Int64,1}()
    fracleftvals = Array{Float64,1}()
    fracbelowvals = Array{Float64,1}()
    for i = 1:Nmax
        if splits[i] >= 4
            c,inds = splitnuminds(histo,i)
            ind = inds[2]
            fl = fracleft[ind]
            fb = fracbelow[i]
            push!(fourinds,ind)
            push!(fracleftvals,fl)
            push!(fracbelowvals,fb)
        end
    end
    if length(fourinds) == 0
        error("Could not find a multimodal slice")
    end
    ind = argmin((fracleftvals.- 0.5).^2 .+ (fracbelowvals .- 0.5).^2)
    tau = histo.edges[1][fourinds[ind]+1]
    return tau
end

function getdividingtimescale(sim,t;taumax=1e6,taumin=1e-18,taures=10.0^0.5,usediag=true)
    Jy = jacobiany(sim.sol(t),sim.p,t,sim.domain,sim.interfaces);
    getdividingtimescale(Jy,taumax=taumax,taumin=taumin,taures=taures,usediag=usediag)
end
export getdividingtimescale

function jacobiany(sol,t,p)
    Jy = zeros(length(sol.prob.u0),length(sol.prob.u0))
    sol.prob.f.jac(Jy,sol(t),p,t)
    return Jy
end

function jacobianp(sol,t,p)
    Jp = zeros(length(sol.prob.u0),length(p))
    sol.prob.f.paramjac(Jp,sol(t),p,t)
    return Jp
end

function jacobianysparse(sol,t,p)
    Jy = spzeros(length(sol.prob.u0),length(sol.prob.u0))
    sol.prob.f.jac(Jy,sol(t),p,t)
    return Jy
end

function jacobianpsparse(sol,t,p)
    Jp = spzeros(length(sol.prob.u0),length(p))
    sol.prob.f.paramjac(Jp,sol(t),p,t)
    return Jp
end

function jacobianp(sol,t,p,ind)
   return ForwardDiff.jacobian(pv->sol.prob.f(sol(t),hcat(p[1:ind-1],[pv],p[ind+1:end]),t),p[ind])
end

function normalizefulltransitorysensitivities!(dSdt,sim::Simulation,t)
    y = sim.sol(t)
    @views ns = y[sim.domain.indexes[1]:sim.domain.indexes[2]]
    dSdt .*= sim.p'
    @views dSdt[sim.domain.indexes[1]:sim.domain.indexes[2],:] ./= ns
    return dSdt
end

function normalizefulltransitorysensitivities!(dSdt,ssys::SystemSimulation,t)
    y = ssys.sol(t)
    dSdt .*= ssys.p'
    for sim in ssys.sims
        @views ns = y[sim.domain.indexes[1]:sim.domain.indexes[2]]
        @views dSdt[sim.domain.indexes[1]:sim.domain.indexes[2],:] ./= ns
    end
    return dSdt
end

function normalizeparamtransitorysensitivities!(dSdt,sim::Simulation,t,ind)
    y = sim.sol(t)
    @views ns = y[sim.domain.indexes[1]:sim.domain.indexes[2]]
    dSdt .*= sim.p[ind]
    @views dSdt[sim.domain.indexes[1]:sim.domain.indexes[2]] ./= ns
    return dSdt
end

function normalizeparamtransitorysensitivities!(dSdt,ssys::SystemSimulation,t,ind)
    y = ssys.sol(t)
    for sim in ssys.sims
        @views ns = y[sim.domain.indexes[1]:sim.domain.indexes[2]]
        dSdt .*= ssys.p[ind]
        @views dSdt[sim.domain.indexes[1]:sim.domain.indexes[2]] ./= ns
    end
    return dSdt
end

function normalizeadjointtransitorysensitivities!(dSdt,sim::Simulation,t,ind)
    y = sim.sol(t)
    dSdt .*= sim.p'
    if ind <= sim.domain.indexes[2]
        dSdt ./= y[ind]
    end
    return dSdt
end

function normalizeadjointtransitorysensitivities!(dSdt,ssys::SystemSimulation,t,ind)
    y = ssys.sol(t)
    for sim in ssys.sims
        @views dSdt[sim.domain.parameterindexes[1]:sim.domain.parameterindexes[2]] .*= ssys.domain.p[sim.domain.parameterindexes[1]:sim.domain.parameterindexes[2]]
        if ind >= sim.domain.indexes[1] && ind <= sim.domain.indexes[2]
            dSdt ./= y[ind]
        end
    end
    return dSdt
end

"""
Compute exact transitory sensitivities using Forward Sensitivity Analysis
"""
function transitorysensitivitiesfullexact(sim::Simulation,t;tau=NaN,
        normalized=true,solver=Sundials.CVODE_BDF(linear_solver=:GMRES),
        abstol=1e-16,reltol=1e-6)

    if isnan(tau)
        Jy = jacobiany(sim.sol,t,sim.p);
        tau = getdividingtimescale(Jy);
    end

    if tau == 0.0
        dSdt = jacobianp(sim.sol,t,sim.p);
    else
        tspan = (0.0,tau)
        react = Reactor(sim.domain,sim.sol(t),tspan,sim.interfaces;p=sim.sol.prob.p,forwardsensitivities=true);
        sol = solve(react.ode,solver,abstol=abstol,reltol=reltol)
        dSdt = reduce(hcat,DiffEqSensitivity.extract_local_sensitivities(sol,tau)[2])./tau
    end

    if normalized
        return normalizefulltransitorysensitivities!(dSdt,sim,t)
    else
        return dSdt
    end
end

"""
Compute exact transitory sensitivities using Forward Sensitivity Analysis
"""
function transitorysensitivitiesfullexact(ssys::SystemSimulation,t;tau=NaN,
        normalized=true,solver=Sundials.CVODE_BDF(linear_solver=:GMRES),
        abstol=1e-16,reltol=1e-6)

    if isnan(tau)
        Jy = jacobiany(ssys.sol,t,ssys.p);
        tau = getdividingtimescale(Jy);
    end

    if tau == 0.0
        dSdt = jacobianp(ssys.sol,t,ssys.p);
    else
        tspan = (0.0,tau)
        react = Reactor(ssys.domains,sim.sol(t),tspan,ssys.interfaces;p=ssys.p,forwardsensitivities=true);
        sol = solve(react.ode,solver,abstol=abstol,reltol=reltol)
        dSdt = reduce(hcat,DiffEqSensitivity.extract_local_sensitivities(sol,tau)[2])./tau
    end

    if normalized
        return normalizefulltransitorysensitivities!(dSdt,sim,t)
    else
        return dSdt
    end
end
export transitorysensitivitiesfullexact

"""
Compute approximate transitory sensitivities
Assumes constant jacobian y and p and then uses trapezoidal rule to
approximate the solution to the resulting forward sensitivity equations
requires one evaluation of both jacobiany and jacobianp
"""
function transitorysensitivitiesfulltrapezoidal(sim,t;tau=NaN,normalized=true)
    if tau == 0.0
        dSdt = jacobianp(sim.sol,t,sim.p);
    elseif isnan(tau)
        Jy = jacobiany(sim.sol,t,sim.p);
        tau = getdividingtimescale(Jy);
        Jp = jacobianp(sim.sol,t,sim.p);
        dSdt = (Jp .+ exp(tau*Jy)*Jp)./2.0
    else
        Jy = jacobiany(sim.sol,t,sim.p);
        Jp = jacobianp(sim.sol,t,sim.p);
        dSdt = (Jp .+ exp(tau*Jy)*Jp)./2.0
    end
    if normalized
        return normalizefulltransitorysensitivities!(dSdt,sim,t)
    else
        return dSdt
    end
end
export transitorysensitivitiesfulltrapezoidal

"""
Compute approximate transitory sensitivities with respect to one parameter
Assumes constant jacobian y and p and then uses trapezoidal rule to
approximate the solution to the resulting forward sensitivity equations
requires one evaluation of both jacobiany and jacobianp
"""
function transitorysensitivitiesparamtrapezoidal(sim,t,ind;tau=NaN,normalized=true)
    if tau == 0.0
        dSdt = jacobianp(sol,t,sim.p,ind);
    elseif isnan(tau)
        Jy = jacobiany(sim.sol,t,sim.p);
        tau = getdividingtimescale(Jy);
        Jp = jacobianp(sim.sol,t,sim.p,ind);
        dSdt = (Jp .+ exp(tau*Jy)*Jp)./2.0
    else
        Jy = jacobiany(sim.sol,t,sim.p);
        Jp = jacobianp(sim.sol,t,sim.p,ind);
        dSdt = (Jp .+ exp(tau*Jy)*Jp)./2.0
    end
    if normalized
        return normalizeparamtransitorysensitivities!(dSdt,sim,t)
    else
        return dSdt
    end
end
export transitorysensitivitiesparamtrapezoidal

"""
Compute exact transitory sensitivities with respect to one parameter
using forward sensitivity analysis
"""
function transitorysensitivitiesparamexact(sim::Simulation,t,ind;tau=NaN,
        normalized=true,solver=Sundials.CVODE_BDF(linear_solver=:GMRES),
        abstol=1e-16,reltol=1e-6)

    if isnan(tau)
        Jy = jacobiany(sim.sol,t,sim.p);
        tau = getdividingtimescale(Jy);
    end

    if tau == 0.0
        dSdt = jacobianp(sim.sol,t,sim.p,ind);
    else
        tspan = (0.0,tau)
        react = Reactor(sim.domain,sim.sol(t),tspan,sim.interfaces;p=sim.sol.prob.p,forwardsensitivities=true);
        fparam(u,p,t) = react.ode.f(u,hcat(sim.p[1:ind-1],[sim.p[ind]],sim.p[ind+1:end]),t)
        jacparam(u,p,tn) = jacobianp(sol,t+tn,p,ind)
        odeparam = remake(react.ode,f=fparam,paramjac=jacparam)
        sol = solve(odeparam,solver,abstol=abstol,reltol=reltol)
        dSdt = reduce(hcat,DiffEqSensitivity.extract_local_sensitivities(sol,tau)[2])./tau
    end

    if normalized
        return normalizeparamtransitorysensitivities!(dSdt,sim,t)
    else
        return dSdt
    end
end

"""
Compute exact transitory sensitivities with respect to one parameter
using forward sensitivity analysis
"""
function transitorysensitivitiesparamexact(ssys::SystemSimulation,t,ind;tau=NaN,
        normalized=true,solver=Sundials.CVODE_BDF(linear_solver=:GMRES),
        abstol=1e-16,reltol=1e-6)

    if isnan(tau)
        Jy = jacobiany(ssys.sol,t,ssys.p);
        tau = getdividingtimescale(Jy);
    end

    if tau == 0.0
        dSdt = jacobianp(ssys.sol,ssys.p,t,ind);
    else
        tspan = (0.0,tau)
        react = Reactor(ssys.domains,sim.sol(t),tspan,ssys.interfaces;p=ssys.p,forwardsensitivities=true);
        fparam(u,p,t) = react.ode.f(u,hcat(sim.sol.prob.p[1:ind-1],[sim.sol.prob.p[ind]],sim.sol.prob.p[ind+1:end]),t)
        jacparam(u,p,tn) = jacobianp(sol,t+tn,p,ind)
        odeparam = remake(react.ode,f=fparam,paramjac=jacparam)
        sol = solve(odeparam,solver,abstol=abstol,reltol=reltol)
        dSdt = reduce(hcat,DiffEqSensitivity.extract_local_sensitivities(sol,tau)[2])./tau
    end

    if normalized
        return normalizeparamtransitorysensitivities!(dSdt,sim,t)
    else
        return dSdt
    end
end
export transitorysensitivitiesparamexact

"""
compute exact transitory sensitivities using adjoint sensitivity analysis
abstol and reltol are for the adjoint simulation
appreciable speed up with minimal accuracy reduction can be achieved by
loosening these tolerances
does not run any jacobian computations if tau != 0, tau != NaN  is supplied and
the solver is jacobian free
"""
function transitorysensitivitiesadjointexact(sim::Simulation,t,name;tau=NaN,
        normalized=true,solver=Sundials.CVODE_BDF(),sensalg=InterpolatingAdjoint(),
        abstol=1e-16,reltol=1e-6)

    @assert name in sim.names || name in ["T","V","P"]
    if name in ["T","V","P"]
        if haskey(sim.domain.thermovariabledict, name)
            ind = sim.domain.thermovariabledict[name]
        else
            throw(error("$(sim.domain) doesn't have $name in its thermovariables"))
        end
    else
        ind = findfirst(isequal(name),sim.names)
    end

    if isnan(tau)
        Jy = jacobiany(sim.sol,t,sim.p);
        tau = getdividingtimescale(Jy);
    end

    if tau == 0.0
        dSdt = jacobianp(sim.sol,t,sim.p)[ind,:];
    else
        prob = remake(sim.sol.prob;u0=sim.sol(t),tspan=(0.0,tau));
        soladj = solve(prob,solver,reltol=reltol,abstol=abstol);
        simadj = Simulation(soladj,sim.domain);
        dSdt = getadjointsensitivities(simadj,name,solver;sensalg=sensalg,abstol=abstol,
            reltol=reltol,normalize=false)./tau;
    end

    if normalized
        return normalizeadjointtransitorysensitivities!(dSdt,sim,t,ind);
    else
        return dSdt
    end
end

"""
compute exact transitory sensitivities using adjoint sensitivity analysis
abstol and reltol are for the adjoint simulation
appreciable speed up with minimal accuracy reduction can be achieved by
loosening these tolerances
does not run any jacobian computations if tau != 0, tau != NaN  is supplied and
the solver is jacobian free
"""
function transitorysensitivitiesadjointexact(ssys::SystemSimulation,t,name;tau=NaN,
        normalized=true,solver=Sundials.CVODE_BDF(),sensalg=InterpolatingAdjoint(),
        abstol=1e-6,reltol=1e-3)

    @assert name in sim.names || name in ["T","V","P"]
    if name in ["T","V","P"]
        if haskey(sim.domain.thermovariabledict, name)
            ind = sim.domain.thermovariabledict[name]
        else
            throw(error("$(sim.domain) doesn't have $name in its thermovariables"))
        end
    else
        ind = findfirst(isequal(name),sim.names)
    end

    if isnan(tau)
        Jy = jacobianysparse(ssys.sol,t,ssys.p);
        tau = getdividingtimescale(Jy);
    end

    if tau == 0.0
        dSdt = jacobianp(ssys.sol,t,ssys.p)[ind,:];
    else
        prob = remake(ssys.sol.prob;u0=sim.sol(t),tspan=(0.0,tau));
        soladj = solve(prob,solver,reltol=reltol,abstol=abstol);
        ssysadj = SystemSimulation(soladj,[sim.domain for sim in ssys.sims],ssys.interfaces,ssys.p);
        dSdt = getadjointsensitivities(ssysadj,name,solver;sensalg=sensalg,abstol=abstol,reltol=reltol)./tau;
    end

    if normalized
        return normalizeadjointtransitorysensitivities!(dSdt,ssys,t,ind)
    else
        return dSdt
    end
end
export transitorysensitivitiesadjointexact

function lmprod(lm,x,p,t,f)
    dy = similar(x,length(x))
    return dot(lm,f(dy,x,p,t))
end

function dlmdt(lm,p,t,f,gu,sol)
    return gu.-ReverseDiff.gradient(x->lmprod(lm,x,p,t,f),sol(t))
end

"""
Compute approximate transitory sensitivities using adjoint sensitivity analysis
solver, abstol and reltol are for the main ode simulation while solveradj,
abstoadj and reltoladj are for the adjoint simulation
Npoints is the number of quadrature points to use
does not run any jacobian computations if tau != 0, tau != NaN is supplied and
the solver is jacobian free
"""
function transitorysensitivitiesadjointapprox(sim::Simulation,t,name;tau=NaN,
        normalized=true,solver=Sundials.CVODE_BDF(),
        solveradj=Sundials.CVODE_BDF(linear_solver=:GMRES),
        abstol=1e-16,reltol=1e-6,abstoladj=1e0,reltoladj=1e0,Npoints=8)

    @assert name in sim.names || name in ["T","V","P"]
    if name in ["T","V","P"]
        if haskey(sim.domain.thermovariabledict, name)
            ind = sim.domain.thermovariabledict[name]
        else
            throw(error("$(sim.domain) doesn't have $name in its thermovariables"))
        end
    else
        ind = findfirst(isequal(name),sim.names)
    end

    if isnan(tau)
        Jy = jacobiany(sim.sol,t,sim.p);
        tau = getdividingtimescale(Jy);
    end

    if tau == 0.0
        dSdt = jacobianp(sim.sol,t,sim.ps)[ind,:];
    else
        prob = remake(sim.sol.prob;u0=sim.sol(t),tspan=(0.0,tau));
        soladj = solve(prob,solver,reltol=reltol,abstol=abstol);
        simadj = Simulation(soladj,sim.domain);

        lm0 = zeros(length(sim.sol.prob.u0))
        function f1(x)
            dy = similar(x,length(x))
            return soladj.prob.f(dy,x,soladj.prob.p,0.0)[ind]
        end
        gu = ReverseDiff.gradient(f1,soladj.prob.u0)
        dlm(lm,p,t) = dlmdt(lm,simadj.sol.prob.p,t,simadj.sol.prob.f,gu,soladj)
        ode = ODEProblem(dlm,lm0,(tau,0.0))
        lmsol = solve(ode,solveradj,reltol=reltoladj,abstola=abstoladj);

        function G1(p,t)
            dy = similar(p,length(soladj.prob.u0))
            return soladj.prob.f(dy,soladj(t),p,t)[ind]
        end
        function G2(p,t)
            dy = similar(p,length(soladj.prob.u0))
            return soladj.prob.f(dy,soladj(t),p,t)
        end
        gp = ReverseDiff.gradient(x->G1(x,0.0),soladj.prob.p)
        G(t) = ReverseDiff.gradient(x->dot(lmsol(t),G2(x,t)),soladj.prob.p) .+ gp
        nodes,weights = gausslegendre(Npoints);
        Gvals = G.((nodes.+1.0).*tau./2.0)
        dSdt = reduce(hcat,Gvals)*weights
    end
    if normalized
        return normalizeadjointtransitorysensitivities!(dSdt',sim,t,ind)
    else
        return dSdt
    end
end

"""
Compute approximate transitory sensitivities using adjoint sensitivity analysis
solver, abstol and reltol are for the main ode simulation while solveradj,
abstoadj and reltoladj are for the adjoint simulation
Npoints is the number of quadrature points to use
does not run any jacobian computations if tau != 0, tau != NaN is supplied and
the solver is jacobian free
"""
function transitorysensitivitiesadjointapprox(ssys::SystemSimulation,t,name;tau=NaN,
        normalized=true,solver=Sundials.CVODE_BDF(),
        solveradj=Sundials.CVODE_BDF(linear_solver=:GMRES),
        abstol=1e-16,reltol=1e-6,abstoladj=1e0,reltoladj=1e0,Nquadpoints=8)

    @assert name in sim.names || name in ["T","V","P"]
    if name in ["T","V","P"]
        if haskey(sim.domain.thermovariabledict, name)
            ind = sim.domain.thermovariabledict[name]
        else
            throw(error("$(sim.domain) doesn't have $name in its thermovariables"))
        end
    else
        ind = findfirst(isequal(name),sim.names)
    end

    if isnan(tau)
        Jy = jacobiany(ssys.sol,t,ssys.p);
        tau = getdividingtimescale(Jy);
    end

    if tau == 0.0
        dSdt = jacobianp(ssys.sol,t,ssys.p)[ind,:];
    else
        prob = remake(ssys.sol.prob;u0=sim.sol(t),tspan=(0.0,tau));
        soladj = solve(prob,solver,reltol=reltol,abstol=abstol);
        ssysadj = SystemSimulation(soladj,[sim.domain for sim in ssys.sims],ssys.interfaces,ssys.p);

        lm0 = zeros(length(sim.sol.prob.u0))
        gu = ReverseDiff.gradient(x->soladj.prob.f(x,soladj.prob.p,0.0)[ind],soladj.prob.u0)
        dlm(lm,p,t) = dlmdt(lm,p,t,simadj.sol.prob.f,gu,soladj)
        ode = ODEProblem(dlm,lm0,(tau,0.0))
        lmsol = solve(ode,solveradj,reltol=reltoladj,abstol=abstoladj);

        gp = ReverseDiff.gradient(p->soladj.prob.f(soladj.prob.u0,p,0.0)[ind],soladj.prob.p)
        G(t) = ReverseDiff.gradient(p->soladj.prob.f(soladj.prob.u0,p,t)[ind],lmsol(t)) .+ gp
        nodes,weights = gausslegendre(Nquadpoints);
        Gvals = G.((nodes.+1.0).*tau./2.0)
        dSdt = reduce(hcat,Gvals)*weights
    end

    if normalized
        return normalizeadjointtransitorysensitivities!(dSdt,ssys,t,ind)
    else
        return dSdt
    end
end


text:
function f() x = 0 ; while x < 10; y = x; x = x + 1; end end
print(y)

max_stars_repo_name: gabriel-doriath-dohler/pjulia

max_issues_repo_name: gabriel-doriath-dohler/pjulia

max_forks_repo_name: gabriel-doriath-dohler/pjulia


text:
const NEWICK_DELIM = ('(', ')', ',', ':', ';')
const COMMENT_DELIM = ('[', ']')
const ALL_DELIMS = union(NEWICK_DELIM, COMMENT_DELIM)

"""
Count the number of nodes in a Newick stream
"""
function _countnodes(io)

    #TODO: Use this function for pre-allocating the node vector of the tree
    #? Could use a `getfreenode` function that returns an unlinked node in the node vector
    #? or a new node created with `addnode!`

    #* The logic below could be simplified using `skipchars`

    i = 0
    chr = peek(io, Char)
    while chr ≠ ';'
        i += ((chr == '(') | (chr == ',')) ? 1 : 0
        chr = read(io, Char)
    end

    return i
end

"""
Read IO until any character in a collection of delimiters is found.
"""
function _readuntilany(io, delimset)

    #TODO: Re-implement this method using `skipchars`

    s = ""
    chr = peek(io, Char)
    while chr ∉ delimset
        s *= read(io, Char)
        chr = peek(io, Char)
    end

    return s
end

"""
Check that `query` is the current character in a IOStram, and read it. Raise exception if
it isn't.
"""
function _readconfirm(io, query)
    chr = read(io, Char)
    if chr != query
        throw(ErrorException("expected \'$(query)\' but got \'$(chr)\'"))
    end

    return nothing
end

function _parsetext(io)
    _whiteslurp(io)
    chr = peek(io, Char)
    if chr == '\"'
        read(io, Char)            # Read opening quote mark
        s = readuntil(io, '\"')   # Read string
    elseif chr == '\''
        read(io, Char)
        s = readuntil(io, '\'')
    else
        s = _readuntilany(io, ALL_DELIMS)
    end

    return s
end

function _parsebranch(io)
    #TODO: Check for a comment here
    _readconfirm(io, ':')
    s = _readuntilany(io, ALL_DELIMS)

    return isempty(s) ? nothing : parse(Float64, s)
end

function _parseouter(io::IO, tree)

    #? All this could just be included in `_parsechild`

    labelstr = _parsetext(io)
    _whiteslurp(io)
    chr = peek(io, Char)
    #TODO: Check for a comment here
    if chr == ':'
        brlen =  _parsebranch(io)
    else
        brlen = nothing
    end

    c = addnode!(tree, 1)
    setlabel!(c, labelstr)
    brlength!(c, brlen)

    return c
end

function _parsechild(io, tree)
    chr = peek(io, Char)
    if chr == '('
        c = _parseinner(io, tree)
    else
        c = _parseouter(io, tree)
    end
    #TODO: Check for a comment here

    return c
end

function _parseinner(io::IO, tree)
    _readconfirm(io, '(')
    _whiteslurp(io)

    p = addnode!(tree, 1)

    chr = peek(io, Char)   #? Not necessary
    c = _parsechild(io, tree)
    link!(p => c)
    _whiteslurp(io)
    chr = peek(io, Char)

    while chr == ','
        _readconfirm(io, ',')    # Eat the comma
        _whiteslurp(io)
        c = _parsechild(io, tree)
        link!(p => c)
        _whiteslurp(io)
        chr = peek(io, Char)
    end

    _readconfirm(io, ')')

    labelstr = _parsetext(io)
    chr = peek(io, Char)
    brlen = chr == ':' ? _parsebranch(io) : nothing

    setlabel!(p, labelstr)
    brlength!(p, brlen)

    return p
end

_whiteslurp(io) = skipchars(isspace, io)

function parsenewick(io::IO)
    tree = createtree(0)

    _whiteslurp(io)
    tree.root = _parseinner(io, tree)
    _whiteslurp(io)

    #TODO: Check for a comment here

    _readconfirm(io, ';')

    return tree
end

parsenewick(s::String) = parsenewick(IOBuffer(s))

max_stars_repo_name: eascarrunz/RootedTrees

max_issues_repo_name: eascarrunz/RootedTrees

max_forks_repo_name: eascarrunz/RootedTrees


text:
"""
    GaussianLikelihood(σ²::T=1e-3) # σ² is the **variance**

Gaussian noise :
```math
    p(y|f) = N(y|f,σ²)
```
There is no augmentation needed for this likelihood which is already conjugate to a Gaussian prior
"""
struct GaussianLikelihood{T<:Real,O,A<:AbstractVector{T}} <: RegressionLikelihood{T}
    σ²::A
    opt_noise::O
    θ::A
    function GaussianLikelihood{T}(σ²::T, opt_noise) where {T<:Real}
        return new{T,typeof(opt_noise),Vector{T}}([σ²], opt_noise)
    end
    function GaussianLikelihood{T}(
        σ²::T, opt_noise, θ::A
    ) where {T<:Real,A<:AbstractVector{T}}
        return new{T,typeof(opt_noise),A}(A([σ²]), opt_noise, θ)
    end
end

function GaussianLikelihood(σ²::T=1e-3; opt_noise=false) where {T<:Real}
    if isa(opt_noise, Bool)
        opt_noise = opt_noise ? ADAM(0.05) : nothing
    end
    return GaussianLikelihood{T}(σ², opt_noise)
end

implemented(::GaussianLikelihood, ::Union{<:AnalyticVI,<:Analytic}) = true

function (l::GaussianLikelihood)(y::Real, f::Real)
    return pdf(Normal(y, sqrt(noise(l))), f)
end

function Distributions.loglikelihood(l::GaussianLikelihood, y::Real, f::Real)
    return logpdf(Normal(y, sqrt(noise(l))), f)
end

noise(l::GaussianLikelihood) = first(l.σ²)

function Base.show(io::IO, l::GaussianLikelihood)
    return print(io, "Gaussian likelihood (σ² = $(noise(l)))")
end

function compute_proba(
    l::GaussianLikelihood{T}, μ::AbstractVector{<:Real}, σ²::AbstractVector{<:Real}
) where {T<:Real}
    return μ, σ² .+ noise(l)
end

function init_likelihood(
    likelihood::GaussianLikelihood{T}, ::AbstractInference, ::Int, nSamplesUsed::Int
) where {T}
    return GaussianLikelihood{T}(
        noise(likelihood), likelihood.opt_noise, fill(inv(noise(likelihood)), nSamplesUsed)
    )
end

function local_updates!(
    l::GaussianLikelihood, y::AbstractVector, μ::AbstractVector, var_f::AbstractVector
)
    if !isnothing(l.opt_noise)
        grad = 0.5 * ((sum(abs2, y - μ) + sum(var_f)) / noise(l) - length(y))
        l.σ² .= exp.(log.(l.σ²) + Optimise.apply!(l.opt_noise, l.σ², [grad]))
    end
    return l.θ .= inv(noise(l))
end

@inline function ∇E_μ(l::GaussianLikelihood{T}, ::AOptimizer, y::AbstractVector) where {T}
    return (y ./ noise(l),)
end

@inline function ∇E_Σ(l::GaussianLikelihood{T}, ::AOptimizer, ::AbstractVector) where {T}
    return (0.5 * l.θ,)
end

function expec_loglikelihood(
    l::GaussianLikelihood,
    ::AnalyticVI,
    y::AbstractVector,
    μ::AbstractVector,
    diagΣ::AbstractVector,
)
    return -0.5 * (
        length(y) * (log(twoπ) + log(noise(l))) + (sum(abs2, y - μ) + sum(diagΣ)) / noise(l)
    )
end

AugmentedKL(::GaussianLikelihood{T}, ::AbstractVector) where {T} = zero(T)


text:
module TestReduce

using Test
using SyncBarriers
using SyncBarriers.Internal: BarrierHandle, acctype

# An example of non-commutative reduction
function larger((i, a), (j, b))
    if a < b
        (j, b)
    else
        (i, a)
    end
end

initof(op, T) = Base.reduce_empty(op, T)
initof(::typeof(xor), T::Type{<:Integer}) = zero(T)
initof(::typeof(larger), ::Type{Tuple{I,X}}) where {I,X} = (zero(I), typemin(X))

samples(::Type{T}) where {T<:Real} = T(100):T(999)
samples(::Type{T}) where {I,X,T<:Tuple{I,X}} = collect(zip(samples(I), samples(X)))

function test_reduce()
    @testset "+" begin
        test_reduce(+, Int)
    end
    @testset "xor" begin
        test_reduce(xor, UInt)
    end
    @testset "larger" begin
        test_reduce(larger, Tuple{Int,Float64})
    end
end

struct TreeBarrier′{NBranches,T,B<:SyncBarriers.TreeBarrier{NBranches,T}} <: SyncBarriers.Barrier
    barrier::B
end

Base.length(barrier::TreeBarrier′) = length(barrier.barrier)

TreeBarrier′{NBranches,T}(op, n) where {NBranches,T} =
    TreeBarrier′(SyncBarriers.TreeBarrier{NBranches,T}(op, n))
SyncBarriers.Internal.acctype(::Type{B}) where {T,B<:TreeBarrier′{<:Any,T}} = T

function SyncBarriers.reduce!(
    handle::BarrierHandle{<:TreeBarrier′},
    value,
    spin::Union{Nothing,Integer} = nothing,
)
    handle = BarrierHandle(handle.barrier.barrier, handle.i)
    SyncBarriers.reduce_arrive!(handle, value)
    SyncBarriers.depart!(handle, spin)
end

function test_reduce(op, T)
    barrier_factories = [
        SyncBarriers.StaticTreeBarrier{2,2,T},
        SyncBarriers.StaticTreeBarrier{3,4,T},
        SyncBarriers.TreeBarrier{3,T},
        SyncBarriers.FlatTreeBarrier{2,T},
        SyncBarriers.FlatTreeBarrier{3,T},
        TreeBarrier′{2,T},
        TreeBarrier′{3,T},
        # ...
    ]
    ntasks_list = unique([1, 2, Threads.nthreads(), 2 * Threads.nthreads()])
    @testset for f in barrier_factories, ntasks in ntasks_list
        @debug "Testing `test_reduce($op, $f, $ntasks)`"
        test_reduce(op, f, ntasks)
    end
end

function test_reduce(op, factory, ntasks)
    output, input = use_barrier(op, factory, ntasks)
    desired =
        repeat(reduce(op, input; init = initof(op, eltype(output)), dims = 2), 1, ntasks)
    @test output[:, 1] == desired[:, 1]
    @test all(output[:, 1] .== output)
    @test output == desired
end

function use_barrier(op, factory, ntasks; kwargs...)
    barrier = factory(op, ntasks)
    if factory isa Type
        @test barrier isa factory
    end
    return use_barrier(barrier; kwargs...)
end

function use_barrier(barrier::SyncBarriers.Barrier; ncycles = 1000)
    ntasks = length(barrier)
    input = rand(samples(acctype(barrier)), ncycles, ntasks)
    output = similar(input)
    return use_barrier!(output, input, barrier)
end

function use_barrier!(
    output::Matrix,
    input::Matrix,
    barrier::SyncBarriers.Barrier,
    spin = nothing,
)
    ncycles, ntasks = size(output)
    @assert ntasks == length(barrier)
    @sync for i in 1:ntasks
        Threads.@spawn try
            for k in 1:ncycles
                output[k, i] = SyncBarriers.reduce!(barrier[i], input[k, i], spin)
            end
        catch err
            @error(
                "`use_barrier` failed",
                exception = (err, catch_backtrace()),
                i,
                current_task()
            )
            # TODO: close(barrier)
            rethrow()
        end
    end
    return output, input
end

end  # module


text:
# RUN_LIST as input
include("run_cutest.jl")

problem_list = default_list()
#problem_list = get_problem_list(20,100)

@show length(problem_list)
# :dynamic,
#RUN_LIST = ARGS
#RUN_LIST = [:corrections, :step_style, :regularizer, :tol]
RUN_LIST = ["Jan4"] #["LowTol_Dec26","test_Dec26_ftb"]#"LowTol_Dec19",

if OnePhase.USE_HSL
    OnePhase.loadHSL("../../src/linear_system_solvers/")
end

if "Jan4" in RUN_LIST

my_par = OnePhase.Class_parameters()
my_par.term.max_time = 60.0 * 60
my_par.term.max_it = 3000

folder_name = "one_phase/Jan4/CUTEst"
if_mkdir("../results/$folder_name")
run_cutest_problems_using_our_solver(problem_list, folder_name, my_par)

my_par = OnePhase.Class_parameters()
my_par.term.max_time = 60.0 * 60
my_par.term.max_it = 3000
my_par.term.tol_opt = 1e-2
my_par.term.tol_inf_1 = 1e-2
my_par.term.tol_inf_2 = 1e-2
#folder_name = "one_phase/test_meh_init2_Dec30"
folder_name = "one_phase/Jan4/LowTol"
if_mkdir("../results/$folder_name")
run_cutest_problems_using_our_solver(problem_list, folder_name, my_par)
end

#=if "Dec26" in RUN_LIST
my_par = Class_parameters()
my_par.term.max_time = 60.0 * 60
my_par.term.max_it = 3000
folder_name = "one_phase/Dec26"
if_mkdir("../results/$folder_name")
run_cutest_problems_using_our_solver(problem_list, folder_name, my_par)
end=#

if "test_Dec26_ftb" in RUN_LIST
my_par = Class_parameters()
my_par.term.max_time = 60.0 * 60
my_par.term.max_it = 3000
my_par.ls.fraction_to_boundary_predict = 0.3
folder_name = "one_phase/test_Dec26_ftb"
if_mkdir("../results/$folder_name")
run_cutest_problems_using_our_solver(problem_list, folder_name, my_par)
end

if "LowTol_Dec26" in RUN_LIST
my_par = Class_parameters()
my_par.term.max_time = 60.0 * 60
my_par.term.max_it = 3000
my_par.term.tol_opt = 1e-2
my_par.term.tol_inf_1 = 1e-2
my_par.term.tol_inf_2 = 1e-2
folder_name = "one_phase/LowTol_Dec26"
if_mkdir("../results/$folder_name")
run_cutest_problems_using_our_solver(problem_list, folder_name, my_par)
end

if "init" in RUN_LIST
my_par.term.max_time = 10.0 * 60.0
my_par.term.max_it = 1000
if_mkdir("../results/one_phase/init")

my_par.init.mehotra_scaling = false

my_par.init.mu_scale = 0.1
folder_name = "one_phase/init/0.1"
if_mkdir("../results/$folder_name")
run_cutest_problems_using_our_solver(problem_list, folder_name, my_par)

my_par.init.mu_scale = 1.0
folder_name = "one_phase/init/1"
if_mkdir("../results/$folder_name")
run_cutest_problems_using_our_solver(problem_list, folder_name, my_par)

my_par.init.mu_scale = 10.0
folder_name = "one_phase/init/10"
if_mkdir("../results/$folder_name")
run_cutest_problems_using_our_solver(problem_list, folder_name, my_par)

my_par.init.mu_scale = 100.0
folder_name = "one_phase/init/100"
if_mkdir("../results/$folder_name")
run_cutest_problems_using_our_solver(problem_list, folder_name, my_par)

my_par.init.mu_scale = 1000.0
folder_name = "one_phase/init/1000"
if_mkdir("../results/$folder_name")
run_cutest_problems_using_our_solver(problem_list, folder_name, my_par)
end

#problem_list = get_problem_list(100,200)

#test_problems(problem_list,1)
if false
problem_list = get_problem_list(6000,7000)
for problem in problem_list
    @show problem
    nlp_raw = CUTEstModel(problem)
    @show nlp_raw.meta.nnzj, nlp_raw.meta.nnzh
    @show nlp_raw.meta.ncon, nlp_raw.meta.nvar
    A = jac(nlp_raw,nlp_raw.meta.x0);
    @show densest_row(A), densest_col(A)

    finalize(nlp_raw)
end
end

if "penalty" in RUN_LIST
if_mkdir("../results/one_phase/Oct22")
folder_name = "one_phase/Oct22/penalty"
my_par.x_norm_penalty = 1e-8
my_par.a_norm_penalty = 1e-4
my_par.use_prox = true
my_par.use_reg = true
if_mkdir("../results/$folder_name")
run_cutest_problems_using_our_solver(problem_list, folder_name, my_par)

folder_name = "one_phase/Oct22/penalty_no_prox"
my_par.x_norm_penalty = 1e-8
my_par.a_norm_penalty = 1e-4
my_par.use_prox = false
my_par.use_reg = true
if_mkdir("../results/$folder_name")
run_cutest_problems_using_our_solver(problem_list, folder_name, my_par)

folder_name = "one_phase/Oct22/a_penalty_only"
my_par.x_norm_penalty = 1e-30
my_par.a_norm_penalty = 1e-4
my_par.use_prox = true
my_par.use_reg = true
if_mkdir("../results/$folder_name")
run_cutest_problems_using_our_solver(problem_list, folder_name, my_par)

folder_name = "one_phase/Oct22/a_penalty_only_no_prox"
my_par.x_norm_penalty = 1e-30
my_par.a_norm_penalty = 1e-4
my_par.use_prox = false
my_par.use_reg = true
if_mkdir("../results/$folder_name")
run_cutest_problems_using_our_solver(problem_list, folder_name, my_par)
end

if "plain" in RUN_LIST
    folder_name = "one_phase/plain"
    if_mkdir("../results/$folder_name")
    run_cutest_problems_using_our_solver(problem_list, folder_name, my_par)
end

if "dynamic" in RUN_LIST
    folder_name = "one_phase/sept_dynamic"
    if_mkdir("../results/$folder_name")
    my_par.adaptive_mu = :paper2
    my_par.primal_bounds_dual_feas = true
    run_cutest_problems_using_our_solver(problem_list, folder_name, my_par)
    my_par.adaptive_mu = :none
    my_par.primal_bounds_dual_feas = false
end

if "corrections" in RUN_LIST
    #folder_name = "one_phase/sept_1_corrections"
    #if_mkdir("../results/$folder_name")
    #my_par.max_it_corrections = 1
    #run_cutest_problems_using_our_solver(problem_list, folder_name, my_par)

    folder_name = "one_phase/sept_2_corrections"
    if_mkdir("../results/$folder_name")
    my_par.max_it_corrections = 2
    run_cutest_problems_using_our_solver(problem_list, folder_name, my_par)

    folder_name = "one_phase/sept_3_corrections"
    if_mkdir("../results/$folder_name")
    my_par.max_it_corrections = 3
    run_cutest_problems_using_our_solver(problem_list, folder_name, my_par)

    folder_name = "one_phase/sept_4_corrections"
    if_mkdir("../results/$folder_name")
    my_par.max_it_corrections = 4
    run_cutest_problems_using_our_solver(problem_list, folder_name, my_par)
end

if "step_style" in RUN_LIST
    folder_name = "one_phase/sept_max_step_stable"
    if_mkdir("../results/$folder_name")
    my_par.adaptive_mu = :paper
    my_par.max_it_corrections = 3
    ls_mode = :accept_aggressive
    my_par.ls_mode_stable_trust = ls_mode
    my_par.ls_mode_stable_delta_zero = ls_mode
    my_par.ls_mode_stable_correction = ls_mode
    run_cutest_problems_using_our_solver(problem_list, folder_name, my_par)

    folder_name = "one_phase/sept_log_barrier_stable"
    if_mkdir("../results/$folder_name")
    my_par.adaptive_mu = :paper
    my_par.max_it_corrections = 3
    ls_mode = :accept_stable
    my_par.ls_mode_stable_trust = ls_mode
    my_par.ls_mode_stable_delta_zero = ls_mode
    my_par.ls_mode_stable_correction = ls_mode
    run_cutest_problems_using_our_solver(problem_list, folder_name, my_par)

    ls_mode = :accept_filter
    my_par.ls_mode_stable_trust = ls_mode
    my_par.ls_mode_stable_delta_zero = ls_mode
    my_par.ls_mode_stable_correction = ls_mode
end

if "regularizer" in RUN_LIST
    folder_name = "one_phase/sept_no_regularizer"
    if_mkdir("../results/$folder_name")
    my_par.use_prox = false
    run_cutest_problems_using_our_solver(problem_list, folder_name, my_par)
    my_par.use_prox = true
end

if "tol" in RUN_LIST
    folder_name = "one_phase/sept_high_tol"
    if_mkdir("../results/$folder_name")
    my_par.tol = 1e-8
    run_cutest_problems_using_our_solver(problem_list, folder_name, my_par)
end

if "dual_ls_style" in RUN_LIST
  folder_name = "one_phase/dual_ls_0"
  if_mkdir("../results/$folder_name")
  my_par.dual_ls = 0
  run_cutest_problems_using_our_solver(problem_list, folder_name, my_par)

  folder_name = "one_phase/dual_ls_1"
  if_mkdir("../results/$folder_name")
  my_par.dual_ls = 1
  run_cutest_problems_using_our_solver(problem_list, folder_name, my_par)

  folder_name = "one_phase/dual_ls_2"
  if_mkdir("../results/$folder_name")
  my_par.dual_ls = 2
  run_cutest_problems_using_our_solver(problem_list, folder_name, my_par)
end

if "eta_strategy" in RUN_LIST
  #folder_name = "one_phase/eta/mehotra"
  #if_mkdir("../results/$folder_name")
  #my_par.agg_eta = :mehrotra
  #run_cutest_problems_using_our_solver(problem_list, folder_name, my_par)

  folder_name = "one_phase/eta/mehrotra_stb"
  if_mkdir("../results/$folder_name")
  my_par.agg_eta = :mehrotra_stb
  run_cutest_problems_using_our_solver(problem_list, folder_name, my_par)

  folder_name = "one_phase/eta/affine"
  if_mkdir("../results/$folder_name")
  my_par.agg_eta = :affine
  run_cutest_problems_using_our_solver(problem_list, folder_name, my_par)
end

max_stars_repo_name: ohinder/OnePhase

max_issues_repo_name: ohinder/OnePhase

max_forks_repo_name: ohinder/OnePhase


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#   Function to run the MGEO Canonical.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# References
#
#   [1] Galski, R.L (2006). Desenvolvimento de versões aprimoradas híbridas,
#       paralela e multiobjetivo do método da otimização extrema generalizada e
#       sua aplicação no projeto de sistemas espaciais (INPE-14795-TDI/1238).
#       PhD Thesis, Instituto Nacional de Pesquisas Espaciais, São José dos
#       Campos, SP, Brazil (in Portuguese).
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

export mgeo_run

"""
    mgeo_run(mgeod::MGEO_Structure{Nv, Nf, Val{:MGEO_Canonical}}, f_obj::Function, show_debug::Bool = false) where {Nv, Nf}

Run the MGEO Canonical configured in `mgeod` using the objective functions
`f_obj`. The algorithm was based in [1].

# Keywords

* `show_debug`: Print debug information (**Default** = `true`).

# Returns

The Pareto frontier, which is an array of `Pareto_Point`.

# Remarks

The objective function must have the following signature

        function f_obj(vars)
            1) Compute f[1], f[2], ..., f[nf] using vars[1], vars[2], ... var[n]
            2) return (valid_point, f)
        end

where `nf` is the number of objective functions, `n` is the number of design
variables, and `valid_point` is a boolean value that indicates if vars yield to
a valid point.

"""
function mgeo_run(mgeod::MGEO_Structure{Nv, Nf, Design_Variable_MGEO_Canonical},
                  f_obj::Function,
                  show_debug::Bool = false) where {Nv, Nf}
    # Pareto frontier.
    pareto = Pareto_Point{Nv,Nf}[]

    # Maximum number of generations per run.
    ngen_max_per_run = floor(mgeod.ngen_max/mgeod.run_max)

    # Get the number of bits the string must have by summing the bits for each
    # design variable. Also, get the maximum number of bits in one variable so
    # that we can compute the `factors` vector to convert the string of bits to
    # a real number.
    num_bits    = 0
    max_bit_num = 0

    @inbounds for k = 1:Nv
        num_bits += mgeod.design_vars[k].bits
        max_bit_num = max(max_bit_num, mgeod.design_vars[k].bits)
    end

    # Compute the vector with the interger factor to convert the string of bits
    # into a real number.
    factors = SVector{max_bit_num, Int64}(2 .^collect(0:1:(max_bit_num-1)))

    # Create the string.
    string = BitArray(undef,num_bits)

    # Initialization of Pareto Frontier
    # =================================

    # Initialize the string with one valid point in the Pareto Frontier.
    @inbounds for run = 1:mgeod.run_max
        # Sample a new string.
        rand!(string)

        # Convert string to real numbers.
        vars = bitstrtonum(mgeod.design_vars, string, factors)

        # Call the objective functions.
        (valid, f) = f_obj(vars)

        # Check if the number of objective functions are valid.
        (length(f) != Nf) && @error("The number of objective functions returned by f_obj is not equal to the one specified in mgeod.")

        if valid
            # Add the point to the Pareto Frontier.
            push!(pareto, Pareto_Point{Nv, Nf}(vars, SVector{Nf, Float64}(f)))
            break
        end
    end

    # Check if a valid point was found.
    if (length(pareto) == 0)
        return pareto
    end

    # Loop - Independent Runs
    # =======================

    # List of all points created after flipping the bits.
    candidate_points = Vector{Pareto_Point{Nv,Nf}}(undef,num_bits)

    # Array to sort the candidate points.
    f_rank = Vector{sRank}(undef, num_bits)

    @inbounds for run = 1:mgeod.run_max
        if show_debug
            println("--------------------------------------------------------------------------------")
            println("RUN = $run")
            println("")
            println("Generations ($ngen_max_per_run):")
            print("    0%")
        end

        # Sample a new string if it is not the first run.
        (run > 1) && rand!(string)

        # Loop - MGEO Generations
        # =======================

        genDebug = 0

        @inbounds for ngen_per_run = 1:ngen_max_per_run
            if (show_debug)
                genDebug += 1

                if (genDebug >= 0.1*ngen_max_per_run)
                    genDebug = 0
                    print("...$(ngen_per_run/ngen_max_per_run*100)%")
                end
            end

            # Choose which objective function will be used to compute the
            # adaptability and to assemble the rank.
            chosen_func = rand(1:Nf)

            # Evaluate the objective functions using parallel computing.
            Threads.@threads for i=1:num_bits
                # Flip the i-th bit.
                string_i              = copy(string)
                @inbounds string_i[i] = !string_i[i]

                # Convert the string into real values.
                vars = bitstrtonum(mgeod.design_vars, string_i, factors)

                # Evaluate the objective functions.
                (valid, f) = f_obj(vars)

                # Check if the solution is valid.
                if valid
                    # Create the candidate point.
                    candidate_point =
                        Pareto_Point{Nv, Nf}(vars, SVector{Nf, Float64}(f))

                    # Add the result to the rank.
                    @inbounds f_rank[i] = sRank(true, i, f[chosen_func])
                    @inbounds candidate_points[i] = candidate_point
                else
                    @inbounds f_rank[i] = sRank(false, i, 0.0)
                end
            end

            # Add the points to the Pareto frontier.
            num_valid_points = 0

            @inbounds for i=1:length(candidate_points)
                if f_rank[i].valid
                    num_valid_points += 1
                    check_dominance!(pareto, candidate_points[i], mgeod)
                end
            end

            # If there are no valid points, an independent run must be called.
            (num_valid_points == 0) && break

            # Ranking
            # ========
            #
            # The points will be ranked according to the selected objective
            # function (`chosen_func`).
            #
            # Notice that the invalid points will be placed after the valid
            # ones.

            sort!(f_rank, lt=(a,b)->begin
                      # Check if data is valid.
                      if (a.valid && b.valid)
                          return a.f < b.f
                      elseif (a.valid)
                          return true
                      else
                          return false
                      end
                  end
                 )

            # Choose a bit to be changed for the new generation.
            bit_accepted = false

            while !bit_accepted
                # Sample a valid bit in the string.
                b_sample = rand(1:num_valid_points)

                # Accept the change with probability r^(-τ), where r is the
                # rank of the bit.
                Pk = b_sample^(-mgeod.τ)

                if rand() <= Pk
                    # The bit is accepted, then exit the loop.
                    @inbounds string[f_rank[b_sample].index] =
                              !string[f_rank[b_sample].index]
                    bit_accepted = true
                end
            end
        end

        if show_debug
            println("")
            println("")
            println("Pareto frontier information:")
            println("    Number of points: $(length(pareto))")
            println("")
        end
    end

    pareto
end


text:
__precompile__()
module EnvelopedArrays

import Base: tail, parent
export EnvelopedArray, envelope

struct EnvelopedArray{T<:Number,N,P} <: AbstractArray{T,N}
   envelope::T
   parent::P
   size::NTuple{N,Int}
end
function EnvelopedArray(envelope::T,parent::P) where {T<:Number,P}
   psize = size(parent)
   esize = (psize[1]+1,tail(psize)...)
   bigtype = promote_type(T,eltype(P))
   pp = convert.(bigtype,parent)
   EnvelopedArray{promote_type(T,eltype(P)),ndims(parent),typeof(pp)}(envelope,pp,esize)
end

# accessors and utilities
parent(A::EnvelopedArray) = A.parent
envelope(A::EnvelopedArray) = A.envelope
Base.size(A::EnvelopedArray) = A.size
Base.similar(A::EnvelopedArray{T,N,P}) where {T,N,P} = EnvelopedArray(A.envelope,zeros(T,size(parent(A))))

# Indexing functions
function Base.getindex(A::EnvelopedArray{T,N,P}, I::Vararg{Int,N}) where {T,N,P}
   if I[1] == 1
      return envelope(A)
   else
      return getindex(parent(A),I[1]-1,tail(I)...)
   end
end

function Base.getindex(A::EnvelopedArray{T,N,P}, i::Int) where {T,N,P}
   @assert i >= 1
   ci = ind2sub(size(A),i)
   if i == 1 || ci[1] == 1
      return envelope(A)
   else
      return getindex(parent(A),ci[1]-1,tail(ci)...)
   end
end

Base.IndexStyle(::Type{<:EnvelopedArray}) = IndexCartesian()

# for now, keep things immutable. to change envelope or parent, just create a new EnvelopedArray
# function Base.setindex(A::EnvelopedArray{T,N,P},x,i::Int) where {T,N,P}
   # @assert i >= 1
   # ci = ind2sub(size(A),i)
   # if i == 1 || ci[1] == 1
      # error("Cannot mutate envelope.")
   # else
      # setindex!(parent(A),x,ci[1]-1,tail(ci)...)
   # end
# end

end


text:
println("===============================================")
println("Testing AFM_importJPK.jl")
println("===============================================")

function _importJPK_data()

    fildir = joinpath(@__DIR__, "testdata", "AFM_sampledata.txt")

    data = importJPK(fildir, AFM(150e-9); sections = ["extend"])

    test1 = (data.t == [5.0e-5, 1.5E-4, 2.5E-4])
    test2 = (data.ϵ == [0.0, 0.016422407656940773, 0.04216370213557422])
    test3 = (data.σ == [0.0, -32711.11111111111, 8000.0])

    test1 && test2 && test3
end

@test _importJPK_data()


text:
using LinearAlgebra
using CompScienceMeshes
using SauterSchwab3D
using StaticArrays

const pI   = point(0,0,0)
const pII  = point(1,0,0)
const pIII = point(0,1,0)
const pIV  = point(0,0,1)

const qIII = point(0,-1,0)
const qIV  = point(0,0,-1)

const P = simplex(pI,pII,pIII,pIV)
const Q = simplex(pI,pII,qIII,qIV)

const sing = SauterSchwab3D.singularity_detection(P,Q)

Accuracy2 = 15
ce_ref = CommonEdge6D(sing,SauterSchwab3D._legendre(Accuracy2,0.0,1.0))

function integrand(x,y)
      return ((x-pI)'*(y-pIV))*exp(-im*1*norm(x-y))/(4pi*norm(x-y))
end

function INTEGRAND(u,v)
   j1 = volume(P) * factorial(dimension(P))
   j2 = volume(Q) * factorial(dimension(Q))
   x = barytocart(P,u)
   y = barytocart(Q,v)
   output = integrand(x,y)*j1*j2

   return output
end

print("Ref: ")
ref = sauterschwab_parameterized(INTEGRAND, ce_ref)
println(ref)
println()

res_tp =[]
res_sp = []
res_gm = []
n1 = []
n2 = []
n3 = []
for i in 2:1:14
   Accuracy = i
   ce = CommonEdge6D(sing,SauterSchwab3D._legendre(Accuracy,0.0,1.0))

   int_tp = sauterschwab_parameterized(INTEGRAND, ce)

   num_pts = 5*length(ce.qps)^6
   push!(n1,num_pts)

   push!(res_tp,int_tp)
end

for i in 2:1:6
   Accuracy = i
   ce_s = CommonEdge6D_S(sing,(SauterSchwab3D._legendre(Accuracy,0.0,1.0),
                          SauterSchwab3D._shunnham2D(Accuracy),
                          SauterSchwab3D._shunnham3D(Accuracy),
                          SauterSchwab3D._shunnham4D(Accuracy)))

   int_sp = sauterschwab_parameterized(INTEGRAND, ce_s)

   num_pts = length(ce_s.qps[2])^3+
             2*length(ce_s.qps[1])*length(ce_s.qps[2])*length(ce_s.qps[3])+
             length(ce_s.qps[2])*length(ce_s.qps[4])+
             length(ce_s.qps[1])^2*length(ce_s.qps[2])^2

   push!(n2,num_pts)
   push!(res_sp,int_sp)
end

for i in 2:1:12
   Accuracy = i
   ce_gm = CommonEdge6D_S(sing,(SauterSchwab3D._legendre(Accuracy,0.0,1.0),
                          SauterSchwab3D._grundmannMoeller2D(2*Accuracy-1),
                          SauterSchwab3D._grundmannMoeller3D(2*Accuracy-1),
                          SauterSchwab3D._grundmannMoeller4D(2*Accuracy-1)))

   int_gm = sauterschwab_parameterized(INTEGRAND, ce_gm)

   num_pts = length(ce_gm.qps[2])^3+
             2*length(ce_gm.qps[1])*length(ce_gm.qps[2])*length(ce_gm.qps[3])+
             length(ce_gm.qps[2])*length(ce_gm.qps[4])+
             length(ce_gm.qps[1])^2*length(ce_gm.qps[2])^2

   push!(n3,num_pts)
   push!(res_gm,int_gm)
end

err_tp = norm.(res_tp.-ref)/norm(ref)
err_sp = norm.(res_sp.-ref)/norm(ref)
err_gm = norm.(res_gm.-ref)/norm(ref)

println(err_tp)
println(err_sp)
println(err_gm)

using Plots
plot( yaxis=:log, xaxis=:log, fontfamily="Times")
plot!(n1,err_tp, label="Gauss Tensor-Product",markershape=:circle)
plot!(n2,err_sp, label="Simplex Tensor-Product",markershape=:rect)
#plot!(n3,err_gm, label="Simplex-Product GM",markershape=:x)
plot!(xlims=(1e2,5e5),ylims=(1e-7,1))
plot!(xlabel="#Quad. pts/Func. evals", ylabel="Rel. Error.", title="Common Edge 6D",legend=:bottomleft)

using BenchmarkTools

ref = -0.0004334592744412087 + 0.0005754325764942029im

ce = CommonEdge6D(sing,SauterSchwab3D._legendre(5,0.0,1.0))

int_tp = sauterschwab_parameterized(INTEGRAND, ce)

num_pts = 5*length(ce.qps)^6
println("#Pts: ",num_pts)
err_tp = norm.(int_tp.-ref)/norm(ref)
println("Rel Err: ",err_tp)
@time for i in 1:100 sauterschwab_parameterized(INTEGRAND, ce) end

ce_s = CommonEdge6D_S(sing,(SauterSchwab3D._legendre(5,0.0,1.0),
                     SauterSchwab3D._shunnham2D(5),
                     SauterSchwab3D._shunnham3D(5),
                     SauterSchwab3D._shunnham4D(5)))

int_sp = sauterschwab_parameterized(INTEGRAND, ce_s)

num_pts = length(ce_s.qps[2])^3+
2*length(ce_s.qps[1])*length(ce_s.qps[2])*length(ce_s.qps[3])+
length(ce_s.qps[2])*length(ce_s.qps[4])+
length(ce_s.qps[1])^2*length(ce_s.qps[2])^2
println("#Pts: ",num_pts)
err_sp = norm.(int_sp.-ref)/norm(ref)
println("Rel Err: ",err_sp)
@time for i in 1:100 sauterschwab_parameterized(INTEGRAND, ce_s) end

#=
ce_gm = CommonEdge6D_S((SauterSchwab3D._legendre(6,0.0,1.0),
                        SauterSchwab3D._grundmannMoeller2D(2*6-1),
                        SauterSchwab3D._grundmannMoeller3D(2*6-1),
                        SauterSchwab3D._grundmannMoeller4D(2*6-1)))

int_gm = sauterschwab_parameterized(INTEGRAND, ce_gm)

num_pts = length(ce_gm.qps[2])^3+
          2*length(ce_gm.qps[1])*length(ce_gm.qps[2])*length(ce_gm.qps[3])+
          length(ce_gm.qps[2])*length(ce_gm.qps[4])+
          length(ce_gm.qps[1])^2*length(ce_gm.qps[2])^2
println("#Pts: ",num_pts)
err_gm = norm.(int_gm.-ref)/norm(ref)
println("Rel Err: ",err_gm)
=#


text:
using Base.Meta: isexpr

function stagem(sig, abst)
    abst = esc(abst)
    if sig isa Symbol
        name = sig
        structdef = :(
            struct $name{T, M} <: $abst
              x::T
              meta::M
            end
        )
    elseif sig isa Expr
        if isexpr(sig, :curly)
            length(sig.args) != 3 && error("invalid TokenStages definition")
            name, Tsig, Msig = sig.args
            T = isexpr(Tsig, :<:) ? Tsig.args[1] : Tsig
            M = isexpr(Msig, :<:) ? Msig.args[1] : Msig
            structdef = :(
                struct $sig <: $abst
                  x::$T
                  meta::$M
                end
            )
        else
            error("invalid TokenStages definition")
        end
    else
        error("invalid TokenStages definition")
    end
    name = esc(name)
    setmeta = esc(:(TextEncodeBase.setmeta))
    setvalue = esc(:(TextEncodeBase.setvalue))
    return quote
        $structdef
        $name(x) = $name(x, nothing)
        $setmeta(x::$name, meta) = $name(x.x, meta)
        $setvalue(x::$name, y) = $name(y, x.meta)
    end
end

"""
    @stage StageName
    @stage StageName{A<:SomeType, B}
    @stage StageName AbstractStage
    @stage StageName{A<:SomeType, B} <: AbstractStage

Define `TokenStages` with two field (`x` and `meta`), it's single arguement constructor,
 and add methods to `setmeta` and `setvalue`.

Equivalent to:

```julia
struct StageName{A<:SomeType, B} <: AbstractStage
    x::A
    meta::B
end

StageName(x) = StageName(x, nothing)
TextEncodeBase.setmeta(x::StageName, meta) = StageName(x.x, meta)
TextEncodeBase.setvalue(x::StageName, y) = StageName(y, x.meta)

```
"""
macro stage(sig, abst=:TokenStages)
    if isexpr(sig, :<:)
        abst != :TokenStages && error("invalid TokenStages definition")
        sig, abst = sig.args
    end
    return stagem(sig, abst)
end


text:
typealias SmallerFloatInt Union{Float32,Float16,Int32,Int16,Int8}
typealias FloatInt Union{Float64,Float32,Float16,Int64,Int32,Int16,Int8}

Float120{T<:SmallerFloatInt}(hi::T) = Float120(hi,zero(T))

function Float120{T<:SmallerFloatInt}(hi::T, lo::T)
    high = Float64(hi)
    low  = Float64(lo)
    high, low = eftAdd(high, low)
    Float120(high,low)
end

Float120{T1<:FloatInt,T2<:FloatInt}(hi::T1, lo::T2) = Float120(Float64(hi), Float64(lo))

function convert(::Type{Float120}, x::BigFloat)
    hi = Float64(x)
    lo = Float64(x-hi)
    Float120(hi,lo)
end

convert(::Type{BigFloat}, x::Float120) = convert(BigFloat,x.hi) + convert(BigFloat,x.lo)

convert(::Type{Float64}, x::Float120) = x.hi
convert(::Type{Float32}, x::Float120) = convert(Float32, x.hi)
convert(::Type{Float16}, x::Float120) = convert(Float16, x.hi)

convert(::Type{Float120}, x::Float64) = Float120(x)
convert(::Type{Float120}, x::Float32) = Float120(convert(Float64,x))
convert(::Type{Float120}, x::Float16) = Float120(convert(Float64,x))

convert(::Type{Int128}, x::Float120) = trunc(Int128,x.hi)+trunc(Int128,x.lo)
convert(::Type{Int64}, x::Float120) = convert(Int64, convert(Int128,x))

function convert(::Type{Int32}, x::Float120)
    if (x.lo == trunc(x.lo)) && (x.hi == trunc(x.hi))
        convert(Int32, trunc(Int64,x.hi)+trunc(Int64,x.lo))
    else
        throw(InexactError())
    end
end

convert(::Type{Int16}, x::Float120) = convert(Int16, convert(Int32,x))
convert(::Type{Int8}, x::Float120) = convert(Int8, convert(Int32,x))

convert(::Type{Float120}, x::Int8) = Float120(convert(Float64,x))
convert(::Type{Float120}, x::Int16) = Float120(convert(Float64,x))
convert(::Type{Float120}, x::Int32) = Float120(convert(Float64,x))

function convert(::Type{Float120}, x::Int64)
    if signbit(x)
        if x < safemin(Float64)
            Float120(convert(BigFloat,convert(BigInt,x)))
        else
            Float120(convert(Float64,x))
        end
    else
        if x > safemax(Float64)
            Float120(convert(BigFloat,convert(BigInt,x)))
        else
            Float120(convert(Float64,x))
        end
    end
end

#const Float120safemax128 = Int128(  207_691_874_341_393_105_096_183_856_895_509_888 )

function convert(::Type{Float120}, x::Int128)
    if signbit(x)
        if x > typemin(Int128)
            Float120(convert(BigFloat,convert(BigInt,x)))
        else
            throw(DomainError())
        end
    else
        if x <= typemax(Int128)
            Float120(convert(BigFloat,convert(BigInt,x)))
        else
            throw(DomainError())
        end
    end
end

# resolve ambiguity
convert(::Type{Bool}, x::Float120) = !(x.hi == 0.0)
convert(::Type{Base.GMP.BigInt}, x::Float120) =
   convert(Base.GMP.BigInt, x.hi) + convert(Base.GMP.BigInt, x.lo)
convert(::Type{Rational{Int128}}, x::Float120) =
   convert(Rational{Int128}, x.hi) + convert(Rational{Int128}, x.lo)
convert(::Type{Rational{Base.GMP.BigInt}}, x::Float120) =
   convert(Rational{Base.GMP.BigInt}, x.hi) + convert(Rational{Base.GMP.BigInt}, x.lo)

function convert{T<:Union{Int128,Base.GMP.BigInt}}(::Type{Float120}, x::Rational{T})
    hi = Float64(x - convert(Rational{T},Float64(x)))
    lo = Float64(x - convert(Rational{T},hi))
    Float120(hi,lo)
end

convert{T<:Integer}(::Type{Rational{T}}, x::Float120) =
    convert(Rational{T},convert(Rational{BigInt},convert(BigFloat,x)))

convert{T<:Integer}(::Type{Float120}, x::Rational{T}) =
   convert(Float120, convert(BigFloat,convert(Rational{BigInt},x)))

# promotions

promote_rule(::Type{Bool}, ::Type{Float120}) = Float120
promote_rule(::Type{Base.GMP.BigInt}, ::Type{Float120}) = Float120
promote_rule(::Type{Rational{Int128}}, ::Type{Float120}) = Rational{Int128}
promote_rule(::Type{Rational{Base.GMP.BigInt}}, ::Type{Float120}) = Rational{Base.GMP.BigInt}
promote_rule(::Type{BigFloat}, ::Type{Float120}) = BigFloat

promote_rule{T<:AbstractFloat}(::Type{T}, ::Type{Float120}) = Float120
promote_rule{T<:Integer}(::Type{T}, ::Type{Float120}) = Float120
promote_rule{T<:Integer}(::Type{Rational{T}}, ::Type{Float120}) = Float120