def test_nsolve():
    # onedimensional
    from sympy import Symbol, sin, pi
    x = Symbol('x')
    assert nsolve(sin(x), 2) - pi.evalf() < 1e-16
    assert nsolve(Eq(2*x, 2), x, -10) == nsolve(2*x - 2, -10)
    # multidimensional
    x1 = Symbol('x1')
    x2 = Symbol('x2')
    f1 = 3 * x1**2 - 2 * x2**2 - 1
    f2 = x1**2 - 2 * x1 + x2**2 + 2 * x2 - 8
    f = Matrix((f1, f2)).T
    F = lambdify((x1, x2), f.T, modules='mpmath')
    for x0 in [(-1, 1), (1, -2), (4, 4), (-4, -4)]:
        x = nsolve(f, (x1, x2), x0, tol=1.e-8)
        assert mnorm(F(*x),1) <= 1.e-10
    # The Chinese mathematician Zhu Shijie was the very first to solve this
    # nonlinear system 700 years ago (z was added to make it 3-dimensional)
    x = Symbol('x')
    y = Symbol('y')
    z = Symbol('z')
    f1 = -x + 2*y
    f2 = (x**2 + x*(y**2 - 2) - 4*y)  /  (x + 4)
    f3 = sqrt(x**2 + y**2)*z
    f = Matrix((f1, f2, f3)).T
    F = lambdify((x,  y,  z), f.T, modules='mpmath')
    def getroot(x0):
        root = nsolve((f1,  f2,  f3), (x,  y,  z), x0)
        assert mnorm(F(*root),1) <= 1.e-8
        return root
    assert map(round,  getroot((1,  1,  1))) == [2.0,  1.0,  0.0]

def test_nsolve_complex():
    x, y = symbols('x y')

    assert nsolve(x**2 + 2, 1j) == sqrt(2.)*I
    assert nsolve(x**2 + 2, I) == sqrt(2.)*I

    assert nsolve([x**2 + 2, y**2 + 2], [x, y], [I, I]) == Matrix([sqrt(2.)*I, sqrt(2.)*I])
    assert nsolve([x**2 + 2, y**2 + 2], [x, y], [I, I]) == Matrix([sqrt(2.)*I, sqrt(2.)*I])

def test_nsolve_dict_kwarg():
    x, y = symbols('x y')
    # one variable
    assert nsolve(x**2 - 2, 1, dict = True) == \
        [{x: sqrt(2.)}]
    # one variable with complex solution
    assert nsolve(x**2 + 2, I, dict = True) == \
        [{x: sqrt(2.)*I}]
    # two variables
    assert nsolve([x**2 + y**2 - 5, x**2 - y**2 + 1], [x, y], [1, 1], dict = True) == \
        [{x: sqrt(2.), y: sqrt(3.)}]

def test_nsolve_precision():
    x, y = symbols('x y')
    sol = nsolve(x**2 - pi, x, 3, prec=128)
    assert abs(sqrt(pi).evalf(128) - sol) < 1e-128
    assert isinstance(sol, Float)

    sols = nsolve((y**2 - x, x**2 - pi), (x, y), (3, 3), prec=128)
    assert isinstance(sols, Matrix)
    assert sols.shape == (2, 1)
    assert abs(sqrt(pi).evalf(128) - sols[0]) < 1e-128
    assert abs(sqrt(sqrt(pi)).evalf(128) - sols[1]) < 1e-128
    assert all(isinstance(i, Float) for i in sols)

	def age_dist(self, age=1.0):

		# Half-lives from Dunham et al. (2014) and age distribution from C2D (Evans et al. 2009) and Sadavoy et al. 2014
		# Note: Not correct half-lives: would need to be recalculated under the assumption of consecutive decay which is why
		# the calculated age distribution does not match the input criteria

		### Relative population fractions observed in Perseus where t = 1 Myr
		age_frac = numpy.asarray([0.068, 0.211, 0.102, 0.545, 0.075])
		age_0 = 1.0

		lmbda_obs = numpy.zeros(5)
		lmbda = numpy.zeros(5)
		for i in range(0,len(age_frac)): lmbda_obs[i] = log(age_frac[i])/(-age_0)

		# Get values of lambda corresponding to observed age distribution
		# Note that the initial guesses are specific to the Class distribution in Perseus at 1 Myr
		x = Symbol('x')
		lmbda[0] = nsolve(exp(-x * age_0) - age_frac[0], x, [lmbda_obs[0]])
		lmbda[1] = nsolve(lmbda[0]/(x-lmbda[0]) * (exp(-lmbda[0]*age_0) - exp(-x*age_0)) - age_frac[1], x, [lmbda_obs[1]])
		lmbda[2] = nsolve(lmbda[0]*lmbda[1] * (exp(-lmbda[0]*age_0)/(lmbda[1]-lmbda[0])/(x-lmbda[0]) + 
			exp(-lmbda[1]*age_0)/(lmbda[0]-lmbda[1])/(x-lmbda[1]) + 
			exp(-x*age_0)/(lmbda[0]-x)/(lmbda[1]-x))-age_frac[2], x, [6.0])
		lmbda[3] = nsolve(lmbda[0]*lmbda[1]*lmbda[2] * (exp(-lmbda[0]*age_0)/(lmbda[1]-lmbda[0])/(lmbda[2]-lmbda[0])/(x-lmbda[0]) +
			exp(-lmbda[1]*age_0)/(lmbda[0]-lmbda[1])/(lmbda[2]-lmbda[1])/(x-lmbda[1]) + 
			exp(-lmbda[2]*age_0)/(lmbda[0]-lmbda[2])/(lmbda[1]-lmbda[2])/(x-lmbda[2]) + 
			exp(-x*age_0)/(lmbda[0]-x)/(lmbda[1]-x)/(lmbda[2]-x)) - age_frac[3], x, [0.3])
		lmbda[4] = nsolve(lmbda[0]*lmbda[1]*lmbda[2]*lmbda[3] * (exp(-lmbda[0]*age_0)/(lmbda[1]-lmbda[0])/(lmbda[2]-lmbda[0])/(lmbda[3]-lmbda[0])/(x-lmbda[0]) +
			exp(-lmbda[1]*age_0)/(lmbda[0]-lmbda[1])/(lmbda[2]-lmbda[1])/(lmbda[3]-lmbda[1])/(x-lmbda[1]) + 
			exp(-lmbda[2]*age_0)/(lmbda[0]-lmbda[2])/(lmbda[1]-lmbda[2])/(lmbda[3]-lmbda[2])/(x-lmbda[2]) + 
			exp(-lmbda[3]*age_0)/(lmbda[0]-lmbda[3])/(lmbda[1]-lmbda[3])/(lmbda[2]-lmbda[3])/(x-lmbda[3]) + 
			exp(-x*age_0)/(lmbda[0]-x)/(lmbda[1]-x)/(lmbda[2]-x)/(lmbda[3]-x)) - age_frac[4], x, [-0.1])

		# Calculate new fractional populations, setting Class III equal to any leftovers
		self.frac = numpy.zeros(5)
		self.frac[0] = exp(-lmbda[0]*age)
		self.frac[1] = lmbda[0]/(lmbda[1]-lmbda[0]) * (exp(-lmbda[0]*age) - exp(-lmbda[1]*age))
		self.frac[2] = lmbda[0]*lmbda[1] * (exp(-lmbda[0]*age)/(lmbda[1]-lmbda[0])/(lmbda[2]-lmbda[0]) + 
			exp(-lmbda[1]*age)/(lmbda[0]-lmbda[1])/(lmbda[2]-lmbda[1]) + 
			exp(-lmbda[2]*age)/(lmbda[0]-lmbda[2])/(lmbda[1]-lmbda[2]))
		self.frac[3] = lmbda[0]*lmbda[1]*lmbda[2] * (exp(-lmbda[0]*age)/(lmbda[1]-lmbda[0])/(lmbda[2]-lmbda[0])/(lmbda[3]-lmbda[0]) +
			exp(-lmbda[1]*age)/(lmbda[0]-lmbda[1])/(lmbda[2]-lmbda[1])/(lmbda[3]-lmbda[1]) + 
			exp(-lmbda[2]*age)/(lmbda[0]-lmbda[2])/(lmbda[1]-lmbda[2])/(lmbda[3]-lmbda[2]) + 
			exp(-lmbda[3]*age)/(lmbda[0]-lmbda[3])/(lmbda[1]-lmbda[3])/(lmbda[2]-lmbda[3]))
		# self.frac[4] = lmbda[0]*lmbda[1]*lmbda[2]*lmbda[3] * (exp(-lmbda[0]*age)/(lmbda[1]-lmbda[0])/(lmbda[2]-lmbda[0])/(lmbda[3]-lmbda[0])/(lmbda[4]-lmbda[0]) +
		#	exp(-lmbda[1]*age)/(lmbda[0]-lmbda[1])/(lmbda[2]-lmbda[1])/(lmbda[3]-lmbda[1])/(lmbda[4]-lmbda[1]) + 
		#	exp(-lmbda[2]*age)/(lmbda[0]-lmbda[2])/(lmbda[1]-lmbda[2])/(lmbda[3]-lmbda[2])/(lmbda[4]-lmbda[2]) + 
		#	exp(-lmbda[3]*age)/(lmbda[0]-lmbda[3])/(lmbda[1]-lmbda[3])/(lmbda[2]-lmbda[3])/(lmbda[4]-lmbda[3]) + 
		#	exp(-lmbda[4]*age)/(lmbda[0]-lmbda[4])/(lmbda[1]-lmbda[4])/(lmbda[2]-lmbda[4])/(lmbda[3]-lmbda[4]))

		# Assume that everything ends up as Class III; no main sequence. An ok assumption when the interest is on Class 0/I sources
		self.frac[4] = 1. - sum(self.frac[:4])

		return self.frac

