// COMPILE_LINE: g++ -o dali22_code dali22_code.cpp -lm -lstdc++ -lgomp -fopenmp


// Supplementary material code for the paper Soler, Molazem, Subr'2022, 
// A Theoretical Analysis of Compactness of the Light Transport Operator
//
// Copyright (C) 2022  Cyril Soler (cyril.soler@inria.fr)
// 
// This program is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
// 
// 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 General Public License
// along with this program.  If not, see <https://www.gnu.org/licenses/>.

#include <cmath>
#include <iomanip>
#include <iostream>
#include <fstream>
#include <vector>
#include <limits.h>

void test_compactness(double p_distance);

int main(int,char *[])
{
	double p = 2.0;

	test_compactness(p);

	return 0;
}

// Computes the form factor of patch to a disk. The disk is tangent to the support plane of
// the patch and orthogonal to it. Parameters are:
//      r  : radius of the disk
//      x,y: coordinates of the patch in its plane, assuming the tangent point of the disk is (0,0),
//			 x is the distance orthogonal to the axis of the disk
//			 y is the distance along the axis of the disk
//
// We use the formula of Eq.2 in 
//     NARAGHI, M. H. N. (1988). 
//     Radiation view factors from differential plane sources to disks - A general formulation. 
//     Journal of Thermophysics and Heat Transfer,  2(3), 271–274. doi:10.2514/3.96
//
// using:
// 		theta=pi/2, h=y, b=x, a=r, phi=pi/2

double point_to_tangent_disk_form_factor(double x,double y,double r)
{
	double alpha = 2*r*r+x*x+y*y;

	return 0.5*r*y/(r*r+x*x) * (alpha/sqrt(alpha*alpha - 4*r*r*(r*r+x*x))-1);
}

// Monte-Calo integration version of the previous function, optionaly used for testing.
//
double point_to_tangent_disk_form_factor_MC(double x,double y,double r)
{
	// uniform sample the disk
	uint32_t NR = 100;
	uint32_t NT = 300;

	double z = 0.0;
	double res = 0.0;

	for(uint32_t i=0;i<NR;++i)
	{
		double R = r*(i+0.5)/NR;

		for(uint32_t j=0;j<NT;++j)
		{
			double theta = (j+0.5)*2.0*M_PI/NT;

			double cx = R*sin(theta);
			double cy = 0.0;
			double cz = r + R*cos(theta);

			double d2 = (cx-x)*(cx-x)+(cy-y)*(cy-y)+(cz-z)*(cz-z);
			double d  = sqrt(d2);

			double ctp_times_d = y;
			double ct_times_d  = cz;

			res += ctp_times_d*ct_times_d / (d2*d2) * R;
		}
	}
	res *= 2*M_PI*r/(NR*NT) / M_PI;

	return res;
}

template<typename T> static std::string number(T x,int decimals=0,bool leading_zeros=false)
{
	std::ostringstream o ;

	if(leading_zeros)
	{
		o << std::setw(decimals);
		o << std::setiosflags(std::ios::fixed);
		o << std::setfill('0') ;
	}
	if(decimals > 0) o << std::setprecision(decimals) ;

	o << x ;
	return o.str() ;
}

bool test_point_to_tangent_disk_form_factor()
{
	double r=1;
	bool succeed = true;

	std::cerr << "Testing analytical point-to-disk form factor formula..." << std::endl;

	for(int i=0;i<10;++i)
	{
		double x = 10*(drand48()-0.5);
		double y = 10*drand48();

		double ff_analytical = point_to_tangent_disk_form_factor(x,y,r);
		double ff_test_MC    = point_to_tangent_disk_form_factor_MC(x,y,r);

		double error = 100.0*fabs(ff_analytical-ff_test_MC)/std::max(ff_test_MC,ff_analytical);
		succeed = succeed && (error<0.01);

		std::cerr  << "   x=" << number(x,8,true) << ",\t y=" << number(y,8,true) << "\t FF_analytical=" << number(ff_analytical,8,true) << ", FF_MC=" << number(ff_test_MC,8,true) << " error=" << number(error,8,true) << "% : " << ((error<0.01)?"OK":"FAILED") << std::endl;
	}

	return succeed;
}

