Skip to content

Commit

Permalink
Merge branch 'bono' into dev
Browse files Browse the repository at this point in the history
  • Loading branch information
@cliu3
cliu3 committed Nov 20, 2015
2 parents ac4d6e7 + b34661b commit 9d88254
Show file tree
Hide file tree
Showing 5 changed files with 130 additions and 103 deletions.
9 changes: 7 additions & 2 deletions backfun/plot_likelihood.m
View file
Expand Up @@ -5,10 +5,14 @@ function plot_likelihood(fish_no,plot_mpt)

filename=['ObsLh',num2str(fish_no),'.mat'];
load(filename)
filename=['mpt',num2str(fish_no),'.mat'];
filename=['tagdata',num2str(fish_no),'.mat'];
load(filename)
days = mpt.time;
if plot_mpt==1
filename=['mpt',num2str(fish_no),'.mat'];
load(filename)
end

days = td.time_plot(td.d24);

global fvcom_tidaldb
load(fvcom_tidaldb)
Expand All @@ -20,6 +24,7 @@ function plot_likelihood(fish_no,plot_mpt)

ndays=size(ObsLh,1);


if plot_mpt==1
[mpt_x,mpt_y] = my_project(mpt.long,mpt.lat,'forward');
end
Expand Down
8 changes: 8 additions & 0 deletions datalik/likelihood_cliu.m
View file
Expand Up @@ -301,6 +301,14 @@ function likelihood_cliu(fish_no,path_to_tags,tagname)


ObsLh(i,:)=ObsLh_dep_total.*ObsLh_temp_total.*AttLh;

% release location treatment
if i==1
[xl,yl]=my_project(tag.release_lon,tag.release_lat,'forward');
dist_rl = ( (fvcom.x-xl).^2+(fvcom.y-yl).^2 ).^0.5;
[~,rel_idx] = find(dist_rl==min(dist_rl));
ObsLh(i,rel_idx) = 1;
end
%
% figure(2);clf;
% patch('Vertices',[fvcom.x,fvcom.y],'Faces',fvcom.tri,'Cdata',ObsLh(i,:),'edgecolor','none','facecolor','interp');
Expand Down
9 changes: 7 additions & 2 deletions filter/hmmfilter.m
View file
@@ -1,4 +1,4 @@
function [phi,normaliser,pred] = hmmfilter(s,db,td,LIK)
function [phi,normaliser,pred,isDtoosmall] = hmmfilter(s,db,td,LIK)
%HMMFILTER Perform the forward sweep of the filtering.
%
% This function is called by hmmgeolocate.m and likelihood.m
Expand All @@ -18,6 +18,7 @@
icalc = length(td.d24);
phi = zeros(row,col,icalc); % Density function
pred = phi;
isDtoosmall = 0;

% Posterior probability for initial position of fish
post = zeros(row,col); post(td.y0,td.x0) = 1; % Dirac delta
Expand Down Expand Up @@ -74,7 +75,11 @@
% pause(0.01)

%if sum(isnan(post(:))) ~= 0, error('NaN found in predicted distribution at day %i',i), end
if sum(post(:)) == 0, error('Zero distribution at day %i, s = [%f,%f]\nIf you are estimating D try changing the bounds (in hmmgeolocate)',i,s(1),s(2)), end
if sum(post(:)) == 0
warning('Zero distribution at day %i, s = [%f,%f]\nIf you are estimating D try changing the bounds (in hmmgeolocate)',i,s(1),s(2));
isDtoosmall=1;
return
end
% Store likelihood function to be optimized later for D
% The normalising constant corresponds to the conditional distribution
% of the depth given the previous measurements
Expand Down
193 changes: 101 additions & 92 deletions filter/hmmgeolocate1.m
View file
@@ -1,5 +1,5 @@
function hmmgeolocate(tagno,mode,viewres,Duser,ext)
%HMMGEOLOCATE Obtain geolocation by filtering preprocessed data
%HMMGEOLOCATE Obtain geolocation by filtering preprocessed data
% HMMGEOLOCATE(TAGNO,MODE,VIEWRES,DUSER,EXT)
%
% - TAGNO indentifier as string for the tag to geolocate.
Expand All @@ -9,7 +9,7 @@ function hmmgeolocate(tagno,mode,viewres,Duser,ext)
% - MODE number of behaviour modes to use (1 or 2).
% default is 2.
% - VIEWRES plots the marginal distributions consecutively
% in an animation when the geolocation has finished
% in an animation when the geolocation has finished
% successfully (by using the viewdistr function).
% default is 'on'.
% - DUSER user defined diffusivity vector e.g DUSER = [10 100].
Expand Down Expand Up @@ -82,8 +82,8 @@ function hmmgeolocate(tagno,mode,viewres,Duser,ext)
disp(sprintf('==Tag #%s==',tagno));

