Context Navigation

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

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

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

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

Download in other formats: