
clear 
clf

% Hypo2002
% a programme to use probabilistic inverse theory 
% to determine hypocenter locations (highly simplified)
% here simplified to a 2D problem, but extendable to 3D

% model and data
ns=4;								% number of stations (and observations)
sx=[5 10 25 30]; 	    			% station coordinates x (km)
sy=[ 0 0 0 0];  	                    % station coordinates y (km)
sz=[ 0 0 0 0];	                    % station coordinates z (km)
% sx=[0 10 12 15 18  20 22 25]; 	    			% station coordinates x (km)
% sy=[0 0 0 0 0 0 0 0];  	                    % station coordinates y (km)
% sz=[3 0 0 3 0 0 0 3];	                    % station coordinates z (km)
ns=length(sx);						% number of stations (and observations)

% Define uncertainties in data using variances (length of vector is number
% of data points
sobs=[ .2 .3 .2 .2];
%sobs=[ .5 .5 1 5 1 .5 .5 .5  ];

if length(sobs) ~= ns, disp(' Please adjust vector length (sobs) '), pause, end
%sobs=[ 0.5 1 0.5 0.5]; 

% TRUE MODEL
v=5.;								% homogeneous velocity (km/s)
xs=20;                              % source location x
ys=0;                               % source location y
zs=8;                               % source location z
t0=20;                              % origin time

% Calculate exact data
for i=1:ns,
    tobs(i)=t0+1/v*sqrt( (sx(i)-xs)^2 + (sy(i)-ys)^2 + (sz(i)-zs)^2);
end
t=tobs*0;         % vector used later



% Define model space and increments for model exploration
xa=10; xe=35; xinc=1.;			% upper and lower bounds in x0
ya=0; ye=0; yinc=.1;            % y0
za=0; ze=15; zinc=1.;           % z0
Ta=15; Te=25; Tinc=.25;         % t0
va=4; ve=6; vinc=.25;        % v
nx=(xe-xa)/xinc+1;              % number of sample in dimension x0
ny=(ye-ya)/yinc+1;              % y0
nz=(ze-za)/zinc+1;              % z0
nT=(Te-Ta)/Tinc+1;              % t0
nv=(ve-va)/vinc+1;              % v
sv=1; v0=4.5;

% Now let s solve the inverse problem
% Some parameters
ntot=nv*nx*ny*nz*nT;               % total number of forward modellings to do
nm=round(ntot/100);                    % info after ntot/100 samples 
nn=0;

disp(sprintf(' Proceed to calculate a posteriori pdf for %i x %i x %i x %i x %i = %i samples ? ',nv,nx,ny,nz,nT,ntot))
disp(' Press any key to continue ')
pause
disp(' Starting exploration of model space ... ')



% We now loop through the model space (x,y,z, and T)

sigD=zeros([nv nx ny nz nT]);      % initilization of pdf to 0

% Loop through model space

iv=0;
for v=va:vinc:ve,
iv=iv+1;

ix=0;
for x=xa:xinc:xe,
   ix=ix+1;
   iy=0;
 %  for y=ya:yinc:ye,
 %     iy=iy+1;
 			iz=0;
 			iy=1;
			 y=0;
      for z=za:zinc:ze,
         iz=iz+1;
         iT=0;
        for T=Ta:Tinc:Te,
            iT=iT+1;
            % calculate forward problem
            
            for is=1:ns,    % loop over number of observables
            t(is)=T+1/v*sqrt((x-sx(is))^2+(y-sy(is))^2+(z-sz(is))^2);
            end
        
            % calculate the associated probability
            
                term=1;
            for is=1:ns,    % loop over number of observables
                term=term*exp(-1/(2*sobs(is)^2)*(t(is)-tobs(is))^2);
            end
            sigM(iv,ix,iy,iz,iT)=term*exp(-1/(2*sv^2)*(v-v0)^2);         % 
          
                                      nn=nn+1;
                                      if rem(nn,nm)==0, 
                                         disp(sprintf(' Calculated %g %%',round(nn/ntot*100)))
                                      end
                                      
                          end
    %                 end
                  end
            end
            
        end
               
               
               
% Calculate the marginal a posteriori probabilities for plotting

% plane x-z
sigMxz=zeros([ix iz]);
for i=1:ix,
   for j=1:iz,
      sigMxz(i,j)=sum(sum(sum(sigM(:,i,:,j,:))));
      end
end

% normalize marginals for plotting
sigMxz=sigMxz/max(max(sigMxz));
x=xa:xinc:xe;
z=za:zinc:ze;

%  marginal a posteriori for T
for i=1:nT,sigMT(i)=sum(sum(sum(sum(sigM(:,:,:,:,i))))); end
s=sum(sigMT);
sigMT=sigMT/s;

%  prior for T
for i=1:nT,roMT(i)=1; end
s=sum(roMT);
roMT=roMT/s;

%  marginal a posteriori for v
for i=1:nv,sigMv(i)=sum(sum(sum(sum(sigM(i,:,:,:,:))))); end
s=sum(sigMv);
sigMv=sigMv/s;

% prior in v
for i=1:nv,roMv(i)=exp(-1/(2*sv^2)*((va+(i-1)*vinc)-v0)^2); end
s=sum(roMv);
roMv=roMv/s;

%  marginal a posteriori for x
for i=1:nx,sigMx(i)=sum(sum(sum(sum(sigM(:,i,:,:,:))))); end
s=sum(sigMx);
sigMx=sigMx/s;

%  prior for x
for i=1:nx,roMx(i)=1; end
s=sum(roMx);
roMx=roMx/s;

%  marginal a posteriori for z
for i=1:nz,sigMz(i)=sum(sum(sum(sum(sigM(:,:,:,i,:))))); end
s=sum(sigMz);
sigMz=sigMz/s;

%  prior for z
for i=1:nz,roMz(i)=1; end
s=sum(roMz);
roMz=roMz/s;


subplot(231)
imagesc(interp(x,4),interp(z,4),interp2(sigMxz',4)), hold on
xlabel(' x(km) ')
ylabel(' z(km) ')
title(' A posteriori probability: hypocenter ')

for i=1:ns
    plot(sx(i),sz(i),'r+'), hold on
end
hold off

subplot(232)
plot(Ta:Tinc:Te,sigMT,'k-',Ta:Tinc:Te,roMT,'r-')
xlabel(' Origin time (s) ')
ylabel(' Probability')
title(' A posteriori (black) - prior (red) probability: origin time ')

subplot(233)
plot(va:vinc:ve,sigMv,'k-',va:vinc:ve,roMv,'r-')
xlabel(' Velocity (km/s) ')
ylabel(' Probability')
title(' A posteriori (black) - prior (red) probability: velocity ')

subplot(234)
plot(xa:xinc:xe,sigMx,'k-',xa:xinc:xe,roMx,'r-')
xlabel(' Source location x (km) ')
ylabel(' Probability')
title(' A posteriori (black) - prior (red) probability: x ')

subplot(235)
plot(za:zinc:ze,sigMz,'k-',za:zinc:ze,roMz,'r-')
xlabel(' Source location z (km) ')
ylabel(' Probability')
title(' A posteriori (black) - prior (red) probability: z ')