package jasext.hist.function;
import java.lang.Math.*;
import jas2.hist.*;
import java.util.*;
import jas2.util.*;
import java.awt.Frame;
import java.awt.Container;
import java.awt.GridLayout;
import java.awt.event.*;
import javax.swing.*;
import javax.swing.event.*;
public class SplineFunction extends Basic1DFunction implements FunctionAdvancedOptions
{
public String[] getParameterNames()
{
final String[] result = {"Smoothness"};
return result;
}
public double[] getParameterValues()
{
final double[] result = {smoothness};
return result;
}
public String getTitle()
{
return "Cubic Spline";
}
public void setParameter(int index, double value)
{
smoothness = value;
// never mind what index was supplied
setChanged();
}
public double valueAt(double x)
{
return spline.get(x);
}
private void init(Rebinnable1DHistogramData data)
{
smoothness = initialSmoothness;
int bins = data.getBins();
double min = data.getMin();
double max = data.getMax();
double bw = (max-min)/bins;
double[][] binnedData = data.rebin(bins,min,max,true,false);
double[] x = new double[bins];
double[] y = binnedData[0];
double[] yerr;
if (binnedData.length > 2 && binnedData[1] != binnedData[2])
throw new RuntimeException("Cannot create spline fit for data with asymetric error bars");
double xx = min + bw/2;
for (int i=0; i 1) yerr= binnedData[1];
else
{
yerr = new double[bins];
for (int i=0; i0) spline.add(x[i], y[i], yerr[i]);
//spline.calcInterpolation(true);
spline.calcApproximation(smoothness);
}
SplineFunction(Rebinnable1DHistogramData data)
{
init(data);
}
private synchronized void setSmoothness(double d)
{
smoothness = d;
spline.calcApproximation(smoothness);
setChanged();
}
public void openAdvancedDialog(Frame f, JASHist hist)
{
new SplineFunctionAdvancedDialog(f, hist);
}
private Spline spline;
private double smoothness;
final double initialSmoothness = 10.0;
private double getSmoothness()
{
return smoothness;
}
private class SplineFunctionAdvancedDialog extends JASDialog
{
private final java.text.NumberFormat format = java.text.NumberFormat.getInstance();
final int sliderMin = 10;
final int sliderMax = 1000;
final double maxSmoothness = 100.0;
// maximum allowed on slider
final double minSmoothness = 1.0;
// minimum allowed
double smoothness;
private int toInt(double smoothness)
{
double fractionIncrease = (smoothness - minSmoothness) /
(maxSmoothness - minSmoothness);
int increase = (int) ( fractionIncrease * (sliderMax-sliderMin) );
return increase + sliderMin;
}
SplineFunctionAdvancedDialog(Frame f, JASHist hist)
{
super(f, "Advanced...", true, JASDialog.OK_BUTTON);
format.setMaximumFractionDigits(5);
// This formatter is used for the number generated from
// the slider, which sometimes can generate very long
// decimal strings.
smoothness = getSmoothness();
final JSlider s = new
JSlider( JSlider.HORIZONTAL, sliderMin, sliderMax,
toInt(smoothness) );
final JASTextField tf = new JASTextField(String.valueOf(smoothness), 5);
tf.addActionListener(new ActionListener()
{
public void actionPerformed(ActionEvent te)
// the enter key was hit
{
try
{
smoothness = ( Double.valueOf(tf.getText()) ).doubleValue();
if (smoothness < 0.0) smoothness = 0.0;
double temp = smoothness;
s.setValue( toInt(smoothness) );
// causes an unwanted ChangeEvent
tf.setText(String.valueOf(temp));
// we hope that by now the ChangeEvent
// has finished
setSmoothness(smoothness = temp);
/*
* Without the temp variable, the method call
* s.setValue() will cause a ChangeEvent and
* set the smoothness to maxSmoothness.
* The variable temp simply stores the real
* smoothness, lets the ChangeEvent mess
* things up, and is used to reset things
* after the ChangeEvent on the slider.
*/
}
catch (NumberFormatException ex) {}
// if the number format is wrong, nothing will change
}
});
s.addChangeListener(new ChangeListener()
{
public void stateChanged(ChangeEvent e)
{
int increase = s.getValue() - sliderMin;
double fractionIncrease =
( (double) increase ) / ( (double) (sliderMax-sliderMin) );
smoothness = minSmoothness + (maxSmoothness - minSmoothness) *
fractionIncrease;
setSmoothness(smoothness);
tf.setText(format.format(smoothness));
}
});
Vector model = new Vector();
Enumeration e = hist.getDataSources();
while (e.hasMoreElements())
{
JASHistData d = (JASHistData) e.nextElement();
if (d.getDataSource() instanceof Rebinnable1DHistogramData) model.addElement(d);
}
JComboBox data = new JComboBox(model);
data.addItemListener(new ItemListener()
{
public void itemStateChanged(ItemEvent ev)
{
JASHistData d = (JASHistData) ev.getItem();
new SplineFunction((Rebinnable1DHistogramData) d.getFittableDataSource() );
setChanged();
}
});
Container c = getContentPane();
c.setLayout( new GridLayout(6, 1) );
c.add( new JLabel("Set spline smoothness.") );
c.add( new JLabel("Slider range is "+ minSmoothness +" to "+ maxSmoothness) );
c.add(s);
c.add(tf);
c.add( new JLabel("Select data source:") );
c.add(data);
pack();
show();
}
}
}
/** This class provides an interpolation/approximation facility to work with data.
* Partly this code is a modification of fortran subroutine smooth (file smooth.f)
* from pppack library at Netlib
* (editor: Eric Grosse,
* master: netlib.bell-labs.com)
* @author Alexander Sokolov,
* Department of Systems Ecology, Stockholm University
* @version July 29, 1997
*/
class Spline
{
/** number of data points */
public int lastElement;
/** array of data abscissae, assumed to be strictly increasing */
public double[] x;
/** array of data ordinates, corresponded to abscissae */
public double[] y;
/** array of smoothed ordinates, corresponded to abscissae */
public double[] f;
double[] d1f, d2f, d3f; // first three derivatives at the left end of the interval
double[] dy; // estimate of uncertainty in data (positive)
double[][] v;
double emptyCode;
/** Constructs the spline class and allocates memory.
* @param size is a max number of data points
* @param emptyCode is a value, which is smaller any of ordinates. (An ordinate
* smaller then emptyCode does not taking into account).
*/
public Spline(int size, double emptyCode)
{ x=new double[size];
y=new double[size];
dy=new double[size];
f=new double[size];
d1f=new double[size];
d2f=new double[size];
d3f=new double[size];
v=new double[size][7];
this.emptyCode=emptyCode;
lastElement=-1;
}
/** Adds a point into the arrays and sets uncertainty to zero.
* @param coord is a value of point's abscissa
* @param func is a value of point's ordinate.
*/
public void add(double coord, double func)
{ if (func > emptyCode)
{ lastElement+=1;
x[lastElement]=coord;
y[lastElement]=func;
dy[lastElement]=0.0;
}
}
/** Adds a point into the arrays.
* @param coord is a value of point's abscissa
* @param func is a value of point's ordinate.
* @param uncertainty is a estimation of uncertainty in data. Must be positive.
*/
public void add(double coord, double func, double uncertainty)
{ if (func > emptyCode)
{ add(coord, func);
dy[lastElement]=uncertainty;
}
}
/** Constructs the cubic smoothing spline to given data (x(),y()).
* This approximation has a second derivative as small as possible
* while sum(((y(i)-f(x(i)))/dy(i))^2) < maxMeanDistance
* @param maxMeanDistance is upper bound on the discrete weighted mean square distance of
* the approximation from the data
*/
public void calcApproximation(double maxMeanDistance)
{ if (lastElement > 0)
{ if (lastElement < 2)
{ for (int i=0; i < lastElement; i++)
{ f[i]=y[i];
d1f[i]=(y[i+1]-y[i])/(x[i+1]-x[i]);
d2f[i]=0.0;
d3f[i]=0.0;
}
f[lastElement]=y[lastElement];
}
else
{ int maxIter=100; // Usually it takes 3 to 5 iterations
double sfq, p=0.0, q=0.0;
SetupQ();
if (maxMeanDistance > 0.0)
{ p=0.0;
Chol1D(p);
sfq=0.0;
for (int i=0; i <= lastElement; i++)
sfq+=Math.pow(f[i]*dy[i], 2);
sfq=sfq*36.0;
if (sfq > maxMeanDistance)
{ double utru=0.0;
for (int i=1; i <= lastElement; i++)
utru+=v[i-1][3]*(d2f[i-1]*(d2f[i-1]+d2f[i])+Math.pow(d2f[i], 2));
double ooss=1.0/Math.sqrt(maxMeanDistance);
double oosf=1.0/Math.sqrt(sfq);
q=-(oosf-ooss)*sfq/(6.0*utru*oosf);
double prevq=0.0;
double prevsf=oosf;
boolean OK;
int iter=0;
do
{ Chol1D(q/(1.0+q));
iter+=1;
sfq=0.0;
for (int i=0; i <= lastElement; i++)
sfq+=Math.pow(f[i]*dy[i], 2);
sfq=sfq*36.0/Math.pow(1.0+q, 2);
if (!(OK=Math.abs(sfq-maxMeanDistance) <= 0.01*maxMeanDistance))
{ oosf=1.0/Math.sqrt(sfq);
double change=(q-prevq)/(oosf-prevsf)*(oosf-ooss);
prevq=q;
q-=change;
prevsf=oosf;
}
}
while (!OK && iter <= maxIter);
p=q/(1.0+q);
}
}
else
{ p=1.0;
Chol1D(p);
sfq=0.0;
}
double temp=6.0/(1.0+q);
for (int i=0; i <= lastElement; i++)
f[i]=y[i]-temp*Math.pow(dy[i], 2)*f[i];
for (int i=0; i <= lastElement; i++)
d2f[i]*=6.0*p;
for (int i=0; i < lastElement; i++)
{ d3f[i]=(d2f[i+1]-d2f[i])/v[i][3];
d1f[i]=(f[i+1]-f[i])/v[i][3]-(d2f[i]+d3f[i]/3.0*v[i][3])/2.0*v[i][3];
}
}
}
}
void SetupQ()
{ v[0][3]=x[1]-x[0];
for (int i=1; i < lastElement; i++)
{ v[i][3]=x[i+1]-x[i];
v[i][0]=dy[i-1]/v[i-1][3];
v[i][1]=-dy[i]/v[i][3]-dy[i]/v[i-1][3];
v[i][2]=dy[i+1]/v[i][3];
}
v[lastElement][0]=0.0;
for (int i=1; i < lastElement; i++)
v[i][4]=Math.pow(v[i][0],2)+Math.pow(v[i][1],2)+Math.pow(v[i][2],2);
for (int i=2; i < lastElement; i++)
v[i-1][5]=v[i-1][1]*v[i][0]+v[i-1][2]*v[i][1];
v[lastElement-1][5]=0.0;
for (int i=3; i < lastElement; i++)
v[i-2][6]=v[i-2][2]*v[i][0];
v[lastElement-2][6]=0.0;
v[lastElement-1][6]=0.0;
// construct q-transp. * y in d3f.
double prev=(y[1]-y[0])/v[0][3];
for (int i=1; i < lastElement; i++)
{ double diff=(y[i+1]-y[i])/v[i][3];
d3f[i]=diff-prev;
prev=diff;
}
}
void Chol1D(double p)
{ double temp=6.0*(1.0-p);
for (int i=1; i < lastElement; i++)
{ v[i][0]=temp*v[i][4]+2.0*p*(v[i-1][3]+v[i][3]);
v[i][1]=temp*v[i][5]+p*v[i][3];
v[i][2]=temp*v[i][6];
}
if (lastElement >= 3)
{ // factorization
for (int i=1; i < lastElement-1; i++)
{ double ratio=v[i][1]/v[i][0];
v[i+1][0]-=ratio*v[i][1];
v[i+1][1]-=ratio*v[i][2];
v[i][1]=ratio;
ratio=v[i][2]/v[i][0];
v[i+2][0]-=ratio*v[i][2];
v[i][2]=ratio;
}
// forward substitution
d2f[0]=0.0;
d2f[1]=d3f[1];
v[0][2]=0.0;
for (int i=1; i < lastElement-1; i++)
d2f[i+1]=d3f[i+1]-v[i][1]*d2f[i]-v[i-1][2]*d2f[i-1];
// back substitution
d2f[lastElement]=0.0;
d2f[lastElement-1]=d2f[lastElement-1]/v[lastElement-1][0];
for (int i=lastElement-2; i > 0; i--)
d2f[i]=d2f[i]/v[i][0]-d2f[i+1]*v[i][1]-d2f[i+2]*v[i][2];
}
else
{ d2f[0]=0.0;
d2f[1]=d3f[1]/v[1][0];
d2f[2]=0.0;
}
// construct q*u
double prev=0.0;
for (int i=1; i <= lastElement; i++)
{ f[i]=(d2f[i]-d2f[i-1])/v[i-1][3];
f[i-1]=f[i]-prev;
prev=f[i];
}
f[lastElement]=-f[lastElement];
}
/**
* Constructs an interpolation function to given data (x(),y()).
* @param isSpline is type of interpolation function. If it is true
* the interpolation function is a cubic spline, otherwise
* a linear function.
*/
public void calcInterpolation(boolean isSpline)
{ if (lastElement > 0)
{ boolean splineInterpolation=(lastElement > 1) && isSpline;
double step, stepB, temp;
double delta=0.0; // What's this? ***
step=x[1]-x[0]; // boundary condition
v[0][3]=-step;
v[0][4]=(step*step)*delta;
for (int i=1; i < lastElement; i++)
{ step=x[i+1]-x[i];
stepB=x[i]-x[i-1];
v[i][3]=2.0*(stepB+step)-(stepB*stepB)/v[i-1][3];
v[i][4]=(y[i+1]-y[i])/step-(y[i]-y[i-1])/stepB-stepB*v[i-1][4]/v[i-1][3];
}
stepB=x[lastElement]-x[lastElement-1]; // boundary condition
v[lastElement][3]=-stepB-(stepB*stepB)/v[lastElement-1][3];
v[lastElement][4]=-(stepB*stepB)*delta-stepB*v[lastElement-1][4]/v[lastElement-1][3];
d2f[lastElement]=splineInterpolation ? v[lastElement][4]/v[lastElement][3] : 0.0;
for (int i=lastElement-1; i >= 0; i--)
d2f[i]=splineInterpolation ? (v[i][4]-(x[i+1]-x[i])*d2f[i+1])/v[i][3] : 0.0 ;
for (int i=0; i < lastElement; i++)
{ step=x[i+1]-x[i];
f[i]=y[i];
d1f[i]=(y[i+1]-y[i])/step-step*(d2f[i+1]+2.0*d2f[i]);
d3f[i]=6.0*(d2f[i+1]-d2f[i])/step;
d2f[i]*=6.0;
}
}
}
/** calculates value of the function.
* @param coord is a value of abscissa for searching function's value.
* @return value of the interpolation/approximation function for coord value of abscissa
* or value less then emptyCode for unappropriate abscissa's value.
*/
public double get(double coord)
{ double temp, result=emptyCode-0.5;
boolean find=false;
if (lastElement > 0)
{ int i;
for (i=0; (i < lastElement) &&
(!(find=(coord >= x[i]) && (coord <= x[i+1]))); i++){};
if (find)
{ temp=coord-x[i];
result=f[i]+temp*(d1f[i]+temp*(d2f[i]+temp*d3f[i]/3.0)/2.0);
}
}
return result;
}
}
Ads help maintain this website.