These are chat archives for symengine/symengine

30th
Apr 2018
Björn Dahlgren
@bjodah
Apr 30 2018 11:06
Thanks! Opened symengine/symengine#1441
Christopher Rackauckas
@ChrisRackauckas
Apr 30 2018 15:30
@isuruf SymEngine doesn't have a simplify routine, or things like analytical forms for ODEs yet?
or have those been added?
Isuru Fernando
@isuruf
Apr 30 2018 15:32
What type of simplification do you want? Got any examples
Christopher Rackauckas
@ChrisRackauckas
Apr 30 2018 15:32
Inverted Jacobian gives a crazy result, hit it with a simplify to make it a smaller equation
for code generation.
Isuru Fernando
@isuruf
Apr 30 2018 15:33
common subexpression elimination is there that can make it simpler
Christopher Rackauckas
@ChrisRackauckas
Apr 30 2018 15:33
any conversions like changing to simpler trig functions and stuff like that
?
Basically, how much of Mathematica FullSimplify is there?
Isuru Fernando
@isuruf
Apr 30 2018 15:34
Nothing like that
Christopher Rackauckas
@ChrisRackauckas
Apr 30 2018 15:34
okay
what about analytical ODEs?
Isuru Fernando
@isuruf
Apr 30 2018 15:34
Do you have an example I can try?
that builds an analytical inverse "W"
Isuru Fernando
@isuruf
Apr 30 2018 15:35
I mean an example you want simplified
Christopher Rackauckas
@ChrisRackauckas
Apr 30 2018 15:35
in the field syminvW
let me see if there's an easy way to "save it"
Isuru Fernando
@isuruf
Apr 30 2018 15:36
Print it. I can parse it.
Christopher Rackauckas
@ChrisRackauckas
Apr 30 2018 15:36
alright
SymEngine.Basic[(1 - b*d*x*y*internal_γ^2/((1- (-c + d*x)*internal_γ + b*d*x*y*internal_γ^2/(1 - (a
 - b*y)*internal_γ))*(1 - (a - b*y)*internal_γ)))/(1 - (a - b*y)*internal_γ) -b*x*internal_γ/((1 - (
-c + d*x)*internal_γ + b*d*x*y*internal_γ^2/(1 - (a - b*y)*internal_γ))*(1 - (a - b*y)*internal_γ));
 d*y*internal_γ/((1 - (-c + d*x)*internal_γ +b*d*x*y*internal_γ^2/(1 - (a - b*y)*internal_γ))*(1 -
(a - b*y)*internal_γ)) (1 - (-c + d*x)*internal_γ + b*d*x*y*internal_γ^2/(1 - (a - b*y)*internal_γ))
^(-1)]
Isuru Fernando
@isuruf
Apr 30 2018 15:46
cse can make this,
x0 = b*y 
x1 = (1 - (a - x0)*internal_γ)**(-1) 
x2 = d*x 
x3 = x0*x2*internal_γ**2 
x4 = (1 + x1*x3 - (-c + x2)*internal_γ)**(-1) 
x5 = x1*x4 
x6 = x5*internal_γ
x7 = x1*(1 - x3*x5)
x8 = -b*x*x6
x9 = d*y*x6
arr = [x7 x8; x9 x4]
Christopher Rackauckas
@ChrisRackauckas
Apr 30 2018 15:50
CSE?
oh
Isuru Fernando
@isuruf
Apr 30 2018 15:50
common subexpression elimination
Christopher Rackauckas
@ChrisRackauckas
Apr 30 2018 15:50
that would be awesome
how would you do that here?
Isuru Fernando
@isuruf
Apr 30 2018 15:51
We don't have a C wrapper for it yet. Just the C++ function
Btw, in SymPy, number of operations = 95, simplified = 77, cse = 26, simplify + cse = 26
so, just the cse is enough
Isuru Fernando
@isuruf
Apr 30 2018 16:04
symengine/symengine#1442
Christopher Rackauckas
@ChrisRackauckas
Apr 30 2018 16:08
why does it need to be wrapped in C instead of using the C++ directly?
Isuru Fernando
@isuruf
Apr 30 2018 16:09
When first writing the julia package, CXX.jl was really hard to install, so wrapped it in C and used ccall
Björn Dahlgren
@bjodah
Apr 30 2018 16:15
Another reason for using a C ABI: you can actually link against a library built by another compiler (whereas clang/libc++ cannot link with a library build using gcc/libstdc++)
Isuru Fernando
@isuruf
Apr 30 2018 19:47
@ChrisRackauckas, looks like julia's internal cse (LLVM's cse pass) does this by default.
function f(x, y, a, b, c, d, internal_γ)
    a1 = (1 - b*d*x*y*internal_γ^2/((1 - (-c + d*x)*internal_γ + b*d*x*y*internal_γ^2/(1 - (a - b*y)*internal_γ))*(1 - (a - b*y)*internal_γ)))/(1 - (a - b*y)*internal_γ)
    a2 = -b*x*internal_γ/((1 - (-c + d*x)*internal_γ + b*d*x*y*internal_γ^2/(1 - (a - b*y)*internal_γ))*(1 - (a - b*y)*internal_γ))
    a3 = d*y*internal_γ/((1 - (-c + d*x)*internal_γ + b*d*x*y*internal_γ^2/(1 - (a - b*y)*internal_γ))*(1 - (a - b*y)*internal_γ))
    a4 = (1 - (-c + d*x)*internal_γ + b*d*x*y*internal_γ^2/(1 - (a - b*y)*internal_γ))^(-1)
    return [a1 a2; a3 a4]
