Sunspot number forecasting

In [1]:

    

import research as r
import scipy
import pandas

def fitness_plots(path, treatment):
    D = r.load_files(r.find_files(path, treatment+"_.*/fitness.dat"))
    D['rmse'] = 100.0 / D['max_fitness'] - 1.0

    figure()
    for i,g in D.groupby('trial'):
        plot(g['update'], g['rmse'])
    ylabel('Sum RMSE, t=1:8')
    xlabel('Update')
    
    figure()
    r.quick_ciplot('update','rmse', D)
    ylabel('Sum RMSE, t=1:8')
    xlabel('Update')

    final = D[D['update']==D['update'].max()]
    print "\nDominant fitness:"
    print final.ix[final['rmse'].idxmin()]
    return D

def rmse_boxplots(path, treatment):
    figure()
    D = r.load_files(r.find_files(path, treatment+"_.*/sunspot_test_rmse.dat"))
    B = D.boxplot(column=['tplus1'])
    ylabel("RMSE per horizon")
    
    print "\nMinimum testing RMSE:"
    print D.ix[D['total_rmse'].idxmin()]
    return D

def predictions(path,R):
    T = R.ix[R['total_rmse'].idxmin()]
    path += "/" + T["treatment"] + "_" + T["trial"]
    
    D = r.load_files(r.find_files(path, "sunspot_test_detail.dat"))
    left, width = 0.1, 0.8
    rect1 = [left, 0.25, width, 0.8]
    rect2 = [left, 0.1, width, 0.1]
    
    for i in range(1,2):
        fig=figure()
        ax1 = fig.add_axes(rect1)
        ax2 = fig.add_axes(rect2, sharex=ax1)

        for label in ax1.get_xticklabels():
            label.set_visible(False)
        
        p1 = ax1.plot(D['observed'], color='blue', label='observed')
        p2 = ax1.plot(D['tplus'+str(i)], color='red', label='t+'+str(i))
        ax1.legend()
        ax1.set_ylabel('SSN')
        
        p3 = ax2.vlines(range(0,len(D)), 0, abs(D['observed']-D['tplus'+str(i)]))
        pyplot.locator_params(axis='y', nbins=2)
        ax2.set_ylabel('Error')
        ax2.set_xlabel('Day')
    return D


    

In [19]:

    

def rmse(D):
    obs = D['observed']
    pre = D['tplus1']
    error = obs[1:] - pre[:-1]
    return sqrt(mean(error**2))


    

In [24]:

    

R = rmse_boxplots('/Users/dk/research/var/sunspot/006-daily','ta0') # 2009-2010
R = rmse_boxplots('/Users/dk/research/var/sunspot/006-daily','tb0') # 2009-2010
R = rmse_boxplots('/Users/dk/research/var/sunspot/006-daily','tc0') # 2009-2010
R = rmse_boxplots('/Users/dk/research/var/sunspot/006-daily','td0') # 2009-2010
#P = predictions('/Users/dk/research/var/sunspot/006-daily', R)
#print rmse(D)
#savefig("tplus1_predictions_td0_18.pdf",pad_inches=0, bbox_inches='tight')


    

    




---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
<ipython-input-24-003928a92121> in <module>()
----> 1 R = rmse_boxplots('/Users/dk/research/var/sunspot/006-daily','ta0') # 2009-2010
      2 R = rmse_boxplots('/Users/dk/research/var/sunspot/006-daily','tb0') # 2009-2010
      3 R = rmse_boxplots('/Users/dk/research/var/sunspot/006-daily','tc0') # 2009-2010
      4 R = rmse_boxplots('/Users/dk/research/var/sunspot/006-daily','td0') # 2009-2010
      5 #P = predictions('/Users/dk/research/var/sunspot/006-daily', R)

<ipython-input-18-46ad9ca2be4c> in rmse_boxplots(path, treatment)
     26     figure()
     27     D = r.load_files(r.find_files(path, treatment+"_.*/sunspot_test_rmse.dat"))
---> 28     B = D.boxplot(column=['tplus1'])
     29     ylabel("RMSE per horizon")
     30 

AttributeError: 'NoneType' object has no attribute 'boxplot'

In [4]:

    

D = fitness_plots('/Users/dk/research/var/sunspot/006-daily','ta0') # 1997-1998


    

    




Dominant fitness:
update               100000
mean_generation    3428.238
min_fitness          1.4465
mean_fitness         8.7176
max_fitness         10.9903
treatment               ta0
trial                    29
rmse               8.098933
Name: 22021






    













    

In [5]:

    

D = fitness_plots('/Users/dk/research/var/sunspot/006-daily','tb0') # 2000-2001


    

    




Dominant fitness:
update               100000
mean_generation     3143.95
min_fitness           0.837
mean_fitness         5.3496
max_fitness          6.9769
treatment               tb0
trial                    10
rmse               13.33301
Name: 2001






    













    

In [6]:

    

D = fitness_plots('/Users/dk/research/var/sunspot/006-daily','tc0') # 2003-2004


    

    




Dominant fitness:
update               100000
mean_generation    3506.626
min_fitness          0.9155
mean_fitness          7.135
max_fitness          9.4208
treatment               tc0
trial                     7
rmse                9.61481
Name: 28027






    













    

In [8]:

    

D = fitness_plots('/Users/dk/research/var/sunspot/006-daily','td0') # 2009-2010


    

    




Dominant fitness:
update               100000
mean_generation     3257.11
min_fitness          5.3134
mean_fitness         16.608
max_fitness         20.2666
treatment               td0
trial                    10
rmse               3.934227
Name: 2001