```size(4).  % Maximum content of a tube
empty(2). % Number of empty tubes

% Tubes are given as ID-contents, where the the first color is the top one.
% In the case of a color of 2 or more units, it is written multiple times.
tubes([1-[blue,orange,cyan,indigo],
2-[red,lime,yellow,purple],
3-[purple,wheat,wheat,grey],
4-[pink,orange,lime,green],
5-[grey,cyan,pink,blue],
6-[indigo,lime,wheat,yellow],
7-[indigo,blue,purple,green],
8-[indigo,lime,red,cyan],
9-[cyan,wheat,purple,red],
10-[red,green,pink,orange],
11-[blue,yellow,yellow,grey],
12-[pink,green,grey,orange],
13-[],
14-[]]).

% Another example, where there seems to be no solution with 2 empty tubes.
% tubes([1-[blue,indigo,orange,red],
%        2-[green,green,orange,yellow],
%        3-[green,wheat,grey,cyan],
%        4-[purple,pink,wheat,yellow],
%        5-[pink,red,indigo,orange],
%        6-[cyan,blue,indigo,lime],
%        7-[pink,purple,wheat,cyan],
%        8-[grey,blue,yellow,grey],
%        9-[pink,red,orange,yellow],
%        10-[lime,purple,purple,grey],
%        11-[cyan,red,lime,blue],
%        12-[wheat,lime,green,indigo],
%        13-[],
%        14-[]]).

%%%%%%%%%%%
% Program %
%%%%%%%%%%%

% top(Xs, Y, Ys, Zs)
%   Xs is the concatenation of Ys = [Y,Y,...] and Zs
top([], _, [], []) :- !.
top([Y|Xs], Y, [Y|C], Xs1) :- !, top(Xs, Y, C, Xs1).
top([X|Xs], Y, [], [X|Xs]) :- X \= Y.

% pour(Ts, I-J, Ts1)
%   Starting from Ts, Ts1 is the state after pouring from I into J
pour(Ts, I-J, [I-Ti1,J-Tj1|Ts2]) :-
select(I-Ti, Ts, Ts1),
Ti = [C0|_], top(Ti, C0, C, Ti1),
select(J-Tj, Ts1, Ts2),
( Tj = [] -> Ti1 = [_|_] ; Tj = [C0|_] ),
length(Tj, L), length(C, Lc),
size(N), L + Lc =< N,
append(C, Tj, Tj1).

one_color([]) :- !.
one_color([_]) :- !.
one_color([X,X|Xs]) :- one_color([X|Xs]).

% A tube is finished if it is full and of one color
finished(T) :- one_color(T), size(N), length(T, N).

%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Solver w/o optimization %
%%%%%%%%%%%%%%%%%%%%%%%%%%%

% solve(Tubes, [I-J|Acc], Steps) :-
%     member(J-T, Tubes), finished(T), !,
%     select(J-T, Tubes, Tubes1),
%     empty(E),
%     ( length(Tubes1, E) -> reverse([I-J|Acc], Steps)
%     ; solve(Tubes1, [I-J|Acc], Steps)
%     ).
% solve(Tubes, Acc, Steps) :-
%     pour(Tubes, Step, Tubes1),
%     solve(Tubes1, [Step|Acc], Steps).
%
% solve(Steps) :- tubes(Tubes), solve(Tubes, [], Steps).

%%%%%%%%%%%%%%%%%%%%
% Optimized solver %
%%%%%%%%%%%%%%%%%%%%

% Idea for optimization:
% Forbid starting tubes that could have been chosen earlier but were skipped

% forbidden(Xs, Y-Z, L)
%   L contains all indices in Xs before Y, except for Z
forbidden([Y-_|_], Y-_, []) :- !.
forbidden([Z-_|Xs], Y-Z, L) :- !, forbidden(Xs, Y-Z, L).
forbidden([X-_|Xs], Y-Z, [X|L]) :- X \= Y, X \= Z, forbidden(Xs, Y-Z, L).

solve(Tubes, Forbidden, [I-J|Acc], Steps) :-
member(J-T, Tubes), finished(T), !,
select(J-T, Tubes, Tubes1),
empty(E),
( length(Tubes1, E) -> reverse([I-J|Acc], Steps)
; solve(Tubes1, Forbidden, [I-J|Acc], Steps)
).
solve(Tubes, Forbidden, Acc, Steps) :-
pour(Tubes, I-J, Tubes1),
( member(I, Forbidden) -> fail
; forbidden(Tubes, I-J, Forbidden1)
),
solve(Tubes1, Forbidden1, [I-J|Acc], Steps).

solve(Steps) :- tubes(Tubes), solve(Tubes, [], [], Steps).
```