https://github.com/cran/BDgraph
Tip revision: 3438405370c4f62bbce4b06797bbeda66307b692 authored by Reza Mohammadi on 07 April 2019, 03:22:42 UTC
version 2.57
version 2.57
Tip revision: 3438405
bdgraph.ts.R
## - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - |
# Copyright (C) 2012 - 2019 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 )
}