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

Rev	Line
[349]	1	// This file is a part of Framsticks SDK. http://www.framsticks.com/
[869]	2	// Copyright (C) 1999-2019 Maciej Komosinski and Szymon Ulatowski.
[349]	3	// See LICENSE.txt for details.
	4
[869]	5	// implementation of the ModelSimil class.
	6
[349]	7	#include "SVD/matrix_tools.h"
[869]	8	#include "hungarian/hungarian.h"
[349]	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>
[492]	15	#ifdef EMSCRIPTEN
[606]	16	#include <cstdlib>
[492]	17	#else
[606]	18	#include <stdlib.h>
[492]	19	#endif
[349]	20	#include <math.h>
	21	#include <string>
	22	#include <limits>
	23	#include <assert.h>
	24	#include <vector>
	25	#include <algorithm>
[666]	26	#include <cstdlib> //required for std::qsort in macos xcode
[349]	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[] = {
[869]	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), "",},
[606]	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), "", },
[818]	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.", },
[606]	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, },
[349]	47	};
	48
	49	#undef FIELDSTRUCT
	50
	51	//////////////////////////////////////////////////////////////////////
	52	// Construction/Destruction
	53	//////////////////////////////////////////////////////////////////////
	54
	55	/** Constructor. Sets default weights. Initializes other fields with zeros.
	56	*/
[356]	57	ModelSimil::ModelSimil() : localpar(MSparam_tab, this), m_iDV(0), m_iDD(0), m_iDN(0), m_dDG(0.0)
[349]	58	{
[606]	59	localpar.setDefault();
[349]	60
[606]	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;
[349]	74
[606]	75	//Determines whether "fuzzy vertex degree" should be used.
	76	//Currently "fuzzy vertex degree" is inactive.
[872]	77	isFuzzy = false;
[606]	78	fuzzyDepth = 10;
[869]	79
[817]	80	//Determines whether weighted MDS should be used.
	81	wMDS = 0;
[869]	82	//Determines whether best matching should be saved using hungarian similarity measure.
[872]	83	saveMatching = false;
[349]	84	}
	85
[869]	86	double ModelSimil::EvaluateDistanceGreedy(const Geno G0, const Geno G1)
[349]	87	{
[606]	88	double dResult = 0;
[349]	89
[606]	90	m_Gen[0] = G0;
	91	m_Gen[1] = G1;
[349]	92
[606]	93	// check whether pointers are not NULL
	94	if (m_Gen[0] == NULL \|\| m_Gen[1] == NULL)
	95	{
[871]	96	logPrintf("ModelSimil", "EvaluateDistanceGreedy", LOG_ERROR, "NULL genotype pointer(s)");
[869]	97	return 0.0;
[606]	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]));
[349]	102
[606]	103	// validate models
	104	if (m_Mod[0] == NULL \|\| m_Mod[1] == NULL \|\| !(m_Mod[0]->isValid()) \|\| !(m_Mod[1]->isValid()))
	105	{
[871]	106	logPrintf("ModelSimil", "EvaluateDistanceGreedy", LOG_ERROR, "NULL or invalid model pointer(s)");
[872]	107	return 0.0;
[606]	108	}
[349]	109
[606]	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	}
[349]	120
[606]	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);
[349]	125
[606]	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);
[349]	131
	132
[606]	133	// assign matching function
	134	int (ModelSimil::* pfMatchingFunction) () = NULL;
[349]	135
[606]	136	pfMatchingFunction = &ModelSimil::MatchPartsGeometry;
[349]	137
[606]	138	// match Parts (vertices of creatures)
	139	if ((this->*pfMatchingFunction)() == 0)
	140	{
[872]	141	logPrintf("ModelSimil", "EvaluateDistanceGreedy", LOG_ERROR, "The matching function returned 0");
	142	return 0.0;
[606]	143	}
[349]	144
[606]	145	// after matching function call we must have full matching
	146	assert(m_pMatching->IsFull() == true);
[349]	147
[606]	148	DB(SaveIntermediateFiles();)
[349]	149
[606]	150	// count differences in matched parts
	151	if (CountPartsDistance() == 0)
	152	{
[872]	153	logPrintf("ModelSimil", "EvaluateDistanceGreedy", LOG_ERROR, "CountPartDistance()==0");
	154	return 0.0;
[606]	155	}
[349]	156
[606]	157	// delete degree arrays created in CreatePartInfoTables
	158	SAFEDELETEARRAY(m_aDegrees[0]);
	159	SAFEDELETEARRAY(m_aDegrees[1]);
