text:
# contains the code for the PCMCM, MCWM and A(DA)-GP-MCMC algorithms with assiciated help functions

################################################################################
######               algorithms                                            #####
################################################################################

# PMCMC/MCWM

doc"""
    mcmc(problem::Problem, store_data::Bool=false, return_cov_matrix::Bool=false)

Runs PMCMC or MCWM with adaptive gaussian random walk.
"""
function mcmc(problem::Problem, store_data::Bool=false, return_cov_matrix::Bool=false)

  # data
  Z = problem.data.Z

  # algorithm parameters
  R = problem.alg_param.R
  N = problem.alg_param.N
  burn_in = problem.alg_param.burn_in
  alg = problem.alg_param.alg
  pf_alg = problem.alg_param.pf_alg
  nbr_of_cores = problem.alg_param.nbr_of_cores
  nbr_x0  = problem.alg_param.nbr_x0
  nbr_x = problem.alg_param.nbr_x
  subsample_interval = problem.alg_param.subsample_int
  dt = problem.alg_param.dt
  dt_U = problem.alg_param.dt_U

  # model parameters
  theta_true = problem.model_param.theta_true
  theta_known = problem.model_param.theta_known
  theta_0 = problem.model_param.theta_0

  # pre-allocate matricies and vectors
  Theta = zeros(length(theta_0),R)
  loglik = zeros(R)
  accept_vec = zeros(R)
  prior_vec = zeros(R)
  theta_star = zeros(length(theta_0),1)
  loglik_current = zero(Float64)
  loglik_star = zero(Float64)
  a_log = zero(Float64)

  if store_data
    Theta_val = zeros(length(theta_0),R-burn_in)
    loglik_val = zeros(R-burn_in)
  end

  # parameters for adaptive update
  adaptive_update_params = set_adaptive_alg_params(problem.adaptive_update, length(theta_0),Theta[:,1], R)

  # parameters for prior dist
  dist_type = problem.prior_dist.dist
  prior_parameters = problem.prior_dist.prior_parameters

  @printf "#####################################################################\n"

  # print information at start of algorithm
  @printf "Starting MCMC with adaptive RW estimating %d parameters\n" length(theta_true)
  @printf "Algorithm: %s\n" alg
  @printf "Adaptation algorithm: %s\n" typeof(problem.adaptive_update)
  @printf "Prior distribution: %s\n" problem.prior_dist.dist
  @printf "Particel filter: %s, on %d cores\n" pf_alg nbr_of_cores
  @printf "Nbr particles for particel filter: %d\n" N

  nbr_of_proc = set_nbr_cores(nbr_of_cores, pf_alg)
  #nbr_of_proc = nbr_of_cores
  loglik_vec = SharedArray{Float64}(nbr_of_proc)

  # print acceptance rate each print_interval:th iteration
  print_interval = 1000

  # first iteration
  @printf "Iteration: %d\n" 1 # print first iteration
  @printf "Covariance:\n"
  print_covariance(problem.adaptive_update,adaptive_update_params, 1)

  Theta[:,1] = theta_0
  loglik[1] = pf_paralell(Z, Theta[:,1],theta_known,N,dt,dt_U,nbr_x0, nbr_x,subsample_interval,false,false, nbr_of_proc,loglik_vec)

  # print start loglik
  @printf "Loglik: %.4f \n" loglik[1]

  for r = 2:R

    # set print_on to false, only print each print_interval:th iteration
    print_on = false

    # Gaussian random walk
    (theta_star,) = gaussian_random_walk(problem.adaptive_update, adaptive_update_params, Theta[:,r-1], r)

    # print acceptance rate for the last print_interval iterations
    if mod(r-1,print_interval) == 0
      # print percentage done
      @printf "Percentage done: %.2f %% \n" 100*(r-1)/R
      print_on = true # print ESS and Nbr resample each print_interval:th iteration
      # print accaptance rate
      @printf "Acceptance rate on iteration %d to %d is %.4f\n" r-print_interval r-1  sum(accept_vec[r-print_interval:r-1])/( r-1 - (r-print_interval) )
      # print covariance matrix
      @printf "Covariance:\n"
      print_covariance(problem.adaptive_update,adaptive_update_params, r)
      # print loglik
      @printf "Loglik: %.4f \n" loglik[r-1]
      # print log-lik vector
      @printf "Loglik values on different cores:\n"
      println(loglik_vec)
    end

    # calc loglik using proposed parameters
    loglik_star = pf_paralell(Z, theta_star,theta_known, N,dt,dt_U,nbr_x0, nbr_x,subsample_interval,print_on,false, nbr_of_proc,loglik_vec)

    if store_data && r > burn_in # store data
      Theta_val[:,r-burn_in] = theta_star
      loglik_val[r-burn_in] = loglik_star
    end

    # run MCWM or PMCMC
    if alg == "MCWM"
      loglik_current =  pf_paralell(Z, Theta[:,r-1],theta_known, N,dt,dt_U,nbr_x0, nbr_x,subsample_interval,print_on, false, nbr_of_proc,loglik_vec)
    else
      loglik_current = loglik[r-1]
    end

    prior_log_star = evaluate_prior(theta_star,prior_parameters, dist_type)
    prior_log_old = evaluate_prior(Theta[:,r-1],prior_parameters, dist_type)

    jacobian_log_star = jacobian(theta_star)
    jacobian_log_old = jacobian(Theta[:,r-1])

    a_log = loglik_star + prior_log_star +  jacobian_log_star - (loglik_current +  prior_log_old + jacobian_log_old)
    # generate log(u)

    u_log = log(rand())
    accept = u_log < a_log # calc accaptace decision

    # update chain
    if accept # the proposal is accapted
      Theta[:,r] = theta_star # update chain with new values
      loglik[r] = loglik_star
      accept_vec[r] = 1
    else
      Theta[:,r] = Theta[:,r-1] # keep old values
      loglik[r] = loglik[r-1]
    end

    # Adaptation
    adaptation(problem.adaptive_update, adaptive_update_params, Theta, r,a_log)

  end

  @printf "Ending MCMC with adaptive RW estimating %d parameters\n" length(theta_true)
  @printf "Algorithm: %s\n" alg
  @printf "Adaptation algorithm: %s\n" typeof(problem.adaptive_update)
  @printf "Prior distribution: %s\n" problem.prior_dist.dist
  @printf "Particel filter: %s, on %d cores\n" pf_alg nbr_of_cores

  @printf "#####################################################################\n"

  if store_data && return_cov_matrix
    cov_prop_kernel = get_covariance(problem.adaptive_update,adaptive_update_params, R)
    return return_results(Theta,loglik,accept_vec,prior_vec, problem,adaptive_update_params), Theta_val, loglik_val, cov_prop_kernel
  elseif store_data
   return return_results(Theta,loglik,accept_vec,prior_vec, problem,adaptive_update_params), Theta_val, loglik_val
  else
    return_results(Theta,loglik,accept_vec,prior_vec, problem,adaptive_update_params)
  end

end

# DA-GP-MCMC

doc"""
    dagpmcmc(problem_traning::Problem, problem::gpProblem, gp::GPModel, cov_matrix::Matrix)

Runs the DA-GP-MCMC algorithm.
"""
function dagpmcmc(problem_traning::Problem, problem::gpProblem, gp::GPModel, cov_matrix::Matrix,return_run_info::Bool=false)

  # data
  Z = problem.data.Z

  # algorithm parameters
  R = problem.alg_param.R # number of iterations
  N = problem.alg_param.N # number of particels
  burn_in = problem.alg_param.burn_in # burn in
  length_training_data = problem.alg_param.length_training_data # length of training data
  independet_sampling = problem.alg_param.independet_sampling # length of training data
  compare_GP_and_PF = problem.alg_param.compare_GP_and_PF # compare GP and PF if true
  noisy_est = problem.alg_param.noisy_est # use noisy estimation
  pred_method = problem.alg_param.pred_method # method used for predictions
  est_method = problem.alg_param.est_method # methods used to estimate the parameters of the GP model
  pf_alg = problem.alg_param.pf_alg # pf algorithm
  alg = problem.alg_param.alg # use PMCMC or MCWM
  #print_interval = problem.alg_param.print_interval # print the accaptance rate every print_interval:th iteration
  lasso = problem.alg_param.lasso # use Lasso
  beta_MH = problem.alg_param.beta_MH

  nbr_of_cores = problem.alg_param.nbr_of_cores
  nbr_x0  = problem.alg_param.nbr_x0
  nbr_x = problem.alg_param.nbr_x
  subsample_interval = problem.alg_param.subsample_int
  dt = problem.alg_param.dt

  # model parameters
  theta_true = problem.model_param.theta_true
  theta_known = problem.model_param.theta_known
  theta_0 = problem.model_param.theta_0

  # pre-allocate matricies and vectors
  Theta = zeros(length(theta_0), R)
  loglik = zeros(R)
  accept_vec = zeros(R)
  prior_vec = zeros(R)
  theta_star = zeros(length(theta_0))
  #compare_GP_PF = zeros(2,R-length_training_data-burn_in)
  compare_GP_PF = zeros(2,R)
  #data_gp_pf = zeros(length(theta_0)+2,R-length_training_data-burn_in)
  data_gp_pf = zeros(length(theta_0)+2,R)
  data_training = zeros(1+length(theta_0), length_training_data)
  #accept_prob_log = zeros(2, R-length_training_data-burn_in) # [gp ; pf]
  accept_prob_log = zeros(2, R) # [gp ; pf]

  loglik_star = zero(Float64)
  loglik_gp = zeros(Float64)
  loglik_gp_old = zero(Float64)
  loglik_gp_new = zero(Float64)
  index_keep_gp_er = zero(Int64)
  nbr_early_rejections = zero(Int64)
  accept = true
  MH_direct = false
  nbr_ordinary_mh = 0
  nbr_ordinary_mh_accapte = 0
  nbr_split_accaptance_region = 0
  nbr_split_accaptance_region_early_accept = 0
  nbr_second_stage_accepted = 0
  nbr_second_stage = 0
  nbr_run_DA = 0
  assumption_list = []
  loglik_list = []
  a_log = zero(Float64)
  loglik_current = zero(Float64)
  nbr_eval_pf = 0
  nbr_eval_pf_secound_stage = 0

  # starting values for times:
  time_pre_er = zero(Float64)
  time_fit_gp = zero(Float64)
  time_er_part = zero(Float64)

  # parameters for prior dist
  dist_type = problem.prior_dist.dist
  prior_parameters = problem.prior_dist.prior_parameters

  # prop kernl for DA-GP-MCMC
  xi = 1.2
  problem.adaptive_update = noAdaptation(xi^2*cov_matrix)

  adaptive_update_params = set_adaptive_alg_params(problem.adaptive_update, length(theta_0),Theta[:,1], R)

  # prop kernl for MH_direct
  kernel_MH_direct = noAdaptation(cov_matrix)
  adaptive_update_params_MH_direct = set_adaptive_alg_params(kernel_MH_direct, length(theta_0),Theta[:,1], R)

  @printf "#####################################################################\n"

  # print information at start of algorithm
  @printf "Starting DA-GP-MCMC estimating %d parameters \n" length(theta_true)
  @printf "MCMC algorithm: %s\n" alg
  @printf "Adaptation algorithm: %s\n" typeof(problem.adaptive_update)
  @printf "Prior distribution: %s\n" problem.prior_dist.dist
  @printf "Particle filter: %s, on %d cores\n" pf_alg nbr_of_cores
  @printf "Nbr particles for particle filter: %d\n" N

  @printf "Covariance - kernel_MH_direct:\n"
  print_covariance(kernel_MH_direct,adaptive_update_params_MH_direct, 1)

  # first iteration
  @printf "Iteration: %d\n" 1 # print first iteration
  @printf "Covariance:\n"
  print_covariance(problem.adaptive_update,adaptive_update_params, 1)

  # set nbr of cores to use for parallel pf
  nbr_of_proc = set_nbr_cores(nbr_of_cores, pf_alg)
  #nbr_of_proc = nbr_of_cores
  loglik_vec = SharedArray{Float64}(nbr_of_cores)

  # print acceptance rate each print_interval:th iteration
  print_interval = 1000

  # first iteration
  @printf "Iteration: %d\n" 1 # print first iteration
  @printf "Covariance:\n"
  print_covariance(problem.adaptive_update,adaptive_update_params, 1)

  # first iteration
  Theta[:,1] = theta_0
  loglik[1] = pf_paralell(Z, Theta[:,1],theta_known,N,dt,dt_U,nbr_x0, nbr_x,subsample_interval,false,false, nbr_of_proc, loglik_vec)

  if alg == "MCWM"
    # do nothing
  else
    loglik_gp_old = predict(Theta[:,1], gp, pred_method,est_method,noisy_est)[1]
  end

  # print start loglik
  @printf "Loglik start: %.4f \n" loglik[1]

  tic()
  for r = 2:R

    # set print_on to false, only print each print_interval:th iteration
    print_on = false

    # print acceptance rate for the last print_interval iterations
    if mod(r-1,print_interval) == 0
      # print percentage done
      @printf "Percentage done: %.2f %% \n" 100*(r-1)/R
      print_on = true # print ESS and Nbr resample each print_interval:th iteration
      # print accaptance rate
      @printf "Acceptance rate on iteration %d to %d is %.4f\n" r-print_interval r-1  sum(accept_vec[r-print_interval:r-1])/( r-1 - (r-print_interval) )
      # print covariance function
      @printf "Covariance:\n"
      print_covariance(problem.adaptive_update,adaptive_update_params, r)
      # print loglik
      @printf "Loglik: %.8f \n" loglik[r-1]
      # print log-lik vector
      @printf "Loglik values on different cores:\n"
      println(loglik_vec)
    end

    MH_direct = rand() < beta_MH # we always run the early-rejection scheme

    if MH_direct

      # secound stage direct

      nbr_ordinary_mh +=  1

      # Gaussian random walk using secound stage direct kernel
      (theta_star, ) = gaussian_random_walk(kernel_MH_direct, adaptive_update_params_MH_direct, Theta[:,r-1], r)

      # calc loglik using proposed parameters
      loglik_star = pf_paralell(Z, theta_star,theta_known, N,dt,dt_U,nbr_x0, nbr_x,subsample_interval,false,false, nbr_of_proc,loglik_vec)

      # run MCWM or PMCMC
      if alg == "MCWM"
        loglik_current =  pf_paralell(Z, Theta[:,r-1],theta_known, N,dt,dt_U,nbr_x0, nbr_x,subsample_interval,false, false, nbr_of_proc,loglik_vec)
      else
        loglik_current = loglik[r-1]
      end

      nbr_eval_pf += 1

      prior_log_star = evaluate_prior(theta_star,prior_parameters, dist_type)
      prior_log_old = evaluate_prior(Theta[:,r-1],prior_parameters, dist_type)

      jacobian_log_star = jacobian(theta_star)
      jacobian_log_old = jacobian(Theta[:,r-1])

      a_log = loglik_star + prior_log_star +  jacobian_log_star - (loglik_current +  prior_log_old + jacobian_log_old)

      accept = log(rand()) < a_log # calc accaptace decision

      if accept # the proposal is accapted
        nbr_ordinary_mh_accapte += 1
        Theta[:,r] = theta_star # update chain with new values
        loglik[r] = loglik_star
        accept_vec[r] = 1
      else
        Theta[:,r] = Theta[:,r-1] # keep old values
        loglik[r] = loglik[r-1]
      end

    else

      # first stage
      nbr_run_DA +=  1

      # set proposal
      (theta_star, ) = gaussian_random_walk(problem.adaptive_update, adaptive_update_params, Theta[:,r-1], r)

      (loglik_gp_pred,loglik_gp_new_std) = predict(theta_star, gp, pred_method,est_method,noisy_est,true)

      loglik_gp_new = loglik_gp_pred[1]
      loglik_gp_new_std = loglik_gp_new_std[1]

      prior_log_star = evaluate_prior(theta_star,prior_parameters,dist_type)
      prior_log = evaluate_prior(Theta[:,r-1],prior_parameters,dist_type)

      prior_log_star = evaluate_prior(theta_star,prior_parameters,dist_type)
      prior_log_old = evaluate_prior(Theta[:,r-1],prior_parameters,dist_type)

      jacobian_log_star = jacobian(theta_star)
      jacobian_log_old = jacobian(Theta[:,r-1])

      # should we recompute the loglik_gp_old value here?
      # we currently recompute loglik_gp_old here!

      if alg == "MCWM"
        loglik_gp_old = predict(Theta[:,r-1], gp, pred_method,est_method,noisy_est)[1]
      else
      end

      a_gp = loglik_gp_new + prior_log_star +  jacobian_log_star -  (loglik_gp_old - prior_log_old - jacobian_log_old)

      accept = log(rand()) < a_gp # calc accept

      if !accept
        # keep old values
        nbr_early_rejections = nbr_early_rejections + 1
        Theta[:,r] = Theta[:,r-1]
        loglik[r] = loglik[r-1]
      else

        # run pf

        nbr_second_stage += 1

        loglik_star = pf_paralell(Z, theta_star,theta_known,N,dt,dt_U,nbr_x0, nbr_x,subsample_interval,false,false, nbr_of_proc,loglik_vec)

        if alg == "MCWM"
          loglik_current =  pf_paralell(Z, Theta[:,r-1],theta_known, N,dt,dt_U,nbr_x0, nbr_x,subsample_interval,false, false, nbr_of_proc,loglik_vec)
        else
          loglik_current = loglik[r-1]
        end

        nbr_eval_pf += 1
        nbr_eval_pf_secound_stage += 1

        a_log = (loglik_star + loglik_gp_old)  -  (loglik_current + loglik_gp_new)

        accept = log(rand()) < a_log # calc accaptance decision

        if accept # the proposal is accapted
          nbr_second_stage_accepted = nbr_second_stage_accepted+1
          Theta[:,r] = theta_star # update chain with proposal
          loglik[r] = loglik_star
          accept_vec[r] = 1
          if alg == "MCWM"
            # do nothing
          else
            loglik_gp_old = loglik_gp_new
          end
        else
          Theta[:,r] = Theta[:,r-1] # keep old values
          loglik[r] = loglik[r-1]
        end
      end
    end
  end

  time_da_part = toc()
  times = [time_pre_er time_fit_gp time_da_part]

  @printf "Ending DA-GP-MCMC with adaptive RW estimating %d parameters\n" length(theta_true)
  @printf "Algorithm: %s\n" alg
  @printf "Adaptation algorithm: %s\n" typeof(problem.adaptive_update)
  @printf "Prior distribution: %s\n" problem.prior_dist.dist
  @printf "Particle filter: %s\n" pf_alg
  @printf "Time pre-er:  %.0f\n" time_pre_er
  @printf "Time fit GP model:  %.0f\n" time_fit_gp
  @printf "Time er-part:  %.0f\n" time_da_part

  @printf "Number early-rejections: %.d\n"  nbr_early_rejections

  @printf "Number of cases directly to ordinary MH: %d\n"  nbr_ordinary_mh
  @printf "Number of cases directly to ordinary MH accepted: %d\n"  nbr_ordinary_mh_accapte

  @printf "Number cases in in DA part: %d\n"  nbr_run_DA
  @printf "Number cases in second stage: %d\n"  nbr_second_stage
  @printf "Number accepted in second stage: %d\n"  nbr_second_stage_accepted

  @printf "Total number of evaluations of the particle filter: %d\n" nbr_eval_pf

  @printf "Acceptance rate for ordinary MH accepted: %.4f\n"  nbr_ordinary_mh_accapte/nbr_ordinary_mh*100
  @printf "Acceptance rate a_1: %.4f\n"  nbr_second_stage/nbr_run_DA*100
  @printf "Acceptance rate a_2: %.4f\n"  nbr_second_stage_accepted/nbr_second_stage*100
  @printf "Acceptance rate a_1*a_2: %.4f\n"  (nbr_second_stage/nbr_run_DA)*(nbr_second_stage_accepted/nbr_second_stage)*100
  @printf "Acceptance rate a (for entier algorithm): %.4f\n" sum(accept_vec)/R*100

  @printf "#####################################################################\n"

  # return resutls
  if return_run_info
    run_info = [nbr_eval_pf;
                nbr_eval_pf_secound_stage;
                nbr_second_stage;
                nbr_ordinary_mh]
    return return_gp_results(gp, Theta,loglik,accept_vec,prior_vec, compare_GP_PF, data_gp_pf,nbr_early_rejections, problem, adaptive_update_params,accept_prob_log,times), run_info
  else
    return return_gp_results(gp, Theta,loglik,accept_vec,prior_vec, compare_GP_PF, data_gp_pf,nbr_early_rejections, problem, adaptive_update_params,accept_prob_log,times)
  end

