Add getindex, setindex! methods for ::Colon

[PerezHz]

This PR is aimed at solving the issue mentioned in #100. It adds getindex and setindex! methods for ::Colon, along with corresponding tests. It also adds a new setindex! method for Taylor1, which allows the user to do:

julia> using TaylorSeries

julia> y = sin(Taylor1(16))
 1.0 t - 0.16666666666666666 t³ + 0.008333333333333333 t⁵ - 0.0001984126984126984 t⁷ + 2.7557319223985893e-6 t⁹ - 2.505210838544172e-8 t¹¹ + 1.6059043836821616e-10 t¹³ - 7.647163731819817e-13 t¹⁵ + 𝒪(t¹⁷)

julia> y[1:5] = rand(5)
5-element Array{Float64,1}:
 0.996061
 0.430487
 0.776109
 0.207677
 0.353624

julia> y
 0.9960609386261976 + 0.43048731930377904 t + 0.7761091106019427 t² + 0.20767709964338787 t³ + 0.3536237457090481 t⁴ + 0.008333333333333333 t⁵ - 0.0001984126984126984 t⁷ + 2.7557319223985893e-6 t⁹ - 2.505210838544172e-8 t¹¹ + 1.6059043836821616e-10 t¹³ - 7.647163731819817e-13 t¹⁵ + 𝒪(t¹⁷)

julia> y == sin(Taylor1(16))
false

Similar new methods for setindex! have been added for TaylorN and HomogeneousPolynomial. Tests are passing in travis. Feedback is welcome!

[dpsanders]

I don't really understand why you would want to be able to modify a Taylor expansion.

[PerezHz]

Well, for example, when solving an ODE using Taylor's method, it is useful to allocate a Taylor1 variable for each dependent variable, and then update its Taylor coefficients each timestep, instead of allocating a new Taylor1, etc...

[coveralls]

Coverage decreased (-0.3%) to 97.745% when pulling 06fe4a6 on PerezHz:jp/getsetindexcolon into 0697deb on JuliaDiff:master.

[lbenet]

