#! /bin/env python
from __future__ import division
import argparse
import os
import functools
from solartherm import postproc
from solartherm import simulation

import numpy as np
import pylab as pl
import random
import time
#import multiprocessing
import scoop
from scoop import futures
import matplotlib.pyplot as plt

t_start = time.time()

try:
	from scipy import optimize as sciopt
except ImportError:
	sciopt = None
try:
	from pyswarm import pso
except ImportError:
	pso = None
try:
	import cma
except ImportError:
	cma = None
try:
	import pyevolve as pe
	from pyevolve import GSimpleGA, G1DList, Initializators, Crossovers, Selectors, Mutators, Scaling, Consts
except ImportError:
	pe = None

try:
	import deap
	from deap import algorithms, base, creator, tools
except ImportError:
	deap = None

#	objfunc = functools.partial(objective_function, sim, mat_fn,
#			(args.start, args.stop, args.step),
#			args.v, scale, offset, perf_i, par_n, resultclass, op_meth, frwk, sign)

def objective_function(sim, mat_fn, stime, initStep, maxStep, integOrder, solver, nls, lv, verb, scale, offset, perf_i, par_n, resultclass, op_meth, frwk, sign, par_val):
	par_v = [str(v*scale[i] + offset[i]) for i, v in enumerate(par_val)]
	sim.update_pars(par_n, par_v)
	sim.simulate(start=stime[0], stop=stime[1], step=stime[2], initStep=initStep, maxStep=maxStep, integOrder=integOrder, solver=solver, nls=nls, lv=lv)

	res = resultclass(mat_fn)
	constr, distance = res.constrained_optimisation() # constr is true if optimisattion is constrained. distance to be added to a constant penalty offset

	if frwk == "soo_min": # Single-objective optimisation to minimise an objective
		obj_min = res.calc_perf()[perf_i[0]] # The objective to be minimised (i.e. lcoe or lcof)

		if constr:
			obj_min = 9.999e+99 + distance
			print "A constraint is violated at this design point!"

		if verb:
			print(par_v)
			print(obj_min)

		if op_meth == 'ga2':
			return obj_min,
		else:
			return obj_min
	elif frwk == "soo_max": # Single-objective optimisation to maximise an objective
		obj_max = res.calc_perf()[perf_i[0]] # The objective to be maximised (i.e. epy, capf or srev)

		if constr:
			obj_max = -9.999e+99 - distance
			print "A constraint is violated at this design point!"

		if verb:
			print(par_v)
			print(obj_max)

		if op_meth == 'ga2':
			return obj_max,
		else:
			return sign * obj_max
	else: # Multi-objective optimisation
		assert(frwk == "moo"), "The framework specified for optimisation must be one of soo_min, soo_max, or moo!"
		if perf_i[0] == 1:
			obj_min = res.calc_perf()[perf_i[0]]
			obj_max = res.calc_perf()[perf_i[1]]
		else:
			obj_min = res.calc_perf()[perf_i[1]]
			obj_max = res.calc_perf()[perf_i[0]]

		if constr:
			obj_min = 9.999e+99 + distance
			obj_max = -9.999e+99 - distance
			print "A constraint is violated at this design point!"

		if verb:
			print(par_v)
			print(obj_min, obj_max)

		return obj_min, \
		obj_max


if __name__ == '__main__':
	"""
	Should make sure parameters are not final (protected), or that other
	derived parameters are final.
	"""
	parser = argparse.ArgumentParser()
	parser.add_argument('file',
			help='model file name')
	parser.add_argument('-v', action='store_true',
			help='verbose')
	parser.add_argument('--nomc', action='store_true',
			help='no model compilation')
	parser.add_argument('--nosc', action='store_true',
			help='no simulation compilation')
	parser.add_argument('--noc', action='store_true',
			help='no compilation at all (--nomc and --nosc)')
	parser.add_argument('--start', type=str, default='0',
			help='simulation start time: <number>[,y,d,m,s]')
	parser.add_argument('--stop', type=str, default='1y',
			help='simulation stop time: <number>[,y,d,m,s]')
	parser.add_argument('--step', type=str, default='5m',
			help='simulation time step: <number>[,y,d,m,s]')
	parser.add_argument('--initStep', type=str, default=None,
			help='simulation initial step size used by dassl and ida methods: <number>[,y,d,m,s]')
	parser.add_argument('--maxStep', type=str, default=None,
			help='simulation maximum absolute step size used by dassl and ida methods: <number>[,y,d,m,s]')
	parser.add_argument('--integOrder', type=str, default='5',
			help='simulation maximum integration order used by dassl and ida methods: <number>')
	parser.add_argument('--solver', type=str, default='dassl',
			help='solver choice for OpenModelica')
	parser.add_argument('--nls', type=str, default='homotopy',
			help='non-linear solver choice for OpenModelica e.g. newton, hybrid, kinsol, mixed, and homotopy')
	parser.add_argument('--lv', type=str, default='-LOG_SUCCESS,-stdout',
			help='a comma-separated String list specifing which logging levels to enable, e.g. LOG_DASSL,LOG_SOLVER etc')
	parser.add_argument('--method', type=str, default='Nelder-Mead',
			help='pso,  cma, ga1, ga2, nsga2 or one of the scipy optimisation methods')
	parser.add_argument('--maxiter', type=int, default=20,
			help='maximum number of iterations (not necessarily number of simulations)')
	parser.add_argument('--objective', type=str, default='lcoe',
			help='quantity to conduct a single (i.e. minimisation or maximisation) or multi-objective optimisation (produced by post processing) in form of objective=lcoe,capf')
	parser.add_argument('par', metavar='P', type=str, nargs='*',
			help='parameters with bounds and optional starting value in form PNAME=LOW,HIGH[,START]')
	parser.add_argument('--fuel', action='store_true',
			help='run post-processing calculations for levelised cost of fuel')
	parser.add_argument('--framework', type=str, default='soo_min',
			help='soo_min, soo_max or moo, being single objective minimisation, single objective maximisation, or multi-objective optimisation respectively')
	parser.add_argument('--outfile_f',  type=str, default=None,
			help='save final optimal solution and objective(s) to outfile_f')
	parser.add_argument('--outfig',  type=str, default=None,
			help='save figure to outfig instead of displaying')
	parser.add_argument('--outfile_p',  type=str, default=None,
			help='save pareto optimal solutions of multi-objective optimisation to outfile_p')
	parser.add_argument('--dm',  type=str, default='linmap',
			help='Decision-making methods for the multi-objective optimisation framework: linmap, topsis')
	args = parser.parse_args()



	sim = simulation.Simulator(args.file);

	if not args.noc:
		if not args.nomc:
			print('Compiling model')
			sim.compile_model()
		if not args.nosc:
			print('Compiling simulator')
			sim.compile_sim(args=([] if args.v else ['-s']))

	if args.fuel:
		resultclass = postproc.SimResultFuel
	else:
		resultclass = postproc.SimResultElec

	op_meth = args.method
	frwk = args.framework

	if args.framework == "soo_max" and op_meth not in ["ga1", "ga2"]:
		sign = -1.0
	else:
		sign = 1.0

	if args.framework == "moo":
		decisionmaker= postproc.DecisionMaker # Decision-maker classs instance
		dm_method = args.dm

	sim.load_init()

	par_n = [] # names
	par_b = [] # bounds
	par_0 = [] # start
	offset = []
	scale = []
	# Don't need bounds (pass None for variable to minimize)
	for pp in args.par:
		k, v = pp.split('=')
		par_n.append(k)
		vals = [simulation.parse_var_val(vv, sim.get_unit(k))
			for vv in v.split(',')]
		assert len(vals) >= 2, 'Expected parameter bounds + optional start value'
		lb = vals[0]
		ub = vals[1]
		assert lb <= ub, 'Lower bound greater than upper bound'
		#par_b.append([lb, ub])
		p0 = (ub + lb)/2
		if len(vals) == 3:
			p0 = vals[2]
		#par_0.append(p0)
		offset.append(lb)
		scale.append((ub - lb))
		par_b.append([0, 1])
		par_0.append((p0 - lb)/(ub - lb))

	mat_fn = sim.model + '_res.mat'

	try:
		obj_n = args.objective.split(",") # A list of objective(s) name in string
		perf_i = [] # perfromance index
		for oo in obj_n:
			perf_i.append(resultclass.perf_n.index(oo))
	except ValueError:
		raise ValueError('Objective(s) value should be in '
				+ str(resultclass.perf_n))

	objfunc = functools.partial(objective_function, sim, mat_fn,
			(args.start, args.stop, args.step), args.initStep,
			args.maxStep, args.integOrder, args.solver, args.nls, args.lv,
			args.v, scale, offset, perf_i, par_n, resultclass, op_meth, frwk, sign)

	print "\n\n\nOptimisation parameter(s): ", par_n, "\n\n\n"

	if args.framework in ['soo_min', 'soo_max']:
		if args.method == 'pso':
			assert pso is not None, 'Library for pso is not installed'
			swarmsize=5
			lb = [v[0] for v in par_b]
			ub = [v[1] for v in par_b]
			res = pso(objfunc, lb, ub, maxiter=args.maxiter, swarmsize=5)
			cand = [scale[i]*v + offset[i] for i, v in enumerate(res[0])]

			t_end = time.time() # NOTE: the clock does not work on the multicored version!
			t_dur = t_end - t_start # Time elapsed to (succesfully) finish the optimisation

			print "\n\nTotal time elapsed: ", t_dur, "seconds."
			if args.outfile_f is not None:
				f = open(args.outfile_f, 'w') # Example: outfile_f = ../examples/result.txt
				f.write("Optimal design parameters: " + str(cand) + '\n')
				if args.framework == 'soo_min':
					f.write("Optimal objective function: " + str(res[1]))
				else:
					f.write("Optimal objective function: " + str(sign * res[1]))
				f.close()
			else:
				if args.framework == 'soo_min':
					print "Optimal design parameters: ", (cand)
					print "Optimal objective function: ", (res[1])
				else:
					print "Optimal design parameters: ", (cand)
					print "Optimal objective function: ", (sign * res[1])
		elif args.method == 'cma':
			assert cma is not None, 'Library for cma is not installed'
			sigma0 = 0.2 # "step size", should be around 1/4 of search domain
			popsize = 5
			lb = [v[0] for v in par_b]
			ub = [v[1] for v in par_b]
			res = cma.fmin(objfunc, par_0, sigma0,
					restarts=0,
					options={
							'bounds': [lb, ub],
							#'maxfevals': args.maxiter,
							'maxiter': args.maxiter,
							'popsize': popsize,
					})
			cand = [scale[i]*v + offset[i] for i, v in enumerate(res[0])]

			t_end = time.time() # NOTE: the clock does not work on the multicored version!
			t_dur = t_end - t_start # Time elapsed to (succesfully) finish the optimisation

			print "\n\n\Total time elapsed: ", t_dur, "seconds."
			if args.outfile_f is not None:
				f = open(args.outfile_f, 'w') # Example: outfile_f = ../examples/result.txt
				f.write("Optimal design parameters: " + str(cand) + '\n')
				if args.framework == 'soo_min':
					f.write("Optimal objective function: " + str(res[1]))
				else:
					f.write("Optimal objective function: " + str(sign * res[1]))
				f.close()
			else:
				if args.framework == 'soo_min':
					print "Optimal design parameters: ", (cand)
					print "Optimal objective function: ", (res[1])
				else:
					print "Optimal design parameters: ", (cand)
					print "Optimal objective function: ", (sign * res[1])
		elif args.method == 'ga1':
			assert pe is not None, 'Library for pyevolve is not installed'

			lb = [v[0] for v in par_b] # Normalised lower bound
			ub = [v[1] for v in par_b] # Normalised upper bound

			ind_size = len(par_n) # Number of design parameters
			pop_size = 100 # Set number of individuals in population
			ngen = 2 # Total number of generation to run
			cxpb = 0.98 # Crossover probability
			mutpb = 0.01 # Mutation probability
			freq_stats = 10 # Frequency of stats
			paral_eval = False # To enable parallel evaluation. Only use it when the fitness function is slow!

			pl.ion()

			genome = G1DList.G1DList(ind_size) # Genome instance
			genome.setParams(rangemin = lb[0], rangemax = ub[0]) # Set the range min and max of the 1D List
			genome.evaluator.set(objfunc) # The evaluator function (fitness/objective function)	
			genome.initializator.set(Initializators.G1DListInitializatorReal) # Real initialization function of G1DList
			genome.crossover.set(Crossovers.G1DListCrossoverUniform) # The G1DList Uniform Crossover
			genome.mutator.set(Mutators.G1DListMutatorRealRange) # Simple real range mutator for G1DList

			ga = GSimpleGA.GSimpleGA(genome) # Genetic algorithm instance
			ga.selector.set(Selectors.GTournamentSelector) # Set the selector method
			if args.framework == 'soo_min':
				ga.setMinimax(Consts.minimaxType["minimize"])
			else:
				ga.setMinimax(Consts.minimaxType["maximize"])
			ga.setPopulationSize(pop_size) # Set the population size for each generation
			ga.setGenerations(ngen) # Set the number of generation
			ga.setCrossoverRate(cxpb) # Set the crossover rate  
			ga.setMutationRate(mutpb) # Set the mutation rate
			ga.terminationCriteria.set(GSimpleGA.ConvergenceCriteria) # Terminate the evolution when the population have converged
			pop = ga.getPopulation() # Return the internal population of GA Engine
			pop.scaleMethod.set(Scaling.SigmaTruncScaling) # Sigma Truncation scaling scheme, allows negative scores
			ga.evolve(freq_stats=freq_stats) # Run the optimisation and print the stats of the ga every n generation
			ga.setMultiProcessing(flag=paral_eval, full_copy=False) # Set the flag to enable/disable the use of python multiprocessing module

			res = ga.bestIndividual() # Best individual in normalised form
			cand = [scale[i]*v + offset[i] for i, v in enumerate(res.genomeList)] # Denormalised best individual
			if args.outfile_f is not None:
				f = open(args.outfile_f, 'w') # Example: outfile_f = ../examples/result.txt
				f.write("Optimal design parameters: " + str(cand) + '\n')
				f.write("Optimal objective function: " + str(res.score))
				f.close()
			else:
				print "Optimal design parameters: ", (cand)
				print "Optimal objective function: ", (res.score)
		elif args.method == 'ga2':
			assert deap is not None, 'Library for deap is not installed'

			lb = [v[0] for v in par_b] # Normalised lower bound
			ub = [v[1] for v in par_b] # Normalised upper bound

			ind_size = len(par_n) # Number of design parameters
			pop_size = 100 # Set number of individuals in population
			ngen = 2 # Total number of generation to run
			cxpb = 0.98 # Crossover probability
			cxpb_in = 0.5 # Probability of crossover within individual
			mutpb = 0.02 # Mutation probability
			mutpb_in = 0.01 # Probability of mutation within individual
			paral_eval = False # To enable parallel evaluation. Only use it when the fitness function is slow!
			n_cores = 4 # Number of cores for parallel evaluation of the fitness function

			if args.framework == 'soo_min':
				creator.create("FitnessMin", base.Fitness, weights=(-1.0,))
			else:
				creator.create("FitnessMin", base.Fitness, weights=(1.0,))
			creator.create("Individual", list, fitness=creator.FitnessMin)

			toolbox = base.Toolbox() # Attribute generator
			toolbox.register("attr_float", random.uniform, lb[0], ub[0]) # Generate a vector of uniform random numbers within the lower and upper bounds of the design variables
			toolbox.register("individual", tools.initRepeat, creator.Individual, toolbox.attr_float, ind_size) # Structure initializer for individuals
			toolbox.register("population", tools.initRepeat, list, toolbox.individual) # Structure initializer for the population

			# Decorator to obtain a constrained domain that is applied to the mutation and crossover
			#def checkBounds(min, max):
				#def decorator(func):
					#def wrappper(*args, **kargs):
						#offspring = func(*args, **kargs)
						#for child in offspring:
							#for i in range(len(child)):
								#if child[i] > max:
									#child[i] = max
								#elif child[i] < min:
									#child[i] = min
						#return offspring
					#return wrappper
				#return decorator

			# NOTE: For penalty function-based constraint handling see: http://deap.readthedocs.io/en/master/tutorials/advanced/constraints.html

			# Register the GA operators
			toolbox.register("evaluate", objfunc) # Set the evaluation function
			toolbox.register("mate", tools.cxUniform, indpb=cxpb_in) # Set the crossover type
			toolbox.register("mutate", tools.mutFlipBit, indpb=mutpb_in) # Set the mutation type
			#toolbox.register("select", tools.selNSGA2) # Set the selector method
			toolbox.register("select", tools.selTournament, tournsize=3) # Set the selector method

			#toolbox.decorate("mate", checkBounds(lb[0], ub[0])) # Check bounds for a constrained domain
			#toolbox.decorate("mutate", checkBounds(lb[0], ub[0])) # Check bounds for a constrained domain

			# Parallel evaluation by running the computation multicore
			if paral_eval:
				#pool = multiprocessing.Pool(processes=n_cores)
				#toolbox.register("map", pool.map) # Change the map functions everywhere to toolbox.map to make the algorithm use a multicored map
				toolbox.register("map", futures.map) # Change the map functions everywhere to toolbox.map to make the algorithm use a multicored map

			pop = toolbox.population(n=pop_size) # Set the size of population (individuals)

			fitnesses = list(toolbox.map(toolbox.evaluate, pop)) # Evaluate the entire population

			for ind, fit in zip(pop, fitnesses):
				ind.fitness.values = fit  # Run the fitness (min mean on each of the individuals)

			fits = [ind.fitness.values[0] for ind in pop] # # Gather all the fitnesses in one list

			# Begin the evolution
			for g in range(ngen):
				print("-- Generation %i --" % g)

				# Select the next generation individuals
				offspring = toolbox.select(pop, len(pop)) # Select the best individuals in the population

				# Clone the selected individuals
				offspring = list(toolbox.map(toolbox.clone, offspring)) #Create the offspring

				# Apply crossover on the offspring
				for child1, child2 in zip(offspring[::2], offspring[1::2]):
					if random.random() < cxpb:
						toolbox.mate(child1, child2)
						del child1.fitness.values
						del child2.fitness.values

				# Apply mutation on the offspring
				for mutant in offspring:
					# mutate an individual with probability mutpb
					if random.random() < mutpb:
						toolbox.mutate(mutant)
						del mutant.fitness.values

				# Evaluate the individuals with an invalid fitness
				invalid_ind = [ind for ind in offspring if not ind.fitness.valid]
				fitnesses = toolbox.map(toolbox.evaluate, invalid_ind)
				for ind, fit in zip(invalid_ind, fitnesses):
					ind.fitness.values = fit

				# Replace the enitre population with the offspring
				pop[:] = offspring
				fits = [ind.fitness.values[0] for ind in pop]

				print(tools.selBest(pop, 1)[0]) # Best individual in the current generation
				print(tools.selBest(pop, 1)[0].fitness.values[0]) # Best fitness in the current generation

			# For parallel evaluation using multiprocessing
			#if paral_eval:
				#pool.close()
				#pool.join()

			res = tools.selBest(pop, 1)[0] # Best individual in normalised form
			cand = [scale[i]*v + offset[i] for i, v in enumerate(res)] # Denormalised best individual

			t_end = time.time() # NOTE: the clock does not work on the multicored version!
			t_dur = t_end - t_start # Time elapsed to (succesfully) finish the optimisation

			print "\n\nTotal time elapsed: ", t_dur, "seconds."
			if args.outfile_f is not None:
				f = open(args.outfile_f, 'w') # Example: outfile_f = ../examples/result.txt
				f.write("Optimal design parameters: " + str(cand) + '\n')
				f.write("Optimal objective function: " + str(res.fitness.values[0]))
				f.close()
			else:
				print "Optimal design parameters: ", (cand)
				print "Optimal objective function: ", (res.fitness.values[0])
		else: # Use one of the scipy optimisation methods. For single-objective maximisation use one of L-BFGS-B, TNC or SLSQP.
			if args.framework == 'soo_max':
				assert sciopt is not None, 'Library for scipy is not installed'
				op_meth = "L-BFGS-B" # Because the Nelder_Mead method does not take bounds, thereby leading to negative domains
			res = sciopt.minimize(objfunc, par_0, method=op_meth, bounds=par_b,
					options={
						#'maxfev': args.maxiter,
						'maxiter': args.maxiter,
						'disp': True,
					})
			print(res)
			cand = [scale[i]*v + offset[i] for i, v in enumerate(res.x)]

			t_end = time.time() # NOTE: the clock does not work on the multicored version!
			t_dur = t_end - t_start # Time elapsed to (succesfully) finish the optimisation

			print "\n\nTotal time elapsed: ", t_dur, "seconds."
			if args.outfile_f is not None:
				f = open(args.outfile_f, 'w') # Example: outfile_f = ../examples/result.txt
				f.write("Optimal design parameters: " + str(cand) + '\n')
				if args.framework == 'soo_min':
					f.write("Optimal objective function: " + str(res.fun))
				else:
					f.write("Optimal objective function: " + str(sign * res.fun))
				f.close()
			else:
				if args.framework == 'soo_min':
					print "Optimal design parameters: ", (cand)
					print "Optimal objective function: ", (res.fun)
				else:
					print "Optimal design parameters: ", (cand)
					print "Optimal objective function: ", (sign * res.fun)
	else:
		assert(frwk == "moo"), "The framework specified for optimisation must be one of soo_min, soo_max, or moo!"
		if args.method == 'nsga2':
			assert deap is not None, 'Library for deap is not installed'

			lb = [v[0] for v in par_b] # Normalised lower bound
			ub = [v[1] for v in par_b] # Normalised upper bound

			ind_size = len(par_n) # Number of design parameters
			pop_size = 100 # Set number of individuals in population
			ngen = 20 # Total number of generation to run
			cxpb = 0.98 # Crossover probability
			#cxpb_in = 0.5 # Probability of crossover within individual
			#mutpb = 0.02 # Mutation probability
			#mutpb_in = 0.01 # Probability of mutation within individual
			paral_eval = False # To enable parallel evaluation. Only use it when the fitness function is slow!
			n_cores = 4 # Number of cores for parallel evaluation of the fitness function

			creator.create("FitnessMulti", base.Fitness, weights=(-1.0, 1.0))
			creator.create("Individual", list, fitness=creator.FitnessMulti)

			toolbox = base.Toolbox() # Attribute generator
			toolbox.register("attr_float", random.uniform, lb[0], ub[0]) # Generate a vector of uniform random numbers within the lower and upper bounds of the design variables
			toolbox.register("individual", tools.initRepeat, creator.Individual, toolbox.attr_float, ind_size) # Structure initializer for individuals
			toolbox.register("population", tools.initRepeat, list, toolbox.individual) # Structure initializer for the population 

						# Decorator to obtain a constrained domain that is applied to the mutation and crossover
			#def checkBounds(min, max):
				#def decorator(func):
					#def wrappper(*args, **kargs):
						#offspring = func(*args, **kargs)
						#for child in offspring:
							#for i in range(len(child)):
								#if child[i] > max:
									#child[i] = max
								#elif child[i] < min:
									#child[i] = min
						#return offspring
					#return wrappper
				#return decorator

			# NOTE: For penalty function-based constraint handling see: http://deap.readthedocs.io/en/master/tutorials/advanced/constraints.html

			# Register the GA operators
			toolbox.register("evaluate", objfunc) # Set the evaluation function
			#toolbox.register("mate", tools.cxUniform, indpb=cxpb_in) # Set the crossover type
			toolbox.register("mate", tools.cxSimulatedBinaryBounded, low=lb[0], up=ub[0], eta=20.0) # Set the crossover type
			#toolbox.register("mutate", tools.mutFlipBit, indpb=mutpb_in) # Set the mutation type
			toolbox.register("mutate", tools.mutPolynomialBounded, low=lb[0], up=ub[0], eta=20.0, indpb=1.0/ind_size) # Set the mutation type
			toolbox.register("select", tools.selNSGA2) # Set the selector method

			#toolbox.decorate("mate", checkBounds(lb[0], ub[0])) # Check bounds for a constrained domain
			#toolbox.decorate("mutate", checkBounds(lb[0], ub[0])) # Check bounds for a constrained domain

			# Parallel evaluation by running the computation multicore
			if paral_eval:
				#pool = multiprocessing.Pool(processes=n_cores)
				#toolbox.register("map", pool.map) # Change the map functions everywhere to toolbox.map to make the algorithm use a multicored map
				toolbox.register("map", futures.map) # Change the map functions everywhere to toolbox.map to make the algorithm use a multicored map

			stats = tools.Statistics(lambda ind: ind.fitness.values) # Set the optimisation statistics
			stats.register("avg", np.mean, axis=0)
			stats.register("std", np.std, axis=0)
			stats.register("min", np.min, axis=0)
			stats.register("max", np.max, axis=0)

			logbook = tools.Logbook() # Logging data
			logbook.header = "gen", "evals", "std", "min", "avg", "max"

			pop = toolbox.population(n=pop_size) # Set the size of population (individuals)
			hof = tools.ParetoFront() # The Pareto front hall of fame containing all the non-dominated individuals that ever lived in the population

			# Evaluate the individuals with an invalid fitness
			invalid_ind = [ind for ind in pop if not ind.fitness.valid]
			fitnesses = toolbox.map(toolbox.evaluate, invalid_ind)
			for ind, fit in zip(invalid_ind, fitnesses):
				ind.fitness.values = fit

			# This is just to assign the crowding distance to the individuals
			# no actual selection is done
			pop = toolbox.select(pop, len(pop))

			record = stats.compile(pop)
			logbook.record(gen=0, evals=len(invalid_ind), **record)
			#print(logbook.stream)

			# Begin the generational process
			print "Starting optimisation..."
			for gen in range(ngen):
				print("-- Generation %i --" % gen)
				# Vary the population
				offspring = tools.selTournamentDCD(pop, len(pop))
				offspring = [toolbox.clone(ind) for ind in offspring]

				for ind1, ind2 in zip(offspring[::2], offspring[1::2]):
					if random.random() <= cxpb:
						toolbox.mate(ind1, ind2)

						toolbox.mutate(ind1)
						toolbox.mutate(ind2)
						del ind1.fitness.values, ind2.fitness.values

				# Evaluate the individuals with an invalid fitness
				invalid_ind = [ind for ind in offspring if not ind.fitness.valid]
				fitnesses = toolbox.map(toolbox.evaluate, invalid_ind)
				for ind, fit in zip(invalid_ind, fitnesses):
					ind.fitness.values = fit

				# Select the next generation population
				pop = toolbox.select(pop + offspring, pop_size)
				record = stats.compile(pop)
				logbook.record(gen=gen, evals=len(invalid_ind), **record)
				#print(logbook.stream)

			# For parallel evaluation using multiprocessing
			#if paral_eval:
				#pool.close()
				#pool.join()

			pop.sort(key=lambda x: x.fitness.values)

			solutions = pop # Normalised non-dominated optimal individuals
			fitnesss = [ind.fitness.values for ind in pop] # optimal values of objective functions

			cand = []
			for ind in pop:
				cand.append([scale[i]*v + offset[i] for i, v in enumerate(ind)]) # Denormalised non-dominated optimal individuals

			cands = np.array(cand)
			front = np.array(fitnesss)

			t_end = time.time() # NOTE: the clock does not work on the multicored version!
			t_dur = t_end - t_start # Time elapsed to (succesfully) finish the optimisation

			# Save the optimal solutions to a text file
			if args.outfile_p is not None:
				f = open(args.outfile_p, 'w') # Example: outfile_p = ../examples/result.txt
				for nn in par_n:
					f.write("%s " % nn)
				for ii in obj_n:
					f.write("%s " % ii)
				f.write("\n")
				np.savetxt(f, np.c_[cands,front],delimiter=' ')
				f.close()

			dm = decisionmaker(cand,fitnesss) # Decision-making class instance
			if dm_method == 'linmap':
				best_ind, best_fitness = dm.linmap()
			else:
				assert(dm_method == 'topsis'), "Decision-making methods must be one of linmap or topsis!"
				best_ind, best_fitness = dm.topsis()

			print "\n\nTotal time elapsed: ", t_dur, "seconds."
			if args.outfile_f is not None:
				f = open(args.outfile_f, 'w') # Example: outfile_f = ../examples/result.txt
				f.write("Final optimal design parameters: " + str(best_ind) + '\n')
				f.write("Final optimal objective functions: " + str(best_fitness))
				f.close()
			else:
				print "Final optimal design parameters: ", (best_ind)
				print "Final optimal objective functions: ", (best_fitness)

			# Plot the Pareto Front
			fig = plt.figure()
			fig.add_subplot(111)

			plt.scatter(front[:,0], front[:,1], c="b", marker='*')
			plt.title('Pareto Front', loc='center')
			plt.axis("tight")
			if perf_i[0] == 1:
				obj_min_u = resultclass.perf_u[perf_i[0]] # objective unit
				obj_max_u = resultclass.perf_u[perf_i[1]] # objective unit
			else:
				obj_min_u = resultclass.perf_u[perf_i[1]] # objective unit
				obj_max_u = resultclass.perf_u[perf_i[0]] # objective unit
			plt.xlabel(obj_n[0] + " (" + obj_min_u + ")")
			plt.ylabel(obj_n[1] + " (" + obj_max_u + ")")
			plt.grid(True,color='0.5', linestyle='--', linewidth=0.5)
			if args.outfig is not None:
				fig.savefig(args.outfig, dpi=600) # Example: outfig = ../examples/ParetoFront.png
			else:
				plt.show(block=True)
		else:
			print "nsga2 is the only optimisation method implemented for the multi-objective optimisation framework"
