优化在值列表中查找最接近的值

Question

我很欣赏下面有人问过类似的问题，但我看到的问题都没有超过数万行。即使代码很慢，这些虽然“大”相对微不足道。就我而言，我正在处理数百万到数亿行，因此即使是小的优化也可能非常有用。

到目前为止，我一直依赖的答案是 here 和 here。我所做的升级并不惊人......

对于一个广泛的概述，我有一个大约 2,000,000 行（格式为easting-northing-value）的参考文件，我从中寻找最接近ProcessedLineData 的值（我不能保证在我的列表中存在完全匹配，所以我需要为每个实例计算整个列表，以找到最接近的可行点）。

我的第一个选项是：

AsciiFile.Sample _closestSample = Modifications.Value_and_Ref_Calc.ReferenceFile.Dataset[0].Data
                            .OrderBy(t => Public.CalculateDistance(t.Easting, t.Northing, ProcessedLineData.Easting, ProcessedLineData.Northing))
                            .First();

在此之后，我认为第二个示例的 Aggregate 会更好，所以我选择了这个：

AsciiFile.Sample _closestSample = Modifications.Value_and_Ref_Calc.ReferenceFile.Dataset[0].Data
                            .Aggregate((x, y) => 
                              Public.CalculateDistance(x.Easting, x.Northing, ProcessedLineData.Easting, ProcessedLineData.Northing) 
                              < 
                              Public.CalculateDistance(y.Easting, y.Northing, ProcessedLineData.Easting, ProcessedLineData.Northing) 
                              ? x : y);

这仍然很慢（与上一个几乎相同）。然后我想到了使用ToLookup()，这也在其中一个示例中进行了描述，但我没有看到它会带来什么具体的好处。

如果相关，Public.CalculateDistance 是：

public static double CalculateDistance(double _x1, double _y1, double _x2, double _y2)
        {
            return Math.Sqrt(Math.Pow(_x1 - _x2, 2) + Math.Pow(_y1 - _y2, 2));
        }

在这些示例中，Modifications.Value_and_Ref_Calc.ReferenceFile.Dataset[0].Data 是我的 2,000,000 个参考点列表。

如果我没有找到最接近的值的数亿个ProcessedLineData 点（嗯，从具有数万个点的小文件到具有数十万个点的大文件），这将不那么重要百万分）...

最终目的：我正在寻找最接近的值，以便能够使用与特定 Modifications.Value_and_Ref_Calc.ReferenceFile.Dataset[0].Data 关联的海拔高度，并使用它来修改我的 ProcessedLineData 中的值。整个序列中最慢的部分肯定是搜索我最接近的值。

有什么明显的方法可以优化我遗漏的代码吗？

Answer 1

在不了解值的范围和精度的任何更多信息的情况下，或者对查找的分布与参考点列表更改进行过多假设的情况下，对 for 循环进行一些简单的优化可以使 100 次查找的速度提高约 30 倍更快OrderBy/First 代码：

将pld 用于您的ProcessedLineData，将data 用于Modifications.Value_and_Ref_Calc.ReferenceFile.Dataset[0].Data，您会得到：

var _closestSample = data[0];
var dist = (_closestSample.Easting - pld.Easting) * (_closestSample.Easting - pld.Easting) + (_closestSample.Northing - pld.Northing) * (_closestSample.Northing - pld.Northing);
for (int j2 = 1; j2 < data.Count; ++j2) {
    var y = data[j2];
    var ydist = (y.Easting - pld.Easting) * (y.Easting - pld.Easting) + (y.Northing - pld.Northing) * (y.Northing - pld.Northing);
    if (ydist < dist) {
        dist = ydist;
        _closestSample = y;
    }
}

我的时间超过 2,000,000 个条目 data 列表和 100 次查找，OrderBy/First 需要 2.22 秒，for 需要 0.06 秒，速度提高了 32 倍。

所以，我确信有比蛮力更好的方法，经过一番研究，我发现了莫顿密码和Hilbert Curves。一些工作使用希尔伯特曲线生成了一个SpatialIndex 类，使用莫顿指数生成了一个SpatialIndexMorton 类。我还将希尔伯特索引调整为仅索引 16-32 位，这提供了每秒的最佳查找。对于我的数据，莫顿曲线现在有点快。

使用相同的随机数据测试，我发现for 方法每秒可以进行 147 次查找，希尔伯特索引每秒可以进行 5634 次查找，而莫顿索引每秒可以进行 7370 次查找，超过 10,000 次查找和 2,000,000 个参考点。请注意，空间索引的设置时间约为 3 秒，因此对于很少的查找，使用for 进行暴力破解会更快 - 我在 468 次查找时获得了收支平衡时间。

为了使这个（有点）通用，我从地球坐标的 (C# 8.0) 接口开始，它提供了一些辅助方法：

public interface ICoordinate {
    double Longitude { get; set; }
    double Latitude { get; set; }

    public ulong MortonCode() {
        float f = (float)Latitude;
        uint ui;
        unsafe { // perform unsafe cast (preserving raw binary)
            float* fRef = &f;
            ui = *((uint*)fRef);
        }
        ulong ixl = ui;

        f = (float)Longitude;
        unsafe { // perform unsafe cast (preserving raw binary)
            float* fRef = &f;
            ui = *((uint*)fRef);
        }
        ulong iyl = ui;

        ixl = (ixl | (ixl << 16)) & 0x0000ffff0000ffffL;
        iyl = (iyl | (iyl << 16)) & 0x0000ffff0000ffffL;

        ixl = (ixl | (ixl << 8)) & 0x00ff00ff00ff00ffL;
        iyl = (iyl | (iyl << 8)) & 0x00ff00ff00ff00ffL;

        ixl = (ixl | (ixl << 4)) & 0x0f0f0f0f0f0f0f0fL;
        iyl = (iyl | (iyl << 4)) & 0x0f0f0f0f0f0f0f0fL;

        ixl = (ixl | (ixl << 2)) & 0x3333333333333333L;
        iyl = (iyl | (iyl << 2)) & 0x3333333333333333L;

        ixl = (ixl | (ixl << 1)) & 0x5555555555555555L;
        iyl = (iyl | (iyl << 1)) & 0x5555555555555555L;

        return ixl | (iyl << 1);
    }

    const int StartBitMinus1 = 31;
    const int EndBit = 16;

    //convert (x,y) to 31-bit Hilbert Index
    public ulong HilbertIndex() {
        float f = (float)Latitude;
        uint x;
        unsafe { // perform unsafe cast (preserving raw binary)
            float* fRef = &f;
            x = *((uint*)fRef);
        }

        f = (float)Longitude;
        uint y;
        unsafe { // perform unsafe cast (preserving raw binary)
            float* fRef = &f;
            y = *((uint*)fRef);
        }

        ulong hi = 0;
        for (int bitpos = StartBitMinus1; bitpos >= EndBit; --bitpos) {
            // extract s'th bit from x & y
            var rx = (x >> bitpos) & 1;
            var ry = (y >> bitpos) & 1;
            hi <<= 2;
            hi += (rx << 1) + (rx ^ ry);

            //rotate/flip a quadrant appropriately
            if (ry == 0) {
                if (rx == 1) {
                    x = ~x;
                    y = ~y;
                }

                //Swap x and y
                uint t = x;
                x = y;
                y = t;
            }
        }
        return hi;
    }
    [MethodImpl(MethodImplOptions.AggressiveInlining)]
    public double DistanceTo(ICoordinate b) =>
        Math.Sqrt((Longitude - b.Longitude) * (Longitude - b.Longitude) + (Latitude - b.Latitude) * (Latitude - b.Latitude));
    [MethodImpl(MethodImplOptions.AggressiveInlining)]
    public double Distance2To(ICoordinate b) => (Longitude - b.Longitude) * (Longitude - b.Longitude) + (Latitude - b.Latitude) * (Latitude - b.Latitude);

    public ICoordinate MakeNew(double plat, double plong);
}

public static class ICoordinateExt {
    public static ICoordinate Minus(this ICoordinate a, ICoordinate b) =>
        a.MakeNew(a.Latitude - b.Latitude, a.Longitude - b.Longitude);
    public static ICoordinate Plus(this ICoordinate a, ICoordinate b) =>
        a.MakeNew(a.Latitude + b.Latitude, a.Longitude + b.Longitude);
}

然后，实现接口的真实类（将替换为您的真实类）：

public class PointOfInterest : ICoordinate {
    public double Longitude { get; set; }
    public double Latitude { get; set; }

    public PointOfInterest(double plat, double plong) {
        Latitude = plat;
        Longitude = plong;
    }

    public ICoordinate MakeNew(double plat, double plong) => new PointOfInterest(plat, plong);
}

还有一个使用希尔伯特曲线将IEnumerable<ICoordinate> 转换为ICoordinate 的空间索引集合的类：

public class SpatialIndex {
    SortedList<ulong, List<ICoordinate>> orderedData;
    List<ulong> orderedIndexes;

    public SpatialIndex(IEnumerable<ICoordinate> data) {
        orderedData = data.GroupBy(d => d.HilbertIndex()).ToSortedList(g => g.Key, g => g.ToList());
        orderedIndexes = orderedData.Keys.ToList();
    }

    public ICoordinate FindNearest(ICoordinate aPoint) {
        var hi = aPoint.HilbertIndex();
        var nearestIndex = orderedIndexes.FindNearestIndex(hi);
        var nearestGuess = orderedData.Values[nearestIndex][0];
        var guessDist = (nearestGuess.Longitude - aPoint.Longitude) * (nearestGuess.Longitude - aPoint.Longitude) + (nearestGuess.Latitude - aPoint.Latitude) * (nearestGuess.Latitude - aPoint.Latitude);
        if (nearestIndex > 0) {
            var tryGuess = orderedData.Values[nearestIndex-1][0];
            var tryDist = (tryGuess.Longitude - aPoint.Longitude) * (tryGuess.Longitude - aPoint.Longitude) + (tryGuess.Latitude - aPoint.Latitude) * (tryGuess.Latitude - aPoint.Latitude);
            if (tryDist < guessDist) {
                nearestGuess = tryGuess;
                guessDist = tryDist;
            }
        }

        var offsetPOI = new PointOfInterest(guessDist, guessDist);
        var minhi = (aPoint.Minus(offsetPOI)).HilbertIndex();
        var minhii = orderedIndexes.FindNearestIndex(minhi);
        if (minhii > 0)
            --minhii;
        var maxhi = (aPoint.Plus(offsetPOI)).HilbertIndex();
        var maxhii = orderedIndexes.FindNearestIndex(maxhi);
        for (int j2 = minhii; j2 < maxhii; ++j2) {
            var tryList = orderedData.Values[j2];
            for (int j3 = 0; j3 < tryList.Count; ++j3) {
                var y = tryList[j3];
                var ydist = (y.Longitude - aPoint.Longitude) * (y.Longitude - aPoint.Longitude) + (y.Latitude - aPoint.Latitude) * (y.Latitude - aPoint.Latitude);
                if (ydist < guessDist) {
                    nearestGuess = y;
                    guessDist = ydist;
                }
            }
        }

        return nearestGuess;
    }
}

还有一个使用莫顿曲线的类似类：

public class SpatialIndexMorton {
    SortedList<ulong, List<ICoordinate>> orderedData;
    List<ulong> orderedIndexes;

    public SpatialIndexMorton(IEnumerable<ICoordinate> data) {
        orderedData = data.GroupBy(d => d.MortonCode()).ToSortedList(g => g.Key, g => g.ToList());
        orderedIndexes = orderedData.Keys.ToList();
    }

    public ICoordinate FindNearest(ICoordinate aPoint) {
        var mc = aPoint.MortonCode();
        var nearestIndex = orderedIndexes.FindNearestIndex(mc);
        var nearestGuess = orderedData.Values[nearestIndex][0];
        var guessDist = (nearestGuess.Longitude - aPoint.Longitude) * (nearestGuess.Longitude - aPoint.Longitude) + (nearestGuess.Latitude - aPoint.Latitude) * (nearestGuess.Latitude - aPoint.Latitude);
        if (nearestIndex > 0) {
            var tryGuess = orderedData.Values[nearestIndex-1][0];
            var tryDist = (tryGuess.Longitude - aPoint.Longitude) * (tryGuess.Longitude - aPoint.Longitude) + (tryGuess.Latitude - aPoint.Latitude) * (tryGuess.Latitude - aPoint.Latitude);
            if (tryDist < guessDist) {
                nearestGuess = tryGuess;
                guessDist = tryDist;
            }
        }

        var offsetPOI = new PointOfInterest(guessDist, guessDist);
        var minmc = (aPoint.Minus(offsetPOI)).MortonCode();
        var minmci = orderedIndexes.FindNearestIndex(minmc);
        if (minmci > 0)
            --minmci;
        var maxmc = (aPoint.Plus(offsetPOI)).MortonCode();
        var maxmci = orderedIndexes.FindNearestIndex(maxmc);
        for (int j2 = minmci; j2 < maxmci; ++j2) {
            var tryList = orderedData.Values[j2];
            for (int j3 = 0; j3 < tryList.Count; ++j3) {
                var y = tryList[j3];
                var ydist = (y.Longitude - aPoint.Longitude) * (y.Longitude - aPoint.Longitude) + (y.Latitude - aPoint.Latitude) * (y.Latitude - aPoint.Latitude);
                if (ydist < guessDist) {
                    nearestGuess = y;
                    guessDist = ydist;
                }
            }
        }

        return nearestGuess;
    }
}

还有一些辅助扩展方法：

public static class ListExt {
    public static int FindNearestIndex<T>(this List<T> l, T possibleKey) {
        var keyIndex = l.BinarySearch(possibleKey);
        if (keyIndex < 0) {
            keyIndex = ~keyIndex;
            if (keyIndex == l.Count)
                keyIndex = l.Count - 1;
        }
        return keyIndex;
    }
}

public static class IEnumerableExt {
    public static SortedList<TKey, TValue> ToSortedList<T, TKey, TValue>(this IEnumerable<T> src, Func<T, TKey> keySelector, Func<T, TValue> valueSelector) =>
        new SortedList<TKey, TValue>(src.ToDictionary(keySelector, valueSelector));    
}

最后，一些使用它的示例代码，您在data 中的引用，以及您在plds 中的查找值：

var hilbertIndex = new SpatialIndex(data);
var ans = new (ICoordinate, ICoordinate)[lookups];
for (int j1 = 0; j1 < lookups; ++j1) {
    ICoordinate pld = plds[j1];
    ans[j1] = (pld, hilbertIndex.FindNearest(pld));
}

更新：我修改了查找最近的算法以获取索引上目标点上方和下方的最接近点，而不是仅尝试上面的那个。这提供了另一个不错的加速。

Answer 2

您可以尝试通过使用网格将平面划分为相等的正方形来缩小每次搜索的范围。每个元素将被存储到适当正方形的桶中。然后你可以使用这个网格来执行搜索，从包含搜索点的正方形开始，向外螺旋，直到你在螺旋的周边找到一个或多个填充的正方形。这是一个幼稚的算法，但在特定条件下可以表现得非常好：

1. 平面中元素的分布必须是随机的（或接近随机的）。如果大多数元素都挤在少数人口稠密的地区，性能就会受到影响。

2. 搜索点将位于地图的一般人口区域。在很远很远的地方或在广阔的无人居住的地区搜索可能需要很长时间。

正确配置网格正方形的大小很重要。太大或太小的方块同样会损害性能。太小的方块也会增加内存消耗，因为每个填充的方块都需要自己的桶。在完美条件下，该算法的执行速度比简单循环快 1000 倍以上。

public class SpatialDictionary<T> : IEnumerable<T>
{
    private readonly Dictionary<(int, int), List<T>> _dictionary;
    private readonly double _squareSize;
    private readonly Func<T, (double, double)> _locationSelector;
    private int _count;

    public int Count => _count;

    public SpatialDictionary(
        double squareSize, Func<T, (double, double)> locationSelector)
    {
        if (squareSize <= 0)
            throw new ArgumentOutOfRangeException(nameof(squareSize));
        _squareSize = squareSize;
        _locationSelector = locationSelector
            ?? throw new ArgumentNullException(nameof(locationSelector));
        _dictionary = new Dictionary<(int, int), List<T>>();
    }

    public void Add(T item)
    {
        var (itemX, itemY) = _locationSelector(item);
        int keyX = checked((int)(itemX / _squareSize));
        int keyY = checked((int)(itemY / _squareSize));
        if (!_dictionary.TryGetValue((keyX, keyY), out var bucket))
        {
            bucket = new List<T>(1);
            _dictionary.Add((keyX, keyY), bucket);
        }
        bucket.Add(item);
        _count++;
    }

    public T FindClosest(double x, double y)
    {
        if (_count == 0) throw new InvalidOperationException();
        int keyX = checked((int)(x / _squareSize));
        int keyY = checked((int)(y / _squareSize));
        double minDistance = Double.PositiveInfinity;
        T minItem = default;
        int radius = 0;
        while (true)
        {
            checked { radius++; }
            foreach (var square in GetSquares(keyX, keyY, radius))
            {
                if (!_dictionary.TryGetValue(square, out var bucket)) continue;
                foreach (var item in bucket)
                {
                    var (itemX, itemY) = _locationSelector(item);
                    var distX = x - itemX;
                    var distY = y - itemY;
                    var distance = Math.Abs(distX * distX + distY * distY);
                    if (distance < minDistance)
                    {
                        minDistance = distance;
                        minItem = item;
                    }
                }
            }
            if (minDistance != Double.PositiveInfinity) return minItem;
        }
    }

    private IEnumerable<(int, int)> GetSquares(int x, int y, int radius)
    {
        if (radius == 1) yield return (x, y);
        for (int i = -radius; i < radius; i++)
        {
            yield return checked((x + i, y + radius));
            yield return checked((x - i, y - radius));
            yield return checked((x + radius, y - i));
            yield return checked((x - radius, y + i));
        }
    }

    public IEnumerator<T> GetEnumerator()
        => _dictionary.Values.SelectMany(b => b).GetEnumerator();

    IEnumerator IEnumerable.GetEnumerator() => GetEnumerator();
}

使用示例：

var spatialData = new SpatialDictionary<DataRow>(squareSize: 10.0,
    dr => (dr.Field<double>("Easting"), dr.Field<double>("Northing")));

foreach (DataRow dataRow in dataTable.Rows)
{
    spatialData.Add(dataRow);
}

DataRow result = spatialData.FindClosest(100.0, 100.0);