def forward_pass(self, xL, xR):

		bs = sp.ones((1,xL.shape[1]), dtype=float)

		# First Layer
		xLb = sp.vstack([xL, bs])
		xRb = sp.vstack([xR, bs])

		a1L = sp.dot(self.w1l, xLb)
		a1R = sp.dot(self.w1r, xRb)
		
		z1L = sp.tanh(a1L)
		z1R = sp.tanh(a1R)

		# Second Layer		
		z1Lb = sp.vstack([z1L, bs])
		z1LRb = sp.vstack([z1L, z1R, bs])
		z1Rb = sp.vstack([z1R, bs])

		a2L = sp.dot(self.w2l, z1Lb)
		a2LR = sp.dot(self.w2lr, z1LRb)
		a2R = sp.dot(self.w2r, z1Rb)

		z2 = a2LR*self.sigmoid(a2L)*self.sigmoid(a2R)

		# Third Layer
		z2b = sp.vstack([z2, bs])
		a3 = sp.dot(self.w3, z2b)

		return a1L, a1R, a2L, a2LR, a2R, a3, z1Lb, z1LRb, z1Rb, z2b, xLb, xRb

    def fgrad_y(self, y, psi, return_precalc = False):
        """
	gradient of f w.r.t to y ([N x 1])
	returns: Nx1 vector of derivatives, unless return_precalc is true,
	then it also returns the precomputed stuff
	"""

	mpsi = psi.copy()
	mpsi[:,0:2] = SP.exp(mpsi[:,0:2])
	s = SP.zeros((len(psi), y.shape[0], y.shape[1]))
	r = SP.zeros((len(psi), y.shape[0], y.shape[1]))	
	d = SP.zeros((len(psi), y.shape[0], y.shape[1]))
	
	grad = 1
	for i in range(len(mpsi)):
	    a,b,c = mpsi[i]
	    s[i] = b*(y+c)
	    r[i] = SP.tanh(s[i])
	    d[i] = 1 - r[i]**2    
	    grad += a*b*d[i]

        #vectorized version
        S = (mpsi[:,1]*(y + mpsi[:,2])).T
        R = SP.tanh(S)
        D = 1-R**2
        GRAD = (1+(mpsi[:,0:1]*mpsi[:,1:2]*D).sum(axis=0))[:,SP.newaxis]

        if return_precalc:
            return GRAD,S,R,D
	    #return grad, s, r, d
	
	return grad

def DP_lib(compuesto, T):
    """Dohrn, R.; Prausnitz, J.M. A simple perturbation term for the Carnahan-Starling equation of state. Fluid Phase Eq. 1990, 61, 53."""
    Tr = T / compuesto.Tc
    if Tr >= 1:
        m = 1
    else:
        m = 0

    a1 = 0.367845 + 0.055966 * compuesto.f_acent
    a2 = -1 ** m * (0.604709 + 0.008477 * compuesto.f_acent)
    ac = 0.550408 * R_atml ** 2 * compuesto.Tc ** 2 / compuesto.Pc.atm
    a = ac * (a1 * tanh(a2 * abs(Tr - 1) ** 0.7) + 1)

    b1 = 0.356983 - 0.190003 * compuesto.f_acent
    b2 = -1 ** m * (1.37 - 1.898981 * compuesto.f_acent)
    bc = 0.187276 * R_atml * compuesto.Tc / compuesto.Pc.atm
    b = bc * (b1 * tanh(b2 * abs(log(Tr)) ** 0.8) + 1)
    sigma = (3 * b / 2 / pi / Avogadro) ** (1.0 / 3)
    return a, b, sigma

    def DP_Z(self, T, P):
        """Factor de compresibilidad según la ecuación de estado de Dohrn-Prausnith"""
        V = self.DP_V(T, P) * self.peso_molecular
        return P * V / R_atml / T

    def DP_V(self, T, P):
        """Volumen según el modelo de Dohrn-Prausnith"""
        a, b, sigma = self.DP_lib(T, P)
        D = sigma
        E = sigma ** 2
        F = sigma ** 3

        if self.Fase(T, P) == "gas":
            V = 25
        else:
            V = 0.5

        def Vm(V):
            nu = b / 4 / V
            Zref = (
                1
                + (3 * D * E / F - 2) * nu
                + (3 * E ** 3 / F ** 2 - 3 * D * E / F + 1) * nu ** 2
                - E ** 3 / F ** 2 * nu ** 3
            ) / (1 - nu) ** 3
            Zpert = a / R_atml / T / V * (1 - 1.41 * b / V / 4 + 5.07 * (b / V / 4) ** 2)
            return V - (Zref + Zpert) * R_atml * T / P

        v = fsolve(Vm, V)

        return unidades.SpecificVolume(v / self.peso_molecular)

    def DP_RhoG(self, T, P):
        """Método para el cálculo de la densidad de gases haciendo uso de la ecuación de estado de Dohrn-Prausnith"""
        z = self.DP_Z(T, P)
        return unidades.SpecificVolume(P / z / R_atml / T * self.peso_molecular)

