Refactor functions #99

lbenet · 2017-03-19T22:46:02Z

This PR aims to build a more consistent API, in the sense that a function f of Taylor1 and functions of TaylorN should be treated as close as possible. Now, most of the functions have more or less the same structure (except maybe for the first coefficient), and call a related function which mutates the coefficients of the resulting function. The functions that mutate the coefficients (which are the recurrence relations on the coefficients) have nicer names and a more homogeneous structure across all function; note that not all functions have been implemented in this way. By doing this, we have extended with (essentially) no effort the inverse trigonometric functions #49 and the hyperbolic functions #64 to TaylorN.

I have also deleted some documentation of functions (e.g., exp(a) for a::Taylor1 or a::TaylorN), which do not really provided anything new but the extension to these types.

There are still some tests to be added, and the example notebooks should also be updated.

Among others, mulHomogCoeff -> mul! and divHomogCoeff -> div!, and internal changes.

Also introduce div! for TaylorN polynomials

fixorder is not needed anymore since getindex was fixed and returns zero of appropriate type whenever the requested coefficient is well inside the maximum order, but not defined as such.

…fficients The following functions for the homogeneous coefficients have been renamed: powerHomogCoef -> pow!, squareHomogCoef -> sqr!, sqrtHomogCoef -> sqrt!

This allows a more smooth approach for Taylor1 and TaylorN functions.

coveralls · 2017-03-19T23:01:27Z

Coverage decreased (-8.9%) to 89.465% when pulling f24f911 on refactor_functions into 487cb8b on master.

dpsanders · 2017-03-20T01:29:28Z

Overall this looks like a very nice improvement, thanks!

I think it can be made even better following some of my comments, trying to unify even more the functions for Taylor1 and TaylorN. I gave some concrete suggestions that I think will help achieve this to a large extent.

Basically, every time there is nearly-duplicated code, it often just needs an extra couple of definitions of extra functions to take care of the small differences, and reduce it to a single function definition.

lbenet · 2017-03-20T01:37:19Z

Did you submit your review observations? I haven't received anything and I am very much looking forward. In any case, thanks a lot for the feedback and the review!

I agree with you that there is indeed a pattern in the mutating functions, and that this should be exploited.

dpsanders · 2017-03-20T00:57:40Z

src/arithmetic.jl

+    @eval function *{T<:$W}(a::$T{T}, b::$T{T})
+        corder = max(a.order, b.order)
+        if $T == Taylor1
+            c = Taylor1(zero(a[1]), corder)


There must be a way of unifying this line and the next, it seems to me.
Is the definition of zero not right?

Just add a method zero(a, order) to make an object of that order filled with zeros, then this should work here and in the next line.

Good point! I'll implement it. This should definitely get rid of some of these kind of ifs.

dpsanders · 2017-03-20T01:01:17Z

src/arithmetic.jl

-    @inbounds for k = 0:a.order
-        coeffs[k+1] = mulHomogCoef(k, a.coeffs, b.coeffs)
+        @inbounds for ord = 0:corder
+            mul!(c, a, b, ord) # updates c[ord+1]


Is there not a nicer way of writing this?

Let me try to explain the "design principle". Essentially, we need the same mutating functions for the coefficients of either Taylor1 or TaylorN, which must include "what to update" (c), the info to update it (a and b here because this is a binary function), and the order of the coefficient to be updated. So, as far as I can see, this is some sort of minimum required. For mul! I tried to be as close as what we already had for HomogeneousPolynomial.

Do you have some suggestion on this?

dpsanders · 2017-03-20T01:02:28Z

src/arithmetic.jl

    order = a.order + b.order

-    order > get_order() && return HomogeneousPolynomial([zero(a[1])], get_order())
+    order > get_order() &&
+        return HomogeneousPolynomial([zero(a[1])], get_order())


This line and the next can be combined with something like min(order, get_order()).
Isn't there a better way of writing this zero?

