--  An example of a TerraME model
--  This is a simple diffusive model of deforestation
--  Potential for deforestation is calculated as an
--  average of the cell's neighbours
--
-- (c) Gilberto Camara, Ana Paula Aguiar, 2007

-- CONSTANTS (MODEL PARAMETERS)


-- GLOBAL VARIABLES
cs = CellularSpace{
	dbType = "ADO",
	host = "localhost",
	database = "c:\\TerraME\\Database\\amazonia.mdb",
	user = "",
	password = "",
	layer = "celulas100",
	theme = "dinamica",
	select = {"defor"}
}
-- RULES

-- CONSTRAINTS 
-- Budget to protect (limit to government costs)
-- Budget to deforest (or area?)
CELL_AREA = 10000;
FINAL_TIME = 1;
DEMAND   = 30000;  -- (km2)
LIMIT      =   30;
LAND_PRICE = 100;   -- (money/km2)
BUDGET     = 500000*LAND_PRICE;  -- (money)
DEFOR_PRICE = 3 * LAND_PRICE; -- (money/km2)
PROT_PRICE = 10 * LAND_PRICE; -- (money/km2)
PROD_PRICE = 20 * LAND_PRICE; -- (money/km2)
WOOD_PRICE = LAND_PRICE;      -- (money/km2)
FINE_VALUE = 100*LAND_PRICE;  -- (money/km2)
PERC_FINE =  10;               

function fine_risk ()
	if math.random (0,100) > PERC_FINE then
	   return FINE_VALUE;
	else
	   return 0; 
	end
end
-- How much it costs to deforest?

function defor_cost_km2 ( perc )
	return LAND_PRICE + fine_risk() + DEFOR_PRICE*(1 - perc);  -- função de custo de desmatamento;
end
-- How much it pays to deforest?

function defor_gain_km2 ( perc )
	return PROD_PRICE*perc + WOOD_PRICE*(1 - perc);
end
-- How much it costs to protect?
function prot_cost_km2 ( perc )
    max = PROT_PRICE*perc; 
	if ( perc < 0.25 ) then return max; 
	else 
		return (max/0.75) * (1 - perc);
	end
end
-- How much it pays to protect? 
function prot_gain_km2 ( perc )
	return 0;
end

cs:load();
CreateMooreNeighbourhood(cs);


for time = 1, FINAL_TIME, 1 do
	cs:synchronize();
	print("t: "..time );
	-- Step 1 = Calcule the deforestation benefit and the protection cost -- in money units
	for i, cell in pairs( cs.cells ) do
         cell.def_payoff  = (defor_gain_km2 (cell.past.defor) - defor_cost_km2 (cell.past.defor))*CELL_AREA;
         cell.prot_cost   = (prot_cost_km2 (cell.past.defor))*CELL_AREA; 
         cell.potential   = 0;
         cell.change = 0;
    end
    
    -- Step 2 = Adjust the protention capacity according to cost and budget
    -- Step 2.1 - Calculate the total protection cost in money units
    
    total_prot_cost = 0;
    for i, cell in pairs( cs.cells ) do
         total_prot_cost = total_prot_cost + cell.prot_cost; 
    end
    print ("total_prot_cost "..total_prot_cost);
    print ("total_prot_budget "..BUDGET);
    -- Step 2.2 - Adjust the protection capacity -- in money units
    rate = 1;
    if ( total_prot_cost > BUDGET ) then  
    	rate = BUDGET/total_prot_cost;
    end
    for i, cell in pairs( cs.cells ) do
           cell.prot_capacity = rate*cell.prot_cost;  -- in money units
    end
    -- Step 3 - Calculate the potential for change

    -- Step 3.1: Calculate the total_potential for change (in money units)
	total_potential = 0;
	for i, cell in pairs( cs.cells ) do
		  -- Fight between deforester and government
		  cell.potential = cell.def_payoff - cell.prot_capacity;  -- in money units
		  if (cell.potential > 0) then 
           	 total_potential  = total_potential + cell.potential; -- in money units
          end
    end
    
    print ("total_potential "..total_potential);
    
    -- Step 3.3: Transform the total_potential for change from money units to land
    defor_km2 = total_potential/LAND_PRICE; -- money/(money/km2) -- km2
    print ("defor_km2  "..defor_km2);
    if (defor_km2 > DEMAND ) then
    	defor_km2 = DEMAND;
    end
    print ("defor  "..defor_km2);

   -- ajust the demand for each cell so that
   -- the maximum demand for change is 100%
   -- adjust the demand so that excess demand is
   -- allocated to the remaining cells 
   -- there is an error limit (30 km2 or 0.1%)

   total_demand = defor_km2; 
   it2 = SpatialIterator {
       cs,
       function(cell) return cell.potential > 0; end,
	   function (c1,c2) return c1.potential > c2.potential; end
	}
   for i, cell in pairs( it2.cells ) do
   		print ("cell.potential "..cell.potential);
   end
	
   while (total_demand > LIMIT ) do  
   	   print("total_demand: "..total_demand );
       for i, cell in pairs( it2.cells ) do
           if (cell.potential > 0 ) then
	           -- calculate the potential to change for each cell
	           prop_cell = cell.potential/total_potential;
	           -- calculate how much area will be deforested in the cell
	           newarea = prop_cell* total_demand;
	           -- calculate the new percentage of deforestation
	           cell.defor = cell.past.defor + newarea/CELL_AREA;
	               
	           if (cell.defor > 1) then
	               total_potential = total_potential - cell.potential;
	               cell.potential  = 0;
	               excess = (cell.defor - 1)*CELL_AREA;
	               cell.defor = 1;
	           else
	               excess = 0;
	           end
	           -- adjust the total demand
	           total_demand = total_demand - (newarea - excess);  -- in km2
           end
           cell.change = cell.defor - cell.past.defor;
        end 
       cs:synchronize(); 
    end
end
cs:save( FINAL_TIME, "defor_t_", {"defor", "change"} );