Home
Manual
Packages
Global Index
Keywords
Quick Reference
|
func zpad (x,pad=,both=){
local s,s2,res;
if(is_void(pad)) pad=1;
s =dimsof(x);
s2= grow(s(1),s(2:)+pad*(is_void(both)?1:2));
res= array(structof(x),s2);
edge=is_void(both)?0:-pad;
if (s(1)==1) res(1+pad:edge)=x;
if (s(1)==2) res(1+pad:edge,1+pad:edge)=x;
if (s(1)==3) res(1+pad:edge,1+pad:edge,1+pad:edge)=x;
if (s(1)==4) res(1+pad:edge,1+pad:edge,1+pad:edge,1+pad:edge)=x;
if (s(1)==5) res(1+pad:edge,1+pad:edge,1+pad:edge,1+pad:edge,1+pad:edge)=x;
if (s(1)==6) res(1+pad:edge,1+pad:edge,1+pad:edge,1+pad:edge,1+pad:edge,1+pad:edge)=x;
if (s(1)==7) res(1+pad:edge,1+pad:edge,1+pad:edge,1+pad:edge,1+pad:edge,1+pad:edge,1+pad:edge)=x;
return res;
}
func bsplit (tens)
/* DOCUMENT bsplit splits tensor tens into vectors defining the
unit sub tensors; works in dimensions 1 to 5.
*/
{
d= dimsof(tens); d0=d(1);
if (d0==1)
{
tens2=transpose([tens,roll(tens,-1)])(,:-1);
}
if (d0==2)
{
n1=dimsof(tens)(2);
n2=dimsof(tens)(3);
tens2=array(structof(tens(1)),2,2,(n1-1)*(n2-1));
for(i1=1;i1<=n1-1;i1++)
for(i2=1;i2<=n2-1;i2++)
tens2(*,i2+(n2-1)*(i1-1))=
tens(1+(i1-1):i1+1,1+(i2-1):i2+1)(*);
}
if (d0==3)
{
n1=dimsof(tens)(2);
n2=dimsof(tens)(3);
n3=dimsof(tens)(4);
tens2=array(structof(tens(1)),2,2,2,(n1-1)*(n2-1)*(n3-1));
for(i1=1;i1<=n1-1;i1++)
for(i2=1;i2<=n2-1;i2++)
for(i3=1;i3<=n3-1;i3++)
tens2(*,i3+(n3-1)*(i2-1)+(n3-1)*(n2-1)*(i1-1))=
tens(1+(i1-1):i1+1,1+(i2-1):i2+1,1+(i3-1):i3+1)(*);
}
if (d0==4)
{
n1=dimsof(tens)(2);
n2=dimsof(tens)(3);
n3=dimsof(tens)(4);
n4=dimsof(tens)(5);
tens2=array(structof(tens(1)),2,2,2,2,(n1-1)*(n2-1)*(n3-1)*(n4-1));
for(i1=1;i1<=n1-1;i1++)
for(i2=1;i2<=n2-1;i2++)
for(i3=1;i3<=n3-1;i3++)
for(i4=1;i4<=n4-1;i4++)
tens2(*,i4+(n4-1)*(i3-1)+(n4-1)*(n3-1)*(i2-1)+(n4-1)*(n3-1)*(n2-1)*(i1-1))=
tens(1+(i1-1):i1+1,1+(i2-1):i2+1,1+(i3-1):i3+1,1+(i4-1):i4+1)(*);
}
if (d0==5)
{
n1=dimsof(tens)(2);
n2=dimsof(tens)(3);
n3=dimsof(tens)(4);
n4=dimsof(tens)(5);
n5=dimsof(tens)(6);
tens2=array(structof(tens(1)),2,2,2,2,2,(n1-1)*(n2-1)*(n3-1)*(n4-1)*(n5-1));
for(i1=1;i1<=n1-1;i1++)
for(i2=1;i2<=n2-1;i2++)
for(i3=1;i3<=n3-1;i3++)
for(i4=1;i4<=n4-1;i4++)
for(i5=1;i5<=n5-1;i5++)
tens2(*,i5+(n5-1)*(i4-1)+(n5-1)*(n4-1)*(i3-1)+
(n5-1)*(n4-1)*(n3-1)*(i2-1)+(n5-1)*(n4-1)*(n3-1)*(n2-1)*(i1-1))=
tens(1+(i1-1):i1+1,
1+(i2-1):i2+1,1+(i3-1):i3+1,1+(i4-1):i4+1,1+(i5-1):i5+1)(*);
}
return tens2;
}
func vsplit3D (x,y,z,binx,biny,binz)
/* DOCUMENT vsplit partitions vector n^3 subregions
it returns for the first componant an array of n^3 pointers corresponding to the different regions
the index points to an index of the componants which fall in
region i and for the second componant an array describing the edge vector
of each region
EXEMPLE x=random(250000); y=random(250000); z=random(250000);
w =vsplit(x,y,z,span(0,1,3),span(0,1,3),span(0,1,3));
pl,x(*w(3,1)),y(*w(3,1)),color=-9;
nn=10; pp=dimsof(w)(2);
u=array(0.,nn,nn,nn,pp);
for(i=1;i<=pp;i++) {
tt=[x(*w(i,1)),y(*w(i,1)),z(*w(i,1))];
indx=(*w(i,2));
u(,,,i)=histo3d(tt,span(indx(1,1),indx(2,1),nn+1),
span(indx(1,2),indx(2,2),nn+1),
span(indx(1,3),indx(2,3),nn+1));
}
v=splitb(u);
pli,v(,avg,);
try also
w =vsplit(x,y,z,span(0,0.5,2),span(0.25,0.5,2),span(0,0.5,2));
to look at small region around a given point.
SEE ALSO split, bsplit, splitb, esplit
*/
{
n=dimsof(binx)(0);
ux=binx(,-:1:n)(,,-:1:n);
uy=binx(-:1:n,)(,,-:1:n);
uz=binx(-:1:n,)(-:1:n,,);
bux= bsplit(ux);
buy= bsplit(uy);
buz= bsplit(uz);
p = dimsof(bux)(0);
w=array(pointer,p);
vec=array(pointer,p);
for(i=1;i<=p;i++) {
tt =((x> bux(1,1,1,i))*(x<= bux(2,1,1,i))*
(y> buy(1,1,1,i))*(y<= buy(1,2,1,i))*
(z> buz(1,1,1,i))*(z<= buz(1,1,2,i)))
if (is_array(where(tt))) w(i)=&(where(tt));
vec(i) = &([[bux(1,1,1,i),bux(2,1,1,i)],[buy(1,1,1,i),buy(1,2,1,i)],[ buz(1,1,1,i),
buz(1,1,2,i)]]);
}
return [w,vec];
}
func vsplit2D (x,y,binx,biny)
/* DOCUMENT vsplit partitions vector n^2 subregions
it returns for the first componant an array of n^2 pointers corresponding to the different regions
the index points to an index of the componants which fall in
region i and for the second componant an array describing the edge vector
of each region
EXEMPLE x=random(25000); y=random(25000);
w =vsplit2D(x,y,span(0,1,3),span(0,1,3));
pl,x(*w(3,1)),y(*w(3,1)),color=-9;
nn=10; pp=dimsof(w)(2);
u=array(0.,nn,nn,pp);
for(i=1;i<=pp;i++) {
tt=[x(*w(i,1)),y(*w(i,1))];
indx=(*w(i,2));
u(,,,i)=histo2d(tt,span(indx(1,1),indx(2,1),nn+1),
span(indx(1,2),indx(2,2),nn+1));
}
v=splitb(u);
pli,v(,avg,);
try also
w =vsplit(x,y,span(0,0.5,2),span(0.25,0.5,2));
to look at small region around a given point.
SEE ALSO split, bsplit, splitb, esplit
*/
{
n=dimsof(binx)(0);
ux=binx(,-:1:n);
uy=binx(-:1:n,);
bux= bsplit(ux);
buy= bsplit(uy);
p = dimsof(bux)(0);
w=array(pointer,p);
vec=array(pointer,p);
for(i=1;i<=p;i++) {
tt =((x> bux(1,1,i))*(x<= bux(2,1,i))*
(y> buy(1,1,i))*(y<= buy(1,2,i)));
if (is_array(where(tt))) w(i)=&(where(tt));
vec(i) = &([[bux(1,1,i),bux(2,1,i)],[buy(1,1,i),buy(1,2,i)]]);
}
return [w,vec];
}
func vsplit1D (x,binx)
/* DOCUMENT vsplit partitions vector n^2 subregions
it returns for the first componant an array of n^2 pointers corresponding to the different regions
the index points to an index of the componants which fall in
region i and for the second componant an array describing the edge vector
of each region
EXEMPLE x=random(25000);
w =vsplit1D(x,span(0,1,5));
pl,x(*w(3,1)),color=-9;
nn=10; pp=dimsof(w)(2);
SEE ALSO split, bsplit, splitb, esplit
*/
{
n=dimsof(binx)(0);
ux=binx;
bux= bsplit(ux);
p = dimsof(bux)(0);
w=array(pointer,p);
vec=array(pointer,p);
for(i=1;i<=p;i++) {
tt =((x> bux(1,i))*(x<= bux(2,i)));
if (is_array(where(tt))) w(i)=&(where(tt));
vec(i) = &([[bux(1,i),bux(2,i)]]);
}
return [w,vec];
}
func split (tens,n)
/* DOCUMENT splits tensor tens into n^p subtensor
where p = dimsof(tens)(1) and returns
a tensor of dimension p+1 where the last index runs alongs
the subtensors; works in dimensions 2 to 5.
*/
{ local d,d0,d1,d2,d3,d4,d5,i1,i2,i3,i4,i5,tens2;
d= dimsof(tens); d0=d(1);
tens2= array(structof(tens(1)),grow(d0+1,grow((d/n)(2:),n^d0)));
if (d0==2)
{
d1=(d/n)(2);
d2=(d/n)(3);
for(i1=1;i1<=n;i1++)
for(i2=1;i2<=n;i2++)
tens2(*,i2+n*(i1-1))=
tens(1+d1*(i1-1):d1*i1,1+d2*(i2-1):d2*i2)(*);
}
if (d0==3)
{
d1=(d/n)(2);
d2=(d/n)(3);
d3=(d/n)(4);
for(i1=1;i1<=n;i1++)
for(i2=1;i2<=n;i2++)
for(i3=1;i3<=n;i3++)
tens2(*,i3+n*(i2-1)+n^2*(i1-1))=
tens(1+d1*(i1-1):d1*i1,1+d2*(i2-1):d2*i2,1+d3*(i3-1):d3*i3)(*);
}
if (d0==4)
{
d1=(d/n)(2);
d2=(d/n)(3);
d3=(d/n)(4);
d4=(d/n)(5);
for(i1=1;i1<=n;i1++)
for(i2=1;i2<=n;i2++)
for(i3=1;i3<=n;i3++)
for(i4=1;i4<=n;i4++)
tens2(*,i4+n*(i3-1)+n^2*(i2-1)+n^3*(i1-1))=
tens(1+d1*(i1-1):d1*i1,1+d2*(i2-1):d2*i2,1+d3*(i3-1):d3*i3,
1+d4*(i4-1):d4*i4)(*);
}
if (d0==5)
{
d1=(d/n)(2);
d2=(d/n)(3);
d3=(d/n)(4);
d4=(d/n)(5);
d5=(d/n)(6);
for(i1=1;i1<=n;i1++)
for(i2=1;i2<=n;i2++)
for(i3=1;i3<=n;i3++)
for(i4=1;i4<=n;i4++)
for(i5=1;i5<=n;i5++)
tens2(*,i5+n*(i4-1)+n^2*(i3-1)+n^3*(i2-1)+n^4*(i1-1))=
tens(1+d1*(i1-1):d1*i1,1+d2*(i2-1):d2*i2,1+d3*(i3-1):d3*i3,
1+d4*(i4-1):d4*i4,1+d5*(i5-1):d5*i5)(*);
}
return tens2;
}
func splitb (tens,pad=)
/* DOCUMENT rebuild tensor tens from n^p subtensor
* SEE ALSO split, bsplit, vsplit, esplit
EXAMPLE
x=eesplit(z,4,pad=1); //adds edges
pli,splitb(x,pad=1)-z; //reconstruct
*/
{ local d,d0,d1,d2,d3,d4,d5,i1,i2,i3,i4,i5,tens2;
d= dimsof(tens);d0=d(1);
if(pad)
{
if (d0==3) tens=tens(pad+1:-pad,pad+1:-pad,);
if (d0==4) tens=tens(pad+1:-pad,pad+1:-pad,pad+1:-pad,);
if (d0==5) tens=tens(pad+1:-pad,pad+1:-pad,pad+1:-pad,pad+1:-pad,);
}
d= dimsof(tens); d0=d(1); n=long(d(0)^(1./(d0-1))+1e-12);
tens2= array(structof(tens(1)),grow(d0-1,(d*n)(2:-1)));
if (d0==3)
{
d1=(d)(2);
d2=(d)(3);
for(i1=1;i1<=n;i1++)
for(i2=1;i2<=n;i2++) {
tt=array(0.,d1,d2);
tt(*)=tens(*,i2+n*(i1-1))(*);
tens2(1+d1*(i1-1):d1*i1,1+d2*(i2-1):d2*i2)=tt ;
}
}
if (d0==4)
{
d1=(d)(2);
d2=(d)(3);
d3=(d)(4);
for(i1=1;i1<=n;i1++)
for(i2=1;i2<=n;i2++)
for(i3=1;i3<=n;i3++) {
tt=array(0.,d1,d2,d3);
tt(*)=tens(*,i3+n*(i2-1)+n^2*(i1-1))(*);
tens2(1+d1*(i1-1):d1*i1,1+d2*(i2-1):d2*i2,1+d3*(i3-1):d3*i3)=tt ;
}
}
if (d0==5)
{
d1=(d)(2);
d2=(d)(3);
d3=(d)(4);
d4=(d)(5);
for(i1=1;i1<=n;i1++)
for(i2=1;i2<=n;i2++)
for(i3=1;i3<=n;i3++)
for(i4=1;i4<=n;i4++){
tt=array(0.,d1,d2,d3,d4);
tt(*)=tens(*,i4+n*(i3-1)+n^2*(i2-1)+n^3*(i1-1))(*);
tens2(1+d1*(i1-1):d1*i1,1+d2*(i2-1):d2*i2,1+d3*(i3-1):d3*i3,1+d4*(i4-1):d4*i4)=tt ;
}
}
return tens2;
}
func esplit (tens,n)
/* DOCUMENT splits tensor tens into n^p subtensor
where p = dimsof(tens)(1) and returns
a tensor of dimension p+1 where the last index runs alongs
the subtensors; works in dimensions 2 to 5.
esplit includes an overlap between subtensors while padding
the edge with zeros.
SEE ALSO split
EXAMPLE
#include "ninteg.i"
x= span(0,15,128)(,-:1:128);
y =span(0,3,128)(-:1:128,);
z =sin(x-y);
x1=esplit(x,8); y1=esplit(y,8); z1=esplit(z,8);
u=array(0.,dimsof(x1)(0));
for(i=1;i<=dimsof(x1)(0);i++){u(i)=ninteg(z1(..,i),y1(..,i),x1(..,i));}
u(sum)
ninteg(z,y,x)
x= span(0,5,28)(,-:1:28)(,,-:1:28);
y =span(0,3,28)(-:1:28,)(,,-:1:28);
z =span(0,2,28)(-:1:28,)(-:1:28,,);
w =sin(x+y+z);
x1=esplit(x,4); y1=esplit(y,4); z1=esplit(z,4); w1=esplit(w,4);
u=array(0.,dimsof(x1)(0));
for(i=1;i<=dimsof(x1)(0);i++){u(i)=ninteg(w1(..,i),z1(..,i),y1(..,i),x1(..,i));}
u(sum)
ninteg(w,z,y,x)
x= span(0,5,28)(,-:1:28)(,,-:1:28)(,,,-:1:28);
y =span(0,3,28)(-:1:28,)(,,-:1:28)(,,,-:1:28);
z =span(0,2,28)(-:1:28,)(-:1:28,,)(,,,-:1:28);
w =span(0,2,28)(-:1:28,)(-:1:28,,)(-:1:28,,,);
v =sin(x+y+z-w);
x1=esplit(x,4); y1=esplit(y,4); z1=esplit(z,4); w1=esplit(w,4);
v1= esplit(v,4);
u=array(0.,dimsof(x1)(0));
for(i=1;i<=dimsof(x1)(0);i++){
u(i)=ninteg(v1(..,i),w1(..,i),z1(..,i),y1(..,i),x1(..,i));}
u(sum)
ninteg(v,w,z,y,x)
*/
{ local d,d0,d1,d2,d3,d4,d5,i1,i2,i3,i4,i5,tens2;
d= dimsof(tens); d0=d(1);
tens2= array(structof(tens(1)),grow(d0+1,(d/n)(2:)+1,n^d0));
if (d0==2)
{
d1=(d/n)(2);
d2=(d/n)(3); tens=zpad(tens);
for(i1=1;i1<=n;i1++)
for(i2=1;i2<=n;i2++)
tens2(*,i2+n*(i1-1))=
tens(d1*(i1-1)+1:d1*i1+1,d2*(i2-1)+1:d2*i2+1)(*);
}
if (d0==3)
{
d1=(d/n)(2);
d2=(d/n)(3);
d3=(d/n)(4); tens=zpad(tens);
for(i1=1;i1<=n;i1++)
for(i2=1;i2<=n;i2++)
for(i3=1;i3<=n;i3++)
tens2(*,i3+n*(i2-1)+n^2*(i1-1))=
tens(1+d1*(i1-1):d1*i1+1,1+d2*(i2-1):d2*i2+1,1+d3*(i3-1):d3*i3+1)(*);
}
if (d0==4)
{
d1=(d/n)(2);
d2=(d/n)(3);
d3=(d/n)(4);
d4=(d/n)(5); tens=zpad(tens);
for(i1=1;i1<=n;i1++)
for(i2=1;i2<=n;i2++)
for(i3=1;i3<=n;i3++)
for(i4=1;i4<=n;i4++)
tens2(*,i4+n*(i3-1)+n^2*(i2-1)+n^3*(i1-1))=
tens(1+d1*(i1-1):d1*i1+1,1+d2*(i2-1):d2*i2+1,1+d3*(i3-1):d3*i3+1,
1+d4*(i4-1):d4*i4+1)(*);
}
if (d0==5)
{
d1=(d/n)(2);
d2=(d/n)(3);
d3=(d/n)(4);
d4=(d/n)(5);
d5=(d/n)(6); tens=zpad(tens);
for(i1=1;i1<=n;i1++)
for(i2=1;i2<=n;i2++)
for(i3=1;i3<=n;i3++)
for(i4=1;i4<=n;i4++)
for(i5=1;i5<=n;i5++)
tens2(*,i5+n*(i4-1)+n^2*(i3-1)+n^3*(i2-1)+n^2*(i1-1))=
tens(1+d1*(i1-1):d1*i1+1,1+d2*(i2-1):d2*i2+1,1+d3*(i3-1):d3*i3+1,
1+d4*(i4-1):d4*i4+1,1+d5*(i5-1):d5*i5+1)(*);
}
return tens2;
}
func eesplit (tens,n,pad=)
/* DOCUMENT
like
esplit but with wider edges
SEE ALSO esplit
EXAMPLE
x= span(0,15,128)(,-:1:128);y =span(0,3,128)(-:1:128,);z =sin(x-y);
pli,splitb(esplit2(z,8,pad=8));
*/
{ local d,d0,d1,d2,d3,d4,d5,i1,i2,i3,i4,i5,tens2;
if(is_void(pad)) pad=1;
d= dimsof(tens); d0=d(1);
tens2= array(structof(tens(1)),grow(d0+1,(d/n)(2:)+2*pad,n^d0));
if (d0==2)
{
d1=(d/n)(2);
d2=(d/n)(3); tens=zpad(tens,pad=pad,both=1);
for(i1=1;i1<=n;i1++)
for(i2=1;i2<=n;i2++)
tens2(*,i2+n*(i1-1))=
tens(d1*(i1-1)+1:d1*i1+2*pad,d2*(i2-1)+1:d2*i2+2*pad)(*);
}
if (d0==3)
{
d1=(d/n)(2);
d2=(d/n)(3);
d3=(d/n)(4); tens=zpad(tens,pad=pad,both=1);
for(i1=1;i1<=n;i1++)
for(i2=1;i2<=n;i2++)
for(i3=1;i3<=n;i3++)
tens2(*,i3+n*(i2-1)+n^2*(i1-1))=
tens(1+d1*(i1-1):d1*i1+2*pad,1+d2*(i2-1):d2*i2+2*pad,1+d3*(i3-1):d3*i3+2*pad)(*);
}
if (d0==4)
{
d1=(d/n)(2);
d2=(d/n)(3);
d3=(d/n)(4);
d4=(d/n)(5); tens=zpad(tens,pad=pad,both=1);
for(i1=1;i1<=n;i1++)
for(i2=1;i2<=n;i2++)
for(i3=1;i3<=n;i3++)
for(i4=1;i4<=n;i4++)
tens2(*,i4+n*(i3-1)+n^2*(i2-1)+n^3*(i1-1))=
tens(1+d1*(i1-1):d1*i1+2*pad,1+d2*(i2-1):d2*i2+2*pad,1+d3*(i3-1):d3*i3+2*pad,
1+d4*(i4-1):d4*i4+2*pad)(*);
}
if (d0==5)
{
d1=(d/n)(2);
d2=(d/n)(3);
d3=(d/n)(4);
d4=(d/n)(5);
d5=(d/n)(6); tens=zpad(tens,pad=pad,both=1);
for(i1=1;i1<=n;i1++)
for(i2=1;i2<=n;i2++)
for(i3=1;i3<=n;i3++)
for(i4=1;i4<=n;i4++)
for(i5=1;i5<=n;i5++)
tens2(*,i5+n*(i4-1)+n^2*(i3-1)+n^3*(i2-1)+n^2*(i1-1))=
tens(1+d1*(i1-1):d1*i1+2*pad,1+d2*(i2-1):d2*i2+2*pad,1+d3*(i3-1):d3*i3+2*pad,
1+d4*(i4-1):d4*i4+2*pad,1+d5*(i5-1):d5*i5+2*pad)(*);
}
return tens2;
}
func recenter (cube,coor)
/* DOCUMENT recenter cube on coordinates coor
EXAMPLE
x=span(0,1,15)(,-:1:16)(,,-:1:17);
y=span(0,1,16)(-:1:15,)(,,-:1:17);
z=span(0,1,17)(-:1:16,)(-:1:15,,);
u=abs(x-0.5,y-0.5,z-0.5);
v=recenter(u,[3,4,5]);
need not be a cube; can be a rectangle or a hypercube
SEE ALSO:
*/
{ dims =dimsof(cube);
if (dims(1)==3){
n1=dims(2);
n2=dims(3);
n3=dims(4);
w1=indgen(n1);
w2=indgen(n2);
w3=indgen(n3);
return cube(((w1+coor(1)-1) % n1)+1,
((w2+coor(2)-1) % n2)+1,
((w3+coor(3)-1) % n3)+1);
}
if (dims(1)==2){
n1=dims(2);
n2=dims(3);
w1=indgen(n1);
w2=indgen(n2);
return cube(((w1+coor(1)-1) % n1)+1,
((w2+coor(2)-1) % n2)+1);
}
if (dims(1)==4){
n1=dims(2);
n2=dims(3);
n3=dims(4);
n4=dims(5);
w1=indgen(n1);
w2=indgen(n2);
w3=indgen(n3);
w4=indgen(n4);
return cube(((w1+coor(1)-1) % n1)+1,
((w2+coor(2)-1) % n2)+1,
((w3+coor(3)-1) % n3)+1,
((w4+coor(4)-1) % n4)+1);
}
}
vsplit= vsplit3D
|