From b598c747b40f9357a5f38e1691c4e4df78765df3 Mon Sep 17 00:00:00 2001
From: Matthieu Schaller <matthieu.schaller@durham.ac.uk>
Date: Tue, 3 Mar 2015 17:46:55 +0000
Subject: [PATCH] Unbiassed correction is now applied on all 3 axis in all 26
configurations.
---
examples/plot_sorted.py | 73 +++++++++++++++++++++------------------
examples/test_bh_sorted.c | 56 ++++++++++++++++++++----------
2 files changed, 77 insertions(+), 52 deletions(-)
diff --git a/examples/plot_sorted.py b/examples/plot_sorted.py
index 2ba8441..b12159a 100644
--- a/examples/plot_sorted.py
+++ b/examples/plot_sorted.py
@@ -22,7 +22,7 @@ params = {'axes.labelsize': 16,
'figure.subplot.wspace' : 0.26 , # the amount of width reserved for blank space between subplots
'figure.subplot.hspace' : 0.22 , # the amount of height reserved for white space between subplots
- 'lines.markersize' : 2,
+ 'lines.markersize' : 4,
'lines.linewidth' : 2,
# 'axes.formatter.limits' : (-2, 1),
@@ -79,14 +79,14 @@ limit_exact = ( dist_cutoff_ratio - 1. )
print limit_exact
#names = ["side", "edge", "corner"]
-#for orientation in range( 26 ):
-for jjj in range(0,3):
- if jjj == 0:
- orientation = 0
- if jjj == 1:
- orientation = 1
- if jjj == 2:
- orientation = 4
+for orientation in range( 26 ):
+# for jjj in range(0,3):
+# if jjj == 0:
+# orientation = 0
+# if jjj == 1:
+# orientation = 1
+# if jjj == 2:
+# orientation = 4
# Read Quickshed accelerations
data=loadtxt( "interaction_dump_%d.dat"%orientation )
@@ -157,30 +157,30 @@ for jjj in range(0,3):
sim = accx_s[(ortho < ortho_max_fit) & (ortho > ortho_min_fit)]
exact = accx_u[(ortho < ortho_max_fit) & (ortho > ortho_min_fit)]
- sim = sim[xx > cell_edge]
- exact = exact[xx > cell_edge]
- xx = xx[xx > cell_edge]
+ #sim = sim[xx > cell_edge]
+ #exact = exact[xx > cell_edge]
+ #xx = xx[xx > cell_edge]
- sim = sim[abs(cell_edge-xx) > limit_exact]
- exact = exact[abs(cell_edge-xx) > limit_exact]
- xx = xx[abs(cell_edge-xx) > limit_exact]
+ sim = sim[xx > limit_exact]
+ exact = exact[xx > limit_exact]
+ xx = xx[xx > limit_exact]
delta = sim / exact
# parabola_param1,pcov = optimize.curve_fit(parabola, xx, delta, p0=[1,1,1], sigma=None, maxfev=100000, xtol=1e-13, ftol=1e-13)
# #print parabola_param1
cubic_param1,pcov = optimize.curve_fit(cubic, xx, delta, p0=[1,1,1,1], sigma=None, maxfev=100000, xtol=1e-13, ftol=1e-13)
print cubic_param1
- # quartic_param1,pcov = optimize.curve_fit(quartic, xx, delta, p0=[1,1,1,1,1], sigma=None, maxfev=100000, xtol=1e-13, ftol=1e-13)
- # #print quartic_param1
+ #quartic_param1,pcov = optimize.curve_fit(quartic, xx, delta, p0=[1,1,1,1,1], sigma=None, maxfev=100000, xtol=1e-13, ftol=1e-13)
+ #print quartic_param1
# quintic_param1,pcov = optimize.curve_fit(quintic, xx, delta, p0=[1,1,1,1,1,1], sigma=None, maxfev=100000, xtol=1e-13, ftol=1e-13)
# #print quintic_param1
-
+
plot(xx, delta, 'bx')
- plot([0,2*cell_edge], [1,1], 'r')
+ plot([min(dist), max(dist)], [1,1], 'r')
text( 0, 1.01, "axis=( %d %d %d )"%(axis[orientation*3 + 0], axis[orientation*3 + 1], axis[orientation*3 + 2]) , ha='center', backgroundcolor='w', fontsize=14)
- xxx = linspace(cell_edge+limit_exact, 2*cell_edge)
+ xxx = linspace(limit_exact, max(dist))
#plot(xxx, parabola(xxx, parabola_param1[0], parabola_param1[1], parabola_param1[2]), 'k-')
- plot(xxx, cubic(xxx, cubic_param1[0], cubic_param1[1], cubic_param1[2], cubic_param1[3]), 'k--')
+ plot(xxx, cubic(xxx, cubic_param1[0], cubic_param1[1], cubic_param1[2], cubic_param1[3]), 'k-')
#plot(xxx, quartic(xxx, quartic_param1[0], quartic_param1[1], quartic_param1[2], quartic_param1[3], quartic_param1[4]), 'k:')
#plot(xxx, quintic(xxx, quintic_param1[0], quintic_param1[1], quintic_param1[2], quintic_param1[3], quintic_param1[4], quintic_param1[5]), 'k-.')
@@ -189,30 +189,35 @@ for jjj in range(0,3):
sim = accx_s[(ortho < ortho_max_fit) & (ortho > ortho_min_fit)]
exact = accx_u[(ortho < ortho_max_fit) & (ortho > ortho_min_fit)]
- sim = sim[xx < cell_edge]
- exact = exact[xx < cell_edge]
- xx = xx[xx < cell_edge]
+ # sim = sim[xx < cell_edge]
+ # exact = exact[xx < cell_edge]
+ # xx = xx[xx < cell_edge]
- sim = sim[abs(cell_edge-xx) > limit_exact]
- exact = exact[abs(cell_edge-xx) > limit_exact]
- xx = xx[abs(cell_edge-xx) > limit_exact]
+ sim = sim[xx < -limit_exact]
+ exact = exact[xx < -limit_exact]
+ xx = xx[xx < -limit_exact]
delta = sim / exact
# parabola_param2,pcov = optimize.curve_fit(parabola, xx, delta, p0=[1,1,1], sigma=None, maxfev=100000, xtol=1e-13, ftol=1e-13)
# #print parabola_param2
cubic_param2,pcov = optimize.curve_fit(cubic, xx, delta, p0=[1,1,1,1], sigma=None, maxfev=100000, xtol=1e-13, ftol=1e-13)
print cubic_param2
- # quartic_param2,pcov = optimize.curve_fit(quartic, xx, delta, p0=[1,1,1,1,1], sigma=None, maxfev=100000, xtol=1e-13, ftol=1e-13)
- # #print quartic_param2
+ #quartic_param2,pcov = optimize.curve_fit(quartic, xx, delta, p0=[1,1,1,1,1], sigma=None, maxfev=100000, xtol=1e-13, ftol=1e-13)
+ #print quartic_param2
# quintic_param2,pcov = optimize.curve_fit(quintic, xx, delta, p0=[1,1,1,1,1,1], sigma=None, maxfev=100000, xtol=1e-13, ftol=1e-13)
# #print quintic_param2
plot(xx, delta, 'bx')
- xxx = linspace(0., cell_edge-limit_exact)
+ xxx = linspace(min(dist), -limit_exact)
#plot(xxx, parabola(xxx, parabola_param2[0], parabola_param2[1], parabola_param2[2]), 'k-')
- plot(xxx, cubic(xxx, cubic_param2[0], cubic_param2[1], cubic_param2[2], cubic_param2[3]), 'k--')
+ plot(xxx, cubic(xxx, cubic_param2[0], cubic_param2[1], cubic_param2[2], cubic_param2[3]), 'k-')
#plot(xxx, quartic(xxx, quartic_param2[0], quartic_param2[1], quartic_param2[2], quartic_param2[3], quartic_param2[4]), 'k:')
#plot(xxx, quintic(xxx, quintic_param2[0], quintic_param2[1], quintic_param2[2], quintic_param2[3], quintic_param2[4], quintic_param2[5]), 'k-.')
+
+ plot([-limit_exact, -limit_exact], [0.9, 1.1], 'k--')
+ plot([limit_exact, limit_exact], [0.9, 1.1], 'k--')
+
+ ylim(0.99, 1.01)
savefig("fit_%d.png"%orientation)
close()
@@ -220,12 +225,12 @@ for jjj in range(0,3):
# Apply the correction ---------------------------------------
# for i in range(size(x)):
- # if cell_edge-dist[i] > limit_exact:
+ # if dist[i] < -limit_exact:
# accx_s[i] /= cubic(dist[i], cubic_param2[0], cubic_param2[1], cubic_param2[2], cubic_param2[3])
- # if dist[i]-cell_edge > limit_exact:
+ # if dist[i] > limit_exact:
# accx_s[i] /= cubic(dist[i], cubic_param1[0], cubic_param1[1], cubic_param1[2], cubic_param1[3])
-
+
# Build error ------------------------------------------------
inv_acc_tot = 1.0 / sqrt(accx_u**2 + accy_u**2 + accz_u**2)
diff --git a/examples/test_bh_sorted.c b/examples/test_bh_sorted.c
index ab5319a..32a1e9b 100644
--- a/examples/test_bh_sorted.c
+++ b/examples/test_bh_sorted.c
@@ -1166,12 +1166,12 @@ static inline void iact_pair_direct_unsorted(struct cell *ci, struct cell *cj) {
float correction_coefs[6*4] =
{
- 0.04730586, -0.37960746, 0.99371996, 0.15444646,
- -0.04397042, 0.10367151, -0.06130978, 1.00754152,
- 0.00536257, -0.06370355, 0.19954405, 0.81828777,
- -0.00677179, -0.01237365, 0.02703505, 0.99483005,
- 0.00701166, -0.0113404, -0.0159732, 1.02330933,
- 0.00297062, 0.01928414, -0.01929454, 1.00190972
+ 0.04730587, -0.13379902, 0.10447393, 0.982606,
+ -0.04397043, -0.12480554, -0.09791498, 0.98388736,
+ 0.00536261, -0.04095222, 0.05153867, 0.98824618,
+ -0.00677177, -0.04110387, -0.04859363, 0.98916257,
+ 0.00425108, 0.01519535, -0.02025942, 1.0035729,
+ -0.00711434, 0.01114372, 0.01865343, 1.00334988
};
@@ -1216,16 +1216,16 @@ static inline void gravity_bias_correction( float d, float axis[3], float corr[3
case 1: /* side */
{
- if ( d > limit_exact)
- for( k = 0; k < 3; k++ )
- corr[k] = 1./cubic( d, correction_coefs[4*4+0], correction_coefs[4*4+1], correction_coefs[4*4+2], correction_coefs[4*4+3]);
-
- else if ( d < -limit_exact)
+
+
+ if ( d < -limit_exact )
for( k = 0; k < 3; k++ )
- corr[k] = 1./cubic( d, correction_coefs[5*4+0], correction_coefs[5*4+1], correction_coefs[5*4+2], correction_coefs[5*4+3]);
+ corr[k] = 1./cubic(d, correction_coefs[5*4+0], correction_coefs[5*4+1], correction_coefs[5*4+2], correction_coefs[5*4+3]);
- message("corr= %f %f %f", corr[0], corr[1], corr[2]);
-
+ else if (d > limit_exact)
+ for( k = 0; k < 3; k++ )
+ corr[k] = 1./cubic(d, correction_coefs[4*4+0], correction_coefs[4*4+1], correction_coefs[4*4+2], correction_coefs[4*4+3]);
+
}
break;
@@ -1233,9 +1233,13 @@ static inline void gravity_bias_correction( float d, float axis[3], float corr[3
case 2: /* edge */
{
- if (d > limit_exact || d < -limit_exact)
+ if ( d < -limit_exact )
+ for( k = 0; k < 3; k++ )
+ corr[k] = 1./cubic(d, correction_coefs[3*4+0], correction_coefs[3*4+1], correction_coefs[3*4+2], correction_coefs[3*4+3]);
+
+ else if (d > limit_exact)
for( k = 0; k < 3; k++ )
- corr[k] = 0.;
+ corr[k] = 1./cubic(d, correction_coefs[2*4+0], correction_coefs[2*4+1], correction_coefs[2*4+2], correction_coefs[2*4+3]);
}
break;
@@ -1243,9 +1247,14 @@ static inline void gravity_bias_correction( float d, float axis[3], float corr[3
case 3: /* corner */
{
- if (d > limit_exact || d < -limit_exact)
+ if ( d < -limit_exact )
for( k = 0; k < 3; k++ )
- corr[k] = 0.;
+ corr[k] = 1./cubic(d, correction_coefs[1*4+0], correction_coefs[1*4+1], correction_coefs[1*4+2], correction_coefs[1*4+3]);
+
+ else if (d > limit_exact)
+ for( k = 0; k < 3; k++ )
+ corr[k] = 1./cubic(d, correction_coefs[0*4+0], correction_coefs[0*4+1], correction_coefs[0*4+2], correction_coefs[0*4+3]);
+
}
break;
@@ -1469,10 +1478,16 @@ static inline void iact_pair_direct_sorted_multipole(struct cell *ci,
/* Get the correction */
gravity_bias_correction( d_norm, axis, corr );
+ //corr[0] = 1.;
+ //corr[1] = 1.;
+ //corr[2] = 1.;
+
+
/* And apply it */
parts_i[i].a[0] *= corr[0];
parts_i[i].a[1] *= corr[1];
parts_i[i].a[2] *= corr[2];
+ parts_i[i].d = d_norm;
}
@@ -1486,10 +1501,15 @@ static inline void iact_pair_direct_sorted_multipole(struct cell *ci,
/* Get the correction */
gravity_bias_correction( d_norm, axis, corr );
+ //corr[0] = 1.;
+ //corr[1] = 1.;
+ //corr[2] = 1.;
+
/* And apply it */
parts_j[j].a[0] *= corr[0];
parts_j[j].a[1] *= corr[1];
parts_j[j].a[2] *= corr[2];
+ parts_j[j].d = d_norm;
}
/* Copy the accelerations back to the original particles. */
--
GitLab