#!/usr/bin/env python
############################################################################
#
# MODULE:	   r.mcda.roughset
# AUTHOR:	   Gianluca Massei - Antonio Boggia
# PURPOSE:	  Generate a MCDA map from several criteria maps using Dominance Rough Set Approach - DRSA 
#					   (DOMLEM algorithm proposed by  (S. Greco, B. Matarazzo, R. Slowinski)	 
# COPYRIGHT:  c) 2010 Gianluca Massei, Antonio Boggia  and the GRASS 
#					   Development Team. This program is free software under the 
#					   GNU General PublicLicense (>=v2). Read the file COPYING 
#					   that comes with GRASS for details.
#
#############################################################################

#%Module
#% description: Generate a MCDA map from several criteria maps using Dominance Rough Set Approach.
#% keywords: raster, Dominance Rough Set Approach 
#% keywords: Multi Criteria Decision Analysis (MCDA)
#%End
#%option
#% key: criteria
#% type: string
#% multiple: yes
#% gisprompt: old,cell,raster
#% key_desc: name
#% description: Name of criteria raster maps 
#% required: yes
#%end
#%option
#% key: preferences
#% type: string
#% key_desc: character
#% description: gain,cost
#% required: yes
#%end
#%option
#% key: decision
#% type: string
#% gisprompt: old,cell,raster
#% key_desc: name
#% description: Name of decision raster map 
#% required: yes
#%end
#%option
#% key: outputMap
#% type: string
#% gisprompt: new,cell,raster
#% description: Output classified raster map
#% required: yes
#%end
#%option
#% key: outputTxt
#% type: string
#% gisprompt: new_file,file,output
#% key_desc: name
#% description: Name for output files (base for *.isf  and *.rls files)
#% answer:infosys
#% required: yes
#%end
#%flag
#% key: l
#% description: Do not remove single rules in vector format
#% answer:false
#%end
#%flag
#% key: n
#% description: Compute null value as zero
#% answer:true
#%end

import sys
import copy
import numpy as np
from time import time, ctime  
import grass.script as grass
import grass.script.array as garray


def BuildFileISF(attributes, preferences, decision, outputMap, outputTxt):
	outputTxt=outputTxt+".isf"
	outf = file(outputTxt,"w")
	outf.write("**ATTRIBUTES\n")
	for i in range(len(attributes)):
		outf.write("+ %s: (continuous)\n" % attributes[i])
	outf.write("+ %s: [" % decision)
	value=[]
	value=grass.read_command("r.describe", flags = "1n", map = decision)
	v=value.split()

	for i in range(len(v)-1):
		outf.write("%s, " % str(v[i]))
	outf.write("%s]\n" % str(v[len(v)-1]))
	outf.write("decision: %s\n" % decision)

	outf.write("\n**PREFERENCES\n")
	for i in range(len(attributes)):
		if(preferences[i]==""):
			preferences[i]="none"
		outf.write("%s: %s\n" % (attributes[i], preferences[i]))
	outf.write("%s: gain\n" % decision)
	
	if flags['n']:
		for i in range(len(attributes)):
			print "%s - convert null to 0" % str(attributes[i])
			grass.run_command("r.null", map=attributes[i], null=0)
			
	outf.write("\n**EXAMPLES\n")
	examples=[]
	MATRIX=[]

	for i in range(len(attributes)):
		grass.mapcalc("rast=if(isnull(${decision})==0,${attribute},null())",
						 rast="rast",
						 decision=decision,
						 attribute=attributes[i])
		tmp=grass.read_command("r.stats", flags = "1n", nv="?", input = "rast")
		example=tmp.split()
		examples.append(example)
	tmp=grass.read_command("r.stats", flags = "1n", nv="?", input = decision)
	example=tmp.split()
	examples.append(example)

	MATRIX=map(list,zip(*examples))
	MATRIX=[r for r in MATRIX if not '?' in r] #remove all rows with almost one "?"
	MATRIX=[list(i) for i in set(tuple(j) for j in MATRIX)] #remove duplicate example 
	
	for r in range(len(MATRIX)):
		for c in range(len(MATRIX[0])):
			outf.write("%s " % (MATRIX[r][c]))
#			outf.write("%s " % round(float(MATRIX[r][c]), 2))
		outf.write("\n")
	   
	outf.write("**END")
	outf.close()
	return outputTxt
	
	

