```python
# goal: have sympy do the mechanical substitutions, to double-check the desired relations
# once this is done, this will also make it easier for a human to check (just double-check the definitions), and easier
# to check for arbitrary splittings

from sympy import *
from sympy import init_printing
init_printing()
```


```python
symbolic_traj = dict()
symbolic_traj['x_0'] = symbols('x_0')
symbolic_traj['v_0'] = symbols('v_0')
dt = symbols('dt')
m = symbols('m')
gamma = symbols('gamma')
f = symbols('f')
kT = symbols('kT')
U = symbols('U')

def count_steps(splitting="OVRVO"):
    n_O = sum([step == "O" for step in splitting])
    n_R = sum([step == "R" for step in splitting])
    n_V = sum([step == "V" for step in splitting])
    return n_O, n_R, n_V

def create_variable_names(i):    
    name_x = 'x_{}'.format(i)
    name_v = 'v_{}'.format(i)
    
    symbolic_traj[name_x] = symbols(name_x)
    symbolic_traj[name_v] = symbols(name_v)
    
    x = symbolic_traj['x_{}'.format(i-1)]
    v = symbolic_traj['v_{}'.format(i-1)]
    
    return x, v, name_x, name_v

def apply_R(i, h):
    x, v, name_x, name_v = create_variable_names(i)
    
    symbolic_traj[name_x] = x + v * h
    symbolic_traj[name_v] = v

def apply_V(i, h):
    x, v, name_x, name_v = create_variable_names(i)
    
    symbolic_traj[name_x] = x
    symbolic_traj[name_v] = v + f(x) * h / m
    
def apply_O(i, h):
    a = exp(-gamma * h)
    b = sqrt(1 - exp(-2 * gamma * h))
    
    x, v, name_x, name_v = create_variable_names(i)
    
    symbolic_traj[name_x] = x
    symbolic_traj[name_v] = a * v + b * sqrt(kT / m) * symbols('xi_{}'.format(i))
    # xi_i is an i.i.d. standard normal r.v.
    
def total_energy(x, v):
    return U(x) + (0.5 * m * v**2)

def get_total_energy_at_step(i):
    x = symbolic_traj['x_{}'.format(i)]
    v = symbolic_traj['v_{}'.format(i)]
    return total_energy(x, v)
    
def apply_integrator(i, splitting="OVRVO"):
    heat, W_shad = 0, 0
    
    n_O, n_R, n_V = count_steps(splitting)
    
    for step in splitting:
        if step == "O":
            apply_O(i, dt / n_O)
            heat += get_total_energy_at_step(i) - get_total_energy_at_step(i-1)
        elif step == "R":
            apply_R(i, dt / n_R)
            W_shad += get_total_energy_at_step(i) - get_total_energy_at_step(i-1)
        elif step == "V":
            apply_V(i, dt / n_V)
            W_shad += get_total_energy_at_step(i) - get_total_energy_at_step(i-1)
        i += 1
        
    return heat, W_shad

splitting = "OVRVO"
heat, W_shad = apply_integrator(1, splitting)
delta_E = get_total_energy_at_step(len(splitting)) - get_total_energy_at_step(0)

delta_E
```


```python
W_shad
```


```python
heat
```


```python
delta_E - (W_shad + heat)
```


```python
# also define path probabilities?


# here, let's define the individual step probabilities
def normal_pdf(x, mu, sigma):
    return (1 / (sqrt(2 * pi * sigma**2))) * exp(-(x - mu)**2 / (2 * sigma**2))

def R_prob(from_x, from_v, to_x, to_v, h):
    '''deterministic position update'''
    if (from_v == to_v) and (to_x == (from_x + from_v * h)):
        return 1
    else:
        return 0

def V_prob(from_x, from_v, to_x, to_v, h):
    if (from_x == to_x) and (to_v == (from_v + f(from_x) * h / m)):
        return 1
    else:
        return 0

def O_prob(from_x, from_v, to_x, to_v, h):
    ''' need to double check this, probably dropped a 2 somewhere...'''
    a = exp(-gamma * h) # okay
    b = sqrt(1 - exp(-2 * gamma * h)) # okay
    
    sigma = b * sqrt(kT / m) # double-check this! is this the definition of sigma or sigma^2?
    mu = a * from_v # okay
    
    if (from_x == to_x):
        return normal_pdf(to_v, mu, sigma) 
    else:
        return 0

# for example, what's the probability of a small v perturbation
O_prob(0, 0, 0, 0.1, dt / 2) 
```


```python
O_prob(symbolic_traj['x_0'], symbolic_traj['v_0'], symbolic_traj['x_1'], symbolic_traj['v_1'], dt / 2)
```


```python
def forward_path_probability(splitting="OVRVO"):
    path_prob = 1
    n_O, n_R, n_V = count_steps(splitting)
    
    
    
    
    for i, step in enumerate(splitting):
        x_0, v_0 = symbolic_traj['x_{}'.format(i)], symbolic_traj['v_{}'.format(i)]
        x_1, v_1 = symbolic_traj['x_{}'.format(i+1)], symbolic_traj['v_{}'.format(i+1)]
        
        if step == "O":
            step_prob = O_prob(x_0, v_0, x_1, v_1,
                               dt / n_O)
        elif step == "R":
            step_prob = R_prob(x_0, v_0, x_1, v_1,
                               dt / n_R)
        elif step == "V":
            step_prob = V_prob(x_0, v_0, x_1, v_1,
                               dt / n_V)

        path_prob = path_prob * step_prob
    return path_prob
```


```python
forward_path_probability("OVRVO")
```


```python
def reverse_path_probability(splitting="OVRVO"):
    '''same as above, just reverse order of trajectory and steps appropriately?'''
    
    path_prob = 1
    n_O, n_R, n_V = count_steps(splitting)

    for i in range(len(splitting))[::-1]:
        step = splitting[i]
        
        
        x_0, v_0 = symbolic_traj['x_{}'.format(i+1)], symbolic_traj['v_{}'.format(i+1)]
        x_1, v_1 = symbolic_traj['x_{}'.format(i)], symbolic_traj['v_{}'.format(i)]
        
        if step == "O":
            step_prob = O_prob(x_0, v_0, x_1, v_1,
                               dt / n_O)
        elif step == "R":
            step_prob = R_prob(x_0, v_0, x_1, v_1,
                               dt / n_R)
        elif step == "V":
            step_prob = V_prob(x_0, v_0, x_1, v_1,
                               dt / n_V)

        path_prob = path_prob * step_prob
    return path_prob
```


```python
def CFT_definition_of_work(splitting="OVRVO"):
    return - ln(reverse_path_probability(splitting) / forward_path_probability(splitting))
```


```python
CFT_definition_of_work("OVRVO")
```


```python
reverse_path_probability()
```


```python
forward_path_probability()
```


```python
# hmm that doesn't look quite right
```

a6ef330526a753e41d7d7a1c09c54974e28e80bc
ipynb
Jupyter Notebook
Using sympy to compute relative path action.ipynb
jchodera/maxentile-notebooks
6e8ca4e3d9dbd1623ea926395d06740a30d9111d
MIT
Using sympy to compute relative path action.ipynb
jchodera/maxentile-notebooks
6e8ca4e3d9dbd1623ea926395d06740a30d9111d
MIT
2018-06-10T12:21:10.000Z
2018-06-10T14:42:45.000Z
Using sympy to compute relative path action.ipynb
jchodera/maxentile-notebooks
6e8ca4e3d9dbd1623ea926395d06740a30d9111d
MIT
2018-06-10T12:14:55.000Z
2018-06-10T12:14:55.000Z
Qwen/Qwen-72B
1. YES
2. YES