source: cpp/frams/genetics/genooperators.cpp @ 1304

// This file is a part of Framsticks SDK.  http://www.framsticks.com/
// Copyright (C) 1999-2023  Maciej Komosinski and Szymon Ulatowski.
// See LICENSE.txt for details.

#include <ctype.h>  //isupper()
#include <algorithm> // std::min, std::max
#include <cmath> // std::floor()
#include "genooperators.h"
#include <common/log.h>
#include <common/nonstd_math.h>
#include <frams/util/rndutil.h>

//
// custom distributions for mutations of various parameters
//
static double distrib_force[] =   // for '!'
{
        3,             // distribution 0 -__/ +1
        0.001, 0.2,    // "slow" neurons
        0.001, 1,
        1, 1,          // "fast" neurons
};
static double distrib_inertia[] =  // for '='
{
        2,             // distribution 0 |..- +1
        0, 0,          // "fast" neurons
        0.7, 0.98,
};
static double distrib_sigmo[] =  // for '/'
{
        5,             // distribution -999 -..-^-..- +999
        -999, -999,    //"perceptron"
        999, 999,
        -5, -1,        // nonlinear
        1, 5,
        -1, 1,         // ~linear
};
/*
static double distrib_weight[] =
{
5,                 // distribution -999 _-^_^-_ +999
-999, 999,         // each weight value may be useful, especially...
-5, -0.3,          // ...little non-zero values
-3, -0.6,
0.6, 3,
0.3, 5,
};
*/

int GenoOperators::roulette(const double *probtab, const int count)
{
        double sum = 0;
        int i;
        for (i = 0; i < count; i++) sum += probtab[i];
        double sel = rndDouble(sum);
        for (sum = 0, i = 0; i < count; i++) { sum += probtab[i]; if (sel < sum) return i; }
        return -1;
}

bool GenoOperators::getMinMaxDef(ParamInterface *p, int i, double &mn, double &mx, double &def)
{
        mn = mx = def = 0;
        int defined = 0;
        if (p->type(i)[0] == 'f')
        {
                double _mn = 0, _mx = 1, _def = 0.5;
                defined = p->getMinMaxDouble(i, _mn, _mx, _def);
                if (defined == 1) _mx = _mn + 1000.0; //only min was defined, so let's set some arbitrary range, just to have some freedom. Assumes _mn is not close to maxdouble...
                if (_mx < _mn && defined == 3) //only default was defined, so let's assume some arbitrary range. Again, no check for min/maxdouble...
                {
                        _mn = _def - 500.0;
                        _mx = _def + 500.0;
                }
                if (defined < 3) _def = (_mn + _mx) / 2.0;
                mn = _mn; mx = _mx; def = _def;
        }
        if (p->type(i)[0] == 'd')
        {
                paInt _mn = 0, _mx = 1, _def = 0;
                defined = p->getMinMaxInt(i, _mn, _mx, _def);
                if (defined == 1) _mx = _mn + 1000; //only min was defined, so let's set some arbitrary range, just to have some freedom. Assumes _mn is not close to maxint...
                if (_mx < _mn && defined == 3) //only default was defined, so let's assume some arbitrary range. Again, no check for min/maxint...
                {
                        _mn = _def - 500;
                        _mx = _def + 500;
                }
                if (defined < 3) _def = (_mn + _mx) / 2;
                mn = _mn; mx = _mx; def = _def;
        }
        return defined == 3;
}

bool GenoOperators::mutateRandomNeuroClassProperty(Neuro* n)
{
        bool mutated = false;
        int prop = selectRandomNeuroClassProperty(n);
        if (prop >= 0)
        {
                if (prop >= GenoOperators::NEUROCLASS_PROP_OFFSET)
                {
                        SyntParam par = n->classProperties();   //commits changes when this object is destroyed
                        mutated = mutateProperty(par, prop - GenoOperators::NEUROCLASS_PROP_OFFSET);
                }
                else
                {
                        Param par = n->extraProperties();
                        mutated = mutateProperty(par, prop);
                }
        }
        return mutated;
}

int GenoOperators::selectRandomNeuroClassProperty(Neuro *n)
{
        int neuext = n->extraProperties().getPropCount(),
                neucls = n->getClass() == NULL ? 0 : n->getClass()->getProperties().getPropCount();
        if (neuext + neucls == 0) return -1; //no properties in this neuron
        int index = rndUint(neuext + neucls);
        if (index >= neuext) index = index - neuext + NEUROCLASS_PROP_OFFSET;
        return index;
}

double GenoOperators::getMutatedNeuroClassProperty(double current, Neuro *n, int i)
{
        if (i == -1)
        {
                logPrintf("GenoOperators", "getMutatedNeuroClassProperty", LOG_WARN, "Deprecated usage in C++ source: to mutate connection weight, use getMutatedNeuronConnectionWeight().");
                return getMutatedNeuronConnectionWeight(current);
        }
        Param p;
        if (i >= NEUROCLASS_PROP_OFFSET) { i -= NEUROCLASS_PROP_OFFSET; p = n->getClass()->getProperties(); }
        else p = n->extraProperties();
        double newval = current;
        /*bool ok=*/getMutatedProperty(p, i, current, newval);
        return newval;
}

double GenoOperators::getMutatedNeuronConnectionWeight(double current)
{
        return mutateCreepNoLimit('f', current, 2, true);
}

bool GenoOperators::mutatePropertyNaive(ParamInterface &p, int i)
{
        double mn, mx, df;
        if (p.type(i)[0] != 'f' && p.type(i)[0] != 'd') return false; //don't know how to mutate
        getMinMaxDef(&p, i, mn, mx, df);

        ExtValue ev;
        p.get(i, ev);
        ev.setDouble(mutateCreep(p.type(i)[0], ev.getDouble(), mn, mx, true));
        p.set(i, ev);
        return true;
}

bool GenoOperators::mutateProperty(ParamInterface &p, int i)
{
        double newval;
        ExtValue ev;
        p.get(i, ev);
        bool ok = getMutatedProperty(p, i, ev.getDouble(), newval);
        if (ok) { ev.setDouble(newval); p.set(i, ev); }
        return ok;
}

bool GenoOperators::getMutatedProperty(ParamInterface &p, int i, double oldval, double &newval)
{
        newval = 0;
        if (p.type(i)[0] != 'f' && p.type(i)[0] != 'd') return false; //don't know how to mutate
        const char *n = p.id(i), *na = p.name(i);
        if (strcmp(n, "si") == 0 && strcmp(na, "Sigmoid") == 0) newval = round(CustomRnd(distrib_sigmo), 3); else
                if (strcmp(n, "in") == 0 && strcmp(na, "Inertia") == 0) newval = round(CustomRnd(distrib_inertia), 3); else
                        if (strcmp(n, "fo") == 0 && strcmp(na, "Force") == 0) newval = round(CustomRnd(distrib_force), 3); else
                        {
                                double mn, mx, df;
                                getMinMaxDef(&p, i, mn, mx, df);
                                newval = mutateCreep(p.type(i)[0], oldval, mn, mx, true);
                        }
        return true;
}

double GenoOperators::mutateCreepNoLimit(char type, double current, double stddev, bool limit_precision_3digits)
{
        double result = RndGen.Gauss(current, stddev);
        if (type == 'd')
        {
                result = int(result + 0.5);
                if (result == current) result += rndUint(2) * 2 - 1; //force some change
        }
        else
        {
                if (limit_precision_3digits)
                        result = round(result, 3);
        }
        return result;
}

double GenoOperators::mutateCreep(char type, double current, double mn, double mx, double stddev, bool limit_precision_3digits)
{
        double result = mutateCreepNoLimit(type, current, stddev, limit_precision_3digits);
        if (result<mn || result>mx) //exceeds boundary, so bring to the allowed range
        {
                //reflect:
                if (result > mx) result = mx - (result - mx); else
                        if (result < mn) result = mn + (mn - result);
                //wrap (just in case 'result' exceeded the allowed range so much that after the reflection above it exceeded the other boundary):
                if (result > mx) result = mn + fmod(result - mx, mx - mn); else
                        if (result < mn) result = mn + fmod(mn - result, mx - mn);
                if (limit_precision_3digits)
                {
                        //reflect and wrap above may have changed the (limited) precision, so try to round again (maybe unnecessarily, because we don't know if reflect+wrap above were triggered)
                        double result_try = round(result, 3);
                        if (mn <= result_try && result_try <= mx) result = result_try; //after rounding still within allowed range, so keep rounded value
                }
        }
        clipNegativeZeroIfNeeded(result, mn); //so we don't get -0.0 when minimum is 0.0
        return result;
}

