Index: cpp/frams/model/similarity/hungarian/hungarian.cpp
===================================================================
--- cpp/frams/model/similarity/hungarian/hungarian.cpp	(revision 869)
+++ cpp/frams/model/similarity/hungarian/hungarian.cpp	(revision 869)
@@ -0,0 +1,380 @@
+//Source: https://github.com/mcximing/hungarian-algorithm-cpp/blob/master/Hungarian.cpp
+//HungarianAlgorithm::Solve() was changed to take arrays instead of vectors as an input
+
+///////////////////////////////////////////////////////////////////////////////
+// Hungarian.cpp: Implementation file for Class HungarianAlgorithm.
+// 
+// This is a C++ wrapper with slight modification of a hungarian algorithm implementation by Markus Buehren.
+// The original implementation is a few mex-functions for use in MATLAB, found here:
+// http://www.mathworks.com/matlabcentral/fileexchange/6543-functions-for-the-rectangular-assignment-problem
+// 
+// Both this code and the orignal code are published under the BSD license.
+// by Cong Ma, 2016
+
+
+#include <stdlib.h>
+#include <cfloat> // for DBL_MAX
+#include <math.h>  // for fabs()
+#include "hungarian.h"
+#include <common/log.h>
+
+
+HungarianAlgorithm::HungarianAlgorithm() {}
+HungarianAlgorithm::~HungarianAlgorithm() {}
+
+
+//********************************************************//
+// A single function wrapper for solving assignment problem.
+//********************************************************//
+double HungarianAlgorithm::Solve(double *&distMatrixIn, int *&assignment, int nRows, int nCols)
+{
+	double cost = 0.0;
+
+	// call solving function
+	assignmentoptimal(assignment, &cost, distMatrixIn, nRows, nCols);
+
+	return cost;
+}
+
+
+//********************************************************//
+// Solve optimal solution for assignment problem using Munkres algorithm, also known as Hungarian Algorithm.
+//********************************************************//
+void HungarianAlgorithm::assignmentoptimal(int *assignment, double *cost, double *distMatrixIn, int nOfRows, int nOfColumns)
+{
+	double *distMatrix, *distMatrixTemp, *distMatrixEnd, *columnEnd, value, minValue;
+	bool *coveredColumns, *coveredRows, *starMatrix, *newStarMatrix, *primeMatrix;
+	int nOfElements, minDim, row, col;
+
+	/* initialization */
+	*cost = 0;
+	for (row = 0; row < nOfRows; row++)
+		assignment[row] = -1;
+
+	/* generate working copy of distance Matrix */
+	/* check if all matrix elements are positive */
+	nOfElements = nOfRows * nOfColumns;
+	distMatrix = (double *)malloc(nOfElements * sizeof(double));
+	distMatrixEnd = distMatrix + nOfElements;
+
+	for (row = 0; row < nOfElements; row++)
+	{
+		value = distMatrixIn[row];
+		if (value < 0)
+			logPrintf("HungarianAlgorithm", "assignmentoptimal", LOG_ERROR, "All matrix elements have to be non-negative, not %g", value);
+		distMatrix[row] = value;
+	}
+
+
+	/* memory allocation */
+	coveredColumns = (bool *)calloc(nOfColumns, sizeof(bool));
+	coveredRows = (bool *)calloc(nOfRows, sizeof(bool));
+	starMatrix = (bool *)calloc(nOfElements, sizeof(bool));
+	primeMatrix = (bool *)calloc(nOfElements, sizeof(bool));
+	newStarMatrix = (bool *)calloc(nOfElements, sizeof(bool)); /* used in step4 */
+
+	/* preliminary steps */
+	if (nOfRows <= nOfColumns)
+	{
+		minDim = nOfRows;
+
+		for (row = 0; row < nOfRows; row++)
+		{
+			/* find the smallest element in the row */
+			distMatrixTemp = distMatrix + row;
+			minValue = *distMatrixTemp;
+			distMatrixTemp += nOfRows;
+			while (distMatrixTemp < distMatrixEnd)
+			{
+				value = *distMatrixTemp;
+				if (value < minValue)
+					minValue = value;
+				distMatrixTemp += nOfRows;
+			}
+
+			/* subtract the smallest element from each element of the row */
+			distMatrixTemp = distMatrix + row;
+			while (distMatrixTemp < distMatrixEnd)
+			{
+				*distMatrixTemp -= minValue;
+				distMatrixTemp += nOfRows;
+			}
+		}
+
+		/* Steps 1 and 2a */
+		for (row = 0; row < nOfRows; row++)
+			for (col = 0; col < nOfColumns; col++)
+				if (fabs(distMatrix[row + nOfRows * col]) < DBL_EPSILON)
+					if (!coveredColumns[col])
+					{
+						starMatrix[row + nOfRows * col] = true;
+						coveredColumns[col] = true;
+						break;
+					}
+	}
+	else /* if(nOfRows > nOfColumns) */
+	{
+		minDim = nOfColumns;
+
+		for (col = 0; col < nOfColumns; col++)
+		{
+			/* find the smallest element in the column */
+			distMatrixTemp = distMatrix + nOfRows * col;
+			columnEnd = distMatrixTemp + nOfRows;
+
+			minValue = *distMatrixTemp++;
+			while (distMatrixTemp < columnEnd)
+			{
+				value = *distMatrixTemp++;
+				if (value < minValue)
+					minValue = value;
+			}
+
+			/* subtract the smallest element from each element of the column */
+			distMatrixTemp = distMatrix + nOfRows * col;
+			while (distMatrixTemp < columnEnd)
+				*distMatrixTemp++ -= minValue;
+		}
+
+		/* Steps 1 and 2a */
+		for (col = 0; col < nOfColumns; col++)
+			for (row = 0; row < nOfRows; row++)
+				if (fabs(distMatrix[row + nOfRows * col]) < DBL_EPSILON)
+					if (!coveredRows[row])
+					{
+						starMatrix[row + nOfRows * col] = true;
+						coveredColumns[col] = true;
+						coveredRows[row] = true;
+						break;
+					}
+		for (row = 0; row < nOfRows; row++)
+			coveredRows[row] = false;
+
+	}
+
+	/* move to step 2b */
+	step2b(assignment, distMatrix, starMatrix, newStarMatrix, primeMatrix, coveredColumns, coveredRows, nOfRows, nOfColumns, minDim);
+
+	/* compute cost and remove invalid assignments */
+	computeassignmentcost(assignment, cost, distMatrixIn, nOfRows);
+
+	/* free allocated memory */
+	free(distMatrix);
+	free(coveredColumns);
+	free(coveredRows);
+	free(starMatrix);
+	free(primeMatrix);
+	free(newStarMatrix);
+
+	return;
+}
+
+/********************************************************/
+void HungarianAlgorithm::buildassignmentvector(int *assignment, bool *starMatrix, int nOfRows, int nOfColumns)
+{
+	int row, col;
+
+	for (row = 0; row < nOfRows; row++)
+		for (col = 0; col < nOfColumns; col++)
+			if (starMatrix[row + nOfRows * col])
+			{
+#ifdef ONE_INDEXING
+				assignment[row] = col + 1; /* MATLAB-Indexing */
+#else
+				assignment[row] = col;
+#endif
+				break;
+			}
+}
+
+/********************************************************/
+void HungarianAlgorithm::computeassignmentcost(int *assignment, double *cost, double *distMatrix, int nOfRows)
+{
+	int row, col;
+
+	for (row = 0; row < nOfRows; row++)
+	{
+		col = assignment[row];
+		if (col >= 0)
+			*cost += distMatrix[row + nOfRows * col];
+	}
+}
+
+/********************************************************/
+void HungarianAlgorithm::step2a(int *assignment, double *distMatrix, bool *starMatrix, bool *newStarMatrix, bool *primeMatrix, bool *coveredColumns, bool *coveredRows, int nOfRows, int nOfColumns, int minDim)
+{
+	bool *starMatrixTemp, *columnEnd;
+	int col;
+
+	/* cover every column containing a starred zero */
+	for (col = 0; col < nOfColumns; col++)
+	{
+		starMatrixTemp = starMatrix + nOfRows * col;
+		columnEnd = starMatrixTemp + nOfRows;
+		while (starMatrixTemp < columnEnd) {
+			if (*starMatrixTemp++)
+			{
+				coveredColumns[col] = true;
+				break;
+			}
+		}
+	}
+
+	/* move to step 3 */
+	step2b(assignment, distMatrix, starMatrix, newStarMatrix, primeMatrix, coveredColumns, coveredRows, nOfRows, nOfColumns, minDim);
+}
+
+/********************************************************/
+void HungarianAlgorithm::step2b(int *assignment, double *distMatrix, bool *starMatrix, bool *newStarMatrix, bool *primeMatrix, bool *coveredColumns, bool *coveredRows, int nOfRows, int nOfColumns, int minDim)
+{
+	int col, nOfCoveredColumns;
+
+	/* count covered columns */
+	nOfCoveredColumns = 0;
+	for (col = 0; col < nOfColumns; col++)
+		if (coveredColumns[col])
+			nOfCoveredColumns++;
+
+	if (nOfCoveredColumns == minDim)
+	{
+		/* algorithm finished */
+		buildassignmentvector(assignment, starMatrix, nOfRows, nOfColumns);
+	}
+	else
+	{
+		/* move to step 3 */
+		step3(assignment, distMatrix, starMatrix, newStarMatrix, primeMatrix, coveredColumns, coveredRows, nOfRows, nOfColumns, minDim);
+	}
+
+}
+
+/********************************************************/
+void HungarianAlgorithm::step3(int *assignment, double *distMatrix, bool *starMatrix, bool *newStarMatrix, bool *primeMatrix, bool *coveredColumns, bool *coveredRows, int nOfRows, int nOfColumns, int minDim)
+{
+	bool zerosFound;
+	int row, col, starCol;
+
+	zerosFound = true;
+	while (zerosFound)
+	{
+		zerosFound = false;
+		for (col = 0; col < nOfColumns; col++)
+			if (!coveredColumns[col])
+				for (row = 0; row < nOfRows; row++)
+					if ((!coveredRows[row]) && (fabs(distMatrix[row + nOfRows * col]) < DBL_EPSILON))
+					{
+						/* prime zero */
+						primeMatrix[row + nOfRows * col] = true;
+
+						/* find starred zero in current row */
+						for (starCol = 0; starCol < nOfColumns; starCol++)
+							if (starMatrix[row + nOfRows * starCol])
+								break;
+
+						if (starCol == nOfColumns) /* no starred zero found */
+						{
+							/* move to step 4 */
+							step4(assignment, distMatrix, starMatrix, newStarMatrix, primeMatrix, coveredColumns, coveredRows, nOfRows, nOfColumns, minDim, row, col);
+							return;
+						}
+						else
+						{
+							coveredRows[row] = true;
+							coveredColumns[starCol] = false;
+							zerosFound = true;
+							break;
+						}
+					}
+	}
+
+	/* move to step 5 */
+	step5(assignment, distMatrix, starMatrix, newStarMatrix, primeMatrix, coveredColumns, coveredRows, nOfRows, nOfColumns, minDim);
+}
+
+/********************************************************/
+void HungarianAlgorithm::step4(int *assignment, double *distMatrix, bool *starMatrix, bool *newStarMatrix, bool *primeMatrix, bool *coveredColumns, bool *coveredRows, int nOfRows, int nOfColumns, int minDim, int row, int col)
+{
+	int n, starRow, starCol, primeRow, primeCol;
+	int nOfElements = nOfRows * nOfColumns;
+
+	/* generate temporary copy of starMatrix */
+	for (n = 0; n < nOfElements; n++)
+		newStarMatrix[n] = starMatrix[n];
+
+	/* star current zero */
+	newStarMatrix[row + nOfRows * col] = true;
+
+	/* find starred zero in current column */
+	starCol = col;
+	for (starRow = 0; starRow < nOfRows; starRow++)
+		if (starMatrix[starRow + nOfRows * starCol])
+			break;
+
+	while (starRow < nOfRows)
+	{
+		/* unstar the starred zero */
+		newStarMatrix[starRow + nOfRows * starCol] = false;
+
+		/* find primed zero in current row */
+		primeRow = starRow;
+		for (primeCol = 0; primeCol < nOfColumns; primeCol++)
+			if (primeMatrix[primeRow + nOfRows * primeCol])
+				break;
+
+		/* star the primed zero */
+		newStarMatrix[primeRow + nOfRows * primeCol] = true;
+
+		/* find starred zero in current column */
+		starCol = primeCol;
+		for (starRow = 0; starRow < nOfRows; starRow++)
+			if (starMatrix[starRow + nOfRows * starCol])
+				break;
+	}
+
+	/* use temporary copy as new starMatrix */
+	/* delete all primes, uncover all rows */
+	for (n = 0; n < nOfElements; n++)
+	{
+		primeMatrix[n] = false;
+		starMatrix[n] = newStarMatrix[n];
+	}
+	for (n = 0; n < nOfRows; n++)
+		coveredRows[n] = false;
+
+	/* move to step 2a */
+	step2a(assignment, distMatrix, starMatrix, newStarMatrix, primeMatrix, coveredColumns, coveredRows, nOfRows, nOfColumns, minDim);
+}
+
+/********************************************************/
+void HungarianAlgorithm::step5(int *assignment, double *distMatrix, bool *starMatrix, bool *newStarMatrix, bool *primeMatrix, bool *coveredColumns, bool *coveredRows, int nOfRows, int nOfColumns, int minDim)
+{
+	double h, value;
+	int row, col;
+
+	/* find smallest uncovered element h */
+	h = DBL_MAX;
+	for (row = 0; row < nOfRows; row++)
+		if (!coveredRows[row])
+			for (col = 0; col < nOfColumns; col++)
+				if (!coveredColumns[col])
+				{
+					value = distMatrix[row + nOfRows * col];
+					if (value < h)
+						h = value;
+				}
+
+	/* add h to each covered row */
+	for (row = 0; row < nOfRows; row++)
+		if (coveredRows[row])
+			for (col = 0; col < nOfColumns; col++)
+				distMatrix[row + nOfRows * col] += h;
+
+	/* subtract h from each uncovered column */
+	for (col = 0; col < nOfColumns; col++)
+		if (!coveredColumns[col])
+			for (row = 0; row < nOfRows; row++)
+				distMatrix[row + nOfRows * col] -= h;
+
+	/* move to step 3 */
+	step3(assignment, distMatrix, starMatrix, newStarMatrix, primeMatrix, coveredColumns, coveredRows, nOfRows, nOfColumns, minDim);
+}
Index: cpp/frams/model/similarity/hungarian/hungarian.h
===================================================================
--- cpp/frams/model/similarity/hungarian/hungarian.h	(revision 869)
+++ cpp/frams/model/similarity/hungarian/hungarian.h	(revision 869)
@@ -0,0 +1,35 @@
+//Source: https://github.com/mcximing/hungarian-algorithm-cpp/blob/master/Hungarian.h
+///////////////////////////////////////////////////////////////////////////////
+// Hungarian.h: Header file for Class HungarianAlgorithm.
+// 
+// This is a C++ wrapper with slight modification of a hungarian algorithm implementation by Markus Buehren.
+// The original implementation is a few mex-functions for use in MATLAB, found here:
+// http://www.mathworks.com/matlabcentral/fileexchange/6543-functions-for-the-rectangular-assignment-problem
+// 
+// Both this code and the orignal code are published under the BSD license.
+// by Cong Ma, 2016
+
+
+#ifndef HUNGARIAN_H
+#define HUNGARIAN_H
+
+class HungarianAlgorithm
+{
+public:
+	HungarianAlgorithm();
+	~HungarianAlgorithm();
+	double Solve(double *&distMatrixIn, int *&assignment, int nRows, int nCols);
+
+private:
+	void assignmentoptimal(int *assignment, double *cost, double *distMatrix, int nOfRows, int nOfColumns);
+	void buildassignmentvector(int *assignment, bool *starMatrix, int nOfRows, int nOfColumns);
+	void computeassignmentcost(int *assignment, double *cost, double *distMatrix, int nOfRows);
+	void step2a(int *assignment, double *distMatrix, bool *starMatrix, bool *newStarMatrix, bool *primeMatrix, bool *coveredColumns, bool *coveredRows, int nOfRows, int nOfColumns, int minDim);
+	void step2b(int *assignment, double *distMatrix, bool *starMatrix, bool *newStarMatrix, bool *primeMatrix, bool *coveredColumns, bool *coveredRows, int nOfRows, int nOfColumns, int minDim);
+	void step3(int *assignment, double *distMatrix, bool *starMatrix, bool *newStarMatrix, bool *primeMatrix, bool *coveredColumns, bool *coveredRows, int nOfRows, int nOfColumns, int minDim);
+	void step4(int *assignment, double *distMatrix, bool *starMatrix, bool *newStarMatrix, bool *primeMatrix, bool *coveredColumns, bool *coveredRows, int nOfRows, int nOfColumns, int minDim, int row, int col);
+	void step5(int *assignment, double *distMatrix, bool *starMatrix, bool *newStarMatrix, bool *primeMatrix, bool *coveredColumns, bool *coveredRows, int nOfRows, int nOfColumns, int minDim);
+};
+
+
+#endif
Index: cpp/frams/model/similarity/simil_match.cpp
===================================================================
--- cpp/frams/model/similarity/simil_match.cpp	(revision 818)
+++ cpp/frams/model/similarity/simil_match.cpp	(revision 869)
@@ -1,4 +1,4 @@
 // This file is a part of Framsticks SDK.  http://www.framsticks.com/
-// Copyright (C) 1999-2015  Maciej Komosinski and Szymon Ulatowski.
+// Copyright (C) 1999-2019  Maciej Komosinski and Szymon Ulatowski.
 // See LICENSE.txt for details.
 
@@ -8,44 +8,44 @@
 
 /** Creates an empty matching for two objects of specified size.
-        @param Obj0Size Size of the first object. Must be positive.
-        @param Obj1Size Size of the second object. Must be positive.
+		@param Obj0Size Size of the first object. Must be positive.
+		@param Obj1Size Size of the second object. Must be positive.
  */
 SimilMatching::SimilMatching(int Obj0Size, int Obj1Size)
 {
-    // assure that sizes of objects are positive
-    assert(Obj0Size > 0);
-    assert(Obj1Size > 0);
-
-    // create necessary vectors
-    m_apvMatched[0] = new std::vector<int>(Obj0Size);
-    m_apvMatched[1] = new std::vector<int>(Obj1Size);
-
-    // assure that vectors are created
-    assert(m_apvMatched[0] != NULL);
-    assert(m_apvMatched[1] != NULL);
-
-    // fill vectors with "unmatched" indicator
+	// assure that sizes of objects are positive
+	assert(Obj0Size > 0);
+	assert(Obj1Size > 0);
+
+	// create necessary vectors
+	m_apvMatched[0] = new std::vector<int>(Obj0Size);
+	m_apvMatched[1] = new std::vector<int>(Obj1Size);
+
+	// assure that vectors are created
+	assert(m_apvMatched[0] != NULL);
+	assert(m_apvMatched[1] != NULL);
+
+	// fill vectors with "unmatched" indicator
 	for (unsigned int i = 0; i < m_apvMatched[0]->size(); i++)
-    {
-        m_apvMatched[0]->operator[](i) = -1;
-    }
+	{
+		m_apvMatched[0]->operator[](i) = -1;
+	}
 	for (unsigned int i = 0; i < m_apvMatched[1]->size(); i++)
-    {
-        m_apvMatched[1]->operator[](i) = -1;
-    }
+	{
+		m_apvMatched[1]->operator[](i) = -1;
+	}
 }
 
 /** A copying constructor.
-        @param Source The object to be copied.
+		@param Source The object to be copied.
  */
 SimilMatching::SimilMatching(const SimilMatching &Source)
 {
-    // copy the vectors of the actual matching
-    m_apvMatched[ 0 ] = new std::vector<int>(* (Source.m_apvMatched[ 0 ]));
-    m_apvMatched[ 1 ] = new std::vector<int>(* (Source.m_apvMatched[ 1 ]));
-
-    // assure that vectors are created
-    assert(m_apvMatched[0] != NULL);
-    assert(m_apvMatched[1] != NULL);
+	// copy the vectors of the actual matching
+	m_apvMatched[0] = new std::vector<int>(*(Source.m_apvMatched[0]));
+	m_apvMatched[1] = new std::vector<int>(*(Source.m_apvMatched[1]));
+
+	// assure that vectors are created
+	assert(m_apvMatched[0] != NULL);
+	assert(m_apvMatched[1] != NULL);
 }
 
@@ -54,147 +54,147 @@
 SimilMatching::~SimilMatching()
 {
-    // delete vectors of matching
-    delete m_apvMatched[0];
-    delete m_apvMatched[1];
+	// delete vectors of matching
+	delete m_apvMatched[0];
+	delete m_apvMatched[1];
 }
 
 /** Gets size of the specified object.
-        @param Index of an object (must be 0 or 1).
-        @return Size of the object (in elements).
+		@param Index of an object (must be 0 or 1).
+		@return Size of the object (in elements).
  */
 int SimilMatching::GetObjectSize(int Obj)
 {
-    // check parameter
-    assert((Obj == 0) || (Obj == 1));
-
-    // return the result
-    return m_apvMatched[ Obj ]->size();
+	// check parameter
+	assert((Obj == 0) || (Obj == 1));
+
+	// return the result
+	return m_apvMatched[Obj]->size();
 }
 
 /** Matches elements given by indices in the given objects.
-        @param Obj0 Index of the first object. Must be 0 or 1.
-        @param index0 Index of element in the first object. Must be a valid index 
-        ( >= 0 and < size of the object).
-        @param Obj1 Index of the second object. Must be 0 or 1 and different from Obj0.
-        @param Index1 index of element in the second object. Must be a valid index 
-        ( >= 0 and < size of the object).
+		@param Obj0 Index of the first object. Must be 0 or 1.
+		@param index0 Index of element in the first object. Must be a valid index
+		( >= 0 and < size of the object).
+		@param Obj1 Index of the second object. Must be 0 or 1 and different from Obj0.
+		@param Index1 index of element in the second object. Must be a valid index
+		( >= 0 and < size of the object).
 
  */
 void SimilMatching::Match(int Obj0, int Index0, int Obj1, int Index1)
 {
-    // check parameters of object 0
-    assert((Obj0 == 0) || (Obj0 == 1));
-    assert((Index0 >= 0) && (Index0 < (int) m_apvMatched[Obj0]->size()));
-    // check parameters of object 1
-    assert(((Obj1 == 0) || (Obj1 == 1)) && (Obj0 != Obj1));
-    assert((Index1 >= 0) && (Index1 < (int) m_apvMatched[Obj1]->size()));
-
-    // match given elements
-    // matching_Obj0(Index0) = Index1
-    m_apvMatched[ Obj0 ]->operator[](Index0) = Index1;
-    // matching_Obj1(Index1) = Index0
-    m_apvMatched[ Obj1 ]->operator[](Index1) = Index0;
+	// check parameters of object 0
+	assert((Obj0 == 0) || (Obj0 == 1));
+	assert((Index0 >= 0) && (Index0 < (int)m_apvMatched[Obj0]->size()));
+	// check parameters of object 1
+	assert(((Obj1 == 0) || (Obj1 == 1)) && (Obj0 != Obj1));
+	assert((Index1 >= 0) && (Index1 < (int)m_apvMatched[Obj1]->size()));
+
+	// match given elements
+	// matching_Obj0(Index0) = Index1
+	m_apvMatched[Obj0]->operator[](Index0) = Index1;
+	// matching_Obj1(Index1) = Index0
+	m_apvMatched[Obj1]->operator[](Index1) = Index0;
 }
 
 /** Checks if the given element in the given object is already matched.
-        @param Obj Index of an object (must be 0 or 1).
-        @param Index Index of an element in the given object. Must be a valid index 
-        ( >=0 and < size of the object).
-        @return true if the given element is matched, false otherwise.
+		@param Obj Index of an object (must be 0 or 1).
+		@param Index Index of an element in the given object. Must be a valid index
+		( >=0 and < size of the object).
+		@return true if the given element is matched, false otherwise.
  */
 bool SimilMatching::IsMatched(int Obj, int Index)
 {
-    // check parameters
-    assert((Obj == 0) || (Obj == 1));
-    assert((Index >= 0) && (Index < (int) m_apvMatched[ Obj ]->size()));
-
-    // check if the element is matched
-    if (m_apvMatched[ Obj ]->operator[](Index) >= 0)
-    {
-        return true;
-    }
-    else
-    {
-        return false;
-    }
+	// check parameters
+	assert((Obj == 0) || (Obj == 1));
+	assert((Index >= 0) && (Index < (int)m_apvMatched[Obj]->size()));
+
+	// check if the element is matched
+	if (m_apvMatched[Obj]->operator[](Index) >= 0)
+	{
+		return true;
+	}
+	else
+	{
+		return false;
+	}
 }
 
 /** Gets index of the element thet is matched in the other object withe the element given
-        by parameters.
-        @param Obj Index of an object (must be 0 or 1).
-        @param Index Index of checked element in the given object.
-        @return Index of the element (in the other organism) that is matched with the given
-        element. WARNING! If the given element is not matched, the result may be smaller than 0 
-        (check IsMatched() before using GetMatchedIndex()).
+		by parameters.
+		@param Obj Index of an object (must be 0 or 1).
+		@param Index Index of checked element in the given object.
+		@return Index of the element (in the other organism) that is matched with the given
+		element. WARNING! If the given element is not matched, the result may be smaller than 0
+		(check IsMatched() before using GetMatchedIndex()).
  */
 int SimilMatching::GetMatchedIndex(int Obj, int Index)
 {
-    // check parameters
-    assert((Obj == 0) || (Obj == 1));
-    assert((Index >= 0) && (Index < (int) m_apvMatched[ Obj ]->size()));
-
-    // return the index of the matched element
-    return m_apvMatched[ Obj ]->operator[](Index);
-}
-
-/** Checks if the matching is already full, i.e. if the smaller object already has all its 
-        elements matched. 
-        @return true if matching is full, false otherwise.
+	// check parameters
+	assert((Obj == 0) || (Obj == 1));
+	assert((Index >= 0) && (Index < (int)m_apvMatched[Obj]->size()));
+
+	// return the index of the matched element
+	return m_apvMatched[Obj]->operator[](Index);
+}
+
+/** Checks if the matching is already full, i.e. if the smaller object already has all its
+		elements matched.
+		@return true if matching is full, false otherwise.
  */
 bool SimilMatching::IsFull()
 {
-    // assume that the matching is full
-    bool bResult = true;
-    // index of the smallest object
-    int nObj;
-
-    // find the smallest object (its index)
-    if (m_apvMatched[ 0 ]->size() < m_apvMatched[ 1 ]->size())
-    {
-        nObj = 0;
-    }
-    else
-    {
-        nObj = 1;
-    }
-
-    // check if all elements of the smallest object are matched
+	// assume that the matching is full
+	bool bResult = true;
+	// index of the smallest object
+	int nObj;
+
+	// find the smallest object (its index)
+	if (m_apvMatched[0]->size() < m_apvMatched[1]->size())
+	{
+		nObj = 0;
+	}
+	else
+	{
+		nObj = 1;
+	}
+
+	// check if all elements of the smallest object are matched
 	for (unsigned int nElem = 0; nElem < m_apvMatched[nObj]->size(); nElem++)
-    {
-        if (m_apvMatched[ nObj ]->operator[](nElem) < 0)
-        {
-            // if any element is not matched, the result is false
-            bResult = false;
-            break;
-        }
-    }
-
-    // return the result
-    return bResult;
+	{
+		if (m_apvMatched[nObj]->operator[](nElem) < 0)
+		{
+			// if any element is not matched, the result is false
+			bResult = false;
+			break;
+		}
+	}
+
+	// return the result
+	return bResult;
 }
 
 /** Checks if the matching is empty (i.e. none of elements is matched).
-    @return true if matching is empty, otherwise - false.
+	@return true if matching is empty, otherwise - false.
  */
 bool SimilMatching::IsEmpty()
 {
-    // result - assume that matching is empty
-    bool bResult = true;
-
-    // matching is empty if either of objects has only unmatched elements
-    // so it may be first object 
-    int nObj = 0;
+	// result - assume that matching is empty
+	bool bResult = true;
+
+	// matching is empty if either of objects has only unmatched elements
+	// so it may be first object 
+	int nObj = 0;
 	for (unsigned int nElem = 0; nElem < m_apvMatched[nObj]->size(); nElem++)
-    {
-        if (m_apvMatched[ nObj ]->operator[](nElem) >= 0)
-        {
-            // if any element of the object is matched (unmatched objects have (-1))
-            bResult = false;
-            break;
-        }
-    }
-
-    // return the result from the loop
-    return bResult;
+	{
+		if (m_apvMatched[nObj]->operator[](nElem) >= 0)
+		{
+			// if any element of the object is matched (unmatched objects have (-1))
+			bResult = false;
+			break;
+		}
+	}
+
+	// return the result from the loop
+	return bResult;
 }
 
