Cube Puzzle solver

% solves this puzzle: https://i0.wp.com/famousartisan.com/wp-content/uploads/2017/10/puzzle-cube.jpg?fit=800%2C470&ssl=1

%% define puzzle pieces
pieces = repmat({false(3,3,3)} ,7,1) ;
pieces{1}(:,:,1) = [ 1 1 1 ; 1 0 0 ; 0 0 0] ;
pieces{2}(:,:,1) = [ 1 1 1 ; 0 1 0 ; 0 0 0] ;
pieces{3}(:,:,1) = [ 1 1 0 ; 0 1 1 ; 0 0 0] ;
pieces{4}(:,:,1) = [ 1 1 0 ; 1 0 0 ; 0 0 0] ;
pieces{5}(:,:,1) = [ 1 1 0 ; 1 0 0 ; 0 0 0] ;
pieces{5}(:,:,2) = [ 1 0 0 ; 0 0 0 ; 0 0 0] ;
pieces{6}(:,:,1) = [ 1 1 0 ; 1 0 0 ; 0 0 0] ;
pieces{6}(:,:,2) = [ 0 0 0 ; 1 0 0 ; 0 0 0] ;
pieces{7}(:,:,1) = [ 1 1 0 ; 1 0 0 ; 0 0 0] ;
pieces{7}(:,:,2) = [ 0 1 0 ; 0 0 0 ; 0 0 0] ;
%% make all possible permutations of transformations of pieces

%rotations and offsets
rotateA = 0: pi/2 : 3 * pi/2 ;
ofsetA = 0 : 2 ;

% precalculate all rotation matrices
rotateIA = 1 : numel(rotateA) ;
RX = [] ;
RY = [] ;
RZ = [] ;
for i = rotateIA
    RX(:,:,i) =  makehgtform('xrotate',rotateA(i))  ;
    RY(:,:,i) =  makehgtform('yrotate',rotateA(i))  ;
    RZ(:,:,i) =  makehgtform('zrotate',rotateA(i))  ;
end
RX = RX(1:3,1:3,:) ;
RY = RY(1:3,1:3,:) ;
RZ = RZ(1:3,1:3,:) ;

rotComb = combvec(rotateIA,rotateIA,rotateIA) ;
RM = zeros(3,3,size(rotComb,2))  ;
for ri =1 : size(rotComb,2)
    RM(:,:,ri) = RX(:,:,rotComb(1,ri)) * RY(:,:,rotComb(2,ri)) * RZ(:,:,rotComb(3,ri))  ;
end

%All combinations of rotation matrices and x,y,z offsets
[allComb] = combvec(ofsetA,ofsetA,ofsetA,1:size(RM,3)) ;

tPieces = cell(7,size(allComb,2)) ;

for ci = 1 : 7
    piece = pieces{ci}  ;

    [xi,yi,zi] = ind2sub([3,3,3],find(piece)) ;

    X = [xi,yi,zi] - 1  ;
    for oi = 1 : size(allComb,2)
        xo = allComb(1,oi) ;
        yo = allComb(2,oi) ;
        zo = allComb(3,oi) ;
        ri = allComb(4,oi) ;

        R = RM(:,:,ri) ;
        newX = X * R+ [xo yo zo] +1 ;

        %Check if in bounds
        if all(newX(:) > 0) && all(newX(:) < 4)
            tPieces{ci,oi} = accumarray(newX,1,[3 3 3]) ;
        end
    end
end

% because each solution (whole cube) can be rotated such that the first
% piece is in any rotation, we only put the first piece in the first
% rotation by discarding the rest.
tPieces(1,allComb(4,:)>1) = cell(size(tPieces(1,allComb(4,:)>1))) ;

