Le wiki de Syracuse

Outils d'utilisateurs

Outils du Site

Piste:
A PCRE internal error occured. This might be caused by a faulty plugin
mc:moindre.pl

====== Différences ====== Cette page vous donne les différences entre la révision choisie et la version actuelle de la page.

Lien vers cette vue

mc:moindre.pl [2008/08/05 20:06]
maxime 		mc:moindre.pl [2010/02/12 02:22] (Version actuelle)
newacct
Ligne 64: Ligne 64:
 [[http://melusine.eu.org/syracuse/journal/2008/07/21/perl/]] [[http://melusine.eu.org/syracuse/journal/2008/07/21/perl/]]
  
 +==== Le code ====
 +
 +<code perl>
 +#!/usr/bin/perl
 +# -----------------------------------------------------------------------------
 +# MOINDRE
 +# Maxime Chupin
 +# -----------------------------------------------------------------------------
 +# version 1.0
 +
 +use Math::Matrix;
 +use Cwd;
 +use Getopt::Std;        # Usage du module standard d'acquisition des options
 +                        # sur la ligne de commande.
 +
 +#---------------- récupération du chemin --------------------------------------
 +my $chemin = cwd();
 +#------------------------------------------------------------------------------
 +
 +#------------------- les options ----------------------------------------------
 +getopts("d:");
 +
 +my $defaut_d = 3;
 +
 +$opt_d and $defaut_d = $opt_d;
 +print ("$defaut_d\n");
 +#------------------------------------------------------------------------------
 +
 +#---------------- subroutine d'un matrice unicolonne --------------------------
 +sub MatCol {
 +    my @list = @_;
 +    my $l = [];
 +    foreach $v (@list){
 + push  @$l, $v;
 +    }
 +    my $mat = new Math::Matrix($l);
 +    return $mat->transpose;
 +}
 +#------------------------------------------------------------------------------
 +
 +#---------------------- listes des points mesures -----------------------------
 +my @xi = ();
 +my @yi = ();
 +#------------------------------------------------------------------------------
 +
 +#---------- récupération des listes de points d'un fichiers externe -----------
 +# Jean-Michel Sarlat, cette opération trie les listes par ordre croissant de x
 +
 +@p = (); # Table de références à des points.
 +
 +# Ouverture en lecture seule du fichier xy.dat
 +open(DAT,"$chemin/$ARGV[0]") or die "Impossible d'ouvrir $ARGV[0]... : $!";
 +while (<DAT>) { # Lecture ligne par ligne du fichier (handle DAT)
 +    chomp; # Suppresion des blancs finaux, y compris le retour à la ligne
 +    s/^\s*//; # Suppression des blancs initiaux (remplacement des blancs par rien...)
 +    if($_){             # Traite la suite si la ligne est non vide
 +        my $r = {}; # Création d'une référence vide...
 +        ($r->{x},$r->{y}) = split /\s+/;  # Découpage de la ligne, le séparateur est 
 +      # une succession de blancs.
 +    push @p, $r;} # Ajout de $r à la table des points.
 +}
 +# Tri de la liste des points selon x.
 +@p = sort { $a->{x} <=> $b->{x} } @p;
 +# Construction des tables des x et y.
 +foreach (@p) {
 +    push @xi, $_->{x};
 +    push @yi, $_->{y};
 +}
 +#------------------------------------------------------------------------------
 +
 +#---------------------- min Max @xi et pas ------------------------------------
 +my $minxi;
 +my $maxxi;
 +$minxi = $xi[0];
 +$maxxi = $xi[-1];
 +$pas = abs($minxi-$maxxi)/50;
 +#------------------------------------------------------------------------------
 +
 +#----------------------------- Les définitions des "outils" -------------------
 +my $taille = @xi; #comme son nom l'indique
 +my @S=(); # tableau des sommes des xi^k, k de 0 à 2*$ordre
 +my @W=(); # tableau des sommes des yi$xi^k, k de 0 à $ordre
 +my $ordre;
 +($prefixe = $ARGV[0]) =~ s/\.dat$//; #récupération du préfixe
 +if("$ARGV[1]" eq "exp"){
 +    $ordre = 1;
 +}
 +else{
 +    $ordre = $ARGV[1];  #ordre de regression
 +}
 +
 +my $matrice = new Math::Matrix([1]) ;# définition de la matrice des Sk
 +#------------------------------------------------------------------------------
 +
 +#-------------------- regression exponentielle --------------------------------
 +my @saveyi = @yi;
 +if("$ARGV[1]" eq "exp"){
 +    for my $i(0..($taille-1)){
 + $yi[$i] = log($yi[$i]);
 +    }
 +}
 +
 +#------------------------ Calculs des S_k et W_k ------------------------------
 +for my $k ( 0 .. (2*$ordre) )
 +  {
 +    for my $i ( 0 .. ($taille-1) )
 +      {
 +   if($k==0){
 +       $xik=1;
 +   }
 +   else{
 +       $xik=1;
 +       $pro=$xi[$i];
 +       for my $l (1..$k){
 +   $xik*=$pro;
 +       }
 +   }
 +   $S[$k]+=$xik;
 +   if($k<=$ordre){
 +       $W[$k]+=$yi[$i]*$xik;
 +   }
 +      } 
 +  }
 +
 +#------------------------------------------------------------------------------
 +
 +$matW = MatCol(@W);
 +
 +for my $i (0 .. $ordre){
 +    for my $j (0 .. $ordre){
 + $matrice->[$i]->[$j] = $S[$i+$j];
 +    }
 +}
 +$sys = $matrice->concat($matW);
 +$sol=$sys->solve;
 +
 +#------------------------------------------------------------------------------
 +
 +#------------------- liste des coordonnées de la régression -------------------
 +my @xr=();
 +my @yr=();
 +my $incr=0;
 +
 +open(DATTH, ">$chemin/regr$ARGV[0]");
 +for(my $i=$minxi; $i<=$maxxi+$pas; $i+=$pas){
 +    my $sum=0;
 +    $xr[$incr]=$i;
 +#si regression exponentielle y=exp(ax+b)
 +    if("$ARGV[1]" eq "exp"){
 + $yr[$incr] = exp(($sol->[0]->[0]))*(exp(($sol->[1]->[0])*$xr[$incr]));
 +    }
 +#regression "normal"
 +    else{
 + for my $j (0 .. $ordre+1){
 +     $sum+=($sol->[$j]->[0])*(($xr[$incr]**($j)));
 + }
 + $yr[$incr]=$sum;
 +    }
 +    print (DATTH "$xr[$incr] $yr[$incr]\n");
 +    $incr++;
 +}
 +close(DATTH);
 +#------------------------------------------------------------------------------
 +
 +#------------------------ on imprime les coefficients -------------------------
 +open(FIC,">$chemin/coef$ARGV[0]");
 +
 +print(FIC "Régression d'ordre $ARGV[1]\n $sol\n");
 +close(FIC);
 +#------------------------------------------------------------------------------
 +#----------------------- On crée le fichier metapost --------------------------
 +my $legende;
 +my $coeffexp;
 +open(FICMP, ">$chemin/MetaRg$prefixe.mp");
 +# Définition de la légende "regression"
 +
 +# Définition de la légende "regression"
 +# JMS - Retouche pour montrer l'utilisation l'instruction sprintf. Elle vient du
 +# C et est très utile pour le formatage des résultats !
 +if("$ARGV[1]" eq "exp"){
 +    $legende = sprintf("%0.${defaut_d}f \\text{e}^{%0.${defaut_d}f x}",
 + exp($sol->[0]->[0]),exp($sol->[1]->[0]));
 +}
 +else {
 +    $legende = sprintf("%0.${defaut_d}f",$sol->[0]->[0]);
 +    for my $i (1 .. $ordre){
 + $legende .= sprintf("%+0.${defaut_d}f x",$sol->[$i]->[0]);
 + $i > 1 and $legende .= "^{$i}";
 +    }
 +}
 +
 +
 +# Calcul des coordonnées de "fixation" de la légende
 +
 +my $xleg = $minxi+1*($maxxi-$minxi)/2;
 +# on trouve le $yi max
 +my $yleg;
 +if("$ARGV[1]" eq "exp"){
 +    @yi = @saveyi;
 +}
 +for(my $i=0; $i<=$taille-1; $i+=1){
 +    if($i==0){
 + $yleg=$yi[$i];
 +    }
 +    else{
 + if($yi[$i]>$yleg){
 +     $yleg=$yi[$i];
 + }
 +    }
 +}
 +
 +print ("$xi[-1] $xleg $yleg \n");
 +print FICMP << "ENDOFMP";
 +verbatimtex
 +%&latex
 +\\documentclass{article}
 +\\usepackage[latin1]{inputenc}
 +\\usepackage{amsmath}
 +\\usepackage[frenchb]{babel}
 +\\begin{document}
 +etex
 + 
 +picture bullet;
 +bullet = image(draw origin withpen pencircle scaled 1.5mm withcolor green);
 + 
 +input graph;
 +
 +beginfig(1);
 +draw begingraph(17cm,12cm);
 +glabel.lft(btex \\Large $ARGV[2] etex rotated 90,OUT);
 +glabel.bot(btex \\Large $ARGV[3] etex,OUT);
 +gdraw "$ARGV[0]" plot bullet ;
 +gdraw "regr$ARGV[0]" withcolor red withpen pencircle scaled 1.3pt;
 +autogrid(grid.bot,grid.lft) withcolor .85white ;
 +setrange(whatever,whatever,whatever,whatever);
 +picture legende;
 +legende = btex Régression  : \$\\boxed{$legende}\$  etex;
 +glabel(image(unfill bbox legende; draw legende) , ($xleg,$yleg));
 +frame;
 +endgraph;
 +endfig;
 +end.
 +ENDOFMP
 +close(FICMP);
 +#------------------------------------------------------------------------------
 +#
 +</code>
mc/moindre.pl.1217959579.txt.gz · Dernière modification: 2008/08/05 20:06 par maxime

Outils de la Page

Sauf mention contraire, le contenu de ce wiki est placé sous la licence suivante : GNU Free Documentation License 1.3
GNU Free Documentation License 1.3 Donate Powered by PHP Valid HTML5 Valid CSS Driven by DokuWiki