Add Gaussian Hypergeometric Function `Hyp2F1` #1288

ColtAllen · 2022-11-09T18:36:14Z

This PR is in reference to #1046 and adds Ops for scipy.special.hyp2f1, scipy.special.poch, and scipy.special.factorial. Currently all three are failing the test_grad tensor broadcasting tests in aesara/tests/tensor/test_math_scipy.py. Here are the results from pytest:

TestHyp2F1Broadcast & TestHyp2F1InplaceBroadcast

test_grad - ValueError: ('abs_err not finite', 'array([[nan, nan, nan],\n       [nan, nan, nan]])', '
Test Elemwise{hyp2f1,no_inplace}::normal: Error occurred while computing the gradient on the following inputs:
[array([[764.16214925, 550.49823533, 542.19109217],\n       [613.93532095, 341.6967123 , 284.38215306]]), 
array([[764.16214925, 550.49823533, 542.19109217],\n       [613.93532095, 341.6967123 , 284.38215306]]), 
array([[764.16214925, 550.49823533, 542.19109217],\n       [613.93532095, 341.6967123 , 284.38215306]]),
array([[0.76416215, 0.55049824, 0.54219109],\n       [0.61393532, 0.34169671, 0.28438215]])]')

TestPochBroadcast & TestPochInplaceBroadcast

test_grad - aesara.gradient.GradientError: GradientError: numeric gradient and analytic gradient exceed tolerance:
E                       At position 1 of argument 0 with shape (2, 1),
E                           val1 = -130.220047      ,  val2 = 172.029567
E                           abs. error = 302.249614,  abs. tolerance = 0.050000
E                           rel. error = 1.000000,  rel. tolerance = 0.000100
E               Exception args: Test Elemwise{poch,no_inplace}::normal: Error occurred while computing the gradient on the following inputs: [array([[0.87809993],
E                      [0.8534127 ]]), array([[2.58322029],
E                      [4.84803489]])]

TestFactorialBroadcast & TestFactorialInplaceBroadcast

test_grad - TypeError: ('verify_grad can work only with floating point inputs, but input 0 has dtype "int64".', 
'Test Elemwise{factorial,no_inplace}::int: 
Error occurred while computing the gradient on the following inputs: 
[array([[ 1,  6,  3],\n       [10,  4,  2]])]')

brandonwillard

Looks good. Seems like we might need to go over the test values being used and/or potentially some subtleties of those gradient implementations.

ColtAllen · 2022-11-10T19:04:17Z

Looks good. Seems like we might need to go over the test values being used and/or potentially some subtleties of those gradient implementations.

Thanks for the quick response. The docs for scipy.special.factorial suggest it only accepts integer arguments, but it works just fine with floats. Changing the testing input types for TestFactorialBroadcast is now throwing a tolerance error similar to that of TestPochBroadcast :

test_grad - aesara.gradient.GradientError: GradientError: numeric gradient and analytic gradient exceed tolerance:
E                       At position 1 of argument 0 with shape (2, 1),
E                           val1 = 14.836027      ,  val2 = 128.113961
E                           abs. error = 113.277934,  abs. tolerance = 0.050000
E                           rel. error = 0.792431,  rel. tolerance = 0.000100
E               Exception args: Test Elemwise{factorial_inplace,inplace}::normal: Error occurred while computing the gradient on the following inputs: [array([[2.50818975],
E                      [4.71063693]])]

What could be causing this GradientError? Perhaps makeBroadcastTester(eps=2e-10) should be adjusted? Does eps even need to be specified at all?

As for TestHyp2F1Broadcast, it seems some combinations of testing input values are shooting to infinity in the summation loop before that gradient expression converges. I'll check to see if the inputs are valid/reasonable and if another loop break criterion can be added.

ColtAllen · 2022-11-10T19:55:03Z

Turns out TestHyp2F1Broadcast also had unreasonable tests values. After making those changes, every Op now has a GradientError:

test_grad - aesara.gradient.GradientError: GradientError: numeric gradient and analytic gradient exceed tolerance:
E                       At position 0 of argument 0 with shape (2, 3),
E                           val1 = -48538346782196181189920139818200480218765603896390449205366505275561749456688053098819858451084244350695136091976910166574438712182302091308166577389179037024256.000000      ,  val2 = 3520008228282594454153617700017716374498340850794903970680924876145639825585099133611446954617413017396337770314707559185864868875633180509958375955496960000000.000000
E                           abs. error = 52058355010478775644073757518218196593263944747185353176047430151707389282273152232431305405701657368091473862291617725760303581057935271818124953344675997024256.000000,  abs. tolerance = 0.050000
E                           rel. error = 1.000000,  rel. tolerance = 0.000100
E               Exception args: Test Elemwise{hyp2f1,no_inplace}::normal: Error occurred while computing the gradient on the following inputs: [array([[764.16214925, 550.49823533, 542.19109217],
E                      [613.93532095, 341.6967123 , 284.38215306]]), array([[764.16214925, 550.49823533, 542.19109217],
E                      [613.93532095, 341.6967123 , 284.38215306]]), array([[764.16214925, 550.49823533, 542.19109217],
E                      [613.93532095, 341.6967123 , 284.38215306]]), array([[0.38208107, 0.27524912, 0.27109555],
E                      [0.30696766, 0.17084836, 0.14219108]])]

If I could get some more clarification on how to resolve these GradientErrors, I think we'll be set.

brandonwillard · 2022-11-10T20:26:09Z

Thanks for the quick response. The docs for scipy.special.factorial suggest it only accepts integer arguments, but it works just fine with floats.

Yeah, some version of the gamma function should be/is likely being used in that case.

brandonwillard · 2022-11-10T20:28:29Z

If I could get some more clarification on how to resolve these GradientErrors, I think we'll be set.

The magnitudes of those printed values might be indicating under/overflow of some sort. A numerical approximation is used for testing, and, if both values are so large, that could mean the test points are in an unstable/unsupported range of values for the given gradient representations (i.e. the exact functional forms of the gradients being used) and/or approximation routine.

brandonwillard · 2022-11-10T20:34:58Z

aesara/scalar/math.py

+class Factorial(UnaryScalarOp):
+    """
+    Factorial function of a scalar or array of numbers.
+
+    """
+
+    nfunc_spec = ("scipy.special.factorial", 1, 1)
+
+    @staticmethod
+    def st_impl(n):
+        return scipy.special.factorial(n)
+
+    def impl(self, n):
+        return Factorial.st_impl(n)
+
+    def grad(self, inputs, grads):
+        (n,) = inputs
+        (gz,) = grads
+        return [gz * gamma(n+1) * tri_gamma(n+1)]
+
+    def c_code(self, *args, **kwargs):
+        raise NotImplementedError()


We should see whether or not the existing Gamma works for this instead. For instance, we could change factorial below to a helper function that constructs Gamma graphs.

ColtAllen · 2022-11-11T13:43:41Z

I don't know if this is exactly kosher, but with makeBroadcastTester(grad=None), all tests are now passing. I've noticed not all tests in aesara/tests/tensor/test_math_scipy.py have this grad parameter set, so I'm unsure when/if it's required.

ColtAllen · 2022-11-11T13:51:35Z

We should see whether or not the existing Gamma works for this instead. For instance, we could change factorial below to a helper function that constructs Gamma graphs.

Like this?

@staticmethod
def st_impl(n):
    return gamma(n+1)

Because if the same is done for Poch, then Hyp2F1 could be defined entirely in terms of existing Ops if that infinite summation loop still works. If none of these Ops require SciPy implementations, would this even still be the appropriate module to add them? Would makeBroadcastTester(expected=scipy.special.hyp2f1) still be valid testing criterion?

brandonwillard · 2022-11-11T17:45:29Z

Like this?
@staticmethod
def st_impl(n):
    return gamma(n+1)
Because if the same is done for Poch, then Hyp2F1 could be defined entirely in terms of existing Ops if that infinite summation loop still works.

Yes, let's try that first. We may need to follow these changes with some rewrites that produce better gradients in certain cases, but that's good, because it better fits the design and intention of the library. Also, no worries if you're not familiar with the rewriting aspects of Aesara; we'll help with or follow-up on that.

If none of these Ops require SciPy implementations, would this even still be the appropriate module to add them? Would makeBroadcastTester(expected=scipy.special.hyp2f1) still be valid testing criterion?

Great question; the answer would be "yes"—to both questions perhaps—if those functions are part of the SciPy interface or have equivalents of some sort. If they aren't then aesara.tensor.math is fine.

brandonwillard · 2022-11-11T18:33:15Z

Sorry, I meant that we should replace Factorial with a helper function that creates gamma(n + 1) graphs using the existing Gamma Op. In other words, we wouldn't need a new Op for the factorial—or anything else we can derive using existing Ops.

ColtAllen · 2022-11-11T19:44:56Z

Sorry, I meant that we should replace Factorial with a helper function that creates gamma(n + 1) graphs using the existing Gamma Op. In other words, we wouldn't need a new Op for the factorial—or anything else we can derive using existing Ops.

So replace the current Factorial Op code with this, in effect?

factorial = aesara.function([n],gamma(n+1))

Doing this for Hyp2F1 is trickier because it requires a scan function for the infinite summation, but here's my first pass at it:

def converge(prev_result, result):
    return prev_result, aesara.scan.utils.until(prev_result==result)
            
results, updates = aesara.scan(
    fn = lambda n, a, b, c, z: if abs(z) >=1 raise NotImplementedError() else ((poch(a,n) * poch(b,n))/poch(c,n)) * ((z**n)/factorial(n)),
        outputs_info=at.constant(1.),
        non_sequences=[a,b,c,z],
        n_steps = 1024
        )
            
        final_result = results.sum()

hyp2f1 = aesara.function(inputs=[a,b,c,z], outputs=final_result , updates=updates)

I want to add that converge early termination function so that it doesn't loop 1,024 times for this 'infinite' summation, but I'm not quite sure how to include it.

brandonwillard · 2022-11-11T23:52:16Z

So replace the current Factorial Op code with this, in effect?
factorial = aesara.function([n],gamma(n+1))

Well, more like the following:

def factorial(x: TensorLike) -> TensorVariable:
    return gamma(x + 1)

A user could then construct a graph with at.factorial(...) and use aesara.function to compile it. The same goes for Poch.

Doing this for Hyp2F1 is trickier because it requires a scan function for the infinite summation, but here's my first pass at it:

We can continue pursuing a custom Op solution for hyp2f1 in this PR—at least until we find a non-series-based approach or improve the Scan while-loop situation more (e.g. it currently allocates large-ish arrays by default for all while loops, the Cython implementation is too inefficient for small numbers of iterations, etc.)

ColtAllen · 2022-11-12T00:29:02Z

Well, more like the following:

def factorial(x: TensorLike) -> TensorVariable:
    return gamma(x + 1)

The typing variables seem to be throwing circular import errors:

_________________________ ERROR collecting tests/tensor/test_math_scipy.py _________________________
tests/tensor/test_math_scipy.py:12: in <module>
    from aesara import function
aesara/__init__.py:120: in <module>
    from aesara import scalar, tensor
aesara/scalar/__init__.py:2: in <module>
    from .math import *
aesara/scalar/math.py:14: in <module>
    from aesara.tensor.var import TensorVariable
aesara/tensor/__init__.py:105: in <module>
    from aesara.tensor import (  # noqa
aesara/tensor/blas.py:165: in <module>
    from aesara.tensor.math import Dot, add, mul, neg, sub
aesara/tensor/math.py:1310: in <module>
    def erf(a):
aesara/tensor/elemwise.py:1813: in scalar_elemwise
    return construct(symbol[0])
aesara/tensor/elemwise.py:1795: in construct
    scalar_op = getattr(scalar, symbolname)
E   AttributeError: partially initialized module 'aesara.scalar' has no attribute 'erf' (most likely due to a circular import)

Dispensing with type hinting is raising a new error in aesara./tensor/elemwise.py:

________________________________ERROR collecting tests/tensor/test_math_scipy.py ________________________________
tests/tensor/test_math_scipy.py:16: in <module>
    from aesara.tensor import inplace
aesara/tensor/inplace.py:406: in <module>
    def factorial_inplace(n):
aesara/tensor/elemwise.py:1813: in scalar_elemwise
    return construct(symbol[0])
aesara/tensor/elemwise.py:1786: in construct
    inplace_scalar_op = scalar_op.__class__(transfer_type(0))
E   TypeError: function() missing required argument 'globals' (pos 2)

We can continue pursuing a custom Op solution for hyp2f1 in this PR—at least until we find a non-series-based approach or improve the Scan while-loop situation more (e.g. it currently allocates large-ish arrays by default for all while loops, the Cython implementation is too inefficient for small numbers of iterations, etc.)

If this is the case, Factorial and Poch can be set aside for now, since they're only included as building blocks for constructing Hyp2F1 natively. However, this is a significant dependency for one of the modules in my BTYD library, and it's in my best interest to optimize Hyp2F1. If you think Factorial and Poch should still be included in this PR, I'm happy to continue working out how to do so.

brandonwillard · 2022-11-12T00:39:11Z

The typing variables seem to be throwing circular import errors:

I was thinking of helper/graph-constructor functions like the following:

from typing import TYPE_CHECKING

import aesara
import aesara.tensor as at


if TYPE_CHECKING:
    from aesara.tensor import TensorLike, TensorVariable


def factorial(n: "TensorLike") -> "TensorVariable":
    return at.gamma(n + 1)


def poch(z: "TensorLike", m: "TensorLike") -> "TensorVariable":
    return at.gamma(z + m) / at.gamma(z)


z = at.scalar("z")
n = at.scalar("n")
res = (factorial(n), poch(z, n))

res_fn = aesara.function([z, n], res)

res_fn(2.0, 7)
# [array(5040.), array(40320.)]

ColtAllen · 2022-11-12T01:29:23Z

Importing aesara.tensor into aesara.scalar.math is still raising the same AttributeError as before. Perhaps these Ops should be defined in aesara.tensor.math instead? The only relation to scipy at this point is an expected value for testing.

brandonwillard · 2022-11-12T01:31:30Z

Importing aesara.tensor into aesara.scalar.math is still raising the same AttributeError as before. Perhaps these Ops should be defined in aesara.tensor.math instead? The only relation to scipy at this point is an expected value for testing.

Ah, yes, aesara.tensor.math would work, but aesara.tensor.special is probably the place those need to be.

ColtAllen · 2022-11-12T01:47:38Z

Ah, yes, aesara.tensor.math would work, but aesara.tensor.special is probably the place those need to be.

Defining factorial in aesara.tensor.special is raising the same TypeError: function() missing required argument 'globals' (pos 2) exception. Removing the inplace definition from aesara.tensor.inplace clears up that error, but raises a new one when running TestFactorialBroadcast::test_good:

AttributeError: ("'function' object has no attribute 'make_node'", 'Test Elemwise{factorial,no_inplace}::normal: Error occurred while making a node with inputs [array([[2.50818975],\n [4.71063693]])]')

All tests will pass if makeBroadcastTester(good = None) , but I'm not sure if this is the best idea.

brandonwillard · 2022-11-12T01:51:27Z

Ah, yes, aesara.tensor.math would work, but aesara.tensor.special is probably the place those need to be.

Defining factorial in aesara.tensor.special is raising the same TypeError: function() missing required argument 'globals' (pos 2) exception. Removing the inplace definition from aesara.tensor.inplace clears up that error, but raises a new one when running TestFactorialBroadcast::test_good:

AttributeError: ("'function' object has no attribute 'make_node'", 'Test Elemwise{factorial,no_inplace}::normal: Error occurred while making a node with inputs [array([[2.50818975],\n [4.71063693]])]')

All tests will pass if makeBroadcastTester(good = None) , but I'm not sure if this is the best idea.

We can't use those test constructors in this case, since they're probably only intended to be used with direct Op implementations. All we need are a few simple tests that make sure each helper function is at least defined correctly (i.e. it computes the quantities of interest). A pytest.mark.parametrized test that checks a few sample points against the corresponding scipy.special functions should be sufficient in each case.

ColtAllen · 2022-11-12T02:06:20Z

All we need are a few simple tests that make sure each helper function is at least defined correctly (i.e. it computes the quantities of interest). A pytest.mark.parametrized test that checks a few sample points against the corresponding scipy.special functions should be sufficient in each case.

Ok, that sounds straightforward. The current tests are also passing if the helper functions call the scalar Ops rather than the tensor equivalents, in which case I suppose they can continue to be defined in aesara.scalar.math.

What about the inplace definition in aesara.tensor.inplace? Tests won't even run unless it is removed.

ColtAllen · 2022-11-12T18:38:34Z

Helper functions and their respective tests have been written, but the Elemwise errors still persist:

AttributeError: ("'function' object has no attribute 'make_node'", 'Test Elemwise{factorial,no_inplace}::normal: Error occurred while making a node with inputs [array([[2.50818975],\n [4.71063693]])]')

These are such simple helper implementations for at.gamma that I've rewritten the derivatives for Hyp2F1 in terms of gamma without issue. If we dispense with the helper factorial and poch functions in this PR, I think it'll be ready to submit for review.

brandonwillard · 2022-11-12T21:52:14Z

aesara/scalar/math.py

+def poch(z: ScalarType, m: ScalarType) -> ScalarVariable:
+    """
+    Pochhammer symbol (rising factorial) function.
+
+    """
+    return gamma(z + m) / gamma(z)
+
+
+def factorial(n: ScalarType) -> ScalarVariable:
+    """
+    Factorial function of a scalar or array of numbers.
+
+    """
+    return gamma(n + 1)


These need to be in aesara.tensor.[special|math] and aesara.tensor.math.gamma needs to be used, because aesara.tensor.math.gamma is an Elemwise Op. That Elemwise version of gamma is created from the scalar aesara.scalar.math.Gamma Op via the scalar_elemwise decorator in aesara.tensor.math, which creates the Elemwise like gamma = Elemwise(Gamma, ...).

In general, Elemwise Ops are the most common way that ScalarOps are used in Aesara, since they add broadcasting to scalar operations.

brandonwillard · 2022-11-12T21:55:28Z

aesara/tensor/math.py

+@scalar_elemwise
+def poch(z, m):
+    """pochhammer symbol (rising factorial) function"""
+
+
+@scalar_elemwise
+def factorial(n):
+    """factorial function"""


The scalar_elemwise decorator will look for a ScalarOp instance with the name poch/factorial, and this is where the errors likely start. Instead, if we remove these and define poch and factorial using the Elemwise version of gamma (i.e. the one created by scalar_elemwise in aesara.tensor.math) then everything should work.

ColtAllen · 2022-11-13T02:13:40Z

Tests are passing! Still need to double-check everything, but this should be about ready for review.

ColtAllen · 2022-11-16T17:23:58Z

If there are no other comments, I'm going to merge the base branch changes and submit this for review.

brandonwillard · 2022-11-16T17:40:57Z

If there are no other comments, I'm going to merge the base branch changes and submit this for review.

Sounds good! N.B. You'll have to rebase onto upstream/main (i.e. no merge commits) for the branch to be mergeable.

brandonwillard

I just rebased, squashed, and removed the type hints. Looks like we need to clean up some other type-related things before we can add type hints to those functions.

codecov · 2022-11-23T03:33:08Z

Codecov Report

Merging #1288 (092fb69) into main (b64cb85) will increase coverage by 0.01%.
The diff coverage is 84.53%.

❗ Current head 092fb69 differs from pull request most recent head c692690. Consider uploading reports for the commit c692690 to get more accurate results

Additional details and impacted files

@@            Coverage Diff             @@
##             main    #1288      +/-   ##
==========================================
+ Coverage   74.69%   74.71%   +0.01%     
==========================================
  Files         194      194              
  Lines       49730    49826      +96     
  Branches    10527    10539      +12     
==========================================
+ Hits        37145    37226      +81     
- Misses      10262    10272      +10     
- Partials     2323     2328       +5

Impacted Files	Coverage Δ
aesara/scalar/math.py	`85.00% <82.55%> (-0.30%)`	⬇️
aesara/tensor/inplace.py	`100.00% <100.00%> (ø)`
aesara/tensor/math.py	`90.42% <100.00%> (+0.03%)`	⬆️
aesara/tensor/special.py	`90.90% <100.00%> (+0.26%)`	⬆️

ColtAllen · 2022-11-23T04:37:58Z

I just rebased, squashed, and removed the type hints. Looks like we need to clean up some other type-related things before we can add type hints to those functions.

Thanks! (Sorry for my delay in getting around to this.)

Looks like the testing coverage report is failing because the derivatives aren't covered for Hyp2F1. Could this be due to TestHyp2F1Broadcast = makeBroadcastTester(grad=None)? These derivatives are contained in the Hyp2F1Der Op, so perhaps a makeBroadcastTester should also be added for it as well?

Note that when/if scan is performant enough for inifinite summations, hyp2f1 can be moved to aesara.tensor.special and refactored in terms of poch and factorial, then Hyp2F1Der can be removed.

brandonwillard · 2022-11-24T18:16:48Z

Could this be due to TestHyp2F1Broadcast = makeBroadcastTester(grad=None)? These derivatives are contained in the Hyp2F1Der Op, so perhaps a makeBroadcastTester should also be added for it as well?

Yeah, that's probably it.

Note that when/if scan is performant enough for inifinite summations, hyp2f1 can be moved to aesara.tensor.special and refactored in terms of poch and factorial, then Hyp2F1Der can be removed.

Definitely. We have quite a few big improvements to Scan that we're waiting to push, so that time will come much sooner than later.

ColtAllen · 2022-11-30T19:53:23Z

Note that when/if scan is performant enough for inifinite summations, hyp2f1 can be moved to aesara.tensor.special and refactored in terms of poch and factorial, then Hyp2F1Der can be removed.

Definitely. We have quite a few big improvements to Scan that we're waiting to push, so that time will come much sooner than later.

Do you think these enhancements will be merged in the next 1-2 months? If so, I'd prefer to wait until that happens. Otherwise I'll be submitting an additional follow-up PR with the aforementioned revisions.

brandonwillard · 2022-11-30T20:03:44Z

Do you think these enhancements will be merged in the next 1-2 months? If so, I'd prefer to wait until that happens. Otherwise I'll be submitting an additional follow-up PR with the aforementioned revisions.

The question is really whether or not we'll have the necessary improvements available for all backends (i.e. C, JAX, and Numba) in that time, and I can say that we're primarily focusing on Numba and JAX at the moment, so 1-2 months might not hold for the C backend.

Regardless, multiple/staged PRs are perfectly fine—and often preferable—on our end, so no worries there.

brandonwillard

Following up on this: there are some Codecov annotations saying that we're missing lines in Hyp2F1 and Hyp2F1Der, so we need to either confirm that those annotations are erroneous and/or add tests that reach those lines.

After that, this is good to merge.

ColtAllen · 2022-12-01T01:24:20Z

Following up on this: there are some Codecov annotations saying that we're missing lines in Hyp2F1 and Hyp2F1Der, so we need to either confirm that those annotations are erroneous and/or add tests that reach those lines.

Tests need to be added for Hyp2F1Der.

The missing lines in Hyp2F1 are all in grad, and are just direct calls to Hyp2F1Der. Do they still require testing?

brandonwillard · 2022-12-01T19:22:59Z

The missing lines in Hyp2F1 are all in grad, and are just direct calls to Hyp2F1Der. Do they still require testing?

Yeah, I think that's how we can get the missing coverage for both classes.

ColtAllen · 2022-12-06T19:21:01Z

I've started writing these additional tests, and could use some help with the gradients and their respective tests.

Here's how the tests are set up for Hyp2F1:

_good_broadcast_quaternary_hyp2f1 = dict(
    normal=(
        random_ranged(0, 1000, (2, 3)),
        random_ranged(0, 1000, (2, 3)),
        random_ranged(0, 1000, (2, 3)),
        random_ranged(0, 0.5, (2, 3)),
    ),
)

TestHyp2F1Broadcast = makeBroadcastTester(
    op=at.hyp2f1,
    expected=expected_hyp2f1,
    good=_good_broadcast_quaternary_hyp2f1,
    grad=_good_broadcast_quaternary_hyp2f1,
    eps=2e-10,
)

TestHyp2F1InplaceBroadcast = makeBroadcastTester(
    op=inplace.hyp2f1_inplace,
    expected=expected_hyp2f1,
    good=_good_broadcast_quaternary_hyp2f1,
    grad=_good_broadcast_quaternary_hyp2f1,
    eps=2e-10,
    inplace=True,
)

They are returning this jumbled error for test_grad:

TypeError: ('float() argument must be a string or a number, not \'ScalarVariable\'\nApply node that caused the error: 
Elemwise{hyp2f1_der}(input 0, input 1, input 2, input 3, TensorConstant{(1, 1) of 0.0})\nToposort index: 0\nInputs types: 
[TensorType(float64, (2, 3)), TensorType(float64, (2, 3)), TensorType(float64, (2, 3)), TensorType(float64, (2, 3)), TensorType(float64, (1, 1))]\n
Inputs shapes: [(2, 3), (2, 3), (2, 3), (2, 3), (1, 1)]\nInputs strides: [(24, 8), (24, 8), (24, 8), (24, 8), (8, 8)]\n
Inputs values: [\'not shown\', \'not shown\', \'not shown\', \'not shown\', array([[0.]])]\nOutputs clients: [[Elemwise{Mul}[(0, 1)](random_projection, Elemwise{hyp2f1_der}.0)]]\n\n
Backtrace when the node is created (use Aesara flag traceback__limit=N to make it longer):\n 
File "/mnt/c/Users/colta/portfolio/aesara/tests/tensor/utils.py", line 585, in test_grad\n    
utt.verify_grad(\n  
File "/mnt/c/Users/colta/portfolio/aesara/tests/unittest_tools.py", line 70, in verify_grad\n   
 orig_verify_grad(op, pt, n_tests, rng, *args, **kwargs)\n 
File "/mnt/c/Users/colta/portfolio/aesara/aesara/gradient.py", line 1854, in verify_grad\n 
symbolic_grad = grad(cost, tensor_pt, disconnected_inputs="ignore")\n  File "/mnt/c/Users/colta/portfolio/aesara/aesara/gradient.py", line 623, in grad\n 
rval: Sequence[Variable] = _populate_grad_dict(\n  File "/mnt/c/Users/colta/portfolio/aesara/aesara/gradient.py", line 1434, in _populate_grad_dict\n 
rval = [access_grad_cache(elem) for elem in wrt]\n  File "/mnt/c/Users/colta/portfolio/aesara/aesara/gradient.py", line 1434, in <listcomp>\n    rval = [access_grad_cache(elem) for elem in wrt]\n
File "/mnt/c/Users/colta/portfolio/aesara/aesara/gradient.py", line 1387, in access_grad_cache\n 
term = access_term_cache(node)[idx]\n  File "/mnt/c/Users/colta/portfolio/aesara/aesara/gradient.py", line 1213, in access_term_cache\n 
input_grads = node.op.L_op(inputs, node.outputs, new_output_grads)\n\n
HINT: Use the Aesara flag `exception_verbosity=high` for a debug print-out and storage map footprint of this Apply node.', 
'Test Elemwise{hyp2f1,no_inplace}::normal: Error occurred while computing the gradient on the following inputs: 
[array([[764.16214925, 550.49823533, 542.19109217],\n       [613.93532095, 341.6967123 , 284.38215306]]), 
array([[764.16214925, 550.49823533, 542.19109217],\n       [613.93532095, 341.6967123 , 284.38215306]]),
array([[764.16214925, 550.49823533, 542.19109217],\n       [613.93532095, 341.6967123 , 284.38215306]]), 
array([[0.38208107, 0.27524912, 0.27109555],\n       [0.30696766, 0.17084836, 0.14219108]])]')

Here's the testing setup for the gradients Op Hyp2F1Der:

_good_broadcast_pentanary_hyp2f1_der = dict(
    normal=(
        random_ranged(0, 1000, (2, 3)),
        random_ranged(0, 1000, (2, 3)),
        random_ranged(0, 1000, (2, 3)),
        random_ranged(0, 0.5, (2, 3)),
        integers_ranged(-1, 3, (2, 3)),
    ),
)

TestHyp2F1DerBroadcast = makeBroadcastTester(
    op=at.hyp2f1_der,
    expected=expected_hyp2f1,
    good=_good_broadcast_pentanary_hyp2f1_der,
    eps=2e-10,
)

The last parameter for this Op is a flag indicating which variable the derivative is taken wrt, but the makeBroadcastTester utility doesn't seem to like it:
numpy.core._exceptions._UFuncOutputCastingError: Cannot cast ufunc 'hyp2f1' output from dtype('float64') to dtype('int64') with casting rule 'same_kind'

brandonwillard · 2023-01-01T23:11:14Z

@ColtAllen, right after your last update on this, it looks like you were asked by the PyMC [Labs] group to move this work to their fork of Aesara where it was later completed.

I don't know if you're aware, but they are downstreaming Aesara, so working on it here means that it will eventually show up there. Also, it should be possible to create a PR for both repos based on the same underlying branch, in which case there would be no delay or disparity.

I apologize if our review efforts weren't to your liking, especially since you seem to have felt the need to take the results elsewhere and leave things in this state. If you had external pressures and—for example—needed things to go faster, we would've preferred that you inform us so that we could have accommodated.

Regardless, we know that the situation created by PyMC [Labs] and their fork can be confusing and divisive, and—because of this—we're more than willing to help however we can.

ColtAllen · 2023-01-02T13:57:24Z

Hey @brandonwillard,

I apologize for my lack of communication in this matter. I spent most of December traveling and have been playing catch-up since I got back last night.

I created this PR because it's a key backend requirement for the btyd library I've been working on, for which I've also been the sole developer. The PyMC [Labs] team are working on a very similar library and reached out to me about merging efforts. I was thrilled to join them because communities are paramount to the success and survival of open-source projects (plus I prefer working with others over going solo).

It was determined in the downstream PR the gradients for the Hyp2F1 Op are best written in C or Stan for performance and numerical stability, and both languages fall outside my present skillset. I'm afraid there isn't much else I can contribute to this PR beyond copy/pasting what's been done downstream by others.

I attempted a PR beyond my abilities and ultimately abandoned one project in favor of another, so if anyone is at fault here, it's me. Your comments were very helpful while I was working on this - in fact, what I enjoyed most about this PR was checking my email and getting a notification you responded.

ColtAllen · 2023-01-16T02:43:26Z

Thought I'd wrap up this PR as a courtesy. Gradients and respective tests have been added.

brandonwillard

We can merge this, but we need to create a follow-up issue for a C/Cython/Numba implementation of Hyp2F1Der.impl.

ColtAllen mentioned this pull request Nov 9, 2022

Request: Add Gaussian Hypergeometric Function #1046

Closed

brandonwillard added the Op implementation Involves the implementation of an Op label Nov 9, 2022

brandonwillard reviewed Nov 9, 2022

View reviewed changes

brandonwillard reviewed Nov 10, 2022

View reviewed changes

brandonwillard linked an issue Nov 10, 2022 that may be closed by this pull request

Request: Add Gaussian Hypergeometric Function #1046

Closed

brandonwillard added the SciPy compatibility label Nov 11, 2022

brandonwillard reviewed Nov 12, 2022

View reviewed changes

brandonwillard changed the title ~~Request: Add Gaussian Hypergeometric Function~~ Add Gaussian Hypergeometric Function Hyp2F1 Nov 13, 2022

brandonwillard force-pushed the scipy_hyp2f1 branch from f66595f to ae44ce6 Compare November 23, 2022 02:08

brandonwillard marked this pull request as ready for review November 23, 2022 02:08

brandonwillard force-pushed the scipy_hyp2f1 branch from ae44ce6 to cd474a5 Compare November 23, 2022 02:17

brandonwillard reviewed Nov 23, 2022

View reviewed changes

brandonwillard force-pushed the scipy_hyp2f1 branch from cd474a5 to db46ad7 Compare November 23, 2022 02:27

brandonwillard reviewed Nov 30, 2022

View reviewed changes

brandonwillard force-pushed the scipy_hyp2f1 branch 3 times, most recently from d86e2e3 to 092fb69 Compare January 18, 2023 23:27

Add Hyp2F1 Op and poch, factorial functions

c692690

brandonwillard force-pushed the scipy_hyp2f1 branch from 092fb69 to c692690 Compare January 30, 2023 20:19

brandonwillard enabled auto-merge (rebase) January 30, 2023 20:19

brandonwillard approved these changes Jan 30, 2023

View reviewed changes

brandonwillard merged commit d789a5f into aesara-devs:main Jan 30, 2023

Add Gaussian Hypergeometric Function Hyp2F1 #1288

Add Gaussian Hypergeometric Function Hyp2F1 #1288

Conversation

ColtAllen commented Nov 9, 2022

TestHyp2F1Broadcast & TestHyp2F1InplaceBroadcast

TestPochBroadcast & TestPochInplaceBroadcast

TestFactorialBroadcast & TestFactorialInplaceBroadcast

brandonwillard left a comment • edited Loading

Choose a reason for hiding this comment

ColtAllen commented Nov 10, 2022

ColtAllen commented Nov 10, 2022

brandonwillard commented Nov 10, 2022

brandonwillard commented Nov 10, 2022 • edited Loading

brandonwillard Nov 10, 2022 • edited Loading

Choose a reason for hiding this comment

ColtAllen commented Nov 11, 2022

ColtAllen commented Nov 11, 2022

brandonwillard commented Nov 11, 2022

brandonwillard commented Nov 11, 2022 • edited Loading

ColtAllen commented Nov 11, 2022 • edited Loading

brandonwillard commented Nov 11, 2022 • edited Loading

ColtAllen commented Nov 12, 2022

brandonwillard commented Nov 12, 2022 • edited Loading

ColtAllen commented Nov 12, 2022

brandonwillard commented Nov 12, 2022

ColtAllen commented Nov 12, 2022

brandonwillard commented Nov 12, 2022 • edited Loading

ColtAllen commented Nov 12, 2022

ColtAllen commented Nov 12, 2022 • edited Loading

brandonwillard Nov 12, 2022 • edited Loading

Choose a reason for hiding this comment

brandonwillard Nov 12, 2022

Choose a reason for hiding this comment

ColtAllen commented Nov 13, 2022

ColtAllen commented Nov 16, 2022

brandonwillard commented Nov 16, 2022

brandonwillard left a comment

Choose a reason for hiding this comment

codecov bot commented Nov 23, 2022 • edited Loading

Codecov Report

ColtAllen commented Nov 23, 2022

brandonwillard commented Nov 24, 2022

ColtAllen commented Nov 30, 2022

brandonwillard commented Nov 30, 2022

brandonwillard left a comment

Choose a reason for hiding this comment

ColtAllen commented Dec 1, 2022

brandonwillard commented Dec 1, 2022

ColtAllen commented Dec 6, 2022

brandonwillard commented Jan 1, 2023 • edited Loading

ColtAllen commented Jan 2, 2023

ColtAllen commented Jan 16, 2023

brandonwillard left a comment

Choose a reason for hiding this comment

Add Gaussian Hypergeometric Function `Hyp2F1` #1288

Add Gaussian Hypergeometric Function `Hyp2F1` #1288

brandonwillard left a comment •

edited

Loading

brandonwillard commented Nov 10, 2022 •

edited

Loading

brandonwillard Nov 10, 2022 •

edited

Loading

brandonwillard commented Nov 11, 2022 •

edited

Loading

ColtAllen commented Nov 11, 2022 •

edited

Loading

brandonwillard commented Nov 11, 2022 •

edited

Loading

brandonwillard commented Nov 12, 2022 •

edited

Loading

brandonwillard commented Nov 12, 2022 •

edited

Loading

ColtAllen commented Nov 12, 2022 •

edited

Loading

brandonwillard Nov 12, 2022 •

edited

Loading

codecov bot commented Nov 23, 2022 •

edited

Loading

brandonwillard commented Jan 1, 2023 •

edited

Loading