def solution(pars):
    (a2,a1,a0),(t1,t2),(k1,k2),(y0,yd0)=pars
    a2,a1,a0 = map(int,[a2,a1,a0])
    den = a2*s**2+a1*s+a0
    V1 = apart(1/s/den).as_ordered_terms()            # const.
    V2 = apart(1/s**2/den).as_ordered_terms()         # lin.
    Iv1,Iv2 = ilt_pfe(V1),ilt_pfe(V2)
    LIC  = (a2*(s*y0+yd0) + a1*y0)/den
    y1 = ILT(LIC, s,t)
    y2 = k1*Iv2                                       #        k1,  t
    y3 = (k2-k1)*Iv2.subs(t,t-t1)*H(t-t1)             #  k2-k1,(t-t1)*H(t-t1)
    y4 = (t1*(k2-k1)-k2*t2)*Iv1.subs(t,t-t1)*H(t-t1)  # t1*(k2-k1)-k2*t2,H(t-t1), nulove pre spojite
    y5 = -k2*Iv2.subs(t,t-t2)*H(t-t2)                 #  -k2, (t-t2)*H(t-t2)
    y = expand_mul(simplify(y1+y2+y3+y4+y5))
    return y

def expand(expression, variable = None):
    expression = expand_(expression)
    if isinstance(expression, Basic) and expression.is_rational_function():
        if variable is None:
           variable = extract_var(expression)
        return apart(expression, variable)
    return expression

 def apply(self, expr, var, evaluation):
     'Apart[expr_, var_Symbol]'
     
     expr_sympy = expr.to_sympy()
     var_sympy = var.to_sympy()
     result = sympy.apart(expr_sympy, var_sympy)
     result = from_sympy(result)
     return result

def expand(expression, variable = None):
    expression = expand_(expression)
    if isinstance(expression, Basic) and expression.is_rational_function():
        # Éviter d'appeler `apart()` pour rien, d'autant que FloatFrac() n'est
        # pas compatible avec apart() (sympy v0.7.2).
        num, den = expression.as_numer_denom()
        if den.free_symbols:
            if variable is None:
               variable = extract_var(expression)
            return apart(expression, variable)
    return expression

 def apply(self, expr, var, evaluation):
     'Apart[expr_, var_Symbol]'
     
     expr_sympy = expr.to_sympy()
     var_sympy = var.to_sympy()
     try:
         result = sympy.apart(expr_sympy, var_sympy)
         result = from_sympy(result)
         return result
     except sympy.PolynomialError:
         # raised e.g. for apart(sin(1/(x**2-y**2)))
         return expr

def test_action_verbs():
    assert nsimplify((1/(exp(3*pi*x/5)+1))) == (1/(exp(3*pi*x/5)+1)).nsimplify()
    assert ratsimp(1/x + 1/y) == (1/x + 1/y).ratsimp()
    assert trigsimp(log(x), deep=True) == (log(x)).trigsimp(deep = True)
    assert radsimp(1/(2+sqrt(2))) == (1/(2+sqrt(2))).radsimp()
    assert powsimp(x**y*x**z*y**z, combine='all') == (x**y*x**z*y**z).powsimp(combine='all')
    assert simplify(x**y*x**z*y**z) == (x**y*x**z*y**z).simplify()
    assert together(1/x + 1/y) == (1/x + 1/y).together()
    assert separate((x*(y*z)**3)**2) == ((x*(y*z)**3)**2).separate()
    assert collect(a*x**2 + b*x**2 + a*x - b*x + c, x) == (a*x**2 + b*x**2 + a*x - b*x + c).collect(x)
    assert apart(y/(y+2)/(y+1), y) == (y/(y+2)/(y+1)).apart(y)
    assert combsimp(y/(x+2)/(x+1)) == (y/(x+2)/(x+1)).combsimp()
    assert factor(x**2+5*x+6) == (x**2+5*x+6).factor()
    assert refine(sqrt(x**2)) == sqrt(x**2).refine()
    assert cancel((x**2+5*x+6)/(x+2)) == ((x**2+5*x+6)/(x+2)).cancel()

    def __init__(self):
        self.num_lines, self.num_marks = 7, 3

        self._qp = {}

        self._qp['function_type'] = random.choice(['exp', 'hyperbola', 'log', 'cubic'])

        if self._qp['function_type'] == 'exp':
            func = all_functions.request_exp(difficulty=3)
            self._qp['domain'], self._qp['range'] = func.domain, func.range

            self._qp['equation'] = func.equation
            self._qp['inverse'] = simplify.canonise_log(all_functions.inverse(self._qp['equation']))

        elif self._qp['function_type'] == 'hyperbola':
            func = all_functions.request_hyperbola(difficulty=3)
            self._qp['domain'], self._qp['range'] = func.domain, func.range

            self._qp['equation'] = func.equation
            self._qp['inverse'] = sympy.apart(all_functions.inverse(self._qp['equation']))

        elif self._qp['function_type'] == 'log':
            func = all_functions.request_log(difficulty=3)
            self._qp['domain'], self._qp['range'] = func.domain, func.range

            self._qp['equation'] = func.equation
            self._qp['inverse'] = all_functions.inverse(self._qp['equation'])

        elif self._qp['function_type'] == 'cubic':  # no base class for cubics so I make it here
            m = not_named_yet.randint_no_zero(-3, 3)
            c = not_named_yet.randint_no_zero(-5, 5)
            self._qp['equation'] = m * x ** 3 + c
            self._qp['inverse'] = all_functions.inverse(self._qp['equation'])

            if random.choice([True, False]):
                self._qp['domain'] = sympy.Interval(-sympy.oo, 0)
                if m > 0:
                    self._qp['range'] = sympy.Interval(-sympy.oo, c)
                else:
                    self._qp['range'] = sympy.Interval(c, sympy.oo)
            else:
                self._qp['domain'] = sympy.Interval(0, sympy.oo)
                if m > 0:
                    self._qp['range'] = sympy.Interval(c, sympy.oo)
                else:
                    self._qp['range'] = sympy.Interval(-sympy.oo, c)

        self._qp['inverse_domain'] = self._qp['range']

def given():
    a, b, c, d, e, f = random.sample(range(-9, -1) + range(1, 9), 6)
    #x = a/b
    rs = a - c
    ls = b * x - c
    rs = d * rs / ls
    ls = 1 - e * x / (d + e * x)
    rs = rs / (d + e * x)
    rs = apart(rs)
    ls = ls / f
    rs = rs / f
    ls1, ls2 = ls.as_two_terms()
    rs1, rs2 = rs.as_two_terms()
    g = Struct()
    g.sls = sstr(S(ls1) + S(ls2))
    g.srs = sstr(S(rs1) + S(rs2))
    return g

 def _make_T_denom_sym(self,):
     r'''
     A method to prepare the symbolic expression T_denom_sym for further
     computations. This expression represents the denominator in terms of
     a symbol x, which represents the shortest time delay in the network.
     '''
     self._make_decimal_delays()
     self.x = sp.symbols('x')
     E_sym = sp.Matrix(np.zeros_like(self.M1))
     for i,delay in enumerate(self.Decimal_delays):
         E_sym[i,i] = self.x**int(delay / self.Decimal_gcd)
     M1_sym = sp.Matrix(self.M1)
     self.T_denom_sym = sp.apart((E_sym - M1_sym).det())
     ## I use apart above because sympy yields a function that is not
     ## completely reduced. Alternatively, can use *.as_numer_denom()
     ## and take the first component for the numerator. However, this could
     ## results in spurious roots if the denominator is nontrivial.
     return

def test_H31():
    f = (x**2 + 2*x + 3)/(x**3 + 4*x**2 + 5*x + 2)
    g = 2 / (x + 1)**2 - 2 / (x + 1) + 3 / (x + 2)
    assert apart(f) == g

def apart_fact(Fs):
    S=0
    Ls = apart(Fs).as_ordered_terms()
    for l in Ls:
        S += numer(l)/factor(denom(l))
    return expand_mul(S)

