Thread: loss of an object between the very end of a loop and the very beginning

Hello,
I spent a lot of times the last week on my code because I get from its execution:

Code:

Erreur de segmentation (core dumped)

It only happend when the parameter nDemes has a value greater than 2.


When I put some "printf" everywhere in the code (through a function named "void afficherPopulation()", showPopulation() in french) I found something very weird.
The main has the following structure (briefly):

Code:

int main(int argc, char* argv[]){
    int i = 0;
    for(i=0; i<x; i++){
        afficherPopulation(population, ... ); // it works if i=0 but not when i=1
        ...
        different functions generating "newPopulation" from "population"
        ...
        copy/paste of values contained within "newPopulation" within the re-created "population"


        afficherPopulation(population, ... ); // it only works when i=0
    }
}

At the very end of the first loop, I correctly have "population" because everything inside are correctly printed when i=0, but at the very beginning of the second iteration of the loop, the exact same function using the exact same objects doesn't work anymore.
Population is an array of "deme" (a created structure) of length nDemes. If nDemes is equal to 1 or 2 (first argument), it works:

Code:

./a.out 2 5 10 5 0.0001 2 0.0001 2 0.2 0.1 1 126

but it's not the case for i=3:

Code:

./a.out 3 5 10 5 0.0001 2 0.0001 2 0.2 0.1 1 126

I really don't understand why.


Do you have an idea on that please ?


Thank you very much,


M.


P.S. I put my code here despite its length...

Code:

// gcc *.c -L/usr/local/lib  -lgsl -lgslcblas -lm
// ./a.out 500 200 10000 20 0.00001 2 0.00001 10 0.001 0.02 2 123
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_statistics.h>
#include <gsl/gsl_permutation.h>
#define MAX_NUMBER_OF_INITIAL_NTRL_ALLELES 100


typedef struct Deme Deme;
struct Deme{
    int nIndividus;
    long* ntrlLoci;    // contains the neutral alleles for nIndividus (number of individuals within the deme) x 2 (diploids) x nNtlrLoci (number of neutral loci)
    double* quantiLoci; // contains the allelic effects for nIndividus (number of individuals within the deme) x 2 (diploids) x nQuantiLoci (number of quantitative loci)
    double* femaleAllocation; // =sum of the allelic effects over the quantitative loci
    double* maleAllocation; // =(1 - femaleAllocation)
    int* nOffsprings; // floor(X) + Binom(n=1, p=X-floor(x)) where X = femaleAllocation x fecundity
};


void initializePopulation(gsl_rng *r, Deme* population, int nDemes, int maxIndPerDem, int nNtrlLoci, int nQuantiLoci, int fecundity);
void libererMemoirePopulation(Deme* population, int nDemes);
void afficherPopulation(Deme* population, int nDemes, int nNtrlLoci, int nQuantiLoci);
void configMetapop(gsl_rng* r, Deme* population, int nDemes, const int maxIndPerDem, double migration, double extinction, int nImmigrants[], int extinctionStatus[], int nProducedSeeds[]);
void setToZero(int nDemes, int nImmigrants[], int extinctionStatus[], int nProducedSeeds[]);
void initializeNewPopulation(Deme* newPopulation, const int nDemes, const int maxIndPerDem, const int nImmigrants[], const int nProducedSeeds[], int extinctionStatus[], const int nNtrlLoci, const int nQuantiLoci, const int recolonization, const int fecundity);
void panmixie(gsl_rng* r, Deme* population, Deme* newPopulation, const int nDemes, const int maxIndPerDem, const int nProducedSeeds[], const int nNtrlLoci, const int nQuantiLoci, long* pointeurSurListOfNtrlAlleles, double ntrlMutation, double quantiMutation, int fecundity);
void weightedSample(gsl_rng* r, const double* liste, const double* weights, double* target, const int sizeOfListe, const int nTrials);
void replacement(gsl_rng* r, Deme* population, Deme* newPopulation, const int nDemes, const maxIndPerDem, const int nNtrlLoci, const int nQuantiLoci, int fecundity, const int nImmigrants[], int nProducedSeeds[], int extinctionStatus[], const int recolonization);


int main(int argc, char *argv[]){
    int i = 0;
    int j = 0;
    int tmp = 0;


    // Get Parameters from comamnd line
    const int nDemes = atoi(argv[1]); // number of demes
    const int maxIndPerDem = atoi(argv[2]);    // carrying capacity per deme
    const int nGeneration = atoi(argv[3]);    // number of generations to simulate
    const int nNtrlLoci = atoi(argv[4]);    // number of neutral loci
    const double ntrlMutation = atof(argv[5]);    // mutation rate of the ntrl loci
    const int nQuantiLoci = atoi(argv[6]);    // number of quantitative loci
    const double quantiMutation = atof(argv[7]);    // mutation rate of the quantative loci
    const int fecundity = atoi(argv[8]);    // max number of offspring when femAlloc=100%
    const double migration = atof(argv[9]);    // immigration rate
    const double extinction = atof(argv[10]);    // extinction rate
    const int recolonization = atoi(argv[11]);    // number of recolonizing individuals
    const int seed = atoi(argv[12]);


    // Random generator
        const gsl_rng_type *T;
        gsl_rng *r;
        gsl_rng_env_setup();


        T=gsl_rng_default;
        r=gsl_rng_alloc(T);
        gsl_rng_set(r, seed);


    // Initializing the metapopulation
    Deme* population = NULL;
    population = malloc(nDemes * sizeof(Deme));
    if(population == NULL){
        exit(0);
    }


    long listOfNtrlAlleles = MAX_NUMBER_OF_INITIAL_NTRL_ALLELES;
    long* pointeurSurListOfNtrlAlleles = &listOfNtrlAlleles;
    // vector containing the most recent alleles for each neutral locus
/*    long* listOfNtrlAlleles = NULL;
    listOfNtrlAlleles = malloc(nNtrlLoci * sizeof(long));
    if(listOfNtrlAlleles == NULL){
        exit(0);
    }
    for(i=0; i<nNtrlLoci; i++){
        listOfNtrlAlleles[i] = MAX_NUMBER_OF_INITIAL_NTRL_ALLELES;
    }
*/
    initializePopulation(r, population, nDemes, maxIndPerDem, nNtrlLoci, nQuantiLoci, fecundity);
    // Evolution of the metapopulation
    for(i=0; i<nGeneration; i++){    // start of the loop 'i' over the generations
        // DEBUG MODE
        printf("%d\n", nDemes);
        afficherPopulation(population, nDemes, nNtrlLoci, nQuantiLoci);
    
        int* nImmigrants = NULL; // stock the number of received immigrants per deme
        int* extinctionStatus = NULL;    // per deme: 0 = non-extincted; 1 = extincted
        int* nProducedSeeds = NULL; // stock the numer of produced seeds per deme
        Deme* newPopulation = NULL;
        nImmigrants = malloc(nDemes * sizeof(int));
        extinctionStatus = malloc(nDemes * sizeof(int));
        nProducedSeeds = malloc(nDemes * sizeof(int));
        newPopulation = malloc(nDemes * sizeof(Deme));


        if(nImmigrants == NULL || extinctionStatus == NULL || nProducedSeeds == NULL || newPopulation == NULL){
            exit(0);
        }
        setToZero(nDemes, nImmigrants, extinctionStatus, nProducedSeeds);    // initialize vectors used to configure the new-population
        
//    printf("nDemes: %d maxIndPerDem: %d migration: %.2lf extinction: %.2lf\n", nDemes, maxIndPerDem, migration, extinction);
        configMetapop(r, population, nDemes, maxIndPerDem, migration, extinction, nImmigrants, extinctionStatus, nProducedSeeds);    // get the parameters to create the new-population


        // DEBUG MODE
        printf("nImmigrants:\n");
        for(j=0; j<nDemes; j++){
            printf("%d ", nImmigrants[j]);
        }
        printf("\nextinctionStatus:\n");
        for(j=0; j<nDemes; j++){
            printf("%d ", extinctionStatus[j]);
        }
        printf("\nnProducedSeeds:\n");
        for(j=0; j<nDemes; j++){
            printf("%d ", nProducedSeeds[j]);
        }
        printf("\n");    


        initializeNewPopulation(newPopulation, nDemes, maxIndPerDem, nImmigrants, nProducedSeeds, extinctionStatus, nNtrlLoci, nQuantiLoci, recolonization, fecundity);


        panmixie(r, population, newPopulation, nDemes, maxIndPerDem, nProducedSeeds, nNtrlLoci, nQuantiLoci, pointeurSurListOfNtrlAlleles, ntrlMutation, quantiMutation, fecundity);


        // replace population by newPopulation
        // free memory
        libererMemoirePopulation(population, nDemes);
        free(population);
        Deme* population = NULL;
        population = malloc(nDemes * sizeof(Deme));
        if(population == NULL){
            exit(0);
        }
        
        replacement(r, population, newPopulation, nDemes, maxIndPerDem, nNtrlLoci, nQuantiLoci, fecundity, nImmigrants, nProducedSeeds, extinctionStatus, recolonization);


        free(nImmigrants);
        free(extinctionStatus);
        free(nProducedSeeds);
        libererMemoirePopulation(newPopulation, nDemes);
        free(newPopulation);


        // DEBUG MODE
        afficherPopulation(population, nDemes, nNtrlLoci, nQuantiLoci);
        printf("\n");
    }    // end of the loop 'i' over the generations
    printf("fin des simuls\n");    // DEBUG MODE
    //free(listOfNtrlAlleles);
    libererMemoirePopulation(population, nDemes);
    free(population);
    return(0);
}




void initializePopulation(gsl_rng* r, Deme* population, int nDemes, int maxIndPerDem, int nNtrlLoci, int nQuantiLoci, int fecundity){
    int i = 0;
    int j = 0;
    int k = 0;
//    unsigned int essais = 1;    // called by gsl_ran_binomial in order to determine whether we floor or ceil the expected number of offsprings
    int valuesNtrlAlleles = MAX_NUMBER_OF_INITIAL_NTRL_ALLELES;
    double minQuanti = 0.0;
    double maxQuanti = 0.0;
    maxQuanti = 1.0/2/nQuantiLoci;
    for(i=0; i<nDemes; i++){    // loop along the demes
        population[i].nIndividus = maxIndPerDem;
        population[i].ntrlLoci = malloc(2 * maxIndPerDem * nNtrlLoci * sizeof(long));
        population[i].quantiLoci = malloc(2 * maxIndPerDem * nQuantiLoci * sizeof(long));
        population[i].femaleAllocation = malloc(maxIndPerDem * sizeof(long));
        population[i].maleAllocation = malloc(maxIndPerDem * sizeof(long));
        population[i].nOffsprings = malloc(maxIndPerDem * sizeof(int));


        if(population[i].ntrlLoci == NULL || population[i].quantiLoci == NULL || population[i].femaleAllocation == NULL || population[i].maleAllocation == NULL || population[i].nOffsprings == NULL){
            exit(0);
        }


        for(j=0; j<maxIndPerDem; j++){    // loop along the individuals
            population[i].femaleAllocation[j] = 0.0;
            population[i].maleAllocation[j] = 1.0;
            for(k=0; k<(2*nNtrlLoci); k++){    // loop along the {2: diploid} x {nNtrlLoci: number of neutral loci} positions
                population[i].ntrlLoci[j*2*nNtrlLoci+k] = gsl_rng_uniform_int(r, valuesNtrlAlleles) + 1;
            }
            for(k=0; k<(2*nQuantiLoci); k++){    // loop along the {2: diploid} x {nQuantiLoci: number of quantitative loci} positions
                population[i].quantiLoci[j*2*nQuantiLoci+k] = gsl_ran_flat(r, minQuanti, maxQuanti);
                population[i].femaleAllocation[j] += population[i].quantiLoci[j*2*nQuantiLoci+k];
                population[i].maleAllocation[j] -= population[i].quantiLoci[j*2*nQuantiLoci+k];
            }
            population[i].nOffsprings[j] = floor(fecundity * population[i].femaleAllocation[j]) + gsl_ran_binomial(r, (fecundity * population[i].femaleAllocation[j]) - floor(fecundity * population[i].femaleAllocation[j]), 1);    // nOffs = floor(fecundity x femaleAllocation) + 1 according to a random Binomial integer
        }    // end of loop along the individuals
    }    // end of loop along the demes
}


void setToZero(int nDemes, int nImmigrants[], int extinctionStatus[], int nProducedSeeds[]){
    // function to set all to demes to zero for the #of immigrants, the extinction status and the #of produced seeds.
    // before configMetapop()
    int i = 0;
    for(i=0; i<nDemes; i++){
        nImmigrants[i] = 0;
        extinctionStatus[i] = 0;
        nProducedSeeds[i] = 0;
    }
}


void configMetapop(gsl_rng* r, Deme* population, int nDemes, const int maxIndPerDem, double migration, double extinction, int nImmigrants[], int extinctionStatus[], int nProducedSeeds[]){
    // function to get per deme the #of immigrants received, the extinction status and the #of produced seeds.
    // after setToZero()
    // modifies nImmigrants[]; extinctionStatus[] and nProducedSeeds[]
    int i = 0;
    int j = 0;
    const unsigned int binomTrials = 1;


    for(i=0; i<nDemes; i++){    // start of the loop over the demes
        nImmigrants[i] = gsl_ran_poisson(r, migration);
        extinctionStatus[i] = gsl_ran_binomial(r, extinction, binomTrials);    // 0: non-extincted; 1: extincted
        for(j=0; j<population[i].nIndividus; j++){    // start of the loop over the individuals
            //printf("deme: %d\tindividu: %d\tfemaleAllocation: %lf\tnOffspring: %d\n", i, j, population[i].femaleAllocation[j], population[i].nOffsprings[j]);
            nProducedSeeds[i] += population[i].nOffsprings[j];
        }    // end of the loop over the individuals
    }    // end of the loop over the demes
}


void initializeNewPopulation(Deme* newPopulation, const int nDemes, const int maxIndPerDem, const int nImmigrants[], const int nProducedSeeds[], int extinctionStatus[], const int nNtrlLoci, const int nQuantiLoci, const int recolonization, const int fecundity){
    // function to initialize the new population: allocate memory and set values to 0
    int i = 0;
    int j = 0;
    int k = 0;
    for(i=0; i<nDemes; i++){ // start of the loop along demes
        int taille = 0;
        taille = nProducedSeeds[i];    // size of the deme = nProducedSeeds within the parental population
        if(taille < fecundity){    // if not enough seeds are produced ==> deme is considered as extincted
            taille = fecundity;
//            extinctionStatus[i] = 1;
        }
        if(taille > maxIndPerDem){    // if too many individuals have to be present in the deme ==> cutoff to the maxIndPerDem (=carrying capacity)
            taille = maxIndPerDem;
        }
    
        newPopulation[i].nIndividus = taille;


        newPopulation[i].ntrlLoci = NULL;
        newPopulation[i].quantiLoci = NULL;
        newPopulation[i].femaleAllocation = NULL;
        newPopulation[i].maleAllocation = NULL;
        newPopulation[i].nOffsprings = NULL;


        newPopulation[i].ntrlLoci = malloc(2 * taille * nNtrlLoci * sizeof(long));
        newPopulation[i].quantiLoci = malloc(2 * taille * nQuantiLoci * sizeof(long));
        newPopulation[i].femaleAllocation = malloc(taille * sizeof(long));
        newPopulation[i].maleAllocation = malloc(taille * sizeof(long));
        newPopulation[i].nOffsprings = malloc(taille * sizeof(int));


        if(newPopulation[i].ntrlLoci == NULL || newPopulation[i].quantiLoci == NULL || newPopulation[i].femaleAllocation == NULL || newPopulation[i].maleAllocation == NULL || newPopulation[i].nOffsprings == NULL){
            exit(0);
        }


        for(j=0; j<taille; j++){ // start the loop along individuals
            newPopulation[i].femaleAllocation[j] = 0.0;
            newPopulation[i].maleAllocation[j] = 0.0;
            newPopulation[i].nOffsprings[j] = 0;
            for(k=0; k<(2*nNtrlLoci); k++){ // loop along the {2: diploid} x {nNtrlLoci: number of neutral loci} positions
                newPopulation[i].ntrlLoci[j*2*nNtrlLoci+k] = 0;
            }
            for(k=0; k<(2*nQuantiLoci); k++){       // loop along the {2: diploid} x {nQuantiLoci: number of quantitative loci} positions
                newPopulation[i].quantiLoci[j*2*nQuantiLoci+k] = 0;
            }
        }    // end of the loop along individuals
    } // end of the loop along demes
}


