Mathematica 模式匹配实验

结论：

利用模式匹配，可以避免循环，完成对数组中特定元素的操作，方法有多种，但是速度有所差别。

1. 在不用到特殊函数的场合，可以考虑使用Compile加速（有C编译器时优先编译为C，另外可加上并行运算指令）; 

2. 不能编译的话，优先使用Map函数；

3. 与Matlab一样，多考虑向量化编程；

4. 尽量减少中间变量；避免引入副作用。

此问题来源于知乎上的一次讨论：

「Mathematica中已知位置列表,如何直接对目标元素进行操作？」

P.S.

Matlab 的对应代码如下：

clear all

N = 40;

tic

vc = zeros(N,N,N,N);

for k1 = 0:N-1        

    for k2 = 0:N-1

        for m = 0:N-1     

            for n = 0:N-1

                c1 = cos (pi/N * (m+1/2) * k1);

                c2 = cos (pi/N * (n+1/2) * k2);

                c = c1 * c2;

                vc(k1+1,k2+1,m+1,n+1) = ((c<0)*10 + (c>0))*c;

            end

        end    

    end

end

vc;

TotalTime = toc

关于 Mathematica 的 Manipulate 函数无法释放CPU的解决办法

最近遇到一个小问题，就是在使用 Manipulate 函数时，一旦涉及到迭代、循环语句，代码就会反复执行（右侧方框加黑——而换用 Animate 函数则没有这种情况），只能手动中止，例如：
   
    Manipulate[
     date = Table[0, {num}];
     For[i = 1, i <= num, i++, date[[i]] = Sin[i]];
     ListPlot[{date}]
     , {{num, 2, "iteration num"}, 1, 20, 1}]

再来一个双层循环的例子：

    Manipulate[           
      buf = Table[0, {10}];
      Mw = Table[0, {10}];
      For[j = 1, j <= 10, j++,
       For[i = 1, i <= 10, i++,
        w = RandomVariate[NormalDistribution[0, var], 10];
        buf[[i]] = Mean[w];];
       Mw[[j]] = Mean@Cos[buf];];     
      ListLinePlot[Mw, PlotRange -> Full, DataRange -> All, Mesh -> Full,
       GridLines -> Automatic,
       GridLinesStyle -> Directive[Orange, Dashed], Frame -> True]     
     ,{{var, 3, "Var"}, 1, 5, 0.5}]
 
这两个例子还不是很严重，但是结构一旦再复杂一点， 就不太好办了，有时甚至会将前端卡死……

去网上搜索、发问，收到了一些具有启发性的回答。有人说到 Manipulate 是实时计算的，而 Animate 是统一计算所有值后再呈现结果，二者不同。

这样一来， 一旦 Manipulate 涉及到迭代、循环一类的语句，就会反复调用，无法退出，这也是逻辑上讲得通的，至少不使用迭代语句，Manipulate 就不会发生这种情况。

如果这是问题的原因，那么如何解决呢？

首先是修改算法，对于简单的问题，将循环语句转化为非循环语句或许行得通，可是对于很大一类迭代算法来说，这不是件轻易可以办到的事情，还是要从工具本身去需找答案。

于是我再去查看 Help文档，终于发现 Mathematica 的开发者早已考虑到这个问题，而且我遇到的问题也不能算作缺陷，只是一种非预设的情况罢了。

以下摘抄自「高级操作（Manipulate）功能」：

    每当步长被从零移开时，内容区会就会连续地更新，一个 CPU 监视器会指示 Mathematica 正在使用 CPU 时间. 你让这种情况持续多久，它就会持续多久.

    然而在某些情况下，持续的再计算是没有意义和不受欢迎的.
   
        例1：   
        Manipulate[
         temp = n;
         temp = temp^3;
         Graphics[{Thickness[0.01], Line[{{0, 0}, {n, temp}}]},
          PlotRange -> 1],
         {n, -1, 1}]
       
        例2：
        Manipulate[
         f[x_] := x^3;
         Graphics[{Thickness[0.01], Line[{{0, 0}, {n, f[n]}}]},
          PlotRange -> 1],
         {n, -1, 1}]
   
    在这两个情况下，这个问题都可以通过把引起问题的那些变量在一个 Module 里设成局部变量来解决. （这无论如何这都是一个好的编程习惯，远不止是为了避免无意义的更新）
   
        例1：
        Manipulate[Module[{temp},
          temp = n;
          temp = temp^3;
          Graphics[{Thickness[0.01], Line[{{0, 0}, {n, temp}}]},
           PlotRange -> 1]],
         {n, -1, 1}]
        
        例2：
        Manipulate[Module[{f},
          f[x_] := x^3;
          Graphics[{Thickness[0.01], Line[{{0, 0}, {n, f[n]}}]},
           PlotRange -> 1]],
         {n, -1, 1}]
    
    不管你对局部 Module 变量做什么都不会造成重新触发，因为这是 Module 定义的一部分，即一次调用后的变量值不会存留到下一次（所以下一轮的结果不会仅仅因为当前一轮运行中对局部变量所做的任何动作而有任何不同）.
   
    另一个解决的办法是使用 Manipulate 的 TrackedSymbols （跟踪的符号）选项来控制哪一个变量能被允许引起更新行为. 默认的值， Full （全部），意味着所有在第一个自变量中明确出现（词汇的）的符号都会被跟踪. （这意味着，除了其它事情之外，在你使用的 Manipulate 例子中函数定义里的临时变量和其它如此的问题将不会引起无限重复计算问题，这是因为它们不在第一个自变量中明显地出现，而只是通过你调用的函数间接地出现.）

    来看第二个例子，如果由于某种原因你不想要 f 成为一个局部的 Module 变量，并且你不能把它的定义移到 Manipulate 之外（在更复杂的例子中，这两种情况有些时候都是很可能会出现的），你能用 TrackedSymbols 来取消由 f 触发的更新：

        Manipulate[
         f[x_] := x^3;
         Graphics[{Thickness[0.01], Line[{{0, 0}, {n, f[n]}}]},
          PlotRange -> 1],
         {n, -1, 1}, TrackedSymbols :> {n}]
    
    这个例子只在移动滑块从而改变 n 的数值时才更新内容区域.

   

以上内容完全解决了我的疑问，我将代码稍作改动，问题搞定，如下：

    Manipulate[Module[{i},
      date = Table[0, {num}];
      For[i = 1, i <= num, i++, date[[i]] = Sin[i]];
      ListPlot[{date}]]
     , {{num, 4, "iteration num"}, 1, 20, 1}]
    
    Manipulate[           
     buf = Table[0, {10}];
     Mw = Table[0, {10}];
     For[j = 1, j <= 10, j++,
      For[i = 1, i <= 10, i++,
       w = RandomVariate[NormalDistribution[0, var], 10];
       buf[[i]] = Mean[w];];
      Mw[[j]] = Mean@Cos[buf];];
     ListLinePlot[Mw, PlotRange -> Full, DataRange -> All, Mesh -> Full,
      GridLines -> Automatic, GridLinesStyle -> Directive[Orange, Dashed],
       Frame -> True]
     , {{var, 3, "Var"}, 1, 5, 0.5}
     , TrackedSymbols :> {var}]

正如高手所说的，工具本身的 Help 是最棒的专家系统。

此外，文档这一节的末尾提到这么一段话：「具体在什么时候一个给定的动态表达式会被更新这个话题是很复杂的，这在 "动态简介" 和 "高级动态功能" 中被提到了. 在阅读那些文件时，始终注意 Manipulate 只是在 Dynamic 中把它的第一个自变量包围起来并且把它的 TrackedSymbols 选项的数值传送个在那里面的 Refresh. 所有与更新有关的事情都是由那个 Dynamic 和 Refresh 处理的.」

于是我又去看了与 Dynamic 和 Refresh 有关的内容，这里摘抄其中一条与上文有些关联的例子：

    一个可能会使人伤脑筋的情况是 RandomReal. 每一次你计算 RandomReal[]（随机实数 []），你都得到一个不同的答案，你也许就会认为 Dynamic[RandomReal[]]（动态[随机实数 []]）因此应该不停地尽快地更新自己. 但是在一般情况下这不会很有用，而且实际上会对一些在内部使用随机性的算法有负面的后果（例如，一个在 Dynamic 里面的 Monte Carlo 积分很可能不应该不停地更新，因为事实上它会每次给出一个稍微有些不同的答案）.

        Dynamic[Refresh[RandomVariate@NormalDistribution[0, 1], UpdateInterval-> 1]]

