← Back to branch summary

~showard314/ubuntu/natty/qtiplot/Python2.7_fix

« back to all changes in this revision

Viewing changes to fitPlugins/fitRational0/fitRational0.cpp

  • Committer: Bazaar Package Importer
  • Author(s): Gudjon I. Gudjonsson
  • Date: 2007-03-25 12:06:27 UTC
  • Revision ID: james.westby@ubuntu.com-20070325120627-5pvdufddr7i0r74x
Tags: upstream-0.9~rc2
Import upstream version 0.9~rc2

Show diffs side-by-side

added added

removed removed

Lines of Context:
fitPlugins/fitRational0/fitRational0.cpp
 
1
#include <math.h>
 
2
#include <stdlib.h>
 
3
#include <stdio.h>
 
4
#include <stddef.h>
 
5
 
 
6
#include <gsl/gsl_vector.h>
 
7
#include <gsl/gsl_blas.h>
 
8
 
 
9
#if defined(_MSC_VER) //MSVC Compiler
 
10
#define MY_EXPORT __declspec(dllexport)
 
11
#else
 
12
#define MY_EXPORT
 
13
#endif
 
14
 
 
15
struct data {
 
16
  size_t n;
 
17
  size_t p;
 
18
  double * X;
 
19
  double * Y;
 
20
  double * sigma;//weighting data
 
21
};
 
22
 
 
23
extern "C" MY_EXPORT char *name()
 
24
{
 
25
return "Rational0";
 
26
}
 
27
 
 
28
extern "C" MY_EXPORT char *function()
 
29
{
 
30
return "(b + c*x)/(1 + a*x)";
 
31
}
 
32
 
 
33
extern "C" MY_EXPORT char *parameters()
 
34
{
 
35
return "a,b,c";
 
36
}
 
37
 
 
38
extern "C" MY_EXPORT double function_eval(double x, double *params)
 
39
{
 
40
return (params[1] + x*params[2])/(1 + x*params[0]);
 
41
}
 
42
 
 
43
extern "C" MY_EXPORT int function_f (const gsl_vector * x, void *params,
 
44
        gsl_vector * f)
 
45
{
 
46
  size_t n = ((struct data *)params)->n;
 
47
  double *X = ((struct data *)params)->X;
 
48
  double *Y = ((struct data *)params)->Y;
 
49
  double *sigma = ((struct data *)params)->sigma;
 
50
 
 
51
  double a = gsl_vector_get (x, 0);
 
52
  double b = gsl_vector_get (x, 1);
 
53
  double c = gsl_vector_get (x, 2);
 
54
 
 
55
  size_t i;
 
56
  for (i = 0; i < n; i++)
 
57
    {
 
58
      double Yi = (b + c*X[i])/(1 + a*X[i]);
 
59
      gsl_vector_set (f, i, (Yi - Y[i])/sigma[i]);
 
60
    }
 
61
 
 
62
  return GSL_SUCCESS;
 
63
}
 
64
extern "C" MY_EXPORT double function_d (const gsl_vector * x, void *params)
 
65
{
 
66
  size_t n = ((struct data *)params)->n;
 
67
  double *X = ((struct data *)params)->X;
 
68
  double *Y = ((struct data *)params)->Y;
 
69
  double *sigma = ((struct data *)params)->sigma;
 
70
 
 
71
  double a = gsl_vector_get (x, 0);
 
72
  double b = gsl_vector_get (x, 1);
 
73
  double c = gsl_vector_get (x, 2);
 
74
 
 
75
  size_t i;
 
76
  double val=0;
 
77
  for (i = 0; i < n; i++)
 
78
    {
 
79
      double dYi = (((b + c*X[i])/(1 + a*X[i]))-Y[i])/sigma[i];
 
80
      val += dYi*dYi;
 
81
    }
 
82
 
 
83
  return val;
 
84
}
 
85
 
 
86
extern "C" MY_EXPORT int function_df (const gsl_vector * x, void *params,
 
87
         gsl_matrix * J)
 
88
{
 
89
  size_t n = ((struct data *)params)->n;
 
90
  double *X = ((struct data *)params)->X;
 
91
  double *sigma = ((struct data *)params)->sigma;
 
92
 
 
93
  double a = gsl_vector_get (x, 0);
 
94
  double b = gsl_vector_get (x, 1);
 
95
  double c = gsl_vector_get (x, 2);
 
96
 
 
97
  size_t i;
 
98
  for (i = 0; i < n; i++)
 
99
    {
 
100
      /* Jacobian matrix J(i,j) = dfi / dxj, 
 
101
         where fi = (Yi - Y[i])/sigma[i],                                   
 
102
         Yi = (b + c*X[i])/(1 + a*X[i])         
 
103
         and the xj are the parameters (a, b, c) */
 
104
 
 
105
          double s = sigma[i];
 
106
      double t = X[i];
 
107
      double e = 1/(1 + a*t);
 
108
      gsl_matrix_set (J, i, 0, -e*e*a*(b+c*t)/s);
 
109
      gsl_matrix_set (J, i, 1, e/s);
 
110
      gsl_matrix_set (J, i, 2, e*t/s);
 
111
 
 
112
    }
 
113
  return GSL_SUCCESS;
 
114
}
 
115
 
 
116
extern "C" MY_EXPORT int function_fdf (const gsl_vector * x, void *params,
 
117
          gsl_vector * f, gsl_matrix * J)
 
118
{
 
119
  function_f (x, params, f);
 
120
  function_df (x, params, J);
 
121
 
 
122
  return GSL_SUCCESS;
 
123
}