def collect_attributes (data):
	"Collects the values of header files isf, puts them in an array of dictionaries"
	header=[]
	attribute=dict()
	j=0
	start=(data.index(['**ATTRIBUTES'])+1)
	end=(data.index(['**PREFERENCES'])-1)
	for r in range(start, end):
		attribute={'name':data[r][1].strip('+:')}
		header.append(attribute)
	decision=data[end-1][1]
	end=(data.index(['**EXAMPLES']))
	
	start=(data.index(['**PREFERENCES'])+1)
	for r in header:
		r['preference']=data[start+j][1]
		j=j+1
	return header	


def collect_examples (data):
	"Collect examples values and put them in a matrix (list of lists) "
	
	matrix=[]
	data=[r for r in data if not '?' in r] #filter objects with " ?"
#	data=[data.remove(r) for r in data if data.count(r)>1] 
	start=(data.index(['**EXAMPLES'])+1)
	end=data.index(['**END'])
	for i in range(start, end):
		data[i]=(map(float, data[i]))
		matrix.append(data[i])
	i=1
	for r in matrix:
		r.insert(0, str(i))
		i=i+1
##	matrix=[list(i) for i in set(tuple(j) for j in matrix)] #remove duplicate example   
	return matrix


def FileToInfoSystem(isf):
	"Read *.isf file and copy it's values in Infosystem dictionary"
	data=[]
	try:
		infile=open(isf,"r")
		rows=infile.readlines()
		for line in rows:
			line=(line.split())
			if (len(line)>0 ):
				data.append(line)
		infile.close()	   
		infosystem={'attributes':collect_attributes(data),'examples':collect_examples(data)}
	except TypeError:
		print "\n\n Computing error or input file %s is not readeable. Exiting gracefully" % isf
		sys.exit(0)
		
	return infosystem


def UnionOfClasses (infosystem):
	"Find upward and downward union for all classes and put it in a dictionary"
	DecisionClass=[]
	AllClasses=[]
	matrix=infosystem['examples']
	for r in matrix:
		DecisionClass.append(int(r[-1]))
	DecisionClass=list(set(DecisionClass))
	for c in range(len(DecisionClass)):
		tmplist=[r for r in matrix if int(r[-1])==DecisionClass[c]]
		AllClasses.append(tmplist)  
	return AllClasses,DecisionClass
	
	
def DownwardUnionsOfClasses (infosystem):
	"For each decision class, downward union corresponding to a decision class\
	is composed of this class and all worse classes (<=)"

	DownwardUnionClass=[]
	DecisionClass=[]
	matrix=infosystem['examples']
	for r in matrix:
		DecisionClass.append(int(r[-1]))
	DecisionClass=list(set(DecisionClass))
	for c in DecisionClass:
		tmplist=[r for r in matrix if int(r[-1])<=c]
		DownwardUnionClass.append(tmplist)
		#label=[row[0] for row in tmplist]
	return DownwardUnionClass


def UpwardUnionsOfClasses (infosystem):
	"For each decision class, upward union corresponding to a decision class \
	is composed of this class and all better classes.(>=)"

	UpwardUnionClass=[]
	DecisionClass=[]
	matrix=infosystem['examples']
	for r in matrix:
		DecisionClass.append(int(r[-1]))
	DecisionClass=list(set(DecisionClass))
	for c in DecisionClass:
		tmplist=[r for r in matrix if int(r[-1])>=c]
		UpwardUnionClass.append(tmplist)
		#label=[row[0] for row in tmplist]
	return UpwardUnionClass
  
	
###############################
def is_better (r1,r2, preference):
	"Check if r1 is better than r2"
	return all((( x >=y and p=='gain') or (x<=y and p=='cost')) for x,y, p in zip(r1,r2, preference) ) 
	

def is_worst (r1,r2, preference): 
	"Check if r1 is worst than r2"
	return all((( x <=y and p=='gain') or (x>=y and p=='cost')) for x,y, p in zip(r1,r2, preference) ) 
 #################################


def DominatingSet (infosystem):
	"Find P-dominating set"
	matrix=infosystem['examples']
	preference=[s['preference'] for s in infosystem['attributes'] ]
	Dominating=[]
	for row in matrix:
		examples=[r  for r in matrix if  is_better(r[1:-1], row[1:-1], preference) ] 
		Dominating.append({'object':row[0], 'dominance':[i[0] for i in examples], 'examples':examples})
##	for dom in Dominating:
##		print  dom['dominance'] ,' dominating ', dom['object'] 
	return Dominating

def DominatedSet (infosystem):
	"Find P-Dominated set"
	matrix=infosystem['examples']
	preference=[s['preference'] for s in infosystem['attributes'] ]
	Dominated=[]
	for row in matrix:
		examples=[r  for r in matrix if  is_worst(r[1:-1], row[1:-1], preference[:-1]) ] 
		Dominated.append({'object':row[0], 'dominance':[i[0] for i in examples], 'examples':examples})
