Homec4science
Diffusion Lenstool-HPC (preparemerge) Browse

R1448/src/chi_computation.cpp
3ca39a32c40bpreparemerge

Last Change
Raw File

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);
//looping through constraints
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);
}
}
//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)
{
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)
{
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]++;
}
}
//____________________________ 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;
}
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];
}
}
c4science · Help