These are chat archives for symengine/symengine

18th
Dec 2016
Isuru Fernando
@isuruf
Dec 18 2016 03:35
@ChrisRackauckas, sure, that would be great. We'll have to come up with an idea for GSoC that can span 3 months. Ideas I have are like symengine/SymEngine.jl#53 which is not suitable for 3 months of work.
Christopher Rackauckas
@ChrisRackauckas
Dec 18 2016 04:03
yeah, that's definitely much quicker
Isuru Fernando
@isuruf
Dec 18 2016 07:21
@ChrisRackauckas, do you know why, v=[BigFloat(123.0)]; unsafe_load(pointer(v)) doesn't work?
Christopher Rackauckas
@ChrisRackauckas
Dec 18 2016 07:21
No, but @ScottPJones might
Isuru Fernando
@isuruf
Dec 18 2016 07:23
if the array is an immutable type, it works
This is probably because BigFloat is not an immutable type and the memory layout is different
Scott P. Jones
@ScottPJones
Dec 18 2016 11:11
Yes, BigFloat is not immutable, it has a pointer into C allocated memory.
Isuru Fernando
@isuruf
Dec 18 2016 11:12
So, if I need to create an array, I have to do it in a C library?
Scott P. Jones
@ScottPJones
Dec 18 2016 11:21
No, and I think this may be a bug in Julia.
I just tried this myself, and doing a unsafe_load on v got me a segfault.
Isuru Fernando
@isuruf
Dec 18 2016 11:22
Oh, thanks for the info. I'll try to find a workaround
Scott P. Jones
@ScottPJones
Dec 18 2016 11:22
I already have one ;-)
The issue is that the pointer to the vector really should be of Ptr{Ptr{BigFloat}}
Isuru Fernando
@isuruf
Dec 18 2016 11:23
Makes sense
Scott P. Jones
@ScottPJones
Dec 18 2016 11:24
(because each element of the vector is really a pointer to an instance of the type BigFloat)
reinterpret(Ptr{Ptr{BigFloat}}, v) will give you a valid pointer
and you can pass that to C
Isuru Fernando
@isuruf
Dec 18 2016 11:25
Thanks
Scott P. Jones
@ScottPJones
Dec 18 2016 11:26
Your welcome! Glad to help!
Fun how @ChrisRackauckas keeps pointing me at interesting little problems ;-)
Isuru Fernando
@isuruf
Dec 18 2016 11:27
:smile:
Scott P. Jones
@ScottPJones
Dec 18 2016 11:28
(of course, I'm in awe of what he's done/is doing with his DiffEq package(s) 🤓)
Scott P. Jones
@ScottPJones
Dec 18 2016 11:45
BTW, why does the description of SymEngine on GitHub point to http://sympy.org? (SymPy says that it's written completely in Python, so I don't get the connection)
Isuru Fernando
@isuruf
Dec 18 2016 11:46
@ScottPJones, SymEngine was named CSymPy before and it made sense earlier
SymEngine is a rewrite of core parts of SymPy in C
Scott P. Jones
@ScottPJones
Dec 18 2016 11:47
Ah! So, hopefully, much better speed than SymPy
Isuru Fernando
@isuruf
Dec 18 2016 11:48
yes, speed is our main goal
Scott P. Jones
@ScottPJones
Dec 18 2016 11:49
(having moved from C/C++ to Julia [after a brief flirtation with Python 3], I wonder how an implementation completely in Julia would perform ;-) )
(I understand that C++ would better serve supporting multiple languages, since you can make a library out of it)
Isuru Fernando
@isuruf
Dec 18 2016 11:52
https://github.com/jcrist/Juniper.jl is one implementation, but I think he was using the same structures as in SymPy.
SymEngine does things a little bit differently, but keeps a higher level API similar to SymPy.
Yes, having it in C++, means we now have wrappers to Python, Ruby and Julia
Scott P. Jones
@ScottPJones
Dec 18 2016 11:52
I was just introduced to SymPy last week (via another one of these "maybe Scott will know" questions, about using the Julia PyCall wrapper to it), and was very impressed. If this gives the same capabilities, but much faster, I'm super impressed ;-)
Isuru Fernando
@isuruf
Dec 18 2016 11:53
There are some speed comparisons at the end here, https://github.com/certik/scipy-2016-symengine-talk/blob/master/talk.pdf
Featurewise though, SymEngine is far behind SymPy and others
Fortunately, @ChrisRackauckas and others need only a little subset of SymPy.
Scott P. Jones
@ScottPJones
Dec 18 2016 11:55
Guess that just means more fun work in the future? 🤓
Isuru Fernando
@isuruf
Dec 18 2016 11:55
:smiley: yes
Scott P. Jones
@ScottPJones
Dec 18 2016 11:56
Is there a RoadMap page for SymEngine?
Isuru Fernando
@isuruf
Dec 18 2016 11:56
Not really. We just make it up as we go
Scott P. Jones
@ScottPJones
Dec 18 2016 11:56
Would be interesting to see the plans for improving SymEngine
Isuru Fernando
@isuruf
Dec 18 2016 11:57
We have this roadmap for using SymEngine in SymPy, https://github.com/symengine/symengine/wiki/SymPy-core-upgrade-to-SymEngine
Scott P. Jones
@ScottPJones
Dec 18 2016 11:57
Nice "meeting" you - I'm going back to bed! ;-) I'll keep lurking here, this is very interesting!
Isuru Fernando
@isuruf
Dec 18 2016 12:01
Good night. Thanks for your time again
Isuru Fernando
@isuruf
Dec 18 2016 15:29
@ChrisRackauckas, I have a WIP branch for symengine/SymEngine.jl#53 if you are interested.
Here's a session,
julia> using SymEngine

