[ create a new paste ] login | about

Link: http://codepad.org/HjhNSk2S    [ raw code | output | fork ]

C++, pasted on Sep 12: 
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
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215

#include <vector>
#include <assert.h>
#include <math.h>
#include <string.h>
#include <iostream>
#include <sstream>

#undef WIN32
#define WINAPI 
#define LPVOID void*
#define DWORD int


template<int DIM, typename TYPE> class Vector;
struct TsFlow {
	int from;		// Feature number in signature 1 
	int to;			// Feature number in signature 2 
	double amount;	// Amount of flow from signature1.features[from] to signature2.features[to]
};

// some example interpolation function that you can use for the particles...
// ...if you're working on a Euclidean space in R^{DIM}
template<int DIM, typename SAMPLESTYPE> double sqrDistLinear(const Vector<DIM,SAMPLESTYPE> &F1, const Vector<DIM,SAMPLESTYPE> &F2);
template<int DIM, typename SAMPLESTYPE> double rbfFuncLinear(const Vector<DIM,SAMPLESTYPE> &center, const Vector<DIM,SAMPLESTYPE> &eval, double r2);
template<int DIM, typename SAMPLESTYPE> Vector<DIM,SAMPLESTYPE> interpolateBinsLinear(const Vector<DIM,SAMPLESTYPE>  &a, const Vector<DIM,SAMPLESTYPE>  &b, double alpha);


// forward declarations - you can skip that, they are not even needed with Visual Studio
template<int DIM, typename SAMPLESTYPE> struct InterpModeData;
template<int DIM, typename SAMPLESTYPE> DWORD WINAPI interpMode(LPVOID lpParam );


// and the interpolation class itself
template<int DIM, typename SAMPLESTYPE=double>
class Interpolator
{
public:
	Interpolator(const std::vector<Vector<DIM,SAMPLESTYPE> > &_A, const std::vector<double> &_wA, 
				 const std::vector<Vector<DIM,SAMPLESTYPE> > &_B, const std::vector<double> &_wB, 
				 double (*_sqrDist)(const Vector<DIM,SAMPLESTYPE> &, const Vector<DIM,SAMPLESTYPE> &), 
				 double (*_kernel)(const Vector<DIM,SAMPLESTYPE> &, const Vector<DIM,SAMPLESTYPE> &, double),
				 Vector<DIM,SAMPLESTYPE>  (*_interpolateBins)(const Vector<DIM,SAMPLESTYPE> &, const Vector<DIM,SAMPLESTYPE> &, double alpha), 
				 double _kNNDist=2.0, unsigned int _NLevels=1) :
		A(_A), B(_B),        // the location in space of each sample point
		wA(_wA), wB(_wB),   // the value of each sample point
		sqrDist(_sqrDist),   // the cost of the particle travelling, typically a squared geodesic distance
		kernel(_kernel),     // the kernel, typically gaussian, to represent each particle
		interpolateBins(_interpolateBins), // the function to move particles along the geodesics
		kNNDist(_kNNDist),   // the standard deviation of the particles kernel, expressed in terms of distance to N'th nearest neighbor (1.5 means the distance halfway between the first and second nearest neighbor)
		NLevels(_NLevels),   // the number of frequency bands
		numberOfNN(3)		// number of nearest neighbors to precompute
	{ };

		void interpolate(double alpha, const std::vector<Vector<DIM,SAMPLESTYPE> > &samples, std::vector<double> &wResult)
		{			
			for (unsigned int numMode=0; numMode<resultFlow.size(); numMode++)
			{
				InterpModeData<DIM, SAMPLESTYPE> data(numMode,alpha,resultFlow, interpolateBins, modes, kernel, variancesA, variancesB, A, B, sumWA, sumWB, samples);
				interpMode<DIM,SAMPLESTYPE>((void*) &data);				
			}
		}

	std::vector<std::vector<double> > modes;

private:

	double (*sqrDist)(const Vector<DIM,SAMPLESTYPE>&, const Vector<DIM,SAMPLESTYPE>&);
	double (*kernel)(const Vector<DIM,SAMPLESTYPE>&, const Vector<DIM,SAMPLESTYPE>&, double);	
	Vector<DIM,SAMPLESTYPE>  (*interpolateBins)(const Vector<DIM,SAMPLESTYPE> &, const Vector<DIM,SAMPLESTYPE> &, double alpha);
	std::vector<Vector<DIM,SAMPLESTYPE> > A, B;
	std::vector<double> wA, wB;
	std::vector<double> sumWA, sumWB;
	unsigned int NLevels;

