00001
00006 #include "variogram.h"
00007 #include "geostat.h"
00008 #include <math.h>
00009 #include <map>
00010 #include <gsl/gsl_blas.h>
00011 #include <gsl/gsl_linalg.h>
00012 #include <algorithm>
00013 #include "matrix_ex.h"
00014
00022 double isoDist(const gsl_vector *c1, const gsl_vector *c2)
00023 {
00024 double norm;
00025 gsl_vector *tmp;
00026
00027 tmp = gsl_vector_alloc(c1->size);
00028 gsl_vector_memcpy(tmp, c1);
00029 gsl_vector_sub(tmp, c2);
00030 norm = gsl_blas_dnrm2(tmp);
00031 gsl_vector_free(tmp);
00032
00033 return norm;
00034 }
00035
00040 CGeostat::CGeostat(void)
00041 {
00042 m_coords.clear();
00043 m_pData = NULL;
00044 m_pTrend = NULL;
00045 m_domain.clear();
00046 m_pVarioCloud = NULL;
00047 m_pVariogram = NULL;
00048 m_pOKSys = NULL;
00049 }
00050
00055 CGeostat::~CGeostat(void)
00056 {
00057 destroy();
00058 }
00059
00068 bool CGeostat::initialize(size_t smpl, size_t dim)
00069 {
00070 destroy();
00071 return allocate(smpl, dim);
00072 }
00073
00082 bool CGeostat::allocate(size_t smpl, size_t dim)
00083 {
00084 destroy();
00085
00086 for (size_t s = 0;s < smpl;s++)
00087 m_coords.push_back(gsl_vector_alloc(dim));
00088
00089 for (size_t d = 0;d < dim;d++)
00090 m_domain.push_back( std::pair<double, double>(0, 0) );
00091
00092 m_pData = gsl_vector_alloc(smpl);
00093 m_pTrend = gsl_vector_calloc(dim + 1);
00094 m_pVarioCloud = new CVariogram;
00095 m_pVariogram = new CVariogram;
00096
00097 m_pOKSys = gsl_matrix_alloc(samples() + 1, samples() + 1);
00098
00099 return true;
00100 }
00101
00106 void CGeostat::destroy(void)
00107 {
00108 while (!m_coords.empty())
00109 {
00110 delete m_coords.back();
00111 m_coords.pop_back();
00112 }
00113 if (m_pData != NULL)
00114 {
00115 delete m_pData;
00116 m_pData = NULL;
00117 }
00118 if (m_pTrend != NULL)
00119 {
00120 delete m_pTrend;
00121 m_pTrend = NULL;
00122 }
00123 if (m_pVariogram != NULL)
00124 {
00125 delete m_pVariogram;
00126 m_pVariogram = NULL;
00127 }
00128 if (m_pVarioCloud != NULL)
00129 {
00130 delete m_pVarioCloud;
00131 m_pVarioCloud = NULL;
00132 }
00133 if (m_pOKSys != NULL)
00134 {
00135 delete m_pOKSys;
00136 m_pOKSys = NULL;
00137 }
00138 m_domain.clear();
00139 }
00140
00149 bool CGeostat::getCoordinate(gsl_vector *c, size_t s) const
00150 {
00151 if (!isActive() || s >= samples())
00152 return false;
00153 gsl_vector_memcpy(c, m_coords[s]);
00154
00155 return true;
00156 }
00157
00164 double CGeostat::getData(size_t s) const
00165 {
00166 if (!isActive() || s >= samples())
00167 return 0.0;
00168
00169 gsl_vector *c = gsl_vector_alloc(dimension());
00170 getCoordinate(c, s);
00171 double d = getResidual(s) + getTrend(c);
00172 gsl_vector_free(c);
00173 return d;
00174 }
00175
00182 double CGeostat::getTrend(const gsl_vector *c) const
00183 {
00184 if (!isActive())
00185 return 0.0;
00186
00187 double trend = gsl_vector_get(m_pTrend, 0);
00188 for (size_t d = 0;d < dimension();d++)
00189 trend += gsl_vector_get(m_pTrend, d + 1) * gsl_vector_get(c, d);
00190
00191 return trend;
00192 }
00193
00200 double CGeostat::getResidual(size_t s) const
00201 {
00202 if (!isActive() || s >= samples())
00203 return 0.0;
00204
00205 return gsl_vector_get(m_pData, s);
00206 }
00207
00216 bool CGeostat::getDomain(gsl_vector *lower, gsl_vector *upper) const
00217 {
00218 if (!isActive())
00219 return false;
00220
00221 for (size_t d = 0;d < dimension();d++)
00222 {
00223 gsl_vector_set(lower, d, m_domain[d].first);
00224 gsl_vector_set(upper, d, m_domain[d].second);
00225 }
00226
00227 return true;
00228 }
00229
00238 bool CGeostat::setDomain(const gsl_vector *lower, const gsl_vector *upper)
00239 {
00240 if (!isActive())
00241 return false;
00242
00243 for (size_t d = 0;d < dimension();d++)
00244 {
00245 m_domain[d].first = gsl_vector_get(lower, d);
00246 m_domain[d].second = gsl_vector_get(upper, d);
00247 }
00248
00249 return true;
00250 }
00251
00260 bool CGeostat::setCoordinate(size_t s, const gsl_vector *c)
00261 {
00262 if (!isActive() || s >= samples())
00263 return false;
00264
00265 gsl_vector_memcpy(m_coords[s], c);
00266
00267 return true;
00268 }
00269
00278 bool CGeostat::setData(size_t s, double d)
00279 {
00280 if (!isActive() || s >= samples())
00281 return false;
00282
00283 gsl_vector_set(m_pData, s, d);
00284 return true;
00285 }
00286
00292 std::vector<CVariogram::VarioItm> CGeostat::computeVariogramCloud(void) const
00293 {
00294 std::vector<CVariogram::VarioItm> vecVario;
00295 CVariogram::VarioItm itm;
00296 gsl_vector *c1, *c2;
00297 double r1, r2;
00298
00299 if (!isActive())
00300 return vecVario;
00301
00302 c1 = gsl_vector_alloc(dimension());
00303 c2 = gsl_vector_alloc(dimension());
00304
00305 for (size_t i = 0;i < samples() - 1;i++)
00306 {
00307 getCoordinate(c1, i);
00308 r1 = getResidual(i);
00309
00310 for (size_t j = i + 1;j < samples();j++)
00311 {
00312 getCoordinate(c2, j);
00313 r2 = getResidual(j);
00314
00315 itm.distance = isoDist(c1, c2);
00316 itm.dissimilarity = (r1 - r2) * (r1 - r2) / 2.0;
00317
00318 vecVario.push_back(itm);
00319 }
00320 }
00321
00322 gsl_vector_free(c1);
00323 gsl_vector_free(c2);
00324
00325 return vecVario;
00326 }
00327
00335 std::vector<CVariogram::VarioItm> CGeostat::computeExperimentalVariogram(const std::vector<CVariogram::VarioItm> &vcloud, double step) const
00336 {
00337 std::vector<CVariogram::VarioItm> vmodel;
00338 std::map<int, CVariogram::VarioItm> variomap;
00339 CVariogram::VarioItm itm;
00340
00341 for (size_t i = 0;i < vcloud.size();i++)
00342 {
00343 int key = static_cast<int>(vcloud[i].distance / step);
00344 if (variomap.find(key) == variomap.end())
00345 {
00346 variomap[key].distance = 1.0;
00347 variomap[key].dissimilarity = vcloud[i].dissimilarity;
00348 }
00349 else
00350 {
00351 variomap[key].distance += 1.0;
00352 variomap[key].dissimilarity += vcloud[i].dissimilarity;
00353 }
00354 }
00355
00356 for (std::map<int, CVariogram::VarioItm>::const_iterator it = variomap.begin();it != variomap.end();++it)
00357 {
00358 itm.distance = it->first * step;
00359 itm.dissimilarity = it->second.dissimilarity / it->second.distance;
00360 vmodel.push_back(itm);
00361 }
00362 return vmodel;
00363 }
00364
00374 bool CGeostat::estimate(int model, double power, double step)
00375 {
00376 if (!isActive())
00377 return false;
00378
00379 std::vector<CVariogram::VarioItm> vcloud;
00380 double n, s, r;
00381
00382
00383 approximateTrendSurface();
00384
00385 vcloud = computeVariogramCloud();
00386 m_pVarioCloud->setSample(vcloud);
00387
00388 m_pVariogram->setSample(computeExperimentalVariogram(vcloud, step));
00389 n = 0;
00390 m_pVariogram->getSample(m_pVariogram->samples() / 2, r, s);
00391 m_pVariogram->estimateModel(model, n, s, r, power, m_pVariogram->maxDistance() / 2.0);
00392
00393 precomputeKrigingSystem();
00394
00395 return true;
00396 }
00397
00408 bool CGeostat::getModelParameters(double &sill, double &range, double &nugget, double &power) const
00409 {
00410 if (!isActive())
00411 return false;
00412
00413 sill = m_pVariogram->sill();
00414 range = m_pVariogram->range();
00415 nugget = m_pVariogram->nugget();
00416 power = m_pVariogram->power();
00417
00418 return false;
00419 }
00420
00427 bool CGeostat::precomputeKrigingSystem(void)
00428 {
00429 gsl_vector *c1, *c2;
00430 gsl_matrix *tmp;
00431
00432 c1 = gsl_vector_alloc(dimension());
00433 c2 = gsl_vector_alloc(dimension());
00434 tmp = gsl_matrix_alloc(samples() + 1, samples() + 1);
00435
00436 for (size_t s1 = 0;s1 < samples();s1++)
00437 {
00438 getCoordinate(c1, s1);
00439
00440 for (size_t s2 = s1;s2 < samples();s2++)
00441 {
00442 getCoordinate(c2, s2);
00443
00444 double dist;
00445 dist = isoDist(c1, c2);
00446
00447 gsl_matrix_set(tmp, s1, s2, m_pVariogram->getModelData(dist));
00448 gsl_matrix_set(tmp, s2, s1, gsl_matrix_get(tmp, s1, s2));
00449 }
00450 }
00451
00452 for (size_t s = 0;s < samples();s++)
00453 {
00454 gsl_matrix_set(tmp, s, samples(), 1.0);
00455 gsl_matrix_set(tmp, samples(), s, 1.0);
00456 }
00457 gsl_matrix_set(tmp, samples(), samples(), 0.0);
00458
00459 MatrixEx::matrix_inverse(m_pOKSys, tmp);
00460
00461 gsl_vector_free(c1);
00462 gsl_vector_free(c2);
00463 gsl_matrix_free(tmp);
00464
00465 return true;
00466 }
00467
00472 double CGeostat::maxDist() const
00473 {
00474 double dist =0.0;
00475
00476 for (size_t d = 0;d < dimension();d++)
00477 dist += pow(m_domain[d].second - m_domain[d].first, 2.0);
00478 return sqrt(dist);
00479 }
00480
00488 bool CGeostat::isDomain(const gsl_vector *c) const
00489 {
00490 for (size_t d = 0;d < dimension();d++)
00491 {
00492 double ct = gsl_vector_get(c, d);
00493 if (ct < m_domain[d].first || ct > m_domain[d].second)
00494 return false;
00495 }
00496 return true;
00497 }
00498
00507 bool CGeostat::getWeightVector(gsl_vector *weight, const gsl_vector *c) const
00508 {
00509 if (!isActive())
00510 return false;
00511
00512 gsl_vector *spos, *vvario;
00513
00514 spos = gsl_vector_alloc(dimension());
00515 vvario = gsl_vector_alloc(samples() + 1);
00516
00517 for (size_t s = 0;s < samples();s++)
00518 {
00519 getCoordinate(spos, s);
00520 double dist = isoDist(c, spos);
00521
00522 gsl_vector_set(vvario, s, m_pVariogram->getModelData(dist));
00523 }
00524 gsl_vector_set(vvario, samples(), 1.0);
00525
00526 gsl_blas_dgemv(CblasNoTrans, 1.0, m_pOKSys, vvario, 0.0, weight);
00527
00528
00529 gsl_vector_free(spos);
00530 gsl_vector_free(vvario);
00531
00532 return true;
00533 }
00534
00544 bool CGeostat::getPredictData(double &pred, double &var, const gsl_vector *c) const
00545 {
00546 if (!isActive())
00547 return false;
00548
00549 gsl_vector *spos;
00550 gsl_vector *vvario, *weight;
00551 spos = gsl_vector_alloc(dimension());
00552 vvario = gsl_vector_alloc(samples() + 1);
00553 weight = gsl_vector_alloc(samples() + 1);
00554
00555 for (size_t s = 0;s < samples();s++)
00556 {
00557 getCoordinate(spos, s);
00558 double dist = isoDist(c, spos);
00559
00560 gsl_vector_set(vvario, s, m_pVariogram->getModelData(dist));
00561 }
00562 gsl_vector_set(vvario, samples(), 1.0);
00563
00564 gsl_blas_dgemv(CblasNoTrans, 1.0, m_pOKSys, vvario, 0.0, weight);
00565
00566 pred = getTrend(c);
00567 for (size_t s = 0;s < samples();s++)
00568 pred += gsl_vector_get(weight, s) * getResidual(s);
00569
00570 gsl_vector_set(vvario, samples(), 0.0);
00571 gsl_blas_ddot(weight, vvario, &var);
00572
00573 gsl_vector_free(spos);
00574 gsl_vector_free(vvario);
00575 gsl_vector_free(weight);
00576
00577 return true;
00578 }
00579
00589 bool CGeostat::getPredictData(double &pred, const gsl_vector *c, const gsl_vector *weight) const
00590 {
00591 if (!isActive())
00592 return false;
00593
00594 pred = getTrend(c);
00595 for (size_t s = 0;s < samples();s++)
00596 pred += gsl_vector_get(weight, s) * getResidual(s);
00597
00598 return true;
00599 }
00600
00607 bool CGeostat::approximateTrendSurface(void)
00608 {
00609 gsl_matrix *mPh, *mPhI;
00610 gsl_vector *c;
00611
00612 mPh = gsl_matrix_alloc(dimension() + 1, samples());
00613 mPhI = gsl_matrix_alloc(samples(), dimension() + 1);
00614 c = gsl_vector_alloc(dimension());
00615
00616 for (size_t s = 0;s < samples();s++)
00617 {
00618 getCoordinate(c, s);
00619 gsl_matrix_set(mPh, 0, s, 1.0);
00620 for (size_t d = 0;d < dimension();d++)
00621 gsl_matrix_set(mPh, d + 1, s, gsl_vector_get(c, d));
00622 }
00623 MatrixEx::matrix_pseudo_inverse(mPhI, mPh);
00624 gsl_blas_dgemv(CblasTrans, 1.0, mPhI, m_pData, 0.0, m_pTrend);
00625
00626 for (size_t s = 0;s < samples();s++)
00627 {
00628 getCoordinate(c, s);
00629 *gsl_vector_ptr(m_pData, s) -= getTrend(c);
00630 }
00631
00632 gsl_vector_free(c);
00633 gsl_matrix_free(mPh);
00634 gsl_matrix_free(mPhI);
00635
00636 return true;
00637 }
00638
00646 bool CGeostat::removeRedundantSample(size_t smpl)
00647 {
00648 if (!isActive() || smpl >= samples())
00649 return false;
00650
00651 std::vector<gsl_vector *> newCoord;
00652 gsl_vector *pcoord, *newData;
00653
00654 newData = gsl_vector_alloc(samples() - 1);
00655
00656 size_t sp1 =0;
00657 for (size_t s = 0;s < samples();s++)
00658 {
00659 if (s == smpl)
00660 continue;
00661 gsl_vector_set(newData, sp1, gsl_vector_get(m_pData, s));
00662
00663 pcoord = gsl_vector_alloc(dimension());
00664 getCoordinate(pcoord, s);
00665 newCoord.push_back(pcoord);
00666
00667 sp1++;
00668 }
00669
00670 while (!m_coords.empty())
00671 {
00672 delete m_coords.back();
00673 m_coords.pop_back();
00674 }
00675 delete m_pData;
00676
00677 m_pData = newData;
00678 m_coords = newCoord;
00679
00680 precomputeKrigingSystem();
00681
00682 return true;
00683 }
00684
00694 bool CGeostat::getPredictionErrorExcluded(double &pred, double &var, size_t s) const
00695 {
00696 if (!isActive() || s >= samples())
00697 return false;
00698
00699 gsl_vector *c, *c1, *c2;
00700 gsl_vector *vvario, *weight;
00701 gsl_matrix *mGamma, *mGammaI;
00702 size_t sp1, sp2;
00703
00704 c = gsl_vector_alloc(dimension());
00705 c1 = gsl_vector_alloc(dimension());
00706 c2 = gsl_vector_alloc(dimension());
00707
00708 vvario = gsl_vector_alloc(samples());
00709 weight = gsl_vector_alloc(samples());
00710 mGamma = gsl_matrix_alloc(samples(), samples());
00711 mGammaI = gsl_matrix_alloc(samples(), samples());
00712
00713 getCoordinate(c, s);
00714
00715 sp1 = 0;
00716 for (size_t s1 = 0;s1 < samples();s1++)
00717 {
00718 if (s1 == s)
00719 continue;
00720
00721 getCoordinate(c1, s1);
00722
00723 sp2 = 0;
00724 for (size_t s2 = 0;s2 < samples();s2++)
00725 {
00726 if (s2 == s)
00727 continue;
00728
00729 getCoordinate(c2, s2);
00730 double dist = isoDist(c1, c2);
00731
00732 gsl_matrix_set(mGamma, sp1, sp2, m_pVariogram->getModelData(dist));
00733
00734 sp2++;
00735 }
00736 sp1++;
00737 }
00738
00739 for (size_t s = 0;s < samples() - 1;s++)
00740 {
00741 gsl_matrix_set(mGamma, s, samples() - 1, 1.0);
00742 gsl_matrix_set(mGamma, samples() - 1, s, 1.0);
00743 }
00744 gsl_matrix_set(mGamma, samples() - 1, samples() - 1, 0.0);
00745 MatrixEx::matrix_inverse(mGammaI, mGamma);
00746
00747
00748 sp1 = 0;
00749 for (size_t s = 0;s < samples();s++)
00750 {
00751 if (s == s)
00752 continue;
00753 getCoordinate(c1, s);
00754 double dist = isoDist(c, c1);
00755 gsl_vector_set(vvario, sp1++, m_pVariogram->getModelData(dist));
00756 }
00757 gsl_vector_set(vvario, samples() - 1, 1.0);
00758 gsl_blas_dgemv(CblasNoTrans, 1.0, mGamma, vvario, 0.0, weight);
00759
00760 pred = getTrend(c);
00761 sp1 = 0;
00762 for (size_t s = 0;s < samples();s++)
00763 {
00764 if (s == s)
00765 continue;
00766 pred += gsl_vector_get(weight, sp1) * getResidual(s);
00767 sp1++;
00768 }
00769 gsl_vector_set(vvario, samples() - 1, 0.0);
00770 gsl_blas_ddot(weight, vvario, &var);
00771
00772
00773 gsl_vector_free(c);
00774 gsl_vector_free(c1);
00775 gsl_vector_free(c2);
00776
00777 gsl_matrix_free(mGamma);
00778 gsl_matrix_free(mGammaI);
00779 gsl_vector_free(vvario);
00780 gsl_vector_free(weight);
00781
00782 return true;
00783 }
00784
00791 double CGeostat::selectEliminatee(size_t &candidate) const
00792 {
00793 if (!isActive())
00794 return 0;
00795
00796 gsl_vector *c, *c1, *c2;
00797 gsl_vector *vvario, *weight;
00798 gsl_matrix *mGamma, *mGammaR, *mGammaI;
00799 size_t sp1, sp2;
00800 double var, varmin;
00801
00802 c = gsl_vector_alloc(dimension());
00803 c1 = gsl_vector_alloc(dimension());
00804 c2 = gsl_vector_alloc(dimension());
00805
00806 mGamma = gsl_matrix_alloc(samples(), samples());
00807 mGammaR = gsl_matrix_alloc(samples(), samples());
00808 mGammaI = gsl_matrix_alloc(samples(), samples());
00809 vvario = gsl_vector_alloc(samples());
00810 weight = gsl_vector_alloc(samples());
00811
00812 for (size_t s1 = 0;s1 < samples();s1++)
00813 {
00814 getCoordinate(c1, s1);
00815 for (size_t s2 = s1;s2 < samples();s2++)
00816 {
00817 getCoordinate(c2, s2);
00818 double dist = isoDist(c1, c2);
00819 gsl_matrix_set(mGamma, s1, s2, m_pVariogram->getModelData(dist));
00820 gsl_matrix_set(mGamma, s2, s1, gsl_matrix_get(mGamma, s1, s2));
00821 }
00822 }
00823 for (size_t s = 0;s < samples() - 1;s++)
00824 {
00825 gsl_matrix_set(mGammaR, s, samples() - 1, 1.0);
00826 gsl_matrix_set(mGammaR, samples() - 1, s, 1.0);
00827 }
00828 gsl_matrix_set(mGammaR, samples() - 1, samples() - 1, 0.0);
00829
00830 varmin = 1.0e6;
00831 for (size_t smpl = 0;smpl < samples();smpl++)
00832 {
00833 getCoordinate(c, smpl);
00834
00835 sp1 = 0;
00836 for (size_t s1 = 0;s1 < samples();s1++)
00837 {
00838 if (s1 != smpl)
00839 {
00840 sp2 = 0;
00841 for (size_t s2 = 0;s2 < samples();s2++)
00842 {
00843 if (s2 != smpl)
00844 gsl_matrix_set(mGammaR, sp1, sp2++, gsl_matrix_get(mGamma, s1, s2));
00845 }
00846 sp1++;
00847 }
00848 }
00849 MatrixEx::matrix_inverse(mGammaI, mGammaR);
00850
00851
00852 sp1 = 0;
00853 for (size_t s = 0;s < samples();s++)
00854 {
00855 if (s != smpl)
00856 {
00857 getCoordinate(c1, s);
00858 double dist = isoDist(c, c1);
00859 gsl_vector_set(vvario, sp1++, m_pVariogram->getModelData(dist));
00860 }
00861 }
00862 gsl_vector_set(vvario, samples() - 1, 1.0);
00863 gsl_blas_dgemv(CblasNoTrans, 1.0, mGammaI, vvario, 0.0, weight);
00864
00865 gsl_vector_set(vvario, samples() - 1, 0.0);
00866 gsl_blas_ddot(weight, vvario, &var);
00867
00868 if (var < varmin)
00869 {
00870 candidate = smpl;
00871 varmin = var;
00872 }
00873 }
00874
00875 gsl_vector_free(c);
00876 gsl_vector_free(c1);
00877 gsl_vector_free(c2);
00878
00879 gsl_matrix_free(mGamma);
00880 gsl_matrix_free(mGammaR);
00881 gsl_matrix_free(mGammaI);
00882
00883 gsl_vector_free(vvario);
00884 gsl_vector_free(weight);
00885
00886
00887 return 1.0 - varmin / m_pVariogram->sill();
00888 }
00889
00898 double CGeostat::evalModelValidity(int model, double power, double step)
00899 {
00900 if (!isActive())
00901 return 0.0;
00902
00903 std::vector<CVariogram*> vecCloud, vecModel;
00904 std::vector<CVariogram::VarioItm> vcloud;
00905 double n =0, s =0, r =0;
00906
00907 approximateTrendSurface();
00908
00909 vcloud = computeVariogramCloud();
00910 m_pVarioCloud = new CVariogram;
00911 m_pVarioCloud->setSample(vcloud);
00912
00913 m_pVariogram = new CVariogram;
00914 m_pVariogram->setSample(computeExperimentalVariogram(vcloud, step));
00915 m_pVariogram->getSample(m_pVariogram->samples() / 2, r, s);
00916
00917 m_pVariogram->estimateModel(model, n, s, r, power, m_pVariogram->maxDistance() / 2.0);
00918
00919
00920 gsl_vector *c, *c1, *c2;
00921 gsl_vector *vvario, *weight;
00922 gsl_matrix *mGamma, *mGammaR, *mGammaI;
00923 size_t sp1, sp2;
00924 double pred, errsum;
00925
00926 c = gsl_vector_alloc(dimension());
00927 c1 = gsl_vector_alloc(dimension());
00928 c2 = gsl_vector_alloc(dimension());
00929
00930 mGamma = gsl_matrix_alloc(samples(), samples());
00931 mGammaR = gsl_matrix_alloc(samples(), samples());
00932 mGammaI = gsl_matrix_alloc(samples(), samples());
00933 vvario = gsl_vector_alloc(samples());
00934 weight = gsl_vector_alloc(samples());
00935
00936 for (size_t s1 = 0;s1 < samples();s1++)
00937 {
00938 getCoordinate(c1, s1);
00939 for (size_t s2 = s1;s2 < samples();s2++)
00940 {
00941 getCoordinate(c2, s2);
00942 double dist = isoDist(c1, c2);
00943 gsl_matrix_set(mGamma, s1, s2, m_pVariogram->getModelData(dist));
00944 gsl_matrix_set(mGamma, s2, s1, gsl_matrix_get(mGamma, s1, s2));
00945 }
00946 }
00947 for (size_t s = 0;s < samples() - 1;s++)
00948 {
00949 gsl_matrix_set(mGammaR, s, samples() - 1, 1.0);
00950 gsl_matrix_set(mGammaR, samples() - 1, s, 1.0);
00951 }
00952 gsl_matrix_set(mGammaR, samples() - 1, samples() - 1, 0.0);
00953
00954
00955 errsum = 0.0;
00956 for (size_t smpl = 0;smpl < samples();smpl++)
00957 {
00958 sp1 = 0;
00959 for (size_t s1 = 0;s1 < samples();s1++)
00960 {
00961 if (s1 != smpl)
00962 {
00963 sp2 = 0;
00964 for (size_t s2 = 0;s2 < samples();s2++)
00965 {
00966 if (s2 != smpl)
00967 {
00968 gsl_matrix_set(mGammaR, sp1, sp2++, gsl_matrix_get(mGamma, s1, s2));
00969 }
00970 }
00971 sp1++;
00972 }
00973 }
00974 MatrixEx::matrix_inverse(mGammaI, mGammaR);
00975
00976 getCoordinate(c, smpl);
00977 sp1 = 0;
00978 for (size_t s = 0;s < samples();s++)
00979 {
00980 if (s != smpl)
00981 {
00982 getCoordinate(c1, s);
00983 double dist = isoDist(c, c1);
00984 gsl_vector_set(vvario, sp1++, m_pVariogram->getModelData(dist));
00985 }
00986 }
00987 gsl_vector_set(vvario, samples() - 1, 1.0);
00988 gsl_blas_dgemv(CblasNoTrans, 1.0, mGammaI, vvario, 0.0, weight);
00989
00990 sp1 = 0;
00991 pred = 0.0;
00992 for (size_t s = 0;s < samples();s++)
00993 {
00994 if (s != smpl)
00995 pred += gsl_vector_get(weight, sp1++) * getResidual(s);
00996 }
00997 errsum += pow(pred - getResidual(smpl), 2.0);
00998 }
00999
01000 gsl_vector_free(c);
01001 gsl_vector_free(c1);
01002 gsl_vector_free(c2);
01003
01004 gsl_matrix_free(mGamma);
01005 gsl_matrix_free(mGammaR);
01006 gsl_matrix_free(mGammaI);
01007
01008 gsl_vector_free(vvario);
01009 gsl_vector_free(weight);
01010
01011 return errsum;
01012 }
