CAS Best Practices

That’s a Python problem. Python has many problem of those lines, like the fact that NumPy even exists at all instead of just using the standard array type. The problem is that Python’s native types are pretty useless for most mathematical compute, so NumPy/TensorFlow/PyTorch/etc. all have their own scalar and array types, and that causes all sorts of weird incompatibilities. Adding to this issue is the fact that the interpreter doesn’t specialize on type information so the conversions can be costly.

Julia’s array and number types are actually usable because of the type-specialization on the parametric typing, so we exploit that and just use the wheel. Julia’s 1 is exact, and 2 is exact. 1/2 does a type promotion to a Float64, but the 1//2 version is the exact Rational{Int} version. That’s all type-generic via the type parameters, so big(1)//big(2) is the Rational{BigInt} version. Right now some of our stuff defaults to ::Real dispatches, but we have ways to allow @variables x::Rational{BigInt}. Rules are then specialized based on the number hierarchy. And then there’s a ton of packages which add alternative arithmetics which composes with, like @JeffreySarnoff’s ArbFloats.jl. By doing numbers like this, we then separate the implementation of numbers from the symbolic arithmetic and just rely on type inference.

The difficulty is just coaxing the dispatch mechanism to treat Sym{T} <: T so that tracing non-symbolic generic code gives behavior specific to the number type. This would be required so that @variables x::Complex then takes the general f(x) dispatch for complex numbers in Julia, which requires subtyping based on the type parameter. This is where overdubbing techniques come into play, i.e.:

and there’s a prototype:

That way, if there’s non-symbolic code which has different behavior for complex than non-complex, it’ll naturally trace through. Then, because the rules system already allows for ::T using dispatch for rule selection, all that’s left is to make type-specific rules (most are ::Real right now), and finally to do things like in solve_for(x::Term,Alg()) allow for specific dispatches where, for example, if all of the coefficients are Float64 just sending it to Flint (which is just the DiffEq dispatching architecture with type-based defaults handling).

6 Likes