source: cpp/frams/model/similarity/simil_model.cpp @ 877

Last change on this file since 877 was 877, checked in by Maciej Komosinski, 16 months ago

Introduced a shared function to avoid code duplication; more informative error message; frees memory in case of error

  • Property svn:eol-style set to native
File size: 68.7 KB
Line 
1// This file is a part of Framsticks SDK.  http://www.framsticks.com/
2// Copyright (C) 1999-2019  Maciej Komosinski and Szymon Ulatowski.
3// See LICENSE.txt for details.
4
5// implementation of the ModelSimil class.
6
7#include "SVD/matrix_tools.h"
8#include "hungarian/hungarian.h"
9#include "simil_model.h"
10#include "simil_match.h"
11#include "frams/model/modelparts.h"
12#include "frams/util/list.h"
13#include "common/nonstd.h"
14#include <frams/vm/classes/genoobj.h>
15#ifdef EMSCRIPTEN
16#include <cstdlib>
17#else
18#include <stdlib.h>
19#endif
20#include <math.h>
21#include <string>
22#include <limits>
23#include <assert.h>
24#include <vector>
25#include <algorithm>
26#include <cstdlib> //required for std::qsort in macos xcode
27
28#define DB(x)  //define as x if you want to print debug information
29
30const int ModelSimil::iNOFactors = 4;
31//depth of the fuzzy neighbourhood
32int fuzDepth = 0;
33
34#define FIELDSTRUCT ModelSimil
35
36static ParamEntry MSparam_tab[] = {
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), "",},
39                { "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.", },
40                { "simil_partdeg", 0, 0, "Weight of parts' degree", "f 0 100 1", FIELD(m_adFactors[1]), "", },
41                { "simil_neuro", 0, 0, "Weight of neurons count", "f 0 100 0.1", FIELD(m_adFactors[2]), "", },
42                { "simil_partgeom", 0, 0, "Weight of parts' geometric distances", "f 0 100 0", FIELD(m_adFactors[3]), "", },
43                { "simil_fixedZaxis", 0, 0, "Fix 'z' (vertical) axis?", "d 0 1 0", FIELD(fixedZaxis), "", },
44                { "simil_weightedMDS", 0, 0, "Should weighted MDS be used?", "d 0 1 0", FIELD(wMDS), "If activated, weighted MDS with vertex (i.e., Part) degrees as weights is used for 3D alignment of body structure.", },
45                { "evaluateDistance", 0, PARAM_DONTSAVE | PARAM_USERHIDDEN, "evaluate model dissimilarity", "p f(oGeno,oGeno)", PROCEDURE(p_evaldistance), "Calculates dissimilarity between two models created from Geno objects.", },
46                { 0, },
47};
48
49#undef FIELDSTRUCT
50
51//////////////////////////////////////////////////////////////////////
52// Construction/Destruction
53//////////////////////////////////////////////////////////////////////
54
55/** Constructor. Sets default weights. Initializes other fields with zeros.
56 */
57ModelSimil::ModelSimil() : localpar(MSparam_tab, this), m_iDV(0), m_iDD(0), m_iDN(0), m_dDG(0.0)
58{
59        localpar.setDefault();
60
61        m_Gen[0] = NULL;
62        m_Gen[1] = NULL;
63        m_Mod[0] = NULL;
64        m_Mod[1] = NULL;
65        m_aDegrees[0] = NULL;
66        m_aDegrees[1] = NULL;
67        m_aPositions[0] = NULL;
68        m_aPositions[1] = NULL;
69        m_fuzzyNeighb[0] = NULL;
70        m_fuzzyNeighb[1] = NULL;
71        m_Neighbours[0] = NULL;
72        m_Neighbours[1] = NULL;
73        m_pMatching = NULL;
74
75        //Determines whether "fuzzy vertex degree" should be used.
76        //Currently "fuzzy vertex degree" is inactive.
77        isFuzzy = false;
78        fuzzyDepth = 10;
79
80        //Determines whether weighted MDS should be used.
81        wMDS = 0;
82        //Determines whether best matching should be saved using hungarian similarity measure.
83        saveMatching = false;
84}
85
86double ModelSimil::EvaluateDistanceGreedy(const Geno *G0, const Geno *G1)
87{
88        double dResult = 0.0;
89
90        m_Gen[0] = G0;
91        m_Gen[1] = G1;
92
93        // create models of objects to compare
94        m_Mod[0] = newModel(m_Gen[0]);
95        m_Mod[1] = newModel(m_Gen[1]);
96
97        if (m_Mod[0] == NULL || m_Mod[1] == NULL)
98                return 0.0;
99
100        // difference in the number of vertices (Parts) - positive
101        // find object that has less parts (m_iSmaller)
102        m_iDV = (m_Mod[0]->getPartCount() - m_Mod[1]->getPartCount());
103        if (m_iDV > 0)
104                m_iSmaller = 1;
105        else
106        {
107                m_iSmaller = 0;
108                m_iDV = -m_iDV;
109        }
110
111        // check if index of the smaller organism is a valid index
112        assert((m_iSmaller == 0) || (m_iSmaller == 1));
113        // validate difference in the parts number
114        assert(m_iDV >= 0);
115
116        // create Parts matching object
117        m_pMatching = new SimilMatching(m_Mod[0]->getPartCount(), m_Mod[1]->getPartCount());
118        // validate matching object
119        assert(m_pMatching != NULL);
120        assert(m_pMatching->IsEmpty() == true);
121
122
123        // assign matching function
124        int (ModelSimil::* pfMatchingFunction) () = NULL;
125
126        pfMatchingFunction = &ModelSimil::MatchPartsGeometry;
127
128        // match Parts (vertices of creatures)
129        if ((this->*pfMatchingFunction)() == 0)
130        {
131                logPrintf("ModelSimil", "EvaluateDistanceGreedy", LOG_ERROR, "The matching function returned 0");
132                return 0.0;
133        }
134
135        // after matching function call we must have full matching
136        assert(m_pMatching->IsFull() == true);
137
138        DB(SaveIntermediateFiles();)
139
140                // count differences in matched parts
141                if (CountPartsDistance() == 0)
142                {
143                        logPrintf("ModelSimil", "EvaluateDistanceGreedy", LOG_ERROR, "CountPartDistance()==0");
144                        return 0.0;
145                }
146
147        // delete degree arrays created in CreatePartInfoTables
148        SAFEDELETEARRAY(m_aDegrees[0]);
149        SAFEDELETEARRAY(m_aDegrees[1]);
150
151        // and position arrays
152        SAFEDELETEARRAY(m_aPositions[0]);
153        SAFEDELETEARRAY(m_aPositions[1]);
154
155        // in fuzzy mode delete arrays of neighbourhood and fuzzy neighbourhood
156        if (isFuzzy)
157        {
158                for (int i = 0; i != 2; ++i)
159                {
160                        for (int j = 0; j != m_Mod[i]->getPartCount(); ++j)
161                        {
162                                delete[] m_Neighbours[i][j];
163                                delete[] m_fuzzyNeighb[i][j];
164                        }
165                        delete[] m_Neighbours[i];
166                        delete[] m_fuzzyNeighb[i];
167                }
168
169        }
170
171        // delete created models
172        SAFEDELETE(m_Mod[0]);
173        SAFEDELETE(m_Mod[1]);
174
175        // delete created matching
176        SAFEDELETE(m_pMatching);
177
178        dResult = m_adFactors[0] * double(m_iDV) +
179                m_adFactors[1] * double(m_iDD) +
180                m_adFactors[2] * double(m_iDN) +
181                m_adFactors[3] * double(m_dDG);
182
183        return dResult;
184}
185
186ModelSimil::~ModelSimil()
187{
188        // matching should have been deleted earlier
189        assert(m_pMatching == NULL);
190}
191
192/**     Creates files matching.txt, org0.txt and org1.txt containing information
193 * about the matching and both organisms for visualization purpose.
194 */
195void ModelSimil::SaveIntermediateFiles()
196{
197        assert(m_pMatching->IsFull() == true);
198        printf("Saving the matching to file 'matching.txt'\n");
199        FILE *pMatchingFile = NULL;
200        // try to open the file
201        pMatchingFile = fopen("matching.txt", "wt");
202        assert(pMatchingFile != NULL);
203
204        int iOrgPart; // original index of a Part
205        int nBigger; // index of the larger organism
206
207        // check which object is bigger
208        if (m_pMatching->GetObjectSize(0) >= m_pMatching->GetObjectSize(1))
209        {
210                nBigger = 0;
211        }
212        else
213        {
214                nBigger = 1;
215        }
216
217        // print first line - original indices of Parts of the bigger organism
218        fprintf(pMatchingFile, "[ ");
219        for (iOrgPart = 0; iOrgPart < m_pMatching->GetObjectSize(nBigger); iOrgPart++)
220        {
221                fprintf(pMatchingFile, "%2i ", iOrgPart);
222        }
223        fprintf(pMatchingFile, "] : ORG[%i]\n", nBigger);
224
225        // print second line - matched original indices of the second organism
226        fprintf(pMatchingFile, "[ ");
227        for (iOrgPart = 0; iOrgPart < m_pMatching->GetObjectSize(nBigger); iOrgPart++)
228        {
229                int iSorted; // index of the iOrgPart after sorting (as used by matching)
230                int iSortedMatched; // index of the matched Part (after sorting)
231                int iOrginalMatched; // index of the matched Part (the original one)
232
233                // find the index of iOrgPart after sorting (in m_aDegrees)
234                for (iSorted = 0; iSorted < m_Mod[nBigger]->getPartCount(); iSorted++)
235                {
236                        // for each iSorted, an index in the sorted m_aDegrees array
237                        if (m_aDegrees[nBigger][iSorted][0] == iOrgPart)
238                        {
239                                // if the iSorted Part is the one with iOrgPart as the orginal index
240                                // remember the index
241                                break;
242                        }
243                }
244                // if the index iSorted was found, then this condition is met
245                assert(iSorted < m_Mod[nBigger]->getPartCount());
246
247                // find the matched sorted index
248                if (m_pMatching->IsMatched(nBigger, iSorted))
249                {
250                        // if Part iOrgPart is matched
251                        // then get the matched Part (sorted) index
252                        iSortedMatched = m_pMatching->GetMatchedIndex(nBigger, iSorted);
253                        assert(iSortedMatched >= 0);
254                        // and find its original index
255                        iOrginalMatched = m_aDegrees[1 - nBigger][iSortedMatched][0];
256                        fprintf(pMatchingFile, "%2i ", iOrginalMatched);
257                }
258                else
259                {
260                        // if the Part iOrgPart is not matched
261                        // just print "X"
262                        fprintf(pMatchingFile, " X ");
263                }
264        } // for ( iOrgPart )
265
266        // now all matched Part indices are printed out, end the line
267        fprintf(pMatchingFile, "] : ORG[%i]\n", 1 - nBigger);
268
269        // close the file
270        fclose(pMatchingFile);
271        // END TEMP
272
273        // TEMP
274        // print out the 2D positions of Parts of both of the organisms
275        // to files "org0.txt" and "org1.txt" using the original indices of Parts
276        int iModel; // index of a model (an organism)
277        FILE *pModelFile;
278        for (iModel = 0; iModel < 2; iModel++)
279        {
280                // for each iModel, a model of a compared organism
281                // write its (only 2D) positions to a file "org<iModel>.txt"
282                // construct the model filename "org<iModel>.txt"
283                std::string sModelFilename("org");
284                //              char *szModelIndex = "0"; // the index of the model (iModel) in the character form
285                char szModelIndex[2];
286                sprintf(szModelIndex, "%i", iModel);
287                sModelFilename += szModelIndex;
288                sModelFilename += ".txt";
289                // open the file for writing
290                pModelFile = fopen(sModelFilename.c_str(), "wt"); //FOPEN_WRITE
291                assert(pModelFile != NULL);
292                // write the 2D positions of iModel to the file
293                int iOriginalPart; // an original index of a Part
294                for (iOriginalPart = 0; iOriginalPart < m_Mod[iModel]->getPartCount(); iOriginalPart++)
295                {
296                        // for each iOriginalPart, a Part of the organism iModel
297                        // get the 2D coordinates of the Part
298                        double dPartX = m_aPositions[iModel][iOriginalPart].x;
299                        double dPartY = m_aPositions[iModel][iOriginalPart].y;
300                        // print the line: <iOriginalPart> <dPartX> <dPartY>
301                        fprintf(pModelFile, "%i %.4lf %.4lf\n", iOriginalPart, dPartX, dPartY);
302                }
303                // close the file
304                fclose(pModelFile);
305        }
306}
307
308/** Comparison function required for qsort() call. Used while sorting groups of
309        Parts with respect to degree. Compares two TDN structures
310        with respect to their [1] field (degree). Highest degree goes first.
311        @param pElem1 Pointer to the TDN structure of some Part.
312        @param pElem2 Pointer to the TDN structure of some Part.
313        @return (-1) - pElem1 should go first, 0 - equal, (1) - pElem2 should go first.
314        */
315int ModelSimil::CompareDegrees(const void *pElem1, const void *pElem2)
316{
317        int *tdn1 = (int *)pElem1;
318        int *tdn2 = (int *)pElem2;
319
320        if (tdn1[1] > tdn2[1])
321        {
322                // when degree - tdn1[1] - is BIGGER
323                return -1;
324        }
325        else
326                if (tdn1[1] < tdn2[1])
327                {
328                        // when degree - tdn2[1] - is BIGGER
329                        return 1;
330                }
331                else
332                {
333                        return 0;
334                }
335}
336
337/** Comparison function required for qsort() call. Used while sorting groups of
338        Parts with respect to fuzzy vertex degree. Compares two TDN structures
339        with respect to their [4] field ( fuzzy vertex degree). Highest degree goes first.
340        @param pElem1 Pointer to the TDN structure of some Part.
341        @param pElem2 Pointer to the TDN structure of some Part.
342        @return (-1) - pElem1 should go first, 0 - equal, (1) - pElem2 should go first.
343        */
344int ModelSimil::CompareFuzzyDegrees(const void *pElem1, const void *pElem2)
345{
346        int *tdn1 = (int *)pElem1;
347        int *tdn2 = (int *)pElem2;
348
349        if (tdn1[4] > tdn2[4])
350        {
351                // when degree - tdn1[4] - is BIGGER
352                return -1;
353        }
354        else
355                if (tdn1[4] < tdn2[4])
356                {
357                        // when degree - tdn2[4] - is BIGGER
358                        return 1;
359                }
360                else
361                {
362                        return 0;
363                }
364}
365
366/** Comparison function required for qsort() call. Used while sorting groups of Parts with
367                the same degree. Firstly, compare NIt. Secondly, compare Neu. If both are equal -
368                compare Parts' original index (they are never equal). So this sorting assures
369                that the order obtained is deterministic.
370                @param pElem1 Pointer to the TDN structure of some Part.
371                @param pElem2 Pointer to the TDN structure of some Part.
372                @return (-1) - pElem1 should go first, 0 - equal, (1) - pElem2 should go first.
373                */
374int ModelSimil::CompareConnsNo(const void *pElem1, const void *pElem2)
375{
376        // pointers to TDN arrays
377        int *tdn1, *tdn2;
378        // definitions of elements being compared
379        tdn1 = (int *)pElem1;
380        tdn2 = (int *)pElem2;
381
382        // comparison according to Neural Connections (to jest TDN[2])
383        if (tdn1[NEURO_CONNS] > tdn2[NEURO_CONNS])
384        {
385                // when number of NConn Elem1 is BIGGER
386                return -1;
387        }
388        else
389                if (tdn1[NEURO_CONNS] < tdn2[NEURO_CONNS])
390                {
391                        // when number of NConn Elem1 is SMALLER
392                        return 1;
393                }
394                else
395                {
396                        // when numbers of NConn are EQUAL
397                        // compare Neu numbers (TDN[3])
398                        if (tdn1[NEURONS] > tdn2[NEURONS])
399                        {
400                                // when number of Neu is BIGGER for Elem1
401                                return -1;
402                        }
403                        else
404                                if (tdn1[NEURONS] < tdn2[NEURONS])
405                                {
406                                        // when number of Neu is SMALLER for Elem1
407                                        return 1;
408                                }
409                                else
410                                {
411                                        // when numbers of Nconn and Neu are equal we check original indices
412                                        // of Parts being compared
413
414                                        // comparison according to OrgIndex
415                                        if (tdn1[ORIG_IND] > tdn2[ORIG_IND])
416                                        {
417                                                // when the number of NIt Deg1 id BIGGER
418                                                return -1;
419                                        }
420                                        else
421                                                if (tdn1[ORIG_IND] < tdn2[ORIG_IND])
422                                                {
423                                                        // when the number of NIt Deg1 id SMALLER
424                                                        return 1;
425                                                }
426                                                else
427                                                {
428                                                        // impossible, indices are alway different
429                                                        return 0;
430                                                }
431                                }
432                }
433}
434
435/** Returns number of factors involved in final distance computation.
436                These factors include differences in numbers of parts, degrees,
437                number of neurons.
438                */
439int ModelSimil::GetNOFactors()
440{
441        return ModelSimil::iNOFactors;
442}
443
444/** Prints the array of degrees for the given organism. Debug method.
445 */
446void ModelSimil::_PrintDegrees(int i)
447{
448        int j;
449        printf("Organizm %i :", i);
450        printf("\n      ");
451        for (j = 0; j < m_Mod[i]->getPartCount(); j++)
452                printf("%3i ", j);
453        printf("\nInd:  ");
454        for (j = 0; j < m_Mod[i]->getPartCount(); j++)
455                printf("%3i ", (int)m_aDegrees[i][j][0]);
456        printf("\nDeg:  ");
457        for (j = 0; j < m_Mod[i]->getPartCount(); j++)
458                printf("%3i ", (int)m_aDegrees[i][j][1]);
459        printf("\nNIt:  ");
460        for (j = 0; j < m_Mod[i]->getPartCount(); j++)
461                printf("%3i ", (int)m_aDegrees[i][j][2]);
462        printf("\nNeu:  ");
463        for (j = 0; j < m_Mod[i]->getPartCount(); j++)
464                printf("%3i ", (int)m_aDegrees[i][j][3]);
465        printf("\n");
466}
467
468/** Prints one array of ints. Debug method.
469        @param array Base pointer of the array.
470        @param base First index of the array's element.
471        @param size Number of elements to print.
472        */
473void ModelSimil::_PrintArray(int *array, int base, int size)
474{
475        int i;
476        for (i = base; i < base + size; i++)
477        {
478                printf("%i ", array[i]);
479        }
480        printf("\n");
481}
482
483void ModelSimil::_PrintArrayDouble(double *array, int base, int size)
484{
485        int i;
486        for (i = base; i < base + size; i++)
487        {
488                printf("%f ", array[i]);
489        }
490        printf("\n");
491}
492
493/** Prints one array of parts neighbourhood.
494        @param index of organism
495        */
496void ModelSimil::_PrintNeighbourhood(int o)
497{
498        assert(m_Neighbours[o] != 0);
499        printf("Neighbourhhod of organism %i\n", o);
500        int size = m_Mod[o]->getPartCount();
501        for (int i = 0; i < size; i++)
502        {
503                for (int j = 0; j < size; j++)
504                {
505                        printf("%i ", m_Neighbours[o][i][j]);
506                }
507                printf("\n");
508        }
509}
510
511/** Prints one array of parts fuzzy neighbourhood.
512        @param index of organism
513        */
514void ModelSimil::_PrintFuzzyNeighbourhood(int o)
515{
516        assert(m_fuzzyNeighb[o] != NULL);
517        printf("Fuzzy neighbourhhod of organism %i\n", o);
518        int size = m_Mod[o]->getPartCount();
519        for (int i = 0; i < size; i++)
520        {
521                for (int j = 0; j < fuzzyDepth; j++)
522                {
523                        printf("%f ", m_fuzzyNeighb[o][i][j]);
524                }
525                printf("\n");
526        }
527}
528
529Model* ModelSimil::newModel(const Geno *g)
530{
531        if (g == NULL)
532        {
533                logPrintf("ModelSimil", "newModel", LOG_ERROR, "NULL genotype pointer");
534                return NULL;
535        }
536        Model *m = new Model(*g);
537        if (!m->isValid())
538        {
539                logPrintf("ModelSimil", "newModel", LOG_ERROR, "Invalid model for the genotype of '%s'", g->getName().c_str());
540                delete m;
541                return NULL;
542        }
543        return m;
544}
545
546
547/** Creates arrays holding information about organisms' Parts (m_aDegrees) andm_Neigh
548        fills them with initial data (original indices and zeros).
549        Assumptions:
550        - Models (m_Mod) are created and available.
551        */
552int ModelSimil::ModelSimil::CreatePartInfoTables()
553{
554        // check assumptions about models
555        assert((m_Mod[0] != NULL) && (m_Mod[1] != NULL));
556        assert(m_Mod[0]->isValid() && m_Mod[1]->isValid());
557
558        int i, j, partCount;
559        // utwórz tablice na informacje o stopniach wierzchołków i liczbie neuroitems
560        for (i = 0; i < 2; i++)
561        {
562                partCount = m_Mod[i]->getPartCount();
563                // utworz i wypelnij tablice dla Parts wartosciami poczatkowymi
564                m_aDegrees[i] = new TDN[partCount];
565
566                if (isFuzzy)
567                {
568                        m_Neighbours[i] = new int*[partCount];
569                        m_fuzzyNeighb[i] = new float*[partCount];
570                }
571
572                if (m_aDegrees[i] != NULL && ((!isFuzzy) || (m_Neighbours[i] != NULL && m_fuzzyNeighb[i] != NULL)))
573                {
574                        // wypelnij tablice zgodnie z sensem TDN[0] - orginalny index
575                        // TDN[1], TDN[2], TDN[3] - zerami
576                        DB(printf("m_aDegrees[%i]: %p\n", i, m_aDegrees[i]);)
577                                for (j = 0; j < partCount; j++)
578                                {
579                                        m_aDegrees[i][j][0] = j;
580                                        m_aDegrees[i][j][1] = 0;
581                                        m_aDegrees[i][j][2] = 0;
582                                        m_aDegrees[i][j][3] = 0;
583                                        m_aDegrees[i][j][4] = 0;
584
585                                        // sprawdz, czy nie piszemy po jakims szalonym miejscu pamieci
586                                        assert(m_aDegrees[i][j] != NULL);
587
588                                        if (isFuzzy)
589                                        {
590                                                m_Neighbours[i][j] = new int[partCount];
591                                                for (int k = 0; k < partCount; k++)
592                                                {
593                                                        m_Neighbours[i][j][k] = 0;
594                                                }
595
596                                                m_fuzzyNeighb[i][j] = new float[fuzzyDepth];
597                                                for (int k = 0; k < fuzzyDepth; k++)
598                                                {
599                                                        m_fuzzyNeighb[i][j][k] = 0;
600                                                }
601
602                                                assert(m_Neighbours[i][j] != NULL);
603                                                assert(m_fuzzyNeighb[i][j] != NULL);
604                                        }
605
606                                }
607                }
608                else
609                {
610                        logPrintf("ModelSimil", "CreatePartInfoTables", LOG_ERROR, "Not enough memory?");
611                        return 0;
612                }
613                // utworz tablice dla pozycji 3D Parts (wielkosc tablicy: liczba Parts organizmu)
614                m_aPositions[i] = new Pt3D[m_Mod[i]->getPartCount()];
615                assert(m_aPositions[i] != NULL);
616                // wypelnij tablice OnJoints i Anywhere wartościami początkowymi
617                // OnJoint
618                m_aOnJoint[i][0] = 0;
619                m_aOnJoint[i][1] = 0;
620                m_aOnJoint[i][2] = 0;
621                m_aOnJoint[i][3] = 0;
622                // Anywhere
623                m_aAnywhere[i][0] = 0;
624                m_aAnywhere[i][1] = 0;
625                m_aAnywhere[i][2] = 0;
626                m_aAnywhere[i][3] = 0;
627        }
628        return 1;
629}
630
631/** Computes degrees of Parts of both organisms. Fills in the m_aDegrees arrays
632        with proper information about degrees.
633        Assumptions:
634        - Models (m_Mod) are created and available.
635        - Arrays m_aDegrees are created.
636        */
637int ModelSimil::CountPartDegrees()
638{
639        // sprawdz zalozenie - o modelach
640        assert((m_Mod[0] != NULL) && (m_Mod[1] != NULL));
641        assert(m_Mod[0]->isValid() && m_Mod[1]->isValid());
642
643        // sprawdz zalozenie - o tablicach
644        assert(m_aDegrees[0] != NULL);
645        assert(m_aDegrees[1] != NULL);
646
647        Part *P1, *P2;
648        int i, j, i1, i2;
649
650        // dla obu stworzen oblicz stopnie wierzcholkow
651        for (i = 0; i < 2; i++)
652        {
653                // dla wszystkich jointow
654                for (j = 0; j < m_Mod[i]->getJointCount(); j++)
655                {
656                        // pobierz kolejny Joint
657                        Joint *J = m_Mod[i]->getJoint(j);
658                        // wez jego konce
659                        P1 = J->part1;
660                        P2 = J->part2;
661                        // znajdz ich indeksy w Modelu (indeksy orginalne)
662                        i1 = m_Mod[i]->findPart(P1);
663                        i2 = m_Mod[i]->findPart(P2);
664                        // zwieksz stopien odpowiednich Parts
665                        m_aDegrees[i][i1][DEGREE]++;
666                        m_aDegrees[i][i2][DEGREE]++;
667                        m_aDegrees[i][i1][FUZZ_DEG]++;
668                        m_aDegrees[i][i2][FUZZ_DEG]++;
669                        if (isFuzzy)
670                        {
671                                m_Neighbours[i][i1][i2] = 1;
672                                m_Neighbours[i][i2][i1] = 1;
673                        }
674                }
675                // dla elementow nie osadzonych na Parts (OnJoint, Anywhere) -
676                // stopnie wierzchołka są już ustalone na zero
677        }
678
679        if (isFuzzy)
680        {
681                CountFuzzyNeighb();
682        }
683
684        return 1;
685}
686
687void ModelSimil::GetNeighbIndexes(int mod, int partInd, std::vector<int> &indexes)
688{
689        indexes.clear();
690        int i, size = m_Mod[mod]->getPartCount();
691
692        for (i = 0; i < size; i++)
693        {
694                if (m_Neighbours[mod][partInd][i] > 0)
695                {
696                        indexes.push_back(i);
697                }
698        }
699}
700
701int cmpFuzzyRows(const void *pa, const void *pb)
702{
703        float **a = (float**)pa;
704        float **b = (float**)pb;
705
706
707        for (int i = 1; i < fuzDepth; i++)
708        {
709                if (a[0][i] > b[0][i])
710                {
711
712                        return -1;
713                }
714                if (a[0][i] < b[0][i])
715                {
716
717                        return 1;
718                }
719        }
720
721        return 0;
722}
723
724void ModelSimil::FuzzyOrder()
725{
726        int i, depth, partInd, prevPartInd, partCount;
727        for (int mod = 0; mod < 2; mod++)
728        {
729                partCount = m_Mod[mod]->getPartCount();
730                partInd = m_fuzzyNeighb[mod][partCount - 1][0];
731                m_aDegrees[mod][partInd][FUZZ_DEG] = 0;
732
733                for (i = (partCount - 2); i >= 0; i--)
734                {
735                        prevPartInd = partInd;
736                        partInd = m_fuzzyNeighb[mod][i][0];
737                        m_aDegrees[mod][partInd][FUZZ_DEG] = m_aDegrees[mod][prevPartInd][FUZZ_DEG];
738                        for (depth = 1; depth < fuzzyDepth; depth++)
739                        {
740                                if (m_fuzzyNeighb[mod][i][depth] != m_fuzzyNeighb[mod][i + 1][depth])
741                                {
742                                        m_aDegrees[mod][partInd][FUZZ_DEG]++;
743                                        break;
744                                }
745                        }
746                }
747        }
748}
749
750//sort according to fuzzy degree
751void ModelSimil::SortFuzzyNeighb()
752{
753        fuzDepth = fuzzyDepth;
754        for (int mod = 0; mod < 2; mod++)
755        {
756                std::qsort(m_fuzzyNeighb[mod], (size_t)m_Mod[mod]->getPartCount(), sizeof(m_fuzzyNeighb[mod][0]), cmpFuzzyRows);
757        }
758}
759
760//computes fuzzy vertex degree
761void ModelSimil::CountFuzzyNeighb()
762{
763        assert(m_aDegrees[0] != NULL);
764        assert(m_aDegrees[1] != NULL);
765
766        assert(m_Neighbours[0] != NULL);
767        assert(m_Neighbours[1] != NULL);
768
769        assert(m_fuzzyNeighb[0] != NULL);
770        assert(m_fuzzyNeighb[1] != NULL);
771
772        std::vector<int> nIndexes;
773        float newDeg = 0;
774
775        for (int mod = 0; mod < 2; mod++)
776        {
777                int partCount = m_Mod[mod]->getPartCount();
778
779                for (int depth = 0; depth < fuzzyDepth; depth++)
780                {
781                        //use first column for storing indices
782                        for (int partInd = 0; partInd < partCount; partInd++)
783                        {
784                                if (depth == 0)
785                                {
786                                        m_fuzzyNeighb[mod][partInd][depth] = partInd;
787                                }
788                                else if (depth == 1)
789                                {
790                                        m_fuzzyNeighb[mod][partInd][depth] = m_aDegrees[mod][partInd][DEGREE];
791                                }
792                                else
793                                {
794                                        GetNeighbIndexes(mod, partInd, nIndexes);
795                                        for (unsigned int k = 0; k < nIndexes.size(); k++)
796                                        {
797                                                newDeg += m_fuzzyNeighb[mod][nIndexes.at(k)][depth - 1];
798                                        }
799                                        newDeg /= nIndexes.size();
800                                        m_fuzzyNeighb[mod][partInd][depth] = newDeg;
801                                        for (int mod = 0; mod < 2; mod++)
802                                        {
803                                                int partCount = m_Mod[mod]->getPartCount();
804                                                for (int i = partCount - 1; i >= 0; i--)
805                                                {
806
807                                                }
808                                        }
809                                        newDeg = 0;
810                                }
811                        }
812                }
813        }
814
815        SortFuzzyNeighb();
816        FuzzyOrder();
817}
818
819/** Gets information about Parts' positions in 3D from models into the arrays
820                m_aPositions.
821                Assumptions:
822                - Models (m_Mod) are created and available.
823                - Arrays m_aPositions are created.
824                @return 1 if success, otherwise 0.
825                */
826int ModelSimil::GetPartPositions()
827{
828        // sprawdz zalozenie - o modelach
829        assert((m_Mod[0] != NULL) && (m_Mod[1] != NULL));
830        assert(m_Mod[0]->isValid() && m_Mod[1]->isValid());
831
832        // sprawdz zalozenie - o tablicach m_aPositions
833        assert(m_aPositions[0] != NULL);
834        assert(m_aPositions[1] != NULL);
835
836        // dla każdego stworzenia skopiuj informację o polozeniu jego Parts
837        // do tablic m_aPositions
838        Part *pPart;
839        int iMod; // licznik modeli (organizmow)
840        int iPart; // licznik Parts
841        for (iMod = 0; iMod < 2; iMod++)
842        {
843                // dla każdego z modeli iMod
844                for (iPart = 0; iPart < m_Mod[iMod]->getPartCount(); iPart++)
845                {
846                        // dla każdego iPart organizmu iMod
847                        // pobierz tego Part
848                        pPart = m_Mod[iMod]->getPart(iPart);
849                        // zapamietaj jego pozycje 3D w tablicy m_aPositions
850                        m_aPositions[iMod][iPart].x = pPart->p.x;
851                        m_aPositions[iMod][iPart].y = pPart->p.y;
852                        m_aPositions[iMod][iPart].z = pPart->p.z;
853                }
854        }
855
856        return 1;
857}
858
859/** Computes numbers of neurons and neurons' inputs for each Part of each
860        organisms and fills in the m_aDegrees array.
861        Assumptions:
862        - Models (m_Mod) are created and available.
863        - Arrays m_aDegrees are created.
864        */
865int ModelSimil::CountPartNeurons()
866{
867        // sprawdz zalozenie - o modelach
868        assert((m_Mod[0] != NULL) && (m_Mod[1] != NULL));
869        assert(m_Mod[0]->isValid() && m_Mod[1]->isValid());
870
871        // sprawdz zalozenie - o tablicach
872        assert(m_aDegrees[0] != NULL);
873        assert(m_aDegrees[1] != NULL);
874
875        Part *P1;
876        Joint *J1;
877        int i, j, i2, neuro_connections;
878
879        // dla obu stworzen oblicz liczbe Neurons + connections dla Parts
880        // a takze dla OnJoints i Anywhere
881        for (i = 0; i < 2; i++)
882        {
883                for (j = 0; j < m_Mod[i]->getNeuroCount(); j++)
884                {
885                        // pobierz kolejny Neuron
886                        Neuro *N = m_Mod[i]->getNeuro(j);
887                        // policz liczbe jego wejść i jego samego tez
888                        // czy warto w ogole liczyc polaczenia...? co to da/spowoduje?
889                        neuro_connections = N->getInputCount() + 1;
890                        // wez Part, na ktorym jest Neuron
891                        P1 = N->getPart();
892                        if (P1)
893                        {
894                                // dla neuronow osadzonych na Partach
895                                i2 = m_Mod[i]->findPart(P1); // znajdz indeks Part w Modelu
896                                m_aDegrees[i][i2][2] += neuro_connections; // zwieksz liczbe Connections+Neurons dla tego Part (TDN[2])
897                                m_aDegrees[i][i2][3]++; // zwieksz liczbe Neurons dla tego Part (TDN[3])
898                        }
899                        else
900                        {
901                                // dla neuronow nie osadzonych na partach
902                                J1 = N->getJoint();
903                                if (J1)
904                                {
905                                        // dla tych na Jointach
906                                        m_aOnJoint[i][2] += neuro_connections; // zwieksz liczbe Connections+Neurons
907                                        m_aOnJoint[i][3]++; // zwieksz liczbe Neurons
908                                }
909                                else
910                                {
911                                        // dla tych "gdziekolwiek"
912                                        m_aAnywhere[i][2] += neuro_connections; // zwieksz liczbe Connections+Neurons
913                                        m_aAnywhere[i][3]++; // zwieksz liczbe Neurons
914                                }
915                        }
916                }
917        }
918        return 1;
919}
920
921/** Sorts arrays m_aDegrees (for each organism) by Part's degree, and then by
922        number of neural connections and neurons in groups of Parts with the same
923        degree.
924        Assumptions:
925        - Models (m_Mod) are created and available.
926        - Arrays m_aDegrees are created.
927        @saeDegrees, CompareItemNo
928        */
929int ModelSimil::SortPartInfoTables()
930{
931        // sprawdz zalozenie - o modelach
932        assert((m_Mod[0] != NULL) && (m_Mod[1] != NULL));
933        assert(m_Mod[0]->isValid() && m_Mod[1]->isValid());
934
935        // sprawdz zalozenie - o tablicach
936        assert(m_aDegrees[0] != NULL);
937        assert(m_aDegrees[1] != NULL);
938
939        int i;
940        int(*pfDegreeFunction) (const void*, const void*) = NULL;
941        pfDegreeFunction = isFuzzy ? &CompareFuzzyDegrees : &CompareDegrees;
942        // sortowanie obu tablic wg stopni punktów - TDN[1]
943        for (i = 0; i < 2; i++)
944        {
945                DB(_PrintDegrees(i));
946                std::qsort(m_aDegrees[i], (size_t)(m_Mod[i]->getPartCount()),
947                        sizeof(TDN), pfDegreeFunction);
948                DB(_PrintDegrees(i));
949        }
950
951        // sprawdzenie wartosci parametru
952        DB(i = sizeof(TDN);)
953                int degreeType = isFuzzy ? FUZZ_DEG : DEGREE;
954
955        // sortowanie obu tablic m_aDegrees wedlug liczby neuronów i
956        // czesci neuronu - ale w obrebie grup o tym samym stopniu
957        for (i = 0; i < 2; i++)
958        {
959                int iPocz = 0;
960                int iDeg, iNewDeg, iPartCount, j;
961                // stopien pierwszego punktu w tablicy Degrees  odniesienie
962                iDeg = m_aDegrees[i][0][degreeType];
963                iPartCount = m_Mod[i]->getPartCount();
964                // po kolei dla kazdego punktu w organizmie
965                for (j = 0; j <= iPartCount; j++)
966                {
967                        // sprawdz stopien punktu (lub nadaj 0 - gdy juz koniec tablicy)
968                        //                      iNewDeg = (j != iPartCount) ? m_aDegrees[i][j][1] : 0;
969                        // usunieto stara wersje porownania!!! wprowadzono znak porownania <
970
971                        iNewDeg = (j < iPartCount) ? m_aDegrees[i][j][degreeType] : 0;
972                        // skoro tablice sa posortowane wg stopni, to mamy na pewno taka kolejnosc
973                        assert(iNewDeg <= iDeg);
974                        if (iNewDeg != iDeg)
975                        {
976                                // gdy znaleziono koniec grupy o tym samym stopniu
977                                // sortuj po liczbie neuronow w obrebie grupy
978                                DB(_PrintDegrees(i));
979                                DB(printf("qsort( poczatek=%i, rozmiar=%i, sizeof(TDN)=%i)\n", iPocz, (j - iPocz), sizeof(TDN));)
980                                        // wyswietlamy z jedna komorka po zakonczeniu tablicy
981                                        DB(_PrintArray(m_aDegrees[i][iPocz], 0, (j - iPocz) * 4);)
982
983                                        std::qsort(m_aDegrees[i][iPocz], (size_t)(j - iPocz),
984                                                sizeof(TDN), ModelSimil::CompareConnsNo);
985                                DB(_PrintDegrees(i));
986                                // wyswietlamy z jedna komorka po zakonczeniu tablicy
987                                DB(_PrintArray(m_aDegrees[i][iPocz], 0, (j - iPocz) * 4);)
988                                        // rozpocznij nowa grupe
989                                        iPocz = j;
990                                iDeg = iNewDeg;
991                        }
992                }
993        }
994        return 1;
995}
996
997
998/** Prints the state of the matching object. Debug method.
999 */
1000void ModelSimil::_PrintPartsMatching()
1001{
1002        // assure that matching exists
1003        assert(m_pMatching != NULL);
1004
1005        printf("Parts matching:\n");
1006        m_pMatching->PrintMatching();
1007}
1008
1009void ModelSimil::ComputeMatching()
1010{
1011        // uniwersalne liczniki
1012        int i, j;
1013
1014        assert(m_pMatching != NULL);
1015        assert(m_pMatching->IsEmpty() == true);
1016
1017        // rozpoczynamy etap dopasowywania Parts w organizmach
1018        // czy dopasowano już wszystkie Parts?
1019        int iCzyDopasowane = 0;
1020        // koniec grupy aktualnie dopasowywanej w każdym organizmie
1021        int aiKoniecGrupyDopasowania[2] = { 0, 0 };
1022        // koniec grupy już w całości dopasowanej
1023        // (Pomiedzy tymi dwoma indeksami znajduja sie Parts w tablicy
1024        // m_aDegrees, ktore moga byc dopasowywane (tam nadal moga
1025        // byc tez dopasowane - ale nie musi to byc w sposob
1026        // ciagly)
1027        int aiKoniecPierwszejGrupy[2] = { 0, 0 };
1028        // Tablica przechowująca odległości poszczególnych Parts z aktualnych
1029        // grup dopasowania. Rozmiar - prostokąt o bokach równych liczbie elementów w
1030        // dopasowywanych aktualnie grupach. Pierwszy wymiar - pierwszy organizm.
1031        // Drugi wymiar - drugi organizm (nie zależy to od tego, który jest mniejszy).
1032        // Wliczane w rozmiar tablicy są nawet już dopasowane elementy - tablice
1033        // paiCzyDopasowany pamiętają stan dopasowania tych elementów.
1034        typedef double *TPDouble;
1035        double **aadOdleglosciParts;
1036        // dwie tablice okreslajace Parts, ktore moga byc do siebie dopasowywane
1037        // rozmiary: [0] - aiRozmiarCalychGrup[1]
1038        //                       [1] - aiRozmiarCalychGrup[0]
1039        std::vector<bool> *apvbCzyMinimalnaOdleglosc[2];
1040        // rozmiar aktualnie dopasowywanej grupy w odpowiednim organizmie (tylko elementy
1041        // jeszcze niedopasowane).
1042        int aiRozmiarGrupy[2];
1043        // rozmiar aktualnie dopasowywanych grup w odpowiednim organizmie (włączone są
1044        // w to również elementy już dopasowane).
1045        int aiRozmiarCalychGrup[2] = { 0, 0 };
1046
1047        // utworzenie tablicy rozmiarow
1048        for (i = 0; i < 2; i++)
1049        {
1050                m_aiPartCount[i] = m_Mod[i]->getPartCount();
1051        }
1052
1053        // DOPASOWYWANIE PARTS
1054        while (!iCzyDopasowane)
1055        {
1056                // znajdz konce obu grup aktualnie dopasowywanych w obu organizmach
1057                for (i = 0; i < 2; i++)
1058                {
1059                        // czyli poszukaj miejsca zmiany stopnia lub konca tablicy
1060                        for (j = aiKoniecPierwszejGrupy[i] + 1; j < m_aiPartCount[i]; j++)
1061                        {
1062                                if (m_aDegrees[i][j][DEGREE] < m_aDegrees[i][j - 1][DEGREE])
1063                                {
1064                                        break;
1065                                }
1066                        }
1067                        aiKoniecGrupyDopasowania[i] = j;
1068
1069                        // sprawdz poprawnosc tego indeksu
1070                        assert((aiKoniecGrupyDopasowania[i] > 0) &&
1071                                (aiKoniecGrupyDopasowania[i] <= m_aiPartCount[i]));
1072
1073                        // oblicz rozmiary całych grup - łącznie z dopasowanymi już elementami
1074                        aiRozmiarCalychGrup[i] = aiKoniecGrupyDopasowania[i] -
1075                                aiKoniecPierwszejGrupy[i];
1076
1077                        // sprawdz teraz rozmiar tej grupy w sensie liczby niedopasowanzch
1078                        // nie moze to byc puste!
1079                        aiRozmiarGrupy[i] = 0;
1080                        for (j = aiKoniecPierwszejGrupy[i]; j < aiKoniecGrupyDopasowania[i]; j++)
1081                        {
1082                                // od poczatku do konca grupy
1083                                if (!m_pMatching->IsMatched(i, j))
1084                                {
1085                                        // jesli niedopasowany, to zwieksz licznik
1086                                        aiRozmiarGrupy[i]++;
1087                                }
1088                        }
1089                        // grupa nie moze byc pusta!
1090                        assert(aiRozmiarGrupy[i] > 0);
1091                }
1092
1093                // DOPASOWYWANIE PARTS Z GRUP
1094
1095                // stworzenie tablicy odległości lokalnych
1096                // stwórz pierwszy wymiar - wg rozmiaru zerowego organizmu
1097                aadOdleglosciParts = new TPDouble[aiRozmiarCalychGrup[0]];
1098                assert(aadOdleglosciParts != NULL);
1099                // stwórz drugi wymiar - wg rozmiaru drugiego organizmu
1100                for (i = 0; i < aiRozmiarCalychGrup[0]; i++)
1101                {
1102                        aadOdleglosciParts[i] = new double[aiRozmiarCalychGrup[1]];
1103                        assert(aadOdleglosciParts[i] != NULL);
1104                }
1105
1106                // stworzenie tablic mozliwosci dopasowania (indykatorow minimalnej odleglosci)
1107                apvbCzyMinimalnaOdleglosc[0] = new std::vector<bool>(aiRozmiarCalychGrup[1], false);
1108                apvbCzyMinimalnaOdleglosc[1] = new std::vector<bool>(aiRozmiarCalychGrup[0], false);
1109                // sprawdz stworzenie tablic
1110                assert(apvbCzyMinimalnaOdleglosc[0] != NULL);
1111                assert(apvbCzyMinimalnaOdleglosc[1] != NULL);
1112
1113                // wypełnienie elementów macierzy (i,j) odległościami pomiędzy
1114                // odpowiednimi Parts: (i) w organizmie 0 i (j) w organizmie 1.
1115                // UWAGA! Uwzględniamy tylko te Parts, które nie są jeszcze dopasowane
1116                // (reszta to byłaby po prostu strata czasu).
1117                int iDeg, iNeu; // ilościowe określenie różnic stopnia, liczby neuronów i połączeń Parts
1118                //int iNIt;
1119                double dGeo; // ilościowe określenie różnic geometrycznych pomiędzy Parts
1120                // indeksy konkretnych Parts - indeksy sa ZERO-BASED, choć właściwy dostep
1121                // do informacji o Part wymaga dodania aiKoniecPierwszejGrupy[]
1122                // tylko aadOdleglosciParts[][] indeksuje sie bezposrednio zawartoscia iIndex[]
1123                int iIndex[2];
1124                int iPartIndex[2] = { -1, -1 }; // at [iModel]: original index of a Part for the specified model (iModel)
1125
1126                // debug - wypisz zakres dopasowywanych indeksow
1127                DB(printf("Organizm 0: grupa: (%2i, %2i)\n", aiKoniecPierwszejGrupy[0],
1128                        aiKoniecGrupyDopasowania[0]);)
1129                        DB(printf("Organizm 1: grupa: (%2i, %2i)\n", aiKoniecPierwszejGrupy[1],
1130                                aiKoniecGrupyDopasowania[1]);)
1131
1132                        for (i = 0; i < aiRozmiarCalychGrup[0]; i++)
1133                        {
1134
1135                                // iterujemy i - Parts organizmu 0
1136                                // (indeks podstawowy - aiKoniecPierwszejGrupy[0])
1137
1138                                if (!(m_pMatching->IsMatched(0, aiKoniecPierwszejGrupy[0] + i)))
1139                                {
1140                                        // interesuja nas tylko te niedopasowane jeszcze (i)
1141                                        for (j = 0; j < aiRozmiarCalychGrup[1]; j++)
1142                                        {
1143
1144                                                // iterujemy j - Parts organizmu 1
1145                                                // (indeks podstawowy - aiKoniecPierwszejGrupy[1])
1146
1147                                                if (!(m_pMatching->IsMatched(1, aiKoniecPierwszejGrupy[1] + j)))
1148                                                {
1149                                                        // interesuja nas tylko te niedopasowane jeszcze (j)
1150                                                        // teraz obliczymy lokalne różnice pomiędzy Parts
1151                                                        iDeg = abs(m_aDegrees[0][aiKoniecPierwszejGrupy[0] + i][1]
1152                                                                - m_aDegrees[1][aiKoniecPierwszejGrupy[1] + j][1]);
1153
1154                                                        //iNit currently is not a component of distance measure           
1155                                                        //iNIt = abs(m_aDegrees[0][ aiKoniecPierwszejGrupy[0] + i ][2]
1156                                                        //           - m_aDegrees[1][ aiKoniecPierwszejGrupy[1] + j ][2]);
1157
1158                                                        iNeu = abs(m_aDegrees[0][aiKoniecPierwszejGrupy[0] + i][3]
1159                                                                - m_aDegrees[1][aiKoniecPierwszejGrupy[1] + j][3]);
1160
1161                                                        // obliczenie także lokalnych różnic geometrycznych pomiędzy Parts
1162                                                        // find original indices of Parts for both of the models
1163                                                        iPartIndex[0] = m_aDegrees[0][aiKoniecPierwszejGrupy[0] + i][0];
1164                                                        iPartIndex[1] = m_aDegrees[1][aiKoniecPierwszejGrupy[1] + j][0];
1165                                                        // now compute the geometrical distance of these Parts (use m_aPositions
1166                                                        // which should be computed by SVD)
1167                                                        Pt3D Part0Pos(m_aPositions[0][iPartIndex[0]]);
1168                                                        Pt3D Part1Pos(m_aPositions[1][iPartIndex[1]]);
1169                                                        dGeo = m_adFactors[3] == 0 ? 0 : Part0Pos.distanceTo(Part1Pos); //no need to compute distane when m_dDG weight is 0
1170
1171                                                        // tutaj prawdopodobnie należy jeszcze dodać sprawdzanie
1172                                                        // identyczności pozostałych własności (oprócz geometrii)
1173                                                        // - żeby móc stwierdzić w ogóle identyczność Parts
1174
1175                                                        // i ostateczna odleglosc indukowana przez te roznice
1176                                                        // (tutaj nie ma różnicy w liczbie wszystkich wierzchołków)
1177                                                        aadOdleglosciParts[i][j] = m_adFactors[1] * double(iDeg) +
1178                                                                m_adFactors[2] * double(iNeu) +
1179                                                                m_adFactors[3] * dGeo;
1180                                                        DB(printf("Parts Distance (%2i,%2i) = %.3lf\n", aiKoniecPierwszejGrupy[0] + i,
1181                                                                aiKoniecPierwszejGrupy[1] + j, aadOdleglosciParts[i][j]);)
1182                                                                DB(printf("Parts geometrical distance = %.20lf\n", dGeo);)
1183                                                                DB(printf("Pos0: (%.3lf %.3lf %.3lf)\n", Part0Pos.x, Part0Pos.y, Part0Pos.z);)
1184                                                                DB(printf("Pos1: (%.3lf %.3lf %.3lf)\n", Part1Pos.x, Part1Pos.y, Part1Pos.z);)
1185                                                }
1186                                        }
1187                                }
1188                        }
1189
1190                // tutaj - sprawdzic tylko, czy tablica odleglosci lokalnych jest poprawnie obliczona
1191
1192                // WYKORZYSTANIE TABLICY ODLEGLOSCI DO BUDOWY DOPASOWANIA
1193
1194                // trzeba raczej iterować aż do zebrania wszystkich możliwych dopasowań w grupie
1195                // dlatego wprowadzamy dodatkowa zmienna - czy skonczyla sie juz grupa
1196                bool bCzyKoniecGrupy = false;
1197                while (!bCzyKoniecGrupy)
1198                {
1199                        for (i = 0; i < 2; i++)
1200                        {
1201                                // iterujemy (i) po organizmach
1202                                // szukamy najpierw jakiegoś niedopasowanego jeszcze Part w organizmach
1203
1204                                // zakładamy, że nie ma takiego Part
1205                                iIndex[i] = -1;
1206
1207                                for (j = 0; j < aiRozmiarCalychGrup[i]; j++)
1208                                {
1209                                        // iterujemy (j) - Parts organizmu (i)
1210                                        // (indeks podstawowy - aiKoniecPierwszejGrupy[0])
1211                                        if (!(m_pMatching->IsMatched(i, aiKoniecPierwszejGrupy[i] + j)))
1212                                        {
1213                                                // gdy mamy w tej grupie jakis niedopasowany element, to ustawiamy
1214                                                // iIndex[i] (chcemy w zasadzie pierwszy taki)
1215                                                iIndex[i] = j;
1216                                                break;
1217                                        }
1218                                }
1219
1220                                // sprawdzamy, czy w ogole znaleziono taki Part
1221                                if (iIndex[i] < 0)
1222                                {
1223                                        // gdy nie znaleziono takiego Part - mamy koniec dopasowywania w
1224                                        // tych grupach
1225                                        bCzyKoniecGrupy = true;
1226                                }
1227                                // sprawdz poprawnosc indeksu niedopasowanego Part - musi byc w aktualnej grupie
1228                                assert((iIndex[i] >= -1) && (iIndex[i] < aiRozmiarCalychGrup[i]));
1229                        }
1230
1231
1232                        // teraz mamy sytuacje:
1233                        // - mamy w iIndex[0] i iIndex[1] indeksy pierwszych niedopasowanych Part
1234                        //              w organizmach, albo
1235                        // - nie ma w ogóle już czego dopasowywać (należy przejść do innej grupy)
1236                        if (!bCzyKoniecGrupy)
1237                        {
1238                                // gdy nie ma jeszcze końca żadnej z grup - możemy dopasowywać
1239                                // najpierw wyszukujemy w tablicy minimum odległości od tych
1240                                // wyznaczonych Parts
1241
1242                                // najpierw wyczyscimy wektory potencjalnych dopasowan
1243                                // dla organizmu 1 (o rozmiarze grupy z 0)
1244                                for (i = 0; i < aiRozmiarCalychGrup[0]; i++)
1245                                {
1246                                        apvbCzyMinimalnaOdleglosc[1]->operator[](i) = false;
1247                                }
1248                                // dla organizmu 0 (o rozmiarze grup z 1)
1249                                for (i = 0; i < aiRozmiarCalychGrup[1]; i++)
1250                                {
1251                                        apvbCzyMinimalnaOdleglosc[0]->operator[](i) = false;
1252                                }
1253
1254                                // szukanie minimum dla Part o indeksie iIndex[0] w organizmie 0
1255                                // wsrod niedopasowanych Parts z organizmu 1
1256                                // zakładamy, że nie znaleŸliśmy jeszcze minimum
1257                                double dMinimum = HUGE_VAL;
1258                                for (i = 0; i < aiRozmiarCalychGrup[1]; i++)
1259                                {
1260                                        if (!(m_pMatching->IsMatched(1, aiKoniecPierwszejGrupy[1] + i)))
1261                                        {
1262
1263                                                // szukamy minimum obliczonej lokalnej odleglosci tylko wsrod
1264                                                // niedopasowanych jeszcze Parts
1265                                                if (aadOdleglosciParts[iIndex[0]][i] < dMinimum)
1266                                                {
1267                                                        dMinimum = aadOdleglosciParts[iIndex[0]][i];
1268                                                }
1269
1270                                                // przy okazji - sprawdz, czy HUGE_VAL jest rzeczywiscie max dla double
1271                                                assert(aadOdleglosciParts[iIndex[0]][i] < HUGE_VAL);
1272                                        }
1273                                }
1274                                // sprawdz, czy minimum znaleziono - musi takie byc, bo jest cos niedopasowanego
1275                                assert((dMinimum >= 0.0) && (dMinimum < HUGE_VAL));
1276
1277                                // teraz zaznaczamy w tablicy te wszystkie Parts, ktore maja
1278                                // rzeczywiscie te minimalna odleglosc do Part iIndex[0] w organizmie 0
1279                                for (i = 0; i < aiRozmiarCalychGrup[1]; i++)
1280                                {
1281                                        if (!(m_pMatching->IsMatched(1, aiKoniecPierwszejGrupy[1] + i)))
1282                                        {
1283                                                if (aadOdleglosciParts[iIndex[0]][i] == dMinimum)
1284                                                {
1285                                                        // jesli w danym miejscu tablicy odleglosci jest faktyczne
1286                                                        // minimum odleglosci, to mamy potencjalna pare dla iIndex[0]
1287                                                        apvbCzyMinimalnaOdleglosc[0]->operator[](i) = true;
1288                                                }
1289
1290                                                // sprawdz poprawnosc znalezionego wczesniej minimum
1291                                                assert(aadOdleglosciParts[iIndex[0]][i] >= dMinimum);
1292                                        }
1293                                }
1294
1295                                // podobnie szukamy minimum dla Part o indeksie iIndex[1] w organizmie 1
1296                                // wsrod niedopasowanych Parts z organizmu 0
1297                                dMinimum = HUGE_VAL;
1298                                for (i = 0; i < aiRozmiarCalychGrup[0]; i++)
1299                                {
1300                                        if (!(m_pMatching->IsMatched(0, aiKoniecPierwszejGrupy[0] + i)))
1301                                        {
1302                                                // szukamy minimum obliczonej lokalnej odleglosci tylko wsrod
1303                                                // niedopasowanych jeszcze Parts
1304                                                if (aadOdleglosciParts[i][iIndex[1]] < dMinimum)
1305                                                {
1306                                                        dMinimum = aadOdleglosciParts[i][iIndex[1]];
1307                                                }
1308                                                // przy okazji - sprawdz, czy HUGE_VAL jest rzeczywiscie max dla double
1309                                                assert(aadOdleglosciParts[i][iIndex[1]] < HUGE_VAL);
1310                                        }
1311                                }
1312                                // sprawdz, czy minimum znaleziono - musi takie byc, bo jest cos niedopasowanego
1313                                assert((dMinimum >= 0.0) && (dMinimum < HUGE_VAL));
1314
1315                                // teraz zaznaczamy w tablicy te wszystkie Parts, ktore maja
1316                                // rzeczywiscie te minimalne odleglosci do Part iIndex[1] w organizmie 1
1317                                for (i = 0; i < aiRozmiarCalychGrup[0]; i++)
1318                                {
1319                                        if (!(m_pMatching->IsMatched(0, aiKoniecPierwszejGrupy[0] + i)))
1320                                        {
1321                                                if (aadOdleglosciParts[i][iIndex[1]] == dMinimum)
1322                                                {
1323                                                        // jesli w danym miejscu tablicy odleglosci jest faktyczne
1324                                                        // minimum odleglosci, to mamy potencjalna pare dla iIndex[1]
1325                                                        apvbCzyMinimalnaOdleglosc[1]->operator[](i) = true;
1326                                                }
1327
1328                                                // sprawdz poprawnosc znalezionego wczesniej minimum
1329                                                assert(aadOdleglosciParts[i][iIndex[1]] >= dMinimum);
1330                                        }
1331                                }
1332
1333                                // teraz mamy juz poszukane potencjalne grupy dopasowania - musimy
1334                                // zdecydowac, co do czego dopasujemy!
1335                                // szukamy Part iIndex[0] posrod mozliwych do dopasowania dla Part iIndex[1]
1336                                // szukamy takze Part iIndex[1] posrod mozliwych do dopasowania dla Part iIndex[0]
1337                                bool PartZ1NaLiscie0 = apvbCzyMinimalnaOdleglosc[0]->operator[](iIndex[1]);
1338                                bool PartZ0NaLiscie1 = apvbCzyMinimalnaOdleglosc[1]->operator[](iIndex[0]);
1339
1340                                if (PartZ1NaLiscie0 && PartZ0NaLiscie1)
1341                                {
1342                                        // PRZYPADEK 1. Oba Parts maja sie wzajemnie na listach mozliwych
1343                                        // dopasowan.
1344                                        // AKCJA. Dopasowanie wzajemne do siebie.
1345
1346                                        m_pMatching->Match(0, aiKoniecPierwszejGrupy[0] + iIndex[0],
1347                                                1, aiKoniecPierwszejGrupy[1] + iIndex[1]);
1348
1349                                        // zmniejsz liczby niedopasowanych elementow w grupach
1350                                        aiRozmiarGrupy[0]--;
1351                                        aiRozmiarGrupy[1]--;
1352                                        // debug - co zostalo dopasowane
1353                                        DB(printf("Przypadek 1.: dopasowane Parts: (%2i, %2i)\n", aiKoniecPierwszejGrupy[0]
1354                                                + iIndex[0], aiKoniecPierwszejGrupy[1] + iIndex[1]);)
1355
1356                                }// PRZYPADEK 1.
1357                                else
1358                                        if (PartZ1NaLiscie0 || PartZ0NaLiscie1)
1359                                        {
1360                                                // PRZYPADEK 2. Tylko jeden z Parts ma drugiego na swojej liscie
1361                                                // mozliwych dopasowan
1362                                                // AKCJA. Dopasowanie jednego jest proste (tego, ktory nie ma
1363                                                // na swojej liscie drugiego). Dla tego drugiego nalezy powtorzyc
1364                                                // duza czesc obliczen (znalezc mu nowa mozliwa pare)
1365
1366                                                // indeks organizmu, ktorego Part nie ma dopasowywanego Part
1367                                                // z drugiego organizmu na swojej liscie
1368                                                int iBezDrugiego;
1369
1370                                                // okreslenie indeksu organizmu z dopasowywanym Part
1371                                                if (!PartZ1NaLiscie0)
1372                                                {
1373                                                        iBezDrugiego = 0;
1374                                                }
1375                                                else
1376                                                {
1377                                                        iBezDrugiego = 1;
1378                                                }
1379                                                // sprawdz indeks organizmu
1380                                                assert((iBezDrugiego == 0) || (iBezDrugiego == 1));
1381
1382                                                int iDopasowywany = -1;
1383                                                // poszukujemy pierwszego z listy dopasowania
1384                                                for (i = 0; i < aiRozmiarCalychGrup[1 - iBezDrugiego]; i++)
1385                                                {
1386                                                        if (apvbCzyMinimalnaOdleglosc[iBezDrugiego]->operator[](i))
1387                                                        {
1388                                                                iDopasowywany = i;
1389                                                                break;
1390                                                        }
1391                                                }
1392                                                // sprawdz poprawnosc indeksu dopasowywanego (musimy cos znalezc!)
1393                                                // nieujemny i w odpowiedniej grupie!
1394                                                assert((iDopasowywany >= 0) &&
1395                                                        (iDopasowywany < aiRozmiarCalychGrup[1 - iBezDrugiego]));
1396
1397                                                // znalezlismy 1. Part z listy dopasowania - dopasowujemy!
1398                                                m_pMatching->Match(
1399                                                        iBezDrugiego,
1400                                                        aiKoniecPierwszejGrupy[iBezDrugiego] + iIndex[iBezDrugiego],
1401                                                        1 - iBezDrugiego,
1402                                                        aiKoniecPierwszejGrupy[1 - iBezDrugiego] + iDopasowywany);
1403                                                DB(printf("Przypadek 2.1.: dopasowane Parts dopasowanie bez drugiego: (%2i, %2i)\n", aiKoniecPierwszejGrupy[iBezDrugiego] + iIndex[iBezDrugiego],
1404                                                        aiKoniecPierwszejGrupy[1 - iBezDrugiego] + iDopasowywany);)
1405
1406                                                        // zmniejsz liczby niedopasowanych elementow w grupach
1407                                                        aiRozmiarGrupy[0]--;
1408                                                aiRozmiarGrupy[1]--;
1409
1410                                                // sprawdz, czy jest szansa dopasowania tego Part z drugiej strony
1411                                                // (ktora miala mozliwosc dopasowania tego Part z poprzedniego organizmu)
1412                                                if ((aiRozmiarGrupy[0] > 0) && (aiRozmiarGrupy[1] > 0))
1413                                                {
1414                                                        // jesli grupy sie jeszcze nie wyczrpaly
1415                                                        // to jest mozliwosc dopasowania w organizmie
1416
1417                                                        int iZDrugim = 1 - iBezDrugiego;
1418                                                        // sprawdz indeks organizmu
1419                                                        assert((iZDrugim == 0) || (iZDrugim == 1));
1420
1421                                                        // bedziemy szukac minimum wsrod niedopasowanych - musimy wykasowac
1422                                                        // poprzednie obliczenia minimum
1423                                                        // dla organizmu 1 (o rozmiarze grupy z 0)
1424                                                        for (i = 0; i < aiRozmiarCalychGrup[0]; i++)
1425                                                        {
1426                                                                apvbCzyMinimalnaOdleglosc[1]->operator[](i) = false;
1427                                                        }
1428                                                        // dla organizmu 0 (o rozmiarze grup z 1)
1429                                                        for (i = 0; i < aiRozmiarCalychGrup[1]; i++)
1430                                                        {
1431                                                                apvbCzyMinimalnaOdleglosc[0]->operator[](i) = false;
1432                                                        }
1433
1434                                                        // szukamy na nowo minimum dla Part o indeksie iIndex[ iZDrugim ] w organizmie iZDrugim
1435                                                        // wsrod niedopasowanych Parts z organizmu 1 - iZDrugim
1436                                                        dMinimum = HUGE_VAL;
1437                                                        for (i = 0; i < aiRozmiarCalychGrup[1 - iZDrugim]; i++)
1438                                                        {
1439                                                                if (!(m_pMatching->IsMatched(
1440                                                                        1 - iZDrugim,
1441                                                                        aiKoniecPierwszejGrupy[1 - iZDrugim] + i)))
1442                                                                {
1443                                                                        // szukamy minimum obliczonej lokalnej odleglosci tylko wsrod
1444                                                                        // niedopasowanych jeszcze Parts
1445                                                                        if (iZDrugim == 0)
1446                                                                        {
1447                                                                                // teraz niestety musimy rozpoznac odpowiedni organizm
1448                                                                                // zeby moc indeksowac niesymetryczna tablice
1449                                                                                if (aadOdleglosciParts[iIndex[0]][i] < dMinimum)
1450                                                                                {
1451                                                                                        dMinimum = aadOdleglosciParts[iIndex[0]][i];
1452                                                                                }
1453                                                                                // przy okazji - sprawdz, czy HUGE_VAL jest rzeczywiscie max dla double
1454                                                                                assert(aadOdleglosciParts[iIndex[0]][i] < HUGE_VAL);
1455
1456                                                                        }
1457                                                                        else
1458                                                                        {
1459                                                                                if (aadOdleglosciParts[i][iIndex[1]] < dMinimum)
1460                                                                                {
1461                                                                                        dMinimum = aadOdleglosciParts[i][iIndex[1]];
1462                                                                                }
1463                                                                                // przy okazji - sprawdz, czy HUGE_VAL jest rzeczywiscie max dla double
1464                                                                                assert(aadOdleglosciParts[i][iIndex[1]] < HUGE_VAL);
1465                                                                        }
1466                                                                }
1467                                                        }
1468                                                        // sprawdz, czy minimum znaleziono - musi takie byc, bo jest cos niedopasowanego
1469                                                        assert((dMinimum >= 0.0) && (dMinimum < HUGE_VAL));
1470
1471                                                        // teraz zaznaczamy w tablicy te wszystkie Parts, ktore maja
1472                                                        // rzeczywiscie te minimalne odleglosci do Part w organizmie iZDrugim
1473                                                        for (i = 0; i < aiRozmiarCalychGrup[1 - iZDrugim]; i++)
1474                                                        {
1475                                                                if (!(m_pMatching->IsMatched(
1476                                                                        1 - iZDrugim,
1477                                                                        aiKoniecPierwszejGrupy[1 - iZDrugim] + i)))
1478                                                                {
1479                                                                        if (iZDrugim == 0)
1480                                                                        {
1481                                                                                // teraz niestety musimy rozpoznac odpowiedni organizm
1482                                                                                // zeby moc indeksowac niesymetryczna tablice
1483                                                                                if (aadOdleglosciParts[iIndex[0]][i] == dMinimum)
1484                                                                                {
1485                                                                                        // jesli w danym miejscu tablicy odleglosci jest faktyczne
1486                                                                                        // minimum odleglosci, to mamy potencjalna pare dla iIndex[1]
1487                                                                                        apvbCzyMinimalnaOdleglosc[0]->operator[](i) = true;
1488                                                                                }
1489                                                                        }
1490                                                                        else
1491                                                                        {
1492                                                                                if (aadOdleglosciParts[i][iIndex[1]] == dMinimum)
1493                                                                                {
1494                                                                                        apvbCzyMinimalnaOdleglosc[1]->operator[](i) = true;
1495                                                                                }
1496                                                                        }
1497                                                                }
1498                                                        }
1499
1500                                                        // a teraz - po znalezieniu potencjalnych elementow do dopasowania
1501                                                        // dopasowujemy pierwszy z potencjalnych
1502                                                        iDopasowywany = -1;
1503                                                        for (i = 0; i < aiRozmiarCalychGrup[1 - iZDrugim]; i++)
1504                                                        {
1505                                                                if (apvbCzyMinimalnaOdleglosc[iZDrugim]->operator[](i))
1506                                                                {
1507                                                                        iDopasowywany = i;
1508                                                                        break;
1509                                                                }
1510                                                        }
1511                                                        // musielismy znalezc jakiegos dopasowywanego
1512                                                        assert((iDopasowywany >= 0) &&
1513                                                                (iDopasowywany < aiRozmiarCalychGrup[1 - iZDrugim]));
1514
1515                                                        // no to juz mozemy dopasowac
1516                                                        m_pMatching->Match(
1517                                                                iZDrugim,
1518                                                                aiKoniecPierwszejGrupy[iZDrugim] + iIndex[iZDrugim],
1519                                                                1 - iZDrugim,
1520                                                                aiKoniecPierwszejGrupy[1 - iZDrugim] + iDopasowywany);
1521                                                        DB(printf("Przypadek 2.1.: dopasowane Parts dopasowaniebz drugim: (%2i, %2i)\n", aiKoniecPierwszejGrupy[iZDrugim] + iIndex[iZDrugim], aiKoniecPierwszejGrupy[1 - iZDrugim] + iDopasowywany);)
1522
1523                                                                //aiKoniecPierwszejGrupy[ 1-iBezDrugiego ] + iDopasowywany ;)
1524                                                                //aiKoniecPierwszejGrupy[ 1-iBezDrugiego ] + iDopasowywany ;)
1525                                                                // zmniejsz liczby niedopasowanych elementow w grupach
1526                                                                aiRozmiarGrupy[0]--;
1527                                                        aiRozmiarGrupy[1]--;
1528
1529                                                }
1530                                                else
1531                                                {
1532                                                        // jedna z grup sie juz wyczerpala
1533                                                        // wiec nie ma mozliwosci dopasowania tego drugiego Partu
1534                                                        /// i trzeba poczekac na zmiane grupy
1535                                                }
1536
1537                                                DB(printf("Przypadek 2.\n");)
1538
1539                                        }// PRZYPADEK 2.
1540                                        else
1541                                        {
1542                                                // PRZYPADEK 3. Zaden z Parts nie ma na liscie drugiego
1543                                                // AKCJA. Niezalezne dopasowanie obu Parts do pierwszych ze swojej listy
1544
1545                                                // najpierw dopasujemy do iIndex[0] w organizmie 0
1546                                                int iDopasowywany = -1;
1547                                                // poszukujemy pierwszego z listy dopasowania - w organizmie 1
1548                                                for (i = 0; i < aiRozmiarCalychGrup[1]; i++)
1549                                                {
1550                                                        if (apvbCzyMinimalnaOdleglosc[0]->operator[](i))
1551                                                        {
1552                                                                iDopasowywany = i;
1553                                                                break;
1554                                                        }
1555                                                }
1556                                                // musielismy znalezc jakiegos dopasowywanego
1557                                                assert((iDopasowywany >= 0) &&
1558                                                        (iDopasowywany < aiRozmiarCalychGrup[1]));
1559
1560                                                // teraz wlasnie dopasowujemy
1561                                                m_pMatching->Match(
1562                                                        0,
1563                                                        aiKoniecPierwszejGrupy[0] + iIndex[0],
1564                                                        1,
1565                                                        aiKoniecPierwszejGrupy[1] + iDopasowywany);
1566
1567                                                // zmniejszamy liczbe niedopasowanych Parts
1568                                                aiRozmiarGrupy[0]--;
1569                                                aiRozmiarGrupy[1]--;
1570
1571                                                // debug - dopasowanie
1572                                                DB(printf("Przypadek 3.: dopasowane Parts: (%2i, %2i)\n", aiKoniecPierwszejGrupy[0]
1573                                                        + iIndex[0], aiKoniecPierwszejGrupy[1] + iDopasowywany);)
1574
1575                                                        // teraz dopasowujemy do iIndex[1] w organizmie 1
1576                                                        iDopasowywany = -1;
1577                                                // poszukujemy pierwszego z listy dopasowania - w organizmie 0
1578                                                for (i = 0; i < aiRozmiarCalychGrup[0]; i++)
1579                                                {
1580                                                        if (apvbCzyMinimalnaOdleglosc[1]->operator[](i))
1581                                                        {
1582                                                                iDopasowywany = i;
1583                                                                break;
1584                                                        }
1585                                                }
1586                                                // musielismy znalezc jakiegos dopasowywanego
1587                                                assert((iDopasowywany >= 0) &&
1588                                                        (iDopasowywany < aiRozmiarCalychGrup[0]));
1589
1590                                                // no i teraz realizujemy dopasowanie
1591                                                m_pMatching->Match(
1592                                                        0,
1593                                                        aiKoniecPierwszejGrupy[0] + iDopasowywany,
1594                                                        1,
1595                                                        aiKoniecPierwszejGrupy[1] + iIndex[1]);
1596
1597                                                // zmniejszamy liczbe niedopasowanych Parts
1598                                                aiRozmiarGrupy[0]--;
1599                                                aiRozmiarGrupy[1]--;
1600
1601                                                // debug - dopasowanie
1602                                                DB(printf("Przypadek 3.: dopasowane Parts: (%2i, %2i)\n", aiKoniecPierwszejGrupy[0]
1603                                                        + iDopasowywany, aiKoniecPierwszejGrupy[1] + iIndex[1]);)
1604
1605
1606                                        } // PRZYPADEK 3.
1607
1608                        }// if (! bCzyKoniecGrupy)
1609                        else
1610                        {
1611                                // gdy mamy juz koniec grup - musimy zlikwidowac tablice aadOdleglosciParts
1612                                // bo za chwilke skonczy sie nam petla
1613                                for (i = 0; i < aiRozmiarCalychGrup[0]; i++)
1614                                {
1615                                        SAFEDELETEARRAY(aadOdleglosciParts[i]);
1616                                }
1617                                SAFEDELETEARRAY(aadOdleglosciParts);
1618
1619                                // musimy tez usunac tablice (wektory) mozliwosci dopasowania
1620                                SAFEDELETE(apvbCzyMinimalnaOdleglosc[0]);
1621                                SAFEDELETE(apvbCzyMinimalnaOdleglosc[1]);
1622                        }
1623                } // while (! bCzyKoniecGrupy)
1624
1625                // PO DOPASOWANIU WSZYSTKIEGO Z GRUP (CO NAJMNIEJ JEDNEJ GRUPY W CALOSCI)
1626
1627                // gdy rozmiar ktorejkolwiek z grup dopasowania spadl do zera
1628                // to musimy przesunac KoniecPierwszejGrupy (wszystkie dopasowane)
1629                for (i = 0; i < 2; i++)
1630                {
1631                        // grupy nie moga miec mniejszego niz 0 rozmiaru
1632                        assert(aiRozmiarGrupy[i] >= 0);
1633                        if (aiRozmiarGrupy[i] == 0)
1634                                aiKoniecPierwszejGrupy[i] = aiKoniecGrupyDopasowania[i];
1635                }
1636
1637                // sprawdzenie warunku konca dopasowywania - gdy nie
1638                // ma juz w jakims organizmie co dopasowywac
1639                if (aiKoniecPierwszejGrupy[0] >= m_aiPartCount[0] ||
1640                        aiKoniecPierwszejGrupy[1] >= m_aiPartCount[1])
1641                {
1642                        iCzyDopasowane = 1;
1643                }
1644        } // koniec WHILE - petli dopasowania
1645}
1646
1647/** Matches Parts in both organisms so that computation of similarity is possible.
1648        New algorithm (assures symmetry of the similarity measure) with geometry
1649        taken into consideration.
1650        Assumptions:
1651        - Models (m_Mod) are created and available.
1652        - Matching (m_pMatching) is created, but empty
1653        Exit conditions:
1654        - Matching (m_pMatching) is full
1655        @return 1 if success, 0 otherwise
1656        */
1657int ModelSimil::MatchPartsGeometry()
1658{
1659        // zaloz, ze sa modele i sa poprawne
1660        assert((m_Mod[0] != NULL) && (m_Mod[1] != NULL));
1661        assert(m_Mod[0]->isValid() && m_Mod[1]->isValid());
1662
1663        if (!CreatePartInfoTables())
1664                return 0;
1665        if (!CountPartDegrees())
1666                return 0;
1667        if (!GetPartPositions())
1668        {
1669                return 0;
1670        }
1671        if (!CountPartNeurons())
1672                return 0;
1673
1674
1675        if (m_adFactors[3] > 0)
1676        {
1677                if (!ComputePartsPositionsBySVD())
1678                {
1679                        return 0;
1680                }
1681        }
1682
1683        DB(printf("Przed sortowaniem:\n");)
1684                DB(_PrintDegrees(0);)
1685                DB(_PrintDegrees(1);)
1686
1687                if (!SortPartInfoTables())
1688                        return 0;
1689
1690        DB(printf("Po sortowaniu:\n");)
1691                DB(_PrintDegrees(0);)
1692                DB(_PrintDegrees(1);)
1693
1694                if (m_adFactors[3] > 0)
1695                {
1696                        // tutaj zacznij pętlę po przekształceniach  geometrycznych
1697                        const int NO_OF_TRANSFORM = 8; // liczba transformacji geometrycznych (na razie tylko ID i O_YZ)
1698                        // tablice transformacji współrzędnych; nie są to dokładnie tablice tranformacji, ale raczej tablice PRZEJŚĆ
1699                        // pomiędzy transformacjami;
1700                        // wartości orginalne transformacji dOrig uzyskuje się przez:
1701                        // for ( iTrans = 0; iTrans <= TRANS_INDEX; iTrans++ ) dOrig *= dMul[ iTrans ];
1702                        //const char *szTransformNames[NO_OF_TRANSFORM] = { "ID", "S_yz", "S_xz", "S_xy", "R180_z", "R180_y", "R180_z", "S_(0,0,0)" };
1703                        const int dMulX[NO_OF_TRANSFORM] = { 1, -1, -1, 1, -1, 1, -1, -1 };
1704                        const int dMulY[NO_OF_TRANSFORM] = { 1, 1, -1, -1, -1, -1, -1, 1 };
1705                        const int dMulZ[NO_OF_TRANSFORM] = { 1, 1, 1, -1, -1, -1, 1, 1 };
1706
1707#ifdef max
1708#undef max //this macro would conflict with line below
1709#endif
1710                        double dMinSimValue = std::numeric_limits<double>::max(); // minimum value of similarity
1711                        int iMinSimTransform = -1; // index of the transformation with the minimum similarity
1712                        SimilMatching *pMinSimMatching = NULL; // matching with the minimum value of similarity
1713
1714                        // remember the original positions of model 0 as computed by SVD in order to restore them later, after
1715                        // all transformations have been computed
1716                        Pt3D *StoredPositions = NULL; // array of positions of Parts, for one (0th) model
1717                        // create the stored positions
1718                        StoredPositions = new Pt3D[m_Mod[0]->getPartCount()];
1719                        assert(StoredPositions != NULL);
1720                        // copy the original positions of model 0 (store them)
1721                        int iPart; // a counter of Parts
1722                        for (iPart = 0; iPart < m_Mod[0]->getPartCount(); iPart++)
1723                        {
1724                                StoredPositions[iPart].x = m_aPositions[0][iPart].x;
1725                                StoredPositions[iPart].y = m_aPositions[0][iPart].y;
1726                                StoredPositions[iPart].z = m_aPositions[0][iPart].z;
1727                        }
1728                        // now the original positions of model 0 are stored
1729
1730
1731                        int iTransform; // a counter of geometric transformations
1732                        for (iTransform = 0; iTransform < NO_OF_TRANSFORM; iTransform++)
1733                        {
1734                                // for each geometric transformation to be done
1735                                // entry conditions:
1736                                // - models (m_Mod) exist and are available
1737                                // - matching (m_pMatching) exists and is empty
1738                                // - all properties are created and available (m_aDegrees and m_aPositions)
1739
1740                                // recompute geometric properties according to the transformation iTransform
1741                                // but only for model 0
1742                                for (iPart = 0; iPart < m_Mod[0]->getPartCount(); iPart++)
1743                                {
1744                                        // for each iPart, a part of the model iMod
1745                                        m_aPositions[0][iPart].x *= dMulX[iTransform];
1746                                        m_aPositions[0][iPart].y *= dMulY[iTransform];
1747                                        m_aPositions[0][iPart].z *= dMulZ[iTransform];
1748                                }
1749                                // now the positions are recomputed
1750                                ComputeMatching();
1751
1752                                // teraz m_pMatching istnieje i jest pełne
1753                                assert(m_pMatching != NULL);
1754                                assert(m_pMatching->IsFull() == true);
1755
1756                                // wykorzystaj to dopasowanie i geometrię do obliczenia tymczasowej wartości miary
1757                                int iTempRes = CountPartsDistance();
1758                                // załóż sukces
1759                                assert(iTempRes == 1);
1760                                double dCurrentSim = m_adFactors[0] * double(m_iDV) +
1761                                        m_adFactors[1] * double(m_iDD) +
1762                                        m_adFactors[2] * double(m_iDN) +
1763                                        m_adFactors[3] * double(m_dDG);
1764                                // załóż poprawną wartość podobieństwa
1765                                assert(dCurrentSim >= 0.0);
1766
1767                                // porównaj wartość obliczoną z dotychczasowym minimum
1768                                if (dCurrentSim < dMinSimValue)
1769                                {
1770                                        // jeśli uzyskano mniejszą wartość dopasowania,
1771                                        // to zapamiętaj to przekształcenie geometryczne i dopasowanie
1772                                        dMinSimValue = dCurrentSim;
1773                                        iMinSimTransform = iTransform;
1774                                        SAFEDELETE(pMinSimMatching);
1775                                        pMinSimMatching = new SimilMatching(*m_pMatching);
1776                                        assert(pMinSimMatching != NULL);
1777                                }
1778
1779                                // teraz już można usunąć stare dopasowanie (dla potrzeb następnego przebiegu pętli)
1780                                m_pMatching->Empty();
1781                        } // for ( iTransform )
1782
1783                        // po pętli przywróć najlepsze dopasowanie
1784                        delete m_pMatching;
1785                        m_pMatching = pMinSimMatching;
1786
1787                        DB(printf("Matching has been chosen!\n");)
1788                                // print the name of the chosen transformation:
1789                                // printf("Chosen transformation: %s\n", szTransformNames[ iMinSimTransform ] );
1790
1791                                // i przywróć jednocześnie pozycje Parts modelu 0 dla tego dopasowania
1792                                // - najpierw przywroc Parts pozycje orginalne, po SVD
1793                                for (iPart = 0; iPart < m_Mod[0]->getPartCount(); iPart++)
1794                                {
1795                                        m_aPositions[0][iPart].x = StoredPositions[iPart].x;
1796                                        m_aPositions[0][iPart].y = StoredPositions[iPart].y;
1797                                        m_aPositions[0][iPart].z = StoredPositions[iPart].z;
1798                                }
1799                        // - usun teraz zapamietane pozycje Parts
1800                        delete[] StoredPositions;
1801                        // - a teraz oblicz na nowo wspolrzedne wlasciwego przeksztalcenia dla model 0
1802                        for (iTransform = 0; iTransform <= iMinSimTransform; iTransform++)
1803                        {
1804                                // for each transformation before and INCLUDING iMinTransform
1805                                // do the transformation (only model 0)
1806                                for (iPart = 0; iPart < m_Mod[0]->getPartCount(); iPart++)
1807                                {
1808                                        m_aPositions[0][iPart].x *= dMulX[iTransform];
1809                                        m_aPositions[0][iPart].y *= dMulY[iTransform];
1810                                        m_aPositions[0][iPart].z *= dMulZ[iTransform];
1811                                }
1812                        }
1813
1814                }
1815                else
1816                {
1817                        ComputeMatching();
1818                }
1819        // teraz dopasowanie musi byc pelne - co najmniej w jednym organizmie musza byc
1820        // wszystkie elementy dopasowane
1821        assert(m_pMatching->IsFull() == true);
1822
1823        //    _PrintDegrees(0);
1824        //    _PrintDegrees(1);
1825
1826        DB(_PrintPartsMatching();)
1827
1828                return 1;
1829}
1830
1831void ModelSimil::_PrintSeamnessTable(std::vector<int> *pTable, int iCount)
1832{
1833        int i;
1834        printf("      ");
1835        for (i = 0; i < iCount; i++)
1836                printf("%3i ", i);
1837        printf("\n      ");
1838        for (i = 0; i < iCount; i++)
1839        {
1840
1841                printf("%3i ", pTable->operator[](i));
1842        }
1843        printf("\n");
1844}
1845
1846/** Computes elements of similarity (m_iDD, m_iDN, m_dDG) based on underlying matching.
1847        Assumptions:
1848        - Matching (m_pMatching) exists and is full.
1849        - Internal arrays m_aDegrees and m_aPositions exist and are properly filled in
1850        Exit conditions:
1851        - Elements of similarity are computed (m)iDD, m_iDN, m_dDG).
1852        @return 1 if success, otherwise 0.
1853        */
1854int ModelSimil::CountPartsDistance()
1855{
1856        int i, temp;
1857
1858        // assume existence of m_pMatching
1859        assert(m_pMatching != NULL);
1860        // musi byc pelne!
1861        assert(m_pMatching->IsFull() == true);
1862
1863        // roznica w stopniach
1864        m_iDD = 0;
1865        // roznica w liczbie neuronów
1866        m_iDN = 0;
1867        // roznica w odleglosci dopasowanych Parts
1868        m_dDG = 0.0;
1869
1870        int iOrgPart, iOrgMatchedPart; // orginalny indeks Part i jego dopasowanego Part
1871        int iMatchedPart; // indeks (wg sortowania) dopasowanego Part
1872
1873        // wykorzystanie dopasowania do zliczenia m_iDD - roznicy w stopniach
1874        // i m_iDN - roznicy w liczbie neuronow
1875        // petla w wiekszej tablicy
1876        for (i = 0; i < m_aiPartCount[1 - m_iSmaller]; i++)
1877        {
1878                // dla kazdego elementu [i] z wiekszego organizmu
1879                // pobierz jego orginalny indeks w organizmie z tablicy TDN
1880                iOrgPart = m_aDegrees[1 - m_iSmaller][i][0];
1881                if (!(m_pMatching->IsMatched(1 - m_iSmaller, i)))
1882                {
1883                        // gdy nie zostal dopasowany
1884                        // dodaj jego stopien do DD
1885                        m_iDD += m_aDegrees[1 - m_iSmaller][i][1];
1886                        // dodaj liczbe neuronow do DN
1887                        m_iDN += m_aDegrees[1 - m_iSmaller][i][3];
1888                        // i oblicz odleglosc tego Part od srodka organizmu (0,0,0)
1889                        // (uzyj orginalnego indeksu)
1890                        //no need to compute distane when m_dDG weight is 0
1891                        m_dDG += m_adFactors[3] == 0 ? 0 : m_aPositions[1 - m_iSmaller][iOrgPart].length();
1892                }
1893                else
1894                {
1895                        // gdy byl dopasowany
1896                        // pobierz indeks (po sortowaniu) tego dopasowanego Part
1897                        iMatchedPart = m_pMatching->GetMatchedIndex(1 - m_iSmaller, i);
1898                        // pobierz indeks orginalny tego dopasowanego Part
1899                        iOrgMatchedPart = m_aDegrees[m_iSmaller][iMatchedPart][0];
1900                        // dodaj ABS roznicy stopni do DD
1901                        temp = m_aDegrees[1 - m_iSmaller][i][1] -
1902                                m_aDegrees[m_iSmaller][iMatchedPart][1];
1903                        m_iDD += abs(temp);
1904                        // dodaj ABS roznicy neuronow do DN
1905                        temp = m_aDegrees[1 - m_iSmaller][i][3] -
1906                                m_aDegrees[m_iSmaller][iMatchedPart][3];
1907                        m_iDN += abs(temp);
1908                        // pobierz polozenie dopasowanego Part
1909                        Pt3D MatchedPartPos(m_aPositions[m_iSmaller][iOrgMatchedPart]);
1910                        // dodaj euklidesowa odleglosc Parts do sumy odleglosci
1911                        //no need to compute distane when m_dDG weight is 0
1912                        m_dDG += m_adFactors[3] == 0 ? 0 : m_aPositions[1 - m_iSmaller][iOrgPart].distanceTo(MatchedPartPos);
1913                }
1914        }
1915
1916        // obliczenie i dodanie różnicy w liczbie neuronów OnJoint...
1917        temp = m_aOnJoint[0][3] - m_aOnJoint[1][3];
1918        m_iDN += abs(temp);
1919        DB(printf("OnJoint DN: %i\n", abs(temp));)
1920                // ... i Anywhere
1921                temp = m_aAnywhere[0][3] - m_aAnywhere[1][3];
1922        m_iDN += abs(temp);
1923        DB(printf("Anywhere DN: %i\n", abs(temp));)
1924
1925                return 1;
1926}
1927
1928/** Computes new positions of Parts of both of organisms stored in the object.
1929                Assumptions:
1930                - models (m_Mod) are created and valid
1931                - positions (m_aPositions) are created and filled with original positions of Parts
1932                - the sorting algorithm was not yet run on the array m_aDegrees
1933                @return true if successful; false otherwise
1934                */
1935bool ModelSimil::ComputePartsPositionsBySVD()
1936{
1937        bool bResult = true; // the result; assume a success
1938
1939        // check assumptions
1940        // the models
1941        assert(m_Mod[0] != NULL && m_Mod[0]->isValid());
1942        assert(m_Mod[1] != NULL && m_Mod[1]->isValid());
1943        // the position arrays
1944        assert(m_aPositions[0] != NULL);
1945        assert(m_aPositions[1] != NULL);
1946
1947        int iMod; // a counter of models
1948        // use SVD to obtain different point of view on organisms
1949        // and store the new positions (currently the original ones are still stored)
1950        for (iMod = 0; iMod < 2; iMod++)
1951        {
1952                // prepare the vector of errors of approximation for the SVD
1953                std::vector<double> vEigenvalues;
1954                int nSize = m_Mod[iMod]->getPartCount();
1955
1956                double *pDistances = new double[nSize * nSize];
1957
1958                for (int i = 0; i < nSize; i++)
1959                {
1960                        pDistances[i] = 0;
1961                }
1962
1963                Model *pModel = m_Mod[iMod]; // use the model of the iMod (current) organism
1964                int iP1, iP2; // indices of Parts in the model
1965                Part *P1, *P2; // pointers to Parts
1966                Pt3D P1Pos, P2Pos; // positions of parts
1967                double dDistance; // the distance between Parts
1968
1969                double *weights = new double[nSize];
1970                for (int i = 0; i < nSize; i++)
1971                {
1972                        if (wMDS == 1)
1973                                weights[i] = 0;
1974                        else
1975                                weights[i] = 1;
1976                }
1977
1978                if (wMDS == 1)
1979                        for (int i = 0; i < pModel->getJointCount(); i++)
1980                        {
1981                                weights[pModel->getJoint(i)->p1_refno]++;
1982                                weights[pModel->getJoint(i)->p2_refno]++;
1983                        }
1984
1985                for (iP1 = 0; iP1 < pModel->getPartCount(); iP1++)
1986                {
1987                        // for each iP1, a Part index in the model of organism iMod
1988                        P1 = pModel->getPart(iP1);
1989                        // get the position of the Part
1990                        P1Pos = P1->p;
1991                        if (fixedZaxis == 1)
1992                        {
1993                                P1Pos.z = 0; //fixed vertical axis, so pretend all points are on the xy plane
1994                        }
1995                        for (iP2 = 0; iP2 < pModel->getPartCount(); iP2++)
1996                        {
1997                                // for each (iP1, iP2), a pair of Parts index in the model
1998                                P2 = pModel->getPart(iP2);
1999                                // get the position of the Part
2000                                P2Pos = P2->p;
2001                                if (fixedZaxis == 1)
2002                                {
2003                                        P2Pos.z = 0; //fixed vertical axis, so pretend all points are on the xy plane
2004                                }
2005                                // compute the geometric (Euclidean) distance between the Parts
2006                                dDistance = P1Pos.distanceTo(P2Pos);
2007                                // store the distance
2008                                pDistances[iP1 * nSize + iP2] = dDistance;
2009                        } // for (iP2)
2010                } // for (iP1)
2011
2012                MatrixTools::weightedMDS(vEigenvalues, nSize, pDistances, m_aPositions[iMod], weights);
2013                if (fixedZaxis == 1) //restore the original vertical coordinate of each Part
2014                {
2015                        for (int part = 0; part < pModel->getPartCount(); part++)
2016                        {
2017                                m_aPositions[iMod][part].z = pModel->getPart(part)->p.z;
2018                        }
2019                }
2020
2021                delete[] pDistances;
2022                delete[] weights;
2023        }
2024
2025        return bResult;
2026}
2027
2028/**     Evaluates distance between two given genotypes. The distance depends strongly
2029        on weights set and the matching algorithm used.
2030        @param G0 Pointer to the first of compared genotypes
2031        @param G1 Pointer to the second of compared genotypes.
2032        @return Distance between two genotypes.
2033        @sa m_adFactors, matching_method
2034        */
2035double ModelSimil::EvaluateDistance(const Geno *G0, const Geno *G1)
2036{
2037        return matching_method == 0 ? EvaluateDistanceHungarian(G0, G1) : EvaluateDistanceGreedy(G0, G1);
2038}
2039
2040void ModelSimil::p_evaldistance(ExtValue *args, ExtValue *ret)
2041{
2042        Geno *g1 = GenoObj::fromObject(args[1]);
2043        Geno *g2 = GenoObj::fromObject(args[0]);
2044        if ((!g1) || (!g2))
2045                ret->setEmpty();
2046        else
2047                ret->setDouble(EvaluateDistance(g1, g2));
2048}
2049
2050void ModelSimil::FillPartsDistances(double*& dist, int bigger, int smaller, bool geo)
2051{
2052        for (int i = 0; i < bigger; i++)
2053        {
2054                for (int j = 0; j < bigger; j++)
2055                {
2056                        // assign penalty for unassignment for vertex from bigger model
2057                        if (j >= smaller)
2058                        {
2059                                if (geo)
2060                                        dist[i*bigger + j] += m_adFactors[3] * m_aPositions[1 - m_iSmaller][i].length();
2061                                else
2062                                        dist[i*bigger + j] = m_adFactors[1] * m_aDegrees[1 - m_iSmaller][i][DEGREE] + m_adFactors[2] * m_aDegrees[1 - m_iSmaller][i][NEURONS];
2063                        }
2064                        // compute distance between parts
2065                        else
2066                        {
2067                                if (geo)
2068                                        dist[i*bigger + j] += m_adFactors[3] * m_aPositions[1 - m_iSmaller][i].distanceTo(m_aPositions[m_iSmaller][j]);
2069                                else
2070                                        dist[i*bigger + j] = m_adFactors[1] * abs(m_aDegrees[1 - m_iSmaller][i][DEGREE] - m_aDegrees[m_iSmaller][j][DEGREE])
2071                                        + m_adFactors[2] * abs(m_aDegrees[1 - m_iSmaller][i][NEURONS] - m_aDegrees[m_iSmaller][j][NEURONS]);
2072                        }
2073
2074                }
2075        }
2076}
2077
2078double ModelSimil::EvaluateDistanceHungarian(const Geno *G0, const Geno *G1)
2079{
2080        double dResult = 0.0;
2081
2082        m_Gen[0] = G0;
2083        m_Gen[1] = G1;
2084
2085        // create models of objects to compare
2086        m_Mod[0] = newModel(m_Gen[0]);
2087        m_Mod[1] = newModel(m_Gen[1]);
2088
2089        if (m_Mod[0] == NULL || m_Mod[1] == NULL)
2090                return 0.0;
2091
2092        //Get information about vertex degrees, neurons and 3D location
2093        if (!CreatePartInfoTables())
2094                return 0.0;
2095        if (!CountPartDegrees())
2096                return 0.0;
2097        if (!GetPartPositions())
2098                return 0.0;
2099        if (!CountPartNeurons())
2100                return 0.0;
2101
2102        m_iSmaller = m_Mod[0]->getPartCount() <= m_Mod[1]->getPartCount() ? 0 : 1;
2103        int nSmaller = m_Mod[m_iSmaller]->getPartCount();
2104        int nBigger = m_Mod[1 - m_iSmaller]->getPartCount();
2105
2106        double* partsDistances = new double[nBigger*nBigger]();
2107        FillPartsDistances(partsDistances, nBigger, nSmaller, false);
2108        int *assignment = new int[nBigger]();
2109
2110        HungarianAlgorithm hungarian;
2111
2112        if (m_adFactors[3] > 0)
2113        {
2114                if (!ComputePartsPositionsBySVD())
2115                {
2116                        return 0.0;
2117                }
2118
2119                // tutaj zacznij pętlę po przekształceniach  geometrycznych
2120                const int NO_OF_TRANSFORM = 8; // liczba transformacji geometrycznych (na razie tylko ID i O_YZ)
2121                // tablice transformacji współrzędnych; nie są to dokładnie tablice tranformacji, ale raczej tablice PRZEJŚĆ
2122                // pomiędzy transformacjami;
2123                const int dMulX[NO_OF_TRANSFORM] = { 1, -1, -1, 1, -1, 1, -1, -1 };
2124                const int dMulY[NO_OF_TRANSFORM] = { 1, 1, -1, -1, -1, -1, -1, 1 };
2125                const int dMulZ[NO_OF_TRANSFORM] = { 1, 1, 1, -1, -1, -1, 1, 1 };
2126
2127                std::vector<int> minAssignment(nBigger);
2128#ifdef max
2129#undef max //this macro would conflict with line below
2130#endif
2131                double dMinSimValue = std::numeric_limits<double>::max(); // minimum value of similarity
2132
2133                int iTransform; // a counter of geometric transformations
2134                for (iTransform = 0; iTransform < NO_OF_TRANSFORM; iTransform++)
2135                {
2136                        // for each geometric transformation to be done
2137                        // entry conditions:
2138                        // - models (m_Mod) exist and are available
2139                        // - all properties are created and available (m_aDegrees and m_aPositions)
2140                        double* tmpPartsDistances = new double[nBigger*nBigger]();
2141                        std::copy(partsDistances, partsDistances + nBigger * nBigger, tmpPartsDistances);
2142                        // recompute geometric properties according to the transformation iTransform
2143                        // but only for model 0
2144                        for (int iPart = 0; iPart < m_Mod[m_iSmaller]->getPartCount(); iPart++)
2145                        {
2146                                // for each iPart, a part of the model iMod
2147                                m_aPositions[m_iSmaller][iPart].x *= dMulX[iTransform];
2148                                m_aPositions[m_iSmaller][iPart].y *= dMulY[iTransform];
2149                                m_aPositions[m_iSmaller][iPart].z *= dMulZ[iTransform];
2150                        }
2151                        // now the positions are recomputed
2152
2153                        FillPartsDistances(tmpPartsDistances, nBigger, nSmaller, true);
2154                        std::fill_n(assignment, nBigger, 0);
2155                        double dCurrentSim = hungarian.Solve(tmpPartsDistances, assignment, nBigger, nBigger);
2156
2157                        delete[] tmpPartsDistances;
2158                        // załóż poprawną wartość podobieństwa
2159                        assert(dCurrentSim >= 0.0);
2160
2161                        // porównaj wartość obliczoną z dotychczasowym minimum
2162                        if (dCurrentSim < dMinSimValue)
2163                        {
2164                                dMinSimValue = dCurrentSim;
2165                                if (saveMatching)
2166                                {
2167                                        minAssignment.clear();
2168                                        minAssignment.insert(minAssignment.begin(), assignment, assignment + nBigger);
2169                                }
2170                        }
2171                }
2172
2173                dResult = dMinSimValue;
2174                if (saveMatching)
2175                        std::copy(minAssignment.begin(), minAssignment.end(), assignment);
2176        }
2177        else
2178        {
2179                dResult = hungarian.Solve(partsDistances, assignment, nBigger, nBigger);
2180        }
2181
2182        //add difference in anywhere and onJoint neurons
2183        dResult += m_adFactors[2] * (abs(m_aOnJoint[0][3] - m_aOnJoint[1][3]) + abs(m_aAnywhere[0][3] - m_aAnywhere[1][3]));
2184        //add difference in part numbers
2185        dResult += (nBigger - nSmaller) * m_adFactors[0];
2186
2187        // delete degree arrays created in CreatePartInfoTables
2188        SAFEDELETEARRAY(m_aDegrees[0]);
2189        SAFEDELETEARRAY(m_aDegrees[1]);
2190
2191        // and position arrays
2192        SAFEDELETEARRAY(m_aPositions[0]);
2193        SAFEDELETEARRAY(m_aPositions[1]);
2194
2195        // delete created models
2196        SAFEDELETE(m_Mod[0]);
2197        SAFEDELETE(m_Mod[1]);
2198
2199        delete[] assignment;
2200        delete[] partsDistances;
2201
2202        return dResult;
2203}
Note: See TracBrowser for help on using the repository browser.