@@ -203,17 +203,17 @@
 void SimilMatching::Empty()
 {
-    for (int iObj = 0; iObj < 2; iObj++) // a counter of objects
-    {
-        // for each object in the matching
+	for (int iObj = 0; iObj < 2; iObj++) // a counter of objects
+	{
+		// for each object in the matching
 		for (unsigned int iElem = 0; iElem < m_apvMatched[iObj]->size(); iElem++) // a counter of objects' elements
-        {
-            // for each element iElem for the object iObj
-            // set it as unmatched (marker: -1)
-            m_apvMatched[ iObj ]->operator[](iElem) = -1;
-        }
-    }
-
-    // the exit condition
-    assert(IsEmpty() == true);
+		{
+			// for each element iElem for the object iObj
+			// set it as unmatched (marker: -1)
+			m_apvMatched[iObj]->operator[](iElem) = -1;
+		}
+	}
+
+	// the exit condition
+	assert(IsEmpty() == true);
 }
 
@@ -222,43 +222,43 @@
 void SimilMatching::PrintMatching()
 {
-    int nBigger;
-
-    // check which object is bigger
-    if (m_apvMatched[ 0 ]->size() >= m_apvMatched[ 1 ]->size())
-    {
-        nBigger = 0;
-    }
-    else
-    {
-        nBigger = 1;
-    }
-
-    // print first line - indices of objects
-    printf("[ ");
+	int nBigger;
+
+	// check which object is bigger
+	if (m_apvMatched[0]->size() >= m_apvMatched[1]->size())
+	{
+		nBigger = 0;
+	}
+	else
+	{
+		nBigger = 1;
+	}
+
+	// print first line - indices of objects
+	printf("[ ");
 	for (unsigned int i = 0; i < m_apvMatched[nBigger]->size(); i++)
-    {
-        printf("%2d ", i);
-    }
-    printf("]\n");
-
-    // print second line and third - indices of elements matched with elements of the objects
-    for (int nObj = 0; nObj < 2; nObj++)
-    {
-        // for both objects - print out lines of matched elements
-        printf("[ ");
+	{
+		printf("%2d ", i);
+	}
+	printf("]\n");
+
+	// print second line and third - indices of elements matched with elements of the objects
+	for (int nObj = 0; nObj < 2; nObj++)
+	{
+		// for both objects - print out lines of matched elements
+		printf("[ ");
 		for (unsigned int i = 0; i < m_apvMatched[nObj]->size(); i++)
-        {
-            if (IsMatched(nObj, i))
-            {
-                // if the element is matched - print the index
-                printf("%2d ", GetMatchedIndex(nObj, i));
-            }
-            else
-            {
-                // if the element is not matched - print "X"
-                printf(" X ");
-            }
-        }
-        printf("]\n");
-    }
-}
+		{
+			if (IsMatched(nObj, i))
+			{
+				// if the element is matched - print the index
+				printf("%2d ", GetMatchedIndex(nObj, i));
+			}
+			else
+			{
+				// if the element is not matched - print "X"
+				printf(" X ");
+			}
+		}
+		printf("]\n");
+	}
+}
Index: cpp/frams/model/similarity/simil_match.h
===================================================================
--- cpp/frams/model/similarity/simil_match.h	(revision 818)
+++ cpp/frams/model/similarity/simil_match.h	(revision 869)
@@ -1,4 +1,4 @@
 // This file is a part of Framsticks SDK.  http://www.framsticks.com/
