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

Last change on this file since 370 was 362, checked in by Maciej Komosinski, 9 years ago

Restored object descriptions/docs that were lost

  • Property svn:eol-style set to native
File size: 77.0 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, iNIt; // ilościowe określenie różnic stopnia, liczby neuronów i połączeń Parts
1138        double dGeo; // ilościowe określenie różnic geometrycznych pomiędzy Parts
1139        // indeksy konkretnych Parts - indeksy sa ZERO-BASED, choć właściwy dostep
1140        // do informacji o Part wymaga dodania aiKoniecPierwszejGrupy[]
1141        // tylko aadOdleglosciParts[][] indeksuje sie bezposrednio zawartoscia iIndex[]
1142        int iIndex[2];
1143        int iPartIndex[ 2 ] = {-1, -1}; // at [iModel]: original index of a Part for the specified model (iModel)
1144
1145        // debug - wypisz zakres dopasowywanych indeksow
1146        DB(printf("Organizm 0: grupa: (%2i, %2i)\n", aiKoniecPierwszejGrupy[0],
1147                  aiKoniecGrupyDopasowania[0]);)
1148                DB(printf("Organizm 1: grupa: (%2i, %2i)\n", aiKoniecPierwszejGrupy[1],
1149                          aiKoniecGrupyDopasowania[1]);)
1150
1151        for (i = 0; i < aiRozmiarCalychGrup[0]; i++)
1152        {
1153
1154            // iterujemy i - Parts organizmu 0
1155            // (indeks podstawowy - aiKoniecPierwszejGrupy[0])
1156
1157            if (!(m_pMatching->IsMatched(0, aiKoniecPierwszejGrupy[0] + i)))
1158            {
1159                // interesuja nas tylko te niedopasowane jeszcze (i)
1160                for (j = 0; j < aiRozmiarCalychGrup[1]; j++)
1161                {
1162
1163                    // iterujemy j - Parts organizmu 1
1164                    // (indeks podstawowy - aiKoniecPierwszejGrupy[1])
1165
1166                    if (!(m_pMatching->IsMatched(1, aiKoniecPierwszejGrupy[1] + j)))
1167                    {
1168                        // interesuja nas tylko te niedopasowane jeszcze (j)
1169                        // teraz obliczymy lokalne różnice pomiędzy Parts
1170                        iDeg = abs(m_aDegrees[0][ aiKoniecPierwszejGrupy[0] + i ][1]
1171                                   - m_aDegrees[1][ aiKoniecPierwszejGrupy[1] + j ][1]);
1172
1173                        iNIt = abs(m_aDegrees[0][ aiKoniecPierwszejGrupy[0] + i ][2]
1174                                   - m_aDegrees[1][ aiKoniecPierwszejGrupy[1] + j ][2]);
1175
1176                        iNeu = abs(m_aDegrees[0][ aiKoniecPierwszejGrupy[0] + i ][3]
1177                                   - m_aDegrees[1][ aiKoniecPierwszejGrupy[1] + j ][3]);
1178
1179                        // obliczenie także lokalnych różnic geometrycznych pomiędzy Parts
1180                        // find original indices of Parts for both of the models
1181                        iPartIndex[ 0 ] = m_aDegrees[0][ aiKoniecPierwszejGrupy[0] + i ][ 0 ];
1182                        iPartIndex[ 1 ] = m_aDegrees[1][ aiKoniecPierwszejGrupy[1] + j ][ 0 ];
1183                        // now compute the geometrical distance of these Parts (use m_aPositions
1184                        // which should be computed by SVD)
1185                        Pt3D Part0Pos(m_aPositions[ 0 ][ iPartIndex[ 0 ] ]);
1186                        Pt3D Part1Pos(m_aPositions[ 1 ][ iPartIndex[ 1 ] ]);
1187                        dGeo = m_adFactors[3] == 0 ? 0 : Part0Pos.distanceTo(Part1Pos); //no need to compute distane when m_dDG weight is 0
1188
1189                        // tutaj prawdopodobnie należy jeszcze dodać sprawdzanie
1190                        // identyczności pozostałych własności (oprócz geometrii)
1191                        // - żeby móc stwierdzić w ogóle identyczność Parts
1192
1193                        // i ostateczna odleglosc indukowana przez te roznice
1194                        // (tutaj nie ma różnicy w liczbie wszystkich wierzchołków)
1195                        aadOdleglosciParts[i][j] = m_adFactors[1] * double(iDeg) +
1196                                m_adFactors[2] * double(iNeu) +
1197                                m_adFactors[3] * dGeo;
1198                        DB(printf("Parts Distance (%2i,%2i) = %.3lf\n", aiKoniecPierwszejGrupy[0] + i,
1199                                  aiKoniecPierwszejGrupy[1] + j, aadOdleglosciParts[i][j]);)
1200                                DB(printf("Parts geometrical distance = %.20lf\n", dGeo);)
1201                                DB(printf("Pos0: (%.3lf %.3lf %.3lf)\n", Part0Pos.x, Part0Pos.y, Part0Pos.z);)
1202                                DB(printf("Pos1: (%.3lf %.3lf %.3lf)\n", Part1Pos.x, Part1Pos.y, Part1Pos.z);)
1203                    }
1204                }
1205            }
1206        }
1207
1208        // tutaj - sprawdzic tylko, czy tablica odleglosci lokalnych jest poprawnie obliczona
1209
1210        // WYKORZYSTANIE TABLICY ODLEGLOSCI DO BUDOWY DOPASOWANIA
1211
1212        // trzeba raczej iterować aż do zebrania wszystkich możliwych dopasowań w grupie
1213        // dlatego wprowadzamy dodatkowa zmienna - czy skonczyla sie juz grupa
1214        bool bCzyKoniecGrupy = false;
1215        while (!bCzyKoniecGrupy)
1216        {
1217            for (i = 0; i < 2; i++)
1218            {
1219                // iterujemy (i) po organizmach
1220                // szukamy najpierw jakiegoś niedopasowanego jeszcze Part w organizmach
1221
1222                // zakładamy, że nie ma takiego Part
1223                iIndex[i] = -1;
1224
1225                for (j = 0; j < aiRozmiarCalychGrup[ i ]; j++)
1226                {
1227                    // iterujemy (j) - Parts organizmu (i)
1228                    // (indeks podstawowy - aiKoniecPierwszejGrupy[0])
1229                    if (!(m_pMatching->IsMatched(i, aiKoniecPierwszejGrupy[i] + j)))
1230                    {
1231                        // gdy mamy w tej grupie jakis niedopasowany element, to ustawiamy
1232                        // iIndex[i] (chcemy w zasadzie pierwszy taki)
1233                        iIndex[i] = j;
1234                        break;
1235                    }
1236                }
1237
1238                // sprawdzamy, czy w ogole znaleziono taki Part
1239                if (iIndex[i] < 0)
1240                {
1241                    // gdy nie znaleziono takiego Part - mamy koniec dopasowywania w
1242                    // tych grupach
1243                    bCzyKoniecGrupy = true;
1244                }
1245                // sprawdz poprawnosc indeksu niedopasowanego Part - musi byc w aktualnej grupie
1246                assert((iIndex[i] >= -1) && (iIndex[i] < aiRozmiarCalychGrup[i]));
1247            }
1248
1249
1250            // teraz mamy sytuacje:
1251            // - mamy w iIndex[0] i iIndex[1] indeksy pierwszych niedopasowanych Part
1252            //          w organizmach, albo
1253            // - nie ma w ogóle już czego dopasowywać (należy przejść do innej grupy)
1254            if (!bCzyKoniecGrupy)
1255            {
1256                // gdy nie ma jeszcze końca żadnej z grup - możemy dopasowywać
1257                // najpierw wyszukujemy w tablicy minimum odległości od tych
1258                // wyznaczonych Parts
1259
1260                // najpierw wyczyscimy wektory potencjalnych dopasowan
1261                // dla organizmu 1 (o rozmiarze grupy z 0)
1262                for (i = 0; i < aiRozmiarCalychGrup[0]; i++)
1263                {
1264                    apvbCzyMinimalnaOdleglosc[1]->operator[](i) = false;
1265                }
1266                // dla organizmu 0 (o rozmiarze grup z 1)
1267                for (i = 0; i < aiRozmiarCalychGrup[1]; i++)
1268                {
1269                    apvbCzyMinimalnaOdleglosc[0]->operator[](i) = false;
1270                }
1271
1272                // szukanie minimum dla Part o indeksie iIndex[0] w organizmie 0
1273                // wsrod niedopasowanych Parts z organizmu 1
1274                // zakładamy, że nie znaleŸliśmy jeszcze minimum
1275                double dMinimum = HUGE_VAL;
1276                for (i = 0; i < aiRozmiarCalychGrup[1]; i++)
1277                {
1278                    if (!(m_pMatching->IsMatched(1, aiKoniecPierwszejGrupy[1] + i)))
1279                    {
1280
1281                        // szukamy minimum obliczonej lokalnej odleglosci tylko wsrod
1282                        // niedopasowanych jeszcze Parts
1283                        if (aadOdleglosciParts[ iIndex[0] ][ i ] < dMinimum)
1284                        {
1285                            dMinimum = aadOdleglosciParts[ iIndex[0] ][ i ];
1286                        }
1287
1288                        // przy okazji - sprawdz, czy HUGE_VAL jest rzeczywiscie max dla double
1289                        assert(aadOdleglosciParts[ iIndex[0] ][ i ] < HUGE_VAL);
1290                    }
1291                }
1292                // sprawdz, czy minimum znaleziono - musi takie byc, bo jest cos niedopasowanego
1293                assert((dMinimum >= 0.0) && (dMinimum < HUGE_VAL));
1294
1295                // teraz zaznaczamy w tablicy te wszystkie Parts, ktore maja
1296                // rzeczywiscie te minimalna odleglosc do Part iIndex[0] w organizmie 0
1297                for (i = 0; i < aiRozmiarCalychGrup[1]; i++)
1298                {
1299                    if (!(m_pMatching->IsMatched(1, aiKoniecPierwszejGrupy[1] + i)))
1300                    {
1301                        if (aadOdleglosciParts[ iIndex[0] ][ i ] == dMinimum)
1302                        {
1303                            // jesli w danym miejscu tablicy odleglosci jest faktyczne
1304                            // minimum odleglosci, to mamy potencjalna pare dla iIndex[0]
1305                            apvbCzyMinimalnaOdleglosc[0]->operator[](i) = true;
1306                        }
1307
1308                        // sprawdz poprawnosc znalezionego wczesniej minimum
1309                        assert(aadOdleglosciParts[ iIndex[0] ][ i ] >= dMinimum);
1310                    }
1311                }
1312
1313                // podobnie szukamy minimum dla Part o indeksie iIndex[1] w organizmie 1
1314                // wsrod niedopasowanych Parts z organizmu 0
1315                dMinimum = HUGE_VAL;
1316                for (i = 0; i < aiRozmiarCalychGrup[0]; i++)
1317                {
1318                    if (!(m_pMatching->IsMatched(0, aiKoniecPierwszejGrupy[0] + i)))
1319                    {
1320                        // szukamy minimum obliczonej lokalnej odleglosci tylko wsrod
1321                        // niedopasowanych jeszcze Parts
1322                        if (aadOdleglosciParts[ i ][ iIndex[1] ] < dMinimum)
1323                        {
1324                            dMinimum = aadOdleglosciParts[ i ][ iIndex[1] ];
1325                        }
1326                        // przy okazji - sprawdz, czy HUGE_VAL jest rzeczywiscie max dla double
1327                        assert(aadOdleglosciParts[ i ][ iIndex[1] ] < HUGE_VAL);
1328                    }
1329                }
1330                // sprawdz, czy minimum znaleziono - musi takie byc, bo jest cos niedopasowanego
1331                assert((dMinimum >= 0.0) && (dMinimum < HUGE_VAL));
1332
1333                // teraz zaznaczamy w tablicy te wszystkie Parts, ktore maja
1334                // rzeczywiscie te minimalne odleglosci do Part iIndex[1] w organizmie 1
1335                for (i = 0; i < aiRozmiarCalychGrup[0]; i++)
1336                {
1337                    if (!(m_pMatching->IsMatched(0, aiKoniecPierwszejGrupy[0] + i)))
1338                    {
1339                        if (aadOdleglosciParts[ i ][ iIndex[1] ] == dMinimum)
1340                        {
1341                            // jesli w danym miejscu tablicy odleglosci jest faktyczne
1342                            // minimum odleglosci, to mamy potencjalna pare dla iIndex[1]
1343                            apvbCzyMinimalnaOdleglosc[1]->operator[](i) = true;
1344                        }
1345
1346                        // sprawdz poprawnosc znalezionego wczesniej minimum
1347                        assert(aadOdleglosciParts[ i ][ iIndex[1] ] >= dMinimum);
1348                    }
1349                }
1350
1351                // teraz mamy juz poszukane potencjalne grupy dopasowania - musimy
1352                // zdecydowac, co do czego dopasujemy!
1353                // szukamy Part iIndex[0] posrod mozliwych do dopasowania dla Part iIndex[1]
1354                // szukamy takze Part iIndex[1] posrod mozliwych do dopasowania dla Part iIndex[0]
1355                bool PartZ1NaLiscie0 = apvbCzyMinimalnaOdleglosc[0]->operator[](iIndex[1]);
1356                bool PartZ0NaLiscie1 = apvbCzyMinimalnaOdleglosc[1]->operator[](iIndex[0]);
1357
1358                if (PartZ1NaLiscie0 && PartZ0NaLiscie1)
1359                {
1360                    // PRZYPADEK 1. Oba Parts maja sie wzajemnie na listach mozliwych
1361                    // dopasowan.
1362                    // AKCJA. Dopasowanie wzajemne do siebie.
1363
1364                    m_pMatching->Match(0, aiKoniecPierwszejGrupy[0] + iIndex[0],
1365                                       1, aiKoniecPierwszejGrupy[1] + iIndex[1]);
1366
1367                    // zmniejsz liczby niedopasowanych elementow w grupach
1368                    aiRozmiarGrupy[0]--;
1369                    aiRozmiarGrupy[1]--;
1370                    // debug - co zostalo dopasowane
1371                    DB(printf("Przypadek 1.: dopasowane Parts: (%2i, %2i)\n", aiKoniecPierwszejGrupy[0]
1372                              + iIndex[0], aiKoniecPierwszejGrupy[1] + iIndex[1]);)
1373
1374                }// PRZYPADEK 1.
1375                else
1376                    if (PartZ1NaLiscie0 || PartZ0NaLiscie1)
1377                {
1378                    // PRZYPADEK 2. Tylko jeden z Parts ma drugiego na swojej liscie
1379                    // mozliwych dopasowan
1380                    // AKCJA. Dopasowanie jednego jest proste (tego, ktory nie ma
1381                    // na swojej liscie drugiego). Dla tego drugiego nalezy powtorzyc
1382                    // duza czesc obliczen (znalezc mu nowa mozliwa pare)
1383
1384                    // indeks organizmu, ktorego Part nie ma dopasowywanego Part
1385                    // z drugiego organizmu na swojej liscie
1386                    int iBezDrugiego;
1387
1388                    // okreslenie indeksu organizmu z dopasowywanym Part
1389                    if (!PartZ1NaLiscie0)
1390                    {
1391                        iBezDrugiego = 0;
1392                    }
1393                    else
1394                    {
1395                        iBezDrugiego = 1;
1396                    }
1397                    // sprawdz indeks organizmu
1398                    assert((iBezDrugiego == 0) || (iBezDrugiego == 1));
1399
1400                    int iDopasowywany = -1;
1401                    // poszukujemy pierwszego z listy dopasowania
1402                    for (i = 0; i < aiRozmiarCalychGrup[ 1 - iBezDrugiego ]; i++)
1403                    {
1404                        if (apvbCzyMinimalnaOdleglosc[iBezDrugiego]->operator[](i))
1405                        {
1406                            iDopasowywany = i;
1407                            break;
1408                        }
1409                    }
1410                    // sprawdz poprawnosc indeksu dopasowywanego (musimy cos znalezc!)
1411                    // nieujemny i w odpowiedniej grupie!
1412                    assert((iDopasowywany >= 0) &&
1413                           (iDopasowywany < aiRozmiarCalychGrup[ 1 - iBezDrugiego ]));
1414
1415                    // znalezlismy 1. Part z listy dopasowania - dopasowujemy!
1416                    m_pMatching->Match(
1417                                       iBezDrugiego,
1418                                       aiKoniecPierwszejGrupy[ iBezDrugiego ] + iIndex[ iBezDrugiego ],
1419                                       1 - iBezDrugiego,
1420                                       aiKoniecPierwszejGrupy[ 1 - iBezDrugiego ] + iDopasowywany);
1421                    DB(printf("Przypadek 2.1.: dopasowane Parts dopasowanie bez drugiego: (%2i, %2i)\n", aiKoniecPierwszejGrupy[ iBezDrugiego ] + iIndex[ iBezDrugiego ],
1422                              aiKoniecPierwszejGrupy[ 1 - iBezDrugiego ] + iDopasowywany);)
1423
1424                            // zmniejsz liczby niedopasowanych elementow w grupach
1425                            aiRozmiarGrupy[0]--;
1426                    aiRozmiarGrupy[1]--;
1427
1428                    // sprawdz, czy jest szansa dopasowania tego Part z drugiej strony
1429                    // (ktora miala mozliwosc dopasowania tego Part z poprzedniego organizmu)
1430                    if ((aiRozmiarGrupy[0] > 0) && (aiRozmiarGrupy[1] > 0))
1431                    {
1432                        // jesli grupy sie jeszcze nie wyczrpaly
1433                        // to jest mozliwosc dopasowania w organizmie
1434
1435                        int iZDrugim = 1 - iBezDrugiego;
1436                        // sprawdz indeks organizmu
1437                        assert((iZDrugim == 0) || (iZDrugim == 1));
1438
1439                        // bedziemy szukac minimum wsrod niedopasowanych - musimy wykasowac
1440                        // poprzednie obliczenia minimum
1441                        // dla organizmu 1 (o rozmiarze grupy z 0)
1442                        for (i = 0; i < aiRozmiarCalychGrup[0]; i++)
1443                        {
1444                            apvbCzyMinimalnaOdleglosc[1]->operator[](i) = false;
1445                        }
1446                        // dla organizmu 0 (o rozmiarze grup z 1)
1447                        for (i = 0; i < aiRozmiarCalychGrup[1]; i++)
1448                        {
1449                            apvbCzyMinimalnaOdleglosc[0]->operator[](i) = false;
1450                        }
1451
1452                        // szukamy na nowo minimum dla Part o indeksie iIndex[ iZDrugim ] w organizmie iZDrugim
1453                        // wsrod niedopasowanych Parts z organizmu 1 - iZDrugim
1454                        dMinimum = HUGE_VAL;
1455                        for (i = 0; i < aiRozmiarCalychGrup[ 1 - iZDrugim ]; i++)
1456                        {
1457                            if (!(m_pMatching->IsMatched(
1458                                                         1 - iZDrugim,
1459                                                         aiKoniecPierwszejGrupy[ 1 - iZDrugim ] + i)))
1460                            {
1461                                // szukamy minimum obliczonej lokalnej odleglosci tylko wsrod
1462                                // niedopasowanych jeszcze Parts
1463                                if (iZDrugim == 0)
1464                                {
1465                                    // teraz niestety musimy rozpoznac odpowiedni organizm
1466                                    // zeby moc indeksowac niesymetryczna tablice
1467                                    if (aadOdleglosciParts[ iIndex[0] ][ i ] < dMinimum)
1468                                    {
1469                                        dMinimum = aadOdleglosciParts[ iIndex[0] ][ i ];
1470                                    }
1471                                    // przy okazji - sprawdz, czy HUGE_VAL jest rzeczywiscie max dla double
1472                                    assert(aadOdleglosciParts[ iIndex[0] ][ i ] < HUGE_VAL);
1473
1474                                }
1475                                else
1476                                {
1477                                    if (aadOdleglosciParts[ i ][ iIndex[1] ] < dMinimum)
1478                                    {
1479                                        dMinimum = aadOdleglosciParts[ i ][ iIndex[1] ];
1480                                    }
1481                                    // przy okazji - sprawdz, czy HUGE_VAL jest rzeczywiscie max dla double
1482                                    assert(aadOdleglosciParts[ i ][ iIndex[1] ] < HUGE_VAL);
1483                                }
1484                            }
1485                        }
1486                        // sprawdz, czy minimum znaleziono - musi takie byc, bo jest cos niedopasowanego
1487                        assert((dMinimum >= 0.0) && (dMinimum < HUGE_VAL));
1488
1489                        // teraz zaznaczamy w tablicy te wszystkie Parts, ktore maja
1490                        // rzeczywiscie te minimalne odleglosci do Part w organizmie iZDrugim
1491                        for (i = 0; i < aiRozmiarCalychGrup[ 1 - iZDrugim ]; i++)
1492                        {
1493                            if (!(m_pMatching->IsMatched(
1494                                                         1 - iZDrugim,
1495                                                         aiKoniecPierwszejGrupy[ 1 - iZDrugim ] + i)))
1496                            {
1497                                if (iZDrugim == 0)
1498                                {
1499                                    // teraz niestety musimy rozpoznac odpowiedni organizm
1500                                    // zeby moc indeksowac niesymetryczna tablice
1501                                    if (aadOdleglosciParts[ iIndex[0] ][ i ] == dMinimum)
1502                                    {
1503                                        // jesli w danym miejscu tablicy odleglosci jest faktyczne
1504                                        // minimum odleglosci, to mamy potencjalna pare dla iIndex[1]
1505                                        apvbCzyMinimalnaOdleglosc[0]->operator[](i) = true;
1506                                    }
1507                                }
1508                                else
1509                                {
1510                                    if (aadOdleglosciParts[ i ][ iIndex[1] ] == dMinimum)
1511                                    {
1512                                        apvbCzyMinimalnaOdleglosc[1]->operator[](i) = true;
1513                                    }
1514                                }
1515                            }
1516                        }
1517
1518                        // a teraz - po znalezieniu potencjalnych elementow do dopasowania
1519                        // dopasowujemy pierwszy z potencjalnych
1520                        iDopasowywany = -1;
1521                        for (i = 0; i < aiRozmiarCalychGrup[ 1 - iZDrugim ]; i++)
1522                        {
1523                            if (apvbCzyMinimalnaOdleglosc[iZDrugim]->operator[](i))
1524                            {
1525                                iDopasowywany = i;
1526                                break;
1527                            }
1528                        }
1529                        // musielismy znalezc jakiegos dopasowywanego
1530                        assert((iDopasowywany >= 0) &&
1531                               (iDopasowywany < aiRozmiarCalychGrup[ 1 - iZDrugim ]));
1532
1533                        // no to juz mozemy dopasowac
1534                        m_pMatching->Match(
1535                                           iZDrugim,
1536                                           aiKoniecPierwszejGrupy[ iZDrugim ] + iIndex[ iZDrugim ],
1537                                           1 - iZDrugim,
1538                                           aiKoniecPierwszejGrupy[ 1 - iZDrugim ] + iDopasowywany);
1539                        DB(printf("Przypadek 2.1.: dopasowane Parts dopasowaniebz drugim: (%2i, %2i)\n", aiKoniecPierwszejGrupy[ iZDrugim ] + iIndex[ iZDrugim ], aiKoniecPierwszejGrupy[ 1 - iZDrugim ] + iDopasowywany);)
1540
1541                        //aiKoniecPierwszejGrupy[ 1-iBezDrugiego ] + iDopasowywany ;)
1542                        //aiKoniecPierwszejGrupy[ 1-iBezDrugiego ] + iDopasowywany ;)
1543                        // zmniejsz liczby niedopasowanych elementow w grupach
1544                        aiRozmiarGrupy[0]--;
1545                        aiRozmiarGrupy[1]--;
1546
1547                    }
1548                    else
1549                    {
1550                        // jedna z grup sie juz wyczerpala
1551                        // wiec nie ma mozliwosci dopasowania tego drugiego Partu
1552                        /// i trzeba poczekac na zmiane grupy
1553                    }
1554
1555                    DB(printf("Przypadek 2.\n");)
1556
1557                }// PRZYPADEK 2.
1558                else
1559                {
1560                    // PRZYPADEK 3. Zaden z Parts nie ma na liscie drugiego
1561                    // AKCJA. Niezalezne dopasowanie obu Parts do pierwszych ze swojej listy
1562
1563                    // najpierw dopasujemy do iIndex[0] w organizmie 0
1564                    int iDopasowywany = -1;
1565                    // poszukujemy pierwszego z listy dopasowania - w organizmie 1
1566                    for (i = 0; i < aiRozmiarCalychGrup[1]; i++)
1567                    {
1568                        if (apvbCzyMinimalnaOdleglosc[0]->operator[](i))
1569                        {
1570                            iDopasowywany = i;
1571                            break;
1572                        }
1573                    }
1574                    // musielismy znalezc jakiegos dopasowywanego
1575                    assert((iDopasowywany >= 0) &&
1576                           (iDopasowywany < aiRozmiarCalychGrup[1]));
1577
1578                    // teraz wlasnie dopasowujemy
1579                    m_pMatching->Match(
1580                                       0,
1581                                       aiKoniecPierwszejGrupy[0] + iIndex[0],
1582                                       1,
1583                                       aiKoniecPierwszejGrupy[1] + iDopasowywany);
1584
1585                    // zmniejszamy liczbe niedopasowanych Parts
1586                    aiRozmiarGrupy[0]--;
1587                    aiRozmiarGrupy[1]--;
1588
1589                    // debug - dopasowanie
1590                    DB(printf("Przypadek 3.: dopasowane Parts: (%2i, %2i)\n", aiKoniecPierwszejGrupy[0]
1591                              + iIndex[0], aiKoniecPierwszejGrupy[1] + iDopasowywany);)
1592
1593                            // teraz dopasowujemy do iIndex[1] w organizmie 1
1594                            iDopasowywany = -1;
1595                    // poszukujemy pierwszego z listy dopasowania - w organizmie 0
1596                    for (i = 0; i < aiRozmiarCalychGrup[0]; i++)
1597                    {
1598                        if (apvbCzyMinimalnaOdleglosc[1]->operator[](i))
1599                        {
1600                            iDopasowywany = i;
1601                            break;
1602                        }
1603                    }
1604                    // musielismy znalezc jakiegos dopasowywanego
1605                    assert((iDopasowywany >= 0) &&
1606                           (iDopasowywany < aiRozmiarCalychGrup[0]));
1607
1608                    // no i teraz realizujemy dopasowanie
1609                    m_pMatching->Match(
1610                                       0,
1611                                       aiKoniecPierwszejGrupy[0] + iDopasowywany,
1612                                       1,
1613                                       aiKoniecPierwszejGrupy[1] + iIndex[1]);
1614
1615                    // zmniejszamy liczbe niedopasowanych Parts
1616                    aiRozmiarGrupy[0]--;
1617                    aiRozmiarGrupy[1]--;
1618
1619                    // debug - dopasowanie
1620                    DB(printf("Przypadek 3.: dopasowane Parts: (%2i, %2i)\n", aiKoniecPierwszejGrupy[0]
1621                              + iDopasowywany, aiKoniecPierwszejGrupy[1] + iIndex[1]);)
1622
1623
1624                } // PRZYPADEK 3.
1625
1626            }// if (! bCzyKoniecGrupy)
1627            else
1628            {
1629                // gdy mamy już koniec grup - musimy zlikwidować tablicę aadOdleglosciParts
1630                // bo za chwilke skonczy sie nam petla
1631                for (i = 0; i < aiRozmiarCalychGrup[0]; i++)
1632                {
1633                    SAFEDELETEARRAY(aadOdleglosciParts[i]);
1634                }
1635                SAFEDELETEARRAY(aadOdleglosciParts);
1636
1637                // musimy tez usunac tablice (wektory) mozliwosci dopasowania
1638                SAFEDELETE(apvbCzyMinimalnaOdleglosc[0]);
1639                SAFEDELETE(apvbCzyMinimalnaOdleglosc[1]);
1640            }
1641        } // while (! bCzyKoniecGrupy)
1642
1643        // PO DOPASOWANIU WSZYSTKIEGO Z GRUP (CO NAJMNIEJ JEDNEJ GRUPY W CAŁOŚCI)
1644
1645        // gdy rozmiar ktorejkolwiek z grup dopasowania spadl do zera
1646        // to musimy przesunac KoniecPierwszejGrupy (wszystkie dopasowane)
1647        for (i = 0; i < 2; i++)
1648        {
1649            // grupy nie moga miec mniejszego niz 0 rozmiaru
1650            assert(aiRozmiarGrupy[i] >= 0);
1651            if (aiRozmiarGrupy[i] == 0)
1652                aiKoniecPierwszejGrupy[i] = aiKoniecGrupyDopasowania[i];
1653        }
1654
1655        // sprawdzenie warunku konca dopasowywania - gdy nie
1656        // ma juz w jakims organizmie co dopasowywac
1657        if (aiKoniecPierwszejGrupy[0] >= m_aiPartCount[0] ||
1658                aiKoniecPierwszejGrupy[1] >= m_aiPartCount[1])
1659        {
1660            iCzyDopasowane = 1;
1661        }
1662    } // koniec WHILE - petli dopasowania
1663}
1664
1665/** Matches Parts in both organisms so that computation of similarity is possible.
1666    New algorithm (assures symmetry of the similarity measure) with geometry
1667    taken into consideration.
1668    Assumptions:
1669    - Models (m_Mod) are created and available.
1670        - Matching (m_pMatching) is created, but empty
1671        Exit conditions:
1672        - Matching (m_pMatching) is full
1673        @return 1 if success, 0 otherwise
1674 */
1675int ModelSimil::MatchPartsGeometry()
1676{
1677    // zaloz, ze sa modele i sa poprawne
1678    assert((m_Mod[0] != NULL) && (m_Mod[1] != NULL));
1679    assert(m_Mod[0]->isValid() && m_Mod[1]->isValid());
1680
1681    if (!CreatePartInfoTables())
1682        return 0;
1683    if (!CountPartDegrees())
1684        return 0;
1685    if (!GetPartPositions())
1686    {
1687        return 0;
1688    }
1689    if (!CountPartNeurons())
1690        return 0;
1691
1692
1693    if (m_adFactors[3] > 0)
1694    {
1695        if (!ComputePartsPositionsBySVD())
1696        {
1697            return 0;
1698        }
1699    }
1700
1701    DB(printf("Przed sortowaniem:\n");)
1702    DB(_PrintDegrees(0);)
1703    DB(_PrintDegrees(1);)
1704
1705    if (!SortPartInfoTables())
1706        return 0;
1707
1708    DB(printf("Po sortowaniu:\n");)
1709    DB(_PrintDegrees(0);)
1710    DB(_PrintDegrees(1);)
1711
1712    if (m_adFactors[3] > 0)
1713    {
1714        // tutaj zacznij pętlę po przekształceniach  geometrycznych
1715        const int NO_OF_TRANSFORM = 8; // liczba transformacji geometrycznych (na razie tylko ID i O_YZ)
1716        // tablice transformacji współrzędnych; nie są to dokładnie tablice tranformacji, ale raczej tablice PRZEJŚĆ
1717        // pomiędzy transformacjami;
1718        // wartości orginalne transformacji dOrig uzyskuje się przez:
1719        // for ( iTrans = 0; iTrans <= TRANS_INDEX; iTrans++ ) dOrig *= dMul[ iTrans ];
1720        //const char *szTransformNames[NO_OF_TRANSFORM] = { "ID", "S_yz", "S_xz", "S_xy", "R180_z", "R180_y", "R180_z", "S_(0,0,0)" };
1721        const int dMulX[ NO_OF_TRANSFORM ] = {1, -1, -1, 1, -1, 1, -1, -1};
1722        const int dMulY[ NO_OF_TRANSFORM ] = {1, 1, -1, -1, -1, -1, -1, 1};
1723        const int dMulZ[ NO_OF_TRANSFORM ] = {1, 1, 1, -1, -1, -1, 1, 1};
1724
1725#ifdef max
1726 #undef max //this macro would conflict with line below
1727#endif
1728        double dMinSimValue = std::numeric_limits<double>::max(); // minimum value of similarity
1729        int iMinSimTransform = -1; // index of the transformation with the minimum similarity
1730        SimilMatching *pMinSimMatching = NULL; // matching with the minimum value of similarity
1731
1732        // remember the original positions of model 0 as computed by SVD in order to restore them later, after
1733        // all transformations have been computed
1734        Pt3D *StoredPositions = NULL; // array of positions of Parts, for one (0th) model
1735        // create the stored positions
1736        StoredPositions = new Pt3D[ m_Mod[ 0 ]->getPartCount() ];
1737        assert(StoredPositions != NULL);
1738        // copy the original positions of model 0 (store them)
1739        int iPart; // a counter of Parts
1740        for (iPart = 0; iPart < m_Mod[ 0 ]->getPartCount(); iPart++)
1741        {
1742            StoredPositions[ iPart ].x = m_aPositions[ 0 ][ iPart ].x;
1743            StoredPositions[ iPart ].y = m_aPositions[ 0 ][ iPart ].y;
1744            StoredPositions[ iPart ].z = m_aPositions[ 0 ][ iPart ].z;
1745        }
1746        // now the original positions of model 0 are stored
1747
1748
1749        int iTransform; // a counter of geometric transformations
1750        for (iTransform = 0; iTransform < NO_OF_TRANSFORM; iTransform++)
1751        {
1752            // for each geometric transformation to be done
1753            // entry conditions:
1754            // - models (m_Mod) exist and are available
1755            // - matching (m_pMatching) exists and is empty
1756            // - all properties are created and available (m_aDegrees and m_aPositions)
1757
1758            // recompute geometric properties according to the transformation iTransform
1759            // but only for model 0
1760            for (iPart = 0; iPart < m_Mod[ 0 ]->getPartCount(); iPart++)
1761            {
1762                // for each iPart, a part of the model iMod
1763                m_aPositions[ 0 ][ iPart ].x *= dMulX[ iTransform ];
1764                m_aPositions[ 0 ][ iPart ].y *= dMulY[ iTransform ];
1765                m_aPositions[ 0 ][ iPart ].z *= dMulZ[ iTransform ];
1766            }
1767            // now the positions are recomputed
1768            ComputeMatching();
1769
1770            // teraz m_pMatching istnieje i jest pełne
1771            assert(m_pMatching != NULL);
1772            assert(m_pMatching->IsFull() == true);
1773
1774            // wykorzystaj to dopasowanie i geometrię do obliczenia tymczasowej wartości miary
1775            int iTempRes = CountPartsDistance();
1776            // załóż sukces
1777            assert(iTempRes == 1);
1778            double dCurrentSim = m_adFactors[0] * double(m_iDV) +
1779                    m_adFactors[1] * double(m_iDD) +
1780                    m_adFactors[2] * double(m_iDN) +
1781                    m_adFactors[3] * double(m_dDG);
1782            // załóż poprawną wartość podobieństwa
1783            assert(dCurrentSim >= 0.0);
1784
1785            // porównaj wartość obliczoną z dotychczasowym minimum
1786            if (dCurrentSim < dMinSimValue)
1787            {
1788                // jeśli uzyskano mniejszą wartość dopasowania,
1789                // to zapamiętaj to przekształcenie geometryczne i dopasowanie
1790                dMinSimValue = dCurrentSim;
1791                iMinSimTransform = iTransform;
1792                SAFEDELETE(pMinSimMatching);
1793                pMinSimMatching = new SimilMatching(*m_pMatching);
1794                assert(pMinSimMatching != NULL);
1795            }
1796
1797            // teraz już można usunąć stare dopasowanie (dla potrzeb następnego przebiegu pętli)
1798            m_pMatching->Empty();
1799        } // for ( iTransform )
1800
1801        // po pętli przywróć najlepsze dopasowanie
1802        delete m_pMatching;
1803        m_pMatching = pMinSimMatching;
1804
1805        DB(printf("Matching has been chosen!\n");)
1806                // print the name of the chosen transformation:
1807                // printf("Chosen transformation: %s\n", szTransformNames[ iMinSimTransform ] );
1808
1809                // i przywróć jednocześnie pozycje Parts modelu 0 dla tego dopasowania
1810                // - najpierw przywroc Parts pozycje orginalne, po SVD
1811        for (iPart = 0; iPart < m_Mod[ 0 ]->getPartCount(); iPart++)
1812        {
1813            m_aPositions[ 0 ][ iPart ].x = StoredPositions[ iPart ].x;
1814            m_aPositions[ 0 ][ iPart ].y = StoredPositions[ iPart ].y;
1815            m_aPositions[ 0 ][ iPart ].z = StoredPositions[ iPart ].z;
1816        }
1817        // - usun teraz zapamietane pozycje Parts
1818        delete [] StoredPositions;
1819        // - a teraz oblicz na nowo wspolrzedne wlasciwego przeksztalcenia dla model 0
1820        for (iTransform = 0; iTransform <= iMinSimTransform; iTransform++)
1821        {
1822            // for each transformation before and INCLUDING iMinTransform
1823            // do the transformation (only model 0)
1824            for (iPart = 0; iPart < m_Mod[ 0 ]->getPartCount(); iPart++)
1825            {
1826                m_aPositions[ 0 ][ iPart ].x *= dMulX[ iTransform ];
1827                m_aPositions[ 0 ][ iPart ].y *= dMulY[ iTransform ];
1828                m_aPositions[ 0 ][ iPart ].z *= dMulZ[ iTransform ];
1829            }
1830        }
1831
1832    }
1833    else
1834    {
1835        ComputeMatching();
1836    }
1837    // teraz dopasowanie musi byc pelne - co najmniej w jednym organizmie musza byc
1838    // wszystkie elementy dopasowane
1839    assert(m_pMatching->IsFull() == true);
1840
1841    //    _PrintDegrees(0);
1842    //    _PrintDegrees(1);
1843
1844    DB(_PrintPartsMatching();)
1845
1846
1847    return 1;
1848}
1849
1850void ModelSimil::_PrintSeamnessTable(std::vector<int> *pTable, int iCount)
1851{
1852    int i;
1853    printf("      ");
1854    for (i = 0; i < iCount; i++)
1855        printf("%3i ", i);
1856    printf("\n      ");
1857    for (i = 0; i < iCount; i++)
1858    {
1859
1860        printf("%3i ", pTable->operator[](i));
1861    }
1862    printf("\n");
1863}
1864
1865/** Computes elements of similarity (m_iDD, m_iDN, m_dDG) based on underlying matching.
1866    Assumptions:
1867    - Matching (m_pMatching) exists and is full.
1868        - Internal arrays m_aDegrees and m_aPositions exist and are properly filled in
1869        Exit conditions:
1870        - Elements of similarity are computed (m)iDD, m_iDN, m_dDG).
1871        @return 1 if success, otherwise 0.
1872 */
1873int ModelSimil::CountPartsDistance()
1874{
1875    int i, temp;
1876
1877    // assume existence of m_pMatching
1878    assert(m_pMatching != NULL);
1879    // musi byc pelne!
1880    assert(m_pMatching->IsFull() == true);
1881
1882    // roznica w stopniach
1883    m_iDD = 0;
1884    // roznica w liczbie neuronów
1885    m_iDN = 0;
1886    // roznica w odleglosci dopasowanych Parts
1887    m_dDG = 0.0;
1888
1889    int iOrgPart, iOrgMatchedPart; // orginalny indeks Part i jego dopasowanego Part
1890    int iMatchedPart; // indeks (wg sortowania) dopasowanego Part
1891
1892    // wykorzystanie dopasowania do zliczenia m_iDD - roznicy w stopniach
1893    // i m_iDN - roznicy w liczbie neuronow
1894    // petla w wiekszej tablicy
1895    for (i = 0; i < m_aiPartCount[1 - m_iSmaller]; i++)
1896    {
1897        // dla kazdego elementu [i] z wiekszego organizmu
1898        // pobierz jego orginalny indeks w organizmie z tablicy TDN
1899        iOrgPart = m_aDegrees[ 1 - m_iSmaller ][ i ][ 0 ];
1900        if (!(m_pMatching->IsMatched(1 - m_iSmaller, i)))
1901        {
1902            // gdy nie zostal dopasowany
1903            // dodaj jego stopien do DD
1904            m_iDD += m_aDegrees[1 - m_iSmaller][i][1];
1905            // dodaj liczbe neuronow do DN
1906            m_iDN += m_aDegrees[1 - m_iSmaller][i][3];
1907            // i oblicz odleglosc tego Part od srodka organizmu (0,0,0)
1908            // (uzyj orginalnego indeksu)
1909            //no need to compute distane when m_dDG weight is 0
1910            m_dDG += m_adFactors[3] == 0 ? 0 : m_aPositions[ 1 - m_iSmaller ][ iOrgPart ].length();
1911        }
1912        else
1913        {
1914            // gdy byl dopasowany
1915            // pobierz indeks (po sortowaniu) tego dopasowanego Part
1916            iMatchedPart = m_pMatching->GetMatchedIndex(1 - m_iSmaller, i);
1917            // pobierz indeks orginalny tego dopasowanego Part
1918            iOrgMatchedPart = m_aDegrees[ m_iSmaller ][ iMatchedPart ][0];
1919            // dodaj ABS roznicy stopni do DD
1920            temp = m_aDegrees[1 - m_iSmaller][i][1] -
1921                    m_aDegrees[ m_iSmaller ][ iMatchedPart ][1];
1922            m_iDD += abs(temp);
1923            // dodaj ABS roznicy neuronow do DN
1924            temp = m_aDegrees[1 - m_iSmaller][i][3] -
1925                    m_aDegrees[ m_iSmaller ][ iMatchedPart ][3];
1926            m_iDN += abs(temp);
1927            // pobierz polozenie dopasowanego Part
1928            Pt3D MatchedPartPos(m_aPositions[ m_iSmaller ][ iOrgMatchedPart ]);
1929            // dodaj euklidesowa odleglosc Parts do sumy odleglosci
1930            //no need to compute distane when m_dDG weight is 0
1931            m_dDG +=m_adFactors[3] == 0 ? 0 :  m_aPositions[ 1 - m_iSmaller ][ iOrgPart ].distanceTo(MatchedPartPos);
1932        }
1933    }
1934
1935    // obliczenie i dodanie różnicy w liczbie neuronów OnJoint...
1936    temp = m_aOnJoint[0][3] - m_aOnJoint[1][3];
1937    m_iDN += abs(temp);
1938    DB(printf("OnJoint DN: %i\n", abs(temp));)
1939    // ... i Anywhere
1940    temp = m_aAnywhere[0][3] - m_aAnywhere[1][3];
1941    m_iDN += abs(temp);
1942    DB(printf("Anywhere DN: %i\n", abs(temp));)
1943
1944    return 1;
1945}
1946
1947/** Computes new positions of Parts of both of organisms stored in the object.
1948        Assumptions:
1949        - models (m_Mod) are created and valid
1950        - positions (m_aPositions) are created and filled with original positions of Parts
1951        - the sorting algorithm was not yet run on the array m_aDegrees
1952        @return true if successful; false otherwise
1953 */
1954bool ModelSimil::ComputePartsPositionsBySVD()
1955{
1956    bool bResult = true; // the result; assume a success
1957
1958    // check assumptions
1959    // the models
1960    assert(m_Mod[0] != NULL && m_Mod[0]->isValid());
1961    assert(m_Mod[1] != NULL && m_Mod[1]->isValid());
1962    // the position arrays
1963    assert(m_aPositions[0] != NULL);
1964    assert(m_aPositions[1] != NULL);
1965
1966    int iMod; // a counter of models
1967    // use SVD to obtain different point of view on organisms
1968    // and store the new positions (currently the original ones are still stored)
1969    for (iMod = 0; iMod < 2; iMod++)
1970    {
1971        // prepare the vector of errors of approximation for the SVD
1972        std::vector<double> vEigenvalues;
1973        int nSize = m_Mod[ iMod ]->getPartCount();
1974
1975        double *pDistances = (double *) malloc(nSize * nSize * sizeof (double));
1976
1977        for (int i = 0; i < nSize; i++)
1978        {
1979            pDistances[i] = 0;
1980        }
1981
1982        Model *pModel = m_Mod[ iMod ]; // use the model of the iMod (current) organism
1983        int iP1, iP2; // indices of Parts in the model
1984        Part *P1, *P2; // pointers to Parts
1985        Pt3D P1Pos, P2Pos; // positions of parts
1986        double dDistance; // the distance between Parts
1987        for (iP1 = 0; iP1 < pModel->getPartCount(); iP1++)
1988        {
1989            // for each iP1, a Part index in the model of organism iMod
1990            P1 = pModel->getPart(iP1);
1991            // get the position of the Part
1992            P1Pos = P1->p;
1993            for (iP2 = 0; iP2 < pModel->getPartCount(); iP2++)
1994            {
1995                // for each (iP1, iP2), a pair of Parts index in the model
1996                P2 = pModel->getPart(iP2);
1997                // get the position of the Part
1998                P2Pos = P2->p;
1999                // compute the geometric (Euclidean) distance between the Parts
2000                dDistance = P1Pos.distanceTo(P2Pos);
2001                // store the distance
2002                pDistances[iP1 * nSize + iP2] = dDistance;
2003            } // for (iP2)
2004        } // for (iP1)
2005
2006        MatrixTools::SVD(vEigenvalues, nSize, pDistances, m_aPositions[ iMod ]);
2007        free(pDistances);
2008    }
2009
2010    return bResult;
2011}
2012
2013void ModelSimil::p_evaldistance(ExtValue *args, ExtValue *ret)
2014{
2015    Geno *g1 = GenoObj::fromObject(args[1]);
2016    Geno *g2 = GenoObj::fromObject(args[0]);
2017    if ((!g1) || (!g2))
2018        ret->setEmpty();
2019    else
2020        ret->setDouble(EvaluateDistance(g1, g2));
2021}
Note: See TracBrowser for help on using the repository browser.