I am in the process of switching all my code that I wrote in python to julia to decrease calculation time. For this I need to translate the following code from python over to Julia but since I am still new to Julia I am unsure what the current capabilities are before I reinvent the wheel.
The python code:
def calculate_2d_probability_over_volume(internal_matrix, px, py):
px = px[px >= 0]
py = py[py >= 0]
integrated_px = [np.trapz(px * row, px) for row in internal_matrix] # \int p_x * f(p_x, p_y) dp_x
return 2 / (2 * np.pi ** 2) * np.trapz(integrated_px[::-1], py) # \int integrated_px(p_y) dp_y
In the first two lines I make sure to only have positive values for px and py, they are both Numpy arrays with dimension N after the first two lines. The matrix is an NxN matrix, it correlates to the px and py values as follow:
let the value of the matrix be defined by a function f(px,py):
f(min(px),min(py)) and the
top right value of the matrix corresponds to f(max(px),max(py)).px, py are in ascending order, so min(px)==px[0]
if the matrix is represented by a two dimensional function f(px, py) then the integral to evaluate is given by:
b * ∫dpy∫dpx px * f(px,py)
with b being a factor that needs no further explanation.
The algorithm uses vectorization of numpy to evaluate it quickly. I now need to reimplement this algorithm within Julia's ecosystem.
Is there a function similar to what np.trapz() does already?
I would prefer if the function comes with an error estimation already too.
EDIT: fixed typo in row and column specification
This isn't a complete answer, but perhaps one which would serve as a starting point. Defining:
trapz(y) = @views sum((y[1:end-1].+y[2:end])/2)
trapz(y,x) = @views sum(((y[1:end-1].+y[2:end])/2).*(x[2:end].-x[1:end-1]))
Gives the same results as the NumPy version of this function:
julia> trapz([1.0,2.0,4.0])
4.5
julia> trapz([1.0, 2.0],[0.0,0.1])
0.15000000000000002