Propagators before doc update

propagators.jl

using QuBase
using Docile

import Base: start, next, done

abstract QuPropagatorMethod

immutable QuPropagator{QPM<:QuPropagatorMethod}
    hamiltonian
    init_state::QuVector
    tlist
    method::QPM
    #options::Dict
    #QuPropagator(hamiltonian, init_state, tlist, method, options) = new(hamiltonian, init_state, tlist, method, options)
    QuPropagator(hamiltonian, init_state, tlist, method) = new(hamiltonian, init_state, tlist, method)
end

# QuPropagator{QPM<:QuPropagatorMethod}(hamiltonian, init_state, tlist, method::QPM, options::Dict) = QuPropagator{QPM}(hamiltonian, init_state, tlist, method, options)
QuPropagator{QPM<:QuPropagatorMethod}(hamiltonian, init_state, tlist, method::QPM) = QuPropagator{QPM}(hamiltonian, init_state, tlist, method)

immutable QuPropagatorState
    psi::QuVector
    t
    t_state
end

function Base.start(prob::QuPropagator)
    t_state = start(prob.tlist)
    t,t_state = next(prob.tlist,t_state)
    return QuPropagatorState(prob.init_state,t,t_state)
end

Base.done(prob::QuPropagator, qustate::QuPropagatorState) = done(prob.tlist, qustate.t_state)

function Base.next{QPM<:QuPropagatorMethod}(prob::QuPropagator{QPM}, qustate::QuPropagatorState)
    op = prob.hamiltonian
    current_qustate = qustate.psi
    current_t = qustate.t
    if !done(prob, qustate)
        t,t_state = next(prob.tlist, qustate.t_state)
        next_qustate = propagate(prob.method, op, t, current_t, current_qustate)
        return (next_qustate, t), QuPropagatorState(next_qustate, t, t_state)
    else
        return QuPropagatorState(current_qustate, current_t, qustate.t_state)
    end
end

immutable QuEuler <: QuPropagatorMethod
end

function propagate(prob::QuEuler, op, t, current_t, current_qustate)
    dt = t - current_t
    return (eye(op)-im*op*dt)*current_qustate
end

immutable QuCN <: QuPropagatorMethod
end

function propagate(prob::QuCN, op, t, current_t, current_qustate)
    dt = t - current_t
    uni = eye(op)-im*op*dt/2
    return \(uni', uni*current_qustate)
end

immutable QuKrylov <: QuPropagatorMethod
    options::Dict
end

function propagate(prob::QuKrylov, op, t, current_t, current_qustate)
    basis_size = prob.options[:NC]
    dt = t-current_t
    # N = min(basis_size, length(qustate.psi))
    N = min(basis_size, length(current_qustate))
    # v = Array(typeof(qustate.psi),0)
    v = Array(typeof(current_qustate),0)
    # push!(v,zeros(qustate.psi))
    push!(v,zeros(current_qustate))
    push!(v,current_qustate)
    # w = Array(typeof(qustate.psi),0)
    w = Array(typeof(current_qustate),0)
    alpha = Array(Complex{Float64},0)
    beta = Array(Complex{Float64},0)
    push!(beta,0.)
    for i=2:N
        # push!(w, prob.hamiltonian*v[i])
        push!(w, op*v[i])
        push!(alpha, w[i-1]'*v[i])
        w[i-1] = w[i-1]-alpha[i-1]*v[i]-beta[i-1]*v[i-1]
        push!(beta, norm(w[i-1]))
        push!(v, w[i-1]/beta[i])
    end
    # push!(w, prob.hamiltonian*v[end])
    push!(w, op*v[end])
    push!(alpha, w[end]'*v[end])
    deleteat!(v,1)
    H_k = full(Tridiagonal(beta[2:end], alpha, beta[2:end]))
    U_k = expm(-im*dt*H_k)
    next_state = zeros(current_qustate)
    for j=1:N
        next_state = next_state + v[j]*U_k[j,1]
    end
    return next_state
end

#@doc """
#Euler Problem
#`hamiltonian` : Hamiltonian of the system
#`init_state`  : initial state of the syste
#`tlist`       : Time step array
#""" ->
#immutable QuEuler{B<:AbstractBasis,T,A} <: QuantumPropagator{B,T,A}
"""immutable QuEuler{B<:AbstractBasis,T,A}
    hamiltonian::QuMatrix{B,T,A}
    init_state::QuVector
    tlist::Array
    QuEuler(hamiltonian, init_state, tlist) = new(hamiltonian, init_state, tlist)
end


QuEuler{B<:AbstractBasis,T,A}(hamiltonian::QuMatrix{B,T,A}, init_state::QuVector, tlist::Array)=QuEuler{B,T,A}(hamiltonian, init_state, tlist)
"""

#@doc """
#Euler State
#`psi`     : State of the system
#`t`       : Current time in time step `tlist` array
#`t_state` : Index of the next time state in the `tlist` array
#""" ->
"""immutable QuEulerState
    psi::QuVector
    t
    t_state
end
"""

#@doc """
#Input Parameters : Euler Problem

#Returns the starting, i.e., the initial state of the system
#""" ->
"""function Base.start(prob::QuEuler)
    t_state = start(prob.tlist)
    t,t_state = next(prob.tlist,t_state)
    return QuEulerState(prob.init_state,t,t_state)
end
"""

#@doc """
#Input Parameters : Euler Problem and Euler State

#Returns the next state using Euler's method.
#""" ->
"""function Base.next(prob::QuEuler, qustate::QuEulerState)
    op = prob.hamiltonian
    t_list = prob.tlist
    current_qustate = qustate.psi
    current_t = qustate.t
    if current_t != t_list[end]
      t,t_state = next(t_list, qustate.t_state)
      dt = t - current_t
      next_state = (eye(op)-im*op*dt)*current_qustate
      return (next_state, t), QuEulerState(next_state,t,t_state)
    else
      return (current_qustate,current_t), QuEuler(current_qustate,current_t,qustate.t_state)
    end
end
"""

#@doc """
#Input Parameters : Euler Problem and Euler State

#Returns true if the current state is final state, else false
#""" ->
# Base.done(prob::QuEuler, qustate::QuEulerState) = done(prob.tlist, qustate.t_state)

@doc """
Internal function to initialize the system
""" ->
function initialize(hamiltonian::QuMatrix, init_state::QuVector, tlist::Array)
    final_state = Array(typeof(init_state),0)
    sizehint(final_state, length(tlist))
    push!(final_state, init_state)
    return final_state
end


@doc """
Input Parameters : Hamiltonian of the system, Initial State of the system, Time steps

Returns the collection of output states at every time step using the Forward Euler Method.
""" ->
function solver_euler(hamiltonian::QuMatrix, init_state::QuVector, tlist::Array)
    final_state = initialize(hamiltonian, init_state, tlist)
    #for i=2:length(tlist)
    #    uni = QuArray(eye(size(hamiltonian,1)))-im*hamiltonian*(tlist[i]-tlist[i-1])
    #    push!(final_state, uni*final_state[i-1])
    #end
    id = eye(hamiltonian)
    next_step = start(tlist)
    while !done(tlist, next_step)
        step, next_step = next(tlist, next_step)
        if next_step == length(tlist)+1
          break
        else
          uni = id-im*hamiltonian*(tlist[next_step]-step)
          push!(final_state, uni*final_state[next_step-1])
        end
    end
    return final_state
end


#@doc """
#Crank Nicolson Problem
#`hamiltonian` : Hamiltonian of the system
#`init_state`  : initial state of the syste
#`tlist`       : Time step array
#""" ->
#immutable QuCN{B<:AbstractBasis,T,A} <: QuantumPropagator{B,T,A}
"""immutable QuCN{B<:AbstractBasis,T,A}
    hamiltonian::QuMatrix{B,T,A}
    init_state::QuVector
    tlist::Array
    QuCN(hamiltonian, init_state, tlist) = new(hamiltonian, init_state, tlist)
end
"""

#@doc """
#Crank Nicolson State
#`psi`     : State of the system
#`t`       : Current time in time step `tlist` array
#`t_state` : Index of the next time state in the `tlist` array
#""" ->
"""
immutable QuCNState
    psi::QuVector
    t
    t_state
end
"""
# QuCN{B<:AbstractBasis,T,A}(hamiltonian::QuMatrix{B,T,A}, init_state::QuVector, tlist::Array)=QuCN{B,T,A}(hamiltonian, init_state, tlist)


#@doc """
#Input Parameters : Crank Nicolson Problem

#Returns the starting, i.e., the initial state of the system
#""" ->
"""function Base.start(prob::QuCN)
    t_state = start(prob.tlist)
    t,t_state = next(prob.tlist,t_state)
    return QuCNState(prob.init_state,t,t_state)
end
"""

#@doc """
#Input Parameters : Crank Nicolson Problem and Crank Nicolson State

#Returns the next state using Euler's method.
#""" ->
"""function Base.next(prob::QuCN, qustate::QuCNState)
    op = prob.hamiltonian
    t_list = prob.tlist
    current_qustate = qustate.psi
    current_t = qustate.t
    if current_t != t_list[end]
        t,t_state = next(t_list,qustate.t_state)
        dt = t - current_t
        uni = eye(op)-im*op*dt/2
        next_state = \(uni', uni*current_qustate)
        return (next_state,t), QuCNState(next_state,t,t_state)
    else
        return (current_qustate,current_t),QuCNState(current_qustate,current_t,qustate.t_state)
    end
end
"""

#@doc """
#Input Parameters : Crank Nicolson Problem and Crank Nicolson State

#Returns true if the current state is final state, else false
#""" ->
# Base.done(prob::QuCN, qustate::QuCNState) = done(prob.tlist,qustate.t_state)


@doc """
Input Parameters : Hamiltonian of the system, Initial State of the system, Time steps

Returns the collection of output states at every time step using Crank Nicolson Method.
""" ->
function solver_CN(hamiltonian::QuMatrix, init_state::QuVector, tlist::Array)
    final_state = initialize(hamiltonian, init_state, tlist)
    #for i=2:length(tlist)
    #    uni = QuArray(eye(size(hamiltonian,1)))-(im*hamiltonian*(tlist[i]-tlist[i-1]))/2
    #    push!(final_state, \(uni',uni*final_state[i-1]))
    #end
    id = eye(hamiltonian)
    next_step = start(tlist)
    while !done(tlist, next_step)
        step, next_step = next(tlist, next_step)
        if next_step == length(tlist)+1
          break
        else
          uni = id-(im*hamiltonian*(tlist[next_step]-step))/2
          push!(final_state, \(uni',uni*final_state[next_step-1]))
        end
    end
    return final_state
end


#@doc """
#Krylov Problem
#`hamiltonian` : Hamiltonian of the system
#`init_state`  : initial state of the syste
#`tlist`       : Time step array
#""" ->
#immutable QuKrylov{B<:AbstractBasis,T,A} <: QuantumPropagator{B,T,A}
"""immutable QuKrylov{B<:AbstractBasis,T,A}
    hamiltonian::QuMatrix{B,T,A}
    init_state::QuVector
    tlist::Array
    options::Dict
    QuKrylov(hamiltonian, init_state, tlist, options) = new(hamiltonian, init_state, tlist, options)
end
"""

# QuKrylov{B<:AbstractBasis,T,A}(hamiltonian::QuMatrix{B,T,A}, init_state::QuVector, tlist::Array, options::Dict)=QuKrylov{B,T,A}(hamiltonian, init_state, tlist,options)


#@doc """
#Krylov State
#`psi`     : State of the system
#`t`       : Current time in time step `tlist` array
#`t_state` : Index of the next time state in the `tlist` array
#""" ->
"""immutable QuKrylovState
    psi::QuVector
    t
    t_state
end
"""

#@doc """
#Input Parameters : QuKrylov Problem

#Returns the starting, i.e., the initial state of the system
#""" ->
"""function Base.start(prob::QuKrylov)
    t_state = start(prob.tlist)
    t,t_state = next(prob.tlist,t_state)
    return QuKrylovState(prob.init_state,t,t_state)
end
"""

#@doc """
#Input Parameters : Krylov Problem and Krylov State

#Returns the next state using Krylov method.
#""" ->
"""function Base.next(prob::QuKrylov, qustate::QuKrylovState)
    op = prob.hamiltonian
    t_list = prob.tlist
    basis_size = prob.options[:NC]
    current_t = qustate.t
    if current_t != t_list[end]
        t,t_state = next(t_list,qustate.t_state)
        dt = t-current_t
        N = min(basis_size, length(qustate.psi))
        v = Array(typeof(qustate.psi),0)
        push!(v,zeros(qustate.psi))
        push!(v,qustate.psi)
        w = Array(typeof(qustate.psi),0)
        alpha = Array(Complex{Float64},0)
        beta = Array(Complex{Float64},0)
        push!(beta,0.)
        for i=2:N
            push!(w, prob.hamiltonian*v[i])
            push!(alpha, w[i-1]'*v[i])
            w[i-1] = w[i-1]-alpha[i-1]*v[i]-beta[i-1]*v[i-1]
            push!(beta, norm(w[i-1]))
            push!(v, w[i-1]/beta[i])
        end
        push!(w, prob.hamiltonian*v[end])
        push!(alpha, w[end]'*v[end])
        deleteat!(v,1)
        H_k = full(Tridiagonal(beta[2:end], alpha, beta[2:end]))
        U_k = expm(-im*dt*H_k)
        next_state = zeros(qustate.psi)
        for j=1:N
            next_state = next_state + v[j]*U_k[j,1]
        end
        return (next_state,t),QuKrylovState(next_state,t,t_state)
    else
        return (qustate.psi,current_t),QuKrylovState(qustate.psi,current_t,qustate.t_state)
    end
end
"""

#@doc """
#Input Parameters : Krylov Problem and Krylov State

#Returns true if the current state is final state, else false
#""" ->
# Base.done(prob::QuKrylov, qustate::QuKrylovState) = done(prob.tlist,qustate.t_state)


export QuEuler,
      QuEulerState,
      solver_euler,
      QuCN,
      QuCNState,
      solver_CN,
      QuKrylov,
      QuKrylovState,
      start,
      next,
      done