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

Last change on this file since 451 was 451, checked in by oriona, 5 years ago

Unused INit commented out. Command line information updated.

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