From 116d976d3b6b46a98e491f0fca4b2aab803f5997 Mon Sep 17 00:00:00 2001 From: Ross Beyer Date: Tue, 31 Jan 2023 16:18:11 -0800 Subject: [PATCH 1/2] feat(crater.determine_production_function()): Added function. --- src/python/synthterrain/crater/__init__.py | 16 ++++++++++------ 1 file changed, 10 insertions(+), 6 deletions(-) diff --git a/src/python/synthterrain/crater/__init__.py b/src/python/synthterrain/crater/__init__.py index b3021ae..e3bfbfe 100644 --- a/src/python/synthterrain/crater/__init__.py +++ b/src/python/synthterrain/crater/__init__.py @@ -42,12 +42,7 @@ def synthesize( synthesized from the input parameters. """ if production_fn is None: - if crater_dist.a >= 10: - production_fn = functions.NPF(a=crater_dist.a, b=crater_dist.b) - elif crater_dist.b <= 2.76: - production_fn = functions.Grun(a=crater_dist.a, b=crater_dist.b) - else: - production_fn = functions.GNPF(a=crater_dist.a, b=crater_dist.b) + production_fn = determine_production_function(crater_dist.a, crater_dist.b) logger.info(f"Production function is {production_fn.__class__}") # Get craters @@ -90,6 +85,15 @@ def synthesize( return df +def determine_production_function(a: float, b: float): + if a >= 10: + return functions.NPF(a, b) + elif b <= 2.76: + return functions.Grun(a, b) + else: + return functions.GNPF(a, b) + + def random_points(poly: Polygon, num_points: int): """Returns two lists, the first being the x coordinates, and the second being the y coordinates, each *num_points* long that represent From 9530e54c409e47980dc2e7197a000c4d9bc71dc2 Mon Sep 17 00:00:00 2001 From: Ross Beyer Date: Tue, 31 Jan 2023 16:52:22 -0800 Subject: [PATCH 2/2] feat(crater.age.estimate_age_by_bin() and cli_convert): Allowed production and equilibrium functions to be None which runs diffusion for the age of the solar system, and also adjusted for the case where there are many craters with only a single diameter. --- src/python/synthterrain/crater/age.py | 54 ++++++++++++------- src/python/synthterrain/crater/cli_convert.py | 26 ++++++--- tests/python/crater/test_age.py | 17 +++++- 3 files changed, 71 insertions(+), 26 deletions(-) diff --git a/src/python/synthterrain/crater/age.py b/src/python/synthterrain/crater/age.py index 54ae3ba..4a9475d 100644 --- a/src/python/synthterrain/crater/age.py +++ b/src/python/synthterrain/crater/age.py @@ -70,9 +70,9 @@ def estimate_age(diameter, dd, max_age): def estimate_age_by_bin( df, - pd_csfd, - eq_csfd, num=50, + pd_csfd=None, + eq_csfd=None, ) -> pd.DataFrame: """ Returns a pandas DataFrame identical to the input *df* but with the @@ -85,25 +85,38 @@ def estimate_age_by_bin( scale boundaries (using the numpy geomspace() function) between the maximum and minimum diameter values. - Then, the center diameter of each bin has its equilibrium age determined via - equilibrium_ages() and diffuse_d_over_D() is run to evaluate the d/D ratio - of that crater over its lifetime. Then, for each crater in the bin, its d/D - is compared to the d/D ratios from the diffusion run, and an estimated age - is assigned. + Then, the center diameter of each bin has diffuse_d_over_D() + run to evaluate the d/D ratio of that crater over the lifetime + of the solar system (if *pd_csfd* and *eq_csfd* are both None), or + the diffusion calculation is run for the appropriate equilibrium age + of each crater diameter. + + Then, for each crater in the bin, its d/D is compared to the d/D ratios from the + diffusion run, and an estimated age is assigned. """ if df.empty: raise ValueError("The provided dataframe has no rows.") logger.info("estimate_age_by_bin start.") - bin_edges = np.geomspace( - df["diameter"].min(), df["diameter"].max(), num=num + 1 - ) - # logger.info(f"{df.shape[0]} craters") - logger.info( - f"Divided the craters into {num} diameter bins (not all bins may have " - "craters)" + if df["diameter"].min() == df["diameter"].max(): + bin_edges = 1 + total_bins = 1 + logger.info("The craters are all the same size.") + else: + bin_edges = np.geomspace( + df["diameter"].min(), df["diameter"].max(), num=num + 1 + ) + total_bins = num + # logger.info(f"{df.shape[0]} craters") + logger.info( + f"Divided the craters into {num} diameter bins (not all bins may have " + "craters)" + ) + + df["diameter_bin"] = pd.cut( + df["diameter"], bins=bin_edges, include_lowest=True, ) - df["diameter_bin"] = pd.cut(df["diameter"], bins=bin_edges, include_lowest=True) + # df["equilibrium_age"] = equilibrium_ages(df["diameter"], pd_csfd, eq_csfd) df["age"] = 0 @@ -111,7 +124,7 @@ def estimate_age_by_bin( df["diameter_bin"].value_counts(sort=False).items() ): logger.info( - f"Processing bin {i}/{num}, interval: {interval}, count: {count}" + f"Processing bin {i + 1}/{total_bins}, interval: {interval}, count: {count}" ) if count == 0: @@ -119,16 +132,19 @@ def estimate_age_by_bin( fresh_dd = stopar_fresh_dd(interval.mid) - eq_age = equilibrium_age(interval.mid, pd_csfd, eq_csfd) + if pd_csfd is not None and eq_csfd is not None: + age = equilibrium_age(interval.mid, pd_csfd, eq_csfd) + else: + age = 4.5e9 dd_rev_list = list(reversed(diffuse_d_over_D( interval.mid, - eq_age, + age, start_dd_adjust=fresh_dd, return_steps=True ))) nsteps = len(dd_rev_list) - years_per_step = eq_age / nsteps + years_per_step = age / nsteps def guess_age(dd): age_step = bisect.bisect_left(dd_rev_list, dd) diff --git a/src/python/synthterrain/crater/cli_convert.py b/src/python/synthterrain/crater/cli_convert.py index 8ff346a..12271c5 100644 --- a/src/python/synthterrain/crater/cli_convert.py +++ b/src/python/synthterrain/crater/cli_convert.py @@ -40,6 +40,15 @@ def arg_parser(): "craters may still yield a zero age if the d/D ratio was large relative " "to the diameter, indicating a very fresh crater." ) + parser.add_argument( + "--full_age", + action="store_true", + help="Ignored unless --estimate_ages is also given. When provided, it will " + "cause the diffusion calculation to run for the age of the solar system " + "instead of just the equilibrium age for each crater size. This may " + "provide improved age estimates, but could also cause longer run times. " + "Please use with caution." + ) parser.add_argument( "infile", type=Path, @@ -77,17 +86,22 @@ def main(): ) if args.estimate_ages: - a = df["diameter"].min() - b = df["diameter"].max() - pd_func = crater_func.GNPF(a=a, b=b) - eq_func = crater_func.VIPER_Env_Spec(a=a, b=b) + if args.full_age: + pd_csfd = None + eq_csfd = None + else: + a = df["diameter"].min() + b = df["diameter"].max() + pd_csfd = crater.determine_production_function(a=a, b=b).csfd + eq_csfd = crater_func.VIPER_Env_Spec(a=a, b=b).csfd + try: if "age" in df.columns: df[df["age"] == 0] = estimate_age_by_bin( - df[df["age"] == 0], pd_func.csfd, eq_func.csfd, num=50 + df[df["age"] == 0], 50, pd_csfd, eq_csfd ) else: - df = estimate_age_by_bin(df, pd_func.csfd, eq_func.csfd, num=50) + df = estimate_age_by_bin(df, 50, pd_csfd, eq_csfd) except ValueError: logger.error( "The provided file has no craters with an age of zero." diff --git a/tests/python/crater/test_age.py b/tests/python/crater/test_age.py index 90b05e3..0311744 100644 --- a/tests/python/crater/test_age.py +++ b/tests/python/crater/test_age.py @@ -50,9 +50,9 @@ def test_estimate_age_by_bin(self): ) df_out = age.estimate_age_by_bin( df, + 50, # the bin size can have a real impact here. pd_func.csfd, eq_func.csfd, - num=50, # the bin size can have a real impact here. ) age_series = pd.Series([ @@ -61,3 +61,18 @@ def test_estimate_age_by_bin(self): 8072614], name="age") pd.testing.assert_series_equal(age_series, df_out["age"]) + + + df2 = pd.DataFrame(data={ + 'diameter': [100., 100., 100., 100.], + 'd/D': [0.01, 0.06, 0.10, 0.17] + }) + + df2_out = age.estimate_age_by_bin( + df2, + 50, # With only one diameter, num is irrelevant + ) + + age2_series = pd.Series([4500000000, 4500000000, 1388888888, 0], name="age") + + pd.testing.assert_series_equal(age2_series, df2_out["age"]) \ No newline at end of file