fct_generateur_fichier.cpp

Aller à la documentation de ce fichier.

//####//------------------------------------------------------------
//####//------------------------------------------------------------
//####//  MAGiC
//####//  Jean Christophe Cuilliere et Vincent FRANCOIS
//####//  Departement de Genie Mecanique - UQTR
//####//------------------------------------------------------------
//####//  MAGIC est un projet de recherche de l equipe ERICCA
//####//  du departement de genie mecanique de l Universite du Quebec a Trois Rivieres
//####//  http://www.uqtr.ca/ericca
//####//  http://www.uqtr.ca/
//####//------------------------------------------------------------
//####//------------------------------------------------------------
//####//
//####//       fct_generateur_fichier.cpp
//####//
//####//------------------------------------------------------------
//####//------------------------------------------------------------
//####//    COPYRIGHT 2000-2024
//####//  jeu 13 jun 2024 11:58:52 EDT
//####//------------------------------------------------------------
//####//------------------------------------------------------------
 
 
#include "gestionversion.h"
 
#pragma hdrstop
 
 
#include <math.h>
#include "fct_generateur_fichier.h"
#include "mg_gestionnaire.h"
#include "ct_face.h"
#include "ct_point.h"
#include "stdlib.h"
#include <fstream>
 
 
FCT_GENERATEUR_FICHIER::FCT_GENERATEUR_FICHIER(std::string nomfich,double dn,double ind1,double ind2):FCT_GENERATEUR_ECHANTILLONAGE_ESTIMATION(dn,ind1,ind2),nomfichier(nomfich)
{
}
 
 
FCT_GENERATEUR_FICHIER::~FCT_GENERATEUR_FICHIER()
{
}
 
 
 
