{-# OPTIONS_HADDOCK hide #-}{-# LANGUAGE FlexibleInstances, FlexibleContexts, MultiParamTypeClasses #-}------------------------------------------------------------------------------- |-- Module : Data.Packed.Array.Internal-- Copyright : (c) Alberto Ruiz 2009-- License : GPL---- Maintainer : Alberto Ruiz <aruiz@um.es>-- Stability : provisional-- Portability : portable---- Multidimensional arrays.---- The arrays provided by this library are immutable, built on top of hmatrix-- structures.-- Operations work on complete structures (indexless), and dimensions have \"names\", -- in order to select the desired contractions in tensor computations.---- This module contains auxiliary functions not required by the end user.-----------------------------------------------------------------------------moduleNumeric.LinearAlgebra.Array.Internal(-- * Data structuresNArray,Idx(..),Name,order,namesR,names,size,sizesR,sizes,typeOf,dims,coords,Compat(..),-- * Array creationscalar,mkNArray,fromVector,fromMatrix,-- * Array manipulationrenameRaw,parts,partsRaw,(|*|),analyzeProduct,smartProduct,zipArray,mapArray,extract,onIndex,-- * UtilitiesseqIdx,reorder,sameStructure,conformable,makeConformant,mapTypes,mapNames,renameSuperRaw,renameExplicit,newIndex,basisOf,common,selDims,mapDims,takeDiagT,atT,firstIdx,fibers,matrixator,matrixatorFree,Coord,asMatrix,asVector,asScalar,debug)whereimportData.PackedimportData.ListimportNumeric.LinearAlgebra((<>),Field)importControl.ApplicativeimportData.Function(on)-- import Control.Parallel.StrategiesimportDebug.Tracedebugmfx=trace(m++show(fx))x-- | Types that can be elements of the multidimensional arrays.class(Num(Vectort),Fieldt)=>CoordtinstanceCoordDoubleinstanceCoord(ComplexDouble)-- | indices are denoted by strings, (frequently single-letter)typeName=String-- | Dimension descriptor.dataIdxi=Idx{iType::i,iDim::Int,iName::Name}deriving(Eq)instanceEqi=>Ord(Idxi)wherecompare=compare`on`iName-- | A multidimensional array with index type i and elements t.dataNArrayit=A{dims::[Idxi]-- ^ Get detailed dimension information about the array.,coords::Vectort-- ^ Get the coordinates of an array as a-- flattened structure (in the order specified by 'dims').}-- | development function not intended for the end usermkNArray::[Idxi]->Vectora->NArrayiamkNArray[]_=error"array with empty dimensions, use scalar"mkNArraydmsvec=Admsvwhereds=mapiDimdmsn=productdsv=ifdimvec==n&&minimumds>0thenvecelseerror$showds++" dimensions and "++show(dimvec)++" coordinates for mkNArray"-- | Create a 0-dimensional structure.scalar::Coordt=>t->NArrayitscalarx=A[](fromList[x])-- | Rename indices (in the internal order). Equal indices are contracted out.renameRaw::(Coordt,Compati)=>NArrayit->[Name]-- ^ new names->NArrayitrenameRawtns=contract(renameSuperRawtns)renameSuperRaw(Adv)l|lengthl==lengthd=Ad'v|otherwise=error$"renameRaw "++showd++" with "++showlwhered'=zipWithfdlfin=i{iName=n}mapDimsf(Adv)=A(mapfd)vmapTypes::(i1->i2)->NArrayi1t->NArrayi2tmapTypesf=mapDims(\i->i{iType=f(iTypei)})mapNames::(Name->Name)->NArrayit->NArrayitmapNamesf=mapDims(\i->i{iName=f(iNamei)})-- | Rename indices using an association list.renameExplicit::(Compati,Coordt)=>[(Name,Name)]->NArrayit->NArrayitrenameExplicital=g.mapNamesfwherefn=maybenid(lookupnal)gt=reorderorig(contractt)whereorig=nub(namesRt)\\common1t-- | Index names (in internal order).namesR::NArrayit->[Name]namesR=mapiName.dims-- | Index names.names::NArrayit->[Name]names=sort.namesR-- | Dimension of given index.size::Name->NArrayit->Intsizent=(iDim.head)(filter((n==).iName)(dimst))sizesR::NArrayit->[Int]sizesR=mapiDim.dimssizes::NArrayit->[Int]sizest=map(flipsizet)(namest)-- | Type of given index.typeOf::Compati=>Name->NArrayit->itypeOfnt=(iType.head)(filter((n==).iName)(dimst))-- | The number of dimensions of a multidimensional array.order::NArrayit->Intorder=length.dimsselDimsds=mapfwherefn=head$filter((n==).iName)ds----------------------------------------------------------common2t1t2=[n1|n1<-namesRt1,n2<-namesRt2,n1==n2]analyzeProduct::(Coordt,Compati)=>NArrayit->NArrayit->Maybe(NArrayit,Int)analyzeProductab=rwherenx=common2abdx1=selDims(dimsa)nxdx2=selDims(dimsb)nxok=and$zipWithcompatdx1dx2(tma,na)=matrixatorFreeanx(mb,nb)=matrixatorFreebnxmc=transtma<>mbda=selDims(dimsa)nadb=selDims(dimsb)nbdc=da++dbc=Adc(flattenmc)sz=product(mapiDimdc)r|ok=Just(c,sz)|otherwise=Nothinginfixl5|*|-- | Tensor product with automatic contraction of repeated indices, following Einstein summation convention.(|*|)::(Coordt,Compati)=>NArrayit->NArrayit->NArrayitt1|*|t2=caseanalyzeProductt1t2ofNothing->error$"wrong contraction2: "++(show$dimst1)++" and "++(show$dimst2)Just(r,_)->r----------------------------------------------------------lastIdxnamet=((d1,d2),m)where(d1,d2)=span(\d->iNamed/=name)(dimst)c=product(mapiDimd2)m=reshapec(coordst)firstIdxnamet=(nd,m')where((d1,d2),m)=lastIdxnametm'=reshapec$flatten$transmnd=d2++d1c=dim(coordst)`div`(iDim$headd2)-- | Obtain a matrix whose columns are the fibers of the array in the given dimension. The column order depends on the selected index (see 'matrixator').fibers::Coordt=>Name->NArrayit->Matrixtfibersn=snd.firstIdxn-- | Reshapes an array as a matrix with the desired dimensions as flattened rows and flattened columns.matrixator::(Coordt)=>NArrayit-- ^ input array->[Name]-- ^ row dimensions->[Name]-- ^ column dimensions->Matrixt-- ^ resultmatrixatortnrnc=reshapes(coordsq)whereq=reorder(nr++nc)ts=product(map(flipsizet)nc)-- | Reshapes an array as a matrix with the desired dimensions as flattened rows and flattened columns. We do not force the order of the columns.matrixatorFree::(Coordt)=>NArrayit-- ^ input array->[Name]-- ^ row dimensions->(Matrixt,[Name])-- ^ (result, column dimensions)matrixatorFreetnr=(reshapes(coordsq),nc)whereq=tridxnrtnc=drop(lengthnr)(mapiName(dimsq))s=product(map(flipsizet)nc)-- | Create a list of the substructures at the given level.parts::(Coordt)=>NArrayit->Name-- ^ index to expand->[NArrayit]partsaname|name`elem`(namesRa)=map(reorderorig)(partsRawaname)|otherwise=error$"parts: "++showname++" is not a dimension of "++(show$namesRa)whereorig=namesRa\\[name]partsRawaname=mapf(toRowsm)where(_:ds,m)=firstIdxnameaft=A{dims=ds,coords=t}tridx[]t=ttridx(name:rest)t=A(d:ds)(joints)whered=caselastIdxnametof((_,d':_),_)->d'_->error"wrong index sequence to reorder"ps=map(tridxrest)(partsRawtname)ts=mapcoordspsds=dims(headps)-- | Change the internal layout of coordinates.-- The array, considered as an abstract object, does not change.reorder::(Coordt)=>[Name]->NArrayit->NArrayitreordernsb|ns==namesRb=b|sortns==sort(namesRb)=tridxnsb|otherwise=error$"wrong index sequence "++showns++" to reorder "++(show$namesRb)------------------------------------------------------------------------ | Apply a function (defined on hmatrix 'Vector's) to all elements of a structure.-- Use @mapArray (mapVector f)@ for general functions.mapArray::Coordb=>(Vectora->Vectorb)->NArrayia->NArrayibmapArrayft|null(dimst)=scalar(f(coordst)@>0)|otherwise=mkNArray(dimst)(f(coordst))liftNA2f(Ad1v1)(A_d2v2)=Ad1(fv1v2)-- | Class of compatible indices for contractions.class(Eqa,Show(Idxa))=>Compatawherecompat::Idxa->Idxa->Boolopos::Idxa->Idxacontract1tname1name2|ok=foldl1'(liftNA2(+))y|otherwise=error$"wrong contraction1: "++(show$dimst)++" "++name1++" "++name2whereok=(compat<$>getNametname1<*>getNametname2)==JustTruex=map(flippartsRawname2)(partsRawtname1)y=maphead$zipWithdrop[0..]xgetNametname=dwherel=filter((==name).iName)(dimst)d=ifnulllthenNothingelseJust(headl)contract1ctn=contract1renamednn'wheren'=" "++n++" "-- forbid spaces in names...renamed=renameSuperRaw(t)auxnamesauxnames=h++(n':r)(h,_:r)=break(==n)(namesRt)common1t=[n1|(a,n1)<-x,(b,n2)<-x,a>b,n1==n2]wherex=zip[0::Int..](namesRt)contractt=foldl'contract1ct(common1t)--------------------------------------------------------------- | Check if two arrays have the same structure.sameStructure::(Eqi)=>NArrayit1->NArrayit2->BoolsameStructureab=sortBy(compare`on`iName)(dimsa)==sortBy(compare`on`iName)(dimsb)--------------------------------------------------------------- | Apply an element-by-element binary function to the coordinates of two arrays. The arguments are automatically made conformant.zipArray::(Coorda,Coordb,Compati)=>(Vectora->Vectorb->Vectorc)-- ^ transformation->NArrayia->NArrayib->NArrayiczipArrayoab=liftNA2oa'b'where(a',b')=makeConformantT(a,b)--------------------------------------------------------- | Create an array from a list of subarrays. (The inverse of 'parts'.)newIndex::(Coordt,Compati)=>i-- ^ index type->Name->[NArrayit]->NArrayitnewIndexinamets=rwhereds=Idxi(lengthts)name:(dims(headcts))cts=makeConformanttsr=mkNArrayds(join$mapcoordscts)--------------------------------------------------------- | Obtain a canonical base for the array.basisOf::Coordt=>NArrayit->[NArrayit]basisOft=map(dimst`mkNArray`)$toRows(ident.dim.coords$t)-------------------------------------------------------------instance(Coordt,Coord(Complext),Compati,ContainerVectort)=>Container(NArrayi)twheretoComplex(r,c)=zipArray(currytoComplex)rcfromComplext=let(r,c)=fromComplex(coordst)in(mapArray(constr)t,mapArray(constc)t)comp=mapArraycompconj=mapArrayconjreal=mapArrayrealcomplex=mapArraycomplex-- instance (NFData t, Element t) => NFData (NArray i t) where-- rnf = rnf . coords------------------------------------------------------------------------ | obtains the common value of a property of a listcommon::(Eqa)=>(b->a)->[b]->Maybeacommonf=commonval.mapfwherecommonval::(Eqa)=>[a]->Maybeacommonval[]=Nothingcommonval[a]=Justacommonval(a:b:xs)=ifa==bthencommonval(b:xs)elseNothing-------------------------------------------------------------------------- | Extract the 'Matrix' corresponding to a two-dimensional array,-- in the rows,cols order.asMatrix::(Coordt)=>NArrayit->MatrixtasMatrixa|ordera==2=reshapec(coordsa')|otherwise=error$"asMatrix requires a 2nd order array."wherec=size(last(namesRa'))a'a'=reorder(sort(namesRa))a-- | Extract the 'Vector' corresponding to a one-dimensional array.asVector::(Coordt)=>NArrayit->VectortasVectora|ordera==1=coordsa|otherwise=error$"asVector requires a 1st order array."-- | Extract the scalar element corresponding to a 0-dimensional array.asScalar::(Coordt)=>NArrayit->tasScalara|ordera==0=coordsa@>0|otherwise=error$"asScalar requires a 0th order array."-------------------------------------------------------------------------- | Create a 1st order array from a 'Vector'.fromVector::Compati=>i->Vectort->NArrayitfromVectoriv=mkNArray[Idxi(dimv)"1"]v-- | Create a 2nd order array from a 'Matrix'.fromMatrix::(Compati,Coordt)=>i->i->Matrixt->NArrayitfromMatrixiricm=mkNArray[Idxir(rowsm)"1",Idxic(colsm)"2"](flattenm)-------------------------------------------------------------------------- | Select some parts of an array, taking into account position and value.extract::(Compati,Coordt)=>(Int->NArrayit->Bool)->Name->NArrayit->NArrayitextractfnamearr=reorder(namesRarr).newIndex(typeOfnamearr)name.mapsnd.filter(uncurryf)$zip[1..](partsarrname)-- | Apply a list function to the parts of an array at a given index.onIndex::(Coorda,Coordb,Compati)=>([NArrayia]->[NArrayib])->Name->NArrayia->NArrayibonIndexfnamet=rwherer=ifsort(namesRx)==sort(namesRt)thenreorder(namesRt)xelsexx=newIndex(typeOfnamet)name(f(partstname))------------------------------------------------------------------------extendalldims(Adv)=reorder(allnames)swhereallnames=mapiNamealldimspref=alldims\\dn=product(mapiDimpref)s=A(pref++d)(join(replicatenv))-- | Obtains most general structure of a list of dimension specificationsconformable::Compati=>[[Idxi]]->Maybe[Idxi]conformableds|ok=Justalldims|otherwise=Nothingwherealldims=nub(concatds)allnames=mapiNamealldimsok=length(allnames)==length(nuballnames)-- | Converts a list of arrays to a common structure.makeConformant::(Coordt,Compati)=>[NArrayit]->[NArrayit]makeConformantts=caseconformable(mapdimsts)ofJustalldims->map(extendalldims)tsNothing->error$"makeConformant with inconsistent dimensions "++show(mapdimsts)-- the same version for tuples with possibly different element typesmakeConformantT(t1,t2)=caseconformable[dimst1,dimst2]ofJustalldims->(extendalldimst1,extendalldimst2)Nothing->error$"makeConformantT with inconsistent dimensions "++show(dimst1,dimst2)---------------------------------------------takeDiagT::(Compati,Coordt)=>NArrayit->[t]takeDiagTt=map(asScalar.atTt)cdswheren=minimum(sizesRt)o=ordertcds=map(replicateo)[0..n-1]atT::(Compati,Coordt)=>NArrayit->[Int]->NArrayitatTtc=atT'ctwhereatT'cs=foldl1'(.)(mapfpartcs)fpartkq=partsq(head(namesRq))!!k------------------------------------------------ not very smart...-- | This is equivalent to the regular 'product', but in the order that minimizes the size of the-- intermediate factors.smartProduct::(Coordt,Compati,Num(NArrayit))=>[NArrayit]->NArrayitsmartProduct[]=1smartProduct[a]=asmartProduct[a,b]=a*bsmartProductts=rwheren=lengthtsks=[0..n-1]xs=zipkstsgab=caseanalyzeProductabofNothing->error$"inconsistent dimensions in smartProduct: "++(show$dimsa)++" and "++(show$dimsb)Just(_,c)->cpairs=[((i,j),gab)|(i,a)<-initxs,(j,b)<-drop(i+1)xs](p,q)=fst$minimumBy(compare`on`snd)pairsr=smartProduct(ts!!p*ts!!q:(dropElemPosp.dropElemPosq)ts)dropElemPoskxs=takekxs++drop(k+1)xs------------------------------------------------ | sequence of n indices with given prefixseqIdx::Int->String->[Name]seqIdxnprefix=[prefix++showk|k<-[1..n]]