%% remove duplicates and invalid solutions (empties) from list
goodPieces = cell(1,7);
for ci = 1 : 7
    thesePieces = tPieces(ci,~cellfun(@isempty,tPieces(ci,:))) ;

    uPieces = unique(cell2mat(cellfun(@(t) t(:)',thesePieces,'uni',0)'),'rows') ;
    goodPieces{ci} = uPieces ;
end

%% verify goodPieces:
%{
[x,y,z] = meshgrid(1:3,1:3,1:3) ;
for ci = 1 %1 : 7

    figure(ci*100)
    plotPlacer(16)
    for oi = 1 : size(goodPieces{ci},1)
        piece = goodPieces{ci}(oi,:) ;

        [xi,yi,zi] = ind2sub([3,3,3],find(piece)) ;
        plotPlacer

        scatter3(xi,yi,zi,'fill')
        hold on  ;
        plot3(xi,yi,zi,'-')
        scatter3(x(:),y(:),z(:))
        axis equal
    end
end
%}

%% now check all combinations of all options wether they overlap?
NaPiece = cellfun(@(g) size(g,1),goodPieces) ;

thisTry = zeros(7,27) ;
altTry = zeros(27,7) ;
indexA = ones(1,7) ;
ci =1 ;
si = 1;
indexS = [] ;
solutions = [] ;
goingDown = false ;
while ci > 0
    %          disp([ci indexA])  % show iterations in
    thisTry(ci,:) = goodPieces{ci}(indexA(ci),:) ;

    indexA(ci) = indexA(ci)+1 ;
    if any(sum(thisTry)>1)
        %if we tried all, reset this piece and go back to the last
        %while ensures you can do that twice in a row
        while indexA(ci) > NaPiece(ci)
            thisTry(ci,:) = 0 ;
            indexA(ci) = 1 ;
            ci = ci - 1 ;
            %when it rolls back all the way to before the first piece  all
            %combinations have been tested
            if ci == 0
                break
            end
        end
    else
        if ci == 7
            %if this try works and its the seventh, a sollution is found
            indexS(:,si) = indexA ; solutions(:,:,si) = thisTry ; si =si + 1; %#ok
            %spoof going to the previous piece, ci will increase again
            ci = ci - 1 ;

            %break ; %if you only want to see one solution, break here
        end
        %if this try works, continue with another piece
        ci = ci + 1 ;
    end
end
%correct index increase before solution storage
indexS(:,end) = indexS(:,end) -1 ;

%% Show solution (there are 480 :O )

[x,y,z] = meshgrid(1:3,1:3,1:3) ;
figure(555) ;
if size(solutions,3) > 0
    plotPlacer(8,10)
else
    plotPlacer(1) ;
end
for si = 1 : size(solutions,3)
    solution = solutions(:,:,si)

    plotPlacer ;
    for ci = 1 : 7
        piece = reshape(solution(ci,:),[3 3 3]) ;
        [xi,yi,zi] = ind2sub([3,3,3],find(piece)) ;
        scatter3(xi,yi,zi,'fill')
        hold on  ;
        %     plot3(xi,yi,zi,'-') %gets messy, might help
        axis equal
    end
end

%% PlotPlacer function i made when i started using matlab bc subplot is anoying
%enjoy the spagetti

function hout = plotPlacer(varargin)
%PLOTPLACER Interface to subplot.
%   Subplot is tedious. plotPlacer only takes arguments the first call, and
%   remembers those for placing all other subplots. Every call without
%   arguments spawns a new subplot, starting anew in the current figure.
%
%   --Input (first call only!)
%      plotPlacer(nPlotx,nPloty).  %plots in this rectangle
%      plotPlacer(nPlot). %plot on a rectangle in 1 figure
%      plotPlacer(...,'vert') draws subsequent axes vertically ('horz' also
%      works. Because of legacy and horizontal and vertical are confusing)

% September 2014

%Todo?
%add support for non regular placement
%Force overwrite new figures as option?
%add option to supply number of plots and number of figures

%now, with output handles :O
hout = [] ;
persistent nPlotx nPloty fn sp nsp next ;
%place a plot
if nargin == 0
    if(sp == nsp)
        fn = fn + 1 ;
        %call figure here for clf ;
        figure(fn) ;
        clf ;
    end
    %try catches when inplots == 0 and do nothing
    try
        %always call figre for more robust plot placement
        figure(fn);
        sp = next(sp,nPlotx,nPloty) ;
        hout = subplot(nPloty,nPlotx,sp) ;
    catch
    end
else
    %process inputs
    if ischar( varargin{end} )
        options = varargin{end} ;
        varargin = varargin(1:end-1) ;
    end
    %initialize
    inSizes = cellfun(@numel,varargin) ;

    %non-vararr input
    if all(inSizes==1) && numel(inSizes) < 3
        %add 'horz' to plot them vertically. Now vert also works
        if exist('options','var') && (strcmp(options,'horz') || strcmp(options,'vert') )
            next = @(n,nx,ny) ((n+nx)>nx*ny)*(-nx*ny+1) + n + nx + (n==(nx*ny))*(-nx) ;
            %O yes he just did
        else
            next =  @(n,nx,ny) mod(n, nx * ny)+1 ;
        end

        switch numel(varargin)  %~= nargin! bc above the options is removed
            case 1
                nPlotx = ceil(sqrt(varargin{1})) ;
                nPloty = ceil(varargin{1}/nPlotx) ;
            case 2
                nPlotx = varargin{2} ;
                nPloty = varargin{1} ;
        end
    else
        %Hiddem feature. To unify the input with that of parPicker (deprec)
        if exist('options','var') && ~strcmp(options,'horz')
            names= strsplit(options) ;
            varnames = arrayfun(@inputname,1:numel(varargin) ,'uni',0) ;

            intmp = cellfun(@(x) strcmp(names,x),varnames,'uni',0) ;
            intmp = cell2mat(reshape(intmp,[numel(intmp),1])) ;

            pioritySizes = arrayfun(@(x) inSizes(intmp(:,x)),1:size(intmp,2)) ;
        else
            pioritySizes = inSizes(find(inSizes>1,2,'first'));
        end

        next =  @(n,nx,ny) mod(n, nx * ny)+1 ;
        if numel(pioritySizes)>0, nPlotx = pioritySizes(1) ; else  nPlotx = 1 ; end
        if numel(pioritySizes)>1, nPloty = pioritySizes(2) ; else  nPloty = 1 ; end
    end
    fig = gcf ;
    fn = get(fig,'Number') ;
    if ischar(fn)
        %Old matlab compatibility
        fn = fig ;
    end
    fn = fn - 1 ;%-1 because we want a call to figure first time
    nsp = nPlotx * nPloty;
    sp = nsp ;
end
end