--  An example of a TerraME model
--  Name: lucc_alocacao_discreto2

--  This is a simple diffusive deforestation allocation model,
--  which takes into consideration some spatial determinants to compute suitability for change
--  Steps:
--  (1) Potential for deforestation is calculated based on the number of deforested neighbours,
--      modified according to the cell attributes
--  (2) Then cells are ordered according to potential 
--  (3) Cells with the highest potential are changed

--  (c) Ana Paula Aguiar, Gilberto Camara, 2007


 -- Initialize the weights of evidence
 


function ave (tab) 
	n = table.getn (tab);
	sum = 0
	for i = 1, n do 
		sum = sum + tab[i];
	end
	return sum/n;
end

function std (tab, ave) 
	n = table.getn (tab);
	sqdiff = 0;
	for i = 1, n do 
		sqdiff = sqdiff + (tab[i] - dist.ave)^2;
	end
	return math.sqrt (sqdiff/n);
end

function median (tab) 
	n = table.getn (tab);
	table.sort (tab);
	i = math.floor (n/2);
	if ( math.mod(n, 2) == 1 ) then -- odd
		return tab[i+1];
	else -- even
		return (tab[i] + tab[i+1])/2;
	end
end
	  

function NormalParams ( tab )
	dist = {}
	dist.ave = ave (tab); 
	dist.std = std (tab, ave);
	return dist;
end

function FuzzyParams ( tab )  -- fuzzy(x) = 1/(1 +alpha*(x - beta)^2) 
	fuz = {}
	fuz.alpha = 1/(median(tab)^2);
	fuz.beta = 0;
	return fuz;
end

function normal (x, dist)
	return 1/(dist.std*math.sqrt(2*math.pi))*math.exp (- ((x - dist.ave)^2)/(dist.std^2))
end

function fuzzy (x, fuz)
	return 1/(1 +fuz.alpha*(x - fuz.beta)^2);
end

function calc_probs (cs1, cs2)
	n = table.getn (cs1.cells);
	dist_vic_change = {};
	dist_vic = {};
	dist_pri_change =  {};
	dist_pri = {};
	dist_sec_change = {};
	dist_sec = {};
	
	for i, cell in pairs (cs2.cells) do
		table.insert (dist_vic,  cs2.cells[i].dist_vicinal ); 
		table.insert (dist_pri,  cs2.cells[i].dist_ramal_princ);
		table.insert (dist_sec,  cs2.cells[i].dist_ramal_sec);
	end
	-- Step 1: Compute the potential for change
	-- The potential of change for each cell is the number of deforested neighbors
	for i = 1,n  do
		if ((cs1.cells[i].past.cover_=="forest") 
	   		and (cs2.cells[i].past.cover_ == "deforested")) then
       		-- assess potential for change according to distance  to roads
		    table.insert (dist_vic_change, cs2.cells[i].dist_vicinal ); 
			table.insert (dist_pri_change, cs2.cells[i].dist_ramal_princ);
			table.insert (dist_sec_change, cs2.cells[i].dist_ramal_sec);
		end
	end
	prob = {}
	prob.norm_vic  = NormalParams (dist_vic);
	prob.norm_pri = NormalParams (dist_pri);
	prob.norm_sec = NormalParams (dist_sec);
	
	prob.fuz_vic  = FuzzyParams (dist_vic_change);
	prob.fuz_pri = FuzzyParams (dist_pri_change);
	prob.fuz_sec = FuzzyParams (dist_sec_change); 
	
	return prob;
end

function potential (cs1, prob)
	for i, cell in pairs (cs1.cells) do
	   cell.potential = 0;
	   if (cell.past.cover_=="forest") then
		    -- modify potential for change according to distance  to roads
		     cell.potential = cell.potential + 
		     	math.log (fuzzy(cell.dist_vicinal,  prob.fuz_vic )/
		     	          normal(cell.dist_vicinal, prob.norm_vic));
		     	          
		     cell.potential = cell.potential + 
		     	math.log  (fuzzy (cell.dist_ramal_princ,  prob.fuz_pri )/
		     	           normal(cell.dist_ramal_princ,  prob.norm_pri));
		     
		     cell.potential = cell.potential + 
		     	math.log  (fuzzy (cell.dist_ramal_sec,  prob.fuz_sec )/
		     	           normal(cell.dist_ramal_sec,  prob.norm_sec));
		     	 
		end
	end
	it = SpatialIterator {
	   cs1,
	   function(cell) return cell.past.cover_== "forest"; end,
	   function (c1,c2) return c1.potential > c2.potential; end
	}
	return it;
end

--- MAIN PROGRAM -- 
-- CONSTANTS (MODEL PARAMETERS)
CELL_AREA = 0.25; -- 500 x 500 meters or 0.25 km2
DEMAND= 200; -- 100 km2
YEARS = 5;
FINAL_TIME = 5;

-- GLOBAL VARIABLES
cs85 = CellularSpace{
	dbType = "ADO",
	host = "localhost",
	database = "c:\\TerraME\\Database\\Rondonia.mdb",
	user = "",
	password = "",
	layer = "celulas_500_Anari_1985",
	theme = "celulas_500_Anari_1985",
	select = {"cover_" }
	}
	
-- RULES
cs85:load(); -- deveria ter o tempo dos atributos (tratar variáveis dinâmicas) 
cs85:synchronize();

-- GLOBAL VARIABLES
cs95 = CellularSpace{
	dbType = "ADO",
	host = "localhost",
	database = "c:\\TerraME\\Database\\Rondonia.mdb",
	user = "",
	password = "",
	layer = "celulas_500_Anari_1995",
	theme = "celulas_500_Anari_1995",
	select = {"cover_","soil_fert", "dist_urban", "protected", "dist_rodovia_BR", "dist_ramal_princ", "dist_ramal_sec", "dist_vicinal"}
	}
	
-- RULES
cs95:load(); -- deveria ter o tempo dos atributos (tratar variáveis dinâmicas) 
cs95:synchronize(); 	 

-- Step 1 - Find Weights of Evidence

prob = calc_probs (cs85, cs95);

-- Step 2 - Implement the Transitions

for time = 1, FINAL_TIME, 1 do
	demandaAtual = 20;
	cs95:synchronize();
	print("t= "..time);
	
	-- Step 2.1 - Calculate the potential for change and return an iterator
	it = potential (cs95, prob)
	
	-- Step 2.2: allocate changes to most suitable cells	 
	for i, cell in pairs( it.cells ) do
		demandaAtual = demandaAtual - CELL_AREA;
		if (demandaAtual < 0) then
			break;
		end;
		cell.cover_ = "deforested";
	end
    cs95:synchronize();
    
    if (time == FINAL_TIME) then
       cs95:save( 2000, "bayes_1995_", {"cover_"});
    end
end -- time step