text:
"""
`module  JAC.MultiPhotonDeExcitation`
    ... a submodel of JAC that contains all methods for computing multi-photon excitation and decay rates.
"""
module MultiPhotonDeExcitation

    using Printf, ..Basics, ..ManyElectron, ..Nuclear, ..Radial, ..TableStrings

    """
    `struct  MultiPhotonDeExcitation.Settings`  ... defines a type for the settings in estimating multi-photon excitation and decay rates.

        + NoPhotons               ::Int64                        ... Number of photons in the multi-photon ionization
        + multipoles              ::Array{EmMultipole}           ... Specifies the multipoles of the radiation field that are to be included.
        + gauges                  ::Array{UseGauge}              ... Specifies the gauges to be included into the computations.
        + printBeforeComputation  ::Bool                         ... True, if all energies and lines are printed before their evaluation.
        + selectLines             ::Bool                         ... True, if lines are selected individually for the computations.
        + selectedLines           ::Array{Tuple{Int64,Int64},1}  ... List of lines, given by tupels (inital-level, final-level).

    """
    struct Settings
        NoPhotons                 ::Int64
        multipoles                ::Array{EmMultipole}
        gauges                    ::Array{UseGauge}
        printBeforeComputation    ::Bool
        selectLines               ::Bool
        selectedLines             ::Array{Tuple{Int64,Int64},1}
    end

    """
    `MultiPhotonDeExcitation.Settings()`  ... constructor for the default values of multi-photon excitation and decay rates.
    """
    function Settings()
        Settings(0, EmMultipole[], UseGauge[], false, false, Tuple{Int64,Int64}[])
    end

    # `Base.show(io::IO, settings::MultiPhotonDeExcitation.Settings)`
    #		... prepares a proper printout of the variable settings::MultiPhotonDeExcitation.Settings.
     function Base.show(io::IO, settings::MultiPhotonDeExcitation.Settings)
        println(io, "NoPhotons:                $(settings.NoPhotons)  ")
        println(io, "multipoles:               $(settings.multipoles)  ")
        println(io, "gauges:                   $(settings.gauges)  ")
        println(io, "printBeforeComputation:   $(settings.printBeforeComputation)  ")
        println(io, "selectLines:              $(settings.selectLines)  ")
        println(io, "selectedLines:            $(settings.selectedLines)  ")
    end

    """
    `struct  MultiPhotonDeExcitation.Channel`
        ... defines a type for a multi-photon channel to help characterize multi-photon (single-electron) ionization with well-defined
            multipolarities.

        + multipoles     ::Array{EmMultipole,1}   ... Multipoles of all N incoming/outgoing photons.
        + gauge          ::UseGauge               ... Gauge for dealing with the (coupled) radiation field.
        + amplitude      ::Complex{Float64}       ... multi-photon ionization amplitude associated with the given channel.
   """
    struct  Channel
        multipoles       ::Array{EmMultipole,1}
        gauge            ::UseGauge
        amplitude        ::Complex{Float64}
    end

    """
    `struct  MultiPhotonDeExcitation.Line`  ... defines a type for a multi-photon ionization line that may include the definition of channels.

        + initialLevel     ::Level                  ... initial-(state) level
        + finalLevel       ::Level                  ... final-(state) level
        + photonEnergy     ::Float64                ... Energy of the incoming photons; all photons are assumed to have equal energy.
        + totalRate        ::EmProperty             ... Cross section for this multi-photon ionization.
        + hasChannels      ::Bool                   ... Determines whether the individual (sub-) channels are defined in terms of their
                                                        multipolarities, etc., or not.
        + channels         ::Array{MultiPhotonDeExcitation.Channel,1}  ... List of MultiPhotonDeExcitation.Channels of this line.
    """
    struct  Line
        initialLevel       ::Level
        finalLevel         ::Level
        photonEnergy       ::Float64
        totalRate          ::EmProperty
        hasChannels        ::Bool
        channels           ::Array{MultiPhotonDeExcitation.Channel,1}
    end

    # `Base.show(io::IO, line::MultiPhotonDeExcitation.Line)`  ... prepares a proper printout of the variable line::MultiPhotonDeExcitation.Line.
    function Base.show(io::IO, line::MultiPhotonDeExcitation.Line)
        println(io, "initialLevel:      $(line.initialLevel)  ")
        println(io, "finalLevel:        $(line.finalLevel)  ")
        println(io, "photonEnergy:      $(line.photonEnergy)  ")
        println(io, "totalRate:         $(line.totalRate)  ")
        println(io, "hasChannels:       $(line.hasChannels)  ")
        println(io, "channels:          $(line.channels)  ")
    end

    """
    `MultiPhotonDeExcitation.computeAmplitudesProperties(line::MultiPhotonDeExcitation.Line, grid::Radial.Grid,
                                                             settings::MultiPhotonDeExcitation.Settings)`
        ... to compute all amplitudes and properties of the given line; a line::Einstein.Line is returned for which the amplitudes
            and properties are now evaluated.
    """
    function  computeAmplitudesProperties(line::MultiPhotonDeExcitation.Line, grid::Radial.Grid, settings::MultiPhotonDeExcitation.Settings)
        global JAC_counter
        newChannels = MultiPhotonDeExcitation.Channel[]
        for channel in line.channels
            # matrix    = MultiPhotonDeExcitation.computeMatrix(channel.multipole, channel.gauge, line.omega, line.finalLevel.basis,
            #                                                                                line.initialLevel.basis, grid, settings)
            # amplitude = line.finalLevel.mc * matrix * line.initialLevel.mc
            JAC_counter = JAC_counter + 1
            if   JAC_counter < 20   println("MultiPhotonDeExcitation.computeAmplitudesProperties-aa: warning ... amplitude always -1.")    end
            amplitude = 1.0
            push!( newChannels, MultiPhotonDeExcitation.Channel( channel.multipoles, channel.gauge, amplitude) )
        end
        # Calculate the photonrate and angular beta if requested
        JAC_counter = JAC_counter + 1
        if   JAC_counter < 20   println("MultiPhotonDeExcitation.computeAmplitudesProperties-ab: warning ... rate always -1.")    end
        totalRate = EmProperty(-1., -1.)
        line = MultiPhotonDeExcitation.Line( line.initialLevel, line.finalLevel, line.photonEnergy, totalRate, true, newChannels)
        return( line )
    end

    """
    `MultiPhotonDeExcitation.computeLines(finalMultiplet::Multiplet, initialMultiplet::Multiplet, grid::Radial.Grid,
                                          settings::MultiPhotonDeExcitation.Settings; output=true)`
        ... to compute the multiphoton transition amplitudes and all properties as requested by the given settings. A list of
            lines::Array{MultiPhotonDeExcitation.Lines} is returned.
    """
    function  computeLines(finalMultiplet::Multiplet, initialMultiplet::Multiplet, grid::Radial.Grid,
                           settings::MultiPhotonDeExcitation.Settings; output=true)
        println("")
        printstyled("MultiPhotonDeExcitation.computeLines(): The computation of multiphoton transition amplitudes starts now ... \n", color=:light_green)
        printstyled("----------------------------------------------------------------------------------------------------------- \n", color=:light_green)
        println("")
        #
        lines = MultiPhotonDeExcitation.determineLines(finalMultiplet, initialMultiplet, settings)
        # Display all selected lines before the computations start
        if  settings.printBeforeComputation    MultiPhotonDeExcitation.displayLines(lines)    end
        # Calculate all amplitudes and requested properties
        newLines = MultiPhotonDeExcitation.Line[]
        for  line in lines
            newLine = MultiPhotonDeExcitation.computeAmplitudesProperties(line, grid, settings)
            push!( newLines, newLine)
        end
        # Print all results to screen
        MultiPhotonDeExcitation.displayRates(stdout, lines)
        printSummary, iostream = Defaults.getDefaults("summary flag/stream")
        if  printSummary    MultiPhotonDeExcitation.displayRates(iostream, lines)   end
        #
        if    output    return( lines )
        else            return( nothing )
        end
    end

    """
    `MultiPhotonDeExcitation.determineChannels(finalLevel::Level, initialLevel::Level, settings::MultiPhotonDeExcitation.Settings)`
        ... to determine a list of MultiPhotonDeExcitation.Channel for a transitions from the initial to final level and by taking into account
            the particular settings of for this computation; an Array{MultiPhotonDeExcitation.Channel,1} is returned.
    """
    function determineChannels(finalLevel::Level, initialLevel::Level, settings::MultiPhotonDeExcitation.Settings)
        global JAC_counter
        channels = MultiPhotonDeExcitation.Channel[];
        symi = LevelSymmetry(initialLevel.J, initialLevel.parity);    symf = LevelSymmetry(finalLevel.J, finalLevel.parity)
        mpListList = Array{Basics.EmMultipole,1}[];    morePhotons = true;    NoPhotons = 1
        for  mp in settings.multipoles    push!( mpListList, [mp] )    end
        if  NoPhotons + 1 > settings.NoPhotons    error("stop a")      end
        while  morePhotons
            NoPhotons     = NoPhotons + 1
            newMpListList = Array{Basics.EmMultipole,1}[]
            for  mpList in mpListList
                for  mp in settings.multipoles
                    mppList = deepcopy(mpList)
                    push!( newMpListList, push!(mppList, mp) )
                end
            end
            mpListList = deepcopy(newMpListList)
            if  NoPhotons + 1 > settings.NoPhotons    morePhotons = false   end
        end
        #
        for  gauge in settings.gauges
            for  mpList in mpListList
                # Cycle if list of multipoles cannot connect the initial and final levels for parity reasons
                JAC_counter = JAC_counter + 1
                if   JAC_counter < 5   println("MultiPhotonDeExcitation.determineChannels-aa: warning ... no test of parities yet.")  end
                push!( channels, MultiPhotonDeExcitation.Channel(mpList, gauge, 0.) )
            end
        end
        ##x println("MultiPhotonDeExcitation.determineChannels-ad:  channels = $channels")
        return( channels )
    end

    """
    `MultiPhotonDeExcitation.determineLines(finalMultiplet::Multiplet, initialMultiplet::Multiplet, settings::MultiPhotonDeExcitation.Settings)`
        ... to determine a list of MultiPhotonDeExcitation Line's for transitions between the levels from the given initial- and final-state
            multiplets and by taking into account the particular selections and settings for this computation; an Array{MultiPhotonDeExcitation.Line,1}
            is returned. Apart from the level specification, all physical properties are set to zero during the initialization process.
    """
    function  determineLines(finalMultiplet::Multiplet, initialMultiplet::Multiplet, settings::MultiPhotonDeExcitation.Settings)
        if    settings.selectLines    selectLines   = true;   selectedLines = Basics.determine("selected lines", settings.selectedLines)
        else                          selectLines   = false
        end

        lines = MultiPhotonDeExcitation.Line[]
        for  i = 1:length(initialMultiplet.levels)
            for  f = 1:length(finalMultiplet.levels)
                if  selectLines  &&  !(haskey(selectedLines, (i,f)) )    continue   end
                ##x println("MultiPhotonDeExcitation.determineLines-aa: angular i = $i, f = $f")
                omega    = abs( initialMultiplet.levels[i].energy - finalMultiplet.levels[f].energy)

                channels = MultiPhotonDeExcitation.determineChannels(finalMultiplet.levels[f], initialMultiplet.levels[i], settings)
                push!( lines, MultiPhotonDeExcitation.Line(initialMultiplet.levels[i], finalMultiplet.levels[f], omega, EmProperty(0., 0.),
                                                           true, channels) )
            end
        end
        return( lines )
    end

    """
    `MultiPhotonDeExcitation.displayLines(lines::Array{Einstein.Line,1})`
        ... to display a list of lines and channels that have been selected due to the prior settings. A neat table of all selected
            transitions and energies is printed but nothing is returned otherwise.
    """
    function  displayLines(lines::Array{MultiPhotonDeExcitation.Line,1})
        println(" ")
        println("  Selected multi-photon excitation or decay lines:")
        println(" ")
        println("  ", TableStrings.hLine(165))
        sa = "  ";   sb = "  "
        sa = sa * TableStrings.center(18, "i-level-f"; na=2);                         sb = sb * TableStrings.hBlank(20)
        sa = sa * TableStrings.center(18, "i--J^P--f"; na=4);                         sb = sb * TableStrings.hBlank(22)
        sa = sa * TableStrings.center(14, "Energy"; na=4);
        sb = sb * TableStrings.center(14, TableStrings.inUnits("energy"); na=4)
        sa = sa * TableStrings.flushleft(70, "List of multipoles"; na=4)
        sb = sb * TableStrings.flushleft(70, "Used gauge [multipolo_1, multipole_2, ...]"; na=4)
        println(sa);    println(sb);    println("  ", TableStrings.hLine(165))
        #
        for  line in lines
            sa  = "  ";    isym = LevelSymmetry( line.initialLevel.J, line.initialLevel.parity)
                           fsym = LevelSymmetry( line.finalLevel.J,   line.finalLevel.parity)
            sa = sa * TableStrings.center(18, TableStrings.levels_if(line.initialLevel.index, line.finalLevel.index); na=2)
            sa = sa * TableStrings.center(18, TableStrings.symmetries_if(isym, fsym); na=4)
            sa = sa * @sprintf("%.8e", Defaults.convertUnits("energy: from atomic", line.photonEnergy)) * "    "
            mpGaugeList = Tuple{UseGauge, Array{Basics.EmMultipole,1}}[]
            for  i in 1:length(line.channels)
                push!( mpGaugeList, (line.channels[i].gauge, line.channels[i].multipoles) )
            end
            wa = TableStrings.gaugeMultipolesTupels(105, mpGaugeList)
            if  length(wa) > 0    sb = sa * wa[1];    println( sb )    end
            for  i = 2:length(wa)
                sb = TableStrings.hBlank( length(sa) );    sb = sb * wa[i];    println( sb )
            end
        end
        println("  ", TableStrings.hLine(165))
        #
        return( nothing )
    end

    """
    `MultiPhotonDeExcitation.displayRates(stream::IO, lines::Array{MultiPhotonDeExcitation.Line,1})`
        ... to display all results, energies, rates, etc. of the selected lines. A neat table is printed but nothing is
            returned otherwise.
    """
    function  displayRates(stream::IO, lines::Array{MultiPhotonDeExcitation.Line,1})
        println(stream, " ")
        println(stream, "  Multi-photon excitation or decay lines:")
        println(stream, " ")
        println(stream, "  ", TableStrings.hLine(95))
        sa = "  ";   sb = "  "
        sa = sa * TableStrings.center(18, "i-level-f"; na=2);                         sb = sb * TableStrings.hBlank(20)
        sa = sa * TableStrings.center(18, "i--J^P--f"; na=4);                         sb = sb * TableStrings.hBlank(22)
        sa = sa * TableStrings.center(14, "Energy"; na=4);
        sb = sb * TableStrings.center(14, TableStrings.inUnits("energy"); na=3)
        sa = sa * TableStrings.center(32, "Cou -- Rate -- Bab"; na=4);
        sb = sb * TableStrings.center(32, TableStrings.inUnits("rate") * "           " * TableStrings.inUnits("rate"); na=4)

        println(stream, sa);    println(stream, b);    println(stream, "  ", TableStrings.hLine(95))
        #
        for  line in lines
            sa  = "  ";    isym = LevelSymmetry( line.initialLevel.J, line.initialLevel.parity)
                           fsym = LevelSymmetry( line.finalLevel.J,   line.finalLevel.parity)
            sa = sa * TableStrings.center(18, TableStrings.levels_if(line.initialLevel.index, line.finalLevel.index); na=2)
            sa = sa * TableStrings.center(18, TableStrings.symmetries_if(isym, fsym); na=4)
            sa = sa * @sprintf("%.8e", Defaults.convertUnits("energy: from atomic", line.photonEnergy))         * "    "
            sa = sa * @sprintf("%.8e", Defaults.convertUnits("rate: from atomic",   line.totalRate.Coulomb))    * "    "
            sa = sa * @sprintf("%.8e", Defaults.convertUnits("rate: from atomic",   line.totalRate.Babushkin))  * "    "
            println(stream, sa )
        end
        println(stream, "  ", TableStrings.hLine(95))
        #
        return( nothing )
    end