Just to make it clear (which it wasn't for me), this PR allows the following to work:

julia> y[1:5] = rand(5)
5-element Array{Float64,1}:
 0.996061
 0.430487
 0.776109
 0.207677
 0.353624

This is pretty useful, as @PerezHz points out, for in-place evaluations (in TaylorIntegration.jl).

[lbenet]

Can you extend this even further, to allow dot-operations that avoid allocations?

Following with your example, it would be nice that the following would work:

julia> y = sin( Taylor1(16) );

julia> y[1:5] .= rand.()
5-element Array{Float64,1}:
 0.654379
 0.412568
 0.460999
 0.806723
 0.885951

julia> y[1:5]
5-element Array{Float64,1}:
  0.0     
  1.0     
  0.0     
 -0.166667
  0.0     

The last line shows that it didn't work. Yet, if I do the same with y.coeffs it does work:

julia> y.coeffs[1:5] .= rand.()
5-element SubArray{Float64,1,Array{Float64,1},Tuple{UnitRange{Int64}},true}:
 0.340261
 0.761345
 0.913692
 0.926201
 0.570581

julia> y[1:5]
5-element Array{Float64,1}:
 0.340261
 0.761345
 0.913692
 0.926201
 0.570581

[coveralls]

Coverage decreased (-0.7%) to 97.311% when pulling d692ded on PerezHz:jp/getsetindexcolon into 0697deb on JuliaDiff:master.

[PerezHz]

Ok, so the problem was that while the syntax y.coeffs[1:5] .= rand.() is equivalent to broadcast!(rand, view(y.coeffs, 1:5)), the syntax y[1:5] .= rand.() was equivalent to broadcast!(rand, y[1:5]), which in turn called getindex(y,1:5). And the way getindex was implemented for UnitRange and Colon, it was creating a independent copy of y[1:5]. The called code for getindex in the above example given by @lbenet was:

getindex(a::Taylor1, n::UnitRange) = a.coeffs[n]

In order to fix this, what I did in the last commit that I just pushed was to modify the getindex methods using view:

getindex(a::Taylor1, n::UnitRange) = view(a.coeffs, n)

and analogously for Colon; I left the Int methods for getindex untouched since for that particular case this change seems to be unnecessary. I also did corresponding changes for HomogeneousPolynomial and TaylorN.

With this change, now everything should be working properly:

julia> y = sin(t)
 1.0 t - 0.16666666666666666 t³ + 0.008333333333333333 t⁵ - 0.0001984126984126984 t⁷ + 2.7557319223985893e-6 t⁹ - 2.505210838544172e-8 t¹¹ + 1.6059043836821616e-10 t¹³ - 7.647163731819817e-13 t¹⁵ + 𝒪(t¹⁷)

julia> y[1:5] .= rand.()
5-element SubArray{Float64,1,Array{Float64,1},Tuple{UnitRange{Int64}},true}:
 0.915906
 0.239288
 0.232181
 0.246147
 0.473498

julia> y
 0.9159064943913973 + 0.239287582828001 t + 0.23218121042262174 t² + 0.24614747149358918 t³ + 0.47349817383751724 t⁴ + 0.008333333333333333 t⁵ - 0.0001984126984126984 t⁷ + 2.7557319223985893e-6 t⁹ - 2.505210838544172e-8 t¹¹ + 1.6059043836821616e-10 t¹³ - 7.647163731819817e-13 t¹⁵ + 𝒪(t¹⁷)

julia> y[1:5]
5-element SubArray{Float64,1,Array{Float64,1},Tuple{UnitRange{Int64}},true}:
 0.915906
 0.239288
 0.232181
 0.246147
 0.473498

4 out of 4 travis tests are passing on my fork (PerezHz/TaylorSeries.jl); but just noticed that 1 out 4 travis tests failed on the JuliaDiff/TaylorSeries.jl repo (only julia 0.6 on OS X failed)

[coveralls]

Coverage decreased (-0.6%) to 97.398% when pulling 6f188d6 on PerezHz:jp/getsetindexcolon into 0697deb on JuliaDiff:master.

[lbenet]

I restarted the tests, and now they all pass.

Coverage indicates that there are few methods of setindex! that have no tests, which involve HomogeneousPolynomial or TaylorN. Can you add them?

Aside from this, LGTM!

[PerezHz]

I restarted the tests, and now they all pass.

Thanks!

Coverage indicates that there are few methods of setindex! that have no tests, which involve HomogeneousPolynomial or TaylorN. Can you add them?

Sure, actually I'm working on it already!

[coveralls]

Coverage decreased (-0.5%) to 97.572% when pulling 346cc11 on PerezHz:jp/getsetindexcolon into 0697deb on JuliaDiff:master.

[coveralls]

Coverage decreased (-0.4%) to 97.658% when pulling 57e8ffb on PerezHz:jp/getsetindexcolon into 0697deb on JuliaDiff:master.

[coveralls]

Coverage increased (+0.05%) to 98.092% when pulling 4cfb124 on PerezHz:jp/getsetindexcolon into 0697deb on JuliaDiff:master.

[coveralls]

Coverage increased (+0.3%) to 98.352% when pulling c2fc172 on PerezHz:jp/getsetindexcolon into 0697deb on JuliaDiff:master.

[coveralls]

Coverage increased (+0.3%) to 98.352% when pulling 037eb38 on PerezHz:jp/getsetindexcolon into 0697deb on JuliaDiff:master.

[PerezHz]

Just finished adding the missing tests for setindex!. While writing those tests, I came across some interesting behaviour on current master:

julia> using TaylorSeries;

julia> p = set_variables("x", numvars=2, order=3);

julia> x = sin(p[1]+p[2]*p[1])+p[2]^2-(1-p[1])^3
 - 1.0 + 4.0 x₁ - 3.0 x₁² + 1.0 x₁ x₂ + 1.0 x₂² + 0.8333333333333334 x₁³ + 𝒪(‖x‖⁴)

julia> x[1]
 - 1.0

julia> x.coeffs[1] = HomogeneousPolynomial(1)^2
 1.0 x₁²

julia> x[1]
 1.0 x₁²

julia> x.coeffs[:]  = HomogeneousPolynomial(1)^2
 1.0 x₁²

julia> x
 1.0 x₁² + 1.0 x₁² + 1.0 x₁² + 1.0 x₁² + 𝒪(‖x‖⁴)

I would think that this is some unwanted behaviour for TaylorN, in the sense that the 0-th order coefficient of a TaylorN should only accept HomogeneousPolynomials of order 0; the 1st order coefficient only HomogeneousPolynomials of order 1; etc... Otherwise, one could do "unsafe" operations with TaylorNs, in the sense that one can destroy their structure. Of course, the new setindex! methods for TaylorN which I have added in this PR allow to do similar things:

julia> using TaylorSeries

julia> p = set_variables("x", numvars=2, order=3)
2-element Array{TaylorSeries.TaylorN{Float64},1}:
  1.0 x₁ + 𝒪(‖x‖⁴)
  1.0 x₂ + 𝒪(‖x‖⁴)

julia> x = sin(p[1]+p[2]*p[1])+p[2]^2-(1-p[1])^3
 - 1.0 + 4.0 x₁ - 3.0 x₁² + 1.0 x₁ x₂ + 1.0 x₂² + 0.8333333333333334 x₁³ + 𝒪(‖x‖⁴)

julia> x[:] = HomogeneousPolynomial(1)^3
 1.0 x₁³

julia> x[1] = HomogeneousPolynomial(1)^2
 1.0 x₁²

julia> x[2:3] = HomogeneousPolynomial(2)^2
 1.0 x₂²

julia> x
 1.0 x₁² + 1.0 x₂² + 1.0 x₂² + 1.0 x₁³ + 𝒪(‖x‖⁴)

Has this been discussed before?

[lbenet]

LGTM. Merging.

[lbenet]

Thanks!