def __init__(self, servIP="127.0.0.1", ownIP="127.0.0.1", port="21560"):
        self.oldScreenValues = None
        self.view = 0
        self.worldRadius = 400
        
        # Start of mousepointer 
        self.lastx = 0
        self.lasty = 15
        self.lastz = 300
        self.zDis = 1
   
        # Start of cube 
        self.cube = [0.0, 0.0, 0.0]
        self.bmpCount = 0
        self.actCount = 0
        self.calcPhysics = 0
        self.newPic = 1
        self.picCount = 0
        self.target = array([80.0, 0.0, 0.0])
      
        self.centerOfGrav = array([0.0, -2.0, 0.0])
        self.points = ones((8, 3), float)
        self.savePics = False
        self.drawCounter = 0
        self.fps = 25
        self.dt = 1.0 / float(self.fps)

        self.client = UDPClient(servIP, ownIP, port)

    def building_loss(self, ci=None, loss_aus_contents=0):
        damage_states = self.get_building_states()
        total_costs = self.structures.cost_breakdown(ci=ci)

        (structure_state, non_structural_state,
            acceleration_sensitive_state) = damage_states
        (structure_cost, non_structural_cost,
            acceleration_cost, contents_cost) = total_costs

        # hardwired loss for each damage state
        f1 = array((0.02, 0.1, 0.5, 1.0))[newaxis, newaxis, :]
        f2 = array((0.02, 0.1, 0.5, 1.0))[newaxis, newaxis, :]
        f3 = array((0.02, 0.1, 0.3, 1.0))[newaxis, newaxis, :]
        f4 = array((0.01, 0.05, 0.25, 0.5))[newaxis, newaxis, :]
        if loss_aus_contents == 1:
            f4 = f4 * 2  # 100% contents loss if building collapses

        structure_ratio = (f1 * structure_state)  # .sum(axis=-1)
        nsd_ratio = (f2 * non_structural_state)  # .sum(axis=-1)
        accel_ratio = (f3 * acceleration_sensitive_state)  # .sum(axis=-1)
        contents_ratio = (f4 * acceleration_sensitive_state)  # .sum(axis=-1)
        loss_ratio = (structure_ratio, nsd_ratio, accel_ratio, contents_ratio)

        structure_loss = structure_ratio * structure_cost[:, newaxis, newaxis]
        nsd_loss = nsd_ratio * non_structural_cost[:, newaxis, newaxis]
        accel_loss = accel_ratio * acceleration_cost[:, newaxis, newaxis]
        contents_loss = contents_ratio * contents_cost[:, newaxis, newaxis]

        total_loss = (structure_loss, nsd_loss, accel_loss, contents_loss)

        return (loss_ratio, total_loss)

