
    Џkh                     p   d Z ddlZddlZddlmZ ddlmZ ddlm	Z	m
Z
mZmZmZmZ ddlmZmZ dd	lmZ dd
lmZ ddlmZ d Zej4                  d        Zd Zd ZddZddddej>                   ej>                  fddddddfdZ d Z!d Z"d Z# G d d      Z$ej>                   ej>                  fddfdZ%y)z'Routines for numerical differentiation.    N)norm)LinearOperator   )issparse
isspmatrixfind	csc_array	csr_array
csr_matrix   )group_densegroup_sparse)array_namespace)
MapWrapper)array_api_extrac                    |dk(  rt        j                  |t              }nA|dk(  r1t        j                  |      }t        j                  |t              }nt        d      t        j                  |t         j                   k(  |t         j                  k(  z        r||fS ||z  }|j                         }| |z
  }	|| z
  }
|dk(  ry| |z   }||k  ||kD  z  }t        j                  |      t        j                  |	|
      k  }|||z  xx   dz  cc<   |
|	k\  | z  }|
|   |z  ||<   |
|	k  | z  }|	|    |z  ||<   ||fS |dk(  r|	|k\  |
|k\  z  }|
|	k\  | z  }t        j                  ||   d|
|   z  |z        ||<   d||<   |
|	k  | z  }t        j                  ||   d|	|   z  |z         ||<   d||<   t        j                  |
|	      |z  }| t        j                  |      |k  z  }||   ||<   d||<   ||fS )	a  Adjust final difference scheme to the presence of bounds.

    Parameters
    ----------
    x0 : ndarray, shape (n,)
        Point at which we wish to estimate derivative.
    h : ndarray, shape (n,)
        Desired absolute finite difference steps.
    num_steps : int
        Number of `h` steps in one direction required to implement finite
        difference scheme. For example, 2 means that we need to evaluate
        f(x0 + 2 * h) or f(x0 - 2 * h)
    scheme : {'1-sided', '2-sided'}
        Whether steps in one or both directions are required. In other
        words '1-sided' applies to forward and backward schemes, '2-sided'
        applies to center schemes.
    lb : ndarray, shape (n,)
        Lower bounds on independent variables.
    ub : ndarray, shape (n,)
        Upper bounds on independent variables.

    Returns
    -------
    h_adjusted : ndarray, shape (n,)
        Adjusted absolute step sizes. Step size decreases only if a sign flip
        or switching to one-sided scheme doesn't allow to take a full step.
    use_one_sided : ndarray of bool, shape (n,)
        Whether to switch to one-sided scheme. Informative only for
        ``scheme='2-sided'``.
    1-sideddtype2-sidedz(`scheme` must be '1-sided' or '2-sided'.      ?TF)np	ones_likeboolabs
zeros_like
ValueErrorallinfcopymaximumminimum)x0h	num_stepsschemelbubuse_one_sidedh_total
h_adjusted
lower_dist
upper_distxviolatedfittingforwardbackwardcentralmin_distadjusted_centrals                      S/var/www/teggl/fontify/venv/lib/python3.12/site-packages/scipy/optimize/_numdiff.py_adjust_scheme_to_boundsr8      sX   > Qd3	9	FF1Iat4CDD	vvrbffW}rvv./-)mGJbJbJLFq2v&&&/RZZ
J%GG8g%&",&+x7(1I=
7+x7 *8 44y@
8& }$$% 
9	(Z7-BC+x7 jjgJj11I=?
7!%g+x7 "

hKz(33i?!A  A
8"&h::j*5	A$Hz(:h(FG'/0@'A
#$*/&'}$$    c                 4   t        j                  t         j                        j                  }d}t        j                  | t         j
                        r@t        j                  |       j                  }t        j                  |       j                  }d}t        j                  |t         j
                        rEt        j                  |      j                  }|r$|k  rt        j                  |      j                  }|dv r|dz  S |dv r|dz  S t        d      )a  
    Calculates relative EPS step to use for a given data type
    and numdiff step method.

    Progressively smaller steps are used for larger floating point types.

    Parameters
    ----------
    f0_dtype: np.dtype
        dtype of function evaluation

    x0_dtype: np.dtype
        dtype of parameter vector

    method: {'2-point', '3-point', 'cs'}

    Returns
    -------
    EPS: float
        relative step size. May be np.float16, np.float32, np.float64

    Notes
    -----
    The default relative step will be np.float64. However, if x0 or f0 are
    smaller floating point types (np.float16, np.float32), then the smallest
    floating point type is chosen.
    FT)2-pointcsr   )3-pointgUUUUUU?zBUnknown step method, should be one of {'2-point', '3-point', 'cs'})	r   finfofloat64eps
issubdtypeinexactr   itemsizeRuntimeError)x0_dtypef0_dtypemethodEPSx0_is_fpx0_itemsizef0_itemsizes          r7   _eps_for_methodrL   ]   s    < ((2::

"
"CH	}}Xrzz*hhx $$hhx(11	}}Xrzz*hhx(11k1((8$((C""Cx	;	Sz : ; 	;r9   c           
         |dk\  j                  t              dz  dz
  }t        |j                  |j                  |      }| 1||z  t	        j
                  dt	        j                  |            z  }|S | |z  t	        j                  |      z  }||z   |z
  }t	        j                  |dk(  ||z  t	        j
                  dt	        j                  |            z  |      }|S )az  
    Computes an absolute step from a relative step for finite difference
    calculation.

    Parameters
    ----------
    rel_step: None or array-like
        Relative step for the finite difference calculation
    x0 : np.ndarray
        Parameter vector
    f0 : np.ndarray or scalar
    method : {'2-point', '3-point', 'cs'}

    Returns
    -------
    h : float
        The absolute step size

    Notes
    -----
    `h` will always be np.float64. However, if `x0` or `f0` are
    smaller floating point dtypes (e.g. np.float32), then the absolute
    step size will be calculated from the smallest floating point size.
    r   r   r         ?)astypefloatrL   r   r   r"   r   where)rel_stepr$   f0rG   sign_x0rstepabs_stepdxs           r7   _compute_absolute_steprX      s    6 Qwu%)A-GBHHbhh7E7?RZZRVVBZ%@@ O g%r
2 H}"88B!G!GObjjbffRj.II$& Or9   c                     d | D        \  }}|j                   dk(  r t        j                  ||j                        }|j                   dk(  r t        j                  ||j                        }||fS )aa  
    Prepares new-style bounds from a two-tuple specifying the lower and upper
    limits for values in x0. If a value is not bound then the lower/upper bound
    will be expected to be -np.inf/np.inf.

    Examples
    --------
    >>> _prepare_bounds([(0, 1, 2), (1, 2, np.inf)], [0.5, 1.5, 2.5])
    (array([0., 1., 2.]), array([ 1.,  2., inf]))
    c              3   R   K   | ]  }t        j                  |t                ! yw)r   N)r   asarrayrP   ).0bs     r7   	<genexpr>z"_prepare_bounds.<locals>.<genexpr>   s     9Qbjj%((9s   %'r   )ndimr   resizeshape)boundsr$   r(   r)   s       r7   _prepare_boundsrc      sY     :&9FB	ww!|YYr288$	ww!|YYr288$r6Mr9   c                    t        |       rt        |       } n7t        j                  |       } | dk7  j	                  t        j
                        } | j                  dk7  rt        d      | j                  \  }}|t        j                  |      r1t        j                  j                  |      }|j                  |      }n0t        j                  |      }|j                  |fk7  rt        d      | dd|f   } t        |       r#t        ||| j                  | j                         }nt#        |||       }|j%                         ||<   |S )a  Group columns of a 2-D matrix for sparse finite differencing [1]_.

    Two columns are in the same group if in each row at least one of them
    has zero. A greedy sequential algorithm is used to construct groups.

    Parameters
    ----------
    A : array_like or sparse array, shape (m, n)
        Matrix of which to group columns.
    order : int, iterable of int with shape (n,) or None
        Permutation array which defines the order of columns enumeration.
        If int or None, a random permutation is used with `order` used as
        a random seed. Default is 0, that is use a random permutation but
        guarantee repeatability.

    Returns
    -------
    groups : ndarray of int, shape (n,)
        Contains values from 0 to n_groups-1, where n_groups is the number
        of found groups. Each value ``groups[i]`` is an index of a group to
        which ith column assigned. The procedure was helpful only if
        n_groups is significantly less than n.

    References
    ----------
    .. [1] A. Curtis, M. J. D. Powell, and J. Reid, "On the estimation of
           sparse Jacobian matrices", Journal of the Institute of Mathematics
           and its Applications, 13 (1974), pp. 117-120.
    r   r   z`A` must be 2-dimensional.Nz`order` has incorrect shape.)r   r	   r   
