Merge pull request #261 from Henry-Best-01/agn_sources

Agn sources

slsim/Sources/agn.py

@@ -128,15 +128,6 @@ def get_mean_mags(self, bands):
         return magnitudes
 
 
-# Include a basic driving light curve to illustrate variability between bands
-# Any light curve can be inserted into the "intrinsic_light_curve" keyword of
-# "kwargs_agn_model". Here we include a placeholder to include (relatively
-# nonsensical) variability without user input.
-basic_light_curve = {
-    "MJD": np.linspace(1, 100, 100),
-    "ps_mag_intrinsic": np.sin(np.linspace(1, 100, 100) * np.pi / 10),
-}
-
 # This dictionary is designed to set the boundaries to draw random parameters from.
 # The bounds of any keys may be redefined using an "input_agn_bounds_dict".
 # :key black_hole_mass_exponent_bounds: mass of SMBH as log_(10)(M_{BH}/M_{sun})
@@ -155,6 +146,8 @@ def get_mean_mags(self, bands):
 #   light curves across all bands. The simplest case os to provide a list of variable
 #   light curves directly, but this will also work with other variability choices using
 #   Source.SourceVariability.variability.Variability(variability_model)
+# :key intrinsic_light_curve: List of light curve objects to randomly choose from. If
+#   None, then a bending power law signal will randomly be generated for a 1000 day period.
 agn_bounds_dict = {
     "black_hole_mass_exponent_bounds": (6.0, 10.0),
     "black_hole_spin_bounds": (-0.997, 0.997),
@@ -163,7 +156,7 @@ def get_mean_mags(self, bands):
     "eddington_ratio_bounds": (0.01, 0.3),
     "supported_disk_models": ["thin_disk"],
     "driving_variability": ["light_curve"],
-    "intrinsic_light_curve": [basic_light_curve],
+    "intrinsic_light_curve": None,
 }
 
 
@@ -237,6 +230,7 @@ def RandomAgn(
     # Check if there was a provided variabilty model.
     # If not, populate driving variability with a simple model
     if agn_driving_variability_model is None:
+
         # Check if other driving variabilities were inserted into bounds dict and randomize
         random_variability_type = random.uniform(
             low=0, high=len(agn_bounds_dict["driving_variability"])
@@ -245,20 +239,48 @@ def RandomAgn(
             "driving_variability"
         ][int(random_variability_type)]
 
-        # Check if other light curves were inserted into bounds dict and randomize
-        random_light_curve_index = random.uniform(
-            low=0, high=len(agn_bounds_dict["intrinsic_light_curve"])
-        )
-        random_light_curve = agn_bounds_dict["intrinsic_light_curve"][
-            int(random_light_curve_index)
-        ]
+        # Check if a list of other light curves were inserted into bounds dict and randomize
+        if input_agn_bounds_dict["intrinsic_light_curve"] is not None:
+
+            random_light_curve_index = random.uniform(
+                low=0, high=len(input_agn_bounds_dict["intrinsic_light_curve"])
+            )
+            random_light_curve = input_agn_bounds_dict["intrinsic_light_curve"][
+                int(random_light_curve_index)
+            ]
 
-        # Set magnitude
-        random_light_curve["ps_mag_intrinsic"] += known_mag
+            # Set magnitude
+            random_light_curve["ps_mag_intrinsic"] += known_mag
+
+            # Define the driving variability model as the light curve to pass into AGN object
+            agn_driving_variability_model = "light_curve"
+            agn_driving_kwargs_variability = random_light_curve
+
+        # If not, generate a bending power law signal from reasonable parameters
+        else:
+            if lightcurve_time is None:
+                length_of_required_light_curve = 1000
+                lightcurve_time = np.linspace(
+                    0,
+                    length_of_required_light_curve - 1,
+                    length_of_required_light_curve,
+                )
+
+            length_of_required_light_curve = np.max(lightcurve_time) - np.min(
+                lightcurve_time
+            )
 
-        # Define the driving variability model as the light curve to pass into AGN object
-        agn_driving_variability_model = "light_curve"
-        agn_driving_kwargs_variability = random_light_curve
+            low_freq_slope = random.uniform(0, 2.0)
+            random_driving_signal_kwargs = {
+                "length_of_light_curve": length_of_required_light_curve,
+                "time_resolution": 1,
+                "log_breakpoint_frequency": random.uniform(-0.5, -2.5),
+                "low_frequency_slope": low_freq_slope,
+                "high_frequency_slope": random.uniform(low_freq_slope, 4.0),
+                "standard_deviation": random.uniform(0.1, 1.0),
+            }
+            agn_driving_variability_model = "bending_power_law"
+            agn_driving_kwargs_variability = random_driving_signal_kwargs
 
     # Define initial speclite filter to be known band
     kwargs_agn_model["speclite_filter"] = known_band

slsim/Util/astro_util.py

@@ -686,7 +686,7 @@ def generate_signal(
         (fourier_transform, fourier_transform[-2:0:-1].conjugate())
     )
     generated_light_curve = ifft(fourier_transform)[
-        : length_of_light_curve // time_resolution
+        : int(length_of_light_curve / time_resolution)
     ]
     if normal_magnitude_variance is False:
         amplitude_baseline = magnitude_to_amplitude(mean_magnitude, zero_point_mag)
@@ -761,7 +761,7 @@ def generate_signal_from_bending_power_law(
     :return: Two arrays, the time_array and the magnitude_array of the variability.
     """
     time_array = np.linspace(
-        0, length_of_light_curve - 1, length_of_light_curve // time_resolution
+        0, length_of_light_curve - 1, int(length_of_light_curve / time_resolution)
     )
     magnitude_array = generate_signal(
         length_of_light_curve,
@@ -816,7 +816,7 @@ def generate_signal_from_generic_psd(
         variability.
     """
     time_array = np.linspace(
-        0, length_of_light_curve - 1, length_of_light_curve // time_resolution
+        0, length_of_light_curve - 1, int(length_of_light_curve / time_resolution)
     )
     magnitude_array = generate_signal(
         length_of_light_curve,

tests/test_Source/test_agn.py

@@ -200,6 +200,28 @@ def test_random_agn():
         input_agn_bounds_dict=input_agn_bounds_dict,
     )
 
+    light_curve_1 = {
+        "MJD": np.asarray([0, 1, 2, 3, 4, 5]),
+        "ps_mag_intrinsic": np.asarray([1, 0, -1, 0, 1, 0]),
+    }
+    light_curve_2 = {
+        "MJD": np.asarray([0, 1, 2, 3, 4, 5]),
+        "ps_mag_intrinsic": np.asarray([1, 0, 1, 0, 1, 0]),
+    }
+    input_agn_bounds_dict["intrinsic_light_curve"] = [light_curve_1, light_curve_2]
+
+    random_agn_3 = RandomAgn(
+        i_band_string,
+        i_band_mag,
+        redshift,
+        random_seed=1,
+        lightcurve_time=lightcurve_time,
+        input_agn_bounds_dict=input_agn_bounds_dict,
+    )
+
+    # test initialization with minimal information
+    RandomAgn(i_band_string, i_band_mag, redshift)
+
     # Test that we raise a warning in RandomAgn when no time axis is input
     with pytest.raises(ValueError):
         RandomAgn(
@@ -214,6 +236,7 @@ def test_random_agn():
 
     # check that a random value from the range [8.0, 8.0) must return 8.0
     assert random_agn_2.kwargs_model["black_hole_mass_exponent"] == 8.0
+    assert random_agn_3.kwargs_model["black_hole_mass_exponent"] == 8.0
 
     # check the inclination is on range [0, 45)
     assert random_agn_2.kwargs_model["inclination_angle"] < 45.1