Context Navigation

← Previous Changeset
Next Changeset →

Changeset 869

Timestamp:

05/02/19 23:50:27 (6 years ago)

Author:

Maciej Komosinski

Message:

Added another, improved way of calculating dissimilarity of two creatures/models. Details and comparisons in https://doi.org/10.1007/978-3-030-16692-2_8

Location:

cpp/frams

Files:

: 3 added
: 5 edited

Makefile-SDK-files (modified) (2 diffs)
model/similarity/hungarian (added)
model/similarity/hungarian/hungarian.cpp (added)
model/similarity/hungarian/hungarian.h (added)
model/similarity/simil_match.cpp (modified) (5 diffs)
model/similarity/simil_match.h (modified) (2 diffs)
model/similarity/simil_model.cpp (modified) (27 diffs)
model/similarity/simil_model.h (modified) (10 diffs)

Legend:

: Unmodified
: Added
: Removed

cpp/frams/Makefile-SDK-files

-                      r864
+                      r869
 # ALL_DIRS is later expanded by the shell, no spaces/newlines allowed, or it breaks
 ALL_DIRS={common,PrintFloat,frams,frams/canvas,frams/config,common/loggers,frams/genetics,frams/genetics/f0,frams/genetics/f1,frams/genetics/f2,frams/genetics/f3,frams/genetics/f4,frams/genetics/f5,frams/genetics/f6,frams/genetics/f7,frams/genetics/f8,frams/genetics/f9,frams/genetics/fn,frams/genetics/fF,frams/genetics/fT,frams/genetics/fB,frams/genetics/fH,frams/genetics/fL,frams/model,frams/neuro,frams/neuro/impl,frams/param,frams/test,frams/util,frams/vm/classes,common/virtfile,frams/_demos,frams/model/geometry,frams/_demos/geometry,frams/model/similarity,frams/model/similarity/SVD}
+ALL_DIRS={common,PrintFloat,frams,frams/canvas,frams/config,common/loggers,frams/genetics,frams/genetics/f0,frams/genetics/f1,frams/genetics/f2,frams/genetics/f3,frams/genetics/f4,frams/genetics/f5,frams/genetics/f6,frams/genetics/f7,frams/genetics/f8,frams/genetics/f9,frams/genetics/fn,frams/genetics/fF,frams/genetics/fT,frams/genetics/fB,frams/genetics/fH,frams/genetics/fL,frams/model,frams/neuro,frams/neuro/impl,frams/param,frams/test,frams/util,frams/vm/classes,common/virtfile,frams/_demos,frams/model/geometry,frams/_demos/geometry,frams/model/similarity,frams/model/similarity/hungarian,frams/model/similarity/SVD}
 GENMANF4=frams/genetics/f4/f4_oper.o
 …
 NEURO_OBJS=frams/neuro/neuroimpl.o frams/neuro/neurofactory.o frams/neuro/impl/neuroimpl-simple.o frams/neuro/impl/neuroimpl-channels.o frams/neuro/impl/neuroimpl-fuzzy.o frams/neuro/impl/neuroimpl-fuzzy-f0.o
 SIMILARITY_OBJS=frams/model/similarity/SVD/lapack.o frams/model/similarity/SVD/matrix_tools.o frams/model/similarity/simil_match.o frams/model/similarity/simil_model.o
+SIMILARITY_OBJS=frams/model/similarity/hungarian/hungarian.o frams/model/similarity/SVD/lapack.o frams/model/similarity/SVD/matrix_tools.o frams/model/similarity/simil_match.o frams/model/similarity/simil_model.o
 NN_LAYOUT_OBJS=frams/canvas/nn_layout_model.o frams/canvas/nn_simple_layout.o frams/canvas/nn_smart_layout.o

cpp/frams/model/similarity/simil_match.cpp

-                      r361
+                      r869
 // 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.
 …
 /** 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);
+}
 …
 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;
+}
 …
 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);
+}
 …
 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");
+        }
+}

cpp/frams/model/similarity/simil_match.h

-                      r357
+                      r869
 // 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.
 …
 #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();
 };

cpp/frams/model/similarity/simil_model.cpp

-                      r818
+                      r869
 // 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"
 …
 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]), "", },
 …
         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;
 …
         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
 …
         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;
+        }
 …
         if ((this->*pfMatchingFunction)() == 0)
+        {
                 DB(printf("ModelSimil::EvaluateDistance - MatchParts() error\n");)
+                DB(printf("ModelSimil::EvaluateDistanceGreedy - MatchParts() error\n");)
                         return 0.0;
+        }
 …
                 if (CountPartsDistance() == 0)
+                {
                 DB(printf("ModelSimil::EvaluateDistance - CountPartDistance() error\n");)
                         return 0.0;
+                        DB(printf("ModelSimil::EvaluateDistanceGreedy - CountPartDistance() error\n");)
+                                return 0.0;
+                }
 …
                 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
 …
                 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
 …
                                 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
 …
                                                 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
 …
+}
+/** 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).
 …
                                 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);
+                                }
+                                }
 …
+}
+//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;
+                                }
 …
         SortFuzzyNeighb();
+        FuzzyOrder();
+}
 …
         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
 …
                                         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
 …
+}
-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.
 …
                         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);)
+                                                }
+                                        }
+                                }
+                        }
+                        }
 …
                                         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],
+- 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(
+- 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(
+- 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;
 …
+                                                        }
+                                                }
+                                                // 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],
+- 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],
+- 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(
+- 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(
+- 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],
+- 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.
 …
                 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];
+                        }
+                }
+                }
 …
         DB(_PrintPartsMatching();)
                 return 1;
 …
                 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++)
 …
                 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++)
+                {
 …
+}
+/**     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)
+{
 …
                 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;
+}

cpp/frams/model/similarity/simil_model.h

-                      r817
+                      r869
 // 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.
 …
 #include "frams/genetics/geno.h"
 #include "frams/model/model.h"
-#include "frams/util/3d.h"
 #include "simil_match.h"
 …
  * Marek Kubiak (concept, implementation)
  * Maciej Komosinski (concept, Framsticks interface)
  * Agnieszka Mensfelt (refactoring)
+ * Agnieszka Mensfelt (refactoring, improvements)
  */
 class ModelSimil
 …
         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();
 …
         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();
 …
         int CountPartDegrees();
-        int SortPartInfoFuzzy();
         void CountFuzzyNeighb();
         void SortFuzzyNeighb();
         void GetNeighbIndexes(int mod, int partInd, std::vector<int> &indexes);
         void CheckFuzzyIdentity();
+        void FuzzyOrder();
         int CreatePartInfoTables();
 …
         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();
 …
 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:
 …
         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
 …
 };
 #endif

Note: See TracChangeset for help on using the changeset viewer.