Changeset 667 for cpp/frams/genetics/fF/conv_fF.cpp

Timestamp:

06/07/17 12:01:16 (8 years ago)

Author:

Maciej Komosinski

Message:

• fF foraminifera encoding supports scaling (radius) of the initial chamber (three new parameters added: radii for x,y,z)
• mutation ignores these three parameters
• float -> double for consistency everywhere

File:

• cpp/frams/genetics/fF/conv_fF.cpp (modified) (12 diffs)

cpp/frams/genetics/fF/conv_fF.cpp

@@ -1 +1 @@
 // This file is a part of Framsticks SDK.  http://www.framsticks.com/
-// Copyright (C) 1999-2015  Maciej Komosinski and Szymon Ulatowski.
+// Copyright (C) 1999-2017  Maciej Komosinski and Szymon Ulatowski.
 // See LICENSE.txt for details.
 
@@ -10 +10 @@
 GenoConv_fF0::GenoConv_fF0()
 {
-        name = "7-value Foraminifera encoding";
+        name = "10-parameter Foraminifera encoding";
         in_format = 'F';
         out_format = '0';
@@ -42 +42 @@
                 return ""; //so we return an invalid f0 genotype
 
-        double div_radius_length = 1; //div_radius_length=1 or kx=ky=kz=1
-        double radius = 1;
-
         Model m;
         m.open();
@@ -52 +49 @@
         fF_chamber3d **chambers = new fF_chamber3d*[gp.number_of_chambers];
         for (int i = 0; i < gp.number_of_chambers; i++)
-                createSphere(i, chambers, radius, div_radius_length, gp.translation, gp.angle1, gp.angle2, gp.scalex, gp.scaley, gp.scalez);
+                createSphere(i, chambers, gp.radius0x, gp.radius0y, gp.radius0z, gp.translation, gp.angle1, gp.angle2, gp.scalex, gp.scaley, gp.scalez);
 
         Part *p1 = addNewPart(&m, chambers[0]);
+        p1->scale = Pt3D(gp.radius0x, gp.radius0y, gp.radius0z); //size of the initial chamber
         for (int i = 1; i < gp.number_of_chambers; i++)
         {
@@ -71 +69 @@
 }
 
-void GenoConv_fF0::createSphere(int which, fF_chamber3d **chambers, double radius_, double div_radius_length_, double div_vector_length_,
-        double alpha_, double gamma_, double kx_, double ky_, double kz_)
-{
-        chambers[which] = new fF_chamber3d(0.0f, 0.0f, 0.0f,
-                (float)radius_, (float)radius_ * (float)kx_, 0.0f, 0.0f,
-                (float)(radius_ * div_vector_length_), 0.0f, 0.0f, 0.0f, 0.0f);
+void GenoConv_fF0::createSphere(int which, fF_chamber3d **chambers, double radius0x, double radius0y, double radius0z, double translation, double alpha_, double gamma_, double kx_, double ky_, double kz_)
+{
+        chambers[which] = new fF_chamber3d(0.0, 0.0, 0.0,
+                radius0x, radius0y, radius0z, radius0x * kx_, 0.0, 0.0,
+                radius0x * translation, 0.0, 0.0, 0.0, 0.0);
         if (which == 0)
-                chambers[which]->points = generate_points(chambers[which], which, kx_, ky_, kz_);
-        if (which > 0) {
-                /* old radius */
-                double radiusOld, radius;
-                radiusOld = chambers[which - 1]->radius;
-                radius = div_radius_length_ * radiusOld;
+                chambers[which]->points = generate_points(chambers[which]);
+        if (which > 0)
+        {
+                chambers[which]->radius_x = get_radius(chambers[which - 1]->radius_x, kx_, chambers[0]->radius_x);
+                chambers[which]->radius_y = get_radius(chambers[which - 1]->radius_y, ky_, chambers[0]->radius_y);
+                chambers[which]->radius_z = get_radius(chambers[which - 1]->radius_z, kz_, chambers[0]->radius_z);
+
                 /* new growth vector length */
-                double len = radius * div_vector_length_;
-                if (radius < fF_TOO_LITTLE) {
-                        radius = fF_TOO_LITTLE;
-                        if (fabs(len) > (fF_TOO_MUCH * radius)) {
-                                len = ((len < 0) ? (-1) : 1) * fF_TOO_MUCH * radius;
-                        }
-                }
-                if (len == 0) {
+                double len = chambers[which]->radius_y * translation;
+                double max_radius = fF_TOO_MUCH * chambers[which]->radius_y;
+                if (fabs(len) > (max_radius))
+                        len = ((len < 0) ? (-1) : 1) * max_radius;
+                if (len == 0)
                         len = -0.0000001;
-                }
 
                 /* aperture of the previous chamber */
@@ -131 +125 @@
                 double alpha = angle - alpha_;
 
-
                 double gamma = chambers[which - 1]->phi + gamma_;
 
-                /* x */
                 double wx = len * cos(alpha);
-                /* y */
                 double wy = len * sin(alpha);
-                /* y */
                 double wz = len * sin(alpha) * sin(gamma);
 
@@ -146 +136 @@
                 double z = pzz + wz;
 
-                chambers[which]->centerX = (float)x;
-                chambers[which]->centerY = (float)y;
-                chambers[which]->centerZ = (float)z;
-                chambers[which]->radius = (float)radius;
-                chambers[which]->vectorTfX = (float)wx;
-                chambers[which]->vectorTfY = (float)wy;
-                chambers[which]->vectorTfZ = (float)wz;
-                chambers[which]->beta = (float)alpha;
-                chambers[which]->phi = (float)gamma;
-
-                chambers[which]->points = generate_points(chambers[which], which, kx_, ky_, kz_);
-                search_hid(which, chambers, kx_, ky_, kz_);
+                chambers[which]->centerX = x;
+                chambers[which]->centerY = y;
+                chambers[which]->centerZ = z;
+                chambers[which]->vectorTfX = wx;
+                chambers[which]->vectorTfY = wy;
+                chambers[which]->vectorTfZ = wz;
+                chambers[which]->beta = alpha;
+                chambers[which]->phi = gamma;
+
+                chambers[which]->points = generate_points(chambers[which]);
+                search_hid(which, chambers);
                 int pun;
-                pun = find_hole(which, pzx, pzy, pzz, chambers, kx_, ky_, kz_);
-
-                chambers[which]->holeX = (float)chambers[which]->points[pun].x;
-                chambers[which]->holeY = (float)chambers[which]->points[pun].y;
-                chambers[which]->holeZ = (float)chambers[which]->points[pun].z;
-        }
+                pun = find_hole(which, pzx, pzy, pzz, chambers);
+                if (pun < 0) //should never happen
+                {
+                        logPrintf("GenoConv_fF0", "createSphere", LOG_ERROR, "find_hole(%d) returned %d", which, pun);
+                        pun = 0;
+                }
+                chambers[which]->holeX = chambers[which]->points[pun].x;
+                chambers[which]->holeY = chambers[which]->points[pun].y;
+                chambers[which]->holeZ = chambers[which]->points[pun].z;
+        }
+}
+
+double GenoConv_fF0::get_radius(double prev_radius, double scale, double start_radius)
+{
+        double radius = prev_radius * scale;
+        double min_radius = fF_TOO_LITTLE*start_radius;
+        if (radius < min_radius) {
+                radius = min_radius;
+        }
+
+        return radius;
 }
 
 void GenoConv_fF0::precompute_cos_and_sin()
 {
-        int i;
-        double pi = acos(-1.0);
-        double angle = pi / (((double)fF_LATITUDE_NUM)*0.5);
-        for (i = 0; i < fF_LATITUDE_NUM; i++)
-        {
-                cosines[i] = cos((double)i * angle);
-                sines[i] = sin((double)i * angle);
-        }
-}
-
-fF_point* GenoConv_fF0::generate_points(fF_chamber3d *chamber, int which, double kx_, double ky_, double kz_)
-{
-        float radius = chamber->radius;
-        float cenx = chamber->centerX;
-        float ceny = chamber->centerY;
-        float cenz = chamber->centerZ;
-
-        double maxX = 0;
-        double maxY = 0;
-        double minX = 0;
-        double minY = 0;
-        double minZ = 0;
-
-        double kx = 1;
-        double ky = 1;
-        double kz = 1;
-
-        if (which > 0)
-        {
-                for (int kt = 1; kt < (which + 1); kt++)
-                {
-                        kx = kx * kx_;
-                        ky = ky * ky_;
-                        kz = kz * kz_;
-                }
-        }
-
-        bool all_k_ones = kx_ == 1 && ky_ == 1 && kz_ == 1;
-
-        double rx = all_k_ones ? radius : kx;
-        double ry = all_k_ones ? radius : ky;
-        double rz = all_k_ones ? radius : kz;
+        double angle = M_PI * 2 / fF_LATITUDE_NUM;
+        for (int i = 0; i < fF_LATITUDE_NUM; i++)
+        {
+                cosines[i] = cos(i * angle);
+                sines[i] = sin(i * angle);
+        }
+}
+
+fF_point* GenoConv_fF0::generate_points(fF_chamber3d *chamber)
+{
+        double cenx = chamber->centerX;
+        double ceny = chamber->centerY;
+        double cenz = chamber->centerZ;
+
+        double rx = chamber->radius_x;
+        double ry = chamber->radius_y;
+        double rz = chamber->radius_z;
 
         fF_point *points = new fF_point[fF_SIZE];
@@ -227 +206 @@
                         p.z = z;
                         p.inside = false;
-
-                        if (x < minX) minX = x;
-                        if (x > maxX) maxX = x;
-                        if (y < minY) minY = y;
-                        if (y > maxY) maxY = y;
-
-                        if (z < minZ) minZ = z;
                 }
         }
@@ -240 +212 @@
 }
 
