extern crate peroxide;
use peroxide::fuga::*;

#[test]
pub fn test_apply() {
    let a = ml_matrix("1 2 3;4 5 6");
    let b = c!(1, 2, 3);
    let c = a.apply(&b);

    assert_eq!(c, c!(14, 32));
}

#[cfg(feature = "O3")]
#[test]
pub fn test_dpotrf() {
    let a = py_matrix(vec![
        vec![
            2.9821852954666,
            2.680083137271475,
            2.317961273223499,
            2.134665492195094,
            3.026598331879419,
            2.283461999828011,
            3.727890627691731,
            1.948290349481972,
            2.509358003244204,
            2.354788397327008,
        ],
        vec![
            2.680083137271475,
            4.018906407367834,
            2.341820593934314,
            2.830765880629095,
            3.552151618608445,
            2.460655598487766,
            3.630225204862198,
            1.796569139726129,
            2.846699793116624,
            2.774282835105069,
        ],
        vec![
            2.317961273223499,
            2.341820593934314,
            2.038910554565652,
            1.966181894408401,
            2.536384937997141,
            1.886915794296248,
            3.136970841370761,
            1.905232489093743,
            2.01726450417587,
            2.238391500871792,
        ],
        vec![
            2.134665492195094,
            2.830765880629095,
            1.966181894408401,
            3.301422256938514,
            2.794302738298506,
            2.205849011689008,
            2.679880725689802,
            1.555759606416891,
            2.475715452431455,
            2.116394736527704,
        ],
        vec![
            3.026598331879419,
            3.552151618608445,
            2.536384937997141,
            2.794302738298506,
            4.337375745514422,
            3.604658556851256,
            4.143687579387594,
            2.385563118734139,
            3.321157832605267,
            2.949560656237908,
        ],
        vec![
            2.283461999828011,
            2.460655598487766,
            1.886915794296248,
            2.205849011689008,
            3.604658556851256,
            3.413575093545912,
            3.168622634363222,
            1.907039710137486,
            2.698561362072113,
            2.249353149769269,
        ],
        vec![
            3.727890627691731,
            3.630225204862198,
            3.136970841370761,
            2.679880725689802,
            4.143687579387594,
            3.168622634363222,
            5.142703863442649,
            2.977549597798733,
            3.210043437463944,
            3.385961680761656,
        ],
        vec![
            1.948290349481972,
            1.796569139726129,
            1.905232489093743,
            1.555759606416891,
            2.385563118734139,
            1.907039710137486,
            2.977549597798733,
            2.230209604622137,
            1.742552591231192,
            2.252632855204552,
        ],
        vec![
            2.509358003244204,
            2.846699793116624,
            2.01726450417587,
            2.475715452431455,
            3.321157832605267,
            2.698561362072113,
            3.210043437463944,
            1.742552591231192,
            3.202516498869207,
            2.228297655595811,
        ],
        vec![
            2.354788397327008,
            2.774282835105069,
            2.238391500871792,
            2.116394736527704,
            2.949560656237908,
            2.249353149769269,
            3.385961680761656,
            2.252632855204552,
            2.228297655595811,
            2.784914722669765,
        ],
    ]);

    let dpotrf_u = lapack_dpotrf(&a, UPLO::Upper).unwrap();
    let dpotrf_l = lapack_dpotrf(&a, UPLO::Lower).unwrap();

    let u = dpotrf_u.get_U().unwrap();
    let l = dpotrf_l.get_L().unwrap();

    let scipy_u = py_matrix(vec![
        vec![
            1.726900488003463,
            1.551961537963327,
            1.34226684706273,
            1.236125362766596,
            1.752618841041957,
            1.322289278213136,
            2.158717687318329,
            1.128200705840593,
            1.453099365410088,
            1.363592409455782,
        ],
        vec![
            0.,
            1.268984551541246,
            0.203843359082252,
            0.718958209858647,
            0.65576415850829,
            0.321921568031328,
            0.220631829192057,
            0.035969734306862,
            0.466156555210059,
            0.518556243550586,
        ],
        vec![
            0.,
            0.,
            0.442355231459957,
            0.362642811042644,
            0.113550376983589,
            0.104958994115387,
            0.439514016124092,
            0.867072131324973,
            -0.063762143110917,
            0.683573627458849,
        ],
        vec![
            0.,
            0.,
            0.,
            1.060662825095517,
            0.108612163548495,
            0.284560665177766,
            -0.289042488015237,
            -0.168890722850109,
            0.346460971292225,
            -0.179029330692296,
        ],
        vec![
            0.,
            0.,
            0.,
            0.,
            0.90054762326669,
            1.14736846650923,
            0.218847271328198,
            0.33819073440549,
            0.486759471799237,
            0.17930981650712,
        ],
        vec![
            0.,
            0.,
            0.,
            0.,
            0.,
            0.391212348671252,
            0.072001696614059,
            -0.069844847295321,
            -0.059587640909243,
            0.135011494209335,
        ],
        vec![
            0.,
            0.,
            0.,
            0.,
            0.,
            0.,
            0.32274897849312,
            0.109238423598569,
            -0.011482954266517,
            -0.226830086097441,
        ],
        vec![
            0.,
            0.,
            0.,
            0.,
            0.,
            0.,
            0.,
            0.211084069199511,
            0.154785258348754,
            0.218861000711332,
        ],
        vec![
            0.,
            0.,
            0.,
            0.,
            0.,
            0.,
            0.,
            0.,
            0.696448259994717,
            -0.007136319580541,
        ],
        vec![0., 0., 0., 0., 0., 0., 0., 0., 0., 0.086731850555121],
    ]);

    let scipy_l = py_matrix(vec![
        vec![1.726900488003463, 0., 0., 0., 0., 0., 0., 0., 0., 0.],
        vec![
            1.551961537963327,
            1.268984551541246,
            0.,
            0.,
            0.,
            0.,
            0.,
            0.,
            0.,
            0.,
        ],
        vec![
            1.34226684706273,
            0.203843359082252,
            0.442355231459957,
            0.,
            0.,
            0.,
            0.,
            0.,
            0.,
            0.,
        ],
        vec![
            1.236125362766596,
            0.718958209858647,
            0.362642811042644,
            1.060662825095517,
            0.,
            0.,
            0.,
            0.,
            0.,
            0.,
        ],
        vec![
            1.752618841041957,
            0.65576415850829,
            0.113550376983589,
            0.108612163548495,
            0.90054762326669,
            0.,
            0.,
            0.,
            0.,
            0.,
        ],
        vec![
            1.322289278213136,
            0.321921568031328,
            0.104958994115387,
            0.284560665177766,
            1.14736846650923,
            0.391212348671252,
            0.,
            0.,
            0.,
            0.,
        ],
        vec![
            2.158717687318329,
            0.220631829192057,
            0.439514016124092,
            -0.289042488015237,
            0.218847271328198,
            0.072001696614058,
            0.32274897849312,
            0.,
            0.,
            0.,
        ],
        vec![
            1.128200705840593,
            0.035969734306862,
            0.867072131324973,
            -0.168890722850108,
            0.33819073440549,
            -0.069844847295322,
            0.109238423598569,
            0.211084069199512,
            0.,
            0.,
        ],
        vec![
            1.453099365410088,
            0.466156555210059,
            -0.063762143110917,
            0.346460971292225,
            0.486759471799238,
            -0.059587640909244,
            -0.011482954266516,
            0.154785258348753,
            0.696448259994717,
            0.,
        ],
        vec![
            1.363592409455782,
            0.518556243550586,
            0.683573627458849,
            -0.179029330692296,
            0.17930981650712,
            0.135011494209336,
            -0.226830086097441,
            0.218861000711332,
            -0.00713631958054,
            0.086731850555118,
        ],
    ]);

    assert_eq!(u, scipy_u);
    assert_eq!(l, scipy_l);
}

