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

Last change on this file since 666 was 666, checked in by Maciej Komosinski, 15 months ago

Added missing #include

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