detect_ignition/new_likelihood/make_spline.m

   1 function [p_like_spline,p_deriv_spline,n_deriv_spline] = make_spline(time_bounds,num_pts)
   2 %function [p_like_spline,p_deriv_spline,n_deriv_spline,n_like_spline] = make_spline(time_bounds,num_pts)
   3 %[like_spline,deriv_spline] = make_spline(time_bounds,num_pts)
   4 % Function returns a splines for the postitive and negative detections
   5 % log-likelihood an their derivatives
   6 % requires curve fitting toolbox
   7 % inputs:
   8 %   time_bounds - hours before and after fire arrival time
   9 %   num pts - number of pts to use in creating the splines, should be an
  10 %   even number
  11 % outputs:
  12 %   p_like_spline - spline which gives the log-likehood of POSITIVE detection as a
  13 %      function of seconds since fire arrival time
  14 %   p_deriv_spline - the derivative of p_like_spline, computed by center
  15 %      difference approximation
  16 %   n_like_spline - spline which gives the log-likehood of NEGATIVE detection as a
  17 %      function of seconds since fire arrival time
  18 %   n_deriv_spline - the derivative of n_like_spline, computed by center
  19 %      difference approximation
  20 %
  21
  22 domain_size = num_pts;
  23 time_start = -time_bounds; time_end = 1.5*time_bounds;
  24 time_strip = linspace(time_start,time_end,num_pts);
  25 ts2 = linspace(0,1,domain_size); %dummy dimension
  26 [x,y] = meshgrid(time_strip,ts2);
  27 times = zeros(domain_size);
  28 for i = 1:domain_size
  29     times(i,:) = time_strip;
  30 end
  31 %figure,quick_mesh(times)
  32 dp = detection_probability(times);
  33 d_strip = dp(domain_size/2,:);
  34 % figure,plot(time_strip,d_strip),title('detection probability'),
  35 % xlabel('time since fire arrival')
  36 % ylabel('probability of detection')
  37 % ylim([0 1]);
  38
  39 like = [];
  40 radius = 30; %size over which to "integrate"
  41 %weight = gauss_weight(radius); %not needed, we shift the curve up
  42 dx = 30; %in meters, this is grid spacing, dy = dx also, maybe compute this by passing
  43          % in some tign data
  44 % make distance matrix one time and then pass into the gaussin computation
  45 d_squared = zeros(radius*2+1,radius*2+1);
  46 [m,n] = size(d_squared);
  47 for i = 1:m
  48     for j = i:n
  49         d_squared(i,j) = i^2+j^2+2*(radius+1)^2-2*(radius+1)*(i+j);
  50         if i == j
  51             d_squared(i,j) = 0.5*d_squared(i,j);
  52         end
  53     end
  54 end
  55
  56 d_squared = d_squared+d_squared';
  57 d_squared = dx^2*d_squared;
  58 sig = 750/3; %1000 meters at sig 3
  59
  60 counter = 1;
  61 for i=radius+1:domain_size-radius
  62     %counter
  63     temp_prob = dp(domain_size/2-radius:domain_size/2+radius,i-radius:i+radius);
  64     %l =  compute_pixel_probability(500,i,heat,radius, weight, detection_probs );
  65     gauss = exp(-d_squared/(2*sig^2));
  66     prob = gauss.*temp_prob;
  67     l = log(sum(prob(:))/(2*pi*sig^2));
  68     like(counter) = l;
  69     counter = counter +1;
  70 end
  71
  72 short_time = time_strip(1,radius+1:end-radius);
  73 %move close to zero
  74 l2 = like-max(like(:))+max(like(:))/1000;
  75 l2 = like-0.9999*max(like(:));
  76
  77
  78 %negative dections likelihood
  79 n_l2 = log(1-exp(l2));
  80
  81
  82 % compute derivative
  83 l2_prime = zeros(1,length(short_time));
  84 n_l2_prime = l2_prime;
  85 h = short_time(2) - short_time(1);
  86 for i = 2:length(short_time)-1
  87     l2_prime(i) = (l2(i+1)-l2(i-1))/(2*h);
  88     n_l2_prime(i) = (n_l2(i+1)-n_l2(i-1))/(2*h);
  89 end
  90 l2_prime(1) = l2_prime(2);
  91 l2_prime(end)=l2_prime(end-1);
  92 n_l2_prime(1) = n_l2_prime(2);
  93 n_l2_prime(end)= n_l2_prime(end-1);
  94
  95 plots = 0;
  96 if plots
  97     figure,plot(short_time,l2),title('pixel log likelihood')
  98     xlabel('Hours since fire arrival')
  99     ylabel('Log likelihood of Pos detection')
 100     figure,plot(short_time,l2_prime),title('Pos Derviative')
 101     figure,plot(short_time,n_l2),title('neg detection log likelihood')
 102     xlabel('Hours since fire arrival')
 103     ylabel('Log likelihood of Negative detection')
 104     figure,plot(short_time,n_l2_prime), title('Neg derivative')
 105 end
 106
 107 %make spline from the data short_time, l2, l2_prime, converting to seconds
 108 [p_like_spline, fit1] = createFit(3600*short_time, l2);
 109 [p_deriv_spline, fit2] = createFit(3600*short_time, l2_prime);
 110 %[n_like_spline, fit3] = createFit(3600*short_time, n_l2);
 111 [n_deriv_spline, fit4] = createFit(3600*short_time, n_l2_prime);
 112
 113 % [p_like_spline, fit1] = createFit(short_time, l2);
 114 % [p_deriv_spline, fit2] = createFit(short_time, l2_prime);
 115 % [n_like_spline, fit3] = createFit(short_time, n_l2);
 116 % [n_deriv_spline, fit4] = createFit(short_time, n_l2_prime);
 117
 118
 119
 120 end
 121