double GenoOperators::mutateCreep(char type, double current, double mn, double mx, bool limit_precision_3digits)
{
        double stddev = (mx - mn) / 2 / 5; // magic arbitrary formula for stddev, which becomes /halfinterval, 5 times narrower
        return mutateCreep(type, current, mn, mx, stddev, limit_precision_3digits);
}

void GenoOperators::setIntFromDoubleWithProbabilisticDithering(ParamInterface &p, int index, double value)
{
        // Deterministic rounding to the closest integer:
        //value += 0.5; // value==2.499 will become int 2 and value==2.5 will become int 3, but we want these cases to be 2 or 3 with almost equal probability (stochastic rounding).

        //stochastic rounding (value==2.1 should turn in most cases to int 2, rarely to int 3; value==-2.1 should become mostly int -2, rarely int -3):
        double lower = std::floor(value);
        value = rndDouble(1) < (value - lower) ? lower + 1 : lower;

        p.setInt(index, (paInt)value);
}

void GenoOperators::linearMix(vector<double> &p1, vector<double> &p2, double proportion)
{
        if (p1.size() != p2.size())
        {
                logPrintf("GenoOperators", "linearMix", LOG_ERROR, "Cannot mix vectors of different length (%d and %d)", p1.size(), p2.size());
                return;
        }
        for (unsigned int i = 0; i < p1.size(); i++)
        {
                double v1 = p1[i];
                double v2 = p2[i];
                p1[i] = v1 * proportion + v2 * (1 - proportion);
                p2[i] = v2 * proportion + v1 * (1 - proportion);
        }
}

void GenoOperators::linearMix(ParamInterface &p1, int i1, ParamInterface &p2, int i2, double proportion)
{
        char type1 = p1.type(i1)[0];
        char type2 = p2.type(i2)[0];
        if (type1 == 'f' && type2 == 'f')
        {
                double v1 = p1.getDouble(i1);
                double v2 = p2.getDouble(i2);
                p1.setDouble(i1, v1 * proportion + v2 * (1 - proportion));
                p2.setDouble(i2, v2 * proportion + v1 * (1 - proportion));
        }
        else
                if (type1 == 'd' && type2 == 'd')
                {
                        int v1 = p1.getInt(i1);
                        int v2 = p2.getInt(i2);
                        setIntFromDoubleWithProbabilisticDithering(p1, i1, v1 * proportion + v2 * (1 - proportion));
                        setIntFromDoubleWithProbabilisticDithering(p2, i2, v2 * proportion + v1 * (1 - proportion));
                }
                else
                        logPrintf("GenoOperators", "linearMix", LOG_WARN, "Cannot mix values of types '%c' and '%c'", type1, type2);
}

int GenoOperators::getActiveNeuroClassCount(Model::ShapeType for_shape_type)
{
        int count = 0;
        for (int i = 0; i < Neuro::getClassCount(); i++)
        {
                NeuroClass *nc = Neuro::getClass(i);
                if (nc->isShapeTypeSupported(for_shape_type) && nc->genactive)
                        count++;
        }
        return count;
}

NeuroClass *GenoOperators::getRandomNeuroClass(Model::ShapeType for_shape_type)
{
        vector<NeuroClass *> active;
        for (int i = 0; i < Neuro::getClassCount(); i++)
        {
                NeuroClass *nc = Neuro::getClass(i);
                if (nc->isShapeTypeSupported(for_shape_type) && nc->genactive)
                        active.push_back(nc);
        }
        if (active.size() == 0) return NULL; else return active[rndUint(active.size())];
}

NeuroClass *GenoOperators::getRandomNeuroClassWithOutput(Model::ShapeType for_shape_type)
{
        vector<NeuroClass *> active;
        for (int i = 0; i < Neuro::getClassCount(); i++)
        {
                NeuroClass *nc = Neuro::getClass(i);
                if (nc->isShapeTypeSupported(for_shape_type) && nc->genactive && nc->getPreferredOutput() != 0)
                        active.push_back(nc);
        }
        if (active.size() == 0) return NULL; else return active[rndUint(active.size())];
}

NeuroClass *GenoOperators::getRandomNeuroClassWithInput(Model::ShapeType for_shape_type)
{
        vector<NeuroClass *> active;
        for (int i = 0; i < Neuro::getClassCount(); i++)
        {
                NeuroClass *nc = Neuro::getClass(i);
                if (nc->isShapeTypeSupported(for_shape_type) && nc->genactive && nc->getPreferredInputs() != 0)
                        active.push_back(nc);
        }
        if (active.size() == 0) return NULL; else return active[rndUint(active.size())];
}

NeuroClass *GenoOperators::getRandomNeuroClassWithOutputAndWantingNoInputs(Model::ShapeType for_shape_type)
{
        vector<NeuroClass *> active;
        for (int i = 0; i < Neuro::getClassCount(); i++)
        {
                NeuroClass *nc = Neuro::getClass(i);
                if (nc->isShapeTypeSupported(for_shape_type) && nc->genactive && nc->getPreferredOutput() != 0 && nc->getPreferredInputs() == 0)
                        active.push_back(nc);
        }
        if (active.size() == 0) return NULL; else return active[rndUint(active.size())];
}

NeuroClass *GenoOperators::getRandomNeuroClassWithOutputAndWantingNoOrAnyInputs(Model::ShapeType for_shape_type)
{
        vector<NeuroClass *> active;
        for (int i = 0; i < Neuro::getClassCount(); i++)
        {
                NeuroClass *nc = Neuro::getClass(i);
                if (nc->isShapeTypeSupported(for_shape_type) && nc->genactive && nc->getPreferredOutput() != 0 && nc->getPreferredInputs() <= 0) // getPreferredInputs() should be 0 or -1 (any)
                        active.push_back(nc);
        }
        if (active.size() == 0) return NULL; else return active[rndUint(active.size())];
}

int GenoOperators::getRandomNeuroClassWithOutput(const vector<NeuroClass *> &NClist)
{
        vector<int> allowed;
        for (int i = 0; i < (int)NClist.size(); i++)
                if (NClist[i]->getPreferredOutput() != 0) //this NeuroClass provides output
                        allowed.push_back(i);
        if (allowed.size() == 0) return -1; else return allowed[rndUint(allowed.size())];
}

int GenoOperators::getRandomNeuroClassWithInput(const vector<NeuroClass *> &NClist)
{
        vector<int> allowed;
        for (int i = 0; i < (int)NClist.size(); i++)
                if (NClist[i]->getPreferredInputs() != 0) //this NeuroClass wants one input connection or more                  
                        allowed.push_back(i);
        if (allowed.size() == 0) return -1; else return allowed[rndUint(allowed.size())];
}

NeuroClass *GenoOperators::parseNeuroClass(char *&s, ModelEnum::ShapeType for_shape_type)
{
        int maxlen = (int)strlen(s);
        int NClen = 0;
        NeuroClass *NC = NULL;
        for (int i = 0; i < Neuro::getClassCount(); i++)
        {
                NeuroClass *nci = Neuro::getClass(i);
                if (!nci->isShapeTypeSupported(for_shape_type))
                        continue;
                const char *nciname = nci->name.c_str();
                int ncinamelen = (int)strlen(nciname);
                if (maxlen >= ncinamelen && ncinamelen > NClen && (strncmp(s, nciname, ncinamelen) == 0))
                {
                        NC = nci;
                        NClen = ncinamelen;
                }
        }
        s += NClen;
        return NC;
}

Neuro *GenoOperators::findNeuro(const Model *m, const NeuroClass *nc)
{
        if (!m) return NULL;
        for (int i = 0; i < m->getNeuroCount(); i++)
                if (m->getNeuro(i)->getClass() == nc) return m->getNeuro(i);
        return NULL; //neuron of class 'nc' was not found
}

