# Copyright (C) 2018 The Regents of the University of California, Michael Ludkovski and Aditya Maheshwari
# All Rights Reserved
# This code is published under the GNU Lesser General Public License (GNU LGPL)

from __future__ import division
import numpy as np
import math
import sys
import matplotlib
import matplotlib.pyplot as plt

import microgridDEA.parameters as bv
import microgridDEA.simulateState_new as simX
import microgridDEA.valueNext_new as fv
if (sys.version_info > (3, 0)):
        import pickle
else:
        import cPickle as pickle

import StOptGrids
import StOptReg
import StOptGeners


def currentCost(param,St, dt):
    return (param.c1*(dt**param.a))*param.dt,  param.c2*np.where(St>0,St,0)*param.dt + param.c3*np.where(St<0,-St,0)*param.dt


def forwardSimulations():

    outSampleSimul = 10000
    maturityOutsample = 10

    param = bv.basicVariables(maturityOutsample)

    condExpFilePath = param.filetoDump

    with open(param.regParamDump, 'rb') as input:
        regParameters = pickle.load(input)

    X0 = np.zeros(outSampleSimul)
    X0[:] = 0  #initial state for residual demand

    Xt = simX.residualDemand(param, param.nstep, outSampleSimul, param.maturity, X0, regParameters[0].rmctype, "outsample")

    r,c = Xt.shape

    It = np.zeros((r,c)) # inventory evolution
    Dt = np.zeros((r,c), dtype=np.int) # regime 
    Dt[:,0] = 0  #initial regime
    
    St = np.zeros((r,c)) # stores Xt - dt - Bt
    Bt = np.zeros((r,c)) # battery output
    dt = np.zeros((r,c)) # diesel output
    costTrue = np.zeros((r,c)) # cost


    It[:,0] = param.I0 #initial state for battery inventory

    archiveToRead =  StOptGeners.BinaryFileArchive(condExpFilePath,"r")


    for iSteps in range(param.nstep):

        contValues = archiveToRead.readGridAndRegressedValue(param.nstep-iSteps,"toStore")

        for sim in range(outSampleSimul):
            _, dt[sim, iSteps], St[sim, iSteps], It[sim, iSteps+1] ,Bt[sim, iSteps] = fv.findOptimalControl(Xt[sim,iSteps], It[sim,iSteps], Dt[sim,iSteps], contValues[0], param, regParameters[iSteps])

        Dt[:, iSteps+1] = np.where(dt[:, iSteps] > 0, 1, 0)

        costTrue[:,iSteps] = np.where((Dt[:,iSteps]==0) & (Dt[:,iSteps+1]==1),param.K,0)


    DieselCost, ArtificalCost = currentCost(param,St, dt) 
    costTrue = costTrue + DieselCost
    costTrue[:,-1] = fv.penalty(param,It[:,-1])



    cumCostTrue = np.cumsum(costTrue,1)

    print("Cost: Total", np.mean(np.sum(costTrue+ArtificalCost,1)), "true", np.mean(np.sum(costTrue,1)), "artifical", np.mean(np.sum(ArtificalCost,1)))
    print("std: Total", np.std(np.sum(costTrue+ArtificalCost,1))/np.sqrt(outSampleSimul), "true", np.std(np.sum(costTrue,1))/np.sqrt(outSampleSimul), "artifical", np.std(np.sum(ArtificalCost,1))/np.sqrt(outSampleSimul))

    # plots to show curtailed energy and proportion of blackouts

    # St1 = np.where(St<-10**-9,1,0)
    # plt.plot(np.sum(St1,0)[:-2]/St1.shape[0])
    # plt.show()


    # St1 = np.where(St>10**-9,1,0)
    # plt.plot(np.sum(St1,0)[:-2]/St1.shape[0])
    # plt.show()
    
if __name__ == "__main__":
    forwardSimulations()