Add QuBase.jl and QuDynamics.jl as Dependencies

LICENSE.md

@@ -1,6 +1,6 @@
 The QuCmp.jl package is licensed under the MIT "Expat" License:
 
-> Copyright (c) 2016: Roger-Luo.
+> Copyright (c) 2016: Contributors to QuCmp.jl.
 >
 > Permission is hereby granted, free of charge, to any person obtaining
 > a copy of this software and associated documentation files (the

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/AQCShrodingerEq.jl

@@ -0,0 +1,82 @@
+import QuDynamics: QuStateEvolution,operator,propagate
+
+function eigenvector(index::Integer,n::Integer)
+    res = 1
+    for i in bin(index,n)
+        if i=='0'
+            res = kron(res,1/sqrt(2)*[1,1])
+        elseif i=='1'
+            res = kron(res,1/sqrt(2)*[1,-1])
+        end
+    end
+    return res
+end
+
+immutable AQCShrodingerEq{H<:QuBase.AbstractQuMatrix} <: QuEquation{1}
+    HB::H
+    HP::H
+    T::Real
+
+    function AQCShrodingerEq(HP::H,maxtime::Real=1)
+        n = size(HP)[1]|>log2|>Int
+        HB = n|>bHamilton
+
+        # Similar Matrix
+        P = eigenvector(0,n)
+        for i=1:2^n-1
+            P = [P eigenvector(i,n)]
+        end
+        invP = inv(P)
+
+        HB = QuArray(P*HB*invP,(ComputBasis(n),ComputBasis(n)))
+
+        # @show HP
+        # @show size(HP)
+        new(HB,HP,maxtime)
+    end
+end
+
+AQCShrodingerEq{H<:QuBase.AbstractQuMatrix}(HP::H,maxtime::Real=1) = AQCShrodingerEq{H}(HP,maxtime)
+
+QuStateEvolution{QPM<:QuDynamics.QuPropagatorMethod, QV<:QuBase.AbstractQuVector}(eq::AQCShrodingerEq, init_state::QV, tlist, method::QPM) = QuStateEvolution{QPM,QV,AQCShrodingerEq}(eq, init_state, tlist, method)
+
+function operator(eq::AQCShrodingerEq,current_t::Real)
+    return (1-current_t/eq.T)*eq.HB+current_t/eq.T*eq.HP
+end
+
+# For ODE solvers
+
+for (qu_ode_type,ode_solver) in QuDynamics.type_to_method_ode
+    @eval begin
+        function propagate(prob::$qu_ode_type, eq::AQCShrodingerEq, t, current_t, current_qustate)
+            op = operator(eq,current_t)
+            dims = size(current_qustate)
+            # Convert the current_qustate to complex as it might result in a Inexact Error. After complex is in QuBase.jl (PR #38)
+            # we could just do a complex(vec(current_qustate)) avoiding the coeffs(coeffs(vec(current_qustate))).
+            next_state = $ode_solver((t,y)-> -im*coeffs(op)*y, complex(coeffs(vec(current_qustate))), [current_t, t], points=:specified,
+                                  reltol = get(prob.options, :reltol, 1.0e-5), abstol = get(prob.options, :abstol, 1.0e-8))[2][end]
+            CQST = QuBase.similar_type(current_qustate)
+            return CQST(reshape(next_state, dims), bases(current_qustate))
+        end
+    end
+end
+
+# For Expm solvers
+
+function propagate(prob::QuExpokit, eq::AQCShrodingerEq, t, current_t, current_qustate)
+    dt = t - current_t
+    dims = size(current_qustate)
+    next_state = Expokit.expmv(dt, -im*coeffs(operator(eq,current_t)), coeffs(vec(current_qustate)), m = get(prob.options, :m, 30), tol = get(prob.options, :tol, 1e-7))
+    CQST = QuBase.similar_type(current_qustate)
+    return CQST(reshape(next_state, dims), bases(current_qustate))
+end
+
+function propagate(prob::QuExpmV, eq::AQCShrodingerEq, t, current_t, current_qustate)
+    dt = t - current_t
+    dims = size(current_qustate)
+    # @show coeffs(operator(eq,current_t))
+    next_state = ExpmV.expmv(im*dt, -coeffs(operator(eq,current_t)), coeffs(vec(current_qustate)), M = get(prob.options, :M, []), prec = get(prob.options, :prec, "double"),
+                            shift = get(prob.options, :shift, false), full_term = get(prob.options, :full_term, false))
+    CQST = QuBase.similar_type(current_qustate)
+    return CQST(reshape(next_state, dims), bases(current_qustate))
+end

src/Adiabatic/Adiabatic.jl

@@ -0,0 +1,46 @@
+################################################################################
+#
+# This file simulates Adiabatic Quantum Computing
+#
+################################################################################
+
+export AQC
+
+include("Hamiltonian.jl")
+include("AQCShrodingerEq.jl")
+
+typealias AbstractQuVecOrMat Union{QuBase.AbstractQuVector,QuBase.AbstractQuMatrix}
+
+type AQC{N,H<:QuBase.AbstractQuMatrix}<:AbstractQC{N}
+    eq::AQCShrodingerEq
+    state::AbstractQuVecOrMat
+    tlist
+    method
+
+    function AQC(HP::H,method,dt,maxtime)
+        eq = AQCShrodingerEq(HP,maxtime)
+        tlist = 0.:dt:maxtime
+        state = QuArray(convert(Array{Complex128,1},[1/sqrt(2^N) for i=1:2^N]),ComputBasis(N))
+        new(eq,state,tlist,method)
+    end
+end
+
+AQC{H<:QuBase.AbstractQuMatrix}(HP::H,n::Int;method = QuODE45(),dt = 1e-2, maxtime = 1.0) = AQC{n,H}(HP,method,dt,maxtime)
+
+
+function Base.start(aqc::AQC)
+    init_state = aqc.state
+    t_state = start(aqc.tlist)
+    t,t_state = next(aqc.tlist,t_state)
+    return QuPropagatorState(init_state,t,t_state)
+end
+
+function Base.next(aqc::AQC, qustate::QuPropagatorState)
+    current_qustate = qustate.state
+    current_t = qustate.t
+    t,t_state = next(aqc.tlist, qustate.t_state)
+    next_qustate = aqc.state = propagate(aqc.method, aqc.eq, t, current_t, current_qustate)
+    return (t, next_qustate), QuPropagatorState(next_qustate, t, t_state)
+end
+
+Base.done(aqc::AQC, qustate::QuPropagatorState) = done(aqc.tlist, qustate.t_state)

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,20 @@
+# Base Hamiltonian
+function bHamilton(bitnum::Int)
+    diagv = Array(Complex128,2^bitnum);
+    index = 1;
+    for i=0:bitnum
+        for j = 1:binomial(bitnum,i)
+            diagv[index] = i
+            index+=1
+        end
+    end
+
+    return spdiagm(diagv)
+end
+
+# problem Hamiltonian for SAT problem
+function pHamilton{M,N}(ins::Instance{M,N},n::Integer)
+    return spdiagm(Complex128[1-ins(Bits(i)) for i=0:2^n-1])
+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/QuCmp.jl

@@ -1,5 +1,18 @@
 module QuCmp
 
-# package code goes here
+using QuBase,QuDynamics
+
+export ComputBasis
+
+abstract AbstractQC{N}
+
+ComputBasis(n::Int) = FiniteBasis(ntuple(x->2,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}<:AbstractQC{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