{ "cells": [ { "cell_type": "code", "execution_count": 1, "id": "e47f1fd3-c09d-43cb-80c8-64f42e8161a4", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "1000" ] }, "execution_count": 1, "metadata": {}, "output_type": "execute_result" } ], "source": [ "using SpiDy\n", "using NPZ\n", "using DataFrames\n", "using CSV\n", "using ProgressMeter\n", "using Random\n", "using Statistics\n", "using LinearAlgebra\n", "using Plots\n", "\n", "########################\n", "########################\n", "\n", "Δt = .2\n", "N = 1000\n", "tspan = (0., N*Δt)\n", "saveat = (0:1:N)*Δt\n", "\n", "matrix = IsoCoupling(1)\n", "\n", "T = .01\n", "noise = ClassicalNoise(T);\n", "\n", "theta = pi * 3/8\n", "\n", "s0 = [0, 1, cos(theta)]\n", "s0 = s0 ./ norm(s0)\n", "\n", "J0 = 1.\n", "JH = Nchain(1, J0)\n", "\n", "nruns = 1000" ] }, { "cell_type": "code", "execution_count": 2, "id": "6e3cdf02-c243-4622-a52a-a6d7f987c68d", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Starting...\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "\u001b[32mProgress: 100%|█████████████████████████████████████████| Time: 0:03:04\u001b[39mm\n" ] } ], "source": [ "α = 10. # 0.16\n", "ω0 = 7. # 1.4\n", "Γ = 5. # 0.5\n", "\n", "J = LorentzianSD(α, ω0, Γ) # coloring the noise\n", "\n", "println(\"Starting...\")\n", "progress = Progress(nruns);\n", "\n", "sols_ohm = zeros(nruns, length(saveat), 3)\n", "\n", "Threads.@threads for i in 1:nruns\n", " bfields = [bfield(N, Δt, J, noise),\n", " bfield(N, Δt, J, noise),\n", " bfield(N, Δt, J, noise)];\n", " sol = diffeqsolver(s0, tspan, J, bfields, matrix; JH=JH, saveat=saveat);\n", " sols_ohm[i, :, :] = transpose(sol[:, :])\n", " next!(progress)\n", "end\n", "\n", "solavg_ohm = mean(sols_ohm, dims=1)[1, :, :];\n", "solstd_ohm = std(sols_ohm, dims=1)[1, :, :];" ] }, { "cell_type": "code", "execution_count": 16, "id": "3422af78-2f32-468c-a8a3-a0732a836fd1", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Starting...\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "\u001b[32mProgress: 100%|█████████████████████████████████████████| Time: 0:00:52\u001b[39m\n" ] } ], "source": [ "α = 0.16\n", "ω0 = 1.4\n", "Γ = 0.5\n", "\n", "J = LorentzianSD(α, ω0, Γ) # coloring the noise\n", "\n", "println(\"Starting...\")\n", "progress = Progress(nruns);\n", "\n", "sols_col = zeros(nruns, length(saveat), 3)\n", "\n", "Threads.@threads for i in 1:nruns\n", " bfields = [bfield(N, Δt, J, noise),\n", " bfield(N, Δt, J, noise),\n", " bfield(N, Δt, J, noise)];\n", " sol = diffeqsolver(s0, tspan, J, bfields, matrix; JH=JH, saveat=saveat);\n", " sols_col[i, :, :] = transpose(sol[:, :])\n", " next!(progress)\n", "end\n", "\n", "solavg_col = mean(sols_col, dims=1)[1, :, :];\n", "solstd_col = std(sols_col, dims=1)[1, :, :];" ] }, { "cell_type": "code", "execution_count": 27, "id": "e4ad51fc-0dc7-4106-99e4-2b29e54b1930", "metadata": {}, "outputs": [], "source": [ "using HypothesisTests\n", "using LaTeXStrings" ] }, { "cell_type": "code", "execution_count": 17, "id": "d16b97dd-b9f5-415c-b388-2b1f3344be45", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "\"/work/benstem/SpiDySims/runs/std_66_5.png\"" ] }, "execution_count": 17, "metadata": {}, "output_type": "execute_result" } ], "source": [ "plot(saveat, solstd_ohm[:, 1], label=L\"\\Delta x\", lw=1)\n", "plot!(saveat, solstd_ohm[:, 2], label=L\"\\Delta y\", lw=1)\n", "plot!(legend=:topright, xlabel=\"Time\", ylabel=\"Standard Deviation\")\n", "\n", "axis = twinx()\n", "plot!(axis, saveat, solavg_col[:, 3], label=L\"s_z\", ls=:dash)\n", "hline!(axis, [0], label=\"\", ls=:dash)\n", "plot!(axis, legend=:bottomright, ylabel=L\"$s_z$ Spin Component\")\n", "\n", "savefig(\"./std_66_5.png\")" ] }, { "cell_type": "code", "execution_count": 18, "id": "a06addc4-47b2-4467-8148-3e68f3eae702", "metadata": {}, "outputs": [ { "data": { "image/svg+xml": [ "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n" ] }, "execution_count": 18, "metadata": {}, "output_type": "execute_result" } ], "source": [ "plot(saveat, solstd_col[:, 1], label=L\"\\Delta s_x\")\n", "plot!(saveat, solstd_col[:, 2], label=L\"\\Delta s_y\")\n", "\n", "axis = twinx()\n", "plot!(axis, saveat, solavg_col[:, 3], label=L\"s_z\", ls=:dash)" ] }, { "cell_type": "code", "execution_count": 19, "id": "bc1686fd-d86d-42b0-9d87-bcb911526d00", "metadata": {}, "outputs": [ { "data": { "image/svg+xml": [ "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n" ] }, "execution_count": 19, "metadata": {}, "output_type": "execute_result" } ], "source": [ "plot(saveat, solavg_ohm[:, 3])\n", "plot!(saveat, solavg_col[:, 3])" ] }, { "cell_type": "code", "execution_count": 20, "id": "59ac92d9-6a78-402b-91ae-3f5fd22cc6ff", "metadata": {}, "outputs": [], "source": [ "projected_ohm = zeros(nruns, length(saveat), 2)\n", "\n", "normalized_avg = [normalize(vec(row)) for row in eachrow(solavg_ohm)]\n", "\n", "Threads.@threads for i in 1:length(saveat)\n", " n = normalized_avg[i, :][1]\n", "\n", " u = normalize(cross(n, [0,0,1]))\n", " v = cross(u, n)\n", "\n", " for j in 1:nruns\n", " b = sols_ohm[j, i, :]\n", " proj = dot(b, n) * n\n", " b_ort = b - proj\n", "\n", " projected_ohm[j, i, 1] = dot(u,b)\n", " projected_ohm[j, i, 2] = dot(v,b)\n", " end\n", "end" ] }, { "cell_type": "code", "execution_count": 21, "id": "f1661b24-42b1-4f59-989a-54dc3ed03bf8", "metadata": {}, "outputs": [], "source": [ "projected_col = zeros(nruns, length(saveat), 2)\n", "\n", "normalized_avg = [normalize(vec(row)) for row in eachrow(solavg_col)]\n", "\n", "Threads.@threads for i in 1:length(saveat)\n", " n = normalized_avg[i, :][1]\n", "\n", " u = normalize(cross(n, [0,0,1]))\n", " v = cross(u, n)\n", "\n", " for j in 1:nruns\n", " b = sols_col[j, i, :]\n", " proj = dot(b, n) * n\n", " b_ort = b - proj\n", "\n", " projected_col[j, i, 1] = dot(u,b)\n", " projected_col[j, i, 2] = dot(v,b)\n", " end\n", "end" ] }, { "cell_type": "code", "execution_count": 28, "id": "e5b9609e-f467-4bd7-9b51-9a4c1f175738", "metadata": {}, "outputs": [], "source": [ "JB_ohm = zeros(N)\n", "for i in 1:N\n", " JB_ohm[i] = HypothesisTests.JarqueBeraTest(projected_ohm[:, i, 1]).JB\n", "end" ] }, { "cell_type": "code", "execution_count": 29, "id": "51ccaeb8-2b3d-4630-afc3-2640142c4d54", "metadata": {}, "outputs": [], "source": [ "JB_col = zeros(N)\n", "for i in 1:N\n", " JB_col[i] = HypothesisTests.JarqueBeraTest(projected_col[:, i, 1]).JB\n", "end" ] }, { "cell_type": "code", "execution_count": 30, "id": "6c63a170-6eff-4c8b-8b80-542eb705666f", "metadata": {}, "outputs": [ { "data": { "image/svg+xml": [ "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n" ] }, "execution_count": 30, "metadata": {}, "output_type": "execute_result" } ], "source": [ "plot(saveat, JB_ohm, title=\"JB Test for Ohmic and Non-Ohmic Coupling\", label=\"Ohmic\", linewidth=2, xlabel=\"Time\", ylabel=\"JB Test\", legend=:topleft, ylimits=(0, 200))\n", "plot!(saveat, JB_col, label=\"Non-Ohmic\", linewidth=2)" ] }, { "cell_type": "code", "execution_count": 25, "id": "ad2674ba-9711-4feb-96ab-a39d88b27ce5", "metadata": {}, "outputs": [ { "ename": "LoadError", "evalue": "BoundsError: attempt to access 1000-element Vector{Float64} at index [1:1001]", "output_type": "error", "traceback": [ "BoundsError: attempt to access 1000-element Vector{Float64} at index [1:1001]", "", "Stacktrace:", " [1] throw_boundserror(A::Vector{Float64}, I::Tuple{UnitRange{Int64}})", " @ Base ./abstractarray.jl:691", " [2] checkbounds", " @ ./abstractarray.jl:656 [inlined]", " [3] getindex(A::Vector{Float64}, I::UnitRange{Int64})", " @ Base ./array.jl:867", " [4] gr_draw_segments(series::Plots.Series, x::StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}, y::Vector{Float64}, z::Nothing, fillrange::Nothing, clims::Tuple{Float64, Float64})", " @ Plots ~/.julia/packages/Plots/ju9dp/src/backends/gr.jl:1825", " [5] gr_add_series(sp::Plots.Subplot{Plots.GRBackend}, series::Plots.Series)", " @ Plots ~/.julia/packages/Plots/ju9dp/src/backends/gr.jl:1742", " [6] gr_display(sp::Plots.Subplot{Plots.GRBackend}, w::Measures.AbsoluteLength, h::Measures.AbsoluteLength, vp_canvas::Plots.GRViewport{Float64})", " @ Plots ~/.julia/packages/Plots/ju9dp/src/backends/gr.jl:964", " [7] (::Plots.var\"#537#538\"{Int64, Int64, Plots.GRViewport{Float64}})(sp::Plots.Subplot{Plots.GRBackend})", " @ Plots ~/.julia/packages/Plots/ju9dp/src/backends/gr.jl:688", " [8] foreach(f::Plots.var\"#537#538\"{Int64, Int64, Plots.GRViewport{Float64}}, itr::Vector{Plots.Subplot})", " @ Base ./abstractarray.jl:2712", " [9] gr_display(plt::Plots.Plot{Plots.GRBackend}, dpi_factor::Float64)", " @ Plots ~/.julia/packages/Plots/ju9dp/src/backends/gr.jl:688", " [10] #578", " @ ~/.julia/packages/Plots/ju9dp/src/backends/gr.jl:2062 [inlined]", " [11] withenv(::Plots.var\"#578#579\"{Plots.Plot{Plots.GRBackend}, Float64}, ::Pair{String, String}, ::Vararg{Pair{String, String}})", " @ Base ./env.jl:172", " [12] _show(io::IOStream, #unused#::MIME{Symbol(\"image/png\")}, plt::Plots.Plot{Plots.GRBackend})", " @ Plots ~/.julia/packages/Plots/ju9dp/src/backends/gr.jl:2057", " [13] #invokelatest#2", " @ ./essentials.jl:716 [inlined]", " [14] invokelatest", " @ ./essentials.jl:714 [inlined]", " [15] show", " @ ~/.julia/packages/Plots/ju9dp/src/output.jl:232 [inlined]", " [16] #345", " @ ~/.julia/packages/Plots/ju9dp/src/output.jl:6 [inlined]", " [17] open(::Plots.var\"#345#346\"{Plots.Plot{Plots.GRBackend}}, ::String, ::Vararg{String}; kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})", " @ Base ./io.jl:330", " [18] open", " @ ./io.jl:328 [inlined]", " [19] png(plt::Plots.Plot{Plots.GRBackend}, fn::String)", " @ Plots ~/.julia/packages/Plots/ju9dp/src/output.jl:6", " [20] savefig(plt::Plots.Plot{Plots.GRBackend}, fn::String)", " @ Plots ~/.julia/packages/Plots/ju9dp/src/output.jl:149", " [21] savefig(fn::String)", " @ Plots ~/.julia/packages/Plots/ju9dp/src/output.jl:154", " [22] top-level scope", " @ In[25]:6" ] } ], "source": [ "right_axis = twinx()\n", "\n", "plot!(right_axis, saveat, solavg_ohm[:, 3], label=\"Ohmic\", ls=:dash, ylabel=L\"S_z\", ylimits=(-1.5, 1.5))\n", "plot!(right_axis, saveat, solavg_col[:, 3], label=\"Non-Ohmic\", ls=:dash)\n", "\n", "savefig(\"./67_5.png\")" ] }, { "cell_type": "code", "execution_count": 26, "id": "a1c77edf-652b-42ea-8bc0-445749a22a2c", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "Jarque-Bera normality test\n", "--------------------------\n", "Population details:\n", " parameter of interest: skewness and kurtosis\n", " value under h_0: \"0 and 3\"\n", " point estimate: \"-0.09472173866793361 and 3.1398692030437387\"\n", "\n", "Test summary:\n", " outcome with 95% confidence: fail to reject h_0\n", " one-sided p-value: 0.3150\n", "\n", "Details:\n", " number of observations: 1000\n", " JB statistic: 2.31051\n" ] }, "execution_count": 26, "metadata": {}, "output_type": "execute_result" } ], "source": [ "HypothesisTests.JarqueBeraTest(projected_col[:, 125, 1])" ] }, { "cell_type": "code", "execution_count": null, "id": "48b2963f-d070-46c9-a077-fe4aad09155a", "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Julia 1.7.3", "language": "julia", "name": "julia-1.7" }, "language_info": { "file_extension": ".jl", "mimetype": "application/julia", "name": "julia", "version": "1.7.3" } }, "nbformat": 4, "nbformat_minor": 5 }