def run_once():
    simplex = Monitor()
    solver = fmin(2)
    solver.SetRandomInitialPoints([0,0],[2,2])
    solver.SetGenerationMonitor(simplex)
    solver.Solve(Corana2, termination=CRT())
    sol = solver.Solution()
    
    for x in simplex.x:
        sam.putarray('x',x)
        sam.eval("plot(x([1,2,3,1],1),x([1,2,3,1],2),'w-')")

    def plot(self,contours=None):
        """Plots the isosurfaces for the field data, at corresponding contour values."""
        import ExcitationSlicer
        path = ExcitationSlicer.__file__
        path=path.strip('__init__.pyc')

        
        if contours is None:
            contours = self._contvals
        
        import sam
        import numpy as np

        try:
            contourvalue = contours[0]
        except:
            raise ValueError, "contourvalues is empty."

        #we first pass the data into Matlab as an array:
        array = np.array(self._getArray())
        (dim0, dim1, dim2) = array.shape
        q0 = np.arange(1, dim0 + 1)
        q1 = np.arange(1, dim1 + 1)
        q2 = np.arange(1, dim2 + 1)

        sam.put('dim0', [dim0])
        sam.put('dim1', [dim1])
        sam.put('dim2', [dim2])

        sam.put('q0', q0)
        sam.put('q1', q1)
        sam.put('q2', q2)
        sam.eval("[q0i, q1i, q2i] = meshgrid(q0, q1, q2);")

        # this causes a segmentation fault for an array of 20x20x20 doubles:
        #sam.putarray('array', array)
        # here is a temporary ugly(!) fix:
        array.shape=(dim0*dim1*dim2,)
        sam.put('array', array)
        sam.eval("earray = reshape(array', dim0, dim1, dim2)")
        # end of fix

        # pass the energy-value for which to draw the isosurface into matlab"
        sam.put('contourvalue', [contourvalue])

        # We use an M-file that gets called when we evaluate the name
        # of the corresponding function in Matlab via SAM

        print "path: ", path
        sam.eval("cd('"+path+"')")
        sam.eval("MlabPhonIsoSurfacePlotter(q0i,q1i,q2i, earray, contourvalue, 'blue')")
        waitforentry=raw_input("Press any key.")
        sam.eval("close")
        return 0

def run_once(x0,x1):
    simplex = Monitor()
    xinit = [x0, x1]

    solver = fmin(len(xinit))
    solver.SetInitialPoints(xinit)
    solver.SetGenerationMonitor(simplex)
    solver.Solve(rosen, termination=CRT())
    sol = solver.Solution()
    
    for x in simplex.x:
        sam.putarray('x',x)
        sam.eval("plot(x([1,2,3,1],1),x([1,2,3,1],2),'w-')")

def run_once_xv():
    simplex = Monitor()
    y1 = y0*random.uniform(0.5,1.5)
    z1 = z0*random.uniform(0.5,1.5)
    xinit = [random.uniform(x0-40,x0+40), y1, z1, random.uniform(v0-0.1,v0+0.1)]

    solver = fmin(len(xinit))
    solver.SetInitialPoints(xinit)
    solver.SetGenerationMonitor(simplex)
    solver.Solve(cost_function, termination=CRT())
    sol = solver.Solution()
    print(sol)

    for x in simplex.x:
        sam.putarray('x',x)
        sam.eval("plot(x([1,2,3,1],1),x([1,2,3,1],2),'w-','LineWidth',2)")
    return sol

def draw_contour_xv():
    import numpy
    x, y = mgrid[-40:40:0.5, -0.1:0.3:.01]
    x = x0 + x
    y = v0 + y
    s,t = x.shape
    c = 0*x
    s,t = x.shape
    for i in range(s):
       for j in range(t):
          xx,yy = x[i,j], y[i,j]
          c[i,j] = cost_function([xx, y0, z0, yy])

    sam.putarray('X',x)
    sam.putarray('Y',y)
    sam.putarray('C',c)

    sam.eval("[c,h]=contourf(X,Y,C,100);set(h,'EdgeColor','none')")
    sam.eval('hold on')

def matlab(scriptdir,title,heading,imagefile,htmlfile='results.html'):
    # use sam module to run a matlab session from python
    import sam
    sam.eval("") # initiate matlab session

    # scriptdir should contain the amroc matlab scripts needed to process
    # the hdf files.  make sure this is in MATLABPATH
    sam.eval("addpath %s" % scriptdir)

    # run plothdf script, using local menu_automation.m script to make plot
    sam.eval("plothdf")
    sam.close() # end matlab session

    # generate HTML web page with results
    makeHTML(title,heading,imagefile,htmlfile)
    return

def draw_contour():
    import numpy

    x, y = numpy.mgrid[0:7.5:0.05,0:7.5:0.05]
    c = 0*x
    s,t = x.shape
    for i in range(s):
       for j in range(t):
          xx,yy = x[i,j], y[i,j]
          c[i,j] = CostFunction([xx,yy])


    sam.putarray('X',x)
    sam.putarray('Y',y)
    sam.putarray('C',c)

    sam.verbose()    
    sam.eval("[c,h]=contourf(X,Y,log(C*20+1)+2,100);set(h,'EdgeColor','none')")
    sam.eval("title('Zimmermann''s Corner. Min at 7,2')")
    sam.eval('hold on')

