"""
Variogram estimation and modelling functions
"""

import numpy as np
from scipy.spatial.distance import pdist
from scipy.optimize import fmin
import random

def pvariance(A):
    """ For data `A`, return the pairwise variance (squared differences). This
    will need to be Cythonized. """
    n = len(A)
    V = np.empty(sum(range(n)))
    ind = 0
    for i in range(len(A)):
        for j in range(i+1, n):
            V[ind] = (A[i] - A[j])**2
            ind += 1
    return V

def fit_model(model, p, lags, variance):
    """ Fit callable *model* with guess *p* to an experimental variogram
    given by *lags*, *variance*.

    *model* should return an estimator function that takes *lags* as an
    argument
    """
    obj = lambda p: model(p)(lags)
    res = fmin(obj, p)
    pbest = res[0]
    return

def estimate_vario(mp, npoints=150, max_dist=None, interval=None):
    """ Given a Multipoint input, estimate a spatial variogram. By default, a
    variogram is quickly estimated. If npoints is None, then the full variogram
    is calculated.

    Parameters:
    -----------
    points : a Multipoint instance
    npoints : number of values use for variogram appromixation. If npoints is
              None, the full variogram is calculated.
    max_dist : the maximum lag
    interval : the lag interval

    Returns:
    --------
    lags : ndarray
    variogram : ndarray
    """
    if npoints is None:
        npoints = len(mp)

    irand = random.sample(np.arange(len(mp)), npoints)
    verts = mp.get_vertices()[irand]
    z = mp.get_data()[irand]
    dist = pdist(verts)
    vari = pvariance(z)

    if max_dist is None:
        max_dist = dist.max()

    if interval is None:
        interval = max_dist / 10.0

    lags = np.arange(0, max_dist+interval, interval)
    sigma_variance = np.empty_like(lags)

    sigma_variance = [np.mean(vari[dist<=lag]) for lag in lags]
    return lags, sigma_variance

def sph_func(a, h):
    """ Spherical model function for range *a* and lag *h*. """
    if h < a:
        val = 1.5*h/a - 0.5 * (h/a)**3
    else:
        val = 1
    return val

def gau_func(a, h):
    """ Gaussian model function for range *a* and lag *h*. """
    val = 1.0 - np.exp(-(h**2 / a**2))
    return val

@@ -134,20 +134,24 @@
  import 
-f
 min
+imize
 %0Aimport 
@@ -744,16 +744,21 @@
 mbda p: 
+sum((
 model(p)
@@ -763,16 +763,32 @@
 p)(lags)
+ - variance)**2)
 %0A    res
@@ -794,12 +794,16 @@
 s = 
-f
 min
+imize
 (obj
@@ -815,33 +815,20 @@
     
-pbest = res%5B0%5D%0A    
 return
+ res.x
 %0A%0Ade
@@ -1769,17 +1769,24 @@
 .arange(
-0
+interval
 , max_di
@@ -1849,17 +1849,82 @@
 e(lags)%0A
-%0A
+    for i, lag in enumerate(lags):%0A        band = dist %3C= lag%0A    
     sigm
@@ -1937,12 +1937,14 @@
 ance
+%5Bi%5D
  = 
-%5B
 np.m
@@ -1956,36 +1956,14 @@
 ari%5B
-dist%3C=lag%5D) for lag in lags%5D
+band%5D)
 %0A   
@@ -2084,31 +2084,29 @@
     
-if h %3C a:%0A        val =
+val = np.where(h %3C a,
  1.5
@@ -2130,34 +2130,14 @@
 )**3
-%0A    else:%0A        val = 1
+, 1.0)
 %0A   

fix divide by zero