ENH: Function Operations

[giovaniceotto]

Pull request type

Please check the type of change your PR introduces:

• Code base additions (bugfix, features)
• Code maintenance (refactoring, formatting, renaming, tests)
• ReadMe, Docs and GitHub maintenance
• Other (please describe):

Pull request checklist

Please check if your PR fulfills the following requirements, depending on the type of PR:

• Code base additions (for bug fixes / features):

  • Tests for the changes have been added
  • Docs have been reviewed and added / updated if needed
  • Lint (black rocketpy) has passed locally and any fixes were made
  • All tests (pytest --runslow) have passed locally

What is the current behavior?

Enter text here...

What is the new behavior?

Enter text here...

Does this introduce a breaking change?

• Yes
• No

Other information

Enter text here...

[MateusStano]

Okay so I changed a few things on this one:

c3c4617 -> created the Function.xArray and Function.yArray which substitutes all the calls of Function.source[:,0] and Function.source[:,1] (there were a lot of them). I think the real right way to do this would be to use @cached_property in a lot of attributes but I think a major restructure in the Function class would be needed to do this.

3550b0f -> Improved and fixed some methods. Mostly change some numpy functions used to others that are faster.

---

Here are some problems:

• The integral method is currently really weird:

RocketPy/rocketpy/Function.py

Lines 2193 to 2283 in 1126751

	def integral(self, a, b, numerical=False):
	"""Evaluate a definite integral of a 1-D Function in the interval
	from a to b.
	Parameters
	----------
	a : float
	Lower limit of integration.
	b : float
	Upper limit of integration.
	numerical : bool
	If True, forces the definite integral to be evaluated numerically.
	The current numerical method used is scipy.integrate.quad.
	If False, try to calculate using interpolation information.
	Currently, only available for spline and linear interpolation. If
	unavailable, calculate numerically anyways.
	Returns
	-------
	ans : float
	Evaluated integral.
	"""
	if self.__interpolation__ == "linear" and not numerical:
	Xs = np.linspace(a, b, int((b - a) * 5))
	Ys = self.getValue(Xs)
	ans = np.trapz(Ys, x=Xs)
	elif self.__interpolation__ == "spline" and not numerical:
	# Integrate using spline coefficients
	xData = self.xArray
	yData = self.yArray
	coeffs = self.__splineCoefficients__
	ans = 0
	# Check to see if interval starts before point data
	if a < xData[0]:
	if self.__extrapolation__ == "constant":
	ans += yData[0] * (xData[0] - a)
	elif self.__extrapolation__ == "natural":
	c = coeffs[:, 0]
	subB = a - xData[0] # subA = 0
	ans -= (
	(c[3] * subB**4) / 4
	+ (c[2] * subB**3 / 3)
	+ (c[1] * subB**2 / 2)
	+ c[0] * subB
	)
	else:
	# self.__extrapolation__ = 'zero'
	pass
	# Integrate in subintervals between Xs of given data up to b
	i = 0
	while i < len(xData) - 1 and xData[i] < b:
	if b < xData[i + 1]:
	subB = b - xData[i] # subA = 0
	else:
	subB = xData[i + 1] - xData[i] # subA = 0
	c = coeffs[:, i]
	subB = xData[i + 1] - xData[i] # subA = 0
	ans += (
	(c[3] * subB**4) / 4
	+ (c[2] * subB**3 / 3)
	+ (c[1] * subB**2 / 2)
	+ c[0] * subB
	)
	i += 1
	# Check to see if interval ends after point data
	if b > xData[-1]:
	if self.__extrapolation__ == "constant":
	ans += yData[-1] * (b - xData[-1])
	elif self.__extrapolation__ == "natural":
	c = coeffs[:, -1]
	subA = xData[-1] - xData[-2]
	subB = b - xData[-2]
	ans -= (
	(c[3] * subA**4) / 4
	+ (c[2] * subA**3 / 3)
	+ (c[1] * subA**2 / 2)
	+ c[0] * subA
	)
	ans += (
	(c[3] * subB**4) / 4
	+ (c[2] * subB**3 / 3)
	+ (c[1] * subB**2 / 2)
	+ c[0] * subB
	)
	else:
	# self.__extrapolation__ = 'zero'
	pass
	else:
	# Integrate numerically
	ans, _ = integrate.quad(self, a, b, epsabs=0.01, limit=10000)
	return ans

First I had to add Xs and Ys to make linearly interpolated functions actually consider the a and b arguments that represent the interval of integration. Now linear interpolation integration is working fine. However, spline interpolated Functions still returns wrong values.

Why not always use scipy. integrate.quad? From what I gathered, this is slower than numpy.trapz but more accurate. However, since we need to chop the array into the intended interval to use numpy.trapz I can't be certain that it is actually faster then scipy.integrate.quad. Also the integration error with numpy.trapz can end up being pretty significant when doing operations on top of itn (like in integralFunction or averageFunction).

• Interpolation is always set to linear when using derivativeFunction is this intended?

• findInput needs to have a good starting guess to work well. Which is as expected, however when we use inverseFunction extra care is needed. Just important to point out

• The averageFunction is a bit problematic because the lower argument can never be equal to x, or in the average calculation there will be a division by 0. Now, I dont think it is worth doing something to deal with exceptions for now but this means that you cant plot these functions without specifing the plot limits. It is important to remember this when using the method

[giovaniceotto]

• Issues with the integral method, which were focused on exact integration of spline interpolated Functions, are now solved.
  • e5eefae
  • 217741e
  • 4f11816
• Thorough testing has been added for the Function.integral method.
  • e7c2d98
  • 2024ab5

[MateusStano]

I think all is good now!

Just one question:

RocketPy/rocketpy/Function.py

Lines 2766 to 2769 in e5eefae

	if __name__ == "__main__":
	import doctest
	doctest.testmod()

What is going on here?

[MateusStano]

Also, I found this in the enh/liquid-motors branch:

RocketPy/rocketpy/Function.py

Lines 2282 to 2304 in 2381a6d

	def differentiate(self, x, dx=1e-6, order=1):
	"""Differentiate a Function object at a given point.
	Parameters
	----------
	x : float
	Point at which to differentiate.
	dx : float
	Step size to use for numerical differentiation.
	order : int
	Order of differentiation.
	Returns
	-------
	ans : float
	Evaluated derivative.
	"""
	if order == 1:
	return (self.getValue(x + dx) - self.getValue(x - dx)) / (2 * dx)
	elif order == 2:
	return (
	self.getValue(x + dx) - 2 * self.getValue(x) + self.getValue(x - dx)
	) / dx**2

Should this be added here too?

[giovaniceotto]

I think all is good now!

Just one question:

RocketPy/rocketpy/Function.py

Lines 2766 to 2769 in e5eefae

	if __name__ == "__main__":
	import doctest
	doctest.testmod()

What is going on here?

This is doctest! It's there so that we can verify if all examples in the docstrings are working as expected.

More info here: https://docs.python.org/3/library/doctest.html