Changeset 34756fd in sasmodels

Timestamp:

Dec 3, 2014 5:04:40 PM (10 years ago)

Author:

Paul Kienzle <pkienzle@…>

Branches:

master, core_shell_microgels, costrafo411, magnetic_model, release_v0.94, release_v0.95, ticket-1257-vesicle-product, ticket_1156, ticket_1265_superball, ticket_822_more_unit_tests

Children:

3699587, 8a3e0af

Parents:

34d49af

Message:

add single to double comparison for the new models

Files:

• compare_many.py (modified) (2 diffs)
• sasmodels/models/capped_cylinder.c (modified) (1 diff)

compare_many.py

@@ -9 +9 @@
                      eval_opencl, eval_ctypes, make_data)
 
-def get_stats(target, model_eval, name, pars, data, index, dtype, cutoff):
-    value, _ = model_eval(name, pars, data, dtype='single', cutoff=cutoff)
+def get_stats(target, value, index):
     resid = abs(value-target)[index]
     relerr = resid/target[index]
@@ -45 +44 @@
 
         env = environment()
-        gpu_single = get_stats(target, eval_opencl, name, pars, data, index, 'single', cutoff)
+        gpu_single_value,_ = eval_opencl(name, pars, data, dtype='single', cutoff=cutoff)
+        gpu_single = get_stats(target, gpu_single_value, index)
         if env.has_double:
-            gpu_double = get_stats(target, eval_opencl, name, pars, data, index, 'double', cutoff)
+            gpu_double_value,_ = eval_opencl(name, pars, data, dtype='double', cutoff=cutoff)
+            gpu_double = get_stats(target, gpu_double_value, index)
         else:
             gpu_double = [0]*len(gpu_single)
-        cpu_double = get_stats(target, eval_ctypes, name, pars, data, index, 'double', cutoff)
+        cpu_double_value,_ =  eval_ctypes(name, pars, data, dtype='double', cutoff=cutoff)
+        cpu_double = get_stats(target, cpu_double_value, index)
+        single_double = get_stats(cpu_double_value, gpu_single_value, index)
 
-        values = list(gpu_single) + list(gpu_double) + list(cpu_double) + [v for _,v in sorted(pars.items())]
+        values = (list(gpu_single) + list(gpu_double) + list(cpu_double)
+                  + list(single_double) + [v for _,v in sorted(pars.items())])
         if gpu_single[0] > 5e-5:
             if first:
-                print_column_headers(pars,'GPU single|GPU double|CPU double'.split('|'))
+                print_column_headers(pars,'GPU single|GPU double|CPU double|single/double'.split('|'))
                 first = False
             print(("%d,"%seed)+','.join("%g"%v for v in values))

sasmodels/models/capped_cylinder.c

@@ -52 +52 @@
 {
     // cap radius should never be less than radius when this is called
-    // Note: cap volume = pi hc/6 * (3 a^2 + hc^2), where a is the cylinder
-    // radius and hc is the height of the cap.  Multiply by two for both ends.
-    // So:
-    //      cap V = pi hc (r^2 + hc^2/3)
-    //      cylinder V = pi r^2 L
-    //      V = cylinder V + cap V
+
+    // Note: volume V = 2*V_cap + V_cyl
+    //
+    // V_cyl = pi r_cyl^2 L
+    // V_cap = 1/6 pi h_c (3 r_cyl^2 + h_c^2) = 1/3 pi h_c^2 (3 r_cap - h_c)
+    //
+    // The docs for capped cylinder give the volume as:
+    //    V = pi r^2 L + 2/3 pi (R-h)^2 (2R + h)
+    // where r_cap=R and h = R - h_c.
+    //
+    // The first part is clearly V_cyl.  The second part requires some work:
+    //    (R-h)^2 => h_c^2
+    //    (2R+h) => 2R+ h_c-h_c + h => 2R + (R-h)-hc + h => 3R-h_c
+    // And so:
+    //    2/3 pi (R-h)^2 (2R + h) => 2/3 pi h_c^2 (3 r_cap - h_c)
+    // which is 2 V_cap, using the second form above.
+    //
+    // In this function we are going to use the first form of V_cap
+    //
+    //      V = V_cyl + 2 V_cap
     //        = pi r^2 L + pi hc (r^2 + hc^2/3)
-    //        = pi * (r^2 (L+hc) + hc^3/3)
+    //        = pi (r^2 (L+hc) + hc^3/3)
     const double hc = cap_radius - sqrt(cap_radius*cap_radius - radius*radius);
     return M_PI*(radius*radius*(length+hc) + 0.333333333333333*hc*hc*hc);