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

Last change on this file since 352 was 352, checked in by oriona, 9 years ago

comparison between signed and unsigned integer expressions fixed

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