%% Sunpack.m
% Soma_osf.m unpacks itself.
    fin = fopen("Soma_osf.m",'r');
    mkdir("Soma");
    L = fgetl(fin);
    while ~isequal(L,'%% EOF')
        F = ['Soma/' L(4:end)];
        disp(F)
        fout = fopen(F,'w');
        L = fgetl(fin);
        while length(L) < 2 || ~isequal(L(1:2),'%%')
            fprintf(fout,'%s\n',L);
            L = fgetl(fin);
        end
        fclose(fout);
    end
    fclose(fin);
%% Soma_Solve.m
function Soma_Solve(idx)
    if nargin == 0
        somafigure
        % McKeeman, MATLAB Central File Exchange #26346, 2010.
        obj = somasol();
        sols = obj.solutions();
        set(gcf,'userdata',sols)
    else
        sols = get(gcf,'userdata');
        idxstr = get(gca,'userdata');
        idxstr.String = int2str(idx);
        somashow(colorvals(sols{idx}),idx); 
    end
end
%% Spack.m
function Spack
    fout = fopen('../Soma_osf.m','w');
    d = dir;
    d = d(2:end);
    d = string({d.name})';
    d(d=="Sunpack.m") = [];
    d(1) = "Sunpack.m";
    for F = d'
        disp(F)
        fprintf(fout,'%%%% %s\n',F);
        fin = fopen(F,'r');
        L = fgetl(fin);
        while ~isequal(L,-1)
            fprintf(fout,'%s\n',L);
            L = fgetl(fin);
        end
        fclose(fin);
    end
    fprintf(fout,'%%%% %s\n','EOF');
    fclose(fout);
end
%% colorvals.m
function cvals = colorvals(sol)
    % somashow-readable solution
    cvals = zeros(3,3,3);
    perm = [ 3 2 4 5 6 7 1];
    for iv = 1:7
        bits = sol(iv);
        for iu = 0:31
            x = mod(iu,3);
            y = mod(floor(iu/3), 3);
            z = mod(floor(iu/9), 3);
            if (bitand(bits, bitshift(uint32(1), iu)))
                cvals(x+1,y+1,z+1) = perm(iv);
            end
        end
    end
end
%% somafigure.m
function somafigure
    % Initialize figure window for Soma cube
    clf
    set(gcf,'name','Soma', ...
        'numbertitle','off', ...
        'color',0.75*[1 1 1])
    uicontrol('Style','pushbutton', ...
        'String','+', ...
        'Units','normalized', ...
        'Position',[.56 .93 .08 .04], ...
        'Fontsize',12, ...
        'Fontweight','bold', ...
        'Callback',@somacb)
    uicontrol('Style','pushbutton', ...
        'String','-', ...
        'Units','normalized', ...
        'Position',[.36 .93 .08 .04], ...
        'Fontweight','bold', ...
        'Fontsize',12, ...
        'Callback',@somacb)
    idxstr = uicontrol( ...
        'units','normalized', ...
        'position',[.46 .93 .08 .04], ...
        'style','pushbutton', ...
        'string',0, ...
        'Fontweight','bold', ...
        'Fontsize',9, ...
        'callback',@somacb);
    set(gca,'position',[0 0 1 1], ...
        'userdata',idxstr, ...
        'clipping','off')
    axis(14*[-1 1 -1 1 -1 1])
    axis('square','off','vis3d')
    rotate3d on
    view(3)
    shg

    function somacb(arg,~)
        idxstr = get(gca,'userdata');
        idx = str2double(idxstr.String);
        switch arg.String
            case '+', idx = mod(idx,240)+1;
            case '-', idx = mod(idx-2,240)+1;
            otherwise, idx = randi(240);
        end
        idxstr.String = num2str(idx);
        Soma_Solve(idx)
    end
