Changeset 869


Ignore:
Timestamp:
05/02/19 23:50:27 (3 weeks 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

Legend:

Unmodified
Added
Removed
  • cpp/frams/Makefile-SDK-files

    r864 r869  
    22
    33# ALL_DIRS is later expanded by the shell, no spaces/newlines allowed, or it breaks
    4 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}
     4ALL_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}
    55
    66GENMANF4=frams/genetics/f4/f4_oper.o
     
    3333NEURO_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
    3434
    35 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
     35SIMILARITY_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
    3636
    3737NN_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  
    11// This file is a part of Framsticks SDK.  http://www.framsticks.com/
    2 // Copyright (C) 1999-2015  Maciej Komosinski and Szymon Ulatowski.
     2// Copyright (C) 1999-2019  Maciej Komosinski and Szymon Ulatowski.
    33// See LICENSE.txt for details.
    44
     
    88
    99/** Creates an empty matching for two objects of specified size.
    10         @param Obj0Size Size of the first object. Must be positive.
    11         @param Obj1Size Size of the second object. Must be positive.
     10                @param Obj0Size Size of the first object. Must be positive.
     11                @param Obj1Size Size of the second object. Must be positive.
    1212 */
    1313SimilMatching::SimilMatching(int Obj0Size, int Obj1Size)
    1414{
    15     // assure that sizes of objects are positive
    16     assert(Obj0Size > 0);
    17     assert(Obj1Size > 0);
    18 
    19     // create necessary vectors
    20     m_apvMatched[0] = new std::vector<int>(Obj0Size);
    21     m_apvMatched[1] = new std::vector<int>(Obj1Size);
    22 
    23     // assure that vectors are created
    24     assert(m_apvMatched[0] != NULL);
    25     assert(m_apvMatched[1] != NULL);
    26 
    27     // fill vectors with "unmatched" indicator
     15        // assure that sizes of objects are positive
     16        assert(Obj0Size > 0);
     17        assert(Obj1Size > 0);
     18
     19        // create necessary vectors
     20        m_apvMatched[0] = new std::vector<int>(Obj0Size);
     21        m_apvMatched[1] = new std::vector<int>(Obj1Size);
     22
     23        // assure that vectors are created
     24        assert(m_apvMatched[0] != NULL);
     25        assert(m_apvMatched[1] != NULL);
     26
     27        // fill vectors with "unmatched" indicator
    2828        for (unsigned int i = 0; i < m_apvMatched[0]->size(); i++)
    29     {
    30         m_apvMatched[0]->operator[](i) = -1;
    31     }
     29        {
     30                m_apvMatched[0]->operator[](i) = -1;
     31        }
    3232        for (unsigned int i = 0; i < m_apvMatched[1]->size(); i++)
    33     {
    34         m_apvMatched[1]->operator[](i) = -1;
    35     }
     33        {
     34                m_apvMatched[1]->operator[](i) = -1;
     35        }
    3636}
    3737
    3838/** A copying constructor.
    39         @param Source The object to be copied.
     39                @param Source The object to be copied.
    4040 */
    4141SimilMatching::SimilMatching(const SimilMatching &Source)
    4242{
    43     // copy the vectors of the actual matching
    44     m_apvMatched[ 0 ] = new std::vector<int>(* (Source.m_apvMatched[ 0 ]));
    45     m_apvMatched[ 1 ] = new std::vector<int>(* (Source.m_apvMatched[ 1 ]));
    46 
    47     // assure that vectors are created
    48     assert(m_apvMatched[0] != NULL);
    49     assert(m_apvMatched[1] != NULL);
     43        // copy the vectors of the actual matching
     44        m_apvMatched[0] = new std::vector<int>(*(Source.m_apvMatched[0]));
     45        m_apvMatched[1] = new std::vector<int>(*(Source.m_apvMatched[1]));
     46
     47        // assure that vectors are created
     48        assert(m_apvMatched[0] != NULL);
     49        assert(m_apvMatched[1] != NULL);
    5050}
    5151
     
    5454SimilMatching::~SimilMatching()
    5555{
    56     // delete vectors of matching
    57     delete m_apvMatched[0];
    58     delete m_apvMatched[1];
     56        // delete vectors of matching
     57        delete m_apvMatched[0];
     58        delete m_apvMatched[1];
    5959}
    6060
    6161/** Gets size of the specified object.
    62         @param Index of an object (must be 0 or 1).
    63         @return Size of the object (in elements).
     62                @param Index of an object (must be 0 or 1).
     63                @return Size of the object (in elements).
    6464 */
    6565int SimilMatching::GetObjectSize(int Obj)
    6666{
    67     // check parameter
    68     assert((Obj == 0) || (Obj == 1));
    69 
    70     // return the result
    71     return m_apvMatched[ Obj ]->size();
     67        // check parameter
     68        assert((Obj == 0) || (Obj == 1));
     69
     70        // return the result
     71        return m_apvMatched[Obj]->size();
    7272}
    7373
    7474/** Matches elements given by indices in the given objects.
    75         @param Obj0 Index of the first object. Must be 0 or 1.
    76         @param index0 Index of element in the first object. Must be a valid index
    77         ( >= 0 and < size of the object).
    78         @param Obj1 Index of the second object. Must be 0 or 1 and different from Obj0.
    79         @param Index1 index of element in the second object. Must be a valid index
    80         ( >= 0 and < size of the object).
     75                @param Obj0 Index of the first object. Must be 0 or 1.
     76                @param index0 Index of element in the first object. Must be a valid index
     77                ( >= 0 and < size of the object).
     78                @param Obj1 Index of the second object. Must be 0 or 1 and different from Obj0.
     79                @param Index1 index of element in the second object. Must be a valid index
     80                ( >= 0 and < size of the object).
    8181
    8282 */
    8383void SimilMatching::Match(int Obj0, int Index0, int Obj1, int Index1)
    8484{
    85     // check parameters of object 0
    86     assert((Obj0 == 0) || (Obj0 == 1));
    87     assert((Index0 >= 0) && (Index0 < (int) m_apvMatched[Obj0]->size()));
    88     // check parameters of object 1
    89     assert(((Obj1 == 0) || (Obj1 == 1)) && (Obj0 != Obj1));
    90     assert((Index1 >= 0) && (Index1 < (int) m_apvMatched[Obj1]->size()));
    91 
    92     // match given elements
    93     // matching_Obj0(Index0) = Index1
    94     m_apvMatched[ Obj0 ]->operator[](Index0) = Index1;
    95     // matching_Obj1(Index1) = Index0
    96     m_apvMatched[ Obj1 ]->operator[](Index1) = Index0;
     85        // check parameters of object 0
     86        assert((Obj0 == 0) || (Obj0 == 1));
     87        assert((Index0 >= 0) && (Index0 < (int)m_apvMatched[Obj0]->size()));
     88        // check parameters of object 1
     89        assert(((Obj1 == 0) || (Obj1 == 1)) && (Obj0 != Obj1));
     90        assert((Index1 >= 0) && (Index1 < (int)m_apvMatched[Obj1]->size()));
     91
     92        // match given elements
     93        // matching_Obj0(Index0) = Index1
     94        m_apvMatched[Obj0]->operator[](Index0) = Index1;
     95        // matching_Obj1(Index1) = Index0
     96        m_apvMatched[Obj1]->operator[](Index1) = Index0;
    9797}
    9898
    9999/** Checks if the given element in the given object is already matched.
    100         @param Obj Index of an object (must be 0 or 1).
    101         @param Index Index of an element in the given object. Must be a valid index
    102         ( >=0 and < size of the object).
    103         @return true if the given element is matched, false otherwise.
     100                @param Obj Index of an object (must be 0 or 1).
     101                @param Index Index of an element in the given object. Must be a valid index
     102                ( >=0 and < size of the object).
     103                @return true if the given element is matched, false otherwise.
    104104 */
    105105bool SimilMatching::IsMatched(int Obj, int Index)
    106106{
    107     // check parameters
    108     assert((Obj == 0) || (Obj == 1));
    109     assert((Index >= 0) && (Index < (int) m_apvMatched[ Obj ]->size()));
    110 
    111     // check if the element is matched
    112     if (m_apvMatched[ Obj ]->operator[](Index) >= 0)
    113     {
    114         return true;
    115     }
    116     else
    117     {
    118         return false;
    119     }
     107        // check parameters
     108        assert((Obj == 0) || (Obj == 1));
     109        assert((Index >= 0) && (Index < (int)m_apvMatched[Obj]->size()));
     110
     111        // check if the element is matched
     112        if (m_apvMatched[Obj]->operator[](Index) >= 0)
     113        {
     114                return true;
     115        }
     116        else
     117        {
     118                return false;
     119        }
    120120}
    121121
    122122/** Gets index of the element thet is matched in the other object withe the element given
    123         by parameters.
    124         @param Obj Index of an object (must be 0 or 1).
    125         @param Index Index of checked element in the given object.
    126         @return Index of the element (in the other organism) that is matched with the given
    127         element. WARNING! If the given element is not matched, the result may be smaller than 0
    128         (check IsMatched() before using GetMatchedIndex()).
     123                by parameters.
     124                @param Obj Index of an object (must be 0 or 1).
     125                @param Index Index of checked element in the given object.
     126                @return Index of the element (in the other organism) that is matched with the given
     127                element. WARNING! If the given element is not matched, the result may be smaller than 0
     128                (check IsMatched() before using GetMatchedIndex()).
    129129 */
    130130int SimilMatching::GetMatchedIndex(int Obj, int Index)
    131131{
    132     // check parameters
    133     assert((Obj == 0) || (Obj == 1));
    134     assert((Index >= 0) && (Index < (int) m_apvMatched[ Obj ]->size()));
    135 
    136     // return the index of the matched element
    137     return m_apvMatched[ Obj ]->operator[](Index);
    138 }
    139 
    140 /** Checks if the matching is already full, i.e. if the smaller object already has all its 
    141         elements matched.
    142         @return true if matching is full, false otherwise.
     132        // check parameters
     133        assert((Obj == 0) || (Obj == 1));
     134        assert((Index >= 0) && (Index < (int)m_apvMatched[Obj]->size()));
     135
     136        // return the index of the matched element
     137        return m_apvMatched[Obj]->operator[](Index);
     138}
     139
     140/** Checks if the matching is already full, i.e. if the smaller object already has all its
     141                elements matched.
     142                @return true if matching is full, false otherwise.
    143143 */
    144144bool SimilMatching::IsFull()
    145145{
    146     // assume that the matching is full
    147     bool bResult = true;
    148     // index of the smallest object
    149     int nObj;
    150 
    151     // find the smallest object (its index)
    152     if (m_apvMatched[ 0 ]->size() < m_apvMatched[ 1 ]->size())
    153     {
    154         nObj = 0;
    155     }
    156     else
    157     {
    158         nObj = 1;
    159     }
    160 
    161     // check if all elements of the smallest object are matched
     146        // assume that the matching is full
     147        bool bResult = true;
     148        // index of the smallest object
     149        int nObj;
     150
     151        // find the smallest object (its index)
     152        if (m_apvMatched[0]->size() < m_apvMatched[1]->size())
     153        {
     154                nObj = 0;
     155        }
     156        else
     157        {
     158                nObj = 1;
     159        }
     160
     161        // check if all elements of the smallest object are matched
    162162        for (unsigned int nElem = 0; nElem < m_apvMatched[nObj]->size(); nElem++)
    163     {
    164         if (m_apvMatched[ nObj ]->operator[](nElem) < 0)
    165         {
    166             // if any element is not matched, the result is false
    167             bResult = false;
    168             break;
    169         }
    170     }
    171 
    172     // return the result
    173     return bResult;
     163        {
     164                if (m_apvMatched[nObj]->operator[](nElem) < 0)
     165                {
     166                        // if any element is not matched, the result is false
     167                        bResult = false;
     168                        break;
     169                }
     170        }
     171
     172        // return the result
     173        return bResult;
    174174}
    175175
    176176/** Checks if the matching is empty (i.e. none of elements is matched).
    177     @return true if matching is empty, otherwise - false.
     177        @return true if matching is empty, otherwise - false.
    178178 */
    179179bool SimilMatching::IsEmpty()
    180180{
    181     // result - assume that matching is empty
    182     bool bResult = true;
    183 
    184     // matching is empty if either of objects has only unmatched elements
    185     // so it may be first object
    186     int nObj = 0;
     181        // result - assume that matching is empty
     182        bool bResult = true;
     183
     184        // matching is empty if either of objects has only unmatched elements
     185        // so it may be first object
     186        int nObj = 0;
    187187        for (unsigned int nElem = 0; nElem < m_apvMatched[nObj]->size(); nElem++)
    188     {
    189         if (m_apvMatched[ nObj ]->operator[](nElem) >= 0)
    190         {
    191             // if any element of the object is matched (unmatched objects have (-1))
    192             bResult = false;
    193             break;
    194         }
    195     }
    196 
    197     // return the result from the loop
    198     return bResult;
     188        {
     189                if (m_apvMatched[nObj]->operator[](nElem) >= 0)
     190                {
     191                        // if any element of the object is matched (unmatched objects have (-1))
     192                        bResult = false;
     193                        break;
     194                }
     195        }
     196
     197        // return the result from the loop
     198        return bResult;
    199199}
    200200
     
    203203void SimilMatching::Empty()
    204204{
    205     for (int iObj = 0; iObj < 2; iObj++) // a counter of objects
    206     {
    207         // for each object in the matching
     205        for (int iObj = 0; iObj < 2; iObj++) // a counter of objects
     206        {
     207                // for each object in the matching
    208208                for (unsigned int iElem = 0; iElem < m_apvMatched[iObj]->size(); iElem++) // a counter of objects' elements
    209         {
    210             // for each element iElem for the object iObj
    211             // set it as unmatched (marker: -1)
    212             m_apvMatched[ iObj ]->operator[](iElem) = -1;
    213         }
    214     }
    215 
    216     // the exit condition
    217     assert(IsEmpty() == true);
     209                {
     210                        // for each element iElem for the object iObj
     211                        // set it as unmatched (marker: -1)
     212                        m_apvMatched[iObj]->operator[](iElem) = -1;
     213                }
     214        }
     215
     216        // the exit condition
     217        assert(IsEmpty() == true);
    218218}
    219219
     
    222222void SimilMatching::PrintMatching()
    223223{
    224     int nBigger;
    225 
    226     // check which object is bigger
    227     if (m_apvMatched[ 0 ]->size() >= m_apvMatched[ 1 ]->size())
    228     {
    229         nBigger = 0;
    230     }
    231     else
    232     {
    233         nBigger = 1;
    234     }
    235 
    236     // print first line - indices of objects
    237     printf("[ ");
     224        int nBigger;
     225
     226        // check which object is bigger
     227        if (m_apvMatched[0]->size() >= m_apvMatched[1]->size())
     228        {
     229                nBigger = 0;
     230        }
     231        else
     232        {
     233                nBigger = 1;
     234        }
     235
     236        // print first line - indices of objects
     237        printf("[ ");
    238238        for (unsigned int i = 0; i < m_apvMatched[nBigger]->size(); i++)
    239     {
    240         printf("%2d ", i);
    241     }
    242     printf("]\n");
    243 
    244     // print second line and third - indices of elements matched with elements of the objects
    245     for (int nObj = 0; nObj < 2; nObj++)
    246     {
    247         // for both objects - print out lines of matched elements
    248         printf("[ ");
     239        {
     240                printf("%2d ", i);
     241        }
     242        printf("]\n");
     243
     244        // print second line and third - indices of elements matched with elements of the objects
     245        for (int nObj = 0; nObj < 2; nObj++)
     246        {
     247                // for both objects - print out lines of matched elements
     248                printf("[ ");
    249249                for (unsigned int i = 0; i < m_apvMatched[nObj]->size(); i++)
    250         {
    251             if (IsMatched(nObj, i))
    252             {
    253                 // if the element is matched - print the index
    254                 printf("%2d ", GetMatchedIndex(nObj, i));
    255             }
    256             else
    257             {
    258                 // if the element is not matched - print "X"
    259                 printf(" X ");
    260             }
    261         }
    262         printf("]\n");
    263     }
    264 }
     250                {
     251                        if (IsMatched(nObj, i))
     252                        {
     253                                // if the element is matched - print the index
     254                                printf("%2d ", GetMatchedIndex(nObj, i));
     255                        }
     256                        else
     257                        {
     258                                // if the element is not matched - print "X"
     259                                printf(" X ");
     260                        }
     261                }
     262                printf("]\n");
     263        }
     264}
  • cpp/frams/model/similarity/simil_match.h

    r357 r869  
    11// This file is a part of Framsticks SDK.  http://www.framsticks.com/
    2 // Copyright (C) 1999-2015  Maciej Komosinski and Szymon Ulatowski.
     2// Copyright (C) 1999-2019  Maciej Komosinski and Szymon Ulatowski.
    33// See LICENSE.txt for details.
    44
     
    88#include <vector>
    99
    10 /** Class describes a mutually single-valued function between two sets of elements
    11         (these sets are called objects). These sets may have different sizes, so this function
    12         is mutually single-valued only for some subset of the bigger set.
    13         This class is used while building a matching function between Parts of two Framsticks'
    14         organisms (similarity measure computation).
     10/** This class describes a mutually single-valued function between two sets of elements
     11        (these sets are called objects). These sets may have different sizes, so this function
     12        is mutually single-valued only for some subset of the bigger set.
     13        This class is used while building a matching function between Parts of two Framsticks'
     14        organisms (similarity measure computation).
    1515 */
    1616class SimilMatching
    1717{
    1818protected:
    19     /** Array of pointers to two vectors. Each one stores indices of matched elements of
    20             the other object. Value <0 means that an element is not matched yet. Sizes of these
    21             vectors are sizes of objects themselves.
    22     */
    23     std::vector<int> *m_apvMatched[2];
     19        /** Array of pointers to two vectors. Each one stores indices of matched elements of
     20                the other object. Value <0 means that an element is not matched yet. Sizes of these
     21                vectors are sizes of objects themselves.
     22        */
     23        std::vector<int> *m_apvMatched[2];
    2424public:
    25     SimilMatching(int Obj0Size, int Obj1Size);
    26     SimilMatching(const SimilMatching &Source);
    27     ~SimilMatching();
    28     int GetObjectSize(int Obj);
    29     void Match(int Obj0, int Index0, int Obj1, int Index1);
    30     bool IsMatched(int Obj, int Index);
    31     int GetMatchedIndex(int Obj, int index);
    32     bool IsFull();
    33     bool IsEmpty();
    34     void Empty();
    35     void PrintMatching();
     25        SimilMatching(int Obj0Size, int Obj1Size);
     26        SimilMatching(const SimilMatching &Source);
     27        ~SimilMatching();
     28        int GetObjectSize(int Obj);
     29        void Match(int Obj0, int Index0, int Obj1, int Index1);
     30        bool IsMatched(int Obj, int Index);
     31        int GetMatchedIndex(int Obj, int index);
     32        bool IsFull();
     33        bool IsEmpty();
     34        void Empty();
     35        void PrintMatching();
    3636};
    3737
  • cpp/frams/model/similarity/simil_model.cpp

    r818 r869  
    11// This file is a part of Framsticks SDK.  http://www.framsticks.com/
    2 // Copyright (C) 1999-2017  Maciej Komosinski and Szymon Ulatowski.
     2// Copyright (C) 1999-2019  Maciej Komosinski and Szymon Ulatowski.
    33// See LICENSE.txt for details.
    44
    5 // simil_model.cpp: implementation of the ModelSimil class.
    6 //
    7 //////////////////////////////////////////////////////////////////////
     5// implementation of the ModelSimil class.
     6
    87#include "SVD/matrix_tools.h"
     8#include "hungarian/hungarian.h"
    99#include "simil_model.h"
    1010#include "simil_match.h"
     
    3535
    3636static ParamEntry MSparam_tab[] = {
    37                 { "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", },
     37                { "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", },
     38                { "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), "",},
    3839                { "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.", },
    3940                { "simil_partdeg", 0, 0, "Weight of parts' degree", "f 0 100 1", FIELD(m_adFactors[1]), "", },
     
    7677        isFuzzy = 0;
    7778        fuzzyDepth = 10;
    78        
     79
    7980        //Determines whether weighted MDS should be used.
    8081        wMDS = 0;
    81 }
    82 
    83 /**     Evaluates distance between two given genotypes. The distance depends strongly
    84         on weights set.
    85         @param G0 Pointer to the first of compared genotypes
    86         @param G1 Pointer to the second of compared genotypes.
    87         @return Distance between two genotypes.
    88         @sa m_adFactors, matching_method
    89         */
    90 double ModelSimil::EvaluateDistance(const Geno *G0, const Geno *G1)
     82        //Determines whether best matching should be saved using hungarian similarity measure.
     83        saveMatching = 0;
     84}
     85
     86double ModelSimil::EvaluateDistanceGreedy(const Geno *G0, const Geno *G1)
    9187{
    9288        double dResult = 0;
     
    9894        if (m_Gen[0] == NULL || m_Gen[1] == NULL)
    9995        {
    100                 DB(printf("ModelSimil::EvaluateDistance - invalid genotypes pointers\n");)
    101                         return 0.0;
     96                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.
     97                        logPrintf("ModelSimil", "EvaluateDistanceGreedy", LOG_ERROR, "NULL genotype pointer(s)");
     98                return 0.0;
    10299        }
    103100        // create models of objects to compare
     
    108105        if (m_Mod[0] == NULL || m_Mod[1] == NULL || !(m_Mod[0]->isValid()) || !(m_Mod[1]->isValid()))
    109106        {
    110                 DB(printf("ModelSimil::EvaluateDistance - invalid models pointers\n");)
     107                DB(printf("ModelSimil::EvaluateDistanceGreedy - invalid model(s) pointers\n");)
    111108                        return 0.0;
    112109        }
     
    143140        if ((this->*pfMatchingFunction)() == 0)
    144141        {
    145                 DB(printf("ModelSimil::EvaluateDistance - MatchParts() error\n");)
     142                DB(printf("ModelSimil::EvaluateDistanceGreedy - MatchParts() error\n");)
    146143                        return 0.0;
    147144        }
     
    155152                if (CountPartsDistance() == 0)
    156153                {
    157                 DB(printf("ModelSimil::EvaluateDistance - CountPartDistance() error\n");)
    158                         return 0.0;
     154                        DB(printf("ModelSimil::EvaluateDistanceGreedy - CountPartDistance() error\n");)
     155                                return 0.0;
    159156                }
    160157
     
    340337                if (tdn1[1] < tdn2[1])
    341338                {
    342                 // when degree - tdn2[1] - is BIGGER
    343                 return 1;
     339                        // when degree - tdn2[1] - is BIGGER
     340                        return 1;
     341                }
     342                else
     343                {
     344                        return 0;
     345                }
     346}
     347
     348/** Comparison function required for qsort() call. Used while sorting groups of
     349        Parts with respect to fuzzy vertex degree. Compares two TDN structures
     350        with respect to their [4] field ( fuzzy vertex degree). Highest degree goes first.
     351        @param pElem1 Pointer to the TDN structure of some Part.
     352        @param pElem2 Pointer to the TDN structure of some Part.
     353        @return (-1) - pElem1 should go first, 0 - equal, (1) - pElem2 should go first.
     354        */
     355int ModelSimil::CompareFuzzyDegrees(const void *pElem1, const void *pElem2)
     356{
     357        int *tdn1 = (int *)pElem1;
     358        int *tdn2 = (int *)pElem2;
     359
     360        if (tdn1[4] > tdn2[4])
     361        {
     362                // when degree - tdn1[4] - is BIGGER
     363                return -1;
     364        }
     365        else
     366                if (tdn1[4] < tdn2[4])
     367                {
     368                        // when degree - tdn2[4] - is BIGGER
     369                        return 1;
    344370                }
    345371                else
     
    374400                if (tdn1[NEURO_CONNS] < tdn2[NEURO_CONNS])
    375401                {
    376                 // when number of NConn Elem1 is SMALLER
    377                 return 1;
     402                        // when number of NConn Elem1 is SMALLER
     403                        return 1;
    378404                }
    379405                else
     
    389415                                if (tdn1[NEURONS] < tdn2[NEURONS])
    390416                                {
    391                                 // when number of Neu is SMALLER for Elem1
    392                                 return 1;
     417                                        // when number of Neu is SMALLER for Elem1
     418                                        return 1;
    393419                                }
    394420                                else
     
    406432                                                if (tdn1[ORIG_IND] < tdn2[ORIG_IND])
    407433                                                {
    408                                                 // when the number of NIt Deg1 id SMALLER
    409                                                 return 1;
     434                                                        // when the number of NIt Deg1 id SMALLER
     435                                                        return 1;
    410436                                                }
    411437                                                else
     
    494520}
    495521
     522/** Prints one array of parts fuzzy neighbourhood.
     523        @param index of organism
     524        */
     525void ModelSimil::_PrintFuzzyNeighbourhood(int o)
     526{
     527        assert(m_fuzzyNeighb[o] != NULL);
     528        printf("Fuzzy neighbourhhod of organism %i\n", o);
     529        int size = m_Mod[o]->getPartCount();
     530        for (int i = 0; i < size; i++)
     531        {
     532                for (int j = 0; j < fuzzyDepth; j++)
     533                {
     534                        printf("%f ", m_fuzzyNeighb[o][i][j]);
     535                }
     536                printf("\n");
     537        }
     538}
     539
    496540/** Creates arrays holding information about organisms' Parts (m_aDegrees) andm_Neigh
    497541        fills them with initial data (original indices and zeros).
     
    526570                                for (j = 0; j < partCount; j++)
    527571                                {
    528                                 m_aDegrees[i][j][0] = j;
    529                                 m_aDegrees[i][j][1] = 0;
    530                                 m_aDegrees[i][j][2] = 0;
    531                                 m_aDegrees[i][j][3] = 0;
    532                                 m_aDegrees[i][j][4] = 0;
    533 
    534                                 // sprawdz, czy nie piszemy po jakims szalonym miejscu pamieci
    535                                 assert(m_aDegrees[i][j] != NULL);
    536 
    537                                 if (isFuzzy)
    538                                 {
    539                                         m_Neighbours[i][j] = new int[partCount];
    540                                         for (int k = 0; k < partCount; k++)
     572                                        m_aDegrees[i][j][0] = j;
     573                                        m_aDegrees[i][j][1] = 0;
     574                                        m_aDegrees[i][j][2] = 0;
     575                                        m_aDegrees[i][j][3] = 0;
     576                                        m_aDegrees[i][j][4] = 0;
     577
     578                                        // sprawdz, czy nie piszemy po jakims szalonym miejscu pamieci
     579                                        assert(m_aDegrees[i][j] != NULL);
     580
     581                                        if (isFuzzy)
    541582                                        {
    542                                                 m_Neighbours[i][j][k] = 0;
     583                                                m_Neighbours[i][j] = new int[partCount];
     584                                                for (int k = 0; k < partCount; k++)
     585                                                {
     586                                                        m_Neighbours[i][j][k] = 0;
     587                                                }
     588
     589                                                m_fuzzyNeighb[i][j] = new float[fuzzyDepth];
     590                                                for (int k = 0; k < fuzzyDepth; k++)
     591                                                {
     592                                                        m_fuzzyNeighb[i][j][k] = 0;
     593                                                }
     594
     595                                                assert(m_Neighbours[i][j] != NULL);
     596                                                assert(m_fuzzyNeighb[i][j] != NULL);
    543597                                        }
    544 
    545                                         m_fuzzyNeighb[i][j] = new float[fuzzyDepth];
    546                                         for (int k = 0; k < fuzzyDepth; k++)
    547                                         {
    548                                                 m_fuzzyNeighb[i][j][k] = 0;
    549                                         }
    550 
    551                                         assert(m_Neighbours[i][j] != NULL);
    552                                         assert(m_fuzzyNeighb[i][j] != NULL);
    553                                 }
    554598
    555599                                }
     
    671715}
    672716
    673 //store information about identity of parts "fuzzy degrees" in the m_aDegrees[4]
    674 void ModelSimil::CheckFuzzyIdentity()
    675 {
    676         int partCount = 0;
     717void ModelSimil::FuzzyOrder()
     718{
     719        int i, depth, partInd, prevPartInd, partCount;
    677720        for (int mod = 0; mod < 2; mod++)
    678721        {
    679                 //for subsequent pairs of parts
    680722                partCount = m_Mod[mod]->getPartCount();
    681                 m_aDegrees[mod][partCount - 1][FUZZ_DEG] = 0;
    682                 for (int partInd = (partCount - 2); partInd >= 0; partInd--)
    683                 {
    684                         m_aDegrees[mod][partInd][FUZZ_DEG] = m_aDegrees[mod][partInd + 1][FUZZ_DEG];
    685                         for (int depth = 1; depth < fuzzyDepth; depth++)
    686                         {
    687                                 if (m_fuzzyNeighb[mod][partInd][depth] != m_fuzzyNeighb[mod][partInd + 1][depth])
    688                                 {
    689                                         m_aDegrees[mod][partInd][FUZZ_DEG] += 1;
     723                partInd = m_fuzzyNeighb[mod][partCount - 1][0];
     724                m_aDegrees[mod][partInd][FUZZ_DEG] = 0;
     725
     726                for (i = (partCount - 2); i >= 0; i--)
     727                {
     728                        prevPartInd = partInd;
     729                        partInd = m_fuzzyNeighb[mod][i][0];
     730                        m_aDegrees[mod][partInd][FUZZ_DEG] = m_aDegrees[mod][prevPartInd][FUZZ_DEG];
     731                        for (depth = 1; depth < fuzzyDepth; depth++)
     732                        {
     733                                if (m_fuzzyNeighb[mod][i][depth] != m_fuzzyNeighb[mod][i + 1][depth])
     734                                {
     735                                        m_aDegrees[mod][partInd][FUZZ_DEG]++;
    690736                                        break;
    691737                                }
     
    761807
    762808        SortFuzzyNeighb();
     809        FuzzyOrder();
    763810}
    764811
     
    884931
    885932        int i;
     933        int(*pfDegreeFunction) (const void*, const void*) = NULL;
     934        pfDegreeFunction = (isFuzzy == 1) ? &CompareFuzzyDegrees : &CompareDegrees;
    886935        // sortowanie obu tablic wg stopni punktów - TDN[1]
    887         if (isFuzzy != 1)
    888         {
    889                 for (i = 0; i < 2; i++)
    890                 {
    891                         DB(_PrintDegrees(i));
    892                         std::qsort(m_aDegrees[i], (size_t)(m_Mod[i]->getPartCount()),
    893                                 sizeof(TDN), ModelSimil::CompareDegrees);
    894                         DB(_PrintDegrees(i));
    895                 }
    896         }//sortowanie wg romzmytego stopnia wierzcholka
    897 
    898         else
    899         {
    900                 SortPartInfoFuzzy();
    901         }
    902 
     936        for (i = 0; i < 2; i++)
     937        {
     938                DB(_PrintDegrees(i));
     939                std::qsort(m_aDegrees[i], (size_t)(m_Mod[i]->getPartCount()),
     940                        sizeof(TDN), pfDegreeFunction);
     941                DB(_PrintDegrees(i));
     942        }
    903943
    904944        // sprawdzenie wartosci parametru
     
    935975
    936976                                        std::qsort(m_aDegrees[i][iPocz], (size_t)(j - iPocz),
    937                                         sizeof(TDN), ModelSimil::CompareConnsNo);
     977                                                sizeof(TDN), ModelSimil::CompareConnsNo);
    938978                                DB(_PrintDegrees(i));
    939979                                // wyswietlamy z jedna komorka po zakonczeniu tablicy
     
    948988}
    949989
    950 int ModelSimil::SortPartInfoFuzzy()
    951 {
    952 
    953         // sprawdz zalozenie - o modelach
    954         assert((m_Mod[0] != NULL) && (m_Mod[1] != NULL));
    955         assert(m_Mod[0]->isValid() && m_Mod[1]->isValid());
    956 
    957         // sprawdz zalozenie - o tablicach
    958         assert(m_aDegrees[0] != NULL);
    959         assert(m_aDegrees[1] != NULL);
    960         // sprawdz zalozenie - o tablicach
    961         assert(m_fuzzyNeighb[0] != NULL);
    962         assert(m_fuzzyNeighb[1] != NULL);
    963 
    964 
    965         TDN * m_aDegreesTmp[2];
    966 
    967         for (int i = 0; i < 2; i++)
    968         {
    969                 int partCount = m_Mod[i]->getPartCount();
    970                 m_aDegreesTmp[i] = new TDN[partCount];
    971 
    972                 for (int j = 0; j < partCount; j++)
    973                 {
    974                         for (int k = 0; k < TDN_SIZE; k++)
    975                         {
    976                                 m_aDegreesTmp[i][j][k] = m_aDegrees[i][j][k];
    977                         }
    978                 }
    979         }
    980 
    981         int newInd = 0;
    982         for (int i = 0; i < 2; i++)
    983         {
    984                 int partCount = m_Mod[i]->getPartCount();
    985                 for (int j = 0; j < partCount; j++)
    986                 {
    987                         newInd = (int)m_fuzzyNeighb[i][j][0];
    988                         for (int k = 0; k < TDN_SIZE; k++)
    989                         {
    990                                 m_aDegrees[i][j][k] = m_aDegreesTmp[i][newInd][k];
    991                         }
    992                 }
    993         }
    994 
    995         SAFEDELETEARRAY(m_aDegreesTmp[0]);
    996         SAFEDELETEARRAY(m_aDegreesTmp[1]);
    997 
    998         CheckFuzzyIdentity();
    999 
    1000         return 1;
    1001 }
    1002 
    1003 /** Checks if given Parts have identical physical and biological properties
    1004         (except for geometry that might differ).
    1005         @param P1 Pointer to first Part.
    1006         @param P2 Pointer to second Part.
    1007         @return 1 - identical properties, 0 - there are differences.
    1008         */
    1009 int ModelSimil::CheckPartsIdentity(Part *P1, Part *P2)
    1010 {
    1011         // sprawdz, czy te Parts chociaz sa w sensownym miejscu pamieci
    1012         assert((P1 != NULL) && (P2 != NULL));
    1013 
    1014         if ((P1->assim != P2->assim) ||
    1015                 (P1->friction != P2->friction) ||
    1016                 (P1->ingest != P2->ingest) ||
    1017                 (P1->mass != P2->mass) ||
    1018                 (P1->size != P2->size) ||
    1019                 (P1->density != P2->density))
    1020                 // gdy znaleziono jakas roznice w parametrach fizycznych i
    1021                 // biologicznych
    1022                 return 0;
    1023         else
    1024                 // gdy nie ma roznic
    1025                 return 1;
    1026 }
    1027990
    1028991/** Prints the state of the matching object. Debug method.
     
    11581121                        aiKoniecGrupyDopasowania[0]);)
    11591122                        DB(printf("Organizm 1: grupa: (%2i, %2i)\n", aiKoniecPierwszejGrupy[1],
    1160                         aiKoniecGrupyDopasowania[1]);)
     1123                                aiKoniecGrupyDopasowania[1]);)
    11611124
    11621125                        for (i = 0; i < aiRozmiarCalychGrup[0]; i++)
    11631126                        {
    11641127
    1165                         // iterujemy i - Parts organizmu 0
    1166                         // (indeks podstawowy - aiKoniecPierwszejGrupy[0])
    1167 
    1168                         if (!(m_pMatching->IsMatched(0, aiKoniecPierwszejGrupy[0] + i)))
    1169                         {
    1170                                 // interesuja nas tylko te niedopasowane jeszcze (i)
    1171                                 for (j = 0; j < aiRozmiarCalychGrup[1]; j++)
    1172                                 {
    1173 
    1174                                         // iterujemy j - Parts organizmu 1
    1175                                         // (indeks podstawowy - aiKoniecPierwszejGrupy[1])
    1176 
    1177                                         if (!(m_pMatching->IsMatched(1, aiKoniecPierwszejGrupy[1] + j)))
     1128                                // iterujemy i - Parts organizmu 0
     1129                                // (indeks podstawowy - aiKoniecPierwszejGrupy[0])
     1130
     1131                                if (!(m_pMatching->IsMatched(0, aiKoniecPierwszejGrupy[0] + i)))
     1132                                {
     1133                                        // interesuja nas tylko te niedopasowane jeszcze (i)
     1134                                        for (j = 0; j < aiRozmiarCalychGrup[1]; j++)
    11781135                                        {
    1179                                                 // interesuja nas tylko te niedopasowane jeszcze (j)
    1180                                                 // teraz obliczymy lokalne różnice pomiędzy Parts
    1181                                                 iDeg = abs(m_aDegrees[0][aiKoniecPierwszejGrupy[0] + i][1]
    1182                                                         - m_aDegrees[1][aiKoniecPierwszejGrupy[1] + j][1]);
    1183 
    1184                                                 //iNit currently is not a component of distance measure           
    1185                                                 //iNIt = abs(m_aDegrees[0][ aiKoniecPierwszejGrupy[0] + i ][2]
    1186                                                 //           - m_aDegrees[1][ aiKoniecPierwszejGrupy[1] + j ][2]);
    1187 
    1188                                                 iNeu = abs(m_aDegrees[0][aiKoniecPierwszejGrupy[0] + i][3]
    1189                                                         - m_aDegrees[1][aiKoniecPierwszejGrupy[1] + j][3]);
    1190 
    1191                                                 // obliczenie także lokalnych różnic geometrycznych pomiędzy Parts
    1192                                                 // find original indices of Parts for both of the models
    1193                                                 iPartIndex[0] = m_aDegrees[0][aiKoniecPierwszejGrupy[0] + i][0];
    1194                                                 iPartIndex[1] = m_aDegrees[1][aiKoniecPierwszejGrupy[1] + j][0];
    1195                                                 // now compute the geometrical distance of these Parts (use m_aPositions
    1196                                                 // which should be computed by SVD)
    1197                                                 Pt3D Part0Pos(m_aPositions[0][iPartIndex[0]]);
    1198                                                 Pt3D Part1Pos(m_aPositions[1][iPartIndex[1]]);
    1199                                                 dGeo = m_adFactors[3] == 0 ? 0 : Part0Pos.distanceTo(Part1Pos); //no need to compute distane when m_dDG weight is 0
    1200 
    1201                                                 // tutaj prawdopodobnie należy jeszcze dodać sprawdzanie
    1202                                                 // identyczności pozostałych własności (oprócz geometrii)
    1203                                                 // - żeby móc stwierdzić w ogóle identyczność Parts
    1204 
    1205                                                 // i ostateczna odleglosc indukowana przez te roznice
    1206                                                 // (tutaj nie ma różnicy w liczbie wszystkich wierzchołków)
    1207                                                 aadOdleglosciParts[i][j] = m_adFactors[1] * double(iDeg) +
    1208                                                         m_adFactors[2] * double(iNeu) +
    1209                                                         m_adFactors[3] * dGeo;
    1210                                                 DB(printf("Parts Distance (%2i,%2i) = %.3lf\n", aiKoniecPierwszejGrupy[0] + i,
    1211                                                         aiKoniecPierwszejGrupy[1] + j, aadOdleglosciParts[i][j]);)
    1212                                                         DB(printf("Parts geometrical distance = %.20lf\n", dGeo);)
    1213                                                         DB(printf("Pos0: (%.3lf %.3lf %.3lf)\n", Part0Pos.x, Part0Pos.y, Part0Pos.z);)
    1214                                                         DB(printf("Pos1: (%.3lf %.3lf %.3lf)\n", Part1Pos.x, Part1Pos.y, Part1Pos.z);)
     1136
     1137                                                // iterujemy j - Parts organizmu 1
     1138                                                // (indeks podstawowy - aiKoniecPierwszejGrupy[1])
     1139
     1140                                                if (!(m_pMatching->IsMatched(1, aiKoniecPierwszejGrupy[1] + j)))
     1141                                                {
     1142                                                        // interesuja nas tylko te niedopasowane jeszcze (j)
     1143                                                        // teraz obliczymy lokalne różnice pomiędzy Parts
     1144                                                        iDeg = abs(m_aDegrees[0][aiKoniecPierwszejGrupy[0] + i][1]
     1145                                                                - m_aDegrees[1][aiKoniecPierwszejGrupy[1] + j][1]);
     1146
     1147                                                        //iNit currently is not a component of distance measure           
     1148                                                        //iNIt = abs(m_aDegrees[0][ aiKoniecPierwszejGrupy[0] + i ][2]
     1149                                                        //           - m_aDegrees[1][ aiKoniecPierwszejGrupy[1] + j ][2]);
     1150
     1151                                                        iNeu = abs(m_aDegrees[0][aiKoniecPierwszejGrupy[0] + i][3]
     1152                                                                - m_aDegrees[1][aiKoniecPierwszejGrupy[1] + j][3]);
     1153
     1154                                                        // obliczenie także lokalnych różnic geometrycznych pomiędzy Parts
     1155                                                        // find original indices of Parts for both of the models
     1156                                                        iPartIndex[0] = m_aDegrees[0][aiKoniecPierwszejGrupy[0] + i][0];
     1157                                                        iPartIndex[1] = m_aDegrees[1][aiKoniecPierwszejGrupy[1] + j][0];
     1158                                                        // now compute the geometrical distance of these Parts (use m_aPositions
     1159                                                        // which should be computed by SVD)
     1160                                                        Pt3D Part0Pos(m_aPositions[0][iPartIndex[0]]);
     1161                                                        Pt3D Part1Pos(m_aPositions[1][iPartIndex[1]]);
     1162                                                        dGeo = m_adFactors[3] == 0 ? 0 : Part0Pos.distanceTo(Part1Pos); //no need to compute distane when m_dDG weight is 0
     1163
     1164                                                        // tutaj prawdopodobnie należy jeszcze dodać sprawdzanie
     1165                                                        // identyczności pozostałych własności (oprócz geometrii)
     1166                                                        // - żeby móc stwierdzić w ogóle identyczność Parts
     1167
     1168                                                        // i ostateczna odleglosc indukowana przez te roznice
     1169                                                        // (tutaj nie ma różnicy w liczbie wszystkich wierzchołków)
     1170                                                        aadOdleglosciParts[i][j] = m_adFactors[1] * double(iDeg) +
     1171                                                                m_adFactors[2] * double(iNeu) +
     1172                                                                m_adFactors[3] * dGeo;
     1173                                                        DB(printf("Parts Distance (%2i,%2i) = %.3lf\n", aiKoniecPierwszejGrupy[0] + i,
     1174                                                                aiKoniecPierwszejGrupy[1] + j, aadOdleglosciParts[i][j]);)
     1175                                                                DB(printf("Parts geometrical distance = %.20lf\n", dGeo);)
     1176                                                                DB(printf("Pos0: (%.3lf %.3lf %.3lf)\n", Part0Pos.x, Part0Pos.y, Part0Pos.z);)
     1177                                                                DB(printf("Pos1: (%.3lf %.3lf %.3lf)\n", Part1Pos.x, Part1Pos.y, Part1Pos.z);)
     1178                                                }
    12151179                                        }
    12161180                                }
    1217                         }
    12181181                        }
    12191182
     
    13881351                                        if (PartZ1NaLiscie0 || PartZ0NaLiscie1)
    13891352                                        {
    1390                                         // PRZYPADEK 2. Tylko jeden z Parts ma drugiego na swojej liscie
    1391                                         // mozliwych dopasowan
    1392                                         // AKCJA. Dopasowanie jednego jest proste (tego, ktory nie ma
    1393                                         // na swojej liscie drugiego). Dla tego drugiego nalezy powtorzyc
    1394                                         // duza czesc obliczen (znalezc mu nowa mozliwa pare)
    1395 
    1396                                         // indeks organizmu, ktorego Part nie ma dopasowywanego Part
    1397                                         // z drugiego organizmu na swojej liscie
    1398                                         int iBezDrugiego;
    1399 
    1400                                         // okreslenie indeksu organizmu z dopasowywanym Part
    1401                                         if (!PartZ1NaLiscie0)
    1402                                         {
    1403                                                 iBezDrugiego = 0;
    1404                                         }
    1405                                         else
    1406                                         {
    1407                                                 iBezDrugiego = 1;
    1408                                         }
    1409                                         // sprawdz indeks organizmu
    1410                                         assert((iBezDrugiego == 0) || (iBezDrugiego == 1));
    1411 
    1412                                         int iDopasowywany = -1;
    1413                                         // poszukujemy pierwszego z listy dopasowania
    1414                                         for (i = 0; i < aiRozmiarCalychGrup[1 - iBezDrugiego]; i++)
    1415                                         {
    1416                                                 if (apvbCzyMinimalnaOdleglosc[iBezDrugiego]->operator[](i))
     1353                                                // PRZYPADEK 2. Tylko jeden z Parts ma drugiego na swojej liscie
     1354                                                // mozliwych dopasowan
     1355                                                // AKCJA. Dopasowanie jednego jest proste (tego, ktory nie ma
     1356                                                // na swojej liscie drugiego). Dla tego drugiego nalezy powtorzyc
     1357                                                // duza czesc obliczen (znalezc mu nowa mozliwa pare)
     1358
     1359                                                // indeks organizmu, ktorego Part nie ma dopasowywanego Part
     1360                                                // z drugiego organizmu na swojej liscie
     1361                                                int iBezDrugiego;
     1362
     1363                                                // okreslenie indeksu organizmu z dopasowywanym Part
     1364                                                if (!PartZ1NaLiscie0)
    14171365                                                {
    1418                                                         iDopasowywany = i;
    1419                                                         break;
     1366                                                        iBezDrugiego = 0;
    14201367                                                }
    1421                                         }
    1422                                         // sprawdz poprawnosc indeksu dopasowywanego (musimy cos znalezc!)
    1423                                         // nieujemny i w odpowiedniej grupie!
    1424                                         assert((iDopasowywany >= 0) &&
    1425                                                 (iDopasowywany < aiRozmiarCalychGrup[1 - iBezDrugiego]));
    1426 
    1427                                         // znalezlismy 1. Part z listy dopasowania - dopasowujemy!
    1428                                         m_pMatching->Match(
    1429                                                 iBezDrugiego,
    1430                                                 aiKoniecPierwszejGrupy[iBezDrugiego] + iIndex[iBezDrugiego],
    1431                                                 1 - iBezDrugiego,
    1432                                                 aiKoniecPierwszejGrupy[1 - iBezDrugiego] + iDopasowywany);
    1433                                         DB(printf("Przypadek 2.1.: dopasowane Parts dopasowanie bez drugiego: (%2i, %2i)\n", aiKoniecPierwszejGrupy[iBezDrugiego] + iIndex[iBezDrugiego],
    1434                                                 aiKoniecPierwszejGrupy[1 - iBezDrugiego] + iDopasowywany);)
    1435 
    1436                                                 // zmniejsz liczby niedopasowanych elementow w grupach
    1437                                                 aiRozmiarGrupy[0]--;
    1438                                         aiRozmiarGrupy[1]--;
    1439 
    1440                                         // sprawdz, czy jest szansa dopasowania tego Part z drugiej strony
    1441                                         // (ktora miala mozliwosc dopasowania tego Part z poprzedniego organizmu)
    1442                                         if ((aiRozmiarGrupy[0] > 0) && (aiRozmiarGrupy[1] > 0))
    1443                                         {
    1444                                                 // jesli grupy sie jeszcze nie wyczrpaly
    1445                                                 // to jest mozliwosc dopasowania w organizmie
    1446 
    1447                                                 int iZDrugim = 1 - iBezDrugiego;
     1368                                                else
     1369                                                {
     1370                                                        iBezDrugiego = 1;
     1371                                                }
    14481372                                                // sprawdz indeks organizmu
    1449                                                 assert((iZDrugim == 0) || (iZDrugim == 1));
    1450 
    1451                                                 // bedziemy szukac minimum wsrod niedopasowanych - musimy wykasowac
    1452                                                 // poprzednie obliczenia minimum
    1453                                                 // dla organizmu 1 (o rozmiarze grupy z 0)
    1454                                                 for (i = 0; i < aiRozmiarCalychGrup[0]; i++)
     1373                                                assert((iBezDrugiego == 0) || (iBezDrugiego == 1));
     1374
     1375                                                int iDopasowywany = -1;
     1376                                                // poszukujemy pierwszego z listy dopasowania
     1377                                                for (i = 0; i < aiRozmiarCalychGrup[1 - iBezDrugiego]; i++)
    14551378                                                {
    1456                                                         apvbCzyMinimalnaOdleglosc[1]->operator[](i) = false;
    1457                                                 }
    1458                                                 // dla organizmu 0 (o rozmiarze grup z 1)
    1459                                                 for (i = 0; i < aiRozmiarCalychGrup[1]; i++)
    1460                                                 {
    1461                                                         apvbCzyMinimalnaOdleglosc[0]->operator[](i) = false;
    1462                                                 }
    1463 
    1464                                                 // szukamy na nowo minimum dla Part o indeksie iIndex[ iZDrugim ] w organizmie iZDrugim
    1465                                                 // wsrod niedopasowanych Parts z organizmu 1 - iZDrugim
    1466                                                 dMinimum = HUGE_VAL;
    1467                                                 for (i = 0; i < aiRozmiarCalychGrup[1 - iZDrugim]; i++)
    1468                                                 {
    1469                                                         if (!(m_pMatching->IsMatched(
    1470                                                                 1 - iZDrugim,
    1471                                                                 aiKoniecPierwszejGrupy[1 - iZDrugim] + i)))
    1472                                                         {
    1473                                                                 // szukamy minimum obliczonej lokalnej odleglosci tylko wsrod
    1474                                                                 // niedopasowanych jeszcze Parts
    1475                                                                 if (iZDrugim == 0)
    1476                                                                 {
    1477                                                                         // teraz niestety musimy rozpoznac odpowiedni organizm
    1478                                                                         // zeby moc indeksowac niesymetryczna tablice
    1479                                                                         if (aadOdleglosciParts[iIndex[0]][i] < dMinimum)
    1480                                                                         {
    1481                                                                                 dMinimum = aadOdleglosciParts[iIndex[0]][i];
    1482                                                                         }
    1483                                                                         // przy okazji - sprawdz, czy HUGE_VAL jest rzeczywiscie max dla double
    1484                                                                         assert(aadOdleglosciParts[iIndex[0]][i] < HUGE_VAL);
    1485 
    1486                                                                 }
    1487                                                                 else
    1488                                                                 {
    1489                                                                         if (aadOdleglosciParts[i][iIndex[1]] < dMinimum)
    1490                                                                         {
    1491                                                                                 dMinimum = aadOdleglosciParts[i][iIndex[1]];
    1492                                                                         }
    1493                                                                         // przy okazji - sprawdz, czy HUGE_VAL jest rzeczywiscie max dla double
    1494                                                                         assert(aadOdleglosciParts[i][iIndex[1]] < HUGE_VAL);
    1495                                                                 }
    1496                                                         }
    1497                                                 }
    1498                                                 // sprawdz, czy minimum znaleziono - musi takie byc, bo jest cos niedopasowanego
    1499                                                 assert((dMinimum >= 0.0) && (dMinimum < HUGE_VAL));
    1500 
    1501                                                 // teraz zaznaczamy w tablicy te wszystkie Parts, ktore maja
    1502                                                 // rzeczywiscie te minimalne odleglosci do Part w organizmie iZDrugim
    1503                                                 for (i = 0; i < aiRozmiarCalychGrup[1 - iZDrugim]; i++)
    1504                                                 {
    1505                                                         if (!(m_pMatching->IsMatched(
    1506                                                                 1 - iZDrugim,
    1507                                                                 aiKoniecPierwszejGrupy[1 - iZDrugim] + i)))
    1508                                                         {
    1509                                                                 if (iZDrugim == 0)
    1510                                                                 {
    1511                                                                         // teraz niestety musimy rozpoznac odpowiedni organizm
    1512                                                                         // zeby moc indeksowac niesymetryczna tablice
    1513                                                                         if (aadOdleglosciParts[iIndex[0]][i] == dMinimum)
    1514                                                                         {
    1515                                                                                 // jesli w danym miejscu tablicy odleglosci jest faktyczne
    1516                                                                                 // minimum odleglosci, to mamy potencjalna pare dla iIndex[1]
    1517                                                                                 apvbCzyMinimalnaOdleglosc[0]->operator[](i) = true;
    1518                                                                         }
    1519                                                                 }
    1520                                                                 else
    1521                                                                 {
    1522                                                                         if (aadOdleglosciParts[i][iIndex[1]] == dMinimum)
    1523                                                                         {
    1524                                                                                 apvbCzyMinimalnaOdleglosc[1]->operator[](i) = true;
    1525                                                                         }
    1526                                                                 }
    1527                                                         }
    1528                                                 }
    1529 
    1530                                                 // a teraz - po znalezieniu potencjalnych elementow do dopasowania
    1531                                                 // dopasowujemy pierwszy z potencjalnych
    1532                                                 iDopasowywany = -1;
    1533                                                 for (i = 0; i < aiRozmiarCalychGrup[1 - iZDrugim]; i++)
    1534                                                 {
    1535                                                         if (apvbCzyMinimalnaOdleglosc[iZDrugim]->operator[](i))
     1379                                                        if (apvbCzyMinimalnaOdleglosc[iBezDrugiego]->operator[](i))
    15361380                                                        {
    15371381                                                                iDopasowywany = i;
     
    15391383                                                        }
    15401384                                                }
    1541                                                 // musielismy znalezc jakiegos dopasowywanego
     1385                                                // sprawdz poprawnosc indeksu dopasowywanego (musimy cos znalezc!)
     1386                                                // nieujemny i w odpowiedniej grupie!
    15421387                                                assert((iDopasowywany >= 0) &&
    1543                                                         (iDopasowywany < aiRozmiarCalychGrup[1 - iZDrugim]));
    1544 
    1545                                                 // no to juz mozemy dopasowac
     1388                                                        (iDopasowywany < aiRozmiarCalychGrup[1 - iBezDrugiego]));
     1389
     1390                                                // znalezlismy 1. Part z listy dopasowania - dopasowujemy!
    15461391                                                m_pMatching->Match(
    1547                                                         iZDrugim,
    1548                                                         aiKoniecPierwszejGrupy[iZDrugim] + iIndex[iZDrugim],
    1549                                                         1 - iZDrugim,
    1550                                                         aiKoniecPierwszejGrupy[1 - iZDrugim] + iDopasowywany);
    1551                                                 DB(printf("Przypadek 2.1.: dopasowane Parts dopasowaniebz drugim: (%2i, %2i)\n", aiKoniecPierwszejGrupy[iZDrugim] + iIndex[iZDrugim], aiKoniecPierwszejGrupy[1 - iZDrugim] + iDopasowywany);)
    1552 
    1553                                                         //aiKoniecPierwszejGrupy[ 1-iBezDrugiego ] + iDopasowywany ;)
    1554                                                         //aiKoniecPierwszejGrupy[ 1-iBezDrugiego ] + iDopasowywany ;)
     1392                                                        iBezDrugiego,
     1393                                                        aiKoniecPierwszejGrupy[iBezDrugiego] + iIndex[iBezDrugiego],
     1394                                                        1 - iBezDrugiego,
     1395                                                        aiKoniecPierwszejGrupy[1 - iBezDrugiego] + iDopasowywany);
     1396                                                DB(printf("Przypadek 2.1.: dopasowane Parts dopasowanie bez drugiego: (%2i, %2i)\n", aiKoniecPierwszejGrupy[iBezDrugiego] + iIndex[iBezDrugiego],
     1397                                                        aiKoniecPierwszejGrupy[1 - iBezDrugiego] + iDopasowywany);)
     1398
    15551399                                                        // zmniejsz liczby niedopasowanych elementow w grupach
    15561400                                                        aiRozmiarGrupy[0]--;
    15571401                                                aiRozmiarGrupy[1]--;
    15581402
    1559                                         }
    1560                                         else
    1561                                         {
    1562                                                 // jedna z grup sie juz wyczerpala
    1563                                                 // wiec nie ma mozliwosci dopasowania tego drugiego Partu
    1564                                                 /// i trzeba poczekac na zmiane grupy
    1565                                         }
    1566 
    1567                                         DB(printf("Przypadek 2.\n");)
     1403                                                // sprawdz, czy jest szansa dopasowania tego Part z drugiej strony
     1404                                                // (ktora miala mozliwosc dopasowania tego Part z poprzedniego organizmu)
     1405                                                if ((aiRozmiarGrupy[0] > 0) && (aiRozmiarGrupy[1] > 0))
     1406                                                {
     1407                                                        // jesli grupy sie jeszcze nie wyczrpaly
     1408                                                        // to jest mozliwosc dopasowania w organizmie
     1409
     1410                                                        int iZDrugim = 1 - iBezDrugiego;
     1411                                                        // sprawdz indeks organizmu
     1412                                                        assert((iZDrugim == 0) || (iZDrugim == 1));
     1413
     1414                                                        // bedziemy szukac minimum wsrod niedopasowanych - musimy wykasowac
     1415                                                        // poprzednie obliczenia minimum
     1416                                                        // dla organizmu 1 (o rozmiarze grupy z 0)
     1417                                                        for (i = 0; i < aiRozmiarCalychGrup[0]; i++)
     1418                                                        {
     1419                                                                apvbCzyMinimalnaOdleglosc[1]->operator[](i) = false;
     1420                                                        }
     1421                                                        // dla organizmu 0 (o rozmiarze grup z 1)
     1422                                                        for (i = 0; i < aiRozmiarCalychGrup[1]; i++)
     1423                                                        {
     1424                                                                apvbCzyMinimalnaOdleglosc[0]->operator[](i) = false;
     1425                                                        }
     1426
     1427                                                        // szukamy na nowo minimum dla Part o indeksie iIndex[ iZDrugim ] w organizmie iZDrugim
     1428                                                        // wsrod niedopasowanych Parts z organizmu 1 - iZDrugim
     1429                                                        dMinimum = HUGE_VAL;
     1430                                                        for (i = 0; i < aiRozmiarCalychGrup[1 - iZDrugim]; i++)
     1431                                                        {
     1432                                                                if (!(m_pMatching->IsMatched(
     1433                                                                        1 - iZDrugim,
     1434                                                                        aiKoniecPierwszejGrupy[1 - iZDrugim] + i)))
     1435                                                                {
     1436                                                                        // szukamy minimum obliczonej lokalnej odleglosci tylko wsrod
     1437                                                                        // niedopasowanych jeszcze Parts
     1438                                                                        if (iZDrugim == 0)
     1439                                                                        {
     1440                                                                                // teraz niestety musimy rozpoznac odpowiedni organizm
     1441                                                                                // zeby moc indeksowac niesymetryczna tablice
     1442                                                                                if (aadOdleglosciParts[iIndex[0]][i] < dMinimum)
     1443                                                                                {
     1444                                                                                        dMinimum = aadOdleglosciParts[iIndex[0]][i];
     1445                                                                                }
     1446                                                                                // przy okazji - sprawdz, czy HUGE_VAL jest rzeczywiscie max dla double
     1447                                                                                assert(aadOdleglosciParts[iIndex[0]][i] < HUGE_VAL);
     1448
     1449                                                                        }
     1450                                                                        else
     1451                                                                        {
     1452                                                                                if (aadOdleglosciParts[i][iIndex[1]] < dMinimum)
     1453                                                                                {
     1454                                                                                        dMinimum = aadOdleglosciParts[i][iIndex[1]];
     1455                                                                                }
     1456                                                                                // przy okazji - sprawdz, czy HUGE_VAL jest rzeczywiscie max dla double
     1457                                                                                assert(aadOdleglosciParts[i][iIndex[1]] < HUGE_VAL);
     1458                                                                        }
     1459                                                                }
     1460                                                        }
     1461                                                        // sprawdz, czy minimum znaleziono - musi takie byc, bo jest cos niedopasowanego
     1462                                                        assert((dMinimum >= 0.0) && (dMinimum < HUGE_VAL));
     1463
     1464                                                        // teraz zaznaczamy w tablicy te wszystkie Parts, ktore maja
     1465                                                        // rzeczywiscie te minimalne odleglosci do Part w organizmie iZDrugim
     1466                                                        for (i = 0; i < aiRozmiarCalychGrup[1 - iZDrugim]; i++)
     1467                                                        {
     1468                                                                if (!(m_pMatching->IsMatched(
     1469                                                                        1 - iZDrugim,
     1470                                                                        aiKoniecPierwszejGrupy[1 - iZDrugim] + i)))
     1471                                                                {
     1472                                                                        if (iZDrugim == 0)
     1473                                                                        {
     1474                                                                                // teraz niestety musimy rozpoznac odpowiedni organizm
     1475                                                                                // zeby moc indeksowac niesymetryczna tablice
     1476                                                                                if (aadOdleglosciParts[iIndex[0]][i] == dMinimum)
     1477                                                                                {
     1478                                                                                        // jesli w danym miejscu tablicy odleglosci jest faktyczne
     1479                                                                                        // minimum odleglosci, to mamy potencjalna pare dla iIndex[1]
     1480                                                                                        apvbCzyMinimalnaOdleglosc[0]->operator[](i) = true;
     1481                                                                                }
     1482                                                                        }
     1483                                                                        else
     1484                                                                        {
     1485                                                                                if (aadOdleglosciParts[i][iIndex[1]] == dMinimum)
     1486                                                                                {
     1487                                                                                        apvbCzyMinimalnaOdleglosc[1]->operator[](i) = true;
     1488                                                                                }
     1489                                                                        }
     1490                                                                }
     1491                                                        }
     1492
     1493                                                        // a teraz - po znalezieniu potencjalnych elementow do dopasowania
     1494                                                        // dopasowujemy pierwszy z potencjalnych
     1495                                                        iDopasowywany = -1;
     1496                                                        for (i = 0; i < aiRozmiarCalychGrup[1 - iZDrugim]; i++)
     1497                                                        {
     1498                                                                if (apvbCzyMinimalnaOdleglosc[iZDrugim]->operator[](i))
     1499                                                                {
     1500                                                                        iDopasowywany = i;
     1501                                                                        break;
     1502                                                                }
     1503                                                        }
     1504                                                        // musielismy znalezc jakiegos dopasowywanego
     1505                                                        assert((iDopasowywany >= 0) &&
     1506                                                                (iDopasowywany < aiRozmiarCalychGrup[1 - iZDrugim]));
     1507
     1508                                                        // no to juz mozemy dopasowac
     1509                                                        m_pMatching->Match(
     1510                                                                iZDrugim,
     1511                                                                aiKoniecPierwszejGrupy[iZDrugim] + iIndex[iZDrugim],
     1512                                                                1 - iZDrugim,
     1513                                                                aiKoniecPierwszejGrupy[1 - iZDrugim] + iDopasowywany);
     1514                                                        DB(printf("Przypadek 2.1.: dopasowane Parts dopasowaniebz drugim: (%2i, %2i)\n", aiKoniecPierwszejGrupy[iZDrugim] + iIndex[iZDrugim], aiKoniecPierwszejGrupy[1 - iZDrugim] + iDopasowywany);)
     1515
     1516                                                                //aiKoniecPierwszejGrupy[ 1-iBezDrugiego ] + iDopasowywany ;)
     1517                                                                //aiKoniecPierwszejGrupy[ 1-iBezDrugiego ] + iDopasowywany ;)
     1518                                                                // zmniejsz liczby niedopasowanych elementow w grupach
     1519                                                                aiRozmiarGrupy[0]--;
     1520                                                        aiRozmiarGrupy[1]--;
     1521
     1522                                                }
     1523                                                else
     1524                                                {
     1525                                                        // jedna z grup sie juz wyczerpala
     1526                                                        // wiec nie ma mozliwosci dopasowania tego drugiego Partu
     1527                                                        /// i trzeba poczekac na zmiane grupy
     1528                                                }
     1529
     1530                                                DB(printf("Przypadek 2.\n");)
    15681531
    15691532                                        }// PRZYPADEK 2.
     
    17241687                if (m_adFactors[3] > 0)
    17251688                {
    1726                 // tutaj zacznij pętlę po przekształceniach  geometrycznych
    1727                 const int NO_OF_TRANSFORM = 8; // liczba transformacji geometrycznych (na razie tylko ID i O_YZ)
    1728                 // tablice transformacji współrzędnych; nie są to dokładnie tablice tranformacji, ale raczej tablice PRZEJŚĆ
    1729                 // pomiędzy transformacjami;
    1730                 // wartości orginalne transformacji dOrig uzyskuje się przez:
    1731                 // for ( iTrans = 0; iTrans <= TRANS_INDEX; iTrans++ ) dOrig *= dMul[ iTrans ];
    1732                 //const char *szTransformNames[NO_OF_TRANSFORM] = { "ID", "S_yz", "S_xz", "S_xy", "R180_z", "R180_y", "R180_z", "S_(0,0,0)" };
    1733                 const int dMulX[NO_OF_TRANSFORM] = { 1, -1, -1, 1, -1, 1, -1, -1 };
    1734                 const int dMulY[NO_OF_TRANSFORM] = { 1, 1, -1, -1, -1, -1, -1, 1 };
    1735                 const int dMulZ[NO_OF_TRANSFORM] = { 1, 1, 1, -1, -1, -1, 1, 1 };
     1689                        // tutaj zacznij pętlę po przekształceniach  geometrycznych
     1690                        const int NO_OF_TRANSFORM = 8; // liczba transformacji geometrycznych (na razie tylko ID i O_YZ)
     1691                        // tablice transformacji współrzędnych; nie są to dokładnie tablice tranformacji, ale raczej tablice PRZEJŚĆ
     1692                        // pomiędzy transformacjami;
     1693                        // wartości orginalne transformacji dOrig uzyskuje się przez:
     1694                        // for ( iTrans = 0; iTrans <= TRANS_INDEX; iTrans++ ) dOrig *= dMul[ iTrans ];
     1695                        //const char *szTransformNames[NO_OF_TRANSFORM] = { "ID", "S_yz", "S_xz", "S_xy", "R180_z", "R180_y", "R180_z", "S_(0,0,0)" };
     1696                        const int dMulX[NO_OF_TRANSFORM] = { 1, -1, -1, 1, -1, 1, -1, -1 };
     1697                        const int dMulY[NO_OF_TRANSFORM] = { 1, 1, -1, -1, -1, -1, -1, 1 };
     1698                        const int dMulZ[NO_OF_TRANSFORM] = { 1, 1, 1, -1, -1, -1, 1, 1 };
    17361699
    17371700#ifdef max
    17381701#undef max //this macro would conflict with line below
    17391702#endif
    1740                 double dMinSimValue = std::numeric_limits<double>::max(); // minimum value of similarity
    1741                 int iMinSimTransform = -1; // index of the transformation with the minimum similarity
    1742                 SimilMatching *pMinSimMatching = NULL; // matching with the minimum value of similarity
    1743 
    1744                 // remember the original positions of model 0 as computed by SVD in order to restore them later, after
    1745                 // all transformations have been computed
    1746                 Pt3D *StoredPositions = NULL; // array of positions of Parts, for one (0th) model
    1747                 // create the stored positions
    1748                 StoredPositions = new Pt3D[m_Mod[0]->getPartCount()];
    1749                 assert(StoredPositions != NULL);
    1750                 // copy the original positions of model 0 (store them)
    1751                 int iPart; // a counter of Parts
    1752                 for (iPart = 0; iPart < m_Mod[0]->getPartCount(); iPart++)
    1753                 {
    1754                         StoredPositions[iPart].x = m_aPositions[0][iPart].x;
    1755                         StoredPositions[iPart].y = m_aPositions[0][iPart].y;
    1756                         StoredPositions[iPart].z = m_aPositions[0][iPart].z;
    1757                 }
    1758                 // now the original positions of model 0 are stored
    1759 
    1760 
    1761                 int iTransform; // a counter of geometric transformations
    1762                 for (iTransform = 0; iTransform < NO_OF_TRANSFORM; iTransform++)
    1763                 {
    1764                         // for each geometric transformation to be done
    1765                         // entry conditions:
    1766                         // - models (m_Mod) exist and are available
    1767                         // - matching (m_pMatching) exists and is empty
    1768                         // - all properties are created and available (m_aDegrees and m_aPositions)
    1769 
    1770                         // recompute geometric properties according to the transformation iTransform
    1771                         // but only for model 0
     1703                        double dMinSimValue = std::numeric_limits<double>::max(); // minimum value of similarity
     1704                        int iMinSimTransform = -1; // index of the transformation with the minimum similarity
     1705                        SimilMatching *pMinSimMatching = NULL; // matching with the minimum value of similarity
     1706
     1707                        // remember the original positions of model 0 as computed by SVD in order to restore them later, after
     1708                        // all transformations have been computed
     1709                        Pt3D *StoredPositions = NULL; // array of positions of Parts, for one (0th) model
     1710                        // create the stored positions
     1711                        StoredPositions = new Pt3D[m_Mod[0]->getPartCount()];
     1712                        assert(StoredPositions != NULL);
     1713                        // copy the original positions of model 0 (store them)
     1714                        int iPart; // a counter of Parts
    17721715                        for (iPart = 0; iPart < m_Mod[0]->getPartCount(); iPart++)
    17731716                        {
    1774                                 // for each iPart, a part of the model iMod
    1775                                 m_aPositions[0][iPart].x *= dMulX[iTransform];
    1776                                 m_aPositions[0][iPart].y *= dMulY[iTransform];
    1777                                 m_aPositions[0][iPart].z *= dMulZ[iTransform];
     1717                                StoredPositions[iPart].x = m_aPositions[0][iPart].x;
     1718                                StoredPositions[iPart].y = m_aPositions[0][iPart].y;
     1719                                StoredPositions[iPart].z = m_aPositions[0][iPart].z;
    17781720                        }
    1779                         // now the positions are recomputed
    1780                         ComputeMatching();
    1781 
    1782                         // teraz m_pMatching istnieje i jest pełne
    1783                         assert(m_pMatching != NULL);
    1784                         assert(m_pMatching->IsFull() == true);
    1785 
    1786                         // wykorzystaj to dopasowanie i geometrię do obliczenia tymczasowej wartości miary
    1787                         int iTempRes = CountPartsDistance();
    1788                         // załóż sukces
    1789                         assert(iTempRes == 1);
    1790                         double dCurrentSim = m_adFactors[0] * double(m_iDV) +
    1791                                 m_adFactors[1] * double(m_iDD) +
    1792                                 m_adFactors[2] * double(m_iDN) +
    1793                                 m_adFactors[3] * double(m_dDG);
    1794                         // załóż poprawną wartość podobieństwa
    1795                         assert(dCurrentSim >= 0.0);
    1796 
    1797                         // porównaj wartość obliczoną z dotychczasowym minimum
    1798                         if (dCurrentSim < dMinSimValue)
    1799                         {
    1800                                 // jeśli uzyskano mniejszą wartość dopasowania,
    1801                                 // to zapamiętaj to przekształcenie geometryczne i dopasowanie
    1802                                 dMinSimValue = dCurrentSim;
    1803                                 iMinSimTransform = iTransform;
    1804                                 SAFEDELETE(pMinSimMatching);
    1805                                 pMinSimMatching = new SimilMatching(*m_pMatching);
    1806                                 assert(pMinSimMatching != NULL);
     1721                        // now the original positions of model 0 are stored
     1722
     1723
     1724                        int iTransform; // a counter of geometric transformations
     1725                        for (iTransform = 0; iTransform < NO_OF_TRANSFORM; iTransform++)
     1726                        {
     1727                                // for each geometric transformation to be done
     1728                                // entry conditions:
     1729                                // - models (m_Mod) exist and are available
     1730                                // - matching (m_pMatching) exists and is empty
     1731                                // - all properties are created and available (m_aDegrees and m_aPositions)
     1732
     1733                                // recompute geometric properties according to the transformation iTransform
     1734                                // but only for model 0
     1735                                for (iPart = 0; iPart < m_Mod[0]->getPartCount(); iPart++)
     1736                                {
     1737                                        // for each iPart, a part of the model iMod
     1738                                        m_aPositions[0][iPart].x *= dMulX[iTransform];
     1739                                        m_aPositions[0][iPart].y *= dMulY[iTransform];
     1740                                        m_aPositions[0][iPart].z *= dMulZ[iTransform];
     1741                                }
     1742                                // now the positions are recomputed
     1743                                ComputeMatching();
     1744
     1745                                // teraz m_pMatching istnieje i jest pełne
     1746                                assert(m_pMatching != NULL);
     1747                                assert(m_pMatching->IsFull() == true);
     1748
     1749                                // wykorzystaj to dopasowanie i geometrię do obliczenia tymczasowej wartości miary
     1750                                int iTempRes = CountPartsDistance();
     1751                                // załóż sukces
     1752                                assert(iTempRes == 1);
     1753                                double dCurrentSim = m_adFactors[0] * double(m_iDV) +
     1754                                        m_adFactors[1] * double(m_iDD) +
     1755                                        m_adFactors[2] * double(m_iDN) +
     1756                                        m_adFactors[3] * double(m_dDG);
     1757                                // załóż poprawną wartość podobieństwa
     1758                                assert(dCurrentSim >= 0.0);
     1759
     1760                                // porównaj wartość obliczoną z dotychczasowym minimum
     1761                                if (dCurrentSim < dMinSimValue)
     1762                                {
     1763                                        // jeśli uzyskano mniejszą wartość dopasowania,
     1764                                        // to zapamiętaj to przekształcenie geometryczne i dopasowanie
     1765                                        dMinSimValue = dCurrentSim;
     1766                                        iMinSimTransform = iTransform;
     1767                                        SAFEDELETE(pMinSimMatching);
     1768                                        pMinSimMatching = new SimilMatching(*m_pMatching);
     1769                                        assert(pMinSimMatching != NULL);
     1770                                }
     1771
     1772                                // teraz już można usunąć stare dopasowanie (dla potrzeb następnego przebiegu pętli)
     1773                                m_pMatching->Empty();
     1774                        } // for ( iTransform )
     1775
     1776                        // po pętli przywróć najlepsze dopasowanie
     1777                        delete m_pMatching;
     1778                        m_pMatching = pMinSimMatching;
     1779
     1780                        DB(printf("Matching has been chosen!\n");)
     1781                                // print the name of the chosen transformation:
     1782                                // printf("Chosen transformation: %s\n", szTransformNames[ iMinSimTransform ] );
     1783
     1784                                // i przywróć jednocześnie pozycje Parts modelu 0 dla tego dopasowania
     1785                                // - najpierw przywroc Parts pozycje orginalne, po SVD
     1786                                for (iPart = 0; iPart < m_Mod[0]->getPartCount(); iPart++)
     1787                                {
     1788                                        m_aPositions[0][iPart].x = StoredPositions[iPart].x;
     1789                                        m_aPositions[0][iPart].y = StoredPositions[iPart].y;
     1790                                        m_aPositions[0][iPart].z = StoredPositions[iPart].z;
     1791                                }
     1792                        // - usun teraz zapamietane pozycje Parts
     1793                        delete[] StoredPositions;
     1794                        // - a teraz oblicz na nowo wspolrzedne wlasciwego przeksztalcenia dla model 0
     1795                        for (iTransform = 0; iTransform <= iMinSimTransform; iTransform++)
     1796                        {
     1797                                // for each transformation before and INCLUDING iMinTransform
     1798                                // do the transformation (only model 0)
     1799                                for (iPart = 0; iPart < m_Mod[0]->getPartCount(); iPart++)
     1800                                {
     1801                                        m_aPositions[0][iPart].x *= dMulX[iTransform];
     1802                                        m_aPositions[0][iPart].y *= dMulY[iTransform];
     1803                                        m_aPositions[0][iPart].z *= dMulZ[iTransform];
     1804                                }
    18071805                        }
    1808 
    1809                         // teraz już można usunąć stare dopasowanie (dla potrzeb następnego przebiegu pętli)
    1810                         m_pMatching->Empty();
    1811                 } // for ( iTransform )
    1812 
    1813                 // po pętli przywróć najlepsze dopasowanie
    1814                 delete m_pMatching;
    1815                 m_pMatching = pMinSimMatching;
    1816 
    1817                 DB(printf("Matching has been chosen!\n");)
    1818                         // print the name of the chosen transformation:
    1819                         // printf("Chosen transformation: %s\n", szTransformNames[ iMinSimTransform ] );
    1820 
    1821                         // i przywróć jednocześnie pozycje Parts modelu 0 dla tego dopasowania
    1822                         // - najpierw przywroc Parts pozycje orginalne, po SVD
    1823                         for (iPart = 0; iPart < m_Mod[0]->getPartCount(); iPart++)
    1824                         {
    1825                         m_aPositions[0][iPart].x = StoredPositions[iPart].x;
    1826                         m_aPositions[0][iPart].y = StoredPositions[iPart].y;
    1827                         m_aPositions[0][iPart].z = StoredPositions[iPart].z;
    1828                         }
    1829                 // - usun teraz zapamietane pozycje Parts
    1830                 delete[] StoredPositions;
    1831                 // - a teraz oblicz na nowo wspolrzedne wlasciwego przeksztalcenia dla model 0
    1832                 for (iTransform = 0; iTransform <= iMinSimTransform; iTransform++)
    1833                 {
    1834                         // for each transformation before and INCLUDING iMinTransform
    1835                         // do the transformation (only model 0)
    1836                         for (iPart = 0; iPart < m_Mod[0]->getPartCount(); iPart++)
    1837                         {
    1838                                 m_aPositions[0][iPart].x *= dMulX[iTransform];
    1839                                 m_aPositions[0][iPart].y *= dMulY[iTransform];
    1840                                 m_aPositions[0][iPart].z *= dMulZ[iTransform];
    1841                         }
    1842                 }
    18431806
    18441807                }
     
    18551818
    18561819        DB(_PrintPartsMatching();)
    1857 
    18581820
    18591821                return 1;
     
    19851947                int nSize = m_Mod[iMod]->getPartCount();
    19861948
    1987                 double *pDistances =  new double[nSize * nSize];
     1949                double *pDistances = new double[nSize * nSize];
    19881950
    19891951                for (int i = 0; i < nSize; i++)
     
    19971959                Pt3D P1Pos, P2Pos; // positions of parts
    19981960                double dDistance; // the distance between Parts
    1999                
     1961
    20001962                double *weights = new double[nSize];
    20011963                for (int i = 0; i < nSize; i++)
    20021964                {
    2003                         if (wMDS==1)
     1965                        if (wMDS == 1)
    20041966                                weights[i] = 0;
    20051967                        else
    20061968                                weights[i] = 1;
    20071969                }
    2008                                
    2009                 if (wMDS==1)
     1970
     1971                if (wMDS == 1)
    20101972                        for (int i = 0; i < pModel->getJointCount(); i++)
    20111973                        {
    20121974                                weights[pModel->getJoint(i)->p1_refno]++;
    2013                                 weights[pModel->getJoint(i)->p2_refno]++; 
     1975                                weights[pModel->getJoint(i)->p2_refno]++;
    20141976                        }
    2015                
     1977
    20161978                for (iP1 = 0; iP1 < pModel->getPartCount(); iP1++)
    20171979                {
     
    20572019}
    20582020
     2021/**     Evaluates distance between two given genotypes. The distance depends strongly
     2022        on weights set and the matching algorithm used.
     2023        @param G0 Pointer to the first of compared genotypes
     2024        @param G1 Pointer to the second of compared genotypes.
     2025        @return Distance between two genotypes.
     2026        @sa m_adFactors, matching_method
     2027        */
     2028double ModelSimil::EvaluateDistance(const Geno *G0, const Geno *G1)
     2029{
     2030        return matching_method == 0 ? EvaluateDistanceHungarian(G0, G1) : EvaluateDistanceGreedy(G0, G1);
     2031}
     2032
    20592033void ModelSimil::p_evaldistance(ExtValue *args, ExtValue *ret)
    20602034{
     
    20662040                ret->setDouble(EvaluateDistance(g1, g2));
    20672041}
     2042
     2043void ModelSimil::FillPartsDistances(double*& dist, int bigger, int smaller, bool geo)
     2044{
     2045        for (int i = 0; i < bigger; i++)
     2046        {
     2047                for (int j = 0; j < bigger; j++)
     2048                {
     2049                        // assign penalty for unassignment for vertex from bigger model
     2050                        if (j >= smaller)
     2051                        {
     2052                                if (geo)
     2053                                        dist[i*bigger + j] += m_adFactors[3] * m_aPositions[1 - m_iSmaller][i].length();
     2054                                else
     2055                                        dist[i*bigger + j] = m_adFactors[1] * m_aDegrees[1 - m_iSmaller][i][DEGREE] + m_adFactors[2] * m_aDegrees[1 - m_iSmaller][i][NEURONS];
     2056                        }
     2057                        // compute distance between parts
     2058                        else
     2059                        {
     2060                                if (geo)
     2061                                        dist[i*bigger + j] += m_adFactors[3] * m_aPositions[1 - m_iSmaller][i].distanceTo(m_aPositions[m_iSmaller][j]);
     2062                                else
     2063                                        dist[i*bigger + j] = m_adFactors[1] * abs(m_aDegrees[1 - m_iSmaller][i][DEGREE] - m_aDegrees[m_iSmaller][j][DEGREE])
     2064                                        + m_adFactors[2] * abs(m_aDegrees[1 - m_iSmaller][i][NEURONS] - m_aDegrees[m_iSmaller][j][NEURONS]);
     2065                        }
     2066
     2067                }
     2068        }
     2069}
     2070
     2071double ModelSimil::EvaluateDistanceHungarian(const Geno *G0, const Geno *G1)
     2072{
     2073        double dResult = 0;
     2074
     2075        m_Gen[0] = G0;
     2076        m_Gen[1] = G1;
     2077
     2078        // check whether pointers are not NULL
     2079        if (m_Gen[0] == NULL || m_Gen[1] == NULL)
     2080        {
     2081                DB(printf("ModelSimil::EvaluateDistanceHungarian - invalid genotypes pointers\n");)
     2082                        return 0.0;
     2083        }
     2084        // create models of objects to compare
     2085        m_Mod[0] = new Model(*(m_Gen[0]));
     2086        m_Mod[1] = new Model(*(m_Gen[1]));
     2087
     2088        // validate models
     2089        if (m_Mod[0] == NULL || m_Mod[1] == NULL || !(m_Mod[0]->isValid()) || !(m_Mod[1]->isValid()))
     2090        {
     2091                DB(printf("ModelSimil::EvaluateDistanceHungarian - invalid models pointers\n");)
     2092                        return 0.0;
     2093        }
     2094
     2095        //Get information about vertex degrees, neurons and 3D location
     2096        if (!CreatePartInfoTables())
     2097                return 0;
     2098        if (!CountPartDegrees())
     2099                return 0;
     2100        if (!GetPartPositions())
     2101                return 0;
     2102        if (!CountPartNeurons())
     2103                return 0;
     2104
     2105        m_iSmaller = m_Mod[0]->getPartCount() <= m_Mod[1]->getPartCount() ? 0 : 1;
     2106        int nSmaller = m_Mod[m_iSmaller]->getPartCount();
     2107        int nBigger = m_Mod[1 - m_iSmaller]->getPartCount();
     2108
     2109        double* partsDistances = new double[nBigger*nBigger]();
     2110        FillPartsDistances(partsDistances, nBigger, nSmaller, false);
     2111        int *assignment = new int[nBigger]();
     2112
     2113        HungarianAlgorithm hungarian;
     2114
     2115        if (m_adFactors[3] > 0)
     2116        {
     2117                if (!ComputePartsPositionsBySVD())
     2118                {
     2119                        return 0;
     2120                }
     2121
     2122                // tutaj zacznij pętlę po przekształceniach  geometrycznych
     2123                const int NO_OF_TRANSFORM = 8; // liczba transformacji geometrycznych (na razie tylko ID i O_YZ)
     2124                // tablice transformacji współrzędnych; nie są to dokładnie tablice tranformacji, ale raczej tablice PRZEJŚĆ
     2125                // pomiędzy transformacjami;
     2126                const int dMulX[NO_OF_TRANSFORM] = { 1, -1, -1, 1, -1, 1, -1, -1 };
     2127                const int dMulY[NO_OF_TRANSFORM] = { 1, 1, -1, -1, -1, -1, -1, 1 };
     2128                const int dMulZ[NO_OF_TRANSFORM] = { 1, 1, 1, -1, -1, -1, 1, 1 };
     2129
     2130                std::vector<int> minAssignment(nBigger);
     2131#ifdef max
     2132#undef max //this macro would conflict with line below
     2133#endif
     2134                double dMinSimValue = std::numeric_limits<double>::max(); // minimum value of similarity
     2135
     2136                int iTransform; // a counter of geometric transformations
     2137                for (iTransform = 0; iTransform < NO_OF_TRANSFORM; iTransform++)
     2138                {
     2139                        // for each geometric transformation to be done
     2140                        // entry conditions:
     2141                        // - models (m_Mod) exist and are available
     2142                        // - all properties are created and available (m_aDegrees and m_aPositions)
     2143                        double* tmpPartsDistances = new double[nBigger*nBigger]();
     2144                        std::copy(partsDistances, partsDistances + nBigger * nBigger, tmpPartsDistances);
     2145                        // recompute geometric properties according to the transformation iTransform
     2146                        // but only for model 0
     2147                        for (int iPart = 0; iPart < m_Mod[m_iSmaller]->getPartCount(); iPart++)
     2148                        {
     2149                                // for each iPart, a part of the model iMod
     2150                                m_aPositions[m_iSmaller][iPart].x *= dMulX[iTransform];
     2151                                m_aPositions[m_iSmaller][iPart].y *= dMulY[iTransform];
     2152                                m_aPositions[m_iSmaller][iPart].z *= dMulZ[iTransform];
     2153                        }
     2154                        // now the positions are recomputed
     2155
     2156                        FillPartsDistances(tmpPartsDistances, nBigger, nSmaller, true);
     2157                        std::fill_n(assignment, nBigger, 0);
     2158                        double dCurrentSim = hungarian.Solve(tmpPartsDistances, assignment, nBigger, nBigger);
     2159
     2160                        delete[] tmpPartsDistances;
     2161                        // załóż poprawną wartość podobieństwa
     2162                        assert(dCurrentSim >= 0.0);
     2163
     2164                        // porównaj wartość obliczoną z dotychczasowym minimum
     2165                        if (dCurrentSim < dMinSimValue)
     2166                        {
     2167                                dMinSimValue = dCurrentSim;
     2168                                if (saveMatching == 1)
     2169                                {
     2170                                        minAssignment.clear();
     2171                                        minAssignment.insert(minAssignment.begin(), assignment, assignment + nBigger);
     2172                                }
     2173                        }
     2174                }
     2175
     2176                dResult = dMinSimValue;
     2177                if (saveMatching == 1)
     2178                        std::copy(minAssignment.begin(), minAssignment.end(), assignment);
     2179        }
     2180
     2181        else
     2182        {
     2183                dResult = hungarian.Solve(partsDistances, assignment, nBigger, nBigger);
     2184        }
     2185
     2186        //add difference in anywhere and onJoint neurons
     2187        dResult += m_adFactors[2] * (abs(m_aOnJoint[0][3] - m_aOnJoint[1][3]) + abs(m_aAnywhere[0][3] - m_aAnywhere[1][3]));
     2188        //add difference in part numbers
     2189        dResult += (nBigger - nSmaller) * m_adFactors[0];
     2190
     2191        // delete degree arrays created in CreatePartInfoTables
     2192        SAFEDELETEARRAY(m_aDegrees[0]);
     2193        SAFEDELETEARRAY(m_aDegrees[1]);
     2194
     2195        // and position arrays
     2196        SAFEDELETEARRAY(m_aPositions[0]);
     2197        SAFEDELETEARRAY(m_aPositions[1]);
     2198
     2199        // delete created models
     2200        SAFEDELETE(m_Mod[0]);
     2201        SAFEDELETE(m_Mod[1]);
     2202
     2203        delete[] assignment;
     2204        delete[] partsDistances;
     2205
     2206        return dResult;
     2207}
  • cpp/frams/model/similarity/simil_model.h

    r817 r869  
    11// This file is a part of Framsticks SDK.  http://www.framsticks.com/
    2 // Copyright (C) 1999-2016  Maciej Komosinski and Szymon Ulatowski.
     2// Copyright (C) 1999-2019  Maciej Komosinski and Szymon Ulatowski.
    33// See LICENSE.txt for details.
    44
     
    99#include "frams/genetics/geno.h"
    1010#include "frams/model/model.h"
    11 #include "frams/util/3d.h"
    1211#include "simil_match.h"
    1312
     
    2726 * Marek Kubiak (concept, implementation)
    2827 * Maciej Komosinski (concept, Framsticks interface)
    29  * Agnieszka Mensfelt (refactoring)
     28 * Agnieszka Mensfelt (refactoring, improvements)
    3029 */
    3130class ModelSimil
     
    3433        ModelSimil();
    3534        virtual ~ModelSimil();
    36         double EvaluateDistance(const Geno *G0, const Geno *G1);
     35        double EvaluateDistance(const Geno *G0, const Geno *G1); //chooses greedy or hungarian
     36        double EvaluateDistanceGreedy(const Geno *G0, const Geno *G1);
     37        double EvaluateDistanceHungarian(const Geno *G0, const Geno *G1);
    3738
    3839        static int CompareDegrees(const void *pElem1, const void *pElem2);
     40        static int CompareFuzzyDegrees(const void *pElem1, const void *pElem2);
    3941        static int CompareConnsNo(const void *pElem1, const void *pElem2);
    4042        static int GetNOFactors();
     
    4850        int MatchPartsGeometry();
    4951        void ComputeMatching();
     52        void FillPartsDistances(double *&dist, int bigger, int smaller, bool geo);
    5053        void _PrintPartsMatching();
    5154        void SaveIntermediateFiles();
    5255
    53         static int CheckPartsIdentity(Part *P1, Part *P2);
     56        static int CheckPartsIdentity(Part *P1, Part *P2); //TODO exists?
    5457        int SortPartInfoTables();
    5558        int CountPartNeurons();
     
    5861        int CountPartDegrees();
    5962
    60         int SortPartInfoFuzzy();
    6163        void CountFuzzyNeighb();
    6264        void SortFuzzyNeighb();
    6365        void GetNeighbIndexes(int mod, int partInd, std::vector<int> &indexes);
    64         void CheckFuzzyIdentity();
     66        void FuzzyOrder();
    6567
    6668        int CreatePartInfoTables();
     
    6870        void _PrintArray(int *array, int base, int size);
    6971        void _PrintNeighbourhood(int i);
     72        void _PrintFuzzyNeighbourhood(int i);
    7073        void _PrintArrayDouble(double *array, int base, int size);
    7174        int CountPartsDistance();
     
    7376
    7477public:
     78        /// 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
     79        /// @sa EvaluateDistance
     80        int matching_method;
     81
    7582        /// Table of weights for weighted distance function.
    7683        /// Weights are for factors in the following order:
     
    8895        int fuzzyDepth;
    8996        int isFuzzy;
    90        
     97
    9198        //For wMDS = 1 weighted MDS with vertex degrees as weights is used for the alignment.
    9299        int wMDS;
     100
     101        //For saveMatching = 1 the best matching found will be saved.
     102        int saveMatching;
    93103
    94104        /// Interface to local parameters
     
    181191};
    182192
    183 
    184193#endif
Note: See TracChangeset for help on using the changeset viewer.