-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathscramblepoints.cpp
158 lines (124 loc) · 3.59 KB
/
scramblepoints.cpp
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
//
// Dale Roberts <[email protected]>
//
#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/io.hpp>
bool randomintzo(void)
{
double u=(rand()+1.0)/static_cast<double>(RAND_MAX+2.0);
//double u=0;
if(u<=0.5)
return(0);
else
return(1);
};
boost::numeric::ublas::vector <bool> getdigits(int n, int m)
{
double nt;
bool ndig;
boost::numeric::ublas::vector <bool> ndigits(m);
nt=((double)n)/(pow((double) 2.0,(double)m));
for(int i=1;i<(m+1);i++)
{
//cout << i << endl;
//cout << 2.0*nt << endl;
ndig= floor(2.0*nt);
//cout << ndig << endl;
nt=2.0*nt-ndig;
//cout << nt << endl;
ndigits[m-i]=ndig;
//cout << ndigits[m-i] << endl;
};
return(ndigits);
};
void getmatrixdimj(boost::numeric::ublas::matrix<int>& mj, boost::numeric::ublas::matrix<bool>& mat,int d,int mact,int m)
{
for(int i=0;i<mact;i++)
{
for(int k=0;k<mact;k++)
{
mj(i,k)=mat(d*m+i,k);
};
};
};
void createscrmatdimj(boost::numeric::ublas::matrix<int>&Ltilde,boost::numeric::ublas::matrix<int>& mj,int d,int mact,int K)
{
//multiply L_j C_j(which is mj), where L_j is K*mact, lower triangular and filled with random 0 ones
boost::numeric::ublas::matrix<int> Ltilde2(K,mact);
for(int i=0;i<mact;i++)
{
for(int k=0;k<i;k++)
{
Ltilde2(k,i)=0;
};
Ltilde2(i,i)=1;
for(int k=i+1;k<K;k++)
{
Ltilde2(k,i)=randomintzo();
};
};
//cout << Ltilde2 << endl;
//now multiply Ltilde with C_j(=mj) and save teh result in Ltilde
for(int i=0;i<K;i++)
{
for(int j=0;j<mact;j++)
{
Ltilde(i+d*K,j)=0;
for(int k=0;k<mact;k++)
{
Ltilde(i+d*K,j)+=Ltilde2(i,k)*mj(k,j);
};
};
};
};
void getejvec(boost::numeric::ublas::vector <bool>& ej,int sact,int K)
{
for(int i=0;i<(sact*K);i++)
{
ej(i)=randomintzo();
};
//cout << ej << endl;
};
void getscrambledpoints(boost::numeric::ublas::matrix<double>& P,boost::numeric::ublas::matrix<bool>& mat,int sact,int mact,int m,int K,int n)
{
boost::numeric::ublas::vector <bool> ndigits(mact); // to contain the result of the matrix vector multiplication
boost::numeric::ublas::vector <int> xdigits(K); // to contain the digits of the point
boost::numeric::ublas::matrix <int> mj(mact,mact); //used to store the matrix of the original point set for a particular dimension
boost::numeric::ublas::vector <bool> ej(sact*K); //used to shift the digits
boost::numeric::ublas::matrix <int> Ltilde(sact*K,mact); //to contain the products
for(int d=0;d<sact;d++)
{
getmatrixdimj(mj, mat,d,mact,m); //mj contains C_j
createscrmatdimj(Ltilde,mj,d,mact,K);
};
getejvec(ej,sact,K); //get the whole vector ej(1:K) for dim 1, ej(K+1:2*k) for dim 2 etc
for(int ni=0;ni<n;ni++)
{
ndigits=getdigits(ni,mact); //contains the b-adic expansion of ni
for(int d=0;d<sact;d++) //for dimension d
{
for(int i=0;i<K;i++) //obtain the K digits for point P(ni,d)
{
xdigits[i]=0;
for(int j=0;j<mact;j++)
{
//cout << mat(d*m+i,j) << endl;
//cout << ndigits[j] << endl;
xdigits[i]+=(Ltilde(d*K+i,j)*ndigits(j));
};
//cout << xdigits[i] << endl;
xdigits[i]=(xdigits(i)+ej(d*K+i))%2;
//cout << xdigits[i] << endl;
//cout << endl;
//cout << xdigits[i] << endl;
};
P(d,ni)=0;
for(int i=0;i<K;i++)
{
//cout << (double) pow(float(2),i+1) << endl;
//cout << (double) xdigits[i] << endl;
P(d,ni)+=(double(xdigits(i))/((double) pow(float(2),i+1)));
};
};
};
};