In [1]:
import numpy as np
import xarray as xr
In [2]:
%load_ext autoreload
%autoreload 2
In [3]:
try:
    import aggutils
except ModuleNotFoundError as e:
    import sys
    sys.path.append('../src')
    import aggutils

not able to load aggregation model: No module named 'aggregation'
In [4]:
from pylab import *
from matplotlib import rc
rc('font',**{'family':'monospace','serif':['Palatino'],'size':10})
fig = lambda **kwargs: plt.figure(figsize=(12,4), **kwargs)
In [5]:
dbloader = lambda f: aggutils.load_agg_db(f)
DS = dict(
    axel = dbloader('../db.axel.nc').isel(uid=slice(None,None,10)),
    pd80mf100 = dbloader('../db.pd80e-6.pdmf100.nc'),
    pd120mf50 = dbloader('../db.pd120e-6.pdmf50.nc'),
)
In [6]:
if 'axel' in DS: # need to hack together variables because they have not been saved in earlier versions
    DS['axel'] = DS['axel'].where(DS['axel'].Nmono!=35, drop=True).where(DS['axel'].Nmono!=350, drop=True)
    if 'l2' not in DS['axel']:
        DS['axel']['l2'] = DS['axel']['aspect']
        DS['axel']['l2'].name = 'aspect'
In [7]:
[ print(k,v) for k,v in DS.items() ]

axel <xarray.Dataset> Size: 695MB
Dimensions:    (uid: 2715057, xyz: 3)
Coordinates:
  * uid        (uid) uint64 22MB 6176477161900751245 ... 17329476830927574502
Dimensions without coordinates: xyz
Data variables: (12/23)
    Nmono      (uid) float64 22MB 2.0 2.0 2.0 2.0 ... 2e+03 2e+03 2e+03 2e+03
    mass       (uid) float64 22MB 8.465e-11 1.194e-09 ... 3.71e-05 3.778e-05
    area2      (uid) float64 22MB 5.8e-09 4.702e-08 ... 0.0001985 0.0001867
    rime       (uid) float64 22MB 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0
    Dmax       (uid) float64 22MB 0.0001084 0.0003199 ... 0.02809 0.02761
    ar0        (uid) float64 22MB 1.154 1.316 1.75 2.8 ... 1.081 1.088 1.02
    ...         ...
    habit      (uid) float64 22MB 0.0 0.0 0.0 0.0 0.0 ... 29.0 29.0 29.0 29.0
    savg_base  (uid) float64 22MB 50.0 50.0 50.0 50.0 ... 490.0 490.0 490.0
    pa_norm    (xyz, uid) float64 65MB 2.629e-05 7.181e-05 ... 0.003098 0.003605
    aspect     (xyz, uid) float64 65MB 1.154 1.316 1.75 ... 1.283 1.369 1.212
    area       (xyz, uid) float64 65MB nan nan nan ... 0.0001985 0.0001867
    l2         (xyz, uid) float64 65MB 1.154 1.316 1.75 ... 1.283 1.369 1.212
pd80mf100 <xarray.Dataset> Size: 963MB
Dimensions:                     (uid: 1501674, xyz: 3, exeyez: 3)
Coordinates:
  * uid                         (uid) uint64 12MB 4718669035105697884 ... 116...
Dimensions without coordinates: xyz, exeyez
Data variables: (12/32)
    savg_base                   (uid) int64 12MB 50 50 50 50 ... 500 500 500 500
    mass                        (uid) float64 12MB 1.708e-10 ... 3.464e-05
    Dmax                        (uid) float64 12MB 0.0001239 ... 0.0326
    area                        (uid, xyz) float64 36MB 3.037e-09 ... 0.0002359
    l1                          (uid, xyz) float64 36MB 0.0001125 ... 0.0318
    l2                          (uid, xyz) float64 36MB 4.25e-05 ... 0.0251
    ...                          ...
    Dmax_parent1                (uid) float64 12MB 3.742e-05 ... 0.0003354
    area_parent1                (uid, xyz) float64 36MB 1.062e-09 ... 2e-08
    l1_parent1                  (uid, xyz) float64 36MB 3.75e-05 ... 0.00015
    l2_parent1                  (uid, xyz) float64 36MB 3.75e-05 ... 0.0003
    pa_norm_parent1             (uid, xyz) float64 36MB 4.868e-05 ... 6.913e-05
    pa_parent1                  (uid, xyz, exeyez) float64 108MB 4.868e-05 .....
Attributes:
    merged_files:  ndb.pd80e-6.pdmf100/aggs.merged.nc
pd120mf50 <xarray.Dataset> Size: 4GB
Dimensions:                     (uid: 5728764, xyz: 3, exeyez: 3)
Coordinates:
  * uid                         (uid) uint64 46MB 10884240525845135471 ... 13...
Dimensions without coordinates: xyz, exeyez
Data variables: (12/32)
    savg_base                   (uid) int64 46MB 50 50 50 50 ... 500 500 500 500
    mass                        (uid) float64 46MB 3.369e-11 ... 4.188e-05
    Dmax                        (uid) float64 46MB 6.135e-05 ... 0.03927
    area                        (uid, xyz) float64 137MB 1.662e-09 ... 0.0002581
    l1                          (uid, xyz) float64 137MB 5e-05 6e-05 ... 0.03928
    l2                          (uid, xyz) float64 137MB 4.5e-05 ... 0.0224
    ...                          ...
    Dmax_parent1                (uid) float64 46MB 4.61e-05 ... 0.0007826
    area_parent1                (uid, xyz) float64 137MB 1.419e-09 ... 1.412e-07
    l1_parent1                  (uid, xyz) float64 137MB 4.5e-05 ... 0.000675
    l2_parent1                  (uid, xyz) float64 137MB 4.25e-05 ... 0.000775
    pa_norm_parent1             (uid, xyz) float64 137MB 1.338e-05 ... 6.426e-05
    pa_parent1                  (uid, xyz, exeyez) float64 412MB 1.338e-05 .....
Attributes:
    merged_files:  ['ndb.pd120e-6.pdmf50/aggs.6459707994.0000.nc', 'ndb.pd120...
Out[7]:
[None, None, None]
In [8]:
from pprint import pprint
mD_mono_db = aggutils.get_monomer_mD_relations(mono_db='../mono*.nc')
pprint(mD_mono_db)
for k,v in mD_mono_db.items():
    plt.plot(*v, label=k, marker='o', linestyle=None)
plt.legend()
mD_mono = aggutils.get_monomer_mD_relations(mono_db='../mono*.nc', interpolate=True)
pprint(mD_mono)
for k,v in mD_mono.items():
    marker = 'x'
    if k in mD_mono_db:
        marker = 'd'
    plt.plot(*v, label=k, marker=marker, linestyle=None, color='.3', alpha=.5)

{'dendrite': array([0.02241413, 2.18473873]),
 'mixDN5': array([0.01313886, 2.06484671]),
 'mixPD5': array([0.08322666, 2.27754602]),
 'mixPN5': array([0.05839223, 2.17948851]),
 'needle': array([0.00605893, 1.91595868]),
 'plate': array([0.47765206, 2.4211984 ])}
{'dendrite': array([0.02241413, 2.18473873]),
 'mixDN1': (0.00690570989060169, 1.9428366848421925),
 'mixDN2': (0.007870837365252613, 1.9697146894629687),
 'mixDN3': (0.008970849023728541, 1.9965926940837448),
 'mixDN4': (0.01022459599556832, 2.0234706987045215),
 'mixDN5': (0.011653564004484923, 2.050348703325298),
 'mixDN6': (0.013282241573700254, 2.0772267079460742),
 'mixDN7': (0.015138539690882224, 2.1041047125668504),
 'mixDN8': (0.01725427012456989, 2.130982717187627),
 'mixDN9': (0.019665690589095102, 2.1578607218084036),
 'mixPD1': (0.030435524829398906, 2.208384694098672),
 'mixPD2': (0.04132756135497623, 2.232030661768164),
 'mixPD3': (0.05611755792361205, 2.255676629437656),
 'mixPD4': (0.07620048713401185, 2.2793225971071482),
 'mixPD5': (0.10347054387798932, 2.3029685647766405),
 'mixPD6': (0.1404998032568778, 2.3266145324461323),
 'mixPD7': (0.1907808152482378, 2.3502605001156245),
 'mixPD8': (0.25905601732577893, 2.3739064677851163),
 'mixPD9': (0.35176503478807797, 2.3975524354546085),
 'mixPN1': (0.00937707336420997, 1.9664826525116845),
 'mixPN2': (0.01451238886885385, 2.0170066248019527),
 'mixPN3': (0.022460038702979405, 2.0675305970922215),
 'mixPN4': (0.03476018614840032, 2.1180545693824895),
 'mixPN5': (0.05379645854800597, 2.1685785416727583),
 'mixPN6': (0.08325786691566694, 2.2191025139630267),
 'mixPN7': (0.12885369391297694, 2.269626486253295),
 'mixPN8': (0.19941988727427823, 2.3201504585435635),
 'mixPN9': (0.3086313650219748, 2.3706744308338323),
 'needle': array([0.00605893, 1.91595868]),
 'plate': array([0.47765206, 2.4211984 ])}
No description has been provided for this image
In [9]:
for i, (k, ds) in enumerate(DS.items()):
    Dmin = 'axel' if k=='axel' else 0
    aggutils.fit_Dnorm(ds, visualize=True, verbose=1, mD_mono=mD_mono, Dmin=Dmin, fignum=(2*i+1,2*i+2))
    plt.suptitle(k)
    plt.tight_layout()

Dmin={'plate': 2e-05, 'dendrite': 2e-05, 'needle': 3e-05, 'mixPN1': 2.9000000000000004e-05, 'mixPN2': 2.8000000000000003e-05, 'mixPN3': 2.7e-05, 'mixPN4': 2.6000000000000002e-05, 'mixPN5': 2.5e-05, 'mixPN6': 2.4e-05, 'mixPN7': 2.3000000000000003e-05, 'mixPN8': 2.2e-05, 'mixPN9': 2.1e-05, 'mixDN1': 2.9000000000000004e-05, 'mixDN2': 2.8000000000000003e-05, 'mixDN3': 2.7e-05, 'mixDN4': 2.6000000000000002e-05, 'mixDN5': 2.5e-05, 'mixDN6': 2.4e-05, 'mixDN7': 2.3000000000000003e-05, 'mixDN8': 2.2e-05, 'mixDN9': 2.1e-05, 'mixPD1': 2e-05, 'mixPD2': 2.0000000000000005e-05, 'mixPD3': 1.9999999999999998e-05, 'mixPD4': 2e-05, 'mixPD5': 2e-05, 'mixPD6': 2e-05, 'mixPD7': 2e-05, 'mixPD8': 2e-05, 'mixPD9': 2e-05}
habit_name=plate           habit_name=plate      eta=[0.47307337] / mu=0.031 sigma=0.155 deviation from lognormal: mae=2.5% p-value: 0.0564
habit_name=needle          habit_name=needle     eta=[0.46215914] / mu=0.033 sigma=0.156 deviation from lognormal: mae=3.9% p-value: 0.0491
habit_name=mixPN1          habit_name=mixPN1     eta=[0.46142706] / mu=0.032 sigma=0.171 deviation from lognormal: mae=2.8% p-value: 0.37
habit_name=mixPN2          habit_name=mixPN2     eta=[0.46091614] / mu=0.032 sigma=0.176 deviation from lognormal: mae=1.7% p-value: 0.505
habit_name=mixPN3          habit_name=mixPN3     eta=[0.46030252] / mu=0.026 sigma=0.171 deviation from lognormal: mae=1.3% p-value: 0.578
habit_name=mixPN5          habit_name=mixPN5     eta=[0.46289552] / mu=0.025 sigma=0.169 deviation from lognormal: mae=1.4% p-value: 0.875
habit_name=mixPN7          habit_name=mixPN7     eta=[0.46655297] / mu=0.026 sigma=0.163 deviation from lognormal: mae=2.1% p-value: 0.147
habit_name=mixPN8          habit_name=mixPN8     eta=[0.47076141] / mu=0.028 sigma=0.167 deviation from lognormal: mae=2.0% p-value: 0.357
habit_name=mixPN9          habit_name=mixPN9     eta=[0.47065738] / mu=0.033 sigma=0.160 deviation from lognormal: mae=2.3% p-value: 0.00336
Dmin={'plate': 0, 'dendrite': 0, 'needle': 0, 'mixPN1': 0.0, 'mixPN2': 0.0, 'mixPN3': 0.0, 'mixPN4': 0.0, 'mixPN5': 0.0, 'mixPN6': 0.0, 'mixPN7': 0.0, 'mixPN8': 0.0, 'mixPN9': 0.0, 'mixDN1': 0.0, 'mixDN2': 0.0, 'mixDN3': 0.0, 'mixDN4': 0.0, 'mixDN5': 0.0, 'mixDN6': 0.0, 'mixDN7': 0.0, 'mixDN8': 0.0, 'mixDN9': 0.0, 'mixPD1': 0.0, 'mixPD2': 0.0, 'mixPD3': 0.0, 'mixPD4': 0.0, 'mixPD5': 0.0, 'mixPD6': 0.0, 'mixPD7': 0.0, 'mixPD8': 0.0, 'mixPD9': 0.0}
habit_name=plate           habit_name=plate      eta=[0.47119623] / mu=0.025 sigma=0.144 deviation from lognormal: mae=3.0% p-value: 0.176
habit_name=needle          habit_name=needle     eta=[0.46345537] / mu=0.038 sigma=0.155 deviation from lognormal: mae=2.2% p-value: 0.0206
habit_name=dendrite        habit_name=dendrite   eta=[0.45915893] / mu=0.008 sigma=0.163 deviation from lognormal: mae=2.6% p-value: 0.252
habit_name=mixPN5          habit_name=mixPN5     eta=[0.46195319] / mu=0.023 sigma=0.168 deviation from lognormal: mae=1.4% p-value: 0.772
habit_name=mixDN5          habit_name=mixDN5     eta=[0.46079297] / mu=0.013 sigma=0.167 deviation from lognormal: mae=2.2% p-value: 0.465
habit_name=mixPD5          habit_name=mixPD5     eta=[0.45831758] / mu=0.003 sigma=0.190 deviation from lognormal: mae=1.9% p-value: 0.0322
Dmin={'plate': 0, 'dendrite': 0, 'needle': 0, 'mixPN1': 0.0, 'mixPN2': 0.0, 'mixPN3': 0.0, 'mixPN4': 0.0, 'mixPN5': 0.0, 'mixPN6': 0.0, 'mixPN7': 0.0, 'mixPN8': 0.0, 'mixPN9': 0.0, 'mixDN1': 0.0, 'mixDN2': 0.0, 'mixDN3': 0.0, 'mixDN4': 0.0, 'mixDN5': 0.0, 'mixDN6': 0.0, 'mixDN7': 0.0, 'mixDN8': 0.0, 'mixDN9': 0.0, 'mixPD1': 0.0, 'mixPD2': 0.0, 'mixPD3': 0.0, 'mixPD4': 0.0, 'mixPD5': 0.0, 'mixPD6': 0.0, 'mixPD7': 0.0, 'mixPD8': 0.0, 'mixPD9': 0.0}
habit_name=plate           habit_name=plate      eta=[0.46121589] / mu=-0.032 sigma=0.131 deviation from lognormal: mae=3.3% p-value: 0.179
habit_name=needle          habit_name=needle     eta=[0.44898227] / mu=-0.033 sigma=0.120 deviation from lognormal: mae=2.0% p-value: 0.033
habit_name=dendrite        habit_name=dendrite   eta=[0.44640013] / mu=-0.058 sigma=0.139 deviation from lognormal: mae=2.7% p-value: 0.014
habit_name=mixPN5          habit_name=mixPN5     eta=[0.4502609] / mu=-0.040 sigma=0.139 deviation from lognormal: mae=1.7% p-value: 0.00308
habit_name=mixDN5          habit_name=mixDN5     eta=[0.44795275] / mu=-0.054 sigma=0.149 deviation from lognormal: mae=1.7% p-value: 0.141
habit_name=mixPD5          habit_name=mixPD5     eta=[0.44785743] / mu=-0.057 sigma=0.172 deviation from lognormal: mae=0.9% p-value: 0.681
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
In [10]:
# Dnorm fit again but this time with Dmax from projected bounding box instead of min bounding sphere
for i, (k, ds) in enumerate(DS.items()):
    if k=='axel': continue # skip because we dont have lengths saved
    dss = ds.copy()
    dss['Dmax'] = aggutils.compute_Dmax_ds(dss)
    Dmin = 'axel' if k=='axel' else 0
    aggutils.fit_Dnorm(dss, visualize=True, verbose=1, mD_mono=mD_mono, Dmin=Dmin, ifit=1, fignum=(2*i+1, 2*i+2))
    plt.suptitle(k)
    plt.tight_layout()

Dmin={'plate': 0, 'dendrite': 0, 'needle': 0, 'mixPN1': 0.0, 'mixPN2': 0.0, 'mixPN3': 0.0, 'mixPN4': 0.0, 'mixPN5': 0.0, 'mixPN6': 0.0, 'mixPN7': 0.0, 'mixPN8': 0.0, 'mixPN9': 0.0, 'mixDN1': 0.0, 'mixDN2': 0.0, 'mixDN3': 0.0, 'mixDN4': 0.0, 'mixDN5': 0.0, 'mixDN6': 0.0, 'mixDN7': 0.0, 'mixDN8': 0.0, 'mixDN9': 0.0, 'mixPD1': 0.0, 'mixPD2': 0.0, 'mixPD3': 0.0, 'mixPD4': 0.0, 'mixPD5': 0.0, 'mixPD6': 0.0, 'mixPD7': 0.0, 'mixPD8': 0.0, 'mixPD9': 0.0}
habit_name=plate           habit_name=plate      eta=[0.4647247] / mu=0.014 sigma=0.161 deviation from lognormal: mae=2.1% p-value: 0.147
habit_name=needle          habit_name=needle     eta=[0.45688864] / mu=0.028 sigma=0.172 deviation from lognormal: mae=1.3% p-value: 0.646
habit_name=dendrite        habit_name=dendrite   eta=[0.45271715] / mu=-0.004 sigma=0.179 deviation from lognormal: mae=1.8% p-value: 0.388
habit_name=mixPN5          habit_name=mixPN5     eta=[0.45534247] / mu=0.012 sigma=0.185 deviation from lognormal: mae=1.0% p-value: 0.24
habit_name=mixDN5          habit_name=mixDN5     eta=[0.45425839] / mu=0.002 sigma=0.182 deviation from lognormal: mae=1.6% p-value: 0.0178
habit_name=mixPD5          habit_name=mixPD5     eta=[0.45177625] / mu=-0.009 sigma=0.203 deviation from lognormal: mae=1.2% p-value: 0.191
Dmin={'plate': 0, 'dendrite': 0, 'needle': 0, 'mixPN1': 0.0, 'mixPN2': 0.0, 'mixPN3': 0.0, 'mixPN4': 0.0, 'mixPN5': 0.0, 'mixPN6': 0.0, 'mixPN7': 0.0, 'mixPN8': 0.0, 'mixPN9': 0.0, 'mixDN1': 0.0, 'mixDN2': 0.0, 'mixDN3': 0.0, 'mixDN4': 0.0, 'mixDN5': 0.0, 'mixDN6': 0.0, 'mixDN7': 0.0, 'mixDN8': 0.0, 'mixDN9': 0.0, 'mixPD1': 0.0, 'mixPD2': 0.0, 'mixPD3': 0.0, 'mixPD4': 0.0, 'mixPD5': 0.0, 'mixPD6': 0.0, 'mixPD7': 0.0, 'mixPD8': 0.0, 'mixPD9': 0.0}
habit_name=plate           habit_name=plate      eta=[0.45456002] / mu=-0.043 sigma=0.148 deviation from lognormal: mae=4.7% p-value: 0.188
habit_name=needle          habit_name=needle     eta=[0.44232075] / mu=-0.046 sigma=0.141 deviation from lognormal: mae=2.8% p-value: 0.00802
habit_name=dendrite        habit_name=dendrite   eta=[0.43980801] / mu=-0.072 sigma=0.155 deviation from lognormal: mae=4.2% p-value: 0.416
habit_name=mixPN5          habit_name=mixPN5     eta=[0.44358062] / mu=-0.053 sigma=0.157 deviation from lognormal: mae=3.1% p-value: 0.356
habit_name=mixDN5          habit_name=mixDN5     eta=[0.4413241] / mu=-0.068 sigma=0.165 deviation from lognormal: mae=3.1% p-value: 0.0629
habit_name=mixPD5          habit_name=mixPD5     eta=[0.44116557] / mu=-0.070 sigma=0.185 deviation from lognormal: mae=1.3% p-value: 0.791
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
In [11]:
ds = DS['pd120mf50'].copy()
ds['Dmax'] = aggutils.compute_Dmax_ds(ds)
aggutils.fit_Dnorm(ds, verbose=1, mD_mono=mD_mono, Dmin=0     , ifit=1, visualize=True, fignum=(1,2))
aggutils.fit_Dnorm(ds, verbose=1, mD_mono=mD_mono, Dmin=120e-6, ifit=1, visualize=True, fignum=(1,3))
aggutils.fit_Dnorm(ds, verbose=1, mD_mono=mD_mono, Dmin=0     , ifit=2, visualize=True, fignum=(1,4))

Dmin={'plate': 0, 'dendrite': 0, 'needle': 0, 'mixPN1': 0.0, 'mixPN2': 0.0, 'mixPN3': 0.0, 'mixPN4': 0.0, 'mixPN5': 0.0, 'mixPN6': 0.0, 'mixPN7': 0.0, 'mixPN8': 0.0, 'mixPN9': 0.0, 'mixDN1': 0.0, 'mixDN2': 0.0, 'mixDN3': 0.0, 'mixDN4': 0.0, 'mixDN5': 0.0, 'mixDN6': 0.0, 'mixDN7': 0.0, 'mixDN8': 0.0, 'mixDN9': 0.0, 'mixPD1': 0.0, 'mixPD2': 0.0, 'mixPD3': 0.0, 'mixPD4': 0.0, 'mixPD5': 0.0, 'mixPD6': 0.0, 'mixPD7': 0.0, 'mixPD8': 0.0, 'mixPD9': 0.0}
habit_name=plate           habit_name=plate      eta=[0.45456002] / mu=-0.043 sigma=0.148 deviation from lognormal: mae=4.7% p-value: 0.343
habit_name=needle          habit_name=needle     eta=[0.44232075] / mu=-0.046 sigma=0.141 deviation from lognormal: mae=2.8% p-value: 0.287
habit_name=dendrite        habit_name=dendrite   eta=[0.43980801] / mu=-0.072 sigma=0.155 deviation from lognormal: mae=4.2% p-value: 0.00991
habit_name=mixPN5          habit_name=mixPN5     eta=[0.44358062] / mu=-0.053 sigma=0.157 deviation from lognormal: mae=3.1% p-value: 0.679
habit_name=mixDN5          habit_name=mixDN5     eta=[0.4413241] / mu=-0.068 sigma=0.165 deviation from lognormal: mae=3.1% p-value: 0.135
habit_name=mixPD5          habit_name=mixPD5     eta=[0.44116557] / mu=-0.070 sigma=0.185 deviation from lognormal: mae=1.3% p-value: 0.841
Dmin={'plate': 0.00012, 'dendrite': 0.00012, 'needle': 0.00012, 'mixPN1': 0.00012000000000000002, 'mixPN2': 0.00012, 'mixPN3': 0.00011999999999999999, 'mixPN4': 0.00012, 'mixPN5': 0.00012, 'mixPN6': 0.00012, 'mixPN7': 0.00012, 'mixPN8': 0.00011999999999999999, 'mixPN9': 0.00012, 'mixDN1': 0.00012000000000000002, 'mixDN2': 0.00012, 'mixDN3': 0.00011999999999999999, 'mixDN4': 0.00012, 'mixDN5': 0.00012, 'mixDN6': 0.00012, 'mixDN7': 0.00012, 'mixDN8': 0.00011999999999999999, 'mixDN9': 0.00012, 'mixPD1': 0.00012000000000000002, 'mixPD2': 0.00012, 'mixPD3': 0.00011999999999999999, 'mixPD4': 0.00012, 'mixPD5': 0.00012, 'mixPD6': 0.00012, 'mixPD7': 0.00012, 'mixPD8': 0.00011999999999999999, 'mixPD9': 0.00012}
habit_name=plate           habit_name=plate      eta=[0.45252379] / mu=-0.035 sigma=0.148 deviation from lognormal: mae=4.8% p-value: 0.259
habit_name=needle          habit_name=needle     eta=[0.4396935] / mu=-0.035 sigma=0.141 deviation from lognormal: mae=2.7% p-value: 0.347
habit_name=dendrite        habit_name=dendrite   eta=[0.43759918] / mu=-0.063 sigma=0.156 deviation from lognormal: mae=4.1% p-value: 0.0434
habit_name=mixPN5          habit_name=mixPN5     eta=[0.44125747] / mu=-0.043 sigma=0.158 deviation from lognormal: mae=2.9% p-value: 0.674
habit_name=mixDN5          habit_name=mixDN5     eta=[0.43899633] / mu=-0.059 sigma=0.166 deviation from lognormal: mae=3.0% p-value: 0.478
habit_name=mixPD5          habit_name=mixPD5     eta=[0.43905412] / mu=-0.061 sigma=0.186 deviation from lognormal: mae=1.3% p-value: 0.0104
Dmin={'plate': 0, 'dendrite': 0, 'needle': 0, 'mixPN1': 0.0, 'mixPN2': 0.0, 'mixPN3': 0.0, 'mixPN4': 0.0, 'mixPN5': 0.0, 'mixPN6': 0.0, 'mixPN7': 0.0, 'mixPN8': 0.0, 'mixPN9': 0.0, 'mixDN1': 0.0, 'mixDN2': 0.0, 'mixDN3': 0.0, 'mixDN4': 0.0, 'mixDN5': 0.0, 'mixDN6': 0.0, 'mixDN7': 0.0, 'mixDN8': 0.0, 'mixDN9': 0.0, 'mixPD1': 0.0, 'mixPD2': 0.0, 'mixPD3': 0.0, 'mixPD4': 0.0, 'mixPD5': 0.0, 'mixPD6': 0.0, 'mixPD7': 0.0, 'mixPD8': 0.0, 'mixPD9': 0.0}
habit_name=plate           habit_name=plate      eta=[1.20399517 0.4665735 ] / mu=-0.005 sigma=0.148 deviation from lognormal: mae=3.7% p-value: 0.59
habit_name=needle          habit_name=needle     eta=[1.34801145 0.4614841 ] / mu=0.015 sigma=0.146 deviation from lognormal: mae=3.9% p-value: 0.148
habit_name=dendrite        habit_name=dendrite   eta=[1.47703118 0.46503256] / mu=0.008 sigma=0.157 deviation from lognormal: mae=4.6% p-value: 0.0787
habit_name=mixPN5          habit_name=mixPN5     eta=[1.34645069 0.46271337] / mu=0.008 sigma=0.160 deviation from lognormal: mae=4.8% p-value: 0.0278
habit_name=mixDN5          habit_name=mixDN5     eta=[1.43684868 0.46474504] / mu=0.006 sigma=0.166 deviation from lognormal: mae=3.9% p-value: 0.235
habit_name=mixPD5          habit_name=mixPD5     eta=[1.43033689 0.46427026] / mu=0.003 sigma=0.186 deviation from lognormal: mae=4.8% p-value: 0.687
Out[11]:
{'plate': array([1.20399517, 0.4665735 ]),
 'needle': array([1.34801145, 0.4614841 ]),
 'dendrite': array([1.47703118, 0.46503256]),
 'mixPN5': array([1.34645069, 0.46271337]),
 'mixDN5': array([1.43684868, 0.46474504]),
 'mixPD5': array([1.43033689, 0.46427026])}
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
In [12]:
def plot_mean_aspects(ds, fignum=None):
    ds['aspect'] = aggutils.compute_aspect_ratio(ds)
    mean_aspects = {}
    for habit in np.unique(ds.habit):
        habit_name = aggutils.habit_codes_inv[habit]
        dsh = ds.where(ds.habit==habit, drop=True)
        mean_aspects[habit_name] = []
        for Nmono in np.unique(dsh.Nmono):
            dss = dsh.where(dsh.Nmono==Nmono, drop=True)
            mean_aspects[habit_name].append([int(Nmono), float(dss.aspect.isel(xyz=2).mean())])
        mean_aspects[habit_name] = np.array(mean_aspects[habit_name])

    fig(num=fignum)
    plt.clf()
    markers = dict(plate='p', needle='X', dendrite='*')
    cmap = matplotlib.colormaps.get_cmap('Spectral').resampled(len(mean_aspects))
    for ihabit, (habit, mean_aspect) in enumerate(mean_aspects.items()):
        pltargs = dict(base=2, marker=markers.get(habit,'.'), color=cmap(ihabit))
        plt.semilogx(mean_aspect[:,0], mean_aspect[:,1], alpha=1., label=f'{habit}', linestyle='-', **pltargs)
        plt.hlines(np.mean(mean_aspect[-3:,1]), 0, mean_aspect[:,0].max(), linestyle='--', color=cmap(ihabit))
    plt.xlabel('Nmono')
    plt.gca().xaxis.set_major_formatter(matplotlib.ticker.ScalarFormatter())
    plt.ylabel('mean aspect ratio')
    plt.legend(loc='lower right')
    plt.ylim(.45,.9)
In [13]:
for k, ds in DS.items():
    print(f"{k} {np.unique(ds.Nmono)=}")
    plot_mean_aspects(ds)
    plt.title(k)

axel np.unique(ds.Nmono)=array([   2.,    3.,    5.,   10.,   20.,   50.,  100.,  200.,  500.,
       1000., 2000.])
pd80mf100 np.unique(ds.Nmono)=array([   2,    3,    4,    6,    8,   11,   16,   23,   32,   45,   64,
         91,  128,  181,  256,  362,  512,  724, 1024, 1448, 2048, 2896])
pd120mf50 np.unique(ds.Nmono)=array([   2,    3,    4,    6,    8,   11,   16,   23,   32,   45,   64,
         91,  128,  181,  256,  362,  512,  724, 1024, 1448, 2048, 2896])
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
In [14]:
def plot_boehms_ellipsoidal_area_ratio_q(ds, fignum=4, aspect_minmax=True):
    Dmax = aggutils.compute_Dmax_ds(ds)
    aspect = aggutils.compute_aspect_ratio(ds, minmax=aspect_minmax).isel(xyz=2)

    ds['boehm_area_ratio'] = ds.area.isel(xyz=2) / (np.pi/4 * Dmax**2)
    ds['boehm_q']          = ds.area.isel(xyz=2) / (np.pi/4 * Dmax**2 * aspect)

    mean_qs = {}
    for habit in np.unique(ds.habit):
        habit_name = aggutils.habit_codes_inv[habit]
        dsh = ds.where(ds.habit==habit, drop=True)
        mean_qs[habit_name] = []
        for Nmono in np.unique(dsh.Nmono):
            dss = dsh.where(dsh.Nmono==Nmono, drop=True)
            mean_qs[habit_name].append([
                int(Nmono),
                float(dss.boehm_q.mean()),
                float(dss.boehm_area_ratio.mean()),
                ])
        mean_qs[habit_name] = np.array(mean_qs[habit_name])

    fig(num=fignum)
    plt.clf()
    cmap = matplotlib.colormaps.get_cmap('Spectral').resampled(len(mean_qs))
    markers = dict(plate='p', needle='X', dendrite='*')
    for ihabit, (habit, mean_q) in enumerate(mean_qs.items()):
        pltargs = dict(base=2, marker=markers.get(habit,'.'), color=cmap(ihabit), alpha=.5)
        plt.semilogx(mean_q[:,0], mean_q[:,2], alpha=.7, label=f"Area ratio: {habit}", linestyle='--', **pltargs)
    for ihabit, (habit, mean_q) in enumerate(mean_qs.items()):
        pltargs = dict(base=2, marker=markers.get(habit,'.'), color=cmap(ihabit))
        plt.semilogx(mean_q[:,0], mean_q[:,1], label=f"Boehms ellipsoidal q: {habit}" , linestyle='-', **pltargs)
    plt.xlabel('Nmono')
    plt.gca().xaxis.set_major_formatter(matplotlib.ticker.ScalarFormatter())
    plt.ylabel(f'mean of area ratio')
    plt.ylim(0.1,.9)
    plt.legend(loc='lower right', fontsize='small')
In [15]:
for k, ds in DS.items():    
    plot_boehms_ellipsoidal_area_ratio_q(ds, fignum=None)
    plt.title(k)
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
In [16]:
def plot_dnorm_vs_aspect_ratio(ds, fignum=7, xmin=.5, xmax=2, cmin=10, aspect_minmax=True):
    eta = aggutils.fit_Dnorm(ds, mD_mono=mD_mono, verbose=0)
    ds['dnorm'] = aggutils.compute_Dnorm(ds, mD_mono=mD_mono, eta_fit=eta)
    aspect = aggutils.compute_aspect_ratio(ds, minmax=aspect_minmax).isel(xyz=2)
    
    fig(num=fignum); plt.clf()

    plt.hist2d(ds.dnorm, aspect, bins=(200,200), cmin=cmin)
    plt.xlabel('Dnorm')
    plt.ylabel(f'horizontal aspect_ratio')
    plt.xlim(xmin, xmax)
    plt.tight_layout()
In [17]:
for k, ds in DS.items():
    plot_dnorm_vs_aspect_ratio(ds, fignum=None, cmin=10)
    plt.title(k)
    plot_dnorm_vs_aspect_ratio(ds, fignum=None, cmin=10, aspect_minmax=False)
    plt.title(f"{k} no aspect_minmax")
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
In [18]:
for k, ds in DS.items():
    if k=='axel': continue # skip because we dont have lengths saved
    fig(); plt.clf()
    maxlen = aggutils.compute_Dmax_ds(ds)
    bx = np.logspace(np.log10(float(maxlen.min())), np.log10(float(maxlen.max())), 100)                                                
    by = np.logspace(np.log10(float(ds.Dmax.min())), np.log10(float(ds.Dmax.max())), 100)
    plt.hist2d(maxlen, ds.Dmax, bins=(bx,by), cmin=1, cmap='Spectral_r',)
    plt.plot(bx,bx,'--',color='black')
    plt.xlabel('Dmax from projected point cloud bounding box')
    plt.ylabel('Dmax from min bounding sphere')
    plt.xscale('log')
    plt.yscale('log')
    plt.title(k)
    plt.colorbar()
No description has been provided for this image
No description has been provided for this image
In [19]:
for k, ds in DS.items():
    fig(); plt.clf()
    by = np.logspace(np.log10(float(ds.Nmono.min())), np.log10(float(ds.Nmono.max())), 100)
    plt.hist2d(ds.dnorm, ds.Nmono, bins=(100,by), cmin=10, cmap='Spectral_r',)
    plt.plot(bx,bx,'--',color='black')
    plt.xlabel('Dnorm')
    plt.ylabel('Nmono')
    plt.yscale('log')
    plt.title(k)
    plt.colorbar()
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image

Joint PDFs¶

In [20]:
mD_mono = aggutils.get_monomer_mD_relations(mono_db='../mono*.nc', normalize_mD_fit_histogram=False, interpolate=True)
for k, ds in DS.items():
    eta_fit = aggutils.fit_Dnorm(ds, verbose=0, mD_mono=mD_mono, Dmin=0e-6)
    ds['Dnorm'] = aggutils.compute_Dnorm(ds=ds, eta_fit=eta_fit, mD_mono=mD_mono)
    ds['aspect'] = aggutils.compute_aspect_ratio(ds, minmax=True)
    Dmax = aggutils.compute_Dmax_ds(ds)
    ds['boehm_area_ratio'] = ds.area.isel(xyz=2) / (np.pi/4 * Dmax**2)
    ds['boehm_q']          = ds.area.isel(xyz=2) / (np.pi/4 * Dmax**2 * ds.aspect.isel(xyz=2))
In [26]:
target_Nmono = 1000

for k, ds in DS.items():
    Nmono = np.unique(ds.Nmono)[abs(np.unique(ds.Nmono)-target_Nmono).argmin()]
    dss = ds.where(ds.Nmono==Nmono, drop=True)
    X, Y = dss.Dnorm, aggutils.compute_aspect_ratio(dss, minmax=False).isel(xyz=2)
    xtransform = [np.log, np.exp] # log-normal for Dnorm
    ytransform = [lambda x:x]*2   # normal distribution for aspect
    data_mean, data_cov, rv = aggutils.fit_joint_pdf(X, Y, visualize=True, xtransform=xtransform, ytransform=ytransform, ylim=(.4,1.1))
    plt.gca().set_title(f'Bivariate Normal Fit: All Habits ({k})')
    print(f'Dataset = {k} r={aggutils.r_value(data_cov)}')
    sampled_pdf = aggutils.random_sample_multivariate_normal_given_x(xtransform[0](X), data_mean, data_cov, size=len(dss.Dnorm))
    plt.figure()
    plt.hist(ytransform[1](sampled_pdf), bins=100, density=True, label=f'sampled bivariate pdf')
    plt.hist(Y, bins=100, density=True, alpha=.5, label=f'data {X.name}')
    plt.xlabel(X.name)
    plt.legend()
    plt.title(f'All Habits ({k})')

    for habit in np.unique(dss.habit):
        habit_name = aggutils.habit_codes_inv[habit]
        dsh = dss.where(dss.habit==habit, drop=True)
        
        X, Y = dsh.Dnorm, aggutils.compute_aspect_ratio(dsh, minmax=False).isel(xyz=2)
        
        data_mean, data_cov, rv = aggutils.fit_joint_pdf(X, Y, visualize=True, xtransform=xtransform, ytransform=ytransform, ylim=(.4,1.1))
        plt.title(f'Bivariate Normal Fit: {habit_name} ({k})')
        
        sampled_pdf = aggutils.random_sample_multivariate_normal_given_x(xtransform[0](X), data_mean, data_cov, size=len(dsh.Dnorm))
        plt.figure()
        plt.hist(ytransform[1](sampled_pdf), bins=100, density=True, label=f'sampled bivariate pdf: {Y.name}')
        plt.hist(Y, bins=100, density=True, alpha=.5, label=f'data: {Y.name}')
        plt.xlabel(X.name)
        plt.legend()
        plt.title(f'{habit_name=} ({k})')
        
        print(f'Dataset = {k} {habit_name=} r-value={aggutils.r_value(data_cov)}')

Dataset = axel r=-0.6181698597173945
Dataset = axel habit_name='plate' r-value=-0.6291586214766198
Dataset = axel habit_name='needle' r-value=-0.595387589132663
Dataset = axel habit_name='mixPN3' r-value=-0.6029742591171319
Dataset = axel habit_name='mixPN5' r-value=-0.6253783918311436
Dataset = axel habit_name='mixPN7' r-value=-0.6380718011814284
Dataset = axel habit_name='mixPN9' r-value=-0.6383519763549302
Dataset = pd80mf100 r=-0.6182458046719742
Dataset = pd80mf100 habit_name='plate' r-value=-0.6500981063342514
Dataset = pd80mf100 habit_name='needle' r-value=-0.5892686168198668

/project/meteo/work/Fabian.Jakub/lib/snowflakes/examples/../src/aggutils.py:1380: RuntimeWarning: More than 20 figures have been opened. Figures created through the pyplot interface (`matplotlib.pyplot.figure`) are retained until explicitly closed and may consume too much memory. (To control this warning, see the rcParam `figure.max_open_warning`). Consider using `matplotlib.pyplot.close()`.
  plt.figure(num=fignum)

Dataset = pd80mf100 habit_name='dendrite' r-value=-0.6271858633274704
Dataset = pd80mf100 habit_name='mixPN5' r-value=-0.5957960454185388
Dataset = pd80mf100 habit_name='mixDN5' r-value=-0.6350392902392707
Dataset = pd80mf100 habit_name='mixPD5' r-value=-0.6171986613846384
Dataset = pd120mf50 r=-0.6241843902202168
Dataset = pd120mf50 habit_name='plate' r-value=-0.6026480792959803
Dataset = pd120mf50 habit_name='needle' r-value=-0.6484936916950024
Dataset = pd120mf50 habit_name='dendrite' r-value=-0.6143262146728602
Dataset = pd120mf50 habit_name='mixPN5' r-value=-0.6402842560353105
Dataset = pd120mf50 habit_name='mixDN5' r-value=-0.6223401694596271
Dataset = pd120mf50 habit_name='mixPD5' r-value=-0.6236862262985707
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
In [27]:
# two lognormal joint pdf fit for Dnorm/boehm_q
target_Nmono = 1000

for k, ds in DS.items():
    Nmono = np.unique(ds.Nmono)[abs(np.unique(ds.Nmono)-target_Nmono).argmin()]
    dss = ds.where(ds.Nmono==Nmono, drop=True)
    xtransform = [np.log, np.exp]
    ytransform = [np.log, np.exp]
    #ytransform = [lambda x:x]*2
    X, Y = dss.Dnorm, dss.boehm_q
    data_mean, data_cov, rv = aggutils.fit_joint_pdf(X, Y, xtransform=xtransform, ytransform=ytransform, visualize=True, ylim=(.1,.7))
    plt.gca().set_title(f'Bivariate Normal Fit: All Habits ({k})')
    print(f'Dataset = {k} r={aggutils.r_value(data_cov)}')
    sampled_pdf = aggutils.random_sample_multivariate_normal_given_x(xtransform[0](X), data_mean, data_cov, size=len(dss.Dnorm))
    plt.figure()
    plt.hist(ytransform[1](sampled_pdf), bins=100, density=True, label=f'sampled bivariate pdf')
    plt.hist(Y, bins=100, density=True, alpha=.5, label=f'data {dss.boehm_q.name}')
    plt.xlabel(X.name)
    plt.legend()
    plt.title(f'All Habits ({k})')

    for habit in np.unique(dss.habit):
        habit_name = aggutils.habit_codes_inv[habit]
        dsh = dss.where(dss.habit==habit, drop=True)
        X, Y = dsh.Dnorm, dsh.boehm_q
        data_mean, data_cov, rv = aggutils.fit_joint_pdf(X, Y, visualize=True, xtransform=xtransform, ytransform=ytransform, ylim=(.1,.7))
        plt.title(f'Bivariate Normal Fit: {habit_name} ({k})')
        sampled_pdf = aggutils.random_sample_multivariate_normal_given_x(xtransform[0](X), data_mean, data_cov, size=len(dsh.Dnorm))
        
        plt.figure()
        plt.hist(ytransform[1](sampled_pdf), bins=100, density=True, label=f'sampled bivariate pdf: {Y.name}')
        plt.hist(Y, bins=100, density=True, alpha=.5, label=f'data: {Y.name}')
        plt.xlabel(X.name)
        plt.legend()
        plt.title(f'{habit_name=} ({k})')
        
        print(f'Dataset = {k} {habit_name=} r-value={aggutils.r_value(data_cov)}')

Dataset = axel r=-0.2982051255368061
Dataset = axel habit_name='plate' r-value=-0.48339308552335686
Dataset = axel habit_name='needle' r-value=-0.15826977302442408
Dataset = axel habit_name='mixPN3' r-value=-0.33604146827635406
Dataset = axel habit_name='mixPN5' r-value=-0.39256343595465754
Dataset = axel habit_name='mixPN7' r-value=-0.44629574863582894
Dataset = axel habit_name='mixPN9' r-value=-0.4704571358511588
Dataset = pd80mf100 r=-0.40635124556230456
Dataset = pd80mf100 habit_name='plate' r-value=-0.6697575105560112
Dataset = pd80mf100 habit_name='needle' r-value=-0.2857208292752462

/project/meteo/work/Fabian.Jakub/lib/snowflakes/examples/../src/aggutils.py:1380: RuntimeWarning: More than 20 figures have been opened. Figures created through the pyplot interface (`matplotlib.pyplot.figure`) are retained until explicitly closed and may consume too much memory. (To control this warning, see the rcParam `figure.max_open_warning`). Consider using `matplotlib.pyplot.close()`.
  plt.figure(num=fignum)

Dataset = pd80mf100 habit_name='dendrite' r-value=-0.49328900996173103
Dataset = pd80mf100 habit_name='mixPN5' r-value=-0.57094347460152
Dataset = pd80mf100 habit_name='mixDN5' r-value=-0.49285129938134997
Dataset = pd80mf100 habit_name='mixPD5' r-value=-0.6134739560180382
Dataset = pd120mf50 r=-0.5698924869550054
Dataset = pd120mf50 habit_name='plate' r-value=-0.7343379816114748
Dataset = pd120mf50 habit_name='needle' r-value=-0.49146455658156124
Dataset = pd120mf50 habit_name='dendrite' r-value=-0.7121287251296111
Dataset = pd120mf50 habit_name='mixPN5' r-value=-0.6253351698516872
Dataset = pd120mf50 habit_name='mixDN5' r-value=-0.6826806785327396
Dataset = pd120mf50 habit_name='mixPD5' r-value=-0.6997510884156257
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image