end
%% somashow.m
function somashow(cvals,idx)
% somashow(cvals)  Display solution colorvals. 

    cla
    Q = Q0;
    F = [1 5 7 3
         3 7 8 4
         1 3 4 2
         2 4 8 6
         1 2 6 5
         5 6 8 7];
    C = get(0,'defaultaxescolororder');
    idxstr = get(gca,'userdata');

    if nargin < 2
        idx = str2double(idxstr.String);
        idx = idx + 1;
        idxstr.String = int2str(idx);
    else
        idxstr.String = int2str(idx);
    end

    % Cube
    for j = 1:27
        patch( ...
            Vertices = Q{j}, ...
            Faces = F, ...
            FaceColor = C(cvals(j),:))
    end

    % Pieces
    for k = 1:7
        switch k
            case 1, jk = [ 1  2     4];  % V
            case 2, jk = [ 1  2  3  4];  % L
            case 3, jk = [ 1  2  3  5];  % T
            case 4, jk = [ 1  2  5  6];  % Z
            case 5, jk = [ 1  2  5 14];  % S
            case 6, jk = [ 1  2 11 14];  % R
            case 7, jk = [ 1  2  4 10];  % Y
        end
        t = 2*pi*k/7;
        shift = 15*[cos(t) sin(t) .15];
        for j = 1:length(jk)
            patch( ...
                Vertices = .75*Q{jk(j)} + shift, ...
                Faces = F, ...
                FaceColor = C(k,:))
        end
    end
    drawnow
    
    function Q = Q0
        % From Qube, the Digital Rubik's Cube.
        v = 0.90*qzero;
        Q = cell(3,3,3);
        for z = [-2 0 2]
            for y = [-2 0 2]
                for x = [-2 0 2]
                    Q{x/2+2,y/2+2,z/2+2} = v + [x y z];
                end
            end 
        end
    end
    
    function q0 = qzero
        % Unit cubelet.
        q0 = [-1 -1 -1
              -1 -1  1
              -1  1 -1
              -1  1  1
               1 -1 -1
               1 -1  1
               1  1 -1 
               1  1  1];
    end