void panmixie(gsl_rng* r, Deme* population, Deme* newPopulation, const int nDemes, const int maxIndPerDem, const int nProducedSeeds[], const int nNtrlLoci, const int nQuantiLoci, long* pointeurSurListOfNtrlAlleles, double ntrlMutation, double quantiMutation, int fecundity){
    // function returning a new deme after a run of panmixia from the old deme 
    // here: only "deme specific" events are simulated (meiosis + mutation)
    int i = 0;
    int j = 0;
    int tmp = 0;
    int nMutation = 0;
        for(i=0; i<nDemes; i++){    // start the loop over the nDemes
            int K = 0;    // #of_individuals in the parental deme.
            int N = 0;    // #of_babies to produce
            double* parentalIndexes = NULL; // contain de the parental indexes (i.e:{0, 1, 2, 3, 4} if K == 5)
            double* mothers = NULL;    // array containing the mothers ID. Ex: {2, 19, 3, 3, 1} means that the first baby has individual#2 as mother, second baby has individual#19 as mother, babies 3 and 4 have individual#3
            double* fathers = NULL;    // array containing the fathers ID.


            K = population[i].nIndividus;    // #of_individuals in the parental deme
            N = newPopulation[i].nIndividus;     // #of_produced_seeds by the K parents


            parentalIndexes = malloc(K * sizeof(double));
            mothers = malloc(N * sizeof(double));    // indexes of mothers of the N autochtones within the deme of size newPopulation[i].nIndividus 
            fathers = malloc(N * sizeof(double));


            if(parentalIndexes == NULL || mothers == NULL || fathers == NULL){
                exit(0);
            }
            
            for(j=0; j<K; j++){
                parentalIndexes[j] = j;
            }
            for(j=0; j<N; j++){
                mothers[j] = 0;
                fathers[j] = 0;
            }


            // Sampling the parents            
            weightedSample(r, parentalIndexes, population[i].femaleAllocation, mothers, K, N); // return the indexes of N mothers (in 'mothers') that will be use to generate the N babies
            weightedSample(r, parentalIndexes, population[i].maleAllocation, fathers, K, N); // return the indexes of N fathers (in 'fathers') that will be use to generate the N babies


            // Meiosis and transmission of gametes
            int pos = 0;    // position in deme.ntrlLoci (or deme.quantiLoci) of offsprings
            int pos2 = 0;    // position in deme.ntrlLoci (or deme.quantiLoci) of parents
            for(j=0; j<N; j++){    // loop along the N sampled parents to fill the new deme with meiose (gsl_ran_binomial) 
                for(tmp=0; tmp<nNtrlLoci; tmp++){    // loop along the neutral loci
                    pos = j*nNtrlLoci*2 + tmp*2 + 0;    // position where to put the allele from the mother
                    pos2 = mothers[j]*nNtrlLoci*2 + tmp*2 + gsl_ran_binomial(r, 0.5, 1); // binomial transmission of the parental alleles with p=0.5 (=recombination at meiosis)
                    newPopulation[i].ntrlLoci[pos] = population[i].ntrlLoci[pos2];


                    pos += 1;    // position where to put the allele from the father = position of mother_allele + 1
                    pos2 = fathers[j]*nNtrlLoci*2 + tmp*2 + gsl_ran_binomial(r, 0.5, 1);
                    newPopulation[i].ntrlLoci[pos] = population[i].ntrlLoci[pos2];
                }    // end of loop along the neutral loci


                for(tmp=0; tmp<nQuantiLoci; tmp++){    // loop along the quantitative loci
                    pos = j*nQuantiLoci*2 + tmp*2 + 0;
                    pos2 = mothers[j]*nQuantiLoci*2 + tmp*2 + gsl_ran_binomial(r, 0.5, 1);
                    newPopulation[i].quantiLoci[pos] = population[i].quantiLoci[pos2];
                    
                    pos +=1;
                    pos2 = fathers[j]*nQuantiLoci*2 + tmp*2 + gsl_ran_binomial(r, 0.5, 1);
                    newPopulation[i].quantiLoci[pos] = population[i].quantiLoci[pos2];
                }    // end of loop along the quantitative loci
            }    // end of loop along the individuals


            // Put mutations
            nMutation = 0;
            nMutation = gsl_ran_binomial(r, ntrlMutation, 2*N*nNtrlLoci);
            if(nMutation > 0 ){
                for(j=0; j<nMutation; j++){
                    *pointeurSurListOfNtrlAlleles += 1;
                    pos = rand()%(2*N*nNtrlLoci);
                    newPopulation[i].ntrlLoci[pos] = *pointeurSurListOfNtrlAlleles;
                }
            }


            nMutation = 0;
            nMutation = gsl_ran_binomial(r, quantiMutation, 2*N*nQuantiLoci);
            if(nMutation > 0){
                for(j=0; j<nMutation; j++){
                    pos = rand()%(2*N*nQuantiLoci);
                    newPopulation[i].quantiLoci[pos] = gsl_ran_flat(r, 0.0, 1.0/2/nQuantiLoci);
                }
            }


            pos = 0;
            for(j=0; j<N; j++){    // loop along the individuals to calculate the new femaleAllocation and the number of offsprings
                tmp = 0;
                do{
                    newPopulation[i].femaleAllocation[j] += newPopulation[i].quantiLoci[pos];    // sum of allelic effects within individuals
                    pos +=1;
                    tmp +=1;
                }while(tmp<2*nQuantiLoci);
                newPopulation[i].maleAllocation[j] = 1 - newPopulation[i].femaleAllocation[j];    // set the male allocation as (1 - femaleAllocation)
                newPopulation[i].nOffsprings[j] = floor(fecundity * newPopulation[i].femaleAllocation[j]) + gsl_ran_binomial(r, (fecundity * newPopulation[i].femaleAllocation[j]) - floor(fecundity * newPopulation[i].femaleAllocation[j]), 1);    // set the #of_offsprings produced for each individual
            }


/*            // debug mode
            int pos3 = 0;
            for(j=0; j<N; j++){
                for(tmp=0; tmp<nNtrlLoci; tmp++){
                    pos = j*nNtrlLoci*2 + tmp*2 + 0;
                    pos2 = mothers[j]*nNtrlLoci*2 + tmp*2;
                    pos3 = fathers[j]*nNtrlLoci*2 + tmp*2;
                    printf("Deme.%d/Ind.%d/Locus.%d\tPos: %d/%d\tMother: %d\tFather: %d\tGenotypeMother: %ld/%ld\tGenotypeFather: %ld/%ld\tAllocF/M: %.1lf/%.1lf\tGenotypeSon: %ld/%ld\tAlloc: %lf\n", i, j, tmp, pos, pos+1, (int)mothers[j], (int)fathers[j], population[i].ntrlLoci[pos2], population[i].ntrlLoci[pos2+1], population[i].ntrlLoci[pos3], population[i].ntrlLoci[pos3+1], population[i].femaleAllocation[(int)mothers[j]], population[i].maleAllocation[(int)fathers[j]], newPopulation[i].ntrlLoci[(int)pos], newPopulation[i].ntrlLoci[pos+1], newPopulation[i].quantiLoci[j]);
                }
            }*/


            free(parentalIndexes);
            free(mothers);
            free(fathers);
    }    // end of loop over the nDemes
}


void replacement(gsl_rng* r, Deme* population, Deme* newPopulation, const int nDemes, const maxIndPerDem, const int nNtrlLoci, const int nQuantiLoci, int fecundity, const int nImmigrants[], int nProducedSeeds[], int extinctionStatus[], const int recolonization){
    int i = 0;
    int j = 0;
    int k = 0;
    int taille = 0;
    int* indexOfDemes = NULL;
    indexOfDemes = malloc(nDemes * sizeof(int));
    if(indexOfDemes == NULL){
        exit(0);
    }
    for(i=0; i<nDemes; i++){
        indexOfDemes[i] = i;
    }
    for(i=0; i<nDemes; i++){
        taille = newPopulation[i].nIndividus;// + nImmigrants[i];
        if(taille > maxIndPerDem){
            taille = maxIndPerDem;
        }
/*        if(extinctionStatus[i] == 1){
            taille = recolonization;
        }*/
        population[i].nIndividus = taille;
        population[i].ntrlLoci = malloc(2 * taille * nNtrlLoci * sizeof(long));
        population[i].quantiLoci = malloc(2 * taille * nQuantiLoci * sizeof(long));
        population[i].femaleAllocation = malloc(taille * sizeof(long));
        population[i].maleAllocation = malloc(taille * sizeof(long));
        population[i].nOffsprings = malloc(taille * sizeof(int));


        if(population[i].ntrlLoci == NULL || population[i].quantiLoci == NULL || population[i].femaleAllocation == NULL || population[i].maleAllocation == NULL || population[i].nOffsprings == NULL){
            exit(0);
        }
        
        for(j=0; j<(2 * population[i].nIndividus * nNtrlLoci); j++){
            population[i].ntrlLoci[j] = newPopulation[i].ntrlLoci[j];
        }
    
        for(j=0; j<(2 * population[i].nIndividus * nQuantiLoci); j++){
            population[i].quantiLoci[j] = newPopulation[i].quantiLoci[j];
        }


        for(j=0; j<population[i].nIndividus; j++){
            population[i].femaleAllocation[j] = newPopulation[i].femaleAllocation[j];
            population[i].maleAllocation[j] = newPopulation[i].maleAllocation[j];
            population[i].nOffsprings[j] = newPopulation[i].nOffsprings[j];
        }
    }    // end of loop over demes


    free(indexOfDemes);
}


void weightedSample(gsl_rng* r, const double* liste, const double* weights, double* target, const int sizeOfListe, const int nTrials){
    // function that fills the vector 'target' of size 'nTrials' containing the weighted-sampled 'sizeOfListe' elements of the 'liste':
    // weightedSample(gsl_rng* r, {2, 4, 6, 8, 10}, {1.2, 0.6, 0.3, 0.15, 0.05}, target, 5, 20)
    // target = {6, 6, 6, 2, 2, 2, 2, 2, 8, 2, 8, 6, 4, 2, 2, 4, 4, 4, 4, 2}
    // but can also be used for boolean sampling (pile ou face) using:
    // weightedSample(gsl_rng* r, {0, 1}, {1, 1}, target, 2, 1)
    int i = 0;
    int* n = NULL;
    int* sampledListe = NULL;
    n = malloc(sizeOfListe * sizeof(double));    // will contain the number of succes after K nTrials for each of the sizeOfListe elements of liste
    sampledListe = malloc(nTrials * sizeof(int));    // if n={0, 3, 1, 1}, sampledListe={1, 1, 1, 4, 5}


    gsl_ran_multinomial(r, sizeOfListe, nTrials, weights, n);    // return in 'n' the number of success for the sizeOfListe elements of liste


    int nValues = 0;
    int tmp = 0;
    int tmp2 = 0;
    for(i=0; i<sizeOfListe; i++){ // loop along the list called 'n' resulting from gsl_ran_multinomial
        nValues = n[i];
        if(nValues != 0){
            tmp2 = 0;
            do{
                sampledListe[tmp] = i;
                tmp++;
                tmp2++;
            }while(tmp2 < nValues);
        }
    }
    
    // shuffle values of the sampledListe
    gsl_permutation* p = gsl_permutation_alloc (nTrials);
    gsl_permutation_init (p);
    gsl_ran_shuffle(r, p -> data, nTrials, sizeof(size_t));


    tmp = 0;    
    for(i=0; i<nTrials; i++){
        tmp=gsl_permutation_get(p, i);
        target[i]=liste[sampledListe[tmp]];
    }
    gsl_permutation_free(p);
    free(n);
    free(sampledListe);
}


