123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735736737738739740741742743744745746747748749750751752753754755756757758759760761762763764765766767768769770771772773774775776777778779780781782783784785786787788789790791792793794795796797798799800801802803804805806807808809810811812813814815816817818819820821822823824825826827828829830831832833834835836837838839840841842843844845846847848849850851852853854855856857858859860861862863864865866867868869870871872873874875876877878879880881882883884885886887888889890891892893894895896897898899900901902903904905906907908909910911912913914915916917918919920921922923924925926927928929930931932933934935936937938939940941942943944945946947948949950951952953954955956957958959960961962963964965966967968969970971972973974975976977978979980981982983984985986987988989990991992993994995996997998999100010011002100310041005100610071008100910101011101210131014101510161017101810191020102110221023102410251026102710281029103010311032103310341035103610371038103910401041104210431044104510461047104810491050105110521053105410551056105710581059106010611062106310641065106610671068106910701071107210731074107510761077107810791080108110821083108410851086108710881089109010911092109310941095109610971098109911001101110211031104110511061107110811091110111111121113111411151116111711181119112011211122112311241125112611271128112911301131113211331134113511361137113811391140114111421143114411451146114711481149115011511152115311541155115611571158115911601161116211631164116511661167116811691170117111721173117411751176117711781179118011811182118311841185118611871188118911901191119211931194119511961197119811991200120112021203120412051206120712081209121012111212121312141215121612171218121912201221122212231224122512261227122812291230123112321233123412351236123712381239124012411242124312441245124612471248124912501251125212531254125512561257125812591260126112621263126412651266126712681269127012711272127312741275127612771278127912801281128212831284128512861287128812891290129112921293129412951296129712981299130013011302130313041305130613071308130913101311131213131314131513161317131813191320132113221323132413251326132713281329133013311332133313341335133613371338133913401341134213431344134513461347134813491350135113521353135413551356135713581359136013611362136313641365136613671368136913701371137213731374137513761377137813791380138113821383138413851386138713881389139013911392139313941395139613971398139914001401140214031404140514061407140814091410141114121413141414151416141714181419142014211422142314241425142614271428142914301431143214331434143514361437143814391440144114421443144414451446144714481449145014511452145314541455145614571458145914601461146214631464146514661467146814691470147114721473147414751476147714781479148014811482148314841485148614871488148914901491149214931494149514961497149814991500150115021503150415051506150715081509151015111512151315141515151615171518151915201521152215231524152515261527152815291530153115321533153415351536153715381539154015411542154315441545154615471548154915501551155215531554155515561557155815591560156115621563156415651566156715681569157015711572157315741575157615771578157915801581158215831584158515861587158815891590159115921593159415951596159715981599160016011602160316041605160616071608160916101611161216131614161516161617161816191620162116221623162416251626162716281629163016311632163316341635163616371638163916401641164216431644164516461647164816491650165116521653165416551656165716581659166016611662166316641665166616671668166916701671167216731674167516761677167816791680168116821683168416851686168716881689169016911692169316941695169616971698169917001701170217031704170517061707170817091710171117121713171417151716171717181719172017211722172317241725172617271728172917301731173217331734173517361737173817391740174117421743174417451746174717481749175017511752175317541755175617571758175917601761176217631764176517661767176817691770177117721773177417751776177717781779178017811782178317841785178617871788178917901791179217931794179517961797179817991800180118021803180418051806180718081809181018111812181318141815181618171818181918201821182218231824182518261827182818291830183118321833183418351836183718381839184018411842184318441845184618471848184918501851185218531854185518561857185818591860186118621863186418651866186718681869187018711872187318741875187618771878187918801881188218831884188518861887188818891890189118921893189418951896189718981899190019011902190319041905190619071908190919101911191219131914191519161917191819191920192119221923192419251926192719281929193019311932193319341935193619371938193919401941194219431944194519461947194819491950195119521953195419551956195719581959196019611962196319641965196619671968196919701971197219731974197519761977197819791980198119821983198419851986198719881989199019911992199319941995199619971998199920002001200220032004200520062007200820092010201120122013201420152016201720182019202020212022202320242025202620272028202920302031203220332034203520362037203820392040204120422043204420452046204720482049205020512052205320542055205620572058205920602061206220632064206520662067206820692070207120722073207420752076207720782079208020812082208320842085208620872088208920902091209220932094209520962097209820992100210121022103210421052106210721082109211021112112211321142115211621172118211921202121212221232124212521262127212821292130213121322133213421352136213721382139214021412142214321442145214621472148214921502151215221532154215521562157215821592160216121622163216421652166216721682169217021712172217321742175217621772178217921802181218221832184218521862187218821892190219121922193219421952196219721982199220022012202220322042205220622072208220922102211221222132214221522162217221822192220222122222223222422252226222722282229223022312232223322342235223622372238223922402241224222432244224522462247224822492250225122522253225422552256225722582259226022612262226322642265226622672268226922702271227222732274227522762277227822792280228122822283228422852286228722882289229022912292229322942295229622972298229923002301230223032304230523062307230823092310231123122313231423152316231723182319232023212322232323242325232623272328232923302331233223332334233523362337233823392340234123422343234423452346234723482349235023512352235323542355235623572358235923602361236223632364236523662367236823692370237123722373237423752376237723782379238023812382238323842385238623872388238923902391239223932394239523962397239823992400240124022403240424052406240724082409241024112412241324142415241624172418241924202421242224232424242524262427242824292430243124322433243424352436243724382439244024412442244324442445244624472448244924502451245224532454245524562457245824592460246124622463246424652466246724682469247024712472247324742475247624772478247924802481248224832484248524862487248824892490249124922493249424952496249724982499250025012502250325042505250625072508250925102511251225132514251525162517251825192520252125222523252425252526252725282529253025312532253325342535253625372538253925402541254225432544254525462547254825492550255125522553255425552556255725582559256025612562256325642565256625672568256925702571257225732574257525762577257825792580258125822583258425852586258725882589259025912592259325942595259625972598259926002601260226032604260526062607260826092610261126122613261426152616261726182619262026212622262326242625262626272628262926302631263226332634263526362637263826392640264126422643264426452646264726482649265026512652265326542655265626572658265926602661266226632664266526662667266826692670267126722673267426752676267726782679268026812682268326842685268626872688268926902691269226932694269526962697269826992700270127022703270427052706270727082709271027112712271327142715271627172718271927202721272227232724272527262727272827292730273127322733273427352736273727382739274027412742274327442745274627472748274927502751275227532754275527562757275827592760276127622763276427652766276727682769277027712772277327742775277627772778277927802781278227832784278527862787278827892790279127922793279427952796279727982799280028012802280328042805280628072808280928102811281228132814281528162817281828192820282128222823282428252826282728282829283028312832283328342835283628372838283928402841284228432844284528462847284828492850285128522853285428552856285728582859286028612862286328642865286628672868286928702871287228732874287528762877287828792880288128822883288428852886288728882889289028912892289328942895289628972898289929002901290229032904290529062907290829092910291129122913291429152916291729182919292029212922292329242925292629272928292929302931293229332934293529362937293829392940294129422943294429452946294729482949295029512952295329542955295629572958295929602961296229632964296529662967296829692970297129722973297429752976297729782979298029812982298329842985298629872988298929902991299229932994299529962997299829993000300130023003300430053006300730083009301030113012301330143015301630173018301930203021302230233024302530263027302830293030303130323033303430353036303730383039304030413042304330443045304630473048304930503051305230533054305530563057305830593060306130623063306430653066306730683069307030713072307330743075307630773078307930803081308230833084308530863087308830893090309130923093309430953096309730983099310031013102310331043105310631073108310931103111311231133114311531163117311831193120312131223123312431253126312731283129313031313132313331343135313631373138313931403141314231433144314531463147314831493150315131523153315431553156315731583159316031613162316331643165316631673168316931703171317231733174317531763177317831793180318131823183318431853186318731883189319031913192319331943195319631973198319932003201320232033204320532063207320832093210321132123213321432153216321732183219322032213222322332243225322632273228322932303231323232333234323532363237323832393240324132423243324432453246324732483249325032513252325332543255325632573258325932603261326232633264326532663267326832693270327132723273327432753276327732783279328032813282328332843285328632873288328932903291329232933294329532963297329832993300330133023303330433053306330733083309331033113312331333143315331633173318331933203321332233233324332533263327332833293330333133323333333433353336333733383339334033413342334333443345334633473348334933503351335233533354335533563357335833593360336133623363336433653366336733683369337033713372337333743375337633773378337933803381338233833384338533863387338833893390339133923393339433953396339733983399340034013402340334043405340634073408340934103411341234133414341534163417341834193420342134223423342434253426342734283429343034313432343334343435343634373438343934403441344234433444344534463447344834493450345134523453345434553456345734583459346034613462346334643465346634673468346934703471347234733474347534763477347834793480348134823483348434853486348734883489349034913492349334943495349634973498349935003501350235033504350535063507350835093510351135123513351435153516351735183519352035213522352335243525352635273528352935303531353235333534353535363537353835393540354135423543354435453546354735483549355035513552355335543555355635573558355935603561356235633564356535663567356835693570357135723573357435753576357735783579358035813582358335843585358635873588358935903591359235933594359535963597359835993600360136023603360436053606360736083609361036113612361336143615361636173618361936203621362236233624362536263627362836293630363136323633363436353636363736383639364036413642364336443645364636473648364936503651365236533654365536563657365836593660366136623663366436653666366736683669367036713672367336743675367636773678367936803681368236833684368536863687368836893690369136923693369436953696369736983699370037013702370337043705370637073708370937103711371237133714371537163717371837193720372137223723372437253726372737283729373037313732373337343735373637373738373937403741374237433744374537463747374837493750375137523753375437553756375737583759376037613762376337643765376637673768376937703771377237733774377537763777377837793780378137823783378437853786378737883789379037913792379337943795379637973798379938003801380238033804380538063807380838093810381138123813381438153816381738183819382038213822382338243825382638273828382938303831383238333834383538363837383838393840384138423843384438453846384738483849385038513852385338543855385638573858385938603861386238633864386538663867386838693870387138723873387438753876387738783879388038813882388338843885388638873888388938903891389238933894389538963897389838993900390139023903390439053906390739083909391039113912391339143915391639173918391939203921392239233924392539263927392839293930393139323933393439353936393739383939394039413942394339443945394639473948394939503951395239533954395539563957395839593960396139623963396439653966396739683969397039713972397339743975397639773978397939803981398239833984398539863987398839893990399139923993399439953996399739983999400040014002400340044005400640074008400940104011401240134014401540164017401840194020402140224023402440254026402740284029403040314032403340344035403640374038403940404041404240434044404540464047404840494050405140524053405440554056405740584059406040614062406340644065406640674068406940704071407240734074407540764077407840794080408140824083408440854086408740884089409040914092409340944095409640974098409941004101410241034104410541064107410841094110411141124113411441154116411741184119412041214122412341244125412641274128412941304131413241334134413541364137413841394140414141424143414441454146414741484149415041514152415341544155415641574158415941604161416241634164416541664167416841694170417141724173417441754176417741784179418041814182418341844185418641874188418941904191419241934194419541964197419841994200420142024203420442054206420742084209421042114212421342144215421642174218421942204221422242234224422542264227422842294230423142324233423442354236423742384239424042414242424342444245424642474248424942504251425242534254425542564257425842594260426142624263426442654266426742684269427042714272427342744275427642774278427942804281428242834284428542864287428842894290429142924293429442954296429742984299430043014302430343044305430643074308430943104311431243134314431543164317431843194320432143224323432443254326432743284329433043314332433343344335433643374338433943404341434243434344434543464347434843494350435143524353435443554356435743584359436043614362436343644365436643674368436943704371437243734374437543764377437843794380438143824383438443854386438743884389439043914392439343944395439643974398439944004401440244034404440544064407440844094410441144124413441444154416441744184419442044214422442344244425442644274428442944304431443244334434443544364437443844394440444144424443444444454446444744484449445044514452445344544455445644574458445944604461446244634464446544664467446844694470447144724473447444754476447744784479448044814482448344844485448644874488448944904491449244934494449544964497449844994500450145024503450445054506450745084509451045114512451345144515451645174518451945204521452245234524452545264527452845294530453145324533453445354536453745384539454045414542454345444545454645474548454945504551455245534554455545564557455845594560456145624563456445654566456745684569457045714572457345744575457645774578457945804581458245834584458545864587458845894590459145924593459445954596459745984599460046014602460346044605460646074608460946104611461246134614461546164617461846194620462146224623462446254626462746284629463046314632463346344635463646374638463946404641464246434644464546464647464846494650465146524653465446554656465746584659466046614662466346644665466646674668466946704671467246734674467546764677467846794680468146824683468446854686468746884689469046914692469346944695469646974698469947004701470247034704470547064707470847094710471147124713471447154716471747184719472047214722472347244725472647274728472947304731473247334734473547364737473847394740474147424743474447454746474747484749475047514752475347544755475647574758475947604761476247634764476547664767476847694770477147724773477447754776477747784779478047814782478347844785478647874788478947904791479247934794479547964797479847994800480148024803480448054806480748084809481048114812481348144815481648174818481948204821482248234824482548264827482848294830483148324833483448354836483748384839484048414842484348444845484648474848484948504851485248534854485548564857485848594860 |
- import collections.abc
- import functools
- import re
- import sys
- import warnings
- import numpy as np
- import numpy.core.numeric as _nx
- from numpy.core import transpose
- from numpy.core.numeric import (
- ones, zeros, arange, concatenate, array, asarray, asanyarray, empty,
- ndarray, around, floor, ceil, take, dot, where, intp,
- integer, isscalar, absolute
- )
- from numpy.core.umath import (
- pi, add, arctan2, frompyfunc, cos, less_equal, sqrt, sin,
- mod, exp, not_equal, subtract
- )
- from numpy.core.fromnumeric import (
- ravel, nonzero, partition, mean, any, sum
- )
- from numpy.core.numerictypes import typecodes
- from numpy.core.overrides import set_module
- from numpy.core import overrides
- from numpy.core.function_base import add_newdoc
- from numpy.lib.twodim_base import diag
- from numpy.core.multiarray import (
- _insert, add_docstring, bincount, normalize_axis_index, _monotonicity,
- interp as compiled_interp, interp_complex as compiled_interp_complex
- )
- from numpy.core.umath import _add_newdoc_ufunc as add_newdoc_ufunc
- import builtins
- # needed in this module for compatibility
- from numpy.lib.histograms import histogram, histogramdd
- array_function_dispatch = functools.partial(
- overrides.array_function_dispatch, module='numpy')
- __all__ = [
- 'select', 'piecewise', 'trim_zeros', 'copy', 'iterable', 'percentile',
- 'diff', 'gradient', 'angle', 'unwrap', 'sort_complex', 'disp', 'flip',
- 'rot90', 'extract', 'place', 'vectorize', 'asarray_chkfinite', 'average',
- 'bincount', 'digitize', 'cov', 'corrcoef',
- 'msort', 'median', 'sinc', 'hamming', 'hanning', 'bartlett',
- 'blackman', 'kaiser', 'trapz', 'i0', 'add_newdoc', 'add_docstring',
- 'meshgrid', 'delete', 'insert', 'append', 'interp', 'add_newdoc_ufunc',
- 'quantile'
- ]
- def _rot90_dispatcher(m, k=None, axes=None):
- return (m,)
- @array_function_dispatch(_rot90_dispatcher)
- def rot90(m, k=1, axes=(0, 1)):
- """
- Rotate an array by 90 degrees in the plane specified by axes.
- Rotation direction is from the first towards the second axis.
- Parameters
- ----------
- m : array_like
- Array of two or more dimensions.
- k : integer
- Number of times the array is rotated by 90 degrees.
- axes: (2,) array_like
- The array is rotated in the plane defined by the axes.
- Axes must be different.
- .. versionadded:: 1.12.0
- Returns
- -------
- y : ndarray
- A rotated view of `m`.
- See Also
- --------
- flip : Reverse the order of elements in an array along the given axis.
- fliplr : Flip an array horizontally.
- flipud : Flip an array vertically.
- Notes
- -----
- rot90(m, k=1, axes=(1,0)) is the reverse of rot90(m, k=1, axes=(0,1))
- rot90(m, k=1, axes=(1,0)) is equivalent to rot90(m, k=-1, axes=(0,1))
- Examples
- --------
- >>> m = np.array([[1,2],[3,4]], int)
- >>> m
- array([[1, 2],
- [3, 4]])
- >>> np.rot90(m)
- array([[2, 4],
- [1, 3]])
- >>> np.rot90(m, 2)
- array([[4, 3],
- [2, 1]])
- >>> m = np.arange(8).reshape((2,2,2))
- >>> np.rot90(m, 1, (1,2))
- array([[[1, 3],
- [0, 2]],
- [[5, 7],
- [4, 6]]])
- """
- axes = tuple(axes)
- if len(axes) != 2:
- raise ValueError("len(axes) must be 2.")
- m = asanyarray(m)
- if axes[0] == axes[1] or absolute(axes[0] - axes[1]) == m.ndim:
- raise ValueError("Axes must be different.")
- if (axes[0] >= m.ndim or axes[0] < -m.ndim
- or axes[1] >= m.ndim or axes[1] < -m.ndim):
- raise ValueError("Axes={} out of range for array of ndim={}."
- .format(axes, m.ndim))
- k %= 4
- if k == 0:
- return m[:]
- if k == 2:
- return flip(flip(m, axes[0]), axes[1])
- axes_list = arange(0, m.ndim)
- (axes_list[axes[0]], axes_list[axes[1]]) = (axes_list[axes[1]],
- axes_list[axes[0]])
- if k == 1:
- return transpose(flip(m, axes[1]), axes_list)
- else:
- # k == 3
- return flip(transpose(m, axes_list), axes[1])
- def _flip_dispatcher(m, axis=None):
- return (m,)
- @array_function_dispatch(_flip_dispatcher)
- def flip(m, axis=None):
- """
- Reverse the order of elements in an array along the given axis.
- The shape of the array is preserved, but the elements are reordered.
- .. versionadded:: 1.12.0
- Parameters
- ----------
- m : array_like
- Input array.
- axis : None or int or tuple of ints, optional
- Axis or axes along which to flip over. The default,
- axis=None, will flip over all of the axes of the input array.
- If axis is negative it counts from the last to the first axis.
- If axis is a tuple of ints, flipping is performed on all of the axes
- specified in the tuple.
- .. versionchanged:: 1.15.0
- None and tuples of axes are supported
- Returns
- -------
- out : array_like
- A view of `m` with the entries of axis reversed. Since a view is
- returned, this operation is done in constant time.
- See Also
- --------
- flipud : Flip an array vertically (axis=0).
- fliplr : Flip an array horizontally (axis=1).
- Notes
- -----
- flip(m, 0) is equivalent to flipud(m).
- flip(m, 1) is equivalent to fliplr(m).
- flip(m, n) corresponds to ``m[...,::-1,...]`` with ``::-1`` at position n.
- flip(m) corresponds to ``m[::-1,::-1,...,::-1]`` with ``::-1`` at all
- positions.
- flip(m, (0, 1)) corresponds to ``m[::-1,::-1,...]`` with ``::-1`` at
- position 0 and position 1.
- Examples
- --------
- >>> A = np.arange(8).reshape((2,2,2))
- >>> A
- array([[[0, 1],
- [2, 3]],
- [[4, 5],
- [6, 7]]])
- >>> np.flip(A, 0)
- array([[[4, 5],
- [6, 7]],
- [[0, 1],
- [2, 3]]])
- >>> np.flip(A, 1)
- array([[[2, 3],
- [0, 1]],
- [[6, 7],
- [4, 5]]])
- >>> np.flip(A)
- array([[[7, 6],
- [5, 4]],
- [[3, 2],
- [1, 0]]])
- >>> np.flip(A, (0, 2))
- array([[[5, 4],
- [7, 6]],
- [[1, 0],
- [3, 2]]])
- >>> A = np.random.randn(3,4,5)
- >>> np.all(np.flip(A,2) == A[:,:,::-1,...])
- True
- """
- if not hasattr(m, 'ndim'):
- m = asarray(m)
- if axis is None:
- indexer = (np.s_[::-1],) * m.ndim
- else:
- axis = _nx.normalize_axis_tuple(axis, m.ndim)
- indexer = [np.s_[:]] * m.ndim
- for ax in axis:
- indexer[ax] = np.s_[::-1]
- indexer = tuple(indexer)
- return m[indexer]
- @set_module('numpy')
- def iterable(y):
- """
- Check whether or not an object can be iterated over.
- Parameters
- ----------
- y : object
- Input object.
- Returns
- -------
- b : bool
- Return ``True`` if the object has an iterator method or is a
- sequence and ``False`` otherwise.
- Examples
- --------
- >>> np.iterable([1, 2, 3])
- True
- >>> np.iterable(2)
- False
- """
- try:
- iter(y)
- except TypeError:
- return False
- return True
- def _average_dispatcher(a, axis=None, weights=None, returned=None):
- return (a, weights)
- @array_function_dispatch(_average_dispatcher)
- def average(a, axis=None, weights=None, returned=False):
- """
- Compute the weighted average along the specified axis.
- Parameters
- ----------
- a : array_like
- Array containing data to be averaged. If `a` is not an array, a
- conversion is attempted.
- axis : None or int or tuple of ints, optional
- Axis or axes along which to average `a`. The default,
- axis=None, will average over all of the elements of the input array.
- If axis is negative it counts from the last to the first axis.
- .. versionadded:: 1.7.0
- If axis is a tuple of ints, averaging is performed on all of the axes
- specified in the tuple instead of a single axis or all the axes as
- before.
- weights : array_like, optional
- An array of weights associated with the values in `a`. Each value in
- `a` contributes to the average according to its associated weight.
- The weights array can either be 1-D (in which case its length must be
- the size of `a` along the given axis) or of the same shape as `a`.
- If `weights=None`, then all data in `a` are assumed to have a
- weight equal to one. The 1-D calculation is::
- avg = sum(a * weights) / sum(weights)
- The only constraint on `weights` is that `sum(weights)` must not be 0.
- returned : bool, optional
- Default is `False`. If `True`, the tuple (`average`, `sum_of_weights`)
- is returned, otherwise only the average is returned.
- If `weights=None`, `sum_of_weights` is equivalent to the number of
- elements over which the average is taken.
- Returns
- -------
- retval, [sum_of_weights] : array_type or double
- Return the average along the specified axis. When `returned` is `True`,
- return a tuple with the average as the first element and the sum
- of the weights as the second element. `sum_of_weights` is of the
- same type as `retval`. The result dtype follows a genereal pattern.
- If `weights` is None, the result dtype will be that of `a` , or ``float64``
- if `a` is integral. Otherwise, if `weights` is not None and `a` is non-
- integral, the result type will be the type of lowest precision capable of
- representing values of both `a` and `weights`. If `a` happens to be
- integral, the previous rules still applies but the result dtype will
- at least be ``float64``.
- Raises
- ------
- ZeroDivisionError
- When all weights along axis are zero. See `numpy.ma.average` for a
- version robust to this type of error.
- TypeError
- When the length of 1D `weights` is not the same as the shape of `a`
- along axis.
- See Also
- --------
- mean
- ma.average : average for masked arrays -- useful if your data contains
- "missing" values
- numpy.result_type : Returns the type that results from applying the
- numpy type promotion rules to the arguments.
- Examples
- --------
- >>> data = np.arange(1, 5)
- >>> data
- array([1, 2, 3, 4])
- >>> np.average(data)
- 2.5
- >>> np.average(np.arange(1, 11), weights=np.arange(10, 0, -1))
- 4.0
- >>> data = np.arange(6).reshape((3,2))
- >>> data
- array([[0, 1],
- [2, 3],
- [4, 5]])
- >>> np.average(data, axis=1, weights=[1./4, 3./4])
- array([0.75, 2.75, 4.75])
- >>> np.average(data, weights=[1./4, 3./4])
- Traceback (most recent call last):
- ...
- TypeError: Axis must be specified when shapes of a and weights differ.
- >>> a = np.ones(5, dtype=np.float128)
- >>> w = np.ones(5, dtype=np.complex64)
- >>> avg = np.average(a, weights=w)
- >>> print(avg.dtype)
- complex256
- """
- a = np.asanyarray(a)
- if weights is None:
- avg = a.mean(axis)
- scl = avg.dtype.type(a.size/avg.size)
- else:
- wgt = np.asanyarray(weights)
- if issubclass(a.dtype.type, (np.integer, np.bool_)):
- result_dtype = np.result_type(a.dtype, wgt.dtype, 'f8')
- else:
- result_dtype = np.result_type(a.dtype, wgt.dtype)
- # Sanity checks
- if a.shape != wgt.shape:
- if axis is None:
- raise TypeError(
- "Axis must be specified when shapes of a and weights "
- "differ.")
- if wgt.ndim != 1:
- raise TypeError(
- "1D weights expected when shapes of a and weights differ.")
- if wgt.shape[0] != a.shape[axis]:
- raise ValueError(
- "Length of weights not compatible with specified axis.")
- # setup wgt to broadcast along axis
- wgt = np.broadcast_to(wgt, (a.ndim-1)*(1,) + wgt.shape)
- wgt = wgt.swapaxes(-1, axis)
- scl = wgt.sum(axis=axis, dtype=result_dtype)
- if np.any(scl == 0.0):
- raise ZeroDivisionError(
- "Weights sum to zero, can't be normalized")
- avg = np.multiply(a, wgt, dtype=result_dtype).sum(axis)/scl
- if returned:
- if scl.shape != avg.shape:
- scl = np.broadcast_to(scl, avg.shape).copy()
- return avg, scl
- else:
- return avg
- @set_module('numpy')
- def asarray_chkfinite(a, dtype=None, order=None):
- """Convert the input to an array, checking for NaNs or Infs.
- Parameters
- ----------
- a : array_like
- Input data, in any form that can be converted to an array. This
- includes lists, lists of tuples, tuples, tuples of tuples, tuples
- of lists and ndarrays. Success requires no NaNs or Infs.
- dtype : data-type, optional
- By default, the data-type is inferred from the input data.
- order : {'C', 'F', 'A', 'K'}, optional
- Memory layout. 'A' and 'K' depend on the order of input array a.
- 'C' row-major (C-style),
- 'F' column-major (Fortran-style) memory representation.
- 'A' (any) means 'F' if `a` is Fortran contiguous, 'C' otherwise
- 'K' (keep) preserve input order
- Defaults to 'C'.
- Returns
- -------
- out : ndarray
- Array interpretation of `a`. No copy is performed if the input
- is already an ndarray. If `a` is a subclass of ndarray, a base
- class ndarray is returned.
- Raises
- ------
- ValueError
- Raises ValueError if `a` contains NaN (Not a Number) or Inf (Infinity).
- See Also
- --------
- asarray : Create and array.
- asanyarray : Similar function which passes through subclasses.
- ascontiguousarray : Convert input to a contiguous array.
- asfarray : Convert input to a floating point ndarray.
- asfortranarray : Convert input to an ndarray with column-major
- memory order.
- fromiter : Create an array from an iterator.
- fromfunction : Construct an array by executing a function on grid
- positions.
- Examples
- --------
- Convert a list into an array. If all elements are finite
- ``asarray_chkfinite`` is identical to ``asarray``.
- >>> a = [1, 2]
- >>> np.asarray_chkfinite(a, dtype=float)
- array([1., 2.])
- Raises ValueError if array_like contains Nans or Infs.
- >>> a = [1, 2, np.inf]
- >>> try:
- ... np.asarray_chkfinite(a)
- ... except ValueError:
- ... print('ValueError')
- ...
- ValueError
- """
- a = asarray(a, dtype=dtype, order=order)
- if a.dtype.char in typecodes['AllFloat'] and not np.isfinite(a).all():
- raise ValueError(
- "array must not contain infs or NaNs")
- return a
- def _piecewise_dispatcher(x, condlist, funclist, *args, **kw):
- yield x
- # support the undocumented behavior of allowing scalars
- if np.iterable(condlist):
- yield from condlist
- @array_function_dispatch(_piecewise_dispatcher)
- def piecewise(x, condlist, funclist, *args, **kw):
- """
- Evaluate a piecewise-defined function.
- Given a set of conditions and corresponding functions, evaluate each
- function on the input data wherever its condition is true.
- Parameters
- ----------
- x : ndarray or scalar
- The input domain.
- condlist : list of bool arrays or bool scalars
- Each boolean array corresponds to a function in `funclist`. Wherever
- `condlist[i]` is True, `funclist[i](x)` is used as the output value.
- Each boolean array in `condlist` selects a piece of `x`,
- and should therefore be of the same shape as `x`.
- The length of `condlist` must correspond to that of `funclist`.
- If one extra function is given, i.e. if
- ``len(funclist) == len(condlist) + 1``, then that extra function
- is the default value, used wherever all conditions are false.
- funclist : list of callables, f(x,*args,**kw), or scalars
- Each function is evaluated over `x` wherever its corresponding
- condition is True. It should take a 1d array as input and give an 1d
- array or a scalar value as output. If, instead of a callable,
- a scalar is provided then a constant function (``lambda x: scalar``) is
- assumed.
- args : tuple, optional
- Any further arguments given to `piecewise` are passed to the functions
- upon execution, i.e., if called ``piecewise(..., ..., 1, 'a')``, then
- each function is called as ``f(x, 1, 'a')``.
- kw : dict, optional
- Keyword arguments used in calling `piecewise` are passed to the
- functions upon execution, i.e., if called
- ``piecewise(..., ..., alpha=1)``, then each function is called as
- ``f(x, alpha=1)``.
- Returns
- -------
- out : ndarray
- The output is the same shape and type as x and is found by
- calling the functions in `funclist` on the appropriate portions of `x`,
- as defined by the boolean arrays in `condlist`. Portions not covered
- by any condition have a default value of 0.
- See Also
- --------
- choose, select, where
- Notes
- -----
- This is similar to choose or select, except that functions are
- evaluated on elements of `x` that satisfy the corresponding condition from
- `condlist`.
- The result is::
- |--
- |funclist[0](x[condlist[0]])
- out = |funclist[1](x[condlist[1]])
- |...
- |funclist[n2](x[condlist[n2]])
- |--
- Examples
- --------
- Define the sigma function, which is -1 for ``x < 0`` and +1 for ``x >= 0``.
- >>> x = np.linspace(-2.5, 2.5, 6)
- >>> np.piecewise(x, [x < 0, x >= 0], [-1, 1])
- array([-1., -1., -1., 1., 1., 1.])
- Define the absolute value, which is ``-x`` for ``x <0`` and ``x`` for
- ``x >= 0``.
- >>> np.piecewise(x, [x < 0, x >= 0], [lambda x: -x, lambda x: x])
- array([2.5, 1.5, 0.5, 0.5, 1.5, 2.5])
- Apply the same function to a scalar value.
- >>> y = -2
- >>> np.piecewise(y, [y < 0, y >= 0], [lambda x: -x, lambda x: x])
- array(2)
- """
- x = asanyarray(x)
- n2 = len(funclist)
- # undocumented: single condition is promoted to a list of one condition
- if isscalar(condlist) or (
- not isinstance(condlist[0], (list, ndarray)) and x.ndim != 0):
- condlist = [condlist]
- condlist = array(condlist, dtype=bool)
- n = len(condlist)
- if n == n2 - 1: # compute the "otherwise" condition.
- condelse = ~np.any(condlist, axis=0, keepdims=True)
- condlist = np.concatenate([condlist, condelse], axis=0)
- n += 1
- elif n != n2:
- raise ValueError(
- "with {} condition(s), either {} or {} functions are expected"
- .format(n, n, n+1)
- )
- y = zeros(x.shape, x.dtype)
- for cond, func in zip(condlist, funclist):
- if not isinstance(func, collections.abc.Callable):
- y[cond] = func
- else:
- vals = x[cond]
- if vals.size > 0:
- y[cond] = func(vals, *args, **kw)
- return y
- def _select_dispatcher(condlist, choicelist, default=None):
- yield from condlist
- yield from choicelist
- @array_function_dispatch(_select_dispatcher)
- def select(condlist, choicelist, default=0):
- """
- Return an array drawn from elements in choicelist, depending on conditions.
- Parameters
- ----------
- condlist : list of bool ndarrays
- The list of conditions which determine from which array in `choicelist`
- the output elements are taken. When multiple conditions are satisfied,
- the first one encountered in `condlist` is used.
- choicelist : list of ndarrays
- The list of arrays from which the output elements are taken. It has
- to be of the same length as `condlist`.
- default : scalar, optional
- The element inserted in `output` when all conditions evaluate to False.
- Returns
- -------
- output : ndarray
- The output at position m is the m-th element of the array in
- `choicelist` where the m-th element of the corresponding array in
- `condlist` is True.
- See Also
- --------
- where : Return elements from one of two arrays depending on condition.
- take, choose, compress, diag, diagonal
- Examples
- --------
- >>> x = np.arange(10)
- >>> condlist = [x<3, x>5]
- >>> choicelist = [x, x**2]
- >>> np.select(condlist, choicelist)
- array([ 0, 1, 2, ..., 49, 64, 81])
- """
- # Check the size of condlist and choicelist are the same, or abort.
- if len(condlist) != len(choicelist):
- raise ValueError(
- 'list of cases must be same length as list of conditions')
- # Now that the dtype is known, handle the deprecated select([], []) case
- if len(condlist) == 0:
- raise ValueError("select with an empty condition list is not possible")
- choicelist = [np.asarray(choice) for choice in choicelist]
- choicelist.append(np.asarray(default))
- # need to get the result type before broadcasting for correct scalar
- # behaviour
- dtype = np.result_type(*choicelist)
- # Convert conditions to arrays and broadcast conditions and choices
- # as the shape is needed for the result. Doing it separately optimizes
- # for example when all choices are scalars.
- condlist = np.broadcast_arrays(*condlist)
- choicelist = np.broadcast_arrays(*choicelist)
- # If cond array is not an ndarray in boolean format or scalar bool, abort.
- for i, cond in enumerate(condlist):
- if cond.dtype.type is not np.bool_:
- raise TypeError(
- 'invalid entry {} in condlist: should be boolean ndarray'.format(i))
- if choicelist[0].ndim == 0:
- # This may be common, so avoid the call.
- result_shape = condlist[0].shape
- else:
- result_shape = np.broadcast_arrays(condlist[0], choicelist[0])[0].shape
- result = np.full(result_shape, choicelist[-1], dtype)
- # Use np.copyto to burn each choicelist array onto result, using the
- # corresponding condlist as a boolean mask. This is done in reverse
- # order since the first choice should take precedence.
- choicelist = choicelist[-2::-1]
- condlist = condlist[::-1]
- for choice, cond in zip(choicelist, condlist):
- np.copyto(result, choice, where=cond)
- return result
- def _copy_dispatcher(a, order=None, subok=None):
- return (a,)
- @array_function_dispatch(_copy_dispatcher)
- def copy(a, order='K', subok=False):
- """
- Return an array copy of the given object.
- Parameters
- ----------
- a : array_like
- Input data.
- order : {'C', 'F', 'A', 'K'}, optional
- Controls the memory layout of the copy. 'C' means C-order,
- 'F' means F-order, 'A' means 'F' if `a` is Fortran contiguous,
- 'C' otherwise. 'K' means match the layout of `a` as closely
- as possible. (Note that this function and :meth:`ndarray.copy` are very
- similar, but have different default values for their order=
- arguments.)
- subok : bool, optional
- If True, then sub-classes will be passed-through, otherwise the
- returned array will be forced to be a base-class array (defaults to False).
- .. versionadded:: 1.19.0
- Returns
- -------
- arr : ndarray
- Array interpretation of `a`.
- See Also
- --------
- ndarray.copy : Preferred method for creating an array copy
- Notes
- -----
- This is equivalent to:
- >>> np.array(a, copy=True) #doctest: +SKIP
- Examples
- --------
- Create an array x, with a reference y and a copy z:
- >>> x = np.array([1, 2, 3])
- >>> y = x
- >>> z = np.copy(x)
- Note that, when we modify x, y changes, but not z:
- >>> x[0] = 10
- >>> x[0] == y[0]
- True
- >>> x[0] == z[0]
- False
- Note that np.copy is a shallow copy and will not copy object
- elements within arrays. This is mainly important for arrays
- containing Python objects. The new array will contain the
- same object which may lead to surprises if that object can
- be modified (is mutable):
- >>> a = np.array([1, 'm', [2, 3, 4]], dtype=object)
- >>> b = np.copy(a)
- >>> b[2][0] = 10
- >>> a
- array([1, 'm', list([10, 3, 4])], dtype=object)
- To ensure all elements within an ``object`` array are copied,
- use `copy.deepcopy`:
- >>> import copy
- >>> a = np.array([1, 'm', [2, 3, 4]], dtype=object)
- >>> c = copy.deepcopy(a)
- >>> c[2][0] = 10
- >>> c
- array([1, 'm', list([10, 3, 4])], dtype=object)
- >>> a
- array([1, 'm', list([2, 3, 4])], dtype=object)
- """
- return array(a, order=order, subok=subok, copy=True)
- # Basic operations
- def _gradient_dispatcher(f, *varargs, axis=None, edge_order=None):
- yield f
- yield from varargs
- @array_function_dispatch(_gradient_dispatcher)
- def gradient(f, *varargs, axis=None, edge_order=1):
- """
- Return the gradient of an N-dimensional array.
- The gradient is computed using second order accurate central differences
- in the interior points and either first or second order accurate one-sides
- (forward or backwards) differences at the boundaries.
- The returned gradient hence has the same shape as the input array.
- Parameters
- ----------
- f : array_like
- An N-dimensional array containing samples of a scalar function.
- varargs : list of scalar or array, optional
- Spacing between f values. Default unitary spacing for all dimensions.
- Spacing can be specified using:
- 1. single scalar to specify a sample distance for all dimensions.
- 2. N scalars to specify a constant sample distance for each dimension.
- i.e. `dx`, `dy`, `dz`, ...
- 3. N arrays to specify the coordinates of the values along each
- dimension of F. The length of the array must match the size of
- the corresponding dimension
- 4. Any combination of N scalars/arrays with the meaning of 2. and 3.
- If `axis` is given, the number of varargs must equal the number of axes.
- Default: 1.
- edge_order : {1, 2}, optional
- Gradient is calculated using N-th order accurate differences
- at the boundaries. Default: 1.
- .. versionadded:: 1.9.1
- axis : None or int or tuple of ints, optional
- Gradient is calculated only along the given axis or axes
- The default (axis = None) is to calculate the gradient for all the axes
- of the input array. axis may be negative, in which case it counts from
- the last to the first axis.
- .. versionadded:: 1.11.0
- Returns
- -------
- gradient : ndarray or list of ndarray
- A set of ndarrays (or a single ndarray if there is only one dimension)
- corresponding to the derivatives of f with respect to each dimension.
- Each derivative has the same shape as f.
- Examples
- --------
- >>> f = np.array([1, 2, 4, 7, 11, 16], dtype=float)
- >>> np.gradient(f)
- array([1. , 1.5, 2.5, 3.5, 4.5, 5. ])
- >>> np.gradient(f, 2)
- array([0.5 , 0.75, 1.25, 1.75, 2.25, 2.5 ])
- Spacing can be also specified with an array that represents the coordinates
- of the values F along the dimensions.
- For instance a uniform spacing:
- >>> x = np.arange(f.size)
- >>> np.gradient(f, x)
- array([1. , 1.5, 2.5, 3.5, 4.5, 5. ])
- Or a non uniform one:
- >>> x = np.array([0., 1., 1.5, 3.5, 4., 6.], dtype=float)
- >>> np.gradient(f, x)
- array([1. , 3. , 3.5, 6.7, 6.9, 2.5])
- For two dimensional arrays, the return will be two arrays ordered by
- axis. In this example the first array stands for the gradient in
- rows and the second one in columns direction:
- >>> np.gradient(np.array([[1, 2, 6], [3, 4, 5]], dtype=float))
- [array([[ 2., 2., -1.],
- [ 2., 2., -1.]]), array([[1. , 2.5, 4. ],
- [1. , 1. , 1. ]])]
- In this example the spacing is also specified:
- uniform for axis=0 and non uniform for axis=1
- >>> dx = 2.
- >>> y = [1., 1.5, 3.5]
- >>> np.gradient(np.array([[1, 2, 6], [3, 4, 5]], dtype=float), dx, y)
- [array([[ 1. , 1. , -0.5],
- [ 1. , 1. , -0.5]]), array([[2. , 2. , 2. ],
- [2. , 1.7, 0.5]])]
- It is possible to specify how boundaries are treated using `edge_order`
- >>> x = np.array([0, 1, 2, 3, 4])
- >>> f = x**2
- >>> np.gradient(f, edge_order=1)
- array([1., 2., 4., 6., 7.])
- >>> np.gradient(f, edge_order=2)
- array([0., 2., 4., 6., 8.])
- The `axis` keyword can be used to specify a subset of axes of which the
- gradient is calculated
- >>> np.gradient(np.array([[1, 2, 6], [3, 4, 5]], dtype=float), axis=0)
- array([[ 2., 2., -1.],
- [ 2., 2., -1.]])
- Notes
- -----
- Assuming that :math:`f\\in C^{3}` (i.e., :math:`f` has at least 3 continuous
- derivatives) and let :math:`h_{*}` be a non-homogeneous stepsize, we
- minimize the "consistency error" :math:`\\eta_{i}` between the true gradient
- and its estimate from a linear combination of the neighboring grid-points:
- .. math::
- \\eta_{i} = f_{i}^{\\left(1\\right)} -
- \\left[ \\alpha f\\left(x_{i}\\right) +
- \\beta f\\left(x_{i} + h_{d}\\right) +
- \\gamma f\\left(x_{i}-h_{s}\\right)
- \\right]
- By substituting :math:`f(x_{i} + h_{d})` and :math:`f(x_{i} - h_{s})`
- with their Taylor series expansion, this translates into solving
- the following the linear system:
- .. math::
- \\left\\{
- \\begin{array}{r}
- \\alpha+\\beta+\\gamma=0 \\\\
- \\beta h_{d}-\\gamma h_{s}=1 \\\\
- \\beta h_{d}^{2}+\\gamma h_{s}^{2}=0
- \\end{array}
- \\right.
- The resulting approximation of :math:`f_{i}^{(1)}` is the following:
- .. math::
- \\hat f_{i}^{(1)} =
- \\frac{
- h_{s}^{2}f\\left(x_{i} + h_{d}\\right)
- + \\left(h_{d}^{2} - h_{s}^{2}\\right)f\\left(x_{i}\\right)
- - h_{d}^{2}f\\left(x_{i}-h_{s}\\right)}
- { h_{s}h_{d}\\left(h_{d} + h_{s}\\right)}
- + \\mathcal{O}\\left(\\frac{h_{d}h_{s}^{2}
- + h_{s}h_{d}^{2}}{h_{d}
- + h_{s}}\\right)
- It is worth noting that if :math:`h_{s}=h_{d}`
- (i.e., data are evenly spaced)
- we find the standard second order approximation:
- .. math::
- \\hat f_{i}^{(1)}=
- \\frac{f\\left(x_{i+1}\\right) - f\\left(x_{i-1}\\right)}{2h}
- + \\mathcal{O}\\left(h^{2}\\right)
- With a similar procedure the forward/backward approximations used for
- boundaries can be derived.
- References
- ----------
- .. [1] Quarteroni A., Sacco R., Saleri F. (2007) Numerical Mathematics
- (Texts in Applied Mathematics). New York: Springer.
- .. [2] Durran D. R. (1999) Numerical Methods for Wave Equations
- in Geophysical Fluid Dynamics. New York: Springer.
- .. [3] Fornberg B. (1988) Generation of Finite Difference Formulas on
- Arbitrarily Spaced Grids,
- Mathematics of Computation 51, no. 184 : 699-706.
- `PDF <http://www.ams.org/journals/mcom/1988-51-184/
- S0025-5718-1988-0935077-0/S0025-5718-1988-0935077-0.pdf>`_.
- """
- f = np.asanyarray(f)
- N = f.ndim # number of dimensions
- if axis is None:
- axes = tuple(range(N))
- else:
- axes = _nx.normalize_axis_tuple(axis, N)
- len_axes = len(axes)
- n = len(varargs)
- if n == 0:
- # no spacing argument - use 1 in all axes
- dx = [1.0] * len_axes
- elif n == 1 and np.ndim(varargs[0]) == 0:
- # single scalar for all axes
- dx = varargs * len_axes
- elif n == len_axes:
- # scalar or 1d array for each axis
- dx = list(varargs)
- for i, distances in enumerate(dx):
- distances = np.asanyarray(distances)
- if distances.ndim == 0:
- continue
- elif distances.ndim != 1:
- raise ValueError("distances must be either scalars or 1d")
- if len(distances) != f.shape[axes[i]]:
- raise ValueError("when 1d, distances must match "
- "the length of the corresponding dimension")
- if np.issubdtype(distances.dtype, np.integer):
- # Convert numpy integer types to float64 to avoid modular
- # arithmetic in np.diff(distances).
- distances = distances.astype(np.float64)
- diffx = np.diff(distances)
- # if distances are constant reduce to the scalar case
- # since it brings a consistent speedup
- if (diffx == diffx[0]).all():
- diffx = diffx[0]
- dx[i] = diffx
- else:
- raise TypeError("invalid number of arguments")
- if edge_order > 2:
- raise ValueError("'edge_order' greater than 2 not supported")
- # use central differences on interior and one-sided differences on the
- # endpoints. This preserves second order-accuracy over the full domain.
- outvals = []
- # create slice objects --- initially all are [:, :, ..., :]
- slice1 = [slice(None)]*N
- slice2 = [slice(None)]*N
- slice3 = [slice(None)]*N
- slice4 = [slice(None)]*N
- otype = f.dtype
- if otype.type is np.datetime64:
- # the timedelta dtype with the same unit information
- otype = np.dtype(otype.name.replace('datetime', 'timedelta'))
- # view as timedelta to allow addition
- f = f.view(otype)
- elif otype.type is np.timedelta64:
- pass
- elif np.issubdtype(otype, np.inexact):
- pass
- else:
- # All other types convert to floating point.
- # First check if f is a numpy integer type; if so, convert f to float64
- # to avoid modular arithmetic when computing the changes in f.
- if np.issubdtype(otype, np.integer):
- f = f.astype(np.float64)
- otype = np.float64
- for axis, ax_dx in zip(axes, dx):
- if f.shape[axis] < edge_order + 1:
- raise ValueError(
- "Shape of array too small to calculate a numerical gradient, "
- "at least (edge_order + 1) elements are required.")
- # result allocation
- out = np.empty_like(f, dtype=otype)
- # spacing for the current axis
- uniform_spacing = np.ndim(ax_dx) == 0
- # Numerical differentiation: 2nd order interior
- slice1[axis] = slice(1, -1)
- slice2[axis] = slice(None, -2)
- slice3[axis] = slice(1, -1)
- slice4[axis] = slice(2, None)
- if uniform_spacing:
- out[tuple(slice1)] = (f[tuple(slice4)] - f[tuple(slice2)]) / (2. * ax_dx)
- else:
- dx1 = ax_dx[0:-1]
- dx2 = ax_dx[1:]
- a = -(dx2)/(dx1 * (dx1 + dx2))
- b = (dx2 - dx1) / (dx1 * dx2)
- c = dx1 / (dx2 * (dx1 + dx2))
- # fix the shape for broadcasting
- shape = np.ones(N, dtype=int)
- shape[axis] = -1
- a.shape = b.shape = c.shape = shape
- # 1D equivalent -- out[1:-1] = a * f[:-2] + b * f[1:-1] + c * f[2:]
- out[tuple(slice1)] = a * f[tuple(slice2)] + b * f[tuple(slice3)] + c * f[tuple(slice4)]
- # Numerical differentiation: 1st order edges
- if edge_order == 1:
- slice1[axis] = 0
- slice2[axis] = 1
- slice3[axis] = 0
- dx_0 = ax_dx if uniform_spacing else ax_dx[0]
- # 1D equivalent -- out[0] = (f[1] - f[0]) / (x[1] - x[0])
- out[tuple(slice1)] = (f[tuple(slice2)] - f[tuple(slice3)]) / dx_0
- slice1[axis] = -1
- slice2[axis] = -1
- slice3[axis] = -2
- dx_n = ax_dx if uniform_spacing else ax_dx[-1]
- # 1D equivalent -- out[-1] = (f[-1] - f[-2]) / (x[-1] - x[-2])
- out[tuple(slice1)] = (f[tuple(slice2)] - f[tuple(slice3)]) / dx_n
- # Numerical differentiation: 2nd order edges
- else:
- slice1[axis] = 0
- slice2[axis] = 0
- slice3[axis] = 1
- slice4[axis] = 2
- if uniform_spacing:
- a = -1.5 / ax_dx
- b = 2. / ax_dx
- c = -0.5 / ax_dx
- else:
- dx1 = ax_dx[0]
- dx2 = ax_dx[1]
- a = -(2. * dx1 + dx2)/(dx1 * (dx1 + dx2))
- b = (dx1 + dx2) / (dx1 * dx2)
- c = - dx1 / (dx2 * (dx1 + dx2))
- # 1D equivalent -- out[0] = a * f[0] + b * f[1] + c * f[2]
- out[tuple(slice1)] = a * f[tuple(slice2)] + b * f[tuple(slice3)] + c * f[tuple(slice4)]
- slice1[axis] = -1
- slice2[axis] = -3
- slice3[axis] = -2
- slice4[axis] = -1
- if uniform_spacing:
- a = 0.5 / ax_dx
- b = -2. / ax_dx
- c = 1.5 / ax_dx
- else:
- dx1 = ax_dx[-2]
- dx2 = ax_dx[-1]
- a = (dx2) / (dx1 * (dx1 + dx2))
- b = - (dx2 + dx1) / (dx1 * dx2)
- c = (2. * dx2 + dx1) / (dx2 * (dx1 + dx2))
- # 1D equivalent -- out[-1] = a * f[-3] + b * f[-2] + c * f[-1]
- out[tuple(slice1)] = a * f[tuple(slice2)] + b * f[tuple(slice3)] + c * f[tuple(slice4)]
- outvals.append(out)
- # reset the slice object in this dimension to ":"
- slice1[axis] = slice(None)
- slice2[axis] = slice(None)
- slice3[axis] = slice(None)
- slice4[axis] = slice(None)
- if len_axes == 1:
- return outvals[0]
- else:
- return outvals
- def _diff_dispatcher(a, n=None, axis=None, prepend=None, append=None):
- return (a, prepend, append)
- @array_function_dispatch(_diff_dispatcher)
- def diff(a, n=1, axis=-1, prepend=np._NoValue, append=np._NoValue):
- """
- Calculate the n-th discrete difference along the given axis.
- The first difference is given by ``out[i] = a[i+1] - a[i]`` along
- the given axis, higher differences are calculated by using `diff`
- recursively.
- Parameters
- ----------
- a : array_like
- Input array
- n : int, optional
- The number of times values are differenced. If zero, the input
- is returned as-is.
- axis : int, optional
- The axis along which the difference is taken, default is the
- last axis.
- prepend, append : array_like, optional
- Values to prepend or append to `a` along axis prior to
- performing the difference. Scalar values are expanded to
- arrays with length 1 in the direction of axis and the shape
- of the input array in along all other axes. Otherwise the
- dimension and shape must match `a` except along axis.
- .. versionadded:: 1.16.0
- Returns
- -------
- diff : ndarray
- The n-th differences. The shape of the output is the same as `a`
- except along `axis` where the dimension is smaller by `n`. The
- type of the output is the same as the type of the difference
- between any two elements of `a`. This is the same as the type of
- `a` in most cases. A notable exception is `datetime64`, which
- results in a `timedelta64` output array.
- See Also
- --------
- gradient, ediff1d, cumsum
- Notes
- -----
- Type is preserved for boolean arrays, so the result will contain
- `False` when consecutive elements are the same and `True` when they
- differ.
- For unsigned integer arrays, the results will also be unsigned. This
- should not be surprising, as the result is consistent with
- calculating the difference directly:
- >>> u8_arr = np.array([1, 0], dtype=np.uint8)
- >>> np.diff(u8_arr)
- array([255], dtype=uint8)
- >>> u8_arr[1,...] - u8_arr[0,...]
- 255
- If this is not desirable, then the array should be cast to a larger
- integer type first:
- >>> i16_arr = u8_arr.astype(np.int16)
- >>> np.diff(i16_arr)
- array([-1], dtype=int16)
- Examples
- --------
- >>> x = np.array([1, 2, 4, 7, 0])
- >>> np.diff(x)
- array([ 1, 2, 3, -7])
- >>> np.diff(x, n=2)
- array([ 1, 1, -10])
- >>> x = np.array([[1, 3, 6, 10], [0, 5, 6, 8]])
- >>> np.diff(x)
- array([[2, 3, 4],
- [5, 1, 2]])
- >>> np.diff(x, axis=0)
- array([[-1, 2, 0, -2]])
- >>> x = np.arange('1066-10-13', '1066-10-16', dtype=np.datetime64)
- >>> np.diff(x)
- array([1, 1], dtype='timedelta64[D]')
- """
- if n == 0:
- return a
- if n < 0:
- raise ValueError(
- "order must be non-negative but got " + repr(n))
- a = asanyarray(a)
- nd = a.ndim
- if nd == 0:
- raise ValueError("diff requires input that is at least one dimensional")
- axis = normalize_axis_index(axis, nd)
- combined = []
- if prepend is not np._NoValue:
- prepend = np.asanyarray(prepend)
- if prepend.ndim == 0:
- shape = list(a.shape)
- shape[axis] = 1
- prepend = np.broadcast_to(prepend, tuple(shape))
- combined.append(prepend)
- combined.append(a)
- if append is not np._NoValue:
- append = np.asanyarray(append)
- if append.ndim == 0:
- shape = list(a.shape)
- shape[axis] = 1
- append = np.broadcast_to(append, tuple(shape))
- combined.append(append)
- if len(combined) > 1:
- a = np.concatenate(combined, axis)
- slice1 = [slice(None)] * nd
- slice2 = [slice(None)] * nd
- slice1[axis] = slice(1, None)
- slice2[axis] = slice(None, -1)
- slice1 = tuple(slice1)
- slice2 = tuple(slice2)
- op = not_equal if a.dtype == np.bool_ else subtract
- for _ in range(n):
- a = op(a[slice1], a[slice2])
- return a
- def _interp_dispatcher(x, xp, fp, left=None, right=None, period=None):
- return (x, xp, fp)
- @array_function_dispatch(_interp_dispatcher)
- def interp(x, xp, fp, left=None, right=None, period=None):
- """
- One-dimensional linear interpolation.
- Returns the one-dimensional piecewise linear interpolant to a function
- with given discrete data points (`xp`, `fp`), evaluated at `x`.
- Parameters
- ----------
- x : array_like
- The x-coordinates at which to evaluate the interpolated values.
- xp : 1-D sequence of floats
- The x-coordinates of the data points, must be increasing if argument
- `period` is not specified. Otherwise, `xp` is internally sorted after
- normalizing the periodic boundaries with ``xp = xp % period``.
- fp : 1-D sequence of float or complex
- The y-coordinates of the data points, same length as `xp`.
- left : optional float or complex corresponding to fp
- Value to return for `x < xp[0]`, default is `fp[0]`.
- right : optional float or complex corresponding to fp
- Value to return for `x > xp[-1]`, default is `fp[-1]`.
- period : None or float, optional
- A period for the x-coordinates. This parameter allows the proper
- interpolation of angular x-coordinates. Parameters `left` and `right`
- are ignored if `period` is specified.
- .. versionadded:: 1.10.0
- Returns
- -------
- y : float or complex (corresponding to fp) or ndarray
- The interpolated values, same shape as `x`.
- Raises
- ------
- ValueError
- If `xp` and `fp` have different length
- If `xp` or `fp` are not 1-D sequences
- If `period == 0`
- See Also
- --------
- scipy.interpolate
- Notes
- -----
- The x-coordinate sequence is expected to be increasing, but this is not
- explicitly enforced. However, if the sequence `xp` is non-increasing,
- interpolation results are meaningless.
- Note that, since NaN is unsortable, `xp` also cannot contain NaNs.
- A simple check for `xp` being strictly increasing is::
- np.all(np.diff(xp) > 0)
- Examples
- --------
- >>> xp = [1, 2, 3]
- >>> fp = [3, 2, 0]
- >>> np.interp(2.5, xp, fp)
- 1.0
- >>> np.interp([0, 1, 1.5, 2.72, 3.14], xp, fp)
- array([3. , 3. , 2.5 , 0.56, 0. ])
- >>> UNDEF = -99.0
- >>> np.interp(3.14, xp, fp, right=UNDEF)
- -99.0
- Plot an interpolant to the sine function:
- >>> x = np.linspace(0, 2*np.pi, 10)
- >>> y = np.sin(x)
- >>> xvals = np.linspace(0, 2*np.pi, 50)
- >>> yinterp = np.interp(xvals, x, y)
- >>> import matplotlib.pyplot as plt
- >>> plt.plot(x, y, 'o')
- [<matplotlib.lines.Line2D object at 0x...>]
- >>> plt.plot(xvals, yinterp, '-x')
- [<matplotlib.lines.Line2D object at 0x...>]
- >>> plt.show()
- Interpolation with periodic x-coordinates:
- >>> x = [-180, -170, -185, 185, -10, -5, 0, 365]
- >>> xp = [190, -190, 350, -350]
- >>> fp = [5, 10, 3, 4]
- >>> np.interp(x, xp, fp, period=360)
- array([7.5 , 5. , 8.75, 6.25, 3. , 3.25, 3.5 , 3.75])
- Complex interpolation:
- >>> x = [1.5, 4.0]
- >>> xp = [2,3,5]
- >>> fp = [1.0j, 0, 2+3j]
- >>> np.interp(x, xp, fp)
- array([0.+1.j , 1.+1.5j])
- """
- fp = np.asarray(fp)
- if np.iscomplexobj(fp):
- interp_func = compiled_interp_complex
- input_dtype = np.complex128
- else:
- interp_func = compiled_interp
- input_dtype = np.float64
- if period is not None:
- if period == 0:
- raise ValueError("period must be a non-zero value")
- period = abs(period)
- left = None
- right = None
- x = np.asarray(x, dtype=np.float64)
- xp = np.asarray(xp, dtype=np.float64)
- fp = np.asarray(fp, dtype=input_dtype)
- if xp.ndim != 1 or fp.ndim != 1:
- raise ValueError("Data points must be 1-D sequences")
- if xp.shape[0] != fp.shape[0]:
- raise ValueError("fp and xp are not of the same length")
- # normalizing periodic boundaries
- x = x % period
- xp = xp % period
- asort_xp = np.argsort(xp)
- xp = xp[asort_xp]
- fp = fp[asort_xp]
- xp = np.concatenate((xp[-1:]-period, xp, xp[0:1]+period))
- fp = np.concatenate((fp[-1:], fp, fp[0:1]))
- return interp_func(x, xp, fp, left, right)
- def _angle_dispatcher(z, deg=None):
- return (z,)
- @array_function_dispatch(_angle_dispatcher)
- def angle(z, deg=False):
- """
- Return the angle of the complex argument.
- Parameters
- ----------
- z : array_like
- A complex number or sequence of complex numbers.
- deg : bool, optional
- Return angle in degrees if True, radians if False (default).
- Returns
- -------
- angle : ndarray or scalar
- The counterclockwise angle from the positive real axis on the complex
- plane in the range ``(-pi, pi]``, with dtype as numpy.float64.
- .. versionchanged:: 1.16.0
- This function works on subclasses of ndarray like `ma.array`.
- See Also
- --------
- arctan2
- absolute
- Notes
- -----
- Although the angle of the complex number 0 is undefined, ``numpy.angle(0)``
- returns the value 0.
- Examples
- --------
- >>> np.angle([1.0, 1.0j, 1+1j]) # in radians
- array([ 0. , 1.57079633, 0.78539816]) # may vary
- >>> np.angle(1+1j, deg=True) # in degrees
- 45.0
- """
- z = asanyarray(z)
- if issubclass(z.dtype.type, _nx.complexfloating):
- zimag = z.imag
- zreal = z.real
- else:
- zimag = 0
- zreal = z
- a = arctan2(zimag, zreal)
- if deg:
- a *= 180/pi
- return a
- def _unwrap_dispatcher(p, discont=None, axis=None):
- return (p,)
- @array_function_dispatch(_unwrap_dispatcher)
- def unwrap(p, discont=pi, axis=-1):
- """
- Unwrap by changing deltas between values to 2*pi complement.
- Unwrap radian phase `p` by changing absolute jumps greater than
- `discont` to their 2*pi complement along the given axis.
- Parameters
- ----------
- p : array_like
- Input array.
- discont : float, optional
- Maximum discontinuity between values, default is ``pi``.
- axis : int, optional
- Axis along which unwrap will operate, default is the last axis.
- Returns
- -------
- out : ndarray
- Output array.
- See Also
- --------
- rad2deg, deg2rad
- Notes
- -----
- If the discontinuity in `p` is smaller than ``pi``, but larger than
- `discont`, no unwrapping is done because taking the 2*pi complement
- would only make the discontinuity larger.
- Examples
- --------
- >>> phase = np.linspace(0, np.pi, num=5)
- >>> phase[3:] += np.pi
- >>> phase
- array([ 0. , 0.78539816, 1.57079633, 5.49778714, 6.28318531]) # may vary
- >>> np.unwrap(phase)
- array([ 0. , 0.78539816, 1.57079633, -0.78539816, 0. ]) # may vary
- """
- p = asarray(p)
- nd = p.ndim
- dd = diff(p, axis=axis)
- slice1 = [slice(None, None)]*nd # full slices
- slice1[axis] = slice(1, None)
- slice1 = tuple(slice1)
- ddmod = mod(dd + pi, 2*pi) - pi
- _nx.copyto(ddmod, pi, where=(ddmod == -pi) & (dd > 0))
- ph_correct = ddmod - dd
- _nx.copyto(ph_correct, 0, where=abs(dd) < discont)
- up = array(p, copy=True, dtype='d')
- up[slice1] = p[slice1] + ph_correct.cumsum(axis)
- return up
- def _sort_complex(a):
- return (a,)
- @array_function_dispatch(_sort_complex)
- def sort_complex(a):
- """
- Sort a complex array using the real part first, then the imaginary part.
- Parameters
- ----------
- a : array_like
- Input array
- Returns
- -------
- out : complex ndarray
- Always returns a sorted complex array.
- Examples
- --------
- >>> np.sort_complex([5, 3, 6, 2, 1])
- array([1.+0.j, 2.+0.j, 3.+0.j, 5.+0.j, 6.+0.j])
- >>> np.sort_complex([1 + 2j, 2 - 1j, 3 - 2j, 3 - 3j, 3 + 5j])
- array([1.+2.j, 2.-1.j, 3.-3.j, 3.-2.j, 3.+5.j])
- """
- b = array(a, copy=True)
- b.sort()
- if not issubclass(b.dtype.type, _nx.complexfloating):
- if b.dtype.char in 'bhBH':
- return b.astype('F')
- elif b.dtype.char == 'g':
- return b.astype('G')
- else:
- return b.astype('D')
- else:
- return b
- def _trim_zeros(filt, trim=None):
- return (filt,)
- @array_function_dispatch(_trim_zeros)
- def trim_zeros(filt, trim='fb'):
- """
- Trim the leading and/or trailing zeros from a 1-D array or sequence.
- Parameters
- ----------
- filt : 1-D array or sequence
- Input array.
- trim : str, optional
- A string with 'f' representing trim from front and 'b' to trim from
- back. Default is 'fb', trim zeros from both front and back of the
- array.
- Returns
- -------
- trimmed : 1-D array or sequence
- The result of trimming the input. The input data type is preserved.
- Examples
- --------
- >>> a = np.array((0, 0, 0, 1, 2, 3, 0, 2, 1, 0))
- >>> np.trim_zeros(a)
- array([1, 2, 3, 0, 2, 1])
- >>> np.trim_zeros(a, 'b')
- array([0, 0, 0, ..., 0, 2, 1])
- The input data type is preserved, list/tuple in means list/tuple out.
- >>> np.trim_zeros([0, 1, 2, 0])
- [1, 2]
- """
- first = 0
- trim = trim.upper()
- if 'F' in trim:
- for i in filt:
- if i != 0.:
- break
- else:
- first = first + 1
- last = len(filt)
- if 'B' in trim:
- for i in filt[::-1]:
- if i != 0.:
- break
- else:
- last = last - 1
- return filt[first:last]
- def _extract_dispatcher(condition, arr):
- return (condition, arr)
- @array_function_dispatch(_extract_dispatcher)
- def extract(condition, arr):
- """
- Return the elements of an array that satisfy some condition.
- This is equivalent to ``np.compress(ravel(condition), ravel(arr))``. If
- `condition` is boolean ``np.extract`` is equivalent to ``arr[condition]``.
- Note that `place` does the exact opposite of `extract`.
- Parameters
- ----------
- condition : array_like
- An array whose nonzero or True entries indicate the elements of `arr`
- to extract.
- arr : array_like
- Input array of the same size as `condition`.
- Returns
- -------
- extract : ndarray
- Rank 1 array of values from `arr` where `condition` is True.
- See Also
- --------
- take, put, copyto, compress, place
- Examples
- --------
- >>> arr = np.arange(12).reshape((3, 4))
- >>> arr
- array([[ 0, 1, 2, 3],
- [ 4, 5, 6, 7],
- [ 8, 9, 10, 11]])
- >>> condition = np.mod(arr, 3)==0
- >>> condition
- array([[ True, False, False, True],
- [False, False, True, False],
- [False, True, False, False]])
- >>> np.extract(condition, arr)
- array([0, 3, 6, 9])
- If `condition` is boolean:
- >>> arr[condition]
- array([0, 3, 6, 9])
- """
- return _nx.take(ravel(arr), nonzero(ravel(condition))[0])
- def _place_dispatcher(arr, mask, vals):
- return (arr, mask, vals)
- @array_function_dispatch(_place_dispatcher)
- def place(arr, mask, vals):
- """
- Change elements of an array based on conditional and input values.
- Similar to ``np.copyto(arr, vals, where=mask)``, the difference is that
- `place` uses the first N elements of `vals`, where N is the number of
- True values in `mask`, while `copyto` uses the elements where `mask`
- is True.
- Note that `extract` does the exact opposite of `place`.
- Parameters
- ----------
- arr : ndarray
- Array to put data into.
- mask : array_like
- Boolean mask array. Must have the same size as `a`.
- vals : 1-D sequence
- Values to put into `a`. Only the first N elements are used, where
- N is the number of True values in `mask`. If `vals` is smaller
- than N, it will be repeated, and if elements of `a` are to be masked,
- this sequence must be non-empty.
- See Also
- --------
- copyto, put, take, extract
- Examples
- --------
- >>> arr = np.arange(6).reshape(2, 3)
- >>> np.place(arr, arr>2, [44, 55])
- >>> arr
- array([[ 0, 1, 2],
- [44, 55, 44]])
- """
- if not isinstance(arr, np.ndarray):
- raise TypeError("argument 1 must be numpy.ndarray, "
- "not {name}".format(name=type(arr).__name__))
- return _insert(arr, mask, vals)
- def disp(mesg, device=None, linefeed=True):
- """
- Display a message on a device.
- Parameters
- ----------
- mesg : str
- Message to display.
- device : object
- Device to write message. If None, defaults to ``sys.stdout`` which is
- very similar to ``print``. `device` needs to have ``write()`` and
- ``flush()`` methods.
- linefeed : bool, optional
- Option whether to print a line feed or not. Defaults to True.
- Raises
- ------
- AttributeError
- If `device` does not have a ``write()`` or ``flush()`` method.
- Examples
- --------
- Besides ``sys.stdout``, a file-like object can also be used as it has
- both required methods:
- >>> from io import StringIO
- >>> buf = StringIO()
- >>> np.disp(u'"Display" in a file', device=buf)
- >>> buf.getvalue()
- '"Display" in a file\\n'
- """
- if device is None:
- device = sys.stdout
- if linefeed:
- device.write('%s\n' % mesg)
- else:
- device.write('%s' % mesg)
- device.flush()
- return
- # See https://docs.scipy.org/doc/numpy/reference/c-api.generalized-ufuncs.html
- _DIMENSION_NAME = r'\w+'
- _CORE_DIMENSION_LIST = '(?:{0:}(?:,{0:})*)?'.format(_DIMENSION_NAME)
- _ARGUMENT = r'\({}\)'.format(_CORE_DIMENSION_LIST)
- _ARGUMENT_LIST = '{0:}(?:,{0:})*'.format(_ARGUMENT)
- _SIGNATURE = '^{0:}->{0:}$'.format(_ARGUMENT_LIST)
- def _parse_gufunc_signature(signature):
- """
- Parse string signatures for a generalized universal function.
- Arguments
- ---------
- signature : string
- Generalized universal function signature, e.g., ``(m,n),(n,p)->(m,p)``
- for ``np.matmul``.
- Returns
- -------
- Tuple of input and output core dimensions parsed from the signature, each
- of the form List[Tuple[str, ...]].
- """
- if not re.match(_SIGNATURE, signature):
- raise ValueError(
- 'not a valid gufunc signature: {}'.format(signature))
- return tuple([tuple(re.findall(_DIMENSION_NAME, arg))
- for arg in re.findall(_ARGUMENT, arg_list)]
- for arg_list in signature.split('->'))
- def _update_dim_sizes(dim_sizes, arg, core_dims):
- """
- Incrementally check and update core dimension sizes for a single argument.
- Arguments
- ---------
- dim_sizes : Dict[str, int]
- Sizes of existing core dimensions. Will be updated in-place.
- arg : ndarray
- Argument to examine.
- core_dims : Tuple[str, ...]
- Core dimensions for this argument.
- """
- if not core_dims:
- return
- num_core_dims = len(core_dims)
- if arg.ndim < num_core_dims:
- raise ValueError(
- '%d-dimensional argument does not have enough '
- 'dimensions for all core dimensions %r'
- % (arg.ndim, core_dims))
- core_shape = arg.shape[-num_core_dims:]
- for dim, size in zip(core_dims, core_shape):
- if dim in dim_sizes:
- if size != dim_sizes[dim]:
- raise ValueError(
- 'inconsistent size for core dimension %r: %r vs %r'
- % (dim, size, dim_sizes[dim]))
- else:
- dim_sizes[dim] = size
- def _parse_input_dimensions(args, input_core_dims):
- """
- Parse broadcast and core dimensions for vectorize with a signature.
- Arguments
- ---------
- args : Tuple[ndarray, ...]
- Tuple of input arguments to examine.
- input_core_dims : List[Tuple[str, ...]]
- List of core dimensions corresponding to each input.
- Returns
- -------
- broadcast_shape : Tuple[int, ...]
- Common shape to broadcast all non-core dimensions to.
- dim_sizes : Dict[str, int]
- Common sizes for named core dimensions.
- """
- broadcast_args = []
- dim_sizes = {}
- for arg, core_dims in zip(args, input_core_dims):
- _update_dim_sizes(dim_sizes, arg, core_dims)
- ndim = arg.ndim - len(core_dims)
- dummy_array = np.lib.stride_tricks.as_strided(0, arg.shape[:ndim])
- broadcast_args.append(dummy_array)
- broadcast_shape = np.lib.stride_tricks._broadcast_shape(*broadcast_args)
- return broadcast_shape, dim_sizes
- def _calculate_shapes(broadcast_shape, dim_sizes, list_of_core_dims):
- """Helper for calculating broadcast shapes with core dimensions."""
- return [broadcast_shape + tuple(dim_sizes[dim] for dim in core_dims)
- for core_dims in list_of_core_dims]
- def _create_arrays(broadcast_shape, dim_sizes, list_of_core_dims, dtypes):
- """Helper for creating output arrays in vectorize."""
- shapes = _calculate_shapes(broadcast_shape, dim_sizes, list_of_core_dims)
- arrays = tuple(np.empty(shape, dtype=dtype)
- for shape, dtype in zip(shapes, dtypes))
- return arrays
- @set_module('numpy')
- class vectorize:
- """
- vectorize(pyfunc, otypes=None, doc=None, excluded=None, cache=False,
- signature=None)
- Generalized function class.
- Define a vectorized function which takes a nested sequence of objects or
- numpy arrays as inputs and returns a single numpy array or a tuple of numpy
- arrays. The vectorized function evaluates `pyfunc` over successive tuples
- of the input arrays like the python map function, except it uses the
- broadcasting rules of numpy.
- The data type of the output of `vectorized` is determined by calling
- the function with the first element of the input. This can be avoided
- by specifying the `otypes` argument.
- Parameters
- ----------
- pyfunc : callable
- A python function or method.
- otypes : str or list of dtypes, optional
- The output data type. It must be specified as either a string of
- typecode characters or a list of data type specifiers. There should
- be one data type specifier for each output.
- doc : str, optional
- The docstring for the function. If None, the docstring will be the
- ``pyfunc.__doc__``.
- excluded : set, optional
- Set of strings or integers representing the positional or keyword
- arguments for which the function will not be vectorized. These will be
- passed directly to `pyfunc` unmodified.
- .. versionadded:: 1.7.0
- cache : bool, optional
- If `True`, then cache the first function call that determines the number
- of outputs if `otypes` is not provided.
- .. versionadded:: 1.7.0
- signature : string, optional
- Generalized universal function signature, e.g., ``(m,n),(n)->(m)`` for
- vectorized matrix-vector multiplication. If provided, ``pyfunc`` will
- be called with (and expected to return) arrays with shapes given by the
- size of corresponding core dimensions. By default, ``pyfunc`` is
- assumed to take scalars as input and output.
- .. versionadded:: 1.12.0
- Returns
- -------
- vectorized : callable
- Vectorized function.
- See Also
- --------
- frompyfunc : Takes an arbitrary Python function and returns a ufunc
- Notes
- -----
- The `vectorize` function is provided primarily for convenience, not for
- performance. The implementation is essentially a for loop.
- If `otypes` is not specified, then a call to the function with the
- first argument will be used to determine the number of outputs. The
- results of this call will be cached if `cache` is `True` to prevent
- calling the function twice. However, to implement the cache, the
- original function must be wrapped which will slow down subsequent
- calls, so only do this if your function is expensive.
- The new keyword argument interface and `excluded` argument support
- further degrades performance.
- References
- ----------
- .. [1] :doc:`/reference/c-api/generalized-ufuncs`
- Examples
- --------
- >>> def myfunc(a, b):
- ... "Return a-b if a>b, otherwise return a+b"
- ... if a > b:
- ... return a - b
- ... else:
- ... return a + b
- >>> vfunc = np.vectorize(myfunc)
- >>> vfunc([1, 2, 3, 4], 2)
- array([3, 4, 1, 2])
- The docstring is taken from the input function to `vectorize` unless it
- is specified:
- >>> vfunc.__doc__
- 'Return a-b if a>b, otherwise return a+b'
- >>> vfunc = np.vectorize(myfunc, doc='Vectorized `myfunc`')
- >>> vfunc.__doc__
- 'Vectorized `myfunc`'
- The output type is determined by evaluating the first element of the input,
- unless it is specified:
- >>> out = vfunc([1, 2, 3, 4], 2)
- >>> type(out[0])
- <class 'numpy.int64'>
- >>> vfunc = np.vectorize(myfunc, otypes=[float])
- >>> out = vfunc([1, 2, 3, 4], 2)
- >>> type(out[0])
- <class 'numpy.float64'>
- The `excluded` argument can be used to prevent vectorizing over certain
- arguments. This can be useful for array-like arguments of a fixed length
- such as the coefficients for a polynomial as in `polyval`:
- >>> def mypolyval(p, x):
- ... _p = list(p)
- ... res = _p.pop(0)
- ... while _p:
- ... res = res*x + _p.pop(0)
- ... return res
- >>> vpolyval = np.vectorize(mypolyval, excluded=['p'])
- >>> vpolyval(p=[1, 2, 3], x=[0, 1])
- array([3, 6])
- Positional arguments may also be excluded by specifying their position:
- >>> vpolyval.excluded.add(0)
- >>> vpolyval([1, 2, 3], x=[0, 1])
- array([3, 6])
- The `signature` argument allows for vectorizing functions that act on
- non-scalar arrays of fixed length. For example, you can use it for a
- vectorized calculation of Pearson correlation coefficient and its p-value:
- >>> import scipy.stats
- >>> pearsonr = np.vectorize(scipy.stats.pearsonr,
- ... signature='(n),(n)->(),()')
- >>> pearsonr([[0, 1, 2, 3]], [[1, 2, 3, 4], [4, 3, 2, 1]])
- (array([ 1., -1.]), array([ 0., 0.]))
- Or for a vectorized convolution:
- >>> convolve = np.vectorize(np.convolve, signature='(n),(m)->(k)')
- >>> convolve(np.eye(4), [1, 2, 1])
- array([[1., 2., 1., 0., 0., 0.],
- [0., 1., 2., 1., 0., 0.],
- [0., 0., 1., 2., 1., 0.],
- [0., 0., 0., 1., 2., 1.]])
- """
- def __init__(self, pyfunc, otypes=None, doc=None, excluded=None,
- cache=False, signature=None):
- self.pyfunc = pyfunc
- self.cache = cache
- self.signature = signature
- self._ufunc = {} # Caching to improve default performance
- if doc is None:
- self.__doc__ = pyfunc.__doc__
- else:
- self.__doc__ = doc
- if isinstance(otypes, str):
- for char in otypes:
- if char not in typecodes['All']:
- raise ValueError("Invalid otype specified: %s" % (char,))
- elif iterable(otypes):
- otypes = ''.join([_nx.dtype(x).char for x in otypes])
- elif otypes is not None:
- raise ValueError("Invalid otype specification")
- self.otypes = otypes
- # Excluded variable support
- if excluded is None:
- excluded = set()
- self.excluded = set(excluded)
- if signature is not None:
- self._in_and_out_core_dims = _parse_gufunc_signature(signature)
- else:
- self._in_and_out_core_dims = None
- def __call__(self, *args, **kwargs):
- """
- Return arrays with the results of `pyfunc` broadcast (vectorized) over
- `args` and `kwargs` not in `excluded`.
- """
- excluded = self.excluded
- if not kwargs and not excluded:
- func = self.pyfunc
- vargs = args
- else:
- # The wrapper accepts only positional arguments: we use `names` and
- # `inds` to mutate `the_args` and `kwargs` to pass to the original
- # function.
- nargs = len(args)
- names = [_n for _n in kwargs if _n not in excluded]
- inds = [_i for _i in range(nargs) if _i not in excluded]
- the_args = list(args)
- def func(*vargs):
- for _n, _i in enumerate(inds):
- the_args[_i] = vargs[_n]
- kwargs.update(zip(names, vargs[len(inds):]))
- return self.pyfunc(*the_args, **kwargs)
- vargs = [args[_i] for _i in inds]
- vargs.extend([kwargs[_n] for _n in names])
- return self._vectorize_call(func=func, args=vargs)
- def _get_ufunc_and_otypes(self, func, args):
- """Return (ufunc, otypes)."""
- # frompyfunc will fail if args is empty
- if not args:
- raise ValueError('args can not be empty')
- if self.otypes is not None:
- otypes = self.otypes
- # self._ufunc is a dictionary whose keys are the number of
- # arguments (i.e. len(args)) and whose values are ufuncs created
- # by frompyfunc. len(args) can be different for different calls if
- # self.pyfunc has parameters with default values. We only use the
- # cache when func is self.pyfunc, which occurs when the call uses
- # only positional arguments and no arguments are excluded.
- nin = len(args)
- nout = len(self.otypes)
- if func is not self.pyfunc or nin not in self._ufunc:
- ufunc = frompyfunc(func, nin, nout)
- else:
- ufunc = None # We'll get it from self._ufunc
- if func is self.pyfunc:
- ufunc = self._ufunc.setdefault(nin, ufunc)
- else:
- # Get number of outputs and output types by calling the function on
- # the first entries of args. We also cache the result to prevent
- # the subsequent call when the ufunc is evaluated.
- # Assumes that ufunc first evaluates the 0th elements in the input
- # arrays (the input values are not checked to ensure this)
- args = [asarray(arg) for arg in args]
- if builtins.any(arg.size == 0 for arg in args):
- raise ValueError('cannot call `vectorize` on size 0 inputs '
- 'unless `otypes` is set')
- inputs = [arg.flat[0] for arg in args]
- outputs = func(*inputs)
- # Performance note: profiling indicates that -- for simple
- # functions at least -- this wrapping can almost double the
- # execution time.
- # Hence we make it optional.
- if self.cache:
- _cache = [outputs]
- def _func(*vargs):
- if _cache:
- return _cache.pop()
- else:
- return func(*vargs)
- else:
- _func = func
- if isinstance(outputs, tuple):
- nout = len(outputs)
- else:
- nout = 1
- outputs = (outputs,)
- otypes = ''.join([asarray(outputs[_k]).dtype.char
- for _k in range(nout)])
- # Performance note: profiling indicates that creating the ufunc is
- # not a significant cost compared with wrapping so it seems not
- # worth trying to cache this.
- ufunc = frompyfunc(_func, len(args), nout)
- return ufunc, otypes
- def _vectorize_call(self, func, args):
- """Vectorized call to `func` over positional `args`."""
- if self.signature is not None:
- res = self._vectorize_call_with_signature(func, args)
- elif not args:
- res = func()
- else:
- ufunc, otypes = self._get_ufunc_and_otypes(func=func, args=args)
- # Convert args to object arrays first
- inputs = [array(a, copy=False, subok=True, dtype=object)
- for a in args]
- outputs = ufunc(*inputs)
- if ufunc.nout == 1:
- res = array(outputs, copy=False, subok=True, dtype=otypes[0])
- else:
- res = tuple([array(x, copy=False, subok=True, dtype=t)
- for x, t in zip(outputs, otypes)])
- return res
- def _vectorize_call_with_signature(self, func, args):
- """Vectorized call over positional arguments with a signature."""
- input_core_dims, output_core_dims = self._in_and_out_core_dims
- if len(args) != len(input_core_dims):
- raise TypeError('wrong number of positional arguments: '
- 'expected %r, got %r'
- % (len(input_core_dims), len(args)))
- args = tuple(asanyarray(arg) for arg in args)
- broadcast_shape, dim_sizes = _parse_input_dimensions(
- args, input_core_dims)
- input_shapes = _calculate_shapes(broadcast_shape, dim_sizes,
- input_core_dims)
- args = [np.broadcast_to(arg, shape, subok=True)
- for arg, shape in zip(args, input_shapes)]
- outputs = None
- otypes = self.otypes
- nout = len(output_core_dims)
- for index in np.ndindex(*broadcast_shape):
- results = func(*(arg[index] for arg in args))
- n_results = len(results) if isinstance(results, tuple) else 1
- if nout != n_results:
- raise ValueError(
- 'wrong number of outputs from pyfunc: expected %r, got %r'
- % (nout, n_results))
- if nout == 1:
- results = (results,)
- if outputs is None:
- for result, core_dims in zip(results, output_core_dims):
- _update_dim_sizes(dim_sizes, result, core_dims)
- if otypes is None:
- otypes = [asarray(result).dtype for result in results]
- outputs = _create_arrays(broadcast_shape, dim_sizes,
- output_core_dims, otypes)
- for output, result in zip(outputs, results):
- output[index] = result
- if outputs is None:
- # did not call the function even once
- if otypes is None:
- raise ValueError('cannot call `vectorize` on size 0 inputs '
- 'unless `otypes` is set')
- if builtins.any(dim not in dim_sizes
- for dims in output_core_dims
- for dim in dims):
- raise ValueError('cannot call `vectorize` with a signature '
- 'including new output dimensions on size 0 '
- 'inputs')
- outputs = _create_arrays(broadcast_shape, dim_sizes,
- output_core_dims, otypes)
- return outputs[0] if nout == 1 else outputs
- def _cov_dispatcher(m, y=None, rowvar=None, bias=None, ddof=None,
- fweights=None, aweights=None, *, dtype=None):
- return (m, y, fweights, aweights)
- @array_function_dispatch(_cov_dispatcher)
- def cov(m, y=None, rowvar=True, bias=False, ddof=None, fweights=None,
- aweights=None, *, dtype=None):
- """
- Estimate a covariance matrix, given data and weights.
- Covariance indicates the level to which two variables vary together.
- If we examine N-dimensional samples, :math:`X = [x_1, x_2, ... x_N]^T`,
- then the covariance matrix element :math:`C_{ij}` is the covariance of
- :math:`x_i` and :math:`x_j`. The element :math:`C_{ii}` is the variance
- of :math:`x_i`.
- See the notes for an outline of the algorithm.
- Parameters
- ----------
- m : array_like
- A 1-D or 2-D array containing multiple variables and observations.
- Each row of `m` represents a variable, and each column a single
- observation of all those variables. Also see `rowvar` below.
- y : array_like, optional
- An additional set of variables and observations. `y` has the same form
- as that of `m`.
- rowvar : bool, optional
- If `rowvar` is True (default), then each row represents a
- variable, with observations in the columns. Otherwise, the relationship
- is transposed: each column represents a variable, while the rows
- contain observations.
- bias : bool, optional
- Default normalization (False) is by ``(N - 1)``, where ``N`` is the
- number of observations given (unbiased estimate). If `bias` is True,
- then normalization is by ``N``. These values can be overridden by using
- the keyword ``ddof`` in numpy versions >= 1.5.
- ddof : int, optional
- If not ``None`` the default value implied by `bias` is overridden.
- Note that ``ddof=1`` will return the unbiased estimate, even if both
- `fweights` and `aweights` are specified, and ``ddof=0`` will return
- the simple average. See the notes for the details. The default value
- is ``None``.
- .. versionadded:: 1.5
- fweights : array_like, int, optional
- 1-D array of integer frequency weights; the number of times each
- observation vector should be repeated.
- .. versionadded:: 1.10
- aweights : array_like, optional
- 1-D array of observation vector weights. These relative weights are
- typically large for observations considered "important" and smaller for
- observations considered less "important". If ``ddof=0`` the array of
- weights can be used to assign probabilities to observation vectors.
- .. versionadded:: 1.10
- dtype : data-type, optional
- Data-type of the result. By default, the return data-type will have
- at least `numpy.float64` precision.
- .. versionadded:: 1.20
- Returns
- -------
- out : ndarray
- The covariance matrix of the variables.
- See Also
- --------
- corrcoef : Normalized covariance matrix
- Notes
- -----
- Assume that the observations are in the columns of the observation
- array `m` and let ``f = fweights`` and ``a = aweights`` for brevity. The
- steps to compute the weighted covariance are as follows::
- >>> m = np.arange(10, dtype=np.float64)
- >>> f = np.arange(10) * 2
- >>> a = np.arange(10) ** 2.
- >>> ddof = 1
- >>> w = f * a
- >>> v1 = np.sum(w)
- >>> v2 = np.sum(w * a)
- >>> m -= np.sum(m * w, axis=None, keepdims=True) / v1
- >>> cov = np.dot(m * w, m.T) * v1 / (v1**2 - ddof * v2)
- Note that when ``a == 1``, the normalization factor
- ``v1 / (v1**2 - ddof * v2)`` goes over to ``1 / (np.sum(f) - ddof)``
- as it should.
- Examples
- --------
- Consider two variables, :math:`x_0` and :math:`x_1`, which
- correlate perfectly, but in opposite directions:
- >>> x = np.array([[0, 2], [1, 1], [2, 0]]).T
- >>> x
- array([[0, 1, 2],
- [2, 1, 0]])
- Note how :math:`x_0` increases while :math:`x_1` decreases. The covariance
- matrix shows this clearly:
- >>> np.cov(x)
- array([[ 1., -1.],
- [-1., 1.]])
- Note that element :math:`C_{0,1}`, which shows the correlation between
- :math:`x_0` and :math:`x_1`, is negative.
- Further, note how `x` and `y` are combined:
- >>> x = [-2.1, -1, 4.3]
- >>> y = [3, 1.1, 0.12]
- >>> X = np.stack((x, y), axis=0)
- >>> np.cov(X)
- array([[11.71 , -4.286 ], # may vary
- [-4.286 , 2.144133]])
- >>> np.cov(x, y)
- array([[11.71 , -4.286 ], # may vary
- [-4.286 , 2.144133]])
- >>> np.cov(x)
- array(11.71)
- """
- # Check inputs
- if ddof is not None and ddof != int(ddof):
- raise ValueError(
- "ddof must be integer")
- # Handles complex arrays too
- m = np.asarray(m)
- if m.ndim > 2:
- raise ValueError("m has more than 2 dimensions")
- if y is not None:
- y = np.asarray(y)
- if y.ndim > 2:
- raise ValueError("y has more than 2 dimensions")
- if dtype is None:
- if y is None:
- dtype = np.result_type(m, np.float64)
- else:
- dtype = np.result_type(m, y, np.float64)
- X = array(m, ndmin=2, dtype=dtype)
- if not rowvar and X.shape[0] != 1:
- X = X.T
- if X.shape[0] == 0:
- return np.array([]).reshape(0, 0)
- if y is not None:
- y = array(y, copy=False, ndmin=2, dtype=dtype)
- if not rowvar and y.shape[0] != 1:
- y = y.T
- X = np.concatenate((X, y), axis=0)
- if ddof is None:
- if bias == 0:
- ddof = 1
- else:
- ddof = 0
- # Get the product of frequencies and weights
- w = None
- if fweights is not None:
- fweights = np.asarray(fweights, dtype=float)
- if not np.all(fweights == np.around(fweights)):
- raise TypeError(
- "fweights must be integer")
- if fweights.ndim > 1:
- raise RuntimeError(
- "cannot handle multidimensional fweights")
- if fweights.shape[0] != X.shape[1]:
- raise RuntimeError(
- "incompatible numbers of samples and fweights")
- if any(fweights < 0):
- raise ValueError(
- "fweights cannot be negative")
- w = fweights
- if aweights is not None:
- aweights = np.asarray(aweights, dtype=float)
- if aweights.ndim > 1:
- raise RuntimeError(
- "cannot handle multidimensional aweights")
- if aweights.shape[0] != X.shape[1]:
- raise RuntimeError(
- "incompatible numbers of samples and aweights")
- if any(aweights < 0):
- raise ValueError(
- "aweights cannot be negative")
- if w is None:
- w = aweights
- else:
- w *= aweights
- avg, w_sum = average(X, axis=1, weights=w, returned=True)
- w_sum = w_sum[0]
- # Determine the normalization
- if w is None:
- fact = X.shape[1] - ddof
- elif ddof == 0:
- fact = w_sum
- elif aweights is None:
- fact = w_sum - ddof
- else:
- fact = w_sum - ddof*sum(w*aweights)/w_sum
- if fact <= 0:
- warnings.warn("Degrees of freedom <= 0 for slice",
- RuntimeWarning, stacklevel=3)
- fact = 0.0
- X -= avg[:, None]
- if w is None:
- X_T = X.T
- else:
- X_T = (X*w).T
- c = dot(X, X_T.conj())
- c *= np.true_divide(1, fact)
- return c.squeeze()
- def _corrcoef_dispatcher(x, y=None, rowvar=None, bias=None, ddof=None, *,
- dtype=None):
- return (x, y)
- @array_function_dispatch(_corrcoef_dispatcher)
- def corrcoef(x, y=None, rowvar=True, bias=np._NoValue, ddof=np._NoValue, *,
- dtype=None):
- """
- Return Pearson product-moment correlation coefficients.
- Please refer to the documentation for `cov` for more detail. The
- relationship between the correlation coefficient matrix, `R`, and the
- covariance matrix, `C`, is
- .. math:: R_{ij} = \\frac{ C_{ij} } { \\sqrt{ C_{ii} * C_{jj} } }
- The values of `R` are between -1 and 1, inclusive.
- Parameters
- ----------
- x : array_like
- A 1-D or 2-D array containing multiple variables and observations.
- Each row of `x` represents a variable, and each column a single
- observation of all those variables. Also see `rowvar` below.
- y : array_like, optional
- An additional set of variables and observations. `y` has the same
- shape as `x`.
- rowvar : bool, optional
- If `rowvar` is True (default), then each row represents a
- variable, with observations in the columns. Otherwise, the relationship
- is transposed: each column represents a variable, while the rows
- contain observations.
- bias : _NoValue, optional
- Has no effect, do not use.
- .. deprecated:: 1.10.0
- ddof : _NoValue, optional
- Has no effect, do not use.
- .. deprecated:: 1.10.0
- dtype : data-type, optional
- Data-type of the result. By default, the return data-type will have
- at least `numpy.float64` precision.
- .. versionadded:: 1.20
- Returns
- -------
- R : ndarray
- The correlation coefficient matrix of the variables.
- See Also
- --------
- cov : Covariance matrix
- Notes
- -----
- Due to floating point rounding the resulting array may not be Hermitian,
- the diagonal elements may not be 1, and the elements may not satisfy the
- inequality abs(a) <= 1. The real and imaginary parts are clipped to the
- interval [-1, 1] in an attempt to improve on that situation but is not
- much help in the complex case.
- This function accepts but discards arguments `bias` and `ddof`. This is
- for backwards compatibility with previous versions of this function. These
- arguments had no effect on the return values of the function and can be
- safely ignored in this and previous versions of numpy.
- Examples
- --------
- In this example we generate two random arrays, ``xarr`` and ``yarr``, and
- compute the row-wise and column-wise Pearson correlation coefficients,
- ``R``. Since ``rowvar`` is true by default, we first find the row-wise
- Pearson correlation coefficients between the variables of ``xarr``.
- >>> import numpy as np
- >>> rng = np.random.default_rng(seed=42)
- >>> xarr = rng.random((3, 3))
- >>> xarr
- array([[0.77395605, 0.43887844, 0.85859792],
- [0.69736803, 0.09417735, 0.97562235],
- [0.7611397 , 0.78606431, 0.12811363]])
- >>> R1 = np.corrcoef(xarr)
- >>> R1
- array([[ 1. , 0.99256089, -0.68080986],
- [ 0.99256089, 1. , -0.76492172],
- [-0.68080986, -0.76492172, 1. ]])
- If we add another set of variables and observations ``yarr``, we can
- compute the row-wise Pearson correlation coefficients between the
- variables in ``xarr`` and ``yarr``.
- >>> yarr = rng.random((3, 3))
- >>> yarr
- array([[0.45038594, 0.37079802, 0.92676499],
- [0.64386512, 0.82276161, 0.4434142 ],
- [0.22723872, 0.55458479, 0.06381726]])
- >>> R2 = np.corrcoef(xarr, yarr)
- >>> R2
- array([[ 1. , 0.99256089, -0.68080986, 0.75008178, -0.934284 ,
- -0.99004057],
- [ 0.99256089, 1. , -0.76492172, 0.82502011, -0.97074098,
- -0.99981569],
- [-0.68080986, -0.76492172, 1. , -0.99507202, 0.89721355,
- 0.77714685],
- [ 0.75008178, 0.82502011, -0.99507202, 1. , -0.93657855,
- -0.83571711],
- [-0.934284 , -0.97074098, 0.89721355, -0.93657855, 1. ,
- 0.97517215],
- [-0.99004057, -0.99981569, 0.77714685, -0.83571711, 0.97517215,
- 1. ]])
- Finally if we use the option ``rowvar=False``, the columns are now
- being treated as the variables and we will find the column-wise Pearson
- correlation coefficients between variables in ``xarr`` and ``yarr``.
- >>> R3 = np.corrcoef(xarr, yarr, rowvar=False)
- >>> R3
- array([[ 1. , 0.77598074, -0.47458546, -0.75078643, -0.9665554 ,
- 0.22423734],
- [ 0.77598074, 1. , -0.92346708, -0.99923895, -0.58826587,
- -0.44069024],
- [-0.47458546, -0.92346708, 1. , 0.93773029, 0.23297648,
- 0.75137473],
- [-0.75078643, -0.99923895, 0.93773029, 1. , 0.55627469,
- 0.47536961],
- [-0.9665554 , -0.58826587, 0.23297648, 0.55627469, 1. ,
- -0.46666491],
- [ 0.22423734, -0.44069024, 0.75137473, 0.47536961, -0.46666491,
- 1. ]])
- """
- if bias is not np._NoValue or ddof is not np._NoValue:
- # 2015-03-15, 1.10
- warnings.warn('bias and ddof have no effect and are deprecated',
- DeprecationWarning, stacklevel=3)
- c = cov(x, y, rowvar, dtype=dtype)
- try:
- d = diag(c)
- except ValueError:
- # scalar covariance
- # nan if incorrect value (nan, inf, 0), 1 otherwise
- return c / c
- stddev = sqrt(d.real)
- c /= stddev[:, None]
- c /= stddev[None, :]
- # Clip real and imaginary parts to [-1, 1]. This does not guarantee
- # abs(a[i,j]) <= 1 for complex arrays, but is the best we can do without
- # excessive work.
- np.clip(c.real, -1, 1, out=c.real)
- if np.iscomplexobj(c):
- np.clip(c.imag, -1, 1, out=c.imag)
- return c
- @set_module('numpy')
- def blackman(M):
- """
- Return the Blackman window.
- The Blackman window is a taper formed by using the first three
- terms of a summation of cosines. It was designed to have close to the
- minimal leakage possible. It is close to optimal, only slightly worse
- than a Kaiser window.
- Parameters
- ----------
- M : int
- Number of points in the output window. If zero or less, an empty
- array is returned.
- Returns
- -------
- out : ndarray
- The window, with the maximum value normalized to one (the value one
- appears only if the number of samples is odd).
- See Also
- --------
- bartlett, hamming, hanning, kaiser
- Notes
- -----
- The Blackman window is defined as
- .. math:: w(n) = 0.42 - 0.5 \\cos(2\\pi n/M) + 0.08 \\cos(4\\pi n/M)
- Most references to the Blackman window come from the signal processing
- literature, where it is used as one of many windowing functions for
- smoothing values. It is also known as an apodization (which means
- "removing the foot", i.e. smoothing discontinuities at the beginning
- and end of the sampled signal) or tapering function. It is known as a
- "near optimal" tapering function, almost as good (by some measures)
- as the kaiser window.
- References
- ----------
- Blackman, R.B. and Tukey, J.W., (1958) The measurement of power spectra,
- Dover Publications, New York.
- Oppenheim, A.V., and R.W. Schafer. Discrete-Time Signal Processing.
- Upper Saddle River, NJ: Prentice-Hall, 1999, pp. 468-471.
- Examples
- --------
- >>> import matplotlib.pyplot as plt
- >>> np.blackman(12)
- array([-1.38777878e-17, 3.26064346e-02, 1.59903635e-01, # may vary
- 4.14397981e-01, 7.36045180e-01, 9.67046769e-01,
- 9.67046769e-01, 7.36045180e-01, 4.14397981e-01,
- 1.59903635e-01, 3.26064346e-02, -1.38777878e-17])
- Plot the window and the frequency response:
- >>> from numpy.fft import fft, fftshift
- >>> window = np.blackman(51)
- >>> plt.plot(window)
- [<matplotlib.lines.Line2D object at 0x...>]
- >>> plt.title("Blackman window")
- Text(0.5, 1.0, 'Blackman window')
- >>> plt.ylabel("Amplitude")
- Text(0, 0.5, 'Amplitude')
- >>> plt.xlabel("Sample")
- Text(0.5, 0, 'Sample')
- >>> plt.show()
- >>> plt.figure()
- <Figure size 640x480 with 0 Axes>
- >>> A = fft(window, 2048) / 25.5
- >>> mag = np.abs(fftshift(A))
- >>> freq = np.linspace(-0.5, 0.5, len(A))
- >>> with np.errstate(divide='ignore', invalid='ignore'):
- ... response = 20 * np.log10(mag)
- ...
- >>> response = np.clip(response, -100, 100)
- >>> plt.plot(freq, response)
- [<matplotlib.lines.Line2D object at 0x...>]
- >>> plt.title("Frequency response of Blackman window")
- Text(0.5, 1.0, 'Frequency response of Blackman window')
- >>> plt.ylabel("Magnitude [dB]")
- Text(0, 0.5, 'Magnitude [dB]')
- >>> plt.xlabel("Normalized frequency [cycles per sample]")
- Text(0.5, 0, 'Normalized frequency [cycles per sample]')
- >>> _ = plt.axis('tight')
- >>> plt.show()
- """
- if M < 1:
- return array([])
- if M == 1:
- return ones(1, float)
- n = arange(1-M, M, 2)
- return 0.42 + 0.5*cos(pi*n/(M-1)) + 0.08*cos(2.0*pi*n/(M-1))
- @set_module('numpy')
- def bartlett(M):
- """
- Return the Bartlett window.
- The Bartlett window is very similar to a triangular window, except
- that the end points are at zero. It is often used in signal
- processing for tapering a signal, without generating too much
- ripple in the frequency domain.
- Parameters
- ----------
- M : int
- Number of points in the output window. If zero or less, an
- empty array is returned.
- Returns
- -------
- out : array
- The triangular window, with the maximum value normalized to one
- (the value one appears only if the number of samples is odd), with
- the first and last samples equal to zero.
- See Also
- --------
- blackman, hamming, hanning, kaiser
- Notes
- -----
- The Bartlett window is defined as
- .. math:: w(n) = \\frac{2}{M-1} \\left(
- \\frac{M-1}{2} - \\left|n - \\frac{M-1}{2}\\right|
- \\right)
- Most references to the Bartlett window come from the signal
- processing literature, where it is used as one of many windowing
- functions for smoothing values. Note that convolution with this
- window produces linear interpolation. It is also known as an
- apodization (which means"removing the foot", i.e. smoothing
- discontinuities at the beginning and end of the sampled signal) or
- tapering function. The fourier transform of the Bartlett is the product
- of two sinc functions.
- Note the excellent discussion in Kanasewich.
- References
- ----------
- .. [1] M.S. Bartlett, "Periodogram Analysis and Continuous Spectra",
- Biometrika 37, 1-16, 1950.
- .. [2] E.R. Kanasewich, "Time Sequence Analysis in Geophysics",
- The University of Alberta Press, 1975, pp. 109-110.
- .. [3] A.V. Oppenheim and R.W. Schafer, "Discrete-Time Signal
- Processing", Prentice-Hall, 1999, pp. 468-471.
- .. [4] Wikipedia, "Window function",
- https://en.wikipedia.org/wiki/Window_function
- .. [5] W.H. Press, B.P. Flannery, S.A. Teukolsky, and W.T. Vetterling,
- "Numerical Recipes", Cambridge University Press, 1986, page 429.
- Examples
- --------
- >>> import matplotlib.pyplot as plt
- >>> np.bartlett(12)
- array([ 0. , 0.18181818, 0.36363636, 0.54545455, 0.72727273, # may vary
- 0.90909091, 0.90909091, 0.72727273, 0.54545455, 0.36363636,
- 0.18181818, 0. ])
- Plot the window and its frequency response (requires SciPy and matplotlib):
- >>> from numpy.fft import fft, fftshift
- >>> window = np.bartlett(51)
- >>> plt.plot(window)
- [<matplotlib.lines.Line2D object at 0x...>]
- >>> plt.title("Bartlett window")
- Text(0.5, 1.0, 'Bartlett window')
- >>> plt.ylabel("Amplitude")
- Text(0, 0.5, 'Amplitude')
- >>> plt.xlabel("Sample")
- Text(0.5, 0, 'Sample')
- >>> plt.show()
- >>> plt.figure()
- <Figure size 640x480 with 0 Axes>
- >>> A = fft(window, 2048) / 25.5
- >>> mag = np.abs(fftshift(A))
- >>> freq = np.linspace(-0.5, 0.5, len(A))
- >>> with np.errstate(divide='ignore', invalid='ignore'):
- ... response = 20 * np.log10(mag)
- ...
- >>> response = np.clip(response, -100, 100)
- >>> plt.plot(freq, response)
- [<matplotlib.lines.Line2D object at 0x...>]
- >>> plt.title("Frequency response of Bartlett window")
- Text(0.5, 1.0, 'Frequency response of Bartlett window')
- >>> plt.ylabel("Magnitude [dB]")
- Text(0, 0.5, 'Magnitude [dB]')
- >>> plt.xlabel("Normalized frequency [cycles per sample]")
- Text(0.5, 0, 'Normalized frequency [cycles per sample]')
- >>> _ = plt.axis('tight')
- >>> plt.show()
- """
- if M < 1:
- return array([])
- if M == 1:
- return ones(1, float)
- n = arange(1-M, M, 2)
- return where(less_equal(n, 0), 1 + n/(M-1), 1 - n/(M-1))
- @set_module('numpy')
- def hanning(M):
- """
- Return the Hanning window.
- The Hanning window is a taper formed by using a weighted cosine.
- Parameters
- ----------
- M : int
- Number of points in the output window. If zero or less, an
- empty array is returned.
- Returns
- -------
- out : ndarray, shape(M,)
- The window, with the maximum value normalized to one (the value
- one appears only if `M` is odd).
- See Also
- --------
- bartlett, blackman, hamming, kaiser
- Notes
- -----
- The Hanning window is defined as
- .. math:: w(n) = 0.5 - 0.5cos\\left(\\frac{2\\pi{n}}{M-1}\\right)
- \\qquad 0 \\leq n \\leq M-1
- The Hanning was named for Julius von Hann, an Austrian meteorologist.
- It is also known as the Cosine Bell. Some authors prefer that it be
- called a Hann window, to help avoid confusion with the very similar
- Hamming window.
- Most references to the Hanning window come from the signal processing
- literature, where it is used as one of many windowing functions for
- smoothing values. It is also known as an apodization (which means
- "removing the foot", i.e. smoothing discontinuities at the beginning
- and end of the sampled signal) or tapering function.
- References
- ----------
- .. [1] Blackman, R.B. and Tukey, J.W., (1958) The measurement of power
- spectra, Dover Publications, New York.
- .. [2] E.R. Kanasewich, "Time Sequence Analysis in Geophysics",
- The University of Alberta Press, 1975, pp. 106-108.
- .. [3] Wikipedia, "Window function",
- https://en.wikipedia.org/wiki/Window_function
- .. [4] W.H. Press, B.P. Flannery, S.A. Teukolsky, and W.T. Vetterling,
- "Numerical Recipes", Cambridge University Press, 1986, page 425.
- Examples
- --------
- >>> np.hanning(12)
- array([0. , 0.07937323, 0.29229249, 0.57115742, 0.82743037,
- 0.97974649, 0.97974649, 0.82743037, 0.57115742, 0.29229249,
- 0.07937323, 0. ])
- Plot the window and its frequency response:
- >>> import matplotlib.pyplot as plt
- >>> from numpy.fft import fft, fftshift
- >>> window = np.hanning(51)
- >>> plt.plot(window)
- [<matplotlib.lines.Line2D object at 0x...>]
- >>> plt.title("Hann window")
- Text(0.5, 1.0, 'Hann window')
- >>> plt.ylabel("Amplitude")
- Text(0, 0.5, 'Amplitude')
- >>> plt.xlabel("Sample")
- Text(0.5, 0, 'Sample')
- >>> plt.show()
- >>> plt.figure()
- <Figure size 640x480 with 0 Axes>
- >>> A = fft(window, 2048) / 25.5
- >>> mag = np.abs(fftshift(A))
- >>> freq = np.linspace(-0.5, 0.5, len(A))
- >>> with np.errstate(divide='ignore', invalid='ignore'):
- ... response = 20 * np.log10(mag)
- ...
- >>> response = np.clip(response, -100, 100)
- >>> plt.plot(freq, response)
- [<matplotlib.lines.Line2D object at 0x...>]
- >>> plt.title("Frequency response of the Hann window")
- Text(0.5, 1.0, 'Frequency response of the Hann window')
- >>> plt.ylabel("Magnitude [dB]")
- Text(0, 0.5, 'Magnitude [dB]')
- >>> plt.xlabel("Normalized frequency [cycles per sample]")
- Text(0.5, 0, 'Normalized frequency [cycles per sample]')
- >>> plt.axis('tight')
- ...
- >>> plt.show()
- """
- if M < 1:
- return array([])
- if M == 1:
- return ones(1, float)
- n = arange(1-M, M, 2)
- return 0.5 + 0.5*cos(pi*n/(M-1))
- @set_module('numpy')
- def hamming(M):
- """
- Return the Hamming window.
- The Hamming window is a taper formed by using a weighted cosine.
- Parameters
- ----------
- M : int
- Number of points in the output window. If zero or less, an
- empty array is returned.
- Returns
- -------
- out : ndarray
- The window, with the maximum value normalized to one (the value
- one appears only if the number of samples is odd).
- See Also
- --------
- bartlett, blackman, hanning, kaiser
- Notes
- -----
- The Hamming window is defined as
- .. math:: w(n) = 0.54 - 0.46cos\\left(\\frac{2\\pi{n}}{M-1}\\right)
- \\qquad 0 \\leq n \\leq M-1
- The Hamming was named for R. W. Hamming, an associate of J. W. Tukey
- and is described in Blackman and Tukey. It was recommended for
- smoothing the truncated autocovariance function in the time domain.
- Most references to the Hamming window come from the signal processing
- literature, where it is used as one of many windowing functions for
- smoothing values. It is also known as an apodization (which means
- "removing the foot", i.e. smoothing discontinuities at the beginning
- and end of the sampled signal) or tapering function.
- References
- ----------
- .. [1] Blackman, R.B. and Tukey, J.W., (1958) The measurement of power
- spectra, Dover Publications, New York.
- .. [2] E.R. Kanasewich, "Time Sequence Analysis in Geophysics", The
- University of Alberta Press, 1975, pp. 109-110.
- .. [3] Wikipedia, "Window function",
- https://en.wikipedia.org/wiki/Window_function
- .. [4] W.H. Press, B.P. Flannery, S.A. Teukolsky, and W.T. Vetterling,
- "Numerical Recipes", Cambridge University Press, 1986, page 425.
- Examples
- --------
- >>> np.hamming(12)
- array([ 0.08 , 0.15302337, 0.34890909, 0.60546483, 0.84123594, # may vary
- 0.98136677, 0.98136677, 0.84123594, 0.60546483, 0.34890909,
- 0.15302337, 0.08 ])
- Plot the window and the frequency response:
- >>> import matplotlib.pyplot as plt
- >>> from numpy.fft import fft, fftshift
- >>> window = np.hamming(51)
- >>> plt.plot(window)
- [<matplotlib.lines.Line2D object at 0x...>]
- >>> plt.title("Hamming window")
- Text(0.5, 1.0, 'Hamming window')
- >>> plt.ylabel("Amplitude")
- Text(0, 0.5, 'Amplitude')
- >>> plt.xlabel("Sample")
- Text(0.5, 0, 'Sample')
- >>> plt.show()
- >>> plt.figure()
- <Figure size 640x480 with 0 Axes>
- >>> A = fft(window, 2048) / 25.5
- >>> mag = np.abs(fftshift(A))
- >>> freq = np.linspace(-0.5, 0.5, len(A))
- >>> response = 20 * np.log10(mag)
- >>> response = np.clip(response, -100, 100)
- >>> plt.plot(freq, response)
- [<matplotlib.lines.Line2D object at 0x...>]
- >>> plt.title("Frequency response of Hamming window")
- Text(0.5, 1.0, 'Frequency response of Hamming window')
- >>> plt.ylabel("Magnitude [dB]")
- Text(0, 0.5, 'Magnitude [dB]')
- >>> plt.xlabel("Normalized frequency [cycles per sample]")
- Text(0.5, 0, 'Normalized frequency [cycles per sample]')
- >>> plt.axis('tight')
- ...
- >>> plt.show()
- """
- if M < 1:
- return array([])
- if M == 1:
- return ones(1, float)
- n = arange(1-M, M, 2)
- return 0.54 + 0.46*cos(pi*n/(M-1))
- ## Code from cephes for i0
- _i0A = [
- -4.41534164647933937950E-18,
- 3.33079451882223809783E-17,
- -2.43127984654795469359E-16,
- 1.71539128555513303061E-15,
- -1.16853328779934516808E-14,
- 7.67618549860493561688E-14,
- -4.85644678311192946090E-13,
- 2.95505266312963983461E-12,
- -1.72682629144155570723E-11,
- 9.67580903537323691224E-11,
- -5.18979560163526290666E-10,
- 2.65982372468238665035E-9,
- -1.30002500998624804212E-8,
- 6.04699502254191894932E-8,
- -2.67079385394061173391E-7,
- 1.11738753912010371815E-6,
- -4.41673835845875056359E-6,
- 1.64484480707288970893E-5,
- -5.75419501008210370398E-5,
- 1.88502885095841655729E-4,
- -5.76375574538582365885E-4,
- 1.63947561694133579842E-3,
- -4.32430999505057594430E-3,
- 1.05464603945949983183E-2,
- -2.37374148058994688156E-2,
- 4.93052842396707084878E-2,
- -9.49010970480476444210E-2,
- 1.71620901522208775349E-1,
- -3.04682672343198398683E-1,
- 6.76795274409476084995E-1
- ]
- _i0B = [
- -7.23318048787475395456E-18,
- -4.83050448594418207126E-18,
- 4.46562142029675999901E-17,
- 3.46122286769746109310E-17,
- -2.82762398051658348494E-16,
- -3.42548561967721913462E-16,
- 1.77256013305652638360E-15,
- 3.81168066935262242075E-15,
- -9.55484669882830764870E-15,
- -4.15056934728722208663E-14,
- 1.54008621752140982691E-14,
- 3.85277838274214270114E-13,
- 7.18012445138366623367E-13,
- -1.79417853150680611778E-12,
- -1.32158118404477131188E-11,
- -3.14991652796324136454E-11,
- 1.18891471078464383424E-11,
- 4.94060238822496958910E-10,
- 3.39623202570838634515E-9,
- 2.26666899049817806459E-8,
- 2.04891858946906374183E-7,
- 2.89137052083475648297E-6,
- 6.88975834691682398426E-5,
- 3.36911647825569408990E-3,
- 8.04490411014108831608E-1
- ]
- def _chbevl(x, vals):
- b0 = vals[0]
- b1 = 0.0
- for i in range(1, len(vals)):
- b2 = b1
- b1 = b0
- b0 = x*b1 - b2 + vals[i]
- return 0.5*(b0 - b2)
- def _i0_1(x):
- return exp(x) * _chbevl(x/2.0-2, _i0A)
- def _i0_2(x):
- return exp(x) * _chbevl(32.0/x - 2.0, _i0B) / sqrt(x)
- def _i0_dispatcher(x):
- return (x,)
- @array_function_dispatch(_i0_dispatcher)
- def i0(x):
- """
- Modified Bessel function of the first kind, order 0.
- Usually denoted :math:`I_0`.
- Parameters
- ----------
- x : array_like of float
- Argument of the Bessel function.
- Returns
- -------
- out : ndarray, shape = x.shape, dtype = float
- The modified Bessel function evaluated at each of the elements of `x`.
- See Also
- --------
- scipy.special.i0, scipy.special.iv, scipy.special.ive
- Notes
- -----
- The scipy implementation is recommended over this function: it is a
- proper ufunc written in C, and more than an order of magnitude faster.
- We use the algorithm published by Clenshaw [1]_ and referenced by
- Abramowitz and Stegun [2]_, for which the function domain is
- partitioned into the two intervals [0,8] and (8,inf), and Chebyshev
- polynomial expansions are employed in each interval. Relative error on
- the domain [0,30] using IEEE arithmetic is documented [3]_ as having a
- peak of 5.8e-16 with an rms of 1.4e-16 (n = 30000).
- References
- ----------
- .. [1] C. W. Clenshaw, "Chebyshev series for mathematical functions", in
- *National Physical Laboratory Mathematical Tables*, vol. 5, London:
- Her Majesty's Stationery Office, 1962.
- .. [2] M. Abramowitz and I. A. Stegun, *Handbook of Mathematical
- Functions*, 10th printing, New York: Dover, 1964, pp. 379.
- http://www.math.sfu.ca/~cbm/aands/page_379.htm
- .. [3] https://metacpan.org/pod/distribution/Math-Cephes/lib/Math/Cephes.pod#i0:-Modified-Bessel-function-of-order-zero
- Examples
- --------
- >>> np.i0(0.)
- array(1.0)
- >>> np.i0([0, 1, 2, 3])
- array([1. , 1.26606588, 2.2795853 , 4.88079259])
- """
- x = np.asanyarray(x)
- if x.dtype.kind == 'c':
- raise TypeError("i0 not supported for complex values")
- if x.dtype.kind != 'f':
- x = x.astype(float)
- x = np.abs(x)
- return piecewise(x, [x <= 8.0], [_i0_1, _i0_2])
- ## End of cephes code for i0
- @set_module('numpy')
- def kaiser(M, beta):
- """
- Return the Kaiser window.
- The Kaiser window is a taper formed by using a Bessel function.
- Parameters
- ----------
- M : int
- Number of points in the output window. If zero or less, an
- empty array is returned.
- beta : float
- Shape parameter for window.
- Returns
- -------
- out : array
- The window, with the maximum value normalized to one (the value
- one appears only if the number of samples is odd).
- See Also
- --------
- bartlett, blackman, hamming, hanning
- Notes
- -----
- The Kaiser window is defined as
- .. math:: w(n) = I_0\\left( \\beta \\sqrt{1-\\frac{4n^2}{(M-1)^2}}
- \\right)/I_0(\\beta)
- with
- .. math:: \\quad -\\frac{M-1}{2} \\leq n \\leq \\frac{M-1}{2},
- where :math:`I_0` is the modified zeroth-order Bessel function.
- The Kaiser was named for Jim Kaiser, who discovered a simple
- approximation to the DPSS window based on Bessel functions. The Kaiser
- window is a very good approximation to the Digital Prolate Spheroidal
- Sequence, or Slepian window, which is the transform which maximizes the
- energy in the main lobe of the window relative to total energy.
- The Kaiser can approximate many other windows by varying the beta
- parameter.
- ==== =======================
- beta Window shape
- ==== =======================
- 0 Rectangular
- 5 Similar to a Hamming
- 6 Similar to a Hanning
- 8.6 Similar to a Blackman
- ==== =======================
- A beta value of 14 is probably a good starting point. Note that as beta
- gets large, the window narrows, and so the number of samples needs to be
- large enough to sample the increasingly narrow spike, otherwise NaNs will
- get returned.
- Most references to the Kaiser window come from the signal processing
- literature, where it is used as one of many windowing functions for
- smoothing values. It is also known as an apodization (which means
- "removing the foot", i.e. smoothing discontinuities at the beginning
- and end of the sampled signal) or tapering function.
- References
- ----------
- .. [1] J. F. Kaiser, "Digital Filters" - Ch 7 in "Systems analysis by
- digital computer", Editors: F.F. Kuo and J.F. Kaiser, p 218-285.
- John Wiley and Sons, New York, (1966).
- .. [2] E.R. Kanasewich, "Time Sequence Analysis in Geophysics", The
- University of Alberta Press, 1975, pp. 177-178.
- .. [3] Wikipedia, "Window function",
- https://en.wikipedia.org/wiki/Window_function
- Examples
- --------
- >>> import matplotlib.pyplot as plt
- >>> np.kaiser(12, 14)
- array([7.72686684e-06, 3.46009194e-03, 4.65200189e-02, # may vary
- 2.29737120e-01, 5.99885316e-01, 9.45674898e-01,
- 9.45674898e-01, 5.99885316e-01, 2.29737120e-01,
- 4.65200189e-02, 3.46009194e-03, 7.72686684e-06])
- Plot the window and the frequency response:
- >>> from numpy.fft import fft, fftshift
- >>> window = np.kaiser(51, 14)
- >>> plt.plot(window)
- [<matplotlib.lines.Line2D object at 0x...>]
- >>> plt.title("Kaiser window")
- Text(0.5, 1.0, 'Kaiser window')
- >>> plt.ylabel("Amplitude")
- Text(0, 0.5, 'Amplitude')
- >>> plt.xlabel("Sample")
- Text(0.5, 0, 'Sample')
- >>> plt.show()
- >>> plt.figure()
- <Figure size 640x480 with 0 Axes>
- >>> A = fft(window, 2048) / 25.5
- >>> mag = np.abs(fftshift(A))
- >>> freq = np.linspace(-0.5, 0.5, len(A))
- >>> response = 20 * np.log10(mag)
- >>> response = np.clip(response, -100, 100)
- >>> plt.plot(freq, response)
- [<matplotlib.lines.Line2D object at 0x...>]
- >>> plt.title("Frequency response of Kaiser window")
- Text(0.5, 1.0, 'Frequency response of Kaiser window')
- >>> plt.ylabel("Magnitude [dB]")
- Text(0, 0.5, 'Magnitude [dB]')
- >>> plt.xlabel("Normalized frequency [cycles per sample]")
- Text(0.5, 0, 'Normalized frequency [cycles per sample]')
- >>> plt.axis('tight')
- (-0.5, 0.5, -100.0, ...) # may vary
- >>> plt.show()
- """
- if M == 1:
- return np.array([1.])
- n = arange(0, M)
- alpha = (M-1)/2.0
- return i0(beta * sqrt(1-((n-alpha)/alpha)**2.0))/i0(float(beta))
- def _sinc_dispatcher(x):
- return (x,)
- @array_function_dispatch(_sinc_dispatcher)
- def sinc(x):
- r"""
- Return the normalized sinc function.
- The sinc function is :math:`\sin(\pi x)/(\pi x)`.
- .. note::
- Note the normalization factor of ``pi`` used in the definition.
- This is the most commonly used definition in signal processing.
- Use ``sinc(x / np.pi)`` to obtain the unnormalized sinc function
- :math:`\sin(x)/(x)` that is more common in mathematics.
- Parameters
- ----------
- x : ndarray
- Array (possibly multi-dimensional) of values for which to to
- calculate ``sinc(x)``.
- Returns
- -------
- out : ndarray
- ``sinc(x)``, which has the same shape as the input.
- Notes
- -----
- ``sinc(0)`` is the limit value 1.
- The name sinc is short for "sine cardinal" or "sinus cardinalis".
- The sinc function is used in various signal processing applications,
- including in anti-aliasing, in the construction of a Lanczos resampling
- filter, and in interpolation.
- For bandlimited interpolation of discrete-time signals, the ideal
- interpolation kernel is proportional to the sinc function.
- References
- ----------
- .. [1] Weisstein, Eric W. "Sinc Function." From MathWorld--A Wolfram Web
- Resource. http://mathworld.wolfram.com/SincFunction.html
- .. [2] Wikipedia, "Sinc function",
- https://en.wikipedia.org/wiki/Sinc_function
- Examples
- --------
- >>> import matplotlib.pyplot as plt
- >>> x = np.linspace(-4, 4, 41)
- >>> np.sinc(x)
- array([-3.89804309e-17, -4.92362781e-02, -8.40918587e-02, # may vary
- -8.90384387e-02, -5.84680802e-02, 3.89804309e-17,
- 6.68206631e-02, 1.16434881e-01, 1.26137788e-01,
- 8.50444803e-02, -3.89804309e-17, -1.03943254e-01,
- -1.89206682e-01, -2.16236208e-01, -1.55914881e-01,
- 3.89804309e-17, 2.33872321e-01, 5.04551152e-01,
- 7.56826729e-01, 9.35489284e-01, 1.00000000e+00,
- 9.35489284e-01, 7.56826729e-01, 5.04551152e-01,
- 2.33872321e-01, 3.89804309e-17, -1.55914881e-01,
- -2.16236208e-01, -1.89206682e-01, -1.03943254e-01,
- -3.89804309e-17, 8.50444803e-02, 1.26137788e-01,
- 1.16434881e-01, 6.68206631e-02, 3.89804309e-17,
- -5.84680802e-02, -8.90384387e-02, -8.40918587e-02,
- -4.92362781e-02, -3.89804309e-17])
- >>> plt.plot(x, np.sinc(x))
- [<matplotlib.lines.Line2D object at 0x...>]
- >>> plt.title("Sinc Function")
- Text(0.5, 1.0, 'Sinc Function')
- >>> plt.ylabel("Amplitude")
- Text(0, 0.5, 'Amplitude')
- >>> plt.xlabel("X")
- Text(0.5, 0, 'X')
- >>> plt.show()
- """
- x = np.asanyarray(x)
- y = pi * where(x == 0, 1.0e-20, x)
- return sin(y)/y
- def _msort_dispatcher(a):
- return (a,)
- @array_function_dispatch(_msort_dispatcher)
- def msort(a):
- """
- Return a copy of an array sorted along the first axis.
- Parameters
- ----------
- a : array_like
- Array to be sorted.
- Returns
- -------
- sorted_array : ndarray
- Array of the same type and shape as `a`.
- See Also
- --------
- sort
- Notes
- -----
- ``np.msort(a)`` is equivalent to ``np.sort(a, axis=0)``.
- """
- b = array(a, subok=True, copy=True)
- b.sort(0)
- return b
- def _ureduce(a, func, **kwargs):
- """
- Internal Function.
- Call `func` with `a` as first argument swapping the axes to use extended
- axis on functions that don't support it natively.
- Returns result and a.shape with axis dims set to 1.
- Parameters
- ----------
- a : array_like
- Input array or object that can be converted to an array.
- func : callable
- Reduction function capable of receiving a single axis argument.
- It is called with `a` as first argument followed by `kwargs`.
- kwargs : keyword arguments
- additional keyword arguments to pass to `func`.
- Returns
- -------
- result : tuple
- Result of func(a, **kwargs) and a.shape with axis dims set to 1
- which can be used to reshape the result to the same shape a ufunc with
- keepdims=True would produce.
- """
- a = np.asanyarray(a)
- axis = kwargs.get('axis', None)
- if axis is not None:
- keepdim = list(a.shape)
- nd = a.ndim
- axis = _nx.normalize_axis_tuple(axis, nd)
- for ax in axis:
- keepdim[ax] = 1
- if len(axis) == 1:
- kwargs['axis'] = axis[0]
- else:
- keep = set(range(nd)) - set(axis)
- nkeep = len(keep)
- # swap axis that should not be reduced to front
- for i, s in enumerate(sorted(keep)):
- a = a.swapaxes(i, s)
- # merge reduced axis
- a = a.reshape(a.shape[:nkeep] + (-1,))
- kwargs['axis'] = -1
- keepdim = tuple(keepdim)
- else:
- keepdim = (1,) * a.ndim
- r = func(a, **kwargs)
- return r, keepdim
- def _median_dispatcher(
- a, axis=None, out=None, overwrite_input=None, keepdims=None):
- return (a, out)
- @array_function_dispatch(_median_dispatcher)
- def median(a, axis=None, out=None, overwrite_input=False, keepdims=False):
- """
- Compute the median along the specified axis.
- Returns the median of the array elements.
- Parameters
- ----------
- a : array_like
- Input array or object that can be converted to an array.
- axis : {int, sequence of int, None}, optional
- Axis or axes along which the medians are computed. The default
- is to compute the median along a flattened version of the array.
- A sequence of axes is supported since version 1.9.0.
- out : ndarray, optional
- Alternative output array in which to place the result. It must
- have the same shape and buffer length as the expected output,
- but the type (of the output) will be cast if necessary.
- overwrite_input : bool, optional
- If True, then allow use of memory of input array `a` for
- calculations. The input array will be modified by the call to
- `median`. This will save memory when you do not need to preserve
- the contents of the input array. Treat the input as undefined,
- but it will probably be fully or partially sorted. Default is
- False. If `overwrite_input` is ``True`` and `a` is not already an
- `ndarray`, an error will be raised.
- keepdims : bool, optional
- If this is set to True, the axes which are reduced are left
- in the result as dimensions with size one. With this option,
- the result will broadcast correctly against the original `arr`.
- .. versionadded:: 1.9.0
- Returns
- -------
- median : ndarray
- A new array holding the result. If the input contains integers
- or floats smaller than ``float64``, then the output data-type is
- ``np.float64``. Otherwise, the data-type of the output is the
- same as that of the input. If `out` is specified, that array is
- returned instead.
- See Also
- --------
- mean, percentile
- Notes
- -----
- Given a vector ``V`` of length ``N``, the median of ``V`` is the
- middle value of a sorted copy of ``V``, ``V_sorted`` - i
- e., ``V_sorted[(N-1)/2]``, when ``N`` is odd, and the average of the
- two middle values of ``V_sorted`` when ``N`` is even.
- Examples
- --------
- >>> a = np.array([[10, 7, 4], [3, 2, 1]])
- >>> a
- array([[10, 7, 4],
- [ 3, 2, 1]])
- >>> np.median(a)
- 3.5
- >>> np.median(a, axis=0)
- array([6.5, 4.5, 2.5])
- >>> np.median(a, axis=1)
- array([7., 2.])
- >>> m = np.median(a, axis=0)
- >>> out = np.zeros_like(m)
- >>> np.median(a, axis=0, out=m)
- array([6.5, 4.5, 2.5])
- >>> m
- array([6.5, 4.5, 2.5])
- >>> b = a.copy()
- >>> np.median(b, axis=1, overwrite_input=True)
- array([7., 2.])
- >>> assert not np.all(a==b)
- >>> b = a.copy()
- >>> np.median(b, axis=None, overwrite_input=True)
- 3.5
- >>> assert not np.all(a==b)
- """
- r, k = _ureduce(a, func=_median, axis=axis, out=out,
- overwrite_input=overwrite_input)
- if keepdims:
- return r.reshape(k)
- else:
- return r
- def _median(a, axis=None, out=None, overwrite_input=False):
- # can't be reasonably be implemented in terms of percentile as we have to
- # call mean to not break astropy
- a = np.asanyarray(a)
- # Set the partition indexes
- if axis is None:
- sz = a.size
- else:
- sz = a.shape[axis]
- if sz % 2 == 0:
- szh = sz // 2
- kth = [szh - 1, szh]
- else:
- kth = [(sz - 1) // 2]
- # Check if the array contains any nan's
- if np.issubdtype(a.dtype, np.inexact):
- kth.append(-1)
- if overwrite_input:
- if axis is None:
- part = a.ravel()
- part.partition(kth)
- else:
- a.partition(kth, axis=axis)
- part = a
- else:
- part = partition(a, kth, axis=axis)
- if part.shape == ():
- # make 0-D arrays work
- return part.item()
- if axis is None:
- axis = 0
- indexer = [slice(None)] * part.ndim
- index = part.shape[axis] // 2
- if part.shape[axis] % 2 == 1:
- # index with slice to allow mean (below) to work
- indexer[axis] = slice(index, index+1)
- else:
- indexer[axis] = slice(index-1, index+1)
- indexer = tuple(indexer)
- # Check if the array contains any nan's
- if np.issubdtype(a.dtype, np.inexact) and sz > 0:
- # warn and return nans like mean would
- rout = mean(part[indexer], axis=axis, out=out)
- return np.lib.utils._median_nancheck(part, rout, axis, out)
- else:
- # if there are no nans
- # Use mean in odd and even case to coerce data type
- # and check, use out array.
- return mean(part[indexer], axis=axis, out=out)
- def _percentile_dispatcher(a, q, axis=None, out=None, overwrite_input=None,
- interpolation=None, keepdims=None):
- return (a, q, out)
- @array_function_dispatch(_percentile_dispatcher)
- def percentile(a, q, axis=None, out=None,
- overwrite_input=False, interpolation='linear', keepdims=False):
- """
- Compute the q-th percentile of the data along the specified axis.
- Returns the q-th percentile(s) of the array elements.
- Parameters
- ----------
- a : array_like
- Input array or object that can be converted to an array.
- q : array_like of float
- Percentile or sequence of percentiles to compute, which must be between
- 0 and 100 inclusive.
- axis : {int, tuple of int, None}, optional
- Axis or axes along which the percentiles are computed. The
- default is to compute the percentile(s) along a flattened
- version of the array.
- .. versionchanged:: 1.9.0
- A tuple of axes is supported
- out : ndarray, optional
- Alternative output array in which to place the result. It must
- have the same shape and buffer length as the expected output,
- but the type (of the output) will be cast if necessary.
- overwrite_input : bool, optional
- If True, then allow the input array `a` to be modified by intermediate
- calculations, to save memory. In this case, the contents of the input
- `a` after this function completes is undefined.
- interpolation : {'linear', 'lower', 'higher', 'midpoint', 'nearest'}
- This optional parameter specifies the interpolation method to
- use when the desired percentile lies between two data points
- ``i < j``:
- * 'linear': ``i + (j - i) * fraction``, where ``fraction``
- is the fractional part of the index surrounded by ``i``
- and ``j``.
- * 'lower': ``i``.
- * 'higher': ``j``.
- * 'nearest': ``i`` or ``j``, whichever is nearest.
- * 'midpoint': ``(i + j) / 2``.
- .. versionadded:: 1.9.0
- keepdims : bool, optional
- If this is set to True, the axes which are reduced are left in
- the result as dimensions with size one. With this option, the
- result will broadcast correctly against the original array `a`.
- .. versionadded:: 1.9.0
- Returns
- -------
- percentile : scalar or ndarray
- If `q` is a single percentile and `axis=None`, then the result
- is a scalar. If multiple percentiles are given, first axis of
- the result corresponds to the percentiles. The other axes are
- the axes that remain after the reduction of `a`. If the input
- contains integers or floats smaller than ``float64``, the output
- data-type is ``float64``. Otherwise, the output data-type is the
- same as that of the input. If `out` is specified, that array is
- returned instead.
- See Also
- --------
- mean
- median : equivalent to ``percentile(..., 50)``
- nanpercentile
- quantile : equivalent to percentile, except with q in the range [0, 1].
- Notes
- -----
- Given a vector ``V`` of length ``N``, the q-th percentile of
- ``V`` is the value ``q/100`` of the way from the minimum to the
- maximum in a sorted copy of ``V``. The values and distances of
- the two nearest neighbors as well as the `interpolation` parameter
- will determine the percentile if the normalized ranking does not
- match the location of ``q`` exactly. This function is the same as
- the median if ``q=50``, the same as the minimum if ``q=0`` and the
- same as the maximum if ``q=100``.
- Examples
- --------
- >>> a = np.array([[10, 7, 4], [3, 2, 1]])
- >>> a
- array([[10, 7, 4],
- [ 3, 2, 1]])
- >>> np.percentile(a, 50)
- 3.5
- >>> np.percentile(a, 50, axis=0)
- array([6.5, 4.5, 2.5])
- >>> np.percentile(a, 50, axis=1)
- array([7., 2.])
- >>> np.percentile(a, 50, axis=1, keepdims=True)
- array([[7.],
- [2.]])
- >>> m = np.percentile(a, 50, axis=0)
- >>> out = np.zeros_like(m)
- >>> np.percentile(a, 50, axis=0, out=out)
- array([6.5, 4.5, 2.5])
- >>> m
- array([6.5, 4.5, 2.5])
- >>> b = a.copy()
- >>> np.percentile(b, 50, axis=1, overwrite_input=True)
- array([7., 2.])
- >>> assert not np.all(a == b)
- The different types of interpolation can be visualized graphically:
- .. plot::
- import matplotlib.pyplot as plt
- a = np.arange(4)
- p = np.linspace(0, 100, 6001)
- ax = plt.gca()
- lines = [
- ('linear', None),
- ('higher', '--'),
- ('lower', '--'),
- ('nearest', '-.'),
- ('midpoint', '-.'),
- ]
- for interpolation, style in lines:
- ax.plot(
- p, np.percentile(a, p, interpolation=interpolation),
- label=interpolation, linestyle=style)
- ax.set(
- title='Interpolation methods for list: ' + str(a),
- xlabel='Percentile',
- ylabel='List item returned',
- yticks=a)
- ax.legend()
- plt.show()
- """
- q = np.true_divide(q, 100)
- q = asanyarray(q) # undo any decay that the ufunc performed (see gh-13105)
- if not _quantile_is_valid(q):
- raise ValueError("Percentiles must be in the range [0, 100]")
- return _quantile_unchecked(
- a, q, axis, out, overwrite_input, interpolation, keepdims)
- def _quantile_dispatcher(a, q, axis=None, out=None, overwrite_input=None,
- interpolation=None, keepdims=None):
- return (a, q, out)
- @array_function_dispatch(_quantile_dispatcher)
- def quantile(a, q, axis=None, out=None,
- overwrite_input=False, interpolation='linear', keepdims=False):
- """
- Compute the q-th quantile of the data along the specified axis.
- .. versionadded:: 1.15.0
- Parameters
- ----------
- a : array_like
- Input array or object that can be converted to an array.
- q : array_like of float
- Quantile or sequence of quantiles to compute, which must be between
- 0 and 1 inclusive.
- axis : {int, tuple of int, None}, optional
- Axis or axes along which the quantiles are computed. The
- default is to compute the quantile(s) along a flattened
- version of the array.
- out : ndarray, optional
- Alternative output array in which to place the result. It must
- have the same shape and buffer length as the expected output,
- but the type (of the output) will be cast if necessary.
- overwrite_input : bool, optional
- If True, then allow the input array `a` to be modified by intermediate
- calculations, to save memory. In this case, the contents of the input
- `a` after this function completes is undefined.
- interpolation : {'linear', 'lower', 'higher', 'midpoint', 'nearest'}
- This optional parameter specifies the interpolation method to
- use when the desired quantile lies between two data points
- ``i < j``:
- * linear: ``i + (j - i) * fraction``, where ``fraction``
- is the fractional part of the index surrounded by ``i``
- and ``j``.
- * lower: ``i``.
- * higher: ``j``.
- * nearest: ``i`` or ``j``, whichever is nearest.
- * midpoint: ``(i + j) / 2``.
- keepdims : bool, optional
- If this is set to True, the axes which are reduced are left in
- the result as dimensions with size one. With this option, the
- result will broadcast correctly against the original array `a`.
- Returns
- -------
- quantile : scalar or ndarray
- If `q` is a single quantile and `axis=None`, then the result
- is a scalar. If multiple quantiles are given, first axis of
- the result corresponds to the quantiles. The other axes are
- the axes that remain after the reduction of `a`. If the input
- contains integers or floats smaller than ``float64``, the output
- data-type is ``float64``. Otherwise, the output data-type is the
- same as that of the input. If `out` is specified, that array is
- returned instead.
- See Also
- --------
- mean
- percentile : equivalent to quantile, but with q in the range [0, 100].
- median : equivalent to ``quantile(..., 0.5)``
- nanquantile
- Notes
- -----
- Given a vector ``V`` of length ``N``, the q-th quantile of
- ``V`` is the value ``q`` of the way from the minimum to the
- maximum in a sorted copy of ``V``. The values and distances of
- the two nearest neighbors as well as the `interpolation` parameter
- will determine the quantile if the normalized ranking does not
- match the location of ``q`` exactly. This function is the same as
- the median if ``q=0.5``, the same as the minimum if ``q=0.0`` and the
- same as the maximum if ``q=1.0``.
- Examples
- --------
- >>> a = np.array([[10, 7, 4], [3, 2, 1]])
- >>> a
- array([[10, 7, 4],
- [ 3, 2, 1]])
- >>> np.quantile(a, 0.5)
- 3.5
- >>> np.quantile(a, 0.5, axis=0)
- array([6.5, 4.5, 2.5])
- >>> np.quantile(a, 0.5, axis=1)
- array([7., 2.])
- >>> np.quantile(a, 0.5, axis=1, keepdims=True)
- array([[7.],
- [2.]])
- >>> m = np.quantile(a, 0.5, axis=0)
- >>> out = np.zeros_like(m)
- >>> np.quantile(a, 0.5, axis=0, out=out)
- array([6.5, 4.5, 2.5])
- >>> m
- array([6.5, 4.5, 2.5])
- >>> b = a.copy()
- >>> np.quantile(b, 0.5, axis=1, overwrite_input=True)
- array([7., 2.])
- >>> assert not np.all(a == b)
- """
- q = np.asanyarray(q)
- if not _quantile_is_valid(q):
- raise ValueError("Quantiles must be in the range [0, 1]")
- return _quantile_unchecked(
- a, q, axis, out, overwrite_input, interpolation, keepdims)
- def _quantile_unchecked(a, q, axis=None, out=None, overwrite_input=False,
- interpolation='linear', keepdims=False):
- """Assumes that q is in [0, 1], and is an ndarray"""
- r, k = _ureduce(a, func=_quantile_ureduce_func, q=q, axis=axis, out=out,
- overwrite_input=overwrite_input,
- interpolation=interpolation)
- if keepdims:
- return r.reshape(q.shape + k)
- else:
- return r
- def _quantile_is_valid(q):
- # avoid expensive reductions, relevant for arrays with < O(1000) elements
- if q.ndim == 1 and q.size < 10:
- for i in range(q.size):
- if q[i] < 0.0 or q[i] > 1.0:
- return False
- else:
- # faster than any()
- if np.count_nonzero(q < 0.0) or np.count_nonzero(q > 1.0):
- return False
- return True
- def _lerp(a, b, t, out=None):
- """ Linearly interpolate from a to b by a factor of t """
- diff_b_a = subtract(b, a)
- # asanyarray is a stop-gap until gh-13105
- lerp_interpolation = asanyarray(add(a, diff_b_a*t, out=out))
- subtract(b, diff_b_a * (1 - t), out=lerp_interpolation, where=t>=0.5)
- if lerp_interpolation.ndim == 0 and out is None:
- lerp_interpolation = lerp_interpolation[()] # unpack 0d arrays
- return lerp_interpolation
- def _quantile_ureduce_func(a, q, axis=None, out=None, overwrite_input=False,
- interpolation='linear', keepdims=False):
- a = asarray(a)
- # ufuncs cause 0d array results to decay to scalars (see gh-13105), which
- # makes them problematic for __setitem__ and attribute access. As a
- # workaround, we call this on the result of every ufunc on a possibly-0d
- # array.
- not_scalar = np.asanyarray
- # prepare a for partitioning
- if overwrite_input:
- if axis is None:
- ap = a.ravel()
- else:
- ap = a
- else:
- if axis is None:
- ap = a.flatten()
- else:
- ap = a.copy()
- if axis is None:
- axis = 0
- if q.ndim > 2:
- # The code below works fine for nd, but it might not have useful
- # semantics. For now, keep the supported dimensions the same as it was
- # before.
- raise ValueError("q must be a scalar or 1d")
- Nx = ap.shape[axis]
- indices = not_scalar(q * (Nx - 1))
- # round fractional indices according to interpolation method
- if interpolation == 'lower':
- indices = floor(indices).astype(intp)
- elif interpolation == 'higher':
- indices = ceil(indices).astype(intp)
- elif interpolation == 'midpoint':
- indices = 0.5 * (floor(indices) + ceil(indices))
- elif interpolation == 'nearest':
- indices = around(indices).astype(intp)
- elif interpolation == 'linear':
- pass # keep index as fraction and interpolate
- else:
- raise ValueError(
- "interpolation can only be 'linear', 'lower' 'higher', "
- "'midpoint', or 'nearest'")
- # The dimensions of `q` are prepended to the output shape, so we need the
- # axis being sampled from `ap` to be first.
- ap = np.moveaxis(ap, axis, 0)
- del axis
- if np.issubdtype(indices.dtype, np.integer):
- # take the points along axis
- if np.issubdtype(a.dtype, np.inexact):
- # may contain nan, which would sort to the end
- ap.partition(concatenate((indices.ravel(), [-1])), axis=0)
- n = np.isnan(ap[-1])
- else:
- # cannot contain nan
- ap.partition(indices.ravel(), axis=0)
- n = np.array(False, dtype=bool)
- r = take(ap, indices, axis=0, out=out)
- else:
- # weight the points above and below the indices
- indices_below = not_scalar(floor(indices)).astype(intp)
- indices_above = not_scalar(indices_below + 1)
- indices_above[indices_above > Nx - 1] = Nx - 1
- if np.issubdtype(a.dtype, np.inexact):
- # may contain nan, which would sort to the end
- ap.partition(concatenate((
- indices_below.ravel(), indices_above.ravel(), [-1]
- )), axis=0)
- n = np.isnan(ap[-1])
- else:
- # cannot contain nan
- ap.partition(concatenate((
- indices_below.ravel(), indices_above.ravel()
- )), axis=0)
- n = np.array(False, dtype=bool)
- weights_shape = indices.shape + (1,) * (ap.ndim - 1)
- weights_above = not_scalar(indices - indices_below).reshape(weights_shape)
- x_below = take(ap, indices_below, axis=0)
- x_above = take(ap, indices_above, axis=0)
- r = _lerp(x_below, x_above, weights_above, out=out)
- # if any slice contained a nan, then all results on that slice are also nan
- if np.any(n):
- if r.ndim == 0 and out is None:
- # can't write to a scalar
- r = a.dtype.type(np.nan)
- else:
- r[..., n] = a.dtype.type(np.nan)
- return r
- def _trapz_dispatcher(y, x=None, dx=None, axis=None):
- return (y, x)
- @array_function_dispatch(_trapz_dispatcher)
- def trapz(y, x=None, dx=1.0, axis=-1):
- """
- Integrate along the given axis using the composite trapezoidal rule.
- Integrate `y` (`x`) along given axis.
- Parameters
- ----------
- y : array_like
- Input array to integrate.
- x : array_like, optional
- The sample points corresponding to the `y` values. If `x` is None,
- the sample points are assumed to be evenly spaced `dx` apart. The
- default is None.
- dx : scalar, optional
- The spacing between sample points when `x` is None. The default is 1.
- axis : int, optional
- The axis along which to integrate.
- Returns
- -------
- trapz : float
- Definite integral as approximated by trapezoidal rule.
- See Also
- --------
- sum, cumsum
- Notes
- -----
- Image [2]_ illustrates trapezoidal rule -- y-axis locations of points
- will be taken from `y` array, by default x-axis distances between
- points will be 1.0, alternatively they can be provided with `x` array
- or with `dx` scalar. Return value will be equal to combined area under
- the red lines.
- References
- ----------
- .. [1] Wikipedia page: https://en.wikipedia.org/wiki/Trapezoidal_rule
- .. [2] Illustration image:
- https://en.wikipedia.org/wiki/File:Composite_trapezoidal_rule_illustration.png
- Examples
- --------
- >>> np.trapz([1,2,3])
- 4.0
- >>> np.trapz([1,2,3], x=[4,6,8])
- 8.0
- >>> np.trapz([1,2,3], dx=2)
- 8.0
- >>> a = np.arange(6).reshape(2, 3)
- >>> a
- array([[0, 1, 2],
- [3, 4, 5]])
- >>> np.trapz(a, axis=0)
- array([1.5, 2.5, 3.5])
- >>> np.trapz(a, axis=1)
- array([2., 8.])
- """
- y = asanyarray(y)
- if x is None:
- d = dx
- else:
- x = asanyarray(x)
- if x.ndim == 1:
- d = diff(x)
- # reshape to correct shape
- shape = [1]*y.ndim
- shape[axis] = d.shape[0]
- d = d.reshape(shape)
- else:
- d = diff(x, axis=axis)
- nd = y.ndim
- slice1 = [slice(None)]*nd
- slice2 = [slice(None)]*nd
- slice1[axis] = slice(1, None)
- slice2[axis] = slice(None, -1)
- try:
- ret = (d * (y[tuple(slice1)] + y[tuple(slice2)]) / 2.0).sum(axis)
- except ValueError:
- # Operations didn't work, cast to ndarray
- d = np.asarray(d)
- y = np.asarray(y)
- ret = add.reduce(d * (y[tuple(slice1)]+y[tuple(slice2)])/2.0, axis)
- return ret
- def _meshgrid_dispatcher(*xi, copy=None, sparse=None, indexing=None):
- return xi
- # Based on scitools meshgrid
- @array_function_dispatch(_meshgrid_dispatcher)
- def meshgrid(*xi, copy=True, sparse=False, indexing='xy'):
- """
- Return coordinate matrices from coordinate vectors.
- Make N-D coordinate arrays for vectorized evaluations of
- N-D scalar/vector fields over N-D grids, given
- one-dimensional coordinate arrays x1, x2,..., xn.
- .. versionchanged:: 1.9
- 1-D and 0-D cases are allowed.
- Parameters
- ----------
- x1, x2,..., xn : array_like
- 1-D arrays representing the coordinates of a grid.
- indexing : {'xy', 'ij'}, optional
- Cartesian ('xy', default) or matrix ('ij') indexing of output.
- See Notes for more details.
- .. versionadded:: 1.7.0
- sparse : bool, optional
- If True a sparse grid is returned in order to conserve memory.
- Default is False.
- .. versionadded:: 1.7.0
- copy : bool, optional
- If False, a view into the original arrays are returned in order to
- conserve memory. Default is True. Please note that
- ``sparse=False, copy=False`` will likely return non-contiguous
- arrays. Furthermore, more than one element of a broadcast array
- may refer to a single memory location. If you need to write to the
- arrays, make copies first.
- .. versionadded:: 1.7.0
- Returns
- -------
- X1, X2,..., XN : ndarray
- For vectors `x1`, `x2`,..., 'xn' with lengths ``Ni=len(xi)`` ,
- return ``(N1, N2, N3,...Nn)`` shaped arrays if indexing='ij'
- or ``(N2, N1, N3,...Nn)`` shaped arrays if indexing='xy'
- with the elements of `xi` repeated to fill the matrix along
- the first dimension for `x1`, the second for `x2` and so on.
- Notes
- -----
- This function supports both indexing conventions through the indexing
- keyword argument. Giving the string 'ij' returns a meshgrid with
- matrix indexing, while 'xy' returns a meshgrid with Cartesian indexing.
- In the 2-D case with inputs of length M and N, the outputs are of shape
- (N, M) for 'xy' indexing and (M, N) for 'ij' indexing. In the 3-D case
- with inputs of length M, N and P, outputs are of shape (N, M, P) for
- 'xy' indexing and (M, N, P) for 'ij' indexing. The difference is
- illustrated by the following code snippet::
- xv, yv = np.meshgrid(x, y, sparse=False, indexing='ij')
- for i in range(nx):
- for j in range(ny):
- # treat xv[i,j], yv[i,j]
- xv, yv = np.meshgrid(x, y, sparse=False, indexing='xy')
- for i in range(nx):
- for j in range(ny):
- # treat xv[j,i], yv[j,i]
- In the 1-D and 0-D case, the indexing and sparse keywords have no effect.
- See Also
- --------
- mgrid : Construct a multi-dimensional "meshgrid" using indexing notation.
- ogrid : Construct an open multi-dimensional "meshgrid" using indexing
- notation.
- Examples
- --------
- >>> nx, ny = (3, 2)
- >>> x = np.linspace(0, 1, nx)
- >>> y = np.linspace(0, 1, ny)
- >>> xv, yv = np.meshgrid(x, y)
- >>> xv
- array([[0. , 0.5, 1. ],
- [0. , 0.5, 1. ]])
- >>> yv
- array([[0., 0., 0.],
- [1., 1., 1.]])
- >>> xv, yv = np.meshgrid(x, y, sparse=True) # make sparse output arrays
- >>> xv
- array([[0. , 0.5, 1. ]])
- >>> yv
- array([[0.],
- [1.]])
- `meshgrid` is very useful to evaluate functions on a grid.
- >>> import matplotlib.pyplot as plt
- >>> x = np.arange(-5, 5, 0.1)
- >>> y = np.arange(-5, 5, 0.1)
- >>> xx, yy = np.meshgrid(x, y, sparse=True)
- >>> z = np.sin(xx**2 + yy**2) / (xx**2 + yy**2)
- >>> h = plt.contourf(x,y,z)
- >>> plt.show()
- """
- ndim = len(xi)
- if indexing not in ['xy', 'ij']:
- raise ValueError(
- "Valid values for `indexing` are 'xy' and 'ij'.")
- s0 = (1,) * ndim
- output = [np.asanyarray(x).reshape(s0[:i] + (-1,) + s0[i + 1:])
- for i, x in enumerate(xi)]
- if indexing == 'xy' and ndim > 1:
- # switch first and second axis
- output[0].shape = (1, -1) + s0[2:]
- output[1].shape = (-1, 1) + s0[2:]
- if not sparse:
- # Return the full N-D matrix (not only the 1-D vector)
- output = np.broadcast_arrays(*output, subok=True)
- if copy:
- output = [x.copy() for x in output]
- return output
- def _delete_dispatcher(arr, obj, axis=None):
- return (arr, obj)
- @array_function_dispatch(_delete_dispatcher)
- def delete(arr, obj, axis=None):
- """
- Return a new array with sub-arrays along an axis deleted. For a one
- dimensional array, this returns those entries not returned by
- `arr[obj]`.
- Parameters
- ----------
- arr : array_like
- Input array.
- obj : slice, int or array of ints
- Indicate indices of sub-arrays to remove along the specified axis.
- .. versionchanged:: 1.19.0
- Boolean indices are now treated as a mask of elements to remove,
- rather than being cast to the integers 0 and 1.
- axis : int, optional
- The axis along which to delete the subarray defined by `obj`.
- If `axis` is None, `obj` is applied to the flattened array.
- Returns
- -------
- out : ndarray
- A copy of `arr` with the elements specified by `obj` removed. Note
- that `delete` does not occur in-place. If `axis` is None, `out` is
- a flattened array.
- See Also
- --------
- insert : Insert elements into an array.
- append : Append elements at the end of an array.
- Notes
- -----
- Often it is preferable to use a boolean mask. For example:
- >>> arr = np.arange(12) + 1
- >>> mask = np.ones(len(arr), dtype=bool)
- >>> mask[[0,2,4]] = False
- >>> result = arr[mask,...]
- Is equivalent to `np.delete(arr, [0,2,4], axis=0)`, but allows further
- use of `mask`.
- Examples
- --------
- >>> arr = np.array([[1,2,3,4], [5,6,7,8], [9,10,11,12]])
- >>> arr
- array([[ 1, 2, 3, 4],
- [ 5, 6, 7, 8],
- [ 9, 10, 11, 12]])
- >>> np.delete(arr, 1, 0)
- array([[ 1, 2, 3, 4],
- [ 9, 10, 11, 12]])
- >>> np.delete(arr, np.s_[::2], 1)
- array([[ 2, 4],
- [ 6, 8],
- [10, 12]])
- >>> np.delete(arr, [1,3,5], None)
- array([ 1, 3, 5, 7, 8, 9, 10, 11, 12])
- """
- wrap = None
- if type(arr) is not ndarray:
- try:
- wrap = arr.__array_wrap__
- except AttributeError:
- pass
- arr = asarray(arr)
- ndim = arr.ndim
- arrorder = 'F' if arr.flags.fnc else 'C'
- if axis is None:
- if ndim != 1:
- arr = arr.ravel()
- # needed for np.matrix, which is still not 1d after being ravelled
- ndim = arr.ndim
- axis = ndim - 1
- else:
- axis = normalize_axis_index(axis, ndim)
- slobj = [slice(None)]*ndim
- N = arr.shape[axis]
- newshape = list(arr.shape)
- if isinstance(obj, slice):
- start, stop, step = obj.indices(N)
- xr = range(start, stop, step)
- numtodel = len(xr)
- if numtodel <= 0:
- if wrap:
- return wrap(arr.copy(order=arrorder))
- else:
- return arr.copy(order=arrorder)
- # Invert if step is negative:
- if step < 0:
- step = -step
- start = xr[-1]
- stop = xr[0] + 1
- newshape[axis] -= numtodel
- new = empty(newshape, arr.dtype, arrorder)
- # copy initial chunk
- if start == 0:
- pass
- else:
- slobj[axis] = slice(None, start)
- new[tuple(slobj)] = arr[tuple(slobj)]
- # copy end chunk
- if stop == N:
- pass
- else:
- slobj[axis] = slice(stop-numtodel, None)
- slobj2 = [slice(None)]*ndim
- slobj2[axis] = slice(stop, None)
- new[tuple(slobj)] = arr[tuple(slobj2)]
- # copy middle pieces
- if step == 1:
- pass
- else: # use array indexing.
- keep = ones(stop-start, dtype=bool)
- keep[:stop-start:step] = False
- slobj[axis] = slice(start, stop-numtodel)
- slobj2 = [slice(None)]*ndim
- slobj2[axis] = slice(start, stop)
- arr = arr[tuple(slobj2)]
- slobj2[axis] = keep
- new[tuple(slobj)] = arr[tuple(slobj2)]
- if wrap:
- return wrap(new)
- else:
- return new
- if isinstance(obj, (int, integer)) and not isinstance(obj, bool):
- # optimization for a single value
- if (obj < -N or obj >= N):
- raise IndexError(
- "index %i is out of bounds for axis %i with "
- "size %i" % (obj, axis, N))
- if (obj < 0):
- obj += N
- newshape[axis] -= 1
- new = empty(newshape, arr.dtype, arrorder)
- slobj[axis] = slice(None, obj)
- new[tuple(slobj)] = arr[tuple(slobj)]
- slobj[axis] = slice(obj, None)
- slobj2 = [slice(None)]*ndim
- slobj2[axis] = slice(obj+1, None)
- new[tuple(slobj)] = arr[tuple(slobj2)]
- else:
- _obj = obj
- obj = np.asarray(obj)
- if obj.size == 0 and not isinstance(_obj, np.ndarray):
- obj = obj.astype(intp)
- if obj.dtype == bool:
- if obj.shape != (N,):
- raise ValueError('boolean array argument obj to delete '
- 'must be one dimensional and match the axis '
- 'length of {}'.format(N))
- # optimization, the other branch is slower
- keep = ~obj
- else:
- keep = ones(N, dtype=bool)
- keep[obj,] = False
- slobj[axis] = keep
- new = arr[tuple(slobj)]
- if wrap:
- return wrap(new)
- else:
- return new
- def _insert_dispatcher(arr, obj, values, axis=None):
- return (arr, obj, values)
- @array_function_dispatch(_insert_dispatcher)
- def insert(arr, obj, values, axis=None):
- """
- Insert values along the given axis before the given indices.
- Parameters
- ----------
- arr : array_like
- Input array.
- obj : int, slice or sequence of ints
- Object that defines the index or indices before which `values` is
- inserted.
- .. versionadded:: 1.8.0
- Support for multiple insertions when `obj` is a single scalar or a
- sequence with one element (similar to calling insert multiple
- times).
- values : array_like
- Values to insert into `arr`. If the type of `values` is different
- from that of `arr`, `values` is converted to the type of `arr`.
- `values` should be shaped so that ``arr[...,obj,...] = values``
- is legal.
- axis : int, optional
- Axis along which to insert `values`. If `axis` is None then `arr`
- is flattened first.
- Returns
- -------
- out : ndarray
- A copy of `arr` with `values` inserted. Note that `insert`
- does not occur in-place: a new array is returned. If
- `axis` is None, `out` is a flattened array.
- See Also
- --------
- append : Append elements at the end of an array.
- concatenate : Join a sequence of arrays along an existing axis.
- delete : Delete elements from an array.
- Notes
- -----
- Note that for higher dimensional inserts `obj=0` behaves very different
- from `obj=[0]` just like `arr[:,0,:] = values` is different from
- `arr[:,[0],:] = values`.
- Examples
- --------
- >>> a = np.array([[1, 1], [2, 2], [3, 3]])
- >>> a
- array([[1, 1],
- [2, 2],
- [3, 3]])
- >>> np.insert(a, 1, 5)
- array([1, 5, 1, ..., 2, 3, 3])
- >>> np.insert(a, 1, 5, axis=1)
- array([[1, 5, 1],
- [2, 5, 2],
- [3, 5, 3]])
- Difference between sequence and scalars:
- >>> np.insert(a, [1], [[1],[2],[3]], axis=1)
- array([[1, 1, 1],
- [2, 2, 2],
- [3, 3, 3]])
- >>> np.array_equal(np.insert(a, 1, [1, 2, 3], axis=1),
- ... np.insert(a, [1], [[1],[2],[3]], axis=1))
- True
- >>> b = a.flatten()
- >>> b
- array([1, 1, 2, 2, 3, 3])
- >>> np.insert(b, [2, 2], [5, 6])
- array([1, 1, 5, ..., 2, 3, 3])
- >>> np.insert(b, slice(2, 4), [5, 6])
- array([1, 1, 5, ..., 2, 3, 3])
- >>> np.insert(b, [2, 2], [7.13, False]) # type casting
- array([1, 1, 7, ..., 2, 3, 3])
- >>> x = np.arange(8).reshape(2, 4)
- >>> idx = (1, 3)
- >>> np.insert(x, idx, 999, axis=1)
- array([[ 0, 999, 1, 2, 999, 3],
- [ 4, 999, 5, 6, 999, 7]])
- """
- wrap = None
- if type(arr) is not ndarray:
- try:
- wrap = arr.__array_wrap__
- except AttributeError:
- pass
- arr = asarray(arr)
- ndim = arr.ndim
- arrorder = 'F' if arr.flags.fnc else 'C'
- if axis is None:
- if ndim != 1:
- arr = arr.ravel()
- # needed for np.matrix, which is still not 1d after being ravelled
- ndim = arr.ndim
- axis = ndim - 1
- else:
- axis = normalize_axis_index(axis, ndim)
- slobj = [slice(None)]*ndim
- N = arr.shape[axis]
- newshape = list(arr.shape)
- if isinstance(obj, slice):
- # turn it into a range object
- indices = arange(*obj.indices(N), dtype=intp)
- else:
- # need to copy obj, because indices will be changed in-place
- indices = np.array(obj)
- if indices.dtype == bool:
- # See also delete
- # 2012-10-11, NumPy 1.8
- warnings.warn(
- "in the future insert will treat boolean arrays and "
- "array-likes as a boolean index instead of casting it to "
- "integer", FutureWarning, stacklevel=3)
- indices = indices.astype(intp)
- # Code after warning period:
- #if obj.ndim != 1:
- # raise ValueError('boolean array argument obj to insert '
- # 'must be one dimensional')
- #indices = np.flatnonzero(obj)
- elif indices.ndim > 1:
- raise ValueError(
- "index array argument obj to insert must be one dimensional "
- "or scalar")
- if indices.size == 1:
- index = indices.item()
- if index < -N or index > N:
- raise IndexError(
- "index %i is out of bounds for axis %i with "
- "size %i" % (obj, axis, N))
- if (index < 0):
- index += N
- # There are some object array corner cases here, but we cannot avoid
- # that:
- values = array(values, copy=False, ndmin=arr.ndim, dtype=arr.dtype)
- if indices.ndim == 0:
- # broadcasting is very different here, since a[:,0,:] = ... behaves
- # very different from a[:,[0],:] = ...! This changes values so that
- # it works likes the second case. (here a[:,0:1,:])
- values = np.moveaxis(values, 0, axis)
- numnew = values.shape[axis]
- newshape[axis] += numnew
- new = empty(newshape, arr.dtype, arrorder)
- slobj[axis] = slice(None, index)
- new[tuple(slobj)] = arr[tuple(slobj)]
- slobj[axis] = slice(index, index+numnew)
- new[tuple(slobj)] = values
- slobj[axis] = slice(index+numnew, None)
- slobj2 = [slice(None)] * ndim
- slobj2[axis] = slice(index, None)
- new[tuple(slobj)] = arr[tuple(slobj2)]
- if wrap:
- return wrap(new)
- return new
- elif indices.size == 0 and not isinstance(obj, np.ndarray):
- # Can safely cast the empty list to intp
- indices = indices.astype(intp)
- indices[indices < 0] += N
- numnew = len(indices)
- order = indices.argsort(kind='mergesort') # stable sort
- indices[order] += np.arange(numnew)
- newshape[axis] += numnew
- old_mask = ones(newshape[axis], dtype=bool)
- old_mask[indices] = False
- new = empty(newshape, arr.dtype, arrorder)
- slobj2 = [slice(None)]*ndim
- slobj[axis] = indices
- slobj2[axis] = old_mask
- new[tuple(slobj)] = values
- new[tuple(slobj2)] = arr
- if wrap:
- return wrap(new)
- return new
- def _append_dispatcher(arr, values, axis=None):
- return (arr, values)
- @array_function_dispatch(_append_dispatcher)
- def append(arr, values, axis=None):
- """
- Append values to the end of an array.
- Parameters
- ----------
- arr : array_like
- Values are appended to a copy of this array.
- values : array_like
- These values are appended to a copy of `arr`. It must be of the
- correct shape (the same shape as `arr`, excluding `axis`). If
- `axis` is not specified, `values` can be any shape and will be
- flattened before use.
- axis : int, optional
- The axis along which `values` are appended. If `axis` is not
- given, both `arr` and `values` are flattened before use.
- Returns
- -------
- append : ndarray
- A copy of `arr` with `values` appended to `axis`. Note that
- `append` does not occur in-place: a new array is allocated and
- filled. If `axis` is None, `out` is a flattened array.
- See Also
- --------
- insert : Insert elements into an array.
- delete : Delete elements from an array.
- Examples
- --------
- >>> np.append([1, 2, 3], [[4, 5, 6], [7, 8, 9]])
- array([1, 2, 3, ..., 7, 8, 9])
- When `axis` is specified, `values` must have the correct shape.
- >>> np.append([[1, 2, 3], [4, 5, 6]], [[7, 8, 9]], axis=0)
- array([[1, 2, 3],
- [4, 5, 6],
- [7, 8, 9]])
- >>> np.append([[1, 2, 3], [4, 5, 6]], [7, 8, 9], axis=0)
- Traceback (most recent call last):
- ...
- ValueError: all the input arrays must have same number of dimensions, but
- the array at index 0 has 2 dimension(s) and the array at index 1 has 1
- dimension(s)
- """
- arr = asanyarray(arr)
- if axis is None:
- if arr.ndim != 1:
- arr = arr.ravel()
- values = ravel(values)
- axis = arr.ndim-1
- return concatenate((arr, values), axis=axis)
- def _digitize_dispatcher(x, bins, right=None):
- return (x, bins)
- @array_function_dispatch(_digitize_dispatcher)
- def digitize(x, bins, right=False):
- """
- Return the indices of the bins to which each value in input array belongs.
- ========= ============= ============================
- `right` order of bins returned index `i` satisfies
- ========= ============= ============================
- ``False`` increasing ``bins[i-1] <= x < bins[i]``
- ``True`` increasing ``bins[i-1] < x <= bins[i]``
- ``False`` decreasing ``bins[i-1] > x >= bins[i]``
- ``True`` decreasing ``bins[i-1] >= x > bins[i]``
- ========= ============= ============================
- If values in `x` are beyond the bounds of `bins`, 0 or ``len(bins)`` is
- returned as appropriate.
- Parameters
- ----------
- x : array_like
- Input array to be binned. Prior to NumPy 1.10.0, this array had to
- be 1-dimensional, but can now have any shape.
- bins : array_like
- Array of bins. It has to be 1-dimensional and monotonic.
- right : bool, optional
- Indicating whether the intervals include the right or the left bin
- edge. Default behavior is (right==False) indicating that the interval
- does not include the right edge. The left bin end is open in this
- case, i.e., bins[i-1] <= x < bins[i] is the default behavior for
- monotonically increasing bins.
- Returns
- -------
- indices : ndarray of ints
- Output array of indices, of same shape as `x`.
- Raises
- ------
- ValueError
- If `bins` is not monotonic.
- TypeError
- If the type of the input is complex.
- See Also
- --------
- bincount, histogram, unique, searchsorted
- Notes
- -----
- If values in `x` are such that they fall outside the bin range,
- attempting to index `bins` with the indices that `digitize` returns
- will result in an IndexError.
- .. versionadded:: 1.10.0
- `np.digitize` is implemented in terms of `np.searchsorted`. This means
- that a binary search is used to bin the values, which scales much better
- for larger number of bins than the previous linear search. It also removes
- the requirement for the input array to be 1-dimensional.
- For monotonically _increasing_ `bins`, the following are equivalent::
- np.digitize(x, bins, right=True)
- np.searchsorted(bins, x, side='left')
- Note that as the order of the arguments are reversed, the side must be too.
- The `searchsorted` call is marginally faster, as it does not do any
- monotonicity checks. Perhaps more importantly, it supports all dtypes.
- Examples
- --------
- >>> x = np.array([0.2, 6.4, 3.0, 1.6])
- >>> bins = np.array([0.0, 1.0, 2.5, 4.0, 10.0])
- >>> inds = np.digitize(x, bins)
- >>> inds
- array([1, 4, 3, 2])
- >>> for n in range(x.size):
- ... print(bins[inds[n]-1], "<=", x[n], "<", bins[inds[n]])
- ...
- 0.0 <= 0.2 < 1.0
- 4.0 <= 6.4 < 10.0
- 2.5 <= 3.0 < 4.0
- 1.0 <= 1.6 < 2.5
- >>> x = np.array([1.2, 10.0, 12.4, 15.5, 20.])
- >>> bins = np.array([0, 5, 10, 15, 20])
- >>> np.digitize(x,bins,right=True)
- array([1, 2, 3, 4, 4])
- >>> np.digitize(x,bins,right=False)
- array([1, 3, 3, 4, 5])
- """
- x = _nx.asarray(x)
- bins = _nx.asarray(bins)
- # here for compatibility, searchsorted below is happy to take this
- if np.issubdtype(x.dtype, _nx.complexfloating):
- raise TypeError("x may not be complex")
- mono = _monotonicity(bins)
- if mono == 0:
- raise ValueError("bins must be monotonically increasing or decreasing")
- # this is backwards because the arguments below are swapped
- side = 'left' if right else 'right'
- if mono == -1:
- # reverse the bins, and invert the results
- return len(bins) - _nx.searchsorted(bins[::-1], x, side=side)
- else:
- return _nx.searchsorted(bins, x, side=side)
|