##	for dom in Dominated:
##		print  dom['dominance'] ,' is dominated by ', dom['object'] 
	return Dominated


def LowerApproximation (UnionClasses,  Dom, DecisionClass):
	"Find Lower approximation and return a dictionaries list"
	c=0
	LowApprox=[]
	single=dict()
	for union in UnionClasses:
		tmp=[]
		UClass=set([row[0] for row in union] )
		for d in Dom:
			if (UClass.issuperset(set(d['dominance']))): #if Union class is a superse of dominating/dominated set, =>single Loer approx.
				tmp.append(d['object'])
		single={'class':DecisionClass[c], 'objects':tmp} #dictionary for lower approximation  -- 
		LowApprox.append(single) #insert all Lower approximation in a list
		c+=1
	return LowApprox


def UpperApproximation (UnionClasses,  Dom, DecisionClass):
	"Find Upper approximation and return a dictionaries list"
	c=0
	UppApprox=[]
	single=dict()
	for union in UnionClasses:
		UnClass=[row[0] for row in union]  #single union class
		s=[]
		for d in Dom:
		   if len(set(d['dominance']) & set(UnClass)) >0:
			   s.append(d['object'])
		single={'class':DecisionClass[c],'objects':list(set(s))}
		UppApprox.append(single)
		c+=1

	return UppApprox


def Boundaries (UppApprox, LowApprox):
	"Find Boundaries like doubtful regions"
	Boundary=[]
	single=dict()

	for i in range(len(UppApprox)):
		single={'class':i, 'objects':list (set(UppApprox[i]['objects'])-set(LowApprox[i]['objects']) )}
		Boundary.append(single)
		
	return Boundary
	

def AccuracyOfApproximation(UppApprox, LowApprox):
	"Define the accuracy of approximation of Upward and downward approximation class"
	return len(LowApprox)/len(UppApprox)
	
	
def QualityOfQpproximation(DownwardBoundary,  infosystem):
	"Defines the quality of approximation of the partition Cl or, briefly, the quality of sorting"
	UnionBoundary=set()
	U=set([i[0] for i in infosystem['examples']])
	for b in DownwardBoundary:
		UnionBoundary=set(UnionBoundary) | set(b['objects'])
	return float(len(U-UnionBoundary)) / float(len(U))
	
	
def FindObjectCovered (rules, selected):
	"Find objects covered by a single rule and return\
	all related examples covered"
	obj=[]
	examples=[]
	
	for rule in rules:
		examples.append(rule['objectsCovered']) 

	if len(examples)>0:
		examples = reduce(set.intersection,map(set,examples))  #functional approach: intersect all lists if example is not empty
		examples = list(set(examples) & set([r[0] for r in selected]))
	return examples #all examples covered from a single rule
	
		
def Evaluate (elem,rules,G,selected,infosystem):
	"Calcolate first and second evaluate index, according with original DOMLEM Algorithm"
	tmpRules=copy.deepcopy(rules)
	tmpElem=copy.deepcopy(elem)
	tmpRules.append(tmpElem)
	Object=[]
	Object=FindObjectCovered(tmpRules,selected)
	if(float(len(Object)))>0:
			 firstEvaluate=float(len(set(G) & set(Object))) / float(len(Object))
			 secondEvaluate=float(len(set(G) & set(Object)))
	else:
			 firstEvaluate=0
			 secondEvaluate=0
			 
	return firstEvaluate,secondEvaluate
	
	
  
def FindBestCondition (best, elem, rules, selected, G, infosystem):
	"Choose the best condition"

	firstElem,secondElem=Evaluate(elem,rules,G,selected,infosystem)
	firstBest,secondBest=Evaluate(best,rules,G,selected,infosystem)
	
	if (firstElem>firstBest) or (firstElem==firstBest and secondElem>=secondBest):
		best=copy.deepcopy(elem)
	else:
		best=best

	return best

	
def Type_one_rule (c,  e,  preference,  matrix):
	elem={'criterion':c,'condition':e, 'sign':preference[c-1],'class':'', \
	'objectsCovered':[r[0] for r in matrix if (((r[c] >= e ) and (preference[c-1] == 'gain')) \
																			  or ((r[c] <= e ) and (preference[c-1] == 'cost' )))],'label':''}
	return elem
	
