18 from __future__
import division, print_function
23 import _SourceXtractorPy
as cpp
24 from .measurement_images
import MeasurementGroup
26 from astropy
import units
as u
27 from astropy.coordinates
import SkyCoord
28 from astropy.coordinates
import Angle
40 Limit, and normalize, the range of values for a model fitting parameter.
45 limits : a tuple (min, max), or a callable that receives a source, and returns a tuple (min, max)
51 Normalized to engine space using a sigmoid function
55 engine = \ln \frac{world - min}{max-world} \\
56 world = min + \frac{max - min}{1 + e^{engine}}
59 Normalized to engine space using an exponential sigmoid function
63 engine = \ln \left( \frac{\ln(world/min)}{\ln(max /world)} \right) \\
64 world = min * e^\frac{ \ln(max / min) }{ (1 + e^{-engine}) }
80 Receives a source, and returns a tuple (min, max)
97 Human readable representation for the object
100 if hasattr(self.
__limits,
'__call__'):
105 RangeType.LINEAR :
'LIN',
106 RangeType.EXPONENTIAL :
'EXP'
108 res +=
',{}]'.format(type_str[self.
__type])
114 Unbounded, but normalize, value of a model fitting parameter
118 normalization_factor: float, or a callable that receives the initial value parameter value and a source,
120 The world value which will be normalized to 1 in engine coordinates
134 Receives the initial parameter value and a source, and returns the world value which will be
135 normalized to 1 in engine coordinates
144 Human readable representation for the object
155 constant_parameter_dict = {}
156 free_parameter_dict = {}
157 dependent_parameter_dict = {}
162 Print a human-readable representation of the configured model fitting parameters.
167 Where to print the representation. Defaults to sys.stderr
169 if constant_parameter_dict:
170 print(
'Constant parameters:', file=file)
171 for n
in constant_parameter_dict:
172 print(
' {}: {}'.format(n, constant_parameter_dict[n]), file=file)
173 if free_parameter_dict:
174 print(
'Free parameters:', file=file)
175 for n
in free_parameter_dict:
176 print(
' {}: {}'.format(n, free_parameter_dict[n]), file=file)
177 if dependent_parameter_dict:
178 print(
'Dependent parameters:', file=file)
179 for n
in dependent_parameter_dict:
180 print(
' {}: {}'.format(n, dependent_parameter_dict[n]), file=file)
185 Base class for all model fitting parameter types.
186 Can not be used directly.
194 Human readable representation for the object
196 return '(ID:{})'.format(self.id)
201 A parameter with a single value that remains constant. It will not be fitted.
205 value : float, or callable that receives a source and returns a float
206 Value for this parameter
213 ParameterBase.__init__(self)
215 constant_parameter_dict[self.id] = self
222 Receives a source and returns a value for the parameter
231 Human readable representation for the object
233 res = ParameterBase.__str__(self)[:-1] +
', value:'
234 if hasattr(self.
__value,
'__call__'):
243 A parameter that will be fitted by the model fitting engine.
247 init_value : float or callable that receives a source, and returns a float
248 Initial value for the parameter.
249 range : instance of Range or Unbounded
250 Defines if this parameter is unbounded or bounded, and how.
259 >>> sersic = FreeParameter(2.0, Range((1.0, 7.0), RangeType.LINEAR))
266 ParameterBase.__init__(self)
269 free_parameter_dict[self.id] = self
276 Receives a source, and returns an initial value for the parameter.
293 Human readable representation for the object
295 res = ParameterBase.__str__(self)[:-1] +
', init:'
301 res +=
', range:' + str(self.
__range)
307 A DependentParameter is not fitted by itself, but its value is derived from another Parameters, whatever their type:
308 FreeParameter, ConstantParameter, or other DependentParameter
313 A callable that will be called with all the parameters specified in this constructor each time a new
314 evaluation is needed.
315 params : list of ParameterBase
316 List of parameters on which this DependentParameter depends.
320 >>> flux = get_flux_parameter()
321 >>> mag = DependentParameter(lambda f: -2.5 * np.log10(f) + args.mag_zeropoint, flux)
322 >>> add_output_column('mf_mag_' + band, mag)
329 ParameterBase.__init__(self)
332 dependent_parameter_dict[self.id] = self
337 Convenience function for the position parameter X and Y.
342 X coordinate, starting at the X coordinate of the centroid and linearly limited to X +/- the object radius.
344 Y coordinate, starting at the Y coordinate of the centroid and linearly limited to Y +/- the object radius.
347 X and Y are fitted on the detection image X and Y coordinates. Internally, these are translated to measurement
348 images using the WCS headers.
351 FreeParameter(
lambda o: o.get_centroid_x(),
Range(
lambda v,o: (v-o.get_radius(), v+o.get_radius()), RangeType.LINEAR)),
352 FreeParameter(
lambda o: o.get_centroid_y(),
Range(
lambda v,o: (v-o.get_radius(), v+o.get_radius()), RangeType.LINEAR))
358 Possible flux types to use as initial value for the flux parameter.
359 Right now, only isophotal is supported.
366 Convenience function for the flux parameter.
371 One of the values defined in FluxParameterType
373 Scaling of the initial flux. Defaults to 1.
378 Flux parameter, starting at the flux defined by `type`, and limited to +/- 1e3 times the initial value.
381 FluxParameterType.ISO :
'get_iso_flux'
383 return FreeParameter(
lambda o: getattr(o, func_map[type])() * scale,
Range(
lambda v,o: (v * 1E-3, v * 1E3), RangeType.EXPONENTIAL))
391 Model a Gaussian prior on a given parameter.
395 param : ParameterBase
396 Model fitting parameter
397 value : float or callable that receives a source and returns a float
399 sigma : float or callable that receives a source and returns a float
400 Standard deviation of the Gaussian
407 cpp.Id.__init__(self)
409 self.
value = value
if hasattr(value,
'__call__')
else lambda o: value
410 self.
sigma = sigma
if hasattr(sigma,
'__call__')
else lambda o: sigma
415 Add a prior to the given parameter.
419 param : ParameterBase
420 value : float or callable that receives a source and returns a float
422 sigma : float or callable that receives a source and returns a float
423 Standard deviation of the Gaussian
425 prior =
Prior(param, value, sigma)
426 prior_dict[prior.id] = prior
429 frame_models_dict = {}
434 if isinstance(x, tuple):
437 if not x.id
in frame_models_dict:
438 frame_models_dict[x.id] = []
439 frame_models_dict[x.id].append(model.id)
444 Add a model to be fitted to the given group.
448 group : MeasurementGroup
451 if not isinstance(group, MeasurementGroup):
452 raise TypeError(
'Models can only be added on MeasurementGroup, got {}'.format(type(group)))
453 if not hasattr(group,
'models'):
455 group.models.append(model)
461 Print a human-readable representation of the configured models.
465 group : MeasurementGroup
466 Print the models for this group.
468 If True, print also the parameters that belong to the model
470 Prefix each line with this string. Used internally for indentation.
472 Where to print the representation. Defaults to sys.stderr
474 if hasattr(group,
'models')
and group.models:
475 print(
'{}Models:'.format(prefix), file=file)
476 for m
in group.models:
477 print(
'{} {}'.format(prefix, m.to_string(show_params)), file=file)
479 if isinstance(x, tuple):
480 print(
'{}{}:'.format(prefix, x[0]), file=file)
484 constant_model_dict = {}
485 point_source_model_dict = {}
486 sersic_model_dict = {}
487 exponential_model_dict = {}
488 de_vaucouleurs_model_dict = {}
489 params_dict = {
"max_iterations": 100,
"modified_chi_squared_scale": 10,
"engine":
""}
497 Max number of iterations for the model fitting.
499 params_dict[
"max_iterations"] = iterations
507 Sets u0, as used by the modified chi squared residual comparator, a function that reduces the effect of large
509 Refer to the SourceXtractor++ documentation for a better explanation of how residuals are computed and how
510 this value affects the model fitting.
512 params_dict[
"modified_chi_squared_scale"] = scale
520 Minimization engine for the model fitting : levmar or gsl
522 params_dict[
"engine"] = engine
527 Base class for all models.
532 class CoordinateModelBase(ModelBase):
534 Base class for positioned models with a flux. It can not be used directly.
538 x_coord : ParameterBase or float
539 X coordinate (in the detection image)
540 y_coord : ParameterBase or float
541 Y coordinate (in the detection image)
542 flux : ParameterBase or float
550 ModelBase.__init__(self)
558 Models a source as a point, spread by the PSF.
562 x_coord : ParameterBase or float
563 X coordinate (in the detection image)
564 y_coord : ParameterBase or float
565 Y coordinate (in the detection image)
566 flux : ParameterBase or float
574 CoordinateModelBase.__init__(self, x_coord, y_coord, flux)
575 global point_source_model_dict
576 point_source_model_dict[self.id] = self
580 Return a human readable representation of the model.
585 If True, include information about the parameters.
592 return 'PointSource[x_coord={}, y_coord={}, flux={}]'.format(
595 return 'PointSource[x_coord={}, y_coord={}, flux={}]'.format(
596 self.x_coord.id, self.y_coord.id, self.flux.id)
601 A model that is constant through all the image.
605 value: ParameterBase or float
606 Value to add to the value of all pixels from the model.
613 ModelBase.__init__(self)
615 global constant_model_dict
616 constant_model_dict[self.id] = self
620 Return a human readable representation of the model.
625 If True, include information about the parameters.
632 return 'ConstantModel[value={}]'.format(self.
value)
634 return 'ConstantModel[value={}]'.format(self.value.id)
639 Base class for the Sersic, Exponential and de Vaucouleurs models. It can not be used directly.
643 x_coord : ParameterBase or float
644 X coordinate (in the detection image)
645 y_coord : ParameterBase or float
646 Y coordinate (in the detection image)
647 flux : ParameterBase or float
649 effective_radius : ParameterBase or float
650 Ellipse semi-major axis, in pixels on the detection image.
651 aspect_ratio : ParameterBase or float
653 angle : ParameterBase or float
654 Ellipse rotation, in radians
657 def __init__(self, x_coord, y_coord, flux, effective_radius, aspect_ratio, angle):
661 CoordinateModelBase.__init__(self, x_coord, y_coord, flux)
669 Model a source with a Sersic profile.
673 x_coord : ParameterBase or float
674 X coordinate (in the detection image)
675 y_coord : ParameterBase or float
676 Y coordinate (in the detection image)
677 flux : ParameterBase or float
679 effective_radius : ParameterBase or float
680 Ellipse semi-major axis, in pixels on the detection image.
681 aspect_ratio : ParameterBase or float
683 angle : ParameterBase or float
684 Ellipse rotation, in radians
685 n : ParameterBase or float
689 def __init__(self, x_coord, y_coord, flux, effective_radius, aspect_ratio, angle, n):
693 SersicModelBase.__init__(self, x_coord, y_coord, flux, effective_radius, aspect_ratio, angle)
695 global sersic_model_dict
696 sersic_model_dict[self.id] = self
700 Return a human readable representation of the model.
705 If True, include information about the parameters.
712 return 'Sersic[x_coord={}, y_coord={}, flux={}, effective_radius={}, aspect_ratio={}, angle={}, n={}]'.format(
715 return 'Sersic[x_coord={}, y_coord={}, flux={}, effective_radius={}, aspect_ratio={}, angle={}, n={}]'.format(
716 self.x_coord.id, self.y_coord.id, self.flux.id, self.effective_radius.id, self.aspect_ratio.id, self.angle.id, self.n.id)
721 Model a source with an exponential profile (Sersic model with an index of 1)
725 x_coord : ParameterBase or float
726 X coordinate (in the detection image)
727 y_coord : ParameterBase or float
728 Y coordinate (in the detection image)
729 flux : ParameterBase or float
731 effective_radius : ParameterBase or float
732 Ellipse semi-major axis, in pixels on the detection image.
733 aspect_ratio : ParameterBase or float
735 angle : ParameterBase or float
736 Ellipse rotation, in radians
739 def __init__(self, x_coord, y_coord, flux, effective_radius, aspect_ratio, angle):
743 SersicModelBase.__init__(self, x_coord, y_coord, flux, effective_radius, aspect_ratio, angle)
744 global exponential_model_dict
745 exponential_model_dict[self.id] = self
749 Return a human readable representation of the model.
754 If True, include information about the parameters.
761 return 'Exponential[x_coord={}, y_coord={}, flux={}, effective_radius={}, aspect_ratio={}, angle={}]'.format(
764 return 'Exponential[x_coord={}, y_coord={}, flux={}, effective_radius={}, aspect_ratio={}, angle={}]'.format(
765 self.x_coord.id, self.y_coord.id, self.flux.id, self.effective_radius.id, self.aspect_ratio.id, self.angle.id)
770 Model a source with a De Vaucouleurs profile (Sersic model with an index of 4)
774 x_coord : ParameterBase or float
775 X coordinate (in the detection image)
776 y_coord : ParameterBase or float
777 Y coordinate (in the detection image)
778 flux : ParameterBase or float
780 effective_radius : ParameterBase or float
781 Ellipse semi-major axis, in pixels on the detection image.
782 aspect_ratio : ParameterBase or float
784 angle : ParameterBase or float
785 Ellipse rotation, in radians
788 def __init__(self, x_coord, y_coord, flux, effective_radius, aspect_ratio, angle):
792 SersicModelBase.__init__(self, x_coord, y_coord, flux, effective_radius, aspect_ratio, angle)
793 global de_vaucouleurs_model_dict
794 de_vaucouleurs_model_dict[self.id] = self
798 Return a human readable representation of the model.
803 If True, include information about the parameters.
810 return 'DeVaucouleurs[x_coord={}, y_coord={}, flux={}, effective_radius={}, aspect_ratio={}, angle={}]'.format(
813 return 'DeVaucouleurs[x_coord={}, y_coord={}, flux={}, effective_radius={}, aspect_ratio={}, angle={}]'.format(
814 self.x_coord.id, self.y_coord.id, self.flux.id, self.effective_radius.id, self.aspect_ratio.id, self.angle.id)
819 Coordinates in right ascension and declination
839 Transform an (X, Y) in pixel coordinates on the detection image to (RA, DEC) in world coordinates.
849 global coordinate_system
850 wc = coordinate_system.image_to_world(cpp.ImageCoordinate(x-1, y-1))
856 Transform an (X, Y) in pixel coordinates on the detection image to astropy SkyCoord.
868 sky_coord = SkyCoord(ra=coord.ra*u.degree, dec=coord.dec*u.degree)
874 Transform a radius in pixels on the detection image to a radius in sky coordinates.
891 Get the separation angle in sky coordinates for two points defined in pixels on the detection image.
902 Separation in degrees
906 return c1.separation(c2).degree
911 Get the position angle in sky coordinates for two points defined in pixels on the detection image.
922 Position angle in degrees, normalized to -/+ 90
926 angle = c1.position_angle(c2).degree
929 return (angle + 90.0) % 180.0 - 90.0
934 Set the global coordinate system. This function is used internally by SourceXtractor++.
936 global coordinate_system
937 coordinate_system = cs
942 Convenience function for generating two dependent parameter with world (alpha, delta) coordinates
943 from image (X, Y) coordinates.
952 ra : DependentParameter
953 dec : DependentParameter
961 >>> x, y = get_pos_parameters()
962 >>> ra, dec = get_world_position_parameters(x, y)
963 >>> add_output_column('mf_ra', ra)
964 >>> add_output_column('mf_dec', dec)
973 Convenience function for generating five dependent parameters, in world coordinates, for the position
974 and shape of a model.
980 radius : ParameterBase
981 angle : ParameterBase
982 ratio : ParameterBase
986 ra : DependentParameter
988 dec : DependentParameter
990 rad : DependentParameter
992 angle : DependentParameter
994 ratio : DependentParameter
995 Aspect ratio. It has to be recomputed as the axis of the ellipse may have different ratios
996 in image coordinates than in world coordinates
1000 >>> flux = get_flux_parameter()
1001 >>> x, y = get_pos_parameters()
1002 >>> radius = FreeParameter(lambda o: o.get_radius(), Range(lambda v, o: (.01 * v, 100 * v), RangeType.EXPONENTIAL))
1003 >>> angle = FreeParameter(lambda o: o.get_angle(), Range((-np.pi, np.pi), RangeType.LINEAR))
1004 >>> ratio = FreeParameter(1, Range((0, 10), RangeType.LINEAR))
1005 >>> add_model(group, ExponentialModel(x, y, flux, radius, ratio, angle))
1006 >>> ra, dec, wc_rad, wc_angle, wc_ratio = get_world_parameters(x, y, radius, angle, ratio)
1007 >>> add_output_column('mf_world_angle', wc_angle)
1012 def get_major_axis(x, y, radius, angle, ratio):
1014 x2 = x + math.cos(angle) * radius
1015 y2 = y + math.sin(angle) * radius
1017 x2 = x + math.sin(angle) * radius * ratio
1018 y2 = y + math.cos(angle) * radius * ratio
1022 def get_minor_axis(x, y, radius, angle, ratio):
1024 x2 = x + math.sin(angle) * radius * ratio
1025 y2 = y + math.cos(angle) * radius * ratio
1027 x2 = x + math.cos(angle) * radius
1028 y2 = y + math.sin(angle) * radius
1032 def wc_rad_func(x, y, radius, angle, ratio):
1033 x2, y2 = get_major_axis(x, y, radius, angle, ratio)
1036 def wc_angle_func(x, y, radius, angle, ratio):
1037 x2, y2 = get_major_axis(x, y, radius, angle, ratio)
1040 def wc_ratio_func(x, y, radius, angle, ratio):
1041 minor_x, minor_y = get_minor_axis(x, y, radius, angle, ratio)
1044 major_x, major_y = get_major_axis(x, y, radius, angle, ratio)
1047 return minor_angle / major_angle
1053 return (ra, dec, wc_rad, wc_angle, wc_ratio)