void FCT_GENERATEUR_FICHIER::remplace_virgule_avec_point(const char * nomfichier)
{
    FILE *in =fopen(nomfichier,"rt");
    if (in==NULL)
        return;
 
    std::string f2;
    char c;
    while (1)
    {
        c = fgetc(in);
        if (c != EOF)
        {
            if (c == ',') c='.';
            f2 += c;
        }
        else
            break;
    }
    fclose(in);
 
    std::ofstream file;
    file.open(nomfichier);
    file << f2;
    file.close();
}
 
 
void FCT_GENERATEUR_FICHIER::echantillonner(class MG_GROUPE_TOPOLOGIQUE* mggt)
{
    remplace_virgule_avec_point(nomfichier.c_str());
    FILE *in =fopen(nomfichier.c_str(),"rt");
    if (in!=NULL)
    {
        char message[500];
        char* res=fgets(message,255,in);
        sscanf(message,"%lf %lf %lf %lf %lf %lf ",&xmin,&ymin,&zmin,&xmax,&ymax,&zmax);
        res=fgets(message,255,in);
        sscanf(message,"%lf ",&dng);
        while (!feof(in))
        {
            res=fgets(message,255,in);
            double x,y,z,coef;
            double lambda=zonezoom;
            double courbure=0;
            int nb=sscanf(message,"%lf %lf %lf %lf %lf %lf ",&x,&y,&z,&coef,&lambda,&courbure);
            if (nb<=3) continue;
            double dist=dng*(1-coef)*lambda;
            BOITE_3D boite(x-dist,y-dist,z-dist,x+dist,y+dist,z+dist);
            double dn=coef*dng;
            DOUBLEN<1> ecart(&dn);
            CT_POINT<1>* pt=new CT_POINT<1>(x,y,z,ecart,boite,courbure);
            lstpoint.ajouter(pt);
 
        }
        fclose(in);
    }
}
 
 
 
 
void FCT_GENERATEUR_FICHIER::construit(int pas_echan,int pas_estima,MG_GROUPE_TOPOLOGIQUE* mggt)
{
    pas_echantillon=pas_echan;
    pas_estimation=pas_estima;
    echantillonner();
    estimation();
}
 
 
bool FCT_GENERATEUR_FICHIER::calibrage_fonction_newton(double eni, double *fonction, double &cr, double &er, double err, int type)
{
    double x, dx, y, z, fx, dc1, dc2, ds, emoins, eplus, er_crt, rm, Eprimmoins, Eprimplus, cmin, cmax = 0;
    int i, div, div1, div2, n_min =0;
    double minim[2][2];  //tableau contenu les deux minimums trouv�s, d�finins par leur c et E
 
    er_crt=err;
 
 
 
 
    double c[4];                 //tableau compos� de {ci,ci-1,cimin,cimin2,}
    double E[5];                 //tableau compos� de {Ei,Ei-1,Ei min,Ei min2,E associ� � Eprim min}
    double Eprim[4];             //tableau compos� de {Eprimi,Eprimi-1, Eprim min, Eprim min2}
    double Esec[2];             //tableau compos� de {Eseci,Eseci-1}
 
    double clarge[100];
    for (i=0;i<100;i++) clarge[i]=0;
    double elarge[100];
    for (i=0;i<100;i++) elarge[i]=0;
 
    dx=0.02;
 
 
 
 
 
 
    i = 0;
 
    dc1 = 5; //pr�cision de recherche des minimums
 
 
 
    for (c[0]=-50-2*dc1;c[0]<50;c[0]+=dc1)          //--------Boucle de recherche de c
    {
 
        eval_erreur(fonction, eni, c[0],er_crt, type, elarge[i], emoins, eplus,dx);
        if (i)
        {
            Eprim[0]=(elarge[i]-elarge[i-1])/dc1;
 
            if (Eprim[0]*Eprim[1]<0 && i>1)
            {
                clarge[n_min]=c[0]-dc1;
                clarge[n_min+1]=c[0];
                n_min+=2;
            }
            Eprim[1]=Eprim[0];
        }
        i++;
    }
 
    if (!n_min)
    {
        return false;
    }
 
    dc2 = 0.001; //pr�cision d'aproximation de la d�riv�e en une fonction du premier ordre
    for (i=0;i<4;i++)
    {
        c[i] = 0;
        E[i] = 0;
        if (i<2) Eprim[i] = 1;
        if (i<2) Esec[i] = 0;
    }
 
 
 
    c[0]=0;
    Eprim[2]=1;
    Eprim[3]=1;
    div=0;
    div1=0;
    div2=0;
    i=0;
 
    for (int nm=0; nm<n_min; nm++)
    {
        cmin=clarge[nm]-dc1;
        cmax=clarge[nm]+dc1;
        c[0]=clarge[nm];
 
        while (div < 2)
        {
 
 
            if (i)
            {
                c[1]=c[0];
                E[1]=E[0];
                Eprim[1]=Eprim[0];
                E[0]=0;
                Eprim[0]=0;
                Esec[1]=Esec[0];
 
                //------calcul de ci
                if (!div)
                {
                    if (div2>10)
                    {
                        c[0]=(E[3]*c[2]+E[2]*c[3])/(E[3]+E[2]);
                        div2++;
                    }
                    else  c[0] = -E[1]/Eprim[1] + c[1];
                }
                else
                {
                    if (div2>5)
                    {
                        c[0]=(Eprim[3]*c[2]+Eprim[2]*c[3])/(Eprim[3]+Eprim[2]);
                        div2=0;
                    }
                    else c[0] = -Eprim[1]/Esec[1] + c[1];
                }
 
                if (div1)
                {
                    c[0]=minim[0][0];
                    div1=0;
                }
            }
 
            if (c[0] > cmax) c[0]=cmax;
            if (c[0] < cmin) c[0]=cmin;
 
            eval_erreur(fonction, eni, c[0],er_crt, type, E[0], emoins, eplus,dx);
 
 
 
            Eprim[0] = (eplus - emoins) / (2 * dc2);
 
 
            Eprimmoins=(E[0]-emoins)/dc2;
            Eprimplus=(eplus-E[0])/dc2;
 
            Esec[0]=(Eprimplus-Eprimmoins)/dc2;
 
            if (!i)
            {
                c[2]=c[0];
                c[3]=c[0];
                E[2]=E[0];
                E[3]=E[0];
                Eprim[2]=Eprim[0];
                Eprim[3]=Eprim[0];
                E[4]=E[0];
            }
            else
            {
                if (!div)
                {
                    if (E[0] < E[2])
                    {
                        c[3]=c[2];
                        E[3]=E[2];
                        c[2]=c[0];
                        E[2]=E[0];
                    }
                    else if (E[0] < E[3] || i==1)
                    {
                        c[3]=c[0];
                        E[3]=E[0];
                    }
                    else div2++;
                }
                else
                {
                    if (fabs(Eprim[0]) < fabs(Eprim[2]))
                    {
                        c[3]=c[2];
                        Eprim[3]=Eprim[2];
                        c[2]=c[0];
                        Eprim[2]=Eprim[0];
                        E[4]=E[0];
                    }
 
                    else if (fabs(Eprim[0]) < fabs(Eprim[3]))
                    {
                        c[3]=c[0];
                        Eprim[3]=Eprim[0];
                    }
                    else div2++;
                }
            }
 
 
 
            if (div2>15)
            {
                minim[0][0]=c[2];
                minim[0][1]=E[2];
                div++;
                div1=1;
                div2=0;
            }
            if (div && (fabs(Eprim[0])<0.000000001 || E[2]<0.000000001 || i>150))
            {
                minim[1][0]=c[2];
                minim[1][1]=E[4];
                div++;
            }
 
            i++;
        }
 
        if (minim[0][1]<minim[1][1]) {
            clarge[nm]=minim[0][0];
            elarge[nm]=minim[0][1];
        }
        else {
            clarge[nm]=minim[1][0];
            elarge[nm]=minim[1][1];
        }
 
        if (!nm) {
            er=elarge[nm];
            cr=clarge[nm];
        }
        else if (elarge[nm]<er) {
            er=elarge[nm];
            cr=clarge[nm];
        }
 
    }
 
 
    if (type)
    {
        er=er*dx/(dng*dng);
    }
    else
    {
        er=er*dx/(2*dng*dng);
    }
 
 
    eval_erreur(fonction, eni, cr, er_crt, type, E[0],emoins, eplus,dx,1);  // stockage pog
    return true;
 
}
 
 
 
 
bool FCT_GENERATEUR_FICHIER::calibrage_fonction_rech(double eni, double *fonction, double &cr, double &er, double err, int type)
{
    double x, dx, y, z, dc,  er_crt, rm, emoins, eplus = 0;
    int i= 0;
    double minim[2][2];  //tableau contenu les deux minimums trouv�s, d�finins par leur c et E
 
    er_crt=err;
 
 
    double c[4];                 //tableau compos� de {ci,ci-1,cimin,cimin2,}
    double E[4];                 //tableau compos� de {Ei,Ei-1,Ei min,Ei min2,E associ� � Eprim min}
 
 
    dx=0.02;
 
 
    dc = 0.01; //pr�cision de recherche
 
    for (i=0;i<4;i++)
    {
        c[i] = 0;
        E[i] = 0;
    }
 
 
 
    c[0]=0;
    i = 0;
 
    for (c[0]=-15;c[0]<15;c[0]+=dc)          //--------Boucle de recherche de c
    {
 
 
        eval_erreur(fonction, eni, c[0],er_crt, type, E[0], emoins, eplus,dx);
 
        if (!i)
        {
            c[2]=c[0];
            c[3]=c[0];
            E[2]=E[0];
            E[3]=E[0];
        }
        else
        {
            if (E[0] < E[2])
            {
                c[3]=c[2];
                E[3]=E[2];
                c[2]=c[0];
                E[2]=E[0];
            }
            else if (E[0] < E[3])
            {
                c[3]=c[0];
                E[3]=E[0];
            }
        }
 
        i++;
    }
 
 
 
    cr=c[2];
    er=E[2];
    if (type)
    {
        er=er*dx/(dng*dng);
    }
    else
    {
        er=er*dx/(2*dng*dng);
    }
 
    eval_erreur(fonction, eni, cr, er_crt, type, E[0],emoins, eplus,dx,1);  // stockage pog
    return true;
 
}
 
