Index: /framspy/FramsticksEvolution.py
===================================================================
--- /framspy/FramsticksEvolution.py	(revision 1310)
+++ /framspy/FramsticksEvolution.py	(revision 1311)
@@ -6,6 +6,8 @@
 from FramsticksLib import FramsticksLib
 
+# Note: this may be less efficient than running the evolution directly in Framsticks, so if performance is key, compare both options.
 
-# Note: this may be less efficient than running the evolution directly in Framsticks, so if performance is key, compare both options.
+
+FITNESS_VALUE_INFEASIBLE_SOLUTION = -999999.0  # DEAP expects fitness to always be a real value (not None), so this special value indicates that a solution is invalid, incorrect, or infeasible. [Related: https://github.com/DEAP/deap/issues/30 ]. Using float('-inf') or -sys.float_info.max here causes DEAP to silently exit. If you are not using DEAP, set this constant to None, float('nan'), or another special/non-float value to avoid clashing with valid real fitness values, and handle such solutions appropriately as a separate case.
 
 
@@ -16,5 +18,5 @@
 		if actual_value > constraint_value:
 			if REPORT_CONSTRAINT_VIOLATIONS:
-				print('Genotype "%s" assigned low fitness because it violates constraint "%s": %s exceeds threshold %s' % (genotype, criterion_name, actual_value, constraint_value))
+				print('Genotype "%s" assigned a special ("infeasible solution") fitness because it violates constraint "%s": %s exceeds the threshold of %s' % (genotype, criterion_name, actual_value, constraint_value))
 			return False
 	return True
@@ -22,6 +24,6 @@
 
 def frams_evaluate(frams_lib, individual):
-	BAD_FITNESS = [-1] * len(OPTIMIZATION_CRITERIA)  # fitness of -1 is intended to discourage further propagation of this genotype via selection ("this genotype is very poor")
-	genotype = individual[0]  # individual[0] because we can't (?) have a simple str as a deap genotype/individual, only list of str.
+	FITNESS_CRITERIA_INFEASIBLE_SOLUTION = [FITNESS_VALUE_INFEASIBLE_SOLUTION] * len(OPTIMIZATION_CRITERIA)  # this special fitness value indicates that the solution should not be propagated via selection ("that genotype is invalid"). The floating point value is only used for compatibility with DEAP. If you implement your own optimization algorithm, instead of a negative value in this constant, use a special value like None to properly distinguish between feasible and infeasible solutions.
+	genotype = individual[0]  # individual[0] because we can't (?) have a simple str as a DEAP genotype/individual, only list of str.
 	data = frams_lib.evaluate([genotype])
 	# print("Evaluated '%s'" % genotype, 'evaluation is:', data)
@@ -34,5 +36,5 @@
 	except (KeyError, TypeError) as e:  # the evaluation may have failed for an invalid genotype (such as X[@][@] with "Don't simulate genotypes with warnings" option), or because the creature failed to stabilize, or for some other reason
 		valid = False
-		print('Problem "%s" so could not evaluate genotype "%s", hence assigned it low fitness: %s' % (str(e), genotype, BAD_FITNESS))
+		print('Problem "%s" so could not evaluate genotype "%s", hence assigned it a special ("infeasible solution") fitness value: %s' % (str(e), genotype, FITNESS_CRITERIA_INFEASIBLE_SOLUTION))
 	if valid:
 		default_evaluation_data['numgenocharacters'] = len(genotype)  # for consistent constraint checking below
@@ -43,11 +45,11 @@
 		valid &= genotype_within_constraint(genotype, default_evaluation_data, 'numgenocharacters', parsed_args.max_numgenochars)
 	if not valid:
-		fitness = BAD_FITNESS
+		fitness = FITNESS_CRITERIA_INFEASIBLE_SOLUTION
 	return fitness
 
 
 def frams_crossover(frams_lib, individual1, individual2):
-	geno1 = individual1[0]  # individual[0] because we can't (?) have a simple str as a deap genotype/individual, only list of str.
-	geno2 = individual2[0]  # individual[0] because we can't (?) have a simple str as a deap genotype/individual, only list of str.
+	geno1 = individual1[0]  # individual[0] because we can't (?) have a simple str as a DEAP genotype/individual, only list of str.
+	geno2 = individual2[0]  # individual[0] because we can't (?) have a simple str as a DEAP genotype/individual, only list of str.
 	individual1[0] = frams_lib.crossOver(geno1, geno2)
 	individual2[0] = frams_lib.crossOver(geno1, geno2)
@@ -56,5 +58,5 @@
 
 def frams_mutate(frams_lib, individual):
