<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 KL-Divergence Filter
In this notebook we will test out a filter for postage stamps that uses the KL-Divergence to determine how alike the object in the postage stamp is to a model Gaussian distribution.
Create base distribution for comparison
Step1:  Run some basic tests
The function we are going to use for this filter is the relative entropy function from scipy
Step2:  Test 2
Step3:  2.2 Wider Gaussian
Here we use a Gaussian with double the sigma.
Step4:  Test 3
Step5:  As expected we see that the KL-divergence when comparing to random noise is much higher than the Gaussian distributions which are similar but still different to our test distribution.
Now that we have some evidence that this test could work, let's apply it to our false positive training set.
Apply to KBMOD False Positives
We load the new training data created by Hayden.
Step6:  Reshape all the arrays into their two-dimensional postage stamp form. Then add an offset to them so they all have a minimum value just above 0. When doing the relative entropy calculation we also normalize the stamps.
Step7:  This looks really promising! There is a clear separation between the peaks for False objects and True sources. There is an area of overlap in the results. To see how well we do when training a classifier to decide True/False when looking at a postage stamps we'll train a logistic regression next.
Optimize separation value
First break into training/test sets.
Step8:  We'll use the LogisticRegression classifier from scikit-learn to learn a boundary on the relative entropy value.
Step9:  What do the false positives look like?
Step10:  Using "real" false positives
This section uses a new set of data where the true sources are verified asteroids found in the Lori Allen data. The false sources are actual false positives from KBMOD searches.
We use the training from the synthetic training set in the first part of the notebook and apply it here to the "real" data.
<ASSISTANT_TASK:>
Python Code:
central_source = np.zeros((21, 21))
central_source[10,10] = 1.
gaussian_test = ndimage.gaussian_filter(central_source, 1.4)
imshow(gaussian_test)
colorbar()
Explanation: KL-Divergence Filter
In this notebook we will test out a filter for postage stamps that uses the KL-Divergence to determine how alike the object in the postage stamp is to a model Gaussian distribution.
Create base distribution for comparison
End of explanation
from scipy.special import rel_entr
kl_div_same_dist = np.sum(rel_entr(gaussian_test, gaussian_test))
print("The KL Divergence of the test Gaussian with itself is: {}".format(kl_div_same_dist))
Explanation: Run some basic tests
The function we are going to use for this filter is the relative entropy function from scipy: scipy.special.rel_entr
Test 1: Same distribution
The KL Divergence should be 0 when the function is compared with itself.
End of explanation
off_center_source = np.zeros((21, 21))
off_center_source[8,10] = 1.
gaussian_off_cent = ndimage.gaussian_filter(off_center_source, 1.4)
plt.imshow(gaussian_off_cent)
# We add a very small offset from 0 to avoid issues in the KL-div calculation
kl_div_off_center_dist = np.nansum(rel_entr(gaussian_test+1e-12, gaussian_off_cent+1e-12))
print("The KL Divergence of the test Gaussian with the offset Gaussian is: {}".format(kl_div_off_center_dist))
Explanation: Test 2: Slightly Different Gaussians
These should be numbers greater than 0.
2.1: Off-center Gaussian
Let's start with a Gaussian with the same width as the test Gaussian but off center by a couple pixels in one direction.
End of explanation
gaussian_doublewide = ndimage.gaussian_filter(central_source, 2.8)
plt.imshow(gaussian_doublewide)
# We add a very small offset from 0 to avoid issues in the KL-div calculation
kl_div_doublewide_dist = np.nansum(rel_entr(gaussian_test+1e-12, gaussian_doublewide+1e-12))
print("The KL Divergence of the test Gaussian with the wider Gaussian is: {}".format(kl_div_doublewide_dist))
Explanation: 2.2 Wider Gaussian
Here we use a Gaussian with double the sigma.
End of explanation
poisson_noise = np.random.RandomState(42).poisson(100, size=(21,21))
poisson_noise = poisson_noise/np.sum(poisson_noise)
plt.imshow(poisson_noise)
plt.colorbar()
# We add a very small offset from 0 to avoid issues in the KL-div calculation
kl_div_poisson_dist = np.nansum(rel_entr(gaussian_test+1e-12, poisson_noise+1e-12))
print("The KL Divergence of the test Gaussian with the Poisson noise is: {}".format(kl_div_poisson_dist))
Explanation: Test 3: Poisson Noise
Above we saw results ~1 or less when comparing different distributions to our test Gaussian. Here we will test against Poisson noise and expect to see that the results are higher than the results in test 2.
End of explanation
full_false = np.load('/epyc/users/smotherh/kbmod_epyc/notebooks/normed_individual_false.npy')
keep_idx = []
idx = 0
for f_im in full_false:
    if np.isnan(f_im[0][0]):
        idx += 1
        continue
    else:
        keep_idx.append(idx)
        idx += 1
full_false = full_false[keep_idx]
full_true = np.load('/epyc/users/smotherh/kbmod_epyc/notebooks/normed_individual_simulated.npy')
imshow(full_false[0].reshape(21,21))
imshow(full_true[0].reshape(21,21))
Explanation: As expected we see that the KL-divergence when comparing to random noise is much higher than the Gaussian distributions which are similar but still different to our test distribution.
Now that we have some evidence that this test could work, let's apply it to our false positive training set.
Apply to KBMOD False Positives
We load the new training data created by Hayden.
End of explanation
false_values = []
for f_im in full_false:
    false_stamp = f_im.reshape(21,21)
    false_stamp -= np.min(false_stamp)
    false_stamp += 1e-12
    #false_kl = special.kl_div(gaussian_test, (false_stamp/np.sum(false_stamp)))
    false_kl = special.rel_entr(gaussian_test, (false_stamp/np.sum(false_stamp)))
    false_values.append(np.nansum(false_kl))
false_values = np.array(false_values)
true_values = []
for t_im in full_true:
    true_stamp = t_im.reshape(21,21)
    true_stamp -= np.min(true_stamp)
    true_stamp += 1e-12
    #true_kl = special.kl_div(gaussian_test, (true_stamp/np.sum(true_stamp)))
    true_kl = special.rel_entr(gaussian_test, (true_stamp/np.sum(true_stamp)))
    true_values.append(np.nansum(true_kl))
true_values = np.array(true_values)
len(full_false)
len(false_values), len(true_values)
fig = plt.figure(figsize=(10,6))
n, bins, _ = hist(false_values, alpha=0.4, label='False', bins=20, range=(0, 5), density=True)
hist(true_values, alpha=0.4, label='True', bins=bins, density=True)
legend()
xlabel('Relative Entropy')
ylabel('Density')
Explanation: Reshape all the arrays into their two-dimensional postage stamp form. Then add an offset to them so they all have a minimum value just above 0. When doing the relative entropy calculation we also normalize the stamps.
End of explanation
rand_state = np.random.RandomState(24)
false_selection = rand_state.choice(np.arange(len(false_values)), size=40000, replace=False)
train_false, test_false = false_values[false_selection[:30000]], false_values[false_selection[30000:]]
true_selection = rand_state.choice(np.arange(len(true_values)), size=40000, replace=False)
train_true, test_true = true_values[true_selection[:30000]], true_values[true_selection[30000:]]
train_data = np.array(train_false)
train_labels = np.zeros(len(train_false))
train_data = np.append(train_data, train_true)
train_labels = np.append(train_labels, np.ones(len(train_true)))
train_idx = np.arange(len(train_data))
rand_state.shuffle(train_idx)
train_data = train_data[train_idx]
train_labels = train_labels[train_idx]
test_data = np.array(test_false)
test_labels = np.zeros(len(test_false))
test_data = np.append(test_data, test_true)
test_labels = np.append(test_labels, np.ones(len(test_true)))
Explanation: This looks really promising! There is a clear separation between the peaks for False objects and True sources. There is an area of overlap in the results. To see how well we do when training a classifier to decide True/False when looking at a postage stamps we'll train a logistic regression next.
Optimize separation value
First break into training/test sets.
End of explanation
from sklearn.linear_model import LogisticRegression
classifier = LogisticRegression()
classifier.fit(train_data.reshape(-1, 1), train_labels)
fig = plt.figure(figsize=(10,6))
plt.plot(np.linspace(0,5,100), classifier.predict_proba(np.linspace(0,5,100).reshape(-1,1))[:,0], lw=12)
plt.scatter(test_data, np.abs(1-test_labels), c='r', s=8, zorder=10)
plt.ylabel('False Object Probability')
plt.xlabel('KL-Divergence Value')
plt.xlim(0,5)
test_results = classifier.predict(test_data.reshape(-1,1))
test_prob = classifier.predict_proba(test_data.reshape(-1,1))
true_positive = 0
false_positive = 0
true_negative = 0
false_negative = 0
true_positive_idx = []
false_positive_idx = []
true_negative_idx = []
false_negative_idx = []
idx = 0
for test_label, test_result in zip(test_labels, test_results):
    if test_label == 1 and test_result == 1:
        true_positive += 1
        true_positive_idx.append(idx)
    elif test_label == 1 and test_result == 0:
        false_negative += 1
        false_negative_idx.append(idx)
    elif test_label == 0 and test_result == 1:
        false_positive += 1
        false_positive_idx.append(idx)
    elif test_label == 0 and test_result == 0:
        true_negative += 1
        true_negative_idx.append(idx)
    else:
        print('error')
        
    idx += 1
matrix = [[true_positive, false_positive], [false_negative, true_negative]]
true_positive
matrix
fig = plt.figure(figsize=(10,10))
plt.imshow(matrix)
plt.text(-0.0,0,'True Positive:\n {}'.format(true_positive), horizontalalignment='center', verticalalignment='center', size=18, color='K')
plt.text(1.,1.,'True Negative:\n {}'.format(true_negative), horizontalalignment='center', verticalalignment='center', size=18, color='K')
plt.text(-0.0,1,'False Negative:\n {}'.format(false_negative), horizontalalignment='center', verticalalignment='center', size=18, color='w')
plt.text(1.,0.,'False Positive:\n {}'.format(false_positive), horizontalalignment='center', verticalalignment='center', size=18, color='w')
plt.colorbar()
Explanation: We'll use the LogisticRegression classifier from scikit-learn to learn a boundary on the relative entropy value.
End of explanation
fig = plt.figure(figsize=(14,14))
for i in range(36):
    fig.add_subplot(6,6,i+1)
    plt.imshow(full_false[false_selection[30000+false_positive_idx[i]]].reshape(21,21))
    plt.title('{:.2f}'.format(test_prob[false_positive_idx[i]][0]))
