import java.io.*;
/**
* Created by IDEA.
* User: Safroshkin
* Date: 18.11.2005
* Time: 10:20:17
* To change this template use File | Settings | File Templates.
*/
// набор утилит для сильного сглаживания на равномерной сетке по "Х"
// степень возможного сглаживания существенно выше, чем при использовании кубических сплайнов
public class ParabolicSplineSmoothing {
public static String read (String Nameinp) throws IOException {
StringBuffer sb = new StringBuffer();
String s;
BufferedReader in = new BufferedReader((new FileReader(Nameinp)));
while ((s = in.readLine()) != null ) {
sb.append(s);
sb.append("\n");
}
in.close();
return sb.toString();
}
public static void chobandet(double d1, double a[], double l[], int n, int m)
{
int d2=0;
int i, j, k, p, q, r, s;
double y;
d1=1.0;
for(i=1; i<=n; i++)
{
if (i>m) p=0;
else p=m-i+1;
r=i-m+p;
for (j=p; j<=m; j++)
{
s=j-1;
q=m-j+p;
y=a[2*i+j];
for (k=p; k<=s; k++)
{
y-=l[i*2+k]*l[r*2+q];
q++;
}
if(j==m)
{
d1*=y;
if (y==0) {d2=0; return;}
while(d1>=1.0)
{
d1/=16.0;
d2+=4;
}
while(Math.abs(d1)<0.0625)
{
d1*=16.0;
d2-=4;
}
if(y<0.0) return;
l[i*2+j]=1.0/Math.sqrt( y );
}
else
{
l[i*2+j]=y*l[r*2+m];
r++;
}
}
}
return;
}
public static void chobandsol(int n, int m, double l[], double b[], double x[])
{
int i, j, k, p, q, s;
double y;
s=m-1;
j=1;
for (i=1; i<=n; i++)
{
p=0;
if (i<=m) p=m-i+1;
q=i;
y=b[i];
for (k=s; k>=p; k--)
{
q--;
y-=l[i*2+k]*x[q];
}
x[i]=y*l[i*2+m];
}
for (i=n; i>0; i--)
{
p=0;
if (n-i<=m) p=m-n+i;
y=x[i];
q=i;
for (k=s; k>=p; k--)
{
q++;
y-=l[q*2+k]*x[q];
}
x[i]=y*l[i*2+m];
}
return;
}
static double kb(int m, int i, int j)
{
double d;
d=(2.0*i-2.0*m-1.0)*(2.0*j-2.0*m-1.0)+1.0/3.0;
return(d/4.0);
}
public static void assa(int nn, double al, double[] a) {
double[] h = new double[4];
double s, c, q;
int i1, i2, fs;
int i, j, n, k, l, m;
n=nn+1;
fs=1;
h[1]=0.5;
h[2]=-1.0;
h[3]=0.5;
for (i=1; i<=n-2; i++)
{
fs=n-2-i;
if(fs>1) fs=1;
for (k=0; k<=fs; k++)
{
s=0.0;
for (j=i+k; j<=i+1; j++)
for (l=j-i+2; l<=3; l++)
for (m=j-i-k+2; m<=3; m++)
s+=h[l] * h[m] * kb(j, i+l-1, i+k+m-1);
a[(i+k)*2+1-k]=s;
}
}
for (i=1; i<=n-2; i++)
{
a[i*2+1]+=al*0.5;
if(i>1)
a[i*2]-=al*0.25;
}
}
public static void rgt(double[] ay, double[] b, int N) {
double[] h = new double[] {0.0, -0.5, 0.5};
double s;
for (int i=1; i<=N-2; i++){
s=0.0;
for (int j=1; j<=2; j++) s+=h[j]*ay[i+j-2];
b[i]=s;
}
}
public static void main (String[] args) throws Exception {
String vasia = read("sample.dat");
String[] file_d = vasia.split("\n");
int nN = file_d.length;
double d1 = 0.0;
double[] x = new double[nN+1];
double[] y = new double[nN+1];
double[] b = new double[nN+1];
double[] af = new double[nN+1];
double[] l = new double[nN*2+2];
double[] a = new double[8191];
String[] sss = new String[2];
// nN - порядка 150
for (int i =0; i < nN; i++) {
sss = file_d[i].split(" ");
x[i] = Double.valueOf(sss[0]).doubleValue();
y[i] = Double.valueOf(sss[1]).doubleValue();
}
// величина "al" определяет степень сглаживания
// для данных с погрешностью свыше 5 %, al нужно брать в диапазоне от 0.2 до 1.0
// для данных с погрешностью от 0.2 до 1 %, al нужно брать в диапазоне от 0.02 до 0.1
double al = 0.3;
assa(nN, al, a);
chobandet(d1, a, l, nN-2, 1);
rgt(y, b, nN);
nN++;
nN--;
chobandsol(nN-2, 1, l, b, af);
double[] h = new double[] {0.0, -0.5, 0.5};
double c = 0.0;
for (int i=1; i<=nN-1; i++)
{
double s=0.0;
for (int j=2; j>0; j--)
{
int m=i-j+1;
if( (1<=m) && (m<=nN-2) ) c=h[j]*af[m];
else c=0.0;
s+=c;
}
y[i-1]-=al*s;
}
nN++;
nN--;
File sample = new File("sample.djj");
sample.delete();
sample.createNewFile();
BufferedWriter out = new BufferedWriter(new OutputStreamWriter(new FileOutputStream("sample.djj", true), "Cp1251"));
try {
for (int i=0; i<nN; i++) {
out.write(Double.toString(x[i]) + " "+ Double.toString(y[i]) + "\n");
}
} finally {
out.close();
}
}
}
