/*
This is a draft program issued under the GPL version 3 or any later version
(https://gnu.org/licenses/gpl.html by John Harold Nixon)
written in D (*) for calculating the probability density estimate based on user-supplied
one-dimensional data so as to minimise the integrated mean square error (IMSE) using the second order diffusion
equation model with variable coefficients. It is based on the paper (http://www.bluesky-home.co.uk/Diffusion_density_paper_4.pdf)
* D (dlang.org) is a general purpose computer programming language that is basically a much improved version of C++ with
a similar syntax and in particular has bounds checking for arrays making it much easier to debug programs without the
need of valgrind or anything similar. It has many more features than C++, and I found it an easy transition to make from
developing in C++, though most of D's features I have never used.
*/
import std.stdio;
import std.math;
import std.algorithm.sorting;
import std.algorithm.comparison;
//import std.file;
//global variables
File file;
double fac1;
double xlo,xhi;//lower and upper limits of interval I of calculation
int n;//number of points in grid (= number of intervals +1).
double delta;//spacing of grid
int[] grid;//grid point nearest to data points
double K;//
int N;//number of data
double[] c,d,integrand,iden,inumd,f,fr,cnew,dnew,integ,term1,term2;
double[][] alpha,beta,G1,G2,dH1ds,H1,H2;
double[] pvi;//final principal value integral.
double[] num1,den1,int_1,int_2,res;
int cp;//index of the centre point for P.V. integral calculation
int j0;//Centre of data (approx. median) as a grid point index.
//******************************************** function pvi_centre *****************************************
//Calculate the principal value integral (res) of the function defined on the grid by num[]/denom[]
//by calculating first the integral of it from a central point in the grid to all other points.
//The central point (counting from zero as the initial point) is determined elsewhere.
void pvi_centre(const ref double num[],const ref double den[],const int centre){
//calculate pvi for bottom half of grid
for(int i=0;i<=centre;++i){
num1[i]=num[centre-i];
den1[i]=den[centre-i];}
for(int i=centre+1;i0?1:-1;
if(i && b_sign_prev*b_sign==-1){sg_change~=i;++k;}//record i where sign change occurs.
//in interval [xlo+(i-1)*delta, xlo+i*delta]
b_sign_prev=b_sign;}
//calculate zero points by linear interpolation of the grid
double[] z,az,bdash;
z.length=az.length=bdash.length=k;
double slope;
int j;
for(int s=0;sxlo+(j-0.5)*delta){
bdash[s]=0.5*(den[j+1]-den[j-1])/delta + (z[s]-(xlo+j*delta))*(den[j+1]-2*den[j]+den[j-1])/(delta*delta);}
else{
bdash[s]=0.5*(den[j]-den[j-2])/delta + (z[s]-xlo-(j-1)*delta)*(den[j]-2*den[j-1]+den[j-2])/(delta*delta);}}
for(int i=0;idata[i])datamin=data[i];
if(datamax