text:
## First pass at math functions. There are more domain restrictions to be implemented.

### SetPrecision

@sjseealso_group(N, SetPrecision, GetPrecison, BigFloatInput, BigIntInput, BI, BF)

@mkapprule SetPrecision :nargs => 1
@sjdoc SetPrecision """
    SetPrecison(n)

Set the precsion of BigFloat numbers to `n` decimal digits. If `N` does not give the result you
want, you can use `SetPrecision`.
"""
@doap function SetPrecision(n::Real)
    setprecision(round(Int,3.322*n))
    return n
end

### GetPrecision

@mkapprule GetPrecision :nargs => 0
@sjdoc GetPrecision """
    GetPrecision(n)

Get the precsion of BigFloat numbers.
"""
@doap function GetPrecision()
    return precision(BigFloat)
end

"""
    evalmath(expr::Expr)

Evaluate `expr`. This is used for writing Symata interfaces
to math functions. Currently, it is equivalent to `eval`.
"""
evalmath(expr::Expr) = Symata.eval(expr::Expr)

## We also need to use these to convert Symata expressions to Julia
## This function writes apprules for math functions. Usally dispatches floating point
## args to Julia functions. Some have a fall through to SymPy
## These tuples have 1,2, or 3 members. Symbols are for Julia, Symata, and SymPy
## If only one member is present, the second is constructed by putting an inital capital on the first.
## If only one or two members are present, we do not fall back to SymPy.
## The last list has tuples with 2 args for which we do not use any Julia function.

# Typical symbols for Julia, Symata, SymPy
function mtr(sym::Symbol)
    s = string(sym)
    sjf = uppercasefirst(s)
    (sym, Symbol(sjf), sym)
end

# We choose not to implement :cis. Julia has it. sympy and mma, not

# List of not yet implemented sympy functions from keys(sympy.functions)
# :E1,:Ei,:FallingFactorial,:Heaviside, :Id,
# :Max,:Min,:Piecewise,
#  :Ynm_c,:Znm,
# :__builtins__,:__doc__,:__file__,:__name__,:__package__,:__path__,
# :acos,:acosh,:acot,:acoth,:acsc,:adjoint,
#  :assoc_laguerre,:assoc_legendre,
#  :bspline_basis, :bspline_basis_set, :ceiling,:chebyshevt_root,
# :chebyshevu_root,:combinatorial
#  :elementary,
#  :exp_polar, :factorial2,:ff
# :floor, :hyper,:im, :jacobi_normalized, :jn, :jn_zeros,
# :periodic_argument,:piecewise_fold,:polar_lift,
# :principal_branch,:re,:real_root,:rf,:root, :special,
# :transpose,:unbranched_argument,:yn

# TODO: FIX: sympy zeta may take two args
# two arg both float or complex : zeta(s,z)  (with domain restrictions)
# BellB should be split into BellB and BellY

# TODO: Implement rewrite
# sympy.airyai(x)[:rewrite](sympy.hyper)
# PyObject -3**(2/3)*x*hyper((), (4/3,), x**3/9)/(3*gamma(1/3)) + 3**(1/3)*hyper((), (2/3,), x**3/9)/(3*gamma(2/3))

# TODO: rename these arrays. Name by number of args plus identifying number.
# ie. two_args1, two_args2. Then explain in comment the particulars. A naming
# scheme to capture all possibilities is to cumbersome.

#  mtr(:exp)  <-- trying with and without this
   const single_arg_float_complex =   # check, some of these can't take complex args
    [ mtr(:sin), mtr(:cos), mtr(:tan), (:sind,:SinD), (:cosd,:CosD),(:tand,:TanD),
         (:sinpi,:SinPi), (:cospi,:CosPi), mtr(:sinh), mtr(:cosh),
         mtr(:tanh), (:acos,:ArcCos,:acos), (:asin,:ArcSin,:asin),
         (:atan,:ArcTan,:atan),(:acosd,:ArcCosD), (:asind,:ArcSinD),
         (:atand,:ArcTanD), mtr(:sec), mtr(:csc), mtr(:cot),(:secd,:SecD),(:csc,:CscD),(:cot,:CotD),
         (:asec,:ArcSec),(:acsc,:ArcCsc),(:acot,:ArcCot),  # (:acotd,:ArcCotD)
         (:csch,), mtr(:coth,),(:asinh,:ASinh,:asinh), (:acosh,:ACosh, :acosh),(:atanh,:ATanh, :atanh),
         (:asech,:ArcSech),(:acsch,:ArcCsch),(:acoth,:ArcCoth, :acoth),
         (:cosc,),
         (:log1p,),(:exp2,),(:exp10,),(:expm1,),(:abs2,),
         mtr(:erfc), mtr(:erfi),(:erfcx,),(:dawson,),(:real,:Re, :re),(:imag,:Im, :im),
         (:angle,:Arg,:arg),  (:lgamma, :LogGamma, :loggamma),
         (:lfact,:LogFactorial), mtr(:digamma), mtr(:trigamma),
         (:airyai,:AiryAi,:airyai),
         (:airybi,:AiryBi,:airybi),(:airyaiprime,:AiryAiPrime,:airyaiprime),(:airybiprime,:AiryBiPrime,:airybiprime),
         (:besselj0,:BesselJ0),(:besselj1,:BesselJ1),(:bessely0,:BesselY0),(:bessely1,:BesselY1),
         (:zeta,:Zeta,:zeta)
#         (:eta,:DirichletEta,:dirichlet_eta)
         ]

# (:log,),   removed from list above, because it must be treated specially (and others probably too!)

const single_arg_float_int_complex = [(:sinc,)]

const single_arg_float = [(:erfcinv,:InverseErfc,:erfcinv),(:invdigamma,:InverseDigamma)]

# FIXME  Ceiling should probably return integer type. So it can't be included in generic code here (or punt to sympy for everything)
    const single_arg_float_int = [(:signbit,:SignBit), (:ceil, :Ceiling, :ceiling), (:floor, :Floor, :floor)]

#  (:factorial,:Factorial, :factorial), <--- we do this by hand

    const single_arg_int = [(:isqrt,:ISqrt),(:ispow2,:IsPow2),(:nextpow2,:NextPow2),(:prevpow2,:PrevPow2),
                      (:isprime,:PrimeQ), (:(~), :BitNot)
                        ]

    const two_arg_int = [(:binomial,:Binomial,:binomial)
                         ]

# Do NDigits by hand for now!
@mkapprule NDigits :nargs => 1:2
@doap NDigits(n::T,b::V) where {T<:Integer,V<:Integer} = ndigits(n, base=b)
@doap NDigits(n::T) where {T<:Integer} = ndigits(n)
@doap NDigits(n) = mx
@doap NDigits(n,b) = mx

# TODO
#  Two or more args
#   (:gcd,:GCD), (:lcm, :LCM), Also sympy does these for polynomials

# Complicated!
    const one_or_two_args1 = [(:polygamma, :PolyGamma, :polygamma)]

    const two_arg_float_and_float_or_complex =
     [
      (:besselj,:BesselJ, :besselj), (:besseljx,:BesselJx), (:bessely,:BesselY,:bessely),
      (:besselyx,:BesselYx), (:hankelh1,:HankelH1,:hankel1), (:hankelh1x,:HankelH1x),
      (:hankelh2,:HankelH2,:hankel2), (:hankelh2x,:HankelH2x), (:besseli,:BesselI, :besseli),
      (:besselix,:BesselIx), (:besselk,:BesselK,:besselk), (:besselkx,:BesselKx),
       (:atan2,:ArcTan2,:atan2)
      ]

    const  two_arg_float = [ (:beta,),(:lbeta,:LogBeta),(:hypot,)]

    ## There are no Julia functions for these (or at least we are not using them).
    ## First symbol is for  Symata, Second is SymPy

# TODO: Migrate functions out this non-specific list
    const no_julia_function = [  (:ExpIntegralE, :expint),
                           (:GegenbauerC, :gegenbauer), (:HermiteH, :hermite),
                           (:BellB, :bell), (:BernoulliB, :bernoulli), (:CatalanNumber, :catalan),
                           (:EulerE, :euler), (:Subfactorial, :subfactorial), (:Factorial2, :factorial2),
                           (:FactorialPower, :FallingFactorial), (:Pochhammer, :RisingFactorial),
 (:LeviCivita, :LeviCivita),
(:KroneckerDelta, :KroneckerDelta),  (:HypergeometricPFQ, :hyper), (:FactorSquareFree, :sqf), (:MellinTransform, :mellin_transform),
(:InverseMellinTransform, :inverse_mellin_transform), (:SineTransform, :sine_transform), (:InverseSineTransform, :inverse_sine_transform),
(:CosineTransform, :cosine_transform), (:InverseCosineTransform, :inverse_cosine_transform),
(:HankelTransform, :hankel_transform), (:InverseHankelTransform, :inverse_hankel_transform)
]
# LeviCivita is different than in mathematica

# Mobius needs filtering
const no_julia_function_one_arg = [ (:EllipticK, :elliptic_k), (:HeavisideTheta,:Heaviside),
                                    (:CosIntegral, :Ci), (:SinIntegral, :Si),
                                    (:FresnelC, :fresnelc), (:FresnelS, :fresnels), (:DiracDelta, :DiracDelta),
                                    (:LogIntegral, :Li), (:Ei, :Ei), (:ExpandFunc, :expand_func), (:Denominator, :denom),
                                    (:MoebiusMu, :mobius), (:EulerPhi, :totient), (:Divisors, :divisors), (:DivisorCount, :divisor_count),
                                     (:SinhIntegral, :Shi),(:CoshIntegral, :Chi) ]

# elliptic_e take more than ints!
const no_julia_function_one_or_two_int = [ (:HarmonicNumber, :harmonic) , (:EllipticE, :elliptic_e) ]

# FIXME: DivisorSigma has arg order reversed in Mma
const no_julia_function_two_args = [(:LegendreP, :legendre_poly), (:EllipticF, :elliptic_f),
                                    (:ChebyshevT, :chebyshevt), (:ChebyshevU, :chebyshevu),
                                    (:Cyclotomic, :cyclotomic_poly), (:SwinnertonDyer, :swinnerton_dyer_poly),
                                    (:SphericalBesselJ, :jn), (:SphericalBesselY, :yn), (:DivisorSigma, :divisor_sigma)
#                                    (:PolyLog, :polylog), (:SphericalBesselJ, :jn), (:SphericalBesselY, :yn), (:DivisorSigma, :divisor_sigma)
                                    ]

const no_julia_function_two_or_three_args = [ (:EllipticPi, :elliptic_pi), (:LaguerreL, :laguerre_poly)]

const no_julia_function_three_args = [ (:LerchPhi, :lerchphi) ]

const no_julia_function_four_args = [ (:JacobiP, :jacobi), (:SphericalHarmonicY, :Ynm) ]

function make_math()

# Note: in Mma and Julia, catalan and Catalan are Catalan's constant. In sympy catalan is the catalan number
# can't find :stirling in sympy

    for x in no_julia_function_one_arg
        sjf = x[1]
        eval(macroexpand(Symata, :( @mkapprule $sjf)))
        write_sympy_apprule(x[1],x[2],1)
    end

    for x in no_julia_function_two_args
        sjf = x[1]
        eval(macroexpand(Symata, :( @mkapprule $sjf :nargs => 2)))
        write_sympy_apprule(x[1],x[2],2)
    end

    for x in no_julia_function_four_args
        nargs = 4
        sjf = x[1]
        eval(macroexpand(Symata, :( @mkapprule $sjf :nargs => $nargs)))
        write_sympy_apprule(x[1],x[2],nargs)
    end

    for x in no_julia_function_three_args
        nargs = 3
        sjf = x[1]
        eval(macroexpand(Symata, :( @mkapprule $sjf :nargs => $nargs)))
        write_sympy_apprule(x[1],x[2],nargs)
    end

    for x in no_julia_function_two_or_three_args
        sjf = x[1]
        eval(macroexpand(Symata, :( @mkapprule $sjf)))
        write_sympy_apprule(x[1],x[2],2)
        write_sympy_apprule(x[1],x[2],3)
    end

    for x in no_julia_function_one_or_two_int
        sjf = x[1]
        eval(macroexpand(Symata, :( @mkapprule $sjf)))
        write_sympy_apprule(x[1],x[2],1)
        write_sympy_apprule(x[1],x[2],2)
    end

    # TODO: update this code
    for x in no_julia_function
        set_up_sympy_default(x...)
        clear_attributes(x[1]) ## FIXME: Listable was set in previous line. Now we unset it. We need a better system for this.
        set_sysattributes(x[1])
    end

    # Ok, this works. We need to clean it up
    for x in single_arg_float_complex
        jf,sjf = only_get_sjstr(x...)
        eval(macroexpand(Symata, :( @mkapprule $sjf :nargs => 1 )))
        write_julia_numeric_rule(jf,sjf,"AbstractFloat")
        write_julia_numeric_rule(jf,sjf,"CAbstractFloat")
        if length(x) == 3
            write_sympy_apprule(sjf,x[3],1)
        end
        set_attribute(Symbol(sjf),:Listable)
    end

    for x in single_arg_float
        jf,sjf = only_get_sjstr(x...)
        eval(macroexpand(Symata, :( @mkapprule $sjf :nargs => 1 )))
        write_julia_numeric_rule(jf,sjf,"AbstractFloat")
        if length(x) == 3
            write_sympy_apprule(sjf,x[3],1)
        end
        set_attribute(Symbol(sjf),:Listable)
    end

    for x in single_arg_float_int
        jf,sjf = only_get_sjstr(x...)
        eval(macroexpand(Symata, :( @mkapprule $sjf :nargs => 1 )))
        write_julia_numeric_rule(jf,sjf,"Real")
        if length(x) == 3
            write_sympy_apprule(sjf,x[3],1)
        end
    end

    # This is all numbers, I suppose
    for x in single_arg_float_int_complex
        jf,sjf = only_get_sjstr(x...)
        eval(macroexpand(Symata, :( @mkapprule $sjf :nargs => 1 )))
        write_julia_numeric_rule(jf,sjf,"Real")
        write_julia_numeric_rule(jf,sjf,"CReal")
        if length(x) == 3
            write_sympy_apprule(sjf,x[3],1)
        end
    end

    for x in single_arg_int
        jf,sjf = only_get_sjstr(x...)
        eval(macroexpand(Symata, :( @mkapprule $sjf :nargs => 1 :argtypes => [Integer] )))
        write_julia_numeric_rule(jf,sjf,"Integer")
        if length(x) == 3
            write_sympy_apprule(sjf,x[3],1)
        end
    end

    for x in two_arg_int
        jf,sjf = only_get_sjstr(x...)
        eval(macroexpand(Symata, :( @mkapprule $sjf :nargs => 2 )))
        write_julia_numeric_rule(jf,sjf,"Integer", "Integer")
        if length(x) == 3
            write_sympy_apprule(sjf,x[3],2)
        end
    end

    # Mma allows one arg, as well
    for x in one_or_two_args1
        jf,sjf = only_get_sjstr(x...)
        eval(macroexpand(Symata, :( @mkapprule $sjf )))
        write_julia_numeric_rule(jf,sjf,"Integer", "AbstractFloat")
        write_julia_numeric_rule(jf,sjf,"Integer", "CAbstractFloat")
        if length(x) == 3
            write_sympy_apprule(sjf,x[3],2)
        end
    end

    for x in two_arg_float
        jf,sjf = only_get_sjstr(x...)
        eval(macroexpand(Symata, :( @mkapprule $sjf :nargs => 2 )))
        write_julia_numeric_rule(jf,sjf,"AbstractFloat", "AbstractFloat")
    end

    for x in two_arg_float_and_float_or_complex
        jf,sjf = only_get_sjstr(x...)
        eval(macroexpand(Symata, :( @mkapprule $sjf :nargs => 2 )))
        write_julia_numeric_rule(jf,sjf,"AbstractFloat", "AbstractFloat")
        write_julia_numeric_rule(jf,sjf,"AbstractFloat", "CAbstractFloat")
        write_julia_numeric_rule(jf,sjf,"Integer", "AbstractFloat")
        write_julia_numeric_rule(jf,sjf,"Integer", "CAbstractFloat")
        if length(x) == 3
            write_sympy_apprule(sjf,x[3],2)
        end
    end
