如何在Mathematica中找到（重复）列表的周期？

Question

如何在Mathematica中找到（重复）列表的周期？

19

什么是在重复列表中找到周期的最佳方法？

例如：

a = {4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2}

有重复项 {4, 5, 1, 2, 3}，剩余部分为 {4, 5, 1, 2}，它们匹配但不完整。

该算法应足够快以处理更长的情况，例如：

b = RandomInteger[10000, {100}];
a = Join[b, b, b, b, Take[b, 27]]

如果没有像上面那样的重复模式，该算法应返回$Failed。

- Arnoud Buzing

我的第一种方法是使用FindSequenceFunction[a]，但是我无法找到一个好的编程方法来确定它的周期。 - Brett Champion

感谢迄今为止的回复。看起来@Szabolcs的解决方案是最好的。我仍在努力开发自己独立的解决方案，但速度仍然慢了一倍。 - Arnoud Buzing

维基百科上有一篇关于“循环检测”的文章，我认为它很相关。 - Mohammad Alaggan

相关链接：http://mathematica.stackexchange.com/q/119907/280 - Alexey Popkov

9个回答

7

如果您的信号没有噪声，上述方法是更好的选择。如果您的信号只是大致的话，傅里叶变换方法可能会有用。我将通过一个“参数化”的设置进行说明，在此设置中，基本信号的长度和重复次数、尾部长度以及噪声扰动的上限都是可以调整的变量。

noise = 20;
extra = 40;
baselen = 103;
base = RandomInteger[10000, {baselen}];
repeat = 5;
signal = Flatten[Join[ConstantArray[base, repeat], Take[base, extra]]];
noisysignal = signal + RandomInteger[{-noise, noise}, Length[signal]];

我们计算FFT的绝对值。我们在两端添加零。目的是通过与相邻值进行比较进行阈值处理。

sigfft = Join[{0.}, Abs[Fourier[noisysignal]], {0}];

现在我们创建两个0-1向量。其中一个是通过对fft中每个元素大于其两个邻居的几何平均值的两倍进行阈值处理，使其为1。另一个向量使用平均值（算术平均值），但我们将大小限制降低到3/4。这是基于一些实验的结果。我们计算每种情况下1的数量。理想情况下，我们应该得到每个100，因为在没有噪声和尾部部分的“完美”情况下，这将是非零元素的数量。

In[419]:= 
thresh1 = 
  Table[If[sigfft[[j]]^2 > 2*sigfft[[j - 1]]*sigfft[[j + 1]], 1, 
    0], {j, 2, Length[sigfft] - 1}];
count1 = Count[thresh1, 1]
thresh2 = 
  Table[If[sigfft[[j]] > 3/4*(sigfft[[j - 1]] + sigfft[[j + 1]]), 1, 
    0], {j, 2, Length[sigfft] - 1}];
count2 = Count[thresh2, 1]

Out[420]= 114

Out[422]= 100

现在我们可以通过将总长度除以计数的平均值并向下取整来得出“repeats”的最佳估计值。

approxrepeats = Floor[2*Length[signal]/(count1 + count2)]
Out[423]= 5

我们发现基本信号重复了5次。这可以帮助我们开始细化估计正确的长度（上面的baselen）。为此，我们可以尝试删除末尾的元素，并观察我们何时得到更接近于实际拥有四个0之间非零值的运行的ffts。

另一种估计重复次数的方法是找到阈值fft的运行长度编码中零的模态数量。虽然我没有真正尝试过，但它看起来可能对如何进行阈值处理的细节的不良选择具有鲁棒性（我的实验似乎有效）。

Daniel Lichtblau

- Daniel Lichtblau

1

这也是我本能想到的做法，但无法可靠地使其工作。所以+1。 - acl

5

以下假设循环从第一个元素开始，提供周期长度和循环信息。

findCyclingList[a_?VectorQ] :=
  Module[{repeats1, repeats2, cl, cLs, vec}, 
  repeats1 = Flatten@Differences[Position[a, First[a]]];
  repeats2 = Flatten[Position[repeats1, First[repeats1]]]; 
  If[Equal @@ Differences[repeats2] && Length[repeats2] > 2(* 
   is potentially cyclic - first element appears cyclically *),
   cl = Plus @@@ Partition[repeats1, First[Differences[repeats2]]];
   cLs = Partition[a, First[cl]];
   If[SameQ @@ cLs  (* candidate cycles all actually the same *), 
    vec = First[cLs];
    {Length[vec], vec}, $Failed], $Failed]  ]

测试

b = RandomInteger[50, {100}];
a = Join[b, b, b, b, Take[b, 27]];

findCyclingList[a]

