Skip to content

Commit

Permalink
Development for optimal macro
Browse files Browse the repository at this point in the history
  • Loading branch information
@theogf
theogf committed Aug 9, 2019
1 parent e1e1744 commit eab9839
Show file tree
Hide file tree
Showing 2 changed files with 196 additions and 15 deletions.
109 changes: 106 additions & 3 deletions dev/test_macro.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -9,7 +9,107 @@ g(y) = 0.0
γ(y) = 1.0
φ(r) = exp(-sqrt(r)/b)
∇φ(r) = -exp(-sqrt(r)/b)/(2*b*sqrt(r))
ll(y,x) = 0.5*exp(0.5*y*x)*sech(0.5*sqrt(x^2))

##
formula = :(p(y|x)=exp(0.5*y*x)*sech(0.5*sqrt(y^2 - 2*y*x + x^2)))
# formula = :(p(y,x)=exp(0.5*y*x)*sech(0.5*sqrt(0.0 - 0.0*x + x^2)))
formula.args[2].args[2].args

topargs = formula.args[2].args[2].args
if topargs[1] == :*
@show topargs[1]
global CC = copy(topargs[2])
popfirst!(topargs)
popfirst!(topargs)
else
global CC = :0
end
args2 = topargs[1]
if args2.args[1] == :exp
gargs = args2.args[2]
if gargs.args[1] == :*
deleteat!(gargs.args,findfirst(isequal(:x),gargs.args))
else
@error "BAD BAD BAD"
end
global GG = copy(gargs)
popfirst!(topargs)
else
global GG = :0
end
args3 = topargs[1]
seq = string(args3)
findh= r"\([^(]*\-.*x.*\+.*x \^ 2[^)]*"
b = occursin(findh,seq)
m = match(findh,seq).match
alphar = r"[^(][^-]*"
malpha = match(alphar,m).match[1:end-1]
betar = r"- [^x]*x"
mbeta = match(betar,m).match[3:end]
gammar = r"\+ [^x]*x \^ 2"
mgamma = match(gammar,m).match[3:end]

AA = :($malpha)
BB = :($(mbeta[1:end-1]))
GG = :($(mgamma == "x ^ 2" ? "1.0" : mgamma[1:end-5]))

loc = findfirst(findh,seq)
newseq = seq[1:loc[1]-1]*"r"*seq[(loc[end]+1):end]
PHI = :($newseq)
##
f_lap = :(p(y|x)=0.5/β * exp(- sqrt((y^2 - 2*y*f + f^2))/β))
display.(AGP.@augmodel NewSuperLaplace Regression (p(y|x)=0.5/β * exp( y * x) * exp(- sqrt((sqrt(y^2) - exp(4.0*y)*x + 2.0*x^2))/β)) β)
pdfstring = "(0.5 / β) * exp(2*y*x) * exp(-(sqrt((y ^ 2 - 2.0 * y * x) + 1.0*x ^ 2)) / β)"

Gstring
##
PhiHstring = match(Regex("(?<=$(AGP.correct_parenthesis(Gstringfull))x\\) \\* ).*"),pdfstring).match
Hstring = match(r"(?<=\().+\-.*x.*\+.+x \^ 2(?=\))",PhiHstring).match
locx = findfirst("x ^ 2",PhiHstring)
count_parenthesis = 1
locf = locx[1]
while count_parenthesis != 0
global locf = locf - 1
println(PhiHstring[locf])
if PhiHstring[locf] == ')'
global count_parenthesis += 1
elseif PhiHstring[locf] == '('
global count_parenthesis -= 1
end
end
Hstring = PhiHstring[(locf+1):locx[end]]

alphar = r"[^(][^-]*"
alpha_string = match(alphar,Hstring).match[1:end-1]
# betar = r"(?>=- )[^x]+(?= * x)"
betar = r"(?<=- )[^x]+(?= * x)"
mbeta = match(betar,Hstring).match
while last(mbeta) == ' ' || last(mbeta) == '*'
global mbeta = mbeta[1:end-1]
end
mbeta
gammar = r"(?<=\+ )[^x]*(?=x \^ 2)"
mgamma = match(gammar,m).match == "" ? "1.0" : match(gammar,m).match
##
findnext(isequal(')'),PhiHstring,locx[end])
code = Meta.parse(PhiHstring)
code.args

S = code.args[2].args[2].args[2].args[2].args[2].args
S = code.args[2].args[2].args[2].args[2].args
for args in S.args
if args == :(x ^ 2)
@show "BLAH"
end
end
S = string(code.args[2].args[2])
Hstring = match(r"(?<=\().*x \^ 2.*\-.*x.*\+.*(?=\))",S)

##

txt = AGP.@augmodel("NewLaplace","Regression",C,g,α,β,γ,φ,∇φ)

# NewLaplaceLikelihood() |> display
N = 500
σ = 1.0
Expand All @@ -19,6 +119,7 @@ L = Matrix(cholesky(K).L)
y_true = rand(MvNormal(K))
y = y_true+randn(length(y_true))*2
p = scatter(X[:],y,lab="data")
NewLaplaceLikelihood() |> display
m = VGP(X,y,RBFKernel(0.5),NewLaplaceLikelihood(),AnalyticVI(),optimizer=false)
train!(m,iterations=100)
y_p, sig_p = proba_y(m,collect(0:0.01:1))
Expand All @@ -30,8 +131,10 @@ y_p2, sig_p2 = proba_y(m2,collect(0:0.01:1))
plot!(X,y_true,lab="truth")

plot!(collect(0:0.01:1),y_p,lab="Auto Laplace")
plot!(collect(0:0.01:1),y_p2,lab="Classic Laplace")
plot!(collect(0:0.01:1),y_p2,lab="Classic Laplace") |> display


# @btime train!($m,iterations=1)
# @btime train!($m2,iterations=1)

@btime train!($m,iterations=1)
@btime train!($m2,iterations=1)
###
102 changes: 90 additions & 12 deletions src/likelihood/generic_likelihood.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3,26 +3,95 @@ macro augmodel()
end

"""
**$(local sname) Likelihood**
** Likelihood**
Template file for likelihood creation
```julia
$(local lname)()
()
```
See all functions you need to implement
---
"""

check_model_name(name::Symbol) = !isnothing(match(r"[a-zA-Z]*",string(name)))
function check_likelihoodtype(ltype::Symbol)
if ltype == :Regression || ltype == :Classification || ltype == :Event
return true
else
return false
end
end

function treat_likelihood(likelihood::Expr)

end

function treat_params(params)

end

function correct_parenthesis(text::AbstractString)
replace(text,r"(?=[()*+-.])"=>"\\")
end

const AGP=AugmentedGaussianProcesses
macro augmodel(name::String,LikelihoodType::String,C::Symbol,g::Symbol,α::Symbol,β::Symbol,γ::Symbol,φ::Symbol,∇φ::Symbol)

macro augmodel(name::Symbol,likelihoodtype::Symbol,likelihood::Expr,params)
latent = :x
@assert occursin(r"p\(y\s*|\s*\w\)",string(likelihood)) "Likelihood should be of the form p(y|x) = C*exp(g(y)*x)*φ(α(y)-β(y)*x+γ(y)*x^2), replacing all functions by 0 if necessary"
pdf_string = string(likelihood.args[2].args[2])
# @show pdfexpr = Meta.parse(pdfstring)
C_string = match(r".*?(?= (\* exp\(.*x\)))",pdf_string).match
G_string_f = match(Regex("(?<=$(AGP.correct_parenthesis(C_string)) \\* exp\\().*(?=x\\) \\*)"),pdf_string).match
G_string = deepcopy(G_string_f)
while last(G_string) == ' ' || last(G_string) == '*'
G_string = G_string[1:end-1]
end
phi_h_string = match(Regex("(?<=$(AGP.correct_parenthesis(G_string_f))x\\) \\* ).*"),pdf_string).match
loc_x² = findfirst("x ^ 2",phi_h_string)
count_parenthesis = 1
loc_start = loc_x²[1]
while count_parenthesis != 0
loc_start = loc_start - 1
if phi_h_string[loc_start] == ')'
count_parenthesis += 1
elseif phi_h_string[loc_start] == '('
count_parenthesis -= 1
end
end
h_string = phi_h_string[(loc_start+1):loc_x²[end]]
phi_string = phi_h_string[1:loc_start]*"r"*phi_h_string[(loc_x²[end]+1):end]
@show alpha_string = match(r"[^(][^-]*",h_string).match[1:end-1]
gamma_string = match(r"(?<=\+ )[^x]*(?=x \^ 2)",h_string).match
gamma_string = gamma_string == "" ? "1.0" : gamma_string
while last(gamma_string) == ' ' || last(gamma_string) == '*'
gamma_string = gamma_string[1:end-1]
end
beta_string = match(Regex("(?<=$(AGP.correct_parenthesis(alpha_string)) - )[^( x )]*(?= x)"),h_string).match
while last(beta_string) == ' ' || last(beta_string) == '*'
beta_string = beta_string[1:end-1]
end
return (C_string,G_string,phi_string,alpha_string,beta_string,gamma_string)
# treat_params(params)
# C,g,α,β,γ,φ = treat_likelihood(likelihood)
end

