B d @sFdZddlZddlZddlmZddlmZddlmZddlm Z ddlmZ ddl m Z m Z ddlmZd d lmZmZmZd d d ddddddddddg ZGdd d ZGdd d ZGdd d ZGdddZGdddZd.ddZd/ddZd d!Zd0d$d%Zd&dZd'dZd(dZ d)d*Z!d+dZ"d,dZ#d-dZ$dS)1a Bounding classes used when proposing new live points, along with a number of useful helper functions. Bounding objects include: UnitCube: The unit N-cube (unconstrained draws from the prior). Ellipsoid: Bounding ellipsoid. MultiEllipsoid: A set of (possibly overlapping) bounding ellipsoids. RadFriends: A set of (possibly overlapping) balls centered on each live point. SupFriends: A set of (possibly overlapping) cubes centered on each live point. N)linalg)cov)spatial)cluster) logsumexpgammaln)kmeans2) unitcheckget_seed_sequenceget_random_generatorUnitCube EllipsoidMultiEllipsoid RadFriends SupFriendslogvol_prefactor randspherebounding_ellipsoidbounding_ellipsoids_bounding_ellipsoids_ellipsoid_bootstrap_expand_friends_bootstrap_radius_friends_leaveoneout_radiusc@s>eZdZdZddZddZdddZdd d Zdd d ZdS)r z An N-dimensional unit cube. Parameters ---------- ndim : int The number of dimensions of the unit cube. cCs||_d|_d|_dS)Ngg?)nlogvolfunit)selfndimr]/work/yifan.wang/ringdown/master-ringdown-env/lib/python3.7/site-packages/dynesty/bounding.py__init__7szUnitCube.__init__cCst|S)z+Checks if unit cube contains the point `x`.)r )rxrrr contains<szUnitCube.containsNcCs|j|jdS)z Draw a sample uniformly distributed within the unit cube. Returns ------- x : `~numpy.ndarray` with shape (ndim,) A coordinate within the unit cube. )size)randomr)rrstaterrr sampleAs zUnitCube.samplecCs|j||jfdS)z Draw `nsamples` samples randomly distributed within the unit cube. Returns ------- x : `~numpy.ndarray` with shape (nsamples, ndim) A collection of coordinates within the unit cube. )r$)r%r)rnsamplesr&rrr samplesNs zUnitCube.samplesrcCsdS)zFiller function.Nr)rpointsr& bootstrappoolrrr update[szUnitCube.update)N)N)NrN) __name__ __module__ __qualname____doc__r!r#r'r)r-rrrr r ,s   c@sjeZdZdZdddZddZddZd d Zd d Zd dZ dddZ dddZ dddZ dddZ dS)ra An N-dimensional ellipsoid defined by:: (x - v)^T A (x - v) = 1 where the vector `v` is the center of the ellipsoid and `A` is a symmetric, positive-definite `N x N` matrix. Parameters ---------- ctr : `~numpy.ndarray` with shape (N,) Coordinates of ellipsoid center. cov : `~numpy.ndarray` with shape (N, N) Covariance matrix describing the axes. NcCst||_t||_t||_tj|jdd\}}t|dkt |@rzt ||_ t |jdt ||_ntd|jd|d|d|dkr||j |_n||_|dkr|d ||j|_n||_||j |_d |_d |_dS) NF) check_finitegg?z2The input precision matrix defining the ellipsoid z is apparently singular with l=z and v=.g?r )lenrnpZasarrayctrrlalgeighallisfinitesqrtaxlensrlogsumr ValueErroraxesTampaxesexpandr)rr6rrBr@lvrrr r!ss$      zEllipsoid.__init__cCsr||j}tt|jd}t|j}||||jkrt||j}|j|d9_|j d|d9_ |j|9_|j |9_ nt |j}|}|j}t j |jdd\} } xPt| dddD]8} tt||| ||d} | || <|| 8}|d8}qWt|} | | d}| || j|_| d|| j|_ |j| 9_|j | |_ ||_dS) z#Scale ellipsoid to a target volume.g?F)r2Nrr )rr5r=r;rr<maxexprrBr@zerosr7r8ZargsortminrA)rrZlogfZ max_log_axlenZ log_axlenfZlogfaxZcurlogfZcurnrErFZcurideltafaxl1rrr scale_to_logvols4       zEllipsoid.scale_to_logvolcCs2t|j}|jdd|f}|j||j|fS)z'Return the endpoints of the major axis.N)r5Zargmaxr<rCr6)rirFrrr major_axis_endpointss zEllipsoid.major_axis_endpointscCs&||j}ttt||j|S)zPCompute the normalized distance to `x` from the center of the ellipsoid.)r6r5r;dotrB)rr"drrr distances zEllipsoid.distancecCs.||jdddf}ttd||j|S)zPCompute the normalized distance to `x` from the center of the ellipsoid.Nz ij,jk,ik->i)r6r5r;einsumrB)rr"rUrrr distance_manyszEllipsoid.distance_manycCs||dkS)z!Checks if ellipsoid contains `x`.g?)rV)rr"rrr r#szEllipsoid.containscCs|jt|jt|j|dS)z Draw a sample uniformly distributed within the ellipsoid. Returns ------- x : `~numpy.ndarray` with shape (ndim,) A coordinate within the ellipsoid. )r&)r6r5rTr@rr)rr&rrr r's zEllipsoid.samplecs"tfddt|D}|S)z Draw `nsamples` samples uniformly distributed within the ellipsoid. Returns ------- x : `~numpy.ndarray` with shape (nsamples, ndim) A collection of coordinates within the ellipsoid. csg|]}jdqS))r&)r').0rR)r&rrr sz%Ellipsoid.samples..)r5arrayrange)rr(r&xsr)r&rr r)s zEllipsoid.samples'cs6fddt|D}tdd|D}d||S)zsUsing `ndraws` Monte Carlo draws, estimate the fraction of overlap between the ellipsoid and the unit cube.csg|]}jdqS))r&)r')rYrR)r&rrr rZsz.Ellipsoid.unitcube_overlap..css|]}t|VqdS)N)r )rYr"rrr sz-Ellipsoid.unitcube_overlap..g?)r\r>)rndrawsr&r)Zninr)r&rr unitcube_overlapszEllipsoid.unitcube_overlaprFcst}|j|_|j|_|j|_|j|_|j|_|j|_|j|_|j|_|j |_ |dkr|dkrft }n|j }ddt |D}fddt |D} t ||} t || | } t|t| } t| } | dkr|j|jt| }|||r|j|d|_dS)a Update the ellipsoid to bound the collection of points. Parameters ---------- points : `~numpy.ndarray` with shape (npoints, ndim) The set of points to bound. rstate : `~numpy.random.Generator`, optional `~numpy.random.Generator` instance. bootstrap : int, optional The number of bootstrapped realizations of the ellipsoid. The maximum distance to the set of points "left out" during each iteration is used to enlarge the resulting volumes. Default is `0`. pool : user-provided pool, optional Use this pool of workers to execute operations in parallel. mc_integrate : bool, optional Whether to use Monte Carlo methods to compute the effective overlap of the final ellipsoid with the unit cube. Default is `False`. rNcSsg|]}dqS)Fr)rYitrrr rZ5sz$Ellipsoid.update..csg|]}qSrr)rYrb)r*rr rZ6sg?)r&)rrr6rrBrr<r@rCrDmapr\r ziplistrrIr5r=rQrar)rr*r&r+r, mc_integrateellMmultispsseedsargsexpandsrDlvr)r*r r-s2"   zEllipsoid.update)NN)N)N)r^N)NrNF)r.r/r0r1r!rQrSrVrXr#r'r)rar-rrrr r`s '$  c@sveZdZdZdddZddZddZd d Zdd d Zdd dZ ddZ d ddZ d!ddZ d"ddZ d#ddZdS)$ra A collection of M N-dimensional ellipsoids. Parameters ---------- ells : list of `Ellipsoid` objects with length M, optional A set of `Ellipsoid` objects that make up the collection of N-ellipsoids. Used to initialize :class:`MultiEllipsoid` if provided. ctrs : `~numpy.ndarray` with shape (M, N), optional Collection of coordinates of ellipsoid centers. Used to initialize :class:`MultiEllipsoid` if :data:`ams` is also provided. covs : `~numpy.ndarray` with shape (M, N, N), optional Collection of matrices describing the axes of the ellipsoids. Used to initialize :class:`MultiEllipsoid` if :data:`ctrs` also provided. Ncs|dk r4dkr*dkr*t||_||_qttdn@dkrNdkrNtdn&t|_fddt|jD|_|t|j|_t |j |_ d|_ d|_ dS)Nz5You cannot specific both `ells` and (`ctrs`, `covs`)!z3You must specify either `ells` or (`ctrs`, `covs`).csg|]}t||qSr)r)rYrR)covsctrsrr rZosz+MultiEllipsoid.__init__..g?r )r4nellsellsr?r\_MultiEllipsoid__update_arraysr5Zonesrmrlogvols logvol_tot expand_totr)rrrrpror)rorpr r!^s     zMultiEllipsoid.__init__cCsdtdd|jD|_tdd|jD|_tdd|jD|_tdd|jD|_dS)zI Update internal arrays to ensure that in sync with ells cSsg|] }|jqSr)r6)rYrgrrr rZ}sz2MultiEllipsoid.__update_arrays..cSsg|] }|jqSr)r)rYrgrrr rZ~scSsg|] }|jqSr)rB)rYrgrrr rZscSsg|] }|jqSr)r)rYrgrrr rZsN)r5r[rrrproamsrt)rrrr Z__update_arraysyszMultiEllipsoid.__update_arrayscs|x&tjD]}j|||q WtfddtjD_t|}j t |j 9_ |_ dS)zJScale ellipoids to a corresponding set of target volumes. csg|]}j|jqSr)rrrD)rYrR)rrr rZsz2MultiEllipsoid.scale_to_logvol..N) r\rqrrrQrsr5r[rmrrvrJru)rrtrRrur)rr rQszMultiEllipsoid.scale_to_logvolcCstdd|jDS)z9Return the endpoints of the major axis of each ellipsoid.cSsg|] }|qSr)rS)rYrgrrr rZsz7MultiEllipsoid.major_axis_endpoints..)r5r[rr)rrrr rSsz#MultiEllipsoid.major_axis_endpointscCsJ|dddf|j}td||j|dk}|dk rar Fr)rpr5rWrwZnonzero)rr"jdeltmaskrrr withins zMultiEllipsoid.withincCst|j||d}|S)zUChecks how many ellipsoid(s) `x` falls within, skipping the `j`-th ellipsoid.)rx)r4r{)rr"rxqrrr overlapszMultiEllipsoid.overlapcCs2|dddf|j}ttd||j|dkS)z-Checks if the set of ellipsoids contains `x`.Nz ai,aij,aj->ar )rpr5anyrWrw)rr"ryrrr r#szMultiEllipsoid.containsFc Cs|jdkr:|jdj|d}d}d}|r2|||fS||fSt|j|j}xt||}|j|j|d}|dddf|j}t d||j |}|dk }|dkr|dk }|dkr| } t d| |r|||fS|dks|d|krN||fSqNWdS)a Sample a point uniformly distributed within the *union* of ellipsoids. Returns ------- x : `~numpy.ndarray` with shape (ndim,) A coordinate within the set of ellipsoids. idx : int The index of the ellipsoid `x` was sampled from. q : int, optional The number of ellipsoids `x` falls within. r r)r&Nz ai,aij,aj->azEllipsoid check failed q=0, g?)rqrrr'r5rJrtru rand_choicerprWrwr>rI RuntimeErrorr%) rr&return_qr"idxr|ZprobsZdeltsZ ell_masksZmax_maskrrr r's0       zMultiEllipsoid.samplecs"tfddt|D}|S)a  Draw `nsamples` samples uniformly distributed within the *union* of ellipsoids. Returns ------- xs : `~numpy.ndarray` with shape (nsamples, ndim) A collection of coordinates within the set of ellipsoids. csg|]}jddqS))r&r)r')rYrR)r&rrr rZsz*MultiEllipsoid.samples..)r5r[r\)rr(r&r]r)r&rr r)s zMultiEllipsoid.samples'Tc slfddt|D}tdd|D}t||j}|rdtdd|D}||}||fS|SdS)zUsing `ndraws` Monte Carlo draws, estimate the log volume of the *union* of ellipsoids. If `return_overlap=True`, also returns the estimated fractional overlap with the unit cube.csg|]}jddqS)T)r&r)r')rYrR)r&rrr rZsz5MultiEllipsoid.monte_carlo_logvol..css|]\}}}d|VqdS)g?Nr)rYr"rr|rrr r_sz4MultiEllipsoid.monte_carlo_logvol..css$|]\}}}d|t|VqdS)g?N)r )rYr"rr|rrr r_ sN)r\r>r5r=ru) rr`r&return_overlapr)qsumrqinr}r)r&rr monte_carlo_logvols z!MultiEllipsoid.monte_carlo_logvolrcshj\}}|dkrtdt}t|} t| _| _tfddDsbtdt j _ t ddjD} j t | } t | } t j | _|dkrH|d krt} n|j} d dt|D}fd dt|D}t||}t|||}t| t|} t| }|d krHj |t |}||rdj|d d\_ _d S)a Update the set of ellipsoids to bound the collection of points. Parameters ---------- points : `~numpy.ndarray` with shape (npoints, ndim) The set of points to bound. rstate : `~numpy.random.Generator`, optional `~numpy.random.Generator` instance. bootstrap : int, optional The number of bootstrapped realizations of the ellipsoids. The maximum distance to the set of points "left out" during each iteration is used to enlarge the resulting volumes. Default is `0`. pool : user-provided pool, optional Use this pool of workers to execute operations in parallel. mc_integrate : bool, optional Whether to use Monte Carlo methods to compute the effective volume and fractional overlap of the final union of ellipsoids with the unit cube. Default is `False`. r z8Cannot compute the bounding ellipsoid of a single point.c3s|]}|VqdS)N)r#)rYp)rrr r_Dsz(MultiEllipsoid.update..z'Rejecting invalid MultiEllipsoid regioncSsg|] }|jqSr)rD)rYrgrrr rZJsz)MultiEllipsoid.update..rNcSsg|]}dqS)Tr)rYrbrrr rZWscsg|]}qSrr)rYrb)r*rr rZXsg?T)r&r)shaperrrr4rqrrrsr9rrtrur5r[r=rJrvrcr\r rdrerrIrQrr)rr*r&r+r,rfnpointsrZfirstellrrrmZ logvols_origZlogvol_tot_origrhrirjrkrlrDZlvsr)r*rr r-s>!         zMultiEllipsoid.update)NNN)N)N)NF)N)rNT)NrNF)r.r/r0r1r!rsrQrSr{r}r#r'r)rr-rrrr rJs"    <  c@sreZdZdZdddZddZddZd d Zd d ZdddZ dddZ d ddZ d!ddZ ddZ ddZdS)"ra0 A collection of N-balls of identical size centered on each live point. Parameters ---------- ndim : int The number of dimensions of each ball. cov : `~numpy.ndarray` with shape `(ndim, ndim)`, optional Covariance structure (correlation and size) of each ball. NcCs||_|dkrt|j}||_t|j|_t|j|_t|j|_ t |j\}}|dksft t |jd||_d|_d|_dS)Nrg?g?r )rr5identityrr7pinvhrBsqrtmr@axes_invrslogdetAssertionErrorr logvol_ballrDr)rrrdetsigndetlnrrr r!{s  zRadFriends.__init__cCsdt||jd|j}|j|d9_|j|d_|j|9_|j|_||_dS)z(Scale ball to encompass a target volume.g?rGN)r5rJrrrrBr@r)rrrMrrr rQs zRadFriends.scale_to_logvolcCs.ttjt|||jdddkd}|S)z#Check which balls `x` falls within.r )axisg?r)r5wherer7normrTr)rr"rpidxsrrr r{s&zRadFriends.withincCst|||}|S)z&Check how many balls `x` falls within.)r4r{)rr"rpr|rrr r}szRadFriends.overlapcCs|||dkS)z'Check if the set of balls contains `x`.r)r})rr"rprrr r#szRadFriends.containsFc Cst|}xt|j|d}t||j}|dkr@d}|d|}n"||} || |}|||}|dks~|s~|d|kr |r||fS|Sq WdS)a$ Sample a point uniformly distributed within the *union* of balls. Returns ------- x : `~numpy.ndarray` with shape (ndim,) A coordinate within the set of balls. q : int, optional The number of balls `x` falls within. )r&r rg?N) r4rrr5rTr@integersr}r%) rrpr&rnctrsdsdxr|r"rrrr r's   zRadFriends.samplecs$tfddt|D}|S)a Draw `nsamples` samples uniformly distributed within the *union* of balls. Returns ------- xs : `~numpy.ndarray` with shape (nsamples, ndim) A collection of coordinates within the set of balls. csg|]}jdqS))r&)r')rYrR)rpr&rrr rZsz&RadFriends.samples..)r5r[r\)rr(rpr&r]r)rpr&rr r)s zRadFriends.samples'Tc sfddt|D}tdd|D}td|}td||tj}|rtdd|D} | |} || fS|SdS)zUsing `ndraws` Monte Carlo draws, estimate the log volume of the *union* of balls. If `return_overlap=True`, also returns the estimated fractional overlap with the unit cube.csg|]}jddqS)T)r&r)r')rYrR)rpr&rrr rZsz1RadFriends.monte_carlo_logvol..cSsg|] }|dqS)r r)rY_rrr rZsg?css"|]\}}d|t|VqdS)g?N)r )rYr"r|rrr r_sz0RadFriends.monte_carlo_logvol..N)r\r5r[r>r=r4r) rrpr`r&rr)qsrrrr}r)rpr&rr rs   zRadFriends.monte_carlo_logvolrcsj|dkrt}n|j}|r&|||_n |||_t|j|_t|j|_t|j|_ t ||j |dkr~t d}nLfddt |D} ddt |D} t||} t| | | } t|t| }t|} |j| d9_|j| d_|j| 9_|j | _ t|j\}}|dks0tt|jd ||_d |_|rf|j|d |d d |_dS)a/ Update the radii of our balls. Parameters ---------- points : `~numpy.ndarray` with shape (npoints, ndim) The set of points to bound. rstate : `~numpy.random.Generator`, optional `~numpy.random.Generator` instance. bootstrap : int, optional The number of bootstrapped realizations of the ellipsoids. The maximum distance to the set of points "left out" during each iteration is used to enlarge the resulting volumes. Default is `0`. pool : user-provided pool, optional Use this pool of workers to execute operations in parallel. mc_integrate : bool, optional Whether to use Monte Carlo methods to compute the effective volume and fractional overlap of the final union of balls with the unit cube. Default is `False`. use_clustering : bool, optional Whether to use clustering to avoid issues with widely-seperated modes. Default is `True`. Ngballscsg|]}qSrr)rYrb)points_trr rZ8sz%RadFriends.update..cSsg|]}dqS)rr)rYrbrrr rZ9srGrg?g?T)rr&r )rc_get_covariance_from_clustersr_get_covariance_from_all_pointsr7rrBrr@rr5rTrr\r rdrerrIrrrrrrrDrr)rr*r&r+r,rfuse_clusteringrhZradiirjftypesrkrlZrmaxrrr)rr r-s<'    zRadFriends.updatecCstj|ddS)z$Compute covariance using all points.F)rowvar)r5r)rr*rrr rTsz*RadFriends._get_covariance_from_all_pointsc Cstjj|d|jd}tj|}tjj|ddd}t |}|dkrN| |Sd}t |}x`t |D]R}|||kddf} | j dd d } |t| } | | ||| ddf<| }qhW| |SdS) z-Compute covariance from re-centered clusters. mahalanobis)metricVIg?rV) criterionr rN)r)r rH)rrVpdistrBr hierarchysinglefclusterr5rIr empty_likeuniquemeanreshaper4) rr* distanceslinkages clusteridxs nclustersrRoverlapped_pointsr group_points group_meanrxrrr rYs&      z(RadFriends._get_covariance_from_clusters)N)NF)N)rNT)NrNFT)r.r/r0r1r!rQr{r}r#r'r)rr-rrrrrr rms$     $   Sc@sreZdZdZdddZddZddZd d Zd d ZdddZ dddZ d ddZ d!ddZ ddZ ddZdS)"ra/ A collection of N-cubes of identical size centered on each live point. Parameters ---------- ndim : int The number of dimensions of the cube. cov : `~numpy.ndarray` with shape `(ndim, ndim)`, optional Covariance structure (correlation and size) of each cube. NcCs||_|dkrt|j}||_t|j|_t|j|_t|j|_ t |j\}}|dksft |jt dd||_d|_d|_dS)Nrg@g?g?r )rr5rrr7rrBrr@rrrrr= logvol_cuberDr)rrrrrrrr r!s  zSupFriends.__init__cCsdt||jd|j}|j|d9_|j|d_|j|9_|j|_||_dS)z(Scale cube to encompass a target volume.g?rGN)r5rJrrrrBr@r)rrrMrrr rQs zSupFriends.scale_to_logvolc Cs4ttjtt|||jdddkd}|S)z$Checks which cubes `x` falls within.r )rg?r)r5rrIabsrTr)rr"rprrrr r{s,zSupFriends.withincCst|||}|S)zIChecks how many cubes `x` falls within, skipping the `j`-th cube.)r4r{)rr"rpr|rrr r}szSupFriends.overlapcCs|||dkS)z(Checks if the set of cubes contains `x`.r)r})rr"rprrr r#szSupFriends.containsFc Cst|}x|jdd|jd}t||j}|dkrD|d|}d}n"||} || |}|||}|dks|s|d|kr |r||fS|Sq WdS)a$ Sample a point uniformly distributed within the *union* of cubes. Returns ------- x : `~numpy.ndarray` with shape (ndim,) A coordinate within the set of cubes. q : int, optional The number of cubes `x` falls within. rHr )r$rg?N) r4uniformrr5rTr@rr}r%) rrpr&rrrrr"r|rrrr r's    zSupFriends.samplecs$tfddt|D}|S)a Draw `nsamples` samples uniformly distributed within the *union* of cubes. Returns ------- xs : `~numpy.ndarray` with shape (nsamples, ndim) A collection of coordinates within the set of cubes. csg|]}jdqS))r&)r')rYrR)rpr&rrr rZsz&SupFriends.samples..)r5r[r\)rr(rpr&r]r)rpr&rr r)s zSupFriends.samples'Tc szfddt|D}tdd|D}td||tj}|rrtdd|D}||} || fS|SdS)zUsing `ndraws` Monte Carlo draws, estimate the log volume of the *union* of cubes. If `return_overlap=True`, also returns the estimated fractional overlap with the unit cube.csg|]}jddqS)T)r&r)r')rYrR)rpr&rrr rZsz1SupFriends.monte_carlo_logvol..css|]\}}d|VqdS)g?Nr)rYr"r|rrr r_sz0SupFriends.monte_carlo_logvol..g?css"|]\}}d|t|VqdS)g?N)r )rYr"r|rrr r_sN)r\r>r5r=r4r) rrpr`r&rr)rrrr}r)rpr&rr rs   zSupFriends.monte_carlo_logvolrcsp|dkrt}n|j}|r&|||_n |||_t|j|_t|j|_t|j|_ t ||j |dkr~t d}nLfddt |D} ddt |D} t||} t| | | } t|t| }t|} |j| d9_|j| d_|j| 9_|j | _ t|j\}}|dks0t|jt d d ||_d |_|rl|j|d |d d|_dS)a; Update the half-side-lengths of our cubes. Parameters ---------- points : `~numpy.ndarray` with shape (npoints, ndim) The set of points to bound. rstate : `~numpy.random.Generator`, optional `~numpy.random.Generator` instance. bootstrap : int, optional The number of bootstrapped realizations of the ellipsoids. The maximum distance to the set of points "left out" during each iteration is used to enlarge the resulting volumes. Default is `0`. pool : user-provided pool, optional Use this pool of workers to execute operations in parallel. mc_integrate : bool, optional Whether to use Monte Carlo methods to compute the effective volume and fractional overlap of the final union of cubes with the unit cube. Default is `False`. use_clustering : bool, optional Whether to use clustering to avoid issues with widely-seperated modes. Default is `True`. Ngcubescsg|]}qSrr)rYrb)rrr rZCsz%SupFriends.update..cSsg|]}dqS)rr)rYrbrrr rZDsrGrg@g?g?T)rr&r )rcrrrr7rrBrr@rr5rTrr\r rdrerrIrrrrr=rrDrr)rr*r&r+r,rfrrhZhsidesrjrrkrlZhsmaxrrr)rr r-s<'    zSupFriends.updatecCstj|ddS)z$Compute covariance using all points.F)r)r5r)rr*rrr r^sz*SupFriends._get_covariance_from_all_pointsc Cstjj|d|jd}tj|}tjj|ddd}t |}|dkrN| |Sd}t |}x`t |D]R}|||kddf} | j dd d } |t| } | | ||| ddf<| }qhW| |SdS) z-Compute covariance from re-centered clusters.r)rrg?rV)rr rN)r)r rH)rrVrrBrrrrr5rIrrrrrr4) rr*rrrrrRrrrrrxrrr rcs&      z(SupFriends._get_covariance_from_clusters)N)NF)N)rNT)NrNFT)r.r/r0r1r!rQr{r}r#r'r)rr-rrrrrr rws$     %   R@cCs>|d9}|td|td|dt||d}|S)a Returns the ln(volume constant) for an `n`-dimensional sphere with an :math:`L^p` norm. The constant is defined as:: lnf = n * ln(2.) + n * LogGamma(1./p + 1) - LogGamma(n/p + 1.) By default the `p=2.` norm is used (i.e. the standard Euclidean norm). g?g@r )r5r=r)rrZlnfrrr rs 2cCs2|j|d}||d|tj|dd}|S)z=Draw a point uniformly within an `n`-dimensional unit sphere.)r$g?F)r2)Zstandard_normalr%r7r)rr&zZxhatrrr rs "cCs,t|}|}tt||t|dS)z Optimized version of numpy's random.choice Return an index of a point selected with the probability pb The pb must sum to 1 r )r5Zcumsumr%rLZ searchsortedr4)Zpbr&p1Zxrrrr rs rdmBc Csx|jd}t|}d}d}xt|D]}d}yhtj|dd\} } | } | } t|  r| dkrpd}q| | |krd}q| | d} Pnd}Wntj k rd}YnX|dkr(|dkrt | || |}| || j }q(|d ||d |d}d |||t |}q(W|dkrRtd t |}|}|} n| d | | j }|dk}|||| fS) a7 Given the covariance matrix improve it, if it is not invertable or eigen values are negative or condition number that is above the limit Returns: a tuple with three elements 1) a boolean flag if a matrix is 'good', so it didn't need adjustments 2) updated matrix 3) its inverse rg|= F)r2rGr g?g?zTFailed to guarantee the ellipsoid axes will be non-singular. Defaulting to a sphere.)rr5r[r\r7r8rIrLr:r9Z LinAlgErrormaximumrAeyewarningswarncopy)Zcovar0ZntriesZmax_condition_numberrcovarZcoeffminZeig_multZtrialfailedZeigvalZeigvecmaxvalminvalr@Z eigval_fixZcoeffrBgood_matrrr improve_covar_matsF         rcCs|j\}}|dkrtdtj|dd}t|dd}||}|dkrNt|}d}d|}xtd D]|}t|\} }} } td || | } |dkr| |kr| |} || 9}| | } | t | 9} |dkr| dkrt d | rdPqdWt ||| | d }|S) a0 Calculate the bounding ellipsoid containing a collection of points. Parameters ---------- points : `~numpy.ndarray` with shape (npoints, ndim) A set of coordinates. Returns ------- ellipsoid : :class:`Ellipsoid` The bounding :class:`Ellipsoid` object. r z6Cannot compute a bounding ellipsoid of a single point.r)rF)rgMbP?g?rGz ij,jk,ik->iz.)rr rrGcSsg|] }t|qSr)r)rYZpoints_jrrr rZyscSsg|] }|jqSr)r)rYerrr rZs)rrSr5Zvstackrcatch_warnings simplefilterrrJr=Z logaddexprrr4r)r*rgrrrp2Z start_ctrsZk2_resZpoints_krrZnparamZvol_dec1outZvol_dec2r)rr*r rGs6    , &$cCst|}t||}t|dS)aa Calculate a set of ellipsoids that bound the collection of points. Parameters ---------- points : `~numpy.ndarray` with shape (npoints, ndim) A set of coordinates. Returns ------- mell : :class:`MultiEllipsoid` object The :class:`MultiEllipsoid` object used to bound the collection of points. )rr)rrr)r*rgrrrrr rs c Cst|}|jd}|j||d}t|}tj|td}d||<|}|dkr\d|dd<||dkrpd|d<||||}} || fS) zl Select the bootstrap set from points. Return: Tuple with selected, and not-selected points r)r$)ZdtypeTrGNr F)r rrr5rrKboolr>) r*rseedr&rrZidx_inZsel_inZn_in points_in points_outrrr _bootstrap_pointss    rc sp|\}}}t||\}t|}|s0|}n,t||}tjtfdd|Ddd}tdt|}|S)a8Internal method used to compute the expansion factor for a bounding ellipsoid or ellipsoids based on bootstrapping. The argument is a tuple: multi: boolean flag if we are doing multiell or single ell decomposition points: 2d array of points rseed: seed to initialize the random generator csg|]}|qSr)rX)rYel)rrr rZsz/_ellipsoid_bootstrap_expand..r)rg?)rrrXrr5rLr[rI) rlmultir*rrrgdistsrrrDr)rr rs   c Csn|\}}}t||\}}t|}|dkrB|j|ddddd}n |dkrb|j|ddtjdd}t|}|S)zInternal method used to compute the radius (half-side-length) for each ball (cube) used in :class:`RadFriends` (:class:`SupFriends`) using bootstrapping.rr rrG)repsrr)rrKDTreequeryr5infrI) rlr*ftyperrrkdtreerdistrrr rs  cCs^t|}|dkr*|j|ddddd}n |dkrJ|j|ddtjdd}|dddf}|S)zInternal method used to compute the radius (half-side-length) for each ball (cube) used in :class:`RadFriends` (:class:`SupFriends`) using leave-one-out (LOO) cross-validation.rrGr)rrrrNr )rrrr5r)r*rrrrrrr rs )r)N)rr)%r1rnumpyr5rrrZscipyrrr7Z scipy.specialrrZscipy.cluster.vqrutilsr r r __all__r rrrrrrrrrrrrrrrrrrr sL       4k%    LMX