diff --git a/qpms/__init__.py b/qpms/__init__.py index e7caea7..1cfe334 100644 --- a/qpms/__init__.py +++ b/qpms/__init__.py @@ -7,7 +7,7 @@ from sys import platform as __platform import warnings as __warnings try: - from .qpms_c import PointGroup, FinitePointGroup, FinitePointGroupElement, Particle, scatsystem_set_nthreads, ScatteringSystem, ScatteringMatrix, pitau, set_gsl_pythonic_error_handling, pgsl_ignore_error, gamma_inc, lll_reduce + from .qpms_c import PointGroup, FinitePointGroup, FinitePointGroupElement, Particle, scatsystem_set_nthreads, ScatteringSystem, ScatteringMatrix, pitau, set_gsl_pythonic_error_handling, pgsl_ignore_error, gamma_inc, lll_reduce, ScatteringSystemCachingMode except ImportError as ex: if __platform == "linux" or __platform == "linux2": if 'LD_LIBRARY_PATH' not in __os.environ.keys(): diff --git a/qpms/qpms_c.pyx b/qpms/qpms_c.pyx index be3d6f5..162b77c 100644 --- a/qpms/qpms_c.pyx +++ b/qpms/qpms_c.pyx @@ -16,7 +16,7 @@ from .cytmatrices cimport CTMatrix, TMatrixFunction, TMatrixGenerator from .cymaterials cimport EpsMuGenerator from libc.stdlib cimport malloc, free, calloc import warnings - +import enum # Set custom GSL error handler. N.B. this is obviously not thread-safe. cdef char *pgsl_err_reason @@ -331,6 +331,10 @@ cpdef void scatsystem_set_nthreads(long n): qpms_scatsystem_set_nthreads(n) return +class ScatteringSystemCachingMode(enum.IntEnum): + NEVER = QPMS_SS_CACHE_NEVER + AUTO = QPMS_SS_CACHE_AUTO + ALWAYS = QPMS_SS_CACHE_ALWAYS cdef class ScatteringSystem: ''' @@ -357,7 +361,9 @@ cdef class ScatteringSystem: #TODO is there a way to disable the constructor outside this module? @staticmethod # We don't have any "standard" constructor for this right now - def create(particles, medium, FinitePointGroup sym, cdouble omega): # TODO tolerances + def create(particles, medium, FinitePointGroup sym, cdouble omega, + caching_mode = QPMS_SS_CACHE_DEFAULT): # TODO tolerances + caching_mode = ScatteringSystemCachingMode(caching_mode) # These we are going to construct cdef ScatteringSystem self cdef _ScatteringSystemAtOmega pyssw @@ -416,6 +422,7 @@ cdef class ScatteringSystem: orig.p[pi].pos = p.cval().pos orig.p[pi].tmatrix_id = tmindices[tm_derived_key] ssw = qpms_scatsys_apply_symmetry(&orig, sym.rawpointer(), omega, &QPMS_TOLERANCE_DEFAULT) + qpms_ss_create_translation_cache(ssw, caching_mode) ss = ssw[0].ss finally: free(orig.tmg) diff --git a/qpms/qpms_cdefs.pxd b/qpms/qpms_cdefs.pxd index a492820..9637655 100644 --- a/qpms/qpms_cdefs.pxd +++ b/qpms/qpms_cdefs.pxd @@ -550,10 +550,16 @@ cdef extern from "scatsystem.h": cdouble omega qpms_epsmu_t medium cdouble wavenumber + ctypedef enum qpms_ss_caching_mode_t: + QPMS_SS_CACHE_NEVER + QPMS_SS_CACHE_ALWAYS + QPMS_SS_CACHE_DEFAULT + QPMS_SS_CACHE_AUTO qpms_scatsys_at_omega_t *qpms_scatsys_apply_symmetry(const qpms_scatsys_t *orig, const qpms_finite_group_t *sym, cdouble omega, const qpms_tolerance_spec_t *tol) qpms_scatsys_at_omega_t *qpms_scatsys_at_omega(const qpms_scatsys_t *ss, cdouble omega) void qpms_scatsys_at_omega_free(qpms_scatsys_at_omega_t *ssw) + int qpms_ss_create_translation_cache(qpms_scatsys_t *ss, qpms_caching_mode_t m) cdouble *qpms_scatsys_irrep_pack_matrix(cdouble *target_packed, const cdouble *orig_full, const qpms_scatsys_t *ss, qpms_iri_t iri) cdouble *qpms_scatsys_irrep_unpack_matrix(cdouble *target_full, diff --git a/qpms/scatsys_translation_booster.c b/qpms/scatsys_translation_booster.c index fce9e7b..de490f1 100644 --- a/qpms/scatsys_translation_booster.c +++ b/qpms/scatsys_translation_booster.c @@ -128,6 +128,30 @@ booster_t *qpms_scatsys_translation_booster_create( return b; } +int qpms_ss_create_translation_cache(qpms_scatsys_t *ss, qpms_ss_caching_mode_t m) { +QPMS_ASSERT(ss); + if (ss->tbooster) { + QPMS_WARN("Translation cache already created?"); + return 0; + } + switch(m) { + case QPMS_SS_CACHE_NEVER: + return 0; + case QPMS_SS_CACHE_AUTO: + QPMS_WARN("Translation operator cache heuristics not implemented, creating the cache"); + case QPMS_SS_CACHE_ALWAYS: + ss->tbooster = qpms_scatsys_translation_booster_create(ss); + if (ss->tbooster) return 0; + else { + QPMS_WARN("Failed to create tranlation operator cache"); + return -1; + } + default: + QPMS_WTF; + } + QPMS_WTF; +} + static qpms_errno_t qpms_scatsys_translation_booster_eval_bessels( const booster_t *b, complex double *target, complex double k // includes ref. ind. effect ) { diff --git a/qpms/scatsystem.c b/qpms/scatsystem.c index 1495051..01ca3cc 100644 --- a/qpms/scatsystem.c +++ b/qpms/scatsystem.c @@ -473,6 +473,8 @@ void qpms_scatsys_free(qpms_scatsys_t *ss) { free(ss->orbit_types); free(ss->saecv_sizes); free(ss->p_by_orbit); + if(ss->tbooster) + qpms_scatsys_translation_booster_free(ss->tbooster); qpms_trans_calculator_free(ss->c); } free(ss); diff --git a/qpms/scatsystem.h b/qpms/scatsystem.h index cb48c96..5826175 100644 --- a/qpms/scatsystem.h +++ b/qpms/scatsystem.h @@ -220,6 +220,7 @@ typedef struct qpms_scatsys_at_omega_t { struct qpms_scatsysw_translation_booster *translation_cache; ///< (private) cache to speedup tranlations } qpms_scatsys_at_omega_t; + /// Creates a new scatsys by applying a symmetry group onto a "proto-scatsys", copying particles if needed. /** In fact, it copies everything except the vswf set specs and qpms_abstract_tmatrix_t instances, * so keep them alive until scatsys is destroyed. @@ -269,6 +270,25 @@ void qpms_scatsys_at_omega_free(qpms_scatsys_at_omega_t *ssw); qpms_scatsys_at_omega_t *qpms_scatsys_at_omega(const qpms_scatsys_t *ss, complex double omega); +/// Determines behaviour of qpms_ss_create_translation_cache(). +typedef enum qpms_ss_caching_mode { + /// Don't create the translation operator cache. + /** Use this if the particle positions are random. */ + QPMS_SS_CACHE_NEVER, + /// Always create the translation operator cache. + /** Use this for highly regular large finite structures. */ + QPMS_SS_CACHE_ALWAYS, + /// Evaluate the need of caching automatically. + QPMS_SS_CACHE_AUTO, + QPMS_SS_CACHE_DEFAULT = QPMS_SS_CACHE_AUTO +} qpms_ss_caching_mode_t; + +/// Creates some data structures for speeding up translation operator calculations. +/** This is mostly useful for "finite lattices", where many pairs of nanoparticles + * share the same relative positions. + */ +int qpms_ss_create_translation_cache(qpms_scatsys_t *ss, qpms_ss_caching_mode_t m); + /// Creates a "full" transformation matrix U that takes a full vector and projects it onto an symmetry adapted basis. /** Mostly as a reference and a debugging tool, as multiplicating these big matrices would be inefficient. * @@ -548,7 +568,9 @@ struct beyn_result_t *qpms_scatsys_finite_find_eigenmodes( double res_tol ///< (default: `0.0`) TODO DOC. ); +/// "private" destructor, called by qpms_scatsys_free() void qpms_scatsys_translation_booster_free(struct qpms_scatsys_translation_booster *); +/// "private" constructor, use qpms_ss_create_translation_cache() instead. struct qpms_scatsys_translation_booster *qpms_scatsys_translation_booster_create( const qpms_scatsys_t *ss);