def setUp(self):
     # setup params for symbolic computation
     th1 = var('th1')
     l1 = var('L1')
     self.sym_gamma = 0
     self.sym_b = 0
     self.sym_alpha = 0
     self.sym_d = l1
     self.sym_theta = th1
     self.sym_r = 0
     self.sym_tmat = Matrix([
         [cos(th1), -sin(th1), 0, l1],
         [sin(th1), cos(th1), 0, 0],
         [0, 0, 1, 0],
         [0, 0, 0, 1]
     ])
     # setup params for numeric computation
     th1 = pi/2
     l1 = 1
     self.num_gamma = 0
     self.num_b = 0
     self.num_alpha = 0
     self.num_d = l1
     self.num_theta = th1
     self.num_r = 0
     self.num_tmat = Matrix([
         [cos(th1), -sin(th1), 0, l1],
         [sin(th1), cos(th1), 0, 0],
         [0, 0, 1, 0],
         [0, 0, 0, 1]
     ])

def test_var():
    var("a")
    assert a  == Symbol("a")

    var("b bb cc zz _x")
    assert b  == Symbol("b")
    assert bb == Symbol("bb")
    assert cc == Symbol("cc")
    assert zz == Symbol("zz")
    assert _x == Symbol("_x")

    v = var(['d','e','fg'])
    assert d  == Symbol('d')
    assert e  == Symbol('e')
    assert fg == Symbol('fg')

    # check return value
    assert v  == (d, e, fg)

    # see if var() really injects into global namespace
    raises(NameError, "z1")
    make_z1()
    assert z1 == Symbol("z1")

    raises(NameError, "z2")
    make_z2()
    assert z2 == Symbol("z2")

    def eq_sympy(input_dim, output_dim, ARD=False):
        """
        Latent force model covariance, exponentiated quadratic with multiple outputs. Derived from a diffusion equation with the initial spatial condition layed down by a Gaussian process with lengthscale given by shared_lengthscale.

        See IEEE Trans Pattern Anal Mach Intell. 2013 Nov;35(11):2693-705. doi: 10.1109/TPAMI.2013.86. Linear latent force models using Gaussian processes. Alvarez MA, Luengo D, Lawrence ND.

        :param input_dim: Dimensionality of the kernel
        :type input_dim: int
        :param output_dim: number of outputs in the covariance function.
        :type output_dim: int
        :param ARD: whether or not to user ARD (default False).
        :type ARD: bool

        """
        real_input_dim = input_dim
        if output_dim>1:
            real_input_dim -= 1
        X = sp.symbols('x_:' + str(real_input_dim))
        Z = sp.symbols('z_:' + str(real_input_dim))
        scale = sp.var('scale_i scale_j',positive=True)
        if ARD:
            lengthscales = [sp.var('lengthscale%i_i lengthscale%i_j' % i, positive=True) for i in range(real_input_dim)]
            shared_lengthscales = [sp.var('shared_lengthscale%i' % i, positive=True) for i in range(real_input_dim)]
            dist_string = ' + '.join(['(x_%i-z_%i)**2/(shared_lengthscale%i**2 + lengthscale%i_i**2 + lengthscale%i_j**2)' % (i, i, i) for i in range(real_input_dim)])
            dist = parse_expr(dist_string)
            f =  variance*sp.exp(-dist/2.)
        else:
            lengthscales = sp.var('lengthscale_i lengthscale_j',positive=True)
            shared_lengthscale = sp.var('shared_lengthscale',positive=True)
            dist_string = ' + '.join(['(x_%i-z_%i)**2' % (i, i) for i in range(real_input_dim)])
            dist = parse_expr(dist_string)
            f =  scale_i*scale_j*sp.exp(-dist/(2*(lengthscale_i**2 + lengthscale_j**2 + shared_lengthscale**2)))
        return kern(input_dim, [spkern(input_dim, f, output_dim=output_dim, name='eq_sympy')])

    def compute_psi_stats(self):
        # define some normal distributions
        mus = [sp.var("mu_%i" % i, real=True) for i in range(self.input_dim)]
        Ss = [sp.var("S_%i" % i, positive=True) for i in range(self.input_dim)]
        normals = [
            (2 * sp.pi * Si) ** (-0.5) * sp.exp(-0.5 * (xi - mui) ** 2 / Si) for xi, mui, Si in zip(self._sp_x, mus, Ss)
        ]

        # do some integration!
        # self._sp_psi0 = ??
        self._sp_psi1 = self._sp_k
        for i in range(self.input_dim):
            print "perfoming integrals %i of %i" % (i + 1, 2 * self.input_dim)
            sys.stdout.flush()
            self._sp_psi1 *= normals[i]
            self._sp_psi1 = sp.integrate(self._sp_psi1, (self._sp_x[i], -sp.oo, sp.oo))
            clear_cache()
        self._sp_psi1 = self._sp_psi1.simplify()

        # and here's psi2 (eek!)
        zprime = [sp.Symbol("zp%i" % i) for i in range(self.input_dim)]
        self._sp_psi2 = self._sp_k.copy() * self._sp_k.copy().subs(zip(self._sp_z, zprime))
        for i in range(self.input_dim):
            print "perfoming integrals %i of %i" % (self.input_dim + i + 1, 2 * self.input_dim)
            sys.stdout.flush()
            self._sp_psi2 *= normals[i]
            self._sp_psi2 = sp.integrate(self._sp_psi2, (self._sp_x[i], -sp.oo, sp.oo))
            clear_cache()
        self._sp_psi2 = self._sp_psi2.simplify()

