1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205

extern crate libc;

use std::f64;
use libc::{c_double, c_int};

mod ffi;

pub trait GLMNet {
    fn is_dim_ok(&self, dim: (usize, usize)) -> bool;
    fn fit(&self, y: &[f64], x: &[f64]) -> GLMNetPath;
}

#[derive(Debug)]
pub enum Family {
    Normal,
    Binomial,
    Multinomial,
    Poisson,
}

impl Default for Family {
    fn default() -> Family {
        Family::Normal
    }
}

#[derive(Debug)]
pub struct GLMNetPath {
    pub family: Family,
    pub a0: Vec<f64>,
    pub betas: Vec<f64>,
    pub rsq: Vec<f64>,
    pub lambda: Vec<f64>,
    pub nlp: i32,
}

impl GLMNetPath {
    fn with_capacity(nftr: usize, nlambda: usize) -> Self {
        GLMNetPath {
            family: Default::default(),
            a0: vec![0f64; nlambda],
            betas: vec![0f64; nftr * nlambda],
            rsq: vec![0f64; nlambda],
            lambda: vec![0f64; nlambda],
            nlp: 0,
        }
    }
}

#[derive(Debug)]
pub struct GLMNetLR {
    alpha: c_double,
    weights: Vec<c_double>,
    cstr: Vec<c_double>,
    penal: Vec<c_double>,
    dfmax: c_int,
    pmax: c_int,
    nlambda: c_int,
    flmin: c_double,
    ulambdas: Vec<c_double>,
    thres: c_double,
    isd: c_int,
    intr: c_int,
    maxit: c_int,
    dim: (usize, usize),
}

impl Default for GLMNetLR {
    fn default() -> GLMNetLR {
        GLMNetLR {
            alpha: 0.5,
            weights: Vec::new(),
            cstr: Vec::new(),
            penal: Vec::new(),
            dfmax: 0,
            pmax: 0,
            nlambda: 100,
            flmin: 1e-4,
            ulambdas: Vec::new(),
            thres: 1e-7,
            isd: 1,
            intr: 0,
            maxit: 1_000_000,
            dim: (0, 0),
        }
    }
}

impl GLMNetLR {
    pub fn new(nobs: usize, nftr: usize) -> Self {
        let mut p: GLMNetLR = Default::default();
        p.pmax(nftr)
            .dim((nobs, nftr))
            .dfmax(nftr)
            .weights(vec![1f64; nobs])
            .penal(vec![1f64; nftr])
            .range((0..(nftr*2)).map(|i| {
                if i % 2 == 0 {
                    f64::NEG_INFINITY
                } else {
                    f64::INFINITY
                }
            }).collect());
        p
    }

    pub fn alpha(&mut self, alpha: f64) -> &mut Self {
        self.alpha = alpha; self
    }

    pub fn dfmax(&mut self, dfmax: usize) -> &mut Self {
        self.dfmax = dfmax as c_int; self
    }

    pub fn pmax(&mut self, pmax: usize) -> &mut Self {
        self.pmax = pmax as c_int ; self
    }

    pub fn weights(&mut self, w: Vec<f64>) -> &mut Self {
        self.weights = w; self
    }

    pub fn range(&mut self, r: Vec<f64>) -> &mut Self {
        self.cstr = r; self
    }

    pub fn dim(&mut self, dim: (usize, usize)) -> &mut Self {
        self.dim = dim; self
    }

    pub fn penal(&mut self, p: Vec<f64>) -> &mut Self {
        self.penal = p; self
    }
}

impl GLMNet for GLMNetLR {
    fn is_dim_ok(&self, dim: (usize, usize)) -> bool { dim == self.dim }

    fn fit(&self, y: &[f64], x: &[f64]) -> GLMNetPath {
        let mut path = GLMNetPath::with_capacity(self.dim.1, self.nlambda as usize);
        let (nr, nc) = (self.dim.0 as c_int, self.dim.1 as c_int);

        let mut temp = vec![0f64; self.pmax as usize * self.nlambda as usize];
        let mut ia = vec![0; self.pmax as usize];
        let mut nin = vec![0; self.nlambda as usize];
        let mut jerr = 0;
        let mut lmu = 0;

        unsafe {
            ffi::elnet_(
                &1,
                &self.alpha,
                &nr,
                &nc,
                x.as_ptr(),
                y.as_ptr(),
                self.weights.as_ptr(),
                &0,
                self.penal.as_ptr(),
                self.cstr.as_ptr(),
                &nc,
                &self.pmax,
                &self.nlambda,
                &1e-6,
                self.ulambdas.as_ptr(),
                &self.thres,
                &self.isd,
                &self.intr,
                &self.maxit,
                &mut lmu,
                path.a0.as_mut_ptr(),
                temp.as_mut_ptr(),
                ia.as_mut_ptr(),
                nin.as_mut_ptr(),
                path.rsq.as_mut_ptr(),
                path.lambda.as_mut_ptr(),
                &mut path.nlp,
                &mut jerr,
            );

            ffi::solns_(
               &nc, &self.pmax, &lmu, temp.as_ptr(), ia.as_ptr(), nin.as_ptr(),
               path.betas.as_mut_ptr(),
            );
        }

        if path.lambda.len() > 2 {
            path.lambda[0] = f64::exp(2f64 * f64::ln(path.lambda[1]) - f64::ln(path.lambda[2]));
        }

        path.a0.truncate(lmu as usize);
        path.betas.truncate(lmu as usize * self.dim.1);
        path.rsq.truncate(lmu as usize);
        path.lambda.truncate(lmu as usize);
        path
    }
}

#[cfg(test)]
mod tests {
    #[test]
    fn it_works() {
        assert_eq!(2 + 2, 4);
    }
}