MCMC in GaussianProcesses.jl

[chris-nemeth]

Hi Theo,

Thanks for fixing the deprecation errors.

We have a world age issue with our default Klara implementation within GaussianProcesses.jl. Our mcmc function can be found here mcmc. As you can see it's a pretty straightforward MALA sampler.

If you run the following as a test you'll see the error

using GaussianProcesses
import Distributions:Normal

srand(112233)
X = rand(20);
X = sort(X);
y = sin.(10*X);
y=convert(Vector{Bool}, y.>0);

#Select mean, kernel and likelihood function
mZero = MeanZero()   #Zero mean function
kern = SE(0.0,0.0)   #Sqaured exponential kernel (note that hyperparameters are on the log scale)
lik = BernLik()      #Bernoulli likelihood

gp = GP(X,vec(y),mZero,kern,lik)     

set_priors!(gp.k,[Normal(0.0,2.0),Normal(0.0,2.0)])

samples = mcmc(gp)

Thanks,

Chris

[papamarkou]

Thanks Chris, leave it with me and I will check it. These world age issues in Julia v0.6 have become my nightmare... if anything else fails, I might have to remove code generation from the package...

[chris-nemeth]

You may have already seen this and not sure if their solution helps world age?

[papamarkou]

I am trying to reproduce your error Chris, to start with, but I get the error

julia> lik = BernLik()      #Bernoulli likelihood
ERROR: UndefVarError: BernLik not defined

Where is BernLik() defined, do I need to import it from somewhere? You may pass on a fully operational script-example (or fully breaking considering the world age issue), and I will work with it.

[chris-nemeth]

Sorry. I should have mentioned that the non-Gaussian likelihoods aren't on the Julia package repository. You'll need to clone the package from Github.

[papamarkou]

I can see the error now Chris. I will ask the community if there is a way to solve the world age issue, or else the package needs re-development.

[papamarkou]

Here is the discussion at discourse:

https://discourse.julialang.org/t/world-age-issue-due-to-code-generation-inside-an-inner-constructor-in-julia-v0-6/5124

[chris-nemeth]

Hopefully there's a simple solution to this and it doesn't create too much extra work for you.

[papamarkou]

Indeed, let's hope Chris, as now nearly everything is code-generated and (world-) aged...

[papamarkou]

@ararslan, @JeffBezanson, @Keno, I know I have been asked before to post on discourse rather than tagging people. Unfortunately discourse hasn't been so helpful, and as you can see here my issue remains unresolved. Without of course wanting you to fix the current Julia issue opened by @chris-nemeth, can you possibly provide a working solution for the issue mentioned at discourse (although I am afraid there isn't one, because some loss of functionality was incurred from Julia v0.5 to v0.6 in relation to this matter)? Your help would be greatly appreciated, as my package is currently broken due to this problem.

[papamarkou]

@chris-nemeth, as you will see below, the following works. I extracted the code from your mcmc() method and it now executes. So, weirdly enough, the problem appears only when the same code is enclosed into a function (your mcmc()):

using Klara
using GaussianProcesses
import Distributions:Normal

srand(112233)
X = rand(20);
X = sort(X);
y = sin.(10*X);
y=convert(Vector{Bool}, y.>0);

#Select mean, kernel and likelihood function
mZero = MeanZero()   #Zero mean function
kern = SE(0.0,0.0)   #Sqaured exponential kernel (note that hyperparameters are on the log scale)
lik = BernLik()      #Bernoulli likelihood

gp = GP(X,vec(y),mZero,kern,lik)     

set_priors!(gp.k,[Normal(0.0,2.0),Normal(0.0,2.0)])

start = GaussianProcesses.get_params(gp)
sampler=Klara.MALA(0.1)
mcrange=Klara.BasicMCRange(nsteps=5000, burnin=1000)

function logtarget(hyp::Vector{Float64})  #log-target
  GaussianProcesses.set_params!(gp, hyp)
  return GaussianProcesses.update_target!(gp)   
end

function dlogtarget(hyp::Vector{Float64}) #gradient of the log-target
  GaussianProcesses.set_params!(gp, hyp)
  return GaussianProcesses.update_target_and_dtarget!(gp)
end

starting = Dict(:p=>start)
q = BasicContMuvParameter(:p, logtarget=logtarget, gradlogtarget=dlogtarget) 
model = likelihood_model(q, false)               #set-up the model
tune = AcceptanceRateMCTuner(0.6, verbose=true)  #set length of tuning (default to burnin length)
job = BasicMCJob(model, sampler, mcrange, starting)   #set-up MCMC job
run(job)
chain = Klara.output(job)

[papamarkou]

P.S. Of course this doesn't solve the problem, just letting you know what happens.

[chris-nemeth]

You are right that it works outside of the mcmc() function although I can't see why that is.

[papamarkou]

I don't know either Chris, it has to do with the new world age counter introduced in Julia v0.6 that messes up Klara. I will try to sort this out by basically dropping code-generation from my package, but this will take me a while, it requires massive amount of work.

[chris-nemeth]

If everything else in Klara.jl works fine then it's not worth modifying the package just for GaussianProcesses.jl.

[papamarkou]

What would happen in this case though with GaussianProcesses? Can you see a workaround for your package? The main difficulty is that the methods in Klara can't be encapsulated in a function body anymore as it seems...

[chris-nemeth]

For the time being, I can implement an HMC algorithm for the mcmc() function.

[papamarkou]

Thanks Chris, sounds good, as it's going to take quite some time before I succeed in removing code generation from the whole codebase.

[papamarkou]

Keno provided a solution on discourse, which solved this world age error. Of course, there are many more world age errors to resolve before GaussianProcesses can operate, but my hope is that the rest of the errors are now soluble in a similar way. I will check this out.

[papamarkou]

I updated discourse with one new issue I encountered while trying to convert all the eval() calls to macro invocations. I replicate below the post from discourse to increase visibility as well as the likelihood of finding a solution.

The following works:

macro codegen(y)
  quote
    function (x::$(typeof(y)))
      x+1
    end
  end
end

f = @codegen(2.)
f(3.)

However, if I try to call the macro from inside a function (function g() below), the following fails:

g(y) = @codegen(y)

f = g(2.)
f(3.)

Is there a way to solve this issue by possibly modifying the macro somehow?

[papamarkou]

Update; I removed several world-age errors associated with the main parameter methods and autodiff. There is still quite some work to do to eliminate all world age issues, so I will leave this ticket open till I am done. I am using your older breaking code @chris-nemeth to check when things will work completely from your end.

[papamarkou]

After my last commit, your initial old GP example script runs for MALA without world age errors. I still need to make the analogous changes for all other samplers, streams and nstates before code-generation is eliminated completely, but perhaps about half of the job is done by now, seeing your initial code is operational.

[chris-nemeth]

Thanks. That's great news! I'm away at the monent but I'll look at this when I return.

[papamarkou]

After the latest upgrade to Klara v0.9.1, there is no code generation anymore in Klara in order to avoid world-age issues. This means that this issue has also been resolved, as there is no way you'd encounter any further world age errors. I will now close the issue.

[chris-nemeth]

Thanks for taking this.