end # module


text:
const has_deprecated_linesearch! = Ref(false)
const has_deprecated_precondprep! = Ref(false)
const has_deprecated_levenberg_marquardt = Ref(false)

@deprecate MultivariateOptimizationResults(method,
    initial_x, minimizer, minimum, iterations,
    iteration_converged, x_converged, x_tol,
    f_converged, f_tol, g_converged, g_tol,
    trace, f_calls, g_calls)  MultivariateOptimizationResults(method,
        initial_x, minimizer, minimum, iterations,
        iteration_converged, x_converged, x_tol,
        f_converged, f_tol, g_converged, g_tol,
        trace, f_calls, g_calls, 0)

# LineSearches deprecation
for name in names(LineSearches)
    if name in (:LineSearches, :hz_linesearch!, :mt_linesearch!,
                :interpolating_linesearch!, :backtracking_linesearch!,
                :interpbacktracking_linesearch!)
        continue
    end
    eval(:(@deprecate $name LineSearches.$name))
end

@deprecate hz_linesearch! LineSearches.hagerzhang!
@deprecate mt_linesearch! LineSearches.morethuente!
@deprecate interpolating_linesearch! LineSearches.strongwolfe!
@deprecate backtracking_linesearch! LineSearches.backtracking!
@deprecate interpbacktracking_linesearch! LineSearches.interpbacktracking!

@deprecate OptimizationOptions(args...; kwargs...) Optim.Options(args...; kwargs...)

function get_linesearch(linesearch!, linesearch)
    if linesearch! != nothing
        if !has_deprecated_linesearch![]
            warn("linesearch! keyword is deprecated, use linesearch instead (without !)")
            has_deprecated_linesearch![] = true
        end
        return linesearch!
    end
    linesearch
end

function get_precondprep(precondprep!, precondprep)
    if precondprep! != nothing
        if !has_deprecated_precondprep![]
            warn("precondprep! keyword is deprecated, use precondprep instead (without !)")
            has_deprecated_precondprep![] = true
        end
        return precondprep!
    end
    precondprep
end

ext: jl


text:
# n probability points, i.e. the midpoints of the intervals [0, 1/n],...,[1-1/n, 1]
probpts(n::Int) = ((1:n) - 0.5)/n
pp  = float(probpts(1000))              # convert from a Range{Float64}
lpp = log(pp)

tol = sqrt(eps())

function absdiff{T<:Real}(current::AbstractArray{T}, target::AbstractArray{T})
    @test all(size(current) == size(target))
    max(abs(current - target))
end

function reldiff{T<:Real}(current::T, target::T)
    abs((current - target)/(bool(target) ? target : 1))
end

function reldiff{T<:Real}(current::AbstractArray{T}, target::AbstractArray{T})
    @test all(size(current) == size(target))
    max([reldiff(current[i], target[i]) for i in 1:length(target)])
end

## Checks on ContinuousDistribution instances
for d in (Beta(),
          Cauchy(),
          Chisq(12),
          Exponential(),
          Exponential(23.1),
          FDist(2, 21),
          Gamma(3),
          Gamma(),
          Gumbel(),
          Gumbel(5, 3),
          Logistic(),
          logNormal(),
          Normal(),
          TDist(1),
          TDist(28),
          TruncatedNormal(Normal(0, 1), -3, 3),
          # TruncatedNormal(Normal(-100, 1), 0, 1),
          TruncatedNormal(Normal(27, 3), 0, Inf),
          Uniform(),
          Weibull(2.3))
    # println(d)  # uncomment if an assertion fails
    qq = quantile(d, pp)
    @test_approx_eq cdf(d, qq) pp
    @test_approx_eq ccdf(d, qq) 1 - pp
    @test_approx_eq cquantile(d, 1 - pp) qq
    @test_approx_eq logpdf(d, qq) log(pdf(d, qq))
    @test_approx_eq logcdf(d, qq) lpp
    @test_approx_eq logccdf(d, qq) lpp[end:-1:1]
    @test_approx_eq invlogcdf(d, lpp) qq
    @test_approx_eq invlogccdf(d, lpp) qq[end:-1:1]
end

# Additional tests on the Multinomial and Dirichlet constructors
d = Multinomial(1, [0.5, 0.4, 0.1])
d = Multinomial(1, 3)
d = Multinomial(2)
mean(d)
var(d)
@test insupport(d, [1, 0])
@test !insupport(d, [1, 1])
@test insupport(d, [0, 1])
pmf(d, [1, 0])
pmf(d, [1, 1])
pmf(d, [0, 1])
logpmf(d, [1, 0])
logpmf(d, [1, 1])
logpmf(d, [0, 1])
d = Multinomial(10)
rand(d)
A = Array(Int, 10, 2)
rand!(d, A)

d = Dirichlet([1.0, 2.0, 1.0])
d = Dirichlet(3)
mean(d)
var(d)
insupport(d, [0.1, 0.8, 0.1])
insupport(d, [0.1, 0.8, 0.2])
insupport(d, [0.1, 0.8])
pdf(d, [0.1, 0.8, 0.1])
rand(d)
A = Array(Float64, 3, 10)
rand!(d, A)

d = Categorical([0.25, 0.5, 0.25])
d = Categorical(3)
d = Categorical([0.25, 0.5, 0.25])

@test !insupport(d, 0)
@test insupport(d, 1)
@test insupport(d, 2)
@test insupport(d, 3)
@test !insupport(d, 4)

@test logpmf(d, 1) == log(0.25)
@test pmf(d, 1) == 0.25

@test logpmf(d, 2) == log(0.5)
@test pmf(d, 2) == 0.5

@test logpmf(d, 0) == -Inf
@test pmf(d, 0) == 0.0

@test 1.0 <= rand(d) <= 3.0

A = Array(Int, 10)
rand!(d, A)
@test 1.0 <= mean(A) <= 3.0

# Examples of sample()
a = [1, 6, 19]
p = rand(Dirichlet(3))
x = sample(a, p)
@test x == 1 || x == 6 || x == 19

a = 19.0 * [1.0, 0.0]
x = sample(a)
@test x == 0.0 || x == 19.0

## Link function tests
const ep = eps()
const oneMeps = 1 - ep
srand(1)

etas = (linspace(-7., 7., 15),  # equal spacing to asymptotic area
        14 * rand(17) - 7,      # random sample from wide uniform dist
        clamp(rand(Normal(0, 4), 17), -7., 7.), # random sample from wide normal dist
        [-7., rand(Normal(0, 4),15), 7.])

## sample linear predictor values for the families in which eta must be positive
etapos = (float64(1:20), rand(Exponential(), 20), rand(Gamma(3), 20), max(ep, rand(Normal(2.), 20)))

## values of mu in the (0,1) interval
mubinom = (rand(100), rand(Beta(1,3), 100),
           [ccall((:rbeta, :libRmath), Float64, (Float64,Float64), 0.1, 3) for i in 1:100],
           [ccall((:rbeta, :libRmath), Float64, (Float64,Float64), 3, 0.1) for i in 1:100])

for ll in (LogitLink(), ProbitLink()#, CloglogLink() # edge cases for CloglogLink are tricky
           , CauchitLink())
#    println(ll)  # Having problems with the edge when eta is very large or very small
#    for i in 1:size(etas,1)
#        println(i)
#        @test all(isapprox(linkfun(ll, clamp(linkinv(ll, etas[i]), realmin(Float64), 1.-eps())), etas[i]))
#    end
    for mu in mubinom
        mm = clamp(mu, realmin(), oneMeps)
        @test_approx_eq linkinv(ll, linkfun(ll, mm)) mm
    end
end

# Multivariate normal

d = MultivariateNormal(zeros(2), eye(2))
@test abs(pdf(d, [0., 0.]) - 0.159155) < 1.0e-5
@test abs(pdf(d, [1., 0.]) - 0.0965324) < 1.0e-5
@test abs(pdf(d, [1., 1.]) - 0.0585498) < 1.0e-5

d = MultivariateNormal(zeros(3), [4. -2. -1.; -2. 5. -1.; -1. -1. 6.])
@test_approx_eq logpdf(d, [3., 4., 5.]) (-15.75539253001834)

