Context Navigation

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

Last change on this file since 870 was 869, checked in by Maciej Komosinski, 5 years ago
Added another, improved way of calculating dissimilarity of two creatures/models. Details and comparisons in https://doi.org/10.1007/978-3-030-16692-2_8
Property svn:eol-style set to `native`
File size: 69.2 KB

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats: