zipfUTest.cxxtest

/** zipfUTest.cxxtest ---
 *
 * Copyright (C) 2019, 2020 Linas Vepstas
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU Affero General Public License v3 as
 * published by the Free Software Foundation and including the exceptions
 * at http://opencog.org/wiki/Licenses
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU Affero General Public License
 * along with this program; if not, write to:
 * Free Software Foundation, Inc.,
 * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
 */

#include <climits>
#include <opencog/util/zipf.h>

using namespace opencog;

class zipfUTest : public CxxTest::TestSuite
{
	std::random_device rd;
	std::mt19937 gen;

public:
	// zipfUTest() : gen(rd()) {}
	zipfUTest() : gen(2) {}

	// Check that the generated distribution is normally distributed.
	void verify(const std::vector<size_t>& pdf, size_t ndraw,
	            int n, double s, double q, double sigma)
	{
		int len = n + 1;

		// Create the expected distribution
		double norm = 0.0;
		std::vector<double> expect;
		expect.resize(len);
		for (int i = len-1; 0<i; i--)
		{
			double term = std::pow (i+q, -s);
			expect[i] = term;
			norm += term;
		}

		// Normalize
		for (int i = 1; i <len; i++)
			expect[i] *= 1.0 / norm;

		// Also normalize the pdf.
		size_t cnt = 0;
		for (int i = 1; i <len; i++)
			cnt += pdf[i];
		TS_ASSERT_EQUALS (cnt, ndraw);

		double fcnt = (double) cnt;

		// Expected standard deviation.  1/(2pi sqrt(cnt))
		double erms = 1.0 / sqrt(2.0 * 3.14159 * fcnt);

		// Take the difference.  We expect the difference
		// to be zero. Also adjust the variance; we expect
		// a normal distribution.
		double mean = 0.0;
		double meansq = 0.0;
		for (int i = 1; i <len; i++)
		{
// #define GRAPH
#ifdef GRAPH
			printf("%d	%g	%g	%g\n", i, expect[i],
			      pdf[i]/fcnt, (expect[i] - pdf[i]/fcnt)/erms);
#endif

			expect[i] -= pdf[i] / fcnt;

			// The tail of the distribution is very accurate.
			// The head, not so much. Multiplying by sqrt seems
			// to turn the distribution into a normal distribution.
			// I don't know the theory; this seems to work, though.
			if (abs(q) < 0.1)
				expect[i] *= pow(i, 0.5*s);
			else
				expect[i] *= pow((i+q)/q, 0.5*s);

			expect[i] /= erms;

			mean += expect[i];
			meansq += expect[i] * expect[i];
		}
		mean *= 1.0 / n;
		meansq *= 1.0 / n;
		double rms = sqrt(meansq);

		printf("Dist: n=%d s=%4.3f q=%4.2f cnt=%zu "
		       "mean=%6.4g rms=%5.3g\n",
			n, s, q, cnt, mean, rms);

		TS_ASSERT_LESS_THAN(mean, 2.0);

		// Hmm. There is some non-trivial expression for nrms
		// that I do not know. Its close to 1.0 but seems to
		// depend on s and depend on .. what ? log(n) ???
		// TS_ASSERT_DELTA(nrms, 1.0, allow);

		// Actually test. We're going to use sigma as the max allowed
		// deviation. So six-sigma allows the unit tests to pass most
		// of the time.
		for (int i = 1; i < len; i++)
		{
			// Should have approx == 1.0, approximately.
			TS_ASSERT_DELTA(expect[i], 0.0, sigma);
		}

// #define PRINT_ERROR_DISTRIBUTION
#ifdef PRINT_ERROR_DISTRIBUTION
		// Dump the error distribution to a file, for graphing.
		std::vector<int> edf;
		edf.resize(200);
		for (int i = 0; i <200; i++)
			edf[i] = 0;

		for (int i = 1; i < len; i++)
		{
			int bin = (100.0 * expect[i] / 6.0) + 100;
			if (0>bin) bin=0;
			if (199 < bin) bin=199;
			edf[bin] += 1;
		}
		printf ("\n#The error distribution:\n#\n");
		for (int i = 0; i <200; i++)
			printf("%d	%f	%d\n", i, 6.0*(i-100.0)/100.0, edf[i]);
#endif // PRINT_ERROR_DISTRIBUTION
	}

	// Draw from the Zipf distribution per given parameters.
	void check_zipf(size_t ndraw, int n, double s, double q, double sigma)
	{
		zipf_distribution<> zipf(n, s, q);
		printf("Init rejection-inversion: n=%d s=%4.3f q=%4.2f draw=%zu\n",
			n, s, q, ndraw);

		std::vector<size_t> pdf;
		pdf.resize(n+1);
		for (int i = 0; i <= n; i++)
			pdf[i] = 0;

		for (size_t s = 0; s < ndraw; s++)
		{
			unsigned int draw = zipf(gen);
			TS_ASSERT_LESS_THAN_EQUALS(1, draw);
			TS_ASSERT_LESS_THAN_EQUALS(draw, n);
			pdf[draw] = pdf[draw] + 1;
		}

		verify(pdf, ndraw, n, s, q, sigma);
		printf("\n");
	}


