Is there an algorithm to compute efficiently the dessin d'enfant from a Belyi pair? Let $(X,f)$ be a Belyi pair, i.e. a Riemann surface $X$ together with a morphism $f: X \to \mathbb{P}^1$, ramified only in $0,1, \infty$. Grothendieck's dessin d'enfant is the pre-image $G$ of the interval $[0,1]$ considered as a graph embedded in $X$ (a dessin d'enfant has more information than that, but let us concentrate on just the graph for the moment). What is the algorithm to compute $G$?
What I mean by an algorithm is some efficient procedure that works in terms of algebraic data used to define $X$ and $f$. We may assume $X$ projective, so let $X$ be given by its homogeneous coordinate ring $R^\bullet$, specified by generators and relations, and $f$ specified by elements of a graded $R^\bullet$-module $M^\bullet$ corresponding to the line bundle $f^*(O(-1))$.
Is there an algorithm that computes the adjacency matrix of $G$ in terms of these data?
update: a dessin d'enfant also contains splitting of the set of vertices into two classes and a cyclic ordering of edges incident to each vertex (which arises from the monodromy action) which can be used to define an action of the free group on two generators on the edges of the dessin. The stabiliser of the action is a subgroup of finite index, and in fact, a conjugacy class of a group of finite index in $F_2$ defines the isomorphism type of a dessin d'enfant.
Therefore my question can also be reformulated accordingly: how to compute efficiently (conjugacy class of) a finite index subgroup of $F_2$ corresponding to a Belyi pair $(X,f)$?
 A: This question is addressed in the recent preprint arXiv:13112529, in section 7. This is a survey of computing Belyi maps from the designs, but section 7 addresses the inverse problem.
A: I'd look into numerical homotopy software, such as Bertini or PHCPack. Numerical homotopy software attempts to solve problems of the following sort: Suppose we have a family of polynomial equations $f_1(t,x_1, \ldots, x_n)=f_2(t,x_1, \ldots, x_n) = \cdots  = f_n(t, x_1, \ldots, x_n)=0$ in $n$ variables $(x_1, \ldots, x_n)$ and one additional parameter $t$. Suppose that these equations have finitely many solutions for some fixed $t_0$, and that we have computed good numerical approximations to them. If we then move $t$ on a specified path from $t_0$ to some other $t_{\infty}$, how do the solutions to the equations move? The idea is very simple: Discretize your path into a number of small steps $t_0$ to $t_1$ to $t_2$ etcetera and use the known solutions at $t=t_i$ as starting points for some sort of Newton's method like solution at $t=t_{i+1}$. In detail, this is a hard numerical analysis problem and you probably should not try to implement it yourself.
This is your problem. Your equations are $f(x)=t$, where $x$ ranges over the Riemann surface $X$. Solve this equation for $t$ near $i$ (say), finding $N$ roots, and compute the monodromy as $t$ travels in a loop around $0$ and in a loop around $1$, each time staying far from $0$ and $1$. This gives you two permutations in $S_N$. The map from the free group to $S_N$ sends the generators of the free group to these permutations.
In short, I am saying that I don't think there is a method better than using the definition of the problem, but there is prepackaged software to implement that definition.
