wolfram-mathematica - 从稀疏定义列表中选择无模式下调值的算法

Question

我有以下问题。

我正在开发一个随机模拟器，该模拟器可以随机采样系统的配置，并存储在特定时间实例中每种配置被访问了多少次的统计信息。大致来说，代码是这样的

f[_Integer][{_Integer..}] :=0
...
someplace later in the code, e.g.,
index = get index;
c = get random configuration (i.e. a tuple of integers, say a pair {n1, n2}); 
f[index][c] = f[index][c] + 1;
which tags that configuration c has occurred once more in the simulation at time instance index.

代码完成后，将出现一个f的定义列表，它看起来像这样(我手动输入只是为了强调最重要的部分)

?f
f[1][{1, 2}] = 112
f[1][{3, 4}] = 114
f[2][{1, 6}] = 216
f[2][{2, 7}] = 227
...
f[index][someconfiguration] = some value
...
f[_Integer][{_Integer..}] :=0

请注意，最先出现的无模式定义可能很少。同样，人们也不知道将选择哪些值和配置。

问题是要有效地提取所需索引的下降值，例如发出类似

result = ExtractConfigurationsAndOccurences[f, 2] 

这应该给出具有结构的列表

result = {list1, list2}

在哪里

list1 = {{1, 6}, {2, 7}} (* the list of configurations that occurred during the simulation*)
list2 = {216, 227} (* how many times each of them occurred *)

问题在于ExtractConfigurationsAndOccurences应该非常快。我唯一能想到的解决方案是使用SubValues [f](提供完整列表)，并使用Cases语句对其进行过滤。我意识到应该不惜一切代价避免这种程序，因为将有成倍地测试许多配置(定义)，这大大降低了代码的速度。

Mathematica中有一种自然的方法可以快速地做到这一点吗？

我希望Mathematica可以将f [2]视为具有多个下降值的单个头部，但是使用DownValues [f [2]]却没有任何效果。同样使用SubValues [f [2]]会导致错误。

Answer 1

这是我先前答案的完整重写。事实证明，在我以前的尝试中，我忽略了基于压缩数组和稀疏数组的组合的简单得多的方法，该方法比所有以前的方法都更快，内存更多-效率更高(至少在样本大小范围内，进行了测试)，而仅以最小的方式更改了基于SubValues的原始方法。由于询问了有关最有效方法的问题，因此我将从答案中删除其他方法(因为它们更加复杂并且占用了大量空间。希望看到它们的人可以查看此方法的以前的修订版)回答)。
原始的基于SubValues的方法
我们首先介绍一个函数来为我们生成配置的测试样本。这里是:

Clear[generateConfigurations];
generateConfigurations[maxIndex_Integer, maxConfX_Integer, maxConfY_Integer, 
  nconfs_Integer] :=
Transpose[{
  RandomInteger[{1, maxIndex}, nconfs],
  Transpose[{
     RandomInteger[{1, maxConfX}, nconfs],
     RandomInteger[{1, maxConfY}, nconfs]
  }]}]; 

我们可以生成一个小样本来说明:

In[3]:= sample  = generateConfigurations[2,2,2,10]
Out[3]= {{2,{2,1}},{2,{1,1}},{1,{2,1}},{1,{1,2}},{1,{1,2}},
          {1,{2,1}},{2,{1,2}},{2,{2,2}},{1,{2,2}},{1,{2,1}}}

在这里，我们只有2个索引，以及“x”和“y”数字都从1变为2的配置-10个此类配置。
以下函数将帮助我们模仿配置的频率累积，因为我们为重复出现的计数器增加了基于SubValues的计数器:

