algorithm हिल्बर्ट वक्र पर एक बिंदु पर मैपिंग एन-आयामी मान




math hilbert-curve (5)

मेरे पास एन-आयामी बिंदुओं का एक विशाल सेट है (लाखों के दस; एन 100 के करीब है)।

स्थानिक इलाके को संरक्षित करते समय मुझे इन बिंदुओं को एक आयाम में मैप करने की आवश्यकता है। मैं इसे करने के लिए हिल्बर्ट स्पेस-भरने वक्र का उपयोग करना चाहता हूं।

प्रत्येक बिंदु के लिए मैं वक्र पर सबसे नज़दीकी बिंदु चुनना चाहता हूं। बिंदु के हिल्बर्ट मूल्य (वक्र की शुरुआत से वक्र की शुरुआत से वक्र लंबाई) एक एकल आयाम मान है जिसे मैं चाहता हूं।

गणना को तत्काल नहीं होना चाहिए, लेकिन मुझे आशा है कि यह सभ्य आधुनिक घरेलू पीसी हार्डवेयर पर कई घंटों से अधिक न हो।

कार्यान्वयन पर कोई सुझाव? क्या कोई पुस्तकालय है जो मेरी मदद करेगा? (भाषा ज्यादा मायने रखती नहीं है।)


N-> 1 और 1-> n से मैपिंग के लिए एल्गोरिदम "हिल्बर्ट स्पेस-भरने वक्र का उपयोग करके एक और एन-आयामी मानों के बीच मैपिंग की गणना" जेके लॉडर

यदि आप "एसएफसी मॉड्यूल और केडलिया ओवरले" के लिए Google को एक समूह पाते हैं जो इसे अपने सिस्टम के हिस्से के रूप में उपयोग करने का दावा करता है। यदि आप स्रोत देखते हैं तो आप शायद प्रासंगिक फ़ंक्शन निकाल सकते हैं।


यह मुझे स्पष्ट नहीं है कि यह वही करेगा जो आप चाहते हैं। इस त्रिभुज 3 डी मामले पर विचार करें:

001 ------ 101
 |\         |\
 | \        | \
 |  011 ------ 111
 |   |      |   |
 |   |      |   |
000 -|---- 100  |
  \  |       \  |
   \ |        \ |
    010 ------ 110

जो निम्नलिखित पथ से "हिल्बर्टाइज्ड" हो सकता है:

001 -----> 101
  \          \
   \          \
    011        111
     ^          |
     |          |
000  |     100  |
  \  |       \  |
   \ |        \ V
    010        110

1 डी आदेश में:

000 -> 010 -> 011 -> 001 -> 101 -> 111 -> 110 -> 100

यहाँ बुरा बिट है। जोड़े और 1 डी दूरी की सूची पर विचार करें:

000 : 100 -> 7
010 : 110 -> 5
011 : 111 -> 3
001 : 101 -> 1

सभी मामलों में, बाएं और दाएं हाथ के मूल्य एक-दूसरे से समान 3 डी दूरी (पहली स्थिति में +/- 1) होते हैं, जो समान "स्थानिक इलाके" को इंगित करते हैं। लेकिन उपरोक्त उदाहरण में आयामी क्रम (वाई, फिर जेड, फिर जेड, किसी भी विकल्प द्वारा रैखिकरण) उस इलाके को तोड़ देता है।

यह कहने का एक और तरीका यह है कि प्रारंभिक बिंदु लेने और शेष बिंदुओं को उस शुरुआती बिंदु से उनकी दूरी से क्रमशः अलग-अलग परिणाम प्रदान करेंगे। शुरुआत के रूप में 000 लेना, उदाहरण के लिए:

1D ordering : distance    3D ordering : distance
----------------------    ----------------------
        010 : 1           001,010,100 : 1
                          011,101,110 : sqrt(2)
                              111     : sqrt(3)
        011 : 2
        001 : 3
        101 : 4
        111 : 5
        110 : 6
        100 : 7