def double_tanh_warp(x, n, lcore, lmid, ledge, la, lb, xa, xb):
    r"""Implements a sum-of-tanh warping function and its derivative.

    .. math::

        l = a\tanh\frac{x-x_a}{l_a} + b\tanh\frac{x-x_b}{l_b}

    Parameters
    ----------
    x : float or array of float
        Locations to evaluate the function at.
    n : int
        Derivative order to take. Used for ALL of the points.
    lcore : float
        Core length scale.
    lmid : float
        Intermediate length scale.
    ledge : float
        Edge length scale.
    la : positive float
        Transition of first tanh.
    lb : positive float
        Transition of second tanh.
    xa : float
        Transition of first tanh.
    xb : float
        Transition of second tanh.

    Returns
    -------
    l : float or array
        Warped length scale at the given locations.

    Raises
    ------
    NotImplementedError
        If `n` > 1.
    """
    a, b, c = scipy.dot([[-0.5, 0, 0.5], [0, 0.5, -0.5], [0.5, 0.5, 0]],
                        [[lcore], [ledge], [lmid]])
    a = a[0]
    b = b[0]
    c = c[0]
    if n == 0:
        return a * scipy.tanh((x - xa) / la) + b * scipy.tanh((x - xb) / lb) + c
    elif n == 1:
        return (a / la * (scipy.cosh((x - xa) / la))**(-2.0) +
                b / lb * (scipy.cosh((x - xb) / lb))**(-2.0))
    else:
        raise NotImplementedError("Only derivatives up to order 1 are supported!")

def plot_neural_net(X, y, clf, segment=False):
    values = X
    pca_plot(X, y, save="00_conn.png", segment=segment)
    counter = 1
    for i, layer in enumerate(clf.nn.modulesSorted):
        name = layer.__class__.__name__
        if name == "BiasUnit":
            continue

        try:
            conn = clf.nn.connections[layer][0]
        except IndexError:
            continue

        if "Linear" not in name:
            if "Sigmoid" in name:
                add = "sigmoid"
                values = sigmoid(values)
            elif "Tanh" in name:
                add = "tanh"
                values = tanh(values)
            pca_plot(values, y, save="%02d_conn_%s.png" % (counter, add), segment=segment)
            counter += 1
        shape = (conn.outdim, conn.indim)
        temp = numpy.dot(numpy.reshape(conn.params, shape), values.T)
        pca_plot(temp.T, y, save="%02d_conn.png" % counter, segment=segment)
        counter += 1
        values = temp.T