atleast_2drO   int32r_   r   ra   isscalarrandomRandomStatepermutationr[   r   indicesindptrr   r!   )Aordermnrnggroupss         r7   group_columnsrs      s   < {aLMM!!VOOBHH%vv{56677DAq}E*ii##E*"

5!;;1$;<<	!U(A{aAIIqxx8Q1%KKMF5MMr9   r=   F c                 t   |dvrt        d| d      ddi}t        |      }t        j                  |j	                  |      d|      }|j
                  }|j                  |j                  d      r|j                  }|j                  ||      }|j                  dkD  rt        d	      t        ||      \  }}|j                  |j                  k7  s|j                  |j                  k7  rt        d
      |r[t        j                  t        j                  |            r(t        j                  t        j                  |            st        d      |
i }
t        | ||	|
      }dx}}| ||      }d}n/t        j                   |      }|j                  dkD  rt        d      t        j"                  ||k  ||kD  z        rt        d      |r7|!t%        |j                  |j                  |      }t'        |||||      \  }}n|t)        ||||      }n|dk\  j                  t*              dz  dz
  }|}||z   |z
  }t        j,                  |dk(  t%        |j                  |j                  |      |z  t        j.                  dt        j0                  |            z  |      }|dk(  rt3        ||dd||      \  }}n |dk(  rt3        ||dd||      \  }}n|dk(  rd}|xs t4        }t7        |      5 }|t9        ||||||      \  }}nt;        |      st=        |      dk(  r|\  }}n|}t?        |      }t;        |      r|jA                         }nt        jB                  |      }t        j                   |      }tE        ||||||||	      \  }}ddd       |r||z  }||d<   |fS S # 1 sw Y   xY w)uZ*  Compute finite difference approximation of the derivatives of a
    vector-valued function.

    If a function maps from R^n to R^m, its derivatives form m-by-n matrix
    called the Jacobian, where an element (i, j) is a partial derivative of
    f[i] with respect to x[j].

    Parameters
    ----------
    fun : callable
        Function of which to estimate the derivatives. The argument x
        passed to this function is ndarray of shape (n,) (never a scalar
        even if n=1). It must return 1-D array_like of shape (m,) or a scalar.
    x0 : array_like of shape (n,) or float
        Point at which to estimate the derivatives. Float will be converted
        to a 1-D array.
    method : {'3-point', '2-point', 'cs'}, optional
        Finite difference method to use:
            - '2-point' - use the first order accuracy forward or backward
                          difference.
            - '3-point' - use central difference in interior points and the
                          second order accuracy forward or backward difference
                          near the boundary.
            - 'cs' - use a complex-step finite difference scheme. This assumes
                     that the user function is real-valued and can be
                     analytically continued to the complex plane. Otherwise,
                     produces bogus results.
    rel_step : None or array_like, optional
        Relative step size to use. If None (default) the absolute step size is
        computed as ``h = rel_step * sign(x0) * max(1, abs(x0))``, with
        `rel_step` being selected automatically, see Notes. Otherwise
        ``h = rel_step * sign(x0) * abs(x0)``. For ``method='3-point'`` the
        sign of `h` is ignored. The calculated step size is possibly adjusted
        to fit into the bounds.
    abs_step : array_like, optional
        Absolute step size to use, possibly adjusted to fit into the bounds.
        For ``method='3-point'`` the sign of `abs_step` is ignored. By default
        relative steps are used, only if ``abs_step is not None`` are absolute
        steps used.
    f0 : None or array_like, optional
        If not None it is assumed to be equal to ``fun(x0)``, in this case
        the ``fun(x0)`` is not called. Default is None.
    bounds : tuple of array_like, optional
        Lower and upper bounds on independent variables. Defaults to no bounds.
        Each bound must match the size of `x0` or be a scalar, in the latter
        case the bound will be the same for all variables. Use it to limit the
        range of function evaluation. Bounds checking is not implemented
        when `as_linear_operator` is True.
    sparsity : {None, array_like, sparse array, 2-tuple}, optional
        Defines a sparsity structure of the Jacobian matrix. If the Jacobian
        matrix is known to have only few non-zero elements in each row, then
        it's possible to estimate its several columns by a single function
        evaluation [3]_. To perform such economic computations two ingredients
        are required:

        * structure : array_like or sparse array of shape (m, n). A zero
          element means that a corresponding element of the Jacobian
          identically equals to zero.
        * groups : array_like of shape (n,). A column grouping for a given
          sparsity structure, use `group_columns` to obtain it.

        A single array or a sparse array is interpreted as a sparsity
        structure, and groups are computed inside the function. A tuple is
        interpreted as (structure, groups). If None (default), a standard
        dense differencing will be used.

        Note, that sparse differencing makes sense only for large Jacobian
        matrices where each row contains few non-zero elements.
    as_linear_operator : bool, optional
        When True the function returns an `scipy.sparse.linalg.LinearOperator`.
        Otherwise it returns a dense array or a sparse array depending on
        `sparsity`. The linear operator provides an efficient way of computing
        ``J.dot(p)`` for any vector ``p`` of shape (n,), but does not allow
        direct access to individual elements of the matrix. By default
        `as_linear_operator` is False.
    args, kwargs : tuple and dict, optional
        Additional arguments passed to `fun`. Both empty by default.
        The calling signature is ``fun(x, *args, **kwargs)``.
    full_output : bool, optional
        If True then the function also returns a dictionary with extra information
        about the calculation.
    workers : int or map-like callable, optional
        Supply a map-like callable, such as
        `multiprocessing.Pool.map` for evaluating the population in parallel.
        This evaluation is carried out as ``workers(fun, iterable)``.
        Alternatively, if `workers` is an int the task is subdivided into `workers`
        sections and the fun evaluated in parallel
        (uses `multiprocessing.Pool <multiprocessing>`).
        Supply -1 to use all available CPU cores.
        It is recommended that a map-like be used instead of int, as repeated
        calls to `approx_derivative` will incur large overhead from setting up
        new processes.

    Returns
    -------
    J : {ndarray, sparse array, LinearOperator}
        Finite difference approximation of the Jacobian matrix.
        If `as_linear_operator` is True returns a LinearOperator
        with shape (m, n). Otherwise it returns a dense array or sparse
        array depending on how `sparsity` is defined. If `sparsity`
        is None then a ndarray with shape (m, n) is returned. If
        `sparsity` is not None returns a csr_array or csr_matrix with
        shape (m, n) following the array/matrix type of the incoming structure.
        For sparse arrays and linear operators it is always returned as
        a 2-D structure. For ndarrays, if m=1 it is returned
        as a 1-D gradient array with shape (n,).

    info_dict : dict
        Dictionary containing extra information about the calculation. The
        keys include:

        - `nfev`, number of function evaluations. If `as_linear_operator` is True
           then `fun` is expected to track the number of evaluations itself.
           This is because multiple calls may be made to the linear operator which
           are not trackable here.

    See Also
    --------
    check_derivative : Check correctness of a function computing derivatives.

    Notes
    -----
    If `rel_step` is not provided, it assigned as ``EPS**(1/s)``, where EPS is
    determined from the smallest floating point dtype of `x0` or `fun(x0)`,
    ``np.finfo(x0.dtype).eps``, s=2 for '2-point' method and
    s=3 for '3-point' method. Such relative step approximately minimizes a sum
    of truncation and round-off errors, see [1]_. Relative steps are used by
    default. However, absolute steps are used when ``abs_step is not None``.
    If any of the absolute or relative steps produces an indistinguishable
    difference from the original `x0`, ``(x0 + dx) - x0 == 0``, then a
    automatic step size is substituted for that particular entry.

    A finite difference scheme for '3-point' method is selected automatically.
    The well-known central difference scheme is used for points sufficiently
    far from the boundary, and 3-point forward or backward scheme is used for
    points near the boundary. Both schemes have the second-order accuracy in
    terms of Taylor expansion. Refer to [2]_ for the formulas of 3-point
    forward and backward difference schemes.

    For dense differencing when m=1 Jacobian is returned with a shape (n,),
    on the other hand when n=1 Jacobian is returned with a shape (m, 1).
    Our motivation is the following: a) It handles a case of gradient
    computation (m=1) in a conventional way. b) It clearly separates these two
    different cases. b) In all cases np.atleast_2d can be called to get 2-D
    Jacobian with correct dimensions.

    References
    ----------
    .. [1] W. H. Press et. al. "Numerical Recipes. The Art of Scientific
           Computing. 3rd edition", sec. 5.7.

    .. [2] A. Curtis, M. J. D. Powell, and J. Reid, "On the estimation of
           sparse Jacobian matrices", Journal of the Institute of Mathematics
           and its Applications, 13 (1974), pp. 117-120.

    .. [3] B. Fornberg, "Generation of Finite Difference Formulas on
           Arbitrarily Spaced Grids", Mathematics of Computation 51, 1988.

    Examples
    --------
    >>> import numpy as np
    >>> from scipy.optimize._numdiff import approx_derivative
    >>>
    >>> def f(x, c1, c2):
    ...     return np.array([x[0] * np.sin(c1 * x[1]),
    ...                      x[0] * np.cos(c2 * x[1])])
    ...
    >>> x0 = np.array([1.0, 0.5 * np.pi])
    >>> approx_derivative(f, x0, args=(1, 2))
    array([[ 1.,  0.],
           [-1.,  0.]])

    Bounds can be used to limit the region of function evaluation.
    In the example below we compute left and right derivative at point 1.0.

    >>> def g(x):
    ...     return x**2 if x >= 1 else x
    ...
    >>> x0 = 1.0
    >>> approx_derivative(g, x0, bounds=(-np.inf, 1.0))
    array([ 1.])
    >>> approx_derivative(g, x0, bounds=(1.0, np.inf))
    array([ 2.])

    We can also parallelize the derivative calculation using the workers
    keyword.

    >>> from multiprocessing import Pool
    >>> import time
    >>> def fun2(x):       # import from an external file for use with multiprocessing
    ...     time.sleep(0.002)
    ...     return rosen(x)

    >>> rng = np.random.default_rng()
    >>> x0 = rng.uniform(high=10, size=(2000,))
    >>> f0 = rosen(x0)

    >>> %timeit approx_derivative(fun2, x0, f0=f0)     # may vary
    10.5 s ± 5.91 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

    >>> elapsed = []
    >>> with Pool() as workers:
    ...     for i in range(10):
    ...         t = time.perf_counter()
    ...         approx_derivative(fun2, x0, workers=workers.map, f0=f0)
    ...         et = time.perf_counter()
    ...         elapsed.append(et - t)
    >>> np.mean(elapsed)    # may vary
    np.float64(1.442545195999901)

    Create a map-like vectorized version. `x` is a generator, so first of all
    a 2-D array, `xx`, is reconstituted. Here `xx` has shape `(Y, N)` where `Y`
    is the number of function evaluations to perform and `N` is the dimensionality
    of the objective function. The underlying objective function is `rosen`, which
    requires `xx` to have shape `(N, Y)`, so a transpose is required.

    >>> def fun(f, x, *args, **kwds):
    ...     xx = np.r_[[xs for xs in x]]
    ...     return f(xx.T)
    >>> %timeit approx_derivative(fun2, x0, workers=fun, f0=f0)    # may vary
    91.8 ms ± 755 μs per loop (mean ± std. dev. of 7 runs, 10 loops each)

    )r;   r=   r<   zUnknown method 'z'. nfevNr   )r_   xpreal floatingz#`x0` must have at most 1 dimension.z,Inconsistent shapes between bounds and `x0`.z7Bounds not supported when `as_linear_operator` is True.r   z&`f0` passed has more than 1 dimension.z `x0` violates bound constraints.r   rN   r;   r   r=   r   r<   F)#r   r   xpx
