Context Navigation

source: branches/HeuristicLab.ALGLIB-2.5.0/ALGLIB-2.5.0/idwint.cs @ 5229

Visit:

Last change on this file since 5229 was 3839, checked in by mkommend, 14 years ago
implemented first version of LR (ticket #1012)
File size: 42.5 KB

Line
1	/*************************************************************************
2	Copyright (c) 2007-2010, Sergey Bochkanov (ALGLIB project).
3
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.
9
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.
14
15	A copy of the GNU General Public License is available at
16	http://www.fsf.org/licensing/licenses
17
18	>>> END OF LICENSE >>>
19	*************************************************************************/
20
21	using System;
22
23	namespace alglib
24	{
25	public class idwint
26	{
27	/*************************************************************************
28	IDW interpolant.
29	*************************************************************************/
30	public struct idwinterpolant
31	{
32	public int n;
33	public int nx;
34	public int d;
35	public double r;
36	public int nw;
37	public nearestneighbor.kdtree tree;
38	public int modeltype;
39	public double[,] q;
40	public double[] xbuf;
41	public int[] tbuf;
42	public double[] rbuf;
43	public double[,] xybuf;
44	public int debugsolverfailures;
45	public double debugworstrcond;
46	public double debugbestrcond;
47	};
48
49
50
51
52	public const double idwqfactor = 1.5;
53	public const int idwkmin = 5;
54
55
56	/*************************************************************************
57	IDW interpolation
58
59	INPUT PARAMETERS:
60	Z - IDW interpolant built with one of model building
61	subroutines.
62	X - array[0..NX-1], interpolation point
63
64	Result:
65	IDW interpolant Z(X)
66
67	-- ALGLIB --
68	Copyright 02.03.2010 by Bochkanov Sergey
69	*************************************************************************/
70	public static double idwcalc(ref idwinterpolant z,
71	ref double[] x)
72	{
73	double result = 0;
74	int nx = 0;
75	int i = 0;
76	int k = 0;
77	double r = 0;
78	double s = 0;
79	double w = 0;
80	double v1 = 0;
81	double v2 = 0;
82	double d0 = 0;
83	double di = 0;
84	double v = 0;
85
86	if( z.modeltype==0 )
87	{
88
89	//
90	// NQ/NW-based model
91	//
92	nx = z.nx;
93	nearestneighbor.kdtreequeryknn(ref z.tree, ref x, z.nw, true);
94	nearestneighbor.kdtreequeryresultsdistances(ref z.tree, ref z.rbuf, ref k);
95	nearestneighbor.kdtreequeryresultstags(ref z.tree, ref z.tbuf, ref k);
96	}
97	if( z.modeltype==1 )
98	{
99
100	//
101	// R-based model
102	//
103	nx = z.nx;
104	nearestneighbor.kdtreequeryrnn(ref z.tree, ref x, z.r, true);
105	nearestneighbor.kdtreequeryresultsdistances(ref z.tree, ref z.rbuf, ref k);
106	nearestneighbor.kdtreequeryresultstags(ref z.tree, ref z.tbuf, ref k);
107	if( k<idwkmin )
108	{
109
110	//
111	// we need at least IDWKMin points
112	//
113	nearestneighbor.kdtreequeryknn(ref z.tree, ref x, idwkmin, true);
114	nearestneighbor.kdtreequeryresultsdistances(ref z.tree, ref z.rbuf, ref k);
115	nearestneighbor.kdtreequeryresultstags(ref z.tree, ref z.tbuf, ref k);
116	}
117	}
118
119	//
120	// initialize weights for linear/quadratic members calculation.
121	//
122	// NOTE 1: weights are calculated using NORMALIZED modified
123	// Shepard's formula. Original formula gives w(i) = sqr((R-di)/(R*di)),
124	// where di is i-th distance, R is max(di). Modified formula have
125	// following form:
126	// w_mod(i) = 1, if di=d0
127	// w_mod(i) = w(i)/w(0), if di<>d0
128	//
129	// NOTE 2: self-match is USED for this query
130	//
131	// NOTE 3: last point almost always gain zero weight, but it MUST
132	// be used for fitting because sometimes it will gain NON-ZERO
133	// weight - for example, when all distances are equal.
134	//
135	r = z.rbuf[k-1];
136	d0 = z.rbuf[0];
137	result = 0;
138	s = 0;
139	for(i=0; i<=k-1; i++)
140	{
141	di = z.rbuf[i];
142	if( (double)(di)==(double)(d0) )
143	{
144
145	//
146	// distance is equal to shortest, set it 1.0
147	// without explicitly calculating (which would give
148	// us same result, but 'll expose us to the risk of
149	// division by zero).
150	//
151	w = 1;
152	}
153	else
154	{
155
156	//
157	// use normalized formula
158	//
159	v1 = (r-di)/(r-d0);
160	v2 = d0/di;
161	w = AP.Math.Sqr(v1*v2);
162	}
163	result = result+w*idwcalcq(ref z, ref x, z.tbuf[i]);
164	s = s+w;
165	}
166	result = result/s;
167	return result;
168	}
169
170
171	/*************************************************************************
172	IDW interpolant using modified Shepard method for uniform point
173	distributions.
174
175	INPUT PARAMETERS:
176	XY - X and Y values, array[0..N-1,0..NX].
177	First NX columns contain X-values, last column contain
178	Y-values.
179	N - number of nodes, N>0.
180	NX - space dimension, NX>=1.
181	D - nodal function type, either:
182	* 0 constant model. Just for demonstration only, worst
183	model ever.
184	* 1 linear model, least squares fitting. Simpe model for
185	datasets too small for quadratic models
186	* 2 quadratic model, least squares fitting. Best model
187	available (if your dataset is large enough).
188	* -1 "fast" linear model, use with caution!!! It is
189	significantly faster than linear/quadratic and better
190	than constant model. But it is less robust (especially
191	in the presence of noise).
192	NQ - number of points used to calculate nodal functions (ignored
193	for constant models). NQ should be LARGER than:
194	* max(1.5*(1+NX),2^NX+1) for linear model,
195	* max(3/4(NX+2)(NX+1),2^NX+1) for quadratic model.
196	Values less than this threshold will be silently increased.
197	NW - number of points used to calculate weights and to interpolate.
198	Required: >=2^NX+1, values less than this threshold will be
199	silently increased.
200	Recommended value: about 2*NQ
201
202	OUTPUT PARAMETERS:
203	Z - IDW interpolant.
204
205	NOTES:
206	* best results are obtained with quadratic models, worst - with constant
207	models
208	* when N is large, NQ and NW must be significantly smaller than N both
209	to obtain optimal performance and to obtain optimal accuracy. In 2 or
210	3-dimensional tasks NQ=15 and NW=25 are good values to start with.
211	* NQ and NW may be greater than N. In such cases they will be
212	automatically decreased.
213	* this subroutine is always succeeds (as long as correct parameters are
214	passed).
215	* see 'Multivariate Interpolation of Large Sets of Scattered Data' by
216	Robert J. Renka for more information on this algorithm.
217	* this subroutine assumes that point distribution is uniform at the small
218	scales. If it isn't - for example, points are concentrated along
219	"lines", but "lines" distribution is uniform at the larger scale - then
220	you should use IDWBuildModifiedShepardR()
221
222
223	-- ALGLIB PROJECT --
224	Copyright 02.03.2010 by Bochkanov Sergey
225	*************************************************************************/
226	public static void idwbuildmodifiedshepard(ref double[,] xy,
227	int n,
228	int nx,
229	int d,
230	int nq,
231	int nw,
232	ref idwinterpolant z)
233	{
234	int i = 0;
235	int j = 0;
236	int k = 0;
237	int j2 = 0;
238	int j3 = 0;
239	double v = 0;
240	double r = 0;
241	double s = 0;
242	double d0 = 0;
243	double di = 0;
244	double v1 = 0;
245	double v2 = 0;
246	int nc = 0;
247	int offs = 0;
248	double[] x = new double[0];
249	double[] qrbuf = new double[0];
250	double[,] qxybuf = new double[0,0];
251	double[] y = new double[0];
252	double[,] fmatrix = new double[0,0];
253	double[] w = new double[0];
254	double[] qsol = new double[0];
255	double[] temp = new double[0];
256	int[] tags = new int[0];
257	int info = 0;
258	double taskrcond = 0;
259	int i_ = 0;
260
261
262	//
263	// assertions
264	//
265	System.Diagnostics.Debug.Assert(n>0, "IDWBuildModifiedShepard: N<=0!");
266	System.Diagnostics.Debug.Assert(nx>=1, "IDWBuildModifiedShepard: NX<1!");
267	System.Diagnostics.Debug.Assert(d>=-1 & d<=2, "IDWBuildModifiedShepard: D<>-1 and D<>0 and D<>1 and D<>2!");
268
269	//
270	// Correct parameters if needed
271	//
272	if( d==1 )
273	{
274	nq = Math.Max(nq, (int)Math.Ceiling(idwqfactor*(1+nx))+1);
275	nq = Math.Max(nq, (int)Math.Round(Math.Pow(2, nx))+1);
276	}
277	if( d==2 )
278	{
279	nq = Math.Max(nq, (int)Math.Ceiling(idwqfactor(nx+2)(nx+1)/2)+1);
280	nq = Math.Max(nq, (int)Math.Round(Math.Pow(2, nx))+1);
281	}
282	nw = Math.Max(nw, (int)Math.Round(Math.Pow(2, nx))+1);
283	nq = Math.Min(nq, n);
284	nw = Math.Min(nw, n);
285
286	//
287	// primary initialization of Z
288	//
289	idwinit1(n, nx, d, nq, nw, ref z);
290	z.modeltype = 0;
291
292	//
293	// Create KD-tree
294	//
295	tags = new int[n];
296	for(i=0; i<=n-1; i++)
297	{
298	tags[i] = i;
299	}
300	nearestneighbor.kdtreebuildtagged(ref xy, ref tags, n, nx, 1, 2, ref z.tree);
301
302	//
303	// build nodal functions
304	//
305	temp = new double[nq+1];
306	x = new double[nx];
307	qrbuf = new double[nq];
308	qxybuf = new double[nq, nx+1];
309	if( d==-1 )
310	{
311	w = new double[nq];
312	}
313	if( d==1 )
314	{
315	y = new double[nq];
316	w = new double[nq];
317	qsol = new double[nx];
318
319	//
320	// NX for linear members,
321	// 1 for temporary storage
322	//
323	fmatrix = new double[nq, nx+1];
324	}
325	if( d==2 )
326	{
327	y = new double[nq];
328	w = new double[nq];
329	qsol = new double[nx+(int)Math.Round(nx(nx+1)0.5)];
330
331	//
332	// NX for linear members,
333	// Round(NX(NX+1)0.5) for quadratic model,
334	// 1 for temporary storage
335	//
336	fmatrix = new double[nq, nx+(int)Math.Round(nx(nx+1)0.5)+1];
337	}
338	for(i=0; i<=n-1; i++)
339	{
340
341	//
342	// Initialize center and function value.
343	// If D=0 it is all what we need
344	//
345	for(i_=0; i_<=nx;i_++)
346	{
347	z.q[i,i_] = xy[i,i_];
348	}
349	if( d==0 )
350	{
351	continue;
352	}
353
354	//
355	// calculate weights for linear/quadratic members calculation.
356	//
357	// NOTE 1: weights are calculated using NORMALIZED modified
358	// Shepard's formula. Original formula is w(i) = sqr((R-di)/(R*di)),
359	// where di is i-th distance, R is max(di). Modified formula have
360	// following form:
361	// w_mod(i) = 1, if di=d0
362	// w_mod(i) = w(i)/w(0), if di<>d0
363	//
364	// NOTE 2: self-match is NOT used for this query
365	//
366	// NOTE 3: last point almost always gain zero weight, but it MUST
367	// be used for fitting because sometimes it will gain NON-ZERO
368	// weight - for example, when all distances are equal.
369	//
370	for(i_=0; i_<=nx-1;i_++)
371	{
372	x[i_] = xy[i,i_];
373	}
374	nearestneighbor.kdtreequeryknn(ref z.tree, ref x, nq, false);
375	nearestneighbor.kdtreequeryresultsxy(ref z.tree, ref qxybuf, ref k);
376	nearestneighbor.kdtreequeryresultsdistances(ref z.tree, ref qrbuf, ref k);
377	r = qrbuf[k-1];
378	d0 = qrbuf[0];
379	for(j=0; j<=k-1; j++)
380	{
381	di = qrbuf[j];
382	if( (double)(di)==(double)(d0) )
383	{
384
385	//
386	// distance is equal to shortest, set it 1.0
387	// without explicitly calculating (which would give
388	// us same result, but 'll expose us to the risk of
389	// division by zero).
390	//
391	w[j] = 1;
392	}
393	else
394	{
395
396	//
397	// use normalized formula
398	//
399	v1 = (r-di)/(r-d0);
400	v2 = d0/di;
401	w[j] = AP.Math.Sqr(v1*v2);
402	}
403	}
404
405	//
406	// calculate linear/quadratic members
407	//
408	if( d==-1 )
409	{
410
411	//
412	// "Fast" linear nodal function calculated using
413	// inverse distance weighting
414	//
415	for(j=0; j<=nx-1; j++)
416	{
417	x[j] = 0;
418	}
419	s = 0;
420	for(j=0; j<=k-1; j++)
421	{
422
423	//
424	// calculate J-th inverse distance weighted gradient:
425	// grad_k = (y_j-y_k)*(x_j-x_k)/sqr(norm(x_j-x_k))
426	// grad = sum(wk*grad_k)/sum(w_k)
427	//
428	v = 0;
429	for(j2=0; j2<=nx-1; j2++)
430	{
431	v = v+AP.Math.Sqr(qxybuf[j,j2]-xy[i,j2]);
432	}
433
434	//
435	// Although x_j<>x_k, sqr(norm(x_j-x_k)) may be zero due to
436	// underflow. If it is, we assume than J-th gradient is zero
437	// (i.e. don't add anything)
438	//
439	if( (double)(v)!=(double)(0) )
440	{
441	for(j2=0; j2<=nx-1; j2++)
442	{
443	x[j2] = x[j2]+w[j](qxybuf[j,nx]-xy[i,nx])(qxybuf[j,j2]-xy[i,j2])/v;
444	}
445	}
446	s = s+w[j];
447	}
448	for(j=0; j<=nx-1; j++)
449	{
450	z.q[i,nx+1+j] = x[j]/s;
451	}
452	}
453	else
454	{
455
456	//
457	// Least squares models: build
458	//
459	if( d==1 )
460	{
461
462	//
463	// Linear nodal function calculated using
464	// least squares fitting to its neighbors
465	//
466	for(j=0; j<=k-1; j++)
467	{
468	for(j2=0; j2<=nx-1; j2++)
469	{
470	fmatrix[j,j2] = qxybuf[j,j2]-xy[i,j2];
471	}
472	y[j] = qxybuf[j,nx]-xy[i,nx];
473	}
474	nc = nx;
475	}
476	if( d==2 )
477	{
478
479	//
480	// Quadratic nodal function calculated using
481	// least squares fitting to its neighbors
482	//
483	for(j=0; j<=k-1; j++)
484	{
485	offs = 0;
486	for(j2=0; j2<=nx-1; j2++)
487	{
488	fmatrix[j,offs] = qxybuf[j,j2]-xy[i,j2];
489	offs = offs+1;
490	}
491	for(j2=0; j2<=nx-1; j2++)
492	{
493	for(j3=j2; j3<=nx-1; j3++)
494	{
495	fmatrix[j,offs] = (qxybuf[j,j2]-xy[i,j2])*(qxybuf[j,j3]-xy[i,j3]);
496	offs = offs+1;
497	}
498	}
499	y[j] = qxybuf[j,nx]-xy[i,nx];
500	}
501	nc = nx+(int)Math.Round(nx(nx+1)0.5);
502	}
503	idwinternalsolver(ref y, ref w, ref fmatrix, ref temp, k, nc, ref info, ref qsol, ref taskrcond);
504
505	//
506	// Least squares models: copy results
507	//
508	if( info>0 )
509	{
510
511	//
512	// LLS task is solved, copy results
513	//
514	z.debugworstrcond = Math.Min(z.debugworstrcond, taskrcond);
515	z.debugbestrcond = Math.Max(z.debugbestrcond, taskrcond);
516	for(j=0; j<=nc-1; j++)
517	{
518	z.q[i,nx+1+j] = qsol[j];
519	}
520	}
521	else
522	{
523
524	//
525	// Solver failure, very strange, but we will use
526	// zero values to handle it.
527	//
528	z.debugsolverfailures = z.debugsolverfailures+1;
529	for(j=0; j<=nc-1; j++)
530	{
531	z.q[i,nx+1+j] = 0;
532	}
533	}
534	}
535	}
536	}
537
538
539	/*************************************************************************
540	IDW interpolant using modified Shepard method for non-uniform datasets.
541
542	This type of model uses constant nodal functions and interpolates using
543	all nodes which are closer than user-specified radius R. It may be used
544	when points distribution is non-uniform at the small scale, but it is at
545	the distances as large as R.
546
547	INPUT PARAMETERS:
548	XY - X and Y values, array[0..N-1,0..NX].
549	First NX columns contain X-values, last column contain
550	Y-values.
551	N - number of nodes, N>0.
552	NX - space dimension, NX>=1.
553	R - radius, R>0
554
555	OUTPUT PARAMETERS:
556	Z - IDW interpolant.
557
558	NOTES:
559	* if there is less than IDWKMin points within R-ball, algorithm selects
560	IDWKMin closest ones, so that continuity properties of interpolant are
561	preserved even far from points.
562
563	-- ALGLIB PROJECT --
564	Copyright 11.04.2010 by Bochkanov Sergey
565	*************************************************************************/
566	public static void idwbuildmodifiedshepardr(ref double[,] xy,
567	int n,
568	int nx,
569	double r,
570	ref idwinterpolant z)
571	{
572	int i = 0;
573	int[] tags = new int[0];
574	int i_ = 0;
575
576
577	//
578	// assertions
579	//
580	System.Diagnostics.Debug.Assert(n>0, "IDWBuildModifiedShepardR: N<=0!");
581	System.Diagnostics.Debug.Assert(nx>=1, "IDWBuildModifiedShepardR: NX<1!");
582	System.Diagnostics.Debug.Assert((double)(r)>(double)(0), "IDWBuildModifiedShepardR: R<=0!");
583
584	//
585	// primary initialization of Z
586	//
587	idwinit1(n, nx, 0, 0, n, ref z);
588	z.modeltype = 1;
589	z.r = r;
590
591	//
592	// Create KD-tree
593	//
594	tags = new int[n];
595	for(i=0; i<=n-1; i++)
596	{
597	tags[i] = i;
598	}
599	nearestneighbor.kdtreebuildtagged(ref xy, ref tags, n, nx, 1, 2, ref z.tree);
600
601	//
602	// build nodal functions
603	//
604	for(i=0; i<=n-1; i++)
605	{
606	for(i_=0; i_<=nx;i_++)
607	{
608	z.q[i,i_] = xy[i,i_];
609	}
610	}
611	}
612
613
614	/*************************************************************************
615	IDW model for noisy data.
616
617	This subroutine may be used to handle noisy data, i.e. data with noise in
618	OUTPUT values. It differs from IDWBuildModifiedShepard() in the following
619	aspects:
620	* nodal functions are not constrained to pass through nodes: Qi(xi)<>yi,
621	i.e. we have fitting instead of interpolation.
622	* weights which are used during least squares fitting stage are all equal
623	to 1.0 (independently of distance)
624	* "fast"-linear or constant nodal functions are not supported (either not
625	robust enough or too rigid)
626
627	This problem require far more complex tuning than interpolation problems.
628	Below you can find some recommendations regarding this problem:
629	* focus on tuning NQ; it controls noise reduction. As for NW, you can just
630	make it equal to 2*NQ.
631	* you can use cross-validation to determine optimal NQ.
632	* optimal NQ is a result of complex tradeoff between noise level (more
633	noise = larger NQ required) and underlying function complexity (given
634	fixed N, larger NQ means smoothing of compex features in the data). For
635	example, NQ=N will reduce noise to the minimum level possible, but you
636	will end up with just constant/linear/quadratic (depending on D) least
637	squares model for the whole dataset.
638
639	INPUT PARAMETERS:
640	XY - X and Y values, array[0..N-1,0..NX].
641	First NX columns contain X-values, last column contain
642	Y-values.
643	N - number of nodes, N>0.
644	NX - space dimension, NX>=1.
645	D - nodal function degree, either:
646	* 1 linear model, least squares fitting. Simpe model for
647	datasets too small for quadratic models (or for very
648	noisy problems).
649	* 2 quadratic model, least squares fitting. Best model
650	available (if your dataset is large enough).
651	NQ - number of points used to calculate nodal functions. NQ should
652	be significantly larger than 1.5 times the number of
653	coefficients in a nodal function to overcome effects of noise:
654	* larger than 1.5*(1+NX) for linear model,
655	* larger than 3/4(NX+2)(NX+1) for quadratic model.
656	Values less than this threshold will be silently increased.
657	NW - number of points used to calculate weights and to interpolate.
658	Required: >=2^NX+1, values less than this threshold will be
659	silently increased.
660	Recommended value: about 2*NQ or larger
661
662	OUTPUT PARAMETERS:
663	Z - IDW interpolant.
664
665	NOTES:
666	* best results are obtained with quadratic models, linear models are not
667	recommended to use unless you are pretty sure that it is what you want
668	* this subroutine is always succeeds (as long as correct parameters are
669	passed).
670	* see 'Multivariate Interpolation of Large Sets of Scattered Data' by
671	Robert J. Renka for more information on this algorithm.
672
673
674	-- ALGLIB PROJECT --
675	Copyright 02.03.2010 by Bochkanov Sergey
676	*************************************************************************/
677	public static void idwbuildnoisy(ref double[,] xy,
678	int n,
679	int nx,
680	int d,
681	int nq,
682	int nw,
683	ref idwinterpolant z)
684	{
685	int i = 0;
686	int j = 0;
687	int k = 0;
688	int j2 = 0;
689	int j3 = 0;
690	double v = 0;
691	int nc = 0;
692	int offs = 0;
693	double taskrcond = 0;
694	double[] x = new double[0];
695	double[] qrbuf = new double[0];
696	double[,] qxybuf = new double[0,0];
697	double[] y = new double[0];
698	double[] w = new double[0];
699	double[,] fmatrix = new double[0,0];
700	double[] qsol = new double[0];
701	int[] tags = new int[0];
702	double[] temp = new double[0];
703	int info = 0;
704	int i_ = 0;
705
706
707	//
708	// assertions
709	//
710	System.Diagnostics.Debug.Assert(n>0, "IDWBuildNoisy: N<=0!");
711	System.Diagnostics.Debug.Assert(nx>=1, "IDWBuildNoisy: NX<1!");
712	System.Diagnostics.Debug.Assert(d>=1 & d<=2, "IDWBuildNoisy: D<>1 and D<>2!");
713
714	//
715	// Correct parameters if needed
716	//
717	if( d==1 )
718	{
719	nq = Math.Max(nq, (int)Math.Ceiling(idwqfactor*(1+nx))+1);
720	}
721	if( d==2 )
722	{
723	nq = Math.Max(nq, (int)Math.Ceiling(idwqfactor(nx+2)(nx+1)/2)+1);
724	}
725	nw = Math.Max(nw, (int)Math.Round(Math.Pow(2, nx))+1);
726	nq = Math.Min(nq, n);
727	nw = Math.Min(nw, n);
728
729	//
730	// primary initialization of Z
731	//
732	idwinit1(n, nx, d, nq, nw, ref z);
733	z.modeltype = 0;
734
735	//
736	// Create KD-tree
737	//
738	tags = new int[n];
739	for(i=0; i<=n-1; i++)
740	{
741	tags[i] = i;
742	}
743	nearestneighbor.kdtreebuildtagged(ref xy, ref tags, n, nx, 1, 2, ref z.tree);
744
745	//
746	// build nodal functions
747	// (special algorithm for noisy data is used)
748	//
749	temp = new double[nq+1];
750	x = new double[nx];
751	qrbuf = new double[nq];
752	qxybuf = new double[nq, nx+1];
753	if( d==1 )
754	{
755	y = new double[nq];
756	w = new double[nq];
757	qsol = new double[1+nx];
758
759	//
760	// 1 for constant member,
761	// NX for linear members,
762	// 1 for temporary storage
763	//
764	fmatrix = new double[nq, 1+nx+1];
765	}
766	if( d==2 )
767	{
768	y = new double[nq];
769	w = new double[nq];
770	qsol = new double[1+nx+(int)Math.Round(nx(nx+1)0.5)];
771
772	//
773	// 1 for constant member,
774	// NX for linear members,
775	// Round(NX(NX+1)0.5) for quadratic model,
776	// 1 for temporary storage
777	//
778	fmatrix = new double[nq, 1+nx+(int)Math.Round(nx(nx+1)0.5)+1];
779	}
780	for(i=0; i<=n-1; i++)
781	{
782
783	//
784	// Initialize center.
785	//
786	for(i_=0; i_<=nx-1;i_++)
787	{
788	z.q[i,i_] = xy[i,i_];
789	}
790
791	//
792	// Calculate linear/quadratic members
793	// using least squares fit
794	// NOTE 1: all weight are equal to 1.0
795	// NOTE 2: self-match is USED for this query
796	//
797	for(i_=0; i_<=nx-1;i_++)
798	{
799	x[i_] = xy[i,i_];
800	}
801	nearestneighbor.kdtreequeryknn(ref z.tree, ref x, nq, true);
802	nearestneighbor.kdtreequeryresultsxy(ref z.tree, ref qxybuf, ref k);
803	nearestneighbor.kdtreequeryresultsdistances(ref z.tree, ref qrbuf, ref k);
804	if( d==1 )
805	{
806
807	//
808	// Linear nodal function calculated using
809	// least squares fitting to its neighbors
810	//
811	for(j=0; j<=k-1; j++)
812	{
813	fmatrix[j,0] = 1.0;
814	for(j2=0; j2<=nx-1; j2++)
815	{
816	fmatrix[j,1+j2] = qxybuf[j,j2]-xy[i,j2];
817	}
818	y[j] = qxybuf[j,nx];
819	w[j] = 1;
820	}
821	nc = 1+nx;
822	}
823	if( d==2 )
824	{
825
826	//
827	// Quadratic nodal function calculated using
828	// least squares fitting to its neighbors
829	//
830	for(j=0; j<=k-1; j++)
831	{
832	fmatrix[j,0] = 1;
833	offs = 1;
834	for(j2=0; j2<=nx-1; j2++)
835	{
836	fmatrix[j,offs] = qxybuf[j,j2]-xy[i,j2];
837	offs = offs+1;
838	}
839	for(j2=0; j2<=nx-1; j2++)
840	{
841	for(j3=j2; j3<=nx-1; j3++)
842	{
843	fmatrix[j,offs] = (qxybuf[j,j2]-xy[i,j2])*(qxybuf[j,j3]-xy[i,j3]);
844	offs = offs+1;
845	}
846	}
847	y[j] = qxybuf[j,nx];
848	w[j] = 1;
849	}
850	nc = 1+nx+(int)Math.Round(nx(nx+1)0.5);
851	}
852	idwinternalsolver(ref y, ref w, ref fmatrix, ref temp, k, nc, ref info, ref qsol, ref taskrcond);
853
854	//
855	// Least squares models: copy results
856	//
857	if( info>0 )
858	{
859
860	//
861	// LLS task is solved, copy results
862	//
863	z.debugworstrcond = Math.Min(z.debugworstrcond, taskrcond);
864	z.debugbestrcond = Math.Max(z.debugbestrcond, taskrcond);
865	for(j=0; j<=nc-1; j++)
866	{
867	z.q[i,nx+j] = qsol[j];
868	}
869	}
870	else
871	{
872
873	//
874	// Solver failure, very strange, but we will use
875	// zero values to handle it.
876	//
877	z.debugsolverfailures = z.debugsolverfailures+1;
878	v = 0;
879	for(j=0; j<=k-1; j++)
880	{
881	v = v+qxybuf[j,nx];
882	}
883	z.q[i,nx] = v/k;
884	for(j=0; j<=nc-2; j++)
885	{
886	z.q[i,nx+1+j] = 0;
887	}
888	}
889	}
890	}
891
892
893	/*************************************************************************
894	Internal subroutine: K-th nodal function calculation
895
896	-- ALGLIB --
897	Copyright 02.03.2010 by Bochkanov Sergey
898	*************************************************************************/
899	private static double idwcalcq(ref idwinterpolant z,
900	ref double[] x,
901	int k)
902	{
903	double result = 0;
904	int nx = 0;
905	int i = 0;
906	int j = 0;
907	int offs = 0;
908
909	nx = z.nx;
910
911	//
912	// constant member
913	//
914	result = z.q[k,nx];
915
916	//
917	// linear members
918	//
919	if( z.d>=1 )
920	{
921	for(i=0; i<=nx-1; i++)
922	{
923	result = result+z.q[k,nx+1+i]*(x[i]-z.q[k,i]);
924	}
925	}
926
927	//
928	// quadratic members
929	//
930	if( z.d>=2 )
931	{
932	offs = nx+1+nx;
933	for(i=0; i<=nx-1; i++)
934	{
935	for(j=i; j<=nx-1; j++)
936	{
937	result = result+z.q[k,offs](x[i]-z.q[k,i])(x[j]-z.q[k,j]);
938	offs = offs+1;
939	}
940	}
941	}
942	return result;
943	}
944
945
946	/*************************************************************************
947	Initialization of internal structures.
948
949	It assumes correctness of all parameters.
950
951	-- ALGLIB --
952	Copyright 02.03.2010 by Bochkanov Sergey
953	*************************************************************************/
954	private static void idwinit1(int n,
955	int nx,
956	int d,
957	int nq,
958	int nw,
959	ref idwinterpolant z)
960	{
961	z.debugsolverfailures = 0;
962	z.debugworstrcond = 1.0;
963	z.debugbestrcond = 0;
964	z.n = n;
965	z.nx = nx;
966	z.d = 0;
967	if( d==1 )
968	{
969	z.d = 1;
970	}
971	if( d==2 )
972	{
973	z.d = 2;
974	}
975	if( d==-1 )
976	{
977	z.d = 1;
978	}
979	z.nw = nw;
980	if( d==-1 )
981	{
982	z.q = new double[n, nx+1+nx];
983	}
984	if( d==0 )
985	{
986	z.q = new double[n, nx+1];
987	}
988	if( d==1 )
989	{
990	z.q = new double[n, nx+1+nx];
991	}
992	if( d==2 )
993	{
994	z.q = new double[n, nx+1+nx+(int)Math.Round(nx(nx+1)0.5)];
995	}
996	z.tbuf = new int[nw];
997	z.rbuf = new double[nw];
998	z.xybuf = new double[nw, nx+1];
999	z.xbuf = new double[nx];
1000	}
1001
1002
1003	/*************************************************************************
1004	Linear least squares solver for small tasks.
1005
1006	Works faster than standard ALGLIB solver in non-degenerate cases (due to
1007	absense of internal allocations and optimized row/colums). In degenerate
1008	cases it calls standard solver, which results in small performance penalty
1009	associated with preliminary steps.
1010
1011	INPUT PARAMETERS:
1012	Y array[0..N-1]
1013	W array[0..N-1]
1014	FMatrix array[0..N-1,0..M], have additional column for temporary
1015	values
1016	Temp array[0..N]
1017	*************************************************************************/
1018	private static void idwinternalsolver(ref double[] y,
1019	ref double[] w,
1020	ref double[,] fmatrix,
1021	ref double[] temp,
1022	int n,
1023	int m,
1024	ref int info,
1025	ref double[] x,
1026	ref double taskrcond)
1027	{
1028	int i = 0;
1029	int j = 0;
1030	double v = 0;
1031	double tau = 0;
1032	double[] b = new double[0];
1033	densesolver.densesolverlsreport srep = new densesolver.densesolverlsreport();
1034	int i_ = 0;
1035	int i1_ = 0;
1036
1037
1038	//
1039	// set up info
1040	//
1041	info = 1;
1042
1043	//
1044	// prepare matrix
1045	//
1046	for(i=0; i<=n-1; i++)
1047	{
1048	fmatrix[i,m] = y[i];
1049	v = w[i];
1050	for(i_=0; i_<=m;i_++)
1051	{
1052	fmatrix[i,i_] = v*fmatrix[i,i_];
1053	}
1054	}
1055
1056	//
1057	// use either fast algorithm or general algorithm
1058	//
1059	if( m<=n )
1060	{
1061
1062	//
1063	// QR decomposition
1064	// We assume that M<=N (we would have called LSFit() otherwise)
1065	//
1066	for(i=0; i<=m-1; i++)
1067	{
1068	if( i<n-1 )
1069	{
1070	i1_ = (i) - (1);
1071	for(i_=1; i_<=n-i;i_++)
1072	{
1073	temp[i_] = fmatrix[i_+i1_,i];
1074	}
1075	reflections.generatereflection(ref temp, n-i, ref tau);
1076	fmatrix[i,i] = temp[1];
1077	temp[1] = 1;
1078	for(j=i+1; j<=m; j++)
1079	{
1080	i1_ = (1)-(i);
1081	v = 0.0;
1082	for(i_=i; i_<=n-1;i_++)
1083	{
1084	v += fmatrix[i_,j]*temp[i_+i1_];
1085	}
1086	v = tau*v;
1087	i1_ = (1) - (i);
1088	for(i_=i; i_<=n-1;i_++)
1089	{
1090	fmatrix[i_,j] = fmatrix[i_,j] - v*temp[i_+i1_];
1091	}
1092	}
1093	}
1094	}
1095
1096	//
1097	// Check condition number
1098	//
1099	taskrcond = rcond.rmatrixtrrcondinf(ref fmatrix, m, true, false);
1100
1101	//
1102	// use either fast algorithm for non-degenerate cases
1103	// or slow algorithm for degenerate cases
1104	//
1105	if( (double)(taskrcond)>(double)(10000nAP.Math.MachineEpsilon) )
1106	{
1107
1108	//
1109	// solve triangular system R*x = FMatrix[0:M-1,M]
1110	// using fast algorithm, then exit
1111	//
1112	x[m-1] = fmatrix[m-1,m]/fmatrix[m-1,m-1];
1113	for(i=m-2; i>=0; i--)
1114	{
1115	v = 0.0;
1116	for(i_=i+1; i_<=m-1;i_++)
1117	{
1118	v += fmatrix[i,i_]*x[i_];
1119	}
1120	x[i] = (fmatrix[i,m]-v)/fmatrix[i,i];
1121	}
1122	}
1123	else
1124	{
1125
1126	//
1127	// use more general algorithm
1128	//
1129	b = new double[m];
1130	for(i=0; i<=m-1; i++)
1131	{
1132	for(j=0; j<=i-1; j++)
1133	{
1134	fmatrix[i,j] = 0.0;
1135	}
1136	b[i] = fmatrix[i,m];
1137	}
1138	densesolver.rmatrixsolvels(ref fmatrix, m, m, ref b, 10000*AP.Math.MachineEpsilon, ref info, ref srep, ref x);
1139	}
1140	}
1141	else
1142	{
1143
1144	//
1145	// use more general algorithm
1146	//
1147	b = new double[n];
1148	for(i=0; i<=n-1; i++)
1149	{
1150	b[i] = fmatrix[i,m];
1151	}
1152	densesolver.rmatrixsolvels(ref fmatrix, n, m, ref b, 10000*AP.Math.MachineEpsilon, ref info, ref srep, ref x);
1153	taskrcond = srep.r2;
1154	}
1155	}
1156	}
1157	}

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

Download in other formats:

Update cookies preferences