text:
using HyperNEAT
using Base.Test

# test cppn functions
genome = HyperNEAT.new_genome(
neurons = [
HyperNEAT.new_neuron(innovation = 1, ntype=:input,  func=:id, parameters=[1.0, 0.0]),
HyperNEAT.new_neuron(innovation = 1, ntype=:output, func=:id, parameters=[1.0, 0.0])],
synapses = [
HyperNEAT.new_synapse(innovation = 1, src=1, dest=2, weight=1.0)])
cppn = HyperNEAT.phenotype(genome)

for i in linspace(-10.0, 10.0, 1001)
  @test [i] == HyperNEAT.update_cppn!(cppn, [i])
end

genome = HyperNEAT.new_genome(
neurons = [
HyperNEAT.new_neuron(innovation = 1, ntype=:input,  func=:id, parameters=[1.0, 0.0]),
HyperNEAT.new_neuron(innovation = 1, ntype=:output, func=:id, parameters=[2.0, 0.0])],
synapses = [
HyperNEAT.new_synapse(innovation = 1, src=1, dest=2, weight=1.0)])
cppn = HyperNEAT.phenotype(genome)

for i in linspace(-10.0, 10.0, 1001)
  @test [2.0 * i] == HyperNEAT.update_cppn!(cppn, [i])
end

genome = HyperNEAT.new_genome(
neurons = [
HyperNEAT.new_neuron(innovation = 1, ntype=:input,  func=:id, parameters=[1.0, 0.0]),
HyperNEAT.new_neuron(innovation = 1, ntype=:output, func=:id, parameters=[1.0, 0.0])],
synapses = [
HyperNEAT.new_synapse(innovation = 1, src=1, dest=2, weight=2.0)])
cppn = HyperNEAT.phenotype(genome)
for i in linspace(-10.0, 10.0, 1001)
  @test [2.0 * i] == HyperNEAT.update_cppn!(cppn, [i])
end

genome = HyperNEAT.new_genome(
neurons = [
HyperNEAT.new_neuron(innovation = 1, ntype=:input,  func=:id,  parameters=[1.0, 0.0]),
HyperNEAT.new_neuron(innovation = 1, ntype=:output, func=:sin, parameters=[1.0, 0.0])],
synapses = [
HyperNEAT.new_synapse(innovation = 1, src=1, dest=2, weight=1.0)])
cppn = HyperNEAT.phenotype(genome)
for i in linspace(-10.0, 10.0, 1001)
  @test [sin(i)] == HyperNEAT.update_cppn!(cppn, [i])
end

#
# sin
#
genome = HyperNEAT.new_genome(
neurons = [
HyperNEAT.new_neuron(innovation = 1, ntype=:input,  func=:id,  parameters=[1.0, 0.0]),
HyperNEAT.new_neuron(innovation = 1, ntype=:output, func=:sin, parameters=[1.0, 1.0])],
synapses = [
HyperNEAT.new_synapse(innovation = 1, src=1, dest=2, weight=1.0)])
cppn = HyperNEAT.phenotype(genome)
for i in linspace(-10.0, 10.0, 1001)
  @test [sin(i + 1.0)] == HyperNEAT.update_cppn!(cppn, [i])
end

genome = HyperNEAT.new_genome(
neurons = [
HyperNEAT.new_neuron(innovation = 1, ntype=:input,  func=:id,  parameters=[1.0, 0.0]),
HyperNEAT.new_neuron(innovation = 1, ntype=:output, func=:sin, parameters=[2.0, 1.0])],
synapses = [
HyperNEAT.new_synapse(innovation = 1, src=1, dest=2, weight=1.0)])
cppn = HyperNEAT.phenotype(genome)
for i in linspace(-10.0, 10.0, 1001)
  @test [sin(2.0 * i + 1.0)] == HyperNEAT.update_cppn!(cppn, [i])
end

#
# cos
#
genome = HyperNEAT.new_genome(
neurons = [
HyperNEAT.new_neuron(innovation = 1, ntype=:input,  func=:id,  parameters=[1.0, 0.0]),
HyperNEAT.new_neuron(innovation = 1, ntype=:output, func=:cos, parameters=[1.0, 1.0])],
synapses = [
HyperNEAT.new_synapse(innovation = 1, src=1, dest=2, weight=1.0)])
cppn = HyperNEAT.phenotype(genome)
for i in linspace(-10.0, 10.0, 1001)
  @test [cos(i + 1.0)] == HyperNEAT.update_cppn!(cppn, [i])
end