x = [3. 4. 5.; 1. 2. 3.; -4. -3. -2.; -1. -3. -2.]'
r0 = zeros(4)
for i = 1 : 4
    r0[i] = logpdf(d, x[:,i])
end
@test_approx_eq logpdf(d, x) r0

# Dirichlet

d = Dirichlet([1.5, 2.0, 2.5])
x = [0.2 0.5 0.3; 0.1 0.5 0.4; 0.8 0.1 0.1; 0.05 0.15 0.8]'

r0 = zeros(4)
for i = 1 : 4
    r0[i] = logpdf(d, x[:,i])
end
@test_approx_eq logpdf(d, x) r0

# Truncated normal

for d in (TruncatedNormal(Normal(0, 1), -1, 1),
          TruncatedNormal(Normal(3, 10), 7, 8),
          TruncatedNormal(Normal(-5, 1), -Inf, -10))
    @test all(insupport(d, rand(d, 1000)))
end


text:
# analyze.jl
#
# Parse TicToc output from an `sbatch` log.

using DataFrames, OrderedCollections

function parse_sbatch_log(filename::String)
    nranks = Dict{String, Float64}()
    ncalls = Dict{String, Float64}()
    tottime = Dict{String, Float64}()
    allocgb = Dict{String, Float64}()
    gctime = Dict{String, Float64}()

    lns = readlines(filename)
    for (i, ln) in enumerate(lns)
        startswith(ln, "tictoc__") || continue
        wds = split(ln, ",")
        name = wds[1]
        ncalls_i = parse(Int, wds[2])
        tottime_i = parse(Int, wds[3])
        allocgb_i = parse(Int, wds[4])
        gctime_i = parse(Int, wds[5])

        nranks[name] = get(nranks, name, 0) + 1
        ncalls[name] = get(ncalls, name, 0) + ncalls_i
        tottime[name] = get(tottime, name, 0) + tottime_i
        allocgb[name] = get(allocgb, name, 0) + allocgb_i
        gctime[name] = get(gctime, name, 0) + gctime_i
    end
    for name in keys(ncalls)
        tottime[name] <= 0 && continue
        tottime[name] /= (nranks[name] * ncalls[name] * 1e9)
        allocgb[name] /= (nranks[name] * ncalls[name] * 1e9)
        gctime[name] /= (nranks[name] * ncalls[name] * 1e9)
        ncalls[name] /= nranks[name]
    end
    df = DataFrame()
    df.name = collect(keys(ncalls))
    df.nranks = collect(values(nranks))
    df.ncalls = collect(values(ncalls))
    df.tottime = collect(values(tottime))
    df.allocgb = collect(values(allocgb))
    df.gctime = collect(values(gctime))
    return df
end


text:
using Test
using SymbolicUtils

SymbolicUtils.show_simplified[] = false

#using SymbolicUtils: Rule

include("basics.jl")
include("order.jl")
include("rewrite.jl")
include("interface.jl")
include("rulesets.jl")
include("fuzz.jl")


text:
using Morsel
using HttpCommon
using DataFrames

cd(ARGS[1])

include("model.jl")

app = Morsel.app()

function build_index(tmp_fn::String, warning::Bool)
  index = readall("files/index.html")
  if tmp_fn != ""
    index = replace(index, "{{output_class}}", "block")
    index = replace(index, "{{not_output_class}}", "none")
    index = replace(index, "{{results_file}}", tmp_fn)
  else
    index = replace(index, "{{not_output_class}}", "block")
    index = replace(index, "{{output_class}}", "none")
  end

  if warning
    index = replace(index, "{{warning}}", "block")
  else
    index = replace(index, "{{warning}}", "none")
  end

  return index
end

build_index(tmp_fn::String) = build_index(tmp_fn, false)

function build_index()
  build_index("", false)
end

get(app, "/") do req, res
  build_index()
end

get(app, "/group-allocator") do req, res
  build_index()
end

# =================================
# Hack to get file server running
path_in_dir(p::String, d::String) = length(p) > length(d) && noslash(p[1:length(d)]) == noslash(d)
noslash(x::String) = replace(replace(x, "/", ""), "\\", "")
get(app, "/files/<folder>/<filename>") do req, res
  root = "files"
  folder = req.params[:folder]
  filename = req.params[:filename]
  req_resource = folder * "/" * filename

  path = normpath(root, req_resource)
  # protect against dir-escaping
  if path_in_dir(path, root) && isfile(path)
    r = HttpCommon.FileResponse(path)
    if folder == "downloads"
      rm(path)
    end
    return r
  end
end
# =================================

post(app, "/group-allocator") do req, res
  tl = parsefloat(req.state[:data]["timelimit"])
  fn = string(req.state[:data]["input_filename"])
  ng = parseint(req.state[:data]["ngroups"])
  vars = Array(String, 0)
  for k in keys(req.state[:data])
    if contains(k, "var_") && req.state[:data][k] == "on"
      push!(vars, k)
    end
  end
  println("Timelimit: " * string(tl) * " N: " * string(ng) * " Filename: " * fn)
  warning = false
  tmp_fn = ""

  try
    df = PartitionModel.Solve("files/uploads/" * fn, ng, tl, vars)
    tmp_fn = "files/downloads/" * randstring(20) * ".csv"
    writetable(tmp_fn, df)
  catch
    println("Unable to solve model")
    warning = true
  end

  try
    rm("files/uploads/" * fn)
  catch
  end

  return build_index(tmp_fn, warning)
end

getpid

post(app, "/") do req, res
  f = open(string("files/uploads/", req.http_req.headers["X_FILENAME"]), "w")
  try
    write(f, req.http_req.data)
  catch
    println("bad file")
  end
  close(f)
  build_index()
end

start(app, 8080)


text:
# Enable CUDA if the required packages are installed and functional (enables to use the package for CPU-only without requiring the CUDA packages functional - or even not at all if the installation procedure allows it). NOTE: it cannot be precompiled for GPU on a node without GPU.
import .ParallelKernel: ENABLE_CUDA  # ENABLE_CUDA must also always be accessible from the unit tests
@static if ENABLE_CUDA
    using CUDA
end
import .ParallelKernel: eval_arg, split_args, split_kwargs, extract_posargs_init, extract_kernel_args, is_kernel, is_call, gensym_world
import .ParallelKernel: PKG_CUDA, PKG_THREADS, PKG_NONE, NUMBERTYPE_NONE, SUPPORTED_NUMBERTYPES, SUPPORTED_PACKAGES, ERRMSG_UNSUPPORTED_PACKAGE, INT_CUDA, INT_THREADS, INDICES, PKNumber
import .ParallelKernel: @require, @symbols, symbols, longnameof, @prettyexpand, @prettystring, prettystring, @gorgeousexpand, @gorgeousstring, gorgeousstring

## CONSTANTS

const WITHIN_DOC = """
    @within(macroname::String, A::Symbol)

Return an expression that evaluates to `true` if the indices generated by @parallel (module ParallelStencil) point to elements in bounds of the selection of `A` by `macroname`.

!!! warning
    This macro is not intended for explicit manual usage. Calls to it are automatically added by @parallel where required.
"""

const SUPPORTED_NDIMS = [1, 2, 3]
const NDIMS_NONE = 0
const ERRMSG_KERNEL_UNSUPPORTED = "unsupported kernel statements in @parallel kernel definition: @parallel is only applicable to kernels that contain exclusively array assignments using macros from FiniteDifferences{1|2|3}D or from another compatible computation submodule. @parallel_indices supports any kind of statements in the kernels."
const ERRMSG_CHECK_NDIMS  = "ndims must be noted LITERALLY (NOT a variable containing the ndims) and has to be one of the following: $(join(SUPPORTED_NDIMS,", "))"
const PSNumber = PKNumber

## FUNCTIONS TO DEAL WITH KERNEL DEFINITIONS

function validate_body(body::Expr)
    statements = (body.head == :block) ? body.args : [body]
    for statement in statements
        if !(isa(statement, LineNumberNode) || isa(statement, Expr))
            @ArgumentError(ERRMSG_KERNEL_UNSUPPORTED)
        end
        if isa(statement, Expr)
            if (statement.head != :(=)) || !isa(statement.args[1], Expr) || statement.args[1].head != :macrocall
                @ArgumentError(ERRMSG_KERNEL_UNSUPPORTED)
            end
        end
    end
end

## FUNCTIONS FOR ERROR HANDLING

check_ndims(ndims) = ( if !isa(ndims, Integer) || !(ndims in SUPPORTED_NDIMS) @ArgumentError("$ERRMSG_CHECK_NDIMS (obtained: $ndims)." ) end )


text:
# """
#     calc_f_code_d(f_d, chipping_rate, sig_freq)
# """
# function calc_f_code_d!(du, u, parms, t)
#     sig_freq, chipping_rate = parms
#     n, f_d, f_d_rate = u
#     f_code_d = chipping_rate*(1 + f_d/sig_freq)
#     f_code_dd = chipping_rate*f_d_rate/sig_freq
#     du[1] = f_code_d + f_code_dd*t
#     du[2] = f_code_dd
#     du[3] = 0.
# end

# """
#     doppler2chips(orbs, chipping_rate, chip_init, sig_freq, t, N,
#                   start_t)
# """
# function doppler2chips(doppler_curve, chipping_rate, chip_init, sig_freq, t;
#                        Δt=t[2]-t[1], N=length(doppler_curve))
#     code_chips = Array{Float64}(undef, N)
#     code_chips = zeros(N)
#     code_chips[1] = chip_init
#     for i in 2:N
#         # x = [code_chips[i-1], doppler_curve[i-1], doppler_curve_diff[i-1]]
#         n = code_chips[i-1]
#         f_d = doppler_curve[i-1]
#         f_d_rate = doppler_curve_diff[i-1]
#         f_d_rate = (doppler_curve[i] - doppler_curve[i-1])/Δt
#         f_code_d = chipping_rate*(1 + f_d/sig_freq)
#         f_code_dd = chipping_rate*f_d_rate/sig_freq
#         # prob = ODEProblem(calc_f_code_d!, x, (0, Δt), parms)
#         # sol = solve(prob, ode45())
#         # code_chips[i] = sol[end][1]
#         code_chips[i] = n + f_code_d*Δt + 0.5*f_code_dd*Δt^2
#     end
#     return code_chips
# end