end

# ADA-GP_MCMC

doc"""
    problem_traning::Problem, problem::gpProblem, gp::GPModel, casemodel::CaseModel, cov_matrix::Matrix)

Runs the ADA-GP-MCMC algorithm.
"""
function adagpmcmc(problem_traning::Problem, problem::gpProblem, gp::GPModel, casemodel::CaseModel, cov_matrix::Matrix, return_run_info::Bool=false)

  # data
  Z = problem.data.Z

  # algorithm parameters
  R = problem.alg_param.R # number of iterations
  N = problem.alg_param.N # number of particels
  burn_in = problem.alg_param.burn_in # burn in
  length_training_data = problem.alg_param.length_training_data # length of training data
  independet_sampling = problem.alg_param.independet_sampling # length of training data
  compare_GP_and_PF = problem.alg_param.compare_GP_and_PF # compare GP and PF if true
  noisy_est = problem.alg_param.noisy_est # use noisy estimation
  pred_method = problem.alg_param.pred_method # method used for predictions
  est_method = problem.alg_param.est_method # methods used to estimate the parameters of the GP model
  pf_alg = problem.alg_param.pf_alg # pf algorithm
  alg = problem.alg_param.alg # use PMCMC or MCWM
  #print_interval = problem.alg_param.print_interval # print the accaptance rate every print_interval:th iteration
  lasso = problem.alg_param.lasso # use Lasso
  beta_MH = problem.alg_param.beta_MH

  nbr_of_cores = problem.alg_param.nbr_of_cores
  nbr_x0  = problem.alg_param.nbr_x0
  nbr_x = problem.alg_param.nbr_x
  subsample_interval = problem.alg_param.subsample_int
  dt = problem.alg_param.dt

  # model parameters
  theta_true = problem.model_param.theta_true
  theta_known = problem.model_param.theta_known
  theta_0 = problem.model_param.theta_0

  # pre-allocate matricies and vectors
  Theta = zeros(length(theta_0), R)
  loglik = zeros(R)
  accept_vec = zeros(R)
  prior_vec = zeros(R)
  theta_star = zeros(length(theta_0))
  #compare_GP_PF = zeros(2,R-length_training_data-burn_in)
  compare_GP_PF = zeros(2,R)

  #data_gp_pf = zeros(length(theta_0)+2,R-length_training_data-burn_in)
  data_gp_pf = zeros(length(theta_0)+2,R)
  data_training = zeros(1+length(theta_0), length_training_data)
  #accept_prob_log = zeros(2, R-length_training_data-burn_in) # [gp ; pf]
  accept_prob_log = zeros(2, R) # [gp ; pf]

  loglik_star = zero(Float64)
  loglik_gp = zeros(Float64)
  loglik_gp_old = zero(Float64)
  loglik_gp_new = zero(Float64)
  index_keep_gp_er = zero(Int64)
  nbr_early_rejections = zero(Int64)
  accept = true
  MH_direct = false
  nbr_ordinary_mh = 0
  nbr_ordinary_mh_accapte = 0
  nbr_run_DA = 0
  nbr_split_accaptance_region_early_accept = 0
  nbr_split_accaptance_region_early_reject = 0
  nbr_second_stage_accepted = 0
  nbr_second_stage = 0
  loglik_list = []
  a_log = zero(Float64)

  # starting values for times:
  time_pre_er = zero(Float64)
  time_fit_gp = zero(Float64)
  time_er_part = zero(Float64)

  nbr_case_1 = zero(Int64)
  nbr_case_2 = zero(Int64)
  nbr_case_3 = zero(Int64)
  nbr_case_4 = zero(Int64)

  nbr_eval_pf = zero(Int64)
  nbr_case_13 = zero(Int64)
  nbr_case_24 = zero(Int64)

  nbr_case_pf_1 = zero(Int64)
  nbr_case_pf_2 = zero(Int64)
  nbr_case_pf_3 = zero(Int64)
  nbr_case_pf_4 = zero(Int64)

  # da new
  loglik_gp_new_std = 0

  #(~,std_loglik_training) = predict(theta_training, gp, pred_method,est_method,noisy_est,true)
  loglik_gp_new_std = 0

  # parameters for prior dist
  dist_type = problem.prior_dist.dist
  prior_parameters = problem.prior_dist.prior_parameters

  # prop kernl for ADA-GP-MCMC
  xi = 1.2
  problem.adaptive_update = noAdaptation(xi^2*cov_matrix)

  adaptive_update_params = set_adaptive_alg_params(problem.adaptive_update, length(theta_0),Theta[:,1], R)

  # prop kernl for MH_direct
  kernel_MH_direct = noAdaptation(cov_matrix)
  adaptive_update_params_MH_direct = set_adaptive_alg_params(kernel_MH_direct, length(theta_0),Theta[:,1], R)

  @printf "#####################################################################\n"

  # print information at start of algorithm
  @printf "Starting ADA-GP-MCMC with adaptive RW estimating %d parameters\n" length(theta_true)
  @printf "MCMC algorithm: %s\n" alg
  @printf "Adaptation algorithm: %s\n" typeof(problem.adaptive_update)
  @printf "Select case model: %s\n" typeof(casemodel)
  @printf "Prior distribution: %s\n" problem.prior_dist.dist
  @printf "Particle filter: %s, on %d cores\n" pf_alg nbr_of_cores
  @printf "Nbr particles for particle filter: %d\n" N

  @printf "Covariance - kernel_MH_direct:\n"
  print_covariance(kernel_MH_direct,adaptive_update_params_MH_direct, 1)

  # first iteration
  @printf "Iteration: %d\n" 1 # print first iteration
  @printf "Covariance:\n"
  print_covariance(problem.adaptive_update,adaptive_update_params, 1)

  # set nbr of cores to use for parallel pf
  nbr_of_proc = set_nbr_cores(nbr_of_cores, pf_alg)
  #nbr_of_proc = nbr_of_cores
  loglik_vec = SharedArray{Float64}(nbr_of_proc)

  # print acceptance rate each print_interval:th iteration
  print_interval = 1000

  # first iteration
  @printf "Iteration: %d\n" 1 # print first iteration
  @printf "Covariance:\n"
  print_covariance(problem.adaptive_update,adaptive_update_params, 1)

  Theta[:,1] = theta_0
  loglik[1] = pf_paralell(Z, Theta[:,1],theta_known,N,dt,dt_U,nbr_x0, nbr_x,subsample_interval,false,false, nbr_of_proc,loglik_vec)

  # print start loglik
  @printf "Loglik: %.4f \n" loglik[1]

  tic()
  for r = 2:R

    # set print_on to false, only print each print_interval:th iteration
    print_on = false

    # print acceptance rate for the last print_interval iterations
    if mod(r-1,print_interval) == 0
      # print percentage done
      @printf "Percentage done: %.2f %% \n" 100*(r-1)/R
      print_on = true # print ESS and Nbr resample each print_interval:th iteration
      # print accaptance rate
      @printf "Acceptance rate on iteration %d to %d is %.4f\n" r-print_interval r-1  sum(accept_vec[r-print_interval:r-1])/( r-1 - (r-print_interval) )
      # print covariance function
      @printf "Covariance:\n"
      print_covariance(problem.adaptive_update,adaptive_update_params, r)
      # print loglik
      @printf "Loglik: %.4f \n" loglik[r-1]
      @printf "Loglik values on different cores:\n"
      println(loglik_vec)
    end

    MH_direct = rand() < beta_MH # we always run the early-rejection scheme

    if MH_direct

      # secound stage direct
      nbr_ordinary_mh += 1

      # Gaussian random walk using secound stage direct kernel
      (theta_star, ) = gaussian_random_walk(kernel_MH_direct, adaptive_update_params_MH_direct, Theta[:,r-1], r)

      # calc loglik using proposed parameters
      loglik_star = pf_paralell(Z, theta_star,theta_known, N,dt,dt_U,nbr_x0, nbr_x,subsample_interval,false,false, nbr_of_proc,loglik_vec)

      # run MCWM or PMCMC
      if alg == "MCWM"
        loglik_current =  pf_paralell(Z, Theta[:,r-1],theta_known, N,dt,dt_U,nbr_x0, nbr_x,subsample_interval,false, false, nbr_of_proc,loglik_vec)
      else
        loglik_current = loglik[r-1]
      end

      nbr_eval_pf += 1

      prior_log_star = evaluate_prior(theta_star,prior_parameters, dist_type)
      prior_log_old = evaluate_prior(Theta[:,r-1],prior_parameters, dist_type)

      jacobian_log_star = jacobian(theta_star)
      jacobian_log_old = jacobian(Theta[:,r-1])

      a_log = loglik_star + prior_log_star +  jacobian_log_star - (loglik_current +  prior_log_old + jacobian_log_old)

      accept = log(rand()) < a_log # calc accaptace decision

      if accept # the proposal is accapted
        nbr_ordinary_mh_accapte += 1
        Theta[:,r] = theta_star # update chain with new values
        loglik[r] = loglik_star
        accept_vec[r] = 1
      else
        Theta[:,r] = Theta[:,r-1] # keep old values
        loglik[r] = loglik[r-1]
      end

    else

      # stage 1
      nbr_run_DA += 1
      # set proposal
      (theta_star, ) = gaussian_random_walk(problem.adaptive_update, adaptive_update_params, Theta[:,r-1], r)

      (loglik_gp_pred,loglik_gp_new_std) = predict(theta_star, gp, pred_method,est_method,noisy_est,true)

      loglik_gp_new = loglik_gp_pred[1]
      loglik_gp_new_std = loglik_gp_new_std[1]

      prior_log_star = evaluate_prior(theta_star,prior_parameters, problem.prior_dist.dist)
      prior_log_old = evaluate_prior(Theta[:,r-1],prior_parameters, problem.prior_dist.dist)

      jacobian_log_star = jacobian(theta_star)
      jacobian_log_old = jacobian(Theta[:,r-1])

      # should we recompute the loglik_gp_old value here?
      # we currently recompute loglik_gp_old here!
      prediction_sample, sigma_pred = predict(Theta[:,r-1], gp, pred_method,est_method,noisy_est, true)
      loglik_gp_old = prediction_sample[1]
      loglik_gp_old_std = sigma_pred[1]

      a_gp = loglik_gp_new + prior_log_star +  jacobian_log_star -  (loglik_gp_old - prior_log_old - jacobian_log_old)

      accept = log(rand()) < a_gp # calc accept

      # store accaptance probability
      #accept_prob_log[1, r] = a_gp

      if !accept # early-accept

        # keep old values
        nbr_early_rejections = nbr_early_rejections + 1
        Theta[:,r] = Theta[:,r-1]
        loglik[r] = loglik[r-1]
        # adaptation of covaraince matrix for the proposal distribution
        # adaptation(problem.adaptive_update, adaptive_update_params, Theta, r,a_gp)

      else

        # stage 2 usign A-DA

        nbr_second_stage += 1

        # A-DA
        u_log_hat = log(rand())

        if loglik_gp_old < loglik_gp_new #&& loglik_gp_new_std < std_limit && loglik_gp_old_std < std_limit

          nbr_case_13 += 1

          # select case 1 or 3
          if selectcase1or3(casemodel, theta_star, loglik_gp_new, loglik_gp_old) == 1

            # case 1

            nbr_case_1 += 1

            if u_log_hat < loglik_gp_old - loglik_gp_new
              nbr_second_stage_accepted = nbr_second_stage_accepted+1
              nbr_split_accaptance_region_early_accept = nbr_split_accaptance_region_early_accept+1
              Theta[:,r] = theta_star # update chain with proposal
              loglik[r] = NaN
              accept_vec[r] = 1
            else

              loglik_star = pf_paralell(Z, theta_star,theta_known,N,dt,dt_U,nbr_x0, nbr_x,subsample_interval,false,false, nbr_of_proc,loglik_vec)
              loglik_old = pf_paralell(Z, Theta[:,r-1],theta_known,N,dt,dt_U,nbr_x0, nbr_x,subsample_interval,false,false, nbr_of_proc,loglik_vec)

              nbr_eval_pf += 1
              nbr_case_pf_1 += 1

              a_log = (loglik_star + loglik_gp_old)  -  (loglik_old + loglik_gp_new)
              accept = u_log_hat < a_log # calc accaptance decision

              if accept # the proposal is accapted
                nbr_second_stage_accepted = nbr_second_stage_accepted+1
                Theta[:,r] = theta_star # update chain with proposal
                loglik[r] = NaN
                accept_vec[r] = 1
              else
                Theta[:,r] = Theta[:,r-1] # keep old values
                loglik[r] = loglik[r-1]
              end
            end
          else

            nbr_case_3 += 1

            # run case 3
            if u_log_hat > loglik_gp_old - loglik_gp_new

              nbr_split_accaptance_region_early_reject = nbr_split_accaptance_region_early_reject+1
              Theta[:,r] = Theta[:,r-1] # keep old values
              loglik[r] = loglik[r-1]
              #accept_vec[r] = 1

            else

              loglik_star = pf_paralell(Z, theta_star,theta_known,N,dt,dt_U,nbr_x0, nbr_x,subsample_interval,false,false, nbr_of_proc,loglik_vec)
              loglik_old = pf_paralell(Z, Theta[:,r-1],theta_known,N,dt,dt_U,nbr_x0, nbr_x,subsample_interval,false,false, nbr_of_proc,loglik_vec)

              nbr_eval_pf += 1
              nbr_case_pf_3 += 1

              a_log = (loglik_star + loglik_gp_old)  -  (loglik_old + loglik_gp_new)
              accept = u_log_hat < a_log # calc accaptance decision

              if accept # the proposal is accapted
                nbr_second_stage_accepted = nbr_second_stage_accepted+1
                Theta[:,r] = theta_star # update chain with proposal
                loglik[r] = NaN
                accept_vec[r] = 1
              else
                Theta[:,r] = Theta[:,r-1] # keep old values
                loglik[r] = loglik[r-1]
              end
            end
          end
        else

          # select case 2 or 4
          nbr_case_24 += 1

          if selectcase2or4(casemodel, theta_star, loglik_gp_new, loglik_gp_old) == 1

            # case 2
            nbr_case_2 += 1

            loglik_star = pf_paralell(Z, theta_star,theta_known,N,dt,dt_U,nbr_x0, nbr_x,subsample_interval,false,false, nbr_of_proc,loglik_vec)
            loglik_old = pf_paralell(Z, Theta[:,r-1],theta_known,N,dt,dt_U,nbr_x0, nbr_x,subsample_interval,false,false, nbr_of_proc,loglik_vec)

            nbr_eval_pf += 1
            nbr_case_pf_2 += 1

            a_log = (loglik_star + loglik_gp_old)  -  (loglik_old + loglik_gp_new)
            accept = u_log_hat < a_log # calc accaptance decision

            if accept # the proposal is accapted
              nbr_second_stage_accepted = nbr_second_stage_accepted+1
              Theta[:,r] = theta_star # update chain with proposal
              loglik[r] = NaN
              accept_vec[r] = 1
            else
              Theta[:,r] = Theta[:,r-1] # keep old values
              loglik[r] = loglik[r-1]
            end

          else

           # case 4
           nbr_case_4 += 1

           # accept directly
           nbr_second_stage_accepted = nbr_second_stage_accepted+1
           Theta[:,r] = theta_star # update chain with proposal
           loglik[r] = NaN
           accept_vec[r] = 1
         end
        end
      end
    end
  end

  time_da_part = toc()
  times = [time_pre_er time_fit_gp time_da_part]

  @printf "Ending ADA-GP-MCMC with adaptive RW estimating %d parameters\n" length(theta_true)
  @printf "Algorithm: %s\n" alg
  @printf "Adaptation algorithm: %s\n" typeof(problem.adaptive_update)
  @printf "Prior distribution: %s\n" problem.prior_dist.dist
  @printf "Particel filter: %s\n" pf_alg
  @printf "Time consumption:"
  @printf "Time pre-er:  %.0f\n" time_pre_er
  @printf "Time fit GP model:  %.0f\n" time_fit_gp
  @printf "Time er-part:  %.0f\n" time_da_part
  @printf "Number early-rejections: %d\n"  nbr_early_rejections
  @printf "Number of cases directly to ordinary MH: %d\n"  nbr_ordinary_mh

  @printf "Number split acceptance region early-accept: %d\n"  nbr_split_accaptance_region_early_accept
  @printf "Number split acceptance region early-reject: %d\n"  nbr_split_accaptance_region_early_reject

  @printf "Number cases in DA part: %d\n"  nbr_run_DA
  @printf "Number cases in second stage: %d\n"  nbr_second_stage
  @printf "Number accepted in second stage: %d\n"  nbr_second_stage_accepted

  @printf "Total number of evaluations of the particle filter: %d\n" nbr_eval_pf

  @printf "Acceptance rate for ordinary MH accepted: %.4f\n"  nbr_ordinary_mh_accapte/nbr_ordinary_mh*100
  @printf "Acceptance rate a_1: %.4f\n"  nbr_second_stage/nbr_run_DA*100
  @printf "Acceptance rate a_2: %.4f\n"  nbr_second_stage_accepted/nbr_second_stage*100
  @printf "Acceptance rate a_1*a_2: %.4f\n"  (nbr_second_stage/nbr_run_DA)*(nbr_second_stage_accepted/nbr_second_stage)*100
  @printf "Acceptance rate a (for entier algorithm): %.4f\n" sum(accept_vec)/R*100

  @printf "Number case 1 or 3: %d\n"  nbr_case_13
  @printf "Number case 2 or 4: %d\n"  nbr_case_24

  @printf "Number case 1: %d, %.4f %% of all cases\n"  nbr_case_1 nbr_case_1/nbr_case_13
  @printf "Number case 2: %d, %.4f %% of all cases\n"  nbr_case_2 nbr_case_2/nbr_case_24
  @printf "Number case 3: %d, %.4f %% of all cases\n"  nbr_case_3 nbr_case_3/nbr_case_13
  @printf "Number case 4: %d, %.4f %% of all cases\n"  nbr_case_4 nbr_case_4/nbr_case_24

  @printf "Number pf runs in case 1: %d, prob pf given case 1 %.4f %%\n"  nbr_case_pf_1 nbr_case_pf_1/nbr_case_1*100
  @printf "Number pf runs in case 2: %d, prob pf given case 2 %.4f %%\n"  nbr_case_pf_2 nbr_case_pf_2/nbr_case_2*100
  @printf "Number pf runs in case 3: %d, prob pf given case 3 %.4f %%\n"  nbr_case_pf_3 nbr_case_pf_3/nbr_case_3*100
  @printf "Number pf runs in case 4: %d, prob pf given case 4 %.4f %%\n"  nbr_case_pf_4 nbr_case_pf_4/nbr_case_4*100

  @printf "#####################################################################\n"

  assumption_list = []
  loglik_list = []

  # return resutls
  if return_run_info
    run_info =  [nbr_eval_pf;
                nbr_ordinary_mh;
                nbr_case_13;
                nbr_case_24;
                nbr_case_1;
                nbr_case_2;
                nbr_case_3;
                nbr_case_4;
                nbr_case_pf_1;
                nbr_case_pf_2;
                nbr_case_pf_3;
                nbr_case_pf_4]
    return return_gp_results(gp, Theta,loglik,accept_vec,prior_vec, compare_GP_PF, data_gp_pf,nbr_early_rejections, problem, adaptive_update_params,accept_prob_log,times), run_info
  else
    return return_gp_results(gp, Theta,loglik,accept_vec,prior_vec, compare_GP_PF, data_gp_pf,nbr_early_rejections, problem, adaptive_update_params,accept_prob_log,times)
  end

