Index: /issm/trunk-jpl/src/c/Container/Observations.cpp =================================================================== --- /issm/trunk-jpl/src/c/Container/Observations.cpp (revision 12288) +++ /issm/trunk-jpl/src/c/Container/Observations.cpp (revision 12289) @@ -162,9 +162,8 @@ /*Output and Intermediaries*/ - int nobs,i,j,k; - double range; + int i,j,k; + double distance; Observation *observation1 = NULL; Observation *observation2 = NULL; - int *indices = NULL; int *counter= (int*)xmalloc(n*sizeof(int)); @@ -175,23 +174,23 @@ observation1=(Observation*)this->GetObjectByOffset(i); - for(j=0;j=0.); + for(j=i+1;jSize();j++){ + observation2=(Observation*)this->GetObjectByOffset(j); - /*Find all observations that are in range*/ - this->quadtree->RangeSearch(&indices,&nobs,observation1->x,observation1->y,range); + distance=sqrt(pow(observation1->x - observation2->x,2.) + pow(observation1->y - observation2->y,2.)); + if(distance>x[n-1]) continue; - for (k=0;kGetObjectByOffset(indices[k]); - gamma[j] += 1./2.*pow(observation1->value - observation2->value,2.); - } + int index = int(distance/(x[1]-x[0])); + if(index>n-1) index = n-1; + if(index<0) index = 0; - counter[j] += nobs; - xfree((void**)&indices); + gamma[index] += 1./2.*pow(observation1->value - observation2->value,2.); + counter[index] += 1; } } /*Normalize semivariogram*/ - for(j=0;jObservationList(&x,&y,&obs,&n_obs,x_interp[idx],y_interp[idx],range); + if(n_obs==0){ + predictions[idx] = -999.; + error[idx] = -999.; + continue; + } /*Allocate intermediary matrix and vectors*/ + Gamma = (double*)xmalloc(n_obs*n_obs*sizeof(double)); gamma0 = (double*)xmalloc(n_obs*sizeof(double)); @@ -147,5 +153,6 @@ for(i=0;iSemiVariogram(x[i]-x[j],y[i]-y[j]); + //Gamma[i*n_obs+j] = variogram->SemiVariogram(x[i]-x[j],y[i]-y[j]); + Gamma[i*n_obs+j] = variogram->Covariance(x[i]-x[j],y[i]-y[j]); Gamma[j*n_obs+i] = Gamma[i*n_obs+j]; } @@ -154,5 +161,6 @@ /*Get semivariogram vector associated to this location*/ - for(i=0;iSemiVariogram(x[i]-x_interp[idx],y[i]-y_interp[idx]); + //for(i=0;iSemiVariogram(x[i]-x_interp[idx],y[i]-y_interp[idx]); + for(i=0;iCovariance(x[i]-x_interp[idx],y[i]-y_interp[idx]); /*Solve the three linear systems*/ Index: /issm/trunk-jpl/src/c/objects/Kriging/ExponentialVariogram.cpp =================================================================== --- /issm/trunk-jpl/src/c/objects/Kriging/ExponentialVariogram.cpp (revision 12288) +++ /issm/trunk-jpl/src/c/objects/Kriging/ExponentialVariogram.cpp (revision 12289) @@ -59,4 +59,22 @@ /*Variogram function*/ +/*FUNCTION ExponentialVariogram::Covariance{{{1*/ +double ExponentialVariogram::Covariance(double deltax,double deltay){ + /*The covariance can be deduced from the variogram from the following + * relationship: + * 2 gamma = C(x,x) + C(y,y) -2 C(x,y) + * so + * C(h) = sill - gamma */ + double h,a,cova; + + /*Calculate length*/ + h=sqrt(pow(deltax,2.)+pow(deltay,2.)); + + /*return covariance*/ + a = 1./3.; + cova = (sill-nugget)*exp(-h/(a*range)); + return cova; +} +/*}}}*/ /*FUNCTION ExponentialVariogram::SemiVariogram{{{1*/ double ExponentialVariogram::SemiVariogram(double deltax,double deltay){ Index: /issm/trunk-jpl/src/c/objects/Kriging/ExponentialVariogram.h =================================================================== --- /issm/trunk-jpl/src/c/objects/Kriging/ExponentialVariogram.h (revision 12288) +++ /issm/trunk-jpl/src/c/objects/Kriging/ExponentialVariogram.h (revision 12289) @@ -31,4 +31,5 @@ /*Variogram functions*/ double SemiVariogram(double deltax,double deltay); + double Covariance(double deltax,double deltay); }; #endif /* _EXPONENTIALVARIOGRAM_H */ Index: /issm/trunk-jpl/src/c/objects/Kriging/GaussianVariogram.cpp =================================================================== --- /issm/trunk-jpl/src/c/objects/Kriging/GaussianVariogram.cpp (revision 12288) +++ /issm/trunk-jpl/src/c/objects/Kriging/GaussianVariogram.cpp (revision 12289) @@ -59,4 +59,25 @@ /*Variogram function*/ +/*FUNCTION GaussianVariogram::Covariance{{{1*/ +double GaussianVariogram::Covariance(double deltax,double deltay){ + /*The covariance can be deduced from the variogram from the following + * relationship: + * 2 gamma = C(x,x) + C(y,y) -2 C(x,y) + * so + * C(h) = sill - gamma */ + double h2,a,cova; + + /*Calculate length square*/ + h2=pow(deltax,2.)+pow(deltay,2.); + + /*return covariance*/ + a = 1./3.; + cova = (sill-nugget)*exp(-h2/(a*range*range)); + + //if(h2<0.0000001) cova = sill; + + return cova; +} +/*}}}*/ /*FUNCTION GaussianVariogram::SemiVariogram{{{1*/ double GaussianVariogram::SemiVariogram(double deltax,double deltay){ @@ -70,4 +91,7 @@ a = 1./3.; gamma = (sill-nugget)*(1.-exp(-h2/(a*range*range))) + nugget; + + //if(h2>1000*1000) printf("gamma = %g h= %g\n",gamma,sqrt(h2)); + printf("h = %g gamma = %g\n",sqrt(h2),gamma); return gamma; } Index: /issm/trunk-jpl/src/c/objects/Kriging/GaussianVariogram.h =================================================================== --- /issm/trunk-jpl/src/c/objects/Kriging/GaussianVariogram.h (revision 12288) +++ /issm/trunk-jpl/src/c/objects/Kriging/GaussianVariogram.h (revision 12289) @@ -31,4 +31,5 @@ /*Variogram functions*/ double SemiVariogram(double deltax,double deltay); + double Covariance(double deltax,double deltay); }; #endif /* _GAUSSIANVARIOGRAM_H */ Index: /issm/trunk-jpl/src/c/objects/Kriging/PowerVariogram.cpp =================================================================== --- /issm/trunk-jpl/src/c/objects/Kriging/PowerVariogram.cpp (revision 12288) +++ /issm/trunk-jpl/src/c/objects/Kriging/PowerVariogram.cpp (revision 12289) @@ -60,4 +60,22 @@ /*Variogram function*/ +/*FUNCTION PowerVariogram::Covariance{{{1*/ +double PowerVariogram::Covariance(double deltax,double deltay){ + /*The covariance can be deduced from the variogram from the following + * relationship: + * 2 gamma = C(x,x) + C(y,y) -2 C(x,y) + * so + * C(h) = sill - gamma */ + double h,cova; + + /*Calculate length square*/ + h=sqrt(pow(deltax,2.)+pow(deltay,2.)); + + /*return covariance*/ + cova = 9999. - this->slope*pow(h,this->power); + + return cova; +} +/*}}}*/ /*FUNCTION PowerVariogram::SemiVariogram{{{1*/ double PowerVariogram::SemiVariogram(double deltax,double deltay){ @@ -71,4 +89,5 @@ gamma = this->nugget + this->slope*pow(h,this->power); + //if(h>1000) printf("gamma = %g h=%g\n",gamma,h); return gamma; } Index: /issm/trunk-jpl/src/c/objects/Kriging/PowerVariogram.h =================================================================== --- /issm/trunk-jpl/src/c/objects/Kriging/PowerVariogram.h (revision 12288) +++ /issm/trunk-jpl/src/c/objects/Kriging/PowerVariogram.h (revision 12289) @@ -31,4 +31,5 @@ /*Variogram functions*/ double SemiVariogram(double deltax,double deltay); + double Covariance(double deltax,double deltay); }; #endif /* _POWERVARIOGRAM_H */ Index: /issm/trunk-jpl/src/c/objects/Kriging/SphericalVariogram.cpp =================================================================== --- /issm/trunk-jpl/src/c/objects/Kriging/SphericalVariogram.cpp (revision 12288) +++ /issm/trunk-jpl/src/c/objects/Kriging/SphericalVariogram.cpp (revision 12289) @@ -59,4 +59,25 @@ /*Variogram function*/ +/*FUNCTION SphericalVariogram::Covariance{{{1*/ +double SphericalVariogram::Covariance(double deltax,double deltay){ + /*The covariance can be deduced from the variogram from the following + * relationship: + * 2 gamma = C(x,x) + C(y,y) -2 C(x,y) + * so + * C(h) = sill - gamma */ + double h,cova; + + /*Calculate length square*/ + h=sqrt(pow(deltax,2.)+pow(deltay,2.)); + + /*return covariance*/ + if(h<=range) + cova = (sill-nugget)*(1 - (3*h)/(2*range) + pow(h,3.)/(2*pow(range,3.)) ); + else + cova = 0.; + + return cova; +} +/*}}}*/ /*FUNCTION SphericalVariogram::SemiVariogram{{{1*/ double SphericalVariogram::SemiVariogram(double deltax,double deltay){ Index: /issm/trunk-jpl/src/c/objects/Kriging/SphericalVariogram.h =================================================================== --- /issm/trunk-jpl/src/c/objects/Kriging/SphericalVariogram.h (revision 12288) +++ /issm/trunk-jpl/src/c/objects/Kriging/SphericalVariogram.h (revision 12289) @@ -31,4 +31,5 @@ /*Variogram functions*/ double SemiVariogram(double deltax,double deltay); + double Covariance(double deltax,double deltay); }; #endif /* _SPHERICALVARIOGRAM_H */ Index: /issm/trunk-jpl/src/c/objects/Kriging/Variogram.h =================================================================== --- /issm/trunk-jpl/src/c/objects/Kriging/Variogram.h (revision 12288) +++ /issm/trunk-jpl/src/c/objects/Kriging/Variogram.h (revision 12289) @@ -13,4 +13,5 @@ virtual ~Variogram(){}; virtual double SemiVariogram(double deltax,double deltay)=0; + virtual double Covariance(double deltax,double deltay)=0; }; Index: sm/trunk-jpl/src/c/shared/Sorting/Quicksort.c =================================================================== --- /issm/trunk-jpl/src/c/shared/Sorting/Quicksort.c (revision 12288) +++ (revision ) @@ -1,19 +1,0 @@ -void swap(int *a, int *b) { - int t=*a; *a=*b; *b=t; -} -void sort(int arr[], int beg, int end) { - - if (end > beg + 1) { - int piv = arr[beg], l = beg + 1, r = end; - while (l < r) - { - if (arr[l] <= piv) - l++; - else - swap(&arr[l], &arr[--r]); - } - swap(&arr[--l], &arr[beg]); - sort(arr, beg, l); - sort(arr, r, end); - } -} Index: /issm/trunk-jpl/src/c/shared/Sorting/binary_search.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Sorting/binary_search.cpp (revision 12288) +++ /issm/trunk-jpl/src/c/shared/Sorting/binary_search.cpp (revision 12289) @@ -11,9 +11,9 @@ #include -int binary_search(int* poffset,int target, int* sorted_integers,int num_integers){ +int binary_search(int* poffset,int target,int* sorted_integers,int num_integers){ /*output: */ int offset; //offset, if found - int found=0; //found=0 if target is not found, 1 otherwise. + int found=0; //found=0 if target is not found, 1 otherwise. /*intermediary: */ @@ -64,3 +64,2 @@ return found; } - Index: /issm/trunk-jpl/src/c/shared/Sorting/sorting.h =================================================================== --- /issm/trunk-jpl/src/c/shared/Sorting/sorting.h (revision 12288) +++ /issm/trunk-jpl/src/c/shared/Sorting/sorting.h (revision 12289) @@ -6,7 +6,5 @@ #define _SORTING_H_ -int binary_search(int* poffset,int target, int* sorted_integers,int num_integers); - +int binary_search(int* poffset,int target,int* sorted_integers,int num_integers); #endif //ifndef _SORTING_H_ -