SpiDySims/runs/Projections.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "f52c5daf-9123-4b36-b83c-d5d2bcae5d91",
   "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 = 1\n",
    "N = 200\n",
    "tspan = (0., N*Δt)\n",
    "saveat = (0:1:N)*Δt\n",
    "α = 0.16\n",
    "ω0 = 1.4\n",
    "Γ = 0.5\n",
    "J = LorentzianSD(α, ω0, Γ) # coloring the noise\n",
    "matrix = AnisoCoupling([-sin(π/4) 0. 0. # coupling to the environment\n",
    "                        0. 0. 0.\n",
    "                        cos(π/4) 0. 0.]);\n",
    "T = .01\n",
    "noise = ClassicalNoise(T);\n",
    "\n",
    "nspin = 1 # number of spins\n",
    "\n",
    "s0 = zeros(3*nspin)\n",
    "for i in 1:nspin\n",
    "    ϵ = 0.1\n",
    "    s0tmp = [ϵ*rand(), ϵ*rand(), -1]\n",
    "    s0[1+(i-1)*3:3+(i-1)*3] = s0tmp./norm(s0tmp)\n",
    "end\n",
    "J0 = 1.\n",
    "JH = Nchain(nspin, J0)\n",
    "\n",
    "nruns = 1000\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6e3cdf02-c243-4622-a52a-a6d7f987c68d",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Starting...\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\u001b[32mProgress:   7%|██▊                                      |  ETA: 0:13:09\u001b[39m"
     ]
    }
   ],
   "source": [
    "println(\"Starting...\")\n",
    "progress = Progress(nruns);\n",
    "\n",
    "sols = zeros(nruns, length(saveat), 3*nspin)\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[i, :, :] = transpose(sol[:, :])\n",
    "    next!(progress)\n",
    "end\n",
    "\n",
    "solavg = mean(sols, dims=1)[1, :, :];\n",
    "solstd = std(sols, dims=1)[1, :, :];"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d16b97dd-b9f5-415c-b388-2b1f3344be45",
   "metadata": {},
   "outputs": [],
   "source": [
    "plot(saveat, solavg[:, 1], ribbon=solstd[:, 1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a06addc4-47b2-4467-8148-3e68f3eae702",
   "metadata": {},
   "outputs": [],
   "source": [
    "plot(saveat, solavg[:, 2], ribbon=solstd[:, 2])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bc1686fd-d86d-42b0-9d87-bcb911526d00",
   "metadata": {},
   "outputs": [],
   "source": [
    "plot(saveat, solavg[:, 3], ribbon=solstd[:, 3])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "59ac92d9-6a78-402b-91ae-3f5fd22cc6ff",
   "metadata": {},
   "outputs": [],
   "source": [
    "projected = zeros(nruns, length(saveat), 2)\n",
    "\n",
    "normalized_avg = eachrow(solavg) |> x -> normalize(x)\n",
    "\n",
    "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[j, i, :]\n",
    "        proj = dot(b, n) * n\n",
    "        b_ort = b - proj\n",
    "\n",
    "        projected[j, i, 1] = dot(u,b)\n",
    "        projected[j, i, 2] = dot(v,b)\n",
    "    end\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "03467eca-3f0c-4bd4-b3f7-5649b68780bc",
   "metadata": {},
   "outputs": [],
   "source": [
    "histogram(projected[:, 197, 1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9c6484d3-3396-47bb-b063-98e8d013d020",
   "metadata": {},
   "outputs": [],
   "source": [
    "histogram(projected[:, 20, 2])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f6df949c-15b2-447f-a264-8ef48ce1725e",
   "metadata": {},
   "outputs": [],
   "source": [
    "function nth_moment(data, N)\n",
    "    return mean(data .^ N)\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e4ad51fc-0dc7-4106-99e4-2b29e54b1930",
   "metadata": {},
   "outputs": [],
   "source": [
    "using HypothesisTests"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e5b9609e-f467-4bd7-9b51-9a4c1f175738",
   "metadata": {},
   "outputs": [],
   "source": [
    "JB = zeros(N)\n",
    "for i in 10:N\n",
    "    JB[i] = HypothesisTests.JarqueBeraTest(projected[:, i, 1]).JB\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3ad3f882-caa1-4dfc-b83c-3b57c214ecea",
   "metadata": {},
   "outputs": [],
   "source": [
    "plot(10:N, JB[10:N])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bb3c6b84-5b21-4d79-b5bc-3345a8f2d893",
   "metadata": {},
   "outputs": [],
   "source": [
    "HypothesisTests.JarqueBeraTest(projected[:, 197, 1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e331f6c7-fb4b-4f21-84d3-6a3f203a95b6",
   "metadata": {},
   "outputs": [],
   "source": [
    "HypothesisTests.JarqueBeraTest(projected[:, 200, 1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "240660fb-d008-4dab-8bc1-aaf530885102",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 1.9.4",
   "language": "julia",
   "name": "julia-1.9"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.9.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}