julia> @vars x y z
(x,y,z)

julia> erf(Basic(1.1))
0.880205069574082

julia> erf(x * 1.1)
erf(1.1*x)

julia> diff(erf(x * 1.1), x)
2.2*exp(-1.21*x^2)/sqrt(pi)
Scott P. Jones
@ScottPJones
Dec 18 2016 15:34
Fascinating stuff, benchmark results look fantastic. You mentioned Juniper.jl, but that hasn't been updated since it was placed there two years ago.
Isuru Fernando
@isuruf
Dec 18 2016 15:36
yes, that's the only native julia CAS I know
Scott P. Jones
@ScottPJones
Dec 18 2016 15:37
SymEngine.jl needs some minor updates (ASCIIString was deprecated in v0.5, I think it will be removed in v0.6)
Would you like a PR?
Isuru Fernando
@isuruf
Dec 18 2016 15:38
Sure
Here's a benchmark with Calculus.jl
julia> using SymEngine

julia> using Calculus

julia> function se_diff(expr, sym, n::Int)
           t = expr
           for i in 1:n
               t = diff(t, sym)
           end
       end
se_diff (generic function with 1 method)

julia> function calc_diff(expr, sym, n::Int)
           t = expr
           for i in 1:n
               t = differentiate(t, sym)
           end
       end
calc_diff (generic function with 1 method)

julia> function calc_diff2(expr, sym, n::Int)
          t = expr
          for i in 1:n
              t = simplify(differentiate(t, sym))
          end
       end
calc_diff2 (generic function with 1 method)

julia> @vars x y z
(x,y,z)

julia> @time calc_diff("erf(x)", :x, 8)
  1.431160 seconds (2.88 M allocations: 204.814 MB, 2.36% gc time)

julia> @time se_diff(erf(x), x, 8)
  0.002720 seconds (38 allocations: 1.313 KB)

julia> @time calc_diff2("erf(x)", :x, 8)
  1.536351 seconds (3.08 M allocations: 216.925 MB, 2.34% gc time)
Scott P. Jones
@ScottPJones
Dec 18 2016 15:46
IIUC, for types (?) in Python that SymEngine doesn't understand, it calls back into SymPy, is that more or less correct?
Isuru Fernando
@isuruf
Dec 18 2016 15:47
In python bindings that's correct for functions
In above, I've used Julia cfunctions for evaluating erf(x) which SymEngine doesn't understand
Isuru Fernando
@isuruf
Dec 18 2016 15:54
I just checked SymPy, and in this case, SymEngine and SymPy are equivalent
I mean diff(erf(x), x, 8)
Christopher Rackauckas
@ChrisRackauckas
Dec 18 2016 17:06