atleast_ndr[   r?   isdtyper   rO   r_   rc   ra   r   r   isinf_Fun_Wrapper
atleast_1danyrL   _linear_operator_differencerX   rP   rQ   r"   r   r8   mapr   _dense_differencer   lenrs   tocscre   _sparse_difference)funr$   rG   rR   rV   rS   rb   sparsityas_linear_operatorargskwargsfull_outputworkers	info_dictrw   _x_dtyper(   r)   fun_wrappedrv   _nfevJ_r%   rT   rW   r*   mf	structurerr   s                                  r7   approx_derivativer     s   F 11+F83788I		B	

2Q2	6BZZF	zz"((O, 
2v	B	ww{>??VR(FB	xx288rxx2883GHH266"((2,#7')vvbhhrl'; 9 : 	: ~sBf5K D5	z_]]277Q;EFF	vvrBw27#$;<<&rxx6BH*;+-xA1 &xR@A Qw&&u-1A5GA 6R-Bq(288VD !#%::c266":#>?A
 Y7Aq)R -A}y 7Aq)R -A}t^!M .S  	AB,["b!)6)+-5  )c(mq.@(0%Iv (I*84FI& ) 1I "i 8Iv.-k2r1-:I-3VRA5#	A*  	&)|5	A 	As   4B N..N7c                      j                   j                   }|dk(  r
 fd}n'|dk(  r	 fd}n|dk(  r	 fd}nt        d      t        |f|      dfS )	Nr;   c                     t        j                  | t        j                  |             rt        j                        S t	        |       z  }|| z  z   } |      z
  }||z  S Nr   array_equalr   zerosr   )	prW   r/   dfrS   r   r%   ro   r$   s	       r7   matvecz+_linear_operator_difference.<locals>.matvecr  sX    ~~aq!12xx{"T!WBRT	AQ"B7Nr9   r=   c                    t        j                  | t        j                  |             rt        j                  	      S dz  t	        |       z  }