#[cfg(feature = "O3")]
#[test]
pub fn test_cholesky() {
    let a = py_matrix(vec![
        vec![
            2.9821852954666,
            2.680083137271475,
            2.317961273223499,
            2.134665492195094,
            3.026598331879419,
            2.283461999828011,
            3.727890627691731,
            1.948290349481972,
            2.509358003244204,
            2.354788397327008,
        ],
        vec![
            2.680083137271475,
            4.018906407367834,
            2.341820593934314,
            2.830765880629095,
            3.552151618608445,
            2.460655598487766,
            3.630225204862198,
            1.796569139726129,
            2.846699793116624,
            2.774282835105069,
        ],
        vec![
            2.317961273223499,
            2.341820593934314,
            2.038910554565652,
            1.966181894408401,
            2.536384937997141,
            1.886915794296248,
            3.136970841370761,
            1.905232489093743,
            2.01726450417587,
            2.238391500871792,
        ],
        vec![
            2.134665492195094,
            2.830765880629095,
            1.966181894408401,
            3.301422256938514,
            2.794302738298506,
            2.205849011689008,
            2.679880725689802,
            1.555759606416891,
            2.475715452431455,
            2.116394736527704,
        ],
        vec![
            3.026598331879419,
            3.552151618608445,
            2.536384937997141,
            2.794302738298506,
            4.337375745514422,
            3.604658556851256,
            4.143687579387594,
            2.385563118734139,
            3.321157832605267,
            2.949560656237908,
        ],
        vec![
            2.283461999828011,
            2.460655598487766,
            1.886915794296248,
            2.205849011689008,
            3.604658556851256,
            3.413575093545912,
            3.168622634363222,
            1.907039710137486,
            2.698561362072113,
            2.249353149769269,
        ],
        vec![
            3.727890627691731,
            3.630225204862198,
            3.136970841370761,
            2.679880725689802,
            4.143687579387594,
            3.168622634363222,
            5.142703863442649,
            2.977549597798733,
            3.210043437463944,
            3.385961680761656,
        ],
        vec![
            1.948290349481972,
            1.796569139726129,
            1.905232489093743,
            1.555759606416891,
            2.385563118734139,
            1.907039710137486,
            2.977549597798733,
            2.230209604622137,
            1.742552591231192,
            2.252632855204552,
        ],
        vec![
            2.509358003244204,
            2.846699793116624,
            2.01726450417587,
            2.475715452431455,
            3.321157832605267,
            2.698561362072113,
            3.210043437463944,
            1.742552591231192,
            3.202516498869207,
            2.228297655595811,
        ],
        vec![
            2.354788397327008,
            2.774282835105069,
            2.238391500871792,
            2.116394736527704,
            2.949560656237908,
            2.249353149769269,
            3.385961680761656,
            2.252632855204552,
            2.228297655595811,
            2.784914722669765,
        ],
    ]);

    let u = a.cholesky(UPLO::Upper);
    let l = a.cholesky(UPLO::Lower);

    let scipy_u = py_matrix(vec![
        vec![
            1.726900488003463,
            1.551961537963327,
            1.34226684706273,
            1.236125362766596,
            1.752618841041957,
            1.322289278213136,
            2.158717687318329,
            1.128200705840593,
            1.453099365410088,
            1.363592409455782,
        ],
        vec![
            0.,
            1.268984551541246,
            0.203843359082252,
            0.718958209858647,
            0.65576415850829,
            0.321921568031328,
            0.220631829192057,
            0.035969734306862,
            0.466156555210059,
            0.518556243550586,
        ],
        vec![
            0.,
            0.,
            0.442355231459957,
            0.362642811042644,
            0.113550376983589,
            0.104958994115387,
            0.439514016124092,
            0.867072131324973,
            -0.063762143110917,
            0.683573627458849,
        ],
        vec![
            0.,
            0.,
            0.,
            1.060662825095517,
            0.108612163548495,
            0.284560665177766,
            -0.289042488015237,
            -0.168890722850109,
            0.346460971292225,
            -0.179029330692296,
        ],
        vec![
            0.,
            0.,
            0.,
            0.,
            0.90054762326669,
            1.14736846650923,
            0.218847271328198,
            0.33819073440549,
            0.486759471799237,
            0.17930981650712,
        ],
        vec![
            0.,
            0.,
            0.,
            0.,
            0.,
            0.391212348671252,
            0.072001696614059,
            -0.069844847295321,
            -0.059587640909243,
            0.135011494209335,
        ],
        vec![
            0.,
            0.,
            0.,
            0.,
            0.,
            0.,
            0.32274897849312,
            0.109238423598569,
            -0.011482954266517,
            -0.226830086097441,
        ],
        vec![
            0.,
            0.,
            0.,
            0.,
            0.,
            0.,
            0.,
            0.211084069199511,
            0.154785258348754,
            0.218861000711332,
        ],
        vec![
            0.,
            0.,
            0.,
            0.,
            0.,
            0.,
            0.,
            0.,
            0.696448259994717,
            -0.007136319580541,
        ],
        vec![0., 0., 0., 0., 0., 0., 0., 0., 0., 0.086731850555121],
    ]);

    let scipy_l = py_matrix(vec![
        vec![1.726900488003463, 0., 0., 0., 0., 0., 0., 0., 0., 0.],
        vec![
            1.551961537963327,
            1.268984551541246,
            0.,
            0.,
            0.,
            0.,
            0.,
            0.,
            0.,
            0.,
        ],
        vec![
            1.34226684706273,
            0.203843359082252,
            0.442355231459957,
            0.,
            0.,
            0.,
            0.,
            0.,
            0.,
            0.,
        ],
        vec![
            1.236125362766596,
            0.718958209858647,
            0.362642811042644,
            1.060662825095517,
            0.,
            0.,
            0.,
            0.,
            0.,
            0.,
        ],
        vec![
            1.752618841041957,
            0.65576415850829,
            0.113550376983589,
            0.108612163548495,
            0.90054762326669,
            0.,
            0.,
            0.,
            0.,
            0.,
        ],
        vec![
            1.322289278213136,
            0.321921568031328,
            0.104958994115387,
            0.284560665177766,
            1.14736846650923,
            0.391212348671252,
            0.,
            0.,
            0.,
            0.,
        ],
        vec![
            2.158717687318329,
            0.220631829192057,
            0.439514016124092,
            -0.289042488015237,
            0.218847271328198,
            0.072001696614058,
            0.32274897849312,
            0.,
            0.,
            0.,
        ],
        vec![
            1.128200705840593,
            0.035969734306862,
            0.867072131324973,
            -0.168890722850108,
            0.33819073440549,
            -0.069844847295322,
            0.109238423598569,
            0.211084069199512,
            0.,
            0.,
        ],
        vec![
            1.453099365410088,
            0.466156555210059,
            -0.063762143110917,
            0.346460971292225,
            0.486759471799238,
            -0.059587640909244,
            -0.011482954266516,
            0.154785258348753,
            0.696448259994717,
            0.,
        ],
        vec![
            1.363592409455782,
            0.518556243550586,
            0.683573627458849,
            -0.179029330692296,
            0.17930981650712,
            0.135011494209336,
            -0.226830086097441,
            0.218861000711332,
            -0.00713631958054,
            0.086731850555118,
        ],
    ]);

    assert_eq!(u, scipy_u);
    assert_eq!(l, scipy_l);
}

#[cfg(feature = "O3")]
#[test]
#[should_panic]
pub fn test_cholesky_panic_1() {
    // Not symmetric
    let a = ml_matrix("1 2;3 10");
    a.cholesky(UPLO::Upper).print();
}

#[cfg(feature = "O3")]
#[test]
#[should_panic]
pub fn test_cholesky_panic_2() {
    // Not Positive Definite
    let a = ml_matrix("1 2;2 1");
    a.cholesky(UPLO::Upper).print();
}