	std::vector<double> variancesA, variancesB;			// the variance actually used for the RBF	
	std::vector<std::vector<int> > allNNidA, allNNidB;  // distance of the k-NN to each sample	

	const unsigned int numberOfNN;

	std::vector<std::vector<TsFlow> > resultFlow;  // EMD flow between A and B
	
	double kNNDist;	
};



/////////////////////////////////////////////////////////////////////
//////////////////////  OTHER CLASSES DEFINITIONS ///////////////////
/////////////////////////////////////////////////////////////////////


template<int DIM, typename TYPE> class Vector
{
public:
	Vector(){for(int i=0; i<DIM; i++) data[i]=(TYPE)0;};
	Vector(TYPE x, TYPE y) {assert(DIM==2); data[0]=x; data[1]=y;};
	Vector(TYPE x, TYPE y, TYPE z) {assert(DIM==2); data[0]=x; data[1]=y; data[2]=z;};
	Vector(const Vector<DIM,TYPE> &rhs){for(int i=0; i<DIM; i++) data[i]=rhs[i];};
	TYPE operator[](int i) const {return data[i];};
	TYPE& operator[](int i) {return data[i];};
	Vector<DIM,TYPE> operator*(double rhs) const
	{
		Vector<DIM,TYPE> result;
		for (int i=0; i<DIM; i++)
			result[i] = data[i]*rhs;
		return result;
	}
	Vector<DIM,TYPE> operator+(const Vector<DIM,TYPE> &rhs) const
	{
		Vector<DIM,TYPE> result;
		for (int i=0; i<DIM; i++)
			result[i] = data[i]+rhs[i];
		return result;
	}

private:
	TYPE data[DIM];
};



////// Datastructures passed to threads

// Data for the nearest neighbor search (for threads)
template<int DIM, typename SAMPLESTYPE>
struct InterpModeData
{
	InterpModeData(int _i, double _alpha, const std::vector<std::vector<TsFlow> > &_resultFlow, Vector<DIM,SAMPLESTYPE>  (*_interpolateBins)(const Vector<DIM,SAMPLESTYPE> &, const Vector<DIM,SAMPLESTYPE> &, double alpha),
		std::vector<std::vector<double> > &_modes, double (*_kernel)(const Vector<DIM,SAMPLESTYPE> &, const Vector<DIM,SAMPLESTYPE> &, double),
		const std::vector<double> &_variancesA, const std::vector<double> &_variancesB,
		const std::vector<Vector<DIM,SAMPLESTYPE> > &_A, const std::vector<Vector<DIM,SAMPLESTYPE> > &_B,
		const std::vector<double> &_sumWA, const std::vector<double> &_sumWB,
		const std::vector<Vector<DIM,SAMPLESTYPE> > &_samples):
	i(_i),alpha(_alpha),resultFlow(_resultFlow), interpolateBins(_interpolateBins),  modes(_modes),
	kernel(_kernel), variancesA(_variancesA), variancesB(_variancesB), A(_A), B(_B), sumWA(_sumWA), sumWB(_sumWB), samples(_samples)
	{
	}

	int i; // mode number
	double alpha;
	const std::vector<std::vector<TsFlow> > &resultFlow ;
	Vector<DIM,SAMPLESTYPE>  (*interpolateBins)(const Vector<DIM,SAMPLESTYPE> &, const Vector<DIM,SAMPLESTYPE> &, double alpha);
	std::vector<std::vector<double> > &modes;
	double (*kernel)(const Vector<DIM,SAMPLESTYPE> &, const Vector<DIM,SAMPLESTYPE> &, double);
	const std::vector<double> &variancesA;
	const std::vector<double> &variancesB;
	const std::vector<Vector<DIM,SAMPLESTYPE> > &A;
	const std::vector<Vector<DIM,SAMPLESTYPE> > &B;
	const std::vector<double> &sumWA;
	const std::vector<double> &sumWB;
	const std::vector<Vector<DIM,SAMPLESTYPE> > &samples;
};


/////////////////////////////////////////////////////////
/////////////////// Interpolation functions /////////////
////////////////////////////////////////////////////////

template<int DIM, typename SAMPLESTYPE>
Vector<DIM,SAMPLESTYPE> interpolateBinsLinear(const Vector<DIM,SAMPLESTYPE>  &a, const Vector<DIM,SAMPLESTYPE>  &b, double alpha)
{
	return a*(1.-alpha) + b*alpha;
}
template<int DIM, typename SAMPLESTYPE>
double sqrDistLinear(const Vector<DIM,SAMPLESTYPE> &F1, const Vector<DIM,SAMPLESTYPE> &F2)
{
	double d = 0;
	for (int i=0; i<DIM; i++)
		d+= (F1[i]-F2[i])*(F1[i]-F2[i]);
	return d;
}
template<int DIM, typename SAMPLESTYPE>
double rbfFuncLinear(const Vector<DIM,SAMPLESTYPE> &center, const Vector<DIM,SAMPLESTYPE> &eval, double r2)
{
	double sqrdist = sqrDistLinear<DIM,SAMPLESTYPE>(center, eval);
	double stddev = sqrt(r2);
	return exp(-sqrdist/(2.*r2))/(pow(2.5066, DIM)*stddev); // 2.5066 = sqrt(2*pi)
}

//--------------------------------------------------------------------------------


template<int DIM, typename SAMPLESTYPE>
DWORD WINAPI interpMode(LPVOID lpParam )
{
	InterpModeData<DIM,SAMPLESTYPE>* argData = (InterpModeData<DIM,SAMPLESTYPE>*) lpParam;
	 //do stuffs
	return 0;
}



int main()
{
	int W=50;

	std::vector<Vector<2,double> > samplesPos(W*W);
	std::vector<double> values1(W*W, 0.);
	std::vector<double> values2(W*W, 0.);
	std::vector<double> valuesR(W*W, 0.);

	Interpolator<2,double> interp(samplesPos, values1,  // source distribution, in R^2 
		samplesPos, values2,							// target distribution
		sqrDistLinear, rbfFuncLinear, interpolateBinsLinear, // how to represent and move the particles
		2,												// particle spread : distance to 2nd nearest neighbor at each sample point
		1);											   // 1 frequency band (standard displacement interpolation)

	int N = 50;
	for (int p=0; p<N; p++)
	{
		double alpha = p/((double)N-1.);
		interp.interpolate(alpha, samplesPos, valuesR);
	}
	return 0;
}


Output: 
1
2
3
4
5
6
7
8
9

cc1plus: warnings being treated as errors
t.cpp: In constructor 'Interpolator<DIM, SAMPLESTYPE>::Interpolator(const __gnu_debug_def::vector<Vector<DIM, SAMPLESTYPE>, std::allocator<Vector<DIM, SAMPLESTYPE> > >&, const __gnu_debug_def::vector<double, std::allocator<double> >&, const __gnu_debug_def::vector<Vector<DIM, SAMPLESTYPE>, std::allocator<Vector<DIM, SAMPLESTYPE> > >&, const __gnu_debug_def::vector<double, std::allocator<double> >&, double (*)(const Vector<DIM, SAMPLESTYPE>&, const Vector<DIM, SAMPLESTYPE>&), double (*)(const Vector<DIM, SAMPLESTYPE>&, const Vector<DIM, SAMPLESTYPE>&, double), Vector<DIM, SAMPLESTYPE> (*)(const Vector<DIM, SAMPLESTYPE>&, const Vector<DIM, SAMPLESTYPE>&, double), double, unsigned int) [with int DIM = 2, SAMPLESTYPE = double]':
t.cpp:207:   instantiated from here
Line 71: warning: 'Interpolator<2, double>::wB' will be initialized after
Line 67: warning:   'double (* Interpolator<2, double>::sqrDist)(const Vector<2, double>&, const Vector<2, double>&)'
Line 43: warning:   when initialized here
Line 82: warning: 'Interpolator<2, double>::kNNDist' will be initialized after
Line 73: warning:   'unsigned int Interpolator<2, double>::NLevels'
Line 43: warning:   when initialized here


Create a new paste based on this one


Comments: