if ( lens[i].type == 5 and fabs(im_position.x - lens[i].position.x) <= dx/2. and fabs(im_position.y - lens[i].position.y) <= 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");
}
}
if(skip_image==0){
//checking whether a closest image has already been found
// 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] && nimagesfound[image_id][i]==0; i++) // we search for an observed image not already linked (nimagesfound=0)
{
if(im_dist[i]<im_dist[second_closest_id])
{
second_closest_id=i;
im_index=i; // im_index value changes only if we found a not linked yet image
tim[image_id][im_index]=temp; // if we found an observed and not already linked image, we allocate the theoretical image temp
}
}
nimagesfound[image_id][im_index]++; // increasing the total number of images found (If we find more than 1 theoretical image linked to 1 real image, these theoretical
// images are included in this number)
}
}
thread_found_image=0; // for next iteration
}
}
}
}
//////////////////////////////////////computing the local chi square//////////////////////////////////////
double chiimage;
for( int iter = 0; iter < nimages_strongLensing[source_id]*nimages_strongLensing[source_id]; iter++){
int i=iter/nimages_strongLensing[source_id];
int j=iter % nimages_strongLensing[source_id];
if(i!=j){
// In the current method, we get the source in the source plane by ray tracing image in nimagesfound[i][i]. If we ray trace back,
// we arrive again at the same position and thus the chi2 from it is 0. Thus we do not calculate the chi2 (-> if i!=j)
if(nimagesfound[i][j]>0){
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
*chi += chiimage;
}
else if(nimagesfound[i][j]==0){
// If we do not find a correpsonding image, we add a big value to the chi2 to disfavor the model
*chi += 100.*nimages_strongLensing[source_id];
}
}
}
/*
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];
}
}
void chi_bruteforce_SOA(double *chi, int *error, runmode_param *runmode, const struct Potential_SOA *lens, const struct grid_param *frame, const int *nimages_strongLensing, galaxy *images){
int bid=0;
int tid=0; // index of the thread within the block
if ( fabs(im_position.x - lens[0].position_x[iterator]) <= dx/2. and fabs(im_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");
}
}
if(skip_image==0){
//checking whether a closest image has already been found
// 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] && nimagesfound[image_id][i]==0; i++) // we search for an observed image not already linked (nimagesfound=0)
{
if(im_dist[i]<im_dist[second_closest_id])
{
second_closest_id=i;
im_index=i; // im_index value changes only if we found a not linked yet image
tim[image_id][im_index]=temp; // if we found an observed and not already linked image, we allocate the theoretical image temp
}
}
nimagesfound[image_id][im_index]++; // increasing the total number of images found (If we find more than 1 theoretical image linked to 1 real image, these theoretical
// images are included in this number)
}
}
thread_found_image=0; // for next iteration
}
}
}
}
//////////////////////////////////////computing the local chi square//////////////////////////////////////
double chiimage;
for( int iter = 0; iter < nimages_strongLensing[source_id]*nimages_strongLensing[source_id]; iter++){
int i=iter/nimages_strongLensing[source_id];
int j=iter % nimages_strongLensing[source_id];
if(i!=j){
// In the current method, we get the source in the source plane by ray tracing image in nimagesfound[i][i]. If we ray trace back,
// we arrive again at the same position and thus the chi2 from it is 0. Thus we do not calculate the chi2 (-> if i!=j)
if(nimagesfound[i][j]>0){
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
*chi += chiimage;
}
else if(nimagesfound[i][j]==0){
// If we do not find a correpsonding image, we add a big value to the chi2 to disfavor the model
*chi += 100.*nimages_strongLensing[source_id];
}
}
}
/*
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];
}
}
void chi_bruteforce_SOA_AVX(double *chi, int *error, runmode_param *runmode, const struct Potential_SOA *lens, const struct grid_param *frame, const int *nimages_strongLensing, galaxy *images){
int bid=0;
int tid=0; // index of the thread within the block
//printf(" %d Found something at x %f y %f by bid %d and tid %d\n",image_id, x_pos,y_pos, bid , tid);
//printf("T ax %f ay %f cx %f cy %f Tsour ax %f ay %f cx %f cy %fby bid %d and tid %d\n", Tinf.a.x, Tinf.a.y , Tinf.c.x, Tinf.c.y ,Tinfsource.a.x, Tinfsource.a.y , Tinfsource.c.x, Tinfsource.c.y , bid , tid);
//printf(" %d Source x %f y %f by bid %d and tid %d\n",source_id, sources[source_id].center.x, sources[source_id].center.y , bid , tid);
//printf(" %d Img x %f y %f im_dist %f and im_index %d\n",image_id, im_position.x, im_position.y , im_dist[im_index] , im_index);
}
}
//
int skip_image = 0;
//
if ( thread_found_image == 1 )
{
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(im_position.x - lens[0].position_x[iterator]) <= dx/2. and fabs(im_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");
}
}
#if 1
if(skip_image==0)
{
//checking whether a closest image has already been found
// 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] && nimagesfound[image_id][i]==0; i++) // we search for an observed image not already linked (nimagesfound=0)
{
if(im_dist[i]<im_dist[second_closest_id])
{
second_closest_id=i;
im_index=i; // im_index value changes only if we found a not linked yet image
tim[image_id][im_index]=temp; // if we found an observed and not already linked image, we allocate the theoretical image temp
}
}
nimagesfound[image_id][im_index]++; // increasing the total number of images found (If we find more than 1 theoretical image linked to 1 real image, these theoretical
// images are included in this number)
}
#endif
}
//thread_found_image=0; // for next iteration
#endif
}
}
}
}
//////////////////////////////////////computing the local chi square//////////////////////////////////////
double chiimage;
for( int iter = 0; iter < nimages_strongLensing[source_id]*nimages_strongLensing[source_id]; iter++){
int i=iter/nimages_strongLensing[source_id];
int j=iter%nimages_strongLensing[source_id];
if(i!=j){
// In the current method, we get the source in the source plane by ray tracing image in nimagesfound[i][i]. If we ray trace back,
// we arrive again at the same position and thus the chi2 from it is 0. Thus we do not calculate the chi2 (-> if i!=j)
if(nimagesfound[i][j]>0){
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
*chi += chiimage;
}
else if(nimagesfound[i][j]==0){
// If we do not find a correpsonding image, we add a big value to the chi2 to disfavor the model
*chi += 100.*nimages_strongLensing[source_id];
}
}
}
/*
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];
}
}
/** @brief Tranform a point from image to source plane. Result stored in sourcepoint argument
*
* Tranform a point from image to source plane using lensequation
*
* @param image_point image position
* @param dlsds dls/ds
* @param nhalos number of halos
* @param potential_param gravitational potential information
* @param source_point address where source information will be stored
*
*
*/
void chi_transformImageToSourcePlane(const runmode_param *runmode, const struct point *image_point, double dlsds, const struct Potential *lens, struct point *source_point)
{ // dlsds is the distance between lens and source divided by the distance observer-source
struct point Grad; // gradient
Grad = module_potentialDerivatives_totalGradient(runmode,image_point, lens);