def test_symbolic_twoport():
    circuit.default_toolkit = symbolic
    cir = SubCircuit()

    k = symbolic.kboltzmann
    var('R1 R0 C1 w T', real=True, positive=True)
    s = 1j*w

    cir['R0'] = R(1, gnd, r=R0)
    cir['R1'] = R(1, 2, r=R1)
#    cir['C1'] = C(2, gnd, c=C1)

    ## Add an AC source to verify that the source will not affect results
#    cir['IS'] = IS(1, gnd, iac=1) 

    ## Run symbolic 2-port analysis
    twoport_ana = TwoPortAnalysis(cir, Node('1'), gnd, Node('2'), gnd,
                                  noise = True, toolkit=symbolic,
                                  noise_outquantity = 'v')
    result = twoport_ana.solve(freqs=s, complexfreq=True)
    
    ABCD = Matrix(result['twoport'].A)
    ABCD.simplify()

    assert_array_equal(ABCD, np.array([[1 + 0*R1*C1*s, R1],
                                    [(1 + 0*R0*C1*s + 0*R1*C1*s) / R0,  (R0 + R1)/R0]]))

    assert_array_equal(simplify(result['Sin'] - (4*k*T/R0 + 4*R1*k*T/R0**2)), 0)
    assert_array_equal(simplify(result['Svn']), 4*k*T*R1)

def test_coeff2_0():
    var('r, kappa')
    psi = Function("psi")
    g = 1/r**2 * (2*r*psi(r).diff(r, 1) + r**2 * psi(r).diff(r, 2))
    g = g.expand()

    assert g.coeff(psi(r).diff(r, 2)) == 1

def test_var_return():
    raises(ValueError, lambda: var(''))
    v2 = var('q')
    v3 = var('q p')

    assert v2 == Symbol('q')
    assert v3 == (Symbol('q'), Symbol('p'))