{100, {47, 15, 42, 10, 14, 29, 12, 29, 11, 37, 6, 19, 14, 50, 4, 38, 
  23, 3, 41, 39, 41, 17, 32, 8, 18, 37, 5, 45, 38, 8, 39, 9, 26, 33, 
  40, 50, 0, 45, 1, 48, 32, 37, 15, 37, 49, 16, 27, 36, 11, 16, 4, 28,
   31, 46, 30, 24, 30, 3, 32, 31, 31, 0, 32, 35, 47, 44, 7, 21, 1, 22,
   43, 13, 44, 35, 29, 38, 31, 31, 17, 37, 49, 22, 15, 28, 21, 8, 31, 
  42, 26, 33, 1, 47, 26, 1, 37, 22, 40, 27, 27, 16}}

b1 = RandomInteger[10000, {100}]; 
a1 = Join[b1, b1, b1, b1, Take[b1, 23]];

findCyclingList[a1]

{100, {1281, 5325, 8435, 7505, 1355, 857, 2597, 8807, 1095, 4203, 
  3718, 3501, 7054, 4620, 6359, 1624, 6115, 8567, 4030, 5029, 6515, 
  5921, 4875, 2677, 6776, 2468, 7983, 4750, 7609, 9471, 1328, 7830, 
  2241, 4859, 9289, 6294, 7259, 4693, 7188, 2038, 3994, 1907, 2389, 
  6622, 4758, 3171, 1746, 2254, 556, 3010, 1814, 4782, 3849, 6695, 
  4316, 1548, 3824, 5094, 8161, 8423, 8765, 1134, 7442, 8218, 5429, 
  7255, 4131, 9474, 6016, 2438, 403, 6783, 4217, 7452, 2418, 9744, 
  6405, 8757, 9666, 4035, 7833, 2657, 7432, 3066, 9081, 9523, 3284, 
  3661, 1947, 3619, 2550, 4950, 1537, 2772, 5432, 6517, 6142, 9774, 
  1289, 6352}}

这个案例应该会失败，因为它不是循环的。

findCyclingList[Join[b, Take[b, 11], b]]

$Failed

我尝试使用Repeated做一些事情，例如a /. Repeated[t__, {2, 100}] -> {t}，但是它对我来说根本不起作用。

- Verbeia

如果第一个If语句中的检查导致一个周期内存在第一个元素的多个出现，则会导致返回$Failed，例如：{1, 2, 3, 1, 2, 1, 2, 3, 1, 2, 1, 2, 3, 1, 2, 1, 2, 3, 1, 2}。我认为你的检查需要是差分列表本身就是周期性的。 - Brett Champion

哦，好观点 - 我现在不确定如何解决这个问题。 - Verbeia

请将以下有关编程的内容从英语翻译成中文。仅返回翻译后的文本：此修订版本应更好。 - Verbeia

@Verbeia a /. {Repeated[k__, {2, Infinity}], Shortest[f___]} -> {{k}} 可以使用，但速度太慢了。 - Dr. belisarius

3

这个对你有效吗？

period[a_] := 
   Quiet[Check[
      First[Cases[
         Table[
            {k, Equal @@ Partition[a, k]}, 
            {k, Floor[Length[a]/2]}], 
         {k_, True} :> k
         ]], 
      $Failed]]

严格来说，这种做法对于一些事情会失败，比如：

a = {1, 2, 3, 1, 2, 3, 1, 2, 3, 4, 5}

尽管可以通过使用以下内容进行修复：

(Equal @@ Partition[a, k]) && (Equal @@ Partition[Reverse[a], k])

（可能是提前计算一次Reverse [a]。）

- Brett Champion

2

我提出如下建议，它借鉴了Verbeia和Brett的答案。

Do[
  If[MatchQ @@ Equal @@ Partition[#, i, i, 1, _], Return @@ i],
  {i, #[[ 2 ;; Floor[Length@#/2] ]] ~Position~ First@#}
] /. Null -> $Failed &

这个函数在长周期上的效率不如Vebeia的函数，但在短周期上更快，而且也更简单。

- Mr.Wizard

+1 是为了创造性的周期测试函数、简单性和（潜在的）速度。你的瓶颈实际上是周期测试。如果使用更快的方法（像我使用的丑陋方法），这也是最快的解决方案。我显然过于复杂化了事情 :-) - Szabolcs

@Szabolcs，优秀的东西不可能丑陋。（别介意其他的，我现在明白了。） - Mr.Wizard

1

我不知道如何在Mathematica中解决它，但以下算法（用Python编写）应该可以工作。它的时间复杂度为O(n)，因此速度不应该是问题。

def period(array):
    if len(array) == 0:
        return False
    else:
        s = array[0]
        match = False
        end = 0
        i = 0

        for k in range(1,len(array)):
            c = array[k]

            if not match:
                if c == s:
                    i = 1
                    match = True
                    end = k
            else:
                if not c == array[i]:
                    match = False
                i += 1

        if match:
            return array[:end]
        else:
            return False

# False         
print(period([4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2,1]))
# [4, 5, 1, 2, 3]            
print(period([4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2]))
# False
print(period([4]))
# [4, 2]
print(period([4,2,4]))
# False
print(period([4,2,1]))
# False
print(period([]))

- Konstantin Weitz

1

好的，这里展示一下我的工作：

ModifiedTortoiseHare[a_List] := Module[{counter, tortoise, hare},
Quiet[
 Check[
  counter = 1;
  tortoise = a[[counter]];
  hare = a[[2 counter]];
  While[(tortoise != hare) || (a[[counter ;; 2 counter - 1]] != a[[2 counter ;; 3 counter - 1]]),
   counter++;
   tortoise = a[[counter]];
   hare = a[[2 counter]];
  ];
 counter,
$Failed]]]

我不确定这是否是100％正确的，特别是在像{pattern，pattern，different，pattern，pattern}这样的情况下，并且当有很多重复元素时，速度会变得越来越慢，就像这样：

{ 1,2,1,1, 1,2,1,1, 1,2,1,1, ...}

因为它进行了太多昂贵的比较。

- Arnoud Buzing

0

#include <iostream>
#include <vector>
using namespace std;


int period(vector<int> v)
{
    int p=0; // period 0

    for(int i=p+1; i<v.size(); i++)
    {

        if(v[i] == v[0])
        {
            p=i; // new potential period


            bool periodical=true;
            for(int i=0; i<v.size()-p; i++)
            {
                if(v[i]!=v[i+p])
                {
                    periodical=false;
                    break;
                }
            }
            if(periodical) return p;


            i=p; // try to find new period
        }
    }

    return 0; // no period

}


int main()
{

    vector<int> v3{1,2,3,1,2,3,1,2,3};
    cout<<"Period is :\t"<<period(v3)<<endl;


    vector<int> v0{1,2,3,1,2,3,1,9,6};
    cout<<"Period is :\t"<<period(v0)<<endl;

    vector<int> v1{1,2,1,1,7,1,2,1,1,7,1,2,1,1};
    cout<<"Period is :\t"<<period(v1)<<endl;
    return 0;
}

- Fahrudin Halilović

并提供一些解释 - HaveNoDisplayName

-1

这听起来可能与序列比对有关。这些算法已经得到了深入的研究，可能已经在Mathematica中实现。

- Gian

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- Szabolcs · Accepted Answer

请参见代码中夹杂的注释，了解它的工作原理。

(* True if a has period p *)
testPeriod[p_, a_] := Drop[a, p] === Drop[a, -p]

(* are all the list elements the same? *)
homogeneousQ[list_List] := Length@Tally[list] === 1
homogeneousQ[{}] := Throw[$Failed] (* yes, it's ugly to put this here ... *)

(* auxiliary for findPeriodOfFirstElement[] *)
reduce[a_] := Differences@Flatten@Position[a, First[a], {1}]

(* the first element occurs every ?th position ? *)
findPeriodOfFirstElement[a_] := Module[{nl},
  nl = NestWhileList[reduce, reduce[a], ! homogeneousQ[#] &];
  Fold[Total@Take[#2, #1] &, 1, Reverse[nl]]
  ]

(* the period must be a multiple of the period of the first element *)
period[a_] := Catch@With[{fp = findPeriodOfFirstElement[a]},
   Do[
    If[testPeriod[p, a], Return[p]],
    {p, fp, Quotient[Length[a], 2], fp}
    ]
   ]

如果findPeriodOfFirstElement[]不清楚，请问。我独立完成了这个（为了好玩！），但现在我看到原则与Verbeia的解决方案相同，除了Brett指出的问题已经修复。

我正在进行测试

b = RandomInteger[100, {1000}];
a = Flatten[{ConstantArray[b, 1000], Take[b, 27]}];

（注意低整数值：同一周期内会有很多重复的元素*）

编辑:根据下面Leonid的评论，可以通过使用专门编译用于整数列表的自定义位置函数再次加速2-3倍（在我的计算机上大约是2.4倍）。

(* Leonid's reduce[] *)

myPosition = Compile[
  {{lst, _Integer, 1}, {val, _Integer}}, 
  Module[{pos = Table[0, {Length[lst]}], i = 1, ctr = 0}, 
    For[i = 1, i <= Length[lst], i++, 
      If[lst[[i]] == val, pos[[++ctr]] = i]
    ]; 
    Take[pos, ctr]
  ], 
  CompilationTarget -> "C", RuntimeOptions -> "Speed"
]

reduce[a_] := Differences@myPosition[a, First[a]]

编译testPeriod在快速测试中可以进一步提高约20%的速度，但我认为这将取决于输入数据：

Clear[testPeriod]
testPeriod = 
 Compile[{{p, _Integer}, {a, _Integer, 1}}, 
  Drop[a, p] === Drop[a, -p]]