end

Typei(i::Int)  =  "T" *string(i)
complex_type(t::AbstractString) = "Complex{" * t * "}"
function write_julia_numeric_rule(jf, sjf, types...)
    annot = join(AbstractString[ "T" * string(i) * "<:" *
                                 (types[i][1] == 'C' ? types[i][2:end] : types[i])  for i in 1:length(types)], ", ")
    protargs = join(AbstractString[ "x" * string(i) * "::" *
                                (types[i][1] == 'C' ? complex_type(Typei(i)) : Typei(i))
                                    for i in 1:length(types)], ", ")
    callargs = join(AbstractString[ "x" * string(i) for i in 1:length(types)], ", ")
#    appstr = "do_$sjf{$annot}(mx::Mxpr{:$sjf},$protargs) = $jf($callargs)"
    appstr = "do_$sjf(mx::Mxpr{:$sjf},$protargs) where {$annot} = $jf($callargs)"
    eval(Meta.parse(appstr))
end

function only_get_sjstr(jf,sjf,args...)
    return jf, sjf
end

function only_get_sjstr(jf)
    st = string(jf)
    sjf = Symbol(uppercasefirst(st))
    return jf, sjf
end

function get_sjstr(jf,sjf)
    do_common(sjf)
    return jf, sjf
end

function get_sjstr(jf)
    st = string(jf)
    sjf = String(uppercasefirst(st))
    do_common(sjf)
    return jf, sjf
end

## FIXME: This is outdated. Some of this is handled in @mkapprule
# Handle functions that do *not* fall back on SymPy
function do_common(sjf)
    aprs = "Symata.apprules(mx::Mxpr{:$sjf}) = do_$sjf(mx,margs(mx)...)"
    aprs1 = "do_$sjf(mx::Mxpr{:$sjf},x...) = mx"
    evalmath(Meta.parse(aprs))
    evalmath(Meta.parse(aprs1))
    set_attribute(Symbol(sjf),:Protected)
    set_attribute(Symbol(sjf),:Listable)
end

# Handle functions that fall back on SymPy
function get_sjstr(jf, sjf, sympyf)
    set_up_sympy_default(sjf, sympyf)
    return jf, sjf
end

# Faster if we don't do interpolation
function write_sympy_apprule(sjf, sympyf, nargs::Int)
    callargs = Array{AbstractString}(undef, 0)
    sympyargs = Array{AbstractString}(undef, 0)
    for i in 1:nargs
        xi = "x" * string(i)
        push!(callargs, xi)
        push!(sympyargs, "sjtopy(" * xi * ")")
        end
    cstr = join(callargs, ", ")
    sstr = join(sympyargs, ", ")
    aprpy = "function do_$sjf(mx::Mxpr{:$sjf},$cstr)
               try
                 (sympy[:$sympyf]($sstr) |> pytosj)
               catch e
                 showerror(stdout, e)
                 mx
               end
            end"
    evalmath(Meta.parse(aprpy))
end

function set_up_sympy_default(sjf, sympyf)
    aprs = "Symata.apprules(mx::Mxpr{:$sjf}) = do_$sjf(mx,margs(mx)...)"
    aprs1 = "function do_$sjf(mx::Mxpr{:$sjf},x...)
               try
                 (sympy[:$sympyf](map(sjtopy,x)...) |> pytosj)
               catch
                   mx
               end
           end"
    evalmath(Meta.parse(aprs))
    evalmath(Meta.parse(aprs1))
    set_attribute(Symbol(sjf),:Protected)
    set_attribute(Symbol(sjf),:Listable)
end

make_math()

#### Gamma function

@mkapprule Gamma :nargs => 1:2

@doap Gamma(a) = sjgamma(a)
@doap Gamma(a,z) = sjgamma(a,z)

set_pytosj(:gamma, :Gamma)
set_pytosj(:uppergamma, :Gamma)
set_sjtopy(:Gamma, :sympy_gamma)

@sjdoc Gamma """
    Gamma(a)

the gamma function.

    Gamma(a,z)

the upper incomplete Gamma function.
"""

### GammaRegularized

## Implementation is not complete
@mkapprule GammaRegularized :nargs => 2:3

@doap function GammaRegularized(a,z::Number)
    if z == 0
       return  _gammaregularizedzero(mx,a)
    end
   return  _gammaregularizednozero(mx,a,z)
end

function _gammaregularizedzero(mx,a::Real)
    a > 0 ? 1 : a == 0 ? 0 : ComplexInfinity
end

function _gammaregularizedzero(mx,a::Complex)
    ra = real(a)
    ra > 0 ? 1 : ra == 0 ? Indeterminate : ComplexInfinity
end

function _gammaregularizedzero(mx,a)
    mx
end

function _gammaregularizednozero(mx,a::Real,z::Number)
    a < 0 && return zero(a) ## need to consider type of z as well
    mxpr(:Gamma,a,z)/mxpr(:Gamma,a)
end

function _gammaregularizednozero(mx,a::Number,z::Number)
    mxpr(:Gamma,a,z)/mxpr(:Gamma,a)
end

function _gammaregularizednozero(mx,a,z)
    mx
end

### Erf

@mkapprule Erf :nargs => 1:2

@doap Erf(x) = sjerf(x)
@doap Erf(x,y) = sjerf(x,y)

set_pytosj(:erf, :Erf)
set_pytosj(:erf2, :Erf)
set_sjtopy(:Erf, :sympy_erf)

### InverseErf

@mkapprule InverseErf :nargs => 1:2

@doap InverseErf(x::T) where {T<:AbstractFloat} = x < 1 && x > -1 ? erfinv(x) : mx
@doap InverseErf(x) = pytosj(sympy[:erfinv](sjtopy(x)))
@doap InverseErf(x,y) = pytosj(sympy[:erf2inv](sjtopy(x),sjtopy(y)))

register_sjfunc_pyfunc(:InverseErf,:erfinv)
register_only_pyfunc_to_sjfunc(:InverseErf,:erf2inv)

### IntegerDigits

@sjdoc IntegerDigits """
    IntegerDigits(n,[, base][, pad])

return an array of the digits of `n` in the given `base`,
optionally padded with zeros to `pad` characters.

In contrast to Julia, more significant digits are at lower indexes.
"""
@mkapprule IntegerDigits  :nargs => 1:3
set_sysattributes(:IntegerDigits, :Listable)

_intdig(args) = MList(reverse!(args)...) |> setfixed

@doap IntegerDigits(n::Integer) = digits(n) |> _intdig
@doap IntegerDigits(n::Integer,b::Integer) = digits(n, base=convert(Int, b)) |> _intdig
@doap IntegerDigits(n::Integer,b::Integer,p::Integer) = digits(n, base=convert(Int, b), convert(Int, p)) |> _intdig

@mkapprule Log :nargs => 1:2

# sjlog defined in wrappers.jl
@doap Log(x,y) = sjlog(x,y)
@doap Log(x) = sjlog(x)

### N

@sjdoc N """
    N(expr)

try to give a the numerical value of `expr`.

    N(expr,p)

try to give `p` decimal digits of precision.

Sometimes `N` does not give the number of digits requested. In this case, you can use `SetPrecision`.
"""

## TODO: call evalf(expr,n) on sympy functions to get arbitrary precision numbers.
## Eg. N(LogIntegral(4),30) does not give correct result.
## ... a list of heads all of which should be converted and then passed to evalf.
## N needs to be rewritten
function apprules(mx::Mxpr{:N})
    outer_N(margs(mx)...)
end

function outer_N(expr)
    do_N(expr)
end

function outer_N(expr, p)
    if p > 16
        pr = precision(BigFloat)
        dig = round(Int,p*3.322)
        set_mpmath_dps(dig)     # This is for some SymPy math. But, it is not clear when this will help us
        setprecision( BigFloat, dig) # new form
        res = meval(do_N(expr,p))
        setprecision(BigFloat, pr)
        restore_mpmath_dps()
        return res
    else
        do_N(expr)
    end
end

do_N(x::Float64,dig) = x  # Mma does this
do_N(x,dig) = x
do_N(x) = x
do_N(n::T) where {T<:AbstractFloat} = n
do_N(n::T) where {T<:Real} = float(n)
do_N(n::Complex{T}) where {T<:Real} = complex(float(real(n)), float(imag(n)))
do_N(n::T,p::V) where {T<:Real,V<:Integer}  = float_with_precision(n,p)

# p is number of decimal digits of precision
# Julia doc says set_bigfloat_precision uses
# binary digits, but it looks more like decimal.

function float_with_precision(x,p)
        res = BigFloat(x)
end

# These rely on fixed-point evaluation to continue with N, this is not efficient
# We need to do it all here.

function _make_N_body(d1,d2)
    quote begin
        get_attribute(mx,:NHoldAll) && return mx;
        len = length(mx)
        args = margs(mx)
        nargs = newargs(len)
        start::Int = 1
        if get_attribute(mx,:NHoldFirst)
            nargs[1] = args[1]
            for i in 2:len
                nargs[i] = $d2
            end
        elseif get_attribute(mx,:NHoldRest)
            nargs[1] = $d1
            for i in 2:len
                nargs[i] = args[i]
            end
        else
            @inbounds for i in 1:len
                nargs[i] = $d2
            end
        end
        res = mxpr(mhead(mx),nargs)
    end
  end
end

function make_Mxpr_N()
    _Nbody = _make_N_body(:(do_N(args[1])), :(do_N(args[i])) )
    @eval begin
        function do_N(mx::Mxpr)
            $_Nbody
            res
        end
    end

    _Nbody = _make_N_body(:(do_N(args[1],p)), :(do_N(args[i],p)) )
    @eval begin
        function do_N(mx::Mxpr,p::T) where T<:Integer
            $_Nbody
            res
        end
    end
end

make_Mxpr_N()

@sjdoc GoldenRatio """
    GoldenRatio

equal to `(1+Sqrt(5))/2`.
"""

# We need to use dispatch as well, not conditionals
function do_N(s::SJSym)
    if s == :Pi || s == :π
        return float(MathConstants.pi)
    elseif s == :E
        return float(MathConstants.e)
    elseif s == :EulerGamma
        return float(MathConstants.eulergamma)
    elseif s == :GoldenRatio
        return float(MathConstants.golden)
    elseif s == :Catalan
        return float(MathConstants.catalan)
    end
    return s
end

function do_N(s::SJSym,pr::T) where T<:Integer
    if s == :Pi || s == :π
        return float_with_precision(MathConstants.pi, pr)
    elseif s == :E
        return float_with_precision(MathConstants.e, pr)
    elseif s == :EulerGamma
        return float_with_precision(MathConstants.eulergamma, pr)
    elseif s == :GoldenRatio
        return float_with_precision(MathConstants.golden, pr)
    elseif s == :Catalan
        return float_with_precision(MathConstants.catalan, pr)
    end
    return s
end

@sjdoc Precision """
    Precision(x)

get the precision of a floating point number `x` in decimal digits

    Precision()

get the current precision of `BigFloat` arithmetic

    Precision(Binary => True)
    Precision(x, Binary => True)

give precision as defined by the effective number of bits in the mantissa.
"""

@mkapprule Precision :nodefault => true, :options => Dict( :Binary => false )

#apprules(mx::Mxpr{:Precision}) = do_Precision(mx,margs(mx)...)

# TODO: mkapprule should write the correct default rule if keywords are present
do_Precision(mx::Mxpr{:Precision},args...; kws...) = mx

function do_Precision(mx::Mxpr{:Precision},x::AbstractFloat; Binary=false)
    Binary ? precision(x) : round(Int,precision(x) * log10(2))
end

function do_Precision(mx::Mxpr{:Precision}; Binary=false)
    pr = precision(BigFloat)
    Binary ? pr : round(Int,pr * log10(2))
