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

Last change on this file since 601 was 601, checked in by oriona, 4 years ago

Possibility of fixing z axis during the alignment added.

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