diff --git a/README.md b/README.md index 8b2411f..f500383 100644 --- a/README.md +++ b/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 diff --git a/src/Adiabatic/Adiabatic.jl b/src/Adiabatic/Adiabatic.jl new file mode 100644 index 0000000..742066f --- /dev/null +++ b/src/Adiabatic/Adiabatic.jl @@ -0,0 +1,10 @@ +################################################################################ +# +# This file simulates an adiabatic computing +# +################################################################################ + + +include("Base.jl") +include("Hamiltonian.jl") +include("Operator.jl") diff --git a/src/Adiabatic/Base.jl b/src/Adiabatic/Base.jl new file mode 100644 index 0000000..bf11256 --- /dev/null +++ b/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 diff --git a/src/Adiabatic/Cooling.jl b/src/Adiabatic/Cooling.jl new file mode 100644 index 0000000..f155401 --- /dev/null +++ b/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 (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 count0 + 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! diff --git a/src/Adiabatic/Hamiltonian.jl b/src/Adiabatic/Hamiltonian.jl new file mode 100644 index 0000000..fe7e3db --- /dev/null +++ b/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 diff --git a/src/Adiabatic/Operator.jl b/src/Adiabatic/Operator.jl new file mode 100644 index 0000000..10bbf16 --- /dev/null +++ b/src/Adiabatic/Operator.jl @@ -0,0 +1,2 @@ +include("timeop.jl") +include("Cooling.jl") diff --git a/src/Adiabatic/README.md b/src/Adiabatic/README.md new file mode 100644 index 0000000..6d9a8c4 --- /dev/null +++ b/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 diff --git a/src/Adiabatic/timeop.jl b/src/Adiabatic/timeop.jl new file mode 100644 index 0000000..ac87450 --- /dev/null +++ b/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! diff --git a/src/QuCmp.jl b/src/QuCmp.jl index 5afc76d..dc7398e 100644 --- a/src/QuCmp.jl +++ b/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 diff --git a/src/circuit/Circuit.jl b/src/circuit/Circuit.jl new file mode 100644 index 0000000..d6bc8c2 --- /dev/null +++ b/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") diff --git a/src/circuit/README.md b/src/circuit/README.md new file mode 100644 index 0000000..e0ea2ba --- /dev/null +++ b/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 diff --git a/src/circuit/chp.jl b/src/circuit/chp.jl new file mode 100644 index 0000000..93349df --- /dev/null +++ b/src/circuit/chp.jl @@ -0,0 +1,74 @@ +# These functions are based on C functions in the chp.c program by Scott Aaronson +# These are pure Julia implementations, and do not link to the chp.c library +# chp.c is Copyright (c), Scott Aaronson, and is not allowed for commercial use + + + +########################### +# Stablizer Circuit +# Circuits with only Hadamard, C-NOT, R (Phase) +# followed by only one bits measurement +########################### + + +export stlzCircuit,HadamardUnit,CNOTUnit,PhaseUnit + +type stlzCircuit{N} <: QuCircuit{N} + gates::Array{GateUnit,1} +end + +stlzCircuit(num::Integer,gates::Array{GateUnit,1}) = stlzCircuit{num}(gates) +stlzCircuit(num::Integer) = stlzCircuit{num}(Array(GateUnit,0)) + +type HadamardUnit<:AbstractGateUnit{1} + pos::Int +end + +function process!{N}(unit::HadamardUnit,input::stlzState{N}) + a = unit.pos + # swap x_{ia} with z_{ia} + vec_temp = Array(Int,Int(2*N)) + vec_temp[:] = input.X[:,a] + input.X[:,a] = input.Z[:,a] + input.Z[:,a] = vec_temp[:] + + for i = 1:2N + input.R[i] $= input.X[i,a]&input.Z[i,a] + end + + return input +end + +type CNOTUnit<:AbstractGateUnit{1} + pos::Int + ctrl::Int +end + +function process!{N}(unit::CNOTUnit,input::stlzState{N}) + a = unit.pos + b = unit.ctrl + + r = input.R + x = input.X + z = input.Z + + for i=1:2N + r[i] $= x[i,a] & z[i,b] & (x[i,b]$z[i,a]$1) + x[i,b] $= x[i,a] + z[i,a] $= z[i,b] + end + return input +end + +type PhaseUnit<:AbstractGateUnit{1} + pos::Int +end + +function process!{N}(unit::PhaseUnit,input::stlzState{N}) + a = unit.pos + for i=1:2N + input.R[i] $= input.X[i,a] & input.Z[i,a] + input.Z[i,a] $= input.X[i,a] + end + return input +end diff --git a/src/circuit/design.md b/src/circuit/design.md new file mode 100644 index 0000000..3e6e9c1 --- /dev/null +++ b/src/circuit/design.md @@ -0,0 +1,22 @@ +# Matrix Operators + + +# Functional Operators + +should be constructed by a function with single input as type `QuState` + +# 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) +# ... +``` + +# QuIDD + +QuIDD methods needs to be done to store the state efficiently. diff --git a/src/circuit/gates.jl b/src/circuit/gates.jl new file mode 100644 index 0000000..2160c1b --- /dev/null +++ b/src/circuit/gates.jl @@ -0,0 +1,75 @@ +########################## +# Gates +########################## + +type Gate{T,N} + name::AbstractString + op::AbstractOp{T,N} + Gate{T,N}(name::AbstractString,op::AbstractOp{T,N}) = new(name,op) +end + +Gate{N}(name::AbstractString,op::MatrixOp{N})=Gate{AbstractMatrix,N}(name,op) + +Hadamard = Gate("Hadamard gate",OP_Hadamard) +X = Gate("Paili X",OP_sigmax) +Y = Gate("Paili Y",OP_sigmay) +Z = Gate("Paili Z",OP_sigmaz) + +function call{T,N}(gate::Gate{T,N}) + return gate.op() +end + + +""" +An AbstractGateUnit should at least contain one element called `pos`, +for definition of `pos`, please refer to the docs of `GateUnit` +""" +abstract AbstractGateUnit{N} + + +""" +GateUnit +--- + +GateUnit is for the gate unit in a circuit + +pos: +the first number in pos is the the column number +the second number in pos is the bits' IDs that is realted to this gate +eg. + +''' + --block 1--|----block 2----| +1 -----------|---------------| +2 ----[A]----|---------------| +3 -----------|---------------| +4 -----------|-----[ ]-----| +5 -----------|-----[ B ]-----| +6 -----------|-----[ ]-----| +''' + +The position of gate B is (2,4,5,6) +""" + +type GateUnit{T,N}<:AbstractGateUnit{N} + gate::Gate{T,N} + pos::Vector{Int} + # TODO : time layer? + # time_layer::Real + + GateUnit{T,N}(gate::Gate{T,N},pos::Vector{Int})=new(gate,pos) +end + +GateUnit{T,N}(gate::Gate{T,N},pos::Vector{Int}) = GateUnit{T,N}(gate,pos) +GateUnit{T,N}(gate::Gate{T,N},pos::Tuple{Int}) = GateUnit{T,N}(gate,collect(pos)) +GateUnit{T,N}(gate::Gate{T,N},pos::Int...) = GateUnit{T,N}(gate,collect(pos)) + +type CtrlGateUnit{T,N}<:AbstractGateUnit{N} + gate::Gate{T,N} + pos::Vector{Int} + ctrl::Vector{Int} +end + +# CtrlGateUnit{T,N}(gate::Gate{T,N},pos::Vector{Int},ctrl::Vector{Int}) = CtrlGateUnit{T,N}(gate,pos,ctrl) + +bitnum{N}(unit::AbstractGateUnit{N}) = N diff --git a/src/circuit/op.jl b/src/circuit/op.jl new file mode 100644 index 0000000..75b26f8 --- /dev/null +++ b/src/circuit/op.jl @@ -0,0 +1,77 @@ +abstract AbstractOp{T,N} + +################################## +# Matrix Operators +################################## + +type MatrixOp{N}<:AbstractOp{AbstractMatrix,N} + name::AbstractString + mat::AbstractMatrix +end + +MatrixOp(num::Integer,name::AbstractString,mat::AbstractMatrix) = MatrixOp{num}(name,mat) + +function call{N}(matop::MatrixOp{N}) + return [matop.mat[:,i] for i=1:size(matop.mat,2)] +end + +function call{N}(matop::MatrixOp{N},state::AbstractVector) + return matop.mat*state +end + +function show{N}(io::IO,matop::MatrixOp{N}) + println("$N bits matrix operator $(matop.name):") + show(matop.mat) +end + + +OP_Hadamard = MatrixOp(1,"Hadamard",hadamard) +# Pauli Groups +OP_sigmax = MatrixOp(1,"Pauli Sigma X",σ₁) +OP_sigmay = MatrixOp(1,"Pauli Sigma Y",σ₂) +OP_sigmaz = MatrixOp(1,"Pauli Sigma Z",σ₃) + +# TODO +# this comes from linalg/uniformscaling.jl +# same operators should be overloaded +immutable IdentityOp{T<:Number}<:AbstractOp{T} + λ::T +end + +OP_I = IdentityOp(1) + +################################## +# Function Operators +################################## + +const FUNCTION_OP_PARA_INF = -1 + +""" +Function operators should be able to accept AbstractSparseVector inputs + +The member `f` should be an one input only function +""" +type FunctionOp{N}<:AbstractOp{Function,N} + name::AbstractString + f::Function +end + +function basis(n::Int) + return [i->sparsevec(Dict[i=>1]) for i=1:2^n] +end + +function call{N}(funcop::FunctionOp{N}) + res = Array(Array{Complex,1},0) + for i in basis(N) + push!(res,funcop.f(i)) + end + return res +end + +function call{N}(funcop::FunctionOp{N},state::AbstractVector) + return funcop(state) +end + +function TimeOp{N}(state::QuState{N};Hamiltonian=I,dt=1e-6) + return expmv(-im*dt,Hamiltonian,state.s) +end diff --git a/src/circuit/process.jl b/src/circuit/process.jl new file mode 100644 index 0000000..45f0822 --- /dev/null +++ b/src/circuit/process.jl @@ -0,0 +1,70 @@ +function bin(x::Integer,pad::Int,id::Int) + return (x>>(pad-id))&1 +end + +function bin(x::Integer,pad::Int,ids::Vector{Int}) + ret = 0 + subid = 1;subpad = length(ids) + for id in ids + ret+=(bin(x,pad,id)<<(subpad-subid)) + subid+=1 + end + return ret +end + +function flip(x::Integer,pad::Int,id::Int) + return x$(1<<(pad-id)) +end + +function flip(x::Integer,pad::Int,ids::Vector{Int},assign::Int) + ret = x + n = length(ids) + for id = 1:n + if bin(assign,n,id)==1 + ret = flip(ret,pad,id) + end + end + return ret +end + +function set(x::Integer,pad::Int,ids::Vector{Int},assign::Integer) + ret = x + n = length(ids) + for i = 1:n + if bin(assign,n,i) != bin(x,n,ids[i]) + ret = flip(ret,pad,ids[i]) + end + end + return ret +end + +function process!{N}(unit::AbstractGateUnit{N},input::AbstractSparseArray) + len = length(unit.pos)-1 + eigens = unit.gate() + ret = spzeros(Complex,size(input)...) + for i in rowvals(input) + subidx = bin(i-1,len,unit.pos[2:end]) + rett = spzeros(Complex,size(input)...) + for j = 0:(2^len-1) + idx = set(i-1,N,unit.pos[2:end],j) + rett[idx+1] = input[i]*eigens[subidx+1][j+1] + end + ret += rett + end + input[:] = ret + return ret +end + +function process!{N}(circuit::Circuit{N},input::AbstractSparseArray) + ret = deepcopy(input) + for unit in circuit.gates + if N==bitnum(unit) + ret = unit.gate(ret) + else + ret = process!(unit,ret) + end + end + return ret +end + +export process! diff --git a/src/circuit/qcircuit.jl b/src/circuit/qcircuit.jl new file mode 100644 index 0000000..38d3e13 --- /dev/null +++ b/src/circuit/qcircuit.jl @@ -0,0 +1,57 @@ +############################ +# Circuits +############################ + +type Circuit{N} <: QuCircuit{N} + gates::Array{GateUnit,1} +end + +Circuit(num::Integer,gates::Array{GateUnit,1}) = Circuit{num}(gates) +Circuit(num::Integer)=Circuit{num}(Array(GateUnit,0)) + +# TODO +# function show(io::IO,circuit::Circuit) +# end + + +############################ +# Circuit Constructor +############################ + +max(a::Int) = a + +# NOTE +# Should we open AbstractGateUnit as API for users? + +function addgate!{N,M}(circuit::QuCircuit{N},gate::AbstractGateUnit{M}) + push!(circuit.gates,gate) + sort!(circuit,gates,alg=QuickSort,by=x->x.pos[1]) +end + +# For gate units +function addgate!{T,N,M}(circuit::QuCircuit{N},gate::Gate{T,M},pos::Vector{Int}) + # Bounds check + @assert length(pos)==M+1 "number of qubits do not match" + @assert maximum(pos[2:end])<=N "bit's id out of range, maximum is $N" + + addgate!(circuit,GateUnit(gate,pos)) +end + +addgate!{T,N,M}(circuit::QuCircuit{N},gate::Gate{T,M},pos::Int...) = addgate!(circuit,gate,collect(pos)) + +# for controled gates +function addgate!{T,N,M}(circuit::QuCircuit{N},gate::Gate{T,M},pos::Vector{Int},ctrl::Vector{Int}) + # Bounds check + @assert length(pos)+length(ctrl) == M+1 + @assert max(maximum(pos[2:end]),maximum(ctrl)) <= N "bit's id out of range, maximum is $N" + + addgate!(circuit,CtrlGateUnit(gate,pos,ctrl)) +end + +function removegate!{N}(circuit::QuCircuit{N},pos::Vector{Int}) + for i in eachindnex(circuits.gates) + if circuits.gates.pos == pos + deleteat!(circuits.gates,i) + end + end +end diff --git a/src/circuit/state.jl b/src/circuit/state.jl new file mode 100644 index 0000000..3d37475 --- /dev/null +++ b/src/circuit/state.jl @@ -0,0 +1,58 @@ + +import Base: +,-,*,/,^ +export stlzState + +# TODO: QuIDD + +# NOTE +# N is the number of bits in a quantum array + +abstract AbstractQuArray{N} + +# TODO +# QuIDD encoding array +# type QuIDDArray{T,N}<:AbstractQuArray{N} +# end + +type QuState{T,N}<:AbstractQuArray{N} + s::AbstractVector{T} +end + +(+)(A::QuState,B::QuState) = A.s+B.s +(-)(A::QuState,B::QuState) = A.s-B.s +(*)(A::QuState,B::QuState) = A.s*B.s +(/)(A::QuState,B::QuState) = A.s/B/s +(^)(A::QuState,b::Real) = A.s^b + + +# TODO +# show(io::IO,state::AbstractQuArray) +# show(io::IO,state::QuState) +# show(io::IO,state::QuIDDArray) + +type stlzState{N}<:AbstractQuArray{N} + X::AbstractMatrix # (2n+1)*n matrix for stabilizer/destabilizer x bits + Z::AbstractMatrix # (2n+1)*n matrix for z bits + R::AbstractVector # phase bits: 0 for +1, 1 for i, 2 for -1, 3 for -i. Normally either 0 or 2 +end + +stlzState(X::BitMatrix,Z::BitMatrix,R::BitVector) = stlzState{length(R)/2}(X,Z,R) + +function stlzState(n::Int) + X = spzeros(Bool,2n,n) + Z = spzeros(Bool,2n,n) + R = spzeros(Bool,2n) + + X[1:n,:] = speye(n) + Z[n+1:2n,:] = speye(n) + + return stlzState{n}(X,Z,R) +end + +function copy!(A::stlzState,B::stlzState) + copy!(A.X,B.X) + copy!(A.Z,B.Z) + copy!(A.R,B.R) + + return A +end diff --git a/src/const.jl b/src/const.jl new file mode 100644 index 0000000..49bbbbf --- /dev/null +++ b/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 diff --git a/src/utils/LogicExpr.jl b/src/utils/LogicExpr.jl new file mode 100644 index 0000000..30c4fce --- /dev/null +++ b/src/utils/LogicExpr.jl @@ -0,0 +1,208 @@ +################################################################### +# +# This file generates SAT's Instance in classical method for tests +# in quantum computing. +# Ref: arxiv:quant-ph/0007071v1 +# +################################################################### + + + +abstract AbstractBits +abstract AbstractClause{N} + +immutable Bits <: AbstractBits + value::UInt32 +end + +function call(b::Bits,n::Integer) + @assert n>0 + return (b.value>>(n-1))&1 +end + +""" +TruthTable Example: + +|bits value | truth value| +|-----------|------------| +| 000 | 0 | +| 001 | 1 | +| 010 | 0 | +| 011 | 1 | +| 100 | 0 | +| 101 | 1 | +| 110 | 0 | +| 111 | 0 | + +TruthTable(0b00101010) +""" +immutable TruthTable + value::UInt32 +end + +function call(truthtable::TruthTable,index::Integer) + return (truthtable.value>>index)&1 +end + +immutable Clause{N} <: AbstractClause{N} + table::TruthTable + ids::AbstractVector{Integer} + + function Clause(table::TruthTable,ids::AbstractVector{Integer}) + @assert length(ids)==N + sort!(ids) + new(table,ids) + end +end + +Clause(table::TruthTable,ids::Integer...) = Clause{length(ids)}(table,[ids...]) +Clause(num::Integer,table::TruthTable,ids::Integer...) = Clause{num}(table,[ids...]) +Clause(num::Integer,table::Integer,ids::Integer...) = Clause{num}(TruthTable(table),[ids...]) + +function call{N}(clause::Clause{N},assign::Integer) + return clause.table(assign) +end + +################################################################################ +# +# Clauses in exact cover problem +# +################################################################################ + +type ECClause{N} <: AbstractClause{N} + ids::Vector{Int} + + function ECClause(ids::Vector{Int}) + @assert length(ids) == N + sort!(ids) + new(ids) + end +end + +ECClause(ids::Vector{Int}) = ECClause{length(ids)}(ids) +ECClause(ids::Integer...) = ECClause{length(ids)}([ids...]) + +function call{N}(c::ECClause{N},assign::Integer) + res = 0 + for i = 1:N + res += assign&1 + assign = assign>>1 + end + return Int(res==1) +end + +function save{N}(io::IO,clause::ECClause{N}) + write(io,clause.ids[1]) + for id in clause.ids[2:end] + write(io,"\t") + write(io,id) + end + write(io,"\n") +end + +""" +Instance is a collection of M clauses: + +~~~TeX +C_1 Λ C_2 Λ... C_M +~~~ + +construct an instance: + +- `num` is the number of bits +- `clauses` is the collection of clause +""" +type Instance{M,N} + c::AbstractVector{AbstractClause{N}} +end + +Instance{N}(num::Integer,clauses::AbstractVector{ECClause{N}}) = Instance{num,N}(clauses) +Instance{N}(num::Integer,clause::ECClause{N},clauses::ECClause{N}...) = Instance(num,[clause,clauses...]) +Instance{N}(num::Integer,clauses::AbstractVector{Clause{N}}) = Instance{num,N}(clauses) +Instance{N}(num::Integer,clause::Clause{N},clauses::Clause{N}...) = Instance(num,[clause,clauses...]) + +function call{M,N}(clauses::Instance{M,N},assign::Bits) + res = 1 + for clause in clauses.c + assignment = 0 + digit = 0 + for id in clause.ids + assignment += assign(id)<shuffle)[1:N] + + clauses = [ECClause(ids)] + ins = Instance(n,clauses) + + for i = 1:maxtry + pre_ids = ids + ids = (list|>shuffle)[1:N] + if sort!(ids) != sort!(pre_ids) + push!(ins,ECClause(ids)) + end + + if AssignNum(ins)<2 + return ins + end + end + + return nothing +end + +function generate(n::Integer,maxiter=1000;maxbits=16,N=3) + for i=1:maxiter + t = engine(n,N) + t === nothing || return t + end + warn("may not have an assignment\n") +end + +export AbstractBits,AbstractClause,Bits,Clause,ECClause,TruthTable,save,show,generate,call,Instance diff --git a/src/utils/math.jl b/src/utils/math.jl new file mode 100644 index 0000000..66a78dc --- /dev/null +++ b/src/utils/math.jl @@ -0,0 +1,5 @@ +include("matrix.jl") + +function normalize!(vec::AbstractVector) + return vec[:]=vec/norm(vec) +end diff --git a/src/utils/matrix.jl b/src/utils/matrix.jl new file mode 100644 index 0000000..673b917 --- /dev/null +++ b/src/utils/matrix.jl @@ -0,0 +1,15 @@ +#trotter expansion +function trotter(A::AbstractMatrix,B::AbstractMatrix,P::Int64) + return (expm(full(A/(2*P)))*expm(full(B/P))*expm(full(A/(2*P))))^P +end + +function (⊗)(A::AbstractMatrix,B::AbstractMatrix) + @assert size(A)[1]==size(A)[2] + @assert size(B)[1]==size(B)[2] + + return kron(A,B) +end + +function (⊕)(A::AbstractMatrix,B::AbstractMatrix) + return full(blkdiag(sparse(A),sparse(B))) +end diff --git a/test/addgate.jl b/test/addgate.jl new file mode 100644 index 0000000..9bc9c74 --- /dev/null +++ b/test/addgate.jl @@ -0,0 +1,8 @@ +using QuCmp + +cc = Circuit(3) + +addgate!(cc,X,2,3) +addgate!(cc,X,1,2) + +@show cc diff --git a/test/chp.jl b/test/chp.jl new file mode 100644 index 0000000..c13fa95 --- /dev/null +++ b/test/chp.jl @@ -0,0 +1,28 @@ +using QuCmp +using Base.Test + +h = HadamardUnit(1) + +input = stlzState(2) + +process!(h,input) + + +@test input.X == [0 0;0 1;1 0;0 0]|>BitMatrix|>sparse +@test input.Z == [1 0;0 0;0 0;0 1]|>BitMatrix|>sparse +@test input.R == [0,0,0,0]|>BitVector|>sparse + +cnot = CNOTUnit(1,2) + +process!(cnot,input) + +@test input.X == [0 0;0 1;1 1;0 0]|>BitMatrix|>sparse +@test input.Z == [1 0;0 0;0 0;1 1]|>BitMatrix|>sparse +@test input.R == [0,0,0,0]|>BitVector|>sparse + +phase = PhaseUnit(2) +process!(phase,input) + +@test input.X == [0 0;0 1;1 1;0 0]|>BitMatrix|>sparse +@test input.Z == [1 0;0 1;0 1;1 1]|>BitMatrix|>sparse +@test input.R == [0,0,0,0]|>BitVector|>sparse diff --git a/test/process.jl b/test/process.jl new file mode 100644 index 0000000..468fe66 --- /dev/null +++ b/test/process.jl @@ -0,0 +1,9 @@ +using QuCmp + +cc = Circuit(2) + +addgate!(cc,X,1,1) + +t = process(cc,sparsevec([1,0,0,0])) +@show t +@show full(t)