end

################################################################################
######               help functions                                        #####
################################################################################

doc"""
    return_results(Theta,loglik,accept_vec,prior_vec, problem,adaptive_update_params)

Constructs the return type of the resutls from the PMCMC and MCWM algorithm.
"""
function return_results(Theta,loglik,accept_vec,prior_vec, problem,adaptive_update_params)
  if typeof(problem.adaptive_update) == AMUpdate_gen
    return (Result(Theta, loglik, accept_vec, prior_vec), adaptive_update_params[6])
  elseif typeof(problem.adaptive_update) == AMUpdate_comp_w
    return (Result(Theta, loglik, accept_vec, prior_vec), adaptive_update_params[7])
  elseif typeof(problem.adaptive_update) == AMUpdate_gen_comp_w
    return (Result(Theta, loglik, accept_vec, prior_vec), adaptive_update_params[7])
  else
    return Result(Theta, loglik, accept_vec, prior_vec)
  end
end

doc"""
    return_gp_results(gp, Theta,loglik,accept_vec,prior_vec, compare_GP_PF, data_gp_pf,nbr_early_rejections, problem, adaptive_update_params,accept_prob_log,times)

Constructs the return type of the resutls from the gpPMCMC algorithm.
"""
function return_gp_results(gp, Theta,loglik,accept_vec,prior_vec, compare_GP_PF, data_gp_pf,nbr_early_rejections, problem, adaptive_update_params,accept_prob_log,times)
  if typeof(problem.adaptive_update) == AMUpdate_gen
    return (gpResult(Theta, loglik, accept_vec, prior_vec,compare_GP_PF,data_gp_pf,nbr_early_rejections), adaptive_update_params[6],gp,accept_prob_log,times)
