text:
export chandra, cyclic, katsura, fourbar, rps10, ipp, ipp2, boon, heart,
    d1, bacillus_subtilis, griewank_osborne, tritangents, cyclooctane

"""
    chandra(n)

The Chandrasekhar H-Equation for `n`.

## References

Laureano Gonzalez-Vega: "Some examples on problem solving by using the
  symbolic viewpoint when dealing with polynomial systems of equations".
  in: "Computer Algebra in Science and Engineering", Editors: J. Fleischer,
  J. Grabmeier, F.W. Hehl and W. Kuechlin.  Pages 102-116
  World Scientific Publishing, 1995.

S. Chandrasekhar: "Radiative Transfer", Dover, NY, 1960.

C.T. Kelley: "Solution of the Chandrasekhar H-equation by Newton's method".
  J. Math. Phys., 21 (1980), pp. 1625-1628.

Jorge J. More': "A collection of nonlinear model problems"
  in: "Computational Solution of Nonlinear Systems of Equations",
  Editors: Eugene L. Allgower and Kurt Georg.  Pages 723-762.
  Lectures in Applied Mathematics, Volume 26, AMS, 1990.

NOTE : the parameter c equals 0.51234.
  In general c can be any number in the interval (0,1]
"""
function chandra(n)
    @polyvar H[1:n]

    return TestSystem([(2n*H[i] - 0.51234H[i]*(1+ sum(i / (j+i) * H[j] for j=1:(n-1))) - 2n) for i=1:n],
        multi_bezout_number=(2^(n-1), [H[1:n-1], [H[n]]]),
        mixed_volume=2^(n-1))
end

"""
    cyclic(n)

The cyclic n-roots problems.

## References

Göran Björck and Ralf Fröberg:
`A faster way to count the solutions of inhomogeneous systems
 of algebraic equations, with applications to cyclic n-roots',
in J. Symbolic Computation (1991) 12, pp 329--336.
"""
function cyclic(n)
    @polyvar z[1:n]

    mv = nothing
    if n == 5
        mv = 70
    elseif n == 6
        mv = 156
    elseif n == 7
        mv = 924
    elseif n == 8
        mv = 2560
    elseif n == 9
        mv = 11016
    elseif n == 10
        mv = 35940
    elseif n == 11
        mv = 184_756
    end

    TestSystem([(sum(prod(z[(k-1) % n + 1] for k in j:j+m) for j in 1:n) for m=0:(n-2))...,prod(z)-1],
        mixed_volume = mv)
end

"""
    katsura(n)

A problem of magnetism in physics.

## References

From the PoSSo test suite.

Shigetoshi Katsura: "Users posing problems to PoSSO",
in the PoSSo Newsletter, no. 2, July 1994,
edited by L. Gonzalez-Vega and T. Recio.
Available at http://janet.dm.unipi.it/

S. Katsura, W. Fukuda, S. Inawashiro, N.M. Fujiki and R. Gebauer,
Cell Biophysics, Vol 11, pages 309--319, 1987.

W. Boege, R. Gebauer, and H. Kredel:
"Some examples for solving systems of algebraic equations by
 calculating Groebner bases", J. Symbolic Computation, 2:83-98, 1986.

Shigetoshi Katsura:
 "Spin Glass Problem by the Method of Integral Equation
  of the Effective Field".  In "New Trends in Magnetism",
edited by Mauricio D. Coutinho-Filho and Sergio M. Resende,
pages 110-121, World Scientific, 1990.
"""
function katsura(n)
  @polyvar x[0:n] # This creates variables x0, x1, ...

  return TestSystem([
    (sum(x[abs(l)+1]*x[abs(m-l)+1] for l=-n:n if abs(m-l)<=n) -
    x[m+1] for m=0:n-1)...,
    x[1] + 2sum(x[i+1] for i=1:n) - 1
  ], mixed_volume=2^n)
end

# Mechanisms

"""
    fourbar()
A four-bar design problem, so-called 5-point problem.

## References

See Morgan, A.P. and Wampler, C.W. :
 `Solving a planar four-bar design problem using continuation'
in Transaction of the ASME, J. of Mechanical Design, Vol. 112
pages 544-550, 1990.

For the coefficients, see Table 2, a = (1,0).
This is the start system, with five random precision points,
that has been used to solve twenty other test systems.
"""
function fourbar()
    @polyvar X1 X2 Y1 Y2

    return TestSystem([
        0.01692601*X1^2*Y1^2 - 0.888509280014*X1^2*Y2^2 + 0.0411717692438*X2^2*Y1^2 - 0.00437457395884*X2^2*Y2^2 + 0.331480641249*X1*X2*Y1^2 - 1.38036964668*X1*X2*Y2^2 - 0.270492270191*X1^2*Y1*Y2 + 1.44135801774*X2^2*Y1*Y2 + 0.859888946812*X1*X2*Y1*Y2 + 0.0791489659197*X1^2*Y1 - 0.00336032777032*X1^2*Y2 - 0.0620826738427*X1*Y1^2 + 0.501879647495*X1*Y2^2 + 0.647156236961*X2^2*Y1 + 0.0926311741907*X2^2*Y2 - 0.255000006226*X2*Y1^2 - 0.0896892386081*X2*Y2^2 - 0.568007271041*X1*X2*Y2 + 0.095991501961*X1*X2*Y1 + 0.165310767618*X1*Y1*Y2 - 0.563962321337*X2*Y1*Y2 - 0.0784871167595*X1*Y1 - 0.0784871167595*X2*Y2 + 0.011807283256*X1*Y2  - 0.011807283256*X2*Y1 + 0.0422876985355*X1^2 + 0.0422876985355*X2^2 + 0.0372427422943*Y1^2 + 0.0372427422943*Y2^2,
        0.518178672335*X1^2*Y1^2 - 0.0414464807343*X1^2*Y2^2 + 2.63600135179*X2^2*Y1^2 - 0.799490472298*X2^2*Y2^2 + 0.29442805494*X1^2*Y1*Y2 + 1.46551534655*X2^2*Y1*Y2 - 0.631878110759*X1*X2*Y1^2 - 1.80296540237*X1*X2*Y2^2 - 2.87586667102*X1*X2*Y1*Y2 - 0.987856648177*X1^2*Y1 - 0.530579106676*X1^2*Y2 - 0.0397576281649*X1*Y1^2 + 0.317719102869*X1*Y2^2 - 1.93710490787*X2^2*Y1  + 0.00127693327315*X2^2*Y2 - 0.581380074072*X2*Y1^2 - 0.0672137066743*X2*Y2^2 + 0.531856039949*X1*X2*Y1 + 0.949248259696*X1*X2*Y2 + 0.514166367398*X1*Y1*Y2 - 0.357476731033*X2*Y1*Y2 + 0.140965913657*X1*Y1  + 0.140965913657*X2*Y2 - 0.153347218606*X1*Y2 + 0.153347218606*X2*Y1 + 0.283274882058*X1^2  + 0.283274882058*X2^2 + 0.0382903330079*Y1^2 + 0.0382903330079*Y2^2,
        0.0233560008057*X1^2*Y1^2 - 0.00428427501149*X1^2*Y1*Y2 - 0.792756311827*X1^2*Y2^2 + 0.0492185850289*X2^2*Y2^2 + 0.0759264856293*X1*X2*Y1^2 + 1.14839711492*X1*X2*Y1*Y2 - 0.283066217262*X2^2*Y1^2 + 0.460041521291*X2^2*Y1*Y2 - 0.388399310674*X1*X2*Y2^2 - 0.0561169736293*X1*Y1^2 + 0.485064247792*X1*Y2^2 + 0.0689567235492*X1^2*Y1 - 0.115620658768*X1^2*Y2 - 0.13286905328*X2*Y1^2 - 0.084375901147*X2*Y2^2 + 0.639964831612*X2^2*Y1 + 0.101386684276*X2^2*Y2 + 0.217007343044*X1*X2*Y1 - 0.571008108063*X1*X2*Y2 + 0.0484931521334*X1*Y1*Y2 - 0.541181221422*X2*Y1*Y2 - 0.00363197918253*X1*Y2 + 0.00363197918253*X2*Y1 - 0.0781302968652*X1*Y1 - 0.0781302968652*X2*Y2 + 0.0471311092612*X1^2 + 0.0471311092612*X2^2 + 0.0324495575052*Y1^2 + 0.0324495575052*Y2^2,
        0.393707415641*X1^2*Y1^2 + 0.59841456862*X1^2*Y2^2 + 0.0735854940135*X2^2*Y1^2 + 0.0548997238169*X2^2*Y2^2 + 0.0116156836985*X1^2*Y1*Y2 + 0.0699694273575*X2^2*Y1*Y2 - 0.305757340849*X1*X2*Y1^2 - 0.364111084508*X1*X2*Y2^2 - 0.223392923175*X1*X2*Y1*Y2 + 0.0996725944534*X1^2*Y1 + 0.0113936468426*X1^2*Y2 - 0.381205205249*X1*Y1^2 - 0.473402150235*X1*Y2^2 - 0.0213613191759*X2^2*Y1 - 0.0372595571271*X2^2*Y2 + 0.148904552394*X2*Y1^2 + 0.142408744984*X2*Y2^2 - 0.0486532039697*X1*X2*Y1 + 0.121033913629*X1*X2*Y2 - 0.00649580741066*X1*Y1*Y2 + 0.092196944986*X2*Y1*Y2 - 0.0483106652705*X1*Y1  - 0.0483106652705*X2*Y2 - 0.00316794272326*X1*Y2 + 0.00316794272326*X2*Y1 + 0.00634952598374*X1^2 + 0.00634952598374*X2^2 + 0.0922886309144*Y1^2  + 0.0922886309144*Y2^2
    ], mixed_volume=80)
end

"""
    ipp()

The six-revolute-joint problem of mechanics.

## References

A. Morgan and A. Sommese
`Computing all solutions to polynomial systems
 using homotopy continuation',
 Appl. Math. Comput., Vol. 24, pp 115-138, 1987.
"""
function ipp()
    @polyvar x1 x2 x3 x4 x5 x6 x7 x8
    TestSystem([x1^2+x2^2-1,
        x3^2+x4^2-1,
        x5^2+x6^2-1,
        x7^2+x8^2-1,
        (-2.4915068E-01*x1*x3+ 1.6091354E+00*x1*x4+ 2.7942343E-01 *x2*x3+ 1.4348016E+00*x2*x4) + (4.0026384E-01*x5*x8-8.0052768E-01*x6*x7+ 7.4052388E-02*x1-8.3050031E-02*x2) - (3.8615961E-01*x3-7.5526603E-01*x4+ 5.0420168E-01*x5 -1.0916287E+00*x6+ 4.0026384E-01*x8) + 4.920729E-02,
        (1.2501635E-01*x1*x3-6.8660736E-01*x1*x4-1.1922812E-01* x2*x3-7.1994047E-01*x2*x4) - (4.3241927E-01*x5*x7-8.6483855E-01*x6*x8-3.715727E-02*x1+ 3.5436896E-02*x2)+ 8.5383482E-02*x3-3.9251967E-02*x5-4.3241927E-01*x7+ 1.387301E-02,
        (-6.3555007E-01*x1*x3-1.1571992E-01*x1*x4-6.6640448E-01 *x2*x3) + (1.1036211E-01*x2*x4+ 2.9070203E-01*x5*x7+ 1.2587767E+00*x5*x8)- (6.2938836E-01*x6*x7+ 5.8140406E-01*x6*x8+ 1.9594662E-01*x1)- (1.2280342E+00*x2-7.9034221E-02*x4+ 2.6387877E-02*x5)- 5.713143E-02*x6-1.1628081E+00*x7+1.2587767E+00*x8+ 2.162575E+00,
        (1.4894773E+00*x1*x3+ 2.3062341E-01*x1*x4+ 1.3281073E+00*x2*x3)-(2.5864503E-01*x2*x4+ 1.165172E+00*x5*x7-2.6908494E-01*x5*x8)+ (5.3816987E-01*x6*x7+ 5.8258598E-01*x6*x8-2.0816985E-01*x1)+(2.686832E+00*x2-6.9910317E-01*x3+ 3.5744413E-01*x4)+ 1.2499117E+00*x5+ 1.467736E+00*x6+ 1.165172E+00*x7+ 1.10763397E+00*x8-6.9686809E-01
        ], mixed_volume=64, nreal_solutions=10, nsolutions=32)
end

"""
    ipp2()

The 6R inverse position problem.

# References

This system occurs as Example 3.3 in a paper by Charles Wampler:
`Bezout Number Calculations for Multi-Homogeneous Polynomial Systems',
Appl. Math. Comput. vol. 51 No. 2--3, pp. 143--157.

For the original formulation of the problem, see
 Charles Wampler and Alexander Morgan:
 `Solving the 6R inverse position problem using a generic-case solution
  methodology', Mech. Mach. Theory, Vol. 26, No. 1, pp. 91-106, 1991.
"""
function ipp2()
    @polyvar z21 z22 z31 z32 z33 z41 z42 z43 z51 z52 z53
    TestSystem([ z21^2 + z22^2 - 1.25830472585 + 1.05384271933im,
        z31^2 + z32^2 + z33^2 - 1,
        z41^2 + z42^2 + z43^2 - 1,
        z51^2 + z52^2 + z53^2 - 1,
        z21*z31 + z22*z32 + (0.642935654806 + 0.819555356316im)*z33- 0.266880023988 - 0.452565255666im,
        z31*z41 + z32*z42 + z33*z43 - 0.26425551745 - 0.342483846503im,
        z41*z51 + z42*z52 + z43*z53 - 0.126010863922 - 0.864590917688im,
        (0.352598136811 + 0.116888144319im)*z51 + (0.539042485525 + 0.687058436892im)*z52 + (0.391154215376 + 0.128900893182im)*z53 - 0.179560356712 - 0.8709166566im,
        (0.984138451804 + 0.414967172346im)*z21 + (0.958341609741 + 0.847442419999im)*z22*z33 + (-0.496254299764 - 0.546020011741im)*z22 + (0.353268870498 + 0.389909226888im)*z31 + (0.964759562277 + 0.71397074519im)*z32*z43 + (0.0783739840935 - 1.33026494666im)*z32 + (-0.964759562277 - 0.71397074519im)*z33*z42 + (0.204379350351 + 0.00374294529684im)*z41 + (0.706319991205 + 0.120097702053im)*z42*z53 + (-0.706319991205 - 0.120097702053im)*z43*z52 + (0.907681632783 + 0.405209293447im)*z51 + (0.148939301127 + 0.182393186752im)*z52 + (-0.0486385369088 - 0.496934768083im)*z53 - 0.437312713588 - 0.914780691357im,
        (-0.958341609741 - 0.847442419999im)*z21*z33 + (0.496254299764 + 0.546020011741im)*z21 + (0.984138451804 + 0.414967172346im)*z22 + (-0.964759562277 - 0.71397074519im)*z31*z43 + (-0.0783739840935 + 1.33026494666im)*z31 + (0.353268870498 + 0.389909226888im)*z32 + ( 0.964759562277 + 0.71397074519im)*z33*z41 + (-0.706319991205 - 0.120097702053im)*z41*z53 + (0.204379350351 + 0.00374294529684im)*z42 + (0.706319991205 + 0.120097702053im)*z43*z51 + (-0.148939301127 - 0.182393186752im)*z51 + (0.907681632783 + 0.405209293447im)*z52 + (0.134045297503 + 0.164748774862im)*z53 - 0.719086796333 - 0.691791591267im,
        (0.958341609741 + 0.847442419999im)*z21*z32 + (-0.958341609741 - 0.847442419999im)*z22*z31 + (0.964759562277 + 0.71397074519im)*z31*z42 + (-0.964759562277 - 0.71397074519im)*z32*z41 + (0.353268870498 + 0.389909226888im)*z33 + (0.706319991205 + 0.120097702053im)*z41*z52 + (-0.706319991205 - 0.120097702053im)*z42*z51 + (0.204379350351 + 0.00374294529684im)*z43 + (0.0486385369088 + 0.496934768083im)*z51  + (-0.134045297503 - 0.164748774862im)*z52 + (0.907681632783 + 0.405209293447im)*z53 - 0.64863071126 + 0.983034576618im
    ],
    multi_bezout_number=(320, [[z21, z22, z41, z42, z43], [z31, z32, z33, z51, z52, z53]]),
    mixed_volume=288)