% Number of days to calculate forward
[row,col]=size(db.depth);
icalc = length(td.d24);
[row,col]=size(db.depth);
icalc = length(td.d24);
if mode == 1
td.behav = ones(1,icalc); disp('Using ONE behavioural mode!')
elseif mode == 2
Expand All @@ -99,87 +99,96 @@ function hmmgeolocate(tagno,mode,viewres,Duser,ext)

%% Set up parameters %%
% k = time step = 1 day
k=1;

%% Initial guess on D %%
if db.h <= 0, error('db.h is less than or equal to zero, ie. the database i degenerate!'), end
if ext
result.D = Duser; % Unit: km^2/d?gn % 2255 18/2
k=1;

while 1
%% Initial guess on D %%
if db.h <= 0, error('db.h is less than or equal to zero, ie. the database i degenerate!'), end
if ext
result.D = Duser; % Unit: km^2/d?gn % 2255 18/2

else
result.D = Duser; % Unit: km^2/d?gn % 2255 18/2
D2s = k/db.h^2; result.D2s = D2s;
sEst = result.D*D2s;
s = kstest(sEst,25); % Calc size of conv kernel, stop if large
result.D = s/D2s;
end
% [s I] = max(result.D)*D2s;
% unc = sqrt(2*s);
% ks = ceil(unc*10+1); ks = ks + mod(ks,2) + 1;
% if ks > 100,
% fac = (25/ks)^2;
% result.D = [result.D(I) result.D(I)]*fac;
% sEst = result.D*D2s;
% warning('Diffusivity too large, overwriting! new D= %f, %f',result.D(1),result.D(2))
% end
% s = max(result.D(unique(td.behav))*D2s);
unc = sqrt(2*s);
ks = ceil(unc*10+1);
ks = ks + mod(ks,2) + 1;

else
result.D = Duser; % Unit: km^2/d?gn % 2255 18/2
D2s = k/db.h^2; result.D2s = D2s;
sEst = result.D*D2s;
s = kstest(sEst,25); % Calc size of conv kernel, stop if large
result.D = s/D2s;
end
% [s I] = max(result.D)*D2s;
% unc = sqrt(2*s);
% ks = ceil(unc*10+1); ks = ks + mod(ks,2) + 1;
% if ks > 100,
% fac = (25/ks)^2;
% result.D = [result.D(I) result.D(I)]*fac;
% sEst = result.D*D2s;
% warning('Diffusivity too large, overwriting! new D= %f, %f',result.D(1),result.D(2))
% end
% s = max(result.D(unique(td.behav))*D2s);
unc = sqrt(2*s);
ks = ceil(unc*10+1);
ks = ks + mod(ks,2) + 1;

%% Plot likelihood function
nevals = 5; % Number of evaluation points on likelihood curve
%LB = [2 2]*D2s; UB = [225 225]*D2s;
LB = [2 2]*D2s; UB = [300 300]*D2s;

%LB = [.2443 .2443]/10000*D2s; UB = [32.8 32.8]/10000*D2s;
%LB = kstest(LB,7); UB = kstest(UB,81);

%plotlikelihoodfunction
%return

%% Find MLE using builtin fmincon (optimisation toolbox)
if estimate == 1
disp('Estimating D...'), tic
ds = 0.00001;

if mode == 1
[sEst,loglik] = fminbnd(@(s) likelihood(s,db,td,LIK),LB(1),UB(1), ...
optimset('TolX',1e-4,'MaxFunEvals',30,'Display','iter'));
sEst = [sEst sEst];
time_fmins = toc; result.D = sEst/(D2s);
hess=(likelihood(sEst(1)+ds,db,td,LIK)+likelihood(sEst(1)-ds,db,td,LIK)-2*loglik)/ds^2;
result.MLvar = [1/(hess*D2s^2) 1/(hess*D2s^2)];
disp(sprintf('\nMLE found!\nDhat = %1.4f,\t stdev = %f\ntime spent: %1.4f sec\n',result.D(1), sqrt(result.MLvar(1)),time_fmins))
%% Plot likelihood function
nevals = 5; % Number of evaluation points on likelihood curve
%LB = [2 2]*D2s; UB = [225 225]*D2s;
LB = [2 2]*D2s; UB = [300 300]*D2s;

%LB = [.2443 .2443]/10000*D2s; UB = [32.8 32.8]/10000*D2s;
%LB = kstest(LB,7); UB = kstest(UB,81);

%plotlikelihoodfunction
%return

%% Find MLE using builtin fmincon (optimisation toolbox)
if estimate == 1
disp('Estimating D...'), tic
ds = 0.00001;

if mode == 1
[sEst,loglik] = fminbnd(@(s) likelihood(s,db,td,LIK),LB(1),UB(1), ...
optimset('TolX',1e-4,'MaxFunEvals',30,'Display','iter'));
sEst = [sEst sEst];
time_fmins = toc; result.D = sEst/(D2s);
hess=(likelihood(sEst(1)+ds,db,td,LIK)+likelihood(sEst(1)-ds,db,td,LIK)-2*loglik)/ds^2;
result.MLvar = [1/(hess*D2s^2) 1/(hess*D2s^2)];
disp(sprintf('\nMLE found!\nDhat = %1.4f,\t stdev = %f\ntime spent: %1.4f sec\n',result.D(1), sqrt(result.MLvar(1)),time_fmins))
else
guess = result.D*D2s;
[sEst,loglik] = fminsearchbnd(@(s) likelihood(s,db,td,LIK),guess,LB,UB, ...
optimset('TolX',1e-4,'MaxFunEvals',30,'Display','iter'));
time_fmins = toc; result.D = sEst/(D2s);
hess1=(likelihood(sEst+[ds 0],db,td,LIK)+likelihood(sEst+[-ds 0],db,td,LIK)-2*loglik)/ds^2;
hess2=(likelihood(sEst+[0 ds],db,td,LIK)+likelihood(sEst+[0 -ds],db,td,LIK)-2*loglik)/ds^2;
result.MLvar(1) = 1/(hess1*D2s^2);
result.MLvar(2) = 1/(hess2*D2s^2);
disp(sprintf('\nMLE found!\nDhat1 = %1.4f,\t stdev1 = %f\nDhat2 = %1.4f,\t stdev2 = %f\ntime spent: %1.4f sec\n',result.D(1), sqrt(result.MLvar(1)),result.D(2), sqrt(result.MLvar(2)),time_fmins))
end
result.loglikval = loglik;
end

if result.D(1)>result.D(2)
fprintf('\nD(1)>D(2), MLE failed... \nUsing default D\n');
result.D=guess./D2s;
end
%% Prediction %%
disp(sprintf('Using D: [%f %f]',result.D(1),result.D(2))),
if ext == 0
disp('Predicting...'),
%tic, [result.phi,normaliser] = hmmfilter(sEst,db,td,LIK); tt=toc;
tic, [result.phi,normaliser,result.pred,isDtoosmall] = hmmfilter(sEst,db,td,LIK); tt=toc;
else
guess = result.D*D2s;
[sEst,loglik] = fminsearchbnd(@(s) likelihood(s,db,td,LIK),guess,LB,UB, ...
optimset('TolX',1e-4,'MaxFunEvals',30,'Display','iter'));
time_fmins = toc; result.D = sEst/(D2s);
hess1=(likelihood(sEst+[ds 0],db,td,LIK)+likelihood(sEst+[-ds 0],db,td,LIK)-2*loglik)/ds^2;
hess2=(likelihood(sEst+[0 ds],db,td,LIK)+likelihood(sEst+[0 -ds],db,td,LIK)-2*loglik)/ds^2;
result.MLvar(1) = 1/(hess1*D2s^2);
result.MLvar(2) = 1/(hess2*D2s^2);
disp(sprintf('\nMLE found!\nDhat1 = %1.4f,\t stdev1 = %f\nDhat2 = %1.4f,\t stdev2 = %f\ntime spent: %1.4f sec\n',result.D(1), sqrt(result.MLvar(1)),result.D(2), sqrt(result.MLvar(2)),time_fmins))
disp('Predicting (extended version)...')
par.s = sEst;
par.mp = mp;
tic, [result.phi,normaliser,result.pred] = hmmfiltermode(par,td,LIK); tt=toc;
end

if isDtoosmall==0
break
else
Duser = Duser.*2;
disp('D is too small, doubled D and trying again...')
end
result.loglikval = loglik;
end

if result.D(1)>result.D(2)
fprintf('\nD(1)>D(2), MLE failed... \nUsing default D\n');
result.D=guess./D2s;
end
%% Prediction %%
disp(sprintf('Using D: [%f %f]',result.D(1),result.D(2))),
if ext == 0
disp('Predicting...'),
%tic, [result.phi,normaliser] = hmmfilter(sEst,db,td,LIK); tt=toc;
tic, [result.phi,normaliser,result.pred] = hmmfilter(sEst,db,td,LIK); tt=toc;
else
disp('Predicting (extended version)...')
par.s = sEst;
par.mp = mp;
tic, [result.phi,normaliser,result.pred] = hmmfiltermode(par,td,LIK); tt=toc;
end
disp(sprintf('\b done in %3.2f sec!',tt))

Expand Down Expand Up @@ -211,32 +220,32 @@ function hmmgeolocate(tagno,mode,viewres,Duser,ext)
%% Creating *.mat file
filename = sprintf('result%s',td.tagno);
disp(sprintf('Saving -> %s.mat <- in\n%s',filename,cd))
try
save(filename,'result')
catch EM
%try
% save(filename,'result')
%catch EM
save(filename,'result','-v7.3'); %for large bathymetric databases / long tags struct is too large for v7 datatype
end
%end
disp(sprintf('\nDONE geolocating!\n'))

%% View smoothed distribution
%if ~strcmp(viewres,'off'), close all, viewdistr(result.smooth_plot); end
if ~strcmp(viewres,'off'), close all,
if ~strcmp(viewres,'off'), close all,
figure, set(gcf,'position',[50 150 900 400])
subplot(121), viewdistr(result.UD_plot), title('Utilisation Distribution')
%subplot(121), viewdistr(result.UD_plot,[],[],[],[],'fancylock',result.land), title('Utilisation Distribution')
subplot(122), viewdistr(result.smooth_plot);
%subplot(122), viewdistr(result.smooth,[],[],[],[],'fancylock',result.land);
subplot(122), viewdistr(result.smooth_plot);
%subplot(122), viewdistr(result.smooth,[],[],[],[],'fancylock',result.land);
end

function sEst=kstest(s,siz)
sEst = s;
[s I] = max(s);
unc = sqrt(2*s);
unc = sqrt(2*s);
ks = ceil(unc*10+1); ks = ks + mod(ks,2) + 1;
if ks > 100,
if ks > 100,
fac = (siz/ks)^2;
sEst = [s s]*fac
warning('Diffusivity too large, overwriting!')
warning('Diffusivity too large, overwriting!')
end


14 changes: 7 additions & 7 deletions tbworkdir/run_tag.m
View file
@@ -1,7 +1,7 @@
%clear all;
%close all;
addpath(genpath('../'));
addpath(genpath('../../preprocess/'));
%addpath(genpath('../../preprocess/'));
%addpath('/opt/matlab/googleearth');

global fvcom_tidaldb % path to fvcom tidal database
Expand All @@ -21,7 +21,7 @@
tag_temp_accu = 0.1; % in degree C

%ptags = [12,22,24,55,56];
ptags = 7;
ptags = 11;

tag_num_range = ptags;

Expand All @@ -44,11 +44,11 @@

%do_parts = 6;

do_parts(1) = 0; %1 generate a new tidaldb, =0 use tidaldb.mat
do_parts(2) = 0; %2 strip
do_parts(3) = 0; %3 behavior
do_parts(4) = 0; %4 likelihood
do_parts(5) = 0; %5 cliu likelihood & tidal constraint
do_parts(1) = 1; %1 generate a new tidaldb, =0 use tidaldb.mat
do_parts(2) = 1; %2 strip
do_parts(3) = 1; %3 behavior
do_parts(4) = 1; %4 likelihood
do_parts(5) = 1; %5 cliu likelihood & tidal constraint
do_parts(6) = 1; %6 geolocate
do_parts(7) = 1; %7 most probable track
do_parts(8) = 0; %8 make a movie
Expand Down

0 comments on commit 9d88254

Please sign in to comment.