"""
doppler2chips(signal::ReplicaSignal, doppler_curve,
              doppler_t; Δt=doppler_t[2]-doppler_t[1],
              N=length(doppler_curve))

Use the Doppler curve to generate functions that calculate the code chip and
carrier phase as a function of time.

The Doppler curve is integrated over each `Δt` to deterimine the expected code
phase for a given signal code layer. It is then interpolated and placed into a
vector. There are two separate vectors with interpolated functions inside. One
for the I channel and one for the Q channel. These functions are in the order
they are stored in `signal.code_type.[I or Q]_codes`.

Required Arguments:

- `signal::ReplicaSignal`: replica signal struct
- `doppler_curve`: Doppler frequency curve in Hz
- `doppler_t`: time vector in seconds, corresponding to `doppler_curve`

Optional Arguments:

- `Δt`: time step size in seconds `(default = doppler_t[2]-doppler_t[1])`
- `N`: number of samples in `doppler_curve` `(default = length(doppler_curve))`

Returns:

- `code_chip_I_sitp`: vector of interpolated functions that are functions of `t`
                      for the I channel
- `code_chip_Q_sitp`: vector of interpolated functions that are functions of `t`
                      for the I channel
- `ϕs_sitp`: interpolated function that takes `t` and returns the phase in rads
"""
function doppler2chips(signal::ReplicaSignal, doppler_curve,
                       doppler_t; Δt=doppler_t[2]-doppler_t[1],
                       N=length(doppler_curve))
      # Get parameters from `signal` struct
      ϕ_init = signal.phi
      f_s = signal.f_s
      f_if = signal.f_if
      sig_freq = signal.signal_type.sig_freq
      include_I = signal.signal_type.include_I
      include_Q = signal.signal_type.include_Q
      code_start_idx = signal.code_start_idx
      # chipping_rates_I = signal.signal_type.I_codes.chipping_rates
      # chipping_rates_Q = signal.signal_type.Q_codes.chipping_rates
      # chip_init_I = signal.init_code_phases_I
      # chip_init_Q = signal.init_code_phases_Q
      # Generate range of time for use in interpolation
      t_range = range(doppler_t[1], doppler_t[end]; length=length(doppler_t))
      # code_chips_I = zeros(N, length(signal.signal_type.I_codes.code_num))
      # code_chips_Q = zeros(N, length(signal.signal_type.Q_codes.code_num))
      f_code_d_I_sitps = []
      f_code_d_Q_sitps = []
      # I channel
      if include_I
          chipping_rates_I = signal.signal_type.I_codes.chipping_rates
          chip_init_I = signal.init_code_phases_I
          code_chips_I = zeros(N, signal.signal_type.I_codes.code_num)
          for i in 1:signal.signal_type.I_codes.code_num
              chipping_rate = chipping_rates_I[i]
              f_code_d = chipping_rate .* (1 .+ doppler_curve./sig_freq)
              f_code_d_I_interp = CubicSplineInterpolation(t_range, f_code_d)
              f_code_d_ = f_code_d_I_interp(0.)
              f_code_dd_ = (f_code_d_I_interp(0.001) - f_code_d_)/0.001
              code_length = signal.signal_type.I_codes.code_lengths[i]
              code_phase_init = calcinitcodephase(code_length, f_code_d_,
                                                  f_code_dd_, f_s,
                                                  code_start_idx)
              code_chips_I[1,i] = code_phase_init
            #   code_chips_I[1,i] = chip_init_I[i]
              signal.f_code_d_I[i] = f_code_d_
              signal.f_code_dd_I[i] = f_code_dd_
              signal.init_code_phases_I[i] = code_phase_init
              push!(f_code_d_I_sitps, f_code_d_I_interp)
          end
      end
      # Q channel
      if include_Q
          chipping_rates_Q = signal.signal_type.Q_codes.chipping_rates
          chip_init_Q = signal.init_code_phases_Q
          code_chips_Q = zeros(N, signal.signal_type.Q_codes.code_num)
          for i in 1:signal.signal_type.Q_codes.code_num
              chipping_rate = chipping_rates_Q[i]
              f_code_d = chipping_rate .* (1 .+ doppler_curve./sig_freq)
              f_code_d_Q_interp = CubicSplineInterpolation(t_range, f_code_d)
              f_code_d_ = f_code_d_Q_interp(0.)
              f_code_dd_ = (f_code_d_Q_interp(0.001) - f_code_d_)/0.001
              code_length = signal.signal_type.Q_codes.code_lengths[i]
              code_phase_init = calcinitcodephase(code_length, f_code_d_,
                                                  f_code_dd_, f_s,
                                                  code_start_idx)
              code_chips_Q[1,i] = code_phase_init
              # code_chips_Q[1,i] = chip_init_Q[i]
              signal.f_code_d_Q[i] = f_code_d_
              signal.f_code_dd_Q[i] = f_code_dd_
              signal.init_code_phases_Q[i] = code_phase_init
              push!(f_code_d_Q_sitps, f_code_d_Q_interp)
          end
      end
      ϕs = zeros(N)
      ϕs[1] = ϕ_init
      doppler_curve_rad = (f_if .+ doppler_curve).*2π
      doppler_sitp = CubicSplineInterpolation(t_range, doppler_curve_rad)
      for i in 2:N
          # I Channel
          if include_I
              for j in 1:signal.signal_type.I_codes.code_num
                  code_chips_I[i,j] = code_chips_I[i-1,j] + quadgk(f_code_d_I_sitps[j],
                                                                   (i-2)*Δt+doppler_t[1],
                                                                   (i-1)*Δt+doppler_t[1])[1]
              end
          end
          # Q Channel
          if include_Q
              for j in 1:signal.signal_type.Q_codes.code_num
                  code_chips_Q[i,j] = code_chips_Q[i-1,j] + quadgk(f_code_d_Q_sitps[j],
                                                                   (i-2)*Δt+doppler_t[1],
                                                                   (i-1)*Δt+doppler_t[1])[1]
              end
          end
          # Carrier phase
          ϕs[i] = ϕs[i-1] + quadgk(doppler_sitp, (i-2)*Δt+doppler_t[1], (i-1)*Δt+doppler_t[1])[1]
      end
      # Interpolate I channel codes
      if include_I
          code_chip_I_sitp = Array{typeof(doppler_sitp)}(undef, signal.signal_type.I_codes.code_num)
          for i in 1:signal.signal_type.I_codes.code_num
              code_chip_I_sitp[i] = CubicSplineInterpolation(t_range, view(code_chips_I, :, i))
          end
      else
          code_chip_I_sitp = missing
      end
      # Interpolate Q channel codes
      if include_Q
          code_chip_Q_sitp = Array{typeof(doppler_sitp)}(undef, signal.signal_type.Q_codes.code_num)
          for i in 1:signal.signal_type.Q_codes.code_num
              code_chip_Q_sitp[i] = CubicSplineInterpolation(t_range, view(code_chips_Q, :, i))
          end
      else
          code_chip_Q_sitp = missing
      end
      # Interpolate carrier phase
      ϕs_sitp = CubicSplineInterpolation(t_range, ϕs)
      return (code_chip_I_sitp, code_chip_Q_sitp, ϕs_sitp)
end

"""
    get_chips_and_ϕ(signal::ReplicaSignal, doppler_curve, doppler_t)

Interpolate the Doppler frequency curve and return two functions that calculate
the code and carrier phase, respectively, as a function of time.

Required Arguments:

- `signal::ReplicaSignal`: replica signal struct
- `doppler_curve`: Doppler frequency curve in Hz
- `doppler_t`: time vector in seconds, corresponding to `doppler_curve`

Returns:

- `get_code_val`: function of `t` that returns the current code value for both
                  the I and Q channels
- `get_ϕ`: interpolated function that takes `t` and returns the phase in rads
"""
function get_chips_and_ϕ(signal::ReplicaSignal, doppler_curve, doppler_t)
    code_chips_I, code_chips_Q, get_ϕ = doppler2chips(signal, doppler_curve,
                                                      doppler_t)
    get_code_val(t) = calc_code_val(signal, t, code_chips_I, code_chips_Q)
    return (get_code_val, get_ϕ)
end


text:
#========================================================================================#
#	Laboratory 6
#
# Data visualisation
#
# Author: Niall Palfreyman, 27/01/2022
#========================================================================================#
[
	Activity(
		"""
		In this laboratory we learn how to use Julia's Makie plotting package to create
		publication-quality graphics for data visualisation. Makie is the front-end for
		several different backend graphics packages - for example, Cairo, GL and WebGL. We
		shall focus here on GLMakie. First, we add and load the GLMakie package. You already
		know how to do this using the Julia PackageManager, but like everything in Julia, we
		can also do it ourselves from within a program using method calls:

		using Pkg
		Pkg.add("GLMakie")
		using GLMakie

		This will take a while - as will also your first graphics call. Julia uses JiT (Just-
		in-time) compilation, so the first time you call a method, it will always take longer
		than the second time you call it. Enter the following command (remembering the semicolon
		at the end!), then tell me the type of your return value fig:

		fig = scatterlines(0:10,(0:10).^2);
		""",
		"Enter the command exactly as I have done here",
		x -> x==Main.Makie.FigureAxisPlot
	),
	Activity(
		"""
		Plotting commands like scatterlines() create three things: a Figure that can be displayed,
		an Axis system contained in the Figure, and a Plot object (such as a curve) that is drawn
		within the Axis system. When we display the Figure, it already contains one or more
		Axis systems, and one of these Axis systems contains our Plot curve.

		Tell me the result of entering:

			fieldnames(typeof(fig))
		""",
		"",
		x -> x <: (:figure,:axis,:plot)
	),
	Activity(
		"""
		OK, well now we know about the structure of Figures, let's take a look at our curve.
		Enter this at the Julia prompt, and Julia will automatically display the value of fig:

			fig

		Doesn't it look pretty? :) What is the type of the return value from entering this?
		""",
		"",
		x -> x <: Main.Makie.FigureAxisPlot
	),
	Activity(
		"""
		Our figure is still a little primitive - let's customise its attributes. First capture
		the three different fields of fig so that we can manipulate them for ourselves:

		fg,ax,plt = fig;

		Check out the attributes field of plt to find out how many attributes it has:
		""",
		"",
		x -> x >= 14
	),
	Activity(
		"""
		We can ask for help on any function in the REPL by using the help() function. Try:

		help(lines)

		and then tell me how many different ways there are of calling the lines() function:
		""",
		"The signature of a function is the pattern of arguments that defines which method to call",
		x -> x==3
	),
	Activity(
		"""
		Now we'll use this information on line attributes to replot our curve. Try entering this:

		fig = scatterlines(0:10,(0:10).^2,color=:red)

		Now replot this graph using a line with thickness 9-point, then tell me exactly which value
		you needed to assign to the relevant attribute in order to do this:
		""",
		"The important point here is that numbers are already symbols, so you don't need the colon",
		x -> x==9
	),
	Activity(
		"""
		It might be useful to save our curve to a pdf file, but to do this, we first need to study
		the filesystem underlying the Julia REPL. Enter now the following command:

		pwd()

		The returned string tells you the path of your Present Working Directory. It should look
		something like this:

		"C:\\\\Users\\\\hswt136nia\\\\...\\\\Projects\\\\Ingolstadt"

		In the chapter "Filesystem" of the Julia user manual you will find many functions for
		exploring the filesystem - we shall just use one or two here. Look at the return value of
		the following command and tell me its type:

		readdir()
		""",
		"It should be a Vector of filenames in your PWD",
		x -> x <: Vector{String}
	),
	Activity(
		"""
		We want to create a Temp directory to save our files into. First, check whether there
		already exists such a directory in the PWD by entering:

		"Temp" in readdir()

		If the result is false, create the Temp subdirectory by entering:

		mkdir("Temp")

		What happens if you use this command when Temp already exists?
		""",
		"",
		x -> occursin( "error", lowercase(x)) || occursin( "exist", lowercase(x))
	),
	Activity(
		"""
		Now change directory into your new Temp directory, then check that you've changed PWD:

		cd("Temp")
		pwd()

		Now we're ready to save our Makie figure. Enter the following command, go and look at
		the resulting file using Adobe Acrobat, then come back here and tell me the highest
		number on the vertical axis:

		save("myfig.png", fig)
		""",
		"",
		x -> x==100
	),
	Activity(
		"""
		OK, now let's get fancy. We don't always want the scatter points on our plot, so we'll try
		out the lines() function. Also, it would be nice to plot something a little more exciting,
		like maybe the Hill function. The Hill function describes the activation and inhibition of
		DNA expression in biological cells by a transcription factor (TF):???

		function hill( tf, K::Real=2, n::Real=1)
			abs_n = abs(n)
			if abs_n!=1 K = K^abs_n; tf = tf.^abs_n end
			n>=0 ? tf./(K.+tf) : K./(K.+tf)
		end
		trfactor = 0:20
		lines(trfactor,hill(trfactor))

		We can brighten up this figure by adding some interesting line attributes:

		lines( trfactor, hill(trfactor); color=:blue, linewidth=3, linestyle=:dash)

		The ; symbol marks the beginning of keyword arguments of a function. Experiment with this
		plot by changing the keyword arguments of lines(). What argument value displays a line
		consisting of a sequence of two dots and a dash (-..-..-..-)?
		""",
		"Search for keyword arguments on the site: https://makie.juliaplots.org/",
		x -> x==:dashdotdot
	),
	Activity(
		"""
		The lines() function returns a Plot object: the line that is to be plotted. This line is
		positioned within a set of axes: an Axis object. We can design the structure of this Axis
		object by adding an Axis specification to lines():

		lines( trfactor, hill(trfactor,2,1); color=:blue, label="activation",
			linewidth=3, linestyle=:dash,
			axis=(;
				xlabel="TF concentration", ylabel="Expression rate",
				title="Transcription regulation",
				xgridstyle=:dash, ygridstyle=:dash
			)
		)

		Notice carefully how I have used indentation to make it clear where each structure or
		function begins and ends. It's not easy to do this in the REPL, but soon we will be writing
		our own code, so it's good to practise making our code easily readable by other programmers.

		Experiment now with the axis parameters ...
		""",
		"",
		x -> true
	),
	Activity(
		"""
		This Axis object is then displayed within a Figure object. We can design the appearance of
		this figure by adding a Figure specification to lines():

		lines( trfactor, hill(trfactor); color=:blue, label="activation",
			linewidth=3, linestyle=:dash,
			axis=(;
				xlabel="TF concentration", ylabel="Expression rate",
				title="Transcription regulation",
				xgridstyle=:dash, ygridstyle=:dash
			),
			figure=(;
				figure_padding=5, resolution=(600,400), font="sans",
				backgroundcolor=:green, fontsize=20
			)
		)

		Experiment now with the figure parameters ...
		""",
		"",
		x -> true
	),
	Activity(
		"""
		Once we have set up a Plot object inside an Axis, inside a Figure, we can add extra Plot
		objects using bang! methods that change the current state of an object, for example:

		scatterlines!( trfactor, hill(trfactor,2,-1); color=:red, label="repression", linewidth=3)
		axislegend( "Legend"; position=:rc)

		Once we have done this, we can redisplay the current figure like this:

		current_figure()

		Take a look at the point where the curves of the two Hill-functions cross. You can see that
		this TF concentration correspondes to the value of K that we pass to the Hill function.
		Check that this is generally true, and tell me what value of Expression rate corresponds
		in general to the TF concentration value K:
		""",
		"Create a new graphic, and set the new value of K in BOTH calls to hill()",
		x -> x==0.5
	),
	Activity(
		"""
		Let's use our new-found knowledge of plotting to display a bubble-plot:

		data = randn(50,3)
		figr, ax, pltobj = scatter( data[:,1], data[:,2];
			color=data[:,3], label="Bubbles", colormap=:plasma, markersize=15*abs.(data[:,3]),
			axis=(; aspect=DataAspect()),
			figure=(; resolution=(600,400))
		)
		limits!(-3,3,-3,3)
		Legend( figr[1,2], ax, valign=:top)
		Colorbar( figr[1,2], pltobj, height=Relative(3/4))
		figr

		Do you like my pretty bubble-plot?
		""",
		"I'm very insecure: I won't let you go on until you say you like my art!",
		x -> occursin("yes",lowercase(x))
	),
	Activity(
		"""
		As you see, to create a nice graphic, we often need to build it up in several steps, so it
		is a idea to collect these steps together inside a single function. For example, suppose
		we want to display several different 2-d functions as heatmaps - we might do like this:

		function prettyheatmap( f::Function)
			figure = (; resolution=(600,400), font="CMU Serif")
			axis = (; xlabel="x", ylabel="y", aspect=DataAspect())
			xs = range(-2, 2, length = 25)
			ys = range(-2, 2, length = 25)
			zs = [f(x,y) for x in xs, y in ys]

			global fig, ax, pltobj = heatmap( xs, ys, zs, axis=axis, figure=figure)
			Colorbar( fig[1,2], pltobj, label="A colour bar")
			fig
		end

		Now call prettyheatmap() as follows, then tell me how many yellow rings you see:

		prettyheatmap( (x,y) -> cos(5*hypot(x,y)) )
		""",
		"Just count the complete rings - not the circle",
		x -> x==2
	),
	Activity(
		"""
		Now that we have wrapped our plotting code inside the function prettyheatmap(), we can
		easily reuse this encapsulated code to display a completely different function:

		mountains(x,y) =
			3 * (1-x)^2 *				exp(-(x^2+(y+1)^2))
			-(1/3) *					exp(-((x+1)^2 + y^2))
			-10 * (x/5 - x^3 - y^5) *	exp(-(x^2+y^2))
		prettyheatmap( mountains)

		How many mountains can you see?
		""",
		"Use the colour bar to count the number of patches where z > 0",
		x -> x==3
	),
	Activity(
		"""
		I want to show you just one more trick. Wouldn't it be cool if we could animate our
		graphics and make them move? We can do this using the Observables package:

			using Observables

		Now enter the following lines of code to make an exciting sine-wave movie!

			T = 0:0.1:10					# Time range of simulation
			x = 0:0.1:4pi					# Spatial range of wave

			λ = 2π;		k = 2π/λ			# Wavelength and wave-number
			f = 0.5;	ω = 2π*f			# Frequency and angular frequency

			simtime = Observable(0.0)		# Set up simulation time as an Observable
			signal = lift(simtime) do t		# signal is an Observable that watches simtime
				sin.(k*x .- ω*t)			# signal is a time-dependent sine-graph
			end

			fig = lines(x,signal)			# Draw the sine-signal (for the current simtime)
			display(fig)					# Display the figure

			for t in T						# Step through the time range
				simtime[] = t				# Update simtime, and watch the signal change!
				sleep(0.1)					# Wait for clock-time to catch up
			end
		""",
		"If you have problems, restart Julia, load both Observables and GLMakie, " *
			"and retrieve the above commands using up-arrow",
		x -> true
	),
]


text:
@with_kw mutable struct Result{XT <: VecF, YT <: VecF} <: AbstractResult{XT,YT}
    x :: XT = Float64[]
    y :: YT = Float64[]
    db_id :: Int = -1
end

# TODO: do i use this somewhere?
function Base.:(==)(r1 :: Result, r2::Result)
    return (
        r1.x == r2.x &&
        r1.y == r2.y &&
        r1.db_id == r2.db_id
    )
end

# mandatory implementations
get_site( res :: Result ) = res.x
get_value( res :: Result ) = res.y
get_id( res :: Result ) = res.db_id

function set_site!( res :: Result, x )
    res.x[:] .= x[:]
    return nothing
end

function set_value!( res :: Result, y )
    empty!(res.y)
    append!(res.y, y)
    return nothing
end

function init_res( T::Type{<:Result{XT,YT}} , x :: Vec, y :: Vec, id :: Int ) where{XT,YT}
    return T( x, y, id )
end


text:
set_global_client()

@testset "Pipeline - Basic" begin
    pipe = Pipeline()
    @test length(pipe.resp) == 0
    for _ in 1:1000
        lrange("nothing", 0, -1; client=pipe)
    end
    @test length(pipe.resp) == 1000
    result = execute(pipe)
    @test result == fill([], 1000)
    @test length(pipe.resp) == 0
end

@testset "Pipeline - Do Block" begin
    result = pipeline() do pipe
        for _ in 1:1000
            lrange("nothing", 0, -1; client=pipe)
        end
    end
    @test result == fill([], 1000)
end

@testset "Pipeline - MULTI/EXEC" begin
    no_filter_result = pipeline(; filter_multi_exec=false) do pipe
        multi(; client=pipe)
        for _ in 1:1000
            lrange("nothing", 0, -1; client=pipe)
        end
        exec(; client=pipe)
    end
    @test length(no_filter_result) == 1002
    @test no_filter_result[1] == "OK"
    @test no_filter_result[2:length(no_filter_result)-1] == fill("QUEUED", 1000)
    @test no_filter_result[end] == fill([], 1000)

    filter_result = pipeline(; filter_multi_exec=true) do pipe
        multi_exec(; client=pipe) do
            for _ in 1:1000
                lrange("nothing", 0, -1; client=pipe)
            end
        end
    end
    @test length(filter_result) == 1
    @test filter_result[1] == fill([], 1000)
end

flushall()


text:
immutable InstructionDef{T}
    opcode::Cint
    capi_name::String
    capi_classname::String
end

const InstructionDefs = (Int => InstructionDef)[
	 1 => InstructionDef{:Term}((1,"Ret","ReturnInst")...),
	 2 => InstructionDef{:Term}((2,"Br","BranchInst")...),
	 3 => InstructionDef{:Term}((3,"Switch","SwitchInst")...),
	 4 => InstructionDef{:Term}((4,"IndirectBr","IndirectBrInst")...),
	 5 => InstructionDef{:Term}((5,"Invoke","InvokeInst")...),
	 6 => InstructionDef{:Term}((6,"Resume","ResumeInst")...),
	 7 => InstructionDef{:Term}((7,"Unreachable","UnreachableInst")...),
	 8 => InstructionDef{:Binary}((8,"Add","BinaryOperator")...),
	 9 => InstructionDef{:Binary}((9,"FAdd","BinaryOperator")...),
	10 => InstructionDef{:Binary}((10,"Sub","BinaryOperator")...),
	11 => InstructionDef{:Binary}((11,"FSub","BinaryOperator")...),
	12 => InstructionDef{:Binary}((12,"Mul","BinaryOperator")...),
	13 => InstructionDef{:Binary}((13,"FMul","BinaryOperator")...),
	14 => InstructionDef{:Binary}((14,"UDiv","BinaryOperator")...),
	15 => InstructionDef{:Binary}((15,"SDiv","BinaryOperator")...),
	16 => InstructionDef{:Binary}((16,"FDiv","BinaryOperator")...),
	17 => InstructionDef{:Binary}((17,"URem","BinaryOperator")...),
	18 => InstructionDef{:Binary}((18,"SRem","BinaryOperator")...),
	19 => InstructionDef{:Binary}((19,"FRem","BinaryOperator")...),
	20 => InstructionDef{:Binary}((20,"Shl","BinaryOperator")...),
	21 => InstructionDef{:Binary}((21,"LShr","BinaryOperator")...),
	22 => InstructionDef{:Binary}((22,"AShr","BinaryOperator")...),
	23 => InstructionDef{:Binary}((23,"And","BinaryOperator")...),
	24 => InstructionDef{:Binary}((24,"Or","BinaryOperator")...),
	25 => InstructionDef{:Binary}((25,"Xor","BinaryOperator")...),
	26 => InstructionDef{:Memory}((26,"Alloca","AllocaInst")...),
	27 => InstructionDef{:Memory}((27,"Load","LoadInst")...),
	28 => InstructionDef{:Memory}((28,"Store","StoreInst")...),
	29 => InstructionDef{:Memory}((29,"GetElementPtr","GetElementPtrInst")...),
	30 => InstructionDef{:Memory}((30,"Fence","FenceInst")...),
	31 => InstructionDef{:Memory}((31,"AtomicCmpXchg","AtomicCmpXchgInst")...),
	32 => InstructionDef{:Memory}((32,"AtomicRMW","AtomicRMWInst")...),
	33 => InstructionDef{:Cast}((33,"Trunc","TruncInst")...),
	34 => InstructionDef{:Cast}((34,"ZExt","ZExtInst")...),
	35 => InstructionDef{:Cast}((35,"SExt","SExtInst")...),
	36 => InstructionDef{:Cast}((36,"FPToUI","FPToUIInst")...),
	37 => InstructionDef{:Cast}((37,"FPToSI","FPToSIInst")...),
	38 => InstructionDef{:Cast}((38,"UIToFP","UIToFPInst")...),
	39 => InstructionDef{:Cast}((39,"SIToFP","SIToFPInst")...),
	40 => InstructionDef{:Cast}((40,"FPTrunc","FPTruncInst")...),
	41 => InstructionDef{:Cast}((41,"FPExt","FPExtInst")...),
	42 => InstructionDef{:Cast}((42,"PtrToInt","PtrToIntInst")...),
	43 => InstructionDef{:Cast}((43,"IntToPtr","IntToPtrInst")...),
	44 => InstructionDef{:Cast}((44,"BitCast","BitCastInst")...),
	45 => InstructionDef{:Cast}((45,"AddrSpaceCast","AddrSpaceCastInst")...),
	46 => InstructionDef{:Other}((46,"ICmp","ICmpInst")...),
	47 => InstructionDef{:Other}((47,"FCmp","FCmpInst")...),
	48 => InstructionDef{:Other}((48,"PHI","PHINode")...),
	49 => InstructionDef{:Other}((49,"Call","CallInst")...),
	50 => InstructionDef{:Other}((50,"Select","SelectInst")...),
	51 => InstructionDef{:Other}((51,"UserOp1","Instruction")...),
	52 => InstructionDef{:Other}((52,"UserOp2","Instruction")...),
	53 => InstructionDef{:Other}((53,"VAArg","VAArgInst")...),
	54 => InstructionDef{:Other}((54,"ExtractElement","ExtractElementInst")...),
	55 => InstructionDef{:Other}((55,"InsertElement","InsertElementInst")...),
	56 => InstructionDef{:Other}((56,"ShuffleVector","ShuffleVectorInst")...),
	57 => InstructionDef{:Other}((57,"ExtractValue","ExtractValueInst")...),
	58 => InstructionDef{:Other}((58,"InsertValue","InsertValueInst")...),
	59 => InstructionDef{:Other}((59,"LandingPad","LandingPadInst")...),
]