void libererMemoirePopulation(Deme* population, int nDemes){
    // free memory taken by the population at the end of each generation.
    int i = 0;
    for(i=0; i<nDemes; i++){
        free(population[i].ntrlLoci);
        free(population[i].quantiLoci);
        free(population[i].femaleAllocation);
        free(population[i].maleAllocation);
        free(population[i].nOffsprings);
    }
}


void afficherPopulation(Deme* population, int nDemes, int nNtrlLoci, int nQuantiLoci){
    // called to print some informations about population in a debug mode
    int i = 0;
    int j = 0;
    int k = 0;
    for(i=0; i<nDemes; i++){
        for(j=0; j<population[i].nIndividus; j++){
            printf("Deme: %d Ind: %d Ntrl: ", i, j);
            for(k=0; k<(2*nNtrlLoci); k++){
                printf("%ld ", population[i].ntrlLoci[2*j*nNtrlLoci+k]);
            }
            printf(" Quanti: ");
            for(k=0; k<(2*nQuantiLoci); k++){
                printf("%.2lf ", population[i].quantiLoci[2*j*nQuantiLoci+k]);
            }
        printf("femAlloc: %.2lf nOffs: %d\n", population[i].femaleAllocation[j], population[i].nOffsprings[j]);
        }
    }
}

Well the short answer is you compile for debug.

Code:

gcc -g *.c -L/usr/local/lib  -lgsl -lgslcblas -lm
gdb ./a.out

At the gdb prompt, you type

Code:

run 500 200 10000 20 0.00001 2 0.00001 10 0.001 0.02 2 123

When it crashes, the first thing to type is 'bt' to find the call stack of where you are in your code.

Other useful commands being 'up' and 'down' to navigate the call stack, and 'print' to print variables.

You can also set breakpoints to stop the code at any point (say where you think the bug might be about to happen), then single step the code to observe its behaviour in detail.

Hi Salem,
and thank you.

I think I know where is the problem of my code, I just don't understand why.

I have something similar to:

Code:

int i = 0;
for(i=0; i<1000; i++){
     printf("%d\n", population[0].nIndividuals);  // only works the first time
     ...
     printf("%d\n", population[0].nIndividuals); // also works the first time only
}

To me, there is no reason to lose "population" when throwing the loop second time since there is no function after the print. So, I probably miss a concept with loops in C

Well you could start by explaining what you think line 140 is meant to achieve.

> Deme* population = NULL;
Because now you have TWO variables called population in overlapping scopes.
One part of the code uses one of them, and another part of the code uses another.

I suggest you add flags such as
-Wall -Wextra -Wshadow
to your compiler command line, and fix everything it complains about before trying to run the code.

OK, I am starting to see.
The global idea was:
1) create an initial population before the loop over generations (variable called "population"), line 69
2) within the loop over generations, create variable called "newPopulation" (containing the next generation) by taking into account informations from "population" (former line 98). Until this step, everything is ok.
3) within the same loop, I want to replace "population" by "newPopulation". Like that, I can let my population evolving through a loop of length equal to nGeneration. To do that, I called:

free(population) // former ligne 139;

in order to release the area dedicated to the former "population" variable.

Then:

Deme* population = NULL; // former Line 140
population = malloc(nDemes * sizeof(Deme));

in order to re-create a new "population" variable that will contain informations stocked within "newPopulation"

Then, I fill "population".

I will try your flags (and will read a tuto about dbg).
Thanks again!

What i see in all your functions and also main, the variable 'population' is not a Deme, it is a pointer to Deme.
If you want to work with the members of Deme, you should dereference the pointer.

If a variable is of type Deme, you can work directly with members.

Code:

…
    Deme population[n];
…
    population[x].nIndividus = y;
…

But if the variable is a pointer to an type of Deme, you must dereference it.

Code:

…
    Deme *population;
    population = malloc(n * sizeof(*population));
…
    (*population)[x].nIndividus = y; // works, but looks ugly
    population[x]->nIndividus = y;    // looks nice
…

Now to your work at line 136 - 153.
It looks for me that you free all memory from population, copy the whole newPopulation to population (new malloc of the needed space) and than free newPopulation.
Both, population and newPupolation are pointers to the same type.
You can assign the pointer from newPopulation to population directly.
In that way, you don't must copy over all datas and free newPopulation, because population will become the newPopulation.
You should only set the pointer newPopulation to NULL.

Code:

        // replace population by newPopulation
        // free memory
        libererMemoirePopulation(population, nDemes);
        free(population);

        population = newPopulation;
        newPopulation = NULL;

        free(nImmigrants);
        free(extinctionStatus);
        free(nProducedSeeds);

His loop uses malloc to create a new block of population structs every loop, why do this if he can use malloc before the loop then realloc during the loop so the array continues to exist for as long as the loop lasts before finally freeing it at the end (though I will admit returning size in my example was bad).

Code:

typedef struct _Demes
{
  int count;
  int _count; // Total elements in allocated block
  size_t size;
  Deme *buff;
} Demes;
Deme* growDemes( Demes *population, int count );
void initDemes( Demes *population, int fromPos );
// main code
Demes population = { 1, 1, sizeof(Deme), malloc( sizeof( Deme ) ) };
if ( !population.buff )
  return failure_code;
for ( etc... )
{
  if ( growDemes( &population, nDemes ) == NULL )
    break;
  initDemes( &population, 0 );
}
free( population.buff );
memset( &population, 0, sizeof( Demes ) );
if ( loop_not_done )
  return failure_code;
// No extra cycles spent on creating whole buffer if not necessary

Like this he has option of expanding what he does with his population arrays beyond what he currently is doing.

Edit: I'm not saying to use the functions the exact same way as I have done here, I was abbreviating for the purpose of showing what I meant.

Originally Posted by awsdert

Code:

typedef struct _Demes
{
  int count;
  int _count; // Total elements in allocated block
  size_t size;
  Deme *buff;
} Demes;

The above is not good because the identifier _Demes is reserved to the implementation. The relevant rules are:

Originally Posted by C99 Clause 7.1.3 Paragraph 1 (part)

— All identifiers that begin with an underscore and either an uppercase letter or another underscore are always reserved for any use.
— All identifiers that begin with an underscore are always reserved for use as identifiers with file scope in both the ordinary and tag name spaces.

_count happens to be okay as an identifier because it is a member name and hence not at file scope and not in the ordinary or tag name spaces, but having members named count and _count for the same struct is likely to cause confusion.

If you want to use realloc, then besides the fix to return the pointer, you also should remember that realloc could return a null pointer.

what do you think the

== NULL

is for? Of course I'm aware realloc can return NULL, I allowed for that in growDemes return values, the _ was just lazy coding on my part, easily rectified with something like structDemes, _count would not cause confusion because it's obviously for private use when there is a more easily referenced member count (which a simple doxygen comment & intelligent IDE can inform a developer of anyway)

Originally Posted by awsdert

Code:

what do you think the
== NULL
is for? Of course I'm aware realloc can return NULL, I allowed for that in growDemes return values

Look at what you wrote earlier:

Code:

oldPopulation = realloc( oldPopulation, size );

The above makes the mistake of overwriting oldPopulation before checking for a null pointer (without a comment to indicate that you're aware of the problem and was just simplifying for a demo), so it is plausible to think that you made the same mistake again. Just because you return a null pointer from growDemes does not mean that you did not also overwrite the member. Furthermore, it is rather unusual that growDemes should return a pointer to Deme instead of a pointer to Demes: you're already providing an interface on the level of Demes, so why return a pointer to Deme?

Originally Posted by awsdert

the _ was just lazy coding on my part, easily rectified with something like structDemes

Agreed, but not everyone is aware of that such identifiers are reserved, so it is worthwhile to point it out. As an alternative, it would be simpler to just place the underscore at the end.

Originally Posted by awsdert

_count would not cause confusion because it's obviously for private use when there is a more easily referenced member count

No, that is not true. It depends on convention, and in fact if you're providing an interface of functions, all the members of the struct could well be "private", even if the opaque pointer idiom is not used.

Originally Posted by awsdert

(which a simple doxygen comment & intelligent IDE can inform a developer of anyway)

Comments can become outdated, perhaps by failing to be updated along with code. Hence, it is better to use good names than to comment poor ones.

Originally Posted by awsdert

as for returning a Deme pointer it was simply for a quick check on whether a buffer exists but can just as easily be an error code instead

Returning a pointer is workable; it is a matter of what pointer. Wrapping the dynamic array of Deme objects in another struct is a good idea: it creates an abstraction for a list of Deme objects, combining the array itself (via a pointer to its first element) and meta information concerning the array. growDemes is good: it creates an abstraction for "growing" a list of Deme objects. In other words, growDemes is part of the Demes interface. Therefore, why should it return a pointer to a Deme? If it should return a pointer, then it should return a pointer to the Demes object that was passed as an argument.

Originally Posted by awsdert

Code:

/// \brief struct that holds an array of Deme, the size of the buffer and a count of the expected elements, any member starting with _ is for private use

Why is there both count and _count to begin with? I might have written:

Code:

typedef struct Demes_
{
    size_t capacity; /* number of elements allocated */
    size_t size; /* number of elements in use */
    Deme *elements;
} Demes;

There is no need to do name decoration to denote a "private" member since I would provide an interface of functions for Demes. A programmer can easily see that this is a dynamic array struct, hence the elements are effectively read-only except by the use of the interface. If I wanted to make the members private, I would use the opaque pointer idiom.