[349]	160
[606]	161	// and position arrays
	162	SAFEDELETEARRAY(m_aPositions[0]);
	163	SAFEDELETEARRAY(m_aPositions[1]);
[349]	164
[606]	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	}
[349]	178
[606]	179	}
[349]	180
[606]	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;
[349]	194	}
	195
	196	ModelSimil::~ModelSimil()
	197	{
[606]	198	// matching should have been deleted earlier
	199	assert(m_pMatching == NULL);
[349]	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	{
[606]	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);
[349]	213
[606]	214	int iOrgPart; // original index of a Part
	215	int nBigger; // index of the larger organism
[349]	216
[606]	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	}
[349]	226
[606]	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);
[349]	234
[606]	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)
[349]	242
[606]	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());
[349]	256
[606]	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 )
[349]	275
[606]	276	// now all matched Part indices are printed out, end the line
	277	fprintf(pMatchingFile, "] : ORG[%i]\n", 1 - nBigger);
[349]	278
[606]	279	// close the file
	280	fclose(pMatchingFile);
	281	// END TEMP
[349]	282
[606]	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	}
[349]	316	}
	317
	318	/** Comparison function required for qsort() call. Used while sorting groups of
[606]	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	*/
[349]	325	int ModelSimil::CompareDegrees(const void pElem1, const void pElem2)
	326	{
[606]	327	int tdn1 = (int )pElem1;
	328	int tdn2 = (int )pElem2;
[349]	329
[606]	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	{
[869]	338	// when degree - tdn2[1] - is BIGGER
	339	return 1;
[606]	340	}
	341	else
	342	{
	343	return 0;
	344	}
[349]	345	}
	346
[869]	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
[349]	376	/** Comparison function required for qsort() call. Used while sorting groups of Parts with
[606]	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	*/
[349]	384	int ModelSimil::CompareConnsNo(const void pElem1, const void pElem2)
	385	{
[606]	386	// pointers to TDN arrays
	387	int tdn1, tdn2;
	388	// definitions of elements being compared
	389	tdn1 = (int *)pElem1;
	390	tdn2 = (int *)pElem2;
[349]	391
[606]	392	// comparison according to Neural Connections (to jest TDN[2])
[782]	393	if (tdn1[NEURO_CONNS] > tdn2[NEURO_CONNS])
[606]	394	{
	395	// when number of NConn Elem1 is BIGGER
	396	return -1;
	397	}
	398	else
[782]	399	if (tdn1[NEURO_CONNS] < tdn2[NEURO_CONNS])
[606]	400	{
[869]	401	// when number of NConn Elem1 is SMALLER
	402	return 1;
[606]	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	{
[869]	416	// when number of Neu is SMALLER for Elem1
	417	return 1;
[606]	418	}
	419	else
	420	{
	421	// when numbers of Nconn and Neu are equal we check original indices
	422	// of Parts being compared
[349]	423
[606]	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	{
[869]	433	// when the number of NIt Deg1 id SMALLER
	434	return 1;
[606]	435	}
	436	else
	437	{
	438	// impossible, indices are alway different
	439	return 0;
	440	}
	441	}
	442	}
[349]	443	}
	444
[606]	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	*/
[349]	449	int ModelSimil::GetNOFactors()
	450	{
[606]	451	return ModelSimil::iNOFactors;
[349]	452	}
	453
	454	/** Prints the array of degrees for the given organism. Debug method.
	455	*/
	456	void ModelSimil::_PrintDegrees(int i)
	457	{
[606]	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");
[349]	476	}
	477
	478	/** Prints one array of ints. Debug method.
[606]	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	*/
[349]	483	void ModelSimil::_PrintArray(int *array, int base, int size)
	484	{
[606]	485	int i;
	486	for (i = base; i < base + size; i++)
	487	{
	488	printf("%i ", array[i]);
	489	}
	490	printf("\n");
[349]	491	}
	492
	493	void ModelSimil::_PrintArrayDouble(double *array, int base, int size)
	494	{
[606]	495	int i;
	496	for (i = base; i < base + size; i++)
	497	{
	498	printf("%f ", array[i]);
	499	}
	500	printf("\n");
[349]	501	}
	502
	503	/** Prints one array of parts neighbourhood.
[606]	504	@param index of organism
	505	*/
[349]	506	void ModelSimil::_PrintNeighbourhood(int o)
	507	{
[606]	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	}
[349]	519	}
	520
