0

This is for a geophysical analysis program I am creating. I already have code to do all this, but I am looking for inspirations and ideas (good datastructures and algorithms).

What I want to model:

  • Velocity as a function of depth (z)
  • The model is built up from multiple layers (<10)
    • Every layer is accessible by an index going from 0 for the top most layer to n for the bottom most layer
  • Every layer has velocity as a linear function of depth (gradient a_k and axis intercept b_k of the kth layer)
  • Every layer has a top and bottom depth (z_k-1 and z_k)
  • The model is complete, there is no space between layers. The point directly between two layers belongs to the lower layer

Requirements:

  • Get velocity at an arbitrary depth within the model. This will be done on the order of 1k to 10k times, so it should be well optimized.
  • Access to the top and bottom depths, gradients and intercepts of a layer by the layer index

What I have so far:
I have working Python code where every layer is saved as a numpy array with the values of z_k (bottom depth), z_k-1 (top depth), a_k (velocity gradient) and b_k (axis intercept). To evaluate the model at a certain depth, I get the layer index (, use that to get the parameters of the layer and pass them to a function that evaluates the linear velocity gradient.

Olaf
  • 586
  • 5
  • 18
  • It seems like you just need a random access ability to the data, every operation you do is already `O(1)`. What else do you want to improve? – Yonlif Jun 27 '19 at 09:07
  • Evaluating the velocity model currently isn't `O(1)` since I need to find the layer index for which I use binary search on the z values to find the place where the depth would fit. I know the number of layers isn't large, but I think there has to be a better approach. – Olaf Jun 27 '19 at 09:13
  • Oh ok so I misunderstood your first requirement. What is the maximum depth? – Yonlif Jun 27 '19 at 09:17
  • There will be multiple different models, but the depth range I am working is <1 km to <10 km. – Olaf Jun 27 '19 at 09:19

1 Answers1

1

So you have piecewise linear dependence, where z-coordinates of pieces ends go irrregular, and want to get function value at given z.

Note that there is no sense to use binary search for 10 pieces (3-4 rounds of BS might be slower than 9 simple comparisons).

But what precision have your depth queries? Note that you can store a table both for 1-meter resolution and for 1 millimeter too - only 10^7 entries provide O(1) access to any precalculated velocity value

For limited number of pieces it is possible to make long formula (involving integer division) but results perhaps should be slower.

Example for arbitrary three-pieces polyline with border points 2 and 4.5:

f = f0 + 0.2*int(z/2.0)*(z-2.0) + 0.03*int(z/4.5)*(z-4.5)
MBo
  • 77,366
  • 5
  • 53
  • 86
  • I also have a version where I use comparisons, but I haven't benchmarked them yet. I can't control the precision for the queries, since one algorithm (ray tracing) solves an IVP for a system of ODEs and the solver (`scipy.integrate.solve_ivp`) uses a root-finding-algorithm around the layer boundaries where the step size varies. But I like the idea, keeping one table, rounding the depth to a given precision and using that to look up the value. Is there a way to string together the piecewise defined linear functions and evaluate without an if? – Olaf Jun 27 '19 at 09:40
  • Yes, it is possible to make single formula - I added info, but doubt in effectiveness – MBo Jun 27 '19 at 09:55