text:
using CuckooFilters, Documenter, Test

@testset "CuckooFilter tests" begin
    @testset "Insertion tests" begin
        filter = CuckooFilter{AbstractVector{UInt8}}()
        @test b"hello, world!" ∉ filter
        @test b"goodbye!" ∉ filter

        @test insert!(filter, b"hello, world!")
        @test insert!(filter, b"goodbye!")

        @test b"hello, world!" ∈ filter
        @test b"hello, world" ∉ filter
        @test b"goodbye!" ∈ filter

        filter = CuckooFilter{Int64}()
        for ii = 1:256
            insert!(filter, ii)
        end

        @test all(x ∈ filter for x = 1:256)
    end

    @testset "Unique insertion tests" begin
        filter = CuckooFilter{AbstractVector{UInt8}}()
        @test insert_unique!(filter, b"hello, world!")
        @test b"hello, world!" ∈ filter

        # insert! returns whether or not the insertion was successful. If we're
        # performing unique insertions, then our second attempt at insertion
        # should return false.
        @test insert_unique!(filter, b"hello, world!") == false
    end

    @testset "Deletion tests" begin
        filter = CuckooFilter{AbstractVector{UInt8}}()
        insert!(filter, b"hello, world!")
        insert!(filter, b"goodbye!")
        delete!(filter, b"hello, world!")

        @test b"hello, world!" ∉ filter
        @test b"goodbye!" ∈ filter
    end
end

@testset "CuckooFilter with custom type" begin
    # Test operations on CuckooFilter using a custom type that supports the
    # fingerprint and fingerprint_type methods
    struct Foo
        x::Int64
        y::Int64
    end

    CuckooFilters.fingerprint(f::Foo) = UInt32((f.x ⊻ f.y) & 0xffff)
    CuckooFilters.fingerprint_type(::Type{Foo}) = UInt32

    @testset "Insertion tests" begin
        filter = CuckooFilter{Foo}()
        @test Foo(1, 2) ∉ filter
        @test Foo(2, 3) ∉ filter
        @test insert!(filter, Foo(1, 2))
        @test insert!(filter, Foo(2, 3))
        @test Foo(1, 2) ∈ filter
        @test Foo(2, 3) ∈ filter

        # Should not be able to insert another type into the CuckooFilter
        @test_throws MethodError insert!(filter, 5)
    end

    @testset "Deletion tests" begin
        filter = CuckooFilter{Foo}()
        @test Foo(-5, 18) ∉ filter
        @test delete!(filter, Foo(-5, 18)) == false
        @test insert!(filter, Foo(-5, 18))
        @test Foo(-5, 18) ∈ filter
        @test delete!(filter, Foo(-5, 18))
        @test Foo(-5, 18) ∉ filter
    end
end

@testset "CuckooFilter doctests" begin
    doctest(CuckooFilters)
end


text:
using Makie

tree_theme = Makie.Theme(
                Lines = (
                    color = :red,
                ),
                Axis = (
                    leftspinevisible = false,
                    rightspinevisible = false,
                    bottomspinevisible = false,
                    topspinevisible = false,
                    xticksvisible = false,
                    yticksvisible = false,
                    xlabelsvisible = false,
                    ylabelsvisible = false,
                    gridvisible = false,
                )
            )

@Makie.recipe(TreePlot, tree) do scene
    # default = tree_theme
    Makie.Attributes(;
        color = :black,
        labels = :all,
        labeloffset = 0,
        branchwidth = 1,
        usebrlength = true,
    )
end

Makie.plottype(::RTree) = TreePlot

function Makie.plot!(plt::TreePlot)
    tree = plt[:tree][]

    A, B = _tree_coords(tree, plt[:usebrlength][])

    x, y = A
    xline, yline = B

    branchwidth = plt[:branchwidth][]
    Makie.lines!(plt, xline, yline, linewidth=branchwidth, color = plt[:color][])
    labelmode = plt[:labels][]
    labeloffset = plt[:labeloffset][]
    if labelmode == :all
        pltlabels = getlabel.(tree.nodes)
        Makie.text!(plt, pltlabels, position = Point.(x, y), align = (:left, :center), space=:screen)
    end
    # hidedecorations!(plt[Axis])

    return plt
end

max_stars_repo_name: eascarrunz/RootedTrees

max_issues_repo_name: eascarrunz/RootedTrees

max_forks_repo_name: eascarrunz/RootedTrees


text:
function tensor_create{order, dim}(::Type{Tensor{order, dim}}, f)
    if order == 1
        ex = Expr(:tuple, [f(i) for i=1:dim]...)
    elseif order == 2
        ex = Expr(:tuple, [f(i,j) for i=1:dim, j=1:dim]...)
    elseif order == 4
        ex = Expr(:tuple, [f(i,j,k,l) for i=1:dim, j=1:dim, k = 1:dim, l = 1:dim]...)
    end
    return ex
end

function tensor_create{order, dim}(::Type{SymmetricTensor{order, dim}}, f)
    ex = Expr(:tuple)
    if order == 2
        for j in 1:dim, i in j:dim
            push!(ex.args, f(i, j))
        end
    elseif order == 4
        for l in 1:dim, k in l:dim, j in 1:dim, i in j:dim
            push!(ex.args, f(i, j, k, l))
        end
    end
    return ex
end

# create recursive muladd exp from two expression arrays
function make_muladd_exp(ex1, ex2)
    N = length(ex1)
    ex = Expr(:call)
    exn = Expr(:call, :*, ex1[1], ex2[1])

    if N == 1 # return only the multiplication
        return exn
    end

    for i in 2:N
        ex = Expr(:call, :muladd)
        push!(ex.args, ex1[i])
        push!(ex.args, ex2[i])
        push!(ex.args, exn)
        exn = ex
    end
    return ex
end


text:
@testset "hostbyteorder" begin
    if Base.ENDIAN_BOM == 0x04030201
        @test MRC.hostbyteorder() == MRC.LittleEndian
    else
        @test MRC.hostbyteorder() == MRC.BigEndian
    end
end

@testset "machstfrombyteorder" begin
    @test MRC.machstfrombyteorder() == [0x44, 0x44, 0x00, 0x00]
    @test MRC.machstfrombyteorder(MRC.LittleEndian) == [0x44, 0x44, 0x00, 0x00]
    @test MRC.machstfrombyteorder(MRC.BigEndian) == [0x11, 0x11, 0x00, 0x00]
end

@testset "byteorderfrommachst" begin
    @test MRC.byteorderfrommachst([0x44, 0x44, 0x00, 0x00]) == MRC.LittleEndian
    @test MRC.byteorderfrommachst([0x44, 0x41, 0x00, 0x00]) == MRC.LittleEndian
    @test MRC.byteorderfrommachst([0x11, 0x11, 0x00, 0x00]) == MRC.BigEndian
    @test (@test_logs (
        :warn,
        "Unrecognized machine stamp $([0x00, 0x00, 0x00, 0x00]). Assuming little endian.",
    ) MRC.byteorderfrommachst([0x00, 0x00, 0x00, 0x00])) == MRC.LittleEndian
end

@testset "bswaptoh" begin
    @test MRC.bswaptoh(MRC.LittleEndian) == MRC.maybeswap(ltoh)
    @test MRC.bswaptoh(MRC.MACHINE_STAMP_LITTLE) == MRC.maybeswap(ltoh)
    @test MRC.bswaptoh(MRC.MACHINE_STAMP_LITTLE_ALT) == MRC.maybeswap(ltoh)
    @test MRC.bswaptoh(MRC.BigEndian) == MRC.maybeswap(ntoh)
    @test MRC.bswaptoh(MRC.MACHINE_STAMP_BIG) == MRC.maybeswap(ntoh)
end

@testset "bswapfromh" begin
    @test MRC.bswapfromh(MRC.LittleEndian) == MRC.maybeswap(htol)
    @test MRC.bswapfromh(MRC.MACHINE_STAMP_LITTLE) == MRC.maybeswap(htol)
    @test MRC.bswapfromh(MRC.MACHINE_STAMP_LITTLE_ALT) == MRC.maybeswap(htol)
    @test MRC.bswapfromh(MRC.BigEndian) == MRC.maybeswap(hton)
    @test MRC.bswapfromh(MRC.MACHINE_STAMP_BIG) == MRC.maybeswap(hton)
end

@testset "maybeswap" begin
    @test MRC.maybeswap(bswap, Int32(5)) === bswap(Int32(5))
    @test MRC.maybeswap(bswap, "foo") === "foo"
    @test MRC.maybeswap(bswap, (Int32(5), Int32(6))) === (bswap(Int32(5)), bswap(Int32(6)))
    @test MRC.maybeswap(bswap, (0x01, 0x02)) === (0x01, 0x02)
    @test MRC.maybeswap(bswap)(Int32(5)) === bswap(Int32(5))
end

@testset "padtruncto!" begin
    x = [1, 2, 3]
    MRC.padtruncto!(x, 4)
    @test x == [1, 2, 3, 0]
    MRC.padtruncto!(x, 6; value=1)
    @test x == [1, 2, 3, 0, 1, 1]
    MRC.padtruncto!(x, 3)
    @test x == [1, 2, 3]
    MRC.padtruncto!(x, 5; value=1.0)
    @test x == [1, 2, 3, 1, 1]
end

@testset "compression: $type" for type in keys(MRC.COMPRESSIONS)
    spec = getproperty(MRC.COMPRESSIONS, type)
    @testset "checkmagic" begin
        io = IOBuffer()
        write(io, spec.magic)
        write(io, [0x01, 0x02, 0x03])
        seek(io, 0)
        type2 = MRC.checkmagic(io)
        @test type2 == type
        close(io)
    end

    @testset "checkextension" begin
        fn = "map.mrc$(spec.extension)"
        type2 = MRC.checkextension(fn)
        @test type2 == type
    end

    @testset "(de)compressstream" begin
        buf = IOBuffer()
        stream = MRC.compressstream(buf, type)
        write(stream, b"foo", TranscodingStreams.TOKEN_END)
        newbuf = IOBuffer(take!(buf))
        @test MRC.checkmagic(newbuf) == type
        newstream = MRC.decompressstream(newbuf, type)
        @test read(newstream) == b"foo"
    end
end


text:
using Knet
include("../infst.jl")
include("../mocommon.jl")

function oldbilstm(model, states, sequence; pdrop=(0, 0))
    total = 0.0
    count = 0
    atype = typeof(AutoGrad.getval(model[:fembed]))

    # forward lstm
    fhiddens = Array(Any, length(sequence))
    sf = copy(states)
    for i=1:length(sequence)-1
        input = convert(atype, sequence[i])
        x = input * model[:fembed]
        h = forward(model[:forw], sf, x)
        fhiddens[i+1] = copy(h)
    end
    fhiddens[1] = convert(atype, zeros(size(fhiddens[2])))

    # backward lstm
    bhiddens = Array(Any, length(sequence))
    sb = copy(states)
    for i=length(sequence):-1:2
        input = convert(atype, sequence[i])
        x = input * model[:bembed]
        bhiddens[i-1] = forward(model[:back], sb, x)
    end
    bhiddens[end] = convert(atype, zeros(size(bhiddens[2])))

    # concatenate layer
    for i=1:length(fhiddens)
        hf = dropout(fhiddens[i], pdrop[1])
        hb = dropout(bhiddens[i], pdrop[2])
        ypred = hcat(hf, hb) * model[:soft][1] .+ model[:soft][2]
        #ypred = hcat(fhiddens[i], bhiddens[i]) * model[:soft][1] .+ model[:soft][2]
        ynorm = logp(ypred, 2)
        ygold = convert(atype, sequence[i])
        count += size(ygold, 1)
        total += sum(ygold .* ynorm)
    end
    return - total / count
end

lg1 = grad(oldbilstm)

function bilstm(model, states, sequence; pdrop=(0, 0))
    total = 0.0
    count = 0
    atype = typeof(AutoGrad.getval(model[:fembed]))

    # forward lstm
    fhiddens = Array(Any, length(sequence))
    sf = copy(states)
    for i=1:length(sequence)-1
        x = model[:fembed][sequence[i], :]
        h = forward(model[:forw], sf, x)
        fhiddens[i+1] = copy(h)
    end
    fhiddens[1] = convert(atype, zeros(size(fhiddens[2])))

    # backward lstm
    bhiddens = Array(Any, length(sequence))
    sb = copy(states)
    for i=length(sequence):-1:2
        x = model[:bembed][sequence[i], :]
        bhiddens[i-1] = forward(model[:back], sb, x)
    end
    bhiddens[end] = convert(atype, zeros(size(bhiddens[2])))

    # concatenate layer
    for i=1:length(fhiddens)
        hf = dropout(fhiddens[i], pdrop[1])
        hb = dropout(bhiddens[i], pdrop[2])
        ypred = hcat(hf, hb) * model[:soft][1] .+ model[:soft][2]
        #ypred = hcat(fhiddens[i], bhiddens[i]) * model[:soft][1] .+ model[:soft][2]
        total += logprob(sequence[i], ypred)
        count += length(sequence[i])
    end
    return - total / count
end

lg2 = grad(bilstm)

function testmodels()
    atype = (gpu() >= 0 ? KnetArray{Float32} : Array{Float64})
    batchsize = 1
    hiddens = [128]
    embedding = 128
    f = open("/mnt/kufs/scratch/okirnap/lexsub/data/testdata.txt")
    global t_vocab = create_vocab("/mnt/kufs/scratch/okirnap/lexsub/data/testvocab")
    sdict = Dict{Int64, Array{Any, 1}}();
    readstream!(f, sdict, t_vocab)
    index_to_word = Array(AbstractString, length(t_vocab))
    for (k, v) in t_vocab; index_to_word[v] = k; end;

    ids = nextbatch(f, sdict, t_vocab, batchsize) # ids of the sequences
    println(ibuild_sentence(index_to_word, ids , 1))

    one2sen(d) = map(i-> (find(x->x==true, i)), d)

    # create minibatch for old model
    data = [ falses(batchsize, length(index_to_word)) for i=1:length(ids)]
    for i=1:length(ids)
        index = ids[i]
        data[i][index] = true
    end
    sen = map(x->index_to_word[x], one2sen(data))

    srand(12)
    m = initmodel(atype, hiddens, embedding, length(index_to_word))
    s = initstate(atype, hiddens, batchsize)
    @assert(bilstm(m, s, ids) == oldbilstm(m, s, data))
    info("forward test passed")
    g1 = lg1(m, s, data)
    g2 = lg2(m, s, ids)
    for (k,v) in g2
        @assert(g1[k] == v)
    end
    info("backward test passed")
end
!isinteractive() && testmodels()

max_stars_repo_name: kirnap/wordsubs

max_issues_repo_name: kirnap/wordsubs

max_forks_repo_name: kirnap/wordsubs


text:
using Memcache
using Base.Test

@test isa(Memcache, Module)

include("test.jl")


text:
#-------------------------------------------------------------------
#* EMSO Model Library (EML) Copyright (C) 2004 - 2007 ALSOC.
#*
#* This LIBRARY is free software; you can distribute it and/or modify
#* it under the therms of the ALSOC FREE LICENSE as available at
#* http://www.enq.ufrgs.br/alsoc.
#*
#* EMSO Copyright (C) 2004 - 2007 ALSOC, original code
#* from http://www.rps.eng.br Copyright (C) 2002-2004.
#* All rights reserved.
#*
#* EMSO is distributed under the therms of the ALSOC LICENSE as
#* available at http://www.enq.ufrgs.br/alsoc.
#*
#*----------------------------------------------------------------------
#* Author: Paula B. Staudt
#* $Id: tray.mso 522 2008-05-21 23:21:12Z arge $
#*--------------------------------------------------------------------
type tray
	tray()=begin
		PP=outers.PP
		NComp=outers.NComp
		new(
			DanaPlugin(Dict{Symbol,Any}(
				:Brief=>"External Physical Properties",
				:Type=>"PP"
			)),
			DanaInteger(Dict{Symbol,Any}(
				:Brief=>"Number of components"
			)),
			fill(molweight(Dict{Symbol,Any}(
				:Brief=>"Component Mol Weight",
				:Hidden=>true
			)),(NComp)),
			acceleration(Dict{Symbol,Any}(
				:Brief=>"Gravity Acceleration",
				:Default=>9.81,
				:Hidden=>true
			)),
			flow_mol(Dict{Symbol,Any}(
				:Brief=>"Stream Flow closed",
				:Default=>0,
				:Hidden=>true
			)),
			flow_mol(Dict{Symbol,Any}(
				:Brief=>"Low stream Flow",
				:Default=>1E-6,
				:Hidden=>true
			)),
			constant(Dict{Symbol,Any}(
				:Brief=>"Pi Number",
				:Default=>3.14159265,
				:Symbol=>"\\pi",
				:Hidden=>true
			)),
			length(Dict{Symbol,Any}(
				:Brief=>"Tray Diameter",
				:Default=>1.600
			)),
			length(Dict{Symbol,Any}(
				:Brief=>"Tray Spacing",
				:Default=>0.600
			)),
			fraction(Dict{Symbol,Any}(
				:Brief=>"Fraction of the active area that is occupied by the holes with respect to the total tray area",
				:Default=>0.10
			)),
			fraction(Dict{Symbol,Any}(
				:Brief=>"Fraction of the downcomer area with respect to the total tray area",
				:Default=>0.20
			)),
			length(Dict{Symbol,Any}(
				:Brief=>"Weir length",
				:Default=>1
			)),
			length(Dict{Symbol,Any}(
				:Brief=>"Weir height",
				:Default=>0.05
			)),
			positive(Dict{Symbol,Any}(
				:Brief=>"Number of liquid passes in the tray",
				:Lower=>1,
				:Default=>1
			)),
			heat_rate(Dict{Symbol,Any}(
				:Brief=>"Rate of heat supply",
				:Default=>0
			)),
			DanaReal(Dict{Symbol,Any}(
				:Brief=>"Aeration fraction",
				:Default=>1
			)),
			DanaReal(Dict{Symbol,Any}(
				:Brief=>"Dry pressure drop coefficient",
				:Default=>0.60
			)),
			DanaReal(Dict{Symbol,Any}(
				:Brief=>"Murphree efficiency for All Trays",
				:Lower=>0.01,
				:Upper=>1
			)),
			area(Dict{Symbol,Any}(
				:Brief=>"Plate area = Atray - Adowncomer",
				:Protected=>true
			)),
			volume(Dict{Symbol,Any}(
				:Brief=>"Total Volume of the tray",
				:Protected=>true
			)),
			area(Dict{Symbol,Any}(
				:Brief=>"Total holes area",
				:Protected=>true
			)),
			DanaReal(Dict{Symbol,Any}(
				:Brief=>"Feeherys correlation coefficient",
				:Unit=>"1/m^4",
				:Default=>1,
				:Hidden=>true
			)),
			DanaReal(Dict{Symbol,Any}(
				:Brief=>"Elgues correlation coefficient",
				:Unit=>"kg/m/mol^2",
				:Default=>1,
				:Hidden=>true
			)),
			DanaReal(Dict{Symbol,Any}(
				:Brief=>"Olsens correlation coefficient",
				:Default=>1,
				:Hidden=>true
			)),
			DanaSwitcher(Dict{Symbol,Any}(
				:Brief=>"Flag for Vapour Flow condition",
				:Valid=>["on", "off"],
				:Default=>"off",
				:Hidden=>true
			)),
			DanaSwitcher(Dict{Symbol,Any}(
				:Brief=>"Flag for Liquid Flow condition",
				:Valid=>["on", "off"],
				:Default=>"off",
				:Hidden=>true
			)),
			stream(Dict{Symbol,Any}(
				:Brief=>"Feed stream",
				:Hidden=>true,
				:PosX=>0,
				:PosY=>0.4932,
				:Symbol=>"_{in}"
			)),
			liquid_stream(Dict{Symbol,Any}(
				:Brief=>"liquid Sidestream",
				:Hidden=>true,
				:Symbol=>"_{outL}"
			)),
			vapour_stream(Dict{Symbol,Any}(
				:Brief=>"vapour Sidestream",
				:Hidden=>true,
				:Symbol=>"_{outV}"
			)),
			stream(Dict{Symbol,Any}(
				:Brief=>"Inlet liquid stream",
				:PosX=>0.5195,
				:PosY=>0,
				:Symbol=>"_{inL}"
			)),
			stream(Dict{Symbol,Any}(
				:Brief=>"Inlet vapour stream",
				:PosX=>0.4994,
				:PosY=>1,
				:Symbol=>"_{inV}"
			)),
			liquid_stream(Dict{Symbol,Any}(
				:Brief=>"Outlet liquid stream",
				:PosX=>0.8277,
				:PosY=>1,
				:Symbol=>"_{outL}"
			)),
			vapour_stream(Dict{Symbol,Any}(
				:Brief=>"Outlet vapour stream",
				:PosX=>0.8043,
				:PosY=>0,
				:Symbol=>"_{outV}"
			)),
			positive(Dict{Symbol,Any}(
				:Brief=>"Flag for Liquid Flow Model",
				:Lower=>0,
				:Default=>1 ,
				:Hidden=>true
			)),
			positive(Dict{Symbol,Any}(
				:Brief=>"Flag for Vapour Flow Model",
				:Lower=>0,
				:Default=>1 ,
				:Hidden=>true
			)),
			fill(mol(Dict{Symbol,Any}(
				:Brief=>"Molar Holdup in the tray"
			)),(NComp)),
			mol(Dict{Symbol,Any}(
				:Brief=>"Molar liquid holdup"
			)),
			mol(Dict{Symbol,Any}(
				:Brief=>"Molar vapour holdup"
			)),
			energy(Dict{Symbol,Any}(
				:Brief=>"Total Energy Holdup on tray"
			)),
			volume_mol(Dict{Symbol,Any}(
				:Brief=>"Liquid Molar Volume"
			)),
			volume_mol(Dict{Symbol,Any}(
				:Brief=>"Vapour Molar volume"
			)),
			length(Dict{Symbol,Any}(
				:Brief=>"Height of clear liquid on plate"
			)),
			fill(fraction()),
			dens_mass(Dict{Symbol,Any}(
				:Brief=>"Mass Density of liquid phase"
			)),
			dens_mass(Dict{Symbol,Any}(
				:Brief=>"Mass Density of vapour phase"
			)),
			[
				:(OutletLiquid.F*vL = 1.84*"1/s"*WeirLength_*((Level-(AerationFraction_*WeirHeight_))/(AerationFraction_))^2),
				:(OutletLiquid.F*vL = 1.84*"m^0.5/s"*WeirLength_*((Level-(AerationFraction_*WeirHeight_))/(AerationFraction_))^1.5),
				:(OutletLiquid.F / "mol/s"= WeirLength_*TrayLiquidPasses_*rhoL/sum(Mw*OutletVapour.z)/(0.665*OlsenCoeff)^1.5 * ((ML*sum(Mw*OutletLiquid.z)/rhoL/PlateArea_)-WeirHeight_)^1.5 * "m^0.5/mol"),
				:(OutletLiquid.F = WeirLength_*rhoL/sum(Mw*OutletLiquid.z) * ((Level-WeirHeight_)/750/"mm")^1.5 * "m^2/s"),
				:(OutletLiquid.F = 2/3*rhoL/sum(Mw*OutletLiquid.z)*WeirLength_*(ML*sum(Mw*OutletLiquid.z)/(PlateArea_*1.3)/rhoL)^1.5*sqrt(2*Gconst/ (2*(1 - 0.3593/"Pa^0.0888545"*abs(OutletVapour.F*sum(Mw*OutletVapour.z)/(PlateArea_*1.3)/sqrt(rhoV))^0.177709)-1))),
				:(OutletLiquid.F = zero_flow),
				:(InletVapour.F*vV = sqrt((InletVapour.P - OutletVapour.P)/(rhoV*DryPdropCoeff_))*HolesArea_),
				:(InletVapour.F = rhoV/PlateArea_/FeeheryCoeff/sum(Mw*OutletVapour.z) * sqrt(((InletVapour.P - OutletVapour.P)-(rhoV*Gconst*ML*vL/PlateArea_))/rhoV)),
				:(InletVapour.F^1.08 * 0.0013 * "kg/m/mol^1.08/s^0.92*1e5" = (InletVapour.P - OutletVapour.P)*1e5 - (AerationFraction_*sum(M*Mw)/(PlateArea_*1.3)*Gconst*1e5) * (rhoV*HolesArea_/sum(Mw*OutletVapour.z))^1.08 * "m^1.08/mol^1.08"),
				:(InletVapour.F * vV = PlateArea_ * sqrt(((InletVapour.P - OutletVapour.P)-rhoL*Gconst*Level)/rhoV)),
				:(InletVapour.F * vV = PlateArea_ * sqrt(((InletVapour.P - OutletVapour.P)-rhoL*Gconst*Level)/rhoV*DryPdropCoeff_)),
				:(InletVapour.F = sqrt((InletVapour.P - OutletVapour.P)/ElgueCoeff)),
				:(InletVapour.F = zero_flow),
				:(OutletVapour.z = MurphreeEff_ * (yideal - InletVapour.z) + InletVapour.z),
				:(diff(E) = ( Inlet.F*Inlet.h + InletLiquid.F*InletLiquid.h + InletVapour.F*InletVapour.h- OutletLiquid.F*OutletLiquid.h - OutletVapour.F*OutletVapour.h - VapourSideStream.F*VapourSideStream.h - LiquidSideStream.F*LiquidSideStream.h + HeatSupply_ )),
				:(E = ML*OutletLiquid.h + MV*OutletVapour.h - OutletLiquid.P*TrayVolume_),
				:(TrayVolume_ = ML* vL + MV*vV),
				:(Level = ML*vL/PlateArea_),
				:(diff(M)=Inlet.F*Inlet.z + InletLiquid.F*InletLiquid.z + InletVapour.F*InletVapour.z- OutletLiquid.F*OutletLiquid.z - OutletVapour.F*OutletVapour.z- LiquidSideStream.F*LiquidSideStream.z-VapourSideStream.F*VapourSideStream.z),
				:(M = ML*OutletLiquid.z + MV*OutletVapour.z),
				:(sum(OutletLiquid.z)= 1.0),
				:(sum(OutletLiquid.z)= sum(OutletVapour.z)),
				:(vL = PP.LiquidVolume(OutletLiquid.T, OutletLiquid.P, OutletLiquid.z)),
				:(vV = PP.VapourVolume(OutletVapour.T, OutletVapour.P, OutletVapour.z)),
				:(rhoL = PP.LiquidDensity(OutletLiquid.T, OutletLiquid.P, OutletLiquid.z)),
				:(rhoV = PP.VapourDensity(InletVapour.T, InletVapour.P, InletVapour.z)),
				:(PP.LiquidFugacityCoefficient(OutletLiquid.T, OutletLiquid.P, OutletLiquid.z)*OutletLiquid.z = PP.VapourFugacityCoefficient(OutletVapour.T, OutletVapour.P, yideal)*yideal),
				:(OutletVapour.T = OutletLiquid.T),
				:(OutletVapour.P = OutletLiquid.P),
				:(OutletVapour.T = VapourSideStream.T),
				:(OutletLiquid.T = LiquidSideStream.T),
				:(OutletVapour.P= VapourSideStream.P),
				:(OutletLiquid.P = LiquidSideStream.P),
				:(OutletLiquid.z= LiquidSideStream.z),
				:(OutletVapour.z= VapourSideStream.z),
			],
			[
				"Francis Equation","Wang_Fl","Olsen","Feehery_Fl","Roffel_Fl","Low level","Reepmeyer","Feehery_Fv","Roffel_Fv","Klingberg","Wang_Fv","Elgue","","Murphree Efficiency","Energy Balance","Energy Holdup","Geometry Constraint","Level of clear liquid over the weir","Component Molar Balance","Molar Holdup","Mol fraction normalisation","Mol fraction constraint","Liquid Volume","Vapour Volume","Liquid Density","Vapour Density","Chemical Equilibrium","Thermal Equilibrium","Mechanical Equilibrium","Thermal Equilibrium Vapour Side Stream","Thermal Equilibrium Liquid Side Stream","Mechanical Equilibrium Vapour Side Stream","Mechanical Equilibrium Liquid Side Stream","Composition Liquid Side Stream","Composition Vapour Side Stream",
			],
			[:PP,:NComp,:Mw,:Gconst,:zero_flow,:low_flow,:Pi,:TrayDiameter_,:TraySpacing_,:Fraction_HoleArea_,:Fraction_DowncomerArea_,:WeirLength_,:WeirHeight_,:TrayLiquidPasses_,:HeatSupply_,:AerationFraction_,:DryPdropCoeff_,:MurphreeEff_,:PlateArea_,:TrayVolume_,:HolesArea_,:FeeheryCoeff,:ElgueCoeff,:OlsenCoeff,:VapourFlow,:LiquidFlow,],
			[:Inlet,:LiquidSideStream,:VapourSideStream,:InletLiquid,:InletVapour,:OutletLiquid,:OutletVapour,:LFlowModel,:VFlowModel,:M,:ML,:MV,:E,:vL,:vV,:Level,:yideal,:rhoL,:rhoV,]
		)
	end
	PP::DanaPlugin
	NComp::DanaInteger
	Mw::Array{molweight}
	Gconst::acceleration
	zero_flow::flow_mol
	low_flow::flow_mol
	Pi::constant
	TrayDiameter_::length
	TraySpacing_::length
	Fraction_HoleArea_::fraction
	Fraction_DowncomerArea_::fraction
	WeirLength_::length
	WeirHeight_::length
	TrayLiquidPasses_::positive
	HeatSupply_::heat_rate
	AerationFraction_::DanaReal
	DryPdropCoeff_::DanaReal
	MurphreeEff_::DanaReal
	PlateArea_::area
	TrayVolume_::volume
	HolesArea_::area
	FeeheryCoeff::DanaReal
	ElgueCoeff::DanaReal
	OlsenCoeff::DanaReal
	VapourFlow::DanaSwitcher
	LiquidFlow::DanaSwitcher
	Inlet::stream
	LiquidSideStream::liquid_stream
	VapourSideStream::vapour_stream
	InletLiquid::stream
	InletVapour::stream
	OutletLiquid::liquid_stream
	OutletVapour::vapour_stream
	LFlowModel::positive
	VFlowModel::positive
	M::Array{mol}
	ML::mol
	MV::mol
	E::energy
	vL::volume_mol
	vV::volume_mol
	Level::length
	yideal::Array{fraction}
	rhoL::dens_mass
	rhoV::dens_mass
	equations::Array{Expr,1}
	equationNames::Array{String,1}
	parameters::Array{Symbol,1}
	variables::Array{Symbol,1}
	attributes::Dict{Symbol,Any}