[869]	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
[349]	539	/** Creates arrays holding information about organisms' Parts (m_aDegrees) andm_Neigh
[606]	540	fills them with initial data (original indices and zeros).
	541	Assumptions:
	542	- Models (m_Mod) are created and available.
	543	*/
[349]	544	int ModelSimil::CreatePartInfoTables()
	545	{
[606]	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());
[349]	549
[606]	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];
[349]	557
[606]	558	if (isFuzzy)
	559	{
	560	m_Neighbours[i] = new int*[partCount];
	561	m_fuzzyNeighb[i] = new float*[partCount];
	562	}
[349]	563
[872]	564	if (m_aDegrees[i] != NULL && ((!isFuzzy) \|\| (m_Neighbours[i] != NULL && m_fuzzyNeighb[i] != NULL)))
[606]	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	{
[869]	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;
[349]	576
[869]	577	// sprawdz, czy nie piszemy po jakims szalonym miejscu pamieci
	578	assert(m_aDegrees[i][j] != NULL);
[349]	579
[869]	580	if (isFuzzy)
[606]	581	{
[869]	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	}
[349]	587
[869]	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);
[606]	596	}
[349]	597
[606]	598	}
	599	}
	600	else
	601	{
[872]	602	logPrintf("ModelSimil", "CreatePartInfoTables", LOG_ERROR, "Not enough memory?");
	603	return 0;
[606]	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;
[349]	621	}
	622
	623	/** Computes degrees of Parts of both organisms. Fills in the m_aDegrees arrays
[606]	624	with proper information about degrees.
	625	Assumptions:
	626	- Models (m_Mod) are created and available.
	627	- Arrays m_aDegrees are created.
	628	*/
[349]	629	int ModelSimil::CountPartDegrees()
	630	{
[606]	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());
[349]	634
[606]	635	// sprawdz zalozenie - o tablicach
	636	assert(m_aDegrees[0] != NULL);
	637	assert(m_aDegrees[1] != NULL);
[349]	638
[606]	639	Part P1, P2;
	640	int i, j, i1, i2;
[349]	641
[606]	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	}
[349]	670
[606]	671	if (isFuzzy)
	672	{
	673	CountFuzzyNeighb();
	674	}
[349]	675
[606]	676	return 1;
[349]	677	}
	678
	679	void ModelSimil::GetNeighbIndexes(int mod, int partInd, std::vector<int> &indexes)
	680	{
[606]	681	indexes.clear();
	682	int i, size = m_Mod[mod]->getPartCount();
[349]	683
[606]	684	for (i = 0; i < size; i++)
	685	{
	686	if (m_Neighbours[mod][partInd][i] > 0)
	687	{
	688	indexes.push_back(i);
	689	}
	690	}
[349]	691	}
	692
	693	int cmpFuzzyRows(const void pa, const void pb)
	694	{
[606]	695	float a = (float)pa;
	696	float b = (float)pb;
[349]	697
	698
[606]	699	for (int i = 1; i < fuzDepth; i++)
	700	{
	701	if (a[0][i] > b[0][i])
	702	{
[349]	703
[606]	704	return -1;
	705	}
	706	if (a[0][i] < b[0][i])
	707	{
[349]	708
[606]	709	return 1;
	710	}
	711	}
[349]	712
[606]	713	return 0;
[349]	714	}
	715
[869]	716	void ModelSimil::FuzzyOrder()
[349]	717	{
[869]	718	int i, depth, partInd, prevPartInd, partCount;
[606]	719	for (int mod = 0; mod < 2; mod++)
	720	{
	721	partCount = m_Mod[mod]->getPartCount();
[869]	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--)
[606]	726	{
[869]	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++)
[606]	731	{
[869]	732	if (m_fuzzyNeighb[mod][i][depth] != m_fuzzyNeighb[mod][i + 1][depth])
[606]	733	{
[869]	734	m_aDegrees[mod][partInd][FUZZ_DEG]++;
[606]	735	break;
	736	}
	737	}
	738	}
	739	}
[349]	740	}
	741
	742	//sort according to fuzzy degree
	743	void ModelSimil::SortFuzzyNeighb()
	744	{
[606]	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	}
[349]	750	}
	751
	752	//computes fuzzy vertex degree
	753	void ModelSimil::CountFuzzyNeighb()
	754	{
[606]	755	assert(m_aDegrees[0] != NULL);
	756	assert(m_aDegrees[1] != NULL);
[349]	757
[606]	758	assert(m_Neighbours[0] != NULL);
	759	assert(m_Neighbours[1] != NULL);
[349]	760
[606]	761	assert(m_fuzzyNeighb[0] != NULL);
	762	assert(m_fuzzyNeighb[1] != NULL);
[349]	763
[606]	764	std::vector<int> nIndexes;
	765	float newDeg = 0;
[349]	766
[606]	767	for (int mod = 0; mod < 2; mod++)
	768	{
	769	int partCount = m_Mod[mod]->getPartCount();
[349]	770
[606]	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);
[361]	787	for (unsigned int k = 0; k < nIndexes.size(); k++)
[606]	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	{
[349]	798
[606]	799	}
	800	}
	801	newDeg = 0;
	802	}
	803	}
	804	}
	805	}
[349]	806
[606]	807	SortFuzzyNeighb();
[869]	808	FuzzyOrder();
[349]	809	}
	810
	811	/** Gets information about Parts' positions in 3D from models into the arrays
[606]	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	*/
[349]	818	int ModelSimil::GetPartPositions()
	819	{
[606]	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());
[349]	823
[606]	824	// sprawdz zalozenie - o tablicach m_aPositions
	825	assert(m_aPositions[0] != NULL);
	826	assert(m_aPositions[1] != NULL);
[349]	827
[606]	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	}
[349]	847
[606]	848	return 1;
[349]	849	}
	850
	851	/** Computes numbers of neurons and neurons' inputs for each Part of each
[606]	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	*/
[349]	857	int ModelSimil::CountPartNeurons()
	858	{
[606]	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());
[349]	862
[606]	863	// sprawdz zalozenie - o tablicach
	864	assert(m_aDegrees[0] != NULL);
	865	assert(m_aDegrees[1] != NULL);
[349]	866
[606]	867	Part *P1;
	868	Joint *J1;
	869	int i, j, i2, neuro_connections;
[349]	870
[606]	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;
[349]	911	}
	912
	913	/** Sorts arrays m_aDegrees (for each organism) by Part's degree, and then by
[606]	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	*/
[349]	921	int ModelSimil::SortPartInfoTables()
	922	{
[606]	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());
[349]	926
[606]	927	// sprawdz zalozenie - o tablicach
	928	assert(m_aDegrees[0] != NULL);
	929	assert(m_aDegrees[1] != NULL);
[349]	930
[606]	931	int i;
[869]	932	int(pfDegreeFunction) (const void, const void*) = NULL;
[872]	933	pfDegreeFunction = isFuzzy ? &CompareFuzzyDegrees : &CompareDegrees;
[606]	934	// sortowanie obu tablic wg stopni punktów - TDN[1]
[869]	935	for (i = 0; i < 2; i++)
[606]	936	{
[869]	937	DB(_PrintDegrees(i));
	938	std::qsort(m_aDegrees[i], (size_t)(m_Mod[i]->getPartCount()),
	939	sizeof(TDN), pfDegreeFunction);
	940	DB(_PrintDegrees(i));
[606]	941	}
[349]	942
[606]	943	// sprawdzenie wartosci parametru
	944	DB(i = sizeof(TDN);)
[872]	945	int degreeType = isFuzzy ? FUZZ_DEG : DEGREE;
[349]	946
[606]	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 <
[349]	962
[606]	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);)
[349]	974
[606]	975	std::qsort(m_aDegrees[i][iPocz], (size_t)(j - iPocz),
[869]	976	sizeof(TDN), ModelSimil::CompareConnsNo);
[606]	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;
[349]	987	}
	988
	989
	990	/** Prints the state of the matching object. Debug method.
	991	*/
	992	void ModelSimil::_PrintPartsMatching()
	993	{
[606]	994	// assure that matching exists
	995	assert(m_pMatching != NULL);
[349]	996
[606]	997	printf("Parts matching:\n");
	998	m_pMatching->PrintMatching();
[349]	999	}
	1000
	1001	void ModelSimil::ComputeMatching()
	1002	{
[606]	1003	// uniwersalne liczniki
	1004	int i, j;
[349]	1005
[606]	1006	assert(m_pMatching != NULL);
	1007	assert(m_pMatching->IsEmpty() == true);
[349]	1008
[606]	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 };
[349]	1038
[606]	1039	// utworzenie tablicy rozmiarow
	1040	for (i = 0; i < 2; i++)
	1041	{
	1042	m_aiPartCount[i] = m_Mod[i]->getPartCount();
	1043	}