macro augmodel(name::Symbol,likelihoodtype::Symbol,likelihood::Expr)

end


macro augmodel(name::Symbol,likelihoodtype::Symbol,C::Symbol,g::Symbol,α::Symbol,β::Symbol,γ::Symbol,φ::Symbol,∇φ::Symbol)
#### Check args here
#Check name has no space
#Check LikelihoodType exists
@assert check_model_name(name) "Please only use alphabetic characters for the name of the likelihood"
@assert checkl_likelihoodtype(likelihoodtype) "Please use a correct likelihood type : Regression, Classification or Event"
#Find gradient with AD if needed
#In a later stage try to find structure automatically
esc(_augmodel(Symbol(name),Symbol(name,"Likelihood"),Symbol(name,"Likelihood{T}"),Symbol(LikelihoodType*"Likelihood"),C,g,α,β,γ,φ,∇φ))
esc(_augmodel(name,Symbol(name,"Likelihood"),Symbol(name,"Likelihood{T}"),Symbol(LikelihoodType*"Likelihood"),C,g,α,β,γ,φ,∇φ))
end
function _augmodel(name,lname,lnameT,ltype,C,g,α,β,γ,φ,∇φ)
quote begin
Expand Down Expand Up @@ -74,19 +143,23 @@ function _augmodel(name,lname,lnameT,ltype,C,g,α,β,γ,φ,∇φ)
end

function φ(l::$(lname),r::T) where {T}
φ.(r)
φ(r)
end

function ∇φ(l::$(lname),r::T) where {T}
∇φ.(r)
∇φ(r)
end

function ∇²φ(l::$(lname),r::T) where {T}
ForwardDiff.gradient(x->∇φ(l,x[1]),r)[1]
end

function pdf(l::$(lname),y::Real,f::Real)
C()*exp(g(y)*f)*φ(α(y)-β(y)*f-α(y)*f)
C()*exp(g(y)*f)*φ(α(y)-β(y)*f+α(y)*f)
end

function Base.show(io::IO,model::$(lname){T}) where {T}
print(io,"$(name) Likelihood")
print(io,"Generic Likelihood")#WARNING TODO, to be fixed!
end

function AGP.compute_proba(l::$(lname){T},μ::AbstractVector{T},σ²::AbstractVector{T}) where {T<:Real}
Expand Down Expand Up @@ -161,10 +234,15 @@ function _augmodel(name,lname,lnameT,ltype,C,g,α,β,γ,φ,∇φ)

### Gradient Section ###

function gradpdf(::$(lname),y::Int,f::T) where {T<:Real}
@inline function grad_log_pdf(l::$(lname){T},y::Real,f::Real) where {T<:Real}
h² = α(l,y) - β(l,y)*f + γ(l,y)*f^2
g(l,y)+(-β(l,y)+2*γ(l,y))*∇φ(l,h²)/φ(l,h²)
end

function hessiandiagpdf(::$(lname),y::Int,f::T) where {T<:Real}
@inline function hessian_log_pdf(l::$(lname){T},y::Real,f::Real) where {T<:Real}
h² = α(l,y) - β(l,y)*f + γ(l,y)*f^2
ϕ = φ(l,h²); ∇ϕ = ∇φ(l,h²); ∇²ϕ = ∇²φ(l,h²)
2*γ(l,y)*∇ϕ/ϕ-(-β(l,y)+2*γ(l,y)*f^2)*∇ϕ^2/ϕ^2+(-β(l,y)+2*γ(l,y)*f^2)^2*∇²ϕ/ϕ
end

end
Expand Down

0 comments on commit eab9839

Please sign in to comment.