#include <cmath>
struct GausND {
GausND( int dim ) :
{}
double operator() (double *x, double *p) {
int k = 0;
for (int i = 0; i<dim; ++i) { X[i] = x[i] - p[k]; k++; }
for (int i = 0; i<dim; ++i) {
CovMat(i,i) = p[k]*p[k];
k++;
}
for (int i = 0; i<dim; ++i) {
for (int j = i+1; j<dim; ++j) {
CovMat(i,j) = p[k]*
sqrt(CovMat(i,i)*CovMat(j,j));
CovMat(j,i) = CovMat(i,j);
k++;
}
}
if (debug) {
}
if (det <= 0) {
Fatal(
"GausND",
"Determinant is <= 0 det = %f",det);
return 0;
}
if (debug) {
std::cout << "det " << det << std::endl;
std::cout << "norm " << norm << std::endl;
std::cout << "fval " << fval << std::endl;
}
return fval;
}
};
void multidimSampling() {
const int N = 10000;
const int DIM = 4;
double xmin[] = {-10,-10,-10, -10};
double xmax[] = { 10, 10, 10, 10};
double par0[] = { 1., -1., 2, 0,
1, 2, 1, 3,
0.5,0.,0.,0.,0.,0.8 };
const int NPAR = DIM + DIM*(DIM+1)/2;
GausND gaus4d(4);
TF1 * f =
new TF1(
"functionND",gaus4d,0,1,14);
double x0[] = {0,0,0,0};
debug = false;
for (int i = 0; i < NPAR; ++i ) {
}
if (sampler == 0) {
Info(
"multidimSampling",
"Default sampler %s is not available try with Foam ",
}
sampler = Factory::CreateDistSampler();
if (sampler == 0) {
Error(
"multidimSampling",
"Foam sampler is not available - exit ");
return;
}
bool ret = sampler->
Init();
std::vector<double> data1(DIM*N);
double v[DIM];
if (!ret) {
Error(
"Sampler::Init",
"Error initializing unuran sampler");
return;
}
for (
int i = 0; i <
N; ++i) {
for (int j = 0; j < DIM; ++j)
data1[N*j + i] = v[j];
}
TFile * file =
new TFile(
"multiDimSampling.root",
"RECREATE");
double x[DIM];
for (
int i = 0; i <
N; ++i) {
for (int j = 0; j < DIM; ++j) {
x[j] = data1[i+N*j];
}
}
t1->
Draw(
"x[0]:x[1]:x[2]:x[3]",
"",
"candle");
int ic=1;
}