I am getting into general relativity, which of course means getting to grips with curvilinear coordinate systems. Repeatedly, the textbook and lectures have emphasized the relationships $$e\_i=g\_{ij}e^j\quad \text{and}\quad g\_{ij}^{-1}=g^{ij}$$

In a problem set, I have been told that, given $$s(r, \theta, \phi) = r\begin{pmatrix}
\sin \theta \cos \phi\\
\sin \theta \sin \phi\\
\cos \theta
\end{pmatrix}$$
to find $e\_r,e\_\theta,e\_\phi$ and then their dual, $e^r,e^\theta,e^\phi$.

I have correctly found $e\_r,e\_\theta,e\_\phi$ using the definition in the text, $e\_j=\frac{\partial s}{\partial q\_j}$, and I have found the correct $g\_{ij}=e\_i\cdot e\_j$, but when I use the above relation, I don't get the correct $e^r,e^\theta,e^\phi$, based on the orthonormality (I get neither orthogonal nor normality) requirement.

Why? As far as I can tell, I should be able to construct the dual basis from the original in this way, but somewhere in the lectures and text I have missed something.

EDIT: to be more explicit on what I have done.

$$e\_r=\frac{\partial s}{\partial r}=\begin{pmatrix}
\sin \theta \cos \phi\\
\sin \theta \sin \phi\\
\cos \theta
\end{pmatrix}$$
$$e\_\theta=\frac{\partial s}{\partial \theta}=r \begin{pmatrix}
\cos\theta\cos\phi\\
\cos\theta\sin\phi\\
-\sin\theta
\end{pmatrix}$$
$$e\_\phi=\frac{\partial s}{\partial\phi}=r \begin{pmatrix}
-\sin\theta\sin\phi\\
\sin \theta \cos \phi\\
0
\end{pmatrix}$$
We note that these vectors are all orthogonal, and so the only nonzero components of the metric are the diagonal ones. $e\_r\cdot e\_r=1$, $e\_\theta\cdot e\_\theta=r^2$, and $e\_\phi\cdot e\_\phi=r^2\sin(\theta)$. This gives
$$g\_{ij}=\begin{pmatrix}
 1 & 0 & 0 \\
 0 & r^2 & 0 \\
 0 & 0 & r^2\sin^2\theta \\
 \end{pmatrix}$$
and we also easily have
$$g^{ij}=\begin{pmatrix}
 1 & 0 & 0 \\
 0 & 1/r^2 & 0 \\
 0 & 0 & 1/(r^2\sin^2\theta) \\
 \end{pmatrix}$$

Finally, based on the first relation, we have that
$$e^r=g^{ij}e\_r=\begin{pmatrix}
\sin \theta \cos \phi\\
\sin \theta \sin \phi/r^2\\
\cos \theta/r^2\sin^2\theta
\end{pmatrix}$$
which already shows the problem. This is not orthonormal with the any of the original basis vectors. I have that the actual answer (constructed another way) is $e^r=e\_r$, $e^\theta=e\_\theta/r^2$, and $e^\phi=e\_\phi/r^2\sin^2\theta$.