Let $\mu$ and $\nu$ be two probability measures on the real line with finite $p$ moment for $p\in [1, \infty)$. I am interested in the functions $f:\mathbb{R}\rightarrow \mathbb{R}$ and $g:\mathbb{R}\rightarrow \mathbb{R}$ that solve the dual of the optimal transport problem $$ OT(\mu, \nu) = \max\left\{ \int f\,d\mu + \int g\,d\nu,\, \Big|\,f\in L^1(\mu), \, g\in L^1(\nu) \, \text{and}\, f(x)+g(y)\leq h(x-y)\, \right\}, $$ for an arbitrary convex function $h:\mathbb{R}\rightarrow \mathbb{R}$. In particular I would be curious to know how $f$ (and $g$) look like if $\mu$ and $\nu$ are discrete measures (or in general not absolutely continuous). Thanks in advance for any help :)
Update: If $\mu$ and $\nu$ are discrete then, if we set $f(x_1)=0$, one of the dual solutions is given by $$ f(x_{i}) = \sum_{t=1}^{i-1} h(G^{-1}(F(x_t)) - x_{t+1}) - h(G^{-1}(F(x_t)), x_{t}) \quad \text{with} \quad i = 2, \dots, n, $$
where $F$ is the cdf of $\mu$ and $G^{-1}$ is the quantile function of $\nu$. This works even if there is no transport map between $\mu$ and $\nu$. Moreover if $\mu$ and $\nu$ are continuous and have a connected support and $h$ is differentialble, then: $$ f(x) = - \int_{-\infty}^x h'(G^{-1}(F(t)) - t)\, dt. $$ (More on the motivation of this will come soon in the form of an answer).
Assume $\mu$ is discrete and has support points $x_1, \dots, x_n\in \mathbb{R}$. Let $F$ be the c.d.f. of $\mu$ and $F^{-1}$and $G^{-1}$ be the quantile functions of $\mu$ and $\nu$ respectively. It is then known that one optimal transport plan is given by the distribution of $(F^{-1}(U), G^{-1}(U))\sim \pi^*$ where $U\sim \mathrm{Unif}[0, 1]$. Consider now the set $$\Gamma = \bigcup_{i=1}^n \left(\{x_i\}\times [G^{-1}\circ F(x_{i-1}), G^{-1}\circ F(x_i)]\right).$$ Since $G^{-1}$ and $F$ are increasing we have that for all $(x_1, y_1), (x_2, y_2)\in \Gamma$, it holds that, if $x_1 < x_2$ then $y_1\leq y_2$. As one can easily prove, this is equivalent to saying that the set is c-cyclically monotone. Moreover, $\Gamma$ is closed and $\mathrm{supp}\, \pi^* = (F^{-1}, G^{-1})[0, 1] \subset \Gamma$.
By Step 3 of Theorem 5.10 in Optimal Transport old and new by Villani, we then know that there exists a c-convex function $f$ such that its c-subdifferential satisfies $\Gamma \subset \partial_c f$. To now constuct $f$ explicitly we can use recursively the characterisation of the subdifferential $$(x, y) \in \partial_c f \quad \text{if and only if}\quad f(x) + f^c(y) = c(x, y),$$ where $f^c$ is the c-transform of $f$. We now set $f(x_1)=0$. Then, since $(x_1, G^{-1}\circ F(x_1))\in \partial_c f$ we find that $$f^c(G^{-1}\circ F(x_1))=c(x_1, y_1).$$ Proceding by induction, let us assume that we do know the values of $f(x_i)$ and $f^c(G^{-1}\circ F(x_i))$. Then we can use that $(x_{i+1}, G^{-1}\circ F(x_i)) \in \Gamma \subset \partial_c f$ to derive: $$f(x_{i+1}) = c(x_{i+1}, G^{-1}\circ F(x_i)) - f^c(G^{-1}\circ F(x_i))$$ and further since $(x_{i+1}, G^{-1}\circ F(x_{i+1})) \in \Gamma \subset \partial_c f$, $$f^c(G^{-1}\circ F (x_{i+1})) = c(x_{i+1}, G^{-1}\circ F(x_{i+1})) - f(x_{i+1}).$$ By induction we then get $$ f(x_{i}) = \sum_{t=1}^{i-1} h(G^{-1}(F(x_t))- x_{t+1}) - h(G^{-1}(F(x_t))- x_{t}) \quad \text{with} \quad i = 2, \dots, n. $$