-	individual[0] = frams_lib.mutate([individual[0]])[0]  # individual[0] because we can't (?) have a simple str as a deap genotype/individual, only list of str.
+	individual[0] = frams_lib.mutate([individual[0]])[0]  # individual[0] because we can't (?) have a simple str as a DEAP genotype/individual, only list of str.
 	return individual,
 
@@ -62,4 +64,35 @@
 def frams_getsimplest(frams_lib, genetic_format, initial_genotype):
 	return initial_genotype if initial_genotype is not None else frams_lib.getSimplest(genetic_format)
+
+
+def is_feasible_fitness_value(fitness_value: float) -> bool:
+	assert isinstance(fitness_value, float), f"feasible_fitness({fitness_value}): argument is not of type float, it is of type {type(fitness_value)}"  # since we are using DEAP, we unfortunately must represent the fitness of an "infeasible solution" as a float...
+	return fitness_value != FITNESS_VALUE_INFEASIBLE_SOLUTION  # ...so if a valid solution happens to have fitness equal to this special value, such a solution will be considered infeasible :/
+
+
+def is_feasible_fitness_criteria(fitness_criteria: tuple) -> bool:
+	return all(is_feasible_fitness_value(fitness_value) for fitness_value in fitness_criteria)
+
+
+def select_feasible(individuals):
+	"""
+	Filters out only feasible individuals (i.e., with fitness different from FITNESS_OF_INFEASIBLE_SOLUTION)
+	"""
+	# for ind in individuals:
+	#	print(ind.fitness.values, ind)
+	feasible_individuals = [ind for ind in individuals if is_feasible_fitness_criteria(ind.fitness.values)]
+	count_all = len(individuals)
+	count_infeasible = count_all - len(feasible_individuals)
+	if count_infeasible != 0:
+		print("Selection: ignoring %d infeasible solution%s in a population of size %d" % (count_infeasible, 's' if count_infeasible > 1 else '', count_all))
+	return feasible_individuals
+
+
+def selTournament_only_feasible(individuals, k, tournsize):
+	return tools.selTournament(select_feasible(individuals), k, tournsize=tournsize)
+
+
+def selNSGA2_only_feasible(individuals, k):
+	return tools.selNSGA2(select_feasible(individuals), k)  # this method (unfortunately) decreases population size permanently each time an infeasible solution is removed
 
 
@@ -82,7 +115,9 @@
 	toolbox.register("mutate", frams_mutate, frams_lib)
 	if len(OPTIMIZATION_CRITERIA) <= 1:
-		toolbox.register("select", tools.selTournament, tournsize=tournament_size)
+		# toolbox.register("select", tools.selTournament, tournsize=tournament_size) # without explicitly filtering out infeasible solutions - eliminating/discriminating infeasible solutions during selection would only rely on their relatively poor fitness value
+		toolbox.register("select", selTournament_only_feasible, tournsize=tournament_size)
 	else:
-		toolbox.register("select", tools.selNSGA2)
+		# toolbox.register("select", selNSGA2) # without explicitly filtering out infeasible solutions - eliminating/discriminating infeasible solutions during selection would only rely on their relatively poor fitness value
+		toolbox.register("select", selNSGA2_only_feasible)
 	return toolbox
 
@@ -147,8 +182,10 @@
 	hof = tools.HallOfFame(parsed_args.hof_size)
 	stats = tools.Statistics(lambda ind: ind.fitness.values)
-	stats.register("avg", np.mean)
-	stats.register("stddev", np.std)
-	stats.register("min", np.min)
-	stats.register("max", np.max)
+	# calculate statistics excluding infeasible solutions (by filtering out their fitness=FITNESS_OF_INFEASIBLE_SOLUTION)
+	filter_feasible_for_function = lambda function, fitness_criteria: function(list(filter(is_feasible_fitness_criteria, fitness_criteria)))
+	stats.register("avg", lambda fitness_criteria: filter_feasible_for_function(np.mean, fitness_criteria))
+	stats.register("stddev", lambda fitness_criteria: filter_feasible_for_function(np.std, fitness_criteria))
+	stats.register("min", lambda fitness_criteria: filter_feasible_for_function(np.min, fitness_criteria))
+	stats.register("max", lambda fitness_criteria: filter_feasible_for_function(np.max, fitness_criteria))
 	pop, log = algorithms.eaSimple(pop, toolbox, cxpb=parsed_args.pxov, mutpb=parsed_args.pmut, ngen=parsed_args.generations, stats=stats, halloffame=hof, verbose=True)
 	print('Best individuals:')
