From 766f98a17aa4ded032bbac8157e91a48977390e1 Mon Sep 17 00:00:00 2001 From: ctwardy Date: Fri, 31 Oct 2014 10:42:40 -0400 Subject: [PATCH 01/50] Fix leaderboard average. (#141) leaderboard() was displaying model_avgrating, but it was stale. Fixed: leaderboard uses its own average *and* updates the database. --- .gitignore | 12 +++++--- Views.py | 85 +++++++++++++++++++++++++++--------------------------- urls.py | 5 +++- 3 files changed, 55 insertions(+), 47 deletions(-) diff --git a/.gitignore b/.gitignore index ab92cef..888cefd 100644 --- a/.gitignore +++ b/.gitignore @@ -1,6 +1,10 @@ *.pyc *.[oa] *~ +*.zip +*.png +.Python +bin/ media/t*png media/profpic*png todo.txt @@ -8,10 +12,10 @@ settings.py temp.txt probmap -*.zip +arc-models/*.pdf -*.png +arc-models/*_flymake.py -arc-models/lognorm_summary.pdf +*.7 -arc-models/motion_model_flymake.py +junk.py \ No newline at end of file diff --git a/Views.py b/Views.py index 046bedd..ce8d101 100644 --- a/Views.py +++ b/Views.py @@ -1,5 +1,5 @@ #!/usr/bin/env python -# -*- coding: UTF-8 -*- +# -*- coding: utf-8 -*- # -*- mode: python -*- # -*- py-indent-offset: 4 -*- @@ -62,8 +62,9 @@ def get_sorted_models(all_models, condition=lambda model: True): ''' Return list of rated models, highest-rated first. Uses model_avgrating attribute and operator.attrgetter method. ''' - rated_models = list(model for model in all_models - if model.model_avgrating != 'unrated' and condition(model)) + rated_models = [model for model in all_models + if model.model_avgrating != 'unrated' + and condition(model)] return sorted(rated_models, key=attrgetter('model_avgrating'), reverse=True) @@ -84,10 +85,10 @@ def confidence_interval(scores): lowerbound = round(avg-halfwidth, 4) upperbound = round(avg+halfwidth, 4) assert lowerbound < avg < upperbound - return (lowerbound, upperbound) + return (lowerbound, upperbound, avg) except: print >> sys.stderr, 'No 95%% CI. N=%d, avg=%6.2f, std=%6.2f' % (N, avg, stdev) - return (0, 0) + return (-1, 1, scores[0]) def check_account_fields(fields, new_user=True): @@ -763,7 +764,7 @@ def model_access(request): request_to_input(request.session, input_dic, 'info', 'error') return render_to_response('ModelScreen.html',input_dic) - # If comming frm account + # If coming frm account else: selection = request.GET['model_in'] if selection == '0': @@ -971,13 +972,19 @@ def completed_test(request): @login_required def leaderboard(request): + '''Generate the data for the main leaderboard: Average, 95% CI, and N. + + Handles filter choice of 'case', 'category', and usual 'model'. + + ''' input_dict = dict(csrf(request)) models = list(Model.objects.filter(Completed_cases__gt=0)) filter_choice = request.POST.get('filter') - instname = lambda model: model.account_set.all()[0].institution_name + if filter_choice == 'case': - model_data, case_name = list(), request.POST.get('case_name') + model_data, case_name = [], request.POST.get('case_name') + input_dict['msg'] = 'of case "%s"' % case_name for model in models: try: test = model.model_tests.get(test_name=case_name, Active=False) @@ -985,28 +992,22 @@ def leaderboard(request): round(float(test.test_rating), 3), '', 1)) except Test.DoesNotExist: pass - input_dict['msg'] = 'of case "%s"' % case_name - elif filter_choice == 'category': - model_data, category = list(), request.POST.get('case_category') - for model in models: - all_tests = model.model_tests.filter(Active=False) - valid_tests = list(test for test in all_tests if test.test_case.subject_category == category) - if len(valid_tests) > 0: - scores = list(float(test.test_rating) for test in valid_tests) - lowerbound, upperbound = (-1, 1) if len(scores) == 1 else confidence_interval(scores) - model_data.append((instname(model), model.model_nameID, - round(sum(scores) / len(scores), 3), '[%5.3f, %5.3f]' % (lowerbound, upperbound), len(scores))) - input_dict['msg'] = 'of category "%s"' % category else: - model_data = list() + if filter_choice == 'category': + category = request.POST.get('case_category') + input_dict['msg'] = 'of "%s"' % category + model_data = [] for model in models: - all_tests = model.model_tests.filter(Active=False) - scores = list(float(test.test_rating) for test in all_tests) - lowerbound, upperbound = (-1, 1) - if len(scores) > 1: - lowerbound, upperbound = confidence_interval(scores) - model_data.append((instname(model), model.model_nameID, - round(float(model.model_avgrating), 3), '[%5.3f, %5.3f]' % (lowerbound, upperbound), len(model.model_tests.all()))) + tests = model.model_tests.filter(Active=False) + if filter_choice == 'category': + tests = [t for t in tests if t.test_case.subject_category == category] + if len(tests) > 0: + scores = [float(test.test_rating) for test in tests] + lowerbound, upperbound, avg = confidence_interval(scores) + model_data.append((instname(model), model.model_nameID, + round(avg, 2), '[%5.2f, %5.2f]' % (lowerbound, upperbound), len(scores))) + model.model_avgrating = '%5.3f' % avg + model.save() input_dict['case_names'] = Case.objects.values_list('case_name', flat=True) input_dict['case_categories'] = set(case.subject_category for case in Case.objects.all()) @@ -1017,7 +1018,11 @@ def leaderboard(request): @login_required def Leader_model(request): - '''Create the leaderboard.''' + '''Duplicate(?) function to create leaderboard. + + TODO: figure out what calls this and why it can't use previous. + + ''' sorted_models = get_sorted_models(Model.objects.all()) # Build Leaderboard inputlist = [] @@ -1027,23 +1032,19 @@ def Leader_model(request): institution = account.institution_name username = account.username name = model.model_nameID - rating = float(model.model_avgrating) + tests = model.model_tests.all() finished_tests = [test for test in tests if not test.Active] - N = len(finished_tests) scores = [float(x.test_rating) for x in finished_tests] + N = len(scores) + lowerbound, upperbound, avg = confidence_interval(scores) # Build case, depending on sample size - case = [institution, name, '%5.3f'%rating, N, username] - if N <= 1: - case.extend([-1, 1, True]) # std=False, sm.sample=True - else: - lowerbound, upperbound = confidence_interval(scores) - if N < 10: # small sample=True - case.extend([lowerbound, upperbound, True]) - else: # small sample = False - case.extend([lowerbound, upperbound, False]) - #print >> sys.stderr, case + case = [institution, name, '%5.3f'%avg, N, username] + if N < 10: # small sample=True + case.extend([lowerbound, upperbound, True]) + else: # small sample = False + case.extend([lowerbound, upperbound, False]) inputlist.append(case) # Prepare variables to send to HTML template @@ -1822,7 +1823,7 @@ def switchboard_toscenario(request): if N < 2: entry.extend([-1, 1, True]) # std=False, sm.sample=True else: - lowerbound, upperbound = confidence_interval(scores) + lowerbound, upperbound, avg = confidence_interval(scores) if N < 10: # small sample=True entry.extend([lowerbound, upperbound, True]) else: # small sample = False diff --git a/urls.py b/urls.py index caf1262..68c07d6 100644 --- a/urls.py +++ b/urls.py @@ -1,4 +1,7 @@ -from django.conf.urls.defaults import patterns, include, url +try: + from django.conf.urls.defaults import patterns, include, url +except: + from django.conf.urls import patterns, include, url from django.conf.urls.static import static from django.http import HttpResponse from django.conf import settings From 71cb7f1bc5350d88a615fb117feadd55f2d7d45d Mon Sep 17 00:00:00 2001 From: ezcawi Date: Thu, 18 Dec 2014 11:23:06 -0500 Subject: [PATCH 02/50] add node adds the best node to the search path. --- senior-design-code/add_node.m | 57 +++++++++++++++++++++++++++++++++++ 1 file changed, 57 insertions(+) create mode 100644 senior-design-code/add_node.m diff --git a/senior-design-code/add_node.m b/senior-design-code/add_node.m new file mode 100644 index 0000000..3188278 --- /dev/null +++ b/senior-design-code/add_node.m @@ -0,0 +1,57 @@ +function [path_sum_new, new_path, probs_update] = add_node(current_path,current_probs, path_sum, m,n,prob_detect) + %this function will insert a point in between the current x,y position + %and the next x,y position in the path or between the previous node and + %current node. Parameters: current_path: the current path of the + %unit, probs: the array of probability distributions + %get neighbors of the different points, path_sum: current path sum, + %m,n: dimensions of array, k: current step on the path + %get the current point on the path and it's immediate neighbors + new_point_choice = zeros(size(current_path)); + best_probs = zeros(length(current_path(1,:)),1); + + for k = 2:(length(current_path(:,1))-1) + current_xy = current_path(k,:); + next_xy = current_path(k+1,:); %next point on path + previous_xy = current_path(k-1,:); %previous point on path + [ncx,ncy] = get_neighbors(current_xy(1), current_xy(2), m,n); + neighbors = [ncx ncy]; + [nnx,nny] = get_neighbors(next_xy(1), next_xy(2),m,n); + next_neighbors = [nnx nny]; + [npx,npy] = get_neighbors(previous_xy(1), previous_xy(2),m,n); + prev_neighbors = [npx npy]; + v1 = intersect(neighbors,next_neighbors,'rows'); + v2 = intersect(neighbors,prev_neighbors,'rows'); + valid_choices = union(v1,v2,'rows'); %the choices for the new point based on the old points + %get intersection of all three arrays as valid points + %look at the probabilities of the + n_probs = zeros(length(valid_choices(:,1)),1); + for p = 1: length(valid_choices(:,1)) + n_probs(p) = current_probs(valid_choices(p,2), valid_choices(p,1)); + end + a = find(n_probs == max(n_probs)); + new_point_choice(k,:) = valid_choices(a,:); + best_probs(k) = current_probs(new_point_choice(k,2), new_point_choice(k,1)); + end + best_point = find(best_probs == max(best_probs)); + if length(best_point)>1 + best_point = best_point(1);%pick the first best choice if there are more than one + end + %update path sum new + current_probs(new_point_choice(best_point,2), new_point_choice(best_point,1)); + path_sum_new = path_sum + prob_detect*current_probs(new_point_choice(best_point,2), new_point_choice(best_point,1)); + probs_update = current_probs; + probs_update(new_point_choice(best_point,2), new_point_choice(best_point,1)) = (1-prob_detect)*probs_update(new_point_choice(best_point,2), new_point_choice(best_point,1));%reflect that the new point has been visited in the probabilities matrix + %update new path + current_xy = current_path(best_point,:); %the current point that yielded the new best point + [ncx,ncy] = get_neighbors(current_xy(1), current_xy(2), m,n); + neighbors = [ncx ncy]; + previous_xy = current_path(best_point-1,:); %previous point on path + [npx,npy] = get_neighbors(previous_xy(1), previous_xy(2),m,n); + prev_neighbors = [npx npy]; + v = intersect(neighbors,prev_neighbors,'rows'); + if ~isempty(intersect(v,new_point_choice(best_point,:), 'rows')) %check to see if it's between the kth point and the previous neighbor + new_path = [current_path(1:best_point-1,:) ; new_point_choice(best_point,:); current_path(best_point:end,:)]; + else %if not,it's between the current point and the next point in the path + new_path = [current_path(1:best_point,:); new_point_choice(best_point,:); current_path(best_point+1:end,:)]; + end +end From eb1c32ce0d3f197c80b17de1c4d88f129d22eb18 Mon Sep 17 00:00:00 2001 From: ezcawi Date: Thu, 18 Dec 2014 11:24:49 -0500 Subject: [PATCH 03/50] finds the best node to replace in the path --- senior-design-code/replace_node.m | 40 +++++++++++++++++++++++++++++++ 1 file changed, 40 insertions(+) create mode 100644 senior-design-code/replace_node.m diff --git a/senior-design-code/replace_node.m b/senior-design-code/replace_node.m new file mode 100644 index 0000000..edcff43 --- /dev/null +++ b/senior-design-code/replace_node.m @@ -0,0 +1,40 @@ +function [path_sum_new, new_path, probs_update] = replace_node(current_path, probs, current_probs,m,n,prob_detect) + %function to replace the kth element of the path. parameters: + %current_path: array representing column, row elements of the current + %path. probs: array of probabilities at the area needed to recompute path sum, current_probs m,n: dimensions of + %array, k, current spot on the path + new_point_choice = zeros(size(current_path)); + best_probs = zeros(length(current_path(1,:)),1); + for k = 2:(length(current_path(:,1))-1) + current_xy = current_path(k,:); + next_xy = current_path(k+1,:); %previous point on path + previous_xy = current_path(k-1,:); %next point on path + [nnx,nny] = get_neighbors(next_xy(1), next_xy(2),m,n); %don't need neighbors of current point because that's being gotten rid of + next_neighbors = [nnx nny]; + [npx,npy] = get_neighbors(previous_xy(1), previous_xy(2),m,n); + prev_neighbors = [npx npy] ; + valid_choices = intersect(next_neighbors, prev_neighbors,'rows');%the idea is that the new in between point should be a neighbor of both points + %choose valid choices as in "add_node" + n_probs = zeros(length(valid_choices(:,1)),1); + for i = 1: length(valid_choices(:,1)) + n_probs(i) = current_probs(valid_choices(i,2), valid_choices(i,1));%have to look at the current probs + end + a = find(n_probs == max(n_probs)); + new_point_choice(k,:) = valid_choices(a,:); + best_probs(k) = probs(new_point_choice(k,2), new_point_choice(k,1)); + end + best_point = find(best_probs == max(best_probs)); + if length(best_point)>1 + best_point = best_point(1);%pick the first best choice if there are more than one + end + %update probability + new_path = [current_path(1:best_point-1,:); new_point_choice(best_point,:); current_path(best_point+1:end,:)]; + %have to update the path sum the same way as in the first loop because + %if a point is replaced it doesn't factor into the whole sum + path_sum_new = 0; + probs_update = probs; + for j = 1:length(current_path(:,1)) + path_sum_new = path_sum_new + prob_detect*probs(new_path(j,2), new_path(j,1)); + probs_update(new_path(j,2), new_path(j,1)) = probs(new_path(j,2), new_path(j,1))*(1-prob_detect); + end +end From ff6a6f3fe7393235c6a53c3de7cace961b6c1f88 Mon Sep 17 00:00:00 2001 From: ezcawi Date: Thu, 18 Dec 2014 11:25:28 -0500 Subject: [PATCH 04/50] gets neighbors in array indeces, checks for edges --- senior-design-code/get_neighbors | 36 ++++++++++++++++++++++++++++++++ 1 file changed, 36 insertions(+) create mode 100644 senior-design-code/get_neighbors diff --git a/senior-design-code/get_neighbors b/senior-design-code/get_neighbors new file mode 100644 index 0000000..563fed5 --- /dev/null +++ b/senior-design-code/get_neighbors @@ -0,0 +1,36 @@ +function [neighbors_x, neighbors_y] = get_neighbors(x,y,m,n) + % x and y are assumed to be column and row indeces, respectively, of + % an mxn matrix. gets all neighbors, including diagonals + + %check for corner cases + if (x == 1) && (y==1) %top left corner + neighbors_x = [1;2;2]; + neighbors_y = [2; 2;1]; + elseif (x == m) && (y==1) %top right corner + neighbors_x = [n; n-1; n-1]; + neighbors_y = [2; 2; 1]; + elseif (x==1) && (y==n) %bottom left corner + neighbors_x = [1;2;2]; + neighbors_y = [n-1; n-1; n]; + elseif (x==m) && (y ==n) %bottom right corner + neighbors_x = [n, n-1, n-1]'; + neighbors_y = [m-1, m-1, m]'; + %check for edge cases + elseif (x~=1) && (y==1) %top edge + neighbors_x = [x-1, x-1, x, x+1, x+1]'; + neighbors_y = [1,2, 2, 2, 1]'; + elseif (x ~=1) && (y == n) %bottom edge + neighbors_x = [x-1, x-1, x, x+1, x+1]'; + neighbors_y = [n, n-1, n-1, n-1, n]'; + elseif (x==1) && (y~=1) %left edge + neighbors_x = [ 1, 2, 2, 2, 1]'; + neighbors_y = [y-1, y-1, y, y+1, y+1]'; + elseif (x == m) && (y~=1) %right edge + neighbors_x = [m, m-1, m-1, m-1, m]'; + neighbors_y = [y-1, y-1, y, y+1, y+1]'; + %middle case + else + neighbors_x = [x-1, x-1, x-1, x, x, x+1, x+1, x+1]'; + neighbors_y = [y-1, y, y+1, y-1, y+1, y-1, y, y+1]'; + end +end From 2a2b6968fd20ccf359a11ee1a25f6af40a37c679 Mon Sep 17 00:00:00 2001 From: ezcawi Date: Thu, 18 Dec 2014 11:25:57 -0500 Subject: [PATCH 05/50] greedy algorithm, no fixed start or end point --- senior-design-code/greedy_algorithm.m | 74 +++++++++++++++++++++++++++ 1 file changed, 74 insertions(+) create mode 100644 senior-design-code/greedy_algorithm.m diff --git a/senior-design-code/greedy_algorithm.m b/senior-design-code/greedy_algorithm.m new file mode 100644 index 0000000..458b63b --- /dev/null +++ b/senior-design-code/greedy_algorithm.m @@ -0,0 +1,74 @@ +%greedy path finding algorithm for search and rescue +%first, load in the array of probabilities, for this example i just picked +%a matrix from a uniform distribution between 0 and 1 and then turn that +%into probabilities that sum to 1 +function best_sum = greedy_algorithm(m,n,probs,prob_detect) +%starting off with 4x4 matrix for simplicity + %turn the example set into probabilities. +%set a cost for the search, here it is arbitrarily 75% of cells that the quad can cover, +%can determine this through flight time, etc +tic +num_cells = ceil(m*n*0.75); + +path_sum = zeros(m,n); +current_x = zeros(m,n,num_cells); +current_y = zeros(m,n,num_cells); %this will keep track of the path for each starting point +%to find the best path moving through only some of the cells,we want to know the best starting point +for x = 1:1:n + for y = 1:1:m %first two loops desginate going through the different starting points + search_value_temp = probs; %going to be changing this matrix so i set the thing + current_x(y,x,1) = x; + current_y(y,x,1) = y; + path_sum(y,x) =prob_detect*search_value_temp(y, x);%initialize path probability sum, using probability of detection because we might miss + search_value_temp(current_y(y,x,1), current_x(y,x,1)) = search_value_temp(current_y(y,x,1), current_x(y,x,1))*(1 - prob_detect); %subtracting some probability from the first cell + + for k = 2:1:num_cells + + [neighbors_x, neighbors_y] = get_neighbors(current_x(y,x,k-1), current_y(y,x,k-1), n, n); + current_neighbors = [neighbors_x neighbors_y]; + + neighbor_prob = zeros(size(neighbors_x)); + for i = 1:length(neighbors_x) + neighbor_prob(i) = search_value_temp(current_neighbors(i,2),current_neighbors(i,1)); + end + best_neighbor = find(neighbor_prob == max(neighbor_prob)); + next_x = current_neighbors(best_neighbor,1); + next_y = current_neighbors(best_neighbor,2); + if length(next_x)>1 %if there are two equal probability cells in the neighbors list just pick one randomly + p = ranperm(length(next_x), 1); + current_x(y,x,k) = next_x(p); + current_y(y,x,k) = next_y(p); + else + current_x(y,x,k) = next_x; + current_y(y,x,k) = next_y; + end; + path_sum(y,x) = path_sum(y,x) + prob_detect*search_value_temp(current_y(y,x,k), current_x(y,x,k)); %update the path sum + search_value_temp(current_y(y,x,k), current_x(y,x,k)) = search_value_temp(current_y(y,x,k), current_x(y,x,k))*(1 - prob_detect); + end + end +end +[best_start_y, best_start_x] = find(path_sum == max(max(path_sum))); +if length(best_start_x)>1 + p = randperm(length(best_start_x),1); + best_start_x = best_start_x(p); + best_start_y = best_start_y(p); +end +best_sum = path_sum(best_start_y, best_start_x); +best_pathx = zeros(num_cells,1); +best_pathy = zeros(num_cells,1); +best_pathtempx = current_x(best_start_y, best_start_x, :); +best_pathtempy = current_y(best_start_y, best_start_x, :); +toc +for i =1:num_cells + best_pathx(i) = best_pathtempx(1,1,i); + best_pathy(i) = best_pathtempy(1,1,i); +end +best_path = [best_pathx best_pathy]; +% display(sprintf('the best greedy path sum is: %d',best_sum)) +% display(sprintf('the best greedy starting point is %d , %d', best_start_x, best_start_y)) +% display(sprintf('the best greedy path found is: ')); +% display(best_path) +toc + +%needed, vectorize for loops to make the program scaleable to a larger array, +%firgure out higher level algorithm to get the effort put into cells From 046b2af53f9c6a77ba5a1b924150996a759a2ec0 Mon Sep 17 00:00:00 2001 From: ezcawi Date: Thu, 18 Dec 2014 11:26:44 -0500 Subject: [PATCH 06/50] greedy hill climbing with random restart --- senior-design-code/hill_climbing.m | 48 ++++++++++++++++++++++++++++++ 1 file changed, 48 insertions(+) create mode 100644 senior-design-code/hill_climbing.m diff --git a/senior-design-code/hill_climbing.m b/senior-design-code/hill_climbing.m new file mode 100644 index 0000000..9b8c9c0 --- /dev/null +++ b/senior-design-code/hill_climbing.m @@ -0,0 +1,48 @@ +function best_path_sum = hill_climbing(m,n,probs,prob_detect,C,spt, ept) + tic + i_best_path = zeros(C,2,10); %for storing the best path + i_path_sum = zeros(1,10); + path_sum_start = zeros(1,10); + for i = 1:10 + probs_current = probs; + %first find random starting path with length less than or equal to C + l = C+1; + while l>C + current_path = gen_rand_path(spt,ept,m,n); + l = length(current_path(:,1)); + end + %compute first path sum + path_sum = 0; + for j = 1:length(current_path(:,1)) + path_sum = path_sum + prob_detect*probs_current(current_path(j,2), current_path(j,1)); + probs_current(current_path(j,2), current_path(j,1)) = probs_current(current_path(j,2), current_path(j,1))*(1-prob_detect); + end + path_sum_start(i) = path_sum; + %hill climbing loop: two rules: add a node, replace a node + for a = 1:100 + if length(current_path(:,1)) 0 %new path is a better sum + path_sum = path_sum_new; %take new path/new path sum + current_path = new_path; + probs_current = probs_update; + end %if the old path was better, don't update anything + end + i_path_sum(i) = path_sum; %store the current path sum + i_best_path(:,:,i) = current_path; %store the current path for access later + end + + best_path_sum = max(i_path_sum); + best_i = find(i_path_sum == best_path_sum); + best_path = i_best_path(:,:,best_i); + best_path_start = path_sum_start(best_i); +% disp('best hill climbing path is:') +% disp(best_path); +% disp('best hill climbing sum is:') +% disp(best_path_sum); + + toc From 5185d56799cd3ae92fdea20d667949f9fa7e2d6b Mon Sep 17 00:00:00 2001 From: ezcawi Date: Thu, 18 Dec 2014 11:28:13 -0500 Subject: [PATCH 07/50] compares greedy to hill climbing for values from uniform distribution --- senior-design-code/comparison.m | 20 ++++++++++++++++++++ 1 file changed, 20 insertions(+) create mode 100644 senior-design-code/comparison.m diff --git a/senior-design-code/comparison.m b/senior-design-code/comparison.m new file mode 100644 index 0000000..574b552 --- /dev/null +++ b/senior-design-code/comparison.m @@ -0,0 +1,20 @@ +%hill climbing +clear all; close all; + +m = 5; +n = 5; +prob_detect = .9; +spt = [1 1]; +ept = [4 5]; + +C = floor(m*n*.75); %total number of cells drone can travel to +for i = 1:30 + search_value = rand(m,n); + probs = search_value/sum(sum(search_value)); %turn the example set into probabilities. + best_hill_sum(i) = hill_climbing(m,n,probs,prob_detect,C,spt, ept); + %compare to greedy algorithm for the same set of probabilities + best_greedy_sum(i) = greedy_algorithm(m,n,probs,prob_detect); + i +end +hill_avg = mean(best_hill_sum); +greedy_avg = mean(best_greedy_sum); From 1330a3968758b8bd3fa1223f73214b6a4f9a52e2 Mon Sep 17 00:00:00 2001 From: ezcawi Date: Thu, 18 Dec 2014 11:29:20 -0500 Subject: [PATCH 08/50] current iteration of big path algorithm needs: to be finished, currently figuring out the good way to get the next search region --- senior-design-code/big_path | 65 +++++++++++++++++++++++++++++++++++++ 1 file changed, 65 insertions(+) create mode 100644 senior-design-code/big_path diff --git a/senior-design-code/big_path b/senior-design-code/big_path new file mode 100644 index 0000000..d25d188 --- /dev/null +++ b/senior-design-code/big_path @@ -0,0 +1,65 @@ +function search_path = big_path(N,probs,C_total,percentages,spt) +%this function will take a large probability distribution and number of +%search areas to create a group of squares - want to extend this to human +%created search regions at some point, possibly with mapsar/translate into +%an array? +%parameters: N should be even to facilitate region creation in the first +%iteration. probs: the array of search probabilities. C: the total cost the +%unit is allowed to travel. percentages: how much should be spent in each +%region, should be length-N vector adding up to 1. Region 1 will be the top left, region +%two will be on it's right, and so on until it reaches the right end of the +%array and then the counting will start again from the left. Spt: user +%defined starting point, length 2 vector of pixel indeces + +m = floor(N/2); +region_width = length(probs(1,:))/m; %for this to work need probs dimensions to be multiple of N +region_height = length(probs(:,1))/m; +num_across = length(probs(1,:))/region_width +num_down = length(probs(:,1))/region_height +current_region = probs(1:region_height,1:region_width); %initialize current search region +current_path = spt; +current_across = 1; +current_down = 1; %location of the current search region +next_across = current_across; +next_down = current_down; +next_width_start = 1; +next_height_start = 1; +for i = 1:N + + %create next search region + + if direction_flag == 0 %going from left to right + if current_across1 + next_across = current_across - 1 + direction_flag = 1; + next_width_start = next_across*region_width + 1; + elseif current_across==1 + next_down = next_down + 1 + direction_flag = 0; %start back the other way + next_height_start = + end + next_region = probs(next_height_start: (next_height_start+region_height -1), next_width_start:(next_width_start + region_width -1)); + + %examine interface between search regions to determine the ending point + %use percentages to determine the cost allowed in the current area + + %run hill climbing algorithm + + %append generated path to current path + + %update starting point and region of interest for the next iteration, unless i = N, then we + %are done + +end + +%compute ending path sum From 8b4aa7c8abe1608cf8f8b324ec6ec044a4ba1d18 Mon Sep 17 00:00:00 2001 From: ezcawi Date: Thu, 8 Jan 2015 11:25:01 -0500 Subject: [PATCH 09/50] Create hand_test_cases.m --- senior-design-code/hand_test_cases.m | 43 ++++++++++++++++++++++++++++ 1 file changed, 43 insertions(+) create mode 100644 senior-design-code/hand_test_cases.m diff --git a/senior-design-code/hand_test_cases.m b/senior-design-code/hand_test_cases.m new file mode 100644 index 0000000..e203866 --- /dev/null +++ b/senior-design-code/hand_test_cases.m @@ -0,0 +1,43 @@ +%test cases for small matrices on greedy and hill climbing, which can be +%verified by hand to ensure algorithms are working +clear all; close all; + +m = 5; +n = 5; +prob_detect = .9; +spt = [1 1]; +ept = [5 5]; + +C = floor(m*n*.75); %total number of cells drone can travel to +best_hill_path = zeros(C,2,3); +best_hill_sum = zeros(1,3); +best_greedy_path = zeros(C,2,3); +best_greed_sum = zeros(1,3); +%%magic square matrix +search_value = magic(m); +probs = search_value/sum(sum(search_value)); %turn the example set into probabilities. +[best_hill_path(:,:,1),best_hill_sum(1)] = hill_climbing(m,n,probs,prob_detect,C,spt, ept); +%compare to greedy algorithm for the same set of probabilities +[best_greedy_path(:,:,1),best_greedy_sum(1)] = greedy_algorithm(m,n,probs,prob_detect); + +search_value = gallery('cauchy',1:m,1:n); +probs = search_value/sum(sum(search_value)); %turn the example set into probabilities. +[best_hill_path(:,:,2),best_hill_sum(2)] = hill_climbing(m,n,probs,prob_detect,C,spt, ept); +%compare to greedy algorithm for the same set of probabilities +[best_greedy_path(:,:,2),best_greedy_sum(2)] = greedy_algorithm(m,n,probs,prob_detect); + +%circular bivariate normal distribution, ro = 0, mu = 0, sigma = 1 for x +%and y +y = linspace(-1,1,m); +x = linspace(-1,1,n); +[xx,yy] = meshgrid(x,y); +mu = [0 0]; +sig = [1 0; 0 1]; +search_value = mvnpdf([xx(:),yy(:)],mu, sig); +search_value = reshape(search_value, length(y), length(x)); +probs = search_value/sum(sum(search_value)); %turn the example set into probabilities. +[best_hill_path(:,:,3),best_hill_sum(3)] = hill_climbing(m,n,probs,prob_detect,C,spt, ept); +%compare to greedy algorithm for the same set of probabilities +[best_greedy_path(:,:,3),best_greedy_sum(3)] = greedy_algorithm(m,n,probs,prob_detect); +hill_avg = mean(best_hill_sum); +greedy_avg = mean(best_greedy_sum); From 1321c0fbd476fba0921893514a5355ef435b14af Mon Sep 17 00:00:00 2001 From: ezcawi Date: Thu, 8 Jan 2015 11:39:15 -0500 Subject: [PATCH 10/50] Create LatLon_to_m.py uses pixel resolution and size of the array to get the meters x and y coordinates of an array. First assuming 0 degree heading, then rotated --- senior-design-code/LatLon_to_m.py | 42 +++++++++++++++++++++++++++++++ 1 file changed, 42 insertions(+) create mode 100644 senior-design-code/LatLon_to_m.py diff --git a/senior-design-code/LatLon_to_m.py b/senior-design-code/LatLon_to_m.py new file mode 100644 index 0000000..34d7082 --- /dev/null +++ b/senior-design-code/LatLon_to_m.py @@ -0,0 +1,42 @@ +import math +import numpy as np +def array_to_meters(probs, center_lat, center_lon, heading,pixel_res) +#probs - numpy array of probabilities, will need size of the array +#center_lat, center_lon - central coordinates of the array +#heading - angle of orientation of the array, measured from north/y axis +#pixel_res - number of meters per pixel + shp = probs.shape + xs = np.zeros(shp) + ys = np.zeros(shp) + #convert central point to meters + xy_c = coord_to_meters(center_lat,center_lon) + c_index = [shp[0]/2 - 1,shp[1]/2 -1] + for i in range(shp[0]) + dy_index = i - c_index[0] + dy = dy_index*pixel_res + ys[i][:] = xy_c[1] + dy + for j in range(shp[1]) + dx_index = j - c_index[1] + dx = dx_index*pixel_res + xs[:][j] = xy_c[0] + dx + #now need to rotate each x,y by the heading + rot_matrix = np.array([[math.cos(-1*heading), -1*math.sin(-1*heading)],[math.sin(-1*heading), math.cos(-1*heading)]]) + xy_rot = np.zeros((shp[0],shp[1],2)) + for i in range(shp[0]) + for j in range(shp[1]) + xy_vec = np.array([[xs[i][j]],[ys[i][j]]) + xy_rot[i][j][:] = np.dot(rot_matrix,xy_vec) + return xy_rot #returns array of x,y coordinates for each element of the array +def coord_to_meters(lat,lon) + lat_degrees = float(lat) + long_degrees = float(lng) + + lat_radians = float(lat_degrees)*(math.pi/180) + long_radians = float(long_degrees)*(math.pi/180) + + R = 6378137 + XY = [0,0] + XY[0] = R*math.cos(lat_radians)*math.cos(long_radians) + XY[1] = R*math.cos(lat_radians)*math.sin(long_radians) + Z = R*math.sin(lat_radians) #not sure if this is right/needed if pixel resolution is known + return XY From 2884d061f50ce78bf229065736e8230829ba673e Mon Sep 17 00:00:00 2001 From: ezcawi Date: Thu, 8 Jan 2015 11:39:39 -0500 Subject: [PATCH 11/50] Update add_node.m --- senior-design-code/add_node.m | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/senior-design-code/add_node.m b/senior-design-code/add_node.m index 3188278..fcc30bb 100644 --- a/senior-design-code/add_node.m +++ b/senior-design-code/add_node.m @@ -29,7 +29,7 @@ n_probs(p) = current_probs(valid_choices(p,2), valid_choices(p,1)); end a = find(n_probs == max(n_probs)); - new_point_choice(k,:) = valid_choices(a,:); + new_point_choice(k,:) = valid_choices(a(1),:); best_probs(k) = current_probs(new_point_choice(k,2), new_point_choice(k,1)); end best_point = find(best_probs == max(best_probs)); @@ -37,7 +37,6 @@ best_point = best_point(1);%pick the first best choice if there are more than one end %update path sum new - current_probs(new_point_choice(best_point,2), new_point_choice(best_point,1)); path_sum_new = path_sum + prob_detect*current_probs(new_point_choice(best_point,2), new_point_choice(best_point,1)); probs_update = current_probs; probs_update(new_point_choice(best_point,2), new_point_choice(best_point,1)) = (1-prob_detect)*probs_update(new_point_choice(best_point,2), new_point_choice(best_point,1));%reflect that the new point has been visited in the probabilities matrix From 2d2db006ef3f90f6284f8544ae2cc4b7dc218f68 Mon Sep 17 00:00:00 2001 From: ezcawi Date: Thu, 8 Jan 2015 11:40:18 -0500 Subject: [PATCH 12/50] Update greedy_algorithm.m --- senior-design-code/greedy_algorithm.m | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/senior-design-code/greedy_algorithm.m b/senior-design-code/greedy_algorithm.m index 458b63b..e601585 100644 --- a/senior-design-code/greedy_algorithm.m +++ b/senior-design-code/greedy_algorithm.m @@ -2,13 +2,13 @@ %first, load in the array of probabilities, for this example i just picked %a matrix from a uniform distribution between 0 and 1 and then turn that %into probabilities that sum to 1 -function best_sum = greedy_algorithm(m,n,probs,prob_detect) +function [best_path,best_sum] = greedy_algorithm(m,n,probs,prob_detect) %starting off with 4x4 matrix for simplicity %turn the example set into probabilities. %set a cost for the search, here it is arbitrarily 75% of cells that the quad can cover, %can determine this through flight time, etc tic -num_cells = ceil(m*n*0.75); +num_cells = floor(m*n*0.75); path_sum = zeros(m,n); current_x = zeros(m,n,num_cells); @@ -35,7 +35,7 @@ next_x = current_neighbors(best_neighbor,1); next_y = current_neighbors(best_neighbor,2); if length(next_x)>1 %if there are two equal probability cells in the neighbors list just pick one randomly - p = ranperm(length(next_x), 1); + p = randperm(length(next_x), 1); current_x(y,x,k) = next_x(p); current_y(y,x,k) = next_y(p); else @@ -58,12 +58,13 @@ best_pathy = zeros(num_cells,1); best_pathtempx = current_x(best_start_y, best_start_x, :); best_pathtempy = current_y(best_start_y, best_start_x, :); -toc + for i =1:num_cells best_pathx(i) = best_pathtempx(1,1,i); best_pathy(i) = best_pathtempy(1,1,i); end best_path = [best_pathx best_pathy]; +size(best_path) % display(sprintf('the best greedy path sum is: %d',best_sum)) % display(sprintf('the best greedy starting point is %d , %d', best_start_x, best_start_y)) % display(sprintf('the best greedy path found is: ')); From 56ee1c19ff770c79a251f56cf1503dc08ffeae89 Mon Sep 17 00:00:00 2001 From: ezcawi Date: Thu, 8 Jan 2015 11:40:46 -0500 Subject: [PATCH 13/50] Update hill_climbing.m --- senior-design-code/hill_climbing.m | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/senior-design-code/hill_climbing.m b/senior-design-code/hill_climbing.m index 9b8c9c0..fc16e3f 100644 --- a/senior-design-code/hill_climbing.m +++ b/senior-design-code/hill_climbing.m @@ -1,4 +1,9 @@ -function best_path_sum = hill_climbing(m,n,probs,prob_detect,C,spt, ept) +function [best_path,best_path_sum] = hill_climbing(m,n,probs,prob_detect,C,spt, ept) + %m,n - dimensions of matrix + %probs - array of probability that the subject is in location + %prob_detect - probability that searcher will see subject + %C - total cost i.e. number of cells that can be visited + %spt, ept - starting and ending points, respectively tic i_best_path = zeros(C,2,10); %for storing the best path i_path_sum = zeros(1,10); From 1c055f81c33b3fd438201e15ca5cbd0e35e41898 Mon Sep 17 00:00:00 2001 From: ezcawi Date: Thu, 8 Jan 2015 11:41:41 -0500 Subject: [PATCH 14/50] Update replace_node.m --- senior-design-code/replace_node.m | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/senior-design-code/replace_node.m b/senior-design-code/replace_node.m index edcff43..34fbedf 100644 --- a/senior-design-code/replace_node.m +++ b/senior-design-code/replace_node.m @@ -20,7 +20,7 @@ n_probs(i) = current_probs(valid_choices(i,2), valid_choices(i,1));%have to look at the current probs end a = find(n_probs == max(n_probs)); - new_point_choice(k,:) = valid_choices(a,:); + new_point_choice(k,:) = valid_choices(a(1),:); best_probs(k) = probs(new_point_choice(k,2), new_point_choice(k,1)); end best_point = find(best_probs == max(best_probs)); From dcdcafe5b935444c6b39864d23d0fd1eee122598 Mon Sep 17 00:00:00 2001 From: ezcawi Date: Thu, 8 Jan 2015 12:00:36 -0500 Subject: [PATCH 15/50] Create gen_rand_path.m --- senior-design-code/gen_rand_path.m | 16 ++++++++++++++++ 1 file changed, 16 insertions(+) create mode 100644 senior-design-code/gen_rand_path.m diff --git a/senior-design-code/gen_rand_path.m b/senior-design-code/gen_rand_path.m new file mode 100644 index 0000000..f9080e3 --- /dev/null +++ b/senior-design-code/gen_rand_path.m @@ -0,0 +1,16 @@ +function rand_path = gen_rand_path(spt, ept,m,n) + %function to generate random starting path + rand_path = spt; + x = spt(1); + y = spt(2); + current_end = spt; + while isempty(intersect(ept,current_end,'rows')) + [x_neighbors,y_neighbors] = get_neighbors(x, y, m,n); + k = randi([1 length(x_neighbors)],1); + x = x_neighbors(k); %update x and y + y = y_neighbors(k); + current_end = [x y]; + rand_path = [rand_path; current_end]; + + end +end From 49a39090e6c7f0b8d43c8afddc2a73348a7cddec Mon Sep 17 00:00:00 2001 From: ezcawi Date: Thu, 8 Jan 2015 14:41:27 -0500 Subject: [PATCH 16/50] Update big_path --- senior-design-code/big_path | 98 +++++++++++++++---------------------- 1 file changed, 39 insertions(+), 59 deletions(-) diff --git a/senior-design-code/big_path b/senior-design-code/big_path index d25d188..bb509fa 100644 --- a/senior-design-code/big_path +++ b/senior-design-code/big_path @@ -1,65 +1,45 @@ -function search_path = big_path(N,probs,C_total,percentages,spt) -%this function will take a large probability distribution and number of -%search areas to create a group of squares - want to extend this to human -%created search regions at some point, possibly with mapsar/translate into -%an array? -%parameters: N should be even to facilitate region creation in the first -%iteration. probs: the array of search probabilities. C: the total cost the -%unit is allowed to travel. percentages: how much should be spent in each -%region, should be length-N vector adding up to 1. Region 1 will be the top left, region -%two will be on it's right, and so on until it reaches the right end of the -%array and then the counting will start again from the left. Spt: user -%defined starting point, length 2 vector of pixel indeces +%stitching paths together +%human planning level +%assumptions: rectangular search region, split up into other rectangular +%regions, search regions numbered columnwise (standard matlab convention)as [1 4 7; 2 5 8; 3 6 9], each region is +%the same size, want to spend time in every search region +%inputs from user: probability files, number of levels to go down - N +% the inputs should be interfaced with the gui somehow -m = floor(N/2); -region_width = length(probs(1,:))/m; %for this to work need probs dimensions to be multiple of N -region_height = length(probs(:,1))/m; -num_across = length(probs(1,:))/region_width -num_down = length(probs(:,1))/region_height -current_region = probs(1:region_height,1:region_width); %initialize current search region -current_path = spt; -current_across = 1; -current_down = 1; %location of the current search region -next_across = current_across; -next_down = current_down; -next_width_start = 1; -next_height_start = 1; -for i = 1:N - - %create next search region - - if direction_flag == 0 %going from left to right - if current_across2 + dx_old = current_path(j,1) - current_path(j-1,1); %previous direction, determines starting point + dy_old = current_path(j,2) - current_path(j-1,2); + else + dx_old = 0; + dy_old = 0; end - else %going from left to right - if current_across>1 - next_across = current_across - 1 - direction_flag = 1; - next_width_start = next_across*region_width + 1; - elseif current_across==1 - next_down = next_down + 1 - direction_flag = 0; %start back the other way - next_height_start = + [spt,ept] = endpts(dx_old, dy_old, dx_new, dy_new,m_new, n_new); + [path_part(:,:,j), path_sum] = hill_climbing(m_new,n_new,current_res_probs,prob_detect,C_area,spt,ept); end - next_region = probs(next_height_start: (next_height_start+region_height -1), next_width_start:(next_width_start + region_width -1)); - - %examine interface between search regions to determine the ending point - %use percentages to determine the cost allowed in the current area - - %run hill climbing algorithm - - %append generated path to current path - - %update starting point and region of interest for the next iteration, unless i = N, then we - %are done + %link the paths together + [current_path,higher_res_probs] = link_paths(path_part,current_path,m_old, n_old, m_new,n_new,current_res_probs, higher_res_probs); end -%compute ending path sum +%compute the final path sum From 5a5744944227e22c40aaf28efd21d204f9025270 Mon Sep 17 00:00:00 2001 From: ezcawi Date: Thu, 8 Jan 2015 14:48:07 -0500 Subject: [PATCH 17/50] Update big_path --- senior-design-code/big_path | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/senior-design-code/big_path b/senior-design-code/big_path index bb509fa..9302daf 100644 --- a/senior-design-code/big_path +++ b/senior-design-code/big_path @@ -13,7 +13,7 @@ C = input('What percentage of the total area can the UAV cover?'); % needs a num starting_path = input('What is the starting path through the top level?') % needs an n x 2 array %need to check inputs for final version current_path = starting_path; -higher_res_probs = top_level_probs; +lower_res_probs = top_level_probs; [m_old, n_old] = size(top_level_probs); for i = 1:N-1 %for every level of resolution for j = 1: length(current_path(:,1))-1 @@ -21,7 +21,7 @@ for i = 1:N-1 %for every level of resolution %resolutions probabilities and segment but segmentation indexing %is proving difficult to wrap my head around. current_res_probs(:,:,j) = load(input('next probability array','s')); - C_area = floor(C*higher_res_probs(current_path(j,2), current_path(j,1))*numel(current_res_probs)); %current cost * probability subject is in an area = proportional time spent in that area * number of cells = number of cells in new resolution we can cover + C_area = floor(C*lower_res_probs(current_path(j,2), current_path(j,1))*numel(current_res_probs)); %current cost * probability subject is in an area = proportional time spent in that area * number of cells = number of cells in new resolution we can cover [m_new,n_new] = size(current_res_probs); %calculate starting and ending points for this level from previous %path resolution @@ -39,7 +39,7 @@ for i = 1:N-1 %for every level of resolution end %link the paths together - [current_path,higher_res_probs] = link_paths(path_part,current_path,m_old, n_old, m_new,n_new,current_res_probs, higher_res_probs); + [current_path,res_probs_update] = link_paths(path_part,current_path,m_old, n_old, m_new,n_new,current_res_probs, lower_res_probs); end %compute the final path sum From 39189208f444f66b8114c7c8034221e73624c90c Mon Sep 17 00:00:00 2001 From: ezcawi Date: Thu, 8 Jan 2015 14:48:30 -0500 Subject: [PATCH 18/50] Create link_paths.m --- senior-design-code/link_paths.m | 9 +++++++++ 1 file changed, 9 insertions(+) create mode 100644 senior-design-code/link_paths.m diff --git a/senior-design-code/link_paths.m b/senior-design-code/link_paths.m new file mode 100644 index 0000000..6c71e14 --- /dev/null +++ b/senior-design-code/link_paths.m @@ -0,0 +1,9 @@ +function [current_path,higher_res_probs] = link_paths(path_part,current_path,m_old, n_old, m_new,n_new,current_res_probs, lower_res_probs) + %path_part - array of all the different path parts + %m_old, n_old - dimensions of higher resolution array + %m_new, n_new - dimensions of newer resolution arrays + %current_res_probs - arrays of probabilities in current resolution + %higher_res_probs - arrays of probabilities in lower resolution + + +end From f7842f8ae4b95ba06e1a5bc081acfb3dc18ab120 Mon Sep 17 00:00:00 2001 From: ezcawi Date: Thu, 8 Jan 2015 14:48:49 -0500 Subject: [PATCH 19/50] Create endpts.m --- senior-design-code/endpts.m | 11 +++++++++++ 1 file changed, 11 insertions(+) create mode 100644 senior-design-code/endpts.m diff --git a/senior-design-code/endpts.m b/senior-design-code/endpts.m new file mode 100644 index 0000000..4c2c3b1 --- /dev/null +++ b/senior-design-code/endpts.m @@ -0,0 +1,11 @@ +function [spt,ept] = endpts(dx_old, dy_old, dx_new, dy_new,m, n) + %dx_old, dy_old - which direction is it coming from + %dx_new, dy_new - which direction is it going + %m, n - dimensions of array + + %spt if statements - depend on old values + + + + %ept if statements - dpend on new values +end From 25dc3ed24d6d33f01f17165f8471a84236e99f69 Mon Sep 17 00:00:00 2001 From: ezcawi Date: Thu, 8 Jan 2015 14:55:30 -0500 Subject: [PATCH 20/50] Update big_path --- senior-design-code/big_path | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/senior-design-code/big_path b/senior-design-code/big_path index 9302daf..6575cb6 100644 --- a/senior-design-code/big_path +++ b/senior-design-code/big_path @@ -40,6 +40,11 @@ for i = 1:N-1 %for every level of resolution %link the paths together [current_path,res_probs_update] = link_paths(path_part,current_path,m_old, n_old, m_new,n_new,current_res_probs, lower_res_probs); + lower_res_probs = res_probs_update; % update the lower resolution for the next iteration end %compute the final path sum +for j = 1:length(current_path(:,1)) + path_sum = path_sum + prob_detect*res_probs_update(current_path(j,2), current_path(j,1)); + res_probs_update(current_path(j,2), current_path(j,1)) = res_probs_update(current_path(j,2), current_path(j,1))*(1-prob_detect); +end From 2d68c8fcc7c6a2f021670fff6e4e9830db1370c2 Mon Sep 17 00:00:00 2001 From: ezcawi Date: Thu, 8 Jan 2015 15:16:19 -0500 Subject: [PATCH 21/50] Update big_path --- senior-design-code/big_path | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/senior-design-code/big_path b/senior-design-code/big_path index 6575cb6..64e8bf4 100644 --- a/senior-design-code/big_path +++ b/senior-design-code/big_path @@ -42,7 +42,7 @@ for i = 1:N-1 %for every level of resolution [current_path,res_probs_update] = link_paths(path_part,current_path,m_old, n_old, m_new,n_new,current_res_probs, lower_res_probs); lower_res_probs = res_probs_update; % update the lower resolution for the next iteration end - +path_sum = 0; %compute the final path sum for j = 1:length(current_path(:,1)) path_sum = path_sum + prob_detect*res_probs_update(current_path(j,2), current_path(j,1)); From 97f30514a28c1d5d8c8c64df80dabf8014e537ed Mon Sep 17 00:00:00 2001 From: ezcawi Date: Mon, 12 Jan 2015 08:44:45 -0500 Subject: [PATCH 22/50] Update endpts.m --- senior-design-code/endpts.m | 42 +++++++++++++++++++++++++++++++++---- 1 file changed, 38 insertions(+), 4 deletions(-) diff --git a/senior-design-code/endpts.m b/senior-design-code/endpts.m index 4c2c3b1..87d1cdc 100644 --- a/senior-design-code/endpts.m +++ b/senior-design-code/endpts.m @@ -4,8 +4,42 @@ %m, n - dimensions of array %spt if statements - depend on old values - - - - %ept if statements - dpend on new values + if (dx_old == 0) && (dy_old ==0) %starting point of path, arbitrarily start at top + spt = [1, floor(n/2)]; + elseif (dx_old > 0) && (dy_old == 0) %coming from the left + spt = [floor(m/2),1]; + elseif (dx_old < 0) && (dy_old == 0) %coming from the right + spt = [floor(m/2),n]; + elseif (dx_old > 0) && (dy_old > 0) %coming from the bottom left + spt = [m,1]; + elseif (dx_old > 0) && (dy_old < 0) %coming from the top left + spt = [1,1]; + elseif (dx_old < 0) && (dy_old > 0) %coming from the bottom right + spt = [m,n]; + elseif (dx_old < 0) && (dy_old < 0) %coming from the top right + spt = [1,n]; + elseif (dx_old == 0) && (dy_old > 0) %coming from the bottom + spt = [m,floor(n/2)]; + elseif (dx_old == 0) && (dy_old < 0) %coming from the top + spt = [1,floor(n/2)]; + end + %ept if statements - depend on new values + + if (dx_new > 0) && (dy_new == 0) %going right + ept = [floor(m/2),n]; + elseif (dx_new < 0) && (dy_new == 0) %going left + ept = [floor(m/2),1]; + elseif (dx_new > 0) && (dy_new > 0) %going top right + ept = [1,n]; + elseif (dx_new > 0) && (dy_new < 0) %going bottom left + ept = [m,1]; + elseif (dx_new < 0) && (dy_new > 0) %going top left + ept = [1,1]; + elseif (dx_new < 0) && (dy_new < 0) %going bottom right + ept = [m,n]; + elseif (dx_new == 0) && (dy_new > 0) %going top + ept = [1,floor(n/2)]; + elseif (dx_new == 0) && (dy_new < 0) %going bottom + ept = [m,floor(n/2)]; + end end From 20fd01410f7635a11b2cf7289f42eaddc27ee197 Mon Sep 17 00:00:00 2001 From: ezcawi Date: Mon, 12 Jan 2015 09:13:14 -0500 Subject: [PATCH 23/50] Update link_paths.m --- senior-design-code/link_paths.m | 23 +++++++++++++++++++++-- 1 file changed, 21 insertions(+), 2 deletions(-) diff --git a/senior-design-code/link_paths.m b/senior-design-code/link_paths.m index 6c71e14..890808e 100644 --- a/senior-design-code/link_paths.m +++ b/senior-design-code/link_paths.m @@ -1,9 +1,28 @@ -function [current_path,higher_res_probs] = link_paths(path_part,current_path,m_old, n_old, m_new,n_new,current_res_probs, lower_res_probs) +function [current_path,res_probs_update] = link_paths(path_part,current_path,m_old, n_old, m_new,n_new,current_res_probs, lower_res_probs) %path_part - array of all the different path parts %m_old, n_old - dimensions of higher resolution array %m_new, n_new - dimensions of newer resolution arrays %current_res_probs - arrays of probabilities in current resolution %higher_res_probs - arrays of probabilities in lower resolution - + %create larger probability matrix + res_probs_update = zeros(m_old*m_new, n_old*n_new); + high_res_path = []; + for i = 1:length(current_path(:,1)) + %put probabilities in right spot in res_probs_update + current_point = current_path(i,:); + scaled_probs = lower_res_probs(current_point(2),current_point(1))*current_res_probs(:,:,i); % scales to make everything probabilities + res_probs_update((current_point(2)-1)*m_old*m_new+1:current_point(2)*m_old*m_new, (current_point(1)-1)*n_old*n_new+1:current_point(1)*n_old*n_new) = scaled_probs; + + %update path + %get path part that isn't zero + [good_rows,good_col] = find(path_part(:,:,i)>0); + good_path = path_part(good_rows,good_col); + %calculate indeces in larger matrix + index_update = ones(size(good_path)); + index_update(:,1) = (current_point(2)-1)*m_old*m_new + 1; + index_update(:,2) = (current_point(1)-1)*n_old*n_new + 1; + high_res_path = [high_res_path; (good_path + index_update)]; + end + current_path = high_res_path; % update current path at the end end From 5d34e3918b5b3b1f046be5de8eded8754dc68fd2 Mon Sep 17 00:00:00 2001 From: ezcawi Date: Mon, 12 Jan 2015 09:14:07 -0500 Subject: [PATCH 24/50] Update big_path --- senior-design-code/big_path | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/senior-design-code/big_path b/senior-design-code/big_path index 64e8bf4..912414c 100644 --- a/senior-design-code/big_path +++ b/senior-design-code/big_path @@ -16,13 +16,15 @@ current_path = starting_path; lower_res_probs = top_level_probs; [m_old, n_old] = size(top_level_probs); for i = 1:N-1 %for every level of resolution + [m_new,n_new] = input('Size of current_resolution blocks'); % expecting a two dimensional array + current_res_probs = zeros(m_new,n_new,length(current_path(:,1))); + path_part = zeros(m_new*n_new, 2,length(current_path(:,1)); for j = 1: length(current_path(:,1))-1 %load this node's next level - would like to load entire %resolutions probabilities and segment but segmentation indexing %is proving difficult to wrap my head around. current_res_probs(:,:,j) = load(input('next probability array','s')); C_area = floor(C*lower_res_probs(current_path(j,2), current_path(j,1))*numel(current_res_probs)); %current cost * probability subject is in an area = proportional time spent in that area * number of cells = number of cells in new resolution we can cover - [m_new,n_new] = size(current_res_probs); %calculate starting and ending points for this level from previous %path resolution dx_new = current_path(j+1,1) - current_path(j,1); %next direction, determines ending point @@ -35,14 +37,14 @@ for i = 1:N-1 %for every level of resolution dy_old = 0; end [spt,ept] = endpts(dx_old, dy_old, dx_new, dy_new,m_new, n_new); - [path_part(:,:,j), path_sum] = hill_climbing(m_new,n_new,current_res_probs,prob_detect,C_area,spt,ept); + [path_part(1:C_area,1:C_area,j), path_sum] = hill_climbing(m_new,n_new,current_res_probs,prob_detect,C_area,spt,ept); end %link the paths together [current_path,res_probs_update] = link_paths(path_part,current_path,m_old, n_old, m_new,n_new,current_res_probs, lower_res_probs); lower_res_probs = res_probs_update; % update the lower resolution for the next iteration end -path_sum = 0; + %compute the final path sum for j = 1:length(current_path(:,1)) path_sum = path_sum + prob_detect*res_probs_update(current_path(j,2), current_path(j,1)); From 0f550c335607dc1f690cd1f5bf7a04820f451eba Mon Sep 17 00:00:00 2001 From: ezcawi Date: Mon, 12 Jan 2015 09:14:36 -0500 Subject: [PATCH 25/50] Update big_path --- senior-design-code/big_path | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/senior-design-code/big_path b/senior-design-code/big_path index 912414c..3e4c748 100644 --- a/senior-design-code/big_path +++ b/senior-design-code/big_path @@ -37,7 +37,7 @@ for i = 1:N-1 %for every level of resolution dy_old = 0; end [spt,ept] = endpts(dx_old, dy_old, dx_new, dy_new,m_new, n_new); - [path_part(1:C_area,1:C_area,j), path_sum] = hill_climbing(m_new,n_new,current_res_probs,prob_detect,C_area,spt,ept); + [path_part(1:C_area,:,j), path_sum] = hill_climbing(m_new,n_new,current_res_probs,prob_detect,C_area,spt,ept); end %link the paths together From 5c082394111d1a4bfc3c94cc1c8423ef050bad0d Mon Sep 17 00:00:00 2001 From: ezcawi Date: Mon, 12 Jan 2015 09:18:47 -0500 Subject: [PATCH 26/50] Update big_path --- senior-design-code/big_path | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/senior-design-code/big_path b/senior-design-code/big_path index 3e4c748..286b59a 100644 --- a/senior-design-code/big_path +++ b/senior-design-code/big_path @@ -10,7 +10,7 @@ N = input('How many layers of Resolution?'); %integer prob_detect = input('What is the UAVs probability of detection?'); % number between 0 and 1 top_level_probs = load(input('Top Level File Name','s')); %filename C = input('What percentage of the total area can the UAV cover?'); % needs a number between 0 and 1 -starting_path = input('What is the starting path through the top level?') % needs an n x 2 array +starting_path = input('What is the starting path through the top level?'); % needs an n x 2 array %need to check inputs for final version current_path = starting_path; lower_res_probs = top_level_probs; @@ -18,7 +18,7 @@ lower_res_probs = top_level_probs; for i = 1:N-1 %for every level of resolution [m_new,n_new] = input('Size of current_resolution blocks'); % expecting a two dimensional array current_res_probs = zeros(m_new,n_new,length(current_path(:,1))); - path_part = zeros(m_new*n_new, 2,length(current_path(:,1)); + path_part = zeros(m_new*n_new, 2,length(current_path(:,1))); for j = 1: length(current_path(:,1))-1 %load this node's next level - would like to load entire %resolutions probabilities and segment but segmentation indexing From fc0e79647fdd51b978e4d5d28640f924c20f1d33 Mon Sep 17 00:00:00 2001 From: ezcawi Date: Mon, 12 Jan 2015 09:23:12 -0500 Subject: [PATCH 27/50] Update big_path --- senior-design-code/big_path | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/senior-design-code/big_path b/senior-design-code/big_path index 286b59a..083f8ce 100644 --- a/senior-design-code/big_path +++ b/senior-design-code/big_path @@ -8,7 +8,7 @@ N = input('How many layers of Resolution?'); %integer prob_detect = input('What is the UAVs probability of detection?'); % number between 0 and 1 -top_level_probs = load(input('Top Level File Name','s')); %filename +top_level_probs = load(input('Top Level File Name','s'),...); %filename C = input('What percentage of the total area can the UAV cover?'); % needs a number between 0 and 1 starting_path = input('What is the starting path through the top level?'); % needs an n x 2 array %need to check inputs for final version @@ -23,7 +23,7 @@ for i = 1:N-1 %for every level of resolution %load this node's next level - would like to load entire %resolutions probabilities and segment but segmentation indexing %is proving difficult to wrap my head around. - current_res_probs(:,:,j) = load(input('next probability array','s')); + current_res_probs(:,:,j) = load(input('next probability array','s'),...); C_area = floor(C*lower_res_probs(current_path(j,2), current_path(j,1))*numel(current_res_probs)); %current cost * probability subject is in an area = proportional time spent in that area * number of cells = number of cells in new resolution we can cover %calculate starting and ending points for this level from previous %path resolution From e880ab0cf179616a08b22f8b533bd72e65aa2e5b Mon Sep 17 00:00:00 2001 From: ezcawi Date: Mon, 12 Jan 2015 09:28:03 -0500 Subject: [PATCH 28/50] Update big_path --- senior-design-code/big_path | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/senior-design-code/big_path b/senior-design-code/big_path index 083f8ce..286b59a 100644 --- a/senior-design-code/big_path +++ b/senior-design-code/big_path @@ -8,7 +8,7 @@ N = input('How many layers of Resolution?'); %integer prob_detect = input('What is the UAVs probability of detection?'); % number between 0 and 1 -top_level_probs = load(input('Top Level File Name','s'),...); %filename +top_level_probs = load(input('Top Level File Name','s')); %filename C = input('What percentage of the total area can the UAV cover?'); % needs a number between 0 and 1 starting_path = input('What is the starting path through the top level?'); % needs an n x 2 array %need to check inputs for final version @@ -23,7 +23,7 @@ for i = 1:N-1 %for every level of resolution %load this node's next level - would like to load entire %resolutions probabilities and segment but segmentation indexing %is proving difficult to wrap my head around. - current_res_probs(:,:,j) = load(input('next probability array','s'),...); + current_res_probs(:,:,j) = load(input('next probability array','s')); C_area = floor(C*lower_res_probs(current_path(j,2), current_path(j,1))*numel(current_res_probs)); %current cost * probability subject is in an area = proportional time spent in that area * number of cells = number of cells in new resolution we can cover %calculate starting and ending points for this level from previous %path resolution From e75223a08ce36dfe17863cc53198906f7fd26cfd Mon Sep 17 00:00:00 2001 From: ezcawi Date: Mon, 12 Jan 2015 09:32:26 -0500 Subject: [PATCH 29/50] Update big_path --- senior-design-code/big_path | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/senior-design-code/big_path b/senior-design-code/big_path index 286b59a..12d32a8 100644 --- a/senior-design-code/big_path +++ b/senior-design-code/big_path @@ -8,7 +8,7 @@ N = input('How many layers of Resolution?'); %integer prob_detect = input('What is the UAVs probability of detection?'); % number between 0 and 1 -top_level_probs = load(input('Top Level File Name','s')); %filename +top_level_probs = struct2cell(load(input('Top Level File Name','s'))); %filename C = input('What percentage of the total area can the UAV cover?'); % needs a number between 0 and 1 starting_path = input('What is the starting path through the top level?'); % needs an n x 2 array %need to check inputs for final version @@ -16,14 +16,15 @@ current_path = starting_path; lower_res_probs = top_level_probs; [m_old, n_old] = size(top_level_probs); for i = 1:N-1 %for every level of resolution - [m_new,n_new] = input('Size of current_resolution blocks'); % expecting a two dimensional array + res_size = input('Size of current_resolution blocks'); % expecting a two dimensional array + m_new = res_size(1); n_new = res_size(2); current_res_probs = zeros(m_new,n_new,length(current_path(:,1))); path_part = zeros(m_new*n_new, 2,length(current_path(:,1))); for j = 1: length(current_path(:,1))-1 %load this node's next level - would like to load entire %resolutions probabilities and segment but segmentation indexing %is proving difficult to wrap my head around. - current_res_probs(:,:,j) = load(input('next probability array','s')); + current_res_probs(:,:,j) = struct2cell(load(input('next probability array','s'))); C_area = floor(C*lower_res_probs(current_path(j,2), current_path(j,1))*numel(current_res_probs)); %current cost * probability subject is in an area = proportional time spent in that area * number of cells = number of cells in new resolution we can cover %calculate starting and ending points for this level from previous %path resolution From 70fba4b94927833d843d45b32c0ab3930752f8e6 Mon Sep 17 00:00:00 2001 From: ezcawi Date: Mon, 12 Jan 2015 09:49:58 -0500 Subject: [PATCH 30/50] Update endpts.m --- senior-design-code/endpts.m | 28 ++++++++++++++-------------- 1 file changed, 14 insertions(+), 14 deletions(-) diff --git a/senior-design-code/endpts.m b/senior-design-code/endpts.m index 87d1cdc..43ab9d5 100644 --- a/senior-design-code/endpts.m +++ b/senior-design-code/endpts.m @@ -5,41 +5,41 @@ %spt if statements - depend on old values if (dx_old == 0) && (dy_old ==0) %starting point of path, arbitrarily start at top - spt = [1, floor(n/2)]; + spt = [floor(n/2),1]; elseif (dx_old > 0) && (dy_old == 0) %coming from the left - spt = [floor(m/2),1]; + spt = [1, floor(m/2)]; elseif (dx_old < 0) && (dy_old == 0) %coming from the right - spt = [floor(m/2),n]; + spt = [n,floor(m/2)]; elseif (dx_old > 0) && (dy_old > 0) %coming from the bottom left - spt = [m,1]; + spt = [1,m]; elseif (dx_old > 0) && (dy_old < 0) %coming from the top left spt = [1,1]; elseif (dx_old < 0) && (dy_old > 0) %coming from the bottom right - spt = [m,n]; + spt = [n,m]; elseif (dx_old < 0) && (dy_old < 0) %coming from the top right - spt = [1,n]; + spt = [n,1]; elseif (dx_old == 0) && (dy_old > 0) %coming from the bottom - spt = [m,floor(n/2)]; + spt = [floor(n/2),m]; elseif (dx_old == 0) && (dy_old < 0) %coming from the top - spt = [1,floor(n/2)]; + spt = [floor(n/2),1]; end %ept if statements - depend on new values if (dx_new > 0) && (dy_new == 0) %going right - ept = [floor(m/2),n]; + ept = [n,floor(m/2)]; elseif (dx_new < 0) && (dy_new == 0) %going left - ept = [floor(m/2),1]; + ept = [1,floor(m/2)]; elseif (dx_new > 0) && (dy_new > 0) %going top right ept = [1,n]; elseif (dx_new > 0) && (dy_new < 0) %going bottom left - ept = [m,1]; + ept = [1,m]; elseif (dx_new < 0) && (dy_new > 0) %going top left ept = [1,1]; elseif (dx_new < 0) && (dy_new < 0) %going bottom right - ept = [m,n]; + ept = [n,m]; elseif (dx_new == 0) && (dy_new > 0) %going top - ept = [1,floor(n/2)]; + ept = [floor(n/2),1]; elseif (dx_new == 0) && (dy_new < 0) %going bottom - ept = [m,floor(n/2)]; + ept = [floor(n/2),m]; end end From 672feb61fe5a85654af3fa53add3f1e65376ed79 Mon Sep 17 00:00:00 2001 From: ezcawi Date: Wed, 14 Jan 2015 18:49:23 -0500 Subject: [PATCH 31/50] Update add_node.m --- senior-design-code/add_node.m | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/senior-design-code/add_node.m b/senior-design-code/add_node.m index fcc30bb..9839bd8 100644 --- a/senior-design-code/add_node.m +++ b/senior-design-code/add_node.m @@ -8,7 +8,6 @@ %get the current point on the path and it's immediate neighbors new_point_choice = zeros(size(current_path)); best_probs = zeros(length(current_path(1,:)),1); - for k = 2:(length(current_path(:,1))-1) current_xy = current_path(k,:); next_xy = current_path(k+1,:); %next point on path @@ -28,6 +27,7 @@ for p = 1: length(valid_choices(:,1)) n_probs(p) = current_probs(valid_choices(p,2), valid_choices(p,1)); end + a = find(n_probs == max(n_probs)); new_point_choice(k,:) = valid_choices(a(1),:); best_probs(k) = current_probs(new_point_choice(k,2), new_point_choice(k,1)); From c68bea7bb6538e49be41bcd9a0426b675620255e Mon Sep 17 00:00:00 2001 From: ezcawi Date: Wed, 14 Jan 2015 18:49:59 -0500 Subject: [PATCH 32/50] Update greedy_algorithm.m --- senior-design-code/greedy_algorithm.m | 108 ++++++++++---------------- 1 file changed, 39 insertions(+), 69 deletions(-) diff --git a/senior-design-code/greedy_algorithm.m b/senior-design-code/greedy_algorithm.m index e601585..6d62511 100644 --- a/senior-design-code/greedy_algorithm.m +++ b/senior-design-code/greedy_algorithm.m @@ -1,75 +1,45 @@ -%greedy path finding algorithm for search and rescue -%first, load in the array of probabilities, for this example i just picked -%a matrix from a uniform distribution between 0 and 1 and then turn that -%into probabilities that sum to 1 -function [best_path,best_sum] = greedy_algorithm(m,n,probs,prob_detect) -%starting off with 4x4 matrix for simplicity - %turn the example set into probabilities. -%set a cost for the search, here it is arbitrarily 75% of cells that the quad can cover, -%can determine this through flight time, etc -tic -num_cells = floor(m*n*0.75); +function greedy_path = greedy_algorithm(probs,prob_detect,spt,C) +%probs - probability array +%prob_detect - probability of detection +%spt - starting point +%C- percentage of cells coverable +%performs a greedy algorithm given starting point and probability of +%detection, no fixed ending point +[m,n] = size(probs); +num_cells = C; -path_sum = zeros(m,n); -current_x = zeros(m,n,num_cells); -current_y = zeros(m,n,num_cells); %this will keep track of the path for each starting point -%to find the best path moving through only some of the cells,we want to know the best starting point -for x = 1:1:n - for y = 1:1:m %first two loops desginate going through the different starting points - search_value_temp = probs; %going to be changing this matrix so i set the thing - current_x(y,x,1) = x; - current_y(y,x,1) = y; - path_sum(y,x) =prob_detect*search_value_temp(y, x);%initialize path probability sum, using probability of detection because we might miss - search_value_temp(current_y(y,x,1), current_x(y,x,1)) = search_value_temp(current_y(y,x,1), current_x(y,x,1))*(1 - prob_detect); %subtracting some probability from the first cell - - for k = 2:1:num_cells - - [neighbors_x, neighbors_y] = get_neighbors(current_x(y,x,k-1), current_y(y,x,k-1), n, n); - current_neighbors = [neighbors_x neighbors_y]; +search_value_temp = probs; %going to be changing this matrix so i set the thing +current_x = zeros(num_cells,1); +current_y = zeros(num_cells,1); +current_x(1) = spt(1); +current_y(1) = spt(2); +search_value_temp(current_y(1), current_x(1)) = search_value_temp(current_y(1), current_x(1))*(1 - prob_detect); %subtracting some probability from the first cell - neighbor_prob = zeros(size(neighbors_x)); - for i = 1:length(neighbors_x) - neighbor_prob(i) = search_value_temp(current_neighbors(i,2),current_neighbors(i,1)); - end - best_neighbor = find(neighbor_prob == max(neighbor_prob)); - next_x = current_neighbors(best_neighbor,1); - next_y = current_neighbors(best_neighbor,2); - if length(next_x)>1 %if there are two equal probability cells in the neighbors list just pick one randomly - p = randperm(length(next_x), 1); - current_x(y,x,k) = next_x(p); - current_y(y,x,k) = next_y(p); - else - current_x(y,x,k) = next_x; - current_y(y,x,k) = next_y; - end; - path_sum(y,x) = path_sum(y,x) + prob_detect*search_value_temp(current_y(y,x,k), current_x(y,x,k)); %update the path sum - search_value_temp(current_y(y,x,k), current_x(y,x,k)) = search_value_temp(current_y(y,x,k), current_x(y,x,k))*(1 - prob_detect); - end +for k = 2:1:num_cells + [neighbors_x, neighbors_y] = get_neighbors(current_x(k-1), current_y(k-1), m, n); + current_neighbors = [neighbors_x neighbors_y]; + neighbor_prob = zeros(size(neighbors_x)); + for i = 1:length(neighbors_x) + neighbor_prob(i) = search_value_temp(current_neighbors(i,2),current_neighbors(i,1)); end + best_neighbor = find(neighbor_prob == max(neighbor_prob)); + next_x = current_neighbors(best_neighbor,1); + next_y = current_neighbors(best_neighbor,2); + if length(next_x)>1 %if there are two equal probability cells in the neighbors list just pick one randomly + p = randperm(length(next_x), 1); + current_x(k) = next_x(p); + current_y(k) = next_y(p); + else + current_x(k) = next_x; + current_y(k) = next_y; + end; + %update search probability by (1-probdetect) * p to reflect searching + %that cell + search_value_temp(current_y(k), current_x(k)) = search_value_temp(current_y(k), current_x(k))*(1 - prob_detect); end -[best_start_y, best_start_x] = find(path_sum == max(max(path_sum))); -if length(best_start_x)>1 - p = randperm(length(best_start_x),1); - best_start_x = best_start_x(p); - best_start_y = best_start_y(p); -end -best_sum = path_sum(best_start_y, best_start_x); -best_pathx = zeros(num_cells,1); -best_pathy = zeros(num_cells,1); -best_pathtempx = current_x(best_start_y, best_start_x, :); -best_pathtempy = current_y(best_start_y, best_start_x, :); -for i =1:num_cells - best_pathx(i) = best_pathtempx(1,1,i); - best_pathy(i) = best_pathtempy(1,1,i); -end -best_path = [best_pathx best_pathy]; -size(best_path) -% display(sprintf('the best greedy path sum is: %d',best_sum)) -% display(sprintf('the best greedy starting point is %d , %d', best_start_x, best_start_y)) -% display(sprintf('the best greedy path found is: ')); -% display(best_path) -toc +greedy_path = [current_x current_y]; + + + -%needed, vectorize for loops to make the program scaleable to a larger array, -%firgure out higher level algorithm to get the effort put into cells From 265336aaa5c1c48f6217afd12cb42a58c8d25229 Mon Sep 17 00:00:00 2001 From: ezcawi Date: Wed, 14 Jan 2015 18:51:13 -0500 Subject: [PATCH 33/50] Update big_path --- senior-design-code/big_path | 83 ++++++++++++++++++++++++------------- 1 file changed, 54 insertions(+), 29 deletions(-) diff --git a/senior-design-code/big_path b/senior-design-code/big_path index 12d32a8..95cca4f 100644 --- a/senior-design-code/big_path +++ b/senior-design-code/big_path @@ -5,49 +5,74 @@ %the same size, want to spend time in every search region %inputs from user: probability files, number of levels to go down - N % the inputs should be interfaced with the gui somehow - -N = input('How many layers of Resolution?'); %integer -prob_detect = input('What is the UAVs probability of detection?'); % number between 0 and 1 -top_level_probs = struct2cell(load(input('Top Level File Name','s'))); %filename -C = input('What percentage of the total area can the UAV cover?'); % needs a number between 0 and 1 -starting_path = input('What is the starting path through the top level?'); % needs an n x 2 array +clear all; close all; clc; +N = input('How many layers of Resolution?: '); %integer +prob_detect = input('What is the UAVs probability of detection?: '); % number between 0 and 1 +top_level_struct = load(input('Top Level File Name: ','s')); %filename +fname = fieldnames(top_level_struct); +top_level_probs = top_level_struct.(fname{1}); % have to do this to convert from struct to array +C = input('What percentage of the total area can the UAV cover?: '); % needs a number between 0 and 1 +starting_path = input('What is the starting path through the top level?: '); % needs an n x 2 array %need to check inputs for final version current_path = starting_path; lower_res_probs = top_level_probs; [m_old, n_old] = size(top_level_probs); -for i = 1:N-1 %for every level of resolution - res_size = input('Size of current_resolution blocks'); % expecting a two dimensional array - m_new = res_size(1); n_new = res_size(2); - current_res_probs = zeros(m_new,n_new,length(current_path(:,1))); +for i = 2:N %for every level of resolution + str = sprintf('Resolution level %d ',i); + disp(str); + current_res_struct = load(input('next resolution probability array: ','s')); + fname = fieldnames(current_res_struct); + current_res_probs = current_res_struct.(fname{1}); + [m_new, n_new] = size(current_res_probs); path_part = zeros(m_new*n_new, 2,length(current_path(:,1))); - for j = 1: length(current_path(:,1))-1 - %load this node's next level - would like to load entire - %resolutions probabilities and segment but segmentation indexing - %is proving difficult to wrap my head around. - current_res_probs(:,:,j) = struct2cell(load(input('next probability array','s'))); + for j = 1: length(current_path(:,1)) + %figure out which segment of the path to use for path planning, + %assumes that m_new is a multiple of m_old, i.e. each cell in the + %old resolution is expanding into an integer number of cells in the + %new resolution + m_current = m_new/m_old; + n_current = n_new/n_old; + current_probs = current_res_probs((current_path(j,2)-1)*m_current+1:current_path(j,2)*m_current, (current_path(j,1)-1)*n_current+1:current_path(j,1)*n_current); + C_area = floor(C*lower_res_probs(current_path(j,2), current_path(j,1))*numel(current_res_probs)); %current cost * probability subject is in an area = proportional time spent in that area * number of cells = number of cells in new resolution we can cover %calculate starting and ending points for this level from previous %path resolution - dx_new = current_path(j+1,1) - current_path(j,1); %next direction, determines ending point - dy_new = current_path(j+1,2) - current_path(j,2); - if j>2 - dx_old = current_path(j,1) - current_path(j-1,1); %previous direction, determines starting point + if j=2 + dx_old = current_path(j,1) - current_path(j-1,1);%previous direction, determines starting point dy_old = current_path(j,2) - current_path(j-1,2); - else + else %j = 1 dx_old = 0; dy_old = 0; end - [spt,ept] = endpts(dx_old, dy_old, dx_new, dy_new,m_new, n_new); - [path_part(1:C_area,:,j), path_sum] = hill_climbing(m_new,n_new,current_res_probs,prob_detect,C_area,spt,ept); + [spt,ept] = endpts(dx_old, dy_old, dx_new, dy_new,m_current, n_current); + if ((ept(1) == 0)&&(ept(2) == 0)) %at the end of the path, do a greedy algorithm + disp('end of current path') + path_part(1:C_area,:,j) = greedy_algorithm(current_probs, prob_detect, spt, C_area); + elseif ((dx_old == 0)&&(dy_old == 0)) %this happens at the beginning, do greedy from the endpoint and then reverse the path + disp('start of current path'); + path_temp = greedy_algorithm(current_probs,prob_detect,ept, C_area); % ending point is starting point of the path + path_part(1:C_area,:,j) = flipud(path_temp);%flips up/down so that the ending point is hte last part of the array + else % the normal case, hill climbing with fixed start and end points + disp('hill climbing'); + path_part(1:C_area,:,j) = hill_climbing(m_current,n_current,current_probs,prob_detect,C_area,spt,ept); + end end - %link the paths together - - [current_path,res_probs_update] = link_paths(path_part,current_path,m_old, n_old, m_new,n_new,current_res_probs, lower_res_probs); - lower_res_probs = res_probs_update; % update the lower resolution for the next iteration + %link the paths together to update current path + current_path = link_paths(path_part,current_path, m_current,n_current); + lower_res_probs = current_res_probs; % update the lower resolution for the next iteration end - +disp('done'); %compute the final path sum +path_sum = 0; for j = 1:length(current_path(:,1)) - path_sum = path_sum + prob_detect*res_probs_update(current_path(j,2), current_path(j,1)); - res_probs_update(current_path(j,2), current_path(j,1)) = res_probs_update(current_path(j,2), current_path(j,1))*(1-prob_detect); + path_sum = path_sum + prob_detect*current_res_probs(current_path(j,2), current_path(j,1)); + current_res_probs(current_path(j,2), current_path(j,1)) = current_res_probs(current_path(j,2), current_path(j,1))*(1-prob_detect); end + From b868d90ce95331c015419e7eac44d56e3e211633 Mon Sep 17 00:00:00 2001 From: ezcawi Date: Wed, 14 Jan 2015 18:52:06 -0500 Subject: [PATCH 34/50] Update hill_climbing.m --- senior-design-code/hill_climbing.m | 17 ++++++++++------- 1 file changed, 10 insertions(+), 7 deletions(-) diff --git a/senior-design-code/hill_climbing.m b/senior-design-code/hill_climbing.m index fc16e3f..bc5cb47 100644 --- a/senior-design-code/hill_climbing.m +++ b/senior-design-code/hill_climbing.m @@ -1,18 +1,19 @@ -function [best_path,best_path_sum] = hill_climbing(m,n,probs,prob_detect,C,spt, ept) +function [hill_path] = hill_climbing(m,n,probs,prob_detect,C,spt, ept) %m,n - dimensions of matrix %probs - array of probability that the subject is in location %prob_detect - probability that searcher will see subject %C - total cost i.e. number of cells that can be visited %spt, ept - starting and ending points, respectively - tic i_best_path = zeros(C,2,10); %for storing the best path i_path_sum = zeros(1,10); path_sum_start = zeros(1,10); + path_start = zeros(C,2,10); for i = 1:10 probs_current = probs; %first find random starting path with length less than or equal to C l = C+1; - while l>C + current_path = spt; + while ((l>C)||(length(current_path(:,1))== 2)) current_path = gen_rand_path(spt,ept,m,n); l = length(current_path(:,1)); end @@ -23,8 +24,10 @@ probs_current(current_path(j,2), current_path(j,1)) = probs_current(current_path(j,2), current_path(j,1))*(1-prob_detect); end path_sum_start(i) = path_sum; + + path_start(1:length(current_path(:,1)),:,i) = current_path; %hill climbing loop: two rules: add a node, replace a node - for a = 1:100 + for a = 1:200 if length(current_path(:,1)) Date: Wed, 14 Jan 2015 18:53:12 -0500 Subject: [PATCH 35/50] Update link_paths.m --- senior-design-code/link_paths.m | 29 +++++++++++------------------ 1 file changed, 11 insertions(+), 18 deletions(-) diff --git a/senior-design-code/link_paths.m b/senior-design-code/link_paths.m index 890808e..350564f 100644 --- a/senior-design-code/link_paths.m +++ b/senior-design-code/link_paths.m @@ -1,28 +1,21 @@ -function [current_path,res_probs_update] = link_paths(path_part,current_path,m_old, n_old, m_new,n_new,current_res_probs, lower_res_probs) +function [linked_path] = link_paths(path_part,current_path, m_current,n_current) %path_part - array of all the different path parts - %m_old, n_old - dimensions of higher resolution array - %m_new, n_new - dimensions of newer resolution arrays - %current_res_probs - arrays of probabilities in current resolution - %higher_res_probs - arrays of probabilities in lower resolution - - %create larger probability matrix - res_probs_update = zeros(m_old*m_new, n_old*n_new); + %current_path - path of the lower resolution + %m_current, n_current - dimensions of newer resolution arrays + %stitches path together from path parts and current path high_res_path = []; for i = 1:length(current_path(:,1)) - %put probabilities in right spot in res_probs_update current_point = current_path(i,:); - scaled_probs = lower_res_probs(current_point(2),current_point(1))*current_res_probs(:,:,i); % scales to make everything probabilities - res_probs_update((current_point(2)-1)*m_old*m_new+1:current_point(2)*m_old*m_new, (current_point(1)-1)*n_old*n_new+1:current_point(1)*n_old*n_new) = scaled_probs; - - %update path %get path part that isn't zero - [good_rows,good_col] = find(path_part(:,:,i)>0); - good_path = path_part(good_rows,good_col); + p = path_part(:,:,i); + [good_rows,good_col] = find(p ~= 0); + good_path = p(1:length(good_rows)/2,:) %calculate indeces in larger matrix index_update = ones(size(good_path)); - index_update(:,1) = (current_point(2)-1)*m_old*m_new + 1; - index_update(:,2) = (current_point(1)-1)*n_old*n_new + 1; + index_update(:,1) = (current_point(1)-1)*m_current; + index_update(:,2) = (current_point(2)-1)*n_current; + good_path + index_update high_res_path = [high_res_path; (good_path + index_update)]; end - current_path = high_res_path; % update current path at the end + linked_path = high_res_path; % update current path at the end end From 348487c2065910c572080fbe22ef678baeb23402 Mon Sep 17 00:00:00 2001 From: ezcawi Date: Wed, 14 Jan 2015 18:59:20 -0500 Subject: [PATCH 36/50] Update endpts.m --- senior-design-code/endpts.m | 21 +++++++++++---------- 1 file changed, 11 insertions(+), 10 deletions(-) diff --git a/senior-design-code/endpts.m b/senior-design-code/endpts.m index 43ab9d5..d89d810 100644 --- a/senior-design-code/endpts.m +++ b/senior-design-code/endpts.m @@ -24,22 +24,23 @@ spt = [floor(n/2),1]; end %ept if statements - depend on new values - if (dx_new > 0) && (dy_new == 0) %going right ept = [n,floor(m/2)]; elseif (dx_new < 0) && (dy_new == 0) %going left ept = [1,floor(m/2)]; - elseif (dx_new > 0) && (dy_new > 0) %going top right - ept = [1,n]; - elseif (dx_new > 0) && (dy_new < 0) %going bottom left - ept = [1,m]; - elseif (dx_new < 0) && (dy_new > 0) %going top left - ept = [1,1]; - elseif (dx_new < 0) && (dy_new < 0) %going bottom right + elseif (dx_new > 0) && (dy_new > 0) %going bottom right ept = [n,m]; - elseif (dx_new == 0) && (dy_new > 0) %going top + elseif (dx_new > 0) && (dy_new < 0) %going top left + ept = [1,1]; + elseif (dx_new < 0) && (dy_new > 0) %going bottom left + ept = [1,m]; + elseif (dx_new < 0) && (dy_new < 0) %going top right + ept = [n,1]; + elseif (dx_new == 0) && (dy_new < 0) %going top ept = [floor(n/2),1]; - elseif (dx_new == 0) && (dy_new < 0) %going bottom + elseif (dx_new == 0) && (dy_new > 0) %going bottom ept = [floor(n/2),m]; + elseif (dx_new == 0) && (dy_new ==0) %hit end of path, will do greedy + ept = [0,0]; end end From d6baf9806ed97c46fde1bdf8b6c18c45cc18b3a6 Mon Sep 17 00:00:00 2001 From: ezcawi Date: Mon, 19 Jan 2015 08:58:42 -0500 Subject: [PATCH 37/50] Update and rename LatLon_to_m.py to array_to_meters_lon.py --- senior-design-code/LatLon_to_m.py | 42 ------------------- senior-design-code/array_to_meters_lon.py | 49 +++++++++++++++++++++++ 2 files changed, 49 insertions(+), 42 deletions(-) delete mode 100644 senior-design-code/LatLon_to_m.py create mode 100644 senior-design-code/array_to_meters_lon.py diff --git a/senior-design-code/LatLon_to_m.py b/senior-design-code/LatLon_to_m.py deleted file mode 100644 index 34d7082..0000000 --- a/senior-design-code/LatLon_to_m.py +++ /dev/null @@ -1,42 +0,0 @@ -import math -import numpy as np -def array_to_meters(probs, center_lat, center_lon, heading,pixel_res) -#probs - numpy array of probabilities, will need size of the array -#center_lat, center_lon - central coordinates of the array -#heading - angle of orientation of the array, measured from north/y axis -#pixel_res - number of meters per pixel - shp = probs.shape - xs = np.zeros(shp) - ys = np.zeros(shp) - #convert central point to meters - xy_c = coord_to_meters(center_lat,center_lon) - c_index = [shp[0]/2 - 1,shp[1]/2 -1] - for i in range(shp[0]) - dy_index = i - c_index[0] - dy = dy_index*pixel_res - ys[i][:] = xy_c[1] + dy - for j in range(shp[1]) - dx_index = j - c_index[1] - dx = dx_index*pixel_res - xs[:][j] = xy_c[0] + dx - #now need to rotate each x,y by the heading - rot_matrix = np.array([[math.cos(-1*heading), -1*math.sin(-1*heading)],[math.sin(-1*heading), math.cos(-1*heading)]]) - xy_rot = np.zeros((shp[0],shp[1],2)) - for i in range(shp[0]) - for j in range(shp[1]) - xy_vec = np.array([[xs[i][j]],[ys[i][j]]) - xy_rot[i][j][:] = np.dot(rot_matrix,xy_vec) - return xy_rot #returns array of x,y coordinates for each element of the array -def coord_to_meters(lat,lon) - lat_degrees = float(lat) - long_degrees = float(lng) - - lat_radians = float(lat_degrees)*(math.pi/180) - long_radians = float(long_degrees)*(math.pi/180) - - R = 6378137 - XY = [0,0] - XY[0] = R*math.cos(lat_radians)*math.cos(long_radians) - XY[1] = R*math.cos(lat_radians)*math.sin(long_radians) - Z = R*math.sin(lat_radians) #not sure if this is right/needed if pixel resolution is known - return XY diff --git a/senior-design-code/array_to_meters_lon.py b/senior-design-code/array_to_meters_lon.py new file mode 100644 index 0000000..7190ec7 --- /dev/null +++ b/senior-design-code/array_to_meters_lon.py @@ -0,0 +1,49 @@ +import math +import numpy as np +def array_to_meters_lon(arr, center_lat, center_lon,HFOV, VFOV, L,theta, alt,met_or_lon) +#arr - numpy array- probabilities or a portion of the picture with possible object in it +#center_lat, center_lon - central coordinates of the array +#hfov,vfov- horizontal and vertical fields of view +#L - focal length of camera +#theta - yaw of quad, in radians +#alt - altitude of quad +#met_or_lon - does the user want global coordinates in meters or longitude, 0 for meters, 1 for lat/lon +#this function takes arrays, along with a central point-for example tagged from a picture, and returns x,y,z or lat lon coordinates +#equations based off of "Robot Navigation", 2011 by Gerald Cook + shp = arr.shape + xs = np.zeros(shp) + ys = np.zeros(shp) + N = shp[0] + M = shp[1] + #convert central point to meters + c_index = [shp[0]/2 - 1,shp[1]/2 -1] + #convert to cartesian + for i in range(shp[0]) #rows of array + ys[i][:] = float(L* (N+1.0-2.0*i)/(N-1.0) * np.tan(VFOV/2)) + for j in range(shp[1]) + xs[:][j] = float(L* (2*j - M -1.0)/(M-1.0) * np.tan(HFOV/2)) + #azimuth and zenith + az = np.atan(-1*xs/L) + zen = np.atan(np.divide(ys, np.sqrt(L*L + np.square(xs)))) + #camera coordinates + xc = Z*np.tan(az) + yc = Z*np.tan(zen) + #world coordinates, 2x2 rotation matrix * [xc,yc], rotates by yaw + xw = xc*np.cos(theta) - yc*np.sin(theta) + yw = xc*np.sin(theta) + yc*np.cos(theta) + #global coordinates in meters + R = 6371000 #average radius of earth, equatorial radius is 6378137 meters + lat_rad = center_lat*math.pi/180.0 + lon_rad = center_lon*math.pi/180.0 + xg = np.add(R*np.cos(lat_rad), xw) + yg = np.add(R*np.cos(lat_rad)*np.sin(lon_rad),yw) + zg = R*np.sin(lat_rad)+Z + if met_or_lon == 0 #users + return xg,yg,zg + elif met_or_lon == 1 + lats = np.atan(np.divide(Zg,np.sqrt(np.add(np.square(xg),np.square(yg))))) + lons = np.atan(np.divide(yg,xg))#in design document it's yw/xw but I think it should be yg/xg + return lats, lons #returns array of x,y coordinates for each element of the array + else + return 'bad input' + From d9c9b5a9b4acd63632472e553acd5ed109ab0af9 Mon Sep 17 00:00:00 2001 From: ezcawi Date: Tue, 20 Jan 2015 17:48:34 -0500 Subject: [PATCH 38/50] fixed starting point logic causing issue in link_path --- senior-design-code/endpts.m | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/senior-design-code/endpts.m b/senior-design-code/endpts.m index d89d810..a699c3e 100644 --- a/senior-design-code/endpts.m +++ b/senior-design-code/endpts.m @@ -4,23 +4,23 @@ %m, n - dimensions of array %spt if statements - depend on old values - if (dx_old == 0) && (dy_old ==0) %starting point of path, arbitrarily start at top + if (dx_old == 0) && (dy_old ==0) %starting point of path, arbitrarily start at top middle spt = [floor(n/2),1]; elseif (dx_old > 0) && (dy_old == 0) %coming from the left spt = [1, floor(m/2)]; elseif (dx_old < 0) && (dy_old == 0) %coming from the right spt = [n,floor(m/2)]; - elseif (dx_old > 0) && (dy_old > 0) %coming from the bottom left + elseif (dx_old > 0) && (dy_old < 0) %coming from the bottom left spt = [1,m]; - elseif (dx_old > 0) && (dy_old < 0) %coming from the top left + elseif (dx_old > 0) && (dy_old > 0) %coming from the top left spt = [1,1]; - elseif (dx_old < 0) && (dy_old > 0) %coming from the bottom right + elseif (dx_old < 0) && (dy_old < 0) %coming from the bottom right spt = [n,m]; - elseif (dx_old < 0) && (dy_old < 0) %coming from the top right + elseif (dx_old < 0) && (dy_old > 0) %coming from the top right spt = [n,1]; - elseif (dx_old == 0) && (dy_old > 0) %coming from the bottom + elseif (dx_old == 0) && (dy_old < 0) %coming from the bottom spt = [floor(n/2),m]; - elseif (dx_old == 0) && (dy_old < 0) %coming from the top + elseif (dx_old == 0) && (dy_old > 0) %coming from the top spt = [floor(n/2),1]; end %ept if statements - depend on new values From 6003f5a4ca8a18ef660c69f17db2adca53898279 Mon Sep 17 00:00:00 2001 From: ezcawi Date: Tue, 10 Feb 2015 10:11:05 -0500 Subject: [PATCH 39/50] Computes Overlook Probability Function Computes overlook probability function presented in kadane 2014 for (L,J) pairs. Simplifying assumptions: No partial searches: The quad either takes a pic of the entire area or doesn't One search technology: Multirotor Overlook probability and cost of search are the same no matter how many times an area is searched (the same in j) (L(i),J(i)) represents searching location L(i) J(i) times --- senior-design-code/S.m | 25 +++++++++++++++++++++++++ 1 file changed, 25 insertions(+) create mode 100644 senior-design-code/S.m diff --git a/senior-design-code/S.m b/senior-design-code/S.m new file mode 100644 index 0000000..9bbc940 --- /dev/null +++ b/senior-design-code/S.m @@ -0,0 +1,25 @@ +function [SLJ,C] = S(p,alpha,c,L,J) + %parameters p - probability arrays, alpha-overlook probabilities for each location, c-cost of each cell, M + %-nlj*clj in this case the values of M are equal as j increases, L- set of locations to do the sum over, J- number of times + %each location is searched. (L(i),J(i)) represents searching location + %L(i) J(i) times + %SLJ = overlook probability + %C = total cost of the path + %reshape input matrices to be column vectors, might be unnecessary but + %not taking that chance + [m,n] = size(p); + p = reshape(p,n*m,1); + alpha = reshape(alpha, n*m,1); + c = reshape(c, n*m,1); + %M = reshape(M, n*m,1); %don't have partial searches so commenting this out + %take the parts that we need + pl = p(L); + alphal = alpha(L); + cl = c(L); + %Ml = M(L); + %compute S((L,J)) assuming that the costs and overlook probabilities are + %the same for every j and only probabilities are different + Il = log(alphal); %since not worrying about M, don't need to divide by c + eIlj = exp(J.*Il);%equivalent to e^(sum_j(Ml,j Il,j)) since the values are equal in J and M = C + SLJ = sum(pl.* eIlj); + C = sum(J.*cl); %total cost incurred by the path assuming it costs the same to search each cell From 486291cb05e040d24c1abb496a9889e7405405ed Mon Sep 17 00:00:00 2001 From: ezcawi Date: Fri, 27 Feb 2015 11:44:06 -0500 Subject: [PATCH 40/50] adding optimized area cost --- senior-design-code/big_path | 159 ++++++++++++++++++++++++------------ 1 file changed, 107 insertions(+), 52 deletions(-) diff --git a/senior-design-code/big_path b/senior-design-code/big_path index 95cca4f..becc2c7 100644 --- a/senior-design-code/big_path +++ b/senior-design-code/big_path @@ -6,73 +6,128 @@ %inputs from user: probability files, number of levels to go down - N % the inputs should be interfaced with the gui somehow clear all; close all; clc; +%declare some global variables +global pl alphal m_current n_current m_new n_new Cp N = input('How many layers of Resolution?: '); %integer prob_detect = input('What is the UAVs probability of detection?: '); % number between 0 and 1 -top_level_struct = load(input('Top Level File Name: ','s')); %filename -fname = fieldnames(top_level_struct); -top_level_probs = top_level_struct.(fname{1}); % have to do this to convert from struct to array -C = input('What percentage of the total area can the UAV cover?: '); % needs a number between 0 and 1 -starting_path = input('What is the starting path through the top level?: '); % needs an n x 2 array -%need to check inputs for final version -current_path = starting_path; -lower_res_probs = top_level_probs; -[m_old, n_old] = size(top_level_probs); -for i = 2:N %for every level of resolution +Cp = input('What percentage of the total area can the UAV cover?: '); % needs a number between 0 and 1 +output_filename = 'C:\Users\Eric Cawi\Documents\senior design\test.txt'; +fid = fopen(output_filename,'w'); +%write some basic information to the test file +strN = sprintf('%d layers of resolution',N); +fprintf(fid,'%s \r\n',strN); +strp = sprintf('probability of detection: %f',prob_detect); +fprintf(fid,'%s \r\n',strp); +strc = sprintf('Total Area that the system can cover: %f',Cp); +fprintf(fid,'%s \r\n',strc); +for i = 1:N %for every level of resolution str = sprintf('Resolution level %d ',i); disp(str); - current_res_struct = load(input('next resolution probability array: ','s')); + current_res_struct = load(input('resolution probability array: ','s')); fname = fieldnames(current_res_struct); current_res_probs = current_res_struct.(fname{1}); - [m_new, n_new] = size(current_res_probs); - path_part = zeros(m_new*n_new, 2,length(current_path(:,1))); - for j = 1: length(current_path(:,1)) - %figure out which segment of the path to use for path planning, - %assumes that m_new is a multiple of m_old, i.e. each cell in the - %old resolution is expanding into an integer number of cells in the - %new resolution - m_current = m_new/m_old; - n_current = n_new/n_old; - current_probs = current_res_probs((current_path(j,2)-1)*m_current+1:current_path(j,2)*m_current, (current_path(j,1)-1)*n_current+1:current_path(j,1)*n_current); - - C_area = floor(C*lower_res_probs(current_path(j,2), current_path(j,1))*numel(current_res_probs)); %current cost * probability subject is in an area = proportional time spent in that area * number of cells = number of cells in new resolution we can cover - %calculate starting and ending points for this level from previous - %path resolution - if j=2 - dx_old = current_path(j,1) - current_path(j-1,1);%previous direction, determines starting point - dy_old = current_path(j,2) - current_path(j-1,2); - else %j = 1 - dx_old = 0; - dy_old = 0; + lower_res_probs = current_res_probs; + else + [m_new, n_new] = size(current_res_probs); + path_part = zeros(m_new*n_new, 2,length(current_path(:,1))); + %prepare to optimize effort + pl = zeros(length(current_path(:,1)),1); + for k =1:length(current_path(:,1)) + pl(k) = lower_res_probs(current_path(k,2),current_path(k,1)); end - [spt,ept] = endpts(dx_old, dy_old, dx_new, dy_new,m_current, n_current); - if ((ept(1) == 0)&&(ept(2) == 0)) %at the end of the path, do a greedy algorithm - disp('end of current path') - path_part(1:C_area,:,j) = greedy_algorithm(current_probs, prob_detect, spt, C_area); - elseif ((dx_old == 0)&&(dy_old == 0)) %this happens at the beginning, do greedy from the endpoint and then reverse the path - disp('start of current path'); - path_temp = greedy_algorithm(current_probs,prob_detect,ept, C_area); % ending point is starting point of the path - path_part(1:C_area,:,j) = flipud(path_temp);%flips up/down so that the ending point is hte last part of the array - else % the normal case, hill climbing with fixed start and end points - disp('hill climbing'); - path_part(1:C_area,:,j) = hill_climbing(m_current,n_current,current_probs,prob_detect,C_area,spt,ept); + alphal = (1-prob_detect)*ones(size(pl)); + m_current = m_new/m_old; + n_current = n_new/n_old; + %determine the effort for the points in the current path with fmincon + n0 = rand(length(current_path(:,1)),1); + n = fmincon(@S,n0,[],[],[],[],zeros(size(n0)),ones(size(n0)),@cnstr); + disp(n) + for j = 1: length(current_path(:,1)) + %figure out which segment of the path to use for path planning, + %assumes that m_new is a multiple of m_old, i.e. each cell in the + %old resolution is expanding into an integer number of cells in the + %new resolution + current_probs = current_res_probs((current_path(j,2)-1)*m_current+1:current_path(j,2)*m_current, (current_path(j,1)-1)*n_current+1:current_path(j,1)*n_current); + C_area = floor(m_current*n_current*n(j)); %current cost * probability subject is in an area = proportional time spent in that area * number of cells = number of cells in new resolution we can cover + disp(C_area) + %calculate starting and ending points for this level from previous + %path resolution + if j=2 + dx_old = current_path(j,1) - current_path(j-1,1);%previous direction, determines starting point + dy_old = current_path(j,2) - current_path(j-1,2); + else %j = 1 + dx_old = 0; + dy_old = 0; + end + [spt,ept] = endpts(dx_old, dy_old, dx_new, dy_new,m_current, n_current); + if ((ept(1) == 0)&&(ept(2) == 0)) %at the end of the path, do a greedy algorithm + disp('end of current path') + path_part(1:C_area,:,j) = greedy_algorithm(current_probs, prob_detect, spt, C_area); + elseif ((dx_old == 0)&&(dy_old == 0)) %this happens at the beginning, do greedy from the endpoint and then reverse the path + disp('start of current path'); + path_temp = greedy_algorithm(current_probs,prob_detect,ept, C_area); % ending point is starting point of the path + path_part(1:C_area,:,j) = flipud(path_temp);%flips up/down so that the ending point is hte last part of the array + else % the normal case, hill climbing with fixed start and end points + disp('hill climbing'); + path_part(1:C_area,:,j) = hill_climbing(m_current,n_current,current_probs,prob_detect,C_area,spt,ept); + end + strn = sprintf('x y - level %d path through level %d point (%d,%d)',i,i-1, current_path(j,1)-1,current_path(j,2)-1); + fprintf(fid,'\r\n%s \r\n',strn); + for k = 1:length(path_part(1:C_area,1,j)) + fprintf(fid,'%d %d \r\n',path_part(k,:,j)-1); + end end + %link the paths together to update current path + current_path = link_paths(path_part,current_path, m_current,n_current,m_new, n_new); + lower_res_probs = current_res_probs; % update the lower resolution for the next iteration + m_old = m_new; + n_old = n_new; end - %link the paths together to update current path - current_path = link_paths(path_part,current_path, m_current,n_current); - lower_res_probs = current_res_probs; % update the lower resolution for the next iteration end -disp('done'); +%close the file at the end %compute the final path sum path_sum = 0; +disp('plotting') +background_probs = imread('C:\Users\Eric Cawi\Documents\senior design\ring_test.png'); +min_x = 1; +max_x = 40; +min_y = 1; +max_y = 40; +figure(); +imagesc([min_x max_x], [min_y max_y], background_probs); +hold on +plot(current_path(:,1),current_path(:,2),'-*','linewidth',2.25) +plot(current_path(1,1),current_path(1,2),'r*','markersize',20) +quiver(current_path(:,1),current_path(:,2),[current_path(2:end,1)-current_path(1:end-1,1);0],[current_path(2:end,2)-current_path(1:end-1,2);0],'linewidth',2.5) +xlabel('x','fontsize',14) +ylabel('y','fontsize',14) +title('Search Path overlaid onto the probability matrix','fontsize',20) for j = 1:length(current_path(:,1)) path_sum = path_sum + prob_detect*current_res_probs(current_path(j,2), current_path(j,1)); current_res_probs(current_path(j,2), current_path(j,1)) = current_res_probs(current_path(j,2), current_path(j,1))*(1-prob_detect); end +strsum = sprintf('total probability covered: %f', path_sum); +fprintf(fid,'\r\n%s \r\n',strsum); +fclose(fid); +disp('done'); +%note, for the output file, the indices are presented as python compatible, +%for further matlab stuff you'll have to add 1 to every array From 61bef02aa6fe1846d14898aab39cfd8a5d1935c4 Mon Sep 17 00:00:00 2001 From: ezcawi Date: Fri, 27 Feb 2015 11:45:23 -0500 Subject: [PATCH 41/50] Update S.m --- senior-design-code/S.m | 37 ++++++++++++------------------------- 1 file changed, 12 insertions(+), 25 deletions(-) diff --git a/senior-design-code/S.m b/senior-design-code/S.m index 9bbc940..21ef94e 100644 --- a/senior-design-code/S.m +++ b/senior-design-code/S.m @@ -1,25 +1,12 @@ -function [SLJ,C] = S(p,alpha,c,L,J) - %parameters p - probability arrays, alpha-overlook probabilities for each location, c-cost of each cell, M - %-nlj*clj in this case the values of M are equal as j increases, L- set of locations to do the sum over, J- number of times - %each location is searched. (L(i),J(i)) represents searching location - %L(i) J(i) times - %SLJ = overlook probability - %C = total cost of the path - %reshape input matrices to be column vectors, might be unnecessary but - %not taking that chance - [m,n] = size(p); - p = reshape(p,n*m,1); - alpha = reshape(alpha, n*m,1); - c = reshape(c, n*m,1); - %M = reshape(M, n*m,1); %don't have partial searches so commenting this out - %take the parts that we need - pl = p(L); - alphal = alpha(L); - cl = c(L); - %Ml = M(L); - %compute S((L,J)) assuming that the costs and overlook probabilities are - %the same for every j and only probabilities are different - Il = log(alphal); %since not worrying about M, don't need to divide by c - eIlj = exp(J.*Il);%equivalent to e^(sum_j(Ml,j Il,j)) since the values are equal in J and M = C - SLJ = sum(pl.* eIlj); - C = sum(J.*cl); %total cost incurred by the path assuming it costs the same to search each cell +function [fx] = S(x) + %f(x) for kadane's probability function + %computes suml(pl*exp(sumj(xlj*log(alphalj)) + %x is the matrix of continuous variables given by [n1,n2,...nL] is the + %effort put into each search area, since with simplifications j isn't + %needed + %pl, alphal are taken as global varables + global alphal pl + eIl = exp(x.*log(alphal));% e^(sumj(n*log(alphal + fx = sum(pl.* eIl); + + From 33839fa630a96c378e143121444f3eddc4e6454e Mon Sep 17 00:00:00 2001 From: ezcawi Date: Fri, 27 Feb 2015 11:46:24 -0500 Subject: [PATCH 42/50] Constraint function for big path --- senior-design-code/cnstr.m | 9 +++++++++ 1 file changed, 9 insertions(+) create mode 100644 senior-design-code/cnstr.m diff --git a/senior-design-code/cnstr.m b/senior-design-code/cnstr.m new file mode 100644 index 0000000..3af2904 --- /dev/null +++ b/senior-design-code/cnstr.m @@ -0,0 +1,9 @@ +%cost constraint function +function [c,ceq] = cnstr(x) + global m_new n_new m_current n_current Cp + %x is the variable representing effort of searching a search region + %m_new/current, n_new/current are global variables from the big_path.m + %code + C = floor(Cp*m_new*n_new); %x should be m_new*n_new long + c =sum(x*(m_current*n_current))-C; %assumes search areas have equal probability + ceq = [];%not using this constraint From 6b7fd0ed424338ae4075e0590a05c76847ea3908 Mon Sep 17 00:00:00 2001 From: ezcawi Date: Tue, 13 Oct 2015 14:28:03 -0500 Subject: [PATCH 43/50] simulates the movement of gas through a non-uniform field --- arc-models/gas_move.m | 38 ++++++++++++++++++++++++++++++++++++++ 1 file changed, 38 insertions(+) create mode 100644 arc-models/gas_move.m diff --git a/arc-models/gas_move.m b/arc-models/gas_move.m new file mode 100644 index 0000000..fb413f2 --- /dev/null +++ b/arc-models/gas_move.m @@ -0,0 +1,38 @@ +function gas = gas_move(con,speed_mat,src_vec, sink_vec,num_it) +%con- initial concentrations of the gas +%speed_mat- matrix of the max percentage transference,takes values between +%0 and .25-each cell can get rid of at most 25% of it's gas +%src_vec = locations of sources in the thing +%sink_vec= locations of sinks in the domain +%num_it- number of iterations to transfer +[m,n] = size(con);%also the size of speed matrix +if ~isempty(src_vec) + src_inds = sub2ind([m,n],src_vec(:,1),src_vec(:,2)); +end +if ~isempty(sink_vec) + sink_inds = sub2ind([m,n],sink_vec(:,1),sink_vec(:,2)); +end +for i = 1:num_it + %create left, right, up, and down transition matrices, keeping in mind + %the borders only have things from 2-3 directions + left = speed_mat(1:m,1:n-1).*(con(1:m,1:n-1)-con(1:m,2:n));%represents the contribution of the cells coming from the left + right = speed_mat(1:m,2:n).*(con(1:m,2:n)-con(1:m,1:n-1));%contribution of cells coming from the right + up = speed_mat(1:m-1,1:n).*(con(1:m-1,1:n)-con(2:m,1:n));%contribution from the cells above + down = speed_mat(2:m,1:n).*(con(2:m,1:n)-con(1:m-1,1:n));%contribution from the below cells + %update concentration + con(1:m,2:n) = con(1:m,2:n)+left; + con(1:m,1:n-1) = con(1:m,1:n-1)+right; + con(2:m,1:n) = con(2:m,1:n)+up; + con(1:m-1,1:n)=con(1:m-1,1:n)+down; + if ~isempty(src_vec) + con(src_inds)=1;%fills up the source(s) + end + if ~isempty(sink_vec) + con(sink_inds)=0;%drains sink + end +end +%at the end number of desired iterations, gas is the final concentration +gas = con; + + + From 5d3c980b49cdf14c433b2ca1b59d8f51a41031ee Mon Sep 17 00:00:00 2001 From: ezcawi Date: Tue, 13 Oct 2015 14:29:30 -0500 Subject: [PATCH 44/50] uses gas_move.m to simulate hikers on real terrain data --- arc-models/walking_speed.m | 33 +++++++++++++++++++++++++++++++++ 1 file changed, 33 insertions(+) create mode 100644 arc-models/walking_speed.m diff --git a/arc-models/walking_speed.m b/arc-models/walking_speed.m new file mode 100644 index 0000000..fc40a68 --- /dev/null +++ b/arc-models/walking_speed.m @@ -0,0 +1,33 @@ +function map = walking_speed(elev,R,num_it,fname) +%creates a map simulating hiker movement as gas +%elev-elevation +%R - reference matrix in terms of [cells/degree,top latitude, west +%longitude +%num_it - number of iterations +%fname - file to which final image is saved + +%first, calculate the slope function +[ASPECT, SLOPE, gradN, gradE] = gradientm(elev, R); +%calculate walking speed using tobler's formula +speed_mat = 6*exp(-3.5*abs(tand(SLOPE)+.05)); +%scale speed matrix down to [0,.25] +speed_mat = .25*speed_mat/max(max(speed_mat)); + +%downsample things to make managable calculations +speed_mat = downsample(downsample(speed_mat,3)',3)'; +%set up last known point +[m,n] = size(speed_mat); +con = zeros(m,n); +con(ceil(m/2),ceil(n/2)) = 1; +%run the iterations +con_2 = gas_move(con,speed_mat,[ceil(m/2) ceil(n/2)],[],num_it); +%upsample +con_2 = con_2/sum(sum(con_2));%normalizing +map = con_2; +figure(); imagesc(map); +%create and resize png to 5000 by 5000 +im_map = imresize(map,[5000,5000]); +%write image to file +im_map = mat2gray(im_map); +imwrite(im_map,fname); +end From 43e90c74f5117ed54dee0fe7870a9dfa3512bbb7 Mon Sep 17 00:00:00 2001 From: ezcawi Date: Tue, 20 Oct 2015 14:01:34 -0500 Subject: [PATCH 45/50] writeup of progress on diffusion model. --- arc-models/diffusionwriteup.tex | 48 +++++++++++++++++++++++++++++++++ 1 file changed, 48 insertions(+) create mode 100644 arc-models/diffusionwriteup.tex diff --git a/arc-models/diffusionwriteup.tex b/arc-models/diffusionwriteup.tex new file mode 100644 index 0000000..6f0d767 --- /dev/null +++ b/arc-models/diffusionwriteup.tex @@ -0,0 +1,48 @@ +\documentclass[]{article} + +%opening +\title{Modeling Lost Person Movement as Diffusion based on elevation and land cover} +\author{Eric Cawi} + +\begin{document} + +\maketitle + +\begin{abstract} +This work creates a model for lost person movement based on gas diffusion. to be filled in +\end{abstract} +\tableofcontents + +\section{Create Working Gas Model} +Burgess and Darken created a simple gas diffusion model to create avenues of approach. \textbf{Insert diagram here}. First, the area is divided into cells. Then the elevation data is transformed onto the interval $[0,.25]$, which represents the maximum transfer of gas from that cell to the another cell. For example, a cell with a value of .25 and a concentration of 1 (units) can transfer up to .25 of it's gas to each of the four neighboring cells. This ensures the conservation of concentration in the simulation. This simulation also allows for sources and sinks to be placed at arbitrary grid cells. Burgess and Darken run this simulation until a relative equilibrium is reached (mathematically taking the limit as time goes to infinity), and then use the gradient of the final concentration to generate possible paths of approach for military applications. Note-using the gradient to generate the avenues of approach is solving the transport equation-hard to solve numerically. + +For search and rescue, I am more interested in the final gas concentration. Modeling the last known point as a source, the gas model simulates the overall movement of many hikers from the last known point (assertion, want to prove this). I use the source because it models a spike of probability at the LKP, and as B/D used originally it represents a continuous stream of people moving around the material + +In this model, I use elevation/land cover data as the speed at which the gas is allowed to propagate. Because the arrays involved are very large, I downsample by a factor of 4 on the NED dataset, then resample up to 5001x5001 pixels using Matlab's default interpolation- i think it's bicubic. Currently I'm running out to a huge number of iterations because I haven't calibrated the transformations to match different categories/mean distances yet. +\section{Transformations from elevation/land cover to impedance} +Elevation- using Tobler's hiking function +Land Cover - using 1-Don's values +\section{Diffusion Equation with source} +The diffusion equation models the overall random walks of particles diffusing through materials. At each point there is some diffusion constant determining the maximum velocity of the concentration. Mathematically, freely diffusing gas like the Burgess Darken model can be modeled with this equation. In 1-D with a constant diffusion coefficient, the solution is a normal distribution, and generalizations give similarly smooth distributions. In our case the diffusion is given by: +\[\frac{dc(x,y,t)}{dt}=\nabla \cdot(D(x,y)\nabla c(x,y,t)=\] +\[(\frac{dD(x,y)}{dx},\frac{dD(x,y)}{dy})\cdot\nabla c(x,y,t) + D(x,y)\Delta c(x,y,t)\] +Current work: finite difference scheme for this equation. +\section{Boundary Considerations} +I'm considering using the boundary of the map as a sink and seeing what happens, or having some other way of modeling the interactions with the outside layer because in the current state we might be getting innaccurate reflection. Boundary conditions will be important for the PDE formulation as well. + +\section{Static Model Preliminary Testing} +Done on 10 arizona cases I had elevation data for. Observations- On flat terrain the concentration behaves like regular diffusion, i.e. a smooth version of the distance rings. Otherwise, it avoids big changes in elevation, as expected. E.g. if a subject doesn't start on a mountain then the mountain will have lower probability. However, a 2 or 3 of the hikers tested were going for the mountain, which makes sense because they were probably trying to hike on the mountain, but this model isn't going to catch that. Fixing this would require more information about the subject, for example "They were heading towards Mt. Awesome when they got lost." Since the odds are they moved in the direction of Mt. Awesome, a sink dropped at the intended location would skew the concentration in that direction. But that isn't available all the time. +Not done for land cover yet +\section{Scaling with Bob's Book} +not done +\section{Testing/fitting on Yosemite Cases} +Prediction- lots of mountains in yosemite, so we should see higher probabilities around ridges/trails that stay at the same elevation. I would like to somehow set up this as a batch + +\section{Comparison with diffusion Models} +Compare results and speed. +\section{Discussion} +Goal: An easy GIS ready application/A diffusion model that responds to terrain and is about as good as the distance model, which means that the search team now has more information. + + +%Need Citations +\end{document} From c176ec46501e56b91c92bd25f57e4a5f8010158f Mon Sep 17 00:00:00 2001 From: ezcawi Date: Tue, 20 Oct 2015 14:02:55 -0500 Subject: [PATCH 46/50] generalized to land cover as well --- arc-models/walking_speed.m | 99 ++++++++++++++++++++++++++++---------- 1 file changed, 74 insertions(+), 25 deletions(-) diff --git a/arc-models/walking_speed.m b/arc-models/walking_speed.m index fc40a68..e52f6b1 100644 --- a/arc-models/walking_speed.m +++ b/arc-models/walking_speed.m @@ -1,33 +1,82 @@ -function map = walking_speed(elev,R,num_it,fname) +function map = walking_speed(data1,data2,R1,R2,num_it,elev_flag,elev_weight,land_flag,land_weight,fname) %creates a map simulating hiker movement as gas -%elev-elevation +%data1 - elevation data +%data2 - land cover data %R - reference matrix in terms of [cells/degree,top latitude, west -%longitude +%longitude] -need one reference matrix for both sets %num_it - number of iterations -%fname - file to which final image is saved - +%elev_flag = doing based on elevation +%land_flag - doing calculation based on land cover +%elev/land weight - weight for combination, if applicable %first, calculate the slope function -[ASPECT, SLOPE, gradN, gradE] = gradientm(elev, R); +if elev_flag==1 && land_flag==0; + figure(); imagesc(data1); title('elevation data') + [ASPECT, SLOPE, gradN, gradE] = gradientm(data1, R1); %calculate walking speed using tobler's formula -speed_mat = 6*exp(-3.5*abs(tand(SLOPE)+.05)); + speed_mat1 = 6*exp(-3.5*abs(tand(SLOPE)+.05)); %scale speed matrix down to [0,.25] -speed_mat = .25*speed_mat/max(max(speed_mat)); + speed_mat1 = .25*speed_mat1/max(max(speed_mat1)); + %downsample + speed_mat1 = downsample(downsample(speed_mat1,12)',12)'; +elseif land_flag ==1&&elev_flag==0; + %convert land cover to impedance + speed_mat1 = land2speed(data2);%using speed mat 1 in case of only using one dataset + figure(); imagesc(speed_mat1);title('impedance data') + %scale down to [0,.25] + speed_mat1=.25*speed_mat1/max(max(speed_mat1)); + %downsample + speed_mat1 = downsample(downsample(speed_mat1,4)',4)'; +elseif land_flag ==1 && elev_flag==1; + %use two speed matrices + figure(); imagesc(data1);title('elevation data') + [ASPECT, SLOPE, gradN, gradE] = gradientm(data1, R1); +%calculate walking speed using tobler's formula + speed_mat1 = 6*exp(-3.5*abs(tand(SLOPE)+.05)); + %scale speed matrix down to [0,.25] + speed_mat1 = .25*speed_mat1/max(max(speed_mat1)); + %calculate second speed matrix, assumed to be land cover + speed_mat2 = land2speed(data2);%using speed mat 1 in case of only using one dataset + figure(); imagesc(speed_mat2);title('impedance data') + %scale down to [0,.25] + speed_mat2=.25*speed_mat2/max(max(speed_mat2)); + %downsample + speed_mat1 = downsample(downsample(speed_mat1,4)',4)'; + speed_mat2 = downsample(downsample(speed_mat2,4)',4)'; +end -%downsample things to make managable calculations -speed_mat = downsample(downsample(speed_mat,3)',3)'; -%set up last known point -[m,n] = size(speed_mat); -con = zeros(m,n); -con(ceil(m/2),ceil(n/2)) = 1; -%run the iterations -con_2 = gas_move(con,speed_mat,[ceil(m/2) ceil(n/2)],[],num_it); -%upsample -con_2 = con_2/sum(sum(con_2));%normalizing -map = con_2; -figure(); imagesc(map); -%create and resize png to 5000 by 5000 -im_map = imresize(map,[5000,5000]); -%write image to file -im_map = mat2gray(im_map); -imwrite(im_map,fname); + +%if we're only doing elevation or land: +if (elev_flag==0)||(land_flag ==0); + %set up last known point + [m,n] = size(speed_mat1); + con = zeros(m,n); + con(ceil(m/2),ceil(n/2)) = 1; + %run the iterations + con_2 = gas_move(con,speed_mat1,[ceil(m/2) ceil(n/2)],[],num_it); + %upsample + con_2 = con_2/sum(sum(con_2));%normalizing + map = con_2; + figure(); imagesc(map); + %create and resize png to 5000 by 5000 + im_map = imresize(map,[5001,5001]); + %write image to file + im_map = mat2gray(im_map); + imwrite(im_map,fname); +else %using both + [m1,n1] = size(speed_mat1); + [m2,n2]= size(speed_mat2); + con1 = zeros(m1,n1); + con1(ceil(m1/2),ceil(n1/2)) = 1; + con2 = zeros(m2,n2); + con2(ceil(m2/2),ceil(n2/2)) = 1; + con_1 = gas_move(con1,speed_mat1,[ceil(m1/2) ceil(n1/2)],[],num_it); + con_2 = gas_move(con2,speed_mat2,[ceil(m2/2) ceil(n2/2)],[],num_it); + con_1 = con_1/sum(sum(con_1)); + con_2 = con_2/sum(sum(con_2)); + im_map1 =imresize(con_1,[5001,5001]); + im_map2 = imresize(con_2,[5001,5001]); + map = elev_weight*im_map1+land_weight*im_map2;%weighted average of the two + im_map = mat2gray(map); + imwrite(im_map,fname); +end end From 0a3ba2fecbe5c59f953d941c6e7f20a79a6570ba Mon Sep 17 00:00:00 2001 From: ezcawi Date: Tue, 27 Oct 2015 12:01:05 -0500 Subject: [PATCH 47/50] cleaned up if statements --- arc-models/walking_speed.m | 71 ++++++++++++++------------------------ 1 file changed, 25 insertions(+), 46 deletions(-) diff --git a/arc-models/walking_speed.m b/arc-models/walking_speed.m index e52f6b1..bae69e0 100644 --- a/arc-models/walking_speed.m +++ b/arc-models/walking_speed.m @@ -1,9 +1,9 @@ -function map = walking_speed(data1,data2,R1,R2,num_it,elev_flag,elev_weight,land_flag,land_weight,fname) +function map = walking_speed(data1,data2,R1,num_it,elev_flag,elev_weight,land_flag,land_weight,fname) %creates a map simulating hiker movement as gas %data1 - elevation data %data2 - land cover data %R - reference matrix in terms of [cells/degree,top latitude, west -%longitude] -need one reference matrix for both sets +%longitude] or georeference object %num_it - number of iterations %elev_flag = doing based on elevation %land_flag - doing calculation based on land cover @@ -17,15 +17,12 @@ %scale speed matrix down to [0,.25] speed_mat1 = .25*speed_mat1/max(max(speed_mat1)); %downsample - speed_mat1 = downsample(downsample(speed_mat1,12)',12)'; + speed_mat1 = downsample(downsample(speed_mat1,4)',4)'; elseif land_flag ==1&&elev_flag==0; %convert land cover to impedance speed_mat1 = land2speed(data2);%using speed mat 1 in case of only using one dataset figure(); imagesc(speed_mat1);title('impedance data') - %scale down to [0,.25] - speed_mat1=.25*speed_mat1/max(max(speed_mat1)); - %downsample - speed_mat1 = downsample(downsample(speed_mat1,4)',4)'; + elseif land_flag ==1 && elev_flag==1; %use two speed matrices figure(); imagesc(data1);title('elevation data') @@ -35,48 +32,30 @@ %scale speed matrix down to [0,.25] speed_mat1 = .25*speed_mat1/max(max(speed_mat1)); %calculate second speed matrix, assumed to be land cover - speed_mat2 = land2speed(data2);%using speed mat 1 in case of only using one dataset + speed_ma%downsample + speed_mat1 = downsample(downsample(speed_mat1,4)',4)';t2 = land2speed(data2);%using speed mat 1 in case of only using one dataset figure(); imagesc(speed_mat2);title('impedance data') - %scale down to [0,.25] - speed_mat2=.25*speed_mat2/max(max(speed_mat2)); + %resize one matrix to the size of the other + speed_mat2 = imresize(speed_mat2,size(speed_mat1)); + %combine speed matrix + speed_mat1 = elev_weight*speed_mat1 + land_weight*speed_mat2; %downsample speed_mat1 = downsample(downsample(speed_mat1,4)',4)'; - speed_mat2 = downsample(downsample(speed_mat2,4)',4)'; end +%set up last known point +[m,n] = size(speed_mat1); +con = zeros(m,n); +con(ceil(m/2),ceil(n/2)) = 1; +%run the iterations +con_2 = gas_move(con,speed_mat1,[ceil(m/2) ceil(n/2)],[],num_it); +%upsample +con_2 = con_2/sum(sum(con_2));%normalizing +map = con_2; +figure(); imagesc(map); +%create and resize png to 5000 by 5000 +im_map = imresize(map,[5001,5001]); +%write image to file +im_map = mat2gray(im_map); +imwrite(im_map,fname); - -%if we're only doing elevation or land: -if (elev_flag==0)||(land_flag ==0); - %set up last known point - [m,n] = size(speed_mat1); - con = zeros(m,n); - con(ceil(m/2),ceil(n/2)) = 1; - %run the iterations - con_2 = gas_move(con,speed_mat1,[ceil(m/2) ceil(n/2)],[],num_it); - %upsample - con_2 = con_2/sum(sum(con_2));%normalizing - map = con_2; - figure(); imagesc(map); - %create and resize png to 5000 by 5000 - im_map = imresize(map,[5001,5001]); - %write image to file - im_map = mat2gray(im_map); - imwrite(im_map,fname); -else %using both - [m1,n1] = size(speed_mat1); - [m2,n2]= size(speed_mat2); - con1 = zeros(m1,n1); - con1(ceil(m1/2),ceil(n1/2)) = 1; - con2 = zeros(m2,n2); - con2(ceil(m2/2),ceil(n2/2)) = 1; - con_1 = gas_move(con1,speed_mat1,[ceil(m1/2) ceil(n1/2)],[],num_it); - con_2 = gas_move(con2,speed_mat2,[ceil(m2/2) ceil(n2/2)],[],num_it); - con_1 = con_1/sum(sum(con_1)); - con_2 = con_2/sum(sum(con_2)); - im_map1 =imresize(con_1,[5001,5001]); - im_map2 = imresize(con_2,[5001,5001]); - map = elev_weight*im_map1+land_weight*im_map2;%weighted average of the two - im_map = mat2gray(map); - imwrite(im_map,fname); -end end From 3ec649f40537194da169e2e279feb83ee5e02524 Mon Sep 17 00:00:00 2001 From: ezcawi Date: Tue, 27 Oct 2015 12:02:40 -0500 Subject: [PATCH 48/50] converts land cover to impedance matrix --- arc-models/land2speed.m | 13 +++++++++++++ 1 file changed, 13 insertions(+) create mode 100644 arc-models/land2speed.m diff --git a/arc-models/land2speed.m b/arc-models/land2speed.m new file mode 100644 index 0000000..849512c --- /dev/null +++ b/arc-models/land2speed.m @@ -0,0 +1,13 @@ +function speed_mat = land2speed(data,land_cover_class) +%converts the land cover data to impedance values +%land cover class has two columns, 1 is the values taken by the raster. The second is impedance values. In my case i'm scaling down to [0,.25] for the gas model +temp = data; +for i = 1:length(land_cover_class(:,1)) + inds = find(data == land_cover_class(i,1));% find all of the points that have a particular classification + temp(inds) = land_cover_class(i,2); %reset all of the data +end +%get rid of any border values by setting unknown things to 0 go +inds = find(temp>.25); +temp(inds) = 0; +speed_mat = temp; +end From 69651f3ebedcbcf23009f434be9e737358df97da Mon Sep 17 00:00:00 2001 From: ezcawi Date: Tue, 27 Oct 2015 12:23:53 -0500 Subject: [PATCH 49/50] better if statements, downsample at end of if block --- arc-models/walking_speed.m | 61 ++++++++++++++++++++++++++++++++++++++ 1 file changed, 61 insertions(+) diff --git a/arc-models/walking_speed.m b/arc-models/walking_speed.m index bae69e0..d27bbad 100644 --- a/arc-models/walking_speed.m +++ b/arc-models/walking_speed.m @@ -7,7 +7,68 @@ %num_it - number of iterations %elev_flag = doing based on elevation %land_flag - doing calculation based on land cover +%elev/land weight - weight for combination, if applicabfunction map = walking_speed(data1,data2,R1,num_it,elev_flag,elev_weight,land_flag,land_weight,landclass,fname) +%creates a map simulating hiker movement as gas +%data1 - elevation data +%data2 - land cover data +%R - reference matrix in terms of [cells/degree,top latitude, west +%longitude] or georeference object +%num_it - number of iterations +%elev_flag = doing based on elevation +%land_flag - doing calculation based on land cover %elev/land weight - weight for combination, if applicable +%fname- desired file location name +%first, calculate the slope function +if elev_flag==1 && land_flag==0; + figure(); imagesc(data1); title('elevation data') + [ASPECT, SLOPE, gradN, gradE] = gradientm(data1, R1); +%calculate walking speed using tobler's formula + speed_mat1 = 6*exp(-3.5*abs(tand(SLOPE)+.05)); +%scale speed matrix down to [0,.25] + speed_mat1 = .25*speed_mat1/max(max(speed_mat1)); +elseif land_flag ==1&&elev_flag==0; + %convert land cover to impedance + speed_mat1 = land2speed(data2,landclass);%using speed mat 1 in case of only using one dataset + figure(); imagesc(speed_mat1);title('impedance data') + +elseif land_flag ==1 && elev_flag==1; + %use two speed matrices + figure(); imagesc(data1);title('elevation data') + [ASPECT, SLOPE, gradN, gradE] = gradientm(data1, R1); +%calculate walking speed using tobler's formula + speed_mat1 = 6*exp(-3.5*abs(tand(SLOPE)+.05)); + %scale speed matrix down to [0,.25] + speed_mat1 = .25*speed_mat1/max(max(speed_mat1)); + %calculate second speed matrix, assumed to be land cover + figure(); imagesc(speed_mat1);title('elevation data') + + speed_mat2 = land2speed(data2,landclass);%using speed mat 1 in case of only using one dataset + figure(); imagesc(speed_mat2);title('impedance data') + %resize one matrix to the size of the other + speed_mat2 = imresize(speed_mat2,size(speed_mat1)); + %combine speed matrix + speed_mat1 = elev_weight*speed_mat1 + land_weight*speed_mat2; + +end +%downsample +speed_mat1 = downsample(downsample(speed_mat1,4)',4)'; +%set up last known point +[m,n] = size(speed_mat1); +con = zeros(m,n); +con(ceil(m/2),ceil(n/2)) = 1; +%run the iterations +con_2 = gas_move(con,speed_mat1,[ceil(m/2) ceil(n/2)],[],num_it); +%upsample +con_2 = con_2/sum(sum(con_2));%normalizing +map = con_2; +figure(); imagesc(map); +%create and resize png to 5000 by 5000 +im_map = imresize(map,[5001,5001]); +%write image to file +im_map = mat2gray(im_map); +imwrite(im_map,fname); + +endle %first, calculate the slope function if elev_flag==1 && land_flag==0; figure(); imagesc(data1); title('elevation data') From 5ba4fdb178f439dc7a28856c7e0dc85d81787f14 Mon Sep 17 00:00:00 2001 From: ctwardy Date: Wed, 23 Dec 2015 02:00:58 -0500 Subject: [PATCH 50/50] minor tweaks Minor tweaks to diffusionwriteup. Update .gitignore. --- .gitignore | 17 ++++++++++++++++- arc-models/diffusionwriteup.tex | 12 ++++++++++-- 2 files changed, 26 insertions(+), 3 deletions(-) diff --git a/.gitignore b/.gitignore index 888cefd..08e8201 100644 --- a/.gitignore +++ b/.gitignore @@ -18,4 +18,19 @@ arc-models/*_flymake.py *.7 -junk.py \ No newline at end of file +junk.py +*.rel + +*.aux + +*.bbl + +*.bib + +*.blg + +*.log + +*.toc + +arc-models/diffusionwriteup.synctex.gz diff --git a/arc-models/diffusionwriteup.tex b/arc-models/diffusionwriteup.tex index 6f0d767..4454fdc 100644 --- a/arc-models/diffusionwriteup.tex +++ b/arc-models/diffusionwriteup.tex @@ -1,7 +1,7 @@ \documentclass[]{article} %opening -\title{Modeling Lost Person Movement as Diffusion based on elevation and land cover} +\title{Modeling Lost Person Movement as Diffusion Using Elevation and Land Cover} \author{Eric Cawi} \begin{document} @@ -13,8 +13,13 @@ \end{abstract} \tableofcontents +We use Burgess and Darken's \ycite{burgess_realistic_2004} gas diffusion model of path planning to create probability maps for lost hikers. Where Burgess and Darken sample the diffusion gradient to find discrete likely paths from source to target, we use the diffusion gradient itself as a probability map, and find it is competitive with the Euclidean ring model in \cite{koester_lost_2008}. + \section{Create Working Gas Model} -Burgess and Darken created a simple gas diffusion model to create avenues of approach. \textbf{Insert diagram here}. First, the area is divided into cells. Then the elevation data is transformed onto the interval $[0,.25]$, which represents the maximum transfer of gas from that cell to the another cell. For example, a cell with a value of .25 and a concentration of 1 (units) can transfer up to .25 of it's gas to each of the four neighboring cells. This ensures the conservation of concentration in the simulation. This simulation also allows for sources and sinks to be placed at arbitrary grid cells. Burgess and Darken run this simulation until a relative equilibrium is reached (mathematically taking the limit as time goes to infinity), and then use the gradient of the final concentration to generate possible paths of approach for military applications. Note-using the gradient to generate the avenues of approach is solving the transport equation-hard to solve numerically. + +\textbf{Insert diagram here}. + +First, the area is divided into cells. Each cell has a concentrationThen the elevation data is transformed onto the interval $[0,.25]$, which represents the maximum transfer of gas from that cell to the another cell. For example, a cell with a value of .25 and a concentration of 1 (units) can transfer up to .25 of its gas to each of the four neighboring cells. This ensures the conservation of concentration in the simulation. This simulation also allows for sources and sinks to be placed at arbitrary grid cells. Burgess and Darken run this simulation until a relative equilibrium is reached (mathematically taking the limit as time goes to infinity), and then use the gradient of the final concentration to generate possible paths of approach for military applications. Note-using the gradient to generate the avenues of approach is solving the transport equation-hard to solve numerically. For search and rescue, I am more interested in the final gas concentration. Modeling the last known point as a source, the gas model simulates the overall movement of many hikers from the last known point (assertion, want to prove this). I use the source because it models a spike of probability at the LKP, and as B/D used originally it represents a continuous stream of people moving around the material @@ -43,6 +48,9 @@ \section{Comparison with diffusion Models} \section{Discussion} Goal: An easy GIS ready application/A diffusion model that responds to terrain and is about as good as the distance model, which means that the search team now has more information. +\section{References} +\bibliographystyle{plain} +\bibliography{diffusionwriteup} %Need Citations \end{document}