Clear[testAccumulate];
testAccumulate[ff_Symbol, data_] :=
  Module[{},
   ClearAll[ff];
   ff[_][_] = 0;
   Do[
     doSomeStuff;
     ff[#1][#2]++ & @@ elem;
     doSomeMoreStaff;
   , {elem, data}]];

此处的doSomeStuff和doSomeMoreStaff符号表示一些可能在计数代码之前或之后的代码。 data参数应该是generateConfigurations产生的形式的列表。例如:

In[6]:= 
testAccumulate[ff,sample];
SubValues[ff]

Out[7]= {HoldPattern[ff[1][{1,2}]]:>2,HoldPattern[ff[1][{2,1}]]:>3,
   HoldPattern[ff[1][{2,2}]]:>1,HoldPattern[ff[2][{1,1}]]:>1,
   HoldPattern[ff[2][{1,2}]]:>1,HoldPattern[ff[2][{2,1}]]:>1,
   HoldPattern[ff[2][{2,2}]]:>1,HoldPattern[ff[_][_]]:>0}

以下函数将从SubValues列表中提取结果数据(索引，配置及其频率):

Clear[getResultingData];
getResultingData[f_Symbol] :=
   Transpose[{#[[All, 1, 1, 0, 1]], #[[All, 1, 1, 1]], #[[All, 2]]}] &@
        Most@SubValues[f, Sort -> False];

例如:

In[10]:= result = getResultingData[ff]
Out[10]= {{2,{2,1},1},{2,{1,1},1},{1,{2,1},3},{1,{1,2},2},{2,{1,2},1},
{2,{2,2},1},{1,{2,2},1}}

为了完成数据处理周期，以下是一个简单的函数，它基于Select为固定索引提取数据:

Clear[getResultsForFixedIndex];
getResultsForFixedIndex[data_, index_] := 
  If[# === {}, {}, Transpose[#]] &[
    Select[data, First@# == index &][[All, {2, 3}]]];

对于我们的测试示例，

In[13]:= getResultsForFixedIndex[result,1]
Out[13]= {{{2,1},{1,2},{2,2}},{3,2,1}}

据推测，这与@zorank在代码中尝试过的内容很接近。
基于压缩数组和稀疏数组的更快解决方案
正如@zorank所指出的，对于具有更多索引和配置的较大样本，这变得很慢。我们现在将生成一个大样本来说明这一点(注意!这需要大约4-5 Gb的RAM，因此，如果超出可用的RAM ，则可能需要减少配置数量):

In[14]:= 
largeSample = generateConfigurations[20,500,500,5000000];
testAccumulate[ff,largeSample];//Timing

Out[15]= {31.89,Null}

现在，我们将从SubValues的ff中提取完整数据:

In[16]:= (largeres = getResultingData[ff]); // Timing
Out[16]= {10.844, Null}

这需要一些时间，但是一个操作只能执行一次。但是，当我们开始为固定索引提取数据时，我们发现它非常慢:

In[24]:= getResultsForFixedIndex[largeres,10]//Short//Timing
Out[24]= {2.687,{{{196,26},{53,36},{360,43},{104,144},<<157674>>,{31,305},{240,291},
 {256,38},{352,469}},{<<1>>}}}

我们将在此处使用的主要思想是加快速度，将单个列表打包在largeres中，这些列表用于索引，组合和频率。虽然无法打包完整列表，但这些部分可以单独:

In[18]:= Timing[
   subIndicesPacked = Developer`ToPackedArray[largeres[[All,1]]];
   subCombsPacked =  Developer`ToPackedArray[largeres[[All,2]]];
   subFreqsPacked =  Developer`ToPackedArray[largeres[[All,3]]];
]
Out[18]= {1.672,Null}

这也需要一些时间，但它又是一次性操作。
然后，将使用以下函数来更有效地提取固定索引的结果:

Clear[extractPositionFromSparseArray];
extractPositionFromSparseArray[HoldPattern[SparseArray[u___]]] := {u}[[4, 2, 2]]

Clear[getCombinationsAndFrequenciesForIndex];
getCombinationsAndFrequenciesForIndex[packedIndices_, packedCombs_, 
    packedFreqs_, index_Integer] :=
With[{positions = 
         extractPositionFromSparseArray[
               SparseArray[1 - Unitize[packedIndices - index]]]},
  {Extract[packedCombs, positions],Extract[packedFreqs, positions]}];

现在，我们有:

In[25]:=  
getCombinationsAndFrequenciesForIndex[subIndicesPacked,subCombsPacked,subFreqsPacked,10]
  //Short//Timing

Out[25]= {0.094,{{{196,26},{53,36},{360,43},{104,144},<<157674>>,{31,305},{240,291},
{256,38},{352,469}},{<<1>>}}}

我们的速度提高了30倍天真的Select方法。
有关复杂性的一些注意事项
请注意，第二种解决方案速度更快，因为它使用了优化的数据结构，但是其复杂性与基于Select的复杂性相同，即基于所有索引的唯一组合的总列表长度是线性的。因此，从理论上讲，先前讨论的基于嵌套哈希表等的解决方案可能在渐近性上更好。问题是，实际上我们可能会在此之前很久就达到内存限制。对于1000万个配置示例，上述代码仍比我之前发布的最快解决方案快2-3倍。
编辑
进行以下修改:

Clear[getCombinationsAndFrequenciesForIndex];
getCombinationsAndFrequenciesForIndex[packedIndices_, packedCombs_, 
    packedFreqs_, index_Integer] :=
 With[{positions =  
          extractPositionFromSparseArray[
             SparseArray[Unitize[packedIndices - index], Automatic, 1]]},
    {Extract[packedCombs, positions], Extract[packedFreqs, positions]}];

使代码仍然快两倍。而且，对于更稀疏的索引(例如，使用诸如generateConfigurations[2000, 500, 500, 5000000]之类的参数调用样本生成函数)，相对于基于Select的函数而言，其提速约为100倍。