diff --git a/backfun/plot_likelihood.m b/backfun/plot_likelihood.m index 4ff8565..3b5c404 100644 --- a/backfun/plot_likelihood.m +++ b/backfun/plot_likelihood.m @@ -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) @@ -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 diff --git a/datalik/likelihood_cliu.m b/datalik/likelihood_cliu.m index 353cdab..5d5cbec 100644 --- a/datalik/likelihood_cliu.m +++ b/datalik/likelihood_cliu.m @@ -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'); diff --git a/filter/hmmfilter.m b/filter/hmmfilter.m index 562744e..04ed653 100755 --- a/filter/hmmfilter.m +++ b/filter/hmmfilter.m @@ -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 @@ -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 @@ -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 diff --git a/filter/hmmgeolocate1.m b/filter/hmmgeolocate1.m index 8b2e921..e6475a3 100755 --- a/filter/hmmgeolocate1.m +++ b/filter/hmmgeolocate1.m @@ -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. @@ -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]. @@ -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 @@ -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)) @@ -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 diff --git a/tbworkdir/run_tag.m b/tbworkdir/run_tag.m index 6fb7764..50258e3 100644 --- a/tbworkdir/run_tag.m +++ b/tbworkdir/run_tag.m @@ -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 @@ -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; @@ -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