Mainly for Raul, as nobody else has commented, at least publicly.
OK - I didn't mean to spend much time on this, but I started wondering
how a
coordinate representation would work.
The 3-d boolean representation is more intuitive, but one advantage of the
coordinate approach is that the representation is less ragged, at least
when
we're dealing with all n-polycubes ; I haven't yet dealt with sets of
polycubes
with different numbers of component cubes. So the set of 29
penta-polycubes
has shape 29 5 3, for example.
The main problem with this approach is that I didn't quite solve the
canonicity
task. "Rotating" (and sorting) the coords does most of what's needed, but
I've resorted to flipping to binary, applying Raul's "canon" routine,
then
flipping back to coordinates. Naturally this adds to the computation
time. It would
be much neater to remove all the isomers at the coordinate level.
So here's quite a long screed copied from my current script, which
includes
Raul's definitions. I'll keep those in for completeness but will remove
his inline
comments; these definitions appear first. Apologies if I've omitted
anything.
NB. https://en.wikipedia.org/wiki/Polycube
rotx=: |."2@(2 1&|:) NB. rotate in the yz plane about the x axis
roty=: |.@(1 0 2&|:) NB. rotate in the xz plane about the y axis
rotz=: |.@(2 1 0&|:) NB. rotate in the xy plane about the z axis
rotall=: {{
'X Y Z'=: y
rotx^:X roty^:Y rotz^:Z x
}}"3 1&(4 4 4#:(i.20),24+i.4)
extend=: {{
Y=. (-2+$y){.(1+$y){.y
new=. ($#:I.@,)Y < (_1&|. +. 1&|. +. 1&|."2 +. _1&|."2 +. 1&|."1 +.
_1&|."1) Y
~.new {{canon 1 (,:x)} y}}"1 3 Y
}}
offset=: {{ 1 i.~ +./+./"2 x |: *y }}
trim=: {{ (0 offset y)}."3 (1 offset y)}."2 (2 offset y)}."1 y }}
rtrim=: trim&.|.&.:(|."2)&.:(|."1)
canon=: {{ rtrim {. (/: +/@p:@I.@,"3) rotall trim y }}^:2
extendall=: {{ ~.;<@extend"3 y }}
NB. #@> extendall&.>^:(i.8)<1 1 1 1$1
NB. 1 1 2 8 29 166 1023 6922
NB. My stuff follows (the copy & paste has removed my indentations!)
NB. I've repaired a couple of line-wraps - others might lurk!
NB. =============================================================
NB. Prepare all 24 general 3x3 rotation matrices
mm =: +/ . * matrix multiply
I3 =: =i.3 3x3 identity matrix
NB. define some helper matrices
Rx =: >1 0 0; 0 0 _1; 0 1 0 NB. pi/2 rotation around x axis
Ry =: _1 |. _1 |."1 Rx NB. y
Rz =: 1 |. 1 |."1 Rx NB. z
Rx3=: Rx mm Rx2=: mm~ Rx
Ry3=: Rx mm Ry2=: mm~ Ry
Rz3=: Rx mm Rz2=: mm~ Rz
Rxs=: I3,Rx,Rx2,:Rx3 NB. rotations by 0, pi/2, pi, 3pi/2 around x axis
Rys=: I3,Ry,Ry2,:Ry3 NB. rotations by 0, pi/2, pi, 3pi/2 around y axis
Rzs=: I3,Rz,Rz2,:Rz3 NB. rotations by 0, pi/2, pi, 3pi/2 around z axis
NB. combine them to generate all 24 rotation matrices
irots =: (4 4 4#:(i.20),24+i.4)
makeRots =: {{
r =. {{1 (0 1 2 ,each y) } 3 3$0}}"1 tap 3 NB. all 6 permutations of 3
axes' labels ==> 3 x 1s in 3x3 matrices
s =. ~./:~,/1 _1,"0 1/,/(,"0)/~1 _1 NB. 8 sets of triples of signs
rs =. ,/s *"1/ r NB. output 48 (!) matrices
}}
NB. Can't see yet how to limit them to 24 except by comparison with
Raul's rotall results on
NB. an irregular set of 3 points
# q=. {{/:~ {{y +"1 <./y}} y -"1 <./ y}}each ($#:I.@:,)each <"3 rotall 1
(0 0 0;1 1 1;1 2 3) } 2 3 4$0 24
#r =. {{ /:~ {{y +"1 <./y}} y -"1 <./ y}}each<"2(0 0 0,1 1 1,:1 2 3) mm
"2 Rots =: makeRots''48
+/ok =. r e. q
24
NB. Global set of rotation matrices to help with finding canonical
coordinate representations
Rots =: (q i.~ ok#r){ ok# Rots NB. !!!!!!!
NB. Demonstration that the 24 Rots can be used for rotation analogously
to Raul's rotall function
NB. compare with Raul's results on all flips of these 3 points - not
face-adjacent cubes as it happens!
NB. but that doesn't matter much
NB. limited to 1st 16 for easier display...
NB. 16{.{{/:~ {{y +"1 <./y}} y -"1 <./ y}}each ($#:I.@:,)each <"3 rotall
1 (0 0 0;1 1 1;1 2 3) } 2 3 4$0
NB.
┌─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┐
NB. │0 0 0│0 2 1│0 1 2│0 0 1│0 1 3│0 0 1│0 0 0│0 3 0│0 0 3│0 2 0│0 2 3│0
0 0│0 1 0│0 0 0│0 0 3│0 3 1│
NB. │1 1 1│2 1 1│0 2 0│1 1 0│1 1 1│1 2 1│1 0 2│1 1 0│0 1 1│1 1 1│1 0 0│2
1 0│1 0 1│1 1 1│1 1 2│1 2 0│
NB. │1 2 3│3 0 0│1 0 3│3 2 0│2 0 0│2 3 0│2 1 3│2 0 1│1 2 0│3 0 1│1 1 2│3
2 1│2 0 3│2 3 1│2 1 0│2 0 0│
NB.
└─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┘
NB. 16{.{{ /:~ {{y +"1 <./y}} y -"1 <./ y}}each<"2(0 0 0,1 1 1,:1 2 3)
mm "2 Rots
NB.
┌─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┐
NB. │0 0 0│0 2 1│0 1 2│0 0 1│0 1 3│0 0 1│0 0 0│0 3 0│0 0 3│0 2 0│0 2 3│0
0 0│0 1 0│0 0 0│0 0 3│0 3 1│
NB. │1 1 1│2 1 1│0 2 0│1 1 0│1 1 1│1 2 1│1 0 2│1 1 0│0 1 1│1 1 1│1 0 0│2
1 0│1 0 1│1 1 1│1 1 2│1 2 0│
NB. │1 2 3│3 0 0│1 0 3│3 2 0│2 0 0│2 3 0│2 1 3│2 0 1│1 2 0│3 0 1│1 1 2│3
2 1│2 0 3│2 3 1│2 1 0│2 0 0│
NB.
└─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┘
NB. Like trim - easier -
NB. eg to cut planes x < 2 , just subtract 2 0 0
tc =: trimcoords =: {{ /:~ y -"1 <./ y}}
NB. convert binary array representation to coord rep
b2c =: ($#:I.@:,) :. c2b
NB. convert coord rep to binary array
c2b =: {{
1 (<"1 y) } 0$~ >: >./ y
}} :. b2c
NB. global for add3 (function def follows)
ADD3 =: _1 0 0,1 0 0,0 _1 0, 0 1 0, 0 0 _1,: 0 0 1
NB. append new cubes by means of their 3-coords
NB. too loopy but ok for now...
add3 =: {{
out =. 0$~ 0,1 0+$y
head=. 0$~0, {:$y
tail=. y
NB. left as loop for now!!!
for_i. i.#y do.
iy =. i{y
new =. iy +"1 ADD3
if. +/ ok =. -. new e. y do.
for_newi. ok#new do.
out =. out, /:~ head, newi, tail
end.
end.
head =. head, iy
tail =. }. tail
end.
NB. "trim" coords to help with checking canonicity
tc"2 ~. out
}}
mycanon =: {{ ~. mycanon1 "2 ~. add3 y }}
mycanon1 =: {{
if. 1 < # q =. Rots ~.@:(tc@:mm"2~)y do.
q{~{./:@:((+/@:p:@:I.@:,)"2) q
else.
{.q
end.
}}
myexpand =: {{ r =. ~. /:~,/ mycanon"2 y
if. 0 = +/, {.r do. r =. }. r end.}}
NB. Performance - starts ok!
NB. Build Raul's extended sets for comparison purposes, and as seeds
NB. uncomment these to run!
NB. e7=.extendall each^:6<1 1 1 1$1
NB. e6=.extendall each^:5<1 1 1 1$1
NB. e5=.extendall each^:4<1 1 1 1$1
NB. e4=.extendall each^:3<1 1 1 1$1
NB. e3=.extendall each^:2<1 1 1 1$1
NB. e2=.extendall each^:1<1 1 1 1$1
NB. Compare results from "myexpand", run on Raul's arrays turned into
coordinate lists, with his own results
NB. timer'$ myexpand b2c"3;e2' NB. b2c"3 turns binary array ;e2 into
sets of coords
NB. ┌─────────┬─────┐
NB. │0.0019989│2 3 3│
NB. └─────────┴─────┘
NB. timer'#extendall >e2'
NB. ┌─────────┬─┐
NB. │0.0059967│2│
NB. └─────────┴─┘
NB. timer'$ myexpand b2c"3;e3'
NB. ┌──────────┬─────┐
NB. │0.00800323│8 4 3│
NB. └──────────┴─────┘
NB. timer'#extendall >e3'
NB. ┌─────────┬─┐
NB. │0.0130005│8│
NB. └─────────┴─┘
NB. timer'$ myexpand b2c"3;e4'
NB. ┌─────────┬──────┐
NB. │0.0279999│29 5 3│
NB. └─────────┴──────┘
NB. timer'#extendall >e4'
NB. ┌─────────┬──┐
NB. │0.0589981│29│
NB. └─────────┴──┘
NB. timer'$q6=: myexpand b2c"3;e5'
NB. ┌────────┬───────┐
NB. │0.129005│169 6 3│ NB. too big - whoops!
NB. └────────┴───────┘
NB. Raul's canon can clean it up!!!
NB. timer'$bq6=:~.>canon@:c2b each <"2 q6 '
NB. ┌─────────┬─────────┐
NB. │0.0830002│166 3 4 6│
NB. └─────────┴─────────┘
NB. timer'$cq6 =: b2c"3 bq6' NB. restoring to coord form
NB. ┌───────────┬───────┐
NB. │0.000999451│166 6 3│
NB. └───────────┴───────┘
NB. timer'#extendall >e5'
NB. ┌────────┬───┐
NB. │0.276001│166│
NB. └────────┴───┘
NB. Wrap up the clean-up ...
NB. timer'${{~.canon&.:c2b "2 myexpand b2c"3 y}};e5 '
NB. ┌────────┬───────┐
NB. │0.230003│166 6 3│
NB. └────────┴───────┘
NB. So incorporate this clean-up in a new routine:
myexpanda =: {{~.canon&.:c2b "2 myexpand b2c"3 y}}
NB. timer'$myexpanda ;e5'
NB. ┌────────┬───────┐
NB. │0.142998│166 6 3│
NB. └────────┴───────┘
NB. timer'$myexpanda ;e6'
NB. ┌────────┬────────┐
NB. │0.976997│1023 7 3│
NB. └────────┴────────┘
NB. timer'#extendall >e6'
NB. ┌─────┬────┐
NB. │2.119│1023│
NB. └─────┴────┘
NB. timer'$myexpanda ;e7'
NB. ┌─────┬────────┐
NB. │6.975│6922 8 3│
NB. └─────┴────────┘
NB. timer'#extendall >e7'
NB. ┌──────┬────┐
NB. │16.231│6922│
NB. └──────┴────┘
Sorry for the length of this offering,
Mike
On 15/07/2023 22:03, Raul Miller wrote:
https://en.wikipedia.org/wiki/Polycube
Today, I've been playing with polycubes (which are polyhedrons formed
from adjacent cubes).
This requires both a representation -- a rank 3 bit array works
nicely, and a mechanism for working with cubic symmetry. For that, I
went with:
rotx=: |."2@(2 1&|:) NB. rotate in the yz plane about the x axis roty=: |.@(1 0
2&|:) NB. rotate in the xz plane about the y axis rotz=: |.@(2 1 0&|:)
NB. rotate in the xy plane about the z axis rotall=: {{ 'X Y Z'=. y
rotx^:X roty^:Y rotz^:Z x }}"3 1&(4 4 4#:(i.20),24+i.4)
Here, rotall gives all 24 rotations of cubic symmetry for a given
collection of 3d cubes represented by a bit array.
I decided to use a brute force mechanism to generate polycubes. In
other words, pad the array representation and add a single cube in
each of the locations adjacent to an existing cube:
extend=: {{
Y=. (-2+$y){.(1+$y){.y
new=. ($#:I.@,)Y < (_1&|. +. 1&|. +. 1&|."2 +. _1&|."2 +. 1&|."1 +.
_1&|."1) Y
~.new {{canon 1 (,:x)} y}}"1 3 Y
}}
Here, Y is the y argument with an extra 0 of padding in all
directions. And 'new' is the list of indices of empty locations
adjacent to cube locations.
I also need a verb (canon) to reduce these generated polycubes to
canonical form. A lot of this is removing unnecessary rows, columns
and planes of zeros. The rest of it is picking arbitrarily (but
consistently) from the 24 possible symmetric rotations:
offset=: {{ 1 i.~ +./+./"2 x |: *y }}
trim=: {{ (0 offset y)}."3 (1 offset y)}."2 (2 offset y)}."1 y }}
rtrim=: trim&.|.&.:(|."2)&.:(|."1)
canon=: {{ rtrim {. (/: +/@p:@I.@,"3) rotall trim y }}^:2
Finally, I want to work with all polycubes of a specific order, and to
get larger batch I can use;
extendall=: {{ ~.;<@extend"3 y }}
And, testing, this works:
#@> extendall&.>^:(i.8)<1 1 1 1$1
1 1 2 8 29 166 1023 6922
But, this is slow -- since I'm brute forcing all possibilities then
discarding duplicates, I need almost seven seconds on the little
laptop I'm currently using, to generate the above sequence.
(Also, for visualization of individual polycubes,
https://code.jsoftware.com/wiki/Scripts/Plot_3D works passibly well.)
Does anyone here have the geometric insight to improve on this approach?
Thanks,
----------------------------------------------------------------------
For information about J forums see http://www.jsoftware.com/forums.htm
