text:
# -----------------------------------------------------------
# -   Model for the Dubins car with free final time example -
# -----------------------------------------------------------

using DifferentialEquations # not necessary is RungeKutta is used for shooting method (default)
using MINPACK
using NLsolve

using Distributions
using Random

export DubinsFree

# Model DubinsFree as a Julia class

mutable struct DubinsFree
    # State (x,y,theta,tf) and control (u) dimensions
    x_dim
    u_dim

    # Dynamics and linearized dynamics
    f
    A
    B

    # Problem settings
    dimLinearConstraintsU
    dimSecondOrderConeConstraintsU
    x_init
    x_final
    tf_guess
    tf
    xMin
    xMax
    v
    k

    B_free_final_time
    B_free_final_angle

    # Cylindrical obstacles (modeled by a center (x,y) and a radius r) and polygon obstacles (not used in this example)
    eps_obstacles
    obstacles
    poly_obstacles
    alpha

    # SCP parameters
    Delta0
    omega0
    omegamax
    # threshold for constraints satisfaction : constraints <= epsilon
    epsilon
    epsilon_xf_constraint
    convergence_threshold # in %
end

# The problem is set in the class constructor

function DubinsFree()
    x_dim = 4 # (last component is final time, only last index is used.)
    u_dim = 1

    v = 1.
    k = 2.
    B_free_final_time  = true # ALWAYS TRUE
    B_free_final_angle = false

    dimLinearConstraintsU = 0
    dimSecondOrderConeConstraintsU = 0
    x_init  = [ 0.; 0; 0]
    x_final = [ 2.; 3; 0]
    tf_guess = 3. # guess
    tf       = 1. # normalized time for dynamics constraints discretization
    myInf = 1.0e6 # Adopted to detect initial and final condition-free state variables
    xMin = [-100., -100., -2*pi]
    xMax = [ 100.,  100.,  2*pi]
    true_cost_weight = 1.

    Delta0 = 3.
    omega0 = 5000.
    omegamax = 1.0e6
    epsilon = 1e-3
    epsilon_xf_constraint = 0.
    convergence_threshold = 0.001

    # Obstacles
    eps_obstacles = 0.05
    obstacles = []
    obs = [[0.8,2],0.5]
    push!(obstacles, obs)
    alpha = 5000.

    poly_obstacles = []

    DubinsFree(x_dim, u_dim,
             [], [], [],
             dimLinearConstraintsU, dimSecondOrderConeConstraintsU,
             x_init, x_final,
             tf_guess, tf,
             xMin, xMax,
             v, k,
             B_free_final_time, B_free_final_angle,
             eps_obstacles,
             obstacles,
             poly_obstacles,
             alpha,
             Delta0,
             omega0,
             omegamax,
             epsilon,
             epsilon_xf_constraint,
             convergence_threshold)
end

# Method that returns the SCP parameters (used for set up)

function get_initial_scp_parameters(m::DubinsFree)
    return m.Delta0, m.omega0, m.omegamax, m.epsilon, m.convergence_threshold
end

# SCP is intialized by zero controls and velocities, and a straight-line in position

function initialize_trajectory(model::DubinsFree, N::Int)
  x_dim,  u_dim   = model.x_dim, model.u_dim
  x_init, x_final = model.x_init, model.x_final

  X        = zeros(x_dim, N)
  X[1:3,:] = hcat(range(x_init, stop=x_final, length=N)...)
  X[end,:] .= model.tf_guess

  U = zeros(u_dim, N-1)

  return X, U
end

# Method that returns the convergence ratio between iterations (in percentage)
# The quantities X, U denote the actual solution over time, whereas Xp, Up denote the solution at the previous step over time

function convergence_metric(model::DubinsFree, X, U, Xp, Up)
    x_dim = model.x_dim
    N = length(X[1,:])

    dt = model.tf/(N-1)
    Delta_uks = 0.
    for k in 1:(N-1)
        Delta_uks = Delta_uks + dt*(norm(U[:,k] - Up[:,k]))^2
    end
    return Delta_uks
end

# Method that returns the original cost

function true_NL_cost(model::DubinsFree, solver_model, X, U, Xp, Up)
    x_dim, u_dim     = model.x_dim, model.u_dim
    N = length(U[1,:])+1

    if model.B_free_final_time
        @NLexpression(solver_model, trueNLcost, (Xp[end,end]/(N-1)) * sum(U[i,k]^2 for k in 1:length(U[1,:]) for i in 1:u_dim))
    else
        @NLexpression(solver_model, trueNLcost, (model.tf/(N-1)) * sum(U[i,k]^2 for k in 1:length(U[1,:]) for i in 1:u_dim))
    end

    return trueNLcost
end

# The following methods return the i-th coordinate at the k-th iteration of the various constraints and their linearized versions (when needed)
# These are returned in the form " g(t,x(t),u(t)) <= 0 "

# Method that gathers all the linear control constraints

function control_linear_constraints(model::DubinsFree, X, U, Xp, Up, k, i)
    return 0.
end

function control_second_order_cone_constraints(model::DubinsFree, X, U, Xp, Up, k, i)
    throw(MethodError("[DubinsFree.jl::control_second_order_cone_constraints] No 2nd order constraints for this system."))
end

# State bounds and trust-region constraints (these are all convex constraints)

function state_max_convex_constraints(model::DubinsFree, X, U, Xp, Up, k, i)
    return ( X[i, k] - model.xMax[i] )
end

function state_min_convex_constraints(model::DubinsFree, X, U, Xp, Up, k, i)
    return ( model.xMin[i] - X[i, k] )
end

function trust_region_max_constraints(model::DubinsFree, X, U, Xp, Up, k, i, Delta)
    return ( (X[i, k] - Xp[i, k]) - Delta )
end

function trust_region_min_constraints(model::DubinsFree, X, U, Xp, Up, k, i, Delta)
    return ( -Delta - (X[i, k] - Xp[i, k]) )
end

# Method that checks whether trust-region constraints are satisifed or not (recall that trust-region constraints are penalized)

function is_in_trust_region(model::DubinsFree, X, U, Xp, Up, Delta)
    B_is_inside = true

    for k = 1:length(X[1,:])
        for i = 1:model.x_dim
            if trust_region_max_constraints(model, X, U, Xp, Up, k, i, Delta) > 0.
                B_is_inside = false
            end
            if trust_region_min_constraints(model, X, U, Xp, Up, k, i, Delta) > 0.
                B_is_inside = false
            end
        end
    end

    return B_is_inside
end

# Initial and final conditions on state variables

function state_initial_constraints(model::DubinsFree, X, U, Xp, Up)
    return ( X[1:3,1] - model.x_init )
end

function state_final_constraints(model::DubinsFree, X, U, Xp, Up)
    if model.B_free_final_angle
      return ( X[1:2,end] - model.x_final[1:2] )
    else
      return ( X[1:3,end] - model.x_final )
    end
end

# Methods that return the cylindrical obstacle-avoidance constraint and its lienarized version
# Here, a merely classical distance function is considered

function obstacle_constraint(model::DubinsFree, X, U, Xp, Up, k, obs_i,
                                               obs_type::String="sphere")
    #  obs_type    : Type of obstacles, can be 'sphere' or 'poly'
    if obs_type=="sphere"
      dim_obs           = 2
      p_obs, obs_radius = model.obstacles[obs_i][1], model.obstacles[obs_i][2]
      p_k               = X[1:dim_obs,k]

      dist       = norm(p_k - p_obs, 2)
      constraint = -( dist - obs_radius )
    else
      print("[DubinsFree.jl::obstacle_constraint_convexified] Unknown obstacle type.")
    end
    return constraint
end
function obstacle_constraint_convexified(model::DubinsFree, X, U, Xp, Up, k, obs_i,
                                                           obs_type::String="sphere")
    #  obs_type    : Type of obstacles, can be 'sphere' or 'poly'
    if obs_type=="sphere"
      dim_obs           = 2
      p_obs, obs_radius = model.obstacles[obs_i][1], model.obstacles[obs_i][2]
      p_k               = X[1:dim_obs,k]
      p_kp              = Xp[1:dim_obs,k]

      dist_prev  = norm(p_kp - p_obs, 2)
      dir_prev   = (p_kp - p_obs)/dist_prev
      constraint = -( dist_prev - obs_radius + sum(dir_prev[i] * (p_k[i] - p_kp[i]) for i=1:dim_obs) )
    else
      print("[DubinsFree.jl::obstacle_constraint_convexified] Unknown obstacle type.")
    end
    return constraint
end

# The following methods return the dynamical constraints and their linearized versions
# These are returned as time-discretized versions of the constraints " x' - f(x,u) = 0 " or " x' - A(t)*(x - xp) - B(t)*(u - up) = 0 ", respectively

# Method that returns the dynamical constraints and their linearized versions all at once

function compute_dynamics(model, Xp, Up)
    N = length(Xp[1,:])

    f_all, A_all, B_all = [], [], []

    for k in 1:N-1
        x_k = Xp[:,k]
        u_k = Up[:,k]

        f_dyn_k, A_dyn_k, B_dyn_k = f_dyn(x_k, u_k, model), A_dyn(x_k, u_k, model), B_dyn(x_k, u_k, model)

        push!(f_all, f_dyn_k)
        push!(A_all, A_dyn_k)
        push!(B_all, B_dyn_k)
    end

    return f_all, A_all, B_all
end

# These methods return the dynamics and its linearizations with respect to the state (matrix A(t)) and the control (matrix B(t)), respectively

function f_dyn(x::Vector, u::Vector, model::DubinsFree)
  x_dim = model.x_dim
  tf    = x[end]

  f = zeros(x_dim)
  f[1] = tf*model.v*cos(x[3])
  f[2] = tf*model.v*sin(x[3])
  f[3] = tf*model.k*u[1]
  f[4] = 0. # final time

  return f
end

function A_dyn(x::Vector, u::Vector, model::DubinsFree)
  x_dim = model.x_dim
  tf    = x[end]

  A      = zeros(x_dim,x_dim)
  A[1,3] = tf * (-model.v*sin(x[3]))
  A[2,3] = tf * ( model.v*cos(x[3]))

  A[1,4] = model.v*cos(x[3])
  A[2,4] = model.v*sin(x[3])
  A[3,4] = model.k*u[1]

  return A
end

function B_dyn(x::Vector, u::Vector, model::DubinsFree)
  x_dim, u_dim = model.x_dim, model.u_dim
  tf           = x[end]

  B      = zeros(x_dim,u_dim)
  B[:,1] = [0; 0; tf*model.k; 0]
  return B
end

function obstacle_potentialSphere_penalty(model::DubinsFree,
                                          pos, pos_obs, radius, alpha)
    pad = model.eps_obstacles
    if norm(pos-pos_obs)^2 < (radius+pad)^2
      return alpha * (norm(pos-pos_obs)^2-(radius+pad)^2)^2
    else
      return 0.
    end
end
function obstacle_potentialSphere_penalty_grad(model::DubinsFree,
                                               g,
                                               pos, pos_obs, radius, alpha)
    pad = model.eps_obstacles
    if norm(pos-pos_obs)^2 < (radius+pad)^2
      g[1] = alpha * 2. * (norm(pos-pos_obs)^2-(radius+pad)^2) * 2. * (pos-pos_obs)[1]
      g[2] = alpha * 2. * (norm(pos-pos_obs)^2-(radius+pad)^2) * 2. * (pos-pos_obs)[2]
    else
      g[1] = 0.
      g[2] = 0.
    end
    return g
end
function define_obs_potential_jump_NL_functions(model::DubinsFree, solver_model,
                                                obs_type::String="sphere")
    if obs_type != "sphere"
      throw(MethodError("[DubinsFree.jl::define_obs_potential_jump_NL_functions] Only spheres are supported."))
    end
    alpha = model.alpha

    # First obstacle
    pos_obs1, radius1 = model.obstacles[1][1], model.obstacles[1][2]
    obsPotFunc1     = (p1,p2)    -> obstacle_potentialSphere_penalty(     model,    [p1,p2], pos_obs1, radius1, alpha)
    obsPotFuncGrad1 = (g, p1,p2) -> obstacle_potentialSphere_penalty_grad(model, g, [p1,p2], pos_obs1, radius1, alpha)
    register(solver_model, :obsPotFunc1, 2, obsPotFunc1, obsPotFuncGrad1)

    # 2nd obstacle
    pos_obs2, radius2 = model.obstacles[2][1], model.obstacles[2][2]
    obsPotFunc2     = (p1,p2)    -> obstacle_potentialSphere_penalty(     model,    [p1,p2], pos_obs2, radius2, alpha)
    obsPotFuncGrad2 = (g, p1,p2) -> obstacle_potentialSphere_penalty_grad(model, g, [p1,p2], pos_obs2, radius2, alpha)
    register(solver_model, :obsPotFunc2, 2, obsPotFunc2, obsPotFuncGrad2)

    # 3-6 obstacles
    pos_obs3, radius3 = model.obstacles[3][1], model.obstacles[3][2]
    obsPotFunc3     = (p1,p2)    -> obstacle_potentialSphere_penalty(     model,    [p1,p2], pos_obs3, radius3, alpha)
    obsPotFuncGrad3 = (g, p1,p2) -> obstacle_potentialSphere_penalty_grad(model, g, [p1,p2], pos_obs3, radius3, alpha)
    register(solver_model, :obsPotFunc3, 2, obsPotFunc3, obsPotFuncGrad3)
    pos_obs4, radius4 = model.obstacles[4][1], model.obstacles[4][2]
    obsPotFunc4     = (p1,p2)    -> obstacle_potentialSphere_penalty(     model,    [p1,p2], pos_obs4, radius4, alpha)
    obsPotFuncGrad4 = (g, p1,p2) -> obstacle_potentialSphere_penalty_grad(model, g, [p1,p2], pos_obs4, radius4, alpha)
    register(solver_model, :obsPotFunc4, 2, obsPotFunc4, obsPotFuncGrad4)
    pos_obs5, radius5 = model.obstacles[5][1], model.obstacles[5][2]
    obsPotFunc5     = (p1,p2)    -> obstacle_potentialSphere_penalty(     model,    [p1,p2], pos_obs5, radius5, alpha)
    obsPotFuncGrad5 = (g, p1,p2) -> obstacle_potentialSphere_penalty_grad(model, g, [p1,p2], pos_obs5, radius5, alpha)
    register(solver_model, :obsPotFunc5, 2, obsPotFunc5, obsPotFuncGrad5)
    pos_obs6, radius6 = model.obstacles[6][1], model.obstacles[6][2]
    obsPotFunc6     = (p1,p2)    -> obstacle_potentialSphere_penalty(     model,    [p1,p2], pos_obs6, radius6, alpha)
    obsPotFuncGrad6 = (g, p1,p2) -> obstacle_potentialSphere_penalty_grad(model, g, [p1,p2], pos_obs6, radius6, alpha)
    register(solver_model, :obsPotFunc6, 2, obsPotFunc6, obsPotFuncGrad6)
