Fix 123 (#124)

* Fixed call to logistic kernel in exp_power

now exponential power kernel is found and logistic kernel works as intended

* correct power law CDF

* swapped shape and scale parameter in power law distribution call. fixes power law distribution in #123

* simplified PDF and removed - in exp.

* simplified log normal cdf using erf function.

* revert cdf for tests

* add negative in lognormal PDF

* revert cdf

* simplified cdf using complimentary error function

* fixed missing semicolon

* Adding/updated tests and changing gamma.cdf

added tests to test_deterministic

Fixing gamma kernel call

reverting log normal for gamma tests

removing line in gamma

fixing problem in deterministic test

fix problem in exponential power

adding tests for log normal

* updated approximation for inverse erf

* updated results for power law tests

* updated power law test results

* updated formula for exponential_power.icdf

updated inverse erf approximation in normal

* updated exponential power test

* fixing error in exponential power test

* updating gamma kernel icdf

* fixed error in gamma kernel

* Update test_deterministic.cpp

* Update radial_kernel.hpp

revert power-law call back to distance then shape.

Co-authored-by: MargaretLawrimore <margaret.lawrimore@gmail.com>

include/pops/exponential_power_kernel.hpp

@@ -69,7 +69,7 @@ class ExponentialPowerKernel
     double pdf(double x)
     {
         return (beta / (2 * alpha * std::tgamma(1.0 / beta)))
-               * pow(exp(-x / alpha), beta);
+               * exp(-pow(x / alpha, beta));
     }
 
     /*!
@@ -83,9 +83,13 @@ class ExponentialPowerKernel
         if (x <= 0 || x >= 1) {
             throw std::invalid_argument("icdf: x must be between 0.0 and 1.0");
         }
+        float sign = ((x - 0.5) > 0) ? 1.0 : ((x - 0.5) < 0 ? -1.0 : 0);
+        if (sign == 0) {
+            return 0;
+        }
         GammaKernel gamma_distribution(1.0 / beta, 1.0 / pow(alpha, beta));
         double gamma = gamma_distribution.icdf(2 * std::abs(x - 0.5));
-        return (x - 0.5) * pow(gamma, 1.0 / beta);
+        return sign * pow(gamma, 1.0 / beta);
     }
 };
 

include/pops/gamma_kernel.hpp

@@ -85,9 +85,9 @@ class GammaKernel
         double beta = 1.0 / theta;
         for (int i = 0; i < alpha; i++) {
             // tgamma = (i-1)! used since c++ has no factorial in std lib
-            sum += pow(beta * x, i) / std::tgamma(i + 1);
+            sum += (1.0 / std::tgamma(i + 1)) * exp(-beta * x) * pow(beta * x, i);
         }
-        return 1 - sum * exp(-beta * x);
+        return 1 - sum;
     }
 
     /*!
@@ -110,24 +110,27 @@ class GammaKernel
         LogNormalKernel lognormal(1);
         double guess = lognormal.icdf(x);
         double check = cdf(guess);
-        double numiterations = 500;  // will need to adjust this
+        double numiterations = 1000;  // will need to adjust this
         double precision = 0.001;  // will need to adjust this
         for (int i = 0; i < numiterations; i++) {
             if (check < (x - precision) || check > (x + precision)) {
                 double dif = check - x;
                 // if dif is positive guess is greater than needed
                 // if dif is negative guess is less than needed
                 double past_guess = guess;
-                double derivative = (check - x) / pdf(guess);
+                double derivative = dif / pdf(guess);
                 // limit size of next guess
                 guess = std::max(guess / 10, std::min(guess * 10, guess - derivative));
                 check = cdf(guess);
                 // Check if we went to far and need to backtrack
-                for (int j = 0; j < 10; j++) {
-                    if (std::abs(dif) < std::abs(check - x)) {
-                        past_guess = guess;
-                        guess = (guess + past_guess) / 2.0;
-                        check = cdf(guess);
+                int count = 0;
+                bool run = true;
+                while ((std::abs(dif) < std::abs(check - x)) && run) {
+                    guess = (guess + past_guess) / 2.0;
+                    check = cdf(guess);
+                    count++;
+                    if (count > 20) {
+                        run = false;
                     }
                 }
             }

include/pops/lognormal_kernel.hpp

@@ -72,7 +72,7 @@ class LogNormalKernel
         if (x == 0) {
             return 0;
         }
-        return (1 / x) * (1 / (sigma * sqrt(2 * M_PI)))
+        return (1 / (x * sigma * sqrt(2 * M_PI)))
                * exp(-(pow(log(x), 2)) / (2 * pow(sigma, 2)));
     }
 
@@ -93,9 +93,14 @@ class LogNormalKernel
         double y = (2 * x) - 1;
         float sign = (y < 0) ? -1.0f : 1.0f;
         //  0.147 used for a relative error of about 2*10^-3
-        float b = 2 / (M_PI * 0.147) + 0.5f * log(1 - pow(y, 2));
+        // float b = 2 / (M_PI * 0.147) + 0.5f * log(1 - pow(y, 2));
+        // double inverf =
+        //    (sign * sqrt(-b + sqrt(pow(b, 2) - (1 / (0.147) * log(1 - pow(y, 2))))));
+        double a = 0.140012;
+        double t = 2.0 / (M_PI * a);
+        double l = log(1 - pow(y, 2));
         double inverf =
-            (sign * sqrt(-b + sqrt(pow(b, 2) - (1 / (0.147) * log(1 - pow(y, 2))))));
+            sign * sqrt(sqrt(pow(t + (l / 2.0), 2) - (l / a)) - (t + (l / 2.0)));
         return exp(sqrt(2 * pow(sigma, 2)) * inverf);
     }
 };

include/pops/normal_kernel.hpp

@@ -87,11 +87,11 @@ class NormalKernel
         //  approximation for inverse error function
         double y = (2 * x) - 1;
         float sign = (y < 0) ? -1.0f : 1.0f;
-        //  0.147 used for a relative error of about 2*10^-3
-        float b = 2.0 / (M_PI * 0.147) + 0.5f * log(1 - pow(y, 2));
+        double a = 0.140012;
+        double t = 2.0 / (M_PI * a);
+        double l = log(1 - pow(y, 2));
         double inverf =
-            (sign * sqrt(-b + sqrt(pow(b, 2) - (1.0 / (0.147) * log(1 - pow(y, 2))))));
-
+            sign * sqrt(sqrt(pow(t + (l / 2.0), 2) - (l / a)) - (t + (l / 2.0)));
         return sigma * std::sqrt(2) * inverf;
     }
 };

include/pops/power_law_kernel.hpp

@@ -31,12 +31,12 @@ using std::pow;
 class PowerLawKernel
 {
 protected:
-    double xmin;
     double alpha;
+    double xmin;
     std::uniform_real_distribution<double> distribution;
 
 public:
-    PowerLawKernel(double xm, double a) : xmin(xm), alpha(a), distribution(0.0, 1.0)
+    PowerLawKernel(double a, double xm) : alpha(a), xmin(xm), distribution(0.0, 1.0)
     {
         //         if (alpha < 1.0) {
         //             throw std::invalid_argument("alpha must be greater than or equal
@@ -92,7 +92,7 @@ class PowerLawKernel
         if (x <= 0 || x >= 1) {
             throw std::invalid_argument("icdf: x must be between 0 and 1.0");
         }
-        return xmin * pow(x, (1.0 / (-alpha + 1.0)));
+        return pow(x / xmin, (-alpha + 1.0));
     }
 };
 

include/pops/radial_kernel.hpp

@@ -184,7 +184,7 @@ class RadialDispersalKernel
         else if (dispersal_kernel_type_ == DispersalKernelType::Gamma) {
             distance = std::abs(gamma_distribution.random(generator));
         }
-        else if (dispersal_kernel_type_ == DispersalKernelType::Logistic) {
+        else if (dispersal_kernel_type_ == DispersalKernelType::ExponentialPower) {
             distance = std::abs(exponential_power_distribution.random(generator));
         }
         else if (dispersal_kernel_type_ == DispersalKernelType::Logistic) {