Vol. III · 2026 · 본문 Entry

저불일치 수열 적용

머리

랜덤하게 나오는 듯 하지만, 잘 겹치지는 않는 그런 방식이 필요했던적이 있었습니다.

근데 찾아보니 저-불일치-순열 이라는 몇몇 정보가 나와서 그냥 작성했습니다.

아참 맨마지막 sobol로 만든 3차원 난수적용 공간을 퍼스펙으로도 찍어봤습니다.

perspectiveView.png
우주같고 이쁘지않나요?(먼지같다는 뜻)

일단 저불일치 순열도 몇개가 있는데

Sequences.gif

그래도 좀 난수같아보이는 halton과 sobol을 대상으로 정했습니다.

Halton sequence

무우우려 유니티에 포함되어있는 코드긴 합니다.

그리고 단순하게 구현할 수 있습니다.

public static float Get(int index, int radix)
{
    float result = 0f;
    float fraction = 1f / radix;

    while (index > 0)
    {
        result += (index % radix) * fraction;

        index /= radix;
        fraction /= radix;
    }

    return result;
}

테스트코드는 다음과 같이 구성했습니다.

using System.Collections;
using System.Collections.Generic;
using UnityEditor;
using UnityEngine;

public class CustomButton : DevMono
{
    private List<GameObject> instances = new List<GameObject>();

    private IEnumerator currentRoutine;

    [Button]
    void AddEditLoop()
    {
        EditorApplication.update += EditorUpdate;
    }

    private void EditorUpdate()
    {
        currentRoutine?.MoveNext();
    }

    [Button]
    private void GenerateCubesViaHalton()
    {
        currentRoutine = GenerateCubes();

        IEnumerator GenerateCubes()
        {
            for (int i = 0; i < 100; ++i)
            {
                var instance = GameObject.CreatePrimitive(PrimitiveType.Cube);

                var pos = new Vector3(Get(i, 2), Get(i, 3), Get(i, 5));
                instance.transform.position = pos;
                instance.transform.localScale = Vector3.one * 0.1f;
                instances.Add(instance);
                yield return null;
            }
        }
    }

    [Button]
    private void ReleaseEditLoop()
    {
        EditorApplication.update -= EditorUpdate;
    }

    [Button]
    private void Delete()
    {
        foreach(var instance in instances)
        {
            DestroyImmediate(instance);
        }

        instances.Clear();
    }

		//Halton
    public static float Get(int index, int radix)
    {
        float result = 0f;
        float fraction = 1f / radix;

        while (index > 0)
        {
            result += (index % radix) * fraction;

            index /= radix;
            fraction /= radix;
        }

        return result;
    }
}

Radix는 기저인데, 소수를 사용합니다.

저차원 저갯수에서는 작은 소수를 사용해도 무방합니다.

GIF 2024-03-13 오후 12-08-04.gif

Sobol Sequence

고차원을 위한 - 귀찮고, 많고, 복잡합니다.

솔직히 처음 만든게 halton이였다면, 소볼은 안건들였을 것 입니다.

internal class SobolSequence : IEnumerable<float>
{
    public class SobolEnumerator<T> : IEnumerator<float>
    {
        protected List<SobolSequence> seq;

        private SobolEnumerator() { }

        public SobolEnumerator(in SobolSequence seqeuence)
        {
            seq = new List<SobolSequence>();

            seq.Add(seqeuence);

            for (int i = 1; i < seqeuence.s; ++i)
            {
                seq.Add(new SobolSequence(seqeuence.s, seqeuence.a, seqeuence.m));
            }
        }

        public IEnumerable<float> Current => seq.Select(s => (float)s.Current);

        object IEnumerator.Current => seq.First().Current;

        float IEnumerator<float>.Current => (float)seq.First().Current;

        public bool MoveNext()
        {
            foreach (var s in seq)
            {
                s.MoveNext();
            }
            return true;
        }

        public void Reset()
        {
            foreach (var s in seq) s.Reset();
        }

        public void Dispose()
        {
            seq = null;
        }
    }

