# Andriy Zhugayevych azh@ukr.net # BasicTools package # created 30.07.2005, modified - see version below # ################################################################################ #root: BasicTools #hfl: BasicTools #toc: _Overview BasicTools:=module() option package; global tex_pm, tex_cdot, tex_degree, tex_ll, tex_gg, tex_lnot, tex_mu, `index/symmetric12`, JP, `simplify/esu`,`simplify/emu`; export ModuleLoad, Setup, Round, Reduce2P, PeriodicMean, PeriodicStandardDeviation, angle2pair, mod1, coeffs2, convert2range, convert2repeat, simplifyradical, GetPath, # Miscellaneous tools SearchPos, Substitute, Dim2, Grid, Grid2, idL, iidL, idI, iidI, idS, iidS, idA, iidA, idSS, iidSS, nextpointS, isleq, Sort, Rort, SortM, SortIdx, MaxVal, MinVal, MaxIdx, MinIdx, SortEE, SortSign, IdentifyPairs, Classify2, BlockSplit, idExch, idSym, idASym, PrintTree, Tree2Seq, Tree2Seq2, Tree2List, # Array tools Squeeze, TrimLeft, TrimRight, Trim, StartNewLine, SearchKeySeq, ReadBody, ReadByMask, Split2, FormattedOperation, FormatFloat, AlignPoints, FormatTable, FormatTime2, PrintColumns, PrintMatrix, PrintUpperMatrix, PrintVector3d, ascii2unicode, # String tools CyclesAndTrees, defistree, CycleBasis2, # Combinatorics and algorithms Circumsphere, # Geometry fullparfrac2, fZero, fDblZero, fZeros, FZERO, NLPSolve2, HessianByGradient, PolynomialInterpolation1der, BezierControlPoints, BezierNodes, ExtrapolateMPE, # Calculus devpsolve, # Differential equations CheckEigenvector, NormalizeVector, NormalizeMatrix, RealJordanBlockMatrix, MatrixFunction2, SparseMM, SparseMHM, SparseMV, SparseMO, SparseMN, SparseMap, SuperimposeM, LocalizeVectors4, # Linear algebra RotationM, ReflectionM, RotationParam, RotationParamEulerX, RotationParamEulerY, RotationMEulerX, RotationMEulerY, EulerY2Rot, Rot2EulerY, RotationMp2q, RotationM2zx, TransSpherCoo, TransSpherCooI, TransTensor, ExtendedTransM, # Linear transformations Phim, NPhim2, Pnma, NPnma2, Ylm, Yxyz, YxyzM, YxyzR, YxyzL, Wignerd, WignerD, Wigner3j, Wigner6j, MatElemYlm, MatElemY, # Spherical harmonics LerchPhi2, polylogR, polylogR_init, # Other special functions ExpandPath, SimplifyPath, FoldersTree, LoadTextFile, UnloadTextFile, SearchFilePos, ReadLines, ReadValue, ReadPaginatedMatrix, WriteLines, ReplaceLines, ReadRecord, WriteRecord, ReadSet, ReadUTF8, WriteUTF8, ReadBIN, WriteBIN, PrintFile, convert2unix, convert2dos, OpenFolder, # File tools GetBox, Spline2, ErrorBar, TextOnLine, Rainbow, WritePlotData, plotps, CorrectEPS, ImageCoreBox, # Plotting and graphics RootMeanSquare, LinearFit2, LinearFitTrig, MultiDimFit, GaussianFit, FindExtremum, ReorderXYdata, Intersection, # Data manipulation SqueezeHTML, SqueezeTable, ReadHref, ResolveURI, GetXMLChild, ReadXML, WriteXML, # HTML/XML tools GetConstantSU, ConvertSU, # Physics tools BuildHelp, AssignModule, ReadModules, Check4Updates, TabulatedFunction, ReduceFloat, ReduceFloat2, ReduceFloat3, ReduceFloatGrid, Timing, InitTiming, ProcessSetupArgs; # System tools local version, MapleVersion, CyrCharMap, complexTrigBasis, conventionY, normalizedY, psplotoptions, phi,theta,psi,alpha,beta,gamma, polylogR_T, polylogR_n, polylogR_N, TextFiles,TextFilesAttr,TextFilesList,nTextFiles,sTextFiles, binformatidentifiers1, cc,hb,NA,kB,q2vu, s_TabulatedFunction, InitialTime,PrintTime,nDone,nLeft; ModuleLoad:=proc() local i,n,V; version:=20231227; # Check version MapleVersion:=sscanf(convert(kernelopts('version'),string),"Maple %f")[1]; if not(type(MapleVersion,numeric)) then error("Cannot determine Maple version: %1",kernelopts('version')) end; if (MapleVersion<10) then WARNING("Version 10 or later is needed for full operability. Detected version is %1",kernelopts('version')) end; # Special symbols tex_pm,tex_cdot,tex_degree,tex_ll,tex_gg,tex_lnot,tex_mu:=op(map(StringTools[Char],[177,183,176,171,187,172,181])); # CyrCharMap CyrCharMap:=table(): for i from 0 to 255 do CyrCharMap[i]:=cat("&#",i,";") end: for i from 192 to 255 do CyrCharMap[i]:=cat("&#",i+848,";") end: CyrCharMap[179]:="і": CyrCharMap[175]:="Ї": #³² CyrCharMap[191]:="ї": CyrCharMap[178]:="І": #¿¯ CyrCharMap[186]:="є": CyrCharMap[170]:="Є": #ºª CyrCharMap[184]:="ё": CyrCharMap[168]:="Ё": #¸¨ CyrCharMap[180]:="ґ": CyrCharMap[165]:="Ґ": #´¥ # Default settings Setup('InitTextFiles', 'complexTrigBasis'=true, 'conventionY'="standard", 'normalizedY'=true, 'nTextFiles'=1000, 'sTextFiles'=5*10^8, 'psplotoptions'="portrait,noborder,width=400,height=400"); # Indexing functions `index/symmetric12`:=proc(id::list(posint),A::Array,v::list) local i; if (nops(id)<2) then error("symmetric12 indexing requires at least 2 indices") end; i:=max(id[1],id[2]),min(id[1],id[2]),op(3..,id); if (nargs=2) then A[i] else A[i]:=op(v) end end; # Legendre and Jacobi polynomials JP:=orthopoly[P]; # LerchPhi patch polylogR_init(20,100); # Binary format identifiers binformatidentifiers1:=[1760568055,-148705176]; # System of units try Units[AddBaseUnit]('esu',context='SU',symbol='esu',dimension='electric_charge') catch: WARNING("Cannot add 'esu' unit because of the following error:"); WARNING(lasterror) end; try Units[AddBaseUnit]('emu',context='SU',symbol='emu',dimension='magnetic_flux') catch: WARNING("Cannot add 'emu' unit because of the following error:"); WARNING(lasterror) end; try Units[AddSystem]('SU',check=true,'cm'['SI'],'g','s','esu','emu','A','K','mol','cd') catch: WARNING("Cannot add 'SU' system because of the following error:"); WARNING(lasterror) end; Units[UseSystem]('SU',strict=true); cc:=ScientificConstants[GetValue](ScientificConstants[Constant]('c',system='SU')); hb:=ScientificConstants[GetValue](ScientificConstants[Constant]('hbar',system='SU')); NA:=ScientificConstants[GetValue](ScientificConstants[Constant]('N[A]',system='SU')); kB:=ScientificConstants[GetValue](ScientificConstants[Constant]('k',system='SU')); `simplify/esu`:=proc(expr) local u,v,e; v,u:=q2vu(simplify(expr)); e:=subs('A'=cc/10*'esu'/'s',parse(convert(op(u),string))); simplify(Units[Unit](v*e)) end; `simplify/emu`:=proc(expr) local u,v,e; v,u:=q2vu(simplify(expr)); e:=subs('A'=1/10*'emu'/'s',parse(convert(op(u),string))); simplify(Units[Unit](v*e)) end; NULL end: #hfl: Setup #toc: _Setup #BasicTools[Setup] Setup:=proc({InitTextFiles::boolean:=false,printout::boolean:=false}) local sq; sq:=ProcessSetupArgs([_rest],['complexTrigBasis','conventionY','normalizedY','polylogR_n','polylogR_N','nTextFiles','sTextFiles','psplotoptions'],['version']); if printout then printf("complexTrigBasis=%a, conventionY=%a, normalizedY=%a\nnTextFiles=%d, sTextFiles=%.0fMB\npsplotoptions=%a\n", complexTrigBasis,conventionY,normalizedY,nTextFiles,sTextFiles*1e-6,psplotoptions) end; if InitTextFiles then TextFiles:=table(); TextFilesAttr:=table(); TextFilesList:=[] end; sq end: ################################################################################ #cat: Miscellaneous tools #hfl: Round Round:=proc(v,n::integer,$) if type(v,realcons) then evalf(round(v*10^n)*10^(-n)) elif type(v,indexable) then map(Round,v,n) elif type(v,complex) then simplify(Round(Re(v),n)+Round(Im(v),n)*I) else error "Unrecognized format",v end end: #hfl: Reduce2P Reduce2P:=proc(x,P:=1,x0:=0) local e,i; e:=indets(x); x-P*floor(`if`(e={},(x-x0)/P,eval((x-x0)/P,{seq(e[i]=exp(-2.07-0.5*i),i=1..nops(e))}))) end: #hfl: Reduce2P PeriodicMean:=proc(V::indexable,i0::list:=[1],P:=1,x0:=0,$) local n,x1,x2,v; n:=add(1,v=V); x1:=V[op(i0)]-P/2; x1:=Reduce2P(add(Reduce2P(v,P,x1),v=V)/n,P,x0); x2:=x1-P/2; x2:=Reduce2P(add(Reduce2P(v,P,x2),v=V)/n,P,x0); if (abs(x1-x2)>n*10^(2-Digits)) then error("Differnt answers (%1<>%2), choose another i0",x1,x2) end; x2 end: #hfl: Reduce2P PeriodicStandardDeviation:=proc(V::indexable,i0::list:=[1],P:=1,x0:=0,$) local x,n,x1,v; x:=PeriodicMean(V,i0,P,x0); n:=add(1,v=V); x1:=x-P/2; evalf(sqrt(add((Reduce2P(v,P,x1)-x)^2,v=V)/(n-1))) end: #hfl: angle2pair angle2pair:=proc(a::numeric,n::integer,$) local a1,a0,da; if (n=0) then a,0 elif (n<0) then a0,da:=angle2pair(a+180/n,-n); Reduce2P(a0-180/n,360,-180+1e-8),da else a1:=Reduce2P(a,360,-180+`if`(type(n,odd),0,180/n)); a0:=round(a1/360*n)*360/n; a0,a1-a0 end end: #hfl: mod1 mod1:=proc(e::integer,m::posint,$) modp(e-1,m)+1 end: #hfl: coeffs2 coeffs2:=proc(e,vars::list,$) local e2,C,M,x,i; if (vars=[]) then return [[e]] end; e2:=collect(e,vars,'distributed'); if not(type(e2,polynom(anything,vars))) then error("Not a polynom of %1: %2",vars,e) end; C:=[coeffs(e2,vars,M)]; M:=[M]; [seq([seq(degree(M[i],x),x=vars),C[i]],i=1..nops(C))] end: #hfl: convert2range convert2range:=proc(ls0::list(integer),minnum::posint:=3,{asstring::boolean:=false},$) local ls,tb,i1,i,v; ls:=ls0; if (nops(ls)>1) then tb:=table(); i1:=1; for i from 2 to nops(ls) do if (ls[i]-ls[i-1]>1) then tb[i1]:=`if`(i=i1+1,ls[i1],[ls[i1],ls[i-1]]); i1:=i elif (ls[i]-ls[i-1]<1) then error("Nonmonotonic at i=%1 list %2",i,ls) end end; tb[i1]:=`if`(i=i1+1,ls[i1],[ls[i1],ls[i-1]]); ls:=map(v->`if`(type(v,list) and v[2]-v[1]+1`if`(type(v,list),sprintf("$%d..%d",op(v)),sprintf("%d",v)),ls))) else ls end end: #hfl: convert2repeat convert2repeat:=proc(ls::list,minnum::posint:=3,format::string:="",$) local tb,i,L,k,v; if (ls=[]) then L:=[] else tb:=table([1=1]); for i from 2 to nops(ls) do if (ls[i]<>ls[i-1]) then tb[i]:=1 end end; L:=sort([indices(tb,'nolist'),nops(ls)+1]); L:=[seq([ls[L[k]],L[k+1]-L[k]],k=1..nops(L)-1)]; L:=map(v->`if`(v[2]`if`(type(v,list),cat(sprintf(format,v[1]),"$",v[2]),sprintf(format,v)),L))) end end: #hfl: simplifyradical simplifyradical:=e->signum(e)*sqrt(expand(rationalize((radnormal(e)^2)))): #hfl: GetPath GetPath:=proc(X::{list,Vector},Y::{list,Vector},N::list(nonnegint),$) local d,o,XY,n,i; d:=nops(N); if not(foldl(`and`,seq(type(X[o]*N[o],integer),o=1..d))) then error("Non-grid X: %1, but N=%2",X,N) end; if not(foldl(`and`,seq(type(Y[o]*N[o],integer),o=1..d))) then error("Non-grid Y: %1, but N=%2",Y,N) end; XY:=Y-X; n:=igcd(seq(XY[o]*N[o],o=1..d)); [seq(X+i/n*XY,i=0..n)] end: ################################################################################ #cat: Array tools #hfl: SearchPos SearchPos:=proc(s::{string,list},xorf,$) local f,y,i; if (type(s,string) and not(type(xorf,procedure))) then return SearchText(xorf,s) end; if type(xorf,procedure) then f:=xorf else f:=y->y=xorf end; for i from 1 to `if`(type(s,string),length(s),nops(s)) do if f(s[i]) then return i end end; 0 end: #hfl: SearchPos Substitute:=proc(s::{string,list},xorf) local f,y,i; if (type(s,string) and not(type(xorf,procedure))) then return StringTools[Substitute](s,xorf,`if`(_nrest>0,_rest,"")) end; if type(xorf,procedure) then f:=xorf else f:=y->y=xorf end; for i from 1 to `if`(type(s,string),length(s),nops(s)) do if f(s[i]) then break end end; if (type(s,list) and i<=nops(s)) then subsop(i=_rest,s) elif (type(s,string) and i<=length(s)) then cat(s[..i-1],_rest,s[i+1..]) else s end end: #hfl: Dim2 Dim2:=(A::{rtable,list})->`if`(type(A,list),nops(A),`if`(type(A,Array),op(map2(op,2,[op(2,A)])),op(1,A))): #hfl: Grid Grid:=proc(rng::complexcons..complexcons,n::integer,f::procedure:=(x->x),avoid::list:=[],$) local x1,x2,x,m,k,t,t1,t2,ls,lsex; x1,x2:=op(rng); m:=`if`(n<2,1,n-1); lsex:=avoid; if (f(x)=x) then ls:=[seq(x1+(x2-x1)*k/m,k=0..m)] else if (type(x1,float) or type(x2,float)) then lsex:=evalf(lsex); t1,t2:=max(fsolve(f(t)=x1,t)),max(fsolve(f(t)=x2,t)); else t1,t2:=max(solve(f(t)=x1,t)),max(solve(f(t)=x2,t)) end; ls:=[seq(f(t1+(t2-t1)*k/m),k=0..m)] end; remove(member,ls,lsex) end: #hfl: Grid Grid2:=proc(lsp::listlist,dx::numeric,c::numeric:=1.1,$) local e,n,ls,v,i; if (nops(lsp)=1) then [map(v->round(v/dx)*dx,lsp[1])] elif (nops(lsp)=2) then n:=round(c*sqrt(add(v^2,v=lsp[2]-lsp[1]))/dx); if (n=0) then Grid2([(lsp[1]+lsp[2])/2],dx,c) else e:=(lsp[2]-lsp[1])/n; ListTools[MakeUnique]([seq(map(v->round(v/dx)*dx,lsp[1]+e*i),i=0..n)]) end else ls:=[seq(Grid2([lsp[i-1],lsp[i]],dx,c),i=2..nops(lsp))]; [seq(op(..-2,v),v=ls[..-2]),op(ls[-1])] end end: #hfl: idL idL:=proc(x::list(integer),L::{posint,list(posint)},$) local d,LL,n,i; d:=nops(x); LL:=`if`(type(L,posint),[L$d],L); n:=x[d] mod LL[d]; for i from d-1 by -1 to 1 do n:=n*LL[i]+(x[i] mod LL[i]) end; n+1 end: #hfl: idL iidL:=proc(n::posint,L::list(posint),$) local d,x,p,i; d:=nops(L); x:=array(1..d); p:=n-1; for i from 1 to d do x[i]:=irem(p,L[i],'p') end; convert(x,list) end: #hfl: idL idI:=m->2*abs(m)+`if`(m>0,0,1): #hfl: idL iidI:=n->`if`(type(n,odd),-1,1)*iquo(n,2): #hfl: idL idS:=proc(x::list(posint),$) local i,j; i:=max(x[1],x[2]); j:=min(x[1],x[2]); i*(i-1)/2+j end: #hfl: idL iidS:=proc(n::posint,$) local i; i:=round(sqrt(2*n)); [i,n-i*(i-1)/2] end: #hfl: idL idA:=proc(x::list(posint),$) local i,j; if (x[1]=x[2]) then 0 else i:=max(x[1],x[2]); j:=min(x[1],x[2]); (i-1)*(i-2)/2+j end end: #hfl: idL iidA:=proc(n::posint,$) local i; i:=round(sqrt(2*n))+1; [i,n-(i-1)*(i-2)/2] end: #hfl: idL idSS:=proc(x::list(posint),$) local n,m; n:=idS(x[1..2]); m:=idS(x[3..4]); [max(n,m),min(m,n)] end: #hfl: idL iidSS:=nm->[op(iidS(nm[1])),op(iidS(nm[2]))]: #hfl: idL nextpointS:=proc(x::list(nonnegint),$) local d,i; d:=nops(x); for i from d by -1 to 2 while (x[i-1]=x[i]) do end; [op(..i-1,x),x[i]+1,0$(d-i)] end: #hfl: isleq isleq:=proc(u,v) local n,i; if type(u,list) then n:=min(nops(u),nops(v)); for i from 1 to n while (u[i]=v[i]) do end; if (i>n) then true else isleq(u[i],v[i]) end else evalb(uu),id::list:=[],$)::list; local F2; if (id=[]) then F2:=`if`(type(F,procedure(anything,anything)),F,(u,v)->evalb(F(u)F(u[op(id)],v[op(id)]),(u,v)->evalb(F(u[op(id)])u),id::list:=[],$)::list; local F2; if (id=[]) then F2:=`if`(type(F,procedure(anything,anything)),(u,v)->F(v,u),(u,v)->evalb(F(v)F(v[op(id)],u[op(id)]),(u,v)->evalb(F(v[op(id)])u),$)::list; #local F; #if type(ForId,list) then F:=(u,v)->evalb(u[op(ForId)]evalb(ForId(u)sort(ls,proc(u,v) local i; for i in id while (u[i]=v[i]) do end; isleq(u[i],v[i]) end): #hfl: SortIdx SortIdx:=proc(ls::indexable, ForId:=(u->u), {nolist::boolean:=false}, $)::list; local tb,i,id,F,ls2,v; if type(ls,list) then tb:=ls; id:=[[i]$i=1..nops(tb)] else tb:=convert(ls,table); id:=[indices(tb)] end; if type(ForId,list) then F:=(u,v)->evalb(u[1][op(ForId)]ForId(u[1],v[1]) elif type(ForId,procedure) then F:=(u,v)->evalb(ForId(u[1])u),$)::list; local V,x,y; V:=`if`(type(ls,list),ls,[entries(convert(ls,table),'nolist')]); if (V=[]) then return [] end; x:=V[1]; if type(ForId,list) then for y in V do if (x[op(ForId)]u),$)::list; local V,x,y; V:=`if`(type(ls,list),ls,[entries(convert(ls,table),'nolist')]); if (V=[]) then return [] end; x:=V[1]; if type(ForId,list) then for y in V do if (y[op(ForId)]u),$)::list; local tb,id,i,j,x; if type(ls,list) then tb:=ls; id:=[[i]$i=1..nops(tb)] else tb:=convert(ls,table); id:=[indices(tb)] end; if (id=[]) then return [] end; i:=id[1]; x:=tb[op(i)]; if type(ForId,list) then for j in id do if (x[op(ForId)]u),$)::list; local tb,id,i,j,x; if type(ls,list) then tb:=ls; id:=[[i]$i=1..nops(tb)] else tb:=convert(ls,table); id:=[indices(tb)] end; if (id=[]) then return [] end; i:=id[1]; x:=tb[op(i)]; if type(ForId,list) then for j in id do if (tb[op(j)][op(ForId)]>, F::procedure:=(u->u),$) local P,v; P:=SortIdx([seq(v,v=ev)],F,'nolist'); ev[P],`if`(evc=<<0>>,NULL,evc[..,P]) end: #hfl: SortIdx SortSign:=proc(ls::{list,Vector},$) local V,v,i,j,s; V:=convert(ls,Vector); s:=0; for i from 2 to op(1,V) do v:=V[i]; for j from i-1 by -1 to 1 while (V[j]>v) do V[j+1]:=V[j] end; s:=s+i-1-j; V[j+1]:=v end; `if`(type(ls,list),convert(V,list),V), `if`(type(s,odd),-1,1) end: #hfl: IdentifyPairs IdentifyPairs:=proc(ls1::indexable, ls2::indexable, idf::procedure, output::string:="l", {nolist::boolean:=false}, $) local id1,id2,idv,i1,i2,n,V,i,k1,k2,R,c; id1,id2:=`if`(type(ls1,list),{seq([i],i=1..nops(ls1))},{indices(ls1)}),`if`(type(ls2,list),{seq([i],i=1..nops(ls2))},{indices(ls2)}); idv:=table(); for i1 in id1 do for i2 in id2 do idv[i1,i2]:=idf(ls1[op(i1)],ls2[op(i2)]) end end; n:=min(nops(id1),nops(id2)); V:=Vector(n); for i from 1 to n do k1,k2:=op(MinIdx(idv)); # V[i]:=[`if`(nolist,op(k1),k1),`if`(nolist,op(k2),k2),idv[k1,k2]]; V[i]:=[k1,k2,idv[k1,k2]]; for i2 in id2 do idv[k1,i2]:=evaln(idv[k1,i2]) end; id1:=`minus`(id1,{k1}); for i1 in id1 do idv[i1,k2]:=evaln(idv[i1,k2]) end; id2:=`minus`(id2,{k2}) end; R:=[op(id1),op(id2)]; seq(`if`(c="l", [seq(`if`(nolist,[op(V[i][1]),op(V[i][2]),V[i][3]],V[i]),i=1..n)], `if`(c="v", [seq([ls1[op(V[i][1])],ls2[op(V[i][2])],V[i][3]],i=1..n)], `if`(c="r", `if`(nolist,map(op,R),R), `if`(c="m", `if`(n=0,0,V[n][3]), NULL)))),c=output) end: #hfl: Classify2 Classify2:=proc(ls,ForId::{procedure,list},{index::boolean:=false,nolist::boolean:=false},$) local F,tb,nolist2,lsi,F2,i,id,v; if type(ForId,list) then F:=v->v[op(ForId)] else F:=v->ForId(v) end; tb:=table(); if index then if type(ls,list) then nolist2:=true; lsi:=[$1..nops(ls)] elif type(ls,table) then nolist2:=nolist; lsi:=[indices(ls,`if`(nolist2,':-nolist',NULL))] elif type(ls,rtable) then nolist2:=evalb(rtable_num_dims(ls)=1); lsi:=`if`(nolist2,map(v->lhs(v),rtable_elems(ls)),map(v->[lhs(v)],rtable_elems(ls))) else error("Indexing is supported only for list,table,rtable") end; if nolist2 then F2:=i->F(ls[i]) else F2:=i->F(ls[op(i)]) end; for i in lsi do id:=F2(i); if assigned(tb[id]) then tb[id]:=[op(tb[id]),i] else tb[id]:=[i] end end; op(tb) elif type(ls,list) then for v in ls do id:=F(v); if assigned(tb[id]) then tb[id]:=[op(tb[id]),v] else tb[id]:=[v] end end; op(tb) else ListTools[Classify](F,ls) end end: #hfl: BlockSplit BlockSplit:=proc(L::list,B::list(nonnegint),$) local N,n,A,i,R,v; N,n:=nops(L),nops(B); A:=Vector(B); if (A[1]>=N) then [L] else for i from 2 to n do A[i]:=A[i-1]+A[i]; if (A[i]>N) then A[i]:=N; n:=i; break end end; R:=[[1,min(N,A[1])],seq([A[i-1]+1,A[i]],i=2..n)]; [seq(L[v[1]..v[2]],v=R)] end end: #hfl: idExch idExch:=proc(e,A::name,id1,id2,$) local s,pre,L,i,j,k,v,x; s:=convert(e,string); pre:=cat(convert(A,string),"["); L:=length(pre); i:=1; do k:=searchtext(pre,s,i..-1); if (k=0) then break end; i:=i+k+L-2; j:=searchtext("]",s,i..-1)+i-1; if (j=i-1) then error "There is no right bracket after left bracket at position %1",i end; v:=parse(s[i..j]); v:=convert(subs(id1=x,id2=id1,x=id2,v),string); s:=cat(s[1..i-1],v,s[j+1..-1]) end; parse(s) end: #hfl: idExch idSym:=proc(e,A::name,i,j) local e2,k; if (nargs=4) then (e+idExch(e,A,i,j))/2 elif (nargs=5) then k:=args[5]; e2:=idExch(e,A,i,j): (e+idExch(e,A,i,k)+idExch(e,A,j,k)+e2+idExch(e2,A,i,k)+idExch(e2,A,j,k))/6 else error "No more than three indices are supported" end end: #hfl: idExch idASym:=proc(e,A::name,i,j,$) (e-idExch(e,A,i,j))/2 end: #hfl: PrintTree PrintTree:=proc( tree, level::nonnegint:=0, str::procedure:=(x->x), { filedesc::{string,integer}:="", prefix::string:="." },$) local s,tree1; s:=cat(prefix$level,str(tree[1]),"\n"); if (filedesc="") then printf(s) else fprintf(filedesc,s) end; for tree1 in tree[2..-1] do PrintTree(tree1,level+1,str,':-filedesc'=filedesc,':-prefix'=prefix) end; NULL end: #hfl: Tree2Seq Tree2Seq:=proc( tree::list, prefix::string:="", str::procedure:=(x->x), {separator::string:="/"},$) local tree1; cat(prefix,str(tree[1])),seq(Tree2Seq(tree1,cat(prefix,str(tree[1]),separator),str,':-separator'=separator),tree1=tree[2..]) end: #hfl: Tree2Seq Tree2Seq2:=proc( tree::list, index::list(posint):=[], size::posint:=1, body::procedure:=(x->x),$) local i; [body(`if`(size=1,tree[1],tree[..size])),index,nops(tree)-size],seq(Tree2Seq2(tree[size+i],[op(index),i],size,body),i=1..nops(tree)-size) end: #hfl: Tree2Seq Tree2List:=proc( tree::list, size::posint:=1, body::procedure:=(x->x),$) local ls,id,i; ls:=[Tree2Seq2(tree,size,body)]; id:=table([seq(op(ls[i][2])=i,i=1..nops(ls))]); map(v->[v[1],`if`(v[2]=[],0,id[op(..-2,v[2])]),[seq(id[op(v[2]),i],i=1..v[3])]],ls) end: ################################################################################ #cat: String tools #hfl: specsymbols #specsymbols #hfl: Trim Squeeze:=proc(s::string,key::string:=" ", multiplicity::posint:=2, $) local dkey,s1,s2,i,j; dkey:=cat(key$multiplicity); s1:=s: s2:="": i:=searchtext(dkey,s1); while (i>0) do s2:=cat(s2,s1[1..i]); j:=i: while (s1[j+1]=key) do j:=j+1 end: s1:=s1[j+1..-1]; i:=searchtext(dkey,s1) end: cat(s2,s1); end: #hfl: Trim TrimLeft:=proc(s::string, ls::list:=[" "]) local i; for i from 1 to length(s) do if not(s[i] in ls) then break end end; s[i..-1] end: #hfl: Trim TrimRight:=proc(s::string, ls::list:=[" "]) local i; for i from length(s) by -1 to 1 do if not(s[i] in ls) then break end end; s[1..i] end: #hfl: Trim Trim:=proc(s::string, ls::list:=[" "]) TrimRight(TrimLeft(s,ls),ls) end: #hfl: StartNewLine StartNewLine:=proc(s::string,ls::list:=[],$) local s1,s2; s1:=s; for s2 in ls do s1:=StringTools[SubstituteAll](s1,s2,cat("\n",s2)) end; s1 end: #hfl: SearchKeySeq SearchKeySeq:=proc(s::string, keyseq::list(string), { start::posint:=1, skip::nonnegint:=0, sep::realcons:=infinity, fulloutput::boolean:=false, noerror::boolean:=false },$) local i,n,j,lsp,k; i:=start; for n from 0 to skip do j:=searchtext(keyseq[1],s,i..-1); if (j=0) then if noerror then return `if`(fulloutput,[],0) else printf(s); error("No key%1 \"%2\" is found after %3/%4 skips",1,keyseq[1],n,skip) end end; i:=i+j-1+length(keyseq[1]) end; lsp:=array(1..nops(keyseq)): lsp[1]:=i; for k from 2 to nops(keyseq) do j:=searchtext(keyseq[k],s,i..-1); if (j=0) then if noerror then return `if`(fulloutput,[],0) else printf(s); error("No key%1 \"%2\" is found",k,keyseq[k]) end end; if (j>sep+1) then if noerror then return `if`(fulloutput,[],0) else printf(s); error("key%1 \"%2\" is separated by %3>%4 positions",k,keyseq[k],j-1,sep) end end; i:=i+j-1+length(keyseq[k]); lsp[k]:=i end; `if`(fulloutput,convert(lsp,list),i) end: #hfl: ReadBody ReadBody:=proc( s::string, key1::string, key2::string:="", trim::list:=[], { noerror::boolean:=false, restoftheline::boolean:=false, format::string:="", wkeys::boolean:=false, pos1::name:=undefined, pos2::name:=undefined},$) local i,j,s2,v; i:=SearchText(key1,s); if (i=0) then if noerror then return "" else error( "No key1 (\"%1\") in %2",key1,s) end else if not(type(pos1,undefined)) then assign(pos1,i-1) end; i:=i+length(key1) end; j:=`if`(key2="",length(s),SearchText(key2,s,i..-1)+i-2); if (j=i-2 and not(key2="")) then if noerror then if restoftheline then j:=length(s) else return "" end else error("No key2 (\"%1\") in %2",key2,s) end end; if not(type(pos2,undefined)) then assign(pos2,j+length(key2)+1) end; s2:=`if`(trim=[],s[i..j],Trim(s[i..j],trim)); if (format="") then `if`(wkeys,cat(key1,s2,key2),s2) else v:=sscanf(s2,format); `if`(v=0 or v=[],s2,op(v)) end end: #hfl: ReadBody ReadByMask:=proc(mask::string,s::string,star::character:="*",$) local i,i1,s1; if (mask[1]="*") then if (mask="*") then s else i:=SearchText(star,mask[2..]); if (i=1) then error("Consecutive stars are not allowed in mask: %1",mask) end; s1:=mask[2..`if`(i=0,-1,i)]; i1:=SearchText(s1,s); if (i1=0) then error("No suffix from mask (%1) in search string (%2)",mask,s) end; s[..i1-1],`if`(i=0,NULL,ReadByMask(mask[i+1..],s[i1+length(s1)..],star)) end else i:=SearchText(star,mask); if (i=0) then if (mask=s) then NULL else error("No match: %1<>%2",mask,s) end else if (mask[..i-1]=s[..i-1]) then ReadByMask(mask[i..],s[i..],star) else error("No prefix from mask (%1) in search string (%2)",mask,s) end end end end: #hfl: Split2 Split2:=proc(s::string,c::string:=" ",{output::string:="s"},$) local i1,n,lsi,i,lss,k,o; if (length(c)<>1) then error("Unrecognized split character: \"%1\"",c) end; i1:=0; n:=length(s); lsi:=table([1=1,n+1=n+1]); do for i1 from i1+1 to n while (s[i1]=c) do end; if (i1>n) then break end; i:=SearchText(c,s,i1..n); if (i>0) then i1:=i1+i-1; lsi[i1]:=i1 else break end end; lsi:=convert(lsi,list); lss:=[seq(s[lsi[k-1]..lsi[k]-1],k=2..nops(lsi))]; seq(`if`(o="s",lss,`if`(o="i",lsi,NULL)),o=output) end: #hfl: FormattedOperation FormattedOperation:=proc(s::string,f::procedure,pos::integer,iorn::integer,$) local i2,i1,s1,ls,S1,n,V1,V2,S2,k,i,j,d,s2,v; i2:=length(s); if (iorn>0) then i2:=min(i2,iorn) end; if (pos>0) then i1,i2,s1:=pos,i2,s[pos..i2] else i1,i2,s1:=1,i2,cat(" "$(-pos),s[..i2]) end; S1:=Split2(s1); if (nops(S1)>1 and Trim(S1[-1])="") then S1:=[op(..-3,S1),cat(S1[-2],S1[-1])] end; if (iorn>0) then n:=nops(S1) elif (iorn=0) then return cat(s1,s[i2+1..]) else n:=-iorn; if (n>nops(S1)) then error("Expected %1 records, detected %2 in %3",n,nops(S1),s) end; S1:=S1[..n]; i2:=i1+add(length(v),v=S1)-1+min(0,pos) end; V1:=Vector(map(op@sscanf,S1,"%f")); V2:=f(V1); S2:=Vector(n); for k from 1 to n do i:=SearchText(".",S1[k]); for j from i+1 to length(S1[k]) while StringTools[IsDigit](S1[k][j]) do end; d:=j-i-1; s2:=sprintf("%.*f%s",d,V2[k],S1[k][j..]); j:=SearchText(".",s2); if (j>i) then WARNING("String %1 is expanded by %2 symbols",S1[k],j-i) end; S2[k]:=cat(" "$max(0,i-j),s2) end; cat(s[..i1-1],seq(v,v=S2),s[i2+1..]) end: #hfl: FormatFloat FormatFloat:=proc( v::realcons, d1::{nonnegint,realcons}, d2::posint:=`if`(type(d1,nonnegint),max(1,d1),9), { maxorder::nonnegint:=3, asfloat::boolean:=false, width::nonnegint:=0, alignsign::boolean:=false, showpoint::boolean:=false, showzero::boolean:=false, hidezero::boolean:=false, showsign::boolean:=false, hidenull::boolean:=false },$) local a,rmin,ppos,r,s; if type(d1,nonnegint) then if (not(asfloat) and type(v,integer)) then r:=sprintf("%d",abs(v)); if showpoint then r:=cat(r,"."); ppos:=0 else ppos:=-1 end elif (not(asfloat) and type(v,rational) and type(v*10^(d1+maxorder),integer)) then for ppos from d1 to d1+maxorder do if type(v*10^ppos,integer) then break end end; r:=FormatFloat(evalf(abs(v)),ppos,d2,':-maxorder'=0,':-showzero'=showzero,':-hidezero'=true) elif (not(asfloat) and type(v,float) and -op(2,v)<=d1+maxorder) then r:=sprintf("%.*f",-op(2,v),abs(v)); if (not(showzero) and r[..2]="0.") then r:=r[2..] end; ppos:=max(0,-op(2,v)) else a:=abs(evalf(v)); rmin:=10^(d2-1); for ppos from d1 to d1+maxorder do r:=round(a*10^ppos); if (r>=rmin) then break end end; if (ppos>d1+maxorder) then ppos:=d1+maxorder end; r:=sprintf("%0*d",ppos,r); if (length(r)>ppos) then if (ppos=0 and d1=0) then if showpoint then r:=cat(r,".") else if (width>0) then r:=cat(r," ") end end else r:=cat(r[..-ppos-1],".",r[-ppos..]) end else r:=cat(`if`(showzero,"0.","."),r) end end else for ppos from 0 to infinity while (round(d1*10^ppos)<=d2) do end; ppos:=ppos-1; r:=sprintf("%.*f(%.0f)",ppos,abs(v),d1*10^ppos); if (not(showzero) and r[..2]="0.") then r:=r[2..] end end; if (evalf(v)<0) then r:=cat("-",r) elif showsign then r:=cat("+",r) end; r:=cat(r," "$min(d1+maxorder-ppos,width-length(r))); s:=cat(""," "$(width-length(r))); s:=`if`(alignsign and (r[1]="+" or r[1]="-"),cat(r[1],s,r[2..]),cat(s,r)); if (hidenull and parse(s)=0) then s:=`if`(width=0,"",cat(" "$width)) end; `if`(hidezero and ppos>0,TrimRight(s,["0"]),s) end: #hfl: FormatFloat AlignPoints:=proc(ls::list(string),$) local n,V,i,s,p,l1,l2,f,v; n:=nops(ls); V:=Vector(n): for i from 1 to n do s:=ls[i]; p:=searchtext(".",s); V[i]:=`if`(p=0,[s," ",""],[s[..p-1],".",s[p+1..]]) end; l1,l2:=max(seq(length(V[i][1]),i=1..n)),max(seq(length(V[i][3]),i=1..n)); f:=cat("%",l1,"s%s%-",l2,"s"); [seq(sprintf(f,op(v)),v=V)] end: #hfl: FormatTable FormatTable:=proc( tb::table, format::string:=" %a", { prefix::string:="", align::boolean:=false, printout::boolean:=false },$) local ls,N,n,S,p,V,c,dp,i,islist,s,v; ls:=[indices(tb)]; N:=nops(ls); n:=min(seq(nops(v),v=ls)); ls:=SortM(ls,[$1..n]); S:=Vector(N,i->convert(ls[i],string)[2..-2]); if align then for i from 1 to N do S[i]:=cat(S[i],",") end; p:=1; do V:=[seq(`if`(p>length(s),0,SearchText(",",s,p..-1)),s=S)]; dp:=max(op(V)); if (dp=0) then break end; for i from 1 to N do if (V[i]"%a"); ls:=[seq(sprintf(cat("( %s )=",format),S[i],`if`(islist,op(tb[op(ls[i])]),tb[op(ls[i])])),i=1..N)]; if printout then for s in ls do printf("%s\n",s) end end; ls end: #hfl: FormatTime2 FormatTime2:=proc(t::numeric,{nodays::boolean:=false},$) local s,m,h,d; s:=irem(round(t),60,'m'); if (m=0) then return sprintf("%d",s) end; m:=irem(m,60,'h'); if (h=0) then return sprintf("%d:%02d",m,s) end; if (nodays or h<24) then return sprintf("%d:%02d:%02d",h,m,s) end; h:=irem(h,24,'d'); sprintf("%d:%02d:%02d:%02d",d,h,m,s) end: #hfl: PrintColumns PrintColumns:=proc(ls::list(string),n0::posint,sep::string:=" ",{before::string:=""},$) local m,n,M,L,i,j; m:=ceil(nops(ls)/n0); n:=ceil(nops(ls)/m); M:=Matrix(m,n,(i,j)->`if`(i+(j-1)*m>nops(ls),"",ls[i+(j-1)*m])); L:=Vector(n,j->max(seq(length(M[i,j]),i=1..m))); for i from 1 to m do printf(before); for j from 1 to n-1 do printf("%-*s%s",L[j],M[i,j],sep) end; printf("%-*s\n",L[n],M[i,n]) end; end: #hfl: PrintMatrix PrintMatrix:=proc( M::Matrix, fmt::{integer,string}:=-interface(displayprecision), corner::string:="", { maxsize::posint:=99, head::{boolean,[indexable,indexable]}:=false, left::boolean:=false, Ichar::character:="i", donotprint::boolean:=false },$) local n,m,d,M1,M2,n1,d1,zero,print0,l,W,C1,R1,w1; n,m:=Dim2(M); if (max(n,m)>maxsize) then error("Max size of %1 is exceeded: [%2,%3]",maxsize,n,m) end; if type(fmt,negint) then d:=-fmt; M1:=map(ReduceFloat2,Matrix(map(Round,M,d),datatype=anything)); n1,d1:=interface(rtablesize,displayprecision); if not(donotprint) then interface(rtablesize=max(n,m),displayprecision=d); print(M1); interface(rtablesize=n1,displayprecision=d1) end else if type(fmt,string) then d:=3; M1:=map2(sprintf,fmt,M) else d:=fmt; zero:=cat("."," "$d); print0:=s->`if`(member(TrimRight(Trim(s),["0"]),["0.","-0."]),zero,s); if (MatrixOptions(M,datatype)=complex[8]) then M1,M2:=seq(map(v->print0(sprintf("%.*f",d,v)),f(M)),f=[Re,Im]); M1:=Matrix(n,m,(i,j)->`if`(M2[i,j]=zero,cat(M1[i,j]," "),`if`(M1[i,j]=zero,cat(M2[i,j],Ichar),cat(M1[i,j],`if`(M2[i,j][1]="-","","+"),M2[i,j],Ichar)))) else M1:=map(v->print0(sprintf("%.*f",d,v)),M) end end; W:=Vector(m,j->max(map(length,M1[..,j])),datatype=integer); if (head=false) then M1:=Matrix(n,m,(i,j)->StringTools[PadLeft](M1[i,j],W[j])) else if type(head,list) then C1,R1:=op(head) else C1:=Vector(n,i->sprintf("%3d",i)); R1:=Vector(m,j->sprintf( "%d",j)) end; C1:=Vector(n+1,i->`if`(i=1,corner,C1[i-1])); w1:=max(map(length,C1)); W:=Vector(m,j->max(W[j],length(R1[j])),datatype=integer); M1:=Matrix(n+1,m+1,(i,j)->`if`(i=1, `if`(j=1,StringTools[PadRight](C1[1],w1),StringTools[PadLeft](StringTools[PadRight](StringTools[Center](R1[j-1],3),d),W[j-1])), `if`(j=1,`if`(left,StringTools[PadRight],StringTools[PadLeft])(C1[i],w1),StringTools[PadLeft](M1[i-1,j-1],W[j-1])))) end; if not(donotprint) then printf("%s\n",M1) end end; M1 end: #hfl: PrintMatrix PrintUpperMatrix:=proc(M::Matrix,fmt::string,width::posint:=16,height::posint:=99,{step::posint:=1},$) local n,m,i1,i; n,m:=op(1,M); for i1 from 1 to step do for i from i1 by step to min(height,n,m) do printf(cat(fmt,"\n"),Vector(min(width,m-i+1),j->M[i,i+j-1])) end end end: #hfl: PrintVector3d PrintVector3d:=proc(V::indexable,$) local U,i1,i2,i3,k,v; U:=[seq(abs(v),v=V)]; i1,i2,i3:=op(SortIdx(U,'nolist')); k:=op(MaxIdx([U[i3],0.70710678*(U[i3]+U[i2]),0.57735027*(U[i3]+U[i2]+U[i1])])); if (k=1) then `if`(V[i3]>0,"XYZ","xyz")[i3] elif (k=2) then sort(cat(`if`(V[i3]>0,"XYZ","xyz")[i3],`if`(V[i2]>0,"XYZ","xyz")[i2])) else cat(`if`(V[1]>0,"X","x"),`if`(V[2]>0,"Y","y"),`if`(V[3]>0,"Z","z")) end end: #hfl: ascii2unicode ascii2unicode:=proc(s::string,CharMap::table:=CyrCharMap,imax::posint:=127,$) local ls,i; ls:=convert(s,'bytes'); for i in ls do if (i>imax) then return cat(seq(CharMap[i],i=ls)) end end; s end: ################################################################################ #cat: Combinatorics and algorithms #hfl: CyclesAndTrees CyclesAndTrees:=proc( G::Graph, istree::procedure:=defistree, extracycles::list(set):=[], {output::{1,2}:=1},$) local nn,cycles,nc,nc1,trees,nt,cbonds,tbonds,p,q,i,j,tb,b; nn:=GraphTheory[Neighbors](G); cycles:=map(convert,GraphTheory[CycleBasis](G),set); cycles:=[op(cycles),op(extracycles)]; nc:=nops(cycles); while (nc<>nc1) do cycles:=map(`union`@op,[ListTools[Categorize]((u,v)->`intersect`(u,v)<>{},cycles)]); nc1,nc:=nc,nops(cycles) end; cycles:=map(convert,cycles,list); trees:=GraphTheory[ConnectedComponents](GraphTheory[DeleteVertex](G,map(op,cycles))); nt:=nops(trees); cycles:=Vector(nc,p->cycles[p]); cbonds,tbonds:=Vector(nc,p->[]),Vector(nt,p->[]): for p from 2 to nc do for i in cycles[p] do for j in nn[i] do if not(member(j,cycles[p])) then for q from 1 to p-1 do if member(j,cycles[q]) then cbonds[p]:=[op(cbonds[p]),[i,j,q]]; cbonds[q]:=[op(cbonds[q]),[j,i,p]]; break end end end end end end; for p from 1 to nt do for i in trees[p] do for j in nn[i] do if not(member(j,trees[p] )) then for q from 1 to nc do if member(j,cycles[q]) then tbonds[p]:=[op(tbonds[p]),[i,j,q]]; break end end end end end end; tb:=table(): for p from 1 to nt do if (nops(tbonds[p])=1 and not(istree(trees[p],op(1..2,tbonds[p][1]),cycles[tbonds[p][1][3]]))) then cycles[tbonds[p][1][3]]:=[op(cycles[tbonds[p][1][3]]),op(trees[p])] else tb[p]:=[trees[p],op(tbonds[p])] end end; trees:=convert(tb,list); nt:=nops(trees); for p from 1 to nt do for b in trees[p][2..] do cbonds[b[3]]:=[op(cbonds[b[3]]),[b[2],b[1],nc+p]] end end; [seq([sort(cycles[p]),op(cbonds[p])],p=1..nc),op(trees)],`if`(output=1,NULL,nc) end: defistree:=((tree,i,j,cycle)->evalb(nops(tree)>1)): #hfl: CycleBasis2 CycleBasis2:=proc(G::Graph,{maxiter::posint:=9},$) local ff,C,B,N,T,iter,p,q,v; ff:=f->select(v->nops(select(member,C[v],f))>2,f); C:=table(zip(`=`,GraphTheory[Vertices](G),GraphTheory[Neighbors](G))); B:=GraphTheory[CycleBasis](G); N:=nops(B); B:=Vector(N,p->{op(B[p])}); T:=map(ff,B); for iter from 0 to maxiter while (max(map(nops,T))>0) do for p from 1 to N do if (T[p]<>{}) then for q from 1 to N do if (q<>p and `subset`(B[q],B[p])) then B[p]:=`minus`(B[p],`minus`(B[q],T[p])); T[p]:=ff(B[p]) end end end end end; B:=[seq([op(v)],v=B)]; if (iter>maxiter) then error("No convergence, B=%1",B) end; B end: ################################################################################ #cat: Geometry #hfl: Circumsphere Circumsphere:=proc() local n,p0,P,i,p,v,o; n:=_nrest-1; if (n<2) then error("At least three vectors must be provided") end; p0:=_rest[1]; P:=Matrix(n); for i from 1 to n do P[i,..]:=_rest[1+i]-p0 end; try p:=LinearAlgebra[LinearSolve](P,Vector(n,i->add(P[i,o]^2,o=1..n)/2)); catch: return undefined,infinity end; p+p0,sqrt(add(v^2,v=p)) end: ################################################################################ #cat: Calculus #hfl: fullparfrac2 fullparfrac2:=proc(e,z::name,$) local e1,e2,a,subex,indet,eq,sol,z1; e1:=convert(e,fullparfrac,z); e2:=0: if not(type(e1,`+`)) then e1:=[e1] end; for a in e1 do if (type(a,`*`) or type(a,`^`)) then e2:=e2+a else subex:=op(1,a); indet:=op(1,op(2,a)); eq:=op(2,op(2,a)); sol:=allvalues(eq); for z1 in sol do e2:=e2+subs(indet=z1,subex) end end end: e2 end: #hfl: fZero fZero:=proc(f::procedure, r::range(realcons), iniguess::realcons:=infinity, {abserror::nonnegative:=0, relerror::nonnegative:=0, maxiter::posint:=100, output::name:=undefined},$) local x1,x2,x0,iter,flag; x1,x2:=op(r); if (x1=-infinity) then x1:=-10^Digits else x1:=evalf(x1) end; if (x2= infinity) then x2:= 10^Digits else x2:=evalf(x2) end; x0:=evalf(iniguess); if not(x0>=x1 and x0<=x2) then x0:=(x1+x2)/2 end; flag,x0,x2,iter:=FZERO(f,x1,x2,x0,relerror,abserror,maxiter); if not(type(output,undefined)) then assign(output,[flag,x0,x2,iter]) end; if (flag=1 or flag=2) then x0 elif (flag=3) then x0 elif (flag=4) then error "No change in sign was found. After %1 iterations x0=%2, x2=%3, dx=%4, f(x0)=%5", iter,x0,x2,x2-x0,evalf(f(x0)) elif (flag=5) then error "Too many iterations: %1. At this moment x0=%2, x2=%3, dx=%4, f(x0)=%5", iter,x0,x2,x2-x0,evalf(f(x0)) else error "Unrecognized flag: %1",flag end end: #hfl: fZero fDblZero:=proc(f::procedure, r::range(realcons), {abserror::nonnegative:=0, relerror::nonnegative:=0, maxiter::posint:=100, noerror::boolean:=false, rootnumber::posint:=0},$) local x1,x2,f1,f2,s,dx,x0,f0,y1,y2,g1,g2; x1,x2:=min(op(r)),max(op(r)); if (x1=-infinity) then x1:=-10^Digits else x1:=evalf(x1) end; if (x2= infinity) then x2:= 10^Digits else x2:=evalf(x2) end; f1:=evalf(f(x1)); f2:=evalf(f(x2)); s:=signum(f1); if (s*f2<0) then if (rootnumber=1 or rootnumber=2) then fZero(f,x1..x2,':-abserror'=abserror,':-relerror'=relerror,':-maxiter'=maxiter) elif noerror then return 'undefined' else error "f(x1) and f(x2) have different signs" end end; dx:=`if`(abserror>0,abserror,(x2-x1)*max(relerror,10^(2-Digits))); x0:=0.5*(x1+x2); f0:=evalf(f(x0)); while (f0=0.) do x0:=max(0.5*(x0+x1),x0-dx); f0:=evalf(f(x0)) end; if (s*f0>0) then do if (x2-x1x), {abserror::nonnegative:=0, relerror::nonnegative:=0, maxiter::posint:=100},$) local x1,x2,x0,iter,i,lsx,lsy,tb,flag; x1,x2:=op(r); if (x1=-infinity) then x1:=-10^Digits else x1:=evalf(x1) end; if (x2= infinity) then x2:= 10^Digits else x2:=evalf(x2) end; lsx:=evalf(Grid(x1..x2,n,GridF)); lsy:=map(evalf@f,lsx); tb:=table(): for i from 2 to n do if lsy[i-1]*lsy[i]<0 then flag,x0,x2,iter:=FZERO(f,lsx[i-1],lsx[i],(lsx[i-1]+lsx[i])/2,relerror,abserror,maxiter); if (flag=1 or flag=2 or flag=3) then tb[i]:=x0 else error "Cannot find root on interval %1",lsx[i-1]..lsx[i] end end end; convert(tb,list) end: #hfl: fZero FZERO:=proc(F,B0,C0,R0,RE,AE,MAXITER) # based on FORTRAN subroutine of SLATEC package # original authors: L. F. Shampine and H. A. Watts local A,ACBS,ACMB,B,C,CMB,FA,FB,FC,FX,FZ,P,Q,R,T,TOL,Z,IC,ITER; B:=evalf(B0); C:=evalf(C0); R:=evalf(R0); TOL:=`if`( AE>0, AE, abs(B-C)*max(RE,10^(2-Digits)) ); Z:=`if`( R<=min(B,C) or R>=max(B,C), C, R ); IC:=0; FZ:=evalf(F(Z)); FC:=FZ; FB:=evalf(F(B)); if (signum(FZ)<>signum(FB)) then C:=Z elif (Z<>C) then FC:=evalf(F(C)); if (signum(FZ)<>signum(FC)) then B:=Z; FB:=FZ end end; A:=C; FA:=FC; ACBS:=abs(B-C); FX:=max(abs(FB),abs(FC)); ITER:=0; ##### MAIN LOOP do # Perform interchange if needed if (abs(FC)FX) then return 3,B,C,ITER else return 1,B,C,ITER end elif (FB=0) then return 2,B,C,ITER elif (ITER>=MAXITER) then return 5,B,C,ITER end; # Calculate new iterate implicitly as B+P/Q, where we arrange P>=0. The implicit form is used to prevent overflow P:=(B-A)*FB; Q:=FA-FB; if (P<0) then P:=-P; Q:=-Q end; # Update A and check for satisfactory reduction in the size of the bracketing interval. If not, perform bisection A:=B; FA:=FB; IC:=IC+1; if (IC>=4 and 8*ACMB>=ACBS) then # Use bisection (C+B)/2 B:=B+CMB else if (IC>=4) then IC:=0; ACBS:=ACMB end; # Test for too small a change if (P>abs(Q)*TOL) then # Root ought to be between B and (C+B)/2 if (P>=CMB*Q) then # Use bisection (C+B)/2 B:=B+CMB else # Use secant rule B:=B+P/Q end else # Increment by TOLerance. B:=B+signum(CMB)*TOL end end; # Have completed computation for new iterate B FB:=evalf(F(B)); ITER:=ITER+1; # Decide whether next step is extrapolation or interpolation if (signum(FB)=signum(FC)) then C:=A; FC:=FA end end do end: #hfl: NLPSolve2 NLPSolve2:=proc( n::posint, f::procedure, g::{0,procedure}:=0, dx::positive:=10^(-iquo(Digits,2)), {onesidegradient::boolean:=false} ) local grad; if type(g,procedure) then grad:=g elif onesidegradient then grad:=proc(x) local f0; f0:=evalf(f(x)); Vector(n,i->evalf((f(Vector(n,j->x[j]+`if`(j=i,dx,0)))-f0)/dx),datatype=float) end else grad:=x->Vector(n,i->evalf((f(Vector(n,j->x[j]+`if`(j=i,dx,0)))-f(Vector(n,j->x[j]-`if`(j=i,dx,0))))/(2*dx)),datatype=float) end; Optimization[NLPSolve]( n, f, _rest, objectivegradient=proc(V,G) local u,i; u:=grad(V); for i from 1 to n do G[i]:=u[i] end end) end: #hfl: HessianByGradient HessianByGradient:=proc( n::posint, g::procedure, x::indexable, dx::positive:=10^(-iquo(Digits,3)), {printout::boolean:=false},$) local x0,H,j,dxj,g1,g2,i; x0:=Vector(n,i->evalf(x[i]),datatype=float); H:=Matrix(n,datatype=float); for j from 1 to n do dxj:=Vector(n,datatype=float); dxj[j]:=dx; g1:=evalf(g(x0+dxj)); g2:=evalf(g(x0-dxj)); for i from 1 to n do H[i,j]:=(g1[i]-g2[i])/(2*dx) end end; if printout then printf("Frobenius asymmetry is %.2g, maximum asymmetry is %.2g\n", sqrt(add(add((H[i,j]-H[j,i])^2,j=i+1..n),i=1..n-1)), max(seq(seq(abs(H[i,j]-H[j,i]),j=i+1..n),i=1..n-1)) ) end; Matrix(n,(i,j)->(H[i,j]+H[j,i])/2,datatype=float,shape=symmetric); end: #hfl: PolynomInterp1der PolynomialInterpolation1der:=proc(xydata::list,x,k::nonnegint,ypk) local n,L,X,Y,i,j; n:=nops(xydata); L:=[seq(i,i=1..n)]; X,Y:=map2(op,1,xydata),map2(op,2,xydata); add(Y[i]*(x-X[k])/(X[i]-X[k])*mul((x-X[j])/(X[i]-X[j]),j=subsop(i=NULL,L)),i=subsop(k=NULL,L)) +(Y[k]+(ypk-Y[k]*add(1/(X[k]-X[j]),j=subsop(k=NULL,L)))*(x-X[k]))*mul((x-X[j])/(X[k]-X[j]),j=subsop(k=NULL,L)) end: #hfl: BezierControlPoints BezierControlPoints:=proc( x1::realcons,y1::realcons,p1::realcons,q1::realcons,k1::realcons, x4::realcons,y4::realcons,p4::realcons,q4::realcons,k4::realcons) local v1,v4,a,b1,b4,e,L; a:=2*(p1*q4-p4*q1); b1:=6*(q1*(x4-x1)-p1*(y4-y1)); b4:=6*(q4*(x4-x1)-p4*(y4-y1)); L:=evalf(sqrt((x4-x1)^2+(y4-y1)^2)); e:=fsolve({k1*v1^2+a*v4=-b1,k4*v4^2+a*v1=b4},{v1=L,v4=L},{v1=0..infinity,v4=0..infinity}); if (nops([e])>1) then WARNING("Multiple solution for velocities: %1. They are set to L=%2. The input was: %3", [e],L,[x1,y1,p1,q1,k1,x4,y4,p4,q4,k4]); e:={v1=L,v4=L} elif (type(e,function)) then WARNING("Cannot solve equations for velocities: %1. They are set to L=%2. The input was: %3", {k1*v1^2+a*v4=-b1,k4*v4^2+a*v1=b4},L,[x1,y1,p1,q1,k1,x4,y4,p4,q4,k4]); e:={v1=L,v4=L} end; evalf(subs(op(e),[[x1,y1],[x1+v1*p1/3,y1+v1*q1/3],[x4-v4*p4/3,y4-v4*q4/3],[x4,y4]])) end: #hfl: BezierNodes BezierNodes:=proc(Points::list, B::list:=[0,0]) local Nodes,n,i,i1,i2,calc,ypv,kv; # Subroutine calc:=proc(k::integer:=0) local x,dx,y,yp,c2,cosa,e; x:=Points[i][1]; e:=`if`(k=0, CurveFitting[PolynomialInterpolation](map(v->[v[1]-x,v[2]],Points[max(i-2,1)..min(i+2,n)]),dx), PolynomialInterpolation1der(map(v->[v[1]-x,v[2]],Points[max(i-2,1)..min(i+2,n)]),dx,k,ypv)); y,yp,c2:=seq(coeff(e,dx,i),i=0..2); cosa:=1/sqrt(1+yp^2); Nodes[i]:=evalf([x,y,cosa,yp*cosa,2*c2*cosa^3]) end; # Begin program n:=nops(Points); Nodes:=Array(1..n); for i from 1 to n do calc() end; i1,i2:=1,n; # Left boundary if (B[1] in [0,1,2]) then i1:=1+B[1] elif (B[1]=3) then ypv:=args[3][1]; for i from 1 to 2 do calc(1) end elif (B[1]=4) then ypv,kv:=args[3][1],args[4][1]; Nodes[1]:=evalf([op(Points[1]),1/sqrt(1+ypv^2),ypv/sqrt(1+ypv^2),kv]); i:=2; calc(1) else error "Unrecognized left boundary condition:",B[1] end; # Right boundary if (B[2] in [0,1,2]) then i2:=n-B[2] elif (B[2]=3) then ypv:=args[3][2]; for i from n-1 to n do calc(n) end elif (B[2]=4) then ypv,kv:=args[3][2],args[4][2]; Nodes[n]:=evalf([op(Points[n]),1/sqrt(1+ypv^2),ypv/sqrt(1+ypv^2),kv]); i:=n-1; calc(n) else error "Unrecognized right boundary condition:",B[2] end; [seq(Nodes[i],i=i1..i2)] end: #hfl: ExtrapolateMPE ExtrapolateMPE:=proc(X::Matrix,$) local n,m,dX,A,i,b; n,m:=op(1,X); m:=m-1; dX:=Matrix(n,m,(i,j)->X[i,j+1]-X[i,j],datatype=float): A:=LinearAlgebra[Transpose](dX).dX; for i from 1 to m do A[m,i]:=1 end; b:=Vector(m,datatype=float): b[m]:=1: X[..,..m].LinearAlgebra[LinearSolve](A,b) end: ################################################################################ #cat: Differential equations #hfl: devpsolve devpsolve:=proc( x1::realcons, x2::realcons, BC, U::procedure(anything), N::posint:=100, { numev::posint:=iquo(N+1,2), evonly::boolean:=false },$) local X,dx,dt,H,i,ev,evc0,evc,p,c,S; X:=Vector(N,i->x1+(x2-x1)*(i-1)/(N-1)); dx:=evalf((x2-x1)/(N-1)); if type(BC,complexcons) then dt:=`if`(BC=1 or BC=-1,float,complex); H:=Matrix(N-1,shape=hermitian); for i from 1 to N-1 do H[i,i]:=evalf(2+U(X[i])*dx^2) end; for i from 1 to N-2 do H[i,i+1]:=-1. end; H[N-1,1]:=-evalf(BC); if evonly then return (1/dx^2)*(LinearAlgebra[Eigenvalues](H))[1..numev] end; ev,evc0:=LinearAlgebra[Eigenvectors](H); evc:=Matrix(N,numev,datatype=dt); evc[1..-2,1..-1]:=evc0[1..-1,1..numev]; evc[N,1..-1]:=evalf(BC).evc0[1,1..numev] elif type(BC,[anything,anything]) then H:=Matrix(N-2,shape=symmetric,datatype=float); for i from 1 to N-2 do H[i,i]:=evalf(2+U(X[i+1])*dx^2) end; for i from 1 to N-3 do H[i,i+1]:=-1. end; for p from 1 to 2 do if (BC[p]=1) then c[p]:=0 elif (BC[p]=2) then c[p]:=1; elif (type(BC[p],list) and type(BC[p][1],positive)) then c[p]:=2/(1+BC[1][1]*dx)-1 else error "Unrecognized %2 boundary condition: %1",BC[p],`if`(p=1,"left","right") end end; H[ 1, 1]:=H[ 1, 1]-c[1]; H[-1,-1]:=H[-1,-1]-c[2]; if evonly then return (1/dx^2)*(LinearAlgebra[Eigenvalues](H))[1..numev] end; ev,evc0:=LinearAlgebra[Eigenvectors](H); evc:=Matrix(N,numev,datatype=float); evc[2..-2,1..-1]:=evc0[1..-1,1..numev]; evc[1,1..-1]:=c[1].evc0[ 1,1..numev]; evc[N,1..-1]:=c[2].evc0[-1,1..numev] else error "Unrecognized boundary condition: %1",BC end; S:=LinearAlgebra[DiagonalMatrix]([seq(`if`(evc[2,i]>=0,1,-1),i=1..numev)]); (1/dx^2)*ev[1..numev],(1/sqrt(dx))*evc.S,evalf(X) end: ################################################################################ #cat: Linear algebra #hfl: CheckEigenvector CheckEigenvector:=proc(M::Matrix,V::Vector,maxdev::numeric:=0,simplifyf::procedure:=simplify,{trueorfalse::boolean:=false},$) local v,R; v:=simplifyf(LinearAlgebra[HermitianTranspose](V).M.V/simplifyf(add(abs(v)^2,v=V))); R:=map(simplifyf,simplify((M-v).V)); if `if`(maxdev=0, verify(R,LinearAlgebra[ZeroVector](op(1,V)),Vector), evalb(LinearAlgebra[Norm](R,2)simplify(u/v),V); break end end end; try u:=ilcm(seq(`if`(type(v,`+`),denom(op(1,v)),denom(v)),v=V)) catch: u:=1 end; if (type(u,realcons) and u<>1) then LinearAlgebra[VectorScalarMultiply](V,u,'inplace') end; try u:=igcd(seq(subs(sbs1,v),v=V)) catch: u:=1 end; if (type(u,realcons) and u<>1) then LinearAlgebra[VectorScalarMultiply](V,1/u,'inplace') end elif (_nrest>0) then LinearAlgebra[Normalize](V,_rest,'inplace') end; if positive then LinearAlgebra[VectorScalarMultiply](V,signum(V[op(MaxIdx(evalf(V),v->abs(v)))]),'inplace') end; V end: #hfl: NormalizeMatrix NormalizeMatrix:=proc(M::Matrix,{positive::boolean:=false,int::boolean:=false,rad::boolean:=false}) local j; for j from 1 to Dim2(M)[2] do M[..,j]:=NormalizeVector(M[..,j],_passed[2..]) end; M end: #hfl: RealJordanBlockMatrix RealJordanBlockMatrix:=proc(lsJB::list,$) local ls,JB,ev,m,b1,b2,i; ls:=[]: for JB in lsJB do ev,m:=op(JB); if Im(ev)=0 then ls:=[op(ls),LinearAlgebra[BandMatrix]([[ev$m],[1$(m-1)]],0)] else b1:=[ Im(ev),seq(op([0, Im(ev)]),i=1..m-1)]; b2:=[-Im(ev),seq(op([0,-Im(ev)]),i=1..m-1)]; ls:=[op(ls),LinearAlgebra[BandMatrix]( [b1,[Re(ev)$(2*m)],b2,[1$(2*m-2)]],1)] end end; LinearAlgebra[DiagonalMatrix](ls) end: #hfl: MatrixFunction2 MatrixFunction2:=proc(M::Matrix,f::procedure,simplifyf::procedure:=(x->x),$) local ev,evc; ev,evc:=LinearAlgebra[Eigenvectors](M); map(simplifyf,evc.LinearAlgebra[DiagonalMatrix](map(f,ev)).evc^(-1)) end: #hfl: SparseMM SparseMM:=proc(M1::Matrix,M2::Matrix,$) local T1,K,T2,js,M,i,j,k,v; T1:=table(op(2,M1)); K:=map(v->map2(op,2,v),Classify2([indices(T1)],[1])); T2:=table(op(2,M2)); js:=map2(op,2,[indices(T2)]); M:=Matrix(Dim2(M1)[1],Dim2(M2)[2]); for i in indices(K,'nolist') do for j in js do M[i,j]:=add(`if`(assigned('T2[k,j]'),T1[i,k]*T2[k,j],0),k=K[i]) end end; M end: #hfl: SparseMM SparseMHM:=proc(M::Matrix,H::{Matrix,undefined}:=undefined,{conjugate::boolean:=false},$) local ff,T,J,js,S,Sshape,pi,i,j,k,l,v; if conjugate then ff:=:-conjugate else ff:=v->v end; T:=table(op(2,M)); J:=map(v->map2(op,1,v),Classify2([indices(T)],[2])); js:=sort([indices(J,'nolist')]); S:=Matrix(Dim2(M)[2]); Sshape:=`if`(type(H,Matrix) and MatrixOptions(H,shape)<>[symmetric] and MatrixOptions(H,shape)<>[hermitian],NULL,`if`(conjugate,hermitian,symmetric)); if (Sshape=NULL) then for i in js do for j in js do S[i,j]:=add(add(ff(T[k,i])*H[k,l]*T[l,j],k=J[i]),l=J[j]) end end else if type(H,Matrix) then for pi from 1 to nops(js) do i:=js[pi]; for j in js[pi..] do S[i,j]:=add(add(ff(T[k,i])*H[k,l]*T[l,j],k=J[i]),l=J[j]) end end else for pi from 1 to nops(js) do i:=js[pi]; for j in js[pi..] do S[i,j]:=add(`if`(assigned('T[k,j]'),ff(T[k,i])*T[k,j],0),k=J[i]) end end end; for pi from 1 to nops(js) do i:=js[pi]; for j in js[..pi-1] do S[i,j]:=ff(S[j,i]) end end end; S end: #hfl: SparseMM SparseMV:=proc(M::Matrix,V::Vector,$) local T,J,U,i,j,v; T:=table(op(2,M)); J:=map(v->map2(op,2,v),Classify2([indices(T)],[1])); U:=Vector(Dim2(M)[1]); for i in indices(J,'nolist') do U[i]:=add(T[i,j]*V[j],j=J[i]) end; U end: #hfl: SparseMM SparseMO:=proc(MM::Matrix,norm::{"","2","int"},sf::procedure:=(x->x),{sortbysize::boolean:=false},$) local ls,tb,j0,J,P,L,U,m,p,U1,U2,J2,i,j,v; ls:=LinearAlgebra[StronglyConnectedBlocks](MM,'outputform'='rowlist'); if sortbysize then ls:=Sort(ls,nops) end; tb:=table(); j0:=0; for J in ls do if (nops(J)=1) then U:=MM[J[1],J[1]]; if (U<>0) then tb[j0]:=[(J[1],j0+1)=`if`(norm="2",1/sqrt(U),1)]; j0:=j0+1 end else P,L,U:=LinearAlgebra[LUDecomposition](Matrix(MM[J,J])); m:=max(seq(lhs(v)[1],v=op(2,U))); #Rank p:=map2(op,1,Sort([seq([lhs(v)],v=op(2,P))],[2]))[..m]; #permutation U1:=U[..m,p]; if (norm="2") then for i from 1 to m do LinearAlgebra[RowOperation](U1,i,1/sqrt(U1[i,i]),'inplace') end end; U2:=map(sf,U1^(-1)); if (norm="int") then U2:=NormalizeMatrix(U2,'int') end; J2:=J[p]; # M2=M[..,J2].U2, HermitianTranspose(U2).MM[J2,J2].U2=`if`(norm="2",1,diagonalmatrix) tb[j0]:=[seq(seq(`if`(U2[i,j]=0,NULL,(J2[i],j0+j)=U2[i,j]),j=i..m),i=1..m)]; j0:=j0+m end end; Matrix(Dim2(MM)[1],j0,map(op,{entries(tb,'nolist')})) end: #hfl: SparseMM SparseMN:=proc(M::Matrix) local T,J,MN,K,V,i,j,k,v; T:=table(op(2,M)); J:=map(v->map2(op,1,v),Classify2([indices(T)],[2])); MN:=Matrix(Dim2(M)); for j in indices(J,'nolist') do K:=J[j]; V:=NormalizeVector(Vector([seq(M[i,j],i=K)]),_rest); for k from 1 to nops(K) do MN[K[k],j]:=V[k] end end; MN end: #hfl: SparseMM SparseMap:=proc(f,A) local v,u; u:=_rest; op(0,A)(op(1,A),map(v->lhs(v)=f(rhs(v),u),op(2,A)),op(3,A)) end: #hfl: SuperimposeM SuperimposeM:=proc( M1::Matrix, M2::Matrix, out::string:="m", {printout::boolean:=false},$) local n,m,U1,S1,Vt1,R,f,d,v,ijs,e,c; n,m:=Dim2(M1); if (Dim2(M1)<>Dim2(M2)) then error("Different matrix sizes: %1<>%2",Dim2(M1)<>Dim2(M2)) end; U1,S1,Vt1:=LinearAlgebra[SingularValues](M1.LinearAlgebra[Transpose](M2),'output'=['U','S','Vt']); R:=U1.Vt1; f:=R.M2-M1; d:=sqrt(add(v^2,v=f)/m); # Frobenius_norm/sqrt(m) ijs:=SortIdx(convert(f,array),v->-abs(v)); ijs:=ijs[..min(5,nops(ijs))]; e:=f[op(ijs[1])]; if printout then printf("RMSD = %.3g\n",d); printf("Deviations:%{c,}s\n",Vector(map(ij->sprintf(" [%d,%d]=%.2g",op(ij),f[op(ij)]),ijs))); printf("Rotation: det=%d, angle=%.3g, axis=<%{c,}.3g>\n",RotationParam(R,'degrees')) end; seq(`if`(c="m",R.M2,`if`(c="d",d,`if`(c="e",e,`if`(c="f",f,`if`(c="r",R,NULL))))),c=out) end: #hfl: LocalizeVectors4 LocalizeVectors4:=proc(M0::Matrix,tol::numeric,{method::{"sequential","parallel"}:="sequential",printout::{boolean,posint}:=false},$) local prn,freq,a0,d,n,m,M,S4,S3,A,i,j,k,l,a,c,s,pair, o,p; if type(printout,posint) then prn,freq:=true,printout else prn,freq:=printout,10 end; a0:=1; d:=1+max(0,floor(-log10(tol))); n,m:=op(1,M0); if (m=1) then return M0 end; M:=copy(M0); if prn then printf("%-*s localization\n",d+2,"angle") end; if (method="sequential") then S4:=Matrix(m,(i,j)->evalhf(add(M[p,i]^2*M[p,j]^2,p=1..n)),shape=symmetric,datatype=float); S3:=Matrix(m,(i,j)->evalhf(add(M[p,i]*M[p,j]*(M[p,i]^2-M[p,j]^2),p=1..n)),shape=antisymmetric,datatype=float); A:=Matrix(m,(i,j)->`if`(iabs(v))); a:=A[i,j]; if (prn and abs(a)|<-s,c>>; for k in [i,j] do for l from 1 to m do S4[k,l]:=evalhf(add(M[p,k]^2*M[p,l]^2,p=1..n)); S3[k,l]:=evalhf(add(M[p,k]*M[p,l]*(M[p,k]^2-M[p,l]^2),p=1..n)) end end; for k in [i,j] do for l from 1 to m do A[op(sort([k,l]))]:=evalhf(signum(l-k)*arctan(4*S3[k,l],(S4[k,k]+S4[l,l]-6*S4[k,l]))/4) end end end else do S4:=Matrix(m,(i,j)->add(M[p,i]^2*M[p,j]^2,p=1..n),shape=symmetric,datatype=float); A:=table([]): for i from 1 to m do for j from i+1 to m do A[i,j]:=arctan(4*add(M[p,i]*M[p,j]*(M[p,i]^2-M[p,j]^2),p=1..n),S4[i,i]+S4[j,j]-6*S4[i,j])/4 end end; a:=abs(MaxVal(A,v->abs(v))); if (prn and abs(a)abs(v))); c,s:=cos(A[i,j]),sin(A[i,j]); M[..,[i,j]]:=M[..,[i,j]].<|<-s,c>>; for k in [i,j] do for l from 1 to m do unassign('A[op(sort([k,l]))]') end end end end end; M end: ################################################################################ #cat: Linear transformations #hfl: RotationM RotationM:=proc(alpha)::Matrix; local n,L,R,theta,phi; if (nargs=1) then R:=<,> else if (nargs=2 and type(args[2],indexable)) then n:=LinearAlgebra[Normalize](convert(args[2],Vector),2,conjugate=false) elif (nargs=3) then if (type(args[2],indexable) and args[3]='normalized') then n:=args[2] else theta,phi:=args[2],args[3]; n:= end else error "Unrecognized arguments: %1",[args] end; L:=Matrix(3,{(3,2)=n[1],(1,3)=n[2],(2,1)=n[3]},shape=antisymmetric); R:=1+sin(alpha)*L+(1-cos(alpha))*L^2 end; `if`(hastype([args],'float'),evalf(R),R) end: #hfl: ReflectionM ReflectionM:=proc()::Matrix; local phi,R; if (nargs=1 and not(type(args[1],indexable))) then phi:=args[1]; R:=<<-cos(2*phi)|-sin(2*phi)>,<-sin(2*phi)|cos(2*phi)>>; `if`(type(phi,'float'),evalf(R),R) else -RotationM(Pi,args) end end: #hfl: RotationParam RotationParam:=proc( R::Matrix(3,3), P::Vector(3):=<0,0,1>, T::name:='undefined', { method::{"Eigenvectors","NullSpace"}:=`if`(hastype(R,'float'),"Eigenvectors","NullSpace"), degrees::boolean:=false, digits::posint:=2 },$) local det,ev,evc,i,n,alpha,e,e1,e2,T1,v; det:=simplify(LinearAlgebra[Determinant](R)); if (evalb(abs(abs(det)-1)>10^(digits-Digits))=true) then error("Impossible determinant: %1",det) else det:=signum(det) end; if (method="Eigenvectors") then ev,evc:=LinearAlgebra[Eigenvectors](det*R); i:=op(MinIdx(ev,v->abs(v-1))); n:=map(Re,evc[..,i]); alpha:=argument(simplify(ev[(i mod 3)+1])) else e:=LinearAlgebra[NullSpace](R-det); if (nops(e)<>1) then if (add(abs(v),v=R-det)=0) then return det,0,P else error "Impossible nullspace",e end end; n:=op(e); e:=simplify((det*LinearAlgebra[Trace](R)-1)/2); if (evalb(abs(e)>1+1e-4)=true) then error "Impossible cosine",e end; alpha:=arccos(e) end; e1:=LinearAlgebra[MatrixNorm](simplify(R-det*RotationM( alpha,n)),'Frobenius'); e2:=LinearAlgebra[MatrixNorm](simplify(R-det*RotationM(-alpha,n)),'Frobenius'); if (evalf(e1)>evalf(e2)) then alpha:=-alpha end; if (min(evalf(e1),evalf(e2))>10^(digits-Digits)) then error("Failed test: %1",min(evalf(e1),evalf(e2))) end; if (LinearAlgebra[DotProduct](n,P)<0) then n,alpha:=-n,-alpha end; if not(T='undefined') then T1:=LinearAlgebra[IsSimilar](RotationM(alpha,0,0),R,output='C'); if hastype(R,'float') then T1:=evalf(T1) end; assign(T,T1) end; if degrees then alpha:=180/Pi*alpha; if hastype(R,'float') then alpha:=evalf(alpha) end end; det,alpha,n end: #hfl: RotationParam RotationParamEulerX:=proc( R::Matrix(3,3), det::{1,-1}:=signum(LinearAlgebra[Determinant](R)), $) `if`(has(R,float),evalf,simplify)([arctan(det*R[1,3],-det*R[2,3]),arccos(det*R[3,3]),arctan(det*R[3,1],det*R[3,2])]) end: #hfl: RotationParam RotationParamEulerY:=proc( R::Matrix(3,3), det::{1,-1}:=signum(LinearAlgebra[Determinant](R)), $) `if`(has(R,float),evalf,simplify)([arctan(det*R[2,3],det*R[1,3]),arccos(det*R[3,3]),arctan(det*R[3,2],-det*R[3,1])]) end: unassign('phi,theta,psi,alpha,beta,gamma'); #hfl: RotationMEulerX RotationMEulerX:=unapply(RotationM(phi,<0,0,1>).RotationM(theta,<1,0,0>).RotationM(psi,<0,0,1>),phi,theta,psi): #hfl: RotationMEulerX RotationMEulerY:=unapply(RotationM(alpha,<0,0,1>).RotationM(beta,<0,1,0>).RotationM(gamma,<0,0,1>),alpha,beta,gamma): #hfl: EulerY2Rot EulerY2Rot:=proc(alpha,beta,gamma,$) simplify(2*arccos(cos(beta/2)*cos((alpha+gamma)/2))), `if`(beta=0,(1-signum(sin((alpha+gamma)/2)))/2*Pi,simplify(arccot(cot(beta/2)*sin((alpha+gamma)/2)))), simplify(Reduce2P((alpha-gamma+Pi)/2,2*Pi)) end: #hfl: EulerY2Rot Rot2EulerY:=proc(alpha,e) local theta,phi,s; if (_npassed=2) then theta:=simplify(arccos(e[3]/sqrt(e[1]^2+e[2]^2+e[3]^2))); phi:=simplify(arctan(e[2],e[1])) elif (_npassed=3) then theta:=e; phi:=_passed[3] else error("Unrecognized input: expected (alpha,e) or (apha,theta,phi), but received %1",_passed) end; s:=simplify(arctan(cos(theta)*sin(alpha/2),cos(alpha/2))); simplify(Reduce2P(s+phi-Pi/2,2*Pi)), simplify(2*arcsin(sin(theta)*sin(alpha/2))), simplify(Reduce2P(s-phi+Pi/2,2*Pi)) end: #hfl: RotationMp2q RotationMp2q:=proc(p::Vector,q::Vector,{normalized::boolean:=false},$) local npnq,L; npnq:=`if`(normalized,1,LinearAlgebra[Norm](p,2)*LinearAlgebra[Norm](q,2)); L:=(q.LinearAlgebra[Transpose](p)-p.LinearAlgebra[Transpose](q))/npnq; 1+L+L^2/(1+LinearAlgebra[Transpose](p).q/npnq) end: #hfl: RotationM2zx RotationM2zx:=proc(p::indexable,q::indexable,axes::{12,13,23,21,32,31}:=31,det::{0,1,-1}:=1,$) local i,j,k,e; i,j:=seq(parse(v),v=convert(axes,string)); k:=6-i-j; e:=Vector(3); e[i]:=LinearAlgebra[Normalize](convert(p,Vector),2); e[k]:=LinearAlgebra[Normalize](LinearAlgebra[CrossProduct](e[i],convert(q,Vector)),2); e[j]:=LinearAlgebra[CrossProduct](e[k],e[i]); if (det=1 and member(axes,[13,32,21]) or det=-1 and member(axes,[12,23,31])) then e[k]:=-e[k] end; Matrix(3,(i,j)->e[i][j]); end: #hfl: TransSpherCoo TransSpherCoo:=(t,p,a,b,g)->( arccos(cos(t)*cos(b)+sin(t)*sin(b)*cos(p-a)), arctan(sin(t)*sin(p-a),sin(t)*cos(b)*cos(p-a)-cos(t)*sin(b))-g ): #hfl: TransSpherCoo TransSpherCooI:=(t,p,a,b,g)->( arccos(cos(t)*cos(b)-sin(t)*sin(b)*cos(p+g)), arctan(sin(t)*sin(p+g),sin(t)*cos(b)*cos(p+g)+cos(t)*sin(b))+a ): #hfl: TransTensor TransTensor:=proc(A::indexable,T::Matrix) local dims,B,V,S,k,i,j; B:=Array(A); dims:=[ArrayDims(B)]; for k from 1 to nops(dims) do V:=Vector(dims[k],i->B[op(subsop(k=i,dims))]); for i in [$dims[k]] do S:=0; for j in [$dims[k]] do S:=S+T[i,j]*V[j] end; B[op(subsop(k=i,dims))]:=S end end; B end: #hfl: ExtendedTransM ExtendedTransM:=proc(T::Matrix) local ET,i,j; ET:=Matrix(6); for i from 1 to 3 do for j from 1 to 3 do ET[i,j]:=T[i,j]^2; ET[i,j+3]:=2*T[i,(j mod 3)+1]*T[i,(j+1 mod 3)+1]; ET[i+3,j]:=T[(i mod 3)+1,j]*T[(i+1 mod 3)+1,j]; ET[i+3,j+3]:=T[(i mod 3)+1,(j mod 3)+1]*T[(i+1 mod 3)+1,(j+1 mod 3)+1] +T[(i mod 3)+1,(j+1 mod 3)+1]*T[(i+1 mod 3)+1,(j mod 3)+1] end end; ET end: ################################################################################ #cat: Spherical harmonics #hfl: Phim Phim:=proc(m::integer,phi,{complex::boolean:=complexTrigBasis},$) if complex then exp(I*m*phi) else if (m=0) then 1 elif (m>0) then cos(m*phi) else sin(abs(m)*phi) end end end: #hfl: Phim NPhim2:=proc(m::integer,{complex::boolean:=complexTrigBasis},$) `if`(complex or m=0,2*Pi,Pi) end: #hfl: Pnma Pnma:=proc(n::nonnegint,m::nonnegint,arg) local x,v,theta,a; if (m>n) then error "m must be not greater than n",m,n end; if (nargs=3) then theta:=args[3]; `if`(m=0, orthopoly['P'](n,cos(theta)), simplify(subs(x=cos(theta),diff(orthopoly['P'](n,x),x$m))*v)/v*sin(theta)^m) else a,theta:=args[3],args[4]; if (evalb(a<=-1/2)=true or evalb(evalf(a)<=-0.5)=true) then error "a must be greater than -1/2",a end; `if`(m=0, orthopoly['G'](n,a+1/2,cos(theta)), simplify(subs(x=cos(theta),diff(orthopoly['G'](n,a+1/2,x),x$m))*v)/v*sin(theta)^m) end end: #hfl: Pnma NPnma2:=(n,m,a)->2*Pi/4^a/GAMMA(a+1/2)^2*GAMMA(2*a+1+n+m)/(2*a+1+2*n)/(n-m)!: #hfl: Ylm Ylm:=proc(l0,m::integer,t0,phi, convention::string:=conventionY, {normalized::boolean:=normalizedY},$) local l,t,e,k,A; if type(l0,list) then l:=ListTools[Reverse](l0) else l:=[l0] end; if type(t0,list) then t:=ListTools[Reverse](t0) else t:=[t0] end; if (nops(l)<>nops(t)) then error "Inconsistent l and theta",l0,t0 end; for k from 1 to nops(l) do if not(type(l[k],nonnegint)) then error "l must be nonnegative integer",l0 end end; for k from 2 to nops(l) do if (l[k]0 and l[1]0,(-1)^m,1)*Phim(m,phi,'complex') elif (convention="complex") then A:=Phim(m,phi,'complex') elif (convention="real") then A:=Phim(m,phi,'complex'=false) else error "Convention can be one of [standard,complex,real] but received %1",convention end; e:=mul(factor(Pnma(l[k],`if`(k=1,abs(m),l[k-1]),(k-1)/2,t[k])),k=1..nops(l))*A; `if`(normalized, e*simplify(1/sqrt(mul(NPnma2(l[k],`if`(k=1,abs(m),l[k-1]),(k-1)/2),k=1..nops(l))*NPhim2(m,'complex'=evalb(convention<>"real")))), e) end: #hfl: Yxyz Yxyz:=proc(l::nonnegint,m::{integer,string},x,y,z,$) local e,r,t,p,m2,i; if type(m,integer) then # e:=algsubs(r^2=x^2+y^2+z^2,algsubs(cos(t)=z/r,simplify(Ylm(l,m,t,p,'normalized')*exp(-m*I*p))*((x+signum(m)*I*y)/sin(t))^abs(m))*r^(l-abs(m))); # simplify(factor(`if`(m=0,1,sqrt(2))*evalc(`if`(m<0,Im,Re)(e)))) e:=algsubs(r^2=x^2+y^2+z^2, simplify( r^(l-abs(m))* algsubs(cos(t)=z/r, algsubs(`if`(m=0,t=t,`if`(m<0,sin(t)^abs(m)*sin(abs(m)*p)=((x+I*y)^abs(m)-(x-I*y)^abs(m))/2/I,sin(t)^m*cos(m*p)=((x+I*y)^m+(x-I*y)^m)/2)), Ylm(l,m,t,p,"real",'normalized'))))) elif (m="trig") then [seq(Yxyz(l,m2,x,y,z),m2=[0,seq(op([i,-i]),i=1..l)])] elif (m="cubic") then [seq(op(Yxyz(l,m2,x,y,z)),m2=YxyzR(l))] elif (l=0) then if (m="A1") then [ 1/sqrt(4*Pi) ] else error("For l=0 m can be one of [A1], but received %1",m) end elif (l=1) then if (m="F1") then [ x*sqrt(3/4/Pi), y*sqrt(3/4/Pi), z*sqrt(3/4/Pi) ] else error("For l=1 m can be one of [F1], but received %1",m) end elif (l=2) then if (m="E" ) then [ (x^2-y^2)*sqrt(15/16/Pi), (2*z^2-x^2-y^2)*sqrt(5/16/Pi) ] elif (m="F2") then [ y*z*sqrt(15/4/Pi), z*x*sqrt(15/4/Pi), x*y*sqrt(15/4/Pi) ] else error("For l=2 m can be one of [E,F2], but received %1",m) end elif (l=3) then if (m="A2") then [ x*y*z*sqrt(105/4/Pi) ] elif (m="F1") then [ x*(2*x^2-3*y^2-3*z^2)*sqrt(7/16/Pi), y*(2*y^2-3*z^2-3*x^2)*sqrt(7/16/Pi), z*(2*z^2-3*x^2-3*y^2)*sqrt(7/16/Pi) ] elif (m="F2") then [ x*(y^2-z^2)*sqrt(105/16/Pi), y*(z^2-x^2)*sqrt(105/16/Pi), z*(x^2-y^2)*sqrt(105/16/Pi) ] else error("For l=3 m can be one of [A2,F1,F2], but received %1",m) end elif (l=4) then if (m="A1") then [ (x^4+y^4+z^4-3*x^2*y^2-3*y^2*z^2-3*z^2*x^2)*sqrt(21/16/Pi) ] elif (m="E" ) then [ (x^2-y^2)*(6*z^2-x^2-y^2)*sqrt(45/64/Pi), (x^4+y^4-2*z^4-12*x^2*y^2+6*x^2*z^2+6*y^2*z^2)*sqrt(15/64/Pi) ] elif (m="F1") then [ y*z*(y^2-z^2)*sqrt(315/16/Pi), z*x*(z^2-x^2)*sqrt(315/16/Pi), x*y*(x^2-y^2)*sqrt(315/16/Pi) ] elif (m="F2") then [ y*z*(6*x^2-y^2-z^2)*sqrt(45/16/Pi), z*x*(6*y^2-z^2-x^2)*sqrt(45/16/Pi), x*y*(6*z^2-x^2-y^2)*sqrt(45/16/Pi) ] else error("For l=4 m can be one of [A1,E,F1,F2], but received %1",m) end else if member(m,["A1","A2","E","F1","F2"]) then error("Max supported l is 4, but received %1",l) else error("m can be one of [trig,cubic,A1,A2,E,F1,F2], but received %1",m) end end end: #hfl: Yxyz YxyzM:=proc(l::nonnegint,m::{integer,string},convention::string:=conventionY,$) option remember; local e,t,p,b,M; e:=Yxyz(l,m,sin(t)*cos(p),sin(t)*sin(p),cos(t)); b:=type(e,list); if not(b) then e:=[e] end; M:=Matrix(nops(e),2*l+1,(i,j)->combine(simplify(int(int(Ylm(l,iidI(j),t,`if`(convention="real",p,-p),convention,'normalized')*e[i],p=0..2*Pi)*sin(t),t=0..Pi)),sqrt)); `if`(b,M,M[1,..]) end: #hfl: Yxyz YxyzR:=proc(l::nonnegint,$) if (l=0) then ["A1"] elif (l=1) then ["F1"] elif (l=2) then ["E","F2"] elif (l=3) then ["A2","F1","F2"] elif (l=4) then ["A1","E","F1","F2"] else error("Max supported l is 4, but received %1",l) end end: #hfl: Yxyz YxyzL:=proc(l::nonnegint,$) if (l=0) then ["S"] elif (l=1) then ["X","Y","Z"] elif (l=2) then ["X2","Z2","YZ","ZX","XY"] elif (l=3) then ["XYZ","X3","Y3","Z3","XY2","YZ2","ZX2"] elif (l=4) then ["S4","X2Z2","Z4","YZY2","ZXZ2","XYX2","YZX2","ZXY2","XYZ2"] else error("Max supported l is 4, but received %1",l) end end: #hfl: WignerD Wignerd:=proc(l,m1,m2,c,s) local k; sqrt((l+m1)!*(l+m2)!*(l-m1)!*(l-m2)!)*add((-1)^(l-k-m2)/k!/(k+m1+m2)!/(l-k-m1)!/(l-k-m2)!*c^(2*k+m1+m2)*s^(2*l-2*k-m1-m2),k=max(0,-m1-m2)..l-max(m1,m2)) end: #hfl: WignerD WignerD:=(l,m1,m2,alpha,beta,gamma)->exp(I*m1*alpha+I*m2*gamma)*Wignerd(l,m1,m2,cos(beta/2),sin(beta/2)): ### Wigner3j is written by Dennis Isbister, dji@adfa.edu.au #hfl: MatElemY Wigner3j:=proc(j1,j2,j3,m1,m2,m3) local s,ss,f,edmonds; if type([args],list(rational)) then if abs(m1)>j1 or abs(m2)>j2 or abs(m3)>j3 then RETURN (0) fi; if m1+m2+m3 <> 0 then RETURN (0) fi; if j3 < abs(j1-j2) or j3 > j1+j2 then RETURN (0) fi; ss:=0: for s from 0 to j3+m3 do if j1-m1-s < 0 or j3+m3-s < 0 or j2-j3+m1+s < 0 then f:=0 else f:=(j1+m1+s)!*(j2+j3-m1-s)!/(s!*(j1-m1-s)!*(j2-j3+m1+s)!*(j3+m3-s)!): fi; ss:=ss+(-1)**s*f; od; f:=(j1+j2-j3)!*(j1-m1)!*(j2-m2)!*(j3-m3)!*(j3+m3)! /((j1+j2+j3+1)!*(j1-j2+j3)!*(-j1+j2+j3)!*(j1+m1)!*(j2+m2)!): edmonds:=(-1)**(-2*j1-m1-j2-m3)*ss*sqrt(f); simplify(edmonds) else 'Wigner3j'(args) fi; end: ### Wigner6j is written by Dennis Isbister, dji@adfa.edu.au #hfl: MatElemY Wigner6j:= proc(j1,j2,j3,l1,l2,l3) local delta,s,smin,smax,ss,pref,f; delta:= proc(j1,j2,j3) local p; if abs(j1 - j2) <= j3 and j3 <= j1 + j2 then p:=(j1+j2-j3)!*(j1-j2+j3)!*(-j1+j2+j3)!/(j1+j2+j3+1)!; else RETURN(0) fi; RETURN(sqrt(p)) end: if type([args],list(rational)) then pref:=delta(j1,j2,j3)*delta(j1,l2,l3) *delta(l1,j2,l3)*delta(l1,l2,j3): if (pref = 0) then RETURN(0) fi; smin:= max(j1+j2+j3,j1+l2+l3,l1+j2+l3,l1+l2+j3); smax:= min(j1+j2+l1+l2,j2+j3+l2+l3,j3+j1+l3+l1); ss:=0: for s from smin to smax do if s-j1-j2-j3<0 or s-j1-l2-l3<0 or s-l1-j2-l3<0 or s-l1-l2-j3<0 or j1+j2+l1+l2-s<0 or j2+j3+l2+l3-s<0 or j3+j1+l3+l1-s<0 then f:= infinity else f:=(s-j1-j2-j3)!*(s-j1-l2-l3)!*(s-l1-j2-l3)! *(s-l1-l2-j3)!*(j1+j2+l1+l2-s)! *(j2+j3+l2+l3-s)!*(j3+j1+l3+l1-s)!: fi; ss:=ss+(-1)**s*(s+1)!/f: od; pref:=delta(j1,j2,j3)*delta(j1,l2,l3) *delta(l1,j2,l3)*delta(l1,l2,j3): simplify(pref*ss) else 'Wigner6j'(args) fi; end: #hfl: MatElemY MatElemYlm:=(l1,m1,l,m,l2,m2)->simplify((-1)^m1*sqrt((2*l1+1)*(2*l+1)*(2*l2+1)/4/Pi)*Wigner3j(l1,l,l2,-m1,m,m2)*Wigner3j(l1,l,l2,0,0,0)): #hfl: MatElemY MatElemY:=proc(l1,m1,f,l2,m2,convention::string:=conventionY,$) simplify(int(int( Ylm(l1,m1,theta,`if`(convention="real",phi,-phi),convention,'normalized') *f(theta,phi) *Ylm(l2,m2,theta,phi,convention,'normalized'),phi=0..2*Pi)*sin(theta),theta=0..Pi)) end: ################################################################################ #cat: Other special functions #hfl: LerchPhi2 LerchPhi2:=proc(z,a,v,vmax::numeric:=10,$) local n; `if`(type(v,numeric) and v>vmax and z<>0 and z<>1, (-ln(z))^(a-1)/z^v*GAMMA(1-a,-v*ln(z))+add(pochhammer(a,n)/n!/v^n*polylogR(z,n),n=0..polylogR_n)/v^a, LerchPhi(z,a,v)) end: #hfl: LerchPhi2 polylogR:=proc(z,n) local v,k; if (z=0) then `if`(n=0,1,0) elif (n>polylogR_n) then undefined elif (n=0) then `if`(z=1,1/2,1/ln(z)+1/(1-z)) elif (n<=10 and z<0.5 or n<=20 and z<0.1) then n!/ln(z)^(n+1)+(-1)^n*polylog(-n,z) else v:=0; for k from polylogR_N by -1 to 0 do v:=polylogR_T[n,k]+v*(1-z) end; v end end: #hfl: LerchPhi2 polylogR_init:=proc(n::posint,N::posint,$) local x,e,k,n2; e:=convert(series(1/ln(1-x)+1/x,x,N+n+3),polynom); polylogR_T:=Array(0..N,0..N+n,datatype=rational); for k from 0 to N+n do polylogR_T[0,k]:=coeff(e,x,k) end; for n2 from 1 to n do for k from 0 to N+n-n2 do polylogR_T[n2,k]:=(k+1)*polylogR_T[n2-1,k+1]-k*polylogR_T[n2-1,k] end end; polylogR_n,polylogR_N:=n,N; NULL end: ################################################################################ #cat: File tools #hfl: ExpandPath ExpandPath:=proc( path::string, output::string:="apnx", { replacedirsepto::string:="/" },$) local i,p,nx,n,x,c,a,u,v; i:=max(StringTools[FirstFromRight]("/",path),StringTools[FirstFromRight]("\\",path)); p:=path[1..i]; nx:=path[i+1..-1]; i:=StringTools[FirstFromRight](".",nx); if (i=0) then n,x:=nx,"" else n,x:=nx[1..i-1],nx[i..-1] end; c:=cat(currentdir(),kernelopts(dirsep)); if (replacedirsepto<>"") then p:=StringTools[SubstituteAll](StringTools[SubstituteAll](p,"/",replacedirsepto),"\\",replacedirsepto); c:=StringTools[SubstituteAll](StringTools[SubstituteAll](c,"/",replacedirsepto),"\\",replacedirsepto) end; a:=`if`(p="/" or StringTools[Has](p,":"),"",c); seq(cat(seq(piecewise(u="a",a,u="c",c,u="p",p,u="n",n,u="x",x,""),u=v)),v=StringTools[Split](output,",")) end: #hfl: ExpandPath SimplifyPath:=proc( path::string, { replacedirsepto::string:="/" },$) local ls,exitdo,j,i,n,v; ls:=StringTools[Split](path,"/\\"); ls:=[seq(`if`(v="",NULL,v),v=ls[..-2]),ls[-1]]; do exitdo:=true; for j from nops(ls) by -1 to 1 do if (ls[j]="..") then for i from j by -1 to 1 while (ls[i]="..") do end; if (i=0) then break end; n:=min(j-i,i); ls:=[op(1..i-n,ls),op(i+1+n..-1,ls)]; exitdo:=false; break end end; if exitdo then break end end; StringTools[Join](ls,`if`(replacedirsepto="",kernelopts('dirsep'),replacedirsepto)) end: #hfl: FoldersTree FoldersTree:=proc(folder::string,depth::{infinity,nonnegint}:=infinity,$) local path; path:=cat(folder,kernelopts('dirsep')); [FileTools[Filename](folder),`if`(depth=0,NULL,op(map(FoldersTree,select(proc(v) local b; try b:=FileTools[IsDirectory](v) catch: WARNING("IsDirectory is failed on %1",v); b:=false end; b end, map2(cat,path,FileTools[ListDirectory](folder))),depth-1)))] end: #hfl: LoadTextFile LoadTextFile:=proc(f::string,keyfun::procedure:=(s->true),nlines::nonnegint:=infinity,{reload::boolean:=false},$) local fd,tb,s,i,size,v; if FileTools[IsOpen](f) then error("Cannot load open file: %1",f) end; if not(FileTools[Exists](f)) then error("File does not exists: %1",f) end; if assigned('TextFiles[f]') then if (reload or FileTools[ModificationTime](f)<>TextFilesAttr[f][1]) then UnloadTextFile(f) else return end end; fd:=fopen(f,READ,TEXT); tb:=table(); s:=readline(fd); for i from 1 to infinity while not(s=0 or keyfun(s)) do tb[i]:=s; s:=readline(fd) end; for i from i to i-1+nlines while not(s=0) do tb[i]:=s; s:=readline(fd) end; fclose(fd); size:=add(length(s),s=entries(tb,'nolist')); if (size>sTextFiles) then error("File size =%1 and is larger than sTextFiles=%2",size,sTextFiles) end; TextFiles[f]:=convert(tb,list); TextFilesAttr[f]:=[ FileTools[ModificationTime](f), size ]; # [modtime,size] TextFilesList:=[op(TextFilesList),f]; # Cleanup while (nops(TextFilesList)>nTextFiles) do UnloadTextFile(TextFilesList[1]) end; while (add(TextFilesAttr[v][2],v=TextFilesList)>sTextFiles) do for i from nops(TextFilesList)-1 by -1 to 2 while (TextFilesAttr[TextFilesList[i]][2]<=add(TextFilesAttr[v][2],v=TextFilesList[..i-1])) do end; UnloadTextFile(TextFilesList[i]) end; NULL end: #hfl: LoadTextFile UnloadTextFile:=proc(f::string,$) if assigned('TextFiles[f]') then unassign('TextFiles[f],TextFilesAttr[f]'); TextFilesList:=Substitute(TextFilesList,f) end; NULL end: #hfl: SearchFilePos SearchFilePos:=proc( f::string, key::string, keypos::range:=1..-1, number::integer:=1, { skiplines::integer:=0, linematch::boolean:=`if`(key="",true,false), reload::boolean:=false, noerror::boolean:=false },$) local j0,L,n,i1,tb,k,i,lsj,j; j0:=op(1,keypos)-1; if (j0<0) then error("No negative lower bound in keypos: %1",keypos) end; LoadTextFile(f,':-reload'=reload); L:=TextFiles[f]; n:=nops(L); i1:=`if`(skiplines<0,max(1,n+skiplines+1),skiplines+1); tb:=table(); k:=0; if linematch then if (number>0) then for i from i1 to n do if (L[i]=key) then k:=k+1; if (k=number) then return [i,1] end end end elif (number<0) then for i from n-i1+1 by -1 to 1 do if (L[i]=key) then k:=k-1; if (k=number) then return [i,1] end end end else return [seq(`if`(L[i]=key,[i,1],NULL),i=i1..n)] end else if (number>0) then for i from i1 to n do lsj:=[StringTools[SearchAll](key,L[i][keypos])]; if (lsj<>[]) then if (nops(lsj)[]) then if (nops(lsj)evalb(s="")), { skiplines::integer:=0, linematch::boolean:=`if`(key="",true,false), shift::integer:=`if`(key=undefined,0,1), nlines::{nonnegint,infinity}:=`if`(key=undefined,10^9,infinity), line::name:=undefined, reload::boolean:=false },$) local i0,p,i2,L,i; if (key=undefined) then LoadTextFile(f,':-reload'=reload); i0:=1 else p:=SearchFilePos(f,key,keypos,`if`(number=0,1,number),':-skiplines'=skiplines,':-linematch'=linematch,':-reload'=reload,'noerror'); if (p=[]) then return [] else i0:=p[1] end end; i2:=min(nops(TextFiles[f]),i0+shift+nlines-1); L:=TextFiles[f][max(1,i0+shift)..i2]; if (nlines=infinity) then for i from 1 to nops(L) do if endkeyfun(L[i]) then i2:=i2+i-nops(L); L:=L[..i-1]; break end end end; if not(line=undefined) then assign(line,i2) end; L end: #hfl: ReadValue ReadValue:=proc( f::string, key::{string,undefined}:=undefined, keypos::range:=1..-1, number::integer:=1, { skiplines::integer:=0, linematch::boolean:=`if`(key="",true,false), shift::integer:=0, startofline::boolean:=false, pos::range:=1..-1, format::string:="", line::name:=undefined, reload::boolean:=false, raiseerror::boolean:=false },$) local p,n,i,j,s,v; if (key=undefined) then LoadTextFile(f,':-reload'=reload); p:=[1,1] else p:=SearchFilePos(f,key,keypos,`if`(number=0,1,number),':-skiplines'=skiplines,':-linematch'=linematch,':-reload'=reload,'noerror'=not(raiseerror)); if (p=[]) then return end end; n:=nops(TextFiles[f]); i:=p[1]+shift; if (i<1 or i>n) then return end; if not(line=undefined) then assign(line,i) end; j:=`if`(startofline or shift<>0,1,p[2]+`if`(key=undefined,0,length(key))); s:=TextFiles[f][i][j..][pos]; if (format="") then s else v:=sscanf(s,format); if (v=0 or v=[]) then if raiseerror then error("Cannot read format=%1 in %2",format,s) else s end else op(v) end end end: #hfl: ReadPaginatedMatrix ReadPaginatedMatrix:=proc( f::string, key::{string,undefined}:=undefined, keypos::range:=1..-1, number::integer:=1, endkeyfun::procedure:=(s->evalb(s="")), { skiplines::integer:=0, linematch::boolean:=`if`(key="",true,false), shift::integer:=`if`(key=undefined,0,1), pos::range:=1..-1, format::string:="%f", sep::string:=" ", sym::boolean:=false, skip::nonnegint:=0, reload::boolean:=false },$) local p,L,k,n,V,m,j,i,M,ls; if (key=undefined) then LoadTextFile(f,':-reload'=reload); p:=[1,1] else p:=SearchFilePos(f,key,keypos,`if`(number=0,1,number),':-skiplines'=skiplines,':-linematch'=linematch,':-reload'=reload) end; L:=TextFiles[f][max(1,p[1]+shift)..]; for k from 1 to nops(L) while not(endkeyfun(L[k])) do end; n:=k-1; V:=Vector(n,i->L[i]); m:=nops(remove(`=`,StringTools[Split](V[n],sep),"")); if (m=0) then error("No entries in %1",V[n]) end; for j from m+1 by m to n do for k from k to nops(L) while endkeyfun(L[k]) do end; k:=k-1+skip; for i from `if`(sym,j,1) to n do k:=k+1; V[i]:=cat(V[i],sep,L[k][pos]) end end; M:=Matrix(n,`if`(sym,shape=symmetric,NULL),`if`(format="%f",datatype=float,`if`(format="%d",datatype=integer,NULL))); for i from 1 to n do ls:=remove(`=`,StringTools[Split](V[i],sep),""); if (nops(ls)=n or sym and nops(ls)=i) then for j from 1 to `if`(sym,i,n) do M[i,j]:=op(sscanf(ls[j],format)) end else error("Wrong number of entries in %1",V[i]) end end; M end: #hfl: WriteLines WriteLines:=proc(f::{string,integer},ls::list,{overwrite::boolean:=false,append::boolean:=false,printout::{boolean,nonnegint,[nonnegint,nonnegint]}:=false},$) local fd,s; if type(f,integer) then fd:=f elif (f="") then fd:='default' else if FileTools[Exists](f) then if (overwrite or append) then fd:=fopen(f,`if`(append,'APPEND','WRITE'),'TEXT') else error("File exists: %1",f) end else fd:=fopen(f,'WRITE','TEXT') end end; for s in ls do if type(s,string) then writeline(fd,s) end end; if (type(f,string) and f<>"") then fclose(fd); if not(printout=false) then PrintFile(f,`if`(printout=true,NULL,`if`(type(printout,posint),printout,op(printout)))) end end; NULL end: #hfl: ReplaceLines ReplaceLines:=proc(f::string,change::list,{nobackup::boolean:=false,forcebackup::boolean:=false},$) local ls; ls:=ReadLines(f); ls:=map(proc(s) local v; for v in change do if SearchText(v[1],s)>0 then return `if`(nops(v)=1,NULL,`if`(nops(v)=2,v[2],StringTools[Substitute](s,v[2],v[3]))) end end; s end,ls); FileTools[Rename](f,cat(f,".bak"),'force'=forcebackup); WriteLines(f,ls,'overwrite'); if nobackup then FileTools[Remove](cat(f,".bak")) end; NULL end: #hfl: ReadRecord ReadRecord:=proc( filename::string, records::string:="", linenum::nonnegint:=0, opt4ReadLines::list:=[], { comment::string:="#", assignment::string:="=", nlines::posint:=1, input::{"file","string","auto"}:="auto", nowarning::boolean:=false },$) local ls,tb,s,l,s1,i,out,dv,t,v; if (input="file" or (input="auto" and FileTools[Exists](filename))) then if (linenum=0) then ls:=remove(s->s[1]=comment,map(Trim,ReadLines(filename,op(opt4ReadLines)))) else ls:=ReadLines(filename,'shift'=linenum-1,':-nlines'=nlines); if (ls<>[]) then ls:=select(s->SearchText(assignment,s)>0,map(Trim,map(op@StringTools[Split],ls),[",",".",";",":","?","!"])) end end elif (input="string" or (input="auto" and SearchText(assignment,filename)>0)) then ls:=StringTools[Split](filename,"\n"); ls:=select(s->SearchText(assignment,s)>0,map(Trim,map(op@StringTools[Split],ls),[",",".",";",":","?","!"])) else if (input="auto" and not(nowarning)) then WARNING("File %1 does not exist and filename does not contain assignment symbol %2",filename,assignment) end; ls:=[] end; tb:=table(); s:=""; for l from nops(ls) by -1 to 1 do s1:=ls[l]; i:=SearchText(comment,s1); if (i>0) then s1:=TrimRight(s1[..i-1]) end; i:=SearchText(assignment,s1); if (i>0) then if (s1[i+1]=assignment) then s:=cat(s1[..i],s1[i+2..],s) else tb[TrimRight(s1[..i-1])]:=cat(TrimLeft(s1[i+1..]),s); s:="" end else s:=cat(s1,s) end end; if (s<>"") then error("Unprocessed string in %1, s=%2",filename,s) end; if (records="") then return op(tb) end; out:=[]; for s in StringTools[Split](records,",") do i:=SearchText(":=",s); if (i>0) then s,dv:=s[..i-1],s[i+2..] else s,dv:=s,"" end; i:=SearchText("::",s); if (i>0) then t:=parse(s[i+2..]); if not(type(t,type)) then error("Unrecognized type: %1",t) end; s:=s[..i-1] else s,t:=s,anything end; if assigned('tb[s]') then if (t=list(string)) then if not(tb[s][1]="[" and tb[s][-1]="]") then error("Unrecognized list(string) record: %1",tb[s]) else v:=map(Trim,StringTools[Split](tb[s][2..-2],",")) end else v:=`if`(t=string,tb[s],eval(parse(tb[s]))) end; if not(type(v,t) or t=table and type(v,list)) then error("Unrecognized value (must be of type %1): %2",t,v) end else if (dv="") then v:=`if`(t=string,"",`if`(t=list or t=list(string) or t=table,[],`if`(subtype(t,extended_numeric),0,undefined))) else v:=`if`(t=string,dv,eval(parse(dv))) end end; # out:=out,`if`(t=table,table(v),v) end; out:=[op(out),[t,v]] end; #`if`(nops([out])=1 and type(out,table),op(out),out) seq(`if`(v[1]=table,table(v[2]),v[2]),v=out) end: #hfl: ReadRecord WriteRecord:=proc(filename::string,records::list({string,[string,anything],[string,anything,string]}),{overwrite::boolean:=false,append::boolean:=false},$) local printdatatype,formatrecord; printdatatype:=v->`if`(rtable_options(v,datatype)=anything,"",sprintf(",datatype=%a",rtable_options(v,datatype))); formatrecord:=proc(v,$) local s; if type(v,string) then v elif type(v,[string,anything]) then `if`(type(v[2],list(string)), cat(v[1],"=[",StringTools[Join](v[2],","),"]"), sprintf("%s=%a",op(v))) else if type(v[2],atomic) then sprintf(cat("%s=",v[3]),v[1],v[2]) elif type(v[2],list(list)) then sprintf("%s=[%{c,}s]",v[1],Vector(nops(v[2]),i->sprintf(cat("[%{c,}",v[3][2..],"]"),Vector(nops(v[2][i]),j->v[2][i][j])))) elif type(v[2],list) then sprintf(cat("%s=[%{c,}",v[3][2..],"]"),v[1],Vector(nops(v[2]),i->v[2][i])) elif type(v[2],Vector) then sprintf(cat("%s=Vector(%d%s,[%{c,}",v[3][2..],"])"),v[1],op(1,v[2]),printdatatype(v[2]),v[2]) elif type(v[2],Matrix) then s:=StringTools[SubstituteAll](sprintf(cat("%{c,}",v[3][2..]),v[2]),"\n",",\n"); sprintf("%s=Matrix(%d,%d%s,[\n%s])",v[1],op(1,v[2]),printdatatype(v[2]),s) elif type(v[2],table) then s:=StringTools[Join](map(StringTools[Substitute],FormatTable(v[2],v[3],'align'),"=","=="),",\n"); sprintf("%s=[\n%s,\nNULL]",v[1],s) else error("Unsupported type in %1",v) end end end; WriteLines(filename,map(formatrecord,records),':-overwrite'=overwrite,':-append'=append) end: #hfl: ReadSet ReadSet:=proc(s::string,f::string,n::integer:=99,$) local ls,ls2,s2; ls:=StringTools[Split](s,","); if (n<0) then [] else if (nops(ls)=1) then ls:=ReadRecord(f,cat(s,"::list(string)"),'input'="file"); if (ls=[]) then [s] else ReadSet(StringTools[Join](ls,","),f,n-1) end else ListTools[MakeUnique]([seq(op(ReadSet(s2,f,n)),s2=ls)]) end end end: #hfl: ReadUTF8 ReadUTF8:=proc(fn::string,{unicodes::boolean:=false},$) local ls,tb,i,k; ls:=FileTools[Binary][Read](fn,integer[1],'byteorder'='native'); fclose(fn); ls:=map(v->`if`(v>0,v,v+256),ls); tb:=table(); for i from 1 to nops(ls) do if (ls[i]<128) then tb[i]:=ls[i] elif (ls[i]<192) then error("Unrecognized byte #%1: %2",i,ls[i]) elif (ls[i]<224) then tb[i]:=64*(ls[i]-192)+ls[i+1]-128; i:=i+1 elif (ls[i]<240) then tb[i]:=64^2*(ls[i]-224)+add(64^(2-k)*(ls[i+k]-128),k=1..2); i:=i+2 else tb[i]:=64^3*(ls[i]-240)+add(64^(3-k)*(ls[i+k]-128),k=1..3); i:=i+3 end end; ls:=convert(tb,list); if unicodes then ls else cat(seq(`if`(k<256,StringTools[Char](k),sprintf("&#%04d;",k)),k=ls)) end end: #hfl: ReadUTF8 WriteUTF8:=proc(fn::string,s::{string,list(nonnegint)},{overwrite::boolean:=false},$) local ls,tb,i,q,k; if (not(overwrite) and FileTools[Exists](fn)) then error("File exists: #1",fn) end; if type(s,list) then ls:=s else tb:=table(); for i from 1 to length(s) do if (s[i]="&" and s[i+1]="#" and s[i+6]=";") then tb[i]:=parse(s[i+2..i+5]); if type(tb[i],nonnegint) then i:=i+6 else tb[i]:=38 end else tb[i]:=StringTools[Ord](s[i]) end end; ls:=convert(tb,list) end; tb:=table(); for i from 1 to nops(ls) do if (ls[i]<128 ) then tb[4*i]:=ls[i] elif (ls[i]<2048 ) then tb[4*i]:=128+irem(ls[i],64,'q'); tb[4*i-1]:=192+q elif (ls[i]<65536) then tb[4*i]:=128+irem(ls[i],64,'q'); tb[4*i-1]:=128+irem(q,64,'q'); tb[4*i-2]:=224+q else tb[4*i]:=128+irem(ls[i],64,'q'); for k from 1 to 2 do tb[4*i-k]:=128+irem(q,64,'q') end; tb[4*i-3]:=240+q end end; ls:=Vector(map(v->`if`(v<128,v,v-256),convert(tb,list)),datatype=integer[1]); FileTools[Binary][Write](fn,ls,'byteorder'='native'); fclose(fn) end: #hfl: ReadBIN ReadBIN:=proc( filename::string, { code::[nonnegint,nonnegint,nonnegint,nonnegint]:=[0,0,0,0], forcecode::boolean:=false, dimensions::list:=[], datapos::integer:=-1, nodescription::boolean:=false, raw::{boolean,nonnegint}:=false, RMO::boolean:=false, defelementsize::boolean:=true, optrtable::list:=[], buffersize::posint:=10000, printout::boolean:=false },$) local fd,cod,pos,id,storedcode,i,dt,bo,dims,ls,j,dim2,dts,size,Vl,tb,V,p,k1,N,R,i1,i2,i3,i4, v,l,k; fd:=fopen(filename,READ,BINARY); try # code and pos if nodescription then id:=undefined; if (datapos<0) then error("Datapos<0, %1",pos) end; cod,pos:=code,datapos; if (dimensions=[]) then if (cod[3]=0) then cod[3]:=1 elif not(member(cod[3],[1,11,12])) then error("Dimensions are required for multidimensional arrays, code[2]=%1",code[2]) end end else id:=op(FileTools[Binary][Read](fd,integer[4],1,byteorder='native')); if not(member(id,binformatidentifiers1)) then error("Wrong file identifier") end; pos:=`if`(datapos<0,32,datapos); storedcode:=readbytes(fd,4); cod:=storedcode; for i from 1 to 4 do if (code[i]>0 and code[i]<>storedcode[i]) then cod[i]:=code[i] else cod[i]:=storedcode[i] end end; if (cod<>storedcode) then if forcecode then WARNING("Format code is changed from %1 to %2",storedcode,cod) else error("Stored format code %1<>%2 provided code",storedcode,cod) end end end; if not(member(cod[1],[1,2]) and member(cod[2],[1,2,4,8]) and member(cod[4],[1,2,3,4])) then error("Unrecognized format code, %1",cod) end; # datatype, byteorder, and dimensions dt:=[integer,float][cod[1]][cod[2]]; bo:=piecewise(cod[4]=1,'native',cod[4]=2,'network',cod[4]=3,'big',cod[4]=4,'little',NULL); dims:=`if`(dimensions=[], `if`(nodescription, [(FileTools[Size](fd)-pos)/cod[2]], FileTools[Binary][Read](fd,integer[4], `if`(member(cod[3],[1,11,12,13,14,15,3,31,32]),1, `if`(member(cod[3],[2,21,22,53]),2, `if`(cod[3]=54,3, cod[3]-40))),'byteorder'=bo)), dimensions); # special bytes if (cod[3]=13) then ls:=readbytes(fd,20); for j from 1 to nops(ls) while (ls[j]>0) do end; dim2:=j-1; dts:=[seq([integer,float][1+iquo(v,128)][irem(v,128)],v=ls[..dim2])] end; # read data size:=`if`(member(cod[3],[3,31,32,53,54]),dims[1]*(dims[1]+1)/2*mul(v,v=dims[2..]),mul(v,v=dims)); if printout then printf("id=%d, code=[%d,%d,%d,%d], pos=%d, size=%d\n",id,op(cod),pos,size) end; if type(raw,nonnegint) then size:=min(size,raw) end; filepos(fd,pos); if (cod[3]=13) then Vl:=FileTools[Binary][Read](fd,dt,dims[1],'byteorder'=bo); size:=add(l,l=Vl); tb:=table(); for j from 1 to dim2 do V:=Vector(size,datatype=dts[j]); for i from 0 by buffersize to size-1 do ls:=FileTools[Binary][Read](fd,dts[j],min(buffersize,size-i),'byteorder'=bo); V[i+1..i+nops(ls)]:=Vector(ls,datatype=dts[j]) end; tb[j]:=copy(V) end else V:=Vector(size,datatype=dt); for i from 0 by buffersize to size-1 do ls:=FileTools[Binary][Read](fd,dt,min(buffersize,size-i),'byteorder'=bo); V[i+1..i+nops(ls)]:=Vector(ls,datatype=dt) end end; catch: error finally fclose(fd) end try; # compose rtable if defelementsize then dt:=op(0,dt) end; if not(raw=false) then [id,cod,dims,`if`(size>buffersize,V,convert(V,list))] elif (cod[3]=1 ) then Vector (dims[1],V,datatype=dt,op(optrtable)) elif (cod[3]=11) then Vector[column](dims[1],V,datatype=dt,op(optrtable)) elif (cod[3]=12) then Vector[row] (dims[1],V,datatype=dt,op(optrtable)) elif (cod[3]=13) then if (dim2=1) then V:=Vector(dims[1]); k1:=0; for i from 1 to dims[1] do V[i]:=[seq(tb[1][k],k=k1+1..k1+Vl[i])]; k1:=k1+Vl[i] end; V else V:=Vector(dims[1]); k1:=0; for i from 1 to dims[1] do V[i]:=[seq([seq(tb[j][k],j=1..dim2)],k=k1+1..k1+Vl[i])]; k1:=k1+Vl[i] end; V end elif (cod[3]=2 or cod[3]=21) then Matrix(op(dims),(i,j)->V[j+dims[2]*(i-1)],datatype=dt,op(optrtable)) elif (cod[3]=22) then Matrix(op(dims),(i,j)->V[i+dims[1]*(j-1)],datatype=dt,op(optrtable)) elif (cod[3]=3 or cod[3]=31) then Matrix(dims[1],(i,j)->V[i*(i-1)/2+j],datatype=dt,shape=symmetric,storage=triangular[lower],op(optrtable)) elif (cod[3]=32) then Matrix(dims[1],(i,j)->V[j*(j-1)/2+i],datatype=dt,shape=symmetric,storage=triangular[lower],op(optrtable)) elif (cod[3]=41) then Array(1..dims[1],(i1)->V[i1],datatype=dt,op(optrtable)) elif (cod[3]=42) then if RMO then Array(1..dims[1],1..dims[2],(i1,i2)->V[i2+dims[2]*(i1-1)],datatype=dt,op(optrtable)) else Array(1..dims[1],1..dims[2],(i1,i2)->V[i1+dims[1]*(i2-1)],datatype=dt,op(optrtable)) end elif (cod[3]=43) then if RMO then Array(1..dims[1],1..dims[2],1..dims[3],(i1,i2,i3)->V[i3+dims[3]*(i2-1+dims[2]*(i1-1))],datatype=dt,op(optrtable)) else Array(1..dims[1],1..dims[2],1..dims[3],(i1,i2,i3)->V[i1+dims[1]*(i2-1+dims[2]*(i3-1))],datatype=dt,op(optrtable)) end elif (cod[3]=44) then if RMO then Array(1..dims[1],1..dims[2],1..dims[3],1..dims[4],(i1,i2,i3,i4)->V[i4+dims[4]*(i3-1+dims[3]*(i2-1+dims[2]*(i1-1)))],datatype=dt,op(optrtable)) else Array(1..dims[1],1..dims[2],1..dims[3],1..dims[4],(i1,i2,i3,i4)->V[i1+dims[1]*(i2-1+dims[2]*(i3-1+dims[3]*(i4-1)))],datatype=dt,op(optrtable)) end elif (cod[3]=53) then N:=dims[1]*(dims[1]+1)/2; R:=Array(symmetric12,1..dims[1],1..dims[1],1..dims[2],datatype=dt,op(optrtable)); for i1 from 1 to dims[1] do for i2 from 1 to i1 do for i3 from 1 to dims[2] do R[i1,i2,i3]:=V[i1*(i1-1)/2+i2+N*(i3-1)] end end end; R elif (cod[3]=54) then N:=dims[1]*(dims[1]+1)/2; R:=Array(symmetric12,1..dims[1],1..dims[1],1..dims[2],1..dims[3],datatype=dt,op(optrtable)); for i1 from 1 to dims[1] do for i2 from 1 to i1 do for i3 from 1 to dims[2] do for i4 from 1 to dims[3] do R[i1,i2,i3,i4]:=V[i1*(i1-1)/2+i2+N*(i3-1+dims[3]*(i4-1))] end end end end; R else error("Unrecognized shape, %1",cod[3]) end end: #hfl: ReadBIN WriteBIN:=proc( filename::string, R::{list,Vector,Matrix,Array}, opt::list:=[], { code::[nonnegint,nonnegint,nonnegint,nonnegint]:=[0,0,0,1], dimensions::list:=[], nodescription::boolean:=false, overwrite::boolean:=false },$) local cod,dt,v,dims,dim2,dts,bo,fd,pos,j, u,i; # determine datatype cod:=code; if (cod[1]*cod[2]>0) then dt:=[integer,float][cod[1]][cod[2]] else dt:=rtable_options(R,'datatype'); v:=op(0,dt); if (v='integer') then cod[1]:=1; cod[2]:=op(1,dt) elif (v='float') then cod[1]:=2; cod[2]:=op(1,dt) elif type(R[1],list) then cod[1]:=1; cod[2]:=4; dt:=integer[4] else error("Unsupported element type %1",v) end end; # determine shape and dimensions if (cod[3]>0 and dimensions<>[]) then dims:=dimensions else if type(R,Vector) then dims:=[op(1,R)]; cod[3]:=`if`(type(R[1],list),13,1) elif type(R,Matrix) then v:=MatrixOptions(R,'shape'); if (v=[]) then cod[3],dims:=2,[op(1,R)] elif (v=['symmetric']) then cod[3],dims:=3,[[op(1,R)][1]] else error("Unsupported shape %1",v) end elif type(R,Array) then dims:=[rtable_dims(R)]; if (nops(dims)>4) then error("Too large number of dimensions: %1",nops(dims)) end; if not(foldl(`and`,seq(evalb(op(1,v)=1),v=dims))) then error("Each dimension of Array must be enumeratted from 1, but received %1",dims) end; cod[3]:=40+nops(dims); dims:=map2(op,2,dims) end end; # process additional arguments if (cod[3]=13) then dim2:=`if`(type(R[1][1],list),nops(R[1][1]),1); if (dim2>20) then error("Too large secondary dimension (max 20), %1",dim2) end; if (nops(opt)=1) then dts:=[opt[1]$dim2] elif (nops(opt)=dim2) then dts:=opt else error("Unrecognized datatypes: %1",opt) end; dts:=map(v->`if`(v='integer',integer[4],`if`(v='float',float[8],v)),dts) end; # prepare writing bo:=['native','network','big','little'][cod[4]]; if (FileTools[Exists](filename) and not(overwrite)) then error("File exists, %1",filename) end; fd:=fopen(filename,WRITE,BINARY); # write description if not(nodescription) then FileTools[Binary][Write](fd,integer[4],[1760568055],'byteorder'=bo); writebytes(fd,cod); FileTools[Binary][Write](fd,integer[4],dims,'byteorder'=bo); pos:=4+4+4*nops(dims); if (cod[3]=13) then writebytes(fd,map(v->`if`(op(0,v)='float',128,0)+op(1,v),dts)); pos:=pos+dim2 end; writebytes(fd,[0$(32-pos)]) end; # write data if type(R,list) then FileTools[Binary][Write](fd,dt,R,'byteorder'=bo) elif (cod[3]=1) then FileTools[Binary][Write](fd, `if`(cod[1]*cod[2]>0, Vector(R,'datatype'=dt), R), 'byteorder'=bo) elif (cod[3]=13) then FileTools[Binary][Write](fd,Vector(dims[1],i->nops(R[i]),'datatype'=dt),'byteorder'=bo); if (dim2=1) then FileTools[Binary][Write](fd,dts[1],[seq(op(v),v=R)],'byteorder'=bo) else for j from 1 to dim2 do FileTools[Binary][Write](fd,dts[j],[seq(seq(u[j],u=v),v=R)],'byteorder'=bo) end end elif (cod[3]=2) then FileTools[Binary][Write](fd, `if`(cod[1]*cod[2]>0, Matrix(LinearAlgebra[Transpose](R),'datatype'=dt), LinearAlgebra[Transpose](R)), 'byteorder'=bo) elif (cod[3]=3) then FileTools[Binary][Write](fd, Vector([seq(seq(R[i,j],i=1..j),j=1..dims[1])],'datatype'=dt), 'byteorder'=bo) else FileTools[Binary][Write](fd,R,'byteorder'=bo) end; fclose(fd); FileTools[Size](filename) end: #hfl: PrintFile PrintFile:=proc(filename::{string,list(string)},nlines::nonnegint:=0,nendlines::nonnegint:=0,$) local ls,s,i; if (nlines=0 and nendlines=0) then if type(filename,string) then printf("%s\n",FileTools[Text][ReadFile](filename)) else seq(printf("%s\n",s),s=filename) end elif (nendlines=0) then ls:=`if`(type(filename,string),ReadLines(filename,':-nlines'='nlines'),filename[..min(nlines,nops(filename))]); for s in ls do printf("%s\n",s) end; if (nops(ls)=nlines) then printf("...\n") end else ls:=`if`(type(filename,string),StringTools[Split](FileTools[Text][ReadFile](filename),"\n"),filename); if (nlines+nendlinesf2 and not(overwrite) and FileTools[Exists](f2)) then error("File exists, %1",f2) end; FileTools[Binary][Open](f1); ls:=FileTools[Binary][Read](f1,integer[1]); FileTools[Binary][Close](f1); tb:=table(): for i from 1 to nops(ls) do if (ls[i]=13 and if2 and not(overwrite) and FileTools[Exists](f2)) then error("File exists, %1",f2) end; if force then WriteLines(f2,ReadLines(f1),':-overwrite'=true); true else fd:=fopen(f1,READ,TEXT); readline(fd); i:=filepos(fd); s:=readline(fd); try filepos(fd,i); s2:=readline(fd) catch: end; fclose(fd); if (s<>s2) then WriteLines(f2,ReadLines(f1),':-overwrite'=true) end; evalb(s<>s2) end end: #hfl: OpenFolder OpenFolder:=proc(filename::string,{noerror::boolean:=false},$) local s; s:=SimplifyPath(ExpandPath(filename),'replacedirsepto'=kernelopts('dirsep')); if not(FileTools[Exists](s)) then if noerror then ssystem(cat("explorer /select,",s)) else error("File or folder does not exist: %1",s) end else if FileTools[IsDirectory](s) then ssystem(cat("explorer /root,",s)) else ssystem(cat("explorer /select,",s)) end end; NULL end: ################################################################################ #cat: Data manipulation #hfl: RootMeanSquare RootMeanSquare:=proc(V::indexable,$) local n,a,v; n:=add(1,v=V); a:=add(v,v=V)/n; sqrt(add((v-a)^2,v=V)/n) end: #hfl: LinearFit2 LinearFit2:=proc( XorXY::{list,Vector,Matrix}, Yor::{list,Vector}:=[], o::posint:=1, ex::list(nonnegint):=[], output::string:="", {printout::boolean:=false, format::string:="%.3g"}, $) local X,Y,lsp,n,tb,x,p,d,std,Y0,v,r2,c,cstd,out,lo,e; if (Yor=[]) then if type(XorXY,list) then X,Y:=map2(op,1,XorXY),map2(op,2,XorXY) else X,Y:=XorXY[1,..],XorXY[2,..] end else X,Y:=XorXY,Yor end; lsp:=remove(member,[$0..o],ex); n:=nops(lsp); tb:=table(Statistics[LinearFit]([seq(x^p,p=lsp)],X,Y,x,':-output'='solutionmodule'):-Results()); d:=tb["degreesoffreedom"]; std:=`if`(assigned('tb["residualstandarddeviation"]'),tb["residualstandarddeviation"],undefined); Y0:=Statistics[Mean](Y); r2:=`if`(assigned('tb["residuals"]'),1-add(v^2,v=tb["residuals"])/add((v-Y0)^2,v=Y),undefined); c:=convert(tb["parametervalues"],list); cstd:=`if`(d=0,[0$n],convert(tb["standarderrors"],list)); if printout then printf("n=%d, o=%d, d=%d, std=%s, r2=%s\n",n,o,d,Trim(sprintf(format,std)),Trim(sprintf(format,r2))); printf("c =[%{c,}s]\ncstd=[%{c,}s]\n",Vector(map2(sprintf,format,c)),Vector(map2(sprintf,format,cstd))) end; out:=StringTools[Split](output,","); lo:=[]: if (out=[""]) then op(tb) else for e in out do if (e="f" ) then lo:=[op(lo),unapply(tb["leastsquaresfunction"],x)] elif (e="std" ) then lo:=[op(lo),std] elif (e="r2" ) then lo:=[op(lo),r2] elif (e="c" ) then lo:=[op(lo),c] elif (e="cstd") then lo:=[op(lo),cstd] elif (e="cint") then lo:=[op(lo),convert(tb["confidenceintervals"],list)] elif (e="res" ) then lo:=[op(lo),`if`(d=0,map(v->0,Y),tb["residuals"])] elif (e="cov" ) then lo:=[op(lo),`if`(d=0,Matrix(n),tb["variancecovariancematrix"])] end end; op(lo) end end: #hfl: LinearFitTrig LinearFitTrig:=proc( m1::nonnegint, m2::nonnegint, P::numeric, XY::list, c1::numeric:=1, XY1::list:=[], c2::numeric:=1, XY2::list:=[], output::string:="c,f,std", { weights::{list,procedure,numeric}:=((x,y,p)->1),maxcondnum::numeric:=1e-4, digits::[posint,nonnegint]:=[6,3], printout::boolean:=false },$) local cc,phi,bas,bas1,bas2,nb,n,n1,n2,m,np,W,B,b,i,XXi,XX,CN,Y,C,R,RMS,RMS0,RMS1,RMS2,f,out, v,s; cc:=evalf(2*Pi/P); bas :=[seq( cos(m*phi),m=0..m1),seq( sin(m*phi),m=1..m2)]; bas1:=[seq(-m *sin(m*phi),m=0..m1),seq( m *cos(m*phi),m=1..m2)]; bas2:=[seq(-m^2*cos(m*phi),m=0..m1),seq(-m^2*sin(m*phi),m=1..m2)]; nb:=nops(bas); n,n1,n2:=nops(XY),nops(XY1),nops(XY2); m:=n+n1; np:=m+n2; if type(weights,list) then if (nops(weights)<>np) then error("nops(weights)=%1<>np=%2",nops(weights),np) end; W:=weights else if type(weights,numeric) then if (np<>n) then error("Boltzmann weights are not supported for derivatives") else W:=[seq(evalf(exp(-weights*v[2])),v=XY)] end else W:=[seq(evalf(weights(op(v),0)),v=XY),seq(evalf(weights(op(v),1)),v=XY1),seq(evalf(weights(op(v),2)),v=XY2)] end; W:=W/max(W) end; W:=LinearAlgebra[DiagonalMatrix](W,datatype=float); B:=Matrix(np,nb,datatype=float); for b from 1 to nb do for i from 1 to n do B[ i,b]:=evalf(eval( bas [b],phi=cc*XY [i][1])) end; for i from 1 to n1 do B[n+i,b]:=evalf(eval(c1*bas1[b],phi=cc*XY1[i][1])) end; for i from 1 to n2 do B[m+i,b]:=evalf(eval(c2*bas2[b],phi=cc*XY2[i][1])) end end; XXi:=Matrix(LinearAlgebra[Transpose](B).W.B,shape=symmetric); XX:=Matrix(LinearAlgebra[MatrixInverse](XXi),shape=symmetric,datatype=float); CN:=LinearAlgebra[Norm](XX.XXi-LinearAlgebra[IdentityMatrix](nb),'Frobenius'); if (CN>maxcondnum) then error("Ill conditioned basis, condnum=%1",CN) end; Y:=Vector(np,datatype=float); for i from 1 to n do Y[ i]:=XY [i][2] end; for i from 1 to n1 do Y[n+i]:=XY1[i][2]*c1/cc end; for i from 1 to n2 do Y[m+i]:=XY2[i][2]*c2/cc^2 end; C:=XX.LinearAlgebra[Transpose](B).W.Y; R:=Y-B.C; RMS:=evalf(LinearAlgebra[Norm](R,2)/sqrt(np)); RMS0:=`if`(n =0 ,undefined,evalf(LinearAlgebra[Norm](R[ ..n],2)/sqrt(n))); RMS1:=`if`(n1=0 or c1=0,undefined,evalf(LinearAlgebra[Norm](R[n+1..m],2)/sqrt(n1)/c1)); RMS2:=`if`(n2=0 or c2=0,undefined,evalf(LinearAlgebra[Norm](R[-n2.. ],2)/sqrt(n2)/c2)); if printout then printf("%*d %*d\n%*.*f\n",digits[1]-digits[2],0,digits[1],Vector(max(m1,m2),m->m),op(digits),C[..1+m1]); if (m2>0) then printf("%*s %*.*f\n",digits[1],"",op(digits),C[-m2..]) end; if (n=np) then printf("RMS=%.3g\n",RMS) else printf("RMS=%.3g (0=%.3g, 1=%.3g, 2=%.3g)\n",RMS,RMS0,RMS1,RMS2) end end; f:=unapply(subs(phi=cc*phi,add(C[b]*bas[b],b=1..nb)),phi); out:=StringTools[Split](output,","); seq(`if`(s="c",C,`if`(s="f",op(f),`if`(s="std",RMS,NULL))),s=out) end: #hfl: MultiDimFit MultiDimFit:=proc(X0::{listlist,Vector,Matrix}, Y::Vector, basis::string, order0::list(nonnegint), rng::list({numeric,numeric..numeric}):=[0$nops(order0)], output::string:="f", { cmin::numeric:=10^(2-Digits), printout::boolean:=false},$) local order,d,bf,l,m,x,nx,d2,X,a,b,xi1,xi2,x1,x2,nb,iidB,j,p,B,C,cmin2,CM, v,e; order:=order0; d:=nops(order); if (length(basis)<>d) then error("Inconsistent sizes of basis=%1 and order=%2",basis,order) end; for l from 1 to d do if (basis[l]="t" and type(order[l],odd)) then order[l]:=order[l]+1 end end; bf:=Array(1..d,0..max(order)); for l from 1 to d do for m from 0 to order[l] do bf[l,m]:=unapply( `if`(basis[l]="t",`if`(type(m,even),cos(m/2*x),sin((m+1)/2*x)), `if`(basis[l]="c",orthopoly['T'](m,x), undefined)),x) end end; if type(X0,Matrix) then nx,d2:=op(1,X0); X:=X0 else nx,d2:=Dim2(X0),Dim2(X0[1]); X:=Matrix(nx,d2,(i,l)->X0[i][l],datatype=float) end; if (d2<>d) then error("Inconsistent sizes of X and order") end; if (nops(rng)<>d) then error("Inconsistent sizes of rng and order") end; a,b:=Vector(d,datatype=float),Vector(d,datatype=float); for l from 1 to d do if (basis[l]="t") then xi1,xi2:=0,2*Pi elif (basis[l]="c") then xi1,xi2:=-1,1 else error("Unrecognized basis: %1",basis[l]) end; if (rng[l]=0) then x1,x2:=min(X[..,l]),max(X[..,l]) elif type(rng[l],numeric) then x1,x2:=0,rng[l] else x1,x2:=op(rng[l]) end; a[l],b[l]:=evalf((xi2-xi1)/(x2-x1)),evalf((x2*xi1-x1*xi2)/(x2-x1)) end; nb:=mul(v+1,v=order); iidB:=Matrix(nb,d,datatype=integer); for j from 1 to nb do p:=j-1; for l from 1 to d do iidB[j,l]:=irem(p,order[l]+1,'p') end end; B:=Matrix(nx,nb,(i,j)->evalf(mul(bf[l,iidB[j,l]](a[l]*X[i,l]+b[l]),l=1..d)),datatype=float); C:=LinearAlgebra[LeastSquares](B,Y); cmin2:=cmin*LinearAlgebra[Norm](C,infinity); for j from 1 to nb do if (abs(C[j])M0[i],datatype=float)); d,n:=Dim2(M); maxlen:=1+floor(log10(n)); maxr2:=`if`(type(maxr20,Vector),maxr20,Vector(n,k->maxr20,datatype=float)); effr2:=Mean(maxr2); cor:=evalf(d/2*(GAMMA(d/2)-GAMMA(d/2,effr2/2))/(GAMMA(d/2+1)-GAMMA(d/2+1,effr2/2))); cor0:=`if`(effr20>0,evalf(d/2*(GAMMA(d/2)-GAMMA(d/2,effr20/2))/(GAMMA(d/2+1)-GAMMA(d/2+1,effr20/2))),1); K2:=K0; for iter from 1 to niter do K:=K2; nk:=nops(K); V:=Vector(d,o->Mean(M[o,K]),datatype=float); dM:=Matrix(d,n,(o,k)->M[o,k]-V[o],datatype=float); dM2:=Matrix(d,(o1,o2)->dM[o1,K].dM[o2,K]/(nk-1)*`if`(iter=1,cor0,cor),shape=symmetric,datatype=float); dM2i:=dM2^(-1); r2:=Vector(n,k->Transpose(dM[..,k]).dM2i.dM[..,k],datatype=float); # add(v,v=r2[K])/(nk-1)=d if printout then printf("nk=%*d, mean=[%+{c,}.3f], sigma=[%{c,}.3f]\n",maxlen,nk,V,map(sqrt,Eigenvalues(dM2))) end; K2:=select(k->r2[k][op(v),v[3]-v[2]],[seq([v,evalf(nops(K2)*GAMMA(d/2,v/2)/(GAMMA(d/2)-GAMMA(d/2,effr2/2))),add(`if`(r2[k]evalb(v[1]>x1 and v[1]0..0 and printout) then printf("After xrng(%a): np=%d, x-range=%g..%g\n",xrng,nops(xy1),min(seq(v[1],v=xy1)),max(seq(v[1],v=xy1))) end; if (cythr<>0.5) then y1,y2:=min(seq(v[2],v=xy1)),max(seq(v[2],v=xy1)); ythr:=y1+cythr*(y2-y1); if (cythr>0 and cythr<0.5) then xy2:=select(v->evalb(v[2]0.5) then xy2:=select(v->evalb(v[2]>ythr),xy1) else error("Expected 00) then xy2:=`if`(cythr<0.5,Sort,Rort)(xy2,[2]); if (nops(xy2)<3) then error("3 points are required for pre-interpolation, %2 points are selected by ythr=%3",3,nops(xy2),ythr) end; if (nops(xy2)evalb(v[1]>x1 and v[1]1) then WARNING("Multiple extrema %1",sol); x0:=sol[op(MinIdx([seq(abs(v-x0),v=sol)]))] else x0:=op(sol) end; if not(type(x0,numeric)) then error("No extremum") end; ypp:=add(n*(n-1)*c[n+1]*x0^(n-2),n=2..o); [x0, sqrt(evalf[2*Digits](add(add(n*m*cov[n+1,m+1]*x0^(n+m-2),n=1..o),m=1..o)))/abs(ypp), add(c[n+1]*x0^n,n=0..o), sqrt(evalf[2*Digits](add(add(cov[n+1,m+1]*x0^(n+m),n=0..o),m=0..o))), ypp, unapply(add(c[n+1]*x^n,n=0..o),x), fitrng] end: #hfl: ReorderXYdata ReorderXYdata:=proc( xydata, n::posint:=Dim2(xydata), m::posint:=Dim2(xydata[1][2]), { order::posint:=2, xindexing::procedure:=((A,i)->A[i][1]), yindexing::procedure:=((A,i,j)->A[i][2][j]) },$) local X,Y,y,i,j,id; X:=Vector(n,i->xindexing(xydata,i)); Y:=Matrix(n,m,(i,j)->yindexing(xydata,i,j)); y:=Vector(m); for i from 1 to n do Y[i,..]:=sort(Y[i,..]) end; for i from order+2 to n do id:=SortIdx(Vector(m,j->evalf(CurveFitting[PolynomialInterpolation](X[i-order-1..i-1],Y[i-order-1..i-1,j],evalf(X[i])))),'nolist'); for j from 1 to m do y[id[j]]:=Y[i,j] end; Y[i,..]:=y end; [seq([X[i],Y[i,..]],i=1..n)] end: #hfl: Intersection Intersection:=proc(xydata, n::posint:=Dim2(xydata), { order::nonnegint:=1, xindexing::procedure:=((A,i)->A[i][1]), yindexing::procedure:=((A,i,j)->A[i][2][j]) },$) local X,Y,i1,i2,s,i,x,e1,e2; X:=Vector(n,i->xindexing(xydata,i)); Y:=Matrix(n,2,(i,j)->yindexing(xydata,i,j)); if (order=0) then i1,i2:=1,n else s:=signum(Y[1][1]-Y[1][2]); for i from 2 to n while s*(Y[i][1]-Y[i][2])>0 do end: if (i>n) then error "No intersection is found" end; i1,i2:=i-iquo(order+2,2),i+iquo(order,2); if (type(order,even) and i1=0 and i2<=n ) then i1:=1 elif (type(order,even) and i1>=1 and i2=n+1) then i2:=n elif (i1<1 or i2>n) then error "For order=%1 provide at least %2 interpolation points",order,2*iquo(order+2,2) elif (type(order,even)) then if (abs(Y[i1][1]-Y[i1][2])0) then dx:=(x2-x1)*rdx; x1,x2:=x1-dx,x2+dx end; if (rdy<>0) then dy:=(y2-y1)*rdy; y1,y2:=y1-dy,y2+dy end; x1,x2,y1,y2 end: #hfl: Spline2 Spline2:=proc(XY,cod::{0,1,2,12}:=0,{degree::posint:=3}) local G,e,x; if (cod>0) then G:=`if`(cod=12,Matrix(2,4*degree+1,{(1,1)=1,(2,2*degree+2)=1}),Matrix(1,4*degree+1,{(1,`if`(cod=1,1,2*degree+2))=1})) end; e:=CurveFitting[Spline](XY,x,':-degree'=degree,`if`(cod>0,'endpoints'=G,NULL),_rest); unapply(simplify(piecewise(xXY[-1][1],'undefined',e)),x) end: #hfl: ErrorBar ErrorBar:=proc(x0,r::range,w::realcons:=0,opt::list:=[],$) local x,y1,y2,dx,e; x,y1,y2:=evalf(x0),op(evalf(r)); e:=plottools[line]([x,y1],[x,y2],op(opt)): if (evalf(w)>0) then dx:=evalf(w/2); CURVES([[x-dx,y1],[x+dx,y1]],[[x-dx,y2],[x+dx,y2]],op(e)) else e end end: #hfl: TextOnLine TextOnLine:=proc(x1::realcons,y1::realcons,x2::realcons,y2::realcons,s::string, h::numeric,h2v0::numeric:=0.55,g0::numeric:=0.03,rxy::list(numeric):=[0.5], {textover::boolean:=false}) local h2v,g,w,rx,ry,x0,y0,Dx,Dy,sx,sy,dx,dy; h2v,g:=max(h2v0,g0),min(h2v0,g0); w:=evalf(length(s)*h*h2v); if (nops(rxy)=1) then rx,ry:=rxy[1],rxy[1] elif (nops(rxy)=2) then rx,ry:=rxy[1],rxy[2] else error("Unrecognized rxy: %1",rxy) end; x0,y0:=x1*(1-rx)+x2*rx,y1*(1-ry)+y2*ry; Dx,Dy:=evalf(x2-x1),evalf(y2-y1); if textover then [plot([[x1+g*Dx,y1+g*Dy],[x2-g*Dx,y2-g*Dy]],_rest), plottools[rectangle]([x0-w/2,y0-h/2],[x0+w/2,y0+h/2],'color'="White",'style'="polygon"), plots[textplot]([x0,y0,s],_rest)] else sx,sy:=`if`(Dx<0,-1,1),`if`(Dy<0,-1,1); if (w*Dy^211) then ls:=[seq(op(ls),i=1..ceil(n/11))][..n] end end; [op(before),op(ls),op(after)] end end: #hfl: WritePlotData WritePlotData:=proc( filename::string, ls::list(list(list)), Xformat::string, Yformat::string, com::string:="", { Xlabel::string:="Xvalue", legend::list(string):=[], emptycell::string:="", sep::character:="\t", overwrite::boolean:=false },$) local X,Y,m,n,i,j,v,fmt,u; n:=nops(ls); X:=Sort(convert({seq(seq(sprintf(Xformat,u[1]),u=v),v=ls)},list),v->parse(v)); m:=nops(X); Y:=Matrix(m,n,(i,j)->emptycell); for j from 1 to n do for v in ls[j] do i:=SearchPos(X,sprintf(Xformat,v[1])); if (i>0) then Y[i,j]:=sprintf(Yformat,v[2]) else error("Cannot find index of x-value for data point %1 in %2",v,ls) end end end; if (nops(legend)>0 and nops(legend)<>n) then error("Inconsistent legend %1 for data %2",legend,ls) end; fmt:=cat("%s",cat(sep,"%s")$n); WriteLines(filename,[ `if`(com="",NULL,cat("# ",com)), `if`(legend=[],NULL,sprintf("%s%s%{c*}s",Xlabel,sep,sep,Vector(legend))), seq(sprintf(fmt,X[i],seq(Y[i,j],j=1..n)),i=1..m)], ':-overwrite'=overwrite) end: #hfl: plotps plotps:=proc( filename0::string, pl, opt::string:=psplotoptions, cbb::[integer,integer,integer,integer]:=[0,0,0,0], { waittime::numeric:=0.1, overwrite::boolean:=false, eps2pdf::string:="eps2pdf.bat" }) local p,filename,e,cd,ans; p,filename,e:=ExpandPath(filename0,"p,pn,x"); filename:=`if`(e="" or e=".pdf",cat(filename,".eps"),filename0); if FileTools[Exists](filename) then if overwrite then fremove(filename) else error("File exists: %1",filename) end end; plotsetup('ps','plotoutput'=filename,'plotoptions'=opt); print(plots[display](pl,_rest)); plotsetup('default'); if (cbb<>[0,0,0,0]) then Threads[Sleep](waittime); CorrectEPS(filename,cbb) end; if (e=".pdf") then if FileTools[Exists](filename0) then if overwrite then fremove(filename0) else error("File exists: %1",filename0) end end; if (p<>"") then cd:=currentdir(); currentdir(p) end; ans:=ssystem(cat(eps2pdf," ",ExpandPath(filename,"nx"))); if (p<>"") then currentdir(cd) end; if (ans[1]=0) then fremove(filename) else error("Converter's error: %1",ans) end end; NULL end: #hfl: CorrectEPS CorrectEPS:=proc(filename::string,cbb::list:=[0,0,0,0],$) local ls,fd,s,bb; ls:=ReadLines(filename); fd:=fopen(filename,WRITE,TEXT): for s in ls do if (s[1..8]="%%Pages:") then elif (s[1..14]="%%BoundingBox:") then bb:=sscanf(s[15..-1],"%d %d %d %d"); fprintf(fd,"%%%%BoundingBox: %d %d %d %d\n",op(bb+cbb)) else writeline(fd,s) end end: fclose(fd); NULL end: #hfl: ImageCoreBox ImageCoreBox:=proc(img::Array,threshold::numeric:=1,$) local nx,ny,v,x1,y1,x2,y2; ny,nx,v:=Dim2(img); for x1 from 1 to nx while (min(img[..,x1,..])>=threshold) do end; for y1 from 1 to ny while (min(img[y1,..,..])>=threshold) do end; for x2 from nx by -1 to x1 while (min(img[..,x2,..])>=threshold) do end; for y2 from ny by -1 to y1 while (min(img[y2,..,..])>=threshold) do end; [x1,y1,x2,y2] end: ################################################################################ #cat: HTML and XML tools #hfl: SqueezeHTML SqueezeHTML:=proc(s::string,ls::list:=[],$) local s1,i,j,ls1,v; s1:=StringTools[Trim](Squeeze(s)); i:=1; while (i<=length(s1)) do j:=searchtext("<",s1,i..-1); if (j=0) then break else i:=i+j end; if (s1[i]=" ") then s1:=cat(s1[1..i-1],s1[i+1..-1]) end end; i:=1; while (i<=length(s1)) do j:=searchtext(">",s1,i..-1); if (j=0) then break else i:=i+j end; if (i>2 and s1[i-2]=" ") then s1:=cat(s1[1..i-3],s1[i-1..-1]) end end; ls1:=[seq([length(v),cat(v,">")],v=ls)]; i:=1; while (i<=length(s1)) do j:=searchtext("<",s1,i..-1); if (j=0) then break else i:=i+j end; if (i>2 and s1[i-2]=" " and foldl(`or`,false,seq(s1[i..i+v[1]]=v[2],v=ls1))) then s1:=cat(s1[1..i-3],s1[i-1..-1]) end end; ls1:=[seq([length(v),cat("<",v)],v=ls)]; i:=1; while (i<=length(s1)) do j:=searchtext(">",s1,i..-1); if (j=0) then break else i:=i+j end; if (s1[i]=" " and foldl(`or`,false,seq(s1[max(1,i-2-v[1])..i-2]=v[2],v=ls1))) then s1:=cat(s1[1..i-1],s1[i+1..-1]) end end; s1 end: #hfl: SqueezeHTML SqueezeTable:=proc(s0::string,$) local s,i,s1,pos,s2; s:=SqueezeHTML(s0,["tr","/tr","td","/td","th","/th"]); i:=SearchKeySeq(s,[""])-4; s1:=Trim(s[1..i-1],[" ","\n"]); do pos:=SearchKeySeq(s,["",""],'start'=i,'fulloutput','noerror'); if (pos=0 or pos=[]) then break end; s2:=Trim(s[i..pos[1]-5],[" ","\n"]); s1:=cat( s1, `if`(s2="","",cat("\n",s2)), "\n", StringTools[SubstituteAll](Trim(s[pos[1]..pos[2]-6],[" ","\n"]),"\n"," "), ""); i:=pos[2] end; s1:=cat(s1,"\n",Trim(s[i..-1],[" ","\n"])); SqueezeHTML(s1,["tr","/tr","td","/td","th","/th"]) end: #hfl: ReadHref ReadHref:=proc(s0::string) local i,s; s:=SqueezeHTML(StringTools[SubstituteAll](s0,"\n",""),["a","/a"]); i:=searchtext("\">",s,10..-5); if (s[1..9]<>""" or i=0) then error "Incorrect href structure: %1",s else [s[i+11..-5],s[10..i+8]] end end: #hfl: ReadHref ResolveURI:=proc(URI::string,thisfile::string,$)::string; local s,n,ls; if (URI[1..8]="file:///") then URI[9..-1] elif (URI[1..7]="http://") then URI else s:=`if`(URI=".","",URI); n:=0; while (s[1..3]="../") do s:=s[4..-1]; n:=n+1 end; ls:=StringTools[Split](thisfile,"/"); if (nops(ls)true),{all::boolean:=false},$) local tb,i,v; if not(XMLTools[IsElement](xml) or XMLTools[IsTree](xml)) then error("First argument must be XML element or tree but received %1",xml) end; tb:=table(): for i from 1 to XMLTools[ContentModelCount](xml) do v:=XMLTools[GetChild](xml,i); if (XMLTools[IsElement](v) and XMLTools[ElementTypeName](v)=nam and f(v)) then if (XMLTools[ContentModelCount](v)=1 and XMLTools[IsText](XMLTools[GetChild](v,1))) then v:=XMLTools[TextText](XMLTools[GetChild](v,1)) end; if not(all) then return v end; tb[i]:=v end end; `if`(all,convert(tb,list),NULL) end: #hfl: ReadXML ReadXML:=proc(filename::string,{reload::boolean:=false},$) local ff,ls,i,XMLdeclaration,DOCTYPEdeclaration,xml; ff:=proc(x) local t,u,v; if type(x,function) then t:=op(0,x); if (t=_XML_Element) then u:=[XMLTools[ElementName](x),op(XMLTools[AttributeTable](x))]; v:=map(ff,op(3,x)); if (v=[]) then [u] elif (nops(v)=1 and type(v[1],string)) then [[op(u),v[1]]] else [u,op(v)] end elif (t=_XML_Text) then op(1,x) else XMLTools[PrintToString](x) end else x end end; ls:=ReadLines(filename,':-reload'=reload); i:=1; if (TrimLeft(ls[i][..5])="", dtd::string:="", overwrite::boolean:=false},$) local ff,XMLdeclaration,DOCTYPEdeclaration; ff:=proc(x) local v,att,u; v:=x[1]; att:=cat(seq(sprintf(" %s='%s'",op(u)),u=v[2])); if (nops(x)=1) then if (nops(v)=3) then sprintf("<%s%s>%s\n",v[1],att,v[-1],v[1]) else sprintf("<%s%s />\n",v[1],att) end else cat(sprintf("<%s%s>\n",v[1],att),seq(ff(u),u=x[2..]),sprintf("\n",v[1])) end end; WriteLines(filename,[`if`(xmld="",NULL,xmld),`if`(dtd="",NULL,dtd),ff(xml)],':-overwrite'=overwrite) end: ################################################################################ #cat: Physics tools #hfl: SU #SU q2vu:=proc(expr) local q,v,u,e; q:=simplify(expr); e:=selectfun(q,Units[Unit]); if nops(e)=0 then u:=1; v:=q elif nops(e)=1 then u:=op(e); v:=simplify(q/u) else error "Unrecognized format in %1",q end; v,u end: #hfl: GetConstantSU GetConstantSU:=proc(const::name) local opts,e,v,u; opts:=[_passed][2..-1]; e:=ScientificConstants[Constant](const,'system'='SU'); e:=ScientificConstants[GetValue](e)*ScientificConstants[GetUnit](e); if ('esu' in opts) then e:=simplify(e,'esu') elif ('emu' in opts) then e:=simplify(e,'emu') else e:=simplify(e) end; v,u:=q2vu(e); if ('value' in opts) then v elif ('unit' in opts) then u else e end end: #hfl: ConvertSU ConvertSU:=proc(expr,sys::name,$) local e,v,u; if (sys='SI') then v,u:=q2vu(expr); e:=subs('esu'=10/cc*'A'*'s','emu'=10*'A'*'s',parse(convert(op(u),string))) elif (sys='CGS') then v,u:=q2vu(expr); e:=parse(convert(op(u),string)); if has(e,'A') then error "A unit is present in %1",expr end; e:=subs( 'esu'=sqrt('g'*'cm'^3/'s'^2), 'emu'=sqrt('g'*'cm'^3/'s'^2), 'K'=kB*'g'*'cm'^2/'s'^2, 'mol'=NA, 'cd'=10^7/683*'g'*'cm'^2/'s'^3, e) elif (sys='cm') then v,u:=q2vu(ConvertSU(expr,'CGS')); e:=subs('s'=cc*'cm','g'=cc/hb/'cm',parse(convert(op(u),string))) else error "Unrecognized system of units: %1",sys end; simplify(Units[Unit](v*e)) end: ################################################################################ #cat: System tools #hfl: BuildHelp BuildHelp:=proc( source::string, hdbfile::string:="", { rootkey::string:="#root:", catkey::string:="#cat:", tockey::string:="#toc:", hflkey::string:="#hfl:", ext::string:=".mw", test::boolean:=false, printout::boolean:=test},$) local v,lib,folder,s,root,category,n1,n2,n3,n4,S,N,l,helpfile,tocentry,topic,i; if (MapleVersion>=18) then error("Maple 18 bug: unable to parse mw-files as help pages. Open required mw-files directly to get help") end; if (hdbfile="") then lib:=NULL else if FileTools[Exists](hdbfile) then WARNING("hdb-file exists and can not be deleted during Maple session. Build over. For sure, remove hdb-file.") end; lib:=hdbfile end; folder,s:=ExpandPath(source,"p,n"); if printout then printf("folder=%s\n",folder) end; root:=cat(s,"/"); category:=""; n1,n2,n3,n4:=length(rootkey),length(catkey),length(tockey),length(hflkey); S:=ReadLines(source); N:=nops(S); for l from 1 to N do s:=S[l]; if (s[..n1]=rootkey) then root:=cat(Trim(s[n1+1..]),"/"); if printout then printf("%s\n",s) end elif (s[..n2]=catkey) then category:=cat(Trim(s[n2+1..]),"/"); if printout then printf("%s\n",s) end elif (s[..n4]=hflkey) then helpfile:=cat(folder,Trim(s[n4+1..]),ext); if printout then printf("%s\n",s) end; if not(FileTools[Exists](helpfile)) then error("Helpfile \"%1\" does not exist",Trim(s[n4+1..])) end; if (l0) then error "Module in module or no module end" end; s1:=s[i..j]: for k from i-2 by -1 to 1 while (s[k]<>"\n") do end; tb[i]:=cat(prefix,Trim(s[k+1..i-2],[" ",":","="])); assign(convert(tb[i],name),s1) end; convert(tb,list) end: #hfl: ReadModules AssignModule:=proc(moduledef::string) local s,m; if (SearchText(" ",moduledef)>0) then s:=moduledef else s:=eval(convert(cat("s_",moduledef),name)); if not(type(s,string)) then ReadModules(moduledef); s:=eval(s) end end; m:=parse(s,'statement'); if member('ModuleInit',m) then m:-ModuleInit(_rest) end; m end: #hfl: Check4Updates Check4Updates:=proc( root::string:="http://cdn.rawgit.com/zhugayevych/zhugayevych.github.io/master/maple/", packages::list(string):=["BasicTools","SSH","FiniteGroups","RandomWalk","MolMod","AstroTools"], {openlinks::boolean:=false},$) local n,p,v1,v2,s,s2; n:=max(map(length,packages)); for p in packages do v1,v2:=0,0; s:=sprintf("%-*s ",n,p); try v1:=eval(parse(cat(p,":-Setup('version')"))); s:=sprintf("%s%d",s,v1) catch: s:=cat(s,"not installed") end; try s2:=HTTP:-Get(cat(root,p,"/",p,".ini"))[2]; v2:=ReadBody(s2,"version:=",";",'format'="%d"); if (v1=v2) then s:=sprintf("%s=%d",s,v2) elif (v1%d - you have newer version",s,v2) end catch: s:=cat(s," - no online information") end; printf("%s\n",s) end; end: #hfl: TabulatedFunction TabulatedFunction:=proc() AssignModule(s_TabulatedFunction,_rest) end: s_TabulatedFunction:="module()\n export ModuleInit,ModuleApply,XV,Remove,Read,Save;\n local f,rho,defdr;\n \n ModuleInit:=proc(ForXVorFilename,dr::positive,rhof::procedure:=((x,y)->`if`(type(x,atomic),abs(x-y),evalf(sqrt(add(v^2,v=x-y))))),{Fdim::posint:=1},$)\n if type(ForXVorFilename,list) then\n XV:=ForXVorFilename;\n f:=x->'undefined'\n elif type(ForXVorFilename,string) then\n XV:=[];\n Read(filename,Fdim);\n f:=x->'undefined'\n else\n XV:=[];\n f:=ForXVorFilename end;\n defdr:=dr;\n rho:=rhof;\n NULL\n end:\n \n ModuleApply:=proc(x,dr::nonnegative:=defdr,$)\n local ls,i,xv;\n if (XV=[]) then\n xv:=[evalf(x),evalf(f(x))]; XV:=[xv]\n else\n ls:=[seq(rho(x,xv[1]),xv=XV)];\n i:=op(MinIdx(ls));\n if (ls[i]>dr) then xv:=[evalf(x),evalf(f(x))]; XV:=[op(XV),xv] else xv:=XV[i] end end;\n xv\n end:\n \n Remove:=proc(x,$)\n local i;\n if (XV<>[]) then\n i:=op(MinIdx([seq(rho(x,xv[1]),xv=XV)])):\n XV:=subsop(i=NULL,XV) end;\n NULL\n end:\n \n Read:=proc(filename::string,Fdim::posint:=1,convert::procedure:=(V->Vector[column](V)),$)\n local M,n,m,Xdim;\n M:=op(fscanf(filename,""%{;h}fm""));\n n,m:=Dim2(M);\n Xdim:=m-Fdim;\n XV:=[op(XV),seq([`if`(Xdim=1,M[i,1],convert(M[i,1..Xdim])),`if`(Fdim=1,M[i,-1],convert(M[i,-Fdim..-1]))],i=1..n)];\n ls\n end:\n \n Save:=proc(filename::string,{digits::posint:=min(10,Digits-2),overwrite::boolean:=false},$)\n local fd,xv;\n if (XV=[]) then\n WARNING(""Nothing to save"")\n else\n if (filename<>"""") then\n if (FileTools[Exists](filename) and not(overwrite)) then error ""File exists, %1"",filename end;\n fd:=fopen(filename,WRITE,TEXT)\n else\n fd:=fopen('terminal',WRITE) end;\n for xv in XV do fprintf(fd,""%.*f %.*f\n"",\n digits, `if`(type(xv[1],float),xv[1],convert(xv[1],Vector)),\n digits, `if`(type(xv[2],float),xv[2],convert(xv[2],Vector))) end;\n if (filename<>"""") then fclose(fd) end end;\n nops(XV)\n end:\n \n end": #hfl: ReduceFloat ReduceFloat:=(f::atomic)->parse(TrimRight(sprintf("%a",evalf(f)),["0"]),'statement'): #hfl: ReduceFloat ReduceFloat2:=proc(f::atomic,$) local s; s:=TrimRight(sprintf("%a",evalf(f)),["0"]); if (s[-1]=".") then s:=s[..-2] end; parse(s,'statement') end: #hfl: ReduceFloat ReduceFloat3:=proc(f::atomic,base::posint,digits::posint:=Digits-2,$) local v; v:=round(f*base); `if`(abs(f-v/base)<10^(-digits),v/base,f) end: #hfl: ReduceFloat ReduceFloatGrid:=proc(G0::list(list),id::{posint,list(posint)},digits::posint,digits2::posint:=digits-2,$) local lso,no,kmins,lsn,n,G,i,o,v; lso:=`if`(type(id,posint),[$1..id],id); no:=nops(lso); kmins:=[seq(MinVal(remove(`<`,[seq(abs(v[o]),v=G0),1],10^(-digits))),o=lso)]; lsn:=map(v->1/v,map(convert,kmins,rational,digits2)); for n in lsn do if not(type(n,integer)) then error("Noninteger result: %1",lsn) end end; G:=map(v->subsop(seq(lso[o]=round(v[lso[o]]*lsn[o])/lsn[o],o=1..no),v),G0); for i from 1 to nops(G) do for o in lso do if (abs(G[i][o]-G0[i][o])>10^(-digits)) then error("Deviation is too large for newG[%1][%2]=%3<>%4",i,o,G[i][o],G0[i][o]) end end end; G,lsn end: #hfl: Timing Timing:=proc(nleft::extended_numeric:=nLeft,{nmin::posint:=3,timestep::positive:=5},$) global InitialTime,PrintTime,nDone,nLeft; local t,tleft; nDone:=nDone+1; nLeft:=nleft-1; t:=time[real](); if (t>PrintTime and nDone>=nmin) then tleft:=(t-InitialTime)*nLeft/nDone; if (tleft<100) then printf("%ds left to finish at %ds\n",round(tleft),round(t+tleft)) elif (tleft<3000) then printf("%dm left to finish at %s\n",round(tleft/60),StringTools[FormatTime]("%R",round(tleft))) elif (tleft<80000) then printf("%dh left to finish at %s\n",round(tleft/3600),StringTools[FormatTime]("%R",round(tleft))) elif (tleft