def tanh_warp_arb(X, l1, l2, lw, x0):
    r"""Warps the `X` coordinate with the tanh model

    .. math::

        l = \frac{l_1 + l_2}{2} - \frac{l_1 - l_2}{2}\tanh\frac{x-x_0}{l_w}

    Parameters
    ----------
    X : :py:class:`Array`, (`M`,) or scalar float
        `M` locations to evaluate length scale at.
    l1 : positive float
        Small-`X` saturation value of the length scale.
    l2 : positive float
        Large-`X` saturation value of the length scale.
    lw : positive float
        Length scale of the transition between the two length scales.
    x0 : float
        Location of the center of the transition between the two length scales.

    Returns
    -------
    l : :py:class:`Array`, (`M`,) or scalar float
        The value of the length scale at the specified point.
    """
    if isinstance(X, scipy.ndarray):
        if isinstance(X, scipy.matrix):
            X = scipy.asarray(X, dtype=float)
        return 0.5 * ((l1 + l2) - (l1 - l2) * scipy.tanh((X - x0) / lw))
    else:
        return 0.5 * ((l1 + l2) - (l1 - l2) * mpmath.tanh((X - x0) / lw))

    def draw_eccentricities(self, eccentricity='flat', emin=0., emax='tidal'):
        """Draw a new eccentricity distribution.

        Either provide an array (with the same length as the number of binaries) or choose from:
        - 'flat': Flat distribution between `emin` and `emax`.
        - 'thermal`: f(e) ~ e^2 between `emin` and `emax`.

        If `emax` is set to 'tidal' it will depend on the period of the binary to mimic tidal circulirization through:
        emax = max(emin, 0.5 * (0.95 + tanh(0.6 * log_10(period in days) - 1.7)))
        If this setting is chosen the eccentricities will have to be redrawn if the periods are redrawn.

        Arguments:
        - `eccentricity`: New eccentricity distribution (array or name of distriution).
        - `emin`: Minimum eccentricity (default: 0.)
        - `emax`: Maximum eccentricity (default: set by tidal circularization)
        """
        nbinaries = self.size
        if emax == 'tidal':
            emax =  .5 * (0.95 + sp.tanh(0.6 * sp.log10(self['period'] * 365.25) - 1.7))
            emax[emax < emin] = emin
        if eccentricity == 'flat':
            self['eccentricity'] = sp.random.rand(nbinaries) * (emax - emin) + emin
        elif eccentricity == 'thermal':
            self['eccentricity'] = sp.sqrt(sp.random.rand(nbinaries) * (emax ** 2. - emin ** 2.) + emin ** 2.)
        elif isinstance(eccentricity, basestring):
            raise ValueError("Eccentricity distribution '%s' not found" % eccentricity)
        else:
            self['eccentricity'] = eccentricity

    def calculate_dwf(self):
        """
        Calculates Debye-Waller factor according to
        Sears and Shelley Acta Cryst. A 47, 441 (1991)
        """
        run = self.vanaws.getRun()
        nhist = self.vanaws.getNumberHistograms()
        thetasort = np.empty(nhist)  # half of the scattering angle, in radians
        for i in range(nhist):
            det = self.vanaws.getDetector(i)
            thetasort[i] = 0.5 * self.vanaws.detectorTwoTheta(det)

        # T in K
        temperature = self.get_temperature()
        # Wavelength, Angstrom
        wlength = float(run.getLogData('wavelength').value)
        # Vanadium mass, kg
        mass_vana = 0.001*self.Mvan/sp.constants.N_A
        temp_ratio = temperature/self.DebyeT

        if temp_ratio < 1.e-3:
            integral = 0.5
        else:
            integral = \
                integrate.quad(lambda x: x/sp.tanh(0.5*x/temp_ratio), 0, 1)[0]

        msd = 3.*sp.constants.hbar**2 / \
            (2.*mass_vana*sp.constants.k * self.DebyeT)*integral*1.e20
        return np.exp(-msd*(4.*sp.pi*sp.sin(thetasort)/wlength)**2)

 def astroOut(self):
     for astro in self.astros:
         astro.act = tanh(astro.act)
         for syn in astro.syns:
             i, j = syn
             self.astroOuts[i][j] = astro.act * self.syn_wi[i][j]  
     return self.astroOuts.copy()

def calculateCPandRT(stim_strengths, k, A, t1_nondt, t2_nondt):
    """
    Given a diffusion to bound model with flat bounds, compute the following
    as a function of stimulus strength:
        1) Probability of t1 choice (correct choice for positive strengths)
        2) Mean reaction times for t1 and t2 choices

    Parameters
    ----------
    stim_strengths : array of unique stimulus strengths

    k : float, proportionality constant between stim_strength and drift rate

    A : float, bound

    t1_nondt: float, non decision time for making t1 choice

    t2_nondt: float, non decision time for making t2 choice

    Returns
    -------
    out : dictionary with keys stim_strength, p_t1, t1_meanrt, t2_meanrt

    """

    # compute probability of t1 response given k and A
    p = 1 / (1 + np.exp(-2 * A * k * stim_strengths))

    # compute mean response times for t1 and t2
    t1_meanrt = np.zeros(len(stim_strengths))
    t2_meanrt = np.zeros(len(stim_strengths))
    for i in range(len(stim_strengths)):
        if stim_strengths[i] == 0:
            t1_meanrt[i] = A ** 2 + t1_nondt
            t2_meanrt[i] = A ** 2 + t2_nondt
        else:
            t1_meanrt[i] = A / (k * stim_strengths[i]) * \
                sp.tanh(A * k * stim_strengths[i]) + t1_nondt
            t2_meanrt[i] = A / (k * stim_strengths[i]) * \
                sp.tanh(A * k * stim_strengths[i]) + t2_nondt

    out = {}
    out['stim_strengths'] = stim_strengths
    out['p_t1'] = p
    out['t1_meanrt'] = t1_meanrt
    out['t2_meanrt'] = t2_meanrt
    return out

