Hi guys,

I work with numerical sampling methods that evolve some vectorial parameter \boldsymbol{x} and associated momentum \boldsymbol{p} via a force function \boldsymbol{F}(\boldsymbol{x}) that comes from an underlying probability density.

The method runs for N iterations, updating parameters, momenta, and forces in each step, and computes the moving average of certain observables \phi(\boldsymbol{x}, \boldsymbol{p}) any n iterations via \langle \phi \rangle = \frac{1}{k} \sum_{i=1}^{k} \phi(\boldsymbol{x}_i, \boldsymbol{p}_i), where k is the current iteration count.

Below is a spaghetti-code example. I was wondering how one could use modules to modularize it properly in order to increase generality.

There are three main components to a simulation like that:

- The choice of the sampling method (different methods evolve the system (\boldsymbol{x}, \boldsymbol{p}) in different ways.
- The choice of the observables to collect, i.e. the functions \phi.
- The force function to use, i.e. the underlying probability model.

Therefore, I thought one could implement three different modules, one for the samplers, one for the observables, one for the forces. The forces and the observables are independent of one another, but the sampler method needs to â€śseeâ€ť both the force function and the observables.

In C++, I had a setup like this using classes and inheritance, but I am not sure about the Julia way.

Any tips would be great!

```
using Random
# force function
function compute_force!(x::Vector{<:Real}, force::Vector{<:Real})
#force[1] = ...
#force[2] = ...
# ....
end
# simulation parameters, eg.
maxIter = Int(1e7) # no. of iterations
dt = 0.01 # step size
init_cond = [0.0, 0.0] # initial parameter, momentum
output_file_name = "output.csv" # print observables to file with that name
t_meas = 5 # collect sample any t_meas iterations
# array initialization
x = [init_cond[1]] # initialize x and p
p = [init_cond[2]]
force = [0.0]
compute_force!(x, force) # get first force
# specify observables to measure, in this case just p^2
p_sq_sum = 0.0
no_samples = 0
res_vec_size = maxIter Ă· t_meas + 1 # size of the results vector
p_sq_results = zeros( res_vec_size ) # going to store the evolving time average of the observable p^2
# define rng and some useful constants
rand_vec = zeros(size(p))
rng = Xoshiro(1234)
# ....other numerical constants....
# run sampler
for i in 0:maxIter
# take a measurement
if i%t_meas == 0
p_sq_sum += p'*p
no_samples += 1
p_sq_results[1+iĂ·t_meas] = p_sq_sum / no_samples
end
### evolve (x,p) via some dynamics governed by force and random number generator
# x=...
# p=...
compute_force!(x, force) # update force
# x=...
# p=...
end
# print result to .csv file, in this case just the array p_sq_results
```