void FCT_GENERATEUR_FICHIER::eval_erreur(double *fct, double eni, double courb, double er_c, int type, double &er, double &erm, double &erp,double dx, int stockage_fichier)
{
 
    double x=0;
    double y=0;
    double z=0;
    double e=0;
    double emoins=0;
    double eplus=0;
    double fx=0;
    double en=0;
    double enm=0;
    double enp=0;
    double de=0;
    double dem=0;
    double dep=0;
    double ds=0;
    double n=0;
    double a=0;
    double b=0;
    double dc;
    double dist;
    double cpd=0;
 
    int cp1=0;
    int cp2=0;
    int N=0;
    int nb_point=0;
 
 
 
    double delta=zonezoom * (dng - eni);
 
    dc = 0.001;  //--------pr�cision d'aproximation de la d�riv�e en une fonction du premier ordre
 
    CT_POINT<1>* pt;
 
 
    xmin=-1.125*delta;
    ymin=-1.125*delta;
    zmin=-1.125*delta;
    xmax=1.125*delta;
    ymax=1.125*delta;
    zmax=1.125*delta;
 
    BOITE_3D boite(xmin,ymin,zmin,xmax,ymax,zmax);
    grille.initialiser(boite.get_xmin(),boite.get_ymin(),boite.get_zmin(),boite.get_xmax(),boite.get_ymax(),boite.get_zmax(),20,20,20);
    BOITE_3D boitep(xmin,ymin,zmin,xmax,ymax,zmax);
 
 
 
    switch (type)
    {
 
    case 0:
    {
        //---- g�n�ration d'un seul point d'�chantillonnage si type=0
        dist=(dng-eni)*zonezoom;
        boitep.reinit(x-dist,y-dist,z-dist,x+dist,y+dist,z+dist);
        DOUBLEN<1> ecart(&eni);
        pt=new CT_POINT<1>(x,y,z,ecart,boitep);
        lstpoint.ajouter(pt);
        grille.inserer(pt);
 
    } break;
 
 
    case 1:
    {
 
 
        ds = distance_max_entre_points(eni, courb, er_c);
 
        n = std::max(1., floor(2 * delta / ds - 0.0001) + 1);
 
 
        N=floor(n/2);
        a=-N+1;
        b=N+1;
        if (N!=n/2) b++;
 
        for (cp1 = a;cp1<b;cp1++)
        {
            y = cp1 * ds - ds/2;
            dist=(dng-eni)*zonezoom;
            boitep.reinit(x-dist,y-dist,z-dist,x+dist,y+dist,z+dist);
            DOUBLEN<1> ecart(&eni);
            pt=new CT_POINT<1>(x,y,z,ecart,boitep);
            lstpoint.ajouter(pt);
            grille.inserer(pt);
        }
 
    }
    break;
 
    case 2:
    {
 
 
        ds = distance_max_entre_points(eni, courb, er_c);
 
        n = std::max(1., floor(2 * delta / ds - 0.0001) + 1);
 
 
        N=floor(n/2);
        a=-N+1;
        b=N+1;
        if (N!=n/2) b++;
 
        for (cp1 = a;cp1<b;cp1++)
        {
            for (cp2 = a;cp2<b;cp2++)
            {
                y = cp1 * ds - ds/2;
                z = cp2 * ds - ds/2;
                dist=(dng-eni)*zonezoom;
                boitep.reinit(x-dist,y-dist,z-dist,x+dist,y+dist,z+dist);
                DOUBLEN<1> ecart(&eni);
                pt=new CT_POINT<1>(x,y,z,ecart,boitep);
                lstpoint.ajouter(pt);
                grille.inserer(pt);
            }
        }
 
    }
    break;
 
    }
 
    if (stockage_fichier)
    {
 
        FILE* in =fopen(nomfichier.c_str(),"w");
        if (in!=NULL)
        {
            double coord[3];
 
            fprintf(in,"%lf %lf %lf %lf %lf %lf\n",xmin,ymin,zmin,xmax,ymax,zmax);
            fprintf(in,"%lf\n",dng);
            double coef= eni/dng;
 
            nb_point=lstpoint.get_nb();
 
            for (cp1=0;cp1<nb_point;cp1++)
            {
                CT_POINT<1> *pt=lstpoint.get(cp1);
                pt->evaluer(coord);
                fprintf(in,"%lf %lf %lf %lf %lf %lf\n",coord[0],coord[1],coord[2],coef,zonezoom,courb);
            }
 
        }
        fclose(in);
 
        return;
    }
 
    x=0;
    y=0;
    z=0;
 
 
    switch (type)
    {
 
    case 0:
    {
 
 
 
        for (x = 0;x<delta;x+=dx * delta)
        {
            fx = fonction_calibrage_theorique(x, eni, fct);
            en = estimer(x,0,0,courb);
            enm = estimer(x,0,0,courb-dc);
            enp = estimer(x,0,0,courb+dc);
 
            de = (fx - en);
            dem = (fx - enm);
            dep = (fx - enp);
 
            e+= de*de;
            emoins += dem*dem;
            eplus += dep*dep;
        }
 
    } break;
 
 
    case 1:
    {
 
 
        for (cpd=0;cpd<0.6;cpd+=0.5)
        {
            for (x = 0;x<delta;x+=dx * delta)
            {
                fx = fonction_calibrage_theorique(x, eni, fct);
                en = estimer(x,cpd*ds,0,courb);
                enm = estimer(x,cpd*ds,0,courb-dc);
                enp = estimer(x,cpd*ds,0,courb+dc);
 
                de = (fx - en);
                dem = (fx - enm);
                dep = (fx - enp);
 
                e+= de * de;
                emoins += dem * dem;
                eplus += dep * dep;
            }
        }
    }
    break;
 
    case 2:
    {
 
 
        for (cpd=0;cpd<0.6;cpd+=0.5)
        {
            for (x = 0;x<delta;x+=dx * delta)
            {
                fx = fonction_calibrage_theorique(x, eni, fct);
                en = estimer(x,cpd*ds,cpd*ds,courb);
                enm = estimer(x,cpd*ds,cpd*ds,courb-dc);
                enp = estimer(x,cpd*ds,cpd*ds,courb+dc);
 
                de = (fx - en);
                dem = (fx - enm);
                dep = (fx - enp);
 
                e+= de * de;
                emoins += dem * dem;
                eplus += dep * dep;
            }
        }
    }
    break;
 
 
    }
 
 
    while (lstpoint.get_nb())
    {
        CT_POINT<1> *pt=lstpoint.get(0);
        lstpoint.supprimer(pt);
        grille.supprimer(pt);
        delete pt;
    }
 
 
 
    er=e;
    erm=emoins;
    erp=eplus;
 
 
}
 