plt.tight_layout()
plt.suptitle('False Positive')
plt.subplots_adjust(top=0.95)
fig = plt.figure(figsize=(14,14))
for i in range(36):
    fig.add_subplot(6,6,i+1)
    plt.imshow(full_false[false_selection[30000+true_negative_idx[i]]].reshape(21,21))
    plt.title('{:.2f}'.format(test_prob[true_negative_idx[i]][0]))
plt.tight_layout()
plt.suptitle('True Negative')
plt.subplots_adjust(top=0.95)
fig = plt.figure(figsize=(14,14))
for i in range(36):
    fig.add_subplot(6,6,i+1)
    plt.imshow(full_true[true_selection[30000+false_negative_idx[i]-10000]].reshape(21,21))
    plt.title('{:.2f}'.format(test_prob[false_negative_idx[i]][0]))
plt.tight_layout()
plt.suptitle('False Negative')
plt.subplots_adjust(top=0.95)
accuracy = (true_positive + true_negative) / (true_positive + true_negative + false_positive + false_negative)
print("Accuracy = {:.2f}%".format(100*accuracy))
Explanation: What do the false positives look like?
End of explanation
real_fake = np.load('/epyc/users/smotherh/kbmod_epyc/notebooks/normed_individual_real_false_positives.npy')
real_true = np.load('/epyc/users/smotherh/kbmod_epyc/notebooks/normed_individual_real.npy')
real_true_values = []
for t_im in real_true:
    true_stamp = t_im.reshape(21,21)
    true_stamp -= np.min(true_stamp)
    true_stamp += 1e-12
    #true_kl = special.kl_div(gaussian_test, (true_stamp/np.sum(true_stamp)))
    true_kl = special.rel_entr(gaussian_test, (true_stamp/np.sum(true_stamp)))
    real_true_values.append(np.nansum(true_kl))
real_true_values = np.array(real_true_values)
%%time
real_fake_values = []
for t_im in real_fake:
    fake_stamp = t_im.reshape(21,21)
    fake_stamp -= np.min(fake_stamp)
    fake_stamp += 1e-12
    #true_kl = special.kl_div(gaussian_test, (true_stamp/np.sum(true_stamp)))
    fake_kl = special.rel_entr(gaussian_test, (fake_stamp/np.sum(fake_stamp)))
    real_fake_values.append(np.nansum(fake_kl))
real_fake_values = np.array(real_fake_values)
real_true_results = classifier.predict(real_true_values.reshape(-1,1))
real_true_prob = classifier.predict_proba(real_true_values.reshape(-1,1))
%%time
real_fake_results = classifier.predict(real_fake_values.reshape(-1,1))
real_fake_prob = classifier.predict_proba(real_fake_values.reshape(-1,1))
fig = figure(figsize=(10, 6))
hist(real_fake_prob[:,1])
xlabel('True Object Probability', size=16)
ylabel('Count', size=16)
# Fraction of false positives correctly classified at 0.5 level
(len(real_fake_results) - np.sum(real_fake_results))/len(real_fake_results)
# Fraction of True Objects correctly classified at 0.5 level
np.sum(real_true_results)/len(real_true_results)
real_true_results[:20]
fig = figure(figsize(12,12))
for i in range(len(real_true)):
    fig.add_subplot(4,4,i+1)
    imshow(real_true[i].reshape(21,21))
    title('Prob. True: {:.2f}'.format(real_true_prob[i][1]))
plt.tight_layout()
plt.suptitle('True Sources')
plt.subplots_adjust(top=0.94)
fig = figure(figsize(16,16))
for i in range(49):
    fig.add_subplot(7,7,i+1)
    imshow(real_fake[i].reshape(21,21))
    title('Prob. True: {:.2f}'.format(real_fake_prob[i][1]))
plt.tight_layout()
plt.suptitle('False positives')
plt.subplots_adjust(top=0.94)
misclass = np.where(real_fake_prob[:,1] > 0.95)
len(misclass[0])
fig = figure(figsize(16,16))
for i, idx in list(enumerate(misclass[0]))[:49]:
    fig.add_subplot(7,7,i+1)
    imshow(real_fake[idx].reshape(21,21))
    title('Prob. True: {:.2f}'.format(real_fake_prob[idx][1]))
plt.tight_layout()
plt.suptitle('False positives with highest probability of true source')
plt.subplots_adjust(top=0.94)
Explanation: Using "real" false positives
This section uses a new set of data where the true sources are verified asteroids found in the Lori Allen data. The false sources are actual false positives from KBMOD searches.
We use the training from the synthetic training set in the first part of the notebook and apply it here to the "real" data.
End of explanation
<END_TASK>