#  elseif typeof(problem.adaptive_update) == AMUpdate_comp_w
#    return (gpResult(Theta, loglik, accept_vec, prior_vec,compare_GP_PF,data_gp_pf,nbr_early_rejections), adaptive_update_params[7],gp,accept_prob_log,times)
  else
    return gpResult(Theta, loglik, accept_vec, prior_vec,compare_GP_PF,data_gp_pf,nbr_early_rejections),gp,accept_prob_log,times
  end
end

doc"""
    set_nbr_cores(nbr_of_cores::Int64, pf_alg::String)

Sets the number of cores to use.
"""
function set_nbr_cores(nbr_of_cores::Int64, pf_alg::String)
  if length(workers()) == 1
    nbr_of_proc = nbr_of_cores
    addprocs(nbr_of_proc)
    if pf_alg == "parallel_bootstrap"
      @everywhere include("run_pf_paralell.jl")
    else
      error("The auxiliary particle filter is not implemented.")
    end
  else
    nbr_of_proc = length(workers())
  end
  return nbr_of_proc
end

doc"""
    evaluate_prior(theta_star, prior_parameters, dist_type)

Calculates the `log-prior` value for the prior distribution for the parameters `theta_star`.
"""
function  evaluate_prior(theta_star, prior_parameters, dist_type)

  # set start value for loglik
  log_prior = 0.

  if dist_type == "Uniform"
    for i = 1:length(theta_star)
      # Update loglik, i.e. add the loglik for each model paramter in theta
      log_prior = log_prior + log_unifpdf( theta_star[i], prior_parameters[i,1], prior_parameters[i,2] )
    end
  elseif dist_type == "Normal"
    for i = 1:length(theta_star)
      # Update loglik, i.e. add the loglik for each model paramter in theta
      log_prior = log_prior + log_normpdf(theta_star[i],prior_parameters[i,1],prior_parameters[i,2])
    end
  elseif dist_type == "nonlog"
    # add code to handle priors on non-log-scale!
    if length(theta_star) == 2
      for i = 1:length(theta_star)
        # the unknown parameters c and d both have normal prior dists
        log_prior = log_prior + log_normpdf( exp(theta_star[i]), prior_parameters[i,1], prior_parameters[i,2] )
      end
    elseif length(theta_star) == 3
      # the unknown parameter A has a inv-gamma prior dist
      log_prior = log_prior + log_invgampdf(exp(theta_star[1]), prior_parameters[1,1], prior_parameters[1,2])
      for i = 2:3
        # The unknown parameters c and d both have normal prior dists
        log_prior = log_prior + log_normpdf( exp(theta_star[i]), prior_parameters[i,1], prior_parameters[i,2] )
      end
    elseif length(theta_star) == 5
      for i in [1 2 5]
        # The unknown parameters Κ,Γ and sigma both have gamma prior dists
        log_prior = log_prior + log_gampdf( exp(theta_star[i]), prior_parameters[i,1], prior_parameters[i,2] )
      end
      for i = 3:4
        # The unknown parameters c and d both have normal prior dists
        log_prior = log_prior + log_normpdf( exp(theta_star[i]), prior_parameters[i,1], prior_parameters[i,2] )
      end
    elseif length(theta_star) == 6
      # The unknown parameter A has a gaminv prior dist
      log_prior = log_prior + log_invgampdf(exp(theta_star[1]), prior_parameters[1,1], prior_parameters[1,2])
      for i = 2:3
        # The unknown parameters c and d both have normal prior dists
        log_prior = log_prior + log_normpdf( exp(theta_star[i]), prior_parameters[i,1], prior_parameters[i,2] )
      end
      for i = 4:length(theta_star)
        # The unknown parameters p1, p2 and sigma both have gamma prior dists
        log_prior = log_prior + log_gampdf( exp(theta_star[i]), prior_parameters[i,1], prior_parameters[i,2] )
      end
    elseif length(theta_star) == 7
      for i = [1,2,5,6,7]
        # The unknown parameters Κ,Γ,power1,power2 and sigma all have gamma prior dists
        log_prior = log_prior + log_gampdf( exp(theta_star[i]), prior_parameters[i,1], prior_parameters[i,2] )
      end
      for i = [3,4]
        # The unknown parameters c and d have normal prior dists
        log_prior = log_prior + log_normpdf( exp(theta_star[i]), prior_parameters[i,1], prior_parameters[i,2] )
      end

    elseif length(theta_star) == 8
      # the unknown parameter A has a inv-gamma prior dist
      log_prior = log_prior + log_invgampdf(exp(theta_star[3]), prior_parameters[3,1], prior_parameters[3,2])
      for i = [1,2,5,6,7]
        # The unknown parameters Κ,Γ,power1,power2 and sigma all have gamma prior dists
        log_prior = log_prior + log_gampdf( exp(theta_star[i]), prior_parameters[i,1], prior_parameters[i,2] )
      end
      for i = [3,4]
        # The unknown parameters c and d have normal prior dists
        log_prior = log_prior + log_normpdf( exp(theta_star[i]), prior_parameters[i,1], prior_parameters[i,2] )
      end
    elseif length(theta_star) == 4
      for i = 1:2
        # The unknown parameters Κ and Γ both have gammma prior dists
        log_prior = log_prior + log_invgampdf( exp(theta_star[i]), prior_parameters[i,1], prior_parameters[i,2] )
      end
      for i = 3:4
        # The unknown parameters c and d both have normal prior dists
        log_prior = log_prior + log_normpdf( exp(theta_star[i]), prior_parameters[i,1], prior_parameters[i,2] )
      end

    else
      for i in [3 6]
        # The unknown parameters A and g have a gaminv prior dist
        log_prior = log_prior + log_invgampdf(exp(theta_star[i]), prior_parameters[i,1], prior_parameters[i,2])
      end
      for i in [4 5]
        # The unknown parameters c and d both have normal prior dists
        log_prior = log_prior + log_normpdf( exp(theta_star[i]), prior_parameters[i,1], prior_parameters[i,2] )
      end
      for i in [1 2 7 8 9]
        # The unknown parameters Κ,Γ,power1,power2 and sigma all have gamma prior dists
        log_prior = log_prior + log_gampdf( exp(theta_star[i]), prior_parameters[i,1], prior_parameters[i,2] )
      end

    end

  end

  return log_prior # return log_lik

end

doc"""
    jacobian(theta::Vector, parameter_transformation::String)

Returnes log-Jacobian for transformation of proposal space.
"""
function jacobian(theta::Vector)

  return sum(theta)

end

doc"""
    log_gampdf(x::Float64, mu::Float64, sigma::Float64)

Computes log(normpdf(x,a,b)) without normalizing constant.

"""
function log_gampdf(x::Float64, a::Float64, b::Float64)
  #return -a*log(b) - log(gamma(a)) - (a-1)*log(x) - x/b
  return -(a-1)*log(x) - x/b

end

doc"""
    log_invgampdf(x::Float64, mu::Float64, sigma::Float64)

Computes log(log_invgampdf(x,a,b)) constant.

"""
function log_invgampdf(x::Float64, a::Float64, b::Float64)
  #return a*log(b) - log(gamma(a)) + (-a-1)*log(x) - b/x
  return (-a-1)*log(x) - b/x
end

doc"""
    log_unifpdf(x::Float64, a::Float64,b::Float64)

Computes log(unifpdf(x,a,b)).

"""
function log_unifpdf(x::Float64, a::Float64,b::Float64)
  if  x >= a && x<= b
    return -log(b-a)
  else
    return log(0)
  end
end

doc"""
    log_normpdf(x::Float64, mu::Float64, sigma::Float64)

Computes log(normpdf(x,a,b)) without normalizing constant.

"""
function log_normpdf(x::Float64, mu::Float64, sigma::Float64)
  #return -log(sigma) - (x-mu)^2/(2*sigma^2)
  return -(x-mu)^2/(2*sigma^2)
end

doc"""
    pf_paralell(Z::Array{Float64},theta::Array{Float64},theta_known::Array{Float64}, N::Int64,dt::Float64, nbr_x0::Int64, nbr_x::Int64,subsample_interval::Int64,print_on::Bool, store_weigths::Bool, nbr_of_proc::Int64)

Runs the bootstrap filter to estimimate the `log-likelihood`
log(P(`Z`|`theta`))
"""
function pf_paralell(Z::Array{Float64},theta::Array{Float64},theta_known::Array{Float64},
  N::Int64,dt::Float64, dt_U::Float64, nbr_x0::Int64, nbr_x::Int64,
  subsample_interval::Int64,print_on::Bool, return_weigths_and_particels::Bool,
  nbr_of_proc::Int64, loglik::SharedArray)

  # set parameters
  (Κ, Γ, A, A_sign,B,c,d,g,f,power1,power2,sigma) = set_parameters(theta, theta_known,length(theta))

  # compute value for A
  A = A*A_sign

  # set value for constant b function
  b_const = sqrt(2.*sigma^2 / 2.)

  # set values needed for calculations in Float64
  (subsample_interval_calc, nbr_x0_calc, nbr_x_calc,N_calc) = map(Float64, (subsample_interval, nbr_x0, nbr_x,N))

  if !return_weigths_and_particels

    # run nbr_of_proc parallel estimations

    #loglik = SharedArray{Float64}(nbr_of_proc)

    @sync begin

    @parallel for i = 1:nbr_of_proc
      loglik[i] = run_pf_paralell(Z,theta,theta_known, N, N_calc, dt, dt_U, nbr_x0, nbr_x0_calc, nbr_x, nbr_x_calc, subsample_interval, subsample_interval_calc, print_on, return_weigths_and_particels, Κ, Γ, A, B, c, d, f, g, power1, power2, b_const)
    end

    end

    #%println(loglik)
    #println(logsumexp(loglik)-log(nbr_of_proc))

    #println(loglik)

    return logsumexp(loglik)-log(nbr_of_proc)

    #=
    loglik = @parallel (+) for i = 1:nbr_of_proc
      run_pf_paralell(Z,theta,theta_known, N, N_calc, dt, dt_U, nbr_x0, nbr_x0_calc, nbr_x, nbr_x_calc, subsample_interval, subsample_interval_calc, print_on, return_weigths_and_particels, Κ, Γ, A, B, c, d, f, g, power1, power2, b_const)
    end

    # return the average of the nbr_of_proc estimations
    return loglik/nbr_of_proc
    =#
  else
    # run nbr_of_proc parallel estimations
    (loglik,weigths, particels) = run_pf_paralell(Z,theta,theta_known, N, N_calc, dt, dt_U, nbr_x0, nbr_x0_calc, nbr_x, nbr_x_calc, subsample_interval, subsample_interval_calc, print_on, return_weigths_and_particels, Κ, Γ, A, B, c, d, f, g, power1, power2, b_const)

    # return the average of the nbr_of_proc estimations
    return loglik, weigths, particels
  end
end

doc"""
    set_parameters(theta::Array{Float64}, theta_known::Array{Float64}, nbr_of_unknown_parameters::Int64)

Sets the model parameters according to the nbr of unknown parameters. Help function for pf.

"""
function set_parameters(theta::Array{Float64}, theta_known::Array{Float64}, nbr_of_unknown_parameters::Int64)

  # set parameter values
  if nbr_of_unknown_parameters == 8
    # estiamte Κ, Γ,A,c,d,power1,power2,sigma
    (Κ,Γ,A,c,d,power1,power2,sigma) = exp.(theta)
    (A_sign,B,f,g) = theta_known
  elseif nbr_of_unknown_parameters == 7
    # estimate Κ,Γ,c,d,power1,power2,sigma
    (Κ,Γ,c,d,power1,power2,sigma) = exp.(theta)
    (A,A_sign,B,f,g) = theta_known
  elseif nbr_of_unknown_parameters == 6
    # estiamte A,c,d,power1,power2,sigma
    (A,c,d,power1,power2,sigma) = exp.(theta)
    (Κ, Γ, A_sign,B,f,g) = theta_known
  elseif nbr_of_unknown_parameters == 5
    # estimate Κ Γ c d sigma
    (Κ,  Γ, c, d, sigma) = exp.(theta)
    (A,A_sign, B, f, g, power1, power2) = theta_known
  elseif nbr_of_unknown_parameters == 4
    # estimate Κ Γ c d
    (Κ, Γ,c,d) = exp.(theta)
    (A,A_sign,B,f,g,power1,power2, sigma) = theta_known
  elseif nbr_of_unknown_parameters == 2
    # estimating c and d
    (c,d) = exp.(theta)
    (Κ, Γ, A, A_sign, B, f, g, power1, power2, sigma) = theta_known
  elseif nbr_of_unknown_parameters == 3
    # estimate A,c,d
    (A, c, d) = exp.(theta)
    (Κ, Γ, A_sign, B, f, g, power1, power2, sigma) = theta_known
  else
    # estimate all parameters. i.e. estimate Κ, Γ,A,c,d,g,power1,power2,sigma
    (Κ, Γ, A,c,d,g,power1,power2,sigma) = exp.(theta)
    (A_sign,B,f) = theta_known
  end

  return (Κ, Γ, A, A_sign, B,c,d,g,f,power1,power2,sigma)

end

################################################################################
###            Functions for particle filter diagnostics                    ####
################################################################################