    /// <summary>
    /// 
    /// </summary>
    /// <param name="s">Dimenstion</param>
    /// <param name="a">coefficients</param>
    /// <param name="m">Seeds</param>
    /// <exception cref="InvalidOperationException"></exception>
    public SobolSequence(int s, ulong a, IList<ulong> m)
    {
        this.m = m.ToArray();
        // store the binary sequence that encodes our direction-generating recurrsion
        // we could do this with a single ulong, but the literature typically specifies
        // the degree s seperately from the coefficients of x^{s-1}, \cdots, x^2 x^1, which
        // are stored as binary digits in a = (a_{s-1} \cdots a_2 a_1)
        this.s = s;
        this.a = a;

        // now we construct the "directions"

        // The ith direction number v_i = u_i / 2^i, where u_i is an whole number with fewer
        // than i bits. 

        // we will store the v's scaled up by a constant factor 2^n, where n is the maximum i, so they
        // can be manipulated as whole binary numbers; then we will divide by 2^n once at the very end
        // when we return a floating point value
        v = new ulong[n];

        // store the scale factor 1 / 2^n
        f = 1.0 / ((double)(((ulong)1) << n));

        // the first few "directions" are given by the initialization seeds
        // the number of seed values must equal the degree of the primitive polynomial
        if (m.Count != s) throw new InvalidOperationException();
        for (int i = 0; i < s; i++)
        {
            // seeds must be odd and the ith must be less than 2^i
            if ((m[i] % 2 == 0) || (m[i] > (((ulong)1) << (i + 1)))) throw new InvalidOperationException();
            // the values are m_i / 2^i, but remember we then scale up by 2^n
            v[i] = m[i] << (n - (i + 1));
        }

        // The remainder of the "directions" are constructed via the recurrance
        //   u_i = u_{i-s} + 2^1 a_1 u_{i-1} + 2^2 a_2 u_{i-2} + \cdots + 2^{s-1} a_{s-1} u_{i-(s-1)} + 2^s u_{i-s}
        // Note + in mod 2 arithmetic is XOR.
        // Use v_i = u_i / 2^(i+1) * 2^n to translate this into a recurrence for the v's
        //   v_i = v_{i-s} 2^{-s} + a_1 v_{i-1} + a_2 v_{i-2} + \cdots + a_{s-1} v_{i-(s-1)} + v_{i-s}
        // where the a's are the binary digits of a = (a_1 a_2 \cdots a_{s-1})
        for (int i = s; i < n; i++)
        {
            // since the first and last coefficients are known to be one, we start with them
            ulong vi = v[i - s] ^ (v[i - s] >> s);
            // proceed through the bits of a, XORing in the appropriate value if the bit is 1
            for (int j = 1; j < s; j++)
            {
                if (((a >> (s - j - 1)) & 1) != 0)
                {
                    vi ^= v[i - j];
                }
            }
            v[i] = vi;
        }

    }

    // degree of primitive polynomial (i.e. degree of recurrsion)
    int s;

    // polynomial coefficients, expressed
    ulong a;

    // the number of "directions"
    // the maximum number of values in the sequence is 2^n
    private const int n = 31;

    // the "directions"
    ulong[] v;
    // the list of "initial direction numbers"
    ulong[] m;

    // the factor to convert a given value to a floating point number on (0,1)
    double f;

    // our positition in the sequence
    // p is the count
    // q stores the "state" of the system
    uint p = 0;
    ulong q = 0;

    public double Current => f * q;

    public bool MoveNext()
    {
        // get index of first zero binary digit of p = (\cdots p_2 p_1 p_0)
        int i = 0;
        while (((p >> i) & 1) != 0) i++;

        q ^= v[i];
        p++;

        return true;
    }

    public void Reset()
    {
        p = 0;
        q = 0;
    }

    IEnumerator IEnumerable.GetEnumerator()
    {
        return new SobolEnumerator<float>(this);
    }

    public IEnumerator<float> GetEnumerator()
    {
        return new SobolEnumerator<float>(this);
    }
}

//other file
public static void Generate(in int dimension, in int count, out List<float> values)
{
    var p = SobolSequenceParameters.sobolParameters[dimension];
    var seq = new SobolSequence(p.Dimension, p.Coefficients, p.Seeds);
    values = seq.Take(count).ToList();
}


//other file
internal class SobolSequenceParameters
{
    private SobolSequenceParameters(int dimension, ulong coefficients, ulong[] seeds)
    {
        this.Dimension = dimension;
        this.Coefficients = coefficients;
        this.Seeds = seeds;
    }
    public int Dimension { get; private set; }
    public ulong Coefficients { get; private set; }
    public ulong[] Seeds { get; private set; }

    //http://web.maths.unsw.edu.au/~fkuo/sobol/

    public static readonly SobolSequenceParameters[] sobolParameters = new SobolSequenceParameters[]
    {
        new SobolSequenceParameters(1, 0, new ulong[] { 1 }),
        new SobolSequenceParameters(2, 1, new ulong[] { 1, 3 }),
        new SobolSequenceParameters(3, 1, new ulong[] { 1, 3, 1 }),
        new SobolSequenceParameters(3, 2, new ulong[] { 1, 1, 1 }),
        new SobolSequenceParameters(4, 1, new ulong[] { 1, 1, 3, 3 }),
        new SobolSequenceParameters(4, 4, new ulong[] { 1, 3, 5, 13 }),
        new SobolSequenceParameters(5, 2, new ulong[] { 1, 1, 5, 5, 17 }),
        new SobolSequenceParameters(5, 4, new ulong[] { 1, 1, 5, 5, 5 }),
        new SobolSequenceParameters(5, 7, new ulong[] { 1, 1, 7, 11, 19 }),
        new SobolSequenceParameters(5, 11, new ulong[] { 1, 1, 5, 1, 1 }),
        new SobolSequenceParameters(5, 13, new ulong[] { 1, 1, 1, 3, 11 }),
        new SobolSequenceParameters(5, 14, new ulong[] { 1, 3, 5, 5, 31 }),
        new SobolSequenceParameters(6, 1, new ulong[] { 1, 3, 3, 9, 7, 49 }),
        new SobolSequenceParameters(6, 13, new ulong[] { 1, 1, 1, 15, 21, 21 }),
        new SobolSequenceParameters(6, 16, new ulong[] { 1, 3, 1, 13, 27, 49 })
    };

}
//TestCode
[Button]
private void GenerateCubesViaSobol()
{
    int SampleCount = 5000;
    SobolSequenceGenerator.Generate(7, SampleCount, out var values_x);
    SobolSequenceGenerator.Generate(9, SampleCount, out var values_y);
    SobolSequenceGenerator.Generate(8, SampleCount, out var values_z);
    currentRoutine = GenerateCubes();

    IEnumerator GenerateCubes()
    {
        for (int i = 0; i < SampleCount; ++i)
        {
            var instance = GameObject.CreatePrimitive(PrimitiveType.Cube);

            var pos = new Vector3(values_x[i], values_y[i], values_z[i]);
            instance.transform.position = pos;
            instance.transform.localScale = Vector3.one * 10f / SampleCount;
            instances.Add(instance);
            yield return null;
            yield return null;
        }
    }
}



소볼의경우 a - coefficients, m- init direction에 따라 값이 달라지는데,

샘플에 사용된 s,a,m는 다음과 같습니다.

new SobolSequenceParameters(5, 4, new ulong[] { 1, 1, 5, 5, 5 }),  //(x)
new SobolSequenceParameters(5, 7, new ulong[] { 1, 1, 7, 11, 19 }),//(z)
new SobolSequenceParameters(5, 11, new ulong[] { 1, 1, 5, 1, 1 }), //(y)
xz → xy → yz 화면으로 회전합니다
xz → xy → yz 화면으로 회전합니다

꼬리

이거 2022년도에 기능으로 구현해서 넣고 2년뒤에서야 gif 파일목록 구경하다가 찾아서 작성합니다 ㅎ….

오랜만에 또 구현하니까 재미있네요