Taichi: [Discussion] Taichi's standard math library

[archibate]

Concisely describe the proposed feature
Currently Taichi's math utility functions are very basic and thus not very rich. And the linear algebra part has some name nit problems, e.g. normalized() is a static method while norm() is not.
I would like to use glsl-alike functions to manuplate vectors so that users from the world of OpenGL can easily get used to our vector system and thus minimize learning cost.

Describe the solution you'd like (if any)
linear algebra:
ti.vec2(x, y) -> ti.Vector([x, y])
u.norm() -> ti.length(u)
u.dot(v) -> ti.dot(u, v)
ti.Vector.normalized(u) -> ti.normalize(u)
(u - v).norm() -> ti.distance(u, v)
ti.reflect(i, d), ti.refract(i, d, eta) (see examples/renderer_util.py)
...
some scalar functions:
step(), smoothstep(), clamp()...

Additional comments
GLSL built-in functions: https://blog.csdn.net/hgl868/article/details/7876257
GLSL-alike C++ header-only math library: https://github.com/g-truc/glm

[archibate]

I found ti.length already exist, could we rename it to ti.sizeof?

[archibate]

Having a math library with well functionality would definitely help people make use of Taichi better and reduce their mind workload.
We want to add functions, but not too much. For example, examples/renderer_utils.py contains some useful functions for renderers, but not so useful for other usages. We should avoid adding these domain specific functions since they are useless in other scenes. If user want them, they should implement by themselves, or import a external package like math.
What we really want is something like clamp, step, distance, which can also be helpful in other usages, not only for renderers.
Once the standard is settled, it's hard to go back. We need discussions to decide what functions should be put in our built-in math library, and what should not.
@yuanming-hu @k-ye @xumingkuan @KLozes @masahi @Psycho7 What do you think, guys? Inputs are welcome!

[yuanming-hu]

Thanks for raising this.

Here are my personal opinions:

• Although glsl's math functions are battle-tested, I don't think many new taichi users are from the world of glsl. So I'm not sure if adding glsl-like functions will lower the learning barrier.
• If the goal is to bring convenience while not introducing too many things to learn, then we should only provide the most commonly used math functions. For example, adding distance and step may be convenient for some people, but I guess for most people they will be confused and have to look up the doc to learn what these functions mean.

Some other things to consider:

• Should scalar functions such as abs take matrices (element-wise operation) as well?
• Should we use ti.dot(a, b) or a.dot(b) ? Same thing to think about for cross

[archibate]

Should scalar functions such as abs take matrices (element-wise operation) as well?

Whether we support or not depends on:

1. whether it's often used.
2. whether it's irreplaceable.

I think element-wise-function-application aren't quite often-used, but irreplaceable because:

1. Write a for-loop to do-element-wise-operations by hand is not easy.
2. Implement FPL-ish m.map(lambda x: abs(x)), which is an overkill.

Should we use ti.dot(a, b) or a.dot(b) ? Same thing to think about for cross

I prefer ti.dot(a, b) than a.dot(b) since sometimes we may have (a - b).dot(c), and ti.dot shows the symmtric between a and b. Also apply to a.norm().

[yuanming-hu]

I think element-wise-function-application aren't quite often-used, but irreplaceable because:

1. Write a for-loop to do-element-wise-operations by hand is not easy.
2. Implement FPL-ish m.map(lambda x: abs(x)), which is an overkill.

Right, I think we should probably implement these scalar functions on matrices/vectors. I think the semantic meaning of these operations (e.g. ti.sin(matrix)) is easy to understand, and sometimes we do need these. I believe @znah also supports this.

[archibate]

Right, I think we should probably implement these scalar functions on matrices/vectors.

To do this, we may modify unary:
https://github.com/taichi-dev/taichi/blob/a03fa13139c8f48eead8af8d41e1518a659be30c/python/taichi/lang/ops.py#L19-L24
to something like:

def unary(foo): 
    def func(expr):
        if ti.is_taichi_class(expr):  # ti.Matrix, or ti.Array in #777
            return expr.element_wise(foo)
        else:
            return foo(Expr(expr)) 

    unary_ops.append(func) 
    return func 

And the same thing could be applied to binary.

[yuanming-hu]

Right, implementing element-wise operators for unary and binary should be straightforward. If anyone is interested please go ahead :-)

[archibate]

As you mentioned in https://github.com/taichi-dev/taichi/pull/937#issuecomment-626209818, ti.inv is for matrix inverse? ti.inverse or ti.Matrix.inverse may be a better name.
Also about rcp and rsqrt, do they really have performance difference with 1/x and 1/sqrt(x)? Or we don't expose rcp/rsqrt, just let the optimizer to create them?

[yuanming-hu]

Sorry about my delayed reply.

As you mentioned in #937 (comment), ti.inv is for matrix inverse? ti.inverse or ti.Matrix.inverse may be a better name.

I guess the debate should be: Do we exposed all matrix operations in the ti scope or in the ti.Matrix scope? We need a standardization here.

