Page MenuHomec4science
Files F76684533

gradient.cpp
No OneTemporary
Actions

File Metadata

Created
Fri, Aug 9, 20:41

gradient.cpp
View Options

#include <iostream>
#include <string.h>
//#include <cuda_runtime.h>
#include <math.h>
#include <sys/time.h>
#include <fstream>
#include <immintrin.h>
//#include "simd_math.h"
#include "gradient.hpp"
#include "structure.h"
//#include "iacaMarks.h"
//#define __INV RCP
//#define __INV RCP_1NR
//#define __INV RCP_2NR
//#define __SQRT _mm256_sqrt_pd
//#define __SQRT SQRT
//#define __SQRT SQRT_1NR
//#define __SQRT SQRT_2NR
void module_readParameters_calculatePotentialparameter_SOA(Potential_SOA *lens, int ii)
{
//std::cout << "module_readParameters_calculatePotentialparameter..." << std::endl;
switch (lens->type[ii])
{
case(5): /*Elliptical Isothermal Sphere*/
//impact parameter b0
lens->b0[ii] = 4* pi_c2 * lens->vdisp[ii] * lens->vdisp[ii] ;
//ellipticity_potential
lens->ellipticity_potential[ii] = lens->ellipticity[ii]/3 ;
break;
case(8): /* PIEMD */
//impact parameter b0
lens->b0[ii] = 6.*pi_c2 * lens->vdisp[ii] * lens->vdisp[ii];
//ellipticity_parameter
if ( lens->ellipticity[ii] == 0. && lens->ellipticity_potential[ii] != 0. ){
// emass is (a2-b2)/(a2+b2)
lens->ellipticity[ii] = 2.*lens->ellipticity_potential[ii] / (1. + lens->ellipticity_potential[ii] * lens->ellipticity_potential[ii]);
//printf("1 : %f %f \n",lens->ellipticity[ii],lens->ellipticity_potential[ii]);
}
else if ( lens->ellipticity[ii] == 0. && lens->ellipticity_potential[ii] == 0. ){
lens->ellipticity_potential[ii] = 0.00001;
//printf("2 : %f %f \n",lens->ellipticity[ii],lens->ellipticity_potential[ii]);
}
else{
// epot is (a-b)/(a+b)
lens->ellipticity_potential[ii] = (1. - sqrt(1 - lens->ellipticity[ii] * lens->ellipticity[ii])) / lens->ellipticity[ii];
//printf("3 : %f %f \n",lens->ellipticity[ii],lens->ellipticity_potential[ii]);
}
break;
default:
std::cout << "ERROR: LENSPARA profil type of clump "<< lens->name << " unknown : "<< lens->type << std::endl;
//printf( "ERROR: LENSPARA profil type of clump %s unknown : %d\n",lens->name, lens->type);
break;
};
}
//
/**@brief Return the gradient of the projected lens potential for one clump
* !!! You have to multiply by dlsds to obtain the true gradient
* for the expressions, see the papers : JP Kneib & P Natarajan, Cluster Lenses, The Astronomy and Astrophysics Review (2011) for 1 and 2
* and JP Kneib PhD (1993) for 3
*
* @param pImage point where the result is computed in the lens plane
* @param lens mass distribution
*/
//
/// Useful functions
//
//static
complex
piemd_1derivatives_ci05(double x, double y, double eps, double rc)
{
double sqe, cx1, cxro, cyro, rem2;
complex zci, znum, zden, zis, zres;
double norm;
//
//std::cout << "piemd_lderivatives" << std::endl;
sqe = sqrt(eps);
cx1 = (1. - eps) / (1. + eps);
cxro = (1. + eps) * (1. + eps);
cyro = (1. - eps) * (1. - eps);
//
rem2 = x * x / cxro + y * y / cyro;
/*zci=cpx(0.,-0.5*(1.-eps*eps)/sqe);
znum=cpx(cx1*x,(2.*sqe*sqrt(rc*rc+rem2)-y/cx1));
zden=cpx(x,(2.*rc*sqe-y));
zis=pcpx(zci,lncpx(dcpx(znum,zden)));
zres=pcpxflt(zis,b0);*/
// --> optimized code
zci.re = 0;
zci.im = -0.5 * (1. - eps * eps) / sqe;
znum.re = cx1 * x;
znum.im = 2.*sqe * sqrt(rc * rc + rem2) - y / cx1;
zden.re = x;
zden.im = 2.*rc * sqe - y;
norm = zden.re * zden.re + zden.im * zden.im; // zis = znum/zden
zis.re = (znum.re * zden.re + znum.im * zden.im) / norm;
zis.im = (znum.im * zden.re - znum.re * zden.im) / norm;
norm = zis.re;
zis.re = log(sqrt(norm * norm + zis.im * zis.im)); // ln(zis) = ln(|zis|)+i.Arg(zis)
zis.im = atan2(zis.im, norm);
// norm = zis.re;
zres.re = zci.re * zis.re - zci.im * zis.im; // Re( zci*ln(zis) )
zres.im = zci.im * zis.re + zis.im * zci.re; // Im( zci*ln(zis) )
//
//zres.re = zis.re*b0;
//zres.im = zis.im*b0;
//
return(zres);
}
//
//// changes the coordinates of point P into a new basis (rotation of angle theta)
//// y' y x'
//// * | /
//// * | / theta
//// * | /
//// *|--------->x
//
static struct point rotateCoordinateSystem(struct point P, double theta)
{
struct point Q;
Q.x = P.x*cos(theta) + P.y*sin(theta);
Q.y = P.y*cos(theta) - P.x*sin(theta);
return(Q);
}
__attribute__((noinline))
static struct point
grad_halo(const struct point *pImage, const struct Potential *lens)
{
struct point true_coord, true_coord_rotation, result;
double R, angular_deviation;
complex zis;
//
//std::cout << "grad_halo..." << lens->type << std::endl;
result.x = result.y = 0.;
//
/*positionning at the potential center*/
// Change the origin of the coordinate system to the center of the clump
true_coord.x = pImage->x - lens->position.x;
true_coord.y = pImage->y - lens->position.y;
//
switch (lens->type)
{
case(5): /*Elliptical Isothermal Sphere*/
/*rotation of the coordiante axes to match the potential axes*/
true_coord_rotation = rotateCoordinateSystem(true_coord, lens->ellipticity_angle);
//
R=sqrt(true_coord_rotation.x*true_coord_rotation.x*(1 - lens->ellipticity/3.) + true_coord_rotation.y*true_coord_rotation.y*(1 + lens->ellipticity/3.)); //ellippot = ellipmass/3
result.x = (1 - lens->ellipticity/3.)*lens->b0*true_coord_rotation.x/(R);
result.y = (1 + lens->ellipticity/3.)*lens->b0*true_coord_rotation.y/(R);
break;
case(8): /* PIEMD */
/*rotation of the coordiante axes to match the potential axes*/
true_coord_rotation = rotateCoordinateSystem(true_coord, lens->ellipticity_angle);
/*Doing something....*/
zis = piemd_1derivatives_ci05(true_coord_rotation.x, true_coord_rotation.y, lens->ellipticity_potential, lens->rcore);
//
result.x = lens->b0*zis.re;
result.y = lens->b0*zis.im;
break;
default:
std::cout << "ERROR: Grad 1 profil type of clump "<< lens->name << " unknown : "<< lens->type << std::endl;
break;
};
return result;
}
//
//
//
struct point module_potentialDerivatives_totalGradient(const runmode_param *runmode, const struct point *pImage, const struct Potential *lens)
{
struct point grad, clumpgrad;
grad.x = 0;
grad.y = 0;
//std::cout << "nhalos = " << runmode->nhalos << std::endl;
for(int i = 0; i < runmode->nhalos; i++)
{
clumpgrad = grad_halo(pImage, &lens[i]); //compute gradient for each clump separately
//nan check
//std::cout << clumpgrad.x << " " << clumpgrad.y << std::endl;
if(clumpgrad.x == clumpgrad.x or clumpgrad.y == clumpgrad.y)
{
// add the gradients
grad.x += clumpgrad.x;
grad.y += clumpgrad.y;
}
}
//
return(grad);
}
struct point grad_halo_piemd_SOA(const struct point *pImage, int iterator, const struct Potential_SOA *lens)
{
struct point clumpgrad;
clumpgrad.x = 0;
clumpgrad.y = 0;
struct point true_coord, true_coord_rotation; //, result;
//double R, angular_deviation;
complex zis;
//
//result.x = result.y = 0.;
//
true_coord.x = pImage->x - lens->position_x[iterator];
true_coord.y = pImage->y - lens->position_y[iterator];
//
//std::cout << true_coord.x <<" "<< true_coord.y << std::endl;
//
// Change the origin of the coordinate system to the center of the clump
true_coord_rotation = rotateCoordinateSystem(true_coord, lens->ellipticity_angle[iterator]);
// std::cout << i << ": " << true_coord.x << " " << true_coord.y << " -> "
// << true_coord_rotation.x << " " << true_coord_rotation.y << std::endl;
//
//double sqe, cx1, cxro, cyro, rem2;
//
double x = true_coord_rotation.x;
double y = true_coord_rotation.y;
double eps = lens->ellipticity_potential[iterator];
double rc = lens->rcore[iterator];
//
//std::cout << "piemd_lderivatives" << std::endl;
//
double sqe = sqrt(eps);
//
double cx1 = (1. - eps)/(1. + eps);
double cxro = (1. + eps)*(1. + eps);
double cyro = (1. - eps)*(1. - eps);
//
double rem2 = x*x/cxro + y*y/cyro;
//
//zci=cpx(0.,-0.5*(1.-eps*eps)/sqe);
//znum=cpx(cx1*x,(2.*sqe*sqrt(rc*rc+rem2)-y/cx1));
//zden=cpx(x,(2.*rc*sqe-y));
//zis=pcpx(zci,lncpx(dcpx(znum,zden)));
//zres=pcpxflt(zis,b0);
// --> optimized code
complex zci, znum, zden, zres;
double norm;
//
zci.re = 0;
zci.im = -0.5*(1. - eps*eps)/sqe;
//
znum.re = cx1*x;
znum.im = 2.*sqe*sqrt(rc*rc + rem2) - y/cx1;
//
zden.re = x;
zden.im = 2.*rc*sqe - y;
norm = (zden.re*zden.re + zden.im*zden.im); // zis = znum/zden
//
zis.re = (znum.re*zden.re + znum.im*zden.im)/norm;
zis.im = (znum.im*zden.re - znum.re*zden.im)/norm;
norm = zis.re;
zis.re = log(sqrt(norm*norm + zis.im*zis.im)); // ln(zis) = ln(|zis|)+i.Arg(zis)
zis.im = atan2(zis.im, norm);
// norm = zis.re;
zres.re = zci.re*zis.re - zci.im*zis.im; // Re( zci*ln(zis) )
zres.im = zci.im*zis.re + zis.im*zci.re; // Im( zci*ln(zis) )
//
zis.re = zres.re;
zis.im = zres.im;
//
//zres.re = zis.re*b0;
//zres.im = zis.im*b0;
//
//
clumpgrad.x = lens->b0[iterator]*zis.re;
clumpgrad.y = lens->b0[iterator]*zis.im;
//nan check
//
return(clumpgrad);
}
struct point grad_halo_sis_SOA(const struct point *pImage, int iterator, const struct Potential_SOA *lens)
{
struct point true_coord, true_coord_rotation, result;
double R, angular_deviation;
complex zis;
result.x = result.y = 0.;
/*positionning at the potential center*/
true_coord.x = pImage->x - lens->position_x[iterator]; // Change the origin of the coordinate system to the center of the clump
true_coord.y = pImage->y - lens->position_y[iterator];
/*Elliptical Isothermal Sphere*/
/*rotation of the coordiante axes to match the potential axes*/
true_coord_rotation = rotateCoordinateSystem(true_coord, lens->ellipticity_angle[iterator]);
R=sqrt(true_coord_rotation.x*true_coord_rotation.x*(1-lens->ellipticity[iterator]/3.)+true_coord_rotation.y*true_coord_rotation.y*(1+lens->ellipticity[iterator]/3.)); //ellippot = ellipmass/3
result.x=(1-lens->ellipticity[iterator]/3.)*lens->b0[iterator]*true_coord_rotation.x/(R);
result.y=(1+lens->ellipticity[iterator]/3.)*lens->b0[iterator]*true_coord_rotation.y/(R);
return result;
}
struct point module_potentialDerivatives_totalGradient_SOA(const int *Nlens, const struct point *pImage, const struct Potential_SOA *lens)
{
struct point grad, clumpgrad;
grad.x=0;
grad.y=0;
//This here could be done with function pointer to better acomodate future ass distributions functions
// However I'm unsure of the time of function pointers -> ask gilles
//for the moment lens and Nlens is organised the following way : 1. SIS, 2. PIEMD
//SIS is the first
for(int i=0; i<Nlens[0]; i++){
clumpgrad=grad_halo_sis_SOA(pImage,i,&lens[0]); //compute gradient for each clump separately
if(clumpgrad.x == clumpgrad.x or clumpgrad.y == clumpgrad.y){ //nan check
grad.x+=clumpgrad.x;
grad.y+=clumpgrad.y;
} // add the gradients
}
//PIEMD is the second
for(int i=0; i<Nlens[1]; i++){
clumpgrad=grad_halo_piemd_SOA(pImage,i,&lens[1]); //compute gradient for each clump separately
if(clumpgrad.x == clumpgrad.x or clumpgrad.y == clumpgrad.y){ //nan check
grad.x+=clumpgrad.x;
grad.y+=clumpgrad.y;
//printf(" %f %f %f %f %f %f \n", grad.x,grad.y, clumpgrad.x,clumpgrad.y,pImage->x,pImage->y);
} // add the gradients,
}
//printf(" %f %f \n", grad.x,grad.y);
return(grad);
}
struct point module_potentialDerivatives_piemd_peelGradient_SOA(const struct point *pImage, const struct Potential_SOA *lens, int shalos, int nhalos)
{
struct point grad, clumpgrad;
grad.x = 0;
grad.y = 0;
//
for(int i = shalos; i < nhalos; i++)
{
//IACA_START;
//
struct point true_coord, true_coord_rot; //, result;
//double R, angular_deviation;
complex zis;
//
//result.x = result.y = 0.;
//
true_coord.x = pImage->x - lens->position_x[i];
true_coord.y = pImage->y - lens->position_y[i];
/*positionning at the potential center*/
// Change the origin of the coordinate system to the center of the clump
true_coord_rot = rotateCoordinateSystem(true_coord, lens->ellipticity_angle[i]);
//
double x = true_coord_rot.x;
double y = true_coord_rot.y;
double eps = lens->ellipticity_potential[i];
double rc = lens->rcore[i];
//
//std::cout << "piemd_lderivatives" << std::endl;
//
double sqe = sqrt(eps);
//
double cx1 = (1. - eps)/(1. + eps);
double cxro = (1. + eps)*(1. + eps);
double cyro = (1. - eps)*(1. - eps);
//
double rem2 = x*x/cxro + y*y/cyro;
//
complex zci, znum, zden, zres;
double norm;
//
zci.re = 0;
zci.im = -0.5*(1. - eps*eps)/sqe;
//
znum.re = cx1*x;
znum.im = 2.*sqe*sqrt(rc*rc + rem2) - y/cx1;
//
zden.re = x;
zden.im = 2.*rc*sqe - y;
norm = (zden.re*zden.re + zden.im*zden.im); // zis = znum/zden
//
zis.re = (znum.re*zden.re + znum.im*zden.im)/norm;
zis.im = (znum.im*zden.re - znum.re*zden.im)/norm;
norm = zis.re;
zis.re = log(sqrt(norm*norm + zis.im*zis.im)); // ln(zis) = ln(|zis|)+i.Arg(zis)
zis.im = atan2(zis.im, norm);
// norm = zis.re;
zres.re = zci.re*zis.re - zci.im*zis.im; // Re( zci*ln(zis) )
zres.im = zci.im*zis.re + zis.im*zci.re; // Im( zci*ln(zis) )
//
zis.re = zres.re;
zis.im = zres.im;
//
//zres.re = zis.re*b0;
//zres.im = zis.im*b0;
// rotation
clumpgrad.x = zis.re;
clumpgrad.y = zis.im;
clumpgrad = rotateCoordinateSystem(clumpgrad, -lens->ellipticity_angle[i]);
//
clumpgrad.x = lens->b0[i]*clumpgrad.x;
clumpgrad.y = lens->b0[i]*clumpgrad.y;
//
//clumpgrad.x = lens->b0[i]*zis.re;
//clumpgrad.y = lens->b0[i]*zis.im;
//nan check
if(clumpgrad.x == clumpgrad.x or clumpgrad.y == clumpgrad.y)
{
// add the gradients
grad.x += clumpgrad.x;
grad.y += clumpgrad.y;
}
}
//IACA_END;
//
return(grad);
}
struct point module_potentialDerivatives_totalGradient_SOA_GPU(const int *Nlens, const struct point *pImage, const struct Potential_SOA *lens)
{
struct point grad, clumpgrad;
grad.x=0;
grad.y=0;
//This here could be done with function pointer to better acomodate future ass distributions functions
// However I'm unsure of the time of function pointers -> ask gilles
//for the moment lens and Nlens is organised the following way : 1. SIS, 2. PIEMD
//SIS is the first
for(int i=0; i<Nlens[0]; i++){
clumpgrad=grad_halo_sis_SOA(pImage,i,&lens[0]); //compute gradient for each clump separately
if(clumpgrad.x == clumpgrad.x or clumpgrad.y == clumpgrad.y){ //nan check
grad.x+=clumpgrad.x;
grad.y+=clumpgrad.y;
} // add the gradients
}
//PIEMD is the second
clumpgrad=grad_halo_piemd_SOA_GPU(pImage,Nlens[1],&lens[1]); //compute gradient for each clump separately
if(clumpgrad.x == clumpgrad.x or clumpgrad.y == clumpgrad.y){ //nan check
grad.x+=clumpgrad.x;
grad.y+=clumpgrad.y;
}
return(grad);
}

Event Timeline

c4science · Help