def scipyOptimize(recipe):
    """Optimize the recipe created above using scipy.
    """
    from scipy.optimize.minpack import leastsq
    print "Fit using scipy's LM optimizer"
    leastsq(recipe.residual, recipe.getValues())
    
    return

 def refine(self):
     '''Optimize the recipe created above using scipy.
     '''
     from scipy.optimize.minpack import leastsq
     leastsq(self.recipe.residual, self.recipe.values)
     self.results = FitResults(self.recipe)
     print("Fit results:\n")
     print(self.results)
     return

def doBestFit(compositeModel, params0, maxfev=None, factor=None):
    if not maxfev:
        maxfev = 500*(len(params0)+1)
    if not factor:
        factor = 100
    residual = compositeModel.residual
    if compositeModel.isAnalyticalDerivs:
        jacobian = compositeModel.jacobian
        full_output = leastsq(residual, params0,\
                            maxfev=maxfev, Dfun=jacobian, col_deriv=True, \
                            factor=factor, full_output=1)
    else:
        full_output = leastsq(residual, params0, maxfev=maxfev, \
                            factor=factor, full_output=1)
    return full_output

 def test_input_untouched(self):
     p0 = array([0, 0, 0], dtype=float64)
     p0_copy = array(p0, copy=True)
     full_output = leastsq(self.residuals, p0, args=(self.y_meas, self.x), full_output=True)
     params_fit, cov_x, infodict, mesg, ier = full_output
     assert_(ier in (1, 2, 3, 4), "solution not found: %s" % mesg)
     assert_array_equal(p0, p0_copy)

 def test_basic(self):
     p0 = array([0,0,0])
     params_fit, ier = leastsq(self.residuals, p0,
                               args=(self.y_meas, self.x))
     assert_(ier in (1,2,3,4), 'solution not found (ier=%d)'%ier)
     # low precision due to random
     assert_array_almost_equal(params_fit, self.abc, decimal=2)

    def minimize1(self, dim, theta, delta):
        # p0Cor = numpy.random.uniform(-1,1,dim**2).reshape(dim, dim)
        # p0Cor = p0Cor - numpy.diag(p0Cor) + numpy.identity(dim)
        
        #p0 = reverse(dim, numpy.identity(dim), numpy.ones(dim)/20)
        p0  = numpy.ones(dim)
        popt, _,infodict,mesg,_ = leastsq(errfunc, p0, args=(theta, delta),full_output=True)
        print(mesg)
         
         
        ss_err=(infodict['fvec']**2).sum()
        ss_tot=((delta-delta.mean())**2).sum()
        rsquared=1-(ss_err/ss_tot)
        print("rsquared", rsquared)
        
        # corrm, var, mu = transform(dim, popt)
        # var = var * self.gbest.position
        # _cov = corr2cov(corrm, var)
        print("found diag_inv_cov:\n", abs(popt))
        
        _cov = numpy.linalg.inv(numpy.diag(abs(popt)))

        #print("used mu:", mu)
        print("found _cov:\n", _cov)

        #sigma = numpy.sqrt(numpy.diag(_cov))
        #print( "=> found sigma:", sigma)
        
        return _cov

 def test_full_output(self):
     p0 = array([0,0,0])
     full_output = leastsq(self.residuals, p0,
                           args=(self.y_meas, self.x),
                           full_output=True)
     params_fit, cov_x, infodict, mesg, ier = full_output
     assert_(ier in (1,2,3,4), 'solution not found: %s'%mesg)

    def summarize(self, report, index):
        
        func1 = lambda v, x, mu: v[0] + v[1] * np.exp(-np.abs((np.array(x) - mu) / v[2]))
        full1 = lambda v, x: v[0] + v[1] * np.exp(-np.abs((np.array(x) - v[2]) / v[3]))

        func2 = lambda v, x, mu: v[0] + v[1] * np.exp(-((x - mu) / v[2]) ** 2) + v[3] * np.exp(-np.abs((np.array(x) - mu) / v[4]))
        full2 = lambda v, x: v[0] + v[1] * np.exp(-((x - v[2]) / v[3]) ** 2) + v[4] * np.exp(-np.abs((np.array(x) - v[2]) / v[5]))

        func3 = lambda v, x, mu: v[0] + v[1] * np.exp(-((np.array(x) - mu) / v[2]) ** 2)
        full3 = lambda v, x: v[0] + v[1] * np.exp(-np.abs((np.array(x) - v[2]) / v[3]) ** 2)
        
        # Force E to go to zero at minimum
        func4 = lambda v, x, mu: v[0] - v[0] * np.exp(-np.abs((np.array(x) - mu) / v[1]))
        full4 = lambda v, x: v[0] - v[0] * np.exp(-np.abs((np.array(x) - v[1]) / v[2]))
        
        funce = lambda v, x, y, mu: (func1(v, x, mu) - y)
        fulle = lambda v, x, y: (full1(v, x) - y)
        use = self.sensels