I'm not sure if I understand your point here.

The case is the following: The multiplication of two HomogeneousPolynomial is another one whose order is the addition of the orders of the polynomials involved. Due to the fact that we have a maximum order for HomogeneousPolynomial (and TaylorN), whenever the resulting polynomial is of larger degree we return a zero of the largest degree. This is what is implemented here. We can certainly return zero with another order, but I'm not sure this helps at all...

I was thinking something like zero(HomogeneousPolynomial{T}, get_order(a))
for whatever is a suitable T. Maybe no better.

dpsanders · 2017-03-20T01:09:45Z

src/arithmetic.jl

+
+    @eval /{T<:NumberNotSeries}(a::$T{T}, b::T) = a * inv(b)
+
+    @eval function /{T<:NumberNotSeries,S<:NumberNotSeries}(a::$T{T}, b::S)


Any particular reason for splitting up the @eval begin...end into separate @eval here?

Not really...

dpsanders · 2017-03-20T01:11:34Z

src/arithmetic.jl

    end
+
+    @eval /{T<:NumberNotSeries}(a::$T, b::T) = a * inv(b)


Why is there this one and the one on line 369?

Good catch! Indeed, line 369 can be deleted.

dpsanders · 2017-03-20T01:20:07Z

src/functions.jl

+                c = $T( aux, order )
+            else
+                order = get_order()
+                @inbounds aux = exp( a[1][1] )


What is aux for here?

Same comments presumably apply to all the following functions.

You are right, this appears all over the place. aux is some sort of initial value which allows initializing the returned Taylor1 or TaylorN, and also defines the initial value for the recurrence relations.

dpsanders · 2017-03-20T01:22:03Z

src/functions.jl

+                @inbounds c = $T( cos(a[1]), order )
+            else
+                order = get_order()
+                @inbounds s = $T( sin(a[1][1]), order )


How about defining a function constant_term (or similar name) that does this correctly for Taylor1 and TaylorN.

i.e.

constant_term(a::Taylor1) = a[1] constant_term(a::TaylorN) = a[1][1]

I think this is the key comment that will eliminate code duplication here!

Very nice! I'll do it!

dpsanders · 2017-03-20T01:23:48Z

src/functions.jl

+# Recursive functions (homogeneous coefficients)
+for T in (:Taylor1, :TaylorN)
+    @eval begin
+        function exp!(c::$T, a::$T, k::Int)