	// Draw from the Zipf distribution per given parameters.
	void check_table(size_t ndraw, int n, double s, double q)
	{
		zipf_table_distribution<> zipf(n, s);
		printf("Init table: n=%d s=%4.3f q=%4.2f draw=%zu\n",
			n, s, q, ndraw);

		std::vector<size_t> pdf;
		pdf.resize(n+1);
		for (int i = 0; i <= n; i++)
			pdf[i] = 0;

		for (size_t s = 0; s < ndraw; s++)
		{
			unsigned int draw = zipf(gen);
			TS_ASSERT_LESS_THAN_EQUALS(1, draw);
			TS_ASSERT_LESS_THAN_EQUALS(draw, n);
			pdf[draw] ++;
		}
		verify(pdf, ndraw, n, s, q, 6.0);
		printf("\n");
	}

	// Test the Zipf distribution, 300 bins, for exponent s=1.
	// This is a pretty small test, as such things go.
	void test_zipf()
	{
		printf("\n");

		// Six sigma should do the trick.
		double sigma = 7.0;

		check_zipf(1623, 30, 1.0, 0.0, sigma);
		check_zipf(1623000, 300, 1.0, 0.0, sigma);
		check_zipf(1623000, 3000, 1.0, 0.0, sigma);
		check_zipf(5623000, 30000, 1.0, 0.0, sigma);

		check_zipf(1623000, 30, 0.2, 0.0, sigma);
		check_zipf(1623000, 30, 0.6, 0.0, sigma);
		check_zipf(1623000, 30, 0.8, 0.0, sigma);
		check_zipf(1623000, 30, 1.2, 0.0, sigma);
		check_zipf(1623000, 30, 1.5, 0.0, sigma);
		check_zipf(1623000, 30, 2.0, 0.0, sigma);
		check_zipf(1623000, 30, 3.1, 0.0, sigma);
		check_zipf(1623000, 30, 4.5, 0.0, sigma);
		check_zipf(1623000, 30, 8.2, 0.0, sigma);

		check_zipf(1623000, 30, 1.0 - 1e-3, 0.0, sigma);
		check_zipf(1623000, 30, 1.0 - 1e-6, 0.0, sigma);
		check_zipf(1623000, 30, 1.0 - 1e-9, 0.0, sigma);
		check_zipf(1623000, 30, 1.0 - 1e-12, 0.0, sigma);

		check_zipf(1623000, 30, 1.0 + 1e-3, 0.0, sigma);
		check_zipf(1623000, 30, 1.0 + 1e-6, 0.0, sigma);
		check_zipf(1623000, 30, 1.0 + 1e-9, 0.0, sigma);
		check_zipf(1623000, 30, 1.0 + 1e-12, 0.0, sigma);

		check_zipf(1623000, 3021, 0.01, 0.0, sigma);
		check_zipf(1623000, 3070, 0.53, 0.0, sigma);
		check_zipf(1623000, 2999, 1.03, 0.0, sigma);
		check_zipf(2623010, 213, 2.31, 0.0, sigma);

		check_zipf(1623000, 30, 1, -0.49, sigma);
		check_zipf(1623000, 30, 1, -0.4, sigma);
		check_zipf(1623000, 30, 1, -0.1, sigma);
		check_zipf(1623000, 30, 1,  0.1, sigma);

		check_zipf(1623000, 30, 1.1, -0.49, sigma);
		check_zipf(1623000, 30, 1.1, -0.4, sigma);
		check_zipf(1623000, 30, 1.1, -0.1, sigma);
		check_zipf(1623000, 30, 1.1,  0.1, sigma);

		// I don't understand the theory of the error
		// distribution for q != 0 but this seems to work.
		check_zipf(1623000, 30, 1.2, 0.4, 1.2*sigma);
		check_zipf(1623000, 30, 1.5, 0.4, 1.5*sigma);
		check_zipf(1623000, 30, 2.0, 0.4, 1.5*sigma);
		check_zipf(1623000, 30, 3.1, 0.4, 3.1*sigma);
		check_zipf(1623000, 30, 4.5, 0.4, 4.5*sigma);
		check_zipf(1623000, 30, 8.2, 0.4, 8.2*sigma);

		check_zipf(1623000, 30, 1.2, 4, sigma);
		check_zipf(1623000, 30, 1.5, 4, sigma);
		check_zipf(1623000, 30, 2.0, 4, sigma);
		check_zipf(1623000, 30, 3.1, 4, sigma);
		check_zipf(1623000, 30, 4.5, 4, 9.5);
		check_zipf(1623000, 30, 8.2, 4, 24.0);

		check_zipf(1623000, 30, 1.2, 44, sigma);
		check_zipf(1623000, 30, 1.5, 44, sigma);
		check_zipf(1623000, 30, 2.0, 44, sigma);
		check_zipf(1623000, 30, 3.1, 44, sigma);
		check_zipf(1623000, 30, 4.5, 44, sigma);
		check_zipf(1623000, 30, 8.2, 44, sigma);

		check_zipf(1623000, 30, 1.2, 44e4, sigma);
		check_zipf(1623000, 30, 1.5, 44e4, sigma);
		check_zipf(1623000, 30, 2.0, 44e4, sigma);
		check_zipf(1623000, 30, 3.1, 44e4, sigma);
		check_zipf(1623000, 30, 4.5, 44e4, sigma);
		check_zipf(1623000, 30, 8.2, 44e4, sigma);

		check_table(1623000, 310, 1.0, 0.0);
		check_table(1623000, 312, 1.2, 0.0);
		check_table(1623000, 217, 0.8, 0.0);
	}
};