Unable to include functions from other file

[itsdfish]

Hi Rob-

I am having trouble including functions from another file, which worked at one point. I'm not sure whether it is my execution or something changed on the CmdStan side. Suppose you have a folder called my_model, which contains three files: Run_model.jl, my_functions.stan, and bernoulli.stan. What I want to do is import reusable functions from my_functions.stan into bernoulli.stan. Here is what the files look like:

my_functions.stan:

my_function(){
    print("I'm working");
}

bernoulli.stan:

functions {
#include ../my_functions.stan
}

data { 
  int<lower=1> N; 
  int<lower=0,upper=1> y[N];
} 
parameters {
  real<lower=0,upper=1> theta;
} 
model {
  my_function();
  theta ~ beta(1,1);
  y ~ bernoulli(theta);
}

Run_model.jl:

using CmdStan,Distributions,Random
ProjDir = @__DIR__
cd(ProjDir)

stream = open("bernoulli.stan","r")
Model = read(stream,String)
close(stream)

stanmodel = Stanmodel(Sample(save_warmup=true, num_warmup=1000,
  num_samples=2000, thin=1), name="bernoulli", model=Model,
  printsummary=false)

observeddata = Dict("N" => 10, "y" => [0, 1, 0, 1, 0, 0, 0, 0, 0, 1])

rc, chn, cnames = stan(stanmodel, observeddata, ProjDir)

I recieve an error along these lines (although the specific directories may differ from the simplified example above):

could not find include file  in the following directories:

make: *** No rule to make target '/home/dfish/.julia/dev/StanACTR/src/tmp/bernoulli.hpp', needed by '/home/dfish/.julia/dev/StanACTR/src/tmp/bernoulli'.  Stop.

I've tried various directories and syntax changes, but to no avail. Do you know how to make this work? Thanks!

[itsdfish]

P.S. Here is a forum thread where including files is discussed.

[ghost]

I cannot replicate this exactly, and it works for me if the last line is modified as:

rc, chn, cnames = stan(stanmodel, observeddata, ProjDir, CmdStanDir=CMDSTAN_HOME)

[itsdfish]

Strange. I am able to replicate this on Ubuntu 18.04 and Mac 10.11. I receive the following error in the build log:

could not find include file ../my_functions.stan in the following directories:

make: *** [/Users/christopher.fisher/.julia/dev/Test/tmp/bernoulli.hpp] Error 255

Changing to:

rc, chn, cnames = stan(stanmodel, observeddata, ProjDir, CmdStanDir=CMDSTAN_HOME)

Did not help in my case. It seems like the problem is that it always checks the tmp folder for my_functions.stan.

[goedman]

Definitely see the same problem as Chris is seeing. Oliver, are you on Windows?

[ghost]

Linux and Julia 1.3.

[goedman]

Thanks Oliver (and also for jumping in!) That makes it even stranger as Chris also sees it on Ubuntu.

Given that the thread above is all based on R I assume Ben's answer relates to R's #include.

the build-log also shown an empty directory:

could not find include file my_functions.stan; in the following directories:

make: *** [/Users/rob/.julia/dev/StanSample/examples/Include_Example/tmp/bernoulli.hpp] Error 255

The last make complaint is logical.

[itsdfish]

That is very strange that it works only on Linux. It is probably not relevant, but my Julia versions are 1.1 on Mac 10.11 and 1.1.1 on Ubuntu 18.04.

[goedman]

I don't think it is the Julia version.

It doesn't quite work, but Chris could you take a look at:

using CmdStan, Distributions, Random

ProjDir = @__DIR__
cd(ProjDir)

stream = open("bernoulli.stan", "r")
initial_model = read(stream, String)
close(stream)

stream = open("my_functions.stan", "r")
funcs = read(stream, String)
close(stream)

model = "functions{\n$(funcs)}\n"*initial_model
println(model
)
stanmodel = Stanmodel(Sample(save_warmup=true, num_warmup=1000,
  num_samples=2000, thin=1), name="bernoulli", model=model,
  printsummary=false)

observeddata = Dict("N" => 10, "y" => [0, 1, 0, 1, 0, 0, 0, 0, 0, 1])