end

@sjdoc Re """
    Re(z)

return the real part of `z`.
"""

@sjdoc Im """
    Im(z)

return the imaginary part of `z`.
"""

# Using Julia would be faster in some cases for ReIm and AbsArg
@sjdoc ReIm """
    ReIm(z)

returns `[Re(z),Im(z)]`.
"""

@mkapprule ReIm
@doap ReIm(z) = mxpr(:List,mxpr(:Re,z), mxpr(:Im,z))

# FIXME Abs(I-3) returns 10^(1/2) instead of (2^(1/2))*(5^(1/2)). i.e. does not evaluate to a fixed point
@sjdoc AbsArg """
    AbsArg(z)

return `[Abs(z),Arg(z)]`.
"""

@mkapprule AbsArg
@doap AbsArg(z) = mxpr(:List,mxpr(:Abs,z), mxpr(:Arg,z))

# Re and Im code below is disabled so it does not interfere with SymPy.
# Mma allows complex numbers of mixed Real type. So does Julia: Complex{Real}(1//3, 3.0)
# Implementation not complete. eg  Im(a + I *b) --> Im(a) + Re(b)

# do_Re{T<:Real}(mx::Mxpr{:Re}, x::Complex{T}) = real(x)
# do_Re(mx::Mxpr{:Re}, x::Real) = x

# function do_Re(mx::Mxpr{:Re}, m::Mxpr{:Times})
#     f = m[1]
#     return is_imaginary_integer(f) ? do_Re_imag_int(m,f) : mx
# end

# # dispatch on type of f. Maybe this is worth something.
# function do_Re_imag_int(m,f)
#     nargs = copy(margs(m))
#     popfirst!(nargs)
#     if length(nargs) == 1
#         return mxpr(:Times,-imag(f),mxpr(:Im,nargs))
#     else
#         return mxpr(:Times,-imag(f),mxpr(:Im,mxpr(:Times,nargs)))
#     end
# end

# do_Im{T<:Real}(mx::Mxpr{:Im}, x::Complex{T}) = imag(x)
# do_Im(mx::Mxpr{:Im}, x::Real) = zero(x)

# function do_Im(mx::Mxpr{:Im}, m::Mxpr{:Times})
#     f = m[1]
#     return is_imaginary_integer(f) ? do_Im_imag_int(m,f) : mx
# end

# function do_Im_imag_int(m,f)
#     nargs = copy(margs(m))
#     popfirst!(nargs)
#     if length(nargs) == 1
#         return mxpr(:Times,imag(f),mxpr(:Re,nargs))
#     else
#         return mxpr(:Times,imag(f),mxpr(:Re,mxpr(:Times,nargs)))
#     end
# end

#### Complex

@sjdoc Complex """
    Complex(a,b)

return a complex number when `a` and `b` are Reals

The conversion is done when the expression is parsed, so it is much faster than `a + I*b`.
"""

# Complex with two numerical arguments is converted at parse time. But, the
# arguments may evaluate to numbers only at run time, so this is needed.
# mkapprule requires that the first parameter do_Complex be annotated with the Mxpr type.

@mkapprule Complex
do_Complex(mx::Mxpr{:Complex},a::T,b::T) where {T<:Number} = complex(a,b)
do_Complex(mx::Mxpr{:Complex},a::T) where {T<:Number} = complex(a)

@sjdoc Rational """
    Rational(a,b), or a/b

return a `Rational` for `Integer`s `a` and `b`.  This is done when the
expression is parsed, so it is much faster than `a/b`.
"""

# Same here. But we need to use mdiv to reduce rationals to ints if possible.
@mkapprule Rational
do_Rational(mx::Mxpr{:Rational},a::Number,b::Number) = mdiv(a,b)

#### Rationalize

@sjdoc Rationalize """
    Rationalize(x)

return a `Rational` approximation of `x`.

    Rationalize(x,tol)

return an approximation differing from `x` by no more than `tol`.
"""

@mkapprule Rationalize
do_Rationalize(mx::Mxpr{:Rationalize},x::AbstractFloat) = rationalize(x)
do_Rationalize(mx::Mxpr{:Rationalize},x::AbstractFloat,tol::Number) = rationalize(x,tol=float(tol))
function do_Rationalize(mx::Mxpr{:Rationalize},x::Symbolic)
    r = doeval(mxpr(:N,x))  # we need to redesign do_N so that we can call it directly. See above
    return isa(r,AbstractFloat) ? do_Rationalize(mx,r) : x
end
function do_Rationalize(mx::Mxpr{:Rationalize},x::Symbolic,tol::Number)
    ndig = round(Int,-log10(tol))      # This is not quite correct.
    r = doeval(mxpr(:N,x,ndig))  # we need to redesign do_N so that we can call it directly. See above.
    return isa(r,AbstractFloat) ? do_Rationalize(mx,r,tol) : x
end
do_Rationalize(mx::Mxpr{:Rationalize},x) = x

#### Numerator

@sjdoc Numerator """
    Numerator(expr)

return the numerator of `expr`.
"""

apprules(mx::Mxpr{:Numerator}) = do_Numerator(mx::Mxpr{:Numerator},margs(mx)...)
do_Numerator(mx::Mxpr{:Numerator},args...) = mx

do_Numerator(mx::Mxpr{:Numerator},x::Rational) = numerator(x)
do_Numerator(mx::Mxpr{:Numerator},x) = x
function do_Numerator(mx::Mxpr{:Numerator},x::Mxpr{:Power})
    find_numerator(x)
end

function do_Numerator(mx::Mxpr{:Numerator},m::Mxpr{:Times})
    nargsn = newargs()
    args = margs(m)
    for i in 1:length(args)
        arg = args[i]
        res = find_numerator(arg)
        if res != 1
            push!(nargsn,res)
        end
    end
    if isempty(nargsn)
        return 1  # which one ?
    end
    if length(nargsn) == 1
        return nargsn[1]
    end
    return mxpr(mhead(m),nargsn)
end

find_numerator(x::Rational) = numerator(x)
find_numerator(x::Mxpr{:Power}) = pow_sign(x,expt(x)) > 0 ? x : 1

function pow_sign(x, texpt::Mxpr{:Times})
    fac = texpt[1]
    return is_Number(fac) && fac < 0 ? -1 : 1
end

pow_sign(x, texpt::Number) = texpt < 0 ? -1 : 1
pow_sign(x, texpt) = 1
find_numerator(x) = x

### Chop

@mkapprule Chop  :nargs => 1:2

@sjdoc Chop """
    Chop(expr)

set small floating point numbers in `expr` to zero.

    Chop(expr,eps)

set floating point numbers smaller in magnitude than `eps` in `expr` to zero.
"""
do_Chop(mx::Mxpr{:Chop}, x) = zchop(x)
do_Chop(mx::Mxpr{:Chop}, x, zeps) = zchop(x,zeps)
const chop_eps = 1e-14
zchop(x::T, eps=chop_eps) where {T<:AbstractFloat} = abs(x) > eps ? x : 0  # don't follow abstract type
zchop(x::T, eps=chop_eps) where {T<:Number} = x
zchop(x::Complex{T}, eps=chop_eps) where {T<:AbstractFloat} = complex(zchop(real(x),eps),zchop(imag(x),eps))
zchop(a::T, eps=chop_eps) where {T<:AbstractArray} = (b = similar(a); for i in 1:length(a) b[i] = zchop(a[i],eps) end ; b)
zchop!(a::T, eps=chop_eps) where {T<:AbstractArray} = (for i in 1:length(a) a[i] = zchop(a[i],eps) end ; a)
zchop(x::Expr,eps=chop_eps) = Expr(x.head,zchop(x.args)...)
zchop(x) = x
zchop(x,eps) = x

function zchop(mx::T) where T<:Mxpr
    nargs = similar(margs(mx))
    for i in 1:length(nargs) nargs[i] = zchop(mx[i]) end
    mxpr(mhead(mx), nargs)
end

function zchop(mx::T,zeps) where T<:Mxpr
    nargs = similar(margs(mx))
    for i in 1:length(nargs) nargs[i] = zchop(mx[i],zeps) end
    mxpr(mhead(mx), nargs)
end

### Conjugate

@mkapprule Conjugate :nargs => 1

@doap Conjugate(z::Number) = conj(z)
@doap Conjugate(z::Mxpr{:Plus}) = mxpr(:Plus, mapmargs( x -> mxpr(:Conjugate, x), margs(z))...)
@doap Conjugate(z::Mxpr) = z |> sjtopy |> sympy[:conjugate] |> pytosj

### Exp

# The parser normally takes care of this,
# But, when converting expressions from Sympy, we get Exp, so we handle it here.
# TODO: don't throw away other args. (but, we hope this is only called with returns from Sympy)
apprules(mx::Mxpr{:Exp}) = mxpr(:Power,:E,mx[1])

# TODO: handle 3rd argument
#### Mod
@mkapprule  Mod  :nargs => 2
do_Mod(mx::Mxpr{:Mod}, x::T, y::V) where {T<:Integer, V<:Integer} = mod1(x,y)

#### DivRem
@mkapprule  DivRem  :nargs => 2
do_DivRem(mx::Mxpr{:DivRem}, x::T, y::V) where {T<:Integer, V<:Integer} = mxprcf(:List,divrem(x,y)...)

#### Abs

@sjdoc Abs """
    Abs(z)

the absolute value of `z`.
"""

@mkapprule Abs :nargs => 1

@doap Abs(n::T) where {T<:Number} = mabs(n)

@doap function Abs(x)
    x |> sjtopy |> sympy[:Abs] |> pytosj
end

# The following code may be useful, but it is not complete
# # Abs(x^n) --> Abs(x)^n  for Real n
function do_Abs(mx::Mxpr{:Abs}, pow::Mxpr{:Power})
     doabs_pow(mx,base(pow),expt(pow))
end

# SymPy does not do this one
doabs_pow(mx,b,e::T) where {T<:Real} = mxpr(:Power,mxpr(:Abs,b),e)

function doabs_pow(mx,b,e)
    mx[1] |> sjtopy |> sympy[:Abs] |> pytosj
end

# do_Abs(mx::Mxpr{:Abs},prod::Mxpr{:Times}) = do_Abs(mx,prod,prod[1])

# #doabs(mx,prod,s::Symbol)

# function do_Abs{T<:Number}(mx::Mxpr{:Abs},prod,f::T)
#     f >=0 && return mx
#     if f == -1
#         return doabsmone(mx,prod,f)
#     end
#     args = copy(margs(prod))
#     args[1] = -args[1]
#     return mxpr(:Abs,mxpr(:Times,args))
# end

# function doabsmone{T<:Integer}(mx,prod,f::T)
#     args = copy(margs(prod))
#     popfirst!(args)
#     if length(args) == 1
#         return mxpr(:Abs,args)
#     else
#         return mxpr(:Abs,mxpr(:Times,args))
#     end
# end

# TODO Fix canonical routines so that 1.0 * a is not simplifed to a
function doabsmone(mx,prod,f::T) where T<:Real
    args = copy(margs(prod))
    popfirst!(args)
    if length(args) == 1
        res = mxpr(:Times,one(f),mxpr(:Abs,args))
    else
        res = mxpr(:Times,one(f),mxpr(:Abs,mxpr(:Times,args)))
    end
    return res
end

#### Sign

@mkapprule Sign :nargs => 1

do_Sign(mx::Mxpr{:Sign}, x::Number) = sign_number(x)

