text:
# the main loop for the trampoline is parameterised over the Config
# type.  this lets us implement different semantics by simply changing
# the provided Config subtype (and implementing the appropriate
# dispatch() functions, etc).

# this works for most configs
parent(k::Config) = k.stack[end][1]

# evaluation without a cache (memoization)

type NoCache{S,I}<:Config{S,I}
    source::S
    @compat stack::Vector{Tuple{Matcher, State}}
    @compat NoCache(source; kargs...) = new(source, Vector{Tuple{Matcher,State}}())
end

function dispatch(k::NoCache, e::Execute)
    push!(k.stack, (e.parent, e.parent_state))
    try
        execute(k, e.child, e.child_state, e.iter)
    catch x
        isa(x, FailureException) ? FAILURE : rethrow()
    end
end

function dispatch(k::NoCache, s::Success)
    (parent, parent_state) = pop!(k.stack)
    try
        success(k, parent, parent_state, s.child_state, s.iter, s.result)
    catch x
        isa(x, FailureException) ? FAILURE : rethrow()
    end
end

function dispatch(k::NoCache, f::Failure)
    (parent, parent_state) = pop!(k.stack)
    try
        failure(k, parent, parent_state)
    catch x
        isa(x, FailureException) ? FAILURE : rethrow()
    end
end

# evaluation with a complete cache (all intermediate results memoized)

@compat typealias Key{I} Tuple{Matcher,State,I}

type Cache{S,I}<:Config{S,I}
    source::S
    @compat stack::Vector{Tuple{Matcher,State,Key{I}}}
    cache::Dict{Key{I},Message}
    @compat Cache(source; kargs...) = new(source, Vector{Tuple{Matcher,State,Key{I}}}(), Dict{Key{I},Message}())
end

function dispatch(k::Cache, e::Execute)
    key = (e.child, e.child_state, e.iter)
    push!(k.stack, (e.parent, e.parent_state, key))
    if haskey(k.cache, key)
        k.cache[key]
    else
        try
            execute(k, e.child, e.child_state, e.iter)
        catch x
            isa(x, FailureException) ? FAILURE : rethrow()
        end
    end
end

function dispatch(k::Cache, s::Success)
    parent, parent_state, key = pop!(k.stack)
    try
        k.cache[key] = s
    catch x
        isa(x, CacheException) ? nothing : rethrow()
    end
    try
        success(k, parent, parent_state, s.child_state, s.iter, s.result)
    catch x
        isa(x, FailureException) ? FAILURE : rethrow()
    end
end

function dispatch(k::Cache, f::Failure)
    parent, parent_state, key = pop!(k.stack)
    try
        k.cache[key] = f
    catch x
        isa(x, CacheException) ? nothing : rethrow()
    end
    try
        failure(k, parent, parent_state)
    catch x
        isa(x, FailureException) ? FAILURE : rethrow()
    end
end

# TODO - some kind of MRU cache?

# a dummy matcher used by the parser

type Root<:Delegate
    name::Symbol
    Root() = new(:Root)
end

immutable RootState<:DelegateState
    state::State
end

success(k::Config, m::Root, s::State, t::State, i, r::Value) = Success(RootState(t), i, r)
failure(k::Config, m::Root, s::State) = FAILURE

# the core loop that drives the parser, calling the appropriate dispatch
# functions (above) depending on which Config was used.
# to modify the behaviour you can create a new Config subtype and then
# add your own dispatch functions.

function producer(k::Config, m::Matcher; debug=false)

    root = Root()
    msg::Message = Execute(root, CLEAN, m, CLEAN, start(k.source))

    try

        while true
            msg = dispatch(k, msg)
            if isempty(k.stack)
                if isa(msg, Execute)
                    error("Unexpected execute message")
                elseif isa(msg, Success)
                    produce(msg.result)
                    # my head hurts
                    msg = Execute(root, CLEAN, m, msg.child_state.state, start(k.source))
                else
                    break
                end
            end
        end

    catch x
        if (debug)
            println("debug was set, so showing error from inside task")
            println(x)
            Base.show_backtrace(STDOUT, catch_backtrace())
        end
        throw(x)
    end

end

# helper functions to generate the parsers from the above

# these assume that any config construct takes a single source argument
# plus optional keyword args

function make{S}(config, source::S, matcher; debug=false, kargs...)
    I = typeof(start(source))
    k = config{S,I}(source; debug=debug, kargs...)
    (k, Task(() -> producer(k, matcher; debug=debug)))
end

function make_all(config; kargs_make...)
    function multiple_results(source, matcher::Matcher; kargs_parse...)
        kargs = vcat(kargs_make, kargs_parse)
        make(config, source, matcher; kargs...)[2]
    end
end

function once(task)
    result = consume(task)
    if task.state == :done
        throw(ParserException("cannot parse"))
    else
        return result
    end
end

function make_one(config; kargs_make...)
    function single_result(source, matcher::Matcher; kargs_parse...)
        kargs = vcat(kargs_make, kargs_parse)
        once(make(config, source, matcher; kargs...)[2])
    end
end

# the default parsers

# we use the cache only when matching multiple results but obviously motivated
# users are free to construct their own...

parse_all_cache = make_all(Cache)
parse_all_nocache = make_all(NoCache)
parse_all = parse_all_cache

parse_one_cache = make_one(Cache)
parse_one_nocache = make_one(NoCache)
parse_one = parse_one_nocache

parse_lines(source, matcher; kargs...) = parse_one(LineSource(source), matcher; kargs...)
parse_lines_cache(source, matcher; kargs...) = parse_one_cache(LineSource(source), matcher; kargs...)


text:
# Find the <k> nearest neighbour documents in <corpus> for <doc>.
function find_k_nearest_neighbours(dist::AbstractDocumentDistance, corpus, doc; k = 10)
    @assert k >= 1
    withdistances = map(cdoc -> (cdoc, evaluate(dist, cdoc, doc)), corpus)
    sort!(withdistances, by = t -> t[2])
    return withdistances[1:min(k, length(withdistances))]
end

ext: jl


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

# RUN: julia --startup-file=no %s %t && llvm-link -S %t/* -o %t/module.ll
# RUN: cat %t/module.ll | FileCheck %s

## Notes:
# This script uses the `emit` function (defined llvmpasses.jl) to emit either
# optimized or unoptimized LLVM IR. Each function is emitted individually and
# `llvm-link` is used to create a single module that can be passed to opt.
# The order in which files are emitted and linked is important since `lit` will
# process the test cases in order.

include(joinpath("..", "testhelpers", "llvmpasses.jl"))

import Base.FastMath

# CHECK: call fast float @llvm.sqrt.f32(float %0)
emit(FastMath.sqrt_fast, Float32)

# Float16 operations should be performed as Float32, unless @fastmath is specified
# TODO: this is not true for platforms that natively support Float16

foo(x::T,y::T) where T = x-y == zero(T)
# LOWER: fsub half %0, %1
# FINAL: %2 = fpext half %0 to float
# FINAL: %3 = fpext half %1 to float
# FINAL: fsub half %2, %3
emit(foo, Float16, Float16)

@fastmath foo(x::T,y::T) where T = x-y == zero(T)
# LOWER: fsub fast half %0, %1
# FINAL: fsub fast half %0, %1
emit(foo, Float16, Float16)

max_stars_repo_name: vanillajonathan/julia

max_issues_repo_name: vanillajonathan/julia

max_forks_repo_name: vanillajonathan/julia


text:
let S = ZZ{BigInt}, P = S[:x], x = P([0, 1])

@testset "construction of quotient ring" begin
    ideal = x^2 + 1
    @test P / ideal <: Quotient{P}
    @test P / (x^3+1) <: Quotient{P}
    Q = P / ideal
    Qp = P / (x^3+1)

    @test basetype(Q) == P
    @test depth(Q) == 3
    @test P / ideal != nothing
    p = 4x + 5
    @test Qp(p) == p
    @test Qp(1) == 1
    @test Qp(ZZ(1)) == 1
    q = Q(p)
    @test q + q == Q(8x + 10)
    @test q - q == 0
    @test q * q == p * p
    @test 3q == Q(12x + 15)
    pp = (4im + 5)^10
    @test q^10 == Q(P([pp.re, pp.im]))

    @test !iszero(q)
    @test iszero(zero(Q))
    @test !isone(q)
    @test isone(one(Q))
    @test isunit(Q(x))
    @test !isunit(Q(x+1))
    @test inv(Q(x)) == -Q(x)
    @test Q(x) * 5 == Q(5x)

    @test hash(q) != 0
    @test sprint(show, q) isa String

    @test value(q + 3 * ideal) == p
end

@testset "quotient of polynomial over GF" begin
    G = GF(2, 2)
    GX = G[:x]
    x = monom(GX)
    p = x^3 + 1
    Q = GX / p
    @test order(Q) == order(G)^3
end

@testset "show" begin
    P = ZZ{Int}[:x]
    x = monom(P)
    Q = P / (x^2 + 17)
    @test sprint(show, Q(1)) == "1"
    @test sprint(show, Q(x)) == "\u23b7-17"
    @test sprint(show, Q(-2x + 10)) == "-2*\u23b7-17 + 10"
    @test sprint(show, -Q(2x + 10)) == "-2*\u23b7-17 - 10"
    R = P / (x^3 + 1)
    @test sprint(show, R(x^2 + x + 1)) == "x^2 + x + 1 mod(x^3 + 1)"
end

end


text:
using Plots
gr()

module SurfaceHopping

export diabatic_potential,H_diabatic,diagonalise,adiabatic_potential,H_adiabatic
export diabatic_state, overlap_phi, diabatic_derivative, adiabatic_states, NACV, a_mn, g_mn
export force, acceleration, classical_propagation, new_a_mn
export riemann

using Distributions
using Calculus
using OrdinaryDiffEq

"""
Numeric integration, taken from http://mth229.github.io/integration.html
"""
function riemann(f::Function, a::Real, b::Real, n::Int; method="right")
    meth(f,l,r) = (1/6) * (f(l) + 4*(f((l+r)/2)) + f(r)) * (r-l)

    xs = a + (0:n) * (b-a)/n
    first = xs[1:end-1]; last = xs[2:end]; join = [first ; last];
    len = length(join)/2
    elements = reshape(join,(Int(len),2))
    pair = [elements[i,:] for i in 1:Int(len)]
    as = [meth(f,l,r) for (l,r) in pair]#meth(f, l, r)
    sum(as)

end

include("Propagators.jl")
include("States.jl")
include("PotentialEnergySurfaces.jl")

end


text:
"""
This example learns to classify hand-written digits from the
[MNIST](http://yann.lecun.com/exdb/mnist) dataset.  There are 60000
training and 10000 test examples. Each input x consists of 784 pixels
representing a 28x28 image.  The pixel values are normalized to
[0,1]. Each output y is converted to a ten-dimensional one-hot vector
(a vector that has a single non-zero component) indicating the correct
class (0-9) for a given image.  10 is used to represent 0.
You can run the demo using `julia mnist-cnn.jl` on the command line or
`julia> include(KnetLayers.dir("examples/mnist-cnn.jl"))` at the Julia prompt.
"""

using Knet, KnetLayers
import Knet: Data
#Data
include(Knet.dir("data","mnist.jl"))
dtrn,dtst = mnistdata(); # dtrn and dtst = [ (x1,y1), (x2,y2), ... ] where xi,yi are

#Model
struct CNN; layers::Chain; end
(m::CNN)(x) = m.layers(x)
(m::CNN)(x,y) = nll(m(x),y)
(m::CNN)(d::Data) = mean(m(x,y) for (x,y) in d)

model = CNN(Chain(Conv(height=5,width=5,inout=1=>12,activation=ReLU()),
                  Pool(window=2, stride=2),
                  Conv(height=5,width=5,inout=12=>12,activation=Tanh()),
                  Pool(window=2, stride=2),
                  x->mat(x),
                  Dense(input=4*4*12,output=13,activation=Tanh()),
                  Linear(input=13,output=10)))
#Train
EPOCH=10
progress!(sgd(model,repeat(dtrn,EPOCH)))

#Test
@show 100accuracy(model, dtst)


text:
using LinearAlgebra

## Test standard basis polynomials with (nearly) the same tests

#   compare upto trailing  zeros
function  upto_tz(as, bs)
    n,m = findlast.(!iszero, (as,bs))
    n == m || return false
    n == nothing &&  return true
    for i in 1:n
        as[i] != bs[i] && return false
    end
    true
end

upto_z(as, bs) = upto_tz(filter(!iszero,as), filter(!iszero,bs))

# compare upto trailing zeros infix operator
==ᵗ⁰(a,b) = upto_tz(a,b)
==ᵗᶻ(a,b) = upto_z(a,b)

Ps = (ImmutablePolynomial, Polynomial, SparsePolynomial, LaurentPolynomial)

@testset "Construction" for coeff in [
    Int64[1, 1, 1, 1],
    Float32[1, -4, 2],
    ComplexF64[1 - 1im, 2 + 3im],
    [3 // 4, -2 // 1, 1 // 1]
]

    for P in Ps
        p = P(coeff)
        @test coeffs(p) ==ᵗ⁰ coeff
        @test degree(p) == length(coeff) - 1
        @test p.var == :x
        P == Polynomial && @test length(p) == length(coeff)
        P == Polynomial && @test size(p) == size(coeff)
        P == Polynomial && @test size(p, 1) == size(coeff, 1)
        P == Polynomial && @test typeof(p).parameters[1] == eltype(coeff)
        @test eltype(p) == eltype(coeff)
        @test all([-200, -0.3, 1, 48.2] .∈ domain(p))
    end
end

@testset "Mapdomain" begin
    for P in Ps
        x = -30:20
        mx = mapdomain(P, x)
        @test mx == x

        x = 0.5:0.01:0.6
        mx = mapdomain(P, x)
        @test mx == x
    end
end

@testset "Other Construction" begin
    for P in (ImmutablePolynomial{10}, Polynomial)

        # Leading 0s
        p = P([1, 2, 0, 0])
        @test coeffs(p) ==ᵗ⁰ [1, 2]
        P == Polynomial && @test length(p) == 2

        # different type
        p = P{Float64}(ones(Int32, 4))
        @test coeffs(p) ==ᵗ⁰ ones(Float64, 4)

        p = P(30)
        @test coeffs(p) ==ᵗ⁰ [30]

        p = zero(P{Int})
        @test coeffs(p) ==ᵗ⁰ [0]

        p = one(P{Int})
        @test coeffs(p) ==ᵗ⁰ [1]

        pNULL = P(Int[])
        @test iszero(pNULL)
        @test degree(pNULL) == -1

        p0 = P([0])
        @test iszero(p0)
        @test degree(p0) == -1

        # variable(), P() to generate `x` in given basis
        @test degree(variable(P)) == 1
        @test variable(P)(1) == 1
        @test degree(P()) == 1
        @test P()(1) == 1
        @test variable(P, :y) == P(:y)

        # test degree, isconstant
        @test degree(zero(P)) ==  -1
        @test degree(one(P)) == 0
        @test degree(P(1))  ==  0
        @test degree(P([1]))  ==  0
        @test degree(P(:x)) ==  1
        @test degree(variable(P)) == 1
        @test degree(Polynomials.basis(P,5)) == 5
        @test Polynomials.isconstant(P(1))
        @test !Polynomials.isconstant(variable(P))
    end
end

@testset "Arithmetic" begin

    for  P in Ps
        pNULL = P(Int[])
        p0 = P([0])
        p1 = P([1,0,0,0,0,0,0,0,0,0,0,0,0,0])
        p2 = P([1,1,0,0])
        p3 = P([1,2,1,0,0,0,0])
        p4 = P([1,3,3,1,0,0])
        p5 = P([1,4,6,4,1,0,0,0,0,0,0,0,0,0,0,0,0,0])
        pN = P([276,3,87,15,24,0])
        pR = P([3 // 4, -2 // 1, 1 // 1])

        @test p3 == P([1,2,1])
        @test pN * 10 == P([2760, 30, 870, 150, 240])
        @test pN / 10.0 == P([27.6, 0.3, 8.7, 1.5, 2.4])
        @test 10 * pNULL + pN == pN
        @test 10 * p0 + pN == pN
        @test p5 + 2 * p1 == P([3,4,6,4,1])
        @test 10 * pNULL - pN == -pN
        @test p0 - pN == -pN
        @test p5 - 2 * p1 == P([-1,4,6,4,1])
        @test p2 * p2 * p2 == p4
        @test p2^4 == p5
        @test pNULL^3 == pNULL
        @test pNULL * pNULL == pNULL

        @test pNULL + 2 == p0 + 2 == 2 + p0 == P([2])
        @test p2 - 2 == -2 + p2 == P([-1,1])
        @test 2 - p2 == P([1,-1])

    end

    for P in Ps # ensure promotion of scalar +,*,/
        p = P([1,2,3])
        @test p + 0.5 == P([1.5, 2.0, 3.0])
        @test p / 2  == P([1/2, 1.0, 3/2])
        @test p * 0.5 == P([1/2, 1.0, 3/2])
    end
end

@testset "Divrem" begin
    for P in  Ps

        p0 = P([0])
        p1 = P([1])
        p2 = P([5, 6, -3, 2 ,4])
        p3 = P([7, -3, 2, 6])
        p4 = p2 * p3
        pN = P([276,3,87,15,24,0])
        pR = P([3 // 4, -2 // 1, 1 // 1])

        @test divrem(p4, p2) == (p3, zero(p3))
        @test p3 % p2 == p3
        @test all((map(abs, coeffs(p2 ÷ p3 - P([1 / 9,2 / 3])))) .< eps())
        @test divrem(p0, p1) == (p0, p0)
        @test divrem(p1, p1) == (p1, p0)
        @test divrem(p2, p2) == (p1, p0)
        @test divrem(pR, pR) == (one(pR), zero(pR))
        @test_throws DivideError p1 ÷ p0
        @test_throws DivideError divrem(p0, p0)
    end
end

@testset "Comparisons" begin
    for P in Ps
        pX = P([1, 2, 3, 4, 5])
        pS1 = P([1, 2, 3, 4, 5], "s")
        pS2 = P([1, 2, 3, 4, 5], 's')
        pS3 = P([1, 2, 3, 4, 5], :s)
        @test pX != pS1
        @test pS1 == pS2
        @test pS1 == pS3
        @test_throws ErrorException pS1 + pX
        @test_throws ErrorException pS1 - pX
        @test_throws ErrorException pS1 * pX
        @test_throws ErrorException pS1 ÷ pX
        @test_throws ErrorException pS1 % pX

        # Testing copying.
        pcpy1 = P([1,2,3,4,5], :y)
        pcpy2 = copy(pcpy1)
        @test pcpy1 == pcpy2

        # Check for isequal
        p1 = P([1.0, -0.0, 5.0, Inf])
        p2 = P([1.0,  0.0, 5.0, Inf])
        p3 = P([0, NaN])

        P != SparsePolynomial && (@test p1 == p2 && !isequal(p1, p2))  # SparsePolynomial doesn't store -0.0,  0.0.
        @test p3 === p3 && p3 ≠ p3 && isequal(p3, p3)

        p = fromroots(P, [1,2,3])
        q = fromroots(P, [1,2,3])
        @test hash(p) == hash(q)

        p1s = P([1,2], :s)
        p1x = P([1,2], :x)
        p2s = P([1], :s)

        @test p1s == p1s
        @test p1s ≠ p1x
        @test p1s ≠ p2s

        @test_throws ErrorException p1s ≈ p1x
        @test p1s ≉ p2s
        @test p1s ≈ P([1,2.], :s)

        @test p2s ≈ 1.0 ≈ p2s
        @test p2s == 1.0 == p2s
        @test p2s ≠ 2.0 ≠ p2s
        @test p1s ≠ 2.0 ≠ p1s

        @test nnz(map(P, sparse(1.0I, 5, 5))) == 5

        @test P([0.5]) + 2 == P([2.5])
        @test 2 - P([0.5]) == P([1.5])
    end
end

@testset "Fitting" begin
    for P in Ps
        xs = range(0, stop = π, length = 10)
        ys = sin.(xs)

        p = fit(P, xs, ys)
        y_fit = p.(xs)
        abs_error = abs.(y_fit .- ys)
        @test maximum(abs_error) <= 0.03

        p = fit(P, xs, ys, 2)
        y_fit = p.(xs)
        abs_error = abs.(y_fit .- ys)
        @test maximum(abs_error) <= 0.03

        # Test weighted
        for W in [1, ones(size(xs)), diagm(0 => ones(size(xs)))]
            p = fit(P, xs, ys, 2, weights = W)
            @test p.(xs) ≈ y_fit
        end

        # Getting error on passing Real arrays to polyfit #146
        xx = Real[20.0, 30.0, 40.0]
        yy = Real[15.7696, 21.4851, 28.2463]
        fit(P, xx, yy, 2)

    end
end

@testset "Values" begin
    for P in Ps
        pNULL = P(Int[])
        p0 = P([0])
        p1 = P([1,0,0,0,0,0,0,0,0,0,0,0,0,0])
        p2 = P([1,1,0,0])
        p3 = P([1,2,1,0,0,0,0])
        p4 = P([1,3,3,1,0,0])
        p5 = P([1,4,6,4,1,0,0,0,0,0,0,0,0,0,0,0,0,0])
        pN = P([276,3,87,15,24,0])
        pR = P([3 // 4, -2 // 1, 1 // 1])

        @test fromroots(P, Int[])(2.) == 1.
        @test pN(-.125) == 276.9609375
        @test pNULL(10) == 0
        @test p0(-10) == 0
        @test fromroots(P, [1 // 2, 3 // 2])(1 // 2) == 0 // 1

        # Check for Inf/NaN operations
        p1 = P([Inf, Inf])
        p2 = P([0, Inf])
        @test p1(Inf) == Inf
        @test isnan(p1(-Inf))
        @test isnan(p1(0))
        @test p2(-Inf) == -Inf

        # issue #189
        p = P([0,1,2,3])
        A = [0 1; 0  0];
        @test  p(A) == A  + 2A^2 + 3A^3

        # issue #209
        ps  = [P([0,1]), P([0,0,1])]
        @test Polynomials.evalpoly.(1/2, ps) ≈ [p(1/2)  for  p  in ps]
    end

    # constant polynomials and type
    Ts = (Int, Float32, Float64, Complex{Int}, Complex{Float64})
    for P in (Polynomial, ImmutablePolynomial{1}, SparsePolynomial)
        for T in Ts
            for S in Ts
                c = 2
                p = P{T}(c)
                x = one(S)
                y = p(x)
                @test y === c * one(T)*one(S)
                q = P{T}([c,c])
                @test typeof(q(x)) == typeof(p(x))
            end
        end
    end

    for P in Ps
        p = P(1)
        x = [1 0; 0 1]
        y = p(x)
        @test y == x

        # Issue #208 and  type of output
        p1=P([1//1])
        p2=P([0, 0.9])
        p3=p1(p2)
        @test isa(p3, P)
        @test eltype(p3) == eltype(p2)
    end

end

@testset "Conversion" begin
    for P in (Polynomial, ImmutablePolynomial{2})
        p = P([0,one(Float64)])
        @test P{Complex{Float64}} == typeof(p + 1im)
        @test P{Complex{Float64}} == typeof(1im - p)
        @test P{Complex{Float64}} == typeof(p * 1im)
    end

    # unnecessary copy in convert #65
    p1 = Polynomial([1,2])
    p2 = convert(Polynomial{Int64}, p1)
    p2[3] = 3
    @test p1[3] == 3

end

@testset "Roots" begin
    for P in Ps

        pNULL = P(Int[])
        p0 = P([0])
        p1 = P([1,0,0,0,0,0,0,0,0,0,0,0,0,0])
        p2 = P([1,1,0,0])
        p3 = P([1,2,1,0,0,0,0])
        p4 = P([1,3,3,1,0,0])
        p5 = P([1,4,6,4,1,0,0,0,0,0,0,0,0,0,0,0,0,0])
        pN = P([276,3,87,15,24,0])
        pR = P([3 // 4, -2 // 1, 1 // 1])

        # From roots
        r = [2, 3]
        @test fromroots(r) == Polynomial([6, -5, 1])
        p = fromroots(P, r)
        @test p == P([6, -5, 1])
        @test sort(roots(p)) ≈ r

        @test roots(p0) == roots(p1) == roots(pNULL) == []
        @test P == LaurentPolynomial ? roots(variable(P)) == [0.0] : roots(P([0,1,0])) == [0.0]

        @test roots(p2) == [-1]
        a_roots = [c for c in coeffs(copy(pN))]
        @test all(map(abs, sort(roots(fromroots(a_roots))) - sort(a_roots)) .< 1e6)
        @test length(roots(p5)) == 4
        @test roots(pNULL) == []
        @test sort(roots(pR)) == [1 // 2, 3 // 2]

        A = [1 0; 0 1]
        @test fromroots(A) == Polynomial(Float64[1, -2, 1])
        p = fromroots(P, A)
        @test p == P(Float64[1, -2, 1])
        @test roots(p) ≈ sort!(eigvals(A), rev = true)

        x = variable()
        plarge = 8.362779449448982e41 - 2.510840694154672e57x + 4.2817430781178795e44x^2 - 1.6225927682921337e31x^3 + 1.0x^4  # #120
        @test length(roots(plarge)) == 4
    end
end

@testset "Integrals and Derivatives" begin
    # Integrals derivatives
    for P in Ps

        pNULL = P(Int[])
        p0 = P([0])
        p1 = P([1,0,0,0,0,0,0,0,0,0,0,0,0,0])
        p2 = P([1,1,0,0])
        p3 = P([1,2,1,0,0,0,0])
        p4 = P([1,3,3,1,0,0])
        p5 = P([1,4,6,4,1,0,0,0,0,0,0,0,0,0,0,0,0,0])
        pN = P([276,3,87,15,24,0])
        pR = P([3 // 4, -2 // 1, 1 // 1])

        c = [1, 2, 3, 4]
        p = P(c)

        der = derivative(p)
        @test coeffs(der) ==ᵗ⁰ [2, 6, 12]
        int = integrate(der, 1)
        @test coeffs(int) ==ᵗ⁰ c

        @test derivative(pR) == P([-2 // 1,2 // 1])
        @test derivative(p3) == P([2,2])
        @test derivative(p1) == derivative(p0) == derivative(pNULL) == pNULL
        @test_throws ErrorException derivative(pR, -1)
        @test integrate(P([1,1,0,0]), 0, 2) == 4.0

        if P ∈ (Polynomial,)
            @test derivative(integrate(pN)) == convert(P{Float64}, pN)
            @test integrate(pNULL, 1) == convert(P{Float64}, p1)
            rc = Rational{Int64}[1,2,3]
            @test integrate(P(rc)) == P{eltype(rc)}([0, 1, 1, 1])
        end

        for i in 1:10
            p = P(rand(1:5, 6))
            @test degree(truncate(p - integrate(derivative(p)), atol=1e-13)) <= 0
            @test degree(truncate(p - derivative(integrate(p)), atol=1e-13)) <= 0
        end

        # Handling of `NaN`s
        p     = P([NaN, 1, 5])
        pder  = derivative(p)
        pint  = integrate(p)

        @test isnan(p(1)) # p(1) evaluates to NaN
        @test isequal(pder, P([NaN]))
        @test isequal(pint, P([NaN]))

        pint  = integrate(p, 0.0im)
        @test isequal(pint, P([NaN]))

        # Issue with overflow and polyder Issue #159
        @test derivative(P(BigInt[0, 1])^100, 100) == P(factorial(big(100)))
    end
end

@testset "Elementwise Operations" begin
    for P in Ps
        p1  = P([1, 2])
        p2  = P([3, 1.])
        p   = [p1, p2]
        q   = [3, p1]
        @test q isa Vector{typeof(p1)}
        @test p isa Vector{typeof(p2)}

        psum  = p .+ 3
        pprod = p .* 3
        pmin  = p .- 3
        @test psum  isa Vector{typeof(p2)}
        @test pprod isa Vector{typeof(p2)}
        @test pmin  isa Vector{typeof(p2)}
    end
end

@testset "Chop and Truncate" begin
    # chop and truncate
    for P in Ps
        if P == Polynomial
            p = P([1, 1, 1, 1])
            coeffs(p)[end] = 0
            @assert coeffs(p) == [1, 1, 1, 0]
            p = chop(p)
        else
            p = P([1, 1, 1, 0])
        end

        @test coeffs(p) ==ᵗ⁰ [1, 1, 1]
        ## truncation
        p1 = P([1,1] / 10)
        p2 = P([1,2] / 10)
        p3 = P([1,3] / 10)
        psum = p1 + p2 - p3
        @test degree(psum) == 1         # will have wrong degree
        @test degree(truncate(psum)) == 0 # the degree should be correct after truncation

        @test truncate(P([2,1]), rtol = 1 / 2, atol = 0) == P([2])
        @test truncate(P([2,1]), rtol = 1, atol = 0)   == P([0])
        @test truncate(P([2,1]), rtol = 0, atol = 1)   == P([2])

        pchop = P([1, 2, 3, 0, 0, 0])
        pchopped = chop(pchop)
        @test roots(pchop) == roots(pchopped)

    end
end

@testset "As matrix elements" begin

    for P in Ps
        p = P([1,2,3], :x)
        A = [1 p; p^2 p^3]
        @test !issymmetric(A)
        @test issymmetric(A*transpose(A))
        diagm(0 => [1, p^3], 1=>[p^2], -1=>[p])
    end

    # issue 206 with mixed variable types and promotion
    for P in (ImmutablePolynomial,)
        λ = P([0,1],:λ)
        A = [1 λ; λ^2 λ^3]
        @test A ==  diagm(0 => [1, λ^3], 1=>[λ], -1=>[λ^2])
        @test all([1 -λ]*[λ^2 λ; λ 1] .== 0)
        @test [λ 1] + [1 λ] == (λ+1) .* [1 1] # (λ+1) not a number, so we broadcast
    end
end

@testset "Linear Algebra" begin
    for P in Ps
        p = P([3, 4])
        @test norm(p) == 5
        p = P([-1, 3, 5, -2])
        @test norm(p) ≈ 6.244997998398398
        p = P([1 - 1im, 2 - 3im])
        p2 = conj(p)
        @test coeffs(p2) ==ᵗ⁰ [1 + 1im, 2 + 3im]
        @test transpose(p) == p
        P != ImmutablePolynomial && @test transpose!(p) == p

        @test norm(P([1., 2.])) == norm([1., 2.])
        @test norm(P([1., 2.]), 1) == norm([1., 2.], 1)
    end
end

@testset "Indexing" begin
    # Indexing
    for P in Ps
        # getindex
        p = P([-1, 3, 5, -2])
        @test p[0] == -1
        @test p[[1, 2]] ==ᵗ⁰ [3, 5]
        @test p[1:2] ==ᵗ⁰ [3, 5]
        @test p[:] ==ᵗ⁰ [-1, 3, 5, -2]

        if P != ImmutablePolynomial
            # setindex
            p1  = P([1,2,1])
            p1[5] = 1
            @test p1[5] == 1
            @test p1 == P([1,2,1,0,0,1])

            @test p[end] == coeffs(p)[end]

            if P != SparsePolynomial
                p[1] = 2
                @test coeffs(p)[2] == 2
                p[2:3] = [1, 2]
                @test coeffs(p)[3:4] == [1, 2]
                p[0:1] = 0
                @test coeffs(p)[1:2] ==ᵗ⁰ [0, 0]

                p[:] = 1
                @test coeffs(p) ==ᵗ⁰ ones(4)
            end

            p[:] = 0
            @test chop(p) ≈ zero(p)
        end

        p1 = P([1,2,0,3])
        for term in p1
            @test isa(term, P)
        end

        @test eltype(p1) == Int
        for Q in (Polynomial, ImmutablePolynomial{5})
            p1 = Q([1,2,0,3])
            @test eltype(collect(p1)) == Q{Int}
            @test eltype(collect(Q{Float64}, p1)) == Q{Float64}
            @test_throws InexactError collect(Q{Int}, Q([1.2]))
        end

        p1 = P([1,2,0,3])
        @test length(collect(p1)) == degree(p1) + 1

        @test [p1[idx] for idx in eachindex(p1)] ==ᵗᶻ [1,2,0,3]
    end
end

@testset "Copying" begin
    for P in Ps
        pcpy1 = P([1,2,3,4,5], :y)
        pcpy2 = copy(pcpy1)
        @test pcpy1 == pcpy2
    end
end

@testset "GCD" begin
    for P in Ps
        p1 = P([2.,5.,1.])
        p2 = P([1.,2.,3.])

        @test degree(gcd(p1, p2)) == 0          # no common roots
        @test degree(gcd(p1, P(5))) == 0          # ditto
        @test degree(gcd(p1, P(eps(0.)))) == 0          # ditto
        @test degree(gcd(p1, P(0))) == degree(p1) # P(0) has the roots of p1
        @test degree(gcd(p1 + p2 * 170.10734737144486, p2)) == 0          # see, c.f., #122

        p1 = fromroots(P, [1.,2.,3.])
        p2 = fromroots(P, [1.,2.,6.])
        res = roots(gcd(p1, p2))
        @test 1. ∈ res
        @test 2. ∈ res
    end
end

@testset "Showing" begin

    p = Polynomial{Rational}([1, 4])
    @test sprint(show, p) == "Polynomial(1//1 + 4//1*x)"

    p = Polynomial{Rational{Int}}([1, 4])
    @test sprint(show, p) == "Polynomial(1//1 + 4//1*x)"

    for P in (Polynomial, ImmutablePolynomial)
        p = P([1, 2, 3])
        @test sprint(show, p) == "$P(1 + 2*x + 3*x^2)"

        p = P([1.0, 2.0, 3.0])
        @test sprint(show, p) == "$P(1.0 + 2.0*x + 3.0*x^2)"

        p = P([1 + 1im, -2im])
        @test sprint(show, p) == "$P((1 + 1im) - 2im*x)"

        p = P([1,2,3,1])  # leading coefficient of 1
        @test repr(p) == "$P(1 + 2*x + 3*x^2 + x^3)"
        p = P([1.0, 2.0, 3.0, 1.0])
        @test repr(p) == "$P(1.0 + 2.0*x + 3.0*x^2 + 1.0*x^3)"
        p = P([1, im])
        @test repr(p) == "$P(1 + im*x)"
        p = P([1 + im, 1 - im, -1 + im, -1 - im])# minus signs
        @test repr(p) == "$P((1 + 1im) + (1 - 1im)*x - (1 - 1im)*x^2 - (1 + 1im)*x^3)"
        p = P([1.0, 0 + NaN * im, NaN, Inf, 0 - Inf * im]) # handle NaN or Inf appropriately
        @test repr(p) == "$P(1.0 + NaN*im*x + NaN*x^2 + Inf*x^3 - Inf*im*x^4)"

        p = P([1,2,3])

        @test repr("text/latex", p) == "\$1 + 2\\cdot x + 3\\cdot x^{2}\$"
        p = P([1 // 2, 2 // 3, 1])
        @test repr("text/latex", p) == "\$\\frac{1}{2} + \\frac{2}{3}\\cdot x + x^{2}\$"

        # customized printing with printpoly
        function printpoly_to_string(args...; kwargs...)
            buf = IOBuffer()
            printpoly(buf, args...; kwargs...)
            return String(take!(buf))
        end
        @test printpoly_to_string(P([1,2,3], "y")) == "1 + 2*y + 3*y^2"
        @test printpoly_to_string(P([1,2,3], "y"), descending_powers = true) == "3*y^2 + 2*y + 1"
        @test printpoly_to_string(P([2, 3, 1], :z), descending_powers = true, offset = -2) == "1 + 3*z^-1 + 2*z^-2"
        @test printpoly_to_string(P([-1, 0, 1], :z), offset = -1, descending_powers = true) == "z - z^-1"
    end
end

@testset "Plotting" begin
    p = fromroots([-1, 1]) # x^2 - 1
    r = -1.4:0.055999999999999994:1.4
    rec = apply_recipe(Dict{Symbol,Any}(), p)
    @test getfield(rec[1], 1) == Dict{Symbol,Any}(:label => "-1 + x^2")
    @test rec[1].args == (r, p.(r))

    r = -1:0.02:1
    rec = apply_recipe(Dict{Symbol,Any}(), p, -1, 1)
    @test getfield(rec[1], 1) == Dict{Symbol,Any}(:label => "-1 + x^2")
    @test rec[1].args == (r, p.(r))

    p = ChebyshevT([1,1,1])
    rec = apply_recipe(Dict{Symbol,Any}(), p)
    r = -1.0:0.02:1.0  # uses domain(p)
    @test rec[1].args == (r, p.(r))

end

@testset "Promotion"  begin

    # Test different types work together
    Ps = (Polynomial,  ImmutablePolynomial, SparsePolynomial, LaurentPolynomial)
    for P₁ in Ps
        for   P₂ in Ps
            p₁, p₂ = P₁(rand(1:5, 4)), P₂(rand(1:5, 5))
            p₁ + p₂
            p₁ * p₂

            p₁, p₂ = P₁(rand(1:5, 4)), P₂(5) # constant
            p₁ + p₂
            p₁ * p₂

            p₁, p₂ = P₁(rand(1:5, 4)), P₂(5, :y) # constant, but wrong variable
            if !(promote_type(P₁, P₂) <: Polynomial || promote_type(P₁, P₂) <: Polynomials.StandardBasisPolynomial)
                p₁ + p₂
                p₁ * p₂
            end
        end
    end

    # P{T}(vector{S}) will promote to P{T}
    for Ts in ((Int32, Int,  BigInt),
               (Int,  Rational{Int}, Float64),
               (Float32, Float64, BigFloat)
               )

        n = length(Ts)
        for i in 1:n-1
            T1,T2 = Ts[i],Ts[i+1]
            for P in (Polynomial,  ImmutablePolynomial{5}, SparsePolynomial, LaurentPolynomial)
                p = P{T2}(T1.(rand(1:3,3)))
                @test typeof(p) == P{T2}
            end
        end
    end

    # test P{T}(...) is P{T}
    for P in (Polynomial,  ImmutablePolynomial{5}, SparsePolynomial, LaurentPolynomial)
        for  T in (Int32, Int64, BigInt)
            p₁ =  P{T}(Float64.(rand(1:3,5)))
            @test typeof(p₁) == P{T} # conversion works
            @test_throws InexactError  P{T}(rand(5))
        end
    end

end


text:
module Fuzzy

    export TriangularMF, GaussianMF, BellMF, TrapezoidalMF, SigmoidMF

    export Rule, FISMamdani, FISSugeno

    export mean_at, eval_fis

    export minimum_value, algebraic_product, bounded_difference, drastic_product, einstein_product, hamacher_product

    export maximum_value, algebraic_sum, bounded_sum, drastic_sum, einstein_sum, hamacher_sum

    # T-Norm
    include("T-norm.jl")

    # S-Norm
    include("S-norm.jl")

    # Membership functions
    include("MF.jl")

    # Membership functions evaluations
    include("EvalMF.jl")

    # FIS
    include("FIS.jl")

    # Evaluation functions
    include("Eval.jl")

end


text:
using ITF1788
using Test

@testset "parse line" begin

    # simple 1 output commands
    @test parse_command("add [10.0, 20.0] [13.0, 17.0] = [23.0, 37.0];") == "@test +(interval(10.0, 20.0), interval(13.0, 17.0)) === Interval(23.0, 37.0)"
    @test parse_command("add [10.0, 20.0] [13.0, 17.0] = [23.0, 38.0];") == "@test_broken +(interval(10.0, 20.0), interval(13.0, 17.0)) === Interval(23.0, 38.0)"
    @test parse_command("add [10.0, 20.0] [13.0, 17.0] = [23.0, 38.0];"; failure=false) == "@test +(interval(10.0, 20.0), interval(13.0, 17.0)) === Interval(23.0, 38.0)"

    # tests with NaN
    @test parse_command("rad [empty] = NaN;") == "@test isnan(radius(emptyinterval()))"

    # two outputs
    @test parse_command("midRad [0.0,2.0] = 1.0 1.0;") == "@test midpoint_radius(interval(0.0,2.0))[1] === 1.0 && midpoint_radius(interval(0.0,2.0))[2] === 1.0"

    # texts
    @test parse_command("""b-textToInterval "[1.e-3, 1.1e-3]" = [0X4.189374BC6A7ECP-12, 0X4.816F0068DB8BCP-12];""") == "@test @interval(\"[1.e-3, 1.1e-3]\") === Interval(0x4.189374BC6A7ECP-12, 0x4.816F0068DB8BCP-12)"
    @test parse_command("""d-textToInterval "[entire]" = [-infinity, +infinity]_dac;""") == "@test @decorated(\"[entire]\") === DecoratedInterval(Interval(-Inf, +Inf), dac)"
    @test parse_command("""d-textToInterval "[nai]" = [nai];""") == "@test isnai(@decorated(\"[nai]\"))"
end


text:
using DaemonMode

runargs()


text:
using Documenter, EvolutionaryPlus

makedocs(;
    modules=[EvolutionaryPlus],
    format=Documenter.HTML(),
    pages=[
        "Home" => "index.md",
        "Manual" => Any[
            "manual/start.md",
            "manual/example.md",
        ],
        "Algorithm" => Any[
            "algorithm/DE.md",
            "algorithm/CMAES.md",
            "algorithm/GPs.md",
        ],
        "Library" => Any[
            "lib/public.md",
            "lib/internal.md",
        ],
        "contributing.md",
    ],
    repo="https://github.com/LDNN97/EvolutionaryPlus.jl/blob/{commit}{path}#L{line}",
    sitename="EvolutionaryPlus.jl",
    authors="LDNN97",
)

deploydocs(;
    repo="github.com/LDNN97/EvolutionaryPlus.jl",
)


text:
using Revise
	using DiffEqOperators, ForwardDiff
	using PseudoArcLengthContinuation, LinearAlgebra, Plots, SparseArrays, Parameters, Setfield
	const PALC = PseudoArcLengthContinuation

norminf = x -> norm(x, Inf)
plotsol!(x, nx = Nx, ny = Ny; kwargs...) = heatmap!(LinRange(0,1,nx), LinRange(0,1,ny), reshape(x, nx, ny)'; color = :viridis, xlabel = "x", ylabel = "y", kwargs...)
plotsol(x, nx = Nx, ny = Ny; kwargs...) = (plot();plotsol!(x, nx, ny; kwargs...))

function Laplacian2D(Nx, Ny, lx, ly, bc = :Neumann)
	hx = 2lx/Nx
	hy = 2ly/Ny
	D2x = CenteredDifference(2, 2, hx, Nx)
	D2y = CenteredDifference(2, 2, hy, Ny)

	Qx = Neumann0BC(hx)
	Qy = Neumann0BC(hy)

	D2xsp = sparse(D2x * Qx)[1]
	D2ysp = sparse(D2y * Qy)[1]
	A = kron(sparse(I, Ny, Ny), D2xsp) + kron(D2ysp, sparse(I, Nx, Nx))
	return A, D2x
end

ϕ(u, λ)  = -10(u-λ*exp(u))
dϕ(u, λ) = -10(1-λ*exp(u))

function NL!(dest, u, p)
	@unpack λ = p
	dest .= ϕ.(u, λ)
	return dest
end

NL(u, p) = NL!(similar(u), u, p)

function Fmit!(f, u, p)
	mul!(f, p.Δ, u)
	f .= f .+ NL(u, p)
	return f
end

Fmit(u, p) = Fmit!(similar(u), u, p)

function dFmit(x, p, dx)
	f = similar(dx)
	mul!(f, p.Δ, dx)
	nl = d1NL(x, p, dx)
	f .= f .+ nl
end

function JFmit(x,p)
	J = p.Δ
	dg = dϕ.(x, p.λ)
	return J + spdiagm(0 => dg)
end

# computation of the derivatives
d1NL(x, p, dx) = ForwardDiff.derivative(t -> NL(x .+ t .* dx, p), 0.)
d1Fmit(x, p, dx) = ForwardDiff.derivative(t -> Fmit(x .+ t .* dx, p), 0.)
d2Fmit(x, p, dx1, dx2) = ForwardDiff.derivative(t2 -> ForwardDiff.derivative( t1 -> Fmit(x .+ t1 .* dx1 .+ t2 .* dx2, p), 0.), 0.)
####################################################################################################
Nx = 100
	Ny = 101
	lx = 0.5
	ly = 0.5

	Δ = Laplacian2D(Nx, Ny, lx, ly)[1]
	par_mit = (λ = .05, Δ = Δ)
	sol0 = 0*ones(Nx, Ny) |> vec
####################################################################################################
eigls = EigArpack(0.5, :LM)
	opt_newton = PALC.NewtonPar(tol = 1e-8, verbose = true, eigsolver = eigls, maxIter = 20)
	out, hist, flag = @time PALC.newton(Fmit, JFmit, sol0, par_mit, opt_newton, normN = norminf)

plotsol(out)

####################################################################################################
opts_br = ContinuationPar(dsmin = 0.001, dsmax = 0.05, ds = 0.01, pMax = 3.5, pMin = 0.025, detectBifurcation = 2, nev = 30, plotEveryNsteps = 10, newtonOptions = (@set opt_newton.verbose = true), maxSteps = 100, precisionStability = 1e-6, nInversion = 4, dsminBisection = 1e-7, maxBisectionSteps = 25)

	br, _ = @time PALC.continuation(
		Fmit, JFmit,
		sol0, par_mit, (@lens _.λ), opts_br;
		printSolution = (x, p) -> norm(x),
		plotSolution = (x, p; kwargs...) -> plotsol!(x ; kwargs...),
		plot = true, verbosity = 3, normC = norminf)

####################################################################################################
# branch switching

D(f, x, p, dx) = ForwardDiff.derivative(t->f(x .+ t .* dx, p), 0.)

d1Fmit(x,p,dx1) = D((z, p0) -> Fmit(z, p0), x, p, dx1)
d2Fmit(x,p,dx1,dx2) = D((z, p0) -> d1Fmit(z, p0, dx1), x, p, dx2)
d3Fmit(x,p,dx1,dx2,dx3) = D((z, p0) -> d2Fmit(z, p0, dx1, dx2), x, p, dx3)

jet = ( Fmit, JFmit, d2Fmit, d3Fmit)

br1, _ = continuation(jet...,
		br, 3,
		setproperties(opts_br;ds = 0.001, maxSteps = 40);
		verbosity = 3, plot = true,
		printSolution = (x, p) -> norm(x),
		plotSolution = (x, p; kwargs...) -> plotsol!(x ; kwargs...),
		normC = norminf)

plot([br,br1],plotfold=false)

br2, _ = continuation(jet...,
		br1, 1, setproperties(opts_br;ds = 0.001, maxSteps = 40);
		verbosity = 3, plot = true,
		printSolution = (x, p) -> norm(x),
		plotSolution = (x, p; kwargs...) -> plotsol!(x ; kwargs...), normC = norminf)

plot([br,br1,br2],plotfold=false)

####################################################################################################
# analyse 2d bifurcation point
bp2d = @time PALC.computeNormalForm(jet..., br, 2;  verbose=true)

PALC.nf(bp2d)[2] |> println

using ProgressMeter
Nd = 200
	L = 0.9
	X = LinRange(-L,L, Nd)
	Y = LinRange(-L,L, Nd)
	P = LinRange(-0.0001,0.0001, Nd+1)

V1a = @showprogress [bp2d(Val(:reducedForm),[x1,y1], p1)[1] for p1 in P, x1 in X, y1 in Y]

Ind1 = findall( abs.(V1a) .<= 9e-4 * maximum(abs.(V1a)))

	V2a = @showprogress [bp2d(Val(:reducedForm),[X[ii[2]],Y[ii[3]]], P[ii[1]])[2] for ii in Ind1]

	Ind2 = findall( abs.(V2a) .<= 3e-3 * maximum(abs.(V2a)))
	@show length(Ind2)

resp = Float64[]
	resx = Vector{Float64}[]
	resnrm = Float64[]
	@showprogress for k in Ind2
		ii = Ind1[k]
		push!(resp, P[ii[1]])
		# push!(resx, max(X[ii[2]],Y[ii[3]]))
		push!(resnrm, sqrt(X[ii[2]]^2+Y[ii[3]]^2))
		push!(resx, [X[ii[2]], Y[ii[3]]])
	end

using LaTeXStrings

plot(
	scatter(1e4resp, map(x->x[1], resx), map(x->x[2], resx); label = "", markerstrokewidth=0, xlabel = L"10^4 \cdot \lambda", ylabel = L"x_1", zlabel = L"x_2", zcolor = resnrm, color = :viridis,colorbar=false),
	scatter(1e4resp, resnrm; label = "", markersize =2, markerstrokewidth=0, xlabel = L"10^4 \cdot \lambda", ylabel = L"\|x\|"))

plotsol(bp2d.ζ[1])
plotsol(bp2d(resx[10], resp[10]))
####################################################################################################
# find isolated branch, see Farrell et al.
deflationOp = DeflationOperator(2.0, (x, y) -> dot(x, y), 1.0, [out])
optdef = setproperties(opt_newton; tol = 1e-8, maxIter = 150)

# eigen-elements close to the second bifurcation point on the branch
# of homogenous solutions
vp, ve, _, _= eigls(JFmit(out, @set par_mit.λ = br.bifpoint[2].param), 5)

for ii=1:size(ve, 2)
		outdef1, _, flag, _ = @time newton(
			Fmit, JFmit,
			# initial guess for newton
			br.bifpoint[2].x .+ 0.01 .* ve[:,ii] .* (1 .+ 0.01 .* rand(Nx*Ny)),
			(@set par_mit.λ = br.bifpoint[2].param + 0.005),
			optdef, deflationOp)
			flag && push!(deflationOp, outdef1)
	end
	length(deflationOp)

l = @layout grid(3,2)
	plot(layout = l)
	for ii=1:length(deflationOp)
		plotsol!(deflationOp[ii], title="$ii", subplot = ii, label = "", xlabel="$ii", colorbar=true)
	end
	title!("")

brdef1, _ = @time PALC.continuation(
	Fmit, JFmit,
	deflationOp[3], (@set par_mit.λ = br.bifpoint[2].param + 0.005), (@lens _.λ),
	# bp2d([0.6,0.6], -0.01), br.bifpoint[2].param - 0.005,
	setproperties(opts_br;ds = 0.001, detectBifurcation = 0, dsmax = 0.01, maxSteps = 500);
	verbosity = 3, plot = true,
	printSolution = (x, p) -> norm(x),
	plotSolution = (x, p; kwargs...) -> plotsol!(x ; kwargs...), normC = norminf)

plot([br,br1,br2, brdef1],plotfold=false)

brdef2, _ = @time PALC.continuation(
	Fmit, JFmit,
	deflationOp[5], (@set par_mit.λ = br.bifpoint[2].param + 0.005), (@lens _.λ),
	setproperties(opts_br;ds = -0.001, detectBifurcation = 0, dsmax = 0.02);
	verbosity = 3, plot = true,
	printSolution = (x, p) -> norm(x),
	plotSolution = (x, p; kwargs...) -> plotsol!(x ; kwargs...), normC = norminf)

plot([br,br1,br2, brdef1, brdef2],plotfold=false, putbifptlegend = false)

plot([brdef1, brdef2],plotfold = false, putbifptlegend = false)


text:
export
    showfull,
    edit_list

"""
    showfull(x)
    showfull(io, x)

Print arguments without truncation.
"""
showfull(io, x) = show(IOContext(io; compact = false, limit = false), x)
showfull(x) = showfull(STDOUT, x)

function _edit_list(ms, n)
    if isempty(ms)
        print_with_color(:red, "no methods.\n")
        return
    end
    ms = collect(ms)
    N = length(ms)
    firstmsg = N > n ? ", only first $n of $N shown" : ""
    print_with_color(:yellow, "select a method (C-d exit$firstmsg)\n")
        for (index, method) in enumerate(Base.Iterators.take(ms, n))
        print_with_color(:green, "$(index)) ")
        print(method)
        println()
    end
    while true
        selection = strip(readline())
        if isempty(selection)
            break
        else
            index = tryparse(Int, selection)
            if isnull(index) || (get(index) ∉ 1:length(ms))
                println_with_color(:red, "invalid selection $(selection), ",
                                   "exit with RET")
            else
                edit(ms[get(index)])
                break
            end
        end
    end
end

"""
    edit_list(f; [max)

List methods of a function, edit the selected one.
"""
edit_list(f::Function; max_methods = 20) =
    _edit_list(methods(f), max_methods)

edit_list(f::Function, t::Type; max_methods = 20, showparents = true) =
    _edit_list(methodswith(t, f, showparents), max_methods)


text:
#!/usr/bin/env julia

"""
# Problem 5: Smallest multiple

2520 is the smallest number that can be divided by each of the numbers from 1
to 10 without any remainder.

What is the smallest positive number that is evenly divisible by all of the
numbers from 1 to `max_num`?
"""
module Problem005

"""
    solve(; max_num::Integer = 20)

Returns the solution for [`Problem005`](@ref) with the given parameters.

## Preconditions
- `max_num > 0`
"""
function solve(; max_num::Integer = 20)
    return lcm((oneunit(max_num):max_num)...)
end

end

max_stars_repo_name: curtislb/ProjectEuler

max_issues_repo_name: curtislb/ProjectEuler

max_forks_repo_name: curtislb/ProjectEuler