def test_geteqsys():
    var('R1 V0 Isat T')
    k = symbolic.kboltzmann
    qelectron = symbolic.qelectron
    

    c = SubCircuit(toolkit=symbolic)
    c['V0'] = VS('net1', gnd, v=V0, toolkit=symbolic)
    c['R1'] = R('net1', 'net2', r=R1)
    c['D1'] = Diode('net2', gnd, IS=Isat, toolkit=symbolic)

    dc = SymbolicDC(c)

    dc.epar.T = T

    eqsys, x = dc.get_eqsys()

    x0, x2, x3 = x

    eqsys_ref = np.array([x3 + x0/R1 - x2/R1, 
                          -Isat*(1 - sympy.exp(qelectron*x2/(T*k))) + x2/R1 - x0/R1, 
                          x0 - V0])

    assert sympy.simplify(eqsys_ref[0] - eqsys_ref[0]) == 0
    assert sympy.simplify(eqsys_ref[0] - eqsys_ref[0]) == 0
    assert sympy.simplify(eqsys_ref[0] - eqsys_ref[0]) == 0

def test_symbolic_noise_kt_over_C():
    pycircuit.circuit.circuit.default_toolkit = symbolic
    cir = SubCircuit(toolkit = symbolic)

    k = symbolic.kboltzmann
    var('r c w w1 V', real=True, positive=True)

    s = I * w

    cir['vs'] = VS(1, gnd, vac=V)
    cir['R'] = R(1, 2, r=r)
    cir['C'] = C(2, gnd, c=c)

    noise = Noise(cir, inputsrc='vs', outputnodes=('2', gnd), 
                  toolkit = symbolic)
    res = noise.solve(s, complexfreq=True)

    T = noise.par.epar.T

    svnout = simplify(res['Svnout'])

    noise_voltage_power = simplify(integrate(svnout, (w, 0, oo)))

    assert_equal(noise_voltage_power, 2*pi*k*T/c)

    assert_equal(simplify(res['gain'] - 1/(1 + s*r*c)), 0)

def S2():
    var("x y z")
    e = (x**sin(x) + y**cos(y) + z**(x + y))**100
    t1 = clock()
    f = e.expand()
    t2 = clock()
    return t2 - t1

def planar2r():
    robo = Robot('planar2r', 2, 2, 3, False, tools.SIMPLE)
    # update geometric params
    params = {
        1: {'sigma': 0, 'mu': 1, 'theta': var('th1')},
        2: {'sigma': 0, 'mu': 1, 'alpha': 0, 'd': var('L1'), 'theta': var('th2')},
        3: {'sigma': 2, 'd': var('L2')}
    }
    robo.update_params('geos', params)
    # update dynamic params
    params = {
        1: {'xx': 0, 'xy': 0, 'xz': 0, 'yy': 0, 'yz': 0, 'ia': 0}
    }
    robo.update_params('dyns', params)
    # update joint params
    params = {
        0: {'qdots': 0, 'qddots': 0, 'torques': 0}
    }
    robo.update_params('misc', params)
    # update gravity vector
    params = {
        0: {'gravity': 0},
        1: {'gravity': 0},
        2: {'gravity': var('G3')},
    }
    robo.update_params('misc', params)
    return robo

def R3():
    var('x y z')
    f = x+y+z
    t1 = clock()
    a = [bool(f==f) for _ in range(10)]
    t2 = clock()
    return t2 - t1

def R7():
    var('x')
    f = x**24+34*x**12+45*x**3+9*x**18 +34*x**10+ 32*x**21
    t1 = clock()
    a = [f.subs({x : random()}) for _ in xrange(10^4)]
    t2 = clock()
    return t2 - t1

