Index: /cpp/frams/model/similarity/measure-distribution.cpp
===================================================================
--- /cpp/frams/model/similarity/measure-distribution.cpp	(revision 1119)
+++ /cpp/frams/model/similarity/measure-distribution.cpp	(revision 1120)
@@ -1,4 +1,4 @@
 // This file is a part of Framsticks SDK.  http://www.framsticks.com/
-// Copyright (C) 1999-2020  Maciej Komosinski and Szymon Ulatowski.
+// Copyright (C) 1999-2021  Maciej Komosinski and Szymon Ulatowski.
 // See LICENSE.txt for details.
 
@@ -15,5 +15,5 @@
 		{ "simil_density", 0, 0, "Density of surface sampling", "f 1 100 10", FIELD(density), "", },
 		{ "simil_bin_num", 0, 0, "Number of bins", "d 1 1000 128", FIELD(bin_num), "", },
-		{ "simil_samples_num", 0, 0, "Number of samples", "d 1 1048576 1048576", FIELD(samples_num), "", },
+		{ "simil_samples_num", 0, 0, "Number of samples", "d 1 1048576 10000", FIELD(samples_num), "", }, //based on experiments, not much accuracy to gain when this is increased above 1000
 		{ "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.", },
 		{ 0, },
@@ -37,11 +37,11 @@
 			funs[i][j] = std::make_pair(0, 0);
 	}
-	
+
 	for (int i = 0; i < 2; i++)
 		sst_models[i] = new SolidsShapeTypeModel((*models[i]));
-	
-	SimilMeasureDistribution::calculateFuns();	
+
+	SimilMeasureDistribution::calculateFuns();
 	dist = SimilMeasureDistribution::compareFuns();
-	
+
 	for (int i = 0; i < 2; i++)
 	{
@@ -60,13 +60,13 @@
 			return -1;
 		}
-	
+
 	density = params.at(0);
 	bin_num = params.at(1);
 	samples_num = params.at(2);
-	
+
 	return 0;
 }
 
-void SimilMeasureDistribution::calculateFun(std::pair<double, float> *fun, Model &sampled)
+void SimilMeasureDistribution::calculateFun(std::pair<double, float> *fun, const Model &sampled)
 {
 	int size = sampled.getPartCount();
@@ -75,9 +75,14 @@
 	//Check if total number of points pairs is smaller than samples number
 	//but first, prevent exceeding int limits
-	if (size < (int) sqrt((double) std::numeric_limits<int>::max()))
-		samples_taken = min(samples_num, size*size);
+	//if (size < (int) sqrt((double) std::numeric_limits<int>::max()))
+	//	samples_taken = min(samples_num, size*size);
+
+	rndgen.seed(55); //For determinism. Otherwise the descriptors (that choose samples pseudo-randomly) for the same Model can yield different values and so the dissimilarity between the object and its copy will not be 0.
+	std::uniform_int_distribution<> uniform_distrib(0, sampled.getPartCount() - 1);
 
 	//Get sampled distribution
-	std::vector<double> dist_vect = (this->*SimilMeasureDistribution::distribution_fun)(samples_taken, &sampled);
+	std::vector<double> dist_vect;
+	dist_vect.reserve(samples_taken); //we will add up to samples_taken elements to this vector
+	(this->*SimilMeasureDistribution::distribution_fun)(samples_taken, uniform_distrib, sampled, dist_vect);
 
 	auto result = std::minmax_element(dist_vect.begin(), dist_vect.end());
@@ -86,14 +91,14 @@
 
 	//Create histogram
-	vector<float> hist(bin_num);
+	vector<int> hist(bin_num);
 	int ind = 0;
 
 	for (unsigned int j = 0; j < dist_vect.size(); j++)
 	{
-		ind = (int) std::floor((dist_vect.at(j)-min)*1/(max-min)*bin_num);
-		if (ind <= (bin_num-1))
-			hist[ind] += 1;
+		ind = (int)std::floor((dist_vect.at(j) - min) * 1 / (max - min) * bin_num);
+		if (ind <= (bin_num - 1))
+			hist[ind]++;
 		else if (ind == bin_num)
-			hist[bin_num-1] += 1;
+			hist[bin_num - 1]++;
 	}
 
@@ -101,5 +106,5 @@
 	for (int j = 0; j < bin_num; j++)
 	{
-		fun[j] = std::make_pair(min+(max-min)/bin_num*(j+0.5), hist[j]);
+		fun[j] = std::make_pair(min + (max - min) / bin_num * (j + 0.5), hist[j]);
 	}
 
@@ -107,7 +112,7 @@
 	float total_mass = 0;
 	for (int j = 0; j < bin_num; j++)
-    {
-        total_mass += fun[j].second;
-    }
+	{
+		total_mass += fun[j].second;
+	}
 
 	for (int j = 0; j < bin_num; j++)
@@ -129,134 +134,117 @@
 }
 
-vector<double> SimilMeasureDistribution::D1(int samples_taken, Model *sampled)
-{
-	vector<double> dist_vect;
-	int size = sampled->getPartCount();
+void SimilMeasureDistribution::D1(int samples_taken, std::uniform_int_distribution<> &distribution, const Model &sampled, std::vector<double> &dist_vect)
+{
+	int size = sampled.getPartCount();
 	double x = 0;
 	double y = 0;
 	double z = 0;
-	
+
 	for (int i = 0; i < size; i++)
 	{
-		Pt3D pos = sampled->getPart(i)->p;
+		Pt3D pos = sampled.getPart(i)->p;
 		x += pos.x;
 		y += pos.y;
 		z += pos.z;
 	}
-	
-	x = x/size;
-	y = y/size;
-	z = z/size;
-	
-	Pt3D centroid = {x, y, z};
-	
-	for (int i = 0; i < samples_taken; i++)
-	{
-		int p1 = rndUint(size);
-		double dist = sampled->getPart(p1)->p.distanceTo(centroid);
-		if (dist > 0)
-			dist_vect.push_back(dist);
-	}
-	
-	return dist_vect;
-}
-
-vector<double> SimilMeasureDistribution::D2(int samples_taken, Model *sampled)
-{
-	vector<double> dist_vect;
-	int size = sampled->getPartCount();
-	for (int i = 0; i < samples_taken; i++)
-	{
-		int p1 = rndUint(size);
-		int p2 = rndUint(size);
-		double dist = sampled->getPart(p1)->p.distanceTo(sampled->getPart(p2)->p);
-		if (dist > 0)
-			dist_vect.push_back(dist);
-	}
-	
-	return dist_vect;
-}
-
-vector<double> SimilMeasureDistribution::D3(int samples_taken, Model *sampled)
-{
-	vector<double> dist_vect;
-	int size = sampled->getPartCount();
-	for (int i = 0; i < samples_taken; i++)
-	{
-		int p1 = rndUint(size);
-		int p2 = rndUint(size);
-		int p3 = rndUint(size);
-		
-		Pt3D v(sampled->getPart(p2)->p);
-		Pt3D w(sampled->getPart(p3)->p);
-		v -= sampled->getPart(p1)->p;
-		w -= sampled->getPart(p1)->p;
+
+	x = x / size;
+	y = y / size;
+	z = z / size;
+
+	Pt3D centroid = { x, y, z };
+
+	for (int i = 0; i < samples_taken; i++)
+	{
+		int p1 = distribution(rndgen);
+		double dist = sampled.getPart(p1)->p.distanceTo(centroid);
+		if (dist > 0)
+			dist_vect.push_back(dist);
+	}
+}
+
+void SimilMeasureDistribution::D2(int samples_taken, std::uniform_int_distribution<> &distribution, const Model &sampled, std::vector<double> &dist_vect)
+{
+	int size = sampled.getPartCount();
+	for (int i = 0; i < samples_taken; i++)
+	{
+		int p1 = distribution(rndgen);
+		int p2 = distribution(rndgen);
+		double dist = sampled.getPart(p1)->p.distanceTo(sampled.getPart(p2)->p);
+		if (dist > 0)
+			dist_vect.push_back(dist);
+	}
+}
+
+void SimilMeasureDistribution::D3(int samples_taken, std::uniform_int_distribution<> &distribution, const Model &sampled, std::vector<double> &dist_vect)
+{
+	for (int i = 0; i < samples_taken; i++)
+	{
+		int p1 = distribution(rndgen);
+		int p2 = distribution(rndgen);
+		int p3 = distribution(rndgen);
+
+		Pt3D v(sampled.getPart(p2)->p);
+		Pt3D w(sampled.getPart(p3)->p);
+		v -= sampled.getPart(p1)->p;
+		w -= sampled.getPart(p1)->p;
 		Pt3D cross_prod(0);
 		cross_prod.vectorProduct(v, w);
-		
+
 		double dist = 0.5 * cross_prod.length();
 		if (dist > 0)
 			dist_vect.push_back(dist);
 	}
-	
-	return dist_vect;
-}
-
-vector<double> SimilMeasureDistribution::D4(int samples_taken, Model *sampled)
-{
-	vector<double> dist_vect;
-	int size = sampled->getPartCount();
-	for (int i = 0; i < samples_taken; i++)
-	{
-		int a = rndUint(size);
-		int b = rndUint(size);
-		int c = rndUint(size);
-		int d = rndUint(size);
-		
-		Pt3D ad(sampled->getPart(a)->p);
-		Pt3D bd(sampled->getPart(b)->p);
-		Pt3D cd(sampled->getPart(c)->p);
-		
-		ad -= sampled->getPart(d)->p;
-		bd -= sampled->getPart(d)->p;
-		cd -= sampled->getPart(d)->p;
-		
+}
+
+void SimilMeasureDistribution::D4(int samples_taken, std::uniform_int_distribution<> &distribution, const Model &sampled, std::vector<double> &dist_vect)
+{
+	for (int i = 0; i < samples_taken; i++)
+	{
+		int a = distribution(rndgen);
+		int b = distribution(rndgen);
+		int c = distribution(rndgen);
+		int d = distribution(rndgen);
+
+		Pt3D ad(sampled.getPart(a)->p);
+		Pt3D bd(sampled.getPart(b)->p);
+		Pt3D cd(sampled.getPart(c)->p);
+
+		ad -= sampled.getPart(d)->p;
+		bd -= sampled.getPart(d)->p;
+		cd -= sampled.getPart(d)->p;
+
 		Pt3D cross_prod(0);
 		cross_prod.vectorProduct(bd, cd);
 		cross_prod.entrywiseProduct(ad);
-		
-		double dist = cross_prod.length()/6;
-		if (dist > 0)
-			dist_vect.push_back(dist);
-	}
-	
-	return dist_vect;
-}
-
-vector<double> SimilMeasureDistribution::A3(int samples_taken, Model *sampled)
-{
-	vector<double> dist_vect;
-	int size = sampled->getPartCount();
-	for (int i = 0; i < samples_taken; i++)
-	{
-		int p1 = rndUint(size);
-		int p2 = rndUint(size);
-		int p3 = rndUint(size);
-		double a = sampled->getPart(p1)->p.distanceTo(sampled->getPart(p3)->p);
-		double b = sampled->getPart(p3)->p.distanceTo(sampled->getPart(p2)->p);
-		double c = sampled->getPart(p1)->p.distanceTo(sampled->getPart(p2)->p);
-		double beta = acos((a*a + b*b - c*c)/(2*a*b));
-	
+
+		double dist = cross_prod.length() / 6;
+		if (dist > 0)
+			dist_vect.push_back(dist);
+	}
+}
+
+void SimilMeasureDistribution::A3(int samples_taken, std::uniform_int_distribution<> &distribution, const Model &sampled, std::vector<double> &dist_vect)
+{
+	int size = sampled.getPartCount();
+	for (int i = 0; i < samples_taken; i++)
+	{
+		int p1 = distribution(rndgen);
+		int p2 = distribution(rndgen);
+		int p3 = distribution(rndgen);
+		double a = sampled.getPart(p1)->p.distanceTo(sampled.getPart(p3)->p);
+		double b = sampled.getPart(p3)->p.distanceTo(sampled.getPart(p2)->p);
+		double c = sampled.getPart(p1)->p.distanceTo(sampled.getPart(p2)->p);
+		double beta = acos((a * a + b * b - c * c) / (2 * a * b));
+
 		if (!std::isnan(beta))
 			dist_vect.push_back(beta);
 	}
-	
-	return dist_vect;
 }
 
 
 float dist(feature_t* F1, feature_t* F2)
-{ 
-	return abs((*F1)-(*F2)); 
+{
+	return abs((*F1) - (*F2));
 }
 
@@ -266,5 +254,5 @@
 	for (int j = 0; j < bin_num; j++)
 	{
-		points[j] = {fun[j].first};
+		points[j] = { fun[j].first };
 		weights[j] = fun[j].second;
 	}
@@ -275,5 +263,5 @@
 	feature_t *points[2];
 	float *weights[2];
-	
+
 	for (int i = 0; i < 2; i++)
 	{
@@ -283,10 +271,10 @@
 	SimilMeasureDistribution::fillPointsWeights(fun1, points[0], weights[0]);
 	SimilMeasureDistribution::fillPointsWeights(fun2, points[1], weights[1]);
-	
-	signature_t sig1 = {bin_num, points[0], weights[0]},
-	sig2 = {bin_num, points[1], weights[1]};
+
+	signature_t sig1 = { bin_num, points[0], weights[0] },
+		sig2 = { bin_num, points[1], weights[1] };
 
 	float e = emd(&sig1, &sig2, dist, 0, 0, bin_num, bin_num);
-	
+
 	for (int i = 0; i < 2; i++)
 	{
Index: /cpp/frams/model/similarity/measure-distribution.h
===================================================================
--- /cpp/frams/model/similarity/measure-distribution.h	(revision 1119)
+++ /cpp/frams/model/similarity/measure-distribution.h	(revision 1120)
@@ -1,4 +1,4 @@
 // This file is a part of Framsticks SDK.  http://www.framsticks.com/
-// Copyright (C) 1999-2020  Maciej Komosinski and Szymon Ulatowski.
+// Copyright (C) 1999-2021  Maciej Komosinski and Szymon Ulatowski.
 // See LICENSE.txt for details.
 
@@ -8,4 +8,5 @@
 #include "simil-measure-base.h"
 #include "EMD/emd.h"
+#include <random>
 
 class SimilMeasureDistribution : public SimilMeasureBase
@@ -13,8 +14,8 @@
 public:
 	SimilMeasureDistribution();
-	~SimilMeasureDistribution(){};
+	~SimilMeasureDistribution() {};
 	int setParams(std::vector<double> params);
 	void calculateFuns();
-	void calculateFun(std::pair<double, float> *fun, Model &sampled);
+	void calculateFun(std::pair<double, float> *fun, const Model &sampled);
 	double compareFuns();
 	double EMD(std::pair<double, float> *fun1, std::pair<double, float> *fun2);
@@ -28,13 +29,16 @@
 	std::pair<double, float> *funs[2];
 	SolidsShapeTypeModel *sst_models[2];
-	std::vector<double> (SimilMeasureDistribution::*distribution_fun)(int samples_taken, Model *sampled);
+	void (SimilMeasureDistribution::*distribution_fun)(int samples_taken, std::uniform_int_distribution<> &distribution, const Model &sampled, std::vector<double> &dist_vect);
 
 	void fillPointsWeights(std::pair<double, float> *fun, feature_t *points, float *weights);
 	double getDistance();
-	std::vector<double> D2(int samples_taken, Model *sampled);
-	std::vector<double> D1(int samples_taken, Model *sampled);
-	std::vector<double> D3(int samples_taken, Model *sampled);
-	std::vector<double> D4(int samples_taken, Model *sampled);
-	std::vector<double> A3(int samples_taken, Model *sampled);
+	void D2(int samples_taken, std::uniform_int_distribution<> &distribution, const Model &sampled, std::vector<double> &dist_vect);
+	void D1(int samples_taken, std::uniform_int_distribution<> &distribution, const Model &sampled, std::vector<double> &dist_vect);
+	void D3(int samples_taken, std::uniform_int_distribution<> &distribution, const Model &sampled, std::vector<double> &dist_vect);
+	void D4(int samples_taken, std::uniform_int_distribution<> &distribution, const Model &sampled, std::vector<double> &dist_vect);
+	void A3(int samples_taken, std::uniform_int_distribution<> &distribution, const Model &sampled, std::vector<double> &dist_vect);
+
+private:
+	std::mt19937 rndgen;
 };
 