int GenoOperators::neuroClassProp(char *&s, NeuroClass *nc, bool also_v1_N_props)
{
        int len = (int)strlen(s);
        int Len = 0, I = -1;
        if (nc)
        {
                Param p = nc->getProperties();
                for (int i = 0; i < p.getPropCount(); i++)
                {
                        const char *n = p.id(i);
                        int l = (int)strlen(n);
                        if (len >= l && l > Len && (strncmp(s, n, l) == 0)) { I = NEUROCLASS_PROP_OFFSET + i; Len = l; }
                        if (also_v1_N_props) //recognize old symbols of properties:  /=!
                        {
                                if (strcmp(n, "si") == 0) n = "/"; else
                                        if (strcmp(n, "in") == 0) n = "="; else
                                                if (strcmp(n, "fo") == 0) n = "!";
                                l = (int)strlen(n);
                                if (len >= l && l > Len && (strncmp(s, n, l) == 0)) { I = NEUROCLASS_PROP_OFFSET + i; Len = l; }
                        }
                }
        }
        Neuro n;
        Param p = n.extraProperties();
        for (int i = 0; i < p.getPropCount(); i++)
        {
                const char *n = p.id(i);
                int l = (int)strlen(n);
                if (len >= l && l > Len && (strncmp(s, n, l) == 0)) { I = i; Len = l; }
        }
        s += Len;
        return I;
}

bool GenoOperators::canStartNeuroClassName(const char firstchar)
{
        return isupper(firstchar) || firstchar == '|' || firstchar == '@' || firstchar == '*';
}

bool GenoOperators::isWS(const char c)
{
        return c == ' ' || c == '\n' || c == '\t' || c == '\r';
}

void GenoOperators::skipWS(char *&s)
{
        if (s == NULL)
                logMessage("GenoOperators", "skipWS", LOG_WARN, "NULL reference!");
        else
                while (isWS(*s)) s++;
}

bool GenoOperators::areAlike(char *g1, char *g2)
{
        while (*g1 || *g2)
        {
                skipWS(g1);
                skipWS(g2);
                if (*g1 != *g2) return false; //when difference
                if (!*g1 && !*g2) break; //both end
                g1++;
                g2++;
        }
        return true; //equal
}

char *GenoOperators::strchrn0(const char *str, char ch)
{
        return ch == 0 ? NULL : strchr((char *)str, ch);
}

int GenoOperators::getRandomChar(const char *choices, const char *excluded)
{
        int allowed_count = 0;
        for (size_t i = 0; i < strlen(choices); i++) if (!strchrn0(excluded, choices[i])) allowed_count++;
        if (allowed_count == 0) return -1; //no char is allowed
        int rnd_index = rndUint(allowed_count) + 1;
        allowed_count = 0;
        for (size_t i = 0; i < strlen(choices); i++)
        {
                if (!strchrn0(excluded, choices[i])) allowed_count++;
                if (allowed_count == rnd_index) return int(i);
        }
        return -1; //never happens
}

string GenoOperators::simplifiedModifiers_rR(const string& str)
{
        int R = 0; //positive means more 'R', negative means more 'r'
        for (char c : str)
        {
                if (c == 'R') R++; else
                        if (c == 'r') R--;
        }
        R %= 8; // 8 * 45 degrees = 360 degrees. After this, we get R=-7..+7

        /* now, simplify homogeneous sequences of rR longer than 4: for example, rrrrr == RRR and RRRRRR == rr
        -7      1
        -6      2
        -5      3
        -4      -4 (or 4; we choose +4 meaning we will never see rrrr)
        -3..3      (no changes)
        4       4 (or -4)
        5       -3
        6       -2
        7       -1
        */
        if (R <= -4) R += 8; //-4 => +4
        else if (R >= 5) R -= 8;

        return R == 0 ? "" : (R > 0 ? string(R, 'R') : string(-R, 'r'));
}