Fortunately, @ChrisRackauckas and others need only a little subset of SymPy.

Yeah, I only need the differentiation and linear algebra stuff, but Julia's generic linear algebra "just works" so it handles that.

yes, that's the only native julia CAS I know

There's Symata.jl. But it moreso uses other CAS software like SymPy as backends, and implements a lot of Mathematica-like pattern matching.

I put an issue on their repo for moving some backend stuff over to SymEngine

In above, I've used Julia cfunctions for evaluating erf(x) which SymEngine doesn't understand

Yes, I have this exact error in a tutorial notebook. Looks like I will have to make the error more complicated now :smile:

Scott P. Jones
@ScottPJones
Dec 18 2016 17:09
@isuruf, would SymEngine be able to handle arbitrary Julia types (like ArbFloats, for example), via that cfunction mechanism?
Isuru Fernando
@isuruf
Dec 18 2016 17:10
No, I've written it only for functions
Do you have a usecase?
In SymEngine, support for Arb is limited, but that'll change in the future.
Scott P. Jones
@ScottPJones
Dec 18 2016 17:15
No particular use case, just one of the great things about @ChrisRackauckas's work, is that (due both to Julia, and his being careful with his generic programming) it supports arbitrary types.
Christopher Rackauckas
@ChrisRackauckas
Dec 18 2016 17:16
I convert SymEngine results to Julia functions.
Scott P. Jones
@ScottPJones
Dec 18 2016 17:16
Ah! That sounds genius. could that be made part of the SymEngine.jl wrapper?
Christopher Rackauckas
@ChrisRackauckas
Dec 18 2016 17:17
I think lambdify if a simple function for it?
I modify it a bit more.
I think a function which goes further and creates an in-place function would be the real deal though.
Isuru Fernando
@isuruf
Dec 18 2016 17:19
@ChrisRackauckas, I don't understand what you mean by that last comment
Christopher Rackauckas
@ChrisRackauckas
Dec 18 2016 17:21

Like, here's a usage. I have this macro:

f = @ode_def LotkaVolterra begin
  dx = a*x - b*x*y
  dy = -c*y + d*x*y
end a=>1.5 b=>1 c=3 d=1

which interprets that differential equation and makes a function function f(t,u,du) equal to

begin  # console, line 2:
    du[1] = p.a * u[1] - p.b * u[1] * u[2] # console, line 3:
    du[2] = -3 * u[2] + 1 * u[1] * u[2]
    nothing
end

SymEngine.jl comes in because I then use it to create an expression for the Jacobian. It's an inplace expression which writes to J, and the f(Val{:jac},t,u,J) is automatically computed as:

begin 
    J[1,1] = p.a - p.b * u[2]
    J[1,2] = -(p.b) * u[1]
    J[2,1] = 1 * u[2]
    J[2,2] = -3 + u[1]
    nothing
