```python
import bpy, math, random
from mathutils import Vector
from math import sin, cos, pi, sqrt

random.seed(18)

# ---------- scene ----------
bpy.ops.object.select_all(action='SELECT')
bpy.ops.object.delete(use_global=False)
for d in (bpy.data.meshes, bpy.data.curves, bpy.data.materials, bpy.data.cameras, bpy.data.lights):
    pass

world = bpy.context.scene.world or bpy.data.worlds.new("Amber Forest World")
bpy.context.scene.world = world
world.use_nodes = True
world.node_tree.nodes["Background"].inputs["Color"].default_value = (0.055, 0.012, 0.006, 1)
world.node_tree.nodes["Background"].inputs["Strength"].default_value = 0.22

scene = bpy.context.scene
scene.render.engine = 'BLENDER_EEVEE_NEXT'
scene.render.resolution_x = 1536
scene.render.resolution_y = 1536
scene.render.resolution_percentage = 100
scene.render.image_settings.file_format = 'PNG'
scene.view_settings.look = 'AgX - Medium High Contrast'
scene.render.film_transparent = False

# ---------- materials ----------
def mat_principled(name, color, rough=0.5, emission=None, estr=0.0, metallic=0.0, subs=0.0):
    m = bpy.data.materials.new(name)
    m.diffuse_color = (*color, 1)
    m.use_nodes = True
    bs = m.node_tree.nodes.get("Principled BSDF")
    bs.inputs["Base Color"].default_value = (*color, 1)
    bs.inputs["Roughness"].default_value = rough
    bs.inputs["Metallic"].default_value = metallic
    if "Subsurface Weight" in bs.inputs:
        bs.inputs["Subsurface Weight"].default_value = subs
    if emission:
        if "Emission Color" in bs.inputs:
            bs.inputs["Emission Color"].default_value = (*emission, 1)
            bs.inputs["Emission Strength"].default_value = estr
        else:
            bs.inputs["Emission"].default_value = (*emission, 1)
            bs.inputs["Emission Strength"].default_value = estr
    return m

def mat_emission(name, color, strength):
    m = bpy.data.materials.new(name)
    m.use_nodes = True
    nt = m.node_tree
    nt.nodes.clear()
    out = nt.nodes.new("ShaderNodeOutputMaterial")
    em = nt.nodes.new("ShaderNodeEmission")
    em.inputs["Color"].default_value = (*color, 1)
    em.inputs["Strength"].default_value = strength
    nt.links.new(em.outputs[0], out.inputs[0])
    return m

capmat = mat_principled("Warm leathery mushroom flesh", (0.33,0.075,0.025), 0.72, (0.10,0.012,0.002), 0.18, subs=0.12)
capedge = mat_principled("Amber cap edge", (0.48,0.13,0.035), 0.62, (0.18,0.025,0.004), 0.22, subs=0.18)
underside = mat_principled("Moist translucent underside", (0.025,0.16,0.16), 0.4, (0.01,0.12,0.14), 0.45, subs=0.22)
teal = mat_emission("Bioluminescent reaction diffusion folds", (0.01,0.72,0.68), 4.2)
teal_dim = mat_principled("Wet teal channel sides", (0.01,0.19,0.20), 0.33, (0.005,0.25,0.25), 1.8)
bark = mat_principled("Ochre old bark", (0.16,0.045,0.012), 0.88)
bark2 = mat_principled("Copper bark highlights", (0.33,0.095,0.018), 0.82)
wood = mat_principled("Treehouse timber", (0.23,0.07,0.018), 0.7)
roofmat = mat_principled("Rust orange roofs", (0.40,0.075,0.018), 0.66)
leafmat = mat_principled("Copper fern leaves", (0.16,0.24,0.035), 0.8)
mossmat = mat_principled("Moss", (0.07,0.16,0.035), 0.95)
soilmat = mat_principled("Dark forest soil", (0.055,0.018,0.008), 1.0)
windowmat = mat_emission("Warm window light", (1.0,0.22,0.035), 8)
rope_mat = mat_principled("Rope", (0.24,0.10,0.035), 0.9)
sporemat = mat_emission("Floating amber spores", (1.0,0.32,0.06), 3.0)

def assign(obj, mat):
    obj.data.materials.append(mat)
    return obj

def cube(name, loc, scale, mat, bevel=0.0):
    bpy.ops.mesh.primitive_cube_add(location=loc)
    o=bpy.context.object; o.name=name; o.dimensions=scale
    bpy.ops.object.transform_apply(location=False, rotation=False, scale=True)
    assign(o,mat)
    if bevel:
        mod=o.modifiers.new("soft worn edges","BEVEL"); mod.width=bevel; mod.segments=2
    return o

def cyl(name, loc, radius, depth, mat, vertices=12, rot=None):
    bpy.ops.mesh.primitive_cylinder_add(vertices=vertices, radius=radius, depth=depth, location=loc, rotation=rot or (0,0,0))
    o=bpy.context.object; o.name=name; assign(o,mat); return o

def curve_obj(name, pts, bevel, mat, cyclic=False):
    cu=bpy.data.curves.new(name,'CURVE'); cu.dimensions='3D'; cu.resolution_u=2
    cu.bevel_depth=bevel; cu.bevel_resolution=2; cu.resolution_u=2
    sp=cu.splines.new('POLY'); sp.points.add(len(pts)-1)
    for p,co in zip(sp.points,pts): p.co=(*co,1)
    sp.use_cyclic_u=cyclic
    ob=bpy.data.objects.new(name,cu); bpy.context.collection.objects.link(ob); assign(ob,mat)
    return ob

# ---------- iterative reaction diffusion field ----------
N=88
u=[[1.0 for x in range(N)] for y in range(N)]
v=[[0.0 for x in range(N)] for y in range(N)]
for _ in range(95):
    nu=[[0.0]*N for _ in range(N)]
    nv=[[0.0]*N for _ in range(N)]
    for y in range(1,N-1):
        for x in range(1,N-1):
            lapu=(u[y][x-1]+u[y][x+1]+u[y-1][x]+u[y+1][x]-4*u[y][x])
            lapv=(v[y][x-1]+v[y][x+1]+v[y-1][x]+v[y+1][x]-4*v[y][x])
            uv=u[y][x]*v[y][x]*v[y][x]
            nu[y][x]=u[y][x]+0.16*lapu-uv+0.055*(1-u[y][x])
            nv[y][x]=v[y][x]+0.08*lapv+uv-0.062*v[y][x]
    u,v=nu,nv
def field(x,y):
    ix=max(1,min(N-2,int((x+8)/16*(N-1))))
    iy=max(1,min(N-2,int((y+8)/16*(N-1))))
    return v[iy][ix]

# ---------- mushroom cap closed irregular shell ----------
rings=9; seg=96
verts=[]; faces=[]
for j in range(rings+1):
    rr=7.5*j/rings
    for i in range(seg):
        a=2*pi*i/seg
        wob=1+0.055*sin(5*a+0.7)+0.035*sin(11*a)
        r=rr*wob
        ztop=5.35 + 1.28*(1-(rr/7.5)**1.55) + 0.10*sin(3*a+rr*.6) + 0.07*sin(8*a)
        verts.append((r*cos(a),r*sin(a),ztop))
# underside rings
under_start=len(verts)
for j in range(rings+1):
    rr=7.5*j/rings
    for i in range(seg):
        a=2*pi*i/seg
        wob=1+0.055*sin(5*a+0.7)+0.035*sin(11*a)
        r=rr*wob
        z=5.13 + 0.11*(1-(rr/7.5)**1.4) + 0.03*sin(7*a+rr)
        verts.append((r*cos(a),r*sin(a),z))
for j in range(rings):
    for i in range(seg):
        n=(i+1)%seg
        faces.append((j*seg+i,j*seg+n,(j+1)*seg+n,(j+1)*seg+i))
        a=under_start+j*seg+i; b=under_start+j*seg+n
        c=under_start+(j+1)*seg+n; d=under_start+(j+1)*seg+i
        faces.append((a,d,c,b))
for i in range(seg):
    n=(i+1)%seg
    faces.append((rings*seg+i,under_start+rings*seg+i,under_start+rings*seg+n,rings*seg+n))
me=bpy.data.meshes.new("single monumental cap mesh"); me.from_pydata(verts,[],faces); me.update()
cap=bpy.data.objects.new("Colossal flattened oyster cap",me); bpy.context.collection.objects.link(cap); assign(cap,capmat)
bev=cap.modifiers.new("soft biological edges","BEVEL"); bev.width=.10; bev.segments=3
# underside base separate for GN-driven procedural structure
bpy.ops.mesh.primitive_circle_add(vertices=96, radius=7.38, fill_type='TRIFAN', location=(0,0,5.12))
base=bpy.context.object; base.name="Reaction diffusion substrate"; assign(base,underside)
# geometry nodes modifier, with procedural named network
ng=bpy.data.node_groups.new("Reaction Diffusion Growth Simulation","GeometryNodeTree")
ng.interface.new_socket(name="Geometry",in_out='INPUT',socket_type='NodeSocketGeometry')
ng.interface.new_socket(name="Geometry",in_out='OUTPUT',socket_type='NodeSocketGeometry')
gi=ng.nodes.new("NodeGroupInput"); go=ng.nodes.new("NodeGroupOutput")
sm=ng.nodes.new("GeometryNodeSetMaterial"); sm.inputs["Material"].default_value=underside
ng.links.new(gi.outputs["Geometry"],sm.inputs["Geometry"]); ng.links.new(sm.outputs["Geometry"],go.inputs["Geometry"])
mod=base.modifiers.new("Iterative reaction-diffusion geometry driver","NODES"); mod.node_group=ng

# ---------- branching folds and channels ----------
# broad radial growth paths, increasingly dense toward edge
for k in range(30):
    a=2*pi*k/30 + random.uniform(-.04,.04)
    pts=[]
    phase=random.random()*8
    for q in range(18):
        r=.18 + 7.05*q/17
        ang=a + .10*sin(q*.55+phase) + .035*sin(q*1.8+phase)
        # field perturbs the biological trajectory
        fx=field(r*cos(ang),r*sin(ang))
        ang += (fx-.12)*.20
        z=5.22 + .035*sin(q*.8+phase)
        pts.append((r*cos(ang),r*sin(ang),z))
    curve_obj("primary branching fold",pts,.075,teal_dim)
    # branched veins
    for q in range(5,16,3):
        if random.random()<.9:
            p=pts[q]; aa=a+random.uniform(-.7,.7)
            sub=[p]
            for t in range(1,6):
                r=math.hypot(p[0],p[1])+t*.22
                ang=aa+random.uniform(-.08,.08)+t*.12
                sub.append((r*cos(ang),r*sin(ang),5.245))
            curve_obj("split fold",sub,.045,teal)
# reconnecting labyrinth rings, much denser at rim
for ring in range(4,15):
    r=ring*.48
    count=8+ring*3
    pts=[]
    for i in range(count+1):
        a=2*pi*i/count + .08*sin(i*1.7+ring)
        rr=r + .10*sin(i*2.3+ring*.8)
        pts.append((rr*cos(a),rr*sin(a),5.255+.018*sin(i)))
    curve_obj("reconnected reaction maze",pts,.035 if ring>8 else .052,teal)
# glowing edge filaments
for k in range(46):
    a=2*pi*k/46
    r0=5.9+random.random()*1.15
    pts=[]
    for i in range(5):
        r=r0+i*.18
        aa=a+.05*sin(i*1.5+k)
        pts.append((r*cos(aa),r*sin(aa),5.29))
    curve_obj("edge luminous branch",pts,.032,teal)

# ---------- stem, roots and ground ----------
cyl("massive fibrous stem",(0,0,2.65),1.18,5.0,underside,20)
for k in range(13):
    a=2*pi*k/13
    pts=[(0.8*cos(a),0.8*sin(a),.15),(1.3*cos(a),1.2*sin(a),.42),(2.8*cos(a),2.0*sin(a),.18),(3.8*cos(a),2.8*sin(a),.08)]
    curve_obj("exposed mycelium root",pts,.11,mossmat)
cube("forest floor",(0,0,-.18),(28,28,.35),soilmat)
for k in range(18):
    a=random.random()*2*pi; r=random.uniform(1.2,5.8)
    cyl("moss mound",(r*cos(a),r*sin(a),.12),random.uniform(.18,.45),.25,mossmat,10)

# ---------- surrounding trunks ----------
trunks=[(-9,2,5.0,1.0), (9,1,5.8,1.15),(-7,-7,6.7,.85),(7,-7,7.4,.95),(-12,-3,8,.7),(12,-4,8,.75)]
for x,y,h,r in trunks:
    cyl("ancient habitat trunk",(x,y,h/2),r,h,bark,14)
    for j in range(4):
        a=random.random()*6.28
        curve_obj("gnarled branch",[(x,y,h*.45+j*.7),(x+random.uniform(-2,2),y+random.uniform(-2,2),h*.65+j*.65),(x+random.uniform(-3,3),y+random.uniform(-3,3),h*.8+j*.5)],.16,bark2)

# ---------- treehouses ----------
houses=[(-8.0,1.4,5.6,.3), (8.0,.8,6.3,-.35),(-6.8,-6.0,6.0,.2),(7.0,-6.3,7.0,-.2)]
for idx,(x,y,z,rot) in enumerate(houses):
    # platform and cabin
    cube("elevated lookout platform",(x,y,z-.35),(3.1,2.5,.28),wood,.12)
    body=cube("small treehouse cabin",(x,y,z+.65),(2.15,1.65,1.65),wood,.12)
    body.rotation_euler[2]=rot
    # roof
    bpy.ops.mesh.primitive_cone_add(vertices=4, radius1=1.65, radius2=0.18, depth=1.05, location=(x,y,z+1.95), rotation=(0,0,pi/4+rot))
    roof=bpy.context.object; roof.name="rust colored treehouse roof"; assign(roof,roofmat)
    # glowing front windows
    cube("golden window",(x-0.5,y-0.84,z+.72),(.48,.035,.55),windowmat,.04)
    cube("golden window",(x+0.5,y-0.84,z+.72),(.48,.035,.55),windowmat,.04)
    cyl("ladder post",(x-1.0,y-.9,z-1.3),.06,2.7,rope_mat,8)
    for j in range(6):
        cube("ladder rung",(x-1.0,y-.9,z-2.5+j*.42),(.65,.10,.08),wood,.02)

# bridges as paired ropes
def bridge(a,b,z):
    a=Vector(a); b=Vector(b)
    for side in (-.18,.18):
        pts=[]
        for i in range(9):
            t=i/8; p=a.lerp(b,t); p.z=z-0.65*sin(pi*t)+.15
            # lateral offset
            d=(b-a).normalized(); p += Vector((-d.y,d.x,0))*side
            pts.append(tuple(p))
        curve_obj("sagging rope bridge",pts,.045,rope_mat)
    for i in range(1,8):
        t=i/8; p=a.lerp(b,t); p.z=z-.65*sin(pi*t)+.15
        cube("bridge plank",tuple(p),(0.65,.42,.10),wood,.025)

for a,b,z in [((-7.5,1.2),(7.4,.8),5.45),((-6.7,-5.8),(0,-3.2),5.5),((7,-6),(0,3.2),6.0)]:
    bridge(a,b,z)

# ---------- foliage and ferns ----------
for k in range(48):
    a=random.random()*2*pi; r=random.uniform(6,13)
    x,y=r*cos(a),r*sin(a)
    z=random.uniform(.15,2.0)
    # fern stem and leaf blades
    curve_obj("fern frond",[(x,y,z),(x+.2*cos(a),y+.2*sin(a),z+.8),(x+.35*cos(a),y+.35*sin(a),z+1.45)],.025,leafmat)
    for j in range(4):
        p=(x+.25*cos(a)*j/3,y+.25*sin(a)*j/3,z+.35*j)
        cube("fern leaflet",(p[0]+.18*cos(a+1.4),p[1]+.18*sin(a+1.4),p[2]),(.38,.08,.06),leafmat,.02)

# ---------- spores ----------
for k in range(100):
    a=random.random()*2*pi; r=random.uniform(5.8,10.5)
    z=random.uniform(.3,5.0)*(1-r/13)+.35
    bpy.ops.mesh.primitive_ico_sphere_add(subdivisions=1,radius=random.uniform(.025,.075),location=(r*cos(a),r*sin(a),z))
    assign(bpy.context.object,sporemat)

# ---------- teal and warm lighting ----------
def light(name, typ, loc, color, energy, size=3):
    data=bpy.data.lights.new(name,typ); data.energy=energy; data.color=color
    if typ=='AREA': data.shape='DISK'; data.size=size
    ob=bpy.data.objects.new(name,data); bpy.context.collection.objects.link(ob); ob.location=loc
    return ob
light("cool undercap raking light",'AREA',(0,-2,4.0),(0.02,.65,.62),850,7).rotation_euler=(0.4,0,0)
light("warm sunset",'AREA',(-8,-10,10),(1.0,.16,.035),1400,8).rotation_euler=(0.55,0,-.55)
for x,y,z in [(-8,0,6),(8,0,7),(-7,-6,6),(7,-6,7)]:
    light("window amber",'POINT',(x,y-.9,z+.7),(1.0,.12,.025),90,1)

# ---------- camera ----------
bpy.ops.object.camera_add(location=(13,-17,8.0))
cam=bpy.context.object; cam.name="Low three-quarter underside presentation camera"
scene.camera=cam
cam.data.lens= forty = 42
target=Vector((0,0,4.9))
cam.rotation_euler=(target-Vector(cam.location)).to_track_quat('-Z','Y').to_euler()
cam.data.dof.use_dof=False

# compositor subtle glare for presentation, while geometry remains physically readable
scene.use_nodes=True
nt=scene.node_tree; nt.nodes.clear()
rl=nt.nodes.new("CompositorNodeRLayers")
gl=nt.nodes.new("CompositorNodeGlare"); gl.glare_type='FOG_GLOW'; gl.quality='HIGH'; gl.threshold=1.0; gl.size=7
comp=nt.nodes.new("CompositorNodeComposite")
nt.links.new(rl.outputs["Image"],gl.inputs["Image"]); nt.links.new(gl.outputs["Image"],comp.inputs["Image"])

# organize
for ob in bpy.context.scene.objects:
    if ob.type=='CURVE' and "reaction" in ob.name.lower():
        ob["procedural_origin"]="iterative Gray-Scott reaction diffusion field"
    ob.select_set(False)

bpy.ops.wm.save_as_mainfile(filepath="/tmp/colossal_reaction_diffusion_mushroom.blend")
```