Context Navigation

source: trunk/sources/ALGLIB/spline3.cs @ 2515

Visit:

Last change on this file since 2515 was 2430, checked in by gkronber, 15 years ago
Moved ALGLIB code into a separate plugin. #783
File size: 40.2 KB

Rev	Line
[2154]	1	/*************************************************************************
	2	Copyright (c) 2007, Sergey Bochkanov (ALGLIB project).
	3
[2430]	4	>>> SOURCE LICENSE >>>
	5	This program is free software; you can redistribute it and/or modify
	6	it under the terms of the GNU General Public License as published by
	7	the Free Software Foundation (www.fsf.org); either version 2 of the
	8	License, or (at your option) any later version.
[2154]	9
[2430]	10	This program is distributed in the hope that it will be useful,
	11	but WITHOUT ANY WARRANTY; without even the implied warranty of
	12	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
	13	GNU General Public License for more details.
[2154]	14
[2430]	15	A copy of the GNU General Public License is available at
	16	http://www.fsf.org/licensing/licenses
[2154]	17
[2430]	18	>>> END OF LICENSE >>>
[2154]	19	*************************************************************************/
	20
	21	using System;
	22
[2430]	23	namespace alglib
[2154]	24	{
[2430]	25	public class spline3
	26	{
	27	/*************************************************************************
	28	This subroutine builds linear spline coefficients table.
[2154]	29
[2430]	30	Input parameters:
	31	X - spline nodes, array[0..N-1]
	32	Y - function values, array[0..N-1]
	33	N - points count, N>=2
	34
	35	Output parameters:
	36	C - coefficients table. Used by SplineInterpolation and other
	37	subroutines from this file.
[2154]	38
[2430]	39	-- ALGLIB PROJECT --
	40	Copyright 24.06.2007 by Bochkanov Sergey
	41	*************************************************************************/
	42	public static void buildlinearspline(double[] x,
	43	double[] y,
	44	int n,
	45	ref double[] c)
	46	{
	47	int i = 0;
	48	int tblsize = 0;
[2154]	49
[2430]	50	x = (double[])x.Clone();
	51	y = (double[])y.Clone();
[2154]	52
[2430]	53	System.Diagnostics.Debug.Assert(n>=2, "BuildLinearSpline: N<2!");
	54
	55	//
	56	// Sort points
	57	//
	58	heapsortpoints(ref x, ref y, n);
	59
	60	//
	61	// Fill C:
	62	// C[0] - length(C)
	63	// C[1] - type(C):
	64	// 3 - general cubic spline
	65	// C[2] - N
	66	// C[3]...C[3+N-1] - x[i], i = 0...N-1
	67	// C[3+N]...C[3+N+(N-1)*4-1] - coefficients table
	68	//
	69	tblsize = 3+n+(n-1)*4;
	70	c = new double[tblsize-1+1];
	71	c[0] = tblsize;
	72	c[1] = 3;
	73	c[2] = n;
	74	for(i=0; i<=n-1; i++)
	75	{
	76	c[3+i] = x[i];
	77	}
	78	for(i=0; i<=n-2; i++)
	79	{
	80	c[3+n+4*i+0] = y[i];
	81	c[3+n+4*i+1] = (y[i+1]-y[i])/(x[i+1]-x[i]);
	82	c[3+n+4*i+2] = 0;
	83	c[3+n+4*i+3] = 0;
	84	}
[2154]	85	}
	86
	87
[2430]	88	/*************************************************************************
	89	This subroutine builds cubic spline coefficients table.
[2154]	90
[2430]	91	Input parameters:
	92	X - spline nodes, array[0..N-1]
	93	Y - function values, array[0..N-1]
	94	N - points count, N>=2
	95	BoundLType - boundary condition type for the left boundary
	96	BoundL - left boundary condition (first or second derivative,
	97	depending on the BoundLType)
	98	BoundRType - boundary condition type for the right boundary
	99	BoundR - right boundary condition (first or second derivative,
	100	depending on the BoundRType)
[2154]	101
[2430]	102	Output parameters:
	103	C - coefficients table. Used by SplineInterpolation and
	104	other subroutines from this file.
	105
	106	The BoundLType/BoundRType parameters can have the following values:
	107	* 0, which corresponds to the parabolically terminated spline
	108	(BoundL/BoundR are ignored).
	109	* 1, which corresponds to the first derivative boundary condition
	110	* 2, which corresponds to the second derivative boundary condition
[2154]	111
[2430]	112	-- ALGLIB PROJECT --
	113	Copyright 23.06.2007 by Bochkanov Sergey
	114	*************************************************************************/
	115	public static void buildcubicspline(double[] x,
	116	double[] y,
	117	int n,
	118	int boundltype,
	119	double boundl,
	120	int boundrtype,
	121	double boundr,
	122	ref double[] c)
	123	{
	124	double[] a1 = new double[0];
	125	double[] a2 = new double[0];
	126	double[] a3 = new double[0];
	127	double[] b = new double[0];
	128	double[] d = new double[0];
	129	int i = 0;
	130	int tblsize = 0;
	131	double delta = 0;
	132	double delta2 = 0;
	133	double delta3 = 0;
[2154]	134
[2430]	135	x = (double[])x.Clone();
	136	y = (double[])y.Clone();
[2154]	137
[2430]	138	System.Diagnostics.Debug.Assert(n>=2, "BuildCubicSpline: N<2!");
	139	System.Diagnostics.Debug.Assert(boundltype==0 \| boundltype==1 \| boundltype==2, "BuildCubicSpline: incorrect BoundLType!");
	140	System.Diagnostics.Debug.Assert(boundrtype==0 \| boundrtype==1 \| boundrtype==2, "BuildCubicSpline: incorrect BoundRType!");
	141	a1 = new double[n-1+1];
	142	a2 = new double[n-1+1];
	143	a3 = new double[n-1+1];
	144	b = new double[n-1+1];
[2154]	145
	146	//
[2430]	147	// Special case:
	148	// * N=2
	149	// * parabolic terminated boundary condition on both ends
[2154]	150	//
[2430]	151	if( n==2 & boundltype==0 & boundrtype==0 )
	152	{
	153
	154	//
	155	// Change task type
	156	//
	157	boundltype = 2;
	158	boundl = 0;
	159	boundrtype = 2;
	160	boundr = 0;
	161	}
	162
	163	//
	164	//
	165	// Sort points
	166	//
	167	heapsortpoints(ref x, ref y, n);
	168
	169	//
	170	// Left boundary conditions
	171	//
	172	if( boundltype==0 )
	173	{
	174	a1[0] = 0;
	175	a2[0] = 1;
	176	a3[0] = 1;
	177	b[0] = 2*(y[1]-y[0])/(x[1]-x[0]);
	178	}
	179	if( boundltype==1 )
	180	{
	181	a1[0] = 0;
	182	a2[0] = 1;
	183	a3[0] = 0;
	184	b[0] = boundl;
	185	}
	186	if( boundltype==2 )
	187	{
	188	a1[0] = 0;
	189	a2[0] = 2;
	190	a3[0] = 1;
	191	b[0] = 3(y[1]-y[0])/(x[1]-x[0])-0.5boundl*(x[1]-x[0]);
	192	}
	193
	194	//
	195	// Central conditions
	196	//
	197	for(i=1; i<=n-2; i++)
	198	{
	199	a1[i] = x[i+1]-x[i];
	200	a2[i] = 2*(x[i+1]-x[i-1]);
	201	a3[i] = x[i]-x[i-1];
	202	b[i] = 3(y[i]-y[i-1])/(x[i]-x[i-1])(x[i+1]-x[i])+3(y[i+1]-y[i])/(x[i+1]-x[i])(x[i]-x[i-1]);
	203	}
	204
	205	//
	206	// Right boundary conditions
	207	//
	208	if( boundrtype==0 )
	209	{
	210	a1[n-1] = 1;
	211	a2[n-1] = 1;
	212	a3[n-1] = 0;
	213	b[n-1] = 2*(y[n-1]-y[n-2])/(x[n-1]-x[n-2]);
	214	}
	215	if( boundrtype==1 )
	216	{
	217	a1[n-1] = 0;
	218	a2[n-1] = 1;
	219	a3[n-1] = 0;
	220	b[n-1] = boundr;
	221	}
	222	if( boundrtype==2 )
	223	{
	224	a1[n-1] = 1;
	225	a2[n-1] = 2;
	226	a3[n-1] = 0;
	227	b[n-1] = 3(y[n-1]-y[n-2])/(x[n-1]-x[n-2])+0.5boundr*(x[n-1]-x[n-2]);
	228	}
	229
	230	//
	231	// Solve
	232	//
	233	solvetridiagonal(a1, a2, a3, b, n, ref d);
	234
	235	//
	236	// Now problem is reduced to the cubic Hermite spline
	237	//
	238	buildhermitespline(x, y, d, n, ref c);
[2154]	239	}
	240
	241
[2430]	242	/*************************************************************************
	243	This subroutine builds cubic Hermite spline coefficients table.
[2154]	244
[2430]	245	Input parameters:
	246	X - spline nodes, array[0..N-1]
	247	Y - function values, array[0..N-1]
	248	D - derivatives, array[0..N-1]
	249	N - points count, N>=2
[2154]	250
[2430]	251	Output parameters:
	252	C - coefficients table. Used by SplineInterpolation and
	253	other subroutines from this file.
[2154]	254
[2430]	255	-- ALGLIB PROJECT --
	256	Copyright 23.06.2007 by Bochkanov Sergey
	257	*************************************************************************/
	258	public static void buildhermitespline(double[] x,
	259	double[] y,
	260	double[] d,
	261	int n,
	262	ref double[] c)
	263	{
	264	int i = 0;
	265	int tblsize = 0;
	266	double delta = 0;
	267	double delta2 = 0;
	268	double delta3 = 0;
[2154]	269
[2430]	270	x = (double[])x.Clone();
	271	y = (double[])y.Clone();
	272	d = (double[])d.Clone();
[2154]	273
[2430]	274	System.Diagnostics.Debug.Assert(n>=2, "BuildHermiteSpline: N<2!");
	275
	276	//
	277	// Sort points
	278	//
	279	heapsortdpoints(ref x, ref y, ref d, n);
	280
	281	//
	282	// Fill C:
	283	// C[0] - length(C)
	284	// C[1] - type(C):
	285	// 3 - general cubic spline
	286	// C[2] - N
	287	// C[3]...C[3+N-1] - x[i], i = 0...N-1
	288	// C[3+N]...C[3+N+(N-1)*4-1] - coefficients table
	289	//
	290	tblsize = 3+n+(n-1)*4;
	291	c = new double[tblsize-1+1];
	292	c[0] = tblsize;
	293	c[1] = 3;
	294	c[2] = n;
	295	for(i=0; i<=n-1; i++)
	296	{
	297	c[3+i] = x[i];
	298	}
	299	for(i=0; i<=n-2; i++)
	300	{
	301	delta = x[i+1]-x[i];
	302	delta2 = AP.Math.Sqr(delta);
	303	delta3 = delta*delta2;
	304	c[3+n+4*i+0] = y[i];
	305	c[3+n+4*i+1] = d[i];
	306	c[3+n+4i+2] = (3(y[i+1]-y[i])-2d[i]delta-d[i+1]*delta)/delta2;
	307	c[3+n+4i+3] = (2(y[i]-y[i+1])+d[i]delta+d[i+1]delta)/delta3;
	308	}
[2154]	309	}
	310
	311
[2430]	312	/*************************************************************************
	313	This subroutine builds Akima spline coefficients table.
[2154]	314
[2430]	315	Input parameters:
	316	X - spline nodes, array[0..N-1]
	317	Y - function values, array[0..N-1]
	318	N - points count, N>=5
[2154]	319
[2430]	320	Output parameters:
	321	C - coefficients table. Used by SplineInterpolation and
	322	other subroutines from this file.
[2154]	323
[2430]	324	-- ALGLIB PROJECT --
	325	Copyright 24.06.2007 by Bochkanov Sergey
	326	*************************************************************************/
	327	public static void buildakimaspline(double[] x,
	328	double[] y,
	329	int n,
	330	ref double[] c)
	331	{
	332	int i = 0;
	333	double[] d = new double[0];
	334	double[] w = new double[0];
	335	double[] diff = new double[0];
[2154]	336
[2430]	337	x = (double[])x.Clone();
	338	y = (double[])y.Clone();
[2154]	339
[2430]	340	System.Diagnostics.Debug.Assert(n>=5, "BuildAkimaSpline: N<5!");
	341
	342	//
	343	// Sort points
	344	//
	345	heapsortpoints(ref x, ref y, n);
	346
	347	//
	348	// Prepare W (weights), Diff (divided differences)
	349	//
	350	w = new double[n-2+1];
	351	diff = new double[n-2+1];
	352	for(i=0; i<=n-2; i++)
[2154]	353	{
[2430]	354	diff[i] = (y[i+1]-y[i])/(x[i+1]-x[i]);
[2154]	355	}
[2430]	356	for(i=1; i<=n-2; i++)
[2154]	357	{
[2430]	358	w[i] = Math.Abs(diff[i]-diff[i-1]);
[2154]	359	}
[2430]	360
	361	//
	362	// Prepare Hermite interpolation scheme
	363	//
	364	d = new double[n-1+1];
	365	for(i=2; i<=n-3; i++)
	366	{
	367	if( Math.Abs(w[i-1])+Math.Abs(w[i+1])!=0 )
	368	{
	369	d[i] = (w[i+1]diff[i-1]+w[i-1]diff[i])/(w[i+1]+w[i-1]);
	370	}
	371	else
	372	{
	373	d[i] = ((x[i+1]-x[i])diff[i-1]+(x[i]-x[i-1])diff[i])/(x[i+1]-x[i-1]);
	374	}
	375	}
	376	d[0] = diffthreepoint(x[0], x[0], y[0], x[1], y[1], x[2], y[2]);
	377	d[1] = diffthreepoint(x[1], x[0], y[0], x[1], y[1], x[2], y[2]);
	378	d[n-2] = diffthreepoint(x[n-2], x[n-3], y[n-3], x[n-2], y[n-2], x[n-1], y[n-1]);
	379	d[n-1] = diffthreepoint(x[n-1], x[n-3], y[n-3], x[n-2], y[n-2], x[n-1], y[n-1]);
	380
	381	//
	382	// Build Akima spline using Hermite interpolation scheme
	383	//
	384	buildhermitespline(x, y, d, n, ref c);
[2154]	385	}
	386
	387
[2430]	388	/*************************************************************************
	389	This subroutine calculates the value of the spline at the given point X.
[2154]	390
[2430]	391	Input parameters:
	392	C - coefficients table. Built by BuildLinearSpline,
	393	BuildHermiteSpline, BuildCubicSpline, BuildAkimaSpline.
	394	X - point
[2154]	395
[2430]	396	Result:
	397	S(x)
[2154]	398
[2430]	399	-- ALGLIB PROJECT --
	400	Copyright 23.06.2007 by Bochkanov Sergey
	401	*************************************************************************/
	402	public static double splineinterpolation(ref double[] c,
	403	double x)
	404	{
	405	double result = 0;
	406	int n = 0;
	407	int l = 0;
	408	int r = 0;
	409	int m = 0;
[2154]	410
[2430]	411	System.Diagnostics.Debug.Assert((int)Math.Round(c[1])==3, "SplineInterpolation: incorrect C!");
	412	n = (int)Math.Round(c[2]);
	413
	414	//
	415	// Binary search in the [ x[0], ..., x[n-2] ] (x[n-1] is not included)
	416	//
	417	l = 3;
	418	r = 3+n-2+1;
	419	while( l!=r-1 )
[2154]	420	{
[2430]	421	m = (l+r)/2;
	422	if( c[m]>=x )
	423	{
	424	r = m;
	425	}
	426	else
	427	{
	428	l = m;
	429	}
[2154]	430	}
[2430]	431
	432	//
	433	// Interpolation
	434	//
	435	x = x-c[l];
	436	m = 3+n+4*(l-3);
	437	result = c[m]+x(c[m+1]+x(c[m+2]+x*c[m+3]));
	438	return result;
[2154]	439	}
	440
	441
[2430]	442	/*************************************************************************
	443	This subroutine differentiates the spline.
[2154]	444
[2430]	445	Input parameters:
	446	C - coefficients table. Built by BuildLinearSpline,
	447	BuildHermiteSpline, BuildCubicSpline, BuildAkimaSpline.
	448	X - point
[2154]	449
[2430]	450	Result:
	451	S - S(x)
	452	DS - S'(x)
	453	D2S - S''(x)
[2154]	454
[2430]	455	-- ALGLIB PROJECT --
	456	Copyright 24.06.2007 by Bochkanov Sergey
	457	*************************************************************************/
	458	public static void splinedifferentiation(ref double[] c,
	459	double x,
	460	ref double s,
	461	ref double ds,
	462	ref double d2s)
	463	{
	464	int n = 0;
	465	int l = 0;
	466	int r = 0;
	467	int m = 0;
[2154]	468
[2430]	469	System.Diagnostics.Debug.Assert((int)Math.Round(c[1])==3, "SplineInterpolation: incorrect C!");
	470	n = (int)Math.Round(c[2]);
	471
	472	//
	473	// Binary search
	474	//
	475	l = 3;
	476	r = 3+n-2+1;
	477	while( l!=r-1 )
[2154]	478	{
[2430]	479	m = (l+r)/2;
	480	if( c[m]>=x )
	481	{
	482	r = m;
	483	}
	484	else
	485	{
	486	l = m;
	487	}
[2154]	488	}
[2430]	489
	490	//
	491	// Differentiation
	492	//
	493	x = x-c[l];
	494	m = 3+n+4*(l-3);
	495	s = c[m]+x(c[m+1]+x(c[m+2]+x*c[m+3]));
	496	ds = c[m+1]+2xc[m+2]+3AP.Math.Sqr(x)c[m+3];
	497	d2s = 2c[m+2]+6x*c[m+3];
[2154]	498	}
	499
	500
[2430]	501	/*************************************************************************
	502	This subroutine makes the copy of the spline.
[2154]	503
[2430]	504	Input parameters:
	505	C - coefficients table. Built by BuildLinearSpline,
	506	BuildHermiteSpline, BuildCubicSpline, BuildAkimaSpline.
[2154]	507
[2430]	508	Result:
	509	CC - spline copy
[2154]	510
[2430]	511	-- ALGLIB PROJECT --
	512	Copyright 29.06.2007 by Bochkanov Sergey
	513	*************************************************************************/
	514	public static void splinecopy(ref double[] c,
	515	ref double[] cc)
	516	{
	517	int s = 0;
	518	int i_ = 0;
[2154]	519
[2430]	520	s = (int)Math.Round(c[0]);
	521	cc = new double[s-1+1];
	522	for(i_=0; i_<=s-1;i_++)
	523	{
	524	cc[i_] = c[i_];
	525	}
[2154]	526	}
	527
	528
[2430]	529	/*************************************************************************
	530	This subroutine unpacks the spline into the coefficients table.
[2154]	531
[2430]	532	Input parameters:
	533	C - coefficients table. Built by BuildLinearSpline,
	534	BuildHermiteSpline, BuildCubicSpline, BuildAkimaSpline.
	535	X - point
[2154]	536
[2430]	537	Result:
	538	Tbl - coefficients table, unpacked format, array[0..N-2, 0..5].
	539	For I = 0...N-2:
	540	Tbl[I,0] = X[i]
	541	Tbl[I,1] = X[i+1]
	542	Tbl[I,2] = C0
	543	Tbl[I,3] = C1
	544	Tbl[I,4] = C2
	545	Tbl[I,5] = C3
	546	On [x[i], x[i+1]] spline is equals to:
	547	S(x) = C0 + C1t + C2t^2 + C3*t^3
	548	t = x-x[i]
[2154]	549
[2430]	550	-- ALGLIB PROJECT --
	551	Copyright 29.06.2007 by Bochkanov Sergey
	552	*************************************************************************/
	553	public static void splineunpack(ref double[] c,
	554	ref int n,
	555	ref double[,] tbl)
	556	{
	557	int i = 0;
[2154]	558
[2430]	559	System.Diagnostics.Debug.Assert((int)Math.Round(c[1])==3, "SplineUnpack: incorrect C!");
	560	n = (int)Math.Round(c[2]);
	561	tbl = new double[n-2+1, 5+1];
	562
	563	//
	564	// Fill
	565	//
	566	for(i=0; i<=n-2; i++)
	567	{
	568	tbl[i,0] = c[3+i];
	569	tbl[i,1] = c[3+i+1];
	570	tbl[i,2] = c[3+n+4*i];
	571	tbl[i,3] = c[3+n+4*i+1];
	572	tbl[i,4] = c[3+n+4*i+2];
	573	tbl[i,5] = c[3+n+4*i+3];
	574	}
[2154]	575	}
	576
	577
[2430]	578	/*************************************************************************
	579	This subroutine performs linear transformation of the spline argument.
[2154]	580
[2430]	581	Input parameters:
	582	C - coefficients table. Built by BuildLinearSpline,
	583	BuildHermiteSpline, BuildCubicSpline, BuildAkimaSpline.
	584	A, B- transformation coefficients: x = A*t + B
	585	Result:
	586	C - transformed spline
[2154]	587
[2430]	588	-- ALGLIB PROJECT --
	589	Copyright 30.06.2007 by Bochkanov Sergey
	590	*************************************************************************/
	591	public static void splinelintransx(ref double[] c,
	592	double a,
	593	double b)
	594	{
	595	int i = 0;
	596	int n = 0;
	597	double v = 0;
	598	double dv = 0;
	599	double d2v = 0;
	600	double[] x = new double[0];
	601	double[] y = new double[0];
	602	double[] d = new double[0];
[2154]	603
[2430]	604	System.Diagnostics.Debug.Assert((int)Math.Round(c[1])==3, "SplineLinTransX: incorrect C!");
	605	n = (int)Math.Round(c[2]);
	606
	607	//
	608	// Special case: A=0
	609	//
	610	if( a==0 )
[2154]	611	{
[2430]	612	v = splineinterpolation(ref c, b);
	613	for(i=0; i<=n-2; i++)
	614	{
	615	c[3+n+4*i] = v;
	616	c[3+n+4*i+1] = 0;
	617	c[3+n+4*i+2] = 0;
	618	c[3+n+4*i+3] = 0;
	619	}
	620	return;
[2154]	621	}
[2430]	622
	623	//
	624	// General case: A<>0.
	625	// Unpack, X, Y, dY/dX.
	626	// Scale and pack again.
	627	//
	628	x = new double[n-1+1];
	629	y = new double[n-1+1];
	630	d = new double[n-1+1];
	631	for(i=0; i<=n-1; i++)
	632	{
	633	x[i] = c[3+i];
	634	splinedifferentiation(ref c, x[i], ref v, ref dv, ref d2v);
	635	x[i] = (x[i]-b)/a;
	636	y[i] = v;
	637	d[i] = a*dv;
	638	}
	639	buildhermitespline(x, y, d, n, ref c);
[2154]	640	}
	641
	642
[2430]	643	/*************************************************************************
	644	This subroutine performs linear transformation of the spline.
[2154]	645
[2430]	646	Input parameters:
	647	C - coefficients table. Built by BuildLinearSpline,
	648	BuildHermiteSpline, BuildCubicSpline, BuildAkimaSpline.
	649	A, B- transformation coefficients: S2(x) = A*S(x) + B
	650	Result:
	651	C - transformed spline
[2154]	652
[2430]	653	-- ALGLIB PROJECT --
	654	Copyright 30.06.2007 by Bochkanov Sergey
	655	*************************************************************************/
	656	public static void splinelintransy(ref double[] c,
	657	double a,
	658	double b)
	659	{
	660	int i = 0;
	661	int n = 0;
	662	double v = 0;
	663	double dv = 0;
	664	double d2v = 0;
	665	double[] x = new double[0];
	666	double[] y = new double[0];
	667	double[] d = new double[0];
[2154]	668
[2430]	669	System.Diagnostics.Debug.Assert((int)Math.Round(c[1])==3, "SplineLinTransX: incorrect C!");
	670	n = (int)Math.Round(c[2]);
	671
	672	//
	673	// Special case: A=0
	674	//
	675	for(i=0; i<=n-2; i++)
	676	{
	677	c[3+n+4i] = ac[3+n+4*i]+b;
	678	c[3+n+4i+1] = ac[3+n+4*i+1];
	679	c[3+n+4i+2] = ac[3+n+4*i+2];
	680	c[3+n+4i+3] = ac[3+n+4*i+3];
	681	}
[2154]	682	}
	683
	684
[2430]	685	/*************************************************************************
	686	This subroutine integrates the spline.
[2154]	687
[2430]	688	Input parameters:
	689	C - coefficients table. Built by BuildLinearSpline,
	690	BuildHermiteSpline, BuildCubicSpline, BuildAkimaSpline.
	691	X - right bound of the integration interval [a, x]
	692	Result:
	693	integral(S(t)dt,a,x)
[2154]	694
[2430]	695	-- ALGLIB PROJECT --
	696	Copyright 23.06.2007 by Bochkanov Sergey
	697	*************************************************************************/
	698	public static double splineintegration(ref double[] c,
	699	double x)
	700	{
	701	double result = 0;
	702	int n = 0;
	703	int i = 0;
	704	int l = 0;
	705	int r = 0;
	706	int m = 0;
	707	double w = 0;
[2154]	708
[2430]	709	System.Diagnostics.Debug.Assert((int)Math.Round(c[1])==3, "SplineIntegration: incorrect C!");
	710	n = (int)Math.Round(c[2]);
	711
	712	//
	713	// Binary search in the [ x[0], ..., x[n-2] ] (x[n-1] is not included)
	714	//
	715	l = 3;
	716	r = 3+n-2+1;
	717	while( l!=r-1 )
[2154]	718	{
[2430]	719	m = (l+r)/2;
	720	if( c[m]>=x )
	721	{
	722	r = m;
	723	}
	724	else
	725	{
	726	l = m;
	727	}
[2154]	728	}
[2430]	729
	730	//
	731	// Integration
	732	//
	733	result = 0;
	734	for(i=3; i<=l-1; i++)
[2154]	735	{
[2430]	736	w = c[i+1]-c[i];
	737	m = 3+n+4*(i-3);
	738	result = result+c[m]*w;
	739	result = result+c[m+1]*AP.Math.Sqr(w)/2;
	740	result = result+c[m+2]AP.Math.Sqr(w)w/3;
	741	result = result+c[m+3]*AP.Math.Sqr(AP.Math.Sqr(w))/4;
[2154]	742	}
[2430]	743	w = x-c[l];
	744	m = 3+n+4*(l-3);
[2154]	745	result = result+c[m]*w;
	746	result = result+c[m+1]*AP.Math.Sqr(w)/2;
	747	result = result+c[m+2]AP.Math.Sqr(w)w/3;
	748	result = result+c[m+3]*AP.Math.Sqr(AP.Math.Sqr(w))/4;
[2430]	749	return result;
[2154]	750	}
	751
	752
[2430]	753	/*************************************************************************
	754	Obsolete subroutine, left for backward compatibility.
	755	*************************************************************************/
	756	public static void spline3buildtable(int n,
	757	int diffn,
	758	double[] x,
	759	double[] y,
	760	double boundl,
	761	double boundr,
	762	ref double[,] ctbl)
	763	{
	764	bool c = new bool();
	765	int e = 0;
	766	int g = 0;
	767	double tmp = 0;
	768	int nxm1 = 0;
	769	int i = 0;
	770	int j = 0;
	771	double dx = 0;
	772	double dxj = 0;
	773	double dyj = 0;
	774	double dxjp1 = 0;
	775	double dyjp1 = 0;
	776	double dxp = 0;
	777	double dyp = 0;
	778	double yppa = 0;
	779	double yppb = 0;
	780	double pj = 0;
	781	double b1 = 0;
	782	double b2 = 0;
	783	double b3 = 0;
	784	double b4 = 0;
[2154]	785
[2430]	786	x = (double[])x.Clone();
	787	y = (double[])y.Clone();
[2154]	788
[2430]	789	n = n-1;
	790	g = (n+1)/2;
[2154]	791	do
	792	{
[2430]	793	i = g;
[2154]	794	do
	795	{
[2430]	796	j = i-g;
	797	c = true;
	798	do
[2154]	799	{
[2430]	800	if( x[j]<=x[j+g] )
	801	{
	802	c = false;
	803	}
	804	else
	805	{
	806	tmp = x[j];
	807	x[j] = x[j+g];
	808	x[j+g] = tmp;
	809	tmp = y[j];
	810	y[j] = y[j+g];
	811	y[j+g] = tmp;
	812	}
	813	j = j-1;
[2154]	814	}
[2430]	815	while( j>=0 & c );
	816	i = i+1;
	817	}
	818	while( i<=n );
	819	g = g/2;
	820	}
	821	while( g>0 );
	822	ctbl = new double[4+1, n+1];
	823	n = n+1;
	824	if( diffn==1 )
	825	{
	826	b1 = 1;
	827	b2 = 6/(x[1]-x[0])*((y[1]-y[0])/(x[1]-x[0])-boundl);
	828	b3 = 1;
	829	b4 = 6/(x[n-1]-x[n-2])*(boundr-(y[n-1]-y[n-2])/(x[n-1]-x[n-2]));
	830	}
	831	else
	832	{
	833	b1 = 0;
	834	b2 = 2*boundl;
	835	b3 = 0;
	836	b4 = 2*boundr;
	837	}
	838	nxm1 = n-1;
	839	if( n>=2 )
	840	{
	841	if( n>2 )
	842	{
	843	dxj = x[1]-x[0];
	844	dyj = y[1]-y[0];
	845	j = 2;
	846	while( j<=nxm1 )
[2154]	847	{
[2430]	848	dxjp1 = x[j]-x[j-1];
	849	dyjp1 = y[j]-y[j-1];
	850	dxp = dxj+dxjp1;
	851	ctbl[1,j-1] = dxjp1/dxp;
	852	ctbl[2,j-1] = 1-ctbl[1,j-1];
	853	ctbl[3,j-1] = 6*(dyjp1/dxjp1-dyj/dxj)/dxp;
	854	dxj = dxjp1;
	855	dyj = dyjp1;
	856	j = j+1;
[2154]	857	}
	858	}
[2430]	859	ctbl[1,0] = -(b1/2);
	860	ctbl[2,0] = b2/2;
	861	if( n!=2 )
	862	{
	863	j = 2;
	864	while( j<=nxm1 )
	865	{
	866	pj = ctbl[2,j-1]*ctbl[1,j-2]+2;
	867	ctbl[1,j-1] = -(ctbl[1,j-1]/pj);
	868	ctbl[2,j-1] = (ctbl[3,j-1]-ctbl[2,j-1]*ctbl[2,j-2])/pj;
	869	j = j+1;
	870	}
	871	}
	872	yppb = (b4-b3ctbl[2,nxm1-1])/(b3ctbl[1,nxm1-1]+2);
	873	i = 1;
	874	while( i<=nxm1 )
	875	{
	876	j = n-i;
	877	yppa = ctbl[1,j-1]*yppb+ctbl[2,j-1];
	878	dx = x[j]-x[j-1];
	879	ctbl[3,j-1] = (yppb-yppa)/dx/6;
	880	ctbl[2,j-1] = yppa/2;
	881	ctbl[1,j-1] = (y[j]-y[j-1])/dx-(ctbl[2,j-1]+ctbl[3,j-1]dx)dx;
	882	yppb = yppa;
	883	i = i+1;
	884	}
	885	for(i=1; i<=n; i++)
	886	{
	887	ctbl[0,i-1] = y[i-1];
	888	ctbl[4,i-1] = x[i-1];
	889	}
[2154]	890	}
	891	}
[2430]	892
	893
	894	/*************************************************************************
	895	Obsolete subroutine, left for backward compatibility.
	896	*************************************************************************/
	897	public static double spline3interpolate(int n,
	898	ref double[,] c,
	899	double x)
[2154]	900	{
[2430]	901	double result = 0;
	902	int i = 0;
	903	int l = 0;
	904	int half = 0;
	905	int first = 0;
	906	int middle = 0;
	907
	908	n = n-1;
	909	l = n;
	910	first = 0;
	911	while( l>0 )
[2154]	912	{
[2430]	913	half = l/2;
	914	middle = first+half;
	915	if( c[4,middle]<x )
[2154]	916	{
[2430]	917	first = middle+1;
	918	l = l-half-1;
[2154]	919	}
[2430]	920	else
[2154]	921	{
[2430]	922	l = half;
[2154]	923	}
	924	}
[2430]	925	i = first-1;
	926	if( i<0 )
[2154]	927	{
[2430]	928	i = 0;
[2154]	929	}
[2430]	930	result = c[0,i]+(x-c[4,i])(c[1,i]+(x-c[4,i])(c[2,i]+c[3,i]*(x-c[4,i])));
	931	return result;
[2154]	932	}
	933
	934
[2430]	935	/*************************************************************************
	936	Internal subroutine. Heap sort.
	937	*************************************************************************/
	938	private static void heapsortpoints(ref double[] x,
	939	ref double[] y,
	940	int n)
	941	{
	942	int i = 0;
	943	int j = 0;
	944	int k = 0;
	945	int t = 0;
	946	double tmp = 0;
	947	bool isascending = new bool();
	948	bool isdescending = new bool();
[2154]	949
[2430]	950
	951	//
	952	// Test for already sorted set
	953	//
	954	isascending = true;
	955	isdescending = true;
	956	for(i=1; i<=n-1; i++)
[2154]	957	{
[2430]	958	isascending = isascending & x[i]>x[i-1];
	959	isdescending = isdescending & x[i]<x[i-1];
[2154]	960	}
[2430]	961	if( isascending )
[2154]	962	{
[2430]	963	return;
[2154]	964	}
[2430]	965	if( isdescending )
[2154]	966	{
[2430]	967	for(i=0; i<=n-1; i++)
[2154]	968	{
[2430]	969	j = n-1-i;
	970	if( j<=i )
	971	{
	972	break;
	973	}
	974	tmp = x[i];
	975	x[i] = x[j];
	976	x[j] = tmp;
	977	tmp = y[i];
	978	y[i] = y[j];
	979	y[j] = tmp;
[2154]	980	}
[2430]	981	return;
[2154]	982	}
[2430]	983
	984	//
	985	// Special case: N=1
	986	//
	987	if( n==1 )
[2154]	988	{
[2430]	989	return;
[2154]	990	}
[2430]	991
	992	//
	993	// General case
	994	//
	995	i = 2;
	996	do
[2154]	997	{
[2430]	998	t = i;
	999	while( t!=1 )
[2154]	1000	{
[2430]	1001	k = t/2;
	1002	if( x[k-1]>=x[t-1] )
[2154]	1003	{
[2430]	1004	t = 1;
[2154]	1005	}
	1006	else
	1007	{
	1008	tmp = x[k-1];
	1009	x[k-1] = x[t-1];
	1010	x[t-1] = tmp;
	1011	tmp = y[k-1];
	1012	y[k-1] = y[t-1];
	1013	y[t-1] = tmp;
	1014	t = k;
	1015	}
	1016	}
[2430]	1017	i = i+1;
[2154]	1018	}
[2430]	1019	while( i<=n );
	1020	i = n-1;
	1021	do
[2154]	1022	{
[2430]	1023	tmp = x[i];
	1024	x[i] = x[0];
	1025	x[0] = tmp;
	1026	tmp = y[i];
	1027	y[i] = y[0];
	1028	y[0] = tmp;
	1029	t = 1;
	1030	while( t!=0 )
[2154]	1031	{
[2430]	1032	k = 2*t;
	1033	if( k>i )
	1034	{
	1035	t = 0;
	1036	}
	1037	else
	1038	{
	1039	if( k<i )
	1040	{
	1041	if( x[k]>x[k-1] )
	1042	{
	1043	k = k+1;
	1044	}
	1045	}
	1046	if( x[t-1]>=x[k-1] )
	1047	{
	1048	t = 0;
	1049	}
	1050	else
	1051	{
	1052	tmp = x[k-1];
	1053	x[k-1] = x[t-1];
	1054	x[t-1] = tmp;
	1055	tmp = y[k-1];
	1056	y[k-1] = y[t-1];
	1057	y[t-1] = tmp;
	1058	t = k;
	1059	}
	1060	}
[2154]	1061	}
[2430]	1062	i = i-1;
[2154]	1063	}
[2430]	1064	while( i>=1 );
[2154]	1065	}
[2430]	1066
	1067
	1068	/*************************************************************************
	1069	Internal subroutine. Heap sort.
	1070	*************************************************************************/
	1071	private static void heapsortdpoints(ref double[] x,
	1072	ref double[] y,
	1073	ref double[] d,
	1074	int n)
[2154]	1075	{
[2430]	1076	int i = 0;
	1077	int j = 0;
	1078	int k = 0;
	1079	int t = 0;
	1080	double tmp = 0;
	1081	bool isascending = new bool();
	1082	bool isdescending = new bool();
	1083
	1084
	1085	//
	1086	// Test for already sorted set
	1087	//
	1088	isascending = true;
	1089	isdescending = true;
	1090	for(i=1; i<=n-1; i++)
[2154]	1091	{
[2430]	1092	isascending = isascending & x[i]>x[i-1];
	1093	isdescending = isdescending & x[i]<x[i-1];
	1094	}
	1095	if( isascending )
	1096	{
	1097	return;
	1098	}
	1099	if( isdescending )
	1100	{
	1101	for(i=0; i<=n-1; i++)
[2154]	1102	{
[2430]	1103	j = n-1-i;
	1104	if( j<=i )
	1105	{
	1106	break;
	1107	}
	1108	tmp = x[i];
	1109	x[i] = x[j];
	1110	x[j] = tmp;
	1111	tmp = y[i];
	1112	y[i] = y[j];
	1113	y[j] = tmp;
	1114	tmp = d[i];
	1115	d[i] = d[j];
	1116	d[j] = tmp;
[2154]	1117	}
[2430]	1118	return;
[2154]	1119	}
[2430]	1120
	1121	//
	1122	// Special case: N=1
	1123	//
	1124	if( n==1 )
[2154]	1125	{
[2430]	1126	return;
	1127	}
	1128
	1129	//
	1130	// General case
	1131	//
	1132	i = 2;
	1133	do
	1134	{
	1135	t = i;
	1136	while( t!=1 )
[2154]	1137	{
[2430]	1138	k = t/2;
	1139	if( x[k-1]>=x[t-1] )
[2154]	1140	{
[2430]	1141	t = 1;
[2154]	1142	}
	1143	else
	1144	{
	1145	tmp = x[k-1];
	1146	x[k-1] = x[t-1];
	1147	x[t-1] = tmp;
	1148	tmp = y[k-1];
	1149	y[k-1] = y[t-1];
	1150	y[t-1] = tmp;
	1151	tmp = d[k-1];
	1152	d[k-1] = d[t-1];
	1153	d[t-1] = tmp;
	1154	t = k;
	1155	}
	1156	}
[2430]	1157	i = i+1;
[2154]	1158	}
[2430]	1159	while( i<=n );
	1160	i = n-1;
	1161	do
	1162	{
	1163	tmp = x[i];
	1164	x[i] = x[0];
	1165	x[0] = tmp;
	1166	tmp = y[i];
	1167	y[i] = y[0];
	1168	y[0] = tmp;
	1169	tmp = d[i];
	1170	d[i] = d[0];
	1171	d[0] = tmp;
	1172	t = 1;
	1173	while( t!=0 )
	1174	{
	1175	k = 2*t;
	1176	if( k>i )
	1177	{
	1178	t = 0;
	1179	}
	1180	else
	1181	{
	1182	if( k<i )
	1183	{
	1184	if( x[k]>x[k-1] )
	1185	{
	1186	k = k+1;
	1187	}
	1188	}
	1189	if( x[t-1]>=x[k-1] )
	1190	{
	1191	t = 0;
	1192	}
	1193	else
	1194	{
	1195	tmp = x[k-1];
	1196	x[k-1] = x[t-1];
	1197	x[t-1] = tmp;
	1198	tmp = y[k-1];
	1199	y[k-1] = y[t-1];
	1200	y[t-1] = tmp;
	1201	tmp = d[k-1];
	1202	d[k-1] = d[t-1];
	1203	d[t-1] = tmp;
	1204	t = k;
	1205	}
	1206	}
	1207	}
	1208	i = i-1;
	1209	}
	1210	while( i>=1 );
[2154]	1211	}
	1212
	1213
[2430]	1214	/*************************************************************************
	1215	Internal subroutine. Tridiagonal solver.
	1216	*************************************************************************/
	1217	private static void solvetridiagonal(double[] a,
	1218	double[] b,
	1219	double[] c,
	1220	double[] d,
	1221	int n,
	1222	ref double[] x)
	1223	{
	1224	int k = 0;
	1225	double t = 0;
[2154]	1226
[2430]	1227	a = (double[])a.Clone();
	1228	b = (double[])b.Clone();
	1229	c = (double[])c.Clone();
	1230	d = (double[])d.Clone();
[2154]	1231
[2430]	1232	x = new double[n-1+1];
	1233	a[0] = 0;
	1234	c[n-1] = 0;
	1235	for(k=1; k<=n-1; k++)
	1236	{
	1237	t = a[k]/b[k-1];
	1238	b[k] = b[k]-t*c[k-1];
	1239	d[k] = d[k]-t*d[k-1];
	1240	}
	1241	x[n-1] = d[n-1]/b[n-1];
	1242	for(k=n-2; k>=0; k--)
	1243	{
	1244	x[k] = (d[k]-c[k]*x[k+1])/b[k];
	1245	}
[2154]	1246	}
	1247
	1248
[2430]	1249	/*************************************************************************
	1250	Internal subroutine. Three-point differentiation
	1251	*************************************************************************/
	1252	private static double diffthreepoint(double t,
	1253	double x0,
	1254	double f0,
	1255	double x1,
	1256	double f1,
	1257	double x2,
	1258	double f2)
	1259	{
	1260	double result = 0;
	1261	double a = 0;
	1262	double b = 0;
[2154]	1263
[2430]	1264	t = t-x0;
	1265	x1 = x1-x0;
	1266	x2 = x2-x0;
	1267	a = (f2-f0-x2/x1(f1-f0))/(AP.Math.Sqr(x2)-x1x2);
	1268	b = (f1-f0-a*AP.Math.Sqr(x1))/x1;
	1269	result = 2at+b;
	1270	return result;
	1271	}
[2154]	1272	}
	1273	}

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Update cookies preferences