flow_calc_hourly.sagews

# %auto %default_mode python3
%auto preparser(False)

import sklearn
import sklearn.svm
sklearn.__version__

'0.18.1'

sklearn.svm.SVC
sklearn.svm.LinearSVC

<class 'sklearn.svm.classes.SVC'> <class 'sklearn.svm.classes.LinearSVC'>

import warnings
warnings.filterwarnings("ignore")
import sys
import os
#sys.path.insert(0, os.path.expanduser('~/.local/lib/python3.4/site-packages/'))
print(sys.executable)
print(sys.version_info)
import datetime
import os
import sklearn
from sklearn import svm, calibration
from sklearn import preprocessing, metrics
import random
import matplotlib.pyplot as plt
import numpy as np
from numpy import cos, sin, pi
# from sage.all import Point`

def get_shift_factors(constraint_name,contingency_name=''):
    sf_dict=dict()
    c_lz=dict()
    f_shift_factor = open(os.path.expanduser("~/shift_factors"),mode='r')
    for line in f_shift_factor:
        line = line.replace("\n","")
        lsplit=  line.split(",")
        if constraint_name in lsplit[3] and contingency_name in lsplit[2]:
            if lsplit[1] not in c_lz:
                c_lz[lsplit[1]] = 0
                sf_dict[lsplit[1]]=0
            c_lz[lsplit[1]]+=1
            sf_dict[lsplit[1]]+=float(lsplit[0])
    for i in sf_dict:
        sf_dict[i]=sf_dict[i]/c_lz[i]
    f_shift_factor.close()
    return sf_dict
    
def get_wind_forecasted():
    f_wind = open(os.path.expanduser("~/wind.csv"),mode='r')
    w_dict = dict()
    hour = datetime.timedelta(hours=int(1))
    w_dict['WEST']=dict()
    w_dict['SOUTH']=dict()
    for line in f_wind:
        line = line.replace("\n","")
        lsplit = line.split(",")
        f_date = datetime.datetime.strptime(lsplit[0].split(" ")[0],"%Y-%m-%d")+int(lsplit[1])*hour
        if 'STWPF_WEST' in lsplit[2]:
            w_dict['WEST'][f_date]=float(lsplit[3])
        if "STWPF_SOUTH" in lsplit[2]:
            w_dict['SOUTH'][f_date]=float(lsplit[3])
    f_wind.close()
    return w_dict

def get_load_forecasted():
    f_load = open(os.path.expanduser("~/load.csv"),mode='r')
    hour = datetime.timedelta(hours=int(1))
    l_dict=dict()
    for line in f_load:
        lsplit = line.split(",")
        f_date = datetime.datetime.strptime(lsplit[0].split(" ")[0],"%Y-%m-%d")+int(lsplit[1].split(":")[0])*hour
        zone = lsplit[2]
        if zone not in l_dict:
            l_dict[zone] = dict()
        l_dict[zone][f_date]=float(lsplit[3])
    f_load.close()
    return l_dict

def flow_hourly(constraint_name,contingency_name=''):
    hour = datetime.timedelta(hours=int(1))
    l_dict = dict()
    f_dict = dict()
    c_lz = dict()
    flow_dict = dict()
    #first do the shift factor dictionary
    sf_dict=get_shift_factors(constraint_name,contingency_name)
    #now create the wind_dictionary
    w_dict=get_wind_forecasted()
    #now create the load_dictionary
    l_dict=get_load_forecasted()
    #now combine using shift_factors to get the hourly flows
    for zone in l_dict:
        for sf_z in sf_dict:
            if zone.lower() in sf_z.lower() and "LZ" in sf_z:
                for d in l_dict[zone]:
                    if d in w_dict['WEST'] and d in w_dict['SOUTH']:
                        if d not in flow_dict:
                            flow_dict[d] = 0
                        flow_dict[d]=flow_dict[d]-l_dict[zone][d]*sf_dict[sf_z]
    for d in w_dict['WEST']:
        if d in w_dict['SOUTH'] and d in flow_dict:
            flow_dict[d]=flow_dict[d]+w_dict['WEST'][d]*sf_dict['LZ_WEST']+w_dict['SOUTH'][d]*sf_dict['LZ_SOUTH']
    return flow_dict


def binding_hours(constraint_name,contingency_name=''):
    b_hours=  dict()
    f=open(os.path.expanduser("~/shadows"))
    for line in f:
        line = line.replace("\n","")
        lsplit = line.split(",")
        f_date = datetime.datetime.strptime(lsplit[0],"%Y-%m-%d %H:%M:%S")
        try:
            sp = float(lsplit[1])
        except:
            sp=0
        if constraint_name in lsplit[4] and contingency_name in lsplit[3] and sp>=0:
            b_hours[f_date]=1
    print(f.readline())
    return b_hours

def con_frequency():
    n_hours=  dict()
    f=open(os.path.expanduser("~/shadows"))
    for line in f:
        line = line.replace("\n","")
        lsplit = line.split(",")
        try:
            sp = float(lsplit[1])
        except:
            sp=0
        if lsplit[4] not in n_hours:
            n_hours[lsplit[4]] = 0
        if sp>0:
            n_hours[lsplit[4]] += 1
    f.close()
    return n_hours


def ML(binding_hours,flow_dict):
    clf = svm.LinearSVC(penalty='l2',class_weight='balanced',dual=False)
    clf = svm.LinearSVC(penalty='l2',dual=False)
    l = list(flow_dict.keys())
    random.shuffle(l)
    x = []
    y = []
    for i in l:
        if i in flow_dict:
            T=i.hour
            x.append([flow_dict[i],np.cos(T*2*np.pi/24),np.sin(T*2*np.pi/24)])
            y.append([binding_hours.get(i,-1)])
    scaler = preprocessing.StandardScaler().fit(x)
    clf2 = calibration.CalibratedClassifierCV(clf,cv=5)
    clf2.fit(scaler.transform(x), y)
    print(clf2.score(scaler.transform(x),y))
    predicted=clf2.predict(scaler.transform(x))
    return clf2,scaler,x,y


cons=[["EMSES",""],["W_DENT",""],["LKPNT",""],["HAMILT","SBRAU"],["ASHERTON-CATARINA",""],["DILLEYSW-COTULAS",""],["LOLITA-FORMOSA",""],["JEWET",""],["SNG-ZEN","DRN"],["BOSQ",""],["COTULAS-DILLEYSW","DDIL"],["BRNSW",""],["LOYOLA",""],["PNHNDL",""],["LISTON",""],["TOKSW-JK_CK","DJEWSNG5"],["HWRDTP-HWRDLN",""],["ODEHV-TROTP",""],["TYLNW",""],["CIBOLO-SCHERT",""],["KLEBERG-KINGSVIL",""],["GILLCR-MCNEIL",""],["BLESSING-LOLITA",""],["BRUNI",""],["192T175",""]]
hour=datetime.timedelta(hours=int(1))
t=datetime.date.today()
d_start = datetime.datetime(t.year,t.month,t.day)+int(1)*datetime.timedelta(days=int(1))
fout = open(os.path.expanduser("~/Activations/activation_"+str(d_start)+".csv"),mode='w')
fout.write("Constraint,Hour,Probability\n")
for con in cons:
    print("test a %s" %con[0])
    b = binding_hours(con[0],con[1])
    a = flow_hourly(con[0],con[1])
    clf,scaler,x,y = ML(b,a)
    probs = list()
    print(d_start)
    for i in range(1,25):
        xx = i
        yy = clf.predict_proba(scaler.transform([[a[d_start+i*hour],cos(i*2*pi/24),sin(i*2*pi/24)]]))[0][1]
        probs.append((xx,yy))
        fout.write(con[0]+","+str(i)+","+str(yy)+"\n")
    Pr=point(probs)    
    Pr.show(title=con[0]+" "+str(d_start.date()))
    s = list()
    pred = clf.predict_proba(scaler.transform(x))
    for i in range(len(x)):
        #print(pred[i])
        xx = x[i][1]
        yy = x[i][2]
        zz = x[i][0]
        if y[i][0]==1:
            s.append(point((xx,zz,yy),color='red'))
        elif pred[i][1]>.12:
            s.append(point((xx,zz,yy),color='green'))
        else:
            s.append(point((xx,zz,yy),color='blue'))
    P = s[0]
    for i in s:
        P+=i
    #P.show(aspect_ratio=[1,1,1])
fout.close()

/projects/sage/sage-7.5/local/bin/python sys.version_info(major=2, minor=7, micro=13, releaselevel='final', serial=0) test a EMSES 0.966890483908 2017-06-10 00:00:00

Error in lines 138-166 Traceback (most recent call last): File "/projects/sage/sage-7.5/local/lib/python2.7/site-packages/smc_sagews/sage_server.py", line 995, in execute exec compile(block+'\n', '', 'single') in namespace, locals File "", line 10, in <module> KeyError: datetime.datetime(2017, 6, 10, 1, 0)