+template<typename T> T Square(T x) { return x * x; }
+
 double GenoConv_fF0::dist(double x1, double y1, double z1, double x2, double y2, double z2)
 {
-        return sqrt((x2 - x1)*(x2 - x1) + (y2 - y1)*(y2 - y1) + (z2 - z1)*(z2 - z1));
-}
-
-void GenoConv_fF0::search_hid(int nr, fF_chamber3d **spheres, double kx_, double ky_, double kz_)
-{
-        double kxsq = kx_*kx_;
-        double kysq = ky_*ky_;
-        double kzsq = kz_*kz_;
-
+        return sqrt(Square(x2 - x1) + Square(y2 - y1) + Square(z2 - z1));
+}
+
+void GenoConv_fF0::search_hid(int nr, fF_chamber3d **chambers)
+{
         for (int i = 0; i < nr; i++)
         {
-                double srX0 = spheres[i]->centerX;
-                double srY0 = spheres[i]->centerY;
-                double srZ0 = spheres[i]->centerZ;
-
-                double radsq = spheres[i]->radius * spheres[i]->radius;
-
-                double a2 = kx_ != 1 ? kxsq : radsq;
-                double b2 = ky_ != 1 ? kysq : radsq;
-                double c2 = kzsq * radsq;
+                fF_chamber3d *chamber = chambers[i];
+
+                double rx_sq = Square(chamber->radius_x);
+                double ry_sq = Square(chamber->radius_y);
+                double rz_sq = Square(chamber->radius_z);
 
                 for (int j = 0; j < fF_AMOUNT; j++)
                 {
-                        fF_point &p = spheres[nr]->points[j];
-
-                        double up1 = (p.x - srX0) * (p.x - srX0);
-                        double up2 = (p.y - srY0) * (p.y - srY0);
-                        double up3 = (p.z - srZ0) * (p.z - srZ0);
-
-                        double exp1 = up1 / a2;
-                        double exp2 = up2 / b2;
-                        double exp3 = up3 / c2;
-
-                        double result = exp1 + exp2 + exp3;
+                        fF_point &p = chambers[nr]->points[j];
+
+                        double upx = Square(p.x - chamber->centerX);
+                        double upy = Square(p.y - chamber->centerY);
+                        double upz = Square(p.z - chamber->centerZ);
+
+                        double expx = upx / rx_sq;
+                        double expy = upy / ry_sq;
+                        double expz = upz / rz_sq;
+
+                        double result = expx + expy + expz;
 
                         if (result < fF_THICK_RATIO)
@@ -285 +251 @@
 }
 
-int GenoConv_fF0::find_hole(int which, double x, double y, double z, fF_chamber3d **chambers, double kx_, double ky_, double kz_)
+int GenoConv_fF0::find_hole(int which, double x, double y, double z, fF_chamber3d **chambers)
 {
         int found = -1;
         double distsq_found;
 
-        double kxsq = kx_*kx_;
-        double kysq = ky_*ky_;
-        double kzsq = kz_*kz_;
-
         for (int i = 0; i < fF_AMOUNT; i++)
         {
@@ -299 +261 @@
                 if (!p.inside) //it is not inside another chamber
                 {
-                        double distancesq = (p.x - x)*(p.x - x) + (p.y - y)*(p.y - y) + (p.z - z)*(p.z - z);
+                        double distancesq = Square(p.x - x) + Square(p.y - y) + Square(p.z - z);
                         if (found < 0)
                         {
@@ -312 +274 @@
                                         for (int j = 0; j < which && good; j++)
                                         {
-                                                double srX0 = chambers[j]->centerX;
-                                                double srY0 = chambers[j]->centerY;
-                                                double srZ0 = chambers[j]->centerZ;
-
-                                                double radsq = chambers[j]->radius * chambers[j]->radius;
-
-                                                double a2 = kxsq * radsq;
-                                                double b2 = kysq * radsq;
-                                                double c2 = kzsq * radsq;
-
-                                                double up1 = (p.x - srX0) * (p.x - srX0);
-                                                double up2 = (p.y - srY0) * (p.y - srY0);
-                                                double up3 = (p.z - srZ0) * (p.z - srZ0);
-
-                                                double exp1 = up1 / a2;
-                                                double exp2 = up2 / b2;
-                                                double exp3 = up3 / c2;
-
-                                                double result = exp1 + exp2 + exp3;
+                                                fF_chamber3d *chamber = chambers[j];
+
+                                                double rx_sq = Square(chamber->radius_x);
+                                                double ry_sq = Square(chamber->radius_y);
+                                                double rz_sq = Square(chamber->radius_z);
+
+                                                double upx = Square(p.x - chamber->centerX);
+                                                double upy = Square(p.y - chamber->centerY);
+                                                double upz = Square(p.z - chamber->centerZ);
+
+                                                double expx = upx / rx_sq;
+                                                double expy = upy / ry_sq;
+                                                double expz = upz / rz_sq;
+
+                                                double result = expx + expy + expz;
                                                 if (result < 1.0)
                                                 {