2009-01-31 9 views
32

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

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

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

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

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

+0

मैंने ओलाप डेटा के बहु आयामी मैपिंग के लिए हिल्बर्ट वक्र का उपयोग किया है, मैंने पाया कि प्रदर्शन के संदर्भ में यह सरल एल्गोरिदम से बेहतर नहीं था। लेकिन मैं आपके से आयामों के छोटे सेटों के साथ परीक्षण कर रहा था। मुझे यकीन नहीं है कि आपका वास्तविक प्रश्न क्या है। – Tom

+0

अच्छा, मेरा सवाल है: "मैं यह कैसे करता हूं?" –

+0

मैंने सप्ताहों में एक ही प्रश्न के उत्तर की तलाश में बिताए हैं। इस विषय पर कागजात या तो अनजान, या पठनीय हैं लेकिन बिना किसी कोड के दिए गए हैं। जब मुझे कोड मिल जाता है, तो मैं इसका पालन नहीं कर सकता, या यह कहता है कि दृष्टिकोण दस आयामों से अधिक पैमाने पर नहीं होगा। फिर भी, विशेषज्ञों का कहना है कि आप जिस दृष्टिकोण का पीछा कर रहे हैं वह ध्वनि है, इसलिए हार मत मानो! –

उत्तर

38

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

आप लेख खरीद करने के लिए चुनते हैं, तो मैं स्रोत कोड में एक त्रुटि की खोज की। एक्स [मैं]^= एक्स [मैं के लिए (;; i> = 0 i-- i = n-1)

:

कोड (समारोह TransposetoAxes में पाया) की निम्न पंक्ति त्रुटि है -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

+0

इसे पोस्ट करने के लिए धन्यवाद। मूल और आपके कोड दोनों के लाइसेंस पर कोई स्पष्टीकरण? –

+1

जहां तक ​​मेरा संबंध है, आप उपरोक्त सी # कोड का स्वतंत्र रूप से उपयोग कर सकते हैं। मैंने जो कुछ किया वह सी से सी # में अनुवाद किया गया था और जर्नल लेख में एक टाइपो को ठीक किया गया था। कृपया स्किलिंग को क्रेडिट दें। –

+1

स्रोत कोड भी यहां पोस्ट किया गया है: http://www.tddft.org/svn/octopus/trunk/src/grid/hilbert.c (जॉन स्किलिंग नाम के अंदर उद्धृत)। –

1

मुझे नहीं लगता कि आप एक आयाम में हिल्बर्ट वक्र का उपयोग कैसे कर सकते हैं।

यदि आप दूरी को संरक्षित करते समय कम आयाम पर अंक मैप करने में रुचि रखते हैं (न्यूनतम त्रुटि के साथ) तो आप "बहुआयामी स्केलिंग" एल्गोरिदम देख सकते हैं।

नकली एनीलिंग एक दृष्टिकोण है।

संपादित करें: टिप्पणी के लिए धन्यवाद। मैं देखता हूं कि अब हिल्बर्ट वक्र दृष्टिकोण से आपका क्या मतलब है। हालांकि, यह एक कठिन समस्या है, और एन = 100 और 10 मिलियन डेटा पॉइंट दिए गए हैं, मुझे नहीं लगता कि कोई दृष्टिकोण इलाके को अच्छी तरह से संरक्षित रखेगा और उचित समय में चलाएगा। मुझे नहीं लगता कि केडी पेड़ यहां काम करेंगे।

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

+0

मैं एन-आयामी अंतरिक्ष में हिल्बर्ट वक्र का उपयोग करना चाहता हूं, न कि एक आयाम में। –

2

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

+0

आयाम के लिए लगभग 100 केडी-पेड़ उपयुक्त समाधान नहीं है। इसकी जटिलता आयामों की संख्या के साथ तेजी से बढ़ती है। – BartoszKP

4

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

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

जो हो सकता है "Hilbertized" निम्न पथ द्वारा:

001 -----> 101 
    \   \ 
    \   \ 
    011  111 
    ^  | 
    |   | 
000 |  100 | 
    \ |  \ | 
    \ |  \ V 
    010  110 
-1 डी क्रम में

:

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

यहाँ बुरा सा इस trival 3 डी मामले पर विचार करें। जोड़े की सूची पर विचार करें और -1 डी नीचे दूरी:

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

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

यह कहने का एक और तरीका यह है कि शुरुआती बिंदु से अपनी दूरी से शेष बिंदुओं को क्रमबद्ध करने और क्रमशः क्रमशः अलग-अलग परिणाम प्रदान करेंगे। 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 

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

+0

ठीक है, अगर मुझे यह सही लगता है, तो इलाके अभी भी संरक्षित है - घन के भीतर सभी बिंदु उस घन के किनारों से और कहीं और नहीं होंगे। कम घन मात्रा है, बेहतर इलाके संरक्षित है। मुझे लगता है कि परिणाम मेरे लिए काफी करीब होना चाहिए। –

+0

"हिल्बर्ट स्पेस-फिलिंग वक्र के क्लस्टरिंग प्रॉपर्टीज का विश्लेषण" लेख भी देखें: http://www.cs.cmu.edu/afs/cs.cmu.edu/user/christos/www/PUBLICATIONS/ieee-tkde -hilbert.ps.gz –

+0

हाँ, लेकिन जैसे ही आपके "क्यूब्स" छोटे हो जाते हैं, आपके हिल्बर्ट वक्र की लंबाई तेजी से बढ़ जाती है। मुझे संदेह है कि आप जिस समस्या को हल करने की कोशिश कर रहे हैं उसकी जटिलता के कारण, आपको अपने अंक पर अच्छा ऑर्डर करने से बहुत लंबा समय लगेगा। – Imran

8

से n-> 1 और 1-> मानचित्रण के लिए एल्गोरिथ्म n यहाँ "Calculation of Mappings Between One and n-dimensional Values Using the Hilbert Space-filling Curve" J K Lawder

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

+0

मैंने आपके द्वारा इंगित लेख को पूरी तरह से चलाया, लेकिन, केवल तालिका 3 को देखकर, मैं एल्गोरिदम को नहीं समझता जो वे 'n' से' 1' तक जाने का सुझाव देते हैं। क्या यह पुनरावृत्त है? वे 'w' से शुरू होते हैं, जो 'tau_tilde' पर निर्भर करता है जिसे अंत में गणना की जाती है। – user3285148

+0

यदि हां, तो मुझे कैसे प्रारंभ करना चाहिए? – user3285148

2

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

ध्यान दें कि मैंने 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; 
    } 
+0

क्या आपने शायद लंबे [] से BigInteger में रूपांतरण भी किया था? (अप्रचलित) मुझे रेंज पूछताछ में रूचि है। –

संबंधित मुद्दे