def test_nsolve_fail():
    x = symbols('x')
    # Sometimes it is better to use the numerator (issue 4829)
    # but sometimes it is not (issue 11768) so leave this to
    # the discretion of the user
    ans = nsolve(x**2/(1 - x)/(1 - 2*x)**2 - 100, x, 0)
    assert ans > 0.46 and ans < 0.47

def test_increased_dps():
    # Issue 8564
    import mpmath
    mpmath.mp.dps = 128
    x = Symbol('x')
    e1 = x**2 - pi
    q = nsolve(e1, x, 3.0)

    assert abs(sqrt(pi).evalf(128) - q) < 1e-128

def test_nsolve():
    # onedimensional
    x = Symbol('x')
    assert nsolve(sin(x), 2) - pi.evalf() < 1e-15
    assert nsolve(Eq(2*x, 2), x, -10) == nsolve(2*x - 2, -10)
    # Testing checks on number of inputs
    raises(TypeError, lambda: nsolve(Eq(2*x, 2)))
    raises(TypeError, lambda: nsolve(Eq(2*x, 2), x, 1, 2))
    # issue 4829
    assert nsolve(x**2/(1 - x)/(1 - 2*x)**2 - 100, x, 0)  # doesn't fail
    # multidimensional
    x1 = Symbol('x1')
    x2 = Symbol('x2')
    f1 = 3 * x1**2 - 2 * x2**2 - 1
    f2 = x1**2 - 2 * x1 + x2**2 + 2 * x2 - 8
    f = Matrix((f1, f2)).T
    F = lambdify((x1, x2), f.T, modules='mpmath')
    for x0 in [(-1, 1), (1, -2), (4, 4), (-4, -4)]:
        x = nsolve(f, (x1, x2), x0, tol=1.e-8)
        assert mnorm(F(*x), 1) <= 1.e-10
    # The Chinese mathematician Zhu Shijie was the very first to solve this
    # nonlinear system 700 years ago (z was added to make it 3-dimensional)
    x = Symbol('x')
    y = Symbol('y')
    z = Symbol('z')
    f1 = -x + 2*y
    f2 = (x**2 + x*(y**2 - 2) - 4*y) / (x + 4)
    f3 = sqrt(x**2 + y**2)*z
    f = Matrix((f1, f2, f3)).T
    F = lambdify((x, y, z), f.T, modules='mpmath')

    def getroot(x0):
        root = nsolve(f, (x, y, z), x0)
        assert mnorm(F(*root), 1) <= 1.e-8
        return root
    assert list(map(round, getroot((1, 1, 1)))) == [2.0, 1.0, 0.0]
    assert nsolve([Eq(
        f1), Eq(f2), Eq(f3)], [x, y, z], (1, 1, 1))  # just see that it works
    a = Symbol('a')
    assert nsolve(1/(0.001 + a)**3 - 6/(0.9 - a)**3, a, 0.3).ae(
        mpf('0.31883011387318591'))

def model2(data, n, price1, price2):
    sample = float(data.shape[0])
    shares = np.histogram(data[:, 1], bins=np.arange(1.0, n + 2.0))[0] / sample
    x = Symbol("x")
    y = Symbol("y")
    eq1 = shares[1] * price1 - (price1 * x) / (1 + x + y)
    eq2 = shares[2] * price2 - (price2 * y) / (1 + x + y)
    sl = nsolve(
        (eq1, eq2),
        (x, y),
        (0.5, 0.5),
        set=True,
        rational=True,
        manual=True,
        implicit=True,
        simplify=True,
        numerical=True,
    )
    # print('prices: 0\t' + str(price1) + '\t' + str(price2))
    # print('shares: ' + str(shares[0]) + '\t' + str(shares[1])  + '\t' + str(shares[2]))
    # print('d = ' + str(sl[0]) + '\t' + str(sl[1]))
    V = np.array([0, sl[0], sl[1]])
    return V

def test_issue_6408_fail():
    x, y = symbols('x y')
    assert nsolve(Integral(x*y, (x, 0, 5)), y, 2) == 0.0

def test_issue_6408():
    x = Symbol('x')
    assert nsolve(Piecewise((x, x < 1), (x**2, True)), x, 2) == 0.0

 def getroot(x0):
     root = nsolve(f, (x, y, z), x0)
     assert mnorm(F(*root), 1) <= 1.e-8
     return root

 def getroot(x0):
     root = nsolve((f1, f2, f3), (x, y, z), x0)
     assert mnorm(F(*root), 1) <= 1.0e-8
     return root

def exactSol(x,t=0.2,x0=0,ql=[1.,1.,0.],qr=[0.25,0.1,0.],gamma=1.4):
	''' Gives exact solution to Sod shock tube problem in order to check for accuracies and compare with computational solution.
	Exact solution is given at x points.
	Algorith taken from http://www.phys.lsu.edu/~tohline/PHYS7412/sod.html
	Ref. Sod, G. A. 1978, Journal of Computational Physics, 27, 1-31.
	'''
	
	#Import stuff
	from sympy.solvers import nsolve
	from sympy import Symbol
	
	#Initiate stuff
	
	shape=x.shape
	x=x.flatten()
	p1 = Symbol('p1')
	[rol,pl,vl]=ql
	[ror,pr,vr]=qr
	
	#Calculate wave velocities and values
	
	cleft=(gamma*pl/rol)**(0.5)
	cright=(gamma*pr/ror)**(0.5)
	m=((gamma-1)/(gamma+1))**0.5
	eq=((2*(gamma**0.5))/(gamma-1)*(1-(p1**((gamma-1)/2/gamma))))-((p1-pr)*(((1-(m**2))**2)*((ror*(p1+(m*m*pr)))**(-1)))**(0.5))
	ppost=float(nsolve(eq,p1,0.))
	rpostByrright=((ppost/pr)+(m*m))/(1+((ppost/pr)*(m*m)))
	vpost=(2*(gamma**0.5))/(gamma-1)*(1-(ppost**((gamma-1)/2/gamma)))
	romid=((ppost/pl)**(1/gamma))*rol
	vshock=vpost*(rpostByrright)/(rpostByrright-1)
	ropost=rpostByrright*ror
	pmid=ppost
	vmid=vpost

	#Calculate locations
	x1=x0-(cleft*t)
	x3=x0+(vpost*t)
	x4=x0+(vshock*t)
	ro=[]
	p=[]
	v=[]
	for i in x:
		csound=((m*m)*(x0-i)/t)+((1-(m*m))*cleft)
		vinst=(1-(m*m))*(((i-x0)/t)+cleft)
		roinst=rol*((csound/cleft)**(2/(gamma-1)))
		pinst=pl*((roinst/rol)**(gamma))
		if i<x1:
			ro.append(rol)
			p.append(pl)
			v.append(vl)
		elif (i>=x4):
			ro.append(ror)
			p.append(pr)
			v.append(vr)
		elif (i<x4) and (i>=x3):
			ro.append(ropost)
			p.append(ppost)
			v.append(vpost)
		elif (i<x3) and (((roinst>rol) and (roinst<romid)) or ((roinst<rol) and (roinst>romid))):
			ro.append(roinst)
			p.append(pinst)
			v.append(vinst)
		else:
			ro.append(romid)
			p.append(pmid)
			v.append(vmid)
			
	#Reshape solutions
	ro=np.array(ro).reshape(shape)
	v=np.array(v).reshape(shape)
	p=np.array(p).reshape(shape)
	
	#calculate conserved variables
	rou = ro*v
	ener=p/(gamma-1)/ro

	return([ro,v,p,ener])

T = 297
gamma = mu*v**2/(c.R*T)
Sg = 2*mu*v/(c.R*T)*Sv

Ck = (7.0-5.0*gamma)/(2*(gamma-1))*c.R
Sck = c.R/(gamma-1)**2*Sg
print "Cк:", Ck
print "Погрешность Cк:", Sck

a = (7.0-5.0*gamma)/(2*(gamma-1))
# Sa = Sg/(gamma-1)**2
print "a:", a
print "a/2:", a/2

y = Symbol('y')
# последнее число это приблизительный хk
# лучше всего посмотреть на картинке приблизительный х,
# которому соответствует Y=a/2
xk = nsolve(2*y**2*exp(-y)/(1-exp(-y))**2-a, 4.5)
print "xk:", xk
Sxk = 0.05

nu = xk*c.k*T/c.h
Snu = c.k*T/c.h*Sxk
print "Частота:", nu
print "Погрешность частоты:", Snu

Th = c.h*nu/c.k
Sth = T*Sxk
print "Характеристическая темп.:", Th
print "Погрешность хар-кой темп.:", Sth

def test_nsolve_denominator():
    x = symbols('x')
    # Test that nsolve uses the full expression (numerator and denominator).
    ans = nsolve((x**2 + 3*x + 2)/(x + 2), -2.1)
    # The root -2 was divided out, so make sure we don't find it.
    assert ans == -1.0

def test_nsolve_rational():
    x = symbols('x')
    assert nsolve(x - Rational(1, 3), 0, prec=100) == Rational(1, 3).evalf(100)

def test_issue_14950():
    x = Matrix(symbols('t s'))
    x0 = Matrix([17, 23])
    eqn = x + x0
    assert nsolve(eqn, x, x0) == -x0
    assert nsolve(eqn.T, x.T, x0.T) == -x0

args = "x,y,z"

# really sloppy handling, use argparse
if len(sys.argv) > 1:
    try:
        N = int(sys.argv[1])
    except:
        raise Exception("argument must be an integer N (partition size in RK4)")

else:
    N = 4

# just set verbosity based on size of N
v = (N < 5)
x0 = np.array([[1,1,1]]).T

x_sol = cm_rk4(case_one, args, x0, N, verbose=v)

# now just do some data
F, X = parse_symbolic(case_one, args)

# get F value at x*
F1 = F.subs(list(zip(X, x_sol)))

err = F1.norm() # wow, these methods are great
print("||F(x*)|| = ", err)

x_check = nsolve(F, X, (1,1,1))
print("solution according to sympy.solvers.nsolve:")
print(x_check.T)