end
function obstacle_potential_penalties(model::DubinsFree, solver_model,
                                    X, U, Xp, Up,
                                    obs_i,
                                    obs_type::String="sphere")
    N, Nb_obstacles = length(X[1,:]), length(model.obstacles)
    if Nb_obstacles != 1 && Nb_obstacles != 2 && Nb_obstacles != 6
      throw(MethodError("[DubinsFree.jl::obstacle_potential_penalties] Incorrect number of obstacles."))
    end
    if !(obs_type=="sphere")
      print("[DubinsFree.jl::obstacle_potential_penalties] Unknown obstacle type.")
    end

    if !(model.B_free_final_time)
      print("[DubinsFree.jl::obstacle_potential_penalties] tf should be free.")
    end

    if obs_i == 1
      @NLexpression(solver_model, obscost1, (X[end,end]/(N-1))*sum(obsPotFunc1(X[1,k],X[2,k]) for k in 1:N))
      return obscost1
    elseif obs_i == 2
      @NLexpression(solver_model, obscost2, (X[end,end]/(N-1))*sum(obsPotFunc2(X[1,k],X[2,k]) for k in 1:N))
      return obscost2
    elseif obs_i == 3
      @NLexpression(solver_model, obscost3, (X[end,end]/(N-1))*sum(obsPotFunc3(X[1,k],X[2,k]) for k in 1:N))
      return obscost3
    elseif obs_i == 4
      @NLexpression(solver_model, obscost4, (X[end,end]/(N-1))*sum(obsPotFunc4(X[1,k],X[2,k]) for k in 1:N))
      return obscost4
    elseif obs_i == 5
      @NLexpression(solver_model, obscost5, (X[end,end]/(N-1))*sum(obsPotFunc5(X[1,k],X[2,k]) for k in 1:N))
      return obscost5
    elseif obs_i == 6
      @NLexpression(solver_model, obscost6, (X[end,end]/(N-1))*sum(obsPotFunc6(X[1,k],X[2,k]) for k in 1:N))
      return obscost6
    else
      print("obs_i=$obs_i")
      throw(MethodError("[DubinsFree.jl::obstacle_potential_penalty] Incorrect number of obstacles."))
    end
end

##### ----------------------------- #####
##### FUNCTIONS FOR SHOOTING METHOD #####

function shooting_solve(p0::Vector, tf_guess::Float64,
                        N::Int, model::DubinsFree;
                        tol_newton=1e-2, tol_hybrd = 1e-6,
                        B_RK4::Bool=false, B_HYBRD::Bool=true, B_verbose=true)
  x_dim, u_dim = model.x_dim, model.u_dim
  x0, xF       = model.x_init, model.x_final

  # Set up shooting function
  if !B_RK4
    shooting_eval! = (F, p_0_tf) -> parameterized_shooting_eval!(F,
                                                            x0, xF, p_0_tf[1:3], p_0_tf[4],
                                                            model)
  else
    N_RK4shooting = N
    shooting_eval! = (F, p_0_tf) -> parameterized_shooting_eval_RK4!(F,
                                                            x0, xF, p_0_tf[1:3], p_0_tf[4],
                                                            model, N_RK4shooting)
  end

  if !B_HYBRD
    # Run Newton method
    sol_newton = nlsolve(shooting_eval!, [p0;tf_guess], iterations=20, ftol=tol_newton)
    # sol_newton = nlsolve(shooting_eval!, [p0;tf_guess], ftol=tol_newton)
    p0_tf_solution = sol_newton.zero
    if B_verbose
      println(sol_newton)
      println(p0_tf_solution)
    end
  else
    # run Advanced modified version of Powell's algorithm
    # sol_minpack = fsolve(shooting_eval!, [p0;tf_guess], iterations=100, tol=tol_hybrd, method=:hybrd)
    sol_minpack = fsolve(shooting_eval!, [p0;tf_guess], tol=tol_hybrd, method=:hybrd)
    p0_tf_solution = sol_minpack.x
    if B_verbose
      println(sol_minpack)
      println("p0_tf=",p0_tf_solution)
    end
  end

  if (!B_HYBRD && sol_newton.f_converged) || (B_HYBRD && sol_minpack.converged)
    xp_0  = [x0; p0_tf_solution[1:3]; p0_tf_solution[4]]

    # check tolerance (again)
    F = zeros(x_dim)
    shooting_eval!(F, p0_tf_solution)
    if (!B_HYBRD && (norm(sum(F)) > 10*x_dim*tol_newton)) ||
       (B_HYBRD  && (norm(sum(F)) > 10*x_dim*tol_hybrd))
      if B_verbose
        println("Shooting method DID NOT converge ! (F=$(F)))")
      end
      return [],[],[],-1,false
    end
    # check that time did not jump too much
    if norm(p0_tf_solution[4]-tf_guess)>3.
      if B_verbose
        println("Shooting method didnt converge - delta tf > 3.")
      end
      return [],[],[],-1,false
    end

    # Solve again to get trajectory
    if !B_RK4
      tspan = (0., p0_tf_solution[4])

      shooting_ode! = (f, u, p, t) -> dynamics_shooting!(f, u[1:3], u[4:6], u[7], model)
      prob = ODEProblem(shooting_ode!, xp_0, tspan)

      dt  = p0_tf_solution[4] / (N-1)

      sol = DifferentialEquations.solve(prob, saveat=dt)#, dtmin=0.1, force_dtmin=true)
      xp   = hcat(sol.u...)

      Xtraj, Ptraj, tf = xp[1:3,:], xp[4:6,:], xp[end,end]
    else
      shooting_ode! = (f, u, t) -> dynamics_shooting!(f, u[1:3], u[4:6], u[7], model)
      t_vec, traj_u = solveRK4(shooting_ode!, model, N, p0_tf_solution[4], xp_0)
      Xtraj, Ptraj, tf = traj_u[1:3,:], traj_u[4:6,:], t_vec[end]
    end
    Xtraj = vcat(Xtraj,tf*ones(1,N))

    Utraj = zeros(u_dim, N-1)
    for k in 1:(N-1)
      Utraj[1,k] = ( Ptraj[3,k]*model.k / 2.)
    end
    return Xtraj, Utraj, Ptraj, tf, true
  else
    if B_verbose
      println("Shooting method DID NOT converge !")
    end
    return [],[],[],-1,false
  end

end

function parameterized_shooting_eval!(F::Vector,
                                      x0::Vector, xF::Vector, p0::Vector, tf_guess::Float64,
                                      model::DubinsFree)
  v, k = model.v, model.k

  xp_0  = [x0; p0; tf_guess]
  tspan = (0., tf_guess)

  shooting_ode! = (f, u, p, t) -> dynamics_shooting!(f, u[1:3], u[4:6], u[7], model)
  prob = ODEProblem(shooting_ode!, xp_0, tspan)
  sol  = DifferentialEquations.solve(prob)

  for i = 1:model.x_dim
    if (i == 1) || (i == 2)
      F[i] = (sol.u[end][i]-xF[i])
    elseif i == 3
      if model.B_free_final_angle
        F[i] = sol.u[end][6] # p_theta(t_f) = 0.
      else
        F[i] = sqrt((cos(xF[i])-cos(sol.u[end][i]))^2 + (sin(xF[i])-sin(sol.u[end][i]))^2)
      end
    elseif i == 4
      tf = tf_guess
      xf = sol.u[end][1:3]
      pf = sol.u[end][4:6]
      uf = (pf[3]*k/2.)
      ff = [v * cos(xf[3]), v * sin(xf[3]), k * uf] #f_dyn(xf, [uf], model)
      Hf = sum(pf[i]*ff[i] for i in 1:3) - uf^2

      # Obstacles
      for obs_i in 1:length(model.obstacles)
        pos_obs, radius = model.obstacles[obs_i][1], model.obstacles[obs_i][2]
        Hf = Hf - obstacle_potentialSphere_penalty(model, [traj_end[1], traj_end[2]],
                                                          pos_obs, radius, model.alpha)
      end
      F[i] = Hf
    else
      println("DubinsFree doesnt have so many dimensions!!!")
    end
  end
end

function parameterized_shooting_eval_RK4!(F::Vector,
                                      x0::Vector, xF::Vector, p0::Vector, tf_guess::Float64,
                                      model::DubinsFree, N::Int)
  v, k = model.v, model.k

  xp_0  = [x0; p0; tf_guess]
  tspan = (0., tf_guess)

  shooting_ode! = (f, u, t) -> dynamics_shooting!(f, u[1:3], u[4:6], u[7], model)
  t_vec, traj_u = solveRK4(shooting_ode!, model, N, tf_guess, xp_0)
  # print("traj_u=", size(traj_u))
  traj_end = traj_u[:,end]

  for i = 1:model.x_dim
    if (i == 1) || (i == 2)
      F[i] = traj_end[i]-xF[i]
    elseif i == 3
      if model.B_free_final_angle
        F[i] = traj_end[6] # p_theta(t_f) = 0.
      else
        F[i] = (cos(xF[i])-cos(traj_end[i]))^2 + (sin(xF[i])-sin(traj_end[i]))^2
      end
    elseif i == 4
      tf = tf_guess
      xf = traj_end[1:3]
      pf = traj_end[4:6]
      uf = (pf[3]*k/2.)
      ff = [v * cos(xf[3]), v * sin(xf[3]), k * uf] #f_dyn(xf, [uf], model)
      Hf = sum(pf[i]*ff[i] for i in 1:3) - uf^2

      # Obstacles
      for obs_i in 1:length(model.obstacles)
        pos_obs, radius = model.obstacles[obs_i][1], model.obstacles[obs_i][2]
        Hf = Hf - obstacle_potentialSphere_penalty(model, [traj_end[1], traj_end[2]],
                                                          pos_obs, radius, model.alpha)
      end

      F[i] = Hf
    else
      println("DubinsFree doesnt have so many dimensions!!!")
    end
  end
end

function solveRK4(dynamics::Function,model::DubinsFree, N, tf, z0, t0=0.0)
    # dynamics : (f,x,t) -> populates f=f(x) (= x_dot)

    x_dim = length(z0)

    h = (tf - t0)/(N - 1)
    t = zeros(N)
    z = zeros(x_dim,N)
    dynTemp = zeros(x_dim)
    zTemp  = zeros(x_dim)
    t[1]   = t0
    z[:,1] = z0

    for i = 1:N-1
        t[i+1] = t[i] + h

        dynamics(dynTemp,z[:,i],t[i])
        z[:,i+1] = z[:,i] + h*dynTemp/6.0

        zTemp = z[:,i] + h*dynTemp/2.0
        dynamics(dynTemp,zTemp,t[i]+h/2.0)
        z[:,i+1] = z[:,i+1] + h*dynTemp/3.0

        zTemp = z[:,i] + h*dynTemp/2.0
        dynamics(dynTemp,zTemp,t[i]+h/2.0)
        z[:,i+1] = z[:,i+1] + h*dynTemp/3.0

        zTemp = z[:,i] + h*dynTemp
        dynamics(dynTemp,zTemp,t[i]+h)
        z[:,i+1] = z[:,i+1] + h*dynTemp/6.0
    end

    return t, z
end

function dynamics_shooting!(f::Vector, x::Vector, p::Vector, tf::Float64,
                            model::DubinsFree)
    v, k = model.v, model.k
    uopt = ( p[3]*k / 2.)

    f[1] = v * cos(x[3])
    f[2] = v * sin(x[3])
    f[3] = k * uopt
    f[4] = 0.
    f[5] = 0.
    f[6] = v * (sin(x[3])*p[1] - cos(x[3])*p[2])
    f[7] = 0.

    # Obstacles
    for obs_i in 1:length(model.obstacles)
      g               = zeros(2)
      pos_obs, radius = model.obstacles[obs_i][1], model.obstacles[obs_i][2]
      grads_obs       = obstacle_potentialSphere_penalty_grad(model, g, [x[1], x[2]],
                                                              pos_obs, radius, model.alpha)
      f[4]            = f[4] + grads_obs[1]
      f[5]            = f[5] + grads_obs[2]
    end
end

