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 and get a circle, but I don't know how to do this in Python.
Or is there any other way to generate it in Python?
Thanks!
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.
You need to use 2 basis vectors for this.
Generate Random unit N-D vector U
so just create array of N random numbers ui = <-1,+1> and then normalize to unit size by dividing them all by the size sqrt(u0^2+u1^2+...) Beware The U must not be zero !!!
Generate N-D vector V which is perpendicular to U
In 2D and 3D is this simple as you can just swap x,y and 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 + ... = 0
So you can create some elements as random and fit the rest (so the sum is zero). To equalize you can use abs higher vi where ui is abs lower... After this normalize V to unit size.
In the C++ example bellow I used this approach:
v[i]=< -1 , +1 > such that partial dot product is lovering in magnitude (so just compute sign so u[i]*v[i] is negative to partial dot product)u[i] and recompute v[i] so dot product is zeroV to unit sizeBeware also V must not be zero !!!
Circle in N-D
use simple parametric equation:
P(t) = center + radius*cos(t)*U + radius*sin(t)*V
Where t = < 0 , 2*Pi >
Here simple C++ example of generating U,V:
//---------------------------------------------------------------------------
const int n=5; // dimension count
//---------------------------------------------------------------------------
double nd_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;
}
//---------------------------------------------------------------------------
void nd_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;
}
//---------------------------------------------------------------------------
double nd_length(double *a) // = |a|
{
int i; double c;
for (c=0.0,i=0;i<n;i++) c+=a[i]*a[i];
return sqrt(c);
}
//---------------------------------------------------------------------------
void nd_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;
}
//---------------------------------------------------------------------------
void nd_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) -> 0
if (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) = 0
if (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;
const int n=5;
//---------------------------------------------------------------------------
double nd_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;
}
//---------------------------------------------------------------------------
void nd_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;
}
//---------------------------------------------------------------------------
double nd_length(double *a) // = |a|
{
int i; double c;
for (c=0.0,i=0;i<n;i++) c+=a[i]*a[i];
return sqrt(c);
}
//---------------------------------------------------------------------------
void nd_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;
}
//---------------------------------------------------------------------------
void nd_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) -> 0
if (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) = 0
if (fabs(b[i])>1e-10)
{
c-=a[i]*b[i];
a[i]=-c/b[i];
break;
}
nd_unit(a);
}
//---------------------------------------------------------------------------
void nd_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 cases
if (n<3) { for (i=0;i<n;i++) a[i]=0.0; return; }
// init
for (i=0;i<n;i++) a[i]=1.0;
ab = nd_dot(a,b);
ac = nd_dot(a,c);
// tune until dot products near zero
for (e=0;(e<1000)&&(fabs(ab)+fabs(ac)>=1e-5);e++) // max 1000 iterations so it will not hang up
for (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];
else if (fabs(b[i])<fabs(c[i])) a[i]=-ac/c[i];
else if (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=new int*[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 circles
for (i=0;i<100;i++)
{
// 3 vector form
nd_rnd(w); // W random unit vector (normal of the plane)
nd_perpendicular(v,w); // V perpendicular to W
nd_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 numbers
if (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.
1st vector (normal)2nd vector (1 dot product to optimize)3th vector (2 dot products to optimize)4th vector (3 dot products to optimize)n-th vector (n-1 dot products to optimize)2 (except the first) vectors as basis (you can select them randomly)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:
//---------------------------------------------------------------------------
void nd_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 b
for (;;)
{
nd_rnd(tmp);
dot=nd_dot(b,tmp);
if (fabs(dot)<=len) break;
}
// make a unit and perpendicular to b
for (i=0;i<n;i++) a[i]=tmp[i]-(dot*b[i]);
nd_unit(a);
}
//---------------------------------------------------------------------------
void nd_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,...
const int M=10; // vi operands max count
int i,j,k,e,m; double dot[M],*v[M]={v0,v1,v2,v3,v4,v5,v6,v7,v8,v9},q;
// spec cases
if (n<3) { for (i=0;i<n;i++) a[i]=0.0; return; }
// init
for (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 zero
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; }
// k = index of abs max dot
for (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 vec
for (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 test
double 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])));
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):
As a complete python implementation :
N=3 # dimensions
R=1 # radius
npoints=1000
from 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()
{var%20f='http://v.t.sina.com.cn/share/share.php?appkey=1515056452',u=z||d.location,p=['&url=',e(u),'&title=',e(t||d.title),'&source=',e(r),'&sourceUrl=',e(l),'&content=',c||'gb2312','&pic=',e(p||'')].join('');function%20a(){if(!window.open([f,p].join(''),'mb',['toolbar=0,status=0,resizable=1,width=440,height=430,left=',(s.width-440)/2,',top=',(s.height-430)/2].join('')))u.href=[f,p].join('');};if(/Firefox/.test(navigator.userAgent))setTimeout(a,0);else%20a();})(screen,document,encodeURIComponent,'','','https://i.stack.imgur.com/uB1EQ.png', '推荐 贪生不怕死 的文章《How to generate a randomly oriented, high dimensio》','https://www.manongdao.com/article-2089961.html','页面编码gb2312|utf-8默认gb2312'));){kind=link}