def g():
        d = len(Mu)
        assert Mu.shape == (d,), "Mu must be a vector"
        assert A.shape == (d, d), "A must be a square matrix"
        assert (A.T == A).all(), "and symmetric"
        assert V.shape == (d, d), "V must be a square matrix"
        assert (V.T == V).all(), "and symmetric"

        a = chol(A)
        v = chol(V)

        B = dot(V, inv(V + A))
        _a2 = V - dot(B, V)
        _a2 = chol(_a2)

        Y, U = array([0.0] * d), array([0.0] * d)

        for i in range(n + burnin):
            for _ in range(thin):  # skipe
                # sample Y | U ~ N(U, V)
                Y = U + dot(v, random.normal(size=d))

                # sample U | Y ~ N(A(A+V)^-1)(Y-Mu) + Mu,
                #                  A - A(A+V)^-1A)
                U = dot(B, (Mu - Y)) + Y + +dot(_a2, random.normal(size=d))

            if i >= burnin:
                yield [U, Y]

 def testRewardFunction(self, x, typ, noise=0.000001):
     if typ == "growSin":
         return (sin((x - self.rangeMin) / 3.0) + 1.5 + x / self.distRange) / 4.0 + random.normal(0, noise)
     if typ == "rastrigin":
         n = x / self.distRange * 10.0
         if abs(n) > 5.0: n = 0.0
         # FIXME: imprecise reimplementation of the Rastrigin function that exists already
         # in rl/environments/functions...
         return (20.0 + n ** 2 - 10.0 * cos(2.0 * 3.1416 * n)) / 55.0 + random.normal(0, noise)
     if typ == "singleGaus":
         return self.getStND(x) + random.normal(0, noise)
     return 0.0

def quadthermo_agent_gen():
	targetT=spr.normal(18,2)
	tolerance=2

	absmax=max(spr.normal(21,1),targetT+tolerance) #no support
	absmin=min(spr.gamma(4,1),targetT-tolerance) #no support

	vec=[0.0,0.21,0.135,0.205,0.115,0.115,0.095,0.065,0.03,0.03]
	r=sp.random.uniform()
	nres=1
	i=0
	
	while i<7 and r>vec[i]:
		r-=vec[i]
		i+=1
	
		
	nres=i+1
	
	occ=active.get_occ_p(nres)
	cons=[]
	for p in occ:
		p.extend([targetT,tolerance])
		cons.append(p)
	q=sp.exp(sp.random.normal(-9.3,2.0))
	
	
	fa=28.39+sp.random.gamma(shape=2.099,scale=28.696) #floor area from cabe dwelling survey
	
	flat=sp.random.uniform(0,1)<0.365 # flat or house
	if flat:
		U = 3.3*sp.sqrt(fa) #insulation in W/K floor area
	else:
		U= 3.6*sp.sqrt(fa)+0.14*fa
	
	
	k=U #insulation in W/K

	cm=1000*sp.exp(sp.random.normal(5.5,0.35)) #thermal capacity in J
	
	P=sp.random.uniform(6000,15000)# power in W

	Prequ=k*20
	
	if Prequ>P:
		print "!!!!!!!!!!!"
	s=str(["quadthermo_agent",absmax,absmin,cons,q,P,cm,k])
	return s

    def newEpisode(self):
        if self.learning:
            params = ravel(self.explorationlayer.module.params)
            target = ravel(sum(self.history.getSequence(self.history.getNumSequences()-1)[2]) / 500)
        
            if target != 0.0:
                self.gp.addSample(params, target)
                if len(self.gp.trainx) > 20:
                    self.gp.trainx = self.gp.trainx[-20:, :]
                    self.gp.trainy = self.gp.trainy[-20:]
                    self.gp.noise = self.gp.noise[-20:]
                    
                self.gp._calculate()
                        
                # get new parameters where mean was highest
                max_cov = diag(self.gp.pred_cov).max()
                indices = where(diag(self.gp.pred_cov) == max_cov)[0]
                pick = indices[random.randint(len(indices))]
                new_param = self.gp.testx[pick]
            
                # check if that one exists already in gp training set
                if len(where(self.gp.trainx == new_param)[0]) > 0:
                    # add some normal noise to it
                    new_param += random.normal(0, 1, len(new_param))

                self.explorationlayer.module._setParameters(new_param)

            else:
                self.explorationlayer.drawRandomWeights()
        
        # don't call StateDependentAgent.newEpisode() because it randomizes the params
        LearningAgent.newEpisode(self)

