Bounding JD in terms of TV

import sys, re, math, os, json, random
import numpy as np

def tvd(x, y):
    return 0.5 * sum(abs(np.array(x) - np.array(y)))

def jpd(x, y):
    ret = 0.0
    for i in range(len(x)):
        if x[i] == 0 or y[i] == 0:
            continue
        ret0 = 0.0
        for j in range(len(x)):
            ret0 += max(x[j] / x[i], y[j] / y[i])
        ret += 1.0 / ret0
    return 1.0 - ret


def simplex_sample(n):
    x = np.random.exponential(scale=10.0, size=n)
    return x / sum(x)


record1 = -100
record2 = 100
for _ in range(1000000):
    x = simplex_sample(3)
    y = simplex_sample(3)
    tv = tvd(x, y)
    jp = jpd(x, y)
    ratio1 = jp / (tv / (1 + tv))
    ratio2 = jp / tv
    if ratio1 > record1:
        print('upper bound: %0.10f' % ratio1)
        record1 = ratio1
    if ratio2 < record2:
        print('lower bound: %0.10f' % ratio2)
        record2 = ratio2