
--
-- Copyright (C) 2021  <fastrgv@gmail.com>
--
-- This program is free software: you can redistribute it and/or modify
-- it under the terms of the GNU General Public License as published by
-- the Free Software Foundation, either version 3 of the License, or
-- (at your option) any later version.
--
-- This program is distributed in the hope that it will be useful,
-- but WITHOUT ANY WARRANTY; without even the implied warranty of
-- MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-- GNU General Public License for more details.
--
-- You may read the full text of the GNU General Public License
-- at <http://www.gnu.org/licenses/>.
--




-- This package provides an interface to the munkres
-- Hungarian algorithm for usage as the heuristic
-- measure in an A* algorithm to solve sokoban puzzles.

-- Herein, the integer type "ushort", defined as 
-- mod 2**16 = 0..65535, is used to fill a cost matrix 
-- that defines the estimated number of moves to slide a 
-- particular box onto a particular goal.  The value 
-- usmx=65535 is used as "infinity" to indicate an 
-- impossible association between a particular box with 
-- a particular goal position,

with munkres;
with utils;

package body hungarian is





-- This is a more specialized application of the method I
-- use in mainstream::utils::dpbox but instead of discarding
-- the actual manhattan costs, here I save them in the
-- hungarian cost matrix. Also this sets individual costs 
-- for each box one-at-a-time so we can use munkres to
-- optimize pairings of boxes to goals.
--
-- Note that this initialization assumes boxes move on 
-- an empty maze, only using initial boxes & goals to 
-- determine all box-valid locations.
--
-- Note also that "bvalid" is NOT a minimal set of locations
-- where boxes reside without deadlock; It is a convenience
-- to define the domain of the hungarian algorithm where
-- manhattan distances define minimal cost of traversal.
procedure inithun(
	ivf,iff: vftype;
	ee:vustype;
	--pvalid: in booltype;
	bvalid: in out booltype;
	numrows, numcols: ushort
	) is --expensive; called once @ readPuzzle
	k: ushort:=0;
	gr,gc,kost, js,ii, up,dn,lf,rt: ushort;
	use text_io;
begin
	nrows:=numrows; --define local copies
	ncols:=numcols; --of dimensions

	hvf:=ivf; --define local copies
	hff:=iff; --of layout arrays

	for r in rowrng range 2..nrows-1 loop
	for c in colrng range 2..ncols-1 loop
		ii:=indx(r,c);
		up:=indx(r-1,c);
		dn:=indx(r+1,c);
		lf:=indx(r,c-1);
		rt:=indx(r,c+1);
		js:=ee(ii);
		if hff(ii)=2 then --goal for puller <==> box for pusher
			k:=k+1;
			hgoalrow(k):=r;
			hgoalcol(k):=c;
		end if;

		if 
			bvalid(ii) and --cell is valid location
			hff(ii)/=1 and --cell is not a wall OR is puller goal (hff=2)
			( 	(hff(dn)/=1 and hff(up)/=1) or --spaces above/below
				(hff(rt)/=1 and hff(lf)/=1) or --spaces left/right
				hvf(ii)=1  --box pos for puller [==pusher goal]
				or hff(ii)=2
			)
		then
			getboxcost(r,c); --defines hunmatrix(r,c)(*,*)
		end if;
		-- a) cannot be in a corner;
		-- b) if on a vertical wall, must have room above & below
		-- c) if on horizontal wall, must have room left & right
		-- d) valid if a box is there already
		--    (goal for pusher <=> box for puller)

	end loop;
	end loop;
	--myassert(k=nboxes,20202,"inithun:20");
	ub := k;
	nb := integer(k);


	nubvalid := (others=>false);
	for br in rowrng range 2..nrows-1 loop
	for bc in colrng range 2..ncols-1 loop
	for gx in 1..ub loop
		gr:=hgoalrow(gx);
		gc:=hgoalcol(gx);

		-- note hunmatrix NOT symmetric wrt boxes/goals
		kost := hunmatrix(br,bc)(gr,gc);

		if kost<usmx/2 then --this cell valid for some puller-goal...
			nubvalid( indx(br,bc) ) := true;
		   --a box @ (br,bc) can be pulled to goal @ (gr,gc)
		end if;

	end loop;
	end loop;
	end loop;

if false then
--debug
put("........inithun............"); new_line;
for br in 2..nrows-1 loop
for bc in 2..ncols-1 loop
ii:=indx(br,bc);
if nubvalid(ii) /= bvalid(ii) then
put("differ @ "&utrim(br)&","&utrim(bc));
	if bvalid(ii) then 
		put(", olbvalid=True");
	else
		put(", olbvalid=False");
	end if;
new_line;
end if;
end loop;
end loop;
end if; -- !different, but perfect results!


	bvalid:=nubvalid;

end inithun;








procedure setboxes(ivf: vftype) is --cheap
	k: ushort := 0;
	ii: ushort;
begin
	hvf:=ivf;

	--define current box layout numbering sequence
	for r in rowrng range 2..nrows-1 loop
	for c in colrng range 2..ncols-1 loop
		ii:=indx(r,c);
		if hvf(ii)=1 then
			k:=k+1;
			hboxrow(k):=r;
			hboxcol(k):=c;
		end if;
	end loop;
	end loop;

end setboxes;



function getBoxNum(br: rowrng; bc: colrng) return integer is
	bnum: integer := 0;
begin
	for k in boxrng loop
		if hboxrow(k)=br and hboxcol(k)=bc then
			bnum:=integer(k);
		end if;
	end loop;
	return bnum;
end getBoxNum;

procedure dpbox(
	r0,c0 : ushort; --puller.pos
	boxcost: in out boxcostype
	); --forward


--called for every B-valid interior point (r0,c0)
procedure getboxcost(r0: rowrng; c0: colrng) is
	boxcost: boxcostype := infinitecost;
begin
	--DP versus walls
	boxcost(r0,c0):=0;

	--now do DP assuming empty maze:
	dpbox(r0,c0,boxcost);
	-- the matrix boxcost represents the costs of 
	-- a) pushing a box from any position onto a goal @ (r0,c0)
	-- OR
	-- b) pulling a box @ (r0,c0) to any possible destination

	hunmatrix(r0,c0):=boxcost;

end getboxcost;











-- NOTE:  I have now specialized this to box-moves;
--        i.e. I allow extra cell for "puller"...
--			 This heuristic must Underestimate true cost.
-- define box domain using relaxation (flood-fill)
-- This defines pull-feasibility.
procedure dpbox(
	r0,c0 : ushort; --assumed current box position
	boxcost: in out boxcostype --manhattan distances to anywhere
	) is

	use text_io;

	ip: constant ushort := indx(r0,c0);
	ndelta : integer;
	cost,irc,ino,iso,iea,iwe: ushort;


--prep to test pull-feasibility
procedure initdpbox is
	use text_io;
	use myintio;
	ic: ushort;
begin
	hfff := (others=>1);
	for row in 1..nrows loop
	for col in 1..ncols loop
		ic:=indx(row,col);

		hfff(ic):=0;
		if hff(ic)=1  then hfff(ic):=1;  end if;

		hbestcost(ic):=usinf; -- usinf=usmx-1
		hviano(ic):=false;
		hviaso(ic):=false;
		hviaea(ic):=false;
		hviawe(ic):=false;

	end loop;
	end loop;

	for row in 2..nrows-1 loop
	for col in 2..ncols-1 loop
		ic:=indx(row,col);

		if 
			hfff(ic)=0                  -- no wall @ cell itself
			and hfff(indx(row-1,col))=0 -- no wall @ pred north
			and hfff(indx(row+1,col))=0 -- room for puller south
		then 
			hviano(ic):=true;
		end if;

		if 
			hfff(ic)=0 
			and hfff(indx(row+1,col))=0 -- no wall @ pred south
			and hfff(indx(row-1,col))=0 -- room for puller north
		then 
			hviaso(ic):=true; 
		end if;

		if 
			hfff(ic)=0                  -- no wall @ cell itself
			and hfff(indx(row,col+1))=0 -- no wall @ pred east
			and hfff(indx(row,col-1))=0 -- room for puller west
		then 
			hviaea(ic):=true; 
		end if;

		if 
			hfff(ic)=0 
			and hfff(indx(row,col-1))=0 -- no wall @ pred west
			and hfff(indx(row,col+1))=0 -- room for puller east
		then 
			hviawe(ic):=true; 
		end if;

	end loop;
	end loop;

end initdpbox;


