Created
January 31, 2010 01:51
-
-
Save Houdini/290833 to your computer and use it in GitHub Desktop.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
import pylab | |
from math import exp | |
def scatter_plot(x, y, r, color='b'): | |
for i in range(len(x)): | |
pylab.scatter(x[i], y[i], s=r, c=color) | |
#x1 = [0, 1., 0.1, 2, 3, 2.3, 3.1, 1.7, 1.9] | |
#y1 = [0.1, 0.2, 0.8, 1, 1.5, 1.9, 2.1, 1.7, 1.9] | |
#x2 = [0.1, 0.5, 2.1, 2.3, 0.8, 2 , 1.2, 1 ,-0.1, 1.5, 0.1] | |
#y2 = [1.5, 3 , 2.9, 2.9, 1.5, 2.5, 3.2, 2.9, 1.8, 1.54, 2.1] | |
x1 = [-3, -4, -5] | |
y1 = [-0.1, 2, 0] | |
x2 = [1, 2, 3] | |
y2 = [0.5, 0, -0.9] | |
eps = 0.01 | |
sigma = 1.3 | |
def get_function_value(x, y, mx, my, sigma_this): | |
value = 1./(2*3.14159265*sigma_this**2)*exp(-0.5/sigma_this**2*((x-mx)**2 + (y-my)**2)) | |
return value | |
x = [i*0.001-3 for i in range(6000)] | |
y = [get_function_value(xi,0,0,0, sigma) for xi in x] | |
def do_one_plot(): | |
scatter_plot(x1, y1, 200) | |
scatter_plot(x2, y2, 200, 'g') | |
pylab.plot(x, y) | |
def get_coef(h, sigma): | |
def find_border(min_or_max, x_or_y): | |
return eval(min_or_max + '([' + min_or_max + '(' + x_or_y + '1),' + min_or_max + '(' + x_or_y + '2)])') | |
left_border = find_border('min', 'x') | |
right_border = find_border('max', 'x') | |
top_border = find_border('max', 'y') | |
bottom_border = find_border('min', 'y') | |
coef = 0 | |
for i in range(int(round((right_border - left_border)/h))): | |
for j in range(int(round((top_border - bottom_border)/h))): | |
center_i = ((i*h + left_border) + h)/2 | |
center_j = ((j*h + bottom_border) + h)/2 | |
sum1 = 0 | |
sum2 = 0 | |
for i_point in range(len(x1)): | |
sum1 += get_function_value(center_i, center_j, x1[i_point], y1[i_point], sigma) | |
for i_point in range(len(x1)): | |
sum2 += get_function_value(center_i, center_j, x2[i_point], y2[i_point], sigma) | |
sum = min([sum1, sum2]) | |
coef += sum*h**2 | |
print "coef = ", coef | |
print "point = ", | |
return coef | |
def plot_full(step_length, all_steps): | |
steps = [i*step_length for i in range(1,all_steps)] | |
coefs = [get_coef(0.01, step) for step in steps] | |
max_coef_i = 0 | |
max_coef = coefs[0] | |
for i in range(len(steps)): | |
if coefs[i] > max_coef: | |
max_coef = coefs[i] | |
max_coef_i = steps[i] | |
print "max_coef_i = ", max_coef_i | |
print "max_coef = ", max_coef | |
pylab.plot(steps, coefs) | |
do_one_plot() | |
#plot_full(0.05, 80) | |
#get_coef(0.01, 0.2) | |
pylab.show() |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment