From a17f2a7da2ec0056fcb9bdfdcb18ad506f9217ab Mon Sep 17 00:00:00 2001 From: Adrian Przybylski Date: Wed, 17 Jul 2019 12:15:12 +0200 Subject: [PATCH 1/9] bugfix initialization of n_iterations_ --- Gpufit/gpu_data.cu | 1 + 1 file changed, 1 insertion(+) diff --git a/Gpufit/gpu_data.cu b/Gpufit/gpu_data.cu index 2ad00bc1..e27750c3 100644 --- a/Gpufit/gpu_data.cu +++ b/Gpufit/gpu_data.cu @@ -135,6 +135,7 @@ void GPUData::init set(scaling_vectors_, 0., chunk_size_ * info_.n_parameters_to_fit_); set(states_, 0, chunk_size_); set(lambdas_, 0.001f, chunk_size_); + set(n_iterations_, 0, chunk_size_); } void GPUData::init_user_info(char const * const user_info) From 0ab8becab0c1716b3e43967e776a1adf58074d02 Mon Sep 17 00:00:00 2001 From: Scinfaxi Date: Mon, 19 Aug 2019 18:21:54 -0700 Subject: [PATCH 2/9] CUDA code and .cpp testing added. --- Gpufit/constants.h | 3 +- Gpufit/examples/Simple_Example.cpp | 268 +++++++++++++++++++++++------ Gpufit/models/models.cuh | 7 +- Gpufit/models/patlak.cuh | 55 ++++++ 4 files changed, 276 insertions(+), 57 deletions(-) create mode 100644 Gpufit/models/patlak.cuh diff --git a/Gpufit/constants.h b/Gpufit/constants.h index 62964ee8..72fe2657 100644 --- a/Gpufit/constants.h +++ b/Gpufit/constants.h @@ -11,7 +11,8 @@ enum ModelID { CAUCHY_2D_ELLIPTIC = 4, LINEAR_1D = 5, FLETCHER_POWELL_HELIX = 6, - BROWN_DENNIS = 7 + BROWN_DENNIS = 7, + PATLAK = 12 }; // estimator ID diff --git a/Gpufit/examples/Simple_Example.cpp b/Gpufit/examples/Simple_Example.cpp index 331c0899..ac9cf091 100644 --- a/Gpufit/examples/Simple_Example.cpp +++ b/Gpufit/examples/Simple_Example.cpp @@ -1,51 +1,129 @@ #include "../gpufit.h" -#include + #include +#include +#include +#include +//using namespace std; -void simple_example() +void patlak_two() { + /* - This example demonstrates a simple, minimal program containing all - of the required parameters for a call to the Gpufit function. The example - can be built and executed within the project environment. Please note that - this code does not actually do anything other than make a single call to - gpufit(). - - In the first section of the code, the *model ID* is set, memory space for - initial parameters and data values is allocated, the *fit tolerance* is set, - the *maximum number of iterations* is set, the *estimator ID* is set, and - the *parameters to fit array* is initialized. Note that in most applications, - the data array will already exist and it will be unnecessary to allocate - additional space for data. In this example, the *parameters to fit* array - is initialized to all ones, indicating that all model parameters should be - adjusted in the fit. + This example generates test data in form of 10000 one dimensional linear + curves with the size of 20 data points per curve. It is noised by normal + distributed noise. The initial guesses were randomized, within a specified + range of the true value. The LINEAR_1D model is fitted to the test data sets + using the LSE estimator. The optional parameter user_info is used to pass + custom x positions of the data sets. The same x position values are used for + every fit. + + The console output shows + - the ratio of converged fits including ratios of not converged fits for + different reasons, + - the values of the true parameters and the mean values of the fitted + parameters including their standard deviation, + - the mean chi square value + - and the mean number of iterations needed. */ - /*************** definition of input and output parameters ***************/ + // number of fits, fit points and parameters + size_t const n_fits = 100; // 10000? 100? + size_t const n_points_per_fit = 60; + size_t const n_model_parameters = 2; + + // custom x positions for the data points of every fit, stored in user info + REAL time[] = { 0.25, 0.5, 0.75, 1, 1.25, 1.5, 1.75, 2, 2.25, 2.5, 2.75, 3, 3.25, 3.5, 3.75, 4, 4.25, 4.5, 4.75, 5, + 5.25, 5.5, 5.75, 6, 6.25, 6.5, 6.75, 7, 7.25, 7.5, 7.75, 8, 8.25, 8.5, 8.75, 9, 9.25, 9.5, 9.75, 10, + 10.25, 10.5, 10.75, 11, 11.25, 11.5, 11.75, 12, 12.25, 12.5, 12.75, 13, 13.25, 13.5, 13.75, 14, 14.25, 14.5, 14.75, 15 }; // in min. up to 15 min + + REAL Cp[] = { 0.0f, 0.0f, 0.0f, 3.01846399851715f, 2.59801604007558f, 2.2361331285733f, 1.92465762011135f, 1.65656816551711f, 1.4258214335524f, + 1.22721588081636f, 1.05627449741415f, 0.909143885218726f, 0.782507393725825f, 0.6735103553914f, 0.579695735090254f, + 0.498948743091769f, 0.429449163006342f, 0.369630320068624f, 0.318143764811612f, 0.273828876023252f, 0.235686697768721f, + 0.20285742070682f, 0.174601000079374f, 0.150280473460109f, 0.12934760220805f, 0.111330512951924f, 0.0958230605172143f, + 0.0824756725126274f, 0.0709874691926393f, 0.0610994809603327f, 0.0525888106179893f, 0.0452636087696102f, 0.0389587491097867f, + 0.033532106110336f, 0.0288613511954976f, 0.0248411951843697f, 0.0213810148391187f, 0.018402810016092f, 0.0158394453694853f, + 0.013633136971665f, 0.0117341497352074f, 0.010099676273659f, 0.00869287192804919f, 0.00748202420651342f, 0.00643983791435146f, + 0.00554281985976681f, 0.00477074926518851f, 0.00410622194607174f, 0.00353425798195557f, 0.00304196403581309f, 0.00261824270962248f, + 0.00225354238438883f, 0.00193964190545541f, 0.00166946525943377f, 0.00143692206515917f, 0.00123677028298367f, 0.00106449804756952f, + 0.000916221960431984f, 0.000788599549519612f, 0.000678753922476738f }; - // number of fits, number of points per fit - size_t const n_fits = 10; - size_t const n_points_per_fit = 10; - // model ID and number of model parameters - int const model_id = GAUSS_1D; - size_t const n_model_parameters = 4; + std::vector< REAL > user_info(2 * n_points_per_fit); + for (size_t i = 0; i < n_points_per_fit; i++) + { + user_info[i] = static_cast(time[i]); + std::cout << user_info[i] << std::endl; + } + std::cout << "Time filled..." << std::endl; + + for (size_t i = n_points_per_fit; i < 2 * n_points_per_fit; i++) + { + user_info[i] = static_cast(Cp[i - n_points_per_fit]); + std::cout << user_info[i] << std::endl; + } + std::cout << "Checking user_info vector..." << std::endl; + for (size_t i = 0; i < 2 * n_points_per_fit; i++) + { + std::cout << user_info[i] << std::endl; + } - // initial parameters + std::cout << user_info.size() << std::endl; + std::cout << user_info.max_size() << std::endl; + + // size of user info in bytes + size_t const user_info_size = 2 * n_points_per_fit * sizeof(REAL); + + // initialize random number generator + std::mt19937 rng; + rng.seed(0); + std::uniform_real_distribution< REAL > uniform_dist(0, 1); + std::normal_distribution< REAL > normal_dist(0, 1); + + // true parameters + std::vector< REAL > true_parameters{ 0.005, 0.02 }; // Ktrans, vp -- ? + + // initial parameters (randomized) std::vector< REAL > initial_parameters(n_fits * n_model_parameters); + for (size_t i = 0; i != n_fits; i++) + { + // random offset + initial_parameters[i * n_model_parameters + 0] = true_parameters[0] * (0.8f + 0.4f * uniform_dist(rng)); + // random slope + initial_parameters[i * n_model_parameters + 1] = true_parameters[0] * (0.8f + 0.4f * uniform_dist(rng)); + } - // data + std::cout << std::endl << "parameters randomized..." << std::endl; + + // generate data -- NEEDS TO BE UPDATED WITH PATLAK MODEL std::vector< REAL > data(n_points_per_fit * n_fits); + for (size_t i = 0; i != data.size(); i++) + { + size_t j = i / n_points_per_fit; // the fit + size_t k = i % n_points_per_fit; // the position within a fit + REAL x = 0; + for (int n = 1; n < k; n++) { // or is point_index i's limit? ya probably + // + REAL spacing = time[n] - time[n - 1]; + x += Cp[n] * spacing + 0.5 * (Cp[n] - Cp[n - 1]); + } + //REAL x = user_info[k]; // how to define w/ 2 x's? integral of Cp(k) to point_index? + REAL y = true_parameters[0] * x + true_parameters[1] * user_info[k + n_points_per_fit]; + data[i] = y + normal_dist(rng); + } // tolerance - REAL const tolerance = 0.001f; + REAL const tolerance = 0.001f; // calculated by Christian to be 10e-15f // maximum number of iterations - int const max_number_iterations = 10; + int const max_number_iterations = 200; // estimator ID int const estimator_id = LSE; + // model ID + int const model_id = PATLAK; + // parameters to fit (all of them) std::vector< int > parameters_to_fit(n_model_parameters, 1); @@ -55,44 +133,124 @@ void simple_example() std::vector< REAL > output_chi_square(n_fits); std::vector< int > output_number_iterations(n_fits); - /***************************** call to gpufit ****************************/ - + // call to gpufit (C interface) int const status = gpufit - ( - n_fits, - n_points_per_fit, - data.data(), - 0, - model_id, - initial_parameters.data(), - tolerance, - max_number_iterations, - parameters_to_fit.data(), - estimator_id, - 0, - 0, - output_parameters.data(), - output_states.data(), - output_chi_square.data(), - output_number_iterations.data() - ); - - /****************************** status check *****************************/ + ( + n_fits, + n_points_per_fit, + data.data(), + 0, + model_id, + initial_parameters.data(), + tolerance, + max_number_iterations, + parameters_to_fit.data(), + estimator_id, + user_info_size, + reinterpret_cast< char* >( user_info.data() ), + output_parameters.data(), + output_states.data(), + output_chi_square.data(), + output_number_iterations.data() + ); + + std::cout << std::endl << "gpufit complete." << std::endl; + + // check status if (status != ReturnState::OK) { throw std::runtime_error(gpufit_get_last_error()); } + + + std::cout << std::endl << "status is ok." << std::endl; + + + // get fit states + std::vector< int > output_states_histogram(5, 0); + for (std::vector< int >::iterator it = output_states.begin(); it != output_states.end(); ++it) + { + output_states_histogram[*it]++; + } + + std::cout << "ratio converged " << (REAL)output_states_histogram[0] / n_fits << "\n"; + std::cout << "ratio max iteration exceeded " << (REAL)output_states_histogram[1] / n_fits << "\n"; + std::cout << "ratio singular hessian " << (REAL)output_states_histogram[2] / n_fits << "\n"; + std::cout << "ratio neg curvature MLE " << (REAL)output_states_histogram[3] / n_fits << "\n"; + std::cout << "ratio gpu not read " << (REAL)output_states_histogram[4] / n_fits << "\n"; + + // compute mean fitted parameters for converged fits + std::vector< REAL > output_parameters_mean(n_model_parameters, 0); + for (size_t i = 0; i != n_fits; i++) + { + if (output_states[i] == FitState::CONVERGED) + { + // add Ktrans + output_parameters_mean[0] += output_parameters[i * n_model_parameters + 0]; + // add vp + output_parameters_mean[1] += output_parameters[i * n_model_parameters + 1]; + } + } + output_parameters_mean[0] /= output_states_histogram[0]; + output_parameters_mean[1] /= output_states_histogram[0]; + + // compute std of fitted parameters for converged fits + std::vector< REAL > output_parameters_std(n_model_parameters, 0); + for (size_t i = 0; i != n_fits; i++) + { + if (output_states[i] == FitState::CONVERGED) + { + // add squared deviation for Ktrans + output_parameters_std[0] += (output_parameters[i * n_model_parameters + 0] - output_parameters_mean[0]) * (output_parameters[i * n_model_parameters + 0] - output_parameters_mean[0]); + // add squared deviation for vp + output_parameters_std[1] += (output_parameters[i * n_model_parameters + 1] - output_parameters_mean[1]) * (output_parameters[i * n_model_parameters + 1] - output_parameters_mean[1]); + } + } + // divide and take square root + output_parameters_std[0] = sqrt(output_parameters_std[0] / output_states_histogram[0]); + output_parameters_std[1] = sqrt(output_parameters_std[1] / output_states_histogram[0]); + + // print mean and std + std::cout << "offset true " << true_parameters[0] << " mean " << output_parameters_mean[0] << " std " << output_parameters_std[0] << "\n"; + std::cout << "slope true " << true_parameters[1] << " mean " << output_parameters_mean[1] << " std " << output_parameters_std[1] << "\n"; + + // compute mean chi-square for those converged + REAL output_chi_square_mean = 0; + for (size_t i = 0; i != n_fits; i++) + { + if (output_states[i] == FitState::CONVERGED) + { + output_chi_square_mean += output_chi_square[i]; + } + } + output_chi_square_mean /= static_cast(output_states_histogram[0]); + std::cout << "mean chi square " << output_chi_square_mean << "\n"; + + // compute mean number of iterations for those converged + REAL output_number_iterations_mean = 0; + for (size_t i = 0; i != n_fits; i++) + { + if (output_states[i] == FitState::CONVERGED) + { + output_number_iterations_mean += static_cast(output_number_iterations[i]); + } + } + + // normalize + output_number_iterations_mean /= static_cast(output_states_histogram[0]); + std::cout << "mean number of iterations " << output_number_iterations_mean << "\n"; } -int main(int argc, char *argv[]) +int main(int argc, char* argv[]) { - simple_example(); + std::cout << std::endl << "Beginning Patlak fit...v2" << std::endl; + patlak_two(); + + std::cout << std::endl << "Patlak fit completed!" << std::endl; + std::cout << "Press ENTER to exit" << std::endl; + std::getchar(); - std::cout << std::endl << "Example completed!" << std::endl; - std::cout << "Press ENTER to exit" << std::endl; - std::getchar(); - return 0; } diff --git a/Gpufit/models/models.cuh b/Gpufit/models/models.cuh index b25360d8..ef45e96b 100644 --- a/Gpufit/models/models.cuh +++ b/Gpufit/models/models.cuh @@ -9,6 +9,7 @@ #include "cauchy_2d_elliptic.cuh" #include "fletcher_powell_helix.cuh" #include "brown_dennis.cuh" +#include "patlak.cuh" __device__ void calculate_model( ModelID const model_id, @@ -49,7 +50,10 @@ __device__ void calculate_model( case BROWN_DENNIS: calculate_brown_dennis(parameters, n_fits, n_points, value, derivative, point_index, fit_index, chunk_index, user_info, user_info_size); break; - default: + case PATLAK: + calculate_patlak(parameters, n_fits, n_points, value, derivative, point_index, fit_index, chunk_index, user_info, user_info_size); + break; + default: break; } } @@ -66,6 +70,7 @@ void configure_model(ModelID const model_id, int & n_parameters, int & n_dimensi case LINEAR_1D: n_parameters = 2; n_dimensions = 1; break; case FLETCHER_POWELL_HELIX: n_parameters = 3; n_dimensions = 1; break; case BROWN_DENNIS: n_parameters = 4; n_dimensions = 1; break; + case PATLAK: n_parameters = 2; n_dimensions = 1; break; default: break; } } diff --git a/Gpufit/models/patlak.cuh b/Gpufit/models/patlak.cuh new file mode 100644 index 00000000..e44fff2e --- /dev/null +++ b/Gpufit/models/patlak.cuh @@ -0,0 +1,55 @@ +#ifndef GPUFIT_PATLAK_CUH_INCLUDED +#define GPUFIT_PATLAK_CUH_INCLUDED +#define REAL float + + +__device__ void calculate_patlak ( // function name + REAL const * parameters, + int const n_fits, + int const n_points, + REAL * value, + REAL * derivative, + int const point_index, // k + int const fit_index, + int const chunk_index, + char * user_info, // likely contains time and Cp values in 1 dimensional array + std::size_t const user_info_size) +{ + // indices + REAL* user_info_float = (REAL*)user_info; + //REAL x = 0; + //if (user_info_size / sizeof(REAL) == n_points) { // unnecessary since this case is always valid for this model? and setting independent variables is below. + // x = user_info_float[point_index]; + //} + + ///////////////////////////// value ////////////////////////////// + + // split user_info array into time and Cp + REAL *T = new REAL[n_points]; + for (int i = 0; i < n_points - 1; i++) + T[i] = user_info_float[i]; + + REAL *Cp = new REAL[n_points]; + for (int i = n_points - 1; i < 2 * n_points - 1; i++) + Cp[i] = user_info_float[i]; + + // integral (trapezoidal rule) + REAL area = 0; + for (int i = 1; i < point_index; i++) { // or is point_index i's limit? ya probably + // + REAL spacing = T[i] - T[i - 1]; + area += Cp[i] * spacing + 0.5 * (Cp[i] - Cp[i - 1]); + } + delete[] T; + delete[] Cp; + + value[point_index] = parameters[0] * area + parameters[1] * Cp[point_index]; // formula calculating fit model values + // C(t) = Ktrans * trapz(Cp(k)) + vp * Cp(k) + + /////////////////////////// derivative /////////////////////////// + REAL * current_derivative = derivative + point_index; + + current_derivative[0 * n_points] = Cp[point_index]; // formula calculating derivative values with respect to parameters[0] (Ktrans) + current_derivative[1 * n_points] = Cp[point_index]; // formula calculating derivative values with respect to parameters[1] (vp) +} +#endif \ No newline at end of file From 8418cdec82bd8479f1cae4f04c4cdaf875ca246b Mon Sep 17 00:00:00 2001 From: lsaca05 Date: Mon, 19 Aug 2019 22:38:39 -0700 Subject: [PATCH 3/9] +comment, -debugging cout line --- Gpufit/examples/Simple_Example.cpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/Gpufit/examples/Simple_Example.cpp b/Gpufit/examples/Simple_Example.cpp index ac9cf091..bef76ffe 100644 --- a/Gpufit/examples/Simple_Example.cpp +++ b/Gpufit/examples/Simple_Example.cpp @@ -33,11 +33,12 @@ void patlak_two() size_t const n_model_parameters = 2; // custom x positions for the data points of every fit, stored in user info + // time independent variable, given in minutes REAL time[] = { 0.25, 0.5, 0.75, 1, 1.25, 1.5, 1.75, 2, 2.25, 2.5, 2.75, 3, 3.25, 3.5, 3.75, 4, 4.25, 4.5, 4.75, 5, 5.25, 5.5, 5.75, 6, 6.25, 6.5, 6.75, 7, 7.25, 7.5, 7.75, 8, 8.25, 8.5, 8.75, 9, 9.25, 9.5, 9.75, 10, - 10.25, 10.5, 10.75, 11, 11.25, 11.5, 11.75, 12, 12.25, 12.5, 12.75, 13, 13.25, 13.5, 13.75, 14, 14.25, 14.5, 14.75, 15 }; // in min. up to 15 min + 10.25, 10.5, 10.75, 11, 11.25, 11.5, 11.75, 12, 12.25, 12.5, 12.75, 13, 13.25, 13.5, 13.75, 14, 14.25, 14.5, 14.75, 15 }; - REAL Cp[] = { 0.0f, 0.0f, 0.0f, 3.01846399851715f, 2.59801604007558f, 2.2361331285733f, 1.92465762011135f, 1.65656816551711f, 1.4258214335524f, + REAL Cp[] = { 0.0f, 0.0f, 0.0f, 3.01846399851715f, 2.59801604007558f, 2.2361331285733f, 1.92465762011135f, 1.65656816551711f, 1.4258214335524f, 1.22721588081636f, 1.05627449741415f, 0.909143885218726f, 0.782507393725825f, 0.6735103553914f, 0.579695735090254f, 0.498948743091769f, 0.429449163006342f, 0.369630320068624f, 0.318143764811612f, 0.273828876023252f, 0.235686697768721f, 0.20285742070682f, 0.174601000079374f, 0.150280473460109f, 0.12934760220805f, 0.111330512951924f, 0.0958230605172143f, @@ -69,7 +70,6 @@ void patlak_two() } std::cout << user_info.size() << std::endl; - std::cout << user_info.max_size() << std::endl; // size of user info in bytes size_t const user_info_size = 2 * n_points_per_fit * sizeof(REAL); From 01da2c6a1b2fae42955b32743c8bbdd0ee7c49ad Mon Sep 17 00:00:00 2001 From: S Barnes Date: Tue, 20 Aug 2019 16:51:31 -0700 Subject: [PATCH 4/9] Fixed bugs in cuda patlak, fitting working --- Gpufit/examples/Linear_Regression_Example.cpp | 4 +- Gpufit/examples/Simple_Example.cpp | 88 +++++++++++-------- Gpufit/models/patlak.cuh | 36 ++++---- 3 files changed, 71 insertions(+), 57 deletions(-) diff --git a/Gpufit/examples/Linear_Regression_Example.cpp b/Gpufit/examples/Linear_Regression_Example.cpp index cb4b6e1a..93af2830 100644 --- a/Gpufit/examples/Linear_Regression_Example.cpp +++ b/Gpufit/examples/Linear_Regression_Example.cpp @@ -27,7 +27,7 @@ void linear_regression_example() // number of fits, fit points and parameters size_t const n_fits = 10000; - size_t const n_points_per_fit = 20; + size_t const n_points_per_fit = 60; size_t const n_model_parameters = 2; // custom x positions for the data points of every fit, stored in user info @@ -38,7 +38,7 @@ void linear_regression_example() } // size of user info in bytes - size_t const user_info_size = n_points_per_fit * sizeof(REAL); + size_t const user_info_size = n_points_per_fit * sizeof(REAL); // initialize random number generator std::mt19937 rng; diff --git a/Gpufit/examples/Simple_Example.cpp b/Gpufit/examples/Simple_Example.cpp index bef76ffe..ac91b00e 100644 --- a/Gpufit/examples/Simple_Example.cpp +++ b/Gpufit/examples/Simple_Example.cpp @@ -1,5 +1,6 @@ #include "../gpufit.h" +#include #include #include #include @@ -27,14 +28,20 @@ void patlak_two() - and the mean number of iterations needed. */ + // start timer + clock_t time_start, time_end; + time_start = clock(); + + // number of fits, fit points and parameters - size_t const n_fits = 100; // 10000? 100? + size_t const n_fits = 10000; // 10000? 100? size_t const n_points_per_fit = 60; size_t const n_model_parameters = 2; + REAL snr = 0.8; // custom x positions for the data points of every fit, stored in user info // time independent variable, given in minutes - REAL time[] = { 0.25, 0.5, 0.75, 1, 1.25, 1.5, 1.75, 2, 2.25, 2.5, 2.75, 3, 3.25, 3.5, 3.75, 4, 4.25, 4.5, 4.75, 5, + REAL timeX[] = { 0.25, 0.5, 0.75, 1, 1.25, 1.5, 1.75, 2, 2.25, 2.5, 2.75, 3, 3.25, 3.5, 3.75, 4, 4.25, 4.5, 4.75, 5, 5.25, 5.5, 5.75, 6, 6.25, 6.5, 6.75, 7, 7.25, 7.5, 7.75, 8, 8.25, 8.5, 8.75, 9, 9.25, 9.5, 9.75, 10, 10.25, 10.5, 10.75, 11, 11.25, 11.5, 11.75, 12, 12.25, 12.5, 12.75, 13, 13.25, 13.5, 13.75, 14, 14.25, 14.5, 14.75, 15 }; @@ -53,67 +60,70 @@ void patlak_two() std::vector< REAL > user_info(2 * n_points_per_fit); for (size_t i = 0; i < n_points_per_fit; i++) { - user_info[i] = static_cast(time[i]); - std::cout << user_info[i] << std::endl; + user_info[i] = static_cast(timeX[i]); } - std::cout << "Time filled..." << std::endl; for (size_t i = n_points_per_fit; i < 2 * n_points_per_fit; i++) - { - user_info[i] = static_cast(Cp[i - n_points_per_fit]); - std::cout << user_info[i] << std::endl; - } - std::cout << "Checking user_info vector..." << std::endl; - for (size_t i = 0; i < 2 * n_points_per_fit; i++) { - std::cout << user_info[i] << std::endl; + user_info[i] = static_cast(Cp[i - n_points_per_fit]); } - std::cout << user_info.size() << std::endl; + //std::cout << user_info.size() << std::endl; // size of user info in bytes size_t const user_info_size = 2 * n_points_per_fit * sizeof(REAL); // initialize random number generator std::mt19937 rng; - rng.seed(0); + rng.seed(time(NULL)); + //rng.seed(0); std::uniform_real_distribution< REAL > uniform_dist(0, 1); std::normal_distribution< REAL > normal_dist(0, 1); // true parameters - std::vector< REAL > true_parameters{ 0.005, 0.02 }; // Ktrans, vp -- ? + std::vector< REAL > true_parameters{ 0.05, 0.03 }; // Ktrans, vp -- ? // initial parameters (randomized) std::vector< REAL > initial_parameters(n_fits * n_model_parameters); for (size_t i = 0; i != n_fits; i++) { // random offset - initial_parameters[i * n_model_parameters + 0] = true_parameters[0] * (0.8f + 0.4f * uniform_dist(rng)); + initial_parameters[i * n_model_parameters + 0] = true_parameters[0] * (0.1f + 1.8f * uniform_dist(rng)); // random slope - initial_parameters[i * n_model_parameters + 1] = true_parameters[0] * (0.8f + 0.4f * uniform_dist(rng)); + initial_parameters[i * n_model_parameters + 1] = true_parameters[0] * (0.1f + 1.8f * uniform_dist(rng)); } - std::cout << std::endl << "parameters randomized..." << std::endl; - - // generate data -- NEEDS TO BE UPDATED WITH PATLAK MODEL + // generate data -- NEEDS TO BE CHECKED WITH PATLAK MODEL std::vector< REAL > data(n_points_per_fit * n_fits); + REAL mean_y = 0; for (size_t i = 0; i != data.size(); i++) { size_t j = i / n_points_per_fit; // the fit size_t k = i % n_points_per_fit; // the position within a fit REAL x = 0; - for (int n = 1; n < k; n++) { // or is point_index i's limit? ya probably + for (int n = 1; n < k; n++) { // possible point of error // - REAL spacing = time[n] - time[n - 1]; - x += Cp[n] * spacing + 0.5 * (Cp[n] - Cp[n - 1]); + REAL spacing = timeX[n] - timeX[n - 1]; + x += (Cp[n - 1] + Cp[n]) / 2 * spacing; } - //REAL x = user_info[k]; // how to define w/ 2 x's? integral of Cp(k) to point_index? - REAL y = true_parameters[0] * x + true_parameters[1] * user_info[k + n_points_per_fit]; - data[i] = y + normal_dist(rng); + REAL y = true_parameters[0] * x + true_parameters[1] * Cp[k]; + //data[i] = y + normal_dist(rng); + //data[i] = y * (0.2f + 1.6f * uniform_dist(rng)); + data[i] = y; + mean_y += y; + //std::cout << data[i] << std::endl; + } + mean_y = mean_y / data.size(); + std::normal_distribution norm_snr(0,mean_y/snr); + for (size_t i = 0; i != data.size(); i++) + { + data[i] = data[i] + norm_snr(rng); } + + // tolerance - REAL const tolerance = 0.001f; // calculated by Christian to be 10e-15f + REAL const tolerance = 10e-8f; // calculated by Christian to be 10e-15f // maximum number of iterations int const max_number_iterations = 200; @@ -154,8 +164,6 @@ void patlak_two() output_number_iterations.data() ); - std::cout << std::endl << "gpufit complete." << std::endl; - // check status if (status != ReturnState::OK) @@ -164,9 +172,6 @@ void patlak_two() } - std::cout << std::endl << "status is ok." << std::endl; - - // get fit states std::vector< int > output_states_histogram(5, 0); for (std::vector< int >::iterator it = output_states.begin(); it != output_states.end(); ++it) @@ -182,6 +187,7 @@ void patlak_two() // compute mean fitted parameters for converged fits std::vector< REAL > output_parameters_mean(n_model_parameters, 0); + std::vector< REAL > output_parameters_mean_error(n_model_parameters, 0); for (size_t i = 0; i != n_fits; i++) { if (output_states[i] == FitState::CONVERGED) @@ -190,6 +196,10 @@ void patlak_two() output_parameters_mean[0] += output_parameters[i * n_model_parameters + 0]; // add vp output_parameters_mean[1] += output_parameters[i * n_model_parameters + 1]; + // add Ktrans + output_parameters_mean_error[0] += abs(output_parameters[i * n_model_parameters + 0]-true_parameters[0]); + // add vp + output_parameters_mean_error[1] += abs(output_parameters[i * n_model_parameters + 1]-true_parameters[1]); } } output_parameters_mean[0] /= output_states_histogram[0]; @@ -212,8 +222,8 @@ void patlak_two() output_parameters_std[1] = sqrt(output_parameters_std[1] / output_states_histogram[0]); // print mean and std - std::cout << "offset true " << true_parameters[0] << " mean " << output_parameters_mean[0] << " std " << output_parameters_std[0] << "\n"; - std::cout << "slope true " << true_parameters[1] << " mean " << output_parameters_mean[1] << " std " << output_parameters_std[1] << "\n"; + std::cout << "Ktrans true " << true_parameters[0] << " mean " << output_parameters_mean[0] << " std " << output_parameters_std[0] << "\n"; + std::cout << "vp true " << true_parameters[1] << " mean " << output_parameters_mean[1] << " std " << output_parameters_std[1] << "\n"; // compute mean chi-square for those converged REAL output_chi_square_mean = 0; @@ -240,17 +250,21 @@ void patlak_two() // normalize output_number_iterations_mean /= static_cast(output_states_histogram[0]); std::cout << "mean number of iterations " << output_number_iterations_mean << "\n"; + + // time + //time(&time_end); + time_end = clock(); + double time_taken_sec = double(time_end-time_start)/double(CLOCKS_PER_SEC); + std::cout << "execution time for " << n_fits << " fits was " << time_taken_sec << " seconds\n"; } int main(int argc, char* argv[]) { - std::cout << std::endl << "Beginning Patlak fit...v2" << std::endl; + std::cout << std::endl << "Beginning Patlak fit..." << std::endl; patlak_two(); std::cout << std::endl << "Patlak fit completed!" << std::endl; - std::cout << "Press ENTER to exit" << std::endl; - std::getchar(); return 0; } diff --git a/Gpufit/models/patlak.cuh b/Gpufit/models/patlak.cuh index e44fff2e..b777892f 100644 --- a/Gpufit/models/patlak.cuh +++ b/Gpufit/models/patlak.cuh @@ -1,7 +1,5 @@ #ifndef GPUFIT_PATLAK_CUH_INCLUDED #define GPUFIT_PATLAK_CUH_INCLUDED -#define REAL float - __device__ void calculate_patlak ( // function name REAL const * parameters, @@ -12,7 +10,7 @@ __device__ void calculate_patlak ( // function name int const point_index, // k int const fit_index, int const chunk_index, - char * user_info, // likely contains time and Cp values in 1 dimensional array + char * user_info, // contains time and Cp values in 1 dimensional array std::size_t const user_info_size) { // indices @@ -25,31 +23,33 @@ __device__ void calculate_patlak ( // function name ///////////////////////////// value ////////////////////////////// // split user_info array into time and Cp - REAL *T = new REAL[n_points]; - for (int i = 0; i < n_points - 1; i++) - T[i] = user_info_float[i]; + REAL* T = user_info_float; +// for (int i = 0; i < n_points; i++) +// T[i] = user_info_float[i]; - REAL *Cp = new REAL[n_points]; - for (int i = n_points - 1; i < 2 * n_points - 1; i++) - Cp[i] = user_info_float[i]; + REAL* Cp = user_info_float + n_points; +// for (int i = n_points; i < 2 * n_points; i++) +// Cp[i-n_points] = user_info_float[i]; // integral (trapezoidal rule) - REAL area = 0; - for (int i = 1; i < point_index; i++) { // or is point_index i's limit? ya probably - // + REAL convCp = 0; + for (int i = 1; i < point_index; i++) { REAL spacing = T[i] - T[i - 1]; - area += Cp[i] * spacing + 0.5 * (Cp[i] - Cp[i - 1]); + convCp += (Cp[i - 1] + Cp[i]) / 2 * spacing; } - delete[] T; - delete[] Cp; - value[point_index] = parameters[0] * area + parameters[1] * Cp[point_index]; // formula calculating fit model values + + value[point_index] = parameters[0] * convCp + parameters[1] * Cp[point_index]; // formula calculating fit model values // C(t) = Ktrans * trapz(Cp(k)) + vp * Cp(k) /////////////////////////// derivative /////////////////////////// REAL * current_derivative = derivative + point_index; - current_derivative[0 * n_points] = Cp[point_index]; // formula calculating derivative values with respect to parameters[0] (Ktrans) + current_derivative[0 * n_points] = convCp; // formula calculating derivative values with respect to parameters[0] (Ktrans) current_derivative[1 * n_points] = Cp[point_index]; // formula calculating derivative values with respect to parameters[1] (vp) + + // deallocate pointers + delete T; + delete Cp; } -#endif \ No newline at end of file +#endif From c143adb5667977a9ac0036487a6c5e6335ef9d37 Mon Sep 17 00:00:00 2001 From: lsaca05 Date: Wed, 21 Aug 2019 23:46:28 -0700 Subject: [PATCH 5/9] Cleaned up completed patlak code, prep'd files for other 4 models. --- Gpufit/constants.h | 6 +- Gpufit/examples/2CXM_fitting.cpp | 267 +++++++++++++++++++++ Gpufit/examples/Patlak_fitting.cpp | 267 +++++++++++++++++++++ Gpufit/examples/Simple_Example.cpp | 11 +- Gpufit/examples/Tissue_Uptake_fitting.cpp | 267 +++++++++++++++++++++ Gpufit/examples/Tofts_Extended_fitting.cpp | 267 +++++++++++++++++++++ Gpufit/examples/Tofts_fitting.cpp | 267 +++++++++++++++++++++ Gpufit/models/models.cuh | 36 ++- Gpufit/models/patlak.cuh | 16 +- Gpufit/models/tissue_uptake.cuh | 45 ++++ Gpufit/models/tofts.cuh | 45 ++++ Gpufit/models/tofts_extended.cuh | 45 ++++ Gpufit/models/two-compartment_exchange.cuh | 45 ++++ 13 files changed, 1553 insertions(+), 31 deletions(-) create mode 100644 Gpufit/examples/2CXM_fitting.cpp create mode 100644 Gpufit/examples/Patlak_fitting.cpp create mode 100644 Gpufit/examples/Tissue_Uptake_fitting.cpp create mode 100644 Gpufit/examples/Tofts_Extended_fitting.cpp create mode 100644 Gpufit/examples/Tofts_fitting.cpp create mode 100644 Gpufit/models/tissue_uptake.cuh create mode 100644 Gpufit/models/tofts.cuh create mode 100644 Gpufit/models/tofts_extended.cuh create mode 100644 Gpufit/models/two-compartment_exchange.cuh diff --git a/Gpufit/constants.h b/Gpufit/constants.h index 72fe2657..980c0308 100644 --- a/Gpufit/constants.h +++ b/Gpufit/constants.h @@ -12,7 +12,11 @@ enum ModelID { LINEAR_1D = 5, FLETCHER_POWELL_HELIX = 6, BROWN_DENNIS = 7, - PATLAK = 12 + PATLAK = 12, + TOFTS = 13, + TOFTS_EXTENDED = 14, + TISSUE_UPTAKE = 15, + TWO_COMPARTMENT_EXCHANGE = 16 }; // estimator ID diff --git a/Gpufit/examples/2CXM_fitting.cpp b/Gpufit/examples/2CXM_fitting.cpp new file mode 100644 index 00000000..7b35fae4 --- /dev/null +++ b/Gpufit/examples/2CXM_fitting.cpp @@ -0,0 +1,267 @@ +#include "../gpufit.h" + +#include +#include +#include +#include +#include + +void patlak_two() +{ + + /* + This example generates test data in form of 10000 one dimensional linear + curves with the size of 20 data points per curve. It is noised by normal + distributed noise. The initial guesses were randomized, within a specified + range of the true value. The LINEAR_1D model is fitted to the test data sets + using the LSE estimator. The optional parameter user_info is used to pass + custom x positions of the data sets. The same x position values are used for + every fit. + + The console output shows + - the ratio of converged fits including ratios of not converged fits for + different reasons, + - the values of the true parameters and the mean values of the fitted + parameters including their standard deviation, + - the mean chi square value + - and the mean number of iterations needed. + */ + + // start timer + clock_t time_start, time_end; + time_start = clock(); + + + // number of fits, fit points and parameters + size_t const n_fits = 10000; + size_t const n_points_per_fit = 60; + size_t const n_model_parameters = 2; + REAL snr = 0.8; + + // custom x positions for the data points of every fit, stored in user info + // time independent variable, given in minutes + REAL timeX[] = { 0.25, 0.5, 0.75, 1, 1.25, 1.5, 1.75, 2, 2.25, 2.5, 2.75, 3, 3.25, 3.5, 3.75, 4, 4.25, 4.5, 4.75, 5, + 5.25, 5.5, 5.75, 6, 6.25, 6.5, 6.75, 7, 7.25, 7.5, 7.75, 8, 8.25, 8.5, 8.75, 9, 9.25, 9.5, 9.75, 10, + 10.25, 10.5, 10.75, 11, 11.25, 11.5, 11.75, 12, 12.25, 12.5, 12.75, 13, 13.25, 13.5, 13.75, 14, 14.25, 14.5, 14.75, 15 }; + + // Concentration of plasma (independent), at 1 min based on equation: Cp(t) = 5.5e^(-.6t) + REAL Cp[] = { 0.0f, 0.0f, 0.0f, 3.01846399851715f, 2.59801604007558f, 2.2361331285733f, 1.92465762011135f, 1.65656816551711f, 1.4258214335524f, + 1.22721588081636f, 1.05627449741415f, 0.909143885218726f, 0.782507393725825f, 0.6735103553914f, 0.579695735090254f, + 0.498948743091769f, 0.429449163006342f, 0.369630320068624f, 0.318143764811612f, 0.273828876023252f, 0.235686697768721f, + 0.20285742070682f, 0.174601000079374f, 0.150280473460109f, 0.12934760220805f, 0.111330512951924f, 0.0958230605172143f, + 0.0824756725126274f, 0.0709874691926393f, 0.0610994809603327f, 0.0525888106179893f, 0.0452636087696102f, 0.0389587491097867f, + 0.033532106110336f, 0.0288613511954976f, 0.0248411951843697f, 0.0213810148391187f, 0.018402810016092f, 0.0158394453694853f, + 0.013633136971665f, 0.0117341497352074f, 0.010099676273659f, 0.00869287192804919f, 0.00748202420651342f, 0.00643983791435146f, + 0.00554281985976681f, 0.00477074926518851f, 0.00410622194607174f, 0.00353425798195557f, 0.00304196403581309f, 0.00261824270962248f, + 0.00225354238438883f, 0.00193964190545541f, 0.00166946525943377f, 0.00143692206515917f, 0.00123677028298367f, 0.00106449804756952f, + 0.000916221960431984f, 0.000788599549519612f, 0.000678753922476738f }; + + + std::vector< REAL > user_info(2 * n_points_per_fit); + for (size_t i = 0; i < n_points_per_fit; i++) + { + user_info[i] = static_cast(timeX[i]); + } + + for (size_t i = n_points_per_fit; i < 2 * n_points_per_fit; i++) + { + user_info[i] = static_cast(Cp[i - n_points_per_fit]); + } + + // size of user info in bytes + size_t const user_info_size = 2 * n_points_per_fit * sizeof(REAL); + + // initialize random number generator + std::mt19937 rng; + rng.seed(time(NULL)); + std::uniform_real_distribution< REAL > uniform_dist(0, 1); + std::normal_distribution< REAL > normal_dist(0, 1); + + // true parameters + std::vector< REAL > true_parameters{ 0.05, 0.03 }; // Ktrans, vp + + // initial parameters (randomized) + std::vector< REAL > initial_parameters(n_fits * n_model_parameters); + for (size_t i = 0; i != n_fits; i++) + { + // random offset + initial_parameters[i * n_model_parameters + 0] = true_parameters[0] * (0.1f + 1.8f * uniform_dist(rng)); + // random slope + initial_parameters[i * n_model_parameters + 1] = true_parameters[0] * (0.1f + 1.8f * uniform_dist(rng)); + } + + // generate data + std::vector< REAL > data(n_points_per_fit * n_fits); + REAL mean_y = 0; + for (size_t i = 0; i != data.size(); i++) + { + size_t j = i / n_points_per_fit; // the fit + size_t k = i % n_points_per_fit; // the position within a fit + REAL x = 0; + for (int n = 1; n < k; n++) { + + REAL spacing = timeX[n] - timeX[n - 1]; + x += (Cp[n - 1] + Cp[n]) / 2 * spacing; + } + REAL y = true_parameters[0] * x + true_parameters[1] * Cp[k]; + //data[i] = y + normal_dist(rng); + //data[i] = y * (0.2f + 1.6f * uniform_dist(rng)); + data[i] = y; + mean_y += y; + //std::cout << data[i] << std::endl; + } + mean_y = mean_y / data.size(); + std::normal_distribution norm_snr(0,mean_y/snr); + for (size_t i = 0; i != data.size(); i++) + { + data[i] = data[i] + norm_snr(rng); + } + + + + // tolerance + REAL const tolerance = 10e-8f; + + // maximum number of iterations + int const max_number_iterations = 200; + + // estimator ID + int const estimator_id = LSE; + + // model ID + int const model_id = PATLAK; + + // parameters to fit (all of them) + std::vector< int > parameters_to_fit(n_model_parameters, 1); + + // output parameters + std::vector< REAL > output_parameters(n_fits * n_model_parameters); + std::vector< int > output_states(n_fits); + std::vector< REAL > output_chi_square(n_fits); + std::vector< int > output_number_iterations(n_fits); + + // call to gpufit (C interface) + int const status = gpufit + ( + n_fits, + n_points_per_fit, + data.data(), + 0, + model_id, + initial_parameters.data(), + tolerance, + max_number_iterations, + parameters_to_fit.data(), + estimator_id, + user_info_size, + reinterpret_cast< char* >( user_info.data() ), + output_parameters.data(), + output_states.data(), + output_chi_square.data(), + output_number_iterations.data() + ); + + + // check status + if (status != ReturnState::OK) + { + throw std::runtime_error(gpufit_get_last_error()); + } + + + // get fit states + std::vector< int > output_states_histogram(5, 0); + for (std::vector< int >::iterator it = output_states.begin(); it != output_states.end(); ++it) + { + output_states_histogram[*it]++; + } + + std::cout << "ratio converged " << (REAL)output_states_histogram[0] / n_fits << "\n"; + std::cout << "ratio max iteration exceeded " << (REAL)output_states_histogram[1] / n_fits << "\n"; + std::cout << "ratio singular hessian " << (REAL)output_states_histogram[2] / n_fits << "\n"; + std::cout << "ratio neg curvature MLE " << (REAL)output_states_histogram[3] / n_fits << "\n"; + std::cout << "ratio gpu not read " << (REAL)output_states_histogram[4] / n_fits << "\n"; + + // compute mean fitted parameters for converged fits + std::vector< REAL > output_parameters_mean(n_model_parameters, 0); + std::vector< REAL > output_parameters_mean_error(n_model_parameters, 0); + for (size_t i = 0; i != n_fits; i++) + { + if (output_states[i] == FitState::CONVERGED) + { + // add Ktrans + output_parameters_mean[0] += output_parameters[i * n_model_parameters + 0]; + // add vp + output_parameters_mean[1] += output_parameters[i * n_model_parameters + 1]; + // add Ktrans + output_parameters_mean_error[0] += abs(output_parameters[i * n_model_parameters + 0]-true_parameters[0]); + // add vp + output_parameters_mean_error[1] += abs(output_parameters[i * n_model_parameters + 1]-true_parameters[1]); + } + } + output_parameters_mean[0] /= output_states_histogram[0]; + output_parameters_mean[1] /= output_states_histogram[0]; + + // compute std of fitted parameters for converged fits + std::vector< REAL > output_parameters_std(n_model_parameters, 0); + for (size_t i = 0; i != n_fits; i++) + { + if (output_states[i] == FitState::CONVERGED) + { + // add squared deviation for Ktrans + output_parameters_std[0] += (output_parameters[i * n_model_parameters + 0] - output_parameters_mean[0]) * (output_parameters[i * n_model_parameters + 0] - output_parameters_mean[0]); + // add squared deviation for vp + output_parameters_std[1] += (output_parameters[i * n_model_parameters + 1] - output_parameters_mean[1]) * (output_parameters[i * n_model_parameters + 1] - output_parameters_mean[1]); + } + } + // divide and take square root + output_parameters_std[0] = sqrt(output_parameters_std[0] / output_states_histogram[0]); + output_parameters_std[1] = sqrt(output_parameters_std[1] / output_states_histogram[0]); + + // print mean and std + std::cout << "Ktrans true " << true_parameters[0] << " mean " << output_parameters_mean[0] << " std " << output_parameters_std[0] << "\n"; + std::cout << "vp true " << true_parameters[1] << " mean " << output_parameters_mean[1] << " std " << output_parameters_std[1] << "\n"; + + // compute mean chi-square for those converged + REAL output_chi_square_mean = 0; + for (size_t i = 0; i != n_fits; i++) + { + if (output_states[i] == FitState::CONVERGED) + { + output_chi_square_mean += output_chi_square[i]; + } + } + output_chi_square_mean /= static_cast(output_states_histogram[0]); + std::cout << "mean chi square " << output_chi_square_mean << "\n"; + + // compute mean number of iterations for those converged + REAL output_number_iterations_mean = 0; + for (size_t i = 0; i != n_fits; i++) + { + if (output_states[i] == FitState::CONVERGED) + { + output_number_iterations_mean += static_cast(output_number_iterations[i]); + } + } + + // normalize + output_number_iterations_mean /= static_cast(output_states_histogram[0]); + std::cout << "mean number of iterations " << output_number_iterations_mean << "\n"; + + // time + //time(&time_end); + time_end = clock(); + double time_taken_sec = double(time_end-time_start)/double(CLOCKS_PER_SEC); + std::cout << "execution time for " << n_fits << " fits was " << time_taken_sec << " seconds\n"; +} + + +int main(int argc, char* argv[]) +{ + std::cout << std::endl << "Beginning Patlak fit..." << std::endl; + patlak_two(); + + std::cout << std::endl << "Patlak fit completed!" << std::endl; + + return 0; +} diff --git a/Gpufit/examples/Patlak_fitting.cpp b/Gpufit/examples/Patlak_fitting.cpp new file mode 100644 index 00000000..7b35fae4 --- /dev/null +++ b/Gpufit/examples/Patlak_fitting.cpp @@ -0,0 +1,267 @@ +#include "../gpufit.h" + +#include +#include +#include +#include +#include + +void patlak_two() +{ + + /* + This example generates test data in form of 10000 one dimensional linear + curves with the size of 20 data points per curve. It is noised by normal + distributed noise. The initial guesses were randomized, within a specified + range of the true value. The LINEAR_1D model is fitted to the test data sets + using the LSE estimator. The optional parameter user_info is used to pass + custom x positions of the data sets. The same x position values are used for + every fit. + + The console output shows + - the ratio of converged fits including ratios of not converged fits for + different reasons, + - the values of the true parameters and the mean values of the fitted + parameters including their standard deviation, + - the mean chi square value + - and the mean number of iterations needed. + */ + + // start timer + clock_t time_start, time_end; + time_start = clock(); + + + // number of fits, fit points and parameters + size_t const n_fits = 10000; + size_t const n_points_per_fit = 60; + size_t const n_model_parameters = 2; + REAL snr = 0.8; + + // custom x positions for the data points of every fit, stored in user info + // time independent variable, given in minutes + REAL timeX[] = { 0.25, 0.5, 0.75, 1, 1.25, 1.5, 1.75, 2, 2.25, 2.5, 2.75, 3, 3.25, 3.5, 3.75, 4, 4.25, 4.5, 4.75, 5, + 5.25, 5.5, 5.75, 6, 6.25, 6.5, 6.75, 7, 7.25, 7.5, 7.75, 8, 8.25, 8.5, 8.75, 9, 9.25, 9.5, 9.75, 10, + 10.25, 10.5, 10.75, 11, 11.25, 11.5, 11.75, 12, 12.25, 12.5, 12.75, 13, 13.25, 13.5, 13.75, 14, 14.25, 14.5, 14.75, 15 }; + + // Concentration of plasma (independent), at 1 min based on equation: Cp(t) = 5.5e^(-.6t) + REAL Cp[] = { 0.0f, 0.0f, 0.0f, 3.01846399851715f, 2.59801604007558f, 2.2361331285733f, 1.92465762011135f, 1.65656816551711f, 1.4258214335524f, + 1.22721588081636f, 1.05627449741415f, 0.909143885218726f, 0.782507393725825f, 0.6735103553914f, 0.579695735090254f, + 0.498948743091769f, 0.429449163006342f, 0.369630320068624f, 0.318143764811612f, 0.273828876023252f, 0.235686697768721f, + 0.20285742070682f, 0.174601000079374f, 0.150280473460109f, 0.12934760220805f, 0.111330512951924f, 0.0958230605172143f, + 0.0824756725126274f, 0.0709874691926393f, 0.0610994809603327f, 0.0525888106179893f, 0.0452636087696102f, 0.0389587491097867f, + 0.033532106110336f, 0.0288613511954976f, 0.0248411951843697f, 0.0213810148391187f, 0.018402810016092f, 0.0158394453694853f, + 0.013633136971665f, 0.0117341497352074f, 0.010099676273659f, 0.00869287192804919f, 0.00748202420651342f, 0.00643983791435146f, + 0.00554281985976681f, 0.00477074926518851f, 0.00410622194607174f, 0.00353425798195557f, 0.00304196403581309f, 0.00261824270962248f, + 0.00225354238438883f, 0.00193964190545541f, 0.00166946525943377f, 0.00143692206515917f, 0.00123677028298367f, 0.00106449804756952f, + 0.000916221960431984f, 0.000788599549519612f, 0.000678753922476738f }; + + + std::vector< REAL > user_info(2 * n_points_per_fit); + for (size_t i = 0; i < n_points_per_fit; i++) + { + user_info[i] = static_cast(timeX[i]); + } + + for (size_t i = n_points_per_fit; i < 2 * n_points_per_fit; i++) + { + user_info[i] = static_cast(Cp[i - n_points_per_fit]); + } + + // size of user info in bytes + size_t const user_info_size = 2 * n_points_per_fit * sizeof(REAL); + + // initialize random number generator + std::mt19937 rng; + rng.seed(time(NULL)); + std::uniform_real_distribution< REAL > uniform_dist(0, 1); + std::normal_distribution< REAL > normal_dist(0, 1); + + // true parameters + std::vector< REAL > true_parameters{ 0.05, 0.03 }; // Ktrans, vp + + // initial parameters (randomized) + std::vector< REAL > initial_parameters(n_fits * n_model_parameters); + for (size_t i = 0; i != n_fits; i++) + { + // random offset + initial_parameters[i * n_model_parameters + 0] = true_parameters[0] * (0.1f + 1.8f * uniform_dist(rng)); + // random slope + initial_parameters[i * n_model_parameters + 1] = true_parameters[0] * (0.1f + 1.8f * uniform_dist(rng)); + } + + // generate data + std::vector< REAL > data(n_points_per_fit * n_fits); + REAL mean_y = 0; + for (size_t i = 0; i != data.size(); i++) + { + size_t j = i / n_points_per_fit; // the fit + size_t k = i % n_points_per_fit; // the position within a fit + REAL x = 0; + for (int n = 1; n < k; n++) { + + REAL spacing = timeX[n] - timeX[n - 1]; + x += (Cp[n - 1] + Cp[n]) / 2 * spacing; + } + REAL y = true_parameters[0] * x + true_parameters[1] * Cp[k]; + //data[i] = y + normal_dist(rng); + //data[i] = y * (0.2f + 1.6f * uniform_dist(rng)); + data[i] = y; + mean_y += y; + //std::cout << data[i] << std::endl; + } + mean_y = mean_y / data.size(); + std::normal_distribution norm_snr(0,mean_y/snr); + for (size_t i = 0; i != data.size(); i++) + { + data[i] = data[i] + norm_snr(rng); + } + + + + // tolerance + REAL const tolerance = 10e-8f; + + // maximum number of iterations + int const max_number_iterations = 200; + + // estimator ID + int const estimator_id = LSE; + + // model ID + int const model_id = PATLAK; + + // parameters to fit (all of them) + std::vector< int > parameters_to_fit(n_model_parameters, 1); + + // output parameters + std::vector< REAL > output_parameters(n_fits * n_model_parameters); + std::vector< int > output_states(n_fits); + std::vector< REAL > output_chi_square(n_fits); + std::vector< int > output_number_iterations(n_fits); + + // call to gpufit (C interface) + int const status = gpufit + ( + n_fits, + n_points_per_fit, + data.data(), + 0, + model_id, + initial_parameters.data(), + tolerance, + max_number_iterations, + parameters_to_fit.data(), + estimator_id, + user_info_size, + reinterpret_cast< char* >( user_info.data() ), + output_parameters.data(), + output_states.data(), + output_chi_square.data(), + output_number_iterations.data() + ); + + + // check status + if (status != ReturnState::OK) + { + throw std::runtime_error(gpufit_get_last_error()); + } + + + // get fit states + std::vector< int > output_states_histogram(5, 0); + for (std::vector< int >::iterator it = output_states.begin(); it != output_states.end(); ++it) + { + output_states_histogram[*it]++; + } + + std::cout << "ratio converged " << (REAL)output_states_histogram[0] / n_fits << "\n"; + std::cout << "ratio max iteration exceeded " << (REAL)output_states_histogram[1] / n_fits << "\n"; + std::cout << "ratio singular hessian " << (REAL)output_states_histogram[2] / n_fits << "\n"; + std::cout << "ratio neg curvature MLE " << (REAL)output_states_histogram[3] / n_fits << "\n"; + std::cout << "ratio gpu not read " << (REAL)output_states_histogram[4] / n_fits << "\n"; + + // compute mean fitted parameters for converged fits + std::vector< REAL > output_parameters_mean(n_model_parameters, 0); + std::vector< REAL > output_parameters_mean_error(n_model_parameters, 0); + for (size_t i = 0; i != n_fits; i++) + { + if (output_states[i] == FitState::CONVERGED) + { + // add Ktrans + output_parameters_mean[0] += output_parameters[i * n_model_parameters + 0]; + // add vp + output_parameters_mean[1] += output_parameters[i * n_model_parameters + 1]; + // add Ktrans + output_parameters_mean_error[0] += abs(output_parameters[i * n_model_parameters + 0]-true_parameters[0]); + // add vp + output_parameters_mean_error[1] += abs(output_parameters[i * n_model_parameters + 1]-true_parameters[1]); + } + } + output_parameters_mean[0] /= output_states_histogram[0]; + output_parameters_mean[1] /= output_states_histogram[0]; + + // compute std of fitted parameters for converged fits + std::vector< REAL > output_parameters_std(n_model_parameters, 0); + for (size_t i = 0; i != n_fits; i++) + { + if (output_states[i] == FitState::CONVERGED) + { + // add squared deviation for Ktrans + output_parameters_std[0] += (output_parameters[i * n_model_parameters + 0] - output_parameters_mean[0]) * (output_parameters[i * n_model_parameters + 0] - output_parameters_mean[0]); + // add squared deviation for vp + output_parameters_std[1] += (output_parameters[i * n_model_parameters + 1] - output_parameters_mean[1]) * (output_parameters[i * n_model_parameters + 1] - output_parameters_mean[1]); + } + } + // divide and take square root + output_parameters_std[0] = sqrt(output_parameters_std[0] / output_states_histogram[0]); + output_parameters_std[1] = sqrt(output_parameters_std[1] / output_states_histogram[0]); + + // print mean and std + std::cout << "Ktrans true " << true_parameters[0] << " mean " << output_parameters_mean[0] << " std " << output_parameters_std[0] << "\n"; + std::cout << "vp true " << true_parameters[1] << " mean " << output_parameters_mean[1] << " std " << output_parameters_std[1] << "\n"; + + // compute mean chi-square for those converged + REAL output_chi_square_mean = 0; + for (size_t i = 0; i != n_fits; i++) + { + if (output_states[i] == FitState::CONVERGED) + { + output_chi_square_mean += output_chi_square[i]; + } + } + output_chi_square_mean /= static_cast(output_states_histogram[0]); + std::cout << "mean chi square " << output_chi_square_mean << "\n"; + + // compute mean number of iterations for those converged + REAL output_number_iterations_mean = 0; + for (size_t i = 0; i != n_fits; i++) + { + if (output_states[i] == FitState::CONVERGED) + { + output_number_iterations_mean += static_cast(output_number_iterations[i]); + } + } + + // normalize + output_number_iterations_mean /= static_cast(output_states_histogram[0]); + std::cout << "mean number of iterations " << output_number_iterations_mean << "\n"; + + // time + //time(&time_end); + time_end = clock(); + double time_taken_sec = double(time_end-time_start)/double(CLOCKS_PER_SEC); + std::cout << "execution time for " << n_fits << " fits was " << time_taken_sec << " seconds\n"; +} + + +int main(int argc, char* argv[]) +{ + std::cout << std::endl << "Beginning Patlak fit..." << std::endl; + patlak_two(); + + std::cout << std::endl << "Patlak fit completed!" << std::endl; + + return 0; +} diff --git a/Gpufit/examples/Simple_Example.cpp b/Gpufit/examples/Simple_Example.cpp index ac91b00e..f943b751 100644 --- a/Gpufit/examples/Simple_Example.cpp +++ b/Gpufit/examples/Simple_Example.cpp @@ -5,7 +5,6 @@ #include #include #include -//using namespace std; void patlak_two() { @@ -34,7 +33,7 @@ void patlak_two() // number of fits, fit points and parameters - size_t const n_fits = 10000; // 10000? 100? + size_t const n_fits = 10000; size_t const n_points_per_fit = 60; size_t const n_model_parameters = 2; REAL snr = 0.8; @@ -45,6 +44,7 @@ void patlak_two() 5.25, 5.5, 5.75, 6, 6.25, 6.5, 6.75, 7, 7.25, 7.5, 7.75, 8, 8.25, 8.5, 8.75, 9, 9.25, 9.5, 9.75, 10, 10.25, 10.5, 10.75, 11, 11.25, 11.5, 11.75, 12, 12.25, 12.5, 12.75, 13, 13.25, 13.5, 13.75, 14, 14.25, 14.5, 14.75, 15 }; + // Concentration of plasma (independent) REAL Cp[] = { 0.0f, 0.0f, 0.0f, 3.01846399851715f, 2.59801604007558f, 2.2361331285733f, 1.92465762011135f, 1.65656816551711f, 1.4258214335524f, 1.22721588081636f, 1.05627449741415f, 0.909143885218726f, 0.782507393725825f, 0.6735103553914f, 0.579695735090254f, 0.498948743091769f, 0.429449163006342f, 0.369630320068624f, 0.318143764811612f, 0.273828876023252f, 0.235686697768721f, @@ -68,15 +68,12 @@ void patlak_two() user_info[i] = static_cast(Cp[i - n_points_per_fit]); } - //std::cout << user_info.size() << std::endl; - // size of user info in bytes size_t const user_info_size = 2 * n_points_per_fit * sizeof(REAL); // initialize random number generator std::mt19937 rng; rng.seed(time(NULL)); - //rng.seed(0); std::uniform_real_distribution< REAL > uniform_dist(0, 1); std::normal_distribution< REAL > normal_dist(0, 1); @@ -101,8 +98,8 @@ void patlak_two() size_t j = i / n_points_per_fit; // the fit size_t k = i % n_points_per_fit; // the position within a fit REAL x = 0; - for (int n = 1; n < k; n++) { // possible point of error - // + for (int n = 1; n < k; n++) { + REAL spacing = timeX[n] - timeX[n - 1]; x += (Cp[n - 1] + Cp[n]) / 2 * spacing; } diff --git a/Gpufit/examples/Tissue_Uptake_fitting.cpp b/Gpufit/examples/Tissue_Uptake_fitting.cpp new file mode 100644 index 00000000..7b35fae4 --- /dev/null +++ b/Gpufit/examples/Tissue_Uptake_fitting.cpp @@ -0,0 +1,267 @@ +#include "../gpufit.h" + +#include +#include +#include +#include +#include + +void patlak_two() +{ + + /* + This example generates test data in form of 10000 one dimensional linear + curves with the size of 20 data points per curve. It is noised by normal + distributed noise. The initial guesses were randomized, within a specified + range of the true value. The LINEAR_1D model is fitted to the test data sets + using the LSE estimator. The optional parameter user_info is used to pass + custom x positions of the data sets. The same x position values are used for + every fit. + + The console output shows + - the ratio of converged fits including ratios of not converged fits for + different reasons, + - the values of the true parameters and the mean values of the fitted + parameters including their standard deviation, + - the mean chi square value + - and the mean number of iterations needed. + */ + + // start timer + clock_t time_start, time_end; + time_start = clock(); + + + // number of fits, fit points and parameters + size_t const n_fits = 10000; + size_t const n_points_per_fit = 60; + size_t const n_model_parameters = 2; + REAL snr = 0.8; + + // custom x positions for the data points of every fit, stored in user info + // time independent variable, given in minutes + REAL timeX[] = { 0.25, 0.5, 0.75, 1, 1.25, 1.5, 1.75, 2, 2.25, 2.5, 2.75, 3, 3.25, 3.5, 3.75, 4, 4.25, 4.5, 4.75, 5, + 5.25, 5.5, 5.75, 6, 6.25, 6.5, 6.75, 7, 7.25, 7.5, 7.75, 8, 8.25, 8.5, 8.75, 9, 9.25, 9.5, 9.75, 10, + 10.25, 10.5, 10.75, 11, 11.25, 11.5, 11.75, 12, 12.25, 12.5, 12.75, 13, 13.25, 13.5, 13.75, 14, 14.25, 14.5, 14.75, 15 }; + + // Concentration of plasma (independent), at 1 min based on equation: Cp(t) = 5.5e^(-.6t) + REAL Cp[] = { 0.0f, 0.0f, 0.0f, 3.01846399851715f, 2.59801604007558f, 2.2361331285733f, 1.92465762011135f, 1.65656816551711f, 1.4258214335524f, + 1.22721588081636f, 1.05627449741415f, 0.909143885218726f, 0.782507393725825f, 0.6735103553914f, 0.579695735090254f, + 0.498948743091769f, 0.429449163006342f, 0.369630320068624f, 0.318143764811612f, 0.273828876023252f, 0.235686697768721f, + 0.20285742070682f, 0.174601000079374f, 0.150280473460109f, 0.12934760220805f, 0.111330512951924f, 0.0958230605172143f, + 0.0824756725126274f, 0.0709874691926393f, 0.0610994809603327f, 0.0525888106179893f, 0.0452636087696102f, 0.0389587491097867f, + 0.033532106110336f, 0.0288613511954976f, 0.0248411951843697f, 0.0213810148391187f, 0.018402810016092f, 0.0158394453694853f, + 0.013633136971665f, 0.0117341497352074f, 0.010099676273659f, 0.00869287192804919f, 0.00748202420651342f, 0.00643983791435146f, + 0.00554281985976681f, 0.00477074926518851f, 0.00410622194607174f, 0.00353425798195557f, 0.00304196403581309f, 0.00261824270962248f, + 0.00225354238438883f, 0.00193964190545541f, 0.00166946525943377f, 0.00143692206515917f, 0.00123677028298367f, 0.00106449804756952f, + 0.000916221960431984f, 0.000788599549519612f, 0.000678753922476738f }; + + + std::vector< REAL > user_info(2 * n_points_per_fit); + for (size_t i = 0; i < n_points_per_fit; i++) + { + user_info[i] = static_cast(timeX[i]); + } + + for (size_t i = n_points_per_fit; i < 2 * n_points_per_fit; i++) + { + user_info[i] = static_cast(Cp[i - n_points_per_fit]); + } + + // size of user info in bytes + size_t const user_info_size = 2 * n_points_per_fit * sizeof(REAL); + + // initialize random number generator + std::mt19937 rng; + rng.seed(time(NULL)); + std::uniform_real_distribution< REAL > uniform_dist(0, 1); + std::normal_distribution< REAL > normal_dist(0, 1); + + // true parameters + std::vector< REAL > true_parameters{ 0.05, 0.03 }; // Ktrans, vp + + // initial parameters (randomized) + std::vector< REAL > initial_parameters(n_fits * n_model_parameters); + for (size_t i = 0; i != n_fits; i++) + { + // random offset + initial_parameters[i * n_model_parameters + 0] = true_parameters[0] * (0.1f + 1.8f * uniform_dist(rng)); + // random slope + initial_parameters[i * n_model_parameters + 1] = true_parameters[0] * (0.1f + 1.8f * uniform_dist(rng)); + } + + // generate data + std::vector< REAL > data(n_points_per_fit * n_fits); + REAL mean_y = 0; + for (size_t i = 0; i != data.size(); i++) + { + size_t j = i / n_points_per_fit; // the fit + size_t k = i % n_points_per_fit; // the position within a fit + REAL x = 0; + for (int n = 1; n < k; n++) { + + REAL spacing = timeX[n] - timeX[n - 1]; + x += (Cp[n - 1] + Cp[n]) / 2 * spacing; + } + REAL y = true_parameters[0] * x + true_parameters[1] * Cp[k]; + //data[i] = y + normal_dist(rng); + //data[i] = y * (0.2f + 1.6f * uniform_dist(rng)); + data[i] = y; + mean_y += y; + //std::cout << data[i] << std::endl; + } + mean_y = mean_y / data.size(); + std::normal_distribution norm_snr(0,mean_y/snr); + for (size_t i = 0; i != data.size(); i++) + { + data[i] = data[i] + norm_snr(rng); + } + + + + // tolerance + REAL const tolerance = 10e-8f; + + // maximum number of iterations + int const max_number_iterations = 200; + + // estimator ID + int const estimator_id = LSE; + + // model ID + int const model_id = PATLAK; + + // parameters to fit (all of them) + std::vector< int > parameters_to_fit(n_model_parameters, 1); + + // output parameters + std::vector< REAL > output_parameters(n_fits * n_model_parameters); + std::vector< int > output_states(n_fits); + std::vector< REAL > output_chi_square(n_fits); + std::vector< int > output_number_iterations(n_fits); + + // call to gpufit (C interface) + int const status = gpufit + ( + n_fits, + n_points_per_fit, + data.data(), + 0, + model_id, + initial_parameters.data(), + tolerance, + max_number_iterations, + parameters_to_fit.data(), + estimator_id, + user_info_size, + reinterpret_cast< char* >( user_info.data() ), + output_parameters.data(), + output_states.data(), + output_chi_square.data(), + output_number_iterations.data() + ); + + + // check status + if (status != ReturnState::OK) + { + throw std::runtime_error(gpufit_get_last_error()); + } + + + // get fit states + std::vector< int > output_states_histogram(5, 0); + for (std::vector< int >::iterator it = output_states.begin(); it != output_states.end(); ++it) + { + output_states_histogram[*it]++; + } + + std::cout << "ratio converged " << (REAL)output_states_histogram[0] / n_fits << "\n"; + std::cout << "ratio max iteration exceeded " << (REAL)output_states_histogram[1] / n_fits << "\n"; + std::cout << "ratio singular hessian " << (REAL)output_states_histogram[2] / n_fits << "\n"; + std::cout << "ratio neg curvature MLE " << (REAL)output_states_histogram[3] / n_fits << "\n"; + std::cout << "ratio gpu not read " << (REAL)output_states_histogram[4] / n_fits << "\n"; + + // compute mean fitted parameters for converged fits + std::vector< REAL > output_parameters_mean(n_model_parameters, 0); + std::vector< REAL > output_parameters_mean_error(n_model_parameters, 0); + for (size_t i = 0; i != n_fits; i++) + { + if (output_states[i] == FitState::CONVERGED) + { + // add Ktrans + output_parameters_mean[0] += output_parameters[i * n_model_parameters + 0]; + // add vp + output_parameters_mean[1] += output_parameters[i * n_model_parameters + 1]; + // add Ktrans + output_parameters_mean_error[0] += abs(output_parameters[i * n_model_parameters + 0]-true_parameters[0]); + // add vp + output_parameters_mean_error[1] += abs(output_parameters[i * n_model_parameters + 1]-true_parameters[1]); + } + } + output_parameters_mean[0] /= output_states_histogram[0]; + output_parameters_mean[1] /= output_states_histogram[0]; + + // compute std of fitted parameters for converged fits + std::vector< REAL > output_parameters_std(n_model_parameters, 0); + for (size_t i = 0; i != n_fits; i++) + { + if (output_states[i] == FitState::CONVERGED) + { + // add squared deviation for Ktrans + output_parameters_std[0] += (output_parameters[i * n_model_parameters + 0] - output_parameters_mean[0]) * (output_parameters[i * n_model_parameters + 0] - output_parameters_mean[0]); + // add squared deviation for vp + output_parameters_std[1] += (output_parameters[i * n_model_parameters + 1] - output_parameters_mean[1]) * (output_parameters[i * n_model_parameters + 1] - output_parameters_mean[1]); + } + } + // divide and take square root + output_parameters_std[0] = sqrt(output_parameters_std[0] / output_states_histogram[0]); + output_parameters_std[1] = sqrt(output_parameters_std[1] / output_states_histogram[0]); + + // print mean and std + std::cout << "Ktrans true " << true_parameters[0] << " mean " << output_parameters_mean[0] << " std " << output_parameters_std[0] << "\n"; + std::cout << "vp true " << true_parameters[1] << " mean " << output_parameters_mean[1] << " std " << output_parameters_std[1] << "\n"; + + // compute mean chi-square for those converged + REAL output_chi_square_mean = 0; + for (size_t i = 0; i != n_fits; i++) + { + if (output_states[i] == FitState::CONVERGED) + { + output_chi_square_mean += output_chi_square[i]; + } + } + output_chi_square_mean /= static_cast(output_states_histogram[0]); + std::cout << "mean chi square " << output_chi_square_mean << "\n"; + + // compute mean number of iterations for those converged + REAL output_number_iterations_mean = 0; + for (size_t i = 0; i != n_fits; i++) + { + if (output_states[i] == FitState::CONVERGED) + { + output_number_iterations_mean += static_cast(output_number_iterations[i]); + } + } + + // normalize + output_number_iterations_mean /= static_cast(output_states_histogram[0]); + std::cout << "mean number of iterations " << output_number_iterations_mean << "\n"; + + // time + //time(&time_end); + time_end = clock(); + double time_taken_sec = double(time_end-time_start)/double(CLOCKS_PER_SEC); + std::cout << "execution time for " << n_fits << " fits was " << time_taken_sec << " seconds\n"; +} + + +int main(int argc, char* argv[]) +{ + std::cout << std::endl << "Beginning Patlak fit..." << std::endl; + patlak_two(); + + std::cout << std::endl << "Patlak fit completed!" << std::endl; + + return 0; +} diff --git a/Gpufit/examples/Tofts_Extended_fitting.cpp b/Gpufit/examples/Tofts_Extended_fitting.cpp new file mode 100644 index 00000000..7b35fae4 --- /dev/null +++ b/Gpufit/examples/Tofts_Extended_fitting.cpp @@ -0,0 +1,267 @@ +#include "../gpufit.h" + +#include +#include +#include +#include +#include + +void patlak_two() +{ + + /* + This example generates test data in form of 10000 one dimensional linear + curves with the size of 20 data points per curve. It is noised by normal + distributed noise. The initial guesses were randomized, within a specified + range of the true value. The LINEAR_1D model is fitted to the test data sets + using the LSE estimator. The optional parameter user_info is used to pass + custom x positions of the data sets. The same x position values are used for + every fit. + + The console output shows + - the ratio of converged fits including ratios of not converged fits for + different reasons, + - the values of the true parameters and the mean values of the fitted + parameters including their standard deviation, + - the mean chi square value + - and the mean number of iterations needed. + */ + + // start timer + clock_t time_start, time_end; + time_start = clock(); + + + // number of fits, fit points and parameters + size_t const n_fits = 10000; + size_t const n_points_per_fit = 60; + size_t const n_model_parameters = 2; + REAL snr = 0.8; + + // custom x positions for the data points of every fit, stored in user info + // time independent variable, given in minutes + REAL timeX[] = { 0.25, 0.5, 0.75, 1, 1.25, 1.5, 1.75, 2, 2.25, 2.5, 2.75, 3, 3.25, 3.5, 3.75, 4, 4.25, 4.5, 4.75, 5, + 5.25, 5.5, 5.75, 6, 6.25, 6.5, 6.75, 7, 7.25, 7.5, 7.75, 8, 8.25, 8.5, 8.75, 9, 9.25, 9.5, 9.75, 10, + 10.25, 10.5, 10.75, 11, 11.25, 11.5, 11.75, 12, 12.25, 12.5, 12.75, 13, 13.25, 13.5, 13.75, 14, 14.25, 14.5, 14.75, 15 }; + + // Concentration of plasma (independent), at 1 min based on equation: Cp(t) = 5.5e^(-.6t) + REAL Cp[] = { 0.0f, 0.0f, 0.0f, 3.01846399851715f, 2.59801604007558f, 2.2361331285733f, 1.92465762011135f, 1.65656816551711f, 1.4258214335524f, + 1.22721588081636f, 1.05627449741415f, 0.909143885218726f, 0.782507393725825f, 0.6735103553914f, 0.579695735090254f, + 0.498948743091769f, 0.429449163006342f, 0.369630320068624f, 0.318143764811612f, 0.273828876023252f, 0.235686697768721f, + 0.20285742070682f, 0.174601000079374f, 0.150280473460109f, 0.12934760220805f, 0.111330512951924f, 0.0958230605172143f, + 0.0824756725126274f, 0.0709874691926393f, 0.0610994809603327f, 0.0525888106179893f, 0.0452636087696102f, 0.0389587491097867f, + 0.033532106110336f, 0.0288613511954976f, 0.0248411951843697f, 0.0213810148391187f, 0.018402810016092f, 0.0158394453694853f, + 0.013633136971665f, 0.0117341497352074f, 0.010099676273659f, 0.00869287192804919f, 0.00748202420651342f, 0.00643983791435146f, + 0.00554281985976681f, 0.00477074926518851f, 0.00410622194607174f, 0.00353425798195557f, 0.00304196403581309f, 0.00261824270962248f, + 0.00225354238438883f, 0.00193964190545541f, 0.00166946525943377f, 0.00143692206515917f, 0.00123677028298367f, 0.00106449804756952f, + 0.000916221960431984f, 0.000788599549519612f, 0.000678753922476738f }; + + + std::vector< REAL > user_info(2 * n_points_per_fit); + for (size_t i = 0; i < n_points_per_fit; i++) + { + user_info[i] = static_cast(timeX[i]); + } + + for (size_t i = n_points_per_fit; i < 2 * n_points_per_fit; i++) + { + user_info[i] = static_cast(Cp[i - n_points_per_fit]); + } + + // size of user info in bytes + size_t const user_info_size = 2 * n_points_per_fit * sizeof(REAL); + + // initialize random number generator + std::mt19937 rng; + rng.seed(time(NULL)); + std::uniform_real_distribution< REAL > uniform_dist(0, 1); + std::normal_distribution< REAL > normal_dist(0, 1); + + // true parameters + std::vector< REAL > true_parameters{ 0.05, 0.03 }; // Ktrans, vp + + // initial parameters (randomized) + std::vector< REAL > initial_parameters(n_fits * n_model_parameters); + for (size_t i = 0; i != n_fits; i++) + { + // random offset + initial_parameters[i * n_model_parameters + 0] = true_parameters[0] * (0.1f + 1.8f * uniform_dist(rng)); + // random slope + initial_parameters[i * n_model_parameters + 1] = true_parameters[0] * (0.1f + 1.8f * uniform_dist(rng)); + } + + // generate data + std::vector< REAL > data(n_points_per_fit * n_fits); + REAL mean_y = 0; + for (size_t i = 0; i != data.size(); i++) + { + size_t j = i / n_points_per_fit; // the fit + size_t k = i % n_points_per_fit; // the position within a fit + REAL x = 0; + for (int n = 1; n < k; n++) { + + REAL spacing = timeX[n] - timeX[n - 1]; + x += (Cp[n - 1] + Cp[n]) / 2 * spacing; + } + REAL y = true_parameters[0] * x + true_parameters[1] * Cp[k]; + //data[i] = y + normal_dist(rng); + //data[i] = y * (0.2f + 1.6f * uniform_dist(rng)); + data[i] = y; + mean_y += y; + //std::cout << data[i] << std::endl; + } + mean_y = mean_y / data.size(); + std::normal_distribution norm_snr(0,mean_y/snr); + for (size_t i = 0; i != data.size(); i++) + { + data[i] = data[i] + norm_snr(rng); + } + + + + // tolerance + REAL const tolerance = 10e-8f; + + // maximum number of iterations + int const max_number_iterations = 200; + + // estimator ID + int const estimator_id = LSE; + + // model ID + int const model_id = PATLAK; + + // parameters to fit (all of them) + std::vector< int > parameters_to_fit(n_model_parameters, 1); + + // output parameters + std::vector< REAL > output_parameters(n_fits * n_model_parameters); + std::vector< int > output_states(n_fits); + std::vector< REAL > output_chi_square(n_fits); + std::vector< int > output_number_iterations(n_fits); + + // call to gpufit (C interface) + int const status = gpufit + ( + n_fits, + n_points_per_fit, + data.data(), + 0, + model_id, + initial_parameters.data(), + tolerance, + max_number_iterations, + parameters_to_fit.data(), + estimator_id, + user_info_size, + reinterpret_cast< char* >( user_info.data() ), + output_parameters.data(), + output_states.data(), + output_chi_square.data(), + output_number_iterations.data() + ); + + + // check status + if (status != ReturnState::OK) + { + throw std::runtime_error(gpufit_get_last_error()); + } + + + // get fit states + std::vector< int > output_states_histogram(5, 0); + for (std::vector< int >::iterator it = output_states.begin(); it != output_states.end(); ++it) + { + output_states_histogram[*it]++; + } + + std::cout << "ratio converged " << (REAL)output_states_histogram[0] / n_fits << "\n"; + std::cout << "ratio max iteration exceeded " << (REAL)output_states_histogram[1] / n_fits << "\n"; + std::cout << "ratio singular hessian " << (REAL)output_states_histogram[2] / n_fits << "\n"; + std::cout << "ratio neg curvature MLE " << (REAL)output_states_histogram[3] / n_fits << "\n"; + std::cout << "ratio gpu not read " << (REAL)output_states_histogram[4] / n_fits << "\n"; + + // compute mean fitted parameters for converged fits + std::vector< REAL > output_parameters_mean(n_model_parameters, 0); + std::vector< REAL > output_parameters_mean_error(n_model_parameters, 0); + for (size_t i = 0; i != n_fits; i++) + { + if (output_states[i] == FitState::CONVERGED) + { + // add Ktrans + output_parameters_mean[0] += output_parameters[i * n_model_parameters + 0]; + // add vp + output_parameters_mean[1] += output_parameters[i * n_model_parameters + 1]; + // add Ktrans + output_parameters_mean_error[0] += abs(output_parameters[i * n_model_parameters + 0]-true_parameters[0]); + // add vp + output_parameters_mean_error[1] += abs(output_parameters[i * n_model_parameters + 1]-true_parameters[1]); + } + } + output_parameters_mean[0] /= output_states_histogram[0]; + output_parameters_mean[1] /= output_states_histogram[0]; + + // compute std of fitted parameters for converged fits + std::vector< REAL > output_parameters_std(n_model_parameters, 0); + for (size_t i = 0; i != n_fits; i++) + { + if (output_states[i] == FitState::CONVERGED) + { + // add squared deviation for Ktrans + output_parameters_std[0] += (output_parameters[i * n_model_parameters + 0] - output_parameters_mean[0]) * (output_parameters[i * n_model_parameters + 0] - output_parameters_mean[0]); + // add squared deviation for vp + output_parameters_std[1] += (output_parameters[i * n_model_parameters + 1] - output_parameters_mean[1]) * (output_parameters[i * n_model_parameters + 1] - output_parameters_mean[1]); + } + } + // divide and take square root + output_parameters_std[0] = sqrt(output_parameters_std[0] / output_states_histogram[0]); + output_parameters_std[1] = sqrt(output_parameters_std[1] / output_states_histogram[0]); + + // print mean and std + std::cout << "Ktrans true " << true_parameters[0] << " mean " << output_parameters_mean[0] << " std " << output_parameters_std[0] << "\n"; + std::cout << "vp true " << true_parameters[1] << " mean " << output_parameters_mean[1] << " std " << output_parameters_std[1] << "\n"; + + // compute mean chi-square for those converged + REAL output_chi_square_mean = 0; + for (size_t i = 0; i != n_fits; i++) + { + if (output_states[i] == FitState::CONVERGED) + { + output_chi_square_mean += output_chi_square[i]; + } + } + output_chi_square_mean /= static_cast(output_states_histogram[0]); + std::cout << "mean chi square " << output_chi_square_mean << "\n"; + + // compute mean number of iterations for those converged + REAL output_number_iterations_mean = 0; + for (size_t i = 0; i != n_fits; i++) + { + if (output_states[i] == FitState::CONVERGED) + { + output_number_iterations_mean += static_cast(output_number_iterations[i]); + } + } + + // normalize + output_number_iterations_mean /= static_cast(output_states_histogram[0]); + std::cout << "mean number of iterations " << output_number_iterations_mean << "\n"; + + // time + //time(&time_end); + time_end = clock(); + double time_taken_sec = double(time_end-time_start)/double(CLOCKS_PER_SEC); + std::cout << "execution time for " << n_fits << " fits was " << time_taken_sec << " seconds\n"; +} + + +int main(int argc, char* argv[]) +{ + std::cout << std::endl << "Beginning Patlak fit..." << std::endl; + patlak_two(); + + std::cout << std::endl << "Patlak fit completed!" << std::endl; + + return 0; +} diff --git a/Gpufit/examples/Tofts_fitting.cpp b/Gpufit/examples/Tofts_fitting.cpp new file mode 100644 index 00000000..7b35fae4 --- /dev/null +++ b/Gpufit/examples/Tofts_fitting.cpp @@ -0,0 +1,267 @@ +#include "../gpufit.h" + +#include +#include +#include +#include +#include + +void patlak_two() +{ + + /* + This example generates test data in form of 10000 one dimensional linear + curves with the size of 20 data points per curve. It is noised by normal + distributed noise. The initial guesses were randomized, within a specified + range of the true value. The LINEAR_1D model is fitted to the test data sets + using the LSE estimator. The optional parameter user_info is used to pass + custom x positions of the data sets. The same x position values are used for + every fit. + + The console output shows + - the ratio of converged fits including ratios of not converged fits for + different reasons, + - the values of the true parameters and the mean values of the fitted + parameters including their standard deviation, + - the mean chi square value + - and the mean number of iterations needed. + */ + + // start timer + clock_t time_start, time_end; + time_start = clock(); + + + // number of fits, fit points and parameters + size_t const n_fits = 10000; + size_t const n_points_per_fit = 60; + size_t const n_model_parameters = 2; + REAL snr = 0.8; + + // custom x positions for the data points of every fit, stored in user info + // time independent variable, given in minutes + REAL timeX[] = { 0.25, 0.5, 0.75, 1, 1.25, 1.5, 1.75, 2, 2.25, 2.5, 2.75, 3, 3.25, 3.5, 3.75, 4, 4.25, 4.5, 4.75, 5, + 5.25, 5.5, 5.75, 6, 6.25, 6.5, 6.75, 7, 7.25, 7.5, 7.75, 8, 8.25, 8.5, 8.75, 9, 9.25, 9.5, 9.75, 10, + 10.25, 10.5, 10.75, 11, 11.25, 11.5, 11.75, 12, 12.25, 12.5, 12.75, 13, 13.25, 13.5, 13.75, 14, 14.25, 14.5, 14.75, 15 }; + + // Concentration of plasma (independent), at 1 min based on equation: Cp(t) = 5.5e^(-.6t) + REAL Cp[] = { 0.0f, 0.0f, 0.0f, 3.01846399851715f, 2.59801604007558f, 2.2361331285733f, 1.92465762011135f, 1.65656816551711f, 1.4258214335524f, + 1.22721588081636f, 1.05627449741415f, 0.909143885218726f, 0.782507393725825f, 0.6735103553914f, 0.579695735090254f, + 0.498948743091769f, 0.429449163006342f, 0.369630320068624f, 0.318143764811612f, 0.273828876023252f, 0.235686697768721f, + 0.20285742070682f, 0.174601000079374f, 0.150280473460109f, 0.12934760220805f, 0.111330512951924f, 0.0958230605172143f, + 0.0824756725126274f, 0.0709874691926393f, 0.0610994809603327f, 0.0525888106179893f, 0.0452636087696102f, 0.0389587491097867f, + 0.033532106110336f, 0.0288613511954976f, 0.0248411951843697f, 0.0213810148391187f, 0.018402810016092f, 0.0158394453694853f, + 0.013633136971665f, 0.0117341497352074f, 0.010099676273659f, 0.00869287192804919f, 0.00748202420651342f, 0.00643983791435146f, + 0.00554281985976681f, 0.00477074926518851f, 0.00410622194607174f, 0.00353425798195557f, 0.00304196403581309f, 0.00261824270962248f, + 0.00225354238438883f, 0.00193964190545541f, 0.00166946525943377f, 0.00143692206515917f, 0.00123677028298367f, 0.00106449804756952f, + 0.000916221960431984f, 0.000788599549519612f, 0.000678753922476738f }; + + + std::vector< REAL > user_info(2 * n_points_per_fit); + for (size_t i = 0; i < n_points_per_fit; i++) + { + user_info[i] = static_cast(timeX[i]); + } + + for (size_t i = n_points_per_fit; i < 2 * n_points_per_fit; i++) + { + user_info[i] = static_cast(Cp[i - n_points_per_fit]); + } + + // size of user info in bytes + size_t const user_info_size = 2 * n_points_per_fit * sizeof(REAL); + + // initialize random number generator + std::mt19937 rng; + rng.seed(time(NULL)); + std::uniform_real_distribution< REAL > uniform_dist(0, 1); + std::normal_distribution< REAL > normal_dist(0, 1); + + // true parameters + std::vector< REAL > true_parameters{ 0.05, 0.03 }; // Ktrans, vp + + // initial parameters (randomized) + std::vector< REAL > initial_parameters(n_fits * n_model_parameters); + for (size_t i = 0; i != n_fits; i++) + { + // random offset + initial_parameters[i * n_model_parameters + 0] = true_parameters[0] * (0.1f + 1.8f * uniform_dist(rng)); + // random slope + initial_parameters[i * n_model_parameters + 1] = true_parameters[0] * (0.1f + 1.8f * uniform_dist(rng)); + } + + // generate data + std::vector< REAL > data(n_points_per_fit * n_fits); + REAL mean_y = 0; + for (size_t i = 0; i != data.size(); i++) + { + size_t j = i / n_points_per_fit; // the fit + size_t k = i % n_points_per_fit; // the position within a fit + REAL x = 0; + for (int n = 1; n < k; n++) { + + REAL spacing = timeX[n] - timeX[n - 1]; + x += (Cp[n - 1] + Cp[n]) / 2 * spacing; + } + REAL y = true_parameters[0] * x + true_parameters[1] * Cp[k]; + //data[i] = y + normal_dist(rng); + //data[i] = y * (0.2f + 1.6f * uniform_dist(rng)); + data[i] = y; + mean_y += y; + //std::cout << data[i] << std::endl; + } + mean_y = mean_y / data.size(); + std::normal_distribution norm_snr(0,mean_y/snr); + for (size_t i = 0; i != data.size(); i++) + { + data[i] = data[i] + norm_snr(rng); + } + + + + // tolerance + REAL const tolerance = 10e-8f; + + // maximum number of iterations + int const max_number_iterations = 200; + + // estimator ID + int const estimator_id = LSE; + + // model ID + int const model_id = PATLAK; + + // parameters to fit (all of them) + std::vector< int > parameters_to_fit(n_model_parameters, 1); + + // output parameters + std::vector< REAL > output_parameters(n_fits * n_model_parameters); + std::vector< int > output_states(n_fits); + std::vector< REAL > output_chi_square(n_fits); + std::vector< int > output_number_iterations(n_fits); + + // call to gpufit (C interface) + int const status = gpufit + ( + n_fits, + n_points_per_fit, + data.data(), + 0, + model_id, + initial_parameters.data(), + tolerance, + max_number_iterations, + parameters_to_fit.data(), + estimator_id, + user_info_size, + reinterpret_cast< char* >( user_info.data() ), + output_parameters.data(), + output_states.data(), + output_chi_square.data(), + output_number_iterations.data() + ); + + + // check status + if (status != ReturnState::OK) + { + throw std::runtime_error(gpufit_get_last_error()); + } + + + // get fit states + std::vector< int > output_states_histogram(5, 0); + for (std::vector< int >::iterator it = output_states.begin(); it != output_states.end(); ++it) + { + output_states_histogram[*it]++; + } + + std::cout << "ratio converged " << (REAL)output_states_histogram[0] / n_fits << "\n"; + std::cout << "ratio max iteration exceeded " << (REAL)output_states_histogram[1] / n_fits << "\n"; + std::cout << "ratio singular hessian " << (REAL)output_states_histogram[2] / n_fits << "\n"; + std::cout << "ratio neg curvature MLE " << (REAL)output_states_histogram[3] / n_fits << "\n"; + std::cout << "ratio gpu not read " << (REAL)output_states_histogram[4] / n_fits << "\n"; + + // compute mean fitted parameters for converged fits + std::vector< REAL > output_parameters_mean(n_model_parameters, 0); + std::vector< REAL > output_parameters_mean_error(n_model_parameters, 0); + for (size_t i = 0; i != n_fits; i++) + { + if (output_states[i] == FitState::CONVERGED) + { + // add Ktrans + output_parameters_mean[0] += output_parameters[i * n_model_parameters + 0]; + // add vp + output_parameters_mean[1] += output_parameters[i * n_model_parameters + 1]; + // add Ktrans + output_parameters_mean_error[0] += abs(output_parameters[i * n_model_parameters + 0]-true_parameters[0]); + // add vp + output_parameters_mean_error[1] += abs(output_parameters[i * n_model_parameters + 1]-true_parameters[1]); + } + } + output_parameters_mean[0] /= output_states_histogram[0]; + output_parameters_mean[1] /= output_states_histogram[0]; + + // compute std of fitted parameters for converged fits + std::vector< REAL > output_parameters_std(n_model_parameters, 0); + for (size_t i = 0; i != n_fits; i++) + { + if (output_states[i] == FitState::CONVERGED) + { + // add squared deviation for Ktrans + output_parameters_std[0] += (output_parameters[i * n_model_parameters + 0] - output_parameters_mean[0]) * (output_parameters[i * n_model_parameters + 0] - output_parameters_mean[0]); + // add squared deviation for vp + output_parameters_std[1] += (output_parameters[i * n_model_parameters + 1] - output_parameters_mean[1]) * (output_parameters[i * n_model_parameters + 1] - output_parameters_mean[1]); + } + } + // divide and take square root + output_parameters_std[0] = sqrt(output_parameters_std[0] / output_states_histogram[0]); + output_parameters_std[1] = sqrt(output_parameters_std[1] / output_states_histogram[0]); + + // print mean and std + std::cout << "Ktrans true " << true_parameters[0] << " mean " << output_parameters_mean[0] << " std " << output_parameters_std[0] << "\n"; + std::cout << "vp true " << true_parameters[1] << " mean " << output_parameters_mean[1] << " std " << output_parameters_std[1] << "\n"; + + // compute mean chi-square for those converged + REAL output_chi_square_mean = 0; + for (size_t i = 0; i != n_fits; i++) + { + if (output_states[i] == FitState::CONVERGED) + { + output_chi_square_mean += output_chi_square[i]; + } + } + output_chi_square_mean /= static_cast(output_states_histogram[0]); + std::cout << "mean chi square " << output_chi_square_mean << "\n"; + + // compute mean number of iterations for those converged + REAL output_number_iterations_mean = 0; + for (size_t i = 0; i != n_fits; i++) + { + if (output_states[i] == FitState::CONVERGED) + { + output_number_iterations_mean += static_cast(output_number_iterations[i]); + } + } + + // normalize + output_number_iterations_mean /= static_cast(output_states_histogram[0]); + std::cout << "mean number of iterations " << output_number_iterations_mean << "\n"; + + // time + //time(&time_end); + time_end = clock(); + double time_taken_sec = double(time_end-time_start)/double(CLOCKS_PER_SEC); + std::cout << "execution time for " << n_fits << " fits was " << time_taken_sec << " seconds\n"; +} + + +int main(int argc, char* argv[]) +{ + std::cout << std::endl << "Beginning Patlak fit..." << std::endl; + patlak_two(); + + std::cout << std::endl << "Patlak fit completed!" << std::endl; + + return 0; +} diff --git a/Gpufit/models/models.cuh b/Gpufit/models/models.cuh index ef45e96b..e7ad1018 100644 --- a/Gpufit/models/models.cuh +++ b/Gpufit/models/models.cuh @@ -53,6 +53,18 @@ __device__ void calculate_model( case PATLAK: calculate_patlak(parameters, n_fits, n_points, value, derivative, point_index, fit_index, chunk_index, user_info, user_info_size); break; + case TOFTS: + calculate_tofts(parameters, n_fits, n_points, value, derivative, point_index, fit_index, chunk_index, user_info, user_info_size); + break; + case TOFTS_EXTENDED: + calculate_tofts_extended(parameters, n_fits, n_points, value, derivative, point_index, fit_index, chunk_index, user_info, user_info_size); + break; + case TISSUE_UPTAKE: + calculate_tissue_uptake(parameters, n_fits, n_points, value, derivative, point_index, fit_index, chunk_index, user_info, user_info_size); + break; + case TWO_COMPARTMENT_EXCHANGE: + calculate_two_compartment_exchange(parameters, n_fits, n_points, value, derivative, point_index, fit_index, chunk_index, user_info, user_info_size); + break; default: break; } @@ -62,16 +74,20 @@ void configure_model(ModelID const model_id, int & n_parameters, int & n_dimensi { switch (model_id) { - case GAUSS_1D: n_parameters = 4; n_dimensions = 1; break; - case GAUSS_2D: n_parameters = 5; n_dimensions = 2; break; - case GAUSS_2D_ELLIPTIC: n_parameters = 6; n_dimensions = 2; break; - case GAUSS_2D_ROTATED: n_parameters = 7; n_dimensions = 2; break; - case CAUCHY_2D_ELLIPTIC: n_parameters = 6; n_dimensions = 2; break; - case LINEAR_1D: n_parameters = 2; n_dimensions = 1; break; - case FLETCHER_POWELL_HELIX: n_parameters = 3; n_dimensions = 1; break; - case BROWN_DENNIS: n_parameters = 4; n_dimensions = 1; break; - case PATLAK: n_parameters = 2; n_dimensions = 1; break; - default: break; + case GAUSS_1D: n_parameters = 4; n_dimensions = 1; break; + case GAUSS_2D: n_parameters = 5; n_dimensions = 2; break; + case GAUSS_2D_ELLIPTIC: n_parameters = 6; n_dimensions = 2; break; + case GAUSS_2D_ROTATED: n_parameters = 7; n_dimensions = 2; break; + case CAUCHY_2D_ELLIPTIC: n_parameters = 6; n_dimensions = 2; break; + case LINEAR_1D: n_parameters = 2; n_dimensions = 1; break; + case FLETCHER_POWELL_HELIX: n_parameters = 3; n_dimensions = 1; break; + case BROWN_DENNIS: n_parameters = 4; n_dimensions = 1; break; + case PATLAK: n_parameters = 2; n_dimensions = 1; break; + case TOFTS: n_parameters = 2; n_dimensions = 1; break; + case TOFTS_EXTENDED: n_parameters = 3; n_dimensions = 1; break; + case TISSUE_UPTAKE: n_parameters = 3; n_dimensions = 1; break; + case TWO_COMPARTMENT_EXCHANGE: n_parameters = 4; n_dimensions = 1; break; + default: break; } } diff --git a/Gpufit/models/patlak.cuh b/Gpufit/models/patlak.cuh index b777892f..c6f6a53d 100644 --- a/Gpufit/models/patlak.cuh +++ b/Gpufit/models/patlak.cuh @@ -7,29 +7,20 @@ __device__ void calculate_patlak ( // function name int const n_points, REAL * value, REAL * derivative, - int const point_index, // k + int const point_index, int const fit_index, int const chunk_index, - char * user_info, // contains time and Cp values in 1 dimensional array + char * user_info, // contains time and Cp values in a 1 dimensional array std::size_t const user_info_size) { // indices REAL* user_info_float = (REAL*)user_info; - //REAL x = 0; - //if (user_info_size / sizeof(REAL) == n_points) { // unnecessary since this case is always valid for this model? and setting independent variables is below. - // x = user_info_float[point_index]; - //} ///////////////////////////// value ////////////////////////////// // split user_info array into time and Cp REAL* T = user_info_float; -// for (int i = 0; i < n_points; i++) -// T[i] = user_info_float[i]; - REAL* Cp = user_info_float + n_points; -// for (int i = n_points; i < 2 * n_points; i++) -// Cp[i-n_points] = user_info_float[i]; // integral (trapezoidal rule) REAL convCp = 0; @@ -38,9 +29,8 @@ __device__ void calculate_patlak ( // function name convCp += (Cp[i - 1] + Cp[i]) / 2 * spacing; } - value[point_index] = parameters[0] * convCp + parameters[1] * Cp[point_index]; // formula calculating fit model values - // C(t) = Ktrans * trapz(Cp(k)) + vp * Cp(k) + // C(t) = Ktrans * trapz(Cp(k)) + vp * Cp(k) /////////////////////////// derivative /////////////////////////// REAL * current_derivative = derivative + point_index; diff --git a/Gpufit/models/tissue_uptake.cuh b/Gpufit/models/tissue_uptake.cuh new file mode 100644 index 00000000..7092396f --- /dev/null +++ b/Gpufit/models/tissue_uptake.cuh @@ -0,0 +1,45 @@ +#ifndef GPUFIT_TISSUE_UPTAKE_CUH_INCLUDED +#define GPUFIT_TISSUE_UPTAKE_CUH_INCLUDED + +__device__ void calculate_patlak ( // function name + REAL const * parameters, + int const n_fits, + int const n_points, + REAL * value, + REAL * derivative, + int const point_index, + int const fit_index, + int const chunk_index, + char * user_info, // contains time and Cp values in 1 dimensional array + std::size_t const user_info_size) +{ + // indices + REAL* user_info_float = (REAL*)user_info; + + ///////////////////////////// value ////////////////////////////// + + // split user_info array into time and Cp + REAL* T = user_info_float; + REAL* Cp = user_info_float + n_points; + + // integral (trapezoidal rule) + REAL convCp = 0; + for (int i = 1; i < point_index; i++) { + REAL spacing = T[i] - T[i - 1]; + convCp += (Cp[i - 1] + Cp[i]) / 2 * spacing; + } + + value[point_index] = parameters[0] * convCp + parameters[1] * Cp[point_index]; // formula calculating fit model values + // C(t) = Ktrans * trapz(Cp(k)) + vp * Cp(k) + + /////////////////////////// derivative /////////////////////////// + REAL * current_derivative = derivative + point_index; + + current_derivative[0 * n_points] = convCp; // formula calculating derivative values with respect to parameters[0] (Ktrans) + current_derivative[1 * n_points] = Cp[point_index]; // formula calculating derivative values with respect to parameters[1] (vp) + + // deallocate pointers + delete T; + delete Cp; +} +#endif diff --git a/Gpufit/models/tofts.cuh b/Gpufit/models/tofts.cuh new file mode 100644 index 00000000..03890d17 --- /dev/null +++ b/Gpufit/models/tofts.cuh @@ -0,0 +1,45 @@ +#ifndef GPUFIT_TOFTS_CUH_INCLUDED +#define GPUFIT_TOFTS_CUH_INCLUDED + +__device__ void calculate_tofts ( // function name + REAL const * parameters, + int const n_fits, + int const n_points, + REAL * value, + REAL * derivative, + int const point_index, + int const fit_index, + int const chunk_index, + char * user_info, // contains time and Cp values in 1 dimensional array + std::size_t const user_info_size) +{ + // indices + REAL* user_info_float = (REAL*)user_info; + + ///////////////////////////// value ////////////////////////////// + + // split user_info array into time and Cp + REAL* T = user_info_float; + REAL* Cp = user_info_float + n_points; + + // integral (trapezoidal rule) + REAL convCp = 0; + for (int i = 1; i < point_index; i++) { + REAL spacing = T[i] - T[i - 1]; + convCp += (Cp[i - 1] + Cp[i]) / 2 * spacing; + } + + value[point_index] = parameters[0] * convCp + parameters[1] * Cp[point_index]; // formula calculating fit model values + // C(t) = Ktrans * trapz(Cp(k)) + vp * Cp(k) + + /////////////////////////// derivative /////////////////////////// + REAL * current_derivative = derivative + point_index; + + current_derivative[0 * n_points] = convCp; // formula calculating derivative values with respect to parameters[0] (Ktrans) + current_derivative[1 * n_points] = Cp[point_index]; // formula calculating derivative values with respect to parameters[1] (vp) + + // deallocate pointers + delete T; + delete Cp; +} +#endif diff --git a/Gpufit/models/tofts_extended.cuh b/Gpufit/models/tofts_extended.cuh new file mode 100644 index 00000000..70502841 --- /dev/null +++ b/Gpufit/models/tofts_extended.cuh @@ -0,0 +1,45 @@ +#ifndef GPUFIT_TOFTS_EXTENDED_CUH_INCLUDED +#define GPUFIT_TOFTS_EXTENDED_CUH_INCLUDED + +__device__ void calculate_patlak ( // function name + REAL const * parameters, + int const n_fits, + int const n_points, + REAL * value, + REAL * derivative, + int const point_index, + int const fit_index, + int const chunk_index, + char * user_info, // contains time and Cp values in 1 dimensional array + std::size_t const user_info_size) +{ + // indices + REAL* user_info_float = (REAL*)user_info; + + ///////////////////////////// value ////////////////////////////// + + // split user_info array into time and Cp + REAL* T = user_info_float; + REAL* Cp = user_info_float + n_points; + + // integral (trapezoidal rule) + REAL convCp = 0; + for (int i = 1; i < point_index; i++) { + REAL spacing = T[i] - T[i - 1]; + convCp += (Cp[i - 1] + Cp[i]) / 2 * spacing; + } + + value[point_index] = parameters[0] * convCp + parameters[1] * Cp[point_index]; // formula calculating fit model values + // C(t) = Ktrans * trapz(Cp(k)) + vp * Cp(k) + + /////////////////////////// derivative /////////////////////////// + REAL * current_derivative = derivative + point_index; + + current_derivative[0 * n_points] = convCp; // formula calculating derivative values with respect to parameters[0] (Ktrans) + current_derivative[1 * n_points] = Cp[point_index]; // formula calculating derivative values with respect to parameters[1] (vp) + + // deallocate pointers + delete T; + delete Cp; +} +#endif diff --git a/Gpufit/models/two-compartment_exchange.cuh b/Gpufit/models/two-compartment_exchange.cuh new file mode 100644 index 00000000..fc233471 --- /dev/null +++ b/Gpufit/models/two-compartment_exchange.cuh @@ -0,0 +1,45 @@ +#ifndef GPUFIT_TWO_COMPARTMENT_EXCHANGE_CUH_INCLUDED +#define GPUFIT_TWO_COMPARTMENT_EXCHANGE_CUH_INCLUDED + +__device__ void calculate_patlak ( // function name + REAL const * parameters, + int const n_fits, + int const n_points, + REAL * value, + REAL * derivative, + int const point_index, + int const fit_index, + int const chunk_index, + char * user_info, // contains time and Cp values in 1 dimensional array + std::size_t const user_info_size) +{ + // indices + REAL* user_info_float = (REAL*)user_info; + + ///////////////////////////// value ////////////////////////////// + + // split user_info array into time and Cp + REAL* T = user_info_float; + REAL* Cp = user_info_float + n_points; + + // integral (trapezoidal rule) + REAL convCp = 0; + for (int i = 1; i < point_index; i++) { + REAL spacing = T[i] - T[i - 1]; + convCp += (Cp[i - 1] + Cp[i]) / 2 * spacing; + } + + value[point_index] = parameters[0] * convCp + parameters[1] * Cp[point_index]; // formula calculating fit model values + // C(t) = Ktrans * trapz(Cp(k)) + vp * Cp(k) + + /////////////////////////// derivative /////////////////////////// + REAL * current_derivative = derivative + point_index; + + current_derivative[0 * n_points] = convCp; // formula calculating derivative values with respect to parameters[0] (Ktrans) + current_derivative[1 * n_points] = Cp[point_index]; // formula calculating derivative values with respect to parameters[1] (vp) + + // deallocate pointers + delete T; + delete Cp; +} +#endif From 221fb4c80d94d67329e96357f7171552fba930ea Mon Sep 17 00:00:00 2001 From: lsaca05 Date: Fri, 23 Aug 2019 00:56:49 -0700 Subject: [PATCH 6/9] Laid further groundwork for the additional models. Began implementation of Tofts model. --- Gpufit/examples/2CXM_fitting.cpp | 2 +- Gpufit/models/models.cuh | 4 ++++ Gpufit/models/tissue_uptake.cuh | 2 +- Gpufit/models/tofts.cuh | 13 ++++++------- Gpufit/models/tofts_extended.cuh | 2 +- Gpufit/models/two-compartment_exchange.cuh | 2 +- 6 files changed, 14 insertions(+), 11 deletions(-) diff --git a/Gpufit/examples/2CXM_fitting.cpp b/Gpufit/examples/2CXM_fitting.cpp index 7b35fae4..ccccdc82 100644 --- a/Gpufit/examples/2CXM_fitting.cpp +++ b/Gpufit/examples/2CXM_fitting.cpp @@ -6,7 +6,7 @@ #include #include -void patlak_two() +void two_compartment_exchange_complicated() { /* diff --git a/Gpufit/models/models.cuh b/Gpufit/models/models.cuh index e7ad1018..338b3ec5 100644 --- a/Gpufit/models/models.cuh +++ b/Gpufit/models/models.cuh @@ -10,6 +10,10 @@ #include "fletcher_powell_helix.cuh" #include "brown_dennis.cuh" #include "patlak.cuh" +#include "tofts.cuh" +#include "tofts_extended.cuh" +#include "tissue_uptake.cuh" +#include "two-compartment_exchange.cuh" __device__ void calculate_model( ModelID const model_id, diff --git a/Gpufit/models/tissue_uptake.cuh b/Gpufit/models/tissue_uptake.cuh index 7092396f..62ed62b7 100644 --- a/Gpufit/models/tissue_uptake.cuh +++ b/Gpufit/models/tissue_uptake.cuh @@ -1,7 +1,7 @@ #ifndef GPUFIT_TISSUE_UPTAKE_CUH_INCLUDED #define GPUFIT_TISSUE_UPTAKE_CUH_INCLUDED -__device__ void calculate_patlak ( // function name +__device__ void calculate_tissue_uptake ( // function name REAL const * parameters, int const n_fits, int const n_points, diff --git a/Gpufit/models/tofts.cuh b/Gpufit/models/tofts.cuh index 03890d17..9cbc31db 100644 --- a/Gpufit/models/tofts.cuh +++ b/Gpufit/models/tofts.cuh @@ -23,20 +23,19 @@ __device__ void calculate_tofts ( // function name REAL* Cp = user_info_float + n_points; // integral (trapezoidal rule) - REAL convCp = 0; + REAL convFunc = 0; for (int i = 1; i < point_index; i++) { - REAL spacing = T[i] - T[i - 1]; - convCp += (Cp[i - 1] + Cp[i]) / 2 * spacing; + } - value[point_index] = parameters[0] * convCp + parameters[1] * Cp[point_index]; // formula calculating fit model values - // C(t) = Ktrans * trapz(Cp(k)) + vp * Cp(k) + value[point_index] = parameters[0] * Cp[point_index] * convFunc; // formula calculating fit model values + // C(t) = Ktrans * Cp * trapz(e^(-Ktrans*t/ve) /////////////////////////// derivative /////////////////////////// REAL * current_derivative = derivative + point_index; - current_derivative[0 * n_points] = convCp; // formula calculating derivative values with respect to parameters[0] (Ktrans) - current_derivative[1 * n_points] = Cp[point_index]; // formula calculating derivative values with respect to parameters[1] (vp) + current_derivative[0 * n_points] = Cp[point_index]; // formula calculating derivative values with respect to parameters[0] (Ktrans) + current_derivative[1 * n_points] = Cp[point_index]; // formula calculating derivative values with respect to parameters[1] (ve) // deallocate pointers delete T; diff --git a/Gpufit/models/tofts_extended.cuh b/Gpufit/models/tofts_extended.cuh index 70502841..7beea4d3 100644 --- a/Gpufit/models/tofts_extended.cuh +++ b/Gpufit/models/tofts_extended.cuh @@ -1,7 +1,7 @@ #ifndef GPUFIT_TOFTS_EXTENDED_CUH_INCLUDED #define GPUFIT_TOFTS_EXTENDED_CUH_INCLUDED -__device__ void calculate_patlak ( // function name +__device__ void calculate_tofts_extended ( // function name REAL const * parameters, int const n_fits, int const n_points, diff --git a/Gpufit/models/two-compartment_exchange.cuh b/Gpufit/models/two-compartment_exchange.cuh index fc233471..3cbc9efd 100644 --- a/Gpufit/models/two-compartment_exchange.cuh +++ b/Gpufit/models/two-compartment_exchange.cuh @@ -1,7 +1,7 @@ #ifndef GPUFIT_TWO_COMPARTMENT_EXCHANGE_CUH_INCLUDED #define GPUFIT_TWO_COMPARTMENT_EXCHANGE_CUH_INCLUDED -__device__ void calculate_patlak ( // function name +__device__ void calculate_two_compartment_exchange ( // function name REAL const * parameters, int const n_fits, int const n_points, From 818106d96002c3c492a0651c995c0edc8289cbae Mon Sep 17 00:00:00 2001 From: S Barnes Date: Fri, 23 Aug 2019 12:00:33 -0700 Subject: [PATCH 7/9] Pushing my (Lucas's) code to other PC. I expect conflicts with GitHub. --- Gpufit/matlab/ModelID.m | 3 ++- Gpufit/models/patlak.cuh | 11 +---------- 2 files changed, 3 insertions(+), 11 deletions(-) diff --git a/Gpufit/matlab/ModelID.m b/Gpufit/matlab/ModelID.m index b7921fd7..861f4f92 100644 --- a/Gpufit/matlab/ModelID.m +++ b/Gpufit/matlab/ModelID.m @@ -8,5 +8,6 @@ LINEAR_1D = 5 FLETCHER_POWELL = 6 BROWN_DENNIS = 7 + PATLAK = 12 end -end \ No newline at end of file +end diff --git a/Gpufit/models/patlak.cuh b/Gpufit/models/patlak.cuh index b777892f..ba6f3acc 100644 --- a/Gpufit/models/patlak.cuh +++ b/Gpufit/models/patlak.cuh @@ -15,21 +15,12 @@ __device__ void calculate_patlak ( // function name { // indices REAL* user_info_float = (REAL*)user_info; - //REAL x = 0; - //if (user_info_size / sizeof(REAL) == n_points) { // unnecessary since this case is always valid for this model? and setting independent variables is below. - // x = user_info_float[point_index]; - //} ///////////////////////////// value ////////////////////////////// // split user_info array into time and Cp REAL* T = user_info_float; -// for (int i = 0; i < n_points; i++) -// T[i] = user_info_float[i]; - REAL* Cp = user_info_float + n_points; -// for (int i = n_points; i < 2 * n_points; i++) -// Cp[i-n_points] = user_info_float[i]; // integral (trapezoidal rule) REAL convCp = 0; @@ -40,7 +31,7 @@ __device__ void calculate_patlak ( // function name value[point_index] = parameters[0] * convCp + parameters[1] * Cp[point_index]; // formula calculating fit model values - // C(t) = Ktrans * trapz(Cp(k)) + vp * Cp(k) + // C(t) = Ktrans * trapz(Cp(k)) + vp * Cp(k) /////////////////////////// derivative /////////////////////////// REAL * current_derivative = derivative + point_index; From 35a7f9defa7fc73f5a10e9fc46013418c1496f89 Mon Sep 17 00:00:00 2001 From: Sam Barnes Date: Fri, 23 Aug 2019 12:52:57 -0700 Subject: [PATCH 8/9] Git scrub here. Penultimate merge for finished Patlak model. --- Gpufit/constants.h | 6 +- Gpufit/examples/2CXM_fitting.cpp | 267 ++++++++++++++++++++ Gpufit/examples/Patlak_fitting.cpp | 267 ++++++++++++++++++++ Gpufit/examples/Simple_Example.cpp | 281 ++++----------------- Gpufit/examples/Tissue_Uptake_fitting.cpp | 267 ++++++++++++++++++++ Gpufit/examples/Tofts_Extended_fitting.cpp | 267 ++++++++++++++++++++ Gpufit/examples/Tofts_fitting.cpp | 267 ++++++++++++++++++++ Gpufit/models/models.cuh | 40 ++- Gpufit/models/patlak.cuh | 9 +- Gpufit/models/tissue_uptake.cuh | 45 ++++ Gpufit/models/tofts.cuh | 44 ++++ Gpufit/models/tofts_extended.cuh | 45 ++++ Gpufit/models/two-compartment_exchange.cuh | 45 ++++ 13 files changed, 1609 insertions(+), 241 deletions(-) create mode 100644 Gpufit/examples/2CXM_fitting.cpp create mode 100644 Gpufit/examples/Patlak_fitting.cpp create mode 100644 Gpufit/examples/Tissue_Uptake_fitting.cpp create mode 100644 Gpufit/examples/Tofts_Extended_fitting.cpp create mode 100644 Gpufit/examples/Tofts_fitting.cpp create mode 100644 Gpufit/models/tissue_uptake.cuh create mode 100644 Gpufit/models/tofts.cuh create mode 100644 Gpufit/models/tofts_extended.cuh create mode 100644 Gpufit/models/two-compartment_exchange.cuh diff --git a/Gpufit/constants.h b/Gpufit/constants.h index 72fe2657..980c0308 100644 --- a/Gpufit/constants.h +++ b/Gpufit/constants.h @@ -12,7 +12,11 @@ enum ModelID { LINEAR_1D = 5, FLETCHER_POWELL_HELIX = 6, BROWN_DENNIS = 7, - PATLAK = 12 + PATLAK = 12, + TOFTS = 13, + TOFTS_EXTENDED = 14, + TISSUE_UPTAKE = 15, + TWO_COMPARTMENT_EXCHANGE = 16 }; // estimator ID diff --git a/Gpufit/examples/2CXM_fitting.cpp b/Gpufit/examples/2CXM_fitting.cpp new file mode 100644 index 00000000..ccccdc82 --- /dev/null +++ b/Gpufit/examples/2CXM_fitting.cpp @@ -0,0 +1,267 @@ +#include "../gpufit.h" + +#include +#include +#include +#include +#include + +void two_compartment_exchange_complicated() +{ + + /* + This example generates test data in form of 10000 one dimensional linear + curves with the size of 20 data points per curve. It is noised by normal + distributed noise. The initial guesses were randomized, within a specified + range of the true value. The LINEAR_1D model is fitted to the test data sets + using the LSE estimator. The optional parameter user_info is used to pass + custom x positions of the data sets. The same x position values are used for + every fit. + + The console output shows + - the ratio of converged fits including ratios of not converged fits for + different reasons, + - the values of the true parameters and the mean values of the fitted + parameters including their standard deviation, + - the mean chi square value + - and the mean number of iterations needed. + */ + + // start timer + clock_t time_start, time_end; + time_start = clock(); + + + // number of fits, fit points and parameters + size_t const n_fits = 10000; + size_t const n_points_per_fit = 60; + size_t const n_model_parameters = 2; + REAL snr = 0.8; + + // custom x positions for the data points of every fit, stored in user info + // time independent variable, given in minutes + REAL timeX[] = { 0.25, 0.5, 0.75, 1, 1.25, 1.5, 1.75, 2, 2.25, 2.5, 2.75, 3, 3.25, 3.5, 3.75, 4, 4.25, 4.5, 4.75, 5, + 5.25, 5.5, 5.75, 6, 6.25, 6.5, 6.75, 7, 7.25, 7.5, 7.75, 8, 8.25, 8.5, 8.75, 9, 9.25, 9.5, 9.75, 10, + 10.25, 10.5, 10.75, 11, 11.25, 11.5, 11.75, 12, 12.25, 12.5, 12.75, 13, 13.25, 13.5, 13.75, 14, 14.25, 14.5, 14.75, 15 }; + + // Concentration of plasma (independent), at 1 min based on equation: Cp(t) = 5.5e^(-.6t) + REAL Cp[] = { 0.0f, 0.0f, 0.0f, 3.01846399851715f, 2.59801604007558f, 2.2361331285733f, 1.92465762011135f, 1.65656816551711f, 1.4258214335524f, + 1.22721588081636f, 1.05627449741415f, 0.909143885218726f, 0.782507393725825f, 0.6735103553914f, 0.579695735090254f, + 0.498948743091769f, 0.429449163006342f, 0.369630320068624f, 0.318143764811612f, 0.273828876023252f, 0.235686697768721f, + 0.20285742070682f, 0.174601000079374f, 0.150280473460109f, 0.12934760220805f, 0.111330512951924f, 0.0958230605172143f, + 0.0824756725126274f, 0.0709874691926393f, 0.0610994809603327f, 0.0525888106179893f, 0.0452636087696102f, 0.0389587491097867f, + 0.033532106110336f, 0.0288613511954976f, 0.0248411951843697f, 0.0213810148391187f, 0.018402810016092f, 0.0158394453694853f, + 0.013633136971665f, 0.0117341497352074f, 0.010099676273659f, 0.00869287192804919f, 0.00748202420651342f, 0.00643983791435146f, + 0.00554281985976681f, 0.00477074926518851f, 0.00410622194607174f, 0.00353425798195557f, 0.00304196403581309f, 0.00261824270962248f, + 0.00225354238438883f, 0.00193964190545541f, 0.00166946525943377f, 0.00143692206515917f, 0.00123677028298367f, 0.00106449804756952f, + 0.000916221960431984f, 0.000788599549519612f, 0.000678753922476738f }; + + + std::vector< REAL > user_info(2 * n_points_per_fit); + for (size_t i = 0; i < n_points_per_fit; i++) + { + user_info[i] = static_cast(timeX[i]); + } + + for (size_t i = n_points_per_fit; i < 2 * n_points_per_fit; i++) + { + user_info[i] = static_cast(Cp[i - n_points_per_fit]); + } + + // size of user info in bytes + size_t const user_info_size = 2 * n_points_per_fit * sizeof(REAL); + + // initialize random number generator + std::mt19937 rng; + rng.seed(time(NULL)); + std::uniform_real_distribution< REAL > uniform_dist(0, 1); + std::normal_distribution< REAL > normal_dist(0, 1); + + // true parameters + std::vector< REAL > true_parameters{ 0.05, 0.03 }; // Ktrans, vp + + // initial parameters (randomized) + std::vector< REAL > initial_parameters(n_fits * n_model_parameters); + for (size_t i = 0; i != n_fits; i++) + { + // random offset + initial_parameters[i * n_model_parameters + 0] = true_parameters[0] * (0.1f + 1.8f * uniform_dist(rng)); + // random slope + initial_parameters[i * n_model_parameters + 1] = true_parameters[0] * (0.1f + 1.8f * uniform_dist(rng)); + } + + // generate data + std::vector< REAL > data(n_points_per_fit * n_fits); + REAL mean_y = 0; + for (size_t i = 0; i != data.size(); i++) + { + size_t j = i / n_points_per_fit; // the fit + size_t k = i % n_points_per_fit; // the position within a fit + REAL x = 0; + for (int n = 1; n < k; n++) { + + REAL spacing = timeX[n] - timeX[n - 1]; + x += (Cp[n - 1] + Cp[n]) / 2 * spacing; + } + REAL y = true_parameters[0] * x + true_parameters[1] * Cp[k]; + //data[i] = y + normal_dist(rng); + //data[i] = y * (0.2f + 1.6f * uniform_dist(rng)); + data[i] = y; + mean_y += y; + //std::cout << data[i] << std::endl; + } + mean_y = mean_y / data.size(); + std::normal_distribution norm_snr(0,mean_y/snr); + for (size_t i = 0; i != data.size(); i++) + { + data[i] = data[i] + norm_snr(rng); + } + + + + // tolerance + REAL const tolerance = 10e-8f; + + // maximum number of iterations + int const max_number_iterations = 200; + + // estimator ID + int const estimator_id = LSE; + + // model ID + int const model_id = PATLAK; + + // parameters to fit (all of them) + std::vector< int > parameters_to_fit(n_model_parameters, 1); + + // output parameters + std::vector< REAL > output_parameters(n_fits * n_model_parameters); + std::vector< int > output_states(n_fits); + std::vector< REAL > output_chi_square(n_fits); + std::vector< int > output_number_iterations(n_fits); + + // call to gpufit (C interface) + int const status = gpufit + ( + n_fits, + n_points_per_fit, + data.data(), + 0, + model_id, + initial_parameters.data(), + tolerance, + max_number_iterations, + parameters_to_fit.data(), + estimator_id, + user_info_size, + reinterpret_cast< char* >( user_info.data() ), + output_parameters.data(), + output_states.data(), + output_chi_square.data(), + output_number_iterations.data() + ); + + + // check status + if (status != ReturnState::OK) + { + throw std::runtime_error(gpufit_get_last_error()); + } + + + // get fit states + std::vector< int > output_states_histogram(5, 0); + for (std::vector< int >::iterator it = output_states.begin(); it != output_states.end(); ++it) + { + output_states_histogram[*it]++; + } + + std::cout << "ratio converged " << (REAL)output_states_histogram[0] / n_fits << "\n"; + std::cout << "ratio max iteration exceeded " << (REAL)output_states_histogram[1] / n_fits << "\n"; + std::cout << "ratio singular hessian " << (REAL)output_states_histogram[2] / n_fits << "\n"; + std::cout << "ratio neg curvature MLE " << (REAL)output_states_histogram[3] / n_fits << "\n"; + std::cout << "ratio gpu not read " << (REAL)output_states_histogram[4] / n_fits << "\n"; + + // compute mean fitted parameters for converged fits + std::vector< REAL > output_parameters_mean(n_model_parameters, 0); + std::vector< REAL > output_parameters_mean_error(n_model_parameters, 0); + for (size_t i = 0; i != n_fits; i++) + { + if (output_states[i] == FitState::CONVERGED) + { + // add Ktrans + output_parameters_mean[0] += output_parameters[i * n_model_parameters + 0]; + // add vp + output_parameters_mean[1] += output_parameters[i * n_model_parameters + 1]; + // add Ktrans + output_parameters_mean_error[0] += abs(output_parameters[i * n_model_parameters + 0]-true_parameters[0]); + // add vp + output_parameters_mean_error[1] += abs(output_parameters[i * n_model_parameters + 1]-true_parameters[1]); + } + } + output_parameters_mean[0] /= output_states_histogram[0]; + output_parameters_mean[1] /= output_states_histogram[0]; + + // compute std of fitted parameters for converged fits + std::vector< REAL > output_parameters_std(n_model_parameters, 0); + for (size_t i = 0; i != n_fits; i++) + { + if (output_states[i] == FitState::CONVERGED) + { + // add squared deviation for Ktrans + output_parameters_std[0] += (output_parameters[i * n_model_parameters + 0] - output_parameters_mean[0]) * (output_parameters[i * n_model_parameters + 0] - output_parameters_mean[0]); + // add squared deviation for vp + output_parameters_std[1] += (output_parameters[i * n_model_parameters + 1] - output_parameters_mean[1]) * (output_parameters[i * n_model_parameters + 1] - output_parameters_mean[1]); + } + } + // divide and take square root + output_parameters_std[0] = sqrt(output_parameters_std[0] / output_states_histogram[0]); + output_parameters_std[1] = sqrt(output_parameters_std[1] / output_states_histogram[0]); + + // print mean and std + std::cout << "Ktrans true " << true_parameters[0] << " mean " << output_parameters_mean[0] << " std " << output_parameters_std[0] << "\n"; + std::cout << "vp true " << true_parameters[1] << " mean " << output_parameters_mean[1] << " std " << output_parameters_std[1] << "\n"; + + // compute mean chi-square for those converged + REAL output_chi_square_mean = 0; + for (size_t i = 0; i != n_fits; i++) + { + if (output_states[i] == FitState::CONVERGED) + { + output_chi_square_mean += output_chi_square[i]; + } + } + output_chi_square_mean /= static_cast(output_states_histogram[0]); + std::cout << "mean chi square " << output_chi_square_mean << "\n"; + + // compute mean number of iterations for those converged + REAL output_number_iterations_mean = 0; + for (size_t i = 0; i != n_fits; i++) + { + if (output_states[i] == FitState::CONVERGED) + { + output_number_iterations_mean += static_cast(output_number_iterations[i]); + } + } + + // normalize + output_number_iterations_mean /= static_cast(output_states_histogram[0]); + std::cout << "mean number of iterations " << output_number_iterations_mean << "\n"; + + // time + //time(&time_end); + time_end = clock(); + double time_taken_sec = double(time_end-time_start)/double(CLOCKS_PER_SEC); + std::cout << "execution time for " << n_fits << " fits was " << time_taken_sec << " seconds\n"; +} + + +int main(int argc, char* argv[]) +{ + std::cout << std::endl << "Beginning Patlak fit..." << std::endl; + patlak_two(); + + std::cout << std::endl << "Patlak fit completed!" << std::endl; + + return 0; +} diff --git a/Gpufit/examples/Patlak_fitting.cpp b/Gpufit/examples/Patlak_fitting.cpp new file mode 100644 index 00000000..7b35fae4 --- /dev/null +++ b/Gpufit/examples/Patlak_fitting.cpp @@ -0,0 +1,267 @@ +#include "../gpufit.h" + +#include +#include +#include +#include +#include + +void patlak_two() +{ + + /* + This example generates test data in form of 10000 one dimensional linear + curves with the size of 20 data points per curve. It is noised by normal + distributed noise. The initial guesses were randomized, within a specified + range of the true value. The LINEAR_1D model is fitted to the test data sets + using the LSE estimator. The optional parameter user_info is used to pass + custom x positions of the data sets. The same x position values are used for + every fit. + + The console output shows + - the ratio of converged fits including ratios of not converged fits for + different reasons, + - the values of the true parameters and the mean values of the fitted + parameters including their standard deviation, + - the mean chi square value + - and the mean number of iterations needed. + */ + + // start timer + clock_t time_start, time_end; + time_start = clock(); + + + // number of fits, fit points and parameters + size_t const n_fits = 10000; + size_t const n_points_per_fit = 60; + size_t const n_model_parameters = 2; + REAL snr = 0.8; + + // custom x positions for the data points of every fit, stored in user info + // time independent variable, given in minutes + REAL timeX[] = { 0.25, 0.5, 0.75, 1, 1.25, 1.5, 1.75, 2, 2.25, 2.5, 2.75, 3, 3.25, 3.5, 3.75, 4, 4.25, 4.5, 4.75, 5, + 5.25, 5.5, 5.75, 6, 6.25, 6.5, 6.75, 7, 7.25, 7.5, 7.75, 8, 8.25, 8.5, 8.75, 9, 9.25, 9.5, 9.75, 10, + 10.25, 10.5, 10.75, 11, 11.25, 11.5, 11.75, 12, 12.25, 12.5, 12.75, 13, 13.25, 13.5, 13.75, 14, 14.25, 14.5, 14.75, 15 }; + + // Concentration of plasma (independent), at 1 min based on equation: Cp(t) = 5.5e^(-.6t) + REAL Cp[] = { 0.0f, 0.0f, 0.0f, 3.01846399851715f, 2.59801604007558f, 2.2361331285733f, 1.92465762011135f, 1.65656816551711f, 1.4258214335524f, + 1.22721588081636f, 1.05627449741415f, 0.909143885218726f, 0.782507393725825f, 0.6735103553914f, 0.579695735090254f, + 0.498948743091769f, 0.429449163006342f, 0.369630320068624f, 0.318143764811612f, 0.273828876023252f, 0.235686697768721f, + 0.20285742070682f, 0.174601000079374f, 0.150280473460109f, 0.12934760220805f, 0.111330512951924f, 0.0958230605172143f, + 0.0824756725126274f, 0.0709874691926393f, 0.0610994809603327f, 0.0525888106179893f, 0.0452636087696102f, 0.0389587491097867f, + 0.033532106110336f, 0.0288613511954976f, 0.0248411951843697f, 0.0213810148391187f, 0.018402810016092f, 0.0158394453694853f, + 0.013633136971665f, 0.0117341497352074f, 0.010099676273659f, 0.00869287192804919f, 0.00748202420651342f, 0.00643983791435146f, + 0.00554281985976681f, 0.00477074926518851f, 0.00410622194607174f, 0.00353425798195557f, 0.00304196403581309f, 0.00261824270962248f, + 0.00225354238438883f, 0.00193964190545541f, 0.00166946525943377f, 0.00143692206515917f, 0.00123677028298367f, 0.00106449804756952f, + 0.000916221960431984f, 0.000788599549519612f, 0.000678753922476738f }; + + + std::vector< REAL > user_info(2 * n_points_per_fit); + for (size_t i = 0; i < n_points_per_fit; i++) + { + user_info[i] = static_cast(timeX[i]); + } + + for (size_t i = n_points_per_fit; i < 2 * n_points_per_fit; i++) + { + user_info[i] = static_cast(Cp[i - n_points_per_fit]); + } + + // size of user info in bytes + size_t const user_info_size = 2 * n_points_per_fit * sizeof(REAL); + + // initialize random number generator + std::mt19937 rng; + rng.seed(time(NULL)); + std::uniform_real_distribution< REAL > uniform_dist(0, 1); + std::normal_distribution< REAL > normal_dist(0, 1); + + // true parameters + std::vector< REAL > true_parameters{ 0.05, 0.03 }; // Ktrans, vp + + // initial parameters (randomized) + std::vector< REAL > initial_parameters(n_fits * n_model_parameters); + for (size_t i = 0; i != n_fits; i++) + { + // random offset + initial_parameters[i * n_model_parameters + 0] = true_parameters[0] * (0.1f + 1.8f * uniform_dist(rng)); + // random slope + initial_parameters[i * n_model_parameters + 1] = true_parameters[0] * (0.1f + 1.8f * uniform_dist(rng)); + } + + // generate data + std::vector< REAL > data(n_points_per_fit * n_fits); + REAL mean_y = 0; + for (size_t i = 0; i != data.size(); i++) + { + size_t j = i / n_points_per_fit; // the fit + size_t k = i % n_points_per_fit; // the position within a fit + REAL x = 0; + for (int n = 1; n < k; n++) { + + REAL spacing = timeX[n] - timeX[n - 1]; + x += (Cp[n - 1] + Cp[n]) / 2 * spacing; + } + REAL y = true_parameters[0] * x + true_parameters[1] * Cp[k]; + //data[i] = y + normal_dist(rng); + //data[i] = y * (0.2f + 1.6f * uniform_dist(rng)); + data[i] = y; + mean_y += y; + //std::cout << data[i] << std::endl; + } + mean_y = mean_y / data.size(); + std::normal_distribution norm_snr(0,mean_y/snr); + for (size_t i = 0; i != data.size(); i++) + { + data[i] = data[i] + norm_snr(rng); + } + + + + // tolerance + REAL const tolerance = 10e-8f; + + // maximum number of iterations + int const max_number_iterations = 200; + + // estimator ID + int const estimator_id = LSE; + + // model ID + int const model_id = PATLAK; + + // parameters to fit (all of them) + std::vector< int > parameters_to_fit(n_model_parameters, 1); + + // output parameters + std::vector< REAL > output_parameters(n_fits * n_model_parameters); + std::vector< int > output_states(n_fits); + std::vector< REAL > output_chi_square(n_fits); + std::vector< int > output_number_iterations(n_fits); + + // call to gpufit (C interface) + int const status = gpufit + ( + n_fits, + n_points_per_fit, + data.data(), + 0, + model_id, + initial_parameters.data(), + tolerance, + max_number_iterations, + parameters_to_fit.data(), + estimator_id, + user_info_size, + reinterpret_cast< char* >( user_info.data() ), + output_parameters.data(), + output_states.data(), + output_chi_square.data(), + output_number_iterations.data() + ); + + + // check status + if (status != ReturnState::OK) + { + throw std::runtime_error(gpufit_get_last_error()); + } + + + // get fit states + std::vector< int > output_states_histogram(5, 0); + for (std::vector< int >::iterator it = output_states.begin(); it != output_states.end(); ++it) + { + output_states_histogram[*it]++; + } + + std::cout << "ratio converged " << (REAL)output_states_histogram[0] / n_fits << "\n"; + std::cout << "ratio max iteration exceeded " << (REAL)output_states_histogram[1] / n_fits << "\n"; + std::cout << "ratio singular hessian " << (REAL)output_states_histogram[2] / n_fits << "\n"; + std::cout << "ratio neg curvature MLE " << (REAL)output_states_histogram[3] / n_fits << "\n"; + std::cout << "ratio gpu not read " << (REAL)output_states_histogram[4] / n_fits << "\n"; + + // compute mean fitted parameters for converged fits + std::vector< REAL > output_parameters_mean(n_model_parameters, 0); + std::vector< REAL > output_parameters_mean_error(n_model_parameters, 0); + for (size_t i = 0; i != n_fits; i++) + { + if (output_states[i] == FitState::CONVERGED) + { + // add Ktrans + output_parameters_mean[0] += output_parameters[i * n_model_parameters + 0]; + // add vp + output_parameters_mean[1] += output_parameters[i * n_model_parameters + 1]; + // add Ktrans + output_parameters_mean_error[0] += abs(output_parameters[i * n_model_parameters + 0]-true_parameters[0]); + // add vp + output_parameters_mean_error[1] += abs(output_parameters[i * n_model_parameters + 1]-true_parameters[1]); + } + } + output_parameters_mean[0] /= output_states_histogram[0]; + output_parameters_mean[1] /= output_states_histogram[0]; + + // compute std of fitted parameters for converged fits + std::vector< REAL > output_parameters_std(n_model_parameters, 0); + for (size_t i = 0; i != n_fits; i++) + { + if (output_states[i] == FitState::CONVERGED) + { + // add squared deviation for Ktrans + output_parameters_std[0] += (output_parameters[i * n_model_parameters + 0] - output_parameters_mean[0]) * (output_parameters[i * n_model_parameters + 0] - output_parameters_mean[0]); + // add squared deviation for vp + output_parameters_std[1] += (output_parameters[i * n_model_parameters + 1] - output_parameters_mean[1]) * (output_parameters[i * n_model_parameters + 1] - output_parameters_mean[1]); + } + } + // divide and take square root + output_parameters_std[0] = sqrt(output_parameters_std[0] / output_states_histogram[0]); + output_parameters_std[1] = sqrt(output_parameters_std[1] / output_states_histogram[0]); + + // print mean and std + std::cout << "Ktrans true " << true_parameters[0] << " mean " << output_parameters_mean[0] << " std " << output_parameters_std[0] << "\n"; + std::cout << "vp true " << true_parameters[1] << " mean " << output_parameters_mean[1] << " std " << output_parameters_std[1] << "\n"; + + // compute mean chi-square for those converged + REAL output_chi_square_mean = 0; + for (size_t i = 0; i != n_fits; i++) + { + if (output_states[i] == FitState::CONVERGED) + { + output_chi_square_mean += output_chi_square[i]; + } + } + output_chi_square_mean /= static_cast(output_states_histogram[0]); + std::cout << "mean chi square " << output_chi_square_mean << "\n"; + + // compute mean number of iterations for those converged + REAL output_number_iterations_mean = 0; + for (size_t i = 0; i != n_fits; i++) + { + if (output_states[i] == FitState::CONVERGED) + { + output_number_iterations_mean += static_cast(output_number_iterations[i]); + } + } + + // normalize + output_number_iterations_mean /= static_cast(output_states_histogram[0]); + std::cout << "mean number of iterations " << output_number_iterations_mean << "\n"; + + // time + //time(&time_end); + time_end = clock(); + double time_taken_sec = double(time_end-time_start)/double(CLOCKS_PER_SEC); + std::cout << "execution time for " << n_fits << " fits was " << time_taken_sec << " seconds\n"; +} + + +int main(int argc, char* argv[]) +{ + std::cout << std::endl << "Beginning Patlak fit..." << std::endl; + patlak_two(); + + std::cout << std::endl << "Patlak fit completed!" << std::endl; + + return 0; +} diff --git a/Gpufit/examples/Simple_Example.cpp b/Gpufit/examples/Simple_Example.cpp index ac91b00e..c91e0121 100644 --- a/Gpufit/examples/Simple_Example.cpp +++ b/Gpufit/examples/Simple_Example.cpp @@ -1,139 +1,50 @@ #include "../gpufit.h" - -#include -#include -#include #include -#include -//using namespace std; +#include -void patlak_two() +void simple_example() { - /* - This example generates test data in form of 10000 one dimensional linear - curves with the size of 20 data points per curve. It is noised by normal - distributed noise. The initial guesses were randomized, within a specified - range of the true value. The LINEAR_1D model is fitted to the test data sets - using the LSE estimator. The optional parameter user_info is used to pass - custom x positions of the data sets. The same x position values are used for - every fit. - - The console output shows - - the ratio of converged fits including ratios of not converged fits for - different reasons, - - the values of the true parameters and the mean values of the fitted - parameters including their standard deviation, - - the mean chi square value - - and the mean number of iterations needed. + This example demonstrates a simple, minimal program containing all + of the required parameters for a call to the Gpufit function. The example + can be built and executed within the project environment. Please note that + this code does not actually do anything other than make a single call to + gpufit(). + In the first section of the code, the *model ID* is set, memory space for + initial parameters and data values is allocated, the *fit tolerance* is set, + the *maximum number of iterations* is set, the *estimator ID* is set, and + the *parameters to fit array* is initialized. Note that in most applications, + the data array will already exist and it will be unnecessary to allocate + additional space for data. In this example, the *parameters to fit* array + is initialized to all ones, indicating that all model parameters should be + adjusted in the fit. */ - // start timer - clock_t time_start, time_end; - time_start = clock(); - - - // number of fits, fit points and parameters - size_t const n_fits = 10000; // 10000? 100? - size_t const n_points_per_fit = 60; - size_t const n_model_parameters = 2; - REAL snr = 0.8; - - // custom x positions for the data points of every fit, stored in user info - // time independent variable, given in minutes - REAL timeX[] = { 0.25, 0.5, 0.75, 1, 1.25, 1.5, 1.75, 2, 2.25, 2.5, 2.75, 3, 3.25, 3.5, 3.75, 4, 4.25, 4.5, 4.75, 5, - 5.25, 5.5, 5.75, 6, 6.25, 6.5, 6.75, 7, 7.25, 7.5, 7.75, 8, 8.25, 8.5, 8.75, 9, 9.25, 9.5, 9.75, 10, - 10.25, 10.5, 10.75, 11, 11.25, 11.5, 11.75, 12, 12.25, 12.5, 12.75, 13, 13.25, 13.5, 13.75, 14, 14.25, 14.5, 14.75, 15 }; - - REAL Cp[] = { 0.0f, 0.0f, 0.0f, 3.01846399851715f, 2.59801604007558f, 2.2361331285733f, 1.92465762011135f, 1.65656816551711f, 1.4258214335524f, - 1.22721588081636f, 1.05627449741415f, 0.909143885218726f, 0.782507393725825f, 0.6735103553914f, 0.579695735090254f, - 0.498948743091769f, 0.429449163006342f, 0.369630320068624f, 0.318143764811612f, 0.273828876023252f, 0.235686697768721f, - 0.20285742070682f, 0.174601000079374f, 0.150280473460109f, 0.12934760220805f, 0.111330512951924f, 0.0958230605172143f, - 0.0824756725126274f, 0.0709874691926393f, 0.0610994809603327f, 0.0525888106179893f, 0.0452636087696102f, 0.0389587491097867f, - 0.033532106110336f, 0.0288613511954976f, 0.0248411951843697f, 0.0213810148391187f, 0.018402810016092f, 0.0158394453694853f, - 0.013633136971665f, 0.0117341497352074f, 0.010099676273659f, 0.00869287192804919f, 0.00748202420651342f, 0.00643983791435146f, - 0.00554281985976681f, 0.00477074926518851f, 0.00410622194607174f, 0.00353425798195557f, 0.00304196403581309f, 0.00261824270962248f, - 0.00225354238438883f, 0.00193964190545541f, 0.00166946525943377f, 0.00143692206515917f, 0.00123677028298367f, 0.00106449804756952f, - 0.000916221960431984f, 0.000788599549519612f, 0.000678753922476738f }; - - - std::vector< REAL > user_info(2 * n_points_per_fit); - for (size_t i = 0; i < n_points_per_fit; i++) - { - user_info[i] = static_cast(timeX[i]); - } - - for (size_t i = n_points_per_fit; i < 2 * n_points_per_fit; i++) - { - user_info[i] = static_cast(Cp[i - n_points_per_fit]); - } - - //std::cout << user_info.size() << std::endl; + /*************** definition of input and output parameters ***************/ - // size of user info in bytes - size_t const user_info_size = 2 * n_points_per_fit * sizeof(REAL); + // number of fits, number of points per fit + size_t const n_fits = 10; + size_t const n_points_per_fit = 10; - // initialize random number generator - std::mt19937 rng; - rng.seed(time(NULL)); - //rng.seed(0); - std::uniform_real_distribution< REAL > uniform_dist(0, 1); - std::normal_distribution< REAL > normal_dist(0, 1); + // model ID and number of model parameters + int const model_id = GAUSS_1D; + size_t const n_model_parameters = 4; - // true parameters - std::vector< REAL > true_parameters{ 0.05, 0.03 }; // Ktrans, vp -- ? - - // initial parameters (randomized) + // initial parameters std::vector< REAL > initial_parameters(n_fits * n_model_parameters); - for (size_t i = 0; i != n_fits; i++) - { - // random offset - initial_parameters[i * n_model_parameters + 0] = true_parameters[0] * (0.1f + 1.8f * uniform_dist(rng)); - // random slope - initial_parameters[i * n_model_parameters + 1] = true_parameters[0] * (0.1f + 1.8f * uniform_dist(rng)); - } - // generate data -- NEEDS TO BE CHECKED WITH PATLAK MODEL + // data std::vector< REAL > data(n_points_per_fit * n_fits); - REAL mean_y = 0; - for (size_t i = 0; i != data.size(); i++) - { - size_t j = i / n_points_per_fit; // the fit - size_t k = i % n_points_per_fit; // the position within a fit - REAL x = 0; - for (int n = 1; n < k; n++) { // possible point of error - // - REAL spacing = timeX[n] - timeX[n - 1]; - x += (Cp[n - 1] + Cp[n]) / 2 * spacing; - } - REAL y = true_parameters[0] * x + true_parameters[1] * Cp[k]; - //data[i] = y + normal_dist(rng); - //data[i] = y * (0.2f + 1.6f * uniform_dist(rng)); - data[i] = y; - mean_y += y; - //std::cout << data[i] << std::endl; - } - mean_y = mean_y / data.size(); - std::normal_distribution norm_snr(0,mean_y/snr); - for (size_t i = 0; i != data.size(); i++) - { - data[i] = data[i] + norm_snr(rng); - } - - // tolerance - REAL const tolerance = 10e-8f; // calculated by Christian to be 10e-15f + REAL const tolerance = 0.001f; // maximum number of iterations - int const max_number_iterations = 200; + int const max_number_iterations = 10; // estimator ID int const estimator_id = LSE; - // model ID - int const model_id = PATLAK; - // parameters to fit (all of them) std::vector< int > parameters_to_fit(n_model_parameters, 1); @@ -143,128 +54,44 @@ void patlak_two() std::vector< REAL > output_chi_square(n_fits); std::vector< int > output_number_iterations(n_fits); - // call to gpufit (C interface) - int const status = gpufit - ( - n_fits, - n_points_per_fit, - data.data(), - 0, - model_id, - initial_parameters.data(), - tolerance, - max_number_iterations, - parameters_to_fit.data(), - estimator_id, - user_info_size, - reinterpret_cast< char* >( user_info.data() ), - output_parameters.data(), - output_states.data(), - output_chi_square.data(), - output_number_iterations.data() - ); + /***************************** call to gpufit ****************************/ + int const status = gpufit + ( + n_fits, + n_points_per_fit, + data.data(), + 0, + model_id, + initial_parameters.data(), + tolerance, + max_number_iterations, + parameters_to_fit.data(), + estimator_id, + 0, + 0, + output_parameters.data(), + output_states.data(), + output_chi_square.data(), + output_number_iterations.data() + ); + + /****************************** status check *****************************/ - // check status if (status != ReturnState::OK) { throw std::runtime_error(gpufit_get_last_error()); } - - - // get fit states - std::vector< int > output_states_histogram(5, 0); - for (std::vector< int >::iterator it = output_states.begin(); it != output_states.end(); ++it) - { - output_states_histogram[*it]++; - } - - std::cout << "ratio converged " << (REAL)output_states_histogram[0] / n_fits << "\n"; - std::cout << "ratio max iteration exceeded " << (REAL)output_states_histogram[1] / n_fits << "\n"; - std::cout << "ratio singular hessian " << (REAL)output_states_histogram[2] / n_fits << "\n"; - std::cout << "ratio neg curvature MLE " << (REAL)output_states_histogram[3] / n_fits << "\n"; - std::cout << "ratio gpu not read " << (REAL)output_states_histogram[4] / n_fits << "\n"; - - // compute mean fitted parameters for converged fits - std::vector< REAL > output_parameters_mean(n_model_parameters, 0); - std::vector< REAL > output_parameters_mean_error(n_model_parameters, 0); - for (size_t i = 0; i != n_fits; i++) - { - if (output_states[i] == FitState::CONVERGED) - { - // add Ktrans - output_parameters_mean[0] += output_parameters[i * n_model_parameters + 0]; - // add vp - output_parameters_mean[1] += output_parameters[i * n_model_parameters + 1]; - // add Ktrans - output_parameters_mean_error[0] += abs(output_parameters[i * n_model_parameters + 0]-true_parameters[0]); - // add vp - output_parameters_mean_error[1] += abs(output_parameters[i * n_model_parameters + 1]-true_parameters[1]); - } - } - output_parameters_mean[0] /= output_states_histogram[0]; - output_parameters_mean[1] /= output_states_histogram[0]; - - // compute std of fitted parameters for converged fits - std::vector< REAL > output_parameters_std(n_model_parameters, 0); - for (size_t i = 0; i != n_fits; i++) - { - if (output_states[i] == FitState::CONVERGED) - { - // add squared deviation for Ktrans - output_parameters_std[0] += (output_parameters[i * n_model_parameters + 0] - output_parameters_mean[0]) * (output_parameters[i * n_model_parameters + 0] - output_parameters_mean[0]); - // add squared deviation for vp - output_parameters_std[1] += (output_parameters[i * n_model_parameters + 1] - output_parameters_mean[1]) * (output_parameters[i * n_model_parameters + 1] - output_parameters_mean[1]); - } - } - // divide and take square root - output_parameters_std[0] = sqrt(output_parameters_std[0] / output_states_histogram[0]); - output_parameters_std[1] = sqrt(output_parameters_std[1] / output_states_histogram[0]); - - // print mean and std - std::cout << "Ktrans true " << true_parameters[0] << " mean " << output_parameters_mean[0] << " std " << output_parameters_std[0] << "\n"; - std::cout << "vp true " << true_parameters[1] << " mean " << output_parameters_mean[1] << " std " << output_parameters_std[1] << "\n"; - - // compute mean chi-square for those converged - REAL output_chi_square_mean = 0; - for (size_t i = 0; i != n_fits; i++) - { - if (output_states[i] == FitState::CONVERGED) - { - output_chi_square_mean += output_chi_square[i]; - } - } - output_chi_square_mean /= static_cast(output_states_histogram[0]); - std::cout << "mean chi square " << output_chi_square_mean << "\n"; - - // compute mean number of iterations for those converged - REAL output_number_iterations_mean = 0; - for (size_t i = 0; i != n_fits; i++) - { - if (output_states[i] == FitState::CONVERGED) - { - output_number_iterations_mean += static_cast(output_number_iterations[i]); - } - } - - // normalize - output_number_iterations_mean /= static_cast(output_states_histogram[0]); - std::cout << "mean number of iterations " << output_number_iterations_mean << "\n"; - - // time - //time(&time_end); - time_end = clock(); - double time_taken_sec = double(time_end-time_start)/double(CLOCKS_PER_SEC); - std::cout << "execution time for " << n_fits << " fits was " << time_taken_sec << " seconds\n"; } -int main(int argc, char* argv[]) +int main(int argc, char *argv[]) { - std::cout << std::endl << "Beginning Patlak fit..." << std::endl; - patlak_two(); - - std::cout << std::endl << "Patlak fit completed!" << std::endl; + simple_example(); + std::cout << std::endl << "Example completed!" << std::endl; + std::cout << "Press ENTER to exit" << std::endl; + std::getchar(); + return 0; } diff --git a/Gpufit/examples/Tissue_Uptake_fitting.cpp b/Gpufit/examples/Tissue_Uptake_fitting.cpp new file mode 100644 index 00000000..7b35fae4 --- /dev/null +++ b/Gpufit/examples/Tissue_Uptake_fitting.cpp @@ -0,0 +1,267 @@ +#include "../gpufit.h" + +#include +#include +#include +#include +#include + +void patlak_two() +{ + + /* + This example generates test data in form of 10000 one dimensional linear + curves with the size of 20 data points per curve. It is noised by normal + distributed noise. The initial guesses were randomized, within a specified + range of the true value. The LINEAR_1D model is fitted to the test data sets + using the LSE estimator. The optional parameter user_info is used to pass + custom x positions of the data sets. The same x position values are used for + every fit. + + The console output shows + - the ratio of converged fits including ratios of not converged fits for + different reasons, + - the values of the true parameters and the mean values of the fitted + parameters including their standard deviation, + - the mean chi square value + - and the mean number of iterations needed. + */ + + // start timer + clock_t time_start, time_end; + time_start = clock(); + + + // number of fits, fit points and parameters + size_t const n_fits = 10000; + size_t const n_points_per_fit = 60; + size_t const n_model_parameters = 2; + REAL snr = 0.8; + + // custom x positions for the data points of every fit, stored in user info + // time independent variable, given in minutes + REAL timeX[] = { 0.25, 0.5, 0.75, 1, 1.25, 1.5, 1.75, 2, 2.25, 2.5, 2.75, 3, 3.25, 3.5, 3.75, 4, 4.25, 4.5, 4.75, 5, + 5.25, 5.5, 5.75, 6, 6.25, 6.5, 6.75, 7, 7.25, 7.5, 7.75, 8, 8.25, 8.5, 8.75, 9, 9.25, 9.5, 9.75, 10, + 10.25, 10.5, 10.75, 11, 11.25, 11.5, 11.75, 12, 12.25, 12.5, 12.75, 13, 13.25, 13.5, 13.75, 14, 14.25, 14.5, 14.75, 15 }; + + // Concentration of plasma (independent), at 1 min based on equation: Cp(t) = 5.5e^(-.6t) + REAL Cp[] = { 0.0f, 0.0f, 0.0f, 3.01846399851715f, 2.59801604007558f, 2.2361331285733f, 1.92465762011135f, 1.65656816551711f, 1.4258214335524f, + 1.22721588081636f, 1.05627449741415f, 0.909143885218726f, 0.782507393725825f, 0.6735103553914f, 0.579695735090254f, + 0.498948743091769f, 0.429449163006342f, 0.369630320068624f, 0.318143764811612f, 0.273828876023252f, 0.235686697768721f, + 0.20285742070682f, 0.174601000079374f, 0.150280473460109f, 0.12934760220805f, 0.111330512951924f, 0.0958230605172143f, + 0.0824756725126274f, 0.0709874691926393f, 0.0610994809603327f, 0.0525888106179893f, 0.0452636087696102f, 0.0389587491097867f, + 0.033532106110336f, 0.0288613511954976f, 0.0248411951843697f, 0.0213810148391187f, 0.018402810016092f, 0.0158394453694853f, + 0.013633136971665f, 0.0117341497352074f, 0.010099676273659f, 0.00869287192804919f, 0.00748202420651342f, 0.00643983791435146f, + 0.00554281985976681f, 0.00477074926518851f, 0.00410622194607174f, 0.00353425798195557f, 0.00304196403581309f, 0.00261824270962248f, + 0.00225354238438883f, 0.00193964190545541f, 0.00166946525943377f, 0.00143692206515917f, 0.00123677028298367f, 0.00106449804756952f, + 0.000916221960431984f, 0.000788599549519612f, 0.000678753922476738f }; + + + std::vector< REAL > user_info(2 * n_points_per_fit); + for (size_t i = 0; i < n_points_per_fit; i++) + { + user_info[i] = static_cast(timeX[i]); + } + + for (size_t i = n_points_per_fit; i < 2 * n_points_per_fit; i++) + { + user_info[i] = static_cast(Cp[i - n_points_per_fit]); + } + + // size of user info in bytes + size_t const user_info_size = 2 * n_points_per_fit * sizeof(REAL); + + // initialize random number generator + std::mt19937 rng; + rng.seed(time(NULL)); + std::uniform_real_distribution< REAL > uniform_dist(0, 1); + std::normal_distribution< REAL > normal_dist(0, 1); + + // true parameters + std::vector< REAL > true_parameters{ 0.05, 0.03 }; // Ktrans, vp + + // initial parameters (randomized) + std::vector< REAL > initial_parameters(n_fits * n_model_parameters); + for (size_t i = 0; i != n_fits; i++) + { + // random offset + initial_parameters[i * n_model_parameters + 0] = true_parameters[0] * (0.1f + 1.8f * uniform_dist(rng)); + // random slope + initial_parameters[i * n_model_parameters + 1] = true_parameters[0] * (0.1f + 1.8f * uniform_dist(rng)); + } + + // generate data + std::vector< REAL > data(n_points_per_fit * n_fits); + REAL mean_y = 0; + for (size_t i = 0; i != data.size(); i++) + { + size_t j = i / n_points_per_fit; // the fit + size_t k = i % n_points_per_fit; // the position within a fit + REAL x = 0; + for (int n = 1; n < k; n++) { + + REAL spacing = timeX[n] - timeX[n - 1]; + x += (Cp[n - 1] + Cp[n]) / 2 * spacing; + } + REAL y = true_parameters[0] * x + true_parameters[1] * Cp[k]; + //data[i] = y + normal_dist(rng); + //data[i] = y * (0.2f + 1.6f * uniform_dist(rng)); + data[i] = y; + mean_y += y; + //std::cout << data[i] << std::endl; + } + mean_y = mean_y / data.size(); + std::normal_distribution norm_snr(0,mean_y/snr); + for (size_t i = 0; i != data.size(); i++) + { + data[i] = data[i] + norm_snr(rng); + } + + + + // tolerance + REAL const tolerance = 10e-8f; + + // maximum number of iterations + int const max_number_iterations = 200; + + // estimator ID + int const estimator_id = LSE; + + // model ID + int const model_id = PATLAK; + + // parameters to fit (all of them) + std::vector< int > parameters_to_fit(n_model_parameters, 1); + + // output parameters + std::vector< REAL > output_parameters(n_fits * n_model_parameters); + std::vector< int > output_states(n_fits); + std::vector< REAL > output_chi_square(n_fits); + std::vector< int > output_number_iterations(n_fits); + + // call to gpufit (C interface) + int const status = gpufit + ( + n_fits, + n_points_per_fit, + data.data(), + 0, + model_id, + initial_parameters.data(), + tolerance, + max_number_iterations, + parameters_to_fit.data(), + estimator_id, + user_info_size, + reinterpret_cast< char* >( user_info.data() ), + output_parameters.data(), + output_states.data(), + output_chi_square.data(), + output_number_iterations.data() + ); + + + // check status + if (status != ReturnState::OK) + { + throw std::runtime_error(gpufit_get_last_error()); + } + + + // get fit states + std::vector< int > output_states_histogram(5, 0); + for (std::vector< int >::iterator it = output_states.begin(); it != output_states.end(); ++it) + { + output_states_histogram[*it]++; + } + + std::cout << "ratio converged " << (REAL)output_states_histogram[0] / n_fits << "\n"; + std::cout << "ratio max iteration exceeded " << (REAL)output_states_histogram[1] / n_fits << "\n"; + std::cout << "ratio singular hessian " << (REAL)output_states_histogram[2] / n_fits << "\n"; + std::cout << "ratio neg curvature MLE " << (REAL)output_states_histogram[3] / n_fits << "\n"; + std::cout << "ratio gpu not read " << (REAL)output_states_histogram[4] / n_fits << "\n"; + + // compute mean fitted parameters for converged fits + std::vector< REAL > output_parameters_mean(n_model_parameters, 0); + std::vector< REAL > output_parameters_mean_error(n_model_parameters, 0); + for (size_t i = 0; i != n_fits; i++) + { + if (output_states[i] == FitState::CONVERGED) + { + // add Ktrans + output_parameters_mean[0] += output_parameters[i * n_model_parameters + 0]; + // add vp + output_parameters_mean[1] += output_parameters[i * n_model_parameters + 1]; + // add Ktrans + output_parameters_mean_error[0] += abs(output_parameters[i * n_model_parameters + 0]-true_parameters[0]); + // add vp + output_parameters_mean_error[1] += abs(output_parameters[i * n_model_parameters + 1]-true_parameters[1]); + } + } + output_parameters_mean[0] /= output_states_histogram[0]; + output_parameters_mean[1] /= output_states_histogram[0]; + + // compute std of fitted parameters for converged fits + std::vector< REAL > output_parameters_std(n_model_parameters, 0); + for (size_t i = 0; i != n_fits; i++) + { + if (output_states[i] == FitState::CONVERGED) + { + // add squared deviation for Ktrans + output_parameters_std[0] += (output_parameters[i * n_model_parameters + 0] - output_parameters_mean[0]) * (output_parameters[i * n_model_parameters + 0] - output_parameters_mean[0]); + // add squared deviation for vp + output_parameters_std[1] += (output_parameters[i * n_model_parameters + 1] - output_parameters_mean[1]) * (output_parameters[i * n_model_parameters + 1] - output_parameters_mean[1]); + } + } + // divide and take square root + output_parameters_std[0] = sqrt(output_parameters_std[0] / output_states_histogram[0]); + output_parameters_std[1] = sqrt(output_parameters_std[1] / output_states_histogram[0]); + + // print mean and std + std::cout << "Ktrans true " << true_parameters[0] << " mean " << output_parameters_mean[0] << " std " << output_parameters_std[0] << "\n"; + std::cout << "vp true " << true_parameters[1] << " mean " << output_parameters_mean[1] << " std " << output_parameters_std[1] << "\n"; + + // compute mean chi-square for those converged + REAL output_chi_square_mean = 0; + for (size_t i = 0; i != n_fits; i++) + { + if (output_states[i] == FitState::CONVERGED) + { + output_chi_square_mean += output_chi_square[i]; + } + } + output_chi_square_mean /= static_cast(output_states_histogram[0]); + std::cout << "mean chi square " << output_chi_square_mean << "\n"; + + // compute mean number of iterations for those converged + REAL output_number_iterations_mean = 0; + for (size_t i = 0; i != n_fits; i++) + { + if (output_states[i] == FitState::CONVERGED) + { + output_number_iterations_mean += static_cast(output_number_iterations[i]); + } + } + + // normalize + output_number_iterations_mean /= static_cast(output_states_histogram[0]); + std::cout << "mean number of iterations " << output_number_iterations_mean << "\n"; + + // time + //time(&time_end); + time_end = clock(); + double time_taken_sec = double(time_end-time_start)/double(CLOCKS_PER_SEC); + std::cout << "execution time for " << n_fits << " fits was " << time_taken_sec << " seconds\n"; +} + + +int main(int argc, char* argv[]) +{ + std::cout << std::endl << "Beginning Patlak fit..." << std::endl; + patlak_two(); + + std::cout << std::endl << "Patlak fit completed!" << std::endl; + + return 0; +} diff --git a/Gpufit/examples/Tofts_Extended_fitting.cpp b/Gpufit/examples/Tofts_Extended_fitting.cpp new file mode 100644 index 00000000..7b35fae4 --- /dev/null +++ b/Gpufit/examples/Tofts_Extended_fitting.cpp @@ -0,0 +1,267 @@ +#include "../gpufit.h" + +#include +#include +#include +#include +#include + +void patlak_two() +{ + + /* + This example generates test data in form of 10000 one dimensional linear + curves with the size of 20 data points per curve. It is noised by normal + distributed noise. The initial guesses were randomized, within a specified + range of the true value. The LINEAR_1D model is fitted to the test data sets + using the LSE estimator. The optional parameter user_info is used to pass + custom x positions of the data sets. The same x position values are used for + every fit. + + The console output shows + - the ratio of converged fits including ratios of not converged fits for + different reasons, + - the values of the true parameters and the mean values of the fitted + parameters including their standard deviation, + - the mean chi square value + - and the mean number of iterations needed. + */ + + // start timer + clock_t time_start, time_end; + time_start = clock(); + + + // number of fits, fit points and parameters + size_t const n_fits = 10000; + size_t const n_points_per_fit = 60; + size_t const n_model_parameters = 2; + REAL snr = 0.8; + + // custom x positions for the data points of every fit, stored in user info + // time independent variable, given in minutes + REAL timeX[] = { 0.25, 0.5, 0.75, 1, 1.25, 1.5, 1.75, 2, 2.25, 2.5, 2.75, 3, 3.25, 3.5, 3.75, 4, 4.25, 4.5, 4.75, 5, + 5.25, 5.5, 5.75, 6, 6.25, 6.5, 6.75, 7, 7.25, 7.5, 7.75, 8, 8.25, 8.5, 8.75, 9, 9.25, 9.5, 9.75, 10, + 10.25, 10.5, 10.75, 11, 11.25, 11.5, 11.75, 12, 12.25, 12.5, 12.75, 13, 13.25, 13.5, 13.75, 14, 14.25, 14.5, 14.75, 15 }; + + // Concentration of plasma (independent), at 1 min based on equation: Cp(t) = 5.5e^(-.6t) + REAL Cp[] = { 0.0f, 0.0f, 0.0f, 3.01846399851715f, 2.59801604007558f, 2.2361331285733f, 1.92465762011135f, 1.65656816551711f, 1.4258214335524f, + 1.22721588081636f, 1.05627449741415f, 0.909143885218726f, 0.782507393725825f, 0.6735103553914f, 0.579695735090254f, + 0.498948743091769f, 0.429449163006342f, 0.369630320068624f, 0.318143764811612f, 0.273828876023252f, 0.235686697768721f, + 0.20285742070682f, 0.174601000079374f, 0.150280473460109f, 0.12934760220805f, 0.111330512951924f, 0.0958230605172143f, + 0.0824756725126274f, 0.0709874691926393f, 0.0610994809603327f, 0.0525888106179893f, 0.0452636087696102f, 0.0389587491097867f, + 0.033532106110336f, 0.0288613511954976f, 0.0248411951843697f, 0.0213810148391187f, 0.018402810016092f, 0.0158394453694853f, + 0.013633136971665f, 0.0117341497352074f, 0.010099676273659f, 0.00869287192804919f, 0.00748202420651342f, 0.00643983791435146f, + 0.00554281985976681f, 0.00477074926518851f, 0.00410622194607174f, 0.00353425798195557f, 0.00304196403581309f, 0.00261824270962248f, + 0.00225354238438883f, 0.00193964190545541f, 0.00166946525943377f, 0.00143692206515917f, 0.00123677028298367f, 0.00106449804756952f, + 0.000916221960431984f, 0.000788599549519612f, 0.000678753922476738f }; + + + std::vector< REAL > user_info(2 * n_points_per_fit); + for (size_t i = 0; i < n_points_per_fit; i++) + { + user_info[i] = static_cast(timeX[i]); + } + + for (size_t i = n_points_per_fit; i < 2 * n_points_per_fit; i++) + { + user_info[i] = static_cast(Cp[i - n_points_per_fit]); + } + + // size of user info in bytes + size_t const user_info_size = 2 * n_points_per_fit * sizeof(REAL); + + // initialize random number generator + std::mt19937 rng; + rng.seed(time(NULL)); + std::uniform_real_distribution< REAL > uniform_dist(0, 1); + std::normal_distribution< REAL > normal_dist(0, 1); + + // true parameters + std::vector< REAL > true_parameters{ 0.05, 0.03 }; // Ktrans, vp + + // initial parameters (randomized) + std::vector< REAL > initial_parameters(n_fits * n_model_parameters); + for (size_t i = 0; i != n_fits; i++) + { + // random offset + initial_parameters[i * n_model_parameters + 0] = true_parameters[0] * (0.1f + 1.8f * uniform_dist(rng)); + // random slope + initial_parameters[i * n_model_parameters + 1] = true_parameters[0] * (0.1f + 1.8f * uniform_dist(rng)); + } + + // generate data + std::vector< REAL > data(n_points_per_fit * n_fits); + REAL mean_y = 0; + for (size_t i = 0; i != data.size(); i++) + { + size_t j = i / n_points_per_fit; // the fit + size_t k = i % n_points_per_fit; // the position within a fit + REAL x = 0; + for (int n = 1; n < k; n++) { + + REAL spacing = timeX[n] - timeX[n - 1]; + x += (Cp[n - 1] + Cp[n]) / 2 * spacing; + } + REAL y = true_parameters[0] * x + true_parameters[1] * Cp[k]; + //data[i] = y + normal_dist(rng); + //data[i] = y * (0.2f + 1.6f * uniform_dist(rng)); + data[i] = y; + mean_y += y; + //std::cout << data[i] << std::endl; + } + mean_y = mean_y / data.size(); + std::normal_distribution norm_snr(0,mean_y/snr); + for (size_t i = 0; i != data.size(); i++) + { + data[i] = data[i] + norm_snr(rng); + } + + + + // tolerance + REAL const tolerance = 10e-8f; + + // maximum number of iterations + int const max_number_iterations = 200; + + // estimator ID + int const estimator_id = LSE; + + // model ID + int const model_id = PATLAK; + + // parameters to fit (all of them) + std::vector< int > parameters_to_fit(n_model_parameters, 1); + + // output parameters + std::vector< REAL > output_parameters(n_fits * n_model_parameters); + std::vector< int > output_states(n_fits); + std::vector< REAL > output_chi_square(n_fits); + std::vector< int > output_number_iterations(n_fits); + + // call to gpufit (C interface) + int const status = gpufit + ( + n_fits, + n_points_per_fit, + data.data(), + 0, + model_id, + initial_parameters.data(), + tolerance, + max_number_iterations, + parameters_to_fit.data(), + estimator_id, + user_info_size, + reinterpret_cast< char* >( user_info.data() ), + output_parameters.data(), + output_states.data(), + output_chi_square.data(), + output_number_iterations.data() + ); + + + // check status + if (status != ReturnState::OK) + { + throw std::runtime_error(gpufit_get_last_error()); + } + + + // get fit states + std::vector< int > output_states_histogram(5, 0); + for (std::vector< int >::iterator it = output_states.begin(); it != output_states.end(); ++it) + { + output_states_histogram[*it]++; + } + + std::cout << "ratio converged " << (REAL)output_states_histogram[0] / n_fits << "\n"; + std::cout << "ratio max iteration exceeded " << (REAL)output_states_histogram[1] / n_fits << "\n"; + std::cout << "ratio singular hessian " << (REAL)output_states_histogram[2] / n_fits << "\n"; + std::cout << "ratio neg curvature MLE " << (REAL)output_states_histogram[3] / n_fits << "\n"; + std::cout << "ratio gpu not read " << (REAL)output_states_histogram[4] / n_fits << "\n"; + + // compute mean fitted parameters for converged fits + std::vector< REAL > output_parameters_mean(n_model_parameters, 0); + std::vector< REAL > output_parameters_mean_error(n_model_parameters, 0); + for (size_t i = 0; i != n_fits; i++) + { + if (output_states[i] == FitState::CONVERGED) + { + // add Ktrans + output_parameters_mean[0] += output_parameters[i * n_model_parameters + 0]; + // add vp + output_parameters_mean[1] += output_parameters[i * n_model_parameters + 1]; + // add Ktrans + output_parameters_mean_error[0] += abs(output_parameters[i * n_model_parameters + 0]-true_parameters[0]); + // add vp + output_parameters_mean_error[1] += abs(output_parameters[i * n_model_parameters + 1]-true_parameters[1]); + } + } + output_parameters_mean[0] /= output_states_histogram[0]; + output_parameters_mean[1] /= output_states_histogram[0]; + + // compute std of fitted parameters for converged fits + std::vector< REAL > output_parameters_std(n_model_parameters, 0); + for (size_t i = 0; i != n_fits; i++) + { + if (output_states[i] == FitState::CONVERGED) + { + // add squared deviation for Ktrans + output_parameters_std[0] += (output_parameters[i * n_model_parameters + 0] - output_parameters_mean[0]) * (output_parameters[i * n_model_parameters + 0] - output_parameters_mean[0]); + // add squared deviation for vp + output_parameters_std[1] += (output_parameters[i * n_model_parameters + 1] - output_parameters_mean[1]) * (output_parameters[i * n_model_parameters + 1] - output_parameters_mean[1]); + } + } + // divide and take square root + output_parameters_std[0] = sqrt(output_parameters_std[0] / output_states_histogram[0]); + output_parameters_std[1] = sqrt(output_parameters_std[1] / output_states_histogram[0]); + + // print mean and std + std::cout << "Ktrans true " << true_parameters[0] << " mean " << output_parameters_mean[0] << " std " << output_parameters_std[0] << "\n"; + std::cout << "vp true " << true_parameters[1] << " mean " << output_parameters_mean[1] << " std " << output_parameters_std[1] << "\n"; + + // compute mean chi-square for those converged + REAL output_chi_square_mean = 0; + for (size_t i = 0; i != n_fits; i++) + { + if (output_states[i] == FitState::CONVERGED) + { + output_chi_square_mean += output_chi_square[i]; + } + } + output_chi_square_mean /= static_cast(output_states_histogram[0]); + std::cout << "mean chi square " << output_chi_square_mean << "\n"; + + // compute mean number of iterations for those converged + REAL output_number_iterations_mean = 0; + for (size_t i = 0; i != n_fits; i++) + { + if (output_states[i] == FitState::CONVERGED) + { + output_number_iterations_mean += static_cast(output_number_iterations[i]); + } + } + + // normalize + output_number_iterations_mean /= static_cast(output_states_histogram[0]); + std::cout << "mean number of iterations " << output_number_iterations_mean << "\n"; + + // time + //time(&time_end); + time_end = clock(); + double time_taken_sec = double(time_end-time_start)/double(CLOCKS_PER_SEC); + std::cout << "execution time for " << n_fits << " fits was " << time_taken_sec << " seconds\n"; +} + + +int main(int argc, char* argv[]) +{ + std::cout << std::endl << "Beginning Patlak fit..." << std::endl; + patlak_two(); + + std::cout << std::endl << "Patlak fit completed!" << std::endl; + + return 0; +} diff --git a/Gpufit/examples/Tofts_fitting.cpp b/Gpufit/examples/Tofts_fitting.cpp new file mode 100644 index 00000000..7b35fae4 --- /dev/null +++ b/Gpufit/examples/Tofts_fitting.cpp @@ -0,0 +1,267 @@ +#include "../gpufit.h" + +#include +#include +#include +#include +#include + +void patlak_two() +{ + + /* + This example generates test data in form of 10000 one dimensional linear + curves with the size of 20 data points per curve. It is noised by normal + distributed noise. The initial guesses were randomized, within a specified + range of the true value. The LINEAR_1D model is fitted to the test data sets + using the LSE estimator. The optional parameter user_info is used to pass + custom x positions of the data sets. The same x position values are used for + every fit. + + The console output shows + - the ratio of converged fits including ratios of not converged fits for + different reasons, + - the values of the true parameters and the mean values of the fitted + parameters including their standard deviation, + - the mean chi square value + - and the mean number of iterations needed. + */ + + // start timer + clock_t time_start, time_end; + time_start = clock(); + + + // number of fits, fit points and parameters + size_t const n_fits = 10000; + size_t const n_points_per_fit = 60; + size_t const n_model_parameters = 2; + REAL snr = 0.8; + + // custom x positions for the data points of every fit, stored in user info + // time independent variable, given in minutes + REAL timeX[] = { 0.25, 0.5, 0.75, 1, 1.25, 1.5, 1.75, 2, 2.25, 2.5, 2.75, 3, 3.25, 3.5, 3.75, 4, 4.25, 4.5, 4.75, 5, + 5.25, 5.5, 5.75, 6, 6.25, 6.5, 6.75, 7, 7.25, 7.5, 7.75, 8, 8.25, 8.5, 8.75, 9, 9.25, 9.5, 9.75, 10, + 10.25, 10.5, 10.75, 11, 11.25, 11.5, 11.75, 12, 12.25, 12.5, 12.75, 13, 13.25, 13.5, 13.75, 14, 14.25, 14.5, 14.75, 15 }; + + // Concentration of plasma (independent), at 1 min based on equation: Cp(t) = 5.5e^(-.6t) + REAL Cp[] = { 0.0f, 0.0f, 0.0f, 3.01846399851715f, 2.59801604007558f, 2.2361331285733f, 1.92465762011135f, 1.65656816551711f, 1.4258214335524f, + 1.22721588081636f, 1.05627449741415f, 0.909143885218726f, 0.782507393725825f, 0.6735103553914f, 0.579695735090254f, + 0.498948743091769f, 0.429449163006342f, 0.369630320068624f, 0.318143764811612f, 0.273828876023252f, 0.235686697768721f, + 0.20285742070682f, 0.174601000079374f, 0.150280473460109f, 0.12934760220805f, 0.111330512951924f, 0.0958230605172143f, + 0.0824756725126274f, 0.0709874691926393f, 0.0610994809603327f, 0.0525888106179893f, 0.0452636087696102f, 0.0389587491097867f, + 0.033532106110336f, 0.0288613511954976f, 0.0248411951843697f, 0.0213810148391187f, 0.018402810016092f, 0.0158394453694853f, + 0.013633136971665f, 0.0117341497352074f, 0.010099676273659f, 0.00869287192804919f, 0.00748202420651342f, 0.00643983791435146f, + 0.00554281985976681f, 0.00477074926518851f, 0.00410622194607174f, 0.00353425798195557f, 0.00304196403581309f, 0.00261824270962248f, + 0.00225354238438883f, 0.00193964190545541f, 0.00166946525943377f, 0.00143692206515917f, 0.00123677028298367f, 0.00106449804756952f, + 0.000916221960431984f, 0.000788599549519612f, 0.000678753922476738f }; + + + std::vector< REAL > user_info(2 * n_points_per_fit); + for (size_t i = 0; i < n_points_per_fit; i++) + { + user_info[i] = static_cast(timeX[i]); + } + + for (size_t i = n_points_per_fit; i < 2 * n_points_per_fit; i++) + { + user_info[i] = static_cast(Cp[i - n_points_per_fit]); + } + + // size of user info in bytes + size_t const user_info_size = 2 * n_points_per_fit * sizeof(REAL); + + // initialize random number generator + std::mt19937 rng; + rng.seed(time(NULL)); + std::uniform_real_distribution< REAL > uniform_dist(0, 1); + std::normal_distribution< REAL > normal_dist(0, 1); + + // true parameters + std::vector< REAL > true_parameters{ 0.05, 0.03 }; // Ktrans, vp + + // initial parameters (randomized) + std::vector< REAL > initial_parameters(n_fits * n_model_parameters); + for (size_t i = 0; i != n_fits; i++) + { + // random offset + initial_parameters[i * n_model_parameters + 0] = true_parameters[0] * (0.1f + 1.8f * uniform_dist(rng)); + // random slope + initial_parameters[i * n_model_parameters + 1] = true_parameters[0] * (0.1f + 1.8f * uniform_dist(rng)); + } + + // generate data + std::vector< REAL > data(n_points_per_fit * n_fits); + REAL mean_y = 0; + for (size_t i = 0; i != data.size(); i++) + { + size_t j = i / n_points_per_fit; // the fit + size_t k = i % n_points_per_fit; // the position within a fit + REAL x = 0; + for (int n = 1; n < k; n++) { + + REAL spacing = timeX[n] - timeX[n - 1]; + x += (Cp[n - 1] + Cp[n]) / 2 * spacing; + } + REAL y = true_parameters[0] * x + true_parameters[1] * Cp[k]; + //data[i] = y + normal_dist(rng); + //data[i] = y * (0.2f + 1.6f * uniform_dist(rng)); + data[i] = y; + mean_y += y; + //std::cout << data[i] << std::endl; + } + mean_y = mean_y / data.size(); + std::normal_distribution norm_snr(0,mean_y/snr); + for (size_t i = 0; i != data.size(); i++) + { + data[i] = data[i] + norm_snr(rng); + } + + + + // tolerance + REAL const tolerance = 10e-8f; + + // maximum number of iterations + int const max_number_iterations = 200; + + // estimator ID + int const estimator_id = LSE; + + // model ID + int const model_id = PATLAK; + + // parameters to fit (all of them) + std::vector< int > parameters_to_fit(n_model_parameters, 1); + + // output parameters + std::vector< REAL > output_parameters(n_fits * n_model_parameters); + std::vector< int > output_states(n_fits); + std::vector< REAL > output_chi_square(n_fits); + std::vector< int > output_number_iterations(n_fits); + + // call to gpufit (C interface) + int const status = gpufit + ( + n_fits, + n_points_per_fit, + data.data(), + 0, + model_id, + initial_parameters.data(), + tolerance, + max_number_iterations, + parameters_to_fit.data(), + estimator_id, + user_info_size, + reinterpret_cast< char* >( user_info.data() ), + output_parameters.data(), + output_states.data(), + output_chi_square.data(), + output_number_iterations.data() + ); + + + // check status + if (status != ReturnState::OK) + { + throw std::runtime_error(gpufit_get_last_error()); + } + + + // get fit states + std::vector< int > output_states_histogram(5, 0); + for (std::vector< int >::iterator it = output_states.begin(); it != output_states.end(); ++it) + { + output_states_histogram[*it]++; + } + + std::cout << "ratio converged " << (REAL)output_states_histogram[0] / n_fits << "\n"; + std::cout << "ratio max iteration exceeded " << (REAL)output_states_histogram[1] / n_fits << "\n"; + std::cout << "ratio singular hessian " << (REAL)output_states_histogram[2] / n_fits << "\n"; + std::cout << "ratio neg curvature MLE " << (REAL)output_states_histogram[3] / n_fits << "\n"; + std::cout << "ratio gpu not read " << (REAL)output_states_histogram[4] / n_fits << "\n"; + + // compute mean fitted parameters for converged fits + std::vector< REAL > output_parameters_mean(n_model_parameters, 0); + std::vector< REAL > output_parameters_mean_error(n_model_parameters, 0); + for (size_t i = 0; i != n_fits; i++) + { + if (output_states[i] == FitState::CONVERGED) + { + // add Ktrans + output_parameters_mean[0] += output_parameters[i * n_model_parameters + 0]; + // add vp + output_parameters_mean[1] += output_parameters[i * n_model_parameters + 1]; + // add Ktrans + output_parameters_mean_error[0] += abs(output_parameters[i * n_model_parameters + 0]-true_parameters[0]); + // add vp + output_parameters_mean_error[1] += abs(output_parameters[i * n_model_parameters + 1]-true_parameters[1]); + } + } + output_parameters_mean[0] /= output_states_histogram[0]; + output_parameters_mean[1] /= output_states_histogram[0]; + + // compute std of fitted parameters for converged fits + std::vector< REAL > output_parameters_std(n_model_parameters, 0); + for (size_t i = 0; i != n_fits; i++) + { + if (output_states[i] == FitState::CONVERGED) + { + // add squared deviation for Ktrans + output_parameters_std[0] += (output_parameters[i * n_model_parameters + 0] - output_parameters_mean[0]) * (output_parameters[i * n_model_parameters + 0] - output_parameters_mean[0]); + // add squared deviation for vp + output_parameters_std[1] += (output_parameters[i * n_model_parameters + 1] - output_parameters_mean[1]) * (output_parameters[i * n_model_parameters + 1] - output_parameters_mean[1]); + } + } + // divide and take square root + output_parameters_std[0] = sqrt(output_parameters_std[0] / output_states_histogram[0]); + output_parameters_std[1] = sqrt(output_parameters_std[1] / output_states_histogram[0]); + + // print mean and std + std::cout << "Ktrans true " << true_parameters[0] << " mean " << output_parameters_mean[0] << " std " << output_parameters_std[0] << "\n"; + std::cout << "vp true " << true_parameters[1] << " mean " << output_parameters_mean[1] << " std " << output_parameters_std[1] << "\n"; + + // compute mean chi-square for those converged + REAL output_chi_square_mean = 0; + for (size_t i = 0; i != n_fits; i++) + { + if (output_states[i] == FitState::CONVERGED) + { + output_chi_square_mean += output_chi_square[i]; + } + } + output_chi_square_mean /= static_cast(output_states_histogram[0]); + std::cout << "mean chi square " << output_chi_square_mean << "\n"; + + // compute mean number of iterations for those converged + REAL output_number_iterations_mean = 0; + for (size_t i = 0; i != n_fits; i++) + { + if (output_states[i] == FitState::CONVERGED) + { + output_number_iterations_mean += static_cast(output_number_iterations[i]); + } + } + + // normalize + output_number_iterations_mean /= static_cast(output_states_histogram[0]); + std::cout << "mean number of iterations " << output_number_iterations_mean << "\n"; + + // time + //time(&time_end); + time_end = clock(); + double time_taken_sec = double(time_end-time_start)/double(CLOCKS_PER_SEC); + std::cout << "execution time for " << n_fits << " fits was " << time_taken_sec << " seconds\n"; +} + + +int main(int argc, char* argv[]) +{ + std::cout << std::endl << "Beginning Patlak fit..." << std::endl; + patlak_two(); + + std::cout << std::endl << "Patlak fit completed!" << std::endl; + + return 0; +} diff --git a/Gpufit/models/models.cuh b/Gpufit/models/models.cuh index ef45e96b..338b3ec5 100644 --- a/Gpufit/models/models.cuh +++ b/Gpufit/models/models.cuh @@ -10,6 +10,10 @@ #include "fletcher_powell_helix.cuh" #include "brown_dennis.cuh" #include "patlak.cuh" +#include "tofts.cuh" +#include "tofts_extended.cuh" +#include "tissue_uptake.cuh" +#include "two-compartment_exchange.cuh" __device__ void calculate_model( ModelID const model_id, @@ -53,6 +57,18 @@ __device__ void calculate_model( case PATLAK: calculate_patlak(parameters, n_fits, n_points, value, derivative, point_index, fit_index, chunk_index, user_info, user_info_size); break; + case TOFTS: + calculate_tofts(parameters, n_fits, n_points, value, derivative, point_index, fit_index, chunk_index, user_info, user_info_size); + break; + case TOFTS_EXTENDED: + calculate_tofts_extended(parameters, n_fits, n_points, value, derivative, point_index, fit_index, chunk_index, user_info, user_info_size); + break; + case TISSUE_UPTAKE: + calculate_tissue_uptake(parameters, n_fits, n_points, value, derivative, point_index, fit_index, chunk_index, user_info, user_info_size); + break; + case TWO_COMPARTMENT_EXCHANGE: + calculate_two_compartment_exchange(parameters, n_fits, n_points, value, derivative, point_index, fit_index, chunk_index, user_info, user_info_size); + break; default: break; } @@ -62,16 +78,20 @@ void configure_model(ModelID const model_id, int & n_parameters, int & n_dimensi { switch (model_id) { - case GAUSS_1D: n_parameters = 4; n_dimensions = 1; break; - case GAUSS_2D: n_parameters = 5; n_dimensions = 2; break; - case GAUSS_2D_ELLIPTIC: n_parameters = 6; n_dimensions = 2; break; - case GAUSS_2D_ROTATED: n_parameters = 7; n_dimensions = 2; break; - case CAUCHY_2D_ELLIPTIC: n_parameters = 6; n_dimensions = 2; break; - case LINEAR_1D: n_parameters = 2; n_dimensions = 1; break; - case FLETCHER_POWELL_HELIX: n_parameters = 3; n_dimensions = 1; break; - case BROWN_DENNIS: n_parameters = 4; n_dimensions = 1; break; - case PATLAK: n_parameters = 2; n_dimensions = 1; break; - default: break; + case GAUSS_1D: n_parameters = 4; n_dimensions = 1; break; + case GAUSS_2D: n_parameters = 5; n_dimensions = 2; break; + case GAUSS_2D_ELLIPTIC: n_parameters = 6; n_dimensions = 2; break; + case GAUSS_2D_ROTATED: n_parameters = 7; n_dimensions = 2; break; + case CAUCHY_2D_ELLIPTIC: n_parameters = 6; n_dimensions = 2; break; + case LINEAR_1D: n_parameters = 2; n_dimensions = 1; break; + case FLETCHER_POWELL_HELIX: n_parameters = 3; n_dimensions = 1; break; + case BROWN_DENNIS: n_parameters = 4; n_dimensions = 1; break; + case PATLAK: n_parameters = 2; n_dimensions = 1; break; + case TOFTS: n_parameters = 2; n_dimensions = 1; break; + case TOFTS_EXTENDED: n_parameters = 3; n_dimensions = 1; break; + case TISSUE_UPTAKE: n_parameters = 3; n_dimensions = 1; break; + case TWO_COMPARTMENT_EXCHANGE: n_parameters = 4; n_dimensions = 1; break; + default: break; } } diff --git a/Gpufit/models/patlak.cuh b/Gpufit/models/patlak.cuh index ba6f3acc..f42bc19b 100644 --- a/Gpufit/models/patlak.cuh +++ b/Gpufit/models/patlak.cuh @@ -7,10 +7,10 @@ __device__ void calculate_patlak ( // function name int const n_points, REAL * value, REAL * derivative, - int const point_index, // k + int const point_index, int const fit_index, int const chunk_index, - char * user_info, // contains time and Cp values in 1 dimensional array + char * user_info, // contains time and Cp values in a 1 dimensional array std::size_t const user_info_size) { // indices @@ -29,9 +29,12 @@ __device__ void calculate_patlak ( // function name convCp += (Cp[i - 1] + Cp[i]) / 2 * spacing; } - value[point_index] = parameters[0] * convCp + parameters[1] * Cp[point_index]; // formula calculating fit model values +<<<<<<< HEAD // C(t) = Ktrans * trapz(Cp(k)) + vp * Cp(k) +======= + // C(t) = Ktrans * trapz(Cp(k)) + vp * Cp(k) +>>>>>>> 221fb4c80d94d67329e96357f7171552fba930ea /////////////////////////// derivative /////////////////////////// REAL * current_derivative = derivative + point_index; diff --git a/Gpufit/models/tissue_uptake.cuh b/Gpufit/models/tissue_uptake.cuh new file mode 100644 index 00000000..62ed62b7 --- /dev/null +++ b/Gpufit/models/tissue_uptake.cuh @@ -0,0 +1,45 @@ +#ifndef GPUFIT_TISSUE_UPTAKE_CUH_INCLUDED +#define GPUFIT_TISSUE_UPTAKE_CUH_INCLUDED + +__device__ void calculate_tissue_uptake ( // function name + REAL const * parameters, + int const n_fits, + int const n_points, + REAL * value, + REAL * derivative, + int const point_index, + int const fit_index, + int const chunk_index, + char * user_info, // contains time and Cp values in 1 dimensional array + std::size_t const user_info_size) +{ + // indices + REAL* user_info_float = (REAL*)user_info; + + ///////////////////////////// value ////////////////////////////// + + // split user_info array into time and Cp + REAL* T = user_info_float; + REAL* Cp = user_info_float + n_points; + + // integral (trapezoidal rule) + REAL convCp = 0; + for (int i = 1; i < point_index; i++) { + REAL spacing = T[i] - T[i - 1]; + convCp += (Cp[i - 1] + Cp[i]) / 2 * spacing; + } + + value[point_index] = parameters[0] * convCp + parameters[1] * Cp[point_index]; // formula calculating fit model values + // C(t) = Ktrans * trapz(Cp(k)) + vp * Cp(k) + + /////////////////////////// derivative /////////////////////////// + REAL * current_derivative = derivative + point_index; + + current_derivative[0 * n_points] = convCp; // formula calculating derivative values with respect to parameters[0] (Ktrans) + current_derivative[1 * n_points] = Cp[point_index]; // formula calculating derivative values with respect to parameters[1] (vp) + + // deallocate pointers + delete T; + delete Cp; +} +#endif diff --git a/Gpufit/models/tofts.cuh b/Gpufit/models/tofts.cuh new file mode 100644 index 00000000..9cbc31db --- /dev/null +++ b/Gpufit/models/tofts.cuh @@ -0,0 +1,44 @@ +#ifndef GPUFIT_TOFTS_CUH_INCLUDED +#define GPUFIT_TOFTS_CUH_INCLUDED + +__device__ void calculate_tofts ( // function name + REAL const * parameters, + int const n_fits, + int const n_points, + REAL * value, + REAL * derivative, + int const point_index, + int const fit_index, + int const chunk_index, + char * user_info, // contains time and Cp values in 1 dimensional array + std::size_t const user_info_size) +{ + // indices + REAL* user_info_float = (REAL*)user_info; + + ///////////////////////////// value ////////////////////////////// + + // split user_info array into time and Cp + REAL* T = user_info_float; + REAL* Cp = user_info_float + n_points; + + // integral (trapezoidal rule) + REAL convFunc = 0; + for (int i = 1; i < point_index; i++) { + + } + + value[point_index] = parameters[0] * Cp[point_index] * convFunc; // formula calculating fit model values + // C(t) = Ktrans * Cp * trapz(e^(-Ktrans*t/ve) + + /////////////////////////// derivative /////////////////////////// + REAL * current_derivative = derivative + point_index; + + current_derivative[0 * n_points] = Cp[point_index]; // formula calculating derivative values with respect to parameters[0] (Ktrans) + current_derivative[1 * n_points] = Cp[point_index]; // formula calculating derivative values with respect to parameters[1] (ve) + + // deallocate pointers + delete T; + delete Cp; +} +#endif diff --git a/Gpufit/models/tofts_extended.cuh b/Gpufit/models/tofts_extended.cuh new file mode 100644 index 00000000..7beea4d3 --- /dev/null +++ b/Gpufit/models/tofts_extended.cuh @@ -0,0 +1,45 @@ +#ifndef GPUFIT_TOFTS_EXTENDED_CUH_INCLUDED +#define GPUFIT_TOFTS_EXTENDED_CUH_INCLUDED + +__device__ void calculate_tofts_extended ( // function name + REAL const * parameters, + int const n_fits, + int const n_points, + REAL * value, + REAL * derivative, + int const point_index, + int const fit_index, + int const chunk_index, + char * user_info, // contains time and Cp values in 1 dimensional array + std::size_t const user_info_size) +{ + // indices + REAL* user_info_float = (REAL*)user_info; + + ///////////////////////////// value ////////////////////////////// + + // split user_info array into time and Cp + REAL* T = user_info_float; + REAL* Cp = user_info_float + n_points; + + // integral (trapezoidal rule) + REAL convCp = 0; + for (int i = 1; i < point_index; i++) { + REAL spacing = T[i] - T[i - 1]; + convCp += (Cp[i - 1] + Cp[i]) / 2 * spacing; + } + + value[point_index] = parameters[0] * convCp + parameters[1] * Cp[point_index]; // formula calculating fit model values + // C(t) = Ktrans * trapz(Cp(k)) + vp * Cp(k) + + /////////////////////////// derivative /////////////////////////// + REAL * current_derivative = derivative + point_index; + + current_derivative[0 * n_points] = convCp; // formula calculating derivative values with respect to parameters[0] (Ktrans) + current_derivative[1 * n_points] = Cp[point_index]; // formula calculating derivative values with respect to parameters[1] (vp) + + // deallocate pointers + delete T; + delete Cp; +} +#endif diff --git a/Gpufit/models/two-compartment_exchange.cuh b/Gpufit/models/two-compartment_exchange.cuh new file mode 100644 index 00000000..3cbc9efd --- /dev/null +++ b/Gpufit/models/two-compartment_exchange.cuh @@ -0,0 +1,45 @@ +#ifndef GPUFIT_TWO_COMPARTMENT_EXCHANGE_CUH_INCLUDED +#define GPUFIT_TWO_COMPARTMENT_EXCHANGE_CUH_INCLUDED + +__device__ void calculate_two_compartment_exchange ( // function name + REAL const * parameters, + int const n_fits, + int const n_points, + REAL * value, + REAL * derivative, + int const point_index, + int const fit_index, + int const chunk_index, + char * user_info, // contains time and Cp values in 1 dimensional array + std::size_t const user_info_size) +{ + // indices + REAL* user_info_float = (REAL*)user_info; + + ///////////////////////////// value ////////////////////////////// + + // split user_info array into time and Cp + REAL* T = user_info_float; + REAL* Cp = user_info_float + n_points; + + // integral (trapezoidal rule) + REAL convCp = 0; + for (int i = 1; i < point_index; i++) { + REAL spacing = T[i] - T[i - 1]; + convCp += (Cp[i - 1] + Cp[i]) / 2 * spacing; + } + + value[point_index] = parameters[0] * convCp + parameters[1] * Cp[point_index]; // formula calculating fit model values + // C(t) = Ktrans * trapz(Cp(k)) + vp * Cp(k) + + /////////////////////////// derivative /////////////////////////// + REAL * current_derivative = derivative + point_index; + + current_derivative[0 * n_points] = convCp; // formula calculating derivative values with respect to parameters[0] (Ktrans) + current_derivative[1 * n_points] = Cp[point_index]; // formula calculating derivative values with respect to parameters[1] (vp) + + // deallocate pointers + delete T; + delete Cp; +} +#endif From 210a4d9e6ebcf7e2d8e8c7f0782c9b2e589c0f70 Mon Sep 17 00:00:00 2001 From: Sam Barnes Date: Fri, 23 Aug 2019 13:02:54 -0700 Subject: [PATCH 9/9] Cleaned up merge conflict gore. --- Gpufit/examples/Simple_Example.cpp | 85 ------------------------------ Gpufit/models/patlak.cuh | 8 --- 2 files changed, 93 deletions(-) diff --git a/Gpufit/examples/Simple_Example.cpp b/Gpufit/examples/Simple_Example.cpp index 5408a0b3..c91e0121 100644 --- a/Gpufit/examples/Simple_Example.cpp +++ b/Gpufit/examples/Simple_Example.cpp @@ -1,10 +1,6 @@ #include "../gpufit.h" #include -<<<<<<< HEAD #include -======= -#include ->>>>>>> 221fb4c80d94d67329e96357f7171552fba930ea void simple_example() { @@ -24,7 +20,6 @@ void simple_example() adjusted in the fit. */ -<<<<<<< HEAD /*************** definition of input and output parameters ***************/ // number of fits, number of points per fit @@ -34,92 +29,12 @@ void simple_example() // model ID and number of model parameters int const model_id = GAUSS_1D; size_t const n_model_parameters = 4; -======= - // start timer - clock_t time_start, time_end; - time_start = clock(); - - - // number of fits, fit points and parameters - size_t const n_fits = 10000; - size_t const n_points_per_fit = 60; - size_t const n_model_parameters = 2; - REAL snr = 0.8; - - // custom x positions for the data points of every fit, stored in user info - // time independent variable, given in minutes - REAL timeX[] = { 0.25, 0.5, 0.75, 1, 1.25, 1.5, 1.75, 2, 2.25, 2.5, 2.75, 3, 3.25, 3.5, 3.75, 4, 4.25, 4.5, 4.75, 5, - 5.25, 5.5, 5.75, 6, 6.25, 6.5, 6.75, 7, 7.25, 7.5, 7.75, 8, 8.25, 8.5, 8.75, 9, 9.25, 9.5, 9.75, 10, - 10.25, 10.5, 10.75, 11, 11.25, 11.5, 11.75, 12, 12.25, 12.5, 12.75, 13, 13.25, 13.5, 13.75, 14, 14.25, 14.5, 14.75, 15 }; - - // Concentration of plasma (independent) - REAL Cp[] = { 0.0f, 0.0f, 0.0f, 3.01846399851715f, 2.59801604007558f, 2.2361331285733f, 1.92465762011135f, 1.65656816551711f, 1.4258214335524f, - 1.22721588081636f, 1.05627449741415f, 0.909143885218726f, 0.782507393725825f, 0.6735103553914f, 0.579695735090254f, - 0.498948743091769f, 0.429449163006342f, 0.369630320068624f, 0.318143764811612f, 0.273828876023252f, 0.235686697768721f, - 0.20285742070682f, 0.174601000079374f, 0.150280473460109f, 0.12934760220805f, 0.111330512951924f, 0.0958230605172143f, - 0.0824756725126274f, 0.0709874691926393f, 0.0610994809603327f, 0.0525888106179893f, 0.0452636087696102f, 0.0389587491097867f, - 0.033532106110336f, 0.0288613511954976f, 0.0248411951843697f, 0.0213810148391187f, 0.018402810016092f, 0.0158394453694853f, - 0.013633136971665f, 0.0117341497352074f, 0.010099676273659f, 0.00869287192804919f, 0.00748202420651342f, 0.00643983791435146f, - 0.00554281985976681f, 0.00477074926518851f, 0.00410622194607174f, 0.00353425798195557f, 0.00304196403581309f, 0.00261824270962248f, - 0.00225354238438883f, 0.00193964190545541f, 0.00166946525943377f, 0.00143692206515917f, 0.00123677028298367f, 0.00106449804756952f, - 0.000916221960431984f, 0.000788599549519612f, 0.000678753922476738f }; - - - std::vector< REAL > user_info(2 * n_points_per_fit); - for (size_t i = 0; i < n_points_per_fit; i++) - { - user_info[i] = static_cast(timeX[i]); - } - - for (size_t i = n_points_per_fit; i < 2 * n_points_per_fit; i++) - { - user_info[i] = static_cast(Cp[i - n_points_per_fit]); - } - - // size of user info in bytes - size_t const user_info_size = 2 * n_points_per_fit * sizeof(REAL); - - // initialize random number generator - std::mt19937 rng; - rng.seed(time(NULL)); - std::uniform_real_distribution< REAL > uniform_dist(0, 1); - std::normal_distribution< REAL > normal_dist(0, 1); ->>>>>>> 221fb4c80d94d67329e96357f7171552fba930ea // initial parameters std::vector< REAL > initial_parameters(n_fits * n_model_parameters); // data std::vector< REAL > data(n_points_per_fit * n_fits); -<<<<<<< HEAD -======= - REAL mean_y = 0; - for (size_t i = 0; i != data.size(); i++) - { - size_t j = i / n_points_per_fit; // the fit - size_t k = i % n_points_per_fit; // the position within a fit - REAL x = 0; - for (int n = 1; n < k; n++) { - - REAL spacing = timeX[n] - timeX[n - 1]; - x += (Cp[n - 1] + Cp[n]) / 2 * spacing; - } - REAL y = true_parameters[0] * x + true_parameters[1] * Cp[k]; - //data[i] = y + normal_dist(rng); - //data[i] = y * (0.2f + 1.6f * uniform_dist(rng)); - data[i] = y; - mean_y += y; - //std::cout << data[i] << std::endl; - } - mean_y = mean_y / data.size(); - std::normal_distribution norm_snr(0,mean_y/snr); - for (size_t i = 0; i != data.size(); i++) - { - data[i] = data[i] + norm_snr(rng); - } - - ->>>>>>> 221fb4c80d94d67329e96357f7171552fba930ea // tolerance REAL const tolerance = 0.001f; diff --git a/Gpufit/models/patlak.cuh b/Gpufit/models/patlak.cuh index c205c0c6..94cb2da9 100644 --- a/Gpufit/models/patlak.cuh +++ b/Gpufit/models/patlak.cuh @@ -30,15 +30,7 @@ __device__ void calculate_patlak ( // function name } value[point_index] = parameters[0] * convCp + parameters[1] * Cp[point_index]; // formula calculating fit model values -<<<<<<< HEAD -<<<<<<< HEAD // C(t) = Ktrans * trapz(Cp(k)) + vp * Cp(k) -======= - // C(t) = Ktrans * trapz(Cp(k)) + vp * Cp(k) ->>>>>>> 221fb4c80d94d67329e96357f7171552fba930ea -======= - // C(t) = Ktrans * trapz(Cp(k)) + vp * Cp(k) ->>>>>>> 221fb4c80d94d67329e96357f7171552fba930ea /////////////////////////// derivative /////////////////////////// REAL * current_derivative = derivative + point_index;