add Julia implementation for chp.c

README.md

@@ -17,3 +17,8 @@ This project will help scientists to design better architecture of quantum compu
 **This package is still under developing, installation is not suggested.**
 
 # Usage
+
+# TODO
+
+- [ ] QuIDD based simulation
+- [ ] reload `show`s

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(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,12 @@
 module QuCmp
 
-# package code goes here
+abstract QuComput
+abstract AbstractModels{N}
+
+include("const.jl")
+include("utils/LogicExpr.jl")
+include("utils/math.jl")
+include("Adiabatic/Adiabatic.jl")
+include("circuit/Circuit.jl")
 
 end # module

src/circuit/Circuit.jl

@@ -0,0 +1,12 @@
+using Expokit
+import Base: show,bin
+export Gate,Hadamard,X,Y,Z,Circuit,addgate!,HadamardUnit
+
+abstract QuCircuit{N}<:AbstractModels{N}
+
+include("state.jl")
+include("op.jl")
+include("gates.jl")
+include("qcircuit.jl")
+include("process.jl")
+include("chp.jl")

src/circuit/README.md

@@ -0,0 +1,72 @@
+# Operators
+
+## Matrix Operators
+
+
+## Functional Operators
+
+should be constructed by a function with single input as type `QuState`
+
+# Gate
+
+## Gate
+
+As interface or wrapper to normally used gates.
+
+## Gate Unit
+
+gate unit is a basic data type for storing gates, the difference between `AbstractGateUnit` and `Gate` is that `Gate` is only a wrapper for quantum operators in quantum information processing. But `AbstractGateUnit` and its subtypes provide a way to store the position and way of processing for a certain gate.
+
+# Circuits
+
+`QuCircuit{N}` is the super-type for all circuit objects, where `N` is the number of qubits involved in this circuit.
+
+## `Circuit{N}`
+
+`Circuit{N}` is the default circuit type. The processing simulation of this circuit type will not be optimized to specific algorithm, which is not recommended if you need high performance.
+
+## `stlzCircuit{N}` (**TODO**)
+
+`stlzCircuit{N}` is for stabilizer circuits, a kind of circuit with only Hadamard, C-NOT, Phase gates, and followed by one bit measurement only. The current plan for this part is the rework of [CHP.c](http://www.scottaaronson.com/chp/) by Aanronson.
+
+# process
+
+This part makes use of Julia's multiple dispatch. If you have a new type of gate unit. You will need to overload the process function if there are special optimized algorithms for these type of gates. eg.
+
+for stablizer circuits, `HadamardUnit`,`PhaseUnit`,`CNOTUnit` are implemented in this package. All of them is implemented with a single process function.
+
+```julia
+function process{N}(unit::CtrlGateUnit{N},input::AbstractSparseArray)
+function process(unit::HadamardUnit,input::AbstractSparseArray)
+# ...
+```
+
+# How to customize your own gate?
+
+If you want to create your own gate and its related simulation algorithm, the following funcitons should be overloaded:
+
+- **define your own gate unit type** :: define your own gate unit type and it should be a subtype of `AbstractGateUnit{N}`
+    - the gate unit should at least contain one element named `pos` for storing the gate's position in a circuit
+- `addgate!`:: provide a method to add gates to subtypes of `QuCircuit{N}`, where `N` is the number of qubits.
+- `process` :: provide a method to process this gate by customized simulation algorithm.
+
+# Roadmap
+
+
+- [x] Basic structure of types in quantum circuit and adiabatic model
+  - [x] C-NOT
+  - [x] Hadamard
+  - [x] Pauli-X, Pauli-Y, Pauli-Z
+  - [x] TimeOp
+  - [x] object: `Circuit`
+
+- Quantum circuits and adiabatic computation
+  - [ ] `addgate!`, `removegate!` -> working on
+  - [x] pHamiltonian
+  - [ ] merge GPU acceleration from AdiaComput.jl
+
+- Visualization
+  - [ ] `plot` (circuit) -> next
+
+- Error correction
+  - [ ] CSS code