[349]	1044
[606]	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;
[349]	1060
[606]	1061	// sprawdz poprawnosc tego indeksu
	1062	assert((aiKoniecGrupyDopasowania[i] > 0) &&
	1063	(aiKoniecGrupyDopasowania[i] <= m_aiPartCount[i]));
[349]	1064
[606]	1065	// oblicz rozmiary całych grup - łącznie z dopasowanymi już elementami
	1066	aiRozmiarCalychGrup[i] = aiKoniecGrupyDopasowania[i] -
	1067	aiKoniecPierwszejGrupy[i];
[349]	1068
[606]	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	}
[349]	1084
[606]	1085	// DOPASOWYWANIE PARTS Z GRUP
[349]	1086
[606]	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	}
[349]	1097
[606]	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);
[349]	1104
[606]	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)
[349]	1117
[606]	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],
[869]	1122	aiKoniecGrupyDopasowania[1]);)
[349]	1123
[606]	1124	for (i = 0; i < aiRozmiarCalychGrup[0]; i++)
	1125	{
[349]	1126
[869]	1127	// iterujemy i - Parts organizmu 0
	1128	// (indeks podstawowy - aiKoniecPierwszejGrupy[0])
[349]	1129
[869]	1130	if (!(m_pMatching->IsMatched(0, aiKoniecPierwszejGrupy[0] + i)))
[606]	1131	{
[869]	1132	// interesuja nas tylko te niedopasowane jeszcze (i)
	1133	for (j = 0; j < aiRozmiarCalychGrup[1]; j++)
	1134	{
[349]	1135
[869]	1136	// iterujemy j - Parts organizmu 1
	1137	// (indeks podstawowy - aiKoniecPierwszejGrupy[1])
[349]	1138
[869]	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]);
[349]	1145
[869]	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]);
[349]	1149
[869]	1150	iNeu = abs(m_aDegrees[0][aiKoniecPierwszejGrupy[0] + i][3]
	1151	- m_aDegrees[1][aiKoniecPierwszejGrupy[1] + j][3]);
[349]	1152
[869]	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
[349]	1162
[869]	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
[349]	1166
[869]	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	}
[606]	1178	}
	1179	}
	1180	}
[349]	1181
[606]	1182	// tutaj - sprawdzic tylko, czy tablica odleglosci lokalnych jest poprawnie obliczona
[349]	1183
[606]	1184	// WYKORZYSTANIE TABLICY ODLEGLOSCI DO BUDOWY DOPASOWANIA
[349]	1185
[606]	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
[349]	1195
[606]	1196	// zakładamy, że nie ma takiego Part
	1197	iIndex[i] = -1;
[349]	1198
[606]	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	}
[349]	1211
[606]	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	}
[349]	1222
	1223
[606]	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
[349]	1233
[606]	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	}
[349]	1245
[606]	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	{
[349]	1254
[606]	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	}
[349]	1261
[606]	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));
[349]	1268
[606]	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	}
[349]	1281
[606]	1282	// sprawdz poprawnosc znalezionego wczesniej minimum
	1283	assert(aadOdleglosciParts[iIndex[0]][i] >= dMinimum);
	1284	}
	1285	}
[349]	1286
[606]	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));
[349]	1306
[606]	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	}
[349]	1319
[606]	1320	// sprawdz poprawnosc znalezionego wczesniej minimum
	1321	assert(aadOdleglosciParts[i][iIndex[1]] >= dMinimum);
	1322	}
	1323	}
[349]	1324
[606]	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]);
[349]	1331
[606]	1332	if (PartZ1NaLiscie0 && PartZ0NaLiscie1)
	1333	{
	1334	// PRZYPADEK 1. Oba Parts maja sie wzajemnie na listach mozliwych
	1335	// dopasowan.
	1336	// AKCJA. Dopasowanie wzajemne do siebie.
[349]	1337
[606]	1338	m_pMatching->Match(0, aiKoniecPierwszejGrupy[0] + iIndex[0],
	1339	1, aiKoniecPierwszejGrupy[1] + iIndex[1]);
[349]	1340
[606]	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]);)
[349]	1347
[606]	1348	}// PRZYPADEK 1.
	1349	else
	1350	if (PartZ1NaLiscie0 \|\| PartZ0NaLiscie1)
	1351	{
[869]	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)
[349]	1357
[869]	1358	// indeks organizmu, ktorego Part nie ma dopasowywanego Part
	1359	// z drugiego organizmu na swojej liscie
	1360	int iBezDrugiego;
