detect_ignition/new_likelihood/test4.m

   1 % test 4
   2
   3 domain_size = 2000;
   4 mask = zeros(domain_size,domain_size);
   5 %strip = ones(domain_size,1);
   6 mask(:,domain_size/2) = 1;
   7
   8 time_start = -72; time_end = 72;
   9 time_strip = linspace(time_start,time_end,domain_size);
  10 ts2 = linspace(0,1,domain_size); %dummy dimension
  11 [x,y] = meshgrid(time_strip,ts2);
  12 times = zeros(domain_size);
  13 for i = 1:domain_size
  14     times(i,:) = time_strip;
  15 end
  16
  17 dp = detection_probability(times);
  18 d_strip = dp(domain_size/2,:);
  19 figure,plot(time_strip,d_strip),title('detection probability'),
  20 xlabel('time since fire arrival')
  21 ylabel('probability of detection')
  22 ylim([0 1]);
  23
  24 % figure,plot(time_strip,log(d_strip)),title('log detection probability'),
  25 % xlabel('time since fire arrival')
  26 % ylabel('log probability of detection')
  27
  28 % sig = 2.2778;
  29 % gauss_strip = 1/(2*pi*sig^2)*exp(-time_strip.^2/(2*sig^2));
  30 %gauss_strip = exp(-time_strip.^2/(2*sig^2));
  31
  32 % figure,plot(time_strip,gauss_strip),title('geolocation stuff');
  33 % xlabel('time since fire arrival')
  34 % ylabel('gaussian')
  35
  36 % log_g = log(gauss_strip);
  37 % figure,plot(time_strip,log_g),title('log geolocation stuff');
  38 % xlabel('time since fire arrival')
  39 % ylabel('log gaussian')
  40
  41 % new = log(gauss_strip)+log(d_strip);
  42 % figure,plot(time_strip,new),title('new')
  43
  44
  45
  46 like = [];
  47 radius = 60;
  48 weight = gauss_weight(radius);
  49 dx = 30; %in meters, this is grid spacing, dy = dx also
  50
  51 % make distance matrix one time and then pass into the gaussin computation
  52 d_squared = zeros(radius*2+1,radius*2+1);
  53 [m,n] = size(d_squared);
  54 for i = 1:m
  55     for j = i:n
  56         d_squared(i,j) = i^2+j^2+2*(radius+1)^2-2*(radius+1)*(i+j);
  57         if i == j
  58             d_squared(i,j) = 0.5*d_squared(i,j);
  59         end
  60     end
  61 end
  62
  63 d_squared = d_squared+d_squared';
  64 d_squared = dx^2*d_squared;
  65 sig = 1000/3; %1000 meters at sig 3
  66 %figure,mesh(exp(-d_squared/(2*pi*sig^2)));
  67
  68 counter = 1;
  69 for i=radius+1:domain_size-radius
  70     %counter
  71     temp_prob = dp(domain_size/2-radius:domain_size/2+radius,i-radius:i+radius);
  72     %l =  compute_pixel_probability(500,i,heat,radius, weight, detection_probs );
  73     gauss = exp(-d_squared/(2*sig^2));
  74     prob = gauss.*temp_prob;
  75     l = log(sum(prob(:))/(2*pi*sig^2));
  76     like(counter) = l;
  77     counter = counter +1;
  78 end
  79
  80 short_time = time_strip(1,radius+1:end-radius);
  81 l2 = like-max(like(:))+max(like(:))/1000;
  82 figure,plot(short_time,l2),title('pixel log likelihood')
  83 xlabel('Hours since fire arrival')
  84 ylabel('Log likelihood of detection')
  85 %xlim([(round(min(short_time))-1) (round(max(short_time))+1)])
  86 %figure,plot(like)
  87
  88 % compute derivative
  89 l2_prime = zeros(1,length(short_time));
  90 h = short_time(2) - short_time(1);
  91 for i = 2:length(short_time)-1
  92     l2_prime(i) = (l2(i-1)-l2(i+1))/(2*h);
  93 end
  94 l2_prime(1) = l2_prime(2);
  95 l2_prime(end)=l2_prime(end-1);
  96 figure,plot(short_time,l2_prime),title('Derviative')