genome = HyperNEAT.new_genome(
neurons = [
HyperNEAT.new_neuron(innovation = 1, ntype=:input,  func=:id,  parameters=[1.0, 0.0]),
HyperNEAT.new_neuron(innovation = 1, ntype=:output, func=:cos, parameters=[2.0, 1.0])],
synapses = [
HyperNEAT.new_synapse(innovation = 1, src=1, dest=2, weight=1.0)])
cppn = HyperNEAT.phenotype(genome)
for i in linspace(-10.0, 10.0, 1001)
  @test [cos(2.0 * i + 1.0)] == HyperNEAT.update_cppn!(cppn, [i])
end

#
# tanh
#
genome = HyperNEAT.new_genome(
neurons = [
HyperNEAT.new_neuron(innovation = 1, ntype=:input,  func=:id,  parameters=[1.0, 0.0]),
HyperNEAT.new_neuron(innovation = 1, ntype=:output, func=:tanh, parameters=[1.0, 1.0])],
synapses = [
HyperNEAT.new_synapse(innovation = 1, src=1, dest=2, weight=1.0)])
cppn = HyperNEAT.phenotype(genome)
for i in linspace(-10.0, 10.0, 1001)
  @test [tanh(i + 1.0)] == HyperNEAT.update_cppn!(cppn, [i])
end

genome = HyperNEAT.new_genome(
neurons = [
HyperNEAT.new_neuron(innovation = 1, ntype=:input,  func=:id,  parameters=[1.0, 0.0]),
HyperNEAT.new_neuron(innovation = 1, ntype=:output, func=:tanh, parameters=[2.0, 1.0])],
synapses = [
HyperNEAT.new_synapse(innovation = 1, src=1, dest=2, weight=1.0)])
cppn = HyperNEAT.phenotype(genome)
for i in linspace(-10.0, 10.0, 1001)
  @test [tanh(2.0 * i + 1.0)] == HyperNEAT.update_cppn!(cppn, [i])
end

#
# gauss
#
genome = HyperNEAT.new_genome(
neurons = [
HyperNEAT.new_neuron(innovation = 1, ntype=:input,  func=:id,    parameters=[1.0, 0.0]),
HyperNEAT.new_neuron(innovation = 1, ntype=:output, func=:gauss, parameters=[1.0, 1.0])],
synapses = [
HyperNEAT.new_synapse(innovation = 1, src=1, dest=2, weight=1.0)])
cppn = HyperNEAT.phenotype(genome)
for i in linspace(-10.0, 10.0, 1001)
  @test [HyperNEAT.gauss(i,1.0,1.0)] == HyperNEAT.update_cppn!(cppn, [i])
end

genome = HyperNEAT.new_genome(
neurons = [
HyperNEAT.new_neuron(innovation = 1, ntype=:input,  func=:id,    parameters=[1.0, 0.0]),
HyperNEAT.new_neuron(innovation = 1, ntype=:output, func=:gauss, parameters=[1.0, 0.0])],
synapses = [
HyperNEAT.new_synapse(innovation = 1, src=1, dest=2, weight=1.0)])
cppn = HyperNEAT.phenotype(genome)
for i in linspace(-10.0, 10.0, 1001)
  @test [HyperNEAT.gauss(i,0.0,1.0)] == HyperNEAT.update_cppn!(cppn, [i])
end

genome = HyperNEAT.new_genome(
neurons = [
HyperNEAT.new_neuron(innovation = 1, ntype=:input,  func=:id,  parameters=[1.0, 0.0]),
HyperNEAT.new_neuron(innovation = 1, ntype=:output, func=:gauss, parameters=[2.0, 1.0])],
synapses = [
HyperNEAT.new_synapse(innovation = 1, src=1, dest=2, weight=1.0)])
cppn = HyperNEAT.phenotype(genome)
for i in linspace(-10.0, 10.0, 1001)
  @test [HyperNEAT.gauss(i, 1.0,  2.0)] == HyperNEAT.update_cppn!(cppn, [i])
end

#
# combination of functions
#
genome = HyperNEAT.new_genome(
neurons = [
HyperNEAT.new_neuron(innovation = 1, ntype=:input,  func=:id,    parameters=[1.0, 0.0]),
HyperNEAT.new_neuron(innovation = 1, ntype=:output, func=:gauss, parameters=[1.0, 0.0]),
HyperNEAT.new_neuron(innovation = 1, ntype=:hidden, func=:sin,   parameters=[1.0, 0.0])],
synapses = [
HyperNEAT.new_synapse(innovation = 1, src=1, dest=3, weight=1.0),
HyperNEAT.new_synapse(innovation = 1, src=3, dest=2, weight=1.0)])
cppn = HyperNEAT.phenotype(genome)
for i in linspace(-12.0, 12.0, 1001)
  value = HyperNEAT.gauss(sin(i), 0.0, 1.0)
  calc  = HyperNEAT.update_cppn!(cppn, [i])
  @test [value] == calc
