Source [part 2]: Adiabatic Computing

src/Adiabatic/Adiabatic.jl

@@ -0,0 +1,10 @@
+################################################################################
+#
+# This file simulates an adiabatic computing
+#
+################################################################################
+
+
+include("Base.jl")
+include("Hamiltonian.jl")
+include("Operator.jl")

src/Adiabatic/Base.jl

@@ -0,0 +1,67 @@
+type AdiaComputer <: QuComput
+    ###########################################
+    #  Eternal
+    ###########################################
+    HB::AbstractMatrix      # Base Hamiltonian
+    HP::AbstractMatrix      # Problem Hamiltonian
+    maxtime::Real           # max evolution time
+    n::Int64           # number of bits
+    dt::Real                # time step
+
+    ###########################################
+    # variable for current state
+    ###########################################
+    location::Real          # time location
+    state::AbstractVector
+    eigens::AbstractVecOrMat
+
+    ###########################################
+    # measure
+    ###########################################
+
+    prob::Real
+
+    ###########################################
+    # other
+    ###########################################
+
+    nev::Int
+
+    # n is the number of bits
+    # maxtime is the max evolution time
+
+    function AdiaComputer{M,N}(pH::AbstractMatrix,n::Int,maxtime::Real;dt=1e-2,nev=6)
+        HP = pH
+        HB = bHamilton(n)
+
+        # initialize time location
+        location = 0
+        # prepare the initial state
+        state = convert(Array{Complex,1},[1/sqrt(2^n) for i=1:2^n])
+
+        # initialize states
+        eigens = eigs(HB,nev=nev,which=:SM)[1].'
+
+        # adjust nev if the number of bits is smaller than 3
+
+        if n<3
+            nev = 2^n
+            warn("adjusting nev to $(nev)\n")
+        end
+
+        # probility
+        prob = norm([x==findmin(diag(HP))[2]?1:0 for x=1:2^n])^2
+
+        new(HB,HP,maxtime,n,dt,location,state,eigens,prob,nev)
+      end
+end
+
+# 3-SAT problem
+
+AdiaComputer{M,N}(ins::Instance{M,N},n::Int,maxtime::Real;dt=1e-2,nev=6) = AdiaComputer{M,N}(pHamilton(ins,n),n,maxtime;dt=dt,nev=nev)
+
+function Hamiltonian(Hs::AdiaComputer)
+    return (1-Hs.location)*Hs.HB+Hs.location*Hs.HP
+end
+
+export AdiaComputer

src/Adiabatic/Cooling.jl

@@ -0,0 +1,69 @@
+function cooler!(Hs::AdiaComputer,gamma::Real,t::Real)
+    # boundscheck(Hs,gamma,t)
+    energy = eig(full(Hamiltonian(Hs)))[1]
+    @assert -pi/2<minimum(energy)*t-gamma<pi/2 "bad cooling parameters"
+    @assert -pi/2<maximum(energy)*t-gamma<pi/2 "bad cooling parameters"
+
+    Hs.state = normalize!((Hs.state-im*exp(im*gamma)*
+            trotter(-im*t*(1-Hs.location)*Hs.HB,
+                -im*t*Hs.location*Hs.HP,3)*Hs.state)/2)
+end
+
+function heater!(Hs::AdiaComputer,gamma::Real,t::Real)
+    # boundscheck(Hs,gamma,t)
+    energy = eig(full(Hamiltonian(Hs)))[1]
+    @assert -pi/2<minimum(energy)*t-gamma<pi/2 "bad cooling parameters"
+    @assert -pi/2<maximum(energy)*t-gamma<pi/2 "bad cooling parameters"
+
+    Hs.state = normalize!((Hs.state+im*exp(im*gamma)*
+                trotter(-im*t*(1-Hs.location)*Hs.HB,
+                    -im*t*Hs.location*Hs.HP,3)*Hs.state)/2)
+end
+
+function daemon!(
+    Hs::AdiaComputer,
+    gamma::Real,
+    t::Real
+    )
+    dice = rand()
+    if dice <= 0.5*(1+sin(gamma))
+        #get |1> (higher energy)
+        heater!(Hs,gamma,t)
+        return 0.5*(1-sin(gamma))
+    else
+        #get |0> (lower energy)
+        cooler!(Hs,gamma,t)
+        return 0.5*(1+sin(gamma))
+    end
+end
+
+function cooling!(
+    Hs::AdiaComputer;
+    n=5)
+
+    count = 0
+    gamma,t = CoolingPara(Hs)
+
+    while count<n
+        Hs.prob *= daemon!(Hs,gamma,t)
+        count += 1
+    end
+    return Hs
+end
+
+function CoolingPara(Hs::AdiaComputer)
+    Eigen = eig(full(Hamiltonian(Hs)))[1]
+    maxEigen = maximum(Eigen)
+    minEigen = minimum(Eigen)
+
+    gamma = (maxEigen+minEigen)/(maxEigen-minEigen) * pi/2 * 0.1
+
+    if (gamma-pi/2)>0
+        t = 0.5*((gamma-pi/2)/minEigen + (gamma+pi/2)/maxEigen)
+    else
+        t = 0.5*((gamma-pi/2)/maxEigen + (gamma+pi/2)/maxEigen)
+    end
+    return gamma,t
+end
+
+export cooling!

src/Adiabatic/Hamiltonian.jl

@@ -0,0 +1,32 @@
+# Base Hamiltonian
+function bHamilton(bitnum::Int)
+    H = sparse(0.5*(I-σ₁))
+    Iden = spdiagm([1,1])
+    for i=2:bitnum
+        H = H⊗Iden
+    end
+
+    for i=2:bitnum
+        H_BC = Iden
+        for j=2:bitnum
+            if i==j
+                H_BC = H_BC⊗sparse(0.5*(I-σ₁))
+            else
+                H_BC = H_BC⊗Iden
+            end
+        end
+        H+=H_BC
+    end
+    return H
+end
+
+# problem Hamiltonian
+function pHamilton{M,N}(ins::Instance{M,N},n::Integer)
+    sum = spzeros(2^n,2^n);
+    for clause in ins.c
+        sum += spdiagm(Int[clause(i) for i=0:2^n-1])
+    end
+    return sum
+end
+
+export bHamilton,pHamilton

src/Adiabatic/Operator.jl

@@ -0,0 +1,2 @@
+include("timeop.jl")
+include("Cooling.jl")

src/Adiabatic/README.md

@@ -0,0 +1,11 @@
+# Adiabatic
+
+This folder is for simulations on adiabatic computing.
+
+# Contents
+
+- **Adiabatic** pack adiabatic related package in this folder
+    - **Base** basic types for adiabatic computing
+    - **Hamiltonian** Hamiltonian generators for adiabatic computing
+    - **Cooling** Daemon-like Cooling accroding to [doi:10.103](http://www.nature.com/nphoton/journal/v8/n2/full/nphoton.2013.354.html)
+    - **Operator** quantum operators realted to adiabatic computing

src/Adiabatic/timeop.jl

@@ -0,0 +1,21 @@
+function realtimeop!(Hs::AdiaComputer)
+    Hs.state = trotter(-im*(1-Hs.location)*Hs.HB,-im*Hs.location*Hs.HP,3)*Hs.state
+    Hs.location += Hs.dt/Hs.maxtime
+end
+
+function next_timestep!(Hs::AdiaComputer;evopercentage::Real=1/3,nev=6)
+    @assert (0<=Hs.location+evopercentage)&&( Hs.location+evopercentage-1<0.1) "evolutoin percentage out of bounds(should be in [0,1])"
+
+    const evotime = Hs.maxtime
+    const dt      = Hs.dt
+
+    # calculate eigen values
+    for i=Hs.location*evotime:dt:(Hs.location+evopercentage)*evotime
+        realtimeop!(Hs)
+        Hs.eigens = [Hs.eigens;eigs(Hamiltonian(Hs),nev=nev,which=:SM)[1].']
+    end
+
+    return Hs
+end
+
+export next_timestep!

src/QuCmp.jl

@@ -1,5 +1,15 @@
 module QuCmp
 
+<<<<<<< HEAD
+abstract QuComput
+
+include("const.jl")
+include("utils/LogicExpr.jl")
+include("utils/math.jl")
+include("Adiabatic/Adiabatic.jl")
+
+=======
 # package code goes here
+>>>>>>> ad742772d3e82260d55f64ba87c31ee32eb94575
 
 end # module

src/const.jl

@@ -0,0 +1,10 @@
+const σ₀=I
+const σ₁=[0 1;1 0]
+const σ₂=[0 -im;im 0]
+const σ₃=[1 0;0 -1]
+const sigmax = [0 1;1 0]
+const sigmay = [0 -im;im 0]
+const sigmaz = [1 0;0 -1]
+const hadamard = [1/sqrt(2) 1/sqrt(2);1/sqrt(2) -1/sqrt(2)]
+
+export σ₀,σ₁,σ₂,σ₃,sigmax,sigmay,sigmaz,hadamard