From bc515d3acdd94847b6e7aa6135bc234b46161db6 Mon Sep 17 00:00:00 2001 From: sotech117 Date: Tue, 7 May 2024 07:00:43 -0400 Subject: add hw9 and hw8 --- .DS_Store | Bin 14340 -> 12292 bytes hw8/10-14-0.0001.png | Bin 0 -> 17051 bytes hw8/10-14-0.001.png | Bin 0 -> 60938 bytes hw8/10-14-0.01.png | Bin 0 -> 19005 bytes hw8/10-14-0.1.png | Bin 0 -> 17313 bytes hw8/10-14-E.png | Bin 0 -> 14507 bytes hw8/10-14-psi.mp4 | Bin 0 -> 3519376 bytes hw8/10-14-psi.png | Bin 0 -> 61965 bytes hw8/10-14.jl | 68 ++++++-- hw8/10-17-1.png | Bin 0 -> 10080 bytes hw8/10-17-2.05.png | Bin 0 -> 24566 bytes hw8/10-17-2.png | Bin 0 -> 30436 bytes hw8/10-17-3.05.png | Bin 0 -> 28064 bytes hw8/10-17-3.png | Bin 0 -> 29569 bytes hw8/10-17-width-.05.png | Bin 0 -> 46547 bytes hw8/10-17-width-.1.png | Bin 0 -> 28921 bytes hw8/10-17-width-.2.png | Bin 0 -> 30538 bytes hw8/10-17.jl | 64 +++++--- hw8/10-3-long.png | Bin 0 -> 27534 bytes hw8/10-3.jl | 39 ++++- hw8/10-3.png | Bin 0 -> 60650 bytes hw8/TimeDependentSchrodingerEquation.jl | 38 ++--- hw9/12-11-all.png | Bin 0 -> 50155 bytes hw9/12-11-all.svg | 131 +++++++++++++++ hw9/12-11.jl | 272 ++++++++++++++++++++++++++++++++ hw9/12-12-all.png | Bin 0 -> 32142 bytes hw9/12-12-test.png | Bin 0 -> 29601 bytes hw9/12-12.jl | 253 +++++++++++++++++++++++++++++ hw9/YaoQuantumComputing.jl | 28 ++++ hw9/qubit00.svg | 39 +++++ hw9/qubit01.svg | 39 +++++ hw9/qubit10.svg | 39 +++++ hw9/qubit11.svg | 39 +++++ hw9/yao-assigment-2.jl | 75 +++++++++ 34 files changed, 1067 insertions(+), 57 deletions(-) create mode 100644 hw8/10-14-0.0001.png create mode 100644 hw8/10-14-0.001.png create mode 100644 hw8/10-14-0.01.png create mode 100644 hw8/10-14-0.1.png create mode 100644 hw8/10-14-E.png create mode 100644 hw8/10-14-psi.mp4 create mode 100644 hw8/10-14-psi.png create mode 100644 hw8/10-17-1.png create mode 100644 hw8/10-17-2.05.png create mode 100644 hw8/10-17-2.png create mode 100644 hw8/10-17-3.05.png create mode 100644 hw8/10-17-3.png create mode 100644 hw8/10-17-width-.05.png create mode 100644 hw8/10-17-width-.1.png create mode 100644 hw8/10-17-width-.2.png create mode 100644 hw8/10-3-long.png create mode 100644 hw8/10-3.png create mode 100644 hw9/12-11-all.png create mode 100644 hw9/12-11-all.svg create mode 100644 hw9/12-11.jl create mode 100644 hw9/12-12-all.png create mode 100644 hw9/12-12-test.png create mode 100644 hw9/12-12.jl create mode 100644 hw9/YaoQuantumComputing.jl create mode 100644 hw9/qubit00.svg create mode 100644 hw9/qubit01.svg create mode 100644 hw9/qubit10.svg create mode 100644 hw9/qubit11.svg create mode 100644 hw9/yao-assigment-2.jl diff --git a/.DS_Store b/.DS_Store index cd84abd..2b5a296 100644 Binary files a/.DS_Store and b/.DS_Store differ diff --git a/hw8/10-14-0.0001.png b/hw8/10-14-0.0001.png new file mode 100644 index 0000000..0dc656a Binary files /dev/null and b/hw8/10-14-0.0001.png differ diff --git a/hw8/10-14-0.001.png b/hw8/10-14-0.001.png new file mode 100644 index 0000000..fcf13d4 Binary files /dev/null and b/hw8/10-14-0.001.png differ diff --git a/hw8/10-14-0.01.png b/hw8/10-14-0.01.png new file mode 100644 index 0000000..8c6d987 Binary files /dev/null and b/hw8/10-14-0.01.png differ diff --git a/hw8/10-14-0.1.png b/hw8/10-14-0.1.png new file mode 100644 index 0000000..74f5045 Binary files /dev/null and b/hw8/10-14-0.1.png differ diff --git a/hw8/10-14-E.png b/hw8/10-14-E.png new file mode 100644 index 0000000..d13b721 Binary files /dev/null and b/hw8/10-14-E.png differ diff --git a/hw8/10-14-psi.mp4 b/hw8/10-14-psi.mp4 new file mode 100644 index 0000000..37c0d10 Binary files /dev/null and b/hw8/10-14-psi.mp4 differ diff --git a/hw8/10-14-psi.png b/hw8/10-14-psi.png new file mode 100644 index 0000000..7911c89 Binary files /dev/null and b/hw8/10-14-psi.png differ diff --git a/hw8/10-14.jl b/hw8/10-14.jl index 27d1b5c..8a6307a 100644 --- a/hw8/10-14.jl +++ b/hw8/10-14.jl @@ -52,9 +52,12 @@ plot(prob(psi0)) # integrate forward in time tf = 1.0 -dt = 0.1 +dt = 0.0007 tspan = (0.0, tf) +println("dx = ", dx) +println("dt = ", dt) + function timeEvolve(psi0, tf, dt) # second order Runge-Kutta algorithm psi = psi0 for t in range(0, stop = tf, step = dt) @@ -64,19 +67,56 @@ function timeEvolve(psi0, tf, dt) # second order Runge-Kutta algorithm return psi end -# tendency(psi, dx, t) = derivative(psi, dx) # use ODE solver -# problem = ODEProblem(tendency, psi0, tspan, dx) # specify ODE -# sol = solve(problem, Tsit5(), reltol = 1e-12, abstol = 1e-12) # solve using Tsit5 algorithm to specified accuracy +# get energy of psi +function energy(psi, dx) + e_complex = dot(psi, H(psi, dx)) + return real(e_complex) +end # compare initial and final wavefunction probabilities -psi = timeEvolve(psi0, tf, dt) -# psi = sol[:, end] -times = sol.t - +psi_1 = timeEvolve(psi0, tf, dt) +psi_2 = timeEvolve(psi0, tf, 0.00088) +psi_3 = timeEvolve(psi0, tf, 0.00089) +psi_4 = timeEvolve(psi0, tf, 0.0009) # check that normalization hasn't deviated too far from 1.0 -println("norm = ", normalization(psi, dx)) - -plot([prob(psi0), prob(psi)]) -savefig("10-14.png") - - +println("norm euler dt=.0007 ->", normalization(psi_1, dx)) +println("norm euler dt=.00088 ->", normalization(psi_2, dx)) +println("norm euler dt=.00089 -> ", normalization(psi_3, dx)) +println("norm euler dt=.0009 -> ", normalization(psi_4, dx)) + + +tendency(psi, dx, t) = derivative(psi, dx) # use ODE solver +problem = ODEProblem(tendency, psi0, tspan, dx) # specify ODE +sol = solve(problem, Tsit5(), reltol = 1e-13, abstol = 1e-13) # solve using Tsit5 algorithm to specified accuracy +psi_5 = sol[:, end] +times_2 = sol.t +println("norm ode :", normalization(psi_5, dx)) + +# cluclate the energy at each timestep +# energies_over_time = [energy(sol[:, i], dx) for i in 1:length(times_2)] +# plot(times_2, energies_over_time, title = "Energy of Wave Packet over Time", xlabel = "t", ylabel = "E", lw = 1.5, ylim = (85, 90), legend = false) +# savefig("hw8/10-14-E.png") + +# make an animation of the wave packet +# anim = @animate for i in 1:length(times_2) +# plot(x, prob(sol[:, i]), title = "Propagation of Wave Packet", xlabel = "x", ylabel = "|ψ(x)|²", lw = 1.5, ylim = (0, 1), label = "t = $(round(times_2[i], digits = 3))", legend = :topright) +# end +# mp4(anim, "hw8/10-14-psi.mp4", fps = 75) + +# plot([prob(psi0), prob(psi_1), prob(psi_2)], label = ["Initial" "Final (Euler dt = $dt) @ t=$tf" "Final (ODE abserr = 10^(-12)) @ t=$tf"], lw = 1.5, title = "Propagation of Wave Packet", xlabel = "x", ylabel = "|ψ(x)|²") +# plot( +# [prob(psi0), prob(psi_1), prob(psi_2), prob(psi_3), prob(psi_4)], +# label = ["Initial" "Final (Euler dt = $dt) @ t=$tf" "Final (Euler dt = 0.00088) @ t=$tf" "Final (Euler dt = 0.00089) @ t=$tf" "Final (ODE abserr = 10^(-12)) @ t=$tf"], +# lw = 1.5, +# title = "Errors on Propagation of Wave Packet k_0 = 1000.0", +# xlabel = "x", +# ) +plot( + [prob(psi0), prob(psi_1), prob(psi_2), prob(psi_3), prob(psi_4), prob(psi_5)], + label = ["Initial" "Final (Euler dt = $dt) @ t=$tf" "Final (Euler dt = 0.00088) @ t=$tf" "Final (Euler dt = 0.00089) @ t=$tf" "Final (Euler dt = 0.0009) @ t=$tf" "Final (ODE abserr = 10^(-12)) @ t=$tf"], + lw = 1.5, + title = "Propagation of Wave Packet (k_0 = 1000.0)", + xlabel = "x", + ylabel = "|ψ(x)|²", +) +savefig("hw8/10-14-psi.png") diff --git a/hw8/10-17-1.png b/hw8/10-17-1.png new file mode 100644 index 0000000..5646660 Binary files /dev/null and b/hw8/10-17-1.png differ diff --git a/hw8/10-17-2.05.png b/hw8/10-17-2.05.png new file mode 100644 index 0000000..885a984 Binary files /dev/null and b/hw8/10-17-2.05.png differ diff --git a/hw8/10-17-2.png b/hw8/10-17-2.png new file mode 100644 index 0000000..7453918 Binary files /dev/null and b/hw8/10-17-2.png differ diff --git a/hw8/10-17-3.05.png b/hw8/10-17-3.05.png new file mode 100644 index 0000000..d18dc6d Binary files /dev/null and b/hw8/10-17-3.05.png differ diff --git a/hw8/10-17-3.png b/hw8/10-17-3.png new file mode 100644 index 0000000..e140e1c Binary files /dev/null and b/hw8/10-17-3.png differ diff --git a/hw8/10-17-width-.05.png b/hw8/10-17-width-.05.png new file mode 100644 index 0000000..61829cc Binary files /dev/null and b/hw8/10-17-width-.05.png differ diff --git a/hw8/10-17-width-.1.png b/hw8/10-17-width-.1.png new file mode 100644 index 0000000..b2e77cb Binary files /dev/null and b/hw8/10-17-width-.1.png differ diff --git a/hw8/10-17-width-.2.png b/hw8/10-17-width-.2.png new file mode 100644 index 0000000..ec49ca7 Binary files /dev/null and b/hw8/10-17-width-.2.png differ diff --git a/hw8/10-17.jl b/hw8/10-17.jl index 6dfff2b..d1d18b7 100644 --- a/hw8/10-17.jl +++ b/hw8/10-17.jl @@ -4,20 +4,29 @@ using Plots using LinearAlgebra using DifferentialEquations +k_0 = 700.0 # useful functions: -function H(psi, dx) # action of Hamiltonian on wavefunction +function H(psi, dx, k = 0, width = 0.1) # action of Hamiltonian on wavefunction Hpsi = zeros(ComplexF64, size(psi)) # -(1/2) * laplacian(psi) (m = hbar = 1) Hpsi[2:end-1] = 0.5 * (2.0 * psi[2:end-1] - psi[3:end] - psi[1:end-2]) / (dx * dx) + # add the potential function V_0 = 2 * k_0^2 for specific region + for i in 1:length(psi) + x = i * dx + if x >= 30.0 && x <= 30.0 + width + Hpsi[i] = 2 * k^2 * psi[i] + end + end + # periodic boundary conditions #Hpsi[1] = 0.5 * (2.0*psi[1] - psi[2] - psi[end])/(dx*dx) #Hpsi[end] = 0.5 * (2.0*psi[end] - psi[end-1] - psi[1])/(dx*dx) return Hpsi end -derivative(psi, dx) = -1.0im * H(psi, dx) +derivative(psi, dx) = -1.0im * H(psi, dx, k_0) function initialWavefunction(x::Vector{Float64}, x0 = 10.0, Delta = 1.0, k = 4.0) Delta2 = Delta^2 @@ -32,15 +41,16 @@ end prob(psi) = real(psi .* conj(psi)) # The actual simulation -N = 400 # number of lattice points -L = 20.0 # x runs from 0 to L +N = 5000 # number of lattice points +L = 40.0 # x runs from 0 to L dx = L / N x = range(0.0, L, N) |> collect # lattice along x-axis -#println(x) +println("dx = ", dx) + # initial wavefunction has position (x0), width (Delta), and momentum (k) -psi0 = initialWavefunction(x, 10.0, 0.05, 700.0) +psi0 = initialWavefunction(x, 10.0, 0.05, k_0) # normalize wavefunction n = normalization(psi0, dx) @@ -51,31 +61,49 @@ println("norm = ", normalization(psi0, dx)) plot(prob(psi0)) # integrate forward in time -tf = 1.0 +tf = 3.0 dt = 5e-7 tspan = (0.0, tf) -function timeEvolve(psi0, tf, dt) # second order Runge-Kutta algorithm - psi = psi0 - for t in range(0, stop = tf, step = dt) - psiMid = psi + 0.5 * dt * derivative(psi, dx) - psi = psi + dt * derivative(psiMid, dx) - end - return psi -end +# function timeEvolve(psi0, tf, dt) # second order Runge-Kutta algorithm +# psi = psi0 +# for t in range(0, stop = tf, step = dt) +# psiMid = psi + 0.5 * dt * derivative(psi, dx) +# psi = psi + dt * derivative(psiMid, dx) + +# # println("t = ", t, " norm = ", real(normalization(psi, dx))) + +# if real(normalization(psi, dx)) < 0.999 || real(normalization(psi, dx)) > 1.010 +# println("dt = ", dt, " t = ", t, " norm = ", normalization(psi, dx)) +# println("Normalization deviated too far from 1.0") +# break +# end +# end +# return psi +# end tendency(psi, dx, t) = derivative(psi, dx) # use ODE solver problem = ODEProblem(tendency, psi0, tspan, dx) # specify ODE -sol = solve(problem, Tsit5(), reltol = 1e-12, abstol = 1e-12) # solve using Tsit5 algorithm to specified accuracy +sol = solve(problem, Tsit5(), reltol = 1e-14, abstol = 1e-14) # solve using Tsit5 algorithm to specified accuracy # compare initial and final wavefunction probabilities -#psi = timeEvolve(psi0, tf, dt) +# psi = timeEvolve(psi0, tf, dt) psi = sol[:, end] times = sol.t # check that normalization hasn't deviated too far from 1.0 println("norm = ", normalization(psi, dx)) -plot([prob(psi0), prob(psi)]) +plot(prob(psi0)) +savefig("10-17-1.png") +# plot a vertical line where the barrier is +# barrier_on_lattice = 6.5 / dx +barrier_on_lattice = 30.0 / dx +plot(prob(psi), label = "final (t=$tf)", xlabel = "x", ylabel = "|ψ|^2", title = "Initial and final probability densities k_0=$k_0", lw = 1.5) +plot!([barrier_on_lattice, barrier_on_lattice], [0.0, 0.3], color = :red, label = "barrier (width .2)", lw = 1.5, linestyle = :dash) +savefig("10-17-width-.1.png") +# plot([prob(psi0), prob(psi)], label = ["initial" "final (t=$tf)"], xlabel = "x", ylabel = "probability density", title = "Initial and final probability densities k_0=$k_0", lw = 1.5) +# plot!([barrier_on_lattice, barrier_on_lattice], [0.0, 2.0], color = :red, label = "barrier (width .2)", lw = 1.5, linestyle = :dash) +# savefig("10-17-width-.2.png") diff --git a/hw8/10-3-long.png b/hw8/10-3-long.png new file mode 100644 index 0000000..8d782a4 Binary files /dev/null and b/hw8/10-3-long.png differ diff --git a/hw8/10-3.jl b/hw8/10-3.jl index 8616a69..2f4a690 100644 --- a/hw8/10-3.jl +++ b/hw8/10-3.jl @@ -60,16 +60,43 @@ function map_n_to_energies(n) return e end -n_max = 18 -n_to_e = [map_n_to_energies(n) for n in 1:n_max] +n_s = collect(0:1:18) +n_to_e = [map_n_to_energies(n) for n in n_s] # plot e[0] for all N -eList = zeros(0) -for i in 1:n_max - push!(eList, n_to_e[i][1]) +ground_state = [] +excited_1 = [] +excited_2 = [] +excited_3 = [] +for i in 1:length(n_to_e) + push!(ground_state, n_to_e[i][1]) + push!(excited_1, n_to_e[i][2]) + push!(excited_2, n_to_e[i][3]) + push!(excited_3, n_to_e[i][4]) end -plot(eList) +plot(ground_state, label = "groud state energy for n", xlabel = "n (level)", ylabel = "energy", title = "excited energy levels for V(n) = abs(x)^n", marker = :circle) +plot!(excited_1, label = "1st excited state", marker = :circle) +plot!(excited_2, label = "2nd excited state", marker = :circle) +plot!(excited_3, label = "3rd excited state", marker = :circle) + +# plot the energies for an inifinite square well as a horizontial line +# function excited_state_to_energy_inf_square_well(n) +# return n^2 * pi^2 / 2 +# end + +# ground_state_inf_square_well = [excited_state_to_energy_inf_square_well(1) for i in 1:length(n_s)] +# excited_1_inf_square_well = [excited_state_to_energy_inf_square_well(2) for i in 1:length(n_s)] +# excited_2_inf_square_well = [excited_state_to_energy_inf_square_well(3) for i in 1:length(n_s)] +# excited_3_inf_square_well = [excited_state_to_energy_inf_square_well(4) for i in 1:length(n_s)] + +# plot!(ground_state_inf_square_well, label = "ground state energy for infinite square well") +# plot!(excited_1_inf_square_well, label = "1st excited state for infinite square well") +# plot!(excited_2_inf_square_well, label = "2nd excited state for infinite square well") +# plot!(excited_3_inf_square_well, label = "3rd excited state for infinite square well") + + +savefig("hw8/10-3.png") # gs(x) = exp(-0.5 * x^2) # Gaussian that is exact ground state of SHO diff --git a/hw8/10-3.png b/hw8/10-3.png new file mode 100644 index 0000000..23c15a9 Binary files /dev/null and b/hw8/10-3.png differ diff --git a/hw8/TimeDependentSchrodingerEquation.jl b/hw8/TimeDependentSchrodingerEquation.jl index dacbced..1d4f315 100644 --- a/hw8/TimeDependentSchrodingerEquation.jl +++ b/hw8/TimeDependentSchrodingerEquation.jl @@ -7,26 +7,26 @@ using DifferentialEquations # useful functions: function H(psi, dx) # action of Hamiltonian on wavefunction - Hpsi = zeros(ComplexF64, size(psi)) - # -(1/2) * laplacian(psi) (m = hbar = 1) - Hpsi[2:end-1] = 0.5 * (2.0*psi[2:end-1] - psi[3:end] - psi[1:end-2])/(dx*dx) - - # periodic boundary conditions - #Hpsi[1] = 0.5 * (2.0*psi[1] - psi[2] - psi[end])/(dx*dx) - #Hpsi[end] = 0.5 * (2.0*psi[end] - psi[end-1] - psi[1])/(dx*dx) - return Hpsi + Hpsi = zeros(ComplexF64, size(psi)) + # -(1/2) * laplacian(psi) (m = hbar = 1) + Hpsi[2:end-1] = 0.5 * (2.0 * psi[2:end-1] - psi[3:end] - psi[1:end-2]) / (dx * dx) + + # periodic boundary conditions + #Hpsi[1] = 0.5 * (2.0*psi[1] - psi[2] - psi[end])/(dx*dx) + #Hpsi[end] = 0.5 * (2.0*psi[end] - psi[end-1] - psi[1])/(dx*dx) + return Hpsi end derivative(psi, dx) = -1.0im * H(psi, dx) function initialWavefunction(x::Vector{Float64}, x0 = 10.0, Delta = 1.0, k = 4.0) - Delta2 = Delta^2 - return exp.(- (x .- x0).^2 / Delta2 ) .* exp.(1.0im * k * x) + Delta2 = Delta^2 + return exp.(-(x .- x0) .^ 2 / Delta2) .* exp.(1.0im * k * x) end function normalization(psi, dx) # normalization of wavefunction - n2 = dot(psi, psi) * dx - return sqrt(n2) + n2 = dot(psi, psi) * dx + return sqrt(n2) end prob(psi) = real(psi .* conj(psi)) @@ -56,17 +56,17 @@ dt = 0.0001 tspan = (0.0, tf) function timeEvolve(psi0, tf, dt) # second order Runge-Kutta algorithm - psi = psi0 - for t in range(0, stop=tf, step=dt) - psiMid = psi + 0.5 * dt * derivative(psi, dx) - psi = psi + dt * derivative(psiMid, dx) - end - return psi + psi = psi0 + for t in range(0, stop = tf, step = dt) + psiMid = psi + 0.5 * dt * derivative(psi, dx) + psi = psi + dt * derivative(psiMid, dx) + end + return psi end tendency(psi, dx, t) = derivative(psi, dx) # use ODE solver problem = ODEProblem(tendency, psi0, tspan, dx) # specify ODE -sol = solve(problem, Tsit5(), reltol=1e-12, abstol=1e-12) # solve using Tsit5 algorithm to specified accuracy +sol = solve(problem, Tsit5(), reltol = 1e-12, abstol = 1e-12) # solve using Tsit5 algorithm to specified accuracy # compare initial and final wavefunction probabilities #psi = timeEvolve(psi0, tf, dt) diff --git a/hw9/12-11-all.png b/hw9/12-11-all.png new file mode 100644 index 0000000..ac5589d Binary files /dev/null and b/hw9/12-11-all.png differ diff --git a/hw9/12-11-all.svg b/hw9/12-11-all.svg new file mode 100644 index 0000000..c6bc3c6 --- /dev/null +++ b/hw9/12-11-all.svg @@ -0,0 +1,131 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + \ No newline at end of file diff --git a/hw9/12-11.jl b/hw9/12-11.jl new file mode 100644 index 0000000..ee7e2b1 --- /dev/null +++ b/hw9/12-11.jl @@ -0,0 +1,272 @@ +# my neural net from 12.11 + +# first, store the memory of A as a lattice of 1, -1, where 1 is the memory and -1 is the absence of memory +A = [ + -1 -1 -1 -1 1 1 -1 -1 -1 -1; + -1 -1 -1 1 1 1 1 -1 -1 -1; + -1 -1 1 1 -1 -1 1 1 -1 -1; + -1 1 1 -1 -1 -1 -1 1 1 -1; + 1 1 -1 -1 -1 -1 -1 -1 1 1; + 1 1 -1 -1 -1 -1 -1 -1 1 1; + -1 1 1 -1 -1 -1 -1 1 1 -1; + -1 -1 1 1 -1 -1 1 1 -1 -1; + -1 -1 -1 1 1 1 1 -1 -1 -1; + -1 -1 -1 -1 1 1 -1 -1 -1 -1 +] + +B = [ + 1 1 1 1 1 -1 -1 -1 -1 -1; + 1 1 -1 -1 1 1 -1 -1 -1 -1; + 1 1 -1 -1 1 1 -1 -1 -1 -1; + 1 1 -1 -1 1 1 -1 -1 -1 -1; + 1 1 1 1 1 -1 -1 -1 -1 -1; + 1 1 1 1 1 -1 -1 -1 -1 -1; + 1 1 -1 -1 1 1 -1 -1 -1 -1; + 1 1 -1 -1 1 1 -1 -1 -1 -1; + 1 1 -1 -1 1 1 -1 -1 -1 -1; + 1 1 1 1 1 -1 -1 -1 -1 -1 +] + +C = [ + -1 1 1 1 1 1 1 -1 -1 -1; + 1 1 1 1 1 1 1 -1 -1 -1; + 1 1 1 -1 -1 -1 -1 -1 -1 -1; + 1 1 -1 -1 -1 -1 -1 -1 -1 -1; + 1 1 -1 -1 -1 -1 -1 -1 -1 -1; + 1 1 -1 -1 -1 -1 -1 -1 -1 -1; + 1 1 -1 -1 -1 -1 -1 -1 -1 -1; + 1 1 1 -1 -1 -1 -1 -1 -1 -1; + 1 1 1 1 1 1 1 -1 -1 -1; + -1 1 1 1 1 1 1 -1 -1 -1 +] + +D = [ + 1 1 1 1 1 -1 -1 -1 -1 -1; + 1 1 -1 -1 1 1 -1 -1 -1 -1; + 1 1 -1 -1 1 1 -1 -1 -1 -1; + 1 1 -1 -1 1 1 -1 -1 -1 -1; + 1 1 -1 -1 1 1 -1 -1 -1 -1; + 1 1 -1 -1 1 1 -1 -1 -1 -1; + 1 1 -1 -1 1 1 -1 -1 -1 -1; + 1 1 -1 -1 1 1 -1 -1 -1 -1; + 1 1 -1 -1 1 1 -1 -1 -1 -1; + 1 1 1 1 1 -1 -1 -1 -1 -1 +] + +E_letter = [ + 1 1 1 1 1 1 1 1 1 1; + 1 1 -1 -1 -1 -1 -1 -1 -1 -1; + 1 1 -1 -1 -1 -1 -1 -1 -1 -1; + 1 1 1 1 1 1 1 1 1 1; + 1 1 -1 -1 -1 -1 -1 -1 -1 -1; + 1 1 -1 -1 -1 -1 -1 -1 -1 -1; + 1 1 1 1 1 1 1 1 1 1; + 1 1 -1 -1 -1 -1 -1 -1 -1 -1; + 1 1 -1 -1 -1 -1 -1 -1 -1 -1; + 1 1 1 1 1 1 1 1 1 1 +] + +F = [ + 1 1 1 1 1 1 1 1 1 1; + 1 1 -1 -1 -1 -1 -1 -1 -1 -1; + 1 1 -1 -1 -1 -1 -1 -1 -1 -1; + 1 1 1 1 1 1 1 1 1 1; + 1 1 -1 -1 -1 -1 -1 -1 -1 -1; + 1 1 -1 -1 -1 -1 -1 -1 -1 -1; + 1 1 -1 -1 -1 -1 -1 -1 -1 -1; + 1 1 -1 -1 -1 -1 -1 -1 -1 -1; + 1 1 -1 -1 -1 -1 -1 -1 -1 -1; + 1 1 -1 -1 -1 -1 -1 -1 -1 -1 +] + +G = [ + 1 1 1 1 1 1 1 1 1 1; + 1 1 -1 -1 -1 -1 -1 -1 -1 -1; + 1 1 -1 -1 -1 -1 -1 -1 -1 -1; + 1 1 -1 -1 -1 -1 -1 -1 -1 -1; + 1 1 -1 -1 -1 -1 -1 1 1 1; + 1 1 -1 -1 -1 -1 -1 -1 1 1; + 1 1 -1 -1 -1 -1 -1 -1 1 1; + 1 1 -1 -1 -1 -1 -1 -1 1 1; + -1 1 1 1 1 1 1 1 1 1; + -1 -1 1 1 1 1 1 1 1 1 +] + +encodings = [A, B, F, G] + +spin_numbers(row, col) = (row - 1) * 10 + col + +# create strength of interactions between ith and jth spinds in Ji,j +J = zeros(100, 100) +for m in 1:100 # row + for n in 1:100 # col + i = (m - 1) % 10 + 1 + j = (n - 1) % 10 + 1 + k = (m - 1) ÷ 10 + 1 + l = (n - 1) ÷ 10 + 1 + + J[m, n] = 0 + for encoding in encodings + J[m, n] += encoding[i, k] * encoding[j, l] + end + J[m, n] /= length(encodings) + end +end + +# the energy function +function energy(s, J) + E = 0 + for i in 1:100 + for j in 1:100 + E += J[i, j] * s[i] * s[j] + end + end + return -E +end + +# the update function (uses monte carlo steps) +function monte_carlo(s, J) + for i in 1:100 # systematically go through each point in the lattice + # calculate the energy of the system + E = energy(s, J) + # randomly flip a spin + s[i] = -s[i] + # calculate the new energy of the system + E_new = energy(s, J) + # calculate the change in energy + dE = E_new - E + # if the change in energy is positive, flip the spin back + if dE > 0 + s[i] = -s[i] + end + end + return s +end + +# create the main function +function main(s, J, nsteps) + E = energy(s, J) + for i in 1:nsteps + s = monte_carlo(s, J) + E = energy(s, J) + end + return s, E +end + +# run the main function + +# randomly change some values in A to see if NN works +function produce_test_arr(enc, prob_change = 0.1) + tmp = zeros(10, 10) + for i in 1:10 + for j in 1:10 + if rand() < prob_change + tmp[i, j] = rand([-1, 1]) # set some random values + else + tmp[i, j] = enc[i, j] + end + end + end + return tmp +end + +function check_if_same(s, enc) + for i in 1:10 + for j in 1:10 + if s[i, j] != enc[i, j] + println("The neural net did not work") + return false + end + end + end + # println("The neural net worked") + return true +end + +function apply_damage_to_J(J, prob_damage = 0.8) + println("Applying damage to J with probability ", prob_damage) + ret = copy(J) + for i in 1:100 + for j in 1:100 + if rand() < prob_damage + ret[i, j] = 0 + else + ret[i, j] = J[i, j] + end + end + end + return ret +end + +function run_tests(num_times, Js, MC_steps = 2) + num_correct = 0 + num_total = 0 + + for i in 1:num_times + # randomly select a memory + enc = encodings[rand(1:length(encodings))] + # produce a test array + s = produce_test_arr(enc) + # run theNN + s, E = main(s, Js, MC_steps) + # check if the NN worked + if check_if_same(s, enc) + num_correct += 1 + end + num_total += 1 + end + + println("Number of correct tests: ", num_correct) + + return num_correct / num_total +end + +function run_tests_for_damage_range(damages, num_times, MC_steps = 2, init_J = J) + res = [] + + for damage in damages + damaged_J = apply_damage_to_J(init_J, damage) + p = run_tests(num_times, damaged_J, MC_steps) + + push!(res, p) + end + + return res +end + + +# J = apply_damage_to_J(J, 0.8) + +# s = produce_test_arr(F) +# s, E = main(s, J, 2) +# check_if_same(s, F) + +s = produce_test_arr(A) +s, E = main(s, J, 2) +check_if_same(s, A) + +# s = produce_test_arr(C, 0.0) +# s, E = main(s, J, 2) +# check_if_same(s, C) + +damage_range = collect(0.1:0.05:1.0) +res = run_tests_for_damage_range(damage_range, 100) + +#plot the results +using Plots +plot(damage_range, res, xlabel = "Damage", ylabel = "Success rate", title = "NN success rate vs damage probability", marker = :circle, label = "MC steps = 2") + +# do the same for different number of MC steps +res = run_tests_for_damage_range(damage_range, 100, 5) +plot!(damage_range, res, xlabel = "Damage", ylabel = "Success rate", title = "NN success rate vs damage probability", marker = :circle, label = "MC steps = 5") + +res = run_tests_for_damage_range(damage_range, 100, 10) +plot!(damage_range, res, xlabel = "Damage", ylabel = "Success rate", title = "NN success rate vs damage probability", marker = :circle, label = "MC steps = 10") + +res = run_tests_for_damage_range(damage_range, 100, 20) +plot!(damage_range, res, xlabel = "Damage", ylabel = "Success rate", title = "NN success rate vs damage probability", marker = :circle, label = "MC steps = 20") + +savefig("12-11-all.png") + + + diff --git a/hw9/12-12-all.png b/hw9/12-12-all.png new file mode 100644 index 0000000..edfb678 Binary files /dev/null and b/hw9/12-12-all.png differ diff --git a/hw9/12-12-test.png b/hw9/12-12-test.png new file mode 100644 index 0000000..f86297e Binary files /dev/null and b/hw9/12-12-test.png differ diff --git a/hw9/12-12.jl b/hw9/12-12.jl new file mode 100644 index 0000000..e08eadc --- /dev/null +++ b/hw9/12-12.jl @@ -0,0 +1,253 @@ +# my neural net from 12.11 + +# first, store the memory of A as a lattice of 1, -1, where 1 is the memory and -1 is the absence of memory +A = [ + -1 -1 -1 -1 1 1 -1 -1 -1 -1; + -1 -1 -1 1 1 1 1 -1 -1 -1; + -1 -1 1 1 -1 -1 1 1 -1 -1; + -1 1 1 -1 -1 -1 -1 1 1 -1; + 1 1 -1 -1 -1 -1 -1 -1 1 1; + 1 1 -1 -1 -1 -1 -1 -1 1 1; + -1 1 1 -1 -1 -1 -1 1 1 -1; + -1 -1 1 1 -1 -1 1 1 -1 -1; + -1 -1 -1 1 1 1 1 -1 -1 -1; + -1 -1 -1 -1 1 1 -1 -1 -1 -1 +] + +B = [ + 1 1 1 1 1 -1 -1 -1 -1 -1; + 1 1 -1 -1 1 1 -1 -1 -1 -1; + 1 1 -1 -1 1 1 -1 -1 -1 -1; + 1 1 -1 -1 1 1 -1 -1 -1 -1; + 1 1 1 1 1 -1 -1 -1 -1 -1; + 1 1 1 1 1 -1 -1 -1 -1 -1; + 1 1 -1 -1 1 1 -1 -1 -1 -1; + 1 1 -1 -1 1 1 -1 -1 -1 -1; + 1 1 -1 -1 1 1 -1 -1 -1 -1; + 1 1 1 1 1 -1 -1 -1 -1 -1 +] + +C = [ + -1 1 1 1 1 1 1 -1 -1 -1; + 1 1 1 1 1 1 1 -1 -1 -1; + 1 1 1 -1 -1 -1 -1 -1 -1 -1; + 1 1 -1 -1 -1 -1 -1 -1 -1 -1; + 1 1 -1 -1 -1 -1 -1 -1 -1 -1; + 1 1 -1 -1 -1 -1 -1 -1 -1 -1; + 1 1 -1 -1 -1 -1 -1 -1 -1 -1; + 1 1 1 -1 -1 -1 -1 -1 -1 -1; + 1 1 1 1 1 1 1 -1 -1 -1; + -1 1 1 1 1 1 1 -1 -1 -1 +] + +D = [ + 1 1 1 1 1 -1 -1 -1 -1 -1; + 1 1 -1 -1 1 1 -1 -1 -1 -1; + 1 1 -1 -1 1 1 -1 -1 -1 -1; + 1 1 -1 -1 1 1 -1 -1 -1 -1; + 1 1 -1 -1 1 1 -1 -1 -1 -1; + 1 1 -1 -1 1 1 -1 -1 -1 -1; + 1 1 -1 -1 1 1 -1 -1 -1 -1; + 1 1 -1 -1 1 1 -1 -1 -1 -1; + 1 1 -1 -1 1 1 -1 -1 -1 -1; + 1 1 1 1 1 -1 -1 -1 -1 -1 +] + +E_letter = [ + 1 1 1 1 1 1 1 1 1 1; + 1 1 -1 -1 -1 -1 -1 -1 -1 -1; + 1 1 -1 -1 -1 -1 -1 -1 -1 -1; + 1 1 1 1 1 1 1 1 1 1; + 1 1 -1 -1 -1 -1 -1 -1 -1 -1; + 1 1 -1 -1 -1 -1 -1 -1 -1 -1; + 1 1 1 1 1 1 1 1 1 1; + 1 1 -1 -1 -1 -1 -1 -1 -1 -1; + 1 1 -1 -1 -1 -1 -1 -1 -1 -1; + 1 1 1 1 1 1 1 1 1 1 +] + +F = [ + 1 1 1 1 1 1 1 1 1 1; + 1 1 -1 -1 -1 -1 -1 -1 -1 -1; + 1 1 -1 -1 -1 -1 -1 -1 -1 -1; + 1 1 1 1 1 1 1 1 1 1; + 1 1 -1 -1 -1 -1 -1 -1 -1 -1; + 1 1 -1 -1 -1 -1 -1 -1 -1 -1; + 1 1 -1 -1 -1 -1 -1 -1 -1 -1; + 1 1 -1 -1 -1 -1 -1 -1 -1 -1; + 1 1 -1 -1 -1 -1 -1 -1 -1 -1; + 1 1 -1 -1 -1 -1 -1 -1 -1 -1 +] + +G = [ + 1 1 1 1 1 1 1 1 1 1; + 1 1 -1 -1 -1 -1 -1 -1 -1 -1; + 1 1 -1 -1 -1 -1 -1 -1 -1 -1; + 1 1 -1 -1 -1 -1 -1 -1 -1 -1; + 1 1 -1 -1 -1 -1 -1 1 1 1; + 1 1 -1 -1 -1 -1 -1 -1 1 1; + 1 1 -1 -1 -1 -1 -1 -1 1 1; + 1 1 -1 -1 -1 -1 -1 -1 1 1; + -1 1 1 1 1 1 1 1 1 1; + -1 -1 1 1 1 1 1 1 1 1 +] + +encodings = [A, B, C, D, E_letter, F, G] + +spin_numbers(row, col) = (row - 1) * 10 + col + +# create strength of interactions between ith and jth spinds in Ji,j +J = zeros(100, 100) +for m in 1:100 # row + for n in 1:100 # col + i = (m - 1) % 10 + 1 + j = (n - 1) % 10 + 1 + k = (m - 1) ÷ 10 + 1 + l = (n - 1) ÷ 10 + 1 + + J[m, n] = 0 + for encoding in encodings + J[m, n] += encoding[i, k] * encoding[j, l] + end + J[m, n] /= length(encodings) + end +end + +# the energy function +function energy(s, J) + E = 0 + for i in 1:100 + for j in 1:100 + E += J[i, j] * s[i] * s[j] + end + end + return -E +end + +# the update function (uses monte carlo steps) +function monte_carlo(s, J) + for i in 1:100 # systematically go through each point in the lattice + # calculate the energy of the system + E = energy(s, J) + # randomly flip a spin + s[i] = -s[i] + # calculate the new energy of the system + E_new = energy(s, J) + # calculate the change in energy + dE = E_new - E + # if the change in energy is positive, flip the spin back + if dE > 0 + s[i] = -s[i] + end + end + return s +end + +# create the main function +function main(s, J, nsteps) + E = energy(s, J) + for i in 1:nsteps + s = monte_carlo(s, J) + E = energy(s, J) + end + return s, E +end + +# run the main function + +# randomly change some values in A to see if NN works +function produce_test_arr(enc, prob_change = 0.1) + tmp = zeros(10, 10) + for i in 1:10 + for j in 1:10 + if rand() < prob_change + tmp[i, j] = rand([-1, 1]) # set some random values + else + tmp[i, j] = enc[i, j] + end + end + end + return tmp +end + +function check_if_same(s, enc) + for i in 1:10 + for j in 1:10 + if s[i, j] != enc[i, j] + println("The neural net did not work") + return false + end + end + end + # println("The neural net worked") + return true +end + +function apply_damage_to_J(J, prob_damage = 0.8) + println("Applying damage to J with probability ", prob_damage) + ret = copy(J) + for i in 1:100 + for j in 1:100 + if rand() < prob_damage + ret[i, j] = 0 + else + ret[i, j] = J[i, j] + end + end + end + return ret +end + +function run_tests(num_times, Js, MC_steps = 2) + num_correct = 0 + num_total = 0 + + for i in 1:num_times + # randomly select a memory + enc = encodings[rand(1:length(encodings))] + # produce a test array + s = produce_test_arr(enc) + # run theNN + s, E = main(s, Js, MC_steps) + # check if the NN worked + if check_if_same(s, enc) + num_correct += 1 + end + num_total += 1 + end + + println("Number of correct tests: ", num_correct) + + return num_correct / num_total +end + +function run_tests_for_damage_range(damages, num_times, MC_steps = 2, init_J = J) + res = [] + + for damage in damages + damaged_J = apply_damage_to_J(init_J, damage) + p = run_tests(num_times, damaged_J, MC_steps) + + push!(res, p) + end + + return res +end + +damage_range = collect(0.0:0.1:0.15) +res = run_tests_for_damage_range(damage_range, 10) + +#plot the results +using Plots +plot(damage_range, res, xlabel = "Damage", ylabel = "Success rate", title = "NN success rate vs damage probability", marker = :circle, label = "MC steps = 2") + +res = run_tests_for_damage_range(damage_range, 10, 10) +plot!(damage_range, res, xlabel = "Damage", ylabel = "Success rate", title = "NN success rate vs damage probability", marker = :circle, label = "MC steps = 10") + +res = run_tests_for_damage_range(damage_range, 10, 100) +plot!(damage_range, res, xlabel = "Damage", ylabel = "Success rate", title = "NN success rate vs damage probability", marker = :circle, label = "MC steps = 100") + +res = run_tests_for_damage_range(damage_range, 10, 1000) +plot!(damage_range, res, xlabel = "Damage", ylabel = "Success rate", title = "NN success rate vs damage probability", marker = :circle, label = "MC steps = 1000") + +savefig("12-12-all.png") \ No newline at end of file diff --git a/hw9/YaoQuantumComputing.jl b/hw9/YaoQuantumComputing.jl new file mode 100644 index 0000000..36c22d6 --- /dev/null +++ b/hw9/YaoQuantumComputing.jl @@ -0,0 +1,28 @@ +using Yao, YaoPlots + +let + circuit = chain(2, put(1=>X), put(2=>X)) + plot(circuit) +end + + +begin + bellcircuit = chain(2, put(1=>H), control(1, 2=>X)) + plot(bellcircuit) +end + + +q1 = ArrayReg(bit"00") #creating the system of two qubits with state |00> +println(state(q1)) +println(q1 |> r->measure(r, nshots=10)) + +a = (q1 |> bellcircuit) +println(state(a)) +println(a |> r->measure(r, nshots=10)) + + +reversebellcircuit = chain(2, control(1,2=>X), put(1=>H)) +plot(reversebellcircuit) + +b = (a |> reversebellcircuit) +println(state(b)) \ No newline at end of file diff --git a/hw9/qubit00.svg b/hw9/qubit00.svg new file mode 100644 index 0000000..677f7e7 --- /dev/null +++ b/hw9/qubit00.svg @@ -0,0 +1,39 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + \ No newline at end of file diff --git a/hw9/qubit01.svg b/hw9/qubit01.svg new file mode 100644 index 0000000..e24d9f4 --- /dev/null +++ b/hw9/qubit01.svg @@ -0,0 +1,39 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + \ No newline at end of file diff --git a/hw9/qubit10.svg b/hw9/qubit10.svg new file mode 100644 index 0000000..8dcead0 --- /dev/null +++ b/hw9/qubit10.svg @@ -0,0 +1,39 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + \ No newline at end of file diff --git a/hw9/qubit11.svg b/hw9/qubit11.svg new file mode 100644 index 0000000..0833dcc --- /dev/null +++ b/hw9/qubit11.svg @@ -0,0 +1,39 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + \ No newline at end of file diff --git a/hw9/yao-assigment-2.jl b/hw9/yao-assigment-2.jl new file mode 100644 index 0000000..7de7987 --- /dev/null +++ b/hw9/yao-assigment-2.jl @@ -0,0 +1,75 @@ +begin + using Pkg + Pkg.activate(mktempdir()) + Pkg.Registry.update() + Pkg.add("Yao") + Pkg.add("YaoPlots") + Pkg.add("StatsBase") + Pkg.add("Plots") + Pkg.add("BitBasis") +end + +using Yao, YaoPlots + +# make the bell circuit +bellcircuit = chain(2, put(1 => H), control(1, 2 => X)) + +# make the reverse bell circuit +reversebellcircuit = chain(2, control(1, 2 => X), put(1 => H)) + +# circuit that takes two qubits and passes it through +# the bell circuit then the reverse bell circuit +circuit = chain(2, bellcircuit, reversebellcircuit) +# plot(circuit) + +# make the qubits |00>, |01>, |10>, |11> +qubit00 = ArrayReg(bit"00") +qubit10 = ArrayReg(bit"10") # circuit reads in reverse order +qubit01 = ArrayReg(bit"01") +qubit11 = ArrayReg(bit"11") + +using StatsBase: Histogram, fit +using Plots: bar, scatter!, gr; +gr(); +using BitBasis +function plotmeasure(x::Array{BitStr{n, Int}, 1}) where n + hist = fit(Histogram, Int.(x), 0:2^n) + x = 0 + if (n <= 3) + s = 8 + elseif (n > 3 && n <= 6) + s = 5 + elseif (n > 6 && n <= 10) + s = 3.2 + elseif (n > 10 && n <= 15) + s = 2 + elseif (n > 15) + s = 1 + end + bar(hist.edges[1] .- 0.5, hist.weights, legend = :none, size = (600 * (2^n) / s, 400), ylims = (0:maximum(hist.weights)), xlims = (0:2^n), grid = :false, ticks = false, border = :none, color = :lightblue, lc = :lightblue) + scatter!(0:2^n-1, ones(2^n, 1), markersize = 0, + series_annotations = "|" .* string.(hist.edges[1]; base = 2, pad = n) .* "⟩") + scatter!(0:2^n-1, zeros(2^n, 1) .+ maximum(hist.weights), markersize = 0, + series_annotations = string.(hist.weights)) +end + +# pass them through the circuit 1024 times, taking measurements +println("Qubit 00") +measured_qubits00 = qubit00 |> circuit -> measure(circuit, nshots = 1024) +println(measured_qubits00) +plotmeasure(measured_qubits00) + +println("Qubit 01") +measured_qubits01 = qubit01 |> circuit -> measure(circuit, nshots = 1024) +println(measured_qubits01) +plotmeasure(measured_qubits01) + +println("Qubit 10") +measured_qubits10 = qubit10 |> circuit -> measure(circuit, nshots = 1024) +println(measured_qubits10) +plotmeasure(measured_qubits10) + +println("Qubit 11") +measured_qubits11 = qubit11 |> circuit -> measure(circuit, nshots = 1024) +println(measured_qubits11) +plotmeasure(measured_qubits11) -- cgit v1.2.3-70-g09d2