यह प्रभाव आयामों की संख्या के साथ तेजी से बढ़ता है (यह मानते हुए कि प्रत्येक आयाम का आकार "आकार" है)।


मैं अंततः टूट गया और कुछ पैसे खोल दिया। एआईपी (अमेरिकन इंस्टीट्यूट ऑफ फिजिक्स) में जॉन स्किलिंग द्वारा सी "प्रोग्रामिंग द हिल्बर्ट वक्र" में स्रोत कोड के साथ एक अच्छा, लघु लेख है (एआईपी कॉन्फ। प्रो। 707, 381 (2004) से) के लिए कोड के साथ एक परिशिष्ट है दोनों दिशाओं में मैपिंग्स। यह किसी भी आयाम> 1 के लिए काम करता है, रिकर्सिव नहीं है, राज्य-संक्रमण लुकअप टेबल का उपयोग नहीं करता है जो बड़ी मात्रा में स्मृति को पकड़ता है, और ज्यादातर बिट ऑपरेशंस का उपयोग करता है। इस प्रकार यह काफी तेज़ है और इसकी अच्छी स्मृति पदचिह्न है।

यदि आप लेख खरीदना चुनते हैं, तो मुझे स्रोत कोड में एक त्रुटि मिली।

कोड की निम्न पंक्ति (फ़ंक्शन TransposetoAxes में मिली) में त्रुटि है:

के लिए (i = n-1; i> = 0; i--) एक्स [i] ^ = एक्स [i-1];

सुधार (>) से अधिक (बराबर) (> =) से अधिक (>) को बदलने के लिए है। इस सुधार के बिना, एक्स सरणी को नकारात्मक सूचकांक का उपयोग करके एक्सेस किया जाता है जब चर "i" शून्य हो जाता है, जिससे प्रोग्राम विफल हो जाता है।

मैं आलेख पढ़ने की सिफारिश करता हूं (जो कोड सहित सात पृष्ठ लंबा है), क्योंकि यह बताता है कि एल्गोरिदम कैसे काम करता है, जो स्पष्ट से बहुत दूर है।

मैंने अपने कोड के लिए अपने कोड का अनुवाद सी # में किया। कोड निम्नानुसार है। स्किलिंग आपके द्वारा पारित वेक्टर को ओवरराइट करने के स्थान पर परिवर्तन करता है। मैंने इनपुट वेक्टर का क्लोन बनाने और एक नई प्रतिलिपि बनाने का विकल्प चुना है। इसके अलावा, मैंने विस्तार विधियों के रूप में विधियों को लागू किया।

स्किलिंग का कोड हिल्बर्ट इंडेक्स को एक सरणी के रूप में संग्रहीत एक ट्रांसपोजर के रूप में दर्शाता है। मुझे बिट्स को अंतःस्थापित करने और एक बिगइन्टर बनाने के लिए और अधिक सुविधाजनक लगता है (शब्दकोशों में अधिक उपयोगी, लूपों में फिर से भरना आसान है), लेकिन मैंने उस ऑपरेशन को अनुकूलित किया और जादू संख्याओं, बिट ऑपरेशंस और इसी तरह के विपरीत, और कोड लंबा है, इसलिए मैंने इसे छोड़ दिया है।