def create_layer_sympy_var_and_maps(layer_prop=['z','kv','kh','et', 'mv',
                                                'surz', 'vacz']):
    """Create sympy variables and maps for use with integrating \
    to generate 1d spectral equations.
    
    Each x in layer prop will get a 'top', and 'bot', suffix
    Each 'xtop' will get mapped to 'xt[layer]', each 'xbot' to 'xb[layer] 
    and added to `prop_map` dict.
    Each s in layer prop will get a linear representation added to the 
    `linear_expression` dict:
        x = (xbot-xtop)/(zbot-ztop)*(z-zbot) +xtop
    `prop_map` will also get a 'mi' to m[i], and mj to m[j] map.
    'z, mi, mj' will become global variables
            
    Parameters
    ----------
    layer_prop : ``list`` of ``str``, optional
        label for properties that vary in a layer
        (default is ['z','kv','kh','et', 'mv', 'surz',' vacz'])
    
         
    Returns
    -------
    prop_map : ``dict``
        maps the string version of a variable to the sympy.tensor.IndexedBase 
        version e.g. prop_map['kvtop'] = kvt[layer]
    linear_expressions : ``dict``
        maps the string version of a variable to an expression describing how 
        that varibale varies linearly within a layer    
        
    Examples
    --------
    >>> prop_map, linear_expressions = create_layer_sympy_var_and_maps(layer_prop=['z','kv'])
    >>> prop_map
    {'kvtop': kvt[layer], 'mi': m[i], 'mj': m[j], 'zbot': zb[layer], 'ztop': zt[layer], 'kvbot': kvb[layer]}
    >>> linear_expressions
    {'z': z, 'kv': kvtop + (kvbot - kvtop)*(z - ztop)/(zbot - ztop)}

    """
    #http://www.daniweb.com/software-development/python/threads/111526/setting-a-string-as-a-variable-name    
    m = sympy.tensor.IndexedBase('m')
    i = sympy.tensor.Idx('i')
    j = sympy.tensor.Idx('j')
    layer = sympy.tensor.Idx('layer')
    
    sympy.var('z, mi, mj')
    suffix={'t':'top','b': 'bot'}        
    prop_map = {}
    linear_expressions ={}
    
    prop_map['mi'] = m[i]
    prop_map['mj'] = m[j]
    
    for prop in layer_prop:            
        for s1, s3 in suffix.iteritems():
            vars()[prop + s1] = string_to_IndexedBase(prop + s1)            
            sympy.var(prop + s3)                        
            prop_map[prop + s3] = vars()[prop + s1][layer]        
        linear_expressions[prop]=linear(z, ztop, eval(prop+suffix['t']), zbot, eval(prop+suffix['b']))
    return (prop_map, linear_expressions)

def test__validate_production_function():
    """Testing validation of production function attribute."""

    # production function must have type sym.Basic
    def invalid_F(x, y, l, r, A, kappa, nu, rho, gamma):
        """Valid F must return a SymPy expression."""
        return r * A * kappa * (nu * x**rho + (1 - nu) * (y * (l / r))**rho)**(gamma / rho)

    with nose.tools.assert_raises(AttributeError):
        models.Model('positive', workers, firms, production=invalid_F,
                     params=valid_F_params)

    # production function must share vars with workers and firms
    m, n = sym.var('m, n')
    invalid_F = r * A * kappa * (nu * m**rho + (1 - nu) * (n * (l / r))**rho)**(gamma / rho)

    with nose.tools.assert_raises(AttributeError):
        models.Model('negative', workers, firms, production=invalid_F,
                     params=valid_F_params)

    # production function must depend on r and l
    m, n = sym.var('m, n')
    invalid_F = r * A * kappa * (nu * x**rho + (1 - nu) * (y * (m / n))**rho)**(gamma / rho)

    with nose.tools.assert_raises(AttributeError):
        models.Model('positive', workers, firms, production=invalid_F,
                     params=valid_F_params)

def symList(string, imin, imax):
    outlist = []
    for i in range(imin, imax+1):
        string2 = (string+('%d' % i))
        sympy.var(string2)
        outlist.append(eval(string2))
    return outlist

