text: import numpy as np
  import numpy.random as npr
  import scipy.stats as ss
  import utilities as ut
  
  
  class Model:
      """
      The models (objects of this class) represent quantum circuits that
      operate on na + nb qubits.
  
      We will call a list1 any list of numpy arrays of shapes given by the
      shapes1, with length given by len1 = len(shapes1). len1 will be referred
      to as the number of layers.
  
      The first version of Quantum Edward considers two models called NbTrols
      and NoNbTrols, but QEdward is more general than those. QEdward was
      written so that it can also analyze other akin models. Akin models
      should have len1 'layers'.
  
      len1=na+nb for both models NbTrols and NoNbTrols, but they have
      different shapes1. In fact, for the NoNbTrols model
  
      self.shapes1 = [(1 << k,) for k in range(na)] + [(powna,)]*nb
  
      where powna = 2^na = 1 << na, whereas for the NbTrols model,
  
      self.shapes1 = [(1 << k,) for k in range(na+nb)]
  
      The Wikipedia article cited below refers to the Beta distribution
      concentrations conc0 and conc1 as alpha and beta, respectively.
              The Bayesian net for the model looks like this
  
               list1_angs[0:na]
             /
          x    list1_angs[na:]
          |  /
          y
  
          with the arrows pointing downward.
  
          y is an int in range(pownb). x is an int in range(powna). list1_angs
          is a list1 of numpy arrays containing qc qubit rotation angles.
  
          y and x are observed variables. nsam measurements of (y, x) are
          given by the training data. list1_angs are hidden variables. We want
          to estimate those values for list1_angs which best fit the training
          data.
  
  
      References
      ----------
      https://en.wikipedia.org/wiki/Beta_distribution
  
      """
  
      def __init__(self, nb, na, list1_conc0_prior=None,
                       list1_conc1_prior=None):
          """
          Constructor
  
          Parameters
          ----------
          nb : int
              Total number of qubits is na + nb. Input goes into first na qubits
              and only last nb measured
          na : int
          list1_conc0_prior : list[np.array]
              a list1 of conc0 for a beta distribution. Used as priors
          list1_conc1_prior : list[np.array]
              a list1 of conc1 for a beta distribution. Used as priors
  
          Returns
          -------
          None
  
          """
          self.na = na
          self.nb = nb
          
          self.shapes1 = self.get_shapes1()
          self.len1 = len(self.shapes1)
          # print("..", self.len1)
  
          self.list1_conc0_prior = list1_conc0_prior
          self.list1_conc1_prior = list1_conc1_prior
          if list1_conc0_prior is None or list1_conc1_prior is None:
              self.list1_conc0_prior = ut.new_uniformly_random_array_list(
                  low=0., high=5., shapes=self.shapes1)
              self.list1_conc1_prior = ut.new_uniformly_random_array_list(
                  low=0., high=5., shapes=self.shapes1)
  
          # with beta function, if two concentrations are equal,
          # mean_x beta(x, conc0=1, conc1=1) = .5
          # ang = prob*dpi
          self.list1_angs_prior = [ss.beta.mean(self.list1_conc0_prior[k],
                                   self.list1_conc1_prior[k])*ut.dpi
                                   for k in range(self.len1)]
          # print(';;', self.list1_angs_prior)
  
      def get_shapes1(self):
          """
          Abstract method. Must be overridden by descendant class. In
          descendant class, this function should return a list of the shapes
          of the elements of a list1.
  
          Returns
          -------
          list[tuple]
  
          """
          assert False, "this function must be overridden"
  
      def prob_x(self, x,
                 list1_angs,
                 verbose=False):
          """
          Abstract method. Must be overridden by descendant class. Returns
          probability of input x, P(x).
  
          Parameters
          ----------
          x : int
              x is an int in range(powna).
          list1_angs : list[np.array]
          verbose : bool
  
          Returns
          -------
          float
  
          """
          assert False, "this function must be overridden"
  
      def prob_y_given_x_and_angs_prior(self, y, x,
                                        list1_angs,
                                        verbose=False):
          """
          Abstract method. Must be overridden by descendant class. In
          descendant class, should return the probability of y given x and
          list1_angs, P(y | x, list1_angs).
  
  
          Parameters
          ----------
          y : int
              y is an int in range(pownb).
          x : int
              x is an int in range(powna).
          list1_angs : list[np.array]
          verbose : bool
  
          Returns
          -------
          float
  
          """
          assert False, "this function must be overridden"
  
      def sum_over_x_of_prob_x(self, list1_angs):
          """
          This function should return 1. It can be used to check that the
          function prob_x() yields a set of probabilities which when summed
          over x gives 1.
  
          Parameters
          ----------
          list1_angs : list[np.array]
  
          Returns
          -------
          float
  
          """
          powna = 1 << self.na
          tot_prob = 0.
          for x in range(powna):
              tot_prob == self.prob_x(x, list1_angs)
          return tot_prob
  
      def sum_over_y_of_prob_y(self, x, list1_angs):
          """
          This function should return 1. It can be used to check that the
          function prob_y_given_x_and_angs_prior() yields a set of
          probabilities which when summed over y gives 1.
  
          Parameters
          ----------
          x : int
          list1_angs : list[np.array]
  
          Returns
          -------
          float
  
          """
          pownb = 1 << self.nb
          tot_prob = 0.
          for y in range(pownb):
              tot_prob == self.prob_y_given_x_and_angs_prior(y, x,
                                                             list1_angs)
          return tot_prob
  
      def gen_toy_data(self, nsam, verbose=False):
          """
          This function generates nsam samples of (y, x). It returns these
          samples as arrays with 0 and 1 entries y_nsam_nb and x_nsam_nb.
          y_nsam_nb.shape = (nsam, nb), x_nsam_na.shape = (nsam, na).
  
          For example, x_nsam_na[0, :] could look like [1, 0, 1, 1], meaning
          that with na=4 and nb=3, the qubit values are
  
          q0=1
          q1=0
          q2=1
          q3=1
          q4=0
          q5=0
          q6=0
  
          Parameters
          ----------
          nsam : int
          verbose : bool
  
          Returns
          -------
          y_nsam_nb, x_nsam_na : tuple(np.array, np.array)
  
          """
  
          powna = 1 << self.na
          # x_nsam entries are integers in range(powna)
          # x_nsam_na entries are integers 0 or 1
          # multinomial with just one draw will give vec with
          # all entries 0 except one entry equal to 1
          pvals = [self.prob_x(x, self.list1_angs_prior,
                              verbose) for x in range(powna)]
          if verbose:
              print('tot_prob=', np.sum(pvals), '\n')
          x_nsam = np.array(
              [list(npr.multinomial(n=1, pvals=pvals)).index(1)
               for sam in range(nsam)])
          x_nsam_na = ut.dec_to_bin_vec(x_nsam, self.na, nsam=nsam)
  
          # p_nsam_pownb entries will be probabilities
          pownb = 1 << self.nb
          p_nsam_pownb = np.zeros((nsam, pownb))
          for sam in range(nsam):
              p_nsam_pownb[sam, :] = np.array([
                  self.prob_y_given_x_and_angs_prior(y, x_nsam[sam],
                                                     self.list1_angs_prior,
                                                     verbose)
                  for y in range(pownb)])
              if verbose:
                  print("sample=", sam,
                        "  total_prob=", np.sum(p_nsam_pownb[sam, :]),
                        '\n')
  
          # multinomial with just one draw will give vec with
          # all entries 0 except one entry equal to 1
          y_nsam = np.array(
              [list(npr.multinomial(n=1, pvals=p_nsam_pownb[sam, :])).index(1)
               for sam in range(nsam)])
  
          y_nsam_nb = ut.dec_to_bin_vec(y_nsam, self.nb, nsam=nsam)
  
          return y_nsam_nb, x_nsam_na
  
  
  if __name__ == "__main__":
      def main():
          print(5)
      main()
  