end


text:
using Test

function func1(arr)
  @test arr[1] == 1.0
  @test arr[2] == 2.0
  @test arr[3] == 3.0
end

module Containers
  include(joinpath(@__DIR__, "testcommon.jl"))
  @wrapmodule CxxWrap.CxxWrapCore.libjlcxx_containers()

  function __init__()
    @initcxx
  end

  export test_tuple, const_ptr, const_ptr_arg, const_vector, const_matrix
end
using Main.Containers

@testset "$(basename(@__FILE__)[1:end-3])" begin

@test test_tuple() == (1,2.0,3.0f0)

cptr = const_ptr()
@test isbitstype(typeof(cptr))
@test const_ptr_arg(cptr) == (1., 2., 3.)

cv = const_vector()
@test size(cv) == (3,)
@test cv == [1.,2.,3.]

cm = const_matrix()
@test size(cm) == (3,2)
@test cm == [[1.,2.,3.] [4.,5.,6.]]
println("Displaying const matrix")
display(cm)
println()

mm = Containers.mutable_array()
println("Displaying mutable matrix")
display(mm)
println()
mm .= 1.0
@test Containers.check_mutable_array(mm)
Containers.do_embedding_test()

@test Containers.array_return() == ["hello", "world"]
@test Containers.tuple_int_pointer() == (C_NULL,1)

let a1 = [UInt8(3)], a2 = [UInt8(5)]
  @test Containers.uint8_ptr(pointer(a1)) == 3
  @test Containers.uint8_ptr(pointer(a2)) == 5
  @test Containers.uint8_arrayref([pointer(a1), pointer(a2)]) == 8
end

end

ext: jl


text:
using GeometricIntegrators
using ChargedParticleDynamics.GuidingCenter4d.TokamakFastSurface
using GeometricExamples

const Δt     = 1.
const ntime  = 25000
const nx     = 200
const ny     = 200
const nsave  = 50
const nplot  = 5
const run_id = "poincare_2nd_fast_firk_dt1"

include("guiding_center_4d_settings_firk.jl")

tableau_list = get_tableau_list_firk()

pinv = guiding_center_4d_ode_poincare_invariant_2nd(Δt, nx, ny, ntime, nsave)

include("guiding_center_4d_fast_poincare_invariant_2nd.jl")


text:
module Plug # module Bukdu

import ..Deps
import ..Assoc

"""
    ApplicationController
"""
abstract type ApplicationController end
abstract type AbstractPlug end
abstract type AbstractRender end

"""
    Render <: AbstractRender
"""
struct Render <: AbstractRender
    content_type::String
    body::Vector{UInt8}
end

include("plugs/logger.jl")
include("plugs/conn.jl")

struct EventStream <: AbstractRender
    content_type::String
    body::Vector{UInt8}
end

include("plugs/Parsers.jl")
include("plugs/static.jl")
include("plugs/websocket.jl")
include("plugs/server_sent_events.jl")

# pipeline plugs
include("plugs/csrf_protection.jl")
include("plugs/auth.jl")

function plug(::Type{T}; kwargs...) where {T <: AbstractPlug}
end

end # module Bukdu.Plug

export Plug, Conn, ApplicationController, Render, EventStream, plug

import .Plug: Conn, ApplicationController, AbstractPlug, AbstractRender, Render, EventStream, plug

# module Bukdu

max_stars_repo_name: wookay/TestBukdu

max_issues_repo_name: wookay/TestBukdu

max_forks_repo_name: wookay/TestBukdu


text:
#nb # %% A slide [markdown] {"slideshow": {"slide_type": "subslide"}}
# ## Resistive Square Patch
# This example will demonstrate the ability of PSSFSS to accurately model finite
# conductivity of FSS metalization.  It consists of a square finitely conducting
# patch in a square lattice.  It is taken from a paper by Alon S. Barlevy and
# Yahya Rahmat-Samii,
# "Fundamental Constraints on the Electrical Characteristics of Frequency Selective
# Surfaces", **Electromagnetics**, vol. 17, 1997, pp. 41-68. This particular example
# is from Section 3.2, Figures 7 and 8.  We will compare PSSFSS results to those digitized
# from the cited figures.
#
# We start by defining a function that creates a patch of the desired sheet resistance:

using Plots, PSSFSS
patch(R) = rectstrip(Nx=10, Ny=10, Px=1, Py=1, Lx=0.5, Ly=0.5, units=cm, Rsheet=R)
plot(patch(0), unitcell=true)
#md savefig("resistive1.png"); nothing  # hide
#-
#md # ![](resistive1.png)

# The patches measure 0.5 cm on a side and lie in a square lattice of period 1 cm.
# Now we perform the analysis, looping over the desired values of sheet resistance.

steering = (ϕ=0, θ=0)
flist = 1:0.5:60
Rs = [0, 10, 30, 100]
calculated = zeros(length(flist), length(Rs)) # preallocate storage
outputs = @outputs s11mag(v,v)
for (i,R) in pairs(Rs)
    strata = [Layer(), patch(R), Layer()]
    results = analyze(strata, flist, steering, showprogress=false,
                      logfile=devnull, resultfile=devnull)
    calculated[:,i] = extract_result(results, outputs)
end

# Looping over the four sheet resistance values, each evaluated at 119 frequencies
# required approximately 20 seconds on my machine.
#
# We plot the results, including those digitized from the paper for comparison:

using DelimitedFiles
markers = (:diamond, :utriangle, :square, :xcross)
colors = (:blue, :red, :green, :black)
p = plot(xlim=(-0.01,60.01), xtick = 0:10:60, ylim=(-0.01,1.01), ytick=0:0.1:1,
         xlabel="Frequency (GHz)", ylabel="Reflection Coefficient Magnitude",
         title = "Resistive Square Patch",
         legend=:topright)
for (i,R) in pairs(Rs)
    scatter!(p, flist, calculated[:,i], label="PSSFSS $R Ω", ms=2, shape=markers[i], color=colors[i])
    data = readdlm("../src/assets/barlevy_patch_digitized_$(R)ohm.csv", ',')
    plot!(p, data[:,1], data[:,2], label="Barlevy $R Ω", ls=:dash, color=colors[i])
end
p
#md savefig("resistive2.png"); nothing  # hide
#-
#md # ![](resistive2.png)

# ### Conclusion
# PSSFSS results are indistinguishable from those reported in the cited paper.
#


text:
using Mea.Transport
using Mea.Periodize

using Base.Test

@testset "Testing Transport.jl" begin

    #setup

        cutoff_test = 200.0
        beta_ = 12.0
        sigmadc_good = 0.7637186262560831
        l11_good = sigmadc_good/beta_

    # cutoff_test = 3.0
    # beta_ = 12.0
    # sigmadc_good = 0.48994580280943345
    # l11_good = sigmadc_good/beta_

    @testset "test dfdw" begin
        @test isapprox(Transport.dfdw(12.0, 0.2), -0.91506 , atol=1.0e-4)
        @test isapprox(Transport.dfdw(12.0, -0.1), -2.13473, atol=1.0e-4)
    end

    @testset "test coefstrans" begin
        modelvec=Periodize.buildmodelvec("./data/self_ctow0.dat", "./data/statsparams0.json")
        transport = Transport.coefstrans(modelvec, beta_, cutoff=cutoff_test, fout_name="dos_green.dat", maxevals=80000, fctper="make_akw2green")
        transport_cum = Transport.coefstrans(modelvec, beta_, cutoff=cutoff_test, fout_name="dos_cum.dat", maxevals=80000, fctper="make_akw2cum")
        transport_trace = Transport.coefstrans(modelvec, beta_, cutoff=cutoff_test, fout_name="dos_trace.dat", maxevals=80000, fctper="make_akw2trace")
        transport_cuba = Transport.coefstrans(modelvec, beta_, cutoff=cutoff_test, fout_name="dos_cuba.dat", maxevals=80000, libintegrator="cuba", fctper="make_akw2green")
        println("transport_cuba = ", transport_cuba)
        println("transport = ", transport)
        println("transport_trace = ", transport_trace)
        println("transport_cum = ", transport_cum)
        @test isapprox(transport["n"], 0.495, atol=1e-2, rtol=1e-2)
        @test isapprox(transport["sigmadc"], sigmadc_good, atol=1.0e-5, rtol=1.0e-5)
        @test isapprox(transport_cuba["sigmadc"], sigmadc_good, atol=1.0e-3, rtol=1.0e-3)
        @test isapprox(transport["l11"], sigmadc_good/beta_, atol=1.0e-5, rtol=1.0e-5)
        @test isapprox(transport_cuba["l11"], sigmadc_good/beta_, atol=1.0e-3, rtol=1.0e-3)
        # @test isapprox(transport["l22"], )
        # @test isapprox(transport["l22"],)
        # @test isapprox(transport["seebeck"],)
    end

end  # testset