def explanation(pars,for_tex=False):
    (a2,a1,a0),(t1,t2),(k1,k2),(y0,yd0)=pars
    a2,a1,a0 = map(int,[a2,a1,a0])
    Chpol = a2*s**2+a1*s+a0
    r12 = solve(Chpol,s)
    if len(r12)==2:
        r1,r2=r12
    else:
        r1=r12[0]; r2=r1
    # nezavisle od parametrov    
    realroots = (complex(r1).imag == 0)
    kps = t1*(k2-k1)-k2*t2
    continuous = (kps == 0)
    
    if not for_tex:
        beg_eqn = r"$$\begin{eqnarray*}"
        end_eqn = r"\end{eqnarray*}$$"
    else:
        beg_eqn = r"\begin{eqnarray*}"
        end_eqn = r"\end{eqnarray*}"        
    Lls = Chpol * Y(s)
    IC = a1*y0+s*y0+yd0
    ict = IC/Chpol
    y0t = simplify(ILT(ict,s,t))
    if realroots:
        y0t = H(t)*expand_mul(simplify(y0t/H(t)))
    if not continuous:
        V1 = apart(1/s/Chpol).as_ordered_terms()
        Iv1 = ilt_pfe(V1)
        y1t = expand_mul(simplify(Iv1/H(t)))
    V2 = apart(1/s**2/Chpol).as_ordered_terms(); Iv2 = ilt_pfe(V2)
    y2t = expand_mul(simplify(Iv2/H(t)))
    if realroots:
        icpfe=latex(apart_fact(ict))
    else:
        icpfe=latex(ict)
    if not continuous:
        Ykstr=r"""%s
        Y_0(s) &=& %s,\\
        Y_1(s) &=& \frac{1}{s (%s)},\\ 
        Y_2(s) &=& \frac{1}{s^2 (%s)}.\\
        %s""" %(beg_eqn,latex(ict),latex(Chpol),latex(Chpol),end_eqn)
        ytstr=r"""%s
        y_0(t) &=& %s,\\ 
        y_1(t) &=& \theta(t) \left({%s}\right),\\
        y_2(t) &=& \theta(t) \left({%s}\right).\\
        %s""" %(beg_eqn,latex(y0t,fold_func_brackets=True), 
                               latex(y1t,fold_func_brackets=True),
                               latex(y2t,fold_func_brackets=True),end_eqn)
        pfestr=r"""%s
        Y_0(s) &=& %s,\\
        Y_1(s) &=& %s,\\ 
        Y_2(s) &=& %s.\\
        %s""" %(beg_eqn,icpfe,latex(apart(1/s/Chpol)),latex(apart(1/s**2/Chpol)),end_eqn) 
    else:
        Ykstr=r"$$Y_0(s) = %s,\ Y_1(s) = \frac{1}{s^2 (%s)}.$$" %(latex(ict),latex(Chpol))
        ytstr=r"""%s
        y_0(t) &=& %s,\\
        y_1(t) &=& \theta(t) \left({%s}\right).\\
        %s""" %(beg_eqn,latex(y0t,fold_func_brackets=True),
                               latex(y2t,fold_func_brackets=True),end_eqn)
        pfestr=r"$$Y_0(s) =%s,\ Y_1(s) = %s.$$" %(icpfe,latex(apart_fact(1/s**2/Chpol)))
    ic1,ic2 = s*Y(s)- y0, s**2*Y(s) - y0*s - yd0
    f = k1*t*(H(t)-H(t-t1))+k2*(t-t2)*(H(t-t1)-H(t-t2))
    ct1,ct2=factor((k2-k1)*(t-t1)),factor(-k2*(t-t2))
    fh =k1*t*H(t)+ct1*H(t-t1)+ct2*H(t-t2)
    if not continuous:
        fh = fh + kps*H(t-t1)
    Lfh = expand_mul(LT(fh,t,s)[0])
    y2t1 = y2t.subs(t,t-t1); y2t2 = y2t.subs(t,t-t2)
    if not continuous:
        y1t1 = y1t.subs(t,t-t1)
    yt =  y0t + k1*y2t + (k2-k1)*y2t1*H(t-t1) - k2*y2t2*H(t-t2)
    if not continuous:
        yt = yt + kps*y1t1*H(t-t1)
    yder = diff(y0t/H(t),t)*H(t) + k1*diff(y2t,t)*H(t)+(k2-k1)*diff(y2t1,t)*H(t-t1)-k2*diff(y2t2,t)*H(t-t2)
    if not continuous:
        yder = yder + kps*diff(y1t1,t)*H(t-t1)
        
    yder2 = diff(y0t/H(t),t,2)*H(t) + k1*diff(y2t,t,2)*H(t)+(k2-k1)*diff(y2t1,t,2)*H(t-t1)-k2*diff(y2t2,t,2)*H(t-t2)
    if not continuous:
        yder2 = yder2 + kps*diff(y1t1,t,2)*H(t-t1)
    if continuous:    
        sol2 = latex(k1 * y1(t) + (k2-k1)*y1(t-t1) - k2 * y1(t-t2))
    else:
        sol2 = latex(k1 * y2(t) + (k2-k1)*y2(t-t1) - k2 * y2(t-t2) + kps*y1(t-t1))
    sol2 = sol2.replace('\\operatorname','')
    if kps < 0:
        SolStr = "y(t) = y_0(t) " + sol2
    else:    
        SolStr = "y(t) = y_0(t) + " + sol2
        
    Fcheck=a2*yder2+a1*yder+a0*yt
    Is_solution=simplify(Fcheck-f)
    Estr=[r"""Nech $y(t) \rightarrow Y(s)$. Pre transformáciu ľavej strany použijeme vzorce:
    %s
    y'(t)  &\rightarrow & s Y(s)- y(0)=%s,\\ 
    y''(t) &\rightarrow & s^2 Y(s) - s y(0) -y'(0)= %s.\\
    %s""" %(beg_eqn,latex(ic1), latex(ic2),end_eqn),
    r"""Pravú stranu prepíšeme pomocou Heavisideovho skoku $\theta(t)\,$ takto:
#.........这里部分代码省略.........

from sympy import together, apart, symbols
from sympy import *

x1, x2, x3, x4 = symbols('x1 x2 x3 x4')
x1/x2 + x3/x4
together(x1/x2 + x3/x4)


apart ((2*x**2+3*x+4)/(x+1))
together(apart ((2*x**2+3*x+4)/(x+1)))

log(4).n()
log(4,4).n()
ln(4).n()
mpmath.log10(4)

from __future__ import division, print_function

from sympy import S, apart, summation


k = S('k')
k_domain = (k, 1, 00)

f_Y = 18/((k+1)*(k+2)*(k+3))


print(
    " = ".join([
        f_Y,
        apart(f_Y)
    ])
)

EY = summation(
    f_Y,
    k_domain
)

print(
    "E[Y] = {}".format(
        EY
    )
)

    if isinstance(e, (sympy.Add, sympy.Mul)):
        t = type(e)
        return t(*[laplace_transform_derivatives(arg) for arg in e.args])

    return e

L_ode_2 = laplace_transform_derivatives(L_ode)
print(L_ode_2)
L_ode_3 = L_ode_2.subs(L_y, Y)
print(L_ode_3)
ics = {y(0): 1, y(t).diff(t).subs(t, 0): 0}
print(ics)
L_ode_4 = L_ode_3.subs(ics)
print(L_ode_4)
Y_sol = sympy.solve(L_ode_4, Y)
print(Y_sol)
sympy.apart(Y_sol[0])
y_sol = sympy.inverse_laplace_transform(Y_sol[0], s, t)
sympy.simplify(y_sol)
y_t = sympy.lambdify(t, y_sol, 'numpy')

fig, ax = plt.subplots(figsize=(8, 4))

tt = np.linspace(0, 10, 500)
ax.plot(tt, y_t(tt).real)
ax.set_xlabel(r"$t$", fontsize=18)
ax.set_ylabel(r"$y(t)$", fontsize=18)
fig.tight_layout()

os.system("pause")