def main():
	

	ITERATIONS = 100
	mc = zeros(ITERATIONS)
	og = zeros(ITERATIONS)

	#farby = 
	QS = 10

	colors = [
		[0, 0, 0],
		[0, 0, 1],
		[0, 1, 0],
		[1, 0, 0],
		[0, 1, 1],
		[1, 0, 1],
		[0, 1, 1],
	]

	for qpre in range(QS):
		q =  qpre + 2
		for it in range(ITERATIONS):
			W = random.normal(0, 0.1, [q, q])
			WI = random.uniform(-.1, .1, [q, 1])
			mc[it] = sum(memory_capacity(W, WI, memory_max=200, runs=1, iterations_coef_measure=5000)[0][:q+2])
			og[it] = matrix_orthogonality(W)
			print(qpre, QS, it, ITERATIONS)
		plt.scatter(og, mc, marker='+', label=q, c=(colors[qpre % len(colors)]))

	plt.xlabel("orthogonality")
	plt.ylabel("memory capacity")
	plt.grid(True)
	plt.legend()
	plt.show()

 def drop_object(self):
     """Drops a random object (box, sphere) into the scene."""
     # choose between boxes and spheres
     if random.uniform() > 0.5:
         (body, geom) = self._create_sphere(self.space, 10, 0.4)
     else:
         (body, geom) = self._create_box(self.space, 10, 0.5, 0.5, 0.5)
     # randomize position slightly
     body.setPosition((random.normal(-6.5, 0.5), 6.0, random.normal(-6.5, 0.5)))
     # body.setPosition( (0.0, 3.0, 0.0) )
     # randomize orientation slightly
     #theta = random.uniform(0,2*pi)
     #ct = cos (theta)
     #st = sin (theta)
     # rotate body and append to (body,geom) tuple list
     # body.setRotation([ct, 0., -st, 0., 1., 0., st, 0., ct])
     self.body_geom.append((body, geom))

def GaussianRandomInitializer(gridShape, sigma=0.2, seed=None, slipSystem=None, slipPlanes=None, slipDirections=None, vacancy=None, smectic=None):

    oldgrid = copy.copy(gridShape)
   
    if len(gridShape) == 1:
	    gridShape = (128,)
    if len(gridShape) == 2:
	    gridShape = (128,128)
    if len(gridShape) == 3:
	    gridShape = (128,128,128)

    """ Returns a random initial set of fields of class type PlasticityState """
    if slipSystem=='gamma':
        state = SlipSystemState.SlipSystemState(gridShape,slipPlanes=slipPlanes,slipDirections=slipDirections)
    elif slipSystem=='betaP':
        state = SlipSystemBetaPState.SlipSystemState(gridShape,slipPlanes=slipPlanes,slipDirections=slipDirections)
    else:
        if vacancy is not None:
            state = VacancyState.VacancyState(gridShape,alpha=vacancy)
        elif smectic is not None:
            state = SmecticState.SmecticState(gridShape)
        else:
            state = PlasticityState.PlasticityState(gridShape)

    field = state.GetOrderParameterField()
    Ksq_prime = FourierSpaceTools.FourierSpaceTools(gridShape).kSq * (-sigma**2/4.)

    if seed is None:
        seed = 0
    n = 0
    random.seed(seed)

    Ksq = FourierSpaceTools.FourierSpaceTools(gridShape).kSq.numpy_array()

    for component in field.components:
        temp = random.normal(scale=gridShape[0],size=gridShape)
        ktemp = fft.rfftn(temp)*(sqrt(pi)*sigma)**len(gridShape)*exp(-Ksq*sigma**2/4.)
        field[component] = numpy.real(fft.irfftn(ktemp))
        #field[component] = GenerateGaussianRandomArray(gridShape, temp ,sigma)
        n += 1

    """
    t, s = LoadState("2dstate32.save", 0)
    for component in field.components:
        for j in range(0,32):
            field[component][:,:,j] = s.betaP[component].numpy_array()
    """

    ## To make seed consistent across grid sizes and convergence comparison
    gridShape = copy.copy(oldgrid)
    if gridShape[0] != 128:
        state = ResizeState(state,gridShape[0],Dim=len(gridShape))

    state = ReformatState(state)
    state.ktools = FourierSpaceTools.FourierSpaceTools(gridShape)
    
    return state 

def keplerSim(tau, e, T0, K, w, sig, tlo, thi, n):
	dt = (thi-tlo)/(n-1)
	data = zeros((n,2), Float)
	data[:,0] = r.uniform(tlo, thi, (n))
#	for i in range(n):
#            data[i,0] = tlo + i*dt
        data[:,1] = v_rad(K, w, tau, e, T0, data[:,0])+r.normal(0.,sig,(n))
        print "Created data."
	return data

    def perturbation(self):
        """ Generate a difference vector with the given standard deviations """
        #print self.sigList
        #print "_*_*_*_*_*_*_*_"
        #raw_input("Press Enter to continue")
        #time.sleep(3)

        
    
        return random.normal(0., self.sigList)

 def drawSample(self):
     sum = 0.0
     rndFakt = random.random()
     for g in range(self.numOGaus):
         sum += self.sigmo(self.alpha[g])
         if rndFakt < sum:
             if self.sigma[g] < self.minSig: self.sigma[g] = self.minSig
             x = random.normal(self.mue[g], self.sigma[g])
             break
     return x

def fill_region(l,r,sigma,v):
    if (l == r or l == r-1):
        pass
    else:
        m = int(round((r+l)*0.5))
        a = v[l] + (v[r]-v[l])*(m - l)/float(r - l)
        s = sigma*sqrt((m-l)*(r-m)/float(r-l))
        v[m] = a + s * random.normal()
        fill_region(l,m,sigma,v)
        fill_region(m,r,sigma,v)

 def add_noise(self, temp):
   """
   Add per-pixel Gaussian random noise.
   """
   self.noise = random.normal(0,temp,[self.npix,self.npix])
   self.Fnoise = ft.fftshift(ft.fft2(self.noise))
   self.Txy = self.Txy + self.noise
   self.Fxy = ft.fftshift(ft.fft2(self.Txy))
   
   self.Clnoise = ((temp*self.mapsize_rad/self.npix)*self.Bl)**-2.e0
   self.Pknoise = np.interp(self.modk, self.k, self.Clnoise)