text:
"""
    $TYPEDEF

Range and theta definition on `(:Euclid, :Circular)` manifold.
"""
struct Polar <: IIF.InferenceVariable
  dims::Int
  manifolds::Tuple{Symbol,Symbol}
  Polar() = new(2,(:Euclid,:Circular))
end

"""
    $TYPEDEF

Prior belief on any Polar related variable.
"""
mutable struct PriorPolar{T1<:IIF.SamplableBelief, T2<:IIF.SamplableBelief} <: IIF.AbstractPrior
  Zrange::T1
  Zangle::T2
  # empty constructor
  PriorPolar{T1,T2}() where {T1,T2} = new{T1,T2}()
  # regular constructor
  PriorPolar{T1,T2}(zr::T1, za::T2) where {T1 <: IIF.SamplableBelief,T2 <: IIF.SamplableBelief} = new{T1,T2}(zr,za)
end
# convenience and default object helper
PriorPolar(zr::T1=Normal(1,1), za::T2=Normal(0,0.1)) where {T1 <: IIF.SamplableBelief,T2 <: IIF.SamplableBelief} = PriorPolar{T1,T2}(zr,za)

function getSample(pp2i::PriorPolar, N::Int=1)
  sps = zeros(2,N)
  sps[1,:] = rand(pp2i.Zrange,N);
  sps[2,:] = rand(pp2i.Zangle,N);
  return (sps, )
end

"""
    $TYPEDEF

Linear offset factor of `IIF.SamplableBelief` between two `Polar` variables.
"""
mutable struct PolarPolar{T1<:IIF.SamplableBelief, T2<:IIF.SamplableBelief} <: IIF.AbstractRelativeRoots
  Zrange::T1
  Zangle::T2
  # empty constructor
  PolarPolar{T1,T2}() where {T1,T2} = new{T1,T2}()
  # regular constructor
  PolarPolar{T1,T2}(zr::T1, za::T2) where {T1 <: IIF.SamplableBelief,T2 <: IIF.SamplableBelief} = new{T1,T2}(zr,za)
end
# convenience and default object helper
PolarPolar(zr::T1=Normal(1,1), za::T2=Normal(0,0.1)) where {T1 <: IIF.SamplableBelief,T2 <: IIF.SamplableBelief} = PolarPolar{T1,T2}(zr,za)

function getSample(pp2i::PolarPolar, N::Int=1)
  sps = zeros(2,N)
  sps[1,:] = rand(pp2i.Zrange,N);
  sps[2,:] = rand(pp2i.Zangle,N);
  return (sps, )
end

function (pp::PolarPolar)(res::Array{Float64},
                          userdata::FactorMetadata,
                          idx::Int,
                          meas::Tuple{Array{Float64,2}},
                          p1::Array{Float64,2},
                          p2::Array{Float64,2} )
  #
  @inbounds res[1:2] = meas[1][1:2,idx] - (p2[1:2,idx] - p1[1:2,idx])
  nothing
end


text:
export samplespace

# each sample space supports a different type of data matrix (dense, sparse, etc.)

type DenseSampleSpace{T}
    data::Array{T,2}
    N::Int64
end
function samplespace{T}(X::Array{T,2}, nsamples::Int)
    DenseSampleSpace(Array(T, size(X)[1], 2*nsamples*size(X)[2]), 0)
end
addsample!(s::DenseSampleSpace, x) = (s.N += 1; s.data[:,s.N] = x)
reset!(s::DenseSampleSpace) = (s.N = 0)
data(s::DenseSampleSpace) = s.data[:,1:s.N]

type SparseSampleSpace{Tv,Ti}
    colptr::Array{Ti,1}
    colptrLen::Ti
    rowval::Array{Ti,1}
    rowvalLen::Ti
    nzval::Array{Tv,1}
    nzvalLen::Ti
    nrows::Ti
end
function samplespace{Tv,Ti}(X::SparseMatrixCSC{Tv,Ti}, nsamples::Int, maxDensity=1.0)
    s = SparseSampleSpace(
        Array(Ti, 2*nsamples*size(X)[2]+1), one(Ti),
        Array(Ti, ceil(Int, maxDensity*size(X)[1]*2*nsamples*size(X)[2]+1)), zero(Ti),
        Array(Tv, ceil(Int, maxDensity*size(X)[1]*2*nsamples*size(X)[2]+1)), zero(Ti), convert(Ti, size(X)[1])
    )
    s.colptr[1] = 1
    s
end
function reset!(s::SparseSampleSpace)
    s.colptrLen = 1
    s.colptr[1] = 1
    s.rowvalLen = 0
    s.nzvalLen = 0
