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

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