# Import Packages
from SimPEG import *
from IPython.html.widgets import interactive
import numpy as np
%pylab inline

M=5              # RecalL is the number of values we have in m 
n=M+1            # Define n as the N+1 dimension of the matrix   
w=n-1            # Define w and the n-1  dimentsion of the matrix   
s = (M,n)        # Store matrix values 
Av = np.zeros(s) # Create a matrix of zeros of the correct dimensions 
                 # and fill in with elements usin the loop below (note the 1/2 is included in here). 
for i in range(M):
    j=i
    k=i+1
    Av[i,j] = 0.5  
    Av[i,k] = 0.5
print Av

g = lambda x, k, p, q: np.exp(-p*k*x)*np.cos(2*np.pi*q*k*x) # create an anonymous function as above
x=np.linspace(0,1,6)                                          # define the nodes of our x-array
#x = np.array([0., 0.2, 0.4, 0.6, 0.8, 1.])
p = 0.01
q = 0.1
j = np.array([1, 2])
Gn = np.zeros((len(x), len(j)))

for i in range(len(j)):
    f = g(x,i,p,q)
    Gn[:,i] = f
    
print Gn

# Make the delta x array

Deltax = 0.2*np.ones(M) # set x-spacings
V = np.diag(Deltax)     # create diagonal matrix     
print V


m = np.array([0.02, 0.05, 0.09, 0.07, 0.04])
print m

Gc=np.transpose(np.dot(Av, Gn)) # cell center matrix Gc
G = np.dot(Gc,V)                # make G     
d = np.dot(G,m)

#Note: in one line our data d, would be:  d = np.dot(np.dot(np.transpose(np.dot(Av, Gn)), V),m) 
print d


clear

# Begin by creating a ficticious set of model data

M=1000                      # Set number of model parameters 
x=np.linspace(0,1,M+1)      # Define 1D domain on nodes
xc = 0.5*(x[1:] + x[0:-1])  # Define 1D domain on cell centers
m = np.zeros(M)             # preallocate m array 

# Define Gaussian function:
gauss = lambda x, a, mean, SD: a*numpy.exp(-(x-mean)**2./(2.*SD**2.)) # create a Gaussian function

# Choose parameters for Gaussian, evaluate, and store in an array, f.
SD=6                           
mean=0
a=1/(SD*sqrt(2*pi))
x2=np.linspace(-20,20,0.2*M)  
f = gauss(x2,a, mean,SD)

# Define a boxcar function:
box = 0.4*np.ones(0.2*M)

# Insert the Gaussian and Boxcar into m:
m[0.2*M:0.4*M]=box
m[0.6*M:0.8*M]=10*f
   
# Plot      
pylab.plot(xc,m)
#pylab.plot(f)
print size(f)
pylab.xlabel('x')
pylab.ylabel('m(x)')
pylab.title('Model, $m(x)$')


print shape(xc)

# Make the set of kernel functions

g = lambda x, i, p, q: np.exp(-p*i*x)*np.cos(2*np.pi*q*i*x) # create an anonymous function as immediately above
p = 0.01      # Set values for p, q
q = 0.15
N = 20        # specify number of output data
Gn = np.zeros((M+1,N))

for i in range(N):
    f = g(x,i,p,q)
    Gn[:,i] = f
    
# Plot  
pylab.plot(x,Gn)
pylab.xlabel('x')
pylab.ylabel('g(x)')
pylab.title('Kernel functions, $g(x)$')

# Make Averaging Matrix
n=M+1            # Define n as the N+1 dimension of the matrix   
w=n-1            # Define w and the n-1  dimentsion of the matrix   
s = (M,n)        # Store matrix values 
Av = np.zeros(s) # Create a matrix of zeros of the correct dimensions 
                 # and fill in with elements usin the loop below (note the 1/2 is included in here). 
for i in range(M):
    j=i
    k=i+1
    Av[i,j] = 0.5  
    Av[i,k] = 0.5
print shape(Av)

# Plot  
pylab.title('Subset of averaging matrix')
pyplot.imshow(Av[:5,:6], cmap='Greys',  interpolation='nearest')
print(Av[:5,:6])

# make the Volume, "delta x" array

Deltax = 0.2*np.ones(M) # set x-spacings
V = np.diag(Deltax)     # create diagonal matrix     
print shape(V)


d = np.dot(np.dot(np.transpose(np.dot(Av, Gn)), V),m)   # Create synthetic data
xd=np.linspace(0,1,len(d))                              # make x array for data 

# Plot

pylab.plot(xd,d)
pylab.xlabel('')
pylab.ylabel('d')
pylab.title('Synthetic Data $d$')

g = lambda x, k, p, q: np.exp(-p*k*x)*np.cos(2*np.pi*q*k*x)  # I redfine our kernel functions, as above.

M = Mesh.TensorMesh([1000])                                  # create a tensor mesh
Glist = [g(M.vectorNx,i+1,p,q) for i in range(20)]           # make kernel functions on nodes with the following lines
Glist = []
for i in range(20):
    gi = g(M.vectorNx,i+1,p,q)
    Glist.append(gi)
Gnodes = np.vstack(Glist).T                            
Gcells = M.aveF2CC * Gnodes                                  # evaluate kernel functions on cell centers directly
V = Utils.sdiag(M.vol)                                       # build volume array  
d = np.dot(Gcells.T, V * m)                                  # create data  
xd=np.linspace(0,1,len(d)) 

# Plot
pylab.plot(xd,d)
pylab.xlabel('')
pylab.ylabel('d')
pylab.title('Synthetic Data $d$')