lasttry01

function [colors,score] = solver(A,GOAL)
colors=[];

[LABELS,VALUES,IDX,R,N,S]=connected4(A);
WALL=~VALUES;
VALUES_GOAL=A(GOAL);
LABELS_GOAL=LABELS(GOAL);
[VALUES_U,VALUES_I,VALUES_J]=unique(VALUES);
VALUES_N=numel(VALUES_U);
VALUES_E=sparse(1:N,VALUES_J,true,N,VALUES_N);
UNITS=accumarray(VALUES_J,ones(N,1),[VALUES_N,1]);
COST=(VALUES_GOAL-VALUES).*S;COST(WALL)=inf;
COST1=sum(A(:));
COST2=sum(COST(~WALL));
COST3=accumarray(VALUES_J,COST);

CONNEC=connectivity;
[DISTANCES_START,DISTANCES_GOAL,DISTANCES_IDX1,DISTANCES_IDX2,DISTANCES_IDX3,VALUES_THR]=distances;
if isinf(DISTANCES_START(LABELS(GOAL))), colors=[];score=max(A(:))*sum(A>0);return;end

Z={}; COLORS=Z;C=[];
NDO=0;MAXDO=min(3,sum(DISTANCES_IDX2<=VALUES_GOAL));
COST1_TYPES=[2,3,4,5,1];
DO_BASELINE=1;
for N1=1:size(DISTANCES_START,2), 
    if ~isinf(DISTANCES_START(LABELS(GOAL),N1)),
        if NDO<MAXDO,
%             [Z{end+1},COLORS{end+1},C(end+1)]=evaluate(DISTANCES_START(:,N1),LABELS(GOAL));
            for COST1_TYPE=COST1_TYPES,[Z{end+1},COLORS{end+1},C(end+1)]=evaluate_stepwise_cost1(LABELS(GOAL),VALUES_THR(:,N1),DISTANCES_GOAL(:,N1));end
            NDO=NDO+1;
        end
    end; 
end
if DO_BASELINE,[COLORS{end+1},C(end+1)]=qcombo3(A,GOAL);Z{end+1}=[];end
[score,idxC]=min(C);
colors=COLORS{idxC};


function [Z,colors,C]=evaluate(d,GOAL)
    % backwards
    n=d(GOAL);
    if ~n||isinf(n),Z=[];colors=[];C=inf;return;end
    idx=zeros(1,n+1);
    idx(1)=find(d==0,1);
    idx(n+1)=GOAL;
    for n1=n+1:-1:3,
        a=find(CONNEC(:,idx(n1)));
        %b=find(d(a)==n1-2,1);
        [nill,b]=min(d(a));
        idx(n1-1)=a(b);
    end
    % forward
    Z=sparse(idx(1),1,true,N,1);
    colors=zeros(n,1);
    for n1=1:n-1,
        colors(n1)=VALUES(idx(n1+1));
        Z=Z|((CONNEC*Z>0)&VALUES_E(:,VALUES_J(idx(n1+1))));
    end
    minC=-inf;n1=0;
    % reconsider end
    while(minC<0),
        V=((CONNEC*Z)>0 & ~Z);
        COLORS_C=accumarray(VALUES_J(V),COST(V),[VALUES_N,1])+VALUES_U;
        COLORS_C(VALUES_J(idx(n+1)))=0;
        [minC,idxC]=min(COLORS_C);
        if minC<0,
            n1=n1+1;
            [nill,idxC2]=max(COLORS_C<0);
            colors(n+n1-1)=VALUES_U(idxC2);
            Z=Z|((CONNEC*Z>0)&VALUES_E(:,idxC2));
        end
    end
    % end
    colors(n+n1)=VALUES(idx(n+1));
    Z=Z|((CONNEC*Z>0)&VALUES_E(:,VALUES_J(idx(n+1))));
    if ~isempty(Z),C=COST1+(VALUES(GOAL)-VALUES(Z))'*S(Z)+sum(colors);else,C=inf;end
end

function [Z,colors,C]=evaluate_stepwise_cost1(GOAL,do,d)
    if nargin>2,units=UNITS-accumarray(VALUES_J,isinf(d),[VALUES_N,1]);else,units=UNITS;end
    % forward
    Y=sparse(1,1,true,N,1);
    W=Y;Z=Y;
    colors=zeros(N,1);
    COLORS_P=zeros(VALUES_N,3);COLORS_P(:,3)=inf;
    COLORS_C=inf+zeros(VALUES_N,1);
    COLORS_I=sparse([],[],[],N,VALUES_N);
    if nargin>2,Z0=d(1);else,Z0=DISTANCES_GOAL(1);end
    Y0=COST2-(VALUES_GOAL-VALUES(Y))*S(Y);X0=COST1-VALUES(Y)*S(Y);n3=0;
    for n1=1:N,
        % update neighborhood
        if nargin<2,V=((CONNEC*Y)>0 & ~W);
        else,V=((CONNEC*Y)>0 & ~W & do);end
        v=find(V);
        for n2=1:numel(v),
            vvalue=VALUES_J(v(n2));
            COLORS_P(vvalue,1)=VALUES(v(n2));
            COLORS_P(vvalue,2)=COLORS_P(vvalue,2)+S(v(n2));
            COLORS_P(vvalue,3)=min(COLORS_P(vvalue,3),d(v(n2)));
            COLORS_I(v(n2),vvalue)=true;
            units(vvalue)=units(vvalue)-1;
        end
        
        idx=find(COLORS_P(:,1));
        minD=min(COLORS_P(idx,3));
        k=((VALUES_GOAL-COLORS_P(idx,1)).*COLORS_P(idx,2) +COLORS_P(idx,1) );
        if COST1_TYPE==1,                   COLORS_C(idx)=COLORS_P(idx,3).*k; 
        elseif COST1_TYPE==2,               COLORS_C(idx)=(COLORS_P(idx,3)>0).*( k+(minD>0)*units(idx).*COLORS_P(idx,1) ); 
        elseif COST1_TYPE==3,               COLORS_C(idx)=units(idx);
        elseif COST1_TYPE==4,               COLORS_C(idx)=COLORS_P(idx,3);
        elseif COST1_TYPE==5,               COLORS_C(idx)=-COLORS_P(idx,3)./(1e-10+units(idx).*COLORS_P(idx,1)); 
        end
        
        [nill,idxbest]=min(COLORS_C);
        if COLORS_I(GOAL,idxbest),
            minC=-inf;n3=0;
            % reconsider end
            while(minC<0),
                V=((CONNEC*Z)>0 & ~Z);
                COLORS_C=accumarray(VALUES_J(V),COST(V),[VALUES_N,1])+VALUES_U;
                COLORS_C(VALUES_J(GOAL))=0;
                [minC,idxC]=min(COLORS_C);
                if minC<0,
                    n3=n3+1;
                    [nill,idxC2]=max(COLORS_C<0);
                    colors(n1+n3-1)=VALUES_U(idxC2);
                    Z=Z|(V&VALUES_E(:,idxC2));
                end
            end
            % end
            colors(n1+n3)=VALUES_U(VALUES_J(GOAL));
            Z=Z|((CONNEC*Z>0)&VALUES_E(:,VALUES_J(GOAL)));
            break;
        end
        
        
        
        colors(n1)=VALUES_U(idxbest);
        Y=COLORS_I(:,idxbest);
        W=W|V;
        Z=Z|Y;
        Y0=Y0-(VALUES_GOAL-COLORS_P(idxbest,1))*COLORS_P(idxbest,2);
        X0=X0-COLORS_P(idxbest,1)*COLORS_P(idxbest,2);
        COLORS_P(idxbest,:)=[0,0,inf];
        COLORS_C(idxbest)=inf;
        COLORS_I(:,idxbest)=false;
        if Z(GOAL),break;end
    end
    colors=colors(1:n1+n3,:);
    if ~isempty(Z),C=COST1+(VALUES(GOAL)-VALUES(Z))'*S(Z)+sum(colors);else,C=inf;end
    
