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

Last change on this file since 854 was 818, checked in by oriona, 6 years ago

Help to wMDS field added.

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