function sign_number(z::Complex{T}) where T<:Real
    av = mpow(real(z*conj(z)),-1//2)
    av == 1 ? z : mxprcf(:Times, z,  av)
end

sign_number(x::Real) = convert(Int,sign(x))

function do_Sign(mx::Mxpr{:Sign}, x)
    x |> sjtopy |> sympy[:sign] |> pytosj
end

# SymPy does not do Sign((I+1)*a), so we catch this case here.
function do_Sign(mx::Mxpr{:Sign}, x::Mxpr{:Times})
    length(x) == 0 && return 1
    x1 = x[1]
    if isa(x1,Number)
        mxpr(:Times, sign_number(x1), mxpr(:Sign , mxpr(:Times,x[2:end]...)))
    else
        x |> sjtopy |> sympy[:sign] |> pytosj
    end
end

# do_Sign(mx::Mxpr{:Sign}, x) = mx

# # Maybe each symbol should be a type. Like irrational or Mxpr
# function do_Sign(mx::Mxpr{:Sign}, x::SJSym)
#     x == :Pi && return 1
#     x == :E && return 1
#     x == :EulerGamma && return 1
#     return mx
# end

#### LowerGamma

@mkapprule LowerGamma

do_LowerGamma(mx::Mxpr{:LowerGamma}, a, z) =  mxpr(:Gamma,a) - mxpr(:Gamma,a,z)

@sjdoc LowerGamma """
    LowerGamma(a,z)

the lower incomplete gamma function.
"""

#### Norm

# TODO: Implement p norm
@mkapprule Norm :nargs => 1

@doap Norm(x::Number) = Abs(x)

@doap function Norm(x::Mxpr{:List})
    vectorq(x) || return mx # TODO error message
    local sum0 = 0
    for i in 1:length(x)
        sum0 += mpow(Abs(x[i]),2)
    end
    mpow(sum0,1//2)
end

@sjdoc ExpandFunc """
    ExpandFunc(expr)

rewrite some multi-parameter special functions using simpler functions.
"""

@sjdoc I """
    I

the imaginary unit.
"""

@sjdoc E """
    E

the base of the natural logarithm.
"""

@sjdoc Pi """
    Pi

the trigonometric constant `π`.
"""

@sjdoc BellB """
    BellB(n)

gives the `n^{th}` Bell number, `B_n`.

    BellB(n, x)

gives the `n^{th}` Bell polynomial, `B_n(x)`.

    BellB(n, k, [x1, x2, ... x_{n-k+1}])
"""

### Identity

@sjdoc Identity """
    Identity(expr)

return `expr`.
"""

@mkapprule Identity

@doap Identity(x) = x

### Total

@mkapprule Total

@sjdoc Total """
    Total(list)

give the sum of elements in `list`
"""

@doap Total(x::AbstractArray) = sum(x)

## We could try some things to speed this up... check that all elements are numbers ?
## with all() it takes about as much time to check as it does to sum. Maybe use try/catch
#@doap Total(x::Mxpr{:List}) = mxpr(:Plus,margs(x))
@doap Total(x::Mxpr{:List}) = _apply_plus(x)

### CosPi, SinPi

## for non-float arguments
## Not sure it is good to do this: N(CosPi(100)) will be less accurate than CosPi(N(100))
@doap CosPi(x) = mxpr(:Cos, mmul(x,Pi))
@doap SinPi(x) = mxpr(:Sin, mmul(x,Pi))

@sjdoc CosPi """
    CosPi(x)

returns `cos(π x)`. This is more accurate than `Cos(Pi*x)` for
floating point `x`.
"""

@sjdoc SinPi """
    SinPi(x)

returns `sin(π x)`. This is more accurate than `Sin(Pi*x)` for
floating point `x`.
"""

### ProductLog

@mkapprule ProductLog :nargs => 1:2

@doap function ProductLog(z::FloatRC)
    LambertW.lambertw(z)
end

@doap function ProductLog(z)
    pytosj(sympy[:LambertW](sjtopy(z)))
end

@doap function ProductLog(z::FloatRC, k::Number)
    LambertW.lambertw(z,k)
end

@doap function ProductLog(z,k)
    pytosj(sympy[:LambertW](sjtopy(z), sjtopy(k)))
end

### DirichletEta

## Note: sympy does not automatically transform eta to a function of zeta.
## However, Mma transforms it immediately
## However, if we do ex = DirichletEta(s) and then substitute 1 for s, we get indeterminate.
## So, one could argue for only converting DirichletEta(s) explicitly

@sjdoc DirichletEta """
    DirichletEta(s)

is the Dirichlet eta function.
"""

@mkapprule DirichletEta  :nargs => 1
@doap DirichletEta(x::FloatRC) = SpecialFunctions.eta(x)
@doap DirichletEta(x::Number) = _dirichlet_eta(x)
@doap DirichletEta(x::Integer) = x == 1 ? Log2 : _dirichlet_eta(x)
@doap DirichletEta(x) = _dirichlet_eta(x)
_dirichlet_eta(x::Number) = (1 - mpow(2, (1 - x))) * Zeta(x)
_dirichlet_eta(x) = mminus(1, mpow(2, mminus(1, x))) * Zeta(x)

### PolyLog

@sjdoc PolyLog """
    PolyLog(s,z)

is the polylogarithm function.
"""

@mkapprule PolyLog :nargs => 2

@doap PolyLog(s, z) = _polylog(mx,s,z)

_polyzero(z) = mmul(z , mpow(mminus(1,z),-1))

function _polylog(mx,s::Integer,z)
    s == 1 && return -Log(mminus(1,z))
    s == 0 && return _polyzero(z)
    if s < 0
        u = gensym()
        r = _polyzero(u)
        for i in 1:(-s)
            r = doeval(mxpr(:Factor, u*mxpr(:D,r,u)))
        end
        return replaceall(r, mxpr(:Rule, u, z))
    end
    if z == 1//2
        if s == 2
            return (Pisq - 6*(mpow(Log2,2)))/12
        elseif s == 3
            return (4*Log2^3  - 2*Pisq*Log2 + 21*Zeta(3))/24
        end
    end
    if s == 2
        z == I && return mmul(I,:Catalan) - Pisq/48
        z == MinusI && return -mmul(I,:Catalan) - Pisq/48
    end
    _polylog1(mx,s,z)
end

_polylog(mx,s,z) = _polylog1(mx,s,z)

function _polylog1(mx,s,z::Integer)
    z == 1 && return Zeta(s)
    z == -1 && return -DirichletEta(s)
    _polylog2(mx,s,z)
end

_polylog1(mx,s,z) = _polylog2(mx,s,z)

function _polylog2(mx,s1::Number,z1::Number)
    (s,z) = promote(s1,z1)
    polylognum(mx,s,z)
end

_polylog2(mx,s,z) = mx

polylognum(mx,s::FloatRC,z::FloatRC) = sympy[:polylog](s,z) |> pytosj
polylognum(mx,s,z) = mx

### ExpPolar

## sympy returns this sometimes. I am not sure what it is, but it usually just hinders
## desired evaluation. We could also handle this as a direct translation in sympy.jl.
## There may be some reason to preserve this.

@mkapprule ExpPolar
@doap ExpPolar(x) = Exp(x)

### Factorial

@sjdoc Factorial """
    Factorial(n)

is the factorial function.
"""
@mkapprule Factorial :nargs => 1
@doap Factorial(n::SJSym) = mx
@doap Factorial(n) = sjfactorial(n)

# function sjfactorial(mx::Mxpr{:DirectedInfinity})
#     _sjfactorial_inf(mx,mx[1])
# end

_sjfactorial_inf(mx, x, v) = mx
_sjfactorial_inf(mx,x,v::Real) = v > 0 ? Infinity : Indeterminate
_sjfactorial_inf(mx,x,v::Complex) = real(v) == 0 ? 0 : mx

@doap function Factorial(x::Mxpr{:DirectedInfinity})
    length(x) == 0 && return mx
    _sjfactorial_inf(mx,x,x[1])
end

### Sqrt

## we must have square root for when this is constructed programatically. eg. Map(Sqrt,[4])
## Sqrt is not always read and translated!
@mkapprule Sqrt :nargs => 1
@doap Sqrt(x) = mpow(x,1//2)

### CubeRoot

@mkapprule CubeRoot :nargs => 1

@doap CubeRoot(x::AbstractFloat) = cbrt(x)
@doap function CubeRoot(x::T) where T<:Union{Integer,Rational}
    x == 0 && return 0
    x > 0 && return mpow(x,1//3)
    return -mpow(-x,1//3)
end

@doap function CubeRoot(x::Complex)
    @symwarn("CubeRoot::preal: The parameter $x should be real valued.")
    return mx
end

@doap function CubeRoot(x::Mxpr)
    return mxpr(:Surd,x,3)
end

@doap function CubeRoot(x::Symbol)
    return mxpr(:Surd,x,3)
end

### Surd

@mkapprule Surd :nargs => 2

@doap function Surd(x,n::Integer)
    if iseven(n)
        @symwarn("Surd::noneg: Surd is not defined for even roots of negative values.")
        return mx
    end
    return _surd(mx,x,n)
end

function _surd(mx,x::Real,n)
    x >= 0 && return mpow(x,1//n)
    mminus(mpow(-x,1//n))
end

function _surd(mx,x::Complex,n)
    @symwarn("Surd::preal: The parameter I should be real valued")
    mx
end

_surd(mx,x,n) = mx

function _surd(mx,x::TimesT,n)
    length(x) < 2 && return mx
    isa(x[1],Real) && return mxpr(:Surd,x[1],n) * mxpr(:Surd, mxpr(:Times, x[2:end]...),n)
    mx
end

### Power

@mkapprule Power :nargs => 2

do_Power(mx::Mxpr{:Power}, b::Symbolic, n::Integer) = n == 1 ? b : n == 0 ? one(n) : mx

@doap Power(b::SJSym, expt) = b == :E ? dopowerE(mx, expt) : mx

@doap Power(b::SJSym, expt::Integer) = expt == 1 ? b : expt == 0 ? one(expt) :
    b == :E ? dopowerE(mx, expt) : mx

dopowerE(mx, expt::T) where {T<:AbstractFloat} = exp(expt)
dopowerE(mx, expt::Complex{T}) where {T<:AbstractFloat} = exp(expt)

function dopowerE(mx, expt)
    syexpt = sjtopy(expt)
    syres = sympy[:exp](syexpt)
    res = pytosj(syres)
    if is_Mxpr(res, :Exp)
        res[1] == expt && return mx
        return mxpr(:Power, :E, margs(res)...)
    end
    res
end

# Don't handle this yet.
do_Power(mx::Mxpr{:Power},   b::Complex{T},expt::Rational{V}) where {T<:Integer,V<:Integer} = mx
do_Power(mx::Mxpr{:Power},   b::Complex{T},expt::Complex{Rational{V}}) where {T<:Integer,V<:Integer} = mx

do_Power(mx::Mxpr{:Power},   b::T,expt::V) where {T<:Number,V<:Number} = mpow(b,expt)

# For some reason, we need this integer rule. For instance for (a^2)^2 --> a^4
do_Power(mx::Mxpr{:Power}, b::Mxpr{:Power}, exp::T) where {T<:Integer} = mpow(base(b), mmul(exp,expt(b)))
do_Power(   mx::Mxpr{:Power}, b::Mxpr{:Power}, exp::T) where {T<:Real} = mpow(base(b), mmul(exp,expt(b)))

do_Power(mx::Mxpr{:Power},   b::Mxpr{:Power}, exp) = is_Number(expt(b)) ? mpow(base(b), mmul(expt(b),exp)) : mx

do_Power(mx::Mxpr{:Power},   b::SJSym,expt::Complex{T}) where {T<:AbstractFloat} = b == :E ? exp(expt) : mx

do_Power(mx::Mxpr{:Power},   b::Mxpr{:DirectedInfinity},expt::T) where {T<:Integer} = mpow(b,expt)
do_Power(mx::Mxpr{:Power},   b::Mxpr{:DirectedInfinity},expt::T) where {T<:Number} = mpow(b,expt)

function do_Power(mx::Mxpr{:Power},b::T,expt::V) where {T<:Integer, V<:Rational}
    mpow(b,expt)
end

do_Power(mx,b,e) = mx

### // (Rational)

apprules(mx::Mxpr{://}) = makerat(mx,mx[1],mx[2])
makerat(mx::Mxpr{://},n::T,d::T) where {T<:Number} = n//d
makerat(mx,n,d) = mx

### BI

@sjdoc BI """
    BI(n)

convert the number `n` to type `BigInt`. Symata neither
detects integer overflow, nor automatically promote integers to `BigInt`.
But, a literal integer will automatically be given a large enough storage type without using
`BI`.
"""

@sjseealso_group(BI,BF,Big,BigIntInput)

### BF

@sjdoc BF """
    BF(n), or BF"n"

convert the number, or string `n` to a BigFloat. Symata neither
detects overflow, nor automatically promotes types from fixed to arbitrary precision. The form
`BF"n"` is more efficient, being a julia macro that converts the string "n" upon parsing.
"""

@sjdoc Big """
    Big(n)

convert `n` a maximum precision representation, typically
`BigInt` or `BigFloat`.
"""

apprules(mx::Mxpr{:BI}) = dobigint(mx,mx[1])
dobigint(mx,x) = mx
dobigint(mx,x::T) where {T<:Number} = BigInt(x)
dobigint(mx,x::T) where {T<:AbstractString} = Meta.parse(BigInt,x)

apprules(mx::Mxpr{:BF}) = dobigfloat(mx,mx[1])
dobigfloat(mx,x) = mx
dobigfloat(mx,x::T) where {T<:Number} = BigFloat(x)
dobigfloat(mx,x::T) where {T<:AbstractString} = Meta.parse(BigFloat,x)

### Big

apprules(mx::Mxpr{:Big}) = do_Big(mx,mx[1])
do_Big(mx,x) = mx
do_Big(mx,x::T) where {T<:Number} = big(x)

### Minus

apprules(mx::Mxpr{:Minus}) = is_Number(mx[1]) ? -mx[1] : mmul(-1, mx[1])


text:
using ForwardDiff

abstract ExtendedKalmanFilter <: KalmanFilter

type NonlinearModel <: Model
    f::Function
    j::Function
    g::Function
    q::Matrix
    function NonlinearModel(f::Function,g::Function,q::Matrix)
        j = forwarddiff_jacobian(f,Float64,fadtype=:typed)
        new(f,j,g,q)
    end
end

function ap(f::NonlinearModel,x::State)
    x1 = f.f(x.x)
    F = f.j(x.x)
    G = f.g(x.x)
    p1 = F*x.p*F' + G*f.q*G'
    State(x1,p1)
end

type NonlinearObservationModel <: ObservationModel
    h::Function
    j::Function
    r::Matrix
    function NonlinearObservationModel(h::Function,r::Matrix)
        j = forwarddiff_jacobian(h,Float64,fadtype=:typed)
        new(h,j,r)
    end
end

type BasicExtendedKalmanFilter <: ExtendedKalmanFilter
    x::State
    f::NonlinearModel
    z::NonlinearObservationModel
end

function covs(kf::BasicExtendedKalmanFilter,y::Observation)
    res = y.y - kf.z.h(kf.x.x)
    H = kf.z.j(kf.x.x)'
    ph = kf.x.p * H'
    s = H * kf.x.p * H' + kf.z.r
    (res,ph,s)
end

ext: jl


text:
module ConformingFESpaces

using Gridap
using Gridap.Helpers
using Gridap.CellValuesGallery
using Gridap.CachedArrays
using Base: @propagate_inbounds
using LinearAlgebra: det
using Base: transpose

using Gridap.BoundaryGrids: _setup_tags

export ConformingFESpace
export H1ConformingFESpace
export DivConformingFESpace
export CurlConformingFESpace
import Gridap: num_free_dofs
import Gridap: num_diri_dofs
import Gridap: diri_tags
import Gridap: apply_constraints
import Gridap: apply_constraints_rows
import Gridap: apply_constraints_cols
import Gridap: celldofids
import Gridap: interpolate_values
import Gridap: interpolate_diri_values
import Gridap: CellField
import Gridap: CellBasis
import Gridap: Triangulation
import Base: size
import Base: getindex

"""
Conforming FE Space, where only one RefFE is possible in the whole mesh
"""
struct ConformingFESpace{D,Z,T} <: FESpace{D,Z,T}
  dim_to_nface_eqclass::Vector{<:IndexCellArray{Int}}
  cell_eqclass::IndexCellArray{Int}
  num_free_dofs::Int
  num_diri_dofs::Int
  diri_tags::Vector{Int}
  reffe::RefFE{Z,T}
  triangulation::Triangulation{D,Z}
  gridgraph::GridGraph
  facelabels::FaceLabels
  cellbasis::CellBasis{Z,T}
end

function ConformingFESpace(
  reffe::RefFE{D,T},
  trian::Triangulation{D,Z},
  graph::GridGraph,
  labels::FaceLabels,
  diri_tags::Vector{Int}) where {D,Z,T}
  args = _setup_conforming_fe_fields(reffe,trian,graph,labels,diri_tags,D)
  ConformingFESpace{D,Z,T}(args...)
end

function ConformingFESpace(
  reffe::RefFE{D,T},
  trian::Triangulation{D,Z},
  graph::GridGraph,
  labels::FaceLabels) where {D,Z,T}
  return ConformingFESpace(reffe, trian, graph, labels, ())
end

num_free_dofs(this::ConformingFESpace) = this.num_free_dofs

num_diri_dofs(this::ConformingFESpace) = this.num_diri_dofs

diri_tags(f::ConformingFESpace) = f.diri_tags

function apply_constraints(
  this::ConformingFESpace, cellvec::CellVector, cellids::CellNumber)
  cellvec
end

function apply_constraints_rows(
  this::ConformingFESpace, cellmat::CellMatrix, cellids::CellNumber)
  cellmat
end

function apply_constraints_cols(
  this::ConformingFESpace, cellmat::CellMatrix, cellids::CellNumber)
  cellmat
end

function celldofids(this::ConformingFESpace)
  this.cell_eqclass
end

function interpolate_values(this::ConformingFESpace,f::Function)
  _interpolate_values(this,f)
end

function interpolate_values(this::ConformingFESpace{D,Z,T},val::T) where {D,Z,T}
  fun(x) = val
  interpolate_values(this,fun)
end

function interpolate_diri_values(this::ConformingFESpace, funs::Vector{<:Function})
  _interpolate_diri_values(this,funs)
end

function interpolate_diri_values(this::ConformingFESpace{D,Z,T}, vals::Vector{T}) where {D,Z,T}
  _interpolate_diri_values(this,vals)
end

function CellField(
  fespace::ConformingFESpace{D,Z,T},
  free_dofs::AbstractVector{E},
  diri_dofs::AbstractVector{E}) where {D,Z,T,E}

  _CellField(fespace, free_dofs, diri_dofs,T,E)
end

CellBasis(this::ConformingFESpace) = this.cellbasis

Triangulation(this::ConformingFESpace) = this.triangulation

function H1ConformingFESpace(
  ::Type{T}, model::DiscreteModel{D}, order::Integer, diri_tags) where {D,T}

  labels = FaceLabels(model)
  H1ConformingFESpace(T,model,labels,order,diri_tags)
end

function H1ConformingFESpace(
  ::Type{T},
  model::DiscreteModel{D},
  labels::FaceLabels,
  order::Integer,
  diri_tags) where {D,T}

  grid = Grid(model,D)
  trian = Triangulation(grid)
  graph = GridGraph(model)
  orders = fill(order,D)
  polytope = _polytope(celltypes(grid))
  fe = LagrangianRefFE(T,polytope, orders)
  _diri_tags = _setup_tags(labels,diri_tags)
  ConformingFESpace(fe,trian,graph,labels,_diri_tags)
end

# @santiagobadia : Create a HDiv conforming constructor

# Helpers

function _CellField(
  fespace, free_dofs, diri_dofs, ::Type{T}, ::Type{E}) where {T,E}

  @assert E == eltype(T)
  @assert num_free_dofs(fespace) == length(free_dofs)
  @assert num_diri_dofs(fespace) == length(diri_dofs)
  reffe = fespace.reffe
  celldofs = celldofids(fespace)
  shb = CellBasis(fespace)
  cdofs = CellVectorFromLocalToGlobalPosAndNeg(
  celldofs, free_dofs, diri_dofs)
  lincomb(shb,cdofs)
  # @santiagobadia: For RT methods, we must add here a local_to_global_dofs
  # or global_to_local_dofs

end

function _setup_conforming_fe_fields(reffe,trian,graph,labels,diri_tags,D)
  dim_to_nface_eqclass, nfree, ndiri  = _generate_dim_to_nface_to_dofs(
  reffe, graph, labels, diri_tags)
  cellvefs_dim = [connections(graph,D,i) for i in 0:D]
  offset = length.(dim_to_nface_eqclass)
  for i in 2:length(offset)
    offset[i] += offset[i-1]
  end
  offset = tuple(offset[1:(end-1)]...)
  cellvefs = local_append(offset, cellvefs_dim...)
  dofs_all = append(dim_to_nface_eqclass...)
  cell_eqclass = CellEqClass(cellvefs, dofs_all, reffe)
  shb = ConstantCellValue(shfbasis(reffe), ncells(trian))
  phi = CellGeomap(trian)
  basis = attachgeomap(shb,phi)
  return dim_to_nface_eqclass, cell_eqclass, nfree, ndiri, diri_tags,
  reffe, trian, graph, labels, basis
end

function _generate_dim_to_nface_to_dofs(
  reffe::RefFE{D},
  graph::GridGraph,
  labels::FaceLabels,
  diri_tags::Vector{Int}) where D

  i_free_dof = 1
  i_diri_dof = -1

  dim_to_nface_to_dofs = IndexCellVector{Int}[]
  tag_to_labels = labels.tag_to_labels

  for d in 0:D

    nface_to_label = labels_on_dim(labels,d)
    cell_to_nfaces = connections(graph,D,d)
    nface_to_cells = connections(graph,d,D)
    icell = 1
    nface_to_cellowner = get_local_item(nface_to_cells,icell)
    nface_to_lnface = find_local_index(nface_to_cellowner, cell_to_nfaces)

    lnface_to_ldofs = nfacedofs(reffe,d)
    lnface_to_nldofs = length.(lnface_to_ldofs)

    num_nfaces = length(nface_to_cells)

    nface_to_dofs_ptrs = zeros(Int, num_nfaces+1)
    nface_to_dofs_data = Int[]

    if any(lnface_to_nldofs != 0)

      i_free_dof, i_diri_dof = _generate_nface_to_dofs!(
      nface_to_dofs_data,
      nface_to_dofs_ptrs,
      nface_to_lnface,
      lnface_to_nldofs,
      nface_to_label,
      diri_tags,
      tag_to_labels,
      i_free_dof,
      i_diri_dof)

    end

    length_to_ptrs!(nface_to_dofs_ptrs)

    nface_to_dofs = CellVectorFromDataAndPtrs(
    nface_to_dofs_data, nface_to_dofs_ptrs)

    push!(dim_to_nface_to_dofs, nface_to_dofs)

  end

  return (dim_to_nface_to_dofs, i_free_dof-1, -i_diri_dof-1)

end

function _generate_nface_to_dofs!(
  nface_to_dofs_data,
  nface_to_dofs_ptrs,
  nface_to_lnface,
  lnface_to_nldofs,
  nface_to_label,
  diri_tags,
  tag_to_labels,
  i_free_dof,
  i_diri_dof)

  for (nface, lnface) in enumerate(nface_to_lnface)

    nldofs = lnface_to_nldofs[lnface]
    label = nface_to_label[nface]
    isdiri = _is_diri(label,diri_tags,tag_to_labels)

    if isdiri
      for i in 1:nldofs
        push!(nface_to_dofs_data,i_diri_dof)
        i_diri_dof += -1
      end
    else
      for i in 1:nldofs
        push!(nface_to_dofs_data,i_free_dof)
        i_free_dof += 1
      end
    end

    nface_to_dofs_ptrs[nface+1] = nldofs
  end

  (i_free_dof, i_diri_dof)
end

_polytope(celltypes) = @notimplemented

function _polytope(celltypes::ConstantCellValue)
  code = celltypes.value
  Polytope(code)
end

function _interpolate_values(fesp::ConformingFESpace{D,Z,T},fun::Function) where {D,Z,T}
  reffe = fesp.reffe
  dofb = dofbasis(reffe)
  trian = fesp.triangulation
  phi = CellGeomap(trian)
  uphys = fun ∘ phi
  celldofs = fesp.cell_eqclass
  nfdofs = fesp.dim_to_nface_eqclass
  E = dof_type(T)
  free_dofs = zeros(E, num_free_dofs(fesp))
  diri_dofs = zeros(E, num_diri_dofs(fesp))
  aux = zeros(E, length(dofb))
  _interpolate_values_kernel!(
  free_dofs,diri_dofs,uphys,celldofs,dofb,aux)
  return free_dofs, diri_dofs
end

function _interpolate_values_kernel!(
  free_dofs,diri_dofs,uphys,celldofs,dofb,aux)

  for (imap,l2g) in zip(uphys,celldofs)
    evaluate!(dofb,imap,aux)
    # @santiagobadia : Here we should add a method for RT that multiplies by -1
    # if the face has this cell as the second one in the grid graph
    # local_to_global_dofs, global_to_local_dofs CellArray
    for (i,gdof) in enumerate(l2g)
      if (gdof > 0)
        free_dofs[gdof] = aux[i]
      else
        diri_dofs[-gdof] = aux[i]
      end
    end
  end
end

function _interpolate_diri_values(fesp::ConformingFESpace{D,Z,T},funs) where {D,Z,T}

  labels = fesp.facelabels
  dim_to_nface_to_label = labels.dim_to_nface_to_label
  dim_to_nface_to_dofs = fesp.dim_to_nface_eqclass
  diri_tags = fesp.diri_tags
  @assert length(diri_tags) == length(funs)
  E = eltype(T)
  diri_dof_to_val = zeros(E, num_diri_dofs(fesp))
  tag_to_labels = labels.tag_to_labels

  for (ifunc,f) in enumerate(funs)

    tag_f = diri_tags[ifunc]
    _ , diri_dof_to_fval = interpolate_values(fesp,f)

    if length(funs) == 1
      return diri_dof_to_fval
    end

    for idim in 0:D
      nface_to_label = dim_to_nface_to_label[idim+1]
      nface_to_dofs  = dim_to_nface_to_dofs[idim+1]
      _interpolate_diri_values_kernel!(
      diri_dof_to_val,
      diri_dof_to_fval,
      nface_to_label,
      nface_to_dofs,
      tag_f,
      tag_to_labels)
    end
  end

  return diri_dof_to_val

end

function _interpolate_diri_values_kernel!(
  diri_dof_to_val,
  diri_dof_to_fval,
  nface_to_label,
  nface_to_dofs,
  tag_f,
  tag_to_labels)

  for (nface, label) in enumerate(nface_to_label)
    dofs = nface_to_dofs[nface]
    if _is_diri(label,tag_f,tag_to_labels)
      for dof in dofs
        i = -dof
        diri_dof_to_val[i] = diri_dof_to_fval[i]
      end
    end
  end
end

@inline function _is_diri(label,diritags,tag_to_labels)
  for tag in diritags
    for dirilabel in tag_to_labels[tag]
      if (label == dirilabel)
        return true
      end
    end
  end
  return false
end

"""
Type encoding the cellwise local to global dof map.
For the moment, for the case of a single RefFE and for oriented nfaces.
"""
struct CellEqClass{T,A,B} <: IndexCellValue{CachedVector{T,Vector{T}},1}
  cell_to_nfaces::A
  nface_to_dofs::B
  lnface_to_ldofs::Vector{Vector{Int}}
  nldofs::Int
  cv::CachedVector{T,Vector{T}}
end

function CellEqClass(
  cell_to_nfaces::IndexCellVector{<:Integer},
  nface_to_dofs::IndexCellVector{T},
  reffe::RefFE) where T<:Integer

  lnface_to_ldofs = nfacedofs(reffe)
  A = typeof(cell_to_nfaces)
  B = typeof(nface_to_dofs)
  nldofs = length(dofbasis(reffe))
  v = zeros(T,nldofs)
  cv = CachedVector(v)
  CellEqClass{T,A,B}(cell_to_nfaces,nface_to_dofs,lnface_to_ldofs,nldofs,cv)
end

size(self::CellEqClass) = (length(self.cell_to_nfaces),)

@propagate_inbounds function getindex(self::CellEqClass,cell::Int)
  nfaces = self.cell_to_nfaces[cell]
  for (lnface,nface) in enumerate(nfaces)
    dofs = self.nface_to_dofs[nface]
    ldofs = self.lnface_to_ldofs[lnface]
    for i in 1:length(dofs)
      dof = dofs[i]
      ldof = ldofs[i]
      self.cv[ldof] = dof
    end
  end
  self.cv
end

function DivConformingFESpace(
  reffe::RefFE{D,T},
  trian::Triangulation{D,Z},
  graph::GridGraph,
  labels::FaceLabels,
  diri_tags::Vector{Int}) where {D,Z,T}

  dim_to_nface_eqclass,
  cell_eqclass,
  nfree,
  ndiri, diri_tags,
  reffe,
  trian,
  graph,
  labels,
  basis = _setup_conforming_fe_fields(reffe,trian,graph,labels,diri_tags,D)

  phi = CellGeomap(trian)
  jac = gradient(phi)
  detjac = det(jac)

  piola_map = (1.0/detjac)*jac

  # physbasis = _attachpiola(basis,piola_map)
  physbasis = basis

  ConformingFESpace{D,Z,T}( dim_to_nface_eqclass,
  cell_eqclass,
  nfree,
  ndiri, diri_tags,
  reffe,
  trian,
  graph,
  labels,
  physbasis )
end

function CurlConformingFESpace(
  reffe::RefFE{D,T},
  trian::Triangulation{D,Z},
  graph::GridGraph,
  labels::FaceLabels,
  diri_tags::Vector{Int}) where {D,Z,T}

  dim_to_nface_eqclass,
  cell_eqclass,
  nfree,
  ndiri, diri_tags,
  reffe,
  trian,
  graph,
  labels,
  basis = _setup_conforming_fe_fields(reffe,trian,graph,labels,diri_tags,D)

  phi = CellGeomap(trian)
  jac = gradient(phi)
  jact = transpose(jac)

  piola_map = inv(jact)

  # physbasis = _attachpiola(basis,piola_map)
  physbasis = piola_map

  ConformingFESpace{D,Z,T}( dim_to_nface_eqclass,
  cell_eqclass,
  nfree,
  ndiri, diri_tags,
  reffe,
  trian,
  graph,
  labels,
  physbasis )
end

# @santiagobadia : It is not ready because I think that the product of a
# integration point values times a
function _attachpiola(a::CellBasis{D},piola::CellFieldLike) where D
  Gridap.CellFieldsOperations._merge_val_and_grad(piola*a.val,piola*a.grad)
  # Gridap.CellFieldsOperations._merge_val_and_grad(a,gradient(a))
end

end # module


text:
using LTC
using GalacticOptim
using Plots
using Functors

function generate_data()
    in_features = 2
    out_features = 1
    N = 48
    data_x = [sin.(range(0,stop=3π,length=N)), cos.(range(0,stop=3π,length=N))]
    data_x = [reshape([Float32(data_x[1][i]),Float32(data_x[2][i])],2,1) for i in 1:N]# |> f32
    data_y = [reshape([Float32(y)],1) for y in sin.(range(0,stop=6π,length=N))]# |> f32

    # data_x = [repeat(x,1,20) for x in data_x]
    # data_y = [repeat(x,1,20) for x in data_y]

    data_x, data_y
end

function data(iter; data_x=nothing, data_y=nothing, short=false, noisy=false)
    #noisy_data = Vector{Tuple{Vector{Matrix{Float32}}, Vector{Vector{Float32}}}}([])
    if data_y === nothing
      data_x, data_y = generate_data()
    end
    noisy_data = []
    for i in 1:iter
        x = data_x
        y = data_y
        if short isa Array
          x = x[short[1]:short[2]]
          y = y[short[1]:short[2]]
        end
        push!(noisy_data, (x , noisy ? add_gauss.(y,0.02) : y))
    end
    noisy_data
end

data_x,data_y = generate_data()

wiring = Wiring(2,1)

display(Plots.heatmap(wiring.sens_mask))
display(Plots.heatmap(wiring.sens_pol))
display(Plots.heatmap(wiring.syn_mask))
display(Plots.heatmap(wiring.syn_pol))

ltc = NCP(wiring)
lower, upper = get_bounds(ltc)

#ltcp, ltcre = Functors.functor(ltc)
ltcp,ltcre = Flux.destructure(ltc)
ltcre(ltcp).([rand(Float32,2,1) for i in 1:50])

function loss_galactic(θ,p,x,y)
  m = ltcre(θ)
  ŷ = m.(x)
  sum(sum([(ŷ[i] .- y[i]) .^ 2 for i in 1:length(y)]))/length(y), ŷ
end
@time loss_galactic(ltcp,[],data_x,data_y,ltcre)

function cbs(θ,l,pred;doplot=true) #callback function to observe training
  display(l)
  if doplot
    #plot([x[1] for x in pred])
    #fig = scatter([(isinf(x[end,1]) || isnan(x[end,1])) ? missing : x[end,1] for x in pred])
    #scatter!(fig, [x[1,1] for x in y])
    #display(fig)
  end
  return false
end
cbs(ltcp,loss_galactic(ltcp,[],data_x,data_y,ltcre)...)

opt = GalacticOptim.Flux.Optimiser(GalacticOptim.Flux.ExpDecay(0.001, 0.1, 500, 1e-4), ADAM(0.001))

optfunc = GalacticOptim.OptimizationFunction((x, p, tx,ty) -> loss_galactic(x,p,tx,ty), GalacticOptim.AutoZygote())
optprob = GalacticOptim.OptimizationProblem(optfunc, ltcp; lb = lower, ub = upper)

res = GalacticOptim.solve(optprob, ADAM(0.1), data(2), cb = cbs, maxiters = 100)
res = GalacticOptim.solve(optprob, Fminbox(BFGS(initial_stepnorm = 0.01)), cb = cbs, allow_f_increases=false)
res = GalacticOptim.solve(optprob, Opt(:LN_BOBYQA, length(ltcp)), maxiters=100, cb = cbs)

sol = GalacticOptim.solve(optprob, NelderMead(), data(3), maxiters=3000, cb = cbs)


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#   Tests of highlighters.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

@testset "Highlighters" begin
    expected = """

    Col. 1
    Col. 2
    Col. 3
    Col. 4

    1
    false
    1.0
    1

    2
    true
    2.0
    2

    3
    false
    3.0
    3

    4
    true
    4.0
    4

    5
    false
    5.0
    5

    6
    true
    6.0
    6

"""

    result = pretty_table(String, data;
                          backend = :html,
                          highlighters = HTMLHighlighter((data, i, j)->data[i,j] > 4,
                                                         HTMLDecoration(color = "red")),
                          standalone = true)
end

@testset "Pre-defined highlighters" begin
    # hl_cell
    # ==========================================================================

    expected = """

    Col. 1
    Col. 2
    Col. 3
    Col. 4

    1
    false
    1.0
    1

    2
    true
    2.0
    2

    3
    false
    3.0
    3

    4
    true
    4.0
    4

    5
    false
    5.0
    5

    6
    true
    6.0
    6

"""

    result = pretty_table(String, data;
                          backend = :html,
                          highlighters = (hl_cell(2, 2, HTMLDecoration(background = "blue")),
                                          hl_cell(3, 3, HTMLDecoration(background = "blue"))),
                          standalone = false)

    @test result == expected

    result = pretty_table(String, data;
                          backend = :html,
                          highlighters = hl_cell([(2, 2), (3, 3)], HTMLDecoration(background = "blue")),
                          standalone = false)

    # hl_col
    # ==========================================================================

    expected = """

    Col. 1
    Col. 2
    Col. 3
    Col. 4

    1
    false
    1.0
    1

    2
    true
    2.0
    2

    3
    false
    3.0
    3

    4
    true
    4.0
    4

    5
    false
    5.0
    5

    6
    true
    6.0
    6

"""

    result = pretty_table(String, data;
                          backend = :html,
                          highlighters = (hl_col(1, HTMLDecoration(font_weight = "bold")),
                                          hl_col(4, HTMLDecoration(font_weight = "bold"))),
                          standalone = false)

    @test result == expected

    result = pretty_table(String, data;
                          backend = :html,
                          highlighters = hl_col([1, 4], HTMLDecoration(font_weight = "bold")),
                          standalone = false)

    @test result == expected

    # hl_row
    # ==========================================================================

    expected = """

    Col. 1
    Col. 2
    Col. 3
    Col. 4

    1
    false
    1.0
    1

    2
    true
    2.0
    2

    3
    false
    3.0
    3

    4
    true
    4.0
    4

    5
    false
    5.0
    5

    6
    true
    6.0
    6

"""

    result = pretty_table(String, data;
                          backend = :html,
                          highlighters = (hl_row(1, HTMLDecoration(font_family = "monospace")),
                                          hl_row(3, HTMLDecoration(font_family = "monospace")),
                                          hl_row(5, HTMLDecoration(font_family = "monospace"))),
                          standalone = false)

    @test result == expected

    result = pretty_table(String, data;
                          backend = :html,
                          highlighters = hl_row([1, 3, 5], HTMLDecoration(font_family = "monospace")),
                          standalone = false)

    @test result == expected

    # hl_lt
    # ==========================================================================

    expected = """

    Col. 1
    Col. 2
    Col. 3
    Col. 4

    1
    false
    1.0
    1

    2
    true
    2.0
    2

    3
    false
    3.0
    3

    4
    true
    4.0
    4

    5
    false
    5.0
    5

    6
    true
    6.0
    6

"""

    result = pretty_table(String, data;
                          backend = :html,
                          highlighters = hl_lt(3, HTMLDecoration(color = "white",
                                                                 background = "black")),
                          standalone = false)

    @test result == expected

    # hl_leq
    # ==========================================================================

    expected = """

    Col. 1
    Col. 2
    Col. 3
    Col. 4

    1
    false
    1.0
    1

    2
    true
    2.0
    2

    3
    false
    3.0
    3

    4
    true
    4.0
    4

    5
    false
    5.0
    5

    6
    true
    6.0
    6

"""

    result = pretty_table(String, data;
                          backend = :html,
                          highlighters = hl_leq(3, HTMLDecoration(color = "white",
                                                                  background = "black")),
                          standalone = false)

    @test result == expected

    # hl_gt
    # ==========================================================================

    expected = """

    Col. 1
    Col. 2
    Col. 3
    Col. 4

    1
    false
    1.0
    1

    2
    true
    2.0
    2

    3
    false
    3.0
    3

    4
    true
    4.0
    4

    5
    false
    5.0
    5

    6
    true
    6.0
    6

"""

    result = pretty_table(String, data;
                          backend = :html,
                          highlighters = hl_gt(3, HTMLDecoration(color = "white",
                                                                 background = "black")),
                          standalone = false)

    @test result == expected

    # hl_geq
    # ==========================================================================

    expected = """

    Col. 1
    Col. 2
    Col. 3
    Col. 4

    1
    false
    1.0
    1

    2
    true
    2.0
    2

    3
    false
    3.0
    3

    4
    true
    4.0
    4

    5
    false
    5.0
    5

    6
    true
    6.0
    6

"""

    result = pretty_table(String, data;
                          backend = :html,
                          highlighters = hl_geq(3, HTMLDecoration(color = "white",
                                                                  background = "black")),
                          standalone = false)

    @test result == expected

    # hl_value
    # ==========================================================================

    expected = """

    Col. 1
    Col. 2
    Col. 3
    Col. 4

    1
    false
    1.0
    1

    2
    true
    2.0
    2

    3
    false
    3.0
    3

    4
    true
    4.0
    4

    5
    false
    5.0
    5

    6
    true
    6.0
    6

"""

    result = pretty_table(String, data;
                          backend = :html,
                          highlighters = hl_value(3, HTMLDecoration(text_decoration = "line-through")),
                          standalone = false)

    @test result == expected
end


text:
using Documenter, DWDataReader

makedocs(;
    modules = [DWDataReader],
    format = Documenter.HTML(prettyurls = false),
    pages = [
        "Home" => "index.md",
    ],
    repo = "https://github.com/fleimgruber/DWDataReader.jl/blob/{commit}{path}#{line}",
    sitename = "DWDataReader.jl",
    authors = "Fabian Leimgruber",
)

deploydocs(; repo = "github.com/fleimgruber/DWDataReader.jl", devbranch = "main")


text:
module MushroomBody

export
run_model,
gif_model,
response_before_after_learning,
response_before_after_learning2,
plotmeanlayer,
meanlayer,
plot_before_after,
doboth

using SparseArrays, Infiltrator, ColorCoding, Random, StaticArrays, Plots
import Base: size, sum, show, getindex, setindex!, getproperty, setproperty!, zero, length, one, fill!
using Random, BenchmarkTools

include("Types.jl")
include("UpdateWeights.jl")
include("UpdateActivation.jl")
include("initialisers.jl")
include("Neurotransmitters.jl")
include("helpers.jl")
include("input.jl")
include("Parameters.jl")
include("Tests2.jl")
include("plotters.jl")
include("SaveData.jl")
include("ImageInput.jl")

"""
Run All the Steps
"""
function run_all_steps(nn, numsteps, m, p, sensory, da; update=true, savevars=nothing, showplot=false, reportvar=nothing, normweights=false)

	returnvariable = !isnothing(reportvar) && initialise_return_variable(numsteps,m,reportvar)

	for t = 1:numsteps

		da = run_step(t, m, p, sensory, da, nn, update=update, normweights=normweights)

		showplot && shownetwork(init_plot(), t, nn, m)

		!isnothing(savevars) && save_variables(m,savevars)

		if !isnothing(reportvar)
			returnvariable[t] = return_variable(m, reportvar)
		end
	end
	return returnvariable
end
function run_all_steps(nn, numsteps, m, p, sensory, da, gf::GifPlot; update=true, normweights=false)
	for t = 1:numsteps
		da = run_step(t, m, p, sensory, da, nn, update=update, normweights=normweights)
		gf = gifnetwork(gf, init_plot(), t, nn, m)
	end
	return gf
end

"""
run single step
"""
function run_step(t, m, p, sensory, da, nn; update=true, normweights=false)
	da = inputandreward!(t, m.input.layers[1], sensory, p.τ[1], p.da_on, da)
	for layer = 1:length(nn)
		run_layer(t, layer, nn, m, p, da, update=update)
	end
	normweights && normalise_layer!(m, p, l=(length(nn)-1)) # normalise the last layer.
	return da
end
"""
run Single layer
"""
function run_layer(t, layer, nn, m, p, da; update=true)
	update_activation!(t, layer, nn, m, p)
	if layer !== length(nn)
		update_pre_layers(t,layer,m,p, da, update=update)
	end
end
"""
weight-changing and output functions for all but last layers.
"""
function update_pre_layers(t, layer, m, p, da; update=true)
	update_ACh!(t, layer, m, p)
	calc_input!(layer, m, p)
	(layer == 2 && update) && update_weights!(t, layer, m, p, da)
end

"""
reset the dynamic aspects of the model, keeping the structur (weights and synapses) intact.
"""
function reset(nn, p, weights, synapses)

	return MatrixTypes(initialise_matrices(nn, p, weights, synapses)...)
end

function reset!(m, p)

	fill!(m.activation, p.vr)
	fill!(m.rec, 0.0)
	fill!(m.spiked, false)
	fill!(m.spt, Int(-1))
	fill!(m.I, 0.0)
	fill!(m.ACh, 0.0)
	fill!(m.input, 0.0)
	fill!(m.γ, 0.0)

end

end


text:
# Based on HttpParser.jl/test/runtests.jl

module HttpParserTest

using HttpParser
using Test
using Compat
import Compat: String

headers = Dict()

function on_message_begin(parser)
    return 0
end

function on_url(parser, at, len)
    return 0
end

function on_status_complete(parser)
    return 0
end

function on_header_field(parser, at, len)
    header = unsafe_string(convert(Ptr{UInt8}, at), Int(len))
    # set the current header
    headers["current_header"] = header
    return 0
end

function on_header_value(parser, at, len)
    s = unsafe_string(convert(Ptr{UInt8}, at), Int(len))
    # once we know we have the header value, that will be the value for current header
    headers[headers["current_header"]] = s
    # reset current_header
    headers["current_header"] = ""
    return 0
end

function on_headers_complete(parser)
    p = unsafe_load(parser)
    # get first two bits of p.type_and_flags

    # The parser type are the bottom two bits
    # 0x03 = 00000011
    ptype = p.type_and_flags & 0x03
    # flags = p.type_and_flags >>> 3
    if ptype == 0
        method = http_method_str(convert(Int, p.method))
    end
    if ptype == 1
        headers["status_code"] = string(convert(Int, p.status_code))
    end
    headers["http_major"] = string(convert(Int, p.http_major))
    headers["http_minor"] = string(convert(Int, p.http_minor))
    headers["Keep-Alive"] = string(http_should_keep_alive(parser))
    return 0
end

function on_body(parser, at, len)
    return 0
end

function on_message_complete(parser)
    return 0
end

function on_chunk_header(parser)
    return 0
end

function on_chunk_complete(parser)
    return 0
end

c_message_begin_cb = cfunction(on_message_begin, HttpParser.HTTP_CB...)
c_url_cb = cfunction(on_url, HttpParser.HTTP_DATA_CB...)
c_status_complete_cb = cfunction(on_status_complete, HttpParser.HTTP_CB...)
c_header_field_cb = cfunction(on_header_field, HttpParser.HTTP_DATA_CB...)
c_header_value_cb = cfunction(on_header_value, HttpParser.HTTP_DATA_CB...)
c_headers_complete_cb = cfunction(on_headers_complete, HttpParser.HTTP_CB...)
c_body_cb = cfunction(on_body, HttpParser.HTTP_DATA_CB...)
c_message_complete_cb = cfunction(on_message_complete, HttpParser.HTTP_CB...)
c_body_cb = cfunction(on_body, HttpParser.HTTP_DATA_CB...)
c_message_complete_cb = cfunction(on_message_complete, HttpParser.HTTP_CB...)
c_chunk_header_cb = cfunction(on_chunk_header, HttpParser.HTTP_CB...)
c_chunk_complete_cb = cfunction(on_chunk_complete, HttpParser.HTTP_CB...)

function parse(bytes)
    # reset request
    # Moved this up for testing purposes
    global headers = Dict()
    parser = Parser()
    http_parser_init(parser, false)
    settings = ParserSettings(c_message_begin_cb, c_url_cb,
                              c_status_complete_cb, c_header_field_cb,
                              c_header_value_cb, c_headers_complete_cb,
                              c_body_cb, c_message_complete_cb,
                              c_chunk_header_cb, c_chunk_complete_cb)

    size = http_parser_execute(parser, settings, bytes)

    return headers
end

end


text:
# # [Sparse linear regression](@id sparse_linreg)
#
# Let's look at a least squares regression problem with L1 regularization:
# we will use the "diabetes dataset" (see [here](https://www4.stat.ncsu.edu/~boos/var.select/)), so let's start by loading the data.

using HTTP

splitlines(s) = split(s, "\n")
splitfields(s) = split(s, "\t")
parsefloat64(s) = parse(Float64, s)

function load_diabetes_dataset()
    res = HTTP.request("GET", "https://www4.stat.ncsu.edu/~boos/var.select/diabetes.tab.txt")
    lines = res.body |> String |> strip |> splitlines
    return hcat((line |> splitfields .|> parsefloat64 for line in lines[2:end])...)'
end

data = load_diabetes_dataset()

training_input = data[1:end-100, 1:end-1]
training_label = data[1:end-100, end]

test_input = data[end-99:end, 1:end-1]
test_label = data[end-99:end, end]

n_training, n_features = size(training_input)

# Now we can set up the optimization problem we want to solve: we will minimize the mean squared error
# for a linear model, appropriately scaled so that the features in the training data are normally distributed
# ("standardization", this is known to help the optimization process).
#
# After some simple manipulation, this standardized linear model can be implemented as follows:

using LinearAlgebra
using Statistics

input_loc = mean(training_input, dims=1) |> vec
input_scale = std(training_input, dims=1) |> vec

linear_model(wb, input) = input * wb[1:end-1] .+ wb[end]

function standardized_linear_model(wb, input)
    w_scaled = wb[1:end-1] ./ input_scale
    wb_scaled = vcat(w_scaled, wb[end] - dot(w_scaled, input_loc))
    return linear_model(wb_scaled, input)
end

# The loss term in the cost is then the following. Note that this is a regular Julia function:
# since the algorithm we will apply requires its gradient, automatic differentiation will
# do the work for us.

mean_squared_error(label, output) = mean((output .- label) .^ 2) / 2

training_loss(wb) = mean_squared_error(training_label, standardized_linear_model(wb, training_input))

# As regularization we will use the L1 norm, implemented in [ProximalOperators](https://github.com/JuliaFirstOrder/ProximalOperators.jl):

using ProximalOperators

reg = ProximalOperators.NormL1(1)

# We want to minimize the sum of `training_loss` and `reg`, and for this task we can use `FastForwardBackward`,
# which implements the fast proximal gradient method (also known as fast forward-backward splitting, or FISTA).
# Therefore we construct the algorithm, then apply it to our problem by providing a starting point,
# and the objective terms `f=training_loss` (smooth) and `g=reg` (non smooth).

using ProximalAlgorithms

ffb = ProximalAlgorithms.FastForwardBackward()
solution, iterations = ffb(x0=zeros(n_features + 1), f=training_loss, g=reg)

# We can now check how well the trained model performs on the test portion of our data.

test_output = standardized_linear_model(solution, test_input)
mean_squared_error(test_label, test_output)


text:
#!/usr/bin/env julia

using EAGOBranchBound

# write your own tests here
println("Testing 1D Interval Optimization...")
t = @elapsed include("D1_Interval_Test.jl")
println("done (took $t seconds).")

println("Testing Option Setting...")
t = @elapsed include("Options_Tests.jl")
println("done (took $t seconds).")

println("Testing Nodes...")
t = @elapsed include("Node_Tests.jl")
println("done (took $t seconds).")

println("Testing Bisection...")
t = @elapsed include("Bisect_Tests.jl")
println("done (took $t seconds).")

println("Testing Branch...")
t = @elapsed include("Branch_Tests.jl")
println("done (took $t seconds).")

println("Testing Access...")
t = @elapsed include("Access_Tests.jl")
println("done (took $t seconds).")

println("Testing Checks...")
t = @elapsed include("Checks_Tests.jl")
println("done (took $t seconds).")

println("Testing Display...")
t = @elapsed include("Display_Tests.jl")
println("done (took $t seconds).")


text:
@testset "slfm" begin
    rng = MersenneTwister(123)
    FDM = FiniteDifferences.central_fdm(5, 1)
    N = 10
    in_dim = 5
    out_dim = 4
    x1 = MOInput([rand(rng, in_dim) for _ in 1:N], out_dim)
    x2 = MOInput([rand(rng, in_dim) for _ in 1:N], out_dim)

    k = LatentFactorMOKernel(
        [Matern32Kernel(), SqExponentialKernel(), FBMKernel()],
        IndependentMOKernel(GaussianKernel()),
        rand(rng, out_dim, 3),
    )
    @test k isa LatentFactorMOKernel
    @test k isa MOKernel
    @test k isa Kernel
    @test k(x1[1], x2[1]) isa Real

    @test kernelmatrix(k, x1, x2) ≈ kernelmatrix(k, collect(x1), collect(x2))
    @test kernelmatrix(k, x1, x1) ≈ kernelmatrix(k, x1)

    @test string(k) == "Semi-parametric Latent Factor Multi-Output Kernel"
    @test repr("text/plain", k) == (
        "Semi-parametric Latent Factor Multi-Output Kernel\n\tgᵢ: " *
        "Matern 3/2 Kernel\n\t\tSquared Exponential Kernel\n" *
        "\t\tFractional Brownian Motion Kernel (h = 0.5)\n\teᵢ: " *
        "Independent Multi-Output Kernel\n\tSquared Exponential Kernel"
    )

    # AD test
    function test_slfm(A::AbstractMatrix, x1, x2)
        k = LatentFactorMOKernel(
            [Matern32Kernel(), SqExponentialKernel(), FBMKernel()],
            IndependentMOKernel(GaussianKernel()),
            A,
        )
        return k((x1, 1), (x2, 1))
    end

    a = rand()
    @test all(
        FiniteDifferences.j′vp(FDM, test_slfm, a, k.A, x1[1][1], x2[1][1]) .≈
        Zygote.pullback(test_slfm, k.A, x1[1][1], x2[1][1])[2](a),
    )
end


text:
using SMLMData
using Documenter

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

makedocs(;
    modules=[SMLMData],
    sitename="SMLMData.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://JuliaSMLM.github.io/SMLMData.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
        "Library.md"
    ],
)

deploydocs(;
    repo="github.com/JuliaSMLM/SMLMData.jl",
    devbranch = "main"
)


text:
# developed with Julia 1.4.2
#
# functions to evaluate a model based on performance metrics
# computed after simulation

function average_cost_per_site(result_model::String,
	result_dummy::String=joinpath(DIR, "metadata", "baseline", "dummy.jld2"),
	result_anticipative::String=joinpath(DIR, "metadata", "baseline", "anticipative.jld2"))

    model = load(result_model)
    dummy = load(result_dummy)
    anticipative = load(result_anticipative)
    scores = DataFrame(site=String[], period=String[], model=Float64[],
    	dummy=Float64[], anticipative=Float64[])

    sites = keys(model)
    ordered_sites = string.(sort!(parse.(Int64, string.(sites))))

    for site in ordered_sites

        simulations_model = model[site]
        simulations_dummy = dummy[site]
        simulations_anticipative = anticipative[site]

        for simulation in zip(simulations_model, simulations_dummy,
        	simulations_anticipative)

            model_cost = sum(simulation[1].result.cost)
            dummy_cost = sum(simulation[2].result.cost)
            anticipative_cost = sum(simulation[3].result.cost)
            push!(scores, [simulation[1].id.site_id, simulation[1].id.period_id,
                    model_cost, dummy_cost, anticipative_cost])

        end

    end

    summary = DataFrame(site=String[], model=Float64[], dummy=Float64[], anticipative=Float64[])

    for site in ordered_sites
        data = scores[scores.site .== site, :]
        model_mean_cost = mean(data[!, :model])
        dummy_mean_cost = mean(data[!, :dummy])
        anticipative_mean_cost = mean(data[!, :anticipative])
        push!(summary, [site, model_mean_cost, dummy_mean_cost, anticipative_mean_cost])
    end

    return summary

end

function evaluate_model(result_model::String,
	result_dummy::String=joinpath(DIR, "metadata", "baseline", "dummy.jld2"),
	result_anticipative::String=joinpath(DIR, "metadata", "baseline", "anticipative.jld2"))

    average_costs = average_cost_per_site(result_model, result_dummy, result_anticipative)

    DataFrame(site=average_costs[!, :site],
    	cost=average_costs[!, :model],
    	gain=average_costs[!, :dummy] - average_costs[!, :model],
    	score=(average_costs[!, :dummy] - average_costs[!, :model]) ./ (average_costs[!, :dummy]
    	- average_costs[!, :anticipative]))

end


text:
#! /usr/bin/env julia

println("Start test 'tests/train.jl'")

include("../src/DCellC.jl")
using DCellC

model = Multiscale3(GreyscaleImage, bn = true, seed = 42)
println("Test train.jl: Creating instance of FCRNA works")

lmgs  = [ synthesize(256, 256, (1500, 2000)) for i in 1:2 ]
println("Test train.jl: Creating synthesized images works")

train!(model, lmgs, patchsize = 128, epochs = 1, patchmode = 2, lr = 1e-3)
train!(model, lmgs, epochs = 1, lr = 1e-3)
println("Test train.jl: Training of model was sucessful")

println("Test train.jl completed")

max_stars_repo_name: tscode/DCellC

max_issues_repo_name: tscode/DCellC


text:
int main()
{
   f(10, 7);
   g(2, 3);
   h(3, 4, 5);
   return 0;
}

void f(int x, int y)
{
   int v, w;
   w = x + y;
   v = x + y;
   int z = v + 2 * w;
   printInt(z);
}

void g(int u, int v)
{
   int z = (v * u) + (u * v) - (3 * ((u * v) + (v * u)));
   printInt(z);
}

void h(int x, int y, int z)
{
   printInt(x * y + x * z);
}

max_stars_repo_name: lukaszcz/javalette

max_issues_repo_name: lukaszcz/javalette

max_forks_repo_name: lukaszcz/javalette


text:
#!/usr/bin/julia
using ExperimentalDesign, StatsModels, GLM, DataFrames, Distributions, Random, CSV, StatsBase, Plots, LinearAlgebra

function plot_df(results)
	mse(y1, y2) = mean((y1 - y2) .^ 2)

	sort!(results, [:expected_insns])

	transform!(results, AsTable(5:15) .=>
						  ByRow.([mean, sqrt∘var]) .=> [:mean, :stdev])

	parallel2 = filter(Row -> Row.parallel == 2, results)
	parallel4 = filter(Row -> Row.parallel == 4, results)
	parallel8 = filter(Row -> Row.parallel == 8, results)
	seq = filter(Row -> Row.parallel == 1, results)

	mse_seq       = mse(seq.mean, seq.predicted)
	mse_parallel2 = mse(parallel2.mean, parallel2.predicted)
	mse_parallel4 = mse(parallel4.mean, parallel4.predicted)
	mse_parallel8 = mse(parallel8.mean, parallel8.predicted)

	mse_seq_str   = string(round(mse_seq, digits=2))
	mse_par2_str  = string(round(mse_parallel2, digits=2))
	mse_par4_str  = string(round(mse_parallel4, digits=2))
	mse_par8_str  = string(round(mse_parallel8, digits=2))

	println("MSE Sequential = " * string(mse_seq))
	println("MSE Parallel 2 = " * string(mse_parallel2))
	println("MSE Parallel 4 = " * string(mse_parallel4))
	println("MSE Parallel 8 = " * string(mse_parallel8))

	pltobj = plot(seq.expected_insns, seq.mean, yerror=1.96*seq.stdev/sqrt(15), label="Sequential", xaxis="Number of Instructions", yaxis="Time (s)")
	plot!(pltobj, parallel8.expected_insns, parallel8.mean, yerror=1.96*seq.stdev/sqrt(15), label="Parallel (8 threads)")
	plot!(pltobj, seq.expected_insns, seq.predicted, label=("Sequential Pred. (mse=" * mse_seq_str * ")"))
	plot!(pltobj, parallel8.expected_insns, parallel8.predicted, label=("Parallel Pred. (mse=" * mse_par8_str * ")"))

	savefig("times-insns.svg")

	sort!(seq, [:functions, :mean])
	sort!(parallel8, [:functions, :mean])

	pltobj = plot(seq.functions, seq.mean, yerror=1.96*seq.stdev/sqrt(15), label="Sequential", xaxis="Number of Functions", yaxis="Time (s)")
	plot!(pltobj, parallel8.functions, parallel8.mean, yerror=1.96*seq.stdev/sqrt(15), label="Parallel (8 threads)")
	plot!(pltobj, seq.functions, seq.predicted, label=("Sequential Pred. (mse=" * mse_seq_str * ")"))
	plot!(pltobj, parallel8.functions, parallel8.predicted, label=("Parallel Pred. (mse=" * mse_par8_str * ")"))

	savefig("times-functions.svg")
end

# Rename labels to lowercase, as @formula seems to not like uppercase strings.
function lowercase_columns!(D)
	for name in names(D)
		h = name => lowercase(name)
		rename!(D, replace.(names(D), h))
	end
end

function kl_prediction(D, formula, experiments)
	selected_rows = kl_exchange(formula,
								D,
								seed_design_size = 2,
								experiments = experiments,
								design_k = experiments,
								candidates_l = 3000)
								#candidates_l = size(D, 1) - experiments)

	return selected_rows.indices[1]
end

# Open Data for analysis
results = CSV.read("output_corei7.csv", delim="; ", DataFrame)

# Rename labels to lowercase
lowercase_columns!(results)

# Filter 0-insns programs
filter!(:expected_insns => x -> x > 0, results)

# Stack data on exectution_i into multiple lines
stacked = stack(results, Not([:filename, :functions, :expected_insns, :parallel, :inlined_functions]))

# Set regression formula
regression_formula = @formula(0 ~ ((expected_insns + expected_insns ^ 2) +
                                            (functions + functions ^ 2) +
                                            (inlined_functions + inlined_functions ^ 2)) *
                                            (parallel + parallel ^ 2))

## O que fazer com isso ?
experiments = 30
train_rows = kl_prediction(select(stacked, Not([:filename, :value, :variable])),  regression_formula, experiments)
##println(stacked[kl])

# Set train set
#train_ratio = 0.2
#train_size = round(Int, 0.2 * nrow(results))
#train_rows = shuffle(1:nrow(results))[1:train_size]
train = stacked[train_rows, :]
println(train)

# Set test set
#test_rows = [x for x in 1:nrow(results) if !(x in train_rows)]
#test  = stacked[test_rows, :]

regression_formula = @formula(value ~ ((expected_insns + expected_insns ^ 2) +
                                            (functions + functions ^ 2) +
                                            (inlined_functions + inlined_functions ^ 2)) *
                                            (parallel + parallel ^ 2))

# Fit data into test set using formula
screening_fit = lm(regression_formula, train)

# Print coefficient and R2 error
println(screening_fit)
println("R² = " * string(adjr2(screening_fit)))

# Predict data based on fitted info
results[!, :predicted] = predict(screening_fit, results)

CSV.write("predicted.csv", results)
plot_df(results)

max_stars_repo_name: giulianobelinassi/LTO-Timing-Analysis

max_issues_repo_name: giulianobelinassi/LTO-Timing-Analysis

max_forks_repo_name: giulianobelinassi/LTO-Timing-Analysis


text:
series = "
73167176531330624919225119674426574742355349194934
96983520312774506326239578318016984801869478851843
85861560789112949495459501737958331952853208805511
12540698747158523863050715693290963295227443043557
66896648950445244523161731856403098711121722383113
62229893423380308135336276614282806444486645238749
30358907296290491560440772390713810515859307960866
70172427121883998797908792274921901699720888093776
65727333001053367881220235421809751254540594752243
52584907711670556013604839586446706324415722155397
53697817977846174064955149290862569321978468622482
83972241375657056057490261407972968652414535100474
82166370484403199890008895243450658541227588666881
16427171479924442928230863465674813919123162824586
17866458359124566529476545682848912883142607690042
24219022671055626321111109370544217506941658960408
07198403850962455444362981230987879927244284909188
84580156166097919133875499200524063689912560717606
05886116467109405077541002256983155200055935729725
71636269561882670428252483600823257530420752963450
"

series = filter(x -> !isspace(x), series)
max = 0
biglist=[]
for i = 1:length(series)-13
    rolling = series[i:i+13]
    products = 1
    for i in rolling
        products *= parse(Int,i)
    end
    if products > max
        global max = products
        global biglist = rolling
    end
end
println(max)
println("List = $biglist")


text:
using Test
import TreeLas.Dual: dual!
import Random: MersenneTwister

@testset "dual: test_dual3               " begin
    root = 1
    local parent = Int[1, 1, 1]
    z = Float64[-3, 1, 2]
    proc_order = Int[2, 3, 1]
    @test sum(z) ≈ 0
    alpha = copy(z)
    dual!(alpha, proc_order, parent, NaN)
    @test isnan(alpha[root])
    @test alpha[2:end] == [1, 2]
end

@testset "dual: cumsum                   " begin
    n = 10
    rseed() = MersenneTwister(42)
    pi = collect(0:(n-1))
    pi[1] = 1
    post_order = n:-1:1 |> collect
    z = randn(rseed(), n)
    z .-= sum(z)/n
    @test sum(z) <= eps() * n
    alpha = copy(z)
    dual!(alpha, post_order, pi)
    @test alpha ≈ reverse!(cumsum(reverse(z)))
end

max_stars_repo_name: EQt/treelas

max_issues_repo_name: EQt/treelas

max_forks_repo_name: EQt/treelas


text:
#========================================================================================#
"""
	QuBits

Module QuBits: My Very Own Quantum Computer

Author: Niall Palfreyman, 03/05/22
"""
module QuBits

# Externally available names:
export Tensor, ⊗, kron, ⊚, kpow, isclose, ishermitian, isunitary
export State, qubit, ampl, phase, prob, maxprob, nbits

# Imports:
import Base:kron

using LinearAlgebra

#========================================================================================#
# Tensor definitions:
#-----------------------------------------------------------------------------------------
"""
	Tensor

Tensor is a synonym for AbstractArray.
"""
Tensor = AbstractArray

#-----------------------------------------------------------------------------------------
"""
	Kronecker product ⊗

Define ⊗ as an infix operator for the Kronecker product of two Tensors
"""
⊗ = kron

#-----------------------------------------------------------------------------------------
# Tensor methods:

#-----------------------------------------------------------------------------------------
"""
	kpow( t::Tensor, n::Int)

Kronecker the Tensor t n times with itself.
"""
function kpow( t::Tensor, n::Int)
	if n==0
		# Base case:
		return 1.0
	end

	# Build Kronecker product:
	tn = t
	for _ in 1:(n-1)
		tn = tn ⊗ t
	end
	tn
end
# Use circled ring as operator for kpow:
⊚ = kpow

#-----------------------------------------------------------------------------------------
"""
	isclose( t::Tensor, u::Tensor, atol::Float64=1e-6)

Check whether all elements of the Tensor t are close to those of the Tensor u.
"""
function isclose( t, u, atol::Float64=1e-6)
	all(abs.(t.-u) .< atol)
end

#-----------------------------------------------------------------------------------------
"""
	ishermitian( t::Tensor)

Check whether the Tensor t is hermitian.
"""
function ishermitian( t::Tensor)
	sz = size(t)
	if length(sz) != 2 || sz[1] != sz[2]
		return false
	end

	isclose(t,t')
end

#-----------------------------------------------------------------------------------------
"""
	isunitary( t::Tensor)

Check whether the Tensor t is unitary.
"""
function isunitary( t::Tensor)
	sz = size(t)
	if length(sz) != 2 || sz[1] != sz[2]
		return false
	end

	isclose(t*t',Matrix(I,size(t)))
end

#========================================================================================#
# State definitions:
#-----------------------------------------------------------------------------------------
"""
	State

A State is a collection of one or more qubits: a wrapper for Vector{Complex}.
"""
struct State <: AbstractVector{Complex}
	amp::Vector{Complex}								# The State amplitudes

	function State(amp::Vector)
		if isclose(amp,0.0)
			error("State vectors must be non-zero")
		end
		new(amp/norm(amp))
	end
end

#-----------------------------------------------------------------------------------------
# Delegated methods:
Base.length(s::State) = length(s.amp)
Base.size(s::State) = size(s.amp)
Base.setindex(s::State,i) = setindex(s.amp,i)
Base.getindex(s::State,i) = getindex(s.amp,i)

#-----------------------------------------------------------------------------------------
"""
	kron( s1::State, s2::State)

Return Kronecker product of two States.
"""
function kron( s1::State, s2::State)
	State(kron(s1.amp,s2.amp))
end

#-----------------------------------------------------------------------------------------
"""
	ampl( s::State, bits::BitVector)

Return amplitude for qubit indexed by bits.
"""
function ampl( s::State, bits::BitVector)
	s.amp[1+bits2dec(bits)]					# Add 1 to index because bits counts from 0
end

ampl( s::State, bits...) = ampl(s,BitVector(bits))

#-----------------------------------------------------------------------------------------
"""
	phase( s::State, bits::BitVector)

Return phase of the amplitude of the qubit indexed by bits.
"""
function phase( s::State, bits::BitVector)
	angle(ampl(s,bits))
end

phase( s::State, bits...) = phase(s,BitVector(bits))

#-----------------------------------------------------------------------------------------
"""
	prob( s::State, bits::BitVector)

Return probability for qubit indexed by bits.
"""
function prob( s::State, bits::BitVector)
	amp = ampl(s,bits)						# Retrieve own amplitude
	real(amp'*amp)							# Calculate its squared norm
end

prob( s::State, bits...) = prob(s,BitVector(bits))

#-----------------------------------------------------------------------------------------
"""
	normalise( s::State)

Normalise the state s (but throw an error if it's a zero state).
"""
function normalise(s::State)
	nrm = norm(s)
	if isclose(nrm,0.0)
		error("Attempting to normaise zero-probability state.")
	end
	s /= nrm
end

#-----------------------------------------------------------------------------------------
"""
	qubit( alpha::Complex, beta::Complex)

Create a single qubit State from the given amplitudes.
"""
function qubit( alpha::Number, beta=nothing)
	if beta===nothing
		beta = sqrt(1.0 - alpha'*alpha)
	end

#=	if !isclose( alpha'*alpha+beta'*beta, 1.0)
		error("Qubit probabilities do not add up to 1.")
	end=#

	State([alpha,beta])
end

function qubit(; alpha=nothing, beta=nothing)
	if alpha === beta === nothing
		# Neither amplitude is given:
		error("alpha, beta or both are required.")
	end

	if alpha === nothing
		# alpha not given:
		alpha = sqrt(1.0 - beta'*beta)
	end

	qubit(alpha,beta)
end

#-----------------------------------------------------------------------------------------
"""
	nbits( stt::State)

Return number of qubits in this state.
"""
function nbits( s::State)
	Int(log2(length(s)))
end

#-----------------------------------------------------------------------------------------
"""
	maxprob( stt::State)

Return tuple (bitindex,prob) of highest probability qubit in this State.
"""
function maxprob( s::State)
	mxbindex, mxprob = BitVector([]), 0.0
	for bindex in bitprod(nbits(s))
		thisprob = prob(s,bindex)
		if thisprob > mxprob
			mxbindex, mxprob = bindex, thisprob
		end
	end

	(mxbindex, mxprob)
end

#-----------------------------------------------------------------------------------------
"""
	Base.show( s::State)

Display the given State.
"""
function Base.show( io::IO, s::State)
	print( io, "State{", nbits(s), "}[")
	for i in 1:length(s)-1
		print( round(s[i],digits=4), ", ")
	end
	print( io, round(s[end],digits=4), "]")
end

#========================================================================================#
# Helper methods:
#-----------------------------------------------------------------------------------------
"""
	bits2dec( bits::BitVector)

Compute decimal representation of a BitVector.
"""
function bits2dec( bits::BitVector)
    s = 0; v = 1
   for i in view(bits,length(bits):-1:1)
        s += v*i
        v <<= 1
    end
    s
end

#-----------------------------------------------------------------------------------------
"""
	dec2bits( dec::Int, nbits::Int)

Compute a hi2lo nbit binary representation of a decimal value.
"""
function dec2bits( dec::Int, nbits::Int)
    BitVector(reverse( digits(dec, base=2, pad=nbits)))
end

#-----------------------------------------------------------------------------------------
"""
	bitprod( nbits::Int)

Construct in numerical order a list of all binary numbers containing nbits.
"""
function bitprod( nbits::Int)
    [dec2bits(i,nbits) for i in 0:2^nbits-1]
end

#-----------------------------------------------------------------------------------------
function demo()
	p1 = qubit(alpha=rand())
	x1 = qubit(beta=rand())
	psi = p1 ⊗ x1

	println( isclose( psi'*psi, 1.0))
	println( isclose( p1'*p1*x1'*x1, 1.0))
	println( ampl(psi,1,0))
	println( prob(psi,1,0))
end

end		# ... of module QuBits