//#include <cassert>
string GenoOperators::simplifiedModifiersFixedOrder(const char *str_of_char_pairs, vector<int> &char_counts)
{
        //      assert(strlen(str_of_char_pairs) == char_counts.size());
        //      assert(char_counts.size() % 2 == 0);
        const int MAX_NUMBER_SAME_TYPE = 8; // max. number of modifiers of each type (case-sensitive) - mainly for rR, even though for rR, 4 would be sufficient if we assume lower or upper can be chosen as required for minimal length just as simplifiedModifiers_rR() does, e.g. rrrrr==RRR, RRRRRR==rr
        string simplified;
        //#define CLUMP_IDENTICAL_MODIFIERS //if GeneProps::normalizeBiol4() is used, this is not good because properties are calculated incrementally, non-linearly, their values are updated after each modifier character and some properties interact with each other due to normalization so they can saturate when clumped, therefore it is better keep the modifiers dispersed to equalize their effects
#ifdef CLUMP_IDENTICAL_MODIFIERS
        for (size_t i = 0; i < strlen(str_of_char_pairs); i++)
                if ((i % 2) == 0) //only even index "i" in str_of_char_pairs
                        for (int j = 0; j < std::min(MAX_NUMBER_SAME_TYPE, abs(char_counts[i] - char_counts[i + 1])); j++) //assume that an even-index char and the following odd-index char have the opposite influence, so they cancel out.
                                simplified += str_of_char_pairs[i + (char_counts[i + 1] > char_counts[i])]; //inner loop adds a sequence of same chars such as rrrrr or QQQ
#else
        for (size_t i = 0; i < strlen(str_of_char_pairs); i++)
                if ((i % 2) == 0) //only even index "i" in str_of_char_pairs
                {
                        char_counts[i] -= char_counts[i + 1]; //from now on, even items in the vector store the difference between antagonistic modifier symbols; odd items are not needed
                        char_counts[i] = std::min(std::max(char_counts[i], -MAX_NUMBER_SAME_TYPE), MAX_NUMBER_SAME_TYPE);
                }
        int remaining;
        do {
                remaining = 0;
                for (size_t i = 0; i < strlen(str_of_char_pairs); i++)
                        if ((i % 2) == 0) //only even index "i" in str_of_char_pairs
                                if (char_counts[i] != 0)
                                {
                                        simplified += str_of_char_pairs[i + (char_counts[i] < 0)];
                                        char_counts[i] += char_counts[i] > 0 ? -1 : +1; //decrease the difference towards zero
                                        remaining += abs(char_counts[i]);
                                }
        } while (remaining > 0);
#endif
        return simplified;
}

string GenoOperators::simplifiedModifiers(const string & original)
{
        const int MAX_NUMBER_SAME_TYPE = 5; // max. number of modifiers of each type (case-insensitive). The more characters, the closer we can get to min and max values of a given property at the expense of the length of evolved genotypes. 5 is "close enough", but how close we get to the extreme also depends on the initial value of a given property, which is not always exactly in the middle of min and max. rR is treated separately in simplification because their influence follows different (i.e., simple additive) logic - so the simplifiedModifiersFixedOrder() logic with cancelling out antagonistic modifiers would be appropriate for rR.
        int counter[256] = {}; //initialize with zeros; 256 is unnecessarily too big and redundant, but enables very fast access (indexed directly by the ascii code)
        string simplified = "";
        for (int i = int(original.size()) - 1; i >= 0; i--) //iterate from end to begin so it is easier to remove "oldest" = first modifiers
        {
                unsigned char c = original[i];
                if (!std::isalpha(c) || c == 'r' || c == 'R') //ignore non-alphabet characters; also, 'r' and 'R' are handled separately by simplifiedModifiers_rR()
                        continue;
                unsigned char lower = std::tolower(c);
                counter[lower]++;
                if (counter[lower] <= MAX_NUMBER_SAME_TYPE) //get rid of modifiers that are too numerous, but get rid of the first ones in the string (="oldest", the last ones looking from the end), because their influence on the parameter value is the smallest
                        simplified += c;
        }
        std::reverse(simplified.begin(), simplified.end()); //"simplified" was built in reverse order, so need to restore the order that corresponds to "original"
        return simplifiedModifiers_rR(original) + simplified;
}