doc"""
    pf_diagnostics(problem::Problem, nbr_iterations::Int64, theta_compute::Array)

Runs the boostrap filter or the auxiliar particle filter for nbr_iterations times.

# Inputs

* `problem`: the type decribing the problem
* `nbr_iterations`: number of iterations of the particel filters

# Output

* `loglik_matrix`: A nbr_iterations x nbr_iterations matrix with the computed log-likelihood values

"""
function pf_diagnostics(problem::Problem, nbr_iterations::Int64, theta_compute::Array)

  # data
  Z = problem.data.Z

  # algorithm parameters
  R = problem.alg_param.R
  N = problem.alg_param.N
  burn_in = problem.alg_param.burn_in
  pf_alg = problem.alg_param.pf_alg
  nbr_x0  = problem.alg_param.nbr_x0
  nbr_x = problem.alg_param.nbr_x
  subsample_interval = problem.alg_param.subsample_int
  dt = problem.alg_param.dt
  nbr_of_cores = problem.alg_param.nbr_of_cores

  # model parameters
  theta_true = problem.model_param.theta_true
  theta_known = problem.model_param.theta_known
  theta_0 = problem.model_param.theta_0

  # pre-allocate matrix for storing computed loglik values
  loglik_vec = zeros(nbr_iterations,1) # pf apf

  @printf "Running the particle filter %d times \n" nbr_iterations

  # set nbr of cores to use
  nbr_of_proc = set_nbr_cores(nbr_of_cores, pf_alg)

  if pf_alg == "parallel_bootstrap"
    for j = 1:nbr_iterations
      loglik_vec[j] = pf_paralell(Z, theta_compute,theta_known,N,dt,nbr_x0, nbr_x,subsample_interval,false,false, nbr_of_proc)
    end

    # run nbr_of_proc parallel estimations
    (loglik, weigths, particels) = pf_paralell(Z, theta_compute,theta_known,N,dt,nbr_x0, nbr_x,subsample_interval,true,true, nbr_of_proc)

  elseif pf_alg == "parallel_apf"

    error("The auxiliary particle filter is not implemented.")

  end

  return loglik_vec, weigths, particels

end

max_stars_repo_path: DWPSDE model/mcmc.jl

max_stars_repo_name: SamuelWiqvist/adamcmcpaper

max_issues_repo_path: DWPSDE model/mcmc.jl

max_issues_repo_name: SamuelWiqvist/adamcmcpaper

max_forks_repo_path: DWPSDE model/mcmc.jl

max_forks_repo_name: SamuelWiqvist/adamcmcpaper


text:
#=
# A parameter to be used in the weak form.
# It is an array of Expr.
# Each represents a function of x,y,z,t, coefficients, or known variables
# Values are in an array to allow vector/tensor valued parameters
=#

struct Parameter
    symbol::Symbol          # symbol used in expressions
    index::Int              # index in the Finch list of parameters
    type::String            # constants for SCALAR, VECTOR, etc.
    value::Array            # An array of either constant values(numbers) or Expr expressions
end

ext: jl

max_stars_repo_name: paralab/Finch

max_issues_repo_name: paralab/Finch

max_forks_repo_name: paralab/Finch


text:
"""
    couple_models(aero::LiftingLine, dyn::RigidBody)

Create a coupled model using a lifting line aerodynamic model and a rigid body
dynamics model.  This model introduces additional parameters corresponding to
the length, position, and orientation of each lifting line element ``L, p_e,
e_1, e_2, e_3``, followed by the freestream air density ``rho_\\infty``, rigid
body inertial properties ``m, Ixx, Iyy, Izz, Ixz, Ixy, Iyz``, and the additional
forces/moments applied on the body ``F_x, F_y, F_z, M_x, M_y, M_z``.

**NOTE: When using this model, the local frame for each lifting line element should be
oriented with the x-axis in the chordwise direction, the y-axis in the spanwise
direction (out the right wing), and the z-axis in the airfoil normal direction**
"""
couple_models(aero::LiftingLine, dyn::RigidBody) = (aero, dyn)

# --- Traits --- #

function number_of_additional_parameters(::Type{<:LiftingLine{NA,TA}}, ::Type{<:RigidBody}) where {NA,TA}
    return 13*NA + 14
end

coupling_inplaceness(::Type{<:LiftingLine}, ::Type{<:RigidBody}) = InPlace()

function coupling_rate_jacobian_type(::Type{LiftingLine{NA,TA}}, ::Type{<:RigidBody}) where {NA,TA}
    model_types = TA.parameters
    if all(isempty.(coupling_rate_jacobian_type.(model_types, Ref(LiftingLineSection))))
        return Empty()
    elseif all(iszero.(coupling_rate_jacobian_type.(model_types, Ref(LiftingLineSection))))
        return Zeros()
    elseif all(isidentity.(coupling_rate_jacobian_type.(model_types, Ref(LiftingLineSection))))
        return Identity()
    elseif all(isinvariant.(coupling_rate_jacobian_type.(model_types, Ref(LiftingLineSection))))
        return Invariant()
    elseif all(isconstant.(coupling_rate_jacobian_type.(model_types, Ref(LiftingLineSection))))
        return Constant()
    elseif all(islinear.(coupling_rate_jacobian_type.(model_types, Ref(LiftingLineSection))))
        return Linear()
    else
        return Nonlinear()
    end
end

function coupling_state_jacobian_type(::Type{LiftingLine{NA,TA}}, ::Type{<:RigidBody}) where {NA,TA}
    model_types = TA.parameters
    if all(isempty.(coupling_state_jacobian_type.(model_types, Ref(LiftingLineSection))))
        return Empty()
    elseif all(iszero.(coupling_state_jacobian_type.(model_types, Ref(LiftingLineSection))))
        return Zeros()
    elseif all(isidentity.(coupling_state_jacobian_type.(model_types, Ref(LiftingLineSection))))
        return Identity()
    elseif all(isinvariant.(coupling_state_jacobian_type.(model_types, Ref(LiftingLineSection))))
        return Invariant()
    elseif all(isconstant.(coupling_state_jacobian_type.(model_types, Ref(LiftingLineSection))))
        return Constant()
    elseif all(islinear.(coupling_state_jacobian_type.(model_types, Ref(LiftingLineSection))))
        return Linear()
    else
        return Nonlinear()
    end
end

coupling_parameter_jacobian_type(::Type{<:LiftingLine}, ::Type{<:RigidBody}) = Nonlinear()

function coupling_time_gradient_type(::Type{LiftingLine{NA,TA}}, ::Type{<:RigidBody}) where {NA,TA}
    model_types = TA.parameters
    if all(isempty.(coupling_time_gradient_type.(model_types, Ref(LiftingLineSection))))
        return Empty()
    elseif all(iszero.(coupling_time_gradient_type.(model_types, Ref(LiftingLineSection))))
        return Zeros()
    elseif all(isinvariant.(coupling_time_gradient_type.(model_types, Ref(LiftingLineSection))))
        return Invariant()
    elseif all(isconstant.(coupling_time_gradient_type.(model_types, Ref(LiftingLineSection))))
        return Constant()
    elseif all(islinear.(coupling_time_gradient_type.(model_types, Ref(LiftingLineSection))))
        return Linear()
    else
        return Nonlinear()
    end
end

# --- Methods --- #

function get_coupling_inputs!(y, aero::LiftingLine{NA,TA}, dyn::RigidBody, dx, x, p, t) where {NA,TA}

    # extract number of state variables, inputs, and parameters
    nua = number_of_states(aero) # number of aerodynamic states
    nud = number_of_states(dyn) # number of rigid body states
    nya = number_of_inputs(aero) # number of aerodynamic inputs
    nyd = number_of_inputs(dyn) # number of rigid body inputs
    npa = number_of_parameters(aero) # number of aerodynamic parameters
    npd = number_of_parameters(dyn) # number of rigid body parameters
    npadd = number_of_additional_parameters(aero, dyn) # number of additional parameters for coupled model

    # get indices for state variables, inputs, and parameters
    iua = 1:nua # indices of aerodynamic states
    iya = 1:nya # indices of aerodynamic inputs
    ipa = 1:npa # indices of aerodynamic parameters
    iud = nua + 1 : nua + nud # indices of rigid body states
    iyd = nya + 1 : nya + nyd # indices of rigid body inputs
    ipd = npa + 1 : npa + npd # indices of rigid body parameters
    ipadd = npa + npd + 1 : npa + npd + npadd # indices of additional parameters
    iuas = state_indices(aero.models) # indices of aerodynamic states for each section
    iyas = input_indices(aero.models) # indices of aerodynamic inputs for each section
    ipas = parameter_indices(aero.models) # indices of aerodynamic inputs for each section

    # separate rates, states, inputs, and parameters
    dua = view(dx, iua) # aerodynamic rates
    ua = view(x, iua) # aerodynamic states
    ya = view(y, iya) # aerodynamic inputs
    pa = view(p, ipa) # aerodynamic parameters
    dud = view(dx, iud) # rigid body rates
    ud = view(x, iud) # rigid body states
    yd = view(y, iyd) # rigid body inputs
    pd = view(p, ipd) # rigid body parameters
    padd = view(p, ipadd) # additional parameters
    duas = view.(Ref(dua), iuas) # aerodynamic rates for each section
    uas = view.(Ref(ua), iuas) # aerodynamic states for each section
    yas = view.(Ref(ya), iyas) # aerodynamic inputs for each section
    pas = view.(Ref(pa), ipas) # aerodynamic parameters for each section

    # rigid body rates and states
    dxr, dyr, dzr, dϕr, dθr, dψr, dur, dvr, dwr, dpr, dqr, drr = dud
    xr, yr, zr, ϕr, θr, ψr, ur, vr, wr, pr, qr, rr = ud

    # global parameters
    ρ = padd[13*NA+1]
    m = padd[13*NA+2]
    Ixx = padd[13*NA+3]
    Iyy = padd[13*NA+4]
    Izz = padd[13*NA+5]
    Ixz = padd[13*NA+6]
    Ixy = padd[13*NA+7]
    Iyz = padd[13*NA+8]
    Fbx = padd[13*NA+9]
    Fby = padd[13*NA+10]
    Fbz = padd[13*NA+11]
    Mbx = padd[13*NA+12]
    Mby = padd[13*NA+13]
    Mbz = padd[13*NA+14]

    # body linear and angular velocity
    V = SVector(ur, vr, wr)
    Ω = SVector(pr, qr, rr)

    # body linear and angular acceleration
    dV = SVector(dur, dvr, dwr)
    dΩ = SVector(dpr, dqr, drr)

    # aerodynamic forces and moments
    Fa = @SVector zeros(3)
    Ma = @SVector zeros(3)
    for i = 1:NA
        # models for this section
        section_aero = aero.models[i]
        section_stru = LiftingLineSection()
        section_models = (section_aero, section_stru)

        # model dimensions for this section
        Nuai = number_of_states(section_aero)
        Nyai = number_of_inputs(section_aero)
        Npai = number_of_parameters(section_aero)
        Nusi = number_of_states(section_stru)
        Nysi = number_of_inputs(section_stru)
        Npsi = number_of_parameters(section_stru)

        # section properties
        offset = 13*(i-1)
        ΔL = padd[offset+1] # element length
        pe = SVector(padd[offset+2], padd[offset+3], padd[offset+4]) # local frame position
        e1 = SVector(padd[offset+5], padd[offset+6], padd[offset+7]) # local frame x-axis
        e2 = SVector(padd[offset+8], padd[offset+9], padd[offset+10]) # local frame y-axis
        e3 = SVector(padd[offset+11], padd[offset+12], padd[offset+13]) # local frame z-axis
        Cab = [e1 e2 e3] # local to body frame transformation matrix
        Cba = Cab' # body to local frame transformation matrix
        vi = -Cba*(V + cross(Ω, pe)) # local linear freestream velocity
        ωi = Cba*Ω # local angular freestream velocity
        dvi =-Cba*(dV + cross(dΩ, pe)) # local linear freestream acceleration
        dωi = Cba*dΩ # local angular freestream acceleration

        # section rates
        duai = SVector{Nuai}(duas[i]) # aerodynamic state variables
        dusi = vcat(dvi, dωi) # rigid body state variables
        dui = vcat(duai, dusi)

        # section states
        uai = SVector{Nuai}(uas[i]) # aerodynamic state variables
        usi = vcat(vi, ωi) # rigid body state variables
        ui = vcat(uai, usi)

        # section parameters
        pai = SVector{Npai}(pas[i]) # aerodynamic parameters
        psi = SVector(ρ) # rigid body parameters
        pi = vcat(pai, psi)

        # section inputs
        yi = get_coupling_inputs(section_models, dui, ui, pi, t)

        # separate inputs
        yai = view(yi, 1:Nyai)
        ysi = view(yi, Nyai + 1 : Nyai + Nysi)

        # save local inputs
        y[iya[iyas[i]]] .= yai # aerodynamic inputs

        # section aerodynamic loads (in body frame)
        Fi = Cab*ΔL*SVector(ysi[1], ysi[2], ysi[3])
        Mi = Cab*ΔL*SVector(ysi[4], ysi[5], ysi[6])

        # add loads to total aerodynamic loads
        Fa += Fi
        Ma += cross(pe, Fi) + Mi
    end

    # save rigid body inputs
    y[iyd[1]] = m
    y[iyd[2]] = Ixx
    y[iyd[3]] = Iyy
    y[iyd[4]] = Izz
    y[iyd[5]] = Ixz
    y[iyd[6]] = Ixy
    y[iyd[7]] = Iyz
    y[iyd[8]] = Fa[1] + Fbx
    y[iyd[9]] = Fa[2] + Fby
    y[iyd[10]] = Fa[3] + Fbz
    y[iyd[11]] = Ma[1] + Mbx
    y[iyd[12]] = Ma[2] + Mby
    y[iyd[13]] = Ma[3] + Mbz

    return y
end

# --- Performance Overloads --- #

# TODO

# --- convenience methods --- #

function set_additional_parameters!(padd, model::LiftingLine{NA,TA},
    dyn::RigidBody; section_parameters, rho, m, Ixx, Iyy, Izz, Ixz, Ixy, Iyz,
    Fbx, Fby, Fbz, Mbx, Mby, Mbz) where {NA,TA}

    for i = 1:NA
        padd[13*(i-1)+1] = section_parameters[i].L
        padd[13*(i-1)+2 : 13*(i-1)+4] = section_parameters[i].pe
        padd[13*(i-1)+5 : 13*(i-1)+7] = section_parameters[i].e1
        padd[13*(i-1)+8 : 13*(i-1)+10] = section_parameters[i].e2
        padd[13*(i-1)+11 : 13*(i-1)+13] = section_parameters[i].e3
    end

    padd[13*NA+1] = rho
    padd[13*NA+2] = m
    padd[13*NA+3] = Ixx
    padd[13*NA+4] = Iyy
    padd[13*NA+5] = Izz
    padd[13*NA+6] = Ixz
    padd[13*NA+7] = Ixy
    padd[13*NA+8] = Iyz
    padd[13*NA+9] = Fbx
    padd[13*NA+10] = Fby
    padd[13*NA+11] = Fbz
    padd[13*NA+12] = Mbx
    padd[13*NA+13] = Mby
    padd[13*NA+14] = Mbz

    return padd