namespace HilbertExtensions
{
    /// <summary>
    /// Convert between Hilbert index and N-dimensional points.
    /// 
    /// The Hilbert index is expressed as an array of transposed bits. 
    /// 
    /// Example: 5 bits for each of n=3 coordinates.
    /// 15-bit Hilbert integer = A B C D E F G H I J K L M N O is stored
    /// as its Transpose                        ^
    /// X[0] = A D G J M                    X[2]|  7
    /// X[1] = B E H K N        <------->       | /X[1]
    /// X[2] = C F I L O                   axes |/
    ///        high low                         0------> X[0]
    ///        
    /// NOTE: This algorithm is derived from work done by John Skilling and published in "Programming the Hilbert curve".
    /// (c) 2004 American Institute of Physics.
    /// 
    /// </summary>
    public static class HilbertCurveTransform
    {
        /// <summary>
        /// Convert the Hilbert index into an N-dimensional point expressed as a vector of uints.
        ///
        /// Note: In Skilling's paper, this function is named TransposetoAxes.
        /// </summary>
        /// <param name="transposedIndex">The Hilbert index stored in transposed form.</param>
        /// <param name="bits">Number of bits per coordinate.</param>
        /// <returns>Coordinate vector.</returns>
        public static uint[] HilbertAxes(this uint[] transposedIndex, int bits)
        {
            var X = (uint[])transposedIndex.Clone();
            int n = X.Length; // n: Number of dimensions
            uint N = 2U << (bits - 1), P, Q, t;
            int i;
            // Gray decode by H ^ (H/2)
            t = X[n - 1] >> 1;
            // Corrected error in Skilling's paper on the following line. The appendix had i >= 0 leading to negative array index.
            for (i = n - 1; i > 0; i--) 
                X[i] ^= X[i - 1];
            X[0] ^= t;
            // Undo excess work
            for (Q = 2; Q != N; Q <<= 1)
            {
                P = Q - 1;
                for (i = n - 1; i >= 0; i--)
                    if ((X[i] & Q) != 0U)
                        X[0] ^= P; // invert
                    else
                    {
                        t = (X[0] ^ X[i]) & P;
                        X[0] ^= t;
                        X[i] ^= t;
                    }
            } // exchange
            return X;
        }

        /// <summary>
        /// Given the axes (coordinates) of a point in N-Dimensional space, find the distance to that point along the Hilbert curve.
        /// That distance will be transposed; broken into pieces and distributed into an array.
        /// 
        /// The number of dimensions is the length of the hilbertAxes array.
        ///
        /// Note: In Skilling's paper, this function is called AxestoTranspose.
        /// </summary>
        /// <param name="hilbertAxes">Point in N-space.</param>
        /// <param name="bits">Depth of the Hilbert curve. If bits is one, this is the top-level Hilbert curve.</param>
        /// <returns>The Hilbert distance (or index) as a transposed Hilbert index.</returns>
        public static uint[] HilbertIndexTransposed(this uint[] hilbertAxes, int bits)
        {
            var X = (uint[])hilbertAxes.Clone();
            var n = hilbertAxes.Length; // n: Number of dimensions
            uint M = 1U << (bits - 1), P, Q, t;
            int i;
            // Inverse undo
            for (Q = M; Q > 1; Q >>= 1)
            {
                P = Q - 1;
                for (i = 0; i < n; i++)
                    if ((X[i] & Q) != 0)
                        X[0] ^= P; // invert
                    else
                    {
                        t = (X[0] ^ X[i]) & P;
                        X[0] ^= t;
                        X[i] ^= t;
                    }
            } // exchange
            // Gray encode
            for (i = 1; i < n; i++)
                X[i] ^= X[i - 1];
            t = 0;
            for (Q = M; Q > 1; Q >>= 1)
                if ((X[n - 1] & Q)!=0)
                    t ^= Q - 1;
            for (i = 0; i < n; i++)
                X[i] ^= t;

            return X;
        }

    }
}

मैंने सी # में गिटूब में कामकाजी कोड पोस्ट किया है।

https://github.com/paulchernoch/HilbertTransformation देखें


मैंने पॉल चेर्नोक के कोड को जावा में अनुवाद करने और इसे साफ करने में थोडा समय बिताया। यह संभव है कि मेरे कोड में एक बग है, विशेष रूप से क्योंकि मेरे पास उस कागज़ तक पहुंच नहीं है जो मूल रूप से है। हालांकि, यह गुजरता है कि मैं कौन से यूनिट परीक्षण लिखने में सक्षम था। यह नीचे है।

