# 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 matplotlib
import matplotlib.pyplot as plt

import microGrid.parameters as bv
import microGrid.simulateState as simX
import microGrid.valueNext as fv


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

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

    Xt, _ = simX.residualDemand(param, param.nstep, outSampleSimul, param.maturity, X0, param.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(39,"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)

        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()