diff --git a/.github/workflows/deploy_release.yml b/.github/workflows/deploy_release.yml index b18c31ff2a..873c5f6e05 100644 --- a/.github/workflows/deploy_release.yml +++ b/.github/workflows/deploy_release.yml @@ -36,11 +36,6 @@ jobs: - run: echo "AMICI_DIR=$(pwd)" >> $GITHUB_ENV - - name: Remove direct dependencies from setup.cfg - # Remove any "git+https"-based dependencies that are not supported - # by PyPI and are only required for testing. - run: sed -i '/git+https/d' python/sdist/pyproject.toml - - name: sdist run: scripts/buildSdist.sh diff --git a/.github/workflows/test_benchmark_collection_models.yml b/.github/workflows/test_benchmark_collection_models.yml index 7b1568b1fc..0b14d73716 100644 --- a/.github/workflows/test_benchmark_collection_models.yml +++ b/.github/workflows/test_benchmark_collection_models.yml @@ -68,7 +68,7 @@ jobs: - name: Run Gradient Checks run: | pip install git+https://github.com/ICB-DCM/fiddy.git \ - && cd tests/benchmark-models && pytest ./test_petab_benchmark.py + && cd tests/benchmark-models && pytest --durations=10 ./test_petab_benchmark.py # upload results - uses: actions/upload-artifact@v4 diff --git a/.github/workflows/test_petab_test_suite.yml b/.github/workflows/test_petab_test_suite.yml index 64277fe267..c16b40bc30 100644 --- a/.github/workflows/test_petab_test_suite.yml +++ b/.github/workflows/test_petab_test_suite.yml @@ -15,7 +15,7 @@ jobs: build: name: PEtab Testsuite - runs-on: ubuntu-22.04 + runs-on: ubuntu-latest env: ENABLE_GCOV_COVERAGE: TRUE @@ -68,6 +68,12 @@ jobs: && source ./venv/bin/activate \ && cd petab_test_suite && pip3 install -e . + - name: Install PEtab benchmark collection + run: | + git clone --depth 1 https://github.com/benchmarking-initiative/Benchmark-Models-PEtab.git \ + && export BENCHMARK_COLLECTION="$(pwd)/Benchmark-Models-PEtab/Benchmark-Models/" \ + && source venv/bin/activate && python -m pip install -e $BENCHMARK_COLLECTION/../src/python + - name: Install petab run: | source ./venv/bin/activate \ diff --git a/.github/workflows/test_python_cplusplus.yml b/.github/workflows/test_python_cplusplus.yml index 5ca8076a69..a74f407d00 100644 --- a/.github/workflows/test_python_cplusplus.yml +++ b/.github/workflows/test_python_cplusplus.yml @@ -338,6 +338,12 @@ jobs: path: ${{ env.pooch-cache }} key: ${{ runner.os }}-py${{ matrix.python-version }}-${{ github.job }} + - name: Install PEtab benchmark collection + run: | + git clone --depth 1 https://github.com/benchmarking-initiative/Benchmark-Models-PEtab.git \ + && export BENCHMARK_COLLECTION="$(pwd)/Benchmark-Models-PEtab/Benchmark-Models/" \ + && source venv/bin/activate && python -m pip install -e $BENCHMARK_COLLECTION/../src/python + - name: Python tests run: | scripts/run-python-tests.sh \ diff --git a/.gitignore b/.gitignore index 453f29a50f..0faae713ab 100644 --- a/.gitignore +++ b/.gitignore @@ -204,3 +204,4 @@ cache_fiddy/* debug/* tests/benchmark-models/cache_fiddy/* venv/* +.coverage diff --git a/CHANGELOG.md b/CHANGELOG.md index 28d5a16f3a..44e57b2f22 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -4,6 +4,29 @@ See also our [versioning policy](https://amici.readthedocs.io/en/latest/versioni ## v0.X Series + +### v0.26.3 (2024-10-03) + +**Fixes** + +* Skip building SuiteSparse shared libraries and build all subprojects together + for slightly faster package installation + + by @dweindl in https://github.com/AMICI-dev/AMICI/pull/2514 + and https://github.com/AMICI-dev/AMICI/pull/2519 + +* Got rid of petab `DeprecationWarnings` when using the `amici_import_petab` + CLI + + by @dweindl in https://github.com/AMICI-dev/AMICI/pull/2517 + +* Now also sundials and suitesparse are built in debug mode when installing + with `ENABLE_AMICI_DEBUGGING=TRUE` + + by @dweindl in https://github.com/AMICI-dev/AMICI/pull/2515 + +**Full Changelog**: https://github.com/AMICI-dev/AMICI/compare/v0.26.2...v0.26.3 + ### v0.26.2 (2024-09-25) **Fixes** diff --git a/CMakeLists.txt b/CMakeLists.txt index 12ecdc3fcd..bdb9d45fb3 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1,10 +1,10 @@ # # Build AMICI library # -cmake_minimum_required(VERSION 3.15) +cmake_minimum_required(VERSION 3.22) # When updating the policy version, please also update it in # src/CMakeLists.template.cmake -cmake_policy(VERSION 3.15...3.30) +cmake_policy(VERSION 3.22...3.30) project(amici) @@ -130,7 +130,8 @@ elseif(AMICI_TRY_ENABLE_HDF5) endif() set(VENDORED_SUNDIALS_DIR ${CMAKE_CURRENT_SOURCE_DIR}/ThirdParty/sundials) -set(SUNDIALS_PRIVATE_INCLUDE_DIRS "${VENDORED_SUNDIALS_DIR}/src") +set(SUNDIALS_PRIVATE_INCLUDE_DIRS + "${VENDORED_SUNDIALS_DIR}/src;${VENDORED_SUNDIALS_DIR}/src/sundials") # Handle different sundials build/install dirs, depending on whether we are # building the Python extension only or the full C++ interface if(AMICI_PYTHON_BUILD_EXT_ONLY) @@ -224,7 +225,14 @@ set(AMICI_CXX_OPTIONS "" CACHE STRING "C++ options for libamici (semicolon-separated)") target_compile_options(${PROJECT_NAME} PRIVATE "${AMICI_CXX_OPTIONS}") - +set_property( + SOURCE src/solver_cvodes.cpp + APPEND + PROPERTY INCLUDE_DIRECTORIES "${SUNDIALS_PRIVATE_INCLUDE_DIRS}") +set_property( + SOURCE src/solver_idas.cpp + APPEND + PROPERTY INCLUDE_DIRECTORIES "${SUNDIALS_PRIVATE_INCLUDE_DIRS}") add_dependencies(${PROJECT_NAME} version) file(GLOB PUBLIC_HEADERS include/amici/*.h) @@ -237,8 +245,7 @@ target_include_directories( $ $ $ - $ - PRIVATE ${SUNDIALS_PRIVATE_INCLUDE_DIRS}) + $) if("$ENV{ENABLE_AMICI_DEBUGGING}" OR "$ENV{ENABLE_GCOV_COVERAGE}" diff --git a/ThirdParty/SuiteSparse/AMD/Source/amd.f b/ThirdParty/SuiteSparse/AMD/Source/amd.f deleted file mode 100644 index 0378073648..0000000000 --- a/ThirdParty/SuiteSparse/AMD/Source/amd.f +++ /dev/null @@ -1,1220 +0,0 @@ -C----------------------------------------------------------------------- -C AMD/Source/amd.f: Fortran version of AMD -C----------------------------------------------------------------------- - -C AMD, Copyright (c) 1996-2022, Timothy A. Davis, Patrick R. Amestoy, -C and C Iain S. Duff. All Rights Reserved. -C SPDX-License-Identifier: BSD-3-clause - -C----------------------------------------------------------------------- - - SUBROUTINE AMD - $ (N, PE, IW, LEN, IWLEN, PFREE, NV, NEXT, - $ LAST, HEAD, ELEN, DEGREE, NCMPA, W) - - INTEGER N, IWLEN, PFREE, NCMPA, IW (IWLEN), PE (N), - $ DEGREE (N), NV (N), NEXT (N), LAST (N), HEAD (N), - $ ELEN (N), W (N), LEN (N) - -C Given a representation of the nonzero pattern of a symmetric matrix, -C A, (excluding the diagonal) perform an approximate minimum -C (UMFPACK/MA38-style) degree ordering to compute a pivot order -C such that the introduction of nonzeros (fill-in) in the Cholesky -C factors A = LL^T are kept low. At each step, the pivot -C selected is the one with the minimum UMFPACK/MA38-style -C upper-bound on the external degree. -C -C Aggresive absorption is used to tighten the bound on the degree. - -C ********************************************************************** -C ***** CAUTION: ARGUMENTS ARE NOT CHECKED FOR ERRORS ON INPUT. ****** -C ********************************************************************** - -C References: -C -C [1] Timothy A. Davis and Iain Duff, "An unsymmetric-pattern -C multifrontal method for sparse LU factorization", SIAM J. -C Matrix Analysis and Applications, vol. 18, no. 1, pp. -C 140-158. Discusses UMFPACK / MA38, which first introduced -C the approximate minimum degree used by this routine. -C -C [2] Patrick Amestoy, Timothy A. Davis, and Iain S. Duff, "An -C approximate degree ordering algorithm," SIAM J. Matrix -C Analysis and Applications, vol. 17, no. 4, pp. 886-905, -C 1996. Discusses AMD, AMDBAR, and MC47B. -C -C [3] Alan George and Joseph Liu, "The evolution of the minimum -C degree ordering algorithm," SIAM Review, vol. 31, no. 1, -C pp. 1-19, 1989. We list below the features mentioned in -C that paper that this code includes: -C -C mass elimination: -C Yes. MA27 relied on supervariable detection for mass -C elimination. -C indistinguishable nodes: -C Yes (we call these "supervariables"). This was also in -C the MA27 code - although we modified the method of -C detecting them (the previous hash was the true degree, -C which we no longer keep track of). A supervariable is -C a set of rows with identical nonzero pattern. All -C variables in a supervariable are eliminated together. -C Each supervariable has as its numerical name that of -C one of its variables (its principal variable). -C quotient graph representation: -C Yes. We use the term "element" for the cliques formed -C during elimination. This was also in the MA27 code. -C The algorithm can operate in place, but it will work -C more efficiently if given some "elbow room." -C element absorption: -C Yes. This was also in the MA27 code. -C external degree: -C Yes. The MA27 code was based on the true degree. -C incomplete degree update and multiple elimination: -C No. This was not in MA27, either. Our method of -C degree update within MC47B/BD is element-based, not -C variable-based. It is thus not well-suited for use -C with incomplete degree update or multiple elimination. - -C----------------------------------------------------------------------- -C Authors, and Copyright (C) 1995 by: -C Timothy A. Davis, Patrick Amestoy, Iain S. Duff, & John K. Reid. -C -C Acknowledgements: -C This work (and the UMFPACK package) was supported by the -C National Science Foundation (ASC-9111263 and DMS-9223088). -C The UMFPACK/MA38 approximate degree update algorithm, the -C unsymmetric analog which forms the basis of MC47B/BD, was -C developed while Tim Davis was supported by CERFACS (Toulouse, -C France) in a post-doctoral position. -C -C Date: September, 1995 -C----------------------------------------------------------------------- - -C----------------------------------------------------------------------- -C INPUT ARGUMENTS (unaltered): -C----------------------------------------------------------------------- - -C n: The matrix order. -C -C Restriction: 1 .le. n .lt. (iovflo/2)-2, where iovflo is -C the largest positive integer that your computer can represent. - -C iwlen: The length of iw (1..iwlen). On input, the matrix is -C stored in iw (1..pfree-1). However, iw (1..iwlen) should be -C slightly larger than what is required to hold the matrix, at -C least iwlen .ge. pfree + n is recommended. Otherwise, -C excessive compressions will take place. -C *** We do not recommend running this algorithm with *** -C *** iwlen .lt. pfree + n. *** -C *** Better performance will be obtained if *** -C *** iwlen .ge. pfree + n *** -C *** or better yet *** -C *** iwlen .gt. 1.2 * pfree *** -C *** (where pfree is its value on input). *** -C The algorithm will not run at all if iwlen .lt. pfree-1. -C -C Restriction: iwlen .ge. pfree-1 - -C----------------------------------------------------------------------- -C INPUT/OUPUT ARGUMENTS: -C----------------------------------------------------------------------- - -C pe: On input, pe (i) is the index in iw of the start of row i, or -C zero if row i has no off-diagonal non-zeros. -C -C During execution, it is used for both supervariables and -C elements: -C -C * Principal supervariable i: index into iw of the -C description of supervariable i. A supervariable -C represents one or more rows of the matrix -C with identical nonzero pattern. -C * Non-principal supervariable i: if i has been absorbed -C into another supervariable j, then pe (i) = -j. -C That is, j has the same pattern as i. -C Note that j might later be absorbed into another -C supervariable j2, in which case pe (i) is still -j, -C and pe (j) = -j2. -C * Unabsorbed element e: the index into iw of the description -C of element e, if e has not yet been absorbed by a -C subsequent element. Element e is created when -C the supervariable of the same name is selected as -C the pivot. -C * Absorbed element e: if element e is absorbed into element -C e2, then pe (e) = -e2. This occurs when the pattern of -C e (that is, Le) is found to be a subset of the pattern -C of e2 (that is, Le2). If element e is "null" (it has -C no nonzeros outside its pivot block), then pe (e) = 0. -C -C On output, pe holds the assembly tree/forest, which implicitly -C represents a pivot order with identical fill-in as the actual -C order (via a depth-first search of the tree). -C -C On output: -C If nv (i) .gt. 0, then i represents a node in the assembly tree, -C and the parent of i is -pe (i), or zero if i is a root. -C If nv (i) = 0, then (i,-pe (i)) represents an edge in a -C subtree, the root of which is a node in the assembly tree. - -C pfree: On input the tail end of the array, iw (pfree..iwlen), -C is empty, and the matrix is stored in iw (1..pfree-1). -C During execution, additional data is placed in iw, and pfree -C is modified so that iw (pfree..iwlen) is always the unused part -C of iw. On output, pfree is set equal to the size of iw that -C would have been needed for no compressions to occur. If -C ncmpa is zero, then pfree (on output) is less than or equal to -C iwlen, and the space iw (pfree+1 ... iwlen) was not used. -C Otherwise, pfree (on output) is greater than iwlen, and all the -C memory in iw was used. - -C----------------------------------------------------------------------- -C INPUT/MODIFIED (undefined on output): -C----------------------------------------------------------------------- - -C len: On input, len (i) holds the number of entries in row i of the -C matrix, excluding the diagonal. The contents of len (1..n) -C are undefined on output. - -C iw: On input, iw (1..pfree-1) holds the description of each row i -C in the matrix. The matrix must be symmetric, and both upper -C and lower triangular parts must be present. The diagonal must -C not be present. Row i is held as follows: -C -C len (i): the length of the row i data structure -C iw (pe (i) ... pe (i) + len (i) - 1): -C the list of column indices for nonzeros -C in row i (simple supervariables), excluding -C the diagonal. All supervariables start with -C one row/column each (supervariable i is just -C row i). -C if len (i) is zero on input, then pe (i) is ignored -C on input. -C -C Note that the rows need not be in any particular order, -C and there may be empty space between the rows. -C -C During execution, the supervariable i experiences fill-in. -C This is represented by placing in i a list of the elements -C that cause fill-in in supervariable i: -C -C len (i): the length of supervariable i -C iw (pe (i) ... pe (i) + elen (i) - 1): -C the list of elements that contain i. This list -C is kept short by removing absorbed elements. -C iw (pe (i) + elen (i) ... pe (i) + len (i) - 1): -C the list of supervariables in i. This list -C is kept short by removing nonprincipal -C variables, and any entry j that is also -C contained in at least one of the elements -C (j in Le) in the list for i (e in row i). -C -C When supervariable i is selected as pivot, we create an -C element e of the same name (e=i): -C -C len (e): the length of element e -C iw (pe (e) ... pe (e) + len (e) - 1): -C the list of supervariables in element e. -C -C An element represents the fill-in that occurs when supervariable -C i is selected as pivot (which represents the selection of row i -C and all non-principal variables whose principal variable is i). -C We use the term Le to denote the set of all supervariables -C in element e. Absorbed supervariables and elements are pruned -C from these lists when computationally convenient. -C -C CAUTION: THE INPUT MATRIX IS OVERWRITTEN DURING COMPUTATION. -C The contents of iw are undefined on output. - -C----------------------------------------------------------------------- -C OUTPUT (need not be set on input): -C----------------------------------------------------------------------- - -C nv: During execution, abs (nv (i)) is equal to the number of rows -C that are represented by the principal supervariable i. If i is -C a nonprincipal variable, then nv (i) = 0. Initially, -C nv (i) = 1 for all i. nv (i) .lt. 0 signifies that i is a -C principal variable in the pattern Lme of the current pivot -C element me. On output, nv (e) holds the true degree of element -C e at the time it was created (including the diagonal part). - -C ncmpa: The number of times iw was compressed. If this is -C excessive, then the execution took longer than what could have -C been. To reduce ncmpa, try increasing iwlen to be 10% or 20% -C larger than the value of pfree on input (or at least -C iwlen .ge. pfree + n). The fastest performance will be -C obtained when ncmpa is returned as zero. If iwlen is set to -C the value returned by pfree on *output*, then no compressions -C will occur. - -C elen: See the description of iw above. At the start of execution, -C elen (i) is set to zero. During execution, elen (i) is the -C number of elements in the list for supervariable i. When e -C becomes an element, elen (e) = -nel is set, where nel is the -C current step of factorization. elen (i) = 0 is done when i -C becomes nonprincipal. -C -C For variables, elen (i) .ge. 0 holds until just before the -C permutation vectors are computed. For elements, -C elen (e) .lt. 0 holds. -C -C On output elen (1..n) holds the inverse permutation (the same -C as the 'INVP' argument in Sparspak). That is, if k = elen (i), -C then row i is the kth pivot row. Row i of A appears as the -C (elen(i))-th row in the permuted matrix, PAP^T. - -C last: In a degree list, last (i) is the supervariable preceding i, -C or zero if i is the head of the list. In a hash bucket, -C last (i) is the hash key for i. last (head (hash)) is also -C used as the head of a hash bucket if head (hash) contains a -C degree list (see head, below). -C -C On output, last (1..n) holds the permutation (the same as the -C 'PERM' argument in Sparspak). That is, if i = last (k), then -C row i is the kth pivot row. Row last (k) of A is the k-th row -C in the permuted matrix, PAP^T. - -C----------------------------------------------------------------------- -C LOCAL (not input or output - used only during execution): -C----------------------------------------------------------------------- - -C degree: If i is a supervariable, then degree (i) holds the -C current approximation of the external degree of row i (an upper -C bound). The external degree is the number of nonzeros in row i, -C minus abs (nv (i)) (the diagonal part). The bound is equal to -C the external degree if elen (i) is less than or equal to two. -C -C We also use the term "external degree" for elements e to refer -C to |Le \ Lme|. If e is an element, then degree (e) holds |Le|, -C which is the degree of the off-diagonal part of the element e -C (not including the diagonal part). - -C head: head is used for degree lists. head (deg) is the first -C supervariable in a degree list (all supervariables i in a -C degree list deg have the same approximate degree, namely, -C deg = degree (i)). If the list deg is empty then -C head (deg) = 0. -C -C During supervariable detection head (hash) also serves as a -C pointer to a hash bucket. -C If head (hash) .gt. 0, there is a degree list of degree hash. -C The hash bucket head pointer is last (head (hash)). -C If head (hash) = 0, then the degree list and hash bucket are -C both empty. -C If head (hash) .lt. 0, then the degree list is empty, and -C -head (hash) is the head of the hash bucket. -C After supervariable detection is complete, all hash buckets -C are empty, and the (last (head (hash)) = 0) condition is -C restored for the non-empty degree lists. - -C next: next (i) is the supervariable following i in a link list, or -C zero if i is the last in the list. Used for two kinds of -C lists: degree lists and hash buckets (a supervariable can be -C in only one kind of list at a time). - -C w: The flag array w determines the status of elements and -C variables, and the external degree of elements. -C -C for elements: -C if w (e) = 0, then the element e is absorbed -C if w (e) .ge. wflg, then w (e) - wflg is the size of -C the set |Le \ Lme|, in terms of nonzeros (the -C sum of abs (nv (i)) for each principal variable i that -C is both in the pattern of element e and NOT in the -C pattern of the current pivot element, me). -C if wflg .gt. w (e) .gt. 0, then e is not absorbed and has -C not yet been seen in the scan of the element lists in -C the computation of |Le\Lme| in loop 150 below. -C -C for variables: -C during supervariable detection, if w (j) .ne. wflg then j is -C not in the pattern of variable i -C -C The w array is initialized by setting w (i) = 1 for all i, -C and by setting wflg = 2. It is reinitialized if wflg becomes -C too large (to ensure that wflg+n does not cause integer -C overflow). - -C----------------------------------------------------------------------- -C LOCAL INTEGERS: -C----------------------------------------------------------------------- - - INTEGER DEG, DEGME, DEXT, DMAX, E, ELENME, ELN, HASH, HMOD, I, - $ ILAST, INEXT, J, JLAST, JNEXT, K, KNT1, KNT2, KNT3, - $ LENJ, LN, MAXMEM, ME, MEM, MINDEG, NEL, NEWMEM, - $ NLEFT, NVI, NVJ, NVPIV, SLENME, WE, WFLG, WNVI, X - -C deg: the degree of a variable or element -C degme: size, |Lme|, of the current element, me (= degree (me)) -C dext: external degree, |Le \ Lme|, of some element e -C dmax: largest |Le| seen so far -C e: an element -C elenme: the length, elen (me), of element list of pivotal var. -C eln: the length, elen (...), of an element list -C hash: the computed value of the hash function -C hmod: the hash function is computed modulo hmod = max (1,n-1) -C i: a supervariable -C ilast: the entry in a link list preceding i -C inext: the entry in a link list following i -C j: a supervariable -C jlast: the entry in a link list preceding j -C jnext: the entry in a link list, or path, following j -C k: the pivot order of an element or variable -C knt1: loop counter used during element construction -C knt2: loop counter used during element construction -C knt3: loop counter used during compression -C lenj: len (j) -C ln: length of a supervariable list -C maxmem: amount of memory needed for no compressions -C me: current supervariable being eliminated, and the -C current element created by eliminating that -C supervariable -C mem: memory in use assuming no compressions have occurred -C mindeg: current minimum degree -C nel: number of pivots selected so far -C newmem: amount of new memory needed for current pivot element -C nleft: n - nel, the number of nonpivotal rows/columns remaining -C nvi: the number of variables in a supervariable i (= nv (i)) -C nvj: the number of variables in a supervariable j (= nv (j)) -C nvpiv: number of pivots in current element -C slenme: number of variables in variable list of pivotal variable -C we: w (e) -C wflg: used for flagging the w array. See description of iw. -C wnvi: wflg - nv (i) -C x: either a supervariable or an element - -C----------------------------------------------------------------------- -C LOCAL POINTERS: -C----------------------------------------------------------------------- - - INTEGER P, P1, P2, P3, PDST, PEND, PJ, PME, PME1, PME2, PN, PSRC - -C Any parameter (pe (...) or pfree) or local variable -C starting with "p" (for Pointer) is an index into iw, -C and all indices into iw use variables starting with -C "p." The only exception to this rule is the iwlen -C input argument. - -C p: pointer into lots of things -C p1: pe (i) for some variable i (start of element list) -C p2: pe (i) + elen (i) - 1 for some var. i (end of el. list) -C p3: index of first supervariable in clean list -C pdst: destination pointer, for compression -C pend: end of memory to compress -C pj: pointer into an element or variable -C pme: pointer into the current element (pme1...pme2) -C pme1: the current element, me, is stored in iw (pme1...pme2) -C pme2: the end of the current element -C pn: pointer into a "clean" variable, also used to compress -C psrc: source pointer, for compression - -C----------------------------------------------------------------------- -C FUNCTIONS CALLED: -C----------------------------------------------------------------------- - - INTRINSIC MAX, MIN, MOD - -C======================================================================= -C INITIALIZATIONS -C======================================================================= - - WFLG = 2 - MINDEG = 1 - NCMPA = 0 - NEL = 0 - HMOD = MAX (1, N-1) - DMAX = 0 - MEM = PFREE - 1 - MAXMEM = MEM - ME = 0 - - DO 10 I = 1, N - LAST (I) = 0 - HEAD (I) = 0 - NV (I) = 1 - W (I) = 1 - ELEN (I) = 0 - DEGREE (I) = LEN (I) -10 CONTINUE - -C ---------------------------------------------------------------- -C initialize degree lists and eliminate rows with no off-diag. nz. -C ---------------------------------------------------------------- - - DO 20 I = 1, N - - DEG = DEGREE (I) - - IF (DEG .GT. 0) THEN - -C ---------------------------------------------------------- -C place i in the degree list corresponding to its degree -C ---------------------------------------------------------- - - INEXT = HEAD (DEG) - IF (INEXT .NE. 0) LAST (INEXT) = I - NEXT (I) = INEXT - HEAD (DEG) = I - - ELSE - -C ---------------------------------------------------------- -C we have a variable that can be eliminated at once because -C there is no off-diagonal non-zero in its row. -C ---------------------------------------------------------- - - NEL = NEL + 1 - ELEN (I) = -NEL - PE (I) = 0 - W (I) = 0 - - ENDIF - -20 CONTINUE - -C======================================================================= -C WHILE (selecting pivots) DO -C======================================================================= - -30 CONTINUE - IF (NEL .LT. N) THEN - -C======================================================================= -C GET PIVOT OF MINIMUM DEGREE -C======================================================================= - -C ------------------------------------------------------------- -C find next supervariable for elimination -C ------------------------------------------------------------- - - DO 40 DEG = MINDEG, N - ME = HEAD (DEG) - IF (ME .GT. 0) GOTO 50 -40 CONTINUE -50 CONTINUE - MINDEG = DEG - -C ------------------------------------------------------------- -C remove chosen variable from link list -C ------------------------------------------------------------- - - INEXT = NEXT (ME) - IF (INEXT .NE. 0) LAST (INEXT) = 0 - HEAD (DEG) = INEXT - -C ------------------------------------------------------------- -C me represents the elimination of pivots nel+1 to nel+nv(me). -C place me itself as the first in this set. It will be moved -C to the nel+nv(me) position when the permutation vectors are -C computed. -C ------------------------------------------------------------- - - ELENME = ELEN (ME) - ELEN (ME) = - (NEL + 1) - NVPIV = NV (ME) - NEL = NEL + NVPIV - -C======================================================================= -C CONSTRUCT NEW ELEMENT -C======================================================================= - -C ------------------------------------------------------------- -C At this point, me is the pivotal supervariable. It will be -C converted into the current element. Scan list of the -C pivotal supervariable, me, setting tree pointers and -C constructing new list of supervariables for the new element, -C me. p is a pointer to the current position in the old list. -C ------------------------------------------------------------- - -C flag the variable "me" as being in Lme by negating nv (me) - NV (ME) = -NVPIV - DEGME = 0 - - IF (ELENME .EQ. 0) THEN - -C ---------------------------------------------------------- -C construct the new element in place -C ---------------------------------------------------------- - - PME1 = PE (ME) - PME2 = PME1 - 1 - - DO 60 P = PME1, PME1 + LEN (ME) - 1 - I = IW (P) - NVI = NV (I) - IF (NVI .GT. 0) THEN - -C ---------------------------------------------------- -C i is a principal variable not yet placed in Lme. -C store i in new list -C ---------------------------------------------------- - - DEGME = DEGME + NVI -C flag i as being in Lme by negating nv (i) - NV (I) = -NVI - PME2 = PME2 + 1 - IW (PME2) = I - -C ---------------------------------------------------- -C remove variable i from degree list. -C ---------------------------------------------------- - - ILAST = LAST (I) - INEXT = NEXT (I) - IF (INEXT .NE. 0) LAST (INEXT) = ILAST - IF (ILAST .NE. 0) THEN - NEXT (ILAST) = INEXT - ELSE -C i is at the head of the degree list - HEAD (DEGREE (I)) = INEXT - ENDIF - - ENDIF -60 CONTINUE -C this element takes no new memory in iw: - NEWMEM = 0 - - ELSE - -C ---------------------------------------------------------- -C construct the new element in empty space, iw (pfree ...) -C ---------------------------------------------------------- - - P = PE (ME) - PME1 = PFREE - SLENME = LEN (ME) - ELENME - - DO 120 KNT1 = 1, ELENME + 1 - - IF (KNT1 .GT. ELENME) THEN -C search the supervariables in me. - E = ME - PJ = P - LN = SLENME - ELSE -C search the elements in me. - E = IW (P) - P = P + 1 - PJ = PE (E) - LN = LEN (E) - ENDIF - -C ------------------------------------------------------- -C search for different supervariables and add them to the -C new list, compressing when necessary. this loop is -C executed once for each element in the list and once for -C all the supervariables in the list. -C ------------------------------------------------------- - - DO 110 KNT2 = 1, LN - I = IW (PJ) - PJ = PJ + 1 - NVI = NV (I) - IF (NVI .GT. 0) THEN - -C ------------------------------------------------- -C compress iw, if necessary -C ------------------------------------------------- - - IF (PFREE .GT. IWLEN) THEN -C prepare for compressing iw by adjusting -C pointers and lengths so that the lists being -C searched in the inner and outer loops contain -C only the remaining entries. - - PE (ME) = P - LEN (ME) = LEN (ME) - KNT1 - IF (LEN (ME) .EQ. 0) THEN -C nothing left of supervariable me - PE (ME) = 0 - ENDIF - PE (E) = PJ - LEN (E) = LN - KNT2 - IF (LEN (E) .EQ. 0) THEN -C nothing left of element e - PE (E) = 0 - ENDIF - - NCMPA = NCMPA + 1 -C store first item in pe -C set first entry to -item - DO 70 J = 1, N - PN = PE (J) - IF (PN .GT. 0) THEN - PE (J) = IW (PN) - IW (PN) = -J - ENDIF -70 CONTINUE - -C psrc/pdst point to source/destination - PDST = 1 - PSRC = 1 - PEND = PME1 - 1 - -C while loop: -80 CONTINUE - IF (PSRC .LE. PEND) THEN -C search for next negative entry - J = -IW (PSRC) - PSRC = PSRC + 1 - IF (J .GT. 0) THEN - IW (PDST) = PE (J) - PE (J) = PDST - PDST = PDST + 1 -C copy from source to destination - LENJ = LEN (J) - DO 90 KNT3 = 0, LENJ - 2 - IW (PDST + KNT3) = IW (PSRC + KNT3) -90 CONTINUE - PDST = PDST + LENJ - 1 - PSRC = PSRC + LENJ - 1 - ENDIF - GOTO 80 - ENDIF - -C move the new partially-constructed element - P1 = PDST - DO 100 PSRC = PME1, PFREE - 1 - IW (PDST) = IW (PSRC) - PDST = PDST + 1 -100 CONTINUE - PME1 = P1 - PFREE = PDST - PJ = PE (E) - P = PE (ME) - ENDIF - -C ------------------------------------------------- -C i is a principal variable not yet placed in Lme -C store i in new list -C ------------------------------------------------- - - DEGME = DEGME + NVI -C flag i as being in Lme by negating nv (i) - NV (I) = -NVI - IW (PFREE) = I - PFREE = PFREE + 1 - -C ------------------------------------------------- -C remove variable i from degree link list -C ------------------------------------------------- - - ILAST = LAST (I) - INEXT = NEXT (I) - IF (INEXT .NE. 0) LAST (INEXT) = ILAST - IF (ILAST .NE. 0) THEN - NEXT (ILAST) = INEXT - ELSE -C i is at the head of the degree list - HEAD (DEGREE (I)) = INEXT - ENDIF - - ENDIF -110 CONTINUE - - IF (E .NE. ME) THEN -C set tree pointer and flag to indicate element e is -C absorbed into new element me (the parent of e is me) - PE (E) = -ME - W (E) = 0 - ENDIF -120 CONTINUE - - PME2 = PFREE - 1 -C this element takes newmem new memory in iw (possibly zero) - NEWMEM = PFREE - PME1 - MEM = MEM + NEWMEM - MAXMEM = MAX (MAXMEM, MEM) - ENDIF - -C ------------------------------------------------------------- -C me has now been converted into an element in iw (pme1..pme2) -C ------------------------------------------------------------- - -C degme holds the external degree of new element - DEGREE (ME) = DEGME - PE (ME) = PME1 - LEN (ME) = PME2 - PME1 + 1 - -C ------------------------------------------------------------- -C make sure that wflg is not too large. With the current -C value of wflg, wflg+n must not cause integer overflow -C ------------------------------------------------------------- - - IF (WFLG + N .LE. WFLG) THEN - DO 130 X = 1, N - IF (W (X) .NE. 0) W (X) = 1 -130 CONTINUE - WFLG = 2 - ENDIF - -C======================================================================= -C COMPUTE (w (e) - wflg) = |Le\Lme| FOR ALL ELEMENTS -C======================================================================= - -C ------------------------------------------------------------- -C Scan 1: compute the external degrees of previous elements -C with respect to the current element. That is: -C (w (e) - wflg) = |Le \ Lme| -C for each element e that appears in any supervariable in Lme. -C The notation Le refers to the pattern (list of -C supervariables) of a previous element e, where e is not yet -C absorbed, stored in iw (pe (e) + 1 ... pe (e) + iw (pe (e))). -C The notation Lme refers to the pattern of the current element -C (stored in iw (pme1..pme2)). If (w (e) - wflg) becomes -C zero, then the element e will be absorbed in scan 2. -C ------------------------------------------------------------- - - DO 150 PME = PME1, PME2 - I = IW (PME) - ELN = ELEN (I) - IF (ELN .GT. 0) THEN -C note that nv (i) has been negated to denote i in Lme: - NVI = -NV (I) - WNVI = WFLG - NVI - DO 140 P = PE (I), PE (I) + ELN - 1 - E = IW (P) - WE = W (E) - IF (WE .GE. WFLG) THEN -C unabsorbed element e has been seen in this loop - WE = WE - NVI - ELSE IF (WE .NE. 0) THEN -C e is an unabsorbed element -C this is the first we have seen e in all of Scan 1 - WE = DEGREE (E) + WNVI - ENDIF - W (E) = WE -140 CONTINUE - ENDIF -150 CONTINUE - -C======================================================================= -C DEGREE UPDATE AND ELEMENT ABSORPTION -C======================================================================= - -C ------------------------------------------------------------- -C Scan 2: for each i in Lme, sum up the degree of Lme (which -C is degme), plus the sum of the external degrees of each Le -C for the elements e appearing within i, plus the -C supervariables in i. Place i in hash list. -C ------------------------------------------------------------- - - DO 180 PME = PME1, PME2 - I = IW (PME) - P1 = PE (I) - P2 = P1 + ELEN (I) - 1 - PN = P1 - HASH = 0 - DEG = 0 - -C ---------------------------------------------------------- -C scan the element list associated with supervariable i -C ---------------------------------------------------------- - - DO 160 P = P1, P2 - E = IW (P) -C dext = | Le \ Lme | - DEXT = W (E) - WFLG - IF (DEXT .GT. 0) THEN - DEG = DEG + DEXT - IW (PN) = E - PN = PN + 1 - HASH = HASH + E - ELSE IF (DEXT .EQ. 0) THEN -C aggressive absorption: e is not adjacent to me, but -C the |Le \ Lme| is 0, so absorb it into me - PE (E) = -ME - W (E) = 0 - ELSE -C element e has already been absorbed, due to -C regular absorption, in do loop 120 above. Ignore it. - CONTINUE - ENDIF -160 CONTINUE - -C count the number of elements in i (including me): - ELEN (I) = PN - P1 + 1 - -C ---------------------------------------------------------- -C scan the supervariables in the list associated with i -C ---------------------------------------------------------- - - P3 = PN - DO 170 P = P2 + 1, P1 + LEN (I) - 1 - J = IW (P) - NVJ = NV (J) - IF (NVJ .GT. 0) THEN -C j is unabsorbed, and not in Lme. -C add to degree and add to new list - DEG = DEG + NVJ - IW (PN) = J - PN = PN + 1 - HASH = HASH + J - ENDIF -170 CONTINUE - -C ---------------------------------------------------------- -C update the degree and check for mass elimination -C ---------------------------------------------------------- - - IF (DEG .EQ. 0) THEN - -C ------------------------------------------------------- -C mass elimination -C ------------------------------------------------------- - -C There is nothing left of this node except for an -C edge to the current pivot element. elen (i) is 1, -C and there are no variables adjacent to node i. -C Absorb i into the current pivot element, me. - - PE (I) = -ME - NVI = -NV (I) - DEGME = DEGME - NVI - NVPIV = NVPIV + NVI - NEL = NEL + NVI - NV (I) = 0 - ELEN (I) = 0 - - ELSE - -C ------------------------------------------------------- -C update the upper-bound degree of i -C ------------------------------------------------------- - -C the following degree does not yet include the size -C of the current element, which is added later: - DEGREE (I) = MIN (DEGREE (I), DEG) - -C ------------------------------------------------------- -C add me to the list for i -C ------------------------------------------------------- - -C move first supervariable to end of list - IW (PN) = IW (P3) -C move first element to end of element part of list - IW (P3) = IW (P1) -C add new element to front of list. - IW (P1) = ME -C store the new length of the list in len (i) - LEN (I) = PN - P1 + 1 - -C ------------------------------------------------------- -C place in hash bucket. Save hash key of i in last (i). -C ------------------------------------------------------- - - HASH = MOD (HASH, HMOD) + 1 - J = HEAD (HASH) - IF (J .LE. 0) THEN -C the degree list is empty, hash head is -j - NEXT (I) = -J - HEAD (HASH) = -I - ELSE -C degree list is not empty -C use last (head (hash)) as hash head - NEXT (I) = LAST (J) - LAST (J) = I - ENDIF - LAST (I) = HASH - ENDIF -180 CONTINUE - - DEGREE (ME) = DEGME - -C ------------------------------------------------------------- -C Clear the counter array, w (...), by incrementing wflg. -C ------------------------------------------------------------- - - DMAX = MAX (DMAX, DEGME) - WFLG = WFLG + DMAX - -C make sure that wflg+n does not cause integer overflow - IF (WFLG + N .LE. WFLG) THEN - DO 190 X = 1, N - IF (W (X) .NE. 0) W (X) = 1 -190 CONTINUE - WFLG = 2 - ENDIF -C at this point, w (1..n) .lt. wflg holds - -C======================================================================= -C SUPERVARIABLE DETECTION -C======================================================================= - - DO 250 PME = PME1, PME2 - I = IW (PME) - IF (NV (I) .LT. 0) THEN -C i is a principal variable in Lme - -C ------------------------------------------------------- -C examine all hash buckets with 2 or more variables. We -C do this by examing all unique hash keys for super- -C variables in the pattern Lme of the current element, me -C ------------------------------------------------------- - - HASH = LAST (I) -C let i = head of hash bucket, and empty the hash bucket - J = HEAD (HASH) - IF (J .EQ. 0) GOTO 250 - IF (J .LT. 0) THEN -C degree list is empty - I = -J - HEAD (HASH) = 0 - ELSE -C degree list is not empty, restore last () of head - I = LAST (J) - LAST (J) = 0 - ENDIF - IF (I .EQ. 0) GOTO 250 - -C while loop: -200 CONTINUE - IF (NEXT (I) .NE. 0) THEN - -C ---------------------------------------------------- -C this bucket has one or more variables following i. -C scan all of them to see if i can absorb any entries -C that follow i in hash bucket. Scatter i into w. -C ---------------------------------------------------- - - LN = LEN (I) - ELN = ELEN (I) -C do not flag the first element in the list (me) - DO 210 P = PE (I) + 1, PE (I) + LN - 1 - W (IW (P)) = WFLG -210 CONTINUE - -C ---------------------------------------------------- -C scan every other entry j following i in bucket -C ---------------------------------------------------- - - JLAST = I - J = NEXT (I) - -C while loop: -220 CONTINUE - IF (J .NE. 0) THEN - -C ------------------------------------------------- -C check if j and i have identical nonzero pattern -C ------------------------------------------------- - - IF (LEN (J) .NE. LN) THEN -C i and j do not have same size data structure - GOTO 240 - ENDIF - IF (ELEN (J) .NE. ELN) THEN -C i and j do not have same number of adjacent el - GOTO 240 - ENDIF -C do not flag the first element in the list (me) - DO 230 P = PE (J) + 1, PE (J) + LN - 1 - IF (W (IW (P)) .NE. WFLG) THEN -C an entry (iw(p)) is in j but not in i - GOTO 240 - ENDIF -230 CONTINUE - -C ------------------------------------------------- -C found it! j can be absorbed into i -C ------------------------------------------------- - - PE (J) = -I -C both nv (i) and nv (j) are negated since they -C are in Lme, and the absolute values of each -C are the number of variables in i and j: - NV (I) = NV (I) + NV (J) - NV (J) = 0 - ELEN (J) = 0 -C delete j from hash bucket - J = NEXT (J) - NEXT (JLAST) = J - GOTO 220 - -C ------------------------------------------------- -240 CONTINUE -C j cannot be absorbed into i -C ------------------------------------------------- - - JLAST = J - J = NEXT (J) - GOTO 220 - ENDIF - -C ---------------------------------------------------- -C no more variables can be absorbed into i -C go to next i in bucket and clear flag array -C ---------------------------------------------------- - - WFLG = WFLG + 1 - I = NEXT (I) - IF (I .NE. 0) GOTO 200 - ENDIF - ENDIF -250 CONTINUE - -C======================================================================= -C RESTORE DEGREE LISTS AND REMOVE NONPRINCIPAL SUPERVAR. FROM ELEMENT -C======================================================================= - - P = PME1 - NLEFT = N - NEL - DO 260 PME = PME1, PME2 - I = IW (PME) - NVI = -NV (I) - IF (NVI .GT. 0) THEN -C i is a principal variable in Lme -C restore nv (i) to signify that i is principal - NV (I) = NVI - -C ------------------------------------------------------- -C compute the external degree (add size of current elem) -C ------------------------------------------------------- - - DEG = MIN (DEGREE (I) + DEGME - NVI, NLEFT - NVI) - -C ------------------------------------------------------- -C place the supervariable at the head of the degree list -C ------------------------------------------------------- - - INEXT = HEAD (DEG) - IF (INEXT .NE. 0) LAST (INEXT) = I - NEXT (I) = INEXT - LAST (I) = 0 - HEAD (DEG) = I - -C ------------------------------------------------------- -C save the new degree, and find the minimum degree -C ------------------------------------------------------- - - MINDEG = MIN (MINDEG, DEG) - DEGREE (I) = DEG - -C ------------------------------------------------------- -C place the supervariable in the element pattern -C ------------------------------------------------------- - - IW (P) = I - P = P + 1 - ENDIF -260 CONTINUE - -C======================================================================= -C FINALIZE THE NEW ELEMENT -C======================================================================= - - NV (ME) = NVPIV + DEGME -C nv (me) is now the degree of pivot (including diagonal part) -C save the length of the list for the new element me - LEN (ME) = P - PME1 - IF (LEN (ME) .EQ. 0) THEN -C there is nothing left of the current pivot element - PE (ME) = 0 - W (ME) = 0 - ENDIF - IF (NEWMEM .NE. 0) THEN -C element was not constructed in place: deallocate part -C of it (final size is less than or equal to newmem, -C since newly nonprincipal variables have been removed). - PFREE = P - MEM = MEM - NEWMEM + LEN (ME) - ENDIF - -C======================================================================= -C END WHILE (selecting pivots) - GOTO 30 - ENDIF -C======================================================================= - -C======================================================================= -C COMPUTE THE PERMUTATION VECTORS -C======================================================================= - -C ---------------------------------------------------------------- -C The time taken by the following code is O(n). At this -C point, elen (e) = -k has been done for all elements e, -C and elen (i) = 0 has been done for all nonprincipal -C variables i. At this point, there are no principal -C supervariables left, and all elements are absorbed. -C ---------------------------------------------------------------- - -C ---------------------------------------------------------------- -C compute the ordering of unordered nonprincipal variables -C ---------------------------------------------------------------- - - DO 290 I = 1, N - IF (ELEN (I) .EQ. 0) THEN - -C ---------------------------------------------------------- -C i is an un-ordered row. Traverse the tree from i until -C reaching an element, e. The element, e, was the -C principal supervariable of i and all nodes in the path -C from i to when e was selected as pivot. -C ---------------------------------------------------------- - - J = -PE (I) -C while (j is a variable) do: -270 CONTINUE - IF (ELEN (J) .GE. 0) THEN - J = -PE (J) - GOTO 270 - ENDIF - E = J - -C ---------------------------------------------------------- -C get the current pivot ordering of e -C ---------------------------------------------------------- - - K = -ELEN (E) - -C ---------------------------------------------------------- -C traverse the path again from i to e, and compress the -C path (all nodes point to e). Path compression allows -C this code to compute in O(n) time. Order the unordered -C nodes in the path, and place the element e at the end. -C ---------------------------------------------------------- - - J = I -C while (j is a variable) do: -280 CONTINUE - IF (ELEN (J) .GE. 0) THEN - JNEXT = -PE (J) - PE (J) = -E - IF (ELEN (J) .EQ. 0) THEN -C j is an unordered row - ELEN (J) = K - K = K + 1 - ENDIF - J = JNEXT - GOTO 280 - ENDIF -C leave elen (e) negative, so we know it is an element - ELEN (E) = -K - ENDIF -290 CONTINUE - -C ---------------------------------------------------------------- -C reset the inverse permutation (elen (1..n)) to be positive, -C and compute the permutation (last (1..n)). -C ---------------------------------------------------------------- - - DO 300 I = 1, N - K = ABS (ELEN (I)) - LAST (K) = I - ELEN (I) = K -300 CONTINUE - -C======================================================================= -C RETURN THE MEMORY USAGE IN IW -C======================================================================= - -C If maxmem is less than or equal to iwlen, then no compressions -C occurred, and iw (maxmem+1 ... iwlen) was unused. Otherwise -C compressions did occur, and iwlen would have had to have been -C greater than or equal to maxmem for no compressions to occur. -C Return the value of maxmem in the pfree argument. - - PFREE = MAXMEM - - RETURN - END - diff --git a/ThirdParty/SuiteSparse/AMD/Source/amdbar.f b/ThirdParty/SuiteSparse/AMD/Source/amdbar.f deleted file mode 100644 index 45b829743f..0000000000 --- a/ThirdParty/SuiteSparse/AMD/Source/amdbar.f +++ /dev/null @@ -1,1212 +0,0 @@ -C----------------------------------------------------------------------- -C AMD/Source/amdbar.f: Fortran version of AMD, no aggress absorption -C----------------------------------------------------------------------- - -C AMD, Copyright (c) 1996-2022, Timothy A. Davis, Patrick R. Amestoy, -C and C Iain S. Duff. All Rights Reserved. -C SPDX-License-Identifier: BSD-3-clause - -C----------------------------------------------------------------------- - - SUBROUTINE AMDBAR - $ (N, PE, IW, LEN, IWLEN, PFREE, NV, NEXT, - $ LAST, HEAD, ELEN, DEGREE, NCMPA, W) - - INTEGER N, IWLEN, PFREE, NCMPA, IW (IWLEN), PE (N), - $ DEGREE (N), NV (N), NEXT (N), LAST (N), HEAD (N), - $ ELEN (N), W (N), LEN (N) - -C Given a representation of the nonzero pattern of a symmetric matrix, -C A, (excluding the diagonal) perform an approximate minimum -C (UMFPACK/MA38-style) degree ordering to compute a pivot order -C such that the introduction of nonzeros (fill-in) in the Cholesky -C factors A = LL^T are kept low. At each step, the pivot -C selected is the one with the minimum UMFPACK/MA38-style -C upper-bound on the external degree. -C -C This routine does not do aggresive absorption (as done by AMD). - -C ********************************************************************** -C ***** CAUTION: ARGUMENTS ARE NOT CHECKED FOR ERRORS ON INPUT. ****** -C ********************************************************************** - -C References: -C -C [1] Timothy A. Davis and Iain Duff, "An unsymmetric-pattern -C multifrontal method for sparse LU factorization", SIAM J. -C Matrix Analysis and Applications, vol. 18, no. 1, pp. -C 140-158. Discusses UMFPACK / MA38, which first introduced -C the approximate minimum degree used by this routine. -C -C [2] Patrick Amestoy, Timothy A. Davis, and Iain S. Duff, "An -C approximate degree ordering algorithm," SIAM J. Matrix -C Analysis and Applications, vol. 17, no. 4, pp. 886-905, -C 1996. Discusses AMD, AMDBAR, and MC47B. -C -C [3] Alan George and Joseph Liu, "The evolution of the minimum -C degree ordering algorithm," SIAM Review, vol. 31, no. 1, -C pp. 1-19, 1989. We list below the features mentioned in -C that paper that this code includes: -C -C mass elimination: -C Yes. MA27 relied on supervariable detection for mass -C elimination. -C indistinguishable nodes: -C Yes (we call these "supervariables"). This was also in -C the MA27 code - although we modified the method of -C detecting them (the previous hash was the true degree, -C which we no longer keep track of). A supervariable is -C a set of rows with identical nonzero pattern. All -C variables in a supervariable are eliminated together. -C Each supervariable has as its numerical name that of -C one of its variables (its principal variable). -C quotient graph representation: -C Yes. We use the term "element" for the cliques formed -C during elimination. This was also in the MA27 code. -C The algorithm can operate in place, but it will work -C more efficiently if given some "elbow room." -C element absorption: -C Yes. This was also in the MA27 code. -C external degree: -C Yes. The MA27 code was based on the true degree. -C incomplete degree update and multiple elimination: -C No. This was not in MA27, either. Our method of -C degree update within MC47B/BD is element-based, not -C variable-based. It is thus not well-suited for use -C with incomplete degree update or multiple elimination. - -C----------------------------------------------------------------------- -C Authors, and Copyright (C) 1995 by: -C Timothy A. Davis, Patrick Amestoy, Iain S. Duff, & John K. Reid. -C -C Acknowledgements: -C This work (and the UMFPACK package) was supported by the -C National Science Foundation (ASC-9111263 and DMS-9223088). -C The UMFPACK/MA38 approximate degree update algorithm, the -C unsymmetric analog which forms the basis of MC47B/BD, was -C developed while Tim Davis was supported by CERFACS (Toulouse, -C France) in a post-doctoral position. -C -C Date: September, 1995 -C----------------------------------------------------------------------- - -C----------------------------------------------------------------------- -C INPUT ARGUMENTS (unaltered): -C----------------------------------------------------------------------- - -C n: The matrix order. -C -C Restriction: 1 .le. n .lt. (iovflo/2)-2, where iovflo is -C the largest positive integer that your computer can represent. - -C iwlen: The length of iw (1..iwlen). On input, the matrix is -C stored in iw (1..pfree-1). However, iw (1..iwlen) should be -C slightly larger than what is required to hold the matrix, at -C least iwlen .ge. pfree + n is recommended. Otherwise, -C excessive compressions will take place. -C *** We do not recommend running this algorithm with *** -C *** iwlen .lt. pfree + n. *** -C *** Better performance will be obtained if *** -C *** iwlen .ge. pfree + n *** -C *** or better yet *** -C *** iwlen .gt. 1.2 * pfree *** -C *** (where pfree is its value on input). *** -C The algorithm will not run at all if iwlen .lt. pfree-1. -C -C Restriction: iwlen .ge. pfree-1 - -C----------------------------------------------------------------------- -C INPUT/OUPUT ARGUMENTS: -C----------------------------------------------------------------------- - -C pe: On input, pe (i) is the index in iw of the start of row i, or -C zero if row i has no off-diagonal non-zeros. -C -C During execution, it is used for both supervariables and -C elements: -C -C * Principal supervariable i: index into iw of the -C description of supervariable i. A supervariable -C represents one or more rows of the matrix -C with identical nonzero pattern. -C * Non-principal supervariable i: if i has been absorbed -C into another supervariable j, then pe (i) = -j. -C That is, j has the same pattern as i. -C Note that j might later be absorbed into another -C supervariable j2, in which case pe (i) is still -j, -C and pe (j) = -j2. -C * Unabsorbed element e: the index into iw of the description -C of element e, if e has not yet been absorbed by a -C subsequent element. Element e is created when -C the supervariable of the same name is selected as -C the pivot. -C * Absorbed element e: if element e is absorbed into element -C e2, then pe (e) = -e2. This occurs when the pattern of -C e (that is, Le) is found to be a subset of the pattern -C of e2 (that is, Le2). If element e is "null" (it has -C no nonzeros outside its pivot block), then pe (e) = 0. -C -C On output, pe holds the assembly tree/forest, which implicitly -C represents a pivot order with identical fill-in as the actual -C order (via a depth-first search of the tree). -C -C On output: -C If nv (i) .gt. 0, then i represents a node in the assembly tree, -C and the parent of i is -pe (i), or zero if i is a root. -C If nv (i) = 0, then (i,-pe (i)) represents an edge in a -C subtree, the root of which is a node in the assembly tree. - -C pfree: On input the tail end of the array, iw (pfree..iwlen), -C is empty, and the matrix is stored in iw (1..pfree-1). -C During execution, additional data is placed in iw, and pfree -C is modified so that iw (pfree..iwlen) is always the unused part -C of iw. On output, pfree is set equal to the size of iw that -C would have been needed for no compressions to occur. If -C ncmpa is zero, then pfree (on output) is less than or equal to -C iwlen, and the space iw (pfree+1 ... iwlen) was not used. -C Otherwise, pfree (on output) is greater than iwlen, and all the -C memory in iw was used. - -C----------------------------------------------------------------------- -C INPUT/MODIFIED (undefined on output): -C----------------------------------------------------------------------- - -C len: On input, len (i) holds the number of entries in row i of the -C matrix, excluding the diagonal. The contents of len (1..n) -C are undefined on output. - -C iw: On input, iw (1..pfree-1) holds the description of each row i -C in the matrix. The matrix must be symmetric, and both upper -C and lower triangular parts must be present. The diagonal must -C not be present. Row i is held as follows: -C -C len (i): the length of the row i data structure -C iw (pe (i) ... pe (i) + len (i) - 1): -C the list of column indices for nonzeros -C in row i (simple supervariables), excluding -C the diagonal. All supervariables start with -C one row/column each (supervariable i is just -C row i). -C if len (i) is zero on input, then pe (i) is ignored -C on input. -C -C Note that the rows need not be in any particular order, -C and there may be empty space between the rows. -C -C During execution, the supervariable i experiences fill-in. -C This is represented by placing in i a list of the elements -C that cause fill-in in supervariable i: -C -C len (i): the length of supervariable i -C iw (pe (i) ... pe (i) + elen (i) - 1): -C the list of elements that contain i. This list -C is kept short by removing absorbed elements. -C iw (pe (i) + elen (i) ... pe (i) + len (i) - 1): -C the list of supervariables in i. This list -C is kept short by removing nonprincipal -C variables, and any entry j that is also -C contained in at least one of the elements -C (j in Le) in the list for i (e in row i). -C -C When supervariable i is selected as pivot, we create an -C element e of the same name (e=i): -C -C len (e): the length of element e -C iw (pe (e) ... pe (e) + len (e) - 1): -C the list of supervariables in element e. -C -C An element represents the fill-in that occurs when supervariable -C i is selected as pivot (which represents the selection of row i -C and all non-principal variables whose principal variable is i). -C We use the term Le to denote the set of all supervariables -C in element e. Absorbed supervariables and elements are pruned -C from these lists when computationally convenient. -C -C CAUTION: THE INPUT MATRIX IS OVERWRITTEN DURING COMPUTATION. -C The contents of iw are undefined on output. - -C----------------------------------------------------------------------- -C OUTPUT (need not be set on input): -C----------------------------------------------------------------------- - -C nv: During execution, abs (nv (i)) is equal to the number of rows -C that are represented by the principal supervariable i. If i is -C a nonprincipal variable, then nv (i) = 0. Initially, -C nv (i) = 1 for all i. nv (i) .lt. 0 signifies that i is a -C principal variable in the pattern Lme of the current pivot -C element me. On output, nv (e) holds the true degree of element -C e at the time it was created (including the diagonal part). - -C ncmpa: The number of times iw was compressed. If this is -C excessive, then the execution took longer than what could have -C been. To reduce ncmpa, try increasing iwlen to be 10% or 20% -C larger than the value of pfree on input (or at least -C iwlen .ge. pfree + n). The fastest performance will be -C obtained when ncmpa is returned as zero. If iwlen is set to -C the value returned by pfree on *output*, then no compressions -C will occur. - -C elen: See the description of iw above. At the start of execution, -C elen (i) is set to zero. During execution, elen (i) is the -C number of elements in the list for supervariable i. When e -C becomes an element, elen (e) = -nel is set, where nel is the -C current step of factorization. elen (i) = 0 is done when i -C becomes nonprincipal. -C -C For variables, elen (i) .ge. 0 holds until just before the -C permutation vectors are computed. For elements, -C elen (e) .lt. 0 holds. -C -C On output elen (1..n) holds the inverse permutation (the same -C as the 'INVP' argument in Sparspak). That is, if k = elen (i), -C then row i is the kth pivot row. Row i of A appears as the -C (elen(i))-th row in the permuted matrix, PAP^T. - -C last: In a degree list, last (i) is the supervariable preceding i, -C or zero if i is the head of the list. In a hash bucket, -C last (i) is the hash key for i. last (head (hash)) is also -C used as the head of a hash bucket if head (hash) contains a -C degree list (see head, below). -C -C On output, last (1..n) holds the permutation (the same as the -C 'PERM' argument in Sparspak). That is, if i = last (k), then -C row i is the kth pivot row. Row last (k) of A is the k-th row -C in the permuted matrix, PAP^T. - -C----------------------------------------------------------------------- -C LOCAL (not input or output - used only during execution): -C----------------------------------------------------------------------- - -C degree: If i is a supervariable, then degree (i) holds the -C current approximation of the external degree of row i (an upper -C bound). The external degree is the number of nonzeros in row i, -C minus abs (nv (i)) (the diagonal part). The bound is equal to -C the external degree if elen (i) is less than or equal to two. -C -C We also use the term "external degree" for elements e to refer -C to |Le \ Lme|. If e is an element, then degree (e) holds |Le|, -C which is the degree of the off-diagonal part of the element e -C (not including the diagonal part). - -C head: head is used for degree lists. head (deg) is the first -C supervariable in a degree list (all supervariables i in a -C degree list deg have the same approximate degree, namely, -C deg = degree (i)). If the list deg is empty then -C head (deg) = 0. -C -C During supervariable detection head (hash) also serves as a -C pointer to a hash bucket. -C If head (hash) .gt. 0, there is a degree list of degree hash. -C The hash bucket head pointer is last (head (hash)). -C If head (hash) = 0, then the degree list and hash bucket are -C both empty. -C If head (hash) .lt. 0, then the degree list is empty, and -C -head (hash) is the head of the hash bucket. -C After supervariable detection is complete, all hash buckets -C are empty, and the (last (head (hash)) = 0) condition is -C restored for the non-empty degree lists. - -C next: next (i) is the supervariable following i in a link list, or -C zero if i is the last in the list. Used for two kinds of -C lists: degree lists and hash buckets (a supervariable can be -C in only one kind of list at a time). - -C w: The flag array w determines the status of elements and -C variables, and the external degree of elements. -C -C for elements: -C if w (e) = 0, then the element e is absorbed -C if w (e) .ge. wflg, then w (e) - wflg is the size of -C the set |Le \ Lme|, in terms of nonzeros (the -C sum of abs (nv (i)) for each principal variable i that -C is both in the pattern of element e and NOT in the -C pattern of the current pivot element, me). -C if wflg .gt. w (e) .gt. 0, then e is not absorbed and has -C not yet been seen in the scan of the element lists in -C the computation of |Le\Lme| in loop 150 below. -C -C for variables: -C during supervariable detection, if w (j) .ne. wflg then j is -C not in the pattern of variable i -C -C The w array is initialized by setting w (i) = 1 for all i, -C and by setting wflg = 2. It is reinitialized if wflg becomes -C too large (to ensure that wflg+n does not cause integer -C overflow). - -C----------------------------------------------------------------------- -C LOCAL INTEGERS: -C----------------------------------------------------------------------- - - INTEGER DEG, DEGME, DMAX, E, ELENME, ELN, HASH, HMOD, I, - $ ILAST, INEXT, J, JLAST, JNEXT, K, KNT1, KNT2, KNT3, - $ LENJ, LN, MAXMEM, ME, MEM, MINDEG, NEL, NEWMEM, - $ NLEFT, NVI, NVJ, NVPIV, SLENME, WE, WFLG, WNVI, X - -C deg: the degree of a variable or element -C degme: size, |Lme|, of the current element, me (= degree (me)) -C dext: external degree, |Le \ Lme|, of some element e -C dmax: largest |Le| seen so far -C e: an element -C elenme: the length, elen (me), of element list of pivotal var. -C eln: the length, elen (...), of an element list -C hash: the computed value of the hash function -C hmod: the hash function is computed modulo hmod = max (1,n-1) -C i: a supervariable -C ilast: the entry in a link list preceding i -C inext: the entry in a link list following i -C j: a supervariable -C jlast: the entry in a link list preceding j -C jnext: the entry in a link list, or path, following j -C k: the pivot order of an element or variable -C knt1: loop counter used during element construction -C knt2: loop counter used during element construction -C knt3: loop counter used during compression -C lenj: len (j) -C ln: length of a supervariable list -C maxmem: amount of memory needed for no compressions -C me: current supervariable being eliminated, and the -C current element created by eliminating that -C supervariable -C mem: memory in use assuming no compressions have occurred -C mindeg: current minimum degree -C nel: number of pivots selected so far -C newmem: amount of new memory needed for current pivot element -C nleft: n - nel, the number of nonpivotal rows/columns remaining -C nvi: the number of variables in a supervariable i (= nv (i)) -C nvj: the number of variables in a supervariable j (= nv (j)) -C nvpiv: number of pivots in current element -C slenme: number of variables in variable list of pivotal variable -C we: w (e) -C wflg: used for flagging the w array. See description of iw. -C wnvi: wflg - nv (i) -C x: either a supervariable or an element - -C----------------------------------------------------------------------- -C LOCAL POINTERS: -C----------------------------------------------------------------------- - - INTEGER P, P1, P2, P3, PDST, PEND, PJ, PME, PME1, PME2, PN, PSRC - -C Any parameter (pe (...) or pfree) or local variable -C starting with "p" (for Pointer) is an index into iw, -C and all indices into iw use variables starting with -C "p." The only exception to this rule is the iwlen -C input argument. - -C p: pointer into lots of things -C p1: pe (i) for some variable i (start of element list) -C p2: pe (i) + elen (i) - 1 for some var. i (end of el. list) -C p3: index of first supervariable in clean list -C pdst: destination pointer, for compression -C pend: end of memory to compress -C pj: pointer into an element or variable -C pme: pointer into the current element (pme1...pme2) -C pme1: the current element, me, is stored in iw (pme1...pme2) -C pme2: the end of the current element -C pn: pointer into a "clean" variable, also used to compress -C psrc: source pointer, for compression - -C----------------------------------------------------------------------- -C FUNCTIONS CALLED: -C----------------------------------------------------------------------- - - INTRINSIC MAX, MIN, MOD - -C======================================================================= -C INITIALIZATIONS -C======================================================================= - - WFLG = 2 - MINDEG = 1 - NCMPA = 0 - NEL = 0 - HMOD = MAX (1, N-1) - DMAX = 0 - MEM = PFREE - 1 - MAXMEM = MEM - ME = 0 - - DO 10 I = 1, N - LAST (I) = 0 - HEAD (I) = 0 - NV (I) = 1 - W (I) = 1 - ELEN (I) = 0 - DEGREE (I) = LEN (I) -10 CONTINUE - -C ---------------------------------------------------------------- -C initialize degree lists and eliminate rows with no off-diag. nz. -C ---------------------------------------------------------------- - - DO 20 I = 1, N - - DEG = DEGREE (I) - - IF (DEG .GT. 0) THEN - -C ---------------------------------------------------------- -C place i in the degree list corresponding to its degree -C ---------------------------------------------------------- - - INEXT = HEAD (DEG) - IF (INEXT .NE. 0) LAST (INEXT) = I - NEXT (I) = INEXT - HEAD (DEG) = I - - ELSE - -C ---------------------------------------------------------- -C we have a variable that can be eliminated at once because -C there is no off-diagonal non-zero in its row. -C ---------------------------------------------------------- - - NEL = NEL + 1 - ELEN (I) = -NEL - PE (I) = 0 - W (I) = 0 - - ENDIF - -20 CONTINUE - -C======================================================================= -C WHILE (selecting pivots) DO -C======================================================================= - -30 CONTINUE - IF (NEL .LT. N) THEN - -C======================================================================= -C GET PIVOT OF MINIMUM DEGREE -C======================================================================= - -C ------------------------------------------------------------- -C find next supervariable for elimination -C ------------------------------------------------------------- - - DO 40 DEG = MINDEG, N - ME = HEAD (DEG) - IF (ME .GT. 0) GOTO 50 -40 CONTINUE -50 CONTINUE - MINDEG = DEG - -C ------------------------------------------------------------- -C remove chosen variable from link list -C ------------------------------------------------------------- - - INEXT = NEXT (ME) - IF (INEXT .NE. 0) LAST (INEXT) = 0 - HEAD (DEG) = INEXT - -C ------------------------------------------------------------- -C me represents the elimination of pivots nel+1 to nel+nv(me). -C place me itself as the first in this set. It will be moved -C to the nel+nv(me) position when the permutation vectors are -C computed. -C ------------------------------------------------------------- - - ELENME = ELEN (ME) - ELEN (ME) = - (NEL + 1) - NVPIV = NV (ME) - NEL = NEL + NVPIV - -C======================================================================= -C CONSTRUCT NEW ELEMENT -C======================================================================= - -C ------------------------------------------------------------- -C At this point, me is the pivotal supervariable. It will be -C converted into the current element. Scan list of the -C pivotal supervariable, me, setting tree pointers and -C constructing new list of supervariables for the new element, -C me. p is a pointer to the current position in the old list. -C ------------------------------------------------------------- - -C flag the variable "me" as being in Lme by negating nv (me) - NV (ME) = -NVPIV - DEGME = 0 - - IF (ELENME .EQ. 0) THEN - -C ---------------------------------------------------------- -C construct the new element in place -C ---------------------------------------------------------- - - PME1 = PE (ME) - PME2 = PME1 - 1 - - DO 60 P = PME1, PME1 + LEN (ME) - 1 - I = IW (P) - NVI = NV (I) - IF (NVI .GT. 0) THEN - -C ---------------------------------------------------- -C i is a principal variable not yet placed in Lme. -C store i in new list -C ---------------------------------------------------- - - DEGME = DEGME + NVI -C flag i as being in Lme by negating nv (i) - NV (I) = -NVI - PME2 = PME2 + 1 - IW (PME2) = I - -C ---------------------------------------------------- -C remove variable i from degree list. -C ---------------------------------------------------- - - ILAST = LAST (I) - INEXT = NEXT (I) - IF (INEXT .NE. 0) LAST (INEXT) = ILAST - IF (ILAST .NE. 0) THEN - NEXT (ILAST) = INEXT - ELSE -C i is at the head of the degree list - HEAD (DEGREE (I)) = INEXT - ENDIF - - ENDIF -60 CONTINUE -C this element takes no new memory in iw: - NEWMEM = 0 - - ELSE - -C ---------------------------------------------------------- -C construct the new element in empty space, iw (pfree ...) -C ---------------------------------------------------------- - - P = PE (ME) - PME1 = PFREE - SLENME = LEN (ME) - ELENME - - DO 120 KNT1 = 1, ELENME + 1 - - IF (KNT1 .GT. ELENME) THEN -C search the supervariables in me. - E = ME - PJ = P - LN = SLENME - ELSE -C search the elements in me. - E = IW (P) - P = P + 1 - PJ = PE (E) - LN = LEN (E) - ENDIF - -C ------------------------------------------------------- -C search for different supervariables and add them to the -C new list, compressing when necessary. this loop is -C executed once for each element in the list and once for -C all the supervariables in the list. -C ------------------------------------------------------- - - DO 110 KNT2 = 1, LN - I = IW (PJ) - PJ = PJ + 1 - NVI = NV (I) - IF (NVI .GT. 0) THEN - -C ------------------------------------------------- -C compress iw, if necessary -C ------------------------------------------------- - - IF (PFREE .GT. IWLEN) THEN -C prepare for compressing iw by adjusting -C pointers and lengths so that the lists being -C searched in the inner and outer loops contain -C only the remaining entries. - - PE (ME) = P - LEN (ME) = LEN (ME) - KNT1 - IF (LEN (ME) .EQ. 0) THEN -C nothing left of supervariable me - PE (ME) = 0 - ENDIF - PE (E) = PJ - LEN (E) = LN - KNT2 - IF (LEN (E) .EQ. 0) THEN -C nothing left of element e - PE (E) = 0 - ENDIF - - NCMPA = NCMPA + 1 -C store first item in pe -C set first entry to -item - DO 70 J = 1, N - PN = PE (J) - IF (PN .GT. 0) THEN - PE (J) = IW (PN) - IW (PN) = -J - ENDIF -70 CONTINUE - -C psrc/pdst point to source/destination - PDST = 1 - PSRC = 1 - PEND = PME1 - 1 - -C while loop: -80 CONTINUE - IF (PSRC .LE. PEND) THEN -C search for next negative entry - J = -IW (PSRC) - PSRC = PSRC + 1 - IF (J .GT. 0) THEN - IW (PDST) = PE (J) - PE (J) = PDST - PDST = PDST + 1 -C copy from source to destination - LENJ = LEN (J) - DO 90 KNT3 = 0, LENJ - 2 - IW (PDST + KNT3) = IW (PSRC + KNT3) -90 CONTINUE - PDST = PDST + LENJ - 1 - PSRC = PSRC + LENJ - 1 - ENDIF - GOTO 80 - ENDIF - -C move the new partially-constructed element - P1 = PDST - DO 100 PSRC = PME1, PFREE - 1 - IW (PDST) = IW (PSRC) - PDST = PDST + 1 -100 CONTINUE - PME1 = P1 - PFREE = PDST - PJ = PE (E) - P = PE (ME) - ENDIF - -C ------------------------------------------------- -C i is a principal variable not yet placed in Lme -C store i in new list -C ------------------------------------------------- - - DEGME = DEGME + NVI -C flag i as being in Lme by negating nv (i) - NV (I) = -NVI - IW (PFREE) = I - PFREE = PFREE + 1 - -C ------------------------------------------------- -C remove variable i from degree link list -C ------------------------------------------------- - - ILAST = LAST (I) - INEXT = NEXT (I) - IF (INEXT .NE. 0) LAST (INEXT) = ILAST - IF (ILAST .NE. 0) THEN - NEXT (ILAST) = INEXT - ELSE -C i is at the head of the degree list - HEAD (DEGREE (I)) = INEXT - ENDIF - - ENDIF -110 CONTINUE - - IF (E .NE. ME) THEN -C set tree pointer and flag to indicate element e is -C absorbed into new element me (the parent of e is me) - PE (E) = -ME - W (E) = 0 - ENDIF -120 CONTINUE - - PME2 = PFREE - 1 -C this element takes newmem new memory in iw (possibly zero) - NEWMEM = PFREE - PME1 - MEM = MEM + NEWMEM - MAXMEM = MAX (MAXMEM, MEM) - ENDIF - -C ------------------------------------------------------------- -C me has now been converted into an element in iw (pme1..pme2) -C ------------------------------------------------------------- - -C degme holds the external degree of new element - DEGREE (ME) = DEGME - PE (ME) = PME1 - LEN (ME) = PME2 - PME1 + 1 - -C ------------------------------------------------------------- -C make sure that wflg is not too large. With the current -C value of wflg, wflg+n must not cause integer overflow -C ------------------------------------------------------------- - - IF (WFLG + N .LE. WFLG) THEN - DO 130 X = 1, N - IF (W (X) .NE. 0) W (X) = 1 -130 CONTINUE - WFLG = 2 - ENDIF - -C======================================================================= -C COMPUTE (w (e) - wflg) = |Le\Lme| FOR ALL ELEMENTS -C======================================================================= - -C ------------------------------------------------------------- -C Scan 1: compute the external degrees of previous elements -C with respect to the current element. That is: -C (w (e) - wflg) = |Le \ Lme| -C for each element e that appears in any supervariable in Lme. -C The notation Le refers to the pattern (list of -C supervariables) of a previous element e, where e is not yet -C absorbed, stored in iw (pe (e) + 1 ... pe (e) + iw (pe (e))). -C The notation Lme refers to the pattern of the current element -C (stored in iw (pme1..pme2)). If (w (e) - wflg) becomes -C zero, then the element e will be absorbed in scan 2. -C ------------------------------------------------------------- - - DO 150 PME = PME1, PME2 - I = IW (PME) - ELN = ELEN (I) - IF (ELN .GT. 0) THEN -C note that nv (i) has been negated to denote i in Lme: - NVI = -NV (I) - WNVI = WFLG - NVI - DO 140 P = PE (I), PE (I) + ELN - 1 - E = IW (P) - WE = W (E) - IF (WE .GE. WFLG) THEN -C unabsorbed element e has been seen in this loop - WE = WE - NVI - ELSE IF (WE .NE. 0) THEN -C e is an unabsorbed element -C this is the first we have seen e in all of Scan 1 - WE = DEGREE (E) + WNVI - ENDIF - W (E) = WE -140 CONTINUE - ENDIF -150 CONTINUE - -C======================================================================= -C DEGREE UPDATE AND ELEMENT ABSORPTION -C======================================================================= - -C ------------------------------------------------------------- -C Scan 2: for each i in Lme, sum up the degree of Lme (which -C is degme), plus the sum of the external degrees of each Le -C for the elements e appearing within i, plus the -C supervariables in i. Place i in hash list. -C ------------------------------------------------------------- - - DO 180 PME = PME1, PME2 - I = IW (PME) - P1 = PE (I) - P2 = P1 + ELEN (I) - 1 - PN = P1 - HASH = 0 - DEG = 0 - -C ---------------------------------------------------------- -C scan the element list associated with supervariable i -C ---------------------------------------------------------- - -C UMFPACK/MA38-style approximate degree: - DO 160 P = P1, P2 - E = IW (P) - WE = W (E) - IF (WE .NE. 0) THEN -C e is an unabsorbed element - DEG = DEG + WE - WFLG - IW (PN) = E - PN = PN + 1 - HASH = HASH + E - ENDIF -160 CONTINUE - -C count the number of elements in i (including me): - ELEN (I) = PN - P1 + 1 - -C ---------------------------------------------------------- -C scan the supervariables in the list associated with i -C ---------------------------------------------------------- - - P3 = PN - DO 170 P = P2 + 1, P1 + LEN (I) - 1 - J = IW (P) - NVJ = NV (J) - IF (NVJ .GT. 0) THEN -C j is unabsorbed, and not in Lme. -C add to degree and add to new list - DEG = DEG + NVJ - IW (PN) = J - PN = PN + 1 - HASH = HASH + J - ENDIF -170 CONTINUE - -C ---------------------------------------------------------- -C update the degree and check for mass elimination -C ---------------------------------------------------------- - - IF (ELEN (I) .EQ. 1 .AND. P3 .EQ. PN) THEN - -C ------------------------------------------------------- -C mass elimination -C ------------------------------------------------------- - -C There is nothing left of this node except for an -C edge to the current pivot element. elen (i) is 1, -C and there are no variables adjacent to node i. -C Absorb i into the current pivot element, me. - - PE (I) = -ME - NVI = -NV (I) - DEGME = DEGME - NVI - NVPIV = NVPIV + NVI - NEL = NEL + NVI - NV (I) = 0 - ELEN (I) = 0 - - ELSE - -C ------------------------------------------------------- -C update the upper-bound degree of i -C ------------------------------------------------------- - -C the following degree does not yet include the size -C of the current element, which is added later: - DEGREE (I) = MIN (DEGREE (I), DEG) - -C ------------------------------------------------------- -C add me to the list for i -C ------------------------------------------------------- - -C move first supervariable to end of list - IW (PN) = IW (P3) -C move first element to end of element part of list - IW (P3) = IW (P1) -C add new element to front of list. - IW (P1) = ME -C store the new length of the list in len (i) - LEN (I) = PN - P1 + 1 - -C ------------------------------------------------------- -C place in hash bucket. Save hash key of i in last (i). -C ------------------------------------------------------- - - HASH = MOD (HASH, HMOD) + 1 - J = HEAD (HASH) - IF (J .LE. 0) THEN -C the degree list is empty, hash head is -j - NEXT (I) = -J - HEAD (HASH) = -I - ELSE -C degree list is not empty -C use last (head (hash)) as hash head - NEXT (I) = LAST (J) - LAST (J) = I - ENDIF - LAST (I) = HASH - ENDIF -180 CONTINUE - - DEGREE (ME) = DEGME - -C ------------------------------------------------------------- -C Clear the counter array, w (...), by incrementing wflg. -C ------------------------------------------------------------- - - DMAX = MAX (DMAX, DEGME) - WFLG = WFLG + DMAX - -C make sure that wflg+n does not cause integer overflow - IF (WFLG + N .LE. WFLG) THEN - DO 190 X = 1, N - IF (W (X) .NE. 0) W (X) = 1 -190 CONTINUE - WFLG = 2 - ENDIF -C at this point, w (1..n) .lt. wflg holds - -C======================================================================= -C SUPERVARIABLE DETECTION -C======================================================================= - - DO 250 PME = PME1, PME2 - I = IW (PME) - IF (NV (I) .LT. 0) THEN -C i is a principal variable in Lme - -C ------------------------------------------------------- -C examine all hash buckets with 2 or more variables. We -C do this by examing all unique hash keys for super- -C variables in the pattern Lme of the current element, me -C ------------------------------------------------------- - - HASH = LAST (I) -C let i = head of hash bucket, and empty the hash bucket - J = HEAD (HASH) - IF (J .EQ. 0) GOTO 250 - IF (J .LT. 0) THEN -C degree list is empty - I = -J - HEAD (HASH) = 0 - ELSE -C degree list is not empty, restore last () of head - I = LAST (J) - LAST (J) = 0 - ENDIF - IF (I .EQ. 0) GOTO 250 - -C while loop: -200 CONTINUE - IF (NEXT (I) .NE. 0) THEN - -C ---------------------------------------------------- -C this bucket has one or more variables following i. -C scan all of them to see if i can absorb any entries -C that follow i in hash bucket. Scatter i into w. -C ---------------------------------------------------- - - LN = LEN (I) - ELN = ELEN (I) -C do not flag the first element in the list (me) - DO 210 P = PE (I) + 1, PE (I) + LN - 1 - W (IW (P)) = WFLG -210 CONTINUE - -C ---------------------------------------------------- -C scan every other entry j following i in bucket -C ---------------------------------------------------- - - JLAST = I - J = NEXT (I) - -C while loop: -220 CONTINUE - IF (J .NE. 0) THEN - -C ------------------------------------------------- -C check if j and i have identical nonzero pattern -C ------------------------------------------------- - - IF (LEN (J) .NE. LN) THEN -C i and j do not have same size data structure - GOTO 240 - ENDIF - IF (ELEN (J) .NE. ELN) THEN -C i and j do not have same number of adjacent el - GOTO 240 - ENDIF -C do not flag the first element in the list (me) - DO 230 P = PE (J) + 1, PE (J) + LN - 1 - IF (W (IW (P)) .NE. WFLG) THEN -C an entry (iw(p)) is in j but not in i - GOTO 240 - ENDIF -230 CONTINUE - -C ------------------------------------------------- -C found it! j can be absorbed into i -C ------------------------------------------------- - - PE (J) = -I -C both nv (i) and nv (j) are negated since they -C are in Lme, and the absolute values of each -C are the number of variables in i and j: - NV (I) = NV (I) + NV (J) - NV (J) = 0 - ELEN (J) = 0 -C delete j from hash bucket - J = NEXT (J) - NEXT (JLAST) = J - GOTO 220 - -C ------------------------------------------------- -240 CONTINUE -C j cannot be absorbed into i -C ------------------------------------------------- - - JLAST = J - J = NEXT (J) - GOTO 220 - ENDIF - -C ---------------------------------------------------- -C no more variables can be absorbed into i -C go to next i in bucket and clear flag array -C ---------------------------------------------------- - - WFLG = WFLG + 1 - I = NEXT (I) - IF (I .NE. 0) GOTO 200 - ENDIF - ENDIF -250 CONTINUE - -C======================================================================= -C RESTORE DEGREE LISTS AND REMOVE NONPRINCIPAL SUPERVAR. FROM ELEMENT -C======================================================================= - - P = PME1 - NLEFT = N - NEL - DO 260 PME = PME1, PME2 - I = IW (PME) - NVI = -NV (I) - IF (NVI .GT. 0) THEN -C i is a principal variable in Lme -C restore nv (i) to signify that i is principal - NV (I) = NVI - -C ------------------------------------------------------- -C compute the external degree (add size of current elem) -C ------------------------------------------------------- - - DEG = MAX (1, MIN (DEGREE (I) + DEGME-NVI, NLEFT-NVI)) - -C ------------------------------------------------------- -C place the supervariable at the head of the degree list -C ------------------------------------------------------- - - INEXT = HEAD (DEG) - IF (INEXT .NE. 0) LAST (INEXT) = I - NEXT (I) = INEXT - LAST (I) = 0 - HEAD (DEG) = I - -C ------------------------------------------------------- -C save the new degree, and find the minimum degree -C ------------------------------------------------------- - - MINDEG = MIN (MINDEG, DEG) - DEGREE (I) = DEG - -C ------------------------------------------------------- -C place the supervariable in the element pattern -C ------------------------------------------------------- - - IW (P) = I - P = P + 1 - ENDIF -260 CONTINUE - -C======================================================================= -C FINALIZE THE NEW ELEMENT -C======================================================================= - - NV (ME) = NVPIV + DEGME -C nv (me) is now the degree of pivot (including diagonal part) -C save the length of the list for the new element me - LEN (ME) = P - PME1 - IF (LEN (ME) .EQ. 0) THEN -C there is nothing left of the current pivot element - PE (ME) = 0 - W (ME) = 0 - ENDIF - IF (NEWMEM .NE. 0) THEN -C element was not constructed in place: deallocate part -C of it (final size is less than or equal to newmem, -C since newly nonprincipal variables have been removed). - PFREE = P - MEM = MEM - NEWMEM + LEN (ME) - ENDIF - -C======================================================================= -C END WHILE (selecting pivots) - GOTO 30 - ENDIF -C======================================================================= - -C======================================================================= -C COMPUTE THE PERMUTATION VECTORS -C======================================================================= - -C ---------------------------------------------------------------- -C The time taken by the following code is O(n). At this -C point, elen (e) = -k has been done for all elements e, -C and elen (i) = 0 has been done for all nonprincipal -C variables i. At this point, there are no principal -C supervariables left, and all elements are absorbed. -C ---------------------------------------------------------------- - -C ---------------------------------------------------------------- -C compute the ordering of unordered nonprincipal variables -C ---------------------------------------------------------------- - - DO 290 I = 1, N - IF (ELEN (I) .EQ. 0) THEN - -C ---------------------------------------------------------- -C i is an un-ordered row. Traverse the tree from i until -C reaching an element, e. The element, e, was the -C principal supervariable of i and all nodes in the path -C from i to when e was selected as pivot. -C ---------------------------------------------------------- - - J = -PE (I) -C while (j is a variable) do: -270 CONTINUE - IF (ELEN (J) .GE. 0) THEN - J = -PE (J) - GOTO 270 - ENDIF - E = J - -C ---------------------------------------------------------- -C get the current pivot ordering of e -C ---------------------------------------------------------- - - K = -ELEN (E) - -C ---------------------------------------------------------- -C traverse the path again from i to e, and compress the -C path (all nodes point to e). Path compression allows -C this code to compute in O(n) time. Order the unordered -C nodes in the path, and place the element e at the end. -C ---------------------------------------------------------- - - J = I -C while (j is a variable) do: -280 CONTINUE - IF (ELEN (J) .GE. 0) THEN - JNEXT = -PE (J) - PE (J) = -E - IF (ELEN (J) .EQ. 0) THEN -C j is an unordered row - ELEN (J) = K - K = K + 1 - ENDIF - J = JNEXT - GOTO 280 - ENDIF -C leave elen (e) negative, so we know it is an element - ELEN (E) = -K - ENDIF -290 CONTINUE - -C ---------------------------------------------------------------- -C reset the inverse permutation (elen (1..n)) to be positive, -C and compute the permutation (last (1..n)). -C ---------------------------------------------------------------- - - DO 300 I = 1, N - K = ABS (ELEN (I)) - LAST (K) = I - ELEN (I) = K -300 CONTINUE - -C======================================================================= -C RETURN THE MEMORY USAGE IN IW -C======================================================================= - -C If maxmem is less than or equal to iwlen, then no compressions -C occurred, and iw (maxmem+1 ... iwlen) was unused. Otherwise -C compressions did occur, and iwlen would have had to have been -C greater than or equal to maxmem for no compressions to occur. -C Return the value of maxmem in the pfree argument. - - PFREE = MAXMEM - - RETURN - END - diff --git a/cmake/AmiciConfig.cmake.in b/cmake/AmiciConfig.cmake.in index d1e86c151e..d760587766 100644 --- a/cmake/AmiciConfig.cmake.in +++ b/cmake/AmiciConfig.cmake.in @@ -7,23 +7,8 @@ find_package(OpenMP) # Current SUNDIALSConfig.cmake doesn't use KLU's FindKLU, but hardcodes paths, # and is not relocatable. This does not work with Python package installation in # tmpdirs. -list(APPEND CMAKE_MODULE_PATH - @CMAKE_SOURCE_DIR@/ThirdParty/SuiteSparse/lib/cmake/SuiteSparse/) -if(NOT DEFINED SUITESPARSE_CONFIG_ROOT) - set(SUITESPARSE_CONFIG_ROOT @CMAKE_SOURCE_DIR@/ThirdParty/SuiteSparse/) -endif() -if(NOT DEFINED AMD_ROOT) - set(AMD_ROOT @CMAKE_SOURCE_DIR@/ThirdParty/SuiteSparse/) -endif() -if(NOT DEFINED BTF_ROOT) - set(BTF_ROOT @CMAKE_SOURCE_DIR@/ThirdParty/SuiteSparse/) -endif() -if(NOT DEFINED COLAMD_ROOT) - set(COLAMD_ROOT @CMAKE_SOURCE_DIR@/ThirdParty/SuiteSparse/) -endif() -if(NOT DEFINED KLU_ROOT) - set(KLU_ROOT @CMAKE_SOURCE_DIR@/ThirdParty/SuiteSparse/) -endif() +list(APPEND CMAKE_PREFIX_PATH + @CMAKE_SOURCE_DIR@/ThirdParty/SuiteSparse/install/) find_dependency(SuiteSparse_config REQUIRED) find_dependency(AMD REQUIRED) find_dependency(BTF REQUIRED) diff --git a/include/amici/model.h b/include/amici/model.h index 1c020400ac..533f0139dc 100644 --- a/include/amici/model.h +++ b/include/amici/model.h @@ -112,10 +112,6 @@ class Model : public AbstractModel, public ModelDimensions { * @param o2mode Second order sensitivity mode * @param idlist Indexes indicating algebraic components (DAE only) * @param z2event Mapping of event outputs to events - * @param pythonGenerated Flag indicating matlab or python wrapping - * @param ndxdotdp_explicit Number of nonzero elements in `dxdotdp_explicit` - * @param ndxdotdx_explicit Number of nonzero elements in `dxdotdx_explicit` - * @param w_recursion_depth Recursion depth of fw * @param state_independent_events Map of events with state-independent * triggers functions, mapping trigger timepoints to event indices. */ @@ -123,9 +119,7 @@ class Model : public AbstractModel, public ModelDimensions { ModelDimensions const& model_dimensions, SimulationParameters simulation_parameters, amici::SecondOrderMode o2mode, std::vector idlist, - std::vector z2event, bool pythonGenerated = false, - int ndxdotdp_explicit = 0, int ndxdotdx_explicit = 0, - int w_recursion_depth = 0, + std::vector z2event, std::map> state_independent_events = {} ); @@ -1427,19 +1421,20 @@ class Model : public AbstractModel, public ModelDimensions { * @param x_solver State variables with conservation laws applied * (solver returns this) */ - void fx_rdata(AmiVector& x_rdata, AmiVector const& x_solver); + void fx_rdata(gsl::span x_rdata, AmiVector const& x_solver); /** * @brief Expand conservation law for state sensitivities. * @param sx_rdata Output buffer for state variables sensitivities with - * conservation laws expanded (stored in `amici::ReturnData`). + * conservation laws expanded + * (stored in `amici::ReturnData` shape `nplist` x `nx`, row-major). * @param sx_solver State variables sensitivities with conservation laws * applied (solver returns this) * @param x_solver State variables with conservation laws * applied (solver returns this) */ void fsx_rdata( - AmiVectorArray& sx_rdata, AmiVectorArray const& sx_solver, + gsl::span sx_rdata, AmiVectorArray const& sx_solver, AmiVector const& x_solver ); @@ -1457,9 +1452,6 @@ class Model : public AbstractModel, public ModelDimensions { */ std::vector const& getReinitializationStateIdxs() const; - /** Flag indicating Matlab- or Python-based model generation */ - bool pythonGenerated = false; - /** * @brief getter for dxdotdp (matlab generated) * @return dxdotdp @@ -1488,17 +1480,7 @@ class Model : public AbstractModel, public ModelDimensions { * * @return Steady-state mask */ - std::vector get_steadystate_mask() const { - return steadystate_mask_.getVector(); - }; - - /** - * @brief Get steady-state mask as AmiVector. - * - * See `set_steadystate_mask` for details. - * @return Steady-state mask - */ - AmiVector const& get_steadystate_mask_av() const { + std::vector get_steadystate_mask() const { return steadystate_mask_; }; @@ -1513,7 +1495,7 @@ class Model : public AbstractModel, public ModelDimensions { * * @param mask Mask of length `nx_solver`. */ - void set_steadystate_mask(std::vector const& mask); + void set_steadystate_mask(std::vector const& mask); /** * Flag indicating whether for @@ -2107,18 +2089,6 @@ class Model : public AbstractModel, public ModelDimensions { realtype min_sigma_{50.0}; private: - /** Sparse dwdp implicit temporary storage (shape `ndwdp`) */ - mutable std::vector dwdp_hierarchical_; - - /** Sparse dwdw temporary storage (shape `ndwdw`) */ - mutable SUNMatrixWrapper dwdw_; - - /** Sparse dwdx implicit temporary storage (shape `ndwdx`) */ - mutable std::vector dwdx_hierarchical_; - - /** Recursion */ - int w_recursion_depth_{0}; - /** Simulation parameters, initial state, etc. */ SimulationParameters simulation_parameters_; @@ -2129,7 +2099,7 @@ class Model : public AbstractModel, public ModelDimensions { * Negative values indicate that the corresponding state variable should * be ignored. */ - AmiVector steadystate_mask_; + std::vector steadystate_mask_; }; bool operator==(Model const& a, Model const& b); diff --git a/include/amici/model_dae.h b/include/amici/model_dae.h index 03ffbf8fb9..a704fcaef2 100644 --- a/include/amici/model_dae.h +++ b/include/amici/model_dae.h @@ -36,10 +36,6 @@ class Model_DAE : public Model { * @param o2mode second order sensitivity mode * @param idlist indexes indicating algebraic components (DAE only) * @param z2event mapping of event outputs to events - * @param pythonGenerated flag indicating matlab or python wrapping - * @param ndxdotdp_explicit number of nonzero elements dxdotdp_explicit - * @param ndxdotdx_explicit number of nonzero elements dxdotdx_explicit - * @param w_recursion_depth Recursion depth of fw * @param state_independent_events Map of events with state-independent * triggers functions, mapping trigger timepoints to event indices. */ @@ -47,15 +43,12 @@ class Model_DAE : public Model { ModelDimensions const& model_dimensions, SimulationParameters simulation_parameters, SecondOrderMode const o2mode, std::vector const& idlist, - std::vector const& z2event, bool const pythonGenerated = false, - int const ndxdotdp_explicit = 0, int const ndxdotdx_explicit = 0, - int const w_recursion_depth = 0, + std::vector const& z2event, std::map> state_independent_events = {} ) : Model( model_dimensions, simulation_parameters, o2mode, idlist, z2event, - pythonGenerated, ndxdotdp_explicit, ndxdotdx_explicit, - w_recursion_depth, state_independent_events + state_independent_events ) { derived_state_.M_ = SUNMatrixWrapper(nx_solver, nx_solver); auto M_nnz = static_cast( diff --git a/include/amici/model_dimensions.h b/include/amici/model_dimensions.h index b5aa1ba21e..7b6ade81c5 100644 --- a/include/amici/model_dimensions.h +++ b/include/amici/model_dimensions.h @@ -9,7 +9,7 @@ namespace amici { /** * @brief Container for model dimensions. * - * Holds number of states, observables, etc. + * Holds number of state variables, observables, etc. */ struct ModelDimensions { /** Default ctor */ @@ -54,6 +54,11 @@ struct ModelDimensions { * @param nnz Number of nonzero elements in Jacobian * @param ubw Upper matrix bandwidth in the Jacobian * @param lbw Lower matrix bandwidth in the Jacobian + * @param pythonGenerated Flag indicating model creation from Matlab or + * Python + * @param ndxdotdp_explicit Number of nonzero elements in `dxdotdp_explicit` + * @param ndxdotdx_explicit Number of nonzero elements in `dxdotdx_explicit` + * @param w_recursion_depth Recursion depth of fw */ ModelDimensions( int const nx_rdata, int const nxtrue_rdata, int const nx_solver, @@ -64,7 +69,8 @@ struct ModelDimensions { int const ndwdw, int const ndxdotdw, std::vector ndJydy, int const ndxrdatadxsolver, int const ndxrdatadtcl, int const ndtotal_cldx_rdata, int const nnz, int const ubw, - int const lbw + int const lbw, bool pythonGenerated = false, int ndxdotdp_explicit = 0, + int ndxdotdx_explicit = 0, int w_recursion_depth = 0 ) : nx_rdata(nx_rdata) , nxtrue_rdata(nxtrue_rdata) @@ -92,7 +98,11 @@ struct ModelDimensions { , nnz(nnz) , nJ(nJ) , ubw(ubw) - , lbw(lbw) { + , lbw(lbw) + , pythonGenerated(pythonGenerated) + , ndxdotdp_explicit(ndxdotdp_explicit) + , ndxdotdx_explicit(ndxdotdx_explicit) + , w_recursion_depth(w_recursion_depth) { Expects(nxtrue_rdata >= 0); Expects(nxtrue_rdata <= nx_rdata); Expects(nxtrue_solver >= 0); @@ -128,6 +138,9 @@ struct ModelDimensions { Expects(nJ >= 0); Expects(ubw >= 0); Expects(lbw >= 0); + Expects(ndxdotdp_explicit >= 0); + Expects(ndxdotdx_explicit >= 0); + Expects(w_recursion_depth >= 0); } /** Number of states */ @@ -229,6 +242,18 @@ struct ModelDimensions { /** Lower bandwidth of the Jacobian */ int lbw{0}; + + /** Flag indicating model creation from Matlab or Python */ + bool pythonGenerated = false; + + /** Number of nonzero elements in `dxdotdx_explicit` */ + int ndxdotdp_explicit = 0; + + /** Number of nonzero elements in `dxdotdp_explicit` */ + int ndxdotdx_explicit = 0; + + /** Recursion depth of fw */ + int w_recursion_depth = 0; }; } // namespace amici diff --git a/include/amici/model_ode.h b/include/amici/model_ode.h index 24d410a33c..8386eb929e 100644 --- a/include/amici/model_ode.h +++ b/include/amici/model_ode.h @@ -35,10 +35,6 @@ class Model_ODE : public Model { * @param o2mode second order sensitivity mode * @param idlist indexes indicating algebraic components (DAE only) * @param z2event mapping of event outputs to events - * @param pythonGenerated flag indicating matlab or python wrapping - * @param ndxdotdp_explicit number of nonzero elements dxdotdp_explicit - * @param ndxdotdx_explicit number of nonzero elements dxdotdx_explicit - * @param w_recursion_depth Recursion depth of fw * @param state_independent_events Map of events with state-independent * triggers functions, mapping trigger timepoints to event indices. */ @@ -46,15 +42,12 @@ class Model_ODE : public Model { ModelDimensions const& model_dimensions, SimulationParameters simulation_parameters, SecondOrderMode const o2mode, std::vector const& idlist, - std::vector const& z2event, bool const pythonGenerated = false, - int const ndxdotdp_explicit = 0, int const ndxdotdx_explicit = 0, - int const w_recursion_depth = 0, + std::vector const& z2event, std::map> state_independent_events = {} ) : Model( model_dimensions, simulation_parameters, o2mode, idlist, z2event, - pythonGenerated, ndxdotdp_explicit, ndxdotdx_explicit, - w_recursion_depth, state_independent_events + state_independent_events ) {} void diff --git a/include/amici/model_state.h b/include/amici/model_state.h index c6c46df3de..6336f8735d 100644 --- a/include/amici/model_state.h +++ b/include/amici/model_state.h @@ -43,9 +43,6 @@ struct ModelState { * (dimension: nplist) */ std::vector plist; - - /** temporary storage for spline values */ - std::vector spl_; }; inline bool operator==(ModelState const& a, ModelState const& b) { @@ -323,6 +320,21 @@ struct ModelStateDerived { /** temporary storage of positified state variables according to * stateIsNonNegative (dimension: `nx_solver`) */ AmiVector x_pos_tmp_{0}; + + /** temporary storage for spline values */ + std::vector spl_; + + /** Sparse dwdp implicit temporary storage (shape `ndwdp`) */ + std::vector dwdp_hierarchical_; + + /** Sparse dwdw temporary storage (shape `ndwdw`) */ + SUNMatrixWrapper dwdw_; + + /** Sparse dwdx implicit temporary storage (shape `ndwdx`) */ + std::vector dwdx_hierarchical_; + + /** Temporary storage for dense dJydy (dimension: `nJ` x `ny`) */ + SUNMatrixWrapper dJydy_dense_; }; /** @@ -332,11 +344,19 @@ struct ModelStateDerived { struct SimulationState { /** timepoint */ realtype t; - /** state variables */ + /** + * partial state vector, excluding states eliminated from conservation laws + */ AmiVector x; - /** state variables */ + /** + * partial time derivative of state vector, excluding states eliminated + * from conservation laws + */ AmiVector dx; - /** state variable sensitivity */ + /** + * partial sensitivity state vector array, excluding states eliminated from + * conservation laws + */ AmiVectorArray sx; /** state of the model that was used for simulation */ ModelState state; diff --git a/include/amici/newton_solver.h b/include/amici/newton_solver.h index 2e8b2f6573..28de4a5f4d 100644 --- a/include/amici/newton_solver.h +++ b/include/amici/newton_solver.h @@ -2,11 +2,8 @@ #define amici_newton_solver_h #include "amici/solver.h" -#include "amici/sundials_matrix_wrapper.h" #include "amici/vector.h" -#include - namespace amici { class Model; @@ -26,26 +23,21 @@ class NewtonSolver { * @brief Initializes solver according to the dimensions in the provided * model * - * @param model pointer to the model object + * @param model the model object + * @param linsol_type type of linear solver to use */ - explicit NewtonSolver(Model const& model); + explicit NewtonSolver(Model const& model, LinearSolver linsol_type); - /** - * @brief Factory method to create a NewtonSolver based on linsolType - * - * @param simulationSolver solver with settings - * @param model pointer to the model instance - * @return solver NewtonSolver according to the specified linsolType - */ - static std::unique_ptr - getSolver(Solver const& simulationSolver, Model const& model); + NewtonSolver(NewtonSolver const&) = delete; + + NewtonSolver& operator=(NewtonSolver const& other) = delete; /** * @brief Computes the solution of one Newton iteration * * @param delta containing the RHS of the linear system, will be * overwritten by solution to the linear system - * @param model pointer to the model instance + * @param model the model instance * @param state current simulation state */ void getStep(AmiVector& delta, Model& model, SimulationState const& state); @@ -53,8 +45,8 @@ class NewtonSolver { /** * @brief Computes steady state sensitivities * - * @param sx pointer to state variable sensitivities - * @param model pointer to the model instance + * @param sx state variable sensitivities + * @param model the model instance * @param state current simulation state */ void computeNewtonSensis( @@ -65,22 +57,19 @@ class NewtonSolver { * @brief Writes the Jacobian for the Newton iteration and passes it to the * linear solver * - * @param model pointer to the model instance + * @param model the model instance * @param state current simulation state */ - virtual void prepareLinearSystem(Model& model, SimulationState const& state) - = 0; + void prepareLinearSystem(Model& model, SimulationState const& state); /** * Writes the Jacobian (JB) for the Newton iteration and passes it to the * linear solver * - * @param model pointer to the model instance + * @param model the model instance * @param state current simulation state */ - virtual void - prepareLinearSystemB(Model& model, SimulationState const& state) - = 0; + void prepareLinearSystemB(Model& model, SimulationState const& state); /** * @brief Solves the linear system for the Newton step @@ -88,28 +77,25 @@ class NewtonSolver { * @param rhs containing the RHS of the linear system, will be * overwritten by solution to the linear system */ - virtual void solveLinearSystem(AmiVector& rhs) = 0; + void solveLinearSystem(AmiVector& rhs); /** * @brief Reinitialize the linear solver * */ - virtual void reinitialize() = 0; + void reinitialize(); /** - * @brief Checks whether linear system is singular + * @brief Checks whether the linear system is singular * - * @param model pointer to the model instance + * @param model the model instance * @param state current simulation state * @return boolean indicating whether the linear system is singular * (condition number < 1/machine precision) */ - virtual bool is_singular(Model& model, SimulationState const& state) const - = 0; - - virtual ~NewtonSolver() = default; + bool is_singular(Model& model, SimulationState const& state) const; - protected: + private: /** dummy rhs, used as dummy argument when computing J and JB */ AmiVector xdot_; /** dummy state, attached to linear solver */ @@ -119,88 +105,9 @@ class NewtonSolver { /** dummy differential adjoint state, used as dummy argument when computing * JB */ AmiVector dxB_; -}; - -/** - * @brief The NewtonSolverDense provides access to the dense linear solver for - * the Newton method. - */ - -class NewtonSolverDense : public NewtonSolver { - - public: - /** - * @brief constructor for sparse solver - * - * @param model model instance that provides problem dimensions - */ - explicit NewtonSolverDense(Model const& model); - - NewtonSolverDense(NewtonSolverDense const&) = delete; - - NewtonSolverDense& operator=(NewtonSolverDense const& other) = delete; - - ~NewtonSolverDense() override; - - void solveLinearSystem(AmiVector& rhs) override; - - void - prepareLinearSystem(Model& model, SimulationState const& state) override; - - void - prepareLinearSystemB(Model& model, SimulationState const& state) override; - - void reinitialize() override; - - bool is_singular(Model& model, SimulationState const& state) const override; - - private: - /** temporary storage of Jacobian */ - SUNMatrixWrapper Jtmp_; - - /** dense linear solver */ - SUNLinearSolver linsol_{nullptr}; -}; - -/** - * @brief The NewtonSolverSparse provides access to the sparse linear solver for - * the Newton method. - */ - -class NewtonSolverSparse : public NewtonSolver { - - public: - /** - * @brief constructor for dense solver - * - * @param model model instance that provides problem dimensions - */ - explicit NewtonSolverSparse(Model const& model); - - NewtonSolverSparse(NewtonSolverSparse const&) = delete; - - NewtonSolverSparse& operator=(NewtonSolverSparse const& other) = delete; - - ~NewtonSolverSparse() override; - - void solveLinearSystem(AmiVector& rhs) override; - - void - prepareLinearSystem(Model& model, SimulationState const& state) override; - - void - prepareLinearSystemB(Model& model, SimulationState const& state) override; - - bool is_singular(Model& model, SimulationState const& state) const override; - - void reinitialize() override; - - private: - /** temporary storage of Jacobian */ - SUNMatrixWrapper Jtmp_; - /** sparse linear solver */ - SUNLinearSolver linsol_{nullptr}; + /** linear solver */ + std::unique_ptr linsol_; }; } // namespace amici diff --git a/include/amici/rdata.h b/include/amici/rdata.h index b5b4c269b4..44f45ceef6 100644 --- a/include/amici/rdata.h +++ b/include/amici/rdata.h @@ -104,12 +104,12 @@ class ReturnData : public ModelDimensions { */ std::vector ts; - /** time derivative (shape `nx`) evaluated at `t_last`. */ + /** time derivative (shape `nx_solver`) evaluated at `t_last`. */ std::vector xdot; /** - * Jacobian of differential equation right hand side (shape `nx` x `nx`, - * row-major) evaluated at `t_last`. + * Jacobian of differential equation right hand side (shape `nx_solver` x + * `nx_solver`, row-major) evaluated at `t_last`. */ std::vector J; @@ -156,11 +156,12 @@ class ReturnData : public ModelDimensions { */ std::vector s2rz; - /** state (shape `nt` x `nx`, row-major) */ + /** state (shape `nt` x `nx_rdata`, row-major) */ std::vector x; /** - * parameter derivative of state (shape `nt` x `nplist` x `nx`, row-major) + * parameter derivative of state (shape `nt` x `nplist` x `nx_rdata`, + * row-major) */ std::vector sx; @@ -358,18 +359,18 @@ class ReturnData : public ModelDimensions { */ realtype posteq_wrms = NAN; - /** initial state (shape `nx`) */ + /** initial state (shape `nx_rdata`) */ std::vector x0; - /** preequilibration steady state (shape `nx`) */ + /** preequilibration steady state (shape `nx_rdata`) */ std::vector x_ss; - /** initial sensitivities (shape `nplist` x `nx`, row-major) */ + /** initial sensitivities (shape `nplist` x `nx_rdata`, row-major) */ std::vector sx0; /** * preequilibration sensitivities - * (shape `nplist` x `nx`, row-major) + * (shape `nplist` x `nx_rdata`, row-major) */ std::vector sx_ss; @@ -463,28 +464,6 @@ class ReturnData : public ModelDimensions { /** offset for sigma_residuals */ realtype sigma_offset; - /** timepoint for model evaluation*/ - realtype t_; - - /** partial state vector, excluding states eliminated from conservation laws - */ - AmiVector x_solver_; - - /** partial time derivative of state vector, excluding states eliminated - * from conservation laws */ - AmiVector dx_solver_; - - /** partial sensitivity state vector array, excluding states eliminated from - * conservation laws */ - AmiVectorArray sx_solver_; - - /** full state vector, including states eliminated from conservation laws */ - AmiVector x_rdata_; - - /** full sensitivity state vector array, including states eliminated from - * conservation laws */ - AmiVectorArray sx_rdata_; - /** array of number of found roots for a certain event type * (shape `ne`) */ std::vector nroots_; @@ -568,18 +547,25 @@ class ReturnData : public ModelDimensions { template void storeJacobianAndDerivativeInReturnData(T const& problem, Model& model) { - readSimulationState(problem.getFinalSimulationState(), model); + auto simulation_state = problem.getFinalSimulationState(); + model.setModelState(simulation_state.state); AmiVector xdot(nx_solver); if (!this->xdot.empty() || !this->J.empty()) - model.fxdot(t_, x_solver_, dx_solver_, xdot); + model.fxdot( + simulation_state.t, simulation_state.x, simulation_state.dx, + xdot + ); if (!this->xdot.empty()) writeSlice(xdot, this->xdot); if (!this->J.empty()) { SUNMatrixWrapper J(nx_solver, nx_solver); - model.fJ(t_, 0.0, x_solver_, dx_solver_, xdot, J); + model.fJ( + simulation_state.t, 0.0, simulation_state.x, + simulation_state.dx, xdot, J + ); // CVODES uses colmajor, so we need to transform to rowmajor for (int ix = 0; ix < model.nx_solver; ix++) for (int jx = 0; jx < model.nx_solver; jx++) @@ -587,21 +573,18 @@ class ReturnData : public ModelDimensions { = J.data()[ix + model.nx_solver * jx]; } } - /** - * @brief sets member variables and model state according to provided - * simulation state - * @param state simulation state provided by Problem - * @param model model that was used for forward/backward simulation - */ - void readSimulationState(SimulationState const& state, Model& model); /** * @brief Residual function * @param it time index * @param model model that was used for forward/backward simulation + * @param simulation_state simulation state the timepoint `it` * @param edata ExpData instance containing observable data */ - void fres(int it, Model& model, ExpData const& edata); + void fres( + int it, Model& model, SimulationState const& simulation_state, + ExpData const& edata + ); /** * @brief Chi-squared function @@ -614,17 +597,25 @@ class ReturnData : public ModelDimensions { * @brief Residual sensitivity function * @param it time index * @param model model that was used for forward/backward simulation + * @param simulation_state simulation state the timepoint `it` * @param edata ExpData instance containing observable data */ - void fsres(int it, Model& model, ExpData const& edata); + void fsres( + int it, Model& model, SimulationState const& simulation_state, + ExpData const& edata + ); /** * @brief Fisher information matrix function * @param it time index * @param model model that was used for forward/backward simulation + * @param simulation_state simulation state the timepoint `it` * @param edata ExpData instance containing observable data */ - void fFIM(int it, Model& model, ExpData const& edata); + void fFIM( + int it, Model& model, SimulationState const& simulation_state, + ExpData const& edata + ); /** * @brief Set likelihood, state variables, outputs and respective @@ -665,46 +656,58 @@ class ReturnData : public ModelDimensions { /** * @brief Extracts output information for data-points, expects that - * x_solver_ and sx_solver_ were set appropriately + * the model state was set appropriately * @param it timepoint index * @param model model that was used in forward solve + * @param simulation_state simulation state the timepoint `it` * @param edata ExpData instance carrying experimental data */ - void getDataOutput(int it, Model& model, ExpData const* edata); + void getDataOutput( + int it, Model& model, SimulationState const& simulation_state, + ExpData const* edata + ); /** * @brief Extracts data information for forward sensitivity analysis, - * expects that x_solver_ and sx_solver_ were set appropriately + * expects that the model state was set appropriately * @param it index of current timepoint * @param model model that was used in forward solve + * @param simulation_state simulation state the timepoint `it` * @param edata ExpData instance carrying experimental data */ - void getDataSensisFSA(int it, Model& model, ExpData const* edata); + void getDataSensisFSA( + int it, Model& model, SimulationState const& simulation_state, + ExpData const* edata + ); /** - * @brief Extracts output information for events, expects that x_solver_ - * and sx_solver_ were set appropriately + * @brief Extracts output information for events, expects that the model + * state was set appropriately * @param t event timepoint * @param rootidx information about which roots fired * (1 indicating fired, 0/-1 for not) * @param model model that was used in forward solve + * @param simulation_state simulation state the timepoint `it` * @param edata ExpData instance carrying experimental data */ void getEventOutput( - realtype t, std::vector const rootidx, Model& model, - ExpData const* edata + realtype t, std::vector const& rootidx, Model& model, + SimulationState const& simulation_state, ExpData const* edata ); /** * @brief Extracts event information for forward sensitivity analysis, - * expects that x_solver_ and sx_solver_ were set appropriately + * expects the model state was set appropriately * @param ie index of event type * @param t event timepoint * @param model model that was used in forward solve + * @param simulation_state simulation state the timepoint `it` * @param edata ExpData instance carrying experimental data */ - void - getEventSensisFSA(int ie, realtype t, Model& model, ExpData const* edata); + void getEventSensisFSA( + int ie, realtype t, Model& model, + SimulationState const& simulation_state, ExpData const* edata + ); /** * @brief Updates contribution to likelihood from quadratures (xQB), @@ -725,12 +728,14 @@ class ReturnData : public ModelDimensions { * (llhS0), if no preequilibration was run or if forward sensitivities were * used * @param model model that was used for forward/backward simulation + * @param simulation_state simulation state the timepoint `it` * @param llhS0 contribution to likelihood for initial state sensitivities * @param xB vector with final adjoint state * (excluding conservation laws) */ void handleSx0Forward( - Model const& model, std::vector& llhS0, AmiVector& xB + Model const& model, SimulationState const& simulation_state, + std::vector& llhS0, AmiVector& xB ) const; }; diff --git a/include/amici/steadystateproblem.h b/include/amici/steadystateproblem.h index 72d248f8bf..f680b6c128 100644 --- a/include/amici/steadystateproblem.h +++ b/include/amici/steadystateproblem.h @@ -400,6 +400,8 @@ class SteadystateProblem { AmiVector xQB_; /** time-derivative of quadrature state vector */ AmiVector xQBdot_; + /** NVector around Model::steadystate_mask_ */ + AmiVector steadystate_mask_; /** maximum number of steps for Newton solver for allocating numlinsteps */ int max_steps_{0}; @@ -450,7 +452,7 @@ class SteadystateProblem { realtype rtol_quad_{NAN}; /** newton solver */ - std::unique_ptr newton_solver_{nullptr}; + NewtonSolver newton_solver_; /** damping factor flag */ NewtonDampingFactorMode damping_factor_mode_{NewtonDampingFactorMode::on}; diff --git a/include/amici/sundials_linsol_wrapper.h b/include/amici/sundials_linsol_wrapper.h index 613eb22156..d9d2105c38 100644 --- a/include/amici/sundials_linsol_wrapper.h +++ b/include/amici/sundials_linsol_wrapper.h @@ -33,10 +33,21 @@ class SUNLinSolWrapper { /** * @brief Wrap existing SUNLinearSolver - * @param linsol + * + * @param linsol SUNLinSolWrapper takes ownership of `linsol`. */ explicit SUNLinSolWrapper(SUNLinearSolver linsol); + /** + * @brief Wrap existing SUNLinearSolver + * + * @param linsol SUNLinSolWrapper takes ownership of `linsol`. + * @param A Matrix + */ + explicit SUNLinSolWrapper( + SUNLinearSolver linsol, SUNMatrixWrapper const& A + ); + virtual ~SUNLinSolWrapper(); /** @@ -80,26 +91,17 @@ class SUNLinSolWrapper { /** * @brief Performs any linear solver setup needed, based on an updated * system matrix A. - * @param A */ - void setup(SUNMatrix A) const; - - /** - * @brief Performs any linear solver setup needed, based on an updated - * system matrix A. - * @param A - */ - void setup(SUNMatrixWrapper const& A) const; + void setup() const; /** * @brief Solves a linear system A*x = b - * @param A * @param x A template for cloning vectors needed within the solver. * @param b * @param tol Tolerance (weighted 2-norm), iterative solvers only * @return error flag */ - int Solve(SUNMatrix A, N_Vector x, N_Vector b, realtype tol) const; + int solve(N_Vector x, N_Vector b, realtype tol) const; /** * @brief Returns the last error flag encountered within the linear solver @@ -119,7 +121,7 @@ class SUNLinSolWrapper { * @brief Get the matrix A (matrix solvers only). * @return A */ - virtual SUNMatrix getMatrix() const; + virtual SUNMatrixWrapper& getMatrix(); protected: /** @@ -131,6 +133,9 @@ class SUNLinSolWrapper { /** Wrapped solver */ SUNLinearSolver solver_{nullptr}; + + /** Matrix A for solver. */ + SUNMatrixWrapper A_; }; /** @@ -139,12 +144,12 @@ class SUNLinSolWrapper { class SUNLinSolBand : public SUNLinSolWrapper { public: /** - * @brief Create solver using existing matrix A without taking ownership of - * A. + * @brief Create solver using existing matrix A + * * @param x A template for cloning vectors needed within the solver. * @param A square matrix */ - SUNLinSolBand(N_Vector x, SUNMatrix A); + SUNLinSolBand(N_Vector x, SUNMatrixWrapper A); /** * @brief Create new band solver and matrix A. @@ -153,12 +158,6 @@ class SUNLinSolBand : public SUNLinSolWrapper { * @param lbw lower bandwidth of band matrix A */ SUNLinSolBand(AmiVector const& x, int ubw, int lbw); - - SUNMatrix getMatrix() const override; - - private: - /** Matrix A for solver, only if created by here. */ - SUNMatrixWrapper A_; }; /** @@ -171,12 +170,6 @@ class SUNLinSolDense : public SUNLinSolWrapper { * @param x A template for cloning vectors needed within the solver. */ explicit SUNLinSolDense(AmiVector const& x); - - SUNMatrix getMatrix() const override; - - private: - /** Matrix A for solver, only if created by here. */ - SUNMatrixWrapper A_; }; /** @@ -192,7 +185,7 @@ class SUNLinSolKLU : public SUNLinSolWrapper { * @param x A template for cloning vectors needed within the solver. * @param A sparse matrix */ - SUNLinSolKLU(N_Vector x, SUNMatrix A); + SUNLinSolKLU(N_Vector x, SUNMatrixWrapper A); /** * @brief Create KLU solver and matrix to operate on @@ -202,11 +195,10 @@ class SUNLinSolKLU : public SUNLinSolWrapper { * @param ordering */ SUNLinSolKLU( - AmiVector const& x, int nnz, int sparsetype, StateOrdering ordering + AmiVector const& x, int nnz, int sparsetype, + StateOrdering ordering = StateOrdering::COLAMD ); - SUNMatrix getMatrix() const override; - /** * @brief Reinitializes memory and flags for a new factorization * (symbolic and numeric) to be conducted at the next solver setup call. @@ -224,9 +216,13 @@ class SUNLinSolKLU : public SUNLinSolWrapper { */ void setOrdering(StateOrdering ordering); - private: - /** Sparse matrix A for solver, only if created by here. */ - SUNMatrixWrapper A_; + /** + * @brief Checks whether the linear system is singular + * + * @return boolean indicating whether the linear system is singular + * (condition number < 1/machine precision) + */ + bool is_singular() const; }; #ifdef SUNDIALS_SUPERLUMT @@ -249,7 +245,7 @@ class SUNLinSolSuperLUMT : public SUNLinSolWrapper { * @param A sparse matrix * @param numThreads Number of threads to be used by SuperLUMT */ - SUNLinSolSuperLUMT(N_Vector x, SUNMatrix A, int numThreads); + SUNLinSolSuperLUMT(N_Vector x, SUNMatrixWrapper A, int numThreads); /** * @brief Create SuperLUMT solver and matrix to operate on @@ -279,18 +275,12 @@ class SUNLinSolSuperLUMT : public SUNLinSolWrapper { int numThreads ); - SUNMatrix getMatrix() const override; - /** * @brief Sets the ordering used by SuperLUMT for reducing fill in the * linear solve. * @param ordering */ void setOrdering(StateOrdering ordering); - - private: - /** Sparse matrix A for solver, only if created by here. */ - SUNMatrixWrapper A; }; #endif diff --git a/models/model_calvetti/CMakeLists.txt b/models/model_calvetti/CMakeLists.txt index 30262f7de3..c685278a4f 100644 --- a/models/model_calvetti/CMakeLists.txt +++ b/models/model_calvetti/CMakeLists.txt @@ -1,11 +1,6 @@ # Build AMICI model -cmake_minimum_required(VERSION 3.15) -cmake_policy(VERSION 3.15...3.27) - -# cmake >=3.27 -if(POLICY CMP0144) - cmake_policy(SET CMP0144 NEW) -endif(POLICY CMP0144) +cmake_minimum_required(VERSION 3.22) +cmake_policy(VERSION 3.22...3.30) project(model_calvetti) @@ -38,6 +33,10 @@ endif() find_package(Amici REQUIRED HINTS ${CMAKE_CURRENT_LIST_DIR}/../../build) message(STATUS "Found AMICI ${Amici_DIR}") +set_target_properties(Upstream::amici PROPERTIES + MAP_IMPORTED_CONFIG_RELWITHDEBINFO RelWithDebInfo;Release; + MAP_IMPORTED_CONFIG_RELEASE Release + MAP_IMPORTED_CONFIG_DEBUG Debug;RelWithDebInfo;) # Debug build? if("$ENV{ENABLE_AMICI_DEBUGGING}" OR "$ENV{ENABLE_GCOV_COVERAGE}") @@ -47,7 +46,6 @@ if("$ENV{ENABLE_AMICI_DEBUGGING}" OR "$ENV{ENABLE_GCOV_COVERAGE}") else() add_compile_options(-O0 -g) endif() - set(CMAKE_BUILD_TYPE "Debug") endif() # coverage options diff --git a/models/model_calvetti/model_calvetti.h b/models/model_calvetti/model_calvetti.h index c189ae8d74..bdc00bdac9 100644 --- a/models/model_calvetti/model_calvetti.h +++ b/models/model_calvetti/model_calvetti.h @@ -1,6 +1,6 @@ #ifndef _amici_model_calvetti_h #define _amici_model_calvetti_h -/* Generated by amiwrap (R2017b) 8b324bca5b796a93094195d22a023e5f8e945ef1 */ +/* Generated by amiwrap (R2017b) d69c026f8f8a89a13f1fa1307548f2b8f7045fe1 */ #include #include #include "amici/defines.h" @@ -72,7 +72,7 @@ class Model_model_calvetti : public amici::Model_DAE { amici::Model* clone() const override { return new Model_model_calvetti(*this); }; - std::string getAmiciCommit() const override { return "8b324bca5b796a93094195d22a023e5f8e945ef1"; }; + std::string getAmiciCommit() const override { return "d69c026f8f8a89a13f1fa1307548f2b8f7045fe1"; }; void fJSparse(SUNMatrixContent_Sparse JSparse, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype cj, const realtype *dx, const realtype *w, const realtype *dwdx) override { JSparse_model_calvetti(JSparse, t, x, p, k, h, cj, dx, w, dwdx); diff --git a/models/model_calvetti/swig/CMakeLists.txt b/models/model_calvetti/swig/CMakeLists.txt index eb9fe681c9..523571c52e 100644 --- a/models/model_calvetti/swig/CMakeLists.txt +++ b/models/model_calvetti/swig/CMakeLists.txt @@ -1,4 +1,15 @@ cmake_minimum_required(VERSION 3.15) +cmake_policy(VERSION 3.15...3.27) + +# cmake >=3.27 +if(POLICY CMP0144) + cmake_policy(SET CMP0144 NEW) +endif(POLICY CMP0144) +# cmake >= 3.30 +if(POLICY CMP0167) + cmake_policy(SET CMP0167 NEW) +endif(POLICY CMP0167) + if(DEFINED ENV{SWIG}) set(SWIG_EXECUTABLE $ENV{SWIG}) diff --git a/models/model_dirac/CMakeLists.txt b/models/model_dirac/CMakeLists.txt index d96169b04a..8c07511c40 100644 --- a/models/model_dirac/CMakeLists.txt +++ b/models/model_dirac/CMakeLists.txt @@ -1,11 +1,6 @@ # Build AMICI model -cmake_minimum_required(VERSION 3.15) -cmake_policy(VERSION 3.15...3.27) - -# cmake >=3.27 -if(POLICY CMP0144) - cmake_policy(SET CMP0144 NEW) -endif(POLICY CMP0144) +cmake_minimum_required(VERSION 3.22) +cmake_policy(VERSION 3.22...3.30) project(model_dirac) @@ -38,6 +33,10 @@ endif() find_package(Amici REQUIRED HINTS ${CMAKE_CURRENT_LIST_DIR}/../../build) message(STATUS "Found AMICI ${Amici_DIR}") +set_target_properties(Upstream::amici PROPERTIES + MAP_IMPORTED_CONFIG_RELWITHDEBINFO RelWithDebInfo;Release; + MAP_IMPORTED_CONFIG_RELEASE Release + MAP_IMPORTED_CONFIG_DEBUG Debug;RelWithDebInfo;) # Debug build? if("$ENV{ENABLE_AMICI_DEBUGGING}" OR "$ENV{ENABLE_GCOV_COVERAGE}") @@ -47,7 +46,6 @@ if("$ENV{ENABLE_AMICI_DEBUGGING}" OR "$ENV{ENABLE_GCOV_COVERAGE}") else() add_compile_options(-O0 -g) endif() - set(CMAKE_BUILD_TYPE "Debug") endif() # coverage options diff --git a/models/model_dirac/model_dirac.h b/models/model_dirac/model_dirac.h index 0425320765..a8803237ca 100644 --- a/models/model_dirac/model_dirac.h +++ b/models/model_dirac/model_dirac.h @@ -1,6 +1,6 @@ #ifndef _amici_model_dirac_h #define _amici_model_dirac_h -/* Generated by amiwrap (R2017b) 8b324bca5b796a93094195d22a023e5f8e945ef1 */ +/* Generated by amiwrap (R2017b) d69c026f8f8a89a13f1fa1307548f2b8f7045fe1 */ #include #include #include "amici/defines.h" @@ -72,7 +72,7 @@ class Model_model_dirac : public amici::Model_ODE { amici::Model* clone() const override { return new Model_model_dirac(*this); }; - std::string getAmiciCommit() const override { return "8b324bca5b796a93094195d22a023e5f8e945ef1"; }; + std::string getAmiciCommit() const override { return "d69c026f8f8a89a13f1fa1307548f2b8f7045fe1"; }; void fJSparse(SUNMatrixContent_Sparse JSparse, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *w, const realtype *dwdx) override { JSparse_model_dirac(JSparse, t, x, p, k, h, w, dwdx); diff --git a/models/model_dirac/swig/CMakeLists.txt b/models/model_dirac/swig/CMakeLists.txt index eb9fe681c9..523571c52e 100644 --- a/models/model_dirac/swig/CMakeLists.txt +++ b/models/model_dirac/swig/CMakeLists.txt @@ -1,4 +1,15 @@ cmake_minimum_required(VERSION 3.15) +cmake_policy(VERSION 3.15...3.27) + +# cmake >=3.27 +if(POLICY CMP0144) + cmake_policy(SET CMP0144 NEW) +endif(POLICY CMP0144) +# cmake >= 3.30 +if(POLICY CMP0167) + cmake_policy(SET CMP0167 NEW) +endif(POLICY CMP0167) + if(DEFINED ENV{SWIG}) set(SWIG_EXECUTABLE $ENV{SWIG}) diff --git a/models/model_events/CMakeLists.txt b/models/model_events/CMakeLists.txt index 53e3a335b8..706de289cf 100644 --- a/models/model_events/CMakeLists.txt +++ b/models/model_events/CMakeLists.txt @@ -1,11 +1,6 @@ # Build AMICI model -cmake_minimum_required(VERSION 3.15) -cmake_policy(VERSION 3.15...3.27) - -# cmake >=3.27 -if(POLICY CMP0144) - cmake_policy(SET CMP0144 NEW) -endif(POLICY CMP0144) +cmake_minimum_required(VERSION 3.22) +cmake_policy(VERSION 3.22...3.30) project(model_events) @@ -38,6 +33,10 @@ endif() find_package(Amici REQUIRED HINTS ${CMAKE_CURRENT_LIST_DIR}/../../build) message(STATUS "Found AMICI ${Amici_DIR}") +set_target_properties(Upstream::amici PROPERTIES + MAP_IMPORTED_CONFIG_RELWITHDEBINFO RelWithDebInfo;Release; + MAP_IMPORTED_CONFIG_RELEASE Release + MAP_IMPORTED_CONFIG_DEBUG Debug;RelWithDebInfo;) # Debug build? if("$ENV{ENABLE_AMICI_DEBUGGING}" OR "$ENV{ENABLE_GCOV_COVERAGE}") @@ -47,7 +46,6 @@ if("$ENV{ENABLE_AMICI_DEBUGGING}" OR "$ENV{ENABLE_GCOV_COVERAGE}") else() add_compile_options(-O0 -g) endif() - set(CMAKE_BUILD_TYPE "Debug") endif() # coverage options diff --git a/models/model_events/model_events.h b/models/model_events/model_events.h index bc2b1b6157..bf68b4529a 100644 --- a/models/model_events/model_events.h +++ b/models/model_events/model_events.h @@ -1,6 +1,6 @@ #ifndef _amici_model_events_h #define _amici_model_events_h -/* Generated by amiwrap (R2017b) 8b324bca5b796a93094195d22a023e5f8e945ef1 */ +/* Generated by amiwrap (R2017b) d69c026f8f8a89a13f1fa1307548f2b8f7045fe1 */ #include #include #include "amici/defines.h" @@ -86,7 +86,7 @@ class Model_model_events : public amici::Model_ODE { amici::Model* clone() const override { return new Model_model_events(*this); }; - std::string getAmiciCommit() const override { return "8b324bca5b796a93094195d22a023e5f8e945ef1"; }; + std::string getAmiciCommit() const override { return "d69c026f8f8a89a13f1fa1307548f2b8f7045fe1"; }; void fJSparse(SUNMatrixContent_Sparse JSparse, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *w, const realtype *dwdx) override { JSparse_model_events(JSparse, t, x, p, k, h, w, dwdx); diff --git a/models/model_events/swig/CMakeLists.txt b/models/model_events/swig/CMakeLists.txt index eb9fe681c9..523571c52e 100644 --- a/models/model_events/swig/CMakeLists.txt +++ b/models/model_events/swig/CMakeLists.txt @@ -1,4 +1,15 @@ cmake_minimum_required(VERSION 3.15) +cmake_policy(VERSION 3.15...3.27) + +# cmake >=3.27 +if(POLICY CMP0144) + cmake_policy(SET CMP0144 NEW) +endif(POLICY CMP0144) +# cmake >= 3.30 +if(POLICY CMP0167) + cmake_policy(SET CMP0167 NEW) +endif(POLICY CMP0167) + if(DEFINED ENV{SWIG}) set(SWIG_EXECUTABLE $ENV{SWIG}) diff --git a/models/model_jakstat_adjoint/CMakeLists.txt b/models/model_jakstat_adjoint/CMakeLists.txt index 75cc527694..4407583796 100644 --- a/models/model_jakstat_adjoint/CMakeLists.txt +++ b/models/model_jakstat_adjoint/CMakeLists.txt @@ -1,11 +1,6 @@ # Build AMICI model -cmake_minimum_required(VERSION 3.15) -cmake_policy(VERSION 3.15...3.27) - -# cmake >=3.27 -if(POLICY CMP0144) - cmake_policy(SET CMP0144 NEW) -endif(POLICY CMP0144) +cmake_minimum_required(VERSION 3.22) +cmake_policy(VERSION 3.22...3.30) project(model_jakstat_adjoint) @@ -38,6 +33,10 @@ endif() find_package(Amici REQUIRED HINTS ${CMAKE_CURRENT_LIST_DIR}/../../build) message(STATUS "Found AMICI ${Amici_DIR}") +set_target_properties(Upstream::amici PROPERTIES + MAP_IMPORTED_CONFIG_RELWITHDEBINFO RelWithDebInfo;Release; + MAP_IMPORTED_CONFIG_RELEASE Release + MAP_IMPORTED_CONFIG_DEBUG Debug;RelWithDebInfo;) # Debug build? if("$ENV{ENABLE_AMICI_DEBUGGING}" OR "$ENV{ENABLE_GCOV_COVERAGE}") @@ -47,7 +46,6 @@ if("$ENV{ENABLE_AMICI_DEBUGGING}" OR "$ENV{ENABLE_GCOV_COVERAGE}") else() add_compile_options(-O0 -g) endif() - set(CMAKE_BUILD_TYPE "Debug") endif() # coverage options diff --git a/models/model_jakstat_adjoint/model_jakstat_adjoint.h b/models/model_jakstat_adjoint/model_jakstat_adjoint.h index aad482c26c..7d3b9290f4 100644 --- a/models/model_jakstat_adjoint/model_jakstat_adjoint.h +++ b/models/model_jakstat_adjoint/model_jakstat_adjoint.h @@ -1,6 +1,6 @@ #ifndef _amici_model_jakstat_adjoint_h #define _amici_model_jakstat_adjoint_h -/* Generated by amiwrap (R2017b) 8b324bca5b796a93094195d22a023e5f8e945ef1 */ +/* Generated by amiwrap (R2017b) d69c026f8f8a89a13f1fa1307548f2b8f7045fe1 */ #include #include #include "amici/defines.h" @@ -75,7 +75,7 @@ class Model_model_jakstat_adjoint : public amici::Model_ODE { amici::Model* clone() const override { return new Model_model_jakstat_adjoint(*this); }; - std::string getAmiciCommit() const override { return "8b324bca5b796a93094195d22a023e5f8e945ef1"; }; + std::string getAmiciCommit() const override { return "d69c026f8f8a89a13f1fa1307548f2b8f7045fe1"; }; void fJSparse(SUNMatrixContent_Sparse JSparse, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *w, const realtype *dwdx) override { JSparse_model_jakstat_adjoint(JSparse, t, x, p, k, h, w, dwdx); diff --git a/models/model_jakstat_adjoint/swig/CMakeLists.txt b/models/model_jakstat_adjoint/swig/CMakeLists.txt index eb9fe681c9..523571c52e 100644 --- a/models/model_jakstat_adjoint/swig/CMakeLists.txt +++ b/models/model_jakstat_adjoint/swig/CMakeLists.txt @@ -1,4 +1,15 @@ cmake_minimum_required(VERSION 3.15) +cmake_policy(VERSION 3.15...3.27) + +# cmake >=3.27 +if(POLICY CMP0144) + cmake_policy(SET CMP0144 NEW) +endif(POLICY CMP0144) +# cmake >= 3.30 +if(POLICY CMP0167) + cmake_policy(SET CMP0167 NEW) +endif(POLICY CMP0167) + if(DEFINED ENV{SWIG}) set(SWIG_EXECUTABLE $ENV{SWIG}) diff --git a/models/model_jakstat_adjoint_o2/CMakeLists.txt b/models/model_jakstat_adjoint_o2/CMakeLists.txt index 4b2b35e223..a8406acf6b 100644 --- a/models/model_jakstat_adjoint_o2/CMakeLists.txt +++ b/models/model_jakstat_adjoint_o2/CMakeLists.txt @@ -1,11 +1,6 @@ # Build AMICI model -cmake_minimum_required(VERSION 3.15) -cmake_policy(VERSION 3.15...3.27) - -# cmake >=3.27 -if(POLICY CMP0144) - cmake_policy(SET CMP0144 NEW) -endif(POLICY CMP0144) +cmake_minimum_required(VERSION 3.22) +cmake_policy(VERSION 3.22...3.30) project(model_jakstat_adjoint_o2) @@ -38,6 +33,10 @@ endif() find_package(Amici REQUIRED HINTS ${CMAKE_CURRENT_LIST_DIR}/../../build) message(STATUS "Found AMICI ${Amici_DIR}") +set_target_properties(Upstream::amici PROPERTIES + MAP_IMPORTED_CONFIG_RELWITHDEBINFO RelWithDebInfo;Release; + MAP_IMPORTED_CONFIG_RELEASE Release + MAP_IMPORTED_CONFIG_DEBUG Debug;RelWithDebInfo;) # Debug build? if("$ENV{ENABLE_AMICI_DEBUGGING}" OR "$ENV{ENABLE_GCOV_COVERAGE}") @@ -47,7 +46,6 @@ if("$ENV{ENABLE_AMICI_DEBUGGING}" OR "$ENV{ENABLE_GCOV_COVERAGE}") else() add_compile_options(-O0 -g) endif() - set(CMAKE_BUILD_TYPE "Debug") endif() # coverage options diff --git a/models/model_jakstat_adjoint_o2/model_jakstat_adjoint_o2.h b/models/model_jakstat_adjoint_o2/model_jakstat_adjoint_o2.h index e44c31b8d5..d7b34417b6 100644 --- a/models/model_jakstat_adjoint_o2/model_jakstat_adjoint_o2.h +++ b/models/model_jakstat_adjoint_o2/model_jakstat_adjoint_o2.h @@ -1,6 +1,6 @@ #ifndef _amici_model_jakstat_adjoint_o2_h #define _amici_model_jakstat_adjoint_o2_h -/* Generated by amiwrap (R2017b) 8b324bca5b796a93094195d22a023e5f8e945ef1 */ +/* Generated by amiwrap (R2017b) d69c026f8f8a89a13f1fa1307548f2b8f7045fe1 */ #include #include #include "amici/defines.h" @@ -75,7 +75,7 @@ class Model_model_jakstat_adjoint_o2 : public amici::Model_ODE { amici::Model* clone() const override { return new Model_model_jakstat_adjoint_o2(*this); }; - std::string getAmiciCommit() const override { return "8b324bca5b796a93094195d22a023e5f8e945ef1"; }; + std::string getAmiciCommit() const override { return "d69c026f8f8a89a13f1fa1307548f2b8f7045fe1"; }; void fJSparse(SUNMatrixContent_Sparse JSparse, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *w, const realtype *dwdx) override { JSparse_model_jakstat_adjoint_o2(JSparse, t, x, p, k, h, w, dwdx); diff --git a/models/model_jakstat_adjoint_o2/swig/CMakeLists.txt b/models/model_jakstat_adjoint_o2/swig/CMakeLists.txt index eb9fe681c9..523571c52e 100644 --- a/models/model_jakstat_adjoint_o2/swig/CMakeLists.txt +++ b/models/model_jakstat_adjoint_o2/swig/CMakeLists.txt @@ -1,4 +1,15 @@ cmake_minimum_required(VERSION 3.15) +cmake_policy(VERSION 3.15...3.27) + +# cmake >=3.27 +if(POLICY CMP0144) + cmake_policy(SET CMP0144 NEW) +endif(POLICY CMP0144) +# cmake >= 3.30 +if(POLICY CMP0167) + cmake_policy(SET CMP0167 NEW) +endif(POLICY CMP0167) + if(DEFINED ENV{SWIG}) set(SWIG_EXECUTABLE $ENV{SWIG}) diff --git a/models/model_nested_events/CMakeLists.txt b/models/model_nested_events/CMakeLists.txt index c609531e1d..0dc2b916bb 100644 --- a/models/model_nested_events/CMakeLists.txt +++ b/models/model_nested_events/CMakeLists.txt @@ -1,11 +1,6 @@ # Build AMICI model -cmake_minimum_required(VERSION 3.15) -cmake_policy(VERSION 3.15...3.27) - -# cmake >=3.27 -if(POLICY CMP0144) - cmake_policy(SET CMP0144 NEW) -endif(POLICY CMP0144) +cmake_minimum_required(VERSION 3.22) +cmake_policy(VERSION 3.22...3.30) project(model_nested_events) @@ -38,6 +33,10 @@ endif() find_package(Amici REQUIRED HINTS ${CMAKE_CURRENT_LIST_DIR}/../../build) message(STATUS "Found AMICI ${Amici_DIR}") +set_target_properties(Upstream::amici PROPERTIES + MAP_IMPORTED_CONFIG_RELWITHDEBINFO RelWithDebInfo;Release; + MAP_IMPORTED_CONFIG_RELEASE Release + MAP_IMPORTED_CONFIG_DEBUG Debug;RelWithDebInfo;) # Debug build? if("$ENV{ENABLE_AMICI_DEBUGGING}" OR "$ENV{ENABLE_GCOV_COVERAGE}") @@ -47,7 +46,6 @@ if("$ENV{ENABLE_AMICI_DEBUGGING}" OR "$ENV{ENABLE_GCOV_COVERAGE}") else() add_compile_options(-O0 -g) endif() - set(CMAKE_BUILD_TYPE "Debug") endif() # coverage options diff --git a/models/model_nested_events/model_nested_events.h b/models/model_nested_events/model_nested_events.h index 5cad2049a9..0ca111adc7 100644 --- a/models/model_nested_events/model_nested_events.h +++ b/models/model_nested_events/model_nested_events.h @@ -1,6 +1,6 @@ #ifndef _amici_model_nested_events_h #define _amici_model_nested_events_h -/* Generated by amiwrap (R2017b) 8b324bca5b796a93094195d22a023e5f8e945ef1 */ +/* Generated by amiwrap (R2017b) d69c026f8f8a89a13f1fa1307548f2b8f7045fe1 */ #include #include #include "amici/defines.h" @@ -75,7 +75,7 @@ class Model_model_nested_events : public amici::Model_ODE { amici::Model* clone() const override { return new Model_model_nested_events(*this); }; - std::string getAmiciCommit() const override { return "8b324bca5b796a93094195d22a023e5f8e945ef1"; }; + std::string getAmiciCommit() const override { return "d69c026f8f8a89a13f1fa1307548f2b8f7045fe1"; }; void fJSparse(SUNMatrixContent_Sparse JSparse, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *w, const realtype *dwdx) override { JSparse_model_nested_events(JSparse, t, x, p, k, h, w, dwdx); diff --git a/models/model_nested_events/swig/CMakeLists.txt b/models/model_nested_events/swig/CMakeLists.txt index eb9fe681c9..523571c52e 100644 --- a/models/model_nested_events/swig/CMakeLists.txt +++ b/models/model_nested_events/swig/CMakeLists.txt @@ -1,4 +1,15 @@ cmake_minimum_required(VERSION 3.15) +cmake_policy(VERSION 3.15...3.27) + +# cmake >=3.27 +if(POLICY CMP0144) + cmake_policy(SET CMP0144 NEW) +endif(POLICY CMP0144) +# cmake >= 3.30 +if(POLICY CMP0167) + cmake_policy(SET CMP0167 NEW) +endif(POLICY CMP0167) + if(DEFINED ENV{SWIG}) set(SWIG_EXECUTABLE $ENV{SWIG}) diff --git a/models/model_neuron/CMakeLists.txt b/models/model_neuron/CMakeLists.txt index 4b580c036a..8de63e6084 100644 --- a/models/model_neuron/CMakeLists.txt +++ b/models/model_neuron/CMakeLists.txt @@ -1,11 +1,6 @@ # Build AMICI model -cmake_minimum_required(VERSION 3.15) -cmake_policy(VERSION 3.15...3.27) - -# cmake >=3.27 -if(POLICY CMP0144) - cmake_policy(SET CMP0144 NEW) -endif(POLICY CMP0144) +cmake_minimum_required(VERSION 3.22) +cmake_policy(VERSION 3.22...3.30) project(model_neuron) @@ -38,6 +33,10 @@ endif() find_package(Amici REQUIRED HINTS ${CMAKE_CURRENT_LIST_DIR}/../../build) message(STATUS "Found AMICI ${Amici_DIR}") +set_target_properties(Upstream::amici PROPERTIES + MAP_IMPORTED_CONFIG_RELWITHDEBINFO RelWithDebInfo;Release; + MAP_IMPORTED_CONFIG_RELEASE Release + MAP_IMPORTED_CONFIG_DEBUG Debug;RelWithDebInfo;) # Debug build? if("$ENV{ENABLE_AMICI_DEBUGGING}" OR "$ENV{ENABLE_GCOV_COVERAGE}") @@ -47,7 +46,6 @@ if("$ENV{ENABLE_AMICI_DEBUGGING}" OR "$ENV{ENABLE_GCOV_COVERAGE}") else() add_compile_options(-O0 -g) endif() - set(CMAKE_BUILD_TYPE "Debug") endif() # coverage options diff --git a/models/model_neuron/model_neuron.h b/models/model_neuron/model_neuron.h index 387a3cadbd..99fb505d60 100644 --- a/models/model_neuron/model_neuron.h +++ b/models/model_neuron/model_neuron.h @@ -1,6 +1,6 @@ #ifndef _amici_model_neuron_h #define _amici_model_neuron_h -/* Generated by amiwrap (R2017b) 8b324bca5b796a93094195d22a023e5f8e945ef1 */ +/* Generated by amiwrap (R2017b) d69c026f8f8a89a13f1fa1307548f2b8f7045fe1 */ #include #include #include "amici/defines.h" @@ -89,7 +89,7 @@ class Model_model_neuron : public amici::Model_ODE { amici::Model* clone() const override { return new Model_model_neuron(*this); }; - std::string getAmiciCommit() const override { return "8b324bca5b796a93094195d22a023e5f8e945ef1"; }; + std::string getAmiciCommit() const override { return "d69c026f8f8a89a13f1fa1307548f2b8f7045fe1"; }; void fJSparse(SUNMatrixContent_Sparse JSparse, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *w, const realtype *dwdx) override { JSparse_model_neuron(JSparse, t, x, p, k, h, w, dwdx); diff --git a/models/model_neuron/swig/CMakeLists.txt b/models/model_neuron/swig/CMakeLists.txt index eb9fe681c9..523571c52e 100644 --- a/models/model_neuron/swig/CMakeLists.txt +++ b/models/model_neuron/swig/CMakeLists.txt @@ -1,4 +1,15 @@ cmake_minimum_required(VERSION 3.15) +cmake_policy(VERSION 3.15...3.27) + +# cmake >=3.27 +if(POLICY CMP0144) + cmake_policy(SET CMP0144 NEW) +endif(POLICY CMP0144) +# cmake >= 3.30 +if(POLICY CMP0167) + cmake_policy(SET CMP0167 NEW) +endif(POLICY CMP0167) + if(DEFINED ENV{SWIG}) set(SWIG_EXECUTABLE $ENV{SWIG}) diff --git a/models/model_neuron_o2/CMakeLists.txt b/models/model_neuron_o2/CMakeLists.txt index 161fd4e9ce..1caad81d23 100644 --- a/models/model_neuron_o2/CMakeLists.txt +++ b/models/model_neuron_o2/CMakeLists.txt @@ -1,11 +1,6 @@ # Build AMICI model -cmake_minimum_required(VERSION 3.15) -cmake_policy(VERSION 3.15...3.27) - -# cmake >=3.27 -if(POLICY CMP0144) - cmake_policy(SET CMP0144 NEW) -endif(POLICY CMP0144) +cmake_minimum_required(VERSION 3.22) +cmake_policy(VERSION 3.22...3.30) project(model_neuron_o2) @@ -38,6 +33,10 @@ endif() find_package(Amici REQUIRED HINTS ${CMAKE_CURRENT_LIST_DIR}/../../build) message(STATUS "Found AMICI ${Amici_DIR}") +set_target_properties(Upstream::amici PROPERTIES + MAP_IMPORTED_CONFIG_RELWITHDEBINFO RelWithDebInfo;Release; + MAP_IMPORTED_CONFIG_RELEASE Release + MAP_IMPORTED_CONFIG_DEBUG Debug;RelWithDebInfo;) # Debug build? if("$ENV{ENABLE_AMICI_DEBUGGING}" OR "$ENV{ENABLE_GCOV_COVERAGE}") @@ -47,7 +46,6 @@ if("$ENV{ENABLE_AMICI_DEBUGGING}" OR "$ENV{ENABLE_GCOV_COVERAGE}") else() add_compile_options(-O0 -g) endif() - set(CMAKE_BUILD_TYPE "Debug") endif() # coverage options diff --git a/models/model_neuron_o2/model_neuron_o2.h b/models/model_neuron_o2/model_neuron_o2.h index 4ae613ee84..f29e1c88f3 100644 --- a/models/model_neuron_o2/model_neuron_o2.h +++ b/models/model_neuron_o2/model_neuron_o2.h @@ -1,6 +1,6 @@ #ifndef _amici_model_neuron_o2_h #define _amici_model_neuron_o2_h -/* Generated by amiwrap (R2017b) 8b324bca5b796a93094195d22a023e5f8e945ef1 */ +/* Generated by amiwrap (R2017b) d69c026f8f8a89a13f1fa1307548f2b8f7045fe1 */ #include #include #include "amici/defines.h" @@ -91,7 +91,7 @@ class Model_model_neuron_o2 : public amici::Model_ODE { amici::Model* clone() const override { return new Model_model_neuron_o2(*this); }; - std::string getAmiciCommit() const override { return "8b324bca5b796a93094195d22a023e5f8e945ef1"; }; + std::string getAmiciCommit() const override { return "d69c026f8f8a89a13f1fa1307548f2b8f7045fe1"; }; void fJSparse(SUNMatrixContent_Sparse JSparse, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *w, const realtype *dwdx) override { JSparse_model_neuron_o2(JSparse, t, x, p, k, h, w, dwdx); diff --git a/models/model_neuron_o2/swig/CMakeLists.txt b/models/model_neuron_o2/swig/CMakeLists.txt index eb9fe681c9..523571c52e 100644 --- a/models/model_neuron_o2/swig/CMakeLists.txt +++ b/models/model_neuron_o2/swig/CMakeLists.txt @@ -1,4 +1,15 @@ cmake_minimum_required(VERSION 3.15) +cmake_policy(VERSION 3.15...3.27) + +# cmake >=3.27 +if(POLICY CMP0144) + cmake_policy(SET CMP0144 NEW) +endif(POLICY CMP0144) +# cmake >= 3.30 +if(POLICY CMP0167) + cmake_policy(SET CMP0167 NEW) +endif(POLICY CMP0167) + if(DEFINED ENV{SWIG}) set(SWIG_EXECUTABLE $ENV{SWIG}) diff --git a/models/model_robertson/CMakeLists.txt b/models/model_robertson/CMakeLists.txt index 1a4c57353a..b9c6c2ecb1 100644 --- a/models/model_robertson/CMakeLists.txt +++ b/models/model_robertson/CMakeLists.txt @@ -1,11 +1,6 @@ # Build AMICI model -cmake_minimum_required(VERSION 3.15) -cmake_policy(VERSION 3.15...3.27) - -# cmake >=3.27 -if(POLICY CMP0144) - cmake_policy(SET CMP0144 NEW) -endif(POLICY CMP0144) +cmake_minimum_required(VERSION 3.22) +cmake_policy(VERSION 3.22...3.30) project(model_robertson) @@ -38,6 +33,10 @@ endif() find_package(Amici REQUIRED HINTS ${CMAKE_CURRENT_LIST_DIR}/../../build) message(STATUS "Found AMICI ${Amici_DIR}") +set_target_properties(Upstream::amici PROPERTIES + MAP_IMPORTED_CONFIG_RELWITHDEBINFO RelWithDebInfo;Release; + MAP_IMPORTED_CONFIG_RELEASE Release + MAP_IMPORTED_CONFIG_DEBUG Debug;RelWithDebInfo;) # Debug build? if("$ENV{ENABLE_AMICI_DEBUGGING}" OR "$ENV{ENABLE_GCOV_COVERAGE}") @@ -47,7 +46,6 @@ if("$ENV{ENABLE_AMICI_DEBUGGING}" OR "$ENV{ENABLE_GCOV_COVERAGE}") else() add_compile_options(-O0 -g) endif() - set(CMAKE_BUILD_TYPE "Debug") endif() # coverage options diff --git a/models/model_robertson/model_robertson.h b/models/model_robertson/model_robertson.h index 4f25f4d9bd..ef524b9239 100644 --- a/models/model_robertson/model_robertson.h +++ b/models/model_robertson/model_robertson.h @@ -1,6 +1,6 @@ #ifndef _amici_model_robertson_h #define _amici_model_robertson_h -/* Generated by amiwrap (R2017b) 8b324bca5b796a93094195d22a023e5f8e945ef1 */ +/* Generated by amiwrap (R2017b) d69c026f8f8a89a13f1fa1307548f2b8f7045fe1 */ #include #include #include "amici/defines.h" @@ -73,7 +73,7 @@ class Model_model_robertson : public amici::Model_DAE { amici::Model* clone() const override { return new Model_model_robertson(*this); }; - std::string getAmiciCommit() const override { return "8b324bca5b796a93094195d22a023e5f8e945ef1"; }; + std::string getAmiciCommit() const override { return "d69c026f8f8a89a13f1fa1307548f2b8f7045fe1"; }; void fJSparse(SUNMatrixContent_Sparse JSparse, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype cj, const realtype *dx, const realtype *w, const realtype *dwdx) override { JSparse_model_robertson(JSparse, t, x, p, k, h, cj, dx, w, dwdx); diff --git a/models/model_robertson/swig/CMakeLists.txt b/models/model_robertson/swig/CMakeLists.txt index eb9fe681c9..523571c52e 100644 --- a/models/model_robertson/swig/CMakeLists.txt +++ b/models/model_robertson/swig/CMakeLists.txt @@ -1,4 +1,15 @@ cmake_minimum_required(VERSION 3.15) +cmake_policy(VERSION 3.15...3.27) + +# cmake >=3.27 +if(POLICY CMP0144) + cmake_policy(SET CMP0144 NEW) +endif(POLICY CMP0144) +# cmake >= 3.30 +if(POLICY CMP0167) + cmake_policy(SET CMP0167 NEW) +endif(POLICY CMP0167) + if(DEFINED ENV{SWIG}) set(SWIG_EXECUTABLE $ENV{SWIG}) diff --git a/models/model_steadystate/CMakeLists.txt b/models/model_steadystate/CMakeLists.txt index 3d0dacaf83..05d04f8081 100644 --- a/models/model_steadystate/CMakeLists.txt +++ b/models/model_steadystate/CMakeLists.txt @@ -1,11 +1,6 @@ # Build AMICI model -cmake_minimum_required(VERSION 3.15) -cmake_policy(VERSION 3.15...3.27) - -# cmake >=3.27 -if(POLICY CMP0144) - cmake_policy(SET CMP0144 NEW) -endif(POLICY CMP0144) +cmake_minimum_required(VERSION 3.22) +cmake_policy(VERSION 3.22...3.30) project(model_steadystate) @@ -38,6 +33,10 @@ endif() find_package(Amici REQUIRED HINTS ${CMAKE_CURRENT_LIST_DIR}/../../build) message(STATUS "Found AMICI ${Amici_DIR}") +set_target_properties(Upstream::amici PROPERTIES + MAP_IMPORTED_CONFIG_RELWITHDEBINFO RelWithDebInfo;Release; + MAP_IMPORTED_CONFIG_RELEASE Release + MAP_IMPORTED_CONFIG_DEBUG Debug;RelWithDebInfo;) # Debug build? if("$ENV{ENABLE_AMICI_DEBUGGING}" OR "$ENV{ENABLE_GCOV_COVERAGE}") @@ -47,7 +46,6 @@ if("$ENV{ENABLE_AMICI_DEBUGGING}" OR "$ENV{ENABLE_GCOV_COVERAGE}") else() add_compile_options(-O0 -g) endif() - set(CMAKE_BUILD_TYPE "Debug") endif() # coverage options diff --git a/models/model_steadystate/model_steadystate.h b/models/model_steadystate/model_steadystate.h index ca02261ef8..a27a5a74df 100644 --- a/models/model_steadystate/model_steadystate.h +++ b/models/model_steadystate/model_steadystate.h @@ -1,6 +1,6 @@ #ifndef _amici_model_steadystate_h #define _amici_model_steadystate_h -/* Generated by amiwrap (R2017b) 8b324bca5b796a93094195d22a023e5f8e945ef1 */ +/* Generated by amiwrap (R2017b) d69c026f8f8a89a13f1fa1307548f2b8f7045fe1 */ #include #include #include "amici/defines.h" @@ -72,7 +72,7 @@ class Model_model_steadystate : public amici::Model_ODE { amici::Model* clone() const override { return new Model_model_steadystate(*this); }; - std::string getAmiciCommit() const override { return "8b324bca5b796a93094195d22a023e5f8e945ef1"; }; + std::string getAmiciCommit() const override { return "d69c026f8f8a89a13f1fa1307548f2b8f7045fe1"; }; void fJSparse(SUNMatrixContent_Sparse JSparse, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *w, const realtype *dwdx) override { JSparse_model_steadystate(JSparse, t, x, p, k, h, w, dwdx); diff --git a/models/model_steadystate/swig/CMakeLists.txt b/models/model_steadystate/swig/CMakeLists.txt index eb9fe681c9..523571c52e 100644 --- a/models/model_steadystate/swig/CMakeLists.txt +++ b/models/model_steadystate/swig/CMakeLists.txt @@ -1,4 +1,15 @@ cmake_minimum_required(VERSION 3.15) +cmake_policy(VERSION 3.15...3.27) + +# cmake >=3.27 +if(POLICY CMP0144) + cmake_policy(SET CMP0144 NEW) +endif(POLICY CMP0144) +# cmake >= 3.30 +if(POLICY CMP0167) + cmake_policy(SET CMP0167 NEW) +endif(POLICY CMP0167) + if(DEFINED ENV{SWIG}) set(SWIG_EXECUTABLE $ENV{SWIG}) diff --git a/python/examples/example_steadystate/ExampleSteadystate.ipynb b/python/examples/example_steadystate/ExampleSteadystate.ipynb index d9f6ae635d..147ae473b6 100644 --- a/python/examples/example_steadystate/ExampleSteadystate.ipynb +++ b/python/examples/example_steadystate/ExampleSteadystate.ipynb @@ -392,10 +392,10 @@ "source": [ "model = model_module.getModel()\n", "\n", - "print(\"Model name: \", model.getName())\n", - "print(\"Model parameters:\", model.getParameterIds())\n", - "print(\"Model outputs: \", model.getObservableIds())\n", - "print(\"Model states: \", model.getStateIds())" + "print(\"Model name: \", model.getName())\n", + "print(\"Model parameters: \", model.getParameterIds())\n", + "print(\"Model outputs: \", model.getObservableIds())\n", + "print(\"Model state variables: \", model.getStateIds())" ] }, { @@ -985,6 +985,21 @@ "print(\"Log-likelihood %f\" % rdata[\"llh\"])" ] }, + { + "cell_type": "markdown", + "metadata": {}, + "source": "The provided measurements can be visualized together with the simulation results by passing the `Expdata` to `amici.plotting.plot_observable_trajectories`:" + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "amici.plotting.plot_observable_trajectories(rdata, edata=edata)\n", + "plt.legend(loc=\"center left\", bbox_to_anchor=(1.04, 0.5))" + ] + }, { "cell_type": "markdown", "metadata": {}, diff --git a/python/sdist/MANIFEST.in b/python/sdist/MANIFEST.in index 5be0b1ebc6..7cd9792a20 100644 --- a/python/sdist/MANIFEST.in +++ b/python/sdist/MANIFEST.in @@ -11,6 +11,7 @@ include LICENSE.md exclude amici/*.so exclude amici/*.dll prune **/build +prune **/install prune amici/share prune amici/lib prune amici/ThirdParty/SuiteSparse/lib diff --git a/python/sdist/amici/petab/cli/import_petab.py b/python/sdist/amici/petab/cli/import_petab.py index b124b4d98b..d61cda5362 100644 --- a/python/sdist/amici/petab/cli/import_petab.py +++ b/python/sdist/amici/petab/cli/import_petab.py @@ -1,6 +1,6 @@ import argparse -import petab +import petab.v1 as petab from ..petab_import import import_model_sbml from petab.v1.models.sbml_model import SbmlModel diff --git a/python/sdist/pyproject.toml b/python/sdist/pyproject.toml index 621fd43952..fd4e7e22b7 100644 --- a/python/sdist/pyproject.toml +++ b/python/sdist/pyproject.toml @@ -61,7 +61,6 @@ classifiers = [ petab = ["petab>=0.4.0"] pysb = ["pysb>=1.13.1"] test = [ - "benchmark_models_petab @ git+https://github.com/Benchmarking-Initiative/Benchmark-Models-PEtab.git@master#subdirectory=src/python", "h5py", "pytest", "pytest-cov", diff --git a/python/sdist/setup.py b/python/sdist/setup.py index d01a0ec75c..dc8372cdc3 100755 --- a/python/sdist/setup.py +++ b/python/sdist/setup.py @@ -45,73 +45,37 @@ def get_extensions(): f"-DCMAKE_MODULE_PATH={prefix_path.as_posix()}", ] - # SuiteSparse Config - suitesparse_config = CMakeExtension( - name="SuiteSparse_config", + debug_build = os.getenv("ENABLE_AMICI_DEBUGGING", "").lower() in [ + "1", + "true", + ] or os.getenv("ENABLE_GCOV_COVERAGE", "").lower() in ["1", "true"] + build_type = "Debug" if debug_build else "Release" + + # SuiteSparse + suitesparse = CMakeExtension( + name="SuiteSparse", install_prefix="amici", - source_dir="amici/ThirdParty/SuiteSparse/SuiteSparse_config", + source_dir="amici/ThirdParty/SuiteSparse", + cmake_build_type=build_type, cmake_configure_options=[ *global_cmake_configure_options, "-DCMAKE_POSITION_INDEPENDENT_CODE=ON", "-DBUILD_SHARED_LIBS=OFF", + "-DBUILD_TESTING=OFF", # Building SuiteSparse_config does not require a BLAS # we just set BLAS_LIBRARIES to skip the search, # the value is not used # "-DBLA_VENDOR=All", "-DBLAS_LIBRARIES=dummy", "-DSUITESPARSE_USE_64BIT_BLAS=ON", + "-DSUITESPARSE_ENABLE_PROJECTS=amd;btf;colamd;klu", "-DSUITESPARSE_USE_CUDA=OFF", "-DSUITESPARSE_USE_FORTRAN=OFF", - ], - ) - # SuiteSparse AMD - amd = CMakeExtension( - name="amd", - install_prefix="amici", - source_dir="amici/ThirdParty/SuiteSparse/AMD", - cmake_configure_options=[ - *global_cmake_configure_options, - "-DBUILD_SHARED_LIBS=OFF", - "-DCMAKE_POSITION_INDEPENDENT_CODE=ON", - "-DSUITESPARSE_USE_FORTRAN=OFF", - ], - ) - # SuiteSparse BTF - btf = CMakeExtension( - name="btf", - install_prefix="amici", - source_dir="amici/ThirdParty/SuiteSparse/BTF", - cmake_configure_options=[ - *global_cmake_configure_options, - "-DSUITESPARSE_USE_FORTRAN=OFF", - "-DBUILD_SHARED_LIBS=OFF", - "-DCMAKE_POSITION_INDEPENDENT_CODE=ON", - ], - ) - # SuiteSparse COLAMD - colamd = CMakeExtension( - name="colamd", - install_prefix="amici", - source_dir="amici/ThirdParty/SuiteSparse/COLAMD", - cmake_configure_options=[ - *global_cmake_configure_options, - "-DSUITESPARSE_USE_FORTRAN=OFF", - "-DBUILD_SHARED_LIBS=OFF", - "-DCMAKE_POSITION_INDEPENDENT_CODE=ON", - ], - ) - # SuiteSparse KLU - klu = CMakeExtension( - name="klu", - install_prefix="amici", - source_dir="amici/ThirdParty/SuiteSparse/KLU", - cmake_configure_options=[ - *global_cmake_configure_options, + "-DSUITESPARSE_USE_PYTHON=OFF", + "-DSUITESPARSE_USE_OPENMP=OFF", + "-DSUITESPARSE_CONFIG_USE_OPENMP=OFF", + "-DCHOLMOD_CAMD=OFF", "-DKLU_USE_CHOLMOD=OFF", - "-DSUITESPARSE_USE_CUDA=OFF", - "-DSUITESPARSE_USE_FORTRAN=OFF", - "-DBUILD_SHARED_LIBS=OFF", - "-DCMAKE_POSITION_INDEPENDENT_CODE=ON", ], ) # SUNDIALS @@ -119,6 +83,7 @@ def get_extensions(): name="sundials", install_prefix="amici", source_dir="amici/ThirdParty/sundials", + cmake_build_type=build_type, cmake_configure_options=[ *global_cmake_configure_options, "-DBUILD_ARKODE=OFF", @@ -141,14 +106,11 @@ def get_extensions(): ], ) # AMICI - debug_build = os.getenv("ENABLE_AMICI_DEBUGGING", "").lower() in [ - "1", - "true", - ] or os.getenv("ENABLE_GCOV_COVERAGE", "").lower() in ["1", "true"] amici_ext = CMakeExtension( name="amici", install_prefix="amici", source_dir="amici", + cmake_build_type=build_type, cmake_configure_options=[ *global_cmake_configure_options, "-Werror=dev" @@ -159,10 +121,9 @@ def get_extensions(): "-DAMICI_PYTHON_BUILD_EXT_ONLY=ON", f"-DPython3_EXECUTABLE={Path(sys.executable).as_posix()}", ], - cmake_build_type="Debug" if debug_build else "Release", ) # Order matters! - return [suitesparse_config, amd, btf, colamd, klu, sundials, amici_ext] + return [suitesparse, sundials, amici_ext] def main(): diff --git a/scripts/buildSuiteSparse.sh b/scripts/buildSuiteSparse.sh index 93341810f8..90ec1d8e4a 100755 --- a/scripts/buildSuiteSparse.sh +++ b/scripts/buildSuiteSparse.sh @@ -8,14 +8,33 @@ script_path=$(dirname "$BASH_SOURCE") amici_path=$(cd "$script_path/.." && pwd) suitesparse_root="${amici_path}/ThirdParty/SuiteSparse" -for subdir in SuiteSparse_config BTF AMD COLAMD KLU; do - export CMAKE_OPTIONS="-DSUITESPARSE_USE_CUDA=OFF -DSUITESPARSE_USE_FORTRAN=OFF" - if [ $subdir = "SuiteSparse_config" ]; then - export CMAKE_OPTIONS="$CMAKE_OPTIONS -DCMAKE_POSITION_INDEPENDENT_CODE=ON -DBLA_VENDOR=All -DSUITESPARSE_USE_64BIT_BLAS=ON -DBLAS_LIBRARIES=dummy" - elif [ $subdir = "KLU" ]; then - export CMAKE_OPTIONS="$CMAKE_OPTIONS -DKLU_USE_CHOLMOD=OFF" - fi +cd "${suitesparse_root}/build" - cd "${suitesparse_root}/${subdir}" && make local install -done +if [[ "${GITHUB_REPOSITORY:-}" = *"/AMICI" ]] || + [ "${ENABLE_AMICI_DEBUGGING:-}" = TRUE ]; then + # Running on CI server + build_type="Debug" +else + build_type="RelWithDebInfo" +fi + +cmake -DSUITESPARSE_ENABLE_PROJECTS="amd;btf;colamd;klu" \ + -DCMAKE_BUILD_TYPE="$build_type" \ + -DCMAKE_INSTALL_PREFIX="${suitesparse_root}/install" \ + -DCHOLMOD_CAMD=OFF \ + -DKLU_USE_CHOLMOD=OFF \ + -DSUITESPARSE_CONFIG_USE_OPENMP=OFF \ + -DSUITESPARSE_USE_CUDA=OFF \ + -DSUITESPARSE_USE_FORTRAN=OFF \ + -DSUITESPARSE_USE_PYTHON=OFF \ + -DSUITESPARSE_USE_OPENMP=OFF \ + -DCMAKE_POSITION_INDEPENDENT_CODE=ON \ + -DSUITESPARSE_USE_64BIT_BLAS=ON \ + -DBLAS_LIBRARIES=dummy \ + -DBUILD_SHARED_LIBS=OFF \ + -DBUILD_TESTING=OFF \ + .. + +cmake --build . +cmake --install . diff --git a/scripts/buildSundials.sh b/scripts/buildSundials.sh index f3410ef31a..642c607cb8 100755 --- a/scripts/buildSundials.sh +++ b/scripts/buildSundials.sh @@ -14,7 +14,8 @@ sundials_build_path="${amici_path}/ThirdParty/sundials/build/" cmake=${CMAKE:-cmake} make=${MAKE:-make} -if [[ $GITHUB_ACTIONS = true ]]; then +if [[ "${GITHUB_REPOSITORY:-}" = *"/AMICI" ]] || + [ "${ENABLE_AMICI_DEBUGGING:-}" = TRUE ]; then # Running on CI server build_type="Debug" else @@ -47,8 +48,8 @@ ${cmake} -DCMAKE_INSTALL_PREFIX="${sundials_build_path}" \ -DEXAMPLES_ENABLE_C=OFF \ -DEXAMPLES_INSTALL=OFF \ -DENABLE_KLU=ON \ - -DKLU_LIBRARY_DIR="${suitesparse_root}/lib" \ - -DKLU_INCLUDE_DIR="${suitesparse_root}/include/suitesparse" \ + -DKLU_LIBRARY_DIR="${suitesparse_root}/install/lib" \ + -DKLU_INCLUDE_DIR="${suitesparse_root}/install/include/suitesparse" \ ${SuperLUMT} \ .. diff --git a/scripts/installAmiciSource.sh b/scripts/installAmiciSource.sh index 1c0a975708..6769ebc100 100755 --- a/scripts/installAmiciSource.sh +++ b/scripts/installAmiciSource.sh @@ -33,7 +33,7 @@ fi export PYTHON_EXECUTABLE="${AMICI_PATH}/venv/bin/python" python -m pip install --upgrade pip wheel -python -m pip install --upgrade pip setuptools cmake_build_extension==0.6.0 numpy +python -m pip install --upgrade pip setuptools cmake_build_extension==0.6.0 numpy petab python -m pip install git+https://github.com/FFroehlich/pysb@fix_pattern_matching # pin to PR for SPM with compartments AMICI_BUILD_TEMP="${AMICI_PATH}/python/sdist/build/temp" \ python -m pip install --verbose -e "${AMICI_PATH}/python/sdist[petab,test,vis]" --no-build-isolation diff --git a/src/CMakeLists.template.cmake b/src/CMakeLists.template.cmake index acffe7e247..4ead82b79d 100644 --- a/src/CMakeLists.template.cmake +++ b/src/CMakeLists.template.cmake @@ -1,11 +1,6 @@ # Build AMICI model -cmake_minimum_required(VERSION 3.15) -cmake_policy(VERSION 3.15...3.30) - -# cmake >=3.27 -if(POLICY CMP0144) - cmake_policy(SET CMP0144 NEW) -endif(POLICY CMP0144) +cmake_minimum_required(VERSION 3.22) +cmake_policy(VERSION 3.22...3.30) project(TPL_MODELNAME) diff --git a/src/model.cpp b/src/model.cpp index c27a5310d9..4658f41008 100644 --- a/src/model.cpp +++ b/src/model.cpp @@ -177,19 +177,15 @@ Model::Model( ModelDimensions const& model_dimensions, SimulationParameters simulation_parameters, SecondOrderMode o2mode, std::vector idlist, std::vector z2event, - bool const pythonGenerated, int const ndxdotdp_explicit, - int const ndxdotdx_explicit, int const w_recursion_depth, std::map> state_independent_events ) : ModelDimensions(model_dimensions) - , pythonGenerated(pythonGenerated) , o2mode(o2mode) , idlist(std::move(idlist)) , state_independent_events_(std::move(state_independent_events)) , derived_state_(model_dimensions) , z2event_(std::move(z2event)) , state_is_non_negative_(nx_solver, false) - , w_recursion_depth_(w_recursion_depth) , simulation_parameters_(std::move(simulation_parameters)) { Expects( model_dimensions.np @@ -217,60 +213,6 @@ Model::Model( root_initial_values_.resize(ne, true); - /* If Matlab wrapped: dxdotdp is a full AmiVector, - if Python wrapped: dxdotdp_explicit and dxdotdp_implicit are CSC matrices - */ - if (pythonGenerated) { - - dwdw_ = SUNMatrixWrapper(nw, nw, ndwdw, CSC_MAT); - // size dynamically adapted for dwdx_ and dwdp_ - derived_state_.dwdx_ = SUNMatrixWrapper(nw, nx_solver, 0, CSC_MAT); - derived_state_.dwdp_ = SUNMatrixWrapper(nw, np(), 0, CSC_MAT); - - for (int irec = 0; irec <= w_recursion_depth_; ++irec) { - /* for the first element we know the exact size, while for all - others we guess the size*/ - dwdp_hierarchical_.emplace_back( - SUNMatrixWrapper(nw, np(), irec * ndwdw + ndwdp, CSC_MAT) - ); - dwdx_hierarchical_.emplace_back( - SUNMatrixWrapper(nw, nx_solver, irec * ndwdw + ndwdx, CSC_MAT) - ); - } - assert( - gsl::narrow(dwdp_hierarchical_.size()) - == w_recursion_depth_ + 1 - ); - assert( - gsl::narrow(dwdx_hierarchical_.size()) - == w_recursion_depth_ + 1 - ); - - derived_state_.dxdotdp_explicit - = SUNMatrixWrapper(nx_solver, np(), ndxdotdp_explicit, CSC_MAT); - // guess size, will be dynamically reallocated - derived_state_.dxdotdp_implicit - = SUNMatrixWrapper(nx_solver, np(), ndwdp + ndxdotdw, CSC_MAT); - derived_state_.dxdotdx_explicit = SUNMatrixWrapper( - nx_solver, nx_solver, ndxdotdx_explicit, CSC_MAT - ); - // guess size, will be dynamically reallocated - derived_state_.dxdotdx_implicit - = SUNMatrixWrapper(nx_solver, nx_solver, ndwdx + ndxdotdw, CSC_MAT); - // dynamically allocate on first call - derived_state_.dxdotdp_full - = SUNMatrixWrapper(nx_solver, np(), 0, CSC_MAT); - - for (int iytrue = 0; iytrue < nytrue; ++iytrue) - derived_state_.dJydy_.emplace_back( - SUNMatrixWrapper(nJ, ny, ndJydy.at(iytrue), CSC_MAT) - ); - } else { - derived_state_.dwdx_ = SUNMatrixWrapper(nw, nx_solver, ndwdx, CSC_MAT); - derived_state_.dwdp_ = SUNMatrixWrapper(nw, np(), ndwdp, CSC_MAT); - derived_state_.dJydy_matlab_ - = std::vector(nJ * nytrue * ny, 0.0); - } requireSensitivitiesForAllParameters(); } @@ -288,8 +230,7 @@ bool operator==(Model const& a, Model const& b) { && (a.state_is_non_negative_ == b.state_is_non_negative_) && (a.sigma_res_ == b.sigma_res_) && (a.min_sigma_ == b.min_sigma_) && (a.state_ == b.state_) - && (a.steadystate_mask_.getVector() - == b.steadystate_mask_.getVector()); + && (a.steadystate_mask_ == b.steadystate_mask_); } bool operator==(ModelDimensions const& a, ModelDimensions const& b) { @@ -381,7 +322,7 @@ void Model::initializeSplines() { splines_ = fcreate_splines( state_.unscaledParameters.data(), state_.fixedParameters.data() ); - state_.spl_.resize(splines_.size(), 0.0); + derived_state_.spl_.resize(splines_.size(), 0.0); for (auto& spline : splines_) { spline.compute_coefficients(); } @@ -1985,20 +1926,20 @@ void Model::fsx0_fixedParameters(AmiVectorArray& sx, AmiVector const& x) { void Model::fsdx0() {} -void Model::fx_rdata(AmiVector& x_rdata, AmiVector const& x) { +void Model::fx_rdata(gsl::span x_rdata, AmiVector const& x) { fx_rdata( x_rdata.data(), computeX_pos(x), state_.total_cl.data(), state_.unscaledParameters.data(), state_.fixedParameters.data() ); if (always_check_finite_) checkFinite( - x_rdata.getVector(), ModelQuantity::x_rdata, + x_rdata, ModelQuantity::x_rdata, std::numeric_limits::quiet_NaN() ); } void Model::fsx_rdata( - AmiVectorArray& sx_rdata, AmiVectorArray const& sx, + gsl::span sx_rdata, AmiVectorArray const& sx, AmiVector const& x_solver ) { realtype* stcl = nullptr; @@ -2006,7 +1947,7 @@ void Model::fsx_rdata( if (ncl() > 0) stcl = &state_.stotal_cl.at(plist(ip) * ncl()); fsx_rdata( - sx_rdata.data(ip), sx.data(ip), stcl, + &sx_rdata[ip * nx_rdata], sx.data(ip), stcl, state_.unscaledParameters.data(), state_.fixedParameters.data(), x_solver.data(), state_.total_cl.data(), plist(ip) ); @@ -2107,7 +2048,7 @@ void Model::fdydp(realtype const t, AmiVector const& x) { state_.unscaledParameters.data(), state_.fixedParameters.data(), state_.h.data(), plist(ip), derived_state_.w_.data(), state_.total_cl.data(), state_.stotal_cl.data(), - state_.spl_.data(), derived_state_.sspl_.data() + derived_state_.spl_.data(), derived_state_.sspl_.data() ); } else { fdydp( @@ -2251,7 +2192,6 @@ void Model::fdJydy(int const it, AmiVector const& x, ExpData const& edata) { if (pythonGenerated) { fdJydsigma(it, x, edata); fdsigmaydy(it, &edata); - SUNMatrixWrapper tmp_dense(nJ, ny); setNaNtoZero(derived_state_.dJydsigma_); setNaNtoZero(derived_state_.dsigmaydy_); @@ -2276,15 +2216,17 @@ void Model::fdJydy(int const it, AmiVector const& x, ExpData const& edata) { // dJydy += dJydsigma * dsigmaydy // C(nJ,ny) A(nJ,ny) * B(ny,ny) // sparse dense dense - tmp_dense.zero(); + derived_state_.dJydy_dense_.zero(); amici_dgemm( BLASLayout::colMajor, BLASTranspose::noTrans, BLASTranspose::noTrans, nJ, ny, ny, 1.0, &derived_state_.dJydsigma_.at(iyt * nJ * ny), nJ, - derived_state_.dsigmaydy_.data(), ny, 1.0, tmp_dense.data(), nJ + derived_state_.dsigmaydy_.data(), ny, 1.0, + derived_state_.dJydy_dense_.data(), nJ ); - auto tmp_sparse = SUNMatrixWrapper(tmp_dense, 0.0, CSC_MAT); + auto tmp_sparse + = SUNMatrixWrapper(derived_state_.dJydy_dense_, 0.0, CSC_MAT); auto ret = SUNMatScaleAdd( 1.0, derived_state_.dJydy_.at(iyt), tmp_sparse ); @@ -2864,7 +2806,7 @@ void Model::fdJrzdsigma( void Model::fspl(realtype const t) { for (int ispl = 0; ispl < nspl; ispl++) - state_.spl_[ispl] = splines_[ispl].get_value(t); + derived_state_.spl_[ispl] = splines_[ispl].get_value(t); } void Model::fsspl(realtype const t) { @@ -2872,8 +2814,9 @@ void Model::fsspl(realtype const t) { realtype* sspl_data = derived_state_.sspl_.data(); for (int ip = 0; ip < nplist(); ip++) { for (int ispl = 0; ispl < nspl; ispl++) - sspl_data[ispl + nspl * plist(ip)] - = splines_[ispl].get_sensitivity(t, ip, state_.spl_[ispl]); + sspl_data[ispl + nspl * plist(ip)] = splines_[ispl].get_sensitivity( + t, ip, derived_state_.spl_[ispl] + ); } } @@ -2884,7 +2827,7 @@ void Model::fw(realtype const t, realtype const* x, bool include_static) { fspl(t); fw(derived_state_.w_.data(), t, x, state_.unscaledParameters.data(), state_.fixedParameters.data(), state_.h.data(), state_.total_cl.data(), - state_.spl_.data(), include_static); + derived_state_.spl_.data(), include_static); if (always_check_finite_) { checkFinite(derived_state_.w_, ModelQuantity::w, t); @@ -2903,26 +2846,26 @@ void Model::fdwdp(realtype const t, realtype const* x, bool include_static) { fsspl(t); fdwdw(t, x, include_static); if (include_static) { - dwdp_hierarchical_.at(0).zero(); - fdwdp_colptrs(dwdp_hierarchical_.at(0)); - fdwdp_rowvals(dwdp_hierarchical_.at(0)); + derived_state_.dwdp_hierarchical_.at(0).zero(); + fdwdp_colptrs(derived_state_.dwdp_hierarchical_.at(0)); + fdwdp_rowvals(derived_state_.dwdp_hierarchical_.at(0)); } fdwdp( - dwdp_hierarchical_.at(0).data(), t, x, + derived_state_.dwdp_hierarchical_.at(0).data(), t, x, state_.unscaledParameters.data(), state_.fixedParameters.data(), state_.h.data(), derived_state_.w_.data(), state_.total_cl.data(), - state_.stotal_cl.data(), state_.spl_.data(), + state_.stotal_cl.data(), derived_state_.spl_.data(), derived_state_.sspl_.data(), include_static ); - for (int irecursion = 1; irecursion <= w_recursion_depth_; + for (int irecursion = 1; irecursion <= w_recursion_depth; irecursion++) { - dwdw_.sparse_multiply( - dwdp_hierarchical_.at(irecursion), - dwdp_hierarchical_.at(irecursion - 1) + derived_state_.dwdw_.sparse_multiply( + derived_state_.dwdp_hierarchical_.at(irecursion), + derived_state_.dwdp_hierarchical_.at(irecursion - 1) ); } - derived_state_.dwdp_.sparse_sum(dwdp_hierarchical_); + derived_state_.dwdp_.sparse_sum(derived_state_.dwdp_hierarchical_); } else { if (!derived_state_.dwdp_.capacity()) @@ -2932,7 +2875,7 @@ void Model::fdwdp(realtype const t, realtype const* x, bool include_static) { derived_state_.dwdp_.data(), t, x, state_.unscaledParameters.data(), state_.fixedParameters.data(), state_.h.data(), derived_state_.w_.data(), state_.total_cl.data(), - state_.stotal_cl.data(), state_.spl_.data(), + state_.stotal_cl.data(), derived_state_.spl_.data(), derived_state_.sspl_.data() ); } @@ -2950,31 +2893,31 @@ void Model::fdwdx(realtype const t, realtype const* x, bool include_static) { derived_state_.dwdx_.zero(); if (pythonGenerated) { - if (!dwdx_hierarchical_.at(0).capacity()) + if (!derived_state_.dwdx_hierarchical_.at(0).capacity()) return; fdwdw(t, x, include_static); if (include_static) { - dwdx_hierarchical_.at(0).zero(); - fdwdx_colptrs(dwdx_hierarchical_.at(0)); - fdwdx_rowvals(dwdx_hierarchical_.at(0)); + derived_state_.dwdx_hierarchical_.at(0).zero(); + fdwdx_colptrs(derived_state_.dwdx_hierarchical_.at(0)); + fdwdx_rowvals(derived_state_.dwdx_hierarchical_.at(0)); } fdwdx( - dwdx_hierarchical_.at(0).data(), t, x, + derived_state_.dwdx_hierarchical_.at(0).data(), t, x, state_.unscaledParameters.data(), state_.fixedParameters.data(), state_.h.data(), derived_state_.w_.data(), state_.total_cl.data(), - state_.spl_.data(), include_static + derived_state_.spl_.data(), include_static ); - for (int irecursion = 1; irecursion <= w_recursion_depth_; + for (int irecursion = 1; irecursion <= w_recursion_depth; irecursion++) { - dwdw_.sparse_multiply( - dwdx_hierarchical_.at(irecursion), - dwdx_hierarchical_.at(irecursion - 1) + derived_state_.dwdw_.sparse_multiply( + derived_state_.dwdx_hierarchical_.at(irecursion), + derived_state_.dwdx_hierarchical_.at(irecursion - 1) ); } - derived_state_.dwdx_.sparse_sum(dwdx_hierarchical_); + derived_state_.dwdx_.sparse_sum(derived_state_.dwdx_hierarchical_); } else { if (!derived_state_.dwdx_.capacity()) @@ -2983,7 +2926,8 @@ void Model::fdwdx(realtype const t, realtype const* x, bool include_static) { fdwdx( derived_state_.dwdx_.data(), t, x, state_.unscaledParameters.data(), state_.fixedParameters.data(), state_.h.data(), - derived_state_.w_.data(), state_.total_cl.data(), state_.spl_.data() + derived_state_.w_.data(), state_.total_cl.data(), + derived_state_.spl_.data() ); } @@ -2997,19 +2941,19 @@ void Model::fdwdw(realtype const t, realtype const* x, bool include_static) { return; if (include_static) { - dwdw_.zero(); - fdwdw_colptrs(dwdw_); - fdwdw_rowvals(dwdw_); + derived_state_.dwdw_.zero(); + fdwdw_colptrs(derived_state_.dwdw_); + fdwdw_rowvals(derived_state_.dwdw_); } fdwdw( - dwdw_.data(), t, x, state_.unscaledParameters.data(), + derived_state_.dwdw_.data(), t, x, state_.unscaledParameters.data(), state_.fixedParameters.data(), state_.h.data(), derived_state_.w_.data(), state_.total_cl.data(), include_static ); if (always_check_finite_) { - checkFinite(dwdw_, ModelQuantity::dwdw, t); + checkFinite(derived_state_.dwdw_, ModelQuantity::dwdw, t); } } @@ -3137,11 +3081,9 @@ std::vector Model::get_trigger_timepoints() const { return trigger_timepoints; } -void Model::set_steadystate_mask(std::vector const& mask) { +void Model::set_steadystate_mask(std::vector const& mask) { if (mask.size() == 0) { - if (steadystate_mask_.getLength() != 0) { - steadystate_mask_ = AmiVector(); - } + steadystate_mask_.clear(); return; } @@ -3151,7 +3093,7 @@ void Model::set_steadystate_mask(std::vector const& mask) { gsl::narrow(mask.size()), nx_solver ); - steadystate_mask_ = AmiVector(mask); + steadystate_mask_ = mask; } const_N_Vector Model::computeX_pos(const_N_Vector x) { diff --git a/src/model_header.template.h b/src/model_header.template.h index 932fdeb1a0..bc87da5a03 100644 --- a/src/model_header.template.h +++ b/src/model_header.template.h @@ -135,7 +135,11 @@ class Model_TPL_MODELNAME : public amici::Model_TPL_MODEL_TYPE_UPPER { TPL_NDTOTALCLDXRDATA, // ndtotal_cldx_rdata 0, // nnz TPL_UBW, // ubw - TPL_LBW // lbw + TPL_LBW, // lbw + true, // pythonGenerated + TPL_NDXDOTDP_EXPLICIT, // ndxdotdp_explicit + TPL_NDXDOTDX_EXPLICIT, // ndxdotdx_explicit + TPL_W_RECURSION_DEPTH // w_recursion_depth ), amici::SimulationParameters( std::vector{TPL_FIXED_PARAMETERS}, // fixedParameters @@ -144,10 +148,6 @@ class Model_TPL_MODELNAME : public amici::Model_TPL_MODEL_TYPE_UPPER { TPL_O2MODE, // o2mode std::vector{TPL_ID}, // idlist std::vector{TPL_Z2EVENT}, // z2events - true, // pythonGenerated - TPL_NDXDOTDP_EXPLICIT, // ndxdotdp_explicit - TPL_NDXDOTDX_EXPLICIT, // ndxdotdx_explicit - TPL_W_RECURSION_DEPTH, // w_recursion_depth {TPL_STATE_INDEPENDENT_EVENTS} // state-independent events ) { root_initial_values_ = std::vector( diff --git a/src/model_state.cpp b/src/model_state.cpp index b6d1d9c850..ade1d4c847 100644 --- a/src/model_state.cpp +++ b/src/model_state.cpp @@ -19,6 +19,65 @@ ModelStateDerived::ModelStateDerived(ModelDimensions const& dim) , w_(dim.nw) , x_rdata_(dim.nx_rdata, 0.0) , sx_rdata_(dim.nx_rdata, 0.0) - , x_pos_tmp_(dim.nx_solver) {} + , x_pos_tmp_(dim.nx_solver) { + + // If Matlab wrapped: dxdotdp is a full AmiVector, + // if Python wrapped: dxdotdp_explicit and dxdotdp_implicit are CSC matrices + if (dim.pythonGenerated) { + + dwdw_ = SUNMatrixWrapper(dim.nw, dim.nw, dim.ndwdw, CSC_MAT); + // size dynamically adapted for dwdx_ and dwdp_ + dwdx_ = SUNMatrixWrapper(dim.nw, dim.nx_solver, 0, CSC_MAT); + dwdp_ = SUNMatrixWrapper(dim.nw, dim.np, 0, CSC_MAT); + + for (int irec = 0; irec <= dim.w_recursion_depth; ++irec) { + /* for the first element we know the exact size, while for all + others we guess the size*/ + dwdp_hierarchical_.emplace_back(SUNMatrixWrapper( + dim.nw, dim.np, irec * dim.ndwdw + dim.ndwdp, CSC_MAT + )); + dwdx_hierarchical_.emplace_back(SUNMatrixWrapper( + dim.nw, dim.nx_solver, irec * dim.ndwdw + dim.ndwdx, CSC_MAT + )); + } + assert( + gsl::narrow(dwdp_hierarchical_.size()) + == dim.w_recursion_depth + 1 + ); + assert( + gsl::narrow(dwdx_hierarchical_.size()) + == dim.w_recursion_depth + 1 + ); + + dxdotdp_explicit = SUNMatrixWrapper( + dim.nx_solver, dim.np, dim.ndxdotdp_explicit, CSC_MAT + ); + // guess size, will be dynamically reallocated + dxdotdp_implicit = SUNMatrixWrapper( + dim.nx_solver, dim.np, dim.ndwdp + dim.ndxdotdw, CSC_MAT + ); + dxdotdx_explicit = SUNMatrixWrapper( + dim.nx_solver, dim.nx_solver, dim.ndxdotdx_explicit, CSC_MAT + ); + // guess size, will be dynamically reallocated + dxdotdx_implicit = SUNMatrixWrapper( + dim.nx_solver, dim.nx_solver, dim.ndwdx + dim.ndxdotdw, CSC_MAT + ); + // dynamically allocate on first call + dxdotdp_full = SUNMatrixWrapper(dim.nx_solver, dim.np, 0, CSC_MAT); + + for (int iytrue = 0; iytrue < dim.nytrue; ++iytrue) + dJydy_.emplace_back( + SUNMatrixWrapper(dim.nJ, dim.ny, dim.ndJydy.at(iytrue), CSC_MAT) + ); + + dJydy_dense_ = SUNMatrixWrapper(dim.nJ, dim.ny); + } else { + dwdx_ = SUNMatrixWrapper(dim.nw, dim.nx_solver, dim.ndwdx, CSC_MAT); + dwdp_ = SUNMatrixWrapper(dim.nw, dim.np, dim.ndwdp, CSC_MAT); + dJydy_matlab_ + = std::vector(dim.nJ * dim.nytrue * dim.ny, 0.0); + } +} } // namespace amici diff --git a/src/newton_solver.cpp b/src/newton_solver.cpp index 72bf0a3d64..bcc189328b 100644 --- a/src/newton_solver.cpp +++ b/src/newton_solver.cpp @@ -10,56 +10,31 @@ namespace amici { -NewtonSolver::NewtonSolver(Model const& model) +NewtonSolver::NewtonSolver(Model const& model, LinearSolver linsol_type) : xdot_(model.nx_solver) , x_(model.nx_solver) , xB_(model.nJ * model.nx_solver) - , dxB_(model.nJ * model.nx_solver) {} - -std::unique_ptr -NewtonSolver::getSolver(Solver const& simulationSolver, Model const& model) { - - std::unique_ptr solver; - - switch (simulationSolver.getLinearSolver()) { - - /* DIRECT SOLVERS */ - case LinearSolver::dense: - solver.reset(new NewtonSolverDense(model)); - break; - - case LinearSolver::band: - throw NewtonFailure(AMICI_NOT_IMPLEMENTED, "getSolver"); - - case LinearSolver::LAPACKDense: - throw NewtonFailure(AMICI_NOT_IMPLEMENTED, "getSolver"); - - case LinearSolver::LAPACKBand: - throw NewtonFailure(AMICI_NOT_IMPLEMENTED, "getSolver"); - - case LinearSolver::diag: - throw NewtonFailure(AMICI_NOT_IMPLEMENTED, "getSolver"); - - /* ITERATIVE SOLVERS */ - case LinearSolver::SPGMR: - throw NewtonFailure(AMICI_NOT_IMPLEMENTED, "getSolver"); - - case LinearSolver::SPBCG: - throw NewtonFailure(AMICI_NOT_IMPLEMENTED, "getSolver"); - - case LinearSolver::SPTFQMR: - throw NewtonFailure(AMICI_NOT_IMPLEMENTED, "getSolver"); - - /* SPARSE SOLVERS */ - case LinearSolver::SuperLUMT: - throw NewtonFailure(AMICI_NOT_IMPLEMENTED, "getSolver"); - case LinearSolver::KLU: - solver.reset(new NewtonSolverSparse(model)); - break; - default: - throw NewtonFailure(AMICI_NOT_IMPLEMENTED, "getSolver"); + , dxB_(model.nJ * model.nx_solver) { + try { + // NOTE: when adding new linear solvers, make sure to also add them to + // NewtonSolver::reinitialize and NewtonSolver::is_singular + switch (linsol_type) { + case LinearSolver::dense: + linsol_.reset(new SUNLinSolDense(x_)); + break; + case LinearSolver::KLU: + linsol_.reset(new SUNLinSolKLU(x_, model.nnz, CSC_MAT)); + break; + default: + throw NewtonFailure( + AMICI_NOT_IMPLEMENTED, "Unknown linear solver type" + ); + } + } catch (NewtonFailure const&) { + throw; + } catch (AmiException const& e) { + throw NewtonFailure(AMICI_ERROR, e.what()); } - return solver; } void NewtonSolver::getStep( @@ -105,144 +80,79 @@ void NewtonSolver::computeNewtonSensis( } } -NewtonSolverDense::NewtonSolverDense(Model const& model) - : NewtonSolver(model) - , Jtmp_(model.nx_solver, model.nx_solver) - , linsol_(SUNLinSol_Dense(x_.getNVector(), Jtmp_)) { - auto status = SUNLinSolInitialize_Dense(linsol_); - if (status != SUNLS_SUCCESS) - throw NewtonFailure(status, "SUNLinSolInitialize_Dense"); -} - -void NewtonSolverDense::prepareLinearSystem( +void NewtonSolver::prepareLinearSystem( Model& model, SimulationState const& state ) { - model.fJ(state.t, 0.0, state.x, state.dx, xdot_, Jtmp_); - Jtmp_.refresh(); - auto status = SUNLinSolSetup_Dense(linsol_, Jtmp_); - if (status != SUNLS_SUCCESS) - throw NewtonFailure(status, "SUNLinSolSetup_Dense"); + auto& J = linsol_->getMatrix(); + if (J.matrix_id() == SUNMATRIX_SPARSE) { + model.fJSparse(state.t, 0.0, state.x, state.dx, xdot_, J); + } else { + model.fJ(state.t, 0.0, state.x, state.dx, xdot_, J); + } + J.refresh(); + try { + linsol_->setup(); + } catch (AmiException const& e) { + throw NewtonFailure(AMICI_SINGULAR_JACOBIAN, e.what()); + } } -void NewtonSolverDense::prepareLinearSystemB( +void NewtonSolver::prepareLinearSystemB( Model& model, SimulationState const& state ) { - model.fJB(state.t, 0.0, state.x, state.dx, xB_, dxB_, xdot_, Jtmp_); - Jtmp_.refresh(); - auto status = SUNLinSolSetup_Dense(linsol_, Jtmp_); - if (status != SUNLS_SUCCESS) - throw NewtonFailure(status, "SUNLinSolSetup_Dense"); + auto& J = linsol_->getMatrix(); + if (J.matrix_id() == SUNMATRIX_SPARSE) { + model.fJSparseB(state.t, 0.0, state.x, state.dx, xB_, dxB_, xdot_, J); + } else { + model.fJB(state.t, 0.0, state.x, state.dx, xB_, dxB_, xdot_, J); + } + J.refresh(); + try { + linsol_->setup(); + } catch (AmiException const& e) { + throw NewtonFailure(AMICI_SINGULAR_JACOBIAN, e.what()); + } } -void NewtonSolverDense::solveLinearSystem(AmiVector& rhs) { - auto status = SUNLinSolSolve_Dense( - linsol_, Jtmp_, rhs.getNVector(), rhs.getNVector(), 0.0 - ); - Jtmp_.refresh(); +void NewtonSolver::solveLinearSystem(AmiVector& rhs) { // last argument is tolerance and does not have any influence on result - + auto status = linsol_->solve(rhs.getNVector(), rhs.getNVector(), 0.0); if (status != SUNLS_SUCCESS) - throw NewtonFailure(status, "SUNLinSolSolve_Dense"); + throw NewtonFailure(status, "SUNLinSolSolve"); } -void NewtonSolverDense::reinitialize() { - /* dense solver does not need reinitialization */ +void NewtonSolver::reinitialize() { + // direct solvers do not need reinitialization + if (dynamic_cast(linsol_.get())) { + return; + } + if (auto s = dynamic_cast(linsol_.get())) { + try { + s->reInit(s->getMatrix().capacity(), SUNKLU_REINIT_PARTIAL); + return; + } catch (AmiException const&) { + throw NewtonFailure( + AMICI_ERROR, "Failed to reinitialize KLU solver" + ); + } + } + throw AmiException( + "Unhandled linear solver type in NewtonSolver::reinitialize." + ); }; -bool NewtonSolverDense::is_singular(Model& model, SimulationState const& state) +bool NewtonSolver::is_singular(Model& model, SimulationState const& state) const { + if (auto s = dynamic_cast(linsol_.get())) { + return s->is_singular(); + } + // dense solver doesn't have any implementation for rcond/condest, so use // sparse solver interface, not the most efficient solution, but who is // concerned about speed and used the dense solver anyways ¯\_(ツ)_/¯ - NewtonSolverSparse sparse_solver(model); + NewtonSolver sparse_solver(model, LinearSolver::KLU); sparse_solver.prepareLinearSystem(model, state); return sparse_solver.is_singular(model, state); } -NewtonSolverDense::~NewtonSolverDense() { - if (linsol_) - SUNLinSolFree_Dense(linsol_); -} - -NewtonSolverSparse::NewtonSolverSparse(Model const& model) - : NewtonSolver(model) - , Jtmp_(model.nx_solver, model.nx_solver, model.nnz, CSC_MAT) - , linsol_(SUNKLU(x_.getNVector(), Jtmp_)) { - auto status = SUNLinSolInitialize_KLU(linsol_); - if (status != SUNLS_SUCCESS) - throw NewtonFailure(status, "SUNLinSolInitialize_KLU"); -} - -void NewtonSolverSparse::prepareLinearSystem( - Model& model, SimulationState const& state -) { - /* Get sparse Jacobian */ - model.fJSparse(state.t, 0.0, state.x, state.dx, xdot_, Jtmp_); - Jtmp_.refresh(); - auto status = SUNLinSolSetup_KLU(linsol_, Jtmp_); - if (status != SUNLS_SUCCESS) - throw NewtonFailure(status, "SUNLinSolSetup_KLU"); -} - -void NewtonSolverSparse::prepareLinearSystemB( - Model& model, SimulationState const& state -) { - /* Get sparse Jacobian */ - model.fJSparseB(state.t, 0.0, state.x, state.dx, xB_, dxB_, xdot_, Jtmp_); - Jtmp_.refresh(); - auto status = SUNLinSolSetup_KLU(linsol_, Jtmp_); - if (status != SUNLS_SUCCESS) - throw NewtonFailure(status, "SUNLinSolSetup_KLU"); -} - -void NewtonSolverSparse::solveLinearSystem(AmiVector& rhs) { - /* Pass pointer to the linear solver */ - auto status = SUNLinSolSolve_KLU( - linsol_, Jtmp_, rhs.getNVector(), rhs.getNVector(), 0.0 - ); - // last argument is tolerance and does not have any influence on result - - if (status != SUNLS_SUCCESS) - throw NewtonFailure(status, "SUNLinSolSolve_KLU"); -} - -void NewtonSolverSparse::reinitialize() { - /* partial reinitialization, don't need to reallocate Jtmp_ */ - auto status = SUNLinSol_KLUReInit( - linsol_, Jtmp_, Jtmp_.capacity(), SUNKLU_REINIT_PARTIAL - ); - if (status != SUNLS_SUCCESS) - throw NewtonFailure(status, "SUNLinSol_KLUReInit"); -} - -bool NewtonSolverSparse:: - is_singular(Model& /*model*/, SimulationState const& /*state*/) const { - // adapted from SUNLinSolSetup_KLU in sunlinsol/klu/sunlinsol_klu.c - auto content = (SUNLinearSolverContent_KLU)(linsol_->content); - // first cheap check via rcond - auto status - = sun_klu_rcond(content->symbolic, content->numeric, &content->common); - if (status == 0) - throw NewtonFailure(content->last_flag, "sun_klu_rcond"); - - auto precision = std::numeric_limits::epsilon(); - - if (content->common.rcond < precision) { - // cheap check indicates singular, expensive check via condest - status = sun_klu_condest( - SM_INDEXPTRS_S(Jtmp_.get()), SM_DATA_S(Jtmp_.get()), - content->symbolic, content->numeric, &content->common - ); - if (status == 0) - throw NewtonFailure(content->last_flag, "sun_klu_rcond"); - return content->common.condest > 1.0 / precision; - } - return false; -} - -NewtonSolverSparse::~NewtonSolverSparse() { - if (linsol_) - SUNLinSolFree_KLU(linsol_); -} - } // namespace amici diff --git a/src/rdata.cpp b/src/rdata.cpp index 649b420ccb..6d016298bd 100644 --- a/src/rdata.cpp +++ b/src/rdata.cpp @@ -46,10 +46,6 @@ ReturnData::ReturnData( , rdata_reporting(rdrm) , sigma_res(sigma_res) , sigma_offset(sigma_offset) - , x_solver_(nx_solver) - , sx_solver_(nx_solver, nplist) - , x_rdata_(nx) - , sx_rdata_(nx, nplist) , nroots_(ne) { switch (rdata_reporting) { case RDataReporting::full: @@ -194,16 +190,14 @@ void ReturnData::processSimulationObjects( void ReturnData::processPreEquilibration( SteadystateProblem const& preeq, Model& model ) { - readSimulationState(preeq.getFinalSimulationState(), model); + auto simulation_state = preeq.getFinalSimulationState(); + model.setModelState(simulation_state.state); if (!x_ss.empty()) { - model.fx_rdata(x_rdata_, x_solver_); - writeSlice(x_rdata_, x_ss); + model.fx_rdata(x_ss, simulation_state.x); } if (!sx_ss.empty() && sensi >= SensitivityOrder::first) { - model.fsx_rdata(sx_rdata_, sx_solver_, x_solver_); - for (int ip = 0; ip < nplist; ip++) - writeSlice(sx_rdata_[ip], slice(sx_ss, ip, nx)); + model.fsx_rdata(sx_ss, simulation_state.sx, simulation_state.x); } /* Get cpu time for pre-equilibration in milliseconds */ preeq_cpu_time = preeq.getCPUTime(); @@ -223,8 +217,9 @@ void ReturnData::processPostEquilibration( for (int it = 0; it < nt; it++) { auto t = model.getTimepoint(it); if (std::isinf(t)) { - readSimulationState(posteq.getFinalSimulationState(), model); - getDataOutput(it, model, edata); + auto simulation_state = posteq.getFinalSimulationState(); + model.setModelState(simulation_state.state); + getDataOutput(it, model, simulation_state, edata); } } /* Get cpu time for Newton solve in milliseconds */ @@ -248,25 +243,24 @@ void ReturnData::processForwardProblem( auto initialState = fwd.getInitialSimulationState(); if (initialState.x.getLength() == 0 && model.nx_solver > 0) return; // if x wasn't set forward problem failed during initialization - readSimulationState(initialState, model); + + model.setModelState(initialState.state); if (!x0.empty()) { - model.fx_rdata(x_rdata_, x_solver_); - writeSlice(x_rdata_, x0); + model.fx_rdata(x0, initialState.x); } if (!sx0.empty()) { - model.fsx_rdata(sx_rdata_, sx_solver_, x_solver_); - for (int ip = 0; ip < nplist; ip++) - writeSlice(sx_rdata_[ip], slice(sx0, ip, nx)); + model.fsx_rdata(sx0, initialState.sx, initialState.x); } // process timepoint data realtype tf = fwd.getFinalTime(); for (int it = 0; it < model.nt(); it++) { if (model.getTimepoint(it) <= tf) { - readSimulationState(fwd.getSimulationStateTimepoint(it), model); - getDataOutput(it, model, edata); + auto simulation_state = fwd.getSimulationStateTimepoint(it); + model.setModelState(simulation_state.state); + getDataOutput(it, model, simulation_state, edata); } else { // check for integration failure but consider postequilibration if (!std::isinf(model.getTimepoint(it))) @@ -278,38 +272,44 @@ void ReturnData::processForwardProblem( if (nz > 0) { auto rootidx = fwd.getRootIndexes(); for (int iroot = 0; iroot <= fwd.getEventCounter(); iroot++) { - readSimulationState(fwd.getSimulationStateEvent(iroot), model); - getEventOutput(t_, rootidx.at(iroot), model, edata); + auto simulation_state = fwd.getSimulationStateEvent(iroot); + model.setModelState(simulation_state.state); + getEventOutput( + simulation_state.t, rootidx.at(iroot), model, simulation_state, + edata + ); } } } -void ReturnData::getDataOutput(int it, Model& model, ExpData const* edata) { +void ReturnData::getDataOutput( + int it, Model& model, SimulationState const& simulation_state, + ExpData const* edata +) { if (!x.empty()) { - model.fx_rdata(x_rdata_, x_solver_); - writeSlice(x_rdata_, slice(x, it, nx)); + model.fx_rdata(slice(x, it, nx), simulation_state.x); } if (!w.empty()) - model.getExpression(slice(w, it, nw), ts[it], x_solver_); + model.getExpression(slice(w, it, nw), ts[it], simulation_state.x); if (!y.empty()) - model.getObservable(slice(y, it, ny), ts[it], x_solver_); + model.getObservable(slice(y, it, ny), ts[it], simulation_state.x); if (!sigmay.empty()) model.getObservableSigma(slice(sigmay, it, ny), it, edata); if (edata) { if (!isNaN(llh)) - model.addObservableObjective(llh, it, x_solver_, *edata); - fres(it, model, *edata); + model.addObservableObjective(llh, it, simulation_state.x, *edata); + fres(it, model, simulation_state, *edata); fchi2(it, *edata); } if (sensi >= SensitivityOrder::first && nplist > 0) { if (sensi_meth == SensitivityMethod::forward) { - getDataSensisFSA(it, model, edata); + getDataSensisFSA(it, model, simulation_state, edata); } else if (edata && !sllh.empty()) { model.addPartialObservableObjectiveSensitivity( - sllh, s2llh, it, x_solver_, *edata + sllh, s2llh, it, simulation_state.x, *edata ); } @@ -321,32 +321,37 @@ void ReturnData::getDataOutput(int it, Model& model, ExpData const* edata) { } } -void ReturnData::getDataSensisFSA(int it, Model& model, ExpData const* edata) { +void ReturnData::getDataSensisFSA( + int it, Model& model, SimulationState const& simulation_state, + ExpData const* edata +) { if (!sx.empty()) { - model.fsx_rdata(sx_rdata_, sx_solver_, x_solver_); - for (int ip = 0; ip < nplist; ip++) { - writeSlice(sx_rdata_[ip], slice(sx, it * nplist + ip, nx)); - } + model.fsx_rdata( + gsl::span(&sx.at(it * nplist * nx), nplist * nx), + simulation_state.sx, simulation_state.x + ); } if (!sy.empty()) { model.getObservableSensitivity( - slice(sy, it, nplist * ny), ts[it], x_solver_, sx_solver_ + slice(sy, it, nplist * ny), ts[it], simulation_state.x, + simulation_state.sx ); } if (edata) { if (!sllh.empty()) model.addObservableObjectiveSensitivity( - sllh, s2llh, it, x_solver_, sx_solver_, *edata + sllh, s2llh, it, simulation_state.x, simulation_state.sx, *edata ); - fsres(it, model, *edata); - fFIM(it, model, *edata); + fsres(it, model, simulation_state, *edata); + fFIM(it, model, simulation_state, *edata); } } void ReturnData::getEventOutput( - realtype t, std::vector rootidx, Model& model, ExpData const* edata + realtype t, std::vector const& rootidx, Model& model, + SimulationState const& simulation_state, ExpData const* edata ) { for (int ie = 0; ie < ne; ie++) { @@ -355,13 +360,15 @@ void ReturnData::getEventOutput( /* get event output */ if (!z.empty()) - model.getEvent(slice(z, nroots_.at(ie), nz), ie, t, x_solver_); + model.getEvent( + slice(z, nroots_.at(ie), nz), ie, t, simulation_state.x + ); /* if called from fillEvent at last timepoint, then also get the root function value */ if (t == model.getTimepoint(nt - 1)) if (!rz.empty()) model.getEventRegularization( - slice(rz, nroots_.at(ie), nz), ie, t, x_solver_ + slice(rz, nroots_.at(ie), nz), ie, t, simulation_state.x ); if (edata) { @@ -372,24 +379,25 @@ void ReturnData::getEventOutput( ); if (!isNaN(llh)) model.addEventObjective( - llh, ie, nroots_.at(ie), t, x_solver_, *edata + llh, ie, nroots_.at(ie), t, simulation_state.x, *edata ); /* if called from fillEvent at last timepoint, add regularization based on rz */ if (t == model.getTimepoint(nt - 1) && !isNaN(llh)) { model.addEventObjectiveRegularization( - llh, ie, nroots_.at(ie), t, x_solver_, *edata + llh, ie, nroots_.at(ie), t, simulation_state.x, *edata ); } } if (sensi >= SensitivityOrder::first) { if (sensi_meth == SensitivityMethod::forward) { - getEventSensisFSA(ie, t, model, edata); + getEventSensisFSA(ie, t, model, simulation_state, edata); } else if (edata && !sllh.empty()) { model.addPartialEventObjectiveSensitivity( - sllh, s2llh, ie, nroots_.at(ie), t, x_solver_, *edata + sllh, s2llh, ie, nroots_.at(ie), t, simulation_state.x, + *edata ); } } @@ -398,7 +406,8 @@ void ReturnData::getEventOutput( } void ReturnData::getEventSensisFSA( - int ie, realtype t, Model& model, ExpData const* edata + int ie, realtype t, Model& model, SimulationState const& simulation_state, + ExpData const* edata ) { if (t == model.getTimepoint(nt - 1)) { // call from fillEvent at last timepoint @@ -408,18 +417,20 @@ void ReturnData::getEventSensisFSA( ); if (!srz.empty()) model.getEventRegularizationSensitivity( - slice(srz, nroots_.at(ie), nz * nplist), ie, t, x_solver_, - sx_solver_ + slice(srz, nroots_.at(ie), nz * nplist), ie, t, + simulation_state.x, simulation_state.sx ); } else if (!sz.empty()) { model.getEventSensitivity( - slice(sz, nroots_.at(ie), nz * nplist), ie, t, x_solver_, sx_solver_ + slice(sz, nroots_.at(ie), nz * nplist), ie, t, simulation_state.x, + simulation_state.sx ); } if (edata && !sllh.empty()) { model.addEventObjectiveSensitivity( - sllh, s2llh, ie, nroots_.at(ie), t, x_solver_, sx_solver_, *edata + sllh, s2llh, ie, nroots_.at(ie), t, simulation_state.x, + simulation_state.sx, *edata ); } } @@ -430,7 +441,9 @@ void ReturnData::processBackwardProblem( ) { if (sllh.empty()) return; - readSimulationState(fwd.getInitialSimulationState(), model); + + auto simulation_state = fwd.getInitialSimulationState(); + model.setModelState(simulation_state.state); std::vector llhS0(model.nJ * model.nplist(), 0.0); auto xB = bwd.getAdjointState(); @@ -439,7 +452,7 @@ void ReturnData::processBackwardProblem( if (preeq && preeq->hasQuadrature()) { handleSx0Backward(model, *preeq, llhS0, xQB); } else { - handleSx0Forward(model, llhS0, xB); + handleSx0Forward(model, simulation_state, llhS0, xB); } for (int iJ = 0; iJ < model.nJ; iJ++) { @@ -486,7 +499,8 @@ void ReturnData::handleSx0Backward( } void ReturnData::handleSx0Forward( - Model const& model, std::vector& llhS0, AmiVector& xB + Model const& model, SimulationState const& simulation_state, + std::vector& llhS0, AmiVector& xB ) const { /* If preequilibration is run in forward mode or is not needed, then adjoint sensitivity analysis still needs the state sensitivities at t=0 (sx0), @@ -497,7 +511,7 @@ void ReturnData::handleSx0Forward( for (int ip = 0; ip < model.nplist(); ++ip) { llhS0[ip] = 0.0; for (int ix = 0; ix < model.nxtrue_solver; ++ix) { - llhS0[ip] += xB[ix] * sx_solver_.at(ix, ip); + llhS0[ip] += xB[ix] * simulation_state.sx.at(ix, ip); } } } else { @@ -506,9 +520,9 @@ void ReturnData::handleSx0Forward( for (int ix = 0; ix < model.nxtrue_solver; ++ix) { llhS0[ip + iJ * model.nplist()] += xB[ix + iJ * model.nxtrue_solver] - * sx_solver_.at(ix, ip) + * simulation_state.sx.at(ix, ip) + xB[ix] - * sx_solver_.at( + * simulation_state.sx.at( ix + iJ * model.nxtrue_solver, ip ); } @@ -575,17 +589,6 @@ void ReturnData::processSolver(Solver const& solver) { } } -void ReturnData::readSimulationState( - SimulationState const& state, Model& model -) { - x_solver_ = state.x; - dx_solver_ = state.dx; - if (computingFSA() || state.t == model.t0()) - sx_solver_ = state.sx; - t_ = state.t; - model.setModelState(state.state); -} - void ReturnData::invalidate(int const it_start) { if (it_start >= nt) return; @@ -817,12 +820,15 @@ static realtype fres_error(realtype sigma_y, realtype sigma_offset) { return sqrt(2 * log(sigma_y) + sigma_offset); } -void ReturnData::fres(int const it, Model& model, ExpData const& edata) { +void ReturnData::fres( + int const it, Model& model, SimulationState const& simulation_state, + ExpData const& edata +) { if (res.empty()) return; std::vector y_it(ny, 0.0); - model.getObservable(y_it, ts[it], x_solver_); + model.getObservable(y_it, ts[it], simulation_state.x); std::vector sigmay_it(ny, 0.0); model.getObservableSigma(sigmay_it, it, &edata); @@ -877,14 +883,19 @@ fsres_error(realtype sigma_y, realtype ssigma_y, realtype sigma_offset) { return ssigma_y / (fres_error(sigma_y, sigma_offset) * sigma_y); } -void ReturnData::fsres(int const it, Model& model, ExpData const& edata) { +void ReturnData::fsres( + int const it, Model& model, SimulationState const& simulation_state, + ExpData const& edata +) { if (sres.empty()) return; std::vector y_it(ny, 0.0); - model.getObservable(y_it, ts[it], x_solver_); + model.getObservable(y_it, ts[it], simulation_state.x); std::vector sy_it(ny * nplist, 0.0); - model.getObservableSensitivity(sy_it, ts[it], x_solver_, sx_solver_); + model.getObservableSensitivity( + sy_it, ts[it], simulation_state.x, simulation_state.sx + ); std::vector sigmay_it(ny, 0.0); model.getObservableSigma(sigmay_it, it, &edata); @@ -917,14 +928,19 @@ void ReturnData::fsres(int const it, Model& model, ExpData const& edata) { } } -void ReturnData::fFIM(int it, Model& model, ExpData const& edata) { +void ReturnData::fFIM( + int it, Model& model, SimulationState const& simulation_state, + ExpData const& edata +) { if (FIM.empty()) return; std::vector y_it(ny, 0.0); - model.getObservable(y_it, ts[it], x_solver_); + model.getObservable(y_it, ts[it], simulation_state.x); std::vector sy_it(ny * nplist, 0.0); - model.getObservableSensitivity(sy_it, ts[it], x_solver_, sx_solver_); + model.getObservableSensitivity( + sy_it, ts[it], simulation_state.x, simulation_state.sx + ); std::vector sigmay_it(ny, 0.0); model.getObservableSigma(sigmay_it, it, &edata); diff --git a/src/steadystateproblem.cpp b/src/steadystateproblem.cpp index d78f9a8705..663d83beba 100644 --- a/src/steadystateproblem.cpp +++ b/src/steadystateproblem.cpp @@ -28,6 +28,7 @@ SteadystateProblem::SteadystateProblem(Solver const& solver, Model const& model) , xQ_(model.nJ * model.nx_solver) , xQB_(model.nplist()) , xQBdot_(model.nplist()) + , steadystate_mask_(AmiVector(model.get_steadystate_mask())) , max_steps_(solver.getNewtonMaxSteps()) , dJydx_(model.nJ * model.nx_solver * model.nt(), 0.0) , state_( @@ -44,7 +45,7 @@ SteadystateProblem::SteadystateProblem(Solver const& solver, Model const& model) , rtol_sensi_(solver.getRelativeToleranceSteadyStateSensi()) , atol_quad_(solver.getAbsoluteToleranceQuadratures()) , rtol_quad_(solver.getRelativeToleranceQuadratures()) - , newton_solver_(NewtonSolver::getSolver(solver, model)) + , newton_solver_(NewtonSolver(model, solver.getLinearSolver())) , damping_factor_mode_(solver.getNewtonDampingFactorMode()) , damping_factor_lower_bound_(solver.getNewtonDampingFactorLowerBound()) , newton_step_conv_(solver.getNewtonStepSteadyStateCheck()) @@ -86,7 +87,7 @@ void SteadystateProblem::workSteadyStateProblem( try { /* this might still fail, if the Jacobian is singular and simulation did not find a steady state */ - newton_solver_->computeNewtonSensis(state_.sx, model, state_); + newton_solver_.computeNewtonSensis(state_.sx, model, state_); } catch (NewtonFailure const&) { throw AmiException( "Steady state sensitivity computation failed due " @@ -249,7 +250,7 @@ void SteadystateProblem::findSteadyStateBySimulation( void SteadystateProblem::initializeForwardProblem( int it, Solver const& solver, Model& model ) { - newton_solver_->reinitialize(); + newton_solver_.reinitialize(); /* process solver handling for pre- or postequilibration */ if (it == -1) { /* solver was not run before, set up everything */ @@ -278,7 +279,7 @@ void SteadystateProblem::initializeForwardProblem( bool SteadystateProblem::initializeBackwardProblem( Solver const& solver, Model& model, BackwardProblem const* bwd ) { - newton_solver_->reinitialize(); + newton_solver_.reinitialize(); /* note that state_ is still set from forward run */ if (bwd) { /* preequilibration */ @@ -358,8 +359,8 @@ void SteadystateProblem::getQuadratureByLinSolve(Model& model) { /* try to solve the linear system */ try { /* compute integral over xB and write to xQ */ - newton_solver_->prepareLinearSystemB(model, state_); - newton_solver_->solveLinearSystem(xQ_); + newton_solver_.prepareLinearSystemB(model, state_); + newton_solver_.solveLinearSystem(xQ_); /* Compute the quadrature as the inner product xQ * dxdotdp */ computeQBfromQ(model, xQ_, xQB_); /* set flag that quadratures is available (for processing in rdata) */ @@ -551,8 +552,7 @@ SteadystateProblem::getWrms(Model& model, SensitivityMethod sensi_method) { "steady state computations. Stopping." ); wrms = getWrmsNorm( - xQB_, xQBdot_, model.get_steadystate_mask_av(), atol_quad_, - rtol_quad_, ewtQB_ + xQB_, xQBdot_, steadystate_mask_, atol_quad_, rtol_quad_, ewtQB_ ); } else { /* If we're doing a forward simulation (with or without sensitivities: @@ -562,8 +562,8 @@ SteadystateProblem::getWrms(Model& model, SensitivityMethod sensi_method) { else updateRightHandSide(model); wrms = getWrmsNorm( - state_.x, newton_step_conv_ ? delta_ : xdot_, - model.get_steadystate_mask_av(), atol_, rtol_, ewt_ + state_.x, newton_step_conv_ ? delta_ : xdot_, steadystate_mask_, + atol_, rtol_, ewt_ ); } return wrms; @@ -583,10 +583,10 @@ realtype SteadystateProblem::getWrmsFSA(Model& model) { state_.t, state_.x, state_.dx, ip, state_.sx[ip], state_.dx, xdot_ ); if (newton_step_conv_) - newton_solver_->solveLinearSystem(xdot_); + newton_solver_.solveLinearSystem(xdot_); wrms = getWrmsNorm( - state_.sx[ip], xdot_, model.get_steadystate_mask_av(), atol_sensi_, - rtol_sensi_, ewt_ + state_.sx[ip], xdot_, steadystate_mask_, atol_sensi_, rtol_sensi_, + ewt_ ); /* ideally this function would report the maximum of all wrms over all ip, but for practical purposes we can just report the wrms for @@ -887,7 +887,7 @@ void SteadystateProblem::getNewtonStep(Model& model) { return; updateRightHandSide(model); delta_.copy(xdot_); - newton_solver_->getStep(delta_, model, state_); + newton_solver_.getStep(delta_, model, state_); delta_updated_ = true; } } // namespace amici diff --git a/src/sundials_linsol_wrapper.cpp b/src/sundials_linsol_wrapper.cpp index 6836949cb8..b94c28f5aa 100644 --- a/src/sundials_linsol_wrapper.cpp +++ b/src/sundials_linsol_wrapper.cpp @@ -9,6 +9,12 @@ namespace amici { SUNLinSolWrapper::SUNLinSolWrapper(SUNLinearSolver linsol) : solver_(linsol) {} +SUNLinSolWrapper::SUNLinSolWrapper( + SUNLinearSolver linsol, SUNMatrixWrapper const& A +) + : solver_(linsol) + , A_(A) {} + SUNLinSolWrapper::~SUNLinSolWrapper() { if (solver_) SUNLinSolFree(solver_); @@ -16,6 +22,14 @@ SUNLinSolWrapper::~SUNLinSolWrapper() { SUNLinSolWrapper::SUNLinSolWrapper(SUNLinSolWrapper&& other) noexcept { std::swap(solver_, other.solver_); + std::swap(A_, other.A_); +} + +SUNLinSolWrapper& SUNLinSolWrapper::operator=(SUNLinSolWrapper&& other +) noexcept { + std::swap(solver_, other.solver_); + std::swap(A_, other.A_); + return *this; } SUNLinearSolver SUNLinSolWrapper::get() const { return solver_; } @@ -27,23 +41,20 @@ SUNLinearSolver_Type SUNLinSolWrapper::getType() const { int SUNLinSolWrapper::initialize() { auto res = SUNLinSolInitialize(solver_); if (res != SUNLS_SUCCESS) - throw AmiException("Solver initialization failed with code %d", res); + throw AmiException( + "Linear solver initialization failed with code %d", res + ); return res; } -void SUNLinSolWrapper::setup(SUNMatrix A) const { - auto res = SUNLinSolSetup(solver_, A); +void SUNLinSolWrapper::setup() const { + auto res = SUNLinSolSetup(solver_, A_.get()); if (res != SUNLS_SUCCESS) - throw AmiException("Solver setup failed with code %d", res); -} - -void SUNLinSolWrapper::setup(SUNMatrixWrapper const& A) const { - return setup(A.get()); + throw AmiException("Linear solver setup failed with code %d", res); } -int SUNLinSolWrapper::Solve(SUNMatrix A, N_Vector x, N_Vector b, realtype tol) - const { - return SUNLinSolSolve(solver_, A, x, b, tol); +int SUNLinSolWrapper::solve(N_Vector x, N_Vector b, realtype tol) const { + return SUNLinSolSolve(solver_, A_.get(), x, b, tol); } long SUNLinSolWrapper::getLastFlag() const { @@ -54,7 +65,7 @@ int SUNLinSolWrapper::space(long* lenrwLS, long* leniwLS) const { return SUNLinSolSpace(solver_, lenrwLS, leniwLS); } -SUNMatrix SUNLinSolWrapper::getMatrix() const { return nullptr; } +SUNMatrixWrapper& SUNLinSolWrapper::getMatrix() { return A_; } SUNNonLinSolWrapper::SUNNonLinSolWrapper(SUNNonlinearSolver sol) : solver(sol) {} @@ -153,31 +164,29 @@ void SUNNonLinSolWrapper::initialize() { ); } -SUNLinSolBand::SUNLinSolBand(N_Vector x, SUNMatrix A) +SUNLinSolBand::SUNLinSolBand(N_Vector x, SUNMatrixWrapper A) : SUNLinSolWrapper(SUNLinSol_Band(x, A)) { if (!solver_) throw AmiException("Failed to create solver."); } SUNLinSolBand::SUNLinSolBand(AmiVector const& x, int ubw, int lbw) - : A_(SUNMatrixWrapper(x.getLength(), ubw, lbw)) { + : SUNLinSolWrapper(nullptr, SUNMatrixWrapper(x.getLength(), ubw, lbw)) { solver_ = SUNLinSol_Band(const_cast(x.getNVector()), A_); if (!solver_) throw AmiException("Failed to create solver."); } -SUNMatrix SUNLinSolBand::getMatrix() const { return A_.get(); } - SUNLinSolDense::SUNLinSolDense(AmiVector const& x) - : A_(SUNMatrixWrapper(x.getLength(), x.getLength())) { + : SUNLinSolWrapper( + nullptr, SUNMatrixWrapper(x.getLength(), x.getLength()) + ) { solver_ = SUNLinSol_Dense(const_cast(x.getNVector()), A_); if (!solver_) throw AmiException("Failed to create solver."); } -SUNMatrix SUNLinSolDense::getMatrix() const { return A_.get(); } - -SUNLinSolKLU::SUNLinSolKLU(N_Vector x, SUNMatrix A) +SUNLinSolKLU::SUNLinSolKLU(N_Vector x, SUNMatrixWrapper A) : SUNLinSolWrapper(SUNLinSol_KLU(x, A)) { if (!solver_) throw AmiException("Failed to create solver."); @@ -186,7 +195,10 @@ SUNLinSolKLU::SUNLinSolKLU(N_Vector x, SUNMatrix A) SUNLinSolKLU::SUNLinSolKLU( AmiVector const& x, int nnz, int sparsetype, StateOrdering ordering ) - : A_(SUNMatrixWrapper(x.getLength(), x.getLength(), nnz, sparsetype)) { + : SUNLinSolWrapper( + nullptr, + SUNMatrixWrapper(x.getLength(), x.getLength(), nnz, sparsetype) + ) { solver_ = SUNLinSol_KLU(const_cast(x.getNVector()), A_); if (!solver_) throw AmiException("Failed to create solver."); @@ -194,8 +206,6 @@ SUNLinSolKLU::SUNLinSolKLU( setOrdering(ordering); } -SUNMatrix SUNLinSolKLU::getMatrix() const { return A_.get(); } - void SUNLinSolKLU::reInit(int nnz, int reinit_type) { int status = SUNLinSol_KLUReInit(solver_, A_, nnz, reinit_type); if (status != SUNLS_SUCCESS) @@ -208,6 +218,30 @@ void SUNLinSolKLU::setOrdering(StateOrdering ordering) { throw AmiException("SUNLinSol_KLUSetOrdering failed with %d", status); } +bool SUNLinSolKLU::is_singular() const { + // adapted from SUNLinSolSetup_KLU in sunlinsol/klu/sunlinsol_klu.c + auto content = (SUNLinearSolverContent_KLU)(solver_->content); + // first cheap check via rcond + auto status + = sun_klu_rcond(content->symbolic, content->numeric, &content->common); + if (status == 0) + throw AmiException("sun_klu_rcond: %d", content->last_flag); + + auto precision = std::numeric_limits::epsilon(); + + if (content->common.rcond < precision) { + // cheap check indicates singular, expensive check via condest + status = sun_klu_condest( + SM_INDEXPTRS_S(A_.get()), SM_DATA_S(A_.get()), content->symbolic, + content->numeric, &content->common + ); + if (status == 0) + throw AmiException("sun_klu_condest: %d", content->last_flag); + return content->common.condest > 1.0 / precision; + } + return false; +} + SUNLinSolPCG::SUNLinSolPCG(N_Vector y, int pretype, int maxl) : SUNLinSolWrapper(SUNLinSol_PCG(y, pretype, maxl)) { if (!solver_) @@ -413,8 +447,10 @@ int SUNNonLinSolFixedPoint::getSysFn(SUNNonlinSolSysFn* SysFn) const { #ifdef SUNDIALS_SUPERLUMT -SUNLinSolSuperLUMT::SUNLinSolSuperLUMT(N_Vector x, SUNMatrix A, int numThreads) - : SUNLinSolWrapper(SUNLinSol_SuperLUMT(x, A, numThreads)) { +SUNLinSolSuperLUMT::SUNLinSolSuperLUMT( + N_Vector x, SUNMatrixWrapper A, int numThreads +) + : SUNLinSolWrapper(SUNLinSol_SuperLUMT(x, A, numThreads), A) { if (!solver) throw AmiException("Failed to create solver."); } @@ -423,7 +459,10 @@ SUNLinSolSuperLUMT::SUNLinSolSuperLUMT( AmiVector const& x, int nnz, int sparsetype, SUNLinSolSuperLUMT::StateOrdering ordering ) - : A(SUNMatrixWrapper(x.getLength(), x.getLength(), nnz, sparsetype)) { + : SUNLinSolWrapper( + nullptr, + SUNMatrixWrapper(x.getLength(), x.getLength(), nnz, sparsetype) + ) { int numThreads = 1; if (auto env = std::getenv("AMICI_SUPERLUMT_NUM_THREADS")) { numThreads = std::max(1, std::stoi(env)); @@ -440,7 +479,10 @@ SUNLinSolSuperLUMT::SUNLinSolSuperLUMT( AmiVector const& x, int nnz, int sparsetype, StateOrdering ordering, int numThreads ) - : A(SUNMatrixWrapper(x.getLength(), x.getLength(), nnz, sparsetype)) { + : SUNLinSolWrapper( + nullptr, + SUNMatrixWrapper(x.getLength(), x.getLength(), nnz, sparsetype) + ) { solver = SUNLinSol_SuperLUMT(x.getNVector(), A.get(), numThreads); if (!solver) throw AmiException("Failed to create solver."); @@ -448,8 +490,6 @@ SUNLinSolSuperLUMT::SUNLinSolSuperLUMT( setOrdering(ordering); } -SUNMatrix SUNLinSolSuperLUMT::getMatrix() const { return A.get(); } - void SUNLinSolSuperLUMT::setOrdering(StateOrdering ordering) { auto status = SUNLinSol_SuperLUMTSetOrdering(solver, static_cast(ordering)); diff --git a/tests/benchmark-models/test_petab_benchmark.py b/tests/benchmark-models/test_petab_benchmark.py index 82d27f85eb..f82787494c 100644 --- a/tests/benchmark-models/test_petab_benchmark.py +++ b/tests/benchmark-models/test_petab_benchmark.py @@ -2,6 +2,8 @@ from pathlib import Path +import fiddy + import amici import numpy as np import pandas as pd @@ -9,11 +11,13 @@ import pytest from amici.petab.petab_import import import_petab_problem import benchmark_models_petab - - -# Absolute and relative tolerances for finite difference gradient checks. -ATOL: float = 1e-3 -RTOL: float = 1e-2 +from collections import defaultdict +from dataclasses import dataclass +from amici import SensitivityMethod +from fiddy import MethodId, get_derivative +from fiddy.derivative_check import NumpyIsCloseDerivativeCheck +from fiddy.extensions.amici import simulate_petab_to_cached_functions +from fiddy.success import Consistency repo_root = Path(__file__).parent.parent.parent @@ -30,10 +34,15 @@ "Beer_MolBioSystems2014", "Alkan_SciSignal2018", "Lang_PLOSComputBiol2024", + "Smith_BMCSystBiol2013", # excluded due to excessive numerical failures "Crauste_CellSystems2017", "Fujita_SciSignal2010", + # FIXME: re-enable + "Raia_CancerResearch2011", + "Zheng_PNAS2012", } +models = list(sorted(models)) debug = False if debug: @@ -41,19 +50,83 @@ debug_path.mkdir(exist_ok=True, parents=True) -# until fiddy is updated +@dataclass +class GradientCheckSettings: + # Absolute and relative tolerances for simulation + atol_sim: float = 1e-16 + rtol_sim: float = 1e-12 + # Absolute and relative tolerances for finite difference gradient checks. + atol_check: float = 1e-3 + rtol_check: float = 1e-2 + # Step sizes for finite difference gradient checks. + step_sizes = [ + 1e-1, + 5e-2, + 1e-2, + 1e-3, + 1e-4, + 1e-5, + ] + rng_seed: int = 0 + + +settings = defaultdict(GradientCheckSettings) +settings["Smith_BMCSystBiol2013"] = GradientCheckSettings( + atol_sim=1e-10, + rtol_sim=1e-10, +) +settings["Oliveira_NatCommun2021"] = GradientCheckSettings( + atol_sim=1e-10, + rtol_sim=1e-10, +) +settings["Okuonghae_ChaosSolitonsFractals2020"] = GradientCheckSettings( + atol_sim=1e-14, + rtol_sim=1e-14, +) +settings["Okuonghae_ChaosSolitonsFractals2020"].step_sizes.insert(0, 0.2) +settings["Giordano_Nature2020"] = GradientCheckSettings( + atol_check=1e-6, rtol_check=1e-3, rng_seed=1 +) + + +def assert_gradient_check_success( + derivative: fiddy.Derivative, + expected_derivative: np.array, + point: np.array, + atol: float, + rtol: float, +) -> None: + if not derivative.df.success.all(): + raise AssertionError( + f"Failed to compute finite differences:\n{derivative.df}" + ) + check = NumpyIsCloseDerivativeCheck( + derivative=derivative, + expectation=expected_derivative, + point=point, + ) + check_result = check(rtol=rtol, atol=atol) + + if check_result.success is True: + return + + raise AssertionError(f"Gradient check failed:\n{check_result.df}") + + @pytest.mark.filterwarnings( - "ignore:Importing amici.petab_objective is deprecated.:DeprecationWarning" + "ignore:divide by zero encountered in log", + # https://github.com/AMICI-dev/AMICI/issues/18 + "ignore:Adjoint sensitivity analysis for models with discontinuous " + "right hand sides .*:UserWarning", +) +@pytest.mark.parametrize("scale", (True, False), ids=["scaled", "unscaled"]) +@pytest.mark.parametrize( + "sensitivity_method", + (SensitivityMethod.forward, SensitivityMethod.adjoint), + ids=["forward", "adjoint"], ) -@pytest.mark.filterwarnings("ignore:divide by zero encountered in log10") -@pytest.mark.parametrize("scale", (True, False)) @pytest.mark.parametrize("model", models) -def test_benchmark_gradient(model, scale): - from fiddy import MethodId, get_derivative - from fiddy.derivative_check import NumpyIsCloseDerivativeCheck - from fiddy.extensions.amici import simulate_petab_to_cached_functions - from fiddy.success import Consistency - +def test_benchmark_gradient(model, scale, sensitivity_method, request): if not scale and model in ( "Smith_BMCSystBiol2013", "Brannmark_JBC2010", @@ -67,6 +140,28 @@ def test_benchmark_gradient(model, scale): # only fail on linear scale pytest.skip() + if ( + model + in ( + "Blasi_CellSystems2016", + "Oliveira_NatCommun2021", + ) + and sensitivity_method == SensitivityMethod.adjoint + ): + # FIXME + pytest.skip() + + if ( + model + in ( + "Weber_BMC2015", + "Sneyd_PNAS2002", + ) + and sensitivity_method == SensitivityMethod.forward + ): + # FIXME + pytest.skip() + petab_problem = benchmark_models_petab.get_problem(model) petab.flatten_timepoint_specific_output_overrides(petab_problem) @@ -75,6 +170,7 @@ def test_benchmark_gradient(model, scale): petab_problem.x_free_ids ] parameter_ids = list(parameter_df_free.index) + cur_settings = settings[model] # Setup AMICI objects. amici_model = import_petab_problem( @@ -82,18 +178,10 @@ def test_benchmark_gradient(model, scale): model_output_dir=benchmark_outdir / model, ) amici_solver = amici_model.getSolver() - amici_solver.setAbsoluteTolerance(1e-12) - amici_solver.setRelativeTolerance(1e-12) - if model in ( - "Smith_BMCSystBiol2013", - "Oliveira_NatCommun2021", - ): - amici_solver.setAbsoluteTolerance(1e-10) - amici_solver.setRelativeTolerance(1e-10) - elif model in ("Okuonghae_ChaosSolitonsFractals2020",): - amici_solver.setAbsoluteTolerance(1e-14) - amici_solver.setRelativeTolerance(1e-14) + amici_solver.setAbsoluteTolerance(cur_settings.atol_sim) + amici_solver.setRelativeTolerance(cur_settings.rtol_sim) amici_solver.setMaxSteps(int(1e5)) + amici_solver.setSensitivityMethod(sensitivity_method) if model in ("Brannmark_JBC2010",): amici_model.setSteadyStateSensitivityMode( @@ -107,78 +195,90 @@ def test_benchmark_gradient(model, scale): solver=amici_solver, scaled_parameters=scale, scaled_gradients=scale, - cache=not debug, + # FIXME: there is some issue with caching in fiddy + # e.g. Elowitz_Nature2000-True fails with cache=True, + # but not with cache=False + # cache=not debug, + cache=False, ) noise_level = 0.1 + np.random.seed(cur_settings.rng_seed) - np.random.seed(0) - if scale: - point = np.asarray( - list( - petab_problem.scale_parameters( - dict(parameter_df_free.nominalValue) - ).values() + # find a point where the derivative can be computed + for _ in range(5): + if scale: + point = np.asarray( + list( + petab_problem.scale_parameters( + dict(parameter_df_free.nominalValue) + ).values() + ) ) - ) - point_noise = np.random.randn(len(point)) * noise_level - else: - point = parameter_df_free.nominalValue.values - point_noise = np.random.randn(len(point)) * point * noise_level - point += point_noise # avoid small gradients at nominal value - - expected_derivative = amici_derivative(point) + point_noise = np.random.randn(len(point)) * noise_level + else: + point = parameter_df_free.nominalValue.values + point_noise = np.random.randn(len(point)) * point * noise_level + point += point_noise # avoid small gradients at nominal value - sizes = [ - 1e-1, - 1e-2, - 1e-3, - 1e-4, - 1e-5, - ] - if model in ("Okuonghae_ChaosSolitonsFractals2020",): - sizes.insert(0, 0.2) + try: + expected_derivative = amici_derivative(point) + break + except RuntimeError as e: + print(e) + continue + else: + raise RuntimeError("Could not compute expected derivative.") derivative = get_derivative( function=amici_function, point=point, - sizes=sizes, + sizes=cur_settings.step_sizes, direction_ids=parameter_ids, method_ids=[MethodId.CENTRAL, MethodId.FORWARD, MethodId.BACKWARD], - success_checker=Consistency(atol=1e-4, rtol=0.2), + success_checker=Consistency(atol=1e-5, rtol=1e-1), expected_result=expected_derivative, relative_sizes=not scale, ) - success = False - if derivative.df.success.all(): - check = NumpyIsCloseDerivativeCheck( - derivative=derivative, - expectation=expected_derivative, - point=point, - ) - success = check(rtol=RTOL, atol=ATOL) if debug: - df = pd.DataFrame( - [ - { - ( - "fd", - r.metadata["size_absolute"], - str(r.method_id), - ): r.value - for c in d.computers - for r in c.results - } - for d in derivative.directional_derivatives - ], - index=parameter_ids, + write_debug_output( + debug_path / f"{request.node.callspec.id}.tsv", + derivative, + expected_derivative, + parameter_ids, ) - df[("fd", "full", "")] = derivative.series.values - df[("amici", "", "")] = expected_derivative - file_name = f"{model}_scale={scale}.tsv" - df.to_csv(debug_path / file_name, sep="\t") + assert_gradient_check_success( + derivative, + expected_derivative, + point, + rtol=cur_settings.rtol_check, + atol=cur_settings.atol_check, + ) + + +def write_debug_output( + file_name, derivative, expected_derivative, parameter_ids +): + df = pd.DataFrame( + [ + { + ( + "fd", + r.metadata["size_absolute"], + str(r.method_id), + ): r.value + for c in d.computers + for r in c.results + } + for d in derivative.directional_derivatives + ], + index=parameter_ids, + ) + df[("fd", "full", "")] = derivative.series.values + df[("amici", "", "")] = expected_derivative + df["abs_diff"] = np.abs(df[("fd", "full", "")] - df[("amici", "", "")]) + df["rel_diff"] = df["abs_diff"] / np.abs(df[("amici", "", "")]) - # The gradients for all parameters are correct. - assert success, derivative.df + df.to_csv(file_name, sep="\t") diff --git a/tests/cpp/CMakeLists.txt b/tests/cpp/CMakeLists.txt index cf8ddece71..7d0096ae39 100644 --- a/tests/cpp/CMakeLists.txt +++ b/tests/cpp/CMakeLists.txt @@ -93,6 +93,7 @@ foreach(MODEL IN ITEMS ${TEST_MODELS}) TEST_COMMAND "" BUILD_ALWAYS 1 DEPENDS amici + CMAKE_ARGS "-DCMAKE_BUILD_TYPE=Debug" BUILD_BYPRODUCTS "${MODEL_LIBRARY}") # Rebuild if amici files are updated ExternalProject_Add_StepDependencies(external_model_${MODEL} build amici) diff --git a/tests/performance/test.py b/tests/performance/test.py index f117371910..677e6e647d 100755 --- a/tests/performance/test.py +++ b/tests/performance/test.py @@ -8,8 +8,8 @@ from pathlib import Path import amici -import petab -from amici.petab_import import import_model +import petab.v1 as petab +from amici.petab.petab_import import import_model def parse_args(): diff --git a/version.txt b/version.txt index 894542aa24..3f45a6442e 100644 --- a/version.txt +++ b/version.txt @@ -1 +1 @@ -0.26.2 +0.26.3