Example of using the BEAGLE library in Julia

hellobeagle.jl

@printf "BEAGLE version %s\n" bytestring(ccall((:beagleGetVersion,"libhmsbeagle"), Ptr{Cchar}, (),))

@printf "%s\n" bytestring(ccall((:beagleGetCitation,"libhmsbeagle"), Ptr{Cchar}, (),))

#typealias Cint Int64

const BEAGLE_OP_COUNT = 7
const BEAGLE_OP_NONE = -1

immutable BeagleResource
    name::Ptr{UInt8}
    description::Ptr{UInt8}
    supportFlags::Clong
    requiredFlags::Clong
end

immutable BeagleResourceList
    list::Ptr{BeagleResource}
    length::Cint
end

immutable BeagleInstanceDetails
    resourceNumber::Cint
    resourceName::Ptr{UInt8}
    implName::Ptr{UInt8}
    implDescription::Ptr{UInt8}
    flags::Clong
end

immutable BeagleOperation
    destinationPartials::Cint
    destinationScaleWrite::Cint
    destinationScaleRead::Cint
    child1Partials::Cint
    child1TransitionMatrix::Cint
    child2Partials::Cint
    child2TransitionMatrix::Cint
end

resourceList = ccall( (:beagleGetResourceList, "libhmsbeagle"),
    Ptr{BeagleResourceList},
    (),
    )
rl = unsafe_load(resourceList)
rl.length
resources = pointer_to_array(rl.list,rl.length)
for i in 1:rl.length
    @printf "BEAGLE Resource %i\n" i-1
    @printf "Name: %s\n" bytestring(resources[i].name)
    @printf "Descriptions: %s\n" bytestring(resources[i].description)
end

mars = "CCGAG-AGCAGCAATGGAT-GAGGCATGGCG"
saturn  = "GCGCGCAGCTGCTGTAGATGGAGGCATGACG"
jupiter = "GCGCGCAGCAGCTGTGGATGGAAGGATGACG"

nPatterns = length(mars)

retInfo = [BeagleInstanceDetails(0,0,0,0,0)]
instance = ccall((:beagleCreateInstance, "libhmsbeagle"),
    Int,
    (Cint,Cint,Cint,Cint,Cint,Cint,Cint,Cint,Cint,Ptr{Int},Cint,Clong,Clong,Ptr{BeagleInstanceDetails}),
    3,2,3,4,nPatterns,1,4,1,0,C_NULL,0,0,0,retInfo)

table = Dict{Char,Int64}(
    'A' => 0,
    'C' => 1,
    'G' => 2,
    'T' => 3,
    'a' => 0,
    'c' => 1,
    'g' => 2,
    't' => 3,
    '-' => 4,
    )

marsStates = Int64[table[b] for b in mars]
saturnStates = Int64[table[b] for b in saturn]
jupiterStates = Int64[table[b] for b in jupiter]

ccall((:beagleSetTipStates, "libhmsbeagle"),
    Int,
(Cint,Cint,Ptr{Int64}),
    instance,0,marsStates)
ccall((:beagleSetTipStates, "libhmsbeagle"),
    Int,
(Cint,Cint,Ptr{Int64}),
    instance,1,saturnStates)
ccall((:beagleSetTipStates, "libhmsbeagle"),
    Int,
(Cint,Cint,Ptr{Int64}),
    instance,2,jupiterStates)

patternWeights = Array(Float64,nPatterns)
fill!(patternWeights,1.0)
ccall((:beagleSetPatternWeights, "libhmsbeagle"),
    Int,
(Cint,Ptr{Cdouble}),
    instance,patternWeights)

freqs=Float64[0.25,0.25,0.25,0.25]
ccall((:beagleSetStateFrequencies, "libhmsbeagle"),
    Int,
(Cint,Cint,Ptr{Cdouble}),
    instance,0,freqs)

weights = Float64[1.0]
rates = Float64[1.0]
ccall((:beagleSetCategoryWeights, "libhmsbeagle"),
    Int,
(Cint,Cint,Ptr{Cdouble}),
instance,0,weights)
ccall((:beagleSetCategoryRates, "libhmsbeagle"),
    Int,
(Cint,Ptr{Cdouble}),
    instance,rates)

evec=Float64[1.0,2.0,0.0,0.5,
    1.0,-2.0,0.5,0.0,
    1.0,2.0,0.0,-0.5,
    1.0,-2.0,-0.5,0.0
]
ivec=Float64[
    0.25,0.25,0.25,0.25,
    0.125,-0.125,0.125,-0.125,
    0.0,1.0,0.0,-1.0,
    1.0 ,0.0,-1.0,0.0
]
ev = Float64[0.0, -1.3333333333333333, -1.3333333333333333, -1.3333333333333333]
ccall((:beagleSetEigenDecomposition, "libhmsbeagle"),
    Int,
(Cint,Cint,Ptr{Cdouble},Ptr{Cdouble},Ptr{Cdouble}),
    instance,0,evec,ivec,ev)

nodeIndices = Int64[0,1,2,3]
edgeLengths = Float64[0.1,0.1,0.2,0.1]

ccall((:beagleUpdateTransitionMatrices, "libhmsbeagle"),
    Cint,
(Cint,Cint,Ptr{Int64},Ptr{Int64},Ptr{Int64},Ptr{Cdouble},Cint),
    instance,0,nodeIndices,C_NULL,C_NULL,edgeLengths,4)

qmat = [Array(Float64,4,4)]
ccall((:beagleGetTransitionMatrix,"libhmsbeagle"),Cint,(Cint,Cint,Ptr{Cdouble}),instance,0,qmat)
qmat

operations = BeagleOperation[
    BeagleOperation(3, BEAGLE_OP_NONE, BEAGLE_OP_NONE, 0, 0, 1, 1),
    BeagleOperation(4, BEAGLE_OP_NONE, BEAGLE_OP_NONE, 2, 2, 3, 3)
        ]
ccall((:beagleUpdatePartials, "libhmsbeagle"),
    Cint,
    (Cint,Ptr{BeagleOperation},Cint,Cint),
instance,operations,2,BEAGLE_OP_NONE)

lenpartial = 4*31*1 # stateCount * patternCount * categoryCount
partials3 = [Array(Float64,lenpartial)]
ccall((:beagleGetPartials, "libhmsbeagle"),
    Cint,
    (Cint,Cint,Cint,Ptr{Cdouble}),
    instance,3,BEAGLE_OP_NONE,partials3
)
partials3

partials4 = [Array(Float64,lenpartial)]
ccall((:beagleGetPartials, "libhmsbeagle"),
    Cint,
    (Cint,Cint,Cint,Ptr{Cdouble}),
    instance,4,BEAGLE_OP_NONE,partials4
)
partials4

logL = Ref{Cdouble}(0.0)

rootIndex = [4]
categoryWeightIndex = [0]
stateFrequencyIndex = [0]
cumulativeScaleIndex = [BEAGLE_OP_NONE]

ccall((:beagleCalculateRootLogLikelihoods, "libhmsbeagle"), Cint, (Cint,Ptr{Int64},Ptr{Int64},Ptr{Int64},Ptr{Int64},Cint,Ref{Cdouble}), instance,rootIndex,categoryWeightIndex,stateFrequencyIndex,cumulativeScaleIndex,1,logL)

logL[]

sitelike=Array(Float64,nPatterns)
ccall((:beagleGetSiteLogLikelihoods,"libhmsbeagle"),Cint,(Cint,Ptr{Cdouble}),instance,sitelike)
sitelike

ccall((:beagleFinalizeInstance,"libhmsbeagle"), Cint, (Cint,),instance)
ccall((:beagleFinalize,"libhmsbeagle"),Cint,(),)