[349]	1361
[869]	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));
[349]	1373
[869]	1374	int iDopasowywany = -1;
	1375	// poszukujemy pierwszego z listy dopasowania
	1376	for (i = 0; i < aiRozmiarCalychGrup[1 - iBezDrugiego]; i++)
[606]	1377	{
[869]	1378	if (apvbCzyMinimalnaOdleglosc[iBezDrugiego]->operator[](i))
	1379	{
	1380	iDopasowywany = i;
	1381	break;
	1382	}
[606]	1383	}
[869]	1384	// sprawdz poprawnosc indeksu dopasowywanego (musimy cos znalezc!)
	1385	// nieujemny i w odpowiedniej grupie!
	1386	assert((iDopasowywany >= 0) &&
	1387	(iDopasowywany < aiRozmiarCalychGrup[1 - iBezDrugiego]));
[349]	1388
[869]	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);)
[349]	1397
[869]	1398	// zmniejsz liczby niedopasowanych elementow w grupach
	1399	aiRozmiarGrupy[0]--;
	1400	aiRozmiarGrupy[1]--;
[349]	1401
[869]	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
[349]	1408
[869]	1409	int iZDrugim = 1 - iBezDrugiego;
	1410	// sprawdz indeks organizmu
	1411	assert((iZDrugim == 0) \|\| (iZDrugim == 1));
[349]	1412
[869]	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	}
[349]	1425
[869]	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++)
[606]	1430	{
[869]	1431	if (!(m_pMatching->IsMatched(
	1432	1 - iZDrugim,
	1433	aiKoniecPierwszejGrupy[1 - iZDrugim] + i)))
[606]	1434	{
[869]	1435	// szukamy minimum obliczonej lokalnej odleglosci tylko wsrod
	1436	// niedopasowanych jeszcze Parts
	1437	if (iZDrugim == 0)
[606]	1438	{
[869]	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
[606]	1448	}
[869]	1449	else
[606]	1450	{
[869]	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);
[606]	1457	}
	1458	}
	1459	}
[869]	1460	// sprawdz, czy minimum znaleziono - musi takie byc, bo jest cos niedopasowanego
	1461	assert((dMinimum >= 0.0) && (dMinimum < HUGE_VAL));
[349]	1462
[869]	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++)
[606]	1466	{
[869]	1467	if (!(m_pMatching->IsMatched(
	1468	1 - iZDrugim,
	1469	aiKoniecPierwszejGrupy[1 - iZDrugim] + i)))
[606]	1470	{
[869]	1471	if (iZDrugim == 0)
[606]	1472	{
[869]	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	}
[606]	1481	}
[869]	1482	else
[606]	1483	{
[869]	1484	if (aadOdleglosciParts[i][iIndex[1]] == dMinimum)
	1485	{
	1486	apvbCzyMinimalnaOdleglosc[1]->operator[](i) = true;
	1487	}
[606]	1488	}
	1489	}
	1490	}
[349]	1491
[869]	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++)
[606]	1496	{
[869]	1497	if (apvbCzyMinimalnaOdleglosc[iZDrugim]->operator[](i))
	1498	{
	1499	iDopasowywany = i;
	1500	break;
	1501	}
[606]	1502	}
[869]	1503	// musielismy znalezc jakiegos dopasowywanego
	1504	assert((iDopasowywany >= 0) &&
	1505	(iDopasowywany < aiRozmiarCalychGrup[1 - iZDrugim]));
[349]	1506
[869]	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);)
[349]	1514
[869]	1515	//aiKoniecPierwszejGrupy[ 1-iBezDrugiego ] + iDopasowywany ;)
	1516	//aiKoniecPierwszejGrupy[ 1-iBezDrugiego ] + iDopasowywany ;)
	1517	// zmniejsz liczby niedopasowanych elementow w grupach
	1518	aiRozmiarGrupy[0]--;
	1519	aiRozmiarGrupy[1]--;
[349]	1520
[869]	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	}
[349]	1528
[869]	1529	DB(printf("Przypadek 2.\n");)
[349]	1530
[606]	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
[349]	1536
[606]	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]));
[349]	1551
[606]	1552	// teraz wlasnie dopasowujemy
	1553	m_pMatching->Match(
	1554	0,
	1555	aiKoniecPierwszejGrupy[0] + iIndex[0],
	1556	1,
	1557	aiKoniecPierwszejGrupy[1] + iDopasowywany);