end

include("rps10.jl")

# Mixed

"""
    boon()

Neurofysiology, posted by Sjirk Boon. Note that the system

## References

The system has been posted to the newsgroup
sci.math.num-analysis by Sjirk Boon.

P. Van Hentenryck, D. McAllester and D. Kapur:
`Solving Polynomial Systems Using a Branch and Prune Approach'
SIAM J. Numerical Analysis, Vol. 34, No. 2, pp 797-827, 1997.

"""
function boon()
    @polyvar s1 g1 s2 g2 C1 C2
    TestSystem([s1^2+g1^2 - 1,
     s2^2+g2^2 - 1,
     C1*g1^3+C2*g2^3 - 1.2,
     C1*s1^3+C2*s2^3 - 1.2,
     C1*g1^2*s1+C2*g2^2*s2 - 0.7,
     C1*g1*s1^2+C2*g2*s2^2 - 0.7], mixed_volume=20, nsolutions=8, nreal_solutions=8)
end

"""
    heart()

The heart-dipole problem.

## References

Nelsen, C.V. and Hodgkin, B.C.:
`Determination of magnitudes, directions, and locations of two independent
dipoles in a circular conducting region from boundary potential measurements'
IEEE Trans. Biomed. Engrg. Vol. BME-28, No. 12, pages 817-823, 1981.

Morgan, A.P. and Sommese, A.J.:
`Coefficient-Parameter Polynomial Continuation'
Appl. Math. Comput. Vol. 29, No. 2, pages 123-160, 1989.
Errata: Appl. Math. Comput. 51:207 (1992)

Morgan, A.P. and Sommese, A. and Watson, L.T.:
`Mathematical reduction of a heart dipole model'
J. Comput. Appl. Math. Vol. 27, pages 407-410, 1989.
"""
function heart()
    @polyvar a b c d t u v w
    TestSystem([a + b - 0.63254;
        c + d + 1.34534;
        t*a + u*b - v*c - w*d + 0.8365348;
        v*a + w*b + t*c + u*d - 1.7345334;
        a*t^2 - a*v^2 - 2*c*t*v + b*u^2 - b*w^2 - 2*d*u*w - 1.352352;
        c*t^2 - c*v^2 + 2*a*t*v + d*u^2 - d*w^2 + 2*b*u*w + 0.843453;
        a*t^3 - 3*a*t*v^2 + c*v^3 - 3*c*v*t^2 + b*u^3 - 3*b*u*w^2 + d*w^3 - 3*d*w*u^2 + 0.9563453;
        c*t^3 - 3*c*t*v^2 - a*v^3 + 3*a*v*t^2 + d*u^3 - 3*d*u*w^2 - b*w^3 + 3*b*w*u^2 - 1.2342523
        ], mixed_volume=121, nsolutions=4, nreal_solutions=2)
end

"""
    d1()

A a sparse system, known as benchmark D1.

## References

H. Hong and V. Stahl:
`Safe Starting Regions by Fixed Points and Tightening',
Computing 53(3-4): 322-335, 1995.
"""
function d1()
    @polyvar x1 x2 x3 x4 x5 x6 x7 x8 x9 x10 x11 x12
    TestSystem([x1^2  + x2^2 - 1;
         x3^2  + x4^2 - 1;
         x5^2  + x6^2 - 1;
         x7^2  + x8^2 - 1;
         x9^2  + x10^2 - 1;
         x11^2 + x12^2 - 1;
         3*x3 + 2*x5 + x7 - 3.9701;
         3*x1*x4 + 2*x1*x6 + x1*x8 - 1.7172;
         3*x2*x4 + 2*x2*x6 + x2*x8 - 4.0616;
         x3*x9 + x5*x9 + x7*x9 - 1.9791;
         x2*x4*x9 + x2*x6*x9 + x2*x8*x9 + x1*x10 - 1.9115;
         - x3*x10*x11 - x5*x10*x11 - x7*x10*x11 + x4*x12 + x6*x12 + x8*x12 - 0.4077],
         mixed_volume=192,
         nreal_solutions=16,
         nsolutions=48)
end

"""

    bacillus_subtilis()

The system comes from biochemical reactions (Published here https://journals.plos.org/ploscompbiol/article?id=10.1371/journal.pcbi.1005267).
Basically it models how the bacteris B. subtilis produces a stress response protein, sigmaB, in response to some stress input. This input to the system is the concentration of a phosphatase (a kind of protein). In this system I call that variable phos
"""
function bacillus_subtilis()
    @polyvar w w2 w2v v w2v2 vP sigmaB w2sigmaB vPp phos
    poly = [(-1 * 0.7 * w + -2 * 3600.0 * (w ^ 2 / 2) + 2 * 18.0 * w2)*(0.2 + sigmaB) + 4.0 * 0.4 * (1 + 30.0sigmaB),
     -1 * 0.7 * w2 + 3600.0 * (w ^ 2 / 2) + -1 * 18.0 * w2 + -1 * 3600.0 * w2 * v + 18.0w2v + 36.0w2v + -1 * 3600.0 * w2 * sigmaB + 18.0w2sigmaB,
     -1 * 0.7 * w2v + 3600.0 * w2 * v + -1 * 18.0 * w2v + -1 * 3600.0 * w2v * v + 18.0w2v2 + -1 * 36.0 * w2v + 36.0w2v2 + 1800.0 * w2sigmaB * v + -1 * 1800.0 * w2v * sigmaB,
     (-1 * 0.7 * v + -1 * 3600.0 * w2 * v + 18.0w2v + -1 * 3600.0 * w2v * v + 18.0w2v2 + -1 * 1800.0 * w2sigmaB * v + 1800.0 * w2v * sigmaB + 180.0vPp)*(0.2 + sigmaB) + 4.5 * 0.4 * (1 + 30.0sigmaB),
     -1 * 0.7 * w2v2 + 3600.0 * w2v * v + -1 * 18.0 * w2v2 + -1 * 36.0 * w2v2,
     -1 * 0.7 * vP + 36.0w2v + 36.0w2v2 + -1 * 3600.0 * vP * phos + 18.0vPp,
     (-1 * 0.7 * sigmaB + -1 * 3600.0 * w2 * sigmaB + 18.0w2sigmaB + 1800.0 * w2sigmaB * v + -1 * 1800.0 * w2v * sigmaB)*(0.2 + sigmaB) + 0.4 * (1 + 30.0sigmaB),
     -1 * 0.7 * w2sigmaB + 3600.0 * w2 * sigmaB + -1 * 18.0 * w2sigmaB + -1 * 1800.0 * w2sigmaB * v + 1800.0 * w2v * sigmaB,
     -1 * 0.7 * vPp + 3600.0 * vP * phos + -1 * 18.0 * vPp + -1 * 180.0 * vPp,
     (phos + vPp) - 2.0]
    TestSystem(poly; nsolutions=44)
end

"""

    tritangents()

See https://www.juliahomotopycontinuation.org/examples/tritangents/.
"""
function tritangents()
    @polyvar h[1:3] # variables for the plane
    @polyvar x[1:3] y[1:3] z[1:3] #variables for the contact points
    @polyvar c[1:20] #variables for the cubic

    #the quadric
    Q = x[3] - x[1] * x[2]
    #the cubic
    C = c ⋅ unique(kron([x;1], [x;1], [x;1]))

    #generate the system P for the contact point x
    P_x = [
      h ⋅ x - 1;
      Q;
      C;
      det([h differentiate(Q, x) differentiate(C, x)])
    ]

    #generate a copy of P for the other contact points y,z
    P_y = [p([h; x; c] => [h; y; c]) for p in P_x]
    P_z = [p([h; x; c] => [h; z; c]) for p in P_x]

    #define F
    F = [P_x; P_y; P_z]

    #create random complex coefficients for C
    # c₁ is result of randn(ComplexF64, 20)
    c₁ = Complex{Float64}[-0.524114+0.798188im, 0.23203+0.203011im, 0.24588+0.265859im, -0.0293688-0.0838781im, 1.29037+1.19918im, 0.225558+0.113139im, 0.757771+0.327947im, 0.716762-1.14998im, -0.721386+0.253249im, 0.22066+0.286244im, 1.06596-1.24645im, -0.287841-0.15352im, 0.980128+0.798647im, -0.918701-1.5926im, -0.325193-0.163909im, 0.0885901-0.437345im, 0.195358-1.07946im, 0.451222-0.677876im, -0.498035+0.368011im, -1.34059+1.35779im]
    #plug in c₁ for c
    poly = [f([h; x; y; z; c] => [h; x; y; z; c₁]) for f in F]
    TestSystem(poly; nsolutions=720)
end

"""

    griewank_osborne()

This is an illustration of a system for which Newton's method fails near an isolated root.
No matter how close one starts to the multiplicity-three isolated root at the origin, (x,y) = (0,0), Newton's method diverges.

Reference: Griewank & Osborne, 1983
"""
function griewank_osborne()
    @polyvar x y
    TestSystem([(29/16)*x^3 - 2*x*y, x^2 - y])
end

"""

    cyclooctane(c::Float64=2.0)

This is a system that describes the confirmation space of cyclooctane, where the squared distance between neighboring carbon atoms is `c`.  See

    juliahomotopycontinuation.org/examples/cyclooctane/

for a detailed description.
"""
function cyclooctane(c::Float64=2.0)
    @polyvar z[1:3, 1:6]
    z_vec = vec(z)[1:17]
    Z = [zeros(3) z[:,1:5] [z[1,6]; z[2,6]; 0] [√c; 0; 0]]

    F1 = [(Z[:, i] - Z[:, i+1]) ⋅ (Z[:, i] - Z[:, i+1]) - c for i in 1:7]
    F2 = [(Z[:, i] - Z[:, i+2]) ⋅ (Z[:, i] - Z[:, i+2]) - 8c/3 for i in 1:6]
    F3 = (Z[:, 7] - Z[:, 1]) ⋅ (Z[:, 7] - Z[:, 1]) - 8c/3
    F4 = (Z[:, 8] - Z[:, 2]) ⋅ (Z[:, 8] - Z[:, 2]) - 8c/3
    f = [F1; F2; F3; F4]
    TestSystem(f)
end


text:
using Rubin
using Tests
using Elliptic
using HypergeometricFunctions
using Polylogarithms
using SpecialFunctions

@test integrate(x^2*SpecialFunctions.expint(1,b*x), x) == :(-1//3*x^3*SpecialFunctions.expint(-2,b*x)+1//3*x^3*SpecialFunctions.expint(1,b*x))
@test integrate(x*SpecialFunctions.expint(1,b*x), x) == :((1/2)*x^2*SpecialFunctions.expint(1,b*x)+-1//2*x^2*SpecialFunctions.expint(-1,b*x))
@test integrate(SpecialFunctions.expint(1,b*x), x) == :(-1*b^-1*SpecialFunctions.expint(2,b*x))
@test integrate(x^-2*SpecialFunctions.expint(1,b*x), x) == :(x^-1*SpecialFunctions.expint(2,b*x)+-1*x^-1*SpecialFunctions.expint(1,b*x))
@test integrate(x^-3*SpecialFunctions.expint(1,b*x), x) == :((1/2)*x^-2*SpecialFunctions.expint(3,b*x)+-1//2*x^-2*SpecialFunctions.expint(1,b*x))
@test integrate(x^-4*SpecialFunctions.expint(1,b*x), x) == :(-1//3*x^-3*SpecialFunctions.expint(1,b*x)+1//3*x^-3*SpecialFunctions.expint(4,b*x))
@test integrate(x^2*SpecialFunctions.expint(2,b*x), x) == :(-1//4*x^3*SpecialFunctions.expint(-2,b*x)+1//4*x^3*SpecialFunctions.expint(2,b*x))
@test integrate(x*SpecialFunctions.expint(2,b*x), x) == :(-1//3*x^2*SpecialFunctions.expint(-1,b*x)+1//3*x^2*SpecialFunctions.expint(2,b*x))
@test integrate(SpecialFunctions.expint(2,b*x), x) == :(-1*b^-1*SpecialFunctions.expint(3,b*x))
@test integrate(x^-1*SpecialFunctions.expint(2,b*x), x) == :(-1*SpecialFunctions.expint(1,b*x)+SpecialFunctions.expint(2,b*x))
@test integrate(x^-3*SpecialFunctions.expint(2,b*x), x) == :(x^-2*SpecialFunctions.expint(3,b*x)+-1*x^-2*SpecialFunctions.expint(2,b*x))
@test integrate(x^-4*SpecialFunctions.expint(2,b*x), x) == :((1/2)*x^-3*SpecialFunctions.expint(4,b*x)+-1//2*x^-3*SpecialFunctions.expint(2,b*x))
@test integrate(x^-5*SpecialFunctions.expint(2,b*x), x) == :(-1//3*x^-4*SpecialFunctions.expint(2,b*x)+1//3*x^-4*SpecialFunctions.expint(5,b*x))
@test integrate(x^2*SpecialFunctions.expint(3,b*x), x) == :(-1//5*x^3*SpecialFunctions.expint(-2,b*x)+1//5*x^3*SpecialFunctions.expint(3,b*x))
@test integrate(x*SpecialFunctions.expint(3,b*x), x) == :(-1//4*x^2*SpecialFunctions.expint(-1,b*x)+1//4*x^2*SpecialFunctions.expint(3,b*x))
@test integrate(SpecialFunctions.expint(3,b*x), x) == :(-1*b^-1*SpecialFunctions.expint(4,b*x))
@test integrate(x^-1*SpecialFunctions.expint(3,b*x), x) == :((1/2)*SpecialFunctions.expint(3,b*x)+-1//2*SpecialFunctions.expint(1,b*x))
@test integrate(x^-2*SpecialFunctions.expint(3,b*x), x) == :(x^-1*SpecialFunctions.expint(3,b*x)+-1*x^-1*SpecialFunctions.expint(2,b*x))
@test integrate(x^-4*SpecialFunctions.expint(3,b*x), x) == :(x^-3*SpecialFunctions.expint(4,b*x)+-1*x^-3*SpecialFunctions.expint(3,b*x))
@test integrate(x^-5*SpecialFunctions.expint(3,b*x), x) == :((1/2)*x^-4*SpecialFunctions.expint(5,b*x)+-1//2*x^-4*SpecialFunctions.expint(3,b*x))
@test integrate(x^-6*SpecialFunctions.expint(3,b*x), x) == :(-1//3*x^-5*SpecialFunctions.expint(3,b*x)+1//3*x^-5*SpecialFunctions.expint(6,b*x))
@test integrate(x^3*SpecialFunctions.expint(-1,b*x), x) == :((1/2)*x^4*SpecialFunctions.expint(-1,b*x)+-1//2*x^4*SpecialFunctions.expint(-3,b*x))
@test integrate(x^2*SpecialFunctions.expint(-1,b*x), x) == :(x^3*SpecialFunctions.expint(-1,b*x)+-1*x^3*SpecialFunctions.expint(-2,b*x))
@test integrate(x*SpecialFunctions.expint(-1,b*x), x) == :(b^-2*expintegralei(-1*b*x)+-1*b^-2*exp(-1*b*x))
@test integrate(SpecialFunctions.expint(-1,b*x), x) == :(-1*b^-2*x^-1*exp(-1*b*x))
@test integrate(x^-1*SpecialFunctions.expint(-1,b*x), x) == :((1/2)*SpecialFunctions.expint(1,b*x)+-1//2*SpecialFunctions.expint(-1,b*x))
@test integrate(x^-2*SpecialFunctions.expint(-1,b*x), x) == :(-1//3*x^-1*SpecialFunctions.expint(-1,b*x)+1//3*x^-1*SpecialFunctions.expint(2,b*x))
@test integrate(x^-3*SpecialFunctions.expint(-1,b*x), x) == :(-1//4*x^-2*SpecialFunctions.expint(-1,b*x)+1//4*x^-2*SpecialFunctions.expint(3,b*x))
@test integrate(x^4*SpecialFunctions.expint(-2,b*x), x) == :((1/2)*x^5*SpecialFunctions.expint(-2,b*x)+-1//2*x^5*SpecialFunctions.expint(-4,b*x))
@test integrate(x^3*SpecialFunctions.expint(-2,b*x), x) == :(x^4*SpecialFunctions.expint(-2,b*x)+-1*x^4*SpecialFunctions.expint(-3,b*x))
@test integrate(x^2*SpecialFunctions.expint(-2,b*x), x) == :(-2*b^-3*exp(-1*b*x)+2*b^-3*expintegralei(-1*b*x)+-1*b^-1*x^2*SpecialFunctions.expint(-1,b*x))
@test integrate(x*SpecialFunctions.expint(-2,b*x), x) == :(x^2*SpecialFunctions.expint(-1,b*x)+-1*x^2*SpecialFunctions.expint(-2,b*x))
@test integrate(SpecialFunctions.expint(-1,b*x), x) == :(-1*b^-2*x^-1*exp(-1*b*x))
@test integrate(x^-1*SpecialFunctions.expint(-2,b*x), x) == :(-1//3*SpecialFunctions.expint(-2,b*x)+1//3*SpecialFunctions.expint(1,b*x))
@test integrate(x^-2*SpecialFunctions.expint(-2,b*x), x) == :(-1//4*x^-1*SpecialFunctions.expint(-2,b*x)+1//4*x^-1*SpecialFunctions.expint(2,b*x))
@test integrate(x^-3*SpecialFunctions.expint(-2,b*x), x) == :(-1//5*x^-2*SpecialFunctions.expint(-2,b*x)+1//5*x^-2*SpecialFunctions.expint(3,b*x))
@test integrate(x^5*SpecialFunctions.expint(-3,b*x), x) == :((1/2)*x^6*SpecialFunctions.expint(-3,b*x)+-1//2*x^6*SpecialFunctions.expint(-5,b*x))
@test integrate(x^4*SpecialFunctions.expint(-3,b*x), x) == :(x^5*SpecialFunctions.expint(-3,b*x)+-1*x^5*SpecialFunctions.expint(-4,b*x))
@test integrate(x^3*SpecialFunctions.expint(-3,b*x), x) == :(-6*b^-4*exp(-1*b*x)+6*b^-4*expintegralei(-1*b*x)+-1*b^-1*x^3*SpecialFunctions.expint(-2,b*x)+-3*b^-2*x^2*SpecialFunctions.expint(-1,b*x))
@test integrate(x^2*SpecialFunctions.expint(-3,b*x), x) == :(x^3*SpecialFunctions.expint(-2,b*x)+-1*x^3*SpecialFunctions.expint(-3,b*x))
@test integrate(x*SpecialFunctions.expint(-3,b*x), x) == :((1/2)*x^2*SpecialFunctions.expint(-1,b*x)+-1//2*x^2*SpecialFunctions.expint(-3,b*x))
@test integrate(SpecialFunctions.expint(-1,b*x), x) == :(-1*b^-2*x^-1*exp(-1*b*x))
@test integrate(x^-1*SpecialFunctions.expint(-3,b*x), x) == :(-1//4*SpecialFunctions.expint(-3,b*x)+1//4*SpecialFunctions.expint(1,b*x))
@test integrate(x^-2*SpecialFunctions.expint(-3,b*x), x) == :(-1//5*x^-1*SpecialFunctions.expint(-3,b*x)+1//5*x^-1*SpecialFunctions.expint(2,b*x))
@test integrate(x^-3*SpecialFunctions.expint(-3,b*x), x) == :(-1//6*x^-2*SpecialFunctions.expint(-3,b*x)+1//6*x^-2*SpecialFunctions.expint(3,b*x))
@test integrate(x^3*SpecialFunctions.expint(-3,b*x), x) == :(-6*b^-4*exp(-1*b*x)+6*b^-4*expintegralei(-1*b*x)+-1*b^-1*x^3*SpecialFunctions.expint(-2,b*x)+-3*b^-2*x^2*SpecialFunctions.expint(-1,b*x))
@test integrate(x^2*SpecialFunctions.expint(-2,b*x), x) == :(-2*b^-3*exp(-1*b*x)+2*b^-3*expintegralei(-1*b*x)+-1*b^-1*x^2*SpecialFunctions.expint(-1,b*x))
@test integrate(x*SpecialFunctions.expint(-1,b*x), x) == :(b^-2*expintegralei(-1*b*x)+-1*b^-2*exp(-1*b*x))
@test integrate(SpecialFunctions.expint(0,b*x), x) == :(b^-1*expintegralei(-1*b*x))
@test integrate(x^m*SpecialFunctions.expint(n,x), x) == :(x^(1+m)*(m+n)^-1*SpecialFunctions.expint(n,x)+-1*x^(1+m)*(m+n)^-1*SpecialFunctions.expint(-1m,x))
@test integrate(x^m*SpecialFunctions.expint(n,b*x), x) == :(x^(1+m)*(m+n)^-1*SpecialFunctions.expint(n,b*x)+-1*x^(1+m)*(m+n)^-1*SpecialFunctions.expint(-1m,b*x))
@test integrate((d*x)^m*SpecialFunctions.expint(n,x), x) == :(d^-1*(d*x)^(1+m)*(m+n)^-1*SpecialFunctions.expint(n,x)+-1*d^-1*(d*x)^(1+m)*(m+n)^-1*SpecialFunctions.expint(-1m,x))
@test integrate((d*x)^m*SpecialFunctions.expint(n,b*x), x) == :(d^-1*(d*x)^(1+m)*(m+n)^-1*SpecialFunctions.expint(n,b*x)+-1*d^-1*(d*x)^(1+m)*(m+n)^-1*SpecialFunctions.expint(-1m,b*x))
@test integrate(x^2*SpecialFunctions.expint(n,b*x), x) == :(x^3*(2+n)^-1*SpecialFunctions.expint(n,b*x)+-1*x^3*(2+n)^-1*SpecialFunctions.expint(-2,b*x))
@test integrate(x*SpecialFunctions.expint(n,b*x), x) == :(x^2*(1+n)^-1*SpecialFunctions.expint(n,b*x)+-1*x^2*(1+n)^-1*SpecialFunctions.expint(-1,b*x))
@test integrate(SpecialFunctions.expint(n,b*x), x) == :(-1*b^-1*SpecialFunctions.expint(1+n,b*x))
@test integrate(x^-1*SpecialFunctions.expint(n,b*x), x) == :((1+-1n)^-1*SpecialFunctions.expint(1,b*x)+-1*(1+-1n)^-1*SpecialFunctions.expint(n,b*x))
@test integrate(x^-2*SpecialFunctions.expint(n,b*x), x) == :(x^-1*(2+-1n)^-1*SpecialFunctions.expint(2,b*x)+-1*x^-1*(2+-1n)^-1*SpecialFunctions.expint(n,b*x))
@test integrate(x^-3*SpecialFunctions.expint(n,b*x), x) == :(x^-2*(3+-1n)^-1*SpecialFunctions.expint(3,b*x)+-1*x^-2*(3+-1n)^-1*SpecialFunctions.expint(n,b*x))
@test integrate((c+d*x)^3*SpecialFunctions.expint(1,a+b*x), x) == :(-1*b^-1*(c+d*x)^3*SpecialFunctions.expint(2,a+b*x)+-6*b^-4*d^3*SpecialFunctions.expint(5,a+b*x)+-6*b^-3*d^2*(c+d*x)*SpecialFunctions.expint(4,a+b*x)+-3*d*b^-2*(c+d*x)^2*SpecialFunctions.expint(3,a+b*x))
@test integrate((c+d*x)^2*SpecialFunctions.expint(1,a+b*x), x) == :(-1*b^-1*(c+d*x)^2*SpecialFunctions.expint(2,a+b*x)+-2*b^-3*d^2*SpecialFunctions.expint(4,a+b*x)+-2*d*b^-2*(c+d*x)*SpecialFunctions.expint(3,a+b*x))
@test integrate((c+d*x)*SpecialFunctions.expint(1,a+b*x), x) == :(-1*d*b^-2*SpecialFunctions.expint(3,a+b*x)+-1*b^-1*(c+d*x)*SpecialFunctions.expint(2,a+b*x))
@test integrate(SpecialFunctions.expint(1,a+b*x), x) == :(-1*b^-1*SpecialFunctions.expint(2,a+b*x))
@test integrate((c+d*x)^-2*SpecialFunctions.expint(1,a+b*x), x) == :(-1*d^-1*(c+d*x)^-1*SpecialFunctions.expint(1,a+b*x)+-1*b*d^-1*(b*c+-1*a*d)^-1*expintegralei(-1a+-1*b*x)+b*d^-1*(b*c+-1*a*d)^-1*exp(-1a+b*c*d^-1)*expintegralei(-1*b*d^-1*(c+d*x)))
@test integrate((c+d*x)^-3*SpecialFunctions.expint(1,a+b*x), x) == :(-1//2*d^-1*(c+d*x)^-2*SpecialFunctions.expint(1,a+b*x)+-1//2*b^2*d^-1*(b*c+-1*a*d)^-2*expintegralei(-1a+-1*b*x)+(1/2)*b^2*d^-1*(b*c+-1*a*d)^-2*exp(-1a+b*c*d^-1)*expintegralei(-1*b*d^-1*(c+d*x))+-1//2*b*d^-1*(c+d*x)^-1*(b*c+-1*a*d)^-1*exp(-1a+-1*b*x)+-1//2*b^2*d^-2*(b*c+-1*a*d)^-1*exp(-1a+b*c*d^-1)*expintegralei(-1*b*d^-1*(c+d*x)))
@test integrate((c+d*x)^-4*SpecialFunctions.expint(1,a+b*x), x) == :(-1//3*d^-1*(c+d*x)^-3*SpecialFunctions.expint(1,a+b*x)+-1//3*b^3*d^-1*(b*c+-1*a*d)^-3*expintegralei(-1a+-1*b*x)+-1//3*b^2*d^-1*(c+d*x)^-1*(b*c+-1*a*d)^-2*exp(-1a+-1*b*x)+-1//3*b^3*d^-2*(b*c+-1*a*d)^-2*exp(-1a+b*c*d^-1)*expintegralei(-1*b*d^-1*(c+d*x))+-1//6*b*d^-1*(c+d*x)^-2*(b*c+-1*a*d)^-1*exp(-1a+-1*b*x)+1//3*b^3*d^-1*(b*c+-1*a*d)^-3*exp(-1a+b*c*d^-1)*expintegralei(-1*b*d^-1*(c+d*x))+1//6*b^2*d^-2*(c+d*x)^-1*(b*c+-1*a*d)^-1*exp(-1a+-1*b*x)+1//6*b^3*d^-3*(b*c+-1*a*d)^-1*exp(-1a+b*c*d^-1)*expintegralei(-1*b*d^-1*(c+d*x)))
@test integrate((c+d*x)^3*SpecialFunctions.expint(2,a+b*x), x) == :(-1*b^-1*(c+d*x)^3*SpecialFunctions.expint(3,a+b*x)+-6*b^-4*d^3*SpecialFunctions.expint(6,a+b*x)+-6*b^-3*d^2*(c+d*x)*SpecialFunctions.expint(5,a+b*x)+-3*d*b^-2*(c+d*x)^2*SpecialFunctions.expint(4,a+b*x))
@test integrate((c+d*x)^2*SpecialFunctions.expint(2,a+b*x), x) == :(-1*b^-1*(c+d*x)^2*SpecialFunctions.expint(3,a+b*x)+-2*b^-3*d^2*SpecialFunctions.expint(5,a+b*x)+-2*d*b^-2*(c+d*x)*SpecialFunctions.expint(4,a+b*x))
@test integrate((c+d*x)*SpecialFunctions.expint(2,a+b*x), x) == :(-1*d*b^-2*SpecialFunctions.expint(4,a+b*x)+-1*b^-1*(c+d*x)*SpecialFunctions.expint(3,a+b*x))
@test integrate(SpecialFunctions.expint(2,a+b*x), x) == :(-1*b^-1*SpecialFunctions.expint(3,a+b*x))
@test integrate((c+d*x)^-3*SpecialFunctions.expint(2,a+b*x), x) == :(-1//2*d^-1*(c+d*x)^-2*SpecialFunctions.expint(2,a+b*x)+(1/2)*b*d^-2*(c+d*x)^-1*SpecialFunctions.expint(1,a+b*x)+(1/2)*b^2*d^-2*(b*c+-1*a*d)^-1*expintegralei(-1a+-1*b*x)+-1//2*b^2*d^-2*(b*c+-1*a*d)^-1*exp(-1a+b*c*d^-1)*expintegralei(-1*b*d^-1*(c+d*x)))
@test integrate((c+d*x)^-4*SpecialFunctions.expint(2,a+b*x), x) == :(-1//3*d^-1*(c+d*x)^-3*SpecialFunctions.expint(2,a+b*x)+1//6*b*d^-2*(c+d*x)^-2*SpecialFunctions.expint(1,a+b*x)+1//6*b^3*d^-2*(b*c+-1*a*d)^-2*expintegralei(-1a+-1*b*x)+-1//6*b^3*d^-2*(b*c+-1*a*d)^-2*exp(-1a+b*c*d^-1)*expintegralei(-1*b*d^-1*(c+d*x))+1//6*b^2*d^-2*(c+d*x)^-1*(b*c+-1*a*d)^-1*exp(-1a+-1*b*x)+1//6*b^3*d^-3*(b*c+-1*a*d)^-1*exp(-1a+b*c*d^-1)*expintegralei(-1*b*d^-1*(c+d*x)))
@test integrate((c+d*x)^3*SpecialFunctions.expint(3,a+b*x), x) == :(-1*b^-1*(c+d*x)^3*SpecialFunctions.expint(4,a+b*x)+-6*b^-4*d^3*SpecialFunctions.expint(7,a+b*x)+-6*b^-3*d^2*(c+d*x)*SpecialFunctions.expint(6,a+b*x)+-3*d*b^-2*(c+d*x)^2*SpecialFunctions.expint(5,a+b*x))
@test integrate((c+d*x)^2*SpecialFunctions.expint(3,a+b*x), x) == :(-1*b^-1*(c+d*x)^2*SpecialFunctions.expint(4,a+b*x)+-2*b^-3*d^2*SpecialFunctions.expint(6,a+b*x)+-2*d*b^-2*(c+d*x)*SpecialFunctions.expint(5,a+b*x))
@test integrate((c+d*x)*SpecialFunctions.expint(3,a+b*x), x) == :(-1*d*b^-2*SpecialFunctions.expint(5,a+b*x)+-1*b^-1*(c+d*x)*SpecialFunctions.expint(4,a+b*x))
@test integrate(SpecialFunctions.expint(3,a+b*x), x) == :(-1*b^-1*SpecialFunctions.expint(4,a+b*x))
@test integrate((c+d*x)^-4*SpecialFunctions.expint(3,a+b*x), x) == :(-1//3*d^-1*(c+d*x)^-3*SpecialFunctions.expint(3,a+b*x)+-1//6*b^2*d^-3*(c+d*x)^-1*SpecialFunctions.expint(1,a+b*x)+-1//6*b^3*d^-3*(b*c+-1*a*d)^-1*expintegralei(-1a+-1*b*x)+1//6*b*d^-2*(c+d*x)^-2*SpecialFunctions.expint(2,a+b*x)+1//6*b^3*d^-3*(b*c+-1*a*d)^-1*exp(-1a+b*c*d^-1)*expintegralei(-1*b*d^-1*(c+d*x)))
@test integrate((c+d*x)^4*SpecialFunctions.expint(-1,a+b*x), x) == :(-8*b^-5*d^4*exp(-1a+-1*b*x)+-1*b^-1*(a+b*x)^-1*(c+d*x)^4*exp(-1a+-1*b*x)+-8*b^-4*d^3*(c+d*x)*exp(-1a+-1*b*x)+-4*b^-5*d^2*(b*c+-1*a*d)^2*exp(-1a+-1*b*x)+-4*b^-5*d^3*(b*c+-1*a*d)*exp(-1a+-1*b*x)+-4*b^-3*d^2*(c+d*x)^2*exp(-1a+-1*b*x)+4*d*b^-5*(b*c+-1*a*d)^3*expintegralei(-1a+-1*b*x)+-4*b^-4*d^2*(c+d*x)*(b*c+-1*a*d)*exp(-1a+-1*b*x))
@test integrate((c+d*x)^3*SpecialFunctions.expint(-1,a+b*x), x) == :(-3*b^-4*d^3*exp(-1a+-1*b*x)+-1*b^-1*(a+b*x)^-1*(c+d*x)^3*exp(-1a+-1*b*x)+-3*b^-4*d^2*(b*c+-1*a*d)*exp(-1a+-1*b*x)+-3*b^-3*d^2*(c+d*x)*exp(-1a+-1*b*x)+3*d*b^-4*(b*c+-1*a*d)^2*expintegralei(-1a+-1*b*x))
@test integrate((c+d*x)^2*SpecialFunctions.expint(-1,a+b*x), x) == :(-2*b^-3*d^2*exp(-1a+-1*b*x)+-1*b^-1*(a+b*x)^-1*(c+d*x)^2*exp(-1a+-1*b*x)+2*d*b^-3*(b*c+-1*a*d)*expintegralei(-1a+-1*b*x))
@test integrate((c+d*x)*SpecialFunctions.expint(-1,a+b*x), x) == :(d*b^-2*expintegralei(-1a+-1*b*x)+-1*b^-1*(a+b*x)^-1*(c+d*x)*exp(-1a+-1*b*x))
@test integrate(SpecialFunctions.expint(-1,a+b*x), x) == :(-1*b^-1*(a+b*x)^-1*exp(-1a+-1*b*x))
@test integrate((c+d*x)^-1*SpecialFunctions.expint(-1,a+b*x), x) == :(-1*d*(b*c+-1*a*d)^-2*expintegralei(-1a+-1*b*x)+-1*(b*c+-1*a*d)^-1*exp(-1a+b*c*d^-1)*expintegralei(-1*b*d^-1*(c+d*x))+d*(b*c+-1*a*d)^-2*exp(-1a+b*c*d^-1)*expintegralei(-1*b*d^-1*(c+d*x))+-1*b^-1*(a+b*x)^-1*(c+d*x)^-1*exp(-1a+-1*b*x)+-1*d*b^-1*(c+d*x)^-1*(b*c+-1*a*d)^-1*exp(-1a+-1*b*x))
@test integrate((c+d*x)^-2*SpecialFunctions.expint(-1,a+b*x), x) == :((c+d*x)^-1*(b*c+-1*a*d)^-1*exp(-1a+-1*b*x)+-1*b^-1*(a+b*x)^-1*(c+d*x)^-2*exp(-1a+-1*b*x)+-2*b*d*(b*c+-1*a*d)^-3*expintegralei(-1a+-1*b*x)+-2*b*(b*c+-1*a*d)^-2*exp(-1a+b*c*d^-1)*expintegralei(-1*b*d^-1*(c+d*x))+-2*d*(c+d*x)^-1*(b*c+-1*a*d)^-2*exp(-1a+-1*b*x)+b*d^-1*(b*c+-1*a*d)^-1*exp(-1a+b*c*d^-1)*expintegralei(-1*b*d^-1*(c+d*x))+-1*d*b^-1*(c+d*x)^-2*(b*c+-1*a*d)^-1*exp(-1a+-1*b*x)+2*b*d*(b*c+-1*a*d)^-3*exp(-1a+b*c*d^-1)*expintegralei(-1*b*d^-1*(c+d*x)))
@test integrate((c+d*x)^-3*SpecialFunctions.expint(-1,a+b*x), x) == :((c+d*x)^-2*(-2*a*d+2*b*c)^-1*exp(-1a+-1*b*x)+-1*b^-1*(a+b*x)^-1*(c+d*x)^-3*exp(-1a+-1*b*x)+-3*d*b^2*(b*c+-1*a*d)^-4*expintegralei(-1a+-1*b*x)+-3*b^2*(b*c+-1*a*d)^-3*exp(-1a+b*c*d^-1)*expintegralei(-1*b*d^-1*(c+d*x))+-3//2*d*(c+d*x)^-2*(b*c+-1*a*d)^-2*exp(-1a+-1*b*x)+3//2*b*(c+d*x)^-1*(b*c+-1*a*d)^-2*exp(-1a+-1*b*x)+-1*d*b^-1*(c+d*x)^-3*(b*c+-1*a*d)^-1*exp(-1a+-1*b*x)+-3*b*d*(c+d*x)^-1*(b*c+-1*a*d)^-3*exp(-1a+-1*b*x)+3*d*b^2*(b*c+-1*a*d)^-4*exp(-1a+b*c*d^-1)*expintegralei(-1*b*d^-1*(c+d*x))+-1//2*b*d^-1*(c+d*x)^-1*(b*c+-1*a*d)^-1*exp(-1a+-1*b*x)+-1//2*b^2*d^-2*(b*c+-1*a*d)^-1*exp(-1a+b*c*d^-1)*expintegralei(-1*b*d^-1*(c+d*x))+3//2*b^2*d^-1*(b*c+-1*a*d)^-2*exp(-1a+b*c*d^-1)*expintegralei(-1*b*d^-1*(c+d*x)))
@test integrate((c+d*x)^4*SpecialFunctions.expint(-2,a+b*x), x) == :(-1*b^-1*(c+d*x)^4*SpecialFunctions.expint(-1,a+b*x)+-12*b^-5*d^4*exp(-1a+-1*b*x)+-12*b^-5*d^3*(b*c+-1*a*d)*exp(-1a+-1*b*x)+-12*b^-4*d^3*(c+d*x)*exp(-1a+-1*b*x)+12*b^-5*d^2*(b*c+-1*a*d)^2*expintegralei(-1a+-1*b*x)+-4*d*b^-2*(a+b*x)^-1*(c+d*x)^3*exp(-1a+-1*b*x))
@test integrate((c+d*x)^3*SpecialFunctions.expint(-2,a+b*x), x) == :(-1*b^-1*(c+d*x)^3*SpecialFunctions.expint(-1,a+b*x)+-6*b^-4*d^3*exp(-1a+-1*b*x)+6*b^-4*d^2*(b*c+-1*a*d)*expintegralei(-1a+-1*b*x)+-3*d*b^-2*(a+b*x)^-1*(c+d*x)^2*exp(-1a+-1*b*x))
@test integrate((c+d*x)^2*SpecialFunctions.expint(-2,a+b*x), x) == :(-1*b^-1*(c+d*x)^2*SpecialFunctions.expint(-1,a+b*x)+2*b^-3*d^2*expintegralei(-1a+-1*b*x)+-2*d*b^-2*(a+b*x)^-1*(c+d*x)*exp(-1a+-1*b*x))
@test integrate((c+d*x)*SpecialFunctions.expint(-2,a+b*x), x) == :(-1*b^-1*(c+d*x)*SpecialFunctions.expint(-1,a+b*x)+-1*d*b^-2*(a+b*x)^-1*exp(-1a+-1*b*x))
@test integrate(SpecialFunctions.expint(-2,a+b*x), x) == :(-1*b^-1*SpecialFunctions.expint(-1,a+b*x))
@test integrate((c+d*x)^-1*SpecialFunctions.expint(-2,a+b*x), x) == :(-1*b^-1*(c+d*x)^-1*SpecialFunctions.expint(-1,a+b*x)+-1*(b*c+-1*a*d)^-1*exp(-1a+b*c*d^-1)*expintegralei(-1*b*d^-1*(c+d*x))+2*d^2*(b*c+-1*a*d)^-3*expintegralei(-1a+-1*b*x)+-2*d^2*(b*c+-1*a*d)^-3*exp(-1a+b*c*d^-1)*expintegralei(-1*b*d^-1*(c+d*x))+2*d*(b*c+-1*a*d)^-2*exp(-1a+b*c*d^-1)*expintegralei(-1*b*d^-1*(c+d*x))+d*b^-2*(a+b*x)^-1*(c+d*x)^-2*exp(-1a+-1*b*x)+b^-2*d^2*(c+d*x)^-2*(b*c+-1*a*d)^-1*exp(-1a+-1*b*x)+-1*d*b^-1*(c+d*x)^-1*(b*c+-1*a*d)^-1*exp(-1a+-1*b*x)+2*b^-1*d^2*(c+d*x)^-1*(b*c+-1*a*d)^-2*exp(-1a+-1*b*x))
@test integrate((c+d*x)^-2*SpecialFunctions.expint(-2,a+b*x), x) == :((c+d*x)^-1*(b*c+-1*a*d)^-1*exp(-1a+-1*b*x)+-1*b^-1*(c+d*x)^-2*SpecialFunctions.expint(-1,a+b*x)+-3*b*(b*c+-1*a*d)^-2*exp(-1a+b*c*d^-1)*expintegralei(-1*b*d^-1*(c+d*x))+-3*d*(c+d*x)^-1*(b*c+-1*a*d)^-2*exp(-1a+-1*b*x)+6*b*d^2*(b*c+-1*a*d)^-4*expintegralei(-1a+-1*b*x)+6*d^2*(c+d*x)^-1*(b*c+-1*a*d)^-3*exp(-1a+-1*b*x)+b*d^-1*(b*c+-1*a*d)^-1*exp(-1a+b*c*d^-1)*expintegralei(-1*b*d^-1*(c+d*x))+-1*d*b^-1*(c+d*x)^-2*(b*c+-1*a*d)^-1*exp(-1a+-1*b*x)+-6*b*d^2*(b*c+-1*a*d)^-4*exp(-1a+b*c*d^-1)*expintegralei(-1*b*d^-1*(c+d*x))+2*d*b^-2*(a+b*x)^-1*(c+d*x)^-3*exp(-1a+-1*b*x)+2*b^-2*d^2*(c+d*x)^-3*(b*c+-1*a*d)^-1*exp(-1a+-1*b*x)+3*b^-1*d^2*(c+d*x)^-2*(b*c+-1*a*d)^-2*exp(-1a+-1*b*x)+6*b*d*(b*c+-1*a*d)^-3*exp(-1a+b*c*d^-1)*expintegralei(-1*b*d^-1*(c+d*x)))
@test integrate((c+d*x)^-3*SpecialFunctions.expint(-2,a+b*x), x) == :((c+d*x)^-2*(-2*a*d+2*b*c)^-1*exp(-1a+-1*b*x)+-1*b^-1*(c+d*x)^-3*SpecialFunctions.expint(-1,a+b*x)+-6*b^2*(b*c+-1*a*d)^-3*exp(-1a+b*c*d^-1)*expintegralei(-1*b*d^-1*(c+d*x))+-2*d*(c+d*x)^-2*(b*c+-1*a*d)^-2*exp(-1a+-1*b*x)+2*b*(c+d*x)^-1*(b*c+-1*a*d)^-2*exp(-1a+-1*b*x)+6*d^2*(c+d*x)^-2*(b*c+-1*a*d)^-3*exp(-1a+-1*b*x)+12*b^2*d^2*(b*c+-1*a*d)^-5*expintegralei(-1a+-1*b*x)+-1*d*b^-1*(c+d*x)^-3*(b*c+-1*a*d)^-1*exp(-1a+-1*b*x)+-12*b^2*d^2*(b*c+-1*a*d)^-5*exp(-1a+b*c*d^-1)*expintegralei(-1*b*d^-1*(c+d*x))+-6*b*d*(c+d*x)^-1*(b*c+-1*a*d)^-3*exp(-1a+-1*b*x)+2*b^2*d^-1*(b*c+-1*a*d)^-2*exp(-1a+b*c*d^-1)*expintegralei(-1*b*d^-1*(c+d*x))+3*d*b^-2*(a+b*x)^-1*(c+d*x)^-4*exp(-1a+-1*b*x)+3*b^-2*d^2*(c+d*x)^-4*(b*c+-1*a*d)^-1*exp(-1a+-1*b*x)+4*b^-1*d^2*(c+d*x)^-3*(b*c+-1*a*d)^-2*exp(-1a+-1*b*x)+12*b*d^2*(c+d*x)^-1*(b*c+-1*a*d)^-4*exp(-1a+-1*b*x)+12*d*b^2*(b*c+-1*a*d)^-4*exp(-1a+b*c*d^-1)*expintegralei(-1*b*d^-1*(c+d*x))+-1//2*b*d^-1*(c+d*x)^-1*(b*c+-1*a*d)^-1*exp(-1a+-1*b*x)+-1//2*b^2*d^-2*(b*c+-1*a*d)^-1*exp(-1a+b*c*d^-1)*expintegralei(-1*b*d^-1*(c+d*x)))
@test integrate((c+d*x)^4*SpecialFunctions.expint(-3,a+b*x), x) == :(-1*b^-1*(c+d*x)^4*SpecialFunctions.expint(-2,a+b*x)+-24*b^-5*d^4*exp(-1a+-1*b*x)+-4*d*b^-2*(c+d*x)^3*SpecialFunctions.expint(-1,a+b*x)+24*b^-5*d^3*(b*c+-1*a*d)*expintegralei(-1a+-1*b*x)+-12*b^-3*d^2*(a+b*x)^-1*(c+d*x)^2*exp(-1a+-1*b*x))
@test integrate((c+d*x)^3*SpecialFunctions.expint(-3,a+b*x), x) == :(-1*b^-1*(c+d*x)^3*SpecialFunctions.expint(-2,a+b*x)+6*b^-4*d^3*expintegralei(-1a+-1*b*x)+-3*d*b^-2*(c+d*x)^2*SpecialFunctions.expint(-1,a+b*x)+-6*b^-3*d^2*(a+b*x)^-1*(c+d*x)*exp(-1a+-1*b*x))
@test integrate((c+d*x)^2*SpecialFunctions.expint(-3,a+b*x), x) == :(-1*b^-1*(c+d*x)^2*SpecialFunctions.expint(-2,a+b*x)+-2*d*b^-2*(c+d*x)*SpecialFunctions.expint(-1,a+b*x)+-2*b^-3*d^2*(a+b*x)^-1*exp(-1a+-1*b*x))
@test integrate((c+d*x)*SpecialFunctions.expint(-3,a+b*x), x) == :(-1*d*b^-2*SpecialFunctions.expint(-1,a+b*x)+-1*b^-1*(c+d*x)*SpecialFunctions.expint(-2,a+b*x))
@test integrate(SpecialFunctions.expint(-3,a+b*x), x) == :(-1*b^-1*SpecialFunctions.expint(-2,a+b*x))
@test integrate((c+d*x)^-1*SpecialFunctions.expint(-3,a+b*x), x) == :(-1*b^-1*(c+d*x)^-1*SpecialFunctions.expint(-2,a+b*x)+-1*(b*c+-1*a*d)^-1*exp(-1a+b*c*d^-1)*expintegralei(-1*b*d^-1*(c+d*x))+-6*d^3*(b*c+-1*a*d)^-4*expintegralei(-1a+-1*b*x)+d*b^-2*(c+d*x)^-2*SpecialFunctions.expint(-1,a+b*x)+-6*d^2*(b*c+-1*a*d)^-3*exp(-1a+b*c*d^-1)*expintegralei(-1*b*d^-1*(c+d*x))+3*d*(b*c+-1*a*d)^-2*exp(-1a+b*c*d^-1)*expintegralei(-1*b*d^-1*(c+d*x))+6*d^3*(b*c+-1*a*d)^-4*exp(-1a+b*c*d^-1)*expintegralei(-1*b*d^-1*(c+d*x))+b^-2*d^2*(c+d*x)^-2*(b*c+-1*a*d)^-1*exp(-1a+-1*b*x)+-1*d*b^-1*(c+d*x)^-1*(b*c+-1*a*d)^-1*exp(-1a+-1*b*x)+-6*b^-1*d^3*(c+d*x)^-1*(b*c+-1*a*d)^-3*exp(-1a+-1*b*x)+-3*b^-2*d^3*(c+d*x)^-2*(b*c+-1*a*d)^-2*exp(-1a+-1*b*x)+-2*b^-3*d^2*(a+b*x)^-1*(c+d*x)^-3*exp(-1a+-1*b*x)+-2*b^-3*d^3*(c+d*x)^-3*(b*c+-1*a*d)^-1*exp(-1a+-1*b*x)+3*b^-1*d^2*(c+d*x)^-1*(b*c+-1*a*d)^-2*exp(-1a+-1*b*x))
@test integrate((c+d*x)^-2*SpecialFunctions.expint(-3,a+b*x), x) == :((c+d*x)^-1*(b*c+-1*a*d)^-1*exp(-1a+-1*b*x)+-1*b^-1*(c+d*x)^-2*SpecialFunctions.expint(-2,a+b*x)+-24*b*d^3*(b*c+-1*a*d)^-5*expintegralei(-1a+-1*b*x)+-24*d^3*(c+d*x)^-1*(b*c+-1*a*d)^-4*exp(-1a+-1*b*x)+-4*b*(b*c+-1*a*d)^-2*exp(-1a+b*c*d^-1)*expintegralei(-1*b*d^-1*(c+d*x))+-4*d*(c+d*x)^-1*(b*c+-1*a*d)^-2*exp(-1a+-1*b*x)+2*d*b^-2*(c+d*x)^-3*SpecialFunctions.expint(-1,a+b*x)+12*d^2*(c+d*x)^-1*(b*c+-1*a*d)^-3*exp(-1a+-1*b*x)+b*d^-1*(b*c+-1*a*d)^-1*exp(-1a+b*c*d^-1)*expintegralei(-1*b*d^-1*(c+d*x))+-1*d*b^-1*(c+d*x)^-2*(b*c+-1*a*d)^-1*exp(-1a+-1*b*x)+-24*b*d^2*(b*c+-1*a*d)^-4*exp(-1a+b*c*d^-1)*expintegralei(-1*b*d^-1*(c+d*x))+-12*b^-1*d^3*(c+d*x)^-2*(b*c+-1*a*d)^-3*exp(-1a+-1*b*x)+-8*b^-2*d^3*(c+d*x)^-3*(b*c+-1*a*d)^-2*exp(-1a+-1*b*x)+-6*b^-3*d^2*(a+b*x)^-1*(c+d*x)^-4*exp(-1a+-1*b*x)+-6*b^-3*d^3*(c+d*x)^-4*(b*c+-1*a*d)^-1*exp(-1a+-1*b*x)+2*b^-2*d^2*(c+d*x)^-3*(b*c+-1*a*d)^-1*exp(-1a+-1*b*x)+4*b^-1*d^2*(c+d*x)^-2*(b*c+-1*a*d)^-2*exp(-1a+-1*b*x)+12*b*d*(b*c+-1*a*d)^-3*exp(-1a+b*c*d^-1)*expintegralei(-1*b*d^-1*(c+d*x))+24*b*d^3*(b*c+-1*a*d)^-5*exp(-1a+b*c*d^-1)*expintegralei(-1*b*d^-1*(c+d*x)))
@test integrate((c+d*x)^3*SpecialFunctions.expint(n,a+b*x), x) == :(-1*b^-1*(c+d*x)^3*SpecialFunctions.expint(1+n,a+b*x)+-6*b^-4*d^3*SpecialFunctions.expint(4+n,a+b*x)+-6*b^-3*d^2*(c+d*x)*SpecialFunctions.expint(3+n,a+b*x)+-3*d*b^-2*(c+d*x)^2*SpecialFunctions.expint(2+n,a+b*x))
@test integrate((c+d*x)^2*SpecialFunctions.expint(n,a+b*x), x) == :(-1*b^-1*(c+d*x)^2*SpecialFunctions.expint(1+n,a+b*x)+-2*b^-3*d^2*SpecialFunctions.expint(3+n,a+b*x)+-2*d*b^-2*(c+d*x)*SpecialFunctions.expint(2+n,a+b*x))
@test integrate((c+d*x)*SpecialFunctions.expint(n,a+b*x), x) == :(-1*d*b^-2*SpecialFunctions.expint(2+n,a+b*x)+-1*b^-1*(c+d*x)*SpecialFunctions.expint(1+n,a+b*x))
@test integrate(SpecialFunctions.expint(n,a+b*x), x) == :(-1*b^-1*SpecialFunctions.expint(1+n,a+b*x))
@test integrate(x^3*expintegralei(b*x), x) == :(1//4*x^4*expintegralei(b*x)+3//2*b^-4*exp(b*x)+-3//2*x*b^-3*exp(b*x)+-1//4*b^-1*x^3*exp(b*x)+3//4*b^-2*x^2*exp(b*x))
@test integrate(x^2*expintegralei(b*x), x) == :(-2//3*b^-3*exp(b*x)+1//3*x^3*expintegralei(b*x)+-1//3*b^-1*x^2*exp(b*x)+2//3*x*b^-2*exp(b*x))
@test integrate(x*expintegralei(b*x), x) == :((1/2)*b^-2*exp(b*x)+(1/2)*x^2*expintegralei(b*x)+-1//2*x*b^-1*exp(b*x))
@test integrate(expintegralei(b*x), x) == :(x*expintegralei(b*x)+-1*b^-1*exp(b*x))
@test integrate(x^-2*expintegralei(b*x), x) == :(b*expintegralei(b*x)+-1*x^-1*exp(b*x)+-1*x^-1*expintegralei(b*x))
@test integrate(x^-3*expintegralei(b*x), x) == :(-1//2*x^-2*expintegralei(b*x)+-1//4*x^-2*exp(b*x)+1//4*b^2*expintegralei(b*x)+-1//4*b*x^-1*exp(b*x))
@test integrate(x^-4*expintegralei(b*x), x) == :(-1//3*x^-3*expintegralei(b*x)+-1//9*x^-3*exp(b*x)+1//18*b^3*expintegralei(b*x)+-1//18*b*x^-2*exp(b*x)+-1//18*b^2*x^-1*exp(b*x))
@test integrate(x^2*expintegralei(b*x)^2, x) == :(-5//6*b^-3*exp(2*b*x)+1//3*x^3*expintegralei(b*x)^2+4//3*b^-3*expintegralei(2*b*x)+-4//3*b^-3*exp(b*x)*expintegralei(b*x)+1//3*x*b^-2*exp(2*b*x)+-2//3*b^-1*x^2*exp(b*x)*expintegralei(b*x)+4//3*x*b^-2*exp(b*x)*expintegralei(b*x))
@test integrate(x*expintegralei(b*x)^2, x) == :((1/2)*b^-2*exp(2*b*x)+(1/2)*x^2*expintegralei(b*x)^2+-1*b^-2*expintegralei(2*b*x)+b^-2*exp(b*x)*expintegralei(b*x)+-1*x*b^-1*exp(b*x)*expintegralei(b*x))
@test integrate(expintegralei(b*x)^2, x) == :(x*expintegralei(b*x)^2+2*b^-1*expintegralei(2*b*x)+-2*b^-1*exp(b*x)*expintegralei(b*x))
@test integrate((d*x)^m*expintegralei(b*x), x) == :(d^-1*(d*x)^(1+m)*(1+m)^-1*expintegralei(b*x)+-1*b^-1*(d*x)^m*(-1*b*x)^(-1m)*(1+m)^-1*SpecialFunctions.gamma(1+m,-1*b*x))
@test integrate(x^3*expintegralei(a+b*x), x) == :(1//4*x^4*expintegralei(a+b*x)+3//2*b^-4*exp(a+b*x)+(1/2)*a*b^-4*exp(a+b*x)+-3//2*x*b^-3*exp(a+b*x)+-1//4*a^4*b^-4*expintegralei(a+b*x)+-1//4*b^-1*x^3*exp(a+b*x)+1//4*a^2*b^-4*exp(a+b*x)+1//4*a^3*b^-4*exp(a+b*x)+3//4*b^-2*x^2*exp(a+b*x)+-1//2*a*x*b^-3*exp(a+b*x)+-1//4*x*a^2*b^-3*exp(a+b*x)+1//4*a*b^-2*x^2*exp(a+b*x))
@test integrate(x^2*expintegralei(a+b*x), x) == :(-2//3*b^-3*exp(a+b*x)+1//3*x^3*expintegralei(a+b*x)+-1//3*a*b^-3*exp(a+b*x)+-1//3*a^2*b^-3*exp(a+b*x)+-1//3*b^-1*x^2*exp(a+b*x)+1//3*a^3*b^-3*expintegralei(a+b*x)+2//3*x*b^-2*exp(a+b*x)+1//3*a*x*b^-2*exp(a+b*x))
@test integrate(x*expintegralei(a+b*x), x) == :((1/2)*b^-2*exp(a+b*x)+(1/2)*x^2*expintegralei(a+b*x)+(1/2)*a*b^-2*exp(a+b*x)+-1//2*x*b^-1*exp(a+b*x)+-1//2*a^2*b^-2*expintegralei(a+b*x))
@test integrate(expintegralei(a+b*x), x) == :(-1*b^-1*exp(a+b*x)+b^-1*(a+b*x)*expintegralei(a+b*x))
@test integrate(x^-2*expintegralei(a+b*x), x) == :(-1*x^-1*expintegralei(a+b*x)+-1*b*a^-1*expintegralei(a+b*x)+b*a^-1*exp(a)*expintegralei(b*x))
@test integrate(x^-3*expintegralei(a+b*x), x) == :(-1//2*x^-2*expintegralei(a+b*x)+(1/2)*a^-2*b^2*expintegralei(a+b*x)+(1/2)*a^-1*b^2*exp(a)*expintegralei(b*x)+-1//2*b*a^-1*x^-1*exp(a+b*x)+-1//2*a^-2*b^2*exp(a)*expintegralei(b*x))
@test integrate(x^2*expintegralei(a+b*x)^2, x) == :(-5//6*b^-3*exp(2a+2*b*x)+1//3*x^3*expintegralei(a+b*x)^2+4//3*b^-3*expintegralei(2a+2*b*x)+2*a*b^-3*expintegralei(2a+2*b*x)+2*a^2*b^-3*expintegralei(2a+2*b*x)+-4//3*b^-3*exp(a+b*x)*expintegralei(a+b*x)+-2//3*a*b^-3*exp(2a+2*b*x)+1//3*x*b^-2*exp(2a+2*b*x)+-2//3*a*b^-3*exp(a+b*x)*expintegralei(a+b*x)+-2//3*a^2*b^-3*exp(a+b*x)*expintegralei(a+b*x)+-2//3*b^-1*x^2*exp(a+b*x)*expintegralei(a+b*x)+-1//3*x*a^2*b^-2*expintegralei(a+b*x)^2+1//3*a^2*b^-3*expintegralei(a+b*x)^2*(a+b*x)+4//3*x*b^-2*exp(a+b*x)*expintegralei(a+b*x)+2//3*a*x*b^-2*exp(a+b*x)*expintegralei(a+b*x))
@test integrate(x*expintegralei(a+b*x)^2, x) == :((1/2)*b^-2*exp(2a+2*b*x)+(1/2)*x^2*expintegralei(a+b*x)^2+-1*b^-2*expintegralei(2a+2*b*x)+b^-2*exp(a+b*x)*expintegralei(a+b*x)+-2*a*b^-2*expintegralei(2a+2*b*x)+a*b^-2*exp(a+b*x)*expintegralei(a+b*x)+(1/2)*a*x*b^-1*expintegralei(a+b*x)^2+-1*x*b^-1*exp(a+b*x)*expintegralei(a+b*x)+-1//2*a*b^-2*expintegralei(a+b*x)^2*(a+b*x))
@test integrate(expintegralei(a+b*x)^2, x) == :(2*b^-1*expintegralei(2a+2*b*x)+b^-1*expintegralei(a+b*x)^2*(a+b*x)+-2*b^-1*exp(a+b*x)*expintegralei(a+b*x))
@test integrate(x^2*expintegralei(d*(a+b*log(c*x^n))), x) == :(1//3*x^3*expintegralei(d*(a+b*log(c*x^n)))+-1//3*x^3*(c*x^n)^(-3*n^-1)*exp(-3*a*b^-1*n^-1)*expintegralei(b^-1*n^-1*(3+b*d*n)*(a+b*log(c*x^n))))
@test integrate(x*expintegralei(d*(a+b*log(c*x^n))), x) == :((1/2)*x^2*expintegralei(d*(a+b*log(c*x^n)))+-1//2*x^2*(c*x^n)^(-2*n^-1)*exp(-2*a*b^-1*n^-1)*expintegralei(b^-1*n^-1*(2+b*d*n)*(a+b*log(c*x^n))))
@test integrate(expintegralei(d*(a+b*log(c*x^n))), x) == :(x*expintegralei(d*(a+b*log(c*x^n)))+-1*x*(c*x^n)^(-1*n^-1)*exp(-1*a*b^-1*n^-1)*expintegralei(b^-1*n^-1*(1+b*d*n)*(a+b*log(c*x^n))))
@test integrate(x^-1*expintegralei(d*(a+b*log(c*x^n))), x) == :(b^-1*n^-1*(a+b*log(c*x^n))*expintegralei(d*(a+b*log(c*x^n)))+-1*b^-1*d^-1*n^-1*(c*x^n)^(b*d)*exp(a*d))
@test integrate(x^-2*expintegralei(d*(a+b*log(c*x^n))), x) == :(-1*x^-1*expintegralei(d*(a+b*log(c*x^n)))+x^-1*(c*x^n)^(n^-1)*exp(a*b^-1*n^-1)*expintegralei(-1*b^-1*n^-1*(1+-1*b*d*n)*(a+b*log(c*x^n))))
@test integrate(x^-3*expintegralei(d*(a+b*log(c*x^n))), x) == :(-1//2*x^-2*expintegralei(d*(a+b*log(c*x^n)))+(1/2)*x^-2*(c*x^n)^(2*n^-1)*exp(2*a*b^-1*n^-1)*expintegralei(-1*b^-1*n^-1*(2+-1*b*d*n)*(a+b*log(c*x^n))))
@test integrate((e*x)^m*expintegralei(d*(a+b*log(c*x^n))), x) == :(e^-1*(e*x)^(1+m)*(1+m)^-1*expintegralei(d*(a+b*log(c*x^n)))+-1*e^-1*(c*x^n)^(-1*n^-1*(1+m))*(e*x)^(1+m)*(1+m)^-1*exp(-1*a*b^-1*n^-1*(1+m))*expintegralei(b^-1*n^-1*(a+b*log(c*x^n))*(1+m+b*d*n)))
@test integrate(x^-3*exp(b*x)*expintegralei(b*x), x) == :(2*b^2*expintegralei(2*b*x)+-1//4*x^-2*exp(2*b*x)+1//4*b^2*expintegralei(b*x)^2+-1*b*x^-1*exp(2*b*x)+-1//2*x^-2*exp(b*x)*expintegralei(b*x)+-1//2*b*x^-1*exp(b*x)*expintegralei(b*x))
@test integrate(x^-2*exp(b*x)*expintegralei(b*x), x) == :((1/2)*b*expintegralei(b*x)^2+-1*x^-1*exp(2*b*x)+2*b*expintegralei(2*b*x)+-1*x^-1*exp(b*x)*expintegralei(b*x))
@test integrate(x^-1*exp(b*x)*expintegralei(b*x), x) == :((1/2)*expintegralei(b*x)^2)
@test integrate(exp(b*x)*expintegralei(b*x), x) == :(-1*b^-1*expintegralei(2*b*x)+b^-1*exp(b*x)*expintegralei(b*x))
@test integrate(x*exp(b*x)*expintegralei(b*x), x) == :(b^-2*expintegralei(2*b*x)+-1//2*b^-2*exp(2*b*x)+-1*b^-2*exp(b*x)*expintegralei(b*x)+x*b^-1*exp(b*x)*expintegralei(b*x))
@test integrate(x^2*exp(b*x)*expintegralei(b*x), x) == :(-2*b^-3*expintegralei(2*b*x)+5//4*b^-3*exp(2*b*x)+2*b^-3*exp(b*x)*expintegralei(b*x)+-1//2*x*b^-2*exp(2*b*x)+b^-1*x^2*exp(b*x)*expintegralei(b*x)+-2*x*b^-2*exp(b*x)*expintegralei(b*x))
@test integrate(x^3*exp(b*x)*expintegralei(b*x), x) == :(-4*b^-4*exp(2*b*x)+6*b^-4*expintegralei(2*b*x)+-6*b^-4*exp(b*x)*expintegralei(b*x)+2*x*b^-3*exp(2*b*x)+-1//2*b^-2*x^2*exp(2*b*x)+b^-1*x^3*exp(b*x)*expintegralei(b*x)+-3*b^-2*x^2*exp(b*x)*expintegralei(b*x)+6*x*b^-3*exp(b*x)*expintegralei(b*x))
@test integrate(x^3*exp(a+b*x)*expintegralei(c+d*x), x) == :(-6*b^-4*exp(a+b*x)*expintegralei(c+d*x)+-6*b^-3*(b+d)^-1*exp(a+c+x*(b+d))+-3*b^-2*(b+d)^-2*exp(a+c+x*(b+d))+-2*b^-1*(b+d)^-3*exp(a+c+x*(b+d))+6*b^-4*exp(a+-1*b*c*d^-1)*expintegralei(d^-1*(b+d)*(c+d*x))+b^-1*x^3*exp(a+b*x)*expintegralei(c+d*x)+-1*b^-1*x^2*(b+d)^-1*exp(a+c+x*(b+d))+-3*b^-2*x^2*exp(a+b*x)*expintegralei(c+d*x)+2*x*b^-1*(b+d)^-2*exp(a+c+x*(b+d))+3*x*b^-2*(b+d)^-1*exp(a+c+x*(b+d))+6*x*b^-3*exp(a+b*x)*expintegralei(c+d*x)+b^-1*c^3*d^-3*exp(a+-1*b*c*d^-1)*expintegralei(d^-1*(b+d)*(c+d*x))+-1*c*b^-1*d^-1*(b+d)^-2*exp(a+c+x*(b+d))+-1*b^-1*c^2*d^-2*(b+d)^-1*exp(a+c+x*(b+d))+-3*c*b^-2*d^-1*(b+d)^-1*exp(a+c+x*(b+d))+3*b^-2*c^2*d^-2*exp(a+-1*b*c*d^-1)*expintegralei(d^-1*(b+d)*(c+d*x))+6*c*b^-3*d^-1*exp(a+-1*b*c*d^-1)*expintegralei(d^-1*(b+d)*(c+d*x))+c*x*b^-1*d^-1*(b+d)^-1*exp(a+c+x*(b+d)))
@test integrate(x^2*exp(a+b*x)*expintegralei(c+d*x), x) == :(b^-1*(b+d)^-2*exp(a+c+x*(b+d))+-2*b^-3*exp(a+-1*b*c*d^-1)*expintegralei(d^-1*(b+d)*(c+d*x))+2*b^-3*exp(a+b*x)*expintegralei(c+d*x)+2*b^-2*(b+d)^-1*exp(a+c+x*(b+d))+b^-1*x^2*exp(a+b*x)*expintegralei(c+d*x)+-1*x*b^-1*(b+d)^-1*exp(a+c+x*(b+d))+-2*x*b^-2*exp(a+b*x)*expintegralei(c+d*x)+c*b^-1*d^-1*(b+d)^-1*exp(a+c+x*(b+d))+-1*b^-1*c^2*d^-2*exp(a+-1*b*c*d^-1)*expintegralei(d^-1*(b+d)*(c+d*x))+-2*c*b^-2*d^-1*exp(a+-1*b*c*d^-1)*expintegralei(d^-1*(b+d)*(c+d*x)))
@test integrate(x*exp(a+b*x)*expintegralei(c+d*x), x) == :(b^-2*exp(a+-1*b*c*d^-1)*expintegralei(d^-1*(b+d)*(c+d*x))+-1*b^-1*(b+d)^-1*exp(a+c+x*(b+d))+-1*b^-2*exp(a+b*x)*expintegralei(c+d*x)+x*b^-1*exp(a+b*x)*expintegralei(c+d*x)+c*b^-1*d^-1*exp(a+-1*b*c*d^-1)*expintegralei(d^-1*(b+d)*(c+d*x)))
@test integrate(exp(a+b*x)*expintegralei(c+d*x), x) == :(b^-1*exp(a+b*x)*expintegralei(c+d*x)+-1*b^-1*exp(a+-1*b*c*d^-1)*expintegralei(d^-1*(b+d)*(c+d*x)))
@test integrate(x^2*logintegral(b*x), x) == :(-1//3*b^-3*expintegralei(4*log(b*x))+1//3*x^3*logintegral(b*x))
@test integrate(x*logintegral(b*x), x) == :((1/2)*x^2*logintegral(b*x)+-1//2*b^-2*expintegralei(3*log(b*x)))
@test integrate(logintegral(b*x), x) == :(x*logintegral(b*x)+-1*b^-1*expintegralei(2*log(b*x)))
@test integrate(x^-1*logintegral(b*x), x) == :(log(b*x)*logintegral(b*x)+-1*b*x)
@test integrate(x^-2*logintegral(b*x), x) == :(b*log(log(b*x))+-1*x^-1*logintegral(b*x))
@test integrate(x^-3*logintegral(b*x), x) == :((1/2)*b^2*expintegralei(-1*log(b*x))+-1//2*x^-2*logintegral(b*x))
@test integrate((d*x)^m*logintegral(b*x), x) == :(d^-1*(d*x)^(1+m)*(1+m)^-1*logintegral(b*x)+-1*b*d^-2*(b*x)^(-2+-1m)*(d*x)^(2+m)*(1+m)^-1*expintegralei((2+m)*log(b*x)))
@test integrate(x^2*logintegral(a+b*x), x) == :(-1//3*b^-3*expintegralei(4*log(a+b*x))+1//3*x^3*logintegral(a+b*x)+a*b^-3*expintegralei(3*log(a+b*x))+-1*a^2*b^-3*expintegralei(2*log(a+b*x))+1//3*a^3*b^-3*logintegral(a+b*x))
@test integrate(x*logintegral(a+b*x), x) == :((1/2)*x^2*logintegral(a+b*x)+-1//2*b^-2*expintegralei(3*log(a+b*x))+a*b^-2*expintegralei(2*log(a+b*x))+-1//2*a^2*b^-2*logintegral(a+b*x))
@test integrate(logintegral(a+b*x), x) == :(-1*b^-1*expintegralei(2*log(a+b*x))+b^-1*(a+b*x)*logintegral(a+b*x))

ext: jl

max_stars_repo_path: test/8 Special functions/8.3 Exponential integral functions.jl


text:
""" [Ackerman2009](@cite) """
Dycoms_RF02_q_tot(::Type{FT}) where {FT} =
    z -> if z <= 795.0
        FT(9.45) / FT(1000.0)
    else
        (FT(5) - FT(3) * (FT(1) - exp(-(z - FT(795)) / FT(500)))) / FT(1000)
    end

""" [Ackerman2009](@cite) """
Dycoms_RF02_θ_liq_ice(::Type{FT}) where {FT} =
    z -> if z <= 795.0
        FT(288.3)
    else
        FT(295.0) + (z - FT(795))^FT(1.0 / 3.0)
    end

""" [Ackerman2009](@cite) """
function Dycoms_RF02_u(::Type{FT}) where {FT}
    return z -> FT(3) + FT(4.3) * z / FT(1000) - FT(5)
end

""" [Ackerman2009](@cite) """
function Dycoms_RF02_v(::Type{FT}) where {FT}
    return z -> FT(-9) + FT(5.6) * z / FT(1000) - FT(-5.5)
end

""" [Ackerman2009](@cite) """
Dycoms_RF02_tke(::Type{FT}) where {FT} =
    z -> if z <= 795.0
        FT(1) - z / FT(1000)
    else
        FT(0)
    end


text:
#
# Correctness Tests
#

using Base.Test
using DataArrays
using Compat

my_tests = ["abstractarray.jl",
            "booleans.jl",
            "constructors.jl",
            "containers.jl",
            "conversions.jl",
            "data.jl",
            "dataarray.jl",
            "datamatrix.jl",
            "indexing.jl",
            "linalg.jl",
            "nas.jl",
            "operators.jl",
            "reduce.jl",
            "reducedim.jl",
            "broadcast.jl",
            "padding.jl",
            "pooleddataarray.jl",
            "extras.jl",
            "sort.jl",
            "statistics.jl",
            "literals.jl",
            "newtests/dataarray.jl",
            "newtests/datamatrix.jl",
            "newtests/datavector.jl"]

println("Running tests:")

for my_test in my_tests
    @printf " * %s\n" my_test
    include(my_test)
end


text:
push!(LOAD_PATH,"../src/")

using TheDancer
using Base.Test

N, Edisorder, Jdisorder, modelJ, B, dipolestrength, r = init!(3, 0.0, 0.0, 0.001, 298, 20)

"""
Test to compare old and new methods of calculating the dipole response.

Test will pass if the two methods are equal to within 1% of each other, else
test will fail.
"""
function New(dipoles, density, dampening)
    Field = FieldFromDensity(density)
    dipoles = UpdateDipole(Field, dipoles, dampening)
end

function main()
    S,E,H,psi,dipoles=prepare_model()
    density = abs.(psi.^2)

    dipoles_new = New(dipoles, density, 0.5)
    dipoles_old = dipoles_by_relaxation(dipoles, density, 0.5)

    diff = sum(abs.(dipoles_old - dipoles_new))
    percentage_diff = diff/sum(dipoles_new)

    ε = 1E-2
    @test  0 ≈ percentage_diff atol = ε

end

main()


text:
#Read/write values with type info to HDF5.
#-------------------------------------------------------------------------------

#=
Useful when reading in values that might be of unknown type. Thus, they are written
to .hdf5 with metadata identifying the data type.

Note that many base types are alreay understood by the hdf5 format, so no need to
write metadata.

# Interface to be implemented

 - _write_typed(grp::HDF5.Group, name::String, v::NEWTYPE): Writes out type data
with type name in attributes.
 - __read(::Type{NEWTYPE}, ds::HDF5.Dataset): Read in data from NEWTYPE (leaf datatypes)
 - __read(::Type{NEWTYPE}, ds::HDF5.Group): Read in data from NEWTYPE (complex structures)

## NOTE:
 - Must register type with MAP_STR2TYPE
=#

#==Useful constants
===============================================================================#

const NOTHING_STRING = "nothing"

#Data types that have read/write functions implemented:
#NOTE: limit supported types so that not all concrete types need separate read/write functions
const MAP_STR2TYPE = Dict{String, Type}(
	"Nothing" => Nothing,
	"Symbol" => Symbol,
	"RGBA" => Colorant, #Write out any Colorant to an #RRGGBBAA string

	#Higher-level types
	"Extents1D" => Extents1D,
	"Axis" => Axis,
	"FoldedAxis" => FoldedAxis,
	"GridCartesian" => GridCartesian,
	"LineAttributes" => LineAttributes,
	"GlyphAttributes" => GlyphAttributes,
	"TextAnnotation" => TextAnnotation,
	"HVMarker" => HVMarker,
	"Waveform" => Waveform,
	"YStrip" => YStrip,
	"DataF1" => DataF1,
	"DataHR" => DataHR,
	"DataRS" => DataRS,
)
const MAP_TYPE2STR = Dict{Type, String}(v=>k for (k,v) in MAP_STR2TYPE)

#Types that already have built-in HDF5 support (just write out natively):
const Types_HDF5Support = Union{Number, String}

#==Type declarations
===============================================================================#
#HDF5 reader will auto-detect type correctly:
struct Type_HDF5AutoDetect; end

#==Read/write type name in attributes
===============================================================================#
function _write_datatype_attr(ds::Union{HDF5.Group, HDF5.Dataset}, ::Type{T}) where T
	typestr = MAP_TYPE2STR[T]
	HDF5.attributes(ds)["TYPE"] = typestr
end

function _read_datatype_attr(ds::Union{HDF5.Group, HDF5.Dataset})
	if !haskey(HDF5.attributes(ds), "TYPE")
		return Type_HDF5AutoDetect
	end

	typestr = HDF5.read(HDF5.attributes(ds)["TYPE"])
	return MAP_STR2TYPE[typestr]
end

#==_write_typed(): Simple (leaf) datatypes. (Labels with type name.)
===============================================================================#
#Default behaviour: Assumes value is supported by HDF5 format
function _write_typed(grp::HDF5.Group, name::String, v::T) where T<:Types_HDF5Support
	grp[name] = v
	return #No need to _write_datatype_attr
end
function _write_typed(grp::HDF5.Group, name::String, v::Nothing)
	grp[name] = NOTHING_STRING #Redundancy check/easier to read HDF5 file
	_write_datatype_attr(grp[name], Nothing)
end
function _write_typed(grp::HDF5.Group, name::String, v::Symbol)
	grp[name] = String(v)
	_write_datatype_attr(grp[name], Symbol)
end
function _write_typed(grp::HDF5.Group, name::String, v::Colorant)
	vstr = "#" * Colors.hex(v, :RRGGBBAA)
	grp[name] = vstr
	_write_datatype_attr(grp[name], Colorant)
end

#==_write_typed(): More complex structures. (Labels with type name.)
===============================================================================#

#==__read(): Read data, but not type information.
===============================================================================#

#Types with built-in HDF5 support:
__read(::Type{Type_HDF5AutoDetect}, ds::HDF5.Dataset) = HDF5.read(ds)
function __read(::Type{Nothing}, ds::HDF5.Dataset)
	v = HDF5.read(ds)
	nstr = NOTHING_STRING
	if nstr != v
		path = HDF5.name(ds)
		throw(Meta.ParseError("__read(::Nothing, ::HDF5.Group): Read $v != $nstr:\n$path"))
	end
	return nothing
end
function __read(::Type{Symbol}, ds::HDF5.Dataset)
	v = HDF5.read(ds)
	return Symbol(v)
end
function __read(::Type{Colorant}, ds::HDF5.Dataset)
	v = HDF5.read(ds)
	return parse(Colorant, v)
end

#==Read/Write functions for high-level objects
===============================================================================#
function _writeobjectdata(grp::HDF5.Group, obj::T) where T
	for fname in fieldnames(T)
		v = getfield(obj, fname)
		_write_typed(grp, String(fname), v)
	end
	return
end

#When type is unknown, _read_typed() figures it out:
function _read_typed(grp::HDF5.Group, name::String)
	ds = grp[name]
	t = _read_datatype_attr(ds)
	return __read(t, ds)
end

#_readobjectdata: Needs object values to be written out with _write_typed():
function _readobjectdata(::Type{T}, grp::HDF5.Group) where T
	vlist = Array{Any}(nothing, fieldcount(T))
	for (i, fname) in enumerate(fieldnames(T))
		vlist[i] = _read_typed(grp, String(fname))
	end
	return T(vlist...)
end

#==Catch-alls (default behaviour)
===============================================================================#

__read(::Type{T}, grp::HDF5.Group) where T = _readobjectdata(T, grp)
function _write_typed(grp::HDF5.Group, name::String, v::T) where T #Catch-all
	#NOTE: need "name" parameter so that call signature is same with built-ins
	#      (a string is not dataset, not a group)
	objgrp = creategrp(grp, name)
	_writeobjectdata(objgrp, v)
	_write_datatype_attr(objgrp, T)
end

#Last line


text:
# using Distributed
# addprocs(5)

using Caesar, RoME, Profile
using RoMEPlotting, Gadfly, KernelDensityEstimatePlotting
using KernelDensityEstimate
using IncrementalInference
using DocStringExtensions
using DelimitedFiles
# using TransformUtils

include(joinpath(@__DIR__,"slamUtils.jl"))

function main(windowstart::Int, windowlen::Int,saswindow::Int,sasstart::Int)

poses = [Symbol("x$i") for i in 1:windowlen]
alldataframes = collect(windowstart:windowstart+windowlen);
sasdataframes = collect(windowstart+sasstart:1:(windowstart+sasstart+saswindow-1))
sasposes = [Symbol("x$i") for i in sasstart:sasstart+saswindow-1]

dataDir = joinpath("/media","data1","data","kayaks","20_gps_pos") #liljon data
cfgFile = joinpath("/media","data1","data","kayaks","SAS2D.yaml");
chirpFile = joinpath("/media","data1","data","kayaks","chirp250.txt");

# dataDir = joinpath(ENV["HOME"],"data","kayaks","20_gps_pos") #local copy
# cfgFile = joinpath(ENV["HOME"],"data","sas","SAS2D.yaml");
# chirpFile = joinpath(ENV["HOME"],"data","sas","chirp250.txt");
cfgd=loadConfigFile(cfgFile)

fg = initfg();
beacon = :l1
addVariable!(fg, beacon, Point2 )

posData = importdata_nav(alldataframes, datadir=dataDir);
navchecked, errorind = sanitycheck_nav(posData)
dposData = deepcopy(posData)
cumulativeDrift!(dposData,[0.0;0],[0.1,0.1])

waveformData = importdata_waveforms(sasdataframes,2, datadir=dataDir);
tcurrent = 1_000_000

for sym in poses
  addVariable!(fg, sym, DynPoint2(ut=tcurrent))
  tcurrent += 1_000_000
end

priors = [1]
#Priors
for i in priors
    xdotp = posData[i+1,1] - posData[i,1];
    ydotp = posData[i+1,2] - posData[i,2];
    dpμ = [posData[i,1];posData[i,2];xdotp;ydotp];
    dpσ = Matrix(Diagonal([0.1;0.1;0.1;0.1].^2))
    pp = DynPoint2VelocityPrior(MvNormal(dpμ,dpσ))
    addFactor!(fg, [poses[i];], pp, autoinit=false)
end

vps = 1:windowlen-1
for i in vps
    xdotp = dposData[i+1,1] - dposData[i,1];
    ydotp = dposData[i+1,2] - dposData[i,2];
    dpμ = [xdotp;ydotp;0;0];
    dpσ = Matrix(Diagonal([0.1;0.1;0.1;0.1].^2))
    vp = VelPoint2VelPoint2(MvNormal(dpμ,dpσ))
    addFactor!(fg, [poses[i];poses[i+1]], vp, autoinit=false)
    # IncrementalInference.doautoinit!(fg,[getVariable(fg,poses[i])])
end

# IncrementalInference.doautoinit!(fg,[getVariable(fg,poses[5])])
sas2d = prepareSAS2DFactor(saswindow, waveformData, rangemodel=:Correlator,
                           cfgd=cfgd, chirpFile=chirpFile)
addFactor!(fg, [beacon;sasposes], sas2d, autoinit=false)

# visualization tools for debugging
writeGraphPdf(fg,viewerapp="", engine="neato", filepath="/media/data1/data/kayaks/testfg.pdf")
# wipeBuildNewTree!(fg, drawpdf=true, show=false)

getSolverParams(fg).drawtree = true
#getSolverParams(fg).showtree = true

## solve the factor graph
tree, smt, hist = solveTree!(fg, recordcliqs=[:x3; :l1])

drawTree(tree, filepath="/media/data1/data/kayaks/testbt.pdf")

plk= [];

for sym in ls(fg) #plotting all syms labeled
    X1 = getKDEMean(getVertKDE(fg,sym))
    push!(plk, layer(x=[X1[1];],y=[X1[2];], Geom.point), Theme(default_color=colorant"green",point_size = 1.5pt,highlight_width = 0pt))
end
#
# for i in 1:pose_counter
#      X1 = getKDEMean(getVertKDE(fg,Symbol("x$i")))
#      push!(plk, layer(x=[X1[1];],y=[X1[2];], Geom.point), Theme(point_size = 1.5pt,highlight_width = 0pt))
# end
#
# push!(plk,layer(x=posData[:,1],y=posData[:,2], Geom.path, Theme(default_color=colorant"green")))
#
igt = [17.0499;1.7832];
push!(plk,layer(x=[igt[1];],y=[igt[2];], label=String["Beacon";],Geom.point,Geom.label(hide_overlaps=false), order=2, Theme(default_color=colorant"red")));
#
L1 = getVal(getVariable(fg, beacon))
K1 = plotKDEContour(getVertKDE(fg,:l1),xlbl="X (m)", ylbl="Y (m)",levels=5,layers=true);
push!(plk,K1...)
push!(plk,Gadfly.Theme(key_position = :none));
push!(plk, Coord.cartesian(xmin=-40, xmax=140, ymin=-150, ymax=75,fixed=true))

plkplot = Gadfly.plot(plk...); plkplot |> PDF("/media/data1/data/kayaks/testplot.pdf");

end


text:
export meanerror, classerror

"""
`meanerror(net::FFNNet,inputs::Vector{Vector{Float64}},
   outputs::Vector{Vector{Float64}}; cost::Function = quaderror)`

Mean error of the network in this sample (consisting of `inputs` and `outputs`). The error is measured using `cost` function.

### Arguments
* `net` (`FFNNet`): Feedforward Neural Network
* `inputs` (`Vector{Vector{Float64}}`): Vector containing input examples (each one is a vector). This have K elements, K being the number of examples in this dataset. Each input vector must have I elements, where I is the input size of `net`.
* `outputs` (`Vector{Vector{Float64}}`): Vector containing input examples (each one is a vector). This must also have K elements, K being the number of examples in this dataset. Each output vector must have O elements, where O is the size of the output layer of `net`.

### Keyword Arguments
* `cost` (`Function`, `quaderror` by default): Cost function to be used as error measure.
"""
function meanerror(net::FFNNet,
                   inputs::Vector{Vector{Float64}},
                   outputs::Vector{Vector{Float64}};
                   cost::Function = quaderror)

    total_error = 0.0
    for ex in eachindex(inputs)
        total_error += cost(propagate(net, inputs[ex]), outputs[ex])
    end
    return total_error/length(inputs)
end

"""
`classerror(net::FFNNet,
                 inputs::Vector{Vector{Float64}},
                 outputs::Vector{Vector{Float64}})`

Classification error of the network in this sample (consisting of `inputs` and
`outputs`). The error is measured counting the ammount of misclassified inputs.

### Arguments
* `net` (`FFNNet`): Feedforward Neural Network
* `inputs` (`Vector{Vector{Float64}}`): Vector containing input examples (each one is a vector). This have K elements, K being the number of examples in this dataset. Each input vector must have I elements, where I is the input size of `net`.
* `outputs` (`Vector{Vector{Float64}}`): Vector containing input examples (each one is a vector). This must also have K elements, K being the number of examples in this dataset. Each output vector must have O elements, where O is the size of the output layer of `net`.

### Details
This function assumes that a classification of a vector is the index of the maximum value inside that vector. For example, the vector `[0,1,0]` is classified as 2. This function takes the classification of the output of the network, and compares to the classification of the example's output, counting the number of mismatches.

Therefore, it's natural to use `softmax` activation on the last layer of `net` and *one-hot* encoding on the examples' outputs, but it's not mandatory.
"""
function classerror(net::FFNNet,
                    inputs::Vector{Vector{Float64}},
                    outputs::Vector{Vector{Float64}})

    return sum([indmax(outputs[i]) != indmax(propagate(net, inputs[i])) for i in 1:length(inputs)])
end


text:
"""
function to_hdf()
"""
function to_hdf(data::DataFrame, file_path::String, name::String)
    h5open(file_path, "w") do file
        g = create_group(file, name) # create a group
        for col in names(data)
            g[col] = data[:,col]
        end
    end
end

"""
function hdf_to_df()
"""
function hdf_to_df(file_path::String, name::String)

    data = DataFrame(h5read(file_path, name))
return data
end


text:
using ModelingToolkit, OrdinaryDiffEq, DiffEqParamEstim, MinimallyDisruptiveCurves, ForwardDiff, LinearAlgebra, Test

function make_model(input)
    @parameters t
    @parameters k, c, m
    D = Differential(t)
    @variables pos(t) vel(t)

    eqs = [D(pos) ~ vel,
        D(vel) ~ (-1 / m) * (c * vel + k * pos - input(t))
  ]

    ps = [k,c,m] .=>  [2.,1.,4.]
    ics = [pos, vel] .=> [1.,0.]
    od = ODESystem(eqs, t, first.(ics), first.(ps)
                                 , defaults=merge(Dict(first.(ics) .=> last.(ics))
                                 , Dict(first.(ps) .=> last.(ps))), name=:mass_spring
                                 )
    tspan = (0., 100.)
  # prob = ODEProblem(od, ics, tspan, ps)
    return od, ics, tspan, ps
end

od, ics, tspan, ps = make_model(t -> 0.)

"""
take a log transform of the od parameter space in two ways: at the level of the ODESystem and the level of the ODEProblem
"""

# to_fix = ["c2c","c2","c2a","c3c", "c1c", "a2"]
# tstrct_fix = fix_params(last.(ps), get_name_ids(ps, to_fix))

p0 = last.(ps)
tr = logabs_transform(p0)
log_od = transform_ODESystem(od, tr)
@test typeof(log_od) == ODESystem

prob1 = ODEProblem(od, [], tspan, [])

log_od2, log_ics2, log_ps2 = transform_problem(prob1, tr; unames=ModelingToolkit.get_states(od), pnames=ModelingToolkit.get_ps(od))

"""
check if the two manners of transforming the ODE system give the same output
"""

@test repr.(ModelingToolkit.get_ps(log_od)) == repr.(ModelingToolkit.get_ps(log_od2))

log_prob1 = ODEProblem(log_od, [], tspan, [])
log_prob2 = ODEProblem(log_od2, [], tspan, [])

sol1 = solve(log_prob1, Tsit5())
sol2 = solve(log_prob2, Tsit5())

@test sol1[end] == sol2[end]

"""
check if  log transforming the cost function on od gives the same result as an untransformed cost function on log_od
"""
tsteps = tspan[1]:1.:tspan[end]
nom_sol = solve(prob1, Tsit5())
lossf(sol) = sum([sum(abs2, el1 - el2) for (el1, el2) in zip(sol(tsteps).u, nom_sol(tsteps).u) ])

nom_cost = build_loss_objective(prob1, Tsit5(), lossf; mpg_autodiff=true)
log_cost = build_loss_objective(log_prob1, Tsit5(), lossf; mpg_autodiff=true)

@test nom_cost(p0) == log_cost(log.(p0))

tr_cost, newp0 = transform_cost(nom_cost, p0, tr)
@test tr_cost(newp0) == log_cost(log.(p0))
grad_holder = deepcopy(p0)
g2 = deepcopy(grad_holder)
"""
test that summing losses works
"""
ll = sum_losses([nom_cost,nom_cost], p0)
@test ll(p0 .+ 1., grad_holder) == 2nom_cost(p0 .+1, g2)
@test grad_holder == 2g2

"""
test gradients of cost functions are zero at minimum as a proxy for correctness of their gradients
"""

for el in (log_cost, tr_cost)
    el(newp0, grad_holder)
    @test norm(grad_holder) < 1e-3 # 0 gradient at minimum
end

"""
test that mdc curve evolves, and listens to mdc_callbacks
"""
H0 = ForwardDiff.hessian(tr_cost, newp0)
mom = 1.
span = (-1., 1.);

newdp0 = (eigen(H0)).vectors[:, 1]

eprob = MDCProblem(log_cost, newp0, newdp0, mom, span);

cb = [
    Verbose([CurveDistance(0.1:1:10), HamiltonianResidual(2.3:4:10)]),
    ParameterBounds([1,3], [-10.,-10.], [10.,10.])
    ]

@time mdc = evolve(eprob, Tsit5; mdc_callback=cb);

"""
check mdc works with mdc_callback vector of subtype T <: CallbackCallable, strict subtype.
"""

cb = [
    Verbose([CurveDistance(0.1:1:10), HamiltonianResidual(2.3:4:10)])
    ]

@time mdc = evolve(eprob, Tsit5; mdc_callback=cb);

"""
 test MDC works and gives reasonable output
"""

@test log_cost(mdc.sol[1][1:3]) < 1e-3
@test log_cost(mdc.sol[end][1:3]) < 1e-3

"""
 test injection loss works OK
"""
l2 = build_injection_loss(prob1, Tsit5(), tsteps)
@test l2(p0, grad_holder) == 0
@test norm(grad_holder) < 1e-5

"""
test fixing parameters (i.e. a transformation that changes the number of parameters) works ok
"""
to_fix = ["c"]
trf = fix_params(last.(ps), get_name_ids(ps, to_fix))
de = MinimallyDisruptiveCurves.transform_ODESystem(od, trf)
@test length(ModelingToolkit.get_ps(de)) == 2

"""
test jumpstarting works
"""
jprob = jumpstart(eprob, 1e-2, true)
mdcj = evolve(jprob, Tsit5; mdc_callback=cb);
@test log_cost(mdcj.sol[1][1:3]) < 1e-3


text:
# Base for implementing all algorithms
"""
    AgentBase

`AgentBase` is the base type for all MAB Algorithms

"""
abstract type AgentBase end

"""
    reset!( agent::AgentBase )

Resets the internal statistics of the bandit algorithm, except for the number of arms and other
algorithm specific parameters. `reset!()` will make `agent` as fresh as it was first created.
"""
reset!( agent::AgentBase ) = nothing

"""
    info_str( agent::AgentBase, latex::Bool = false )

Return a information string about the agent and it's parameters.  If flag `latex` is set to
`true`, then the returned string will be a compactable latex string.
"""
info_str( agent::AgentBase, latex::Bool = false ) = nothing

"""
    show( agent::AgentBase )

Prints information about agent into the specified I/O.
```julia-repl
```
"""
function Base.show( io::IO, ::MIME"text/plain", agent::AgentBase )
    print( io, @sprintf("Agent: %s",info_str(agent)) )
    for param in fieldnames(agent)
        print( @sprintf("\n    %-16s: ",param), getfield(agent,param) )
    end
end

# Functions to make instances of types
function call( ::Type{T}, params... ) where T <: AgentBase
    return T( params... )
end

"""
    make_agents_with_k( K::Int64, agent_list::Vector{} )

Returns a vector of agents specified in `agent_list` with `K` arms.
```julia-repl

```
"""
function make_agents_with_k( K::Int64, agent_list::Vector{} )
    agents = Vector{AgentBase}()
    for (agent,params) ∈ agent_list
        if params != nothing
            push!( agents, agent(K,params...) )
        else
            push!( agents, agent(K) )
        end
    end
    return agents
end

import Base: ==
"""
    ==

Can be used to compare two agents.
```julia-repl

```
"""
function ==( agent1::AgentBase, agent2::AgentBase )
    if typeof(agent1) != typeof(agent2)
        return false
    else
        for param in fieldnames(typeof(agent1))
            if getfield(agent1,param) != getfield(agent2,param)
                # info( "Failed at ", param )
                return false
            end
        end
        return true
    end
    false
end


text:
@ms include("solvers/greedy.jl")

function loop_greedy(args)
    instances = [
        "01",
        "02",
        "03",
        "04",
        "05",
        "06",
        "07",
        "08",
        "09",
        "10",
        "11",
        "12",
        "13"
    ]
    println("Instance n° & Cost & Time")
    for instance_name in instances
        instance_path = "data/$instance_name.alp"
        instance = Instance(instance_path)
        solution = Solution(instance)

        ms_start = ms()
        # initial_sort!(solution)
        greedy!(solution)
        ms_stop = ms()
        nb_sec = round(ms_stop - ms_start, digits=3)
        # println("$instance_name & $(solution.cost) & $nb_sec")
        println(solution.cost)
    end
end

loop_greedy(Args.args)


text:
using MCMCBenchmarks,CSV

ProjDir = @__DIR__
cd(ProjDir)

path = joinpath(pathof(MCMCBenchmarks),"../../Examples/Linear_Regression/results/")

folder = "2019_06_15T14_20_00/"
results = CSV.read(path*folder*"results.csv")

# Make plots
include("plotting.jl")


text:
# common Base functionality

function Base.repeat(a::AbstractGPUVecOrMat, m::Int, n::Int = 1)
    o, p = size(a, 1), size(a, 2)
    b = similar(a, o*m, p*n)
    if length(b) == 0
        return b
    end
    gpu_call(b, a, o, p, m, n; elements=n) do ctx, b, a, o, p, m, n
        j = linear_index(ctx)
        j > n && return
        d = (j - 1) * p + 1
        @inbounds for i in 1:m
            c = (i - 1) * o + 1
            for r in 1:p
                for k in 1:o
                    b[k - 1 + c, r - 1 + d] = a[k, r]
                end
            end
        end
        return
    end
    return b
end

function Base.repeat(a::AbstractGPUVector, m::Int)
    o = length(a)
    b = similar(a, o*m)
    if length(b) == 0
        return b
    end
    gpu_call(b, a, o, m; elements=m) do ctx, b, a, o, m
        i = linear_index(ctx)
        i > m && return
        c = (i - 1)*o + 1
        @inbounds for i in 1:o
            b[c + i - 1] = a[i]
        end
        return
    end
    return b
end

## PermutedDimsArrays

using Base: PermutedDimsArrays

# PermutedDimsArrays' custom copyto! doesn't know how to deal with GPU arrays
function PermutedDimsArrays._copy!(dest::PermutedDimsArray{T,N,<:Any,<:Any,<:AbstractGPUArray}, src) where {T,N}
    dest .= src
    dest
end

## concatenation

# hacky overloads to make simple vcat and hcat with numbers work as expected.
# we can't really make this work in general without Base providing
# a dispatch mechanism for output container type.
@inline Base.vcat(a::Number, b::AbstractGPUArray) =
    vcat(fill!(similar(b, typeof(a), (1,size(b)[2:end]...)), a), b)
@inline Base.hcat(a::Number, b::AbstractGPUArray) =
    hcat(fill!(similar(b, typeof(a), (size(b)[1:end-1]...,1)), a), b)


text:
module PSDMatrices

import Base: \, /, size, inv, copy, copy!, ==, show, similar, Matrix
using LinearAlgebra
import LinearAlgebra: det, logabsdet, diag

struct PSDMatrix{T,FactorType} <: AbstractMatrix{T}
    R::FactorType
    PSDMatrix(R::AbstractMatrix{T}) where {T} = new{T,typeof(R)}(R)
end

# Base overloads
Matrix(M::PSDMatrix) = M.R' * M.R
size(M::PSDMatrix) = (size(M.R, 2), size(M.R, 2))
inv(M::PSDMatrix) = PSDMatrix(inv(M.R'))
\(A::PSDMatrix, B::AbstractVecOrMat) = A.R \ (A.R' \ B)
/(B::AbstractVecOrMat, A::PSDMatrix) = B / A.R / A.R'
copy(M::PSDMatrix) = PSDMatrix(copy(M.R))
similar(M::PSDMatrix, element_type::Type=eltype(M)) = PSDMatrix(similar(M.R, element_type))
copy!(dst::PSDMatrix, src::PSDMatrix) = (copy!(dst.R, src.R); dst)
==(M1::PSDMatrix, M2::PSDMatrix) = M1.R == M2.R  # todo: same as isequal()?!
function show(io::IO, M::PSDMatrix)
    print(io, "$(size(M,1))x$(size(M,2)) $(typeof(M)); R=")
    show(io, M.R)
end
function show(io::IO, m::MIME"text/plain", M::PSDMatrix)
    println(io, "$(size(M,1))x$(size(M,2)) $(typeof(M)) ")
    print(io, " Right square root: R=")
    show(io, m, M.R)
end

# LinearAlgebra overloads

function det(M::PSDMatrix)
    confirm_factor_is_square(M)
    return det(M.R)^2
end

function logabsdet(M::PSDMatrix)
    confirm_factor_is_square(M)
    _logabsdet, _sign = logabsdet(M.R)
    return 2 * _logabsdet, _sign^2
end

function diag(M::PSDMatrix)
    out = similar(M.R, size(M.R, 2))
    sum!(abs2, out', M.R)
    return out
end

function confirm_factor_is_square(M::PSDMatrix)
    if size(M.R, 1) != size(M.R, 2)
        msg = (
            "The requested operation is not available for a PSDMatrix with a non-square factor." *
            "The factor of the received PSDMatrix has dimensions ($(size(M.R,1)), $(size(M.R,2))). " *
            "Try turning the PSDMatrix into a dense matrix first."
        )
        throw(MethodError(msg))
    end
end

# Custom functions

X_A_Xt(; A::PSDMatrix, X::AbstractMatrix) = PSDMatrix(A.R * X')
X_A_Xt(A::PSDMatrix, X::AbstractMatrix) = X_A_Xt(A=A, X=X)
function X_A_Xt!(out::PSDMatrix; A::PSDMatrix, X::AbstractMatrix)
    mul!(out.R, A.R, X')
    return out
end
X_A_Xt!(out::PSDMatrix, A::PSDMatrix, X::AbstractMatrix) = X_A_Xt!(out, A=A, X=X)

function add_cholesky(A::PSDMatrix, B::PSDMatrix)
    sum_dense = Matrix(A) + Matrix(B)
    factor = cholesky(sum_dense).U
    return PSDMatrix(factor)
end

function add_qr(A::PSDMatrix, B::PSDMatrix)
    stack = vcat(A.R, B.R)
    matrix = PSDMatrix(stack)
    return triangularize_factor(matrix)
end

function triangularize_factor(M::PSDMatrix)
    R = qr(M.R).R
    return PSDMatrix(UpperTriangular(R))
end

export PSDMatrix
export add_cholesky
export add_qr
export triangularize_factor
export X_A_Xt
export X_A_Xt!

end