Announcement

Collapse
No announcement yet.
X
  • Filter
  • Time
  • Show
Clear All
new posts

  • Can Mata find the moves to solve a small cube, and can Stata show them?

    While experimenting with coordinates and 3D graphs, I wondered whether Stata/Mata could do more than draw a cube. Could it find a solution to a small 2x2 cube problem and then show the moves?

    This version starts with 24 sticker colors. It knows the legal turns and the target arrangement, but receives neither the scramble history nor a solution sequence.

    Mata explores states up to three moves from each end, removes duplicates, and finds the meeting with the shortest combined path. Stata then draws that verified path, using two fixed views so that all six faces can be seen. The panels show the chosen solution, not every branch explored by the search.

    For the supplied state, Stata/MP 16.1 reports:

    Code:
    FOUND BY SEARCH: F2 U R2 F' U' R'
    Length: 6 face turns; states: 385 start + 385 goal
    R, U and F mean the right, upper and front faces. A prime means counterclockwise when viewed from outside that face; 2 means a half turn. Each counts as one move.
    Click image for larger version

Name:	cube_search_2x2.png
Views:	1
Size:	758.4 KB
ID:	1787793


    The scope is deliberately small: only these three faces are turned, the target orientation is fixed, and the search covers solutions of at most six moves. Within those rules the returned path is shortest. Failure to find one would mean "not found within the limit", not "this cube is impossible". The final arrangement is checked against all 24 target stickers.

    I found the distinction between turning a face and changing the camera particularly useful: one changes the puzzle state, while the other changes how I see it.

    The complete code follows, including the graphics. Paste the whole block into one do-file. No downloads, packages or data files are needed. It preserves the current dataset, uses sb2_/SB2_ names, and writes cube_search_2x2.png in the current folder (replacing that file if it exists). The graphics are static snapshots after each turn.

    Code:
    version 16.1
    set more off
    preserve
    capture quietly mata: sb2_goal()
    if _rc {
    mata:
    struct sb2_tree {
        real matrix states
        real colvector depth, parent, move
        transmorphic scalar seen
    }
    struct sb2_result {
        real scalar found, length, startnodes, goalnodes, checks
        real rowvector path
        real matrix layers
    }
    
    real matrix sb2_geometry()
    {
        real matrix N,G
        real rowvector n,a,b,c
        real scalar f,i,j
        N=(1,0,0\-1,0,0\0,0,1\0,0,-1\0,-1,0\0,1,0)
        G=J(0,7,.)
        for(f=1;f<=6;f++) {
            n=N[f,]
            if(abs(n[1])) { a=(0,1,0); b=(0,0,1); }
            else if(abs(n[2])) { a=(1,0,0); b=(0,0,1); }
            else { a=(1,0,0); b=(0,1,0); }
            for(i=-1;i<=1;i=i+2) for(j=-1;j<=1;j=j+2) {
                c=n+.5*i*a+.5*j*b
                G=G\(c,n,f)
           }
       }
        return(G)
    }
    
    real matrix sb2_moves(real matrix G)
    {
        real matrix P,Q,H
        real rowvector axis,ix,base
        real scalar f,i,k
        P=J(9,24,.)
        // Row vectors multiply Q. Clockwise seen from outside R, U, F.
        for(f=1;f<=3;f++) {
            if(f==1) { Q=(1,0,0\0,0,-1\0,1,0); axis=(1,0,0); }
            if(f==2) { Q=(0,-1,0\1,0,0\0,0,1); axis=(0,0,1); }
            if(f==3) { Q=(0,0,-1\0,1,0\1,0,0); axis=(0,-1,0); }
            H=G
            for(i=1;i<=24;i++) if(G[i,1..3]*axis'>0) {
                H[i,1..3]=G[i,1..3]*Q
                H[i,4..6]=G[i,4..6]*Q
           }
            k=3*f-2
            for(i=1;i<=24;i++) {
                ix=selectindex(rowsum(abs(G[,1..6]:-H[i,1..6])):==0)'
                assert(cols(ix)==1)
                P[k,ix[1]]=i
           }
            base=P[k,]
            P[k+1,]=base[base]
            P[k+2,]=base[base[base]]
            assert(all(sort(base',1):==(1::24)))
            assert(all(base[base[base[base]]]:==(1..24)))
       }
        return(P)
    }
    
    real rowvector sb2_goal()
    {
        return((1,1,1,1,2,2,2,2,3,3,3,3,4,4,4,4,5,5,5,5,6,6,6,6))
    }
    string scalar sb2_key(real rowvector s)
    {
        return(invtokens(strofreal(s,"%1.0f"),""))
    }
    real scalar sb2_inverse(real scalar m)
    {
        return(3*floor((m-1)/3)+4-mod(m-1,3)-1)
    }
    string scalar sb2_names(real rowvector path)
    {
        string rowvector names
        names=("R","R2","R'","U","U2","U'","F","F2","F'")
        if(cols(path)==0) return("(none)")
        return(invtokens(names[path]))
    }
    real rowvector sb2_parse(string scalar text)
    {
        string rowvector words,names
        real rowvector path,ix
        real scalar i
        words=tokens(text)
        names=("R","R2","R'","U","U2","U'","F","F2","F'")
        path=J(1,0,.)
        for(i=1;i<=cols(words);i++) {
            ix=selectindex(names:==words[i])
            assert(cols(ix)==1)
            path=path,ix[1]
       }
        return(path)
    }
    real rowvector sb2_apply(real rowvector s,real rowvector path,real matrix P)
    {
        real scalar i
        real rowvector out
        out=s
        for(i=1;i<=cols(path);i++) out=out[P[path[i],]]
        return(out)
    }
    
    struct sb2_tree scalar sb2_bfs(real rowvector root,real matrix P,real scalar radius)
    {
        struct sb2_tree scalar T
        real scalar cap,head,n,d,m
        real rowvector s
        string scalar key
        // Upper bound on all 9-branch paths. This demo limits radius to 4.
        cap=(9^(radius+1)-1)/8
        T.states=J(cap,24,.)
        T.depth=J(cap,1,.)
        T.parent=J(cap,1,0)
        T.move=J(cap,1,0)
        T.seen=asarray_create("string",1)
        T.states[1,]=root; T.depth[1]=0
        asarray(T.seen,sb2_key(root),1)
        head=1; n=1
        while(head<=n) {
            d=T.depth[head]
            if(d<radius) for(m=1;m<=9;m++) {
                s=T.states[head,P[m,]]
                key=sb2_key(s)
                if(!asarray_contains(T.seen,key)) {
                    n++
                    T.states[n,]=s; T.depth[n]=d+1
                    T.parent[n]=head; T.move[n]=m
                    asarray(T.seen,key,n)
               }
           }
            head++
       }
        T.states=T.states[1..n,]
        T.depth=T.depth[1..n]; T.parent=T.parent[1..n]; T.move=T.move[1..n]
        return(T)
    }
    real rowvector sb2_trace(struct sb2_tree scalar T, real scalar i)
    {
        real rowvector path
        path=J(1,0,.)
        while(T.parent[i]>0) {
            path=T.move[i],path
            i=T.parent[i]
       }
        return(path)
    }
    
    struct sb2_result scalar sb2_solve(real rowvector start,real matrix P,real scalar radius)
    {
        struct sb2_tree scalar A,B
        struct sb2_result scalar R
        real scalar i,j,d,best,ai,bi,c
        real rowvector back
        string scalar key
        assert(cols(start)==24 & rows(start)==1)
        assert(all(start:==floor(start)) & min(start)>=1 & max(start)<=6)
        for(c=1;c<=6;c++) assert(sum(start:==c)==4)
        assert(radius>=0 & radius<=4 & radius==floor(radius))
        // Only start, legal moves and depth limit enter. No scramble history.
        A=sb2_bfs(start,P,radius)
        B=sb2_bfs(sb2_goal(),P,radius)
        R.startnodes=rows(A.states); R.goalnodes=rows(B.states)
        R.checks=rows(A.states); R.found=0; R.length=.
        R.path=J(1,0,.); R.layers=J(radius+1,3,0)
        for(d=0;d<=radius;d++) R.layers[d+1,]=(d,sum(A.depth:==d),sum(B.depth:==d))
        best=.; ai=0; bi=0
        // Check ALL meetings: first meeting need not give the shortest combined path.
        for(i=1;i<=rows(A.states);i++) {
            key=sb2_key(A.states[i,])
            if(asarray_contains(B.seen,key)) {
                j=asarray(B.seen,key)
                d=A.depth[i]+B.depth[j]
                if(d<best) { best=d; ai=i; bi=j; }
           }
       }
        if(ai==0) return(R)
        R.path=sb2_trace(A,ai)
        back=sb2_trace(B,bi)
        for(i=cols(back);i>=1;i--) R.path=R.path,sb2_inverse(back[i])
        R.found=1; R.length=cols(R.path)
        assert(R.length==best)
        assert(all(sb2_apply(start,R.path,P):==sb2_goal()))
        return(R)
    }
    
    // Static orthographic rendering. At quarter-turn endpoints the cube is convex;
    // outward-facing stickers are sufficient, without a general z-buffer renderer.
    void sb2_render(real rowvector state,real matrix G)
    {
        real matrix out,pts
        real rowvector right,up,cam,n,a,b,c
        real scalar view,i,poly,shift
        out=J(0,5,.); poly=0
        for(view=1;view<=2;view++) {
            right=(.8,.6,0); up=(-.36,.48,.8); cam=(.48,-.64,.6)
            if(view==2) { right=-right; cam=-cam; }
            shift=(view==1 ? 0 : 3.8)
            for(i=1;i<=24;i++) {
                n=G[i,4..6]
                if(n*cam'<=0) continue
                if(abs(n[1])) { a=(0,1,0); b=(0,0,1); }
                else if(abs(n[2])) { a=(1,0,0); b=(0,0,1); }
                else { a=(1,0,0); b=(0,1,0); }
                // First sticker of each face also draws its dark backing square.
                if(mod(i-1,4)==0) {
                    pts=(n-a-b\n+a-b\n+a+b\n-a+b\n-a-b)
                    poly++
                    out=out\(J(5,1,poly),pts*right':+shift,pts*up',J(5,1,0),J(5,1,view))
               }
                c=G[i,1..3]
                pts=(c-.46*a-.46*b\c+.46*a-.46*b\c+.46*a+.46*b\c-.46*a+.46*b\c-.46*a-.46*b)
                poly++
                out=out\(J(5,1,poly),pts*right':+shift,pts*up',J(5,1,state[i]),J(5,1,view))
           }
       }
        st_addobs(rows(out))
        st_addvar("double",("polygon","px","py","color","view"))
        st_store(.,.,out)
    }
    end
    }
    clear
    mata:
    sb2_G=sb2_geometry()
    sb2_P=sb2_moves(sb2_G)
    // Current colors only: no scramble history or solution is supplied.
    // Slots: R(1-4), L(5-8), U(9-12), D(13-16), F(17-20), B(21-24).
    // Colors: 1 red, 2 orange, 3 white, 4 yellow, 5 green, 6 blue.
    sb2_start=(2,6,5,3,5,6,2,1,3,3,1,6,2,4,3,4,4,2,5,4,6,5,1,1)
    sb2_answer=sb2_solve(sb2_start,sb2_P,3)
    if(!sb2_answer.found) {
        printf("No solution within 6 face turns; this does not prove impossibility.\n")
        _error(3498)
    }
    printf("INPUT STATE: %s\n",sb2_key(sb2_start))
    printf("FOUND BY SEARCH: %s\n",sb2_names(sb2_answer.path))
    printf("Length: %g face turns; states: %g start + %g goal\n",sb2_answer.length,sb2_answer.startnodes,sb2_answer.goalnodes)
    st_numscalar("SB2_length",sb2_answer.length)
    st_numscalar("SB2_startnodes",sb2_answer.startnodes)
    st_numscalar("SB2_goalnodes",sb2_answer.goalnodes)
    st_matrix("SB2_layers",sb2_answer.layers)
    sb2_states=sb2_start
    for(sb2_i=1;sb2_i<=sb2_answer.length;sb2_i++) sb2_states=sb2_states\sb2_apply(sb2_states[rows(sb2_states),],sb2_answer.path[sb2_i],sb2_P)
    assert(all(sb2_states[rows(sb2_states),]:==sb2_goal()))
    end
    local panels
    forvalues step=0/`=scalar(SB2_length)' {
    clear
    mata: sb2_render(sb2_states[`step'+1,],sb2_G)
    local applied "Start"
    if `step'>0 {
    mata: st_local("applied",sb2_names(sb2_answer.path[`step']))
    }
    local next "SOLVED: all 24 facelets checked"
    if `step'<scalar(SB2_length) {
    mata: st_local("next","Next: "+sb2_names(sb2_answer.path[`step'+1]))
    }
    local colors `""40 45 53" "204 57 57" "239 145 50" "244 244 235" "246 211 64" "49 151 113" "66 111 185""'
    quietly levelsof polygon, local(polys)
    local plots
    foreach p of local polys {
    quietly summarize color if polygon==`p', meanonly
    local ci=r(min)+1
    local rgb : word `ci' of `colors'
    local plots `"`plots' (area py px if polygon==`p', color("`rgb'") lcolor("40 45 53") lwidth(vthin))"'
    }
    twoway `plots', aspectratio(.5070422535) xscale(range(-1.65 5.45) noextend off) yscale(range(-1.6 2.0) noextend off) xlabel(none) ylabel(none) legend(off) text(1.75 0 "Front / right / up", size(small) color(gs6)) text(1.75 3.8 "Back / left / down", size(small) color(gs6)) title("Step `step': `applied'", size(medium)) subtitle("`next'", size(small)) graphregion(color(white)) plotregion(color(white)) xsize(9) ysize(5.8) name(sb2_step`step', replace)
    local panels `panels' sb2_step`step'
    }
    twoway scatteri 0 0, msymbol(none) xscale(range(0 10) off) yscale(range(0 10) off) xlabel(none) ylabel(none) legend(off) text(8.7 5 "Input: 24 colors, no move history", size(medium)) text(6.9 5 "Start side: `=scalar(SB2_startnodes)' distinct states", size(medium) color("40 104 150")) text(5.8 5 "Goal side: `=scalar(SB2_goalnodes)' distinct states", size(medium) color("40 104 150")) text(4.0 5 "Depth:       0     1     2      3", size(medium)) text(3.0 5 "Start: `=SB2_layers[1,2]' / `=SB2_layers[2,2]' / `=SB2_layers[3,2]' / `=SB2_layers[4,2]'", size(medium)) text(1.1 5 "Shortest path found: `=scalar(SB2_length)' face turns", size(medium)) title("The search behind the pictures", size(medium)) graphregion(color(white)) plotregion(color(white)) xsize(9) ysize(5.8) name(sb2_summary, replace)
    graph combine `panels' sb2_summary, cols(2) imargin(tiny) graphregion(color(white)) title("A 2x2 cube solved by search in Mata", size(medium)) note("Mata finds the moves; Stata draws the verified states. Both fixed views together show all six faces." "R/right, U/up, F/front. Prime: counterclockwise viewed from outside. 2: a half turn.", size(vsmall)) xsize(12) ysize(16) name(sb2_overview, replace)
    graph export "cube_search_2x2.png", width(1800) replace
    restore
    display "CUBE_SANDBOX_PASS"
    I would be interested in a simpler way to show the search itself, or a clearer visual cue for the next turn. This seems a useful small playground for coordinates, search and graphics together.
    One question this leaves me with: does a cube with more stickers in their target positions necessarily need fewer moves to solve? I would like to compare “looking closer” with the shortest-path distance.
    Last edited by Choonjoo Lee; 27 Oct 2026, 07:34.

  • #2
    More correct colors, fewer moves? A small cube experiment in Mata

    After letting Mata find a solution to a 2x2 cube, I wondered: if more stickers already match the target, should fewer moves remain?

    Here, "correct" means the sticker color matches that position in a fixed target, not merely that several stickers share a color. The picture shows the target and both sides of each example.

    Mata starts from the target and searches all states within three moves: 1 + 9 + 54 + 321 = 385 distinct states. Breadth-first search gives their exact minimum distances under the specified moves.

    Code:
                  Correct colors     Minimum moves
    State A           16/24                3
    State B           10/24                2
    These are deliberately selected examples: A maximizes correct colors among states three moves away; B minimizes them among states two moves away. Ties use search order. They are not a random sample.
    Click image for larger version

Name:	cube_closeness.png
Views:	1
Size:	528.1 KB
ID:	1787797

    One shortest solution for A is R2 U2 R2. Its correct-color count goes 16 -> 12 -> 16 -> 24. For B, U' R' gives 10 -> 16 -> 24. More correct colors did not mean fewer moves; even a shortest route can initially make that score worse.

    R, U and F turn the right, upper and front faces. A prime reverses the turn; 2 means a half turn. Each counts as one move. The target orientation is fixed, and only these three faces are turned.

    This leaves me with another question for the next experiment: if I rotate the view, the cube may look different, but should its distance from the target change? Suggestions for other small comparisons are welcome.

    The complete code below runs in Stata 16.1. Paste the whole block into one do-file. It needs no packages or input files, preserves the current dataset, and uses cq_/CQ_ names. It writes cube_closeness.png/.gph, closeness_results.txt, closeness_states.csv and closeness_counts.csv in the current folder, replacing those files if present. All drawing and searching use Stata/Mata.

    Code:
    version 16.1
    set more off
    preserve
    capture quietly mata: cq_goal()
    if _rc {
    mata:
    struct cq_tree {
        real matrix states
        real colvector depth, parent, move
        transmorphic scalar seen
    }
    real matrix cq_geometry()
    {
        real matrix N,G
        real rowvector n,a,b,c
        real scalar f,i,j
        N=(1,0,0\-1,0,0\0,0,1\0,0,-1\0,-1,0\0,1,0)
        G=J(0,7,.)
        for(f=1;f<=6;f++) {
            n=N[f,]
            if(abs(n[1])) { a=(0,1,0); b=(0,0,1); }
            else if(abs(n[2])) { a=(1,0,0); b=(0,0,1); }
            else { a=(1,0,0); b=(0,1,0); }
            for(i=-1;i<=1;i=i+2) for(j=-1;j<=1;j=j+2) {
                c=n+.5*i*a+.5*j*b
                G=G\(c,n,f)
           }
       }
        return(G)
    }
    
    real matrix cq_moves(real matrix G)
    {
        real matrix P,Q,H
        real rowvector axis,ix,base
        real scalar f,i,j,t,k
        P=J(9,24,.)
        // Row vectors multiply Q. Clockwise seen from outside R, U, F.
        for(f=1;f<=3;f++) {
            if(f==1) { Q=(1,0,0\0,0,-1\0,1,0); axis=(1,0,0); }
            if(f==2) { Q=(0,-1,0\1,0,0\0,0,1); axis=(0,0,1); }
            if(f==3) { Q=(0,0,-1\0,1,0\1,0,0); axis=(0,-1,0); }
            H=G
            for(i=1;i<=24;i++) if(G[i,1..3]*axis'>0) {
                H[i,1..3]=G[i,1..3]*Q
                H[i,4..6]=G[i,4..6]*Q
           }
            k=3*f-2
            for(i=1;i<=24;i++) {
                ix=selectindex(rowsum(abs(G[,1..6]:-H[i,1..6])):==0)'
                assert(cols(ix)==1)
                P[k,ix[1]]=i
           }
            base=P[k,]
            P[k+1,]=base[base]
            P[k+2,]=base[base[base]]
            assert(all(sort(base',1):==(1::24)))
            assert(all(base[base[base[base]]]:==(1..24)))
       }
        return(P)
    }
    
    real rowvector cq_goal()
    {
        return((1,1,1,1,2,2,2,2,3,3,3,3,4,4,4,4,5,5,5,5,6,6,6,6))
    }
    string scalar cq_key(real rowvector s)
    {
        return(invtokens(strofreal(s,"%1.0f"),""))
    }
    real scalar cq_inverse(real scalar m)
    {
        return(3*floor((m-1)/3)+4-mod(m-1,3)-1)
    }
    string scalar cq_names(real rowvector path)
    {
        string rowvector names
        names=("R","R2","R'","U","U2","U'","F","F2","F'")
        if(cols(path)==0) return("(none)")
        return(invtokens(names[path]))
    }
    real rowvector cq_apply(real rowvector s,real rowvector path,real matrix P)
    {
        real scalar i
        real rowvector out
        out=s
        for(i=1;i<=cols(path);i++) out=out[P[path[i],]]
        return(out)
    }
    
    struct cq_tree scalar cq_bfs(real rowvector root,real matrix P,real scalar radius)
    {
        struct cq_tree scalar T
        real scalar cap,head,n,d,m
        real rowvector s
        string scalar key
        // Upper bound on all 9-branch paths. This demo limits radius to 4.
        cap=(9^(radius+1)-1)/8
        T.states=J(cap,24,.)
        T.depth=J(cap,1,.)
        T.parent=J(cap,1,0)
        T.move=J(cap,1,0)
        T.seen=asarray_create("string",1)
        T.states[1,]=root; T.depth[1]=0
        asarray(T.seen,cq_key(root),1)
        head=1; n=1
        while(head<=n) {
            d=T.depth[head]
            if(d<radius) for(m=1;m<=9;m++) {
                s=T.states[head,P[m,]]
                key=cq_key(s)
                if(!asarray_contains(T.seen,key)) {
                    n++
                    T.states[n,]=s; T.depth[n]=d+1
                    T.parent[n]=head; T.move[n]=m
                    asarray(T.seen,key,n)
               }
           }
            head++
       }
        T.states=T.states[1..n,]
        T.depth=T.depth[1..n]; T.parent=T.parent[1..n]; T.move=T.move[1..n]
        return(T)
    }
    real rowvector cq_trace(struct cq_tree scalar T, real scalar i)
    {
        real rowvector path
        real scalar node
        node=i
        path=J(1,0,.)
        while(T.parent[node]>0) {
            path=T.move[node],path
            node=T.parent[node]
       }
        return(path)
    }
    
    // Static orthographic rendering. At quarter-turn endpoints the cube is convex;
    // outward-facing stickers are sufficient, without a general z-buffer renderer.
    void cq_render(real rowvector state,real matrix G)
    {
        real matrix out,pts
        real rowvector right,up,cam,n,a,b,c,px,py
        real scalar view,f,i,j,poly,s,shift
        out=J(0,5,.); poly=0
        for(view=1;view<=2;view++) {
            right=(.8,.6,0); up=(-.36,.48,.8); cam=(.48,-.64,.6)
            if(view==2) { right=-right; cam=-cam; }
            shift=(view==1 ? 0 : 3.8)
            for(i=1;i<=24;i++) {
                n=G[i,4..6]
                if(n*cam'<=0) continue
                if(abs(n[1])) { a=(0,1,0); b=(0,0,1); }
                else if(abs(n[2])) { a=(1,0,0); b=(0,0,1); }
                else { a=(1,0,0); b=(0,1,0); }
                // First sticker of each face also draws its dark backing square.
                if(mod(i-1,4)==0) {
                    pts=(n-a-b\n+a-b\n+a+b\n-a+b\n-a-b)
                    poly++
                    out=out\(J(5,1,poly),pts*right':+shift,pts*up',J(5,1,0),J(5,1,view))
               }
                c=G[i,1..3]
                pts=(c-.46*a-.46*b\c+.46*a-.46*b\c+.46*a+.46*b\c-.46*a+.46*b\c-.46*a-.46*b)
                poly++
                out=out\(J(5,1,poly),pts*right':+shift,pts*up',J(5,1,state[i]),J(5,1,view))
           }
       }
        st_addobs(rows(out))
        st_addvar("double",("polygon","px","py","color","view"))
        st_store(.,.,out)
    }
    void cq_experiment()
    {
        struct cq_tree scalar T
        real matrix P, chosen, census, route
        real colvector matches, ix
        real rowvector back, solution, s, neighbors
        real scalar a,b,j,k,d,fh
        string scalar label
        P=cq_moves(cq_geometry())
        T=cq_bfs(cq_goal(),P,3)
        matches=rowsum(T.states:==cq_goal())
        // Select extremes by a declared rule, not by visual preference.
        ix=selectindex(T.depth:==3)
        a=ix[selectindex(matches[ix]:==max(matches[ix]))[1]]
        ix=selectindex(T.depth:==2)
        b=ix[selectindex(matches[ix]:==min(matches[ix]))[1]]
        assert(matches[a]>matches[b] & T.depth[a]>T.depth[b])
        chosen=cq_goal()\T.states[a,]\T.states[b,]
        st_matrix("CQ_chosen",chosen)
        census=(T.depth,matches,T.states)
        st_matrix("CQ_census",census)
        printf("All %g distinct states within three moves:\n",rows(T.states))
        for(d=0;d<=3;d++) {
            ix=selectindex(T.depth:==d)
            printf("Depth %g: %g states; correct colors %g..%g\n",d,length(ix),min(matches[ix]),max(matches[ix]))
        }
        fh=fopen("closeness_results.txt","w")
        fput(fh,"Fixed target; moves R R2 R' U U2 U' F F2 F'; each costs one.")
        fput(fh,"Complete BFS ball of radius 3: "+strofreal(rows(T.states))+" distinct states.")
        for(j=1;j<=2;j++) {
            k=(j==1 ? a : b)
            label=(j==1 ? "A" : "B")
            back=cq_trace(T,k)
            solution=J(1,0,.)
            for(d=cols(back);d>=1;d--) solution=solution,cq_inverse(back[d])
            s=T.states[k,]
            route=(0,sum(s:==cq_goal()))
            for(d=1;d<=cols(solution);d++) {
                s=cq_apply(s,solution[d],P)
                route=route\(d,sum(s:==cq_goal()))
            }
            assert(all(s:==cq_goal()) & cols(solution)==T.depth[k])
            neighbors=J(1,9,.)
            for(d=1;d<=9;d++) neighbors[d]=sum(cq_apply(T.states[k,],d,P):==cq_goal())
            printf("%s: %g/24 correct; minimum %g moves; solve: %s\n",label,matches[k],T.depth[k],cq_names(solution))
            printf("Correct-color counts along this path: %s\n",invtokens(strofreal(route[,2]',"%2.0f")," -> "))
            printf("All nine next moves: correct-color counts %s\n",invtokens(strofreal(neighbors,"%2.0f")))
            fput(fh,label+": "+cq_key(T.states[k,])+"; correct="+strofreal(matches[k])+"; distance="+strofreal(T.depth[k]))
            fput(fh,"One shortest solution: "+cq_names(solution))
            fput(fh,"Correct colors along that solution: "+invtokens(strofreal(route[,2]',"%2.0f")," -> "))
            fput(fh,"All nine neighbor scores: "+invtokens(strofreal(neighbors,"%2.0f")))
            st_local("cq_title"+strofreal(j+1),label+": "+strofreal(matches[k])+"/24 correct; minimum "+strofreal(T.depth[k])+" moves")
            st_local("cq_sub"+strofreal(j+1),"One shortest solution: "+cq_names(solution)+"   |   Correct: "+invtokens(strofreal(route[,2]',"%2.0f")," -> "))
        }
        fclose(fh)
    }
    
    end
    }
    
    capture erase "closeness_results.txt"
    mata: cq_experiment()
    local cq_title1 "Target: 24/24 correct; 0 moves"
    local cq_sub1 "The same fixed target and the same two camera views for every row"
    local panels
    forvalues j=1/3 {
        clear
        mata: cq_render(st_matrix("CQ_chosen")[`j',],cq_geometry())
        local colors `""40 45 53" "204 57 57" "239 145 50" "244 244 235" "246 211 64" "49 151 113" "66 111 185""'
        quietly levelsof polygon, local(polys)
        local plots
        foreach p of local polys {
            quietly summarize color if polygon==`p', meanonly
            local ci=r(min)+1
            local rgb : word `ci' of `colors'
            local plots `"`plots' (area py px if polygon==`p', color("`rgb'") lcolor("40 45 53") lwidth(vthin))"'
        }
        twoway `plots', aspectratio(.5070422535) xscale(range(-1.65 5.45) noextend off) yscale(range(-1.6 2.0) noextend off) xlabel(none) ylabel(none) legend(off) text(1.75 0 "Front / right / up", size(small) color(gs6)) text(1.75 3.8 "Back / left / down", size(small) color(gs6)) title("`cq_title`j''", size(medium)) subtitle("`cq_sub`j''", size(small)) graphregion(color(white)) plotregion(color(white)) xsize(9) ysize(5) name(cq_panel`j', replace)
        local panels `panels' cq_panel`j'
    }
    graph combine `panels', cols(1) imargin(tiny) graphregion(color(white)) title("More correct colors, fewer moves?", size(medium)) note("R/right, U/up, F/front; prime = counterclockwise from outside; 2 = half turn." "Each turn costs one. Distance is exact under these moves and this fixed target.", size(vsmall)) xsize(9) ysize(14) name(cq_comparison, replace)
    graph export "cube_closeness.png", width(1500) replace
    graph save "cube_closeness.gph", replace
    clear
    svmat double CQ_census
    rename CQ_census1 distance
    rename CQ_census2 correct
    forvalues j=1/24 {
        rename CQ_census`=`j'+2' s`j'
    }
    assert _N==385
    export delimited using "closeness_states.csv", replace
    contract distance correct, freq(states)
    sort distance correct
    list, noobs sepby(distance)
    export delimited using "closeness_counts.csv", replace
    restore
    display "CUBE_CLOSENESS_PASS"

    Comment


    • #3
      Did the cube change, or just my view? A Stata/Mata experiment

      In the previous experiment, a cube with more correct colors needed more moves. My next question was simpler: if I change the viewing direction, does anything about the puzzle itself change?

      I kept state A, moved the camera 90 degrees around the vertical axis, and compared that with actually turning the right face by 180 degrees. Each row of the picture uses two opposite views to show all six faces.

      Code:
                                    Correct colors   Minimum moves
      Original state A                   16/24             3
      Same A, different camera            16/24             3
      A after a right-face half turn      12/24             2
      The first two rows contain exactly the same 24 colors at the same model positions. Only their projection onto the screen changes. The third row changes the color arrangement. Mata checks distances by breadth-first search from the fixed target, covering all 385 states within three moves.

      Here is a small coordinate check from the same code. A point stays at (1,-1,1) in the model, while its projected coordinates in the main view change:

      Code:
      Camera angle    Model coordinates     Screen coordinates
           0             (1,-1,1)              ( 0.20,-0.04)
          90             (1,-1,1)              (-1.40, 0.68)
      These are plotting coordinates, not pixels. The code obtains them by taking dot products with the camera's right and up directions; it uses orthographic projection, without perspective. Camera parameters never enter the search.

      R, U and F continue to name the same model faces even when the view changes. Only those faces can turn; a quarter turn, inverse quarter turn or half turn each costs one move. The target orientation stays fixed. Thus "distance" here means moves to that target, not distance on the screen.

      This makes me curious about the next distinction: could moving the object and moving the coordinate origin produce the same-looking picture? Other small examples are welcome.
      Click image for larger version

