Context Navigation

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

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

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

Download in other formats: