How To Generate A Randomly Oriented, High Dimension Circle In Python?
I want to generate a randomly oriented circle in R^n. I have been able to successfully generate points on the surface of an n-sphere. I read that you can intersect it with a plane
Solution 1:
Since you already know how to generate a random point on the surface of an n-sphere, just generate two such points, call them P1 and P2. These will determine the plane in which the circle will lie.
I am assuming that both these points are a distance of 1 from the origin, (which will be true if you picked 1 as the radius of your n-sphere). If not, divide each point by its length.
Now we want to make P2 perpendicular to P1. This can be done by
P2 = P2 - dot(P1, P2) P1
P2 = P2 / || P2 ||
Then for each point, you can generate a random angle between 0 and 2pi. Convert this angle to a point on the circle by:
cos(angle) * P1 + sin(angle) * P2.
Solution 2:
You need to use 2 basis vectors for this.
Generate Random unit N-D vector
Uso just create array of N random numbers
ui = <-1,+1>and then normalize to unit size by dividing them all by the sizesqrt(u0^2+u1^2+...)Beware TheUmust not be zero !!!Generate N-D vector
Vwhich is perpendicular toUIn 2D and 3D is this simple as you can just swap
x,yand negate one in 2D or use cross product in 3D. Unfortunately The cross product is not defined in N-D (at least not to my knowledge) So you need to generate non zero vector for which dot product is zero:(U.V) = dot(U,V) = 0 u0.v0 + u1.v1 . u2.v2 + ... = 0So you can create some elements as random and fit the rest (so the sum is zero). To equalize you can use abs higher
viwhereuiis abs lower... After this normalizeVto unit size.In the C++ example bellow I used this approach:
- compute random
v[i]=< -1 , +1 >such that partial dot product is lovering in magnitude (so just compute sign sou[i]*v[i]is negative to partial dot product) - find any non zero
u[i]and recomputev[i]so dot product is zero - normalize
Vto unit size
Beware also
Vmust not be zero !!!- compute random
Circle in N-D
use simple parametric equation:
P(t) = center + radius*cos(t)*U + radius*sin(t)*VWhere
t = < 0 , 2*Pi >
Here simple C++ example of generating U,V:
//---------------------------------------------------------------------------constint n=5; // dimension count//---------------------------------------------------------------------------doublend_dot(double *a,double *b)// = (a.b){
int i; double c;
for (c=0.0,i=0;i<n;i++) c+=a[i]*b[i];
return c;
}
//---------------------------------------------------------------------------voidnd_unit(double *a)// a = a / |a|{
int i; double c;
for (c=0.0,i=0;i<n;i++) c+=a[i]*a[i];
c=sqrt(c); if (c>1e-10) c=1.0/c; else c=0.0;
for (i=0;i<n;i++) a[i]*=c;
}
//---------------------------------------------------------------------------doublend_length(double *a)// = |a|{
int i; double c;
for (c=0.0,i=0;i<n;i++) c+=a[i]*a[i];
returnsqrt(c);
}
//---------------------------------------------------------------------------voidnd_rnd(double *a)// a = (? ? ? ... ?) , |a| = 1.0{
int i; double c;
for (;;)
{
for (i=0;i<n;i++) a[i]=Random()-0.5;
for (c=0.0,i=0;i<n;i++) c+=a[i]*a[i];
if (c>1e-20) break;
}
c=1.0/sqrt(c);
for (i=0;i<n;i++) a[i]*=c;
}
//---------------------------------------------------------------------------voidnd_perpendicular(double *a,double *b)// a = (? ? ? ... ?) , |a| = 1.0 , (a.b) = 0.0{
if (n==1) // 1D case
{
if (fabs(b[0]>1e-10)) a[0]=1.0; else a[0]=1.0;
return;
}
int i; double c1,s;
for (c1=0.0,i=0;i<n;i++) // c1 = dot(a,b)
{
s=1.0; // s = signum of a[i] so (a.b) -> 0if (b[i]<0.0) s=-1.0;
if (c1*s>0.0) s=-s;
a[i]=s*Random();
c1+=a[i]*b[i];
}
for (i=0;i<n;i++) // correct single element so (a.b) = 0if (fabs(b[i])>1e-10)
{
c1-=a[i]*b[i];
a[i]=-c1/b[i];
break;
}
nd_unit(a);
}
//---------------------------------------------------------------------------AnsiString nd_prn(double *a)// this is for printing you can ignore this{
int i; AnsiString s="( ";
for (i=0;i<n;i++) s+=AnsiString().sprintf("%6.3lf ",a[i]);
s+=" )"; return s;
}
//---------------------------------------------------------------------------// This is usage example in VCL you can ignore this:__fastcall TForm1::TForm1(TComponent* Owner):TForm(Owner)
{double u[n],v[n];
Randomize();
for (int i=0;i<10;i++)
{
nd_rnd(u);
nd_perpendicular(v,u);
mm_log->Lines->Add("U = "+nd_prn(u));
mm_log->Lines->Add("V = "+nd_prn(v));
mm_log->Lines->Add("(U.V) = "+AnsiString().sprintf("%6.3lf",nd_dot(u,v)));
}
}
//-------------------------------------------------------------------------And the result of this for confirmation of approach:
U = ( -0.526 -0.623 -0.535 -0.215 -0.044 )
V = ( 0.620 -0.023 -0.287 -0.682 -0.260 )
(U.V) = -0.000
U = ( 0.444 0.282 0.517 -0.203 0.644 )
V = ( 0.072 -0.564 -0.499 0.134 0.640 )
(U.V) = -0.000
U = ( 0.636 0.339 0.357 -0.022 -0.594 )
V = ( -0.559 -0.015 -0.108 -0.497 -0.655 )
(U.V) = 0.000
U = ( -0.626 -0.195 0.491 -0.374 -0.435 )
V = ( -0.722 -0.063 -0.541 0.039 0.424 )
(U.V) = 0.000
U = ( -0.532 -0.481 0.467 -0.517 0.019 )
V = ( 0.536 -0.716 -0.344 -0.205 -0.199 )
(U.V) = -0.000
U = ( 0.696 -0.588 0.018 -0.078 0.405 )
V = ( -0.106 -0.514 -0.645 -0.079 -0.550 )
(U.V) = 0.000
U = ( 0.072 0.679 0.124 -0.204 -0.690 )
V = ( 0.995 -0.058 0.041 0.060 0.037 )
(U.V) = -0.000
U = ( -0.742 -0.283 0.579 -0.150 0.109 )
V = ( -0.043 -0.798 -0.512 -0.308 -0.066 )
(U.V) = -0.000
U = ( 0.606 0.389 -0.551 0.041 -0.420 )
V = ( -0.457 -0.153 -0.489 -0.691 -0.228 )
(U.V) = 0.000
U = ( 0.947 -0.225 0.156 -0.075 0.151 )
V = ( 0.125 -0.153 -0.043 -0.034 -0.979 )
(U.V) = -0.000
[Edit1] randomness issue
As @RoobieNuby pointed out method above can have randomness issues (the planes will not have uniform distribution) After some testing I confirmed it:
So I tested 3 and 2 basis vectors form and the difference can be visually seen. On the right is the above method in 5D (truncated to 2D view) and on the left is new form obtained like this:
//$$---- Form CPP ----//---------------------------------------------------------------------------#include<vcl.h>#include<math.h>#pragma hdrstop#include"Unit1.h"//---------------------------------------------------------------------------#pragma package(smart_init)#pragma resource "*.dfm"
TForm1 *Form1;
Graphics::TBitmap *bmp=NULL;
int xs=0,ys=0,**pyx=NULL;
constint n=5;
//---------------------------------------------------------------------------doublend_dot(double *a,double *b)// = (a.b){
int i; double c;
for (c=0.0,i=0;i<n;i++) c+=a[i]*b[i];
return c;
}
//---------------------------------------------------------------------------voidnd_unit(double *a)// a = a / |a|{
int i; double c;
for (c=0.0,i=0;i<n;i++) c+=a[i]*a[i];
c=sqrt(c); if (c>1e-10) c=1.0/c; else c=0.0;
for (i=0;i<n;i++) a[i]*=c;
}
//---------------------------------------------------------------------------doublend_length(double *a)// = |a|{
int i; double c;
for (c=0.0,i=0;i<n;i++) c+=a[i]*a[i];
returnsqrt(c);
}
//---------------------------------------------------------------------------voidnd_rnd(double *a)// a = (? ? ? ... ?) , |a| = 1.0{
int i; double c;
for (;;)
{
for (i=0;i<n;i++) a[i]=Random()-0.5;
for (c=0.0,i=0;i<n;i++) c+=a[i]*a[i];
if (c>1e-20) break;
}
c=1.0/sqrt(c);
for (i=0;i<n;i++) a[i]*=c;
}
//---------------------------------------------------------------------------voidnd_perpendicular(double *a,double *b)// a = (? ? ? ... ?) , |a| = 1.0 , (a.b) = 0.0{
if (n==1) // 1D case
{
if (fabs(b[0]>1e-10)) a[0]=1.0; else a[0]=1.0;
return;
}
int i; double c,s;
for (c=0.0,i=0;i<n;i++) // c = dot(a,b)
{
s=1.0; // s = signum of a[i] so (a.b) -> 0if (b[i]<0.0) s=-1.0;
if (c*s>0.0) s=-s;
a[i]=s*Random();
c+=a[i]*b[i];
}
for (i=0;i<n;i++) // correct single element so (a.b) = 0if (fabs(b[i])>1e-10)
{
c-=a[i]*b[i];
a[i]=-c/b[i];
break;
}
nd_unit(a);
}
//---------------------------------------------------------------------------voidnd_perpendicular(double *a,double *b,double *c)// a = (? ? ? ... ?) , |a| = 1.0 , (a.b) = 0.0 , (a.c) = 0.0{ // this is not in-place so (&a != &b) and (&a != &c)int i,e; double ab,ac;
// spec casesif (n<3) { for (i=0;i<n;i++) a[i]=0.0; return; }
// initfor (i=0;i<n;i++) a[i]=1.0;
ab = nd_dot(a,b);
ac = nd_dot(a,c);
// tune until dot products near zerofor (e=0;(e<1000)&&(fabs(ab)+fabs(ac)>=1e-5);e++) // max 1000 iterations so it will not hang upfor (i=0;i<n;i++)
{
ab-=a[i]*b[i];
ac-=a[i]*c[i];
if (fabs(b[i])>fabs(c[i])) a[i]=-ab/b[i];
elseif (fabs(b[i])<fabs(c[i])) a[i]=-ac/c[i];
elseif (fabs(ab)>=fabs(ac)) a[i]=-ab/b[i];
else a[i]=-ac/c[i];
ab+=a[i]*b[i];
ac+=a[i]*c[i];
}
nd_unit(a);
}
//---------------------------------------------------------------------------AnsiString nd_prn(double *a){
int i; AnsiString s="( ";
for (i=0;i<n;i++) s+=AnsiString().sprintf("%6.3lf ",a[i]);
s+=" )"; return s;
}
//---------------------------------------------------------------------------// VCL application init (can ignore this)__fastcall TForm1::TForm1(TComponent* Owner):TForm(Owner)
{
bmp=new Graphics::TBitmap;
bmp->HandleType=bmDIB;
bmp->PixelFormat=pf32bit;
if (bmp==NULL) Application->Terminate();
}
//-------------------------------------------------------------------------// VCL application resize event (can ignore this)void __fastcall TForm1::FormResize(TObject *Sender){
if (bmp==NULL) return;
xs=ClientWidth-mm_log->Width;
ys=ClientHeight;
bmp->SetSize(xs,ys);
xs=bmp->Width;
ys=bmp->Height;
if (pyx) delete[] pyx;
pyx=newint*[ys];
for (int y=0;y<ys;y++) pyx[y]=(int*)bmp->ScanLine[y];
Paint();
}
//---------------------------------------------------------------------------// VCL application exit (can ignore this)void __fastcall TForm1::FormDestroy(TObject *Sender){
if (bmp) delete bmp; bmp=NULL;
if (pyx) delete[] pyx; pyx=NULL;
}
//---------------------------------------------------------------------------// VCL application repaint (consider this void main()...)void __fastcall TForm1::FormPaint(TObject *Sender){
if (bmp==NULL) return;
if (pyx==NULL) return;
int i,e,x,y;
double u[n],v[n],w[n];
double a,da,x0,y0,r,c,s;
Randomize();
// clear screen (xs,ys is resolution, pyx is direct pixel access to bitmap bmp)for (y=0;y<ys;y++)
for (x=0;x<xs;x++)
pyx[y][x]=0x00000000;
da=1.0*M_PI/180.0; // angle step
x0=0.5*xs; // center
y0=0.5*ys;
r=150.0; // radius// 100 random circlesfor (i=0;i<100;i++)
{
// 3 vector formnd_rnd(w); // W random unit vector (normal of the plane)nd_perpendicular(v,w); // V perpendicular to Wnd_perpendicular(u,v,w); // U perpendicular to V and W// old 2 vector form// nd_rnd(u);// nd_perpendicular(v,u);for (e=1,a=0.0;e;a+=da) // circle points loop
{
if (a>=2.0*M_PI) { e=0; a=0.0; }
c=r*cos(a);
s=r*sin(a);
x=double(x0+(c*u[0])+(s*v[0])); // use just x,y for render
y=double(y0+(c*u[1])+(s*v[1]));
if ((x>=0)&&(x<xs)&&(y>=0)&&(y<ys)) // if inside screen
pyx[y][x]=0x00FFFFFF; // render it
}
// debug info log to see the U,V,W as numbersif (i==0)
{
mm_log->Text="";
mm_log->Lines->Add("U = "+nd_prn(u));
mm_log->Lines->Add("V = "+nd_prn(v));
mm_log->Lines->Add("W = "+nd_prn(w));
mm_log->Lines->Add("(U.V) = "+AnsiString().sprintf("%6.3lf",nd_dot(u,v)));
mm_log->Lines->Add("(U.W) = "+AnsiString().sprintf("%6.3lf",nd_dot(u,w)));
mm_log->Lines->Add("(W.V) = "+AnsiString().sprintf("%6.3lf",nd_dot(w,v)));
}
}
// refresh Application window with bitmap (bacbuffering)
Canvas->Draw(0,0,bmp);
}
//---------------------------------------------------------------------------// VCLo mouse wheel event ... just force repaint on mouse wheel to check the randomness in time (can ignore this)void __fastcall TForm1::FormMouseWheel(TObject *Sender, TShiftState Shift, int WheelDelta, TPoint &MousePos, bool &Handled){
Paint();
}
//---------------------------------------------------------------------------The important stuff here is void nd_perpendicular(double *a,double *b,double *c); function which returns perpendicular vector to both b,c. To make this robust (In case randomness matters) in N-D You should have n vectors not just 3.
The operation of function is similar to the previous one. The difference is only in that we optimize more that one dot product at once. So you need to choose which dot product to optimize per axis (based on abs coordinate value of the perpendicular vector). So you should rewrite this to optimize n-1 dot products at once.
- create
1stvector (normal) - create
2ndvector (1dot product to optimize) - create
3thvector (2dot products to optimize) - create
4thvector (3dot products to optimize) - ...
- create
n-thvector (n-1dot products to optimize) - use any
2(except the first) vectors as basis (you can select them randomly) - produce your circle
It looks like just 3 vectors are enough for n>=3 (so you can ignore #4,#5) but to confirm that for sure some sort of statistical analysis would be needed which I am too lazy to do.
[Edit2] vector perpendicular to more vectors ...
Have implemented the dot product perpendicular vector computation from the oher answer. And computation of perpendicular vector to multiple vectors:
//---------------------------------------------------------------------------voidnd_perpendicular(double *a,double *b)// a = (? ? ? ... ?) , |a| = 1.0 , (a.b) = 0.0{
int i; double tmp[n],dot,len;
// |dot| limit to consider as non parallel
len=0.7*nd_length(b);
// tmp = random unit vector not parallel to bfor (;;)
{
nd_rnd(tmp);
dot=nd_dot(b,tmp);
if (fabs(dot)<=len) break;
}
// make a unit and perpendicular to bfor (i=0;i<n;i++) a[i]=tmp[i]-(dot*b[i]);
nd_unit(a);
}
//---------------------------------------------------------------------------voidnd_perpendicular(double *a,double *v0,double *v1,double *v2=NULL,double *v3=NULL,double *v4=NULL,double *v5=NULL,double *v6=NULL,double *v7=NULL,double *v8=NULL,double *v9=NULL){
// this is not in-place so (&a != &b) and (&a != &c)// a = unit vector perpendicular to all not NULL vectors v0,v1,v2,v3,...constint M=10; // vi operands max countint i,j,k,e,m; double dot[M],*v[M]={v0,v1,v2,v3,v4,v5,v6,v7,v8,v9},q;
// spec casesif (n<3) { for (i=0;i<n;i++) a[i]=0.0; return; }
// initfor (i=0;i<n;i++) a[i]=1.0;
nd_rnd(a);
for (m=0,j=0;j<M;j++) if (v[j]) { m=j+1; dot[j]=nd_dot(a,v[j]); } else dot[j]=0.0;
// tune until dot products near zerofor (e=0;e<1000;e++) // max 1000 iterations so it will not hang up
{
for (q=0.0,j=0;j<m;j++) q+=fabs(dot[j]);
if (q<1e-3) { e=-1; break; }
// k = index of abs max dotfor (k=0,j=0;j<m;j++) if (fabs(dot[j])>=fabs(dot[k])) k=j;
// i = index of abs max coordinate of v[k]for (i=0,j=0;j<n;j++) if (fabs(v[k][j])>=fabs(v[k][i])) i=j;
// update dots and a[i]for (j=0;j<m;j++) dot[j]-=a[i]*v[j][i];
a[i]=-dot[k]/v[k][i];
for (j=0;j<m;j++) dot[j]+=a[i]*v[j][i];
}
if (e>=0)
for (e=0;e<1000;e++) // max 1000 iterations so it will not hang up
{
for (q=0.0,j=0;j<m;j++) q+=fabs(dot[j]);
if (q<1e-3) { e=-1; break; }
for (i=0;i<n;i++)
{
// k = index of abs max vecfor (k=0,j=0;j<m;j++) if (fabs(v[j][i])>=fabs(v[k][i])) k=j;
// update dots and a[i]for (j=0;j<m;j++) dot[j]-=a[i]*v[j][i];
a[i]=-dot[k]/v[k][i];
for (j=0;j<m;j++) dot[j]+=a[i]*v[j][i];
}
}
nd_unit(a);
}
//---------------------------------------------------------------------------Too boost robustness It uses 2 different iterative approaches. Here example of usage (n=5):
// nd_perpendicular testdouble V[n][n]; int x,y;
nd_rnd(V[0]);
nd_perpendicular(V[1],V[0]);
nd_perpendicular(V[2],V[1],V[0]);
nd_perpendicular(V[3],V[2],V[1],V[0]);
nd_perpendicular(V[4],V[3],V[2],V[1],V[0]);
for (x=0;x<n;x++)
for (y=x+1;y<n;y++)
mm_log->Lines->Add(AnsiString().sprintf("(V%i.V%i) = %6.3lf",x,y,nd_dot(V[x],V[y])));
Solution 3:
The following is assuming that by high dimensional circle you mean the intersection of an N-sphere with an N-hyperplane containing the sphere's centre (i.e. an N-1-sphere in N+1-space).
You can do the following (assuming you are in N+1-space, considering the N sphere):
- pick a normal vector n, your circle will be the subset of the sphere that is orthogonal to n
- use SVD to get an orthonormal basis B for the orthogonal complement of n
- use your method to create points p on the N-1 sphere and multiply them by B (B is N+1 x N, so Bp will be in N+1-space)
Solution 4:
As a complete python implementation :
N=3# dimensions
R=1# radius
npoints=1000from numpy.random import rand
from numpy.linalg import norm
from numpy import dot,sin,cos,outer,pi
V1,V2,origin=rand(3,N)
u1=V1/norm(V1)
u2=V2/norm(V2)
V3=u1-dot(u1,u2)*u2
u3=V3/norm(V3)
print(norm(u2),norm(u3),dot(u2,u3))
#1.0 1.0 0.0#u2,u3 is orthonormed
theta=np.arange(0,2*pi,2*pi/npoints)
circle=origin+R*(outer(cos(theta),u2)+outer(sin(theta),u3))
import matplotlib as mpl
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
fig = plt.figure()
ax = fig.gca(projection='3d')
ax.plot(*circle.T)
plt.show()
Post a Comment for "How To Generate A Randomly Oriented, High Dimension Circle In Python?"