def test_FunctionWrapper():
    import sympy
    sympy.var("n, m, theta, phi")
    r = sympy.Ynm(n, m, theta, phi)
    s = Integer(2)*r
    assert isinstance(s, Mul)
    assert isinstance(s.args[1]._sympy_(), sympy.Ynm)

    def test_2d5pt_ECMCPU(self):
        store_file = os.path.join(self.temp_dir, 'test_2d5pt_ECMCPU.pickle')
        output_stream = StringIO()

        parser = kc.create_parser()
        args = parser.parse_args(['-m', self._find_file('phinally_gcc.yaml'),
                                  '-p', 'ECMCPU',
                                  self._find_file('2d-5pt.c'),
                                  '-D', 'N', '2000',
                                  '-D', 'M', '1000',
                                  '-vvv',
                                  '--unit=cy/CL',
                                  '--store', store_file])
        kc.check_arguments(args, parser)
        kc.run(parser, args, output_file=output_stream)

        results = pickle.load(open(store_file, 'rb'))

        # Check if results contains correct kernel
        self.assertEqual(list(results), ['2d-5pt.c'])

        # Check for correct variations of constants
        six.assertCountEqual(self,
            [sorted(map(str, r)) for r in results['2d-5pt.c']],
            [sorted(map(str, r)) for r in [((sympy.var('M'), 1000), (sympy.var('N'), 2000))]])

        # Output of first result:
        result = list(results['2d-5pt.c'].values())[0]

        six.assertCountEqual(self, result, ['ECMCPU'])

        ecmd = result['ECMCPU']
        self.assertAlmostEqual(ecmd['T_OL'], 24.8, places=1)
        self.assertAlmostEqual(ecmd['T_nOL'], 20, places=1)

def main():
    var("mu nu rho theta phi sigma tau a t x y z w")

    transform("polar", Matrix([rho*cos(phi), rho*sin(phi)]), [rho, phi])

    transform("cylindrical", Matrix([rho*cos(phi), rho*sin(phi), z]),
            [rho, phi, z])

    transform("spherical",
            Matrix([rho*sin(theta)*cos(phi), rho*sin(theta)*sin(phi),
                rho*cos(theta)]),
            [rho, theta, phi],
            recursive=True
            )

    transform("rotating disk",
            Matrix([t, x*cos(w*t)-y*sin(w*t), x*sin(w*t)+y*cos(w*t), z]),
            [t, x, y, z])

    transform("parabolic",
            Matrix([sigma*tau, (tau**2-sigma**2)/2]),
            [sigma, tau])

    # too complex:
    #transform("bipolar",
    #        Matrix([a*sinh(tau)/(cosh(tau)-cos(sigma)),
    #            a*sin(sigma)/(cosh(tau)-cos(sigma))]),
    #        [sigma, tau]
    #        )

    transform("elliptic",
            Matrix([a*cosh(mu)*cos(nu), a*sinh(mu)*sin(nu)]),
            [mu, nu]
            )

    def test_2d5pt_RooflineIACA(self):
        store_file = os.path.join(self.temp_dir, 'test_2d5pt_RooflineIACA.pickle')
        output_stream = StringIO()

        parser = kc.create_parser()
        args = parser.parse_args(['-m', self._find_file('phinally_gcc.yaml'),
                                  '-p', 'RooflineIACA',
                                  self._find_file('2d-5pt.c'),
                                  '-D', 'N', '4000',
                                  '-D', 'M', '1000',
                                  '-vvv',
                                  '--unit=FLOP/s',
                                  '--store', store_file])
        kc.check_arguments(args, parser)
        kc.run(parser, args, output_file=output_stream)

        results = pickle.load(open(store_file, 'rb'))

        # Check if results contains correct kernel
        self.assertEqual(list(results), ['2d-5pt.c'])

        # Check for correct variations of constants
        six.assertCountEqual(self,
            [sorted(map(str, r)) for r in results['2d-5pt.c']],
            [sorted(map(str, r)) for r in [((sympy.var('M'), 1000), (sympy.var('N'), 4000))]])

        # Output of first result:
        result = list(results['2d-5pt.c'].values())[0]

        six.assertCountEqual(self, result, ['RooflineIACA'])

        roofline = result['RooflineIACA']
        self.assertAlmostEqual(roofline['min performance'], 2900000000.0, places=0)
        self.assertEqual(roofline['bottleneck level'], 3)