Context Navigation

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

Last change on this file since 978 was 972, checked in by Maciej Komosinski, 4 years ago
separate "0" and "0s" formats (for SHAPE_BALL_AND_STICK and SHAPE_SOLIDS, respectively) converting to format list (Geno::F0_FORMAT_LIST = "0,0s") (optional) declaring Model as SHAPE_BALL_AND_STICK or SHAPE_SOLIDS (or SHAPE_UNKNOWN)
Property svn:eol-style set to `native`
File size: 68.8 KB

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

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

Download in other formats: