CODEWITHSUNDEEP
pythonnumpymultidimensional-arrayjuliaintegration
derdotte
derdotte

Reputation: 91

Write 2D integration over discrete points from python to julia

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):

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

Upvotes: 1

Views: 144

Answers (1)

Dan Getz
Dan Getz

Reputation: 18227

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

Upvotes: 1

Related Questions