lm

function d = solver(map,n)

[x,y]=size(map);
d = zeros(x,y);
s=sum(map(:))/n;
ovidx=1;
dir=0;
ndir=0;
%Find single districts
[h,l]=find(map>=s);
idx=1;
lh=length(h);
if lh>0,
  for idx=1:lh
    d(h(idx),l(idx))=idx;
  end
  idx=idx+1;
end

%Find high density districts
[h,l]=ind2sub(size(map),find(map(find(d==0))>=s*.90));
lh=length(h);
if lh>0,
 for idxc=1:lh
  nn=1;
  hh=h(idxc);
  ll=l(idxc);
  d(hh,ll)=idx;
  pos=[];
  ss=sum(map(find(d==idx)));
  if (hh > 1) & d(hh-1,ll)==0    %up
      pos=[pos [hh-1 ll ss+map(hh-1,ll)-s]'];
  end
  if (ll > 1) & d(hh,ll-1)==0    %left
      pos=[pos [hh ll-1 ss+map(hh,ll-1)-s]'];      
  end
  if(ll < y) & d(hh,ll+1)==0     %right
      pos=[pos [hh ll+1 ss+map(hh,ll+1)-s]'];      
  end
  if (hh < x) & d(hh+1,ll)==0    %down
      pos=[pos [hh+1 ll ss+map(hh+1,ll)-s]'];      
  end
  if ~isempty(pos)
      [hh,ll]=min(abs(pos(3,:)));
      d(pos(1,ll),pos(2,ll))=idx;
  end  
  idx=idx+1;
 end

end

while(idx<=n)
  [a,b] =find(~d);
  j=1;
  if isempty(a)
    %Could not allocate all
    %For now just take away from other districts
    n1=n-idx+1;
    for idxc=1:min(n1,ovidx-1)    
      [a,b]=find(d==ov(idxc));
      d(a(1),b(1))=idx;
      idx=idx+1;
    end
    if idx<=n
     [a,b] =find(~d);
      if isempty(a)
        %Could not allocate all
        %For now just take away from other districts
        n1=n-idx+1;
        for idxc=1:n1
          d(idxc) = idx;
          idx=idx+1;
        end  
      end
     end
else    

 while(j~= 0)
   aa = a(j);
   bb = b(j);
   j = j-1;
   d(aa,bb) = idx;
   if dir~=1
   %move right
   if (bb < y) & d(aa,bb+1)==0
      j = j+1;
      a(j) = aa;
      b(j) = bb+1;
   end
   
   %move up
   if (aa > 1) & d(aa-1,bb)==0
      j = j+1;
      a(j) = aa-1;
      b(j) = bb;
   end
   end
   

   %move down
   if (aa < x) & d(aa+1,bb)==0
      j = j+1;
      a(j) = aa+1;
      b(j) = bb;
   end
   %move left
   if (bb > 1) & d(aa,bb-1)==0
      j = j+1;
      a(j) = aa;
      b(j) = bb-1;
      if bb==1
        ndir=1;
      end
   end
   
   if dir==1
   %move right
   if (bb < y) & d(aa,bb+1)==0
      j = j+1;
      a(j) = aa;
      b(j) = bb+1;
   end
   %move up
   if (aa > 1) & d(aa-1,bb)==0
      j = j+1;
      a(j) = aa-1;
      b(j) = bb;
      if aa==1
        ndir=0;
      end    
    end
   end
   dir=ndir;

   if idx<n
     if sum(map(find(d==idx)))>=s
       if sum(map(find(d==idx)))-s>map(aa,bb)/1.95
          d(aa,bb)=0;
       else
          ov(ovidx)=idx;
          ovidx=ovidx+1;
       end
       break
     end
   end
 end
end
if j==0
    ild=find(d==idx);
    if sum(map(ild))<.5*s
        d(ild)=-1;
    else
        idx=idx+1;
    end
else
idx=idx+1;
end
end

d(find(d==-1))=0;

[a,b] =find(~d);
if ~isempty(a)
    %Could not allocate all
    %For now just attach unresolved locations to nearest district
 for j=1:length(a)
   aa = a(j);
   bb = b(j);
      
   if (aa > 1) & d(aa-1,bb)~=0        %get up
      idx=d(aa-1,bb);
   elseif (bb > 1) & d(aa,bb-1)~=0    %get left
      idx=d(aa,bb-1);
   elseif(bb < y) & d(aa,bb+1)~=0     %get right
      idx=d(aa,bb+1);      
   elseif (aa < x) & d(aa+1,bb)~=0    %get down
      idx=d(aa+1,bb);
   end
   d(aa,bb)=idx;
 end
end