# randomized Monte-Carlo experiments
function get_randomized_conditions(model::DubinsFree; N_MC::Int=10, NB_obs::Int=1, seed=12345)
  rng = MersenneTwister(seed)

  # Boundary Conditions
  min_pos, max_pos = -1., 1.
  min_θ0f, max_θ0f = -pi, pi

  x0s = [transpose(min_pos .+ (max_pos-min_pos) .* rand(rng, N_MC));
         transpose(min_pos .+ (max_pos-min_pos) .* rand(rng, N_MC));
         transpose(min_θ0f .+ (max_θ0f-min_θ0f) .* rand(rng, N_MC))]
  x0s = transpose(x0s) # (N_MC x x_dim)

  # Sample final states so they are far enough from each other
  angles    = -pi .+ 2*pi * rand(rng, N_MC)
  dist_x0fs = 4.  .+ 3    * rand(rng, N_MC)
  xfs = [transpose(dist_x0fs .* cos.(angles));
        transpose(dist_x0fs .* sin.(angles));
        transpose(min_θ0f .+ (max_θ0f-min_θ0f) .* rand(rng, N_MC))]
  xfs = transpose(xfs) + x0s

  # set final and initial angles within a cone
  delta_x0fs = xfs-x0s
  for i = 1:N_MC
    theta = atan(delta_x0fs[i,2], delta_x0fs[i,1])
    min_θ0f, max_θ0f = -pi/4+theta, pi/4+theta
    x0s[i,3] = min_θ0f + (max_θ0f-min_θ0f) .* rand(rng)
    xfs[i,3] = min_θ0f + (max_θ0f-min_θ0f) .* rand(rng)
  end

  # Obstacles
  rads = 0.4*ones(N_MC)

  obstacles_vec = []
  for i=1:N_MC
    obss = []
    for j=1:NB_obs
      r    = rads[i]

      # place obs randomly in between
      min_xb, max_xb = min(x0s[i,1], xfs[i,1])+2*r, max(x0s[i,1], xfs[i,1])-2*r
      min_yb, max_yb = min(x0s[i,2], xfs[i,2])+2*r, max(x0s[i,2], xfs[i,2])-2*r
      if min_xb>max_xb
          min_xb_old, max_xb_old = min_xb, max_xb
          min_xb = (min_xb_old+max_xb_old)/2
          max_xb = (min_xb_old+max_xb_old)/2
      end
      if min_yb>max_yb
          min_yb_old, max_yb_old = min_yb, max_yb
          min_yb = (min_yb_old+max_yb_old)/2
          max_yb = (min_yb_old+max_yb_old)/2
      end
      o_x = min_xb+(max_xb-min_xb).*rand(rng)
      o_y = min_yb+(max_yb-min_yb).*rand(rng)
      o   = [o_x,o_y]

      # Check that obstacle is some distance appart from other obstacles
      if j>1
        op = obss[end][1]
        if norm(o-op)<5*r
          o = op + (o-op)/norm(o-op) * 5*r
          o_x,o_y = o[1], o[2]
        end
      end

      # push obstacles so they are 8r from x0 and xf
      if norm(o-x0s[i,1:2]) < 8*r
        o = x0s[i,1:2] + (o-x0s[i,1:2])/norm(o-x0s[i,1:2]) * 8*r
        o_x,o_y = o[1], o[2]
      end
      if norm(o-xfs[i,1:2]) < 8*r
        o = xfs[i,1:2] + (o-xfs[i,1:2])/norm(o-xfs[i,1:2]) * 8*r
        o_x,o_y = o[1], o[2]
      end

      # Again: Check that obstacle is some distance appart from other obstacles
      if j>1
        op = obss[end][1]
        if norm(o-op)<3*r
          o = op + (o-op)/norm(o-op) * 3*r
          o_x,o_y = o[1], o[2]
        end
      end

      # save a single obstacle
      obs = [[o_x, o_y], rads[i]]
      push!(obss, obs)
    end
    push!(obstacles_vec, obss)
  end

  # final times
  # min_tf, max_tf = 0.5, 2.
  min_tf, max_tf = 4, 6.
  tfs = min_tf .+ (max_tf-min_tf) .* rand(rng, N_MC)

  return x0s, xfs, obstacles_vec, tfs
end

max_stars_repo_name: StanfordASL/jlSCP

max_issues_repo_name: StanfordASL/jlSCP

max_forks_repo_name: StanfordASL/jlSCP


text:
#################### Random Walk Metropolis ####################

#################### Types and Constructors ####################

mutable struct RWMTune <: SamplerTune
  logf::Union{Function, Missing}
  scale::Union{Float64, Vector{Float64}}
  proposal::SymDistributionType

  RWMTune() = new()

  function RWMTune(x::Vector, scale::Real, logf::Union{Function, Missing};
                   proposal::SymDistributionType=Normal)
    new(logf, Float64(scale), proposal)
  end

  function RWMTune(x::Vector, scale::Vector{T},
                  logf::Union{Function, Missing};
                  proposal::SymDistributionType=Normal) where {T<:Real}
    new(logf, convert(Vector{Float64}, scale), proposal)
  end
end

RWMTune(x::Vector, scale::ElementOrVector{T}; args...) where {T<:Real} =
  RWMTune(x, scale, missing; args...)

const RWMVariate = SamplerVariate{RWMTune}

validate(v::RWMVariate) = validate(v, v.tune.scale)

validate(v::RWMVariate, scale::Float64) = v

function validate(v::RWMVariate, scale::Vector)
  n = length(v)
  length(scale) == n ||
    throw(ArgumentError("length(scale) differs from variate length $n"))
  v
end

#################### Sampler Constructor ####################

function RWM(params::ElementOrVector{Symbol},
              scale::ElementOrVector{T}; args...) where {T<:Real}
  samplerfx = function(model::Model, block::Integer)
    block = SamplingBlock(model, block, true)
    v = SamplerVariate(block, scale; args...)
    sample!(v, x -> logpdf!(block, x))
    relist(block, v)
  end
  Sampler(params, samplerfx, RWMTune())
end

#################### Sampling Functions ####################

sample!(v::RWMVariate) = sample!(v, v.tune.logf)

function sample!(v::RWMVariate, logf::Function)
  x = v + v.tune.scale .* rand(v.tune.proposal(0.0, 1.0), length(v))
  if rand() < exp(logf(x) - logf(v.value))
    v[:] = x
  end
  v
end

ext: jl


text:
# Main SimString search algorithm

"""
    search(measure::AbstractSimilarityMeasure, db_collection::AbstractSimStringDB, query::AbstractString;
        α=0.7, ranked=true)

Search for strings in a string collection using the SimString algorithm and a
similarity measure.

# Arguments:
* `measure`::AbstractSimilarityMeasure - The similarity measure to use.
* `db_collection`::AbstractSimStringDB - The database collection to search.
* `query`::AbstractString - The query string to search for.
* `α`::float - The α parameter for the SimString algorithm.
* `ranked`::Boolean - Whether to return the results in ranked order.

# Example
```julia
db = DictDB(CharacterNGrams(2, " "));
append!(db, ["foo", "bar", "fooo"]);

search(Dice(), db, "foo"; α=0.8, ranked=true)
# 2-element Vector{Tuple{String, Float64}}:
#  ("foo", 1.0)
#  ("fooo", 0.8888888888888888)
```

# Returns
* A Vector of results, where each element is a Tuple of the form (`string`, `similarity measure score`).
"""
function search(measure::AbstractSimilarityMeasure, db_collection::AbstractSimStringDB, query::AbstractString; α=0.7, ranked=true)
    return search!(measure, db_collection, query; α=α, ranked=ranked)
end

"""
Internal function which ranks the results of a search using the specified similarity measure.
"""
function rank_search_results(measure::AbstractSimilarityMeasure, db_collection::DictDB, query, results; ranked=true)
    features = extract_features(db_collection.feature_extractor, query)

    # Compute similarity scores for each result
    ranked_results = map(results) do x
        x, similarity_score(measure, features, extract_features(db_collection.feature_extractor, x) )
    end

    # Sort by similarity score and return
    return ranked ? sort(ranked_results, by = i -> i[2], rev=true) : ranked_results
end

"""
Internal function which performs the overlap join
"""
function overlap_join(db_collection::AbstractSimStringDB, features, τ, candidate_size)
    # length of features
    query_feature_length = length(features)

    # Sort features from the most uncommon and the most common
    features = sort(features, by = i -> length(lookup_feature_set_by_size_feature(db_collection, candidate_size, i) ) )

    # Count the occurrences of each feature
    candidate_match_counts = DefaultDict{String, Int}(0)
    feature_slice_index = query_feature_length - τ + 1
    idx = query_feature_length - τ
    focus_features = feature_slice_index < 0 ? (@view features[0:end + feature_slice_index]) : (@view features[0:idx])

    @inbounds @views for i in focus_features
        for s in lookup_feature_set_by_size_feature(db_collection, candidate_size, i)
            candidate_match_counts[s] += 1
        end
    end

    results = String[]

    for (candidate, match_count) in candidate_match_counts
        for i in (query_feature_length - τ + 1) : query_feature_length
            if candidate in lookup_feature_set_by_size_feature(db_collection, candidate_size, features[i])
                match_count += 1
            end

            if match_count >= τ
                append!(results, [candidate])
                break
            end

            remaining_count = query_feature_length - i - 1
            if (match_count + remaining_count) < τ
                break
            end
        end
    end
    return results
end

"""
Search for strings in custom DictDB string collection using the SimString algorithm
and a similarity measure.
"""
function search!(measure::AbstractSimilarityMeasure, db_collection::DictDB, query::AbstractString; α=0.7, ranked=true)
    # Generate features from query string
    features = extract_features(db_collection.feature_extractor, query)

    results = String[]

    # Generate and return results from the potential candidate size pool
    @inbounds for candidate_size in minimum_feature_size(measure, length(features), α) : maximum_feature_size(measure, db_collection, length(features), α)
        # Minimum overlap
        τ = minimum_overlap(measure, length(features), candidate_size, α)

        # Generate approximate candidates from the overlap join
        append!(results, overlap_join(db_collection, features, τ, candidate_size))
    end

    # Rank search results
    return rank_search_results(measure, db_collection, query, results; ranked=ranked)
end


text:
###########################################################################
export VecScatter

# describes a scatter operation context (input/output index sets etc.)
 mutable struct VecScatter{T}
  p::C.VecScatter{T}

  function VecScatter{T}(p::C.VecScatter{T}) where {T}
    o = new{T}(p)
    if (!deactivate_finalizers)
      finalizer(VecScatterDestroy, o)
    end
    return o
  end
end

comm(a::VecScatter{T}) where {T} = MPI.Comm(C.PetscObjectComm(T, a.p.pobj))

function VecScatterDestroy(o::VecScatter{T}) where {T}
  PetscFinalized(T) || C.VecScatterDestroy(Ref(o.p))
end

function VecScatter(x::Vec{T}, ix::IS{T}, y::Vec{T}, iy::IS{T}) where {T}
  scatter_c = Ref{C.VecScatter{T}}()
  chk(C.VecScatterCreate(x.p, ix.p, y.p, iy.p, scatter_c))
  return VecScatter{T}(scatter_c[])
end

function Base.copy(i::VecScatter{T}) where {T}
  vs_c = Ref{C.VecScatter{T}}()
  chk(C.VecScatterCopy(i.p, vs_c))
  return VecScatter{T}(vs_c[])
end

###########################################################################
export scatter!

function scatter!(scatter::VecScatter{T}, x::Vec{T}, y::Vec{T}; imode=C.INSERT_VALUES, smode=C.SCATTER_FORWARD) where {T}
  chk(C.VecScatterBegin(scatter.p, x.p, y.p, imode, smode))
  yield() # do async computations while messages are in transit
  chk(C.VecScatterEnd(scatter.p, x.p, y.p, imode, smode))
  return y
end

function scatter!(x::Vec{T}, ix::AbstractVector{I1},
                           y::Vec{T}, iy::AbstractVector{I2};
                          imode=C.INSERT_VALUES, smode=C.SCATTER_FORWARD) where {T,I1,I2}
  scatter = VecScatter(x, IS(T, ix, comm=comm(x)),
                       y, IS(T, iy, comm=comm(y)))
  scatter!(scatter, x, y; imode=imode, smode=smode)
end


text:
module WhitenedRBMs

import Random
import Statistics
import LinearAlgebra
import ValueHistories
import Flux
import RestrictedBoltzmannMachines as RBMs

using LinearAlgebra: Diagonal, cholesky, diagm, Symmetric, I
using ValueHistories: MVHistory
using RestrictedBoltzmannMachines: RBM, AbstractLayer
using RestrictedBoltzmannMachines: Binary, Spin, Potts, Gaussian, ReLU, dReLU, pReLU, xReLU
using RestrictedBoltzmannMachines: visible, hidden, weights, flatten
using RestrictedBoltzmannMachines: inputs_v_to_h, inputs_h_to_v
using RestrictedBoltzmannMachines: batchmean, batchvar

include("affine.jl")
include("whiterbm.jl")
include("whiten.jl")
include("data.jl")
include("binary_white_rbm.jl")
include("layers.jl")
include("train/pcd.jl")

end # module


text:
module util_convexopt

using LinearAlgebra
export LinOp, compute_opnorm, proj_l1!, grad!, div!, prox_nuclear!

# See http://www.ima.umn.edu/materials/2019-2020/SW10.14-18.19/28302/talk.pdf
# struct LinOp
#     forward::Function
#     adjoint::Function
# end
# Base.:(*)(A::LinOp, x) = A.forward(x) # for A*x
# Base.adjoint(A::LinOp) = LinOp(A.adjoint, A.forward) # for A'

# D_  = f -> grad(f)
# Dt_ = u -> -div(u)

# "D: gradient operator"
# D = LinOp(D_, Dt_)

# export D

"Compute opeartor norm based on the Power method"
function compute_opnorm(A, niter=3)
    lamb = 0
    x = rand(size(A,2))
    x_next = similar(x)
    Ax = zeros(size(A, 1))
    At = A'

    for i=1:niter
        mul!(Ax, A, x) # A'*A*x is very slow (WHY?) Todo
        mul!(x_next, At, Ax)
        lamb_prev = lamb
        lamb = sqrt(sum(x_next.^2))

        if (abs(lamb_prev-lamb) < 1e-9)
            println("break at iteration: ", i)
            break
        end
        x = x_next / lamb
    end
    #eig_max = x'*A'*A*x / x*x
    return sqrt(lamb)
end

"Compute forward gradient [HxWx2] of an image `u` with Neumann boundary condition."
function grad(u::Array{T, 2}) where {T<:AbstractFloat}
    ux = circshift(u, [0, -1]) - u
    uy = circshift(u, [-1, 0]) - u
    ux[:, end] .= 0.0
    uy[end, :] .= 0.0

    return cat(ux, uy, dims=3)
end

function grad!(du::Array{T, 3}, u::Array{T, 2}) where {T<:AbstractFloat}
    du1, du2 = view(du, :, :, 1), view(du, :, :, 2)
    du1 .= circshift!(du1, u, [0, -1]) .- u
    du2 .= circshift!(du2, u, [-1, 0]) .- u
    du1[:, end] .= 0.0
    du2[end, :] .= 0.0
end

function grad!(du::Array{T, 4}, u::Array{T, 3}) where {T<:AbstractFloat}
    du1, du2, du3 = view(du,:,:,:,1), view(du,:,:,:,2), view(du,:,:,:,3)
    circshift!(du1, u, [0, -1, 0])
    du1 .= circshift!(du1, u, [0, -1, 0]) .- u
    du2 .= circshift!(du2, u, [-1, 0, 0]) .- u
    du3 .= circshift!(du3, u, [0, 0, -1]) .- u
    du1[:, end, :] .= 0.0
    du2[end, :, :] .= 0.0
    du3[:, :, end] .= 0.0
end

function div(p::Array{T, 3}) where {T<:AbstractFloat}
    return div2d(view(p,:,:,1), view(p,:,:,2))
end

function div!(divp::Array{T, 4}, p::Array{T, 4}) where {T<:AbstractFloat}
    p1, p2, p3 = view(p,:,:,:,1), view(p,:,:,:,2), view(p,:,:,:,3)
    p1_x, p2_y, p3_z = view(divp,:,:,:,1), view(divp,:,:,:,2), view(divp,:,:,:,3)

    # p3_z temp variable to save memory
    p1_x .= p1 .- circshift!(p3_z, p1, [0, 1, 0])
    p1_x[:, end, :] .= -p1[:, end-1, :]
    p1_x[:,   1, :] .=  p1[:, 1, :]

    p2_y .= p2 .- circshift!(p3_z, p2, [1, 0, 0])
    p2_y[end, :, :] .= -p2[end-1, :, :]
    p2_y[1,   :, :] .=  p2[1, :, :]

    p1_x .+= p2_y

    # p2_y: temporary variable to save memory
    p3_z .= p3 .- circshift!(p2_y, p3, [0, 0, 1])
    p3_z[:, :, end] .= -p3[:, :, end-1]
    p3_z[:, :,   1] .=  p3[:, :, 1]

    p1_x .+= p3_z
    return p1_x
end

function div!(divp::Array{T, 3}, p::Array{T, 3}) where {T<:AbstractFloat}
    p1, p2 = view(p,:,:,1), view(p,:,:,2)
    p1_x, p2_y = view(divp,:,:,1), view(divp,:,:,2)

    # p2_y is temp variable to save memory
    p1_x         .=  p1 - circshift!(p2_y, p1, [0, 1])
    p1_x[:, end] .= -p1[:, end-1]
    p1_x[:,   1] .=  p1[:, 1]

    circshift!(p2_y, p2, [1, 0])
    p2_y         .= p2 - p2_y
    p2_y[end, :] .= -p2[end-1, :]
    p2_y[1,   :] .=  p2[1, :]

    p1_x .+= p2_y
    return p1_x
end

"Compute backward divergence dim:[height, width]"
function div2d(p1, p2)
    p1_x         =  p1 - circshift(p1, [0, 1])
    p1_x[:, end] .= -p1[:, end-1]
    p1_x[:,   1] .=  p1[:, 1]

    p2_y = p2 - circshift(p2, [1, 0])
    p2_y[end, :] .= -p2[end-1, :]
    p2_y[1,   :] .=  p2[1, :]

    div_p = p1_x + p2_y
    return div_p
end

"Project p[H,W,2] to dual isonorm"
function proj_dual_iso!(p::Array{T, 3}, temp, weight) where {T<:AbstractFloat}
    "TODO: Avoid copying"
    p1, p2 = view(p,:,:,1), view(p,:,:,2)
    temp .= 1 ./ max.(1.0, sqrt.( p1 .^2 + p2 .^ 2 ) ./ (weight+1e-8) )
    p1 .*= temp
    p2 .*= temp
end

"Project p[H,W,Z,3] to dual isonorm"
function proj_dual_iso!(p::Array{T, 4}, temp, weight) where {T<:AbstractFloat}
    "TODO: Avoid copying"
    p1, p2, p3 = view(p,:,:,:,1), view(p,:,:,:,2), view(p,:,:,:,3)
    temp .= 1 ./ max.(1.0, sqrt.( p1 .^2 + p2 .^2 + p3 .^2 )./ (weight+1e-8))
    p1 .*= temp
    p2 .*= temp
    p3 .*= temp
end

"Project l1 norm (soft thresholding)"
function proj_l1!(x::AbstractArray{T}, weight::T) where {T<:AbstractFloat}
    x .= sign.(x) .* max.(abs.(x) .- weight, 0.0)
end

"Prox to dual l1 norm = Project to unit norm ball"
function proj_dual_l1!(x::AbstractArray{T}, weight::T) where {T<:AbstractFloat}
    x ./= max.(1.0, abs.(x) / weight)
end

"""

Proximal operator of nuclear norm to x. y=prox_{s f}(x)

"""
function prox_nuclear!(y, x, stepsize)
    H, W, M, C = size(x) # Dim : 2 for 2D, 3 for 3D
    x_HxWxCx2 = PermutedDimsArray(x, [1,2,4,3])
    x_ = reshape(x_HxWxCx2, H*W, C, 2)
    y_ = reshape(y, H*W, 2, C)

    EPS = eps()

    x_svd = mapslices(svd!, x_, dims=[2,3])
    # x_svd = dropdims(x_svd_, dims=4)
    for (hw, svd_slice) in enumerate(x_svd)
        for i=1:2
            # svd_slice.S[i] svd_slice.S .- stepsize
            if (svd_slice.S[i] > EPS)
                svd_slice.S[i] = max(svd_slice.S[i] - stepsize, 0.0) / svd_slice.S[i]
            else
                svd_slice.S[i] = max(svd_slice.S[i] - stepsize, 0.0)
            end
        end
        x_hw = view(x_, hw, :, :) # [C x 2]
        y_hw = view(y_, hw, :, :) # [2 x C]

        svd_slice.Vt[1,1] *= svd_slice.S[1]
        svd_slice.Vt[2,2] *= svd_slice.S[2]

        mul!(x_hw, x_hw, svd_slice.Vt')
        mul!(x_hw, x_hw, svd_slice.Vt)

        y_hw[1,1] = x_hw[1,1]
        y_hw[2,2] = x_hw[2,2]
        y_hw[2,1] = x_hw[1,2]
        y_hw[1,2] = x_hw[2,1]
    end
    return y
end

end

max_stars_repo_name: QPImaging/xfromProjections

max_issues_repo_name: QPImaging/xfromProjections

max_forks_repo_name: QPImaging/xfromProjections


text:
struct LocalSystem{T}
    model   ::Model{T, Triangle{T}}
    Ξ       ::PositionVector{T}
    elements::TriangleVector{T}
    umol    ::Vector{T}
    qmol    ::Vector{T}
end

function LocalSystem(model::Model{T, Triangle{T}}) where T
    cuΞ    = PositionVector([e.center for e in model.elements])
    cuelms = TriangleVector(model.elements)
    umol   = model.params.εΩ .\ φmol(model, tolerance=_etol(T))
    qmol   = model.params.εΩ .\ ∂ₙφmol(model)
    LocalSystem(model, cuΞ, cuelms, umol, qmol)
end

struct LocalSystemMatrix{T} <: PotentialMatrix{T}
    K::LaplacePotentialMatrix{T, DoubleLayer}
    α::T
    β::T
end

@inline LocalSystemMatrix(
    K     ::LaplacePotentialMatrix{T, DoubleLayer},
    params::Option{T}
) where T = LocalSystemMatrix(K, (1 + params.εΩ / params.εΣ) * T(2π), params.εΩ / params.εΣ - 1)

@inline Base.size(A::LocalSystemMatrix{T}) where T = Base.size(A.K)

@inline Base.:*(
    A::LocalSystemMatrix{T},
    x::AbstractArray{T, 1}
) where T = A.α .* x .+ (A.β .* (A.K * x))

@inline function LinearAlgebra.diag(
    A::LocalSystemMatrix{T},
    k::Int = 0
) where T
    k != 0 && error("diag not defined for k != 0 on ", typeof(A))
    A.α .* ones(T, size(A, 1))
end

function solve(sys::LocalSystem{T}) where T
    εΩ = sys.model.params.εΩ
    εΣ = sys.model.params.εΣ
    K  = LaplacePotentialMatrix{DoubleLayer}(sys.Ξ, sys.elements)
    V  = LaplacePotentialMatrix{SingleLayer}(sys.Ξ, sys.elements)
    A  = LocalSystemMatrix(K, sys.model.params)

    u = _solve_linear_system(A, K * sys.umol .- (T(2π) .* sys.umol) .- (εΩ / εΣ .* (V * sys.qmol)))
    q = _solve_linear_system(V, K * u .+ (T(2π) .* u))

    LocalBEMResult(sys.model, u, q, sys.umol, sys.qmol)
end

@inline solve(
         ::Type{LocalES},
    model::Model{T, Triangle{T}}
) where T = solve(LocalSystem(model))


text:
include("electrostatics.jl")
include("Partridge_Schwenke.jl")
include("utilities.jl")

mutable struct TTM3F <: TTM_Potential
    num_waters::Int

    M_site_coords::Vector{SVector{3, Float64}}
    constants::TTM3_Constants
    elec_data::Electrostatics
end

TTM3F(num_waters::Int) = TTM3F(num_waters, [@SVector zeros(3) for _ in 1:(num_waters > 0 ? 4*num_waters : 1)], TTM3_Constants(), Electrostatics(zeros(Int(num_waters * 4)), TTM3_Constants()))
TTM3F() = TTM3F(0)

function evaluate!(ttm3f::TTM3F, coords::AbstractMatrix{Float64}, grads::Union{Matrix{Float64}, Nothing}=nothing)
    """
    Evaluates the ttm3f energy of the system containing water called coords.
    """
    @assert isinteger(size(coords, 2) / 3) "Number of coordinates not divisble by three. Is this water?"
    if (ttm3f.num_waters != Int(size(coords, 2) / 3))
        temp_ttm = TTM3F(Int(size(coords, 2) / 3))
        ttm3f.num_waters = temp_ttm.num_waters
        ttm3f.M_site_coords = temp_ttm.M_site_coords
        ttm3f.elec_data = temp_ttm.elec_data
    end

    compute_M_site_crd!(ttm3f, coords)
    #Rij_vectors = distance_vectors(M_site_coords)
    # @SPEED Ignore all non M-Site calculations in the electrostatics calculation !!!!!
    # which is only executed if the gradients are actually requested
    grads_E = [@MVector zeros(3) for _ in 1:3 * 4 * ttm3f.num_waters] # O H H M
    grads_q = [@MArray zeros(3,3,3) for _ in 1:ttm3f.num_waters]
    # calculate the INTRA-molecular energy and Dipole Moment
    # Surface (Partridge-Schwenke)
    E_int = PS_energies_and_gradients!(ttm3f, coords, grads_q, grads, true)

    ttm3f.elec_data.q *= chargecon()
    if grads !== nothing
        grads_q *= chargecon()
    end

    #---------------------------------------------------------------!
    # Calculate the vdW interactions for all atoms                  !
    #---------------------------------------------------------------!

    E_vdw::Float64 = 0.0
    # Loop through all of the waters for van der Waal's interactions
    for i in 1:ttm3f.num_waters-1
        for j in i+1:ttm3f.num_waters
            # Calculate Oxygen-Oxygen interactions
            @inbounds @views Rij = coords[:, 3*(i-1)+1] - coords[:, 3*(j-1)+1]
            expon = ttm3f.constants.vdwD * exp(-ttm3f.constants.vdwE * norm(Rij))

            # Add Evdw van der Waal energy for O-O interaction
            E_vdw += ttm3f.constants.vdwC / norm(Rij)^6 + expon

            if (grads !== nothing)
                tmp::Float64 = (-(6.0 * ttm3f.constants.vdwC / norm(Rij)^6) / (Rij ⋅ Rij) - ttm3f.constants.vdwE * expon / norm(Rij))

                # add forces to indices for the M site
                @inbounds @views grads[:, 3*(i-1)+1] += tmp * Rij
                @inbounds @views grads[:, 3*(j-1)+1] -= tmp * Rij
            end
        end
    end

    # if no gradients then just return energy here.
    if (grads === nothing)
        return E_int + E_vdw + electrostatics(ttm3f.elec_data, ttm3f.M_site_coords, nothing, true)
    end

    #-------------------------------------------------------------------------!
    # Calculate the remaining part of the derivatives                         !
    #-------------------------------------------------------------------------!

    E_elec::Float64 = electrostatics(ttm3f.elec_data, ttm3f.M_site_coords, grads_E, true)

    add_electrostatic_and_dms_derivatives!(ttm3f, grads, grads_E, grads_q)

    return E_int + E_vdw + E_elec
end

max_stars_repo_name: heindelj/MolecularPotentials

max_issues_repo_name: heindelj/MolecularPotentials

max_forks_repo_name: heindelj/MolecularPotentials


text:
include(joinpath(@__DIR__, "..", "CoordTrans", "CoordTrans.jl"))
include(joinpath(@__DIR__, "..", "CoordTrans", "CoordTrans_ESMF.jl"))

using .CoordTrans
using .CoordTrans_ESMF
using ArgParse
using JSON

function parse_commandline()

    s = ArgParseSettings()
    @add_arg_table s begin

        "--w-file"
            help = "Weighting file."
            arg_type = String
            required = true

        "--s-file"
            help = "Input source file. Filenames separated by commas."
            arg_type = String
            required = true

        "--d-file"
            help = "Output destination file. Filenames separated by commas."
            arg_type = String
            required = true

        "--vars"
            help = "Variable names list. They should by dimension 2 (x, y) or 3 (x, y, z) with or without record (time) dimension. Ex: --vars=Ts,Ss,MLD"
            arg_type = String

        "--x-dim"
            help = "Variable name of x-dimension."
            arg_type = String
            required = true

        "--y-dim"
            help = "Variable name of y-dimension."
            arg_type = String
            required = true

        "--z-dim"
            help = "Variable name of z-dimension."
            arg_type = String

        "--t-dim"
            help = "Variable name of time-dimension."
            arg_type = String
            default = "time"

        "--t-len"
            help = "If this is not -1 then time will be a non-record dimension with length specified here."
            arg_type = Int64
            default = -1

        "--s-dir"
            help = "Source files directory. If unset then this variable will not be applied. Convinent when working with multiple files"
            arg_type = String
            default = ""

        "--d-dir"
            help = "Destination files directory. If unset then this variable will not be applied. Convinent when working with multiple files"
            arg_type = String
            default = ""

        "--algo"
            help = "Type of algorithm. Valid values: `XTT`, `ESMF`."
            arg_type = String
            default = "ESMF"

    end

    return parse_args(ARGS, s)
end

println("Running ", @__FILE__)

parsed = parse_commandline()
print(json(parsed, 4))

if parsed["vars"] != nothing
    varnames = collect(split(parsed["vars"], ","; keepempty=false))
else
    varnames = nothing
end

s_files = collect(split(parsed["s-file"], ","; keepempty=false))
d_files = collect(split(parsed["d-file"], ","; keepempty=false))

if length(s_files) != length(d_files)
    throw(ErrorException("Numbers of source files and destination files do not match."))

end

if parsed["algo"] == "XTT"
    wi = CoordTrans.readWeightInfo(parsed["w-file"])
elseif parsed["algo"] == "ESMF"
    wi = CoordTrans_ESMF.readWeightInfo(parsed["w-file"])
end

for i in 1:length(s_files)

    s_file = ( parsed["s-dir"] == "" ) ? s_files[i] : joinpath(parsed["s-dir"], s_files[i])
    d_file = ( parsed["d-dir"] == "" ) ? d_files[i] : joinpath(parsed["d-dir"], d_files[i])
    println("Trasform: ", s_file, " => ", d_file)

    if parsed["algo"] == "XTT"

        println("### Algo: XTT ###")
        CoordTrans.convertFile(
            s_file,
            d_file,
            wi,
            varnames=varnames;
            xdim = parsed["x-dim"],
            ydim = parsed["y-dim"],
            zdim = parsed["z-dim"],
            tdim = parsed["t-dim"],
            tlen = parsed["t-len"]
        )

    elseif parsed["algo"] == "ESMF"
        println("### Algo: ESMF ###")

        CoordTrans_ESMF.convertFile(
            s_file,
            d_file,
            wi,
            varnames=varnames;
            xdim = parsed["x-dim"],
            ydim = parsed["y-dim"],
            zdim = parsed["z-dim"],
            tdim = parsed["t-dim"],
            tlen = parsed["t-len"]
        )

    end
end

max_stars_repo_name: meteorologytoday/IOM

max_issues_repo_name: meteorologytoday/IOM

max_forks_repo_name: meteorologytoday/IOM


text:
using PyCall

export xr, np, cfg, tasks, workflow, MBsandbox

# Import OGGM sub-libraries in Julia
netCDF4 = pyimport("netCDF4")
cfg = pyimport("oggm.cfg")
utils = pyimport("oggm.utils")
workflow = pyimport("oggm.workflow")
tasks = pyimport("oggm.tasks")
graphics = pyimport("oggm.graphics")
bedtopo = pyimport("oggm.shop.bedtopo")
MBsandbox = pyimport("MBsandbox.mbmod_daily_oneflowline")

# Essential Python libraries
np = pyimport("numpy")
xr = pyimport("xarray")
# matplotlib = pyimport("matplotlib")
# matplotlib.use("Qt5Agg")


text:
# importance sampling

struct Metropolis end
struct MetropolisHastings end # todo
struct HeatBath end # todo
struct Glauber end # todo
struct RealMicrocanonical end # todo

#TODO: generalize in the following way
#eval(log_difference::Float64)::Bool
#@inline Base.@propagate_inbounds function _evaluate(d::UnionMetrics, a::Union{Array, ArraySlice}, b::Union{Array, ArraySlice})
#@inline function _evaluate(d::UnionMetrics, a::AbstractArray, b::AbstractArray)

@inline function _evaluate(rng::AbstractRNG, log_difference::Float64)::Bool
    if log_difference > 0
        return true
    elseif rand(rng) < exp(log_difference)
        return true
    else
        return false
    end
end

@doc """
    accept(log_weight::Function, args_new::Tuple{Number, N}, args_old::Tuple{Number, N}, rng::AbstractRNG)::Bool where N

Evaluate most general acceptance probability for imporance sampling of ``P(x) \\propto e^{log\\_weight(x)}``.

# Arguments
- `log_weight(args)`: logarithmic ensemble weight function, e.g., canomical ensemble ``log\\_weight(x) = -\\beta x``
- `args_new`: arguments (can be Number or Tuple) for new (proposed) state
- `args_old`: arguments (can be Number or Tuple) for old                        state
- `rng`: random number generator, e.g. MersenneTwister

# Specializations
- accept(alg::Metropolis(), rng::AbstractRNG, x_new::T, x_old::T) where T (@ref)
"""
function accept(rng::AbstractRNG, log_weight::Function, args_new::NTuple{N,T}, args_old::NTuple{N,T})::Bool where {N,T}
    log_difference = log_weight(args_new...) - log_weight(args_old...)
    return _evaluate(rng, log_difference)
end

function accept(rng::AbstractRNG, log_weight::Function, args_new::T, args_old::T)::Bool where T
    log_difference = log_weight(args_new) - log_weight(args_old)
    return _evaluate(rng, log_difference)
end

@doc """
    accept(alg::Metropolis, rng::AbstractRNG, beta::Float64, dx::T)::Bool where T

Standard metropolis algorithm with
```p(x\\to x^') = \\text{min}\\left(1, e^{-\\beta \\Delta x}\\right)```
where `dx = x^' - x`
"""
function accept(alg::Metropolis, rng::AbstractRNG, beta::Float64, dx::T)::Bool where T
    if dx < 0
        return true
    elseif rand(rng) < exp(-beta * dx)
        return true
    else
        return false
    end
end

"""
    accept(alg::MetropolisHastings, rng::AbstractRNG, beta::Float64, dx::T, q_old::Float64, q_new::Float64)::Bool where T
"""
function accept(alg::MetropolisHastings, rng::AbstractRNG, beta::Float64, dx::T, q_old::Float64, q_new::Float64)::Bool where T
    println("not implemented yet because API uncleaer")
    return false
end

# sweep
@doc """
    sweep(list_updates, list_weights::AbstractWeights, rng::AbstractRNG; number_updates::Int=1) where T<:AbstractFloat

Randomly pick und run update (has to check acceptance by itself!) from
`list_updates` with probability specified in `list_probabilities` and repeat
this `number_updates` times.

# Remarks
* Walker's alias method?
"""
function sweep(list_updates, list_weights::AbstractWeights, rng::AbstractRNG; number_updates::Int = 1)
    @assert length(list_updates) == length(list_weights)

    for i in 1:number_updates
        id = StatsBase.sample(rng, list_weights)
        # update is requred to call Metropolis.accept() itself
        list_updates[id]()
    end
end

# TODO: Should we keep this specialization?
function sweep(update::Function, rng::AbstractRNG; number_updates::Int = 1)
    for i in 1:number_updates
        # update is requred to call Metropolis.accept() itself
        update()
    end
end


text:
using Plots, MAT

nprocs = (2, 2) # nprocs (x, y) dim
nt     = 100

@views function vizme2D_mpi(nprocs, nt)
    dname = "tmp_gif"; if isdir("./out2D/$(dname)")==false mkdir("./out2D/$(dname)") end; loadpath = "./out2D/$(dname)"; anim = Animation(loadpath,String[])
    for it = 1:nt
        T  = []
        ip = 1
        for ipx = 1:nprocs[1]
            for ipy = 1:nprocs[2]
                file = matopen("./out2D/T_$(ip-1)_$(it).mat"); T_loc = read(file, "T"); close(file)
                nx_i, ny_i = size(T_loc,1)-2, size(T_loc,2)-2
                ix1, iy1   = 1+(ipx-1)*nx_i, 1+(ipy-1)*ny_i
                if (ip==1)  T = zeros(nprocs[1]*nx_i, nprocs[2]*ny_i)  end
                T[ix1:ix1+nx_i-1,iy1:iy1+ny_i-1] .= T_loc[2:end-1,2:end-1]
                ip += 1
            end
        end
        heatmap(T', framestyle=:box, aspect_ratio=1, clims=(0.0,1.0), xlims=(1, size(T,1)), ylims=(1, size(T,2)), c=:hot, title="heat diffusion step=$(it)"); frame(anim)
    end
    gif(anim, "heat_2D_mpi_$(nprocs[1]*nprocs[2])procs.gif", fps = 15)
    return
end

vizme2D_mpi(nprocs, nt)

max_stars_repo_name: luraess/geo-hpc-course

max_issues_repo_name: luraess/geo-hpc-course

max_forks_repo_name: luraess/geo-hpc-course


text:
# integration segment (a,b), estimated integral I, and estimated error E
struct Segment{TX,TI,TE}
    a::TX
    b::TX
    I::TI
    E::TE
end
Base.@pure Base.promote_rule(::Type{Segment{TX,TI,TE}}, ::Type{Segment{TX′,TI′,TE′}}) where {TX,TI,TE,TX′,TI′,TE′} =
    Segment{promote_type(TX,TX′), promote_type(TI,TI′), promote_type(TE,TE′)}
Base.convert(::Type{T}, s::Segment) where {T<:Segment} = T(s.a,s.b,s.I,s.E)
Base.isless(i::Segment, j::Segment) = isless(i.E, j.E)

# Internal routine: approximately integrate f(x) over the interval (a,b)
# by evaluating the integration rule (x,w,gw). Return a Segment.
function evalrule(f, a,b, x,w,gw, nrm)
    # Ik and Ig are integrals via Kronrod and Gauss rules, respectively
    s = convert(eltype(x), 0.5) * (b-a)
    n1 = 1 - (length(x) & 1) # 0 if even order, 1 if odd order
    # unroll first iterationof loop to get correct type of Ik and Ig
    fg = f(a + (1+x[2])*s) + f(a + (1-x[2])*s)
    fk = f(a + (1+x[1])*s) + f(a + (1-x[1])*s)
    Ig = fg * gw[1]
    Ik = fg * w[2] + fk * w[1]
    for i = 2:length(gw)-n1
        fg = f(a + (1+x[2i])*s) + f(a + (1-x[2i])*s)
        fk = f(a + (1+x[2i-1])*s) + f(a + (1-x[2i-1])*s)
        Ig += fg * gw[i]
        Ik += fg * w[2i] + fk * w[2i-1]
    end
    if n1 == 0 # even: Gauss rule does not include x == 0
        Ik += f(a + s) * w[end]
    else # odd: don't count x==0 twice in Gauss rule
        f0 = f(a + s)
        Ig += f0 * gw[end]
        Ik += f0 * w[end] +
            (f(a + (1+x[end-1])*s) + f(a + (1-x[end-1])*s)) * w[end-1]
    end
    Ik_s, Ig_s = Ik * s, Ig * s # new variable since this may change the type
    E = nrm(Ik_s - Ig_s)
    if isnan(E) || isinf(E)
        throw(DomainError(a+s, "integrand produced $E in the interval ($a, $b)"))
    end
    return Segment(oftype(s, a), oftype(s, b), Ik_s, E)
end

# compute result = f(x1) + f(x2) in-place
function eval2x!(result, f::InplaceIntegrand{F}, x1, x2) where {F}
    f.f!(result, x1)
    f.f!(f.fx, x2)
    result .+= f.fx
end

# as above, but call assume a mutable result type (e.g. an array) and
# act in-place using `f!(result, x)`.
function evalrule(f::InplaceIntegrand{F}, a,b, x,w,gw, nrm) where {F}
    # Ik and Ig are integrals via Kronrod and Gauss rules, respectively
    s = convert(eltype(x), 0.5) * (b-a)
    n1 = 1 - (length(x) & 1) # 0 if even order, 1 if odd order
    fg, fk, Ig, Ik = f.fg, f.fk, f.Ig, f.Ik # pre-allocated temporary arrays
    # unroll first iterationof loop to get correct type of Ik and Ig
    eval2x!(fg, f, a + (1+x[2])*s, a + (1-x[2])*s)
    eval2x!(fk, f, a + (1+x[1])*s, a + (1-x[1])*s)
    Ig .= fg .* gw[1]
    Ik .= fg .* w[2] .+ fk .* w[1]
    for i = 2:length(gw)-n1
        eval2x!(fg, f, a + (1+x[2i])*s, a + (1-x[2i])*s)
        eval2x!(fk, f, a + (1+x[2i-1])*s, a + (1-x[2i-1])*s)
        Ig .+= fg .* gw[i]
        Ik .+= fg .* w[2i] .+ fk .* w[2i-1]
    end
    f.f!(f.fx, a + s)
    if n1 == 0 # even: Gauss rule does not include x == 0
        Ik .+= f.fx .* w[end]
    else # odd: don't count x==0 twice in Gauss rule
        Ig .+= f.fx .* gw[end]
        f.f!(f.fg, a + (1+x[end-1])*s)
        f.f!(f.fk, a + (1-x[end-1])*s)
        Ik .+= f.fx .* w[end] .+ (f.fg + f.fk) .* w[end-1]
    end
    Ik_s = Ik * s # new variable since this may change the type
    f.Idiff .= Ik_s .- Ig .* s
    E = nrm(f.Idiff)
    if isnan(E) || isinf(E)
        throw(DomainError(a+s, "integrand produced $E in the interval ($a, $b)"))
    end
    return Segment(oftype(s, a), oftype(s, b), Ik_s, E)
end


text:
using HTTP: HTTP, URI, request
using JSON3: JSON3
using DataFrames: DataFrames, DataFrame, combine, groupby, ByRow, select!, transform!, innerjoin
using Statistics: Statistics, mean
using CSV: CSV

const API_BLS = get(ENV, "api_bls_token", "")

response = request("POST",
                   URI(scheme = "https",
                       host = "api.bls.gov",
                       path = "/publicAPI/v2/timeseries/data",
                       query = ["registrationkey" => API_BLS,
                                "startyear" => 2008,
                                "endyear" => 2019,
                                "annualaverage" => true,
                                "seriesid" => "CMU2020000120000D,CMU2030000120000D"]),
                   ["Content-Type" => "application/x-www-form-urlencoded"])
json = JSON3.read(response.body)
function parse_series(series)
    id = series.seriesID
    data = DataFrame((year = elem.year, value = parse(Float64, elem.value)) for elem in series.data)
    data = combine(groupby(data, :year), :value => mean => id)
end
x = parse_series(json.Results.series[1])
y = parse_series(json.Results.series[2])
z = innerjoin(x, y, on = :year)
transform!(z,
           [:CMU2020000120000D, :CMU2030000120000D] => ByRow((sw, tb) -> (sw + tb) / sw) =>
           :salary_to_total_compensation)
sort!(select!(z, [:year, :salary_to_total_compensation]))
CSV.write(joinpath("data", "bls_salary_wages_to_total_compensation.csv"), z)

max_stars_repo_name: Nosferican/Estimating-Resource-Costs

max_issues_repo_name: Nosferican/Estimating-Resource-Costs

max_forks_repo_name: Nosferican/Estimating-Resource-Costs


text:
filename = "q1_q17_dr25_koi.csv"
#filename = "q1_q17_dr25_stellar.csv"
(data,header) = readcsv(filename,header=true)

ncols = length(header)
data_dict = Dict{String,Any}()
for i in 1:ncols
  data_vec = data[:,i]
  try
    data_vec = convert(Vector{Float64},data[:,i])
    #dat_dict[header[i]] = data_vec
    println("# Converted col:", header[i])
  catch
    #println("# Couldn't convert colm: ",header[i])
  end
end


text:
@testset "diffusion setup object" begin
    param = :complexConjug
    θ₀ = [10.0, -5.0, 5.0, 0.0, 3.0]
    P˟ = BSI.FitzhughDiffusion(param, θ₀...)
    obs = [[0.0, 0.0], [-0.3, -0.2], [0.4, 0.2]]
    tt = [0.0, 1.0, 1.5]
    P̃ = [BSI.FitzhughDiffusionAux(param, θ₀..., tt[1], obs[1], tt[2], obs[2]),
         BSI.FitzhughDiffusionAux(param, θ₀..., tt[2], obs[2], tt[3], obs[3])]
    setup = BSI.DiffusionSetup(P˟, P̃, PartObs())

    @testset "initialisation" begin
        @test setup.P˟ == P˟
        @test setup.P̃ == P̃
        @test typeof(setup.adaptive_prop) <: Adaptation{Val{false}}
        @test setup.skip_for_save == 1
        @test !any(values(setup.setup_completion))
        @test @suppress !BSI.check_if_complete(setup, [:obs])
        @test @suppress !BSI.check_if_complete(setup, [:imput])
        @test @suppress !BSI.check_if_complete(setup, [:prior])
    end

    L = [1. 0.; 0. 1.]
    Σ = [0.5 0.0; 0.0 1.0]
    BSI.set_observations!(setup, [L, L], [Σ, 2*Σ], obs, tt)

    @testset "setting observations" begin
        @test setup.Ls == [L, L]
        @test setup.Σs == [Σ, 2*Σ]
        @test setup.obs == obs
        @test setup.obs_times == tt
        @test setup.fpt == [nothing, nothing]
        @test @suppress !BSI.check_if_complete(setup, [:imput])
        @test @suppress !BSI.check_if_complete(setup, [:prior])
        @test @suppress BSI.check_if_complete(setup, [:obs])
    end

    dt = 0.01
    τ(t₀,T) = (x) ->  t₀ + (x-t₀) * (2-(x-t₀)/(T-t₀))
    BSI.set_imputation_grid!(setup, dt)

    @testset "setting imputation grid" begin
        @test setup.dt == dt
        @test setup.τ(tt[1], tt[2])(0.5*(tt[1]+tt[2])) == τ(tt[1], tt[2])(0.5*(tt[1]+tt[2]))
        @test @suppress !BSI.check_if_complete(setup, [:prior])
        @test @suppress BSI.check_if_complete(setup, [:obs, :imput])
    end

    x0_prior = BSI.KnownStartingPt(obs[1])
    BSI.set_x0_prior!(setup, x0_prior)
    @testset "setting priors" begin
        @test setup.x0_prior == x0_prior
        @test @suppress BSI.check_if_complete(setup, [:obs, :imput, :prior])
    end

    @testset "determining data type" begin
        @test BSI.determine_data_type(setup) == (SArray{Tuple{2},Float64,1,2},
                                             Float64)
    end

    BSI.prepare_containers!(setup)
    @testset "setting internal containers" begin
        @test setup.Wnr == Wiener{Float64}()
        @test eltype(setup.WW) == SamplePath{Float64}
        @test length(setup.WW) == 2
        @test eltype(setup.XX) == SamplePath{SArray{Tuple{2},Float64,1,2}}
        @test length(setup.XX) == 2
        @test setup.Ls == [SMatrix{2,2}(L), SMatrix{2,2}(L)]
        @test setup.Σs == [SMatrix{2,2}(Σ), SMatrix{2,2}(2*Σ)]
        @test setup.obs == map(x->SVector{2}(x), obs)
    end

    BSI.initialise!(Float64, setup, BSI.Vern7(), false, BSI.NoChangePt())
    @testset "initialisation of proposal law" begin
        @test length(setup.P) == 2
        @test typeof(setup.P[1].Pt) == typeof(setup.P̃[1])
        @test setup.P[1].Target == setup.P˟
    end
end


text:
# MIT License
#
# Copyright (c) 2020 Joshua E Gentry

# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.

# =================================================================================================
"""
    pix_read(
        filename::AbstractString;
        jpgLuminance::Bool = false,
        jpgFailOnBadData::Bool = false
    )::Union{Pix, Nothing}

Load an image from from disk of an unspecified type.  Returns `nothing` if the file could not be
loaded.

__Arguments:__

| T | Name             | Default | Description
|:--| :--------------- | :------ | :----------
| R | filename         |         | The filename of the image to load.
| O | jpgLuminance     | `false` | If the image is a JPEG should we only load luminance data?
| O | jpgFailOnBadData | `false` | If the image is a JPEG should we ignore bit errors?

__Details:__

The two optional parameters are obviously ignored if the image is not a JPEG.
"""
function pix_read(
            filename::AbstractString;
            jpgLuminance::Bool = false,
            jpgFailOnBadData::Bool = false
        )::Union{Pix, Nothing}
    local hint = Int32(0)

    if jpgLuminance == true
        hint |= 0x1
    end
    if jpgFailOnBadData == true
        hint |= 0x2
    end

    local ptr = @threadcall(
        (:pixReadWithHint, LEPTONICA),
        Ptr{Cvoid},
        (Cstring, Cint),
        filename,
        hint
    )

    if ptr == C_NULL
        return nothing
    end

    return Pix(ptr)
end

# =================================================================================================
"""
    pix_read(
        stream::IO
    )::Union{Pix, Nothing}

Load an image from from a stream.  Returns `nothing` if the file could not be loaded.

__Arguments:__

| T | Name   | Default | Description
|:--| :----- | :------ | :----------
| R | stream |         | The stream to read the image from.

__Details:__

This implementation mirrors the API provided by Leptonica when you pass in a FILE pointer. This
assumes that the remainder of the stream contains an image, so the remainder of the stream will be
read and decoded as an image.
"""
function pix_read(
            stream::IO
        )::Union{Pix, Nothing}
    local data = read(stream)
    return pix_read(data)
end

# =================================================================================================
"""
    pix_read(
        data::AbstractArray{UInt8}
    )::Union{Pix, Nothing}

Load an image from from a byte array.  Returns `nothing` if the file could not be loaded.

__Arguments:__

| T | Name | Default | Description
|:--| :--- | :------ | :----------
| R | data |         | The byte array to read the image from.
"""
function pix_read(
            data::AbstractArray{UInt8}
        )::Union{Pix, Nothing}
    local ptr = ccall(
        (:pixReadMem, LEPTONICA),
        Ptr{Cvoid},
        (Ptr{Cvoid}, Csize_t),
        data,
        length(data)
    )

    if ptr == C_NULL
        return nothing
    end

    return Pix(ptr)
end

# =================================================================================================
"""
    pix_write_i(
        pix::Pix,
        format::IFF
    )::Union{Pix, Nothing}

Write an image to memory in the specified format.  If there is an error `(C_NULL, 0)` is returned.

__Arguments:__

| T | Name   | Default  | Description
|:--| :----- | :------- | :----------
| R | pix    |          | The image to write to memory.
| R | format |          | The image format to write the image out as.

__Details:__

**This method is meant for internal use only.**  It returns the memory pointer and number of bytes
allocated by the Leptonica library to write out the image.  The caller is responsible for freeing
this memory with a call to `lept_free()`.
"""
function pix_write_i(
            pix::Pix,
            format::IFF
        )::Tuple{Ptr{UInt8}, Csize_t}
    local temp = Ref(Ptr{UInt8}(C_NULL))
    local size = Ref(Csize_t(0))

    if is_valid(pix) == false
        @error "Pix has been freed."
        return (C_NULL, 0)
    end

    local retval = ccall(
        (:pixWriteMem, LEPTONICA),
        Cint,
        (Ptr{Ptr{UInt8}}, Ptr{Csize_t}, Ptr{Cvoid}, Cint),
        temp,
        size,
        pix,
        format
    )

    if retval == 1
        return (C_NULL, 0)
    end

    return (temp[], size[])
end

# =================================================================================================
"""
    pix_write(
        filename::String,
        pix::Pix,
        format::IFF = IFF_DEFAULT
    )::Bool

Write an image to disk in the specified format.  If there is an error `false` is returned.

__Arguments:__

| T | Name     | Default       | Description
|:--| :------- | :------------ | :----------
| R | filename |               | The filename to write the image to.
| R | pix      |               | The image to write to disk.
| O | format   | `IFF_DEFAULT` | The image format to write the image out as.

__Details:__

If the file exists it will be overwritten.

If format isn't specified Leptonica will choose the correct format to write the image as based on
the input format.
"""
function pix_write(
            filename::String,
            pix::Pix,
            format::IFF = IFF_DEFAULT
        )::Bool

    if is_valid(pix) == false
        @error "Pix has been freed."
        return false
    end

    local retval = @threadcall(
        (:pixWrite, LEPTONICA),
        Cint,
        (Cstring, Ptr{Cvoid}, Cint),
        filename,
        pix,
        format
    )

    return retval == 0
end

# =================================================================================================
"""
    pix_write(
        stream::IO,
        pix::Pix,
        format::IFF = IFF_DEFAULT
    )::Bool

Write an image to the IO stream in the specified format.  If there is an error `false` is returned.

__Arguments:__

| T | Name   | Default       | Description
|:--| :----- | :------------ | :----------
| R | stream |               | The stream to write the image to.
| R | pix    |               | The image to write to the stream.
| O | format | `IFF_DEFAULT` | The image format to write the image out as.

__Details:__

If format isn't specified Leptonica will choose the correct format to write the image as based on
the input format.
"""
function pix_write(
            stream::IO,
            pix::Pix,
            format::IFF = IFF_DEFAULT
        )::Bool
    local result = false
    local data, size = pix_write_i(pix, format)

    if data != C_NULL
        try
            if size > 0
                unsafe_write(stream, data, size)
                result = true
            end
        finally
            lept_free(data)
        end
    end

    return result
end

# =================================================================================================
"""
    pix_write(
        pix::Pix,
        format::IFF = IFF_DEFAULT
    ):::Union{Vector{UInt8}, Nothing}

Write an image to a byte array in the specified format.  If there is an error `nothing` is returned.

__Arguments:__

| T | Name   | Default       | Description
|:--| :----- | :------------ | :----------
| R | pix    |               | The image to write to the stream.
| O | format | `IFF_DEFAULT` | The image format to write the image out as.

__Details:__

If format isn't specified Leptonica will choose the correct format to write the image as based on
the input format.
"""
function pix_write(
            pix::Pix,
            format::IFF = IFF_DEFAULT
        )::Union{Vector{UInt8}, Nothing}
    local output     = nothing
    local data, size = pix_write_i(pix, format)

    if data != C_NULL
        try
            if size > 0
                output = Vector{UInt8}(undef, size)
                unsafe_copyto!(pointer(output), data, size)
            end
        finally
            lept_free(data)
        end
    end

    return output
end

# =================================================================================================
"""
    pix_write_implied_format(
        filename::AbstractString,     # The name of the file to write to.
        pix::Pix;                     # The image to write.
        quality::Integer  = Int32(75), # If the image is JPEG the quality to encode as.
        progressive::Bool = false     # If the image is JEPG use progressive encoding?
    )::Bool

Write the image to a file based on the file extension.  The quality and progressive parameters are
only used if the image type is JPEG.  Returns true on success or false on failure.
"""
function pix_write_implied_format(
            filename::AbstractString,
            pix::Pix;
            quality::Integer = Int32(75),
            progressive::Bool = false
        )::Bool
    if is_valid(pix) == false
        @error "Pix has been freed."
        return false
    end
    if quality < 1 || quality > 100
        @error "Quality must be between 1 and 100 inclusive."
        return false
    end

    local retval = @threadcall(
        (:pixWriteImpliedFormat, LEPTONICA),
        Cint,
        (Cstring, Ptr{Cvoid}, Cint, Cint),
        filename,
        pix,
        quality,
        progressive ? 1 : 0
    )

    return retval == 0
end


text:
include(joinpath(dirname(@__FILE__), "ExtensionTests.jl"))

Docile.Cache.register_module(joinpath(dirname(@__FILE__), "ExtensionTests.jl"))

import ExtensionTests
import Docile: Cache, Formats

Cache.getmeta(ExtensionTests)[:format] = Formats.PlaintextFormatter

facts("Extensions.") do

    context("!!set and !!get") do

        @fact Cache.getparsed(qs(ExtensionTests, :set_and_get)...)      --> "test"
        @fact Cache.getmeta(qs(ExtensionTests, :set_and_get)...)[:name] --> "test"

    end

    context("!!var") do

        @fact Cache.getparsed(qs(ExtensionTests, :var)...)      --> "test"
        @fact Cache.getmeta(qs(ExtensionTests, :var)...)[:name] --> "test"

    end

    context("!!summary") do

        @fact Cache.getparsed(qs(ExtensionTests, :summary)...)         --> "summary"
        @fact Cache.getmeta(qs(ExtensionTests, :summary)...)[:summary] --> "summary"

    end

    context("!!longform") do

        @fact Cache.getparsed(qs(ExtensionTests, :longform)...) --> "..."

    end

    context("!!include") do

        @fact Cache.getparsed(qs(ExtensionTests, :includes)...)         --> "one\n"
        @fact Cache.getparsed(qs(ExtensionTests, :set_includes_get)...) --> " one\n test"

    end

    context("!!include with metas") do

        @fact Cache.getparsed(qs(ExtensionTests, :includes_nested)...) -->
            string("Text with Outer\nLevel1\nLevel2\nLevel3\nLevel4\nLevel5\nLevel6\n\n",
                   "Some unicode 所以不多说了 бежал мета\n Get the inner again: бежал мета")

        @fact Cache.getmeta(qs(ExtensionTests, :includes_nested)...)[:level6] --> "Level6"

    end

end


text:
include("_test_1_lti_primitives.jl")
include("_test_2_lti_systems.jl")
include("_test_3_lti_kernel.jl")


text:
"""
## fmplat!

This subroutine forms the 2nd derivatives for rectangular
plate bending fin_els.

### Method
```julia
fmplat!(d2x, d2y, d2xy, points, aa, bb, i)
```

### Arguments
```julia
* d2x::Vector{Float64}       : x derivative term (Updated)
* d2y::Vector{Float64}       : y derivative term (Updated)
* d2xy::Vector{Float64}      : x,y derivative term (Updated)
* points::Matrix{Float64}    : Location of Gauss points
* aa::Float64                : Dimension of plate
* bb::Float64                : Dimension of plate
* i::Int                   : Gauss point to use
```
"""
function fmplat!(d2x::Vector{Float64}, d2y::Vector{Float64}, d2xy::Vector{Float64},
  points::Matrix{Float64}, aa::Float64, bb::Float64, i::Int)
  #
  # This subroutine forms the 2nd derivatives for rectangular
  # plate bending fin_els.
  #
  # Updates d2x, d2y and d2xy vectors
  #
 x = points[i,1]
 e = points[i,2]
 xp1 = x+1.0
 xp12 = xp1*xp1
 xp13 = xp12*xp1
 ep1 = e+1.0
 ep12 = ep1*ep1
 ep13 = ep12*ep1
 p1 = 1.0-0.75*xp12+0.25*xp13
 q1 = 1.0-0.75*ep12+0.25*ep13
 p2 = 0.5*aa*xp1*(1.0-xp1+0.25*xp12)
 q2 = 0.5*bb*ep1*(1.0-ep1+0.25*ep12)
 p3 = 0.25*xp12*(3.0-xp1)
 q3 = 0.25*ep12*(3.0-ep1)
 p4 = 0.25*aa*xp12*(0.5*xp1-1.0)
 q4 = 0.25*bb*ep12*(0.5*ep1-1.0)
 dp1 = 1.5*xp1*(0.5*xp1-1.0)
 dq1 = 1.5*ep1*(0.5*ep1-1.0)
 dp2 = aa*(0.5-xp1+0.375*xp12)
 dq2 = bb*(0.5-ep1+0.375*ep12)
 dp3 = 1.5*xp1*(1.0-0.5*xp1)
 dq3 = 1.5*ep1*(1.0-0.5*ep1)
 dp4 = 0.5*aa*xp1*(0.75*xp1-1.0)
 dq4 = 0.5*bb*ep1*(0.75*ep1-1.0)
 d2p1 = 1.5*x
 d2p2 = 0.25*aa*(3.0*x-1.0)
 d2p3 = -d2p1
 d2p4 = 0.25*aa*(3.0*x+1.0)
 d2q1 = 1.5*e
 d2q2 = 0.25*bb*(3.0*e-1.0)
 d2q3 = -d2q1
 d2q4 = 0.25*bb*(3.0*e+1.0)
 d2x[1] = d2p1*q1
 d2x[2] = d2p2*q1
 d2x[3] = d2p1*q2
 d2x[4] = d2p2*q2
 d2x[5] = d2p1*q3
 d2x[6] = d2p2*q3
 d2x[7] = d2p1*q4
 d2x[8] = d2p2*q4
 d2x[9] = d2p3*q3
 d2x[10] = d2p4*q3
 d2x[11] = d2p3*q4
 d2x[12] = d2p4*q4
 d2x[13] = d2p3*q1
 d2x[14] = d2p4*q1
 d2x[15] = d2p3*q2
 d2x[16] = d2p4*q2
 d2y[1] = p1*d2q1
 d2y[2] = p2*d2q1
 d2y[3] = p1*d2q2
 d2y[4] = p2*d2q2
 d2y[5] = p1*d2q3
 d2y[6] = p2*d2q3
 d2y[7] = p1*d2q4
 d2y[8] = p2*d2q4
 d2y[9] = p3*d2q3
 d2y[10] = p4*d2q3
 d2y[11] = p3*d2q4
 d2y[12] = p4*d2q4
 d2y[13] = p3*d2q1
 d2y[14] = p4*d2q1
 d2y[15] = p3*d2q2
 d2y[16] = p4*d2q2
 d2xy[1] = dp1*dq1
 d2xy[2] = dp2*dq1
 d2xy[3] = dp1*dq2
 d2xy[4] = dp2*dq2
 d2xy[5] = dp1*dq3
 d2xy[6] = dp2*dq3
 d2xy[7] = dp1*dq4
 d2xy[8] = dp2*dq4
 d2xy[9] = dp3*dq3
 d2xy[10] = dp4*dq3
 d2xy[11] = dp3*dq4
 d2xy[12] = dp4*dq4
 d2xy[13] = dp3*dq1
 d2xy[14] = dp4*dq1
 d2xy[15] = dp3*dq2
 d2xy[16] = dp4*dq2
 #(d2x, d2y, d2xy)
end


text:
"""
  dataupdat(X::Matrix{T}, Xplus::Matrix{T}) where T<:AbstractFloat

Returns Matrix{Float} of size(X), first u = size(Xup, 1) rows of X are removed and
at the end the updat Xplus is appended.

```jldocstests

julia> a = ones(4,4);

julia> b = zeros(2,4);

julia> dataupdat(a,b)
4×4 Array{Float64,2}:
 1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0
 0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0

```
"""
dataupdat(X::Matrix{T}, Xplus::Matrix{T}) where T<:AbstractFloat =
  vcat(X,Xplus)[1+size(Xplus, 1):end, :]

"""

    momentupdat(M::SymmetricTensor{Float, N}, X::Matrix, Xplus::Matrix)

Returns SymmetricTensor{Float, N} updated moment, given original moment, original data and update
of data - dataup

```jldocstests
julia> x = ones(6, 2);

julia> m = moment(x, 3);

julia> y = 2*ones(2,2);

julia> momentupdat(m, x, y)
SymmetricTensors.SymmetricTensor{Float64,3}(Union{Array{Float64,3}, Void}[[3.33333 3.33333; 3.33333 3.33333]
[3.33333 3.33333; 3.33333 3.33333]], 2, 1, 2, true)

```
"""
function momentupdat(M::SymmetricTensor{T, d}, X::Matrix{T}, Xplus::Matrix{T}) where {T<:AbstractFloat, d}
  tup = size(Xplus,1)
  if tup == 0
    return M
  else
    return M + tup/size(X, 1)*(moment(Xplus, d, M.bls) - moment(X[1:tup,:], d, M.bls))
  end
end

"""

  momentupdat(M::Vector{SymmetricTensor{T}}, X::Matrix{T}, Xplus::Matrix{T})

  Returns Vector{SymmetricTensor} of updated moments

"""
momentupdat(M::Vector{SymmetricTensor{T}}, X::Matrix{T}, Xplus::Matrix{T}) where T <: AbstractFloat =
    [momentupdat(M[i], X, Xplus) for i in 1:length(M)]

"""
  momentarray(X::Matrix{Float}, m::Int, b::Int)

Returns an array of Symmetric Tensors of moments given data and maximum moment order - d
"""
momentarray(X::Matrix{T}, d::Int = 4, b::Int = 4) where T <: AbstractFloat =
    [moment(X, i, b) for i in 1:d]

"""

  moms2cums!(M::Vector{SymmetricTensor})

Changes vector of Symmetric Tensors of moments to vector of Symmetric Tensors of cumulants
```jldocstests

julia> m = momentarray(ones(20,3), 3);

julia> moms2cums!(m)

julia> m[3]

SymmetricTensors.SymmetricTensor{Float64,3}(Union{Array{Float64,3}, Void}[[0.0 0.0; 0.0 0.0]
[0.0 0.0; 0.0 0.0] #undef; #undef #undef]
```

"""
function moms2cums!(M::Vector{SymmetricTensor{T}}) where T <: AbstractFloat
  d = length(M)
  for i in 1:d
    f(σ::Int) = outerprodcum(i, σ, M...; exclpartlen = 0)
    prods = pmap(f, [(2:d)...])
    for k in 2:d
      @inbounds M[i] -= prods[k-1]
    end
  end
end

"""

  cums2moms(cum::Vector{SymmetricTensor})

Returns vector of Symmetric Tensors of moments given vector of Symmetric Tensors
of cumulants
"""
function cums2moms(cum::Vector{SymmetricTensor{T}}) where T <: AbstractFloat
  m = length(cum)
  Mvec = Array{SymmetricTensor{T}}(undef, m)
  for i in 1:m
    f(σ::Int) = outerprodcum(i, σ, cum...; exclpartlen = 0)
    @inbounds Mvec[i] = cum[i]
    prods = pmap(f, [(2:m)...])
    for k in 2:m
      @inbounds Mvec[i] += prods[k-1]
    end
  end
  Mvec
end

"""
  mutable struct DataMoments{T <: AbstractFloat}

structure that stores data (X), array of moments (M) and parameters:
d - maximal moment order
b - a size of a block
"""
mutable struct DataMoments{T <: AbstractFloat}
    X::Matrix{T}
    d::Int
    b::Int
    M::Vector{SymmetricTensor{T}}
end

"""
  DataMoments(X, d, b)

a constructor, claculates an Array of moments given data and parameters
"""
DataMoments(X::Matrix{T}, d::Int, b::Int) where T <: AbstractFloat =
  DataMoments(X, d, b, momentarray(X, d, b))

"""
  function cumulantsupdate!(dm::DataMoments{T}, Xplus::Matrix{T}) where T <: AbstractFloat

updates the DataMoments structure in a sliding window, given an updat Xplus
and returns cumululants of updated data.

```jldocstests

julia> x = ones(10,2);

julia> s = DataMoments(x, 4, 2);

julia> y = zeros(4,2);

julia> cumulantsupdate!(s,y)[4]
SymmetricTensors{Float64,4}(Union{Array{Float64,4}, Void}[[0.0064 0.0064; 0.0064 0.0064]

[0.0064 0.0064; 0.0064 0.0064]

[0.0064 0.0064; 0.0064 0.0064]

[0.0064 0.0064; 0.0064 0.0064]], 2, 1, 2, true)

```
"""
function cumulantsupdate!(dm::DataMoments{T}, Xplus::Matrix{T}) where T <: AbstractFloat
  Mup = momentupdat(dm.M, dm.X, Xplus)
  Xup = dataupdat(dm.X, Xplus)
  dm.M = Mup
  dm.X = Xup
  c = copy(Mup)
  moms2cums!(c)
  c
end

"""
  savedm(dm::DataMoments, str::String)

saves a DataMoment structure at a given direcory
"""
savedm(dm::DataMoments, dir::String) = save(dir, Dict{String, DataMoments}("dm" => dm))

"""
  loaddm(str::String)

loads a DataMoment structure from a given direcory
"""
loaddm(dir::String) = load(dir)["dm"]


text:
function check_bin_1()
    nb_var_1 = 0
    nb_files = 0
    for file in readdir("solutions_mp_data")
        nb_files +=1
        lines = readlines(joinpath("solutions_mp_data", file))
        for line in lines
            var, value = split(line, ';')
            if var == "Shunt_9"
                if parse(Float64, value) == 1.0
                    nb_var_1 +=1
                else
                    println(file)
                end
            end
        end
    end
    return nb_files, nb_var_1
end

 nb_files, nb_var_1 = check_bin_1()
println("nb files = $nb_files")
println("nb var 1 = $nb_var_1")


text:
export BatchedSampler

mutable struct BatchedObsDMSampler{BN, P, S, Sc, Sv, Pv, LC, Lv, Pd} <: AbstractIterativeSampler
    observables::Dict

    bnet::BN
    hilb_do::P

    sampler::S
    sampler_cache::Sc
    samples::Sv

    ψvals::Pv

    accum::LC

    local_vals::Lv

    results::Dict

    parallel_cache::Pd
end

function BatchedObsDMSampler(bnet,
                        sampl,
                        hilb_doubled;
                        batch_sz=2^4,
                        local_batch_sz=batch_sz,
                        par_type=automatic_parallel_type())
    hilb_ph        = physical(hilb_doubled)

    par_cache      = parallel_execution_cache(par_type)

    v              = state(hilb_ph, bnet)
    sampler_cache  = cache(sampl, hilb_ph, bnet, par_cache)

    ch_len         = chain_length(sampl, sampler_cache)

    samples        = state(hilb_ph, bnet, ch_len)
    ψvals          = similar(trainable_first(bnet), out_type(bnet), 1, batch_sz, ch_len)

    local_acc      = AccumulatorObsScalar(bnet, hilb_doubled, v, local_batch_sz)
    Llocal_vals    = similar(ψvals, size(ψvals)[2:end]...)

    nq = BatchedObsDMSampler(Dict(), bnet, hilb_doubled,
            sampl, sampler_cache, samples,
            ψvals,
            local_acc, Llocal_vals,
            Dict(),
            par_cache)

    return nq
end

function add_observable!(is::BatchedObsDMSampler, name::String, obs::AbsLinearOperator)
    is.observables[name] = obs
end

function compute_observables(is::BatchedObsDMSampler)
    ch_len         = chain_length(is.sampler, is.sampler_cache)
    batch_sz       = size(is.local_vals, 1)

    # skip if no observables
    isempty(is.observables) && return nothing

    # Sample phase
    _sample_state!(is)

    # Compute logψ and ∇logψ
    logψ!(is.ψvals, is.bnet, is.samples)

    # Compute LdagL
    for (name, Ô) = is.observables
        is.results[name] = _local_obs_eval(is, Ô)
    end

    return is.results
end

function compute_observable(is::BatchedObsDMSampler, Ô)
    # Sample phase
    _sample_state!(is)

    # Compute logψ and ∇logψ
    logψ!(is.ψvals, is.bnet, is.samples)

    return _local_obs_eval(is, Ô)
end

function _local_obs_eval(is::BatchedObsDMSampler, Ô)
    ch_len         = chain_length(is.sampler, is.sampler_cache)
    batch_sz       = size(is.local_vals, 1)

    for i=1:ch_len
        for j = 1:batch_sz
            σv = state_uview(is.samples, j, i)
            init!(is.accum, σv, is.ψvals[1,j,i])
            accumulate_connections!(is.accum, Ô, σv)
            O_loc = NeuralQuantum.finalize!(is.accum)
            is.local_vals[j, i] = O_loc
        end
    end

    return stat_analysis(is.local_vals, is.parallel_cache)
end

Base.show(io::IO, is::BatchedObsDMSampler) = print(io,
    "BatchedObsDMSampler for :"*
    "\n\tnet\t\t: $(is.bnet)"*
    "\n\tsampler\t: $(is.sampler)")


text:
module CryptoRandom

@windows_only cryptprov = Array(Ptr{Void}, 1)

function __init__()
    @windows_only begin
        retval = ccall((:CryptAcquireContextA, "Advapi32"), Cint,
                          (Ptr{Ptr{Void}}, Ptr{UInt8}, Ptr{UInt8},
                           UInt32, UInt32),
                          cryptprov, C_NULL, C_NULL, 1 ,0) # 1 = PROV_RSA_FULL

        if retval == 0
            error("Failed to acquire crypt context.")
        end
    end

end

function random(len)
    str = Array(UInt8, len)
    retval = ccall((:CryptGenRandom, "Advapi32"), Cint,
                   (Ptr{Void}, Cuint, Ptr{UInt8}),
                   cryptprov[1], len, str)
    if retval == 0
        error("Failed to get random data.")
    else
        str
    end
end

end # module

max_stars_repo_name: LachlanGunn/CryptoRandom

max_issues_repo_name: LachlanGunn/CryptoRandom

max_forks_repo_name: LachlanGunn/CryptoRandom


text:
using ProjAssistant
@quickactivate

# ---------------------------------------------------------------------------------
@time begin
    import Chemostat_EColi
    const ChE = Chemostat_EColi

    import Chemostat_Kayser2005
    const ChK = Chemostat_Kayser2005
    const KiJR = ChK.iJR904
    const Kd = ChK.KayserData

    import Chemostat_Nanchen2006
    const ChN = Chemostat_Nanchen2006
    const NiJR = ChN.iJR904
    const Nd = ChN.NanchenData

    import Chemostat_Folsom2014
    const ChF = Chemostat_Folsom2014
    const FiJR = ChF.iJR904
    const Fd = ChF.FolsomData

    import ChemostatPlots
    const ChP = ChemostatPlots

    import Chemostat
    import Chemostat.LP.MathProgBase
    const Ch = Chemostat
    const ChU = Ch.Utils
    const ChLP = Ch.LP

    import ImgTools
    import ImgTools: make_grid, WHITE_PIX

    using ProgressMeter
    using Base.Threads
    using SparseArrays
    using MagicDicts
    using Serialization
    using Statistics
    using Plots
    using Plots.Measures

    using LaTeXStrings

    import GR
    !isinteractive() && GR.inline("png")

end

## ---------------------------------------------------------------------------------
ider_colors = Dict(
    "GLC" => :black,
    "CO2" => :blue,
    "O2" => :red,
    "AC" => :pink,
    "NH4" => :orange,
    "FORM" => :purple,
    "SUCC" => :yellow,
)

exp_colors = let
    EXPS = 1:20
    len = length(EXPS)
    colors = Plots.distinguishable_colors(len)
    Dict(ider => color for (ider, color) in zip(EXPS, colors))
end

subs_markers = Dict(
    "glycolysis" => :circle,
    "krebs" => :diamond,
    "pentose phosphate" => :utriangle,
    "others" => :star,
    "glyoxylate shunt" => :square,
    "glc tranport" => :hex
)

source_markers = Dict(
    Kd => :square,
    Fd => :circle,
    Nd => :utriangle,
)

## ------------------------------------------------------
function _textbf(strs...)
   strs = replace.(string.(strs), " " => "~")
   latexstring(string("\\mathbf{", strs..., "}"))
end

## ---------------------------------------------------------------------------------
# figure 8
include("1.1_pol_box_size.jl")
include("1.2_Xexp_Xmax_study.jl")
include("1.3_residual_Glc_study.jl")

let
    # panel 1
    # TODO: Add colorbar_title
    Dlim = [0.01, 0.5]
    cgD_Xlim = [0.0, 40.0]
    bins = 60
    panel1 = plot_pol_box_size(
        ChN.iJR904, ChN.NanchenData;
        Dlim, cgD_Xlim, bins
    )

    # panel 2
    panel2 = plot_Xexp_Xmax_study()

    # panel 3
    panel3 = plot_residual_Glc_study()

    for (title, p) in [
            ("A", panel1),
            ("B", panel2),
            ("C", panel3)
        ]
        plot!(p;
            title,
            titlefont = 22,
            legendfont = 18,
            guidefont = 22,
            xtickfont = 21,
            ytickfont = 21,
            thickness_scaling = 1.6,
            size = (1320, 940)
        )
    end

    sfig(ChE, [panel1, panel2, panel3],
        "figure_8", ".png";
        layout = (1, 3)
    )

end

## ---------------------------------------------------------------------------------
# figure 9
include("1.4_inner_flx_corr.jl")
include("1.5_exchs_flx_correlations.jl")

let

    METHODS = [
        :FBA_MAX_Z_MIN_COST,
        :FBA_Z_FIX_MAX_VG_MIN_COST,
        :FBA_Z_FIX_MIN_VG_MIN_COST,
        :FBA_Z_FIX_MAX_ATP_MAX_COST,
        :ME_MAX_POL
    ]

    METHODS_LABELS = [
        "FBA max. z",
        "FBA min. ug",
        "FBA max. ug",
        "FBA max. ATP",
        "ME",
    ]

    exch_ps = plot_exch_corrs(METHODS, METHODS_LABELS)
    inner_ps = plot_inner_plots(METHODS, METHODS_LABELS)

    ps = []
    push!(ps, exch_ps...)
    push!(ps, inner_ps[0.05]...)
    push!(ps, inner_ps[0.1]...)
    push!(ps, inner_ps[0.2]...)
    push!(ps, inner_ps[0.3]...)
    push!(ps, inner_ps[0.4]...)
    push!(ps, inner_ps[:tot]...)

    params = (;
       titlefont = 26,
       axisfont = 26,
       guidefont = 26,
       xtickfont = 22,
       ytickfont = 22,
       legendfont = 12,
       thickness_scaling = 1.6,
       size = (1220, 940)
    )

    for p in ps
        plot!(p; params...)
    end

    # add margin
    for ps in [exch_ps, values(inner_ps)...]
        plot!(first(ps); leftmargin = 5mm)
        plot!(last(ps); rightmargin = 7mm)
    end

    sfig(ChE, ps,
        "figure_9", ".png";
        layout = (7, 5)
    )

end;

## ---------------------------------------------------------------------------------
# figure 11
include("1.6_var_study_and_marginals.jl")

let
    corr_av_p, corr_va_p = plot_ME_methods_comparizon()
    pGLC, pAC, pz = plot_me_marginals()

    size0 = (1220, 940)
    for p in [pGLC, pAC, pz]
        plot!(p;
            titlefont = 26,
            axisfont = 26,
            xtickfont = 22,
            ytickfont = 22,
            legendfont = 22,
            guidefont = 26,
            thickness_scaling = 1.6,
            size = size0
        )
    end

    margs_file = sfig(ChE, [pGLC, pAC, pz],
        "margs.png";
        layout = (3, 1)
    )

    for p in [corr_av_p, corr_va_p]
        plot!(p;
            titlefont = 26,
            axisfont = 26,
            guidefont = 26,
            xtickfont = 22,
            ytickfont = 22,
            legendfont = 22,
            thickness_scaling = 1.6,
            size = 1.5 .* size0
        )
    end

    corrs_file = sfig(ChE, [corr_av_p, corr_va_p],
        "corrs.png";
        layout = (2, 1)
    )

    sfig(ChE, lfig.([margs_file, corrs_file]),
        "figure_11", ".png";
        layout = (1, 2)
    )

    rm.([margs_file, corrs_file]; force = true)

end;


text:
function Base.show(io::IO, ::MIME"text/plain", ok::Ok)
    printstyled(io, "Try.Ok"; color = :green, bold = true)
    print(io, ": ")
    show(io, MIME"text/plain"(), Try.unwrap(ok))
end

function Base.show(io::IO, ::MIME"text/plain", err::Err)
    printstyled(io, "Try.Err"; color = :red, bold = true)
    print(io, ": ")
    ex = Try.unwrap_err(err)
    backtrace = err.backtrace
    if backtrace === nothing
        showerror(io, ex)
    else
        showerror(io, ex, simplify_backtrace(err.backtrace))
    end
end


text:
sampletype(::Dists.DiscreteUnivariateDistribution) = Int
sampletype(::Dists.DiscreteMultivariateDistribution) = Vector{Int}
sampletype(::Dists.ContinuousUnivariateDistribution) = Float64
sampletype(::Dists.ContinuousMultivariateDistribution) = Vector{Float64}


text:
int main()
{
   double a[33];
   double x;
   a[x]++;
   return 0;
}

max_stars_repo_name: lukaszcz/javalette

max_issues_repo_name: lukaszcz/javalette

max_forks_repo_name: lukaszcz/javalette


text:
# utilities for network flow
using LightGraphs: inneighbors, outneighbors

export uniq_flow, flow_path_decomp

"Create dense demand vector for a topology"
function dense_demand(inst::Instance, topo::Topology)
    if length(inst.nodes) == length(topo.nodes)
        # just stop if there are no extra Steiner nodes
        return inst.demand
    end
    if sum(topo.nodes .== topo.nodes') > length(topo.nodes)
        # there must be duplicate nodes, leading to ambiguity!
        throw(ArgumentError("Duplicate nodes in topology!"))
    end

    # allow for additional steiner nodes (demand = 0)
    termidx = termindex(topo.nodes, inst.nodes)
    demand = fill(0.0, length(topo.nodes))
    demand[termidx] = inst.demand
    demand
end

"Compute unique flow on arcs for given topology and demand vector."
function uniq_flow(topo::Topology, demand::Vector{Float64})
    length(topo.nodes) == length(demand) ||
        throw(ArgumentError("Topology and demand incompatible!"))
    is_tree(topo) || throw(ArgumentError("Topology is not a tree!"))
    sum(demand) ≈ 0.0 || throw(ArgumentError("Demand is not balanced!"))

    N = incidence(topo)
    N \ demand
end

function uniq_flow(inst::Instance, topo::Topology)
    uniq_flow(topo, dense_demand(inst, topo))
end

"Reorient the arcs of a tree so that flow goes forward (positive)"
function reorient_fwdflow(topo::Topology, demand::Vector{Float64})
    arcflow = uniq_flow(topo, demand)
    arcs = Arc[]
    for (arc, flow) in zip(topo.arcs, arcflow)
        if flow ≥ 0
            push!(arcs, arc)
        else
            push!(arcs, Arc(arc.head, arc.tail))
        end
    end
    Topology(topo.nodes, arcs)
end

function reorient_fwdflow(inst::Instance, topo::Topology)
    reorient_fwdflow(topo, dense_demand(inst, topo))
end

"Decompose arc flow into paths from sources to sinks."
function flow_path_decomp(topo::Topology, arcflow::Vector{Float64})
    dg = digraph_from_topology(topo)
    nnodes, narcs = length(topo.nodes), length(topo.arcs)
    @assert length(arcflow) == narcs
    @assert all(arcflow .≥ 0.0)

    arcidx = arcindex(topo)
    inarcs = [[arcidx[Arc(u,v)] for u in inneighbors(dg,v)] for v in 1:nnodes]
    outarcs = [[arcidx[Arc(v,w)] for w in outneighbors(dg,v)] for v in 1:nnodes]

    # initialize residual flow
    resflow = copy(arcflow)
    resdem = fill(0.0, nnodes)
    for v = 1:nnodes
        resdem[v] += sum(resflow[inarcs[v]]) - sum(resflow[outarcs[v]])
    end
    @assert sum(resdem) ≈ 0.0

    paths = Vector{Arc}[]
    pathflows = Float64[]

    # iteratively move flow on paths
    while any(nonzero, resdem)
        path = []

        # start from largest remaining source
        current = argmin(resdem)
        pflow = -resdem[current]
        @assert pflow > 0.0

        while resdem[current] <= 0.0
            # pick largest outgoing arc
            maxarc = argmax([resflow[a] for a in outarcs[current]])
            curarc = outarcs[current][maxarc]
            pflow = min(pflow, resflow[curarc])
            @assert pflow > 0.0

            push!(path, topo.arcs[curarc])
            current = path[end].head
        end

        push!(paths, path)
        push!(pathflows, pflow)

        resdem[path[1].tail] += pflow
        for arc in path
            resflow[arcidx[arc]] -= pflow
        end
        resdem[path[end].head] -= pflow
    end

    @assert !any(nonzero, resdem)
    @assert !any(nonzero, resflow)

    paths, pathflows
end


text:
abstract type AbstractAdvection <: AbstractPhysicsComponent end

struct NonLinearAdvection <: AbstractAdvection end

struct ESDGLinearAdvection <: AbstractAdvection end
Base.@kwdef struct ESDGNonLinearAdvection{𝒯} <: AbstractAdvection
    eos::𝒯
end

@inline calc_advective_flux!(flux, ::Nothing, _...) = nothing
@inline calc_advective_flux!(flux, ::Nothing, _...) = nothing
@inline calc_advective_flux!(flux, ::AbstractAdvection, _...) = nothing

@inline function calc_advective_flux!(flux, ::NonLinearAdvection, state, aux, t)
    ρ = state.ρ
    ρu = state.ρu
    ρθ = state.ρθ

    flux.ρ += state.ρu
    flux.ρu += ρu ⊗ ρu / ρ
    flux.ρθ += ρu * ρθ / ρ

    return nothing
end

@inline function calc_flux!(flux, ::ESDGLinearAdvection, state, aux, t)
    ρu = state.ρu
    ρᵣ = aux.ref_state.ρ
    pᵣ = aux.ref_state.p
    ρeᵣ = aux.ref_state.ρe

    flux.ρ += ρu
    #flux.ρu += -0
    flux.ρe += (ρeᵣ + pᵣ) / ρᵣ * ρu

    return nothing
end

@inline function calc_flux!(flux, advection::ESDGNonLinearAdvection, state, aux, t)
    ρ = state.ρ
    ρu = state.ρu
    ρe = state.ρe
    eos = advection.eos

    p = calc_pressure(eos, state, aux)

    flux.ρ += ρu
    flux.ρu += ρu ⊗ ρu / ρ
    flux.ρe += ρu * (ρe + p) / ρ

    return nothing
end

max_stars_repo_name: CliMA/CouplerMachine

max_issues_repo_name: CliMA/CouplerMachine