|dz  | z  z
  }
|dz  | z  z   } |      } |      }||z
  }||z  S Nr   r   )r   rW   x1x2f1f2r   r   r%   ro   r$   s          r7   r   z+_linear_operator_difference.<locals>.matvec|  s    ~~aq!12xx{"1tAwBr!tQhBr!tQhBRBRBbB7Nr9   r<   c                     t        j                  | t        j                  |             rt        j                        S t	        |       z  }|| z  dz  z   } |      }|j
                  }||z  S N              ?)r   r   r   r   r   imag)	r   rW   r/   r   r   r   r%   ro   r$   s	        r7   r   z+_linear_operator_difference.<locals>.matvec  sa    ~~aq!12xx{"T!WBRT#XAQBB7Nr9   Never be here.r   )sizerD   r   )r   r$   rS   r%   rG   rp   r   ro   s   ````   @r7   r   r   l  sr    
A
A	 	 
9			 		 
4	 	 +,,1a&&)1,,r9   c           
      ,   |j                   }|j                   t        j                  |f      }d}	|dk(  rfd}
 ||  |
||            }t              D cg c]  }||   ||   z   ||   z
   }}|D cg c]  }||z
  	 }}t	        ||      D cg c]
  \  }}||z   }}}|	t        |      z  }	n|dk(  rd }t         ||  ||||                  } ||||      }t               }t               }t        |      D ]  \  }}t        |      }t        |      }t        |      }t        |      }|r8|j                  ||   ||   z
         |j                  d|z  d|z  z   |z
         l|j                  ||   ||   z
         |j                  ||z
          t	        ||      D cg c]
  \  }}||z   }}}|	dt        |      z  z  }	nh|d	k(  rXfd
}t         ||  |||                  }t	        ||      D cg c]  \  }}|j                  |z   }}}|	t        |      z  }	nt        d      t        |      D ]
  \  }}|||<    |dk(  rt        j                  |      }|j                  |	fS c c}w c c}w c c}}w c c}}w c c}}w )Nr   r;   c              3   |   K   t              D ])  }t        j                  |       }| |   ||   z   ||<   | + y wr   )ranger   r!   )r$   r%   ir   rp   s       r7   x_generator2z'_dense_difference.<locals>.x_generator2  sC     1X 	 WWR[1!1	s   9<r=   c              3     K   t        |      D ]u  \  }}t        j                  |       }t        j                  |       }|r | |   ||   z   ||<   | |   d||   z  z   ||<   n| |   ||   z
  ||<   | |   ||   z   ||<   | | w y wr   )	enumerater   r!   )r$   r%   r*   r   	one_sidedr   r   s          r7   x_generator3z'_dense_difference.<locals>.x_generator3  s      )- 8 
9WWR[WWR[qEAaDLBqEqEAadFNBqEqEAaDLBqEqEAaDLBqE
s   BBg         r   r<   c              3      K   t              D ]0  }| j                  t        d      }||xx   ||   dz  z  cc<   | 2 y w)NTr!   r   )r   rO   complex)r$   r%   r   xcrp   s       r7   x_generator_csz)_dense_difference.<locals>.x_generator_cs  sG     1X YYwTY211#s   A Ar   r   )r   r   emptyr   zipr   iterlistr   nextappendr   rD   ravelT)r   r$   rS   r%   r*   rG   r   ro   J_transposedrv   r   f_evalsr   rW   f_evalr   delfdelxdf_dxr   genr   lur   r   r   hivrp   s                                @r7   r   r     s   
A
A88QF#LD
	 #|B23.3Ah7r!uqt|r!u$77(/0ffrk00/22r{;t;;E
	9		 wsLQ$FGH2q-0VV%m4 	#LAyS	AS	AgBgB		!A$A,'		$)a"f,r12		!A$1+&		"r'"	# 032r{;t;;CJ	4	 wsN2q$9:;,/O<&"b2<<E
+,,%  1Q 	Avxx->>4u 80;F < =s   I:8I?J.J
Jc	                   #$ |j                   }	j                   }
g }g }g }t        j                        dz   $d}$fd##fd}#fd}#fd}|dk(  rt         ||  |                   } |       }n@|dk(  rt         ||  |                   } |       }n|d	k(  rt         ||  |                   } #       D ]  }t        j                  |      \  }t        |d d |f         \  }}}||   }|dk(  r"t              z
  }t              |z
  }|dz  }n|dk(  rt              }t        |      }|z  } |z  }t        j                  |
      }||   |   z
  ||<   ||   ||   z
  ||<   t              }t        |      } |d
z  }|   }!t        j                  |	      }||!   }"d||"   z  d||"   z  z   | |"   z
  ||"<   ||!    }"| |"   ||"   z
  ||"<   n2|d	k(  r"t              }|dz  }|j                  }|z  }nt        d      |j                  |       |j                  |       |j                  ||   ||   z          t        j                  |      }t        j                  |      }t        j                  |      }t        |      rt        |||ff|	|
f      |fS t        |||ff|	|
f      |fS )Nr   r   c               3   ^   K   t              D ]  } t        j                  |         y wr   )r   r   equal)grouprr   n_groupss    r7   e_generatorz'_sparse_difference.<locals>.e_generator  s+     8_ 	*E((5&))	*s   *-c               3   F   K           } | D ]  }|z  }|z   }|  y wr   rt   )e_geneh_vecr/   r   r%   r$   s       r7   r   z(_sparse_difference.<locals>.x_generator2  s4      	AEEU
AG	s   !c               3      K           } | D ]}  }|z  }
j                         }
j                         }	|z  }||xx   ||   z  cc<   ||xx   d||   z  z  cc<   	 |z  }||xx   ||   z  cc<   ||xx   ||   z  cc<   | |  y wr   r   )r   r   r   r   r   mask_1mask_2r   r%   r*   r$   s          r7   r   z(_sparse_difference.<locals>.x_generator3  s      	AEEBB"Q&FvJ%-'JvJ!eFm++J#^a'FvJ%-'JvJ%-'JHH	s   BBc               3   H   K           } | D ]  }|z  }|dz  z     y wr   rt   )r   r   r   r   r%   r$   s      r7   r   z*_sparse_difference.<locals>.x_generator_cs  s5      	#AEEus{""	#s   "r;   r=   r<   r   r   r   )ra   )r   r   maxr   nonzeror   r   r   r   r   r   r   hstackr   r   r
   )%r   r$   rS   r%   r*   r   rr   rG   r   ro   rp   row_indicescol_indices	fractionsrv   r   r   r   r   xsr   colsr   jr   rW   r   r   r   r   r   r   r   maskrowsr   r   s%    ` `` `                            @@r7   r   r     s   
A
AKKIvvf~!HD*
"# wsLN34^	9	wsLN34^	4wsN$456] 2( 