这次关于 Mathematica 动态功能的探索就到此为止了，最大的获益就是：以后再遇到什么问题，首先去搜强大的 Help 文档。

Mathematica Tips – 2

1.Module[]中注意加分号, 否则在计算时会出现意想不到的错误;而且格式化代码时会乱掉.

2.提取列表中的某一项元素出现的位置
Position[{a, b, a, a, b, c, b}, b]

3.Reap 和 Sow 用法：
In[3]:= Reap[Sow[1, {x, x}]; Sow[2, y]; Sow[4, x], {x, x, y}]
Out[3]= {4, {{{1, 1, 4}}, {{1, 1, 4}}, {{2}}}}

关于Reap和Sow以及Scan的一个巧妙用法:
In[18]:= partition[l_, v_, comp_] :=
Flatten /@
  Reap[Scan[
     Which[comp[#1, v], Sow[#1, less], comp[v, #1], Sow[#1, large],
       True, Sow[#1, equ]] &, l], {large, equ, less}][[2]]

In[19]:= partition[{3, 5, 7, 9, 2, 4, 6, 8, 3, 4}, 4, Less]
Out[19]= {{5,7,9,6,8},{4,4},{3,2,3}}

4.不要轻易加脚标:
Subscript[s, k] = Array[0 &, nK];
Subscript[s, k][[1]] = 1;

5.Map(/@) 和 Apply(@@) 对比:

Map不是替换,而是把函数映射到列表的每一个元素, 即把元素分别传递给函数的变量:
In[295]:= f /@ {1, 2, 3, 4, 5}
Out[295]= {f[1], f[2], f[3], f[4], f[5]}

Apply 自动把列表中的元素传递给多变量函数(也就是它仅仅替换函数头),@@@形式自动作用于第一层
In[52]:= Mod[#1, #2] & @@@ {{10, 4}, {5, 2}}
Out[52]= {2,1}

同样的效果, 二者结合方式(借助纯函数):
In[41]:= Apply[Mod, #] & /@ {{10, 4}, {5, 2}}
Out[41]= {2,1}

6.Evaluate函数很有用:
    In[1]:= ch = ChebyshevT[5, x]
    Out[1]= 16 x^5-20 x^3+5 x

    In[4]:= Function[x, Evaluate[ch]]
    Out[4]= x\[Function]16 x^5-20 x^3+5 x

    In[5]:= %[10]
    Out[5]= 1580050
   
   
7.Compile 函数可以提高数值运算速度! (当然采用内部函数是最快的)
它不但可以处理数学表达式，还可以处理各种简单的 Mathematica 程序. 例如，Compile 可以处理条件和控制流结构.
对比:
1.使用Compile:
In[40]:= newtonIteration :=
  Compile[{x, {n, _Integer}}, Module[{t}, t = x;
    Do[t = (t + x/t)/2, {n}]; t]];
newtonIteration[2.4, 666666] // Timing

Out[41]= {0.094,1.54919}

2.不使用Compile:
In[42]:= newtonIteration2[x_, n_] := Module[{t}, t = x;
   Do[t = (t + x/t)/2, {n}]; t];
newtonIteration2[2.4, 666666] // Timing

Out[43]= {2.137,1.54919}

运行时间相差20倍.

8.编写高效代码最重要的方式之一就是要避免显式部分引用，特别是在内部循环中.

如果将要进行实数的操作，一定要确保使用实数进行初始化.

混合的符号/数值矩阵通常比操作数值矩阵要慢.

一个整数矩阵将使用符号计算技术，其速度较慢, 但可以给出精确解。

9.条件迭代：
ff = # + 1 &;
NestWhileList[ff, 2, # < 5 &]

双层迭代的一个例子:
Mp = Table[1, {1}, {4}]
Module[{i, k},
For[i = 1, i < 5, i++,
   Mpp = {};
   Do[AppendTo[Mpp, Sin[k + i]];
    , {k, 4}];
   AppendTo[Mp, Mpp]];]
Mp

 

10.迭代例子2:

intN = 10;
intK = 200;
(*注意是矩阵*)pki = RandomInteger[{1, 2 intK}, {1, intK}];
mA = RandomInteger[{1, 2 intK}, {intK, intN}];
ak = RandomInteger[{1, 2 intK}, {intK, intN}];
y = RandomInteger[{1, 2 intK}, intN];
(*mP[list_]:=DiagonalMatrix[list];*)
(*Subscript[R, y](i)更新*)
Ryi[i_] := mA\[HermitianConjugate].(DiagonalMatrix[pki[[i]]]).mA;
Ryi[1] = RandomInteger[{1, 2 intN}, {intN, intN}];
(*Subscript[p, k](i)更新*)
pkki[k_, i_] := (Abs[ak[[k]]\[Conjugate].Inverse[Ryi[i - 1]].y])^2/(ak[[k]]\[Conjugate].Inverse[Ryi[i - 1]].ak[[k]])^2;

For[i = 1; endCond = 1, endCond > 10^-4(*i<=11*), i++,
bufMp = {};
Do[AppendTo[bufMp,(*N[ pkki[k,i+1] ]*)pkki[k, i + 1]];
  , {k, intK}];
(*AppendTo使用需要注意层次*)
AppendTo[pki, bufMp];
endCond = Norm[pki[[i + 1]] - pki[[i]]]/Norm[pki[[i]]];
]
vectorPk = pki[[i]]

Mathematica Tips - 1

    注意空格，养成书写严格形式的习惯！
   
    (*不要轻易使用MatrixForm命令*)
    (*使用Take需注意层数*)
   
    局部变量（染成绿色）：
        In[2]:= m = i^2
        Out[2]= i^2

        在块内 i+m 的计算过程中，i 用了局部值.
        In[3]:= Block[{i = a}, i + m]
        Out[3]= a + a^2

        这里仅明显出现在 i+m 中的 i 被当作局部变量处理.
        In[4]:= Module[{i = a}, i + m]
        Out[4]= a + i^2

       
    ->和:>区别：
        -> 在首次输入时计算；:> 在使用时计算：
        In[6]:= {x, x, x, x} /. x -> RandomReal[]
        Out[6]= {0.339388, 0.339388, 0.339388, 0.339388}

        In[5]:= {x, x, x, x} /. x :> RandomReal[]
        Out[5]= {0.715386, 0.330938, 0.828099, 0.762862}

   
    Mathematica将向量、矩阵均视为张量处理；与MatLab不同。
       
    数值计算的结果精度不比在变量精度的基础上判断的结果高！
   
    随机数设定种子！
   
   
    @@用来替换函数头部：
        In[2]:= f @@ {a, b, c, d}
        Out[2]= f[a, b, c, d]
       
    /@用来映射函数体
        In[1]:= Reverse /@ {{a, b}, {c, d}, {e, f}}
        Out[1]= {{b, a}, {d, c}, {f, e}}
   
    @,@@,/@区别:
        In[4]:= f@g[a, b, c, d]
        Out[4]= f[g[a, b, c, d]]
       
        In[5]:= f @@ g[a, b, c, d]
        Out[5]= f[a, b, c, d]
       
        In[6]:= f /@ g[a, b, c, d]
        Out[6]= g[f[a], f[b], f[c], f[d]]
       
   
    /.用来替换表达式中的某个变量
        In[2]:= {x, x^2, y, z} /. x -> {a, b}
        Out[2]= {{a, b}, {a^2, b^2}, y, z}
   
    f[##]&    纯函数，##：所有参量
        In[5]:= f[##] & /@ {a, b, c}等价于In[5]:= f[##] & [Sequence[{a, b, c}]]
        Out[5]= {f[a], f[b], f[c]}

        In[4]:= f[##] &[a, b, c]
        Out[4]= f[a, b, c]
   
    MapIndexed传递2个参数：[list,{list num}]
   
    表达式中可以用Print 函数输出：
    Do[If[PrimeQ[2^n - 1], Print[n]], {n, 50, 100}]
   
    Show可以显示多个图形
        Show[g1,g2,...]
       
    Nest可以重复调用函数：
        Nest[Framed, x, 6]
       
    [[]]与[[{}]]的区别，后者有「头部」
    要特别注意「层数」的概念！