[349]	1558
[606]	1559	// zmniejszamy liczbe niedopasowanych Parts
	1560	aiRozmiarGrupy[0]--;
	1561	aiRozmiarGrupy[1]--;
[349]	1562
[606]	1563	// debug - dopasowanie
	1564	DB(printf("Przypadek 3.: dopasowane Parts: (%2i, %2i)\n", aiKoniecPierwszejGrupy[0]
	1565	+ iIndex[0], aiKoniecPierwszejGrupy[1] + iDopasowywany);)
[349]	1566
[606]	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]));
[349]	1581
[606]	1582	// no i teraz realizujemy dopasowanie
	1583	m_pMatching->Match(
	1584	0,
	1585	aiKoniecPierwszejGrupy[0] + iDopasowywany,
	1586	1,
	1587	aiKoniecPierwszejGrupy[1] + iIndex[1]);
[349]	1588
[606]	1589	// zmniejszamy liczbe niedopasowanych Parts
	1590	aiRozmiarGrupy[0]--;
	1591	aiRozmiarGrupy[1]--;
[349]	1592
[606]	1593	// debug - dopasowanie
	1594	DB(printf("Przypadek 3.: dopasowane Parts: (%2i, %2i)\n", aiKoniecPierwszejGrupy[0]
	1595	+ iDopasowywany, aiKoniecPierwszejGrupy[1] + iIndex[1]);)
[349]	1596
	1597
[606]	1598	} // PRZYPADEK 3.
[349]	1599
[606]	1600	}// if (! bCzyKoniecGrupy)
	1601	else
	1602	{
[647]	1603	// gdy mamy juz koniec grup - musimy zlikwidowac tablice aadOdleglosciParts
[606]	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);
[349]	1610
[606]	1611	// musimy tez usunac tablice (wektory) mozliwosci dopasowania
	1612	SAFEDELETE(apvbCzyMinimalnaOdleglosc[0]);
	1613	SAFEDELETE(apvbCzyMinimalnaOdleglosc[1]);
	1614	}
	1615	} // while (! bCzyKoniecGrupy)
[349]	1616
[647]	1617	// PO DOPASOWANIU WSZYSTKIEGO Z GRUP (CO NAJMNIEJ JEDNEJ GRUPY W CALOSCI)
[349]	1618
[606]	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	}
[349]	1628
[606]	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
[349]	1637	}
	1638
	1639	/** Matches Parts in both organisms so that computation of similarity is possible.
[606]	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	*/
[349]	1649	int ModelSimil::MatchPartsGeometry()
	1650	{
[606]	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());
[349]	1654
[606]	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;
[349]	1665
	1666
[606]	1667	if (m_adFactors[3] > 0)
	1668	{
	1669	if (!ComputePartsPositionsBySVD())
	1670	{
	1671	return 0;
	1672	}
	1673	}
[349]	1674
[606]	1675	DB(printf("Przed sortowaniem:\n");)
	1676	DB(_PrintDegrees(0);)
	1677	DB(_PrintDegrees(1);)
[349]	1678
[606]	1679	if (!SortPartInfoTables())
	1680	return 0;
[349]	1681
[606]	1682	DB(printf("Po sortowaniu:\n");)
	1683	DB(_PrintDegrees(0);)
	1684	DB(_PrintDegrees(1);)
[349]	1685
[606]	1686	if (m_adFactors[3] > 0)
	1687	{
[869]	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 };
[349]	1698
[361]	1699	#ifdef max
[606]	1700	#undef max //this macro would conflict with line below
[361]	1701	#endif
[869]	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
[349]	1705
[869]	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
[606]	1714	for (iPart = 0; iPart < m_Mod[0]->getPartCount(); iPart++)
	1715	{
[869]	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;
[606]	1719	}
[869]	1720	// now the original positions of model 0 are stored
[349]	1721
	1722
[869]	1723	int iTransform; // a counter of geometric transformations
	1724	for (iTransform = 0; iTransform < NO_OF_TRANSFORM; iTransform++)
[606]	1725	{
[869]	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)
[349]	1731
[869]	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();
[349]	1743
[869]	1744	// teraz m_pMatching istnieje i jest pełne
	1745	assert(m_pMatching != NULL);
	1746	assert(m_pMatching->IsFull() == true);
[349]	1747
[869]	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);
[349]	1758
[869]	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++)
[606]	1795	{
[869]	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	}
[606]	1804	}
[349]	1805
[606]	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);
[349]	1814
[606]	1815	// _PrintDegrees(0);
	1816	// _PrintDegrees(1);