-// Copyright (C) 1999-2015  Maciej Komosinski and Szymon Ulatowski.
+// Copyright (C) 1999-2019  Maciej Komosinski and Szymon Ulatowski.
 // See LICENSE.txt for details.
 
@@ -8,30 +8,30 @@
 #include <vector>
 
-/** Class describes a mutually single-valued function between two sets of elements 
-        (these sets are called objects). These sets may have different sizes, so this function
-        is mutually single-valued only for some subset of the bigger set.
-        This class is used while building a matching function between Parts of two Framsticks' 
-        organisms (similarity measure computation).
+/** This class describes a mutually single-valued function between two sets of elements
+	(these sets are called objects). These sets may have different sizes, so this function
+	is mutually single-valued only for some subset of the bigger set.
+	This class is used while building a matching function between Parts of two Framsticks'
+	organisms (similarity measure computation).
  */
 class SimilMatching
 {
 protected:
-    /** Array of pointers to two vectors. Each one stores indices of matched elements of
-            the other object. Value <0 means that an element is not matched yet. Sizes of these 
-            vectors are sizes of objects themselves.
-     */
-    std::vector<int> *m_apvMatched[2];
+	/** Array of pointers to two vectors. Each one stores indices of matched elements of
+		the other object. Value <0 means that an element is not matched yet. Sizes of these
+		vectors are sizes of objects themselves.
+	 */
+	std::vector<int> *m_apvMatched[2];
 public:
-    SimilMatching(int Obj0Size, int Obj1Size);
-    SimilMatching(const SimilMatching &Source);
-    ~SimilMatching();
-    int GetObjectSize(int Obj);
-    void Match(int Obj0, int Index0, int Obj1, int Index1);
-    bool IsMatched(int Obj, int Index);
-    int GetMatchedIndex(int Obj, int index);
-    bool IsFull();
-    bool IsEmpty();
-    void Empty();
-    void PrintMatching();
+	SimilMatching(int Obj0Size, int Obj1Size);
+	SimilMatching(const SimilMatching &Source);
+	~SimilMatching();
+	int GetObjectSize(int Obj);
+	void Match(int Obj0, int Index0, int Obj1, int Index1);
+	bool IsMatched(int Obj, int Index);
+	int GetMatchedIndex(int Obj, int index);
+	bool IsFull();
+	bool IsEmpty();
+	void Empty();
+	void PrintMatching();
 };
 
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;
+}
Index: cpp/frams/model/similarity/simil_model.h
===================================================================
--- cpp/frams/model/similarity/simil_model.h	(revision 818)
+++ cpp/frams/model/similarity/simil_model.h	(revision 869)
@@ -1,4 +1,4 @@
 // This file is a part of Framsticks SDK.  http://www.framsticks.com/
-// Copyright (C) 1999-2016  Maciej Komosinski and Szymon Ulatowski.
+// Copyright (C) 1999-2019  Maciej Komosinski and Szymon Ulatowski.
 // See LICENSE.txt for details.
 
@@ -9,5 +9,4 @@
 #include "frams/genetics/geno.h"
 #include "frams/model/model.h"
-#include "frams/util/3d.h"
 #include "simil_match.h"
 
@@ -27,5 +26,5 @@
  * Marek Kubiak (concept, implementation)
  * Maciej Komosinski (concept, Framsticks interface)
- * Agnieszka Mensfelt (refactoring)
+ * Agnieszka Mensfelt (refactoring, improvements)
  */
 class ModelSimil
@@ -34,7 +33,10 @@
 	ModelSimil();
 	virtual ~ModelSimil();
-	double EvaluateDistance(const Geno *G0, const Geno *G1);
+	double EvaluateDistance(const Geno *G0, const Geno *G1); //chooses greedy or hungarian
+	double EvaluateDistanceGreedy(const Geno *G0, const Geno *G1);
+	double EvaluateDistanceHungarian(const Geno *G0, const Geno *G1);
 
 	static int CompareDegrees(const void *pElem1, const void *pElem2);
+	static int CompareFuzzyDegrees(const void *pElem1, const void *pElem2);
 	static int CompareConnsNo(const void *pElem1, const void *pElem2);
 	static int GetNOFactors();
@@ -48,8 +50,9 @@
 	int MatchPartsGeometry();
 	void ComputeMatching();
+	void FillPartsDistances(double *&dist, int bigger, int smaller, bool geo);
 	void _PrintPartsMatching();
 	void SaveIntermediateFiles();
 
-	static int CheckPartsIdentity(Part *P1, Part *P2);
+	static int CheckPartsIdentity(Part *P1, Part *P2); //TODO exists?
 	int SortPartInfoTables();
 	int CountPartNeurons();
@@ -58,9 +61,8 @@
 	int CountPartDegrees();
 
-	int SortPartInfoFuzzy();
 	void CountFuzzyNeighb();
 	void SortFuzzyNeighb();
 	void GetNeighbIndexes(int mod, int partInd, std::vector<int> &indexes);
-	void CheckFuzzyIdentity();
+	void FuzzyOrder();
 
 	int CreatePartInfoTables();
@@ -68,4 +70,5 @@
 	void _PrintArray(int *array, int base, int size);
 	void _PrintNeighbourhood(int i);
+	void _PrintFuzzyNeighbourhood(int i);
 	void _PrintArrayDouble(double *array, int base, int size);
 	int CountPartsDistance();
@@ -73,4 +76,8 @@
 
 public:
+	/// Currently selected matching algorithm. Allowed values: 0 (more exact, slower), 1 (more greedy, faster). Details in https://doi.org/10.1007/978-3-030-16692-2_8
+	/// @sa EvaluateDistance
+	int matching_method;
+
 	/// Table of weights for weighted distance function.
 	/// Weights are for factors in the following order:
@@ -88,7 +95,10 @@
 	int fuzzyDepth;
 	int isFuzzy;
-        
+
 	//For wMDS = 1 weighted MDS with vertex degrees as weights is used for the alignment.
 	int wMDS;
+
+	//For saveMatching = 1 the best matching found will be saved.
+	int saveMatching;
 
 	/// Interface to local parameters
@@ -181,4 +191,3 @@
 };
 
-
 #endif