Name:	cube_camera.png
Views:	1
Size:	496.8 KB
ID:	1787799

      Paste the whole block below into one do-file (Stata 16.1). No packages or input files are required. It preserves the current dataset, uses cv_/CV_ names, and writes cube_camera.png/.gph, camera_results.txt, camera_metrics.csv and camera_point_coordinates.csv to the current folder, replacing those files if present. All searching and drawing use Stata/Mata.

      Code:
      version 16.1
      set more off
      preserve
      capture quietly mata: cv_goal()
      if _rc {
      mata:
      struct cv_tree {
          real matrix states
          real colvector depth, parent, move
          transmorphic scalar seen
      }
      real matrix cv_geometry()
      {
          real matrix N,G
          real rowvector n,a,b,c
          real scalar f,i,j
          N=(1,0,0\-1,0,0\0,0,1\0,0,-1\0,-1,0\0,1,0)
          G=J(0,7,.)
          for(f=1;f<=6;f++) {
              n=N[f,]
              if(abs(n[1])) { a=(0,1,0); b=(0,0,1); }
              else if(abs(n[2])) { a=(1,0,0); b=(0,0,1); }
              else { a=(1,0,0); b=(0,1,0); }
              for(i=-1;i<=1;i=i+2) for(j=-1;j<=1;j=j+2) {
                  c=n+.5*i*a+.5*j*b
                  G=G\(c,n,f)
             }
         }
          return(G)
      }
      
      real matrix cv_moves(real matrix G)
      {
          real matrix P,Q,H
          real rowvector axis,ix,base
          real scalar f,i,j,t,k
          P=J(9,24,.)
          // Row vectors multiply Q. Clockwise seen from outside R, U, F.
          for(f=1;f<=3;f++) {
              if(f==1) { Q=(1,0,0\0,0,-1\0,1,0); axis=(1,0,0); }
              if(f==2) { Q=(0,-1,0\1,0,0\0,0,1); axis=(0,0,1); }
              if(f==3) { Q=(0,0,-1\0,1,0\1,0,0); axis=(0,-1,0); }
              H=G
              for(i=1;i<=24;i++) if(G[i,1..3]*axis'>0) {
                  H[i,1..3]=G[i,1..3]*Q
                  H[i,4..6]=G[i,4..6]*Q
             }
              k=3*f-2
              for(i=1;i<=24;i++) {
                  ix=selectindex(rowsum(abs(G[,1..6]:-H[i,1..6])):==0)'
                  assert(cols(ix)==1)
                  P[k,ix[1]]=i
             }
              base=P[k,]
              P[k+1,]=base[base]
              P[k+2,]=base[base[base]]
              assert(all(sort(base',1):==(1::24)))
              assert(all(base[base[base[base]]]:==(1..24)))
         }
          return(P)
      }
      
      real rowvector cv_goal()
      {
          return((1,1,1,1,2,2,2,2,3,3,3,3,4,4,4,4,5,5,5,5,6,6,6,6))
      }
      string scalar cv_key(real rowvector s)
      {
          return(invtokens(strofreal(s,"%1.0f"),""))
      }
      real scalar cv_inverse(real scalar m)
      {
          return(3*floor((m-1)/3)+4-mod(m-1,3)-1)
      }
      string scalar cv_names(real rowvector path)
      {
          string rowvector names
          names=("R","R2","R'","U","U2","U'","F","F2","F'")
          if(cols(path)==0) return("(none)")
          return(invtokens(names[path]))
      }
      real rowvector cv_apply(real rowvector s,real rowvector path,real matrix P)
      {
          real scalar i
          real rowvector out
          out=s
          for(i=1;i<=cols(path);i++) out=out[P[path[i],]]
          return(out)
      }
      
      struct cv_tree scalar cv_bfs(real rowvector root,real matrix P,real scalar radius)
      {
          struct cv_tree scalar T
          real scalar cap,head,n,d,m
          real rowvector s
          string scalar key
          // Upper bound on all 9-branch paths. This demo limits radius to 4.
          cap=(9^(radius+1)-1)/8
          T.states=J(cap,24,.)
          T.depth=J(cap,1,.)
          T.parent=J(cap,1,0)
          T.move=J(cap,1,0)
          T.seen=asarray_create("string",1)
          T.states[1,]=root; T.depth[1]=0
          asarray(T.seen,cv_key(root),1)
          head=1; n=1
          while(head<=n) {
              d=T.depth[head]
              if(d<radius) for(m=1;m<=9;m++) {
                  s=T.states[head,P[m,]]
                  key=cv_key(s)
                  if(!asarray_contains(T.seen,key)) {
                      n++
                      T.states[n,]=s; T.depth[n]=d+1
                      T.parent[n]=head; T.move[n]=m
                      asarray(T.seen,key,n)
                 }
             }
              head++
         }
          T.states=T.states[1..n,]
          T.depth=T.depth[1..n]; T.parent=T.parent[1..n]; T.move=T.move[1..n]
          return(T)
      }
      real rowvector cv_trace(struct cv_tree scalar T, real scalar i)
      {
          real rowvector path
          real scalar node
          node=i
          path=J(1,0,.)
          while(T.parent[node]>0) {
              path=T.move[node],path
              node=T.parent[node]
         }
          return(path)
      }
      
      
      // Rows are screen-right, screen-up and the direction toward the camera.
      real matrix cv_camera(real scalar degrees)
      {
          real scalar a
          real matrix Q,C
          a=degrees*pi()/180
          Q=(cos(a),sin(a),0\-sin(a),cos(a),0\0,0,1)
          C=(.8,.6,0\-.36,.48,.8\.48,-.64,.6)*Q
          assert(max(abs(C*C'-I(3)))<1e-12)
          return(C)
      }
      
      // Orthographic projection; two opposite views show all 24 stickers.
      real matrix cv_polygons(real rowvector state,real matrix G,real scalar angle)
      {
          real matrix out,pts,C
          real rowvector right,up,cam,n,a,b,c
          real scalar view,i,poly,shift
          out=J(0,5,.); poly=0
          C=cv_camera(angle)
          for(view=1;view<=2;view++) {
              right=C[1,]; up=C[2,]; cam=C[3,]
              if(view==2) { right=-right; cam=-cam; }
              shift=(view==1 ? 0 : 3.8)
              for(i=1;i<=24;i++) {
                  n=G[i,4..6]
                  if(n*cam'<=0) continue
                  if(abs(n[1])) { a=(0,1,0); b=(0,0,1); }
                  else if(abs(n[2])) { a=(1,0,0); b=(0,0,1); }
                  else { a=(1,0,0); b=(0,1,0); }
                  if(mod(i-1,4)==0) {
                      pts=(n-a-b\n+a-b\n+a+b\n-a+b\n-a-b)
                      poly++
                      out=out\(J(5,1,poly),pts*right':+shift,pts*up',J(5,1,0),J(5,1,view))
                  }
                  c=G[i,1..3]
                  pts=(c-.46*a-.46*b\c+.46*a-.46*b\c+.46*a+.46*b\c-.46*a+.46*b\c-.46*a-.46*b)
                  poly++
                  out=out\(J(5,1,poly),pts*right':+shift,pts*up',J(5,1,state[i]),J(5,1,view))
              }
          }
          return(out)
      }
      void cv_render(real rowvector state,real matrix G,real scalar angle)
      {
          real matrix out
          out=cv_polygons(state,G,angle)
          st_addobs(rows(out))
          st_addvar("double",("polygon","px","py","color","view"))
          st_store(.,.,out)
      }
      void cv_experiment()
      {
          struct cv_tree scalar T
          real matrix G,P,states,metrics,point,C,D,E,G0
          real rowvector A,turned,angles,back,path,world,A0
          real scalar j,k,i,dist,correct,fh
          G=cv_geometry(); P=cv_moves(G)
          // Same state A as the preceding experiment; only its colors are input.
          A=(2,1,2,1,2,1,2,1,4,4,3,3,4,4,3,3,5,5,5,5,6,6,6,6)
          turned=cv_apply(A,2,P) // R2: a physical half turn of the right face.
          states=A\A\turned
          angles=(0,90,0)
          T=cv_bfs(cv_goal(),P,3)
          metrics=J(3,5,.)
          fh=fopen("camera_results.txt","w")
          fput(fh,"Fixed target and moves R R2 R' U U2 U' F F2 F'; each costs one.")
          for(j=1;j<=3;j++) {
              assert(asarray_contains(T.seen,cv_key(states[j,])))
              k=asarray(T.seen,cv_key(states[j,]))
              dist=T.depth[k]; correct=sum(states[j,]:==cv_goal())
              metrics[j,]=(j,angles[j],sum(states[j,]:!=A),correct,dist)
              back=cv_trace(T,k); path=J(1,0,.)
              for(i=cols(back);i>=1;i--) path=path,cv_inverse(back[i])
              assert(all(cv_apply(states[j,],path,P):==cv_goal()))
              printf("Panel %g: camera=%g; changed colors=%g; correct=%g/24; minimum=%g\n",j,angles[j],metrics[j,3],correct,dist)
              fput(fh,"Panel "+strofreal(j)+": "+cv_key(states[j,])+"; camera="+strofreal(angles[j])+"; correct="+strofreal(correct)+"; distance="+strofreal(dist))
              fput(fh,"One shortest solution: "+cv_names(path))
              st_local("cv_sub"+strofreal(j),strofreal(correct)+"/24 correct; minimum "+strofreal(dist)+" moves | Solution: "+cv_names(path))
          }
          // Rendering must neither permute the colors nor alter world geometry.
          A0=A; G0=G
          D=cv_polygons(A,G,0); E=cv_polygons(A,G,90)
          assert(all(A:==A0) & all(G:==G0))
          assert(max(abs(D[,2..3]-E[,2..3]))>0.1)
          assert(all(states[1,]:==states[2,]) & !all(states[1,]:==states[3,]))
          assert(all(metrics[,4..5]:==(16,3\16,3\12,2)))
          // A fixed world point; coordinates before the display's second-view offset.
          world=(1,-1,1)
          point=J(2,6,.)
          for(j=1;j<=2;j++) {
              C=cv_camera(angles[j])
              point[j,]=(angles[j],world,world*C[1,]',world*C[2,]')
          }
          st_matrix("CV_states",states)
          st_matrix("CV_metrics",metrics)
          st_matrix("CV_point",point)
          fput(fh,"World point P=(1,-1,1); main-view screen coordinates:")
          for(j=1;j<=2;j++) fput(fh,strofreal(angles[j])+" degrees: "+invtokens(strofreal(point[j,5..6],"%7.3f")))
          fclose(fh)
      }
      
      end
      }
      
      capture erase "camera_results.txt"
      mata: cv_experiment()
      local cv_title1 "1. State A: original camera"
      local cv_title2 "2. Same state A: camera moved 90 degrees around the cube"
      local cv_title3 "3. Right face turned 180 degrees: original camera"
      local panels
      forvalues j=1/3 {
          local angle=cond(`j'==2,90,0)
          local left "Front / right / up"
          local right "Back / left / down"
          if `j'==2 {
              local left "Right / back / up"
              local right "Left / front / down"
          }
          clear
          mata: cv_render(st_matrix("CV_states")[`j',],cv_geometry(),`angle')
          local colors `""40 45 53" "204 57 57" "239 145 50" "244 244 235" "246 211 64" "49 151 113" "66 111 185""'
          quietly levelsof polygon, local(polys)
          local plots
          foreach p of local polys {
              quietly summarize color if polygon==`p', meanonly
              local ci=r(min)+1
              local rgb : word `ci' of `colors'
              local plots `"`plots' (area py px if polygon==`p', nodropbase color("`rgb'") lcolor("40 45 53") lwidth(vthin))"'
          }
          twoway `plots', aspectratio(.5070422535) xscale(range(-1.65 5.45) noextend off) yscale(range(-1.6 2.0) noextend off) xlabel(none) ylabel(none) legend(off) text(1.75 0 "`left'", size(small) color(gs6)) text(1.75 3.8 "`right'", size(small) color(gs6)) title("`cv_title`j''", size(medium)) subtitle("`cv_sub`j''", size(small)) graphregion(color(white)) plotregion(color(white)) xsize(9) ysize(5) name(cv_panel`j', replace)
          local panels `panels' cv_panel`j'
      }
      graph combine `panels', cols(1) imargin(tiny) graphregion(color(white)) title("Did the cube change, or just my view?", size(medium)) note("Each row shows all six faces. Counts use all 24 positions in the same fixed target." "R/right, U/up, F/front. 2 = half turn; each turn costs one. Face names stay fixed to the model.", size(vsmall)) xsize(10) ysize(14) name(cv_comparison, replace)
      graph export "cube_camera.png", width(1500) replace
      graph save "cube_camera.gph", replace
      clear
      matrix colnames CV_metrics = panel camera_degrees changed_colors correct minimum_moves
      svmat double CV_metrics, names(col)
      list, noobs
      export delimited using "camera_metrics.csv", replace
      clear
      matrix colnames CV_point = camera_degrees world_x world_y world_z screen_x screen_y
      svmat double CV_point, names(col)
      format screen* %6.2f
      list, noobs
      export delimited using "camera_point_coordinates.csv", replace
      restore
      display "CUBE_CAMERA_PASS"

      Comment


      • #4
        Move the cube, or move the origin? Another small Stata/Mata experiment

        After changing the viewing direction, I wondered whether moving an object and moving the coordinate origin could produce the same-looking picture.

        I added a fixed landmark L to the cube. C denotes its center and K the camera position. The figure compares four cases:

        A. Original scene.
        B. Move only the cube by +2 along x.
        C. Move the coordinate origin by -2 along x, converting the object, landmark AND camera coordinates.
        D. Use C's new object and landmark coordinates, but retain the old camera numbers.

        Mata compares every projected polygon vertex. A and C give the same scene. B and D give the same cube image, but the landmark appears in different places:

        Code:
        Case    3D distance from cube center to landmark
         A                         3
         B                         5
         C                         3
         D                         3
        The distinction becomes clearer in one line. With origin displacement o, a point P becomes P-o and the camera K becomes K-o:

        Code:
        (P-o) - (K-o) = P-K
        The camera-relative position is unchanged. Moving only the object instead gives (P+t)-K. In D, keeping the camera's old numbers in the new frame actually changes its physical position. It is no longer just a relabeling of the original scene.

        So the cube image alone was not enough to tell what moved. The landmark made the difference visible. No face turns occur here: all four cases retain the same sticker arrangement. This experiment uses orthographic projection with fixed camera orientation and plotting limits.
        Click image for larger version

Name:	cube_origin.png
Views:	1
Size:	342.7 KB
ID:	1787801

        I am now curious about changing coordinate units: if the axes are rescaled, what must happen to the distance calculation? Other small examples are welcome.

        Paste the entire block into one do-file (Stata 16.1). No packages or input files are needed. It preserves the current dataset, uses co_/CO_ names, and writes cube_origin.png/.gph and origin_metrics.csv in the current folder, replacing existing files with those names. All geometry and drawing use Stata/Mata.

        Code:
        version 16.1
        set more off
        preserve
        capture quietly mata: co_basis()
        if _rc {
        mata:
        real matrix co_geometry()
        {
            real matrix N,G
            real rowvector n,a,b,c
            real scalar f,i,j
            N=(1,0,0\-1,0,0\0,0,1\0,0,-1\0,-1,0\0,1,0)
            G=J(0,7,.)
            for(f=1;f<=6;f++) {
                n=N[f,]
                if(abs(n[1])) { a=(0,1,0); b=(0,0,1); }
                else if(abs(n[2])) { a=(1,0,0); b=(0,0,1); }
                else { a=(1,0,0); b=(0,1,0); }
                for(i=-1;i<=1;i=i+2) for(j=-1;j<=1;j=j+2) {
                    c=n+.5*i*a+.5*j*b
                    G=G\(c,n,f)
               }
           }
            return(G)
        }
        
        
        real matrix co_basis()
        {
            return((.8,.6,0\-.36,.48,.8\.48,-.64,.6))
        }
        real matrix co_project(real matrix points, real rowvector camera)
        {
            real matrix B
            B=co_basis()
            return((points:-camera)*B[1..2,]')
        }
        real matrix co_polygons(real rowvector state, real rowvector center, real rowvector camera)
        {
            real matrix G,B,out,pts
            real rowvector n,a,b,c
            real scalar i,poly
            G=co_geometry(); B=co_basis()
            out=J(0,4,.); poly=0
            for(i=1;i<=24;i++) {
                n=G[i,4..6]
                if(n*B[3,]'<=0) continue
                if(abs(n[1])) { a=(0,1,0); b=(0,0,1); }
                else if(abs(n[2])) { a=(1,0,0); b=(0,0,1); }
                else { a=(1,0,0); b=(0,1,0); }
                if(mod(i-1,4)==0) {
                    pts=(n-a-b\n+a-b\n+a+b\n-a+b\n-a-b):+center
                    poly++
                    out=out\(J(5,1,poly),co_project(pts,camera),J(5,1,0))
                }
                c=G[i,1..3]
                pts=(c-.46*a-.46*b\c+.46*a-.46*b\c+.46*a+.46*b\c-.46*a+.46*b\c-.46*a-.46*b):+center
                poly++
                out=out\(J(5,1,poly),co_project(pts,camera),J(5,1,state[i]))
            }
            return(out)
        }
        void co_render(real rowvector state,real rowvector center,real rowvector camera)
        {
            real matrix out
            out=co_polygons(state,center,camera)
            st_addobs(rows(out))
            st_addvar("double",("polygon","px","py","color"))
            st_store(.,.,out)
        }
        void co_experiment()
        {
            real matrix B,centers,cameras,landmarks,origins,metrics,xy,scene1,scene2,scene3,scene4
            real rowvector K,L,t,o,state,backup
            real scalar j,dist
            B=co_basis()
            assert(max(abs(B*B'-I(3)))<1e-12)
            K=8*B[3,]; L=(-3,0,0); t=(2,0,0); o=-t
            // A baseline; B translate object; C change origin consistently;
            // D intentionally retain the old camera numbers in the new coordinates.
            centers=(0,0,0)\t\(-o)\(-o)
            cameras=K\K\(K-o)\K
            landmarks=L\L\(L-o)\(L-o)
            origins=(0,0,0)\(0,0,0)\o\o
            state=(2,1,2,1,2,1,2,1,4,4,3,3,4,4,3,3,5,5,5,5,6,6,6,6)
            backup=state
            scene1=co_polygons(state,centers[1,],cameras[1,])
            scene2=co_polygons(state,centers[2,],cameras[2,])
            scene3=co_polygons(state,centers[3,],cameras[3,])
            scene4=co_polygons(state,centers[4,],cameras[4,])
            assert(all(state:==backup))
            // Compare ALL projected polygon vertices, not just the center.
            assert(max(abs(scene1-scene3))<1e-12)
            assert(max(abs(scene2-scene4))<1e-12)
            assert(max(abs(scene1[,2..3]-scene2[,2..3]))>1)
            metrics=J(4,10,.)
            for(j=1;j<=4;j++) {
                xy=co_project(centers[j,]\landmarks[j,],cameras[j,])
                dist=sqrt(sum((centers[j,]-landmarks[j,]):^2))
                metrics[j,]=(j,origins[j,1],centers[j,1],landmarks[j,1],cameras[j,1],xy[1,],xy[2,],dist)
                st_local("co_sub"+strofreal(j),"Cx = "+strofreal(centers[j,1])+"; Lx = "+strofreal(landmarks[j,1])+"; Kx = "+strofreal(cameras[j,1],"%4.2f"))
                st_local("co_gap"+strofreal(j),"3D distance C-L = "+strofreal(dist))
            }
            assert(max(abs(metrics[,10]-(3\5\3\3)))<1e-12)
            assert(max(abs(metrics[1,6..9]-metrics[3,6..9]))<1e-12)
            assert(max(abs(metrics[2,6..7]-metrics[4,6..7]))<1e-12)
            assert(max(abs(metrics[2,8..9]-metrics[4,8..9]))>1)
            st_matrix("CO_centers",centers)
            st_matrix("CO_cameras",cameras)
            st_matrix("CO_landmarks",landmarks)
            st_matrix("CO_origins",origins)
            st_matrix("CO_state",state)
            st_matrix("CO_metrics",metrics)
        }
        
        end
        }
        
        mata: co_experiment()
        local co_title1 "A. Original scene"
        local co_title2 "B. Move only the cube: +2 in x"
        local co_title3 "C. Origin -2: convert all coordinates"
        local co_title4 "D. Origin -2: keep old camera numbers"
        local panels
        forvalues j=1/4 {
            clear
            mata: co_render(st_matrix("CO_state"),st_matrix("CO_centers")[`j',],st_matrix("CO_cameras")[`j',])
            local cx=el(CO_metrics,`j',6)
            local cy=el(CO_metrics,`j',7)
            local lx=el(CO_metrics,`j',8)
            local ly=el(CO_metrics,`j',9)
            local colors `""40 45 53" "204 57 57" "239 145 50" "244 244 235" "246 211 64" "49 151 113" "66 111 185""'
            quietly levelsof polygon, local(polys)
            local plots `"(pci `ly' `lx' `cy' `cx', lcolor(gs7) lpattern(dash))"'
            foreach p of local polys {
                quietly summarize color if polygon==`p', meanonly
                local ci=r(min)+1
                local rgb : word `ci' of `colors'
                local plots `"`plots' (area py px if polygon==`p', nodropbase color("`rgb'") lcolor("40 45 53") lwidth(vthin))"'
            }
            local plots `"`plots' (scatteri `ly' `lx' (12) "L", msymbol(O) mcolor("180 55 50") msize(medium) mlabcolor("180 55 50")) (scatteri `cy' `cx' (3) "C", msymbol(+) mcolor(black) msize(medium) mlabcolor(black))"'
            twoway `plots', aspectratio(.625) xscale(range(-4 4) noextend) yscale(range(-2.5 2.5) noextend) xlabel(-4(2)4, labsize(small)) ylabel(-2(1)2, labsize(small) angle(0)) xtitle("Projected horizontal coordinate", size(small)) ytitle("Projected vertical coordinate", size(small)) legend(off) title("`co_title`j''", size(medium)) subtitle("`co_sub`j''", size(small)) text(2.1 0 "`co_gap`j''", size(small) color(gs6)) graphregion(color(white)) plotregion(color(white)) xsize(8) ysize(6) name(co_panel`j', replace)
            local panels `panels' co_panel`j'
        }
        graph combine `panels', cols(2) imargin(small) graphregion(color(white)) title("Move the cube, or move the origin?", size(medium)) note("C = cube center; L = fixed world landmark; K = camera. Subscripts x denote coordinates in each panel's frame." "Camera orientation, plotting limits and cube colors stay fixed. Ky = -5.12, Kz = 4.80 in every panel." "A and C have the same image. B and D have the same cube image, but different landmark positions.", size(vsmall)) xsize(14) ysize(11) name(co_comparison, replace)
        graph export "cube_origin.png", width(1800) replace
        graph save "cube_origin.gph", replace
        clear
        matrix colnames CO_metrics = panel origin_x center_x landmark_x camera_x screen_cx screen_cy screen_lx screen_ly distance_CL
        svmat double CO_metrics, names(col)
        format screen* camera_x %6.2f
        list, noobs
        export delimited using "origin_metrics.csv", replace
        restore
        display "CUBE_ORIGIN_PASS"

        Comment


        • #5
          Did the distance change, or just the numbers? A cube experiment in Stata/Mata

          Changing the origin left the previous scene unchanged. What if I change the units on just one axis?

          I used a cube with sides of 2 m. Between opposite vertices P and Q, the coordinate differences are (2,2,2), so the body diagonal is sqrt(12) m.

          The picture compares four cases. For B, I set u=2*x, leaving the other coordinates unchanged: one u-unit now represents 0.5 m, while one unit on either other axis represents 1 m.

          Code:
          Case                         Coordinate differences   Physical PQ (m)
          A. Original cube                    (2,2,2)                3.464102
          B. Change x units                   (4,2,2)                3.464102
          C. Physically double x length       (4,2,2)                4.898979
          D. Stretch only the display         (2,2,2)                3.464102
          B and C have the same coordinate differences but different physical lengths. For B, simply calculating sqrt(4^2+2^2+2^2) mixes unequal axis units. Converting each component to meters gives:

          Code:
          sqrt((4/2)^2 + 2^2 + 2^2) = sqrt(12)
          More generally, with u=a*x, v=b*y and w=c*z on these perpendicular axes, the distance in the original units is sqrt((du/a)^2+(dv/b)^2+(dw/c)^2), for positive scale factors a, b and c.

          Mata also converts the camera mapping consistently, so B looks like A. C changes the object itself. D doubles the final horizontal display coordinate without changing the stored 3D geometry. The dashed PQ is drawn through the solid; measuring its 2D projection would not give the 3D length, even in A.
          Click image for larger version

Name:	cube_scale.png
Views:	1
Size:	322.8 KB
ID:	1787803

          This leaves me wondering what changes if the coordinate axes are tilted relative to one another, rather than simply rescaled. Suggestions for another small comparison are welcome.

          Paste the whole block below into one do-file (Stata 16.1). No packages or input files are needed. It preserves the current dataset, uses cs_/CS_ names, and writes cube_scale.png/.gph and scale_metrics.csv in the current folder, replacing existing files with those names. All geometry, calculations and drawing use Stata/Mata.

          Code:
          version 16.1
          set more off
          preserve
          capture quietly mata: cs_basis()
          if _rc {
          mata:
          real matrix cs_geometry()
          {
              real matrix N,G
              real rowvector n,a,b,c
              real scalar f,i,j
              N=(1,0,0\-1,0,0\0,0,1\0,0,-1\0,-1,0\0,1,0)
              G=J(0,7,.)
              for(f=1;f<=6;f++) {
                  n=N[f,]
                  if(abs(n[1])) { a=(0,1,0); b=(0,0,1); }
                  else if(abs(n[2])) { a=(1,0,0); b=(0,0,1); }
                  else { a=(1,0,0); b=(0,1,0); }
                  for(i=-1;i<=1;i=i+2) for(j=-1;j<=1;j=j+2) {
                      c=n+.5*i*a+.5*j*b
                      G=G\(c,n,f)
                 }
             }
              return(G)
          }
          
          
          real matrix cs_basis()
          {
              return((.8,.6,0\-.36,.48,.8\.48,-.64,.6))
          }
          
          // physical: local cube -> physical meters (actual object deformation).
          // units: physical meters -> stored coordinate numbers.
          // display: projected coordinates -> display coordinates.
          real matrix cs_project(real matrix localpoints,real rowvector camera,
              real rowvector physical,real rowvector units,real rowvector display)
          {
              real matrix coords,B
              real rowvector kc
              B=cs_basis()
              coords=localpoints*diag(physical)*diag(units)
              kc=camera*diag(units)
              return((coords:-kc)*diag(1:/units)*B[1..2,]'*diag(display))
          }
          real scalar cs_distance(real rowvector difference,real rowvector units)
          {
              // One unit on axis j equals 1/units[j] meters.
              return(sqrt(sum((difference:/units):^2)))
          }
          real matrix cs_polygons(real rowvector state,real rowvector physical,
              real rowvector units,real rowvector display)
          {
              real matrix G,B,out,pts
              real rowvector n,a,b,c,K
              real scalar i,poly
              assert(all(physical:>0) & all(units:>0) & all(display:>0))
              G=cs_geometry(); B=cs_basis(); K=8*B[3,]
              out=J(0,4,.); poly=0
              // Only positive, axis-aligned scales: face-normal directions stay unchanged.
              for(i=1;i<=24;i++) {
                  n=G[i,4..6]
                  if(n*B[3,]'<=0) continue
                  if(abs(n[1])) { a=(0,1,0); b=(0,0,1); }
                  else if(abs(n[2])) { a=(1,0,0); b=(0,0,1); }
                  else { a=(1,0,0); b=(0,1,0); }
                  if(mod(i-1,4)==0) {
                      pts=(n-a-b\n+a-b\n+a+b\n-a+b\n-a-b)
                      poly++
                      out=out\(J(5,1,poly),cs_project(pts,K,physical,units,display),J(5,1,0))
                  }
                  c=G[i,1..3]
                  pts=(c-.46*a-.46*b\c+.46*a-.46*b\c+.46*a+.46*b\c-.46*a+.46*b\c-.46*a-.46*b)
                  poly++
                  out=out\(J(5,1,poly),cs_project(pts,K,physical,units,display),J(5,1,state[i]))
              }
              return(out)
          }
          void cs_render(real rowvector state,real rowvector physical,
              real rowvector units,real rowvector display)
          {
              real matrix out
              out=cs_polygons(state,physical,units,display)
              st_addobs(rows(out))
              st_addvar("double",("polygon","px","py","color"))
              st_store(.,.,out)
          }
          void cs_experiment()
          {
              real matrix W,U,D,metrics,endpoints,xy,world,coords,B,s1,s2,s3,s4
              real rowvector state,P,Q,delta,K
              real scalar j,raw,physical,screen
              W=(1,1,1\1,1,1\2,1,1\1,1,1)
              U=(1,1,1\2,1,1\1,1,1\1,1,1)
              D=(1,1\1,1\1,1\2,1)
              state=(2,1,2,1,2,1,2,1,4,4,3,3,4,4,3,3,5,5,5,5,6,6,6,6)
              P=(-1,-1,-1); Q=(1,1,1)
              B=cs_basis(); K=8*B[3,]
              metrics=J(4,11,.); endpoints=J(4,4,.)
              for(j=1;j<=4;j++) {
                  world=(P\Q)*diag(W[j,])
                  coords=world*diag(U[j,])
                  delta=coords[2,]-coords[1,]
                  raw=sqrt(sum(delta:^2))
                  physical=cs_distance(delta,U[j,])
                  assert(abs(physical-sqrt(sum((world[2,]-world[1,]):^2)))<1e-12)
                  xy=cs_project(P\Q,K,W[j,],U[j,],D[j,])
                  screen=sqrt(sum((xy[2,]-xy[1,]):^2))
                  metrics[j,]=(j,W[j,1],U[j,1],D[j,1],delta,raw,physical,screen,1/U[j,1])
                  endpoints[j,]=(xy[1,],xy[2,])
                  st_local("cs_delta"+strofreal(j),"Coordinate difference: ("+invtokens(strofreal(delta,"%2.0f"),", ")+")")
                  st_local("cs_length"+strofreal(j),"Physical PQ = "+strofreal(physical,"%8.6f")+" m")
              }
              s1=cs_polygons(state,W[1,],U[1,],D[1,])
              s2=cs_polygons(state,W[2,],U[2,],D[2,])
              s3=cs_polygons(state,W[3,],U[3,],D[3,])
              s4=cs_polygons(state,W[4,],U[4,],D[4,])
              assert(max(abs(s1-s2))<1e-12)
              assert(max(abs(s4[,2]-2*s1[,2]))<1e-12)
              assert(max(abs(s4[,3]-s1[,3]))<1e-12)
              assert(max(abs(s3[,2..3]-s4[,2..3]))>.1)
              assert(max(abs(metrics[,9]-(sqrt(12)\sqrt(12)\sqrt(24)\sqrt(12))))<1e-12)
              st_matrix("CS_physical",W)
              st_matrix("CS_units",U)
              st_matrix("CS_display",D)
              st_matrix("CS_state",state)
              st_matrix("CS_metrics",metrics)
              st_matrix("CS_endpoints",endpoints)
          }
          
          end
          }
          
          mata: cs_experiment()
          local cs_title1 "A. Original cube: sides 2 m"
          local cs_title2 "B. Change x units: 1 unit = 0.5 m"
          local cs_title3 "C. Physically double the x length"
          local cs_title4 "D. Stretch only the display horizontally"
          local panels
          forvalues j=1/4 {
              clear
              mata: cs_render(st_matrix("CS_state"),st_matrix("CS_physical")[`j',],st_matrix("CS_units")[`j',],st_matrix("CS_display")[`j',])
              local px=el(CS_endpoints,`j',1)
              local py=el(CS_endpoints,`j',2)
              local qx=el(CS_endpoints,`j',3)
              local qy=el(CS_endpoints,`j',4)
              local colors `""40 45 53" "204 57 57" "239 145 50" "244 244 235" "246 211 64" "49 151 113" "66 111 185""'
              quietly levelsof polygon, local(polys)
              local plots
              foreach p of local polys {
                  quietly summarize color if polygon==`p', meanonly
                  local ci=r(min)+1
                  local rgb : word `ci' of `colors'
                  local plots `"`plots' (area py px if polygon==`p', nodropbase color("`rgb'") lcolor("40 45 53") lwidth(vthin))"'
              }
              local plots `"`plots' (pci `py' `px' `qy' `qx', lcolor(black) lpattern(dash) lwidth(medthin)) (scatteri `py' `px' (9) "P" `qy' `qx' (3) "Q", msymbol(O) mcolor(black) msize(small) mlabcolor(black) mlabsize(medium))"'
              twoway `plots', aspectratio(.7878787879) xscale(range(-3.3 3.3) noextend off) yscale(range(-2.6 2.6) noextend off) xlabel(none) ylabel(none) legend(off) title("`cs_title`j''", size(medium)) subtitle("`cs_delta`j''" "`cs_length`j''", size(small)) graphregion(color(white)) plotregion(color(white)) xsize(8) ysize(7) name(cs_panel`j', replace)
              local panels `panels' cs_panel`j'
          }
          graph combine `panels', cols(2) imargin(small) graphregion(color(white)) title("Did the distance change, or just the numbers?", size(medium)) note("PQ is the body diagonal, drawn through the solid. Its 2D projection is not its 3D length." "B changes units consistently; C deforms the object; D stretches the final display." "All panels keep the same camera orientation and sticker colors.", size(vsmall)) xsize(14) ysize(12) name(cs_comparison, replace)
          graph export "cube_scale.png", width(1800) replace
          graph save "cube_scale.gph", replace
          clear
          matrix colnames CS_metrics = panel physical_x coordinate_x display_x du dv dw raw_norm distance_m display_distance x_unit_m
          svmat double CS_metrics, names(col)
          format raw_norm distance_m display_distance %9.6f
          list panel du dv dw raw_norm distance_m, noobs
          export delimited using "scale_metrics.csv", replace
          restore
          display "CUBE_SCALE_PASS"

          Comment


          • #6
            Same cube, tilted axes: is the usual sum of squares enough?

            After changing coordinate units, I wondered what would happen if the axes were no longer perpendicular.

            I kept the same cube, with sides of 2 m, and measured the body diagonal PQ. This time every coordinate basis vector has length 1 m: only their directions change.

            A uses the original perpendicular axes. B rotates the entire coordinate frame by 45 degrees, preserving right angles. C uses u and v axes separated by 60 degrees, with w perpendicular to both. D instead deforms the actual object so its two horizontal edge directions meet at 60 degrees, retaining perpendicular coordinates.

            Code:
            Case                           Raw coordinate norm   Physical PQ (m)
            A. Original axes                     3.464102             3.464102
            B. Rotate the whole frame            3.464102             3.464102
            C. Axes separated by 60 degrees       3.169837             3.464102
            D. Deform the object                 4.000000             4.000000
            A-C are the same physical cube with the same camera. For C, the coordinate difference becomes approximately (0.845299, 2.309401, 2), although the physical displacement is still (2,2,2) m in the original frame.

            The ordinary sum of squared coordinate differences assumes perpendicular unit basis vectors. For the unit-length axes in C, the missing contribution is a cross term:

            Code:
            d^2 = du^2 + dv^2 + dw^2 + 2*cos(60 degrees)*du*dv
                = du^2 + dv^2 + dw^2 + du*dv
            This restores sqrt(12) m. The extra term comes from the angle between the coordinate directions; it is not an adjustment to the observations.

            In the Mata code, rows of E contain the basis vectors in physical space. A coordinate row vector q represents q*E. With delta the coordinate difference, the distance is sqrt(delta*(E*E')*delta'). For B, E*E' is the identity; for C, its two off-diagonal u-v entries are 0.5.

            The blue insets show the coordinate directions in the xy-plane, not their camera projections. The dashed diagonal is drawn through the solid. Other small checks are welcome.
            Click image for larger version

Name:	cube_oblique.png
Views:	1
Size:	388.1 KB
ID:	1787805

            Paste the entire block into one do-file (Stata 16.1). No packages or input files are needed. It preserves the current dataset, uses ct_/CT_ names, and writes cube_oblique.png/.gph and oblique_metrics.csv in the current folder, replacing files with those names. All calculations and drawing use Stata/Mata.

            Code:
            version 16.1
            set more off
            preserve
            capture quietly mata: ct_oblique()
            if _rc {
            mata:
            real matrix ct_geometry()
            {
                real matrix N,G
                real rowvector n,a,b,c
                real scalar f,i,j
                N=(1,0,0\-1,0,0\0,0,1\0,0,-1\0,-1,0\0,1,0)
                G=J(0,7,.)
                for(f=1;f<=6;f++) {
                    n=N[f,]
                    if(abs(n[1])) { a=(0,1,0); b=(0,0,1); }
                    else if(abs(n[2])) { a=(1,0,0); b=(0,0,1); }
                    else { a=(1,0,0); b=(0,1,0); }
                    for(i=-1;i<=1;i=i+2) for(j=-1;j<=1;j=j+2) {
                        c=n+.5*i*a+.5*j*b
                        G=G\(c,n,f)
                   }
               }
                return(G)
            }
            
            
            real matrix ct_basis()
            {
                return((.8,.6,0\-.36,.48,.8\.48,-.64,.6))
            }
            
            
            // Rows of E are the coordinate basis vectors in physical meters.
            // Coordinate row vector q represents the physical point q*E.
            real matrix ct_oblique()
            {
                return((1,0,0\.5,sqrt(3)/2,0\0,0,1))
            }
            real matrix ct_frame(real scalar panel)
            {
                real scalar a
                a=sqrt(.5)
                if(panel==2) return((a,a,0\-a,a,0\0,0,1))
                if(panel==3) return(ct_oblique())
                return(I(3))
            }
            real matrix ct_deform(real scalar panel)
            {
                if(panel==4) return(ct_oblique())
                return(I(3))
            }
            real matrix ct_project(real matrix localpoints,real scalar panel)
            {
                real matrix E,W,B,q
                real rowvector K,kq
                E=ct_frame(panel); W=ct_deform(panel); B=ct_basis()
                K=8*B[3,]
                q=localpoints*W*luinv(E)
                kq=K*luinv(E)
                return((q:-kq)*E*B[1..2,]')
            }
            real scalar ct_distance(real rowvector delta,real matrix E)
            {
                real matrix metric
                metric=E*E'
                return(sqrt(delta*metric*delta'))
            }
            real matrix ct_polygons(real rowvector state,real scalar panel)
            {
                real matrix G,B,W,out,pts
                real rowvector n,nworld,a,b,c
                real scalar i,poly
                G=ct_geometry(); B=ct_basis(); W=ct_deform(panel)
                out=J(0,4,.); poly=0
                for(i=1;i<=24;i++) {
                    n=G[i,4..6]
                    // Transform face normals for the ACTUAL deformation in panel D.
                    nworld=n*luinv(W)'
                    if(nworld*B[3,]'<=0) continue
                    if(abs(n[1])) { a=(0,1,0); b=(0,0,1); }
                    else if(abs(n[2])) { a=(1,0,0); b=(0,0,1); }
                    else { a=(1,0,0); b=(0,1,0); }
                    if(mod(i-1,4)==0) {
                        pts=(n-a-b\n+a-b\n+a+b\n-a+b\n-a-b)
                        poly++
                        out=out\(J(5,1,poly),ct_project(pts,panel),J(5,1,0))
                    }
                    c=G[i,1..3]
                    pts=(c-.46*a-.46*b\c+.46*a-.46*b\c+.46*a+.46*b\c-.46*a+.46*b\c-.46*a-.46*b)
                    poly++
                    out=out\(J(5,1,poly),ct_project(pts,panel),J(5,1,state[i]))
                }
                return(out)
            }
            void ct_render(real rowvector state,real scalar panel)
            {
                real matrix out
                out=ct_polygons(state,panel)
                st_addobs(rows(out))
                st_addvar("double",("polygon","px","py","color"))
                st_store(.,.,out)
            }
            void ct_experiment()
            {
                real matrix metrics,endpoints,axes,E,W,world,q,xy,m,s1,sj
                real rowvector state,P,Q,delta
                real scalar j,raw,physical,cross
                state=(2,1,2,1,2,1,2,1,4,4,3,3,4,4,3,3,5,5,5,5,6,6,6,6)
                P=(-1,-1,-1); Q=(1,1,1)
                metrics=J(4,9,.); endpoints=J(4,4,.); axes=J(4,4,.)
                for(j=1;j<=4;j++) {
                    E=ct_frame(j); W=ct_deform(j)
                    world=(P\Q)*W
                    q=world*luinv(E); delta=q[2,]-q[1,]
                    m=E*E'
                    raw=sqrt(sum(delta:^2))
                    physical=ct_distance(delta,E)
                    cross=2*m[1,2]*delta[1]*delta[2]
                    assert(abs(physical-sqrt(sum((world[2,]-world[1,]):^2)))<1e-12)
                    metrics[j,]=(j,delta,raw,physical,m[1,2],cross,physical^2-raw^2)
                    xy=ct_project(P\Q,j)
                    endpoints[j,]=(xy[1,],xy[2,])
                    axes[j,]=(E[1,1..2],E[2,1..2])
                    st_local("ct_delta"+strofreal(j),"Coordinate difference: ("+invtokens(strtrim(strofreal(delta,"%6.3f")),", ")+")")
                    st_local("ct_length"+strofreal(j),"Raw norm = "+strofreal(raw,"%8.6f")+"; physical PQ = "+strofreal(physical,"%8.6f")+" m")
                }
                s1=ct_polygons(state,1)
                for(j=2;j<=3;j++) {
                    sj=ct_polygons(state,j)
                    assert(max(abs(s1-sj))<1e-12)
                }
                sj=ct_polygons(state,4)
                assert(max(abs(s1[,2..3]-sj[,2..3]))>.1)
                assert(max(abs(metrics[,6]-(sqrt(12)\sqrt(12)\sqrt(12)\4)))<1e-12)
                assert(abs(metrics[3,8]-metrics[3,9])<1e-12)
                st_matrix("CT_state",state)
                st_matrix("CT_metrics",metrics)
                st_matrix("CT_endpoints",endpoints)
                st_matrix("CT_axes",axes)
            }
            
            end
            }
            
            mata: ct_experiment()
            local ct_title1 "A. Original perpendicular axes"
            local ct_title2 "B. Rotate the coordinate frame 45 degrees"
            local ct_title3 "C. Coordinate axes at 60 degrees"
            local ct_title4 "D. Deform the object to 60-degree edges"
            local panels
            forvalues j=1/4 {
                clear
                mata: ct_render(st_matrix("CT_state"),`j')
                local px=el(CT_endpoints,`j',1)
                local py=el(CT_endpoints,`j',2)
                local qx=el(CT_endpoints,`j',3)
                local qy=el(CT_endpoints,`j',4)
                local ux=-2.5+.5*el(CT_axes,`j',1)
                local uy=-2.25+.5*el(CT_axes,`j',2)
                local vx=-2.5+.5*el(CT_axes,`j',3)
                local vy=-2.25+.5*el(CT_axes,`j',4)
                local colors `""40 45 53" "204 57 57" "239 145 50" "244 244 235" "246 211 64" "49 151 113" "66 111 185""'
                quietly levelsof polygon, local(polys)
                local plots
                foreach p of local polys {
                    quietly summarize color if polygon==`p', meanonly
                    local ci=r(min)+1
                    local rgb : word `ci' of `colors'
                    local plots `"`plots' (area py px if polygon==`p', nodropbase color("`rgb'") lcolor("40 45 53") lwidth(vthin))"'
                }
                local plots `"`plots' (pci `py' `px' `qy' `qx', lcolor(black) lpattern(dash)) (scatteri `py' `px' (9) "P" `qy' `qx' (3) "Q", msymbol(O) mcolor(black) msize(small) mlabcolor(black) mlabsize(medium)) (pci -2.25 -2.5 `uy' `ux' -2.25 -2.5 `vy' `vx', lcolor("40 104 150") lwidth(medthick)) (scatteri `uy' `ux' (3) "u" `vy' `vx' (12) "v", msymbol(none) mlabcolor("40 104 150") mlabsize(small))"'
                twoway `plots', aspectratio(.875) xscale(range(-3.2 3.2) noextend off) yscale(range(-2.8 2.8) noextend off) xlabel(none) ylabel(none) legend(off) title("`ct_title`j''", size(medium)) subtitle("`ct_delta`j''" "`ct_length`j''", size(small)) text(-2.25 1 "Axes in the xy-plane" "(w perpendicular)", size(vsmall) color(gs6)) graphregion(color(white)) plotregion(color(white)) xsize(8) ysize(7) name(ct_panel`j', replace)
                local panels `panels' ct_panel`j'
            }
            graph combine `panels', cols(2) imargin(small) graphregion(color(white)) title("If the axes are not perpendicular, how long is the diagonal?", size(medium)) note("A-C describe the same physical cube with the same camera. Every coordinate basis vector has length 1 m." "C needs a cross term in the distance calculation. D changes the object while keeping perpendicular coordinates." "Blue insets show basis directions in the xy-plane, not the camera projection. PQ is drawn through the solid.", size(vsmall)) xsize(14) ysize(12) name(ct_comparison, replace)
            graph export "cube_oblique.png", width(1800) replace
            graph save "cube_oblique.gph", replace
            clear
            matrix colnames CT_metrics = panel du dv dw raw_norm distance_m uv_inner cross_term d2_minus_raw2
            svmat double CT_metrics, names(col)
            format du dv dw raw_norm distance_m cross_term d2_minus_raw2 %9.6f
            list panel du dv dw raw_norm distance_m cross_term, noobs
            export delimited using "oblique_metrics.csv", replace
            restore
            display "CUBE_OBLIQUE_PASS"

            Comment


            • #7
              Can a deformed cube always be recovered? A final Stata/Mata experiment

              After changing origins, units and axis angles, I wanted to try the reverse journey: deform the cube, then recover its original coordinates.

              The transformation is x=u+v*cos(theta), y=v*sin(theta), z=w. For positive angles up to 90 degrees it has an inverse. Mata recovers coordinates by solving the corresponding linear system with lusolve(); the recovery function receives the observed coordinates and the known angle, not the original cube.

              To test sensitivity, I deliberately add +1e-8 m to deformed y coordinates on the original v=+1 side, and -1e-8 m on the v=-1 side. This controlled input error is the same at every positive angle. It is not an error automatically introduced by Stata.

              For the eight vertices of a cube with sides of 2 m, Stata/MP 16.1 gives:

              Code:
              Angle (degrees)    Maximum vertex error after recovery (m)
                                   Original coordinates    Perturbed coordinates
              60                       1.11e-16                 1.29e-8
              0.0001                   1.11e-16                 0.008103
              0.000001                 0                        0.810285
              0                        No unique inverse; recovery not attempted
              The zero in the third row is the result for these particular unperturbed vertices in this run, not a guarantee of exact arithmetic. The experiment does not say that a nearly singular transformation must always fail.

              For positive theta, the inverse contains v=y/sin(theta) and u=x-v*cos(theta). A small error in y is therefore strongly amplified as sin(theta) approaches zero. For our chosen perturbation pattern, the vertex error is epsilon*sqrt(1+cos(theta)^2)/sin(theta), matching the numerical results.

              At exactly zero degrees, the situation is different: (1,-1,1) and (-1,1,1) both become (0,0,1). The collapsed coordinates alone cannot identify which original point was supplied. More numerical precision cannot restore that lost distinction; additional assumptions would be needed.

              I found this a useful place to end the cube experiments: a changed picture, a changed coordinate description, an actual deformation, and a sensitive or nonunique inverse are different things. Corrections or simpler examples are welcome.
              Click image for larger version

Name:	cube_recovery.png
Views:	1
Size:	406.3 KB
ID:	1787807

              Paste the whole block into one do-file (Stata 16.1). No packages or input files are needed. It preserves the current dataset, uses cr_/CR_ names, and writes cube_recovery.png/.gph, recovery_metrics.csv and recovery_vertices.csv in the current folder, replacing existing files with those names. The picture shows 60 degrees, 0.000001 degrees and zero; the CSV also includes intermediate angles. All calculation and drawing use Stata/Mata.

              Code:
              version 16.1
              set more off
              preserve
              capture quietly mata: cr_vertices()
              if _rc {
              mata:
              real matrix cr_geometry()
              {
                  real matrix N,G
                  real rowvector n,a,b,c
                  real scalar f,i,j
                  N=(1,0,0\-1,0,0\0,0,1\0,0,-1\0,-1,0\0,1,0)
                  G=J(0,7,.)
                  for(f=1;f<=6;f++) {
                      n=N[f,]
                      if(abs(n[1])) { a=(0,1,0); b=(0,0,1); }
                      else if(abs(n[2])) { a=(1,0,0); b=(0,0,1); }
                      else { a=(1,0,0); b=(0,1,0); }
                      for(i=-1;i<=1;i=i+2) for(j=-1;j<=1;j=j+2) {
                          c=n+.5*i*a+.5*j*b
                          G=G\(c,n,f)
                     }
                 }
                  return(G)
              }
              
              
              real matrix cr_basis()
              {
                  return((.8,.6,0\-.36,.48,.8\.48,-.64,.6))
              }
              
              
              
              real matrix cr_vertices()
              {
                  return((-1,-1,-1\-1,-1,1\-1,1,-1\-1,1,1\1,-1,-1\1,-1,1\1,1,-1\1,1,1))
              }
              real matrix cr_transform(real scalar degrees)
              {
                  real scalar a
                  assert(degrees>=0 & degrees<=90)
                  a=degrees*pi()/180
                  return((1,0,0\cos(a),sin(a),0\0,0,1))
              }
              real matrix cr_recover(real matrix observed,real scalar degrees)
              {
                  real matrix E
                  assert(degrees>0)
                  E=cr_transform(degrees)
                  // The recovery function receives only observed coordinates and the map.
                  return(lusolve(E',observed')')
              }
              real matrix cr_warp(real matrix points,real scalar degrees,
                  real scalar stage,real scalar epsilon)
              {
                  real matrix Y
                  if(stage==1) return(points)
                  Y=points*cr_transform(degrees)
                  if(stage==2) return(Y)
                  // Controlled benchmark: +/- epsilon at vertices with original v=+/-1.
                  // Interpolate that perturbation linearly across the displayed faces.
                  Y[,2]=Y[,2]+epsilon*points[,2]
                  return(cr_recover(Y,degrees))
              }
              real matrix cr_project(real matrix X)
              {
                  real matrix B
                  B=cr_basis()
                  return(X*B[1..2,]')
              }
              real matrix cr_polygons(real rowvector state,real scalar degrees,
                  real scalar stage,real scalar epsilon)
              {
                  real matrix G,B,W,out,pts,xy
                  real rowvector n,nworld,a,b,c
                  real scalar i,poly
                  // At zero degrees the cube collapses to this plane. Do not invert it.
                  if(degrees==0 & stage==2) {
                      pts=(-2,0,-1\2,0,-1\2,0,1\-2,0,1\-2,0,-1)
                      return((J(5,1,1),cr_project(pts),J(5,1,7)))
                  }
                  G=cr_geometry(); B=cr_basis()
                  W=cr_warp(I(3),degrees,stage,epsilon)
                  out=J(0,4,.); poly=0
                  for(i=1;i<=24;i++) {
                      n=G[i,4..6]; nworld=n*luinv(W)'
                      if(nworld*B[3,]'<=0) continue
                      if(abs(n[1])) { a=(0,1,0); b=(0,0,1); }
                      else if(abs(n[2])) { a=(1,0,0); b=(0,0,1); }
                      else { a=(1,0,0); b=(0,1,0); }
                      if(mod(i-1,4)==0) {
                          pts=(n-a-b\n+a-b\n+a+b\n-a+b\n-a-b)
                          poly++
                          xy=cr_project(cr_warp(pts,degrees,stage,epsilon))
                          out=out\(J(5,1,poly),xy,J(5,1,0))
                      }
                      c=G[i,1..3]
                      pts=(c-.46*a-.46*b\c+.46*a-.46*b\c+.46*a+.46*b\c-.46*a+.46*b\c-.46*a-.46*b)
                      poly++
                      xy=cr_project(cr_warp(pts,degrees,stage,epsilon))
                      out=out\(J(5,1,poly),xy,J(5,1,state[i]))
                  }
                  return(out)
              }
              void cr_render(real scalar degrees,real scalar stage,real scalar epsilon)
              {
                  real rowvector state
                  real matrix out
                  state=(2,1,2,1,2,1,2,1,4,4,3,3,4,4,3,3,5,5,5,5,6,6,6,6)
                  out=cr_polygons(state,degrees,stage,epsilon)
                  st_addobs(rows(out))
                  st_addvar("double",("polygon","px","py","color"))
                  st_store(.,.,out)
              }
              void cr_experiment()
              {
                  real matrix V,E,Y,N,clean,recovered,metrics,audit,collision
                  real rowvector angles
                  real scalar i,theta,a,epsilon,cleanerr,noisyerr,expected,condition
                  V=cr_vertices(); epsilon=1e-8
                  angles=(60,10,1,.01,.0001,.000001,0)
                  metrics=J(7,8,.); audit=J(0,17,.)
                  for(i=1;i<=cols(angles);i++) {
                      theta=angles[i]; E=cr_transform(theta); Y=V*E
                      if(theta==0) {
                          metrics[i,]=(theta,0,0,.,.,.,.,.)
                          audit=audit\(J(8,1,theta),(1::8),V,Y,J(8,9,.))
                          continue
                      }
                      N=Y; N[,2]=N[,2]+epsilon*V[,2]
                      clean=cr_recover(Y,theta)
                      recovered=cr_recover(N,theta)
                      cleanerr=max(sqrt(rowsum((clean-V):^2)))
                      noisyerr=max(sqrt(rowsum((recovered-V):^2)))
                      a=theta*pi()/180
                      expected=epsilon*sqrt(1+cos(a)^2)/sin(a)
                      condition=(1+cos(a))/sin(a)
                      assert(abs(noisyerr-expected)<1e-12)
                      assert(cleanerr<1e-12)
                      metrics[i,]=(theta,1,sin(a),condition,epsilon,cleanerr,noisyerr,expected)
                      audit=audit\(J(8,1,theta),(1::8),V,Y,N,clean,recovered)
                  }
                  collision=(1,-1,1\-1,1,1)*cr_transform(0)
                  assert(all(collision[1,]:==collision[2,]))
                  assert(all(collision[1,]:==(0,0,1)))
                  st_matrix("CR_metrics",metrics)
                  st_matrix("CR_audit",audit)
                  st_numscalar("CR_epsilon",epsilon)
                  st_local("cr_cleanmax",strtrim(strofreal(max(metrics[1..6,6]),"%9.3g")))
                  st_local("cr_err1",strtrim(strofreal(metrics[1,7],"%9.3g")))
                  st_local("cr_err2",strtrim(strofreal(metrics[6,7],"%9.6f")))
              }
              
              end
              }
              
              mata: cr_experiment()
              local panels
              forvalues row=1/3 {
                  local angle=cond(`row'==1,60,cond(`row'==2,.000001,0))
                  local label "60 degrees"
                  if `row'==2 local label "0.000001 degrees"
                  if `row'==3 local label "0 degrees"
                  forvalues stage=1/3 {
                      local panel=3*(`row'-1)+`stage'
                      if `row'==3 & `stage'==3 {
                          twoway scatteri 0 0, msymbol(none) xscale(range(-3 3) off) yscale(range(-2.6 2.6) off) xlabel(none) ylabel(none) legend(off) aspectratio(.8666666667) title("No unique inverse", size(medium)) subtitle("Distinct inputs have merged", size(small)) text(1.3 0 "(1, -1, 1)", size(medium)) text(.65 0 "and (-1, 1, 1)", size(medium)) text(-.15 0 "both become", size(medium) color(gs6)) text(-.95 0 "(0, 0, 1)", size(large) color("180 55 50")) graphregion(color(white)) plotregion(color(white)) xsize(6) ysize(5.4) name(cr_panel`panel', replace)
                      }
                      else {
                          clear
                          mata: cr_render(`angle',`stage',st_numscalar("CR_epsilon"))
                          local colors `""40 45 53" "204 57 57" "239 145 50" "244 244 235" "246 211 64" "49 151 113" "66 111 185" "197 210 221""'
                          quietly levelsof polygon, local(polys)
                          local plots
                          foreach p of local polys {
                              quietly summarize color if polygon==`p', meanonly
                              local ci=r(min)+1
                              local rgb : word `ci' of `colors'
                              local plots `"`plots' (area py px if polygon==`p', nodropbase color("`rgb'") lcolor("40 45 53") lwidth(vthin))"'
                          }
                          local title "Original cube"
                          local sub "Test angle: `label'"
                          if `stage'==2 {
                              local title "Deform: `label'"
                              local sub "Before adding coordinate error"
                          }
                          if `stage'==3 {
                              local title "Recovery with input error"
                              local sub "Max vertex error: `cr_err`row'' m"
                          }
                          if `row'==3 & `stage'==2 {
                              local sub "Collapsed onto a plane"
                              local plots `"`plots' (scatteri .8 0 (12) "Same point", msymbol(O) mcolor("180 55 50") msize(medium) mlabcolor("180 55 50"))"'
                          }
                          twoway `plots', aspectratio(.8666666667) xscale(range(-3 3) noextend off) yscale(range(-2.6 2.6) noextend off) xlabel(none) ylabel(none) legend(off) title("`title'", size(medium)) subtitle("`sub'", size(small)) graphregion(color(white)) plotregion(color(white)) xsize(6) ysize(5.4) name(cr_panel`panel', replace)
                      }
                      local panels `panels' cr_panel`panel'
                  }
              }
              graph combine `panels', cols(3) imargin(small) graphregion(color(white)) title("Can a deformed cube always be recovered?", size(medium)) note("At positive angles, add y errors of +/- 1e-8 m on the original v = +/-1 sides." "Input errors are supplied deliberately, not generated by Stata. Faces interpolate the vertices." "Unperturbed max error in this run: `cr_cleanmax' m. At zero, the collapsed coordinates have no unique inverse.", size(vsmall)) xsize(15) ysize(14) name(cr_comparison, replace)
              graph export "cube_recovery.png", width(2000) replace
              graph save "cube_recovery.gph", replace
              clear
              matrix colnames CR_metrics = degrees invertible determinant condition2 input_error clean_error perturbed_error predicted_error
              svmat double CR_metrics, names(col)
              format degrees %12.6g
              format clean_error perturbed_error predicted_error %12.6g
              list degrees invertible input_error clean_error perturbed_error, noobs
              export delimited using "recovery_metrics.csv", replace
              clear
              matrix colnames CR_audit = degrees vertex u v w x y z noisy_x noisy_y noisy_z clean_u clean_v clean_w recovered_u recovered_v recovered_w
              svmat double CR_audit, names(col)
              export delimited using "recovery_vertices.csv", replace
              restore
              display "CUBE_RECOVERY_PASS"

              Comment


              • #8
                This is great work, Choonjoo Lee. Thanks for sharing! I'll share #1 with my daughter in case she gets stuck solving a \(2 \times 2\) Rubik's cube.I've added the highlighted code to help her specify the initial color vector (sb2_start).

                Code:
                preserve
                capture quietly mata: sb2_goal()
                if _rc {
                mata:
                struct sb2_tree {
                    real matrix states
                    real colvector depth, parent, move
                    transmorphic scalar seen
                }
                struct sb2_result {
                    real scalar found, length, startnodes, goalnodes, checks
                    real rowvector path
                    real matrix layers
                }
                
                real matrix sb2_geometry()
                {
                    real matrix N,G
                    real rowvector n,a,b,c
                    real scalar f,i,j
                    N=(1,0,0\-1,0,0\0,0,1\0,0,-1\0,-1,0\0,1,0)
                    G=J(0,7,.)
                    for(f=1;f<=6;f++) {
                        n=N[f,]
                        if(abs(n[1])) { a=(0,1,0); b=(0,0,1); }
                        else if(abs(n[2])) { a=(1,0,0); b=(0,0,1); }
                        else { a=(1,0,0); b=(0,1,0); }
                        for(i=-1;i<=1;i=i+2) for(j=-1;j<=1;j=j+2) {
                            c=n+.5*i*a+.5*j*b
                            G=G\(c,n,f)
                       }
                   }
                    return(G)
                }
                
                real matrix sb2_moves(real matrix G)
                {
                    real matrix P,Q,H
                    real rowvector axis,ix,base
                    real scalar f,i,k
                    P=J(9,24,.)
                    // Row vectors multiply Q. Clockwise seen from outside R, U, F.
                    for(f=1;f<=3;f++) {
                        if(f==1) { Q=(1,0,0\0,0,-1\0,1,0); axis=(1,0,0); }
                        if(f==2) { Q=(0,-1,0\1,0,0\0,0,1); axis=(0,0,1); }
                        if(f==3) { Q=(0,0,-1\0,1,0\1,0,0); axis=(0,-1,0); }
                        H=G
                        for(i=1;i<=24;i++) if(G[i,1..3]*axis'>0) {
                            H[i,1..3]=G[i,1..3]*Q
                            H[i,4..6]=G[i,4..6]*Q
                       }
                        k=3*f-2
                        for(i=1;i<=24;i++) {
                            ix=selectindex(rowsum(abs(G[,1..6]:-H[i,1..6])):==0)'
                            assert(cols(ix)==1)
                            P[k,ix[1]]=i
                       }
                        base=P[k,]
                        P[k+1,]=base[base]
                        P[k+2,]=base[base[base]]
                        assert(all(sort(base',1):==(1::24)))
                        assert(all(base[base[base[base]]]:==(1..24)))
                   }
                    return(P)
                }
                
                real rowvector sb2_goal()
                {
                    return((1,1,1,1,2,2,2,2,3,3,3,3,4,4,4,4,5,5,5,5,6,6,6,6))
                }
                string scalar sb2_key(real rowvector s)
                {
                    return(invtokens(strofreal(s,"%1.0f"),""))
                }
                real scalar sb2_inverse(real scalar m)
                {
                    return(3*floor((m-1)/3)+4-mod(m-1,3)-1)
                }
                string scalar sb2_names(real rowvector path)
                {
                    string rowvector names
                    names=("R","R2","R'","U","U2","U'","F","F2","F'")
                    if(cols(path)==0) return("(none)")
                    return(invtokens(names[path]))
                }
                real rowvector sb2_parse(string scalar text)
                {
                    string rowvector words,names
                    real rowvector path,ix
                    real scalar i
                    words=tokens(text)
                    names=("R","R2","R'","U","U2","U'","F","F2","F'")
                    path=J(1,0,.)
                    for(i=1;i<=cols(words);i++) {
                        ix=selectindex(names:==words[i])
                        assert(cols(ix)==1)
                        path=path,ix[1]
                   }
                    return(path)
                }
                real rowvector sb2_apply(real rowvector s,real rowvector path,real matrix P)
                {
                    real scalar i
                    real rowvector out
                    out=s
                    for(i=1;i<=cols(path);i++) out=out[P[path[i],]]
                    return(out)
                }
                
                struct sb2_tree scalar sb2_bfs(real rowvector root,real matrix P,real scalar radius)
                {
                    struct sb2_tree scalar T
                    real scalar cap,head,n,d,m
                    real rowvector s
                    string scalar key
                    // Upper bound on all 9-branch paths. This demo limits radius to 4.
                    cap=(9^(radius+1)-1)/8
                    T.states=J(cap,24,.)
                    T.depth=J(cap,1,.)
                    T.parent=J(cap,1,0)
                    T.move=J(cap,1,0)
                    T.seen=asarray_create("string",1)
                    T.states[1,]=root; T.depth[1]=0
                    asarray(T.seen,sb2_key(root),1)
                    head=1; n=1
                    while(head<=n) {
                        d=T.depth[head]
                        if(d<radius) for(m=1;m<=9;m++) {
                            s=T.states[head,P[m,]]
                            key=sb2_key(s)
                            if(!asarray_contains(T.seen,key)) {
                                n++
                                T.states[n,]=s; T.depth[n]=d+1
                                T.parent[n]=head; T.move[n]=m
                                asarray(T.seen,key,n)
                           }
                       }
                        head++
                   }
                    T.states=T.states[1..n,]
                    T.depth=T.depth[1..n]; T.parent=T.parent[1..n]; T.move=T.move[1..n]
                    return(T)
                }
                real rowvector sb2_trace(struct sb2_tree scalar T, real scalar i)
                {
                    real rowvector path
                    path=J(1,0,.)
                    while(T.parent[i]>0) {
                        path=T.move[i],path
                        i=T.parent[i]
                   }
                    return(path)
                }
                
                struct sb2_result scalar sb2_solve(real rowvector start,real matrix P,real scalar radius)
                {
                    struct sb2_tree scalar A,B
                    struct sb2_result scalar R
                    real scalar i,j,d,best,ai,bi,c
                    real rowvector back
                    string scalar key
                    assert(cols(start)==24 & rows(start)==1)
                    assert(all(start:==floor(start)) & min(start)>=1 & max(start)<=6)
                    for(c=1;c<=6;c++) assert(sum(start:==c)==4)
                    assert(radius>=0 & radius<=4 & radius==floor(radius))
                    // Only start, legal moves and depth limit enter. No scramble history.
                    A=sb2_bfs(start,P,radius)
                    B=sb2_bfs(sb2_goal(),P,radius)
                    R.startnodes=rows(A.states); R.goalnodes=rows(B.states)
                    R.checks=rows(A.states); R.found=0; R.length=.
                    R.path=J(1,0,.); R.layers=J(radius+1,3,0)
                    for(d=0;d<=radius;d++) R.layers[d+1,]=(d,sum(A.depth:==d),sum(B.depth:==d))
                    best=.; ai=0; bi=0
                    // Check ALL meetings: first meeting need not give the shortest combined path.
                    for(i=1;i<=rows(A.states);i++) {
                        key=sb2_key(A.states[i,])
                        if(asarray_contains(B.seen,key)) {
                            j=asarray(B.seen,key)
                            d=A.depth[i]+B.depth[j]
                            if(d<best) { best=d; ai=i; bi=j; }
                       }
                   }
                    if(ai==0) return(R)
                    R.path=sb2_trace(A,ai)
                    back=sb2_trace(B,bi)
                    for(i=cols(back);i>=1;i--) R.path=R.path,sb2_inverse(back[i])
                    R.found=1; R.length=cols(R.path)
                    assert(R.length==best)
                    assert(all(sb2_apply(start,R.path,P):==sb2_goal()))
                    return(R)
                }
                
                // Static orthographic rendering. At quarter-turn endpoints the cube is convex;
                // outward-facing stickers are sufficient, without a general z-buffer renderer.
                void sb2_render(real rowvector state,real matrix G)
                {
                    real matrix out,pts
                    real rowvector right,up,cam,n,a,b,c
                    real scalar view,i,poly,shift
                    out=J(0,5,.); poly=0
                    for(view=1;view<=2;view++) {
                        right=(.8,.6,0); up=(-.36,.48,.8); cam=(.48,-.64,.6)
                        if(view==2) { right=-right; cam=-cam; }
                        shift=(view==1 ? 0 : 3.8)
                        for(i=1;i<=24;i++) {
                            n=G[i,4..6]
                            if(n*cam'<=0) continue
                            if(abs(n[1])) { a=(0,1,0); b=(0,0,1); }
                            else if(abs(n[2])) { a=(1,0,0); b=(0,0,1); }
                            else { a=(1,0,0); b=(0,1,0); }
                            // First sticker of each face also draws its dark backing square.
                            if(mod(i-1,4)==0) {
                                pts=(n-a-b\n+a-b\n+a+b\n-a+b\n-a-b)
                                poly++
                                out=out\(J(5,1,poly),pts*right':+shift,pts*up',J(5,1,0),J(5,1,view))
                           }
                            c=G[i,1..3]
                            pts=(c-.46*a-.46*b\c+.46*a-.46*b\c+.46*a+.46*b\c-.46*a+.46*b\c-.46*a-.46*b)
                            poly++
                            out=out\(J(5,1,poly),pts*right':+shift,pts*up',J(5,1,state[i]),J(5,1,view))
                       }
                   }
                    st_addobs(rows(out))
                    st_addvar("double",("polygon","px","py","color","view"))
                    st_store(.,.,out)
                    // Coordinates and labels for sticker positions 1–24
                st_addvar("double",("labelx","labely","slot"))
                
                poly=0
                for(view=1;view<=2;view++) {
                    right=(.8,.6,0)
                    up=(-.36,.48,.8)
                    cam=(.48,-.64,.6)
                
                    if(view==2) {
                        right=-right
                        cam=-cam
                    }
                
                    shift=(view==1 ? 0 : 3.8)
                
                    for(i=1;i<=24;i++) {
                        n=G[i,4..6]
                        if(n*cam'<=0) continue
                
                        c=G[i,1..3]
                        poly++
                        st_store(poly,("labelx","labely","slot"),
                            (c*right'+shift,c*up',i))
                    }
                }
                    
                    
                }
                end
                }
                clear
                mata:
                sb2_G=sb2_geometry()
                sb2_P=sb2_moves(sb2_G)
                // Current colors only: no scramble history or solution is supplied.
                // Slots: R(1-4), L(5-8), U(9-12), D(13-16), F(17-20), B(21-24).
                // Colors: 1 red, 2 orange, 3 white, 4 yellow, 5 green, 6 blue.
                sb2_start=(2,6,5,3,5,6,2,1,3,3,1,6,2,4,3,4,4,2,5,4,6,5,1,1)
                sb2_answer=sb2_solve(sb2_start,sb2_P,3)
                if(!sb2_answer.found) {
                    printf("No solution within 6 face turns; this does not prove impossibility.\n")
                    _error(3498)
                }
                printf("INPUT STATE: %s\n",sb2_key(sb2_start))
                printf("FOUND BY SEARCH: %s\n",sb2_names(sb2_answer.path))
                printf("Length: %g face turns; states: %g start + %g goal\n",sb2_answer.length,sb2_answer.startnodes,sb2_answer.goalnodes)
                st_numscalar("SB2_length",sb2_answer.length)
                st_numscalar("SB2_startnodes",sb2_answer.startnodes)
                st_numscalar("SB2_goalnodes",sb2_answer.goalnodes)
                st_matrix("SB2_layers",sb2_answer.layers)
                sb2_states=sb2_start
                for(sb2_i=1;sb2_i<=sb2_answer.length;sb2_i++) sb2_states=sb2_states\sb2_apply(sb2_states[rows(sb2_states),],sb2_answer.path[sb2_i],sb2_P)
                assert(all(sb2_states[rows(sb2_states),]:==sb2_goal()))
                end
                local panels
                forvalues step=0/`=scalar(SB2_length)' {
                clear
                mata: sb2_render(sb2_states[`step'+1,],sb2_G)
                local applied "Start"
                if `step'>0 {
                mata: st_local("applied",sb2_names(sb2_answer.path[`step']))
                }
                local next "SOLVED: all 24 facelets checked"
                if `step'<scalar(SB2_length) {
                mata: st_local("next","Next: "+sb2_names(sb2_answer.path[`step'+1]))
                }
                local colors `""40 45 53" "204 57 57" "239 145 50" "244 244 235" "246 211 64" "49 151 113" "66 111 185""'
                quietly levelsof polygon, local(polys)
                local plots
                foreach p of local polys {
                quietly summarize color if polygon==`p', meanonly
                local ci=r(min)+1
                local rgb : word `ci' of `colors'
                local plots `"`plots' (area py px if polygon==`p', color("`rgb'") lcolor("40 45 53") lwidth(vthin))"'
                local plots `"`plots' (scatter labely labelx, msymbol(none) mlabel(slot) mlabposition(0) mlabcolor(black) mlabsize(large))"'
                }
                twoway `plots', aspectratio(.5070422535) xscale(range(-1.65 5.45) noextend off) yscale(range(-1.6 2.0) noextend off) xlabel(none) ylabel(none) legend(off) text(1.75 0 "Front / right / up", size(small) color(gs6)) text(1.75 3.8 "Back / left / down", size(small) color(gs6)) title("Step `step': `applied'", size(medium)) subtitle("`next'", size(small)) graphregion(color(white)) plotregion(color(white)) xsize(9) ysize(5.8) name(sb2_step`step', replace)
                local panels `panels' sb2_step`step'
                }
                twoway scatteri 0 0, msymbol(none) xscale(range(0 10) off) yscale(range(0 10) off) xlabel(none) ylabel(none) legend(off) text(8.7 5 "Input: 24 colors, no move history", size(medium)) text(6.9 5 "Start side: `=scalar(SB2_startnodes)' distinct states", size(medium) color("40 104 150")) text(5.8 5 "Goal side: `=scalar(SB2_goalnodes)' distinct states", size(medium) color("40 104 150")) text(4.0 5 "Depth:       0     1     2      3", size(medium)) text(3.0 5 "Start: `=SB2_layers[1,2]' / `=SB2_layers[2,2]' / `=SB2_layers[3,2]' / `=SB2_layers[4,2]'", size(medium)) text(1.1 5 "Shortest path found: `=scalar(SB2_length)' face turns", size(medium)) title("The search behind the pictures", size(medium)) graphregion(color(white)) plotregion(color(white)) xsize(9) ysize(5.8) name(sb2_summary, replace)
                graph combine `panels' sb2_summary, cols(2) imargin(tiny) graphregion(color(white)) title("A 2x2 cube solved by search in Mata", size(medium)) note("Mata finds the moves; Stata draws the verified states. Both fixed views together show all six faces." "R/right, U/up, F/front. Prime: counterclockwise viewed from outside. 2: a half turn.", size(vsmall)) xsize(12) ysize(16) name(sb2_overview, replace)
                graph export "cube_search_2x2.png", width(1800) replace
                restore
                display "CUBE_SANDBOX_PASS"
                Click image for larger version

Name:	sb2_step0.png
Views:	1
Size:	153.7 KB
ID:	1787935


                Comment

                Working...
                X