1yD)*1aGYbBBg#BAIDy  bBbB"Q&F#^a'F!BFbj0BvJFbj0BvJgBgBAID #D!BT7DBtH}q2d8|3bh>BtHdU8D$x"T(*BtHt^gBAIDBQB-.. 	11AA'e2(h ))K(K))K(K		)$I)9{K&@A!QPRVVVi+{!;<QFKTQQr9   c                       e Zd Zd Zd Zy)r}   c                 <    || _         || _        || _        || _        y r   )r   r$   r   r   )selfr   r$   r   r   s        r7   __init__z_Fun_Wrapper.__init__a  s    	r9   c                 b   t        | j                        }|j                  |j                  d      r&|j	                  || j                  j                        }t        j                   | j                  |g| j                  i | j                        }|j                  dkD  rt        d      |S )Nrx   r   z-`fun` return value has more than 1 dimension.)r   r$   r{   r   rO   r   r~   r   r   r   r_   rD   )r   r/   rw   fs       r7   __call__z_Fun_Wrapper.__call__g  s     TWW%::agg/		!TWW]]+AMM($((1@tyy@DKK@A66A:  8 9 9r9   N)__name__
__module____qualname__r   r   rt   r9   r7   r}   r}   _  s    r9   r}   c           	      Z   |i } ||g|i |}t        |      rt        | |||||      }t        |      }||z
  }t        |      \  }	}
}t	        j
                  ||	|
f         j                         }t	        j                  t	        j                  |      t	        j                  dt	        j                  |            z        S t        | ||||      }t	        j                  ||z
        }t	        j                  |t	        j                  dt	        j                  |            z        S )aS	  Check correctness of a function computing derivatives (Jacobian or
    gradient) by comparison with a finite difference approximation.

    Parameters
    ----------
    fun : callable
        Function of which to estimate the derivatives. The argument x
        passed to this function is ndarray of shape (n,) (never a scalar
        even if n=1). It must return 1-D array_like of shape (m,) or a scalar.
    jac : callable
        Function which computes Jacobian matrix of `fun`. It must work with
        argument x the same way as `fun`. The return value must be array_like
        or sparse array with an appropriate shape.
    x0 : array_like of shape (n,) or float
        Point at which to estimate the derivatives. Float will be converted
        to 1-D array.
    bounds : 2-tuple of array_like, optional
        Lower and upper bounds on independent variables. Defaults to no bounds.
        Each bound must match the size of `x0` or be a scalar, in the latter
        case the bound will be the same for all variables. Use it to limit the
        range of function evaluation.
    args, kwargs : tuple and dict, optional
        Additional arguments passed to `fun` and `jac`. Both empty by default.
        The calling signature is ``fun(x, *args, **kwargs)`` and the same
        for `jac`.

    Returns
    -------
    accuracy : float
        The maximum among all relative errors for elements with absolute values
        higher than 1 and absolute errors for elements with absolute values
        less or equal than 1. If `accuracy` is on the order of 1e-6 or lower,
        then it is likely that your `jac` implementation is correct.

    See Also
    --------
    approx_derivative : Compute finite difference approximation of derivative.

    Examples
    --------
    >>> import numpy as np
    >>> from scipy.optimize._numdiff import check_derivative
    >>>
    >>>
    >>> def f(x, c1, c2):
    ...     return np.array([x[0] * np.sin(c1 * x[1]),
    ...                      x[0] * np.cos(c2 * x[1])])
    ...
    >>> def jac(x, c1, c2):
    ...     return np.array([
    ...         [np.sin(c1 * x[1]),  c1 * x[0] * np.cos(c1 * x[1])],
    ...         [np.cos(c2 * x[1]), -c2 * x[0] * np.sin(c2 * x[1])]
    ...     ])
    ...
    >>>
    >>> x0 = np.array([1.0, 0.5 * np.pi])
    >>> check_derivative(f, jac, x0, args=(1, 2))
    2.4492935982947064e-16
    )rb   r   r   r   r   )rb   r   r   )
r   r   r
   r   r   r[   r   r   r   r"   )r   jacr$   rb   r   r   	J_to_testJ_diffabs_errr   r   abs_err_dataJ_diff_datas                r7   check_derivativer   v  s   z ~B(((I	"36I(,V=i(	f$!']1ljj1.446vvbff\*jjBFF;$789 : 	: #36(,V=&&V+,vvg

1bffVn ==>>r9   )r   )&__doc__	functoolsnumpyr   numpy.linalgr   scipy.sparse.linalgr   sparser   r   r   r	   r
   r   _group_columnsr   r   scipy._lib._array_apir   scipy._lib._utilr   
scipy._libr   ry   r8   	lru_cacherL   rX   rc   rs   r    r   r   r   r   r}   r   rt   r9   r7   <module>r     s    -    . Q Q 5 1 ' -L%^ 2; 2;j.b*:z '0$w&7$).R"'Sl
(-VP frRj . -/FF7BFF*;" M?r9   