Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F120598397
chi_computation.cpp
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Sat, Jul 5, 12:30
Size
5 KB
Mime Type
text/x-c
Expires
Mon, Jul 7, 12:30 (2 d)
Engine
blob
Format
Raw Data
Handle
27162398
Attached To
R1448 Lenstool-HPC
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];
}
}
Event Timeline
Log In to Comment