In [ ]:

    

from rpy2.rinterface import R_VERSION_BUILD
print(R_VERSION_BUILD)


    

In [ ]:

    

# Load in the r magic
%load_ext rpy2.ipython


    

In [ ]:

    

# This magic command makes only this line to an R interpreted cell!
%R X=c(1,4,5,7); sd(X); mean(X)


    

In [ ]:

    

%%R #This magic command make this cell an R interpreted cell!
X=c(1,4,5,7);
sd(X); 
mean(X)


    

In [ ]:

    

# We need ggplot2 (in R)
%R require(ggplot2)
# Load in the pandas library
import pandas as pd 
# Make a pandas DataFrame
df = pd.DataFrame({'Alphabet': ['a', 'b', 'c', 'd','e', 'f', 'g', 'h','i'],
                   'A': [4, 3, 5, 2, 1, 7, 7, 5, 9],
                   'B': [0, 4, 3, 6, 7, 10,11, 9, 13],
                   'C': [1, 2, 3, 1, 2, 3, 1, 2, 3]})


    

In [ ]:

    

%%R -i df # %%R MUST be the first line 
# Take the name of input variable df and assign it to an R variable of the same name
print(df)
# Plot the DataFrame df
p <- ggplot(data=df) + geom_point(aes(x=A, y=B, color=C))
print(p)


    

In [ ]:

    

import numpy as np
import matplotlib.pyplot as plt
import control as con
%pylab inline

K = 1
d = 0.5
T = 10
delay = 10

a0 = 1
a1 = (2 * d * T) #16
a2 = (T**2) #100
b0 = K

tf_1 = con.matlab.tf(K, [a2, a1, a0])
#print tf_1
ss_1a = con.matlab.tf2ss(tf_1)
#print ss_1a

d_num, d_den = con.pade(delay, 1)
tf_delay = con.tf(d_num, d_den)
ss_delay = con.matlab.tf2ss(con.series(tf_delay, tf_1))
print ss_delay

#print con.matlab.tf2ss(ss_delay)
d_yout, d_T = con.matlab.step(ss_delay)
yout, T = con.matlab.step(tf_1) # step without delay
xin = np.ones(len(yout))



plt.plot(d_T, d_yout, 'r-', label='poly_est')
plt.plot(np.add(d_T, delay), yout, 'g-', label='idealized') #delay in timeaxis!


    

In [ ]:

    

# Don't know how to import more than one per line...
%R -i T
%R -i d_yout
%R -i xin


    

In [ ]:

    

len(xin) == len(d_yout) == len(T)


    

In [ ]:

    

%%R -o sys
library(sysid)
require(ggplot2)

data_tf <- idframe(output = d_yout, input = xin)
sys = arx(data_tf, c(3,3,5))
#ls()
#names(sys)
#sys['fitted.values']
plot(T, d_yout)
lines(T, predict(sys, data_tf, nahead=1))
#print(sys)
names(sys)
print(sys)


    

In [ ]:

    

um = np.matrix(xin) # input measurements
ym = np.matrix(d_yout) # output measurements

Rexpt = None
fexpt = None
K= 4

for k in xrange(1,K):
#    print k
    #print ym[k-1]
    vphi = np.matrix([-ym.item(k-1), -ym.item(k-2), -ym.item(k-3), um.item(k-2), um.item(k-3)]).transpose();
    N = float(len(vphi))
    if Rexpt is not None:
        Rexpt = Rexpt + (1/N)*vphi*vphi.transpose();
        fexpt = fexpt + (1/N)*vphi*ym.item(k);
    else:
        Rexpt = (1/N)*vphi*vphi.transpose();
        fexpt = (1/N)*vphi*ym.item(k);

#theta = Rexpt/fexpt # least-squares fit
theta = np.linalg.lstsq(Rexpt, fexpt)
#print theta[0]
Bvec = [0, 0, theta[0].item(3), theta[0].item(4)];
Avec = [1, theta[0].item(0), theta[0].item(1), theta[0].item(2)];

Ghat = con.tf(Bvec,Avec);
est_yout, est_T = con.matlab.step(Ghat)
plt.plot(est_T, est_yout, 'r-', label='tf_est')


    

In [ ]:

    

# Y = B/A * U
# B = 0.4559*q - 0.41-0.051
# A = 1 * q - 0.9521-0.0212
num = [0.4559, -0.41+0.051]
den = [1, -0.9521+0.0212]
est_tf = con.tf(num, den)

est_yout, est_T = con.matlab.step(est_tf)
plt.plot(est_T, est_yout, 'r-', label='tf_est')


    

In [ ]:

    




    

In [ ]:

    




    

In [ ]:

    




    

In [ ]: