We derive a graph expansion for the thermal partition function of solvable two-dimensional models with boundaries. This expansion of the integration measure over the virtual particles winding around the time cycle is obtained with the help of the matrix-tree theorem. The free energy is a sum over all connected graphs, which can be either trees or trees with one loop. The generating function for the connected trees satisfies a non-linear integral equation, which is equivalent to the TBA equation. The sum over connected graphs gives the bulk free energy as well as the exact g-functions for the two boundaries. We reproduced the integral formula conjectured by Dorey, Fioravanti, Rim and Tateo, and proved subsequently by Pozsgay. Our method can be extended to the case of non-diagonal bulk scattering and diagonal reflection matrices with proper regularization.