percolation.sd7

(********************************************************************)
(*                                                                  *)
(*  percolation.sd7 Given a composite systems comprised of randomly *)
(*    distributed insulating and metallic materials: what fraction  *)
(*    of the materials need to be metallic so that the composite    *)
(*    system is an electrical conductor? Given a porous landscape   *)
(*    with water on the surface (or oil below) under what conditions*)
(*    will the water be able to drain through to the bottom (or the *)
(*    oil to gush through to the surface)? Scientists have defined  *)
(*    an abstract process known as <i>percolation</i> to model this *)
(*                                                                  *)
(*  Purpose: implement percolation model for Algorithms Part I      *)
(*           by Robert Sedgewick and Kevin Wayne                    *)
(*                                                                  *)
(*  Copyright (C) 2014, 2015  Arkady Kuleshov                       *)
(*                                                                  *)
(*  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 2 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 should have received a copy of the GNU General Public       *)
(*  License along with this program; if not, write to the           *)
(*  Free Software Foundation, Inc., 51 Franklin Street,             *)
(*  Fifth Floor, Boston, MA  02110-1301, USA.                       *)
(*                                                                  *)
(********************************************************************)

$ include "seed7_05.s7i";
  include "time.s7i";
  include "duration.s7i";
  include "float.s7i";
  include "draw.s7i";
  include "stdfont9.s7i";
  include "pixmap_file.s7i";
  include "keybd.s7i";
  include "unionfnd.s7i";

const type: percolation is new interface;

const func boolean: open(inout percolation: model, in integer: i, in integer: j) is DYNAMIC;

const func integer: dimensions(in percolation: model) is DYNAMIC;

const func boolean: percolates(in percolation: model) is DYNAMIC;

const func boolean: isOpen(in percolation: model, in integer: i, in integer: j) is DYNAMIC;

const func boolean: isFull(in percolation: model, in integer: i, in integer: j) is DYNAMIC;

const type: virtualGrid is new struct
		var integer: N is 0;
		var integer: size is 0;
		const integer: virtualIndex is 0;
		var integer: percolationIndex is 0;
		var array boolean: status is 0 times FALSE;
		var unionfind: unionVirtualTop is unionfind.value;
		var unionfind: unionVirtualBottom is unionfind.value;
	end struct;

type_implements_interface(virtualGrid, percolation);

const proc: initVirtualGrid(inout virtualGrid: vgrid, in integer: N) is func
	local
		var integer: rowIndex is 0;
		var integer: rowIndex2 is 0;
	begin
		vgrid.N := N;
		vgrid.size := N * N + N;
		init_unionfind(vgrid.unionVirtualTop, vgrid.size);
		init_unionfind(vgrid.unionVirtualBottom, vgrid.size);
		vgrid.status := vgrid.size times FALSE;
        vgrid.virtualIndex := 1;
        vgrid.percolationIndex := 0;

        rowIndex := 2 * N;
        rowIndex2 := vgrid.size;
        while rowIndex > N do
            union(vgrid.unionVirtualTop, vgrid.virtualIndex, rowIndex);
            union(vgrid.unionVirtualBottom, vgrid.virtualIndex, rowIndex2);
            decr(rowIndex);
            decr(rowIndex2);
        end while;

	end func;

const proc: connect(inout virtualGrid: vgrid, in integer: idx1, in integer: idx2) is func
	begin
		union(vgrid.unionVirtualTop, idx1, idx2);
		union(vgrid.unionVirtualBottom, idx1, idx2);
	end func;

const func integer: coordsTo1D(in integer: N, in integer: i, in integer: j) is
	return i * N + j - 1;


(**
 * open site (row i, column j) if it is not already
 *
 * @param i 1-based row index
 * @param j 1-based column index
 *)
const func boolean: open(inout virtualGrid: vgrid, in integer: i, in integer: j) is func
  result
    var boolean: have_opened is FALSE;
	local
		var integer: index is 0;
		var integer: offset is 0;
		var integer: N is 0;
	begin
        #validateInput(i, j);
        N := vgrid.N;
        index := coordsTo1D(N, i, j);
        if not vgrid.status[index] then
	        # connect to upper site
	        offset := index - N;
	        if i > 1 and vgrid.status[offset] then
	            connect(vgrid, index, offset);
	        end if;
	        # connect to left site
	        offset := index - 1;
	        if j > 1 and vgrid.status[offset] then
	            connect(vgrid, index, offset);
	        end if;
	        # connect to right site
	        offset := index + 1;
	        if j < N and vgrid.status[offset] then
	            connect(vgrid, index, offset);
	        end if;
	        # connect to lower/bottom site
	        offset := index + N;
	        if i < N and vgrid.status[offset] then
	            connect(vgrid, index, offset);
	        end if;
	        vgrid.status[index] := TRUE;
	        if is_connected(vgrid.unionVirtualTop, index, vgrid.virtualIndex) and
			       is_connected(vgrid.unionVirtualBottom, index, vgrid.virtualIndex) then
	            vgrid.percolationIndex := max(vgrid.percolationIndex, index);
	        end if;
          have_opened := TRUE;
        end if;
	end func;

const func integer: dimensions(in virtualGrid: vgrid) is func
	result
		var integer: dim is 0;
	begin
		dim := vgrid.N;
	end func;

const func boolean: percolates(in virtualGrid: vgrid) is func
	result
		var boolean: res is FALSE;
	begin
		res := vgrid.percolationIndex > 0;
	end func;

const func boolean: isOpen(in virtualGrid: vgrid, in integer: i, in integer: j) is func
	result
		var boolean: res is FALSE;
	local
		var integer: index is 0;
	begin
		index := coordsTo1D(vgrid.N, i, j);
		res := vgrid.status[index];
	end func;

const func boolean: isFull(in virtualGrid: vgrid, in integer: i, in integer: j) is func
	result
		var boolean: res is FALSE;
	local
		var integer: index is 0;
	begin
		index := coordsTo1D(vgrid.N, i, j);
		if vgrid.status[index] = TRUE and vgrid.percolationIndex > 0
			and is_connected(vgrid.unionVirtualTop, vgrid.percolationIndex, index) then
			res := TRUE;
		end if;
	end func;

const integer: AREA_LENGTH is 480;
const integer: WINDOW_WIDTH is AREA_LENGTH;
const integer: TEXT_AREA_HEIGHT is 100;
const integer: WINDOW_HEIGHT is AREA_LENGTH + TEXT_AREA_HEIGHT;

const proc: delay (in float: secs) is func
	local
		var time: start_time is time.value;
		var integer: seconds is 0;
		var integer: micro_seconds is 0;
	begin
		start_time := time(NOW);
		seconds := trunc(secs);
		micro_seconds := round((secs - flt(seconds)) * 1000000.0);
		await(start_time + seconds . SECONDS + micro_seconds . MICRO_SECONDS);
	end func;

const proc: draw_update(in percolation: model) is func
	local
		var float: tile_size is 0.0;
    var integer: N is 0;
		var integer: x1 is 0;
		var integer: y1 is 0;
		var integer: x2 is 0;
		var integer: y2 is 0;
		var integer: row is 0;
		var integer: col is 0;
	begin
    N := dimensions(model);
		(* set X scale [0..N] and Y scale [0..N] *)
		tile_size := flt(AREA_LENGTH) / flt(N);
		(* loop row from 1 to N *)
		for row range 1 to N do
			(* loop col from 1 to N *)
			for col range 1 to N do
				x1 := trunc(tile_size * flt(col - 1)) + 1;
				x2 := trunc(tile_size * flt(col)) - 1;
				y1 := trunc(tile_size * flt(row - 1)) + 1;
				y2 := trunc(tile_size * flt(row)) - 1;
        (* if isFull draw light blue rectangle in this position *)
				if isFull(model, row, col) then
					rect(x1, y1, x2 - x1 + 1, y2 - y1 + 1, dark_blue);
				(* else draw white for isOpen rectangle *)
				elsif isOpen(model, row, col) then
					rect(x1, y1, x2 - x1 + 1, y2 - y1 + 1, white);
        else
          rect(x1, y1, x2 - x1 + 1, y2 - y1 + 1, dark_gray);
        end if;
			end for;
		end for;
		flushGraphic;
	end func;

(* generate random number between 1 and dimensions(model) ^ 2 *)
const proc: open_random(inout percolation: model, inout integer: counter) is func
	local
		var integer: dim is 0;
		var integer: col is 0;
		var integer: row is 0;
		var integer: randNum is 0;
	begin
		dim := dimensions(model);
		randNum := rand(0, dim * dim - 1);
		row := randNum div dim;
		col := randNum mod dim;
		if open(model, row + 1, col + 1) then
      incr(counter);
    end if;
	end func;


const func string: build_stats(in array float: data) is func
  result
    var string: stri is "";
  local
    var float: total is 0.0;
    var float: k is 0.0;
  begin
    if length(data) > 0 then
      for k range data do
        total +:= k;
      end for;
      stri &:= "" <& (total/flt(length(data))) digits 2 <& "% avg for "
          <& length(data) <& " round(s)";
    end if;
  end func;


const proc: main is func
	local
		var char: ch is ' ';
    var integer: num_open is 0;
		var integer: num_blocks is 10;
		const float: DELAY is 0.08; (* sec *)
    var array float: all_simulations is 0 times 0.0;
    var float: pct_open is 0.0;
		var text: info_sheet is STD_NULL;
		var boolean: cont_simulation is TRUE;
		var virtualGrid: grid is virtualGrid.value;
	begin
		screen(WINDOW_WIDTH, WINDOW_HEIGHT);
		clear(curr_win, white);
		color(white, black);
		KEYBOARD := GRAPH_KEYBOARD;
		info_sheet := openPixmapFontFile(curr_win, 3, AREA_LENGTH + 1);
		setFont(info_sheet, stdFont9);
    if length(argv(PROGRAM)) = 1 then
      num_blocks := integer(argv(PROGRAM)[1]);
    end if;
		initVirtualGrid(grid, num_blocks);

		while cont_simulation do
      (* draw black background *)
      rect(1, 1, AREA_LENGTH, AREA_LENGTH, black);
      color(info_sheet, black, white);
			clear(info_sheet);
			writeln(info_sheet, "Percolation demo");
			writeln(info_sheet, "Current simulation:");
			flush(info_sheet);
			repeat
				open_random(grid, num_open);
				draw_update(grid);
        setPos(info_sheet, 2, 12);
        pct_open := flt(num_open)*100.0/flt(num_blocks*num_blocks);
        write(info_sheet, "" <& pct_open digits 2 <& "% is open ");
        delay(DELAY);
			until percolates(grid);
      all_simulations &:= pct_open;
			setPos(info_sheet, 3, 1);
      writeln(info_sheet, "Statistics: " <& build_stats(all_simulations));
      write(info_sheet, "Continue (Y or N) ? ");
			repeat
				ch := upper(getc(KEYBOARD));
			until ch = 'Y' or ch = 'N';
			if ch = 'N' then
				cont_simulation := FALSE;
      else
        initVirtualGrid(grid, num_blocks);
        num_open := 0;
      end if;
			setPos(info_sheet, 4, 1);
			write(info_sheet, "        ");
		end while;
	end func;