Millet new

doc/example/millet/SR_diameters.csv

@@ -1,19 +1,19 @@
-position(cm);root_diameter(µm)
-6.846;504.0315
-27.5;507.2155
-14.686;378.7185
-26.8;615.9995
-1.271;327.473
-5.524;388.558
-9.963;465.4815
-14.631;455.2715
-14.631;455.2715
-18.87;491.699
-18.87;491.699
-18.87;491.699
-5.638;336.655
-10.556;286.85
-3.624;347.4265
-6.444;459.23975
-7.985;505.649
-29.287;540.9115
+position(cm),root_diameter(µm)
+6.846,504.0315
+27.5,507.2155
+14.686,378.7185
+26.8,615.9995
+1.271,327.473
+5.524,388.558
+9.963,465.4815
+14.631,455.2715
+14.631,455.2715
+18.87,491.699
+18.87,491.699
+18.87,491.699
+5.638,336.655
+10.556,286.85
+3.624,347.4265
+6.444,459.23975
+7.985,505.649
+29.287,540.9115

src/openalea/hydroroot/conductance.py

@@ -156,7 +156,7 @@ def poiseuille(radius, length, viscosity=1e-3):  # DEPRECATED
 
 
 def compute_k(g, k0 = 300.):
-    """Compute the radial conductance k (:math:`10^{-9}\ m.s^{-1}.MPa^{-1}`) of each vertex of the MTG.
+    r"""Compute the radial conductance k (:math:`10^{-9}\ m.s^{-1}.MPa^{-1}`) of each vertex of the MTG.
 
     .. math::
         k = 2 \\pi r l k0

src/openalea/hydroroot/millet/architecture.py

@@ -247,3 +247,200 @@ def create_randomized_delayed_axis(nid, n, anchors=anchors, **kwds):
     return g
 
 
+#####################################
+# New implementation
+####################################
+import random
+import numpy as np
+
+from openalea.mtg import MTG
+from openalea.mtg.algo import orders
+from openalea.mtg.mtg import fat_mtg
+
+
+def millet_mtg2(
+    g=None,
+    vid=None,
+    seed=None,
+    # seminal
+    primary_length=0.2,          # [m]
+    nude_tip_length=0.03,        # [m]
+    order_max=2,
+    branching_stability=0.8,     # 1=regular spacing, 0=random spacing
+    # NEW: explicit number of laterals on primary
+    n_laterals_primary=5,       # [count]
+    # fixed lateral length
+    lateral_length=0.89,         # [m]
+    # crown
+    nb_crown=2,
+    crown_length=0.06,           # [m]
+    # discretization
+    segment_length=1e-4,         # [m]
+    **kwargs
+):
+    """
+    Millet RSA MTG with an explicit number of laterals (no branching_variability).
+
+    - Primary axis vertex count = int(primary_length/segment_length)
+    - Crown axis vertex count   = int(crown_length/segment_length)
+    - Each lateral axis length  = lateral_length (fixed)
+    - Exactly n_laterals_primary laterals are initiated on the primary axis
+      (unless branchable zone is too short, then it is clipped).
+    """
+
+    # -------------------------
+    # Helper functions that will be used later to build the MTG
+    # -------------------------
+    def steps_from_length(length_m: float):
+        return max(1, int(length_m / segment_length))
+
+    def pick_lateral_indices(n_axis_steps: int, nude_steps: int, n_laterals: int):
+        """
+        Choose exactly n_laterals indices i in [0 .. max_i] where laterals initiate,
+        where max_i excludes the nude tip region near the apex.
+        """
+        max_i = n_axis_steps - 1 - nude_steps
+        if max_i <= 0 or n_laterals <= 0:
+            return []
+
+        # cannot place more laterals than available slots
+        n_laterals = min(int(n_laterals), max_i + 1)
+
+        if branching_stability >= 0.5:
+            # quasi-regular spacing with limited jitter
+            grid = np.linspace(0, max_i, num=n_laterals, endpoint=True)
+            # jitter amplitude in steps (smaller if stability high)
+            jitter_amp = max(1, int((1.0 - branching_stability) * (max_i / max(n_laterals, 1))))
+            idx = []
+            occupied = set()
+            for x in grid:
+                j = int(round(x + np.random.randint(-jitter_amp, jitter_amp + 1)))
+                j = max(0, min(max_i, j))
+                # resolve collisions locally
+                if j in occupied:
+                    for dj in range(1, 50):
+                        for cand in (j - dj, j + dj):
+                            if 0 <= cand <= max_i and cand not in occupied:
+                                j = cand
+                                break
+                        if j not in occupied:
+                            break
+                occupied.add(j)
+                idx.append(j)
+            return sorted(idx)
+        else:
+            # fully random positions without replacement
+            return sorted(random.sample(range(0, max_i + 1), k=n_laterals))
+
+    def create_linear_axis(root_node, axis_length_m: float, label: str, order: int, lateral_idx=None, anchors=None):
+        """
+        Create a successor axis (<). Optionally mark certain indices as anchors.
+        Store position_index = distance-to-tip in meters.
+        """
+        n_steps = steps_from_length(axis_length_m)
+        lateral_idx = set(lateral_idx or [])
+        anchors = anchors if anchors is not None else []
+
+        nid = root_node
+        for i in range(n_steps):
+            nid = nid.add_child(edge_type='<', label=label, order=order, **kwargs)
+            dist_to_tip_steps = (n_steps - 1 - i)
+            nid.position_index = dist_to_tip_steps * segment_length
+            if i in lateral_idx:
+                anchors.append(nid)
+        return nid
+
+    def add_fixed_length_lateral(parent_node, lateral_length_m: float, order: int):
+        """Create a lateral root: '+' then '<' chain of fixed length."""
+        n_steps = steps_from_length(lateral_length_m)
+        cid = parent_node.add_child(edge_type='+', label='Lateral', order=order, **kwargs)
+        nid = cid
+        for i in range(n_steps):
+            nid = nid.add_child(edge_type='<', label='Lateral', order=order, **kwargs)
+            dist_to_tip_steps = (n_steps - 1 - i)
+            nid.position_index = dist_to_tip_steps * segment_length
+        return cid
+
+    def add_crown_axes(collet_ids):
+        for i in range(min(nb_crown, len(collet_ids))):
+            parent = g.node(collet_ids[i])
+            n_steps = steps_from_length(crown_length)
+
+            cid = parent.add_child(edge_type='+', label='Crown', order=0, **kwargs)
+            cid.position_index = (n_steps - 1) * segment_length
+
+            nid = cid
+            for k in range(1, n_steps):
+                nid = nid.add_child(edge_type='<', label='Crown', order=0, **kwargs)
+                nid.position_index = (n_steps - 1 - k) * segment_length
+
+    # -------------------------
+    # Seeding
+    # -------------------------
+    if seed is not None:
+        random.seed(seed)
+        np.random.seed(seed)
+
+    # -------------------------
+    # MTG init + collets: Crown and and primary roots are built from the collet
+    # -------------------------
+    if g is None:
+        g = MTG()
+
+    if vid is None:
+        vid = g.add_component(g.root, label='collet', order=0, **kwargs)
+
+    collet_ids = [vid]
+    for _ in range(nb_crown):
+        vid = g.add_child(vid, label='collet', edge_type='<', order=0, **kwargs)
+        collet_ids.append(vid)
+
+    # -------------------------
+    # Primary axis with explicit laterals: rather than using lateral root density, we reported the number of laterals observed experimentally for the different millet lines
+    # -------------------------
+    n_primary_steps = steps_from_length(primary_length)
+    nude_steps = max(0, steps_from_length(nude_tip_length))
+
+    lateral_indices = pick_lateral_indices(
+        n_axis_steps=n_primary_steps,
+        nude_steps=nude_steps,
+        n_laterals=n_laterals_primary
+    )
+
+    anchors = []
+    seminal_root = g.node(collet_ids[-1])
+    create_linear_axis(
+        root_node=seminal_root,
+        axis_length_m=primary_length,
+        label='Seminal',
+        order=0,
+        lateral_idx=lateral_indices,
+        anchors=anchors
+    )
+
+    # -------------------------
+    # Create laterals (fixed length)
+    # -------------------------
+    hard_vertex_cap = max(5000, n_primary_steps * 50)
+    while anchors and len(g) < hard_vertex_cap:
+        nid = anchors.pop(0)
+        current_order = getattr(nid, "order", 0)
+        if current_order >= order_max:
+            continue
+        add_fixed_length_lateral(
+            parent_node=nid,
+            lateral_length_m=float(lateral_length),
+            order=current_order + 1
+        )
+
+    # -------------------------
+    # Crown roots
+    # -------------------------
+    add_crown_axes(collet_ids)
+
+    # -------------------------
+    # Finalize
+    # -------------------------
+    g = fat_mtg(g)
+    g.properties()['order'] = orders(g, scale=g.max_scale()) # add order property
+    return g

src/openalea/hydroroot/millet/conductance.py

@@ -5,7 +5,7 @@
 """
 
 import numpy as np
-from scipy.interpolate import UnivariateSpline
+#from scipy.interpolate import UnivariateSpline
 
 def fit_property_from_spline(g, spline, prop_in, prop_out): 
     """ compute a property from another one using a spline transformation.
@@ -24,3 +24,32 @@ def fit_property_from_spline(g, spline, prop_in, prop_out):
 
     return g
 
+
+def compute_K_from_laws(g, seminal_axial_conductivity_law, crown_axial_conductivity_law, lateral_axial_conductivity_law):
+    """
+
+    compute the axial conductance K versus the vertex position according to some laws and to the root types: crown, seminal,
+    laterals
+
+    :param g: (MTG)
+
+    """
+    K={}
+
+
+    positions= g.property('relative_position')
+    orders = g.property('order')
+
+    for vid in g.vertices_iter(g.max_scale()):
+            if g.label(vid)=='Crown':
+                K[vid] = crown_axial_conductivity_law(positions[vid])
+            elif g.label(vid) == 'Seminal' :
+                K[vid] = seminal_axial_conductivity_law(positions[vid])
+            elif g.label(vid) == 'Lateral' :
+                    K[vid] = lateral_axial_conductivity_law(positions[vid])
+            else: # for the collet and other unknown tpes
+                K[vid] = seminal_axial_conductivity_law(positions[vid])
+
+    g.properties()['K'] = K
+
+    return g

src/openalea/hydroroot/millet/display_millet.py

@@ -19,7 +19,7 @@
 from openalea.plantgl import *
 import openalea.plantgl.all as pgl
 
-from openalea.hydroroot.millet import turtle, law
+from openalea.hydroroot.millet import turtle_millet as turtle, law
 from openalea.hydroroot.display import my_colormap
 ########################################################################################