rc, chn, cnames = stan(stanmodel, observeddata, ProjDir, CmdStanDir=CMDSTAN_HOME)

where bernoulli.stan is:

data { 
  int<lower=1> N; 
  int<lower=0,upper=1> y[N];
} 
parameters {
  real<lower=0,upper=1> theta;
} 
model {
  my_function();
  theta ~ beta(1,1);
  y ~ bernoulli(theta);
}

and my_functions.stan contains:

  my_function(){
      print("I'm working");
  }

Maybe you see the error I'm making here?

[itsdfish]

Let me think about this. You are adding a function block twice. What it needs is a way to insert my_function into the function block (and possibly remove the include statement).

[goedman]

Don't think I'm adding a second function block. Took that out of the starting bernoulli.stan.

[goedman]

Isn’t that what I am achieving by using string interpolation?

[itsdfish]

Here is what I find:

model = "functions{\n$(funcs)}\n"*initial_model
println(model)

functions{
my_function(){
    print("I'm working");
}
}
functions {
#include ../../my_functions.stan
}
...

What I think it needs is (1) string split into three parts: "functions {", "include path", " remaining model body (2) insert the functions between the first and third parts from the first step, and (3) some sort of logic that makes the splitting robust to white space.

Unfortunately, I got stuck at (1):

str=split(initial_model,"functions {")
 str
  2-element Array{SubString{String},1}:
 ""                                                                                                                                                                                                                                 
 "#include ../../my_functions.stan\n}\n\ndata {\n  int<lower=1> N;\n  int<lower=0,upper=1> y[N];\n}\nparameters {\n  real<lower=0,upper=1> theta;\n}\nmodel {\n  my_function();\n  theta ~ beta(1,1);\n  y ~ bernoulli(theta);\n}\n"

For some reason its not splitting properly. I'll see if I can help more later today.

[goedman]

The 2nd functions block must be part of your initial-model to get your final model.

[itsdfish]

Just for my edification, is the process of including a file something that you handle with your own parser or is that something on the Stan side? I guess I'm confused because it seems to be working via CmdStan on Oliver's system.

[ghost]

julia> using CmdStan, Distributions, Random

julia> ProjDir = @__DIR__
"/home/oliver"

julia> cd(ProjDir)

julia> stream = open("bernoulli.stan", "r")
IOStream(<file bernoulli.stan>)

julia> initial_model = read(stream, String)
"functions {\n  #include ../my_functions.stan\n}\ndata {\n  int<lower=1> N;\n  int<lower=0,upper=1> y[N];\n}\nparameters {\n  real<lower=0,upper=1> theta;\n}\nmodel {\n  my_function();\n  theta ~ Beta(1,1);\n  y ~ Bernoulli(theta);\n}\n"

julia> close(stream)

julia> stream = open("my_functions.stan", "r")
IOStream(<file my_functions.stan>)

julia> funcs = read(stream, String)
"my_function(){\n print(\"I'm working\");\n}\n"

julia> close(stream)

julia> model = "functions{\n$(funcs)}\n"*initial_model
"functions{\nmy_function(){\n print(\"I'm working\");\n}\n}\nfunctions {\n  #include ../my_functions.stan\n}\ndata {\n  int<lower=1> N;\n  int<lower=0,upper=1> y[N];\n}\nparameters {\n  real<lower=0,upper=1> theta;\n}\nmodel {\n  my_function();\n  theta ~ Beta(1,1);\n  y ~ Bernoulli(theta);\n}\n"

julia> println(model
       )
functions{
my_function(){
 print("I'm working");
}
}
functions {
  #include ../my_functions.stan
}
data {
  int<lower=1> N;
  int<lower=0,upper=1> y[N];
}
parameters {
  real<lower=0,upper=1> theta;
}
model {
  my_function();
  theta ~ Beta(1,1);
  y ~ Bernoulli(theta);
}

julia> stanmodel = Stanmodel(Sample(save_warmup=true, num_warmup=1000,
         num_samples=2000, thin=1), name="bernoulli", model=model,
         printsummary=false)

File /home/oliver/tmp/bernoulli.stan will be updated.

  name =                    "bernoulli"
  nchains =                 4
  num_samples =             1000
  num_warmup =                1000
  thin =                    1
  monitors =                String[]
  model_file =              "bernoulli.stan"
  data_file =               ""
  output =                  Output()
    file =                    ""
    diagnostics_file =        ""
    refresh =                 100
  pdir =                   "/home/oliver"
  tmpdir =                 "/home/oliver/tmp"
  output_format =           :mcmcchains
  method =                  Sample()
    num_samples =             2000
    num_warmup =              1000
    save_warmup =             true
    thin =                    1
    algorithm =               HMC()
      engine =                  NUTS()
        max_depth =               10
      metric =                  CmdStan.diag_e
      stepsize =                1.0
      stepsize_jitter =         1.0
    adapt =                   Adapt()
      gamma =                   0.05
      delta =                   0.8
      kappa =                   0.75
      t0 =                      10.0
      init_buffer =             75
      term_buffer =             50
      window =                  25

julia> observeddata = Dict("N" => 10, "y" => [0, 1, 0, 1, 0, 0, 0, 0, 0, 1])
Dict{String,Any} with 2 entries:
  "N" => 10
  "y" => [0, 1, 0, 1, 0, 0, 0, 0, 0, 1]

julia> rc, chn, cnames = stan(stanmodel, observeddata, ProjDir, CmdStanDir=CMDSTAN_HOME)
(0, Object of type Chains, with data of type 3000×8×4 Array{Float64,3}

Iterations        = 1:3000
Thinning interval = 1
Chains            = 1, 2, 3, 4
Samples per chain = 3000
internals         = lp__, accept_stat__, stepsize__, treedepth__, n_leapfrog__, divergent__, energy__
parameters        = theta

2-element Array{ChainDataFrame,1}

Summary Statistics

│ Row │ parameters │ mean     │ std     │ naive_se   │ mcse       │ ess     │ r_hat   │
│     │ Symbol     │ Float64  │ Float64 │ Float64    │ Float64    │ Any     │ Any     │
├─────┼────────────┼──────────┼─────────┼────────────┼────────────┼─────────┼─────────┤
│ 1   │ theta      │ 0.332135 │ 0.12908 │ 0.00117833 │ 0.00201116 │ 3844.34 │ 1.00021 │

Quantiles

│ Row │ parameters │ 2.5%     │ 25.0%    │ 50.0%    │ 75.0%    │ 97.5%   │
│     │ Symbol     │ Float64  │ Float64  │ Float64  │ Float64  │ Float64 │
├─────┼────────────┼──────────┼──────────┼──────────┼──────────┼─────────┤
│ 1   │ theta      │ 0.108581 │ 0.237781 │ 0.324441 │ 0.416432 │ 0.60805 │
, ["lp__", "accept_stat__", "stepsize__", "treedepth__", "n_leapfrog__", "divergent__", "energy__", "theta"])

[ghost]

julia> versioninfo()
Julia Version 1.3.0-alpha.0
Commit 6c11e7c2c4 (2019-07-23 01:46 UTC)
Platform Info:
  OS: Linux (x86_64-pc-linux-gnu)
  CPU: Intel(R) Core(TM) i5-9600K CPU @ 3.70GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-6.0.1 (ORCJIT, skylake)
Environment:
  JULIA_CMDSTAN_HOME = /home/oliver/Library/cmdstan

[ghost]

oliver@debian:~/Library/cmdstan$ make build

--- CmdStan v2.20.0 built ---

[goedman]

In your original question, Stan’s cmdstan binary tried to handle that. With string interpolation Julia handles it and passes the generated bernoulli.stan in tmp to Stan.

[goedman]

@ode33 uses 2 function blocks. Need to try that. I’m also on CmdStan 2.20.0 which helps big time with Stan’s compilation time.

[itsdfish]

Although it runs on Oliver's machine, I still get the same error. Oliver, can you verify that the bernoulli.stan file in tmp updated properly. Sometimes it does not always update.

One question I have is whether it will run code defined in both function blocks?

[ghost]

Stan files are not supported, here is what's in /tmp/: bernoulli.txt

[itsdfish]