def generate_data(N=100, true_params=secret_true_params,
                  seed = 42):
  x = np.linspace(-2.5, 2.5, N)
  y1 = my_model(x, *true_params)
  y2 = 1.0 * random.normal(size=N)
  # Create the data
  data = np.array([x,y1+y2]).T
  # Shuffle the data
  permuted_data = random.permutation(data)
  # Save the data
  np.savetxt("dataN%d.txt"%N, data)
  return data

    def __init__(self, statedim, actiondim, sigma= -2.):
        Explorer.__init__(self, actiondim, actiondim)
        self.statedim = statedim
        self.actiondim = actiondim

        # initialize parameters to sigma
        ParameterContainer.__init__(self, actiondim, stdParams=0)
        self.sigma = [sigma] * actiondim

        # exploration matrix (linear function)
        self.explmatrix = random.normal(0., expln(self.sigma), (statedim, actiondim))

        # store last state
        self.state = None

 def drawSample(self, dm):
     sum = 0.0
     rndFakt = random.random()
     if dm == "max":
         for g in range(self.numOGaus):
             sum += self.sigmo(self.alpha[g])
             if rndFakt < sum:
                 if self.sigma[g] < self.minSig: self.sigma[g] = self.minSig
                 x = random.normal(self.mue[g], self.sigma[g])
                 break
         return x
     if dm == "dist":
         return rndFakt * self.distRange + self.rangeMin
     return 0.0

        def model(times):
            t_fold, t_fold_model = self.period_folding(times, available, m, m_err, out_dict)
            data = empty(0)
            rms = empty(0)
            for time in t_fold_model:
                # we're going to create a window around the desired time and sample a gaussian distribution around that time
                period = 1.0 / f
                assert (
                    period < available.ptp() * 1.5
                ), (
                    "period is greater than ####SEE VARIABLE CONTSTRAINT#### of the duration of available data"
                )  # alterring this.  originally 1/3
                # window is 2% of the period
                passed = False
                for x in arange(0.01, 0.1, 0.01):
                    t_min = time - x * period
                    t_max = time + x * period
                    window = logical_and(
                        (t_fold < t_max), (t_fold > t_min)
                    )  # picks the available times that are within that window
                    try:
                        # there must be more than # points in the window for this to work:
                        assert window.sum() >= 2, str(time)  # jhiggins changed sum from 5 to 2
                    except AssertionError:
                        continue
                    else:
                        passed = True
                        break
                assert passed, "No adequate window found"
                m_window = m[window]
                mean_window = mean(m_window)
                std_window = std(m_window)

                # now we're ready to sample that distribution and create our point
                new = (random.normal(loc=mean_window, scale=std_window, size=1))[0]
                data = append(data, new)
                rms = append(rms, std_window)
            period_folded_model_file = file("period_folded_model.txt", "w")
            #             model_file = file("model.txt", "w")
            for n in range(len(t_fold_model)):
                period_folded_model_file.write("%f\t%f\t%f\n" % (t_fold_model[n], data[n], rms[n]))
            #                 model_file.write("%f\t%f\t%f\n" % (available[n], data[n], rms[n]))
            #             model_file.close()
            period_folded_model_file.close()
            return {"flux": data, "rms": rms}

    def __init__(self, mapsize=10.e0, pixels=1024, cosm=Cosmology()):
        """
    Constructor.
    Default will create a 10x10 degree FOV with 1024 pixels and
    a WMAP7 Cosmology.
    """
        self.cosm = cosm

        self.mapsize_deg = mapsize  # map size in np.real domain
        self.mapsize_rad = np.deg2rad(self.mapsize_deg)
        self.fsky = (mapsize**2.e0) / 41253.e0
        self.npix = pixels

        self.Fmapsize = 1.e0 / self.mapsize_rad  # map size in Fourier domain

        self.pixsize = self.mapsize_rad / self.npix  # pixel size in radians
        self.pixsize_deg = self.mapsize_deg / self.npix  # pixel size in degrees

        self.mapaxis = (
            (self.mapsize_deg / self.npix) *  # range of axes
            np.arange(-self.mapsize_deg / 2.e0, self.mapsize_deg / 2.e0, 1))
        self.Fmapaxis = 1.e0 / self.mapaxis

        self.krange = (
            self.Fmapsize *  # define k-space
            np.arange(-self.npix / 2.e0, self.npix / 2.e0, 1))
        self.kx, self.ky = np.meshgrid(self.krange, self.krange)
        self.modk = sqrt(self.kx**2.e0 + self.ky**2.e0)

        self.Txy = random.normal(
            0, 1, [self.npix, self.npix])  # Gaussian random field
        self.Txy = self.Txy - np.mean(self.Txy)
        self.Fxy = ft.fftshift(ft.fft2(self.Txy))  # Fourier domain GRF
        self.Fxy = self.Fxy / sqrt(np.var(self.Fxy))

        self.build_Pk(self.cosm)  # Get flat-sky P(k) for cosmology
        self.Fxy = self.Fxy * self.Pk  # Apply the power spectrum
        self.Txy = np.real(ft.ifft2(ft.fftshift(self.Fxy)))

        self.ymap = np.zeros([self.npix, self.npix])

 def step(self):
     """ integrate state using simple rectangle rule """
     thrust = float(self.action[0])
     rudder = float(self.action[1])
     h, hdot, v = self.sensors
     rnd = random.normal(0,1.0, size=3)
     
     thrust = min(max(thrust,-1),+2)
     rudder = min(max(rudder,-90),+90)
     drag = 5*h + (rudder**2 + rnd[0])
     force = 30.0*thrust - 2.0*v - 0.02*v*drag + rnd[1]*3.0
     v = v + self.dt*force/self.mass
     v = min(max(v,-10),+40)
     torque = -v*(rudder + h + 1.0*hdot + rnd[2]*10.)
     last_hdot = hdot
     hdot += torque / self.I
     hdot = min(max(hdot,-180),180)
     h += (hdot + last_hdot) / 2.0
     if h>180.: 
         h -= 360.
     elif h<-180.: 
         h += 360.
     self.sensors = (h,hdot,v)

def SmecticInitializer(gridShape, sigma=0.2, seed=None):
    if seed is None:
        seed = 0
    random.seed(seed)

    state = SmecticState.SmecticState(gridShape)
    field = state.GetOrderParameterField()

    Ksq = FourierSpaceTools.FourierSpaceTools(gridShape).kSq.numpy_array()

    for component in field.components:
        temp = random.normal(scale=gridShape[0],size=gridShape)
        ktemp = fft.rfftn(temp)*(sqrt(pi)*sigma)**len(gridShape)*exp(-Ksq*sigma**2/4.)
        field[component] = numpy.real(fft.irfftn(ktemp))

    ## To make seed consistent across grid sizes and convergence comparison
    gridShape = copy.copy(oldgrid)
    if gridShape[0] != 128:
        state = ResizeState(state,gridShape[0],Dim=len(gridShape))

    state = ReformatState(state)
    state.ktools = FourierSpaceTools.FourierSpaceTools(gridShape)
    
    return state 

 def drawRandomWeights(self):
     self.module._setParameters(random.normal(0, expln(self.params), self.module.paramdim))