end

function g(x, y, a, b, c, d, internal_γ)
   x0 = b*y 
   x1 = (1 - (a - x0)*internal_γ)^(-1) 
   x2 = d*x 
   x3 = x0*x2*internal_γ^2 
   x4 = (1 + x1*x3 - (-c + x2)*internal_γ)^(-1) 
   x5 = x1*x4 
   x6 = x5*internal_γ
   x7 = x1*(1 - x3*x5)
   x8 = -b*x*x6
   x9 = d*y*x6
   return  [x7 x8; x9 x4]
end


@btime f(1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0)
@btime g(1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0)
Christopher Rackauckas
@ChrisRackauckas
Apr 30 2018 19:47
oh interesting.
Isuru Fernando
@isuruf
Apr 30 2018 19:48
There's a catch though.
function h(x, y, a, b, c, d, internal_γ)
           return [[(1 + -1 * (-1 * b * y + a) * internal_γ) ^ -1 * (1 + -1 * b * d * x * y * (1 + b * d * x * y * (1 + -1 * (-1 * b * y + a) * internal_γ) ^ -1 * internal_γ ^ 2 + -1 * (d * x + -1c) * internal_γ) ^ -1 * (1 + -1 * (-1 * b * y + a) * internal_γ) ^ -1 * internal_γ ^ 2) -1 * b * x * (1 + b * d * x * y * (1 + -1 * (-1 * b * y + a) * internal_γ) ^ -1 * internal_γ ^ 2 + -1 * (d * x + -1c) * internal_γ) ^ -1 * (1 + -1 * (-1 * b * y + a) * internal_γ) ^ -1 * internal_γ]; [d * y * (1 + b * d * x * y * (1 + -1 * (-1 * b * y + a) * internal_γ) ^ -1 * internal_γ ^ 2 + -1 * (d * x + -1c) * internal_γ) ^ -1 * (1 + -1 * (-1 * b * y + a) * internal_γ) ^ -1 * internal_γ (1 + b * d * x * y * (1 + -1 * (-1 * b * y + a) * internal_γ) ^ -1 * internal_γ ^ 2 + -1 * (d * x + -1c) * internal_γ) ^ -1]]
       end
is 6 times slower than f and g
Christopher Rackauckas
@ChrisRackauckas
Apr 30 2018 19:49
well because it's doing extra array creation and concatenation?
that's a memory wasteful way?
Isuru Fernando
@isuruf
Apr 30 2018 19:49
there is?
Christopher Rackauckas
@ChrisRackauckas
Apr 30 2018 19:50
yeah, ; is a vcat IIRC
oh wait nevermind
but there is extra array overhead
that's two arrays vcat
not a matrix
Isuru Fernando
@isuruf
Apr 30 2018 19:51
got it
function h2(x, y, a, b, c, d, internal_γ)
          return [(1 - b*d*x*y*internal_γ^2/((1 - (-c + d*x)*internal_γ + b*d*x*y*internal_γ^2/(1 - (a - b*y)*internal_γ))*(1 - (a - b*y)*internal_γ)))/(1 - (a - b*y)*internal_γ) -b*x*internal_γ/((1 - (-c + d*x)*internal_γ + b*d*x*y*internal_γ^2/(1 - (a - b*y)*internal_γ))*(1 - (a - b*y)*internal_γ)); d*y*internal_γ/((1 - (-c + d*x)*internal_γ + b*d*x*y*internal_γ^2/(1 - (a - b*y)*internal_γ))*(1 - (a - b*y)*internal_γ)) (1 - (-c + d*x)*internal_γ + b*d*x*y*internal_γ^2/(1 - (a - b*y)*internal_γ))^(-1)]
end
is the same as f and g
Christopher Rackauckas
@ChrisRackauckas
Apr 30 2018 19:52
looks like it's an array of tuples?
Isuru Fernando
@isuruf
Apr 30 2018 19:53
gives a 2*2 array
Christopher Rackauckas
@ChrisRackauckas
Apr 30 2018 19:53
oh, then it's just hard to read :P
I assumed you wrote
Isuru Fernando
@isuruf
Apr 30 2018 19:53
it's like [a b; c d]
Christopher Rackauckas
@ChrisRackauckas
Apr 30 2018 19:53
[(x7,x8);(x9,x4)]
so then yeah that's fine.
[a b; c d]
it does it by expressions
I wouldn't use the short form syntax in many of these cases though, but huge expressions like this are definitely the edge case :-/
Isuru Fernando
@isuruf
Apr 30 2018 19:56
So, the takeaway is that LLVM's cse does the same thing as symengine's cse which is equivalent to sympy's simplify+cse
Christopher Rackauckas
@ChrisRackauckas
Apr 30 2018 19:56
yup
that's good to know
does SymEngine have analytical ODE solving on the radar?
Isuru Fernando
@isuruf
Apr 30 2018 19:57
Nope
Btw, symengine's cse makes some assumptions that you can turn on in LLVM with @fastmath
Christopher Rackauckas
@ChrisRackauckas
Apr 30 2018 19:58
yeah
I assumed it would need @fastmath
Isuru Fernando
@isuruf
Apr 30 2018 19:59
In this case, there was no need, but might be needed in other cases.