end
export tray
function set(in::tray)
	Mw = PP.MolecularWeight()
	 zero_flow = 0 * "kmol/h"
	 low_flow = 1E-6 * "kmol/h"
	 PlateArea_ = 0.25*Pi*(TrayDiameter_^2)*(1-Fraction_DowncomerArea_)
	 TrayVolume_ = 0.25*Pi*(TrayDiameter_^2)*TraySpacing_
	 HolesArea_ = 0.25*Pi*(TrayDiameter_^2)*Fraction_HoleArea_

end
function setEquationFlow(in::tray)
	# LiquidFlow and VapourFlow equations need to be linerized to avoid indetermination !
	let switch=LiquidFlow
		if LiquidFlow==Level < (AerationFraction_ *WeirHeight_)
			set(switch,"off")
		end
		if LiquidFlow==Level > (AerationFraction_ * WeirHeight_)
			set(switch,"on")
		end
		if switch=="on"
			if LFlowModel == 1
				addEquation(1)
			else
				if LFlowModel == 2
					addEquation(2)
				else
					if LFlowModel == 3
						addEquation(3)
					else
						if LFlowModel == 4
							addEquation(4)
						else
							addEquation(5)
							#/'(kg/m)^0.0888545/s^0.177709';

						end
					end
				end
			end
		elseif switch=="off"
			addEquation(6)
		end
	end
	let switch=VapourFlow
		if VapourFlow==InletVapour.F < low_flow
			set(switch,"off")
		end
		if VapourFlow==InletVapour.P > OutletVapour.P
			set(switch,"on")
		end
		if switch=="on"
			if VFlowModel == 1
				addEquation(7)
			else
				if VFlowModel == 2
					addEquation(8)
				else
					if VFlowModel == 3
						addEquation(9)
					else
						if VFlowModel == 4
							addEquation(10)
						else
							if VFlowModel == 5
								addEquation(11)
							else
								addEquation(12)
							end
						end
					end
				end
			end
		elseif switch=="off"
			addEquation(13)
		end
	end
	addEquation(14)
	addEquation(15)
	addEquation(16)
	addEquation(17)
	addEquation(18)
	addEquation(19)
	addEquation(20)
	addEquation(21)
	addEquation(22)
	addEquation(23)
	addEquation(24)
	addEquation(25)
	addEquation(26)
	addEquation(27)
	addEquation(28)
	addEquation(29)
	addEquation(30)
	addEquation(31)
	addEquation(32)
	addEquation(33)
	addEquation(34)
	addEquation(35)
end
function atributes(in::tray,_::Dict{Symbol,Any})
	fields::Dict{Symbol,Any}=Dict{Symbol,Any}()
	fields[:Pallete]=false
	fields[:Icon]="icon/Tray"
	fields[:Brief]="Complete model of a column tray."
	fields[:Info]="== Assumptions ==
* both phases (liquid and vapour) exists all the time;
* thermodymanic equilibrium with Murphree plate efficiency;
* no entrainment of liquid or vapour phase;
* no weeping;
* the dymanics in the downcomer are neglected.

== Options ==
You can choose the equation for the liquid outlet flow and the vapour
inlet flow calculation through the VapourFlowModel and LiquidFlowModel
switchers.

== References ==
* ELGUE, S.; PRAT, L.; CABASSUD, M.; LANN, J. L.; C虅ZERAC, J. Dynamic models for start-up operations of batch distillation columns with experimental validation. Computers and Chemical Engineering, v. 28, p. 2735-2747, 2004.
* FEEHERY, W. F. Dynamic Optimization with Path Constraints. Tese (Doutorado) - Massachusetts Institute of Technology, June 1998.
* KLINGBERG, A. Modeling and Optimization of Batch Distillation. Disserta袯⭠o (Mestrado) - Department of Automatic Control, Lund Institute of Technology, Lund, Sweden, fev. 2000.
* OLSEN, I.; ENDRESTOL, G. O.; SIRA, T. A rigorous and efficient distillation column model for engineering and training simulators. Computers and Chemical Engineering,v. 21, n. Suppl, p. S193-S198, 1997.
* REEPMEYER, F.; REPKE, J.-U.; WOZNY, G. Analysis of the start-up process for reactive distillation. Chemical Engineering Technology, v. 26, p. 81-86, 2003.
* ROFFEL, B.; BETLEM, B.; RUIJTER, J. de. First principles dynamic modeling and multivariable control of a cryogenic distillation column process. Computers and Chemical Engineering, v. 24, p. 111-123, 2000.
* WANG, L.; LI, P.; WOZNY, G.; WANG, S. A start-up model for simulation of batch distillation starting from a cold state. Computers and Chemical Engineering, v. 27, p.1485-1497, 2003.
"
	drive!(fields,_)
	return fields
end
tray(_::Dict{Symbol,Any})=begin
	newModel=tray()
	newModel.attributes=atributes(newModel,_)
	newModel
end


text:
using BenchmarkTools, MINLPLib

const SUITE = BenchmarkGroup()


text:
time_period_conversion(time_period::Union{Dates.TimePeriod, Dates.DatePeriod}) = convert(Dates.Millisecond, time_period)
time_period_conversion(time_periods::Dict{String, <:Dates.Period}) = convert(Dict{String, Dates.Millisecond}, time_periods)