double FCT_GENERATEUR_FICHIER::distance_max_entre_points(double eni, double c,double err)
{
    double r=0;
    double d=0;
    double enmax=(1+err)*eni;
 
 
    if (c > 0.0001) r = 0.05 * (sqrt((zonetrans*(dng-enmax)-eni+enmax)*(zonetrans*(dng-enmax)-eni+enmax)-40*c*zonetrans*(eni-enmax)*(dng-enmax))-zonetrans*(dng-enmax)+eni-enmax)/(c*zonetrans*(dng-enmax));
    if (c<=0.0001) r = (enmax-eni)/(zonetrans*(dng-enmax)-eni+enmax);
 
 
    d = 2 * r * zonezoom * (dng -eni);
 
    return d;
 
}
 
 
double FCT_GENERATEUR_FICHIER::fonction_calibrage_theorique(double x, double eni, double *f)
{
    double fx=0;
    double ci=0;
    int i=0;
    int j=0;
    double delta=zonezoom * (dng - eni);
    if (f[0])
    {
        for (i=0;i<f[0]+1;i++)
        {
            ci=f[i+1];
            j=0;
            while (j<i)
            {
                ci*=x;
                j++;
            }
            fx+=ci;
        }
    }
    else    fx = eni + (dng - eni) * (x/delta);
 
 
    return fx;
 
}
 
 
 
 
 
 
#pragma package(smart_init)