How does exp! differ from exp?
The usual trick is to define exp! correctly, and then exp is just exp(x) = (y = copy(x); exp!(y)

exp(a) returns a polynomial of type typeof(a).

exp!(c, a, k) mutates c[k+1] with the k-th order coefficient of c = exp(a). The recursion relations are actually coded in exp!.

Maybe it would be useful to use https://github.com/alsam/OffsetArrays.jl
to do away with all these annoying +1 everywhere.

We should indeed consider it, but let's leave it (again) for another PR.

dpsanders · 2017-03-20T01:24:49Z

src/functions.jl

 doc"""
-    tanHomogCoef(kcoef, ac, coeffst2)
+    reverse(f)


Or can we overload inv?

What is the correct term for this function?

I guess inverse or inversion would be ok; it implements Lagrange series inversion.

Cc: @MikaelSlevinsky

That's why I was wondering if overloading inv makes sense here?

I prefer not to use inv since it usually refers to the inverse multiplication, i.e. 1/a.

I'm ok with whatever term you choose. For the case that is implemented (which is numerically stable), I believe inversion and reversion are equivalent.

Ok. So, I'll change it into inverse.

dpsanders · 2017-03-20T01:26:01Z

src/other_functions.jl

-end
+    ## real, imag, conj and ctranspose ##
+    for f in (:real, :imag, :conj)
+        fdot = Symbol(f, ".")


You don't seem to use fdot.

Good catch! Indeed, that was somehow problematic, and at the end $f.(...) worked fine. I'll erase it.

dpsanders · 2017-03-20T01:38:30Z

Sorry, done. Not used to the new interface!

lbenet · 2017-03-20T01:42:48Z

Thanks! I'll read it now!

Also, eliminate an unused function (fdot) in src/other_functions

lbenet · 2017-03-20T05:46:20Z

Thank you very much for the good suggestions and the whole review.

I just pushed few commits addressing your observations. There are some I haven't yet implemented (get_order be replaced by order), but I got rid of lots of code duplication and, actually, of the nasty if-block inside some @eval blocks. I introduced a couple of new functions (whose name we may indeed change).

Lots of tests are still to be implemented, but this is left for tomorrow.

Once again, thanks a lot for the review and ideas!

coveralls · 2017-03-20T06:03:34Z

Coverage decreased (-5.2%) to 93.145% when pulling 6f51e98 on refactor_functions into 487cb8b on master.

dpsanders

Very nice work! Just a few more small comments here, but feel free to merge!

dpsanders · 2017-03-20T13:30:29Z

src/functions.jl

+# Functions
+for T in (:Taylor1, :TaylorN)
+    @eval begin
+        ## exp ##


These comments can now be removed ;)

dpsanders · 2017-03-20T13:30:54Z

src/functions.jl

+        ## exp ##
+        function exp(a::$T)
+            order = max_order(a)
+            @inbounds aux = exp( constant_term(a) )


The @inbounds can be moved inside the constant_term function.

You are right; I thought I removed them from here...

dpsanders · 2017-03-20T13:34:04Z

src/functions.jl

+        function exp(a::$T)
+            order = max_order(a)
+            @inbounds aux = exp( constant_term(a) )
+            c = $T( aux, order )


Maybe get rid of aux everywhere and just write

c = $T( exp(constant_term(a)), order )

Can we think of a better name for constant_term?

I agree. There is at least one exception: tan. There, there is an auxiliary function which is the square, so we save one (expensive?) calculation...

Actually, there are other cases where it is convenient to have aux or a0... sorry for the randomness in naming this.

dpsanders · 2017-03-20T13:35:46Z

src/functions.jl

+                """))
+
+            order = max_order(a)
+            c = $T( asin(a0), order )


Do you mean the definition of r? I'm not sure...

I mean that this is in the acos function, but the definition of c has asin in it.

It's actually a nice trick that @blas-ko came up with: the recursion formula is essentially the same as for asin except for a sign and the constant term.

OK, thanks. Maybe this should be mentioned in a comment somewhere.

Good point.

dpsanders · 2017-03-20T13:37:18Z

src/functions.jl


-Return the Taylor expansion of $\sin(a)$, of order `a.order`, for
-`a::Taylor1` polynomial.
+            c[k+1] = zero_korder(c, k)


What does zero_korder do?

