Issue with threading matrices

Question

I like this answer, which provides an "autoThread" function for more intelligent threading, but it doesn't work in this example and was hoping for a fix.

So for this code:

Format[\[Gamma][i_, j_]] := Subscript[\[Gamma], {i, j}]
dim = 5;
decayDecoherenceMatrix = 
  Table[ 1/2 ([Gamma][i, i] + [Gamma][j, j]), {i, dim}, {j, dim}];
dephasingMatrix = ConstantArray[0, {5, 5}];
dephasingMatrix[[1, 5]] = [Gamma]d;
dephasingMatrix[[5, 1]] = [Gamma]d;
(the ground states are 1 and 5 in this case!)
decayRules = {[Gamma][1, 1] -> 
   0, [Gamma][2, 2] -> [CapitalGamma], [Gamma][3, 
    3] -> [CapitalGamma], [Gamma][4, 
    4] -> [CapitalGamma], [Gamma][5, 5] -> 0} 
decayResult = (decayDecoherenceMatrix + dephasingMatrix) /. decayRules
decayMatrix = Table[ [Gamma][i, j], {i, dim}, {j, dim}];
autoThread[decayMatrix -> decayResult]

Using the autoThread function provided in the link:

SetAttributes[{step, $stepHold}, HoldAll]
step[expr_] := 
 Module[{P}, P = (P = Return[$stepHold @@ #, TraceScan] &) &;
  TraceScan[P, expr, TraceDepth -> 1]]
Attributes[autoThread] = HoldFirst;
autoThread[body_] := autoThread[body, List]
autoThread[body : [___, h[___], ___], h_, seq_: All] := 
 Thread[Unevaluated@body, h, seq]
autoThread[body : f_[arg___], h_, seq_: All] := 
 With[{new = 
    Replace[MapAll[step, Unevaluated[body], 
      Heads -> True], (step | $stepHold)[x_] :> x, -1, 
     Heads -> True]}, (new /. $stepHold[eval_] :> 
      autoThread[eval, h, seq]) /; new =!= $stepHold[body]]
autoThread[else_, x___] := 
 step[else] /. $stepHold[eval_] :> autoThread[eval, x]
complex[expr_] := 
 Block[{w, num, dem}, w = Together /@ ComplexExpand@expr;
  num = Numerator /@ w;
  dem = Simplify@Denominator /@ w[[1]];
  num/dem]

Produces an output of the form:

{{
\!\(\*SubscriptBox[\(\[Gamma]\), \({1, 1}\)]\), 
\!\(\*SubscriptBox[\(\[Gamma]\), \({1, 2}\)]\), 
\!\(\*SubscriptBox[\(\[Gamma]\), \({1, 3}\)]\), 
\!\(\*SubscriptBox[\(\[Gamma]\), \({1, 4}\)]\), 
\!\(\*SubscriptBox[\(\[Gamma]\), \({1, 
     5}\)]\)} -> {0, \[CapitalGamma]/2, \[CapitalGamma]/
   2, \[CapitalGamma]/2, \[Gamma]d}, {
\!\(\*SubscriptBox[\(\[Gamma]\), \({2, 1}\)]\), 
\!\(\*SubscriptBox[\(\[Gamma]\), \({2, 2}\)]\), 
\!\(\*SubscriptBox[\(\[Gamma]\), \({2, 3}\)]\), 
\!\(\*SubscriptBox[\(\[Gamma]\), \({2, 4}\)]\), 
\!\(\*SubscriptBox[\(\[Gamma]\), \({2, 5}\)]\)} -> {\[CapitalGamma]/
   2, \[CapitalGamma], \[CapitalGamma], \[CapitalGamma], \
\[CapitalGamma]/2}, {
\!\(\*SubscriptBox[\(\[Gamma]\), \({3, 1}\)]\), 
\!\(\*SubscriptBox[\(\[Gamma]\), \({3, 2}\)]\), 
\!\(\*SubscriptBox[\(\[Gamma]\), \({3, 3}\)]\), 
\!\(\*SubscriptBox[\(\[Gamma]\), \({3, 4}\)]\), 
\!\(\*SubscriptBox[\(\[Gamma]\), \({3, 5}\)]\)} -> {\[CapitalGamma]/
   2, \[CapitalGamma], \[CapitalGamma], \[CapitalGamma], \
\[CapitalGamma]/2}, {
\!\(\*SubscriptBox[\(\[Gamma]\), \({4, 1}\)]\), 
\!\(\*SubscriptBox[\(\[Gamma]\), \({4, 2}\)]\), 
\!\(\*SubscriptBox[\(\[Gamma]\), \({4, 3}\)]\), 
\!\(\*SubscriptBox[\(\[Gamma]\), \({4, 4}\)]\), 
\!\(\*SubscriptBox[\(\[Gamma]\), \({4, 5}\)]\)} -> {\[CapitalGamma]/
   2, \[CapitalGamma], \[CapitalGamma], \[CapitalGamma], \
\[CapitalGamma]/2}, {
\!\(\*SubscriptBox[\(\[Gamma]\), \({5, 1}\)]\), 
\!\(\*SubscriptBox[\(\[Gamma]\), \({5, 2}\)]\), 
\!\(\*SubscriptBox[\(\[Gamma]\), \({5, 3}\)]\), 
\!\(\*SubscriptBox[\(\[Gamma]\), \({5, 4}\)]\), 
\!\(\*SubscriptBox[\(\[Gamma]\), \({5, 
     5}\)]\)} -> {\[Gamma]d, \[CapitalGamma]/2, \[CapitalGamma]/
   2, \[CapitalGamma]/2, 0}}

(Sorry I can't remember the proper way to format this output for Mathematica!)

I thought that this autoThread function was supposed to automatically thread through at any depth, but maybe I was wrong? Can autoThread be modified to be able to cover this kind of case automatically?

I know I can fix my code by hand:

decoherenceRules = 
  Thread[autoThread[decayMatrix -> decayResult][[#]]] & /@ Range[dim];

But I was hoping for something that could handle this automatically.

Answer 1

If I am understanding correctly, to thread tensors at a given level you can use MapThread. For your example this would be MapThread[#1 -> #2 &, {decayMatrix, decayResult}, 2]. I don't know that this is necessarily related to the issue you linked.