def BowTieH_paramSet( paramSet ):
    """
    The Hamiltonian for a single bowtie, (2*S + 1)^5 states, for a set of d
    values -- this will be a generator equation
    """
    
    A = BowTieAdjacencyDic()
    
    #   Break the paramater set into sub-sets based on the spin value.
    paramSetDic = {}
    paramSet = scipy.array( paramSet )
    
    for row in paramSet:
        sNow = row[0]
        if sNow not in paramSetDic:
            paramSetDic[ sNow ] = []
        
        paramSetDic[ sNow ].append( row )
    
    for S in paramSetDic:
        paramSetDic[S] = scipy.array( paramSetDic[S] )
    
        N = (2*S + 1)**5
        
        sPlusMinus, sMinusPlus, sZZ, sZ2 = HParts( S, A )
        
        for params in paramSetDic[S]:
            S, Jx, Jz, d = params
            H = Jx * .5 * ( sPlusMinus + sMinusPlus ) + Jz * sZZ + d * sZ2
            yield (H, params)

 def objfn_data_to_mesh_project(self, x0, args):
     mesh, Xd, Td = args[0], args[1], args[2]
     mesh.set_variables(x0)
     err = scipy.zeros(Xd.shape[0])
     ind = 0
     for xd in Xd:
         xi1 = mesh.elements[1].project(xd)
         xi2 = mesh.elements[2].project(xd)
         if 0<=xi1<=1:
             xi = xi1
         elif 0<=xi2<=1:
             xi = xi2
         else:
             Xi = scipy.array([xi1, xi1-1, xi2, xi2-1])
             Xi2 = Xi*Xi
             ii = Xi2.argmin()
             xi = Xi[ii]
             if ii < 2:
                 elem = 1
             else:
                 elem = 2
         dx = mesh.elements[elem].evaluate(scipy.array([xi]))[0] - xd
         err[ind] = scipy.sum(dx * dx)
         ind += 1
     return err

	def _initParams_fast(self):
		""" 
		initialize the gp parameters
			1) project Y on the known factor X0 -> Y0
				average variance of Y0 is used to initialize the variance explained by X0
			2) considers the residual Y1 = Y-Y0 (this equivals to regress out X0)
			3) perform PCA on cov(Y1) and considers the first k PC for initializing X
			4) the variance of all other PCs is used to initialize the noise
			5) the variance explained by interaction is set to a small random number 
		"""
		Xd = LA.pinv(self.X0)
		Y0 = self.X0.dot(Xd.dot(self.Y))
		Y1 = self.Y-Y0
		YY = SP.cov(Y1)
		S,U = LA.eigh(YY)
		X = U[:,-self.k:]*SP.sqrt(S[-self.k:])
		a = SP.array([SP.sqrt(Y0.var(0).mean())])
		b = 1e-3*SP.randn(1)
		c = SP.array([SP.sqrt((YY-SP.dot(X,X.T)).diagonal().mean())])
		# gp hyper params
		params = limix.CGPHyperParams()
		if self.interaction:
			params['covar'] = SP.concatenate([a,X.reshape(self.N*self.k,order='F'),SP.ones(1),b])
		else:
			params['covar'] = SP.concatenate([a,X.reshape(self.N*self.k,order='F')])
		params['lik'] = c
		return params

    def __init__(self, x, y, z, bbox=[None] *4, kx=3, ky=3, s=0, bounds_error=True, fill_value=scipy.nan):

        super(RectBivariateSpline, self).__init__( x, y, z, bbox=bbox, kx=kx, ky=ky, s=s)
        self._xlim = scipy.array((x.min(), x.max()))
        self._ylim = scipy.array((y.min(), y.max()))
        self.bounds_error = bounds_error
        self.fill_value = fill_value

def pos2Ray(pos, tokamak, angle=None, eps=1e-6):
    r"""Take in GENIE pos vectors and convert it into TRIPPy rays
    
    Args:
        pos: 4 element tuple or 4x scipy-array
            Each pos is assembled into points of (R1,Z1,RT,phi)

        tokamak: 
            Tokamak object in which the pos vectors are defined.
            
    Returns:
        Ray: Ray object or typle of ray objects.
        
    """

    r1 = scipy.array(pos[0])
    z1 = scipy.array(pos[1])
    rt = scipy.array(pos[2])
    phi = scipy.array(pos[3])

    zt = z1 - scipy.tan(phi)*scipy.sqrt(r1**2 - rt**2)
    angle2  = scipy.arccos(rt/r1)

    if angle is None:
        angle = scipy.zeros(r1.shape)

    pt1 = geometry.Point((r1,angle,z1),tokamak)
    pt2 = geometry.Point((rt,angle+angle2,zt),tokamak)

    output = Ray(pt1,pt2)
    output.norm.s[-1] = eps
    tokamak.trace(output)
    return output

    def __init__(self, x, y, z, a, g, h):
        """
		Construct a Scatterer object, encapsulating the shape and material
		properties of a deformed-cylindrical object with sound speed and
		density similar to the surrounding fluid medium.

		Parameters
		----------
		x, y, z : array-like
			Posiions delimiting the central axis of the scatterer.
		a : array-like
			Array of radii along the centerline of the scatterer.
		g : array-like
			Array of sound speed contrasts (sound speed inside the scatterer
			divided by sound speed in the surrounding medium)
		h : array-like
			Array of density contrasts (density inside the scatterer
			divided by density in the surrounding medium)

		"""
        super(Scatterer, self).__init__()
        self.r = sp.matrix([x, y, z])
        self.a = sp.array(a)
        self.g = sp.array(g)
        self.h = sp.array(h)
        self.cum_rotation = sp.matrix(sp.eye(3))

    def _setnorm(self, input = None, target = None):
        """
        Retrieves normalization info from training data and normalizes data.
        """
        numi = len(self.inno); numo = len(self.outno)
        if input is None and target is None:
            self.inlimits  = array( [[0.15, 0.85]]*numi ) #informative only
            self.outlimits = array( [[0.15, 0.85]]*numo ) #informative only
            self.eni = self.dei = array( [[1., 0.]] * numi )
            self.eno = self.deo = array( [[1., 0.]] * numo )
            self.ded = ones((numo, numi), 'd')
        else:
            input, target = self._testdata(input, target)
            
            # Warn if any input or target node takes a one single value
            # I'm still not sure where to put this check....
            for i, col in enumerate(input.transpose()):
                if max(col) == min(col):
                    print "Warning: %ith input node takes always a single value of %f." %(i+1, max(col))

            for i, col in enumerate(target.transpose()):
                if max(col) == min(col):
                    print "Warning: %ith target node takes always a single value of %f." %(i+1, max(col))
            
            #limits are informative only, eni,dei/eno,deo are input/output coding-decoding
            self.inlimits, self.eni, self.dei = _norms(input, lower=0.15, upper=0.85)
            self.outlimits, self.eno, self.deo = _norms(target, lower=0.15, upper=0.85)
            self.ded = zeros((numo,numi), 'd')
            for o in xrange(numo):
                for i in xrange(numi):
                    self.ded[o,i] = self.eni[i,0] * self.deo[o,0]
            return _normarray(input, self.eni), _normarray(target, self.eno)

    def test_gpkronprod(self):
       # initialize
       covar_c = linear.LinearCF(n_dimensions=self.n_latent)
       covar_r = linear.LinearCF(n_dimensions=self.n_dimensions)
       X0_c = SP.random.randn(self.n_tasks,self.n_latent)
       
       lik = likelihood_base.GaussIsoLik()
       gp = gp_kronprod.KronProdGP(covar_c=covar_c, covar_r=covar_r, likelihood=lik)
       gp.setData(Y=self.Ykronprod['train'],X_r=self.X['train'])
       hyperparams = {'lik':SP.array([0.5]), 'X_c':X0_c, 'covar_r':SP.array([0.5]), 'covar_c':SP.array([0.5]), 'X_r':self.X['train']}
       # check predictions, likelihood and gradients
       gp.predict(hyperparams,Xstar_r=self.X['test'],debugging=True)

       gp._LML_covar(hyperparams,debugging=True)
       gp._LMLgrad_covar(hyperparams,debugging=True)
       gp._LMLgrad_lik(hyperparams,debugging=True)
       gp._LMLgrad_x(hyperparams,debugging=True)
       
       # optimize
       hyperparams = {'lik':SP.array([0.5]), 'X_c':X0_c, 'covar_r':SP.array([0.5]), 'covar_c':SP.array([0.5])}
       opts = {'gradcheck':True}
       hyperparams_opt,lml_opt = optimize_base.opt_hyper(gp,hyperparams,opts=opts)
       Kest = covar_c.K(hyperparams_opt['covar_c'])

       # check predictions, likelihood and gradients
       gp._invalidate_cache() # otherwise debugging parameters are not up to date!
       gp.predict(hyperparams_opt,debugging=True,Xstar_r=self.X['test'])
       gp._LML_covar(hyperparams_opt,debugging=True)
       gp._LMLgrad_covar(hyperparams_opt,debugging=True)
       gp._LMLgrad_lik(hyperparams_opt,debugging=True)
       gp._LMLgrad_x(hyperparams_opt,debugging=True)

def getStepWindow(t, v):
  # return time and voltage vectors during the stimulus period only
  
  # find the point of maximum voltage, and cut off everything afterwards
  maxInd, maxV = max(enumerate(v), key=lambda x: x[1])
  minInd, minV = min(enumerate(v), key=lambda x: x[1])
  if maxV - v[0] > v[0] - minV:
    # this is a positive step
    t = t[:maxInd]
    v = scipy.array(v[:maxInd])
  else:
    # this is a negative step, flip it for now
    t = t[:minInd]
    v = v[0] - scipy.array(v[:minInd])
  
  # re-center time to start at the point of maximum voltage change
  diffV = diff(v)
  dVInd, maxDV = max(enumerate(diffV), key=lambda x: x[1])
  dVInd -= 1
  while diffV[dVInd] > 0:
    dVInd -= 1
  dVInd += 1
  
  t -= t[dVInd]
  v -= v[dVInd]
  
  return t, v, dVInd

