Index: cpp/frams/model/similarity/simil_model.cpp
===================================================================
--- cpp/frams/model/similarity/simil_model.cpp	(revision 818)
+++ cpp/frams/model/similarity/simil_model.cpp	(revision 869)
@@ -1,10 +1,10 @@
 // This file is a part of Framsticks SDK.  http://www.framsticks.com/
-// Copyright (C) 1999-2017  Maciej Komosinski and Szymon Ulatowski.
+// Copyright (C) 1999-2019  Maciej Komosinski and Szymon Ulatowski.
 // See LICENSE.txt for details.
 
-// simil_model.cpp: implementation of the ModelSimil class.
-//
-//////////////////////////////////////////////////////////////////////
+// implementation of the ModelSimil class.
+
 #include "SVD/matrix_tools.h"
+#include "hungarian/hungarian.h"
 #include "simil_model.h"
 #include "simil_match.h"
@@ -35,5 +35,6 @@
 
 static ParamEntry MSparam_tab[] = {
-		{ "Creature: Similarity", 1, 7, "ModelSimilarity", "Evaluates morphological dissimilarity. More information:\nhttp://www.framsticks.com/bib/Komosinski-et-al-2001\nhttp://www.framsticks.com/bib/Komosinski-and-Kubiak-2011\nhttp://www.framsticks.com/bib/Komosinski-2016", },
+		{ "Creature: Similarity", 1, 8, "ModelSimilarity", "Evaluates morphological dissimilarity. More information:\nhttp://www.framsticks.com/bib/Komosinski-et-al-2001\nhttp://www.framsticks.com/bib/Komosinski-and-Kubiak-2011\nhttp://www.framsticks.com/bib/Komosinski-2016\nhttps://doi.org/10.1007/978-3-030-16692-2_8", },
+		{ "simil_method", 0, 0, "Similarity algorithm", "d 0 1 0 ~New (flexible criteria order, optimal matching)~Old (vertex degree order, greedy matching)", FIELD(matching_method), "",},
 		{ "simil_parts", 0, 0, "Weight of parts count", "f 0 100 0", FIELD(m_adFactors[0]), "Differing number of parts is also handled by the 'part degree' similarity component.", },
 		{ "simil_partdeg", 0, 0, "Weight of parts' degree", "f 0 100 1", FIELD(m_adFactors[1]), "", },
@@ -76,17 +77,12 @@
 	isFuzzy = 0;
 	fuzzyDepth = 10;
-	
+
 	//Determines whether weighted MDS should be used.
 	wMDS = 0;
-}
-
-/**	Evaluates distance between two given genotypes. The distance depends strongly
-	on weights set.
-	@param G0 Pointer to the first of compared genotypes
-	@param G1 Pointer to the second of compared genotypes.
-	@return Distance between two genotypes.
-	@sa m_adFactors, matching_method
-	*/
-double ModelSimil::EvaluateDistance(const Geno *G0, const Geno *G1)
+	//Determines whether best matching should be saved using hungarian similarity measure.
+	saveMatching = 0;
+}
+
+double ModelSimil::EvaluateDistanceGreedy(const Geno *G0, const Geno *G1)
 {
 	double dResult = 0;
@@ -98,6 +94,7 @@
 	if (m_Gen[0] == NULL || m_Gen[1] == NULL)
 	{
-		DB(printf("ModelSimil::EvaluateDistance - invalid genotypes pointers\n");)
-			return 0.0;
+		DB(printf("ModelSimil::EvaluateDistanceGreedy - invalid genotype(s) pointers\n");) //TODO convert all such error printfs to legacy error messages, since if's are not in DB(). Example below.
+			logPrintf("ModelSimil", "EvaluateDistanceGreedy", LOG_ERROR, "NULL genotype pointer(s)");
+		return 0.0;
 	}
 	// create models of objects to compare
@@ -108,5 +105,5 @@
 	if (m_Mod[0] == NULL || m_Mod[1] == NULL || !(m_Mod[0]->isValid()) || !(m_Mod[1]->isValid()))
 	{
-		DB(printf("ModelSimil::EvaluateDistance - invalid models pointers\n");)
+		DB(printf("ModelSimil::EvaluateDistanceGreedy - invalid model(s) pointers\n");)
 			return 0.0;
 	}
@@ -143,5 +140,5 @@
 	if ((this->*pfMatchingFunction)() == 0)
 	{
-		DB(printf("ModelSimil::EvaluateDistance - MatchParts() error\n");)
+		DB(printf("ModelSimil::EvaluateDistanceGreedy - MatchParts() error\n");)
 			return 0.0;
 	}
@@ -155,6 +152,6 @@
 		if (CountPartsDistance() == 0)
 		{
-		DB(printf("ModelSimil::EvaluateDistance - CountPartDistance() error\n");)
-			return 0.0;
+			DB(printf("ModelSimil::EvaluateDistanceGreedy - CountPartDistance() error\n");)
+				return 0.0;
 		}
 
@@ -340,6 +337,35 @@
 		if (tdn1[1] < tdn2[1])
 		{
-		// when degree - tdn2[1] - is BIGGER
-		return 1;
+			// when degree - tdn2[1] - is BIGGER
+			return 1;
+		}
+		else
+		{
+			return 0;
+		}
+}
+
+/** Comparison function required for qsort() call. Used while sorting groups of
+	Parts with respect to fuzzy vertex degree. Compares two TDN structures
+	with respect to their [4] field ( fuzzy vertex degree). Highest degree goes first.
+	@param pElem1 Pointer to the TDN structure of some Part.
+	@param pElem2 Pointer to the TDN structure of some Part.
+	@return (-1) - pElem1 should go first, 0 - equal, (1) - pElem2 should go first.
+	*/
+int ModelSimil::CompareFuzzyDegrees(const void *pElem1, const void *pElem2)
+{
+	int *tdn1 = (int *)pElem1;
+	int *tdn2 = (int *)pElem2;
+
+	if (tdn1[4] > tdn2[4])
+	{
+		// when degree - tdn1[4] - is BIGGER
+		return -1;
+	}
+	else
+		if (tdn1[4] < tdn2[4])
+		{
+			// when degree - tdn2[4] - is BIGGER
+			return 1;
 		}
 		else
@@ -374,6 +400,6 @@
 		if (tdn1[NEURO_CONNS] < tdn2[NEURO_CONNS])
 		{
-		// when number of NConn Elem1 is SMALLER
-		return 1;
+			// when number of NConn Elem1 is SMALLER
+			return 1;
 		}
 		else
@@ -389,6 +415,6 @@
 				if (tdn1[NEURONS] < tdn2[NEURONS])
 				{
-				// when number of Neu is SMALLER for Elem1 
-				return 1;
+					// when number of Neu is SMALLER for Elem1 
+					return 1;
 				}
 				else
@@ -406,6 +432,6 @@
 						if (tdn1[ORIG_IND] < tdn2[ORIG_IND])
 						{
-						// when the number of NIt Deg1 id SMALLER
-						return 1;
+							// when the number of NIt Deg1 id SMALLER
+							return 1;
 						}
 						else
@@ -494,4 +520,22 @@
 }
 
+/** Prints one array of parts fuzzy neighbourhood.
+	@param index of organism
+	*/
+void ModelSimil::_PrintFuzzyNeighbourhood(int o)
+{
+	assert(m_fuzzyNeighb[o] != NULL);
+	printf("Fuzzy neighbourhhod of organism %i\n", o);
+	int size = m_Mod[o]->getPartCount();
+	for (int i = 0; i < size; i++)
+	{
+		for (int j = 0; j < fuzzyDepth; j++)
+		{
+			printf("%f ", m_fuzzyNeighb[o][i][j]);
+		}
+		printf("\n");
+	}
+}
+
 /** Creates arrays holding information about organisms' Parts (m_aDegrees) andm_Neigh
 	fills them with initial data (original indices and zeros).
@@ -526,30 +570,30 @@
 				for (j = 0; j < partCount; j++)
 				{
-				m_aDegrees[i][j][0] = j;
-				m_aDegrees[i][j][1] = 0;
-				m_aDegrees[i][j][2] = 0;
-				m_aDegrees[i][j][3] = 0;
-				m_aDegrees[i][j][4] = 0;
-
-				// sprawdz, czy nie piszemy po jakims szalonym miejscu pamieci
-				assert(m_aDegrees[i][j] != NULL);
-
-				if (isFuzzy)
-				{
-					m_Neighbours[i][j] = new int[partCount];
-					for (int k = 0; k < partCount; k++)
+					m_aDegrees[i][j][0] = j;
+					m_aDegrees[i][j][1] = 0;
+					m_aDegrees[i][j][2] = 0;
+					m_aDegrees[i][j][3] = 0;
+					m_aDegrees[i][j][4] = 0;
+
+					// sprawdz, czy nie piszemy po jakims szalonym miejscu pamieci
+					assert(m_aDegrees[i][j] != NULL);
+
+					if (isFuzzy)
 					{
-						m_Neighbours[i][j][k] = 0;
+						m_Neighbours[i][j] = new int[partCount];
+						for (int k = 0; k < partCount; k++)
+						{
+							m_Neighbours[i][j][k] = 0;
+						}
+
+						m_fuzzyNeighb[i][j] = new float[fuzzyDepth];
+						for (int k = 0; k < fuzzyDepth; k++)
+						{
+							m_fuzzyNeighb[i][j][k] = 0;
+						}
+
+						assert(m_Neighbours[i][j] != NULL);
+						assert(m_fuzzyNeighb[i][j] != NULL);
 					}
-
-					m_fuzzyNeighb[i][j] = new float[fuzzyDepth];
-					for (int k = 0; k < fuzzyDepth; k++)
-					{
-						m_fuzzyNeighb[i][j][k] = 0;
-					}
-
-					assert(m_Neighbours[i][j] != NULL);
-					assert(m_fuzzyNeighb[i][j] != NULL);
-				}
 
 				}
@@ -671,21 +715,23 @@
 }
 
-//store information about identity of parts "fuzzy degrees" in the m_aDegrees[4]
-void ModelSimil::CheckFuzzyIdentity()
-{
-	int partCount = 0;
+void ModelSimil::FuzzyOrder()
+{
+	int i, depth, partInd, prevPartInd, partCount;
 	for (int mod = 0; mod < 2; mod++)
 	{
-		//for subsequent pairs of parts
 		partCount = m_Mod[mod]->getPartCount();
-		m_aDegrees[mod][partCount - 1][FUZZ_DEG] = 0;
-		for (int partInd = (partCount - 2); partInd >= 0; partInd--)
-		{
-			m_aDegrees[mod][partInd][FUZZ_DEG] = m_aDegrees[mod][partInd + 1][FUZZ_DEG];
-			for (int depth = 1; depth < fuzzyDepth; depth++)
-			{
-				if (m_fuzzyNeighb[mod][partInd][depth] != m_fuzzyNeighb[mod][partInd + 1][depth])
-				{
-					m_aDegrees[mod][partInd][FUZZ_DEG] += 1;
+		partInd = m_fuzzyNeighb[mod][partCount - 1][0];
+		m_aDegrees[mod][partInd][FUZZ_DEG] = 0;
+
+		for (i = (partCount - 2); i >= 0; i--)
+		{
+			prevPartInd = partInd;
+			partInd = m_fuzzyNeighb[mod][i][0];
+			m_aDegrees[mod][partInd][FUZZ_DEG] = m_aDegrees[mod][prevPartInd][FUZZ_DEG];
+			for (depth = 1; depth < fuzzyDepth; depth++)
+			{
+				if (m_fuzzyNeighb[mod][i][depth] != m_fuzzyNeighb[mod][i + 1][depth])
+				{
+					m_aDegrees[mod][partInd][FUZZ_DEG]++;
 					break;
 				}
@@ -761,4 +807,5 @@
 
 	SortFuzzyNeighb();
+	FuzzyOrder();
 }
 
@@ -884,21 +931,14 @@
 
 	int i;
+	int(*pfDegreeFunction) (const void*, const void*) = NULL;
+	pfDegreeFunction = (isFuzzy == 1) ? &CompareFuzzyDegrees : &CompareDegrees;
 	// sortowanie obu tablic wg stopni punktów - TDN[1]
-	if (isFuzzy != 1)
-	{
-		for (i = 0; i < 2; i++)
-		{
-			DB(_PrintDegrees(i));
-			std::qsort(m_aDegrees[i], (size_t)(m_Mod[i]->getPartCount()),
-				sizeof(TDN), ModelSimil::CompareDegrees);
-			DB(_PrintDegrees(i));
-		}
-	}//sortowanie wg romzmytego stopnia wierzcholka
-
-	else
-	{
-		SortPartInfoFuzzy();
-	}
-
+	for (i = 0; i < 2; i++)
+	{
+		DB(_PrintDegrees(i));
+		std::qsort(m_aDegrees[i], (size_t)(m_Mod[i]->getPartCount()),
+			sizeof(TDN), pfDegreeFunction);
+		DB(_PrintDegrees(i));
+	}
 
 	// sprawdzenie wartosci parametru
@@ -935,5 +975,5 @@
 
 					std::qsort(m_aDegrees[i][iPocz], (size_t)(j - iPocz),
-					sizeof(TDN), ModelSimil::CompareConnsNo);
+						sizeof(TDN), ModelSimil::CompareConnsNo);
 				DB(_PrintDegrees(i));
 				// wyswietlamy z jedna komorka po zakonczeniu tablicy
@@ -948,81 +988,4 @@
 }
 
-int ModelSimil::SortPartInfoFuzzy()
-{
-
-	// sprawdz zalozenie - o modelach
-	assert((m_Mod[0] != NULL) && (m_Mod[1] != NULL));
-	assert(m_Mod[0]->isValid() && m_Mod[1]->isValid());
-
-	// sprawdz zalozenie - o tablicach
-	assert(m_aDegrees[0] != NULL);
-	assert(m_aDegrees[1] != NULL);
-	// sprawdz zalozenie - o tablicach
-	assert(m_fuzzyNeighb[0] != NULL);
-	assert(m_fuzzyNeighb[1] != NULL);
-
-
-	TDN * m_aDegreesTmp[2];
-
-	for (int i = 0; i < 2; i++)
-	{
-		int partCount = m_Mod[i]->getPartCount();
-		m_aDegreesTmp[i] = new TDN[partCount];
-
-		for (int j = 0; j < partCount; j++)
-		{
-			for (int k = 0; k < TDN_SIZE; k++)
-			{
-				m_aDegreesTmp[i][j][k] = m_aDegrees[i][j][k];
-			}
-		}
-	}
-
-	int newInd = 0;
-	for (int i = 0; i < 2; i++)
-	{
-		int partCount = m_Mod[i]->getPartCount();
-		for (int j = 0; j < partCount; j++)
-		{
-			newInd = (int)m_fuzzyNeighb[i][j][0];
-			for (int k = 0; k < TDN_SIZE; k++)
-			{
-				m_aDegrees[i][j][k] = m_aDegreesTmp[i][newInd][k];
-			}
-		}
-	}
-
-	SAFEDELETEARRAY(m_aDegreesTmp[0]);
-	SAFEDELETEARRAY(m_aDegreesTmp[1]);
-
-	CheckFuzzyIdentity();
-
-	return 1;
-}
-
-/** Checks if given Parts have identical physical and biological properties
-	(except for geometry that might differ).
-	@param P1 Pointer to first Part.
-	@param P2 Pointer to second Part.
-	@return 1 - identical properties, 0 - there are differences.
-	*/
-int ModelSimil::CheckPartsIdentity(Part *P1, Part *P2)
-{
-	// sprawdz, czy te Parts chociaz sa w sensownym miejscu pamieci
-	assert((P1 != NULL) && (P2 != NULL));
-
-	if ((P1->assim != P2->assim) ||
-		(P1->friction != P2->friction) ||
-		(P1->ingest != P2->ingest) ||
-		(P1->mass != P2->mass) ||
-		(P1->size != P2->size) ||
-		(P1->density != P2->density))
-		// gdy znaleziono jakas roznice w parametrach fizycznych i
-		// biologicznych
-		return 0;
-	else
-		// gdy nie ma roznic
-		return 1;
-}
 
 /** Prints the state of the matching object. Debug method.
@@ -1158,62 +1121,62 @@
 			aiKoniecGrupyDopasowania[0]);)
 			DB(printf("Organizm 1: grupa: (%2i, %2i)\n", aiKoniecPierwszejGrupy[1],
-			aiKoniecGrupyDopasowania[1]);)
+				aiKoniecGrupyDopasowania[1]);)
 
 			for (i = 0; i < aiRozmiarCalychGrup[0]; i++)
 			{
 
-			// iterujemy i - Parts organizmu 0
-			// (indeks podstawowy - aiKoniecPierwszejGrupy[0])
-
-			if (!(m_pMatching->IsMatched(0, aiKoniecPierwszejGrupy[0] + i)))
-			{
-				// interesuja nas tylko te niedopasowane jeszcze (i)
-				for (j = 0; j < aiRozmiarCalychGrup[1]; j++)
-				{
-
-					// iterujemy j - Parts organizmu 1
-					// (indeks podstawowy - aiKoniecPierwszejGrupy[1])
-
-					if (!(m_pMatching->IsMatched(1, aiKoniecPierwszejGrupy[1] + j)))
+				// iterujemy i - Parts organizmu 0
+				// (indeks podstawowy - aiKoniecPierwszejGrupy[0])
+
+				if (!(m_pMatching->IsMatched(0, aiKoniecPierwszejGrupy[0] + i)))
+				{
+					// interesuja nas tylko te niedopasowane jeszcze (i)
+					for (j = 0; j < aiRozmiarCalychGrup[1]; j++)
 					{
-						// interesuja nas tylko te niedopasowane jeszcze (j)
-						// teraz obliczymy lokalne różnice pomiędzy Parts
-						iDeg = abs(m_aDegrees[0][aiKoniecPierwszejGrupy[0] + i][1]
-							- m_aDegrees[1][aiKoniecPierwszejGrupy[1] + j][1]);
-
-						//iNit currently is not a component of distance measure            
-						//iNIt = abs(m_aDegrees[0][ aiKoniecPierwszejGrupy[0] + i ][2]
-						//           - m_aDegrees[1][ aiKoniecPierwszejGrupy[1] + j ][2]);
-
-						iNeu = abs(m_aDegrees[0][aiKoniecPierwszejGrupy[0] + i][3]
-							- m_aDegrees[1][aiKoniecPierwszejGrupy[1] + j][3]);
-
-						// obliczenie także lokalnych różnic geometrycznych pomiędzy Parts
-						// find original indices of Parts for both of the models
-						iPartIndex[0] = m_aDegrees[0][aiKoniecPierwszejGrupy[0] + i][0];
-						iPartIndex[1] = m_aDegrees[1][aiKoniecPierwszejGrupy[1] + j][0];
-						// now compute the geometrical distance of these Parts (use m_aPositions
-						// which should be computed by SVD)
-						Pt3D Part0Pos(m_aPositions[0][iPartIndex[0]]);
-						Pt3D Part1Pos(m_aPositions[1][iPartIndex[1]]);
-						dGeo = m_adFactors[3] == 0 ? 0 : Part0Pos.distanceTo(Part1Pos); //no need to compute distane when m_dDG weight is 0
-
-						// tutaj prawdopodobnie należy jeszcze dodać sprawdzanie
-						// identyczności pozostałych własności (oprócz geometrii)
-						// - żeby móc stwierdzić w ogóle identyczność Parts
-
-						// i ostateczna odleglosc indukowana przez te roznice
-						// (tutaj nie ma różnicy w liczbie wszystkich wierzchołków)
-						aadOdleglosciParts[i][j] = m_adFactors[1] * double(iDeg) +
-							m_adFactors[2] * double(iNeu) +
-							m_adFactors[3] * dGeo;
-						DB(printf("Parts Distance (%2i,%2i) = %.3lf\n", aiKoniecPierwszejGrupy[0] + i,
-							aiKoniecPierwszejGrupy[1] + j, aadOdleglosciParts[i][j]);)
-							DB(printf("Parts geometrical distance = %.20lf\n", dGeo);)
-							DB(printf("Pos0: (%.3lf %.3lf %.3lf)\n", Part0Pos.x, Part0Pos.y, Part0Pos.z);)
-							DB(printf("Pos1: (%.3lf %.3lf %.3lf)\n", Part1Pos.x, Part1Pos.y, Part1Pos.z);)
+
+						// iterujemy j - Parts organizmu 1
+						// (indeks podstawowy - aiKoniecPierwszejGrupy[1])
+
+						if (!(m_pMatching->IsMatched(1, aiKoniecPierwszejGrupy[1] + j)))
+						{
+							// interesuja nas tylko te niedopasowane jeszcze (j)
+							// teraz obliczymy lokalne różnice pomiędzy Parts
+							iDeg = abs(m_aDegrees[0][aiKoniecPierwszejGrupy[0] + i][1]
+								- m_aDegrees[1][aiKoniecPierwszejGrupy[1] + j][1]);
+
+							//iNit currently is not a component of distance measure            
+							//iNIt = abs(m_aDegrees[0][ aiKoniecPierwszejGrupy[0] + i ][2]
+							//           - m_aDegrees[1][ aiKoniecPierwszejGrupy[1] + j ][2]);
+
+							iNeu = abs(m_aDegrees[0][aiKoniecPierwszejGrupy[0] + i][3]
+								- m_aDegrees[1][aiKoniecPierwszejGrupy[1] + j][3]);
+
+							// obliczenie także lokalnych różnic geometrycznych pomiędzy Parts
+							// find original indices of Parts for both of the models
+							iPartIndex[0] = m_aDegrees[0][aiKoniecPierwszejGrupy[0] + i][0];
+							iPartIndex[1] = m_aDegrees[1][aiKoniecPierwszejGrupy[1] + j][0];
+							// now compute the geometrical distance of these Parts (use m_aPositions
+							// which should be computed by SVD)
+							Pt3D Part0Pos(m_aPositions[0][iPartIndex[0]]);
+							Pt3D Part1Pos(m_aPositions[1][iPartIndex[1]]);
+							dGeo = m_adFactors[3] == 0 ? 0 : Part0Pos.distanceTo(Part1Pos); //no need to compute distane when m_dDG weight is 0
+
+							// tutaj prawdopodobnie należy jeszcze dodać sprawdzanie
+							// identyczności pozostałych własności (oprócz geometrii)
+							// - żeby móc stwierdzić w ogóle identyczność Parts
+
+							// i ostateczna odleglosc indukowana przez te roznice
+							// (tutaj nie ma różnicy w liczbie wszystkich wierzchołków)
+							aadOdleglosciParts[i][j] = m_adFactors[1] * double(iDeg) +
+								m_adFactors[2] * double(iNeu) +
+								m_adFactors[3] * dGeo;
+							DB(printf("Parts Distance (%2i,%2i) = %.3lf\n", aiKoniecPierwszejGrupy[0] + i,
+								aiKoniecPierwszejGrupy[1] + j, aadOdleglosciParts[i][j]);)
+								DB(printf("Parts geometrical distance = %.20lf\n", dGeo);)
+								DB(printf("Pos0: (%.3lf %.3lf %.3lf)\n", Part0Pos.x, Part0Pos.y, Part0Pos.z);)
+								DB(printf("Pos1: (%.3lf %.3lf %.3lf)\n", Part1Pos.x, Part1Pos.y, Part1Pos.z);)
+						}
 					}
 				}
-			}
 			}
 
@@ -1388,150 +1351,31 @@
 					if (PartZ1NaLiscie0 || PartZ0NaLiscie1)
 					{
-					// PRZYPADEK 2. Tylko jeden z Parts ma drugiego na swojej liscie
-					// mozliwych dopasowan
-					// AKCJA. Dopasowanie jednego jest proste (tego, ktory nie ma 
-					// na swojej liscie drugiego). Dla tego drugiego nalezy powtorzyc
-					// duza czesc obliczen (znalezc mu nowa mozliwa pare)
-
-					// indeks organizmu, ktorego Part nie ma dopasowywanego Part
-					// z drugiego organizmu na swojej liscie
-					int iBezDrugiego;
-
-					// okreslenie indeksu organizmu z dopasowywanym Part
-					if (!PartZ1NaLiscie0)
-					{
-						iBezDrugiego = 0;
-					}
-					else
-					{
-						iBezDrugiego = 1;
-					}
-					// sprawdz indeks organizmu
-					assert((iBezDrugiego == 0) || (iBezDrugiego == 1));
-
-					int iDopasowywany = -1;
-					// poszukujemy pierwszego z listy dopasowania
-					for (i = 0; i < aiRozmiarCalychGrup[1 - iBezDrugiego]; i++)
-					{
-						if (apvbCzyMinimalnaOdleglosc[iBezDrugiego]->operator[](i))
+						// PRZYPADEK 2. Tylko jeden z Parts ma drugiego na swojej liscie
+						// mozliwych dopasowan
+						// AKCJA. Dopasowanie jednego jest proste (tego, ktory nie ma 
+						// na swojej liscie drugiego). Dla tego drugiego nalezy powtorzyc
+						// duza czesc obliczen (znalezc mu nowa mozliwa pare)
+
+						// indeks organizmu, ktorego Part nie ma dopasowywanego Part
+						// z drugiego organizmu na swojej liscie
+						int iBezDrugiego;
+
+						// okreslenie indeksu organizmu z dopasowywanym Part
+						if (!PartZ1NaLiscie0)
 						{
-							iDopasowywany = i;
-							break;
+							iBezDrugiego = 0;
 						}
-					}
-					// sprawdz poprawnosc indeksu dopasowywanego (musimy cos znalezc!)
-					// nieujemny i w odpowiedniej grupie!
-					assert((iDopasowywany >= 0) &&
-						(iDopasowywany < aiRozmiarCalychGrup[1 - iBezDrugiego]));
-
-					// znalezlismy 1. Part z listy dopasowania - dopasowujemy!
-					m_pMatching->Match(
-						iBezDrugiego,
-						aiKoniecPierwszejGrupy[iBezDrugiego] + iIndex[iBezDrugiego],
-						1 - iBezDrugiego,
-						aiKoniecPierwszejGrupy[1 - iBezDrugiego] + iDopasowywany);
-					DB(printf("Przypadek 2.1.: dopasowane Parts dopasowanie bez drugiego: (%2i, %2i)\n", aiKoniecPierwszejGrupy[iBezDrugiego] + iIndex[iBezDrugiego],
-						aiKoniecPierwszejGrupy[1 - iBezDrugiego] + iDopasowywany);)
-
-						// zmniejsz liczby niedopasowanych elementow w grupach
-						aiRozmiarGrupy[0]--;
-					aiRozmiarGrupy[1]--;
-
-					// sprawdz, czy jest szansa dopasowania tego Part z drugiej strony
-					// (ktora miala mozliwosc dopasowania tego Part z poprzedniego organizmu)
-					if ((aiRozmiarGrupy[0] > 0) && (aiRozmiarGrupy[1] > 0))
-					{
-						// jesli grupy sie jeszcze nie wyczrpaly
-						// to jest mozliwosc dopasowania w organizmie
-
-						int iZDrugim = 1 - iBezDrugiego;
+						else
+						{
+							iBezDrugiego = 1;
+						}
 						// sprawdz indeks organizmu
-						assert((iZDrugim == 0) || (iZDrugim == 1));
-
-						// bedziemy szukac minimum wsrod niedopasowanych - musimy wykasowac
-						// poprzednie obliczenia minimum
-						// dla organizmu 1 (o rozmiarze grupy z 0)
-						for (i = 0; i < aiRozmiarCalychGrup[0]; i++)
+						assert((iBezDrugiego == 0) || (iBezDrugiego == 1));
+
+						int iDopasowywany = -1;
+						// poszukujemy pierwszego z listy dopasowania
+						for (i = 0; i < aiRozmiarCalychGrup[1 - iBezDrugiego]; i++)
 						{
-							apvbCzyMinimalnaOdleglosc[1]->operator[](i) = false;
-						}
-						// dla organizmu 0 (o rozmiarze grup z 1)
-						for (i = 0; i < aiRozmiarCalychGrup[1]; i++)
-						{
-							apvbCzyMinimalnaOdleglosc[0]->operator[](i) = false;
-						}
-
-						// szukamy na nowo minimum dla Part o indeksie iIndex[ iZDrugim ] w organizmie iZDrugim
-						// wsrod niedopasowanych Parts z organizmu 1 - iZDrugim
-						dMinimum = HUGE_VAL;
-						for (i = 0; i < aiRozmiarCalychGrup[1 - iZDrugim]; i++)
-						{
-							if (!(m_pMatching->IsMatched(
-								1 - iZDrugim,
-								aiKoniecPierwszejGrupy[1 - iZDrugim] + i)))
-							{
-								// szukamy minimum obliczonej lokalnej odleglosci tylko wsrod
-								// niedopasowanych jeszcze Parts
-								if (iZDrugim == 0)
-								{
-									// teraz niestety musimy rozpoznac odpowiedni organizm
-									// zeby moc indeksowac niesymetryczna tablice
-									if (aadOdleglosciParts[iIndex[0]][i] < dMinimum)
-									{
-										dMinimum = aadOdleglosciParts[iIndex[0]][i];
-									}
-									// przy okazji - sprawdz, czy HUGE_VAL jest rzeczywiscie max dla double
-									assert(aadOdleglosciParts[iIndex[0]][i] < HUGE_VAL);
-
-								}
-								else
-								{
-									if (aadOdleglosciParts[i][iIndex[1]] < dMinimum)
-									{
-										dMinimum = aadOdleglosciParts[i][iIndex[1]];
-									}
-									// przy okazji - sprawdz, czy HUGE_VAL jest rzeczywiscie max dla double
-									assert(aadOdleglosciParts[i][iIndex[1]] < HUGE_VAL);
-								}
-							}
-						}
-						// sprawdz, czy minimum znaleziono - musi takie byc, bo jest cos niedopasowanego
-						assert((dMinimum >= 0.0) && (dMinimum < HUGE_VAL));
-
-						// teraz zaznaczamy w tablicy te wszystkie Parts, ktore maja
-						// rzeczywiscie te minimalne odleglosci do Part w organizmie iZDrugim
-						for (i = 0; i < aiRozmiarCalychGrup[1 - iZDrugim]; i++)
-						{
-							if (!(m_pMatching->IsMatched(
-								1 - iZDrugim,
-								aiKoniecPierwszejGrupy[1 - iZDrugim] + i)))
-							{
-								if (iZDrugim == 0)
-								{
-									// teraz niestety musimy rozpoznac odpowiedni organizm
-									// zeby moc indeksowac niesymetryczna tablice
-									if (aadOdleglosciParts[iIndex[0]][i] == dMinimum)
-									{
-										// jesli w danym miejscu tablicy odleglosci jest faktyczne 
-										// minimum odleglosci, to mamy potencjalna pare dla iIndex[1]
-										apvbCzyMinimalnaOdleglosc[0]->operator[](i) = true;
-									}
-								}
-								else
-								{
-									if (aadOdleglosciParts[i][iIndex[1]] == dMinimum)
-									{
-										apvbCzyMinimalnaOdleglosc[1]->operator[](i) = true;
-									}
-								}
-							}
-						}
-
-						// a teraz - po znalezieniu potencjalnych elementow do dopasowania
-						// dopasowujemy pierwszy z potencjalnych
-						iDopasowywany = -1;
-						for (i = 0; i < aiRozmiarCalychGrup[1 - iZDrugim]; i++)
-						{
-							if (apvbCzyMinimalnaOdleglosc[iZDrugim]->operator[](i))
+							if (apvbCzyMinimalnaOdleglosc[iBezDrugiego]->operator[](i))
 							{
 								iDopasowywany = i;
@@ -1539,31 +1383,150 @@
 							}
 						}
-						// musielismy znalezc jakiegos dopasowywanego
+						// sprawdz poprawnosc indeksu dopasowywanego (musimy cos znalezc!)
+						// nieujemny i w odpowiedniej grupie!
 						assert((iDopasowywany >= 0) &&
-							(iDopasowywany < aiRozmiarCalychGrup[1 - iZDrugim]));
-
-						// no to juz mozemy dopasowac
+							(iDopasowywany < aiRozmiarCalychGrup[1 - iBezDrugiego]));
+
+						// znalezlismy 1. Part z listy dopasowania - dopasowujemy!
 						m_pMatching->Match(
-							iZDrugim,
-							aiKoniecPierwszejGrupy[iZDrugim] + iIndex[iZDrugim],
-							1 - iZDrugim,
-							aiKoniecPierwszejGrupy[1 - iZDrugim] + iDopasowywany);
-						DB(printf("Przypadek 2.1.: dopasowane Parts dopasowaniebz drugim: (%2i, %2i)\n", aiKoniecPierwszejGrupy[iZDrugim] + iIndex[iZDrugim], aiKoniecPierwszejGrupy[1 - iZDrugim] + iDopasowywany);)
-
-							//aiKoniecPierwszejGrupy[ 1-iBezDrugiego ] + iDopasowywany ;)
-							//aiKoniecPierwszejGrupy[ 1-iBezDrugiego ] + iDopasowywany ;)
+							iBezDrugiego,
+							aiKoniecPierwszejGrupy[iBezDrugiego] + iIndex[iBezDrugiego],
+							1 - iBezDrugiego,
+							aiKoniecPierwszejGrupy[1 - iBezDrugiego] + iDopasowywany);
+						DB(printf("Przypadek 2.1.: dopasowane Parts dopasowanie bez drugiego: (%2i, %2i)\n", aiKoniecPierwszejGrupy[iBezDrugiego] + iIndex[iBezDrugiego],
+							aiKoniecPierwszejGrupy[1 - iBezDrugiego] + iDopasowywany);)
+
 							// zmniejsz liczby niedopasowanych elementow w grupach
 							aiRozmiarGrupy[0]--;
 						aiRozmiarGrupy[1]--;
 
-					}
-					else
-					{
-						// jedna z grup sie juz wyczerpala
-						// wiec nie ma mozliwosci dopasowania tego drugiego Partu
-						/// i trzeba poczekac na zmiane grupy
-					}
-
-					DB(printf("Przypadek 2.\n");)
+						// sprawdz, czy jest szansa dopasowania tego Part z drugiej strony
+						// (ktora miala mozliwosc dopasowania tego Part z poprzedniego organizmu)
+						if ((aiRozmiarGrupy[0] > 0) && (aiRozmiarGrupy[1] > 0))
+						{
+							// jesli grupy sie jeszcze nie wyczrpaly
+							// to jest mozliwosc dopasowania w organizmie
+
+							int iZDrugim = 1 - iBezDrugiego;
+							// sprawdz indeks organizmu
+							assert((iZDrugim == 0) || (iZDrugim == 1));
+
+							// bedziemy szukac minimum wsrod niedopasowanych - musimy wykasowac
+							// poprzednie obliczenia minimum
+							// dla organizmu 1 (o rozmiarze grupy z 0)
+							for (i = 0; i < aiRozmiarCalychGrup[0]; i++)
+							{
+								apvbCzyMinimalnaOdleglosc[1]->operator[](i) = false;
+							}
+							// dla organizmu 0 (o rozmiarze grup z 1)
+							for (i = 0; i < aiRozmiarCalychGrup[1]; i++)
+							{
+								apvbCzyMinimalnaOdleglosc[0]->operator[](i) = false;
+							}
+
+							// szukamy na nowo minimum dla Part o indeksie iIndex[ iZDrugim ] w organizmie iZDrugim
+							// wsrod niedopasowanych Parts z organizmu 1 - iZDrugim
+							dMinimum = HUGE_VAL;
+							for (i = 0; i < aiRozmiarCalychGrup[1 - iZDrugim]; i++)
+							{
+								if (!(m_pMatching->IsMatched(
+									1 - iZDrugim,
+									aiKoniecPierwszejGrupy[1 - iZDrugim] + i)))
+								{
+									// szukamy minimum obliczonej lokalnej odleglosci tylko wsrod
+									// niedopasowanych jeszcze Parts
+									if (iZDrugim == 0)
+									{
+										// teraz niestety musimy rozpoznac odpowiedni organizm
+										// zeby moc indeksowac niesymetryczna tablice
+										if (aadOdleglosciParts[iIndex[0]][i] < dMinimum)
+										{
+											dMinimum = aadOdleglosciParts[iIndex[0]][i];
+										}
+										// przy okazji - sprawdz, czy HUGE_VAL jest rzeczywiscie max dla double
+										assert(aadOdleglosciParts[iIndex[0]][i] < HUGE_VAL);
+
+									}
+									else
+									{
+										if (aadOdleglosciParts[i][iIndex[1]] < dMinimum)
+										{
+											dMinimum = aadOdleglosciParts[i][iIndex[1]];
+										}
+										// przy okazji - sprawdz, czy HUGE_VAL jest rzeczywiscie max dla double
+										assert(aadOdleglosciParts[i][iIndex[1]] < HUGE_VAL);
+									}
+								}
+							}
+							// sprawdz, czy minimum znaleziono - musi takie byc, bo jest cos niedopasowanego
+							assert((dMinimum >= 0.0) && (dMinimum < HUGE_VAL));
+
+							// teraz zaznaczamy w tablicy te wszystkie Parts, ktore maja
+							// rzeczywiscie te minimalne odleglosci do Part w organizmie iZDrugim
+							for (i = 0; i < aiRozmiarCalychGrup[1 - iZDrugim]; i++)
+							{
+								if (!(m_pMatching->IsMatched(
+									1 - iZDrugim,
+									aiKoniecPierwszejGrupy[1 - iZDrugim] + i)))
+								{
+									if (iZDrugim == 0)
+									{
+										// teraz niestety musimy rozpoznac odpowiedni organizm
+										// zeby moc indeksowac niesymetryczna tablice
+										if (aadOdleglosciParts[iIndex[0]][i] == dMinimum)
+										{
+											// jesli w danym miejscu tablicy odleglosci jest faktyczne 
+											// minimum odleglosci, to mamy potencjalna pare dla iIndex[1]
+											apvbCzyMinimalnaOdleglosc[0]->operator[](i) = true;
+										}
+									}
+									else
+									{
+										if (aadOdleglosciParts[i][iIndex[1]] == dMinimum)
+										{
+											apvbCzyMinimalnaOdleglosc[1]->operator[](i) = true;
+										}
+									}
+								}
+							}
+
+							// a teraz - po znalezieniu potencjalnych elementow do dopasowania
+							// dopasowujemy pierwszy z potencjalnych
+							iDopasowywany = -1;
+							for (i = 0; i < aiRozmiarCalychGrup[1 - iZDrugim]; i++)
+							{
+								if (apvbCzyMinimalnaOdleglosc[iZDrugim]->operator[](i))
+								{
+									iDopasowywany = i;
+									break;
+								}
+							}
+							// musielismy znalezc jakiegos dopasowywanego
+							assert((iDopasowywany >= 0) &&
+								(iDopasowywany < aiRozmiarCalychGrup[1 - iZDrugim]));
+
+							// no to juz mozemy dopasowac
+							m_pMatching->Match(
+								iZDrugim,
+								aiKoniecPierwszejGrupy[iZDrugim] + iIndex[iZDrugim],
+								1 - iZDrugim,
+								aiKoniecPierwszejGrupy[1 - iZDrugim] + iDopasowywany);
+							DB(printf("Przypadek 2.1.: dopasowane Parts dopasowaniebz drugim: (%2i, %2i)\n", aiKoniecPierwszejGrupy[iZDrugim] + iIndex[iZDrugim], aiKoniecPierwszejGrupy[1 - iZDrugim] + iDopasowywany);)
+
+								//aiKoniecPierwszejGrupy[ 1-iBezDrugiego ] + iDopasowywany ;)
+								//aiKoniecPierwszejGrupy[ 1-iBezDrugiego ] + iDopasowywany ;)
+								// zmniejsz liczby niedopasowanych elementow w grupach
+								aiRozmiarGrupy[0]--;
+							aiRozmiarGrupy[1]--;
+
+						}
+						else
+						{
+							// jedna z grup sie juz wyczerpala
+							// wiec nie ma mozliwosci dopasowania tego drugiego Partu
+							/// i trzeba poczekac na zmiane grupy
+						}
+
+						DB(printf("Przypadek 2.\n");)
 
 					}// PRZYPADEK 2.
@@ -1724,121 +1687,121 @@
 		if (m_adFactors[3] > 0)
 		{
-		// tutaj zacznij pętlę po przekształceniach  geometrycznych
-		const int NO_OF_TRANSFORM = 8; // liczba transformacji geometrycznych (na razie tylko ID i O_YZ)
-		// tablice transformacji współrzędnych; nie są to dokładnie tablice tranformacji, ale raczej tablice PRZEJŚĆ
-		// pomiędzy transformacjami; 
-		// wartości orginalne transformacji dOrig uzyskuje się przez:
-		// for ( iTrans = 0; iTrans <= TRANS_INDEX; iTrans++ ) dOrig *= dMul[ iTrans ];
-		//const char *szTransformNames[NO_OF_TRANSFORM] = { "ID", "S_yz", "S_xz", "S_xy", "R180_z", "R180_y", "R180_z", "S_(0,0,0)" };
-		const int dMulX[NO_OF_TRANSFORM] = { 1, -1, -1, 1, -1, 1, -1, -1 };
-		const int dMulY[NO_OF_TRANSFORM] = { 1, 1, -1, -1, -1, -1, -1, 1 };
-		const int dMulZ[NO_OF_TRANSFORM] = { 1, 1, 1, -1, -1, -1, 1, 1 };
+			// tutaj zacznij pętlę po przekształceniach  geometrycznych
+			const int NO_OF_TRANSFORM = 8; // liczba transformacji geometrycznych (na razie tylko ID i O_YZ)
+			// tablice transformacji współrzędnych; nie są to dokładnie tablice tranformacji, ale raczej tablice PRZEJŚĆ
+			// pomiędzy transformacjami; 
+			// wartości orginalne transformacji dOrig uzyskuje się przez:
+			// for ( iTrans = 0; iTrans <= TRANS_INDEX; iTrans++ ) dOrig *= dMul[ iTrans ];
+			//const char *szTransformNames[NO_OF_TRANSFORM] = { "ID", "S_yz", "S_xz", "S_xy", "R180_z", "R180_y", "R180_z", "S_(0,0,0)" };
+			const int dMulX[NO_OF_TRANSFORM] = { 1, -1, -1, 1, -1, 1, -1, -1 };
+			const int dMulY[NO_OF_TRANSFORM] = { 1, 1, -1, -1, -1, -1, -1, 1 };
+			const int dMulZ[NO_OF_TRANSFORM] = { 1, 1, 1, -1, -1, -1, 1, 1 };
 
 #ifdef max
 #undef max //this macro would conflict with line below
 #endif
-		double dMinSimValue = std::numeric_limits<double>::max(); // minimum value of similarity
-		int iMinSimTransform = -1; // index of the transformation with the minimum similarity
-		SimilMatching *pMinSimMatching = NULL; // matching with the minimum value of similarity
-
-		// remember the original positions of model 0 as computed by SVD in order to restore them later, after 
-		// all transformations have been computed
-		Pt3D *StoredPositions = NULL; // array of positions of Parts, for one (0th) model
-		// create the stored positions
-		StoredPositions = new Pt3D[m_Mod[0]->getPartCount()];
-		assert(StoredPositions != NULL);
-		// copy the original positions of model 0 (store them)
-		int iPart; // a counter of Parts
-		for (iPart = 0; iPart < m_Mod[0]->getPartCount(); iPart++)
-		{
-			StoredPositions[iPart].x = m_aPositions[0][iPart].x;
-			StoredPositions[iPart].y = m_aPositions[0][iPart].y;
-			StoredPositions[iPart].z = m_aPositions[0][iPart].z;
-		}
-		// now the original positions of model 0 are stored
-
-
-		int iTransform; // a counter of geometric transformations
-		for (iTransform = 0; iTransform < NO_OF_TRANSFORM; iTransform++)
-		{
-			// for each geometric transformation to be done
-			// entry conditions:
-			// - models (m_Mod) exist and are available
-			// - matching (m_pMatching) exists and is empty
-			// - all properties are created and available (m_aDegrees and m_aPositions)
-
-			// recompute geometric properties according to the transformation iTransform
-			// but only for model 0
+			double dMinSimValue = std::numeric_limits<double>::max(); // minimum value of similarity
+			int iMinSimTransform = -1; // index of the transformation with the minimum similarity
+			SimilMatching *pMinSimMatching = NULL; // matching with the minimum value of similarity
+
+			// remember the original positions of model 0 as computed by SVD in order to restore them later, after 
+			// all transformations have been computed
+			Pt3D *StoredPositions = NULL; // array of positions of Parts, for one (0th) model
+			// create the stored positions
+			StoredPositions = new Pt3D[m_Mod[0]->getPartCount()];
+			assert(StoredPositions != NULL);
+			// copy the original positions of model 0 (store them)
+			int iPart; // a counter of Parts
 			for (iPart = 0; iPart < m_Mod[0]->getPartCount(); iPart++)
 			{
-				// for each iPart, a part of the model iMod
-				m_aPositions[0][iPart].x *= dMulX[iTransform];
-				m_aPositions[0][iPart].y *= dMulY[iTransform];
-				m_aPositions[0][iPart].z *= dMulZ[iTransform];
+				StoredPositions[iPart].x = m_aPositions[0][iPart].x;
+				StoredPositions[iPart].y = m_aPositions[0][iPart].y;
+				StoredPositions[iPart].z = m_aPositions[0][iPart].z;
 			}
-			// now the positions are recomputed
-			ComputeMatching();
-
-			// teraz m_pMatching istnieje i jest pełne
-			assert(m_pMatching != NULL);
-			assert(m_pMatching->IsFull() == true);
-
-			// wykorzystaj to dopasowanie i geometrię do obliczenia tymczasowej wartości miary
-			int iTempRes = CountPartsDistance();
-			// załóż sukces
-			assert(iTempRes == 1);
-			double dCurrentSim = m_adFactors[0] * double(m_iDV) +
-				m_adFactors[1] * double(m_iDD) +
-				m_adFactors[2] * double(m_iDN) +
-				m_adFactors[3] * double(m_dDG);
-			// załóż poprawną wartość podobieństwa
-			assert(dCurrentSim >= 0.0);
-
-			// porównaj wartość obliczoną z dotychczasowym minimum
-			if (dCurrentSim < dMinSimValue)
-			{
-				// jeśli uzyskano mniejszą wartość dopasowania, 
-				// to zapamiętaj to przekształcenie geometryczne i dopasowanie
-				dMinSimValue = dCurrentSim;
-				iMinSimTransform = iTransform;
-				SAFEDELETE(pMinSimMatching);
-				pMinSimMatching = new SimilMatching(*m_pMatching);
-				assert(pMinSimMatching != NULL);
+			// now the original positions of model 0 are stored
+
+
+			int iTransform; // a counter of geometric transformations
+			for (iTransform = 0; iTransform < NO_OF_TRANSFORM; iTransform++)
+			{
+				// for each geometric transformation to be done
+				// entry conditions:
+				// - models (m_Mod) exist and are available
+				// - matching (m_pMatching) exists and is empty
+				// - all properties are created and available (m_aDegrees and m_aPositions)
+
+				// recompute geometric properties according to the transformation iTransform
+				// but only for model 0
+				for (iPart = 0; iPart < m_Mod[0]->getPartCount(); iPart++)
+				{
+					// for each iPart, a part of the model iMod
+					m_aPositions[0][iPart].x *= dMulX[iTransform];
+					m_aPositions[0][iPart].y *= dMulY[iTransform];
+					m_aPositions[0][iPart].z *= dMulZ[iTransform];
+				}
+				// now the positions are recomputed
+				ComputeMatching();
+
+				// teraz m_pMatching istnieje i jest pełne
+				assert(m_pMatching != NULL);
+				assert(m_pMatching->IsFull() == true);
+
+				// wykorzystaj to dopasowanie i geometrię do obliczenia tymczasowej wartości miary
+				int iTempRes = CountPartsDistance();
+				// załóż sukces
+				assert(iTempRes == 1);
+				double dCurrentSim = m_adFactors[0] * double(m_iDV) +
+					m_adFactors[1] * double(m_iDD) +
+					m_adFactors[2] * double(m_iDN) +
+					m_adFactors[3] * double(m_dDG);
+				// załóż poprawną wartość podobieństwa
+				assert(dCurrentSim >= 0.0);
+
+				// porównaj wartość obliczoną z dotychczasowym minimum
+				if (dCurrentSim < dMinSimValue)
+				{
+					// jeśli uzyskano mniejszą wartość dopasowania, 
+					// to zapamiętaj to przekształcenie geometryczne i dopasowanie
+					dMinSimValue = dCurrentSim;
+					iMinSimTransform = iTransform;
+					SAFEDELETE(pMinSimMatching);
+					pMinSimMatching = new SimilMatching(*m_pMatching);
+					assert(pMinSimMatching != NULL);
+				}
+
+				// teraz już można usunąć stare dopasowanie (dla potrzeb następnego przebiegu pętli)
+				m_pMatching->Empty();
+			} // for ( iTransform )
+
+			// po pętli przywróć najlepsze dopasowanie
+			delete m_pMatching;
+			m_pMatching = pMinSimMatching;
+
+			DB(printf("Matching has been chosen!\n");)
+				// print the name of the chosen transformation:
+				// printf("Chosen transformation: %s\n", szTransformNames[ iMinSimTransform ] );
+
+				// i przywróć jednocześnie pozycje Parts modelu 0 dla tego dopasowania
+				// - najpierw przywroc Parts pozycje orginalne, po SVD
+				for (iPart = 0; iPart < m_Mod[0]->getPartCount(); iPart++)
+				{
+					m_aPositions[0][iPart].x = StoredPositions[iPart].x;
+					m_aPositions[0][iPart].y = StoredPositions[iPart].y;
+					m_aPositions[0][iPart].z = StoredPositions[iPart].z;
+				}
+			// - usun teraz zapamietane pozycje Parts
+			delete[] StoredPositions;
+			// - a teraz oblicz na nowo wspolrzedne wlasciwego przeksztalcenia dla model 0
+			for (iTransform = 0; iTransform <= iMinSimTransform; iTransform++)
+			{
+				// for each transformation before and INCLUDING iMinTransform
+				// do the transformation (only model 0)
+				for (iPart = 0; iPart < m_Mod[0]->getPartCount(); iPart++)
+				{
+					m_aPositions[0][iPart].x *= dMulX[iTransform];
+					m_aPositions[0][iPart].y *= dMulY[iTransform];
+					m_aPositions[0][iPart].z *= dMulZ[iTransform];
+				}
 			}
-
-			// teraz już można usunąć stare dopasowanie (dla potrzeb następnego przebiegu pętli)
-			m_pMatching->Empty();
-		} // for ( iTransform )
-
-		// po pętli przywróć najlepsze dopasowanie
-		delete m_pMatching;
-		m_pMatching = pMinSimMatching;
-
-		DB(printf("Matching has been chosen!\n");)
-			// print the name of the chosen transformation:
-			// printf("Chosen transformation: %s\n", szTransformNames[ iMinSimTransform ] );
-
-			// i przywróć jednocześnie pozycje Parts modelu 0 dla tego dopasowania
-			// - najpierw przywroc Parts pozycje orginalne, po SVD
-			for (iPart = 0; iPart < m_Mod[0]->getPartCount(); iPart++)
-			{
-			m_aPositions[0][iPart].x = StoredPositions[iPart].x;
-			m_aPositions[0][iPart].y = StoredPositions[iPart].y;
-			m_aPositions[0][iPart].z = StoredPositions[iPart].z;
-			}
-		// - usun teraz zapamietane pozycje Parts
-		delete[] StoredPositions;
-		// - a teraz oblicz na nowo wspolrzedne wlasciwego przeksztalcenia dla model 0
-		for (iTransform = 0; iTransform <= iMinSimTransform; iTransform++)
-		{
-			// for each transformation before and INCLUDING iMinTransform
-			// do the transformation (only model 0)
-			for (iPart = 0; iPart < m_Mod[0]->getPartCount(); iPart++)
-			{
-				m_aPositions[0][iPart].x *= dMulX[iTransform];
-				m_aPositions[0][iPart].y *= dMulY[iTransform];
-				m_aPositions[0][iPart].z *= dMulZ[iTransform];
-			}
-		}
 
 		}
@@ -1855,5 +1818,4 @@
 
 	DB(_PrintPartsMatching();)
-
 
 		return 1;
@@ -1985,5 +1947,5 @@
 		int nSize = m_Mod[iMod]->getPartCount();
 
-		double *pDistances =  new double[nSize * nSize];
+		double *pDistances = new double[nSize * nSize];
 
 		for (int i = 0; i < nSize; i++)
@@ -1997,21 +1959,21 @@
 		Pt3D P1Pos, P2Pos; // positions of parts 
 		double dDistance; // the distance between Parts
-		
+
 		double *weights = new double[nSize];
 		for (int i = 0; i < nSize; i++)
 		{
-			if (wMDS==1)
+			if (wMDS == 1)
 				weights[i] = 0;
 			else
 				weights[i] = 1;
 		}
-				
-		if (wMDS==1)
+
+		if (wMDS == 1)
 			for (int i = 0; i < pModel->getJointCount(); i++)
 			{
 				weights[pModel->getJoint(i)->p1_refno]++;
-				weights[pModel->getJoint(i)->p2_refno]++; 
+				weights[pModel->getJoint(i)->p2_refno]++;
 			}
-		
+
 		for (iP1 = 0; iP1 < pModel->getPartCount(); iP1++)
 		{
@@ -2057,4 +2019,16 @@
 }
 
+/**	Evaluates distance between two given genotypes. The distance depends strongly
+	on weights set and the matching algorithm used.
+	@param G0 Pointer to the first of compared genotypes
+	@param G1 Pointer to the second of compared genotypes.
+	@return Distance between two genotypes.
+	@sa m_adFactors, matching_method
+	*/
+double ModelSimil::EvaluateDistance(const Geno *G0, const Geno *G1)
+{
+	return matching_method == 0 ? EvaluateDistanceHungarian(G0, G1) : EvaluateDistanceGreedy(G0, G1);
+}
+
 void ModelSimil::p_evaldistance(ExtValue *args, ExtValue *ret)
 {
@@ -2066,2 +2040,168 @@
 		ret->setDouble(EvaluateDistance(g1, g2));
 }
