@ariddell, I’ve got most of it working! I’m on this branch:

`feature/refactor-services`

Here’s how I’m building it:

```
rm pystan/*.so
python3 setup.py build_ext --inplace
```

And here’s what I’m using to test:

```
import pystan
schools_code = """
data {
int<lower=0> J; // number of schools
real y[J]; // estimated treatment effects
real<lower=0> sigma[J]; // s.e. of effect estimates
}
parameters {
real mu;
real<lower=0> tau;
real eta[J];
}
transformed parameters {
real theta[J];
for (j in 1:J)
theta[j] = mu + tau * eta[j];
}
model {
eta ~ normal(0, 1);
y ~ normal(theta, sigma);
}
"""
schools_dat = {'J': 8,
'y': [28, 8, -3, 7, -1, 1, 18, 12],
'sigma': [15, 10, 16, 11, 9, 11, 10, 18]}
fit = pystan.stan(model_code=schools_code, data=schools_dat, verbose=True, seed=1, chains=1, algorithm = "NUTS", control = dict(adapt_engaged = True, metric = "dense_e"))
## 1. diagnose: FIXME
pystan.stan(fit = fit, data=schools_dat, test_grad = True, chains = 1)
## 2: optimize newton
sm = pystan.StanModel(model_code=schools_code)
op = sm.optimizing(data=schools_dat, algorithm = "Newton")
## 3: optimize bfgs
op = sm.optimizing(data=schools_dat, algorithm = "BFGS")
## 4: optimize lbfgs
op = sm.optimizing(data=schools_dat, algorithm = "LBFGS")
## 5: sample fixed_param: FIXME
fit = pystan.stan(fit = fit, data=schools_dat, algorithm = "Fixed_param")
## 6: sample nuts dense_e adapt_false
fit = pystan.stan(fit = fit, data=schools_dat, algorithm = "NUTS", control = dict(metric = "dense_e", adapt_engaged = False))
## 7: sample nuts dense_e adapt_true
fit = pystan.stan(fit = fit, data=schools_dat, algorithm = "NUTS", control = dict(metric = "dense_e", adapt_engaged = True))
## 8: sample nuts diag_e adapt_false
fit = pystan.stan(fit = fit, data=schools_dat, algorithm = "NUTS", control = dict(metric = "diag_e", adapt_engaged = False))
## 9: sample nuts diag_e adapt_true
fit = pystan.stan(fit = fit, data=schools_dat, algorithm = "NUTS", control = dict(metric = "diag_e", adapt_engaged = True))
## 10: sample nuts unit_e adapt_false
fit = pystan.stan(fit = fit, data=schools_dat, algorithm = "NUTS", control = dict(metric = "unit_e", adapt_engaged = False))
## 11: sample nuts unit_e adapt_true
fit = pystan.stan(fit = fit, data=schools_dat, algorithm = "NUTS", control = dict(metric = "unit_e", adapt_engaged = True))
## 12: sample static dense_e adapt_false
fit = pystan.stan(fit = fit, data=schools_dat, algorithm = "HMC", control = dict(metric = "dense_e", adapt_engaged = False))
## 13: sample static dense_e adapt_true
fit = pystan.stan(fit = fit, data=schools_dat, algorithm = "HMC", control = dict(metric = "dense_e", adapt_engaged = True))
## 14: sample static diag_e adapt_false
fit = pystan.stan(fit = fit, data=schools_dat, algorithm = "HMC", control = dict(metric = "diag_e", adapt_engaged = False))
## 15: sample static diag_e adapt_true
fit = pystan.stan(fit = fit, data=schools_dat, algorithm = "HMC", control = dict(metric = "diag_e", adapt_engaged = True))
## 16: sample static unit_e adapt_false
fit = pystan.stan(fit = fit, data=schools_dat, algorithm = "HMC", control = dict(metric = "unit_e", adapt_engaged = False))
## 17: sample static unit_e adapt_true
fit = pystan.stan(fit = fit, data=schools_dat, algorithm = "HMC", control = dict(metric = "unit_e", adapt_engaged = True))
## 18: variational fullrank
sm.vb(data=schools_dat, algorithm = "fullrank")
## 19: variational meanfield
sm.vb(data=schools_dat, algorithm = "meanfield")
```

I’m having trouble with:

- test_grad. The default arguments aren’t propagating into the c++. I think we saw this before. It runs the code, but the
`epsilon`

and`error`

variables are both 0 instead of 1e-6. - use of sample_file. I tried setting it, but it doesn’t get set. It’s always writing to the temp directory. I don’t think this has anything to do with what I’ve done, but I can’t be sure.

I haven’t tested as thoroughly as RStan, but that’s cause I’m really fumbling around with Python. @ariddell, mind taking a look? Most of the functionality looks ok as far as I can tell.