end
function addsample!(s::SparseSampleSpace, x)
    s.colptrLen += 1
    s.colptr[s.colptrLen] = s.colptr[s.colptrLen-1]
    for i in 1:length(x)
        if x[i] != 0
            s.colptr[s.colptrLen] += 1
            s.rowvalLen += 1
            s.rowval[s.rowvalLen] = i
            s.nzvalLen += 1
            s.nzval[s.nzvalLen] = x[i]
        end
    end
end
function data(s::SparseSampleSpace)
    SparseMatrixCSC(
        s.nrows, s.colptrLen-1,
        s.colptr[1:s.colptrLen],
        s.rowval[1:s.rowvalLen],
        s.nzval[1:s.nzvalLen]
    )
end


text:
"""
A series of `Solution`s. One can `push`, `unshift`, `pop`, and `shift` new/old `Solution`s
to/from it.
"""
type Branch
	solutions::Vector{Solution}	# ordered
	parent # ::Project
	Branch(V::Vector{Solution}) = new(V, nothing)
	Branch(P) = Branch(Solution[], P) #TODO runtime typeassert bco circ type?
	Branch(S,P) = new(S,P)
end
@relay(Branch, :solutions)

Base.show(io::Base.IO, P::Branch) = write(io, "***BRANCH***")

#TODO IMPORTANT: handle case of activeSolution::Solution

function Base.push!(B::Branch, sols...)
	P = B.parent
	for s in sols
		push!(B.solutions, s)
		mbObserve.notify(P, :pushSolution, B, s)
	end
	return P
end

function Base.unshift!(B::Branch, sols...)
	P = B.parent
	for s in sols
		unshift!(B.solutions, s)
		mbObserve.notify(P, :unshiftSolution, B, s)
	end
	return P
end

function Base.pop!(B::Branch)
	P = B.parent
	s = pop!(B.solutions)
	s == P.activeSolution && (P.activeSolution=s.data)
	if length(B) == 1
		deleteat!(P, findfirst(P, B))
		B[1] == P.activeSolution && setActiveSolution(P, B[1].data)
		mbObserve.notify(P, :delBranch, B)
	else
		mbObserve.notify(P, :popSolution, B, s)
	end
	return s
end

function Base.shift!(B::Branch)
	P = B.parent
	s = shift!(B.solutions)
	s == P.activeSolution && (P.activeSolution=s.data)
	if length(B.solutions) == 1
		deleteat!(P, findfirst(P, B))
		B[1] == P.activeSolution && setActiveSolution(P, B[1].data)
		mbObserve.notify(P, :delBranch, B)
	else
		mbObserve.notify(P, :shiftSolution, B, s)
	end
	return s
end