def grate(k,L):
    l2 = (np.pi**2)/(4*(L**2))
    mu = np.sqrt(k**2 + l2) 
    th = sp.tanh(mu/2.)
    cth = 1./th
    sigma = (k/mu)*np.sqrt( (mu/2.-cth)*(th-mu/2.)  )
    sigma[np.isnan(sigma)] = 0.
    return sigma

 def activateHid(self):
     # self.ah[:] = tanh(sum(self.wi.T * self.ai, axis=1))
     for j in range(self.hiddim):
         s = 0
         for i in range(self.indim):
             s += sum(self.wi[i][j] * (self.ai + self.astroOuts[i][j]))
             # print self.astroOuts
         self.ah[j] = tanh(s)

def TempProfile(z,T0=1000.,z0=100.):
    """This function creates a tempreture profile for test purposes."""
    zall = (z-z0)*2.*sp.exp(1)/400. -sp.exp(1)
    atanshp = (sp.tanh(zall)+1.)/2
    Te = 1700*atanshp+T0
    Ti = 500*atanshp+T0

    return (Te,Ti)

 def performAction(self, action):
     #Filtered mapping towards performAction of the underlying environment
     #The standard Johnnie task uses a PID controller to controll directly angles instead of forces
     #This makes most tasks much simpler to learn
     isJoints=self.env.getSensorByName('JointSensor') #The joint angles
     isSpeeds=self.env.getSensorByName('JointVelocitySensor') #The joint angular velocitys
     act=(action+1.0)/2.0*(self.env.cHighList-self.env.cLowList)+self.env.cLowList #norm output to action intervall
     action=tanh((act-isJoints-isSpeeds)*16.0)*self.maxPower*self.env.tourqueList #simple PID
     EpisodicTask.performAction(self, action)

