Where communities thrive


  • Join over 1.5M+ people
  • Join over 100K+ communities
  • Free without limits
  • Create your own community
People
Repo info
Activity
  • Jun 24 23:54
    nkanazawa1989 starred symengine/symengine
  • Jun 24 21:13
    isuruf commented #1914
  • Jun 24 21:13
    isuruf commented #1914
  • Jun 24 21:06
    jmanthony3 commented #1914
  • Jun 24 21:02
    isuruf commented #1914
  • Jun 24 20:54
    jmanthony3 commented #1914
  • Jun 24 20:52
    isuruf review_requested #1913
  • Jun 24 20:52
    isuruf commented #1914
  • Jun 24 20:50
    jmanthony3 edited #1914
  • Jun 24 20:42
    jmanthony3 opened #1914
  • Jun 24 20:00
    isuruf reopened #1913
  • Jun 24 20:00
    isuruf synchronize #1913
  • Jun 24 19:32
    isuruf synchronize #1913
  • Jun 24 19:32
    isuruf closed #1913
  • Jun 24 15:36
    isuruf synchronize #1913
  • Jun 24 15:05
    isuruf synchronize #1913
  • Jun 24 14:54
    isuruf synchronize #1913
  • Jun 24 14:24
    isuruf opened #1913
  • Jun 24 11:07
    griels starred symengine/symengine
  • Jun 23 00:02
    MWalbecq opened #1912
Christopher Rackauckas
@ChrisRackauckas
okay
what about analytical ODEs?
Isuru Fernando
@isuruf
Do you have an example I can try?
that builds an analytical inverse "W"
Isuru Fernando
@isuruf
I mean an example you want simplified
Christopher Rackauckas
@ChrisRackauckas
in the field syminvW
let me see if there's an easy way to "save it"
Isuru Fernando
@isuruf
Print it. I can parse it.
Christopher Rackauckas
@ChrisRackauckas
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
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
CSE?
oh
Isuru Fernando
@isuruf
common subexpression elimination
Christopher Rackauckas
@ChrisRackauckas
that would be awesome
how would you do that here?
Isuru Fernando
@isuruf
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
symengine/symengine#1442
Christopher Rackauckas
@ChrisRackauckas
why does it need to be wrapped in C instead of using the C++ directly?
Isuru Fernando
@isuruf
When first writing the julia package, CXX.jl was really hard to install, so wrapped it in C and used ccall
Bjorn
@bjodah
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
@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
oh interesting.
Isuru Fernando
@isuruf
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
well because it's doing extra array creation and concatenation?
that's a memory wasteful way?
Isuru Fernando
@isuruf
there is?
Christopher Rackauckas
@ChrisRackauckas
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
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
looks like it's an array of tuples?
Isuru Fernando
@isuruf
gives a 2*2 array
Christopher Rackauckas
@ChrisRackauckas
oh, then it's just hard to read :P
I assumed you wrote
Isuru Fernando
@isuruf
it's like [a b; c d]
Christopher Rackauckas
@ChrisRackauckas
[(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
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
yup