+
+void ModelSimil::FillPartsDistances(double*& dist, int bigger, int smaller, bool geo)
+{
+	for (int i = 0; i < bigger; i++)
+	{
+		for (int j = 0; j < bigger; j++)
+		{
+			// assign penalty for unassignment for vertex from bigger model
+			if (j >= smaller)
+			{
+				if (geo)
+					dist[i*bigger + j] += m_adFactors[3] * m_aPositions[1 - m_iSmaller][i].length();
+				else
+					dist[i*bigger + j] = m_adFactors[1] * m_aDegrees[1 - m_iSmaller][i][DEGREE] + m_adFactors[2] * m_aDegrees[1 - m_iSmaller][i][NEURONS];
+			}
+			// compute distance between parts
+			else
+			{
+				if (geo)
+					dist[i*bigger + j] += m_adFactors[3] * m_aPositions[1 - m_iSmaller][i].distanceTo(m_aPositions[m_iSmaller][j]);
+				else
+					dist[i*bigger + j] = m_adFactors[1] * abs(m_aDegrees[1 - m_iSmaller][i][DEGREE] - m_aDegrees[m_iSmaller][j][DEGREE])
+					+ m_adFactors[2] * abs(m_aDegrees[1 - m_iSmaller][i][NEURONS] - m_aDegrees[m_iSmaller][j][NEURONS]);
+			}
+
+		}
+	}
+}
+
+double ModelSimil::EvaluateDistanceHungarian(const Geno *G0, const Geno *G1)
+{
+	double dResult = 0;
+
+	m_Gen[0] = G0;
+	m_Gen[1] = G1;
+
+	// check whether pointers are not NULL
+	if (m_Gen[0] == NULL || m_Gen[1] == NULL)
+	{
+		DB(printf("ModelSimil::EvaluateDistanceHungarian - invalid genotypes pointers\n");)
+			return 0.0;
+	}
+	// create models of objects to compare
+	m_Mod[0] = new Model(*(m_Gen[0]));
+	m_Mod[1] = new Model(*(m_Gen[1]));
+
+	// validate models
+	if (m_Mod[0] == NULL || m_Mod[1] == NULL || !(m_Mod[0]->isValid()) || !(m_Mod[1]->isValid()))
+	{
+		DB(printf("ModelSimil::EvaluateDistanceHungarian - invalid models pointers\n");)
+			return 0.0;
+	}
+
+	//Get information about vertex degrees, neurons and 3D location
+	if (!CreatePartInfoTables())
+		return 0;
+	if (!CountPartDegrees())
+		return 0;
+	if (!GetPartPositions())
+		return 0;
+	if (!CountPartNeurons())
+		return 0;
+
+	m_iSmaller = m_Mod[0]->getPartCount() <= m_Mod[1]->getPartCount() ? 0 : 1;
+	int nSmaller = m_Mod[m_iSmaller]->getPartCount();
+	int nBigger = m_Mod[1 - m_iSmaller]->getPartCount();
+
+	double* partsDistances = new double[nBigger*nBigger]();
+	FillPartsDistances(partsDistances, nBigger, nSmaller, false);
+	int *assignment = new int[nBigger]();
+
+	HungarianAlgorithm hungarian;
+
+	if (m_adFactors[3] > 0)
+	{
+		if (!ComputePartsPositionsBySVD())
+		{
+			return 0;
+		}
+
+		// tutaj zacznij pętlę po przekształceniach  geometrycznych
+		const int NO_OF_TRANSFORM = 8; // liczba transformacji geometrycznych (na razie tylko ID i O_YZ)
+		// tablice transformacji współrzędnych; nie są to dokładnie tablice tranformacji, ale raczej tablice PRZEJŚĆ
+		// pomiędzy transformacjami; 
+		const int dMulX[NO_OF_TRANSFORM] = { 1, -1, -1, 1, -1, 1, -1, -1 };
+		const int dMulY[NO_OF_TRANSFORM] = { 1, 1, -1, -1, -1, -1, -1, 1 };
+		const int dMulZ[NO_OF_TRANSFORM] = { 1, 1, 1, -1, -1, -1, 1, 1 };
+
+		std::vector<int> minAssignment(nBigger);
+#ifdef max
+#undef max //this macro would conflict with line below
+#endif
+		double dMinSimValue = std::numeric_limits<double>::max(); // minimum value of similarity
+
+		int iTransform; // a counter of geometric transformations
+		for (iTransform = 0; iTransform < NO_OF_TRANSFORM; iTransform++)
+		{
+			// for each geometric transformation to be done
+			// entry conditions:
+			// - models (m_Mod) exist and are available
+			// - all properties are created and available (m_aDegrees and m_aPositions)
+			double* tmpPartsDistances = new double[nBigger*nBigger]();
+			std::copy(partsDistances, partsDistances + nBigger * nBigger, tmpPartsDistances);
+			// recompute geometric properties according to the transformation iTransform
+			// but only for model 0
+			for (int iPart = 0; iPart < m_Mod[m_iSmaller]->getPartCount(); iPart++)
+			{
+				// for each iPart, a part of the model iMod
+				m_aPositions[m_iSmaller][iPart].x *= dMulX[iTransform];
+				m_aPositions[m_iSmaller][iPart].y *= dMulY[iTransform];
+				m_aPositions[m_iSmaller][iPart].z *= dMulZ[iTransform];
+			}
+			// now the positions are recomputed
+
+			FillPartsDistances(tmpPartsDistances, nBigger, nSmaller, true);
+			std::fill_n(assignment, nBigger, 0);
+			double dCurrentSim = hungarian.Solve(tmpPartsDistances, assignment, nBigger, nBigger);
+
+			delete[] tmpPartsDistances;
+			// załóż poprawną wartość podobieństwa
+			assert(dCurrentSim >= 0.0);
+
+			// porównaj wartość obliczoną z dotychczasowym minimum
+			if (dCurrentSim < dMinSimValue)
+			{
+				dMinSimValue = dCurrentSim;
+				if (saveMatching == 1)
+				{
+					minAssignment.clear();
+					minAssignment.insert(minAssignment.begin(), assignment, assignment + nBigger);
+				}
+			}
+		}
+
+		dResult = dMinSimValue;
+		if (saveMatching == 1)
+			std::copy(minAssignment.begin(), minAssignment.end(), assignment);
+	}
+
+	else
+	{
+		dResult = hungarian.Solve(partsDistances, assignment, nBigger, nBigger);
+	}
+
+	//add difference in anywhere and onJoint neurons
+	dResult += m_adFactors[2] * (abs(m_aOnJoint[0][3] - m_aOnJoint[1][3]) + abs(m_aAnywhere[0][3] - m_aAnywhere[1][3]));
+	//add difference in part numbers
+	dResult += (nBigger - nSmaller) * m_adFactors[0];
+
+	// delete degree arrays created in CreatePartInfoTables 
+	SAFEDELETEARRAY(m_aDegrees[0]);
+	SAFEDELETEARRAY(m_aDegrees[1]);
+
+	// and position arrays
+	SAFEDELETEARRAY(m_aPositions[0]);
+	SAFEDELETEARRAY(m_aPositions[1]);
+
+	// delete created models
+	SAFEDELETE(m_Mod[0]);
+	SAFEDELETE(m_Mod[1]);
+
+	delete[] assignment;
+	delete[] partsDistances;
+
+	return dResult;
+}
