function_base.py 155 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735736737738739740741742743744745746747748749750751752753754755756757758759760761762763764765766767768769770771772773774775776777778779780781782783784785786787788789790791792793794795796797798799800801802803804805806807808809810811812813814815816817818819820821822823824825826827828829830831832833834835836837838839840841842843844845846847848849850851852853854855856857858859860861862863864865866867868869870871872873874875876877878879880881882883884885886887888889890891892893894895896897898899900901902903904905906907908909910911912913914915916917918919920921922923924925926927928929930931932933934935936937938939940941942943944945946947948949950951952953954955956957958959960961962963964965966967968969970971972973974975976977978979980981982983984985986987988989990991992993994995996997998999100010011002100310041005100610071008100910101011101210131014101510161017101810191020102110221023102410251026102710281029103010311032103310341035103610371038103910401041104210431044104510461047104810491050105110521053105410551056105710581059106010611062106310641065106610671068106910701071107210731074107510761077107810791080108110821083108410851086108710881089109010911092109310941095109610971098109911001101110211031104110511061107110811091110111111121113111411151116111711181119112011211122112311241125112611271128112911301131113211331134113511361137113811391140114111421143114411451146114711481149115011511152115311541155115611571158115911601161116211631164116511661167116811691170117111721173117411751176117711781179118011811182118311841185118611871188118911901191119211931194119511961197119811991200120112021203120412051206120712081209121012111212121312141215121612171218121912201221122212231224122512261227122812291230123112321233123412351236123712381239124012411242124312441245124612471248124912501251125212531254125512561257125812591260126112621263126412651266126712681269127012711272127312741275127612771278127912801281128212831284128512861287128812891290129112921293129412951296129712981299130013011302130313041305130613071308130913101311131213131314131513161317131813191320132113221323132413251326132713281329133013311332133313341335133613371338133913401341134213431344134513461347134813491350135113521353135413551356135713581359136013611362136313641365136613671368136913701371137213731374137513761377137813791380138113821383138413851386138713881389139013911392139313941395139613971398139914001401140214031404140514061407140814091410141114121413141414151416141714181419142014211422142314241425142614271428142914301431143214331434143514361437143814391440144114421443144414451446144714481449145014511452145314541455145614571458145914601461146214631464146514661467146814691470147114721473147414751476147714781479148014811482148314841485148614871488148914901491149214931494149514961497149814991500150115021503150415051506150715081509151015111512151315141515151615171518151915201521152215231524152515261527152815291530153115321533153415351536153715381539154015411542154315441545154615471548154915501551155215531554155515561557155815591560156115621563156415651566156715681569157015711572157315741575157615771578157915801581158215831584158515861587158815891590159115921593159415951596159715981599160016011602160316041605160616071608160916101611161216131614161516161617161816191620162116221623162416251626162716281629163016311632163316341635163616371638163916401641164216431644164516461647164816491650165116521653165416551656165716581659166016611662166316641665166616671668166916701671167216731674167516761677167816791680168116821683168416851686168716881689169016911692169316941695169616971698169917001701170217031704170517061707170817091710171117121713171417151716171717181719172017211722172317241725172617271728172917301731173217331734173517361737173817391740174117421743174417451746174717481749175017511752175317541755175617571758175917601761176217631764176517661767176817691770177117721773177417751776177717781779178017811782178317841785178617871788178917901791179217931794179517961797179817991800180118021803180418051806180718081809181018111812181318141815181618171818181918201821182218231824182518261827182818291830183118321833183418351836183718381839184018411842184318441845184618471848184918501851185218531854185518561857185818591860186118621863186418651866186718681869187018711872187318741875187618771878187918801881188218831884188518861887188818891890189118921893189418951896189718981899190019011902190319041905190619071908190919101911191219131914191519161917191819191920192119221923192419251926192719281929193019311932193319341935193619371938193919401941194219431944194519461947194819491950195119521953195419551956195719581959196019611962196319641965196619671968196919701971197219731974197519761977197819791980198119821983198419851986198719881989199019911992199319941995199619971998199920002001200220032004200520062007200820092010201120122013201420152016201720182019202020212022202320242025202620272028202920302031203220332034203520362037203820392040204120422043204420452046204720482049205020512052205320542055205620572058205920602061206220632064206520662067206820692070207120722073207420752076207720782079208020812082208320842085208620872088208920902091209220932094209520962097209820992100210121022103210421052106210721082109211021112112211321142115211621172118211921202121212221232124212521262127212821292130213121322133213421352136213721382139214021412142214321442145214621472148214921502151215221532154215521562157215821592160216121622163216421652166216721682169217021712172217321742175217621772178217921802181218221832184218521862187218821892190219121922193219421952196219721982199220022012202220322042205220622072208220922102211221222132214221522162217221822192220222122222223222422252226222722282229223022312232223322342235223622372238223922402241224222432244224522462247224822492250225122522253225422552256225722582259226022612262226322642265226622672268226922702271227222732274227522762277227822792280228122822283228422852286228722882289229022912292229322942295229622972298229923002301230223032304230523062307230823092310231123122313231423152316231723182319232023212322232323242325232623272328232923302331233223332334233523362337233823392340234123422343234423452346234723482349235023512352235323542355235623572358235923602361236223632364236523662367236823692370237123722373237423752376237723782379238023812382238323842385238623872388238923902391239223932394239523962397239823992400240124022403240424052406240724082409241024112412241324142415241624172418241924202421242224232424242524262427242824292430243124322433243424352436243724382439244024412442244324442445244624472448244924502451245224532454245524562457245824592460246124622463246424652466246724682469247024712472247324742475247624772478247924802481248224832484248524862487248824892490249124922493249424952496249724982499250025012502250325042505250625072508250925102511251225132514251525162517251825192520252125222523252425252526252725282529253025312532253325342535253625372538253925402541254225432544254525462547254825492550255125522553255425552556255725582559256025612562256325642565256625672568256925702571257225732574257525762577257825792580258125822583258425852586258725882589259025912592259325942595259625972598259926002601260226032604260526062607260826092610261126122613261426152616261726182619262026212622262326242625262626272628262926302631263226332634263526362637263826392640264126422643264426452646264726482649265026512652265326542655265626572658265926602661266226632664266526662667266826692670267126722673267426752676267726782679268026812682268326842685268626872688268926902691269226932694269526962697269826992700270127022703270427052706270727082709271027112712271327142715271627172718271927202721272227232724272527262727272827292730273127322733273427352736273727382739274027412742274327442745274627472748274927502751275227532754275527562757275827592760276127622763276427652766276727682769277027712772277327742775277627772778277927802781278227832784278527862787278827892790279127922793279427952796279727982799280028012802280328042805280628072808280928102811281228132814281528162817281828192820282128222823282428252826282728282829283028312832283328342835283628372838283928402841284228432844284528462847284828492850285128522853285428552856285728582859286028612862286328642865286628672868286928702871287228732874287528762877287828792880288128822883288428852886288728882889289028912892289328942895289628972898289929002901290229032904290529062907290829092910291129122913291429152916291729182919292029212922292329242925292629272928292929302931293229332934293529362937293829392940294129422943294429452946294729482949295029512952295329542955295629572958295929602961296229632964296529662967296829692970297129722973297429752976297729782979298029812982298329842985298629872988298929902991299229932994299529962997299829993000300130023003300430053006300730083009301030113012301330143015301630173018301930203021302230233024302530263027302830293030303130323033303430353036303730383039304030413042304330443045304630473048304930503051305230533054305530563057305830593060306130623063306430653066306730683069307030713072307330743075307630773078307930803081308230833084308530863087308830893090309130923093309430953096309730983099310031013102310331043105310631073108310931103111311231133114311531163117311831193120312131223123312431253126312731283129313031313132313331343135313631373138313931403141314231433144314531463147314831493150315131523153315431553156315731583159316031613162316331643165316631673168316931703171317231733174317531763177317831793180318131823183318431853186318731883189319031913192319331943195319631973198319932003201320232033204320532063207320832093210321132123213321432153216321732183219322032213222322332243225322632273228322932303231323232333234323532363237323832393240324132423243324432453246324732483249325032513252325332543255325632573258325932603261326232633264326532663267326832693270327132723273327432753276327732783279328032813282328332843285328632873288328932903291329232933294329532963297329832993300330133023303330433053306330733083309331033113312331333143315331633173318331933203321332233233324332533263327332833293330333133323333333433353336333733383339334033413342334333443345334633473348334933503351335233533354335533563357335833593360336133623363336433653366336733683369337033713372337333743375337633773378337933803381338233833384338533863387338833893390339133923393339433953396339733983399340034013402340334043405340634073408340934103411341234133414341534163417341834193420342134223423342434253426342734283429343034313432343334343435343634373438343934403441344234433444344534463447344834493450345134523453345434553456345734583459346034613462346334643465346634673468346934703471347234733474347534763477347834793480348134823483348434853486348734883489349034913492349334943495349634973498349935003501350235033504350535063507350835093510351135123513351435153516351735183519352035213522352335243525352635273528352935303531353235333534353535363537353835393540354135423543354435453546354735483549355035513552355335543555355635573558355935603561356235633564356535663567356835693570357135723573357435753576357735783579358035813582358335843585358635873588358935903591359235933594359535963597359835993600360136023603360436053606360736083609361036113612361336143615361636173618361936203621362236233624362536263627362836293630363136323633363436353636363736383639364036413642364336443645364636473648364936503651365236533654365536563657365836593660366136623663366436653666366736683669367036713672367336743675367636773678367936803681368236833684368536863687368836893690369136923693369436953696369736983699370037013702370337043705370637073708370937103711371237133714371537163717371837193720372137223723372437253726372737283729373037313732373337343735373637373738373937403741374237433744374537463747374837493750375137523753375437553756375737583759376037613762376337643765376637673768376937703771377237733774377537763777377837793780378137823783378437853786378737883789379037913792379337943795379637973798379938003801380238033804380538063807380838093810381138123813381438153816381738183819382038213822382338243825382638273828382938303831383238333834383538363837383838393840384138423843384438453846384738483849385038513852385338543855385638573858385938603861386238633864386538663867386838693870387138723873387438753876387738783879388038813882388338843885388638873888388938903891389238933894389538963897389838993900390139023903390439053906390739083909391039113912391339143915391639173918391939203921392239233924392539263927392839293930393139323933393439353936393739383939394039413942394339443945394639473948394939503951395239533954395539563957395839593960396139623963396439653966396739683969397039713972397339743975397639773978397939803981398239833984398539863987398839893990399139923993399439953996399739983999400040014002400340044005400640074008400940104011401240134014401540164017401840194020402140224023402440254026402740284029403040314032403340344035403640374038403940404041404240434044404540464047404840494050405140524053405440554056405740584059406040614062406340644065406640674068406940704071407240734074407540764077407840794080408140824083408440854086408740884089409040914092409340944095409640974098409941004101410241034104410541064107410841094110411141124113411441154116411741184119412041214122412341244125412641274128412941304131413241334134413541364137413841394140414141424143414441454146414741484149415041514152415341544155415641574158415941604161416241634164416541664167416841694170417141724173417441754176417741784179418041814182418341844185418641874188418941904191419241934194419541964197419841994200420142024203420442054206420742084209421042114212421342144215421642174218421942204221422242234224422542264227422842294230423142324233423442354236423742384239424042414242424342444245424642474248424942504251425242534254425542564257425842594260426142624263426442654266426742684269427042714272427342744275427642774278427942804281428242834284428542864287428842894290429142924293429442954296429742984299430043014302430343044305430643074308430943104311431243134314431543164317431843194320432143224323432443254326432743284329433043314332433343344335433643374338433943404341434243434344434543464347434843494350435143524353435443554356435743584359436043614362436343644365436643674368436943704371437243734374437543764377437843794380438143824383438443854386438743884389439043914392439343944395439643974398439944004401440244034404440544064407440844094410441144124413441444154416441744184419442044214422442344244425442644274428442944304431443244334434443544364437443844394440444144424443444444454446444744484449445044514452445344544455445644574458445944604461446244634464446544664467446844694470447144724473447444754476447744784479448044814482448344844485448644874488448944904491449244934494449544964497449844994500450145024503450445054506450745084509451045114512451345144515451645174518451945204521452245234524452545264527452845294530453145324533453445354536453745384539454045414542454345444545454645474548454945504551455245534554455545564557455845594560456145624563456445654566456745684569457045714572457345744575457645774578457945804581458245834584458545864587458845894590459145924593459445954596459745984599460046014602460346044605460646074608460946104611461246134614461546164617461846194620462146224623462446254626462746284629463046314632463346344635463646374638463946404641464246434644464546464647464846494650465146524653465446554656465746584659466046614662466346644665466646674668466946704671467246734674467546764677467846794680468146824683468446854686468746884689469046914692469346944695469646974698469947004701470247034704470547064707470847094710471147124713471447154716471747184719472047214722472347244725472647274728472947304731473247334734473547364737473847394740474147424743474447454746474747484749475047514752475347544755475647574758475947604761476247634764476547664767476847694770477147724773477447754776477747784779478047814782478347844785478647874788478947904791479247934794479547964797479847994800480148024803480448054806480748084809481048114812481348144815481648174818481948204821482248234824482548264827482848294830483148324833483448354836483748384839484048414842484348444845484648474848484948504851485248534854485548564857485848594860
  1. import collections.abc
  2. import functools
  3. import re
  4. import sys
  5. import warnings
  6. import numpy as np
  7. import numpy.core.numeric as _nx
  8. from numpy.core import transpose
  9. from numpy.core.numeric import (
  10. ones, zeros, arange, concatenate, array, asarray, asanyarray, empty,
  11. ndarray, around, floor, ceil, take, dot, where, intp,
  12. integer, isscalar, absolute
  13. )
  14. from numpy.core.umath import (
  15. pi, add, arctan2, frompyfunc, cos, less_equal, sqrt, sin,
  16. mod, exp, not_equal, subtract
  17. )
  18. from numpy.core.fromnumeric import (
  19. ravel, nonzero, partition, mean, any, sum
  20. )
  21. from numpy.core.numerictypes import typecodes
  22. from numpy.core.overrides import set_module
  23. from numpy.core import overrides
  24. from numpy.core.function_base import add_newdoc
  25. from numpy.lib.twodim_base import diag
  26. from numpy.core.multiarray import (
  27. _insert, add_docstring, bincount, normalize_axis_index, _monotonicity,
  28. interp as compiled_interp, interp_complex as compiled_interp_complex
  29. )
  30. from numpy.core.umath import _add_newdoc_ufunc as add_newdoc_ufunc
  31. import builtins
  32. # needed in this module for compatibility
  33. from numpy.lib.histograms import histogram, histogramdd
  34. array_function_dispatch = functools.partial(
  35. overrides.array_function_dispatch, module='numpy')
  36. __all__ = [
  37. 'select', 'piecewise', 'trim_zeros', 'copy', 'iterable', 'percentile',
  38. 'diff', 'gradient', 'angle', 'unwrap', 'sort_complex', 'disp', 'flip',
  39. 'rot90', 'extract', 'place', 'vectorize', 'asarray_chkfinite', 'average',
  40. 'bincount', 'digitize', 'cov', 'corrcoef',
  41. 'msort', 'median', 'sinc', 'hamming', 'hanning', 'bartlett',
  42. 'blackman', 'kaiser', 'trapz', 'i0', 'add_newdoc', 'add_docstring',
  43. 'meshgrid', 'delete', 'insert', 'append', 'interp', 'add_newdoc_ufunc',
  44. 'quantile'
  45. ]
  46. def _rot90_dispatcher(m, k=None, axes=None):
  47. return (m,)
  48. @array_function_dispatch(_rot90_dispatcher)
  49. def rot90(m, k=1, axes=(0, 1)):
  50. """
  51. Rotate an array by 90 degrees in the plane specified by axes.
  52. Rotation direction is from the first towards the second axis.
  53. Parameters
  54. ----------
  55. m : array_like
  56. Array of two or more dimensions.
  57. k : integer
  58. Number of times the array is rotated by 90 degrees.
  59. axes: (2,) array_like
  60. The array is rotated in the plane defined by the axes.
  61. Axes must be different.
  62. .. versionadded:: 1.12.0
  63. Returns
  64. -------
  65. y : ndarray
  66. A rotated view of `m`.
  67. See Also
  68. --------
  69. flip : Reverse the order of elements in an array along the given axis.
  70. fliplr : Flip an array horizontally.
  71. flipud : Flip an array vertically.
  72. Notes
  73. -----
  74. rot90(m, k=1, axes=(1,0)) is the reverse of rot90(m, k=1, axes=(0,1))
  75. rot90(m, k=1, axes=(1,0)) is equivalent to rot90(m, k=-1, axes=(0,1))
  76. Examples
  77. --------
  78. >>> m = np.array([[1,2],[3,4]], int)
  79. >>> m
  80. array([[1, 2],
  81. [3, 4]])
  82. >>> np.rot90(m)
  83. array([[2, 4],
  84. [1, 3]])
  85. >>> np.rot90(m, 2)
  86. array([[4, 3],
  87. [2, 1]])
  88. >>> m = np.arange(8).reshape((2,2,2))
  89. >>> np.rot90(m, 1, (1,2))
  90. array([[[1, 3],
  91. [0, 2]],
  92. [[5, 7],
  93. [4, 6]]])
  94. """
  95. axes = tuple(axes)
  96. if len(axes) != 2:
  97. raise ValueError("len(axes) must be 2.")
  98. m = asanyarray(m)
  99. if axes[0] == axes[1] or absolute(axes[0] - axes[1]) == m.ndim:
  100. raise ValueError("Axes must be different.")
  101. if (axes[0] >= m.ndim or axes[0] < -m.ndim
  102. or axes[1] >= m.ndim or axes[1] < -m.ndim):
  103. raise ValueError("Axes={} out of range for array of ndim={}."
  104. .format(axes, m.ndim))
  105. k %= 4
  106. if k == 0:
  107. return m[:]
  108. if k == 2:
  109. return flip(flip(m, axes[0]), axes[1])
  110. axes_list = arange(0, m.ndim)
  111. (axes_list[axes[0]], axes_list[axes[1]]) = (axes_list[axes[1]],
  112. axes_list[axes[0]])
  113. if k == 1:
  114. return transpose(flip(m, axes[1]), axes_list)
  115. else:
  116. # k == 3
  117. return flip(transpose(m, axes_list), axes[1])
  118. def _flip_dispatcher(m, axis=None):
  119. return (m,)
  120. @array_function_dispatch(_flip_dispatcher)
  121. def flip(m, axis=None):
  122. """
  123. Reverse the order of elements in an array along the given axis.
  124. The shape of the array is preserved, but the elements are reordered.
  125. .. versionadded:: 1.12.0
  126. Parameters
  127. ----------
  128. m : array_like
  129. Input array.
  130. axis : None or int or tuple of ints, optional
  131. Axis or axes along which to flip over. The default,
  132. axis=None, will flip over all of the axes of the input array.
  133. If axis is negative it counts from the last to the first axis.
  134. If axis is a tuple of ints, flipping is performed on all of the axes
  135. specified in the tuple.
  136. .. versionchanged:: 1.15.0
  137. None and tuples of axes are supported
  138. Returns
  139. -------
  140. out : array_like
  141. A view of `m` with the entries of axis reversed. Since a view is
  142. returned, this operation is done in constant time.
  143. See Also
  144. --------
  145. flipud : Flip an array vertically (axis=0).
  146. fliplr : Flip an array horizontally (axis=1).
  147. Notes
  148. -----
  149. flip(m, 0) is equivalent to flipud(m).
  150. flip(m, 1) is equivalent to fliplr(m).
  151. flip(m, n) corresponds to ``m[...,::-1,...]`` with ``::-1`` at position n.
  152. flip(m) corresponds to ``m[::-1,::-1,...,::-1]`` with ``::-1`` at all
  153. positions.
  154. flip(m, (0, 1)) corresponds to ``m[::-1,::-1,...]`` with ``::-1`` at
  155. position 0 and position 1.
  156. Examples
  157. --------
  158. >>> A = np.arange(8).reshape((2,2,2))
  159. >>> A
  160. array([[[0, 1],
  161. [2, 3]],
  162. [[4, 5],
  163. [6, 7]]])
  164. >>> np.flip(A, 0)
  165. array([[[4, 5],
  166. [6, 7]],
  167. [[0, 1],
  168. [2, 3]]])
  169. >>> np.flip(A, 1)
  170. array([[[2, 3],
  171. [0, 1]],
  172. [[6, 7],
  173. [4, 5]]])
  174. >>> np.flip(A)
  175. array([[[7, 6],
  176. [5, 4]],
  177. [[3, 2],
  178. [1, 0]]])
  179. >>> np.flip(A, (0, 2))
  180. array([[[5, 4],
  181. [7, 6]],
  182. [[1, 0],
  183. [3, 2]]])
  184. >>> A = np.random.randn(3,4,5)
  185. >>> np.all(np.flip(A,2) == A[:,:,::-1,...])
  186. True
  187. """
  188. if not hasattr(m, 'ndim'):
  189. m = asarray(m)
  190. if axis is None:
  191. indexer = (np.s_[::-1],) * m.ndim
  192. else:
  193. axis = _nx.normalize_axis_tuple(axis, m.ndim)
  194. indexer = [np.s_[:]] * m.ndim
  195. for ax in axis:
  196. indexer[ax] = np.s_[::-1]
  197. indexer = tuple(indexer)
  198. return m[indexer]
  199. @set_module('numpy')
  200. def iterable(y):
  201. """
  202. Check whether or not an object can be iterated over.
  203. Parameters
  204. ----------
  205. y : object
  206. Input object.
  207. Returns
  208. -------
  209. b : bool
  210. Return ``True`` if the object has an iterator method or is a
  211. sequence and ``False`` otherwise.
  212. Examples
  213. --------
  214. >>> np.iterable([1, 2, 3])
  215. True
  216. >>> np.iterable(2)
  217. False
  218. """
  219. try:
  220. iter(y)
  221. except TypeError:
  222. return False
  223. return True
  224. def _average_dispatcher(a, axis=None, weights=None, returned=None):
  225. return (a, weights)
  226. @array_function_dispatch(_average_dispatcher)
  227. def average(a, axis=None, weights=None, returned=False):
  228. """
  229. Compute the weighted average along the specified axis.
  230. Parameters
  231. ----------
  232. a : array_like
  233. Array containing data to be averaged. If `a` is not an array, a
  234. conversion is attempted.
  235. axis : None or int or tuple of ints, optional
  236. Axis or axes along which to average `a`. The default,
  237. axis=None, will average over all of the elements of the input array.
  238. If axis is negative it counts from the last to the first axis.
  239. .. versionadded:: 1.7.0
  240. If axis is a tuple of ints, averaging is performed on all of the axes
  241. specified in the tuple instead of a single axis or all the axes as
  242. before.
  243. weights : array_like, optional
  244. An array of weights associated with the values in `a`. Each value in
  245. `a` contributes to the average according to its associated weight.
  246. The weights array can either be 1-D (in which case its length must be
  247. the size of `a` along the given axis) or of the same shape as `a`.
  248. If `weights=None`, then all data in `a` are assumed to have a
  249. weight equal to one. The 1-D calculation is::
  250. avg = sum(a * weights) / sum(weights)
  251. The only constraint on `weights` is that `sum(weights)` must not be 0.
  252. returned : bool, optional
  253. Default is `False`. If `True`, the tuple (`average`, `sum_of_weights`)
  254. is returned, otherwise only the average is returned.
  255. If `weights=None`, `sum_of_weights` is equivalent to the number of
  256. elements over which the average is taken.
  257. Returns
  258. -------
  259. retval, [sum_of_weights] : array_type or double
  260. Return the average along the specified axis. When `returned` is `True`,
  261. return a tuple with the average as the first element and the sum
  262. of the weights as the second element. `sum_of_weights` is of the
  263. same type as `retval`. The result dtype follows a genereal pattern.
  264. If `weights` is None, the result dtype will be that of `a` , or ``float64``
  265. if `a` is integral. Otherwise, if `weights` is not None and `a` is non-
  266. integral, the result type will be the type of lowest precision capable of
  267. representing values of both `a` and `weights`. If `a` happens to be
  268. integral, the previous rules still applies but the result dtype will
  269. at least be ``float64``.
  270. Raises
  271. ------
  272. ZeroDivisionError
  273. When all weights along axis are zero. See `numpy.ma.average` for a
  274. version robust to this type of error.
  275. TypeError
  276. When the length of 1D `weights` is not the same as the shape of `a`
  277. along axis.
  278. See Also
  279. --------
  280. mean
  281. ma.average : average for masked arrays -- useful if your data contains
  282. "missing" values
  283. numpy.result_type : Returns the type that results from applying the
  284. numpy type promotion rules to the arguments.
  285. Examples
  286. --------
  287. >>> data = np.arange(1, 5)
  288. >>> data
  289. array([1, 2, 3, 4])
  290. >>> np.average(data)
  291. 2.5
  292. >>> np.average(np.arange(1, 11), weights=np.arange(10, 0, -1))
  293. 4.0
  294. >>> data = np.arange(6).reshape((3,2))
  295. >>> data
  296. array([[0, 1],
  297. [2, 3],
  298. [4, 5]])
  299. >>> np.average(data, axis=1, weights=[1./4, 3./4])
  300. array([0.75, 2.75, 4.75])
  301. >>> np.average(data, weights=[1./4, 3./4])
  302. Traceback (most recent call last):
  303. ...
  304. TypeError: Axis must be specified when shapes of a and weights differ.
  305. >>> a = np.ones(5, dtype=np.float128)
  306. >>> w = np.ones(5, dtype=np.complex64)
  307. >>> avg = np.average(a, weights=w)
  308. >>> print(avg.dtype)
  309. complex256
  310. """
  311. a = np.asanyarray(a)
  312. if weights is None:
  313. avg = a.mean(axis)
  314. scl = avg.dtype.type(a.size/avg.size)
  315. else:
  316. wgt = np.asanyarray(weights)
  317. if issubclass(a.dtype.type, (np.integer, np.bool_)):
  318. result_dtype = np.result_type(a.dtype, wgt.dtype, 'f8')
  319. else:
  320. result_dtype = np.result_type(a.dtype, wgt.dtype)
  321. # Sanity checks
  322. if a.shape != wgt.shape:
  323. if axis is None:
  324. raise TypeError(
  325. "Axis must be specified when shapes of a and weights "
  326. "differ.")
  327. if wgt.ndim != 1:
  328. raise TypeError(
  329. "1D weights expected when shapes of a and weights differ.")
  330. if wgt.shape[0] != a.shape[axis]:
  331. raise ValueError(
  332. "Length of weights not compatible with specified axis.")
  333. # setup wgt to broadcast along axis
  334. wgt = np.broadcast_to(wgt, (a.ndim-1)*(1,) + wgt.shape)
  335. wgt = wgt.swapaxes(-1, axis)
  336. scl = wgt.sum(axis=axis, dtype=result_dtype)
  337. if np.any(scl == 0.0):
  338. raise ZeroDivisionError(
  339. "Weights sum to zero, can't be normalized")
  340. avg = np.multiply(a, wgt, dtype=result_dtype).sum(axis)/scl
  341. if returned:
  342. if scl.shape != avg.shape:
  343. scl = np.broadcast_to(scl, avg.shape).copy()
  344. return avg, scl
  345. else:
  346. return avg
  347. @set_module('numpy')
  348. def asarray_chkfinite(a, dtype=None, order=None):
  349. """Convert the input to an array, checking for NaNs or Infs.
  350. Parameters
  351. ----------
  352. a : array_like
  353. Input data, in any form that can be converted to an array. This
  354. includes lists, lists of tuples, tuples, tuples of tuples, tuples
  355. of lists and ndarrays. Success requires no NaNs or Infs.
  356. dtype : data-type, optional
  357. By default, the data-type is inferred from the input data.
  358. order : {'C', 'F', 'A', 'K'}, optional
  359. Memory layout. 'A' and 'K' depend on the order of input array a.
  360. 'C' row-major (C-style),
  361. 'F' column-major (Fortran-style) memory representation.
  362. 'A' (any) means 'F' if `a` is Fortran contiguous, 'C' otherwise
  363. 'K' (keep) preserve input order
  364. Defaults to 'C'.
  365. Returns
  366. -------
  367. out : ndarray
  368. Array interpretation of `a`. No copy is performed if the input
  369. is already an ndarray. If `a` is a subclass of ndarray, a base
  370. class ndarray is returned.
  371. Raises
  372. ------
  373. ValueError
  374. Raises ValueError if `a` contains NaN (Not a Number) or Inf (Infinity).
  375. See Also
  376. --------
  377. asarray : Create and array.
  378. asanyarray : Similar function which passes through subclasses.
  379. ascontiguousarray : Convert input to a contiguous array.
  380. asfarray : Convert input to a floating point ndarray.
  381. asfortranarray : Convert input to an ndarray with column-major
  382. memory order.
  383. fromiter : Create an array from an iterator.
  384. fromfunction : Construct an array by executing a function on grid
  385. positions.
  386. Examples
  387. --------
  388. Convert a list into an array. If all elements are finite
  389. ``asarray_chkfinite`` is identical to ``asarray``.
  390. >>> a = [1, 2]
  391. >>> np.asarray_chkfinite(a, dtype=float)
  392. array([1., 2.])
  393. Raises ValueError if array_like contains Nans or Infs.
  394. >>> a = [1, 2, np.inf]
  395. >>> try:
  396. ... np.asarray_chkfinite(a)
  397. ... except ValueError:
  398. ... print('ValueError')
  399. ...
  400. ValueError
  401. """
  402. a = asarray(a, dtype=dtype, order=order)
  403. if a.dtype.char in typecodes['AllFloat'] and not np.isfinite(a).all():
  404. raise ValueError(
  405. "array must not contain infs or NaNs")
  406. return a
  407. def _piecewise_dispatcher(x, condlist, funclist, *args, **kw):
  408. yield x
  409. # support the undocumented behavior of allowing scalars
  410. if np.iterable(condlist):
  411. yield from condlist
  412. @array_function_dispatch(_piecewise_dispatcher)
  413. def piecewise(x, condlist, funclist, *args, **kw):
  414. """
  415. Evaluate a piecewise-defined function.
  416. Given a set of conditions and corresponding functions, evaluate each
  417. function on the input data wherever its condition is true.
  418. Parameters
  419. ----------
  420. x : ndarray or scalar
  421. The input domain.
  422. condlist : list of bool arrays or bool scalars
  423. Each boolean array corresponds to a function in `funclist`. Wherever
  424. `condlist[i]` is True, `funclist[i](x)` is used as the output value.
  425. Each boolean array in `condlist` selects a piece of `x`,
  426. and should therefore be of the same shape as `x`.
  427. The length of `condlist` must correspond to that of `funclist`.
  428. If one extra function is given, i.e. if
  429. ``len(funclist) == len(condlist) + 1``, then that extra function
  430. is the default value, used wherever all conditions are false.
  431. funclist : list of callables, f(x,*args,**kw), or scalars
  432. Each function is evaluated over `x` wherever its corresponding
  433. condition is True. It should take a 1d array as input and give an 1d
  434. array or a scalar value as output. If, instead of a callable,
  435. a scalar is provided then a constant function (``lambda x: scalar``) is
  436. assumed.
  437. args : tuple, optional
  438. Any further arguments given to `piecewise` are passed to the functions
  439. upon execution, i.e., if called ``piecewise(..., ..., 1, 'a')``, then
  440. each function is called as ``f(x, 1, 'a')``.
  441. kw : dict, optional
  442. Keyword arguments used in calling `piecewise` are passed to the
  443. functions upon execution, i.e., if called
  444. ``piecewise(..., ..., alpha=1)``, then each function is called as
  445. ``f(x, alpha=1)``.
  446. Returns
  447. -------
  448. out : ndarray
  449. The output is the same shape and type as x and is found by
  450. calling the functions in `funclist` on the appropriate portions of `x`,
  451. as defined by the boolean arrays in `condlist`. Portions not covered
  452. by any condition have a default value of 0.
  453. See Also
  454. --------
  455. choose, select, where
  456. Notes
  457. -----
  458. This is similar to choose or select, except that functions are
  459. evaluated on elements of `x` that satisfy the corresponding condition from
  460. `condlist`.
  461. The result is::
  462. |--
  463. |funclist[0](x[condlist[0]])
  464. out = |funclist[1](x[condlist[1]])
  465. |...
  466. |funclist[n2](x[condlist[n2]])
  467. |--
  468. Examples
  469. --------
  470. Define the sigma function, which is -1 for ``x < 0`` and +1 for ``x >= 0``.
  471. >>> x = np.linspace(-2.5, 2.5, 6)
  472. >>> np.piecewise(x, [x < 0, x >= 0], [-1, 1])
  473. array([-1., -1., -1., 1., 1., 1.])
  474. Define the absolute value, which is ``-x`` for ``x <0`` and ``x`` for
  475. ``x >= 0``.
  476. >>> np.piecewise(x, [x < 0, x >= 0], [lambda x: -x, lambda x: x])
  477. array([2.5, 1.5, 0.5, 0.5, 1.5, 2.5])
  478. Apply the same function to a scalar value.
  479. >>> y = -2
  480. >>> np.piecewise(y, [y < 0, y >= 0], [lambda x: -x, lambda x: x])
  481. array(2)
  482. """
  483. x = asanyarray(x)
  484. n2 = len(funclist)
  485. # undocumented: single condition is promoted to a list of one condition
  486. if isscalar(condlist) or (
  487. not isinstance(condlist[0], (list, ndarray)) and x.ndim != 0):
  488. condlist = [condlist]
  489. condlist = array(condlist, dtype=bool)
  490. n = len(condlist)
  491. if n == n2 - 1: # compute the "otherwise" condition.
  492. condelse = ~np.any(condlist, axis=0, keepdims=True)
  493. condlist = np.concatenate([condlist, condelse], axis=0)
  494. n += 1
  495. elif n != n2:
  496. raise ValueError(
  497. "with {} condition(s), either {} or {} functions are expected"
  498. .format(n, n, n+1)
  499. )
  500. y = zeros(x.shape, x.dtype)
  501. for cond, func in zip(condlist, funclist):
  502. if not isinstance(func, collections.abc.Callable):
  503. y[cond] = func
  504. else:
  505. vals = x[cond]
  506. if vals.size > 0:
  507. y[cond] = func(vals, *args, **kw)
  508. return y
  509. def _select_dispatcher(condlist, choicelist, default=None):
  510. yield from condlist
  511. yield from choicelist
  512. @array_function_dispatch(_select_dispatcher)
  513. def select(condlist, choicelist, default=0):
  514. """
  515. Return an array drawn from elements in choicelist, depending on conditions.
  516. Parameters
  517. ----------
  518. condlist : list of bool ndarrays
  519. The list of conditions which determine from which array in `choicelist`
  520. the output elements are taken. When multiple conditions are satisfied,
  521. the first one encountered in `condlist` is used.
  522. choicelist : list of ndarrays
  523. The list of arrays from which the output elements are taken. It has
  524. to be of the same length as `condlist`.
  525. default : scalar, optional
  526. The element inserted in `output` when all conditions evaluate to False.
  527. Returns
  528. -------
  529. output : ndarray
  530. The output at position m is the m-th element of the array in
  531. `choicelist` where the m-th element of the corresponding array in
  532. `condlist` is True.
  533. See Also
  534. --------
  535. where : Return elements from one of two arrays depending on condition.
  536. take, choose, compress, diag, diagonal
  537. Examples
  538. --------
  539. >>> x = np.arange(10)
  540. >>> condlist = [x<3, x>5]
  541. >>> choicelist = [x, x**2]
  542. >>> np.select(condlist, choicelist)
  543. array([ 0, 1, 2, ..., 49, 64, 81])
  544. """
  545. # Check the size of condlist and choicelist are the same, or abort.
  546. if len(condlist) != len(choicelist):
  547. raise ValueError(
  548. 'list of cases must be same length as list of conditions')
  549. # Now that the dtype is known, handle the deprecated select([], []) case
  550. if len(condlist) == 0:
  551. raise ValueError("select with an empty condition list is not possible")
  552. choicelist = [np.asarray(choice) for choice in choicelist]
  553. choicelist.append(np.asarray(default))
  554. # need to get the result type before broadcasting for correct scalar
  555. # behaviour
  556. dtype = np.result_type(*choicelist)
  557. # Convert conditions to arrays and broadcast conditions and choices
  558. # as the shape is needed for the result. Doing it separately optimizes
  559. # for example when all choices are scalars.
  560. condlist = np.broadcast_arrays(*condlist)
  561. choicelist = np.broadcast_arrays(*choicelist)
  562. # If cond array is not an ndarray in boolean format or scalar bool, abort.
  563. for i, cond in enumerate(condlist):
  564. if cond.dtype.type is not np.bool_:
  565. raise TypeError(
  566. 'invalid entry {} in condlist: should be boolean ndarray'.format(i))
  567. if choicelist[0].ndim == 0:
  568. # This may be common, so avoid the call.
  569. result_shape = condlist[0].shape
  570. else:
  571. result_shape = np.broadcast_arrays(condlist[0], choicelist[0])[0].shape
  572. result = np.full(result_shape, choicelist[-1], dtype)
  573. # Use np.copyto to burn each choicelist array onto result, using the
  574. # corresponding condlist as a boolean mask. This is done in reverse
  575. # order since the first choice should take precedence.
  576. choicelist = choicelist[-2::-1]
  577. condlist = condlist[::-1]
  578. for choice, cond in zip(choicelist, condlist):
  579. np.copyto(result, choice, where=cond)
  580. return result
  581. def _copy_dispatcher(a, order=None, subok=None):
  582. return (a,)
  583. @array_function_dispatch(_copy_dispatcher)
  584. def copy(a, order='K', subok=False):
  585. """
  586. Return an array copy of the given object.
  587. Parameters
  588. ----------
  589. a : array_like
  590. Input data.
  591. order : {'C', 'F', 'A', 'K'}, optional
  592. Controls the memory layout of the copy. 'C' means C-order,
  593. 'F' means F-order, 'A' means 'F' if `a` is Fortran contiguous,
  594. 'C' otherwise. 'K' means match the layout of `a` as closely
  595. as possible. (Note that this function and :meth:`ndarray.copy` are very
  596. similar, but have different default values for their order=
  597. arguments.)
  598. subok : bool, optional
  599. If True, then sub-classes will be passed-through, otherwise the
  600. returned array will be forced to be a base-class array (defaults to False).
  601. .. versionadded:: 1.19.0
  602. Returns
  603. -------
  604. arr : ndarray
  605. Array interpretation of `a`.
  606. See Also
  607. --------
  608. ndarray.copy : Preferred method for creating an array copy
  609. Notes
  610. -----
  611. This is equivalent to:
  612. >>> np.array(a, copy=True) #doctest: +SKIP
  613. Examples
  614. --------
  615. Create an array x, with a reference y and a copy z:
  616. >>> x = np.array([1, 2, 3])
  617. >>> y = x
  618. >>> z = np.copy(x)
  619. Note that, when we modify x, y changes, but not z:
  620. >>> x[0] = 10
  621. >>> x[0] == y[0]
  622. True
  623. >>> x[0] == z[0]
  624. False
  625. Note that np.copy is a shallow copy and will not copy object
  626. elements within arrays. This is mainly important for arrays
  627. containing Python objects. The new array will contain the
  628. same object which may lead to surprises if that object can
  629. be modified (is mutable):
  630. >>> a = np.array([1, 'm', [2, 3, 4]], dtype=object)
  631. >>> b = np.copy(a)
  632. >>> b[2][0] = 10
  633. >>> a
  634. array([1, 'm', list([10, 3, 4])], dtype=object)
  635. To ensure all elements within an ``object`` array are copied,
  636. use `copy.deepcopy`:
  637. >>> import copy
  638. >>> a = np.array([1, 'm', [2, 3, 4]], dtype=object)
  639. >>> c = copy.deepcopy(a)
  640. >>> c[2][0] = 10
  641. >>> c
  642. array([1, 'm', list([10, 3, 4])], dtype=object)
  643. >>> a
  644. array([1, 'm', list([2, 3, 4])], dtype=object)
  645. """
  646. return array(a, order=order, subok=subok, copy=True)
  647. # Basic operations
  648. def _gradient_dispatcher(f, *varargs, axis=None, edge_order=None):
  649. yield f
  650. yield from varargs
  651. @array_function_dispatch(_gradient_dispatcher)
  652. def gradient(f, *varargs, axis=None, edge_order=1):
  653. """
  654. Return the gradient of an N-dimensional array.
  655. The gradient is computed using second order accurate central differences
  656. in the interior points and either first or second order accurate one-sides
  657. (forward or backwards) differences at the boundaries.
  658. The returned gradient hence has the same shape as the input array.
  659. Parameters
  660. ----------
  661. f : array_like
  662. An N-dimensional array containing samples of a scalar function.
  663. varargs : list of scalar or array, optional
  664. Spacing between f values. Default unitary spacing for all dimensions.
  665. Spacing can be specified using:
  666. 1. single scalar to specify a sample distance for all dimensions.
  667. 2. N scalars to specify a constant sample distance for each dimension.
  668. i.e. `dx`, `dy`, `dz`, ...
  669. 3. N arrays to specify the coordinates of the values along each
  670. dimension of F. The length of the array must match the size of
  671. the corresponding dimension
  672. 4. Any combination of N scalars/arrays with the meaning of 2. and 3.
  673. If `axis` is given, the number of varargs must equal the number of axes.
  674. Default: 1.
  675. edge_order : {1, 2}, optional
  676. Gradient is calculated using N-th order accurate differences
  677. at the boundaries. Default: 1.
  678. .. versionadded:: 1.9.1
  679. axis : None or int or tuple of ints, optional
  680. Gradient is calculated only along the given axis or axes
  681. The default (axis = None) is to calculate the gradient for all the axes
  682. of the input array. axis may be negative, in which case it counts from
  683. the last to the first axis.
  684. .. versionadded:: 1.11.0
  685. Returns
  686. -------
  687. gradient : ndarray or list of ndarray
  688. A set of ndarrays (or a single ndarray if there is only one dimension)
  689. corresponding to the derivatives of f with respect to each dimension.
  690. Each derivative has the same shape as f.
  691. Examples
  692. --------
  693. >>> f = np.array([1, 2, 4, 7, 11, 16], dtype=float)
  694. >>> np.gradient(f)
  695. array([1. , 1.5, 2.5, 3.5, 4.5, 5. ])
  696. >>> np.gradient(f, 2)
  697. array([0.5 , 0.75, 1.25, 1.75, 2.25, 2.5 ])
  698. Spacing can be also specified with an array that represents the coordinates
  699. of the values F along the dimensions.
  700. For instance a uniform spacing:
  701. >>> x = np.arange(f.size)
  702. >>> np.gradient(f, x)
  703. array([1. , 1.5, 2.5, 3.5, 4.5, 5. ])
  704. Or a non uniform one:
  705. >>> x = np.array([0., 1., 1.5, 3.5, 4., 6.], dtype=float)
  706. >>> np.gradient(f, x)
  707. array([1. , 3. , 3.5, 6.7, 6.9, 2.5])
  708. For two dimensional arrays, the return will be two arrays ordered by
  709. axis. In this example the first array stands for the gradient in
  710. rows and the second one in columns direction:
  711. >>> np.gradient(np.array([[1, 2, 6], [3, 4, 5]], dtype=float))
  712. [array([[ 2., 2., -1.],
  713. [ 2., 2., -1.]]), array([[1. , 2.5, 4. ],
  714. [1. , 1. , 1. ]])]
  715. In this example the spacing is also specified:
  716. uniform for axis=0 and non uniform for axis=1
  717. >>> dx = 2.
  718. >>> y = [1., 1.5, 3.5]
  719. >>> np.gradient(np.array([[1, 2, 6], [3, 4, 5]], dtype=float), dx, y)
  720. [array([[ 1. , 1. , -0.5],
  721. [ 1. , 1. , -0.5]]), array([[2. , 2. , 2. ],
  722. [2. , 1.7, 0.5]])]
  723. It is possible to specify how boundaries are treated using `edge_order`
  724. >>> x = np.array([0, 1, 2, 3, 4])
  725. >>> f = x**2
  726. >>> np.gradient(f, edge_order=1)
  727. array([1., 2., 4., 6., 7.])
  728. >>> np.gradient(f, edge_order=2)
  729. array([0., 2., 4., 6., 8.])
  730. The `axis` keyword can be used to specify a subset of axes of which the
  731. gradient is calculated
  732. >>> np.gradient(np.array([[1, 2, 6], [3, 4, 5]], dtype=float), axis=0)
  733. array([[ 2., 2., -1.],
  734. [ 2., 2., -1.]])
  735. Notes
  736. -----
  737. Assuming that :math:`f\\in C^{3}` (i.e., :math:`f` has at least 3 continuous
  738. derivatives) and let :math:`h_{*}` be a non-homogeneous stepsize, we
  739. minimize the "consistency error" :math:`\\eta_{i}` between the true gradient
  740. and its estimate from a linear combination of the neighboring grid-points:
  741. .. math::
  742. \\eta_{i} = f_{i}^{\\left(1\\right)} -
  743. \\left[ \\alpha f\\left(x_{i}\\right) +
  744. \\beta f\\left(x_{i} + h_{d}\\right) +
  745. \\gamma f\\left(x_{i}-h_{s}\\right)
  746. \\right]
  747. By substituting :math:`f(x_{i} + h_{d})` and :math:`f(x_{i} - h_{s})`
  748. with their Taylor series expansion, this translates into solving
  749. the following the linear system:
  750. .. math::
  751. \\left\\{
  752. \\begin{array}{r}
  753. \\alpha+\\beta+\\gamma=0 \\\\
  754. \\beta h_{d}-\\gamma h_{s}=1 \\\\
  755. \\beta h_{d}^{2}+\\gamma h_{s}^{2}=0
  756. \\end{array}
  757. \\right.
  758. The resulting approximation of :math:`f_{i}^{(1)}` is the following:
  759. .. math::
  760. \\hat f_{i}^{(1)} =
  761. \\frac{
  762. h_{s}^{2}f\\left(x_{i} + h_{d}\\right)
  763. + \\left(h_{d}^{2} - h_{s}^{2}\\right)f\\left(x_{i}\\right)
  764. - h_{d}^{2}f\\left(x_{i}-h_{s}\\right)}
  765. { h_{s}h_{d}\\left(h_{d} + h_{s}\\right)}
  766. + \\mathcal{O}\\left(\\frac{h_{d}h_{s}^{2}
  767. + h_{s}h_{d}^{2}}{h_{d}
  768. + h_{s}}\\right)
  769. It is worth noting that if :math:`h_{s}=h_{d}`
  770. (i.e., data are evenly spaced)
  771. we find the standard second order approximation:
  772. .. math::
  773. \\hat f_{i}^{(1)}=
  774. \\frac{f\\left(x_{i+1}\\right) - f\\left(x_{i-1}\\right)}{2h}
  775. + \\mathcal{O}\\left(h^{2}\\right)
  776. With a similar procedure the forward/backward approximations used for
  777. boundaries can be derived.
  778. References
  779. ----------
  780. .. [1] Quarteroni A., Sacco R., Saleri F. (2007) Numerical Mathematics
  781. (Texts in Applied Mathematics). New York: Springer.
  782. .. [2] Durran D. R. (1999) Numerical Methods for Wave Equations
  783. in Geophysical Fluid Dynamics. New York: Springer.
  784. .. [3] Fornberg B. (1988) Generation of Finite Difference Formulas on
  785. Arbitrarily Spaced Grids,
  786. Mathematics of Computation 51, no. 184 : 699-706.
  787. `PDF <http://www.ams.org/journals/mcom/1988-51-184/
  788. S0025-5718-1988-0935077-0/S0025-5718-1988-0935077-0.pdf>`_.
  789. """
  790. f = np.asanyarray(f)
  791. N = f.ndim # number of dimensions
  792. if axis is None:
  793. axes = tuple(range(N))
  794. else:
  795. axes = _nx.normalize_axis_tuple(axis, N)
  796. len_axes = len(axes)
  797. n = len(varargs)
  798. if n == 0:
  799. # no spacing argument - use 1 in all axes
  800. dx = [1.0] * len_axes
  801. elif n == 1 and np.ndim(varargs[0]) == 0:
  802. # single scalar for all axes
  803. dx = varargs * len_axes
  804. elif n == len_axes:
  805. # scalar or 1d array for each axis
  806. dx = list(varargs)
  807. for i, distances in enumerate(dx):
  808. distances = np.asanyarray(distances)
  809. if distances.ndim == 0:
  810. continue
  811. elif distances.ndim != 1:
  812. raise ValueError("distances must be either scalars or 1d")
  813. if len(distances) != f.shape[axes[i]]:
  814. raise ValueError("when 1d, distances must match "
  815. "the length of the corresponding dimension")
  816. if np.issubdtype(distances.dtype, np.integer):
  817. # Convert numpy integer types to float64 to avoid modular
  818. # arithmetic in np.diff(distances).
  819. distances = distances.astype(np.float64)
  820. diffx = np.diff(distances)
  821. # if distances are constant reduce to the scalar case
  822. # since it brings a consistent speedup
  823. if (diffx == diffx[0]).all():
  824. diffx = diffx[0]
  825. dx[i] = diffx
  826. else:
  827. raise TypeError("invalid number of arguments")
  828. if edge_order > 2:
  829. raise ValueError("'edge_order' greater than 2 not supported")
  830. # use central differences on interior and one-sided differences on the
  831. # endpoints. This preserves second order-accuracy over the full domain.
  832. outvals = []
  833. # create slice objects --- initially all are [:, :, ..., :]
  834. slice1 = [slice(None)]*N
  835. slice2 = [slice(None)]*N
  836. slice3 = [slice(None)]*N
  837. slice4 = [slice(None)]*N
  838. otype = f.dtype
  839. if otype.type is np.datetime64:
  840. # the timedelta dtype with the same unit information
  841. otype = np.dtype(otype.name.replace('datetime', 'timedelta'))
  842. # view as timedelta to allow addition
  843. f = f.view(otype)
  844. elif otype.type is np.timedelta64:
  845. pass
  846. elif np.issubdtype(otype, np.inexact):
  847. pass
  848. else:
  849. # All other types convert to floating point.
  850. # First check if f is a numpy integer type; if so, convert f to float64
  851. # to avoid modular arithmetic when computing the changes in f.
  852. if np.issubdtype(otype, np.integer):
  853. f = f.astype(np.float64)
  854. otype = np.float64
  855. for axis, ax_dx in zip(axes, dx):
  856. if f.shape[axis] < edge_order + 1:
  857. raise ValueError(
  858. "Shape of array too small to calculate a numerical gradient, "
  859. "at least (edge_order + 1) elements are required.")
  860. # result allocation
  861. out = np.empty_like(f, dtype=otype)
  862. # spacing for the current axis
  863. uniform_spacing = np.ndim(ax_dx) == 0
  864. # Numerical differentiation: 2nd order interior
  865. slice1[axis] = slice(1, -1)
  866. slice2[axis] = slice(None, -2)
  867. slice3[axis] = slice(1, -1)
  868. slice4[axis] = slice(2, None)
  869. if uniform_spacing:
  870. out[tuple(slice1)] = (f[tuple(slice4)] - f[tuple(slice2)]) / (2. * ax_dx)
  871. else:
  872. dx1 = ax_dx[0:-1]
  873. dx2 = ax_dx[1:]
  874. a = -(dx2)/(dx1 * (dx1 + dx2))
  875. b = (dx2 - dx1) / (dx1 * dx2)
  876. c = dx1 / (dx2 * (dx1 + dx2))
  877. # fix the shape for broadcasting
  878. shape = np.ones(N, dtype=int)
  879. shape[axis] = -1
  880. a.shape = b.shape = c.shape = shape
  881. # 1D equivalent -- out[1:-1] = a * f[:-2] + b * f[1:-1] + c * f[2:]
  882. out[tuple(slice1)] = a * f[tuple(slice2)] + b * f[tuple(slice3)] + c * f[tuple(slice4)]
  883. # Numerical differentiation: 1st order edges
  884. if edge_order == 1:
  885. slice1[axis] = 0
  886. slice2[axis] = 1
  887. slice3[axis] = 0
  888. dx_0 = ax_dx if uniform_spacing else ax_dx[0]
  889. # 1D equivalent -- out[0] = (f[1] - f[0]) / (x[1] - x[0])
  890. out[tuple(slice1)] = (f[tuple(slice2)] - f[tuple(slice3)]) / dx_0
  891. slice1[axis] = -1
  892. slice2[axis] = -1
  893. slice3[axis] = -2
  894. dx_n = ax_dx if uniform_spacing else ax_dx[-1]
  895. # 1D equivalent -- out[-1] = (f[-1] - f[-2]) / (x[-1] - x[-2])
  896. out[tuple(slice1)] = (f[tuple(slice2)] - f[tuple(slice3)]) / dx_n
  897. # Numerical differentiation: 2nd order edges
  898. else:
  899. slice1[axis] = 0
  900. slice2[axis] = 0
  901. slice3[axis] = 1
  902. slice4[axis] = 2
  903. if uniform_spacing:
  904. a = -1.5 / ax_dx
  905. b = 2. / ax_dx
  906. c = -0.5 / ax_dx
  907. else:
  908. dx1 = ax_dx[0]
  909. dx2 = ax_dx[1]
  910. a = -(2. * dx1 + dx2)/(dx1 * (dx1 + dx2))
  911. b = (dx1 + dx2) / (dx1 * dx2)
  912. c = - dx1 / (dx2 * (dx1 + dx2))
  913. # 1D equivalent -- out[0] = a * f[0] + b * f[1] + c * f[2]
  914. out[tuple(slice1)] = a * f[tuple(slice2)] + b * f[tuple(slice3)] + c * f[tuple(slice4)]
  915. slice1[axis] = -1
  916. slice2[axis] = -3
  917. slice3[axis] = -2
  918. slice4[axis] = -1
  919. if uniform_spacing:
  920. a = 0.5 / ax_dx
  921. b = -2. / ax_dx
  922. c = 1.5 / ax_dx
  923. else:
  924. dx1 = ax_dx[-2]
  925. dx2 = ax_dx[-1]
  926. a = (dx2) / (dx1 * (dx1 + dx2))
  927. b = - (dx2 + dx1) / (dx1 * dx2)
  928. c = (2. * dx2 + dx1) / (dx2 * (dx1 + dx2))
  929. # 1D equivalent -- out[-1] = a * f[-3] + b * f[-2] + c * f[-1]
  930. out[tuple(slice1)] = a * f[tuple(slice2)] + b * f[tuple(slice3)] + c * f[tuple(slice4)]
  931. outvals.append(out)
  932. # reset the slice object in this dimension to ":"
  933. slice1[axis] = slice(None)
  934. slice2[axis] = slice(None)
  935. slice3[axis] = slice(None)
  936. slice4[axis] = slice(None)
  937. if len_axes == 1:
  938. return outvals[0]
  939. else:
  940. return outvals
  941. def _diff_dispatcher(a, n=None, axis=None, prepend=None, append=None):
  942. return (a, prepend, append)
  943. @array_function_dispatch(_diff_dispatcher)
  944. def diff(a, n=1, axis=-1, prepend=np._NoValue, append=np._NoValue):
  945. """
  946. Calculate the n-th discrete difference along the given axis.
  947. The first difference is given by ``out[i] = a[i+1] - a[i]`` along
  948. the given axis, higher differences are calculated by using `diff`
  949. recursively.
  950. Parameters
  951. ----------
  952. a : array_like
  953. Input array
  954. n : int, optional
  955. The number of times values are differenced. If zero, the input
  956. is returned as-is.
  957. axis : int, optional
  958. The axis along which the difference is taken, default is the
  959. last axis.
  960. prepend, append : array_like, optional
  961. Values to prepend or append to `a` along axis prior to
  962. performing the difference. Scalar values are expanded to
  963. arrays with length 1 in the direction of axis and the shape
  964. of the input array in along all other axes. Otherwise the
  965. dimension and shape must match `a` except along axis.
  966. .. versionadded:: 1.16.0
  967. Returns
  968. -------
  969. diff : ndarray
  970. The n-th differences. The shape of the output is the same as `a`
  971. except along `axis` where the dimension is smaller by `n`. The
  972. type of the output is the same as the type of the difference
  973. between any two elements of `a`. This is the same as the type of
  974. `a` in most cases. A notable exception is `datetime64`, which
  975. results in a `timedelta64` output array.
  976. See Also
  977. --------
  978. gradient, ediff1d, cumsum
  979. Notes
  980. -----
  981. Type is preserved for boolean arrays, so the result will contain
  982. `False` when consecutive elements are the same and `True` when they
  983. differ.
  984. For unsigned integer arrays, the results will also be unsigned. This
  985. should not be surprising, as the result is consistent with
  986. calculating the difference directly:
  987. >>> u8_arr = np.array([1, 0], dtype=np.uint8)
  988. >>> np.diff(u8_arr)
  989. array([255], dtype=uint8)
  990. >>> u8_arr[1,...] - u8_arr[0,...]
  991. 255
  992. If this is not desirable, then the array should be cast to a larger
  993. integer type first:
  994. >>> i16_arr = u8_arr.astype(np.int16)
  995. >>> np.diff(i16_arr)
  996. array([-1], dtype=int16)
  997. Examples
  998. --------
  999. >>> x = np.array([1, 2, 4, 7, 0])
  1000. >>> np.diff(x)
  1001. array([ 1, 2, 3, -7])
  1002. >>> np.diff(x, n=2)
  1003. array([ 1, 1, -10])
  1004. >>> x = np.array([[1, 3, 6, 10], [0, 5, 6, 8]])
  1005. >>> np.diff(x)
  1006. array([[2, 3, 4],
  1007. [5, 1, 2]])
  1008. >>> np.diff(x, axis=0)
  1009. array([[-1, 2, 0, -2]])
  1010. >>> x = np.arange('1066-10-13', '1066-10-16', dtype=np.datetime64)
  1011. >>> np.diff(x)
  1012. array([1, 1], dtype='timedelta64[D]')
  1013. """
  1014. if n == 0:
  1015. return a
  1016. if n < 0:
  1017. raise ValueError(
  1018. "order must be non-negative but got " + repr(n))
  1019. a = asanyarray(a)
  1020. nd = a.ndim
  1021. if nd == 0:
  1022. raise ValueError("diff requires input that is at least one dimensional")
  1023. axis = normalize_axis_index(axis, nd)
  1024. combined = []
  1025. if prepend is not np._NoValue:
  1026. prepend = np.asanyarray(prepend)
  1027. if prepend.ndim == 0:
  1028. shape = list(a.shape)
  1029. shape[axis] = 1
  1030. prepend = np.broadcast_to(prepend, tuple(shape))
  1031. combined.append(prepend)
  1032. combined.append(a)
  1033. if append is not np._NoValue:
  1034. append = np.asanyarray(append)
  1035. if append.ndim == 0:
  1036. shape = list(a.shape)
  1037. shape[axis] = 1
  1038. append = np.broadcast_to(append, tuple(shape))
  1039. combined.append(append)
  1040. if len(combined) > 1:
  1041. a = np.concatenate(combined, axis)
  1042. slice1 = [slice(None)] * nd
  1043. slice2 = [slice(None)] * nd
  1044. slice1[axis] = slice(1, None)
  1045. slice2[axis] = slice(None, -1)
  1046. slice1 = tuple(slice1)
  1047. slice2 = tuple(slice2)
  1048. op = not_equal if a.dtype == np.bool_ else subtract
  1049. for _ in range(n):
  1050. a = op(a[slice1], a[slice2])
  1051. return a
  1052. def _interp_dispatcher(x, xp, fp, left=None, right=None, period=None):
  1053. return (x, xp, fp)
  1054. @array_function_dispatch(_interp_dispatcher)
  1055. def interp(x, xp, fp, left=None, right=None, period=None):
  1056. """
  1057. One-dimensional linear interpolation.
  1058. Returns the one-dimensional piecewise linear interpolant to a function
  1059. with given discrete data points (`xp`, `fp`), evaluated at `x`.
  1060. Parameters
  1061. ----------
  1062. x : array_like
  1063. The x-coordinates at which to evaluate the interpolated values.
  1064. xp : 1-D sequence of floats
  1065. The x-coordinates of the data points, must be increasing if argument
  1066. `period` is not specified. Otherwise, `xp` is internally sorted after
  1067. normalizing the periodic boundaries with ``xp = xp % period``.
  1068. fp : 1-D sequence of float or complex
  1069. The y-coordinates of the data points, same length as `xp`.
  1070. left : optional float or complex corresponding to fp
  1071. Value to return for `x < xp[0]`, default is `fp[0]`.
  1072. right : optional float or complex corresponding to fp
  1073. Value to return for `x > xp[-1]`, default is `fp[-1]`.
  1074. period : None or float, optional
  1075. A period for the x-coordinates. This parameter allows the proper
  1076. interpolation of angular x-coordinates. Parameters `left` and `right`
  1077. are ignored if `period` is specified.
  1078. .. versionadded:: 1.10.0
  1079. Returns
  1080. -------
  1081. y : float or complex (corresponding to fp) or ndarray
  1082. The interpolated values, same shape as `x`.
  1083. Raises
  1084. ------
  1085. ValueError
  1086. If `xp` and `fp` have different length
  1087. If `xp` or `fp` are not 1-D sequences
  1088. If `period == 0`
  1089. See Also
  1090. --------
  1091. scipy.interpolate
  1092. Notes
  1093. -----
  1094. The x-coordinate sequence is expected to be increasing, but this is not
  1095. explicitly enforced. However, if the sequence `xp` is non-increasing,
  1096. interpolation results are meaningless.
  1097. Note that, since NaN is unsortable, `xp` also cannot contain NaNs.
  1098. A simple check for `xp` being strictly increasing is::
  1099. np.all(np.diff(xp) > 0)
  1100. Examples
  1101. --------
  1102. >>> xp = [1, 2, 3]
  1103. >>> fp = [3, 2, 0]
  1104. >>> np.interp(2.5, xp, fp)
  1105. 1.0
  1106. >>> np.interp([0, 1, 1.5, 2.72, 3.14], xp, fp)
  1107. array([3. , 3. , 2.5 , 0.56, 0. ])
  1108. >>> UNDEF = -99.0
  1109. >>> np.interp(3.14, xp, fp, right=UNDEF)
  1110. -99.0
  1111. Plot an interpolant to the sine function:
  1112. >>> x = np.linspace(0, 2*np.pi, 10)
  1113. >>> y = np.sin(x)
  1114. >>> xvals = np.linspace(0, 2*np.pi, 50)
  1115. >>> yinterp = np.interp(xvals, x, y)
  1116. >>> import matplotlib.pyplot as plt
  1117. >>> plt.plot(x, y, 'o')
  1118. [<matplotlib.lines.Line2D object at 0x...>]
  1119. >>> plt.plot(xvals, yinterp, '-x')
  1120. [<matplotlib.lines.Line2D object at 0x...>]
  1121. >>> plt.show()
  1122. Interpolation with periodic x-coordinates:
  1123. >>> x = [-180, -170, -185, 185, -10, -5, 0, 365]
  1124. >>> xp = [190, -190, 350, -350]
  1125. >>> fp = [5, 10, 3, 4]
  1126. >>> np.interp(x, xp, fp, period=360)
  1127. array([7.5 , 5. , 8.75, 6.25, 3. , 3.25, 3.5 , 3.75])
  1128. Complex interpolation:
  1129. >>> x = [1.5, 4.0]
  1130. >>> xp = [2,3,5]
  1131. >>> fp = [1.0j, 0, 2+3j]
  1132. >>> np.interp(x, xp, fp)
  1133. array([0.+1.j , 1.+1.5j])
  1134. """
  1135. fp = np.asarray(fp)
  1136. if np.iscomplexobj(fp):
  1137. interp_func = compiled_interp_complex
  1138. input_dtype = np.complex128
  1139. else:
  1140. interp_func = compiled_interp
  1141. input_dtype = np.float64
  1142. if period is not None:
  1143. if period == 0:
  1144. raise ValueError("period must be a non-zero value")
  1145. period = abs(period)
  1146. left = None
  1147. right = None
  1148. x = np.asarray(x, dtype=np.float64)
  1149. xp = np.asarray(xp, dtype=np.float64)
  1150. fp = np.asarray(fp, dtype=input_dtype)
  1151. if xp.ndim != 1 or fp.ndim != 1:
  1152. raise ValueError("Data points must be 1-D sequences")
  1153. if xp.shape[0] != fp.shape[0]:
  1154. raise ValueError("fp and xp are not of the same length")
  1155. # normalizing periodic boundaries
  1156. x = x % period
  1157. xp = xp % period
  1158. asort_xp = np.argsort(xp)
  1159. xp = xp[asort_xp]
  1160. fp = fp[asort_xp]
  1161. xp = np.concatenate((xp[-1:]-period, xp, xp[0:1]+period))
  1162. fp = np.concatenate((fp[-1:], fp, fp[0:1]))
  1163. return interp_func(x, xp, fp, left, right)
  1164. def _angle_dispatcher(z, deg=None):
  1165. return (z,)
  1166. @array_function_dispatch(_angle_dispatcher)
  1167. def angle(z, deg=False):
  1168. """
  1169. Return the angle of the complex argument.
  1170. Parameters
  1171. ----------
  1172. z : array_like
  1173. A complex number or sequence of complex numbers.
  1174. deg : bool, optional
  1175. Return angle in degrees if True, radians if False (default).
  1176. Returns
  1177. -------
  1178. angle : ndarray or scalar
  1179. The counterclockwise angle from the positive real axis on the complex
  1180. plane in the range ``(-pi, pi]``, with dtype as numpy.float64.
  1181. .. versionchanged:: 1.16.0
  1182. This function works on subclasses of ndarray like `ma.array`.
  1183. See Also
  1184. --------
  1185. arctan2
  1186. absolute
  1187. Notes
  1188. -----
  1189. Although the angle of the complex number 0 is undefined, ``numpy.angle(0)``
  1190. returns the value 0.
  1191. Examples
  1192. --------
  1193. >>> np.angle([1.0, 1.0j, 1+1j]) # in radians
  1194. array([ 0. , 1.57079633, 0.78539816]) # may vary
  1195. >>> np.angle(1+1j, deg=True) # in degrees
  1196. 45.0
  1197. """
  1198. z = asanyarray(z)
  1199. if issubclass(z.dtype.type, _nx.complexfloating):
  1200. zimag = z.imag
  1201. zreal = z.real
  1202. else:
  1203. zimag = 0
  1204. zreal = z
  1205. a = arctan2(zimag, zreal)
  1206. if deg:
  1207. a *= 180/pi
  1208. return a
  1209. def _unwrap_dispatcher(p, discont=None, axis=None):
  1210. return (p,)
  1211. @array_function_dispatch(_unwrap_dispatcher)
  1212. def unwrap(p, discont=pi, axis=-1):
  1213. """
  1214. Unwrap by changing deltas between values to 2*pi complement.
  1215. Unwrap radian phase `p` by changing absolute jumps greater than
  1216. `discont` to their 2*pi complement along the given axis.
  1217. Parameters
  1218. ----------
  1219. p : array_like
  1220. Input array.
  1221. discont : float, optional
  1222. Maximum discontinuity between values, default is ``pi``.
  1223. axis : int, optional
  1224. Axis along which unwrap will operate, default is the last axis.
  1225. Returns
  1226. -------
  1227. out : ndarray
  1228. Output array.
  1229. See Also
  1230. --------
  1231. rad2deg, deg2rad
  1232. Notes
  1233. -----
  1234. If the discontinuity in `p` is smaller than ``pi``, but larger than
  1235. `discont`, no unwrapping is done because taking the 2*pi complement
  1236. would only make the discontinuity larger.
  1237. Examples
  1238. --------
  1239. >>> phase = np.linspace(0, np.pi, num=5)
  1240. >>> phase[3:] += np.pi
  1241. >>> phase
  1242. array([ 0. , 0.78539816, 1.57079633, 5.49778714, 6.28318531]) # may vary
  1243. >>> np.unwrap(phase)
  1244. array([ 0. , 0.78539816, 1.57079633, -0.78539816, 0. ]) # may vary
  1245. """
  1246. p = asarray(p)
  1247. nd = p.ndim
  1248. dd = diff(p, axis=axis)
  1249. slice1 = [slice(None, None)]*nd # full slices
  1250. slice1[axis] = slice(1, None)
  1251. slice1 = tuple(slice1)
  1252. ddmod = mod(dd + pi, 2*pi) - pi
  1253. _nx.copyto(ddmod, pi, where=(ddmod == -pi) & (dd > 0))
  1254. ph_correct = ddmod - dd
  1255. _nx.copyto(ph_correct, 0, where=abs(dd) < discont)
  1256. up = array(p, copy=True, dtype='d')
  1257. up[slice1] = p[slice1] + ph_correct.cumsum(axis)
  1258. return up
  1259. def _sort_complex(a):
  1260. return (a,)
  1261. @array_function_dispatch(_sort_complex)
  1262. def sort_complex(a):
  1263. """
  1264. Sort a complex array using the real part first, then the imaginary part.
  1265. Parameters
  1266. ----------
  1267. a : array_like
  1268. Input array
  1269. Returns
  1270. -------
  1271. out : complex ndarray
  1272. Always returns a sorted complex array.
  1273. Examples
  1274. --------
  1275. >>> np.sort_complex([5, 3, 6, 2, 1])
  1276. array([1.+0.j, 2.+0.j, 3.+0.j, 5.+0.j, 6.+0.j])
  1277. >>> np.sort_complex([1 + 2j, 2 - 1j, 3 - 2j, 3 - 3j, 3 + 5j])
  1278. array([1.+2.j, 2.-1.j, 3.-3.j, 3.-2.j, 3.+5.j])
  1279. """
  1280. b = array(a, copy=True)
  1281. b.sort()
  1282. if not issubclass(b.dtype.type, _nx.complexfloating):
  1283. if b.dtype.char in 'bhBH':
  1284. return b.astype('F')
  1285. elif b.dtype.char == 'g':
  1286. return b.astype('G')
  1287. else:
  1288. return b.astype('D')
  1289. else:
  1290. return b
  1291. def _trim_zeros(filt, trim=None):
  1292. return (filt,)
  1293. @array_function_dispatch(_trim_zeros)
  1294. def trim_zeros(filt, trim='fb'):
  1295. """
  1296. Trim the leading and/or trailing zeros from a 1-D array or sequence.
  1297. Parameters
  1298. ----------
  1299. filt : 1-D array or sequence
  1300. Input array.
  1301. trim : str, optional
  1302. A string with 'f' representing trim from front and 'b' to trim from
  1303. back. Default is 'fb', trim zeros from both front and back of the
  1304. array.
  1305. Returns
  1306. -------
  1307. trimmed : 1-D array or sequence
  1308. The result of trimming the input. The input data type is preserved.
  1309. Examples
  1310. --------
  1311. >>> a = np.array((0, 0, 0, 1, 2, 3, 0, 2, 1, 0))
  1312. >>> np.trim_zeros(a)
  1313. array([1, 2, 3, 0, 2, 1])
  1314. >>> np.trim_zeros(a, 'b')
  1315. array([0, 0, 0, ..., 0, 2, 1])
  1316. The input data type is preserved, list/tuple in means list/tuple out.
  1317. >>> np.trim_zeros([0, 1, 2, 0])
  1318. [1, 2]
  1319. """
  1320. first = 0
  1321. trim = trim.upper()
  1322. if 'F' in trim:
  1323. for i in filt:
  1324. if i != 0.:
  1325. break
  1326. else:
  1327. first = first + 1
  1328. last = len(filt)
  1329. if 'B' in trim:
  1330. for i in filt[::-1]:
  1331. if i != 0.:
  1332. break
  1333. else:
  1334. last = last - 1
  1335. return filt[first:last]
  1336. def _extract_dispatcher(condition, arr):
  1337. return (condition, arr)
  1338. @array_function_dispatch(_extract_dispatcher)
  1339. def extract(condition, arr):
  1340. """
  1341. Return the elements of an array that satisfy some condition.
  1342. This is equivalent to ``np.compress(ravel(condition), ravel(arr))``. If
  1343. `condition` is boolean ``np.extract`` is equivalent to ``arr[condition]``.
  1344. Note that `place` does the exact opposite of `extract`.
  1345. Parameters
  1346. ----------
  1347. condition : array_like
  1348. An array whose nonzero or True entries indicate the elements of `arr`
  1349. to extract.
  1350. arr : array_like
  1351. Input array of the same size as `condition`.
  1352. Returns
  1353. -------
  1354. extract : ndarray
  1355. Rank 1 array of values from `arr` where `condition` is True.
  1356. See Also
  1357. --------
  1358. take, put, copyto, compress, place
  1359. Examples
  1360. --------
  1361. >>> arr = np.arange(12).reshape((3, 4))
  1362. >>> arr
  1363. array([[ 0, 1, 2, 3],
  1364. [ 4, 5, 6, 7],
  1365. [ 8, 9, 10, 11]])
  1366. >>> condition = np.mod(arr, 3)==0
  1367. >>> condition
  1368. array([[ True, False, False, True],
  1369. [False, False, True, False],
  1370. [False, True, False, False]])
  1371. >>> np.extract(condition, arr)
  1372. array([0, 3, 6, 9])
  1373. If `condition` is boolean:
  1374. >>> arr[condition]
  1375. array([0, 3, 6, 9])
  1376. """
  1377. return _nx.take(ravel(arr), nonzero(ravel(condition))[0])
  1378. def _place_dispatcher(arr, mask, vals):
  1379. return (arr, mask, vals)
  1380. @array_function_dispatch(_place_dispatcher)
  1381. def place(arr, mask, vals):
  1382. """
  1383. Change elements of an array based on conditional and input values.
  1384. Similar to ``np.copyto(arr, vals, where=mask)``, the difference is that
  1385. `place` uses the first N elements of `vals`, where N is the number of
  1386. True values in `mask`, while `copyto` uses the elements where `mask`
  1387. is True.
  1388. Note that `extract` does the exact opposite of `place`.
  1389. Parameters
  1390. ----------
  1391. arr : ndarray
  1392. Array to put data into.
  1393. mask : array_like
  1394. Boolean mask array. Must have the same size as `a`.
  1395. vals : 1-D sequence
  1396. Values to put into `a`. Only the first N elements are used, where
  1397. N is the number of True values in `mask`. If `vals` is smaller
  1398. than N, it will be repeated, and if elements of `a` are to be masked,
  1399. this sequence must be non-empty.
  1400. See Also
  1401. --------
  1402. copyto, put, take, extract
  1403. Examples
  1404. --------
  1405. >>> arr = np.arange(6).reshape(2, 3)
  1406. >>> np.place(arr, arr>2, [44, 55])
  1407. >>> arr
  1408. array([[ 0, 1, 2],
  1409. [44, 55, 44]])
  1410. """
  1411. if not isinstance(arr, np.ndarray):
  1412. raise TypeError("argument 1 must be numpy.ndarray, "
  1413. "not {name}".format(name=type(arr).__name__))
  1414. return _insert(arr, mask, vals)
  1415. def disp(mesg, device=None, linefeed=True):
  1416. """
  1417. Display a message on a device.
  1418. Parameters
  1419. ----------
  1420. mesg : str
  1421. Message to display.
  1422. device : object
  1423. Device to write message. If None, defaults to ``sys.stdout`` which is
  1424. very similar to ``print``. `device` needs to have ``write()`` and
  1425. ``flush()`` methods.
  1426. linefeed : bool, optional
  1427. Option whether to print a line feed or not. Defaults to True.
  1428. Raises
  1429. ------
  1430. AttributeError
  1431. If `device` does not have a ``write()`` or ``flush()`` method.
  1432. Examples
  1433. --------
  1434. Besides ``sys.stdout``, a file-like object can also be used as it has
  1435. both required methods:
  1436. >>> from io import StringIO
  1437. >>> buf = StringIO()
  1438. >>> np.disp(u'"Display" in a file', device=buf)
  1439. >>> buf.getvalue()
  1440. '"Display" in a file\\n'
  1441. """
  1442. if device is None:
  1443. device = sys.stdout
  1444. if linefeed:
  1445. device.write('%s\n' % mesg)
  1446. else:
  1447. device.write('%s' % mesg)
  1448. device.flush()
  1449. return
  1450. # See https://docs.scipy.org/doc/numpy/reference/c-api.generalized-ufuncs.html
  1451. _DIMENSION_NAME = r'\w+'
  1452. _CORE_DIMENSION_LIST = '(?:{0:}(?:,{0:})*)?'.format(_DIMENSION_NAME)
  1453. _ARGUMENT = r'\({}\)'.format(_CORE_DIMENSION_LIST)
  1454. _ARGUMENT_LIST = '{0:}(?:,{0:})*'.format(_ARGUMENT)
  1455. _SIGNATURE = '^{0:}->{0:}$'.format(_ARGUMENT_LIST)
  1456. def _parse_gufunc_signature(signature):
  1457. """
  1458. Parse string signatures for a generalized universal function.
  1459. Arguments
  1460. ---------
  1461. signature : string
  1462. Generalized universal function signature, e.g., ``(m,n),(n,p)->(m,p)``
  1463. for ``np.matmul``.
  1464. Returns
  1465. -------
  1466. Tuple of input and output core dimensions parsed from the signature, each
  1467. of the form List[Tuple[str, ...]].
  1468. """
  1469. if not re.match(_SIGNATURE, signature):
  1470. raise ValueError(
  1471. 'not a valid gufunc signature: {}'.format(signature))
  1472. return tuple([tuple(re.findall(_DIMENSION_NAME, arg))
  1473. for arg in re.findall(_ARGUMENT, arg_list)]
  1474. for arg_list in signature.split('->'))
  1475. def _update_dim_sizes(dim_sizes, arg, core_dims):
  1476. """
  1477. Incrementally check and update core dimension sizes for a single argument.
  1478. Arguments
  1479. ---------
  1480. dim_sizes : Dict[str, int]
  1481. Sizes of existing core dimensions. Will be updated in-place.
  1482. arg : ndarray
  1483. Argument to examine.
  1484. core_dims : Tuple[str, ...]
  1485. Core dimensions for this argument.
  1486. """
  1487. if not core_dims:
  1488. return
  1489. num_core_dims = len(core_dims)
  1490. if arg.ndim < num_core_dims:
  1491. raise ValueError(
  1492. '%d-dimensional argument does not have enough '
  1493. 'dimensions for all core dimensions %r'
  1494. % (arg.ndim, core_dims))
  1495. core_shape = arg.shape[-num_core_dims:]
  1496. for dim, size in zip(core_dims, core_shape):
  1497. if dim in dim_sizes:
  1498. if size != dim_sizes[dim]:
  1499. raise ValueError(
  1500. 'inconsistent size for core dimension %r: %r vs %r'
  1501. % (dim, size, dim_sizes[dim]))
  1502. else:
  1503. dim_sizes[dim] = size
  1504. def _parse_input_dimensions(args, input_core_dims):
  1505. """
  1506. Parse broadcast and core dimensions for vectorize with a signature.
  1507. Arguments
  1508. ---------
  1509. args : Tuple[ndarray, ...]
  1510. Tuple of input arguments to examine.
  1511. input_core_dims : List[Tuple[str, ...]]
  1512. List of core dimensions corresponding to each input.
  1513. Returns
  1514. -------
  1515. broadcast_shape : Tuple[int, ...]
  1516. Common shape to broadcast all non-core dimensions to.
  1517. dim_sizes : Dict[str, int]
  1518. Common sizes for named core dimensions.
  1519. """
  1520. broadcast_args = []
  1521. dim_sizes = {}
  1522. for arg, core_dims in zip(args, input_core_dims):
  1523. _update_dim_sizes(dim_sizes, arg, core_dims)
  1524. ndim = arg.ndim - len(core_dims)
  1525. dummy_array = np.lib.stride_tricks.as_strided(0, arg.shape[:ndim])
  1526. broadcast_args.append(dummy_array)
  1527. broadcast_shape = np.lib.stride_tricks._broadcast_shape(*broadcast_args)
  1528. return broadcast_shape, dim_sizes
  1529. def _calculate_shapes(broadcast_shape, dim_sizes, list_of_core_dims):
  1530. """Helper for calculating broadcast shapes with core dimensions."""
  1531. return [broadcast_shape + tuple(dim_sizes[dim] for dim in core_dims)
  1532. for core_dims in list_of_core_dims]
  1533. def _create_arrays(broadcast_shape, dim_sizes, list_of_core_dims, dtypes):
  1534. """Helper for creating output arrays in vectorize."""
  1535. shapes = _calculate_shapes(broadcast_shape, dim_sizes, list_of_core_dims)
  1536. arrays = tuple(np.empty(shape, dtype=dtype)
  1537. for shape, dtype in zip(shapes, dtypes))
  1538. return arrays
  1539. @set_module('numpy')
  1540. class vectorize:
  1541. """
  1542. vectorize(pyfunc, otypes=None, doc=None, excluded=None, cache=False,
  1543. signature=None)
  1544. Generalized function class.
  1545. Define a vectorized function which takes a nested sequence of objects or
  1546. numpy arrays as inputs and returns a single numpy array or a tuple of numpy
  1547. arrays. The vectorized function evaluates `pyfunc` over successive tuples
  1548. of the input arrays like the python map function, except it uses the
  1549. broadcasting rules of numpy.
  1550. The data type of the output of `vectorized` is determined by calling
  1551. the function with the first element of the input. This can be avoided
  1552. by specifying the `otypes` argument.
  1553. Parameters
  1554. ----------
  1555. pyfunc : callable
  1556. A python function or method.
  1557. otypes : str or list of dtypes, optional
  1558. The output data type. It must be specified as either a string of
  1559. typecode characters or a list of data type specifiers. There should
  1560. be one data type specifier for each output.
  1561. doc : str, optional
  1562. The docstring for the function. If None, the docstring will be the
  1563. ``pyfunc.__doc__``.
  1564. excluded : set, optional
  1565. Set of strings or integers representing the positional or keyword
  1566. arguments for which the function will not be vectorized. These will be
  1567. passed directly to `pyfunc` unmodified.
  1568. .. versionadded:: 1.7.0
  1569. cache : bool, optional
  1570. If `True`, then cache the first function call that determines the number
  1571. of outputs if `otypes` is not provided.
  1572. .. versionadded:: 1.7.0
  1573. signature : string, optional
  1574. Generalized universal function signature, e.g., ``(m,n),(n)->(m)`` for
  1575. vectorized matrix-vector multiplication. If provided, ``pyfunc`` will
  1576. be called with (and expected to return) arrays with shapes given by the
  1577. size of corresponding core dimensions. By default, ``pyfunc`` is
  1578. assumed to take scalars as input and output.
  1579. .. versionadded:: 1.12.0
  1580. Returns
  1581. -------
  1582. vectorized : callable
  1583. Vectorized function.
  1584. See Also
  1585. --------
  1586. frompyfunc : Takes an arbitrary Python function and returns a ufunc
  1587. Notes
  1588. -----
  1589. The `vectorize` function is provided primarily for convenience, not for
  1590. performance. The implementation is essentially a for loop.
  1591. If `otypes` is not specified, then a call to the function with the
  1592. first argument will be used to determine the number of outputs. The
  1593. results of this call will be cached if `cache` is `True` to prevent
  1594. calling the function twice. However, to implement the cache, the
  1595. original function must be wrapped which will slow down subsequent
  1596. calls, so only do this if your function is expensive.
  1597. The new keyword argument interface and `excluded` argument support
  1598. further degrades performance.
  1599. References
  1600. ----------
  1601. .. [1] :doc:`/reference/c-api/generalized-ufuncs`
  1602. Examples
  1603. --------
  1604. >>> def myfunc(a, b):
  1605. ... "Return a-b if a>b, otherwise return a+b"
  1606. ... if a > b:
  1607. ... return a - b
  1608. ... else:
  1609. ... return a + b
  1610. >>> vfunc = np.vectorize(myfunc)
  1611. >>> vfunc([1, 2, 3, 4], 2)
  1612. array([3, 4, 1, 2])
  1613. The docstring is taken from the input function to `vectorize` unless it
  1614. is specified:
  1615. >>> vfunc.__doc__
  1616. 'Return a-b if a>b, otherwise return a+b'
  1617. >>> vfunc = np.vectorize(myfunc, doc='Vectorized `myfunc`')
  1618. >>> vfunc.__doc__
  1619. 'Vectorized `myfunc`'
  1620. The output type is determined by evaluating the first element of the input,
  1621. unless it is specified:
  1622. >>> out = vfunc([1, 2, 3, 4], 2)
  1623. >>> type(out[0])
  1624. <class 'numpy.int64'>
  1625. >>> vfunc = np.vectorize(myfunc, otypes=[float])
  1626. >>> out = vfunc([1, 2, 3, 4], 2)
  1627. >>> type(out[0])
  1628. <class 'numpy.float64'>
  1629. The `excluded` argument can be used to prevent vectorizing over certain
  1630. arguments. This can be useful for array-like arguments of a fixed length
  1631. such as the coefficients for a polynomial as in `polyval`:
  1632. >>> def mypolyval(p, x):
  1633. ... _p = list(p)
  1634. ... res = _p.pop(0)
  1635. ... while _p:
  1636. ... res = res*x + _p.pop(0)
  1637. ... return res
  1638. >>> vpolyval = np.vectorize(mypolyval, excluded=['p'])
  1639. >>> vpolyval(p=[1, 2, 3], x=[0, 1])
  1640. array([3, 6])
  1641. Positional arguments may also be excluded by specifying their position:
  1642. >>> vpolyval.excluded.add(0)
  1643. >>> vpolyval([1, 2, 3], x=[0, 1])
  1644. array([3, 6])
  1645. The `signature` argument allows for vectorizing functions that act on
  1646. non-scalar arrays of fixed length. For example, you can use it for a
  1647. vectorized calculation of Pearson correlation coefficient and its p-value:
  1648. >>> import scipy.stats
  1649. >>> pearsonr = np.vectorize(scipy.stats.pearsonr,
  1650. ... signature='(n),(n)->(),()')
  1651. >>> pearsonr([[0, 1, 2, 3]], [[1, 2, 3, 4], [4, 3, 2, 1]])
  1652. (array([ 1., -1.]), array([ 0., 0.]))
  1653. Or for a vectorized convolution:
  1654. >>> convolve = np.vectorize(np.convolve, signature='(n),(m)->(k)')
  1655. >>> convolve(np.eye(4), [1, 2, 1])
  1656. array([[1., 2., 1., 0., 0., 0.],
  1657. [0., 1., 2., 1., 0., 0.],
  1658. [0., 0., 1., 2., 1., 0.],
  1659. [0., 0., 0., 1., 2., 1.]])
  1660. """
  1661. def __init__(self, pyfunc, otypes=None, doc=None, excluded=None,
  1662. cache=False, signature=None):
  1663. self.pyfunc = pyfunc
  1664. self.cache = cache
  1665. self.signature = signature
  1666. self._ufunc = {} # Caching to improve default performance
  1667. if doc is None:
  1668. self.__doc__ = pyfunc.__doc__
  1669. else:
  1670. self.__doc__ = doc
  1671. if isinstance(otypes, str):
  1672. for char in otypes:
  1673. if char not in typecodes['All']:
  1674. raise ValueError("Invalid otype specified: %s" % (char,))
  1675. elif iterable(otypes):
  1676. otypes = ''.join([_nx.dtype(x).char for x in otypes])
  1677. elif otypes is not None:
  1678. raise ValueError("Invalid otype specification")
  1679. self.otypes = otypes
  1680. # Excluded variable support
  1681. if excluded is None:
  1682. excluded = set()
  1683. self.excluded = set(excluded)
  1684. if signature is not None:
  1685. self._in_and_out_core_dims = _parse_gufunc_signature(signature)
  1686. else:
  1687. self._in_and_out_core_dims = None
  1688. def __call__(self, *args, **kwargs):
  1689. """
  1690. Return arrays with the results of `pyfunc` broadcast (vectorized) over
  1691. `args` and `kwargs` not in `excluded`.
  1692. """
  1693. excluded = self.excluded
  1694. if not kwargs and not excluded:
  1695. func = self.pyfunc
  1696. vargs = args
  1697. else:
  1698. # The wrapper accepts only positional arguments: we use `names` and
  1699. # `inds` to mutate `the_args` and `kwargs` to pass to the original
  1700. # function.
  1701. nargs = len(args)
  1702. names = [_n for _n in kwargs if _n not in excluded]
  1703. inds = [_i for _i in range(nargs) if _i not in excluded]
  1704. the_args = list(args)
  1705. def func(*vargs):
  1706. for _n, _i in enumerate(inds):
  1707. the_args[_i] = vargs[_n]
  1708. kwargs.update(zip(names, vargs[len(inds):]))
  1709. return self.pyfunc(*the_args, **kwargs)
  1710. vargs = [args[_i] for _i in inds]
  1711. vargs.extend([kwargs[_n] for _n in names])
  1712. return self._vectorize_call(func=func, args=vargs)
  1713. def _get_ufunc_and_otypes(self, func, args):
  1714. """Return (ufunc, otypes)."""
  1715. # frompyfunc will fail if args is empty
  1716. if not args:
  1717. raise ValueError('args can not be empty')
  1718. if self.otypes is not None:
  1719. otypes = self.otypes
  1720. # self._ufunc is a dictionary whose keys are the number of
  1721. # arguments (i.e. len(args)) and whose values are ufuncs created
  1722. # by frompyfunc. len(args) can be different for different calls if
  1723. # self.pyfunc has parameters with default values. We only use the
  1724. # cache when func is self.pyfunc, which occurs when the call uses
  1725. # only positional arguments and no arguments are excluded.
  1726. nin = len(args)
  1727. nout = len(self.otypes)
  1728. if func is not self.pyfunc or nin not in self._ufunc:
  1729. ufunc = frompyfunc(func, nin, nout)
  1730. else:
  1731. ufunc = None # We'll get it from self._ufunc
  1732. if func is self.pyfunc:
  1733. ufunc = self._ufunc.setdefault(nin, ufunc)
  1734. else:
  1735. # Get number of outputs and output types by calling the function on
  1736. # the first entries of args. We also cache the result to prevent
  1737. # the subsequent call when the ufunc is evaluated.
  1738. # Assumes that ufunc first evaluates the 0th elements in the input
  1739. # arrays (the input values are not checked to ensure this)
  1740. args = [asarray(arg) for arg in args]
  1741. if builtins.any(arg.size == 0 for arg in args):
  1742. raise ValueError('cannot call `vectorize` on size 0 inputs '
  1743. 'unless `otypes` is set')
  1744. inputs = [arg.flat[0] for arg in args]
  1745. outputs = func(*inputs)
  1746. # Performance note: profiling indicates that -- for simple
  1747. # functions at least -- this wrapping can almost double the
  1748. # execution time.
  1749. # Hence we make it optional.
  1750. if self.cache:
  1751. _cache = [outputs]
  1752. def _func(*vargs):
  1753. if _cache:
  1754. return _cache.pop()
  1755. else:
  1756. return func(*vargs)
  1757. else:
  1758. _func = func
  1759. if isinstance(outputs, tuple):
  1760. nout = len(outputs)
  1761. else:
  1762. nout = 1
  1763. outputs = (outputs,)
  1764. otypes = ''.join([asarray(outputs[_k]).dtype.char
  1765. for _k in range(nout)])
  1766. # Performance note: profiling indicates that creating the ufunc is
  1767. # not a significant cost compared with wrapping so it seems not
  1768. # worth trying to cache this.
  1769. ufunc = frompyfunc(_func, len(args), nout)
  1770. return ufunc, otypes
  1771. def _vectorize_call(self, func, args):
  1772. """Vectorized call to `func` over positional `args`."""
  1773. if self.signature is not None:
  1774. res = self._vectorize_call_with_signature(func, args)
  1775. elif not args:
  1776. res = func()
  1777. else:
  1778. ufunc, otypes = self._get_ufunc_and_otypes(func=func, args=args)
  1779. # Convert args to object arrays first
  1780. inputs = [array(a, copy=False, subok=True, dtype=object)
  1781. for a in args]
  1782. outputs = ufunc(*inputs)
  1783. if ufunc.nout == 1:
  1784. res = array(outputs, copy=False, subok=True, dtype=otypes[0])
  1785. else:
  1786. res = tuple([array(x, copy=False, subok=True, dtype=t)
  1787. for x, t in zip(outputs, otypes)])
  1788. return res
  1789. def _vectorize_call_with_signature(self, func, args):
  1790. """Vectorized call over positional arguments with a signature."""
  1791. input_core_dims, output_core_dims = self._in_and_out_core_dims
  1792. if len(args) != len(input_core_dims):
  1793. raise TypeError('wrong number of positional arguments: '
  1794. 'expected %r, got %r'
  1795. % (len(input_core_dims), len(args)))
  1796. args = tuple(asanyarray(arg) for arg in args)
  1797. broadcast_shape, dim_sizes = _parse_input_dimensions(
  1798. args, input_core_dims)
  1799. input_shapes = _calculate_shapes(broadcast_shape, dim_sizes,
  1800. input_core_dims)
  1801. args = [np.broadcast_to(arg, shape, subok=True)
  1802. for arg, shape in zip(args, input_shapes)]
  1803. outputs = None
  1804. otypes = self.otypes
  1805. nout = len(output_core_dims)
  1806. for index in np.ndindex(*broadcast_shape):
  1807. results = func(*(arg[index] for arg in args))
  1808. n_results = len(results) if isinstance(results, tuple) else 1
  1809. if nout != n_results:
  1810. raise ValueError(
  1811. 'wrong number of outputs from pyfunc: expected %r, got %r'
  1812. % (nout, n_results))
  1813. if nout == 1:
  1814. results = (results,)
  1815. if outputs is None:
  1816. for result, core_dims in zip(results, output_core_dims):
  1817. _update_dim_sizes(dim_sizes, result, core_dims)
  1818. if otypes is None:
  1819. otypes = [asarray(result).dtype for result in results]
  1820. outputs = _create_arrays(broadcast_shape, dim_sizes,
  1821. output_core_dims, otypes)
  1822. for output, result in zip(outputs, results):
  1823. output[index] = result
  1824. if outputs is None:
  1825. # did not call the function even once
  1826. if otypes is None:
  1827. raise ValueError('cannot call `vectorize` on size 0 inputs '
  1828. 'unless `otypes` is set')
  1829. if builtins.any(dim not in dim_sizes
  1830. for dims in output_core_dims
  1831. for dim in dims):
  1832. raise ValueError('cannot call `vectorize` with a signature '
  1833. 'including new output dimensions on size 0 '
  1834. 'inputs')
  1835. outputs = _create_arrays(broadcast_shape, dim_sizes,
  1836. output_core_dims, otypes)
  1837. return outputs[0] if nout == 1 else outputs
  1838. def _cov_dispatcher(m, y=None, rowvar=None, bias=None, ddof=None,
  1839. fweights=None, aweights=None, *, dtype=None):
  1840. return (m, y, fweights, aweights)
  1841. @array_function_dispatch(_cov_dispatcher)
  1842. def cov(m, y=None, rowvar=True, bias=False, ddof=None, fweights=None,
  1843. aweights=None, *, dtype=None):
  1844. """
  1845. Estimate a covariance matrix, given data and weights.
  1846. Covariance indicates the level to which two variables vary together.
  1847. If we examine N-dimensional samples, :math:`X = [x_1, x_2, ... x_N]^T`,
  1848. then the covariance matrix element :math:`C_{ij}` is the covariance of
  1849. :math:`x_i` and :math:`x_j`. The element :math:`C_{ii}` is the variance
  1850. of :math:`x_i`.
  1851. See the notes for an outline of the algorithm.
  1852. Parameters
  1853. ----------
  1854. m : array_like
  1855. A 1-D or 2-D array containing multiple variables and observations.
  1856. Each row of `m` represents a variable, and each column a single
  1857. observation of all those variables. Also see `rowvar` below.
  1858. y : array_like, optional
  1859. An additional set of variables and observations. `y` has the same form
  1860. as that of `m`.
  1861. rowvar : bool, optional
  1862. If `rowvar` is True (default), then each row represents a
  1863. variable, with observations in the columns. Otherwise, the relationship
  1864. is transposed: each column represents a variable, while the rows
  1865. contain observations.
  1866. bias : bool, optional
  1867. Default normalization (False) is by ``(N - 1)``, where ``N`` is the
  1868. number of observations given (unbiased estimate). If `bias` is True,
  1869. then normalization is by ``N``. These values can be overridden by using
  1870. the keyword ``ddof`` in numpy versions >= 1.5.
  1871. ddof : int, optional
  1872. If not ``None`` the default value implied by `bias` is overridden.
  1873. Note that ``ddof=1`` will return the unbiased estimate, even if both
  1874. `fweights` and `aweights` are specified, and ``ddof=0`` will return
  1875. the simple average. See the notes for the details. The default value
  1876. is ``None``.
  1877. .. versionadded:: 1.5
  1878. fweights : array_like, int, optional
  1879. 1-D array of integer frequency weights; the number of times each
  1880. observation vector should be repeated.
  1881. .. versionadded:: 1.10
  1882. aweights : array_like, optional
  1883. 1-D array of observation vector weights. These relative weights are
  1884. typically large for observations considered "important" and smaller for
  1885. observations considered less "important". If ``ddof=0`` the array of
  1886. weights can be used to assign probabilities to observation vectors.
  1887. .. versionadded:: 1.10
  1888. dtype : data-type, optional
  1889. Data-type of the result. By default, the return data-type will have
  1890. at least `numpy.float64` precision.
  1891. .. versionadded:: 1.20
  1892. Returns
  1893. -------
  1894. out : ndarray
  1895. The covariance matrix of the variables.
  1896. See Also
  1897. --------
  1898. corrcoef : Normalized covariance matrix
  1899. Notes
  1900. -----
  1901. Assume that the observations are in the columns of the observation
  1902. array `m` and let ``f = fweights`` and ``a = aweights`` for brevity. The
  1903. steps to compute the weighted covariance are as follows::
  1904. >>> m = np.arange(10, dtype=np.float64)
  1905. >>> f = np.arange(10) * 2
  1906. >>> a = np.arange(10) ** 2.
  1907. >>> ddof = 1
  1908. >>> w = f * a
  1909. >>> v1 = np.sum(w)
  1910. >>> v2 = np.sum(w * a)
  1911. >>> m -= np.sum(m * w, axis=None, keepdims=True) / v1
  1912. >>> cov = np.dot(m * w, m.T) * v1 / (v1**2 - ddof * v2)
  1913. Note that when ``a == 1``, the normalization factor
  1914. ``v1 / (v1**2 - ddof * v2)`` goes over to ``1 / (np.sum(f) - ddof)``
  1915. as it should.
  1916. Examples
  1917. --------
  1918. Consider two variables, :math:`x_0` and :math:`x_1`, which
  1919. correlate perfectly, but in opposite directions:
  1920. >>> x = np.array([[0, 2], [1, 1], [2, 0]]).T
  1921. >>> x
  1922. array([[0, 1, 2],
  1923. [2, 1, 0]])
  1924. Note how :math:`x_0` increases while :math:`x_1` decreases. The covariance
  1925. matrix shows this clearly:
  1926. >>> np.cov(x)
  1927. array([[ 1., -1.],
  1928. [-1., 1.]])
  1929. Note that element :math:`C_{0,1}`, which shows the correlation between
  1930. :math:`x_0` and :math:`x_1`, is negative.
  1931. Further, note how `x` and `y` are combined:
  1932. >>> x = [-2.1, -1, 4.3]
  1933. >>> y = [3, 1.1, 0.12]
  1934. >>> X = np.stack((x, y), axis=0)
  1935. >>> np.cov(X)
  1936. array([[11.71 , -4.286 ], # may vary
  1937. [-4.286 , 2.144133]])
  1938. >>> np.cov(x, y)
  1939. array([[11.71 , -4.286 ], # may vary
  1940. [-4.286 , 2.144133]])
  1941. >>> np.cov(x)
  1942. array(11.71)
  1943. """
  1944. # Check inputs
  1945. if ddof is not None and ddof != int(ddof):
  1946. raise ValueError(
  1947. "ddof must be integer")
  1948. # Handles complex arrays too
  1949. m = np.asarray(m)
  1950. if m.ndim > 2:
  1951. raise ValueError("m has more than 2 dimensions")
  1952. if y is not None:
  1953. y = np.asarray(y)
  1954. if y.ndim > 2:
  1955. raise ValueError("y has more than 2 dimensions")
  1956. if dtype is None:
  1957. if y is None:
  1958. dtype = np.result_type(m, np.float64)
  1959. else:
  1960. dtype = np.result_type(m, y, np.float64)
  1961. X = array(m, ndmin=2, dtype=dtype)
  1962. if not rowvar and X.shape[0] != 1:
  1963. X = X.T
  1964. if X.shape[0] == 0:
  1965. return np.array([]).reshape(0, 0)
  1966. if y is not None:
  1967. y = array(y, copy=False, ndmin=2, dtype=dtype)
  1968. if not rowvar and y.shape[0] != 1:
  1969. y = y.T
  1970. X = np.concatenate((X, y), axis=0)
  1971. if ddof is None:
  1972. if bias == 0:
  1973. ddof = 1
  1974. else:
  1975. ddof = 0
  1976. # Get the product of frequencies and weights
  1977. w = None
  1978. if fweights is not None:
  1979. fweights = np.asarray(fweights, dtype=float)
  1980. if not np.all(fweights == np.around(fweights)):
  1981. raise TypeError(
  1982. "fweights must be integer")
  1983. if fweights.ndim > 1:
  1984. raise RuntimeError(
  1985. "cannot handle multidimensional fweights")
  1986. if fweights.shape[0] != X.shape[1]:
  1987. raise RuntimeError(
  1988. "incompatible numbers of samples and fweights")
  1989. if any(fweights < 0):
  1990. raise ValueError(
  1991. "fweights cannot be negative")
  1992. w = fweights
  1993. if aweights is not None:
  1994. aweights = np.asarray(aweights, dtype=float)
  1995. if aweights.ndim > 1:
  1996. raise RuntimeError(
  1997. "cannot handle multidimensional aweights")
  1998. if aweights.shape[0] != X.shape[1]:
  1999. raise RuntimeError(
  2000. "incompatible numbers of samples and aweights")
  2001. if any(aweights < 0):
  2002. raise ValueError(
  2003. "aweights cannot be negative")
  2004. if w is None:
  2005. w = aweights
  2006. else:
  2007. w *= aweights
  2008. avg, w_sum = average(X, axis=1, weights=w, returned=True)
  2009. w_sum = w_sum[0]
  2010. # Determine the normalization
  2011. if w is None:
  2012. fact = X.shape[1] - ddof
  2013. elif ddof == 0:
  2014. fact = w_sum
  2015. elif aweights is None:
  2016. fact = w_sum - ddof
  2017. else:
  2018. fact = w_sum - ddof*sum(w*aweights)/w_sum
  2019. if fact <= 0:
  2020. warnings.warn("Degrees of freedom <= 0 for slice",
  2021. RuntimeWarning, stacklevel=3)
  2022. fact = 0.0
  2023. X -= avg[:, None]
  2024. if w is None:
  2025. X_T = X.T
  2026. else:
  2027. X_T = (X*w).T
  2028. c = dot(X, X_T.conj())
  2029. c *= np.true_divide(1, fact)
  2030. return c.squeeze()
  2031. def _corrcoef_dispatcher(x, y=None, rowvar=None, bias=None, ddof=None, *,
  2032. dtype=None):
  2033. return (x, y)
  2034. @array_function_dispatch(_corrcoef_dispatcher)
  2035. def corrcoef(x, y=None, rowvar=True, bias=np._NoValue, ddof=np._NoValue, *,
  2036. dtype=None):
  2037. """
  2038. Return Pearson product-moment correlation coefficients.
  2039. Please refer to the documentation for `cov` for more detail. The
  2040. relationship between the correlation coefficient matrix, `R`, and the
  2041. covariance matrix, `C`, is
  2042. .. math:: R_{ij} = \\frac{ C_{ij} } { \\sqrt{ C_{ii} * C_{jj} } }
  2043. The values of `R` are between -1 and 1, inclusive.
  2044. Parameters
  2045. ----------
  2046. x : array_like
  2047. A 1-D or 2-D array containing multiple variables and observations.
  2048. Each row of `x` represents a variable, and each column a single
  2049. observation of all those variables. Also see `rowvar` below.
  2050. y : array_like, optional
  2051. An additional set of variables and observations. `y` has the same
  2052. shape as `x`.
  2053. rowvar : bool, optional
  2054. If `rowvar` is True (default), then each row represents a
  2055. variable, with observations in the columns. Otherwise, the relationship
  2056. is transposed: each column represents a variable, while the rows
  2057. contain observations.
  2058. bias : _NoValue, optional
  2059. Has no effect, do not use.
  2060. .. deprecated:: 1.10.0
  2061. ddof : _NoValue, optional
  2062. Has no effect, do not use.
  2063. .. deprecated:: 1.10.0
  2064. dtype : data-type, optional
  2065. Data-type of the result. By default, the return data-type will have
  2066. at least `numpy.float64` precision.
  2067. .. versionadded:: 1.20
  2068. Returns
  2069. -------
  2070. R : ndarray
  2071. The correlation coefficient matrix of the variables.
  2072. See Also
  2073. --------
  2074. cov : Covariance matrix
  2075. Notes
  2076. -----
  2077. Due to floating point rounding the resulting array may not be Hermitian,
  2078. the diagonal elements may not be 1, and the elements may not satisfy the
  2079. inequality abs(a) <= 1. The real and imaginary parts are clipped to the
  2080. interval [-1, 1] in an attempt to improve on that situation but is not
  2081. much help in the complex case.
  2082. This function accepts but discards arguments `bias` and `ddof`. This is
  2083. for backwards compatibility with previous versions of this function. These
  2084. arguments had no effect on the return values of the function and can be
  2085. safely ignored in this and previous versions of numpy.
  2086. Examples
  2087. --------
  2088. In this example we generate two random arrays, ``xarr`` and ``yarr``, and
  2089. compute the row-wise and column-wise Pearson correlation coefficients,
  2090. ``R``. Since ``rowvar`` is true by default, we first find the row-wise
  2091. Pearson correlation coefficients between the variables of ``xarr``.
  2092. >>> import numpy as np
  2093. >>> rng = np.random.default_rng(seed=42)
  2094. >>> xarr = rng.random((3, 3))
  2095. >>> xarr
  2096. array([[0.77395605, 0.43887844, 0.85859792],
  2097. [0.69736803, 0.09417735, 0.97562235],
  2098. [0.7611397 , 0.78606431, 0.12811363]])
  2099. >>> R1 = np.corrcoef(xarr)
  2100. >>> R1
  2101. array([[ 1. , 0.99256089, -0.68080986],
  2102. [ 0.99256089, 1. , -0.76492172],
  2103. [-0.68080986, -0.76492172, 1. ]])
  2104. If we add another set of variables and observations ``yarr``, we can
  2105. compute the row-wise Pearson correlation coefficients between the
  2106. variables in ``xarr`` and ``yarr``.
  2107. >>> yarr = rng.random((3, 3))
  2108. >>> yarr
  2109. array([[0.45038594, 0.37079802, 0.92676499],
  2110. [0.64386512, 0.82276161, 0.4434142 ],
  2111. [0.22723872, 0.55458479, 0.06381726]])
  2112. >>> R2 = np.corrcoef(xarr, yarr)
  2113. >>> R2
  2114. array([[ 1. , 0.99256089, -0.68080986, 0.75008178, -0.934284 ,
  2115. -0.99004057],
  2116. [ 0.99256089, 1. , -0.76492172, 0.82502011, -0.97074098,
  2117. -0.99981569],
  2118. [-0.68080986, -0.76492172, 1. , -0.99507202, 0.89721355,
  2119. 0.77714685],
  2120. [ 0.75008178, 0.82502011, -0.99507202, 1. , -0.93657855,
  2121. -0.83571711],
  2122. [-0.934284 , -0.97074098, 0.89721355, -0.93657855, 1. ,
  2123. 0.97517215],
  2124. [-0.99004057, -0.99981569, 0.77714685, -0.83571711, 0.97517215,
  2125. 1. ]])
  2126. Finally if we use the option ``rowvar=False``, the columns are now
  2127. being treated as the variables and we will find the column-wise Pearson
  2128. correlation coefficients between variables in ``xarr`` and ``yarr``.
  2129. >>> R3 = np.corrcoef(xarr, yarr, rowvar=False)
  2130. >>> R3
  2131. array([[ 1. , 0.77598074, -0.47458546, -0.75078643, -0.9665554 ,
  2132. 0.22423734],
  2133. [ 0.77598074, 1. , -0.92346708, -0.99923895, -0.58826587,
  2134. -0.44069024],
  2135. [-0.47458546, -0.92346708, 1. , 0.93773029, 0.23297648,
  2136. 0.75137473],
  2137. [-0.75078643, -0.99923895, 0.93773029, 1. , 0.55627469,
  2138. 0.47536961],
  2139. [-0.9665554 , -0.58826587, 0.23297648, 0.55627469, 1. ,
  2140. -0.46666491],
  2141. [ 0.22423734, -0.44069024, 0.75137473, 0.47536961, -0.46666491,
  2142. 1. ]])
  2143. """
  2144. if bias is not np._NoValue or ddof is not np._NoValue:
  2145. # 2015-03-15, 1.10
  2146. warnings.warn('bias and ddof have no effect and are deprecated',
  2147. DeprecationWarning, stacklevel=3)
  2148. c = cov(x, y, rowvar, dtype=dtype)
  2149. try:
  2150. d = diag(c)
  2151. except ValueError:
  2152. # scalar covariance
  2153. # nan if incorrect value (nan, inf, 0), 1 otherwise
  2154. return c / c
  2155. stddev = sqrt(d.real)
  2156. c /= stddev[:, None]
  2157. c /= stddev[None, :]
  2158. # Clip real and imaginary parts to [-1, 1]. This does not guarantee
  2159. # abs(a[i,j]) <= 1 for complex arrays, but is the best we can do without
  2160. # excessive work.
  2161. np.clip(c.real, -1, 1, out=c.real)
  2162. if np.iscomplexobj(c):
  2163. np.clip(c.imag, -1, 1, out=c.imag)
  2164. return c
  2165. @set_module('numpy')
  2166. def blackman(M):
  2167. """
  2168. Return the Blackman window.
  2169. The Blackman window is a taper formed by using the first three
  2170. terms of a summation of cosines. It was designed to have close to the
  2171. minimal leakage possible. It is close to optimal, only slightly worse
  2172. than a Kaiser window.
  2173. Parameters
  2174. ----------
  2175. M : int
  2176. Number of points in the output window. If zero or less, an empty
  2177. array is returned.
  2178. Returns
  2179. -------
  2180. out : ndarray
  2181. The window, with the maximum value normalized to one (the value one
  2182. appears only if the number of samples is odd).
  2183. See Also
  2184. --------
  2185. bartlett, hamming, hanning, kaiser
  2186. Notes
  2187. -----
  2188. The Blackman window is defined as
  2189. .. math:: w(n) = 0.42 - 0.5 \\cos(2\\pi n/M) + 0.08 \\cos(4\\pi n/M)
  2190. Most references to the Blackman window come from the signal processing
  2191. literature, where it is used as one of many windowing functions for
  2192. smoothing values. It is also known as an apodization (which means
  2193. "removing the foot", i.e. smoothing discontinuities at the beginning
  2194. and end of the sampled signal) or tapering function. It is known as a
  2195. "near optimal" tapering function, almost as good (by some measures)
  2196. as the kaiser window.
  2197. References
  2198. ----------
  2199. Blackman, R.B. and Tukey, J.W., (1958) The measurement of power spectra,
  2200. Dover Publications, New York.
  2201. Oppenheim, A.V., and R.W. Schafer. Discrete-Time Signal Processing.
  2202. Upper Saddle River, NJ: Prentice-Hall, 1999, pp. 468-471.
  2203. Examples
  2204. --------
  2205. >>> import matplotlib.pyplot as plt
  2206. >>> np.blackman(12)
  2207. array([-1.38777878e-17, 3.26064346e-02, 1.59903635e-01, # may vary
  2208. 4.14397981e-01, 7.36045180e-01, 9.67046769e-01,
  2209. 9.67046769e-01, 7.36045180e-01, 4.14397981e-01,
  2210. 1.59903635e-01, 3.26064346e-02, -1.38777878e-17])
  2211. Plot the window and the frequency response:
  2212. >>> from numpy.fft import fft, fftshift
  2213. >>> window = np.blackman(51)
  2214. >>> plt.plot(window)
  2215. [<matplotlib.lines.Line2D object at 0x...>]
  2216. >>> plt.title("Blackman window")
  2217. Text(0.5, 1.0, 'Blackman window')
  2218. >>> plt.ylabel("Amplitude")
  2219. Text(0, 0.5, 'Amplitude')
  2220. >>> plt.xlabel("Sample")
  2221. Text(0.5, 0, 'Sample')
  2222. >>> plt.show()
  2223. >>> plt.figure()
  2224. <Figure size 640x480 with 0 Axes>
  2225. >>> A = fft(window, 2048) / 25.5
  2226. >>> mag = np.abs(fftshift(A))
  2227. >>> freq = np.linspace(-0.5, 0.5, len(A))
  2228. >>> with np.errstate(divide='ignore', invalid='ignore'):
  2229. ... response = 20 * np.log10(mag)
  2230. ...
  2231. >>> response = np.clip(response, -100, 100)
  2232. >>> plt.plot(freq, response)
  2233. [<matplotlib.lines.Line2D object at 0x...>]
  2234. >>> plt.title("Frequency response of Blackman window")
  2235. Text(0.5, 1.0, 'Frequency response of Blackman window')
  2236. >>> plt.ylabel("Magnitude [dB]")
  2237. Text(0, 0.5, 'Magnitude [dB]')
  2238. >>> plt.xlabel("Normalized frequency [cycles per sample]")
  2239. Text(0.5, 0, 'Normalized frequency [cycles per sample]')
  2240. >>> _ = plt.axis('tight')
  2241. >>> plt.show()
  2242. """
  2243. if M < 1:
  2244. return array([])
  2245. if M == 1:
  2246. return ones(1, float)
  2247. n = arange(1-M, M, 2)
  2248. return 0.42 + 0.5*cos(pi*n/(M-1)) + 0.08*cos(2.0*pi*n/(M-1))
  2249. @set_module('numpy')
  2250. def bartlett(M):
  2251. """
  2252. Return the Bartlett window.
  2253. The Bartlett window is very similar to a triangular window, except
  2254. that the end points are at zero. It is often used in signal
  2255. processing for tapering a signal, without generating too much
  2256. ripple in the frequency domain.
  2257. Parameters
  2258. ----------
  2259. M : int
  2260. Number of points in the output window. If zero or less, an
  2261. empty array is returned.
  2262. Returns
  2263. -------
  2264. out : array
  2265. The triangular window, with the maximum value normalized to one
  2266. (the value one appears only if the number of samples is odd), with
  2267. the first and last samples equal to zero.
  2268. See Also
  2269. --------
  2270. blackman, hamming, hanning, kaiser
  2271. Notes
  2272. -----
  2273. The Bartlett window is defined as
  2274. .. math:: w(n) = \\frac{2}{M-1} \\left(
  2275. \\frac{M-1}{2} - \\left|n - \\frac{M-1}{2}\\right|
  2276. \\right)
  2277. Most references to the Bartlett window come from the signal
  2278. processing literature, where it is used as one of many windowing
  2279. functions for smoothing values. Note that convolution with this
  2280. window produces linear interpolation. It is also known as an
  2281. apodization (which means"removing the foot", i.e. smoothing
  2282. discontinuities at the beginning and end of the sampled signal) or
  2283. tapering function. The fourier transform of the Bartlett is the product
  2284. of two sinc functions.
  2285. Note the excellent discussion in Kanasewich.
  2286. References
  2287. ----------
  2288. .. [1] M.S. Bartlett, "Periodogram Analysis and Continuous Spectra",
  2289. Biometrika 37, 1-16, 1950.
  2290. .. [2] E.R. Kanasewich, "Time Sequence Analysis in Geophysics",
  2291. The University of Alberta Press, 1975, pp. 109-110.
  2292. .. [3] A.V. Oppenheim and R.W. Schafer, "Discrete-Time Signal
  2293. Processing", Prentice-Hall, 1999, pp. 468-471.
  2294. .. [4] Wikipedia, "Window function",
  2295. https://en.wikipedia.org/wiki/Window_function
  2296. .. [5] W.H. Press, B.P. Flannery, S.A. Teukolsky, and W.T. Vetterling,
  2297. "Numerical Recipes", Cambridge University Press, 1986, page 429.
  2298. Examples
  2299. --------
  2300. >>> import matplotlib.pyplot as plt
  2301. >>> np.bartlett(12)
  2302. array([ 0. , 0.18181818, 0.36363636, 0.54545455, 0.72727273, # may vary
  2303. 0.90909091, 0.90909091, 0.72727273, 0.54545455, 0.36363636,
  2304. 0.18181818, 0. ])
  2305. Plot the window and its frequency response (requires SciPy and matplotlib):
  2306. >>> from numpy.fft import fft, fftshift
  2307. >>> window = np.bartlett(51)
  2308. >>> plt.plot(window)
  2309. [<matplotlib.lines.Line2D object at 0x...>]
  2310. >>> plt.title("Bartlett window")
  2311. Text(0.5, 1.0, 'Bartlett window')
  2312. >>> plt.ylabel("Amplitude")
  2313. Text(0, 0.5, 'Amplitude')
  2314. >>> plt.xlabel("Sample")
  2315. Text(0.5, 0, 'Sample')
  2316. >>> plt.show()
  2317. >>> plt.figure()
  2318. <Figure size 640x480 with 0 Axes>
  2319. >>> A = fft(window, 2048) / 25.5
  2320. >>> mag = np.abs(fftshift(A))
  2321. >>> freq = np.linspace(-0.5, 0.5, len(A))
  2322. >>> with np.errstate(divide='ignore', invalid='ignore'):
  2323. ... response = 20 * np.log10(mag)
  2324. ...
  2325. >>> response = np.clip(response, -100, 100)
  2326. >>> plt.plot(freq, response)
  2327. [<matplotlib.lines.Line2D object at 0x...>]
  2328. >>> plt.title("Frequency response of Bartlett window")
  2329. Text(0.5, 1.0, 'Frequency response of Bartlett window')
  2330. >>> plt.ylabel("Magnitude [dB]")
  2331. Text(0, 0.5, 'Magnitude [dB]')
  2332. >>> plt.xlabel("Normalized frequency [cycles per sample]")
  2333. Text(0.5, 0, 'Normalized frequency [cycles per sample]')
  2334. >>> _ = plt.axis('tight')
  2335. >>> plt.show()
  2336. """
  2337. if M < 1:
  2338. return array([])
  2339. if M == 1:
  2340. return ones(1, float)
  2341. n = arange(1-M, M, 2)
  2342. return where(less_equal(n, 0), 1 + n/(M-1), 1 - n/(M-1))
  2343. @set_module('numpy')
  2344. def hanning(M):
  2345. """
  2346. Return the Hanning window.
  2347. The Hanning window is a taper formed by using a weighted cosine.
  2348. Parameters
  2349. ----------
  2350. M : int
  2351. Number of points in the output window. If zero or less, an
  2352. empty array is returned.
  2353. Returns
  2354. -------
  2355. out : ndarray, shape(M,)
  2356. The window, with the maximum value normalized to one (the value
  2357. one appears only if `M` is odd).
  2358. See Also
  2359. --------
  2360. bartlett, blackman, hamming, kaiser
  2361. Notes
  2362. -----
  2363. The Hanning window is defined as
  2364. .. math:: w(n) = 0.5 - 0.5cos\\left(\\frac{2\\pi{n}}{M-1}\\right)
  2365. \\qquad 0 \\leq n \\leq M-1
  2366. The Hanning was named for Julius von Hann, an Austrian meteorologist.
  2367. It is also known as the Cosine Bell. Some authors prefer that it be
  2368. called a Hann window, to help avoid confusion with the very similar
  2369. Hamming window.
  2370. Most references to the Hanning window come from the signal processing
  2371. literature, where it is used as one of many windowing functions for
  2372. smoothing values. It is also known as an apodization (which means
  2373. "removing the foot", i.e. smoothing discontinuities at the beginning
  2374. and end of the sampled signal) or tapering function.
  2375. References
  2376. ----------
  2377. .. [1] Blackman, R.B. and Tukey, J.W., (1958) The measurement of power
  2378. spectra, Dover Publications, New York.
  2379. .. [2] E.R. Kanasewich, "Time Sequence Analysis in Geophysics",
  2380. The University of Alberta Press, 1975, pp. 106-108.
  2381. .. [3] Wikipedia, "Window function",
  2382. https://en.wikipedia.org/wiki/Window_function
  2383. .. [4] W.H. Press, B.P. Flannery, S.A. Teukolsky, and W.T. Vetterling,
  2384. "Numerical Recipes", Cambridge University Press, 1986, page 425.
  2385. Examples
  2386. --------
  2387. >>> np.hanning(12)
  2388. array([0. , 0.07937323, 0.29229249, 0.57115742, 0.82743037,
  2389. 0.97974649, 0.97974649, 0.82743037, 0.57115742, 0.29229249,
  2390. 0.07937323, 0. ])
  2391. Plot the window and its frequency response:
  2392. >>> import matplotlib.pyplot as plt
  2393. >>> from numpy.fft import fft, fftshift
  2394. >>> window = np.hanning(51)
  2395. >>> plt.plot(window)
  2396. [<matplotlib.lines.Line2D object at 0x...>]
  2397. >>> plt.title("Hann window")
  2398. Text(0.5, 1.0, 'Hann window')
  2399. >>> plt.ylabel("Amplitude")
  2400. Text(0, 0.5, 'Amplitude')
  2401. >>> plt.xlabel("Sample")
  2402. Text(0.5, 0, 'Sample')
  2403. >>> plt.show()
  2404. >>> plt.figure()
  2405. <Figure size 640x480 with 0 Axes>
  2406. >>> A = fft(window, 2048) / 25.5
  2407. >>> mag = np.abs(fftshift(A))
  2408. >>> freq = np.linspace(-0.5, 0.5, len(A))
  2409. >>> with np.errstate(divide='ignore', invalid='ignore'):
  2410. ... response = 20 * np.log10(mag)
  2411. ...
  2412. >>> response = np.clip(response, -100, 100)
  2413. >>> plt.plot(freq, response)
  2414. [<matplotlib.lines.Line2D object at 0x...>]
  2415. >>> plt.title("Frequency response of the Hann window")
  2416. Text(0.5, 1.0, 'Frequency response of the Hann window')
  2417. >>> plt.ylabel("Magnitude [dB]")
  2418. Text(0, 0.5, 'Magnitude [dB]')
  2419. >>> plt.xlabel("Normalized frequency [cycles per sample]")
  2420. Text(0.5, 0, 'Normalized frequency [cycles per sample]')
  2421. >>> plt.axis('tight')
  2422. ...
  2423. >>> plt.show()
  2424. """
  2425. if M < 1:
  2426. return array([])
  2427. if M == 1:
  2428. return ones(1, float)
  2429. n = arange(1-M, M, 2)
  2430. return 0.5 + 0.5*cos(pi*n/(M-1))
  2431. @set_module('numpy')
  2432. def hamming(M):
  2433. """
  2434. Return the Hamming window.
  2435. The Hamming window is a taper formed by using a weighted cosine.
  2436. Parameters
  2437. ----------
  2438. M : int
  2439. Number of points in the output window. If zero or less, an
  2440. empty array is returned.
  2441. Returns
  2442. -------
  2443. out : ndarray
  2444. The window, with the maximum value normalized to one (the value
  2445. one appears only if the number of samples is odd).
  2446. See Also
  2447. --------
  2448. bartlett, blackman, hanning, kaiser
  2449. Notes
  2450. -----
  2451. The Hamming window is defined as
  2452. .. math:: w(n) = 0.54 - 0.46cos\\left(\\frac{2\\pi{n}}{M-1}\\right)
  2453. \\qquad 0 \\leq n \\leq M-1
  2454. The Hamming was named for R. W. Hamming, an associate of J. W. Tukey
  2455. and is described in Blackman and Tukey. It was recommended for
  2456. smoothing the truncated autocovariance function in the time domain.
  2457. Most references to the Hamming window come from the signal processing
  2458. literature, where it is used as one of many windowing functions for
  2459. smoothing values. It is also known as an apodization (which means
  2460. "removing the foot", i.e. smoothing discontinuities at the beginning
  2461. and end of the sampled signal) or tapering function.
  2462. References
  2463. ----------
  2464. .. [1] Blackman, R.B. and Tukey, J.W., (1958) The measurement of power
  2465. spectra, Dover Publications, New York.
  2466. .. [2] E.R. Kanasewich, "Time Sequence Analysis in Geophysics", The
  2467. University of Alberta Press, 1975, pp. 109-110.
  2468. .. [3] Wikipedia, "Window function",
  2469. https://en.wikipedia.org/wiki/Window_function
  2470. .. [4] W.H. Press, B.P. Flannery, S.A. Teukolsky, and W.T. Vetterling,
  2471. "Numerical Recipes", Cambridge University Press, 1986, page 425.
  2472. Examples
  2473. --------
  2474. >>> np.hamming(12)
  2475. array([ 0.08 , 0.15302337, 0.34890909, 0.60546483, 0.84123594, # may vary
  2476. 0.98136677, 0.98136677, 0.84123594, 0.60546483, 0.34890909,
  2477. 0.15302337, 0.08 ])
  2478. Plot the window and the frequency response:
  2479. >>> import matplotlib.pyplot as plt
  2480. >>> from numpy.fft import fft, fftshift
  2481. >>> window = np.hamming(51)
  2482. >>> plt.plot(window)
  2483. [<matplotlib.lines.Line2D object at 0x...>]
  2484. >>> plt.title("Hamming window")
  2485. Text(0.5, 1.0, 'Hamming window')
  2486. >>> plt.ylabel("Amplitude")
  2487. Text(0, 0.5, 'Amplitude')
  2488. >>> plt.xlabel("Sample")
  2489. Text(0.5, 0, 'Sample')
  2490. >>> plt.show()
  2491. >>> plt.figure()
  2492. <Figure size 640x480 with 0 Axes>
  2493. >>> A = fft(window, 2048) / 25.5
  2494. >>> mag = np.abs(fftshift(A))
  2495. >>> freq = np.linspace(-0.5, 0.5, len(A))
  2496. >>> response = 20 * np.log10(mag)
  2497. >>> response = np.clip(response, -100, 100)
  2498. >>> plt.plot(freq, response)
  2499. [<matplotlib.lines.Line2D object at 0x...>]
  2500. >>> plt.title("Frequency response of Hamming window")
  2501. Text(0.5, 1.0, 'Frequency response of Hamming window')
  2502. >>> plt.ylabel("Magnitude [dB]")
  2503. Text(0, 0.5, 'Magnitude [dB]')
  2504. >>> plt.xlabel("Normalized frequency [cycles per sample]")
  2505. Text(0.5, 0, 'Normalized frequency [cycles per sample]')
  2506. >>> plt.axis('tight')
  2507. ...
  2508. >>> plt.show()
  2509. """
  2510. if M < 1:
  2511. return array([])
  2512. if M == 1:
  2513. return ones(1, float)
  2514. n = arange(1-M, M, 2)
  2515. return 0.54 + 0.46*cos(pi*n/(M-1))
  2516. ## Code from cephes for i0
  2517. _i0A = [
  2518. -4.41534164647933937950E-18,
  2519. 3.33079451882223809783E-17,
  2520. -2.43127984654795469359E-16,
  2521. 1.71539128555513303061E-15,
  2522. -1.16853328779934516808E-14,
  2523. 7.67618549860493561688E-14,
  2524. -4.85644678311192946090E-13,
  2525. 2.95505266312963983461E-12,
  2526. -1.72682629144155570723E-11,
  2527. 9.67580903537323691224E-11,
  2528. -5.18979560163526290666E-10,
  2529. 2.65982372468238665035E-9,
  2530. -1.30002500998624804212E-8,
  2531. 6.04699502254191894932E-8,
  2532. -2.67079385394061173391E-7,
  2533. 1.11738753912010371815E-6,
  2534. -4.41673835845875056359E-6,
  2535. 1.64484480707288970893E-5,
  2536. -5.75419501008210370398E-5,
  2537. 1.88502885095841655729E-4,
  2538. -5.76375574538582365885E-4,
  2539. 1.63947561694133579842E-3,
  2540. -4.32430999505057594430E-3,
  2541. 1.05464603945949983183E-2,
  2542. -2.37374148058994688156E-2,
  2543. 4.93052842396707084878E-2,
  2544. -9.49010970480476444210E-2,
  2545. 1.71620901522208775349E-1,
  2546. -3.04682672343198398683E-1,
  2547. 6.76795274409476084995E-1
  2548. ]
  2549. _i0B = [
  2550. -7.23318048787475395456E-18,
  2551. -4.83050448594418207126E-18,
  2552. 4.46562142029675999901E-17,
  2553. 3.46122286769746109310E-17,
  2554. -2.82762398051658348494E-16,
  2555. -3.42548561967721913462E-16,
  2556. 1.77256013305652638360E-15,
  2557. 3.81168066935262242075E-15,
  2558. -9.55484669882830764870E-15,
  2559. -4.15056934728722208663E-14,
  2560. 1.54008621752140982691E-14,
  2561. 3.85277838274214270114E-13,
  2562. 7.18012445138366623367E-13,
  2563. -1.79417853150680611778E-12,
  2564. -1.32158118404477131188E-11,
  2565. -3.14991652796324136454E-11,
  2566. 1.18891471078464383424E-11,
  2567. 4.94060238822496958910E-10,
  2568. 3.39623202570838634515E-9,
  2569. 2.26666899049817806459E-8,
  2570. 2.04891858946906374183E-7,
  2571. 2.89137052083475648297E-6,
  2572. 6.88975834691682398426E-5,
  2573. 3.36911647825569408990E-3,
  2574. 8.04490411014108831608E-1
  2575. ]
  2576. def _chbevl(x, vals):
  2577. b0 = vals[0]
  2578. b1 = 0.0
  2579. for i in range(1, len(vals)):
  2580. b2 = b1
  2581. b1 = b0
  2582. b0 = x*b1 - b2 + vals[i]
  2583. return 0.5*(b0 - b2)
  2584. def _i0_1(x):
  2585. return exp(x) * _chbevl(x/2.0-2, _i0A)
  2586. def _i0_2(x):
  2587. return exp(x) * _chbevl(32.0/x - 2.0, _i0B) / sqrt(x)
  2588. def _i0_dispatcher(x):
  2589. return (x,)
  2590. @array_function_dispatch(_i0_dispatcher)
  2591. def i0(x):
  2592. """
  2593. Modified Bessel function of the first kind, order 0.
  2594. Usually denoted :math:`I_0`.
  2595. Parameters
  2596. ----------
  2597. x : array_like of float
  2598. Argument of the Bessel function.
  2599. Returns
  2600. -------
  2601. out : ndarray, shape = x.shape, dtype = float
  2602. The modified Bessel function evaluated at each of the elements of `x`.
  2603. See Also
  2604. --------
  2605. scipy.special.i0, scipy.special.iv, scipy.special.ive
  2606. Notes
  2607. -----
  2608. The scipy implementation is recommended over this function: it is a
  2609. proper ufunc written in C, and more than an order of magnitude faster.
  2610. We use the algorithm published by Clenshaw [1]_ and referenced by
  2611. Abramowitz and Stegun [2]_, for which the function domain is
  2612. partitioned into the two intervals [0,8] and (8,inf), and Chebyshev
  2613. polynomial expansions are employed in each interval. Relative error on
  2614. the domain [0,30] using IEEE arithmetic is documented [3]_ as having a
  2615. peak of 5.8e-16 with an rms of 1.4e-16 (n = 30000).
  2616. References
  2617. ----------
  2618. .. [1] C. W. Clenshaw, "Chebyshev series for mathematical functions", in
  2619. *National Physical Laboratory Mathematical Tables*, vol. 5, London:
  2620. Her Majesty's Stationery Office, 1962.
  2621. .. [2] M. Abramowitz and I. A. Stegun, *Handbook of Mathematical
  2622. Functions*, 10th printing, New York: Dover, 1964, pp. 379.
  2623. http://www.math.sfu.ca/~cbm/aands/page_379.htm
  2624. .. [3] https://metacpan.org/pod/distribution/Math-Cephes/lib/Math/Cephes.pod#i0:-Modified-Bessel-function-of-order-zero
  2625. Examples
  2626. --------
  2627. >>> np.i0(0.)
  2628. array(1.0)
  2629. >>> np.i0([0, 1, 2, 3])
  2630. array([1. , 1.26606588, 2.2795853 , 4.88079259])
  2631. """
  2632. x = np.asanyarray(x)
  2633. if x.dtype.kind == 'c':
  2634. raise TypeError("i0 not supported for complex values")
  2635. if x.dtype.kind != 'f':
  2636. x = x.astype(float)
  2637. x = np.abs(x)
  2638. return piecewise(x, [x <= 8.0], [_i0_1, _i0_2])
  2639. ## End of cephes code for i0
  2640. @set_module('numpy')
  2641. def kaiser(M, beta):
  2642. """
  2643. Return the Kaiser window.
  2644. The Kaiser window is a taper formed by using a Bessel function.
  2645. Parameters
  2646. ----------
  2647. M : int
  2648. Number of points in the output window. If zero or less, an
  2649. empty array is returned.
  2650. beta : float
  2651. Shape parameter for window.
  2652. Returns
  2653. -------
  2654. out : array
  2655. The window, with the maximum value normalized to one (the value
  2656. one appears only if the number of samples is odd).
  2657. See Also
  2658. --------
  2659. bartlett, blackman, hamming, hanning
  2660. Notes
  2661. -----
  2662. The Kaiser window is defined as
  2663. .. math:: w(n) = I_0\\left( \\beta \\sqrt{1-\\frac{4n^2}{(M-1)^2}}
  2664. \\right)/I_0(\\beta)
  2665. with
  2666. .. math:: \\quad -\\frac{M-1}{2} \\leq n \\leq \\frac{M-1}{2},
  2667. where :math:`I_0` is the modified zeroth-order Bessel function.
  2668. The Kaiser was named for Jim Kaiser, who discovered a simple
  2669. approximation to the DPSS window based on Bessel functions. The Kaiser
  2670. window is a very good approximation to the Digital Prolate Spheroidal
  2671. Sequence, or Slepian window, which is the transform which maximizes the
  2672. energy in the main lobe of the window relative to total energy.
  2673. The Kaiser can approximate many other windows by varying the beta
  2674. parameter.
  2675. ==== =======================
  2676. beta Window shape
  2677. ==== =======================
  2678. 0 Rectangular
  2679. 5 Similar to a Hamming
  2680. 6 Similar to a Hanning
  2681. 8.6 Similar to a Blackman
  2682. ==== =======================
  2683. A beta value of 14 is probably a good starting point. Note that as beta
  2684. gets large, the window narrows, and so the number of samples needs to be
  2685. large enough to sample the increasingly narrow spike, otherwise NaNs will
  2686. get returned.
  2687. Most references to the Kaiser window come from the signal processing
  2688. literature, where it is used as one of many windowing functions for
  2689. smoothing values. It is also known as an apodization (which means
  2690. "removing the foot", i.e. smoothing discontinuities at the beginning
  2691. and end of the sampled signal) or tapering function.
  2692. References
  2693. ----------
  2694. .. [1] J. F. Kaiser, "Digital Filters" - Ch 7 in "Systems analysis by
  2695. digital computer", Editors: F.F. Kuo and J.F. Kaiser, p 218-285.
  2696. John Wiley and Sons, New York, (1966).
  2697. .. [2] E.R. Kanasewich, "Time Sequence Analysis in Geophysics", The
  2698. University of Alberta Press, 1975, pp. 177-178.
  2699. .. [3] Wikipedia, "Window function",
  2700. https://en.wikipedia.org/wiki/Window_function
  2701. Examples
  2702. --------
  2703. >>> import matplotlib.pyplot as plt
  2704. >>> np.kaiser(12, 14)
  2705. array([7.72686684e-06, 3.46009194e-03, 4.65200189e-02, # may vary
  2706. 2.29737120e-01, 5.99885316e-01, 9.45674898e-01,
  2707. 9.45674898e-01, 5.99885316e-01, 2.29737120e-01,
  2708. 4.65200189e-02, 3.46009194e-03, 7.72686684e-06])
  2709. Plot the window and the frequency response:
  2710. >>> from numpy.fft import fft, fftshift
  2711. >>> window = np.kaiser(51, 14)
  2712. >>> plt.plot(window)
  2713. [<matplotlib.lines.Line2D object at 0x...>]
  2714. >>> plt.title("Kaiser window")
  2715. Text(0.5, 1.0, 'Kaiser window')
  2716. >>> plt.ylabel("Amplitude")
  2717. Text(0, 0.5, 'Amplitude')
  2718. >>> plt.xlabel("Sample")
  2719. Text(0.5, 0, 'Sample')
  2720. >>> plt.show()
  2721. >>> plt.figure()
  2722. <Figure size 640x480 with 0 Axes>
  2723. >>> A = fft(window, 2048) / 25.5
  2724. >>> mag = np.abs(fftshift(A))
  2725. >>> freq = np.linspace(-0.5, 0.5, len(A))
  2726. >>> response = 20 * np.log10(mag)
  2727. >>> response = np.clip(response, -100, 100)
  2728. >>> plt.plot(freq, response)
  2729. [<matplotlib.lines.Line2D object at 0x...>]
  2730. >>> plt.title("Frequency response of Kaiser window")
  2731. Text(0.5, 1.0, 'Frequency response of Kaiser window')
  2732. >>> plt.ylabel("Magnitude [dB]")
  2733. Text(0, 0.5, 'Magnitude [dB]')
  2734. >>> plt.xlabel("Normalized frequency [cycles per sample]")
  2735. Text(0.5, 0, 'Normalized frequency [cycles per sample]')
  2736. >>> plt.axis('tight')
  2737. (-0.5, 0.5, -100.0, ...) # may vary
  2738. >>> plt.show()
  2739. """
  2740. if M == 1:
  2741. return np.array([1.])
  2742. n = arange(0, M)
  2743. alpha = (M-1)/2.0
  2744. return i0(beta * sqrt(1-((n-alpha)/alpha)**2.0))/i0(float(beta))
  2745. def _sinc_dispatcher(x):
  2746. return (x,)
  2747. @array_function_dispatch(_sinc_dispatcher)
  2748. def sinc(x):
  2749. r"""
  2750. Return the normalized sinc function.
  2751. The sinc function is :math:`\sin(\pi x)/(\pi x)`.
  2752. .. note::
  2753. Note the normalization factor of ``pi`` used in the definition.
  2754. This is the most commonly used definition in signal processing.
  2755. Use ``sinc(x / np.pi)`` to obtain the unnormalized sinc function
  2756. :math:`\sin(x)/(x)` that is more common in mathematics.
  2757. Parameters
  2758. ----------
  2759. x : ndarray
  2760. Array (possibly multi-dimensional) of values for which to to
  2761. calculate ``sinc(x)``.
  2762. Returns
  2763. -------
  2764. out : ndarray
  2765. ``sinc(x)``, which has the same shape as the input.
  2766. Notes
  2767. -----
  2768. ``sinc(0)`` is the limit value 1.
  2769. The name sinc is short for "sine cardinal" or "sinus cardinalis".
  2770. The sinc function is used in various signal processing applications,
  2771. including in anti-aliasing, in the construction of a Lanczos resampling
  2772. filter, and in interpolation.
  2773. For bandlimited interpolation of discrete-time signals, the ideal
  2774. interpolation kernel is proportional to the sinc function.
  2775. References
  2776. ----------
  2777. .. [1] Weisstein, Eric W. "Sinc Function." From MathWorld--A Wolfram Web
  2778. Resource. http://mathworld.wolfram.com/SincFunction.html
  2779. .. [2] Wikipedia, "Sinc function",
  2780. https://en.wikipedia.org/wiki/Sinc_function
  2781. Examples
  2782. --------
  2783. >>> import matplotlib.pyplot as plt
  2784. >>> x = np.linspace(-4, 4, 41)
  2785. >>> np.sinc(x)
  2786. array([-3.89804309e-17, -4.92362781e-02, -8.40918587e-02, # may vary
  2787. -8.90384387e-02, -5.84680802e-02, 3.89804309e-17,
  2788. 6.68206631e-02, 1.16434881e-01, 1.26137788e-01,
  2789. 8.50444803e-02, -3.89804309e-17, -1.03943254e-01,
  2790. -1.89206682e-01, -2.16236208e-01, -1.55914881e-01,
  2791. 3.89804309e-17, 2.33872321e-01, 5.04551152e-01,
  2792. 7.56826729e-01, 9.35489284e-01, 1.00000000e+00,
  2793. 9.35489284e-01, 7.56826729e-01, 5.04551152e-01,
  2794. 2.33872321e-01, 3.89804309e-17, -1.55914881e-01,
  2795. -2.16236208e-01, -1.89206682e-01, -1.03943254e-01,
  2796. -3.89804309e-17, 8.50444803e-02, 1.26137788e-01,
  2797. 1.16434881e-01, 6.68206631e-02, 3.89804309e-17,
  2798. -5.84680802e-02, -8.90384387e-02, -8.40918587e-02,
  2799. -4.92362781e-02, -3.89804309e-17])
  2800. >>> plt.plot(x, np.sinc(x))
  2801. [<matplotlib.lines.Line2D object at 0x...>]
  2802. >>> plt.title("Sinc Function")
  2803. Text(0.5, 1.0, 'Sinc Function')
  2804. >>> plt.ylabel("Amplitude")
  2805. Text(0, 0.5, 'Amplitude')
  2806. >>> plt.xlabel("X")
  2807. Text(0.5, 0, 'X')
  2808. >>> plt.show()
  2809. """
  2810. x = np.asanyarray(x)
  2811. y = pi * where(x == 0, 1.0e-20, x)
  2812. return sin(y)/y
  2813. def _msort_dispatcher(a):
  2814. return (a,)
  2815. @array_function_dispatch(_msort_dispatcher)
  2816. def msort(a):
  2817. """
  2818. Return a copy of an array sorted along the first axis.
  2819. Parameters
  2820. ----------
  2821. a : array_like
  2822. Array to be sorted.
  2823. Returns
  2824. -------
  2825. sorted_array : ndarray
  2826. Array of the same type and shape as `a`.
  2827. See Also
  2828. --------
  2829. sort
  2830. Notes
  2831. -----
  2832. ``np.msort(a)`` is equivalent to ``np.sort(a, axis=0)``.
  2833. """
  2834. b = array(a, subok=True, copy=True)
  2835. b.sort(0)
  2836. return b
  2837. def _ureduce(a, func, **kwargs):
  2838. """
  2839. Internal Function.
  2840. Call `func` with `a` as first argument swapping the axes to use extended
  2841. axis on functions that don't support it natively.
  2842. Returns result and a.shape with axis dims set to 1.
  2843. Parameters
  2844. ----------
  2845. a : array_like
  2846. Input array or object that can be converted to an array.
  2847. func : callable
  2848. Reduction function capable of receiving a single axis argument.
  2849. It is called with `a` as first argument followed by `kwargs`.
  2850. kwargs : keyword arguments
  2851. additional keyword arguments to pass to `func`.
  2852. Returns
  2853. -------
  2854. result : tuple
  2855. Result of func(a, **kwargs) and a.shape with axis dims set to 1
  2856. which can be used to reshape the result to the same shape a ufunc with
  2857. keepdims=True would produce.
  2858. """
  2859. a = np.asanyarray(a)
  2860. axis = kwargs.get('axis', None)
  2861. if axis is not None:
  2862. keepdim = list(a.shape)
  2863. nd = a.ndim
  2864. axis = _nx.normalize_axis_tuple(axis, nd)
  2865. for ax in axis:
  2866. keepdim[ax] = 1
  2867. if len(axis) == 1:
  2868. kwargs['axis'] = axis[0]
  2869. else:
  2870. keep = set(range(nd)) - set(axis)
  2871. nkeep = len(keep)
  2872. # swap axis that should not be reduced to front
  2873. for i, s in enumerate(sorted(keep)):
  2874. a = a.swapaxes(i, s)
  2875. # merge reduced axis
  2876. a = a.reshape(a.shape[:nkeep] + (-1,))
  2877. kwargs['axis'] = -1
  2878. keepdim = tuple(keepdim)
  2879. else:
  2880. keepdim = (1,) * a.ndim
  2881. r = func(a, **kwargs)
  2882. return r, keepdim
  2883. def _median_dispatcher(
  2884. a, axis=None, out=None, overwrite_input=None, keepdims=None):
  2885. return (a, out)
  2886. @array_function_dispatch(_median_dispatcher)
  2887. def median(a, axis=None, out=None, overwrite_input=False, keepdims=False):
  2888. """
  2889. Compute the median along the specified axis.
  2890. Returns the median of the array elements.
  2891. Parameters
  2892. ----------
  2893. a : array_like
  2894. Input array or object that can be converted to an array.
  2895. axis : {int, sequence of int, None}, optional
  2896. Axis or axes along which the medians are computed. The default
  2897. is to compute the median along a flattened version of the array.
  2898. A sequence of axes is supported since version 1.9.0.
  2899. out : ndarray, optional
  2900. Alternative output array in which to place the result. It must
  2901. have the same shape and buffer length as the expected output,
  2902. but the type (of the output) will be cast if necessary.
  2903. overwrite_input : bool, optional
  2904. If True, then allow use of memory of input array `a` for
  2905. calculations. The input array will be modified by the call to
  2906. `median`. This will save memory when you do not need to preserve
  2907. the contents of the input array. Treat the input as undefined,
  2908. but it will probably be fully or partially sorted. Default is
  2909. False. If `overwrite_input` is ``True`` and `a` is not already an
  2910. `ndarray`, an error will be raised.
  2911. keepdims : bool, optional
  2912. If this is set to True, the axes which are reduced are left
  2913. in the result as dimensions with size one. With this option,
  2914. the result will broadcast correctly against the original `arr`.
  2915. .. versionadded:: 1.9.0
  2916. Returns
  2917. -------
  2918. median : ndarray
  2919. A new array holding the result. If the input contains integers
  2920. or floats smaller than ``float64``, then the output data-type is
  2921. ``np.float64``. Otherwise, the data-type of the output is the
  2922. same as that of the input. If `out` is specified, that array is
  2923. returned instead.
  2924. See Also
  2925. --------
  2926. mean, percentile
  2927. Notes
  2928. -----
  2929. Given a vector ``V`` of length ``N``, the median of ``V`` is the
  2930. middle value of a sorted copy of ``V``, ``V_sorted`` - i
  2931. e., ``V_sorted[(N-1)/2]``, when ``N`` is odd, and the average of the
  2932. two middle values of ``V_sorted`` when ``N`` is even.
  2933. Examples
  2934. --------
  2935. >>> a = np.array([[10, 7, 4], [3, 2, 1]])
  2936. >>> a
  2937. array([[10, 7, 4],
  2938. [ 3, 2, 1]])
  2939. >>> np.median(a)
  2940. 3.5
  2941. >>> np.median(a, axis=0)
  2942. array([6.5, 4.5, 2.5])
  2943. >>> np.median(a, axis=1)
  2944. array([7., 2.])
  2945. >>> m = np.median(a, axis=0)
  2946. >>> out = np.zeros_like(m)
  2947. >>> np.median(a, axis=0, out=m)
  2948. array([6.5, 4.5, 2.5])
  2949. >>> m
  2950. array([6.5, 4.5, 2.5])
  2951. >>> b = a.copy()
  2952. >>> np.median(b, axis=1, overwrite_input=True)
  2953. array([7., 2.])
  2954. >>> assert not np.all(a==b)
  2955. >>> b = a.copy()
  2956. >>> np.median(b, axis=None, overwrite_input=True)
  2957. 3.5
  2958. >>> assert not np.all(a==b)
  2959. """
  2960. r, k = _ureduce(a, func=_median, axis=axis, out=out,
  2961. overwrite_input=overwrite_input)
  2962. if keepdims:
  2963. return r.reshape(k)
  2964. else:
  2965. return r
  2966. def _median(a, axis=None, out=None, overwrite_input=False):
  2967. # can't be reasonably be implemented in terms of percentile as we have to
  2968. # call mean to not break astropy
  2969. a = np.asanyarray(a)
  2970. # Set the partition indexes
  2971. if axis is None:
  2972. sz = a.size
  2973. else:
  2974. sz = a.shape[axis]
  2975. if sz % 2 == 0:
  2976. szh = sz // 2
  2977. kth = [szh - 1, szh]
  2978. else:
  2979. kth = [(sz - 1) // 2]
  2980. # Check if the array contains any nan's
  2981. if np.issubdtype(a.dtype, np.inexact):
  2982. kth.append(-1)
  2983. if overwrite_input:
  2984. if axis is None:
  2985. part = a.ravel()
  2986. part.partition(kth)
  2987. else:
  2988. a.partition(kth, axis=axis)
  2989. part = a
  2990. else:
  2991. part = partition(a, kth, axis=axis)
  2992. if part.shape == ():
  2993. # make 0-D arrays work
  2994. return part.item()
  2995. if axis is None:
  2996. axis = 0
  2997. indexer = [slice(None)] * part.ndim
  2998. index = part.shape[axis] // 2
  2999. if part.shape[axis] % 2 == 1:
  3000. # index with slice to allow mean (below) to work
  3001. indexer[axis] = slice(index, index+1)
  3002. else:
  3003. indexer[axis] = slice(index-1, index+1)
  3004. indexer = tuple(indexer)
  3005. # Check if the array contains any nan's
  3006. if np.issubdtype(a.dtype, np.inexact) and sz > 0:
  3007. # warn and return nans like mean would
  3008. rout = mean(part[indexer], axis=axis, out=out)
  3009. return np.lib.utils._median_nancheck(part, rout, axis, out)
  3010. else:
  3011. # if there are no nans
  3012. # Use mean in odd and even case to coerce data type
  3013. # and check, use out array.
  3014. return mean(part[indexer], axis=axis, out=out)
  3015. def _percentile_dispatcher(a, q, axis=None, out=None, overwrite_input=None,
  3016. interpolation=None, keepdims=None):
  3017. return (a, q, out)
  3018. @array_function_dispatch(_percentile_dispatcher)
  3019. def percentile(a, q, axis=None, out=None,
  3020. overwrite_input=False, interpolation='linear', keepdims=False):
  3021. """
  3022. Compute the q-th percentile of the data along the specified axis.
  3023. Returns the q-th percentile(s) of the array elements.
  3024. Parameters
  3025. ----------
  3026. a : array_like
  3027. Input array or object that can be converted to an array.
  3028. q : array_like of float
  3029. Percentile or sequence of percentiles to compute, which must be between
  3030. 0 and 100 inclusive.
  3031. axis : {int, tuple of int, None}, optional
  3032. Axis or axes along which the percentiles are computed. The
  3033. default is to compute the percentile(s) along a flattened
  3034. version of the array.
  3035. .. versionchanged:: 1.9.0
  3036. A tuple of axes is supported
  3037. out : ndarray, optional
  3038. Alternative output array in which to place the result. It must
  3039. have the same shape and buffer length as the expected output,
  3040. but the type (of the output) will be cast if necessary.
  3041. overwrite_input : bool, optional
  3042. If True, then allow the input array `a` to be modified by intermediate
  3043. calculations, to save memory. In this case, the contents of the input
  3044. `a` after this function completes is undefined.
  3045. interpolation : {'linear', 'lower', 'higher', 'midpoint', 'nearest'}
  3046. This optional parameter specifies the interpolation method to
  3047. use when the desired percentile lies between two data points
  3048. ``i < j``:
  3049. * 'linear': ``i + (j - i) * fraction``, where ``fraction``
  3050. is the fractional part of the index surrounded by ``i``
  3051. and ``j``.
  3052. * 'lower': ``i``.
  3053. * 'higher': ``j``.
  3054. * 'nearest': ``i`` or ``j``, whichever is nearest.
  3055. * 'midpoint': ``(i + j) / 2``.
  3056. .. versionadded:: 1.9.0
  3057. keepdims : bool, optional
  3058. If this is set to True, the axes which are reduced are left in
  3059. the result as dimensions with size one. With this option, the
  3060. result will broadcast correctly against the original array `a`.
  3061. .. versionadded:: 1.9.0
  3062. Returns
  3063. -------
  3064. percentile : scalar or ndarray
  3065. If `q` is a single percentile and `axis=None`, then the result
  3066. is a scalar. If multiple percentiles are given, first axis of
  3067. the result corresponds to the percentiles. The other axes are
  3068. the axes that remain after the reduction of `a`. If the input
  3069. contains integers or floats smaller than ``float64``, the output
  3070. data-type is ``float64``. Otherwise, the output data-type is the
  3071. same as that of the input. If `out` is specified, that array is
  3072. returned instead.
  3073. See Also
  3074. --------
  3075. mean
  3076. median : equivalent to ``percentile(..., 50)``
  3077. nanpercentile
  3078. quantile : equivalent to percentile, except with q in the range [0, 1].
  3079. Notes
  3080. -----
  3081. Given a vector ``V`` of length ``N``, the q-th percentile of
  3082. ``V`` is the value ``q/100`` of the way from the minimum to the
  3083. maximum in a sorted copy of ``V``. The values and distances of
  3084. the two nearest neighbors as well as the `interpolation` parameter
  3085. will determine the percentile if the normalized ranking does not
  3086. match the location of ``q`` exactly. This function is the same as
  3087. the median if ``q=50``, the same as the minimum if ``q=0`` and the
  3088. same as the maximum if ``q=100``.
  3089. Examples
  3090. --------
  3091. >>> a = np.array([[10, 7, 4], [3, 2, 1]])
  3092. >>> a
  3093. array([[10, 7, 4],
  3094. [ 3, 2, 1]])
  3095. >>> np.percentile(a, 50)
  3096. 3.5
  3097. >>> np.percentile(a, 50, axis=0)
  3098. array([6.5, 4.5, 2.5])
  3099. >>> np.percentile(a, 50, axis=1)
  3100. array([7., 2.])
  3101. >>> np.percentile(a, 50, axis=1, keepdims=True)
  3102. array([[7.],
  3103. [2.]])
  3104. >>> m = np.percentile(a, 50, axis=0)
  3105. >>> out = np.zeros_like(m)
  3106. >>> np.percentile(a, 50, axis=0, out=out)
  3107. array([6.5, 4.5, 2.5])
  3108. >>> m
  3109. array([6.5, 4.5, 2.5])
  3110. >>> b = a.copy()
  3111. >>> np.percentile(b, 50, axis=1, overwrite_input=True)
  3112. array([7., 2.])
  3113. >>> assert not np.all(a == b)
  3114. The different types of interpolation can be visualized graphically:
  3115. .. plot::
  3116. import matplotlib.pyplot as plt
  3117. a = np.arange(4)
  3118. p = np.linspace(0, 100, 6001)
  3119. ax = plt.gca()
  3120. lines = [
  3121. ('linear', None),
  3122. ('higher', '--'),
  3123. ('lower', '--'),
  3124. ('nearest', '-.'),
  3125. ('midpoint', '-.'),
  3126. ]
  3127. for interpolation, style in lines:
  3128. ax.plot(
  3129. p, np.percentile(a, p, interpolation=interpolation),
  3130. label=interpolation, linestyle=style)
  3131. ax.set(
  3132. title='Interpolation methods for list: ' + str(a),
  3133. xlabel='Percentile',
  3134. ylabel='List item returned',
  3135. yticks=a)
  3136. ax.legend()
  3137. plt.show()
  3138. """
  3139. q = np.true_divide(q, 100)
  3140. q = asanyarray(q) # undo any decay that the ufunc performed (see gh-13105)
  3141. if not _quantile_is_valid(q):
  3142. raise ValueError("Percentiles must be in the range [0, 100]")
  3143. return _quantile_unchecked(
  3144. a, q, axis, out, overwrite_input, interpolation, keepdims)
  3145. def _quantile_dispatcher(a, q, axis=None, out=None, overwrite_input=None,
  3146. interpolation=None, keepdims=None):
  3147. return (a, q, out)
  3148. @array_function_dispatch(_quantile_dispatcher)
  3149. def quantile(a, q, axis=None, out=None,
  3150. overwrite_input=False, interpolation='linear', keepdims=False):
  3151. """
  3152. Compute the q-th quantile of the data along the specified axis.
  3153. .. versionadded:: 1.15.0
  3154. Parameters
  3155. ----------
  3156. a : array_like
  3157. Input array or object that can be converted to an array.
  3158. q : array_like of float
  3159. Quantile or sequence of quantiles to compute, which must be between
  3160. 0 and 1 inclusive.
  3161. axis : {int, tuple of int, None}, optional
  3162. Axis or axes along which the quantiles are computed. The
  3163. default is to compute the quantile(s) along a flattened
  3164. version of the array.
  3165. out : ndarray, optional
  3166. Alternative output array in which to place the result. It must
  3167. have the same shape and buffer length as the expected output,
  3168. but the type (of the output) will be cast if necessary.
  3169. overwrite_input : bool, optional
  3170. If True, then allow the input array `a` to be modified by intermediate
  3171. calculations, to save memory. In this case, the contents of the input
  3172. `a` after this function completes is undefined.
  3173. interpolation : {'linear', 'lower', 'higher', 'midpoint', 'nearest'}
  3174. This optional parameter specifies the interpolation method to
  3175. use when the desired quantile lies between two data points
  3176. ``i < j``:
  3177. * linear: ``i + (j - i) * fraction``, where ``fraction``
  3178. is the fractional part of the index surrounded by ``i``
  3179. and ``j``.
  3180. * lower: ``i``.
  3181. * higher: ``j``.
  3182. * nearest: ``i`` or ``j``, whichever is nearest.
  3183. * midpoint: ``(i + j) / 2``.
  3184. keepdims : bool, optional
  3185. If this is set to True, the axes which are reduced are left in
  3186. the result as dimensions with size one. With this option, the
  3187. result will broadcast correctly against the original array `a`.
  3188. Returns
  3189. -------
  3190. quantile : scalar or ndarray
  3191. If `q` is a single quantile and `axis=None`, then the result
  3192. is a scalar. If multiple quantiles are given, first axis of
  3193. the result corresponds to the quantiles. The other axes are
  3194. the axes that remain after the reduction of `a`. If the input
  3195. contains integers or floats smaller than ``float64``, the output
  3196. data-type is ``float64``. Otherwise, the output data-type is the
  3197. same as that of the input. If `out` is specified, that array is
  3198. returned instead.
  3199. See Also
  3200. --------
  3201. mean
  3202. percentile : equivalent to quantile, but with q in the range [0, 100].
  3203. median : equivalent to ``quantile(..., 0.5)``
  3204. nanquantile
  3205. Notes
  3206. -----
  3207. Given a vector ``V`` of length ``N``, the q-th quantile of
  3208. ``V`` is the value ``q`` of the way from the minimum to the
  3209. maximum in a sorted copy of ``V``. The values and distances of
  3210. the two nearest neighbors as well as the `interpolation` parameter
  3211. will determine the quantile if the normalized ranking does not
  3212. match the location of ``q`` exactly. This function is the same as
  3213. the median if ``q=0.5``, the same as the minimum if ``q=0.0`` and the
  3214. same as the maximum if ``q=1.0``.
  3215. Examples
  3216. --------
  3217. >>> a = np.array([[10, 7, 4], [3, 2, 1]])
  3218. >>> a
  3219. array([[10, 7, 4],
  3220. [ 3, 2, 1]])
  3221. >>> np.quantile(a, 0.5)
  3222. 3.5
  3223. >>> np.quantile(a, 0.5, axis=0)
  3224. array([6.5, 4.5, 2.5])
  3225. >>> np.quantile(a, 0.5, axis=1)
  3226. array([7., 2.])
  3227. >>> np.quantile(a, 0.5, axis=1, keepdims=True)
  3228. array([[7.],
  3229. [2.]])
  3230. >>> m = np.quantile(a, 0.5, axis=0)
  3231. >>> out = np.zeros_like(m)
  3232. >>> np.quantile(a, 0.5, axis=0, out=out)
  3233. array([6.5, 4.5, 2.5])
  3234. >>> m
  3235. array([6.5, 4.5, 2.5])
  3236. >>> b = a.copy()
  3237. >>> np.quantile(b, 0.5, axis=1, overwrite_input=True)
  3238. array([7., 2.])
  3239. >>> assert not np.all(a == b)
  3240. """
  3241. q = np.asanyarray(q)
  3242. if not _quantile_is_valid(q):
  3243. raise ValueError("Quantiles must be in the range [0, 1]")
  3244. return _quantile_unchecked(
  3245. a, q, axis, out, overwrite_input, interpolation, keepdims)
  3246. def _quantile_unchecked(a, q, axis=None, out=None, overwrite_input=False,
  3247. interpolation='linear', keepdims=False):
  3248. """Assumes that q is in [0, 1], and is an ndarray"""
  3249. r, k = _ureduce(a, func=_quantile_ureduce_func, q=q, axis=axis, out=out,
  3250. overwrite_input=overwrite_input,
  3251. interpolation=interpolation)
  3252. if keepdims:
  3253. return r.reshape(q.shape + k)
  3254. else:
  3255. return r
  3256. def _quantile_is_valid(q):
  3257. # avoid expensive reductions, relevant for arrays with < O(1000) elements
  3258. if q.ndim == 1 and q.size < 10:
  3259. for i in range(q.size):
  3260. if q[i] < 0.0 or q[i] > 1.0:
  3261. return False
  3262. else:
  3263. # faster than any()
  3264. if np.count_nonzero(q < 0.0) or np.count_nonzero(q > 1.0):
  3265. return False
  3266. return True
  3267. def _lerp(a, b, t, out=None):
  3268. """ Linearly interpolate from a to b by a factor of t """
  3269. diff_b_a = subtract(b, a)
  3270. # asanyarray is a stop-gap until gh-13105
  3271. lerp_interpolation = asanyarray(add(a, diff_b_a*t, out=out))
  3272. subtract(b, diff_b_a * (1 - t), out=lerp_interpolation, where=t>=0.5)
  3273. if lerp_interpolation.ndim == 0 and out is None:
  3274. lerp_interpolation = lerp_interpolation[()] # unpack 0d arrays
  3275. return lerp_interpolation
  3276. def _quantile_ureduce_func(a, q, axis=None, out=None, overwrite_input=False,
  3277. interpolation='linear', keepdims=False):
  3278. a = asarray(a)
  3279. # ufuncs cause 0d array results to decay to scalars (see gh-13105), which
  3280. # makes them problematic for __setitem__ and attribute access. As a
  3281. # workaround, we call this on the result of every ufunc on a possibly-0d
  3282. # array.
  3283. not_scalar = np.asanyarray
  3284. # prepare a for partitioning
  3285. if overwrite_input:
  3286. if axis is None:
  3287. ap = a.ravel()
  3288. else:
  3289. ap = a
  3290. else:
  3291. if axis is None:
  3292. ap = a.flatten()
  3293. else:
  3294. ap = a.copy()
  3295. if axis is None:
  3296. axis = 0
  3297. if q.ndim > 2:
  3298. # The code below works fine for nd, but it might not have useful
  3299. # semantics. For now, keep the supported dimensions the same as it was
  3300. # before.
  3301. raise ValueError("q must be a scalar or 1d")
  3302. Nx = ap.shape[axis]
  3303. indices = not_scalar(q * (Nx - 1))
  3304. # round fractional indices according to interpolation method
  3305. if interpolation == 'lower':
  3306. indices = floor(indices).astype(intp)
  3307. elif interpolation == 'higher':
  3308. indices = ceil(indices).astype(intp)
  3309. elif interpolation == 'midpoint':
  3310. indices = 0.5 * (floor(indices) + ceil(indices))
  3311. elif interpolation == 'nearest':
  3312. indices = around(indices).astype(intp)
  3313. elif interpolation == 'linear':
  3314. pass # keep index as fraction and interpolate
  3315. else:
  3316. raise ValueError(
  3317. "interpolation can only be 'linear', 'lower' 'higher', "
  3318. "'midpoint', or 'nearest'")
  3319. # The dimensions of `q` are prepended to the output shape, so we need the
  3320. # axis being sampled from `ap` to be first.
  3321. ap = np.moveaxis(ap, axis, 0)
  3322. del axis
  3323. if np.issubdtype(indices.dtype, np.integer):
  3324. # take the points along axis
  3325. if np.issubdtype(a.dtype, np.inexact):
  3326. # may contain nan, which would sort to the end
  3327. ap.partition(concatenate((indices.ravel(), [-1])), axis=0)
  3328. n = np.isnan(ap[-1])
  3329. else:
  3330. # cannot contain nan
  3331. ap.partition(indices.ravel(), axis=0)
  3332. n = np.array(False, dtype=bool)
  3333. r = take(ap, indices, axis=0, out=out)
  3334. else:
  3335. # weight the points above and below the indices
  3336. indices_below = not_scalar(floor(indices)).astype(intp)
  3337. indices_above = not_scalar(indices_below + 1)
  3338. indices_above[indices_above > Nx - 1] = Nx - 1
  3339. if np.issubdtype(a.dtype, np.inexact):
  3340. # may contain nan, which would sort to the end
  3341. ap.partition(concatenate((
  3342. indices_below.ravel(), indices_above.ravel(), [-1]
  3343. )), axis=0)
  3344. n = np.isnan(ap[-1])
  3345. else:
  3346. # cannot contain nan
  3347. ap.partition(concatenate((
  3348. indices_below.ravel(), indices_above.ravel()
  3349. )), axis=0)
  3350. n = np.array(False, dtype=bool)
  3351. weights_shape = indices.shape + (1,) * (ap.ndim - 1)
  3352. weights_above = not_scalar(indices - indices_below).reshape(weights_shape)
  3353. x_below = take(ap, indices_below, axis=0)
  3354. x_above = take(ap, indices_above, axis=0)
  3355. r = _lerp(x_below, x_above, weights_above, out=out)
  3356. # if any slice contained a nan, then all results on that slice are also nan
  3357. if np.any(n):
  3358. if r.ndim == 0 and out is None:
  3359. # can't write to a scalar
  3360. r = a.dtype.type(np.nan)
  3361. else:
  3362. r[..., n] = a.dtype.type(np.nan)
  3363. return r
  3364. def _trapz_dispatcher(y, x=None, dx=None, axis=None):
  3365. return (y, x)
  3366. @array_function_dispatch(_trapz_dispatcher)
  3367. def trapz(y, x=None, dx=1.0, axis=-1):
  3368. """
  3369. Integrate along the given axis using the composite trapezoidal rule.
  3370. Integrate `y` (`x`) along given axis.
  3371. Parameters
  3372. ----------
  3373. y : array_like
  3374. Input array to integrate.
  3375. x : array_like, optional
  3376. The sample points corresponding to the `y` values. If `x` is None,
  3377. the sample points are assumed to be evenly spaced `dx` apart. The
  3378. default is None.
  3379. dx : scalar, optional
  3380. The spacing between sample points when `x` is None. The default is 1.
  3381. axis : int, optional
  3382. The axis along which to integrate.
  3383. Returns
  3384. -------
  3385. trapz : float
  3386. Definite integral as approximated by trapezoidal rule.
  3387. See Also
  3388. --------
  3389. sum, cumsum
  3390. Notes
  3391. -----
  3392. Image [2]_ illustrates trapezoidal rule -- y-axis locations of points
  3393. will be taken from `y` array, by default x-axis distances between
  3394. points will be 1.0, alternatively they can be provided with `x` array
  3395. or with `dx` scalar. Return value will be equal to combined area under
  3396. the red lines.
  3397. References
  3398. ----------
  3399. .. [1] Wikipedia page: https://en.wikipedia.org/wiki/Trapezoidal_rule
  3400. .. [2] Illustration image:
  3401. https://en.wikipedia.org/wiki/File:Composite_trapezoidal_rule_illustration.png
  3402. Examples
  3403. --------
  3404. >>> np.trapz([1,2,3])
  3405. 4.0
  3406. >>> np.trapz([1,2,3], x=[4,6,8])
  3407. 8.0
  3408. >>> np.trapz([1,2,3], dx=2)
  3409. 8.0
  3410. >>> a = np.arange(6).reshape(2, 3)
  3411. >>> a
  3412. array([[0, 1, 2],
  3413. [3, 4, 5]])
  3414. >>> np.trapz(a, axis=0)
  3415. array([1.5, 2.5, 3.5])
  3416. >>> np.trapz(a, axis=1)
  3417. array([2., 8.])
  3418. """
  3419. y = asanyarray(y)
  3420. if x is None:
  3421. d = dx
  3422. else:
  3423. x = asanyarray(x)
  3424. if x.ndim == 1:
  3425. d = diff(x)
  3426. # reshape to correct shape
  3427. shape = [1]*y.ndim
  3428. shape[axis] = d.shape[0]
  3429. d = d.reshape(shape)
  3430. else:
  3431. d = diff(x, axis=axis)
  3432. nd = y.ndim
  3433. slice1 = [slice(None)]*nd
  3434. slice2 = [slice(None)]*nd
  3435. slice1[axis] = slice(1, None)
  3436. slice2[axis] = slice(None, -1)
  3437. try:
  3438. ret = (d * (y[tuple(slice1)] + y[tuple(slice2)]) / 2.0).sum(axis)
  3439. except ValueError:
  3440. # Operations didn't work, cast to ndarray
  3441. d = np.asarray(d)
  3442. y = np.asarray(y)
  3443. ret = add.reduce(d * (y[tuple(slice1)]+y[tuple(slice2)])/2.0, axis)
  3444. return ret
  3445. def _meshgrid_dispatcher(*xi, copy=None, sparse=None, indexing=None):
  3446. return xi
  3447. # Based on scitools meshgrid
  3448. @array_function_dispatch(_meshgrid_dispatcher)
  3449. def meshgrid(*xi, copy=True, sparse=False, indexing='xy'):
  3450. """
  3451. Return coordinate matrices from coordinate vectors.
  3452. Make N-D coordinate arrays for vectorized evaluations of
  3453. N-D scalar/vector fields over N-D grids, given
  3454. one-dimensional coordinate arrays x1, x2,..., xn.
  3455. .. versionchanged:: 1.9
  3456. 1-D and 0-D cases are allowed.
  3457. Parameters
  3458. ----------
  3459. x1, x2,..., xn : array_like
  3460. 1-D arrays representing the coordinates of a grid.
  3461. indexing : {'xy', 'ij'}, optional
  3462. Cartesian ('xy', default) or matrix ('ij') indexing of output.
  3463. See Notes for more details.
  3464. .. versionadded:: 1.7.0
  3465. sparse : bool, optional
  3466. If True a sparse grid is returned in order to conserve memory.
  3467. Default is False.
  3468. .. versionadded:: 1.7.0
  3469. copy : bool, optional
  3470. If False, a view into the original arrays are returned in order to
  3471. conserve memory. Default is True. Please note that
  3472. ``sparse=False, copy=False`` will likely return non-contiguous
  3473. arrays. Furthermore, more than one element of a broadcast array
  3474. may refer to a single memory location. If you need to write to the
  3475. arrays, make copies first.
  3476. .. versionadded:: 1.7.0
  3477. Returns
  3478. -------
  3479. X1, X2,..., XN : ndarray
  3480. For vectors `x1`, `x2`,..., 'xn' with lengths ``Ni=len(xi)`` ,
  3481. return ``(N1, N2, N3,...Nn)`` shaped arrays if indexing='ij'
  3482. or ``(N2, N1, N3,...Nn)`` shaped arrays if indexing='xy'
  3483. with the elements of `xi` repeated to fill the matrix along
  3484. the first dimension for `x1`, the second for `x2` and so on.
  3485. Notes
  3486. -----
  3487. This function supports both indexing conventions through the indexing
  3488. keyword argument. Giving the string 'ij' returns a meshgrid with
  3489. matrix indexing, while 'xy' returns a meshgrid with Cartesian indexing.
  3490. In the 2-D case with inputs of length M and N, the outputs are of shape
  3491. (N, M) for 'xy' indexing and (M, N) for 'ij' indexing. In the 3-D case
  3492. with inputs of length M, N and P, outputs are of shape (N, M, P) for
  3493. 'xy' indexing and (M, N, P) for 'ij' indexing. The difference is
  3494. illustrated by the following code snippet::
  3495. xv, yv = np.meshgrid(x, y, sparse=False, indexing='ij')
  3496. for i in range(nx):
  3497. for j in range(ny):
  3498. # treat xv[i,j], yv[i,j]
  3499. xv, yv = np.meshgrid(x, y, sparse=False, indexing='xy')
  3500. for i in range(nx):
  3501. for j in range(ny):
  3502. # treat xv[j,i], yv[j,i]
  3503. In the 1-D and 0-D case, the indexing and sparse keywords have no effect.
  3504. See Also
  3505. --------
  3506. mgrid : Construct a multi-dimensional "meshgrid" using indexing notation.
  3507. ogrid : Construct an open multi-dimensional "meshgrid" using indexing
  3508. notation.
  3509. Examples
  3510. --------
  3511. >>> nx, ny = (3, 2)
  3512. >>> x = np.linspace(0, 1, nx)
  3513. >>> y = np.linspace(0, 1, ny)
  3514. >>> xv, yv = np.meshgrid(x, y)
  3515. >>> xv
  3516. array([[0. , 0.5, 1. ],
  3517. [0. , 0.5, 1. ]])
  3518. >>> yv
  3519. array([[0., 0., 0.],
  3520. [1., 1., 1.]])
  3521. >>> xv, yv = np.meshgrid(x, y, sparse=True) # make sparse output arrays
  3522. >>> xv
  3523. array([[0. , 0.5, 1. ]])
  3524. >>> yv
  3525. array([[0.],
  3526. [1.]])
  3527. `meshgrid` is very useful to evaluate functions on a grid.
  3528. >>> import matplotlib.pyplot as plt
  3529. >>> x = np.arange(-5, 5, 0.1)
  3530. >>> y = np.arange(-5, 5, 0.1)
  3531. >>> xx, yy = np.meshgrid(x, y, sparse=True)
  3532. >>> z = np.sin(xx**2 + yy**2) / (xx**2 + yy**2)
  3533. >>> h = plt.contourf(x,y,z)
  3534. >>> plt.show()
  3535. """
  3536. ndim = len(xi)
  3537. if indexing not in ['xy', 'ij']:
  3538. raise ValueError(
  3539. "Valid values for `indexing` are 'xy' and 'ij'.")
  3540. s0 = (1,) * ndim
  3541. output = [np.asanyarray(x).reshape(s0[:i] + (-1,) + s0[i + 1:])
  3542. for i, x in enumerate(xi)]
  3543. if indexing == 'xy' and ndim > 1:
  3544. # switch first and second axis
  3545. output[0].shape = (1, -1) + s0[2:]
  3546. output[1].shape = (-1, 1) + s0[2:]
  3547. if not sparse:
  3548. # Return the full N-D matrix (not only the 1-D vector)
  3549. output = np.broadcast_arrays(*output, subok=True)
  3550. if copy:
  3551. output = [x.copy() for x in output]
  3552. return output
  3553. def _delete_dispatcher(arr, obj, axis=None):
  3554. return (arr, obj)
  3555. @array_function_dispatch(_delete_dispatcher)
  3556. def delete(arr, obj, axis=None):
  3557. """
  3558. Return a new array with sub-arrays along an axis deleted. For a one
  3559. dimensional array, this returns those entries not returned by
  3560. `arr[obj]`.
  3561. Parameters
  3562. ----------
  3563. arr : array_like
  3564. Input array.
  3565. obj : slice, int or array of ints
  3566. Indicate indices of sub-arrays to remove along the specified axis.
  3567. .. versionchanged:: 1.19.0
  3568. Boolean indices are now treated as a mask of elements to remove,
  3569. rather than being cast to the integers 0 and 1.
  3570. axis : int, optional
  3571. The axis along which to delete the subarray defined by `obj`.
  3572. If `axis` is None, `obj` is applied to the flattened array.
  3573. Returns
  3574. -------
  3575. out : ndarray
  3576. A copy of `arr` with the elements specified by `obj` removed. Note
  3577. that `delete` does not occur in-place. If `axis` is None, `out` is
  3578. a flattened array.
  3579. See Also
  3580. --------
  3581. insert : Insert elements into an array.
  3582. append : Append elements at the end of an array.
  3583. Notes
  3584. -----
  3585. Often it is preferable to use a boolean mask. For example:
  3586. >>> arr = np.arange(12) + 1
  3587. >>> mask = np.ones(len(arr), dtype=bool)
  3588. >>> mask[[0,2,4]] = False
  3589. >>> result = arr[mask,...]
  3590. Is equivalent to `np.delete(arr, [0,2,4], axis=0)`, but allows further
  3591. use of `mask`.
  3592. Examples
  3593. --------
  3594. >>> arr = np.array([[1,2,3,4], [5,6,7,8], [9,10,11,12]])
  3595. >>> arr
  3596. array([[ 1, 2, 3, 4],
  3597. [ 5, 6, 7, 8],
  3598. [ 9, 10, 11, 12]])
  3599. >>> np.delete(arr, 1, 0)
  3600. array([[ 1, 2, 3, 4],
  3601. [ 9, 10, 11, 12]])
  3602. >>> np.delete(arr, np.s_[::2], 1)
  3603. array([[ 2, 4],
  3604. [ 6, 8],
  3605. [10, 12]])
  3606. >>> np.delete(arr, [1,3,5], None)
  3607. array([ 1, 3, 5, 7, 8, 9, 10, 11, 12])
  3608. """
  3609. wrap = None
  3610. if type(arr) is not ndarray:
  3611. try:
  3612. wrap = arr.__array_wrap__
  3613. except AttributeError:
  3614. pass
  3615. arr = asarray(arr)
  3616. ndim = arr.ndim
  3617. arrorder = 'F' if arr.flags.fnc else 'C'
  3618. if axis is None:
  3619. if ndim != 1:
  3620. arr = arr.ravel()
  3621. # needed for np.matrix, which is still not 1d after being ravelled
  3622. ndim = arr.ndim
  3623. axis = ndim - 1
  3624. else:
  3625. axis = normalize_axis_index(axis, ndim)
  3626. slobj = [slice(None)]*ndim
  3627. N = arr.shape[axis]
  3628. newshape = list(arr.shape)
  3629. if isinstance(obj, slice):
  3630. start, stop, step = obj.indices(N)
  3631. xr = range(start, stop, step)
  3632. numtodel = len(xr)
  3633. if numtodel <= 0:
  3634. if wrap:
  3635. return wrap(arr.copy(order=arrorder))
  3636. else:
  3637. return arr.copy(order=arrorder)
  3638. # Invert if step is negative:
  3639. if step < 0:
  3640. step = -step
  3641. start = xr[-1]
  3642. stop = xr[0] + 1
  3643. newshape[axis] -= numtodel
  3644. new = empty(newshape, arr.dtype, arrorder)
  3645. # copy initial chunk
  3646. if start == 0:
  3647. pass
  3648. else:
  3649. slobj[axis] = slice(None, start)
  3650. new[tuple(slobj)] = arr[tuple(slobj)]
  3651. # copy end chunk
  3652. if stop == N:
  3653. pass
  3654. else:
  3655. slobj[axis] = slice(stop-numtodel, None)
  3656. slobj2 = [slice(None)]*ndim
  3657. slobj2[axis] = slice(stop, None)
  3658. new[tuple(slobj)] = arr[tuple(slobj2)]
  3659. # copy middle pieces
  3660. if step == 1:
  3661. pass
  3662. else: # use array indexing.
  3663. keep = ones(stop-start, dtype=bool)
  3664. keep[:stop-start:step] = False
  3665. slobj[axis] = slice(start, stop-numtodel)
  3666. slobj2 = [slice(None)]*ndim
  3667. slobj2[axis] = slice(start, stop)
  3668. arr = arr[tuple(slobj2)]
  3669. slobj2[axis] = keep
  3670. new[tuple(slobj)] = arr[tuple(slobj2)]
  3671. if wrap:
  3672. return wrap(new)
  3673. else:
  3674. return new
  3675. if isinstance(obj, (int, integer)) and not isinstance(obj, bool):
  3676. # optimization for a single value
  3677. if (obj < -N or obj >= N):
  3678. raise IndexError(
  3679. "index %i is out of bounds for axis %i with "
  3680. "size %i" % (obj, axis, N))
  3681. if (obj < 0):
  3682. obj += N
  3683. newshape[axis] -= 1
  3684. new = empty(newshape, arr.dtype, arrorder)
  3685. slobj[axis] = slice(None, obj)
  3686. new[tuple(slobj)] = arr[tuple(slobj)]
  3687. slobj[axis] = slice(obj, None)
  3688. slobj2 = [slice(None)]*ndim
  3689. slobj2[axis] = slice(obj+1, None)
  3690. new[tuple(slobj)] = arr[tuple(slobj2)]
  3691. else:
  3692. _obj = obj
  3693. obj = np.asarray(obj)
  3694. if obj.size == 0 and not isinstance(_obj, np.ndarray):
  3695. obj = obj.astype(intp)
  3696. if obj.dtype == bool:
  3697. if obj.shape != (N,):
  3698. raise ValueError('boolean array argument obj to delete '
  3699. 'must be one dimensional and match the axis '
  3700. 'length of {}'.format(N))
  3701. # optimization, the other branch is slower
  3702. keep = ~obj
  3703. else:
  3704. keep = ones(N, dtype=bool)
  3705. keep[obj,] = False
  3706. slobj[axis] = keep
  3707. new = arr[tuple(slobj)]
  3708. if wrap:
  3709. return wrap(new)
  3710. else:
  3711. return new
  3712. def _insert_dispatcher(arr, obj, values, axis=None):
  3713. return (arr, obj, values)
  3714. @array_function_dispatch(_insert_dispatcher)
  3715. def insert(arr, obj, values, axis=None):
  3716. """
  3717. Insert values along the given axis before the given indices.
  3718. Parameters
  3719. ----------
  3720. arr : array_like
  3721. Input array.
  3722. obj : int, slice or sequence of ints
  3723. Object that defines the index or indices before which `values` is
  3724. inserted.
  3725. .. versionadded:: 1.8.0
  3726. Support for multiple insertions when `obj` is a single scalar or a
  3727. sequence with one element (similar to calling insert multiple
  3728. times).
  3729. values : array_like
  3730. Values to insert into `arr`. If the type of `values` is different
  3731. from that of `arr`, `values` is converted to the type of `arr`.
  3732. `values` should be shaped so that ``arr[...,obj,...] = values``
  3733. is legal.
  3734. axis : int, optional
  3735. Axis along which to insert `values`. If `axis` is None then `arr`
  3736. is flattened first.
  3737. Returns
  3738. -------
  3739. out : ndarray
  3740. A copy of `arr` with `values` inserted. Note that `insert`
  3741. does not occur in-place: a new array is returned. If
  3742. `axis` is None, `out` is a flattened array.
  3743. See Also
  3744. --------
  3745. append : Append elements at the end of an array.
  3746. concatenate : Join a sequence of arrays along an existing axis.
  3747. delete : Delete elements from an array.
  3748. Notes
  3749. -----
  3750. Note that for higher dimensional inserts `obj=0` behaves very different
  3751. from `obj=[0]` just like `arr[:,0,:] = values` is different from
  3752. `arr[:,[0],:] = values`.
  3753. Examples
  3754. --------
  3755. >>> a = np.array([[1, 1], [2, 2], [3, 3]])
  3756. >>> a
  3757. array([[1, 1],
  3758. [2, 2],
  3759. [3, 3]])
  3760. >>> np.insert(a, 1, 5)
  3761. array([1, 5, 1, ..., 2, 3, 3])
  3762. >>> np.insert(a, 1, 5, axis=1)
  3763. array([[1, 5, 1],
  3764. [2, 5, 2],
  3765. [3, 5, 3]])
  3766. Difference between sequence and scalars:
  3767. >>> np.insert(a, [1], [[1],[2],[3]], axis=1)
  3768. array([[1, 1, 1],
  3769. [2, 2, 2],
  3770. [3, 3, 3]])
  3771. >>> np.array_equal(np.insert(a, 1, [1, 2, 3], axis=1),
  3772. ... np.insert(a, [1], [[1],[2],[3]], axis=1))
  3773. True
  3774. >>> b = a.flatten()
  3775. >>> b
  3776. array([1, 1, 2, 2, 3, 3])
  3777. >>> np.insert(b, [2, 2], [5, 6])
  3778. array([1, 1, 5, ..., 2, 3, 3])
  3779. >>> np.insert(b, slice(2, 4), [5, 6])
  3780. array([1, 1, 5, ..., 2, 3, 3])
  3781. >>> np.insert(b, [2, 2], [7.13, False]) # type casting
  3782. array([1, 1, 7, ..., 2, 3, 3])
  3783. >>> x = np.arange(8).reshape(2, 4)
  3784. >>> idx = (1, 3)
  3785. >>> np.insert(x, idx, 999, axis=1)
  3786. array([[ 0, 999, 1, 2, 999, 3],
  3787. [ 4, 999, 5, 6, 999, 7]])
  3788. """
  3789. wrap = None
  3790. if type(arr) is not ndarray:
  3791. try:
  3792. wrap = arr.__array_wrap__
  3793. except AttributeError:
  3794. pass
  3795. arr = asarray(arr)
  3796. ndim = arr.ndim
  3797. arrorder = 'F' if arr.flags.fnc else 'C'
  3798. if axis is None:
  3799. if ndim != 1:
  3800. arr = arr.ravel()
  3801. # needed for np.matrix, which is still not 1d after being ravelled
  3802. ndim = arr.ndim
  3803. axis = ndim - 1
  3804. else:
  3805. axis = normalize_axis_index(axis, ndim)
  3806. slobj = [slice(None)]*ndim
  3807. N = arr.shape[axis]
  3808. newshape = list(arr.shape)
  3809. if isinstance(obj, slice):
  3810. # turn it into a range object
  3811. indices = arange(*obj.indices(N), dtype=intp)
  3812. else:
  3813. # need to copy obj, because indices will be changed in-place
  3814. indices = np.array(obj)
  3815. if indices.dtype == bool:
  3816. # See also delete
  3817. # 2012-10-11, NumPy 1.8
  3818. warnings.warn(
  3819. "in the future insert will treat boolean arrays and "
  3820. "array-likes as a boolean index instead of casting it to "
  3821. "integer", FutureWarning, stacklevel=3)
  3822. indices = indices.astype(intp)
  3823. # Code after warning period:
  3824. #if obj.ndim != 1:
  3825. # raise ValueError('boolean array argument obj to insert '
  3826. # 'must be one dimensional')
  3827. #indices = np.flatnonzero(obj)
  3828. elif indices.ndim > 1:
  3829. raise ValueError(
  3830. "index array argument obj to insert must be one dimensional "
  3831. "or scalar")
  3832. if indices.size == 1:
  3833. index = indices.item()
  3834. if index < -N or index > N:
  3835. raise IndexError(
  3836. "index %i is out of bounds for axis %i with "
  3837. "size %i" % (obj, axis, N))
  3838. if (index < 0):
  3839. index += N
  3840. # There are some object array corner cases here, but we cannot avoid
  3841. # that:
  3842. values = array(values, copy=False, ndmin=arr.ndim, dtype=arr.dtype)
  3843. if indices.ndim == 0:
  3844. # broadcasting is very different here, since a[:,0,:] = ... behaves
  3845. # very different from a[:,[0],:] = ...! This changes values so that
  3846. # it works likes the second case. (here a[:,0:1,:])
  3847. values = np.moveaxis(values, 0, axis)
  3848. numnew = values.shape[axis]
  3849. newshape[axis] += numnew
  3850. new = empty(newshape, arr.dtype, arrorder)
  3851. slobj[axis] = slice(None, index)
  3852. new[tuple(slobj)] = arr[tuple(slobj)]
  3853. slobj[axis] = slice(index, index+numnew)
  3854. new[tuple(slobj)] = values
  3855. slobj[axis] = slice(index+numnew, None)
  3856. slobj2 = [slice(None)] * ndim
  3857. slobj2[axis] = slice(index, None)
  3858. new[tuple(slobj)] = arr[tuple(slobj2)]
  3859. if wrap:
  3860. return wrap(new)
  3861. return new
  3862. elif indices.size == 0 and not isinstance(obj, np.ndarray):
  3863. # Can safely cast the empty list to intp
  3864. indices = indices.astype(intp)
  3865. indices[indices < 0] += N
  3866. numnew = len(indices)
  3867. order = indices.argsort(kind='mergesort') # stable sort
  3868. indices[order] += np.arange(numnew)
  3869. newshape[axis] += numnew
  3870. old_mask = ones(newshape[axis], dtype=bool)
  3871. old_mask[indices] = False
  3872. new = empty(newshape, arr.dtype, arrorder)
  3873. slobj2 = [slice(None)]*ndim
  3874. slobj[axis] = indices
  3875. slobj2[axis] = old_mask
  3876. new[tuple(slobj)] = values
  3877. new[tuple(slobj2)] = arr
  3878. if wrap:
  3879. return wrap(new)
  3880. return new
  3881. def _append_dispatcher(arr, values, axis=None):
  3882. return (arr, values)
  3883. @array_function_dispatch(_append_dispatcher)
  3884. def append(arr, values, axis=None):
  3885. """
  3886. Append values to the end of an array.
  3887. Parameters
  3888. ----------
  3889. arr : array_like
  3890. Values are appended to a copy of this array.
  3891. values : array_like
  3892. These values are appended to a copy of `arr`. It must be of the
  3893. correct shape (the same shape as `arr`, excluding `axis`). If
  3894. `axis` is not specified, `values` can be any shape and will be
  3895. flattened before use.
  3896. axis : int, optional
  3897. The axis along which `values` are appended. If `axis` is not
  3898. given, both `arr` and `values` are flattened before use.
  3899. Returns
  3900. -------
  3901. append : ndarray
  3902. A copy of `arr` with `values` appended to `axis`. Note that
  3903. `append` does not occur in-place: a new array is allocated and
  3904. filled. If `axis` is None, `out` is a flattened array.
  3905. See Also
  3906. --------
  3907. insert : Insert elements into an array.
  3908. delete : Delete elements from an array.
  3909. Examples
  3910. --------
  3911. >>> np.append([1, 2, 3], [[4, 5, 6], [7, 8, 9]])
  3912. array([1, 2, 3, ..., 7, 8, 9])
  3913. When `axis` is specified, `values` must have the correct shape.
  3914. >>> np.append([[1, 2, 3], [4, 5, 6]], [[7, 8, 9]], axis=0)
  3915. array([[1, 2, 3],
  3916. [4, 5, 6],
  3917. [7, 8, 9]])
  3918. >>> np.append([[1, 2, 3], [4, 5, 6]], [7, 8, 9], axis=0)
  3919. Traceback (most recent call last):
  3920. ...
  3921. ValueError: all the input arrays must have same number of dimensions, but
  3922. the array at index 0 has 2 dimension(s) and the array at index 1 has 1
  3923. dimension(s)
  3924. """
  3925. arr = asanyarray(arr)
  3926. if axis is None:
  3927. if arr.ndim != 1:
  3928. arr = arr.ravel()
  3929. values = ravel(values)
  3930. axis = arr.ndim-1
  3931. return concatenate((arr, values), axis=axis)
  3932. def _digitize_dispatcher(x, bins, right=None):
  3933. return (x, bins)
  3934. @array_function_dispatch(_digitize_dispatcher)
  3935. def digitize(x, bins, right=False):
  3936. """
  3937. Return the indices of the bins to which each value in input array belongs.
  3938. ========= ============= ============================
  3939. `right` order of bins returned index `i` satisfies
  3940. ========= ============= ============================
  3941. ``False`` increasing ``bins[i-1] <= x < bins[i]``
  3942. ``True`` increasing ``bins[i-1] < x <= bins[i]``
  3943. ``False`` decreasing ``bins[i-1] > x >= bins[i]``
  3944. ``True`` decreasing ``bins[i-1] >= x > bins[i]``
  3945. ========= ============= ============================
  3946. If values in `x` are beyond the bounds of `bins`, 0 or ``len(bins)`` is
  3947. returned as appropriate.
  3948. Parameters
  3949. ----------
  3950. x : array_like
  3951. Input array to be binned. Prior to NumPy 1.10.0, this array had to
  3952. be 1-dimensional, but can now have any shape.
  3953. bins : array_like
  3954. Array of bins. It has to be 1-dimensional and monotonic.
  3955. right : bool, optional
  3956. Indicating whether the intervals include the right or the left bin
  3957. edge. Default behavior is (right==False) indicating that the interval
  3958. does not include the right edge. The left bin end is open in this
  3959. case, i.e., bins[i-1] <= x < bins[i] is the default behavior for
  3960. monotonically increasing bins.
  3961. Returns
  3962. -------
  3963. indices : ndarray of ints
  3964. Output array of indices, of same shape as `x`.
  3965. Raises
  3966. ------
  3967. ValueError
  3968. If `bins` is not monotonic.
  3969. TypeError
  3970. If the type of the input is complex.
  3971. See Also
  3972. --------
  3973. bincount, histogram, unique, searchsorted
  3974. Notes
  3975. -----
  3976. If values in `x` are such that they fall outside the bin range,
  3977. attempting to index `bins` with the indices that `digitize` returns
  3978. will result in an IndexError.
  3979. .. versionadded:: 1.10.0
  3980. `np.digitize` is implemented in terms of `np.searchsorted`. This means
  3981. that a binary search is used to bin the values, which scales much better
  3982. for larger number of bins than the previous linear search. It also removes
  3983. the requirement for the input array to be 1-dimensional.
  3984. For monotonically _increasing_ `bins`, the following are equivalent::
  3985. np.digitize(x, bins, right=True)
  3986. np.searchsorted(bins, x, side='left')
  3987. Note that as the order of the arguments are reversed, the side must be too.
  3988. The `searchsorted` call is marginally faster, as it does not do any
  3989. monotonicity checks. Perhaps more importantly, it supports all dtypes.
  3990. Examples
  3991. --------
  3992. >>> x = np.array([0.2, 6.4, 3.0, 1.6])
  3993. >>> bins = np.array([0.0, 1.0, 2.5, 4.0, 10.0])
  3994. >>> inds = np.digitize(x, bins)
  3995. >>> inds
  3996. array([1, 4, 3, 2])
  3997. >>> for n in range(x.size):
  3998. ... print(bins[inds[n]-1], "<=", x[n], "<", bins[inds[n]])
  3999. ...
  4000. 0.0 <= 0.2 < 1.0
  4001. 4.0 <= 6.4 < 10.0
  4002. 2.5 <= 3.0 < 4.0
  4003. 1.0 <= 1.6 < 2.5
  4004. >>> x = np.array([1.2, 10.0, 12.4, 15.5, 20.])
  4005. >>> bins = np.array([0, 5, 10, 15, 20])
  4006. >>> np.digitize(x,bins,right=True)
  4007. array([1, 2, 3, 4, 4])
  4008. >>> np.digitize(x,bins,right=False)
  4009. array([1, 3, 3, 4, 5])
  4010. """
  4011. x = _nx.asarray(x)
  4012. bins = _nx.asarray(bins)
  4013. # here for compatibility, searchsorted below is happy to take this
  4014. if np.issubdtype(x.dtype, _nx.complexfloating):
  4015. raise TypeError("x may not be complex")
  4016. mono = _monotonicity(bins)
  4017. if mono == 0:
  4018. raise ValueError("bins must be monotonically increasing or decreasing")
  4019. # this is backwards because the arguments below are swapped
  4020. side = 'left' if right else 'right'
  4021. if mono == -1:
  4022. # reverse the bins, and invert the results
  4023. return len(bins) - _nx.searchsorted(bins[::-1], x, side=side)
  4024. else:
  4025. return _nx.searchsorted(bins, x, side=side)