diff --git a/week7/img/animation.mp4 b/week7/img/animation.mp4 new file mode 100644 index 0000000000000000000000000000000000000000..ec6bca7e480c6a61797cffc1cc089e9214caadbf Binary files /dev/null and b/week7/img/animation.mp4 differ diff --git a/week7/week7_9_2.py b/week7/week7_9_2.py index 8ba9556723814beb9df91ee57baf7f0763553426..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391 100644 --- a/week7/week7_9_2.py +++ b/week7/week7_9_2.py @@ -1,163 +0,0 @@ -import scipy.spatial as spatial -from scipy import integrate as integrate -import numpy as np -import matplotlib.pyplot as plt -import matplotlib.animation as animation -import time - -# import jit from numba -# from numba import jit - -# @jit(nopython=True) -def calc_force(P1, P2): - # calculate the force between points P1 and P2 (both 2D numpy arrays) - # the neutral distance is the variable neutral_distance - dist = np.linalg.norm(P1 - P2) # we should have this from scipy already but recalculating for now... - force = spring_constant * (dist - neutral_distance) * (P1 - P2) / dist - # repulsive forces only - # force = np.where(dist < neutral_distance, force, force*0.1) - return force - -# @jit(nopython=True) -def solve_forces(particles, velocities): - particles_tree = spatial.cKDTree(particles) - # For each point, calculate the force due to the nearest neighbors - F = np.zeros((x_num*y_num, 2)) - for particle in range(x_num*y_num): - # return the index of all particles within radius of radius_check - neighbors = particles_tree.query_ball_point(particles[particle, :], radius_check) - # Moving local damping here to use average velocity to avoid overapplying it - avg_vel = np.mean(velocities[neighbors, :], axis=0) # axis = 0 should average the columns - F_damping = damping_constant * (velocities[particle, :] - avg_vel) - # calculate the force on the particle due its neighbors - for neighbor in neighbors: - if neighbor != particle: - F[particle, :] -= calc_force(particles[particle, :], particles[neighbor, :]) - F[particle, :] -= F_damping # leaving this out here for now - # Add some global damping - F -= damping_constant_global * velocities - F[:, 1] -= gravity # gravity - np.clip(F, -force_threshold, force_threshold, out=F) # clip the force to be within limits of the threshold - return F - -# @jit(nopython=True) -def solve_velocities(F, particles, dt): - # update the velocity of each particle from V_prev according to the force acting on it - V = dt * F / particle_mass * 0.95 # some damping here as well - # if we contact the ground, reverse the y velocity and dampen it - V[:, 1] = np.where((particles[:, 1] < 0), -V[:, 1]*coef_of_restitution, V[:, 1]) # if we cross into the floor, reverse the y velocity and apply a damping factor - # apply friction to the ground - V[:, 0] = np.where((particles[:, 1] < 0), V[:, 0]*0.9, V[:, 0]) - return V - -# @jit(nopython=True) -def solve_positions(particles, V): - # update the position of each particle from particles according to the velocity acting on it - particles_new = particles + dt * V - particles_new[:, 1] = np.where((particles_new[:, 1] < 0), 0, particles_new[:, 1]) # if we cross into the floor, reverse the y velocity and apply a damping factor - # CANTILEVERED BEAM - set all particles in the first column to have a position of zero - ''' - for particle_num in range(len(particles)): - # if the particle is in the first row: - if particle_num <= y_num: - # keep that particle in place - particles_new[particle_num, 0] = particles[particle_num, 0] - particles_new[particle_num, 1] = particles[particle_num, 1] - ''' - return particles_new - -# Constants: -x_num = 40 # 10 or 40 -y_num = 40 # 10 or 40 -t_final = 1.0 -dt = 0.005 -sub_increments = 20 # how many timesteps to calculate between each frame -mesh_spacing = 1 -neutral_distance = mesh_spacing -height_offset = 8 # 3 for 10x10, 8 for 40x40 - -# Physical Constants (these alter the particle physics) -radius_check = neutral_distance*1.35 -spring_constant = 5000 # 1000 -damping_constant = 25 # 10 -damping_constant_global = 0.1 # air resistance -coef_of_restitution = .75 # 0.5 -gravity = 1 -particle_mass = 0.5 -force_threshold_distance = 0.75 # springs shouldn't get any stronger past this distance -force_threshold = force_threshold_distance*neutral_distance * spring_constant - - -# Create our triangular mesh of points -particles = np.zeros((x_num*y_num, 2)) -velocities = np.zeros((x_num*y_num, 2)) -forces = np.zeros((x_num*y_num, 2)) - -# create a mesh of particles, indexes start at point 0,0 and grow up the y-axis, then move to the next column and grow up that one -for i in range(x_num): - for j in range(y_num): - particles[i*y_num+j, 0] = i - particles[i*y_num+j, 1] = j * np.sqrt(3) / 2 # height of equilateral triangle - if j % 2 == 1: # add an offset to every other point - particles[i*y_num+j, 0] += 0.5 - -# particles[-1, 0] -= mesh_spacing*.5 # offset a particle for chaos -particles *= mesh_spacing # scale the mesh to the desired spacing -particles[:,1] += mesh_spacing*height_offset -centroid = np.mean(particles, axis=0) # find the centroid of the mesh - -# rotate the mesh of particles about the origin -theta = np.radians(30) -rot_matrix = np.array([[np.cos(theta), -np.sin(theta)], [np.sin(theta), np.cos(theta)]]) -particles = np.dot(particles-centroid, rot_matrix) + centroid - -# Create a scatterplot of X -my_dpi = 96 -fig = plt.figure(figsize=(800/my_dpi, 400/my_dpi), dpi=my_dpi) -ax = fig.add_axes([0, 0, 1, 1], frameon=False) -# numpy array of colors for each particle in shades from red to blue -colors = np.zeros((x_num*y_num, 4)) -colors[:, 0] = np.linspace(1, 0, x_num*y_num) -colors[:, 2] = np.linspace(0, 1, x_num*y_num) -colors[:, 3] = 1 -plt.axis('equal') # ensure axes are equally scaled -plt.plot([-mesh_spacing*10, x_num*mesh_spacing+mesh_spacing*10], [0, 0], color = 'k') # Draw a line representing the ground at y = 0; -plot = ax.scatter(particles[:, 0], particles[:, 1], - marker = 'o', - s = 20, - color = colors, -) - -time1 = time.time() - - -# @jit(nopython=True) -def animate(w): - global particles, velocities, forces - for i in range(sub_increments): - # Semi-implicit Euler integration -# velocities = velocities + solve_velocities(forces, particles) -# particles = solve_positions(particles, velocities) -# forces = solve_forces(particles, velocities) - # Velocity Verlet integration - velocities_half = velocities + solve_velocities(forces, particles, dt/2) - particles = solve_positions(particles, velocities_half) - forces = solve_forces(particles, velocities) - velocities = velocities_half + solve_velocities(forces, particles, dt/2) - plot.set_offsets(np.c_[particles[:, 0],particles[:, 1]]) - -anim = animation.FuncAnimation(fig, animate, frames=int(2*t_final/dt/2)) -# plt.show() -solve_time = time.time() - time1 - - -time2 = time.time() -import matplotlib as mpl -mpl.rcParams['animation.ffmpeg_path'] = r'C:\ffmpeg_2023_03_30\bin\ffmpeg.exe' -path = r"c://Users/david/Desktop/NMM/week7/img/animation.mp4" -writervideo = animation.FFMpegWriter(fps=30, extra_args=['-vcodec', 'libx264']) -anim.save(path, writer=writervideo, dpi=my_dpi) -mp4_time = time.time() - time2 - -print("Solve time: ", solve_time) -print("MP4 time: ", mp4_time) \ No newline at end of file diff --git a/week9/README.MD b/week9/README.MD new file mode 100644 index 0000000000000000000000000000000000000000..c0966a213ac5262e495b5178882f1c2c07d2734d --- /dev/null +++ b/week9/README.MD @@ -0,0 +1,23 @@ +# Week 9 - Function Fitting + +## 11.1 + +<img src="img/11_1.png" height="600"> + +Fit line: y = 1.648 + 3.772x +Known error: slope error: 0.3401 +Known error: intercept error: 0.1034 +Bootstrap: slope error: 0.0278 +Bootstrap: intercept error: 0.0090 +Ensemble: slope error: 0.0543 +Ensemble: intercept error: 0.0172 + +## 11.2 + +The top plot is our noisy data + fit and the bottom are the a and b coefficients converging over 100 steps of Levenberg-Marquardt iteration. + +<img src="img/11_2.png" height="600"> + +Here's an animation of the coefficients actually updating and converging to fit. + +<img src="img/animation.gif" height="600"> \ No newline at end of file diff --git a/week9/img/11_1.png b/week9/img/11_1.png new file mode 100644 index 0000000000000000000000000000000000000000..33382d65a65ed66ea062ddb36a94da8439eacbbc Binary files /dev/null and b/week9/img/11_1.png differ diff --git a/week9/img/11_2.png b/week9/img/11_2.png new file mode 100644 index 0000000000000000000000000000000000000000..0f21d2ee6a6a1bb8128c7e12918ecc76fdb2f731 Binary files /dev/null and b/week9/img/11_2.png differ diff --git a/week9/img/animation.gif b/week9/img/animation.gif new file mode 100644 index 0000000000000000000000000000000000000000..1bb36167dcbb22ba2abbd3e995018064008d9866 Binary files /dev/null and b/week9/img/animation.gif differ diff --git a/week9/week9_11_1.py b/week9/week9_11_1.py new file mode 100644 index 0000000000000000000000000000000000000000..c8428173ad12794edbca2817c278e9d6281b4b82 --- /dev/null +++ b/week9/week9_11_1.py @@ -0,0 +1,80 @@ +import numpy as np +import matplotlib.pyplot as plt + +def generate_data(size, sigma): + # Generate 100 points uniformly distributed between 0 and 1 + x = np.random.rand(100) + # Ley y - 2+3x + C, where C is a gaussian random variable with standard deviation of 0.5 + C = np.random.normal(0, sigma, size) + y = 2 + 3*x + C + return x, y + +# Use an SVD to fit y = a + bx to this data set, finding a and b. +# https://ltetrel.github.io/data-science/2018/06/08/line_svd.html + + +def fit_line(x, y): + # Make a matrix A with the x values in the first column and 1's in the second column + A = np.vstack([x, y]).T + avg = np.mean(A, axis = 0) + u, s, V = np.linalg.svd(A - avg) # V is a 2x2 matrix with the eigenvectors as columns, with leftmost columns providing the best fit + # get a rise over run slope value from the direction vector V + slope = V[0,1]/V[0,0] + # find the y intercept by plugging in the x value of the average + intercept = avg[1] - slope*avg[0] + return slope, intercept, V, s + +# Generate the data +size = 100 +sigma = 0.5 +x, y = generate_data(size, sigma) +slope, intercept, V, s = fit_line(x, y) + +slope = round(slope, 3) +intercept = round(intercept, 3) +print("Fit line: y = {0} + {1}x".format(intercept, slope)) + +# Plot the data and the fit +plt.title("Fit line: y = {0} + {1}x".format(intercept, slope)) +plt.plot(x, y, 'o', label='Original data', markersize=3) +plt.plot(x, slope*x + intercept, 'r', label='Fitted line') +# set the limits of the plot +plt.show() + +''' +Evaluate the errors in a and b +''' + +# Part A) According to Equation 12.34 +errors = [] +for i in range(2): + error = 0 + for j in range(2): + error += (V[j, i] / s[j]) ** 2 + errors.append(error) +print("Known error: slope error: %.4f" % np.sqrt(errors[0] * (sigma ** 2))) +print("Known error: intercept error: %.4f" % np.sqrt(errors[1] * (sigma ** 2))) + +# Part B) Bootstrap Sampling +space = np.linspace(0, 99, 100) # set up a linespace for each x value +A, B = [], [] +for i in range(100): + sample = np.random.choice(space, 100) + x_b, y_b = x[sample.astype(int)], y[sample.astype(int)] + a, b = np.polyfit(x_b, y_b, 1) + A.append(a) + B.append(b) +print("Bootstrap: slope error: %.4f" % np.var(A)) +print("Bootstrap: intercept error: %.4f" % np.var(B)) + +# Part C) From fitting an ensemble of 100 independent data sets +# Here we are literaly just going to run the same code as part B 100 times +# and then take the variance of the results +A, B = [], [] +for i in range(100): + x, y = generate_data(size, sigma) + a,b,V,s = fit_line(x, y) + A.append(a) + B.append(b) +print("Ensemble: slope error: %.4f" % np.var(A)) +print("Ensemble: intercept error: %.4f" % np.var(B)) \ No newline at end of file diff --git a/week9/week9_11_2.py b/week9/week9_11_2.py new file mode 100644 index 0000000000000000000000000000000000000000..e4ce973340c129c6a96983f1ca99fa58d3ef25bd --- /dev/null +++ b/week9/week9_11_2.py @@ -0,0 +1,91 @@ +import numpy as np +import matplotlib.pyplot as plt +import matplotlib.animation as animation + +''' +Generate 100 points x uniformly distributed between 0 and 1, and let y = sin(2+ 3x) + ζ, where ζ is a Gaussian random variable +with a standard deviation of 0.1. Write a Levenberg-Marquardt routine to fit y = sin(a+bx) to this data set starting from a = b = 1 +(remembering that the second-derivative term can be dropped in the Hessian), and investigate the convergence for both fixed and adaptively +adjusted λ values. +''' + +def generate_data(size, sigma): + # Generate 100 points uniformly distributed between 0 and 1 + x = np.random.rand(100) + # Ley y - 2+3x + C, where C is a gaussian random variable with standard deviation of 0.5 + C = np.random.normal(0, sigma, size) + y = np.sin(2 + 3*x) + C + return x, y + +# Levenberg-Marquardt routine to fit y = sin(a+bx) to this data set starting from a = b = 1 +def levenberg_step(x, y, a, b, l): + # Calculate the Jacobian + J = np.array([np.ones(len(x)), x]).T + # Calculate the Hessian + H = np.dot(J.T, J) + # Calculate the gradient + g = np.dot(J.T, y - np.sin(a + b*x)) + # update = (H + l*I)^-1 * g (where I is the identity matrix np.eye(2) + coefficients = np.dot(np.linalg.inv(H + l*np.eye(2)), g) + return coefficients + +a = 1 +b = 1 +l = 0.1 +steps = 100 + +size = 100 +sigma = 0.1 +x, y = generate_data(size, sigma) + +A = [] +B = [] + +for i in range(steps): + coefficients = levenberg_step(x, y, a, b, l) + A.append(a) + B.append(b) + a += coefficients[0] + b += coefficients[1] + # print("a: %.4f, b: %.4f" % (a, b)) + +''' +# Here we plot the data and the fit line +fig, (ax1, ax2) = plt.subplots(2, 1) +# slim plot +fig.subplots_adjust(hspace=0.5) +ax1.plot(x, y, 'o') +x_sorted = np.sort(x) +ax1.plot(x_sorted, np.sin(a + b*x_sorted)) +# print inital parameters +ax1.text(0,-1, 'Fit line: y = sin(a+bx)\ny = sin(%.2f + %.2fx)' % (a, b), fontsize=10) +ax1.set_title('Levenberg-Marquardt Fit\nsteps = %.2f, sigma = %.2f' % (steps, sigma)) +# Plot A and B as an XY scatter plot with color associated with the value of the iteration number +ax2.scatter(A, B, c=range(steps)) +ax2.set_title('Converging a and b values') +plt.show() +''' + +# Here we animate the fit converging over the scatter plot +my_dpi = 96 +fig = plt.figure(figsize=(800/my_dpi, 400/my_dpi), dpi=my_dpi) +ax = fig.add_subplot(111) +# Scatter plot +scatter = ax.plot(x, y, 'o') +# Fit line +x_sorted = np.sort(x) +fit = ax.plot(x_sorted, np.sin(a + b*x_sorted)) + +# Animate the fit converging over the scatter plot +def animate(w): + global A, B + x = A[w] + y = B[w] + fit[0].set_data(x_sorted, np.sin(x + y*x_sorted)) + +anim = animation.FuncAnimation(fig, animate, frames=int(20)) +# plt.show() + +path = r"c://Users/david/Desktop/NMM/week9/img/animation.gif" +# save as a gif +anim.save(path, writer='imagemagick', fps=8, dpi=my_dpi)