primer3 0.1.0

//! Thermodynamic analysis of DNA sequences.
//!
//! Provides functions for calculating hairpin, homodimer, heterodimer, and
//! 3' end stability thermodynamics using the nearest-neighbor model.

use std::ffi::{CStr, c_void};
use std::fmt;

use crate::conditions::SolutionConditions;
use crate::error::{Primer3Error, Result};
use crate::init::ensure_initialized;

/// Result of a thermodynamic calculation (hairpin, dimer, or end stability).
#[derive(Debug, Clone)]
#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))]
pub struct ThermoResult {
    /// Whether a thermodynamically stable structure was found.
    structure_found: bool,
    /// Melting temperature in degrees Celsius.
    tm: f64,
    /// Delta G (Gibbs free energy) in cal/mol.
    dg: f64,
    /// Delta H (enthalpy) in cal/mol.
    dh: f64,
    /// Delta S (entropy) in cal/(K*mol).
    ds: f64,
    /// ASCII representation of the secondary structure, if requested.
    ascii_structure: Option<String>,
}

impl ThermoResult {
    /// Whether a thermodynamically stable structure was found.
    pub fn structure_found(&self) -> bool {
        self.structure_found
    }

    /// Melting temperature in degrees Celsius.
    pub fn tm(&self) -> f64 {
        self.tm
    }

    /// Delta G (Gibbs free energy) in cal/mol.
    pub fn dg(&self) -> f64 {
        self.dg
    }

    /// Delta H (enthalpy) in cal/mol.
    pub fn dh(&self) -> f64 {
        self.dh
    }

    /// Delta S (entropy) in cal/(K*mol).
    pub fn ds(&self) -> f64 {
        self.ds
    }

    /// ASCII representation of the secondary structure.
    ///
    /// Only available if `output_structure` was set to `true` in the
    /// calculation parameters.
    pub fn ascii_structure(&self) -> Option<&str> {
        self.ascii_structure.as_deref()
    }

    /// Compares two `ThermoResult` values for approximate equality.
    ///
    /// Returns `true` if all `f64` fields are within `epsilon` of each other
    /// and all non-float fields are equal.
    pub fn approx_eq(&self, other: &Self, epsilon: f64) -> bool {
        self.structure_found == other.structure_found
            && (self.tm - other.tm).abs() < epsilon
            && (self.dg - other.dg).abs() < epsilon
            && (self.dh - other.dh).abs() < epsilon
            && (self.ds - other.ds).abs() < epsilon
            && self.ascii_structure == other.ascii_structure
    }
}

impl fmt::Display for ThermoResult {
    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
        write!(
            f,
            "ThermoResult(structure_found={}, tm={:.2}, dg={:.0}, dh={:.0}, ds={:.1})",
            self.structure_found, self.tm, self.dg, self.dh, self.ds,
        )
    }
}

/// Parameters for thermodynamic calculations (hairpin, dimer, end stability).
///
/// Defaults match primer3-py's standard PCR conditions.
///
/// Note: The C library's `dimer` flag is handled automatically: enabled for
/// homodimer/heterodimer calculations, disabled for hairpin.
#[derive(Debug, Clone)]
#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))]
pub struct ThermoArgs {
    /// Salt and buffer concentrations.
    pub conditions: SolutionConditions,
    /// Simulation temperature in Celsius for dG calculation (default: 37.0).
    pub temp_c: f64,
    /// Maximum loop size for hairpin/dimer calculations (default: 30).
    pub max_loop: usize,
    /// Whether to compute and return the ASCII structure representation.
    pub output_structure: bool,
}

impl Default for ThermoArgs {
    fn default() -> Self {
        Self {
            conditions: SolutionConditions::default(),
            temp_c: 37.0,
            max_loop: 30,
            output_structure: false,
        }
    }
}

/// Maximum sequence length for the short oligo in thermodynamic alignment.
/// One sequence must be <= this length for the two-state NN model.
const THAL_MAX_ALIGN: usize = 60;

/// Maximum sequence length for the long sequence in thermodynamic alignment.
/// Used for template mispriming checks.
const THAL_MAX_SEQ: usize = 10_000;

/// Calls the C `thal()` function and converts the result to a `ThermoResult`.
fn call_thal(
    seq1: &str,
    seq2: &str,
    alignment_type: u32,
    is_dimer: bool,
    args: &ThermoArgs,
) -> Result<ThermoResult> {
    let mode = if args.output_structure {
        primer3_sys::thal_mode_THL_STRUCT
    } else {
        primer3_sys::thal_mode_THL_FAST
    };

    let thal_args = primer3_sys::thal_args {
        type_: alignment_type,
        maxLoop: args.max_loop as std::os::raw::c_int,
        mv: args.conditions.mv_conc,
        dv: args.conditions.dv_conc,
        dntp: args.conditions.dntp_conc,
        dna_conc: args.conditions.dna_conc,
        temp: args.temp_c + 273.15, // thal() expects Kelvin
        dimer: i32::from(is_dimer),
    };

    // Null-terminated uppercase byte sequences for C
    let mut s1_c: Vec<u8> = seq1.bytes().map(|b| b.to_ascii_uppercase()).collect();
    s1_c.push(0);
    let mut s2_c: Vec<u8> = seq2.bytes().map(|b| b.to_ascii_uppercase()).collect();
    s2_c.push(0);

    let mut results: primer3_sys::thal_results = unsafe { std::mem::zeroed() };

    unsafe {
        primer3_sys::thal(s1_c.as_ptr(), s2_c.as_ptr(), &thal_args, mode, &mut results);
    }

    // Check for error message
    if results.msg[0] != 0 {
        let msg = unsafe { CStr::from_ptr(results.msg.as_ptr()) }.to_string_lossy().into_owned();
        // Clean up sec_struct if allocated
        if !results.sec_struct.is_null() {
            unsafe { libc::free(results.sec_struct.cast::<c_void>()) };
        }
        return Err(Primer3Error::Library(msg.into_boxed_str()));
    }

    let ascii_structure = if results.sec_struct.is_null() {
        None
    } else {
        let s = unsafe { CStr::from_ptr(results.sec_struct) }.to_string_lossy().into_owned();
        unsafe { libc::free(results.sec_struct.cast::<c_void>()) };
        Some(s)
    };

    // thal() returns temp=0.0 when no stable structure is found
    let structure_found = results.temp > 0.0;

    Ok(ThermoResult {
        structure_found,
        tm: results.temp,
        dg: results.dg,
        dh: results.dh,
        ds: results.ds,
        ascii_structure,
    })
}

fn validate_single_seq(seq: &str, operation: &'static str) -> Result<()> {
    if seq.is_empty() {
        return Err(Primer3Error::InvalidSequence("sequence is empty".into()));
    }
    if seq.len() > THAL_MAX_ALIGN {
        return Err(Primer3Error::SequenceTooLong {
            length: seq.len(),
            max: THAL_MAX_ALIGN,
            operation,
        });
    }
    Ok(())
}

fn validate_pair_seqs(seq1: &str, seq2: &str, operation: &'static str) -> Result<()> {
    if seq1.is_empty() || seq2.is_empty() {
        return Err(Primer3Error::InvalidSequence("sequence is empty".into()));
    }
    if seq1.len() > THAL_MAX_ALIGN && seq2.len() > THAL_MAX_ALIGN {
        return Err(Primer3Error::InvalidSequence(
            format!("at least one sequence must be <= {THAL_MAX_ALIGN} bp for {operation}")
                .into_boxed_str(),
        ));
    }
    let longer = seq1.len().max(seq2.len());
    if longer > THAL_MAX_SEQ {
        return Err(Primer3Error::SequenceTooLong { length: longer, max: THAL_MAX_SEQ, operation });
    }
    Ok(())
}

/// Calculates hairpin formation thermodynamics for a DNA sequence.
///
/// The maximum sequence length is 60 bp (the limit for reliable two-state
/// nearest-neighbor models).
///
/// # Errors
///
/// Returns an error if the sequence is empty, too long, or invalid.
///
/// # Example
///
/// ```no_run
/// let result = primer3::calc_hairpin("CCCCCATCCGATCAGGGGG").unwrap();
/// println!("Hairpin Tm = {:.1}°C", result.tm());
/// ```
pub fn calc_hairpin(seq: &str) -> Result<ThermoResult> {
    calc_hairpin_with(seq, &ThermoArgs::default())
}

/// Calculates hairpin formation thermodynamics with custom parameters.
///
/// # Errors
///
/// Returns an error if the sequence is empty, too long, or the C library reports an error.
pub fn calc_hairpin_with(seq: &str, args: &ThermoArgs) -> Result<ThermoResult> {
    ensure_initialized()?;
    validate_single_seq(seq, "hairpin")?;

    call_thal(seq, seq, primer3_sys::thal_alignment_type_thal_hairpin, false, args)
}

/// Calculates homodimer formation thermodynamics for a DNA sequence.
///
/// The maximum sequence length is 60 bp.
///
/// # Errors
///
/// Returns an error if the sequence is empty, too long, or invalid.
///
/// # Example
///
/// ```no_run
/// let result = primer3::calc_homodimer("AGTCTAGTCTATCGATCG").unwrap();
/// println!("Homodimer dG = {:.0} cal/mol", result.dg());
/// ```
pub fn calc_homodimer(seq: &str) -> Result<ThermoResult> {
    calc_homodimer_with(seq, &ThermoArgs::default())
}

/// Calculates homodimer formation thermodynamics with custom parameters.
///
/// # Errors
///
/// Returns an error if the sequence is empty, too long, or the C library reports an error.
pub fn calc_homodimer_with(seq: &str, args: &ThermoArgs) -> Result<ThermoResult> {
    ensure_initialized()?;
    validate_single_seq(seq, "homodimer")?;

    call_thal(seq, seq, primer3_sys::thal_alignment_type_thal_any, true, args)
}

/// Calculates heterodimer formation thermodynamics between two DNA sequences.
///
/// At least one sequence must be <= 60 bp. The other may be up to 10,000 bp
/// (useful for checking primer-template mispriming).
///
/// # Errors
///
/// Returns an error if sequences are empty, both too long, or invalid.
///
/// # Example
///
/// ```no_run
/// let result = primer3::calc_heterodimer("AAAAAAAAAA", "TTTTTTTTTT").unwrap();
/// println!("Heterodimer Tm = {:.1}°C", result.tm());
/// ```
pub fn calc_heterodimer(seq1: &str, seq2: &str) -> Result<ThermoResult> {
    calc_heterodimer_with(seq1, seq2, &ThermoArgs::default())
}

/// Calculates heterodimer formation thermodynamics with custom parameters.
///
/// # Errors
///
/// Returns an error if either sequence is empty, too long, or the C library reports an error.
pub fn calc_heterodimer_with(seq1: &str, seq2: &str, args: &ThermoArgs) -> Result<ThermoResult> {
    ensure_initialized()?;
    validate_pair_seqs(seq1, seq2, "heterodimer")?;

    call_thal(seq1, seq2, primer3_sys::thal_alignment_type_thal_any, true, args)
}

/// Calculates 3' end stability of `seq1` hybridized against `seq2`.
///
/// Uses `thal_end1` alignment type: aligns seq1's 3' end against seq2.
/// Typically seq1 is the primer and seq2 is the template.
///
/// At least one sequence must be <= 60 bp.
///
/// # Errors
///
/// Returns an error if sequences are empty, both too long, or invalid.
///
/// # Example
///
/// ```no_run
/// let result = primer3::calc_end_stability(
///     "AGCTATTTTTTTTTTT",
///     "CATGATTTTTTTTTTTTTT",
/// ).unwrap();
/// println!("End stability dG = {:.0} cal/mol", result.dg());
/// ```
pub fn calc_end_stability(seq1: &str, seq2: &str) -> Result<ThermoResult> {
    calc_end_stability_with(seq1, seq2, &ThermoArgs::default())
}

/// Calculates 3' end stability with custom parameters.
///
/// Uses `thal_end1` alignment type: aligns seq1's 3' end against seq2.
/// Typically seq1 is the primer and seq2 is the template.
///
/// # Errors
///
/// Returns an error if either sequence is empty, too long, or the C library reports an error.
pub fn calc_end_stability_with(seq1: &str, seq2: &str, args: &ThermoArgs) -> Result<ThermoResult> {
    ensure_initialized()?;
    validate_pair_seqs(seq1, seq2, "end_stability")?;

    call_thal(seq1, seq2, primer3_sys::thal_alignment_type_thal_end1, true, args)
}