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

Last change on this file since 870 was 869, checked in by Maciej Komosinski, 5 years ago

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

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