end
%% somasol.m
function obj = somasol()
% Bill McKeeman, MATLAB Central File Exchange #26346, 2010.
% CBM, insert line 364, June 12, 2022.
%
% SOMASOL   Compute all the unique solutions to the soma cube.
% Examples:
%  >> obj     = somasol();
%  >> sols    = obj.solutions();
%  >> solct   = numel(sols);
%  >> as_vec  = obj.vecsol(sols{1}(1));  % it will be the T
%  >> as_bits = obj.bitsol(sols{1}(1));  % T as bit vector
%  >> howto   = obj.layers(sols{1});     % how to build solution 1
%
%  bot mid top
%  TLL RVS RRS    The first solution
%  TTL RVV YSS
%  TZL YZZ YYZ
%
% The soma cube puzzle consists of 7 pieces (named V L T Z R S Y) which can
% be assembled into a 3x3x3 cube (consisting of 27 labeled cubelets).  See
% the soma demo in MATLAB(TM) for more details.
%
% Historical note:  This computation was originally done in V4 MATLAB.  The
% generated tables were used in the MATLAB soma demo but the table
% generating code itself was lost.  I have recreated the table computation
% here. In V4 on a mac it took several hours.  In R2009b on a PC it took 3
% seconds.
%
% The mapping of cubelet labels to bit positions is rectangular, the lower
% left corner being 0, increasing by 1 in the x direction, by 3 in the y
% direction, and 9 in the z direction.  Converting a label to base-3
% notation gives the xyz coordinates of the cubelet.
%{
    bottom  middle   top
     6 7 8  151617  242526  
  y  3 4 5  121314  212223
     0 1 2   91011  181920
       x          
%}
% somasol computes 240 solutions, each of which is a vector of 7 piece
% positions.  
%
% The computational strategy is to consider the cube as a set, represented
% by the leftmost 27 1-bits in a uint32 value.  A piece position is another
% uint32 value (containing only 3 or 4 1-bits).  Two piece positions collide
% when their intersection (bitand) is non-zero.  Any collection consisting
% of one each of the seven pieces, the union (bitor) of which fills the 27
% bits is a solution.
%
% The first task is computing all 688 piece positions.  The second task is
% classifying them into mutually compatible subsets (theory needed here).
% Ref: McKeeman, W.M., A Formal Model for Space-filling Puzzles, 
% Machine Intelligence vol.8, Elcock&Michie (ed.), Wiley (1977) pp 86-93
% Conway, Berlekamp and Guy, Winning Ways, vol. 2.  See also
% http://www.cs.dartmouth.edu/~mckeeman
%
% The third task is to remove equivalent solutions (rotation and
% reflection) before doing the computation (special knowledge of the soma
% cube needed here). Finally backtracking is applied to the remaining
% search space.

  % 3x3x3 rotation permutation vectors.
  r = uint8([0  3  6]);  r = [r+2 r+1 r];    r = [r r+9 r+18]; xy = r;
  r = uint8([6  7  8]);  r = [r   r+9 r+18]; r = [r r-3 r-6];  yz = r;  
  r = uint8([18 9  0]);  r = [r   r+3 r+6];  r = [r r+1 r+2];  zx = r;
  
  % vefc categorized pieces positions and vefc categories
  [V,vs] = Vsets();   if size(vs,1) ~= 6; error('bad V'); end
  [L,ls] = Lsets();   if size(ls,1) ~= 4; error('bad L'); end
  %[T,ts] = Tsets();  % use T piece to remove rotations
  T = {v2b(uint8([0 1 2 4]))}; ts = [2 1 1 0];
  [Z,zs] = Zsets();   if size(zs,1) ~= 2; error('bad Z'); end
  [R,rs] = Rsets();   if size(rs,1) ~= 4; error('bad R'); end
  [S,ss] = Ssets();   if size(ss,1) ~= 4; error('bad S'); end
  [Y,ys] = Ysets();   if size(ys,1) ~= 4; error('bad Y'); end
  
  % remove reflections using Y piece positions
  ed = uint8([0 3 6]); leftface = v2b([ed; ed+9; ed+18]);
  for ih = 1:numel(Y); Y{ih} = remove(Y{ih}, leftface); end
    
  % compute mutually consistent sets of piece positions
  compat = makeTrials(vs, ls, ts, zs, rs, ss, ys);
 
  % solution strategy:  
  %   place T in its fixed position first.
  %   place V last (it is easiest).
  sols = {};
  for iz=1:size(compat,1)
    idx = compat(iz,:);         % idx order V L T Z S R Y
    trial = {T{1}, L{idx(2)}, Z{idx(4)}, ...
             R{idx(5)}, S{idx(6)}, Y{idx(7)}, V{idx(1)}};
    sols = [sols, solve(trial)];
  end
  
  obj = public();
  return;                       % that's all folks
  
  %          --- nested function definitions ---
  
  % the public interface
  function o = public()
    o = struct;
    o.solutions = sols;
    o.vecsol    = @b2v;
    o.bitsol    = @b2str;
    o.layers    = @b2layers;
  end

  % all piece positions in terms of labels for V L R Z R S Y
  function res = Vvecs() 
    p = uint8([0 1 3]);    % init V
    p = [p; p+1]; p = [p; p+3]; p = [p; p+9; p+18];
    p = rXYZ(p);
    res = unique(sort(p,2), 'rows');
  end

  function res = Lvecs() 
    p = uint8([0 1 2 3]);  % init L
    p = [p; p+3]; p = [p; p+9; p+18];
    p = rXYZ(p);
    res = unique(sort(p,2), 'rows');
  end

  function res = Tvecs() 
    p = uint8([0 1 2 4]);  % init T
    p = [p; p+3]; p = [p; p+9; p+18];
    p = rXYZ(p);
    res = unique(sort(p,2), 'rows');
  end

  function res = Zvecs() 
    p = uint8([0 1 4 5]);  % init Z
    p = [p; p+3]; p = [p; p+9; p+18];
    p = rXYZ(p);
    res = unique(sort(p,2), 'rows');
  end

  function res = Rvecs() 
    p = uint8([0 1 3 12]); % init R
    p = [p; p+1]; p = [p; p+3]; p = [p; p+9];
    p = rXYZ(p);
    res = unique(sort(p,2), 'rows');
  end

  function res = Svecs() 
    p = uint8([0 1 3 10]); % init S
    p = [p; p+1]; p = [p; p+3]; p = [p; p+9];
    p = rXYZ(p);
    res = unique(sort(p,2), 'rows');
  end

  function res = Yvecs() 
    p = uint8([0 1 3 9]);  % init Y
    p = [p; p+1]; p = [p; p+3]; p = [p; p+9];
    p = rXYZ(p);
    res = unique(sort(p,2), 'rows');
  end

  % get all piece positions as bit vectors
  function bits = Vbits(); bits = bitRep(Vvecs()); end
  function bits = Lbits(); bits = bitRep(Lvecs()); end
  function bits = Tbits(); bits = bitRep(Tvecs()); end
  function bits = Zbits(); bits = bitRep(Zvecs()); end
  function bits = Rbits(); bits = bitRep(Rvecs()); end
  function bits = Sbits(); bits = bitRep(Svecs()); end
  function bits = Ybits(); bits = bitRep(Yvecs()); end

  % get signatures for each piece position
  function sigs = Vsigs(); sigs = makeSigs(Vbits()); end
  function sigs = Lsigs(); sigs = makeSigs(Lbits()); end
  function sigs = Tsigs(); sigs = makeSigs(Tbits()); end
  function sigs = Zsigs(); sigs = makeSigs(Zbits()); end
  function sigs = Rsigs(); sigs = makeSigs(Rbits()); end
  function sigs = Ssigs(); sigs = makeSigs(Sbits()); end
  function sigs = Ysigs(); sigs = makeSigs(Ybits()); end

  % separate piece positions by signature
  function [sets, s] = Vsets(); [sets, s] = makeSets(Vbits(), Vsigs()); end
  function [sets, s] = Lsets(); [sets, s] = makeSets(Lbits(), Lsigs()); end
  function [sets, s] = Tsets(); [sets, s] = makeSets(Tbits(), Tsigs()); end
  function [sets, s] = Zsets(); [sets, s] = makeSets(Zbits(), Zsigs()); end
  function [sets, s] = Rsets(); [sets, s] = makeSets(Rbits(), Rsigs()); end
  function [sets, s] = Ssets(); [sets, s] = makeSets(Sbits(), Ssigs()); end
  function [sets, s] = Ysets(); [sets, s] = makeSets(Ybits(), Ysigs()); end

  % make vector of signatures parallel to vector of piece positions
  function sigs = makeSigs(pp)
    vc = v2b(vCat()); ec = v2b(eCat()); 
    fc = v2b(fCat()); cc = v2b(cCat());
    sigs = zeros(size(pp,1), 4);
    for ig=1:size(pp,1)
      t(1) = bitsum(bitand(pp(ig,:), vc));
      t(2) = bitsum(bitand(pp(ig,:), ec));
      t(3) = bitsum(bitand(pp(ig,:), fc));
      t(4) = bitsum(bitand(pp(ig,:), cc));
      sigs(ig,:) = t;                % signature for this piece position
    end
  end

  % given piece position & signature, separate by signature
  function [sets, s] = makeSets(pp, sigs)
    s = unique(sigs, 'rows');
    sets = cell(size(s,1), 1);
    for jm=1:size(s,1); sets{jm} = uint32([]); end
    for im=1:size(pp,1)
      for jm=1:size(s,1)
        if all(s(jm,:)==sigs(im,:))
          sets{jm} = [sets{jm}; pp(im)];
        end
      end
    end
  end

  % solve integer linear equation
  function compat = makeTrials(v1,v2,v3,v4,v5,v6,v7)
    compat = zeros(0,7,'uint8');
    for i1=1:size(v1,1)
      for i2=1:size(v2,1)
        for i3=1:size(v3,1)
          for i4=1:size(v4,1)
           for i5=1:size(v5,1)
             for i6=1:size(v6,1)
               for i7=1:size(v7,1)
                 s = sum([v1(i1,:);v2(i2,:);v3(i3,:);v4(i4,:);
                          v5(i5,:);v6(i6,:);v7(i7,:)]);
                          % v  e f c
                 if all(s==[8 12 6 1])   % a solution
                   compat = [compat; [i1,i2,i3,i4,i5,i6,i7]];
                 end
               end
             end
           end
          end
        end
      end
    end
  end

  % rotationally symmetric catgories
  function res = vCat()     % vertex
    res = rXYZ(uint8(0));   % 0 2 6 8 18 20 24 26
  end

  function res = eCat()     % edge
    res = rXYZ(uint8(1));   % 1 3 5 7 9 11 15 17 19 21 23 25
  end

  function res = fCat()     % face
    res = rXYZ(uint8(4));   % 4 10 12 14 16 22
  end

  function res = cCat()     % center
    res = uint8(13);        % 13
  end

  % rotation actions
  function res = rXY(p)
    res = xy(p+1);      % 90 degrees X into Y
    if size(p,2)==1; res = res'; end
  end

  function res = rYZ(p)
    res = yz(p+1);      % 90 degrees Y into Z
    if size(p,2)==1; res = res'; end
  end

  function res = rZX(p)
    res = zx(p+1);      % 90 degrees Z into X
    if size(p,2)==1; res = res'; end
  end

  function res = rXYZ(p)  % all 24 rotations
    p = unique([p; rXY(p); rXY(rXY(p)); rXY(rXY(rXY(p)))], 'rows');
    p = unique([p; rYZ(p); rYZ(rYZ(p)); rYZ(rYZ(rYZ(p)))], 'rows');
    p = unique([p; rZX(p); rZX(rZX(p)); rZX(rZX(rZX(p)))], 'rows');
    res = p;
  end

  % vector to bits
  function bits = v2b(vec)
    bits = uint32(0);
    for ib=1:numel(vec)
      bits = bitor(bits, bitshift(uint32(1),vec(ib))); 
    end
  end

  % bits to vector
  function vec = b2v(bits)
    vec = zeros(1,0, 'uint8');
    for iv=0:31
      if (bitand(bits, bitshift(uint32(1), iv)))
        vec(end+1) = iv;
      end
    end
  end

  % for printing
  function str = b2str(bits)
    str = char(zeros(1,32));
    for iv=0:31
      if (bitand(bits, bitshift(uint32(1), iv)))
        str(iv+1) = '1';
      else
        str(iv+1) = '0';
      end
    end
  end 

  % human-readable solution
  function str = b2layers(sol)
    str = char(zeros(3,3));
    names = 'TLZSRYV';  % the order in the solve input
    for iv=1:numel(sol)
      bits = sol(iv);
      for iu=0:31
        x = mod(iu,3); 
        y = mod(floor(iu/3), 3);  
        z = mod(floor(iu/9), 3); 
        if (bitand(bits, bitshift(uint32(1), iu)))
          str(x+1,y+1,z+1) = names(iv);
        end
      end
    end
  end 

  % add up bits
  function n = bitsum(bits)
    n = 0;
    for i=0:31
      if bitand(bits, uint32(1)); n = n + 1; end
      bits = bitshift(bits, -1);
    end
  end

  % change vector sets into bit sets
  function bits = bitRep(matrix)
    bits = zeros(size(matrix, 1), 1, 'uint32');
    for i=1:size(matrix,1)
      bits(i) = v2b(matrix(i,:));
    end
  end

  % Q is a vector of piece positions
  % F is any part of the 3x3x3 cube.  
  % Discard all Q that intersect F.
  function Q = remove(Q, F)
    for ij = size(Q,1):-1:1   % each piece
      if bitand(Q(ij), F) ~= 0
        Q(ij) = [];           % discard pieces that hit forbid
      end
    end
  end

  % given signature-compatible inputs, get solutions
  function sols = solve(ppsets)
    sols = {};
    cube = v2b(0:26);
    recur(uint32(0), uint32([]), 1);   % nothing so far...
    
    function recur(sofar, partial, n)
      if n > numel(ppsets)             % cannot try more
        if sofar == cube               % found solution
          sols = [sols, partial];
          somashow(colorvals(partial))  % cbm, 2022
        end
        return;
      end
      pps = ppsets{n};
      for ir = 1:size(pps,1)
        pp = pps(ir);
        if bitand(sofar,pp) == 0       % no collisions
          %fprintf('%d %d %s\n', n, ir, b2str(sofar))
          recur(bitor(sofar,pp), [partial, pp], n+1);
        end
      end
    end
  end
end
%% EOF