Merge pull request #6 from hz-b/dev/feature/math-comb

[TASK] binom, added test, throwing exception on overflow

python/thor_scsi/utils/accelerator.py

@@ -101,10 +101,11 @@ def extract_orbit_from_standard_observers(
  tps_tmp = [ob.get_truncated_power_series_a() for ob in ob_sub]
  logger.debug(tps_tmp)
  tps = xr.DataArray(
- data=[t.jacobian() for t in tps_tmp],
+ data=np.array([t.jacobian() for t in tps_tmp], float),
  name="tps_jacobian",
  dims=["index", "phase_coordinate_row", "phase_coordinate_col"],
- coords=[indices, phase_space_coords_names, phase_space_coords_names],
+# todo: don't understand where this comes from
+ coords=[indices, phase_space_coords_names, phase_space_coords_names + ["unknown"]],
  )
  ps = xr.DataArray(
  data=[np.array(t.cst().iloc) for t in tps_tmp],

src/thor_scsi/CMakeLists.txt

@@ -179,55 +179,12 @@ target_link_libraries(test_transform_phase_space
 
 add_test(transform_phase_space test_transform_phase_space)
 
-## add_executable(a_kick elements/a_kick.cc)
-## target_include_directories(a_kick
-## PUBLIC
-## "$<BUILD_INTERFACE:${CMAKE_CURRENT_SOURCE_DIR}/../>"
-## "$<INSTALL_INTERFACE:${CMAKE_INSTALL_INCLUDEDIR}>"
-## )
-## target_link_libraries(a_kick
-## thor_scsi
-## thor_scsi_core
-## tpsa_lin
-## # PRIVATE Boost::boost
-## # quadmath
-## ${Boost_PRG_EXEC_MONITOR_LIBRARY}
-## ${Boost_UNIT_TEST_FRAMEWORK_LIBRARY}
-## )
-
-#
-# Non linear part needs still to be ported to new interface
-## set(tpsa_nlin_FILES
-## tpsa_for_pm.cc
-## # tpsa_for.cc
-## )
-##
-## add_library(tpsa_nlin SHARED
-## ${tpsa_common_FILES}
-## ${tpsa_nlin_FILES}
-## ${tpsa_HEADERS}
-## )
-
-
-
-## set(thor_scsi_CORE_HEADERS
-## thor/core/util.h
-## flame/core/base.h
-## flame/core/config.h
-## )
-##
-##
-##
-## add_library(thor_scsi SHARED
-## ${thor_scsi_CORE_files}
-## ${thor_scsi_CORE_HEADERS}
-## )
-
 
 add_subdirectory(elements)
 add_subdirectory(std_machine)
 add_subdirectory(examples)
 add_subdirectory(custom)
+add_subdirectory(core)
 
 # ---- install helpers --
 

src/thor_scsi/core/CMakeLists.txt

@@ -1,23 +1,39 @@
 # could be added to file above
 # but not yet required
 # lets leave it here
+message(STATUS "thor_scsi::core gtpsa cpp include dir ${gtpsa_cpp_INCLUDE_DIR}")
+
+add_executable(test_math_comb
+ test_math_comb.cc
+)
+target_include_directories(test_math_comb
+ PUBLIC
+ "$<BUILD_INTERFACE:${CMAKE_CURRENT_SOURCE_DIR}/../../>"
+)
+target_link_libraries(test_math_comb
+ thor_scsi_core
+ ${Boost_PRG_EXEC_MONITOR_LIBRARY}
+ ${Boost_UNIT_TEST_FRAMEWORK_LIBRARY}
+)
+
+add_test(math_comb test_math_comb)
 
 add_executable(test_transform
- test_transform.cc
+  test_transform.cc
 )
 
-message(STATUS "thor_scsi::core gtpsa cpp include dir ${gtpsa_cpp_INCLUDE_DIR}")
+
 target_include_directories(test_transform
- PUBLIC
- "$<BUILD_INTERFACE:${CMAKE_CURRENT_SOURCE_DIR}/../../>"
- "$<BUILD_INTERFACE:${gtpsa_cpp_INCLUDE_DIR}>"
- "$<BUILD_INTERFACE:${gtpsa_mad_ng_INCLUDE_DIR}>"
- "$<INSTALL_INTERFACE:${CMAKE_INSTALL_INCLUDEDIR}>"
+  PUBLIC
+  "$<BUILD_INTERFACE:${CMAKE_CURRENT_SOURCE_DIR}/../../>"
+  "$<BUILD_INTERFACE:${gtpsa_cpp_INCLUDE_DIR}>"
+  "$<BUILD_INTERFACE:${gtpsa_mad_ng_INCLUDE_DIR}>"
+  "$<INSTALL_INTERFACE:${CMAKE_INSTALL_INCLUDEDIR}>"
 )
 
 target_link_libraries(test_transform
- ${Boost_PRG_EXEC_MONITOR_LIBRARY}
- ${Boost_UNIT_TEST_FRAMEWORK_LIBRARY}
+  ${Boost_PRG_EXEC_MONITOR_LIBRARY}
+  ${Boost_UNIT_TEST_FRAMEWORK_LIBRARY}
 )
 
 add_test(transform test_transform)
@@ -38,12 +54,8 @@ target_include_directories(test_two_dimensional_multipoles
  "$<INSTALL_INTERFACE:${CMAKE_INSTALL_INCLUDEDIR}>"
 )
 target_link_libraries(test_two_dimensional_multipoles
- # PRIVATE Boost::boost
- # quadmath
- # gsl
- # gslcblas
- gtpsa::c++
- gtpsa
- ${Boost_PRG_EXEC_MONITOR_LIBRARY}
- ${Boost_UNIT_TEST_FRAMEWORK_LIBRARY}
+ gtpsa::gtpsa-c++
+ gtpsa::gtpsa
+ ${Boost_PRG_EXEC_MONITOR_LIBRARY}
+ ${Boost_UNIT_TEST_FRAMEWORK_LIBRARY}
 )

src/thor_scsi/core/math_comb.cc

@@ -1,56 +1,60 @@
 #include <thor_scsi/core/math_comb.h>
 #include <limits.h>
-#include <cassert>
+#include <stdexcept>
+#include <sstream>
 
 namespace tsc = thor_scsi;
 
 unsigned long thor_scsi::core::binom(unsigned long n, unsigned long k) {
- unsigned long c = 1, i;
-
- if (k > n-k) // take advantage of symmetry
- k = n-k;
-
- for (i = 1; i <= k; i++, n--) {
- if (c/i >= ULONG_MAX/n){
- // raise an exception
- assert(0);
- return 0;
- }
- // split c * n / i into (c / i * i + c % i) * n / i
- c = c / i * n + c % i * n / i;
- }
-
- return c;
+ unsigned long c = 1, i;
+
+ if (k > n-k) // take advantage of symmetry
+ k = n-k;
+
+ for (i = 1; i <= k; i++, n--) {
+ if (c/i >= ULONG_MAX/n){
+ std::stringstream strm;
+ strm << "(c/i = "<< c/i << ") >= (ULONG_MAX/n " << ULONG_MAX/n << ")."
+ << " n = " << n << ", k = " << k << ", c = " << c
+ ;
+ throw std::overflow_error(strm.str());
+ return 0;
+ }
+ // split c * n / i into (c / i * i + c % i) * n / i
+ c = c / i * n + c % i * n / i;
+ }
+
+ return c;
 }
 
 #if 0
 double thor_scsi::binom (int const n, int const k)
 {
- int i;
- double num, den;
+ int i;
+ double num, den;
 
- if (k<0 || n<0 || k>n){
- return -1.0;
- }
+ if (k<0 || n<0 || k>n){
+ return -1.0;
+ }
 
- if (k==0 || k==n){
- return 1.0;
- }
+ if (k==0 || k==n){
+ return 1.0;
+ }
 
- if (k==1 || k==n-1){
- return n;
- }
+ if (k==1 || k==n-1){
+ return n;
+ }
 
- num = den = 1.0;
+ num = den = 1.0;
 
- for(i=n-k+1; i<=n; i++){
- num *= i;
- }
- for(i=1; i<=k; i++){
- den *= i;
- }
+ for(i=n-k+1; i<=n; i++){
+ num *= i;
+ }
+ for(i=1; i<=k; i++){
+ den *= i;
+ }
 
- return num/den;
+ return num/den;
 }
 #endif
 /*

src/thor_scsi/core/math_comb.h

@@ -1,12 +1,12 @@
 #ifndef _THOR_SCSI_MATH_COMB_H_
 #define _THOR_SCSI_MATH_COMB_H_
 namespace thor_scsi::core {
- /**
-  * @brief calculate binomial coefficient
-  *
- * @ todo: check if in standard library or in GSL
-  */
- unsigned long binom (unsigned long const n, unsigned long const k);
+ /**
+  * @brief calculate binomial coefficient
+  *
+ * @todo: check if in standard library or in GSL
+  */
+ unsigned long binom (unsigned long const n, unsigned long const k);
 }
 #endif /* THOR_SCSI_MATH_COMB_H */
 /*

src/thor_scsi/core/test_math_comb.cc

@@ -0,0 +1,57 @@
+#define BOOST_TEST_MODULE transform
+#define BOOST_TEST_DYN_LINK
+
+#include <boost/test/unit_test.hpp>
+#include <thor_scsi/core/math_comb.h>
+
+auto binom = thor_scsi::core::binom;
+
+BOOST_AUTO_TEST_CASE(test10_binom_pascal_triangle)
+{
+ //BOOST_TEST(binom(0, 1) == 0);
+ //BOOST_TEST(binom(0, 2) == 0);
+ //BOOST_TEST(binom(0, 3) == 0);
+
+ // the top
+ BOOST_TEST(binom(1, 1) == 1);
+
+ // quadratic expansion
+ BOOST_TEST(binom(2, 0) == 1);
+ BOOST_TEST(binom(2, 1) == 2);
+ BOOST_TEST(binom(2, 2) == 1);
+
+ // cubic expansion
+ BOOST_TEST(binom(3, 0) == 1);
+ BOOST_TEST(binom(3, 1) == 3);
+ BOOST_TEST(binom(3, 2) == 3);
+ BOOST_TEST(binom(3, 3) == 1);
+
+ // 4th order expansion
+ BOOST_TEST(binom(4, 0) == 1);
+ BOOST_TEST(binom(4, 1) == 4);
+ BOOST_TEST(binom(4, 2) == 6);
+ BOOST_TEST(binom(4, 3) == 4);
+ BOOST_TEST(binom(4, 4) == 1);
+
+ // 5th order expansion
+ BOOST_TEST(binom(5, 0) == 1);
+ BOOST_TEST(binom(5, 1) == 5);
+ BOOST_TEST(binom(5, 2) == 10);
+ BOOST_TEST(binom(5, 3) == 10);
+ BOOST_TEST(binom(5, 4) == 5);
+ BOOST_TEST(binom(5, 5) == 1);
+
+ // 6th order expansion
+ BOOST_TEST(binom(6, 0) == 1);
+ BOOST_TEST(binom(6, 1) == 6);
+ BOOST_TEST(binom(6, 2) == 15);
+ BOOST_TEST(binom(6, 3) == 20);
+ BOOST_TEST(binom(6, 4) == 15);
+ BOOST_TEST(binom(6, 5) == 6);
+ BOOST_TEST(binom(6, 6) == 1);
+}
+
+BOOST_AUTO_TEST_CASE(test10_binom_too_large)
+{
+ BOOST_CHECK_THROW(binom(2000 * 1000, 3000), std::overflow_error);
+}

src/thor_scsi/core/test_two_dimensional_multipoles.cc

@@ -3,7 +3,6 @@
 #include <boost/test/unit_test.hpp>
 #include <boost/test/tools/output_test_stream.hpp>
 #include <thor_scsi/core/multipoles.h>
-#include <gtpsa/utils_tps.hpp>
 #include <cmath>
 #include <array>
 
@@ -266,18 +265,6 @@ BOOST_AUTO_TEST_CASE(test20_set_harmonic_quadrupole)
  << "By\n " << By << "\n";
 
  }
- {
-
- auto x = tps(0, 0+1);
- auto y = tps(0, 2+1);
- tps Bx, By;
-
- std::cout << "x\n " << x << "\n"
- << "y\n " << y << "\n";
- h.field(x, y, &Bx, &By);
- std::cout << "Bx\n " << Bx << "\n"
- << "By\n " << By << "\n";
- }
  {
  auto x = t_ref.clone(), y = t_ref.clone(), Bx = t_ref.clone(), By = t_ref.clone();
  x = 1;