def Type_three_rule (c,  e,  preference,  matrix):
	elem={'criterion':c,'condition':e, 'sign':preference[c-1],'class':'', \
	'objectsCovered':[r[0] for r in matrix if (((r[c] <= e ) and (preference[c-1] == 'gain')) \
														or ((r[c] >= e ) and (preference[c-1] == 'cost' )))],'label':''}
	return elem


def Find_rules (B, infosystem, type_rule):
	"Search rule from a family of lower approximation of upward unions \
	of decision classes"
	start=time()
	matrix=copy.deepcopy(infosystem['examples'])
	criteria_num=len(infosystem['attributes'])
	criteria=[r[1:-1] for r in matrix]
	preference=[s['preference'] for s in infosystem['attributes'] ] #extract preference label
	num_rules=0 #total rules number for each lower approximation
	G=copy.deepcopy(B)		#a set of objects from the given approximation
	E=[]		#a set  of rules covering set B (is a list of dictionary)
	all_obj_cov_by_rules=[] #all objects covered by all rules in E
	selected=copy.deepcopy(matrix) #storage reduct matrix by single elementary condition
	while (len(G)!=0 ):
		rules=[]	 #starting comples (single rule built from elementary conditions  )
		S=copy.deepcopy(G)		 #set of objects currently covered by rule
		control=0
		while (len(rules)==0 or set(obj_cov_by_rules).issubset(B)==False):
			obj_cov_by_rules=[] #set covered by rules
			best={'criterion':'','condition':'','sign':'','class':'','objectsCovered':'','label':'', 'type':''} #best candidate for elementary condition - start as empty
			for c in range(1, criteria_num):
				Cond=[r[c] for r in selected if  r[0] in S] #for each positive object from S create an elementary condition
				for e in Cond:
					if type_rule=='one':
						elem= Type_one_rule (c,  e,  preference,  matrix)
					elif type_rule=='three':
						elem= Type_three_rule (c,  e,  preference,  matrix)
					else:
						elem={'criterion':'','condition':'','sign':'','class':'','objectsCovered':'','label':'', 'type':''}
					best=FindBestCondition(best, elem, rules, selected, G, infosystem)
			if best not in rules:
				rules.append(best)   #add the best condition to the complex 

			for r in rules:
				obj_cov_by_rules.append(r['objectsCovered'])
			obj_cov_by_rules=list((reduce(set.intersection,map(set,obj_cov_by_rules)))) #reduce():Apply function of two arguments cumulatively to the items of iterable, from left to right, so as to reduce the iterable to a single value.

			S=list(set(S) & set(best['objectsCovered'] ))
			control+=1

#		rules=CheckMinimalCondition (rules,B,matrix)
		
		if rules not in E:
			E.append(rules) #add the induced rule
			num_rules+=1
		all_obj_cov_by_rules=list(set(all_obj_cov_by_rules) | set(obj_cov_by_rules))
		
		G=list(set(B)-set(all_obj_cov_by_rules)) #remove example coverred by all finded rule -- this operation is a set difference
		selected=[o  for o in selected if  not  o[0]  in all_obj_cov_by_rules] #reduct matrix, remove object coverred by all finded rule
		num_rules+=1
		
	return E
	
	
  
def Domlem(Lu,Ld, infosystem):
	"DOMLEM algoritm \
	(An algorithm for induction of decision rules consistent with the dominance\
	principle - Greco S., Matarazzo, B., Slowinski R., Stefanowski J.)"
	attributes=infosystem['attributes']
	
	RULES=[]

##  *** AT MOST {<= Class} - Type 3 rules ***"		
	for b in Ld[:-1]:
		B=b['objects']
		E=Find_rules(B, infosystem, 'three')
		for e in E:
			for i in e:
				i['class']=b['class']
				i['label']=attributes[i['criterion']-1]['name']
				i['type']='at_most'
				if (attributes[i['criterion']-1]['preference']=='gain'):
					i['sign']='<='
				else:
					i['sign']='>='
			RULES.append(e)
			
## *** AT LEAST {>= Class} - Type 1 rules *** "
	for a in Lu[1:]:
		B=a['objects']
		E=Find_rules (B, infosystem, 'one')
		for e in E:
			for i in e:
				i['class']=a['class']
				i['label']=attributes[i['criterion']-1]['name']
				i['type']='at_least'
				if (attributes[i['criterion']-1]['preference']=='gain'):
					i['sign']='>='
				else:
					i['sign']='<='
			RULES.append(e)	   

	return RULES
	
	
def Print_rules(RULES, outputTxt):
	"Print rls output file"
	i=1
	outfile=open(outputTxt+".rls","w")
	outfile.write('[RULES]\n')
	
	for R in RULES:
		outfile.write("%d: " % i, )
		for e in R:
				outfile.write("( %s %s %.3f )" % (e['label'], e['sign'],e['condition'] ))
		outfile.write("=> ( class %s , %s )\n" % ( e['type'], e['class'] ))
		i+=1
	outfile.close()
	return 0

def Parser_mapcalc(RULES, outputMap):
	"Parser to build a formula to be included  in mapcalc command"
	i=1
	category=[]
	maps=[]
	stringa=[]
	
	for R in RULES: 
		formula="if("
		for e in R[:-1]: #build a mapcalc formula
			formula+= "(%s %s %.4f ) && " % (e['label'],  e['sign'] ,  e['condition'] )
		formula+= "(%s %s %.4f ),%d,null())" % (R[-1]['label'],R[-1]['sign'], R[-1]['condition'] ,i )
		mappa="r%d_%s_%d" % ( i, R[0]['type'], R[0]['class'] ) #build map name for mapcalc output
		category.append({'id':i, 'type': R[0]['type'], 'class':R[0]['class']}) #extract category name
		maps.append(mappa) #extract maps name
		grass.mapcalc(mappa +"=" +formula)
		i+=1
	mapstring=",".join(maps)

	
	#make one layer for each label rule
	labels=["_".join(m.split('_')[1:]) for m in maps] 
	labels=list(set(labels))
	for l in labels:
		print "mapping %s rule" % str(l)
		map_synth=[]
		for m in maps:
			if l == "_".join(m.split('_')[1:]):
				map_synth.append(m)
		if len(map_synth)>1:
			grass.run_command("r.patch", overwrite='True', input=(",".join(map_synth)), output=l )
		else:
			grass.run_command("g.copy",rast=("".join(map_synth),l))
		grass.run_command("r.to.vect", overwrite='True', flags='s', input=l, output=l, feature='area')
		grass.run_command("v.db.addcol", map=l, columns='rule varchar(25)')
		grass.run_command("v.db.update", map=l, column='rule', value=l)
		grass.run_command("v.db.addcol", map=l, columns='label varchar(25)')
		grass.run_command("v.db.update", map=l, column='label', value=l)
	mapslabels=",".join(labels)
	

	if len(maps)>1:
		grass.run_command("v.patch", overwrite='True', flags='e', input=mapslabels, output=outputMap)
	else:
		grass.run_command("g.copy",vect=(mapslabels,outputMap))
			
	if not flags['l']:
		grass.run_command("g.remove",  rast=mapstring)
		grass.run_command("g.remove",  vect=mapstring)
		
	return 0
	  
	
def main():
	"main function for DOMLEM algorithm"
	try:
		start=time()
		attributes = options['criteria'].split(',')
		preferences=options['preferences'].split(',')
		decision=options['decision']
		outputMap= options['outputMap']
		outputTxt= options['outputTxt']
		out=BuildFileISF(attributes, preferences, decision, outputMap, outputTxt)
		infosystem=FileToInfoSystem(out)
		
		AllClasses,DecisionClass=UnionOfClasses(infosystem)
		DownwardUnionClass=DownwardUnionsOfClasses(infosystem)
		UpwardUnionClass=UpwardUnionsOfClasses(infosystem)
		Dominating=DominatingSet(infosystem)
		Dominated=DominatedSet(infosystem)
	##	upward union class
		print "elaborate upward union"
		Lu=LowerApproximation(UpwardUnionClass, Dominating, DecisionClass) #lower approximation of upward union for type 1 rules
		Uu=UpperApproximation(UpwardUnionClass,Dominated, DecisionClass ) #upper approximation of upward union
		UpwardBoundary=Boundaries(Uu, Lu)
	##	downward union class
		print "elaborate downward union"
		Ld=LowerApproximation(DownwardUnionClass, Dominated, DecisionClass) # lower approximation of  downward union for type 3 rules
		Ud=UpperApproximation(DownwardUnionClass,Dominating, DecisionClass ) # upper approximation of  downward union 
		DownwardBoundary=Boundaries(Ud, Ld)
		QualityOfQpproximation(DownwardBoundary,  infosystem)
		print "RULES extraction (*)" 
		RULES=Domlem(Lu,Ld, infosystem)
		Parser_mapcalc(RULES, outputMap)		   
		Print_rules(RULES, outputTxt)
		end=time()
		print "Time computing-> %.4f s" % (end-start)
		return 0
	except:
		print "ERROR! Rules does not generated!"
		sys.exit()

if __name__ == "__main__":
	options, flags = grass.parser()
	sys.exit(main())