text:
function regulateall!(system::CVSystem,n::Int64,terms::Vector{Int64})
    # changes to HR, ventricular contractility
    system.cns.H[n+2] = ((1 - system.solverparams.h/system.cns.tauH)*
        system.cns.H[n+1]/minTos + system.solverparams.h/system.cns.tauH*
        (system.cns.alphaH*system.cns.ns[system.solverparams.numbeats+1]-
        system.cns.betaH*system.cns.np[system.solverparams.numbeats+1]+
        system.cns.gammaH))*minTos;
    system.cns.Emaxlv[n+2] = (1 - system.solverparams.h/system.cns.tauEmax)*
        system.cns.Emaxlv[n+1] + system.solverparams.h/system.cns.tauEmax*(
        system.cns.alphalv*system.cns.ns[system.solverparams.numbeats+1]+
        system.cns.gammalv);
    system.cns.Emaxrv[n+2] = (1 - system.solverparams.h/system.cns.tauEmax)*
        system.cns.Emaxrv[n+1] + system.solverparams.h/system.cns.tauEmax*(
        system.cns.alpharv*system.cns.ns[system.solverparams.numbeats+1]+
        system.cns.gammarv);

    # changes to peripheral impedance
    system.cns.R2L[n+2] = (1 - system.solverparams.h/system.cns.tauR)*
        system.cns.R2L[n+1] + system.solverparams.h/system.cns.tauR*(
        system.cns.alphaR2L*system.cns.ns[system.solverparams.numbeats+1]+
        system.cns.gammaR2L);
    system.cns.R3L[n+2] = (1 - system.solverparams.h/system.cns.tauR)*
        system.cns.R3L[n+1] + system.solverparams.h/system.cns.tauR*(
        system.cns.alphaR3L*system.cns.ns[system.solverparams.numbeats+1]+
        system.cns.gammaR3L);
    system.cns.R2U[n+2] = (1 - system.solverparams.h/system.cns.tauR)*
        system.cns.R2U[n+1] + system.solverparams.h/system.cns.tauR*(
        system.cns.alphaR2U*system.cns.ns[system.solverparams.numbeats+1]+
        system.cns.gammaR2U);
    system.cns.R3U[n+2] = (1 - system.solverparams.h/system.cns.tauR)*
        system.cns.R3U[n+1] + system.solverparams.h/system.cns.tauR*(
        system.cns.alphaR3U*system.cns.ns[system.solverparams.numbeats+1]+
        system.cns.gammaR3U);
    system.cns.C4L[n+2] = (1 - system.solverparams.h/system.cns.tauC)*
        system.cns.C4L[n+1] + system.solverparams.h/system.cns.tauC*(
        -system.cns.alphaC4L*system.cns.ns[system.solverparams.numbeats+1]+
        system.cns.gammaC4L);
    system.cns.C5L[n+2] = (1 - system.solverparams.h/system.cns.tauC)*
        system.cns.C5L[n+1] + system.solverparams.h/system.cns.tauC*(
        -system.cns.alphaC5L*system.cns.ns[system.solverparams.numbeats+1]+
        system.cns.gammaC5L);
    system.cns.C4U[n+2] = (1 - system.solverparams.h/system.cns.tauC)*
        system.cns.C4U[n+1] + system.solverparams.h/system.cns.tauC*(
        -system.cns.alphaC4U*system.cns.ns[system.solverparams.numbeats+1]+
        system.cns.gammaC4U);
    system.cns.C5U[n+2] = (1 - system.solverparams.h/system.cns.tauC)*
        system.cns.C5U[n+1] + system.solverparams.h/system.cns.tauC*(
        -system.cns.alphaC5U*system.cns.ns[system.solverparams.numbeats+1]+
        system.cns.gammaC5U);
    system.cns.V4L[n+2] = (1 - system.solverparams.h/system.cns.tauC)*
        system.cns.V4L[n+1] + system.solverparams.h/system.cns.tauC*(
        -system.cns.alphaV4L*system.cns.ns[system.solverparams.numbeats+1]+
        system.cns.gammaV4L);
    system.cns.V5L[n+2] = (1 - system.solverparams.h/system.cns.tauC)*
        system.cns.V5L[n+1] + system.solverparams.h/system.cns.tauC*(
        -system.cns.alphaV5L*system.cns.ns[system.solverparams.numbeats+1]+
        system.cns.gammaV5L);
    system.cns.V4U[n+2] = (1 - system.solverparams.h/system.cns.tauC)*
        system.cns.V4U[n+1] + system.solverparams.h/system.cns.tauC*(
        -system.cns.alphaV4U*system.cns.ns[system.solverparams.numbeats+1]+
        system.cns.gammaV4U);
    system.cns.V5U[n+2] = (1 - system.solverparams.h/system.cns.tauC)*
        system.cns.V5U[n+1] + system.solverparams.h/system.cns.tauC*(
        -system.cns.alphaV5U*system.cns.ns[system.solverparams.numbeats+1]+
        system.cns.gammaV5U);

    CVModule.regulateperiphery!(system,n,terms);

    # check for start of new cardiac cycle
    oldnumbeats = system.solverparams.numbeats;
    CVModule.setnumbeats!(system,n);
    if oldnumbeats < system.solverparams.numbeats
        # baroreflex pressure
        CVModule.reflexpressure!(system,n);
        # SNS/PSNS activations
        CVModule.cnsactivations!(system,n);
        CVModule.regulateheart!(system,n);
        CVModule.correctvolume!(system,n);
        # update time discretization
        if system.solverparams.numbeats < system.solverparams.numbeatstotal
            CVModule.updatediscretization!(system);
        end
    end
end

max_stars_repo_name: dcanuto/CardioModeling

max_issues_repo_name: dcanuto/CardioModeling

max_forks_repo_name: dcanuto/CardioModeling


text:
export compute_interactions!

