<SYSTEM_TASK:>
Solve the following problem using Python, implementing the functions described below, one line at a time
<END_TASK>
<USER_TASK:>
Description:
 def conjugate_gradient_normal(op, x, rhs, niter=1, callback=None):
 """Optimized implementation of CG for the normal equation. This method solves the inverse problem (of the first kind) :: A(x) == rhs with a linear `Operator` ``A`` by looking at the normal equation :: A.adjoint(A(x)) == A.adjoint(rhs) It uses a minimum amount of memory copies by applying re-usable temporaries and in-place evaluation. The method is described (for linear systems) in a `Wikipedia article <https://en.wikipedia.org/wiki/Conjugate_gradient_method#\ Conjugate_gradient_on_the_normal_equations>`_. Parameters op : `Operator` Operator in the inverse problem. If not linear, it must have an implementation of `Operator.derivative`, which in turn must implement `Operator.adjoint`, i.e. the call ``op.derivative(x).adjoint`` must be valid. x : ``op.domain`` element Element to which the result is written. Its initial value is used as starting point of the iteration, and its values are updated in each iteration step. rhs : ``op.range`` element Right-hand side of the equation defining the inverse problem niter : int Number of iterations. callback : callable, optional Object executing code per iteration, e.g. plotting each iterate. See Also -------- conjugate_gradient : Optimized solver for symmetric matrices odl.solvers.smooth.nonlinear_cg.conjugate_gradient_nonlinear : Equivalent solver for the nonlinear case """

    # TODO: add a book reference
    # TODO: update doc

    if x not in op.domain:
        raise TypeError('`x` {!r} is not in the domain of `op` {!r}'
                        ''.format(x, op.domain))

    d = op(x)
    d.lincomb(1, rhs, -1, d)               # d = rhs - A x
    p = op.derivative(x).adjoint(d)
    s = p.copy()
    q = op.range.element()
    sqnorm_s_old = s.norm() ** 2  # Only recalculate norm after update

    for _ in range(niter):
        op(p, out=q)                       # q = A p
        sqnorm_q = q.norm() ** 2
        if sqnorm_q == 0.0:  # Return if residual is 0
            return

        a = sqnorm_s_old / sqnorm_q
        x.lincomb(1, x, a, p)               # x = x + a*p
        d.lincomb(1, d, -a, q)              # d = d - a*Ap
        op.derivative(p).adjoint(d, out=s)  # s = A^T d

        sqnorm_s_new = s.norm() ** 2
        b = sqnorm_s_new / sqnorm_s_old
        sqnorm_s_old = sqnorm_s_new

        p.lincomb(1, s, b, p)               # p = s + b * p

        if callback is not None:
            callback(x)