Home
Manual
Packages
Global Index
Keywords
Quick Reference
|
#include "Chris/ftools.i"
/*
defines 2pt and 3pt correlation functions in 1,2 & 3D
its quite slow for 3pt correlation
*/
// #include "randfield.i"
// ui=genrandfield3D(50); u=fft(ui,[-1,-1,-1]).re;
// mp_include,"Yorick/threptcor.i"
// res =twoptcor(u,pp=15,iso=1)
/*
zrange=span(1e-4,5e-1,50);
pp=1500; u=random(3,pp); tt=Dtwoptcor3D(u,zrange,period=0);
PL,tt/(zrange(zcen)^2)/dif(zrange)(2),zrange(zcen),incolor=-(long(random()*7)+5);
zrange=span(1e-4,5e-1,50);
pp=1000; u=random(2,pp);
tt=Dtwoptcor2D(u,zrange,period=1);
PL,tt/(zrange(zcen))//(zrange(dif)(2)),zrange(zcen),incolor=-(long(random()*7)+5);
#endif 1
*/
func invcor (cor,cov)
/* DOCUMENT
SEE ALSO:
*/
{
ymwrite,"~/.tmp.tmp",cor;
pp=popen("math ",1);
write,pp,"<<.init.m\n"; fflush,pp;
write,pp,"<<BSinv.m\n"; fflush,pp;
write,pp,"a=YorickImport[\"~/.tmp.tmp\"];\n";fflush,pp;
write,pp,"a//Dimensions\n";fflush,pp;
write,pp,"Quit[]\n"; fflush,pp;
}
func cov (u,nmax=,wgts=)
/* DOCUMENT
returns covariance matrix of field u
(averaged over the second dim of u)
note that
res(indgen(n1)+n1*indgen(0:n1-1))-res(,rms)^2==0;
SEE ALSO:
*/
{
dims =dimsof(u);
n1=dims(2);
if(nmax) n1=min(nmax,n1);
u=u(1:n1,);
w1=indgen(n1);
if(is_array(wgts)) wgts /=double(sum(wgts));
res= array(0.,n1,n1);
if(is_void(wgts)) um=u(,avg); else um=(u*wgts(-,))(,sum);
u-=um(,-);
for(i1=1;i1<=n1;i1++)
for(j1=i1;j1<=n1;j1++)
{
if(is_void(wgts)) res(i1,j1)= (u(i1,)*u(j1,))(avg); else
res(i1,j1)=(u(i1,)*u(j1,)*wgts)(sum);
res(j1,i1)=res(i1,j1);
}
return res;
}
func redshiftDist (x1,x2)
{
return (x1-x2)/(1.+(x1+x2)/2.);
}
func __tst_fun_wght1D (x,x1,fx,fx1)
{
return (abs(x(,-)-x1(-,))*fx(,-)*fx1(-,))(*);
}
func Dtwoptcor1D (x,&rrange,cross,fx=,fcross=,period=,fun=)
/* DOCUMENT
returns
the 1+xi, the corrected number of pairs found in the 1D cloud of points
dist= gives the oportunity to specify the distance to be used.
period = allows for periodic boundary conditions (warning period *is* the period)
fun is a function of the vector distance and fx which
returns a scalar weight which is applied to all pairs: the weighting is fun(x,x1,fx,fx1)
// random noise
rrange=span(1e-4,1e-1,50);
cc=[]; for(i=1;i<=5;i++){ pp=2^7; tt=Dtwoptcor1D(random(pp),rrange,period=1); grow,cc,[tt]; }pler,cc-1,mean=1
// GRF
require, "Chris/randfield.i"
cc=[]; for(i=1;i<=15;i++){ tt=Dtwoptcor1D(Dgenrandfield1D(2^10),rrange,period=1); grow,cc,[tt]; }
// check that both ways of measuring a correlation are equivalent
cc=dd=[]; for(i=1;i<=15;i++){ rrange=span(1e-4,0.5,2^7); x=Dgenrandfield1D(2^10,u,pp=2^8);tt=Dtwoptcor1D(x,rrange,period=1); grow,cc,[tt]; grow,dd,[twoptcor1D(u/avg(u)-1,ft=1,iso=1)]; } pler,cc-1,mean=1; pler,dd,mean=1,color=-1;
// test weighting
pp=64; tt=Dtwoptcor1D(random(pp),rrange,fun=__tst_fun_wght1D,fx=random_normal(pp))
SEE ALSO:
*/
{
require,"Chris/split.i";
pp=dimsof(x)(0);
if(is_void(rrange)) rrange=span(0.001,0.1,12);
norm =0.5/rrange(dif)(1);
if(is_void(period)) norm *= 1./(1-rrange(zcen) ); // normalise
if (pp> 2000)
{ // splitting not yet implemented
u2=[];
nn=max(long(pp/2000.),1);
w=vsplit1D(x,span(0,1,nn+1));
qq=dimsof(w)(2);
for(i=1;i<=qq;i++)
for(j=1;j<=qq;j++)
{
x1=x(*w(i,1));
cross1=x(*w(j,1));
if(!is_void(fun))
{
fx1 =fx(..,*w(i,1));
if(is_void(cross)) fcross1 =fx(..,*w(j,1)); else fcross1 =fcross(..,*w(j,1));
}
u1=double(Dtwoptcor1D(x1,rrange,
cross1,fx=fx1,fcross=fcross1,period=period,fun=fun));
grow,u2,[u1];
}
return u2(,avg);
}
else // no splitting
{
if(is_void(cross))
{
if(is_void(dist))
{
if(is_void(period)) tt=abs(x(,-)-x(-,));
else tt=abs(min(abs(x(,-)-x(-,)),period-abs(x(,-)-x(-,))));
}
else tt=dist(x(,-),x(-,));
}
else
{
if(is_void(dist))
{
if(is_void(period)) tt=abs(x(,-)-cross(-,));
else tt=abs(min(abs(x(,-)-cross(-,)),period-abs(x(,-)-cross(-,))));
}
else tt=dist(x(,-),cross(-,));
}
if(!is_void(fun))
{
if(is_void(cross)) wght=fun(x,x,fx,fx); else wght=fun(x,cross,fx,fcross);
}
return histo1d(tt(*),rrange,wght=wght)*norm/double(pp*(pp-1));
}
}
func __tst_fun_wght2D (xy,xy2,fxy,fxy2)
{
tt=abs(xy(1,)(-,)-xy2(1,)(,-),
xy(2,)(-,)-xy2(2,)(,-));
tt *= fxy(-,)*fxy2(,-);
return tt(*);
}
func Dtwoptcor2D (xy,&rrange,cross,fxy=,fcross=,period=,fun=)
/* DOCUMENT
returns
the 1+xi, the corrected excess number of pairs found in the 2D cloud of points
dist= gives the oportunity to specify the distance to be used.
period = allows for periodic boundary conditions (warning period *is* the period)
fun = is a function of the vector distance and which
returns a scalar weight which is applied to all pairs: the weighting is fun(xy,xy2,fxy,fxy2))
EXAMPLE
random noise;
cc=[];pp=1000; for(i=1;i<=15;i++){ u=random(2,pp);tt=Dtwoptcor2D(u,rrange,period=1); grow,cc,[tt];}
gaussian field
require, "Chris/randfield.i";
cc=[]; for(i=1;i<=25;i++){ pp=2^10; xy=Dgenrandfield2D(pp); tt=Dtwoptcor2D(transpose(xy),rrange); grow,cc,[tt]; }
pler,cc-1,mean=1;
// this doesn't seem to work fully well at this stage...
ws;cc=dd=[]; for(i=1;i<=50;i++){ rrange=span(1e-4,0.5,2^7); xy=Dgenrandfield2D(2^12,u,pp=2^8);tt=Dtwoptcor2D(transpose(xy),rrange,period=1); grow,cc,[tt]; grow,dd,[twoptcor2D(u/avg(u)-1,ft=1,iso=1)(1:2^7)]; } pler,cc-1,mean=1; pler,dd,mean=1,color=-1;
pp=64; tt=Dtwoptcor2D(random(2,pp),rrange,fun=__tst_fun_wght2D,fxy=random_normal(pp))
SEE ALSO:
*/
{
require,"Chris/split.i";
if(is_void(rrange)) rrange=span(0.001,0.1,12);
norm= 1./pi/2./(rrange(zcen))/(rrange(dif)(2));
if(is_void(period)) norm *=1./(1+(rrange(zcen)^2-4*rrange(zcen))/pi ); // normalise
pp=dimsof(xy)(0);
if (pp> 2^11)
{
u2=[];
// nn=max(long((pp)^(1./2)/20.),1); write,nn;
nn= max(long(sqrt(pp/2^11))+1,1);
x=xy(1,);y=xy(2,);
w=vsplit2D(x,y,span(0,1,nn+1),span(0,1,nn+1));
qq=dimsof(w)(2);
for(i=1;i<=qq;i++)
for(j=1;j<=qq;j++)
{ if( !((i*qq+j) %10) ) write,pr01(100*((i-1)*qq+j)/double(qq^2)), "% done";
xy1=transpose([x(*w(i,1)),y(*w(i,1))]);
cross1=transpose([x(*w(j,1)),y(*w(j,1))]);
if(!is_void(fun))
{
fxy1 =fxy(..,*w(i,1));
if(is_void(cross)) fcross1 =fxy(..,*w(j,1)); else fcross1 =fcross(..,*w(j,1));
}
u1=double(Dtwoptcor2D(xy1,rrange,
cross1,fxy=fxy1,fcross=fcross1,period=period,fun=fun));
grow,u2,[u1];
}
return u2(,avg);
}
else // no splitting
{
if(is_void(cross)){
if(is_void(period)) tt=abs(xy(1,)(-,)-xy(1,)(,-),xy(2,)(-,)-xy(2,)(,-));
else tt=abs(min(abs(xy(1,)(-,)-xy(1,)(,-)),
period-abs(xy(1,)(-,)-xy(1,)(,-))),
min(abs(xy(2,)(-,)-xy(2,)(,-)),
period-abs(xy(2,)(-,)-xy(2,)(,-))));
}
else
{
if(is_void(period)) tt=abs(xy(1,)(-,)-cross(1,)(,-),
xy(2,)(-,)-cross(2,)(,-));
else tt=abs(min(abs(xy(1,)(-,)-cross(1,)(,-)),
period-abs(xy(1,)(-,)-cross(1,)(,-))),
min(abs(xy(2,)(-,)-cross(2,)(,-)),
period-abs(xy(2,)(-,)-cross(2,)(,-))));
}
if(!is_void(fun))
{
if(is_void(cross)) wght=fun(xy,xy,fxy,fxy); else wght=fun(xy,cross,fxy,fcross);
}
return histo1d(tt(*),rrange,wght=wght)*norm/double(pp*(pp-1));
}
}
func __tst_fun_wght3D (xyz,xyz2,fxyz,fxyz2)
{
tt = fxyz(-,)*fxyz2(,-);
return tt(*);
}
func Dtwoptcor3D (xyz,&rrange,cross,fxyz=,fcross=,period=,fun=)
/* DOCUMENT
returns
the 1+xi, the corrected excess number of pairs found in the 2D cloud of points
dist= gives the oportunity to specify the distance to be used.
period = allows for periodic boundary conditions (warning period *is* the period)
fun = is a function of the positions xyz and xyz2 and a (possibly vector) field fxyz and fxyz2 which
returns a scalar weight which is applied to all pairs: the weighting is fun(xyz,xyz2,fxyz,fxyz2)
TODO : deal with log rrange
deal with proper correction for
EXAMPLE
cc=[]; for(i=1;i<=15;i++){ pp=1500; u=random(3,pp);tt=Dtwoptcor3D(u,rrange,period=1); grow,cc,[tt];}
pler,cc,mean=1,color=-1,marker=-1;
// need to account for window;
// might want to consider linear sampling ...
SEE ALSO:
*/
{
require,"Chris/split.i";
if(is_void(rrange)) rrange=span(0.001,0.1,12);
pp=dimsof(xyz)(0);
norm =1./4./pi/(rrange(zcen))^2/(rrange(dif)(2));
if(is_void(period))
norm *=1./(1+2*rrange(zcen)^2/pi-3*rrange(zcen)/2. -rrange(zcen)^3/4./pi ); // normalise
if (pp> 1500)
{
u2=[];
nn=max(long((pp)^(1./3)/3.5),1);
x=xyz(1,);y=xyz(2,); z=xyz(3,);
w=vsplit(x,y,z,span(0,1,nn+1),span(0,1,nn+1),span(0,1,nn+1));
qq=dimsof(w)(2);
for(i=1;i<=qq;i++)
for(j=1;j<=qq;j++) {
xyz1=transpose([x(*w(i,1)),y(*w(i,1)),z(*w(i,1))]);
cross1=transpose([x(*w(j,1)),y(*w(j,1)),z(*w(j,1))]);
if(!is_void(fun))
{
fxyz1 =fxyz(..,*w(i,1));
if(is_void(cross)) fcross1 =fxyz(..,*w(j,1)); else fcross1 =fcross(..,*w(j,1));
}
u1=double(Dtwoptcor3D(xyz1,rrange, cross1,fxyz=fxyz1,fcross=fcross1,period=period,fun=fun));
grow,u2,[u1];
}
return u2(,avg);
}
else
{ // no splitting
if(is_void(cross))
{
if(is_void(period)) tt=abs(xyz(1,)(-,)-xyz(1,)(,-),
xyz(2,)(-,)-xyz(2,)(,-),
xyz(3,)(-,)-xyz(3,)(,-));
else tt=abs(min(abs(xyz(1,)(-,)-xyz(1,)(,-)),
period-abs(xyz(1,)(-,)-xyz(1,)(,-))),
min(abs(xyz(2,)(-,)-xyz(2,)(,-)),
period-abs(xyz(2,)(-,)-xyz(2,)(,-))),
min(abs(xyz(3,)(-,)-xyz(3,)(,-)),
period-abs(xyz(3,)(-,)-xyz(3,)(,-))));
}
else
{
if(is_void(period)) tt=abs(xyz(1,)(-,)-cross(1,)(,-),
xyz(2,)(-,)-cross(2,)(,-),
xyz(3,)(-,)-cross(3,)(,-));
else tt=abs(min(abs(xyz(1,)(-,)-cross(1,)(,-)),
period-abs(xyz(1,)(-,)-cross(1,)(,-))),
min(abs(xyz(2,)(-,)-cross(2,)(,-)),
period-abs(xyz(2,)(-,)-cross(2,)(,-))),
min(abs(xyz(3,)(-,)-cross(3,)(,-)),
period-abs(xyz(3,)(-,)-cross(3,)(,-))));
} // end cross correl
if(!is_void(fun))
{
if(is_void(cross)) wght=fun(xyz,xyz,fxyz,fxyz); else wght=fun(xyz,cross,fxyz,fcross);
}
return histo1d(tt(*),rrange,wght=wght)*norm/double(pp*(pp-1));
}
}
func twoptcor3D (u,&px,pp=,iso=,norm=,ft=,cross=){
/* DOCUMENT compute 2 pt correlation function in 3D
flag pp = defines the number of pts over which to average
EXAMPLE
#include "Chris/randfield.i"
ui=genrandfield3D(64); u=fft(ui,[-1,-1,-1]).re;
res =twoptcor3D(u)
iso= 1 returns the azimuthally averaged corelation function
ft= 1 does the corelation via fft
SEE ALSO: threeptcor
*/
require, "Eric/histo.i";
u=float(u);
if(!is_void(ft))
{
ns=dimsof(u)(0);
x1=(indgen(ns)-1)*2*pi/float(ns);
y1=(indgen(ns)-1)*2*pi/float(ns);
z1=(indgen(ns)-1)*2*pi/float(ns);
r1= abs(x1, y1(-,),z1(-,-,));
if(is_void(cross)) zc=fft(abs2(fft(u,1)),-1).re; else zc=fft(fft(u,1)*fft(u,-1),-1).re;
zc /= dimsof(u)(2)^3;
if(is_void(iso)) return zc/numberof(u);
py= histo2(r1, px, weight=zc, average=1, interp=1,binsize=r1(dif)(min));
px /= 2*pi;
return py/numberof(u);
if(is_void(norm)) return py;
return py/double(py(1));
}
if(is_void(pp)) pp=5;
dims =dimsof(u);
w=array(0,dims);
w(*)= indgen(dims(2)*dims(3)*dims(4));
n1=dims(2);
n2=dims(3);
n3=dims(4);
w1=indgen(n1);
w2=indgen(n2);
w3=indgen(n3);
res= array(0.,pp,pp,pp);
for(i1=1;i1<=pp;i1++)
for(i2=1;i2<=pp;i2++)
for(i3=1;i3<=pp;i3++)
{
a=[i1-1,i2-1,i3-1];
if(is_void(cross)) u1= u(((w1+a(1)-1) % n1)+1,((w2+a(2)-1) % n2)+1,((w3+a(3)-1) % n3)+1); else
u1= cross(((w1+a(1)-1) % n1)+1,((w2+a(2)-1) % n2)+1,((w3+a(3)-1) % n3)+1);
if(!is_void(norm)) { res(i1,i2,i3)= sum((u*u1))/double(sum((u^2)));}
else res(i1,i2,i3)= sum((u*u1));
}
if (!is_void(iso)){
ns=dimsof(res)(0);
x1=(indgen(ns)-1)/float(ns);
y1=(indgen(ns)-1)/float(ns);
z1=(indgen(ns)-1)/float(ns);
r1= abs(x1, y1(-,),z1(-,-,));
py= histo2(r1, px, weight=res, average=1, interp=1,binsize=r1(dif)(min));
return py/numberof(u);
}
return res/numberof(u);
}
func twoptcor2D (u,&px,pp=,iso=,norm=,ft=,cross=){
/* DOCUMENT compute 2 pt correlation function in 2D
flag pp = defines the number of pts over which to average
EXAMPLE
#include "randfield.i"
ui=genrandfield2D(50); u=fft(ui,[-1,-1]).re;
res =twoptcor2D(u)
iso= 1 returns the azimuthally averaged corelation function
ft= 1 does the corelation via fft
SEE ALSO: threeptcor
*/
require, "Eric/histo.i";
u=float(u);
if(dimsof(u)(1)==3)
{
if(is_void(pp)) res=u*0.;
//=== add case iso !
for(i=1;i<=dimsof(u)(0);i++)
{
if(is_void(cross))
res(,,i)=twoptcor2D(u(,,i),ft=ft,pp=pp,iso=iso,norm=norm); else
res(,,i)=twoptcor2D(u(,,i),ft=ft,pp=pp,iso=iso,norm=norm,cross=cross(,,i));
}
return res;
}
if(!is_void(ft))
{
ns=dimsof(u)(0);
x1=(indgen(ns)-1)*2*pi/float(ns);
y1=(indgen(ns)-1)*2*pi/float(ns);
r1= abs(x1, y1(-,));
if(is_void(cross)) zc=fft(abs2(fft(u,1)),-1).re; else zc=fft(fft(u,1)*fft(cross,-1),-1).re;
zc /= dimsof(u)(2)^2;
if(is_void(iso)) return zc;
py= histo2(r1, px, weight=zc, average=1, interp=1,binsize=r1(dif)(min));
px /= 2*pi;
if(is_void(norm)) return py*1./numberof(u);
return py/double(py(1));
}
if(is_void(pp)) pp=5;
dims =dimsof(u);
w=array(0,dims);
w(*)= indgen(dims(2)*dims(3));
n1=dims(2);
n2=dims(3);
w1=indgen(n1);
w2=indgen(n2);
res= array(0.,pp,pp);
for(i1=1;i1<=pp;i1++)
for(i2=1;i2<=pp;i2++)
{
a=[i1-1,i2-1];
if(is_void(cross)) u1= u(((w1+a(1)-1) % n1)+1,((w2+a(2)-1) % n2)+1); else
u1= cross(((w1+a(1)-1) % n1)+1,((w2+a(2)-1) % n2)+1);
if(!is_void(norm)) { res(i1,i2)= sum((u*u1))/double(sum((u^2)));}
else res(i1,i2)= sum((u*u1));
}
if (!is_void(iso)){
ns=dimsof(res)(0);
x1=(indgen(ns)-1)/float(ns);
y1=(indgen(ns)-1)/float(ns);
r1= abs(x1, y1(-,));
py= histo2(r1, px, weight=res, average=1, interp=1,binsize=r1(dif)(min));
if( pp != n1) px *= pp*1./float(n1);
return py/numberof(u);
}
return res/numberof(u);
}
func twoptcor1D (u,&px,pp=,iso=,norm=,ft=,cross=){
/* DOCUMENT compute 2 pt correlation function in 1D
flag pp = defines the number of pts over which to average
EXAMPLE
#include "randfield.i"
ui=genrandfield1D(1024); u=fft(ui,[-1]).re;
res =twoptcor1D(u)
iso= 1 returns the azimuthally averaged corelation function
ft= 1 does the corelation via fft
SEE ALSO: threeptcor
*/
require, "Eric/histo.i";
u=float(u);
if(dimsof(u)(1)==2)
{
if(is_void(pp)) res=u*0.; else res=array(0.,[2,pp,dimsof(u)(0)]);
for(i=1;i<=dimsof(u)(0);i++)
{
if(is_void(cross)) res(,i)=twoptcor1D(u(,i),px,ft=ft,pp=pp,iso=iso,norm=norm);
else res(,i)=twoptcor1D(u(,i),px,ft=ft,pp=pp,iso=iso,norm=norm,cross=cross(,i));
}
return res;
}
if(!is_void(ft))
{
ns=dimsof(u)(0);
if(is_void(cross)) py=double(fft(abs2(fft(u,1)),-1));
else py =double(fft(fft(u,1)*fft(cross,-1),-1));
px=1./ns*indgen(0:ns-1);
if(norm) return 1./double(py(1))*py;
return 1./ns^2*py;
}
if(is_void(pp)) pp=25;
dims =dimsof(u);
//w=array(0,dims);
// w(*)= indgen(dims(2));
n1=dims(2);
w1=indgen(n1);
res= array(0.,pp);
for(i1=1;i1<=pp;i1++)
{
a=[i1-1];
if(is_void(cross)) u1= u(((w1+a(1)-1) % n1)+1); else u1= cross(((w1+a(1)-1) % n1)+1);
if(!is_void(norm)) { res(i1)= sum((u*u1))/double(sum((u^2)));}
else res(i1)= sum((u*u1));
}
if (!is_void(iso)){
ns=dimsof(res)(0);
x1=(indgen(ns)-1)/float(ns);
r1= abs(x1);
py= histo2(r1, px, weight=res, average=1, interp=1,binsize=r1(dif)(min));
py=py(1:-1);px=px(1:-1); //==== uhmm its odd we need to do this ...
if( pp != n1) px *= pp/float(n1);
return py/numberof(u);
}
return res/numberof(u);
}
func genRt (v,qq=,eps=){
/* DOCUMENT for a given 6d hypercube of 3 pts corelations
compute the residual 3 pt cor versus R and theta
EXAMPLE
//==== this can't work !!! n1 is not defined ????
SEE ALSO:
*/
pp =dimsof(v)(0);
x1=x2=x3=x4=x5=x6=indgen(pp);
n1=dimsof(v)(2); n2=dimsof(v)(3); n3=dimsof(v)(4);
n4=dimsof(v)(5); n5=dimsof(v)(6); n6=dimsof(v)(7);
a1 = x1(,-:1:n2)(,,-:1:n3)(,,,-:1:n1)(,,,,-:1:n2)(,,,,,-:1:n3);
a2 = x2(-:1:n1,)(,,-:1:n3)(,,,-:1:n1)(,,,,-:1:n2)(,,,,,-:1:n3);
a3 = x3(-:1:n2,)(-:1:n1,,)(,,,-:1:n1)(,,,,-:1:n2)(,,,,,-:1:n3);
b1 = x4(-:1:n3,)(-:1:n2,,)(-:1:n1,,,)(,,,,-:1:n2)(,,,,,-:1:n3);
b2 = x5(-:1:n1,)(-:1:n3,,)(-:1:n2,,,)(-:1:n1,,,,)(,,,,,-:1:n3);
b3 = x6(-:1:n2,)(-:1:n1,,)(-:1:n3,,,)(-:1:n2,,,,)(-:1:n1,,,,,);
if(is_void(qq)) qq=10;
if(is_void(eps)) eps=1;
RRs= span(1,pp,qq);
ths= span(0,pi/2,qq);
QQ= array(0.,qq,qq);
for(i=1;i<=qq;i++)
for(j=1;j<=qq;j++)
{
RR=RRs(i);
theta=ths(j);
w =(abs(a1^2+a2^2+a3^2-RR^2)<eps)*(abs(b1^2+b2^2+b3^2-RR^2)<eps);
w *= (abs(cos(theta)- (a1*b1+a2*b2+a3*b3)/RR^2)<eps/5.);
if(is_array(where(w))){
QQ(i,j)=v(where(w))(sum)/dimsof(w(*))(0);}
else QQ(i,j)=0.;
}
return QQ;
}
func genRt2D (v,qq=,eps=){
/* DOCUMENT for a given 4d hypercube of 3 pts corelations
compute the residual 3 pt cor versus R and theta
EXAMPLE
SEE ALSO:
*/
pp =dimsof(v)(0);
x1=x2=x3=x4=indgen(pp);
n1=n2=n3=n4=pp;
a1 = x1(,-:1:n2)(,,-:1:n3)(,,,-:1:n4);
a2 = x2(-:1:n1,)(,,-:1:n3)(,,,-:1:n4);
b1 = x3(-:1:n4,)(-:1:n2,,)(-:1:n1,,,);
b2 = x4(,-:1:n3)(-:1:n2,,)(-:1:n1,,,);
if(is_void(qq)) qq=10;
if(is_void(eps)) eps=1;
RRs= span(1,pp,qq);
ths= span(0,pi/2,qq);
QQ= array(0.,qq,qq);
for(i=1;i<=qq;i++)
for(j=1;j<=qq;j++)
{
RR=RRs(i);
theta=ths(j);
w =(abs(a1^2+a2^2-RR^2)<eps)*(abs(b1^2+b2^2-RR^2)<eps);
w *= (abs(cos(theta) - (a1*b1+a2*b2)/RR^2)<eps/5.);
if(is_array(where(w))){
QQ(i,j)=1./dimsof(w(*))(0)*sum(v(where(w)));}
else QQ(i,j)=0.;
}
return QQ;
}
func threeptcor3D (u,&px,pp=,ft=,cross=){
/* DOCUMENT compute 3 pt correlation function in 3D
flag pp = defines the number of pts over which to average
EXAMPLE
#include "randfield.i"
ui=genrandfield3D(50); u=fft(ui,[-1,-1,-1]).re;
res =threeptcor(u)
SEE ALSO:
*/
if(is_void(pp)) pp=2;
if(dimsof(u)(1)==4)
{
dd= dimsof(u);
if(ft) res=array(0.,dd(2),dd(3),dd(4),dd(2),dd(3),dd(4),dd(5));
else res=array(0.,pp,pp,pp,pp,pp,pp,dd(4));
for(i=1;i<=dimsof(u)(0);i++)
res(,,,,,,i)=threeptcor(u(,,,i),ft=ft,pp=pp,cross=cross(,,,i));
return res;
}
if(ft)
{
uk=fft(u,[1,1,1]);
if (! is_array(uk) || (d = dimsof(uk))(1) != 3 || (n = d(2)) != d(3))
error, "bad dimension list for UK";
i = indgen(0:n-1); /* easier to work with zero-based indices then add one
at the end of computation... */
/* compute 6-D indices in UK for 1st, 2nd and 3rd term in bispectrum */
j = (1 + i + n*i(-,)+n*n*i(-,-,)); /* temporary 3D indice in UK */
i1 = j( , , ,-:1:n,-:1:n,-:1:n);
i2 = j(-:1:n,-:1:n,-:1:n, , , );
j = (n + i - i(-,))%n; /* modulo-N distance between 2 coordinates */
i3 = 1 + j( ,-, ) + n*j(-, ,-, ); /* 4-D index in UK */
/* bispectrum */
u1 = uk(i1);
u2 = uk(i2);
if(is_void(cross)) u3 = uk(i3); else u3=fft(cross,[1,1,1])(i3);
bsp = u1*u2*u3; /* bispectrum */
return double(fft(bsp,[-1,-1,-1,-1,-1,-1]));
}
dims =dimsof(u);
w=array(0,dims);
w(*)= indgen(dims(2)*dims(3)*dims(4));
n1=dims(2);
n2=dims(3);
n3=dims(4);
w1=indgen(n1);
w2=indgen(n2);
w3=indgen(n3);
res= array(0.,pp,pp,pp,pp,pp,pp);
for(i1=1;i1<=pp;i1++)
for(i2=1;i2<=pp;i2++)
for(i3=1;i3<=pp;i3++)
for(j1=1;j1<=pp;j1++)
for(j2=1;j2<=pp;j2++)
for(j3=1;j3<=pp;j3++)
{
a=[i1-1,i2-1,i3-1];
b=[j1-1,j2-1,j3-1];
u1= u(((w1+a(1)-1) % n1)+1,((w2+a(2)-1) % n2)+1,((w3+a(3)-1) % n3)+1);
if(is_void(cross)) u2= u(((w1+b(1)-1) % n1)+1,((w2+b(2)-1) % n2)+1,((w3+b(3)-1) % n3)+1);
else u2= cross(((w1+b(1)-1) % n1)+1,((w2+b(2)-1) % n2)+1,((w3+b(3)-1) % n3)+1);
res(i1,i2,i3,j1,j2,j3)= sum((u*u1*u2))/sum((u^3));
}
return res;
}
func threeptcor2D (u,&px,pp=,ft=,cross=){
/* DOCUMENT compute 3 pt correlation function in 2D
flag pp = defines the number of pts over which to average
EXAMPLE
#include "randfield.i"
ui=genrandfield2D(20); u=fft(ui,[-1,-1]).re;
res =threeptcor2D(u)
SEE ALSO:
*/
if(is_void(pp)) pp=2;
if(dimsof(u)(1)==3)
{
dd= dimsof(u);
if(ft) res=array(0.,dd(2),dd(3),dd(2),dd(3),dd(4));
else res=array(0.,pp,pp,pp,pp,dd(4));
for(i=1;i<=dimsof(u)(0);i++)
res(,,,,i)=threeptcor2D(u(,,i),ft=ft,pp=pp,cross=cross(,,i));
return res;
}
if(ft)
{
uk=fft(u,[1,1]);
if (! is_array(uk) || (d = dimsof(uk))(1) != 2 || (n = d(2)) != d(3))
error, "bad dimension list for UK";
i = indgen(0:n-1); /* easier to work with zero-based indices then add one
at the end of computation... */
/* compute 4-D indices in UK for 1st, 2nd and 3rd term in bispectrum */
j = (1 + i + n*i(-,)); /* temporary 2D indice in UK */
i1 = j( , ,-:1:n,-:1:n);
i2 = j(-:1:n,-:1:n, , );
j = (n + i - i(-,))%n; /* modulo-N distance between 2 coordinates */
i3 = 1 + j( ,-, ) + n*j(-, ,-, ); /* 4-D index in UK */
/* bispectrum */
u1 = uk(i1);
u2 = uk(i2);
if(is_void(cross)) u3 = uk(i3); else u3= cross(i3);
bsp = u1*u2*u3; /* bispectrum */
return double(fft(bsp,[-1,-1,-1,-1]));
}
dims =dimsof(u);
w=array(0,dims);
w(*)= indgen(dims(2)*dims(3));
n1=dims(2);
n2=dims(3);
w1=indgen(n1);
w2=indgen(n2);
res= array(0.,pp,pp,pp,pp);
for(i1=1;i1<=pp;i1++)
for(i2=1;i2<=pp;i2++)
for(j1=1;j1<=pp;j1++)
for(j2=1;j2<=pp;j2++)
{
a=[i1-1,i2-1];
b=[j1-1,j2-1];
u1= u(((w1+a(1)-1) % n1)+1,((w2+a(2)-1) % n2)+1)
if(is_void(cross)) u2= u(((w1+b(1)-1) % n1)+1,((w2+b(2)-1) % n2)+1); else
u2= cross(((w1+b(1)-1) % n1)+1,((w2+b(2)-1) % n2)+1);
res(i1,i2,j1,j2)= sum((u*u1*u2))/sum((u^3));
}
return res;
}
func threeptcor1D (u,&px,pp=,ft=,cross=){
/* DOCUMENT compute 3 pt correlation function in 1D
flag pp = defines the number of pts over which to average
EXAMPLE
#include "randfield.i"
ui=genrandfield1D(20); u=fft(ui,[-1]).re;
res =threeptcor1D(u)
SEE ALSO:
*/
if(is_void(pp)) pp=2;
if(numberof(pp)==1) pp=indgen(pp);
if(ft) pp=dimsof(u)(2);
dims =dimsof(u);
if(dims(1)==2)
{
dd= dimsof(u);
res=array(0.,[3,pp,pp,dd(3)]);
for(i=1;i<=dimsof(u)(0);i++)
res(,,i)=threeptcor1D(u(,i),pp=pp,ft=ft);
return res;
}
if(ft)
{
n=numberof(u);
uk=fft(u,1);
// k1=indgen(n)(,-);
// k2=indgen(n)(-,);
// k3=indgen(0:n-1);k3=1 + (k3 + k3(-,))%n;
// b=uk*uk(-,)*conj(uk(k3));
(k3=indgen(n:1:-1))(1)=0; k3=( (k3+k3(-,)) %n )+1;
if(is_void(cross)) b=uk*uk(-,)*uk(k3); else b=uk*uk(-,)*(fft(cross,1))(k3);
res=double(fft(b,[-1,-1]));
return res;
}
w=array(0,dims);
w(*)= indgen(dims(2));
n1=dims(2);
w1=indgen(n1);
res= array(0.,numberof(pp),numberof(pp));
for(i1=1;i1<=numberof(pp);i1++)
for(j1=1;j1<=numberof(pp);j1++)
{
a=[pp(i1)-1];
b=[pp(j1)-1];
u1= u(((w1+a(1)-1) % n1)+1);
if(is_void(cross)) u2= u(((w1+b(1)-1) % n1)+1); else u2= cross(((w1+b(1)-1) % n1)+1);
if(is_void(norm)) res(i1,j1)= sum((u*u1*u2));
else res(i1,j1)= sum((u*u1*u2))/double(sum((u^3)));
}
return res;
}
func fourptcor1D (u,&px,pp=,ft=,cross=){
/* DOCUMENT compute 4 pt correlation function in 1D
flag pp = defines the number of pts over which to average
EXAMPLE
#include "randfield.i"
ui=genrandfield1D(20); u=fft(ui,[-1]).re;
res =fourptcor1D(u)
SEE ALSO:
*/
if(is_void(pp)) pp=2;
if(numberof(pp)==1) pp=indgen(pp);
if(ft) pp=dimsof(u)(2);
dims =dimsof(u);
if(dims(1)==2)
{
dd= dimsof(u);
res=array(0.,[4,numberof(pp),numberof(pp),numberof(pp),dd(3)]);
for(i=1;i<=dimsof(u)(0);i++)
res(,,,i)=fourptcor1D(u(,i),pp=pp,ft=ft,cross=cross(,i));
return res;
}
if(ft)
{
n=numberof(u);
uk=fft(u,1);
(k3=indgen(n:1:-1))(1)=0; k3=( (k3+k3(-,)+k3(-,-,)) %n )+1;
if(is_void(cross)) b=uk*uk(-,)*uk(-,-,)*uk(k3); else b=uk*uk(-,)*uk(-,-,)*(fft(cross,1))(k3);
res=double(fft(b,[-1,-1,-1]));
return res;
}
w=array(0,dims);
w(*)= indgen(dims(2));
n1=dims(2);
w1=indgen(n1);
res= array(0.,numberof(pp),numberof(pp),numberof(pp));
for(i1=1;i1<=numberof(pp);i1++)
for(j1=1;j1<=numberof(pp);j1++)
for(k1=1;k1<=numberof(pp);k1++)
{
a=[pp(i1)-1];
b=[pp(j1)-1];
c=[pp(k1)-1];
u1= u(((w1+a(1)-1) % n1)+1);
u2= u(((w1+b(1)-1) % n1)+1);
if(is_void(cross)) u3= u(((w1+c(1)-1) % n1)+1); else u3= cross(((w1+c(1)-1) % n1)+1);
if(is_void(norm)) res(i1,j1,k1)= sum((u*u1*u2*u3));
else res(i1,j1,k1)= sum((u*u1*u2*u3))/double(sum((u^4)));
}
return res;
}
func mp_threeptcor (u,pp=,ntrips=){
/* DOCUMENT compute 3 pt correlation function in 3D in //
flag pp = defines the number of pts over which to average
EXAMPLE
#include "Chris/randfield.i"
ui=genrandfield3D(25); u=fft(ui,[-1,-1,-1]).re;
res =mp_threeptcor(u)
SEE ALSO:
*/
if(is_void(pp)) pp=2;
if (is_void(ntrips)) ntrips= 2;
if( is_void(u)) error,"the field must be defined !";
save,createb(".crap.dat"),u;
n1=n2=n3=pp;
x1=x2=x3=x4=x5=x6=indgen(pp);
t1 = x1(,-:1:n2)(,,-:1:n3)(,,,-:1:n1)(,,,,-:1:n2)(,,,,,-:1:n3);
t2 = x2(-:1:n1,)(,,-:1:n3)(,,,-:1:n1)(,,,,-:1:n2)(,,,,,-:1:n3);
t3 = x3(-:1:n2,)(-:1:n1,,)(,,,-:1:n1)(,,,,-:1:n2)(,,,,,-:1:n3);
t4 = x4(-:1:n3,)(-:1:n2,,)(-:1:n1,,,)(,,,,-:1:n2)(,,,,,-:1:n3);
t5 = x5(-:1:n1,)(-:1:n3,,)(-:1:n2,,,)(-:1:n1,,,,)(,,,,,-:1:n3);
t6 = x6(-:1:n2,)(-:1:n1,,)(-:1:n3,,,)(-:1:n2,,,,)(-:1:n1,,,,,);
x1 = t1(*);
x2 = t2(*);
x3 = t3(*);
x4 = t4(*);
x5 = t5(*);
x6 = t6(*);
result= array(0.,pp^6);
njobs=dimsof(x1)(0);
ntasks= mp_partition(njobs, ntrips, 0);
mp_pool, ntasks, threeptcor_sow, threeptcor_work, threeptcor_reap;
result= reform(result,[6,pp,pp,pp,pp,pp,pp]);
return result;
}
func threeptcor_sow (to, i)
{
range= mp_prange(i, ntasks, njobs);
mp_send, to, x1(range);
mp_send, to, x2(range);
mp_send, to, x3(range);
mp_send, to, x4(range);
mp_send, to, x5(range);
mp_send, to, x6(range);
}
func threeptcor_work
{
if( is_void(u)) upload,".crap.dat",s=1;
rgi1= mp_recv();
rgi2= mp_recv();
rgi3= mp_recv();
rgi4= mp_recv();
rgi5= mp_recv();
rgi6= mp_recv();
len =dimsof(rgi1)(0);
res= rgi1*0.;
dims =dimsof(u);
n1=dims(2);
n2=dims(3);
n3=dims(4);
w1=indgen(n1);
w2=indgen(n2);
w3=indgen(n3);
for(i=1;i<=len;i++)
{
i1=rgi1(i);
i2=rgi2(i);
i3=rgi3(i);
j1=rgi4(i);
j2=rgi5(i);
j3=rgi6(i);
a=[i1-1,i2-1,i3-1];
b=[j1-1,j2-1,j3-1];
u1= u(((w1+a(1)-1) % n1)+1,((w2+a(2)-1) % n2)+1,((w3+a(3)-1) % n3)+1);
u2= u(((w1+b(1)-1) % n1)+1,((w2+b(2)-1) % n2)+1,((w3+b(3)-1) % n3)+1);
res(i)= sum((u*u1*u2))/sum((u^3));
}
mp_send, 0, res;
mp_send, 0, 1;
}
func threeptcor_reap (i, m)
{
if (m==1) {
msg= mp_recv();
rng=mp_prange(i, ntasks, njobs);
result(rng)=msg;
} else if (m==2) {
if (mp_recv()!=1) error, "garbled transmission from int_work";
}
return (m==2);
}
|