Take the 2-minute tour ×
Computational Science Stack Exchange is a question and answer site for scientists using computers to solve scientific problems. It's 100% free, no registration required.

I'm looking into using artificial neural networks (ANN) to predict the reaction rates in my fluid instead of solving the full system of stiff ODEs. Some people from my lab have already done some work on that so I don't start from scratch but I am having problems with my applications. One of them I think relates to the quality of my dataset for training. We usually extract training data from CFD simulations that are either 1D/2D/3D. No matter what, we end up with a multidimensional array of data to feed to the neural network. To give you an idea of the size of the problem, I am looking into training 8 nets with 10 inputs and 1 output for each. I feel like a training set of about 100,000 points would be reasonable but the problem is that these 100,000 points need to cover a specific region of my multi-dimensional space. Simply extracting the whole flowfield from my simulations is not very good for 2 reasons:

  • For each snapshot, only a small portion of the points lie in the region where I need a high sampling to make sure my training is accurate
  • As I compile snapshots together, I end up with many near-duplicate points which (I believe) have a negative effect on my ANN training by a) biasing the training by putting more weight on these regions b) adding unnecessary points.

So I've been trying to filter the points I record before including them in my training set. As I see it, that involves checking whether a new point is within a certain n-dimensional radius of every point of my dataset. This brute force approach, which barring a few tricks scales like n^2, works so-so for extracting 10,000 points out of 100,000 (takes 30 min say) but breaks down as I increase the sizes and numbers of the snapshots... Clearly, there must be a more clever way of doing this, but I am not sure in which direction to start looking. I first tried with python and could move to FORTRAN to speed things up but I feel like I should look for a better strategy first. Is my only hope some kind of k-d tree? I have little to no experience with them and the problem that I see is that my tree will grow as I build my dataset and this can only increase the complexity. Would a python k-d tree library suit my need? Should I move to FORTRAN given the size of my problem? Any advice is appreciated, thank you !

share|improve this question
add comment

2 Answers

up vote 3 down vote accepted
+150

In molecular dynamics simulations we have the same problem: given a cutoff radius $r_c$, find all particles within at most $r_c$ of each other. The simplest $\mathcal O(n)$ approach is to divide the space into cells of edge length of at least $r_c$ and to compare every particle in every cell to all particles in the 26 neighboring cells (in three dimensions, for $d$ dimensions, this is $3^d-1$). By construction, particles in other cells will be farther away than $r_c$.

In your case where you're only testing a single position for a fixed set of points, divide the space into cells of edge length $\ge r_c$, sort the points into the cells (can be done in $\mathcal O(n)$), and then for each trial point, find the cell in which that point lies (can be done in $\mathcal O(1)$) and check the distance to the points in that cell and the surrounding cells (can be done in $\mathcal O(1)$, which is dependent on the point density and $r_c$, but not on the total number of points).

The whole procedure is described here. If you want more details, try googling for "cell linked list" or "cell list".

The k-d tree is also a good approach, but may be more difficult to implement on your own (there seems to be a Python implementation here, but it doesn't allow adding points). Don't worry about the complexity of the tree, however, since, for a more or less uniform point distribution, the search depth will behave as $\mathcal O(\log_2 n)$ for $n$ points. That is, doubling the amount of points will make the search effort grow by a constant factor. You also only have to construct the tree once and can update it quickly when adding new points.

share|improve this answer
    
Thanks I'm gonna look into this, although I'm afraid my high-dimensionality might be a hindrance. Considering d=10, and about 100 pts per "side" (and I'm not sure this is not too low of a number), that's a grid with 10^20 cells no? Even with 10 points per "side" that's still 10^10 cells and I have 59000 neighbors for each cell :) But I will look at your links. Thank you ! –  FrenchKheldar Feb 4 '12 at 7:39
    
Actually, the number of cells depends on your point density relative to the cutoff distance. You should never end up with more cells than points! The number of neighbors is an issue in high dimensions, but then I would move to k-d trees. –  Pedro Feb 4 '12 at 10:47
    
I meant 59000 neighboring cells sorry. I've found a nicely documented kdtree nearest-neighbour search algorithm cs.umd.edu/~mount/ANN I will give it a shot. –  FrenchKheldar Feb 4 '12 at 16:42
add comment

The problem you have stated is a classic problem in Computational Geometry: Range Queries. That is:

Input: a subset S of n-dimensional euclidean space, and a set of points in that space P.
Output: the subset of P that intersects S.

The book Algorithms in a Nutshell describes how to solve a similar problem on page 292. It describes an algorithm for a rectangular region S (not an n-dimensional sphere as in your case). You can obtain a solution which is $O(n^{1-1/d}+r)$, where n is the number of points, d is the dimension of the space, and r is the number of points reported by the query. If d is very large, the asymptotic performance is practically $O(n)$, thus is what's normally referred to as the "curse of dimensionality". That is, if the dimension is high, then performance suffers greatly! Still, $O(n^{1-\frac{1}{d}}+r)$ is still asymptotically faster than $O(n)$, just not as much as we wish it could be.

The algorithm involves a divide and conquer strategy (recursion) and a special data structure (kd-tree). Here is an outline of the algorithm:

Begin MainProgram

    results = a new SET   (a kd-tree)
    search(space,root,results)
    return results
endMainProgram

Subroutine search(space,node,results)
    if (space contains node.region) then
        add node.point to the results
        for each descendant of node
            add d.point to results
        return
    if (space contains node.point) then 
        add node.point to results
    if (space extends below node.coordinate) then
        search(space,node.below,results)
    if (space extends above node.coordinate) then
        search(space,node.above,results)
endsubroutine

-source: Algorithms in a Nutshell, pg 298

The key to the rectangular case is that we can easily define the data structure such that we can include an entire subset all at once... hence the KD-tree. The KD tree simply splits your n-dimensional space by cutting it successively by hyperplanes, much the same way a 3D space can be split by planes.

For your particular problem involving range queries in the radial direction (not rectangular boxes), you may be able to find a similar n-sphere based recursive space splitting... There is a similar data-structure called a vp-tree which is designed for implementation of space partitioning in hyper-spherical coordinates. You may want to see this publication for more details on the theory of vp-trees.

Of course, coding can be a hassle if your goal is really to use the algorithm as a tool to conduct science. In which case, I would suggest looking into libraries that already implement these data structures that you can use. A computational geometry library would be very useful in this situation. The CGAL library has subroutines for d-dimensional range searches which may be of interest to you. Here is another list of libraries with range query subroutines.

Alternatively, if you are ok with obtaining most of the points in the spherical range that you are searching (but not exactly all), you may want to consider using an approximation algorithm such the one in this paper.

share|improve this answer
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.