end
Basically, if you have an array of SymEngine.Basic, it would be nice to have it write a function which results in that form in a mutating output, since that can get you quite a bit more speed if the function is repeatedly called
Isuru Fernando
@isuruf
Dec 18 2016 17:26
I still have trouble understanding. You have a jacobian J and you have to evaluate it for different u ?
Christopher Rackauckas
@ChrisRackauckas
Dec 18 2016 17:27
J is "an empty 2x2 array` or a cache array.
I want to evaluate the Jacobian of the system at time (t,u) into the cache J
instead of making a new matrix
(this case is actually a little more subtle since there's call overloading going on here, but that inplace updating setup should actually be very common)
Isuru Fernando
@isuruf
Dec 18 2016 17:29
do you have like a simple example?
Christopher Rackauckas
@ChrisRackauckas
Dec 18 2016 17:33
using SymEngine
x = symbols("x")
A = Array{SymEngine.Basic}(2,2)
A[1,1] = x
A[1,2] = x^2
A[2,1] = 2x
A[2,2] = 2x^2
# iA = inv(A) # Not working?
# I would like inplace_lambdify(A) to produce:
f! = function (out,x) 
  out[1,1] = x
  out[1,2] = x^2
  out[2,1] = 2x
  out[2,2] = 2x^2
  nothing
end
Isuru Fernando
@isuruf
Dec 18 2016 17:34
Thanks, I now understand it.
these xs are doubles?
Scott P. Jones
@ScottPJones
Dec 18 2016 17:36
They could be any type that has ^ and * operations
Christopher Rackauckas
@ChrisRackauckas
Dec 18 2016 17:36
yeah
since it's left open, it could be BigFloat, ArbFloat, etc.
Scott P. Jones
@ScottPJones
Dec 18 2016 17:37
(scarily [and a pun I dislike in Julia], they could even be strings!)
Christopher Rackauckas
@ChrisRackauckas
Dec 18 2016 17:37
oh yeah, in this case they can be strings :smile:
^ on strings is repeat, right?
Scott P. Jones
@ScottPJones
Dec 18 2016 17:38
(I do so wish * and ^ on strings would get deprecated!)
Isuru Fernando
@isuruf
Dec 18 2016 17:38
There's no * for strings right?
Christopher Rackauckas
@ChrisRackauckas
Dec 18 2016 17:38
concatenation
Scott P. Jones
@ScottPJones
Dec 18 2016 17:38
Unfortunately, so
Christopher Rackauckas
@ChrisRackauckas
Dec 18 2016 17:38
Yeah, so it's a fluke that this works :smile:
Scott P. Jones
@ScottPJones
Dec 18 2016 17:38
I've advocated ++ (i.e. like Haskell)
Isuru Fernando
@isuruf
Dec 18 2016 17:38
2"asd" doesn't work
Scott P. Jones
@ScottPJones
Dec 18 2016 17:38
for a totally generic cat operator.
Christopher Rackauckas
@ChrisRackauckas
Dec 18 2016 17:39
it probably reads that as a string macro?
MethodError: no method matching *(::Int64, ::String)
Closest candidates are:
  *(::Any, ::Any, !Matched::Any, !Matched::Any...) at operators.jl:287
  *{T<:Union{Int128,Int16,Int32,Int64,Int8,UInt128,UInt16,UInt32,UInt64,UInt8}}(::T<:Union{Int128,Int16,Int32,Int64,Int8,UInt128,UInt16,UInt32,UInt64,UInt8}, !Matched::T<:Union{Int128,Int16,Int32,Int64,Int8,UInt128,UInt16,UInt32,UInt64,UInt8}) at int.jl:33
  *(::Real, !Matched::Complex{Bool}) at complex.jl:215
  ...
 in include_string(::String, ::String) at .\loading.jl:478
 in eval(::Module, ::Any) at .\boot.jl:236
 in (::Atom.##67#70)() at C:\Users\Chris\.julia\v0.6\Atom\src\eval.jl:40
 in withpath(::Atom.##67#70, ::Void) at C:\Users\Chris\.julia\v0.6\CodeTools\src\utils.jl:30
 in withpath(::Function, ::Void) at C:\Users\Chris\.julia\v0.6\Atom\src\eval.jl:46
 in macro expansion at C:\Users\Chris\.julia\v0.6\Atom\src\eval.jl:109 [inlined]
 in (::Atom.##66#69)() at .\task.jl:60
Ahh it's missing a dispatch for that to work
f! = function {T}(out,x::T) 
  out[1,1] = x
  out[1,2] = x^2
  out[2,1] = T(2)*x
  out[2,2] = T(2)*x^2
end
is totally unnecessary but would work with strings then.
Isuru Fernando
@isuruf
Dec 18 2016 17:42
@ChrisRackauckas, inplace_lambdify should be implemented. In Python, we've implemented both and have a unified interface
Scott P. Jones
@ScottPJones
Dec 18 2016 17:42
Funny, because I'd thought *(a,b) was equivalent to string(a,b), if either a or b were an AbstractString
but apparently not. string(2,"foo") will give you "2foo", because it stringifies all arguments.
Christopher Rackauckas
@ChrisRackauckas
Dec 18 2016 17:43
They must've gotten rid of that dispatch. Makes sense because it would probably only lead to errors.
Hmm.. it's still something that's pretty simple to implement. Can't think of GSoC projects for SymEngine.jl