diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml index e64cc07..9d9ee80 100644 --- a/.github/workflows/test.yml +++ b/.github/workflows/test.yml @@ -7,11 +7,17 @@ on: - main jobs: + test: - runs-on: ubuntu-latest + strategy: matrix: + os: ['ubuntu-latest', 'macos-latest', 'windows-latest'] python-version: ["3.10"] + runs-on: ${{ matrix.os }} + defaults: + run: + shell: bash -el {0} steps: - uses: actions/checkout@v3 @@ -20,11 +26,10 @@ jobs: with: python-version: ${{ matrix.python-version }} - name: Install dependencies - uses: mamba-org/setup-micromamba@v1 + uses: mamba-org/setup-micromamba@v1.4.3 with: - generate-run-shell: true environment-file: binder/environment.yml + init-shell: bash powershell - name: Test with pytest run: | - pytest --disable-pytest-warnings - shell: micromamba-shell {0} \ No newline at end of file + pytest --disable-pytest-warnings \ No newline at end of file diff --git a/binder/environment.yml b/binder/environment.yml index 7bfd5ba..a825317 100644 --- a/binder/environment.yml +++ b/binder/environment.yml @@ -10,9 +10,6 @@ dependencies: - toml - matplotlib - scipy - - boost - - boost-cpp - - cppyy=2.4 - openalea.lpy - openalea.plantgl=3.20.1 - openalea.mtg diff --git a/notebooks/simple_simulation.ipynb b/notebooks/simple_simulation.ipynb index 909b3d2..00c74d2 100644 --- a/notebooks/simple_simulation.ipynb +++ b/notebooks/simple_simulation.ipynb @@ -2,7 +2,7 @@ "cells": [ { "cell_type": "code", - "execution_count": null, + "execution_count": 1, "id": "c53a8105-14b5-437b-91ce-604152bb9c61", "metadata": { "tags": [] @@ -17,12 +17,28 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 2, "id": "b74cb80a-cac4-4d6f-ae18-165a52517072", "metadata": { "tags": [] }, - "outputs": [], + "outputs": [ + { + "data": { + "application/vnd.jupyter.widget-view+json": { + "model_id": "47503da9c36d4478af4b36964c741add", + "version_major": 2, + "version_minor": 0 + }, + "text/plain": [ + "SceneWidget(scenes=[{'id': 'XaKImpeVjzSwKnKEkcfVxz8VN', 'data': b'x\\xdaSLrw\\xf5\\xf7e`Pp\\xe0\\xe5RPVVd\\x00\\x020\\…" + ] + }, + "execution_count": 2, + "metadata": {}, + "output_type": "execute_result" + } + ], "source": [ "widget = SceneWidget()\n", "widget" @@ -30,7 +46,7 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 3, "id": "f4e4ddad-ba3d-40c8-9ecc-28740c010452", "metadata": { "tags": [] @@ -44,7 +60,7 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 4, "id": "89991d0d-2539-478d-8bb5-fe63a9a8095a", "metadata": { "tags": [] @@ -57,7 +73,7 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 5, "id": "8f50dd49", "metadata": { "tags": [] @@ -70,7 +86,7 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 6, "id": "b243f0f4-cf7c-4c35-8fa3-a6e894c36e8b", "metadata": { "tags": [] @@ -84,12 +100,53 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 7, "id": "64f4118d-93ff-4d43-99ed-f919d820003c", "metadata": { "tags": [] }, - "outputs": [], + "outputs": [ + { + "data": { + "text/plain": [ + "['edge_type',\n", + " 'label',\n", + " '_axial_id',\n", + " 'growth_units',\n", + " 'year',\n", + " 'inflorescence',\n", + " 'number',\n", + " 'closest_apex',\n", + " 'farthest_apex',\n", + " 'sons_nb',\n", + " 'observation',\n", + " 'parent_observation',\n", + " 'parent_unit_id',\n", + " 'parent_fbr_id',\n", + " 'parent_tree_id',\n", + " 'zone',\n", + " 'cumulated_mass',\n", + " 'radius',\n", + " 'offset',\n", + " 'developped',\n", + " 'phyllotactic_angle',\n", + " 'branching_angle',\n", + " 'rigidity',\n", + " 'age',\n", + " 'length',\n", + " 'trunk',\n", + " 'fruit_age',\n", + " 'fruit_mass',\n", + " 'leaf_age',\n", + " 'leaf_mass',\n", + " 'leaf_area']" + ] + }, + "execution_count": 7, + "metadata": {}, + "output_type": "execute_result" + } + ], "source": [ "tree.property_names()" ] diff --git a/test/notebooks/README.md b/test/notebooks/README.md deleted file mode 100644 index 6b405e9..0000000 --- a/test/notebooks/README.md +++ /dev/null @@ -1 +0,0 @@ -Some tests to verify that loading native libraries and C/C++ optimizations via [cppyy](https://github.com/wlav/cppyy) actually works. Simplifies the model setup & install. diff --git a/test/notebooks/optimization.ipynb b/test/notebooks/optimization.ipynb deleted file mode 100644 index 95894c4..0000000 --- a/test/notebooks/optimization.ipynb +++ /dev/null @@ -1,204 +0,0 @@ -{ - "cells": [ - { - "cell_type": "code", - "execution_count": 1, - "id": "3032b678-2cd6-4657-bc32-486b187140ef", - "metadata": {}, - "outputs": [], - "source": [ - "import os\n", - "import cppyy" - ] - }, - { - "cell_type": "code", - "execution_count": 2, - "id": "3205ce27-8792-45d5-9e37-3bdbd703b314", - "metadata": {}, - "outputs": [ - { - "data": { - "text/plain": [ - "True" - ] - }, - "execution_count": 2, - "metadata": {}, - "output_type": "execute_result" - } - ], - "source": [ - "cppyy.add_include_path(f'{os.environ[\"CONDA_PREFIX\"]}/include')\n", - "cppyy.add_library_path(f'{os.environ[\"CONDA_PREFIX\"]}/lib')\n", - "cppyy.load_library('pglmath')\n", - "cppyy.include('../../vmapplet/optimization.h')" - ] - }, - { - "cell_type": "code", - "execution_count": 3, - "id": "dce27e54-0a67-434b-9b1c-c8bdbf6355a8", - "metadata": {}, - "outputs": [], - "source": [ - "c_second_moment_of_area_annular_section = cppyy.gbl.optimization.second_moment_of_area_annular_section\n", - "c_rotate = cppyy.gbl.optimization.rotate" - ] - }, - { - "cell_type": "code", - "execution_count": 4, - "id": "d0812580-ce55-4579-a8dd-c6808167ae4d", - "metadata": {}, - "outputs": [], - "source": [ - "from math import sin, cos\n", - "from openalea.plantgl.all import Vector3\n", - "def py_second_moment_of_area_annular_section(inner_radius, thickness, section):\n", - " rt = inner_radius+thickness\n", - " rt2 = rt*rt\n", - " rt4 = rt2*rt2\n", - " r = inner_radius\n", - " r2 = r*r\n", - " r4=r2*r2\n", - " return 0.125 * (rt4 - r4)*(section + sin(section))\n", - "\n", - "def py_rotate(v3x, v3y, v3z, angle, vx, vy, vz):\n", - " c = cos(angle)\n", - " t2 = 1 - c\n", - " t6 = t2*v3x\n", - " t7 = t6*v3y\n", - " s = sin(angle)\n", - " t9 = s*v3z\n", - " t11 = t6*v3z\n", - " t12 = s*v3y\n", - " t19 = t2*v3y*v3z\n", - " t20 = s*v3x\n", - " t24 = v3z*v3z\n", - " R00 = c + t2*v3x*v3x\n", - " R01 = t7 - t9\n", - " R02 = t11 + t12\n", - " R10 = t7 + t9\n", - " R11 = c + t2*v3y*v3y\n", - " R12 = t19 - t20\n", - " R20 = t11 - t12\n", - " R21 = t19 + t20\n", - " R22 = c + t2*t24\n", - " return Vector3(R00*vx+R01*vy+R02*vz, R10*vx+R11*vy+R12*vz, R20*vx+R21*vy+R22*vz)" - ] - }, - { - "cell_type": "code", - "execution_count": 5, - "id": "1ee2930d-5e35-425f-ab5c-e6a99298b69c", - "metadata": {}, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "223 ns ± 6.79 ns per loop (mean ± std. dev. of 7 runs, 1,000,000 loops each)\n" - ] - } - ], - "source": [ - "%%timeit\n", - "py_second_moment_of_area_annular_section(1,2,3)" - ] - }, - { - "cell_type": "code", - "execution_count": 6, - "id": "35a1f638-235e-4171-8a97-7a68718459c5", - "metadata": {}, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "98.4 ns ± 1.09 ns per loop (mean ± std. dev. of 7 runs, 10,000,000 loops each)\n" - ] - } - ], - "source": [ - "%%timeit\n", - "c_second_moment_of_area_annular_section(1,2,3)" - ] - }, - { - "cell_type": "code", - "execution_count": 7, - "id": "e7f63991-27a2-4797-9d31-e8a98d6ad1de", - "metadata": {}, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "1.61 µs ± 8.89 ns per loop (mean ± std. dev. of 7 runs, 1,000,000 loops each)\n" - ] - } - ], - "source": [ - "%%timeit\n", - "py_rotate(1,2,3,4,5,6,7)" - ] - }, - { - "cell_type": "code", - "execution_count": 8, - "id": "47d9534a-dc1a-46d2-a9e4-899ee994f095", - "metadata": {}, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "622 ns ± 5.99 ns per loop (mean ± std. dev. of 7 runs, 1,000,000 loops each)\n" - ] - } - ], - "source": [ - "%%timeit\n", - "c_rotate(1,2,3,4,5,6,7)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "c87a4529-2b28-4645-9500-d79f2b1102b2", - "metadata": {}, - "outputs": [], - "source": [] - } - ], - "metadata": { - "kernelspec": { - "display_name": "Python 3 (ipykernel)", - "language": "python", - "name": "python3" - }, - "language_info": { - "codemirror_mode": { - "name": "ipython", - "version": 3 - }, - "file_extension": ".py", - "mimetype": "text/x-python", - "name": "python", - "nbconvert_exporter": "python", - "pygments_lexer": "ipython3", - "version": "3.9.15" - }, - "widgets": { - "application/vnd.jupyter.widget-state+json": { - "state": {}, - "version_major": 2, - "version_minor": 0 - } - } - }, - "nbformat": 4, - "nbformat_minor": 5 -} diff --git a/test/test_metamer.py b/test/test_metamer.py index 5442172..cba2dc0 100644 --- a/test/test_metamer.py +++ b/test/test_metamer.py @@ -2,12 +2,11 @@ from vmapplet import Simulation, Options from vmapplet.organs.metamer import MetamerData, CambialLayer -from vmapplet.optimisation import reaction_wood_target from vmapplet.organs.wood import Wood from vmapplet.organs.fruit import AppleFruit from vmapplet.organs.leaf import AppleLeaf from vmapplet.organs.internode import Internode -from vmapplet.physics import Frame +from vmapplet.physics import Frame, reaction_wood_target from vmapplet.enums import FruitState diff --git a/test/test_physics.py b/test/test_physics.py index b0f1deb..8c320f8 100644 --- a/test/test_physics.py +++ b/test/test_physics.py @@ -8,7 +8,6 @@ second_moment_of_area_circular_section, Frame, rotate_frame_at_branch, - reorient_frame, ) @@ -52,13 +51,3 @@ def test_rotate_frame_at_branch(): v3 = Vector3(1, 1, 2) frame = Frame(v1, v2, v3) rotate_frame_at_branch(frame, 45, 45) - - -def test_reorient_frame(): - v1 = Vector3(2, 1, 1) - v2 = Vector3(1, 2, 1) - v3 = Vector3(1, 1, 2) - rotation_velocity = Vector3(1, 1, 1) - length = 1 - frame = Frame(v1, v2, v3) - reorient_frame(frame, rotation_velocity, length) diff --git a/vmapplet/__init__.py b/vmapplet/__init__.py index aa2eb89..7a403d7 100644 --- a/vmapplet/__init__.py +++ b/vmapplet/__init__.py @@ -1,13 +1,13 @@ -from typing import Optional +from typing import Optional, TYPE_CHECKING -from pgljupyter import SceneWidget +if TYPE_CHECKING: + from pgljupyter import SceneWidget -from ._cppyy import * # noqa from .simulation import Simulation from .options import Options -def run(simulation: Simulation, scene_widget: Optional[SceneWidget] = None): +def run(simulation: Simulation, scene_widget: Optional["SceneWidget"] = None): """ """ date_start = simulation.options.general.date_start diff --git a/vmapplet/_cppyy.py b/vmapplet/_cppyy.py deleted file mode 100644 index 08d299a..0000000 --- a/vmapplet/_cppyy.py +++ /dev/null @@ -1,10 +0,0 @@ -import os -import pathlib - -import cppyy - -# initalize cppyy -cppyy.add_include_path(f'{os.environ["CONDA_PREFIX"]}/include') -cppyy.add_library_path(f'{os.environ["CONDA_PREFIX"]}/lib') -cppyy.include(pathlib.Path(__file__).parent.resolve() / "optimization.h") -cppyy.load_library("pglmath") diff --git a/vmapplet/lpy/mechanics_backward.lpy b/vmapplet/lpy/mechanics_backward.lpy index 2e04923..f4fd29a 100644 --- a/vmapplet/lpy/mechanics_backward.lpy +++ b/vmapplet/lpy/mechanics_backward.lpy @@ -1,10 +1,10 @@ import openalea.plantgl.all as pgl from vmapplet import ( - constants, - optimisation + constants ) from vmapplet.enums import LeafState +from vmapplet.physics import * gravity = pgl.Vector3(0.0, 0.0, -9.81) @@ -28,10 +28,10 @@ root() >> metamer(m): tree.fruit_load = tree.fruits / tree.trunk_cross_sectional_area metamer(m) >> SB() branch() metamer(mb) EB() metamer(mr): - radius = optimisation.get_new_radius(mb.radius, mr.radius) + radius = get_new_radius(mb.radius, mr.radius) if m.leaf.state == LeafState.GROWING: - radius = optimisation.get_new_radius(radius, m.leaf.petiole_radius) - m.radius = optimisation.max(radius, m.radius); + radius = get_new_radius(radius, m.leaf.petiole_radius) + m.radius = max(radius, m.radius); # update last layer thickness m.layers[-1].thickness = m.radius - m.layers[-1].radius m.compute_mass(mr, mb) @@ -49,8 +49,8 @@ metamer(m) >> SB() branch() metamer(mb) EB() metamer(mr): metamer(m) >> metamer(mr): radius = mr.radius if m.leaf.state == LeafState.GROWING: - radius = optimisation.get_new_radius(mr.radius, m.leaf.petiole_radius) - m.radius = optimisation.max(radius, m.radius) + radius = get_new_radius(mr.radius, m.leaf.petiole_radius) + m.radius = max(radius, m.radius) m.layers[-1].thickness = m.radius - m.layers[-1].radius m.compute_mass(mr) if options.general.mechanics: @@ -66,8 +66,8 @@ metamer(m) >> apex(a): # wood.density, m.fruit_mass are units objects radius = a.radius if m.leaf.state == LeafState.GROWING: - radius = optimisation.get_new_radius(a.radius, m.leaf.petiole_radius) - m.radius = optimisation.max(radius, m.radius) + radius = get_new_radius(a.radius, m.leaf.petiole_radius) + m.radius = max(radius, m.radius) m.layers[-1].thickness = m.radius - m.layers[-1].radius m.compute_mass() m.cumulated_torque = cross( m.hlu.heading * m.length, tree.tropism) diff --git a/vmapplet/lpy/mechanics_forward.lpy b/vmapplet/lpy/mechanics_forward.lpy index dc62a8d..3f0f54f 100644 --- a/vmapplet/lpy/mechanics_forward.lpy +++ b/vmapplet/lpy/mechanics_forward.lpy @@ -1,5 +1,4 @@ -from vmapplet import optimisation -from vmapplet.physics import rotate_frame_at_branch +from vmapplet.physics import * module apex(apex_data): scale=2 module branch(): scale=1 @@ -14,12 +13,12 @@ production: metamer(ml) branch() << metamer(m): m.hlu = rotate_frame_at_branch(ml.hlu, ml.branching_angle, ml.phyllotactic_angle); - m.hlu = optimisation.reorient_frame(m.hlu, m.rotation_velocity, m.rv_norm, m.length) + m.hlu = reorient_frame(m.hlu, m.rotation_velocity, m.rv_norm, m.length) m.update_position(ml.position) produce metamer(m) metamer(ml) << metamer(m): - m.hlu = optimisation.reorient_frame(ml.hlu, m.rotation_velocity, m.rv_norm, m.length) + m.hlu = reorient_frame(ml.hlu, m.rotation_velocity, m.rv_norm, m.length) m.update_position(ml.position) produce metamer(m) diff --git a/vmapplet/lpy/parameters.lpy b/vmapplet/lpy/parameters.lpy index 046f1d5..dacd68b 100644 --- a/vmapplet/lpy/parameters.lpy +++ b/vmapplet/lpy/parameters.lpy @@ -1,6 +1,3 @@ -from vmapplet import optimisation -from vmapplet.physics import rotate_frame_at_branch - module apex(apex_data): scale=2 module branch(): scale=1 module growth_unit(growth_unit_data): scale=1 diff --git a/vmapplet/optimisation.py b/vmapplet/optimisation.py deleted file mode 100644 index d32d292..0000000 --- a/vmapplet/optimisation.py +++ /dev/null @@ -1,101 +0,0 @@ -from math import fabs, acos - -import cppyy - -from .frame import Frame - -second_moment_of_area_annular_section = ( - cppyy.gbl.optimization.second_moment_of_area_annular_section -) -second_moment_of_area_circle = cppyy.gbl.optimization.second_moment_of_area_circle -get_new_radius = cppyy.gbl.optimization.get_new_radius -# _reaction_wood_target = cppyy.gbl.optimization.reaction_wood_target - - -def rotate( - v3x: float, v3y: float, v3z: float, angle: float, vx: float, vy: float, vz: float -): - # pgl can not handle the vector: wrap in list - return list(cppyy.gbl.optimization.rotate(v3x, v3y, v3z, angle, vx, vy, vz)) - - -import openalea.plantgl.all as pgl - - -def reorient_frame(initial_hlu, rotation_velocity, rv_norm, length): - h = pgl.Vector3(initial_hlu.heading) - h.normalize() - l = pgl.Vector3(initial_hlu.left) - l.normalize() - - # vl = rotation_velocity.normalize() #_ look at v3d length definition - # vl is replaced by rv_norm - - if fabs(rv_norm * length) >= 0.01: - h = pgl.Vector3( - rotate( - rotation_velocity.x, - rotation_velocity.y, - rotation_velocity.z, - rv_norm * length, - h.x, - h.y, - h.z, - ) - ) - l = pgl.Vector3( - rotate( - rotation_velocity.x, - rotation_velocity.y, - rotation_velocity.z, - rv_norm * length, - l.x, - l.y, - l.z, - ) - ) - h.normalize() - l.normalize() - return Frame(h, l, pgl.cross(h, l)) - - -# cppyy complains maybe clash with boost PyObjects and plantgl types -# def reaction_wood_target(a, b, c): -# # print('reaction_wood_target') -# return _reaction_wood_target(list(a), list(b), list(c)) - - -def reaction_wood_target(up, heading, previous_heading): - cos_gh = pgl.Vector3(0.0, 0.0, 1.0) * heading - cos_pu = previous_heading * up - cos_ph = previous_heading * heading - - inclination = 0 - if cos_pu * cos_ph >= 0.0: - try: - inclination = acos(cos_ph) - except Exception: - pass - # print(str(cos_ph)) - else: - try: - inclination = -acos(cos_ph) - except Exception: - pass - # print(str(cos_ph)) - percentage = 0.1635 * (1.0 - cos_gh) - 0.1778 * inclination - r = 3.14159 * 2.0 * percentage - - if r < 0.0: - r = 0.0 - elif r > 3.14159: - r = 3.141459 - - return r - - -def max(a, b): - if a > b: - return a - else: - return b diff --git a/vmapplet/optimization.h b/vmapplet/optimization.h deleted file mode 100644 index 41479a7..0000000 --- a/vmapplet/optimization.h +++ /dev/null @@ -1,86 +0,0 @@ -#include -#include -#include "plantgl/math/util_vector.h" - -namespace optimization { - -typedef std::array vec3; - - -vec3 rotate(float v3x, float v3y, float v3z, float angle, float vx, float vy, float vz) -{ - float c = cosf(angle); - float t2 = 1 - c; - float t6 = t2*v3x; - float t7 = t6*v3y; - float s = sinf(angle); - float t9 = s*v3z; - float t11 = t6*v3z; - float t12 = s*v3y; - float t19 = t2*v3y*v3z; - float t20 = s*v3x; - float t24 = v3z*v3z; - float R00 = c + t2*v3x*v3x; - float R01 = t7 - t9; - float R02 = t11 + t12; - float R10 = t7 + t9; - float R11 = c + t2*v3y*v3y; - float R12 = t19 - t20; - float R20 = t11 - t12; - float R21 = t19 + t20; - float R22 = c + t2*t24; - return {R00*vx+R01*vy+R02*vz, R10*vx+R11*vy+R12*vz, R20*vx+R21*vy+R22*vz}; -} - -double second_moment_of_area_annular_section(double inner_radius, double thickness, double section) -{ - double rt = inner_radius + thickness; - double rt2 = rt * rt; - double rt4 = rt2 * rt2; - double r = inner_radius; - double r2 = r * r; - double r4= r2 * r2; - return 0.125 * (rt4 - r4) * (section + sinf(section)); -} - -double second_moment_of_area_circle(double radius) -{ - return 0.78539816339744828 * radius * radius * radius * radius; -} - -double get_new_radius(double ra,double rb, double exponent = 2.49, double previous_rt = -1.) -{ - double rap = powf(ra, exponent); - double rbp = powf(rb, exponent); - double newrt = powf(rap+rbp, 1./exponent); - return newrt; -} - -// cppyy runns but complains: IncrementalExecutor::executeFunction: symbol '_ZN3PGLmlERKNS_7Vector3ES2_' unresolved while linking symbol '__cf_11'! -// double reaction_wood_target(vec3 _up, vec3 _heading, vec3 _previous_heading) -// { -// PGL::Vector3 up = PGL::Vector3(_up.at(0), _up.at(1), _up.at(2)); -// PGL::Vector3 heading = PGL::Vector3(_heading.at(0), _heading.at(1), _heading.at(2)); -// PGL::Vector3 previous_heading = PGL::Vector3(_previous_heading.at(0), _previous_heading.at(1), _previous_heading.at(2)); -// double cos_gh = PGL::Vector3(0.0, 0.0, 1.0) * heading; -// double cos_pu = previous_heading * up; -// double cos_ph = previous_heading * heading; -// double inclination, percentage, r; - -// if (cos_pu*cos_ph >= 0.0) -// inclination = acosf(cos_ph); -// else -// inclination = -acosf(cos_ph); - -// percentage = 0.1635 * (1.0 - cos_gh) - 0.1778 * inclination; -// r = 3.14159*2. * percentage; - -// if (r < 0.0) -// r = 0.0; -// else if (r > 3.14159) -// r = 3.141459; - -// return r; -// } - -} diff --git a/vmapplet/organs/metamer.py b/vmapplet/organs/metamer.py index f2e2820..00cd6b6 100644 --- a/vmapplet/organs/metamer.py +++ b/vmapplet/organs/metamer.py @@ -4,8 +4,13 @@ from openalea.plantgl.all import Vector3, dot -from .. import constants, optimisation -from ..physics import rotate_frame_at_branch +from .. import constants +from ..physics import ( + rotate_frame_at_branch, + second_moment_of_area_annular_section, + reaction_wood_target, + second_moment_of_area_circle, +) from ..srandom import boolean_event from ..enums import Zone, FruitState, LeafState @@ -312,7 +317,7 @@ def update_metamer_parameters(self, simulation, cambial=None): # if we create a new layer, then computation on previous one will # be redundant. Compute them once for all second_moment_of_area = ( - optimisation.second_moment_of_area_annular_section( + second_moment_of_area_annular_section( self.layers[-1].radius, self.layers[-1].thickness, self.layers[-1].reaction_wood, @@ -329,7 +334,7 @@ def update_metamer_parameters(self, simulation, cambial=None): # TEST if not the central layer # TODO days or seconds ? if self.nlayers >= 2: - r = optimisation.reaction_wood_target( + r = reaction_wood_target( self.hlu.up, self.hlu.heading, self.season_initial_heading ) if r > self.layers[-1].reaction_wood: @@ -347,11 +352,10 @@ def update_metamer_parameters(self, simulation, cambial=None): # Updating second_moment_of_area second_moment_of_area = ( - self.total_second_moment_of_area - + optimisation.second_moment_of_area_circle(self.radius) + self.total_second_moment_of_area + second_moment_of_area_circle(self.radius) ) second_moment_of_area += ( - optimisation.second_moment_of_area_annular_section( + second_moment_of_area_annular_section( self.layers[-1].radius, self.layers[-1].thickness, self.layers[-1].reaction_wood, @@ -514,67 +518,6 @@ def leaf_mass(self): return self.leaf.mass -def reaction_wood_target(up, heading, previous_heading): - r"""Reaction wood target - - The reaction wood is proportional to the change in inclination - over a season (Almeras, 2001). Hypothesis: The reaction wood is - also proportional to the inclination from gravity. - - .. math:: - - P_r = 0.1635 - 0.1778 \theta_s - - where :math:`P_r` is the radial portion of the outermost cambial layer that - become reaction wood and :math:`\theta_s` is the change in inclination of the - internode over the season (i.e., the cahnge in :math:`\vec{H}`, the heading - vector of the HLU frame since the start of the spring. - - In order to also take into account the gravity, The firs term in the equation above - must be multiply by a coefficient that varies with the angle between the internode and - :math:`\vec{u}` a unit vector opposite to :math:`-\vec{g}` - - :param Vector3 up: - :param heading: vector3 - :param previous_heading: vector3 - :returns: reaction_wood (double) - """ - - # multiply by gravity normalised - cos_gh = Vector3(0.0, 0.0, 1.0) * heading - cos_pu = previous_heading * up - cos_ph = previous_heading * heading - - if cos_pu * cos_ph >= 0.0: - try: - inclination = acos(cos_ph / 1.0001) - except Exception: - print("Problem with acos(cos_ph) where cos_ph=%f" % cos_ph) - inclination = 0.0 - else: - try: - inclination = -acos(cos_ph) - except Exception: - tol = 1e-6 - try: - inclination = -acos(cos_ph - tol) - ValueError("try again with tol set %s %s" % (cos_ph, cos_ph - 1.0)) - except Exception: - print("cos+_pu=", cos_pu, " cos_ph=", cos_ph, "cosgh=", cos_gh) - print(cos_ph - 1.0) - raise ValueError("Problem with acos(cos_ph) where cos_ph=%f" % cos_ph) - - percentage = 0.1635 * (1.0 - cos_gh) - 0.1778 * inclination - r = constants.two_pi * percentage - - if r < 0.0: - r = 0.0 - elif r > constants.pi: - r = constants.pi - - return r - - def _clamp_if_near_zero(data, tol=0.0001): if data < tol and data > -tol: return 0 diff --git a/vmapplet/physics.py b/vmapplet/physics.py index 81b22a0..6c70e4a 100644 --- a/vmapplet/physics.py +++ b/vmapplet/physics.py @@ -19,22 +19,114 @@ from openalea.stocatree.physics import * -.. note:: many functions have been reimplemented in optimisation.pyx - """ -from math import sin +from math import acos, cos, sin, pow, fabs import openalea.plantgl.all as pgl from . import constants -from . import optimisation from .frame import Frame error_tolerance = 0.0001 -def reorient_frame(initial_hlu, rotation_velocity, length): +def reaction_wood_target(up, heading, previous_heading): + r"""Reaction wood target + + The reaction wood is proportional to the change in inclination + over a season (Almeras, 2001). Hypothesis: The reaction wood is + also proportional to the inclination from gravity. + + .. math:: + + P_r = 0.1635 - 0.1778 \theta_s + + where :math:`P_r` is the radial portion of the outermost cambial layer that + become reaction wood and :math:`\theta_s` is the change in inclination of the + internode over the season (i.e., the cahnge in :math:`\vec{H}`, the heading + vector of the HLU frame since the start of the spring. + + In order to also take into account the gravity, The firs term in the equation above + must be multiply by a coefficient that varies with the angle between the internode and + :math:`\vec{u}` a unit vector opposite to :math:`-\vec{g}` + + :param Vector3 up: + :param heading: vector3 + :param previous_heading: vector3 + :returns: reaction_wood (double) + """ + + # multiply by gravity normalised + cos_gh = pgl.Vector3(0.0, 0.0, 1.0) * heading + cos_pu = previous_heading * up + cos_ph = previous_heading * heading + inclination = 0.0 + + if cos_pu * cos_ph >= 0.0: + try: + inclination = acos(cos_ph / 1.0001) + except Exception: + print("Problem with acos(cos_ph) where cos_ph=%f" % cos_ph) + inclination = 0.0 + else: + try: + inclination = -acos(cos_ph) + except Exception: + tol = 1e-6 + try: + inclination = -acos(cos_ph - tol) + ValueError("try again with tol set %s %s" % (cos_ph, cos_ph - 1.0)) + except Exception: + print("cos+_pu=", cos_pu, " cos_ph=", cos_ph, "cosgh=", cos_gh) + print(cos_ph - 1.0) + raise ValueError("Problem with acos(cos_ph) where cos_ph=%f" % cos_ph) + + percentage = 0.1635 * (1.0 - cos_gh) - 0.1778 * inclination + r = constants.two_pi * percentage + + if r < 0.0: + r = 0.0 + elif r > constants.pi: + r = constants.pi + + return r + + +def rotate( + v3x: float, v3y: float, v3z: float, angle: float, vx: float, vy: float, vz: float +) -> pgl.Vector3: + c = cos(angle) + t2 = 1 - c + t6 = t2 * v3x + t7 = t6 * v3y + s = sin(angle) + t9 = s * v3z + t11 = t6 * v3z + t12 = s * v3y + t19 = t2 * v3y * v3z + t20 = s * v3x + t24 = v3z * v3z + R00 = c + t2 * v3x * v3x + R01 = t7 - t9 + R02 = t11 + t12 + R10 = t7 + t9 + R11 = c + t2 * v3y * v3y + R12 = t19 - t20 + R20 = t11 - t12 + R21 = t19 + t20 + R22 = c + t2 * t24 + + return pgl.Vector3( + [ + R00 * vx + R01 * vy + R02 * vz, + R10 * vx + R11 * vy + R12 * vz, + R20 * vx + R21 * vy + R22 * vz, + ] + ) + + +def reorient_frame(initial_hlu, rotation_velocity, rv_norm, length): """Reorient frame The initial frame is an HLU (Heading Up Left, turtle axes) frame, made of tree orthogonal vectors that indicate the heading, @@ -53,50 +145,39 @@ def reorient_frame(initial_hlu, rotation_velocity, length): .. todo:: describe the algorithm - .. note:: this function has been optimised to used the rotate function from - optimisation.pyx that computes the rotation of the AxisAngle around a - given vector. It replaces the AxisAngle./quaternion of the origninal code of - MappleT. - - :: - - import openalea.stocatree.optimisation as optimisation - optimisation.rotate(hlu, pgl.Vector3, length) - """ - heading = pgl.Vector3(initial_hlu.heading) - heading.normalize() - left = pgl.Vector3(initial_hlu.left) - left.normalize() - velocity = rotation_velocity.normalize() - if abs(velocity * length) >= 0.01: - heading = pgl.Vector3( - optimisation.rotate( - rotation_velocity.x, - rotation_velocity.y, - rotation_velocity.z, - velocity * length, - heading.x, - heading.y, - heading.z, - ) + h = pgl.Vector3(initial_hlu.heading) + h.normalize() + l = pgl.Vector3(initial_hlu.left) + l.normalize() + + # vl = rotation_velocity.normalize() #_ look at v3d length definition + # vl is replaced by rv_norm + + if fabs(rv_norm * length) >= 0.01: + h = rotate( + rotation_velocity.x, + rotation_velocity.y, + rotation_velocity.z, + rv_norm * length, + h.x, + h.y, + h.z, ) - left = pgl.Vector3( - optimisation.rotate( - rotation_velocity.x, - rotation_velocity.y, - rotation_velocity.z, - velocity * length, - left.x, - left.y, - left.z, - ) + l = rotate( + rotation_velocity.x, + rotation_velocity.y, + rotation_velocity.z, + rv_norm * length, + l.x, + l.y, + l.z, ) - heading.normalize() - left.normalize() - return Frame(heading, left, pgl.cross(heading, left)) + h.normalize() + l.normalize() + return Frame(h, l, pgl.cross(h, l)) def second_moment_of_area_circle(radius): @@ -133,11 +214,6 @@ def second_moment_of_area_circle(radius): :References: MappleT - .. note:: this function is duplicated in a cython version inside optimisation.pyx - - >>> import openalea.stocatree.optimisation as optimisation - >>> a = optimisation.second_moment_of_area_circle(1.) - """ @@ -181,10 +257,6 @@ def second_moment_of_area_annular_section(inner_radius, thickness, section): :param thickness: the thickness in meters :param section: the section angle in radians - .. warning:: this function is duplicated in a cython version inside optimisation.pyx - - >>> import openalea.stocatree.optimisation as optimisation - >>> a = optimisation.second_moment_of_area_annular_section(1., 0.1, 0.78) """ assert section <= 2 * constants.pi @@ -206,62 +278,54 @@ def rotate_frame_at_branch(initial_hlu, branching_angle, phyllotactic_angle): :param phyllotactic_angle: the angle related to phyllotaxy .. note:: this function has been optimised to used the rotate function from - optimisation.pyx that computes the rotation of the AxisAngle around a + pyx that computes the rotation of the AxisAngle around a given vector. It replaces the AxisAngle./quaternion of the origninal code of MappleT. - .. note:: this function is not part of optimisation.pyx since it already use the + .. note:: this function is not part of pyx since it already use the rotate function and is not called as much as the reorient_frame or rotate function itself. """ hlu = Frame(initial_hlu.heading, initial_hlu.up, initial_hlu.left) - hlu.heading = pgl.Vector3( - optimisation.rotate( - initial_hlu.left.x, - initial_hlu.left.y, - initial_hlu.left.z, - branching_angle, - initial_hlu.heading.x, - initial_hlu.heading.y, - initial_hlu.heading.z, - ) + hlu.heading = rotate( + initial_hlu.left.x, + initial_hlu.left.y, + initial_hlu.left.z, + branching_angle, + initial_hlu.heading.x, + initial_hlu.heading.y, + initial_hlu.heading.z, ) - hlu.up = pgl.Vector3( - optimisation.rotate( - initial_hlu.left.x, - initial_hlu.left.y, - initial_hlu.left.z, - branching_angle, - initial_hlu.up.x, - initial_hlu.up.y, - initial_hlu.up.z, - ) + hlu.up = rotate( + initial_hlu.left.x, + initial_hlu.left.y, + initial_hlu.left.z, + branching_angle, + initial_hlu.up.x, + initial_hlu.up.y, + initial_hlu.up.z, ) hlu.heading.normalize() hlu.up.normalize() - hlu.heading = pgl.Vector3( - optimisation.rotate( - initial_hlu.heading.x, - initial_hlu.heading.y, - initial_hlu.heading.z, - phyllotactic_angle, - hlu.heading.x, - hlu.heading.y, - hlu.heading.z, - ) + hlu.heading = rotate( + initial_hlu.heading.x, + initial_hlu.heading.y, + initial_hlu.heading.z, + phyllotactic_angle, + hlu.heading.x, + hlu.heading.y, + hlu.heading.z, ) - hlu.up = pgl.Vector3( - optimisation.rotate( - initial_hlu.heading.x, - initial_hlu.heading.y, - initial_hlu.heading.z, - phyllotactic_angle, - hlu.up.x, - hlu.up.y, - hlu.up.z, - ) + hlu.up = rotate( + initial_hlu.heading.x, + initial_hlu.heading.y, + initial_hlu.heading.z, + phyllotactic_angle, + hlu.up.x, + hlu.up.y, + hlu.up.z, ) hlu.heading.normalize() @@ -289,3 +353,12 @@ def rupture(torque, radius, modulus_of_rupture=50e6): .. todo: not used for the moment """ return stress(torque, radius) > modulus_of_rupture + + +def get_new_radius( + ra: float, rb: float, exponent: float = 2.49, previous_rt: float = -1.0 +): + rap = pow(ra, exponent) + rbp = pow(rb, exponent) + newrt = pow(rap + rbp, 1.0 / exponent) + return newrt