begin -- dbbox

	initdpbox;
	hbestcost(ip):=0;
	ndelta:=5;



	-- we assume that any reachable position has a
	-- manhattan distance bounded by 254...

	while ndelta>0 loop
		ndelta:=0;

		for row in 2..nrows-1 loop --downward

			for col in 2..ncols-1 loop --rightward
				irc:=indx(row,col);
				ino:=indx(row-1,col);
				iso:=indx(row+1,col);
				iea:=indx(row,col+1);
				iwe:=indx(row,col-1);
				if hviano(irc) and hbestcost(irc)>hbestcost(ino)+1 then
					hbestcost(irc):=hbestcost(ino)+1;
					ndelta:=ndelta+1;
				end if;
				if hviaso(irc) and hbestcost(irc)>hbestcost(iso)+1 then
					hbestcost(irc):=hbestcost(iso)+1;
					ndelta:=ndelta+1;
				end if;
				if hviaea(irc) and hbestcost(irc)>hbestcost(iea)+1 then
					hbestcost(irc):=hbestcost(iea)+1;
					ndelta:=ndelta+1;
				end if;
				if hviawe(irc) and hbestcost(irc)>hbestcost(iwe)+1 then
					hbestcost(irc):=hbestcost(iwe)+1;
					ndelta:=ndelta+1;
				end if;
			end loop; --rightward

		end loop; --downward


		for row in reverse 2..nrows-1 loop --upward

			for col in reverse 2..ncols-1 loop --leftward
			irc:=indx(row,col);
				ino:=indx(row-1,col);
				iso:=indx(row+1,col);
				iea:=indx(row,col+1);
				iwe:=indx(row,col-1);
				if hviano(irc) and hbestcost(irc)>hbestcost(ino)+1 then
					hbestcost(irc):=hbestcost(ino)+1;
					ndelta:=ndelta+1;
				end if;
				if hviaso(irc) and hbestcost(irc)>hbestcost(iso)+1 then
					hbestcost(irc):=hbestcost(iso)+1;
					ndelta:=ndelta+1;
				end if;
				if hviaea(irc) and hbestcost(irc)>hbestcost(iea)+1 then
					hbestcost(irc):=hbestcost(iea)+1;
					ndelta:=ndelta+1;
				end if;
				if hviawe(irc) and hbestcost(irc)>hbestcost(iwe)+1 then
					hbestcost(irc):=hbestcost(iwe)+1;
					ndelta:=ndelta+1;
				end if;
			end loop; --leftward

		end loop; --upward


	end loop; --while ndelta


	for row in 1..nrows loop
		for col in 1..ncols loop
			irc:=indx(row,col);
			cost:=ushort( hbestcost(irc) );
			boxcost(row,col):=cost;
		end loop;
	end loop;


end dpbox;





function ndx(r,c,d: ushort) return integer is
begin
	return integer((r-1)*d+c);
end ndx;







--this is the core of the Hungarian algorithm...
function evalpairings(
	ontgt: out integer; 
	ok: out boolean
	) return ushort is
	use text_io;
	ib,j,i,totcost : integer := 0;
	bx,gx: ushort := 0;
	br,bc,gr,gc: ushort;
	lassn: munkres.iatype(1..nb);    --local assn array
	lcost: munkres.iatype(1..nb*nb); --local cost matrix
	kost: integer;
begin
	ontgt:=0;
	--prepare data for munkres:
	for bx in 1..ub loop

		br:=hboxrow(bx);
		bc:=hboxcol(bx);
	for gx in 1..ub loop
		i:=ndx(bx,gx,ub); --linearized 2D array
		gr:=hgoalrow(gx);
		gc:=hgoalcol(gx);

		-- note hunmatrix NOT symmetric wrt boxes/goals
		kost := integer( hunmatrix(br,bc)(gr,gc) );
		--cost of PULLing a box @ (br,bc) to goal @ (gr,gc)

		lcost(i):= kost; --local [properly-sized]
		gcost(i):= kost; --global [oversized]

	end loop; --gx
	end loop; --bx

	munkres.hungarian(lcost,lassn,Ok); --HungarianAlgo

	-- interpret resulting data;  
	-- use "usmx" as "infinity"=>impossible
	if Ok then

		for i in 1..nb loop
			gassn(i):=lassn(i); --Xfer local output to global ASSN array
			bx:=ushort(i); 
			gx:=ushort( lassn(i) ); --bogus
			j := ndx(bx,gx,ub);
			totcost:=totcost + lcost(j);
			if lcost(j)=0 then ontgt:=ontgt+1; end if;
		end loop;
		if totcost>integer(usmx) then totcost:=integer(usmx); end if;

	else
		totcost:=integer(usmx);
	end if;

	return ushort(totcost);
	
end evalpairings;







------- below here are Utility routines ----------

function getcost( bx, gx : boxrng ) return ushort is
begin
	return ushort( gcost( ndx(bx,gx,ub) ) );
end getcost;


function getcost(  
	r0: rowrng; c0: colrng;
	r1: rowrng; c1: colrng
) return ushort is
begin
	return hunmatrix(r0,c0)(r1,c1);
end getcost;





function getmatch(bx: boxrng) return boxrng is
	gx : boxrng;
begin
	gx := boxrng( gassn( integer(bx) ) );
	return gx;
end getmatch;

procedure getBoxRc(bx: boxrng; r: out rowrng; c: out colrng) is
begin
	r:=hboxrow(bx);
	c:=hboxcol(bx);
end getboxrc;

procedure getGoalRc(gx: boxrng; r: out rowrng; c: out colrng) is
begin
	r:=hgoalrow(gx);
	c:=hgoalcol(gx);
end getgoalrc;






end hungarian; --package