#        pylab.figure()
        for p in self.s:
            if p in use:
                f = report.figure(cols=2)
                e1_list = self.e1_lists[use.index(p)]
                for c, clist in enumerate(e1_list):
                    E1 = np.sum(np.array(clist), axis=0) / self.n_sampels[c]
                    V = np.std(np.array(clist), axis=0)
                    
                    v0 = np.array([1, -1, self.area / 2, 1])
                    mu = np.argmin(V)
                    v, _ = leastsq(funce, v0, args=(range(self.area), V, mu), maxfev=10000)
                    
                    v0 = [v[0], v[1], mu , v[3]]
                    v, _ = leastsq(fulle, v0, args=(range(self.area), V), maxfev=10000)
                    
                    x = np.linspace(0, self.area - 1, 100)
                    d = x[np.argmin(full1(v, x))] - self.area / 2
                    
#                    pylab.subplot(len(self.sensels), len(self.commands), use.index(p) * 2 + c + 1)
                    with f.plot('p%g_c%g' % (p, c), caption='Pixel: %g  Command: %g  Samples: %g' % (p, c, self.n_sampels[c])) as pylab:
                        pylab.hold(True)
                        pylab.errorbar(range(self.area), E1 * 0, V)
                        pylab.xlim((0, self.area))
                        
                        pylab.plot(x, full1(v, x), color='g')
                        
                    report.text('test%g_cmd_p%g_c%g' % (index, p, c), str(d))

def fitGaussianNLLS(xp,yp,patch):
	#print("patch: ",patch)
	paramsGuess = array([0,patch[hw,hw],hw,hw,1])
	tA = time.time()
	(paramsOut, cov, infodict, msg, ier) = minpack.leastsq(residualFunction,paramsGuess,(xp,yp,patch), full_output=1, ftol = .1, xtol = .1)
	#print("leastsq")
	#print("nfev = %d" % infodict['nfev'])
	#print(time.time()-tA)
	return paramsOut

def scipyOptimize(recipe):
    """Optimize the recipe created above using scipy.

    The FitRecipe we created in makeRecipe has a 'residual' method that we can
    be minimized using a scipy optimizer. The details are described in the
    source.

    """

    # We're going to use the least-squares (Levenberg-Marquardt) optimizer from
    # scipy. We simply have to give it the function to minimize
    # (recipe.residual) and the starting values of the Variables
    # (recipe.getValues()).
    from scipy.optimize.minpack import leastsq
    print "Fit using scipy's LM optimizer"
    leastsq(recipe.residual, recipe.getValues())

    return

    def test_reentrant_func(self):
        def func(*args):
            self.test_basic()
            return self.residuals(*args)

        p0 = array([0,0,0])
        params_fit, ier = leastsq(func, p0,
                                  args=(self.y_meas, self.x))
        assert_(ier in (1,2,3,4), 'solution not found (ier=%d)' % ier)
        # low precision due to random
        assert_array_almost_equal(params_fit, self.abc, decimal=2)

    def simulateSteadyState(self):
        """
        Perform a stady state simulation.

        Call self.clear() before using this method!

        This function computes one steady state solution, of the system of
        differential equations. Which solution of the potentially many solutions
        is found, depends on the initial guess. A steady state solution is a
        vector of all variables. This vector is appended to the array of results.
        Usually one will compute a series of stady state solutions, each with
        slightly different parameters.

        Initial guess: When there are no prior results, the initial values are
        (ab)used as an initial guess; otherwise the latest results are used as
        the initial guess.
        In the time array the count of current simulation is stored. This way the
        graph function still produces useful graphs with steady state simulations.

        The results can be displayed with the graph(...) function and stored
        with the store function. The funcion getAttributes(...) returns the
        simulation result of a speciffic attribute.
        """

        raise Exception('This method is currently broken!')

        #TODO: Use flag self._isSteadyStateStart instead of self.time
        if not hasattr(self, 'time'):
            #this is the first call in a row of steady state simulations - setup everything
            lastResult = -1
            self.resultArray = array([[0]], 'float64')
            self.time = array([0], 'float64')
            self.initialize()    #initial guess for root finder: initial values abused
            x0 = self.initialValues
            t0 = -1
        else:
            lastResult = shape(self.resultArray)[0]-1
            x0 = self.resultArray[lastResult, 0:self.stateVectorLen] #initial guess for root finder: last result
            t0 = self.time[lastResult]

        #compute the state variables of the steady state solution
        (xmin, msg) = minpack.leastsq(self.dynamic, x0, (t0)) #funcion will also report local minima that are no roots. Caution!
##        xmin = optimize.fsolve(self.dynamic, x0, (0)) #function is always stuck in one (the trivial) minimum
        #also compute the algebraic variables
        currRes = ones(4) #dummy to make pydev show no error
        #currRes = self.outputEquations(xmin)
        #expand the storage and save the results
        self.resultArray = resize(self.resultArray,
                                  (lastResult+2,
                                   self.stateVectorLen + self.algVectorLen))
        self.resultArray[lastResult+1,:] = currRes[0,:]
        self.time = resize(self.time, (lastResult+2,))
        self.time[lastResult+1] = t0 + 1

    def coherence(self, referenceStream, referenceAzimuth, NFFT=None, noverlap=None, Fs=None):
        """
        Determine angle of a north and east channel based on a reference north channel.

        @param referenceStream stream containing a reference North channel.
        @param referenceAzimuth reference angle of North channel.
        @param NFFT number of points to use in FFT.
        @param noverlap number of overlapping points between FFTs.
        @param Fs sampling frequency. If omitted, uses referenceStream sampling_rate.
        @return relative angle from unknown north to reference north
                (referenceAngle-angle = actual angle of unknown channel).
        """
        referenceExtent = referenceStream.getTimeExtent()
        unknownExtent = self.getTimeExtent()
        if referenceExtent['start'] != unknownExtent['start'] or referenceExtent['end'] != unknownExtent['end']:
            raise StreamsException("reference and unknown streams must have same time extent")
        if Fs is None:
            Fs = referenceStream.data.traces[0].stats.sampling_rate
        if noverlap is None:
            noverlap = NFFT / 4
        # internal function to determine the coherence of unrotated and rotated data
        def cohere1(theta):
            theta_r = math.radians(theta)
            rotated = (self.data[0].data)*math.cos(theta_r) + (self.data[1].data)*math.sin(theta_r)
            coh,fre = mlab.cohere(referenceStream.data[0].data, rotated,
                    NFFT=NFFT, noverlap=noverlap, Fs=Fs)
            return (coh - 1).sum()
        # most coherent angle of rotation
        theta1 = solv.leastsq(cohere1, 0)
        theta1 = normalize360(theta1[0][0])

        # rotate data and compare against reference stream
        rotated = self.rotate(theta1)
        rotatedData1 = rotated.data[0].data.astype('Float64')
        referenceData1 = referenceStream.data[0].data.astype('Float64')
        scale1 = sum(abs(rotatedData1)) / sum(abs(referenceData1))
        residual1 = sum(referenceData1**2-rotatedData1*scale1)**2
        rotatedData2 = rotated.data[1].data.astype('Float64')
        referenceData2 = referenceStream.data[1].data.astype('Float64')
        scale2 = sum(abs(rotatedData2)) / sum(abs(referenceData2))
        residual2 = sum(referenceData2**2-rotatedData2*scale2)**2

        return {
            'theta': theta1,
            'azimuth': referenceAzimuth - theta1,
            'rotated': rotated,
            'scale1': scale1,
            'residual1': residual1,
            'scale2': scale2,
            'residual2': residual2
        }

def least_example():
    print('least example')
    def f(x, b0, b1):
        return b0*sin(b1*x)


    def res(params, xdata, ydata, function):
        return function(xdata, *params) - ydata

    x = linspace(20,40, 100)
    y = f(x,2.4,1.14)
    args = (x,y,f)
    out = leastsq(res, [2.2, 1.2], args=args)
    print(out)

def aproximate_lineal(t_exp, y_exp, test=False):
    # define cost function - adapt to your usage
    #
    # single exponential
    function = single
    x0 = [0, y_exp[-1], t_exp.shape[0]]
    param = (t_exp, y_exp, function)

    # perform least squares fit
    A_final, cov_x, infodict, mesg, ier = leastsq(objective, x0, args=param, full_output=True)
    if ier != 1:
        print "No fit! %s " % mesg 
        return None
    y_final = function(A_final, t_exp)
    chi2 = sum((y_exp-y_final)**2 / y_final)
    if not test:
        return y_final
    else:
        return y_final,chi2

    def __init__(self, *args):
        apply(QWidget.__init__, (self,) + args)

	# make a QwtPlot widget
	self.plot = QwtPlot('A PyQwt and MinPack Demonstration', self)

	# initialize the noisy data
        scatter = 0.05
        x = arrayrange(-5.0, 5.0, 0.1)
        y = RandomArray.uniform(1.0-scatter, 1.0+scatter, shape(x)) * \
            function([1.0, 1.0, -2.0, 2.0], x)

        # fit from a reasonable initial guess
        guess = asarray([0.5, 1.5, -1.0, 3.0])
        yGuess = function(guess, x)
        solution = leastsq(function, guess, args=(x, y))
        yFit = function(solution[0], x)
        print solution

	# insert a few curves
	c1 = self.plot.insertCurve('data')
	c2 = self.plot.insertCurve('guess')
        c3 = self.plot.insertCurve('fit')
        
	# set curve styles
	self.plot.setCurvePen(c1, QPen(Qt.black))
        self.plot.setCurvePen(c2, QPen(Qt.red))
	self.plot.setCurvePen(c3, QPen(Qt.green))

	# copy the data
	self.plot.setCurveData(c1, x, y)
        self.plot.setCurveData(c2, x, yGuess)
        self.plot.setCurveData(c3, x, yFit)
	
	# set axis titles
	self.plot.setAxisTitle(QwtPlot.xBottom, 'x -->')
	self.plot.setAxisTitle(QwtPlot.yLeft, 'y -->')

        self.plot.enableLegend(1)
        self.plot.replot()

    def minimize1(self, dim, theta, delta):

        p0 = reverse(dim, numpy.identity(dim), numpy.ones(dim) / 20)
        popt, _, infodict, mesg, _ = leastsq(errfunc, p0, args=(theta, delta, self.gbest.position), full_output=True)
        print(mesg)

        ss_err = (infodict["fvec"] ** 2).sum()
        ss_tot = ((delta - delta.mean()) ** 2).sum()
        rsquared = 1 - (ss_err / ss_tot)
        print("rsquared", rsquared)

        corrm, var, mu = transform(dim, popt)
        var = var * self.gbest.position
        _cov = corr2cov(corrm, var)

        print("used mu:", mu)
        print("found _cov:\n", _cov)

        sigma = numpy.sqrt(numpy.diag(_cov))
        print("=> found sigma:", sigma)

        return _cov

def fit(x_array, y_array, function, A_start):
    """
    used to fit things
    20101209/RB: started
    
    INPUT:
    x_array: the array with time or something
    y-array: the array with the values that have to be fitted
    function: one of the functions, in the format as in the file "functions"
    A_start: a starting point for the fitting
    
    OUTPUT:
    A_final: the final parameters of the fitting
    
    WARNING:
    Always check the result, it might sometimes be sensitive to a good starting point.

    """
    param = (x_array, y_array, function)

    A_final, cov_x, infodict, mesg, ier = leastsq(minimize, A_start, args=param, full_output=True)  # , warning=True)

    return A_final

    int0 = numpy.argsort(tau)
    tau = tau[int0]; delflx2 = delflx2[int0]; delflxerr = delflxerr[int0]


    # Fitting function; exponential model; sf = sf_inf*(1-e^(-t/Tau))^0.5 
    fp = lambda v, x: v[0]*numpy.power((1.-numpy.exp(-numpy.abs(x)/v[1])),0.5)

    # Error function
    e  = lambda v, x, y, dy: (fp(v,x)-y)/dy  # error function for the leastsq 

    # Initial guess
    v0 = [0.2, 100.]    # initial guess for exp

    # leastsq
    v_whole, cov, info, mesg, ier = sc.leastsq(e, v0, args=(tau,numpy.sqrt(delflx2),delflxerr),full_output=True) 
    if ier != 1: print 'Fail in fitting SF for ' + str(inFiles[i])
    if cov == None:
        t0 = 0.0; t1 = 0.0
    else:
        t0 = numpy.sqrt(cov[0][0]); t1 = numpy.sqrt(cov[1][1])

    expected = v_whole[0]*numpy.power((1.-numpy.exp(-numpy.abs(tau)/v_whole[1])),0.5) # Expected Value
    chi2 = numpy.sum((numpy.sqrt(delflx2)-expected)**2/expected)

    WriteColumns(outDir+'tau%s.dat' % (file),((reshape(tau,(len(tau),1))),), header='# SF')
    WriteColumns(outDir+'sf%s.dat' % (file),((reshape(numpy.sqrt(delflx2),(len(tau),1))),), header='# SF')  
    outf.write('%.5f %.5f %.5f %.5f %.5f %i\n' % (v_whole[0],t0,v_whole[1],t1,chi2,len(tau))) 

outf.close()

# main program
random.seed(12345)       # seed or the random number generator

# loads data from file
data = loadtxt(os.path.join(datadir,'count.txt'), skiprows = 1)    
yvec = data[:, 0]
xmat = data[:, 1:data.shape[1]]
xmat = hstack([ones((data.shape[0], 1)), xmat])

data ={'yvec':yvec, 'xmat':xmat}

# use bayesian regression to initialise
bayesreg = BayesRegression(yvec, xmat)     
sig, beta0 = bayesreg.posterior_mean()

init_beta, info = leastsq(minfunc, beta0, args = (yvec, xmat))
data['betaprec'] =-llhessian(data, init_beta)
scale = linalg.inv(data['betaprec'])

# Initialise the random walk MH algorithm
samplebeta = RWMH(posterior, scale, init_beta, 'beta')

ms = MCMC(20000, 4000, data, [samplebeta],
          loglike = (logl, xmat.shape[1], 'yvec'))
ms.sampler()

ms.output()
ms.plot('beta')