void compute_pairwise_distances(double r1,double r2,int N,double& Lp_dist,double p,double& Lp_v)
{
	double res1 = 0.0;
	double res3 = 0.0;

	double S_n   = M_PI*r1*r1;
	double S_np1 = M_PI*r2*r2;

	// make || L_n ||_p = 1
	double L_n   = 1.0/pow(M_PI*S_n  ,1.0/p);
	double L_np1 = 1.0/pow(M_PI*S_np1,1.0/p);

	for(int i=0;i<N;++i)
	{
		double x = (i+0.5)/(double)N - 0.5;// [-0.5,0.5]

		for(int j=0;j<N;++j)
		{
			double z = (j+1)/(double)N; // [0,1]

			double F_n   = point_to_tangent_disk_form_factor(x,z,r1);
			double F_np1 = point_to_tangent_disk_form_factor(x,z,r2);

			double TL_n   = L_n*F_n/M_PI;
			double TL_np1 = L_np1*F_np1/M_PI;

			res1 += pow(TL_n,p);
			res3 += pow(fabs(TL_n - TL_np1),p);
		}
	}
	Lp_dist = pow(res3/(N*N),1.0/p);
	Lp_v = pow(res1,1.0/p);
}

std::ostream& print(const std::vector<double>& a,int N,int P,std::ostream& o,int max_size=0)
{
	if(max_size == 0)
		max_size = INT_MAX;

	for(size_t i=0;i < max_size && i<(size_t)N;++i)
	{
		o << "[ " ;
		for(size_t j=0;j < max_size && j<(size_t)P;++j)
		{
			o.width(12) ;
			o << a[i+N*j] ;
			o.width(1) ;
			o << "," ;
		}
		if(max_size < P)
			o << "  ..." ;

		o << " ]" << std::endl;
	}
	if(max_size < N)
		o << "..." << std::endl;

	return o ;
}

void test_compactness(double p_distance)
{
	// This test considers a scene made of two squares connected by an edge, making an angle of pi/2
	// We form a bounded sequence of functions {u_n}_{n>0} defined by
	//
	//         u_n (x) = (\pi r^2)^{-1/p} delta(C_n)(x)
	//
	//  ...where C_n is the circle of radius r=1/1.2^n included in the support plan of one of the two squares,
	//  and tangent to the common edge at the middle point. For all n, we have ||u_n||_p = 1, so the series is bounded.
	//
	// We compute the illumination function on the other square (i.e. T u_n) using explicit form-factor formulas
	// and verify whether or not the series is Cauchy, meaning that the distance between e.g. consecutive values converges to
	// 0. If not, the series is not convergent since the space of square-integrable functions over the scene is Hilbert and therefore
	// complete (meaning all Cauchy series are convergent).

	test_point_to_tangent_disk_form_factor();

	// First, draw the shape of the transported radiant exitance on the receiver (Fig.3, middle)

	std::cerr << std::endl;
	std::cerr << "Generating data for Figure 3.b to file fig3b_data.gnu" << std::endl;
	std::cerr << "Plot with gnuplot, using:  splot 'fig3b_data.gnu' matrix w li" << std::endl;

	std::ofstream of("fig3b_data.gnu");

	const int N = 80; // number of points on the receiving surface
	double r_n=1.0/pow(1.2,7);
	double S_n=M_PI*r_n*r_n;

	for(int i=0;i<N;++i)
	{
		double x = 5*((i+0.5)/N - 0.5);

		for(int j=0;j<N;++j)
		{
			double y = 5*(j+0.5)/N;
			double tbn = point_to_tangent_disk_form_factor(x,y,r_n);

			double L_n   = 1.0/pow(S_n,0.5);

			of << L_n*tbn << " ";
		}
		of << std::endl;
	}
	of.close();
	std::cerr << "Done."<< std::endl;
	std::cerr << std::endl;

	// Then, compute pairwise distances between T_b B_n and T_b B_{n+p}, starting from n=2 (arbitrary choice)

	std::cerr << "Computing pairwise distances between T_b B_n and T_b B__{n+p} for n=2..20 and p=1..20, for Fig3c." << std::endl;
	uint32_t n_layers = 20;

	std::vector<double> Lp_distances(n_layers*n_layers,0.0) ;
	std::vector<double> Lp_v(n_layers,0.0);

	for(uint32_t n=2;n<n_layers;n+=1)
	{
		double r_n = pow(1.2,-(int)n);

#pragma omp parallel for
		for(uint32_t p=1;p<n_layers;p+=1)
		{
			double r_np = pow(1.2,-(int)(n+p));

			compute_pairwise_distances(r_n,r_np,500,Lp_distances[n+n_layers*p],2.0,Lp_v[n]);
		}
		std::cerr << "n=" << n << "\r" ;
		std::cerr.flush();
	}

	std::cerr << "Pairwise distances matrix (Lp distance with p=" << p_distance << "):" << std::endl;
	print(Lp_distances,n_layers,n_layers,std::cerr,15);
}