def draw_contour():
    import numpy

    x, y = numpy.mgrid[0:2.1:0.05,0:2.1:0.05]
    c = 0*x
    s,t = x.shape
    for i in range(s):
       for j in range(t):
          xx,yy = x[i,j], y[i,j]
          c[i,j] = Corana2([xx,yy])


    sam.putarray('X',x)
    sam.putarray('Y',y)
    sam.putarray('C',c)

    sam.verbose()    
    sam.eval("[c,h]=contourf(X,Y,C,100);set(h,'EdgeColor','none')")
    #sam.eval("[c,h]=contourf(X,Y,log(C*20+1)+2,100);set(h,'EdgeColor','none')")
    sam.eval("title('Corana''s Parabola in 2D. Min at 0,0')")
    sam.eval('hold on')

def draw_contour_xy():
    import numpy
    x, y = mgrid[-40:40:0.5, -40:40:0.5]
    x = x0 + x
    y = y0 + y
    s,t = x.shape
    c = 0*x
    s,t = x.shape
    for i in range(s):
       for j in range(t):
          xx,yy = x[i,j], y[i,j]
          c[i,j] = cost_function([xx,yy, z0, v0])


    sam.putarray('X',x)
    sam.putarray('Y',y)
    sam.putarray('C',c)

    sam.verbose()    
    sam.eval("[c,h]=contourf(X,Y,C,100);set(h,'EdgeColor','none')")
    sam.eval("title('Mogi Fitting')")
    sam.eval('hold on')

# The dual problem verification begins here.
npt = xy.shape[0]
Q = dot(xy, transpose(xy))
#Q = zeros((npt,npt))
#for i in range(npt):
#   for j in range(npt):
##       Q[i,j] = dot(xy[i,:],xy[j,:])
#Q = array(Q)
f = -diag(Q)+10
H = Q*2
A = ones((1,npt))
b = ones(1)
sam.putarray('H',H);
sam.putarray('f',f);
sam.eval("npt = %d;" % npt);
sam.eval("al = quadprog(H,f,[],[],ones(1,npt),1,zeros(npt,1),ones(npt,1));")
alpha = sam.getarray('al').flatten()


def getobj(H,f, x):
    return 0.5 * dot(dot(x,H),x) + dot(f,x)

def chop(x):
    if abs(x) > 1e-6:
        return x
    else:
        return 0

x = array([chop(y) for y in alpha])
center = dot(alpha,xy)

    solver = fmin(len(xinit))
    solver.SetInitialPoints(xinit)
    solver.SetGenerationMonitor(simplex)
    solver.Solve(cost_function, termination=CRT())
    sol = solver.Solution()
    print(sol)

    for x in simplex.x:
        sam.putarray('x',x)
        sam.eval("plot(x([1,2,3,1],1),x([1,2,3,1],2),'w-','LineWidth',2)")
    return sol

draw_contour_xv()
xysol = run_once_xy()

sam.eval("figure(2)")

draw_contour_xy()
xvsol = run_once_xv()


sam.eval("figure(3)")
plot_noisy_data()
sam.eval("hold on")
plot_sol(xysol)
plot_sol(xvsol)

getch()

# end of file

    def renderOnSurface(self, surface):
        """Renders the scattering intensity on given detector surface in reciprocal space."""

        import sam
        import numpy as np

        ########################
        ###
        ### this is to be implemented using the 'slice()' method of Matlab.
        ###
        ########################

        #we first pass the data into Matlab as an array:
        array = np.array(self._getArray())
        (dim0, dim1, dim2) = array.shape
        q0 = np.arange(1, dim0 + 1)
        q1 = np.arange(1, dim1 + 1)
        q2 = np.arange(1, dim2 + 1)
        # set up the ranges for the grid in reciprocal space
        # these are some dummy ranges for now ([-0.5,0.5])
        q0 = (q0 - (dim0+1)/2.0) / dim0 
        q1 = (q1 - (dim1+1)/2.0) / dim1
        q2 = (q2 - (dim2+1)/2.0) / dim2

        sam.put('dim0', [dim0])
        sam.put('dim1', [dim1])
        sam.put('dim2', [dim2])

        sam.put('q0', q0)
        sam.put('q1', q1)
        sam.put('q2', q2)
        sam.eval("[qxi, qyi, qzi] = meshgrid(q0, q1, q2);")

        # this causes a segmentation fault for an array of 20x20x20 doubles:
        #sam.putarray('array', array)
        # here is a temporary ugly(!) fix:
        array.shape=(dim0*dim1*dim2,)
        sam.put('array', array)
        sam.eval("scattarray = reshape(array', dim0, dim1, dim2)")
        # end of fix

        # we need to pass the detector surface into Matlab:
        sam.eval("[sxi,syi] = meshgrid(q0,q1);")
        sam.eval("szi = sqrt(1.0-(sxi.^2))-0.75 ;")

        # We use an M-file that gets called when we evaluate the name
        # of the corresponding function in Matlab via SAM
        sam.eval("MlabScatteringIntensitySurfaceSlicer(qxi,qyi,qzi,scattarray, sxi, syi,szi)")

        waitforentry=raw_input("Press any key.")
        sam.eval("close")
        return 0