I personally vote for the ti scope (i.e. ti.inverse(mat) instead of ti.Matrix.inverse(mat), for convenience.

Also about rcp and rsqrt, do they really have performance difference with 1/x and 1/sqrt(x)?

Yes, the performance difference is huge. For example, on Intel CPUs

• rcp: https://software.intel.com/sites/landingpage/IntrinsicsGuide/#text=_mm256_rcp_ps&expand=4483
• div: https://software.intel.com/sites/landingpage/IntrinsicsGuide/#text=_mm256_div_ps&expand=4483,4804,2159
• rsqrt: https://software.intel.com/sites/landingpage/IntrinsicsGuide/#text=_mm256_rsqrt_ps&expand=4483,4804
• sqrt: https://software.intel.com/sites/landingpage/IntrinsicsGuide/#text=_mm256_sqrt_ps&expand=4483,4804,2159,5367

We probably need to expose these for performance-aware people.

Or we don't expose rcp/rsqrt, just let the optimizer to create them?

We can also add a compiler flag for these: See https://docs.nvidia.com/cuda/cuda-c-programming-guide/index.html

In general, code compiled with -ftz=true (denormalized numbers are flushed to zero) tends to have higher performance than code compiled with -ftz=false. Similarly, code compiled with -prec div=false (less precise division) tends to have higher performance code than code compiled with -prec div=true, and code compiled with -prec-sqrt=false (less precise square root) tends to have higher performance than code compiled with -prec-sqrt=true. The nvcc user manual describes these compilation flags in more details.

[yuanming-hu]

Thanks for all the valuable discussions. I think it's time to make a decision here.

High-level guidelines:

• Be Pythonic: "There should be one --- and preferably only one --- obvious way to do it."
• The API should be clear at first glance
  
  
  
  • Avoid unclear abbreviations
  
  
  • No in-place matrix operations should be provided (to avoid confusion)
  
  
• Try not to put too many things in the ti namespace
• No operations should be done via ti.Matrix.method(mat)
• Convenience is important

Based on the inputs above, I suggest the following set of Matrix API:

• [x] a.cross(b), a.dot(b) instead of ti.cross(a, b), ti.dot(a, b)
• [ ] U, sig, V = A.svd(dt) instead of U, sig, V = ti.svd(A, dt)
• [ ] R, S = A.polar_decompose(dt) instead of R, S = ti.polar_decompse(A, dt)
• [x] A.transpose() instead of A.T()/A.transposed()/ti.transpose(A)/ti.Matrix.transpose(A)
• [x] A.inverse() instead of A.inv()/A.inversed()/ti.inversed(A)/ti.Matrix.inversed(A)
  
  
  
  • Note that A**-1 computes the element-wise inverse of A instead of the inverse matrix of A. Maybe we should forbid this usage to avoid confusion. For example, ti.Matrix([[3, 4]]) ** -1 gives you [[1/3, 1/4]]
  
  
• [x] A.normalized() instead of A.normalize()/ti.normalize(A)/ti.Matrix.normalize(A)
• [x] A.trace() instead of A.tr()/ti.tr(A)/ti.trace(A)
• [x] A.determinant() instead of A.det()/ti.det(A)/ti.determinant(A)
• [x] A.norm(...) instead of ti.Matrix.norm(A)
• Use existing unary/binary operations for element-wise matrix computations. For example,
  
  
  
  • [x] ti.sin(A) instead of A.sin()
  
  
  • [x] abs(A) instead of ti.abs(A)/A.abs()
  
  
  • [x] max(A, B) instead of ti.max(A, B)/A.max(B) (element-wise max)
  
  
  • [ ] Exception: A.cast(data_type) instead of ti.cast(A, data_type)
  
  

[archibate]

A.norm()

I prefer A.length(), not sure if norm is understood to all.

Exception: A.cast(data_type) instead of ti.cast(A, data_type)

So do we have Expr.cast()? Also Expr.bit_cast()?

[yuanming-hu]

A.norm()

I prefer A.length(), not sure if norm is understood to all.

Eigen also uses norm: https://eigen.tuxfamily.org/dox/group__QuickRefPage.html

Exception: A.cast(data_type) instead of ti.cast(A, data_type)

So do we have Expr.cast()? Also Expr.bit_cast()?

Good point. I think cast needs more discussions. It's not a simple unary operator. Maybe you are right, we should use Expr/Matrix.cast() instead of ti.cast().

[archibate]

Btw, I think we should also consider:
len(a) instead of ti.length(a), which may be confused with vector length / norm.
a.append(b) instead of ti.append(a, b).
For dynamic SNodes, wdyt?

[yuanming-hu]

Btw, I think we should also consider:
len(a) instead of ti.length(a), which may be confused with vector length / norm.
a.append(b) instead of ti.append(a, b).
For dynamic SNodes, wdyt?

Great idea! These belong to sparse computation API (instead of math API) though.

[Rullec]

Hi all,
Sorry for my disturbuting!
IMAO my previous PR maybe is nearly finished, so I'm hunting for my next target to contribute.
Do you think if I can do something here?I'm very glad to help in this issue.:)

Bests,
Xudong

[TH3CHARLie]

so I'm hunting for my next target to contribute.

I am actually looking for places to work with as well since my progress on obsoleting ti.cfg is a little bit hung(but I am confident finishing it on time). I think it would be great to record https://github.com/taichi-dev/taichi/issues/833#issuecomment-636476353 as TODO lists so we can pick them to implement and track progress. @yuanming-hu

[yuanming-hu]

Thank you both for offering your help!

@Rullec For this issue, I think we need a bit more discussion before moving forward. Meanwhile, could you surf the forum and QQ groups, and collect common installation issues and solutions (such as https://forum.taichi.graphics/t/import-taichi/365/2) and add them to https://github.com/taichi-dev/taichi/blob/master/docs/install.rst#troubleshooting? That will be very helpful to the students. #1116 is another TODO-item which will be very helpful.

@TH3CHARLie could you help with #938 (Add bound checks for tensors with coordinate offsets)? You'll need to modify https://github.com/taichi-dev/taichi/blob/master/taichi/transforms/check_out_of_bound.cpp This needs a bit of understanding of Taichi's IR.