def vpd_calc(airtemp= scipy.array([]),\
             rh= scipy.array([])):
    '''
    Function to calculate vapour pressure deficit.

    Parameters:
        - airtemp: measured air temperatures [Celsius].
        - rh: (array of) rRelative humidity [%].
        
    Returns:
        - vpd: (array of) vapour pressure deficits [Pa].
        
    Examples
    --------
    
        >>> vpd_calc(30,60)
        1697.090397862653
        >>> T=[20,25]
        >>> RH=[50,100]
        >>> vpd_calc(T,RH)
        array([ 1168.54009896,     0.        ])
        
    '''
    
    # Test input array/value
    airtemp,rh = _arraytest(airtemp, rh)
    
    # Calculate saturation vapour pressures
    es = es_calc(airtemp)
    eact = ea_calc(airtemp, rh) 
    # Calculate vapour pressure deficit
    vpd = es - eact
    return vpd # in hPa

def ea_calc(airtemp= scipy.array([]),\
            rh= scipy.array([])):
    '''
    Function to calculate actual vapour pressure from relative humidity:
    
    .. math::    
        e_a = \\frac{rh \\cdot e_s}{100}
        
    where es is the saturated vapour pressure at temperature T.

    Parameters:
        - airtemp: array of measured air temperatures [Celsius].
        - rh: Relative humidity [%].

    Returns:
        - ea: array of actual vapour pressure [Pa].

    Examples
    --------
    
        >>> ea_calc(25,60)
        1900.0946514729308

    '''
    
    # Test input array/value
    airtemp,rh = _arraytest(airtemp, rh)

    # Calculate saturation vapour pressures
    es = es_calc(airtemp)
    # Calculate actual vapour pressure
    eact = rh / 100.0 * es
    return eact # in Pa

    def run(self, phase=None, throats=None):
        logger.warning('This algorithm can take some time...')
        conduit_lengths = sp.sum(misc.conduit_lengths(network=self._net,
                                 mode='centroid'), axis=1)
        graph = self._net.create_adjacency_matrix(data=conduit_lengths,
                                                  sprsfmt='csr')

        if phase is not None:
            self._phase = phase
            if 'throat.occupancy' in self._phase.props():
                temp = conduit_lengths*(self._phase['throat.occupancy'] == 1)
                graph = self._net.create_adjacency_matrix(data=temp,
                                                          sprsfmt='csr',
                                                          prop='temp')
        path = spgr.shortest_path(csgraph=graph, method='D', directed=False)

        Px = sp.array(self._net['pore.coords'][:, 0], ndmin=2)
        Py = sp.array(self._net['pore.coords'][:, 1], ndmin=2)
        Pz = sp.array(self._net['pore.coords'][:, 2], ndmin=2)

        Cx = sp.square(Px.T - Px)
        Cy = sp.square(Py.T - Py)
        Cz = sp.square(Pz.T - Pz)
        Ds = sp.sqrt(Cx + Cy + Cz)

        temp = path / Ds

        temp[sp.isnan(temp)] = 0
        temp[sp.isinf(temp)] = 0

        return temp

def GenerateLabels(n):
    " Get proper labeling for output states. "
    # Generate bitstrings
    bitstring = []
    for i in range(0,n+1): 
        bitstring.append(kbits(n, i))
    # Flatten
    bitstring = list(itertools.chain.from_iterable(bitstring))
    # Generate unit vectors
    statelist = []
    poslist = []
    pos = 0
    unit0 = sp.array([1,0])
    unit1 = sp.array([0,1])
    for i in range(len(bitstring)):
        # Construct unit vector corresponding to bitstring
        state = unit1 if (bitstring[i][0] == '1') else unit0
        for j in range(n-1):
            state = sp.kron(state, 
                            unit1 if (bitstring[i][j+1] == '1') else unit0)
        statelist.append(state)
        # Record orientation of unit vector (using position of 1 value)
        for j in range(2**n):
            if (statelist[-1][j]):
                pos = j
                break
        poslist.append(pos)
    # Sort the states
    sortperm = sp.array(poslist).argsort()
    bitstring = [ bitstring[i] for i in sortperm ]

    return bitstring

def ensure_numeric(A, typecode=None):
    """Ensure that sequence is a Numeric array.
    Inputs:
        A: Sequence. If A is already a Numeric array it will be returned
                     unaltered
                     If not, an attempt is made to convert it to a Numeric
                     array
        A: Scalar.   Return 0-dimensional array of length 1, containing that
                     value
        A: String.   Array of ASCII values
        typecode: Numeric type. If specified, use this in the conversion.
                                If not, let Numeric decide

    This function is necessary as array(A) can cause memory overflow.
    """

    if typecode is None:
        if isinstance(A, ndarray):
            return A
        else:
            return array(A)
    else:
        if isinstance(A, ndarray):
            if A.dtype == typecode:
                return array(A)  # FIXME: Shouldn't this just return A?
            else:
                return array(A, typecode)
        else:
            import types
            if isinstance(A, types.StringType):
                return array(A, dtype=int)
            return array(A, typecode)