"""
    compute_interations!(lj; period_x, period_y)

Computes a Lennard-Jones force for particles stored under `lj.particles` given
neighbors stored in the cell list `lj.cell_list`.  If `period_x > 0` and
`period_y > 0`, then apply periodic boundary conditions to the pair interaction.
"""
function compute_interactions!(lj::LennardJones; period_x::Float64 = -1.0, period_y::Float64 = -1.0)
    @use_threads lj.multithreaded for particle in lj.particles
        x, y = particle.x, particle.y
        i = trunc(Int64, x / lj.cell_list.cell_spacing_x)
        j = trunc(Int64, y / lj.cell_list.cell_spacing_y)
        f_x, f_y = 0.0, 0.0
        @inbounds for dj = -1 : 1, di = -1 : 1
            idi = mod(i + di, lj.cell_list.num_cells_x) + 1
            jdj = mod(j + dj, lj.cell_list.num_cells_y) + 1
            pid = lj.cell_list.start_pid[idi, jdj]
            while pid > 0
                neighbor = lj.cell_list.particles[pid]
                Δx = wrap_displacement(x - neighbor.x; period = period_x)
                Δy = wrap_displacement(y - neighbor.y; period = period_y)
                Δr² = Δx^2 + Δy^2

                σ² = (lj.σ < 0.0 ? particle.R + neighbor.R : lj.σ)^2
                cutoff² = (lj.cutoff < 0.0 ? 1.2599210498948732 * σ² : lj.cutoff^2)
                if 0.0 < Δr² < cutoff²
                    inv² = σ² / Δr²
                    inv⁶ = inv² * inv² * inv²
                    coef = lj.ϵ * (48.0 * inv⁶ - 24.0) * inv⁶ / Δr²

                    f_x += (_f_x = coef * Δx)
                    f_y += (_f_y = coef * Δy)
                    if lj.use_newton_3rd
                        neighbor.f_x -= _f_x
                        neighbor.f_y -= _f_y
                    end
                end
                pid = lj.cell_list.next_pid[pid]
            end
        end
        particle.f_x += f_x
        particle.f_y += f_y
    end
end

"""
    compute_interations!(hbond; period_x, period_y)

Compute a harmonic force between pairs of particles stored under
`hb.particle_pairs`.  Apply periodic boundary conditions if `period_x > 0` or
`period_y > 0`.
"""
function compute_interactions!(hbond::HarmonicBond; period_x::Float64 = -1.0, period_y::Float64 = -1.0)
    @use_threads hbond.multithreaded for (particle_1, particle_2) in hbond.particle_pairs
        Δx = wrap_displacement(particle_1.x - particle_2.x; period = period_x)
        Δy = wrap_displacement(particle_1.y - particle_2.y; period = period_y)
        Δr = sqrt(Δx^2 + Δy^2)

        # this check isn't really necessary unless user starts the particles in the same place
        if Δr > 0.0
            coef = -hbond.k_bond * (1 - hbond.l_rest / Δr)
            f_x = coef * Δx
            f_y = coef * Δy

            particle_1.f_x += f_x
            particle_1.f_y += f_y
            particle_2.f_x -= f_x
            particle_2.f_y -= f_y
        end
    end
end

"""
    compute_interations!(hcosθ; period_x, period_y)

Compute `HarmonicCosingAngle` forces for triplets of particles stored in
`hcosθ.particle_triplets`.  Apply periodic boundary conditions if
`period_x > 0.0` or `period_y > 0.0`.
"""
function compute_interactions!(hcosθ::HarmonicCosineAngle; period_x::Float64 = -1.0, period_y::Float64 = -1.0)
    @use_threads hcosθ.multithreaded for (particle_left, particle_mid, particle_right) in hcosθ.particle_triplets
        x_l, y_l = particle_left.x, particle_left.y
        x_m, y_m = particle_mid.x, particle_mid.y
        x_r, y_r = particle_right.x, particle_right.y
        Δx_lm = wrap_displacement(x_l - x_m; period = period_x)
        Δy_lm = wrap_displacement(y_l - y_m; period = period_y)
        Δr_lm² = Δx_lm^2 + Δy_lm^2
        Δx_rm = wrap_displacement(x_r - x_m; period = period_x)
        Δy_rm = wrap_displacement(y_r - y_m; period = period_y)
        Δr_rm² = Δx_rm^2 + Δy_rm^2

        norm = 1.0 / sqrt(Δr_lm² * Δr_rm²)
        cosθ = (Δx_lm * Δx_rm + Δy_lm * Δy_rm) * norm
        sinθ = (Δx_lm * Δy_rm - Δx_rm * Δy_lm) * norm
        coef = -hcosθ.k_cosθ * (cosθ - hcosθ.cosθ_rest) * sinθ

        # particle_left force
        coef_l = coef / Δr_lm²
        particle_left.f_x += (f_x_l = -coef_l * Δy_lm)
        particle_left.f_y += (f_y_l = coef_l * Δx_lm)

        # particle_right force
        coef_r = coef / Δr_rm²
        particle_right.f_x += (f_x_r = coef_r * Δy_rm)
        particle_right.f_y += (f_y_r = -coef_r * Δx_rm)

        # particle_mid force
        particle_mid.f_x -= f_x_l + f_x_r
        particle_mid.f_y -= f_y_l + f_y_r
    end
end