source: cpp/frams/genetics/fF/conv_fF.cpp @ 736

Last change on this file since 736 was 736, checked in by Maciej Komosinski, 7 years ago

Added the new "using_checkpoints" argument to genetic converters so they can now call Model.checkpoint() when desired, see conv_f1.cpp for an example

  • Property svn:eol-style set to native
File size: 8.0 KB
Line 
1// This file is a part of Framsticks SDK.  http://www.framsticks.com/
2// Copyright (C) 1999-2018  Maciej Komosinski and Szymon Ulatowski.
3// See LICENSE.txt for details.
4
5#include "conv_fF.h"
6#include "fF_genotype.h"
7#include <common/nonstd_stl.h>
8#include <common/Convert.h>
9
10GenoConv_fF0::GenoConv_fF0()
11{
12        name = "10-parameter Foraminifera encoding";
13        in_format = 'F';
14        out_format = '0';
15        mapsupport = 0;
16        cosines = new double[fF_LATITUDE_NUM];
17        sines = new double[fF_LATITUDE_NUM];
18        precompute_cos_and_sin();
19}
20
21GenoConv_fF0::~GenoConv_fF0()
22{
23        delete[] cosines;
24        delete[] sines;
25}
26
27Part *GenoConv_fF0::addNewPart(Model *m, const fF_chamber3d* c)
28{
29        Part *part = m->addNewPart(Part::SHAPE_ELLIPSOID);
30        part->p = Pt3D(c->centerX, c->centerY, c->centerZ);
31        Pt3D hole = Pt3D(c->holeX, c->holeY, c->holeZ);
32        Orient o;
33        o.lookAt(part->p - hole);
34        part->setOrient(o);
35        return part;
36}
37
38SString GenoConv_fF0::convert(SString &in, MultiMap *map, bool using_checkpoints)
39{
40        fF_growth_params gp;
41        if (!gp.load(in.c_str())) //invalid input genotype?
42                return ""; //so we return an invalid f0 genotype
43
44        Model m;
45        m.open();
46
47        m.vis_style = "foram"; //dedicated visual look for Foraminifera
48
49        fF_chamber3d **chambers = new fF_chamber3d*[gp.number_of_chambers];
50        for (int i = 0; i < gp.number_of_chambers; i++)
51                createSphere(i, chambers, gp.radius0x, gp.radius0y, gp.radius0z, gp.translation, gp.angle1, gp.angle2, gp.scalex, gp.scaley, gp.scalez);
52
53        Part *p1 = addNewPart(&m, chambers[0]);
54        p1->scale = Pt3D(gp.radius0x, gp.radius0y, gp.radius0z); //size of the initial chamber
55        for (int i = 1; i < gp.number_of_chambers; i++)
56        {
57                Part *p2 = addNewPart(&m, chambers[i]);
58                p2->scale = p1->scale.entrywiseProduct(Pt3D(gp.scalex, gp.scaley, gp.scalez)); //each part's scale is its predecessor's scale * scaling
59                m.addNewJoint(p1, p2, Joint::SHAPE_FIXED); //all parts must be connected
60                p1 = p2;
61        }
62
63        for (int i = 0; i < gp.number_of_chambers; i++)
64                delete chambers[i];
65        delete[] chambers;
66
67        m.close();
68        return m.getF0Geno().getGenes();
69}
70
71void 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_)
72{
73        chambers[which] = new fF_chamber3d(0.0, 0.0, 0.0,
74                radius0x, radius0y, radius0z, radius0x * kx_, 0.0, 0.0,
75                radius0x * translation, 0.0, 0.0, 0.0, 0.0);
76        if (which == 0)
77                chambers[which]->points = generate_points(chambers[which]);
78        if (which > 0)
79        {
80                chambers[which]->radius_x = get_radius(chambers[which - 1]->radius_x, kx_, chambers[0]->radius_x);
81                chambers[which]->radius_y = get_radius(chambers[which - 1]->radius_y, ky_, chambers[0]->radius_y);
82                chambers[which]->radius_z = get_radius(chambers[which - 1]->radius_z, kz_, chambers[0]->radius_z);
83
84                /* new growth vector length */
85                double len = chambers[which]->radius_y * translation;
86                double max_radius = fF_TOO_MUCH * chambers[which]->radius_y;
87                if (fabs(len) > (max_radius))
88                        len = ((len < 0) ? (-1) : 1) * max_radius;
89                if (len == 0)
90                        len = -0.0000001;
91
92                /* aperture of the previous chamber */
93                double pzx = chambers[which - 1]->holeX;
94                double pzy = chambers[which - 1]->holeY;
95                double pzz = chambers[which - 1]->holeZ;
96
97                //center of the previous chamber
98                double pcx = chambers[which - 1]->centerX;
99                double pcy = chambers[which - 1]->centerY;
100                //double pcz = chambers[which - 1]->centerZ; //not used
101
102                /* aperture of the next to last chamber */
103                double ppx;
104                double ppy;
105                //double ppz; //not used
106
107                if (which == 1)
108                {
109                        ppx = pcx;
110                        ppy = pcy;
111                        //ppz = pcz;
112                }
113                else
114                {
115                        ppx = chambers[which - 2]->holeX;
116                        ppy = chambers[which - 2]->holeY;
117                        //ppz = chambers[which - 2]->holeZ;
118                }
119
120                double pzxprim = pzx - ppx;
121                double pzyprim = pzy - ppy;
122                double angle;
123
124                angle = Convert::atan_2(pzyprim, pzxprim);
125                double alpha = angle - alpha_;
126
127                double gamma = chambers[which - 1]->phi + gamma_;
128
129                double wx = len * cos(alpha);
130                double wy = len * sin(alpha);
131                double wz = len * sin(alpha) * sin(gamma);
132
133                /*center of the new sphere*/
134                double x = pzx + wx;
135                double y = pzy + wy;
136                double z = pzz + wz;
137
138                chambers[which]->centerX = x;
139                chambers[which]->centerY = y;
140                chambers[which]->centerZ = z;
141                chambers[which]->vectorTfX = wx;
142                chambers[which]->vectorTfY = wy;
143                chambers[which]->vectorTfZ = wz;
144                chambers[which]->beta = alpha;
145                chambers[which]->phi = gamma;
146
147                chambers[which]->points = generate_points(chambers[which]);
148                search_hid(which, chambers);
149                int pun;
150                pun = find_hole(which, pzx, pzy, pzz, chambers);
151                if (pun < 0) //should never happen
152                {
153                        logPrintf("GenoConv_fF0", "createSphere", LOG_ERROR, "find_hole(%d) returned %d", which, pun);
154                        pun = 0;
155                }
156                chambers[which]->holeX = chambers[which]->points[pun].x;
157                chambers[which]->holeY = chambers[which]->points[pun].y;
158                chambers[which]->holeZ = chambers[which]->points[pun].z;
159        }
160}
161
162double GenoConv_fF0::get_radius(double prev_radius, double scale, double start_radius)
163{
164        double radius = prev_radius * scale;
165        double min_radius = fF_TOO_LITTLE*start_radius;
166        if (radius < min_radius) {
167                radius = min_radius;
168        }
169
170        return radius;
171}
172
173void GenoConv_fF0::precompute_cos_and_sin()
174{
175        double angle = M_PI * 2 / fF_LATITUDE_NUM;
176        for (int i = 0; i < fF_LATITUDE_NUM; i++)
177        {
178                cosines[i] = cos(i * angle);
179                sines[i] = sin(i * angle);
180        }
181}
182
183fF_point* GenoConv_fF0::generate_points(fF_chamber3d *chamber)
184{
185        double cenx = chamber->centerX;
186        double ceny = chamber->centerY;
187        double cenz = chamber->centerZ;
188
189        double rx = chamber->radius_x;
190        double ry = chamber->radius_y;
191        double rz = chamber->radius_z;
192
193        fF_point *points = new fF_point[fF_SIZE];
194
195        for (int i = 0; i < fF_LONGITUDE_NUM; i++)
196        {
197                double y = ceny + ry * cosines[i];
198
199                for (int j = 0; j < fF_LATITUDE_NUM; j++)
200                {
201                        double x = cenx + rx * cosines[j] * sines[i];
202                        double z = cenz + rz * sines[j] * sines[i];
203                        fF_point &p = points[(i * fF_LATITUDE_NUM) + j];
204                        p.x = x;
205                        p.y = y;
206                        p.z = z;
207                        p.inside = false;
208                }
209        }
210        return points;
211
212}
213
214template<typename T> T Square(T x) { return x * x; }
215
216double GenoConv_fF0::dist(double x1, double y1, double z1, double x2, double y2, double z2)
217{
218        return sqrt(Square(x2 - x1) + Square(y2 - y1) + Square(z2 - z1));
219}
220
221void GenoConv_fF0::search_hid(int nr, fF_chamber3d **chambers)
222{
223        for (int i = 0; i < nr; i++)
224        {
225                fF_chamber3d *chamber = chambers[i];
226
227                double rx_sq = Square(chamber->radius_x);
228                double ry_sq = Square(chamber->radius_y);
229                double rz_sq = Square(chamber->radius_z);
230
231                for (int j = 0; j < fF_AMOUNT; j++)
232                {
233                        fF_point &p = chambers[nr]->points[j];
234
235                        double upx = Square(p.x - chamber->centerX);
236                        double upy = Square(p.y - chamber->centerY);
237                        double upz = Square(p.z - chamber->centerZ);
238
239                        double expx = upx / rx_sq;
240                        double expy = upy / ry_sq;
241                        double expz = upz / rz_sq;
242
243                        double result = expx + expy + expz;
244
245                        if (result < fF_THICK_RATIO)
246                        {
247                                p.inside = true;
248                        }
249                }
250        }
251}
252
253int GenoConv_fF0::find_hole(int which, double x, double y, double z, fF_chamber3d **chambers)
254{
255        int found = -1;
256        double distsq_found;
257
258        for (int i = 0; i < fF_AMOUNT; i++)
259        {
260                fF_point &p = chambers[which]->points[i];
261                if (!p.inside) //it is not inside another chamber
262                {
263                        double distancesq = Square(p.x - x) + Square(p.y - y) + Square(p.z - z);
264                        if (found < 0)
265                        {
266                                found = i;
267                                distsq_found = distancesq;
268                        }
269                        if (distancesq < distsq_found)
270                        {
271                                if (which != 0)
272                                {
273                                        bool good = true;
274                                        for (int j = 0; j < which && good; j++)
275                                        {
276                                                fF_chamber3d *chamber = chambers[j];
277
278                                                double rx_sq = Square(chamber->radius_x);
279                                                double ry_sq = Square(chamber->radius_y);
280                                                double rz_sq = Square(chamber->radius_z);
281
282                                                double upx = Square(p.x - chamber->centerX);
283                                                double upy = Square(p.y - chamber->centerY);
284                                                double upz = Square(p.z - chamber->centerZ);
285
286                                                double expx = upx / rx_sq;
287                                                double expy = upy / ry_sq;
288                                                double expz = upz / rz_sq;
289
290                                                double result = expx + expy + expz;
291                                                if (result < 1.0)
292                                                {
293                                                        good = false;
294                                                }
295                                        }
296                                        if (good)
297                                        {
298                                                found = i;
299                                                distsq_found = distancesq;
300                                        }
301                                }
302                        }
303                }
304        }
305
306        return found;
307}
Note: See TracBrowser for help on using the repository browser.