ध्यान दें कि मैंने बड़े पैमाने पर डेटा सेट पर स्थानिक इंडेक्सिंग के लिए Z-Order और हिल्बर्ट वक्र दोनों का मूल्यांकन किया है। मुझे कहना है कि जेड-ऑर्डर बहुत बेहतर गुणवत्ता प्रदान करता है। लेकिन अपने लिए कोशिश करने के लिए स्वतंत्र महसूस करें।

    /**
     * Convert the Hilbert index into an N-dimensional point expressed as a vector of uints.
     *
     * Note: In Skilling's paper, this function is named TransposetoAxes.
     * @param transposedIndex The Hilbert index stored in transposed form.
     * @param bits Number of bits per coordinate.
     * @return Point in N-space.
     */
    static long[] HilbertAxes(final long[] transposedIndex, final int bits) {
        final long[] result = transposedIndex.clone();
        final int dims = result.length;
        grayDecode(result, dims);
        undoExcessWork(result, dims, bits);
        return result;
    }

    static void grayDecode(final long[] result, final int dims) {
        final long swap = result[dims - 1] >>> 1;
        // Corrected error in Skilling's paper on the following line. The appendix had i >= 0 leading to negative array index.
        for (int i = dims - 1; i > 0; i--)
            result[i] ^= result[i - 1];
        result[0] ^= swap;
    }

    static void undoExcessWork(final long[] result, final int dims, final int bits) {
        for (long bit = 2, n = 1; n != bits; bit <<= 1, ++n) {
            final long mask = bit - 1;
            for (int i = dims - 1; i >= 0; i--)
                if ((result[i] & bit) != 0)
                    result[0] ^= mask; // invert
                else
                    swapBits(result, mask, i);
        }
    }

    /**
     * Given the axes (coordinates) of a point in N-Dimensional space, find the distance to that point along the Hilbert curve.
     * That distance will be transposed; broken into pieces and distributed into an array.
     *
     * The number of dimensions is the length of the hilbertAxes array.
     *
     * Note: In Skilling's paper, this function is called AxestoTranspose.
     * @param hilbertAxes Point in N-space.
     * @param bits Depth of the Hilbert curve. If bits is one, this is the top-level Hilbert curve.
     * @return The Hilbert distance (or index) as a transposed Hilbert index.
     */
    static long[] HilbertIndexTransposed(final long[] hilbertAxes, final int bits) {
        final long[] result = hilbertAxes.clone();
        final int dims = hilbertAxes.length;
        final long maxBit = 1L << (bits - 1);
        inverseUndo(result, dims, maxBit);
        grayEncode(result, dims, maxBit);
        return result;
    }

    static void inverseUndo(final long[] result, final int dims, final long maxBit) {
        for (long bit = maxBit; bit != 0; bit >>>= 1) {
            final long mask = bit - 1;
            for (int i = 0; i < dims; i++)
                if ((result[i] & bit) != 0)
                    result[0] ^= mask; // invert
                else
                    swapBits(result, mask, i);
        } // exchange
    }

    static void grayEncode(final long[] result, final int dims, final long maxBit) {
        for (int i = 1; i < dims; i++)
            result[i] ^= result[i - 1];
        long mask = 0;
        for (long bit = maxBit; bit != 0; bit >>>= 1)
            if ((result[dims - 1] & bit) != 0)
                mask ^= bit - 1;
        for (int i = 0; i < dims; i++)
            result[i] ^= mask;
    }

    static void swapBits(final long[] array, final long mask, final int index) {
        final long swap = (array[0] ^ array[index]) & mask;
        array[0] ^= swap;
        array[index] ^= swap;
    }

एक और संभावना आपके डेटा पर एक kd-tree बनाने के लिए, और फिर आदेश प्राप्त करने के लिए पेड़ के इन-ऑर्डर ट्रैवर्सल को बनाना होगा। केडी-पेड़ का निर्माण करने के लिए केवल आपको एक अच्छा औसत खोज करने वाला एल्गोरिदम होना चाहिए, जिसमें से कई हैं।





dimension-reduction