end



function [Y1,Y2,idx1,idx2,idx3,VALUES_THR]=distances
[nill,sort3]=sort(-COST3');
%sort3=1:VALUES_N;
if min(VALUES)>0,
    VALUES_THR=repmat(cat(2,true(N,1),repmat(VALUES_J,[1,VALUES_N])~=repmat(sort3,[N,1])),[1,2]); %1:VALUES_N
    k=VALUES_N+1;
    idx1=1+sort3;%2:VALUES_N+1;
    idx2=[0;VALUES_U(sort3)];
    idx3=[0,sort3];%1:VALUES_N;
else,
    VALUES_THR=repmat(repmat(VALUES_J,[1,VALUES_N])~=repmat(sort3,[N,1]),[1,2]);%1:VALUES_N
    k=VALUES_N;
    idx1=sort3;%1:VALUES_N;
    idx2=VALUES_U(sort3);
    idx3=sort3;%1:VALUES_N;
end
X=sparse([ones(1,k),N+zeros(1,k)],1:2*k,true(1,2*k),N,2*k);
Y=inf(N,2*k);
Y(X)=0;
for n=1:N,
    X=(CONNEC*X)>0;
    X=X&(Y>n)&VALUES_THR;
    if ~any(X(:)),break;end
    Y(X)=n;
end
Y1=Y(:,1:k);Y2=Y(:,k+1:end);
end

function [G] = connectivity
[m n] = size(LABELS);
V = LABELS(:,1:n-1) ~= LABELS(:,2:n);
H = LABELS(1:m-1,:) ~= LABELS(2:m,:);
G = sparse([reshape(LABELS([V false(m,1)]),[],1); reshape(LABELS([H; false(1,n)]),[],1)], [reshape(LABELS([false(m,1) V]),[],1); reshape(LABELS([false(1,n); H]),[],1)], 1, N, N);
G = G + G' + speye(N);
G(WALL,:)=false;G(:,WALL)=false;
end

function [B,C,p,r,nc,s] = connected4(A)
% [K,a,idx,r,nc,s]=connected4(A);
% finds 4-connected components of identical elements in matrix A
%   A [m,n]
%   K [m,n] components
%   a [1,nc] values of A for each component
%   idx [1,m*n] r [1,nc+1] listed component indexes
%   such that:
%    A(idx(r(k):r(k+1)-1)=a(k)
%    K(idx(r(k):r(k+1)-1)=k
%

[m n] = size(A);
N = m*n ;
K = reshape (1:N, m, n) ;
V = A(:,1:n-1) == A(:,2:n);
H = A(1:m-1,:) == A(2:m,:);
G = sparse([reshape(K([V false(m,1)]),[],1); reshape(K([H; false(1,n)]),[],1)], [reshape(K([false(m,1) V]),[],1); reshape(K([false(1,n); H]),[],1)], 1, N, N);
G = G + G' + speye(N);
[p q r] = dmperm(G);p=p';q=q';r=r';
nc = numel(r) - 1;
C = A(p(r(1:nc)));
B = ones(m, n);
for a = 2:nc
    B(p(r(a):r(a+1)-1)) = a;
end
s=diff(r);
end

end











function [c,s] = qcombo3(A,goal)
[cc vcc ncc] = group(A);
if (max(A(:))>40)
    %[c,s] = special(A,goal,cc,vcc,ncc);
    c=[];s=inf;
else
    c = lateorearly(A,goal);
    [c,s] = glean(c,A,cc,vcc,ncc,goal);
end;
end

function [colors,score] = glean(colors,A,cc,vcc,ncc,goal)
ecc = vcc-A(goal);
hcc = hist(cc(:),1:ncc);
cben = hcc.*ecc;  % benefit of each cc
[nrow,ncol] = size(A);
flood = cc==cc(1);
for c = colors(1:end-1)'
    xflood = ~flood&([false(1,ncol);flood(1:end-1,:)]|[false(nrow,1),flood(:,1:end-1)]|[flood(2:end,:);false(1,ncol)]|[flood(:,2:end),false(nrow,1)]);
    vcc = cc(xflood&(A==c));
    flood = flood|ismember(cc,vcc);
end;
while true
    xflood = ~flood&([false(1,ncol);flood(1:end-1,:)]|[false(nrow,1),flood(:,1:end-1)]|[flood(2:end,:);false(1,ncol)]|[flood(:,2:end),false(nrow,1)]);
    cval = nonzeros(unique(A(xflood)));
    ncval = numel(cval);
    scr = zeros(1,ncval);
    vcc = cell(1,ncval);
    for i = 1:ncval
        if (cval(i)==colors(end))
            scr(i) = -inf;
        else
            vcc{i} = cc(xflood&(A==cval(i)));
            scr(i) = sum(cben(vcc{i}))-cval(i);
        end;
    end;
    [sch,ich] = max(scr);
    if (sch < 0)
        break
    end;
    colors = [colors(1:end-1);cval(ich);colors(end)];
    flood = flood|ismember(cc,vcc{ich});
end;
A(flood) = colors(end);
score = sum(A(:))+sum(colors);
end

function [colors,score] = special(A,goal,cc,vcc,ncc)
excess = A-A(goal);
ecc = vcc-A(goal);
[colors,score] = solvetaboo(A,cc,ecc,ncc,goal,[],excess);
cval = unique(A);
dear = cval < A(goal);
for i = dear
    [icolors,iscore] = solvetaboo(A,cc,ecc,ncc,goal,i,excess);
    if (iscore < score)
        score = iscore;
        colors = icolors;
    end;
end;
end


function [colors,score] = solvetaboo(A,cc,ecc,ncc,goal,taboo,excess)
[nrow,ncol] = size(A);
trow = mod(goal-1,nrow)+1;
tcol = ceil(goal/nrow);
colors = zeros(1,ncc);
flood = cc==cc(1);
move = 0;
while ~flood(goal)
    move = move+1;
    xflood = ~flood&([false(1,ncol);flood(1:end-1,:)]|[false(nrow,1),flood(:,1:end-1)]|[flood(2:end,:);false(1,ncol)]|[flood(:,2:end),false(nrow,1)]);
    cval = nonzeros(unique(A(xflood)));
    ncval = numel(cval);
    scr = rand(1,ncval);
    scr(cval > 40) = 2;
    scr(ismember(cval,taboo)) = 1.5;
    [sch,ich] = min(scr);
    colors(move) = cval(ich);
    chvcc = cc(xflood&(A==cval(ich)));
    flood = flood|ismember(cc,chvcc);
end;
colors = [colors(1:move-1) cval(cval > 40) colors(move)];
Ascr = A;
Ascr(flood) = A(goal);
score = sum(Ascr(:))+sum(colors);
end

function c = lateorearly(A, t)
B = uint8(A);
nc = max(max(A));
[G D ng] = group(B);
C = comp_cost(B, G, t, ng);
T = connections(G, D, ng);
C = push_channels(T, C);
c = dijkstra(T, nc, C, 1, G(t));
end

% Based on Tim Davis' find_components
function [B C nc] = group(A)
[m n] = size(A);
N = m*n;
K = reshape(1:N, m, n);
V = A(:,1:n-1) == A(:,2:n);
H = A(1:m-1,:) == A(2:m,:);
G = sparse([reshape(K([V false(m,1)]),[],1); reshape(K([H; false(1,n)]),[],1)], [reshape(K([false(m,1) V]),[],1); reshape(K([false(1,n); H]),[],1)], 1, N, N);
G = G + G' + speye(N);
[p q r] = dmperm(G);
nc = numel(r) - 1;
C = A(p(r(1:nc)));
B = ones(m, n, 'uint16');
for a = 2:nc
    B(p(r(a):r(a+1)-1)) = uint16(a);
end
end

% Similar idea to Dijkstra's shortest path (but not optimal here)
function [path loss Tcurr] = dijkstra(T, nC, cost, from, to)
nG = size(T, 1);
nT = ceil(nG * 1.5);
parent = ones(nG, nT, class(T));
k = zeros(nT, 1, 'uint16');
pointers = nT + ones(nG, 1, 'uint16');
pointers(from) = 1;
loss = zeros(nT+1, 1);
loss(:) = 1e9;
loss(1) = cost(from);
Tcurr = ones(nG, nT, class(T));
Tcurr(:,1) = T(:,from);
Tcurr(from,1) = 2;
sz = [nC+3 1];
noff = [0 0 0 1:nC]' * 2;
changed = false(nT, 1);
changed(1) = true;
b = 0;
while 1
    b = b + 1;
    if b > nT
        if ~any(changed)
            break;
        end
        b = 1;
    end
    if ~changed(b)
        continue;
    end
    changed(b) = false;
    if Tcurr(to,b) == 2
        continue;
    end
    ncost = accumarray(Tcurr(:,b), cost, sz);
    ncost(1:2) = 1e5;
    ncost = loss(b) + (ncost + noff) + 1;
    I = ncost(Tcurr(:,b)) < loss(pointers);
    if ~any(I)
        continue;
    end
    L = find(I)';
    [cl I] = sort(Tcurr(L,b));
    L = L(I);
    D = [0; find(diff(cl) ~= 0); numel(cl)];
    cl = cl(D(2:end));
    pointers(L) = numel(loss);
    I = true(nT+1, 1);
    I(pointers) = false;
    n = numel(cl);
    C = find(I(b+1:end-1), n) + b;
    n = n - numel(C);
    if n
        C = [C; find(I(1:b), n)];
    end
    ncost = ncost(cl);
    loss(C) = ncost;
    changed(C) = true;
    n = k(b) + 1;
    k(C) = n;
    K = Tcurr(:,b) == 1;
    for d = 1:numel(C)
        c = C(d);
        Tcurr(:,c) = Tcurr(:,b);
        J = Tcurr(:,c) == cl(d);
        Tcurr(J,c) = 2;
        Tcurr(K,c) = max(T(K,J), [], 2);
        parent(1:n,c) = parent(1:n,b);
        parent(n,c) = cl(d);
        pointers(L(D(d)+1:D(d+1))) = c;
    end
end
to_ = pointers(to);
path = parent(1:k(to_),to_) - 3;
%loss = loss(to_);
end

function G = connections(A, B, ng)
C = uint32(A);
[h w] = size(A);
T = false(ng);
I = A(1:h-1,:) ~= A(2:h,:);
X = C([I; false(1,w)]);
Y = C([false(1,w); I]);
T(X+(Y-1)*ng) = true;
T(Y+(X-1)*ng) = true;
I = A(:,1:w-1) ~= A(:,2:w);
X = C([I false(h,1)]);
Y = C([false(h,1) I]);
T(X+(Y-1)*ng) = true;
T(Y+(X-1)*ng) = true;
[y x] = find(T);
B = B + 3;
B(B==3) = 1;
G = ones(ng, class(B));
G(T) = B(y);
end

function C = comp_cost(A, B, t, ng)
C = int32(A);
C = C(t) - C(:);
C(A==0) = 1e6;
C = int32(accumarray(B(:), C, [ng 1]));
end

% Heuristic to deal with large regions connected by thin channels
function C = push_channels(T, C)
M = T > 3;
Q = find(sum(M, 1) == 1);
while ~isempty(Q)
    q = Q(1);
    r = find(M(:,q));
    if isscalar(r)
        C(r) = C(r) + C(q);
        C(q) = 0;
        M(q,r) = false;
        Q(1) = r;
    else
        Q = Q(2:end);
    end
end
end