diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index 591a0187..609c003e 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -2,14 +2,14 @@ exclude: 'doc/conf.py' repos: - repo: https://github.com/asottile/pyupgrade - rev: v2.20.0 + rev: v2.37.3 hooks: - id: pyupgrade # for now don't force to change from %-operator to {} args: [--keep-percent-format, --py37-plus] - repo: https://github.com/pre-commit/pre-commit-hooks - rev: v4.0.1 + rev: v4.3.0 hooks: - id: check-ast - id: check-builtin-literals @@ -22,18 +22,18 @@ repos: args: [--remove] - repo: https://github.com/pre-commit/mirrors-autopep8 - rev: 'v1.5.7' + rev: 'v1.7.0' hooks: - id: autopep8 args: [--max-line-length=120] - repo: https://github.com/PyCQA/isort/ - rev: 5.9.1 + rev: 5.10.1 hooks: - id: isort - repo: https://github.com/nbQA-dev/nbQA - rev: 0.13.1 + rev: 1.4.0 hooks: - id: nbqa-black additional_dependencies: [black==20.8b1] diff --git a/.travis.yml b/.travis.yml index ce1af5ca..c8f774dd 100644 --- a/.travis.yml +++ b/.travis.yml @@ -1,11 +1,14 @@ # C++ build language: cpp -compiler: - - gcc - - clang -os: -- osx -- linux + +jobs: + include: + - os: linux + compiler: gcc + - os: linux + compiler: clang + - os: osx + compiler: gcc #before_install: # - sudo apt-get install cmake @@ -35,5 +38,3 @@ script: - ./oaconference -N 12 -k 6 -v 0 - ./oapareto result-20.2-2-2.oa -f B -o dummy.oa - - diff --git a/CMakeLists.txt b/CMakeLists.txt index 5c676ef6..ea394e98 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1,5 +1,5 @@ -CMAKE_MINIMUM_REQUIRED(VERSION 2.8) +CMAKE_MINIMUM_REQUIRED(VERSION 3.0) PROJECT(OApackage) @@ -8,9 +8,7 @@ PROJECT(OApackage) # ---------------------------------------------------------------------------- if(WIN32) - #message(STATUS "Set build type...") set(CMAKE_CONFIGURATION_TYPES Release CACHE TYPE INTERNAL FORCE) # no Debug - #set(CMAKE_BUILD_TYPE Release STRING) endif() if (NOT WIN32) @@ -27,8 +25,6 @@ if(ZLIB_FOUND) endif() message(STATUS "ZLIB_FOUND: ${ZLIB_FOUND}") if(ZLIB_FOUND) - #message(STATUS "ZLIB_INCLUDE_DIRS: ${ZLIB_INCLUDE_DIRS}") - #message(STATUS "ZLIB_LIBRARIES: ${ZLIB_LIBRARIES}") include_directories(${ZLIB_INCLUDE_DIRS}) endif() @@ -38,12 +34,13 @@ if (NumPy_FOUND) message(STATUS "NumPy was found") endif() -# find_package(MPI) # Needed for legacy code and OpenMP functionality # ---------------------------------------------------------------------------- # Code # ---------------------------------------------------------------------------- +set (CMAKE_CXX_STANDARD 11) + include_directories(.) include_directories(src) include_directories(nauty) @@ -151,10 +148,6 @@ if("$ENV{OADEBUG}" STREQUAL "1") set(COPTS "-g -O0 -DOADEBUG" ) # -std=c++0x -Weffc++") endif() -if(WIN32) -else() - set(LINKOPTS "${LINKOPTS}") # or use -s option -endif() if(WIN32) set(COPTS "/MT /wd4018 /wd4996") @@ -167,7 +160,7 @@ endif() set(COPTS "${COPTS} -DFULLPACKAGE") -if(ZLIB_FOUND AND 1) +if(ZLIB_FOUND) set(COPTS "${COPTS} -DUSEZLIB") include_directories(ZLIB_INCLUDE_DIRS) else() @@ -208,7 +201,7 @@ set(progs oacat oajoin oapareto oasplit oaanalyse oainfo oaunittest oaconvert oa set(progsextend ) list(APPEND progs oacheck oastreaming oaranktest oatest oaclustergather) -list(APPEND progs oa_depth_extend) +list(APPEND progs oa_depth_extend oa_select_maxj) if (WIN32) add_library(oalib STATIC ${srcsextend} ${headersextend} ${nautysrc} ${graphsrc}) @@ -240,39 +233,6 @@ endif() TARGET_LINK_LIBRARIES(oaextendsingle oalib) set(extendprogs oaextendsingle ) - if (MPI_FOUND AND NOT WIN32 AND 1) - message(STATUS "Found MPI package: adding oaextendmpi") - if (VERBOSE) - message(STATUS "MPI executable: ${MPIEXEC}, libs ${MPI_LIBRARIES}") - endif() - #message(STATUS " libs ${MPI_LIBRARIES}; flags ${MPI_COMPILE_FLAGS}; libs ${MPI_LIBRARIES}") - include_directories( ${MPI_INCLUDE_PATH} ) - include_directories( src/mpitools ) - - add_executable(oaextendmpi EXCLUDE_FROM_ALL "utils/oaextend.cpp" ${headersextend} ${srcsextend} "src/mpitools/mpitools.cpp") - target_link_libraries(oaextendmpi ${MPI_LIBRARIES}) - if(MPI_COMPILE_FLAGS) - set_target_properties(oaextendmpi PROPERTIES COMPILE_FLAGS "${MPI_COMPILE_FLAGS}") - endif() - if(MPI_LINK_FLAGS) - if (WIN32) - set_target_properties(oaextendmpi PROPERTIES LINK_FLAGS "${MPI_LINK_FLAGS} /MT") - endif() - set_target_properties(oaextendmpi PROPERTIES LINK_FLAGS "${MPI_LINK_FLAGS}") - endif() - - set_target_properties(oaextendmpi PROPERTIES COMPILE_DEFINITIONS "OAEXTEND_MULTICORE;NEWINTERFACE" ) - - if (WIN32) - if(ZLIB_FOUND) - TARGET_LINK_LIBRARIES(oaextendmpi ${ZLIB_LIBRARIES}) - endif() - else() - TARGET_LINK_LIBRARIES(oaextendmpi m ${ZLIB_LIBRARIES}) - endif() - - endif(MPI_FOUND AND NOT WIN32 AND 1) - message(STATUS "Extend progs: ${extendprogs}") message(STATUS "Progs: ${progs}") diff --git a/CONTRIBUTORS.md b/CONTRIBUTORS.md index 655a14a0..29d37bdf 100644 --- a/CONTRIBUTORS.md +++ b/CONTRIBUTORS.md @@ -3,5 +3,4 @@ * Pieter Eendebak * Eric Schoen * Alan Vazquez-Alcocer - - +* Alexandre Bohyn diff --git a/appveyor.yml b/appveyor.yml index 6eb8b96a..dcf21a78 100644 --- a/appveyor.yml +++ b/appveyor.yml @@ -11,6 +11,10 @@ environment: matrix: + - PYTHON: "C:\\Python310-x64" + PYTHON_VERSION: "3.10.4" + PYTHON_ARCH: "64" + - PYTHON: "C:\\Python39-x64" PYTHON_VERSION: "3.9.1" PYTHON_ARCH: "64" @@ -27,10 +31,6 @@ environment: PYTHON_VERSION: "3.7.x" PYTHON_ARCH: "64" - - PYTHON: "C:\\Python37" - PYTHON_VERSION: "3.7.x" - PYTHON_ARCH: "32" - matrix: fast_finish: true @@ -100,7 +100,7 @@ install: - "%CMD_IN_ENV% pip install mock" - "%CMD_IN_ENV% python -c \"import mock\"" - - cinst -y swig --version 4.0.1 + - cinst -y swig # --version 3.0.12 build: off diff --git a/docs/examples/example_mixed_level_gwlp.ipynb b/docs/examples/example_mixed_level_gwlp.ipynb index b33da2d9..8a46a7a9 100644 --- a/docs/examples/example_mixed_level_gwlp.ipynb +++ b/docs/examples/example_mixed_level_gwlp.ipynb @@ -11,13 +11,43 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "References [xx]" + "This notebook contains example code to compute the type-specific generalized word-length pattern (GWLP) for mixed-level designs, more specifically for regular four-and-two-level designs.\n", + "\n", + "In four-and-two-level designs, the four-level factors are constructed us the grouping scheme of [Wu & Zhang (1989)](https://doi.org/10.1093/biomet/80.1.203) where the levels of a four-level factor $A$ are based on the levels of three two-level factors $a_1$, $a_2$ and $a_3$, called the pseudo-factors, where $I=a_1a_2a_3$:\n", + "$$\n", + "\\begin{array}{ccccc}\n", + " a_1 & a_2 & a_3 & & A \\\\\n", + " 1 & 1 & 0 & \\rightarrow & 0 \\\\\n", + " 1 & 0 & 1 & \\rightarrow & 1 \\\\\n", + " 0 & 1 & 1 & \\rightarrow & 2 \\\\\n", + " 0 & 0 & 0 & \\rightarrow & 3 \\\\\n", + " \\end{array}\n", + "$$\n", + "The two-level factors can either be main factors or be entirely aliased with a combination of the main factors. Such aliased factors are called added factors. If one of the main factors, used in an added factor, is also used as pseudo-factor in a four-level factor, then the added factor has type $I$. If it is used as pseudo-factor in two distinct four-level factor, it has type $II$, etc $\\ldots$." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "\n", + "import oapackage" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Create a $4^{1}2^{7}$ four-and-two-level regular design in 32 runs with 1 four-level factor and 7 two-level factors." + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, "outputs": [ { "name": "stdout", @@ -60,16 +90,22 @@ } ], "source": [ - "import numpy as np\n", - "\n", - "import oapackage\n", - "\n", "array = oapackage.exampleArray(56, 1)\n", "array = array.selectFirstColumns(7)\n", "arrayclass = oapackage.arraylink2arraydata(array)\n", "array.showarraycompact()" ] }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The GWLP for mixed-level designs is computed using the adapted forms of equations (7) and (8) in [Xu and Wu (2001)](https://www.doi.org/10.1214/aos/1009210552).\n", + "\n", + "## Distance distribution\n", + "First, the distance distribution is first computed for all combinations of $(i,j)$ with $i \\in \\{0,1\\}$ and $j=0,\\ldots,7$ and $i\\neq j$. In this example with a $4^12^{7}$ design, that is a $2 \\times 8$ matrix $B$, where $B_{i,j}$ is the number of rows that have $(1-i)$ different elements if their four-level parts and $(7-j)$ different elements in their two-level parts, divided by the number of runs." + ] + }, { "cell_type": "code", "execution_count": 3, @@ -79,7 +115,8 @@ "name": "stdout", "output_type": "stream", "text": [ - "distance distribution for mixed-level design [[1 0 1 4 1 0 1]\n", + "distance distribution for mixed-level design\n", + "[[1 0 1 4 1 0 1]\n", " [0 2 6 8 6 2 0]]\n" ] } @@ -87,7 +124,22 @@ "source": [ "Dm = oapackage.distance_distribution_mixed(array, 0)\n", "D = np.array(Dm).astype(int)\n", - "print(f\"distance distribution for mixed-level design {D}\")" + "print(f\"distance distribution for mixed-level design\\n{D}\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "This is verified since the sum of the matrix is 32 which is equal to $\\binom{32}{2} 32^{-1}$.\n", + "\n", + "## McWilliams Transforms\n", + "\n", + "Now, the second step in the computation of the GWLP, is to use the McWilliams Transforms to obtain the GWLP from the distance distribution matrix. The McWilliams Transform will create another $2 \\times 8$ matrix $B^{\\prime}$, obtained using the following formula\n", + "$$\n", + "B_{j_{1}, j_{2}}^{\\prime}=N^{-1} \\sum_{i_{1}=0}^{m} \\sum_{i_{2}=0}^{n} B_{i_{1}, i_{2}} P_{j_{1}}\\left(i_{1} ; 1, 4\\right) P_{j_{2}}\\left(i_{2} ; 7, 2\\right)\n", + "$$\n", + "where $P_j(x,n,s)$ is the Krawtchouck polynomials for a total of $n$ factors with $s$ levels." ] }, { @@ -108,17 +160,24 @@ "source": [ "N = array.n_rows\n", "factor_levels_for_groups = arrayclass.factor_levels_column_groups()\n", - "\n", - "Bprime = oapackage.macwilliams_transform_mixed(\n", - " Dm, N, factor_levels_for_groups=arrayclass.factor_levels_column_groups(), verbose=0\n", - ")\n", + "Bprime = oapackage.macwilliams_transform_mixed(Dm, N, factor_levels_for_groups, verbose=0)\n", "print(f\"MacWilliams transform:\")\n", - "print(np.array(Bprime))" + "print(np.array(Bprime).astype(int))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "This matrix is equivalent to the type specific generalized word-length pattern, where the row index indicates the type of words (0 to 1) and the column index indicates the length of the words (0 to 7).\n", + "\n", + "The generalized word-length pattern ($A$) can be obtained by summing the rows of the $B^{\\prime}$ matrix anti-diagonally. That is:\n", + "$$A_j(D) = \\sum_{i^{\\prime}+j^{\\prime}=j} B^{\\prime}_{i^{\\prime},j^{\\prime}}$$" ] }, { "cell_type": "code", - "execution_count": 5, + "execution_count": 4, "metadata": {}, "outputs": [ { @@ -136,6 +195,9 @@ } ], "metadata": { + "interpreter": { + "hash": "3a893c4079a4533c9db639b850e04910700d284bba4afe56b82dceecf522f906" + }, "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", diff --git a/oalib.i b/oalib.i index d974c3a0..1e86a70b 100644 --- a/oalib.i +++ b/oalib.i @@ -516,6 +516,13 @@ public: } } +%extend jstructconference_t { +public: + std::string __repr__() { + return $self->showstr(); + } +} + #ifdef SWIGPYTHON // Add module docstring %pythoncode diff --git a/oapackage/oahelper.py b/oapackage/oahelper.py index 0cdc8428..225bf26e 100644 --- a/oapackage/oahelper.py +++ b/oapackage/oahelper.py @@ -578,6 +578,15 @@ def write_text_arrayfile(filename: str, designs: List[Any], comment: str = None) afile.append_arrays(designs) afile.closefile() +def arrayfile_generator(afile): + """ Return generator to read all files in the array file """ + af = oapackage.arrayfile_t(afile) + + na=af.narrays + for ii in range(na): + yield af.readnext() + af.closefile() + def runcommand(cmd: str, dryrun=0, idstr: Optional[None] = None, verbose: int = 1, logfile: Optional[str] = None, shell: bool = True): """ Run specified command in external environment diff --git a/setup.py b/setup.py index 0d1e6490..35a24821 100644 --- a/setup.py +++ b/setup.py @@ -14,7 +14,6 @@ from os import path import setuptools.command.build_ext -# %% Load packages from setuptools import Extension, find_packages, setup from setuptools.command.install import install as setuptools_install from setuptools.command.test import test as TestCommand @@ -384,6 +383,8 @@ def readme(): 'Programming Language :: Python :: 3.7', "Programming Language :: Python :: 3.8", "Programming Language :: Python :: 3.9", + "Programming Language :: Python :: 3.10", + "Programming Language :: Python :: 3.11", 'License :: OSI Approved :: BSD License' ] ) diff --git a/src/arrayproperties.cpp b/src/arrayproperties.cpp index 67d24013..f383889b 100644 --- a/src/arrayproperties.cpp +++ b/src/arrayproperties.cpp @@ -40,7 +40,7 @@ mvalue_t< long > A3A4 (const array_link &al) { w.push_back (w3); w.push_back (w4); - mvalue_t< long > wm (w, mvalue_t< long >::LOW); + mvalue_t< long > wm (w, mvalue_t< long >::direction_t::LOW); return wm; } @@ -520,16 +520,16 @@ std::vector< double > GWLP (const array_link &al, int verbose, int truncate) { int N = al.n_rows; int n = al.n_columns; - int domixed = al.is_mixed_level(); + int is2level = al.is2level(); if (verbose) - myprintf ("GWLP: N %d, domixed %d\n", N, domixed); + myprintf ("GWLP: N %d, is2level %d\n", N, is2level); - if (domixed) { - std::vector< double > gma = GWLPmixed (al, verbose, truncate); + if (is2level) { + std::vector< double > gma = GWLP_two_level_design(al, verbose, truncate); return gma; } else { - std::vector< double > gma = GWLP_two_level_design(al, verbose, truncate); + std::vector< double > gma = GWLPmixed (al, verbose, truncate); return gma; } } diff --git a/src/arrayproperties.h b/src/arrayproperties.h index ca1818cd..a7eb009f 100644 --- a/src/arrayproperties.h +++ b/src/arrayproperties.h @@ -604,7 +604,7 @@ inline mvalue_t< long > F4 (const array_link &al, int verbose = 1) { myprintf ("\n"); } - mvalue_t< long > v (FF, mvalue_t< long >::LOW); + mvalue_t< long > v (FF, mvalue_t< long >::direction_t::LOW); return v; } diff --git a/src/arraytools.cpp b/src/arraytools.cpp index e06cf17c..5fb2e514 100644 --- a/src/arraytools.cpp +++ b/src/arraytools.cpp @@ -633,19 +633,8 @@ std::vector< int > getJcounts (arraylist_t *arraylist, int N, int k, int verbose int jl = 5; while (jl <= al.n_columns) { jstruct_t js (al.n_rows, al.n_columns, jl); - if (0) { - arraylist_t *all = new arraylist_t (); - all->push_back (al); - std::vector< jstruct_t > xx = analyseArrays (*all, verbose, jl); - js = jstruct_t (xx[0]); - if (verbose >= 2) { - myprintf (" old: "); - xx[0].show (); - } - delete all; - } else { - foldtest (js, al, jl, verbose); - } + + foldtest (js, al, jl, verbose); if (!js.allzero ()) { if (verbose >= 3) { myprintf (" new: "); @@ -1093,6 +1082,7 @@ int compareLMC(const array_link &lhs, const array_link &rhs) { * \return Selected example array */ array_link exampleArray (int idx, int verbose) { + /* Note: to generate C code from python, use oapackage.oahelper.formatC(A) */ if (idx == -1) { for (int i = 0; i < 1000; i++) { try { array_link al = exampleArray(i, verbose); } @@ -2072,7 +2062,15 @@ array_link exampleArray (int idx, int verbose) { array.setarraydata(array_data_tmp, array.n_rows* array.n_columns); return array; } - + case 57: { + dstr = "design in OA(18, 3^{6})"; + if (verbose) { + myprintf("exampleArray %d: %s\n", idx, dstr.c_str()); + } + array_link array ( 18,6, 0 ); + int array_data_tmp[] = {0,0,0,0,0,0,1,1,1,1,1,1,2,2,2,2,2,2,0,0,1,1,2,2,0,0,1,1,2,2,0,0,1,1,2,2,0,0,1,1,2,2,1,2,0,2,0,1,1,2,0,2,0,1,0,1,0,2,1,2,1,0,2,0,2,1,2,2,1,1,0,0,0,1,0,2,1,2,2,2,1,1,0,0,1,0,2,0,2,1,0,1,2,1,2,0,0,2,0,1,2,1,2,1,2,0,1,0}; array.setarraydata(array_data_tmp, array.n_rows * array.n_columns); + return array; + } } // end of switch return array_link (); @@ -2217,10 +2215,7 @@ bool array_link::is2level () const { int N = this->n_rows; for (int r = 0; r < this->n_rows; r++) { for (int c = 0; c < this->n_columns; c++) { - if (this->array[r + c * N] < 0) { - return false; - } - if (this->array[r + c * N] > 1) { + if (this->array[r + c * N] < 0 || this->array[r + c * N] > 1) { return false; } } @@ -2852,6 +2847,8 @@ std::vector< int > array_link::Jcharacteristics (int jj) const { } } +/** Calculate J_k for the specified k-tuple of columns +*/ int jvalue_conference (const array_link &ar, const int J, const int *column_indices) { int jval = 0; @@ -3299,11 +3296,9 @@ array_link arraydata_t::create_root (int n_columns, int fill_value) const { array_link al (this->N, n_columns, -1); al.create_root (*this); - for (int i = this->strength; i < al.n_columns; i++) { - for (int r = 0; r < this->N; r++) { - al.at (r, i) = fill_value; - } - } + const int N = this->N; + std::fill(al.array+this->strength*N, al.array+ al.n_columns*N, fill_value); + return al; } @@ -3656,6 +3651,7 @@ void jstruct_t::calcj4 (const array_link &al) { void jstruct_t::calcj5 (const array_link &al) { myassert (jj == 5, "jj should be 5"); + myassert( k==al.n_columns, "number of columns of input array should match jstruct_t"); int *pp = new_perm_init< int > (jj); int ncolcombs = ncombs (k, jj); const int nr = al.n_rows; @@ -3790,8 +3786,9 @@ void jstruct_t::showdata () { } std::string jstructbase_t::showstr () { - std::string s = "jstruct_t: " + printfstring ("jj %d, values ", jj); - return s; + std::string values_string = permutation2string(this->values, 16); + std::string s = "jstruct_t: " + printfstring ("jj %d, values %s", this->jj, values_string.c_str()); + return s; } void jstructbase_t::show () { std::cout << "jstruct_t: " << printfstring ("jj %d, values ", jj); @@ -4604,7 +4601,7 @@ void arrayfile_t::writeheader () { if (this->mode == arrayfile::ABINARY_DIFFZERO) { // NOTE: needed because in ABINARY_DIFFZERO we diff modulo 2 if (this->nbits != 1) - myprintf ("not implemented...!\n"); + myprintf ("arrayfile_t::writeheader: not implemented for nbits %d!\n", this->nbits); assert (this->nbits == 1); } fwrite ((const void *)&magic, sizeof (int), 1, this->nfid); @@ -4654,18 +4651,6 @@ void update_array_link(array_link &al, long* pymatinput, int number_of_rows, int #endif -std::string printfstring(const char *message, ...) { - char buf[32 * 1024]; - - va_list va; - va_start(va, message); - vsprintf(buf, message, va); - va_end(va); - - std::string str(buf); - return str; -} - /** * @brief Read file with design of OA * @param file @@ -5364,8 +5349,7 @@ void arrayfile_t::write_array_binary_diffzero (const array_link &A) { this->write_array_binary (z); } else { - array_link z = rest; - this->write_array_binary (z); + this->write_array_binary (rest); } // update with previous array @@ -5414,7 +5398,7 @@ void arrayfile_t::write_array_binary (carray_t *array, const int nrows, const in if (array[i]) { bit_array_set_bit (bitarr, i); } else { - bit_array_clear_bit (bitarr, i); + bit_array_clear_bit_fast (bitarr, i); } } word_addr_t num_of_words = diff --git a/src/arraytools.h b/src/arraytools.h index ead3eeb0..3985b77f 100644 --- a/src/arraytools.h +++ b/src/arraytools.h @@ -116,7 +116,6 @@ typedef const short int carray_t; /* change definition below together with array_t !!!! */ #define MPI_ARRAY_T MPI_SHORT -/*other options for MPI_ARRAY_T are: char: MPI_CHAR, short: MPI_SHORT, int: MPI_INT, long: MPI_LONG */ typedef short int rowindex_t; /** type used for row indexing */ typedef int colindex_t; /** type used for column indexing */ @@ -340,12 +339,6 @@ struct arraydata_t { /// Read array configuration from file arraydata_t *readConfigFile (const char *file); -/** - * @brief Function similar to printf returning C++ style string - * @param message - * @return - */ -std::string printfstring(const char *message, ...); /** * @brief Make a copy of an array @@ -1017,6 +1010,9 @@ class jstruct_t { /// return 1 if all J values are zero, otherwise return 0 int allzero() const; + /// calculate J-characteristics of a 2-level array, special function for jj=5 + void calcj5(const array_link &al); + private: /// init data structures void init(int N, int k, int jj); @@ -1024,8 +1020,6 @@ class jstruct_t { void calc(const array_link &al); /// calculate J-characteristics of a 2-level array, special function for jj=4 void calcj4(const array_link &al); - /// calculate J-characteristics of a 2-level array, special function for jj=5 - void calcj5(const array_link &al); }; diff --git a/src/bit_array.cpp b/src/bit_array.cpp deleted file mode 100644 index ab06e898..00000000 --- a/src/bit_array.cpp +++ /dev/null @@ -1,559 +0,0 @@ -/* - bit_array.c - project: bit array C library - url: https://github.com/noporpoise/BitArray/ - Adapted from: http://stackoverflow.com/a/2633584/431087 - author: Isaac Turner - - Copyright (c) 2011, Isaac Turner - All rights reserved. - - Redistribution and use in source and binary forms, with or without - modification, are permitted provided that the following conditions are met: - * Redistributions of source code must retain the above copyright - notice, this list of conditions and the following disclaimer. - - THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND - ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED - WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE - DISCLAIMED. IN NO EVENT SHALL BE LIABLE FOR ANY - DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES - (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; - LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND - ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT - (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS - SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. -*/ - -#include -#include // ULONG_MAX -#include -#include -#include // memset - -#ifdef _WIN32 -#define inline __inline -#endif - -#include "bit_array.h" - -#define MIN(a, b) (((a) <= (b)) ? (a) : (b)) -#define MAX(a, b) (((a) >= (b)) ? (a) : (b)) - -// -// For internal use -// - -// sizeof gives size in bytes (8 bits per byte) -const int WORD_SIZE = sizeof (word_t) * 8; - -// Number of words required to store so many bits -word_addr_t nwords (bit_index_t b) { return (b + WORD_SIZE - 1) / WORD_SIZE; } - -// Index of word -inline word_addr_t bindex (bit_index_t b) { return b / WORD_SIZE; } - -// Offset within a word (values up to 64 most likely) -inline unsigned int boffset (bit_index_t b) { return b % WORD_SIZE; } - -//#define BIT_MASK(length) (((word_t)1 << (length))-1) // overflows -#define BIT_MASK(length) ((word_t)ULONG_MAX >> (WORD_SIZE - (length))) - -// -// Constructor -// -BIT_ARRAY *bit_array_create (bit_index_t nbits) { - BIT_ARRAY *bitarr = (BIT_ARRAY *)malloc (sizeof (BIT_ARRAY)); - - if (bitarr == NULL) { - // error - could not allocate enough memory - errno = ENOMEM; - return NULL; - } - - word_addr_t num_of_words = nwords (nbits); - -#ifdef DEBUG - printf ("Creating BIT_ARRAY (bits: %lu; words: %lu; WORD_SIZE: %i)\n", (unsigned long)nbits, - (unsigned long)num_of_words, WORD_SIZE); -#endif - - // bitarr->num_of_words = num_of_words; - bitarr->num_of_bits = nbits; - bitarr->words = (word_t *)malloc (sizeof (word_t) * num_of_words); - - if (bitarr->words == NULL) { - // error - could not allocate enough memory - free (bitarr); - errno = ENOMEM; - return NULL; - } - - // Initialise to zero - bit_array_fill_zeros (bitarr); - - return bitarr; -} - -// -// Destructor -// -void bit_array_free (BIT_ARRAY *bitarr) { - free (bitarr->words); - free (bitarr); -} - -// -// Methods -// -void bit_array_set_bit (BIT_ARRAY *bitarr, bit_index_t b) { - if (b < 0 || b >= bitarr->num_of_bits) { - // out of bounds error - fprintf (stderr, "bit_array.c: bit_array_set_bit() - " - "out of bounds error (index: %lu; length: %lu)\n", - b, bitarr->num_of_bits); - - errno = EDOM; - - throw; - exit (EXIT_FAILURE); - } - - bitarr->words[bindex (b)] |= ((word_t)1 << (boffset (b))); -} - -void bit_array_clear_bit (BIT_ARRAY *bitarr, bit_index_t b) { - if (b < 0 || b >= bitarr->num_of_bits) { - // out of bounds error - fprintf (stderr, "bit_array.c: bit_array_clear_bit() - " - "out of bounds error (index: %lu; length: %lu)\n", - b, bitarr->num_of_bits); - - errno = EDOM; - throw; - exit (EXIT_FAILURE); - } - - bitarr->words[bindex (b)] &= ~((word_t)1 << (boffset (b))); -} - -char bit_array_get_bit_nocheck (BIT_ARRAY *bitarr, bit_index_t b) { - return (bitarr->words[bindex (b)] >> (boffset (b))) & 0x1; -} - -char bit_array_get_bit (BIT_ARRAY *bitarr, bit_index_t b) { - if (b < 0 || b >= bitarr->num_of_bits) { - // out of bounds error - fprintf (stderr, "bit_array.c: bit_array_get_bit() - " - "out of bounds error (index: %lu; length: %lu)\n", - b, bitarr->num_of_bits); - - errno = EDOM; - throw; - exit (EXIT_FAILURE); - } - - return (bitarr->words[bindex (b)] >> (boffset (b))) & 0x1; -} - -/* set all elements of data to zero */ -void bit_array_fill_zeros (BIT_ARRAY *bitarr) { - // size_t num_of_bytes = (bitarr->num_of_bits / 8) + 1; - size_t num_of_bytes = nwords (bitarr->num_of_bits) * sizeof (word_t); - memset (bitarr->words, 0, num_of_bytes); -} - -/* set all elements of data to one */ -void bit_array_fill_ones (BIT_ARRAY *bitarr) { - size_t num_of_bytes = nwords (bitarr->num_of_bits) * sizeof (word_t); - memset (bitarr->words, 0xFF, num_of_bytes); -} - -// To string method (remember to free the result!) -char *bit_array_to_string (BIT_ARRAY *bitarr) { - char *str = (char *)malloc (sizeof (char) * (bitarr->num_of_bits + 1)); - - bit_index_t i; - - for (i = 0; i < bitarr->num_of_bits; i++) { - str[i] = bit_array_get_bit (bitarr, i) ? '1' : '0'; - } - - str[bitarr->num_of_bits] = '\0'; - - return str; -} - -BIT_ARRAY *bit_array_clone (BIT_ARRAY *bitarr) { - BIT_ARRAY *cpy = (BIT_ARRAY *)malloc (sizeof (BIT_ARRAY)); - - word_addr_t num_of_words = nwords (bitarr->num_of_bits); - - cpy->num_of_bits = bitarr->num_of_bits; - cpy->words = (word_t *)malloc (sizeof (word_t) * num_of_words); - - // Copy across bits - memcpy (cpy->words, bitarr->words, num_of_words * sizeof (word_t)); - - return cpy; -} - -/* -void bit_array_copy(BIT_ARRAY* dest, bit_index_t dstindx, - BIT_ARRAY* src, bit_index_t srcindx, bit_index_t length) -{ - -} -*/ - -// Enlarge or shrink the size of a bit array -// Shrinking will free some memory if it is large -// Enlarging an array will add zeros to the end of it -// returns 1 on success, 0 on failure -char bit_array_resize (BIT_ARRAY *bitarr, bit_index_t new_num_of_bits) { - bit_index_t old_num_of_bits = bitarr->num_of_bits; - bitarr->num_of_bits = new_num_of_bits; - - word_addr_t old_num_of_words = nwords (old_num_of_bits); - word_addr_t new_num_of_words = nwords (new_num_of_bits); - - if (new_num_of_words != old_num_of_words) { - // Need to change the amount of memory used - bitarr->words = (word_t *)realloc ((void *)bitarr->words, new_num_of_words * sizeof (word_t)); - - if (bitarr->words == NULL) { - // error - could not allocate enough memory - errno = ENOMEM; - return 0; - } - } - - // if we are growing - need to zero new bits - if (new_num_of_bits > old_num_of_bits) { - // zero entire words first - if (new_num_of_words > old_num_of_words) { - memset (bitarr->words + old_num_of_words, 0, - (new_num_of_words - old_num_of_words) * sizeof (word_t)); - } - - // zero bits on the end of what used to be the last word - unsigned int bits_on_last_word = boffset (old_num_of_bits); - - if (bits_on_last_word > 0) { - bitarr->words[old_num_of_words - 1] &= BIT_MASK (bits_on_last_word); - } - } - - return 1; -} - -// -// Logic operators -// - -void bit_array_and (BIT_ARRAY *dest, BIT_ARRAY *src1, BIT_ARRAY *src2) { - if (dest->num_of_bits != src1->num_of_bits || src1->num_of_bits != src2->num_of_bits) { - // error - fprintf (stderr, "bit_array.c: bit_array_and() : " - "dest, src1 and src2 must be of the same length\n"); - throw; - exit (EXIT_FAILURE); - } - - word_addr_t num_of_words = nwords (src1->num_of_bits); - - word_addr_t i; - for (i = 0; i < num_of_words; i++) { - dest->words[i] = src1->words[i] & src2->words[i]; - } -} - -void bit_array_or (BIT_ARRAY *dest, BIT_ARRAY *src1, BIT_ARRAY *src2) { - if (dest->num_of_bits != src1->num_of_bits || src1->num_of_bits != src2->num_of_bits) { - // error - fprintf (stderr, "bit_array.c: bit_array_and() : " - "dest, src1 and src2 must be of the same length\n"); - throw; - exit (EXIT_FAILURE); - } - - word_addr_t num_of_words = nwords (src1->num_of_bits); - - word_addr_t i; - for (i = 0; i < num_of_words; i++) { - dest->words[i] = src1->words[i] | src2->words[i]; - } -} - -void bit_array_xor (BIT_ARRAY *dest, BIT_ARRAY *src1, BIT_ARRAY *src2) { - if (dest->num_of_bits != src1->num_of_bits || src1->num_of_bits != src2->num_of_bits) { - // error - fprintf (stderr, "bit_array.c: bit_array_and() : " - "dest, src1 and src2 must be of the same length\n"); - throw; - exit (EXIT_FAILURE); - } - - word_addr_t num_of_words = nwords (src1->num_of_bits); - - word_addr_t i; - for (i = 0; i < num_of_words; i++) { - dest->words[i] = src1->words[i] ^ src2->words[i]; - } -} - -void bit_array_not (BIT_ARRAY *dest, BIT_ARRAY *src) { - if (dest->num_of_bits != src->num_of_bits) { - // error - fprintf (stderr, "bit_array.c: bit_array_and() : " - "dest and src1 must be of the same length\n"); - throw; - exit (EXIT_FAILURE); - } - - word_addr_t num_of_words = nwords (dest->num_of_bits); - - word_addr_t i; - for (i = 0; i < num_of_words; i++) { - dest->words[i] = ~(src->words[i]); - } -} - -inline word_t get_word (BIT_ARRAY *bitarr, word_addr_t word_index, word_addr_t nwords) { - if (word_index >= nwords) { - return 0; - } - - unsigned int offset; - - if (word_index == nwords - 1 && (offset = boffset (bitarr->num_of_bits)) != 0) { - // get masked - return bitarr->words[word_index] & BIT_MASK (offset); - } else { - return bitarr->words[word_index]; - } -} - -// Compare two bit arrays by value stored -// arrays do not have to be the same length (e.g. 101 (5) > 00000011 (3)) -int bit_array_compare (BIT_ARRAY *bitarr1, BIT_ARRAY *bitarr2) { - word_addr_t nwords1 = nwords (bitarr1->num_of_bits); - word_addr_t nwords2 = nwords (bitarr2->num_of_bits); - - word_addr_t max_words = MAX (nwords1, nwords2); - - word_addr_t i; - word_t word1, word2; - - for (i = max_words - 1; i >= 0; i--) { - word1 = get_word (bitarr1, i, nwords1); - word2 = get_word (bitarr2, i, nwords2); - - if (word1 > word2) { - return 1; - } else if (word1 < word2) { - return -1; - } - } - - return 0; -} - -// Return 0 if there was an overflow error, 1 otherwise -char bit_array_add (BIT_ARRAY *dest, BIT_ARRAY *src1, BIT_ARRAY *src2) { - word_addr_t nwords1 = nwords (src1->num_of_bits); - word_addr_t nwords2 = nwords (src2->num_of_bits); - - word_addr_t max_words = MAX (nwords1, nwords2); - - word_addr_t dest_words = nwords (dest->num_of_bits); - - char carry = 0; - - word_addr_t i; - word_t word1, word2; - - for (i = 0; i < max_words; i++) { - word1 = get_word (src1, i, nwords1); - word2 = get_word (src2, i, nwords2); - - word_t result = word1 + word2 + carry; - carry = (result < src1->words[i] || result < src2->words[i]) ? 1 : 0; - - // Check we can store this result - if (i < dest_words - 1) { - dest->words[i] = result; - } else if (i >= dest_words || carry) { - // overflow error - return 0; - } else { - // Check last word (i == dest_words-1) - unsigned int bits_on_last_word = boffset (dest->num_of_bits); - - if (bits_on_last_word > 0 && BIT_MASK (bits_on_last_word) < result) { - // overflow error - return 0; - } - } - } - - if (carry) { - dest->words[max_words] = 1; - } - - // Zero the rest of dest - for (i = max_words + 1; i < dest_words; i++) { - dest->words[i] = 0; - } - - return 1; -} - -bit_index_t bit_array_length (BIT_ARRAY *bit_arr) { return bit_arr->num_of_bits; } - -// If there is an overflow, bit array will be set to all 1s and 0 is returned -// Returns 0 if there was an overflow, 1 otherwise -char bit_array_increment (BIT_ARRAY *bitarr) { - word_addr_t num_of_words = nwords (bitarr->num_of_bits); - - char carry = 1; - - word_addr_t i; - for (i = 0; i < num_of_words - 1; i++) { - if (bitarr->words[i] + 1 < bitarr->words[i]) { - // Carry continues - bitarr->words[i] = 0; - } else { - // Carry is absorbed - bitarr->words[i]++; - carry = 0; - break; - } - } - - // Deal with last word - if (carry) { - unsigned int bits_in_last_word = boffset (bitarr->num_of_bits - 1) + 1; - word_t mask = BIT_MASK (bits_in_last_word); - word_t prev_last_word = bitarr->words[num_of_words - 1] & mask; - - // Increment - bitarr->words[num_of_words - 1] = (prev_last_word + 1) & mask; - - if (prev_last_word == mask) { - // overflow - return 0; - } - } - - return 1; -} - -// If there is an underflow, bit array will be set to all 0s and 0 is returned -// Returns 0 if there was an underflow, 1 otherwise -char bit_array_decrement (BIT_ARRAY *bitarr) { - word_addr_t num_of_words = nwords (bitarr->num_of_bits); - - word_addr_t i; - for (i = 0; i < num_of_words - 1; i++) { - if (bitarr->words[i] > 0) { - bitarr->words[i]--; - - i--; - while (i >= 0) { - bitarr->words[i--] = ULONG_MAX; - } - - return 1; - } - } - - // Must subtract from last word - unsigned int bits_in_last_word = boffset (bitarr->num_of_bits - 1) + 1; - word_t mask = BIT_MASK (bits_in_last_word); - word_t prev_last_word = bitarr->words[num_of_words - 1] & mask; - - if (prev_last_word == 0) { - // underflow - // number unchanged - return 0; - } else { - bitarr->words[num_of_words - 1] = prev_last_word - 1; - return 1; - } -} - -long bit_array_get_long (BIT_ARRAY *bitarr, bit_index_t start) { - // Bounds checking - if (start >= bitarr->num_of_bits) { - fprintf (stderr, "bit_array.c: bit_array_get_long() - out of bounds error " - "(index: %lu, length: %lu)\n", - start, bitarr->num_of_bits); - throw; - exit (EXIT_FAILURE); - } - - word_addr_t num_of_words = nwords (bitarr->num_of_bits); - - word_addr_t word_index = bindex (start); - unsigned int start_offset = boffset (start); - - long result = get_word (bitarr, word_index++, num_of_words) >> start_offset; - - unsigned int offset = WORD_SIZE - start_offset; - - // 64 bits in a long - while (offset < 64) { - result |= get_word (bitarr, word_index++, num_of_words) << offset; - offset += WORD_SIZE; - } - - return result; -} - -int bit_array_get_int (BIT_ARRAY *bitarr, bit_index_t start) { - // Bounds checking - if (start >= bitarr->num_of_bits) { - fprintf (stderr, "bit_array.c: bit_array_get_long() - out of bounds error " - "(index: %lu, length: %lu)\n", - start, bitarr->num_of_bits); - exit (EXIT_FAILURE); - } - - word_addr_t num_of_words = nwords (bitarr->num_of_bits); - - word_addr_t word_index = bindex (start); - unsigned int start_offset = boffset (start); - - int result = ((get_word (bitarr, word_index++, num_of_words) >> start_offset) & UINT_MAX); - - unsigned int offset = WORD_SIZE - start_offset; - - // 32 bits in an int - while (offset < 32) { - result |= (get_word (bitarr, word_index++, num_of_words) << offset) & UINT_MAX; - offset += WORD_SIZE; - } - - return result; -} - -char bit_array_get_char (BIT_ARRAY *bitarr, bit_index_t start) { - // Bounds checking - if (start >= bitarr->num_of_bits) { - fprintf (stderr, "bit_array.c: bit_array_get_long() - out of bounds error " - "(index: %lu, length: %lu)\n", - start, bitarr->num_of_bits); - exit (EXIT_FAILURE); - } - - word_addr_t num_of_words = nwords (bitarr->num_of_bits); - - word_addr_t word_index = bindex (start); - unsigned int start_offset = boffset (start); - - char result = get_word (bitarr, word_index, num_of_words) >> start_offset; - result |= get_word (bitarr, word_index + 1, num_of_words) << (WORD_SIZE - start_offset); - - return result; -} diff --git a/src/bitarray/bit_array.cpp b/src/bitarray/bit_array.cpp index e153633e..12572b05 100644 --- a/src/bitarray/bit_array.cpp +++ b/src/bitarray/bit_array.cpp @@ -136,6 +136,12 @@ void bit_array_set_bit(BIT_ARRAY* bitarr, bit_index_t b) bitarr->words[bindex(b)] |= ((word_t)1 << (boffset(b))); } + +void bit_array_clear_bit_fast(BIT_ARRAY* bitarr, bit_index_t b) +{ + bitarr->words[bindex(b)] &= ~((word_t)1 << (boffset(b))); +} + void bit_array_clear_bit(BIT_ARRAY* bitarr, bit_index_t b) { if(b < 0 || b >= bitarr->num_of_bits) diff --git a/src/bitarray/bit_array.h b/src/bitarray/bit_array.h index 777ed59e..90781160 100644 --- a/src/bitarray/bit_array.h +++ b/src/bitarray/bit_array.h @@ -62,6 +62,9 @@ void bit_array_set_bit(BIT_ARRAY* bitarr, bit_index_t b); // clear a bit (to 0) at position b void bit_array_clear_bit(BIT_ARRAY* bitarr, bit_index_t b); +// clear a bit (to 0) at position b (no error checking) +void bit_array_clear_bit_fast(BIT_ARRAY* bitarr, bit_index_t b); + // Get the value of a bit (returns 0 or 1) char bit_array_get_bit(BIT_ARRAY* bitarr, bit_index_t b); char bit_array_get_bit_nocheck(BIT_ARRAY* bitarr, bit_index_t b); diff --git a/src/conference.cpp b/src/conference.cpp index 876dcab6..b8c9f853 100644 --- a/src/conference.cpp +++ b/src/conference.cpp @@ -2088,7 +2088,7 @@ int selectZmax (int maxzpos, const conference_t::conference_type &ctype, const a maxzpos = al.n_rows - 1; break; default: - printfd ("not implemented...\n"); + printfd ("selectZmax: not implemented for case %d\n", ctype); maxzpos = al.n_rows - 1; } } diff --git a/src/depth_extend.cpp b/src/depth_extend.cpp index 0e856ba8..f02f98bb 100644 --- a/src/depth_extend.cpp +++ b/src/depth_extend.cpp @@ -57,12 +57,12 @@ int extend_array_dynamic ( const arraylist_t &alist, arraydata_t &adx, OAextend std::vector > edata0; extend_array_dynamic ( al, adx, oaextend, earrays0, edata0, dextend, kfinal, Dfinal, rs, verbose ); - #pragma omp critical(extenddynamic) + #pragma omp critical(extenddynamic) { appendArrays ( earrays0, earrays ); edata.insert ( edata.end(), edata0.begin(), edata0.end() ); - + ntotal += dextend.ntotal; nlmc += dextend.nlmc; n += dextend.n; @@ -162,7 +162,7 @@ int extend_array_dynamic ( const array_link &input_array, arraydata_t &adx, OAex oaextendx.use_row_symmetry=rs; oaextendx.nLMC=500000; oaextendx.setAlgorithmAuto ( &adx ); - oaextendx.extendarraymode=OAextend::APPENDEXTENSION; // save memory by only appending the partial array + oaextendx.extendarraymode=OAextend::extendarray_mode_t::APPENDEXTENSION; // save memory by only appending the partial array if ( kfinal=4 ) printf ( "lmc_t %d, col %d\n", v, reduction.lastcol ); @@ -240,12 +240,12 @@ int extend_array_dynamic ( const array_link &input_array, arraydata_t &adx, OAex tmparray.setcolumn ( lastcol, earrays[ii] ); dextend.Deff[ii] = dextend_t::NO_VALUE; - if ( dextend.Dcheck==DCALC_ALWAYS || dextend.lmctype[ii]==LMC_MORE ) { + if ( dextend.Dcheck== dcalc_mode::DCALC_ALWAYS || dextend.lmctype[ii]==LMC_MORE ) { if (adx.is2level() ) dextend.Deff[ii] = Defficiency ( tmparray, verbose>=2 ); else { array_link altmp = hstack(tmparray, earrays[ii] ); - dextend.Deff[ii] =altmp.Defficiencies ()[0]; + dextend.Deff[ii] =altmp.Defficiencies ()[0]; } } @@ -280,7 +280,7 @@ int extend_array_dynamic ( const array_link &input_array, arraydata_t &adx, OAex dextend.DefficiencyFilter ( Dfinal, k, kfinal, Lmax, verbose ); - //int ngood =std::accumulate ( dextend.filter.begin(),dextend.filter.end(),0 ); + //int ngood =std::accumulate ( dextend.filter.begin(),dextend.filter.end(),0 ); int ngoodlmc = countOccurences ( dextend.lmctype, LMC_MORE ); std::vector ctype = dextend.filterArrays ( input_array, earrays, earraysout, edata,verbose ); diff --git a/src/evenodd.cpp b/src/evenodd.cpp index b28d7634..826e7d97 100644 --- a/src/evenodd.cpp +++ b/src/evenodd.cpp @@ -63,7 +63,7 @@ void counter_t::addNumberFound (const counter_t &de) { } } } - + void counter_t::showcountscompact () const { #ifdef DOOPENMP #pragma omp critical @@ -78,14 +78,14 @@ void counter_t::showcountscompact () const { void counter_t::showcounts (const arraydata_t &ad) const { myprintf ("--results--\n"); for (size_t i = ad.strength; i <= (size_t)ad.ncols; i++) { - myprintf ("depth_extend: column %ld: found %d\n", i, this->nfound[i]); + myprintf ("depth_extend: column %ld: found %ld\n", (long)i, (long)this->nfound[i]); } } void counter_t::showcounts (const char *str, int first, int last) const { myprintf ("--results--\n"); for (size_t i = first; i <= (size_t)last; i++) { - myprintf ("%s: column %ld: found %d\n", str, i, this->nfound[i]); + myprintf ("%s: column %ld: found %ld\n", str, (long)i, (long)this->nfound[i]); } } @@ -182,9 +182,6 @@ arraylist_t depth_extend_sub_t::initialize (const arraylist_t &alist, const arra valididx[i] = i; } - if (verbose >= 3) { - myprintf (" v %ld\n", v.size ()); - } return v; } @@ -225,7 +222,7 @@ void processDepth (const arraylist_t &goodarrays, depth_alg_t depthalg, depth_ex OAextend oaextendDirect = dextend.oaextend; oaextendDirect.use_row_symmetry = 1; - oaextendDirect.extendarraymode = OAextend::APPENDFULL; + oaextendDirect.extendarraymode = OAextend::extendarray_mode_t::APPENDFULL; oaextendDirect.checkarrays = 1; depth_extend_hybrid (goodarrays, dextend, extensioncol, oaextendDirect, 1); @@ -355,18 +352,20 @@ void depth_extend_array (const array_link &al, depth_extend_t &dextend, const ar int64_t ps = al.row_symmetry_group ().permsize (); + al.row_symmetry_group ().show(); + depth_alg_t depthalg = DEPTH_DIRECT; depth_extend_sub_t dextendsub; - if (ps > 20000 || ps > 1999999) { + if (ps > 20000 || extensioncol < 9) { // IDEA: enable caching of extensions at later stage! - printfd ("depth_extend_array: warning: direct extension due to large symmetry group!\n"); + printfd ("depth_extend_array: warning: direct extension due to large symmetry group! (size %ld)\n", ps); depthalg = DEPTH_DIRECT; OAextend oaextendDirect = oaextendx; oaextendDirect.use_row_symmetry = 1; - oaextendDirect.extendarraymode = OAextend::APPENDFULL; + oaextendDirect.extendarraymode = OAextend::extendarray_mode_t::APPENDFULL; oaextendDirect.checkarrays = 1; extend_array (al, &adfull, extensioncol, goodarrays, oaextendDirect); @@ -514,7 +513,7 @@ void depth_extend_omp (const arraylist_t &alist, depth_extend_t &dextend, depth_ if (b) { // do lmc check - LMCreduction_t reduction (&adlocal); + LMCreduction_t reduction (&adlocal); LMCreduction_t tmp = reduction; // make sure code is thread safe @@ -769,7 +768,7 @@ void Jcounter::showcompact () const { } myprintf ("; total %ld\n", nt); } - + std::vector< long > Jcounter::getTotalsJvalue (int jval) const { int nmax = maxCols (); std::vector< long > k (nmax + 1); @@ -801,7 +800,7 @@ std::vector< long > Jcounter::getTotals () const { } return k; } - + /// show statistics of the object void Jcounter::show () const { @@ -998,7 +997,7 @@ int compareJ54(const array_link &lhs, const array_link &rhs) { else if(j5l < j5r ) { return 1; } - + std::vector lj4 = lhs5.Jcharacteristics(4); std::vector rj4 = rhs5.Jcharacteristics(4); diff --git a/src/extend.cpp b/src/extend.cpp index 806a1d19..941df61d 100644 --- a/src/extend.cpp +++ b/src/extend.cpp @@ -1,6 +1,4 @@ -#ifdef OAEXTEND_MULTICORE -#include -#endif + #include #include #include @@ -90,13 +88,13 @@ void dextend_t::DefficiencyFilter (double Dfinal, int k, int kfinal, double Lmax int chk = 1; switch (dextend.filtermode) { - case DFILTER_NONE: + case dfilter_t::DFILTER_NONE: chk = 1; break; - case DFILTER_BASIC: + case dfilter_t::DFILTER_BASIC: chk = Ci >= Cfinal; break; - case DFILTER_MULTI: + case dfilter_t::DFILTER_MULTI: // chk= Ci >= Cfinalmulti; chk = Lmaxmulti * Ci >= Cfinal; break; @@ -393,7 +391,7 @@ int check_branch (extend_data_t *es, carray_t *array, extendpos *p, split *stack #else if (1) { #endif - p->value = i; + p->value = i; /* strength t check */ if (valid_element (es, p, array)) { stack->valid[stack->count][npos] = p->value; @@ -457,7 +455,7 @@ int check_branch_2level (extend_data_t *es, carray_t *array_colstart, extendpos #else if (1) { #endif - p->value = i; + p->value = i; /* strength t check */ if (valid_element_2level (es, p)) { stack->valid[stack->count][npos] = p->value; @@ -740,28 +738,18 @@ bool return_stack (split *stack, extendpos *p, array_t *array, int col_offset) { inline void showLoopProgress (array_t *array, const int col_offset, const rowindex_t N, const int node_rank = 0, int nlmcarrays = -1) { -#ifdef _OPENMPX -#define SAFEOMP -#ifdef SAFEOMP - static long _nloops[16] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}; - int tid = omp_get_thread_num (); - _nloops[tid % 16]++; - long nloops = _nloops[tid % 16]; -#else - static long nloops = 0; - nloops++; -#endif -#else + static long nloops = 0; nloops++; -#endif - if (nloops % 10000 == 0) { + if (nloops % 20000 == 0) { fflush (stdout); if (nloops % (500 * 1000 * 1000) == 0) { + if (checkloglevel(QUIET)) { std::cout << "node [" << node_rank << "]: extend loop " << nloops / (1000 * 1000); cout << "m, "; print_perm (array + col_offset, N, 28); + } } } } @@ -786,12 +774,10 @@ int extend_array (const array_link &input_array, const arraydata_t *fullad, cons array_t *array = create_array (fullad->N, ncolsextension); copy_array (origarray, array, fullad->N, extensioncol); -#ifdef OACHECK if (fullad->strength < 1) { log_print (SYSTEM, " extend_array: error: function not defined for strength < 1\n"); throw_runtime_exception("extend_array: strength should be >=1"); } -#endif /* array data */ arraydata_t *ad = new arraydata_t (fullad, ncolsextension); @@ -820,7 +806,6 @@ int extend_array (const array_link &input_array, const arraydata_t *fullad, cons } else { if (oaextend.init_column_previous == INITCOLUMN_J5) { if (extensioncol > 5) { - // printf("extensioncol %d: here\n", extensioncol); init_column_previous (array, p, col_offset, stack, es, oaextend); } else init_column_full (array, p, col_offset, stack, es); @@ -847,11 +832,7 @@ int extend_array (const array_link &input_array, const arraydata_t *fullad, cons #endif const rowindex_t N = p->ad->N; -#ifdef OAEXTEND_MULTICORE - const int node_rank = MPI::COMM_WORLD.Get_rank (); -#else const int node_rank = 0; -#endif LMCreduction_t reduction (ad); reduction.init_state = COPY; @@ -966,15 +947,15 @@ int extend_array (const array_link &input_array, const arraydata_t *fullad, cons } /* the extension found is LMC */ switch (oaextend.extendarraymode) { - case OAextend::APPENDFULL: { + case OAextend::extendarray_mode_t::APPENDFULL: { array_link tmp_extension (array, N, p->col + 1, nlmcarrays); extensions.push_back (tmp_extension); } break; - case OAextend::APPENDEXTENSION: { + case OAextend::extendarray_mode_t::APPENDEXTENSION: { array_link tmp_extension (array + p->col * N, N, 1, nlmcarrays); extensions.push_back (tmp_extension); } break; - case OAextend::STOREARRAY: { + case OAextend::extendarray_mode_t::STOREARRAY: { array_link tmp_extension (array, N, p->col + 1, nlmcarrays); arrayfile_t *storefile = (arrayfile_t *)&oaextend.storefile; // trick to prevent const warnings @@ -986,7 +967,7 @@ int extend_array (const array_link &input_array, const arraydata_t *fullad, cons */ storefile->append_array (tmp_extension); } break; - case OAextend::NONE: { + case OAextend::extendarray_mode_t::NONE: { // do nothing } break; diff --git a/src/extend.h b/src/extend.h index a173a95c..51110ed8 100644 --- a/src/extend.h +++ b/src/extend.h @@ -6,9 +6,6 @@ #ifndef EXTEND_H #define EXTEND_H -#ifdef OAEXTEND_MULTICORE -#include -#endif #include "arraytools.h" #include "lmc.h" #include "oaoptions.h" @@ -43,7 +40,7 @@ class OAextend { int init_column_previous; /// Specification of how to use the generated extensions - enum extendarray_mode_t { + enum class extendarray_mode_t { /// append extension column to extension list APPENDEXTENSION, /// append full array to extension list @@ -56,23 +53,23 @@ class OAextend { /// determines how the extension arrays are stored extendarray_mode_t extendarraymode; - arrayfile_t storefile; + arrayfile_t storefile; // special cases j5structure_t j5structure; private: /// Algorithm mode - algorithm_t algmode; + algorithm_t algmode; public: /** Options for the extension algorithm - * - * + * + * */ OAextend () : singleExtendTime (10.0), nLMC (40000), checkarrays (1), check_maximal (0), use_row_symmetry (1), - init_column_previous (1), extendarraymode (APPENDFULL), j5structure (J5_45), algmode (MODE_AUTOSELECT){}; + init_column_previous (1), extendarraymode (extendarray_mode_t::APPENDFULL), j5structure (J5_45), algmode (MODE_AUTOSELECT){}; /// @copydoc OAextend() OAextend (const OAextend &o) : singleExtendTime (o.singleExtendTime) { this->nLMC = o.nLMC; @@ -87,13 +84,13 @@ class OAextend { // we do not copy the storefile: this->storefile = o.storefile; }; /** @copydoc OAextend() - * + * * The algorithm is automatically determined from the specified arrayclass. - * + * */ OAextend (arraydata_t &arrayclass) : singleExtendTime (10.0), nLMC (40000), checkarrays (1), check_maximal (0), use_row_symmetry (1), - init_column_previous (1), extendarraymode (APPENDFULL), j5structure (J5_45), algmode (MODE_AUTOSELECT) { + init_column_previous (1), extendarraymode (extendarray_mode_t::APPENDFULL), j5structure (J5_45), algmode (MODE_AUTOSELECT) { setAlgorithmAuto (&arrayclass); }; /// Set the algorithm to use for LMC checks @@ -162,7 +159,7 @@ struct extendpos { * \param array_class Class of arrays to generate * \param oaextend_options Parameters for the extension algorithm * \return List of all generated arrays -* +* * @see extend_array(const array_link &, arraydata_t &, OAextend const &) */ arraylist_t extend_arraylist (const arraylist_t &array_list, arraydata_t &array_class, OAextend const &oaextend_options); @@ -174,7 +171,7 @@ arraylist_t extend_arraylist (const arraylist_t &array_list, arraydata_t &array_ arraylist_t extend_arraylist (const arraylist_t &array_list, const arraydata_t &arrayclass); /** @copydoc extend_arraylist(const arraylist_t &, arraydata_t &, OAextend const &) - * + * * \param extensioncol Index of column to be added to the designs * \param extensions List to append generated designs to * \return Number of candidate arrays generated @@ -211,14 +208,14 @@ arraylist_t extend_array (const array_link &array, arraydata_t &arrayclass); int extend_array (const array_link &array, const arraydata_t *arrayclass, const colindex_t extension_column, arraylist_t &extensions, OAextend const &oaextend); -/** Run the LMC extension algorithm starting with the root array +/** Run the LMC extension algorithm starting with the root array * * @see extend_array(const array_link &, arraydata_t &, OAextend const &) */ arraylist_t runExtendRoot (arraydata_t arrayclass, int max_number_columns, int verbose = 0); -enum dfilter_t { +enum class dfilter_t { /// no filtering on D-efficiency DFILTER_NONE, /// filtering on D-efficiency @@ -226,7 +223,7 @@ enum dfilter_t { /// filtering on D-efficiency with multi column prediction DFILTER_MULTI }; -enum dcalc_mode { +enum class dcalc_mode { /// always calculate efficiency DCALC_ALWAYS, /// only calculate efficiency for LMC_LESS @@ -257,7 +254,7 @@ struct dextend_t { /// perform immediate LMC check in extension int directcheck; - dextend_t () : filtermode (DFILTER_MULTI), Dcheck (DCALC_COND), directcheck (1){}; + dextend_t () : filtermode (dfilter_t::DFILTER_MULTI), Dcheck (dcalc_mode::DCALC_COND), directcheck (1){}; void resize (int nn) { this->lmctype.resize (nn); diff --git a/src/lmc.cpp b/src/lmc.cpp index 9018d6b3..6eb9845d 100644 --- a/src/lmc.cpp +++ b/src/lmc.cpp @@ -1,5 +1,6 @@ #include #include +#include #include "arraytools.h" #include "extend.h" @@ -7,7 +8,6 @@ #include "mathtools.h" #include "oaoptions.h" #include "tools.h" - #include "nonroot.h" #ifndef DOOPENMP @@ -21,7 +21,8 @@ bool operator!= (symmdataPointer const &ptr, int x) { // A friendly reminder to not test pointers against any values except 0 (NULL) myassert (!x, "invalid pointer in comparison"); - return ptr; + bool nonzero = valid_ptr(ptr); + return nonzero; } #else @@ -30,6 +31,14 @@ bool operator!= (symmdataPointer const &ptr, int x) { #endif #endif +bool valid_ptr(const symmdataPointer& sd) { +#ifdef SDSMART + return sd.get() != 0; +#else + return sd != 0; +#endif +} + #ifdef DOOPENMP #include #endif @@ -282,7 +291,7 @@ LMCreduction_t::LMCreduction_t (const LMCreduction_t &at) { array = create_array (transformation->ad); } -LMCreduction_t &LMCreduction_t::operator= (const LMCreduction_t &at) +LMCreduction_t &LMCreduction_t::operator= (const LMCreduction_t &at) { mode = at.mode; state = at.state; @@ -296,7 +305,7 @@ LMCreduction_t &LMCreduction_t::operator= (const LMCreduction_t &at) ncols = at.ncols; nrows = at.nrows; - staticdata = at.staticdata; + staticdata = at.staticdata; free (); @@ -489,9 +498,9 @@ void deallocate_rowsort(rowsort_t *& rowsort) { rowsorter_t::rowsorter_t(int number_of_rows) { this->number_of_rows = number_of_rows; - this->rowsort = allocate_rowsort(number_of_rows); + this->rowsort = allocate_rowsort(number_of_rows); this->reset_rowsort(); - } + } void rowsorter_t::reset_rowsort() { @@ -502,19 +511,19 @@ void rowsorter_t::reset_rowsort() } rowsorter_t::~rowsorter_t() { - deallocate_rowsort(this->rowsort); + deallocate_rowsort(this->rowsort); } dyndata_t::dyndata_t (int N_, int col_) { this->N = N_; this->col = col_; this->rowsort = allocate_rowsort(this->N); - this->colperm = new_perm_init< colindex_t > (this->N); + this->colperm = new_perm_init< colindex_t > (this->N); this->rowsortl = 0; for (int i = 0; i < N; i++) { this->rowsort[i].r = i; - this->rowsort[i].val = 0; + this->rowsort[i].val = 0; } } @@ -555,7 +564,7 @@ void dyndata_t::initdata (const dyndata_t &dd) { memcpy (this->rowsort, dd.rowsort, N * sizeof (rowsort_t)); } if (dd.rowsortl != 0) { - this->rowsortl = new_perm< rowindex_t > (this->N); + this->rowsortl = new_perm< rowindex_t > (this->N); copy_perm (dd.rowsortl, this->rowsortl, this->N); } } @@ -590,7 +599,7 @@ void dyndata_t::copydata (const dyndata_t &dd) { } if (dd.rowsortl != 0) { deleterowsortl (); - this->rowsortl = new_perm< rowindex_t > (this->N); + this->rowsortl = new_perm< rowindex_t > (this->N); copy_perm (dd.rowsortl, this->rowsortl, this->N); } } @@ -707,12 +716,18 @@ void create_root_permutations_index_helper (rowperm_t *rperms, levelperm_t *lper * @return */ rowperm_t *create_root_permutations_index (const arraydata_t *ad, int &totalpermsr) { - /* OPTIMIZE: can this function be made simpler by noting that all permutations of N/oaindex occur ??*/ + /* OPTIMIZE: can this function be made simpler by noting that all permutations of N/oaindex occur */ totalpermsr = 1; - for (int i = 0; i < ad->strength; i++) - totalpermsr *= factorial< int > (ad->s[i]); - + for (int i = 0; i < ad->strength; i++) { + int f = factorial< int > (ad->s[i]); + if (totalpermsr>std::numeric_limits::max()/f) { + myprintf("create_root_permutations_index: design specification: %s\n", ad->fullidstr().c_str()); + throw_runtime_exception(printfstring("create_root_permutations_index: integer overflow for specified factor levels")); + } + totalpermsr *= f; + } + log_print (DEBUG, "create_root_permutations_index: allocating rootrowperms table of size " "%d*%d*sizeof(rowsort_t)=%ld=%.1f MB (sizeof(rowsort_t)=%d)\n", totalpermsr, ad->N, (long)totalpermsr * (long)ad->N * (long)sizeof (rowsort_t), @@ -765,7 +780,7 @@ rowperm_t *create_root_permutations_index_full (const arraydata_t *ad, int &tota // loop over all columns permutations for (int j = 0; j < pow (double(2), ad->strength); j++) { // loop over all possible level permutations - root_row_permutation_from_index (j, ad, lperms); + root_row_permutation_from_index (j, ad, lperms); cperm_lperms_to_rowsort (tmprperm, lperms, rperms[permcounter], cp, ad); permcounter++; @@ -1173,14 +1188,14 @@ lmc_t LMCreduce_root_level_perm (array_t const *original, const arraydata_t *ad, /* pass to non-root stage */ if (oaextend.getAlgorithm () == MODE_LMC_2LEVEL || - (oaextend.getAlgorithm () == MODE_J5ORDERX && reduction->sd != 0) || - (oaextend.getAlgorithm () == MODE_J5ORDER_2LEVEL && reduction->sd != 0)) { + (oaextend.getAlgorithm () == MODE_J5ORDERX && valid_ptr(reduction->sd)) || + (oaextend.getAlgorithm () == MODE_J5ORDER_2LEVEL && valid_ptr(reduction->sd))) { dyndatatmp.initrowsortl (); ret = LMCreduce_non_root_2level (original, ad, &dyndatatmp, reduction, oaextend, - tmpStatic); + tmpStatic); } else ret = LMCreduce_non_root (original, ad, &dyndatatmp, reduction, oaextend, - tmpStatic); + tmpStatic); if (ret == LMC_LESS) { break; @@ -1432,7 +1447,10 @@ std::vector< int > symmetrygroup2splits (const symmetry_group &sg, int ncols, in return splits; } -void _calculate_j4_values(std::vector &j4_values, carray_t *array, const int N, const colperm_t comb) +/** Calculate J4 values for all delete-one-factor combinations of the specified comb 5-column combination + * +*/ +void _calculate_dof_j4_values(std::vector &j4_values, carray_t *array, const int N, const colperm_t comb) { colindex_t lc[4]; init_perm(lc, 4); @@ -1462,8 +1480,7 @@ int jj45split (carray_t *array, rowindex_t N, int jj, const colperm_t comb, cons // allocate buffer to hold the values std::vector< int > j4_values (5); - /* calculate the J4 values */ - _calculate_j4_values(j4_values, array, N, comb); + _calculate_dof_j4_values(j4_values, array, N, comb); indexsort s (5); s.sort (j4_values); @@ -1480,12 +1497,7 @@ int jj45split (carray_t *array, rowindex_t N, int jj, const colperm_t comb, cons // create the sorted column combination perform_inv_perm (comb, comby, jj, s.indices); - if (verbose >= 2) { - myprintf ("initial comb: "); - print_perm (comb, jj); - myprintf (" result comb: "); - print_perm (comby, jj); - } + // make splits init_perm (perm, ad.ncols); create_perm_from_comb2< colindex_t > (perm, comby, jj, ad.ncols); @@ -1496,7 +1508,7 @@ int jj45split (carray_t *array, rowindex_t N, int jj, const colperm_t comb, cons } dyndata_t dyndata (ad.N); - dyndata.col = 0; + dyndata.col = 0; dyndata.setColperm (perm, ad.ncols); if (oaextend.getAlgorithm () == MODE_J5ORDERX || oaextend.getAlgorithm () == MODE_J5ORDER_2LEVEL) { @@ -1606,7 +1618,7 @@ jj45_t jj45val (carray_t *array, rowindex_t N, const colperm_t comb, int dosort int ii = 5 - i - 1; fastJupdate (array, N, 1, comb + ii, tmpval); ww[i + 1] = abs (fastJupdateValue (N, tmpval)); - fastJupdate (array, N, 1, comb + ii, tmpval); + fastJupdate (array, N, 1, comb + ii, tmpval); } if (dosort) { @@ -1680,7 +1692,7 @@ lmc_t LMCcheckj5 (array_link const &al, arraydata_t const &adin, LMCreduction_t lmc_t ret = LMC_EQUAL; - /* loop over all possible column combinations for the first jj columns */ + /* calculate invariants for the first jj columns */ int jbase = abs (jvaluefast (array, ad.N, jj, firstcolcomb)); jj45_t j54_base = jj45val (array, ad.N, firstcolcomb, 0); @@ -1691,7 +1703,7 @@ lmc_t LMCcheckj5 (array_link const &al, arraydata_t const &adin, LMCreduction_t const int orig = oaextend.j5structure == J5_ORIGINAL; int ncolsfirst = adin.colgroupsize[0]; - int nc = ncombs (ncolsfirst, jj); + int nc = ncombs (ncolsfirst, jj); if (dverbose) { myprintf ("LMCcheckj5: selected ncolsfirst %d, nc %d, jbase %d, wbase %f\n", ncolsfirst, nc, jbase, @@ -1772,7 +1784,7 @@ lmc_t LMCcheckj5 (array_link const &al, arraydata_t const &adin, LMCreduction_t break; } - next_combination (combroot, ad.strength, jj); + next_combination (combroot, ad.strength, jj); } } int retx = ret; @@ -1787,7 +1799,7 @@ lmc_t LMCcheckj5 (array_link const &al, arraydata_t const &adin, LMCreduction_t break; } - next_combination (firstcolcomb, jj, ncolsfirst); + next_combination (firstcolcomb, jj, ncolsfirst); } delete_perm (perm); delete_perm (pp); @@ -2071,7 +2083,7 @@ lmc_t LMCreduction_train(const array_link &al, const arraydata_t *ad, LMCreducti /// full reduction, no root-trick lmc_t LMCreduceFull (carray_t *original, const array_t *array, const arraydata_t *adx, const dyndata_t *dyndata, LMCreduction_t *reduction, const OAextend &oaextend, LMCreduction_helper_t &tmpStatic) { - arraydata_t *ad = new arraydata_t (*adx); + arraydata_t *ad = new arraydata_t (*adx); ad->oaindex = ad->N; // NOTE: this is to prevent processing on blocks in LMC_check_col if (dyndata->col == 0) { @@ -2302,7 +2314,7 @@ inline int predictJrowsort (const array_t *array, const int N, const rowsort_t * /// select the unique arrays in a list, the original list is sorted in place void selectUniqueArrays (arraylist_t &input_arrays, arraylist_t &output_arrays, int verbose) { - sort (input_arrays.begin (), input_arrays.end ()); + sort (input_arrays.begin (), input_arrays.end ()); if (input_arrays.size () > 0) { std::vector< int > vv (input_arrays.size ()); vv[0] = 1; @@ -2359,8 +2371,8 @@ array_link reduceDOPform (const array_link &al, int verbose) { } -/** Calculate mixed-level projection GWLP values from normal GWLP values. - * +/** Calculate mixed-level projection GWLP values from normal GWLP values. + * * These are the normal projection values, with added the factor level of the removed column. * */ @@ -2387,10 +2399,10 @@ std::vector< GWLPvalue > mixedProjGWLP (const std::vector< GWLPvalue > dopgwp, c } std::vector< GWLPvalue > projectionDOFvalues (const array_link &array, int verbose ) { - + arraydata_t arrayclass=arraylink2arraydata(array); std::vector< GWLPvalue > projection_dof_values = projectionGWLPs(array); - + if (arrayclass.ismixed() ) { projection_dof_values = mixedProjGWLP(projection_dof_values, arrayclass, verbose); } @@ -2420,7 +2432,7 @@ else { } } -bool _check_dof_order_minimal(std::vector &dopgwp, int ncols, int verbose) +bool _check_dof_order_minimal(std::vector &dopgwp, int ncols, int verbose) { GWLPvalue x = *(min_element(dopgwp.begin(), dopgwp.begin() + ncols - 1)); if (verbose >= 2) { @@ -2490,7 +2502,7 @@ array_transformation_t reductionDOP (const array_link &input_array, int verbose) reduction.init_state = INIT; reduction.setArray (alf); int changed = check_root_update (alf.array, ad, reduction.array); - copy_array (alf.array, reduction.array, input_array.n_rows, input_array.n_columns); + copy_array (alf.array, reduction.array, input_array.n_rows, input_array.n_columns); } lmc_t ret = LMCcheck (alf, ad, oaextend, reduction); @@ -2578,7 +2590,7 @@ void reduceArraysGWLP (const arraylist_t &input_arrays, arraylist_t &output_arra printfd ("reduceArraysGWLP: ret %d (LMC_MORE %d, strength %d)\n", ret, LMC_MORE, ad.strength); helper_compare_dof_reductions(alf, lm, verbose, ret); - + xlist.push_back (lm); } if (verbose) { @@ -2590,7 +2602,7 @@ void reduceArraysGWLP (const arraylist_t &input_arrays, arraylist_t &output_arra selectUniqueArrays (xlist, output_arrays, verbose); } else { - sort (xlist.begin (), xlist.end ()); + sort (xlist.begin (), xlist.end ()); for (size_t i = 0; i < xlist.size (); i++) { output_arrays.push_back (xlist[i]); diff --git a/src/lmc.h b/src/lmc.h index 32e4d8e7..174b459e 100644 --- a/src/lmc.h +++ b/src/lmc.h @@ -87,10 +87,10 @@ enum algorithm_t { const algorithm_t MODE_ORIGINAL = MODE_LMC; -/// Return string representation of available algorithm modes +/// Return string representation of available algorithm modes inline std::string algorithm_t_list () { std::string ss = - printfstring ("%d (automatic), %d (original), %d (check j4), %d (j5 order), %d (j5 order dominant), %d " + printfstring ("%d (automatic), %d (original), %d (check j4), %d (j5 order), %d (MODE_J5ORDERX, j5 order dominant), %d " "(MODE_J5ORDERXFAST)", MODE_AUTOSELECT, MODE_ORIGINAL, MODE_J4, MODE_J5ORDER, MODE_J5ORDERX, MODE_J5ORDER_2LEVEL); ss += printfstring (", %d (MODE_LMC_SYMMETRY), %d (MODE_LMC_2LEVEL)", MODE_LMC_SYMMETRY, MODE_LMC_2LEVEL); @@ -110,7 +110,7 @@ enum initcolumn_t { enum j5structure_t { /// Ordering based in J5 in succesive columns J5_ORIGINAL, - /// Ordering based on J5 and the 5-tuple of J4 values + /// Ordering based on J5 and the 5-tuple of J4 values. Also called the L5 ordering J5_45 }; /// return name of the algorithm @@ -221,7 +221,7 @@ struct LMCreduction_helper_t { rootrowperms = this->rootrowperms; lperm_p = this->current_trans->lperms; - this->LMC_root_rowperms_init = 1; + this->LMC_root_rowperms_init = 1; } /** @brief Static initialization of root row permutations (full group) @@ -248,14 +248,14 @@ void clear_LMCreduction_pool (); enum REDUCTION_STATE { /// the reduction is equal to the initial REDUCTION_INITIAL, - /// the reduction was changed + /// the reduction was changed REDUCTION_CHANGED }; //! main mode for the LMC routine: test, reduce or reduce with initialization enum OA_MODE { /// test for minimal form OA_TEST, /// reduce to minimal form - OA_REDUCE, + OA_REDUCE, /// reduce to partial minimal form OA_REDUCE_PARTIAL }; @@ -285,17 +285,6 @@ typedef std::shared_ptr< symmdata > symmdataPointer; #include typedef std::shared_ptr< symmdata > symmdataPointer; -/** @brief An operator to test whether a pointer holds a null value or not - * - * It is recomended to not compare pointers to integers and use explicit - * conversion `ptr to bool` instead or as another approach to compare pointers against - * nullptr which has type of std::nullptr. But such methods would force migration to C++11. Hence, - * the purpose of this method is to avoid ether rewritting the whole project or changing the code - * that already works on other platforms (aka Windows) by adding a missing in header - * comparison `bool operator!=(smart_ptr, int)` - */ -bool operator!= (symmdataPointer const &ptr, int x); - #else #include typedef std::tr1::shared_ptr< symmdata > symmdataPointer; @@ -305,6 +294,9 @@ typedef std::tr1::shared_ptr< symmdata > symmdataPointer; typedef symmdata *symmdataPointer; #endif +/** Return true if the (smart) symmdataPointer pointer is allocated */ +bool valid_ptr(const symmdataPointer& sd); + /// initial state for reduction algorithm enum INIT_STATE { // invalid state @@ -338,7 +330,7 @@ struct LMCreduction_t { OA_MODE mode; REDUCTION_STATE state; - INIT_STATE init_state; + INIT_STATE init_state; //! maximum depth for search tree int maxdepth; @@ -349,8 +341,8 @@ struct LMCreduction_t { //! counter for number of reductions made long nred; - int targetcol; - int mincol; + int targetcol; + int mincol; int nrows, ncols; @@ -386,7 +378,7 @@ struct LMCreduction_t { #else symmdata *sdp = sd; #endif - if (sdp != 0 && cache) { + if (sdp != 0 && cache) { } else { // update symmetry data this->sd = symmdataPointer (new symmdata (al, 1)); @@ -441,7 +433,7 @@ class rowsorter_t public: int number_of_rows; rowsort_t *rowsort; - + rowsorter_t(int number_of_rows); ~rowsorter_t(); @@ -549,7 +541,7 @@ struct dyndata_t { void initdata (const dyndata_t &dd); }; -/** Return True if the array is in root form +/** Return True if the array is in root form * * \param array Array to check * \param strength Strength to use @@ -591,10 +583,10 @@ lmc_t LMCcheckOriginal (const array_link &array); void reduceArraysGWLP (const arraylist_t &input_arrays, arraylist_t &reduced_arrays, int verbose, int dopruning = 1, int strength = 2, int dolmc = 1); -/** Caculate the transformation reducing an array to delete-on-factor normal - * +/** Caculate the transformation reducing an array to delete-on-factor normal + * * The normal form is described in "A canonical form for non-regular arrays based on generalized wordlength pattern values of delete-one-factor projections", Eendebak, 2014 - * + * * \param array Orthogonal array * \param verbose Verbosity level * \returns The transformation that reduces the array to normal form @@ -602,9 +594,9 @@ void reduceArraysGWLP (const arraylist_t &input_arrays, arraylist_t &reduced_arr array_transformation_t reductionDOP (const array_link &array, int verbose = 0); /** Reduce an array to canonical form using delete-1-factor ordering - * + * * The normal form is described in "A canonical form for non-regular arrays based on generalized wordlength pattern values of delete-one-factor projections", Eendebak, 2014 - * + * * \param array Orthogonal array * \param verbose Verbosity level * \returns The array transformed to normal form @@ -632,7 +624,7 @@ lmc_t LMCcheckLex(array_link const &array, arraydata_t const &arrayclass); lmc_t LMCcheckj4 (array_link const &array, arraydata_t const &arrayclass, LMCreduction_t &reduction, const OAextend &oaextend, int jj = 4); -/// Perform minimal form check for J5 ordering +/// Perform minimal form check for J5 ordering (in the paper this is called the L5 ordering) lmc_t LMCcheckj5 (array_link const &array, arraydata_t const &arrayclass, LMCreduction_t &reduction, const OAextend &oaextend); @@ -643,5 +635,3 @@ lmc_t LMCcheckj5 (array_link const &array, arraydata_t const &arrayclass, LMCred */ void print_rowsort (rowsort_t *rowsort, int N); void print_column_rowsort (const array_t *arraycol, rowsort_t *rowsort, int N); - - diff --git a/src/mathtools.cpp b/src/mathtools.cpp index 7aeadf8e..90d4b3e8 100644 --- a/src/mathtools.cpp +++ b/src/mathtools.cpp @@ -1,7 +1,20 @@ #include +#include #include "mathtools.h" +std::string printfstring(const char* message, ...) { + char buf[64 * 1024]; + + va_list va; + va_start(va, message); + vsprintf(buf, message, va); + va_end(va); + + std::string str(buf); + return str; +} + void set_srand (unsigned int s) { srand (s); } template < class Type > void symmetry_group::init (const std::vector< Type > vals, bool ascendingx, int verbose) { diff --git a/src/mathtools.h b/src/mathtools.h index 7eb8ea6b..0df42dce 100644 --- a/src/mathtools.h +++ b/src/mathtools.h @@ -41,6 +41,12 @@ inline double round (double x) { return floor (x + 0.5); } #include +/** + * @brief Function similar to printf returning C++ style string + * @param message + * @return + */ +std::string printfstring(const char* message, ...); /** Class to make a pool of objects that can be re-used * @@ -244,7 +250,7 @@ struct mvalue_t { /// vector containing the values std::vector< NumericType > values; // int direction; - enum direction_t { + enum class direction_t { /// Order from high to low values HIGH, /// Order from low to high values @@ -256,7 +262,7 @@ struct mvalue_t { * * The object consists of a vector of elements. */ - mvalue_t () : ordering (HIGH){}; + mvalue_t () : ordering (direction_t::HIGH){}; ~mvalue_t (){}; /** @copydoc mvalue_t::mvalue_t() @@ -264,7 +270,7 @@ struct mvalue_t { * \param element Single element to add to the vector * \param dd Ordering to use */ - mvalue_t (NumericType element, direction_t dd = HIGH) { + mvalue_t (NumericType element, direction_t dd = direction_t::HIGH) { values.push_back (element); ordering = dd; } @@ -273,7 +279,7 @@ struct mvalue_t { * \param elements Vector to use for initalization of the object * \param dd Ordering to use */ - mvalue_t (std::vector< NumericType > elements, direction_t dd = HIGH) { + mvalue_t (std::vector< NumericType > elements, direction_t dd = direction_t::HIGH) { ordering = dd; this->values = elements; } @@ -283,7 +289,7 @@ struct mvalue_t { * \param elements Vector to use for initalization of the object * \param dd Ordering to use */ - template < class T > mvalue_t(std::vector< T > elements, direction_t dd = HIGH) { + template < class T > mvalue_t(std::vector< T > elements, direction_t dd = direction_t::HIGH) { ordering = dd; values.clear(); values.resize(elements.size()); @@ -331,7 +337,7 @@ struct mvalue_t { bool operator< (const mvalue_t &rhs) const { bool val = 0; - if (ordering == HIGH) + if (ordering == direction_t::HIGH) val = (bool)worse (rhs); else val = (bool)better (rhs); @@ -339,7 +345,7 @@ struct mvalue_t { } bool operator> (const mvalue_t &rhs) const { bool val = 0; - if (ordering == HIGH) + if (ordering == direction_t::HIGH) val = (bool)better (rhs); else val = (bool)worse (rhs); @@ -539,29 +545,42 @@ static void print_perm (std::ostream &out, const larray< permutationType > s, co } } - -/// Print permutation to output stream +/// Print permutation to string template < class permutationType > /* permtype should be a numeric type, i.e. int or long */ -static void print_perm_int(const permutationType*s, const int len, const int maxlen = 256, const bool ret = true) { +std::string permutation2string(const permutationType* s, const int len, const int maxlen = 256, const bool add_return = true) { int plen = std::min(len, maxlen); - myprintf("{"); + std::string str = printfstring("{"); for (int i = 0; i < plen - 1; i++) - myprintf("%d,", s[i]); + str += printfstring("%d,", s[i]); if (len == 0) { - myprintf("}"); + str += printfstring("}"); } else { if (plen < len) - myprintf("%d,...", s[plen - 1]); + str +=printfstring("%d,...}", s[plen - 1]); else - myprintf("%d}", s[plen - 1]); + str +=printfstring("%d}", s[plen - 1]); } - if (ret) { - myprintf("\n"); + if (add_return) { + str +=printfstring("\n"); } + return str; +} + +/// Print permutation to string +template < class permutationType > /* permtype should be a numeric type, i.e. int or long */ +std::string permutation2string(const std::vector< permutationType > s, const int maxlen = 256, const bool add_return = true) { + return permutation2string(s.data(), s.size(), maxlen, add_return); +} + +/// Print permutation to output stream +template < class permutationType > /* permtype should be a numeric type, i.e. int or long */ +static void print_perm_int(const permutationType*s, const int len, const int maxlen = 256, const bool ret = true) { + std::string str = permutation2string(s, len, maxlen, ret); + myprintf("%s", str.c_str()); } /// Print permutation to output stream diff --git a/src/nonroot.cpp b/src/nonroot.cpp index 57d3fda5..41b0b0c2 100644 --- a/src/nonroot.cpp +++ b/src/nonroot.cpp @@ -396,7 +396,7 @@ lmc_t LMC_check_col_ft_2level (const array_t *originalcol, const array_t *arrayc const int scol = dd->col - 1; rowsortlight_t *rowsort = dd->rowsortl; - int nb = sd.ft.atfast (sd.ft.n_rows - 1, scol); + int nb = sd.ft.atfast (sd.ft.n_rows - 1, scol); array_t *symmdata_column_pointer = sd.ft.array + scol * sd.ft.n_rows; /* we check in blocks determined by the ft */ @@ -612,7 +612,7 @@ inline lmc_t LMC_check_col_j5order (const array_t *original, const array_t *arra assert (ad->order == ORDER_J5); const int cpoffset = ad->N * dd->colperm[dd->col]; - if (dd->col < 4) { + if (dd->col < 4) { lmc_t ret = LMC_check_col (original + dd->col * +ad->N, array + cpoffset, lperm, ad, dd); return ret; } @@ -630,7 +630,7 @@ inline lmc_t LMC_check_col_j5order (const array_t *original, const array_t *arra for (int x = 0; x < 5; x++) pp[x] = dd->colperm[x]; pp[4] = dd->colperm[dd->col]; - int jcol = abs (jvaluefast (array, ad->N, 5, pp)); + int jcol = abs (jvaluefast (array, ad->N, 5, pp)); if ( dd->col > 4) { if (log_print (DEBUG, "")) { myprintf (" xxx col %d, jbase %d, jcol %d: colperm: ", dd->col, jbase, jcol); @@ -674,13 +674,13 @@ lmc_t LMCreduce_non_root_j4 (const array_t *original, const arraydata_t *ad, con helper_structure.init_nonroot_stage (lperm_p, colperm_p, localcolperm_p, dynd_p, dynd_p_nelem, colbuffer, ad); - const int number_remaining = ad->ncols - dyndata->col; - const int number_levels = ad->s[dyndata->col]; + const int number_remaining = ad->ncols - dyndata->col; + const int number_levels = ad->s[dyndata->col]; levelperm_t lperm = lperm_p[dyndata->col]; colperm_t colpermloop = localcolperm_p[dyndata->col]; dyndata_t *dyndatacpy = dynd_p[dyndata->col]; - const int current_column_group = ad->get_col_group (dyndata->col); + const int current_column_group = ad->get_col_group (dyndata->col); /* number of columns remaining in this group */ const int ncolsremgroup = ad->colgroupsize[current_column_group] + ad->colgroupindex[current_column_group] - dyndata->col; @@ -806,7 +806,7 @@ lmc_t LMCreduce_non_root (const array_t *original, const arraydata_t *arrayclass helper_structure.init_nonroot_stage (lperm_p, colperm_p, localcolperm_p, dynd_p, dynd_p_nelem, colbuffer, arrayclass); - const int number_remaining = arrayclass->ncols - dyndata->col; + const int number_remaining = arrayclass->ncols - dyndata->col; const int nlevels = arrayclass->s[dyndata->col]; /* number of levels at this column */ levelperm_t lperm = lperm_p[dyndata->col]; @@ -838,7 +838,7 @@ lmc_t LMCreduce_non_root (const array_t *original, const arraydata_t *arrayclass for (int j = 0; j < nlevelperms; j++) { /* copy dyndata rowsort data */ cpy_dyndata_rowsort (dyndata, dyndatacpy); - dyndatacpy->col = dyndata->col; + dyndatacpy->col = dyndata->col; /* LMC_check_col performs level permutations, updates the sorting structure and compares with * the original array on blocks of oaindex */ @@ -863,16 +863,16 @@ lmc_t LMCreduce_non_root (const array_t *original, const arraydata_t *arrayclass if ((oaextend.getAlgorithm () == MODE_LMC_2LEVEL || oaextend.getAlgorithm () == MODE_J5ORDER_2LEVEL) && - reduction->sd != 0) { - myassert (reduction->sd != 0, "LMC_check_col_ft"); + valid_ptr(reduction->sd)) { + myassert (valid_ptr(reduction->sd), "LMC_check_col_ft: no valid pointer"); ret = LMC_check_col_ft_2level ( reduction->array + dyndata->col * +arrayclass->N, original + cpoffset, - lperm, arrayclass, dyndatacpy, *(reduction->sd)); + lperm, arrayclass, dyndatacpy, *(reduction->sd)); } else { ret = LMC_check_col (reduction->array + dyndata->col * +arrayclass->N, original + cpoffset, lperm, arrayclass, dyndatacpy); } - dyndatacpy->col = dyndata->col; + dyndatacpy->col = dyndata->col; } else { ret = @@ -970,7 +970,7 @@ lmc_t LMCreduce_non_root_2level (const array_t *original, const arraydata_t *ad, dyndatacpy->allocate_rowsortl (); dyndatacpy->col = dyndata->col; - const int current_column_group = ad->get_col_group (dyndata->col); + const int current_column_group = ad->get_col_group (dyndata->col); /* number of columns remaining in this group */ const int ncolsremgroup = ad->colgroupsize[current_column_group] + ad->colgroupindex[current_column_group] - dyndata->col; @@ -992,9 +992,9 @@ lmc_t LMCreduce_non_root_2level (const array_t *original, const arraydata_t *ad, /* loop over all level permutations */ for (int j = 0; j < nlevelperms; j++) { /* copy dyndata rowsort data */ - cpy_dyndata_rowsortl (dyndata, dyndatacpy); + cpy_dyndata_rowsortl (dyndata, dyndatacpy); - dyndatacpy->col = dyndata->col; + dyndatacpy->col = dyndata->col; /* LMC_check_col performs level permutations, updates the sorting structure and compares with * the original array on blocks of oaindex */ @@ -1014,7 +1014,7 @@ lmc_t LMCreduce_non_root_2level (const array_t *original, const arraydata_t *ad, original + cpoffset, lperm, ad, dyndatacpy, *(reduction->sd), 0); } - dyndatacpy->col = dyndata->col; + dyndatacpy->col = dyndata->col; } else { ret = @@ -1041,7 +1041,7 @@ lmc_t LMCreduce_non_root_2level (const array_t *original, const arraydata_t *ad, oaextend, tmpStatic); } } - + } if (ret == LMC_LESS) { break; // break from levelperms loop diff --git a/src/tools.cpp b/src/tools.cpp index 40f25230..422a4f87 100644 --- a/src/tools.cpp +++ b/src/tools.cpp @@ -197,7 +197,24 @@ double get_time_ms () { #include #include -double get_time_ms (double t0) { +double get_time_ms (double t0) +{ +#ifdef WIN32 + struct timeb tb; + ftime (&tb); + return (double)tb.time + ((double)tb.millitm / 1000.0f) - t0; +#else + struct timespec ts; + + if (clock_gettime (CLOCK_MONOTONIC, &ts) == 0) { + return (double) (double(ts.tv_sec) + double(ts.tv_nsec ) / 1000000000) - t0; + } + else + return 0; +#endif +} + +double get_time_ms2 (double t0) { const time_t epoch = 0; time_t now; time(&now); diff --git a/src/tools.h b/src/tools.h index 74d8ec15..1f6f1593 100644 --- a/src/tools.h +++ b/src/tools.h @@ -229,7 +229,6 @@ DataType **malloc2d_irr (const int nrows, const rtype *rowsizes, int &nelements) template < class DataType, class numtype > DataType **malloc2d (const numtype nrows, const int rowsize) { DataType **data; - data = new DataType *[nrows]; if (data == 0) { throw_runtime_exception ("malloc2d: error with memory allocation"); diff --git a/src/version.h b/src/version.h index 364a4e9c..20e5c409 100644 --- a/src/version.h +++ b/src/version.h @@ -1 +1 @@ -const char *__version__ = "2.7.1"; +const char *__version__ = "2.7.2"; diff --git a/tests/test_cpplibrary.py b/tests/test_cpplibrary.py index 21fb4520..2840e7dd 100644 --- a/tests/test_cpplibrary.py +++ b/tests/test_cpplibrary.py @@ -6,6 +6,7 @@ import os import sys import tempfile +import time import unittest import unittest.mock as mock from contextlib import redirect_stdout @@ -400,6 +401,12 @@ def setUp(self): if getattr(self, 'assertRaisesRegex', None) is None: self.assertRaisesRegex = self.assertRaisesRegexp + def test_get_time_ms(self): + t0=oapackage.get_time_ms() + time.sleep(.2) + dt=oapackage.get_time_ms(t0) + np.testing.assert_almost_equal(dt, 0.2, decimal=1) + def test_root_form(self): array = oapackage.exampleArray(1, 0) self.assertTrue(oapackage.is_root_form(array, 2)) diff --git a/tests/test_cpplibrary_gwlp.py b/tests/test_cpplibrary_gwlp.py index 55ae9447..4cfb2224 100644 --- a/tests/test_cpplibrary_gwlp.py +++ b/tests/test_cpplibrary_gwlp.py @@ -24,6 +24,12 @@ def test_gwlp_mixed2(self): al = oapackage.exampleArray(4, 0) self.assertEqual(al.GWLP(), (1.0, 0.0, 0.0, 3.5, 2.5, 0.5, 0.5, 0.0)) + def test_regression_gwlp_non_mixed(self): + al = oapackage.exampleArray(57, 0) + gwlp = al.GWLP() + self.assertEqual(gwlp, (1.0, 0.0, 0.0, 13.0, 13.5, 9.0, 4.0)) + self.assertEqual(oapackage.GWLPmixed(al), (1.0, 0.0, 0.0, 13.0, 13.5, 9.0, 4.0)) + if __name__ == '__main__': """ Test code """ diff --git a/utils/oa_depth_extend.cpp b/utils/oa_depth_extend.cpp index 9cb6fee8..f30c6b78 100644 --- a/utils/oa_depth_extend.cpp +++ b/utils/oa_depth_extend.cpp @@ -72,7 +72,7 @@ int main (int argc, char *argv[]) { opt.addUsage (" -Q [time] Logging time (in seconds)"); opt.addUsage ( printfstring ( - " --j5structure [VALUE] Integer, can be J5_ORIGINAL (%d) or J5_45 (%d)", + " -j [VALUE] --j5structure [VALUE] Integer, can be J5_ORIGINAL (%d) or J5_45 (%d)", J5_ORIGINAL, J5_45) .c_str ()); std::string ss = printfstring (" -m [MODE] Algorithm (") + algorithm_t_list () + ")\n"; @@ -178,7 +178,7 @@ int main (int argc, char *argv[]) { oaextendx.checkarrays = 0; oaextendx.use_row_symmetry = 0; - oaextendx.extendarraymode = OAextend::APPENDFULL; + oaextendx.extendarraymode = OAextend::extendarray_mode_t::APPENDFULL; oaextendx.init_column_previous = INITCOLUMN_J5; oaextendx.nLMC = 500000; oaextendx.info (); @@ -284,7 +284,7 @@ int main (int argc, char *argv[]) { delete arraylist; clear_LMCreduction_pool (); - + logstream (SYSTEM) << "#time end: " << currenttime () << std::endl; if (verbose) std::cout << "#time total: " << printfstring ("%.1f", get_time_ms () - time0) << " [s]" << std::endl; diff --git a/utils/oa_select_maxj.cpp b/utils/oa_select_maxj.cpp new file mode 100644 index 00000000..3c695637 --- /dev/null +++ b/utils/oa_select_maxj.cpp @@ -0,0 +1,171 @@ +/** \file oa_select_maxj.cpp + + C++ program: oa_select_maxj + + Select designs from a file with maxj value and write to outputfile + + Author: Pieter Eendebak + Copyright: See LICENSE.txt file that comes with this distribution +*/ + +#include +#include +#include + +#include "anyoption.h" +#include "arraytools.h" +#include "tools.h" + +bool strings_ends_with(const string& a, const string& b) { + if (b.size() > a.size()) return false; + return std::equal(a.begin() + a.size() - b.size(), a.end(), b.begin()); +} + +bool replace(std::string& str, const std::string& from, const std::string& to) { + size_t start_pos = str.find(from); + if(start_pos == std::string::npos) + return false; + str.replace(start_pos, from.length(), to); + return true; +} + +void write_map(const char *filename, const std::map &m) { + map::const_iterator it; + + FILE *fid = fopen(filename,"wt"); + fprintf(fid, "{"); + for (it = m.begin(); it != m.end(); it++) { + fprintf(fid, " \"%ld\": %ld", it->first, it->second); + if (std::next(it)!=m.end()) + fprintf(fid, ","); + } + fprintf(fid, "}\n"); + fclose(fid); +} + +int main (int argc, char *argv[]) { + AnyOption opt; + + /* parse command line options */ + opt.setFlag ("help", 'h'); /* a flag (takes no argument), supporting long and short form */ + opt.setOption ("input", 'i'); + opt.setOption ("output", 'o'); + opt.setOption ("verbose", 'v'); + opt.setOption ("format", 'f'); + + opt.addUsage ("oa_select_maxj: Select designs with max(J5)"); + opt.addUsage ("Usage: oa_select_maxj -i [FILE] -o file [OPTIONS]"); + opt.addUsage (""); + opt.addUsage (" -h --help Prints this help "); + opt.addUsage (" -i [FILE] --input [FILE] Input file "); + opt.addUsage (" -o [STR] --output [FILE] Output file "); + opt.addUsage (" -f [FORMAT] Output format (default: TEXT, or BINARY; B, or DIFF; D, or DIFFZERO; Z) "); + opt.addUsage (" -v [INTEGER] Verbose (default: 2) "); + opt.processCommandArgs (argc, argv); + + print_copyright (); + + /* parse options */ + if (opt.getFlag ("help") || opt.getFlag ('h')) { + opt.printUsage (); + exit (0); + } + + if (opt.getValue ("input") == NULL && opt.getValue ('i') == 0) { + printf ("No input file specified, use oa_select_maxj -h to get help\n"); + exit (1); + } + + const char *inputfile; + inputfile = opt.getValue ('i'); + + const char *outputfile = opt.getValue ('o'); + + if (!strings_ends_with(outputfile, ".oa")) { + throw_runtime_exception("outputfile should end with .oa"); + } + + std::string outputjson = std::string(outputfile); + replace(outputjson, std::string(".oa"), std::string(".json")); + std::map maxj_values; + + int verbose = opt.getIntValue ("verbose", 2); + + std::string format = opt.getStringValue ('f', "BINARY"); + arrayfile::arrayfilemode_t mode = arrayfile_t::parseModeString (format); + + int nb = 8; + if (mode == ABINARY_DIFFZERO) { + nb = 1; + } + /* open the file containing the arrays */ + arrayfile_t *afile = new arrayfile_t (inputfile); + const int rows = afile->nrows; + const int cols = afile->ncols; + const int narrays = afile->narrays; + + if (!afile->isopen ()) { + printf ("oa_select_maxj: problem opening file %s\n", inputfile); + exit (1); + } + if (verbose) + printf ("oa_select_maxj: input %s, output %s\n", inputfile, outputfile); + + /* prepare the output files */ + arrayfile_t *outfid = 0; + + outfid = new arrayfile_t(outputfile, rows, cols, -1, mode, nb); + + jstruct_t *js = 0; + array_link al (rows, cols, 0); + + array_link active_al5; + int active_value = -1; + int active_maxj = - 1; + for (int i = 0; i < narrays; i++) { + if (i % 500000 == 0 && verbose) + printf (" oa_select_maxj: scanning array %d/%d\n", i, narrays); + + + afile->read_array (al); + + array_link al5 = al.selectFirstColumns(5); + + if (js==0) { + // initial initialization + js = new jstruct_t(al5, 5); + std::vector v= js->Fval(); + for(std::vector::iterator it = v.begin(); it != v.end(); ++it) { + maxj_values[(*it)]= 0; + } + } + + if (al5==active_al5) { + // nothing to do + } + else { + js->calcj5(al5); + + active_al5= al5; + active_maxj = js->maxJ (); + active_value = js->calculateF()[0]; + } + + maxj_values[active_maxj]++; + + if (active_value) + outfid->append_array (al); + } + + write_map(outputjson.c_str(), maxj_values); + + outfid->finisharrayfile (); + + /* clean up memory */ + afile->closefile (); + outfid->closefile(); + delete afile; + delete outfid; + + return 0; +} diff --git a/utils/oaextend.cpp b/utils/oaextend.cpp index 522a877d..926f7df3 100644 --- a/utils/oaextend.cpp +++ b/utils/oaextend.cpp @@ -1,18 +1,14 @@ /** \file oaextend.cpp - C++ program: oaextendmpi / oaextendsingle + C++ program: oaextendsingle - oaextendmpi and oaextendsingle can calculate series of non-isomorphic orthogonal arrays + oaextendsingle can calculate series of non-isomorphic orthogonal arrays Author: Pieter Eendebak , (C) 2014 Copyright: See COPYING file that comes with this distribution */ -#ifdef OAEXTEND_MULTICORE -// this must be included before stdio.h -#include -#endif #include #include @@ -34,18 +30,12 @@ #include "strength.h" #include "tools.h" -#ifdef OAEXTEND_MULTICORE - -#include "mpi.h" -#include "mpitools.h" -#else /** number of master processor */ const int MASTER = 0; /* OAEXTEND_SINGLECORE */ int extend_slave_code (const int this_rank, OAextend const &oaextend) { return 0; } -#endif /*! Save_arrays writes all the arrays from solutions to a file. The filename is obtained from the number of factors and the number of columns so far. Then the file header contains the number of columns in the @@ -57,7 +47,7 @@ the number of runs and the number of arrays in the file. \param resultprefix Prefix string for filename \param mode Mode for output file */ -int save_arrays(arraylist_t &solutions, const arraydata_t *ad, +int save_arrays(arraylist_t &solutions, const arraydata_t *ad, const char *resultprefix, arrayfile::arrayfilemode_t mode) { string fname = resultprefix; @@ -103,11 +93,7 @@ AnyOption *parseOptions (int argc, char *argv[], algorithm_t &algorithm) { opt->setOption ("initcolprev", 'I'); /* algorithm method */ opt->addUsage ("Orthonal Arrays: extend orthogonal arrays"); -#ifdef OAEXTEND_SINGLECORE - opt->addUsage ("Usage: oaextendmpi [OPTIONS]"); -#else opt->addUsage ("Usage: oaextendsingle [OPTIONS]"); -#endif opt->addUsage (""); opt->addUsage (" -h --help Prints this help "); @@ -168,41 +154,22 @@ int init_restart (const char *fname, colindex_t &cols, arraylist_t &solutions) { } /** - * @brief Main function for oaextendmpi and oaextendsingle + * @brief Main function for oaextendsingle * @param argc * @param argv[] * @return */ int main (int argc, char *argv[]) { -#ifdef OAEXTEND_SINGLECORE const int n_processors = 1; const int this_rank = 0; -#else - MPI::Init (argc, argv); - - const int n_processors = MPI::COMM_WORLD.Get_size (); - const int this_rank = MPI::COMM_WORLD.Get_rank (); -#endif - - // printf("MPI: this_rank %d\n", this_rank); /*----------SET STARTING ARRAY(S)-----------*/ if (this_rank != MASTER) { -#ifdef OAEXTEND_MULTICORE - slave_print (QUIET, "M: running core %d/%d\n", this_rank, n_processors); -#endif algorithm_t algorithm = MODE_INVALID; AnyOption *opt = parseOptions (argc, argv, algorithm); OAextend oaextend; oaextend.setAlgorithm (algorithm); -#ifdef OAEXTEND_MULTICORE - log_print (NORMAL, "slave: receiving algorithm int\n"); - algorithm = (algorithm_t)receive_int_slave (); - oaextend.setAlgorithm (algorithm); -// printf("slave %d: receive algorithm %d\n", this_rank, algorithm); -#endif - extend_slave_code (this_rank, oaextend); delete opt; @@ -241,7 +208,7 @@ int main (int argc, char *argv[]) { if (streaming) { logstream (SYSTEM) << "operating in streaming mode, sorting of arrays will not work " << std::endl; - oaextend.extendarraymode = OAextend::STOREARRAY; + oaextend.extendarraymode = OAextend::extendarray_mode_t::STOREARRAY; } // J5_45 @@ -284,23 +251,11 @@ int main (int argc, char *argv[]) { oaextend.setAlgorithmAuto (ad); } -#ifdef OAEXTEND_SINGLECORE if (initcolprev == 0) { log_print (NORMAL, "setting oaextend.init_column_previous to %d\n", INITCOLUMN_ZERO); oaextend.init_column_previous = INITCOLUMN_ZERO; } -#endif -#ifdef OAEXTEND_MULTICORE - int alg = oaextend.getAlgorithm (); - for (int i = 0; i < n_processors; i++) { - if (i == MASTER) { - continue; - } - log_print (NORMAL, "MASTER: sending algorithm %d to slave %d\n", alg, i); - send_int (alg, i); - } -#endif colindex_t col_start; log_print (SYSTEM, "using algorithm %d (%s)\n", oaextend.getAlgorithm (), @@ -315,9 +270,6 @@ int main (int argc, char *argv[]) { logstream (SYSTEM) << "Problem with restart from " << opt->getValue ('r') << "" << endl; -#ifdef OAEXTEND_MPI - MPI_Abort (MPI_COMM_WORLD, 1); -#endif exit (1); } @@ -340,12 +292,9 @@ int main (int argc, char *argv[]) { } else { // starting with root - if (check_divisibility (ad) == false) { + if (check_divisibility (ad) == false) { log_print (SYSTEM, "ERROR: Failed divisibility test!\n"); -#ifdef OAEXTEND_MPI - MPI_Abort (MPI_COMM_WORLD, 1); -#endif exit (1); } create_root (ad, solutions); @@ -376,7 +325,7 @@ int main (int argc, char *argv[]) { } log_print (SYSTEM, "Starting with column %d (%d, total time: %.2f [s])\n", - current_col + 1, (int)solutions.size (), get_time_ms () - Tstart); + current_col + 1, (int)solutions.size (), get_time_ms (Tstart)); nr_extensions = 0; arraylist_t::const_iterator cur_extension; @@ -390,29 +339,11 @@ int main (int argc, char *argv[]) { nr_extensions += extend_array (*cur_extension, adcol, current_col, extensions, oaextend); } else { -#ifdef OAEXTEND_MPI - throw_runtime_exception("mpi version of oextend is no longer supported"); - double Ttmp = get_time_ms (); - int slave = collect_solutions_single (extensions, adcol); - log_print (DEBUG, " time: %.2f, collect time %.2f\n", - (get_time_ms () - Tstart), (get_time_ms () - Ttmp)); - - /* OPTIMIZE: send multiple arrays to slave 1 once step */ - extend_array_mpi (cur_extension->array, 1, adcol, current_col, slave); -#endif } -#ifdef OAEXTEND_MPI - if ((csol + 1) % nsolsync == 0) { - collect_solutions_all (extensions, adcol); - } -#endif // OPTIMIZE: periodically write part of the solutions to disk csol++; /* increase current solution */ } -#ifdef OAEXTEND_MPI - collect_solutions_all (extensions, adcol); -#endif if (checkloglevel (NORMAL)) { csol = 0; @@ -438,10 +369,6 @@ int main (int argc, char *argv[]) { delete adcol; -#ifdef OAANALYZE_DISCR - logstream (NORMAL) << "Discriminant: " << endl; - print_discriminant (ad->N, current_col + 1); -#endif } /*------------------------*/ time (&seconds); @@ -454,19 +381,11 @@ int main (int argc, char *argv[]) { solutions.clear (); delete ad; - -#ifdef OAEXTEND_MPI - stop_slaves (); -#endif } delete opt; } /* end of master code */ -#ifdef OAEXTEND_MPI - MPI_Finalize (); -#endif - return 0; } diff --git a/utils/oasplit.cpp b/utils/oasplit.cpp index de43bdba..6a7384da 100644 --- a/utils/oasplit.cpp +++ b/utils/oasplit.cpp @@ -82,12 +82,11 @@ int main (int argc, char *argv[]) { if (mode == ABINARY_DIFFZERO) { nb = 1; } - /* open the file containint the arrays */ - int rows, cols, narrays; + /* open the file containing the arrays */ arrayfile_t *afile = new arrayfile_t (inputfile); - rows = afile->nrows; - cols = afile->ncols; - narrays = afile->narrays; + const int rows = afile->nrows; + const int cols = afile->ncols; + const int narrays = afile->narrays; if (!afile->isopen ()) { printf ("oasplit: problem opening file %s\n", inputfile); diff --git a/utils/oastreaming.cpp b/utils/oastreaming.cpp index 135c8d7b..d468055a 100644 --- a/utils/oastreaming.cpp +++ b/utils/oastreaming.cpp @@ -142,7 +142,7 @@ int main (int argc, char *argv[]) { OAextend oaextend; oaextend.setAlgorithm (algorithm); - oaextend.extendarraymode = OAextend::STOREARRAY; + oaextend.extendarraymode = OAextend::extendarray_mode_t::STOREARRAY; // J5_45 int xx = opt->getFlag ('x'); diff --git a/utils/oatest.cpp b/utils/oatest.cpp index 619e5b39..03bfeca6 100644 --- a/utils/oatest.cpp +++ b/utils/oatest.cpp @@ -121,6 +121,23 @@ int main(int argc, char* argv[]) { if (input == 0) input = "test.oa"; + int run_size = 100; +int strength = 2; +int number_of_factors = 4; + int factor_levels = 8; + arraydata_t adata(8, run_size, strength, number_of_factors); + adata.show(); + + arraylist_t ll; + printf("create root:\n"); + ll.push_back(adata.create_root()); + + extend_arraylist(ll, adata); + + printf("done\n"); + exit(0); +//array_list_3columns=oapackage.extend_arraylist(array_list, arrayclass) + array_link G(4, 4, 0); G.at(0, 0) = 1; G.at(1, 0) = 1; @@ -140,9 +157,9 @@ int main(int argc, char* argv[]) { myprintf("perm: "); print_perm(perm); return 0; - printf("test! %ld\n", choose(6,4)); + printf("test! %ld\n", (long)choose(6,4)); for (int i = 4; i < 12; i++) - printf("choose(%d, %d): %ld\n", i, i-2, choose(i, i - 2) - ncombs(i, i - 2)); + printf("choose(%d, %d): %ld\n", i, i-2, (long)(choose(i, i - 2) - ncombs(i, i - 2)) ); fflush(0); array_link A = exampleArray(56, 1); diff --git a/utils/oaunittest.cpp b/utils/oaunittest.cpp index bd55d29f..2c62bc00 100644 --- a/utils/oaunittest.cpp +++ b/utils/oaunittest.cpp @@ -114,18 +114,15 @@ int oaunittest (int verbose, int writetests = 0, int randval = 0) { } { - cprintf (verbose, "%s: conference matrix Fvalues\n", bstr); + cprintf(verbose, "%s: conference matrix Fvalues\n", bstr); + std::vector< int > jj5 = Jcharacteristics_conference(exampleArray(18, 0), 5); + myassert(jj5 == std::vector{ -1, -1, -1, -5, 11, -3, 3, 3, -1, -1, 1, -1, 3, -1, 1, 7, 3, 5, -1, 1, 1}, "Jcharacteristics_conference result incorrect"); + array_link al = exampleArray (22, 0); if (verbose >= 2) al.show (); - if (0) { - std::vector< int > f3 = al.FvaluesConference (3); - if (verbose >= 2) { - printf ("F3: "); - display_vector (f3); - printf ("\n"); - } - } + std::vector< int > jj3 = Jcharacteristics_conference(al, 3); + myassert(jj3 == std::vector { 0, 0, 0, 0 }, "Jcharacteristics_conference result incorrect"); const int N = al.n_rows; jstructconference_t js (N, 4); @@ -203,16 +200,6 @@ int oaunittest (int verbose, int writetests = 0, int randval = 0) { conference_transformation_t Ti = T.inverse (); array_link alx = Ti.apply (T.apply (al)); - if (0) { - printf ("input array:\n"); - al.showarray (); - T.show (); - printf ("transformed array:\n"); - T.apply (al).showarray (); - Ti.show (); - alx.showarray (); - } - myassert (alx == al, "transformation of conference matrix\n"); } @@ -567,12 +554,7 @@ int oaunittest (int verbose, int writetests = 0, int randval = 0) { array_link al2 = reduction.transformation->apply (al); array_link alr = al2.reduceLMC (); - if (0) { - printf ("\n reduction complete:\n"); - al2.showarray (); - printf (" --->\n"); - alr.showarray (); - } + bool c = (al == alr); if (!c) { printf ("oaunittest: error: reduction of randomized array failed!\n"); @@ -647,7 +629,7 @@ int oaunittest (int verbose, int writetests = 0, int randval = 0) { exit (1); } - for (int ii = 11; ii < 11; ii++) { + for (int ii = 11; ii < 12; ii++) { printf ("ii %d: ", ii); al = exampleArray (ii, vb); al.showarray (); @@ -687,7 +669,7 @@ int oaunittest (int verbose, int writetests = 0, int randval = 0) { boost::filesystem::path tmpdir = boost::filesystem::temp_directory_path (); boost::filesystem::path temp = boost::filesystem::unique_path ("test-%%%%%%%.oa"); - const std::string tempstr = (tmpdir / temp).native (); + const std::string tempstr = (tmpdir / temp).native (); if (verbose >= 2) printf ("generate text OA file: %s\n", tempstr.c_str ()); @@ -697,7 +679,7 @@ int oaunittest (int verbose, int writetests = 0, int randval = 0) { int narrays = 10; arrayfile_t afile (tempstr.c_str (), nrows, ncols, narrays, ATEXT); for (int i = 0; i < narrays; i++) { - array_link al (nrows, ncols, array_link::INDEX_DEFAULT); + array_link al (nrows, ncols, array_link::INDEX_DEFAULT); afile.append_array (al); } afile.closefile ();