Skip to content

Commit

Permalink
GSL-2.x support.
Browse files Browse the repository at this point in the history
  • Loading branch information
D. V. Wiebe committed Mar 10, 2016
1 parent 66e2719 commit a9d24f9
Show file tree
Hide file tree
Showing 2 changed files with 82 additions and 51 deletions.
67 changes: 41 additions & 26 deletions src/plugins/fits/non_linear.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -18,6 +18,7 @@
#include <gsl/gsl_blas.h>
#include <gsl/gsl_multifit_nlin.h>
#include <gsl/gsl_statistics.h>
#include <gsl/gsl_version.h>
#include "common.h"

struct data {
Expand Down Expand Up @@ -100,6 +101,7 @@ bool kstfit_nonlinear(
gsl_multifit_function_fdf function;
gsl_vector_view vectorViewInitial;
gsl_matrix* pMatrixCovariance;
gsl_matrix *pMatrixJacobian;
struct data d;
double dXInitial[NUM_PARAMS];
double* pInputX;
Expand Down Expand Up @@ -181,37 +183,50 @@ bool kstfit_nonlinear(
}
iIterations++;
} while( iStatus == GSL_CONTINUE && iIterations < MAX_NUM_ITERATIONS );
gsl_multifit_covar( pSolver->J, 0.0, pMatrixCovariance );

//
// determine the fitted values...
//
for( i=0; i<NUM_PARAMS; i++ ) {
dXInitial[i] = gsl_vector_get( pSolver->x, i );
}

for( i=0; i<iLength; i++ ) {
vectorOutYFitted->raw_V_ptr()[i] = function_calculate( pInputX[i], dXInitial );
vectorOutYResiduals->raw_V_ptr()[i] = pInputY[i] - vectorOutYFitted->raw_V_ptr()[i];
}
#if GSL_MAJOR_VERSION >= 2
pMatrixJacobian = gsl_matrix_alloc( iLength, NUM_PARAMS );
#else
pMatrixJacobian = pSolver->J;
#endif
if ( pMatrixJacobian != NULL) {
#if GSL_MAJOR_VERSION >= 2
gsl_multifit_fdfsolver_jac( pSolver, pMatrixJacobian );
#endif
gsl_multifit_covar( pMatrixJacobian, 0.0, pMatrixCovariance );

//
// determine the fitted values...
//
for( i=0; i<NUM_PARAMS; i++ ) {
dXInitial[i] = gsl_vector_get( pSolver->x, i );
}

//
// fill in the parameter values and covariance matrix...
//
for( i=0; i<NUM_PARAMS; i++ ) {
vectorOutYParameters->raw_V_ptr()[i] = gsl_vector_get( pSolver->x, i );
for( j=0; j<NUM_PARAMS; j++ ) {
vectorOutYCovariance->raw_V_ptr()[(i*NUM_PARAMS)+j] = gsl_matrix_get( pMatrixCovariance, i, j );
for( i=0; i<iLength; i++ ) {
vectorOutYFitted->raw_V_ptr()[i] = function_calculate( pInputX[i], dXInitial );
vectorOutYResiduals->raw_V_ptr()[i] = pInputY[i] - vectorOutYFitted->raw_V_ptr()[i];
}
}

//
// determine the value of chi^2/nu
//
scalarOutChi->setValue(gsl_blas_dnrm2( pSolver->f ));
//
// fill in the parameter values and covariance matrix...
//
for( i=0; i<NUM_PARAMS; i++ ) {
vectorOutYParameters->raw_V_ptr()[i] = gsl_vector_get( pSolver->x, i );
for( j=0; j<NUM_PARAMS; j++ ) {
vectorOutYCovariance->raw_V_ptr()[(i*NUM_PARAMS)+j] = gsl_matrix_get( pMatrixCovariance, i, j );
}
}

bReturn = true;
//
// determine the value of chi^2/nu
//
scalarOutChi->setValue(gsl_blas_dnrm2( pSolver->f ));

bReturn = true;

#if GSL_MAJOR_VERSION >= 2
gsl_matrix_free( pMatrixJacobian );
#endif
}
gsl_matrix_free( pMatrixCovariance );
}
gsl_multifit_fdfsolver_free( pSolver );
Expand Down
66 changes: 41 additions & 25 deletions src/plugins/fits/non_linear_weighted.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -18,6 +18,7 @@
#include <gsl/gsl_blas.h>
#include <gsl/gsl_multifit_nlin.h>
#include <gsl/gsl_statistics.h>
#include <gsl/gsl_version.h>
#include "common.h"

struct data {
Expand Down Expand Up @@ -101,6 +102,7 @@ bool kstfit_nonlinear_weighted(
gsl_multifit_function_fdf function;
gsl_vector_view vectorViewInitial;
gsl_matrix* pMatrixCovariance;
gsl_matrix *pMatrixJacobian;
struct data d;
double dXInitial[NUM_PARAMS];
double* pInputs[3];
Expand Down Expand Up @@ -197,37 +199,51 @@ bool kstfit_nonlinear_weighted(
}
while( iStatus == GSL_CONTINUE && iIterations < MAX_NUM_ITERATIONS );

gsl_multifit_covar( pSolver->J, 0.0, pMatrixCovariance );

//
// determine the fitted values...
//
for( i=0; i<NUM_PARAMS; i++ ) {
dXInitial[i] = gsl_vector_get( pSolver->x, i );
}
#if GSL_MAJOR_VERSION >= 2
pMatrixJacobian = gsl_matrix_alloc( iLength, NUM_PARAMS );
#else
pMatrixJacobian = pSolver-J;
#endif

if ( pMatrixJacobian != NULL) {
#if GSL_MAJOR_VERSION >= 2
gsl_multifit_fdfsolver_jac( pSolver, pMatrixJacobian );
#endif
gsl_multifit_covar( pMatrixJacobian, 0.0, pMatrixCovariance );

//
// determine the fitted values...
//
for( i=0; i<NUM_PARAMS; i++ ) {
dXInitial[i] = gsl_vector_get( pSolver->x, i );
}

for( i=0; i<iLength; i++ ) {
vectorOutYFitted->raw_V_ptr()[i] = function_calculate( pInputs[XVALUES][i], dXInitial );
vectorOutYResiduals->raw_V_ptr()[i] = pInputs[YVALUES][i] - vectorOutYFitted->raw_V_ptr()[i];
}
for( i=0; i<iLength; i++ ) {
vectorOutYFitted->raw_V_ptr()[i] = function_calculate( pInputs[XVALUES][i], dXInitial );
vectorOutYResiduals->raw_V_ptr()[i] = pInputs[YVALUES][i] - vectorOutYFitted->raw_V_ptr()[i];
}

//
// fill in the parameter values and covariance matrix...
//
for( i=0; i<NUM_PARAMS; i++ ) {
vectorOutYParameters->raw_V_ptr()[i] = gsl_vector_get( pSolver->x, i );
for( j=0; j<NUM_PARAMS; j++ ) {
vectorOutYCovariance->raw_V_ptr()[(i*NUM_PARAMS)+j] = gsl_matrix_get( pMatrixCovariance, i, j );
//
// fill in the parameter values and covariance matrix...
//
for( i=0; i<NUM_PARAMS; i++ ) {
vectorOutYParameters->raw_V_ptr()[i] = gsl_vector_get( pSolver->x, i );
for( j=0; j<NUM_PARAMS; j++ ) {
vectorOutYCovariance->raw_V_ptr()[(i*NUM_PARAMS)+j] = gsl_matrix_get( pMatrixCovariance, i, j );
}
}
}

//
// determine the value of chi^2/nu
//
scalarOutChi->setValue(gsl_blas_dnrm2( pSolver->f ));
//
// determine the value of chi^2/nu
//
scalarOutChi->setValue(gsl_blas_dnrm2( pSolver->f ));

bReturn = true;
bReturn = true;

#if GSL_MAJOR_VERSION >= 2
gsl_matrix_free( pMatrixJacobian );
#endif
}
gsl_matrix_free( pMatrixCovariance );
}
gsl_multifit_fdfsolver_free( pSolver );
Expand Down

0 comments on commit a9d24f9

Please sign in to comment.