/*
    * Licensed to the Apache Software Foundation (ASF) under one or more
    * contributor license agreements.  See the NOTICE file distributed with
    * this work for additional information regarding copyright ownership.
    * The ASF licenses this file to You under the Apache License, Version 2.0
    * (the "License"); you may not use this file except in compliance with
    * the License.  You may obtain a copy of the License at
    *
    *      http://www.apache.org/licenses/LICENSE-2.0
   *
   * Unless required by applicable law or agreed to in writing, software
   * distributed under the License is distributed on an "AS IS" BASIS,
   * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
   * See the License for the specific language governing permissions and
   * limitations under the License.
   */
  package org.apache.commons.math4.legacy.fitting.leastsquares;
  
  import org.apache.commons.math4.legacy.analysis.MultivariateMatrixFunction;
  import org.apache.commons.math4.legacy.analysis.MultivariateVectorFunction;
  import org.apache.commons.math4.core.jdkmath.JdkMath;
  
  import java.util.ArrayList;
  import java.util.List;
  
  class BevingtonProblem {
      private List<Double> time;
      private List<Double> count;
  
      BevingtonProblem() {
          time = new ArrayList<>();
          count = new ArrayList<>();
      }
  
      public void addPoint(double t, double c) {
          time.add(t);
          count.add(c);
      }
  
      public MultivariateVectorFunction getModelFunction() {
          return new MultivariateVectorFunction() {
              @Override
              public double[] value(double[] params) {
                  double[] values = new double[time.size()];
                  for (int i = 0; i < values.length; ++i) {
                      final double t = time.get(i);
                      values[i] = params[0] +
                          params[1] * JdkMath.exp(-t / params[3]) +
                          params[2] * JdkMath.exp(-t / params[4]);
                  }
                  return values;
              }
          };
      }
  
      public MultivariateMatrixFunction getModelFunctionJacobian() {
          return new MultivariateMatrixFunction() {
              @Override
              public double[][] value(double[] params) {
                  double[][] jacobian = new double[time.size()][5];
  
                  for (int i = 0; i < jacobian.length; ++i) {
                      final double t = time.get(i);
                      jacobian[i][0] = 1;
  
                      final double p3 =  params[3];
                      final double p4 =  params[4];
                      final double tOp3 = t / p3;
                      final double tOp4 = t / p4;
                      jacobian[i][1] = JdkMath.exp(-tOp3);
                      jacobian[i][2] = JdkMath.exp(-tOp4);
                      jacobian[i][3] = params[1] * JdkMath.exp(-tOp3) * tOp3 / p3;
                      jacobian[i][4] = params[2] * JdkMath.exp(-tOp4) * tOp4 / p4;
                  }
                  return jacobian;
              }
          };
      }
  }