The infinite product representation for $\vartheta_{1}(z, q)$ given by $$\vartheta_{1}(z, q) = 2q^{1/4}\sin z\prod_{n = 1}^{\infty}(1 - q^{2n})(1 - 2q^{2n}\cos 2z + q^{4n})\tag{1}$$ I hope you are familiar with the above result. Taking logs we get
\begin{align}\log\vartheta_{1}(z, q) &= \log(2q^{1/4}) + \log\sin z + \log\prod_{n = 1}^{\infty}(1 - q^{2n})\notag\\
&\,\,\,\,\,\,\,\, + \sum_{n = 1}^{\infty}\log(1 - 2q^{2n}\cos 2z + q^{4n})\tag{2}
\end{align}
and differentiating with respect to $z$ we get
\begin{align}
\frac{\vartheta'_{1}(z, q)}{\vartheta_{1}(z, q)} &= \cot z + 4\sum_{n = 1}^{\infty}\frac{q^{2n}\sin 2z}{1 - 2q^{2n}\cos 2z + q^{4n}}\notag\\
&= \cot z + \frac{4}{2i}\sum_{n = 1}^{\infty}\frac{q^{2n}(e^{2iz} - e^{-2iz})}{(1 - q^{2n}e^{2iz})(1 - q^{2n}e^{-2iz})}\notag\\
&= \cot z + \frac{4}{2i}\sum_{n = 1}^{\infty}q^{2n}\left(\frac{e^{2iz}}{1 - q^{2n}e^{2iz}} - \frac{e^{-2iz}}{1 - q^{2n}e^{-2iz}}\right)\notag\\
&= \cot z + \frac{4}{2i}\left(\sum_{n = 1}^{\infty}\frac{q^{2n}e^{2iz}}{1 - q^{2n}e^{2iz}} - \sum_{n = 1}^{\infty}\frac{q^{2n}e^{-2iz}}{1 - q^{2n}e^{-2iz}}\right)\notag\\
&= \cot z + \frac{4}{2i}\left(\sum_{n = 1}^{\infty}\sum_{m = 1}^{\infty}\{q^{2n}e^{2iz}\}^{m} - \sum_{n = 1}^{\infty}\sum_{m = 1}^{\infty}\{q^{2n}e^{-2iz}\}^{m}\right)\notag\\
&= \cot z + \frac{4}{2i}\left(\sum_{n = 1}^{\infty}\sum_{m = 1}^{\infty}q^{2m}e^{2imz}q^{2m(n - 1)}\right. \notag\\
&\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\, \left. - \sum_{n = 1}^{\infty}\sum_{m = 1}^{\infty}q^{2m}e^{-2imz}q^{2m(n - 1)}\right)\notag\\
&= \cot z + \frac{4}{2i}\left(\sum_{m = 1}^{\infty}q^{2m}e^{2imz}\sum_{n = 1}^{\infty}q^{2m(n - 1)}\right.\notag\\
&\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\, \left. - \sum_{m = 1}^{\infty}q^{2m}e^{-2imz}\sum_{n = 1}^{\infty}q^{2m(n - 1)}\right)\notag\\
&= \cot z + \frac{4}{2i}\left(\sum_{m = 1}^{\infty}\frac{q^{2m}e^{2imz}}{1 - q^{2m}} - \sum_{m = 1}^{\infty}\frac{q^{2m}e^{-2imz}}{1 - q^{2m}}\right)\notag\\
&= \cot z + 4\sum_{n = 1}^{\infty}\frac{q^{2n}}{1 - q^{2n}}\sin 2nz\tag{3}\end{align}
The above proof is taken from one of my blog posts (see equation $(25)$ of that post).
I do not know a closed form for
$$I_n=\int_{-1}^{+1} (x;x)_{n+1}\,dx$$
For even values of $n$, this generates the sequence
$$\left\{2,\frac{152}{105},\frac{4672}{3465},\frac{333280768}{253161909},\frac{141841
28217088}{10883356008525},\frac{257964528813440303104}{198899418641149766925}\right\}$$ For odd values of $n$
$$\left\{\frac{4}{3},\frac{128}{99},\frac{328448}{255255},\frac{34963456}{27188525},
\frac{530068713377923072}{412106981947605945}\right\}$$ which both decrese when $n$ increase. As you wrote, I do not see any clear pattern.
I did not find these terms in $OEIS$.
Concerning
$$I_\infty=\int_{-1}^{+1} (x;x)_\infty\,dx$$ it is just a number $\color{red}{\text{(this is WRONG})}$
$$I_\infty=1.28830088867392123018090\cdots$$ which is not recognized by inverse symbolic calculators.
Edit (after @Brian Constantinescu's comment)
In Quora, using residues, @Brian Constantinescu provided a very elegant closed form of the result.
Simplifying his last result
$$I_\infty=\int_{-1}^{+1} (x;x)_\infty\,dx=4\pi \sqrt{\frac{6}{23}}\,\, \frac{\sinh \left(\frac{\sqrt{23}
\pi }{4}\right)+\sqrt{2} \sinh \left(\frac{\sqrt{23} \pi }{6}\right)}{2 \cosh \left(\frac{\sqrt{23} \pi }{3}\right)-1}$$
Using the same steps and notations as @Brian Constantinescu, computing the sums in terms of polylogarithms and simplifying, we have
$$A=\sum_{-\infty}^{\infty}\frac{(-1)^n}{3 n^2-n+2}$$ $$A=-\frac{i \pi \left(\tan \left(\frac{1}{12} \left(5+i \sqrt{23}\right) \pi \right)+\cot
\left(\frac{1}{12} \left(5+i \sqrt{23}\right) \pi \right)-2 \csc \left(\frac{1}{6}
\left(1+i \sqrt{23}\right) \pi \right)\right)}{2 \sqrt{23}}$$ Expanding the complex numbers
$$\color{blue}{A=4 \pi \sqrt{\frac{3}{23}}\,\,\frac{\sinh \left(\frac{\sqrt{23} \pi }{6}\right)}{2 \cosh \left(\frac{\sqrt{23} \pi}{3}\right)-1}} $$
$$B_x=\sum_{-\infty}^{\infty}\frac{(-1)^x}{12 x^2-2 x+2}$$
$$B_x=\frac{i \pi \left(\csc \left(\frac{1}{12} \left(1+i \sqrt{23}\right) \pi
\right)-\text{sech}\left(\frac{1}{12} \left(\sqrt{23}-5 i\right) \pi \right)\right)}{2\sqrt{23}}$$
$$\color{blue}{B_x=-\pi\sqrt{\frac{2}{23}} \left(\sqrt{3}+1\right)\,\,\frac{\sinh \left(\frac{\sqrt{23} \pi }{12}\right)}{\sqrt{3}-2 \cosh
\left(\frac{\sqrt{23} \pi }{6}\right)}}$$
$$B_y=\sum_{-\infty}^{\infty}\frac{(-1)^y}{12 y^2+10 y+4}$$
$$B_y=\frac{i \pi \left(\csc \left(\frac{1}{12} \left(5+i \sqrt{23}\right) \pi
\right)-\text{sech}\left(\frac{1}{12} \left(\sqrt{23}-i\right) \pi \right)\right)}{2 \sqrt{23}}$$
$$\color{blue}{B_y=\pi\sqrt{\frac{2}{23}} \left(\sqrt{3}-1\right)\,\,\frac{\sinh \left(\frac{\sqrt{23} \pi }{12}\right)}{\sqrt{3}+2 \cosh
\left(\frac{\sqrt{23} \pi }{6}\right)}}$$
$$\color{red}{I_\infty=2(A+B_x+B_y)=4\pi \sqrt{\frac{6}{23}}\,\, \frac{\sinh \left(\frac{\sqrt{23}
\pi }{4}\right)+\sqrt{2} \sinh \left(\frac{\sqrt{23} \pi }{6}\right)}{2 \cosh \left(\frac{\sqrt{23} \pi }{3}\right)-1}}$$
Best Answer
A possible approach could be to search an approximated estimate of the summation. Applying the Euler-Maclaurin formula to $f(n)=a^nq^{n^2}$ in the interval between $0$ and $N$, we have
$$\displaystyle \sum_0^N a^nq^{n^2}=\int_0^N a^nq^{n^2}+ \frac{a^Nq^{N^2}+1}{2} \\ + \sum_{k=1}^{\infty }\,{\frac {B_{2k}}{(2k)!}}\left(f^{(2k-1)}(N)-f^{(2k-1)}(0)\right)$$
where $B_{2k}$ denotes the $2k^{th}$ Bernoulli number and $f^{(2k-1)}$ is the $(2k-1)^{th}$ derivative.
To estimate the integral, we can make the substitution $$n=\frac{t}{\sqrt{\log q}}-\frac{\log a}{2 \log q}$$ We have
$$\displaystyle a^nq^{n^2}=a^{t/\sqrt{\log q} -\log a/(2\log q) } q^ {[t/\sqrt{\log q} -\log a/(2\log q) ]^2} \\ = e^{-\log^2 a/(4 \log q)}e^{t^2} $$
Also, we have
$$t=\frac{\log a +2n \log q}{2\sqrt{\log q}}$$ $$dt=\sqrt{\log q } \, dn$$ $$dn=\frac{1}{\sqrt{\log q}} dt$$
Then our integral can be rewritten as
$$\displaystyle \int_0^N a^nq^{n^2}dn\\= \frac{e^{-\log^2 a/(4 \log q)}}{\sqrt{\log q}} \int_{\frac{\log a}{2\sqrt{\log q}}}^{\frac{\log a +2N\log q}{2\sqrt{\log q}}} e^{t^2}$$
Setting by simplicity $$K_1=\log a/(2\sqrt{\log q})$$ and $$K_2=(\log a +2N\log q)/(2\sqrt{\log q})$$ we can further rewrite the integral as
$$\displaystyle \int_0^N a^nq^{n^2}dn =\frac{e^{-K_1^2 }}{\sqrt{\log q}} \int_{K_1}^{K_2} q^{t^2}\\ =\frac{e^{-K_1^2}}{\sqrt{\log q}} \left(\int_{0}^{K_2} q^{t^2} -\int_{0}^{K_1} q^{t^2} \right)\\ = \frac{\sqrt{\pi}\,e^{-K_1^2 }}{2\sqrt{\log q}} \bigg[\operatorname{erfi}(K_2) - \operatorname{erfi}(K_1)\bigg]\\$$
where $\operatorname{erfi}$ indicates the imaginary error function, defined as
$$\displaystyle \operatorname {erfi} (x)={\frac {2}{\sqrt {\pi }}}\int_{0}^{x}e^{t^{2}}\,dt$$
Alternatively, our integral could be written in terms of the Dawson function. In fact
$$\displaystyle \int_0^N a^nq^{n^2}dn\\ =\frac{e^{-\frac{\log^2 a}{4\log q} }}{\sqrt{\log q}} \left(\int_{0}^{K_2} q^{t^2} -\int_{0}^{K_1} q^{t^2} \right)\\ = \frac{e^{-K_2^2 +n \log a +2n^2\log q }}{\sqrt{\log q}} \int_{0}^{K_2} q^{t^2} \\ - \frac{e^{-K_1^2 }}{\sqrt{\log q}} \int_{0}^{K_1} q^{t^2} \\ =\frac{a^Nq^{N^2}}{\sqrt{\log q}} \, e^{-K_2^2} \int_{0}^{K_2} e^{t^2}\\ - \frac{e^{-K_1^1}}{ \sqrt{\log q}} \int_{0}^{K_1} e^{t^2} \\ =\frac{{a^Nq^{N^2} \,F_+(K_2) }- {F_+(K_1)}}{\sqrt{\log q}} $$
where $F_+(K_2)$ and $F_+(K_1)$ indicate the Dawson function, i.e. the one-sided Fourier–Laplace sine transform of the Gaussian function, usually expressed in the following two forms:
$${\displaystyle F_{+}(x)=e^{-x^{2}}\int _{0}^{x}e^{t^{2}}\,dt}$$
$${\displaystyle F_{-}(x)=e^{x^{2}}\int _{0}^{x}e^{-t^{2}}\,dt.\!}$$
Note that the Dawson function is linked to the error function and the imaginary error function by
$${\displaystyle \operatorname {erf} (x)={\frac {2}{\sqrt {\pi }}e^{-x^{2}}F_-(x)}}$$
$${\displaystyle \operatorname {erfi} (x)={\frac {2}{\sqrt {\pi }}e^{-x^{2}}F_+(x)}}$$
To proceed further, we have to note that the initial sum $\sum_0^\infty a^n q^{n^2}$ converges for $q<1$ and diverges for $q>1$ (the case $q=1$ reduces to the trivial sum $\sum_0^\infty a^n$, whose convergence or divergence depends on $a$). These two cases determine the presence of real or imaginary terms in the last equations. Since the OP asks for the result of infinite sums, I will focus only on the first case. However, the same following approach could be used to estimate partial sums in the second case.
If $q<1$, then $K_2$, $K_1$, and the term $\sqrt{ \log q}$ in the denominators of the RHS of the equations above are all imaginary. So if we write
$$K_2=\frac{\log a +2N\log q}{2\sqrt{\log q}}= \frac{\log a +2N\log q}{2\sqrt{-\log 1/q}}= -i\frac{\log a +2N\log q}{2\sqrt{\log 1/q}} =i\hat{K_2} $$
and similarly
$$K_1=\frac{\log a }{2\sqrt{\log q}}= \frac{\log a }{2\sqrt{-\log 1/q}}= -i\frac{\log a }{2\sqrt{\log 1/q}} =i\hat{K_1} $$
reminding that $\operatorname{erfi}(z i)= i\,\operatorname{erf}(z)$, we have
$$\displaystyle \int_0^N a^nq^{n^2}dn\\ = \frac{\sqrt{\pi} \, e^{\hat{K_1}^2} }{2\sqrt{\log q}}\left[ { \operatorname{erfi}(i\hat{K_2}) }- { \operatorname{erfi}(i\hat{K_1}) }\right] \\ = \frac{\sqrt{\pi} \, e^{\hat{K_1}^2} } {2\sqrt{\log (1/q)}}\left[ { \operatorname{erf}(\hat{K_2}) }- { \operatorname{erf}(\hat{K_1}) }\right] $$
Now we can consider that, as $N$ increases, the error function shows a fast convergence to $1$, satisfying $\operatorname{erf}(z)=1+O(e^{-z^2}z^{-1})$. Therefore, making the inverse substitution $\hat{K_1}= - (\log a )/(2\sqrt{\log 1/q}) $, we get
$$\displaystyle \int_0^\infty a^nq^{n^2}dn\\ = \frac{\sqrt{\pi} \, e^{ \frac{\log^2 a}{4\log 1/q}} } {2\sqrt{\log (1/q)}}\left[ 1- { \operatorname{erf}\left(-\frac{ \log a }{2\sqrt{\log 1/q}}\right) }\right] $$
Turning back to the initial Euler-Maclaurin expansion, the degree of approximation for the sum $\sum_0^\infty a^nq^{n^2}$ depends on the number of terms considered after the integral. Among these, it is not difficult to note that the leading terms result from the values of the function and the derivatives in $n=0$. For most values of $a$ and $q$, a good level of approximation can be obtained considering these terms up to the third derivative. Because $$f^{(1)}(0)=\log a$$ $$f^{(3)}(0)=\log^2 a + 6\log a \log q$$
and $$B_2=\frac{1}{6}$$ $$B_4=-\frac{1}{30}$$
we have
$$\displaystyle \sum_0^N a^nq^{n^2}\approx \frac{\sqrt{\pi} \, e^{ \frac{\log^2 a}{4\log 1/q}} } {2\sqrt{\log (1/q)}}\left[ 1- { \operatorname{erf}\left(-\frac{ \log a }{2\sqrt{\log 1/q}}\right) }\right] + \frac{1}{2} -\frac{\log a}{12}+\frac{\log a(\log^2 a +6 \log q)}{720} $$
For example, setting $a=1/2$ and $q=3/4$, we have
$$\sum_0^\infty a^nq^{n^2}=1.46413758...$$
as shown by WA here. Our approximated formula yields
$$\displaystyle \sum_0^N a^nq^{n^2}\approx \frac{\sqrt{\pi} \, e^{ \frac{\log^2 1/2}{4\log 4/3}} } {2\sqrt{\log (4/3)}}\left[ 1- { \operatorname{erf}\left(-\frac{ \log 1/2 }{2\sqrt{\log 4/3}}\right) }\right] + \frac{1}{2} -\frac{\log 1/2}{12}+\frac{\log 1/2 (\log^2 1/2 +6 \log 4/3)}{720} \\ =1.46407819...$$
as confirmed by WA here. The relative error, expressed in percentage, is $-0.004\%$.
Including further terms of the Euler-Maclaurin espansion, the approximation can be improved. For example, considering the terms up to the fifth derivative, taking into account that
$$f^{(5)}(0)\\=\log^5 a + 20 \log^3 a \log q \\ + 60 \log a \log^2 q $$
and $$B_6=-\frac{1}{42}$$
the approximated formula above would include the additional term
$$-\frac{1}{30240} \bigg[\log^5 (1/2) + 20 \log^3 (1/2) \log (3/4) + 60 \log(1/2) \log^2 (3/4) \bigg]=0.00005574... $$
as shown here. In this way, the estimated sum would become
$$=1.46407819 + 0.00005574 = 1.46413393 $$
with a relative error of $-0.00025\%$.
Lastly, to facilitate the computation of the approximated formula (e.g. in Python), the value of the $\operatorname{erf}$ function can be calculated using several numerical approximations. These allow to calculate the value of $$\operatorname{erf}\left(-\frac{ \log a }{2\sqrt{\log 1/q}}\right) $$ with good precision using elementary functions.