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

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

Fixed error messages, changed int to paInt (for 32/64-bit safety), changed two int variables to bool as they should be

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