It avoids the nasty if $T == Taylor1 ... end block that were there. For Taylor1 it simply returns zero(a[1]), while for TaylorN it returns HomogeneousPolynomial(zero(a[1][1]),k.

dpsanders · 2017-03-20T13:52:06Z

src/power.jl

+        n == 1 && return a
+        n == 2 && return square(a)
+        n < 0 && return inv( a^(-n) )
+        return Base.power_by_squaring(a, n)


Maybe use some kind of binomial expansion instead?
Since after normalising we have something of the form [1 + f(t)]^n.
Hmm, sounds complicated.

Or just use exp(n*log(a)) here too -- have you compared efficiency?

The advantage of power_by_squaring is that, if you have a::Taylor1{Int}, the result is of the same type, as it would be expected. This is not preserved by exp(n*log(a)), which would return Taylor1{Float64}.

Long time ago I did some comparisons with what we have now and the same type of implementation for ^(a,r::Real); power_by_squaring was better, at least for small n`.

The same subtlety about the type of the result appears here. The nice thing of preserving the type is that you get exact answers for Fateman-like tests.

dpsanders · 2017-03-20T13:53:56Z

src/power.jl


-with $a$ a `Taylor1` polynomial.
+        if $T == Taylor1
+            @inbounds c[k+1] = zero(c[1])


Can these two lines be unified?

I forgot completly this file; I'll implement the same changes as in functions.jl here.

dpsanders · 2017-03-20T13:54:31Z

src/power.jl

+    @inbounds posTb = pos_table[c.order+1]
+
+    for na = 1:num_coeffs_a
+        @inbounds ca = a[na]


I think you can just put the @inbounds outside the for loop.

dpsanders · 2017-03-20T13:55:47Z

src/power.jl

+    for na = 1:num_coeffs_a
+        @inbounds ca = a[na]
+        ca == zero(T) && continue
+        @inbounds inda = index_table[a.order+1][na]


Try given a name to index_table[a.order+1] outside the loop. I expect that would be a bit more efficient. (Just like with posTb.)

Good point; done!

dpsanders · 2017-03-20T13:56:41Z

src/power.jl

+    return c
+end
+
+function sqrt(a::TaylorN)


Can this function now be unified with the previous one?

The catch here is that it is simple to factorize Taylor1 polynomials of the form a2 t^2 + a3 t^3 +... into something like t^2 (a2 + a3 t +...) so then taking the sqrt makes it simpler. The same trick doesn't apply in a straight forward way to TaylorN polynomials. (Yet, Taylor1{Taylor1{T}} behaves nicely again!)

…verse Some tests have been modified to reflect current functionality.

lbenet · 2017-03-20T16:04:45Z

I just pushed the changes we discussed. I'll add some more tests related to the new functions, perhaps some docstrings, before merging. Thanks again for the code review!

coveralls · 2017-03-20T16:34:36Z

Coverage decreased (-4.3%) to 93.998% when pulling 957d7be on refactor_functions into 487cb8b on master.

coveralls · 2017-03-20T16:34:36Z

Coverage decreased (-4.3%) to 93.998% when pulling 957d7be on refactor_functions into 487cb8b on master.

coveralls · 2017-03-20T16:34:36Z

Coverage decreased (-4.3%) to 93.998% when pulling 957d7be on refactor_functions into 487cb8b on master.

dpsanders · 2017-03-20T18:20:08Z

Looks good and tests pass! Shall I merge?

lbenet · 2017-03-20T18:21:36Z

Hold on... I'm checking performance, and it seems that it is worst.

lbenet · 2017-03-20T18:23:03Z

Memory allocation has jump ~10x... I'll check this later

coveralls · 2017-03-20T18:54:10Z

Coverage increased (+0.2%) to 98.502% when pulling 4c24a05 on refactor_functions into 487cb8b on master.

lbenet · 2017-03-23T02:07:15Z

I just pushed another commit which addresses some optimizations in memory allocation I observed using perf/fateman.jl. The problem was that expressions like c[i] += a[j] * b[k] for TaylorNs, where creating temporary HomogeneousPolynomials and therefore allocating memory. So I replaced them by mul!(c[i], a[j], b[k]) which does not allocate such temporaries.

This achieves much better memory allocation and performs better than current master for perf/fateman.jl.

coveralls · 2017-03-23T02:14:36Z

Coverage decreased (-0.3%) to 98.068% when pulling 9e77bd8 on refactor_functions into 487cb8b on master.

coveralls · 2017-03-23T02:14:36Z

Coverage decreased (-0.3%) to 98.068% when pulling 9e77bd8 on refactor_functions into 487cb8b on master.

coveralls · 2017-03-25T01:21:03Z

Coverage decreased (-0.7%) to 97.591% when pulling e98f3fd on refactor_functions into 487cb8b on master.

coveralls · 2017-03-25T02:23:08Z

Coverage decreased (-1.4%) to 96.924% when pulling e5122ca on refactor_functions into 487cb8b on master.

lbenet · 2017-03-25T03:53:31Z

I'm merging... What ever is left, will be the subject of a new PR.

Luis Benet added 10 commits March 16, 2017 19:06

More homegeneous structure for * and / for Taylor1 and TaylorN

9434d57

Among others, mulHomogCoeff -> mul! and divHomogCoeff -> div!, and internal changes.

Fix getindex for Taylor1 and TaylorN and get rid of fixorder in * and /

cd84749

Also introduce div! for TaylorN polynomials

Fix an evaluate method and remove fixorder

7e4cf7b

fixorder is not needed anymore since getindex was fixed and returns zero of appropriate type whenever the requested coefficient is well inside the maximum order, but not defined as such.

Use Base.power_by_squaring instead of defining more methods

873991e

Metaprogramming for * and mul!, and cleanup code

ceff4bf

Refactor ^, square and sqrt and related functions for homogeneous coe…

4b5e47a

…fficients The following functions for the homogeneous coefficients have been renamed: powerHomogCoef -> pow!, squareHomogCoef -> sqr!, sqrtHomogCoef -> sqrt!

Metaprogramming for the elementary functions and recurrence relations

d1ab529

This allows a more smooth approach for Taylor1 and TaylorN functions.

Some metaprogramming in other_functions.jl

b15d4df

Adapt docs/src/api.md to changes in the naming of functions

9ea7ac0

Fix two tests

f24f911

dpsanders requested changes Mar 20, 2017

View reviewed changes

Luis Benet added 4 commits March 19, 2017 19:50

Add add! and subst! for completeness, and tests

6ca4f07

Add methods for zero and one and metaprogramming for auxiliary functions

4b56d32

Add constant_term and max_order, and reduce repeated code

773ac9a

Add zero_korder and reduce code duplication

6f51e98

Also, eliminate an unused function (fdot) in src/other_functions

dpsanders approved these changes Mar 20, 2017

View reviewed changes

Address comments pointed out in code review, and rename reverse -> in…

957d7be

…verse Some tests have been modified to reflect current functionality.

Add tests related to new functions, some docstrings, and update api.md

4c24a05

Improve memory allocation mainly for TaylorN using mul!

9e77bd8

Two more small changes

e98f3fd

Add more tests for TaylorN

e5122ca

lbenet merged commit 0697deb into master Mar 25, 2017

lbenet deleted the refactor_functions branch March 29, 2017 18:40

PerezHz mentioned this pull request Mar 31, 2017

Notable performance slowdown after PR #99 in julia-0.6.0-pre.alpha #102

Closed


		@eval /{T<:NumberNotSeries}(a::$T{T}, b::T) = a * inv(b)

		@eval function /{T<:NumberNotSeries,S<:NumberNotSeries}(a::$T{T}, b::S)

Refactor functions #99

Refactor functions #99

Conversation

lbenet commented Mar 19, 2017

coveralls commented Mar 19, 2017 • edited Loading

dpsanders commented Mar 20, 2017

lbenet commented Mar 20, 2017

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

dpsanders commented Mar 20, 2017

lbenet commented Mar 20, 2017

lbenet commented Mar 20, 2017

coveralls commented Mar 20, 2017 • edited Loading

dpsanders left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

lbenet commented Mar 20, 2017

coveralls commented Mar 20, 2017

coveralls commented Mar 20, 2017

coveralls commented Mar 20, 2017 • edited Loading

dpsanders commented Mar 20, 2017

lbenet commented Mar 20, 2017

lbenet commented Mar 20, 2017

coveralls commented Mar 20, 2017 • edited Loading

coveralls commented Mar 19, 2017 •

edited

Loading

coveralls commented Mar 20, 2017 •

edited

Loading

coveralls commented Mar 20, 2017 •

edited

Loading

coveralls commented Mar 20, 2017 •

edited

Loading