crossBridge(T)

function [schedule,timings,total,iterations] = crossBridge(T)


% [schedule,timings,total,iterations] = crossBridge(T)
% 
% Dynamic Programming solver for The bridge crossing problem
% 
% The bridge-crossing problem is a mathematical puzzle where
% a group of N persons have to cross a bridge at night. 
% It is dark and they can only cross the bridge if they carry a lamp. 
% Only one lamp is available and at most two persons can cross at the same time. 
% It is not possible to cross from a side if the lamp is not on that side.
% The time of crossing is the time of the slowest person crossing.
% 
% 
% In this simple exercise, the bridge crossing problem is modelled as
% undiscounted Dynamic Programming problem with termination state.
% I know very little about this problem, and I made this script just for fun
% I guess other approaches exist which are more computationally efficient than this.
% 
% Anyway, this one works and can it be used to solve problems up to 10~15
% persons in reasonable time.
% 
% The input T is a row vector, one per person. T(i) is the time taken by the i-th
% person to cross.
% 
% the i-th row of 'schedule' matrix contains the index of the persons
% making the i-th crossing. If there is a zero then only one person is crossing
% The i-th element of 'timings' vector is the time taken for the i-th crossing,
% and 'total' is the total time taken.
% 'iterations' is the number of DP iterations to solve the problem
% 
% example: 
% 
% T = [ 1 2 5 10];
% [schedule,timings,total,iterations] = crossRiver(T)
%
% schedule =
%      1     2
%      1     0
%      3     4
%      2     0
%      1     2
% timings =
%      2
%      1
%     10
%      2
%      2
% total =
%     17
% iterations =
%     19



% number of persons
N = length(T);

pow2 = uint32( 2.^[0:N] );

% number of states
nX = 2^(N+1) - 2;

% reset the Bellman function
v0 = ones(nX,1);
v1 = zeros(nX,1);

iterations = 0;
i = uint32(0);
x = uint32(0);
while ~isequal(v1 , v0)

    iterations = iterations + 1;

    v0 = v1;
    
    fprintf('\n iteration %i',iterations);

    % explore all possible states except for the last one whose cost is zero
    for i = 1:nX-1

        % determine the i-th state
        x = bitand(i * uint32(ones(1,N+1)),pow2) ~= 0;
        x(1) = ~x(1);

        % x(1) is the position of the flag: 1 means Left and 0 means Right
        % x(i), i>1 is the position of the i-1 person: 0 means Left and 1 means Right

        % S is the set of person on the side with the flag
        if x(1) S = find(x) - 1;
        else S = find(~x) - 1;
        end
        % this is just to eliminate the zeros
        S = S(find(S));

        % determine the possible ways to cross Ax
        % determine the corresponding timinings t
        if length(S) > 1
            I = nchoosek(S,2)';
            Ax = [ I [S ; zeros(size(S))] ];
            t = [ max( T( I ) ) T(S)];
        else
            Ax = [S ; 0];
            t = T(S);
        end

        Q = [];
        
        for j = 1:size(Ax,2)

            % determine the next state if the j-th action is taken
            y = x;
            y(Ax(1,j) + 1) = ~y(Ax(1,j) + 1);
            if Ax(2,j) y(Ax(2,j) + 1) = ~y(Ax(2,j) + 1); end

            % determine the position of y in v
            yIndex = sum(pow2(y ~= 0));
            Q(j) = t(j) + v0( yIndex );
        end

        aBest = find(Q == min(Q));
        aBest = aBest(1);

        v1(i) = Q(aBest);

    end

end

% convergence achieved

schedule = [];
timings = [];

% obtain the result
i = 1;
while i < nX


    x = bitand(i * uint32(ones(1,N+1)),pow2) ~= 0;
    x(1) = ~x(1);

    % x(1) is the position of the flag: 1 means Left and 0 means Right
    % x(i), i>1 is the position of the i-1 person: 0 means Left and 1 means Right

    % S is the set of person on the side with the flag
    if x(1) S = find(x) - 1;
    else S = find(~x) - 1;
    end
    % this is just to eliminate the zeros
    S = S(find(S));

    % determine the possible ways to cross Ax
    % determine the corresponding timinings t
    if length(S) > 1
        I = nchoosek(S,2)';
        Ax = [ I [S ; zeros(size(S))] ];
        t = [ max( T( I ) ) T(S)];
    else
        Ax = [S ; 0];
        t = T(S);
    end

    Q = [];
    yIndex = [];
    for j = 1:size(Ax,2)

        % determine the next state if the j-th action is taken
        y = x;
        % the flag goes to the other side
        y(Ax(1,j) + 1) = ~y(Ax(1,j) + 1);
        if Ax(2,j) y(Ax(2,j) + 1) = ~y(Ax(2,j) + 1); end

        % determine the position of y in v
        yIndex(j) = sum(pow2(y ~= 0));

        Q(j) = t(j) + v0( yIndex(j) );
    end

    aBest = find(Q == min(Q));
    aBest = aBest(1);

    schedule(end+1 , :) = Ax(:,aBest)';
    timings(end+1 , 1) = t(aBest);
    i = yIndex(aBest);


end

total = sum(timings);