## - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - |
#     Copyright (C) 2012 - 2018  Reza Mohammadi                                                    |
#                                                                                                  |
#     This file is part of BDgraph package.                                                        |
#                                                                                                  |
#     BDgraph is free software: you can redistribute it and/or modify it under                     |
#     the terms of the GNU General Public License as published by the Free                         |
#     Software Foundation; see <https://cran.r-project.org/web/licenses/GPL-3>.                    |
#                                                                                                  |
#     Maintainer: Reza Mohammadi <a.mohammadi@uva.nl>                                              |
## - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - |
#     BDMCMC algorithm for time series                                                             |
## - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - |
#     Nlength is the length of the time series                                                     |
#     data is the aggregate periodogram Pk, which is arranged as                                   |
#     a large p x ( Nlength * p ) matrix [ P1, P2, ... , PNlength ]                                |
## - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - |

bdgraph.ts = function( data, Nlength = NULL, n, iter = 1000, burnin = iter / 2, 
                       g.prior = 0.5, df.prior = rep( 3, Nlength ), g.start = "empty", 
                       save = FALSE, print = 500, cores = NULL )
{
    num_machine_cores = BDgraph::detect_cores()
    if( is.null( cores ) ) cores = num_machine_cores - 1
    if( cores == "all" )   cores = num_machine_cores
    
    .C( "omp_set_num_cores", as.integer( cores ), PACKAGE = "BDgraph" )
    
    burnin = floor( burnin )
    
    if( !is.matrix( data ) & !is.data.frame( data ) ) stop( "Data should be a matrix or dataframe" )
    if( is.data.frame( data ) ) data <- data.matrix( data )
    if( iter <= burnin )   stop( "Number of iteration must be more than number of burn-in" )
    
    if( any( is.na( data ) ) ) stop( "This method does not deal with missing value" )
    
    dimd <- dim( data )
    p    <- dimd[ 1 ]
    if( is.null( Nlength ) ) Nlength <- dimd[ 2 ] / p
    
    if( !is.matrix( g.prior ) )
    {
        if( ( g.prior <= 0 ) | ( g.prior >= 1 ) ) stop( " 'g.prior' must be between 0 and 1" )
        g.prior = matrix( g.prior, p, p )
    }else{
        if( ( nrow( g.prior ) != p ) | ( ncol( g.prior ) != p ) ) stop( " 'g.prior' and 'data' have non-conforming size" )
        if( any( g.prior < 0 ) || any( g.prior > 1 ) ) stop( " Element of 'g.prior', as a matrix, must be between 0 and 1" )
    }
    g_prior = g.prior
    
    P = data 
    W = matrix( 0, p, p * Nlength )  # I( p * p )
    
    for( t in 1 : Nlength )
    diag( W[ , ( t * p - p + 1 ):( t * p ) ] ) = 1
    
    K      = W
    sigma  = K
    b      = df.prior
    b_star = b + n
    Ws     = W + P
    
    Ls = matrix( 0, 2 * p, 2 * p * Nlength ) # Ls = t(chol([Re(inv_Pk), -Im[inv_Pk]; Re(inv_Pk), Im[inv_Pk]]))
    
    for( k in 1:Nlength )
    {
        inv_Ws = solve( Ws[, ( k * p - p + 1 ):( k * p )] )
        row1   = cbind( Re( inv_Ws ), -Im( inv_Ws ) )
        row2   = cbind( Im( inv_Ws ), Re( inv_Ws ) )
        sig    = rbind( row1, row2 ) / 2
        Ls[ , ( k * 2 * p - 2 * p + 1 ):( k * 2 * p ) ] = t( chol( sig ) )
    }
    
    if( class( g.start ) == "bdgraph" ) 
    {
        G <- g.start $ last_graph
        K <- g.start $ last_K
        if( dim( K )[ 2 ] != p * Nlength ) stop( " K should be a p x ( p x Nlength ) matrix" )
    } 
    
    if( class( g.start ) == "sim" ) 
    {
        G <- as.matrix( g.start $ G )
        K <- as.matrix( g.start $ K )
        if( dim( K )[ 2 ] != p * Nlength ) stop( " K should be a p x (p x Nlength) matrix" )
    } 
    
    if( class( g.start ) == "character" && g.start == "empty" )
    {
        G = matrix( 0, p, p )
        
        for( t in 1:Nlength )
        {
            result = .C( "rgcwish_c", as.integer(G), as.double(Ls), K = as.complex(K), as.integer(b_star), as.integer(p), PACKAGE = "BDgraph" )
            K[, ( t * p - p + 1 ):( t * p )] = matrix ( result $ K, p, p )
            sigma[ , ( t * p - p + 1 ):( t * p ) ] = solve( K[ , ( t * p - p + 1 ):( t * p ) ] )
        }
    }
    
    if( class( g.start ) == "character" && g.start == "full" )
    {
        G         = matrix( 1, p, p )
        diag( G ) = 0
        
        for( t in 1:Nlength )
        {
            result = .C( "rcwish_c", as.double(Ls), K = as.complex(K), as.integer(b_star), as.integer(p), PACKAGE = "BDgraph" )
            K[ , ( t * p - p + 1 ):( t * p ) ] = matrix ( result $ K, p, p ) 
            sigma[ , ( t * p - p + 1 ):( t * p ) ] = solve( K[ , ( t * p - p + 1 ):( t * p ) ] )
        }	
    }
    
    if( is.matrix( g.start ) )
    {
        G         = g.start
        diag( G ) = 0
        
        for( t in 1:Nlength )
        {
            result = .C( "rgcwish_c", as.integer(G), as.double(Ls), K = as.complex(K), as.integer(b_star), as.integer(p), PACKAGE = "BDgraph" )
            K[ , ( t * p - p + 1 ):( t * p ) ] = matrix ( result $ K, p, p )
            sigma[ , ( t * p - p + 1 ):( t * p ) ] = solve( K[ , ( t * p - p + 1 ):( t * p ) ] )
        }	
    }
    
    if( save == TRUE )
    {
        qp1           = ( p * ( p - 1 ) / 2 ) + 1
        string_g      = paste( rep( 0, qp1 ), collapse = '' )
        sample_graphs = rep ( string_g, iter - burnin )  # vector of numbers like "10100" 
        graph_weights = rep ( 0, iter - burnin )         # waiting time for every state
        all_graphs    = rep ( 0, iter - burnin )         # vector of numbers like "10100"
        all_weights   = rep ( 1, iter - burnin )         # waiting time for every state		
        size_sample_g = 0
        exit          = 0
    }else{
        p_links = 0 * G
    }
    
    if( ( save == TRUE ) && ( p > 50 & iter > 20000 ) )
    {
        cat( "  WARNING: Memory needs to run this function is around " )
        print( ( iter - burnin ) * utils::object.size( string_g ), units = "auto" ) 
    } 
    
    r_K        = Re( K )
    i_K        = Im( K )
    r_K_hat    = 0 * K
    i_K_hat    = 0 * K
    r_sigma    = Re( sigma )
    i_sigma    = Im( sigma )
    last_graph = G
    
    mes <- paste( c( iter, " iteration is started.                    " ), collapse = "" )
    #cat( mes, "\r" )
    
    if ( save == TRUE )
    {
        result = .C( "bdmcmc_map_for_multi_dim", as.integer(iter), as.integer(burnin), G = as.integer(G), as.double(Ls), r_K = as.double(r_K), 
        i_K = as.double(i_K), as.integer(p), as.integer(Nlength), r_sigma = as.double(r_sigma), i_sigma = as.double(i_sigma), 
        all_graphs = as.integer(all_graphs), all_weights = as.double(all_weights), r_K_hat = as.double(r_K_hat), i_K_hat = as.double(i_K_hat), 
        sample_graphs = as.character(sample_graphs), graph_weights = as.double(graph_weights), size_sample_g = as.integer(size_sample_g), 
        exit = as.integer(exit), as.integer(b), as.integer(b_star), r_Ds = as.double(Re(Ws)), i_Ds = as.double(Im(Ws)), PACKAGE = "BDgraph" )
    }else{
        result = .C( "bdmcmc_for_multi_dim", as.integer(iter), as.integer(burnin), G = as.integer(G), as.double(g_prior), as.double(Ls), 
        r_K = as.double(r_K), i_K = as.double(i_K), as.integer(p), as.integer(Nlength), r_sigma = as.double(r_sigma), i_sigma = as.double(i_sigma), 
        r_K_hat = as.double(r_K_hat), i_K_hat = as.double(i_K_hat), p_links = as.double(p_links), as.integer(b), as.integer(b_star), 
        r_Ds = as.double(Re(Ws)), i_Ds = as.double(Im(Ws)), print = as.integer(print), PACKAGE = "BDgraph" )  
    }
    
    r_K_hat = result $ r_K_hat
    i_K_hat = result $ i_K_hat
    K_hat   = matrix( complex( 1, r_K_hat, i_K_hat ), p, p * Nlength )
    
    last_graph = matrix( result $ G, p, p )
    
    last_r_K = result $ r_K
    last_i_K = result $ i_K
    last_K   = matrix( complex( 1, last_r_K, last_i_K ), p, p * Nlength )
    
    if( save == TRUE )
    {
        status = result $ exit
        
        if( status )
        {
            mes <- paste( c( "Exit at iteration ", status ), collapse = "" )
            cat( mes, "\r" )
        }
        
        if ( status == 0 | status >= burnin - 1 )
        {
            size_sample_g = result $ size_sample_g
            sample_graphs = result $ sample_graphs[ 1 : size_sample_g ]
            graph_weights = result $ graph_weights[ 1 : size_sample_g ]
            all_graphs    = result $ all_graphs + 1
            all_weights   = result $ all_weights
            
            output = list( sample_graphs = sample_graphs, graph_weights = graph_weights, K_hat = K_hat, 
            all_graphs = all_graphs, all_weights = all_weights, last_graph = last_graph, last_K = last_K, status = status)
        }else{
            p_links   = matrix( result $ G, p, p ) 
            p_links[ lower.tri( p_links ) ] = 0
            output = list( p_links = p_links, K_hat = K_hat, last_graph = last_graph, last_K = last_K )
        }
    }else{
        p_links = matrix( result $ p_links, p, p ) 
        p_links[ lower.tri( p_links ) ] = 0
        
        output = list( p_links = p_links, K_hat = K_hat, last_graph = last_graph, last_K = last_K )
    }
    
    class( output ) = "bdgraph"
    return( output )   
}