end

function separate_additional_parameters(model::LiftingLine{NA,TA},
    dyn::RigidBody, padd) where {NA,TA}

    # section parameters
    return (
        section_parameters = ntuple(i -> (
            L = padd[13*(i-1)+1], # element length
            pe = view(padd, 13*(i-1)+2 : 13*(i-1)+4), # local frame position
            e1 = view(padd, 13*(i-1)+5 : 13*(i-1)+7), # local frame x-axis
            e2 = view(padd, 13*(i-1)+8 : 13*(i-1)+10), # local frame y-axis
            e3 = view(padd, 13*(i-1)+11 : 13*(i-1)+13), # local frame z-axis
            ), NA),
        rho = padd[13*NA+1],
        m = padd[13*NA+2],
        Ixx = padd[13*NA+3],
        Iyy = padd[13*NA+4],
        Izz = padd[13*NA+5],
        Ixz = padd[13*NA+6],
        Ixy = padd[13*NA+7],
        Iyz = padd[13*NA+8],
        Fbx = padd[13*NA+9],
        Fby = padd[13*NA+10],
        Fbz = padd[13*NA+11],
        Mbx = padd[13*NA+12],
        Mby = padd[13*NA+13],
        Mbz = padd[13*NA+14],
        )
end


text:
using Test
using KernelFunctions
using Random: seed!

seed!(42)

dims = (10,5)
X = rand(dims...)
x = rand(dims[1])
s = 3.0
v1 = vcat(3.0,4.0*ones(dims[2]-1))
v2 = vcat(3.0,4.0*ones(dims[1]-1))
P = rand(5,10)
sdims = [1,2,3]
f(x) = sin.(x)

@testset "Transform Test" begin
    ## Test Scale Transform
    @testset "ScaleTransform" begin
        t = ScaleTransform(s)
        vt1 = ScaleTransform(v1)
        vt2 = ScaleTransform(v2)
        @test all(KernelFunctions.transform(t,X).==s*X)
        @test all(KernelFunctions.transform(vt1,X,1).==v1'.*X)
        @test all(KernelFunctions.transform(vt2,X,2).==v2.*X)
    end
    ## Test LowRankTransform
    @testset "LowRankTransform" begin
        tp = LowRankTransform(P)
        @test all(KernelFunctions.transform(tp,X,2).==P*X)
        @test all(KernelFunctions.transform(tp,x).==P*x)
        @test all(KernelFunctions.params(tp).==P)
        P2 = rand(5,10)
        KernelFunctions.set!(tp,P2)
        @test all(tp.proj.==P2)
    end
    ## Test FunctionTransform
    @testset "FunctionTransform" begin
        tf = FunctionTransform(f)
        KernelFunctions.transform(tf,X,1)
        @test all(KernelFunctions.transform(tf,X,1).==f(X))
    end
    ## Test SelectTransform
    @testset "SelectTransform" begin
        ts = SelectTransform(sdims)
        @test all(KernelFunctions.transform(ts,X,2).==X[sdims,:])
        @test all(KernelFunctions.transform(ts,x).==x[sdims])
        @test all(KernelFunctions.params(ts).==sdims)
        sdims2 = [2,3,5]
        KernelFunctions.set!(ts,sdims2)
        @test all(ts.select.==sdims2)
    end
    ## Test ChainTransform
    @testset "ChainTransform" begin
        t = ScaleTransform(s)
        tp = LowRankTransform(P)
        tf = FunctionTransform(f)
        tchain = ChainTransform([t,tp,tf])
        @test all(KernelFunctions.transform(tchain,X,2).==f(P*(s*X)))
        @test all(KernelFunctions.transform(tchain,X,2).==
                    KernelFunctions.transform(tf∘tp∘t,X,2))
    end
end


text:
include("test_RotationMatrix.jl")
include("test_Quaternion.jl")
include("test_Interpolation.jl")


text:
using Weave
using Test

cleanup = true

#Test hold and term options
weave("documents/test_hold.mdw", doctype="pandoc", plotlib="Gadfly")
result =  read("documents/test_hold.md", String)
ref =  read("documents/test_hold_ref.md", String)
@test result == ref
cleanup && rm("documents/test_hold.md")

#Test setting and restoring chunk options
Weave.weave("documents/default_opts.noweb", doctype = "tex")
result =  read("documents/default_opts.tex", String)
ref =  read("documents/default_opts_ref.tex", String)
@test result == ref
cleanup && rm("documents/default_opts.tex")


text:
struct ASTNode
	func::Symbol
	children::Vector{ASTNode}
	args::Dict{Symbol,Any}
	source::Vector{AbstractString}
	refs::Vector{ASTNode}
	warnings::Vector{AbstractString}
	reachable::Bool
	rulename::Symbol
end

function ASTNode(func::Symbol, children::Array{ASTNode} = (ASTNode)[], args::Dict{Symbol,Any} = Dict{Symbol,Any}())
	ASTNode(func, children, args, (AbstractString)[], (ASTNode)[], (AbstractString)[], false, :unnamed)
end

describe_ast_node(node::ASTNode) = "$(node.func)" * (haskey(node.args,:name) ? " " * "$(node.args[:name])" : "")

import Base.show
function show(io::IO, astnode::ASTNode, indent=0)
  indentstr = " "^indent
  print(io, "$(indentstr)$(astnode.func)")
  for arg in astnode.args
    if typeof(arg[2]) <: Array
      print(io, " $(arg[1])=>$(join(arg[2],","))")
    elseif typeof(arg[2]) <: OrdinalRange
      range = arg[2]
      print(io, " $(arg[1])=>$(range[1]):$(range[end])")
    else
      print(io, " $(arg[1])=>$(arg[2])")
    end
  end
  for sourceline in astnode.source
    print(io, " SOURCE: $(sourceline)")
  end
  for ref in astnode.refs
    print(io, " REF: $(describe_ast_node(ref))")  # avoids issues of circular references, and lengthy output
  end
  for warning in astnode.warnings
    print(io, " WARNING: $(warning)")
  end
  if (!astnode.reachable)
    print(io, " UNREACHABLE")
  end
  if astnode.rulename != :unnamed
    print(io, " RULENAME: $(astnode.rulename)")
  end
  println(io)
  for child in astnode.children
    show(io, child, indent+2)
  end
end


text:
function reorderPhylo(args...; x = Phy,
    order::ASCIIString = "postorder", index_only = false)

    ORDER = ["cladewise", "postorder", "pruningwise"]
    io = [1:length(ORDER)][ORDER .== order]
    nb_edge = size(x["edge"])[1]

    if isempty(io)
        error("reorderPhylo: Wrong Order\n", "ambiguous order!")
    end

    if io == 1
        if index_only
            return [1:1:nb_edge]
        else
            return x
        end
    end

    nb_node = x["Nnode"]
    if nb_node == 1
        if index_only
            return [1:1:nb_edge]
        else
            return x
        end
    end

    nb_tip = length(x["tip_label"])
    if nb_node >= nb_tip
        error("reorderPhylo:badtree\n", "tree apparently badly conformed!")
    end

    if io == 3 # need to revise
        # x = reorder(x)
        # neworder = ccall((:neworder_pruningwise,"/Users/JINGZHAI/MkDocs/GitHub/phylogeneticTree/reorder_phylo.so"),
        #     RetType, (Int,Int,Array{Int64,1},Array{Int64,1},Int,Int), Int(nb_tip), Int(nb_node),Int(x["edge"][:,1],
        #         Int(x["edge"][:,2]),Int(nb_edge),Int(nb_edge))
    else
        neworder = zeros(Int, nb_edge)
        neworder = neworder_phylo(n=Int(nb_tip), e1=int(x["edge"][:,1]),
                 e2=int(x["edge"][:,2]), N=Int(nb_edge),
                 neworder=neworder,order=Int(io[1]))
    end

    if index_only
        return neworder
    else
        x["edge"] = x["edge"][neworder, :]

        if ! isempty(x["edge_length"])
            x["edge_length"] = x["edge"][neworder, :]
        end

        return x
    end
end


text:
# MIT license
# Copyright (c) Microsoft Corporation. All rights reserved.
# See LICENSE in the project root for full license information.

@testset "General" begin
    @testset "QuadraticRoots" begin
        Random.seed!(SEED)
        similarroots(r1, r2, x1, x2) = isapprox([r1, r2], [x1, x2], rtol = 1e-9) || isapprox([r2, r1], [x1, x2], rtol = 1e-9)

        for i in 1:10000
            r1, r2, scale = rand(3) .- 0.5
            a, b, c = scale .* (1, r1 + r2, r1 * r2)
            x1, x2 = quadraticroots(a, b, c)
            @test similarroots(-r1, -r2, x1, x2)
        end

        a, b, c = 1, 0, -1
        x1, x2 = quadraticroots(a, b, c)
        @test similarroots(1, -1, x1, x2)

        a, b, c = 1, 2, 1 # (x+1)(x+1)= x^2 + 2x + 1, double root at one
        x1, x2 = quadraticroots(a, b, c)
        @test similarroots(-1, -1, x1, x2)
    end # testset QuadraticRoots

    @testset "Transform" begin
        Random.seed!(SEED)
        @test isapprox(rotmatd(180, 0, 0), [1.0 0.0 0.0; 0.0 -1.0 0.0; 0.0 0.0 -1.0], rtol = RTOLERANCE, atol = ATOLERANCE)
        @test isapprox(rotmatd(0.0, 180.0, 0.0), [-1.0 0.0 0.0; 0.0 1.0 0.0; 0.0 0.0 -1.0], rtol = RTOLERANCE, atol = ATOLERANCE)
        @test isapprox(rotmatd(0, 0, 180), [-1.0 0.0 0.0; 0.0 -1.0 0.0; 0.0 0.0 1.0], rtol = RTOLERANCE, atol = ATOLERANCE)
        @test isapprox(rotmatd(0, 90, 0), [0.0 0.0 1.0; 0.0 1.0 0.0; -1.0 0.0 0.0], rtol = RTOLERANCE, atol = ATOLERANCE)
        @test isapprox(rotmatd(45, -45, 45), [0.5 -0.8535533905932737 0.1464466094067261; 0.5 0.14644660940672644 -0.8535533905932737; 0.7071067811865475 0.5 0.5], rtol = RTOLERANCE, atol = ATOLERANCE)
        x, y, z = rand(3)
        @test isapprox(rotmatd(x * 180 / π, y * 180 / π, z * 180 / π), rotmat(x, y, z), rtol = RTOLERANCE, atol = ATOLERANCE)

        @test isapprox(Transform(rotmatd(0, 90, 0), SVector(0.0, 0.0, 1.0)) * SVector(1.0, 0.0, 0.0), [0.0, 0.0, 0.0], rtol = RTOLERANCE, atol = ATOLERANCE)

        ta = Transform(rotmatd(0, 90, 0), SVector(0.0, 0.0, 1.0))
        tb = Transform(rotmatd(90, 0, 0), SVector(1.0, 0.0, 0.0))
        @test isapprox(collect(ta * tb), collect(Transform(rotmatd(90, 90, 0), SVector(0.0, 0.0, 0.0))), rtol = RTOLERANCE, atol = ATOLERANCE)

        @test isapprox(collect(ta * inv(ta)), collect(identitytransform()), rtol = RTOLERANCE, atol = ATOLERANCE)
        @test isapprox(collect(tb * inv(tb)), collect(identitytransform()), rtol = RTOLERANCE, atol = ATOLERANCE)
        @test isapprox(collect(inv(ta)), collect(Transform(rotmatd(0, -90, 0), SVector(1.0, 0.0, 0.0))), rtol = RTOLERANCE, atol = ATOLERANCE)

        #verify that matrix and vector{vector} forms of point representations are transformed the same way
        pts = SVector(SVector(1.0,1.0,1.0),SVector(2.0,2.0,2.0),SVector(3.0,3.0,3.0))
        ptsmat = SMatrix{3,3}(1.0,1.0,1.0,2.0,2.0,2.0,3.0,3.0,3.0)
        rottrans = rotationd(20,30,40)
        for i in 1:3
            @test isapprox((rottrans*pts[i]),(rottrans*ptsmat)[:,i])
        end
    end # testset Transform

    @testset "Interval" begin
        pt1 = Intersection(0.3, [4.0, 5.0, 6.0], normalize(rand(3)), 0.4, 0.5, NullInterface())
        pt2 = Intersection(0.5, [1.0, 2.0, 3.0], normalize(rand(3)), 0.2, 0.3, NullInterface())
        @test pt1 < pt2 && pt1 <= pt2 && pt2 > pt1 && pt2 >= pt1 && pt1 != pt2 && pt1 <= pt1
        intvl2 = positivehalfspace(pt1)
        @test α(halfspaceintersection(intvl2)) == α(pt1)

        ### interval intersection
        intersectionat = TestData.intersectionat
        ## interval/interval
        # one in another
        a = Interval(intersectionat(0.1), intersectionat(0.2))
        b = Interval(intersectionat(0.0), intersectionat(0.3))
        @test intervalintersection(a, b) == intervalintersection(b, a) == a
        # overlap
        a = Interval(intersectionat(0.1), intersectionat(0.3))
        b = Interval(intersectionat(0.2), intersectionat(0.4))
        @test intervalintersection(a, b) == intervalintersection(b, a) == Interval(intersectionat(0.2), intersectionat(0.3))
        # no overlap
        a = Interval(intersectionat(0.1), intersectionat(0.2))
        b = Interval(intersectionat(0.3), intersectionat(0.4))
        @test intervalintersection(a, b) == intervalintersection(b, a) isa EmptyInterval
        # start == end
        a = Interval(intersectionat(0.1), intersectionat(0.2))
        b = Interval(intersectionat(0.2), intersectionat(0.3))
        @test intervalintersection(a, b) == intervalintersection(b, a) isa EmptyInterval
        ## interval/du
        # encompass one
        a = Interval(intersectionat(0.0), intersectionat(0.3))
        b = DisjointUnion(Interval(intersectionat(0.1), intersectionat(0.2)), Interval(intersectionat(0.5), intersectionat(0.6)))
        @test intervalintersection(a, b) == intervalintersection(b, a) == b[1]
        # encompass two
        a = Interval(intersectionat(0.0), intersectionat(0.8))
        b = DisjointUnion(Interval(intersectionat(0.1), intersectionat(0.2)), Interval(intersectionat(0.5), intersectionat(0.6)))
        @test intervalintersection(a, b) == intervalintersection(b, a) == b
        # overlap one
        a = Interval(intersectionat(0.0), intersectionat(0.2))
        b = DisjointUnion(Interval(intersectionat(0.1), intersectionat(0.3)), Interval(intersectionat(0.5), intersectionat(0.6)))
        @test intervalintersection(a, b) == intervalintersection(b, a) == Interval(intersectionat(0.1), intersectionat(0.2))
        # overlap two
        a = Interval(intersectionat(0.2), intersectionat(0.5))
        b = DisjointUnion(Interval(intersectionat(0.1), intersectionat(0.3)), Interval(intersectionat(0.4), intersectionat(0.6)))
        res = intervalintersection(a, b)
        @test res == intervalintersection(b, a) && res[1] == Interval(intersectionat(0.2), intersectionat(0.3)) && res[2] == Interval(intersectionat(0.4), intersectionat(0.5))
        # no overlap
        a = Interval(intersectionat(0.1), intersectionat(0.2))
        b = DisjointUnion(Interval(intersectionat(0.3), intersectionat(0.4)), Interval(intersectionat(0.5), intersectionat(0.6)))
        @test intervalintersection(a, b) == intervalintersection(b, a) isa EmptyInterval
        ## du/du
        # no overlap
        a = DisjointUnion(Interval(intersectionat(0.0), intersectionat(0.1)), Interval(intersectionat(0.4), intersectionat(0.5)))
        b = DisjointUnion(Interval(intersectionat(0.2), intersectionat(0.3)), Interval(intersectionat(0.6), intersectionat(0.7)))
        @test intervalintersection(a, b) == intervalintersection(b, a) isa EmptyInterval
        # one in a overlaps one in b
        a = DisjointUnion(Interval(intersectionat(0.0), intersectionat(0.2)), Interval(intersectionat(0.4), intersectionat(0.5)))
        b = DisjointUnion(Interval(intersectionat(0.1), intersectionat(0.3)), Interval(intersectionat(0.6), intersectionat(0.7)))
        @test intervalintersection(a, b) == intervalintersection(b, a) == Interval(intersectionat(0.1), intersectionat(0.2))
        # one in a encompasses one in b
        a = DisjointUnion(Interval(intersectionat(0.0), intersectionat(0.3)), Interval(intersectionat(0.4), intersectionat(0.5)))
        b = DisjointUnion(Interval(intersectionat(0.1), intersectionat(0.2)), Interval(intersectionat(0.6), intersectionat(0.7)))
        @test intervalintersection(a, b) == intervalintersection(b, a) == b[1]
        # one in a overlaps two in b
        a = DisjointUnion(Interval(intersectionat(0.2), intersectionat(0.5)), Interval(intersectionat(0.7), intersectionat(0.8)))
        b = DisjointUnion(Interval(intersectionat(0.1), intersectionat(0.3)), Interval(intersectionat(0.4), intersectionat(0.6)))
        res = intervalintersection(a, b)
        @test res == intervalintersection(b, a) && res[1] == Interval(intersectionat(0.2), intersectionat(0.3)) && res[2] == Interval(intersectionat(0.4), intersectionat(0.5))
        # one in a encompasses two in b
        a = DisjointUnion(Interval(intersectionat(0.0), intersectionat(0.5)), Interval(intersectionat(0.7), intersectionat(0.8)))
        b = DisjointUnion(Interval(intersectionat(0.1), intersectionat(0.2)), Interval(intersectionat(0.3), intersectionat(0.4)))
        @test intervalintersection(a, b) == intervalintersection(b, a) == b
        # each in a overlap one in b
        a = DisjointUnion(Interval(intersectionat(0.0), intersectionat(0.2)), Interval(intersectionat(0.4), intersectionat(0.6)))
        b = DisjointUnion(Interval(intersectionat(0.1), intersectionat(0.3)), Interval(intersectionat(0.5), intersectionat(0.7)))
        res = intervalintersection(a, b)
        @test res == intervalintersection(b, a) && res[1] == Interval(intersectionat(0.1), intersectionat(0.2)) && res[2] == Interval(intersectionat(0.5), intersectionat(0.6))
        # each in a encompass one in b
        a = DisjointUnion(Interval(intersectionat(0.0), intersectionat(0.3)), Interval(intersectionat(0.4), intersectionat(0.7)))
        b = DisjointUnion(Interval(intersectionat(0.1), intersectionat(0.2)), Interval(intersectionat(0.5), intersectionat(0.6)))
        @test intervalintersection(a, b) == intervalintersection(b, a) == b
        ## empties
        intvl = positivehalfspace(pt1)
        du = intervalcomplement(Interval(pt1, pt2))
        @test intervalintersection(EmptyInterval(), EmptyInterval()) isa EmptyInterval
        @test intervalintersection(EmptyInterval(), intvl) isa EmptyInterval
        @test intervalintersection(intvl, EmptyInterval()) isa EmptyInterval
        @test intervalintersection(EmptyInterval(), du) isa EmptyInterval
        @test intervalintersection(du, EmptyInterval()) isa EmptyInterval

        ### interval union
        ## interval/interval
        # one in another
        a = Interval(intersectionat(0.1), intersectionat(0.2))
        b = Interval(intersectionat(0.0), intersectionat(0.3))
        @test intervalunion(a, b) == intervalunion(b, a) == b
        # overlap
        a = Interval(intersectionat(0.1), intersectionat(0.3))
        b = Interval(intersectionat(0.2), intersectionat(0.4))
        @test intervalunion(a, b) == intervalunion(b, a) == Interval(intersectionat(0.1), intersectionat(0.4))
        # no overlap
        a = Interval(intersectionat(0.1), intersectionat(0.2))
        b = Interval(intersectionat(0.3), intersectionat(0.4))
        res = intervalunion(a, b)
        @test res == intervalunion(b, a) && res[1] == a && res[2] == b
        # start == end
        a = Interval(intersectionat(0.1), intersectionat(0.2))
        b = Interval(intersectionat(0.2), intersectionat(0.3))
        @test intervalunion(a, b) == intervalunion(b, a) == Interval(intersectionat(0.1), intersectionat(0.3))
        ## interval/du
        # encompass one
        a = Interval(intersectionat(0.0), intersectionat(0.3))
        b = DisjointUnion(Interval(intersectionat(0.1), intersectionat(0.2)), Interval(intersectionat(0.5), intersectionat(0.6)))
        res = intervalunion(a, b)
        @test res == intervalunion(b, a) && res[1] == a && res[2] == b[2]
        # encompass two
        a = Interval(intersectionat(0.0), intersectionat(0.8))
        b = DisjointUnion(Interval(intersectionat(0.1), intersectionat(0.2)), Interval(intersectionat(0.5), intersectionat(0.6)))
        @test intervalunion(a, b) == intervalunion(b, a) == a
        # overlap one
        a = Interval(intersectionat(0.0), intersectionat(0.2))
        b = DisjointUnion(Interval(intersectionat(0.1), intersectionat(0.3)), Interval(intersectionat(0.5), intersectionat(0.6)))
        res = intervalunion(a, b)
        @test res == intervalunion(b, a) && res[1] == Interval(intersectionat(0.0), intersectionat(0.3)) && res[2] == b[2]
        # overlap two
        a = Interval(intersectionat(0.2), intersectionat(0.5))
        b = DisjointUnion(Interval(intersectionat(0.1), intersectionat(0.3)), Interval(intersectionat(0.4), intersectionat(0.6)))
        @test intervalunion(a, b) == intervalunion(b, a) == Interval(intersectionat(0.1), intersectionat(0.6))
        # no overlap
        a = Interval(intersectionat(0.1), intersectionat(0.2))
        b = DisjointUnion(Interval(intersectionat(0.3), intersectionat(0.4)), Interval(intersectionat(0.5), intersectionat(0.6)))
        res = intervalunion(a, b)
        @test res == intervalunion(b, a) && res[1] == a && res[2] == b[1] && res[3] == b[2]
        ## du/du
        # no overlap
        a = DisjointUnion(Interval(intersectionat(0.0), intersectionat(0.1)), Interval(intersectionat(0.4), intersectionat(0.5)))
        b = DisjointUnion(Interval(intersectionat(0.2), intersectionat(0.3)), Interval(intersectionat(0.6), intersectionat(0.7)))
        res = intervalunion(a, b)
        @test res == intervalunion(b, a) && res[1] == a[1] && res[2] == b[1] && res[3] == a[2] && res[4] == b[2]
        # one in a overlaps one in b
        a = DisjointUnion(Interval(intersectionat(0.0), intersectionat(0.2)), Interval(intersectionat(0.4), intersectionat(0.5)))
        b = DisjointUnion(Interval(intersectionat(0.1), intersectionat(0.3)), Interval(intersectionat(0.6), intersectionat(0.7)))
        res = intervalunion(a, b)
        @test res == intervalunion(b, a) && res[1] == Interval(intersectionat(0.0), intersectionat(0.3)) && res[2] == a[2] && res[3] == b[2]
        # one in a encompasses one in b
        a = DisjointUnion(Interval(intersectionat(0.0), intersectionat(0.3)), Interval(intersectionat(0.4), intersectionat(0.5)))
        b = DisjointUnion(Interval(intersectionat(0.1), intersectionat(0.2)), Interval(intersectionat(0.6), intersectionat(0.7)))
        res = intervalunion(a, b)
        @test res == intervalunion(b, a) && res[1] == a[1] && res[2] == a[2] && res[3] == b[2]
        # one in a overlaps two in b
        a = DisjointUnion(Interval(intersectionat(0.2), intersectionat(0.5)), Interval(intersectionat(0.7), intersectionat(0.8)))
        b = DisjointUnion(Interval(intersectionat(0.1), intersectionat(0.3)), Interval(intersectionat(0.4), intersectionat(0.6)))
        res = intervalunion(a, b)
        @test res == intervalunion(b, a) && res[1] == Interval(intersectionat(0.1), intersectionat(0.6)) && res[2] == a[2]
        # one in a encompasses two in b
        a = DisjointUnion(Interval(intersectionat(0.0), intersectionat(0.5)), Interval(intersectionat(0.7), intersectionat(0.8)))
        b = DisjointUnion(Interval(intersectionat(0.1), intersectionat(0.2)), Interval(intersectionat(0.3), intersectionat(0.4)))
        @test intervalunion(a, b) == intervalunion(b, a) == a
        # each in a overlap one in b
        a = DisjointUnion(Interval(intersectionat(0.0), intersectionat(0.2)), Interval(intersectionat(0.4), intersectionat(0.6)))
        b = DisjointUnion(Interval(intersectionat(0.1), intersectionat(0.3)), Interval(intersectionat(0.5), intersectionat(0.7)))
        res = intervalunion(a, b)
        @test res == intervalunion(b, a) && res[1] == Interval(intersectionat(0.0), intersectionat(0.3)) && res[2] == Interval(intersectionat(0.4), intersectionat(0.7))
        # each in a encompass one in b
        a = DisjointUnion(Interval(intersectionat(0.0), intersectionat(0.3)), Interval(intersectionat(0.4), intersectionat(0.7)))
        b = DisjointUnion(Interval(intersectionat(0.1), intersectionat(0.2)), Interval(intersectionat(0.5), intersectionat(0.6)))
        @test intervalunion(a, b) == intervalunion(b, a) == a
        ## empties
        @test intervalunion(EmptyInterval(), EmptyInterval()) isa EmptyInterval
        @test intervalunion(EmptyInterval(), intvl) == intvl
        @test intervalunion(intvl, EmptyInterval()) == intvl
        @test intervalunion(EmptyInterval(), du) == du
        @test intervalunion(du, EmptyInterval()) == du

        # interval complement
        int = intervalcomplement(EmptyInterval())
        @test lower(int) isa RayOrigin && upper(int) isa Infinity

        int = intervalcomplement(rayorigininterval(Infinity()))
        @test int isa EmptyInterval

        int = intervalcomplement(rayorigininterval(pt1))
        @test lower(int) == pt1 && upper(int) isa Infinity

        int = intervalcomplement(positivehalfspace(pt1))
        @test lower(int) isa RayOrigin && upper(int) == pt1

        int = intervalcomplement(Interval(pt1, pt2))
        @test lower(int[1]) isa RayOrigin && upper(int[1]) == pt1 && lower(int[2]) == pt2 && upper(int[2]) isa Infinity

        a = DisjointUnion(Interval(intersectionat(0.1), intersectionat(0.3)), Interval(intersectionat(0.4), intersectionat(0.7)))
        res = intervalcomplement(a)
        @test res[1] == rayorigininterval(intersectionat(0.1)) && res[2] == Interval(intersectionat(0.3), intersectionat(0.4)) && res[3] == positivehalfspace(intersectionat(0.7))

        a = DisjointUnion(rayorigininterval(intersectionat(0.2)), Interval(intersectionat(0.4), intersectionat(0.7)))
        res = intervalcomplement(a)
        @test res[1] == Interval(intersectionat(0.2), intersectionat(0.4)) && res[2] == positivehalfspace(intersectionat(0.7))
    end # testset interval
end # testset General


text:
using Test, Unicode, Pkg
using BinDeps

if VERSION >= v"0.7.0-DEV.3382"
    using Libdl
end

Pkg.build("Cairo")  # Tests apt-get code paths
using Cairo
Pkg.build("GSL")
using GSL

Pkg.build("Gumbo") # Test Autotools code paths
using Gumbo

# PR 171
@test BinDeps.lower(nothing, nothing) === nothing

# PR 271
BinDeps.debug("Cairo")

let gv = glibc_version()
    if Sys.islinux()
        lddv = lowercase(readchomp(`ldd --version`))
        if occursin("gnu", lddv) || occursin("glibc", lddv)
            @test isa(gv, VersionNumber)
            @test gv >= v"1.0.0"
        else
            # Assume non-glibc
            @test gv === nothing
        end
    else
        @test gv === nothing
    end
end


text:
module EMD

include("EMDUtil.jl")
using Splines

export findExtrema
export IMF
export IF

#Function to calculate intrinsic mode functions
#N is the maximum number of modes to find
function IMF(y, t, toldev=0.01, tolzero = 0.01, order=4, N=5,window=0)

	if window==0
		window = ones(length(t))
	end

	n = length(y)
	f = zeros(n,N)
	tempy = copy(y)
	eps = 0.00001

	n_modes = 0;

	for i = 1:N
	    avg = zeros(n,1)+1
	    sd = 2*toldev

	    while(mean(abs(avg))>tolzero && sd > toldev)

	    	# Interpolate a spline through the maxima and minima
	        max_ar, min_ar, tmax, tmin = findExtrema(tempy, t)

	        # # Don't use too high an order to interpolate, restrict it if there are not many extrema
	        # p_max = min(order, length(max_ar)-2)
	        # p_min = min(order, length(min_ar)-2)

	        # # Make even order
	        # p_max = Int(2*floor(p_max/2))
	        # p_min = Int(2*floor(p_min/2))

	        # # At least linear
	        # p_max = max(p_max,2)
	        # p_min = max(p_min,2)
	        p_max = (length(max_ar) >= order)? order : 4
	        p_min = (length(min_ar) >= order)? order : 4

	        S1 = Spline(max_ar, tmax, p_max)
	        S2 = Spline(min_ar, tmin, p_min)

	        # Find mean of envelope
	        avg = (S1(t) + S2(t)) / 2
	        avg = avg.*window

	        tempy = tempy-avg

	        sd = mean( (avg.^2)./((y-f[:,i]).^2 + eps) )

	        f[:,i] = f[:,i] + avg

	        #println(sd)
	    end

	    tempy = copy(f[:,i])

	    # Check to see if it's worth continuing or if the remainder is monotone
	    c = mean(abs(tempy))
	    d = diff(tempy)
	    n_modes = n_modes + 1
	    if all(d+c*tolzero .> 0) || all(d-c*tolzero .< 0)
	    	break
	    end
	end

	C = zeros(n,N)
	C[:,1] = y - f[:,1]
	for i = 2:n_modes
		C[:,i] = f[:,i-1]-f[:,i]
	end

	return C
end

#function to calculate instantaneous frequencies
function IF(C, t, x=0, order=4, window=0)

	n,N = size(C)
	n = length(t)

	if x == 0
		x = t+1e-6
	end

	if window == 0
		window = ones(length(x))
	end

	nx = length(x)
	Phi = zeros(nx,N)

	for i =1:N
		H1 = Spline(C[:,i], t, order)

		Phi[:,i] = 1/(2pi)*(H1(x).*H1(x,1,true) - H1(x, 1).*H1(x,0,true))./(H1(x).^2 + H1(x,0,true).^2)
		Phi[:,i] = Phi[:,i].*window
	end

	return Phi

end

end # module


text:
# deps.jl - Animated Particles

using Luxor, FileIO, ImageView, LinearAlgebra, Combinatorics

mutable struct particle
    center::Point
    radius::Float64
    velocity::Float64 # number of pixels to move per screen refresh
    angle::Float64
end

import Base.mod
function Base.:mod(a::Point, b::Array{Int64, 1})
    return Point(mod(a.x, b[1]), mod(a.y, b[2]))
end

import Base.+
function Base.:+(a::particle, b::particle)
    thisCenter = a.center + b.center
    thisRadius = a.radius + b.radius
    thisVelocity = a.velocity + b.velocity
    thisAng = a.angle + b.angle

    return particle(thisCenter, thisRadius, thisVelocity, thisAng)
end

# Helper functions
function genScreen(num_parts::Int, bounds::Array{Int64})
    # Returns a vector of particles
    displayed = Vector{particle}(undef, num_parts)
    scrnDiag = norm(bounds)
    for i = 1:num_parts
        thisCenter = randompoint(-div(bounds[1],2), -div(bounds[2], 2), div(bounds[1], 2), div(bounds[2], 2))
        thisRadius = rand(1:.1:sqrt(scrnDiag))
        thisVelocity = rand(0:.1:50)
        thisAng = rand(0:.05:2*π) # granularity of .05 radians
        displayed[i] = particle(thisCenter, thisRadius, thisVelocity, thisAng)
    end
    return displayed
end

function revScreen(partsList::Vector{particle})
    # Called once for each frame to be rendered
    numParts = length(partsList)
    lineList = []
    # Need to calculate any and all overlapping circles
    for connection in arraybyindex(partsList, combinations(1:numParts, 2))
        a = connection[1]
        b = connection[2]
        if intersection2circles(a.center, a.radius, b.center, b.radius) > 0
            # add the line between the centers
            append!(lineList, [a.center b.center])
        end
    end

    return lineList
end

function drawScreen(screen::Scene, framenumber)
    # Actually handle the Luxor interface to generate an image
    partsList = screen.opts
    bounds = convert.(Int, [screen.movie.width, screen.movie.height])
    # refresh particles
    revParticles!(partsList, div.(bounds, 2))
    # refresh screen
    currLines = revScreen(partsList)
    # draw
    sethue("white")
    setline(2.5)
    setdash("solid")
    for i in 1:2:length(currLines)
        line(currLines[i], currLines[i+1], :stroke)
    end

    sethue("blue")
    setline(.5)
    for i in 1:length(partsList)
        circle(partsList[i].center, partsList[i].radius, :stroke)
    end

end

function revParticle!(a::particle, bounds::Array{Int64})
# update single particle
# NOTE: Assumes bounds define the max x, y of a point.
    # updates position given velocity/angle
    # grows/shrinks circle randomly
    delCenter = a.velocity * Point(cos(a.angle), sin(a.angle))
    tmpCenter = a.center + delCenter
    # a.center = mod(tmpCenter, 2 .* bounds) - Point(bounds[1], bounds[2])
    a.center = tmpCenter
    a.radius = rand(.9:.01:1.1) * a.radius
    a.velocity = rand(.8:.01:1.1) * a.velocity
    a.angle = rand(.5:.01:1.5) * a.angle

    return a
end

function revParticles!(partsList::Vector, bounds::Array{Int64})
    # Changes Particles between screen refreshes
    for i = 1:length(partsList)
        partsList[i] = revParticle!(partsList[i], bounds)
    end

    return partsList
end

# Animation Functions
function backdrop(scene::Scene, framenumber)
    background("black")
end

function arraybyindex(a, idx)
    res = []
    for e in idx
        push!(res, a[e])
    end
    return res
end

max_stars_repo_name: mcbaron/procedural-art

max_issues_repo_name: mcbaron/procedural-art

max_forks_repo_name: mcbaron/procedural-art


text:
#=
    Specifies the rounding mode for the result of an approximate operation.
ARF_RND_NEAR   (4)
    Round to the nearest representable number, rounding to even if there is a tie.
ARF_RND_FLOOR  (2)
    Round to the nearest representable number in the direction towards minus infinity.
ARF_RND_CEIL   (3)
    Round to the nearest representable number in the direction towards plus infinity.
ARF_RND_DOWN   (0)
    Round to the nearest representable number in the direction towards zero.
ARF_RND_UP     (1)
    Round to the nearest representable number in the direction away from zero.
=#

const ArbRoundNearest  = Cint(4) # RoundingMode{:Nearest}()
const ArbRoundDown     = Cint(2) # RoundingMode{:Down}()
const ArbRoundUp       = Cint(3) # RoundingMode{:Up}()
const ArbRoundToZero   = Cint(0) # RoundingMode{:ToZero}()
const ArbRoundFromZero = Cint(1) # RoundingMode{:FromZero}()

const MpfrRoundNearest  = Cint(0) # RoundingMode{:Nearest}()
const MpfrRoundDown     = Cint(3) # RoundingMode{:Down}()
const MpfrRoundUp       = Cint(2) # RoundingMode{:Up}()
const MpfrRoundToZero   = Cint(1) # RoundingMode{:ToZero}()
const MpfrRoundFromZero = Cint(4) # RoundingMode{:FromZero}()

# matched with Julia's terminology

match_rounding_mode(::Type{Val{ArbRoundNearest}})      = RoundNearest
match_rounding_mode(::Type{Val{ArbRoundDown}})         = RoundDown
match_rounding_mode(::Type{Val{ArbRoundUp}})           = RoundUp
match_rounding_mode(::Type{Val{ArbRoundToZero}})       = RoundToZero
match_rounding_mode(::Type{Val{ArbRoundFromZero}})     = RoundFromZero

match_rounding_mode(::Type{Val{RoundNearest}})         = ArbRoundNearest
match_rounding_mode(::Type{Val{RoundDown}})            = ArbRoundDown
match_rounding_mode(::Type{Val{RoundUp}})              = ArbRoundUp
match_rounding_mode(::Type{Val{RoundToZero}})          = ArbRoundToZero
match_rounding_mode(::Type{Val{RoundFromZero}})        = ArbRoundFromZero
match_rounding_mode(::Type{Val{RoundNearestTiesAway}}) = ArbRoundNearest
match_rounding_mode(::Type{Val{RoundNearestTiesUp}})   = throw(ErrorException("RoundNearestTiesUp is not supported"))

@inline match_rounding_mode(mode::Cint) =
    match_rounding_mode(Val{mode})

@inline match_rounding_mode(mode::RoundingMode{S}) where {S} =
    match_rounding_mode(Val{mode})


text:
path = dirname(@__FILE__)

### Model
model = Model1002("ss10")

### Parameters

@testset "Test parameter bounds checking" begin
    for θ in model.parameters
        @test isa(θ, AbstractParameter)
        if !θ.fixed
            (left, right) = θ.valuebounds
            @test left < θ.value < right
        end
    end
end

### Model indices

@testset "Checking model indices" begin
    # Endogenous states
    endo = model.endogenous_states
    @test length(endo) == 68
    @test endo[:Ez_t] == 60

    # Exogenous shocks
    exo = model.exogenous_shocks
    @test length(exo) == 24
    @test exo[:corepce_sh] == 16

    # Expectation shocks
    ex = model.expected_shocks
    @test length(ex) == 13
    @test ex[:EL_f_sh] == 12

    # Equations
    eq = model.equilibrium_conditions
    @test length(eq) == 68
    @test eq[:eq_Ez] == 59

    # Additional states
    endo_new = model.endogenous_states_augmented
    @test length(endo_new) == 16
    @test endo_new[:y_t1] == 69

    # Observables
    obs = model.observables
    @test length(obs) == 19
    @test obs[:obs_tfp] == 12
end

### Equilibrium conditions
Γ0, Γ1, C, Ψ, Π = eqcond(model)

# Transition and measurement equations
TTT, RRR, CCC = solve(model)
meas = measurement(model, TTT, RRR, CCC)

### Pseudo-measurement equation
pseudo_meas = pseudo_measurement(model, TTT, RRR, CCC)

@testset "Check eqcond and state-space dimensions" begin
    # Matrices are of expected dimensions
    @test size(Γ0) == (68, 68)
    @test size(Γ1) == (68, 68)
    @test size(C) == (68,)
    @test size(Ψ) == (68, 24)
    @test size(Π) == (68, 13)

    @test size(meas[:ZZ]) == (19,84)
    @test size(meas[:DD]) == (19,)
    @test size(meas[:QQ]) == (24,24)
    @test size(meas[:EE]) == (19,19)

    @test size(TTT) == (84,84)
    @test size(RRR) == (84,24)
    @test size(CCC) == (84,)
end

G0_exp, G1_exp, C_exp, PSI_exp, PIE_exp = h5open("$path/eqcond.h5") do file
    read(file, "G0"), read(file, "G1"), read(file, "C"),
    read(file, "PSI"), read(file, "PIE")
end
TTT_exp, RRR_exp, CCC_exp = h5open("$path/transition.h5") do file
    read(file, "TTT"), read(file, "RRR"), read(file, "CCC")
end
Z_exp, D_exp, Q_exp, E_exp = h5open("$path/measurement.h5") do file
    read(file, "ZZ"), read(file, "DD"), read(file, "QQ"), read(file, "EE")
end
Z_pseudo_exp, D_pseudo_exp = h5open("$path/pseudo_measurement.h5") do file
    read(file, "ZZ_pseudo"), read(file, "DD_pseudo")
end

@testset "Compare eqcond, transition, measurement, and pseudo-measurement matrices against reference matrices" begin
    @test @test_matrix_approx_eq Γ0 G0_exp
    @test @test_matrix_approx_eq Γ1 G1_exp
    @test @test_matrix_approx_eq C C_exp
    @test @test_matrix_approx_eq Ψ PSI_exp
    @test @test_matrix_approx_eq Π PIE_exp

    @test @test_matrix_approx_eq TTT TTT_exp
    @test @test_matrix_approx_eq RRR RRR_exp
    @test @test_matrix_approx_eq CCC CCC_exp

    @test @test_matrix_approx_eq meas[:ZZ] Z_exp
    @test @test_matrix_approx_eq meas[:DD] D_exp
    @test @test_matrix_approx_eq meas[:QQ] Q_exp
    @test @test_matrix_approx_eq meas[:EE] E_exp

    @test @test_matrix_approx_eq pseudo_meas[:ZZ_pseudo] Z_pseudo_exp[1:21, :] # Test files had been generated w/extra pseudo-obs
    @test @test_matrix_approx_eq pseudo_meas[:DD_pseudo] D_pseudo_exp[1:21, :] # Test files had been generated w/extra pseudo-obs
end