Kriging algorithms now suports 2D anisotropy, and ND is range covaria…

…nce models. (Nd kriging without rotation is now OK)

edist.m

@@ -1,90 +1,80 @@
 % edist : Euclidean distance
 %
-% Call : 
+% Call :
 %   D=edist(p1,p2,transform,isorange)
-% 
+%
 % p1,p2 : vectors
 %
-% transform : GSTAT anisotropy and/or range information 
+% transform : GSTAT anisotropy and/or range information
 %
-% isorange : [0] (default), transform is the usual GSTAT-anisotropy setting 
+% isorange : [0] (default), transform is the usual GSTAT-anisotropy setting
 % isorange : [1] means that transform simply lists the range in
-% each dimensions, and that no rotation is performed 
+% each dimensions, and that no rotation is performed
 
 function D=edist(p1,p2,transform,isorange)
-
-
-
-  if nargin<4
+
+if nargin<4
     isorange=0;
-  end
-  
-  if nargin==1;
+end
+
+if nargin==1;
     p2=0.*p1;
-  end
-
-  n_dim=size(p1,2);
-
-  if size(p1,1)==1
-    dp=(p1-p2)';
-  else
-    dp=(p1-p2);
-  end
-
-  if isorange==1
+end
+
+n_dim=size(p1,2);
+
+dp=(p1-p2);
+
+if isorange==1
+    %mgstat_verbose(sprintf('%s : isorange',mfilename))
     % ONLY SACLING, no transformation
 
-    %rescale=transform
-    %RescaleMat=eye(length(rescale));
-    %for i=1:length(rescale)
-    %  RescaleMat(i,i)=rescale(i);
-    %end
-    %dp=RescaleMat*dp;
+    if length(transform)==1;
+        transform=ones(1,n_dim).*transform;
+    end
 
     if transform==0
-      % Do not transfrom since this is likely a nugget
+        % Do not transfrom since this is likely a nugget
     else
-      dp=transform(1).*dp./transform(:);
+        for j=1:n_dim;
+            dp(:,j)=dp(:,j)./transform(j);
+        end
+        dp=transform(1).*dp;
     end
-    D=sqrt(dp'*dp);
-    return    
-  end
-
-
+else
 
-
-  % 2D COORDINATE TRANSFORMATION
-  if (nargin>2)&(length(p1)==2);
 
-    if length(transform)>1
-
-      rescale=transform(1:2);
-      rotate=transform(3);
-
-
-      RotMat=[cos(rotate) -sin(rotate);sin(rotate) cos(rotate)];
-      dp=RotMat*dp;
-
-      RescaleMat=eye(length(rescale));
-      for i=1:length(rescale)
-        RescaleMat(i,i)=rescale(i);
-      end
-      dp=RescaleMat*dp;
-    end
-  end
+    % 2D COORDINATE TRANSFORMATION
+    if (nargin>2)&(n_dim==2);
+        mgstat_verbose(sprintf('%s : 2D coordinate transform'),0);
+        if length(transform)>1
+            rescale=transform(1:2);
+            rotate=transform(3)*pi/180;
 
-  % 3D COORDINATE TRANSFORMATION
-  if (nargin>2)&(length(p1)==3);
-    % NOT YET IMPLEMENTED, SEE GSLIB BOOK
-  end
+            dp=dp';
 
+            RotMat=[cos(rotate) -sin(rotate);sin(rotate) cos(rotate)];
+            RescaleMat=eye(length(rescale));
+            RescaleMat(1,1)=1;
+            RescaleMat(2,2)=rescale(2);
+
+            %dp=RotMat*dp;%dp=RescaleMat*dp;
+            dp=RescaleMat*RotMat*dp;
+            dp=dp./rescale(2);
+
+            dp=dp';
+        end
+    end
+
+    % 3D COORDINATE TRANSFORMATION
+    if (nargin>2)&(n_dim==3);
+        mgstat_verbose(sprintf('%s : 3D anisotropy not yet implemented',mfilename),-1)
+        % NOT YET IMPLEMENTED, SEE GSLIB BOOK
+    end
+end
 
-  if n_dim==1
-      D=sqrt(dp.^2+dp.^2);
-  elseif size(p1,1)==2
-      D=sqrt(dp(1,:).^2+dp(2,:).^2);
-      %  elseif size(p1,1)==3
-      % D=sqrt(dp(1,:).^2+dp(2,:).^2+dp(3,:).^2);
-  else
-      D=sqrt(dp'*dp);
-  end
+if n_dim==1
+    D=sqrt(dp.^2+dp.^2);
+else
+    D=transpose(sqrt(sum(transpose(dp.^2))));
+end

krig.m

@@ -157,14 +157,14 @@
     mgstat_verbose('Kriging with a trend',1);
   elseif  (any(strcmp(fieldnames(options),'mean')) | any(strcmp(fieldnames(options),'sk_mean')))
     ktype=0; % SK
-    mgstat_verbose('Simple Kriging',-1);
+    mgstat_verbose('Simple Kriging',1);
   else 
     if size(val_known,1)==1
         ktype=0;
         mgstat_verbose('Forcing simple kriging (only one data point)',20);
     else
         ktype=1; % OK
-        mgstat_verbose('Ordinary Kriging',-1);
+        mgstat_verbose('Ordinary Kriging',1);
     end
   end
 
@@ -175,7 +175,7 @@
     mgstat_verbose('Warning : you called krig with more than one',10)
     mgstat_verbose('          unknown data location',10)
     mgstat_verbose(sprintf('%s --- Calling krig_npoint',mfilename),0)
-    [d_est,d_var,d2d,d2u]=krig_npoint(pos_known,val_known,pos_est,V,options);
+    [d_est,d_var]=krig_npoint(pos_known,val_known,pos_est,V,options);
     lambda_sk=[];K=[];k=[];inhood=[];
     return
   end
@@ -223,6 +223,7 @@
       catch
           keyboard
       end
+      % MAKE USE OF DIRECT CALL TO PRECAL_COV
     end    
     K=gvar-K;
   end
@@ -231,6 +232,7 @@
   for i=1:nknown
     K(i,i)=K(i,i)+unc_known(i);
   end
+
   % Data to Unknown matrix
   if any(strcmp(fieldnames(options),'d2u')); 
     k=options.d2u(inhood,:);
@@ -242,7 +244,12 @@
       for i=1:nknown;
         d(i)=edist(pos_known(i,:),pos_est(1,:),V(iV).par2,isorange);      
       end      
-      k=k+semivar_synth(V(iV),d);
+      try
+          V(iV).par2=V(iV).par2(1); % SCALE TO FIRST RANGE
+          k=k+semivar_synth(V(iV),d);
+      catch
+          keyboard
+      end
     end
     k=gvar-k;
   end

krig_covar_lik.m

@@ -34,8 +34,8 @@
 if isfield(options,'Cm');
     Cm=options.Cm;
 else
-    %Cm=precal_cov(pos_known,pos_known,V,options);
-    Cm=precal_cov(pos_known,pos_known,V);
+    Cm=precal_cov(pos_known,pos_known,V,options);
+    %Cm=precal_cov(pos_known,pos_known,V);
 end
 % TMH 05/05/2009 :
 % CM IS NOT RIGHT IN CASE OF A NUGGET
@@ -53,6 +53,9 @@
     sample_size=nknown;
     sill=sum([V.par1]);
 
+    Cd=eye(size(Cm))+0.00001*sill; % FOR BETTER PERFORMANCE
+    Cm=Cm+Cd;
+
     Q=Cm./sill;
     d_val=val_known(:,1)-mean(val_known(:,1));
     iQ=inv(Q);

krig_npoint.m

@@ -1,12 +1,12 @@
 % krig_npoint : as 'krig' butfor multiple estimation position.
 %
-% [d_est,d_var,d2d,d2u]=krig_npoint(pos_known,val_known,pos_est,V,options);
+% [d_est,d_var,options]=krig_npoint(pos_known,val_known,pos_est,V,options);
 %
 % As krig, but allowing size(pos_known,1)>1
 %
 % See also : krig
 % 
-function [d_est,d_var,d2d,d2u]=krig_npoint(pos_known,val_known,pos_est,V,options);
+function [d_est,d_var,options]=krig_npoint(pos_known,val_known,pos_est,V,options);
 
 if nargin<5
   options.null=0;
@@ -27,29 +27,32 @@
 
 %% %% %% BUG BUG IN SETTUP ING d2d table!!!
 
-% if isfield(options,'d2d')
-%   d2d=options.d2d;
-% elseif ~isfield(options,'noprecalc_d2d')
-%     if size(pos_known,2)>1; % SOMETHING WRONG WHEN USING 1D and options.d2d BUG
-%         d2d=precal_cov(pos_known,pos_known,V,options);
-%         options.d2d=d2d;
-%         mgstat_verbose(sprintf('%s : calulated d2d',mfilename),12)
-%     else        
-%         %%% BUG
-%         mgstat_verbose(sprintf('%s : COULD NOT calulated d2d for 1D data set',mfilename),12)
-%     end
-% end
+if isfield(options,'d2d')
+    d2d=options.d2d;
+elseif ~isfield(options,'noprecalc_d2d')
+    mgstat_verbose(sprintf('%s : precalculating data2data covariance matrix',mfilename),-1)
+    options.d2d=precal_cov(pos_known,pos_known,V,options);
+end
 
 % 
 if isfield(options,'d2u')
   d2u=options.d2u;
 elseif isfield(options,'precalc_d2u')
-  d2u=precal_cov(pos_known,pos_est,V);
-  options.d2u=d2u;
+    if options.precalc_d2u==1;
+    mgstat_verbose(sprintf('%s : precalculating data2unknown covariance matrix',mfilename),-1)
+    for j=1:size(pos_est,1)
+        if (j/500)==round(j/500), 
+            progress_txt(j,n_est,sprintf('%s : precal d2u',mfilename));
+        end
+        d2u(:,j)=precal_cov(pos_est(j,:),pos_known,V,options);
+    end
+    options.d2u=d2u;
+    end
 end
 d_est=zeros(n_est,1);
 d_var=zeros(n_est,1);
 
+
 %options=rmfield(options,'d2d');
 if isfield(options,'d2u');
   % MAKE USE OF PRECALCULATED DATA 2 UNKNOWN
@@ -62,8 +65,8 @@
   end
 else
   for i=1:n_est
-%    if (i/50)==round(i/50), 
-    if (i/5)==round(i/5), 
+    if (i/100)==round(i/100), 
+%    if (i/5)==round(i/5), 
       progress_txt(i,n_est,sprintf('%s : kriging',mfilename));
     end
     % SOMETHING WRONG WHEN USING 1D and options.d2d

krig_optim_ml.m

@@ -28,7 +28,7 @@
 % FIRST SAMPLE THE ATTRIBUTE SPACE
 %maxit=options.maxit;
 %options.maxit=10;
-figure(5);
+figure(5);clf;set_paper('landscape');
 [V,be,L,par2,nugfrac,Vall,options]=krig_optim_mcmc(pos_known,val_known,V,options);
 %options.maxit=maxit;
 iop=find(L==max(L));iop=iop(1);
@@ -44,7 +44,9 @@
 
 %% NEXT LINE TO SKIP LOCAL SEARCH
 %Vop2=Vop1;return
-figure(4);
+
+
+figure(4);clf;set_paper('landscape');
 % THEN LOOK FOR THE LOCAL MAX LIKELIHOOD MODEL
 options.step_nugfrac=options.step_nugfrac.*.01;
 options.step_range=options.step_range.*.01;