code: from __future__ import print_function
  
  """
  Wrapper for the C++ drizzle overlap code.
  
  History
  -------
  
  Created by Jon Nielsen in 2012
  Updated for the new cubing algorithm by Francesco D'Eugenio 16/02/2017
  
  Notes
  -----
  This module contains a testing function. At the moment it requires that the libraries path be hardcoded (FDE).
  """
  
  import ctypes as C
  import numpy as np
  import os.path
  
  # Load the shared library
  try:
      libcm = C.CDLL(os.path.join(os.path.dirname(__file__), "libcCirc.so"))
  except:
      raise ImportError
      pass
      #libcm = C.CDLL(os.path.join('/home/franz/software/dev/sami-software-dev/dr0.10/utils', 'libcCirc.so'))
  
  # Specify the arguments our function takes:
  #  First argument is a 1D array of doubles.  It must be contiguous in memory.
  #  Second argument is a regular C long.
  #  Third argument is a pointer to a double.
  libcm.weight_map.argtypes = [
      C.c_long, C.c_long, C.c_double, C.c_double, C.c_double,
      np.ctypeslib.ndpointer(dtype='d', ndim=1, flags='C_CONTIGUOUS')]
  
  
  
  def resample_circle(nx, ny, xc, yc, r, *args):
    output = np.zeros(nx*ny)
    libcm.weight_map(nx, ny, xc, yc, r, output)
    return output.reshape(ny, nx)
  
  
  
  
  # Specify the arguments our function takes:
  #  First argument is a 1D array of doubles.  It must be contiguous in memory.
  #  Second argument is a regular C long.
  #  Third argument is a pointer to a double.
  libcm.weight_map_Gaussian.argtypes = [
       C.c_long, C.c_long, C.c_double, C.c_double, C.c_double, C.c_double,
       np.ctypeslib.ndpointer(dtype='d', ndim=1, flags='C_CONTIGUOUS')]
  
  
  
  def inteGrauss2d(nx, ny, xc, yc, sigma, n_sigma):
    output = np.zeros(nx*ny)
    libcm.weight_map_Gaussian(nx, ny, xc, yc, sigma, n_sigma, output)
    return output.reshape(ny, nx)
  
  
  
  if __name__ == "__main__":
  
      import time
      import pylab as plt
      from mpl_toolkits.axes_grid1 import ImageGrid
      plt.switch_backend('pdf')
  
      try:
          from texttable import Texttable
  
          def print_table(headers, times):
              table = Texttable()
              headers = ['Method'] + headers
              times = [u'Elapsed time (\u03bcs)'] + times
              table.add_rows([[h, t] for h, t in zip(headers, times)])
              table.set_cols_align(["l", "c"])
  
              print(table.draw())
  
      except ImportError:
          warnings.warn('To print formatted output, please install `texttable`')
  
          def print_table(headers, times):
              for h, e in zip(headers, elapsed):
                  print(h, '{:8.2f}'.format(e/float(n_iterations)*1.e6),
                        u'\u03bc'+'s')
  
  
      from circ import resample_circle as rc_py
      from circ import inteGrauss2d as ig_py
      rc_C = resample_circle
      ig_C = inteGrauss2d
  
      # This circle show that the python implementation of cCirc wraps around,
      # which is dumb
      # (50, 50, 27.527636662069785, 2.716882503265406, 4.9423403572267581,
      #  2.0454945513296701))
  
      # Define some representative cases.
      # Each has (xpix, ypix, xc, yc, [r/sigma], [clip])
      tests = ((10, 10, 5, 5, 1, 3),
               (10, 10, 3, 5, 1, 3),
               (10, 10, 3.1, 5.7, 1, 3))
      tests = ((50, 50, 27.527636662069785, 23.716882503265406, 4.9423403572267581, 2.0454945513296701),
               (50, 50, 27.527636662069785, 2.716882503265406, 4.9423403572267581, 2.0454945513296701)) # <=
  
      # Define some circles to plot. Notice that for `resample_circle` the last
      # element of each tuple is ignored, while the penultimate is the radius.
      # for `inteGrauss2d` the last element defines the truncation and the
      # penultimate element defines the standard deviation of the Gaussian.
      tests = (
          (10, 10, 5., 1.0, 2., 2.),
          (10, 10, 5.65, 1.2, 2.31, 2.001))
      tests = ((50, 50, 9., 9., 1.6, 5.),)
  
      # Number of iterations for the benchmark.
      n_iterations = 10000
  
      n_tests = len(tests)
  
      plt.clf()
      fig = plt.figure(figsize=(18, 3 * n_tests))
      grid = ImageGrid(fig, 111, nrows_ncols=(n_tests, 6),
                       cbar_location='top', cbar_mode='edge',
                       cbar_pad='0%', cbar_size='5%', direction='row',
                       axes_pad=(0, 0))
   
      wmaps = [[rc_py(*t), rc_C(*t), ig_py(*t), ig_C(*t)] for t in tests]
      wmaps = [[w[0], w[1], w[0]-w[1], w[2], w[3], w[2]-w[3]] for w in wmaps]
      wmaps = sum(wmaps, []) # Flatten the list.
   
      # Reformat `tests` so that we can iterate over it for each subplot.
      tests = ((t,)*6 for t in tests) # Repeat each test for every plot that uses it.
      tests = sum(tests, ())          # Flatten this tuple.
      
  
      y_labels = ['$\mathrm{Test\;'+ '{}'.format(n) + '}$'
                  for n in range(1, n_tests+1)]
      x_labels = ['$\mathrm{resample\_circ \; (python)}$',
                  '$\mathrm{resample\_circ \; (C++)}$',
                  '$\mathrm{residuals}$',
                  '$\mathrm{inteGrauss2d \; (python)}$',
                  '$\mathrm{inteGrauss2d \; (C++)}$',
                  '$\mathrm{residuals}$']
  
      for n, (wm, t, ax) in enumerate(zip(wmaps, tests, grid)):
          img = ax.imshow(wm, origin='lower', interpolation='none')
          ax.plot([t[2]-.5], [t[3]-.5], 'mx')
          if n >= (n_tests - 1) * 6: # Trigger the colourbar.
              cbar = ax.cax.colorbar(img)
              cbar.solids.set_edgecolor('face')
  
          if n % 6 == 0:              # We are at the first column.
              ax.set_ylabel(y_labels[n // 6])
          if n // 6 >= (n_tests - 1): # We are at the bottom row.
              ax.set_xlabel(x_labels[n % 6])
   
      #plt.subplots_adjust(wspace=0., hspace=0.)
      plt.savefig('test_ccircs.pdf')
   
   
   
      # Benchmarking.
      
      # Create random tests to prevent system caching of interemdiate results in
      # the C implementation - that would be unfair.
   
      grid_size = np.repeat(50, n_iterations)
      xc, yc = np.random.uniform(10, 40, (2, n_iterations))
      radius = np.random.uniform(1, 5, n_iterations)
      n_sigma = np.random.uniform(1, 5, n_iterations)
   
      tests = [
          (xpix, ypix, x, y, r, ns)
          for xpix, ypix, x, y, r, ns in zip(
              grid_size, grid_size, xc, yc, radius, n_sigma)]
   
      start_time, elapsed, results = [], [], []
   
      print('Benchmarking function `resample_circle` (Top hat cubing, python)...')
      start_time.append(time.time())
      results.append(np.array([rc_py(*t) for t in tests]))
      elapsed.append(time.time() - start_time[-1])
   
      print('Benchmarking function `resample_circle` (Top hat cubing, C++)...')
      start_time.append(time.time())
      results.append(np.array([rc_C(*t) for t in tests]))
      elapsed.append(time.time() - start_time[-1])
   
      print('Benchmarking function `inteGrauss2d` (Gaussian cubing, python)...')
      start_time.append(time.time())
      results.append(np.array([ig_py(*t) for t in tests]))
      elapsed.append(time.time() - start_time[-1])
   
      print('Benchmarking function `inteGrauss2d` (Gaussian cubing, C++)...')
      start_time.append(time.time())
      results.append(np.array([ig_C(*t) for t in tests]))
      elapsed.append(time.time() - start_time[-1])
   
      print('Summary:')
      headers = ['`resample_circle` (Top hat cubing, python)',
                 '`resample_circle` (Top hat cubing, C++)',
                 '`inteGrauss2d` (Gaussian cubing, python)',
                 '`inteGrauss2d` (Gaussian cubing, C++)']
      print_table(headers, elapsed)
   """
  This file contains various fitting functions for general use with SAMI codes.
  
  Currently included:
  
  GaussFitter - Gaussian Fitter (1d)
  GaussHermiteFitter - Fits a truncated Gauss-Hermite expansion (1d)
  TwoDGaussFitter - Gaussian Fitter (2d, optionally w/ PA and different widths)
  
  Would be nice:
  
  Exponential Fitter?
  Others?
  
  Example of the class format is below. Should have a list of things that can
  be accessed in all class definitions (chi-2 etc.)
  
  To use these fitting classes, initialise them using the initial guesses of the
  parameters, along with the coordinates, data and (optionally) weights. Then
  call the fit function to perform the fit. The best fit parameters are then
  stored in p. For example:
  
  my_fitter = TwoDGaussFitter(initial_p, x, y, data, weights)
  my_fitter.fit()
  best_fit_p = my_fitter.p
  
  If you want to integrate over each fibre, use the fibre_integrator function
  *before* performing the fit. The diameter must be provided in the same units
  that x and y will be in. For example:
  
  my_fitter = TwoDGaussFitter(initial_p, x, y, data, weights)
  fibre_integrator(my_fitter, 1.6)
  my_fitter.fit()
  best_fit_p = my_fitter.p
  
  For integration over square pixels rather than round fibres, use:
  fibre_integrator(my_fitter, 0.7, pixel=True)
  
  Calling an instance of a fitter will return the model values at the provided
  coordinates. So, after either of the above examples:
  
  my_fitter(x, y)
  
  would return the best-fit model values at the coordinates (x, y).
  
  TODO: Make a BaseFitter class containing the basic functionality that other
  classes can inherit from.
  
  TODO: Implement limits in a better way, i.e. use a minimisation function
  that incorporates limits rather than the "return 1e99" method used below.
  
  TODO: Rename this module to fitting, rather than samifitting.
  """
  from __future__ import absolute_import, division, print_function, unicode_literals
  
  from scipy.optimize import leastsq
  import scipy as sp
  import numpy as np
  
  class FittingException(Exception):
      """Could I make this do something useful?"""
      pass
  
  class GaussFitter:
      """ Fits a 1d Gaussian to data. Params in form list p (amplitude, mean, sigma, offset). Offset is optional."""
  
      def __init__(self, p, x, y, weights=None):
          self.p_start = p
          self.p = p
          self.x = x
          self.y = y
          if weights == None:
              self.weights = sp.ones(len(self.y))
          else:
              self.weights = weights
  
          self.perr = 0.
          self.var_fit = 0.
  
          if len(p) == 4 and p[0]>0.:
              self.fitfunc = self.f1
          elif len(p) == 3:
              # no base
              self.fitfunc = self.f2
          elif len(p) == 4 and p[0] < 0.:
              self.p[0] = abs(self.p[0])
              self.fitfunc = self.f3
      
      def f1(self, p, x): 
          return p[0]*sp.exp(-(p[1]-x)**2/(2*p[2]**2)) + p[3]
      
      def f2(self, p, x): 
          return p[0]*sp.exp(-(p[1]-x)**2/(2*p[2]**2)) + 0.
      
      def f3(self, p, x): 
          return -p[0]*sp.exp(-(p[1]-x)**2/(2*p[2]**2)) + p[3]
  
      def errfunc(self, p, x, y, weights):
          # if width < 0 return input
          if p[2] < 0. or p[0] < 0.:
              # if we get a negative sigma value then penalise massively because
              # that is silly.
              # This method isn't great, as it can make leastsq confused about
              # the convergence criteria. Would be better to use something like
              # scipy.optimize.minimize that properly incorporates limits. Same
              # goes for the other classes in this module.
              return 1e99
          else:
              return weights*(self.fitfunc(p, x) - y)
  
      def fit(self):
  
          self.p, self.cov_x, self.infodict, self.mesg, self.success = \
          leastsq(self.errfunc, self.p, \
                  args=(self.x, self.y, self.weights), full_output=1)
  
          var_fit = (self.errfunc(self.p, self.x, \
              self.y, self.weights)**2).sum()/(len(self.y)-len(self.p))
  
          self.var_fit = var_fit
  
          if self.cov_x is not None:
              self.perr = sp.sqrt(self.cov_x.diagonal())*self.var_fit
  
          if not self.success in [1,2,3,4]:
              print("Fit Failed")
              #raise FittingException("Fit failed") # This does nothing.
              
          self.linestr=self.p[0]*self.p[2]*sp.sqrt(2*sp.pi)
          #self.line_err=S.sqrt(self.linestr*self.linestr*((self.perr[0]/self.p[0])**2+(self.perr[2]/self.p[2])**2))
  
      def __call__(self, x):
          return self.fitfunc(self.p, x)
  
  class GaussHermiteFitter:
      """Parameters list p contains, in order, amplitude, mean, sigma, h3, h4, bias, where the bias is optional"""
  
      def __init__(self, p, x, y, weights=None):
          self.p_start = p
          self.p = p
          self.x = x
          self.y = y
          if weights == None:
              self.weights = S.ones(len(self.y))
          else:
              self.weights = weights
  
          self.perr = 0.
          self.var_fit = 0.
  
          # Which function to use depends on the input parameters 
          if len(p) == 5 and p[0]>0.:
              # Fit a truncated Gauss-Hermite sequence.
              self.fitfunc = self.f1
          elif len(p) == 6 and p[0]>0.:
              # Fit a truncated Gauss-Hermite sequence with a bias.
              self.fitfunc = self.f2
          else:
              raise Exception
      
      def f1(self, p, x):
          w=(p[1]-x)/(p[2])
          H3=(p[3]*sp.sqrt(2)/sp.sqrt(6))*((2*w**3)-(3*w))
          H4=(p[4]/sp.sqrt(24))*(4*w**4-12*w**2+3)
          gauss=p[0]*sp.exp(-w**2/2)
          
          gh=gauss*(1+H3+H4)
          return gh
  
      def f2(self, p, x):
          w=(p[1]-x)/(p[2])
          H3=(p[3]*sp.sqrt(2)/sp.sqrt(6))*((2*w**3)-(3*w))
          H4=(p[4]/sp.sqrt(24))*(4*w**4-12*w**2+3)
          gauss=p[0]*sp.exp(-w**2/2)
          
          gh2=gauss*(1+H3+H4)+p[5]
          return gh2
  
      def errfunc(self, p, x, y, weights):
          if p[2] < 0. or p[0] < 0.:
              # if we get a negative sigma value then penalise massively because
              # that is silly.
              return 1e99
          else:
              return weights*(self.fitfunc(p, x) - y)
  
      def fit(self):
  
          self.p, self.cov_x, self.infodict, self.mesg, self.success = \
          leastsq(self.errfunc, self.p, \
                  args=(self.x, self.y, self.weights), full_output=1)
  
          var_fit = (self.errfunc(self.p, self.x, \
              self.y, self.weights)**2).sum()/(len(self.y)-len(self.p))
  
          self.var_fit = var_fit
  
          if self.cov_x is not None:
              self.perr = S.sqrt(self.cov_x.diagonal())*self.var_fit
  
          # Would like to return the linestrength and associated error
          gamma=self.p[0]*self.p[2]*S.sqrt(2*S.pi)
          gamma_err=S.sqrt(gamma*gamma*((self.perr[0]/self.p[0])**2+(self.perr[1]/self.p[1])**2))
          self.linestr=gamma*(1+S.sqrt(6)*self.p[4]/4)
          self.line_err=S.sqrt(gamma_err**2*(1+S.sqrt(6)*self.p[4]/4)**2+self.perr[4]**2*(S.sqrt(6)*gamma_err/4)**2)
  
          if not self.success in [1,2,3,4]:
              print("Fit Failed...")
              #raise FittingException("Fit failed")
  
      def __call__(self, x):
          return self.fitfunc(self.p, x)
  
  class TwoDGaussFitter:
      """ Fits a 2d Gaussian with PA and ellipticity. Params in form (amplitude, mean_x, mean_y, sigma_x, sigma_y,
      rotation, offset). Offset is optional. To fit a circular Gaussian use (amplitude, mean_x, mean_y, sigma, offset),
      again offset is optional."""
  
      def __init__(self, p, x, y, z, weights=None):
          self.p_start = p
          self.p = p
          self.x = x
          self.y = y
          self.z = z
          
          if weights == None:
              self.weights = sp.ones(len(self.z))
          else:
              self.weights = weights
  
          self.perr = 0.
          self.var_fit = 0.
  
          if len(p) == 7:
              # 2d elliptical Gaussian with offset.
              self.p[0] = abs(self.p[0]) # amplitude should be positive.
              self.fitfunc = self.f1
              
          elif len(p) == 6:
              # 2d elliptical Gaussian witout offset.
              self.p[0] = abs(self.p[0]) # amplitude should be positive.
              self.fitfunc = self.f2
              
          elif len(p) == 5:
              # 2d circular Gaussian with offset.
              self.p[0] = abs(self.p[0]) # amplitude should be positive.
              self.fitfunc = self.f3
              
          elif len(p) == 4:
              self.p[0] = abs(self.p[0])
              self.fitfunc = self.f4
  
          else:
              raise Exception
  
      def f1(self, p, x, y):
          # f1 is an elliptical Gaussian with PA and a bias level.
  
          rot_rad=p[5]*sp.pi/180 # convert rotation into radians.
  
          rc_x=p[1]*sp.cos(rot_rad)-p[2]*sp.sin(rot_rad)
          rc_y=p[1]*sp.sin(rot_rad)+p[2]*sp.cos(rot_rad)
      
          return p[0]*sp.exp(-(((rc_x-(x*sp.cos(rot_rad)-y*sp.sin(rot_rad)))/p[3])**2\
                                      +((rc_y-(x*sp.sin(rot_rad)+y*sp.cos(rot_rad)))/p[4])**2)/2)+p[6]
  
      def f2(self, p, x, y):
          # f2 is an elliptical Gaussian with PA and no bias level.
  
          rot_rad=p[5]*sp.pi/180 # convert rotation into radians.
  
          rc_x=p[1]*sp.cos(rot_rad)-p[2]*sp.sin(rot_rad)
          rc_y=p[1]*sp.sin(rot_rad)+p[2]*sp.cos(rot_rad)
      
          return p[0]*sp.exp(-(((rc_x-(x*sp.cos(rot_rad)-y*sp.sin(rot_rad)))/p[3])**2\
                                      +((rc_y-(x*sp.sin(rot_rad)+y*sp.cos(rot_rad)))/p[4])**2)/2)
  
      def f3(self, p, x, y):
          # f3 is a circular Gaussian, p in form (amplitude, mean_x, mean_y, sigma, offset).
          return p[0]*sp.exp(-(((p[1]-x)/p[3])**2+((p[2]-y)/p[3])**2)/2)+p[4]
  
      def f4(self, p, x, y):
          # f4 is a circular Gaussian as f3 but without an offset
          return p[0]*sp.exp(-(((p[1]-x)/p[3])**2+((p[2]-y)/p[3])**2)/2)
  
      def errfunc(self, p, x, y, z, weights):
          # if width < 0 return input
          #if p[3] < 0. or p[4] < 0. or p[0] < 0.: # This isn't valid for the circular case!
              # if we get negative sigma values then penalise massively because
              # that is silly.
              #return 1e99
          #else:
          return weights*(self.fitfunc(p, x, y) - z)
  
      def fit(self):
  
          #print(np.shape(self.x), np.shape(self.y), np.shape(self.z))
  
          self.p, self.cov_x, self.infodict, self.mesg, self.success = \
          leastsq(self.errfunc, self.p, \
                  args=(self.x, self.y, self.z, self.weights), full_output=1)
  
          var_fit = (self.errfunc(self.p, self.x, \
              self.y, self.z, self.weights)**2).sum()/(len(self.z)-len(self.p))
  
          self.var_fit = var_fit
  
          if self.cov_x is not None:
              self.perr = sp.sqrt(self.cov_x.diagonal())*self.var_fit
  
          if not self.success in [1,2,3,4]:
              print("Fit Failed...")
              #raise ExpFittingException("Fit failed")
  
      def __call__(self, x, y):
          return self.fitfunc(self.p, x, y)
      
  def fibre_integrator(fitter, diameter, pixel=False):
      """Edits a fitter's fitfunc so that it integrates over each SAMI fibre."""
  
      # Save the diameter; not used here but could be useful later
      fitter.diameter = diameter
  
      # Define the subsampling points to use
      n_pix = 5       # Number of sampling points across the fibre
      # First make a 1d array of subsample points
      delta_x = np.linspace(-0.5 * (diameter * (1 - 1.0/n_pix)),
                            0.5 * (diameter * (1 - 1.0/n_pix)),
                            num=n_pix)
      delta_y = delta_x
      # Then turn that into a 2d grid of (delta_x, delta_y) centred on (0, 0)
      delta_x = np.ravel(np.outer(delta_x, np.ones(n_pix)))
      delta_y = np.ravel(np.outer(np.ones(n_pix), delta_y))
      if pixel:
          # Square pixels; keep everything
          n_keep = n_pix**2
      else:
          # Round fibres; only keep the points within one radius
          keep = np.where(delta_x**2 + delta_y**2 < (0.5 * diameter)**2)[0]
          n_keep = np.size(keep)
          delta_x = delta_x[keep]
          delta_y = delta_y[keep]
  
      old_fitfunc = fitter.fitfunc
  
      def integrated_fitfunc(p, x, y):
          # The fitter's fitfunc will be replaced by this one
          n_fib = np.size(x)
          x_sub = (np.outer(delta_x, np.ones(n_fib)) +
                   np.outer(np.ones(n_keep), x))
          y_sub = (np.outer(delta_y, np.ones(n_fib)) +
                   np.outer(np.ones(n_keep), y))
          return np.mean(old_fitfunc(p, x_sub, y_sub), 0)
  
      # Replace the fitter's fitfunc
      fitter.fitfunc = integrated_fitfunc
  
      return
  
  
apis: numpy.random.uniform numpy.ctypeslib.ndpointer numpy.repeat numpy.zeros
apis: scipy.sqrt numpy.ones scipy.cos scipy.exp scipy.optimize.leastsq numpy.size scipy.sin numpy.where numpy.linspace