lbfgsbrs 0.1.2

Rust port of L-BFGS-B-C
Documentation

use log::{warn, info, debug, trace};

pub type ftnlen = i32;

pub const ERROR: i32 = 200 as i32;
pub const ERROR_END: i32 = 240 as i32;
pub const ERROR_N0: i32 = 209 as i32;
pub const ERROR_M0: i32 = 210 as i32;
pub const ERROR_FACTR: i32 = 211 as i32;
pub const ERROR_NBD: i32 = 212 as i32;
pub const ERROR_FEAS: i32 = 213 as i32;
pub const WORD_DEFAULT: i32 = 0 as i32;
pub const WORD_CON: i32 = 1 as i32;
pub const WORD_BND: i32 = 2 as i32;
pub const WORD_TNT: i32 = 3 as i32;

// Safe implementation of prn1lb
pub fn prn1lb(
    n: i32,
    m: i32,
    l: &[f64],
    u: &[f64],
    x: &[f64],
    epsmch: f64,
) -> i32 {
    info!("           * * *");
    info!("        RUNNING THE L-BFGS-B CODE");
    info!("           * * *");
    info!("Machine precision = {:.2e}", epsmch);
    info!(" N = {:10}   M = {:10}", n, m);
    
    let mut l_values = String::new();
    for i in 0..n as usize {
        l_values.push_str(&format!("{:.2e} ", l[i]));
    }
    trace!("L  = {}", l_values);
    
    let mut x0_values = String::new();
    for i in 0..n as usize {
        x0_values.push_str(&format!("{:.2e} ", x[i]));
    }
    trace!("X0 = {}", x0_values);
    
    let mut u_values = String::new();
    for i in 0..n as usize {
        u_values.push_str(&format!("{:.2e} ", u[i]));
    }
    trace!("U  = {}", u_values);

    0
}

// Safe implementation of prn2lb
pub fn prn2lb(
    n: i32,
    x: &[f64],
    f: f64,
    g: &[f64],
    iter: i32,
    sbgnrm: f64,
    word: &mut i32,
    iword: i32,
    iback: i32,
    xstep: f64,
) -> i32 {
    if iword == 0 {
        *word = WORD_CON as i32;
    } else if iword == 1 {
        *word = WORD_BND as i32;
    } else if iword == 5 {
        *word = WORD_TNT as i32;
    } else {
        *word = WORD_DEFAULT as i32;
    }

    debug!("LINE SEARCH {} times; norm of step = {:.2e}", iback, xstep);
    debug!("At iterate {:5}, f(x)= {:5.2e}, ||proj grad||_infty = {:.2e}", iter, f, sbgnrm);
    
    let mut x_values = String::new();
    for i in 0..n as usize {
        x_values.push_str(&format!("{:.2e} ", x[i]));
    }
    trace!("X = {}", x_values);
    
    let mut g_values = String::new();
    for i in 0..n as usize {
        g_values.push_str(&format!("{:.2e} ", g[i]));
    }
    trace!("G = {}", g_values);

    info!("At iterate {:5}, f(x)= {:5.2e}, ||proj grad||_infty = {:.2e}", iter, f, sbgnrm);

    0
}

// Safe implementation of prn3lb
pub fn prn3lb(
    n: i32,
    x: &[f64],
    f: f64,
    task: &mut i32,
    info: i32,
    iter: i32,
    nfgv: i32,
    nintol: i32,
    nskip: i32,
    nact: i32,
    sbgnrm: f64,
    time: f64,
    k: i32,
    cachyt: f64,
    sbtime: f64,
    lnscht: f64,
) -> i32 {
    let is_error = *task >= ERROR as i32 && *task <= ERROR_END as i32;
    
    if !is_error {
        info!("           * * * ");
        info!("Tit   = total number of iterations");
        info!("Tnf   = total number of function evaluations");
        info!("Tnint = total number of segments explored during Cauchy searches");
        info!("Skip  = number of BFGS updates skipped");
        info!("Nact  = number of active bounds at final generalized Cauchy point");
        info!("Projg = norm of the final projected gradient");
        info!("F     = final function value");
        info!("           * * * ");
        info!("   N    Tit   Tnf  Tnint  Skip  Nact      Projg        F");
        info!("{:5} {:5} {:5} {:5} {:5} {:5}\t{:6.2e} {:9.5e}", 
            n, iter, nfgv, nintol, nskip, nact, sbgnrm, f);
            
        let mut x_values = String::new();
        for i in 0..n as usize {
            x_values.push_str(&format!(" {:.2e}", x[i]));
        }
        trace!("X = {}", x_values);

        info!("F(x) = {:.9e}", f);
    }

    info!("{}", *task);
    if info != 0 {
        if info == -1 {
            warn!("Matrix in 1st Cholesky factorization in formk is not Pos. Def.");
        }
        if info == -2 {
            warn!("Matrix in 2nd Cholesky factorization in formk is not Pos. Def.");
        }
        if info == -3 {
            warn!("Matrix in the Cholesky factorization in formt is not Pos. Def.");
        }
        if info == -4 {
            warn!("Derivative >= 0, backtracking line search impossible.");
            warn!("Previous x, f and g restored.");
            warn!("Possible causes: 1 error in function or gradient evaluation;");
            warn!("                 2 rounding errors dominate computation.");
        }
        if info == -5 {
            warn!("Warning:  more than 10 function and gradient");
            warn!("  evaluations in the last line search.  Termination");
            warn!("  may possibly be caused by a bad search direction.");
        }
        if info == -6 {
            warn!("Input nbd({}) is invalid", k);
        }
        if info == -7 {
            warn!("l({}) > u({}). No feasible solution.", k, k);
        }
        if info == -8 {
            warn!("The triangular system is singular.");
        }
        if info == -9 {
            warn!("Line search cannot locate an adequate point after 20 function");
            warn!("and gradient evaluations.  Previous x, f and g restored.");
            warn!("Possible causes: 1 error in function or gradient evaluation;");
            warn!("                 2 rounding error dominate computation.");
        }
    }

    info!("Cauchy                time {:.3e} seconds.", cachyt);
    info!("Subspace minimization time {:.3e} seconds.", sbtime);
    info!("Line search           time {:.3e} seconds.", lnscht);

    info!("Total User time {:.3e} seconds.", time);
    
    0 as i32
}

// Safe implementation of errclb
pub fn errclb(
    n: i32,
    m: i32,
    factr: f64,
    l: &[f64],
    u: &[f64],
    nbd: &[i32],
    task: &mut i32,
    info: &mut i32,
    k: &mut i32,
) -> i32 {
    if n <= 0 {
        *task = ERROR_N0 as i32;
    }
    if m <= 0 {
        *task = ERROR_M0 as i32;
    }
    if factr < 0.0f64 {
        *task = ERROR_FACTR as i32;
    }
    
    for i in 0..n as usize {
        if nbd[i] < 0 || nbd[i] > 3 {
            *task = ERROR_NBD as i32;
            *info = -6;
            *k = (i + 1) as i32; // Convert back to 1-based indexing for error reporting
        }
        if nbd[i] == 2 {
            if l[i] > u[i] {
                *task = ERROR_FEAS as i32;
                *info = -7;
                *k = (i + 1) as i32; // Convert back to 1-based indexing for error reporting
            }
        }
    }
    
    0
}