diff --git a/.github/workflows/pr.yaml b/.github/workflows/pr.yaml
index 8467d20..e9c3c90 100644
--- a/.github/workflows/pr.yaml
+++ b/.github/workflows/pr.yaml
@@ -81,6 +81,35 @@ jobs:
--python 3.14 \
--no-anaconda-upload
+ pre-commit:
+ name: Pre-commit
+ needs: unit
+ runs-on: ubuntu-26.04
+ timeout-minutes: 15
+
+ steps:
+ - name: Checkout
+ uses: actions/checkout@8e8c483db84b4bee98b60c0593521ed34d9990e8
+
+ - name: Set up Python 3.14
+ uses: actions/setup-python@a309ff8b426b58ec0e2a45f0f869d46889d02405
+ with:
+ python-version: "3.14"
+
+ - name: Install pre-commit dependencies
+ run: |
+ python -m pip install --upgrade pip
+ python -m pip install -e .[pre-commit]
+
+ - name: Run pre-commit
+ run: |
+ pre-commit install
+ pre-commit run --all-files || {
+ git status --short
+ git diff
+ exit 1
+ }
+
docs:
name: Docs
needs: unit
@@ -154,4 +183,4 @@ jobs:
with:
github-token: ${{ secrets.GITHUB_TOKEN }}
file: coverage.xml
- fail-on-error: false
\ No newline at end of file
+ fail-on-error: false
diff --git a/.github/workflows/weekly-regression.yaml b/.github/workflows/weekly-regression.yaml
index 167870a..8f6f3c1 100644
--- a/.github/workflows/weekly-regression.yaml
+++ b/.github/workflows/weekly-regression.yaml
@@ -53,4 +53,4 @@ jobs:
-c conda-forge \
-c bioconda \
--python ${{ matrix.python-version }} \
- --no-anaconda-upload
\ No newline at end of file
+ --no-anaconda-upload
diff --git a/.gitignore b/.gitignore
index 4753c57..654b986 100644
--- a/.gitignore
+++ b/.gitignore
@@ -86,4 +86,4 @@ coverage.xml
# Regression test artefacts
.testdata/
-test-results/
\ No newline at end of file
+test-results/
diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml
new file mode 100644
index 0000000..427ad86
--- /dev/null
+++ b/.pre-commit-config.yaml
@@ -0,0 +1,46 @@
+repos:
+ - repo: https://github.com/astral-sh/ruff-pre-commit
+ rev: v0.15.2
+ hooks:
+ - id: ruff-check
+ args: [--fix]
+ - id: ruff-format
+
+ - repo: https://github.com/pocc/pre-commit-hooks
+ rev: v1.3.5
+ hooks:
+ - id: clang-format
+ files: \.(c|cc|cpp|cxx|h|hh|hpp|hxx)$
+ - id: clang-tidy
+ files: \.(c|cc|cpp|cxx)$
+ args: ["-checks=-*,clang-analyzer-*,bugprone-*"]
+
+ - repo: https://github.com/pre-commit/pre-commit-hooks
+ rev: v6.0.0
+ hooks:
+ - id: check-added-large-files
+ - id: check-merge-conflict
+ - id: check-yaml
+ exclude: ^conda-recipe/meta\.yaml$
+ - id: check-toml
+ - id: check-case-conflict
+ - id: check-ast
+ - id: end-of-file-fixer
+ - id: trailing-whitespace
+
+ - repo: https://github.com/rstcheck/rstcheck
+ rev: v6.2.0
+ hooks:
+ - id: rstcheck
+ name: rstcheck (docstrings only)
+ files: \.py$
+ args: ["--report-level", "warning"]
+
+ - repo: local
+ hooks:
+ - id: sphinx-docs
+ name: sphinx-build (pre-push, warnings as errors)
+ entry: bash -lc 'make -C docs clean && make -C docs html SPHINXOPTS="-W --keep-going"'
+ language: system
+ pass_filenames: false
+ stages: [pre-push]
diff --git a/.readthedocs.yaml b/.readthedocs.yaml
index 87e396c..08801c7 100644
--- a/.readthedocs.yaml
+++ b/.readthedocs.yaml
@@ -20,5 +20,3 @@ sphinx:
python:
install:
- requirements: docs/requirements.txt
-
-
diff --git a/CITATION.cff b/CITATION.cff
index 3c3d94c..9459d76 100644
--- a/CITATION.cff
+++ b/CITATION.cff
@@ -72,4 +72,4 @@ keywords:
license: MIT
version: 1.0.0
-date-released: '2026-09-14'
\ No newline at end of file
+date-released: '2026-09-14'
diff --git a/README.md b/README.md
index 217c1b8..33548de 100644
--- a/README.md
+++ b/README.md
@@ -1,4 +1,4 @@
-SIST: Stress-Induced Structural Transitions in superhelical DNA
+SIST: Stress-Induced Structural Transitions in superhelical DNA
==============================
| Category | Badges |
@@ -12,7 +12,7 @@ SIST: Stress-Induced Structural Transitions in superhelical DNA
## Purpose of SIST
-The codes in this repository are for analyzing three types of structural transitions in superhelical DNA molecules of specified base sequences and kilobase lengths. These are strand separation, BZ transitions and cruciform extrusion. More types of transitions may be added as their energetics become known. The statistical mechanical methods and algorithms used in these analyses are described in the papers cited below.
+The codes in this repository are for analyzing three types of structural transitions in superhelical DNA molecules of specified base sequences and kilobase lengths. These are strand separation, BZ transitions and cruciform extrusion. More types of transitions may be added as their energetics become known. The statistical mechanical methods and algorithms used in these analyses are described in the papers cited below.
diff --git a/docs/source/_static/examples/one_line.pbr322.toy.fa b/docs/source/_static/examples/one_line.pbr322.toy.fa
index 0ca67ed..8786cd9 100644
--- a/docs/source/_static/examples/one_line.pbr322.toy.fa
+++ b/docs/source/_static/examples/one_line.pbr322.toy.fa
@@ -1,2 +1,2 @@
>one_line.pbr322.toy.fa
-AGTCAGGCACCGTGTATGAAATCTAACAATGCGCTCATCGTCATCCTCGGCACCGTCACCCTGGATGCTGTAGGCATAGGCTTGGTTATGCCGGTACTGCCGGGCCTCTTGCGGGATATCGTCCATTCCGACAGCATCGCCAGTCACTATGGCGTGCTGCTAGCGCTATATGCGTTGATGCAATTTCTATGCGCACCCGTTCTCGGAGCACTGTCCGAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAACGGAGCCACTATCGACTACGCGATCATGGCGACCACACCCGTCCTGTGGATCCTCTACGCCGGACGCATCGTGGCCGGCATCACCGGCGCCACAGGTGCGGTTGCTGGCGCCTATATCGCCGACATCACCGATGGGGAAGATCGGGCTCGCCACTTCGGGCTCATGAGCGCTTGTTTCGGCGTGGGTATGGTGGCAGGCCCCGTGGCCGGGGGACTGTTGGGCGCCATCTCCTTGCATGCACCATTCCTTGCGGCGGCGGTGCTCAACGGCCTCAACCTACTACTGGGCTGCTTCCTAATGCAGGAGTCGCATAAGGGAGAGCGTCGACCGATGCCCTTGAGAGCCTTCAACCCAGTCAGCTCCTTCCGGTGGGCGCGGGGCATGACTATCGTCGCCGCACTTATGACTGTCTTCTTTATCATGCAACTCGTAGGACAGGTGCCGGCAGCGCTCTGGGTCATTTTCGGCGAGGACCGCTTTCGCTGGAGCGCGACGATGATCGGCCTGTCGCTTGCGGTATTCGGAATCTTGCACGCCCTCGCTCAAGCCTTCGTCACTGGTCCCGCCACCAAACGTTTCGGCGAGAAGCAGGCCATTATCGCCGGCATGGCGGCCGACGCGCTGGGCTACGTCTTGCTGGCGTTCGCGACGCGAGGCTGGATGGCCTTCCCCATTATGATTCTTCTCGCTTCCGGCGGCATCGGGATGCCCGCGTTGCAGGCCATGCTGTCCAGGCAGGTCGCGCGCGCGCGCGCGCGCAGATGACAAGGATCGCTCGCGGCTCTTACCAGCCTAACTTCGATCACTGGACCGCTGATCGTCACGGCGATTTATGCCGCCTCGGCGAGCACATGGAACGGGTTGGCATGGATTGTAGGCGCCGCCCTATACCTTGTCTGCCTCCCCGCGTTGCGTCGCGGTGCATGGAGCCGGGCCACCTCGACCTGAATGGAAGCCGGCGGCACCTCGCTAACGGATTCACCACTCCAAGAATTGGAGCCAATCAATTCTTGCGGAGAACTGTGAATGCGCAAACCAACCCTTGGCAGAACATATCCATCGCGTCCGCCATCTCCAGCAGCCGCACGCGGCGCATCTCGGGCAGCGTTGGGTCCTGGCCACGGGTGCGCATGATCGTGCTCCTGTCGTTGAGGACCCGGCTAGGCTGGCGGGGTTGCCTTACTGGTTAGCAGAATGAATCACCGATACGCGAGCGAACGTGAAGCGACTGCTGCTGCAAAACGTCTGCGACCTGAGCAACAACATGAATGGTCTTCGGTTTCCGTGTTTCGTAAAGTCTGGAAACGCGGAAGTCAGCGCCCTGCACCATTATGTTCCGGATCTGCATCGCAGGATGCTGCTGGCTACCCTGTGGAACACCTACATCTGTATTAACGAAGCGCTGGCATTGACCCTGAGTGATTTTTCTCTGGTCCCGCCGCATCCATACCGCCAGTTGTTTACCCTCACAACGTTCCAGTAACCGGGCATGTTCATCATCAGTAACCCGTATCGTGAGCATCCTCTCTCGTTTCATCGGTATCATTACCCCCATGAACAGAAATCCCCCTTACACGGAGGCATCAGTGACCAAACAGGAAAAAACCGCCCTTAACATGGCCCGCTTTATCAGAAGCCAGACATTAACGCTTCTGGAGAAACTCAACGAGCTGGACGCGGATGAACAGGCAGACATCTGTGAATCGCTTCACGACCACGCTGATGAGCTTTACCGCAGCTGCCTCGCGCGTTTCGGTGATGACGGTGAAAACCTCTGACACATGCAGCTCCCGGAGACGGTCACAGCTTGTCTGTAAGCGGATGCCGGGAGCAGACAAGCCCGTCAGGGCGCGTCAGCGGGTGTTGGCGGGTGTCGGGGCGCAGCCATGACCCAGTCACGTAGCGATAGCGGAGTGTATACTGGCTTAACTATGCGGCATCAGAGCAGATTGTACTGAGAGTGCACCATATGCGGTGTGAAATACCGCACAGATGCGTAAGGAGAAAATACCGCATCAGGCGCTCTTCCGCTTCCTCGCTCACTGACTCGCTGCGCTCGGTCGTTCGGCTGCGGCGAGCGGTATCAGCTCACTCAAAGGCGGTAATACGGTTATCCACAGAATCAGGGGATAACGCAGGAAAGAACATGTGAGCAAAAGGCCAGCAAAAGGCCAGGAACCGTAAAAAGGCCGCGTTGCTGGCGTTTTTCCATAGGCTCCGCCCCCCTGACGAGCATCACAAAAATCGACGCTCAAGTCAGAGGTGGCGAAACCCGACAGGACTATAAAGATACCAGGCGTTTCCCCCTGGAAGCTCCCTCGTGCGCTCTCCTGTTCCGACCCTGCCGCTTACCGGATACCTGTCCGCCTTTCTCCCTTCGGGAAGCGTGGCGCTTTCTCATAGCTCACGCTGTAGGTATCTCAGTTCGGTGTAGGTCGTTCGCTCCAAGCTGGGCTGTGTGCACGAACCCCCCGTTCAGCCCGACCGCTGCGCCTTATCCGGTAACTATCGTCTTGAGTCCAACCCGGTAAGACACGACTTATCGCCACTGGCAGCAGCCACTGGTAACAGGATTAGCAGAGCGAGGTATGTAGGCGGTGCTACAGAGTTCTTGAAGTGGTGGCCTAACTACGGCTACACTAGAAGGACAGTATTTGGTATCTGCGCTCTGCTGAAGCCAGTTACCTTCGGAAAAAGAGTTGGTAGCTCTTGATCCGGCAAACAAACCACCGCTGGTAGCGGTGGTTTTTTTGTTTGCAAGCAGCAGATTACGCGCAGAAAAAAAGGATCTCAAGAAGATCCTTTGATCTTTTCTACGGGGTGCTCAGTGAACCAATTGGCCAACCGGAAGGAAAACCTTCCGGTTGGCCAATTGGTTGAACGAAAACTATCCTAGATCCTTTTAAATTAAAAATGAAGTTTTAAATCAATCTAAAGTATATATGAGTAAACTTGGTCTGACAGTTACCAATGCTTAATCAGTGAGGCACCTATCTCAGCGATCTGTCTATTTCGTTCATCCATAGTTGCCTGACTCCCCGTCGTGTAGATAACTACGATACGGGAGGGCTTACCATCTGGCCCCAGTGCTGCAATGATACCGCGAGACCCACGCTCACCGGCTCCAGATTTATCAGCAATAAACCAGCCAGCCGGAAGGGCCGAGCGCAGAAGTGGTCCTGCAACTTTATCCGCCTCCATCCAGTCTATTAATTGTTGCCGGGAAGCTAGAGTAAGTAGTTCGCCAGTTAATAGTTTGCGCAACGTTGTTGCCATTGCTGCAGGCATCGTGGTGTCACGCTCGTCGTTTGGTATGGCTTCATTCAGCTCCGGTTCCCAACGATCAAGGCGAGTTACATGATCCCCCATGTTGTGCAAAAAAGCGGTTAGCTCCTTCGGTCCTCCGATCGTTGTCAGAAGTAAGTTGGCCGCAGTGTTATCACTCATGGTTATGGCAGCACTGCATAATTCTCTTACTGTCATGCCATCCGTAAGATGCTTTTCTGTGACTGGTGAGTACTCAACCAAGTCATTCTGAGAATAGTGTATGCGGCGACCGAGTTGCTCTTGCCCGGCGTCAACACGGGATAATACCGCGCCACATAGCAGAACTTTAAAAGTGCTCATCATTGGAAAACGTTCTTCGGGGCGAAAACTCTCAAGGATCTTACCGCTGTTGAGATCCAGTTCGATGTAACCCACTCGTGCACCCAACTGATCTTCAGCATCTTTTACTTTCACCAGCGTTTCTGGGTGAGCAAAAACAGGAAGGCAAAATGCCGCAAAAAAGGGAATAAGGGCGACACGGAAATGTTGAATACTCATACTCTTCCTTTTTCAATATTATTGAAGCATTTATCAGGGTTATTGTCTCATGAGCGGATACATATTTGAATGTATTTAGAAAAATAAACAAATAGGGGTTCCGCGCACATTTCCCCGAAAAGTGCCACCTGACGTCTAAGAAACCATTATTATCATGACATTAACCTATAAAAATAGGCGTATCACGAGGCCCTTTCGTCTTCAAGAA
\ No newline at end of file
+AGTCAGGCACCGTGTATGAAATCTAACAATGCGCTCATCGTCATCCTCGGCACCGTCACCCTGGATGCTGTAGGCATAGGCTTGGTTATGCCGGTACTGCCGGGCCTCTTGCGGGATATCGTCCATTCCGACAGCATCGCCAGTCACTATGGCGTGCTGCTAGCGCTATATGCGTTGATGCAATTTCTATGCGCACCCGTTCTCGGAGCACTGTCCGAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAACGGAGCCACTATCGACTACGCGATCATGGCGACCACACCCGTCCTGTGGATCCTCTACGCCGGACGCATCGTGGCCGGCATCACCGGCGCCACAGGTGCGGTTGCTGGCGCCTATATCGCCGACATCACCGATGGGGAAGATCGGGCTCGCCACTTCGGGCTCATGAGCGCTTGTTTCGGCGTGGGTATGGTGGCAGGCCCCGTGGCCGGGGGACTGTTGGGCGCCATCTCCTTGCATGCACCATTCCTTGCGGCGGCGGTGCTCAACGGCCTCAACCTACTACTGGGCTGCTTCCTAATGCAGGAGTCGCATAAGGGAGAGCGTCGACCGATGCCCTTGAGAGCCTTCAACCCAGTCAGCTCCTTCCGGTGGGCGCGGGGCATGACTATCGTCGCCGCACTTATGACTGTCTTCTTTATCATGCAACTCGTAGGACAGGTGCCGGCAGCGCTCTGGGTCATTTTCGGCGAGGACCGCTTTCGCTGGAGCGCGACGATGATCGGCCTGTCGCTTGCGGTATTCGGAATCTTGCACGCCCTCGCTCAAGCCTTCGTCACTGGTCCCGCCACCAAACGTTTCGGCGAGAAGCAGGCCATTATCGCCGGCATGGCGGCCGACGCGCTGGGCTACGTCTTGCTGGCGTTCGCGACGCGAGGCTGGATGGCCTTCCCCATTATGATTCTTCTCGCTTCCGGCGGCATCGGGATGCCCGCGTTGCAGGCCATGCTGTCCAGGCAGGTCGCGCGCGCGCGCGCGCGCAGATGACAAGGATCGCTCGCGGCTCTTACCAGCCTAACTTCGATCACTGGACCGCTGATCGTCACGGCGATTTATGCCGCCTCGGCGAGCACATGGAACGGGTTGGCATGGATTGTAGGCGCCGCCCTATACCTTGTCTGCCTCCCCGCGTTGCGTCGCGGTGCATGGAGCCGGGCCACCTCGACCTGAATGGAAGCCGGCGGCACCTCGCTAACGGATTCACCACTCCAAGAATTGGAGCCAATCAATTCTTGCGGAGAACTGTGAATGCGCAAACCAACCCTTGGCAGAACATATCCATCGCGTCCGCCATCTCCAGCAGCCGCACGCGGCGCATCTCGGGCAGCGTTGGGTCCTGGCCACGGGTGCGCATGATCGTGCTCCTGTCGTTGAGGACCCGGCTAGGCTGGCGGGGTTGCCTTACTGGTTAGCAGAATGAATCACCGATACGCGAGCGAACGTGAAGCGACTGCTGCTGCAAAACGTCTGCGACCTGAGCAACAACATGAATGGTCTTCGGTTTCCGTGTTTCGTAAAGTCTGGAAACGCGGAAGTCAGCGCCCTGCACCATTATGTTCCGGATCTGCATCGCAGGATGCTGCTGGCTACCCTGTGGAACACCTACATCTGTATTAACGAAGCGCTGGCATTGACCCTGAGTGATTTTTCTCTGGTCCCGCCGCATCCATACCGCCAGTTGTTTACCCTCACAACGTTCCAGTAACCGGGCATGTTCATCATCAGTAACCCGTATCGTGAGCATCCTCTCTCGTTTCATCGGTATCATTACCCCCATGAACAGAAATCCCCCTTACACGGAGGCATCAGTGACCAAACAGGAAAAAACCGCCCTTAACATGGCCCGCTTTATCAGAAGCCAGACATTAACGCTTCTGGAGAAACTCAACGAGCTGGACGCGGATGAACAGGCAGACATCTGTGAATCGCTTCACGACCACGCTGATGAGCTTTACCGCAGCTGCCTCGCGCGTTTCGGTGATGACGGTGAAAACCTCTGACACATGCAGCTCCCGGAGACGGTCACAGCTTGTCTGTAAGCGGATGCCGGGAGCAGACAAGCCCGTCAGGGCGCGTCAGCGGGTGTTGGCGGGTGTCGGGGCGCAGCCATGACCCAGTCACGTAGCGATAGCGGAGTGTATACTGGCTTAACTATGCGGCATCAGAGCAGATTGTACTGAGAGTGCACCATATGCGGTGTGAAATACCGCACAGATGCGTAAGGAGAAAATACCGCATCAGGCGCTCTTCCGCTTCCTCGCTCACTGACTCGCTGCGCTCGGTCGTTCGGCTGCGGCGAGCGGTATCAGCTCACTCAAAGGCGGTAATACGGTTATCCACAGAATCAGGGGATAACGCAGGAAAGAACATGTGAGCAAAAGGCCAGCAAAAGGCCAGGAACCGTAAAAAGGCCGCGTTGCTGGCGTTTTTCCATAGGCTCCGCCCCCCTGACGAGCATCACAAAAATCGACGCTCAAGTCAGAGGTGGCGAAACCCGACAGGACTATAAAGATACCAGGCGTTTCCCCCTGGAAGCTCCCTCGTGCGCTCTCCTGTTCCGACCCTGCCGCTTACCGGATACCTGTCCGCCTTTCTCCCTTCGGGAAGCGTGGCGCTTTCTCATAGCTCACGCTGTAGGTATCTCAGTTCGGTGTAGGTCGTTCGCTCCAAGCTGGGCTGTGTGCACGAACCCCCCGTTCAGCCCGACCGCTGCGCCTTATCCGGTAACTATCGTCTTGAGTCCAACCCGGTAAGACACGACTTATCGCCACTGGCAGCAGCCACTGGTAACAGGATTAGCAGAGCGAGGTATGTAGGCGGTGCTACAGAGTTCTTGAAGTGGTGGCCTAACTACGGCTACACTAGAAGGACAGTATTTGGTATCTGCGCTCTGCTGAAGCCAGTTACCTTCGGAAAAAGAGTTGGTAGCTCTTGATCCGGCAAACAAACCACCGCTGGTAGCGGTGGTTTTTTTGTTTGCAAGCAGCAGATTACGCGCAGAAAAAAAGGATCTCAAGAAGATCCTTTGATCTTTTCTACGGGGTGCTCAGTGAACCAATTGGCCAACCGGAAGGAAAACCTTCCGGTTGGCCAATTGGTTGAACGAAAACTATCCTAGATCCTTTTAAATTAAAAATGAAGTTTTAAATCAATCTAAAGTATATATGAGTAAACTTGGTCTGACAGTTACCAATGCTTAATCAGTGAGGCACCTATCTCAGCGATCTGTCTATTTCGTTCATCCATAGTTGCCTGACTCCCCGTCGTGTAGATAACTACGATACGGGAGGGCTTACCATCTGGCCCCAGTGCTGCAATGATACCGCGAGACCCACGCTCACCGGCTCCAGATTTATCAGCAATAAACCAGCCAGCCGGAAGGGCCGAGCGCAGAAGTGGTCCTGCAACTTTATCCGCCTCCATCCAGTCTATTAATTGTTGCCGGGAAGCTAGAGTAAGTAGTTCGCCAGTTAATAGTTTGCGCAACGTTGTTGCCATTGCTGCAGGCATCGTGGTGTCACGCTCGTCGTTTGGTATGGCTTCATTCAGCTCCGGTTCCCAACGATCAAGGCGAGTTACATGATCCCCCATGTTGTGCAAAAAAGCGGTTAGCTCCTTCGGTCCTCCGATCGTTGTCAGAAGTAAGTTGGCCGCAGTGTTATCACTCATGGTTATGGCAGCACTGCATAATTCTCTTACTGTCATGCCATCCGTAAGATGCTTTTCTGTGACTGGTGAGTACTCAACCAAGTCATTCTGAGAATAGTGTATGCGGCGACCGAGTTGCTCTTGCCCGGCGTCAACACGGGATAATACCGCGCCACATAGCAGAACTTTAAAAGTGCTCATCATTGGAAAACGTTCTTCGGGGCGAAAACTCTCAAGGATCTTACCGCTGTTGAGATCCAGTTCGATGTAACCCACTCGTGCACCCAACTGATCTTCAGCATCTTTTACTTTCACCAGCGTTTCTGGGTGAGCAAAAACAGGAAGGCAAAATGCCGCAAAAAAGGGAATAAGGGCGACACGGAAATGTTGAATACTCATACTCTTCCTTTTTCAATATTATTGAAGCATTTATCAGGGTTATTGTCTCATGAGCGGATACATATTTGAATGTATTTAGAAAAATAAACAAATAGGGGTTCCGCGCACATTTCCCCGAAAAGTGCCACCTGACGTCTAAGAAACCATTATTATCATGACATTAACCTATAAAAATAGGCGTATCACGAGGCCCTTTCGTCTTCAAGAA
diff --git a/docs/source/_static/examples/one_line.pbr322.toy.fa.2.10.10.80.10.20.10000.100.1.html b/docs/source/_static/examples/one_line.pbr322.toy.fa.2.10.10.80.10.20.10000.100.1.html
index 922ccd8..1853342 100644
--- a/docs/source/_static/examples/one_line.pbr322.toy.fa.2.10.10.80.10.20.10000.100.1.html
+++ b/docs/source/_static/examples/one_line.pbr322.toy.fa.2.10.10.80.10.20.10000.100.1.html
@@ -1,9 +1,9 @@
one_line.pbr322.toy.fa.2.10.10.80.10.20.10000.100.1.html
Inverted Repeats Finder Program written by:
Gary Benson
Department of Biomathematical Sciences
Mount Sinai School of Medicine
Version 3.05
Sequence: one_line.pbr322.toy.fa
-Parameters: 2 10 10 80 10 20 10000 100
+Parameters: 2 10 10 80 10 20 10000 100
Length: 4291
-Tables: 1
+Tables: 1
This is table 1 of 1 ( 2 repeats found )
@@ -16,7 +16,7 @@
-Tables: 1
+Tables: 1
The End!
diff --git a/docs/source/_static/examples/one_line.pbr322.toy.fa.2.10.10.80.10.20.10000.100.1.txt.html b/docs/source/_static/examples/one_line.pbr322.toy.fa.2.10.10.80.10.20.10000.100.1.txt.html
index bc3b275..319ac65 100644
--- a/docs/source/_static/examples/one_line.pbr322.toy.fa.2.10.10.80.10.20.10000.100.1.txt.html
+++ b/docs/source/_static/examples/one_line.pbr322.toy.fa.2.10.10.80.10.20.10000.100.1.txt.html
@@ -8,7 +8,7 @@
Version 3.05
Sequence: one_line.pbr322.toy.fa
-Parameters: 2 10 10 80 10 20 10000 100
+Parameters: 2 10 10 80 10 20 10000 100
Length: 4291
ACGTcount: A:0.23, C:0.28, G:0.26, T:0.23, N:0.00
@@ -25,7 +25,7 @@
2972 >> (LF) TCCGGCAAAC
3016 << (RF) CGTTTGTTTT
-
+
2982 >> AAACCACCGCT >> 2992
3006 << *********** << 2996
@@ -60,7 +60,7 @@
3079 >> (LF) GTGCTCAGTG
3146 << (RF) CAAAAGCAAG
-
+
3089 >> AACCAATTGGCCAACCGGAAGG >> 3110
3136 << ********************** << 3115
diff --git a/docs/source/_static/examples/pbr322.toy.fa b/docs/source/_static/examples/pbr322.toy.fa
index 23324fb..ab43845 100644
--- a/docs/source/_static/examples/pbr322.toy.fa
+++ b/docs/source/_static/examples/pbr322.toy.fa
@@ -61,4 +61,3 @@ ACACGGAAATGTTGAATACTCATACTCTTCCTTTTTCAATATTATTGAAGCATTTATCAGGGTTATTGTC
TCATGAGCGGATACATATTTGAATGTATTTAGAAAAATAAACAAATAGGGGTTCCGCGCACATTTCCCCG
AAAAGTGCCACCTGACGTCTAAGAAACCATTATTATCATGACATTAACCTATAAAAATAGGCGTATCACG
AGGCCCTTTCGTCTTCAAGAA
-
diff --git a/docs/source/_static/logos/SIST-logo-black-text.svg b/docs/source/_static/logos/SIST-logo-black-text.svg
index 9974caa..0052e2b 100644
--- a/docs/source/_static/logos/SIST-logo-black-text.svg
+++ b/docs/source/_static/logos/SIST-logo-black-text.svg
@@ -1 +1 @@
-
\ No newline at end of file
+
diff --git a/docs/source/_static/logos/SIST-logo-white-text.svg b/docs/source/_static/logos/SIST-logo-white-text.svg
index 66ae614..865f9b2 100644
--- a/docs/source/_static/logos/SIST-logo-white-text.svg
+++ b/docs/source/_static/logos/SIST-logo-white-text.svg
@@ -1 +1 @@
-
\ No newline at end of file
+
diff --git a/docs/source/citation.rst b/docs/source/citation.rst
index d3d390e..18ef49f 100644
--- a/docs/source/citation.rst
+++ b/docs/source/citation.rst
@@ -54,4 +54,4 @@ The original SIST documentation listed the following contacts:
* Dina Zhabinskaya - dzhabinskaya@ucdavis.edu
* Craig Benham - cjbenham@ucdavis.edu
-* Sally Madden - sallymadden@gmail.com
\ No newline at end of file
+* Sally Madden - sallymadden@gmail.com
diff --git a/docs/source/conf.py b/docs/source/conf.py
index 0d04b6b..8c46a5a 100644
--- a/docs/source/conf.py
+++ b/docs/source/conf.py
@@ -6,8 +6,8 @@
# -- Project information -----------------------------------------------------
# https://www.sphinx-doc.org/en/master/usage/configuration.html#project-information
-project = 'SIST'
-copyright = '2026, CCPBioSim'
+project = "SIST"
+copyright = "2026, CCPBioSim"
# -- General configuration ---------------------------------------------------
# https://www.sphinx-doc.org/en/master/usage/configuration.html#general-configuration
@@ -29,18 +29,17 @@
napoleon_use_param = False
napoleon_use_ivar = True
-templates_path = ['_templates']
+templates_path = ["_templates"]
exclude_patterns = []
-
# -- Options for HTML output -------------------------------------------------
# https://www.sphinx-doc.org/en/master/usage/configuration.html#options-for-html-output
-html_theme = 'furo'
+html_theme = "furo"
html_theme_options = {
"dark_logo": "logos/SIST-logo-white-text.svg",
"light_logo": "logos/SIST-logo-black-text.svg",
}
-html_static_path = ['_static']
+html_static_path = ["_static"]
diff --git a/docs/source/development.rst b/docs/source/development.rst
index 88a495a..16ff24f 100644
--- a/docs/source/development.rst
+++ b/docs/source/development.rst
@@ -48,6 +48,36 @@ During a normal source test run, the test fixtures create a temporary copy of
the repository, build the C++ executables, and run the supported calculations
from that working copy.
+Pre-commit hooks
+----------------
+
+SIST uses **pre-commit hooks** to maintain code quality and consistent style
+across the Python and C++ code.
+
+Install the pre-commit dependencies and enable the hooks:
+
+.. code-block:: bash
+
+ python -m pip install -e '.[pre-commit]'
+ pre-commit install
+
+Our tooling stack:
+
+* **Python linting and formatting** via ``ruff``
+* **C++ formatting and static analysis** via ``clang-format`` and ``clang-tidy``
+* **Basic repository checks** via ``pre-commit-hooks``
+* **Docstring RST validation** via ``rstcheck``
+
+Run the checks manually against the whole repository:
+
+.. code-block:: bash
+
+ pre-commit run --all-files
+
+.. note::
+
+ Pull requests must pass all pre-commit checks before being merged.
+
Scientific regression baselines
--------------------------------
diff --git a/pyproject.toml b/pyproject.toml
index 579aa2e..020d352 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -34,7 +34,7 @@ Repository = "https://github.com/CCPBioSim/SIST"
[project.scripts]
sist = "sist.cli:main"
# Deprecated aliases for scripted workflows that still invoke the old Perl
-# script names directly.
+# script names directly.
"master.pl" = "sist.cli:main"
"IR_finder.pl" = "sist.ir_finder_cli:main"
@@ -45,6 +45,15 @@ testing = [
"mypy>=1.14,<2.0",
]
+pre-commit = [
+ "pre-commit>=4.5,<5.0",
+ "ruff>=0.16,<0.17",
+ "pylint>=4.0,<5.0",
+ "rstcheck>=6.2,<7.0",
+ "clang-format>=23.0,<24.0",
+ "clang-tidy>=22.0,<23.0",
+]
+
[tool.setuptools.packages.find]
where = ["src"]
include = ["sist*"]
diff --git a/src/trans_compete/G_x.cpp b/src/trans_compete/G_x.cpp
index 8d2a047..824fbac 100644
--- a/src/trans_compete/G_x.cpp
+++ b/src/trans_compete/G_x.cpp
@@ -22,54 +22,47 @@
// Construction/Destruction
//////////////////////////////////////////////////////////////////////
-G_x::G_x()
-{
+G_x::G_x() {
- sum_xG = 0.0;
- sum_xB = 0.0;
- ave_Gx = 0.0;
- px = 0.0;
+ sum_xG = 0.0;
+ sum_xB = 0.0;
+ ave_Gx = 0.0;
+ px = 0.0;
}
-void G_x::reset()
-{
- sum_xG = 0.0;
- sum_xB = 0.0;
- ave_Gx = 0.0;
- px = 0.0;
+void G_x::reset() {
+ sum_xG = 0.0;
+ sum_xB = 0.0;
+ ave_Gx = 0.0;
+ px = 0.0;
}
-void G_x::add(double gs, double exponent, double rt, double lastG, double lastB)
-{
- if(gs == -1.0 && rt == -1.0){
- sum_xG = sum_xG + lastG;
- sum_xB = sum_xB + lastB;
- }
- else{
- sum_xG = sum_xG + gs*exponent + lastG;
- sum_xB = sum_xB + exponent + lastB;
- }
+void G_x::add(double gs, double exponent, double rt, double lastG,
+ double lastB) {
+ if (gs == -1.0 && rt == -1.0) {
+ sum_xG = sum_xG + lastG;
+ sum_xB = sum_xB + lastB;
+ } else {
+ sum_xG = sum_xG + gs * exponent + lastG;
+ sum_xB = sum_xB + exponent + lastB;
+ }
}
-//intput: average energy and partition function
-void G_x::calc(double ave_Gs, double ZsumB)
-{
- if (ZsumB == 0.0)
- px = 1.0;
- else
- // V: This is the calculation for the probability
- px = sum_xB / ZsumB;
-
- if(sum_xB == 0.0)
- ave_Gx = INFINITE_Gx;
- else
- // V: This is the calculation of the energy required to transition
- ave_Gx = sum_xG / sum_xB - ave_Gs;
-}
-
-G_x::~G_x()
-{
+// intput: average energy and partition function
+void G_x::calc(double ave_Gs, double ZsumB) {
+ if (ZsumB == 0.0)
+ px = 1.0;
+ else
+ // V: This is the calculation for the probability
+ px = sum_xB / ZsumB;
+ if (sum_xB == 0.0)
+ ave_Gx = INFINITE_Gx;
+ else
+ // V: This is the calculation of the energy required to transition
+ ave_Gx = sum_xG / sum_xB - ave_Gs;
}
+G_x::~G_x() {}
+
#endif
diff --git a/src/trans_compete/G_x.h b/src/trans_compete/G_x.h
index e0cc1e8..3c83215 100644
--- a/src/trans_compete/G_x.h
+++ b/src/trans_compete/G_x.h
@@ -18,26 +18,25 @@
#include
const double INFINITE_Gx = -10000.0;
-class G_x
-{
+class G_x {
private:
- double sum_xG; // total energy with states of x opening
- double sum_xB; // sum of Boltzman's factors with states of x open
- double ave_Gx; // average free energy of opening base x
- double px; // probability of opening base x
+ double sum_xG; // total energy with states of x opening
+ double sum_xB; // sum of Boltzman's factors with states of x open
+ double ave_Gx; // average free energy of opening base x
+ double px; // probability of opening base x
public:
- G_x();
- void reset();
- void add(double gs, double exponent, double rt, double lastG = 0.0, double lastB = 0.0); // add one state
- void calc(double ave_Gs, double zsumb); // computing p(x) and G(x)
- double get_sum_xG(){return sum_xG;}
- double get_ave_Gx(){return ave_Gx;}
- void set_ave_Gx(double x){ave_Gx = x;}
- double get_sum_xB(){return sum_xB;}
- double get_px(){return px;}
- virtual ~G_x();
-
+ G_x();
+ void reset();
+ void add(double gs, double exponent, double rt, double lastG = 0.0,
+ double lastB = 0.0); // add one state
+ void calc(double ave_Gs, double zsumb); // computing p(x) and G(x)
+ double get_sum_xG() { return sum_xG; }
+ double get_ave_Gx() { return ave_Gx; }
+ void set_ave_Gx(double x) { ave_Gx = x; }
+ double get_sum_xB() { return sum_xB; }
+ double get_px() { return px; }
+ virtual ~G_x();
};
#endif // !defined(AFX_G_X_H__024B41EB_E8B4_4F19_BA76_1276D508094E__INCLUDED_)
diff --git a/src/trans_compete/Makefile b/src/trans_compete/Makefile
index 93b2108..80609b4 100644
--- a/src/trans_compete/Makefile
+++ b/src/trans_compete/Makefile
@@ -36,7 +36,7 @@ depend: $(OBJECTS:.o=.cpp)
$(CXX) -MM $^ > $@
test:
- ./qsidd example.fasta
+ ./qsidd example.fasta
tar:
rm -rf /tmp/$(APP)
@@ -70,4 +70,3 @@ coverage:
################
include depend
-
diff --git a/src/trans_compete/SIDD_1R.cpp b/src/trans_compete/SIDD_1R.cpp
index f99cfcd..33a23f1 100644
--- a/src/trans_compete/SIDD_1R.cpp
+++ b/src/trans_compete/SIDD_1R.cpp
@@ -22,136 +22,154 @@
// Construction/Destruction
//////////////////////////////////////////////////////////////////////
-SIDD_1R::SIDD_1R():SIDD_Base() {
- for (int t=0; t<=Ns; t++)
- min_RE[t] = min_E; //minimun energy from the continuous part
- flag_minE1 = false;
- StatesInOneRun = 0;
-}
-SIDD_1R::~SIDD_1R()
-{
+SIDD_1R::SIDD_1R() : SIDD_Base() {
+ for (int t = 0; t <= Ns; t++)
+ min_RE[t] = min_E; // minimun energy from the continuous part
+ flag_minE1 = false;
+ StatesInOneRun = 0;
}
+SIDD_1R::~SIDD_1R() {}
+
+void SIDD_1R::gen_OpenBaseEnergy() {
+ double e;
+ // t=0 for melting, t=1 for Z-DNA, and t=2 for cruciforms
+ for (int t = 0; t <= Ns; t++) {
+ for (int i = MinWindowSize; i <= MaxWindowSize;
+ i++) { // window size: 1 - MaxWindowSizez
+ for (int j = 0; j < length_seq; j++) { // start position to length_seq - 1
+ // V: a - is the nucleation energy to initiate run of transition t
+ // V: calc_OPenBasesEnergy. The function is in SIDD_base. It calculates
+ // the energy associated
+ // with transition t given the run at position startp:j and with size
+ // n:i
+ e = a[t] + calc_OPenBasesEnergy(j, i, t);
+ // V: It creates object of type stat_1R. It seems it just stores the
+ // values of position j,
+ // energy e and RT
+ stat_1R s1r(j, e,
+ RT); // input j=starting position, e=energy, RT=constant
+ // V: It then appends this object the Lst_OBE array/list.
+ // This Lst_OBE seems to be a list made for each run/window of size i
+ // with transition t.
+ Lst_OBE[i][t].push_back(s1r);
+ if (e < min_RE[t])
+ min_RE[t] = e;
+ }
+ }
+ }
-void SIDD_1R::gen_OpenBaseEnergy()
-{
- double e;
- //t=0 for melting, t=1 for Z-DNA, and t=2 for cruciforms
- for (int t=0; t<=Ns; t++) {
- for(int i = MinWindowSize; i <= MaxWindowSize; i++){ // window size: 1 - MaxWindowSizez
- for(int j = 0; j < length_seq; j++){ // start position to length_seq - 1
- // V: a - is the nucleation energy to initiate run of transition t
- // V: calc_OPenBasesEnergy. The function is in SIDD_base. It calculates the energy associated
- // with transition t given the run at position startp:j and with size n:i
- e = a[t] + calc_OPenBasesEnergy(j, i, t);
- // V: It creates object of type stat_1R. It seems it just stores the values of position j,
- // energy e and RT
- stat_1R s1r(j, e, RT); //input j=starting position, e=energy, RT=constant
- // V: It then appends this object the Lst_OBE array/list.
- // This Lst_OBE seems to be a list made for each run/window of size i with transition t.
- Lst_OBE[i][t].push_back(s1r);
- if(e < min_RE[t]) min_RE[t] = e;
- }
- }
- }
-
- // sorting each list of OBE
- // V: I see, for the list of windows/runs of size k of transition t=0-melting, t=1-Z, t=2 cruciform, it sorts
- // from low to high. I think this is useful for filtering high energy states...
- for (int t=0; t<=Ns; t++) {
- for(int k = MinWindowSize; k <= MaxWindowSize; k++){
- Lst_OBE[k][t].sort();
- }
- }
+ // sorting each list of OBE
+ // V: I see, for the list of windows/runs of size k of transition t=0-melting,
+ // t=1-Z, t=2 cruciform, it sorts
+ // from low to high. I think this is useful for filtering high energy
+ // states...
+ for (int t = 0; t <= Ns; t++) {
+ for (int k = MinWindowSize; k <= MaxWindowSize; k++) {
+ Lst_OBE[k][t].sort();
+ }
+ }
}
+bool SIDD_1R::Search_Low1RE() {
+ long count = 0;
+ flag_minE1 = false;
+ sum_1RG = sum_1RB = 0.0;
+ for (int t = 0; t <= Ns; t++) {
+ exp_one[t] = 0;
+ runs[t] = 0;
+ Prob[t] = 0;
+ }
-bool SIDD_1R::Search_Low1RE()
-{
- long count = 0;
- flag_minE1 = false;
- sum_1RG = sum_1RB = 0.0;
- for (int t=0; t<=Ns; t++) {
- exp_one[t] = 0;
- runs[t]=0;
- Prob[t]=0;
+ LstState::iterator j;
+ // collecting one run states withing the energy threshold
+ // V: t is the index indicating transition, t=0 melting, t=1 Z-DNA and t=2
+ // cruciform
+ for (int t = 0; t <= Ns; t++) {
+ for (int i = MinWindowSize; i <= MaxWindowSize;
+ i++) { // window size: 1 - MaxWindowSize
+ // V: This iterates over the list which is composed of states with run of
+ // size i with transition t.
+ for (j = Lst_OBE[i][t].begin(); j != Lst_OBE[i][t].end(); ++j) {
+ // V: This line retrieves the energy from the iterator j.
+ double e = j->get_energy();
+ if (e >= 10000)
+ break;
+ int pos = j->get_pos1();
+ // V: Calculates the superhelical energy? And adds it. I think it is the
+ // residual energy.
+ // I think it makes sense because each transition absorbs certain
+ // amount of superhelicity / free energy
+ e += calc_Gres(i * delta_fnc(t, 0), i * delta_fnc(t, 1),
+ i * delta_fnc(t, 2), t * delta_fnc(t, 1));
+ // V: This pretty much does what it says. it doesn't count high energy
+ // states
+ if (e >= max_E)
+ break; // outside of threshold
+ count++; // counting the number of one run states
+ // V: If minimum energy found, register it.
+ if (e < min_E) {
+ min_E = e; // update min_E
+ // V: theta = threshold. You always look for states that are within
+ // the minimum energy state
+ // plus the threshold.
+ max_E = e + theta;
+ min_WS = i; // the window size corresponding to the minimum energy
+ flag_minE1 = true;
+ }
+ // V: This is where it stores some values
+ // V: So, basically what this does is that the exponent of the energy e
+ // retrieved from
+ // Lst_OBE (see above) is calculated.
+ // Then, these exponents are added to the quantitiy exp_one[t]
+ // according transition t. The number of runs[t] per transition t
+ // starts collecting this exponents as well. Then the matrix
+ // update_promatrix is updated - I still don't know what promatrix
+ // is, but I think it is a vector list with coordinates
+ // promatrix[window_size][structural][position] and contains
+ // startp:position, n:window_size, struc:structural_transition,
+ // x:free_energy, bzfactor:Boltzman factor I believe each entry of
+ // promatrix has the form of G_x, but I'm not that sure sum_1RG sums
+ // all the states within one run. sum_1RG sums all Boltzman factors
+ // within one run
+ if (write_profile) {
+ // double exponent= exp(-e/RT);
+ double exponent =
+ exp(-e / RT + scaling_factor); // V: I added this scaling factor
+ // double exponent_test = exp(-e/RT + 200.0);
+ // V: exp_one[] is a vector of dimension 3, that sums all the
+ // exponents of runs 1.
+ // it is used for calculating the probability of transition (see
+ // below). that probability sums all these calculated exponents.
+ exp_one[t] += exponent;
+ // V: It doesn't make a difference if runs is inside the loop?
+ // But I guess it counts the number of runs.
+ runs[t] = exp_one[t];
+ update_promatrix(pos, i, t, e, exponent);
+ sum_1RB += exponent;
+ // sum_1RB += exponent_test;
+ sum_1RG = sum_1RG + e * exponent;
+ }
+ }
}
+ }
- LstState::iterator j;
- //collecting one run states withing the energy threshold
- // V: t is the index indicating transition, t=0 melting, t=1 Z-DNA and t=2 cruciform
- for (int t=0; t<=Ns; t++) {
- for(int i = MinWindowSize; i <= MaxWindowSize; i++){ // window size: 1 - MaxWindowSize
- // V: This iterates over the list which is composed of states with run of size i with transition t.
- for(j = Lst_OBE[i][t].begin(); j != Lst_OBE[i][t].end(); ++j){
- // V: This line retrieves the energy from the iterator j.
- double e = j->get_energy();
- if (e >= 10000)
- break;
- int pos =j->get_pos1();
- // V: Calculates the superhelical energy? And adds it. I think it is the residual energy.
- // I think it makes sense because each transition absorbs certain amount of superhelicity /
- // free energy
- e += calc_Gres(i*delta_fnc(t,0),i*delta_fnc(t,1),i*delta_fnc(t,2),t*delta_fnc(t,1));
- // V: This pretty much does what it says. it doesn't count high energy states
- if(e >= max_E) break; //outside of threshold
- count++; //counting the number of one run states
- // V: If minimum energy found, register it.
- if(e < min_E){
- min_E = e; //update min_E
- // V: theta = threshold. You always look for states that are within the minimum energy state
- // plus the threshold.
- max_E = e + theta;
- min_WS = i; // the window size corresponding to the minimum energy
- flag_minE1 = true;
- }
- // V: This is where it stores some values
- // V: So, basically what this does is that the exponent of the energy e retrieved from
- // Lst_OBE (see above) is calculated.
- // Then, these exponents are added to the quantitiy exp_one[t] according transition t.
- // The number of runs[t] per transition t starts collecting this exponents as well.
- // Then the matrix update_promatrix is updated - I still don't know what promatrix is, but I think
- // it is a vector list with coordinates promatrix[window_size][structural][position] and contains
- // startp:position, n:window_size, struc:structural_transition, x:free_energy,
- // bzfactor:Boltzman factor
- // I believe each entry of promatrix has the form of G_x, but I'm not that sure
- // sum_1RG sums all the states within one run.
- // sum_1RG sums all Boltzman factors within one run
- if(write_profile){
- // double exponent= exp(-e/RT);
- double exponent= exp(-e/RT + scaling_factor); // V: I added this scaling factor
- // double exponent_test = exp(-e/RT + 200.0);
- // V: exp_one[] is a vector of dimension 3, that sums all the exponents of runs 1.
- // it is used for calculating the probability of transition (see below).
- // that probability sums all these calculated exponents.
- exp_one[t] += exponent;
- // V: It doesn't make a difference if runs is inside the loop?
- // But I guess it counts the number of runs.
- runs[t]=exp_one[t];
- update_promatrix(pos, i, t, e, exponent);
- sum_1RB += exponent;
- // sum_1RB += exponent_test;
- sum_1RG = sum_1RG + e*exponent;
- }
- }
- }
- }
-
- StatesInOneRun = count;
- // adding the zero open bases term to the partition function
- double closed = exp(-alpha*alpha*K/2/RT + scaling_factor);
- // V: Note that the zero runs state is added to the partition function, but not to the probabilities.
- ZsumB = sum_1RB + closed;
- // V: Note that ZsumG has the form like sum_1RG (up there). Which is sum_1RG + e*exponent, where e = alpha*alpha*K/2
- ZsumG = sum_1RG+alpha*alpha*K/2*closed;
+ StatesInOneRun = count;
+ // adding the zero open bases term to the partition function
+ double closed = exp(-alpha * alpha * K / 2 / RT + scaling_factor);
+ // V: Note that the zero runs state is added to the partition function, but
+ // not to the probabilities.
+ ZsumB = sum_1RB + closed;
+ // V: Note that ZsumG has the form like sum_1RG (up there). Which is sum_1RG +
+ // e*exponent, where e = alpha*alpha*K/2
+ ZsumG = sum_1RG + alpha * alpha * K / 2 * closed;
- if(write_profile && results) {
- cout << "Number of one-run states = " << count << endl;
- for (int t=0; t<=Ns; t++) {
- Prob[t] += exp_one[t];
- }
+ if (write_profile && results) {
+ cout << "Number of one-run states = " << count << endl;
+ for (int t = 0; t <= Ns; t++) {
+ Prob[t] += exp_one[t];
}
- return flag_minE1;
+ }
+ return flag_minE1;
}
#endif
-
diff --git a/src/trans_compete/SIDD_1R.h b/src/trans_compete/SIDD_1R.h
index a1be4c8..74dcfea 100644
--- a/src/trans_compete/SIDD_1R.h
+++ b/src/trans_compete/SIDD_1R.h
@@ -20,47 +20,47 @@
using namespace std;
-typedef list LstState;
+typedef list LstState;
-class SIDD_1R : public SIDD_Base
-{
+class SIDD_1R : public SIDD_Base {
protected:
- LstState Lst_OBE[MaxInitialWindowSize+1][3]; // lists of open base energy
- bool flag_minE1; // flagging if new min_E found in one run
- double min_RE[3]; // store minimum run energy (a + NI)
- long StatesInOneRun;
+ LstState Lst_OBE[MaxInitialWindowSize + 1][3]; // lists of open base energy
+ bool flag_minE1; // flagging if new min_E found in one run
+ double min_RE[3]; // store minimum run energy (a + NI)
+ long StatesInOneRun;
public:
- SIDD_1R();
- virtual ~SIDD_1R();
- void gen_OpenBaseEnergy();
- bool Search_Low1RE();
-
- void Update_Low1RE();
- void Show_Low1RE();
- double runs[3];
- double Prob[3];
- double exp_one[3];
- double exp_two[3][3];
- double exp_three[3][3][3];
- double exp_four[3][3][3][3];
+ SIDD_1R();
+ virtual ~SIDD_1R();
+ void gen_OpenBaseEnergy();
+ bool Search_Low1RE();
- double exp_type_B;
+ void Update_Low1RE();
+ void Show_Low1RE();
+ double runs[3];
+ double Prob[3];
+ double exp_one[3];
+ double exp_two[3][3];
+ double exp_three[3][3][3];
+ double exp_four[3][3][3][3];
- // V: This Prob_compete parameter was added for the competition branch
- // t=0 for melting, t=1 for Z-DNA, and t=2 for cruciforms
- double Prob_compete[3][3]; // V: This should be the probability of intersection?
- // When [0][1], means you have runs of t1=0 and t2=1, melting and Z-DNA
- // When [0][2], means you have t1=0 and t2=2, melting and cruciforms
- // When [1][2], means you have t1=1 and t2=2, Z-DNA and cruciforms.
- // Cases when [t][t] means having multiple runs of the same transition
-
- // V: Also these arrays were added for the competition branch
- // t=0 for melting, t=1 for Z-DNA, and t=2 for cruciforms
- double Prob_conditional[3][3]; // V: This should be the conditional probability
- double Prob_conditional_not[3][3]; // V: This should be the conditional probability of having [A] but not [B]
+ double exp_type_B;
+ // V: This Prob_compete parameter was added for the competition branch
+ // t=0 for melting, t=1 for Z-DNA, and t=2 for cruciforms
+ double Prob_compete[3]
+ [3]; // V: This should be the probability of intersection?
+ // When [0][1], means you have runs of t1=0 and t2=1, melting and Z-DNA
+ // When [0][2], means you have t1=0 and t2=2, melting and cruciforms
+ // When [1][2], means you have t1=1 and t2=2, Z-DNA and cruciforms.
+ // Cases when [t][t] means having multiple runs of the same transition
+ // V: Also these arrays were added for the competition branch
+ // t=0 for melting, t=1 for Z-DNA, and t=2 for cruciforms
+ double Prob_conditional[3]
+ [3]; // V: This should be the conditional probability
+ double Prob_conditional_not[3][3]; // V: This should be the conditional
+ // probability of having [A] but not [B]
};
#endif // !defined(AFX_SIDD_1R_H__2688B961_9B45_44D5_8F6F_F45636E04460__INCLUDED_)
diff --git a/src/trans_compete/SIDD_2R.cpp b/src/trans_compete/SIDD_2R.cpp
index b642eb1..c91c620 100644
--- a/src/trans_compete/SIDD_2R.cpp
+++ b/src/trans_compete/SIDD_2R.cpp
@@ -18,184 +18,188 @@
#include "SIDD_2R.h"
-#include
#include
+#include
//////////////////////////////////////////////////////////////////////
// Construction/Destruction
//////////////////////////////////////////////////////////////////////
-SIDD_2R::SIDD_2R():SIDD_1R()
-{
- flag_minE2 = false;
- StatesInTwoRuns = 0;
- sum_2RB = sum_2RG = 0.0;
-
+SIDD_2R::SIDD_2R() : SIDD_1R() {
+ flag_minE2 = false;
+ StatesInTwoRuns = 0;
+ sum_2RB = sum_2RG = 0.0;
}
-SIDD_2R::~SIDD_2R()
-{
-}
+SIDD_2R::~SIDD_2R() {}
// search all lowest states for two runs
-bool SIDD_2R::Search_Low2RE()
-{
- int par[3];
- long count = 0;
- for (int t1=0; t1<=Ns; t1++) {
- for (int t2=0; t2<=Ns; t2++) {
- exp_two[t1][t2] = 0;
-
- // V: This was added for the probability of having at least two runs of two transitions
- Prob_compete[t1][t2]=0;
- // V: These will be the conditionals
- Prob_conditional[t1][t2]=0;
- Prob_conditional_not[t1][t2]=0;
-
- }
+bool SIDD_2R::Search_Low2RE() {
+ int par[3];
+ long count = 0;
+ for (int t1 = 0; t1 <= Ns; t1++) {
+ for (int t2 = 0; t2 <= Ns; t2++) {
+ exp_two[t1][t2] = 0;
+
+ // V: This was added for the probability of having at least two runs of
+ // two transitions
+ Prob_compete[t1][t2] = 0;
+ // V: These will be the conditionals
+ Prob_conditional[t1][t2] = 0;
+ Prob_conditional_not[t1][t2] = 0;
}
- LstState::iterator i;
- LstState::iterator j;
- flag_minE2 = false; // default value
-
- sum_2RB = sum_2RG = 0.0;
- for (int t1=0; t1<=Ns; t1++) {
- for (int t2=0; t2<=Ns; t2++) {
- for(int n1 = MinWindowSize; n1 <= MaxWindowSize/2; n1++){
- if(Lst_OBE[n1][t1].size() < 1)
- continue;
- if(Lst_OBE[n1][t1].begin()->get_energy() + min_RE[t2] + minGres >= max_E)
- continue;
- for(int n2 = n1; n2 <= MaxWindowSize - n1; n2++){
- if (n1==n2 && t1==1 && t2==0) continue;
- if (n1==n2 && t1==2 && t2==0) continue;
- if (n1==n2 && t1==2 && t2==1) continue;
- if(Lst_OBE[n2][t2].size() < 1) continue;
- for (int p=0; p<=2; p++)
- par[p]=n1*delta_fnc(t1,p)+n2*delta_fnc(t2,p);
- double Gr=calc_Gres(par[0],par[1],par[2],t1*delta_fnc(t1,1)+t2*delta_fnc(t2,1));
- if(Lst_OBE[n1][t1].begin()->get_energy() + Lst_OBE[n2][t2].begin()->get_energy() + Gr >= max_E)
- continue;
- for(i = Lst_OBE[n1][t1].begin(); i != Lst_OBE[n1][t1].end(); ++i){
- int p1 = i->get_pos1();
- double e1 = i->get_energy();
- if (e1 >= 10000) break;
- j = Lst_OBE[n2][t2].begin();
- if(n1 == n2 && t1==t2) ++j; // no repeat of the same group
- if (e1 + j->get_energy() + Gr >= max_E)
- break;
- for(; j != Lst_OBE[n2][t2].end(); ++j){
- int p2 = j->get_pos1();
- double e2 = j->get_energy();
- if (e2 >= 10000) break;
- double e = e1 + e2 + Gr; //total energy for two runs
- if (e >= max_E) break;
- if(!overlap(p1, p2, n1, n2)){
- if(e < min_E){
- min_E = e;
- max_E = min_E + theta;
- flag_minE2 = true;
- }
- if(e < max_E){
- if(write_profile){
- // double exponent= exp(-e/RT);
- // V: the scaling factor is added
- double exponent= exp(-e/RT + scaling_factor);
-
- exp_two[t1][t2] += exponent;
- for (int t=0; t<=Ns; t++)
- runs[t] +=exponent*(delta_fnc(t,t1)+delta_fnc(t,t2));
- update_promatrix(p1, n1, t1, e, exponent);
- update_promatrix(p2, n2, t2, e, exponent);
- sum_2RB += exponent;
- sum_2RG = sum_2RG + e*exponent;
- }
- count++;
-
- }
- }
- }
- }
- }
- }
- }
- }
-
- StatesInTwoRuns = count;
- ZsumB += sum_2RB;
- ZsumG +=sum_2RG;
- if (results)
- cout << "Number of two-run states = " << count << endl;
- for (int t=0; t<=Ns; t++) {
- for (int t1=0; t1<=Ns; t1++) {
- for (int t2=0; t2<=Ns; t2++) {
- if (t1==t || t2==t)
- Prob[t] += exp_two[t1][t2];
+ }
+ LstState::iterator i;
+ LstState::iterator j;
+ flag_minE2 = false; // default value
+
+ sum_2RB = sum_2RG = 0.0;
+ for (int t1 = 0; t1 <= Ns; t1++) {
+ for (int t2 = 0; t2 <= Ns; t2++) {
+ for (int n1 = MinWindowSize; n1 <= MaxWindowSize / 2; n1++) {
+ if (Lst_OBE[n1][t1].size() < 1)
+ continue;
+ if (Lst_OBE[n1][t1].begin()->get_energy() + min_RE[t2] + minGres >=
+ max_E)
+ continue;
+ for (int n2 = n1; n2 <= MaxWindowSize - n1; n2++) {
+ if (n1 == n2 && t1 == 1 && t2 == 0)
+ continue;
+ if (n1 == n2 && t1 == 2 && t2 == 0)
+ continue;
+ if (n1 == n2 && t1 == 2 && t2 == 1)
+ continue;
+ if (Lst_OBE[n2][t2].size() < 1)
+ continue;
+ for (int p = 0; p <= 2; p++)
+ par[p] = n1 * delta_fnc(t1, p) + n2 * delta_fnc(t2, p);
+ double Gr = calc_Gres(par[0], par[1], par[2],
+ t1 * delta_fnc(t1, 1) + t2 * delta_fnc(t2, 1));
+ if (Lst_OBE[n1][t1].begin()->get_energy() +
+ Lst_OBE[n2][t2].begin()->get_energy() + Gr >=
+ max_E)
+ continue;
+ for (i = Lst_OBE[n1][t1].begin(); i != Lst_OBE[n1][t1].end(); ++i) {
+ int p1 = i->get_pos1();
+ double e1 = i->get_energy();
+ if (e1 >= 10000)
+ break;
+ j = Lst_OBE[n2][t2].begin();
+ if (n1 == n2 && t1 == t2)
+ ++j; // no repeat of the same group
+ if (e1 + j->get_energy() + Gr >= max_E)
+ break;
+ for (; j != Lst_OBE[n2][t2].end(); ++j) {
+ int p2 = j->get_pos1();
+ double e2 = j->get_energy();
+ if (e2 >= 10000)
+ break;
+ double e = e1 + e2 + Gr; // total energy for two runs
+ if (e >= max_E)
+ break;
+ if (!overlap(p1, p2, n1, n2)) {
+ if (e < min_E) {
+ min_E = e;
+ max_E = min_E + theta;
+ flag_minE2 = true;
+ }
+ if (e < max_E) {
+ if (write_profile) {
+ // double exponent= exp(-e/RT);
+ // V: the scaling factor is added
+ double exponent = exp(-e / RT + scaling_factor);
+
+ exp_two[t1][t2] += exponent;
+ for (int t = 0; t <= Ns; t++)
+ runs[t] +=
+ exponent * (delta_fnc(t, t1) + delta_fnc(t, t2));
+ update_promatrix(p1, n1, t1, e, exponent);
+ update_promatrix(p2, n2, t2, e, exponent);
+ sum_2RB += exponent;
+ sum_2RG = sum_2RG + e * exponent;
+ }
+ count++;
+ }
+ }
}
+ }
}
+ }
}
- // V: This counts for states where you have two different structural transitions, ti and tj.
- // V: Watch out for the diagonal, it doesn't tell you the prob of having one run of another of the same transition. It doesn't even give you p(M)
- for (int ti=0; ti<=Ns; ti++) {
- for (int tj=0; tj<=Ns; tj++) {
- for (int t1=0; t1<=Ns; t1++) {
- for (int t2=0; t2<=Ns; t2++) {
-
- // Let's build vectors of type of transitions, then check if these transitions are included
- std::vector firstVec = { t1, t2};
- std::vector secondVec = { ti, tj};
-
- // Sort first vector
- std::sort(firstVec.begin(), firstVec.end());
- // Sort second vector
- std::sort(secondVec.begin(), secondVec.end());
-
- // Check if all elements of a second vector exists in first vector
- bool my_condition = std::includes(firstVec.begin(), firstVec.end(), secondVec.begin(), secondVec.end());
-
- if (my_condition) {
- Prob_compete[ti][tj] += exp_two[t1][t2];
- }
- }
- }
+ }
+
+ StatesInTwoRuns = count;
+ ZsumB += sum_2RB;
+ ZsumG += sum_2RG;
+ if (results)
+ cout << "Number of two-run states = " << count << endl;
+ for (int t = 0; t <= Ns; t++) {
+ for (int t1 = 0; t1 <= Ns; t1++) {
+ for (int t2 = 0; t2 <= Ns; t2++) {
+ if (t1 == t || t2 == t)
+ Prob[t] += exp_two[t1][t2];
+ }
+ }
+ }
+ // V: This counts for states where you have two different structural
+ // transitions, ti and tj. V: Watch out for the diagonal, it doesn't tell you
+ // the prob of having one run of another of the same transition. It doesn't
+ // even give you p(M)
+ for (int ti = 0; ti <= Ns; ti++) {
+ for (int tj = 0; tj <= Ns; tj++) {
+ for (int t1 = 0; t1 <= Ns; t1++) {
+ for (int t2 = 0; t2 <= Ns; t2++) {
+
+ // Let's build vectors of type of transitions, then check if these
+ // transitions are included
+ std::vector firstVec = {t1, t2};
+ std::vector secondVec = {ti, tj};
+
+ // Sort first vector
+ std::sort(firstVec.begin(), firstVec.end());
+ // Sort second vector
+ std::sort(secondVec.begin(), secondVec.end());
+
+ // Check if all elements of a second vector exists in first vector
+ bool my_condition = std::includes(firstVec.begin(), firstVec.end(),
+ secondVec.begin(), secondVec.end());
+
+ if (my_condition) {
+ Prob_compete[ti][tj] += exp_two[t1][t2];
+ }
}
+ }
}
+ }
- return flag_minE2;
+ return flag_minE2;
}
-
-
-
// checking if two regions are overlapping
-bool SIDD_2R::overlap(int p1, int p2, int r1, int r2)
-{
-
- if(r1 + r2 > MaxWindowSize){
- cerr << "out of window size.\n";
- return false;
- }
-
- if(p2 + r2 <= length_seq - 1 && p1 + r1 <= length_seq - 1){
- if(p1 + r1 < p2 - 1 || p2 + r2 < p1 - 1) // at least one gap
- return false; // no overlapping
- else
- return true;
- }
- else if(p2 + r2 > length_seq - 1 && p1 + r1 <= length_seq - 1){
- if(p1 + r1 < p2 - 1 && r2 - (length_seq - 1 - p2) < p1 - 1)
- return false; // no overlapping
- else
- return true;
- }
- else if(p1 + r1 > length_seq - 1 && p2 + r2 <= length_seq - 1){
- if(p2 + r2 < p1 - 1 && r1 - (length_seq - 1 - p1) < p2 - 1)
- return false; // no overlapping
- else
- return true;
- }
- else
- return true; // overlapping
+bool SIDD_2R::overlap(int p1, int p2, int r1, int r2) {
+
+ if (r1 + r2 > MaxWindowSize) {
+ cerr << "out of window size.\n";
+ return false;
+ }
+
+ if (p2 + r2 <= length_seq - 1 && p1 + r1 <= length_seq - 1) {
+ if (p1 + r1 < p2 - 1 || p2 + r2 < p1 - 1) // at least one gap
+ return false; // no overlapping
+ else
+ return true;
+ } else if (p2 + r2 > length_seq - 1 && p1 + r1 <= length_seq - 1) {
+ if (p1 + r1 < p2 - 1 && r2 - (length_seq - 1 - p2) < p1 - 1)
+ return false; // no overlapping
+ else
+ return true;
+ } else if (p1 + r1 > length_seq - 1 && p2 + r2 <= length_seq - 1) {
+ if (p2 + r2 < p1 - 1 && r1 - (length_seq - 1 - p1) < p2 - 1)
+ return false; // no overlapping
+ else
+ return true;
+ } else
+ return true; // overlapping
}
#endif
diff --git a/src/trans_compete/SIDD_2R.h b/src/trans_compete/SIDD_2R.h
index 863cb02..38473ff 100644
--- a/src/trans_compete/SIDD_2R.h
+++ b/src/trans_compete/SIDD_2R.h
@@ -17,22 +17,18 @@
#include "SIDD_1R.h"
-
-class SIDD_2R : public SIDD_1R
-{
+class SIDD_2R : public SIDD_1R {
protected:
-
- bool flag_minE2;
- long StatesInTwoRuns;
+ bool flag_minE2;
+ long StatesInTwoRuns;
public:
- SIDD_2R();
- virtual ~SIDD_2R();
-
- // two-run states
- bool overlap(int p1, int p2, int r1, int r2);
- bool Search_Low2RE();
+ SIDD_2R();
+ virtual ~SIDD_2R();
+ // two-run states
+ bool overlap(int p1, int p2, int r1, int r2);
+ bool Search_Low2RE();
};
#endif // !defined(AFX_SIDD_2R_H__79785266_E18C_42AD_8C80_6BC6C04A20DD__INCLUDED_)
diff --git a/src/trans_compete/SIDD_3R.cpp b/src/trans_compete/SIDD_3R.cpp
index 986147a..9f6c4db 100644
--- a/src/trans_compete/SIDD_3R.cpp
+++ b/src/trans_compete/SIDD_3R.cpp
@@ -18,160 +18,188 @@
#include "SIDD_3R.h"
-#include
#include
-
+#include
//////////////////////////////////////////////////////////////////////
// Construction/Destruction
//////////////////////////////////////////////////////////////////////
-SIDD_3R::SIDD_3R():SIDD_2R()
-{
-
- flag_minE3 = false;
- StatesInThreeRuns = 0;
- sum_3RB = sum_3RG = 0.0;
-}
+SIDD_3R::SIDD_3R() : SIDD_2R() {
-SIDD_3R::~SIDD_3R()
-{
+ flag_minE3 = false;
+ StatesInThreeRuns = 0;
+ sum_3RB = sum_3RG = 0.0;
}
+SIDD_3R::~SIDD_3R() {}
+
// search all lowest states for three runs
-bool SIDD_3R::Search_Low3RE()
-{
-
- double count = 0.0;
- int par[3];
- for (int t1=0; t1<=Ns; t1++) {
- for (int t2=0; t2<=Ns; t2++) {
- for (int t3=0; t3<=Ns; t3++) {
- exp_three[t1][t2][t3] = 0;
- }
- }
+bool SIDD_3R::Search_Low3RE() {
+
+ double count = 0.0;
+ int par[3];
+ for (int t1 = 0; t1 <= Ns; t1++) {
+ for (int t2 = 0; t2 <= Ns; t2++) {
+ for (int t3 = 0; t3 <= Ns; t3++) {
+ exp_three[t1][t2][t3] = 0;
+ }
}
- LstState::iterator i;
- LstState::iterator k;
- LstState::iterator j;
- flag_minE3 = false;
- sum_3RB = sum_3RG = 0.0;
- int window_size = 200;
- double Epsilon = 0.0;
- for (int t1=0; t1<=2; t1++) {
- for (int t2=0; t2<=2; t2++) {
- for (int t3=0; t3<=2; t3++) {
- for(int n1 = 1; n1 <= window_size/3; n1++){
- if(Lst_OBE[n1][t1].begin()->get_energy() + min_RE[t2]+min_RE[t3] + minGres >= max_E)
- continue;
- for(int n2 = n1; n2 <= (window_size - n1)/2 ; n2++){
- if (n1==n2 && t1==1 && t2==0) continue;
- if(Lst_OBE[n1][t1].begin()->get_energy() + Lst_OBE[n2][t2].begin()->get_energy()+ min_RE[t3] + minGres >= max_E)
- continue;
- for(int n3 = n2; n3 <= window_size - n1 - n2; n3++){
- if (n1==n3 && t1==1 && t3==0) continue;
- if (n2==n3 && t2==1 && t3==0) continue;
- for (int p=0; p<=2; p++)
- par[p]=n1*delta_fnc(t1,p)+n2*delta_fnc(t2,p)+n3*delta_fnc(t3,p);
- double Gr=calc_Gres(par[0],par[1],par[2],t1*delta_fnc(t1,1)+t2*delta_fnc(t2,1)+t3*delta_fnc(t3,1));
- if(Lst_OBE[n1][t1].begin()->get_energy() + Lst_OBE[n2][t2].begin()->get_energy() + Lst_OBE[n3][t3].begin()->get_energy() + Gr>= max_E)
- continue;
- for(i = Lst_OBE[n1][t1].begin(); i != Lst_OBE[n1][t1].end(); ++i){
- int p1 = i->get_pos1();
- double e1 = i->get_energy();
- if (e1 >= 10000) break;
- j = Lst_OBE[n2][t2].begin();
- if(n1 == n2 && t1==t2) ++j; // no repeat of the same group
- if(e1 + j->get_energy() + Lst_OBE[n3][t3].begin()->get_energy()+Gr >= max_E)
- break;
- for(; j != Lst_OBE[n2][t2].end(); ++j){
- int p2 = j->get_pos1();
- double e2 = j->get_energy();
- if (e2 >= 10000) break;
- k = Lst_OBE[n3][t3].begin();
- if(n2 == n3 && t2==t3) ++k;
- if(e1 + e2 + k->get_energy() + Gr >= max_E)
- break;
- for(; k != Lst_OBE[n3][t3].end(); ++k){
- int p3 = k->get_pos1();
- double e3 = k->get_energy();
- if (e3 >= 10000) break;
- double e = e1 + e2 + e3 + Gr;
- if(e >= max_E) break;
- if(!overlap(p1, p2, n1, n2) && !overlap(p1, p3, n1, n3) && !overlap(p2, p3, n2, n3)){
- if(min_E - e >=Epsilon){
- min_E = e;
- max_E = e + theta;
- }
- count++;
- if(write_profile){
- // double exponent = exp(-e/RT);
- double exponent= exp(-e/RT + scaling_factor); // V: scaling factor added
- exp_three[t1][t2][t3] +=exponent;
- for (int t=0; t<=Ns; t++)
- runs[t] +=exponent*(delta_fnc(t,t1)+delta_fnc(t,t2)+delta_fnc(t,t3));
- update_promatrix(p1, n1, t1, e, exponent);
- update_promatrix(p2, n2, t2, e, exponent);
- update_promatrix(p3, n3, t3, e, exponent);
- sum_3RB += exponent;
- sum_3RG = sum_3RG + e*exponent;
- }
- }
- }
- }
- }
- }
- }
- }
- }
- }
- }
-
- if (results)
- cout << "Number of three-run states = " << count << endl;
- StatesInThreeRuns = (long)count;
- ZsumB += sum_3RB;
- ZsumG += sum_3RG;
- for (int t=0; t<=Ns; t++) {
- for (int t1=0; t1<=Ns; t1++) {
- for (int t2=0; t2<=Ns; t2++) {
- for (int t3=0; t3<=Ns; t3++) {
- if (t1==t || t2==t || t3==t)
- Prob[t]+=exp_three[t1][t2][t3];
+ }
+ LstState::iterator i;
+ LstState::iterator k;
+ LstState::iterator j;
+ flag_minE3 = false;
+ sum_3RB = sum_3RG = 0.0;
+ int window_size = 200;
+ double Epsilon = 0.0;
+ for (int t1 = 0; t1 <= 2; t1++) {
+ for (int t2 = 0; t2 <= 2; t2++) {
+ for (int t3 = 0; t3 <= 2; t3++) {
+ for (int n1 = 1; n1 <= window_size / 3; n1++) {
+ if (Lst_OBE[n1][t1].begin()->get_energy() + min_RE[t2] + min_RE[t3] +
+ minGres >=
+ max_E)
+ continue;
+ for (int n2 = n1; n2 <= (window_size - n1) / 2; n2++) {
+ if (n1 == n2 && t1 == 1 && t2 == 0)
+ continue;
+ if (Lst_OBE[n1][t1].begin()->get_energy() +
+ Lst_OBE[n2][t2].begin()->get_energy() + min_RE[t3] +
+ minGres >=
+ max_E)
+ continue;
+ for (int n3 = n2; n3 <= window_size - n1 - n2; n3++) {
+ if (n1 == n3 && t1 == 1 && t3 == 0)
+ continue;
+ if (n2 == n3 && t2 == 1 && t3 == 0)
+ continue;
+ for (int p = 0; p <= 2; p++)
+ par[p] = n1 * delta_fnc(t1, p) + n2 * delta_fnc(t2, p) +
+ n3 * delta_fnc(t3, p);
+ double Gr =
+ calc_Gres(par[0], par[1], par[2],
+ t1 * delta_fnc(t1, 1) + t2 * delta_fnc(t2, 1) +
+ t3 * delta_fnc(t3, 1));
+ if (Lst_OBE[n1][t1].begin()->get_energy() +
+ Lst_OBE[n2][t2].begin()->get_energy() +
+ Lst_OBE[n3][t3].begin()->get_energy() + Gr >=
+ max_E)
+ continue;
+ for (i = Lst_OBE[n1][t1].begin(); i != Lst_OBE[n1][t1].end();
+ ++i) {
+ int p1 = i->get_pos1();
+ double e1 = i->get_energy();
+ if (e1 >= 10000)
+ break;
+ j = Lst_OBE[n2][t2].begin();
+ if (n1 == n2 && t1 == t2)
+ ++j; // no repeat of the same group
+ if (e1 + j->get_energy() +
+ Lst_OBE[n3][t3].begin()->get_energy() + Gr >=
+ max_E)
+ break;
+ for (; j != Lst_OBE[n2][t2].end(); ++j) {
+ int p2 = j->get_pos1();
+ double e2 = j->get_energy();
+ if (e2 >= 10000)
+ break;
+ k = Lst_OBE[n3][t3].begin();
+ if (n2 == n3 && t2 == t3)
+ ++k;
+ if (e1 + e2 + k->get_energy() + Gr >= max_E)
+ break;
+ for (; k != Lst_OBE[n3][t3].end(); ++k) {
+ int p3 = k->get_pos1();
+ double e3 = k->get_energy();
+ if (e3 >= 10000)
+ break;
+ double e = e1 + e2 + e3 + Gr;
+ if (e >= max_E)
+ break;
+ if (!overlap(p1, p2, n1, n2) && !overlap(p1, p3, n1, n3) &&
+ !overlap(p2, p3, n2, n3)) {
+ if (min_E - e >= Epsilon) {
+ min_E = e;
+ max_E = e + theta;
+ }
+ count++;
+ if (write_profile) {
+ // double exponent = exp(-e/RT);
+ double exponent =
+ exp(-e / RT +
+ scaling_factor); // V: scaling factor added
+ exp_three[t1][t2][t3] += exponent;
+ for (int t = 0; t <= Ns; t++)
+ runs[t] +=
+ exponent * (delta_fnc(t, t1) + delta_fnc(t, t2) +
+ delta_fnc(t, t3));
+ update_promatrix(p1, n1, t1, e, exponent);
+ update_promatrix(p2, n2, t2, e, exponent);
+ update_promatrix(p3, n3, t3, e, exponent);
+ sum_3RB += exponent;
+ sum_3RG = sum_3RG + e * exponent;
+ }
+ }
+ }
}
+ }
}
+ }
+ }
+ }
+ }
+ }
+
+ if (results)
+ cout << "Number of three-run states = " << count << endl;
+ StatesInThreeRuns = (long)count;
+ ZsumB += sum_3RB;
+ ZsumG += sum_3RG;
+ for (int t = 0; t <= Ns; t++) {
+ for (int t1 = 0; t1 <= Ns; t1++) {
+ for (int t2 = 0; t2 <= Ns; t2++) {
+ for (int t3 = 0; t3 <= Ns; t3++) {
+ if (t1 == t || t2 == t || t3 == t)
+ Prob[t] += exp_three[t1][t2][t3];
}
+ }
}
- // V: This was added for the competition of transitions:
- for (int ti=0; ti<=Ns; ti++) {
- for (int tj=0; tj<=Ns; tj++) {
- for (int t1=0; t1<=Ns; t1++) {
- for (int t2=0; t2<=Ns; t2++) {
- for (int t3=0; t3<=Ns; t3++) {
+ }
+ // V: This was added for the competition of transitions:
+ for (int ti = 0; ti <= Ns; ti++) {
+ for (int tj = 0; tj <= Ns; tj++) {
+ for (int t1 = 0; t1 <= Ns; t1++) {
+ for (int t2 = 0; t2 <= Ns; t2++) {
+ for (int t3 = 0; t3 <= Ns; t3++) {
- // Let's build vectors of type of transitions, then check if these transitions are included
- std::vector firstVec = { t1, t2, t3};
- std::vector secondVec = { ti, tj};
+ // Let's build vectors of type of transitions, then check if these
+ // transitions are included
+ std::vector firstVec = {t1, t2, t3};
+ std::vector secondVec = {ti, tj};
- // Sort first vector
- std::sort(firstVec.begin(), firstVec.end());
- // Sort second vector
- std::sort(secondVec.begin(), secondVec.end());
+ // Sort first vector
+ std::sort(firstVec.begin(), firstVec.end());
+ // Sort second vector
+ std::sort(secondVec.begin(), secondVec.end());
- // Check if all elements of a second vector exists in first vector
- bool my_condition = std::includes(firstVec.begin(), firstVec.end(), secondVec.begin(), secondVec.end());
+ // Check if all elements of a second vector exists in first vector
+ bool my_condition =
+ std::includes(firstVec.begin(), firstVec.end(),
+ secondVec.begin(), secondVec.end());
- if (my_condition) {
- Prob_compete[ti][tj] += exp_three[t1][t2][t3];
- //cout << " ti "<
#include
#include
+#include
//////////////////////////////////////////////////////////////////////
// Construction/Destruction
//////////////////////////////////////////////////////////////////////
-SIDD_4R::SIDD_4R():SIDD_3R()
-{
-
- flag_minE4 = false;
- StatesInFourRuns = 0;
- sum_4RB = sum_4RG = 0.0;
-}
+SIDD_4R::SIDD_4R() : SIDD_3R() {
-SIDD_4R::~SIDD_4R()
-{
+ flag_minE4 = false;
+ StatesInFourRuns = 0;
+ sum_4RB = sum_4RG = 0.0;
}
+SIDD_4R::~SIDD_4R() {}
+
// search all lowest states for four runs
-bool SIDD_4R::Search_Low4RE()
-{
- double count = 0.0;
- int par[3];
- for (int t1=0; t1<=Ns; t1++) {
- for (int t2=0; t2<=Ns; t2++) {
- for (int t3=0; t3<=Ns; t3++) {
- for (int t4=0; t4<=Ns; t4++) {
- exp_four[t1][t2][t3][t4] = 0;
- }
- }
+bool SIDD_4R::Search_Low4RE() {
+ double count = 0.0;
+ int par[3];
+ for (int t1 = 0; t1 <= Ns; t1++) {
+ for (int t2 = 0; t2 <= Ns; t2++) {
+ for (int t3 = 0; t3 <= Ns; t3++) {
+ for (int t4 = 0; t4 <= Ns; t4++) {
+ exp_four[t1][t2][t3][t4] = 0;
}
+ }
}
- LstState::iterator i;
- LstState::iterator k;
- LstState::iterator j;
- LstState::iterator l;
- flag_minE4 = false;
- sum_4RB = sum_4RG = 0.0;
- int window_size = 200;
- double Epsilon = 0.0;
- for (int t1=0; t1<=2; t1++) {
- for (int t2=0; t2<=2; t2++) {
- for (int t3=0; t3<=2; t3++) {
- for (int t4=0; t4<=1; t4++) {
- for(int n1 = 1; n1 <= window_size/4; n1++){
- //start here
- if(Lst_OBE[n1][t1].begin()->get_energy() + min_RE[t2] + min_RE[t3]+ min_RE[t4] + minGres >= max_E)
- continue;
- for(int n2 = n1; n2 <= (window_size - n1)/3 ; n2++){
- if (n1==n2 && t1==1 && t2==0)
- continue;
- if(Lst_OBE[n1][t1].begin()->get_energy() + Lst_OBE[n2][t2].begin()->get_energy()+ min_RE[t3] + min_RE[t4] + minGres >= max_E)
- continue;
- for(int n3 = n2; n3 <= (window_size - n1 - n2)/2; n3++){
- if (n1==n3 && t1==1 && t3==0) continue;
- if (n2==n3 && t2==1 && t3==0) continue;
- if(Lst_OBE[n1][t1].begin()->get_energy() + Lst_OBE[n2][t2].begin()->get_energy() + Lst_OBE[n3][t3].begin()->get_energy() + min_RE[t4] + minGres >= max_E)
- continue;
- for(int n4 = n3; n4 <= (window_size - n1 - n2-n3); n4++){
- if (n1==n4 && t1==1 && t4==0) continue;
- if (n2==n4 && t2==1 && t4==0) continue;
- if (n3==n4 && t3==1 && t4==0) continue;
- for (int p=0; p<=2; p++)
- par[p]=n1*delta_fnc(t1,p)+n2*delta_fnc(t2,p)+n3*delta_fnc(t3,p)+n4*delta_fnc(t4,p);
- double Gr=calc_Gres(par[0],par[1],par[2],t1*delta_fnc(t1,1)+t2*delta_fnc(t2,1)+t3*delta_fnc(t3,1)+t4*delta_fnc(t4,1));
- if(Lst_OBE[n1][t1].begin()->get_energy() +Lst_OBE[n2][t2].begin()->get_energy() + Lst_OBE[n3][t3].begin()->get_energy() + Gr>= max_E)
- continue;
- for(i = Lst_OBE[n1][t1].begin(); i != Lst_OBE[n1][t1].end(); ++i){
- int p1 = i->get_pos1();
- double e1 = i->get_energy();
- if (e1 >= 10000)
- break;
- j = Lst_OBE[n2][t2].begin();
- if(n1 == n2 && t1==t2) ++j; // no repeat of the same group
- if(e1 + j->get_energy() + Lst_OBE[n3][t3].begin()->get_energy() + Lst_OBE[n4][t4].begin()->get_energy()+ Gr>= max_E)
- break;
- for(; j != Lst_OBE[n2][t3].end(); ++j){
- int p2 = j->get_pos1();
- double e2 = j->get_energy();
- if (e2 >= 10000) break;
- k = Lst_OBE[n3][t3].begin();
- if(n2 == n3 && t2==t3) ++k;
- if(e1 + e2 + k->get_energy() + Lst_OBE[n4][t4].begin()->get_energy()+Gr >= max_E)
- break;
- for(; k != Lst_OBE[n3][t3].end(); ++k){
- int p3 = k->get_pos1();
- double e3 = k->get_energy();
- if (e3 >= 10000) break;
- l = Lst_OBE[n4][t4].begin();
- if(n3 == n4 && t3==t4) ++l;
- if(e1 + e2 + e3 + l->get_energy() +Gr >= max_E)
- break;
- for(; l != Lst_OBE[n4][t4].end(); ++l){
- int p4 = l->get_pos1();
- double e4 = l->get_energy();
- if (e4 >= 10000) break;
- double e = e1 + e2 + e3 + e4 +Gr;
- if(e >= max_E) break;
- if(!overlap(p1, p2, n1, n2) && !overlap(p1, p3, n1, n3) && !overlap(p1, p4, n1, n4) && !overlap(p2, p3, n2, n3) && !overlap(p2, p4, n2, n4) && !overlap(p3, p4, n3, n4)){
- if(min_E - e >=Epsilon){
- min_E = e;
- max_E = e + theta;
- }
- count++;
-
- if(write_profile){
- //double exponent = exp(-e/RT);
- double exponent= exp(-e/RT + scaling_factor); // V: scaling factor added
-
- exp_four[t1][t2][t3][t4]+=exponent;
- for (int t=0; t<=Ns; t++)
- runs[t] +=exponent*(delta_fnc(t,t1)+delta_fnc(t,t2)+delta_fnc(t,t3)+delta_fnc(t,t4));
- update_promatrix(p1, n1, t1, e, exponent);
- update_promatrix(p2, n2, t2, e, exponent);
- update_promatrix(p3, n3, t3, e, exponent);
- update_promatrix(p4, n4, t4, e, exponent);
- sum_4RB += exponent;
- sum_4RG = sum_4RG + e*exponent;
-
- }
- }
- }
- }
- }
- }
- }
- }
- }
- }
- }
- }
- }
- }
-
- StatesInFourRuns = (long)count;
- ZsumB += sum_4RB;
- ZsumG += sum_4RG;
- totalFreq = StatesInOneRun + StatesInTwoRuns + StatesInThreeRuns+StatesInFourRuns;
- if (results) {
- cout << "Number of four-run states = " << count << endl;
- cout << "Total number of states = " << totalFreq << endl;
-
- // cout << "Total Partition Function without exp= " << ZsumB << endl;
- // V: The scaling factor needs to be removed from the partition function
- cout << "Total Partition Function = " << exp(-scaling_factor)*ZsumB << endl;
- cout << "Scaling factor " << scaling_factor << endl;
-
- cout << "Number of M_runs = " << runs[0]/ZsumB << endl;
- cout<< "Number of Z_runs = " << runs[1]/ZsumB << endl;
- cout<< "Number of C_runs = " << runs[2]/ZsumB << endl;
- double runs= (sum_1RB+2*sum_2RB+3*sum_3RB+4*sum_4RB)/ZsumB;
- cout << "Average number of runs = " << runs << endl;
- }
+ }
+ LstState::iterator i;
+ LstState::iterator k;
+ LstState::iterator j;
+ LstState::iterator l;
+ flag_minE4 = false;
+ sum_4RB = sum_4RG = 0.0;
+ int window_size = 200;
+ double Epsilon = 0.0;
+ for (int t1 = 0; t1 <= 2; t1++) {
+ for (int t2 = 0; t2 <= 2; t2++) {
+ for (int t3 = 0; t3 <= 2; t3++) {
+ for (int t4 = 0; t4 <= 1; t4++) {
+ for (int n1 = 1; n1 <= window_size / 4; n1++) {
+ // start here
+ if (Lst_OBE[n1][t1].begin()->get_energy() + min_RE[t2] +
+ min_RE[t3] + min_RE[t4] + minGres >=
+ max_E)
+ continue;
+ for (int n2 = n1; n2 <= (window_size - n1) / 3; n2++) {
+ if (n1 == n2 && t1 == 1 && t2 == 0)
+ continue;
+ if (Lst_OBE[n1][t1].begin()->get_energy() +
+ Lst_OBE[n2][t2].begin()->get_energy() + min_RE[t3] +
+ min_RE[t4] + minGres >=
+ max_E)
+ continue;
+ for (int n3 = n2; n3 <= (window_size - n1 - n2) / 2; n3++) {
+ if (n1 == n3 && t1 == 1 && t3 == 0)
+ continue;
+ if (n2 == n3 && t2 == 1 && t3 == 0)
+ continue;
+ if (Lst_OBE[n1][t1].begin()->get_energy() +
+ Lst_OBE[n2][t2].begin()->get_energy() +
+ Lst_OBE[n3][t3].begin()->get_energy() + min_RE[t4] +
+ minGres >=
+ max_E)
+ continue;
+ for (int n4 = n3; n4 <= (window_size - n1 - n2 - n3); n4++) {
+ if (n1 == n4 && t1 == 1 && t4 == 0)
+ continue;
+ if (n2 == n4 && t2 == 1 && t4 == 0)
+ continue;
+ if (n3 == n4 && t3 == 1 && t4 == 0)
+ continue;
+ for (int p = 0; p <= 2; p++)
+ par[p] = n1 * delta_fnc(t1, p) + n2 * delta_fnc(t2, p) +
+ n3 * delta_fnc(t3, p) + n4 * delta_fnc(t4, p);
+ double Gr = calc_Gres(
+ par[0], par[1], par[2],
+ t1 * delta_fnc(t1, 1) + t2 * delta_fnc(t2, 1) +
+ t3 * delta_fnc(t3, 1) + t4 * delta_fnc(t4, 1));
+ if (Lst_OBE[n1][t1].begin()->get_energy() +
+ Lst_OBE[n2][t2].begin()->get_energy() +
+ Lst_OBE[n3][t3].begin()->get_energy() + Gr >=
+ max_E)
+ continue;
+ for (i = Lst_OBE[n1][t1].begin(); i != Lst_OBE[n1][t1].end();
+ ++i) {
+ int p1 = i->get_pos1();
+ double e1 = i->get_energy();
+ if (e1 >= 10000)
+ break;
+ j = Lst_OBE[n2][t2].begin();
+ if (n1 == n2 && t1 == t2)
+ ++j; // no repeat of the same group
+ if (e1 + j->get_energy() +
+ Lst_OBE[n3][t3].begin()->get_energy() +
+ Lst_OBE[n4][t4].begin()->get_energy() + Gr >=
+ max_E)
+ break;
+ for (; j != Lst_OBE[n2][t3].end(); ++j) {
+ int p2 = j->get_pos1();
+ double e2 = j->get_energy();
+ if (e2 >= 10000)
+ break;
+ k = Lst_OBE[n3][t3].begin();
+ if (n2 == n3 && t2 == t3)
+ ++k;
+ if (e1 + e2 + k->get_energy() +
+ Lst_OBE[n4][t4].begin()->get_energy() + Gr >=
+ max_E)
+ break;
+ for (; k != Lst_OBE[n3][t3].end(); ++k) {
+ int p3 = k->get_pos1();
+ double e3 = k->get_energy();
+ if (e3 >= 10000)
+ break;
+ l = Lst_OBE[n4][t4].begin();
+ if (n3 == n4 && t3 == t4)
+ ++l;
+ if (e1 + e2 + e3 + l->get_energy() + Gr >= max_E)
+ break;
+ for (; l != Lst_OBE[n4][t4].end(); ++l) {
+ int p4 = l->get_pos1();
+ double e4 = l->get_energy();
+ if (e4 >= 10000)
+ break;
+ double e = e1 + e2 + e3 + e4 + Gr;
+ if (e >= max_E)
+ break;
+ if (!overlap(p1, p2, n1, n2) &&
+ !overlap(p1, p3, n1, n3) &&
+ !overlap(p1, p4, n1, n4) &&
+ !overlap(p2, p3, n2, n3) &&
+ !overlap(p2, p4, n2, n4) &&
+ !overlap(p3, p4, n3, n4)) {
+ if (min_E - e >= Epsilon) {
+ min_E = e;
+ max_E = e + theta;
+ }
+ count++;
- // V: Here, they calculate the probabilities of having a particular transition.
- for (int t=0; t<=Ns; t++) {
- for (int t1=0; t1<=Ns; t1++) {
- for (int t2=0; t2<=Ns; t2++) {
- for (int t3=0; t3<=Ns; t3++) {
- for (int t4=0; t4<=Ns; t4++) {
- if (t1==t || t2==t || t3==t || t4==t)
- Prob[t]+=exp_four[t1][t2][t3][t4];
- }
- }
- }
- }
- }
- // V: This was added for the competition of transitions:
- for (int ti=0; ti<=Ns; ti++) {
- for (int tj=0; tj<=Ns; tj++) {
- for (int t1=0; t1<=Ns; t1++) {
- for (int t2=0; t2<=Ns; t2++) {
- for (int t3=0; t3<=Ns; t3++) {
- for (int t4=0; t4<=Ns; t4++) {
-
- // Let's build vectors of type of transitions, then check if these transitions are included
- std::vector firstVec = { t1, t2, t3, t4};
- std::vector secondVec = { ti, tj};
-
- // Sort first vector
- std::sort(firstVec.begin(), firstVec.end());
- // Sort second vector
- std::sort(secondVec.begin(), secondVec.end());
-
- // Check if all elements of a second vector exists in first vector
- bool my_condition = std::includes(firstVec.begin(), firstVec.end(), secondVec.begin(), secondVec.end());
-
- if (my_condition) {
- Prob_compete[ti][tj] += exp_four[t1][t2][t3][t4];
+ if (write_profile) {
+ // double exponent = exp(-e/RT);
+ double exponent = exp(
+ -e / RT +
+ scaling_factor); // V: scaling factor added
+
+ exp_four[t1][t2][t3][t4] += exponent;
+ for (int t = 0; t <= Ns; t++)
+ runs[t] +=
+ exponent *
+ (delta_fnc(t, t1) + delta_fnc(t, t2) +
+ delta_fnc(t, t3) + delta_fnc(t, t4));
+ update_promatrix(p1, n1, t1, e, exponent);
+ update_promatrix(p2, n2, t2, e, exponent);
+ update_promatrix(p3, n3, t3, e, exponent);
+ update_promatrix(p4, n4, t4, e, exponent);
+ sum_4RB += exponent;
+ sum_4RG = sum_4RG + e * exponent;
}
+ }
}
+ }
}
+ }
}
+ }
}
+ }
}
+ }
}
+ }
+
+ StatesInFourRuns = (long)count;
+ ZsumB += sum_4RB;
+ ZsumG += sum_4RG;
+ totalFreq =
+ StatesInOneRun + StatesInTwoRuns + StatesInThreeRuns + StatesInFourRuns;
+ if (results) {
+ cout << "Number of four-run states = " << count << endl;
+ cout << "Total number of states = " << totalFreq << endl;
+
+ // cout << "Total Partition Function without exp= " << ZsumB << endl;
+ // V: The scaling factor needs to be removed from the partition function
+ cout << "Total Partition Function = " << exp(-scaling_factor) * ZsumB
+ << endl;
+ cout << "Scaling factor " << scaling_factor << endl;
- // V: Now, let's calculate the conditionals
- for (int t1=0; t1<=Ns; t1++) {
- for (int t2=0; t2<=Ns; t2++) {
- // V: The diagonal for Prob_conditional_not should be 0, right?
- Prob_conditional[t1][t2] = Prob_compete[t1][t2]/Prob[t2];
- // Probability of t1 given that we don't have t2
- double aa = Prob[t1]/ZsumB;
- double bb = Prob_compete[t1][t2]/ZsumB;
- double cc = Prob[t2]/ZsumB;
-
- //if (cc >= 1.0) Prob_conditional_not[t1][t2] = std::nan("NaN");
- // If the denominator is 0, then let's make the conditional probability 0.
- if (cc >= 1.0) Prob_conditional_not[t1][t2] = 0.0;
- else Prob_conditional_not[t1][t2] = (aa-bb)/(1.0-cc); // If it isn't 0, then calculate it.
+ cout << "Number of M_runs = " << runs[0] / ZsumB << endl;
+ cout << "Number of Z_runs = " << runs[1] / ZsumB << endl;
+ cout << "Number of C_runs = " << runs[2] / ZsumB << endl;
+ double runs = (sum_1RB + 2 * sum_2RB + 3 * sum_3RB + 4 * sum_4RB) / ZsumB;
+ cout << "Average number of runs = " << runs << endl;
+ }
+
+ // V: Here, they calculate the probabilities of having a particular
+ // transition.
+ for (int t = 0; t <= Ns; t++) {
+ for (int t1 = 0; t1 <= Ns; t1++) {
+ for (int t2 = 0; t2 <= Ns; t2++) {
+ for (int t3 = 0; t3 <= Ns; t3++) {
+ for (int t4 = 0; t4 <= Ns; t4++) {
+ if (t1 == t || t2 == t || t3 == t || t4 == t)
+ Prob[t] += exp_four[t1][t2][t3][t4];
+ }
}
+ }
}
+ }
+ // V: This was added for the competition of transitions:
+ for (int ti = 0; ti <= Ns; ti++) {
+ for (int tj = 0; tj <= Ns; tj++) {
+ for (int t1 = 0; t1 <= Ns; t1++) {
+ for (int t2 = 0; t2 <= Ns; t2++) {
+ for (int t3 = 0; t3 <= Ns; t3++) {
+ for (int t4 = 0; t4 <= Ns; t4++) {
+
+ // Let's build vectors of type of transitions, then check if these
+ // transitions are included
+ std::vector firstVec = {t1, t2, t3, t4};
+ std::vector secondVec = {ti, tj};
+
+ // Sort first vector
+ std::sort(firstVec.begin(), firstVec.end());
+ // Sort second vector
+ std::sort(secondVec.begin(), secondVec.end());
+
+ // Check if all elements of a second vector exists in first
+ // vector
+ bool my_condition =
+ std::includes(firstVec.begin(), firstVec.end(),
+ secondVec.begin(), secondVec.end());
- if (results) {
- cout << "Prob_M = " << Prob[0]/ZsumB << endl;
- cout<< "Prob_Z = " << Prob[1]/ZsumB << endl;
- cout << "Prob_C = " << Prob[2]/ZsumB << endl;
+ if (my_condition) {
+ Prob_compete[ti][tj] += exp_four[t1][t2][t3][t4];
+ }
+ }
+ }
+ }
+ }
}
+ }
- flag_minE4 = 0;
- return flag_minE4;
-}
+ // V: Now, let's calculate the conditionals
+ for (int t1 = 0; t1 <= Ns; t1++) {
+ for (int t2 = 0; t2 <= Ns; t2++) {
+ // V: The diagonal for Prob_conditional_not should be 0, right?
+ Prob_conditional[t1][t2] = Prob_compete[t1][t2] / Prob[t2];
+ // Probability of t1 given that we don't have t2
+ double aa = Prob[t1] / ZsumB;
+ double bb = Prob_compete[t1][t2] / ZsumB;
+ double cc = Prob[t2] / ZsumB;
+ // if (cc >= 1.0) Prob_conditional_not[t1][t2] = std::nan("NaN");
+ // If the denominator is 0, then let's make the conditional probability
+ // 0.
+ if (cc >= 1.0)
+ Prob_conditional_not[t1][t2] = 0.0;
+ else
+ Prob_conditional_not[t1][t2] =
+ (aa - bb) / (1.0 - cc); // If it isn't 0, then calculate it.
+ }
+ }
+
+ if (results) {
+ cout << "Prob_M = " << Prob[0] / ZsumB << endl;
+ cout << "Prob_Z = " << Prob[1] / ZsumB << endl;
+ cout << "Prob_C = " << Prob[2] / ZsumB << endl;
+ }
+
+ flag_minE4 = 0;
+ return flag_minE4;
+}
#endif
diff --git a/src/trans_compete/SIDD_4R.h b/src/trans_compete/SIDD_4R.h
index f7fd243..a4689e2 100644
--- a/src/trans_compete/SIDD_4R.h
+++ b/src/trans_compete/SIDD_4R.h
@@ -17,17 +17,15 @@
#include "SIDD_3R.h"
-class SIDD_4R : public SIDD_3R
-{
+class SIDD_4R : public SIDD_3R {
private:
- bool flag_minE4;
- long StatesInFourRuns;
-
+ bool flag_minE4;
+ long StatesInFourRuns;
+
public:
- SIDD_4R();
- virtual ~SIDD_4R();
- bool Search_Low4RE();
-
+ SIDD_4R();
+ virtual ~SIDD_4R();
+ bool Search_Low4RE();
};
#endif // !defined(AFX_SIDD_3R_H__B5CD0095_10DE_4C22_A817_256ABBAD9425__INCLUDED_)
diff --git a/src/trans_compete/SIDD_Base.cpp b/src/trans_compete/SIDD_Base.cpp
index 5d34422..a6a6c62 100644
--- a/src/trans_compete/SIDD_Base.cpp
+++ b/src/trans_compete/SIDD_Base.cpp
@@ -1,9 +1,9 @@
// SIDD_Base.cpp: implementation of the SIDD_Base class.
-// This class is the very base one defining all the parameters and all
+// This class is the very base one defining all the parameters and all
// the common data and function members
//
// program to implement algorithm developed by Craig Benham
-//
+//
// author: Chengpeng Bi
// modifiers: Dina Zhabinskaya, Sally Madden, Ian Korf
// compiler: g++
@@ -26,968 +26,892 @@ using namespace std;
// Construction/Destruction
//////////////////////////////////////////////////////////////////////
-SIDD_Base::SIDD_Base()
-{
- plasmid_seq = 0; // null pointer
- en_cruciforms = 0;
- for (int t=0; t<=Ns; t++) {
- profile[t] = 0;
- for(int m = 0; m <= MaxInitialWindowSize; m++) {
- promatrix[m][t] = 0;
- }
- }
-
- R = 8.314/(4.2*1000.0); // gas constant
- C = 3.6; // stiffness coefficient
- A = 10.4; // bps per turn in B-DNA
- Az = -12.0; //bps per turn in Z-DNA
- tz = 0.4; //undertwist at the two junctions of Z-DNA
- alpha = 0.0; // linking difference
- K0 = 2220.0; // linking energy coefficient
- pea = 0; // null pointer
- Flag_PEA = false;
- EnergyType = Copolymeric; // default value
- MoleculeType = Linear; // default value
-
- min_E = 0.0; // minimum energy
- max_E = 0.0;
- min_WS = 0; // the window size corresponding to the minimum energy
- write_profile = false; // flagging if writing starts
-
- ZsumB = 0.0; // summation of Boltzman frequence
- ZsumG = 0.0; // store partition function value
- sum_1RG = sum_2RG = sum_3RG = 0.0;
- sum_1RB = sum_2RB = sum_3RB = 0.0;
- // MEEE: I ADDED THIS FOR CALCULATING NUUMBER OF RUNS ACCURATELY
- sum_1RB = sum_2RB = sum_3RB = 0.0;
- set_showbase(1);
-}
-
-SIDD_Base::~SIDD_Base()
-{
-
- delete [] plasmid_seq;
- delete [] sa;
- delete [] pea;
- for(int i = 0; i < length_seq; i++) {
- delete [] en_cruciforms[i];
+SIDD_Base::SIDD_Base() {
+ plasmid_seq = 0; // null pointer
+ en_cruciforms = 0;
+ for (int t = 0; t <= Ns; t++) {
+ profile[t] = 0;
+ for (int m = 0; m <= MaxInitialWindowSize; m++) {
+ promatrix[m][t] = 0;
}
- for (int t=0; t<=Ns; t++) {
- delete [] profile[t];
- for(int i = MinWindowSize; i <= MaxWindowSize; i++) {
- delete [] promatrix[i][t];
- }
+ }
+
+ R = 8.314 / (4.2 * 1000.0); // gas constant
+ C = 3.6; // stiffness coefficient
+ A = 10.4; // bps per turn in B-DNA
+ Az = -12.0; // bps per turn in Z-DNA
+ tz = 0.4; // undertwist at the two junctions of Z-DNA
+ alpha = 0.0; // linking difference
+ K0 = 2220.0; // linking energy coefficient
+ pea = 0; // null pointer
+ Flag_PEA = false;
+ EnergyType = Copolymeric; // default value
+ MoleculeType = Linear; // default value
+
+ min_E = 0.0; // minimum energy
+ max_E = 0.0;
+ min_WS = 0; // the window size corresponding to the minimum energy
+ write_profile = false; // flagging if writing starts
+
+ ZsumB = 0.0; // summation of Boltzman frequence
+ ZsumG = 0.0; // store partition function value
+ sum_1RG = sum_2RG = sum_3RG = 0.0;
+ sum_1RB = sum_2RB = sum_3RB = 0.0;
+ // MEEE: I ADDED THIS FOR CALCULATING NUUMBER OF RUNS ACCURATELY
+ sum_1RB = sum_2RB = sum_3RB = 0.0;
+ set_showbase(1);
+}
+
+SIDD_Base::~SIDD_Base() {
+
+ delete[] plasmid_seq;
+ delete[] sa;
+ delete[] pea;
+ for (int i = 0; i < length_seq; i++) {
+ delete[] en_cruciforms[i];
+ }
+ for (int t = 0; t <= Ns; t++) {
+ delete[] profile[t];
+ for (int i = MinWindowSize; i <= MaxWindowSize; i++) {
+ delete[] promatrix[i][t];
}
+ }
}
-void SIDD_Base::set_stress_level(double stress_level)
-{
- if(stress_level > 0)stress_level = -stress_level;
- supdensity = stress_level;
+void SIDD_Base::set_stress_level(double stress_level) {
+ if (stress_level > 0)
+ stress_level = -stress_level;
+ supdensity = stress_level;
}
-void SIDD_Base::set_threshold(double threshold)
-{
- theta = threshold;
-}
+void SIDD_Base::set_threshold(double threshold) { theta = threshold; }
-//energy parameters
-void SIDD_Base::set_salt_conc(double salt_conc)
-{
- Salt_Conc = salt_conc;
- TMAT = 354.65 + 16.6*log10(Salt_Conc);
- TMGC = TMAT + 41.0;
+// energy parameters
+void SIDD_Base::set_salt_conc(double salt_conc) {
+ Salt_Conc = salt_conc;
+ TMAT = 354.65 + 16.6 * log10(Salt_Conc);
+ TMGC = TMAT + 41.0;
}
-//energy parameters
-void SIDD_Base::set_temperature(double temperature)
-{
- Temperature = temperature;
- RT = 1.9872*Temperature/1000.0; // constant
- BAT = 7.2464*(1.0 - Temperature / TMAT); // coefficient (kcal)
- BGC = 9.0172*(1.0 - Temperature / TMGC); // coefficient (kcal)
+// energy parameters
+void SIDD_Base::set_temperature(double temperature) {
+ Temperature = temperature;
+ RT = 1.9872 * Temperature / 1000.0; // constant
+ BAT = 7.2464 * (1.0 - Temperature / TMAT); // coefficient (kcal)
+ BGC = 9.0172 * (1.0 - Temperature / TMGC); // coefficient (kcal)
}
-void SIDD_Base::set_Cinitiation()
-{
- Ecr = 192.5-Temperature*0.565-4*BAT-2*2.44*RT*log(4);
+void SIDD_Base::set_Cinitiation() {
+ Ecr = 192.5 - Temperature * 0.565 - 4 * BAT - 2 * 2.44 * RT * log(4);
}
+void SIDD_Base::set_EnergyType(Energetics et) { EnergyType = et; }
-void SIDD_Base::set_EnergyType(Energetics et)
-{
- EnergyType = et;
+void SIDD_Base::set_min_e(double min_e) {
+ min_E = min_e;
+ max_E = min_E + theta;
}
+void SIDD_Base::set_MoleculeType(Molecule mt) { MoleculeType = mt; }
-void SIDD_Base::set_min_e(double min_e)
-{
- min_E = min_e;
- max_E = min_E + theta;
-}
-void SIDD_Base::set_MoleculeType(Molecule mt)
-{
- MoleculeType = mt;
+void SIDD_Base::set_MaxWindowSize() {
+ MaxWindowSize = 350;
+ if (length_seq < MaxWindowSize) {
+ MaxWindowSize = length_seq;
+ }
}
-void SIDD_Base::set_MaxWindowSize()
-{
- MaxWindowSize = 350;
- if (length_seq < MaxWindowSize){
- MaxWindowSize = length_seq;
- }
-}
-
-void SIDD_Base::set_MinWindowSize()
-{
- MinWindowSize = 1;
-}
+void SIDD_Base::set_MinWindowSize() { MinWindowSize = 1; }
// linking energy coefficient
-void SIDD_Base::set_K()
-{
- K = K0*RT / length_seq;
-}
-
-void SIDD_Base::set_cruciform_string(std::string cruciform_string)
-{
- cr_string = cruciform_string;
-}
+void SIDD_Base::set_K() { K = K0 * RT / length_seq; }
-//nucleation energies
-void SIDD_Base::set_junction()
-{
- a[0]=10.84; // Melt junction energy
- a[1]=10.0; // Z-DNA junction energy
- a[2]=0.0; //zero for cruciform since the energy is included in IR input string
+void SIDD_Base::set_cruciform_string(std::string cruciform_string) {
+ cr_string = cruciform_string;
}
-//linking difference
-void SIDD_Base::set_Alpha()
-{
- alpha = supdensity*length_seq / A;
-
+// nucleation energies
+void SIDD_Base::set_junction() {
+ a[0] = 10.84; // Melt junction energy
+ a[1] = 10.0; // Z-DNA junction energy
+ a[2] =
+ 0.0; // zero for cruciform since the energy is included in IR input string
}
-//sequence length
-int SIDD_Base::get_sequence_length(){
- if(MoleculeType == Linear)
- return (length_seq - Len_Gap_Seq);
- return length_seq;
-}
+// linking difference
+void SIDD_Base::set_Alpha() { alpha = supdensity * length_seq / A; }
-void SIDD_Base::set_sequence_length(int len){
- length_seq = len;
- if(MoleculeType != Circular)
- length_seq+=Len_Gap_Seq;
+// sequence length
+int SIDD_Base::get_sequence_length() {
+ if (MoleculeType == Linear)
+ return (length_seq - Len_Gap_Seq);
+ return length_seq;
}
-void SIDD_Base::set_showres(int showres)
-{
- results = showres;
+void SIDD_Base::set_sequence_length(int len) {
+ length_seq = len;
+ if (MoleculeType != Circular)
+ length_seq += Len_Gap_Seq;
}
-//read sequence
+void SIDD_Base::set_showres(int showres) { results = showres; }
+
+// read sequence
double SIDD_Base::prepare_sequence(std::string sequence) {
- plasmid_seq = new int[length_seq+1];
- for(int j = 0; j < length_seq; j++){
- plasmid_seq[j] = 0;
- }
- int read = 1;
- int i = 0;
- string::iterator iter;
- for(iter=sequence.begin();iter !=sequence.end();iter++)
- {
- char ch = *iter;
- if(ch == '>') {
- read = 0;
- }
- if(!read) {
- if(ch == '\n')
- read = 1;
- continue;
- }
- if(read) {
- ch = toupper(ch);
- switch(ch)
- {
- case 'A':
- plasmid_seq[i++] = 0;
- break;
- case 'C':
- plasmid_seq[i++] = 1;
- break;
- case 'G':
- plasmid_seq[i++] = 2;
- break;
- case 'T':
- plasmid_seq[i++] = 3;
- break;
- //turn Ns into Gs
- case 'N':
- plasmid_seq[i++] = 2;
- break;
- default:
- break;
- }
- }
+ plasmid_seq = new int[length_seq + 1];
+ for (int j = 0; j < length_seq; j++) {
+ plasmid_seq[j] = 0;
+ }
+ int read = 1;
+ int i = 0;
+ string::iterator iter;
+ for (iter = sequence.begin(); iter != sequence.end(); iter++) {
+ char ch = *iter;
+ if (ch == '>') {
+ read = 0;
}
- length_seq = i;
- if(!length_seq) {
- cerr << "sequence file not available.\n";
- return false;
+ if (!read) {
+ if (ch == '\n')
+ read = 1;
+ continue;
}
-
- if(MoleculeType != Circular)
- length_seq+=Len_Gap_Seq;
- return true;
+ if (read) {
+ ch = toupper(ch);
+ switch (ch) {
+ case 'A':
+ plasmid_seq[i++] = 0;
+ break;
+ case 'C':
+ plasmid_seq[i++] = 1;
+ break;
+ case 'G':
+ plasmid_seq[i++] = 2;
+ break;
+ case 'T':
+ plasmid_seq[i++] = 3;
+ break;
+ // turn Ns into Gs
+ case 'N':
+ plasmid_seq[i++] = 2;
+ break;
+ default:
+ break;
+ }
+ }
+ }
+ length_seq = i;
+ if (!length_seq) {
+ cerr << "sequence file not available.\n";
+ return false;
+ }
+
+ if (MoleculeType != Circular)
+ length_seq += Len_Gap_Seq;
+ return true;
}
-//read IR input string
+// read IR input string
void SIDD_Base::prepare_cruciforms() {
- string str1 = cr_string;
- size_t found1 = 0;
- int IR[2] = {0, 0};
- double energy = 0.0;
- int pos1 = 0;
- int pos2 = 0;
- string str2;
- string str3;
-
- en_cruciforms = new double* [length_seq+1];
- for(int i = 0; i < length_seq; i++){
- en_cruciforms[i] = new double[MaxWindowSize+1];
- }
-
- for(int i = 0; i < length_seq; i++){
- for(int j = 1; j <= MaxWindowSize; j++){
- en_cruciforms[i][j] = 10000;
- }
+ string str1 = cr_string;
+ size_t found1 = 0;
+ int IR[2] = {0, 0};
+ double energy = 0.0;
+ int pos1 = 0;
+ int pos2 = 0;
+ string str2;
+ string str3;
+
+ en_cruciforms = new double *[length_seq + 1];
+ for (int i = 0; i < length_seq; i++) {
+ en_cruciforms[i] = new double[MaxWindowSize + 1];
+ }
+
+ for (int i = 0; i < length_seq; i++) {
+ for (int j = 1; j <= MaxWindowSize; j++) {
+ en_cruciforms[i][j] = 10000;
}
-
- while (found1!=string::npos) {
- found1=str1.find("|",pos1);
- if(found1!=string::npos) {
- str2 = str1.substr(pos1,found1-pos1);
- pos2 = 0;
- size_t found2 = 0;
- int i = 0;
- while (found2!=string::npos) {
- found2=str2.find(",",pos2);
- if(found2!=string::npos) {
- str3 = str2.substr(pos2,found2-pos2);
- pos2 = int(found2)+1;
- if (i == 2)
- energy = atof(str3.c_str());
- else if (i < 2)
- IR[i] = atoi(str3.c_str());
- i++;
- }
- }
- if (i == 3 && IR[0] > 0 && IR[0] <= length_seq &&
- IR[1] > 0 && IR[1] <= MaxWindowSize)
- en_cruciforms[IR[0]-1][IR[1]] = energy;
- pos1 = int(found1)+1;
+ }
+
+ while (found1 != string::npos) {
+ found1 = str1.find("|", pos1);
+ if (found1 != string::npos) {
+ str2 = str1.substr(pos1, found1 - pos1);
+ pos2 = 0;
+ size_t found2 = 0;
+ int i = 0;
+ while (found2 != string::npos) {
+ found2 = str2.find(",", pos2);
+ if (found2 != string::npos) {
+ str3 = str2.substr(pos2, found2 - pos2);
+ pos2 = int(found2) + 1;
+ if (i == 2)
+ energy = atof(str3.c_str());
+ else if (i < 2)
+ IR[i] = atoi(str3.c_str());
+ i++;
}
+ }
+ if (i == 3 && IR[0] > 0 && IR[0] <= length_seq && IR[1] > 0 &&
+ IR[1] <= MaxWindowSize)
+ en_cruciforms[IR[0] - 1][IR[1]] = energy;
+ pos1 = int(found1) + 1;
}
- }
-
- //TODO: This is where the energetics are stored. But what does it actually mean each of the temrs...?
-//nearest neighbor energetics
-void SIDD_Base::set_Delta_G()
-{
- double Delta_H[4][4] = {
- {32.3, 35.8, 36.0, 31.2},
- {35.8, 35.1, 39.8, 36.0},
- {36.0, 39.8, 35.1, 35.8},
- {31.2, 36.0, 35.8, 32.3}
- };
-
- double Delta_S[4][4] ={ // J per mol per Kelvin degree
- {95.6, 100.4, 102.4, 95.3},
- {100.4, 93.9, 102.4, 102.4},
- {102.4, 102.4, 93.9, 100.4},
- {95.3, 102.4, 100.4, 95.6}
- };
-
-
- // unit conversion: KJ to Kcal
- // 1 kcal = 4.184 KJ
- int i, j;
- for(i = 0; i < 4; i++){
- for(j = 0; j < 4; j++){
- Delta_H[i][j] /= 4.184;
- Delta_S[i][j] /= 4184.0;
- }
- }
-
- for(i = 0; i < 4; i++){
- for(int j = 0; j < 4; j++){
- Delta_S[i][j] = 1.0/(16.6*log10(Salt_Conc/0.1)/Delta_H[i][j] + 1.0 / Delta_S[i][j]);
- }
- }
-
- for(i = 0; i < 4; i++){
- for(int j = 0; j < 4; j++){
- Delta_G[i][j] = Delta_H[i][j] - Temperature*Delta_S[i][j]; // kcal per mole
- }
- }
-}
-
-//superhelical energy
-double SIDD_Base::calc_Gres(int n_m, int n_z, int n_c, int nr)
-{
- // V: This has the form of 2013 paper competitive transitions DIna and Craig.
- double t1 = (alpha + n_m/A+ n_z/A - n_z/Az + 2*tz*nr + n_c/A) * (alpha + n_m/A+ n_z/A - n_z/Az + 2*tz*nr + n_c/A);
- double t2 = 2*PI*PI*C*K*t1 /(4*PI*PI*C + K*n_m);
- return t2;
-}
-
-void SIDD_Base::set_E_limit()
-{
- min_E = 0.5*K*alpha*alpha;
- max_E = min_E + theta;
-}
-
-//minimum residual energy
-void SIDD_Base::set_minGres()
-{
-minGres = calc_Gres(1,0,0,0);
- for (int r = 0; r <= 4; r++) {
- for(int w1 = 0; w1 <= MaxWindowSize; w1++){
- for(int w2 = 0; w2 <= MaxWindowSize; w2++){
- for(int w3 = 0; w3 <= MaxWindowSize; w3++){
- if(calc_Gres(w1,w2,w3,r) < minGres)
- minGres=calc_Gres(w1,w2,w3,r);
- }
- }
- }
+ }
+}
+
+// TODO: This is where the energetics are stored. But what does it actually mean
+// each of the temrs...?
+// nearest neighbor energetics
+void SIDD_Base::set_Delta_G() {
+ double Delta_H[4][4] = {{32.3, 35.8, 36.0, 31.2},
+ {35.8, 35.1, 39.8, 36.0},
+ {36.0, 39.8, 35.1, 35.8},
+ {31.2, 36.0, 35.8, 32.3}};
+
+ double Delta_S[4][4] = {// J per mol per Kelvin degree
+ {95.6, 100.4, 102.4, 95.3},
+ {100.4, 93.9, 102.4, 102.4},
+ {102.4, 102.4, 93.9, 100.4},
+ {95.3, 102.4, 100.4, 95.6}};
+
+ // unit conversion: KJ to Kcal
+ // 1 kcal = 4.184 KJ
+ int i, j;
+ for (i = 0; i < 4; i++) {
+ for (j = 0; j < 4; j++) {
+ Delta_H[i][j] /= 4.184;
+ Delta_S[i][j] /= 4184.0;
+ }
+ }
+
+ for (i = 0; i < 4; i++) {
+ for (int j = 0; j < 4; j++) {
+ Delta_S[i][j] = 1.0 / (16.6 * log10(Salt_Conc / 0.1) / Delta_H[i][j] +
+ 1.0 / Delta_S[i][j]);
+ }
+ }
+
+ for (i = 0; i < 4; i++) {
+ for (int j = 0; j < 4; j++) {
+ Delta_G[i][j] =
+ Delta_H[i][j] - Temperature * Delta_S[i][j]; // kcal per mole
}
+ }
}
-int SIDD_Base::delta_fnc(int x, int y) {
- return (x == y);
+// superhelical energy
+double SIDD_Base::calc_Gres(int n_m, int n_z, int n_c, int nr) {
+ // V: This has the form of 2013 paper competitive transitions DIna and Craig.
+ double t1 = (alpha + n_m / A + n_z / A - n_z / Az + 2 * tz * nr + n_c / A) *
+ (alpha + n_m / A + n_z / A - n_z / Az + 2 * tz * nr + n_c / A);
+ double t2 = 2 * PI * PI * C * K * t1 / (4 * PI * PI * C + K * n_m);
+ return t2;
}
-// computing transition energy: Z-DNA
-void SIDD_Base::calc_Z(int pos)
-{
-
- double energy_Z_AS[4][4] = {
- {3.9,4.6,3.4,5.9},
- {1.3,2.4,0.7,3.4},
- {3.4,4.0,2.4,4.6},
- {2.5,3.4,1.3,3.9}
- };
-
- double energy_Z_SA[4][4] = {
- {3.9,1.3,3.4,2.5},
- {4.6,2.4,4.0,3.4},
- {3.4,0.7,2.4,1.3},
- {5.9,3.4,4.6,3.9}
- };
-
- double energy_Z_ZZ[4][4] = {
- {7.4,4.5,6.3,5.6},
- {4.5,4.0,4.0,6.3},
- {6.3,4.0,4.0,4.5},
- {5.6,6.3,4.5,7.4}
- };
-
- // Base code before pos
- int p1 = plasmid_seq[pos-1 < 0 ? length_seq - 1 : pos - 1];
- // Base code after pos
- int p2 = plasmid_seq[pos+1 == length_seq ? 0 : pos + 1];
- // Base code of pos
- int p = plasmid_seq[pos];
- // Config code (anti or syn) before pos
- int c1 = sa[pos-1 < 0 ? length_seq - 1 : pos - 1];
- // Config code (anti or syn) after pos
- int c2 = sa[pos+1 == length_seq ? 0 : pos + 1];
- // Config code (anti or syn) of pos
- int c = sa[pos];
- e1=0;
- e2=0;
- zz=0;
- if (c != c2)
- {
- if (c == 'a') {
- e1 = energy_Z_AS[p][p2];
- }
- else {
- e1 = energy_Z_SA[p][p2];
- }
- }
-
- if (c != c1)
- {
- if (c == 'a') {
- e2 = energy_Z_AS[p][p1];
- }
- else {
- e2 = energy_Z_SA[p][p1];
- }
- }
- if (c == c1) {
- zz =energy_Z_ZZ[p][p1];
- int p0 = sa[pos-2 < 0 ? length_seq - 2 : pos - 2];
- if (energy_Z_ZZ[p][p2] >= energy_Z_ZZ[p0][p1]) {
- if (c=='a')
- e2 = energy_Z_SA[p1][p];
- else
- e2 = energy_Z_AS[p1][p];
- }
- else {
- if (c1=='a')
- e2 = energy_Z_AS[p1][p];
- else
- e2 = energy_Z_SA[p1][p];
- }
-
- }
- if (c == c2) {
- int p3 = plasmid_seq[pos+2 == length_seq ? 1 : pos + 2];
- if (energy_Z_ZZ[p2][p3] >= energy_Z_ZZ[p][p1]) {
- zz = energy_Z_ZZ[p][p1];
- if (c2=='a')
- e1 = energy_Z_SA[p][p2];
- else
- e1 = energy_Z_AS[p][p2];
- }
- else {
- zz = energy_Z_ZZ[p2][p3];
- if (c=='a')
- e1 = energy_Z_AS[p][p2];
- else
- e1 = energy_Z_SA[p][p2];
- }
- }
+void SIDD_Base::set_E_limit() {
+ min_E = 0.5 * K * alpha * alpha;
+ max_E = min_E + theta;
+}
+
+// minimum residual energy
+void SIDD_Base::set_minGres() {
+ minGres = calc_Gres(1, 0, 0, 0);
+ for (int r = 0; r <= 4; r++) {
+ for (int w1 = 0; w1 <= MaxWindowSize; w1++) {
+ for (int w2 = 0; w2 <= MaxWindowSize; w2++) {
+ for (int w3 = 0; w3 <= MaxWindowSize; w3++) {
+ if (calc_Gres(w1, w2, w3, r) < minGres)
+ minGres = calc_Gres(w1, w2, w3, r);
+ }
+ }
+ }
+ }
}
+int SIDD_Base::delta_fnc(int x, int y) { return (x == y); }
+
+// computing transition energy: Z-DNA
+void SIDD_Base::calc_Z(int pos) {
+
+ double energy_Z_AS[4][4] = {{3.9, 4.6, 3.4, 5.9},
+ {1.3, 2.4, 0.7, 3.4},
+ {3.4, 4.0, 2.4, 4.6},
+ {2.5, 3.4, 1.3, 3.9}};
+
+ double energy_Z_SA[4][4] = {{3.9, 1.3, 3.4, 2.5},
+ {4.6, 2.4, 4.0, 3.4},
+ {3.4, 0.7, 2.4, 1.3},
+ {5.9, 3.4, 4.6, 3.9}};
+
+ double energy_Z_ZZ[4][4] = {{7.4, 4.5, 6.3, 5.6},
+ {4.5, 4.0, 4.0, 6.3},
+ {6.3, 4.0, 4.0, 4.5},
+ {5.6, 6.3, 4.5, 7.4}};
+
+ // Base code before pos
+ int p1 = plasmid_seq[pos - 1 < 0 ? length_seq - 1 : pos - 1];
+ // Base code after pos
+ int p2 = plasmid_seq[pos + 1 == length_seq ? 0 : pos + 1];
+ // Base code of pos
+ int p = plasmid_seq[pos];
+ // Config code (anti or syn) before pos
+ int c1 = sa[pos - 1 < 0 ? length_seq - 1 : pos - 1];
+ // Config code (anti or syn) after pos
+ int c2 = sa[pos + 1 == length_seq ? 0 : pos + 1];
+ // Config code (anti or syn) of pos
+ int c = sa[pos];
+ e1 = 0;
+ e2 = 0;
+ zz = 0;
+ if (c != c2) {
+ if (c == 'a') {
+ e1 = energy_Z_AS[p][p2];
+ } else {
+ e1 = energy_Z_SA[p][p2];
+ }
+ }
+
+ if (c != c1) {
+ if (c == 'a') {
+ e2 = energy_Z_AS[p][p1];
+ } else {
+ e2 = energy_Z_SA[p][p1];
+ }
+ }
+ if (c == c1) {
+ zz = energy_Z_ZZ[p][p1];
+ int p0 = sa[pos - 2 < 0 ? length_seq - 2 : pos - 2];
+ if (energy_Z_ZZ[p][p2] >= energy_Z_ZZ[p0][p1]) {
+ if (c == 'a')
+ e2 = energy_Z_SA[p1][p];
+ else
+ e2 = energy_Z_AS[p1][p];
+ } else {
+ if (c1 == 'a')
+ e2 = energy_Z_AS[p1][p];
+ else
+ e2 = energy_Z_SA[p1][p];
+ }
+ }
+ if (c == c2) {
+ int p3 = plasmid_seq[pos + 2 == length_seq ? 1 : pos + 2];
+ if (energy_Z_ZZ[p2][p3] >= energy_Z_ZZ[p][p1]) {
+ zz = energy_Z_ZZ[p][p1];
+ if (c2 == 'a')
+ e1 = energy_Z_SA[p][p2];
+ else
+ e1 = energy_Z_AS[p][p2];
+ } else {
+ zz = energy_Z_ZZ[p2][p3];
+ if (c == 'a')
+ e1 = energy_Z_AS[p][p2];
+ else
+ e1 = energy_Z_SA[p][p2];
+ }
+ }
+}
// start position - startp
// window size - n
// return - sum of window energy
-double SIDD_Base::sum_WindowZ(int startp, int n)
-{
- double s = 0.0;
- if (n%2==1 || n==2 || n==4 || n==6) {
- s=10000;
- return s;
- }
- else {
- if(startp + n - 1 < length_seq){
- for(int i = startp; i < startp + n; i++) {
- calc_Z(i);
- if (i==startp)
- s+=e1/2;
- if ((i-startp)%2==0 && i!=startp)
- s+=e1/2+zz;
- if ((i-startp)%2==1)
- s+=e2/2;
- }
- return s;
- }
- else{
- for(int i = startp; i < length_seq; i++) {
- calc_Z(i);
- if (i==startp)
- s+=e1/2;
- if ((i-startp)%2==0 && i!=startp)
- s+=e1/2+zz;
- if ((i-startp)%2==1)
- s+=e2/2;
- }
- for(int j = 0; j < startp + n - length_seq; j++) {
- calc_Z(j);
- if (j==startp)
- s+=e1/2;
- if ((length_seq-startp+j)%2==0 && j!=startp)
- s+=e1/2+zz;
- if ((length_seq-startp+j)%2==1)
- s+=e2/2;
- }
- return s;
- }
- }
+double SIDD_Base::sum_WindowZ(int startp, int n) {
+ double s = 0.0;
+ if (n % 2 == 1 || n == 2 || n == 4 || n == 6) {
+ s = 10000;
+ return s;
+ } else {
+ if (startp + n - 1 < length_seq) {
+ for (int i = startp; i < startp + n; i++) {
+ calc_Z(i);
+ if (i == startp)
+ s += e1 / 2;
+ if ((i - startp) % 2 == 0 && i != startp)
+ s += e1 / 2 + zz;
+ if ((i - startp) % 2 == 1)
+ s += e2 / 2;
+ }
+ return s;
+ } else {
+ for (int i = startp; i < length_seq; i++) {
+ calc_Z(i);
+ if (i == startp)
+ s += e1 / 2;
+ if ((i - startp) % 2 == 0 && i != startp)
+ s += e1 / 2 + zz;
+ if ((i - startp) % 2 == 1)
+ s += e2 / 2;
+ }
+ for (int j = 0; j < startp + n - length_seq; j++) {
+ calc_Z(j);
+ if (j == startp)
+ s += e1 / 2;
+ if ((length_seq - startp + j) % 2 == 0 && j != startp)
+ s += e1 / 2 + zz;
+ if ((length_seq - startp + j) % 2 == 1)
+ s += e2 / 2;
+ }
+ return s;
+ }
+ }
}
-
// computing transition energy: neighbor interation
-double SIDD_Base::calc_NI(int pos)
-{
- if(pos < length_seq && pos >= 0){
-
- if(Flag_PEA && pea[pos] != 0.0)
- return pea[pos]; // assign specified energy
- int p1 = plasmid_seq[pos-1 < 0 ? length_seq - 1 : pos - 1];
- int p2 = plasmid_seq[pos+1 == length_seq ? 0 : pos + 1];
- int p = plasmid_seq[pos];
- return (Delta_G[p1][p] + Delta_G[p][p2]) / 2.0;
- }
- else{
- cerr << "position: out of range.\n";
- return 0;
- }
+double SIDD_Base::calc_NI(int pos) {
+ if (pos < length_seq && pos >= 0) {
+
+ if (Flag_PEA && pea[pos] != 0.0)
+ return pea[pos]; // assign specified energy
+ int p1 = plasmid_seq[pos - 1 < 0 ? length_seq - 1 : pos - 1];
+ int p2 = plasmid_seq[pos + 1 == length_seq ? 0 : pos + 1];
+ int p = plasmid_seq[pos];
+ return (Delta_G[p1][p] + Delta_G[p][p2]) / 2.0;
+ } else {
+ cerr << "position: out of range.\n";
+ return 0;
+ }
}
// start position - startp
// window size - n
// return - sum of window energy
-double SIDD_Base::sum_WindowNI(int startp, int n)
-{
- double s = 0.0;
- if(startp + n - 1 < length_seq){
- for(int i = startp; i < startp + n; i++)
- s += calc_NI(i);
- return s;
- }
- else{
- for(int i = startp; i < length_seq; i++)
- s += calc_NI(i);
- for(int j = 0; j < startp + n - length_seq; j++)
- s += calc_NI(j);
- }
- return s;
-}
-
-double SIDD_Base::sum_WindowNIbyBases(int startp, int n)
-{
- int at, gc;
- at = gc = 0;
- double s = 0.0;
- double tempE = 0.0;
- count_AT_GC(startp, n, at, gc, tempE);
- s = BAT*at + BGC*gc + tempE;
+double SIDD_Base::sum_WindowNI(int startp, int n) {
+ double s = 0.0;
+ if (startp + n - 1 < length_seq) {
+ for (int i = startp; i < startp + n; i++)
+ s += calc_NI(i);
return s;
+ } else {
+ for (int i = startp; i < length_seq; i++)
+ s += calc_NI(i);
+ for (int j = 0; j < startp + n - length_seq; j++)
+ s += calc_NI(j);
+ }
+ return s;
}
-double SIDD_Base::sum_cruciform(int startp, int n)
-{
- return en_cruciforms[startp][n];
+double SIDD_Base::sum_WindowNIbyBases(int startp, int n) {
+ int at, gc;
+ at = gc = 0;
+ double s = 0.0;
+ double tempE = 0.0;
+ count_AT_GC(startp, n, at, gc, tempE);
+ s = BAT * at + BGC * gc + tempE;
+ return s;
}
-// counting A/T bases
-void SIDD_Base::count_AT_GC(int startp, int n, int& c_AT, int& c_GC, double& tempE)
-{
- c_AT = c_GC = 0;
- tempE = 0.0;
- if(startp + n - 1 < length_seq){
- for(int i = startp; i < startp + n; i++){
- if(Flag_PEA && pea[i] != 0.0)
- tempE += pea[i];
- else{
- if(plasmid_seq[i] == 0 || plasmid_seq[i] == 3)
- c_AT++;
- else
- c_GC++;
- }
- }
- }
-
- else{
- for(int i = startp; i < length_seq; i++){
- if(Flag_PEA && pea[i] != 0.0)
- tempE += pea[i];
- else{
-
- if(plasmid_seq[i] == 0 || plasmid_seq[i] == 3)
- c_AT++;
- else
- c_GC++;
- }
- }
-
- for(int j = 0; j < startp + n - length_seq; j++){
- if(Flag_PEA && pea[j] != 0.0)
- tempE += pea[j];
- else{
- if(plasmid_seq[j] == 0 || plasmid_seq[j] == 3)
- c_AT++;
- else
- c_GC++;
- }
- }
- }
+double SIDD_Base::sum_cruciform(int startp, int n) {
+ return en_cruciforms[startp][n];
}
+// counting A/T bases
+void SIDD_Base::count_AT_GC(int startp, int n, int &c_AT, int &c_GC,
+ double &tempE) {
+ c_AT = c_GC = 0;
+ tempE = 0.0;
+ if (startp + n - 1 < length_seq) {
+ for (int i = startp; i < startp + n; i++) {
+ if (Flag_PEA && pea[i] != 0.0)
+ tempE += pea[i];
+ else {
+ if (plasmid_seq[i] == 0 || plasmid_seq[i] == 3)
+ c_AT++;
+ else
+ c_GC++;
+ }
+ }
+ }
+
+ else {
+ for (int i = startp; i < length_seq; i++) {
+ if (Flag_PEA && pea[i] != 0.0)
+ tempE += pea[i];
+ else {
+
+ if (plasmid_seq[i] == 0 || plasmid_seq[i] == 3)
+ c_AT++;
+ else
+ c_GC++;
+ }
+ }
+
+ for (int j = 0; j < startp + n - length_seq; j++) {
+ if (Flag_PEA && pea[j] != 0.0)
+ tempE += pea[j];
+ else {
+ if (plasmid_seq[j] == 0 || plasmid_seq[j] == 3)
+ c_AT++;
+ else
+ c_GC++;
+ }
+ }
+ }
+}
// energy of a transition for a sequence segment
// starting at starp of length n for transition defined by type
-double SIDD_Base::calc_OPenBasesEnergy(int startp, int n, int type)
-{
- if (type == 0) {
- switch(EnergyType){
- case Near_Neighbor:
- return sum_WindowNI(startp, n);
- case Copolymeric:
- return sum_WindowNIbyBases(startp, n);
- default:
- return sum_WindowNIbyBases(startp, n);
- }
+double SIDD_Base::calc_OPenBasesEnergy(int startp, int n, int type) {
+ if (type == 0) {
+ switch (EnergyType) {
+ case Near_Neighbor:
+ return sum_WindowNI(startp, n);
+ case Copolymeric:
+ return sum_WindowNIbyBases(startp, n);
+ default:
+ return sum_WindowNIbyBases(startp, n);
}
- else if (type == 1)
- return sum_WindowZ(startp,n);
- else
- return sum_cruciform(startp,n);
+ } else if (type == 1)
+ return sum_WindowZ(startp, n);
+ else
+ return sum_cruciform(startp, n);
}
-
// reading sequence and initializing parameters
-bool SIDD_Base::initializer(std::string sequence,int len)
-{
- set_sequence_length(len);
- if(!prepare_sequence(sequence)) {
- //sequence either 0 or too long
- return 0;
- }
- set_MaxWindowSize();
- set_MinWindowSize();
- set_K();
- set_Alpha();
- set_junction();
- set_Delta_G();
- set_E_limit();
- prepare_cruciforms();
-
- sa = new char[length_seq+1];
- pea = new double[length_seq+1];
- for(int i = 0; i < length_seq+1; i++)
- pea[i] = 0.0;
-
- for (int t=0; t<=Ns; t++) {
- profile[t] = new G_x[length_seq+1];
- for(int m = MinWindowSize; m <= MaxWindowSize; m++) {
- promatrix[m][t] = new G_x[length_seq+1];
- }
- }
-
- if(MoleculeType == Linear){
- for(int k = length_seq - Len_Gap_Seq; k < length_seq; k++)
- plasmid_seq[k] = 2; //add a tail of G's
- }
-
- set_minGres();
- reset_profile();
-
- for(int i = 0; i < length_seq; i++){
- switch(plasmid_seq[i]){
- case 0:
- sa[i]='s';
- break;
- case 1:
- sa[i]='a';
- break;
- case 2:
- sa[i]='s';
- break;
- case 3:
- sa[i]='a';
- break;
- }
- }
- for(int i = 0; i < length_seq; i++){
- // Config code (anti or syn) before pos
- int c1 = sa[i-1 < 0 ? length_seq - 1 : i - 1];
- // Config code (anti or syn) after pos
- int c2 = sa[i+1 == length_seq ? 0 : i + 1];
- // Config code (anti or syn) of pos
- int c = sa[i];
- if (c==c1 && c==c2) {
- if (c=='s')
- sa[i]='a';
- else
- sa[i]='s';
- }
- }
-
- for(int i = 0; i < length_seq; i++){
- // Config code (anti or syn) after pos
- int c2 = sa[i+1 == length_seq ? 0 : i + 1];
- // Config code (anti or syn) of pos
- int c = sa[i];
- int c3 = sa[i+2 == length_seq ? 1 : i + 2];
- int c4 = sa[i+3 == length_seq ? 2 : i + 3];
- int c5 = sa[i+4 == length_seq ? 3 : i + 4];
- if (c==c2 && c3==c4 && c!=c3) {
- if (c5=='a') {
- sa[i]='a';
- sa[i+1]='s';
- sa[i+2]='a';
- sa[i+3]='s';
- }
- else {
- sa[i]='s';
- sa[i+1]='a';
- sa[i+2]='s';
- sa[i+3]='a';
- }
- }
- }
-
- return true;
+bool SIDD_Base::initializer(std::string sequence, int len) {
+ set_sequence_length(len);
+ if (!prepare_sequence(sequence)) {
+ // sequence either 0 or too long
+ return 0;
+ }
+ set_MaxWindowSize();
+ set_MinWindowSize();
+ set_K();
+ set_Alpha();
+ set_junction();
+ set_Delta_G();
+ set_E_limit();
+ prepare_cruciforms();
+
+ sa = new char[length_seq + 1];
+ pea = new double[length_seq + 1];
+ for (int i = 0; i < length_seq + 1; i++)
+ pea[i] = 0.0;
+
+ for (int t = 0; t <= Ns; t++) {
+ profile[t] = new G_x[length_seq + 1];
+ for (int m = MinWindowSize; m <= MaxWindowSize; m++) {
+ promatrix[m][t] = new G_x[length_seq + 1];
+ }
+ }
+
+ if (MoleculeType == Linear) {
+ for (int k = length_seq - Len_Gap_Seq; k < length_seq; k++)
+ plasmid_seq[k] = 2; // add a tail of G's
+ }
+
+ set_minGres();
+ reset_profile();
+
+ for (int i = 0; i < length_seq; i++) {
+ switch (plasmid_seq[i]) {
+ case 0:
+ sa[i] = 's';
+ break;
+ case 1:
+ sa[i] = 'a';
+ break;
+ case 2:
+ sa[i] = 's';
+ break;
+ case 3:
+ sa[i] = 'a';
+ break;
+ }
+ }
+ for (int i = 0; i < length_seq; i++) {
+ // Config code (anti or syn) before pos
+ int c1 = sa[i - 1 < 0 ? length_seq - 1 : i - 1];
+ // Config code (anti or syn) after pos
+ int c2 = sa[i + 1 == length_seq ? 0 : i + 1];
+ // Config code (anti or syn) of pos
+ int c = sa[i];
+ if (c == c1 && c == c2) {
+ if (c == 's')
+ sa[i] = 'a';
+ else
+ sa[i] = 's';
+ }
+ }
+
+ for (int i = 0; i < length_seq; i++) {
+ // Config code (anti or syn) after pos
+ int c2 = sa[i + 1 == length_seq ? 0 : i + 1];
+ // Config code (anti or syn) of pos
+ int c = sa[i];
+ int c3 = sa[i + 2 == length_seq ? 1 : i + 2];
+ int c4 = sa[i + 3 == length_seq ? 2 : i + 3];
+ int c5 = sa[i + 4 == length_seq ? 3 : i + 4];
+ if (c == c2 && c3 == c4 && c != c3) {
+ if (c5 == 'a') {
+ sa[i] = 'a';
+ sa[i + 1] = 's';
+ sa[i + 2] = 'a';
+ sa[i + 3] = 's';
+ } else {
+ sa[i] = 's';
+ sa[i + 1] = 'a';
+ sa[i + 2] = 's';
+ sa[i + 3] = 'a';
+ }
+ }
+ }
+
+ return true;
}
// assign energy position-wide
-void SIDD_Base::assign_pea(char* efile)
-{
- ifstream ins(efile);
-
- if(ins.is_open()){
- while(!ins.eof()){
- double e0;
- int p;
- ins >> p >> e0;
- if(p >=0 && p < length_seq && e0 != 0.0){
- pea[p] = e0;
- cout << "assign energy at " << p << ": " << pea[p] << endl;
- }
- }
- }
- else{
- cerr << "fail to open the file: " << efile << endl;
- Flag_PEA = false;
- }
- ins.close();
-}
-
-void SIDD_Base::reset_promatrix()
-{
- ZsumB = 0.0;
- ZsumG = 0.0;
- totalFreq = 0;
- for (int t=0; t<=Ns; t++) {
- for(int k = MinWindowSize; k <= MaxWindowSize; k++){
- for(int i = 0; i < length_seq; i++){
- promatrix[k][t][i].reset();
- }
- }
- }
+void SIDD_Base::assign_pea(char *efile) {
+ ifstream ins(efile);
+
+ if (ins.is_open()) {
+ while (!ins.eof()) {
+ double e0;
+ int p;
+ ins >> p >> e0;
+ if (p >= 0 && p < length_seq && e0 != 0.0) {
+ pea[p] = e0;
+ cout << "assign energy at " << p << ": " << pea[p] << endl;
+ }
+ }
+ } else {
+ cerr << "fail to open the file: " << efile << endl;
+ Flag_PEA = false;
+ }
+ ins.close();
+}
+
+void SIDD_Base::reset_promatrix() {
+ ZsumB = 0.0;
+ ZsumG = 0.0;
+ totalFreq = 0;
+ for (int t = 0; t <= Ns; t++) {
+ for (int k = MinWindowSize; k <= MaxWindowSize; k++) {
+ for (int i = 0; i < length_seq; i++) {
+ promatrix[k][t][i].reset();
+ }
+ }
+ }
}
// startp - start position
// n - window size
// x - free energy
-bool SIDD_Base::update_promatrix(int startp, int n, int struc, double x, double bzfactor)
-{
- if(startp < 0 || startp >= length_seq || n > MaxWindowSize)
- return false;
-
- promatrix[n][struc][startp].add(x, bzfactor, RT);
- return true;
-}
-
-void SIDD_Base::fill_profile()
-{
- for(int t=0; t<=Ns; t++) {
- for(int n = MinWindowSize; n <= MaxWindowSize; n++){
- for(int startp = 0; startp < length_seq; startp++){
- double lastG = promatrix[n][t][startp].get_sum_xG();
- double lastB = promatrix[n][t][startp].get_sum_xB();
- if(startp + n - 1 < length_seq){
- for(int i = startp; i < startp + n; i++){
- profile[t][i].add(-1.0, 0.0, -1.0, lastG, lastB);
- }
- }
- else{
- for(int i = startp; i < length_seq; i++)
- profile[t][i].add(-1.0, 0.0, -1.0, lastG, lastB);
- for(int k = 0; k < startp + n - length_seq; k++)
- profile[t][k].add(-1.0, 0.0, -1.0 , lastG, lastB);
- }
- }
- }
- }
-}
-
-void SIDD_Base::reset_profile()
-{
- ZsumB = 0.0;
- ZsumG = 0.0;
- totalFreq = 0;
- for(int t=0; t<=Ns; t++) {
- for(int i = 0; i < length_seq; i++){
- profile[t][i].reset();
- }
- }
-
- reset_promatrix();
+bool SIDD_Base::update_promatrix(int startp, int n, int struc, double x,
+ double bzfactor) {
+ if (startp < 0 || startp >= length_seq || n > MaxWindowSize)
+ return false;
+
+ promatrix[n][struc][startp].add(x, bzfactor, RT);
+ return true;
+}
+
+void SIDD_Base::fill_profile() {
+ for (int t = 0; t <= Ns; t++) {
+ for (int n = MinWindowSize; n <= MaxWindowSize; n++) {
+ for (int startp = 0; startp < length_seq; startp++) {
+ double lastG = promatrix[n][t][startp].get_sum_xG();
+ double lastB = promatrix[n][t][startp].get_sum_xB();
+ if (startp + n - 1 < length_seq) {
+ for (int i = startp; i < startp + n; i++) {
+ profile[t][i].add(-1.0, 0.0, -1.0, lastG, lastB);
+ }
+ } else {
+ for (int i = startp; i < length_seq; i++)
+ profile[t][i].add(-1.0, 0.0, -1.0, lastG, lastB);
+ for (int k = 0; k < startp + n - length_seq; k++)
+ profile[t][k].add(-1.0, 0.0, -1.0, lastG, lastB);
+ }
+ }
+ }
+ }
}
+void SIDD_Base::reset_profile() {
+ ZsumB = 0.0;
+ ZsumG = 0.0;
+ totalFreq = 0;
+ for (int t = 0; t <= Ns; t++) {
+ for (int i = 0; i < length_seq; i++) {
+ profile[t][i].reset();
+ }
+ }
+
+ reset_promatrix();
+}
// startp - start position
// n - window size
// x - free energy
-void SIDD_Base::calc_profile()
-{
- // V: avg_Gs My guess is that it is the average energy.
- double ave_Gs = 0.0;
- if(ZsumB == 0.0)
- ave_Gs = 0.0;
- else
- ave_Gs = ZsumG / ZsumB;
- // V: Goes through all the bases for each transition, and calculates the free energy required for transition
- // and the probability of transitioning.
- // profile is of class G_x, so check those functions (calc for example).
- for(int t=0; t<=Ns; t++) {
- for(int i = 0; i < length_seq; i++){
- profile[t][i].calc(ave_Gs, ZsumB);
- }
- }
- // V: The following lines is just a search of the maximum Gx for each transition, to substitute the Infinite.
- maxGx = profile[1][0].get_ave_Gx();
- for(int t=0; t<=Ns; t++) {
- for(int j = 1; j < length_seq; j++){
- if(maxGx < profile[t][j].get_ave_Gx()){
- maxGx = profile[t][j].get_ave_Gx(); // find max Gx
- }
- }
- }
- // V: This is for the substitution.
- for(int t=0; t<=Ns; t++) {
- for(int k = 0; k < length_seq; k++){
- if(profile[t][k].get_ave_Gx() == INFINITE_Gx)
- profile[t][k].set_ave_Gx(maxGx); // replace INF Gx with maxGx
- }
- }
-}
-
-//average number of bp in each transition
-void SIDD_Base::sum_open()
-{
- for(int t=0; t<=Ns; t++) {
- sumn[t]=0.0;
- for(int i = 0; i < length_seq; i++) {
- sumn[t]+=profile[t][i].get_px();
- }
- }
- if (results) {
- cout << "N_Melting = " << sumn[0] << endl;
- cout << "N_Z = " << sumn[1] << endl;
- cout << "N_Cruciform = " << sumn[2] << endl;
- }
+void SIDD_Base::calc_profile() {
+ // V: avg_Gs My guess is that it is the average energy.
+ double ave_Gs = 0.0;
+ if (ZsumB == 0.0)
+ ave_Gs = 0.0;
+ else
+ ave_Gs = ZsumG / ZsumB;
+ // V: Goes through all the bases for each transition, and calculates the free
+ // energy required for transition
+ // and the probability of transitioning.
+ // profile is of class G_x, so check those functions (calc for example).
+ for (int t = 0; t <= Ns; t++) {
+ for (int i = 0; i < length_seq; i++) {
+ profile[t][i].calc(ave_Gs, ZsumB);
+ }
+ }
+ // V: The following lines is just a search of the maximum Gx for each
+ // transition, to substitute the Infinite.
+ maxGx = profile[1][0].get_ave_Gx();
+ for (int t = 0; t <= Ns; t++) {
+ for (int j = 1; j < length_seq; j++) {
+ if (maxGx < profile[t][j].get_ave_Gx()) {
+ maxGx = profile[t][j].get_ave_Gx(); // find max Gx
+ }
+ }
+ }
+ // V: This is for the substitution.
+ for (int t = 0; t <= Ns; t++) {
+ for (int k = 0; k < length_seq; k++) {
+ if (profile[t][k].get_ave_Gx() == INFINITE_Gx)
+ profile[t][k].set_ave_Gx(maxGx); // replace INF Gx with maxGx
+ }
+ }
+}
+
+// average number of bp in each transition
+void SIDD_Base::sum_open() {
+ for (int t = 0; t <= Ns; t++) {
+ sumn[t] = 0.0;
+ for (int i = 0; i < length_seq; i++) {
+ sumn[t] += profile[t][i].get_px();
+ }
+ }
+ if (results) {
+ cout << "N_Melting = " << sumn[0] << endl;
+ cout << "N_Z = " << sumn[1] << endl;
+ cout << "N_Cruciform = " << sumn[2] << endl;
+ }
}
void SIDD_Base::get_column_header() {
-if(get_showbase()) {
- cout << "Position" << "\tBase" << "\tP_melt" << "\tP_Z"<<"\tP_cruciform"<
#include
-#include
-#include
#include
+#include
+#include
#include
+#include
#include
#include "G_x.h"
@@ -32,155 +32,156 @@ const int MaxInitialWindowSize = 350;
const double GMAX_ADJ = 10.22;
const int Len_Gap_Seq = 50;
const int Ns = 2;
-// V: This is where the scaling factor is declared. It might have to be modified for higher superhelical values.
+// V: This is where the scaling factor is declared. It might have to be modified
+// for higher superhelical values.
const double scaling_factor = 300.0;
-enum Energetics{Copolymeric, Near_Neighbor, Z_DNA};
-enum Molecule{Circular, Linear};
-enum Transition{on, off};
+enum Energetics { Copolymeric, Near_Neighbor, Z_DNA };
+enum Molecule { Circular, Linear };
+enum Transition { on, off };
-class SIDD_Base
-{
+class SIDD_Base {
protected:
- int* plasmid_seq; // hold encoded sequence
- double** en_cruciforms; //array of cruciform pos,length,energies
- int length_seq;
- double Delta_G[4][4]; // free energy
- double Temperature;
- double Salt_Conc; // salt concentration
- double R; // constant
- double C; // torsional stiffness
- double a[3]; // initial energy for melting, Z-DNA
-
- double A; // bases per turn
- double Az;
- double tz;
- double stresslevel;
- double supdensity;
- double alpha; // linking difference
- double theta; // threshold
- double K0;
- double K;
- double RT;
- double Ecr;
- double TMAT;
- double TMGC;
- double BAT;
- double BGC;
- double Zmin;
- int MaxWindowSize;
- int MinWindowSize;
-
- std::string cr_string;
-
- double min_E;
- double max_E;
- int min_WS;
-
- char* sa;
- double e1;
- double e2;
- double zz;
-
- Energetics EnergyType;
- Molecule MoleculeType;
- Transition TransitionState;
-
- // profile information
- G_x* profile[3]; // hold profile
- G_x* promatrix[MaxInitialWindowSize+1][3];
-
- double minGres;
- double sumn[3];
-
-
- double* pea; // position-wide energy assignment
- bool Flag_PEA;
- long totalFreq;
- bool write_profile;
-
- double ZsumG; // total free energy
- double ZsumB; // value by partition function Z
- double sum_1RG; // summation of all states within one run
- double sum_1RB; // summation of all Boltzman factors within one run
- double sum_2RG; // summation of all states within two runs
- double sum_2RB; // summation of all Boltzman factors within two runs
- double sum_3RG;
- double sum_3RB; // summation of all Boltzman factors within three runs
- double sum_4RG;
- double sum_4RB; // summation of all Boltzman factors within three runs
-
- double Prob_1R;
- double Prob_2R;
- double Prob_3R;
- double maxGx;
- double minGx;
- int showbase;
- int results;
+ int *plasmid_seq; // hold encoded sequence
+ double **en_cruciforms; // array of cruciform pos,length,energies
+ int length_seq;
+ double Delta_G[4][4]; // free energy
+ double Temperature;
+ double Salt_Conc; // salt concentration
+ double R; // constant
+ double C; // torsional stiffness
+ double a[3]; // initial energy for melting, Z-DNA
+
+ double A; // bases per turn
+ double Az;
+ double tz;
+ double stresslevel;
+ double supdensity;
+ double alpha; // linking difference
+ double theta; // threshold
+ double K0;
+ double K;
+ double RT;
+ double Ecr;
+ double TMAT;
+ double TMGC;
+ double BAT;
+ double BGC;
+ double Zmin;
+ int MaxWindowSize;
+ int MinWindowSize;
+
+ std::string cr_string;
+
+ double min_E;
+ double max_E;
+ int min_WS;
+
+ char *sa;
+ double e1;
+ double e2;
+ double zz;
+
+ Energetics EnergyType;
+ Molecule MoleculeType;
+ Transition TransitionState;
+
+ // profile information
+ G_x *profile[3]; // hold profile
+ G_x *promatrix[MaxInitialWindowSize + 1][3];
+
+ double minGres;
+ double sumn[3];
+
+ double *pea; // position-wide energy assignment
+ bool Flag_PEA;
+ long totalFreq;
+ bool write_profile;
+
+ double ZsumG; // total free energy
+ double ZsumB; // value by partition function Z
+ double sum_1RG; // summation of all states within one run
+ double sum_1RB; // summation of all Boltzman factors within one run
+ double sum_2RG; // summation of all states within two runs
+ double sum_2RB; // summation of all Boltzman factors within two runs
+ double sum_3RG;
+ double sum_3RB; // summation of all Boltzman factors within three runs
+ double sum_4RG;
+ double sum_4RB; // summation of all Boltzman factors within three runs
+
+ double Prob_1R;
+ double Prob_2R;
+ double Prob_3R;
+ double maxGx;
+ double minGx;
+ int showbase;
+ int results;
private:
- void set_K();
- void set_Delta_G();
- void set_Alpha();
- void set_junction();
- void set_E_limit();
- void set_minGres();
- void set_MinWindowSize();
- //sally
- void set_MaxWindowSize();
+ void set_K();
+ void set_Delta_G();
+ void set_Alpha();
+ void set_junction();
+ void set_E_limit();
+ void set_minGres();
+ void set_MinWindowSize();
+ // sally
+ void set_MaxWindowSize();
public:
- SIDD_Base();
- virtual ~SIDD_Base();
- bool initializer(std::string filename,int len);
- void set_EnergyType(Energetics);
- void set_temperature(double temperature);
- void set_showres(int showres);
- void set_threshold(double threshold);
- void set_Cinitiation();
- void set_cruciform_string(std::string cruciform_string);
- void set_MoleculeType(Molecule);
- void show_seq();
- void show_deltaG();
- void show_parameter();
- double get_Delta_G(int i, int j){return (i< 4 && j < 4)?Delta_G[i][j]:0.0;};
- double calc_Gres(int,int,int,int);
- double calc_NI(int);
- void calc_Z(int);
- double sum_WindowNI(int startp, int n);
- double sum_WindowZ(int startp, int n);
- void count_AT_GC(int startp, int n, int& c_AT, int& c_GC, double& tempE);
- double sum_WindowNIbyBases(int startp, int n);
- double sum_cruciform(int startp, int n);
- double calc_OPenBasesEnergy(int startp, int n, int type);
- bool update_promatrix(int startp, int n, int struc, double x, double bzfactor);
- void reset_promatrix();
- void get_column_header();
- void fill_profile();
- void reset_profile();
- void sum_open();
- void show_profile();
- void calc_profile();
- int delta_fnc(int,int);
- void write_close(){write_profile = false;}
- void write_open(){write_profile = true;}
- void set_Flag_PEA(bool f){Flag_PEA = f;}
- bool get_Flag_PEA(){return Flag_PEA;}
- void assign_pea(char*);
- double prepare_sequence(std::string dnasequence);
- void prepare_cruciforms();
-
- void set_salt_conc(double salt_conc);
- void set_min_e(double min_e);
- void set_stress_level(double stress_level);
- double get_zsumb(){return ZsumB;}
-
- int get_sequence_length();
- void set_sequence_length(int len);
- char decode_base(int base);
- void set_showbase(int show) {showbase = show;}
- int get_showbase() {return showbase;}
-
+ SIDD_Base();
+ virtual ~SIDD_Base();
+ bool initializer(std::string filename, int len);
+ void set_EnergyType(Energetics);
+ void set_temperature(double temperature);
+ void set_showres(int showres);
+ void set_threshold(double threshold);
+ void set_Cinitiation();
+ void set_cruciform_string(std::string cruciform_string);
+ void set_MoleculeType(Molecule);
+ void show_seq();
+ void show_deltaG();
+ void show_parameter();
+ double get_Delta_G(int i, int j) {
+ return (i < 4 && j < 4) ? Delta_G[i][j] : 0.0;
+ };
+ double calc_Gres(int, int, int, int);
+ double calc_NI(int);
+ void calc_Z(int);
+ double sum_WindowNI(int startp, int n);
+ double sum_WindowZ(int startp, int n);
+ void count_AT_GC(int startp, int n, int &c_AT, int &c_GC, double &tempE);
+ double sum_WindowNIbyBases(int startp, int n);
+ double sum_cruciform(int startp, int n);
+ double calc_OPenBasesEnergy(int startp, int n, int type);
+ bool update_promatrix(int startp, int n, int struc, double x,
+ double bzfactor);
+ void reset_promatrix();
+ void get_column_header();
+ void fill_profile();
+ void reset_profile();
+ void sum_open();
+ void show_profile();
+ void calc_profile();
+ int delta_fnc(int, int);
+ void write_close() { write_profile = false; }
+ void write_open() { write_profile = true; }
+ void set_Flag_PEA(bool f) { Flag_PEA = f; }
+ bool get_Flag_PEA() { return Flag_PEA; }
+ void assign_pea(char *);
+ double prepare_sequence(std::string dnasequence);
+ void prepare_cruciforms();
+
+ void set_salt_conc(double salt_conc);
+ void set_min_e(double min_e);
+ void set_stress_level(double stress_level);
+ double get_zsumb() { return ZsumB; }
+
+ int get_sequence_length();
+ void set_sequence_length(int len);
+ char decode_base(int base);
+ void set_showbase(int show) { showbase = show; }
+ int get_showbase() { return showbase; }
};
#endif // !defined(AFX_SIDD_BASE_H__FA2336AB_FC74_40DB_AD46_8B6DE9CDFA25__INCLUDED_)
diff --git a/src/trans_compete/qsidd.cpp b/src/trans_compete/qsidd.cpp
index 00b5c66..0b5fc13 100644
--- a/src/trans_compete/qsidd.cpp
+++ b/src/trans_compete/qsidd.cpp
@@ -1,6 +1,6 @@
// qsidd.cpp: main() function
// program to implement algorithm developed by Craig Benham
-//
+//
// author: Chengpeng Bi
// modifiers: Dina Zhabinskaya, Sally Madden, Ian Korf
// compiler: g++
@@ -12,17 +12,16 @@
//////////////////////////////////////////////////////////////////////////////
#include "SIDD_4R.h"
-#include
#include
-#include
#include
-#include
#include
+#include
+#include
+#include
// #include
// #include
-
static char help[] = "\
QSIDD Help\n\
A program for computing the competition of melting, Z-DNA, and cruciform transitions in superhelical DNA.\n\n\
@@ -70,154 +69,190 @@ options:\n\
-r print ensemble average results\n\
";
+int main(int argc, char *argv[]) {
+ SIDD_4R sidd;
+ time_t time_1, time_2;
+ int c;
+ extern int optind;
+
+ // default parameters
+ Energetics et = Copolymeric;
+ Molecule mt = Linear;
+ char *energy = NULL;
+ char *dnafile = NULL;
-int main(int argc, char* argv[])
-{
- SIDD_4R sidd;
- time_t time_1, time_2;
- int c;
- extern int optind;
-
- // default parameters
- Energetics et = Copolymeric;
- Molecule mt = Linear;
- char *energy = NULL;
- char *dnafile = NULL;
-
- stringstream dnasequence;
- stringstream cruciform_string;
- int showbase = 0;
- int verbose = 0;
- int showpar = 0;
- int showres = 0;
- double salt_conc = 0.01;
- double temperature = 310.00;
- double stress_level=0.06;
- double threshold = 12;
- int usefile = 0;
-
- // option processing
- while ((c = getopt(argc, argv, "hcfnbsvpriTXte:")) != -1) {
- switch (c) {
- case 'c': mt = Circular; break;
- case 'n': et = Near_Neighbor; break;
- case 'f': usefile = 1; break;
- case 'e': energy = argv[optind++]; break;
- case 'b': showbase = 1; break;
- case 'v': verbose = 1; break;
- case 'p': showpar = 1; break;
- case 'r': showres = 1; break;
- case 'T':temperature=atof(argv[optind++]);break;
- case 's':stress_level=atof(argv[optind++]);break;
- case 'i': salt_conc=atof(argv[optind++]);break;
- case 't': threshold=atof(argv[optind++]);break;
- case 'X': cruciform_string << argv[optind++]; break;
- case 'h': cout< 15)
- cout << "WARNING: threshold is too high, execution time may be very long"<< endl;
- if (stress_level > 0.15 || stress_level < -0.15)
- cout << "WARNING: superhelical density is outside of physiological range"<< endl;
- if (temperature < 220 || temperature > 320)
- cout << "WARNING: temperature is outside of physiological range"<< endl;
- if (salt_conc < 0.0001)
- cout << "WARNING: salt concentration is outside of physiological range"<< endl;
- if (length < 1500)
- cout << "WARNING: sequence length is too short"<< endl;
- if (length > 10000)
- cout << "WARNING: sequence length is too long"<< endl;
-
- //sidd.show_seq(); // show sequence for debugging
-
- // V: In the code, they use the word Open very often, but it doesn't always refer to melting, it is more like the
- // calculation of the energy at a given state (melt, Z, cruci...). It seems that it was recycled from the
- // strand-separation code SIDD.
- // V: gen_OpenBaseEnergy() actually calculates the energies of states with 1 run for each type of transition.
- // I think the objective is to obtain the 2-dimensional vector/list Lst_OBE[window_size][transition].
- // Each entry of this vector, is a list containing objects of type stat_1R with associated position and energy.
-
- sidd.gen_OpenBaseEnergy(); // SIDD_1R (calculate opening energies and sort w/ increasing energy for each window size 1 to 250)
- sidd.write_close(); //don't write profile yet, just find minE
- sidd.Search_Low1RE(); // this step is here to update minE (from the zero-run state) if it's found in one-run states
- sidd.write_open(); // start to store info
- if(verbose)
- cout << "writing profile...\n";
- if (showpar)
- sidd.show_parameter();
- sidd.reset_profile(); // initialize profile
- sidd.Search_Low1RE(); // searching for one-run states and now storing info
- sidd.Search_Low2RE(); // searching states for two-run
- sidd.Search_Low3RE(); // searching states for three-run
- sidd.Search_Low4RE(); // searching states for four-run
-
- // V: I think for all these Search_LowNRE, global quantities such as the number of states, partition function,
- // energies Boltzmann factors and everything is calculated (the method is applied).
- // V: But the following lines actually calculate parameters for each region on the DNA
- // (e.g., the probability of transitioning each bp).
-
- // V: I am not sure what fill_profile() does, but I think it prepares the array with the profiles...
- sidd.fill_profile(); // computing profile
- // V: calc_profile does what it says, it calculates the profiles...
- sidd.calc_profile(); // calculate profile
- // V: This one sums probabilites, which basically gives you the number of runs in a particular transition.
- sidd.sum_open();
- time_2 = time(NULL); //end time
+ }
+ // start time
+ time_1 = time(NULL);
+ sidd.set_showres(showres);
+ sidd.set_salt_conc(salt_conc);
+ sidd.set_stress_level(stress_level);
+ sidd.set_temperature(temperature);
+ sidd.set_MoleculeType(mt);
+ sidd.set_EnergyType(et);
+ sidd.set_threshold(threshold);
+ sidd.set_cruciform_string(cruciform_string.str().c_str());
+ sidd.set_showbase(showbase);
+
+ if (argc - optind != 1) {
+ cerr << usage << endl;
+ exit(1);
+ }
+
+ if (usefile) {
+ dnafile = argv[argc - 1];
if (showpar)
- cout << "Run time = " << (time_2-time_1) << " sec" << endl;
- // V: This is where they print everything.
- sidd.show_profile(); // send output to screen or a disk file
-
- return 0; // end of program
-
+ cout << "DNA sequence file: " << dnafile << endl;
+ ifstream f(dnafile);
+ if (f) {
+ dnasequence << f.rdbuf();
+ f.close();
+ }
+ } else {
+ dnasequence << argv[argc - 1];
+ }
+ int length = dnasequence.str().size();
+
+ if (verbose) {
+ if (energy) {
+ cout << "Energy assignment file: " << energy << endl;
+ }
+ sidd.set_Flag_PEA(energy);
+ }
+
+ if (!sidd.initializer(dnasequence.str().c_str(), length)) {
+ cerr << "sidd.initializer failed" << endl;
+ exit(1);
+ }
+ if (sidd.get_Flag_PEA())
+ sidd.assign_pea(energy);
+
+ if (threshold < 9)
+ cout << "WARNING: threshold is too small, results may be inaccurate"
+ << endl;
+ if (threshold > 15)
+ cout << "WARNING: threshold is too high, execution time may be very long"
+ << endl;
+ if (stress_level > 0.15 || stress_level < -0.15)
+ cout << "WARNING: superhelical density is outside of physiological range"
+ << endl;
+ if (temperature < 220 || temperature > 320)
+ cout << "WARNING: temperature is outside of physiological range" << endl;
+ if (salt_conc < 0.0001)
+ cout << "WARNING: salt concentration is outside of physiological range"
+ << endl;
+ if (length < 1500)
+ cout << "WARNING: sequence length is too short" << endl;
+ if (length > 10000)
+ cout << "WARNING: sequence length is too long" << endl;
+
+ // sidd.show_seq(); // show sequence for debugging
+
+ // V: In the code, they use the word Open very often, but it doesn't always
+ // refer to melting, it is more like the
+ // calculation of the energy at a given state (melt, Z, cruci...). It seems
+ // that it was recycled from the strand-separation code SIDD.
+ // V: gen_OpenBaseEnergy() actually calculates the energies of states with 1
+ // run for each type of transition.
+ // I think the objective is to obtain the 2-dimensional vector/list
+ // Lst_OBE[window_size][transition]. Each entry of this vector, is a list
+ // containing objects of type stat_1R with associated position and energy.
+
+ sidd.gen_OpenBaseEnergy(); // SIDD_1R (calculate opening energies and sort w/
+ // increasing energy for each window size 1 to 250)
+ sidd.write_close(); // don't write profile yet, just find minE
+ sidd.Search_Low1RE(); // this step is here to update minE (from the zero-run
+ // state) if it's found in one-run states
+ sidd.write_open(); // start to store info
+ if (verbose)
+ cout << "writing profile...\n";
+ if (showpar)
+ sidd.show_parameter();
+ sidd.reset_profile(); // initialize profile
+ sidd.Search_Low1RE(); // searching for one-run states and now storing info
+ sidd.Search_Low2RE(); // searching states for two-run
+ sidd.Search_Low3RE(); // searching states for three-run
+ sidd.Search_Low4RE(); // searching states for four-run
+
+ // V: I think for all these Search_LowNRE, global quantities such as the
+ // number of states, partition function, energies Boltzmann factors and
+ // everything is calculated (the method is applied). V: But the following
+ // lines actually calculate parameters for each region on the DNA (e.g., the
+ // probability of transitioning each bp).
+
+ // V: I am not sure what fill_profile() does, but I think it prepares the
+ // array with the profiles...
+ sidd.fill_profile(); // computing profile
+ // V: calc_profile does what it says, it calculates the profiles...
+ sidd.calc_profile(); // calculate profile
+ // V: This one sums probabilites, which basically gives you the number of runs
+ // in a particular transition.
+ sidd.sum_open();
+ time_2 = time(NULL); // end time
+ if (showpar)
+ cout << "Run time = " << (time_2 - time_1) << " sec" << endl;
+ // V: This is where they print everything.
+ sidd.show_profile(); // send output to screen or a disk file
+
+ return 0; // end of program
}
diff --git a/src/trans_compete/stat_1R.cpp b/src/trans_compete/stat_1R.cpp
index 87520f3..c8f29cd 100644
--- a/src/trans_compete/stat_1R.cpp
+++ b/src/trans_compete/stat_1R.cpp
@@ -21,38 +21,28 @@
// Construction/Destruction
//////////////////////////////////////////////////////////////////////
-stat_1R::stat_1R(int s, double e, double RT) //constructor
+stat_1R::stat_1R(int s, double e, double RT) // constructor
{
- start_pos1 = s; //starting position
- energy = e;
+ start_pos1 = s; // starting position
+ energy = e;
}
-stat_1R::~stat_1R() //destructor ???
-{
-
-}
+stat_1R::~stat_1R() // destructor ???
+{}
+bool stat_1R::operator<(const stat_1R &x1) const {
-
-bool stat_1R::operator<(const stat_1R& x1)const
-{
-
- return (energy < x1.get_energy());
-
+ return (energy < x1.get_energy());
}
-bool stat_1R::operator==(const stat_1R& x1)const
-{
-
- return (energy == x1.get_energy());
+bool stat_1R::operator==(const stat_1R &x1) const {
+ return (energy == x1.get_energy());
}
-bool stat_1R::operator>(const stat_1R& x1)const
-{
-
- return (energy > x1.get_energy());
+bool stat_1R::operator>(const stat_1R &x1) const {
+ return (energy > x1.get_energy());
}
#endif
diff --git a/src/trans_compete/stat_1R.h b/src/trans_compete/stat_1R.h
index 952b2d1..27ee50a 100644
--- a/src/trans_compete/stat_1R.h
+++ b/src/trans_compete/stat_1R.h
@@ -16,19 +16,19 @@
#ifndef STAT_1R_H_
#define STAT_1R_H_
-class stat_1R
-{
+class stat_1R {
protected:
- int start_pos1;
- double energy;
+ int start_pos1;
+ double energy;
+
public:
- stat_1R(int start_position, double energy, double rt);
- virtual ~stat_1R();
- int get_pos1() const{return start_pos1;}
- bool operator <(const stat_1R&) const;
- bool operator ==(const stat_1R&)const;
- bool operator >(const stat_1R&) const;
- double get_energy()const {return energy;}
+ stat_1R(int start_position, double energy, double rt);
+ virtual ~stat_1R();
+ int get_pos1() const { return start_pos1; }
+ bool operator<(const stat_1R &) const;
+ bool operator==(const stat_1R &) const;
+ bool operator>(const stat_1R &) const;
+ double get_energy() const { return energy; }
};
#endif // !defined(AFX_STAT_1R_H__ACD4635F_79AA_463A_A66C_6B1CA4476EE0__INCLUDED_)
diff --git a/src/trans_three/G_x.cpp b/src/trans_three/G_x.cpp
index 373f7fa..8c805ab 100644
--- a/src/trans_three/G_x.cpp
+++ b/src/trans_three/G_x.cpp
@@ -2,7 +2,7 @@
// This class is designed for storing profiles e.g. p(x), G(x), A(x)
//
// program to implement algorithm developed by Craig Benham
-//
+//
// author: Chengpeng Bi
// modifiers: Dina Zhabinskaya, Sally Madden, Ian Korf
// compiler: g++
@@ -13,8 +13,6 @@
// UC Davis Genome Center
//////////////////////////////////////////////////////////////////////////////
-
-
#ifndef G_x_CPP
#define G_x_CPP
@@ -24,53 +22,47 @@
// Construction/Destruction
//////////////////////////////////////////////////////////////////////
-G_x::G_x()
-{
- sum_xG = 0.0;
- sum_xB = 0.0;
- ave_Gx = 0.0;
- px = 0.0;
-}
-
-void G_x::reset()
-{
- sum_xG = 0.0;
- sum_xB = 0.0;
- ave_Gx = 0.0;
- px = 0.0;
+G_x::G_x() {
+ sum_xG = 0.0;
+ sum_xB = 0.0;
+ ave_Gx = 0.0;
+ px = 0.0;
}
-
-void G_x::add(double gs, double exponent, double rt, double lastG, double lastB)
-{
- if(gs == -1.0 && rt == -1.0){
- sum_xG = sum_xG + lastG;
- sum_xB = sum_xB + lastB;
- }
- else{
- sum_xG = sum_xG + gs*exponent + lastG;
- sum_xB = sum_xB + exponent + lastB;
- }
+void G_x::reset() {
+ sum_xG = 0.0;
+ sum_xB = 0.0;
+ ave_Gx = 0.0;
+ px = 0.0;
}
-
-// HERE IS THE DIVISION THAT MAKES IT P=1.0, when ZsumB == 0.0. We need another way to figure out this.
-void G_x::calc(double ave_Gs, double ZsumB) //input: average energy and partition function
-{
- if (ZsumB == 0.0)
- px = 1.0;
- else
- px = sum_xB / ZsumB;
-
- if(sum_xB == 0.0)
- ave_Gx = INFINITE_Gx;
- else
- ave_Gx = sum_xG / sum_xB - ave_Gs;
+void G_x::add(double gs, double exponent, double rt, double lastG,
+ double lastB) {
+ if (gs == -1.0 && rt == -1.0) {
+ sum_xG = sum_xG + lastG;
+ sum_xB = sum_xB + lastB;
+ } else {
+ sum_xG = sum_xG + gs * exponent + lastG;
+ sum_xB = sum_xB + exponent + lastB;
+ }
}
-G_x::~G_x()
+// HERE IS THE DIVISION THAT MAKES IT P=1.0, when ZsumB == 0.0. We need another
+// way to figure out this.
+void G_x::calc(double ave_Gs,
+ double ZsumB) // input: average energy and partition function
{
+ if (ZsumB == 0.0)
+ px = 1.0;
+ else
+ px = sum_xB / ZsumB;
+ if (sum_xB == 0.0)
+ ave_Gx = INFINITE_Gx;
+ else
+ ave_Gx = sum_xG / sum_xB - ave_Gs;
}
+G_x::~G_x() {}
+
#endif
diff --git a/src/trans_three/G_x.h b/src/trans_three/G_x.h
index 4132658..d940d28 100644
--- a/src/trans_three/G_x.h
+++ b/src/trans_three/G_x.h
@@ -2,7 +2,7 @@
// This class is designed for storing profiles e.g. p(x), G(x), A(x)
//
// program to implement algorithm developed by Craig Benham
-//
+//
// author: Chengpeng Bi
// modifiers: Dina Zhabinskaya, Sally Madden, Ian Korf
// compiler: g++
@@ -19,26 +19,25 @@
#include
const double INFINITE_Gx = -10000.0;
-class G_x
-{
+class G_x {
private:
- double sum_xG; // total energy with states of x opening
- double sum_xB; // sum of Boltzman's factors with states of x open
- double ave_Gx; // average free energy of opening base x
- double px; // probability of opening base x
+ double sum_xG; // total energy with states of x opening
+ double sum_xB; // sum of Boltzman's factors with states of x open
+ double ave_Gx; // average free energy of opening base x
+ double px; // probability of opening base x
public:
- G_x();
- void reset();
- void add(double gs, double exponent, double rt, double lastG = 0.0, double lastB = 0.0); // add one state
- void calc(double ave_Gs, double zsumb); // computing p(x) and G(x)
- double get_sum_xG(){return sum_xG;}
- double get_ave_Gx(){return ave_Gx;}
- void set_ave_Gx(double x){ave_Gx = x;}
- double get_sum_xB(){return sum_xB;}
- double get_px(){return px;}
- virtual ~G_x();
-
+ G_x();
+ void reset();
+ void add(double gs, double exponent, double rt, double lastG = 0.0,
+ double lastB = 0.0); // add one state
+ void calc(double ave_Gs, double zsumb); // computing p(x) and G(x)
+ double get_sum_xG() { return sum_xG; }
+ double get_ave_Gx() { return ave_Gx; }
+ void set_ave_Gx(double x) { ave_Gx = x; }
+ double get_sum_xB() { return sum_xB; }
+ double get_px() { return px; }
+ virtual ~G_x();
};
-#endif
+#endif
diff --git a/src/trans_three/Makefile b/src/trans_three/Makefile
index 1908f1d..98fa52a 100644
--- a/src/trans_three/Makefile
+++ b/src/trans_three/Makefile
@@ -36,7 +36,7 @@ depend: $(OBJECTS:.o=.cpp)
$(CXX) -MM $^ > $@
test:
- ./qsidd example.fasta
+ ./qsidd example.fasta
tar:
rm -rf /tmp/$(APP)
@@ -70,4 +70,3 @@ coverage:
################
include depend
-
diff --git a/src/trans_three/SIDD_1R.cpp b/src/trans_three/SIDD_1R.cpp
index 3c3f954..046af1b 100644
--- a/src/trans_three/SIDD_1R.cpp
+++ b/src/trans_three/SIDD_1R.cpp
@@ -2,7 +2,7 @@
// This class is designed for one run derived from SIDD_Base
//
// program to implement algorithm developed by Craig Benham
-//
+//
// author: Chengpeng Bi
// modifiers: Dina Zhabinskaya, Sally Madden, Ian Korf
// compiler: g++
@@ -22,77 +22,78 @@
// Construction/Destruction
//////////////////////////////////////////////////////////////////////
-SIDD_1R::SIDD_1R():SIDD_Base() {
- min_RE = min_E; //minimun energy from the continuous part
- flag_minE1 = false;
- StatesInOneRun = 0;
-}
-SIDD_1R::~SIDD_1R()
-{
-
+SIDD_1R::SIDD_1R() : SIDD_Base() {
+ min_RE = min_E; // minimun energy from the continuous part
+ flag_minE1 = false;
+ StatesInOneRun = 0;
}
+SIDD_1R::~SIDD_1R() {}
-void SIDD_1R::gen_OpenBaseEnergy()
-{
- for(int i = MinWindowSize; i <= MaxWindowSize; i++){ // window size: 1 - MaxWindowSize
- for(int j = 0; j < length_seq; j++){ // start position to length_seq - 1
- double e = a + calc_OPenBasesEnergy(j, i);
- if(e < min_RE) min_RE = e;
- stat_1R s1r(j, e, RT); //input j=starting position, e=energy, RT=constant
- Lst_OBE[i].push_back(s1r);
- }
- }
- // sorting each list of OBE
- for(int k = MinWindowSize; k <= MaxWindowSize; k++){
- Lst_OBE[k].sort();
- }
+void SIDD_1R::gen_OpenBaseEnergy() {
+ for (int i = MinWindowSize; i <= MaxWindowSize;
+ i++) { // window size: 1 - MaxWindowSize
+ for (int j = 0; j < length_seq; j++) { // start position to length_seq - 1
+ double e = a + calc_OPenBasesEnergy(j, i);
+ if (e < min_RE)
+ min_RE = e;
+ stat_1R s1r(j, e, RT); // input j=starting position, e=energy, RT=constant
+ Lst_OBE[i].push_back(s1r);
+ }
+ }
+ // sorting each list of OBE
+ for (int k = MinWindowSize; k <= MaxWindowSize; k++) {
+ Lst_OBE[k].sort();
+ }
}
-bool SIDD_1R::Search_Low1RE()
-{
- long count = 0;
- flag_minE1 = false;
- sum_1RG = sum_1RB = 0.0;
- LstState::iterator j;
- //collecting one run states withing the energy threshold
- for(int i = MinWindowSize; i <= MaxWindowSize; i++){
- for(j = Lst_OBE[i].begin(); j != Lst_OBE[i].end(); ++j){
- double e = j->get_energy();
- if (e >= 10000)
- break;
- if(e +minGres >= max_E) break; //outside of threshold
- e += Gres[i][0];
- if(e < max_E) { //below threshold
- count++; //counting the number of one run states
- if(e < min_E){
- min_E = e; //update min_E
- max_E = e + theta;
- min_WS = i; // the window size corresponding to the minimum energy
- flag_minE1 = true;
- }
- if(write_profile){
- // double exponent= exp(-e/RT);
- double exponent= exp(-e/RT + scaling_factor); // V: I added this
+bool SIDD_1R::Search_Low1RE() {
+ long count = 0;
+ flag_minE1 = false;
+ sum_1RG = sum_1RB = 0.0;
+ LstState::iterator j;
+ // collecting one run states withing the energy threshold
+ for (int i = MinWindowSize; i <= MaxWindowSize; i++) {
+ for (j = Lst_OBE[i].begin(); j != Lst_OBE[i].end(); ++j) {
+ double e = j->get_energy();
+ if (e >= 10000)
+ break;
+ if (e + minGres >= max_E)
+ break; // outside of threshold
+ e += Gres[i][0];
+ if (e < max_E) { // below threshold
+ count++; // counting the number of one run states
+ if (e < min_E) {
+ min_E = e; // update min_E
+ max_E = e + theta;
+ min_WS = i; // the window size corresponding to the minimum energy
+ flag_minE1 = true;
+ }
+ if (write_profile) {
+ // double exponent= exp(-e/RT);
+ double exponent = exp(-e / RT + scaling_factor); // V: I added this
- // V: I'll include here the residual superhelical density
- update_promatrix(j->get_pos1(), i, e, exponent);
- // update_promatrix(j->get_pos1(), i, e, exponent);
- sum_1RB += exponent;
- sum_1RG = sum_1RG + e*exponent;
- }
- }
- }
- }
- StatesInOneRun = count;
- // adding the zero open bases term to the partition function
- // double closed = exp(-alpha*alpha*K/2/RT);
- double closed = exp(-alpha*alpha*K/2/RT + scaling_factor); // V: Note that the zero runs state is added to the partition function, but not to the probabilities.
+ // V: I'll include here the residual superhelical density
+ update_promatrix(j->get_pos1(), i, e, exponent);
+ // update_promatrix(j->get_pos1(), i, e, exponent);
+ sum_1RB += exponent;
+ sum_1RG = sum_1RG + e * exponent;
+ }
+ }
+ }
+ }
+ StatesInOneRun = count;
+ // adding the zero open bases term to the partition function
+ // double closed = exp(-alpha*alpha*K/2/RT);
+ double closed =
+ exp(-alpha * alpha * K / 2 / RT +
+ scaling_factor); // V: Note that the zero runs state is added to the
+ // partition function, but not to the probabilities.
- ZsumB = sum_1RB + closed;
- ZsumG = sum_1RG+alpha*alpha*K/2*closed;
- if (write_profile && results)
- cout << "Number of one-run states = " << count << endl;
- return flag_minE1;
+ ZsumB = sum_1RB + closed;
+ ZsumG = sum_1RG + alpha * alpha * K / 2 * closed;
+ if (write_profile && results)
+ cout << "Number of one-run states = " << count << endl;
+ return flag_minE1;
}
#endif
diff --git a/src/trans_three/SIDD_1R.h b/src/trans_three/SIDD_1R.h
index f5d0c78..0530d16 100644
--- a/src/trans_three/SIDD_1R.h
+++ b/src/trans_three/SIDD_1R.h
@@ -2,7 +2,7 @@
// This class is designed for one run derived from SIDD_Base
//
// program to implement algorithm developed by Craig Benham
-//
+//
// author: Chengpeng Bi
// modifiers: Dina Zhabinskaya, Sally Madden, Ian Korf
// compiler: g++
@@ -13,7 +13,6 @@
// UC Davis Genome Center
//////////////////////////////////////////////////////////////////////////////
-
#ifndef SIDD_1R_H_
#define SIDD_1R_H_
@@ -22,24 +21,23 @@
using namespace std;
-typedef list LstState;
+typedef list LstState;
-class SIDD_1R : public SIDD_Base
-{
+class SIDD_1R : public SIDD_Base {
protected:
- LstState Lst_OBE[MaxInitialWindowSize+1]; // lists of open base energy
+ LstState Lst_OBE[MaxInitialWindowSize + 1]; // lists of open base energy
- bool flag_minE1; // flagging if new min_E found in one run
- double min_RE; // store minimum run energy (a + NI)
- long StatesInOneRun;
+ bool flag_minE1; // flagging if new min_E found in one run
+ double min_RE; // store minimum run energy (a + NI)
+ long StatesInOneRun;
public:
- SIDD_1R();
- virtual ~SIDD_1R();
- void gen_OpenBaseEnergy();
- bool Search_Low1RE();
- void Update_Low1RE();
- void Show_Low1RE();
+ SIDD_1R();
+ virtual ~SIDD_1R();
+ void gen_OpenBaseEnergy();
+ bool Search_Low1RE();
+ void Update_Low1RE();
+ void Show_Low1RE();
};
#endif // !defined(AFX_SIDD_1R_H__2688B961_9B45_44D5_8F6F_F45636E04460__INCLUDED_)
diff --git a/src/trans_three/SIDD_2R.cpp b/src/trans_three/SIDD_2R.cpp
index 1e82c1a..492774c 100644
--- a/src/trans_three/SIDD_2R.cpp
+++ b/src/trans_three/SIDD_2R.cpp
@@ -2,7 +2,7 @@
// This class is designed for two runs derived from SIDD_1R
//
// program to implement algorithm developed by Craig Benham
-//
+//
// author: Chengpeng Bi
// modifiers: Dina Zhabinskaya, Sally Madden, Ian Korf
// compiler: g++
@@ -22,115 +22,112 @@
// Construction/Destruction
//////////////////////////////////////////////////////////////////////
-SIDD_2R::SIDD_2R():SIDD_1R()
-{
- flag_minE2 = false;
- StatesInTwoRuns = 0;
- sum_2RB = sum_2RG = 0.0;
-
+SIDD_2R::SIDD_2R() : SIDD_1R() {
+ flag_minE2 = false;
+ StatesInTwoRuns = 0;
+ sum_2RB = sum_2RG = 0.0;
}
-SIDD_2R::~SIDD_2R()
-{
-
-}
+SIDD_2R::~SIDD_2R() {}
// search all lowest states for two runs
-bool SIDD_2R::Search_Low2RE()
-{
- long count = 0;
- LstState::iterator i;
- LstState::iterator j;
- flag_minE2 = false; // default value
- sum_2RB = sum_2RG = 0.0;
+bool SIDD_2R::Search_Low2RE() {
+ long count = 0;
+ LstState::iterator i;
+ LstState::iterator j;
+ flag_minE2 = false; // default value
+ sum_2RB = sum_2RG = 0.0;
- for(int n1 = MinWindowSize; n1 <= MaxWindowSize/2; n1++){
- if(Lst_OBE[n1].size() < 1)
- continue;
- if(Lst_OBE[n1].begin()->get_energy() + min_RE + minGres >= max_E)
- continue;
- for(int n2 = n1; n2 <= MaxWindowSize - n1; n2++){
- if(Lst_OBE[n2].size() < 1)
- continue;
- if(Lst_OBE[n1].begin()->get_energy() + Lst_OBE[n2].begin()->get_energy() + Gres[n1+n2][1] >= max_E)
- continue;
+ for (int n1 = MinWindowSize; n1 <= MaxWindowSize / 2; n1++) {
+ if (Lst_OBE[n1].size() < 1)
+ continue;
+ if (Lst_OBE[n1].begin()->get_energy() + min_RE + minGres >= max_E)
+ continue;
+ for (int n2 = n1; n2 <= MaxWindowSize - n1; n2++) {
+ if (Lst_OBE[n2].size() < 1)
+ continue;
+ if (Lst_OBE[n1].begin()->get_energy() +
+ Lst_OBE[n2].begin()->get_energy() + Gres[n1 + n2][1] >=
+ max_E)
+ continue;
- // V: I added this line to calculate the residual superhelical density
- //double a_residual = calc_alpha_res(n1+n2);
+ // V: I added this line to calculate the residual superhelical density
+ // double a_residual = calc_alpha_res(n1+n2);
- for(i = Lst_OBE[n1].begin(); i != Lst_OBE[n1].end(); ++i){
- int p1 = i->get_pos1();
- double e1 = i->get_energy();
+ for (i = Lst_OBE[n1].begin(); i != Lst_OBE[n1].end(); ++i) {
+ int p1 = i->get_pos1();
+ double e1 = i->get_energy();
- if (e1 >= 10000) break;
- j = Lst_OBE[n2].begin();
- if(n1 == n2) ++j; // no repeat of the same group
- if(e1 + j->get_energy() + Gres[n1+n2][1]>= max_E) break;
- for(; j != Lst_OBE[n2].end(); ++j){
- int p2 = j->get_pos1();
- double e2 = j->get_energy();
- if (e2 >= 10000) break;
- double e = e1 + e2 + Gres[n1+n2][1]; //total energy for two runs
- if(e >= max_E) break;
- if(!overlap(p1, p2, n1, n2)){
- if(e < min_E){
- min_E = e;
- max_E = min_E + theta;
- flag_minE2 = true;
- }
- if(e < max_E){
- if(write_profile){
- double exponent= exp(-e/RT + scaling_factor); // V: I added this
- // double exponent= exp(-e/RT);
- update_promatrix(p1, n1, e, exponent);
- update_promatrix(p2, n2, e, exponent);
- sum_2RB += exponent;
- sum_2RG = sum_2RG + e*exponent;
- }
- count++;
- }
- }
- }
- }
- }
- }
- StatesInTwoRuns = count;
- ZsumB += sum_2RB;
- ZsumG +=sum_2RG;
-
- if (results)
- cout << "Number of two-run states = " << count << endl;
- return flag_minE2;
-}
+ if (e1 >= 10000)
+ break;
+ j = Lst_OBE[n2].begin();
+ if (n1 == n2)
+ ++j; // no repeat of the same group
+ if (e1 + j->get_energy() + Gres[n1 + n2][1] >= max_E)
+ break;
+ for (; j != Lst_OBE[n2].end(); ++j) {
+ int p2 = j->get_pos1();
+ double e2 = j->get_energy();
+ if (e2 >= 10000)
+ break;
+ double e = e1 + e2 + Gres[n1 + n2][1]; // total energy for two runs
+ if (e >= max_E)
+ break;
+ if (!overlap(p1, p2, n1, n2)) {
+ if (e < min_E) {
+ min_E = e;
+ max_E = min_E + theta;
+ flag_minE2 = true;
+ }
+ if (e < max_E) {
+ if (write_profile) {
+ double exponent =
+ exp(-e / RT + scaling_factor); // V: I added this
+ // double exponent= exp(-e/RT);
+ update_promatrix(p1, n1, e, exponent);
+ update_promatrix(p2, n2, e, exponent);
+ sum_2RB += exponent;
+ sum_2RG = sum_2RG + e * exponent;
+ }
+ count++;
+ }
+ }
+ }
+ }
+ }
+ }
+ StatesInTwoRuns = count;
+ ZsumB += sum_2RB;
+ ZsumG += sum_2RG;
+ if (results)
+ cout << "Number of two-run states = " << count << endl;
+ return flag_minE2;
+}
// checking if two regions are overlapping
-bool SIDD_2R::overlap(int p1, int p2, int r1, int r2)
-{
- if(r1 + r2 > MaxWindowSize){
- cerr << "out of window size.\n";
- return false;
- }
- if(p2 + r2 <= length_seq - 1 && p1 + r1 <= length_seq - 1){
- if(p1 + r1 < p2 - 1 || p2 + r2 < p1 - 1) // at least one gap
- return false; // no overlapping
- else
- return true;
- }
- else if(p2 + r2 > length_seq - 1 && p1 + r1 <= length_seq - 1){
- if(p1 + r1 < p2 - 1 && r2 - (length_seq - 1 - p2) < p1 - 1)
- return false; // no overlapping
- else
- return true;
- }
- else if(p1 + r1 > length_seq - 1 && p2 + r2 <= length_seq - 1){
- if(p2 + r2 < p1 - 1 && r1 - (length_seq - 1 - p1) < p2 - 1)
- return false; // no overlapping
- else
- return true;
- }
- else
- return true; // overlapping
+bool SIDD_2R::overlap(int p1, int p2, int r1, int r2) {
+ if (r1 + r2 > MaxWindowSize) {
+ cerr << "out of window size.\n";
+ return false;
+ }
+ if (p2 + r2 <= length_seq - 1 && p1 + r1 <= length_seq - 1) {
+ if (p1 + r1 < p2 - 1 || p2 + r2 < p1 - 1) // at least one gap
+ return false; // no overlapping
+ else
+ return true;
+ } else if (p2 + r2 > length_seq - 1 && p1 + r1 <= length_seq - 1) {
+ if (p1 + r1 < p2 - 1 && r2 - (length_seq - 1 - p2) < p1 - 1)
+ return false; // no overlapping
+ else
+ return true;
+ } else if (p1 + r1 > length_seq - 1 && p2 + r2 <= length_seq - 1) {
+ if (p2 + r2 < p1 - 1 && r1 - (length_seq - 1 - p1) < p2 - 1)
+ return false; // no overlapping
+ else
+ return true;
+ } else
+ return true; // overlapping
}
#endif
diff --git a/src/trans_three/SIDD_2R.h b/src/trans_three/SIDD_2R.h
index eb89c9f..8d752c4 100644
--- a/src/trans_three/SIDD_2R.h
+++ b/src/trans_three/SIDD_2R.h
@@ -2,7 +2,7 @@
// This class is designed for two runs derived from SIDD_1R
//
// program to implement algorithm developed by Craig Benham
-//
+//
// author: Chengpeng Bi
// modifiers: Dina Zhabinskaya, Sally Madden, Ian Korf
// compiler: g++
@@ -18,20 +18,16 @@
#include "SIDD_1R.h"
-
-class SIDD_2R : public SIDD_1R
-{
+class SIDD_2R : public SIDD_1R {
protected:
-
- bool flag_minE2;
- long StatesInTwoRuns;
+ bool flag_minE2;
+ long StatesInTwoRuns;
public:
- SIDD_2R();
- virtual ~SIDD_2R();
- bool overlap(int p1, int p2, int r1, int r2);
- bool Search_Low2RE();
-
+ SIDD_2R();
+ virtual ~SIDD_2R();
+ bool overlap(int p1, int p2, int r1, int r2);
+ bool Search_Low2RE();
};
-#endif // !defined(AFX_SIDD_2R_H__B5CD0095_10DE_4C22_A817_256ABBAD9425__INCLUDED_)
+#endif // !defined(AFX_SIDD_2R_H__B5CD0095_10DE_4C22_A817_256ABBAD9425__INCLUDED_)
diff --git a/src/trans_three/SIDD_3R.cpp b/src/trans_three/SIDD_3R.cpp
index ee4cebc..d4c110c 100644
--- a/src/trans_three/SIDD_3R.cpp
+++ b/src/trans_three/SIDD_3R.cpp
@@ -2,7 +2,7 @@
// This class is designed for three runs derived from SIDD_2R
//
// program to implement algorithm developed by Craig Benham
-//
+//
// author: Chengpeng Bi
// modifiers: Dina Zhabinskaya, Sally Madden, Ian Korf
// compiler: g++
@@ -22,90 +22,105 @@
// Construction/Destruction
//////////////////////////////////////////////////////////////////////
-SIDD_3R::SIDD_3R():SIDD_2R()
-{
- flag_minE3 = false;
- StatesInThreeRuns = 0;
- sum_3RB = sum_3RG = 0.0;
+SIDD_3R::SIDD_3R() : SIDD_2R() {
+ flag_minE3 = false;
+ StatesInThreeRuns = 0;
+ sum_3RB = sum_3RG = 0.0;
}
-SIDD_3R::~SIDD_3R()
-{
-
-}
+SIDD_3R::~SIDD_3R() {}
// search all lowest states for three runs
-bool SIDD_3R::Search_Low3RE()
-{
- double count = 0.0;
- LstState::iterator i;
- LstState::iterator k;
- LstState::iterator j;
- flag_minE3 = false;
- sum_3RB = sum_3RG = 0.0;
+bool SIDD_3R::Search_Low3RE() {
+ double count = 0.0;
+ LstState::iterator i;
+ LstState::iterator k;
+ LstState::iterator j;
+ flag_minE3 = false;
+ sum_3RB = sum_3RG = 0.0;
- int window_size = 200;
- double Epsilon = 0.0;
- for(int n1 = 1; n1 <= window_size/3; n1++){
- if(Lst_OBE[n1].begin()->get_energy() + 2*min_RE + minGres >= max_E) continue;
- for(int n2 = n1; n2 <= (window_size - n1)/2 ; n2++){
- if(Lst_OBE[n1].begin()->get_energy() + Lst_OBE[n2].begin()->get_energy()+ min_RE + minGres >= max_E) continue;
- for(int n3 = n2; n3 <= window_size - n1 - n2; n3++){
- if(Lst_OBE[n1].begin()->get_energy() + Lst_OBE[n2].begin()->get_energy() + Lst_OBE[n3].begin()->get_energy() + Gres[n1+n2+n3][2]>= max_E) continue;
+ int window_size = 200;
+ double Epsilon = 0.0;
+ for (int n1 = 1; n1 <= window_size / 3; n1++) {
+ if (Lst_OBE[n1].begin()->get_energy() + 2 * min_RE + minGres >= max_E)
+ continue;
+ for (int n2 = n1; n2 <= (window_size - n1) / 2; n2++) {
+ if (Lst_OBE[n1].begin()->get_energy() +
+ Lst_OBE[n2].begin()->get_energy() + min_RE + minGres >=
+ max_E)
+ continue;
+ for (int n3 = n2; n3 <= window_size - n1 - n2; n3++) {
+ if (Lst_OBE[n1].begin()->get_energy() +
+ Lst_OBE[n2].begin()->get_energy() +
+ Lst_OBE[n3].begin()->get_energy() + Gres[n1 + n2 + n3][2] >=
+ max_E)
+ continue;
- // V: I added this line to calculate the residual superhelical density
- //double a_residual = calc_alpha_res(n1+n2+n3);
+ // V: I added this line to calculate the residual superhelical density
+ // double a_residual = calc_alpha_res(n1+n2+n3);
- for(i = Lst_OBE[n1].begin(); i != Lst_OBE[n1].end(); ++i){
- int p1 = i->get_pos1();
- double e1 = i->get_energy();
- if (e1 >= 10000) break;
- j = Lst_OBE[n2].begin();
- if(n1 == n2) ++j; // no repeat of the same group
- if(e1 + j->get_energy() + Lst_OBE[n3].begin()->get_energy() + Gres[n1+n2+n3][2] >= max_E) break;
- for(; j != Lst_OBE[n2].end(); ++j){
- int p2 = j->get_pos1();
- double e2 = j->get_energy();
- if (e2 >= 10000) break;
- k = Lst_OBE[n3].begin();
- if(n2 == n3) ++k;
- if(e1 + e2 + k->get_energy() + Gres[n1+n2+n3][2] >= max_E) break;
- for(; k != Lst_OBE[n3].end(); ++k){
- int p3 = k->get_pos1();
- double e3 = k->get_energy();
- if (e3 >= 10000) break;
- double e = e1 + e2 + e3 +Gres[n1+n2+n3][2];
- if(e >= max_E) break;
- //checking for overlaps
- if(!overlap(p1, p2, n1, n2) && !overlap(p1, p3, n1, n3) && !overlap(p2, p3, n2, n3)){
- if(min_E - e >=Epsilon){
- min_E = e;
- max_E = e + theta;
- }
- count++;
- if(write_profile){
- // double exponent = exp(-e/RT);
- double exponent= exp(-e/RT + scaling_factor); // V: I added this
+ for (i = Lst_OBE[n1].begin(); i != Lst_OBE[n1].end(); ++i) {
+ int p1 = i->get_pos1();
+ double e1 = i->get_energy();
+ if (e1 >= 10000)
+ break;
+ j = Lst_OBE[n2].begin();
+ if (n1 == n2)
+ ++j; // no repeat of the same group
+ if (e1 + j->get_energy() + Lst_OBE[n3].begin()->get_energy() +
+ Gres[n1 + n2 + n3][2] >=
+ max_E)
+ break;
+ for (; j != Lst_OBE[n2].end(); ++j) {
+ int p2 = j->get_pos1();
+ double e2 = j->get_energy();
+ if (e2 >= 10000)
+ break;
+ k = Lst_OBE[n3].begin();
+ if (n2 == n3)
+ ++k;
+ if (e1 + e2 + k->get_energy() + Gres[n1 + n2 + n3][2] >= max_E)
+ break;
+ for (; k != Lst_OBE[n3].end(); ++k) {
+ int p3 = k->get_pos1();
+ double e3 = k->get_energy();
+ if (e3 >= 10000)
+ break;
+ double e = e1 + e2 + e3 + Gres[n1 + n2 + n3][2];
+ if (e >= max_E)
+ break;
+ // checking for overlaps
+ if (!overlap(p1, p2, n1, n2) && !overlap(p1, p3, n1, n3) &&
+ !overlap(p2, p3, n2, n3)) {
+ if (min_E - e >= Epsilon) {
+ min_E = e;
+ max_E = e + theta;
+ }
+ count++;
+ if (write_profile) {
+ // double exponent = exp(-e/RT);
+ double exponent =
+ exp(-e / RT + scaling_factor); // V: I added this
- update_promatrix(p1, n1, e, exponent);
- update_promatrix(p2, n2, e, exponent);
- update_promatrix(p3, n3, e, exponent);
- sum_3RB += exponent;
- sum_3RG = sum_3RG + e*exponent;
- }
- }
- }
- }
- }
- }
- }
- }
- if (results)
- cout << "Number of three-run states = " << count << endl;
- StatesInThreeRuns = (long)count;
- ZsumB += sum_3RB;
- ZsumG += sum_3RG;
- return flag_minE3;
+ update_promatrix(p1, n1, e, exponent);
+ update_promatrix(p2, n2, e, exponent);
+ update_promatrix(p3, n3, e, exponent);
+ sum_3RB += exponent;
+ sum_3RG = sum_3RG + e * exponent;
+ }
+ }
+ }
+ }
+ }
+ }
+ }
+ }
+ if (results)
+ cout << "Number of three-run states = " << count << endl;
+ StatesInThreeRuns = (long)count;
+ ZsumB += sum_3RB;
+ ZsumG += sum_3RG;
+ return flag_minE3;
}
#endif
diff --git a/src/trans_three/SIDD_3R.h b/src/trans_three/SIDD_3R.h
index 4fad84b..c8ce00f 100644
--- a/src/trans_three/SIDD_3R.h
+++ b/src/trans_three/SIDD_3R.h
@@ -2,7 +2,7 @@
// This class is designed for three runs derived from SIDD_2R
//
// program to implement algorithm developed by Craig Benham
-//
+//
// author: Chengpeng Bi
// modifiers: Dina Zhabinskaya, Sally Madden, Ian Korf
// compiler: g++
@@ -18,16 +18,15 @@
#include "SIDD_2R.h"
-class SIDD_3R : public SIDD_2R
-{
+class SIDD_3R : public SIDD_2R {
protected:
- bool flag_minE3;
- long StatesInThreeRuns;
+ bool flag_minE3;
+ long StatesInThreeRuns;
public:
- SIDD_3R();
- virtual ~SIDD_3R();
- bool Search_Low3RE();
+ SIDD_3R();
+ virtual ~SIDD_3R();
+ bool Search_Low3RE();
};
#endif // !defined(AFX_SIDD_3R_H__B5CD0095_10DE_4C22_A817_256ABBAD9425__INCLUDED_)
diff --git a/src/trans_three/SIDD_4R.cpp b/src/trans_three/SIDD_4R.cpp
index 729d8aa..9187a93 100644
--- a/src/trans_three/SIDD_4R.cpp
+++ b/src/trans_three/SIDD_4R.cpp
@@ -2,7 +2,7 @@
// This class is designed for four runs derived from SIDD_3R
//
// program to implement algorithm developed by Craig Benham
-//
+//
// author: Chengpeng Bi
// modifiers: Dina Zhabinskaya, Sally Madden, Ian Korf
// compiler: g++
@@ -22,124 +22,155 @@
// Construction/Destruction
//////////////////////////////////////////////////////////////////////
-SIDD_4R::SIDD_4R():SIDD_3R()
-{
- flag_minE4 = false;
- StatesInFourRuns = 0;
- sum_4RB = sum_4RG = 0.0;
+SIDD_4R::SIDD_4R() : SIDD_3R() {
+ flag_minE4 = false;
+ StatesInFourRuns = 0;
+ sum_4RB = sum_4RG = 0.0;
}
-SIDD_4R::~SIDD_4R()
-{
-
-}
+SIDD_4R::~SIDD_4R() {}
// search all lowest states for three runs
-bool SIDD_4R::Search_Low4RE()
-{
- double count = 0.0;
- LstState::iterator i;
- LstState::iterator k;
- LstState::iterator j;
- LstState::iterator l;
- flag_minE4 = false;
- sum_4RB = sum_4RG = 0.0;
- int window_size = 200;
- double Epsilon = 0.0;
- for(int n1 = 1; n1 <= window_size/4; n1++){
- if(Lst_OBE[n1].begin()->get_energy() + 3*min_RE + minGres >= max_E)
- continue;
- for(int n2 = n1; n2 <= (window_size - n1)/3 ; n2++){
- if(Lst_OBE[n1].begin()->get_energy() + Lst_OBE[n2].begin()->get_energy()+ 2*min_RE + minGres >= max_E)
- continue;
- for(int n3 = n2; n3 <= (window_size - n1 - n2)/2; n3++){
- if(Lst_OBE[n1].begin()->get_energy() + Lst_OBE[n2].begin()->get_energy() + Lst_OBE[n3].begin()->get_energy() + min_RE + minGres >= max_E)
- continue;
- for(int n4 = n3; n4 <= (window_size - n1 - n2-n3); n4++){
- if(Lst_OBE[n1].begin()->get_energy() +Lst_OBE[n2].begin()->get_energy() + Lst_OBE[n3].begin()->get_energy() + Lst_OBE[n4].begin()->get_energy() + Gres[n1+n2+n3+n4][3]>= max_E)
- continue;
+bool SIDD_4R::Search_Low4RE() {
+ double count = 0.0;
+ LstState::iterator i;
+ LstState::iterator k;
+ LstState::iterator j;
+ LstState::iterator l;
+ flag_minE4 = false;
+ sum_4RB = sum_4RG = 0.0;
+ int window_size = 200;
+ double Epsilon = 0.0;
+ for (int n1 = 1; n1 <= window_size / 4; n1++) {
+ if (Lst_OBE[n1].begin()->get_energy() + 3 * min_RE + minGres >= max_E)
+ continue;
+ for (int n2 = n1; n2 <= (window_size - n1) / 3; n2++) {
+ if (Lst_OBE[n1].begin()->get_energy() +
+ Lst_OBE[n2].begin()->get_energy() + 2 * min_RE + minGres >=
+ max_E)
+ continue;
+ for (int n3 = n2; n3 <= (window_size - n1 - n2) / 2; n3++) {
+ if (Lst_OBE[n1].begin()->get_energy() +
+ Lst_OBE[n2].begin()->get_energy() +
+ Lst_OBE[n3].begin()->get_energy() + min_RE + minGres >=
+ max_E)
+ continue;
+ for (int n4 = n3; n4 <= (window_size - n1 - n2 - n3); n4++) {
+ if (Lst_OBE[n1].begin()->get_energy() +
+ Lst_OBE[n2].begin()->get_energy() +
+ Lst_OBE[n3].begin()->get_energy() +
+ Lst_OBE[n4].begin()->get_energy() +
+ Gres[n1 + n2 + n3 + n4][3] >=
+ max_E)
+ continue;
- // V: I added this line to calculate the residual superhelical density
- // double a_residual = calc_alpha_res(n1+n2+n3+n4);
+ // V: I added this line to calculate the residual superhelical density
+ // double a_residual = calc_alpha_res(n1+n2+n3+n4);
- for(i = Lst_OBE[n1].begin(); i != Lst_OBE[n1].end(); ++i){
- int p1 = i->get_pos1();
- double e1 = i->get_energy();
- if (e1 >= 10000) break;
- j = Lst_OBE[n2].begin();
- if(n1 == n2) ++j; // no repeat of the same group
- if(e1 + j->get_energy() + Lst_OBE[n3].begin()->get_energy() + Lst_OBE[n4].begin()->get_energy() + Gres[n1+n2+n3+n4][3] >= max_E) break;
- for(; j != Lst_OBE[n2].end(); ++j){
- int p2 = j->get_pos1();
- double e2 = j->get_energy();
- if (e2 >= 10000) break;
- k = Lst_OBE[n3].begin();
- if(n2 == n3) ++k;
- if(e1 + e2 + k->get_energy() + Lst_OBE[n4].begin()->get_energy() + Gres[n1+n2+n3+n4][3] >= max_E) break;
- for(; k != Lst_OBE[n3].end(); ++k){
- int p3 = k->get_pos1();
- double e3 = k->get_energy();
- if (e3 >= 10000) break;
- l = Lst_OBE[n4].begin();
- if(n3 == n4) ++l;
- if(e1 + e2 + e3 + l->get_energy() + Gres[n1+n2+n3+n4][3] >= max_E) break;
- for(; l != Lst_OBE[n4].end(); ++l){
- int p4 = l->get_pos1();
- double e4 = l->get_energy();
- if (e4 >= 10000) break;
- double e = e1 + e2 + e3 + e4 +Gres[n1+n2+n3+n4][3];
- if(e >= max_E) break;
- // checking for overlaps
- if(!overlap(p1, p2, n1, n2) && !overlap(p1, p3, n1, n3) && !overlap(p1, p4, n1, n4) && !overlap(p2, p3, n2, n3) && !overlap(p2, p4, n2, n4) && !overlap(p3, p4, n3, n4)){
- if(min_E - e >=Epsilon){
- min_E = e;
- max_E = e + theta;
- }
- count++;
- if(write_profile){
- // double exponent = exp(-e/RT);
- double exponent= exp(-e/RT + scaling_factor); // V: I added this
- update_promatrix(p1, n1, e, exponent);
- update_promatrix(p2, n2, e, exponent);
- update_promatrix(p3, n3, e, exponent);
- update_promatrix(p4, n4, e, exponent);
- sum_4RB += exponent;
- sum_4RG = sum_4RG + e*exponent;
- }
- }
- }
- }
- }
+ for (i = Lst_OBE[n1].begin(); i != Lst_OBE[n1].end(); ++i) {
+ int p1 = i->get_pos1();
+ double e1 = i->get_energy();
+ if (e1 >= 10000)
+ break;
+ j = Lst_OBE[n2].begin();
+ if (n1 == n2)
+ ++j; // no repeat of the same group
+ if (e1 + j->get_energy() + Lst_OBE[n3].begin()->get_energy() +
+ Lst_OBE[n4].begin()->get_energy() +
+ Gres[n1 + n2 + n3 + n4][3] >=
+ max_E)
+ break;
+ for (; j != Lst_OBE[n2].end(); ++j) {
+ int p2 = j->get_pos1();
+ double e2 = j->get_energy();
+ if (e2 >= 10000)
+ break;
+ k = Lst_OBE[n3].begin();
+ if (n2 == n3)
+ ++k;
+ if (e1 + e2 + k->get_energy() +
+ Lst_OBE[n4].begin()->get_energy() +
+ Gres[n1 + n2 + n3 + n4][3] >=
+ max_E)
+ break;
+ for (; k != Lst_OBE[n3].end(); ++k) {
+ int p3 = k->get_pos1();
+ double e3 = k->get_energy();
+ if (e3 >= 10000)
+ break;
+ l = Lst_OBE[n4].begin();
+ if (n3 == n4)
+ ++l;
+ if (e1 + e2 + e3 + l->get_energy() +
+ Gres[n1 + n2 + n3 + n4][3] >=
+ max_E)
+ break;
+ for (; l != Lst_OBE[n4].end(); ++l) {
+ int p4 = l->get_pos1();
+ double e4 = l->get_energy();
+ if (e4 >= 10000)
+ break;
+ double e = e1 + e2 + e3 + e4 + Gres[n1 + n2 + n3 + n4][3];
+ if (e >= max_E)
+ break;
+ // checking for overlaps
+ if (!overlap(p1, p2, n1, n2) && !overlap(p1, p3, n1, n3) &&
+ !overlap(p1, p4, n1, n4) && !overlap(p2, p3, n2, n3) &&
+ !overlap(p2, p4, n2, n4) && !overlap(p3, p4, n3, n4)) {
+ if (min_E - e >= Epsilon) {
+ min_E = e;
+ max_E = e + theta;
}
+ count++;
+ if (write_profile) {
+ // double exponent = exp(-e/RT);
+ double exponent =
+ exp(-e / RT + scaling_factor); // V: I added this
+ update_promatrix(p1, n1, e, exponent);
+ update_promatrix(p2, n2, e, exponent);
+ update_promatrix(p3, n3, e, exponent);
+ update_promatrix(p4, n4, e, exponent);
+ sum_4RB += exponent;
+ sum_4RG = sum_4RG + e * exponent;
+ }
+ }
}
+ }
}
+ }
}
+ }
}
- StatesInFourRuns = (long)count;
- ZsumB += sum_4RB;
- ZsumG += sum_4RG;
+ }
+ StatesInFourRuns = (long)count;
+ ZsumB += sum_4RB;
+ ZsumG += sum_4RG;
- // V: Scaling factor included
- double prob = (ZsumB-exp(-alpha*alpha*K/2/RT + scaling_factor))/ZsumB;
-
- StatesInFourRuns = (long)count;
- totalFreq = StatesInOneRun + StatesInTwoRuns + StatesInThreeRuns+StatesInFourRuns;
- double runs= (sum_1RB+2*sum_2RB+3*sum_3RB+4*sum_4RB)/ZsumB;
- if (results) {
- cout << "Number of four-run states = " << count << endl;
- cout << "Total number of states = " << totalFreq << endl;
- // cout << "Total Partition Function = " << ZsumB << endl;
- // V: We need to multiply by the scaling factor to get the actual partition funciton.
- cout << "Total Partition Function = " << exp(-scaling_factor)*ZsumB << endl;
- // V: And just in case I'll print the scaling factor.
- cout << "Scaling factor " << scaling_factor << endl;
- cout << "Average number of runs = " << runs << endl;
- cout << "Transition Probability = " << prob << endl;
+ // V: Scaling factor included
+ double prob =
+ (ZsumB - exp(-alpha * alpha * K / 2 / RT + scaling_factor)) / ZsumB;
- }
- totalFreq = StatesInOneRun + StatesInTwoRuns + StatesInThreeRuns+StatesInFourRuns;
- flag_minE4 = 0;
- return flag_minE4;
+ StatesInFourRuns = (long)count;
+ totalFreq =
+ StatesInOneRun + StatesInTwoRuns + StatesInThreeRuns + StatesInFourRuns;
+ double runs = (sum_1RB + 2 * sum_2RB + 3 * sum_3RB + 4 * sum_4RB) / ZsumB;
+ if (results) {
+ cout << "Number of four-run states = " << count << endl;
+ cout << "Total number of states = " << totalFreq << endl;
+ // cout << "Total Partition Function = " << ZsumB << endl;
+ // V: We need to multiply by the scaling factor to get the actual partition
+ // funciton.
+ cout << "Total Partition Function = " << exp(-scaling_factor) * ZsumB
+ << endl;
+ // V: And just in case I'll print the scaling factor.
+ cout << "Scaling factor " << scaling_factor << endl;
+ cout << "Average number of runs = " << runs << endl;
+ cout << "Transition Probability = " << prob << endl;
+ }
+ totalFreq =
+ StatesInOneRun + StatesInTwoRuns + StatesInThreeRuns + StatesInFourRuns;
+ flag_minE4 = 0;
+ return flag_minE4;
}
#endif
diff --git a/src/trans_three/SIDD_4R.h b/src/trans_three/SIDD_4R.h
index a397b87..b665cf1 100644
--- a/src/trans_three/SIDD_4R.h
+++ b/src/trans_three/SIDD_4R.h
@@ -2,7 +2,7 @@
// This class is designed for four runs derived from SIDD_3R
//
// program to implement algorithm developed by Craig Benham
-//
+//
// author: Chengpeng Bi
// modifiers: Dina Zhabinskaya, Sally Madden, Ian Korf
// compiler: g++
@@ -13,22 +13,20 @@
// UC Davis Genome Center
//////////////////////////////////////////////////////////////////////////////
-
#ifndef SIDD_4R_H_
#define SIDD_4R_H_
#include "SIDD_3R.h"
-class SIDD_4R : public SIDD_3R
-{
+class SIDD_4R : public SIDD_3R {
private:
- bool flag_minE4;
- long StatesInFourRuns;
+ bool flag_minE4;
+ long StatesInFourRuns;
public:
- SIDD_4R();
- virtual ~SIDD_4R();
- bool Search_Low4RE();
+ SIDD_4R();
+ virtual ~SIDD_4R();
+ bool Search_Low4RE();
};
#endif // !defined(AFX_SIDD_3R_H__B5CD0095_10DE_4C22_A817_256ABBAD9425__INCLUDED_)
diff --git a/src/trans_three/SIDD_Base.cpp b/src/trans_three/SIDD_Base.cpp
index b84ee1d..7600108 100644
--- a/src/trans_three/SIDD_Base.cpp
+++ b/src/trans_three/SIDD_Base.cpp
@@ -1,9 +1,9 @@
// SIDD_Base.cpp: implementation of the SIDD_Base class.
-// This class is the very base one defining all the parameters and all
+// This class is the very base one defining all the parameters and all
// the common data and function members
//
// program to implement algorithm developed by Craig Benham
-//
+//
// author: Chengpeng Bi
// modifiers: Dina Zhabinskaya, Sally Madden, Ian Korf
// compiler: g++
@@ -14,7 +14,6 @@
// UC Davis Genome Center
//////////////////////////////////////////////////////////////////////////////
-
#ifndef SIDD_BASE_CPP
#define SIDD_BASE_CPP
@@ -27,984 +26,886 @@ using namespace std;
// Construction/Destruction
//////////////////////////////////////////////////////////////////////
-SIDD_Base::SIDD_Base()
-{
- plasmid_seq = 0;
- en_cruciforms = 0;
- profile = 0;
- for(int m = 0; m <= MaxInitialWindowSize; m++)
- promatrix[m] = 0;
- R = 8.314/(4.2*1000.0); // gas constant
- C = 3.6; // stiffness constant
- A = 10.4; // bps per turn in B-DNA
- Az = -12.0; //bps per turn in Z-DNA
- tz = 0.4; //undertwist at the two junctions of Z-DNA
- alpha = 0.0; // linking difference
- K0 = 2220.0; // linking energy coefficient
- pea = 0; // null pointer
- Flag_PEA = false;
- EnergyType = Copolymeric; // default value
- MoleculeType = Linear; // default value
- min_E = 0.0; // minimum energy
- max_E = 0.0;
- min_WS = 0; // the window size corresponding to the minimum energy
- write_profile = false; // flagging if writing starts
- ZsumB = 0.0; // summation of Boltzman frequencies
- ZsumG = 0.0; // store partition function value
- sum_1RG = sum_2RG = sum_3RG = 0.0;
- sum_1RB = sum_2RB = sum_3RB = 0.0;
- set_showbase(1);
-
-}
-
-SIDD_Base::~SIDD_Base()
-{
- delete [] plasmid_seq;
- delete [] sa;
- delete [] profile;
- delete [] pea;
- for(int i = 0; i < length_seq; i++) {
- delete [] en_cruciforms[i];
- }
- for(int i = MinWindowSize; i <= MaxWindowSize; i++)
- delete [] promatrix[i];
+SIDD_Base::SIDD_Base() {
+ plasmid_seq = 0;
+ en_cruciforms = 0;
+ profile = 0;
+ for (int m = 0; m <= MaxInitialWindowSize; m++)
+ promatrix[m] = 0;
+ R = 8.314 / (4.2 * 1000.0); // gas constant
+ C = 3.6; // stiffness constant
+ A = 10.4; // bps per turn in B-DNA
+ Az = -12.0; // bps per turn in Z-DNA
+ tz = 0.4; // undertwist at the two junctions of Z-DNA
+ alpha = 0.0; // linking difference
+ K0 = 2220.0; // linking energy coefficient
+ pea = 0; // null pointer
+ Flag_PEA = false;
+ EnergyType = Copolymeric; // default value
+ MoleculeType = Linear; // default value
+ min_E = 0.0; // minimum energy
+ max_E = 0.0;
+ min_WS = 0; // the window size corresponding to the minimum energy
+ write_profile = false; // flagging if writing starts
+ ZsumB = 0.0; // summation of Boltzman frequencies
+ ZsumG = 0.0; // store partition function value
+ sum_1RG = sum_2RG = sum_3RG = 0.0;
+ sum_1RB = sum_2RB = sum_3RB = 0.0;
+ set_showbase(1);
}
-void SIDD_Base::set_stress_level(double stress_level)
-{
- if(stress_level > 0) stress_level = -stress_level;
- supdensity = stress_level;
+SIDD_Base::~SIDD_Base() {
+ delete[] plasmid_seq;
+ delete[] sa;
+ delete[] profile;
+ delete[] pea;
+ for (int i = 0; i < length_seq; i++) {
+ delete[] en_cruciforms[i];
+ }
+ for (int i = MinWindowSize; i <= MaxWindowSize; i++)
+ delete[] promatrix[i];
}
-void SIDD_Base::set_threshold(double threshold)
-{
- theta = threshold;
+void SIDD_Base::set_stress_level(double stress_level) {
+ if (stress_level > 0)
+ stress_level = -stress_level;
+ supdensity = stress_level;
}
-void SIDD_Base::set_salt_conc(double salt_conc)
-{
- Salt_Conc = salt_conc;
- TMAT = 354.65 + 16.6*log10(Salt_Conc);
- TMGC = TMAT + 41.0;
-}
+void SIDD_Base::set_threshold(double threshold) { theta = threshold; }
-void SIDD_Base::set_temperature(double temperature)
-{
- Temperature = temperature;
- RT = 1.9872*Temperature/1000.0; // constant
- BAT = 7.2464*(1.0 - Temperature / TMAT); // coefficient (kcal)
- BGC = 9.0172*(1.0 - Temperature / TMGC); // coefficient (kcal)
+void SIDD_Base::set_salt_conc(double salt_conc) {
+ Salt_Conc = salt_conc;
+ TMAT = 354.65 + 16.6 * log10(Salt_Conc);
+ TMGC = TMAT + 41.0;
}
-void SIDD_Base::set_Cinitiation()
-{
- Ecr = 192.5-Temperature*0.565-4*BAT-2*2.44*RT*log(4);
+void SIDD_Base::set_temperature(double temperature) {
+ Temperature = temperature;
+ RT = 1.9872 * Temperature / 1000.0; // constant
+ BAT = 7.2464 * (1.0 - Temperature / TMAT); // coefficient (kcal)
+ BGC = 9.0172 * (1.0 - Temperature / TMGC); // coefficient (kcal)
}
-
-void SIDD_Base::set_EnergyType(Energetics et)
-{
- EnergyType = et;
+void SIDD_Base::set_Cinitiation() {
+ Ecr = 192.5 - Temperature * 0.565 - 4 * BAT - 2 * 2.44 * RT * log(4);
}
-void SIDD_Base::set_min_e(double min_e)
-{
- min_E = min_e;
- max_E = min_E + theta;
-}
+void SIDD_Base::set_EnergyType(Energetics et) { EnergyType = et; }
-void SIDD_Base::set_MoleculeType(Molecule mt)
-{
- MoleculeType = mt;
+void SIDD_Base::set_min_e(double min_e) {
+ min_E = min_e;
+ max_E = min_E + theta;
}
-void SIDD_Base::set_MaxWindowSize()
-{
- MaxWindowSize = 250;
- if (length_seq < MaxWindowSize){
- MaxWindowSize = length_seq;
- }
+void SIDD_Base::set_MoleculeType(Molecule mt) { MoleculeType = mt; }
+
+void SIDD_Base::set_MaxWindowSize() {
+ MaxWindowSize = 250;
+ if (length_seq < MaxWindowSize) {
+ MaxWindowSize = length_seq;
+ }
}
-void SIDD_Base::set_MinWindowSize()
-{
- if (EnergyType == Z_DNA)
- MinWindowSize = 8;
- else
- MinWindowSize = 1;
+void SIDD_Base::set_MinWindowSize() {
+ if (EnergyType == Z_DNA)
+ MinWindowSize = 8;
+ else
+ MinWindowSize = 1;
}
-void SIDD_Base::set_K()
-{
- K = K0*RT / length_seq; // linking energy coefficient
+void SIDD_Base::set_K() {
+ K = K0 * RT / length_seq; // linking energy coefficient
}
// junction energies
-void SIDD_Base::set_junction()
-{
- if (EnergyType==Z_DNA)
- a=10.0;
- else if (EnergyType==Cruciform)
- a=0; //in string IR
- else
- a = 10.84;
+void SIDD_Base::set_junction() {
+ if (EnergyType == Z_DNA)
+ a = 10.0;
+ else if (EnergyType == Cruciform)
+ a = 0; // in string IR
+ else
+ a = 10.84;
}
-void SIDD_Base::set_Alpha()
-{
- alpha = supdensity*length_seq/A;
-}
+void SIDD_Base::set_Alpha() { alpha = supdensity * length_seq / A; }
-int SIDD_Base::get_sequence_length(){
- if(MoleculeType == Linear)
- return (length_seq - Len_Gap_Seq);
- return length_seq;
+int SIDD_Base::get_sequence_length() {
+ if (MoleculeType == Linear)
+ return (length_seq - Len_Gap_Seq);
+ return length_seq;
}
-void SIDD_Base::set_sequence_length(int len){
- length_seq = len;
- if(MoleculeType != Circular)
- length_seq+=Len_Gap_Seq;
+void SIDD_Base::set_sequence_length(int len) {
+ length_seq = len;
+ if (MoleculeType != Circular)
+ length_seq += Len_Gap_Seq;
}
-void SIDD_Base::set_cruciform_string(std::string cruciform_string)
-{
- cr_string = cruciform_string;
+void SIDD_Base::set_cruciform_string(std::string cruciform_string) {
+ cr_string = cruciform_string;
}
-void SIDD_Base::set_showres(int showres)
-{
- results = showres;
-}
+void SIDD_Base::set_showres(int showres) { results = showres; }
+// read sequence
+double SIDD_Base::prepare_sequence(std::string sequence) {
+ plasmid_seq = new int[length_seq + 1];
+ for (int j = 0; j < length_seq; j++) {
+ plasmid_seq[j] = 0;
+ }
+ int read = 1;
+ int i = 0;
+ string::iterator iter;
+ for (iter = sequence.begin(); iter != sequence.end(); iter++) {
+ char ch = *iter;
+ if (ch == '>') {
+ read = 0;
+ }
+ if (!read) {
+ if (ch == '\n')
+ read = 1;
+ continue;
+ }
+ if (read) {
+ ch = toupper(ch);
+ switch (ch) {
+ case 'A':
+ plasmid_seq[i++] = 0;
+ break;
+
+ case 'C':
+ plasmid_seq[i++] = 1;
+ break;
+ case 'G':
+ plasmid_seq[i++] = 2;
+ break;
+ case 'T':
+ plasmid_seq[i++] = 3;
+ break;
+
+ case 'N':
+ plasmid_seq[i++] = 2;
+ break;
+ default:
+ break;
+ }
+ }
+ }
+ length_seq = i;
+ if (!length_seq) {
+ cerr << "sequence file not available.\n";
+ return false;
+ }
+ if (MoleculeType != Circular)
+ length_seq += Len_Gap_Seq;
-//read sequence
-double SIDD_Base::prepare_sequence(std::string sequence) {
- plasmid_seq = new int[length_seq+1];
- for(int j = 0; j < length_seq; j++){
- plasmid_seq[j] = 0;
- }
- int read = 1;
- int i = 0;
- string::iterator iter;
- for(iter=sequence.begin();iter !=sequence.end();iter++)
- {
- char ch = *iter;
- if(ch == '>') {
- read = 0;
- }
- if(!read) {
- if(ch == '\n')
- read = 1;
- continue;
- }
- if(read) {
- ch = toupper(ch);
- switch(ch)
- {
- case 'A':
- plasmid_seq[i++] = 0;
- break;
-
- case 'C':
- plasmid_seq[i++] = 1;
- break;
- case 'G':
- plasmid_seq[i++] = 2;
- break;
- case 'T':
- plasmid_seq[i++] = 3;
- break;
-
- case 'N':
- plasmid_seq[i++] = 2;
- break;
- default:
- break;
- }
- }
- }
- length_seq = i;
- if(!length_seq) {
- cerr << "sequence file not available.\n";
- return false;
- }
-
- if(MoleculeType != Circular)
- length_seq+=Len_Gap_Seq;
-
- return true;
+ return true;
}
// reading IR string
void SIDD_Base::prepare_cruciforms() {
- string str1 = cr_string;
- size_t found1 = 0;
- int IR[2] = {0, 0};
- double energy = 0.0;
- int pos1 = 0;
- int pos2 = 0;
- string str2;
- string str3;
- en_cruciforms = new double* [length_seq+1];
- for(int i = 0; i < length_seq; i++){
- en_cruciforms[i] = new double[MaxWindowSize+1];
+ string str1 = cr_string;
+ size_t found1 = 0;
+ int IR[2] = {0, 0};
+ double energy = 0.0;
+ int pos1 = 0;
+ int pos2 = 0;
+ string str2;
+ string str3;
+ en_cruciforms = new double *[length_seq + 1];
+ for (int i = 0; i < length_seq; i++) {
+ en_cruciforms[i] = new double[MaxWindowSize + 1];
+ }
+ for (int i = 0; i < length_seq; i++) {
+ for (int j = 0; j <= MaxWindowSize; j++) {
+ en_cruciforms[i][j] = 10000;
}
- for(int i = 0; i < length_seq; i++){
- for(int j = 0; j <= MaxWindowSize; j++){
- en_cruciforms[i][j] = 10000;
+ }
+ while (found1 != string::npos) {
+ found1 = str1.find("|", pos1);
+ if (found1 != string::npos) {
+ str2 = str1.substr(pos1, found1 - pos1);
+ pos2 = 0;
+ size_t found2 = 0;
+ int i = 0;
+ while (found2 != string::npos) {
+ found2 = str2.find(",", pos2);
+ if (found2 != string::npos) {
+ str3 = str2.substr(pos2, found2 - pos2);
+ pos2 = int(found2) + 1;
+ if (i == 2)
+ energy = atof(str3.c_str());
+ else if (i < 2)
+ IR[i] = atoi(str3.c_str());
+ i++;
}
+ }
+ if (i == 3 && IR[0] > 0 && IR[0] <= length_seq && IR[1] > 0 &&
+ IR[1] <= MaxWindowSize)
+ en_cruciforms[IR[0] - 1][IR[1]] = energy;
+ pos1 = int(found1) + 1;
}
- while (found1!=string::npos) {
- found1=str1.find("|",pos1);
- if(found1!=string::npos) {
- str2 = str1.substr(pos1,found1-pos1);
- pos2 = 0;
- size_t found2 = 0;
- int i = 0;
- while (found2!=string::npos) {
- found2=str2.find(",",pos2);
- if(found2!=string::npos) {
- str3 = str2.substr(pos2,found2-pos2);
- pos2 = int(found2)+1;
- if (i == 2)
- energy = atof(str3.c_str());
- else if (i < 2)
- IR[i] = atoi(str3.c_str());
- i++;
- }
- }
- if (i == 3 && IR[0] > 0 && IR[0] <= length_seq &&
- IR[1] > 0 && IR[1] <= MaxWindowSize)
- en_cruciforms[IR[0]-1][IR[1]] = energy;
- pos1 = int(found1)+1;
- }
+ }
+}
+
+// nearest neighbor energetics
+void SIDD_Base::set_Delta_G() {
+ double Delta_H[4][4] = {{32.3, 35.8, 36.0, 31.2},
+ {35.8, 35.1, 39.8, 36.0},
+ {36.0, 39.8, 35.1, 35.8},
+ {31.2, 36.0, 35.8, 32.3}};
+
+ double Delta_S[4][4] = {// J per mol per Kelvin degree
+ {95.6, 100.4, 102.4, 95.3},
+ {100.4, 93.9, 102.4, 102.4},
+ {102.4, 102.4, 93.9, 100.4},
+ {95.3, 102.4, 100.4, 95.6}};
+ // unit conversion: KJ to Kcal
+ // 1 kcal = 4.184 KJ
+ int i, j;
+ for (i = 0; i < 4; i++) {
+ for (j = 0; j < 4; j++) {
+ Delta_H[i][j] /= 4.184;
+ Delta_S[i][j] /= 4184.0;
}
-}
+ }
+ for (i = 0; i < 4; i++) {
+ for (int j = 0; j < 4; j++) {
+ Delta_S[i][j] = 1.0 / (16.6 * log10(Salt_Conc / 0.1) / Delta_H[i][j] +
+ 1.0 / Delta_S[i][j]);
+ }
+ }
-//nearest neighbor energetics
-void SIDD_Base::set_Delta_G()
-{
- double Delta_H[4][4] = {
- {32.3, 35.8, 36.0, 31.2},
- {35.8, 35.1, 39.8, 36.0},
- {36.0, 39.8, 35.1, 35.8},
- {31.2, 36.0, 35.8, 32.3}
- };
-
- double Delta_S[4][4] ={ // J per mol per Kelvin degree
- {95.6, 100.4, 102.4, 95.3},
- {100.4, 93.9, 102.4, 102.4},
- {102.4, 102.4, 93.9, 100.4},
- {95.3, 102.4, 100.4, 95.6}
- };
- // unit conversion: KJ to Kcal
- // 1 kcal = 4.184 KJ
- int i, j;
- for(i = 0; i < 4; i++){
- for(j = 0; j < 4; j++){
- Delta_H[i][j] /= 4.184;
- Delta_S[i][j] /= 4184.0;
- }
- }
-
- for(i = 0; i < 4; i++){
- for(int j = 0; j < 4; j++){
- Delta_S[i][j] = 1.0/(16.6*log10(Salt_Conc/0.1)/Delta_H[i][j] + 1.0 / Delta_S[i][j]);
- }
- }
-
- for(i = 0; i < 4; i++){
- for(int j = 0; j < 4; j++){
- Delta_G[i][j] = Delta_H[i][j] - Temperature*Delta_S[i][j]; // kcal per mole
- }
- }
+ for (i = 0; i < 4; i++) {
+ for (int j = 0; j < 4; j++) {
+ Delta_G[i][j] =
+ Delta_H[i][j] - Temperature * Delta_S[i][j]; // kcal per mole
+ }
+ }
}
-
// reading sequence and initializing parameters
-bool SIDD_Base::initializer(std::string sequence,int len)
-{
- set_sequence_length(len);
- if(!prepare_sequence(sequence)) {
- //sequence either 0 or too long
- return 0;
- }
- set_MaxWindowSize();
- set_MinWindowSize();
- set_K();
- set_Alpha();
- set_junction();
- set_Delta_G();
- set_E_limit();
- prepare_cruciforms();
- sa = new char[length_seq+1];
- pea = new double[length_seq+1];
- for(int i = 0; i < length_seq+1; i++)
- pea[i] = 0.0;
- profile = new G_x[length_seq+1];
- for(int m = MinWindowSize; m <= MaxWindowSize; m++)
- promatrix[m] = new G_x[length_seq+1];
- if(MoleculeType == Linear){
- for(int k = length_seq - Len_Gap_Seq; k < length_seq; k++)
- plasmid_seq[k] = 2; // add 50 Gs to tail when linear
- }
- set_minGres();
- reset_profile();
- if (EnergyType == Z_DNA)
- {
- for(int i = 0; i < length_seq; i++){
- switch(plasmid_seq[i]){
- case 0:
- sa[i]='s';
- break;
- case 1:
- sa[i]='a';
- break;
- case 2:
- sa[i]='s';
- break;
- case 3:
- sa[i]='a';
- break;
- }
- }
- for(int i = 0; i < length_seq; i++){
- int c1 = sa[i-1 < 0 ? length_seq - 1 : i - 1];
- int c2 = sa[i+1 == length_seq ? 0 : i + 1];
- int c = sa[i];
- if (c==c1 && c==c2) {
- if (c=='s')
- sa[i]='a';
- else
- sa[i]='s';
- }
- }
- for(int i = 0; i < length_seq; i++){
- int c2 = sa[i+1 == length_seq ? 0 : i + 1];
- int c = sa[i];
- int c3 = sa[i+2 == length_seq ? 1 : i + 2];
- int c4 = sa[i+3 == length_seq ? 2 : i + 3];
- int c5 = sa[i+4 == length_seq ? 3 : i + 4];
- if (c==c2 && c3==c4 && c!=c3) {
- if (c5=='a') {
- sa[i]='a';
- sa[i+1]='s';
- sa[i+2]='a';
- sa[i+3]='s';
- }
- else {
- sa[i]='s';
- sa[i+1]='a';
- sa[i+2]='s';
- sa[i+3]='a';
- }
- }
- }
- }
- return true;
+bool SIDD_Base::initializer(std::string sequence, int len) {
+ set_sequence_length(len);
+ if (!prepare_sequence(sequence)) {
+ // sequence either 0 or too long
+ return 0;
+ }
+ set_MaxWindowSize();
+ set_MinWindowSize();
+ set_K();
+ set_Alpha();
+ set_junction();
+ set_Delta_G();
+ set_E_limit();
+ prepare_cruciforms();
+ sa = new char[length_seq + 1];
+ pea = new double[length_seq + 1];
+ for (int i = 0; i < length_seq + 1; i++)
+ pea[i] = 0.0;
+ profile = new G_x[length_seq + 1];
+ for (int m = MinWindowSize; m <= MaxWindowSize; m++)
+ promatrix[m] = new G_x[length_seq + 1];
+ if (MoleculeType == Linear) {
+ for (int k = length_seq - Len_Gap_Seq; k < length_seq; k++)
+ plasmid_seq[k] = 2; // add 50 Gs to tail when linear
+ }
+ set_minGres();
+ reset_profile();
+ if (EnergyType == Z_DNA) {
+ for (int i = 0; i < length_seq; i++) {
+ switch (plasmid_seq[i]) {
+ case 0:
+ sa[i] = 's';
+ break;
+ case 1:
+ sa[i] = 'a';
+ break;
+ case 2:
+ sa[i] = 's';
+ break;
+ case 3:
+ sa[i] = 'a';
+ break;
+ }
+ }
+ for (int i = 0; i < length_seq; i++) {
+ int c1 = sa[i - 1 < 0 ? length_seq - 1 : i - 1];
+ int c2 = sa[i + 1 == length_seq ? 0 : i + 1];
+ int c = sa[i];
+ if (c == c1 && c == c2) {
+ if (c == 's')
+ sa[i] = 'a';
+ else
+ sa[i] = 's';
+ }
+ }
+ for (int i = 0; i < length_seq; i++) {
+ int c2 = sa[i + 1 == length_seq ? 0 : i + 1];
+ int c = sa[i];
+ int c3 = sa[i + 2 == length_seq ? 1 : i + 2];
+ int c4 = sa[i + 3 == length_seq ? 2 : i + 3];
+ int c5 = sa[i + 4 == length_seq ? 3 : i + 4];
+ if (c == c2 && c3 == c4 && c != c3) {
+ if (c5 == 'a') {
+ sa[i] = 'a';
+ sa[i + 1] = 's';
+ sa[i + 2] = 'a';
+ sa[i + 3] = 's';
+ } else {
+ sa[i] = 's';
+ sa[i + 1] = 'a';
+ sa[i + 2] = 's';
+ sa[i + 3] = 'a';
+ }
+ }
+ }
+ }
+ return true;
}
-
// assign energy position-wide
-void SIDD_Base::assign_pea(char* efile)
-{
- ifstream ins(efile);
-
- if(ins.is_open()){
- while(!ins.eof()){
- double e0;
- int p;
- ins >> p >> e0;
- if(p >=0 && p < length_seq && e0 != 0.0){
- pea[p] = e0;
- cout << "assign energy at " << p << ": " << pea[p] << endl;
- }
- }
- }
- else{
- cerr << "fail to open the file: " << efile << endl;
- Flag_PEA = false;
- }
- ins.close();
-}
-
-//superhelical energy for each transition
-double SIDD_Base::calc_Gres(int n, int nr)
-{
- if (EnergyType == Z_DNA) {
- double t2 = 0.5*K*(alpha+n/A-n/Az+2*tz*nr)*(alpha+n/A-n/Az+2*tz*nr);
- return t2;
- }
- else if (EnergyType == Cruciform) {
- double t2 = 0.5*K*(alpha + n/A) * (alpha + n/A);
- return t2;
- }
- else{
- double t1 = (alpha + n / A) * (alpha + n / A) ;
- double t2 = 2*PI*PI*C*K*t1;
- return t2 / (4*PI*PI*C + K*n);
+void SIDD_Base::assign_pea(char *efile) {
+ ifstream ins(efile);
+
+ if (ins.is_open()) {
+ while (!ins.eof()) {
+ double e0;
+ int p;
+ ins >> p >> e0;
+ if (p >= 0 && p < length_seq && e0 != 0.0) {
+ pea[p] = e0;
+ cout << "assign energy at " << p << ": " << pea[p] << endl;
+ }
}
+ } else {
+ cerr << "fail to open the file: " << efile << endl;
+ Flag_PEA = false;
+ }
+ ins.close();
+}
+
+// superhelical energy for each transition
+double SIDD_Base::calc_Gres(int n, int nr) {
+ if (EnergyType == Z_DNA) {
+ double t2 = 0.5 * K * (alpha + n / A - n / Az + 2 * tz * nr) *
+ (alpha + n / A - n / Az + 2 * tz * nr);
+ return t2;
+ } else if (EnergyType == Cruciform) {
+ double t2 = 0.5 * K * (alpha + n / A) * (alpha + n / A);
+ return t2;
+ } else {
+ double t1 = (alpha + n / A) * (alpha + n / A);
+ double t2 = 2 * PI * PI * C * K * t1;
+ return t2 / (4 * PI * PI * C + K * n);
+ }
}
// V: I added this function to calculate the residual superhelical density
// n = number of denatured bases
-//double SIDD_Base::calc_alpha_res(int n)
+// double SIDD_Base::calc_alpha_res(int n)
//{/
// if (EnergyType == Z_DNA) {
- // return 0.0;
- // }
- // else if (EnergyType == Cruciform) {
- // return 0.0;
+// return 0.0;
+// }
+// else if (EnergyType == Cruciform) {
+// return 0.0;
// }
// else{
- // double t1 = 4.0*PI*PI*C*(A*alpha + n) ;
- // double t2 = A*(4.0*PI*PI*C + n*K) ;
- // return t1/t2 ;
- // }
+// double t1 = 4.0*PI*PI*C*(A*alpha + n) ;
+// double t2 = A*(4.0*PI*PI*C + n*K) ;
+// return t1/t2 ;
+// }
//}
-
-void SIDD_Base::set_E_limit()
-{
- min_E = 0.5*K*alpha*alpha;
- max_E = min_E + theta;
+void SIDD_Base::set_E_limit() {
+ min_E = 0.5 * K * alpha * alpha;
+ max_E = min_E + theta;
}
// set minimum residual energy
-void SIDD_Base::set_minGres()
-{
- minGres = Gres[MinWindowSize][0] = calc_Gres(MinWindowSize,1);
- for (int r = 0; r < 4; r++) {
- for(int w = MinWindowSize+1; w <= MaxWindowSize; w++){
- Gres[w][r] = calc_Gres(w,r+1);
- if (r == 0) {
- if(Gres[w][r] < minGres)
- minGres = Gres[w][r];
- }
- }
+void SIDD_Base::set_minGres() {
+ minGres = Gres[MinWindowSize][0] = calc_Gres(MinWindowSize, 1);
+ for (int r = 0; r < 4; r++) {
+ for (int w = MinWindowSize + 1; w <= MaxWindowSize; w++) {
+ Gres[w][r] = calc_Gres(w, r + 1);
+ if (r == 0) {
+ if (Gres[w][r] < minGres)
+ minGres = Gres[w][r];
+ }
}
+ }
}
// computing transition energy: Z-DNA
-void SIDD_Base::calc_Z(int pos)
-{
- //Z-DNA energetics
- double energy_Z_AS[4][4] = {
- {3.9,4.6,3.4,5.9},
- {1.3,2.4,0.7,3.4},
- {3.4,4.0,2.4,4.6},
- {2.5,3.4,1.3,3.9}
- };
-
- double energy_Z_SA[4][4] = {
- {3.9,1.3,3.4,2.5},
- {4.6,2.4,4.0,3.4},
- {3.4,0.7,2.4,1.3},
- {5.9,3.4,4.6,3.9}
- };
-
- double energy_Z_ZZ[4][4] = {
- {7.4,4.5,6.3,5.6},
- {4.5,4.0,4.0,6.3},
- {6.3,4.0,4.0,4.5},
- {5.6,6.3,4.5,7.4}
- };
- // Base code before pos
- int p1 = plasmid_seq[pos-1 < 0 ? length_seq - 1 : pos - 1];
- // Base code after pos
- int p2 = plasmid_seq[pos+1 == length_seq ? 0 : pos + 1];
- // Base code of pos
- int p = plasmid_seq[pos];
- // Config code (anti or syn) before pos
- int c1 = sa[pos-1 < 0 ? length_seq - 1 : pos - 1];
- // Config code (anti or syn) after pos
- int c2 = sa[pos+1 == length_seq ? 0 : pos + 1];
- // Config code (anti or syn) of pos
- int c = sa[pos];
- e1=0;
- e2=0;
- zz=0;
- if (c != c2)
- {
- if (c == 'a') {
- e1 = energy_Z_AS[p][p2];
- }
- else {
- e1 = energy_Z_SA[p][p2];
- }
- }
- if (c != c1)
- {
- if (c == 'a') {
- e2 = energy_Z_AS[p][p1];
- }
- else {
- e2 = energy_Z_SA[p][p1];
- }
- }
- if (c == c1) {
- zz =energy_Z_ZZ[p][p1];
- int p0 = sa[pos-2 < 0 ? length_seq - 2 : pos - 2];
- if (energy_Z_ZZ[p][p2] >= energy_Z_ZZ[p0][p1]) {
- if (c=='a')
- e2 = energy_Z_SA[p1][p];
- else
- e2 = energy_Z_AS[p1][p];
- }
- else {
- if (c1=='a')
- e2 = energy_Z_AS[p1][p];
- else
- e2 = energy_Z_SA[p1][p];
- }
-
- }
- if (c == c2) {
- int p3 = plasmid_seq[pos+2 == length_seq ? 1 : pos + 2];
- if (energy_Z_ZZ[p2][p3] >= energy_Z_ZZ[p][p1]) {
- zz = energy_Z_ZZ[p][p1];
- if (c2=='a')
- e1 = energy_Z_SA[p][p2];
- else
- e1 = energy_Z_AS[p][p2];
- }
- else {
- zz = energy_Z_ZZ[p2][p3];
- if (c=='a')
- e1 = energy_Z_AS[p][p2];
- else
- e1 = energy_Z_SA[p][p2];
- }
- }
+void SIDD_Base::calc_Z(int pos) {
+ // Z-DNA energetics
+ double energy_Z_AS[4][4] = {{3.9, 4.6, 3.4, 5.9},
+ {1.3, 2.4, 0.7, 3.4},
+ {3.4, 4.0, 2.4, 4.6},
+ {2.5, 3.4, 1.3, 3.9}};
+
+ double energy_Z_SA[4][4] = {{3.9, 1.3, 3.4, 2.5},
+ {4.6, 2.4, 4.0, 3.4},
+ {3.4, 0.7, 2.4, 1.3},
+ {5.9, 3.4, 4.6, 3.9}};
+
+ double energy_Z_ZZ[4][4] = {{7.4, 4.5, 6.3, 5.6},
+ {4.5, 4.0, 4.0, 6.3},
+ {6.3, 4.0, 4.0, 4.5},
+ {5.6, 6.3, 4.5, 7.4}};
+ // Base code before pos
+ int p1 = plasmid_seq[pos - 1 < 0 ? length_seq - 1 : pos - 1];
+ // Base code after pos
+ int p2 = plasmid_seq[pos + 1 == length_seq ? 0 : pos + 1];
+ // Base code of pos
+ int p = plasmid_seq[pos];
+ // Config code (anti or syn) before pos
+ int c1 = sa[pos - 1 < 0 ? length_seq - 1 : pos - 1];
+ // Config code (anti or syn) after pos
+ int c2 = sa[pos + 1 == length_seq ? 0 : pos + 1];
+ // Config code (anti or syn) of pos
+ int c = sa[pos];
+ e1 = 0;
+ e2 = 0;
+ zz = 0;
+ if (c != c2) {
+ if (c == 'a') {
+ e1 = energy_Z_AS[p][p2];
+ } else {
+ e1 = energy_Z_SA[p][p2];
+ }
+ }
+ if (c != c1) {
+ if (c == 'a') {
+ e2 = energy_Z_AS[p][p1];
+ } else {
+ e2 = energy_Z_SA[p][p1];
+ }
+ }
+ if (c == c1) {
+ zz = energy_Z_ZZ[p][p1];
+ int p0 = sa[pos - 2 < 0 ? length_seq - 2 : pos - 2];
+ if (energy_Z_ZZ[p][p2] >= energy_Z_ZZ[p0][p1]) {
+ if (c == 'a')
+ e2 = energy_Z_SA[p1][p];
+ else
+ e2 = energy_Z_AS[p1][p];
+ } else {
+ if (c1 == 'a')
+ e2 = energy_Z_AS[p1][p];
+ else
+ e2 = energy_Z_SA[p1][p];
+ }
+ }
+ if (c == c2) {
+ int p3 = plasmid_seq[pos + 2 == length_seq ? 1 : pos + 2];
+ if (energy_Z_ZZ[p2][p3] >= energy_Z_ZZ[p][p1]) {
+ zz = energy_Z_ZZ[p][p1];
+ if (c2 == 'a')
+ e1 = energy_Z_SA[p][p2];
+ else
+ e1 = energy_Z_AS[p][p2];
+ } else {
+ zz = energy_Z_ZZ[p2][p3];
+ if (c == 'a')
+ e1 = energy_Z_AS[p][p2];
+ else
+ e1 = energy_Z_SA[p][p2];
+ }
+ }
}
-
// start position - startp
// window size - n
// return - sum of window energy
-double SIDD_Base::sum_WindowZ(int startp, int n)
-{
- double s = 0.0;
- if (n%2==1) {
- s=10000;
- return s;
- }
- else {
- if(startp + n - 1 < length_seq){
- for(int i = startp; i < startp + n; i++) {
- calc_Z(i);
- if (i==startp)
- s+=e1/2;
- if ((i-startp)%2==0 && i!=startp)
- s+=e1/2+zz;
- if ((i-startp)%2==1)
- s+=e2/2;
- }
- return s;
- }
- else{
- for(int i = startp; i < length_seq; i++) {
- calc_Z(i);
- if (i==startp)
- s+=e1/2;
- if ((i-startp)%2==0 && i!=startp)
- s+=e1/2+zz;
- if ((i-startp)%2==1)
- s+=e2/2;
- }
- for(int j = 0; j < startp + n - length_seq; j++) {
- calc_Z(j);
- if (j==startp)
- s+=e1/2;
- if ((length_seq-startp+j)%2==0 && j!=startp)
- s+=e1/2+zz;
- if ((length_seq-startp+j)%2==1)
- s+=e2/2;
- }
- return s;
- }
- }
+double SIDD_Base::sum_WindowZ(int startp, int n) {
+ double s = 0.0;
+ if (n % 2 == 1) {
+ s = 10000;
+ return s;
+ } else {
+ if (startp + n - 1 < length_seq) {
+ for (int i = startp; i < startp + n; i++) {
+ calc_Z(i);
+ if (i == startp)
+ s += e1 / 2;
+ if ((i - startp) % 2 == 0 && i != startp)
+ s += e1 / 2 + zz;
+ if ((i - startp) % 2 == 1)
+ s += e2 / 2;
+ }
+ return s;
+ } else {
+ for (int i = startp; i < length_seq; i++) {
+ calc_Z(i);
+ if (i == startp)
+ s += e1 / 2;
+ if ((i - startp) % 2 == 0 && i != startp)
+ s += e1 / 2 + zz;
+ if ((i - startp) % 2 == 1)
+ s += e2 / 2;
+ }
+ for (int j = 0; j < startp + n - length_seq; j++) {
+ calc_Z(j);
+ if (j == startp)
+ s += e1 / 2;
+ if ((length_seq - startp + j) % 2 == 0 && j != startp)
+ s += e1 / 2 + zz;
+ if ((length_seq - startp + j) % 2 == 1)
+ s += e2 / 2;
+ }
+ return s;
+ }
+ }
}
-
// computing transition energy: neighbor interation
-double SIDD_Base::calc_NI(int pos)
-{
- if(pos < length_seq && pos >= 0){
-
- if(Flag_PEA && pea[pos] != 0.0)
- return pea[pos]; // assign specified energy
- int p1 = plasmid_seq[pos-1 < 0 ? length_seq - 1 : pos - 1];
- int p2 = plasmid_seq[pos+1 == length_seq ? 0 : pos + 1];
- int p = plasmid_seq[pos];
- return (Delta_G[p1][p] + Delta_G[p][p2]) / 2.0;
- }
- else{
- cerr << "position: out of range.\n";
- return 0;
- }
+double SIDD_Base::calc_NI(int pos) {
+ if (pos < length_seq && pos >= 0) {
+
+ if (Flag_PEA && pea[pos] != 0.0)
+ return pea[pos]; // assign specified energy
+ int p1 = plasmid_seq[pos - 1 < 0 ? length_seq - 1 : pos - 1];
+ int p2 = plasmid_seq[pos + 1 == length_seq ? 0 : pos + 1];
+ int p = plasmid_seq[pos];
+ return (Delta_G[p1][p] + Delta_G[p][p2]) / 2.0;
+ } else {
+ cerr << "position: out of range.\n";
+ return 0;
+ }
}
// start position - startp
// window size - n
// return - sum of window energy
-double SIDD_Base::sum_WindowNI(int startp, int n)
-{
- double s = 0.0;
- if(startp + n - 1 < length_seq){
- for(int i = startp; i < startp + n; i++)
- s += calc_NI(i);
- return s;
- }
- else{
- for(int i = startp; i < length_seq; i++)
- s += calc_NI(i);
- for(int j = 0; j < startp + n - length_seq; j++)
- s += calc_NI(j);
- }
+double SIDD_Base::sum_WindowNI(int startp, int n) {
+ double s = 0.0;
+ if (startp + n - 1 < length_seq) {
+ for (int i = startp; i < startp + n; i++)
+ s += calc_NI(i);
+ return s;
+ } else {
+ for (int i = startp; i < length_seq; i++)
+ s += calc_NI(i);
+ for (int j = 0; j < startp + n - length_seq; j++)
+ s += calc_NI(j);
+ }
- return s;
+ return s;
}
-double SIDD_Base::sum_cruciform(int startp, int n)
-{
- return en_cruciforms[startp][n];
+double SIDD_Base::sum_cruciform(int startp, int n) {
+ return en_cruciforms[startp][n];
}
-double SIDD_Base::sum_WindowNIbyBases(int startp, int n)
-{
- int at, gc;
- at = gc = 0;
- double tempE = 0.0;
- count_AT_GC(startp, n, at, gc, tempE);
- return (BAT*at + BGC*gc + tempE);
-
+double SIDD_Base::sum_WindowNIbyBases(int startp, int n) {
+ int at, gc;
+ at = gc = 0;
+ double tempE = 0.0;
+ count_AT_GC(startp, n, at, gc, tempE);
+ return (BAT * at + BGC * gc + tempE);
}
// counting A/T bases
-void SIDD_Base::count_AT_GC(int startp, int n, int& c_AT, int& c_GC, double& tempE)
-{
- c_AT = c_GC = 0;
- tempE = 0.0;
- if(startp + n - 1 < length_seq){
- for(int i = startp; i < startp + n; i++){
- if(Flag_PEA && pea[i] != 0.0)
- tempE += pea[i];
- else{
- if(plasmid_seq[i] == 0 || plasmid_seq[i] == 3)
- c_AT++;
- else
- c_GC++;
- }
- }
- }
-
- else{
- for(int i = startp; i < length_seq; i++){
- if(Flag_PEA && pea[i] != 0.0)
- tempE += pea[i];
- else{
-
- if(plasmid_seq[i] == 0 || plasmid_seq[i] == 3)
- c_AT++;
- else
- c_GC++;
- }
- }
-
- for(int j = 0; j < startp + n - length_seq; j++){
- if(Flag_PEA && pea[j] != 0.0)
- tempE += pea[j];
- else{
- if(plasmid_seq[j] == 0 || plasmid_seq[j] == 3)
- c_AT++;
- else
- c_GC++;
- }
- }
- }
-}
-
-
-double SIDD_Base::calc_OPenBasesEnergy(int startp, int n)
-{
- // base pair energy type
- switch(EnergyType){
- case Z_DNA:
- return sum_WindowZ(startp,n);
- case Cruciform:
- return sum_cruciform(startp,n);
- case Near_Neighbor:
- return sum_WindowNI(startp, n);
- case Copolymeric:
- return sum_WindowNIbyBases(startp, n);
- default:
- return sum_WindowNIbyBases(startp, n);
- }
-
-}
-
-
-void SIDD_Base::reset_promatrix()
-{
- ZsumB = 0.0;
- ZsumG = 0.0;
- totalFreq = 0;
- for(int k = MinWindowSize; k <= MaxWindowSize; k++){
- for(int i = 0; i < length_seq; i++){
- promatrix[k][i].reset();
- }
- }
+void SIDD_Base::count_AT_GC(int startp, int n, int &c_AT, int &c_GC,
+ double &tempE) {
+ c_AT = c_GC = 0;
+ tempE = 0.0;
+ if (startp + n - 1 < length_seq) {
+ for (int i = startp; i < startp + n; i++) {
+ if (Flag_PEA && pea[i] != 0.0)
+ tempE += pea[i];
+ else {
+ if (plasmid_seq[i] == 0 || plasmid_seq[i] == 3)
+ c_AT++;
+ else
+ c_GC++;
+ }
+ }
+ }
+
+ else {
+ for (int i = startp; i < length_seq; i++) {
+ if (Flag_PEA && pea[i] != 0.0)
+ tempE += pea[i];
+ else {
+
+ if (plasmid_seq[i] == 0 || plasmid_seq[i] == 3)
+ c_AT++;
+ else
+ c_GC++;
+ }
+ }
+
+ for (int j = 0; j < startp + n - length_seq; j++) {
+ if (Flag_PEA && pea[j] != 0.0)
+ tempE += pea[j];
+ else {
+ if (plasmid_seq[j] == 0 || plasmid_seq[j] == 3)
+ c_AT++;
+ else
+ c_GC++;
+ }
+ }
+ }
+}
+
+double SIDD_Base::calc_OPenBasesEnergy(int startp, int n) {
+ // base pair energy type
+ switch (EnergyType) {
+ case Z_DNA:
+ return sum_WindowZ(startp, n);
+ case Cruciform:
+ return sum_cruciform(startp, n);
+ case Near_Neighbor:
+ return sum_WindowNI(startp, n);
+ case Copolymeric:
+ return sum_WindowNIbyBases(startp, n);
+ default:
+ return sum_WindowNIbyBases(startp, n);
+ }
+}
+
+void SIDD_Base::reset_promatrix() {
+ ZsumB = 0.0;
+ ZsumG = 0.0;
+ totalFreq = 0;
+ for (int k = MinWindowSize; k <= MaxWindowSize; k++) {
+ for (int i = 0; i < length_seq; i++) {
+ promatrix[k][i].reset();
+ }
+ }
}
// startp - start position
// n - window size
// x - free energy
-bool SIDD_Base::update_promatrix(int startp, int n, double x, double bzfactor)
-{
- if(startp < 0 || startp >= length_seq || n > MaxWindowSize)
- return false;
-
- promatrix[n][startp].add(x, bzfactor, RT);
-// return promatrix[n][startp].get_exponent();
- return true;
-}
-
-void SIDD_Base::fill_profile()
-{
- for(int n = MinWindowSize; n <= MaxWindowSize; n++){
- for(int startp = 0; startp < length_seq; startp++){
- double lastG = promatrix[n][startp].get_sum_xG();
- double lastB = promatrix[n][startp].get_sum_xB();
- if(startp + n - 1 < length_seq){
- for(int i = startp; i < startp + n; i++){
- profile[i].add(-1.0, 0.0, -1.0, lastG, lastB);
- }
- }
- else{
- for(int i = startp; i < length_seq; i++)
- profile[i].add(-1.0, 0.0, -1.0, lastG, lastB);
- for(int k = 0; k < startp + n - length_seq; k++)
- profile[k].add(-1.0, 0.0, -1.0 , lastG, lastB);
- }
- }
- }
-}
-
-
-void SIDD_Base::reset_profile()
-{
- ZsumB = 0.0;
- ZsumG = 0.0;
- totalFreq = 0;
- for(int i = 0; i < length_seq; i++){
- profile[i].reset();
- }
-
- reset_promatrix();
-}
-
-void SIDD_Base::calc_profile()
-{
- double ave_Gs = 0.0;
- if(ZsumB == 0.0) {
- ave_Gs = 0.0;
- }
- else{
- ave_Gs = ZsumG / ZsumB;
+bool SIDD_Base::update_promatrix(int startp, int n, double x, double bzfactor) {
+ if (startp < 0 || startp >= length_seq || n > MaxWindowSize)
+ return false;
+
+ promatrix[n][startp].add(x, bzfactor, RT);
+ // return promatrix[n][startp].get_exponent();
+ return true;
+}
+
+void SIDD_Base::fill_profile() {
+ for (int n = MinWindowSize; n <= MaxWindowSize; n++) {
+ for (int startp = 0; startp < length_seq; startp++) {
+ double lastG = promatrix[n][startp].get_sum_xG();
+ double lastB = promatrix[n][startp].get_sum_xB();
+ if (startp + n - 1 < length_seq) {
+ for (int i = startp; i < startp + n; i++) {
+ profile[i].add(-1.0, 0.0, -1.0, lastG, lastB);
+ }
+ } else {
+ for (int i = startp; i < length_seq; i++)
+ profile[i].add(-1.0, 0.0, -1.0, lastG, lastB);
+ for (int k = 0; k < startp + n - length_seq; k++)
+ profile[k].add(-1.0, 0.0, -1.0, lastG, lastB);
+ }
}
+ }
+}
+
+void SIDD_Base::reset_profile() {
+ ZsumB = 0.0;
+ ZsumG = 0.0;
+ totalFreq = 0;
+ for (int i = 0; i < length_seq; i++) {
+ profile[i].reset();
+ }
- for(int i = 0; i < length_seq; i++){
- profile[i].calc(ave_Gs, ZsumB);
- }
+ reset_promatrix();
+}
- maxGx = profile[0].get_ave_Gx();
- for(int j = 1; j < length_seq; j++){
- if(maxGx < profile[j].get_ave_Gx()){
- maxGx = profile[j].get_ave_Gx(); // find max Gx
- }
- }
+void SIDD_Base::calc_profile() {
+ double ave_Gs = 0.0;
+ if (ZsumB == 0.0) {
+ ave_Gs = 0.0;
+ } else {
+ ave_Gs = ZsumG / ZsumB;
+ }
- for(int k = 0; k < length_seq; k++){
- if(profile[k].get_ave_Gx() == INFINITE_Gx) {
- profile[k].set_ave_Gx(maxGx); // replace INF Gx with maxGx
- }
- }
-}
-
-void SIDD_Base::sum_open()
-{
- double sumn=0.0;
- for(int i = 0; i < length_seq; i++)
- sumn+=profile[i].get_px();
- if (results) {
- if (EnergyType == Z_DNA)
- cout << "Number of Z-DNA bases = " << sumn << endl;
- else if (EnergyType == Cruciform)
- cout << "Number of cruciform bases = " << sumn << endl;
- else
- cout << "Number of melted bases = " << sumn << endl;
- }
-}
-
-void SIDD_Base::get_column_header()
-{
- if(EnergyType==Z_DNA || EnergyType == Cruciform){
- if(get_showbase()) {
- cout << "Position" << "\tBase" << "\tP(x)"<
-#include
-#include
#include
+#include
+#include
#include
+#include
#include
-#include "G_x.h"
const double PI = 3.1415926535897932;
const int MaxInitialWindowSize = 250;
@@ -32,137 +32,137 @@ const int Len_Gap_Seq = 50;
// V: I Added this
const double scaling_factor = 300.0;
+enum Energetics { Copolymeric, Near_Neighbor, Z_DNA, Cruciform };
+enum Molecule { Circular, Linear };
-enum Energetics{Copolymeric, Near_Neighbor, Z_DNA, Cruciform};
-enum Molecule{Circular, Linear};
-
-class SIDD_Base
-{
+class SIDD_Base {
protected:
- int showbase; //include base pairs in output
- int* plasmid_seq; // hold encoded sequence
- double** en_cruciforms; //array of cruciform pos,length,energies
- int length_seq;
- double Delta_G[4][4]; // free energy
- double Temperature;
- double Salt_Conc; // salt concentration
- double R; // constant
- double C; // torsional stiffness
- double a; // initial energy
- double A; // bases per turn
- double Az;
- double tz;
- double stresslevel;
- double supdensity;
- double alpha; // linking difference
- double theta; // threshold
- double K0;
- double K;
- double RT;
- double TMAT;
- double TMGC;
- double BAT;
- double BGC;
- double Ecr;
- double Zmin;
- int MinWindowSize;
- int MaxWindowSize;
- int results;
- std::string cr_string;
- double min_E;
- double max_E;
- int min_WS;
- char* sa;
- double e1;
- double e2;
- double zz;
- Energetics EnergyType;
- Molecule MoleculeType;
- // profile information
- G_x* profile; // hold profile
- G_x* promatrix[MaxInitialWindowSize+1];
- double minGres;
- double Gres[MaxInitialWindowSize+1][4];
- double* pea; // position-wide energy assignment
- bool Flag_PEA;
- long totalFreq;
- bool write_profile;
- double ZsumG; // total free energy
- double ZsumB; // value by partition function Z
- double sum_1RG; // summation of all states within one run
- double sum_1RB; // summation of all Boltzman factors within one run
- double sum_2RG; // summation of all states within two runs
- double sum_2RB; // summation of all Boltzman factors within two runs
- double sum_3RG;
- double sum_3RB; // summation of all Boltzman factors within three runs
- double sum_4RG;
- double sum_4RB; // summation of all Boltzman factors within three runs
- double Prob_1R;
- double Prob_2R;
- double Prob_3R;
- double maxGx;
- double minGx;
+ int showbase; // include base pairs in output
+ int *plasmid_seq; // hold encoded sequence
+ double **en_cruciforms; // array of cruciform pos,length,energies
+ int length_seq;
+ double Delta_G[4][4]; // free energy
+ double Temperature;
+ double Salt_Conc; // salt concentration
+ double R; // constant
+ double C; // torsional stiffness
+ double a; // initial energy
+ double A; // bases per turn
+ double Az;
+ double tz;
+ double stresslevel;
+ double supdensity;
+ double alpha; // linking difference
+ double theta; // threshold
+ double K0;
+ double K;
+ double RT;
+ double TMAT;
+ double TMGC;
+ double BAT;
+ double BGC;
+ double Ecr;
+ double Zmin;
+ int MinWindowSize;
+ int MaxWindowSize;
+ int results;
+ std::string cr_string;
+ double min_E;
+ double max_E;
+ int min_WS;
+ char *sa;
+ double e1;
+ double e2;
+ double zz;
+ Energetics EnergyType;
+ Molecule MoleculeType;
+ // profile information
+ G_x *profile; // hold profile
+ G_x *promatrix[MaxInitialWindowSize + 1];
+ double minGres;
+ double Gres[MaxInitialWindowSize + 1][4];
+ double *pea; // position-wide energy assignment
+ bool Flag_PEA;
+ long totalFreq;
+ bool write_profile;
+ double ZsumG; // total free energy
+ double ZsumB; // value by partition function Z
+ double sum_1RG; // summation of all states within one run
+ double sum_1RB; // summation of all Boltzman factors within one run
+ double sum_2RG; // summation of all states within two runs
+ double sum_2RB; // summation of all Boltzman factors within two runs
+ double sum_3RG;
+ double sum_3RB; // summation of all Boltzman factors within three runs
+ double sum_4RG;
+ double sum_4RB; // summation of all Boltzman factors within three runs
+ double Prob_1R;
+ double Prob_2R;
+ double Prob_3R;
+ double maxGx;
+ double minGx;
private:
- void set_K();
- void set_Delta_G();
- void set_Alpha();
- void set_junction();
- void set_E_limit();
- void set_minGres();
- void set_MinWindowSize();
- void set_MaxWindowSize();
+ void set_K();
+ void set_Delta_G();
+ void set_Alpha();
+ void set_junction();
+ void set_E_limit();
+ void set_minGres();
+ void set_MinWindowSize();
+ void set_MaxWindowSize();
public:
- SIDD_Base();
- virtual ~SIDD_Base();
- bool initializer(std::string filename,int len);
- void set_EnergyType(Energetics);
- void set_temperature(double temperature);
- void set_showres(int showres);
- void set_MoleculeType(Molecule);
- void set_Cinitiation();
- void show_seq();
- void show_deltaG();
- void show_parameter();
- double get_Delta_G(int i, int j){return (i< 4 && j < 4)?Delta_G[i][j]:0.0;};
- double calc_Gres(int,int);
- double calc_NI(int);
- void calc_Z(int);
- double sum_WindowNI(int startp, int n);
- double sum_WindowZ(int startp, int n);
- void count_AT_GC(int startp, int n, int& c_AT, int& c_GC, double& tempE);
- double sum_WindowNIbyBases(int startp, int n);
- double sum_cruciform(int startp, int n);
- double calc_OPenBasesEnergy(int startp, int n);
- bool update_promatrix(int startp, int n, double x, double bzfactor);
- void reset_promatrix();
- void fill_profile();
- void reset_profile();
- void get_column_header();
- void show_profile();
- void calc_profile();
- void sum_open();
- void show_avg_runs();
- void write_close(){write_profile = false;}
- void write_open(){write_profile = true;}
- void set_Flag_PEA(bool f){Flag_PEA = f;}
- bool get_Flag_PEA(){return Flag_PEA;}
- void assign_pea(char*);
- double prepare_sequence(std::string dnasequence);
- void set_cruciform_string(std::string cruciform_string);
- void prepare_cruciforms();
+ SIDD_Base();
+ virtual ~SIDD_Base();
+ bool initializer(std::string filename, int len);
+ void set_EnergyType(Energetics);
+ void set_temperature(double temperature);
+ void set_showres(int showres);
+ void set_MoleculeType(Molecule);
+ void set_Cinitiation();
+ void show_seq();
+ void show_deltaG();
+ void show_parameter();
+ double get_Delta_G(int i, int j) {
+ return (i < 4 && j < 4) ? Delta_G[i][j] : 0.0;
+ };
+ double calc_Gres(int, int);
+ double calc_NI(int);
+ void calc_Z(int);
+ double sum_WindowNI(int startp, int n);
+ double sum_WindowZ(int startp, int n);
+ void count_AT_GC(int startp, int n, int &c_AT, int &c_GC, double &tempE);
+ double sum_WindowNIbyBases(int startp, int n);
+ double sum_cruciform(int startp, int n);
+ double calc_OPenBasesEnergy(int startp, int n);
+ bool update_promatrix(int startp, int n, double x, double bzfactor);
+ void reset_promatrix();
+ void fill_profile();
+ void reset_profile();
+ void get_column_header();
+ void show_profile();
+ void calc_profile();
+ void sum_open();
+ void show_avg_runs();
+ void write_close() { write_profile = false; }
+ void write_open() { write_profile = true; }
+ void set_Flag_PEA(bool f) { Flag_PEA = f; }
+ bool get_Flag_PEA() { return Flag_PEA; }
+ void assign_pea(char *);
+ double prepare_sequence(std::string dnasequence);
+ void set_cruciform_string(std::string cruciform_string);
+ void prepare_cruciforms();
- void set_threshold(double threshold);
- void set_salt_conc(double salt_conc);
- void set_min_e(double min_e);
- void set_stress_level(double stress_level);
- double get_zsumb(){return ZsumB;}
- int get_sequence_length();
- void set_sequence_length(int len);
- char decode_base(int base);
- void set_showbase(int show) {showbase = show;}
- int get_showbase() {return showbase;}
+ void set_threshold(double threshold);
+ void set_salt_conc(double salt_conc);
+ void set_min_e(double min_e);
+ void set_stress_level(double stress_level);
+ double get_zsumb() { return ZsumB; }
+ int get_sequence_length();
+ void set_sequence_length(int len);
+ char decode_base(int base);
+ void set_showbase(int show) { showbase = show; }
+ int get_showbase() { return showbase; }
};
-#endif
+#endif
diff --git a/src/trans_three/qsidd.cpp b/src/trans_three/qsidd.cpp
index 3f3be2b..0924def 100644
--- a/src/trans_three/qsidd.cpp
+++ b/src/trans_three/qsidd.cpp
@@ -1,6 +1,6 @@
// qsidd.cpp: main() function
// program to implement algorithm developed by Craig Benham
-//
+//
// author: Chengpeng Bi
// modifiers: Dina Zhabinskaya, Sally Madden, Ian Korf
// compiler: g++
@@ -12,12 +12,12 @@
//////////////////////////////////////////////////////////////////////////////
#include "SIDD_4R.h"
-#include
#include
-#include
#include
-#include
#include
+#include
+#include
+#include
static char help[] = "\
QSIDD Help\n\
@@ -51,7 +51,6 @@ Position P(x) G(x)\n\
5 7.51455e-07 11.4479\n\
";
-
static char usage[] = "\
usage: qsidd [options] \n\
options:\n\
@@ -73,144 +72,181 @@ options:\n\
-r print ensemble average results\n\
";
-int main(int argc, char* argv[])
-{
- SIDD_4R sidd;
- time_t time_1, time_2;
- int c;
- extern int optind;
-
- // default parameters
- Energetics et = Copolymeric;
- Molecule mt = Linear;
- char *energy = NULL;
- char *dnafile = NULL;
-
- stringstream dnasequence;
- stringstream cruciform_string;
+int main(int argc, char *argv[]) {
+ SIDD_4R sidd;
+ time_t time_1, time_2;
+ int c;
+ extern int optind;
+
+ // default parameters
+ Energetics et = Copolymeric;
+ Molecule mt = Linear;
+ char *energy = NULL;
+ char *dnafile = NULL;
- int verbose = 0;
- int showpar = 0;
- int showres = 0;
- int showbase = 0;
- double salt_conc = 0.01;
- double temperature = 310.00;
- double stress_level=0.06;
- double threshold = 12;
- int usefile = 0;
- // option processing
- while ((c = getopt(argc, argv, "hcfnbvpraietTZCXsm:")) != -1) {
- switch (c) {
- case 'c': mt = Circular; break;
- case 'n': et = Near_Neighbor; break;
- case 'b': showbase = 1; break;
- case 'f': usefile = 1; break;
- case 'Z': et = Z_DNA; break;
- case 'C': et = Cruciform; break;
- case 'X': cruciform_string << argv[optind++]; break;
- case 'e': energy = argv[optind++]; break;
- case 'v': verbose = 1; break;
- case 'p': showpar = 1; break;
- case 'r': showres = 1; break;
- case 'T': temperature = atof(argv[optind++]); break;
- case 's': stress_level = atof(argv[optind++]); break;
- case 'i': salt_conc=atof(argv[optind++]); break;
- case 't': threshold=atof(argv[optind++]); break;
- case 'h': cout< 15)
- cout << "WARNING: threshold is too high, execution time may be very long"<< endl;
- if (stress_level > 0.15 || stress_level < -0.15)
- cout << "WARNING: superhelical density is outside of physiological range"<< endl;
- if (temperature < 220 || temperature > 320)
- cout << "WARNING: temperature is outside of physiological range"<< endl;
- if (salt_conc < 0.0001)
- cout << "WARNING: salt concentration is outside of physiological range"<< endl;
- if (length < 1500)
- cout << "WARNING: sequence length is too short"<< endl;
- if (length > 10000)
- cout << "WARNING: sequence length is too long"<< endl;
+ if (!sidd.initializer(dnasequence.str().c_str(), length)) {
+ cerr << "sidd.initializer failed" << endl;
+ exit(1);
+ }
+ if (sidd.get_Flag_PEA())
+ sidd.assign_pea(energy);
-// sidd.show_seq(); // show sequence
- sidd.gen_OpenBaseEnergy(); // SIDD_1R (calculate opening energies and sort w/ increasing energy for each window size 1 to 250)
- sidd.write_close(); //don't write profile yet, just find minE
- sidd.Search_Low1RE(); // this step is here to update minE (from the zero-run state) if it's found in one-run states
- sidd.write_open(); // start to store info
- if(verbose)
- cout << "writing profile...\n";
- if (showpar)
- sidd.show_parameter();
- //cout << "AAAA1"<< endl;
- sidd.reset_profile(); // initialize profile
- //cout << "AAAA2"<< endl;
- sidd.Search_Low1RE(); // searching for one-run states and now storing info
- //cout << "AAAA3"<< endl;
- sidd.Search_Low2RE(); // searching states for two-run
- //cout << "AAAA4"<< endl;
- sidd.Search_Low3RE(); // searching states for three-run
- //cout << "AAAA5"<< endl;
- sidd.Search_Low4RE(); // searching states for four-run
- //cout << "AAAA6"<< endl;
- sidd.fill_profile(); // computing profile
- //cout << "AAAA7"<< endl;
- sidd.calc_profile(); // calculate profile
- time_2 = time(NULL); //end time
- sidd.sum_open();
- time_2 = time(NULL); //end time
- if (showpar)
- cout << "Run time = " << (time_2-time_1) << " sec" << endl;
- sidd.show_profile(); // send output to screen or a disk file
- return 0; // end of program
+ if (threshold < 9)
+ cout << "WARNING: threshold is too small, results may be inaccurate"
+ << endl;
+ if (threshold > 15)
+ cout << "WARNING: threshold is too high, execution time may be very long"
+ << endl;
+ if (stress_level > 0.15 || stress_level < -0.15)
+ cout << "WARNING: superhelical density is outside of physiological range"
+ << endl;
+ if (temperature < 220 || temperature > 320)
+ cout << "WARNING: temperature is outside of physiological range" << endl;
+ if (salt_conc < 0.0001)
+ cout << "WARNING: salt concentration is outside of physiological range"
+ << endl;
+ if (length < 1500)
+ cout << "WARNING: sequence length is too short" << endl;
+ if (length > 10000)
+ cout << "WARNING: sequence length is too long" << endl;
+ // sidd.show_seq(); // show sequence
+ sidd.gen_OpenBaseEnergy(); // SIDD_1R (calculate opening energies and sort w/
+ // increasing energy for each window size 1 to 250)
+ sidd.write_close(); // don't write profile yet, just find minE
+ sidd.Search_Low1RE(); // this step is here to update minE (from the zero-run
+ // state) if it's found in one-run states
+ sidd.write_open(); // start to store info
+ if (verbose)
+ cout << "writing profile...\n";
+ if (showpar)
+ sidd.show_parameter();
+ // cout << "AAAA1"<< endl;
+ sidd.reset_profile(); // initialize profile
+ // cout << "AAAA2"<< endl;
+ sidd.Search_Low1RE(); // searching for one-run states and now storing info
+ // cout << "AAAA3"<< endl;
+ sidd.Search_Low2RE(); // searching states for two-run
+ // cout << "AAAA4"<< endl;
+ sidd.Search_Low3RE(); // searching states for three-run
+ // cout << "AAAA5"<< endl;
+ sidd.Search_Low4RE(); // searching states for four-run
+ // cout << "AAAA6"<< endl;
+ sidd.fill_profile(); // computing profile
+ // cout << "AAAA7"<< endl;
+ sidd.calc_profile(); // calculate profile
+ time_2 = time(NULL); // end time
+ sidd.sum_open();
+ time_2 = time(NULL); // end time
+ if (showpar)
+ cout << "Run time = " << (time_2 - time_1) << " sec" << endl;
+ sidd.show_profile(); // send output to screen or a disk file
+ return 0; // end of program
}
diff --git a/src/trans_three/stat_1R.cpp b/src/trans_three/stat_1R.cpp
index 87520f3..c8f29cd 100644
--- a/src/trans_three/stat_1R.cpp
+++ b/src/trans_three/stat_1R.cpp
@@ -21,38 +21,28 @@
// Construction/Destruction
//////////////////////////////////////////////////////////////////////
-stat_1R::stat_1R(int s, double e, double RT) //constructor
+stat_1R::stat_1R(int s, double e, double RT) // constructor
{
- start_pos1 = s; //starting position
- energy = e;
+ start_pos1 = s; // starting position
+ energy = e;
}
-stat_1R::~stat_1R() //destructor ???
-{
-
-}
+stat_1R::~stat_1R() // destructor ???
+{}
+bool stat_1R::operator<(const stat_1R &x1) const {
-
-bool stat_1R::operator<(const stat_1R& x1)const
-{
-
- return (energy < x1.get_energy());
-
+ return (energy < x1.get_energy());
}
-bool stat_1R::operator==(const stat_1R& x1)const
-{
-
- return (energy == x1.get_energy());
+bool stat_1R::operator==(const stat_1R &x1) const {
+ return (energy == x1.get_energy());
}
-bool stat_1R::operator>(const stat_1R& x1)const
-{
-
- return (energy > x1.get_energy());
+bool stat_1R::operator>(const stat_1R &x1) const {
+ return (energy > x1.get_energy());
}
#endif
diff --git a/src/trans_three/stat_1R.h b/src/trans_three/stat_1R.h
index 952b2d1..27ee50a 100644
--- a/src/trans_three/stat_1R.h
+++ b/src/trans_three/stat_1R.h
@@ -16,19 +16,19 @@
#ifndef STAT_1R_H_
#define STAT_1R_H_
-class stat_1R
-{
+class stat_1R {
protected:
- int start_pos1;
- double energy;
+ int start_pos1;
+ double energy;
+
public:
- stat_1R(int start_position, double energy, double rt);
- virtual ~stat_1R();
- int get_pos1() const{return start_pos1;}
- bool operator <(const stat_1R&) const;
- bool operator ==(const stat_1R&)const;
- bool operator >(const stat_1R&) const;
- double get_energy()const {return energy;}
+ stat_1R(int start_position, double energy, double rt);
+ virtual ~stat_1R();
+ int get_pos1() const { return start_pos1; }
+ bool operator<(const stat_1R &) const;
+ bool operator==(const stat_1R &) const;
+ bool operator>(const stat_1R &) const;
+ double get_energy() const { return energy; }
};
#endif // !defined(AFX_STAT_1R_H__ACD4635F_79AA_463A_A66C_6B1CA4476EE0__INCLUDED_)
diff --git a/tests/regression/data/pbr322.toy.fa b/tests/regression/data/pbr322.toy.fa
index 23324fb..ab43845 100644
--- a/tests/regression/data/pbr322.toy.fa
+++ b/tests/regression/data/pbr322.toy.fa
@@ -61,4 +61,3 @@ ACACGGAAATGTTGAATACTCATACTCTTCCTTTTTCAATATTATTGAAGCATTTATCAGGGTTATTGTC
TCATGAGCGGATACATATTTGAATGTATTTAGAAAAATAAACAAATAGGGGTTCCGCGCACATTTCCCCG
AAAAGTGCCACCTGACGTCTAAGAAACCATTATTATCATGACATTAACCTATAAAAATAGGCGTATCACG
AGGCCCTTTCGTCTTCAAGAA
-
diff --git a/tests/regression/test_cli.py b/tests/regression/test_cli.py
index 1c92844..b745d11 100644
--- a/tests/regression/test_cli.py
+++ b/tests/regression/test_cli.py
@@ -1,10 +1,8 @@
from __future__ import annotations
import pytest
-
from conftest import SistRun
-
pytestmark = pytest.mark.regression
@@ -26,13 +24,9 @@ def assert_output_is_created(sist_run: SistRun) -> None:
output_path = sist_run.output_path
- assert output_path.is_file(), (
- f"Expected output file was not created: {output_path}"
- )
+ assert output_path.is_file(), f"Expected output file was not created: {output_path}"
- assert output_path.stat().st_size > 0, (
- f"Output file is empty: {output_path}"
- )
+ assert output_path.stat().st_size > 0, f"Output file is empty: {output_path}"
def test_competition_command_succeeds(
@@ -71,9 +65,7 @@ def test_competition_output_contains_expected_sections(
)
for section in expected_sections:
- assert section in output, (
- f"Expected output section was not found: {section!r}"
- )
+ assert section in output, f"Expected output section was not found: {section!r}"
def test_transition_command_succeeds(
@@ -109,6 +101,5 @@ def test_transition_output_contains_expected_sections(
for section in expected_sections:
assert section in output, (
- f"{transition_run.name}: expected output section "
- f"was not found: {section!r}"
+ f"{transition_run.name}: expected output section was not found: {section!r}"
)
diff --git a/tests/regression/test_regression.py b/tests/regression/test_regression.py
index 29d1e0a..c146d10 100644
--- a/tests/regression/test_regression.py
+++ b/tests/regression/test_regression.py
@@ -4,10 +4,8 @@
from pathlib import Path
import pytest
-
from conftest import SistRun
-
pytestmark = pytest.mark.regression
@@ -15,10 +13,7 @@
REFERENCE_DIRECTORY = Path(__file__).resolve().parent / "reference" / REFERENCE_VERSION
-COMPETITION_REFERENCE = (
- REFERENCE_DIRECTORY
- / "competition.rebuilt.txt"
-)
+COMPETITION_REFERENCE = REFERENCE_DIRECTORY / "competition.rebuilt.txt"
# The baselines reproduced all shared printed values exactly.
RELATIVE_TOLERANCE = 0.0
@@ -144,9 +139,7 @@ def parse_competition_output(path: Path) -> ParsedCompetitionOutput:
continue
if line.startswith("Scaling factor "):
- value = parse_numeric_value(
- line.removeprefix("Scaling factor ")
- )
+ value = parse_numeric_value(line.removeprefix("Scaling factor "))
if value is None:
raise AssertionError(
@@ -244,11 +237,7 @@ def parse_transition_output(path: Path) -> ParsedTransitionOutput:
base = columns[1]
probability = float(columns[2])
- energy = (
- float(columns[3])
- if profile_has_energy
- else None
- )
+ energy = float(columns[3]) if profile_has_energy else None
profile[position] = TransitionProfileRow(
base=base,
@@ -259,9 +248,7 @@ def parse_transition_output(path: Path) -> ParsedTransitionOutput:
continue
if line.startswith("Scaling factor "):
- value = parse_numeric_value(
- line.removeprefix("Scaling factor ")
- )
+ value = parse_numeric_value(line.removeprefix("Scaling factor "))
if value is None:
raise AssertionError(
@@ -312,10 +299,7 @@ def assert_number_matches(
expected,
rel=RELATIVE_TOLERANCE,
abs=ABSOLUTE_TOLERANCE,
- ), (
- f"{name} changed: expected {expected!r}, "
- f"actual {actual!r}"
- )
+ ), f"{name} changed: expected {expected!r}, actual {actual!r}"
def assert_metadata_matches(
@@ -430,10 +414,7 @@ def test_transition_scientific_results_match_baseline(
f"stderr:\n{process.stderr}"
)
- reference_output = (
- REFERENCE_DIRECTORY
- / f"{transition_run.name}.txt"
- )
+ reference_output = REFERENCE_DIRECTORY / f"{transition_run.name}.txt"
assert reference_output.is_file(), (
f"Reference output does not exist: {reference_output}"
@@ -473,11 +454,8 @@ def test_transition_scientific_results_match_baseline(
actual=actual_row.probability,
)
- assert (actual_row.energy is None) == (
- expected_row.energy is None
- ), (
- f"{transition_run.name} energy output changed at "
- f"position {position}"
+ assert (actual_row.energy is None) == (expected_row.energy is None), (
+ f"{transition_run.name} energy output changed at position {position}"
)
if expected_row.energy is not None: