<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Make X not baD!
Load X and Y from getting satellite channels 1 and 6 data, see satellite-to-sensor-model.ipynb for more info.
Step1:  Random Forest seems to be giving the best results, so we'll stick with that for now
Step2:  Maybe I should only use the AOD values since the sensor values are sporadic and don't add as much to the data anyway.  Also, I could remove the rows were the sensor readings are zero.  I could look through the Y rows and make a mask to apply to both X and Y rows.
Step3:  Let's try Ridge again.
Step4:  Not too good. Or at least not an improvement over itself even.
I may have to do some feature engineering for the images.  Averaging, filtering or both. I wonder if I do some ratios of channels?  Sort of like the "greenness" parameter that farmers use.  Also, maybe +/- 3 hours is too long for searching for a satellite image. (Reduce to +/- 1 hour?)
Step5:  So, from the plot above, there are values above 0, even though my Y_train and Y_test are all below zero.
Step6:  Try random forest with new X_hist
Step7:  The score has improved (0.2) but the most remarkable is that all the predicted values are now below 0, as they should be.
Let's try something Extra Trees!
Step8:  Not much better than random forest regressor, but does offer multicore support which speeds things up significantly. Edit
Step9:  Wow, even better, but is it good enough?  If I include ratios between channels, spatial and other information will be implicitly collected.
Downsample CH1 and divide by CH6
Step10:  Try it all again with just X, not X_reduced?
Step11:  After pickling the sensor to power model, use it to predict power from satellite data!
Step12:  Seems like dimensions are working out.  Probably will need to pickle the sat to sensor model too and just feed in data start to finish and see how it does.
Step13:  Nice, now we have predictions for every datetime!  Now we have to check this.
Step14:  Tres beau, call it a day!  That's a pretty graph!
Step15:  I shut down the nb, so restart!
<ASSISTANT_TASK:>
Python Code:
import numpy as np
import matplotlib.pyplot as plt
from data_helper_functions import *
from IPython.display import display
pd.options.display.max_columns = 999
%matplotlib inline
with np.load('data/X.npz') as data: #old X, don't use, start at "Now with all channels..."
    X = data['X'] 
with np.load('data/Y.npz') as data: #old Y, don't use
    Y = data['Y']
print X.shape
print Y.shape
Explanation: Make X not baD!
Load X and Y from getting satellite channels 1 and 6 data, see satellite-to-sensor-model.ipynb for more info.
End of explanation
from sklearn.ensemble import RandomForestRegressor
from sklearn.cross_validation import train_test_split
rfr = RandomForestRegressor(oob_score=True)
X_train, X_test, Y_train, Y_test = train_test_split(X, Y, test_size = 0.2, random_state = 14)
rfr.fit(X_train,Y_train)
rfr.score(X_test,Y_test)
Y_pred = rfr.predict(X_test)
from random import randint
val = randint(0,508)
print Y_pred[val]
print Y_test[val]
Explanation: Random Forest seems to be giving the best results, so we'll stick with that for now
End of explanation
mask = []
for i,row in enumerate(Y):
    if row[0] == 0:
        mask.append(False)
    else:
        mask.append(True)    
mask = np.array(mask)
X_reduced = X[mask]
Y_reduced = Y[mask][:,-7::]
print X_reduced.shape
print Y_reduced.shape
rfr = RandomForestRegressor(oob_score=True)
X_train, X_test, Y_train, Y_test = train_test_split(X_reduced, Y_reduced, test_size = 0.3, random_state = 14)
rfr.fit(X_train,Y_train)
print rfr.score(X_test,Y_test)
print rfr.oob_score_
Y_pred = rfr.predict(X_test)
from random import randint
val = randint(0,Y_pred.shape[0])
print Y_pred[val]
print Y_test[val]
Explanation: Maybe I should only use the AOD values since the sensor values are sporadic and don't add as much to the data anyway.  Also, I could remove the rows were the sensor readings are zero.  I could look through the Y rows and make a mask to apply to both X and Y rows.
End of explanation
from sklearn.linear_model import Ridge
ridge = Ridge(solver = 'auto')
X_train, X_test, Y_train, Y_test = train_test_split(X_reduced, Y_reduced, test_size = 0.5, random_state = 14)
ridge.fit(X_train,Y_train)
ridge.score(X_test,Y_test)
Y_pred = ridge.predict(X_test)
from random import randint
val = randint(0,Y_pred.shape[0]-1)
print Y_pred[val]
print Y_test[val]
Explanation: Let's try Ridge again.
End of explanation
plt.plot(Y_pred,'go');
Explanation: Not too good. Or at least not an improvement over itself even.
I may have to do some feature engineering for the images.  Averaging, filtering or both. I wonder if I do some ratios of channels?  Sort of like the "greenness" parameter that farmers use.  Also, maybe +/- 3 hours is too long for searching for a satellite image. (Reduce to +/- 1 hour?)
End of explanation
val = randint(0,2301)
_ = plt.hist(X_reduced[:,0:1972][val], alpha = 0.5, normed=True, bins=25, label='visible',range=(0,25000)) #visible light histogram
_ = plt.hist(X_reduced[:,1973:2476][val], alpha = 0.5, normed=True, bins=25, label='IR',range=(0,25000))
plt.legend(loc='upper right');
X_hist = []
for i in xrange(X_reduced.shape[0]):
    hist1, _ = np.histogram(X_reduced[:,0:1972][i], density=True, bins=25, range=(0,25000)) #vis
    hist2, _ = np.histogram(X_reduced[:,1973:2476][i], density=True, bins=25, range=(0,25000)) #IR
    X_hist.append(np.hstack((hist1,hist2)))
X_hist = np.array(X_hist)
Explanation: So, from the plot above, there are values above 0, even though my Y_train and Y_test are all below zero.
End of explanation
rfr = RandomForestRegressor(oob_score=True)
X_train, X_test, Y_train, Y_test = train_test_split(X_hist, Y_reduced, test_size = 0.3, random_state = 14)
rfr.fit(X_train,Y_train)
print rfr.score(X_test,Y_test)
print rfr.oob_score_
from random import randint
val = randint(0,Y_test.shape[0])
print Y_pred[val]
print Y_test[val]
plt.plot(Y_pred,'go');
Explanation: Try random forest with new X_hist:
End of explanation
X_hist = []
for i in xrange(X_reduced.shape[0]):
    hist1, _ = np.histogram(X_reduced[:,0:1972][i], density=True, bins=30, range=(0,25000)) #vis
    hist2, _ = np.histogram(X_reduced[:,1973:2476][i], density=True, bins=30, range=(0,25000)) #IR
    X_hist.append(np.hstack((hist1,hist2)))
X_hist = np.array(X_hist)
from sklearn.ensemble import ExtraTreesRegressor
etr = ExtraTreesRegressor(oob_score=True, bootstrap=True, n_jobs=-1, n_estimators=100) #njobs uses all cores!
X_train, X_test, Y_train, Y_test = train_test_split(X_hist, Y_reduced, test_size = 0.3, random_state = 14)
etr.fit(X_train,Y_train)
print etr.score(X_test,Y_test)
print etr.oob_score_
from random import randint
val = randint(0,Y_test.shape[0])
print Y_pred[val]
print Y_test[val]
plt.plot(Y_pred,'go');
Explanation: The score has improved (0.2) but the most remarkable is that all the predicted values are now below 0, as they should be.
Let's try something Extra Trees!
End of explanation
with np.load('data/X_all_channels.npz') as data:
    X = data['X']
    
with np.load('data/Y_all_channels.npz') as data:
    Y = data['Y']
print X.shape
print Y.shape
mask = []
for i,row in enumerate(Y):
    if row[0] == 0:
        mask.append(False)
    else:
        mask.append(True)    
mask = np.array(mask)
X_reduced = X[mask]
Y_reduced = Y[mask][:,-7::]
print X_reduced.shape
print Y_reduced.shape
val = randint(0,2301)
_ = plt.hist(X_reduced[:,0:1972][val], alpha = 0.2, normed=True, bins=30, label='CH1',range=(0,25000)) #visible light histogram
_ = plt.hist(X_reduced[:,1973:2476][val], alpha = 0.2, normed=True, bins=30, label='CH2',range=(0,25000))
_ = plt.hist(X_reduced[:,2477:2980][val], alpha = 0.2, normed=True, bins=30, label='CH3',range=(0,25000))
_ = plt.hist(X_reduced[:,2981:3484][val], alpha = 0.2, normed=True, bins=30, label='CH4',range=(0,25000))
_ = plt.hist(X_reduced[:,3484:3988][val], alpha = 0.2, normed=True, bins=30, label='CH6',range=(0,25000))
plt.legend(loc='upper right');
X_hist = []
bins = 20
for i in xrange(X_reduced.shape[0]):
    hist1, _ = np.histogram(X_reduced[:,0:1972][i], density=True, bins=bins, range=(0,25000))
    hist2, _ = np.histogram(X_reduced[:,1973:2476][i], density=True, bins=bins, range=(0,25000))
    hist3, _ = np.histogram(X_reduced[:,2477:2980][i], density=True, bins=bins, range=(0,25000))
    hist4, _ = np.histogram(X_reduced[:,2981:3484][i], density=True, bins=bins, range=(0,25000))
    hist5, _ = np.histogram(X_reduced[:,3484:3988][i], density=True, bins=bins, range=(0,25000))
    X_hist.append(np.hstack((hist1,hist2,hist3,hist4,hist5)))
X_hist = np.array(X_hist)
from sklearn.ensemble import ExtraTreesRegressor
etr = ExtraTreesRegressor(
    oob_score=True, bootstrap=True, n_jobs=-1, n_estimators=1000
                         ) #nj_obs uses all cores!
X_train, X_test, Y_train, Y_test = train_test_split(X_hist, Y_reduced, test_size = 0.25, random_state = 12)
etr.fit(X_train,Y_train)
print etr.score(X_test,Y_test)
print etr.oob_score_
from random import randint
val = randint(0,Y_test.shape[0])
print Y_pred[val]
print Y_test[val]
plt.plot(Y_pred,'go');
Explanation: Not much better than random forest regressor, but does offer multicore support which speeds things up significantly. Edit: actually, this is better when increasing the number of n_estimators.
Now with all channels! (Skip to "try it again..." for model used)
End of explanation
from scipy.ndimage import zoom
from __future__ import division
X_reduced_ratio_1_2 = []
for i in xrange(X_reduced.shape[0]):
    CH1 = zoom(X_reduced[:,0:1972][i].reshape((29,68)),zoom=(0.48, 0.53), order=5)
    CH2 = X_reduced[:,1972:2476][i].reshape((14,36))
    X_reduced_ratio_1_2.append(25000* CH2  / (CH1 + CH2) )
X_reduced_ratio_1_2 = np.array(X_reduced_ratio_1_2)
X_reduced_ratio_1_6 = []
for i in xrange(X_reduced.shape[0]):
    CH1 = zoom(X_reduced[:,0:1972][i].reshape((29,68)),zoom=(0.48, 0.53), order=5)
    CH6 = X_reduced[:,3484:3988][i].reshape((14,36))
    X_reduced_ratio_1_6.append(25000* CH6 / (CH1 + CH6) )
X_reduced_ratio_1_6 = np.array(X_reduced_ratio_1_6)
X_reduced_ratio_2_6 = []
for i in xrange(X_reduced.shape[0]):
    CH2 = X_reduced[:,1972:2476][i].reshape((14,36))
    CH6 = X_reduced[:,3484:3988][i].reshape((14,36))
    X_reduced_ratio_2_6.append(25000* CH6 / (CH2 + CH6) )
X_reduced_ratio_2_6 = np.array(X_reduced_ratio_2_6)
val = randint(0,2301)
_ = plt.hist(X_reduced[:,0:1972][val], alpha = 0.2, normed=True, bins=30, label='CH1',range=(0,25000)) #visible light histogram
_ = plt.hist(X_reduced[:,1973:2476][val], alpha = 0.2, normed=True, bins=30, label='CH2',range=(0,25000))
_ = plt.hist(X_reduced[:,2477:2980][val], alpha = 0.2, normed=True, bins=30, label='CH3',range=(0,25000))
_ = plt.hist(X_reduced[:,2981:3484][val], alpha = 0.2, normed=True, bins=30, label='CH4',range=(0,25000))
_ = plt.hist(X_reduced[:,3484:3988][val], alpha = 0.2, normed=True, bins=30, label='CH6',range=(0,25000))
_ = plt.hist(np.ravel(X_reduced_ratio_1_2[val]), alpha = 0.2, normed=True, bins=30, label='CH2/CH1',range=(0,25000))
_ = plt.hist(np.ravel(X_reduced_ratio_1_6[val]), alpha = 0.2, normed=True, bins=30, label='CH6/CH1',range=(0,25000))
_ = plt.hist(np.ravel(X_reduced_ratio_2_6[val]), alpha = 0.2, normed=True, bins=30, label='CH6/CH2',range=(0,25000))
plt.legend(loc='upper right');
X_hist = []
bins = 30
for i in xrange(X_reduced.shape[0]):
    hist1, _ = np.histogram(X_reduced[:,0:1972][i], density=True, bins=bins, range=(0,25000))
    hist2, _ = np.histogram(X_reduced[:,1972:2476][i], density=True, bins=bins, range=(0,25000))
    hist3, _ = np.histogram(X_reduced[:,2476:2980][i], density=True, bins=bins, range=(0,25000))
    hist4, _ = np.histogram(X_reduced[:,2980:3484][i], density=True, bins=bins, range=(0,25000))
    hist5, _ = np.histogram(X_reduced[:,3484:3988][i], density=True, bins=bins, range=(0,25000))
    hist6, _ = np.histogram(np.ravel(X_reduced_ratio_1_2[i]), density=True, bins=bins, range=(0,25000))
    hist7, _ = np.histogram(np.ravel(X_reduced_ratio_1_6[i]), density=True, bins=bins, range=(0,25000))
    hist8, _ = np.histogram(np.ravel(X_reduced_ratio_2_6[i]), density=True, bins=bins, range=(0,25000))
    X_hist.append(np.hstack((hist1,hist2,hist3,hist4,hist5,hist6,hist7,hist8)))
X_hist = np.array(X_hist)
from sklearn.ensemble import ExtraTreesRegressor
etr = ExtraTreesRegressor(oob_score=True, bootstrap=True,
                          n_jobs=-1, n_estimators=1000) #nj_obs uses all cores!
X_train, X_test, Y_train, Y_test = train_test_split(X_hist, Y_reduced, test_size = 0.2, random_state = 12)
etr.fit(X_train,Y_train)
print etr.score(X_test,Y_test)
print etr.oob_score_
from random import randint
val = randint(0,Y_test.shape[0])
print Y_pred[val]
print Y_test[val]
Explanation: Wow, even better, but is it good enough?  If I include ratios between channels, spatial and other information will be implicitly collected.
Downsample CH1 and divide by CH6:
End of explanation
from scipy.ndimage import zoom
from __future__ import division
from random import randint
X_ratio_1_2 = []
for i in xrange(X.shape[0]):
    CH1 = zoom(X[:,0:1972][i].reshape((29,68)),zoom=(0.48, 0.53), order=5)
    CH2 = X[:,1972:2476][i].reshape((14,36))
    X_ratio_1_2.append(25000* (CH2)  / (CH1 + CH2+1.0) )
X_ratio_1_2 = np.array(X_ratio_1_2)
X_ratio_1_6 = []
for i in xrange(X.shape[0]):
    CH1 = zoom(X[:,0:1972][i].reshape((29,68)),zoom=(0.48, 0.53), order=5)
    CH6 = X[:,3484:3988][i].reshape((14,36))
    X_ratio_1_6.append(25000* CH6 / (CH1 + CH6 + 0.1) )
X_ratio_1_6 = np.array(X_ratio_1_6)
X_ratio_2_6 = []
for i in xrange(X.shape[0]):
    CH2 = X[:,1972:2476][i].reshape((14,36))
    CH6 = X[:,3484:3988][i].reshape((14,36))
    X_ratio_2_6.append(25000* CH6 / (CH2 + CH6 + 0.1) )
X_ratio_2_6 = np.array(X_ratio_2_6)
val = randint(0,2301)
_ = plt.hist(X[:,0:1972][val], alpha = 0.2, normed=True, bins=30, label='CH1',range=(0,25000)) #visible light histogram
_ = plt.hist(X[:,1973:2476][val], alpha = 0.2, normed=True, bins=30, label='CH2',range=(0,25000))
_ = plt.hist(X[:,2477:2980][val], alpha = 0.2, normed=True, bins=30, label='CH3',range=(0,25000))
_ = plt.hist(X[:,2981:3484][val], alpha = 0.2, normed=True, bins=30, label='CH4',range=(0,25000))
_ = plt.hist(X[:,3484:3988][val], alpha = 0.2, normed=True, bins=30, label='CH6',range=(0,25000))
_ = plt.hist(np.ravel(X_ratio_1_2[val]), alpha = 0.2, normed=True, bins=30, label='CH2/CH1',range=(0,25000))
_ = plt.hist(np.ravel(X_ratio_1_6[val]), alpha = 0.2, normed=True, bins=30, label='CH6/CH1',range=(0,25000))
_ = plt.hist(np.ravel(X_ratio_2_6[val]), alpha = 0.2, normed=True, bins=30, label='CH6/CH2',range=(0,25000))
plt.legend(loc='upper right');
import pandas as pd
X_hist = []
bins = 25
for i in xrange(X.shape[0]):
    myval1 = pd.DataFrame(np.ravel(X_ratio_1_2[i])).fillna(np.mean).values.flatten();
    myval2 = pd.DataFrame(np.ravel(X_ratio_1_6[i])).fillna(np.mean).values.flatten();
    myval3 = pd.DataFrame(np.ravel(X_ratio_2_6[i])).fillna(np.mean).values.flatten();
    
    hist1, _ = np.histogram(X[:,0:1972][i], density=True, bins=bins, range=(0,25000))
    hist2, _ = np.histogram(X[:,1972:2476][i], density=True, bins=bins, range=(0,25000))
    hist3, _ = np.histogram(X[:,2476:2980][i], density=True, bins=bins, range=(0,25000))
    hist4, _ = np.histogram(X[:,2980:3484][i], density=True, bins=bins, range=(0,25000))
    hist5, _ = np.histogram(X[:,3484:3988][i], density=True, bins=bins, range=(0,25000))
    hist6, _ = np.histogram( myval1 , density=True, bins=bins, range=(0,25000) )
    hist7, _ = np.histogram( myval2 , density=True, bins=bins, range=(0,25000))
    hist8, _ = np.histogram( myval3, density=True, bins=bins, range=(0,25000))
    X_hist.append(np.hstack((hist1,hist2,hist3,hist4,hist5,hist6,hist7,hist8)))
X_hist = np.array(X_hist)
from sklearn.ensemble import ExtraTreesRegressor
from sklearn.cross_validation import train_test_split
etr = ExtraTreesRegressor(oob_score=True, bootstrap=True,
                          n_jobs=-1, n_estimators=500) #nj_obs uses all cores!
X_train, X_test, Y_train, Y_test = train_test_split(X_hist, Y, test_size = 0.2, random_state = 12)
etr.fit(X_train,Y_train) #pickle it!
from sklearn.externals import joblib
joblib.dump(etr, 'webapp/solarApp/models/sat-to-sensor-model/sat-to-sensor-model.pkl')
sat_to_sensor_model = joblib.load('webapp/solarApp/models/sat-to-sensor-model/sat-to-sensor-model.pkl')
sat_to_sensor_model.predict(X_hist).shape
X_hist.shape
print etr.score(X_test,Y_test)
print etr.oob_score_
Y_pred = etr.predict(X_test)
from random import randint
val = randint(0,Y_test.shape[0])
print Y_pred[val]
print Y_test[val]
Explanation: Try it all again with just X, not X_reduced?
End of explanation
from sklearn.externals import joblib
etr2 = joblib.load('data/sensor-to-power-model/sensor-to-power-model.pkl')
etr2.predict(Y).shape #power predictions
with np.load('data/y.npz') as data:
    y = data['y']
y_pred2 = etr2.predict(Y_pred) #the predicted
y_pred2.shape
y.shape
Explanation: After pickling the sensor to power model, use it to predict power from satellite data!
End of explanation
with np.load('data/good_times.npz') as data:
    good_times = data['good_times']
print good_times.shape
print etr2.predict(Y).shape
Explanation: Seems like dimensions are working out.  Probably will need to pickle the sat to sensor model too and just feed in data start to finish and see how it does.
End of explanation
from datetime import datetime, timedelta, time
pvoutput_filefolder = 'data/pvoutput/pvoutput6months/'
datetime_index = 1747
print 'Predicted power for ' + str(good_times[datetime_index]) + \
      ' is ' + str(etr2.predict(Y[datetime_index])[0]) + ' W'
#pvoutput data
desired_datetime = good_times[datetime_index]
desired_date = (desired_datetime - timedelta(hours=6)).date() #make sure correct date
desired_date = datetime.combine(desired_date, time.min) #get into datetime format
pvoutput_filename = find_file_from_date(desired_date, pvoutput_filefolder)
df_pvoutput = return_pvoutput_data(pvoutput_filename, pvoutput_filefolder)
try:
    print "True power: " + df_pvoutput[df_pvoutput.index == desired_datetime].values[0][0].astype(str)
except:
    print "0"
from datetime import datetime, timedelta, time
pvoutput_filefolder = 'data/pvoutput/pvoutput6months/'
predicted_powers = []
true_powers = []
for datetime_index in range(len(good_times)):
    predicted_powers.append(str(etr2.predict(Y[datetime_index])[0]) )
    #pvoutput data
    desired_datetime = good_times[datetime_index]
    desired_date = (desired_datetime - timedelta(hours=6)).date() #make sure correct date
    desired_date = datetime.combine(desired_date, time.min) #get into datetime format
    pvoutput_filename = find_file_from_date(desired_date, pvoutput_filefolder)
    df_pvoutput = return_pvoutput_data(pvoutput_filename, pvoutput_filefolder)
    try:
        true_powers.append(df_pvoutput[df_pvoutput.index == desired_datetime].values[0][0].astype(str))
    except:
        true_powers.append(0)
true_powers = np.array(true_powers).astype(float)
predicted_powers = np.array(predicted_powers).astype(float)
error_add = []
for datetime_index in range(len(good_times)):
        error_add.append(
            np.abs(true_powers[datetime_index]-predicted_powers[datetime_index]) \
            / true_powers[datetime_index]
            )
error_sum = 0
for error in error_add:
    if (np.isnan(error) or np.isinf(error)):
        pass
    else:
        error_sum += error
        
error_sum/len(good_times)
x = np.linspace(0,13000,100)
y = 1.0 * x
plt.figure(figsize=(8,8))
plt.plot(true_powers,predicted_powers,'go', alpha=0.25)
plt.plot(x,y,'r',linewidth = 5, alpha = 0.5)
plt.xlabel('True Power (W)', fontsize = 16)
plt.xticks(size=16)
plt.ylabel('Predicted Power (W)', fontsize = 16)
plt.yticks(size=16)
Explanation: Nice, now we have predictions for every datetime!  Now we have to check this.
End of explanation
np.savez_compressed('data/true_powers.npz',true_powers=true_powers) #save
np.savez_compressed('data/predicted_powers.npz',predicted_powers=predicted_powers) #save
Explanation: Tres beau, call it a day!  That's a pretty graph!
End of explanation
import numpy as np
with np.load('data/true_powers.npz') as data:
    true_powers = data['true_powers']
with np.load('data/predicted_powers.npz') as data:
    predicted_powers = data['predicted_powers']
import matplotlib.pyplot as plt
%matplotlib inline
import statsmodels.api as sm
# Fit and summarize OLS model
mod = sm.OLS(predicted_powers,true_powers)
res = mod.fit()
print res.summary()
Explanation: I shut down the nb, so restart!
End of explanation
<END_TASK>