diff --git a/LICENSE.md b/LICENSE.md index b35570c..b561db5 100644 --- a/LICENSE.md +++ b/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 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/AQCShrodingerEq.jl b/src/Adiabatic/AQCShrodingerEq.jl new file mode 100644 index 0000000..b6f3719 --- /dev/null +++ b/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 diff --git a/src/Adiabatic/Adiabatic.jl b/src/Adiabatic/Adiabatic.jl new file mode 100644 index 0000000..541dd0d --- /dev/null +++ b/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) diff --git a/src/Adiabatic/Cooling.jl b/src/Adiabatic/Cooling.jl new file mode 100644 index 0000000..9d1b9fb --- /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..11a86ca --- /dev/null +++ b/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 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/QuCmp.jl b/src/QuCmp.jl index 5afc76d..c56cb70 100644 --- a/src/QuCmp.jl +++ b/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 diff --git a/src/circuit/Circuit.jl b/src/circuit/Circuit.jl new file mode 100644 index 0000000..dcf79ff --- /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}<:AbstractQC{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..800171f --- /dev/null +++ b/src/utils/LogicExpr.jl @@ -0,0 +1,226 @@ +################################################################### +# +# This file generates SAT's Instance in classical method for tests +# in quantum computing. +# Ref: arxiv:quant-ph/0007071v1 +# +################################################################### + + +import Base: show,call +export AbstractBits, AbstractClause, Bits, + Clause, ECClause, TruthTable, save, + generate, Instance, AssignNum, readins + + +abstract AbstractBits +abstract AbstractClause{N} + +immutable Bits <: AbstractBits + value::UInt32 + + Bits(n::Integer) = new(n) +end + +function call(b::Bits,n::Integer) + @assert n>0 + return (b.value>>(n-1))&1 +end + +# function convert(Bits,Integer) +# return Bits() + +""" +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 + + asign = AssignNum(ins) + if asign[1]==1 + return ins,asign[2] + 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") + return nothing,nothing +end 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/adiabatic.jl b/test/adiabatic.jl new file mode 100644 index 0000000..34ffb87 --- /dev/null +++ b/test/adiabatic.jl @@ -0,0 +1,20 @@ +using QuCmp +using QuBase + +n = 5 + +ins,ans = generate(n) +pH = QuArray(pHamilton(ins,n),(ComputBasis(n),ComputBasis(n))) + +@show typeof(pH)<: QuBase.AbstractQuArray + +aqc = AQC(pH,n;maxtime=100) + +state = QuArray(convert(Array{Complex128,1},[1/sqrt(2^n) for i=1:2^n]),ComputBasis(n)) + +@time for (t,psi) in aqc + state = psi +end + + +print("p: $(norm(aqc.state[ans[1]+1])^2)\n") 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)