def main():
	a = array([1, 1])
	b = array([0, 0])
	x = array([-0.5, 1])
	y = tanh( a * x + b)

	dbma = [-(-0.46 * 2.79), -(0.76*2.42)]


	print("ma byt:", dbma[0], dbma[1])
	print("je:", ipgauss(x, y, a, b)[0])

 def showFor(self, segment, multiplier=1.0):
     nhi = self.takeNumber(self.getSymbol("hi"))
     nlo = self.takeNumber(self.getSymbol("lo"))
     hi = max(nlo + 1.0, nhi)
     lo = min(nhi - 1.0, nlo)
     self.modSymbol("hi", self.makeNumber(hi))
     self.modSymbol("lo", self.makeNumber(lo))
     duration = abs(hi + lo) / 2.0
     duration += (abs(hi - lo) / 2.0) * scipy.tanh(self.tparam)
     print segment + "\r",
     sys.stdout.flush()
     time.sleep(duration * multiplier)

    def ungray(self):
        self.R = self.source[:,:,0]
        self.G = self.source[:,:,1]
        self.B = self.source[:,:,2]

        self.R = filter2D(self.R, -1, self.kernel)
        self.G = filter2D(self.G, -1, self.kernel)
        self.B = filter2D(self.B, -1, self.kernel)

        slope = self.control[ord('s')].val
        self.R = (1.0 + tanh(slope*self.R)) / 2.0
        self.G = (1.0 + tanh(slope*self.G)) / 2.0
        self.B = (1.0 + tanh(slope*self.B)) / 2.0

        self.R = filter2D(self.R, -1, self.gauss)
        self.G = filter2D(self.G, -1, self.gauss)
        self.B = filter2D(self.B, -1, self.gauss)

        self.target[:,:,0] = self.R * self.scale
        self.target[:,:,1] = self.G * self.scale
        self.target[:,:,2] = self.B * self.scale

 def performAction(self, action):
     #Filtered mapping towards performAction of the underlying environment
     #The standard CCRL task uses a PID controller to controll directly angles instead of forces
     #This makes most tasks much simpler to learn
     self.oldAction = action
     #Grasping as reflex depending on the distance to target - comment in for more easy grasping
     if abs(abs(self.dist[:3]).sum())<2.0: action[15]=1.0 #self.grepRew=action[15]*.01
     else: action[15]=-1.0 #self.grepRew=action[15]*-.03
     isJoints=array(self.env.getSensorByName('JointSensor')) #The joint angles
     isSpeeds=array(self.env.getSensorByName('JointVelocitySensor')) #The joint angular velocitys
     act=(action+1.0)/2.0*(self.env.cHighList-self.env.cLowList)+self.env.cLowList #norm output to action intervall
     action=tanh((act-isJoints-0.9*isSpeeds*self.env.tourqueList)*16.0)*self.maxPower*self.env.tourqueList #simple PID
     EpisodicTask.performAction(self, action)

    def dNeurons(self, statevec, t):
    
        # Extract relevant parameters from
        train = training > 0 and t > training and t < training + train_dur

        x_i = statevec[0:Ng]
        w_i = statevec[Ng:2*Ng]

        # generate noise patterns
        exp_noise_amp = 0.1
        if train:
            zeta = np.random.uniform(-exp_noise_amp, exp_noise_amp, 1)
        else:
            zeta = 0.0
            
        # Compute Firing Rates and feedback signals
        r_i = sp.tanh(x_i) + zeta_state(t)
        z_i = np.dot(w_i, r_i) + zeta
        
        # Compute next timestep depending on if training or not
        if train:
        
            dxidt = (-x_i + lambd * np.dot(W_rec, r_i) + np.dot(W_in_left, self.uleft[t]) + np.dot(W_in_right, self.uright[t]) + z_i*W_fb )/tau
            x_new = x_i + dxidt*dt
        
            P = -1.0*sp.power(z_i - self.f_target[t], 2)
            M = 1.0*(P > self.P_avg)
        
            dwdt = eta(t) * (z_i - self.z_avg) * M * r_i
            w_new = w_i + dwdt*dt

            self.P_avg = (1 - (dt/tau_avg)) * self.P_avg + (dt/tau_avg) * P
            self.z_avg = (1 - (dt/tau_avg)) * self.z_avg + (dt/tau_avg) * z_i
        
        else:
        
            dxidt = (-x_i + lambd * np.dot(W_rec, r_i) + np.dot(W_in_left, self.uleft[t]) + np.dot(W_in_right, self.uright[t]) + z_i*W_fb )/tau
                    
            x_new = x_i + dxidt*dt
            dwdt = np.zeros(np.shape(w_i))
            w_new = w_i
            
        #weight change magnitude:
        dwmag = LA.norm(dwdt)
        rmag = LA.norm(r_i)
        if t%10000 == 0: print(str(t) + '    ' + str(z_i) + '   ' + str(train))
        self.zsave[t] = z_i
        self.rsave[t] = rmag
        return np.concatenate((x_new, w_new))

	def Calculate_Sigma(self, data):
		for temperature in xrange(self.n_temp):
			tmp = 1/(sp.tanh(data.energy/(2*K_B*(temperature+1))))
			self.IP_cluster_qu[temperature] = np.sqrt(np.sum((data.IP_cluster_a1*data.IP_cluster_a1)/2 * tmp))
			self.EA_cluster_qu[temperature] = np.sqrt(np.sum((data.EA_cluster_a1*data.EA_cluster_a1)/2 * tmp))
			self.IP_alone_qu[temperature] = np.sqrt(np.sum((data.IP_alone_a1*data.IP_alone_a1)/2 * tmp))
			self.EA_alone_qu[temperature] = np.sqrt(np.sum((data.EA_alone_a1*data.EA_alone_a1)/2 * tmp))
			self.P_plus_qu[temperature] = np.sqrt(np.sum((data.P_plus_a1*data.P_plus_a1)/2 * tmp))
			self.P_minus_qu[temperature] = np.sqrt(np.sum((data.P_minus_a1*data.P_minus_a1)/2 * tmp))

			if temperature+1 == 300:
				self.IP_cluster_qu_300K = 100*((data.IP_cluster_a1*data.IP_cluster_a1)/2 * tmp)/(self.IP_cluster_qu[temperature]*self.IP_cluster_qu[temperature])
				self.EA_cluster_qu_300K = 100*((data.EA_cluster_a1*data.EA_cluster_a1)/2 * tmp)/(self.EA_cluster_qu[temperature]*self.EA_cluster_qu[temperature])
				self.IP_alone_qu_300K = 100*((data.IP_alone_a1*data.IP_alone_a1)/2 * tmp)/(self.IP_alone_qu[temperature]*self.IP_alone_qu[temperature])
				self.EA_alone_qu_300K = 100*((data.EA_alone_a1*data.EA_alone_a1)/2 * tmp)/(self.EA_alone_qu[temperature]*self.EA_alone_qu[temperature])
				self.P_plus_qu_300K = 100*((data.P_plus_a1*data.P_plus_a1)/2 * tmp)/(self.P_plus_qu[temperature]*self.P_plus_qu[temperature])
				self.P_minus_qu_300K = 100*((data.P_minus_a1*data.P_minus_a1)/2 * tmp)/(self.P_minus_qu[temperature]*self.P_minus_qu[temperature])
				
			tmp = 2*K_B*(temperature+1)/data.energy
			self.IP_cluster_cl[temperature] = np.sqrt(np.sum((data.IP_cluster_a1*data.IP_cluster_a1)/2 * tmp))
			self.EA_cluster_cl[temperature] = np.sqrt(np.sum((data.EA_cluster_a1*data.EA_cluster_a1)/2 * tmp))
			self.IP_alone_cl[temperature] = np.sqrt(np.sum((data.IP_alone_a1*data.IP_alone_a1)/2 * tmp))
			self.EA_alone_cl[temperature] = np.sqrt(np.sum((data.EA_alone_a1*data.EA_alone_a1)/2 * tmp))
			self.P_plus_cl[temperature] = np.sqrt(np.sum((data.P_plus_a1*data.P_plus_a1)/2 * tmp))
			self.P_minus_cl[temperature] = np.sqrt(np.sum((data.P_minus_a1*data.P_minus_a1)/2 * tmp))
			
			if temperature+1 == 300:
				self.IP_cluster_cl_300K = 100*((data.IP_cluster_a1*data.IP_cluster_a1)/2 * tmp)/(self.IP_cluster_cl[temperature]*self.IP_cluster_cl[temperature])
				self.EA_cluster_cl_300K = 100*((data.EA_cluster_a1*data.EA_cluster_a1)/2 * tmp)/(self.EA_cluster_cl[temperature]*self.EA_cluster_cl[temperature])
				self.IP_alone_cl_300K = 100*((data.IP_alone_a1*data.IP_alone_a1)/2 * tmp)/(self.IP_alone_cl[temperature]*self.IP_alone_cl[temperature])
				self.EA_alone_cl_300K = 100*((data.EA_alone_a1*data.EA_alone_a1)/2 * tmp)/(self.EA_alone_cl[temperature]*self.EA_alone_cl[temperature])
				self.P_plus_cl_300K = 100*((data.P_plus_a1*data.P_plus_a1)/2 * tmp)/(self.P_plus_cl[temperature]*self.P_plus_cl[temperature])
				self.P_minus_cl_300K = 100*((data.P_minus_a1*data.P_minus_a1)/2 * tmp)/(self.P_minus_cl[temperature]*self.P_minus_cl[temperature])
				
		self.IP_cluster_L = np.sum((data.IP_cluster_a1*data.IP_cluster_a1)/(2*data.energy))
		self.EA_cluster_L = np.sum((data.EA_cluster_a1*data.EA_cluster_a1)/(2*data.energy))
		self.IP_alone_L = np.sum((data.IP_alone_a1*data.IP_alone_a1)/(2*data.energy))
		self.EA_alone_L = np.sum((data.EA_alone_a1*data.EA_alone_a1)/(2*data.energy))
		self.P_plus_L = np.sum((data.P_plus_a1*data.P_plus_a1)/(2*data.energy))
		self.P_minus_L = np.sum((data.P_minus_a1*data.P_minus_a1)/(2*data.energy))
		
		self.IP_cluster_G2 = np.sum((data.IP_cluster_a1*data.IP_cluster_a1)/(2))
		self.EA_cluster_G2 = np.sum((data.EA_cluster_a1*data.EA_cluster_a1)/(2))
		self.IP_alone_G2 = np.sum((data.IP_alone_a1*data.IP_alone_a1)/(2))
		self.EA_alone_G2 = np.sum((data.EA_alone_a1*data.EA_alone_a1)/(2))
		self.P_plus_G2 = np.sum((data.P_plus_a1*data.P_plus_a1)/(2))
		self.P_minus_G2 = np.sum((data.P_minus_a1*data.P_minus_a1)/(2))