def buildSharedCrossedNetwork():
    """ build a network with shared connections. Two hidden modules are
    symmetrically linked, but to a different input neuron than the
    output neuron. The weights are random. """
    N = FeedForwardNetwork('shared-crossed')
    h = 1
    a = LinearLayer(2, name = 'a')
    b = LinearLayer(h, name = 'b')
    c = LinearLayer(h, name = 'c')
    d = LinearLayer(2, name = 'd')
    N.addInputModule(a)
    N.addModule(b)
    N.addModule(c)
    N.addOutputModule(d)

    m1 = MotherConnection(h)
    m1.params[:] = scipy.array((1,))

    m2 = MotherConnection(h)
    m2.params[:] = scipy.array((2,))

    N.addConnection(SharedFullConnection(m1, a, b, inSliceTo = 1))
    N.addConnection(SharedFullConnection(m1, a, c, inSliceFrom = 1))
    N.addConnection(SharedFullConnection(m2, b, d, outSliceFrom = 1))
    N.addConnection(SharedFullConnection(m2, c, d, outSliceTo = 1))
    N.sortModules()
    return N

def makesumrule(ptype,plen,ts,lagtype='centered'):
    """ This function will return the sum rule.
        Inputs
            ptype - The type of pulse.
            plen - Length of the pulse in seconds.
            ts - Sample time in seconds.
            lagtype -  Can be centered forward or backward.
        Output
            sumrule - A 2 x nlags numpy array that holds the summation rule.
    """
    nlags = sp.round_(plen/ts)
    if ptype.lower()=='long':
        if lagtype=='forward':
            arback=-sp.arange(nlags,dtype=int)
            arforward = sp.zeros(nlags,dtype=int)
        elif lagtype=='backward':
            arback = sp.zeros(nlags,dtype=int)
            arforward=sp.arange(nlags,dtype=int)
        else:
            arback = -sp.ceil(sp.arange(0,nlags/2.0,0.5)).astype(int)
            arforward = sp.floor(sp.arange(0,nlags/2.0,0.5)).astype(int)
        sumrule = sp.array([arback,arforward])
    elif ptype.lower()=='barker':
        sumrule = sp.array([[0],[0]])
    return sumrule

def parse_text( fname, n ):
    """Parse a text file containing a sentence on each line. 
    Output the file as a list of arrays with integer word indices and
    the corresponding dictionary."""

    dictionary = {}
    inv_dictionary = []
    def get_idx( w ):
        if not w in dictionary:
            dictionary[w] = len(dictionary)
            inv_dictionary.append( w )
        return dictionary[w]

    corpus = []

    f = open( fname )
    for line in f.xreadlines():
        # Read at most n documents
        if len(corpus) > n:
            break
        words = line.split()
        words = array( map( get_idx, words ) )
        if len(words) > 0:
            corpus.append( words )
    f.close()

    return array(corpus), array(inv_dictionary)

def read_as_array(filename):
    """reads audio file as scipy array using gstreamer framework

    return:
        data as scipy array
        duration in seconds
        channels as int
        samplerate
    """
    path = os.path.abspath(os.path.expanduser(filename))
    with GstAudioFile(path) as f:
        samplerate = f.samplerate
        duration = float(f.duration) / 1000000000 # in seconds
        channels = f.channels

        data_left = []
        data_right = []
        data_mono = []
        for s in f:
            # http://docs.python.org/library/struct.html
            # little or big endian is choosen automatically by python
            # if its stereo (2 channels): first one is left channel, second is right channel, third is left channel...
            # every short (h) is 2 bytes long and padding on stereo with two x (xx) for other channel
            if channels == 1:
                data_mono += list(struct.unpack( ("h"*(len(s)/2)), s))
            elif channels == 2:
                data_left += list(struct.unpack( ("hxx"*(len(s)/4)), s))
                data_right += list(struct.unpack( ("xxh"*(len(s)/4)), s))

        if channels == 1:
            data = scipy.array(data_mono)
        elif channels == 2:
            data = scipy.array([data_left, data_right])

    return data, duration, channels, samplerate