Tricky selection of elements based on their position

I need to make a list of some elements from a 4D array with dimensions 4x4x4x4. I need to select the elements based on their position in the following way: name an element e(x,y,z,w); I need all the elements where x≥y and z≥w. This basically excludes all the upper triangular part of each submatrix and all the submatrices in the upper triangular part of the tensor. The array in question is this one:

b={{{{12, 0, 0, 0}, {-18 Sqrt[2], 0, 0, 0}, {12 Sqrt[3], 0, 0, 0}, {-6,
0, 0, 0}}, {{0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0,
0}}, {{0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0,
0}}, {{0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0,
0}}}, {{{-18 Sqrt[2], 0, 0, 0}, {84, 0, 0, 0}, {-48 Sqrt[6], 0, 0,
0}, {54 Sqrt[2], 0, 0, 0}}, {{0, 0, 0, 0}, {0, 60, 0,
0}, {0, -90, 0, 0}, {0, 18 Sqrt[10], 0, 0}}, {{0, 0, 0, 0}, {0, 0,
0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}}, {{0, 0, 0, 0}, {0, 0, 0,
0}, {0, 0, 0, 0}, {0, 0, 0, 0}}}, {{{12 Sqrt[3], 0, 0,
0}, {-48 Sqrt[6], 0, 0, 0}, {276, 0, 0, 0}, {-186 Sqrt[3], 0, 0,
0}}, {{0, 0, 0, 0}, {0, -90, 0, 0}, {0, 240, 0,
0}, {0, -90 Sqrt[10], 0, 0}}, {{0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0,
168, 0}, {0, 0, -84 Sqrt[6], 0}}, {{0, 0, 0, 0}, {0, 0, 0, 0}, {0,
0, 0, 0}, {0, 0, 0, 0}}}, {{{-6, 0, 0, 0}, {54 Sqrt[2], 0, 0,
0}, {-186 Sqrt[3], 0, 0, 0}, {648, 0, 0, 0}}, {{0, 0, 0, 0}, {0,
18 Sqrt[10], 0, 0}, {0, -90 Sqrt[10], 0, 0}, {0, 600, 0, 0}}, {{0,
0, 0, 0}, {0, 0, 0, 0}, {0, 0, -84 Sqrt[6], 0}, {0, 0, 504,
0}}, {{0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 360}}}}


I tried in many ways but I can't figure it out. Thanks in advance.

• You can try MapIndexed for level {-1} in combination with Reap/Sow. May 16, 2013 at 23:06

This basically excludes all the upper triangular part of each submatrix and all the submatrices in the upper triangular part of the tensor.

With[{lt = LowerTriangularize, ca = ConstantArray},
Pick[b, # * ca[#, {4, 4}] & @ lt @ ca[1, {4, 4}], 1]
] ~Flatten~ 3

• I think this is what Im going to use, because very soon I will have to do it for a 20x20x20x20 array and this seems pretty quick. May 17, 2013 at 10:56

Do these give what you want?

rv1[b_] :=
Module[{tag},
SparseArray[{x_, y_, z_, w_} /; x >= y && z >= w :>
b[[x, y, z, w]], Dimensions@b, tag]["NonzeroValues"]];

rv2[b_] :=
Extract[b,
Cases[Range@4~Tuples~{4}, {x_, y_, z_, w_} /; x >= y && z >= w]];

rv3[b_] :=
Extract[b,
Table[{x, y, z, w}, {y, 4}, {x, y, 4}, {w, 4}, {z, w, 4}]~Flatten~
3];

rv4[b_] :=
With[{extract =
With[{parts = Table[{x, y}, {y, 4}, {x, y, 4}]~Flatten~1},
Extract[#, parts] &]}, Flatten[extract /@ extract@b, 1]];

(*This is @Spawn1701D's suggestion*)
spawn[b_] :=
Reap[MapIndexed[Sow[#1, #2[[1]] >= #2[[2]] && #2[[3]] >= #2[[4]]] &,
b, {4}], True][[-1, 1]]

(* from @MrWizard *)
wiz1[b_] :=
With[{lt = LowerTriangularize, ca = ConstantArray},
Pick[b, #*ca[#, {4, 4}] &@lt@ca[1, {4, 4}], 1]]~Flatten~3


Let's check the results

contenders = {rv1, rv2, rv3, rv4, spawn, wiz1};

Equal @@ Sort /@ Through@contenders[b]
(* True *)


Some timings

Through@(Composition @@@
Tuples@{{First}, {AbsoluteTiming}, contenders})[b]

(* {0.019003, 0., 0.00100, 0., 0.002001, 0.} *)


Using this function TimingAverage that I took from somewhere sometime and never looked into

TimingAverage[expr_, time_ : 1.] :=
Module[{t = Timing[expr;][[1]], tries = 1},
While[t < time, tries *= 2;
t = Timing[Do[expr, {tries}];][[1]];];
t/tries]


You get

Through@(Composition @@@
Tuples@{{TimingAverage}, contenders})[b]

(* {0.0185251, 0.000716020, 0.000154249, 0.0000656987, 0.00237658, 0.000148536} *)


So, from top to bottom, rv4, wiz1, rv3, rv2, spawn, rv1

• Thanks. I'm curious, do you personally use TimingAverage over the timeAvg function (originally from Timo) I've posted many times willfully, or was this an ad hoc function? May 16, 2013 at 23:40
• @Mr.Wizard because I never even looked into any of them to see pros and cons and I have it packaged up from a long time ago. If you think the other one is better, edit it in and I'll trust you and change my packaging ;)
– Rojo
May 16, 2013 at 23:41
• No need; it does the same thing I believe, only I chose a multiplier of five. May 16, 2013 at 23:48
• Could you add timings for the method I just posted, please? May 16, 2013 at 23:50
• @Mr.Wizard already done. +1 btw
– Rojo
May 16, 2013 at 23:53

A flexible solution :

Reap[ReplacePart[
b,
({x_, y_, z_, w_} /; x >= y && z >= w) :> Sow[Extract[b, {x, y, z, w}]]
];][[2, 1]]


So far I know, ReplacePart[] is the only instruction that accepts a pattern of indexes as argument. That's the reason why I use ReplacePart[]. In fact we don't care of the result of the replacements.

This solution is very flexible but not efficient in terms of quickness and memory occupation.

Two Pick's; first picks out matrices, second elements:

With[{s = ConstantArray[True, {3, 3}] // UpperTriangularize},
Map[
Pick[#, s] &,
Pick[m, s], {-3}]] // Flatten


Pick also works with SparseArray objects, so:

Pick[b, SparseArray[{x_, y_, z_, w_} /; x >= y && z >= w -> True,
Dimensions@b]] // Flatten
`
• I like this. +1 May 21, 2013 at 13:39