/** \file
 *
 *  Contains the HFRadarCompactModelForecaster class implementation.
 *
 *  Copyright (c) 2013, 2014 MBARI
 *  MBARI Proprietary Information.  All Rights Reserved
 */

#include "HFRadarCompactModelForecaster.h"
#include "HFRadarCompactModelForecasterIF.h"

#include "data/BlobWriter.h"
#include "data/ConfigReader.h"
#include "data/Location.h"
#include "data/UniversalDataReader.h"
#include "utils/Timestamp.h"

#include <cerrno>
#include <cstddef>
#include <sys/stat.h>
#include <sys/types.h>
#include <dirent.h>

HFRadarCompactModelForecaster::HFRadarCompactModelForecaster( const Module* module )
    : SyncEstimationComponent( HFRadarCompactModelForecasterIF::NAME, module ), // TODO: consider making this an asynchronous component
      verbosity_( 0 ),
      // load the neural net parameters
      xlags_( "Resources/HFRadarModel/net.inputDataHistory.mtx", logger_ ),
      ulags_( "Resources/HFRadarModel/net.forcingsHistory.mtx", logger_ ),
      netInpMean_( "Resources/HFRadarModel/net.inpMean.mtx", logger_ ),
      netTpca_( "Resources/HFRadarModel/net.Tpca.mtx", logger_ ),
      netW1_( "Resources/HFRadarModel/net.w1.mtx", logger_ ),
      netB1_( "Resources/HFRadarModel/net.b1.mtx", logger_ ),
      // initialize other member variables
      historyHours_( ( int )( xlags_( 0, -1 ) - xlags_( 0, 0 ) + 1 ) ),
      projectionHours_( HFRadarCompactModelForecasterIF::PROJECTION_HOURS ),
      expansionCoefficients_( logger_, HFRadarCompactModelForecasterIF::MAX_EOFS, historyHours_ + projectionHours_ ),
      forecastExpansionCoefficients_( logger_, HFRadarCompactModelForecasterIF::MAX_EOFS, projectionHours_ ),
      forecastStartTime_( Timestamp::NOT_SET_TIME )
{
    ignoreECsMoreRecentThanCfgReader_ = newConfigReader( HFRadarCompactModelForecasterIF::IGNORE_ECS_MORE_RECENT_THAN_CFG );
    platformCommunicationsReader_ = newUniversalReader( UniversalURI::PLATFORM_COMMUNICATIONS );
    forecastExpansionCoefficientsWriter_ = newBlobWriter( HFRadarCompactModelForecasterIF::FORECAST_EXPANSION_COEFFICIENTS );
    forecastTimesWriter_ = newBlobWriter( HFRadarCompactModelForecasterIF::FORECAST_TIMES );
}

HFRadarCompactModelForecaster::~HFRadarCompactModelForecaster()
{}

void HFRadarCompactModelForecaster::initialize( void )
{
    if( loadHistory() ) // always try to load the history, compute a forecast, and publish it on the first cycle
    {
        forecast(); // propogate the historical expansion coefficients through the neural net to produce the forecast
        publishForecast(); // write the expansion coefficients and associated times out to the slate
    }
    run();
}

void HFRadarCompactModelForecaster::run( void )
{
    if( platformCommunicationsReader_->wasTouchedSinceLastRun( this ) ) // only check for new history if we have actually had comms
    {
        if( loadHistory() ) // Only continue if you can load new history
        {
            forecast(); // propogate the historical expansion coefficients through the neural net to produce the forecast
            publishForecast(); // write the expansion coefficients and associated times out to the slate
        }
    }
}

void HFRadarCompactModelForecaster::forecast( void )
{
    //
    Mtx x( logger_, 1, HFRadarCompactModelForecasterIF::MAX_EOFS * xlags_.getN() + ulags_.getN(), 0.0f );
    for( int i = 0; i < projectionHours_; ++i )
    {
        for( int j = 0; j < xlags_.getN(); ++j )
        {
            int xInd = historyHours_ + ( int )xlags_( 0, xlags_.getN() - j - 1 ) + i - 1;
            for( int k = 0; k < HFRadarCompactModelForecasterIF::MAX_EOFS; ++k )
            {
                x[0][j * HFRadarCompactModelForecasterIF::MAX_EOFS + k] = expansionCoefficients_( k, xInd );
            }
        }
        Mtx ec = glm2fwd( x );
        for( int j = 0; j < HFRadarCompactModelForecasterIF::MAX_EOFS; ++j )
        {
            expansionCoefficients_[j][historyHours_ + i] = ec( 0, j );
        }
    }

    // separate out the forecast ECs from the history ECs
    forecastExpansionCoefficients_ = forecastExpansionCoefficients_.subset( 0, -1, historyHours_, -1 );

    // generate the timestamps for the different sets of ECs
    double forecastStartEpoch = forecastStartTime_.asDouble();
    for( int j = 0; j < projectionHours_ + historyHours_; j++ )
    {
        forecastEpochSeconds_[j] = forecastStartEpoch + ( j - historyHours_ + 1 ) * 3600.0;
    }
}

Mtx HFRadarCompactModelForecaster::glm2fwd( const Mtx& buf )
{
    Mtx mtx( buf );
    mtx -= netInpMean_;
    mtx = mtx * netTpca_; // (1x463)*(463x428)= 1x428
    mtx = mtx * netW1_; // (1x428)*(428xMAX_EOFS)=1xMAX_EOFS
    mtx += netB1_;
    return mtx;
}

void HFRadarCompactModelForecaster::publishForecast( void )
{
    // TODO: Consider writing this out as the transpose (hours in rows, modes in columns) for consistency with the shore-side netCDF4 outputs
    Timestamp calculationCompleteTime = Timestamp::Now(); // TODO: move this to right after the glm2fwd calls
    // TODO: evaluate whether to use expansionCoefficients_ or forecastExpansionCoefficients_
    logger_.syslog( "forecast time " + forecastStartTime_.toSmallerString() + " published", Syslog::IMPORTANT );
    if( verbosity_ > 1 )
    {
        char msg[256];
        sprintf( msg, "ec[0][0] is %g, ec[0][1] is %g, ec[1][0] is %g", expansionCoefficients_.get( 0, 0 ), expansionCoefficients_.get( 0, 1 ), expansionCoefficients_.get( 1, 0 ) );
        logger_.syslog( msg, Syslog::DEBUG );
    }
    forecastExpansionCoefficientsWriter_->write2DClass( Units::NONE, expansionCoefficients_, calculationCompleteTime );
    forecastTimesWriter_->write1DArray( Units::EPOCH_SECOND, forecastEpochSeconds_, calculationCompleteTime );
    if( verbosity_ > 1 ) logger_.syslog( "expansionCoefficients_[0][" + Str( HFRadarCompactModelForecasterIF::HISTORY_HOURS - 1 ) + "] = " + Str( expansionCoefficients_[0][HFRadarCompactModelForecasterIF::HISTORY_HOURS - 1] ) + " for " + Timestamp( forecastEpochSeconds_[HFRadarCompactModelForecasterIF::HISTORY_HOURS - 1] ).toString(), Syslog::IMPORTANT );
}

bool HFRadarCompactModelForecaster::loadHistory( void )
{
    // Test existence of the directory that contains the historical ECs.
    Str hfRadarPath( "Data/HFRadarModel/" );
    int error = mkdir( hfRadarPath.cStr(), ACCESSPERMS );
    if( error && EEXIST != errno )
    {
        logger_.syslog( "Unable to access " + hfRadarPath, Syslog::ERROR );
        return false;
    }

    // Test access to the directory that contains the historical ECs.
    DIR* hfRadarDir( opendir( hfRadarPath.cStr() ) );
    if( NULL == hfRadarDir )
    {
        logger_.syslog( "Unable to open " + hfRadarPath, Syslog::ERROR );
        return false;
    }
    union // To use thread-safe readdir_r make sure dirent struct is big enough
    {
        struct dirent ent_;
        char b_[offsetof( struct dirent, d_name ) + NAME_MAX + 1];
    } dir;
    struct dirent* entry;

    Timestamp maxExpansionCoefficientTimestamp( Timestamp::NOT_SET_TIME );
    Timestamp thisExpansionCoefficientTimestamp( Timestamp::NOT_SET_TIME );
    float ignoreECsMoreRecentThanHoursCfg( 0.0 );
    ignoreECsMoreRecentThanCfgReader_->read( Units::HOUR, ignoreECsMoreRecentThanHoursCfg );

    // find the most recent file of historical ECs
    while( 0 == readdir_r( hfRadarDir, &dir.ent_, &entry ) && entry != NULL )
    {
        int year( 0 ), month( 0 ), day( 0 ), hour( 0 );
        if( 4 != sscanf( entry->d_name, "%04d%02d%02d%02d", &year, &month, &day, &hour ) )
        {
            continue;
        }
        thisExpansionCoefficientTimestamp = Timestamp( year, month, day, hour, 0, 0 );
        if( thisExpansionCoefficientTimestamp > Timestamp::Now() - Timespan::Hours( ignoreECsMoreRecentThanHoursCfg ) )
        {
            if( verbosity_ > 0 ) logger_.syslog( "Found ECs from " + thisExpansionCoefficientTimestamp.toSmallString() + ", more recent than " + Str( ignoreECsMoreRecentThanHoursCfg ) + " hours -- perhaps running a simulation in the past?", Syslog::INFO );
            // TODO: Can't just use a break here because readdir_r doesn't appear to read a *sorted* list of the directory contents
            // break; // do not try to read _very_ recent ECs or those in the future (useful for simulations)
        }
        else if( thisExpansionCoefficientTimestamp > maxExpansionCoefficientTimestamp )
        {
            if( verbosity_ > 0 ) logger_.syslog( "found more recent ECs from " + thisExpansionCoefficientTimestamp.toSmallString(), Syslog::INFO );
            maxExpansionCoefficientTimestamp = thisExpansionCoefficientTimestamp;
        }
    }
    closedir( hfRadarDir ); // Avoid a big memory leak

    if( maxExpansionCoefficientTimestamp > forecastStartTime_ ) // actually read the ECs into expansionCoefficient_
    {
        logger_.syslog( "Found new ECs!", Syslog::DEBUG );
        forecastStartTime_ = maxExpansionCoefficientTimestamp;
        if( forecastStartTime_ < Timestamp::Now() - Timespan::Hours( 6.0 ) )
        {
            logger_.syslog( "Most recent ECs are more than six hours old.", Syslog::ERROR );
        }
        for( int j = 0; j < historyHours_; ++j )
        {
            int historyHour( j - historyHours_ + 1 );
            Timespan dt( historyHour * 3600.0 );
            Timestamp inTime( forecastStartTime_ + dt );
            inTime.round( 3600 ); // just to be safe? shouldn't this already be incrementing by 3600?
            char inName[] = "Data/HFRadarModel/yyyymmddhh";
            struct tm tm( inTime.asStructTm() );
            strftime( inName + sizeof( inName ) - 11, 11, "%Y%m%d%H", &tm );
            Mtx temp( ( char* )inName, logger_ ); // read in the historical ECs from the file
            for( int i = 0; i < HFRadarCompactModelForecasterIF::MAX_EOFS; ++i )
            {
                expansionCoefficients_[i][j] = temp.getM() > i && temp.getN() > 0 ? temp[i][0] : 0.0f;
            }
        }
        return true;
    }
    else
    {
        logger_.syslog( "Did not find new historical expansion coefficients.", Syslog::DEBUG );
        return false; // could not load any NEW historical data
    }
}