[349]	1817
[606]	1818	DB(_PrintPartsMatching();)
[349]	1819
[606]	1820	return 1;
[349]	1821	}
	1822
	1823	void ModelSimil::_PrintSeamnessTable(std::vector<int> *pTable, int iCount)
	1824	{
[606]	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	{
[349]	1832
[606]	1833	printf("%3i ", pTable->operator[](i));
	1834	}
	1835	printf("\n");
[349]	1836	}
	1837
	1838	/** Computes elements of similarity (m_iDD, m_iDN, m_dDG) based on underlying matching.
[606]	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	*/
[349]	1846	int ModelSimil::CountPartsDistance()
	1847	{
[606]	1848	int i, temp;
[349]	1849
[606]	1850	// assume existence of m_pMatching
	1851	assert(m_pMatching != NULL);
	1852	// musi byc pelne!
	1853	assert(m_pMatching->IsFull() == true);
[349]	1854
[606]	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;
[349]	1861
[606]	1862	int iOrgPart, iOrgMatchedPart; // orginalny indeks Part i jego dopasowanego Part
	1863	int iMatchedPart; // indeks (wg sortowania) dopasowanego Part
[349]	1864
[606]	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	}
[349]	1907
[606]	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));)
[349]	1916
[606]	1917	return 1;
[349]	1918	}
	1919
	1920	/** Computes new positions of Parts of both of organisms stored in the object.
[606]	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	*/
[349]	1927	bool ModelSimil::ComputePartsPositionsBySVD()
	1928	{
[606]	1929	bool bResult = true; // the result; assume a success
[349]	1930
[606]	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);
[349]	1938
[606]	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();
[349]	1947
[869]	1948	double pDistances = new double[nSize nSize];
[349]	1949
[606]	1950	for (int i = 0; i < nSize; i++)
	1951	{
	1952	pDistances[i] = 0;
	1953	}
[349]	1954
[606]	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
[869]	1960
[817]	1961	double *weights = new double[nSize];
	1962	for (int i = 0; i < nSize; i++)
	1963	{
[869]	1964	if (wMDS == 1)
[817]	1965	weights[i] = 0;
	1966	else
	1967	weights[i] = 1;
	1968	}
[869]	1969
	1970	if (wMDS == 1)
[817]	1971	for (int i = 0; i < pModel->getJointCount(); i++)
	1972	{
	1973	weights[pModel->getJoint(i)->p1_refno]++;
[869]	1974	weights[pModel->getJoint(i)->p2_refno]++;
[817]	1975	}
[869]	1976
[606]	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;
[605]	1983	if (fixedZaxis == 1)
[606]	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;
[605]	1993	if (fixedZaxis == 1)
[606]	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)
[349]	2003
[817]	2004	MatrixTools::weightedMDS(vEigenvalues, nSize, pDistances, m_aPositions[iMod], weights);
[605]	2005	if (fixedZaxis == 1) //restore the original vertical coordinate of each Part
[601]	2006	{
[606]	2007	for (int part = 0; part < pModel->getPartCount(); part++)
	2008	{
	2009	m_aPositions[iMod][part].z = pModel->getPart(part)->p.z;
	2010	}
[601]	2011	}
[606]	2012
[817]	2013	delete[] pDistances;
	2014	delete[] weights;
[601]	2015	}
[349]	2016
[606]	2017	return bResult;
[349]	2018	}
	2019
[869]	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
[349]	2032	void ModelSimil::p_evaldistance(ExtValue args, ExtValue ret)
	2033	{
[606]	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));
[356]	2040	}
[869]	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	{
[871]	2080	logPrintf("ModelSimil", "EvaluateDistanceHungarian", LOG_ERROR, "NULL genotype pointer(s)");
[872]	2081	return 0.0;
[869]	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	{
[871]	2090	logPrintf("ModelSimil", "EvaluateDistanceHungarian", LOG_ERROR, "NULL or invalid model pointer(s)");
[872]	2091	return 0.0;
[869]	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;
[872]	2167	if (saveMatching)
[869]	2168	{
	2169	minAssignment.clear();
	2170	minAssignment.insert(minAssignment.begin(), assignment, assignment + nBigger);
	2171	}
	2172	}
	2173	}
	2174
	2175	dResult = dMinSimValue;
[872]	2176	if (saveMatching)
[869]	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: