The distribution for $\theta$ conditional on $s$ (and $n$) is
\begin{equation}
p(\theta|s) = \textsf{Beta}(\theta|s+1,n-s+1) . 
\end{equation}
The problem is that we don't observe $s$. 

Instead we are given $\pi = (\pi_1, \ldots, \pi_n)$, where 
\begin{equation}
p(x_i|\pi_i) = \textsf{Bernoulli}(x_i|\pi_i)  
\end{equation}
and 
\begin{equation}
p(x|\pi) = \prod_{i=1}^n p(x_i|\pi_i) .
\end{equation}
We can compute $p(s|\pi)$ from $p(x|\pi)$ and then
\begin{equation}
p(\theta|\pi) = \sum_{s=0}^n p(\theta|s)\,p(s|\pi) ,
\end{equation}
which is a mixture of beta distributions. 

The problem now is that $p(s|\pi)$ is the Poisson-binomial distribution which is difficult to compute analytically unless $n$ is quite small. One solution is to use simulation to approximate the distribution. In particular, let
\begin{equation}
p(\theta|\pi) \approx \frac{1}{R} \sum_{r=1}^R p(\theta|s^{(r)}) ,
\end{equation}
where $s^{(r)} \sim p(s|\pi)$. 

Let me provide more detail. Let
\begin{equation}
s^{(r)} = \sum_{i=1}^n x_i^{(r)} ,
\end{equation}
where
\begin{equation}
x_i^{(r)} \sim p(x_i|\pi_i) .
\end{equation}
Now define the approximate mixture weights:
\begin{equation}
\hat w_j = \frac{1}{R}\sum_{r=1}^R 1(s^{(r)} = j) ,
\end{equation}
where the indicator function is
\begin{equation}
1(x) = \begin{cases}
1 & \text{$x$ is true} \\
0 & \text{$x$ is false}
\end{cases} .
\end{equation}
Then 
\begin{equation}
p(\theta|\pi) \approx \sum_{j=0}^n \hat w_j\,p(\theta|s=j) . 
\end{equation}
You will need to experiment to see how large $R$ needs to be so that you get a sufficiently accurate answer. Depending on $n$ and $\pi$, you might need $R = 10^5$ or more.