Thanks for confirming Oliver. Very odd. That doesn't work on my system. One reason might be that it still has the include statement in the original function block.

Does it work if you change your original bernoulli file (not in tmp) to the following:

functions {
#include ../../my_functions.stan
my_other_function(){
  print("printing from original function block");
  }
}

data {
  int<lower=1> N;
  int<lower=0,upper=1> y[N];
}
parameters {
  real<lower=0,upper=1> theta;
}
model {
  my_function();
  my_other_function();
  theta ~ beta(1,1);
  y ~ bernoulli(theta);
}

I'm just wondering whether Stan uses code from both function blocks.

[ghost]

The updated model using the code from the original example works for me.

[itsdfish]

Thanks for looking into that. My best guess is that #include ../../my_functions.stan needs to be removed in order for it to work on my system and Rob's. I'll look into this once I have a free moment.

[itsdfish]

Rob-

I am unable to run code with two function blocks:

bernoulli.stan:

functions {
  int my_other_function(){
    print("my other function");
    return 0;
    }
}

functions {
  int my_function(){
    print("my function");
    return 0;
    }
}

data {
  int<lower=1> N;
  int<lower=0,upper=1> y[N];
}
parameters {
  real<lower=0,upper=1> theta;
}
model {
  //my_function();
  //my_other_function();
  theta ~ beta(1,1);
  y ~ bernoulli(theta);
}

Error:

PARSER FAILED TO PARSE INPUT COMPLETELY
STOPPED AT LINE 8: 
functions {
  int my_function(){
    print("my function");
    return 0;
    }
}

data {
  int<lower=1> N;
  int<lower=0,upper=1> y[N];
}
parameters {
  real<lower=0,upper=1> theta;
}
model {
  //my_function();
  //my_other_function();
  theta ~ beta(1,1);
  y ~ bernoulli(theta);
}

make: *** [/Users/christopher.fisher/.julia/dev/Test/tmp/bernoulli.hpp] Error 253

[goedman]

Chris, below include_example.jl now works. I've changed my_functions.jl.

include_example.jl:

using CmdStan, Distributions, Random

ProjDir = @__DIR__
cd(ProjDir)

stream = open("bernoulli.stan", "r")
initial_model = read(stream, String)
close(stream)

stream = open("my_functions.stan", "r")
funcs = read(stream, String)
close(stream)

model = "functions{\n$(funcs)}\n"*initial_model
println(model)

stanmodel = Stanmodel(Sample(save_warmup=true, num_warmup=1000,
  num_samples=2000, thin=1), name="bernoulli", model=model,
  printsummary=false)

observeddata = Dict("N" => 10, "y" => [0, 1, 0, 1, 0, 0, 0, 0, 0, 1])

rc, chn, cnames = stan(stanmodel, observeddata, ProjDir, CmdStanDir=CMDSTAN_HOME)

bernoulli.stan:

data { 
  int<lower=1> N; 
  int<lower=0,upper=1> y[N];
} 
parameters {
  real<lower=0,upper=1> theta;
} 
model {
  ;
  theta ~ beta(my_function(),1);
  y ~ bernoulli(theta);
}

my_functions.stan:

  real my_function(){
      real x = 1.0;
      return x;
  }

[goedman]

With 2 function blocks I get the same error. Right now I am trying to understand if there are restrictions on functions Stan can include or if we have to escape the "s in print("I'm working");.

[itsdfish]

Thanks Rob. That code works on my system. It looks like the missing piece is integrating functions defined in my_functions.stan with the functions defined in the function block of bernoulli.stan. Do you have any clever ideas? String split was not working in my initial attempts.

[goedman]

Isn't that what I'm doing with the string interpolation followed by the string concatenation?

[goedman]

I'm envisioning something like:

function full_model(base_model, functions)
  stream = open(base_model, "r")
  initial_model = read(stream, String)
  close(stream)

  stream = open(functions, "r")
  funcs = read(stream, String)
  close(stream)

  model = "functions{\n$(funcs)}\n"*initial_model
  model
end

model = full_model("bernoulli.stan", "my_functions.stan")

[itsdfish]

I'm afraid not. My apologies if I failed to articulate the problem clearly. Below, I attached a working example based on your most recent code. The file called intended_model.stan shows what the model should look like when everything is integrated properly. What your code currently produces is something like this, with two function blocks:

functions{
//Sharable functions across models
real general_function(){
    real x = 1.0;
    return x;
}
}
functions {
  real model_specific_function(){
      real x = 1.0;
      return x;
  }
}

data {
  int<lower=1> N;
  int<lower=0,upper=1> y[N];
}
parameters {
  real<lower=0,upper=1> theta;
}
model {
  model_specific_function();
  theta ~ beta(general_function(),1);
  y ~ bernoulli(theta);
}

Test.zip

[goedman]

Aah, you need to mix in shared and 'local' function?

That could be done in a similar way if it is ok not to include them in your initial_model (bernoulli.stan in Test.zip), something like:

function full_model(base_model, functions)
  initial_model = open(f->read(f, String), base_model)
  funcs = open(f->read(f, String), functions)

  model = "functions{\n$(funcs)}\n"*initial_model
  model
end

function full_model(base_model, shared_functions, local_functions)
  initial_model = open(f->read(f, String), base_model)
  shared_funcs = open(f->read(f, String), shared_functions)
  local_funcs = open(f->read(f, String), local_functions)

  model = "functions{\n$(shared_funcs)\n$(local_funcs)}\n"*initial_model
  model
end

If you really want them in the base_model we could use String.split.

[goedman]

An example:

using CmdStan, Distributions, Random

ProjDir = @__DIR__
cd(ProjDir)

function full_model(base_model, shared_functions, local_functions)
  initial_model = open(f->read(f, String), base_model)
  shared_funcs = open(f->read(f, String), shared_functions)
  local_funcs = open(f->read(f, String), local_functions)

  model = "functions{\n$(shared_funcs)$(local_funcs)}\n"*initial_model
  model
end

model2 = full_model("bernoulli2.stan", "shared_funcs.stan", "local_funcs.stan")

stanmodel = Stanmodel(Sample(save_warmup=true, num_warmup=1000,
  num_samples=2000, thin=1), name="bernoulli2", model=model2,
  printsummary=false)

observeddata = Dict("N" => 10, "y" => [0, 1, 0, 1, 0, 0, 0, 0, 0, 1])

rc, chn, cnames = stan(stanmodel, observeddata, ProjDir)

if rc == 0
  read_summary(stanmodel)
end

where bernoulli.stan:

data {
  int<lower=1> N;
  int<lower=0,upper=1> y[N];
}
parameters {
  real<lower=0,upper=1> theta;
}
model {
  model_specific_function();
  theta ~ beta(my_function(),1);
  y ~ bernoulli(theta);
}

and local_funcs.stan:

  void model_specific_function(){
      real x = 1.0;
      return;
  }

and shared_funs.stan:

  real my_function(){
      real x = 1.0;
      return x;
  }

[goedman]

In this setup full_model used multiple dispatch with 2 or 3 arguments.

[itsdfish]

Aah, you need to mix in shared and 'local' function?

Yup. Sorry about the miscommunication.

If you really want them in the base_model we could use String.split.

I don't think that is necessary, but it might be convenient to have the local functions with the other blocks and to define models in the typical fashion. I will see if I can modify your example with string split. I think the biggest problem might be dealing white space. It might be more trouble than its worth. Thanks for your help and sorry for the miscommunication!

[goedman]

There are ways of course, but it might prove error prone.

Every time you ask a question I learn something (how you and other folks use the package!).

Note: I'm adding the example to both CmdStan and StanSample. Now I'll look in why Oliver's seem to work!

[itsdfish]

Yeah. # include does not appear to be a well-documented feature. I imagine that most people do not use it too often because they are doing standard GLMs. I have some ideas for parsing the model text. If there is no way to solve the problem on your and and I get the parser to work, would you be interested in adding to CmdStan?

[goedman]

Absolutely.

[goedman]

What would be pretty cool it to parse for an #include line in base_stan and string interpolate the file to be included when creating base_stan in tmp.

[itsdfish]

Awesome. Well, I have a working prototype. It seems to work under a variety of conditions. It simply inserts the functions at the location of the #include myfile statement. If this seems suitable for you, feel free to integrate it with your code, presumably immediately after receiving the model string.

Parse_Example.zip

[itsdfish]

Hi Rob. I integrated my code with yours and added tests. If you can add me as a collaborator I can add my branch and open a pull request when I get home today.

[goedman]

Hi Chris, I realized what we discussed and came up with is probably not the best way of including this feature.

I worked yesterday a bit more on an approach that would do this transparently when creating the .stan file in the tmp directory and is conform the Stan language spec. Base part is working, need to add one little part (removing comments at the end of the #include line).

Once that is finished I’ll send you what I came up with and we can discuss what we like best. The tests certainly need to be included.

Did send you an invite!

[itsdfish]

Thanks for sending the invite. I don't have the tests with me, but here is what I came up with for the stanmodel.jl file. It seems to work pretty well. Feel free to choose which ever approach works best.

stanmodel.jl.zip

[itsdfish]

I've been using it all morning and it seems to work well. The main limitation is that it only works for a single include statement. One simple way to generalize the code would be to wrap addFunctions in a while loop and iteratively include files until no more include statements are found. This will work because addFunctions finds the first instance of # include

[goedman]

This is what I have right now:

function parse_and_interpolate(model::AbstractString)
  newmodel = ""
  lines = split(model, "\n")
  for l in lines
    ls = strip(l)
    replace_strings = findall("#include", ls)
    if length(replace_strings) == 1
      for r in replace_strings
        ls = ls[r[end]+1:end]
        func = open(f -> read(f, String), strip(ls))
        newmodel *= "   "*func*"\n"
      end
    else
      if length(replace_strings) > 1
        error("Inproper #includes")
      else
        newmodel *= l*"\n"
      end
    end
  end
  newmodel
end

model = parse_and_interpolate(base_model)
println(model)

This handles multiple include lines. No yet comments at the end of #include lines (// comment).

new model is what update_model_file() will use.

[itsdfish]

Excellent. Its good to know we have some viable solutions.

[goedman]

HI Chris, I've now updated CmdStan master branch. It contains 2 new example directories vs. the last released version. The Parse_and_Interpolate_Example is now handled in update_model_file(). A couple of cases I have tested are in example.jl (some commented out, can't have same function defined twice). If the include file includes itself, cmdstan should catch that. I thing Julia's string interpolation in the include line is nice as well.

Let me know if all your test cases are covered.

For now I have also left the example Include_Example in the repo. It is another, more manual way of doing it and not according to the Stan language spec.

We could have a 3rd example along your work provided we can dispatch to the right update_model_file().

Haven't included the tests in runtests.jl yet, first wanted to know what route seems best.

[itsdfish]

Thanks Rob. I pushed my changes to the branch include_file.

Unfortunately, I receive the following error when I try to run your parse example:

MethodError: no method matching findall(::String, ::String)
Closest candidates are:
  findall(::String, !Matched::LibGit2.GitIndex) at /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.1/LibGit2/src/index.jl:192
  findall(!Matched::Function, ::Any) at array.jl:1966
  findall(::Any) at array.jl:2007
parse_and_interpolate(::SubString{String}) at update_model_file.jl:48
update_model_file(::String, ::SubString{String}) at update_model_file.jl:27
#Stanmodel#3(::String, ::Int64, ::Int64, ::Int64, ::Int64, ::String, ::Array{String,1}, ::Array{Dict{String,Any},1}, ::CmdStan.Random, ::Array{Dict{String,Any},1}, ::CmdStan.Output, ::Bool, ::String, ::String, ::Symbol, ::Type, ::Sample) at stanmodel.jl:164
(::getfield(Core, Symbol("#kw#Type")))(::NamedTuple{(:name, :model, :printsummary, :tmpdir),Tuple{String,String,Bool,String}}, ::Type{Stanmodel}, ::Sample) at none:0
top-level scope at none:0

I wonder if it is due to differences in versions of Julia?

[goedman]

Yes, to my BIG surprise that extension of findall() is only available in Julia 1.3!!!

[itsdfish]

The only other thing I might consider adding to the test is a second include statement to ensure that it works with multiple files.