Page MenuHomec4science
Files F103425825

chi_computation.cpp
No OneTemporary
Actions

File Metadata

Created
Sat, Mar 1, 23:45

chi_computation.cpp
View Options

//
#include <string.h>
#include "structure_hpc.hpp"
#include "delense_CPU_utils.hpp"
#ifdef __WITH_MPI
#include <mpi.h>
#include "mpi_check.h"
#endif
#define MAX(a,b) (((a)>(b))?(a):(b))
void
chi_computation (double* chi, int* images_found, int* numimagesfound, struct point* imagesposition, const int MAXIMPERSOURCE, const int *nimages_strongLensing,
const galaxy* images, const int nsets)
{
/*
int MAXIMPERSOURCE;
int numImages = 0;
int maxImagesPerSource = 0;
for( int source_id = 0; source_id < nsets; source_id ++)
{
numImages += nimages_strongLensing[source_id];
maxImagesPerSource = MAX(nimages_strongLensing[source_id], maxImagesPerSource);
}
MAXIMPERSOURCE = maxImagesPerSource;
*/
int nimagesfound[nsets][MAXIMPERSOURCE]; // number of images found from the theoretical sources
int prediction_assigned[MAXIMPERSOURCE]; // if prediction has already been assigned
struct point tim[MAXIMPERSOURCE]; // theoretical images (computed from sources)
int index = 0;
//
for (int source_id = 0; source_id < nsets; source_id++)
{
unsigned short int nimages = nimages_strongLensing[source_id];
//
//______________________Initialisation______________________
//
//for (unsigned short int i = 0; i < nimages; ++i)
// nimagesfound[source_id][i] = 0;
memset (&nimagesfound[source_id][0], 0, nimages * sizeof(int));
memset (&prediction_assigned[0], 0, MAXIMPERSOURCE * sizeof(int));
//
//@todo if(numimagesfound[source_id] > 0) ---> check against no img needed
// Creating image_distance matrix
double image_dist[MAXIMPERSOURCE][MAXIMPERSOURCE]; //[constraints][predictions]
struct point prediction_position;
int image_index;
//looping through predictions
for (int ii = 0; ii < numimagesfound[source_id]; ++ii)
{
(*images_found)++;
prediction_assigned[ii] = 0;
//
prediction_position = imagesposition[source_id * MAXIMPERSOURCE + ii];
image_dist[0][ii] = mychi_dist (prediction_position, images[index + 0].center);
printf ("\nsource = %d n images found = %d, n images = %d,\n --> pos x y = %f %f\n", source_id, numimagesfound[source_id],
nimages_strongLensing[source_id], prediction_position.x, prediction_position.y);
fflush (stdout);
//looping through constraints
printf (" ix iy = %f %f, dist = %f\n", images[index + 0].center.x, images[index + 0].center.y, image_dist[0][ii]);
for (int jj = 1; jj < nimages; jj++)
{
// get the distance to each real image and keep the index of the closest real image
image_dist[jj][ii] = mychi_dist (prediction_position, images[index + jj].center);
printf (" ix iy = %f %f, dist = %f\n", images[index + jj].center.x, images[index + jj].center.y, image_dist[jj][ii]);
}
}
//assigning preditions to contraints
//looping through constraints
for (int jj = 0; jj < nimages; jj++)
{
int prediction_index = -1;
//find first non assigned prediction
for (int ii = 0; ii < numimagesfound[source_id]; ++ii)
{
printf(" cont %d predic %d : prediction_index %d %d\n", jj, ii, prediction_index,prediction_assigned[ii]);
if (prediction_assigned[ii] == 0)
{
prediction_index = ii;
prediction_assigned[ii] = 1;
break;
}
}
//If all predictions have not already been assigned
if (prediction_index != -1)
{
//looping through predictions
for (int ii = 1; ii < numimagesfound[source_id]; ++ii)
{
printf(" cont %d predic %d : prediction_index %d %d %f %f\n", jj, ii, prediction_index,prediction_assigned[ii] , image_dist[jj][ii], image_dist[jj][prediction_index]);
if (image_dist[jj][ii] < image_dist[jj][prediction_index] and prediction_assigned[ii] == 0)
{
/* Uncertain about if that should be done ...
//checking if prediction is indeed the closest // wrong should be if not already assigned
bool none_closer = true;
for (int kk = 0; kk < nimages; kk++)
{
printf(" cont %d predic %d %d : %f %f\n", jj, ii, kk, image_dist[jj][ii], image_dist[kk][ii]);
if (image_dist[jj][ii] > image_dist[kk][ii])
{
none_closer = false;
}
}
if(none_closer == true){
prediction_index = ii ;
prediction_assigned[ii] +=1;
}*/
prediction_assigned[prediction_index] = 0;
prediction_index = ii;
prediction_assigned[ii] = 1;
}
}
//write position
tim[jj] = imagesposition[source_id * MAXIMPERSOURCE + prediction_index];
nimagesfound[source_id][jj]++;
printf (" --> tim%d = %f %f\n", jj, tim[jj].x, tim[jj].y);
}
}
#if 0
struct point image_position;
int image_index;
//
for (int ii = 0; ii < /*loc*/numimagesfound[source_id]; ++ii)
{
//
//
int image_index = 0;
//
image_position = imagesposition[source_id*MAXIMPERSOURCE + ii];
image_dist[0] = mychi_dist(image_position, images[index + 0].center);// get the distance to the real image
printf("\nsource = %d n images found = %d, n images = %d,\n --> pos x y = %f %f\n", source_id, numimagesfound[source_id], nimages_strongLensing[source_id], image_position.x, image_position.y); fflush(stdout);
printf(" ix iy = %f %f, dist = %f\n", images[index].center.x, images[index].center.y, mychi_dist(image_position, images[index].center));
for(int i = 1; i < nimages_strongLensing[source_id]; i++)
{
// get the distance to each real image and keep the index of the closest real image
image_dist[i] = mychi_dist(image_position, images[index + i].center);
printf(" ix iy = %f %f, dist = %f\n", images[index+i].center.x, images[index+i].center.y, image_dist[i]);
if (image_dist[i] < image_dist[image_index])
{
image_index = i;
}
}
//@@printf(" --> image_index = %d\n", image_index); fflush(stdout);
//
// we should exit loops here
//
// p1_time += myseconds();
//
int skip_image = 0;
// Sometimes due to the numerical errors at the centerpoint,
// for SIE potentials an additional image will appear at the center of the Potential.
// This is due to the fact that it is not possible to simulate an infinity value
// at the center correctly, Check that sis correspond to Nlens[0]
/*
for (int iterator = 0; iterator < runmode->Nlens[0]; ++iterator)
{
if ( fabs(image_position.x - lens[0].position_x[iterator]) <= dx/2. and fabs(image_position.y - lens[0].position_y[iterator]) <= dx/2.)
{
skip_image = 1;
printf("WARNING: You are using SIE potentials. An image to close to one of the potential centers has been classified as numerical error and removed \n");
}
}*/
//printf("%d %d %d %d\n", x_id, y_id,thread_found_image, skip_image);
if (!skip_image)
{
//#pragma omp atomic
(*images_found)++;
struct point temp;
printf(" nimagesfound %d\n", source_id, image_index, nimagesfound[source_id][image_index]);
//checking whether a closest image has already been found
if (nimagesfound[source_id][image_index] == 0)
{ // if no image found up to now
//image position is allocated to theoretical image
//#pragma omp critical
tim[image_index] = image_position;
//#pragma omp atomic
nimagesfound[source_id][image_index]++;
}
else if (nimagesfound[source_id][image_index] > 0)
{ // if we have already found an image
// If the new image we found is closer than the previous image
printf(" tim2: %f %f\n", image_dist[image_index], mychi_dist(images[index + image_index].center, tim[image_index]));
if (image_dist[image_index] < mychi_dist(images[index + image_index].center, tim[image_index]))
{
temp = tim[image_index]; // we store the position of the old image in temp
//#pragma omp critical
tim[image_index] = image_position;// we link the observed image with the image we just found
printf("on remplace la vieille img %d = %f %f\n", image_index, temp.x, temp.y);
//il faut recalculer les distances pr la vieille image
image_dist[0] = mychi_dist(temp, images[index + 0].center);// get the distance to the real image
printf(" ix iy = %f %f, dist = %f\n", images[index].center.x, images[index].center.y, image_dist[0]);
for(int i = 1; i < nimages_strongLensing[source_id]; i++)
{
// get the distance to each real image and keep the index of the closest real image
image_dist[i] = mychi_dist(temp, images[index + i].center);
printf(" ix iy = %f %f, dist = %f\n", images[index+i].center.x, images[index+i].center.y, image_dist[i]);
}
}
else
{
printf("on garde la vieille img %d = %f %f\n", image_index, tim[image_index].x, tim[image_index].y);
temp = image_position; // we store the position of the image we just found in temp
}
// initialising second_closest_id to the highest value
// Loop over all images in the set except the closest one
// and initialize to the furthest away image
int second_closest_id = 0;
///////////////////////////////////////////////////////////////
// Loop over all images in the set that are not yet allocated to a theoretical image
// and allocate the closest one
for(int i = 0; i < nimages_strongLensing[source_id]; i++)
{
printf("second_closest_id %d de %f %f %d %d %d %f %f\n", second_closest_id, temp.x, temp.y, i, nimages_strongLensing[source_id],image_index,image_dist[i] , image_dist[second_closest_id]);
if(image_dist[i] <= image_dist[second_closest_id] and i != image_index)
{
second_closest_id = i;
// im_index value changes only if we found a not linked yet image
image_index = i;
//#pragma omp critical
tim[image_index] = temp;// if we found an observed and not already linked image, we allocate the theoretical image temp
}
}
// increasing the total number of images found (If we find more than 1 theoretical image linked to 1 real image,
// these theoretical
//#pragma omp atomic
nimagesfound[source_id][image_index]++;
// images are included in this number)
}
}
//#pragma omp atomic
//loc_images_found++;
//thread_found_image = 0; // for next iteration
}
#endif
//
//}
//#pragma omp barrier
//____________________________ end of image loop
//
//____________________________ computing the local chi square
//
double chiimage;
//
int _nimages = nimages_strongLensing[source_id];
//
for (int iter = 0; iter < _nimages; iter++)
{
int i = iter;
//int i=iter/nimages_strongLensing[source_id];
//int j=iter%nimages_strongLensing[source_id];
//if(i != j)
//{
if (nimagesfound[source_id][i] > 0)
{
double pow1 = images[index + i].center.x - tim[i].x;
double pow2 = images[index + i].center.y - tim[i].y;
//
//chiimage = pow(images[index + j].center.x - tim[i][j].x, 2) + pow(images[index + j].center.y - tim[i][j].y, 2); // compute the chi2
chiimage = pow1 * pow1 + pow2 * pow2; // compute the chi2
//printf("%d %d = %.15f\n", i, j, chiimage);
(*chi) += chiimage;
printf (" -> chi = %f\n", *chi);
}
else
//if(nimagesfound[source_id][i] == 0)
{
// If we do not find a correpsonding image, we add a big value to the chi2 to disfavor the model
(*chi) += 1000. * nimages_strongLensing[source_id];
printf (" source_id = %d, image number = %d has not found an image, chi = %f\n", source_id, i, *chi);
}
//}
}
//
//____________________________ end of computing the local chi square
//
//printf("%d: chi = %.15f\n", source_id, *chi);
/*
for (int i=0; i < nimages_strongLensing[source_id]; ++i){
for (int j=0; j < nimages_strongLensing[source_id]; ++j){
printf(" %d",nimagesfound[i][j]);
}
printf("\n");
}*/
//Incrementing Index: Images already treated by previous source_id
index += nimages_strongLensing[source_id];
}
}

Event Timeline

c4science · Help