text:
#functions and extensiosn for 2SLS
#WARNING: Probably broken, hasn't been maintained in a while

######2SLS constructor
#Holds the modelling componenets for a 2SLS model (except errors)
#Variable descriptions are in the struct
# Runs standard 2SLS using QR optimization
mutable struct FM2SLS <: FMModel
    zaqr::FMQR #QR decomposition of Za (Z|W)
    xaqr::FMQR #QR decomposition of fitted 1st stage values (X̃|W)

    X::Matrix{Float64} #Endogeneous covariates
    W::Matrix{Float64} #Exogeneous covariates
    Y::Vector{Float64} #Dependent var
    Z::Matrix{Float64} #IV data

    N::Int #number of data points
    K::Int # of endogeneous covariates
    L::Int # of instrumental variables
    KW::Int # of exogeneous covariates
    dof::Int

    Π1::Matrix{Float64} #first stage results
    Ξ1::Matrix{Float64} #first stage residuals

    δ2::Vector{Float64} #2SLS coefficients
    ξ2::Vector{Float64} #2SLS residuals

    XNames::Vector{Symbol} #names of the X variables
    WNames::Vector{Symbol} #names of the W variables
    YName::Symbol #names of the Y variables
    ZNames::Vector{Symbol} #names of the Z variables

    #FM2SLS()
end

#= Convenience constructor for FM2SLS that allocates a 0 element array for the
case of no exogeneous covariates. NOTE: no exog covariates implies no intercept,
so there probably should be  some exogeneous covariates
#IN: X Matrix, Y Vector, and the Z matrix of Instrumental variables
#OUT: FM2SLS object=#
FM2SLS(X::Matrix{Float64},Y::Vector{Float64},Z::Matrix{Float64};
    XNames::Vector{Symbol} = Symbol.(:X, 1:size(X,2)),
    YName::Symbol = :Y,
    ZNames::Vector{Symbol} = Symbol.(:Z, 1:size(Z,2)))::FM2SLS =
        FM2SLS(X,Matrix{Float64}(undef, Size(X,1),0),Y,Z)

#=Main constructor for FM2SLS object
IN: Endogeneous matrix X, Exogeneous matrix W, Dependent vector Y,
instrumental zariable matrix Z
OUT: A 2SLS object

#NOTE: not sure about the dof
=#
function FM2SLS(X::Matrix{Float64}, W::Matrix{Float64},
  Y::Vector{Float64},Z::Matrix{Float64};
  XNames::Vector{Symbol} = Symbol.(:X, 1:size(X,2)),
  WNames::Vector{Symbol} = Symbol.(:W, 1:size(W,2)),
  YName::Symbol = :Y,
  ZNames::Vector{Symbol} = Symbol.(:Z, 1:size(Z,2)),
  dof::Int = size(X,1) - size(X,2) - size(W,2))::FM2SLS

  #memory pre-allocaiton
  N::Int = size(X,1) #number of data points
  K::Int = size(X,2) # of endogeneous variables
  L::Int = size(Z,2) # of instrumental variables
  KW::Int = size(W,2)# of exogeneous covariates

  Π1::Matrix{Float64} = Matrix{Float64}(undef, L+KW,K)
  X̂::Matrix{Float64} = Matrix{Float64}(N,K) #holds fitted first stage values
  Ξ1::Matrix{Float64} = similar(X̂) #holds first stage residuals
  δ2::Vector{Float64} = Vector{Float64}(undef, K+KW)
  ξ2::Vector{Float64} = similar(Y)
  #Xa::Matrix{Float64} = Matrix{Float64}(N,K+KW) #fitted 1st stage plus exog covariates

  #start with setting up and running the first stage
  ZW::Matrix{Float64} = [Z W]
  zaqr::FMQR = FMQR(ZW, keepX=false) #augment the matrices to ZW and get the QR
  getCoeff!(zaqr, X, Π1) #get the first stage coefficients
  BLAS.gemm!('N','N',1.0,ZW,Π1,0.0,X̂) #get the fitted first stage values
  getResid!(ZW, X, Π1, Ξ1) #get the residuals

  #now run the second stage
  xaqr::FMQR = FMQR([X̂ W], keepX = false) #augment X̂W and get the QR
  getCoeff!(xaqr, Y, δ2) #get the 2SLS coefficients
  getResid!([X W], Y, δ2, ξ2) #get the 2SLS residuals

  #Finally, return the constructed model
  return FM2SLS(zaqr, xaqr, X, W, Y, Z, N, K, L, KW, dof, Π1, Ξ1, δ2, ξ2,
    XNames, WNames, YName, ZNames)
end

#Constructor and helper method which gets required info from DataFrame
#IN: The source dataframe, a formula expression/symbol for X, W, Y, and Z
# the dependent variable as a symbol
#OUT: A FMLM object
function FM2SLS(df::DataFrame, XExpr::T, WExpr::U, YSym::Symbol, ZExpr::V;
    XNames::Vector{Symbol} = Vector{Symbol}(),
    WNames::Vector{Symbol} = Vector{Symbol}([:Intercept]),
    YName::Symbol = :Y,
    ZNames::Vector{Symbol} = Vector{Symbol}(),
    containsmissings=true)::FM2SLS where {T <: FMExpr, U <: FMExpr, V <: FMExpr}

    #get the list of symbols for X, W and Z
    XSym::Vector{Symbol} = getSymbolsFromExpr(XExpr)
    WSym::Vector{Symbol} = getSymbolsFromExpr(WExpr)
    ZSym::Vector{Symbol} = getSymbolsFromExpr(ZExpr)

    #make a view so if we drop nulls it doesn't affect the original
    dfSub::SubDataFrame = view(df[ [XSym; WSym; YSym; ZSym]],1:size(df,1),:)

    # if we are going to drop nulls
    if containsmissings
        dfSub = cloakmissings(dfSub)
    end

    #get the model matrices, which will expand all factors
    XModelMatrix::Matrix{Float64} = getModelMatrix(dfSub, XExpr)
    WModelMatrix::Matrix{Float64} = getModelMatrix(dfSub, WExpr)
    ZModelMatrix::Matrix{Float64} = getModelMatrix(dfSub, ZExpr)

    #this is a check to make sure the factor expansion follows other columns
    #We must repeat the check on each variable matrix
    for symVec::Vector{Symbol} ∈ [XSym, WSym, ZSym]
        if length(symVec) > 1
            for i ∈ 2:length(symVec) #all types
                if typeof(dfSub[symVec[i-1]])<:CategoricalArray &&
                    !(typeof(dfSub[symVec[i]])<:CategoricalArray)
                    @warn "Type of $(typeof(dfSub[i])) is preceeded by type of
                        factor. Column names are likely MISALIGNED. Put factors
                        after numerical variables in RHS expressions."
                end
            end
        end
    end

    #now create the namesfor the expanded matrix
    XNamesFull::Vector{Symbol} = Vector{Symbol}(undef, size(XModelMatrix,2))
    for i ∈ 1:length(XNamesFull)
        XNamesFull[i] = i≤length(XNames) ? XNames[i] : Symbol(:X,i)
    end

    WNamesFull::Vector{Symbol} = Vector{Symbol}(undef, size(WModelMatrix,2))
    for i ∈ 1:length(WNamesFull)
        WNamesFull[i] = i≤length(WNames) ? WNames[i] : Symbol(:W,i)
    end

    ZNamesFull::Vector{Symbol} = Vector{Symbol}(undef, size(ZModelMatrix,2))
    for i ∈ 1:length(ZNamesFull)
        ZNamesFull[i] = i ≤ length(ZNames) ? ZNames[i] : Symbol(:Z,i)
    end

    #get the appropriate matrices and construct the FM2SLS object
    return FM2SLS(XModelMatrix, WModelMatrix, Vector{Float64}(undef, dfSub[YSym]),
        ZModelMatrix, XNames=XNamesFull, WNames=WNamesFull,
        YName=YName, ZNames=ZNamesFull)

end

#default constructor to help with memory management
FM2SLS() = FM2SLS(FMQR(), FMQR(), Matrix{Float64}(), Matrix{Float64}(),
    Vector{Float64}(), Matrix{Float64}(), 0, 0, 0, 0, Matrix{Float64}(),
    Matrix{Float64}(), Vector{Float64}(), Vector{Float64}(),
    Vector{Symbol}(), Vector{Symbol}(), Symbol(), Vector{Symbol}())::FM2SLS

###########################2SLS Helper functions

#=function get1stStage
#The purpose is to get the first stage of the IV as a linear model
#Unfortunately, we need to copy some of the data, so this can be expensive,
#although Z and W should only be copied once
#IN: an IV object
#OUT: a vector of FMLM objects, one for each focal (X) variable in the
#initial regression=#
function get1stStage(iv::FM2SLS)::Vector{FMLM}
    ZW::Matrix{Float64} =  [iv.Z iv.W]
    iv1st::Vector = Vector{FMLM}()

    for i ∈ 1:length(iv.XNames)
        push!(iv1st, FMLM(iv.zaqr, [iv.Z iv.W], iv.X[:,i],
            iv.zaqr.N, iv.zaqr.K, iv.Π1[:,i], iv.Ξ1[:,i],
            [iv.ZNames; iv.WNames], iv.XNames[i]))
    end

    return iv1st
end

#=gets the residuals in the 2SLS framework
Note this is a wrapper function which can accept multiple independent
and multiple dependent variables. NOTE the change in the naming convention
relative to other verisons of this method (due to the 2SLS context)
IN: QR decomposition of the RHS variables, LHS dependent variable matrix,
a coefficient matrix, optionally the memory for the residual matrix
OUT: Writes and returns the residual matrix=#
function getResid!(ZW::Matrix{Float64}, X::Matrix{Float64}, Π1::Matrix{Float64},
  Ξ1::Matrix{Float64}=similar(X))::Matrix{Float64}

  #get the residuals by reference
  for i ∈ 1:size(X,2) #get the vector of reisduals multiple times
    getResid!(ZW, view(X,:,i), view(Π1,:,i), view(Ξ1,:,i))
  end
  return Ξ1
end

#=This is a wrapper function which gets a matrix of coefficients for use in
# the first stage of 2SLS. NOTE the change in the naming convention
relative to other verisons of this method (due to the 2SLS context)
IN: QR decomposition of the RHS variables, LHS dependent variable matrix,
optionally the memory for the coefficient matrix
OUT: Writes and returns the coefficient matrix=#
function getCoeff!(zaqr::FMQR, X::Matrix{Float64},
  Π1::Matrix{Float64} = Matrix{Float64}(undef, zaqr.K, size(X,2)))::Matrix{Float64}

  for i ∈ 1:size(X,2) # for each dependent variable vector
    getCoeff!(zaqr, view(X,:,i), view(Π1,:,i)) #run the regression as needed
  end

  return Π1
end

################################2SLS Standard Errors

#Convenience method for above: IN: A FM2SLS object #OUT: A covariance matrix
getHomoskedΣSlow(iv::FM2SLS)::Matrix{Float64} =
  getHomoskedΣSlow([[iv.Z iv.W]*iv.Π1 iv.W], iv.Y.-[iv.X iv.W] * iv.δ2)

#=helper method for the above function which extracts the QR decomposition
from a FM2SLS object
IN: A linear model, optionally memory for the covariance matrix
OUT: WRites and returns the covariance matrix=#
function getHomoskedΣ!(iv::FM2SLS,
  Σ::Matrix{Float64}=Matrix{Float64}(undef, iv.xaqr.K, iv.xaqr.K))::Matrix{Float64}

  return getHomoskedΣ!(iv.xaqr, iv.ξ2, Σ, iv.N / iv.dof)
end

#=helper method which extracts the required components from a 2SLS model
IN: A 2SLS model, optionally memory for the covariance matrix
OUT: Writes and allocates the covariance matrix=#
function getModWhiteΣ!(iv::FM2SLS,
  Σ::Matrix{Float64}=Matrix{Float64}(undef, iv.xaqr.K, iv.xaqr.K))::Matrix{Float64}

  return getModWhiteΣ!(iv.xaqr, iv.ξ2, Σ, iv.N/iv.dof)
end

#Convenience method for above: IN: A 2SLS object #OUT: A covariance matrix
getModWhiteΣSlow(iv::FM2SLS)::Matrix{Float64} =
  getModWhiteΣSlow([[iv.Z iv.W]*iv.Π1 iv.W], iv.Y.-[iv.X iv.W] * iv.δ2)

###########################IO
#if we get serious about 2SLS again these can go into the IO file

function texTable(models::Vector{FM2SLS}, getΣ::T where T<:Union{Function, Vector{Function}},
    rows::Vector{Symbol};
    titleCaption::String = "",
    caption::String = "",
    colNames::Vector{Vector{String}} = [["C$i" for i ∈ 1:length(models)]],
    contentRowNames::Vector{String} = String.(rows),
    descRowNames::Vector{String}=Vector{String}(),
    descContent::Vector{Vector{String}}=Vector{Vector{String}}(),
    notes::Vector{String} = Vector{String}(),
    arrayStretch::Float64 = 1.5,
    lineSpacer::String = "\\\\",
    summaryMathMode::Bool = true,
    widthColNames::Vector{Vector{Int}} =
        broadcast((i::Int)->ones(Int,length(colNames[i])),1:length(colNames)),
    alignmentColNames::Vector{Vector{String}} = #contains the number of columns for each entry
      broadcast((i::Int)->["r" for i ∈ 1:length(colNames[i])],1:length(colNames)),
    widthDescContent::Vector{Vector{Int}} =
        broadcast((i::Int)->ones(Int,length(descContent[i])),1:length(descContent)),
    stars::Bool=true,
    starLvls::Vector{Float64} = [.9, .95, .99],
    starLegend::String = stars ? DEFAULT_STAR_LEGEND : "",
    starStrings::Vector{String} =
      ["\\ensuremath{^\\text{*}}","\\ensuremath{^\\text{**}}","\\ensuremath{^\\text{***}}"],
    scaling::Vector{Float64}=ones(length(rows)),
    decimalDigits::Int = 2,
    colHeaderName::Vector{String} = ["" for i ∈ 1:length(colNames)],
    nakedTable=nakedTable,
    alignmentstring::String = string(" l | ", join(["r" for i ∈ 1:length(descContent[1])])))

    numCols = length(models)
    numContentRows = length(rows)

    #Pre-allocate the vector of SE errors
    modelsσ::Vector{Vector{Float64}} =
    [Vector{Float64}(models[i].K+models[i].KW) for i∈1:numCols]

    #pull out the β coefficients and N
    modelsβ::Vector{Vector{Float64}} = [models[i].δ2 for i ∈ 1:numCols]
    modelsN::Vector{Int} = [models[i].N for i ∈ 1:numCols]
    modelsXWNames::Vector{Vector{Symbol}} = [[models[i].XNames; models[i].WNames]  for i ∈ 1:numCols]

    #get the standard errors
    for c ∈ 1:numCols
        if typeof(getΣ) <: Function
          modelsσ[c] .= sqrt.(diag(getΣ(models[c])))
        else
          modelsσ[c] .= sqrt.(diag(getΣ[c](models[c])))
        end
    end

    if length(starLegend) > 0
        notes = [starLegend; notes]
    end

    #get the content matrices
    content::Vector{Matrix{String}} = getContentMatrices!( modelsβ, modelsσ, modelsXWNames, modelsN, rows,
      stars=stars, starLvls=starLvls, scaling=scaling,
      decimalDigits=decimalDigits, starStrings=starStrings)

    texTable(titleCaption = titleCaption,caption=caption, colNames=colNames,
      contentRowNames=contentRowNames, content=content, descRowNames=descRowNames,
      descContent=descContent, notes=notes, arrayStretch=arrayStretch, lineSpacer=lineSpacer,
      summaryMathMode=summaryMathMode, widthColNames=widthColNames,
      widthDescContent=widthDescContent, colHeaderName=colHeaderName,
      alignmentColNames=alignmentColNames, nakedTable=nakedTable, alignmentstring=alignmentstring)
end

function IVTest()

  srand(11)

  N::Int = 200 #Number of Observations
  K::Int = 3 # of endogneous covariates
  L::Int = 4 # of instruments
  KW::Int = 5# of exog covariates

  #Allocate
  X::Matrix{Float64} = Matrix{Float64}(undef, N,K)
  Z::Matrix{Float64} = Matrix{Float64}(undef, N,L)
  W::Matrix{Float64} = Matrix{Float64}(undef, N,KW)
  Y::Vector{Float64} = Vector{Float64}(undef, N)
  A::Vector{Float64} = Vector{Float64}(undef, N)

  Π1::Matrix{Float64} = Matrix{Float64}(undef, L+KW,K)
  π2::Vector{Float64} = Vector{Float64}(undef, L+KW)
  Ξ1::Matrix{Float64} = similar(X)
  e2::Vector{Float64} = similar(Y)

  #parameters for the simulation
  σ2::Vector{Float64} = [2.0^2 for i::Int ∈ 1:N]
  σ2[1:ceil(Int, N/2)] .= 0.5^2.0

  #generate π2 coefficients (pattern is 1 times the index, ie 1 2 3 ...)
  π2 .= [(10*i) for i::Int ∈ 1:(L+KW)]
  π2[L+1] .= 1.0 #intercept term
  #println("π2: ",π2)

  #generate Π1, pattern is δ2 times 1/10 the index, ie ie .1δ2 .2δ2 .3δ2 ...)
  for i ∈ 1:(L+KW), j∈1:K
    Π1[i,j] .= i + j^2
  end
  Π1[L+1,:] .= 1.0 #set the intercept term
  #display(Π1)

  #generate the exogeneous variables
  for i ∈ 2:KW
    W[:,i] .= rand(Uniform(-10,10),N)
  end
  W[:,1] .= 1.0 #intercept

  #generate the instrumental variables
  for i ∈ 1:L
    Z[:,i] .= rand(Uniform(-10,10),N)
  end

  #map((s2::Float64)->rand(Normal(0.0,s2^0.5)),σ2)
  #generate the errors
  for i ∈ 1:N
    for j ∈ 1:K
      Ξ1[i,j] = rand(Normal(0.0,σ2[i]^0.5))
    end
    e2[i] = rand(Normal(0.0,σ2[i]^0.5))
  end

  #println(Ξ1)
  #populate Y and X
  X .= [Z W]*Π1 .+ Ξ1
  Y .= [Z W]*π2 .+ e2

  #add some endogeneity
  for i ∈ 1:K
    A .= rand(Uniform(-i,i),N)
    X[:,i] .+= A
    Y .+= A .* 10
  end

  #probably is a more graceful way to do this
  XNames::Vector{Symbol} = [Symbol("X$i") for i = 1:K]
  WNames::Vector{Symbol} = [Symbol("W$i") for i = 1:KW]
  ZNames::Vector{Symbol} = [Symbol("Z$i") for i = 1:L]

  iv::FM2SLS = FM2SLS(X,W,Y,Z, XNames=XNames, YName=:Y, WNames=WNames, ZNames=ZNames)
  println("coefficients: \n", [XNames; WNames] ,"\n", iv.δ2)

  #get the homoskedastic SEs
  ΣHomosked::Matrix{Float64} = getHomoskedΣ!(iv)
  ΣHomoskedSlow::Matrix{Float64} = getHomoskedΣSlow(iv)

  #print the coefficients
  println("Homoskedastic Errors: ", diag(ΣHomosked).^.5)
  println("Check: ", diag(ΣHomoskedSlow).^.5)

  #get the modified white SEs
  ΣWhite::Matrix{Float64} = getModWhiteΣ!(iv)
  ΣWhiteSlow::Matrix{Float64} = getModWhiteΣSlow(iv)

  #print the coefficients
  println("Modified White Errors: ",diag(ΣWhite).^.5)
  println("Check: ", diag(ΣWhiteSlow).^.5)

  iv1st::Vector{FMLM} = get1stStage(iv)
  println("1st Stage results: \n $(iv1st[1].XNames)")
  for iv1 ∈ iv1st
      println("coefficients $(iv1.YName):  $(iv1.β )")
      println("Error (Homosked $(iv1.YName)): $(diag(getHomoskedΣ!(iv1)).^.5)")
  end

  oStream::IOStream = open(
    "C:\\Users\\Clinton\\Dropbox\\Projects\\Summer17RTester\\Summer17RTester\\IVTestOut.csv","w+")

  write(oStream, [string.(iv.XNames); string.(iv.WNames)]::Vector{String})
  writecsv(oStream, [X W Y Z])
  close(oStream)
end

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
#!/usr/bin/env julia

const DEBUG = false
if !DEBUG
    const IN = 32
    const OUT = 64
    const NUM_ITS = 10 ^ 6
else
    const IN = 4
    const OUT = 2
    const NUM_ITS = 1
end
const INIT_VAL = 0.4711

function forward(W, x, a)
    # a = tanh(W * x)
    for row in 1:size(W, 1)
        a[row] = 0
        for col in 1:size(W, 2)
            @inbounds a[row] += W[row, col] * x[col]
        end
        a[row] = tanh(a[row])
    end
end

function tanh_prime(x)
    return 1 - x .^ 2
end

function backward(x, a, d, b)
    # b = return (tanh_prime(a) .* d) * x'
    for row in 1:size(b, 1)
        @inbounds ad = (1.0 - a[row] ^ 2) * d[row]
        for col in 1:size(b, 2)
            @inbounds b[row, col] = ad * x[col]
        end
    end
end

W = fill!(Array(Float64, OUT, IN), INIT_VAL)
x = fill!(Array(Float64, IN), INIT_VAL)
d = fill!(Array(Float64, OUT), INIT_VAL)

a = Array(Float64, OUT)
b = Array(Float64, OUT, IN)

tic()
for _ in 1:NUM_ITS
    forward(W, x, a)
    backward(x, a, d, b)

    # Avoids scoping issues and has a minimal, if any, effect on performance.
    if DEBUG
        print(a')
        print(b)
    end
end
toc = toq()

println(toc)

ext: jl

max_stars_repo_name: ninjin/ppod

max_issues_repo_name: ninjin/ppod

max_forks_repo_name: ninjin/ppod


text:
# Autogenerated wrapper script for GR_jll for aarch64-linux-gnu-cxx03
export gksqt, libGKS, libGR, libGR3, libGRM

using Bzip2_jll
using Cairo_jll
using FFMPEG_jll
using Fontconfig_jll
using GLFW_jll
using JpegTurbo_jll
using libpng_jll
using Libtiff_jll
using Pixman_jll
using Qt5Base_jll
using Zlib_jll
JLLWrappers.@generate_wrapper_header("GR")
JLLWrappers.@declare_library_product(libGKS, "libGKS.so")
JLLWrappers.@declare_library_product(libGR, "libGR.so")
JLLWrappers.@declare_library_product(libGR3, "libGR3.so")
JLLWrappers.@declare_library_product(libGRM, "libGRM.so")
JLLWrappers.@declare_executable_product(gksqt)
function __init__()
    JLLWrappers.@generate_init_header(Bzip2_jll, Cairo_jll, FFMPEG_jll, Fontconfig_jll, GLFW_jll, JpegTurbo_jll, libpng_jll, Libtiff_jll, Pixman_jll, Qt5Base_jll, Zlib_jll)
    JLLWrappers.@init_library_product(
        libGKS,
        "lib/libGKS.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libGR,
        "lib/libGR.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libGR3,
        "lib/libGR3.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libGRM,
        "lib/libGRM.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_executable_product(
        gksqt,
        "bin/gksqt",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using Revise # reduce recompilation time
using Plots
# using Documenter
using LinearAlgebra
using SparseArrays
using BenchmarkTools
using UnPack

push!(LOAD_PATH, "./examples/EntropyStableEuler.jl/src")
using EntropyStableEuler
using EntropyStableEuler.Fluxes1D
const γ = 1.4

function wavespeed_1D(rho,rhou,E)
    p = pfun_nd(rho,(rhou,),E)
    cvel = @. sqrt(γ*p/rho)
    return @. abs(rhou/rho) + cvel
end
unorm(U) = sum(map((x->x.^2),U))
function pfun_nd(rho, rhoU, E)
    rhoUnorm2 = unorm(rhoU)./rho
    return @. (γ-1)*(E - .5*rhoUnorm2)
end

K = 200
T = 6.0

# Bl = -0.5
# Br = 0.5
Bl = 0.0
Br = 9.0
VX = LinRange(Bl,Br,K+1)
# Assume uniform intervals
m_i = (Br-Bl)/K # mass
C = spdiagm(-1 => -1/2*ones(K-1), 1 => 1/2*ones(K-1))# stiffness
C[1,end] = -1/2
C[end,1] = 1/2

x = zeros(K,1)
for i = 1:K
    x[i] = (VX[i+1]+VX[i])/2
end
"""Initial condition"""
rho_x(x) = (x <= 0.0) ? 1.0 : 0.125
u_x(x) = 0.0
p_x(x) = (x <= 0.0) ? 1.0 : 0.1
# rho_x(x) = (x <= 0.0) ? 1.0 : 0.5
# u_x(x) = 1.0
# p_x(x) = (x <= 0.0) ? 1.0 : 1.0
rho_x(x) = (x <= 3.0) ? 1.0 : 0.001
u_x(x) = 0.0
p_x(x) = (x <= 3.0) ? 0.1 : 1e-7

rho = @. rho_x(x)
u = @. u_x(x)
p = @. p_x(x)
Q = primitive_to_conservative(rho,u,p)

function rhs_IDP(Q,K,m_i,C)
    p = pfun_nd.(Q[1],Q[2],Q[3])
    flux = zero.(Q)
    @. flux[1] = Q[2]
    @. flux[2] = Q[2]^2/Q[1]+p
    @. flux[3] = Q[3]*Q[2]/Q[1]+p*Q[2]/Q[1]

    dfdx = (zeros(K),zeros(K),zeros(K))
    for i = 1:K
        for c = 1:3
            dfdx[c][i] += C[i,mod1(i-1,K)]*(flux[c][mod1(i-1,K)]+flux[c][i])+
                          C[i,mod1(i+1,K)]*(flux[c][mod1(i+1,K)]+flux[c][i])
        end
    end

    visc = (zeros(K),zeros(K),zeros(K))
    for i = 1:K
        wavespd_curr = abs(wavespeed_1D(Q[1][i],Q[2][i],Q[3][i]))
        wavespd_R = abs(wavespeed_1D(Q[1][mod1(i+1,K)],Q[2][mod1(i+1,K)],Q[3][mod1(i+1,K)]))
        wavespd_L = abs(wavespeed_1D(Q[1][mod1(i-1,K)],Q[2][mod1(i-1,K)],Q[3][mod1(i-1,K)]))
        wavespd_curr_m = abs(wavespeed_1D(Q[1][i],-Q[2][i],Q[3][i]))
        wavespd_R_m = abs(wavespeed_1D(Q[1][mod1(i+1,K)],-Q[2][mod1(i+1,K)],Q[3][mod1(i+1,K)]))
        wavespd_L_m = abs(wavespeed_1D(Q[1][mod1(i-1,K)],-Q[2][mod1(i-1,K)],Q[3][mod1(i-1,K)]))
        dL = 1/2*max(wavespd_curr,wavespd_L)
        dR = 1/2*max(wavespd_curr,wavespd_R)
        for c = 1:3
            visc[c][i] = dL*(Q[c][mod1(i-1,K)]-Q[c][i]) + dR*(Q[c][mod1(i+1,K)]-Q[c][i])
        end
    end

    return -1/m_i.*(dfdx.-visc)
end

gr(size=(300,300),ylims=(0,1.2),legend=false,markerstrokewidth=0,markersize=1)
plt = plot()

dt = 0.0001
t = 0.0
Nsteps = Int(T/dt)
Q = collect(Q)
@gif for i = 1:Nsteps
    rhsQ = rhs_IDP(Q,K,m_i,C)
    global @. Q =  Q + dt*rhsQ

    if i%10==0 || i==Nsteps
        println("Number of time steps $i out of $Nsteps")

        scatter(x,Q[1])
        #for k = 1:K
            #xp = LinRange(VX[k], VX[k+1], 3)
            #scatter!(xp,Q[1][k]*ones(length(xp)))
        #end
    end
end every 1000

# p = pfun_nd.(Q[1],Q[2],Q[3])
# for k = 1:K
#     xp = LinRange(VX[k], VX[k+1], 3)
#     scatter!(xp,Q[1][k]*ones(length(xp)))

# end
# plt

max_stars_repo_name: yiminllin/ESDG-PosLimit

max_issues_repo_name: yiminllin/ESDG-PosLimit

max_forks_repo_name: yiminllin/ESDG-PosLimit


text:
module Symmetry

import LinearAlgebra

import MutableArithmetics
const MA = MutableArithmetics
import MultivariatePolynomials
const MP = MultivariatePolynomials
import MultivariateBases
const MB = MultivariateBases

import SymbolicWedderburn

import SumOfSquares

struct Pattern{GT,AT<:SymbolicWedderburn.Action}
    group::GT
    action::AT
end

include("utils.jl")
include("wedderburn.jl")
include("block_diag.jl")

end

max_stars_repo_name: blegat/SumOfSquares

max_issues_repo_name: blegat/SumOfSquares

max_forks_repo_name: blegat/SumOfSquares


text:
module MultivariateMoments

using LinearAlgebra

using MultivariatePolynomials
const MP = MultivariatePolynomials

export AbstractMeasureLike, AbstractMomentLike
abstract type AbstractMeasureLike{T} end
abstract type AbstractMomentLike{T} <: AbstractMeasureLike{T} end

export AbstractMeasure, AbstractMoment
abstract type AbstractMoment{T} <: AbstractMomentLike{T} end
abstract type AbstractMeasure{T} <: AbstractMeasureLike{T} end

include("moment.jl")
include("measure.jl")
include("expectation.jl")
include("symmatrix.jl")
include("moment_matrix.jl")
include("atomic.jl")
include("extract.jl")

end # module


text:
using DataStructures

function generateHypergraph()
    h = Hypergraph{Bool}(5,3)
    h[[1,2,5],1] .= true;
    h[:, 2] .= true;
    h[[1,3,5],3] .= true;
    metaV = [2,2,1,1,2]
    metaE = [2,2,2]
    actV = zeros(Bool, nhv(h));
    actE = zeros(Bool, nhe(h));
    actV[[1,4]] .= true;
    (h, actV, actE, metaV, metaE)
end

function simulate!(h::Hypergraph{Bool},
                  actV::Vector{Bool}, actE::Vector{Bool},
                  metaV::Vector{Int}, metaE::Vector{Int};
                  printme = true, max_step=1_000_000 )
    step = 0
    while true && step < max_step
        step += 1
        actE_cp = deepcopy(actE)
        #deb = DataFrame(step=Int[], aSum)
        for e in 1:nhe(h)
            aSum = sum(actV[ collect(keys(getvertices(h,e))) ])
            if aSum >= metaE[e]
                actE_cp[e] = true
            end
            printme && println("$step e=$e $aSum $(metaE[e]) $(actE_cp[e])")
        end
        #sum(actE_cp) == sum(actE) && return (actvs = sum(actV), step=step-1)
        actE .= actE_cp
        actV_cp = deepcopy(actV)
        for v in 1:nhv(h)
            aSum = sum(actE[ collect(keys(gethyperedges(h,v))) ])
            if aSum >= metaV[v]
                actV_cp[v] = true
            end
            printme && println("$step v=$v $aSum $(metaV[v]) $(actV_cp[v])")
        end
        sum(actV_cp) == sum(actV) && return (actvs = sum(actV), step=step-1)
        actV .= actV_cp
    end
    step
end

function randMetaV(h)
    metaV = Vector{Int}(undef, nhv(h))
    for v in 1:nhv(h)
        metaV[v] = rand(1:length(h.v2he[v]))
    end
    metaV
end

function randMetaE(h)
    metaE = Vector{Int}(undef, nhe(h))
    for e in 1:nhe(h)
        metaE[e] = rand(1:length(h.he2v[e]))
    end
    metaE
end

function proportionalMetaV(h::Hypergraph,prop)
    @assert 0.0 < prop <= 1.0
    metaV = Vector{Int}(undef, nhv(h))
    for v in 1:nhv(h)
        metaV[v] = ceil(length(h.v2he[v])*prop)
    end
    metaV
end

function proportionalMetaE(h::Hypergraph,prop)
    @assert 0.0 < prop <= 1.0
    metaE = Vector{Int}(undef, nhe(h))
    for e in 1:nhe(h)
        metaE[e] = ceil(length(h.he2v[e])*prop)
    end
    metaE
end


text:
"""
TopDistance(args...)

Estimate the distance between points given. Plot either distance between each point or distance with respect to specific point

``D=\\sqrt { { \\left( { X }_{ 1 }-{ X }_{ 2 } \\right)  }^{ 2 }+{ \\left( { Y }_{ 1 }-{ Y }_{ 2 } \\right)  }^{ 2 }{ +\\left( { Z }_{ 1 }-{ Z }_{ 2 } \\right)  }^{ 2 } } ``

The next table show the list of variables allowed:

|PropertyName|Args|Default|Input|Description
|---|---|---|---|---|
|x|Mandatory|--|Array{number,1}|X coordinate of points|
|y|Mandatory|--|Array{number,1}|Y coordinate of points|
|z|Optional|z=false|Array{number,1}|z coordinate of points. If not set, a zeros array is established|
|Names|Optional|Names=false|Array{String,1}|Name of each point. If not set a numbered list is created as a name|
|Show|Optional|Show=false|Show=Number|Refers to a numbered point which distance between other points and this are shown <br> if not set, all distances are shown|
"""
@userplot topdistance

@recipe function f(h::topdistance;
            Names=false,Show=false)

    if length(h.args)==2
        x, y = h.args
    elseif length(h.args)==3
        x, y, z = h.args
    end

    n=size(x,1)  #Number of inputs
    list=1:n    #Numbered list of elements

    if length(h.args)==2      #If height not provided, zero vector is created
        z=zeros(n)
    end

    if Names==false  #If names are not provided, it is created a numbered list
        Names=1:n
    end

        WellSym=map(x->Symbol(x),Names)     #DataFrame COlumns Name
        Distance=DataFrame()

        for i=1:n
            Distance[WellSym[i]]=zeros(n)
            for j=1:n
              Distance[WellSym[i]][j]=sqrt( (x[i]-x[j])^2  + (y[i]-y[j])^2 + (z[i]-z[j])^2)
            end
        end

    if Show==false                         # If Show==false plot each distance of each well, otherwise plot distance Well choosed
        m=binomial(n,2)
        AvrD=DataFrame(X=zeros(m),Y=zeros(m),Dist=zeros(m))
        k=0
           for i=1:n-1
               for j=i+1:n
                k=k+1
                AvrD.X[k]=mean([x[i],x[j]])
                AvrD.Y[k]=mean([y[i],y[j]])
                AvrD.Dist[k]=round(Distance[i,j],digits=2)
                end
            end

    else
        m=n-1
        AvrD=DataFrame(X=zeros(m),Y=zeros(m),Dist=zeros(m))
        c=list[list.!=Show]
        k=0

        for i in c
            k=k+1
            AvrD.X[k]=mean([x[Show],x[i]])
            AvrD.Y[k]=mean([y[Show],y[i]])
            AvrD.Dist[k]=round(Distance[Show,i],digits=2)
        end
    end

# Plot the Distances

@series begin
    seriestype := :scatter
    seriescolor := :black
    series_annotations := map(x->Plots.text(x,:left,:top,8),Names)
    x, y
end
#
@series begin
        seriestype := :scatter
        seriescolor := :white
        series_annotations := map(x->Plots.text("$x m" ,:right,:top,8),AvrD.Dist)
        AvrD.X,AvrD.Y
end

    ## Plot properties for union lines
    seriestype := :path
    linestyle := :dot
    arrow := Plots.arrow(:open, :both )
    seriescolor := :black
    legend --> false

if Show==false  # If Show==false plot each distance of each well, otherwise plot distance Well choosed
    for i=1:n-1
        for j=i+1:n

            @series begin
                        x[[i,j]],y[[i,j]]
                    end
        end
    end
else
    for i in c
                     @series begin
                        x[[Show,i]],y[[Show,i]]
                    end
    end

end

end


text:
include("StructDomain.jl")
include("GenerateGS.jl")
include("FEA.jl")
include("TrussTop.jl")

using BenchmarkTools

function TrussScript4()

    # ----- INPUT DATA -----

    Nx, Ny = 12, 12  	 	 		            # Number of cells
    Lx, Ly = 12., 12.  	 		                # Domain size
    Supp = [[0.,Ly/2., 1,1],[Lx,Ly/2., 1,1]]    # Support list: position and type
    Load = [[Lx/2.,Ly/2., 0.,-1.]]              # Load list: position and value
    GSLvl = 12 		 	 		                # Ground structure level
    GSColTol = 0.999999  		                # Ground structure colinear tolerance
    VolMax = Lx*Ly/9000. 		                # Volume constraint
    Tol = 1e-8			 		                # Optimization convergence tolerance
    MaxIter = 4000		 		                # Max number of optimization iterations
    E = 7e+7 			 		                # Young's modulus
    Filter = 0.01;                              # Max area fraction threshold (for vieweing)

    # ----- PRE PROCESSING -----

    # Create structured domain
    println("Creating structured domain ..")
    @time Node, Elem, Supp, Load = StructDomain(Nx, Ny, Lx, Ly, Supp, Load)

    # Generate Ground-Structure
    println("Generating Ground structure ..")
    @time Bars = GenerateGS(Node, Elem, GSLvl, GSColTol)

    # Init finite element analysis
    println("Initializing finite element analysis ..")
    @time FEA = InitFEA(Node, Bars, Supp, Load, E)

    # Calc initial areas
    Area = VolMax/sum(FEA.ElemLen)
    XInitial = Area*ones(FEA.NElem)

    # Init optimization
    println("Initializing optimization procedure ..")
    @time Opt = InitOPT(XInitial, VolMax, Tol, MaxIter)

    # ----- OPTIMIZATION -----

    # Run TrussTOP
    println("Running optimization procedure, TrussTOP ..")
    @time XHist, ObjFcnHist, FEA = TrussTop(FEA, Opt)

    # ----- POST PROCESSING -----

    # Show general results
    TrussTopView(XHist, ObjFcnHist, FEA, Filter)
end

TrussScript4()

max_stars_repo_name: ricardoaf/jutop

max_issues_repo_name: ricardoaf/jutop

max_forks_repo_name: ricardoaf/jutop


text:
@testset "source" begin

nX, nY, nZ = 1, 2, 3

@testset "distweights" begin
    ∆lg = 2:2:22  # make sure ∆l's are varying to prevent false positive
    l′ = cumsum([-1; ∆lg])  # [-1,1,5,11,19,29,41,55,71,89,109,131]: spacing is even

    lprim_g = StaggeredGridCalculus.movingavg(l′)  # [0,3,8,15,24,35,48,63,80,99,120]: spacing is odd, 11 entries (l has 10 entry)
    ldual_g = l′[1:end-1]  # [-1,1,5,11,19,29,41,55,71,89,109]: spacing is even

    domain = [lprim_g[1], lprim_g[end]]

    lprim, ∆lprim = lprim_g[1:end-1], diff(ldual_g)
    ldual, ∆ldual = ldual_g[2:end], diff(lprim_g)

    # Test for sources to assign at the primal grid points [0,3,8,15,24,35,48,63,80,99 (,120)].
    c = 8  # at grid point
    @test all(MaxwellFDFD.distweights(c, PRIM, domain, lprim, ∆lprim, true) .≈ ([3,3], [1/∆lprim[3], 0]))
    @test all(MaxwellFDFD.distweights(c, PRIM, domain, lprim, ∆lprim, false) .≈ ([3,3], [1/∆lprim[3], 0]))

    c = 11  # between grid points
    r = (c-8) / (15-8)
    @test all(MaxwellFDFD.distweights(c, PRIM, domain, lprim, ∆lprim, true) .≈ ([3,4], [1/∆lprim[3]*(1-r), 1/∆lprim[4]*r]))
    @test all(MaxwellFDFD.distweights(c, PRIM, domain, lprim, ∆lprim, false) .≈ ([3,4], [1/∆lprim[3]*(1-r), 1/∆lprim[4]*r]))

    c = 3  # farthest point affected by negative-end boundary
    @test all(MaxwellFDFD.distweights(c, PRIM, domain, lprim, ∆lprim, true) .≈ ([2,2], [1/∆lprim[2], 0]))
    @test all(MaxwellFDFD.distweights(c, PRIM, domain, lprim, ∆lprim, false) .≈ ([2,2], [1/∆lprim[2], 0]))

    c = 1  # within range affected by negative-end boundary
    r = (3-c) / (3-0)
    @test all(MaxwellFDFD.distweights(c, PRIM, domain, lprim, ∆lprim, true) .≈ ([1,2], [1/∆lprim[1]*r, 1/∆lprim[2]*(1-r)]))
    @test all(MaxwellFDFD.distweights(c, PRIM, domain, lprim, ∆lprim, false) .≈ ([2,2], [1/∆lprim[2]*(1-r), 0]))

    c = 0  # negative-end boundary
    @test all(MaxwellFDFD.distweights(c, PRIM, domain, lprim, ∆lprim, true) .≈ ([1,1], [1/∆lprim[1], 0]))
    @test all(MaxwellFDFD.distweights(c, PRIM, domain, lprim, ∆lprim, false) .≈ ([2,2], [0, 0]))

    c = -0.5  # within range affected by negative-end PDC, but outside Bloch and PPC domain (starting at l = 0)
    @test_throws ArgumentError MaxwellFDFD.distweights(c, PRIM, domain, lprim, ∆lprim, true)
    @test_throws ArgumentError MaxwellFDFD.distweights(c, PRIM, domain, lprim, ∆lprim, false)

    c = -1  # behind negative-end boundary
    @test_throws ArgumentError MaxwellFDFD.distweights(c, PRIM, domain, lprim, ∆lprim, true)
    @test_throws ArgumentError MaxwellFDFD.distweights(c, PRIM, domain, lprim, ∆lprim, false)

    c = 99  # farthest point affected by positive-end boundary
    @test all(MaxwellFDFD.distweights(c, PRIM, domain, lprim, ∆lprim, true) .≈ ([10,10], [1/∆lprim[10], 0]))
    @test all(MaxwellFDFD.distweights(c, PRIM, domain, lprim, ∆lprim, false) .≈ ([10,10], [1/∆lprim[10], 0]))

    c = 105  # within range affected by positive-end boundary
    r = (c-99) / (120-99)
    @test all(MaxwellFDFD.distweights(c, PRIM, domain, lprim, ∆lprim, true) .≈ ([10,1], [1/∆lprim[10]*(1-r), 1/∆lprim[1]*r]))
    @test all(MaxwellFDFD.distweights(c, PRIM, domain, lprim, ∆lprim, false) .≈ ([10,10], [1/∆lprim[10]*(1-r), 0]))

    c = 120  # positive-end boundary
    r = (c-99) / (120-99)
    @test all(MaxwellFDFD.distweights(c, PRIM, domain, lprim, ∆lprim, true) .≈ ([10,1], [1/∆lprim[10]*(1-r), 1/∆lprim[1]*r]))
    @test all(MaxwellFDFD.distweights(c, PRIM, domain, lprim, ∆lprim, false) .≈ ([10,10], [1/∆lprim[10]*(1-r), 0]))

    # Test for sources to assign at the dual grid points [(-1,) 1,5,11,19,29,41,55,71,89,109].
    # Note the domain boundaries are [0, 120].
    c = 11  # at grid point
    @test all(MaxwellFDFD.distweights(c, DUAL, domain, ldual, ∆ldual, true) .≈ ([3,3], [1/∆ldual[3], 0]))
    @test all(MaxwellFDFD.distweights(c, DUAL, domain, ldual, ∆ldual, false) .≈ ([3,3], [1/∆ldual[3], 0]))

    c = 15  # between grid points
    r = (c-11) / (19-11)
    @test all(MaxwellFDFD.distweights(c, DUAL, domain, ldual, ∆ldual, true) .≈ ([3,4], [1/∆ldual[3]*(1-r), 1/∆ldual[4]*r]))
    @test all(MaxwellFDFD.distweights(c, DUAL, domain, ldual, ∆ldual, false) .≈ ([3,4], [1/∆ldual[3]*(1-r), 1/∆ldual[4]*r]))

    c = 1  # farthest point affected by negative-end boundary
    @test all(MaxwellFDFD.distweights(c, DUAL, domain, ldual, ∆ldual, true) .≈ ([1,1], [1/∆ldual[1], 0]))
    @test all(MaxwellFDFD.distweights(c, DUAL, domain, ldual, ∆ldual, false) .≈ ([1,1], [1/∆ldual[1], 0]))

    c = 0.5  # within range affected by negative-end boundary
    r = (1-c) / ((1-0) + (120-109))
    @test all(MaxwellFDFD.distweights(c, DUAL, domain, ldual, ∆ldual, true) .≈ ([1,10], [1/∆ldual[1]*(1-r), 1/∆ldual[10]*r]))
    @test all(MaxwellFDFD.distweights(c, DUAL, domain, ldual, ∆ldual, false) .≈ ([1,1], [1/∆ldual[1], 0]))

    c = 0  # negative-end boundary
    r = (1-c) / ((1-0) + (120-109))
    @test all(MaxwellFDFD.distweights(c, DUAL, domain, ldual, ∆ldual, true) .≈ ([1,10], [1/∆ldual[1]*(1-r), 1/∆ldual[10]*r]))
    @test all(MaxwellFDFD.distweights(c, DUAL, domain, ldual, ∆ldual, false) .≈ ([1,1], [1/∆ldual[1], 0]))

    c = -0.5  # behind negative-end boundary
    @test_throws ArgumentError MaxwellFDFD.distweights(c, DUAL, domain, ldual, ∆ldual, true)
    @test_throws ArgumentError MaxwellFDFD.distweights(c, DUAL, domain, ldual, ∆ldual, false)

    c = 109  # farthest point affected by positive-end boundary
    @test all(MaxwellFDFD.distweights(c, DUAL, domain, ldual, ∆ldual, true) .≈ ([10,10], [1/∆ldual[10], 0]))
    @test all(MaxwellFDFD.distweights(c, DUAL, domain, ldual, ∆ldual, false) .≈ ([10,10], [1/∆ldual[10], 0]))

    c = 115  # within range affected by positive-end boundary
    r = (c-109) / ((120-109) + (1-0))
    @test all(MaxwellFDFD.distweights(c, DUAL, domain, ldual, ∆ldual, true) .≈ ([10,1], [1/∆ldual[10]*(1-r), 1/∆ldual[1]*r]))
    @test all(MaxwellFDFD.distweights(c, DUAL, domain, ldual, ∆ldual, false) .≈ ([10,10], [1/∆ldual[10], 0]))

    c = 120  # positive-end boundary
    r = (c-109) / ((120-109) + (1-0))
    @test all(MaxwellFDFD.distweights(c, DUAL, domain, ldual, ∆ldual, true) .≈ ([10,1], [1/∆ldual[10]*(1-r), 1/∆ldual[1]*r]))
    @test all(MaxwellFDFD.distweights(c, DUAL, domain, ldual, ∆ldual, false) .≈ ([10,10], [1/∆ldual[10], 0]))

    # Test for N = 1
    domain = [0, 2]
    lprim, ldual = [0], [1]
    ∆lprim, ∆ldual = [9], [10]  # this doesn't make sense geometrically, but for test purpose (they are used to get current density from current)
    c₁, c₂ = 2rand(), 2rand()  # within domain

    ind_c₁, wt_c₁ = MaxwellFDFD.distweights(c₁, PRIM, domain, lprim, ∆lprim, true)
    ind_c₂, wt_c₂ = MaxwellFDFD.distweights(c₂, PRIM, domain, lprim, ∆lprim, true)
    @test ind_c₁ == ind_c₂ == [1,1]
    @test sum(wt_c₁) ≈ sum(wt_c₂) ≈ 1 / ∆lprim[1]

    ind_c₁, wt_c₁ = MaxwellFDFD.distweights(c₁, DUAL, domain, ldual, ∆ldual, true)
    ind_c₂, wt_c₂ = MaxwellFDFD.distweights(c₂, DUAL, domain, ldual, ∆ldual, true)
    @test ind_c₁ == ind_c₂ == [1,1]
    @test sum(wt_c₁) ≈ sum(wt_c₂) ≈ 1 / ∆ldual[1]
end  # @testset "distweights"

@testset "PlaneSrc" begin
    src = PlaneSrc([0,0,1], 0, [1,0,0])  # x-polarized source in z-normal plane

    isbloch = [true, true, true]

    # Coarse grid
    lprim = (-10:10, -10:10, -10:10)
    g3 = Grid(lprim, isbloch)
    ∆a = 1.0^2  # area element

    ft = EE
    boundft = SVec(EE,EE,EE)

    j3d = create_field_array(g3.N)
    add_src!(j3d, ft2gt.(ft,boundft), g3.bounds, g3.l, g3.∆l, g3.isbloch, src)

    @test maximum(abs, j3d) == 1.0  # J = K / ∆z
    @test all(j3d[:,:,:,nY] .== 0)  # Jy = 0
    @test all(j3d[:,:,:,nZ] .== 0)  # Jz = 0

    # Fine grid
    lprim = (-10:0.5:10, -10:0.5:10, -10:0.5:10)
    g3_fine = Grid(lprim, isbloch)

    j3d_fine = create_field_array(g3_fine.N)
    add_src!(j3d_fine, ft2gt.(ft,boundft), g3_fine.bounds, g3_fine.l, g3_fine.∆l, g3_fine.isbloch, src)
    ∆a_fine = 0.5^2  # area element

    @test sum(j3d[1,:,:,nX]) * ∆a ≈ sum(j3d_fine[1,:,:,nX]) * ∆a_fine  # total current through x-normal cross section is independent of grid resolution
    @test maximum(abs, j3d_fine) == 2.0  # J = K / ∆z

    # Nonuniform grid in z
    zprim = sort(rand(21)) * 20
    z_avg = mean(zprim)
    zprim .-= z_avg  # z = 0 is within the z-range
    lprim = (-10:10, -10:10, zprim)
    g3_nu = Grid(lprim, isbloch)

    ∆yprim = g3_nu.∆l[nPR][nY]
    ∆zprim = g3_nu.∆l[nPR][nZ]

    j3d_nu = create_field_array(g3_nu.N)
    add_src!(j3d_nu, ft2gt.(ft,boundft), g3_nu.bounds, g3_nu.l, g3_nu.∆l, g3_nu.isbloch, src)

    @test sum(j3d[1,:,:,nX]) * ∆a ≈ sum(j3d_nu[1,:,:,nX] .* (∆yprim * transpose(∆zprim)))  # total current through x-normal cross section is independent of grid resolution

    # Verify if the total current does not change with the location of the intercept in 3D for Bloch.
end  # @testset "PlaneSrc"

@testset "PointSrc" begin
    src = PointSrc([0.7,0.7,0.7], [1,1,1])  # polarized in [1,1,1] direction

    isbloch = [true, true, true]

    # Coarse grid
    lprim = (-10:10, -10:10, -10:10)
    g3 = Grid(lprim, isbloch)
    ∆v = 1.0^3  # volume element (dipole-normal area element * dipole length)

    ft = EE
    boundft = SVec(EE,EE,EE)

    j3d = create_field_array(g3.N)
    add_src!(j3d, ft2gt.(ft,boundft), g3.bounds, g3.l, g3.∆l, g3.isbloch, src)

    @test sum(j3d[:,:,:,nX] .!= 0) == sum(j3d[:,:,:,nY] .!= 0) == sum(j3d[:,:,:,nZ] .!= 0) == 8
    @test sum(j3d[:,:,:,nX]) * ∆v == src.I∆r * src.p[nX]
    @test sum(j3d[:,:,:,nY]) * ∆v == src.I∆r * src.p[nY]
    @test sum(j3d[:,:,:,nZ]) * ∆v == src.I∆r * src.p[nZ]

    # Fine grid
    lprim = (-10:0.5:10, -10:0.5:10, -10:0.5:10)
    g3_fine = Grid(lprim, isbloch)

    j3d_fine = create_field_array(g3_fine.N)
    add_src!(j3d_fine, ft2gt.(ft,boundft), g3_fine.bounds, g3_fine.l, g3_fine.∆l, g3_fine.isbloch, src)
    ∆v_fine = 0.5^3  # volume element (dipole-normal area element * dipole length)

    @test sum(j3d_fine[:,:,:,nX] .!= 0) == sum(j3d_fine[:,:,:,nY] .!= 0) == sum(j3d_fine[:,:,:,nZ] .!= 0) == 8
    @test sum(j3d_fine[:,:,:,nX]) * ∆v_fine == src.I∆r * src.p[nX]
    @test sum(j3d_fine[:,:,:,nY]) * ∆v_fine == src.I∆r * src.p[nY]
    @test sum(j3d_fine[:,:,:,nZ]) * ∆v_fine == src.I∆r * src.p[nZ]

    # Nonuniform grid in x, y, z
    xprim = sort(rand(21)) * 20; x_avg = mean(xprim); xprim .-= (x_avg-0.7)  # x = 0.7 is within the x-range (average x location is moved to 0.7)
    yprim = sort(rand(21)) * 20; y_avg = mean(yprim); yprim .-= (y_avg-0.7)  # y = 0.7 is within the y-range (average y location is moved to 0.7)
    zprim = sort(rand(21)) * 20; z_avg = mean(zprim); zprim .-= (z_avg-0.7)  # z = 0.7 is within the z-range (average z location is moved to 0.7)
    lprim = (xprim, yprim, zprim)
    g3_nu = Grid(lprim, isbloch)

    ∆xprim = g3_nu.∆l[nPR][nX]
    ∆yprim = g3_nu.∆l[nPR][nY]
    ∆zdual = g3_nu.∆l[nDL][nZ]

    j3d_nu = create_field_array(g3_nu.N)
    add_src!(j3d_nu, ft2gt.(ft,boundft), g3_nu.bounds, g3_nu.l, g3_nu.∆l, g3_nu.isbloch, src)

    Nx, Ny, Nz = g3_nu.N.data
    @test sum(j3d[:,:,:,nZ]) * ∆v ≈ sum(j3d_nu[:,:,:,nZ] .* (reshape(∆xprim,Nx,1,1) .* reshape(∆yprim,1,Ny,1) .* reshape(∆zdual,1,1,Nz)))  # total current through x-normal cross section is independent of grid resolution

    # Verify if the total current dipole strength does not change with the location of the point
    # source in 3D for Bloch.

    # Verify if the total current dipole strength does not change with the z-location in 2D
    # assignment for Bloch (because the current assigned at the positive boundary is added to the
    # negative boundary)
end  # @testset "PlaneSrc"

end  # @testset "source"


text:
module OpenSlide

include("../deps/deps.jl")

export
    OpenSlideImage,
    AssociatedImage,
    open_slide,
    level_dims,
    read_slide,
    properties,
    associated_images,
    read_associated

immutable openslide_t end
typealias OSt Ptr{openslide_t}

type OpenSlideImage
    file::String
    levels::Int
    leveldims
    osptr::Ptr{openslide_t}

    OpenSlideImage(file,levels,dims,osptr) = new(file,int(levels),dims,osptr)
end
# TODO finalizer

type AssociatedImage
    name::String
    dims::Array{Int,1}
    primary::OpenSlideImage
end

################################################################################
include("libopenslide.jl")
################################################################################

function open_slide(file::ASCIIString)
    osptr = openslide_open(file)
    (osptr == C_NULL) && error("Unrecognized slide file: ", file)

    levels = get_level_count(osptr)
    dims = Array{Int,1}[]
    for level in 1:levels
        push!(dims, get_level_dimensions(osptr,level-1))
    end

    return OpenSlideImage(file,levels,dims,osptr)
end

function level_dims(img::OpenSlideImage, level::Int)
    get_level_dimensions(img.osptr, level)
end

function read_slide(img::OpenSlideImage; level = -1, origin = [0,0], extent = None)
    if (level < 0)
        level = 1
    end
    if (extent == None)
        extent = img.leveldims[level]
    end
    read_region(img.osptr, origin..., level, extent...)
end

function properties(img::OpenSlideImage)
    rv = Dict{ASCIIString,ASCIIString}()
    for name in get_property_names(img.osptr)
        rv[name] = get_property_value(img.osptr,name)
    end
    return rv
end

function associated_images(img::OpenSlideImage)
    rv = AssociatedImage[]
    for name in associated_names(s)
        dims = associated_dimensions(name)
        push!(rv, AssociatedImage(name, dims))
    end
end

function read_associated(img::AssociatedImage)
    read_associated_image(img.primary.osptr, img.name, img.dims)
end

end # module


text:
function Mm(subName::String)
    try
        return molarMassData[subName]
    catch
        println("Cannot find $subName")
    end
end

function Mm()
    println(keys(molarMassData))
end


text:
# Duffing_orbit_diagram_vs_period_single_v3.jl
# ============================================

using DrWatson
@quickactivate "NonlinearDynamicsTextbook"
include(srcdir("style.jl"))
using DynamicalSystems, PyPlot, OrdinaryDiffEq
using StaticArrays
using Random
using Statistics

function Duffing_ODE(du,u,p,t)
    d, a, omega = p
    du[1] =  u[2]
    du[2] =  - u[1] - u[1]*u[1]*u[1] - d*u[2] + a*sin(omega*t)
end

function Poincare_orbit(d,a,omega,y_init,ttrans,nattr,pcol)

    p = [d a omega]
    period = 2*pi / omega
    ntrans = ceil( ttrans / period )   # no. of periods required for transient

 # transient to attractor
 # ----------------------
    tspan = (0., ntrans*period)
    prob = ODEProblem(Duffing_ODE, y_init, tspan, p)
    sol = solve(prob, Tsit5(), reltol=1e-8, abstol=1e-8)

  # attractor  with root of Poincare map
  # ---------
    tspan = (0., nattr*period)
    y_init = sol.u[end]
    prob = ODEProblem(Duffing_ODE, y_init, tspan, p)
    sol = solve(prob, Tsit5(), reltol=1e-8, abstol=1e-8, saveat = period/2)

    u1vec = [u[1] for u in sol.u]
    npts = length(u1vec)
    for j = 2:2:npts
        u1vec[j] = -u1vec[j]
    end

  # plot
  # ----
    plt.scatter(period*ones(npts),u1vec,s=2, c = pcol)

    return y_init

end # of function

function plot_orbit(a,d,T_orbit,y_init,delta_T,delta_u1p,color_symbol,color_orbit,linlen,dd)

color_line = COLORS[3]

omega  = 2. * pi / T_orbit
period = T_orbit

# transient to attractor
# ----------------------
ttrans = 100.
ntrans = ceil( ttrans / period )   # no. of periods required for transient
p = [d a omega]
tspan = (0., ntrans*period)
prob = ODEProblem(Duffing_ODE, y_init, tspan, p)
sol = solve(prob, Tsit5(), reltol=1e-8, abstol=1e-8)

# orbit
# -----
nsamp = 300
tspan = (0., nattr*period)
y_init = sol.u[end]
prob = ODEProblem(Duffing_ODE, y_init, tspan, p)
sol = solve(prob, Tsit5(), reltol=1e-8, abstol=1e-8, saveat = period/nsamp)

u1vec = [u[1] for u in sol.u]
u2vec = [u[2] for u in sol.u]
u1p = u1vec[end]

# plot line pointing to orbit
# ---------------------------
xx = [ T_orbit + 0.15delta_T  , T_orbit + linlen*delta_T ]
yy = [ u1p + 0.15delta_u1p , u1p + linlen*delta_u1p ]
plt.plot(xx, yy, c = color_line, linewidth = 1)

# plot symbol indicating current solution
# ---------------------------------------
u1p = u1vec[end]
u2p = u2vec[end]
plt.scatter(T_orbit,u1p,s=55, c = color_symbol, zorder = 3 )

# plot orbit
# ----------
cc = [T_orbit + delta_T, u1p + delta_u1p] # center of orbit
# dd = [0.8 0.2]    # scaling

plt.plot(cc[1].+ dd[1] .* u1vec,  cc[2].+ dd[2] .* u2vec, c = color_orbit, linewidth=1.5)
end

# ----------------------- main program -------------------------

@time begin

# parameters and initial values
# -----------------------------
d = 0.1
a = 7

T_min = 1 # 1
T_max = 29 # 32
n_T = 2000  # 1500
T_vec = [LinRange(T_min, T_max, n_T); LinRange(T_max, T_min, n_T)]

ttrans = 1000 # 150  # minimal transient time
nattr = 4     # no. of points on the attractor
y_init = [0.1 0.1]

# open plot
# ---------
fig = figure(figsize=(0.65*figx,figy))

pcol = COLORS[1]

#  control parameter loop (period)
#  ----------------------
for i_T = 1:length(T_vec)
    period = T_vec[i_T]
    omega = 2*pi / period
    y_init =  Poincare_orbit(d,a,omega,y_init,ttrans,nattr,pcol)
end # of period loop

# tight_layout()
ax = gca()
ax.set_ylim([-2.,1.7])
ax.set_xlim([T_min,T_max])

# Add orbits
# ==========
color_symbol = COLORS[4]
color_orbit = COLORS[4]

# Asymmetric orbit Ia
# ------------------
T_orbit = 6.
y_init = [0.5 -1.]
delta_T = 1.7
delta_u1p = 0.6
linlen = 0.3
dd = [0.8 0.2]  * 0.6
plot_orbit(a,d,T_orbit,y_init,delta_T,delta_u1p,color_symbol,color_orbit,linlen,dd)

# Asymmetric orbit Ib
# ------------------
T_orbit = 6.
y_init = [-1. 1.]
delta_T = 1.2
delta_u1p = -0.8
linlen = 0.6
dd = [0.8 0.2] * 0.6
plot_orbit(a,d,T_orbit,y_init,delta_T,delta_u1p,color_symbol,color_orbit,linlen,dd)

# Asymmetric orbit IIa
# --------------------
T_orbit = 15.
y_init = [1. 1.]
delta_T = -1.5
delta_u1p = 0.8
linlen = 0.6
dd = [0.8 0.2]
plot_orbit(a,d,T_orbit,y_init,delta_T,delta_u1p,color_symbol,color_orbit,linlen,dd)

# Asymmetric orbit IIb
# --------------------
T_orbit = 15.
y_init = [-1. 1.]
delta_T = -1.8
delta_u1p = -0.8
linlen = 0.3
dd = [0.8 0.2]
plot_orbit(a,d,T_orbit,y_init,delta_T,delta_u1p,color_symbol,color_orbit,linlen,dd)

# Symmetric orbit I
# -----------------
color_symbol = COLORS[6]
color_orbit = COLORS[6]

T_orbit = 17.
y_init = [-1. 1.]
delta_T = 1.
delta_u1p = 1.
linlen = 0.5
dd = [0.8 0.2]
plot_orbit(a,d,T_orbit,y_init,delta_T,delta_u1p,color_symbol,color_orbit,linlen,dd)

# Symmetric orbit II
# ------------------
T_orbit = 25.
y_init = [-1. 1.]
delta_T = 0.
delta_u1p = 0.8
linlen = 0.6
dd = [0.8 0.2] * 1.7
plot_orbit(a,d,T_orbit,y_init,delta_T,delta_u1p,color_symbol,color_orbit,linlen,dd)

# Symmetric orbit III
# ------------------
T_orbit = 3.
y_init = [-1. 1.]
delta_T = 0.
delta_u1p = 1.
linlen = 0.4
dd = [0.8 0.2] * 0.5
plot_orbit(a,d,T_orbit,y_init,delta_T,delta_u1p,color_symbol,color_orbit,linlen,dd)

end # of time

xlabel(L"T = 2\pi / \omega")
ylim(-2,2)
yticks([-2,2])
ylabel(L"u_1"; labelpad = -30)

fig.tight_layout(pad=0.3)

wsave(plotsdir("9", "duffing_orbitdiagram"), fig)

max_stars_repo_name: JuliaDynamics/NonlinearDynamicsTextbook

max_issues_repo_name: JuliaDynamics/NonlinearDynamicsTextbook

max_forks_repo_name: JuliaDynamics/NonlinearDynamicsTextbook


text:
################################################################################
#
#   FUNCTIONS FOR PLOTTING OF RECIPROCAL POINTS
#   (overwriting plotSite function for AbstractSite)
#
################################################################################

# PLOTTING RECIPROCAL POINTS IN 2D
function LatPhysPlottingPyPlot.plotSite(
            site    :: S,
            radius  :: Real,
            color   :: Vector{<:Integer}
            ;
            site_labels :: Bool = true,
            use_LaTeX :: Bool = true,
            site_label_fontsize :: Real = 12,
            site_label_offset :: Vector{<:Real} = [0.1, 0.0],
            kwargs...
        ) where {S<:AbstractReciprocalPoint{2}}

    # get the coordinates
    x = point(site)[1]
    y = point(site)[2]
    # scatter the point
    PyPlot.scatter([x,],[y,], color=color./255)
    # maybe annotate the label as text
    if site_labels
        if use_LaTeX
            PyPlot.text(x+site_label_offset[1],y+site_label_offset[2], LaTeXString(L"$"*labelLaTeX(site)*L"$"), fontsize=site_label_fontsize, color=color./255)
        else
            PyPlot.text(x+site_label_offset[1],y+site_label_offset[2], String(label(site)), fontsize=site_label_fontsize, color=color./255)
        end
    end
end

# PLOTTING RECIPROCAL POINTS IN 3D
function LatPhysPlottingPyPlot.plotSite(
            site    :: S,
            radius  :: Real,
            color   :: Vector{<:Integer}
            ;
            site_labels :: Bool = true,
            use_LaTeX :: Bool = true,
            site_label_fontsize :: Real = 12,
            site_label_offset :: Vector{<:Real} = [0.1, 0.0, 0.0],
            kwargs...
        ) where {S<:AbstractReciprocalPoint{3}}

    # get the coordinates
    x = point(site)[1]
    y = point(site)[2]
    z = point(site)[3]
    # scatter the point
    PyPlot.scatter3D([x,],[y,],[z,], color=color./255)
    # maybe annotate the label as text
    if site_labels
        if use_LaTeX
            PyPlot.text3D(x+site_label_offset[1],y+site_label_offset[2],z+site_label_offset[3],  LaTeXString(L"$"*labelLaTeX(site)*L"$"), fontsize=site_label_fontsize, color=color./255)
        else
            PyPlot.text3D(x+site_label_offset[1],y+site_label_offset[2],z+site_label_offset[3],  String(label(site)), fontsize=site_label_fontsize, color=color./255)
        end
    end
end

# export the function
export plotSite


text:
# Axis-aligned multivariate normal (i.e., independent entries, i.e., diagonal covariance matrix)
# using conjugate prior, but not integrating it out. This is just to enable comparison with MVNaaC.
module MVNaaN

module MVNaaNmodel # submodule for component family definitions
export Theta, Data, log_likelihood, log_prior, prior_sample!, new_theta, Theta_clear!, Theta_adjoin!, Theta_remove!,
       Hyperparameters, construct_hyperparameters, update_hyperparameters!, update_parameter!

include("Random.jl")
using .Random

const Data = Array{Float64,1}

type Theta
    mu::Array{Float64,1}     # means
    lambda::Array{Float64,1} # precisions
    d::Int64                 # dimension
    n::Int64                 # number of data points assigned to this cluster
    sum_x::Array{Float64,1}  # sum of the data points x assigned to this cluster
    sum_xx::Array{Float64,1} # sum of x.*x for the data points assigned to this cluster
    Theta(mu,lambda) = (p=new(); p.mu=mu; p.lambda=lambda; p.d=length(mu); p.n=0; p.sum_x=zeros(p.d); p.sum_xx=zeros(p.d); p)
end
new_theta(H) = Theta(zeros(H.d),ones(H.d))

const constant = -0.5*log(2*pi)
Normal_logpdf(x,mu,lambda) = -0.5*lambda*(x-mu)*(x-mu) + 0.5*log(lambda) + constant
MVNaaN_logpdf(x,mu,lambda) = (r=0.0; for i=1:length(x); r+=Normal_logpdf(x[i],mu[i],lambda[i]); end; r)
Theta_clear!(p) = (for i=1:p.d; p.sum_x[i] = 0.0; p.sum_xx[i] = 0.0; end; p.n = 0)
Theta_adjoin!(p,x) = (for i=1:p.d; p.sum_x[i] += x[i]; p.sum_xx[i] += x[i]*x[i]; end; p.n += 1)
Theta_remove!(p,x) = (for i=1:p.d; p.sum_x[i] -= x[i]; p.sum_xx[i] -= x[i]*x[i]; end; p.n -= 1)

type Hyperparameters
    d::Int64    # dimension
    m::Float64  # prior mean of mu's
    c::Float64  # prior precision multiplier for mu's
    a::Float64  # prior shape of lambda's
    b::Float64  # prior rate of lambda's
end

Gamma_logpdf(x,a,b) = (a-1)*log(x) - b*x + a*log(b) - lgamma(a)

log_likelihood(x,p) = MVNaaN_logpdf(x,p.mu,p.lambda)
# prior: Normal(mu[i]|H.m,1/(H.c*lambda[i])) * Gamma(lambda[i]|H.a,H.b)
log_prior(p,m,c,a,b) = (r=0.0; for i=1:p.d; r += Normal_logpdf(p.mu[i],m,c*p.lambda[i]) + Gamma_logpdf(p.lambda[i],a,b); end; r)
log_prior(p,H) = log_prior(p,H.m,H.c,H.a,H.b)
prior_sample!(p,H) = (for i=1:H.d; p.lambda[i] = Random.gamma(H.a,H.b); p.mu[i] = randn()/sqrt(H.c*p.lambda[i]) + H.m; end)

function construct_hyperparameters(options)
    x = options.x
    d = length(x[1])
    mu = mean(x)
    v = mean([xi.*xi for xi in x]) - mu.*mu  # sample variance
    @assert(all(abs.(mu) .< 1e-10) && all(abs.(v - 1.0) .< 1e-10), "Data must be normalized to zero mean, unit variance.")
    m = 0.0
    c = 1.0
    a = 1.0
    b = 1.0
    return Hyperparameters(d,m,c,a,b)
end

function update_parameter!(theta_a,theta_b,H,active,density)
    n,sum_x,sum_xx = theta_b.n,theta_b.sum_x,theta_b.sum_xx
    logp = 0.0
    A = H.a + 0.5*(n+1)
    for i = 1:H.d
        # update mean
        L = theta_a.lambda[i]*(n + H.c)
        M = (sum_x[i] + H.c*H.m) / (n + H.c)
        if active; theta_b.mu[i] = randn()/sqrt(L) + M; end
        if density; logp += Normal_logpdf(theta_b.mu[i],M,L); end

        # update precision
        B = H.b + 0.5*(sum_xx[i] - 2*theta_b.mu[i]*sum_x[i] + n*theta_b.mu[i]^2) + 0.5*H.c*(theta_b.mu[i] - H.m)^2
        if active; theta_b.lambda[i] = Random.gamma(A,B); end
        if density; logp += Gamma_logpdf(theta_b.lambda[i],A,B); end
    end
    return (density? logp : NaN)
end

function update_hyperparameters!(H,theta,list,t,x,z)
    #error("update_hyperparameters is not yet implemented.")
end

end # module MVNaaNmodel
using .MVNaaNmodel

# Include generic code
include("generic.jl")

# Include core sampler code
include("coreNonconjugate.jl")

end # module MVNaaN


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct ApplicationGatewayPathRule <: SwaggerModel
    id::Any # spec type: Union{ Nothing, String } # spec name: id

    function ApplicationGatewayPathRule(;id=nothing)
        o = new()
        validate_property(ApplicationGatewayPathRule, Symbol("id"), id)
        setfield!(o, Symbol("id"), id)
        o
    end
end # type ApplicationGatewayPathRule

const _property_map_ApplicationGatewayPathRule = Dict{Symbol,Symbol}(Symbol("id")=>Symbol("id"))
const _property_types_ApplicationGatewayPathRule = Dict{Symbol,String}(Symbol("id")=>"String")
Base.propertynames(::Type{ ApplicationGatewayPathRule }) = collect(keys(_property_map_ApplicationGatewayPathRule))
Swagger.property_type(::Type{ ApplicationGatewayPathRule }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_ApplicationGatewayPathRule[name]))}
Swagger.field_name(::Type{ ApplicationGatewayPathRule }, property_name::Symbol) =  _property_map_ApplicationGatewayPathRule[property_name]

function check_required(o::ApplicationGatewayPathRule)
    true
end

function validate_property(::Type{ ApplicationGatewayPathRule }, name::Symbol, val)
end


text:
import TimeZones
const TIME_ZONE = TimeZones.TimeZone("America/New_York")


text:
import LinearAlgebra: lmul!, rmul!
import Base: Matrix

function compute_rotate(f, g)
    r = hypot(f, g)
    cs = f / r
    sn = g / r
    if cs < 0 && abs(f) > abs(g)
        cs = -cs
        sn = -sn
        r = -r
    end
    return cs, sn, r
end

# Some utility to materialize a rotation to a matrix.
function Matrix(r::Rotation2{Tc,Ts}, n::Int) where {Tc,Ts}
    r.i < n || throw(ArgumentError("Matrix should have order $(r.i+1) or larger"))
    G = Matrix{promote_type(Tc,Ts)}(I, n, n)
    G[r.i+0,r.i+0] = r.c
    G[r.i+1,r.i+0] = -conj(r.s)
    G[r.i+0,r.i+1] = r.s
    G[r.i+1,r.i+1] = r.c
    return G
end

function Matrix(r::Rotation3{Tc,Ts}, n::Int) where {Tc,Ts}
    G₁ = Matrix(Rotation2(r.c₁, r.s₁, r.i + 1), n)
    G₂ = Matrix(Rotation2(r.c₂, r.s₂, r.i), n)
    return G₂ * G₁
end

"""
Get a rotation that maps [p₁, p₂] to a multiple of [1, 0]
"""
function get_rotation(p₁, p₂, i::Int)
    c, s, nrm = compute_rotate(p₁, p₂)
    Rotation2(c, s, i), nrm
end

"""
Get a rotation that maps [p₁, p₂, p₃] to a multiple of [1, 0, 0]
"""
function get_rotation(p₁, p₂, p₃, i::Int)
    c₁, s₁, nrm₁ = compute_rotate(p₂, p₃)
    c₂, s₂, nrm₂ = compute_rotate(p₁, nrm₁)
    Rotation3(c₁, s₁, c₂, s₂, i), nrm₂
end

lmul!(::SmallRotation, ::NotWanted, args...) = nothing
rmul!(::NotWanted, ::SmallRotation, args...) = nothing

@inline function lmul!(G::Rotation3, A::AbstractMatrix)
    @inbounds for j = axes(A, 2)
        a₁ = A[G.i+0,j]
        a₂ = A[G.i+1,j]
        a₃ = A[G.i+2,j]

        a₂′ = G.c₁ * a₂ + G.s₁ * a₃
        a₃′ = -G.s₁' * a₂ + G.c₁ * a₃

        a₁′′ = G.c₂ * a₁ + G.s₂ * a₂′
        a₂′′ = -G.s₂' * a₁ + G.c₂ * a₂′

        A[G.i+0,j] = a₁′′
        A[G.i+1,j] = a₂′′
        A[G.i+2,j] = a₃′
    end

    A
end

@inline function rmul!(A::AbstractMatrix, G::Rotation3)
    @inbounds for j = axes(A, 1)
        a₁ = A[j,G.i+0]
        a₂ = A[j,G.i+1]
        a₃ = A[j,G.i+2]

        a₂′ = a₂ * G.c₁ + a₃ * G.s₁'
        a₃′ = a₂ * -G.s₁ + a₃ * G.c₁

        a₁′′ = a₁ * G.c₂ + a₂′ * G.s₂'
        a₂′′ = a₁ * -G.s₂ + a₂′ * G.c₂

        A[j,G.i+0] = a₁′′
        A[j,G.i+1] = a₂′′
        A[j,G.i+2] = a₃′
    end
    A
end

@inline function lmul!(G::Rotation2, A::AbstractMatrix)
    @inbounds for j = axes(A, 2)
        a₁ = A[G.i+0,j]
        a₂ = A[G.i+1,j]

        a₁′ = G.c * a₁ + G.s * a₂
        a₂′ = -G.s' * a₁ + G.c * a₂

        A[G.i+0,j] = a₁′
        A[G.i+1,j] = a₂′
    end

    A
end

@inline function rmul!(A::AbstractMatrix, G::Rotation2)
    @inbounds for j = axes(A, 1)
        a₁ = A[j,G.i+0]
        a₂ = A[j,G.i+1]

        a₁′ = a₁ * G.c + a₂ * G.s'
        a₂′ = a₁ * -G.s + a₂ * G.c

        A[j,G.i+0] = a₁′
        A[j,G.i+1] = a₂′
    end
    A
end


text:
struct Input{T}
    year::Vector{Int}
    month::Vector{Int}
    day::Vector{Int}
    hour::Vector{Int}
    SW::Vector{T}
    LW::Vector{T}
    Sf::Vector{T}
    Rf::Vector{T}
    Ta::Vector{T}
    RH::Vector{T}
    Ua::Vector{T}
    Ps::Vector{T}
end

Base.@kwdef struct Constants{T}

    cp::T = 1005          # Specific heat capacity of dry air (J/K/kg)
    eps::T = 0.622       # Ratio of molecular weights of water and dry air    BEFORE eps
    e0::T = 611.213       # Saturation vapour pressure at Tm (Pa)
    g::T = 9.81           # Acceleration due to gravity (m/s^2)
    hcap_ice::T = 2100.   # Specific heat capacity of ice (J/K/kg)
    hcap_wat::T = 4180.   # Specific heat capacity of water (J/K/kg)
    hcon_air::T = 0.025   # Thermal conductivity of air (W/m/K)
    hcon_clay::T = 1.16   # Thermal conductivity of clay (W/m/K)
    hcon_ice::T = 2.24    # Thermal conducivity of ice (W/m/K)
    hcon_sand::T = 1.57   # Thermal conductivity of sand (W/m/K)
    hcon_wat::T = 0.56    # Thermal conductivity of water (W/m/K)
    Lc::T = 2.501e6       # Latent heat of condensation (J/kg)
    Lf::T = 0.334e6       # Latent heat of fusion (J/kg)
    Ls::T = Lc + Lf       # Latent heat of sublimation (J/kg)
    R::T = 8.3145         # Universal gas constant (J/K/mol)
    Rgas::T = 287         # Gas constant for dry air (J/K/kg)
    Rwat::T = 462         # Gas constant for water vapour (J/K/kg)
    rho_ice::T = 917.     # Density of ice (kg/m^3)
    rho_wat::T = 1000.    # Density of water (kg/m^3)
    sb::T = 5.67e-8       # Stefan-Boltzmann constant (W/m^2/K^4)
    Tm::T = 273.15        # Melting point (K)
    vkman::T = 0.4        # Von Karman constant

end

Base.@kwdef mutable struct EBM{T}

    ### Settings ##############################################################

    # Maximum number of snow layers
    Nsmax::Int = 3

    # Number of soil layers
    Nsoil::Int = 4

    # Time step (s)
    dt::T = 3600

    # Minimum snow layer thicknesses (m)
    Dzsnow::Vector{T} = [0.1, 0.2, 0.4]

    # Soil layer thicknesses (m)
    Dzsoil::Vector{T} = [0.1, 0.2, 0.4, 0.8]

    # Temperature measurement height (m)
    zT::T = 2

    # Wind measurement height (m)
    zU::T = 10

    # Subtract snow depth from measurement height
    zvar::Bool = true

    ### Snow parameters #######################################################

    # Maximum albedo for fresh snow
    asmx::T = 0.8

    # Minimum albedo for melting snow
    asmn::T = 0.5

    # Stability slope parameter
    bstb::T = 5

    # Snow thermal conductivity exponent
    bthr::T = 2

    # Snow cover fraction depth scale (m)
    hfsn::T = 0.1

    # Fixed thermal conductivity of snow (W/m/K)
    kfix::T = 0.24

    # Fixed snow density (kg/m^3)
    rho0::T = 300

    # Fresh snow density (kg/m^3)
    rhof::T = 100

    # Maximum density for cold snow (kg/m^3)
    rcld::T = 300

    # Maximum density for melting snow (kg/m^3)
    rmlt::T = 500

    # Snowfall to refresh albedo (kg/m^2)
    Salb::T = 10

    # Albedo decay temperature threshold (C)
    Talb::T = -2

    # Cold snow albedo decay timescale (h)
    tcld::T = 1000

    # Melting snow albedo decay timescale (h)
    tmlt::T = 100

    # Snow compaction time scale (h)
    trho::T = 200

    # Irreducible liquid water content of snow
    Wirr::T = 0.03

    # Snow roughness length (m)
    z0sn::T = 0.01

    ### Surface parameters ####################################################

    # Snow-free ground albedo
    alb0::T = 0.2

    # Surface conductance for saturated soil (m/s)
    gsat::T = 0.01

    # Snow-free roughness length (m)
    z0sf::T = 0.1

    ### Soil properties #######################################################

    # Soil clay fraction
    fcly::T = 0.3

    # Soil sand fraction
    fsnd::T = 0.6

    # Clapp-Hornberger exponent
    b::T = 7.630000000000001

    # Volumetric heat capacity of dry soil (J/K/m^3)
    hcap_soil::T = 2.2993333333333335e6

    # Saturated soil water pressure (m)
    sathh::T = 0.10789467222298288

    # Volumetric soil moisture concentration at saturation
    Vsat::T = 0.4087

    # Volumetric soil moisture concentration at critical point
    Vcrit::T = 0.2603859120641063

    # Thermal conductivity of dry soil (W/m/K)
    hcon_soil::T = 0.2740041303112452

    ### State variables #######################################################

    # Snow albedo
    albs::T = 0.8

    # Snow layer thicknesses (m)
    Ds::Vector{T} = fill(0, Nsmax)

    # Number of snow layers
    Nsnow::Int = 0

    # Ice content of snow layers (kg/m^2)
    Sice::Vector{T} = fill(0, Nsmax)

    # Liquid content of snow layers (kg/m^2)
    Sliq::Vector{T} = fill(0, Nsmax)

    # Volumetric moisture content of soil layers
    theta::Vector{T} = fill(0.5 * Vsat, Nsoil)

    # Snow layer temperatures (K)
    Tsnow::Vector{T} = fill(273.15, Nsmax)

    # Soil layer temperatures (K)
    Tsoil::Vector{T} = fill(285, Nsoil)

    # Surface skin temperature (K)
    Tsurf::T = Tsoil[1]

    ### Configurations ########################################################

    # Snow albedo model
    am::Int = 0

    # Snow conductivity model
    cm::Int = 0

    # Snow density model
    dm::Int = 0

    # Surface exchange model
    em::Int = 0

    # Snow hydraulics model
    hm::Int = 0

    ### Internal variables ####################################################

    # Thermal conductivity of snow (W/m/K)
    ksnow::Vector{T} = fill(0, Nsmax)

    # Thermal conductivity of soil (W/m/K)
    ksoil::Vector{T} = fill(0, Nsoil)

    # Areal heat capacity of soil (J/K/m^2)
    csoil::Vector{T} = fill(0, Nsoil)

    # Surface moisture conductance (m/s)
    gs::T = 0

    # Effective albedo
    alb::T = 0

    # Fresh snow density (kg/m^3)
    rfs::T = 0

    # Snow cover fraction
    fsnow::T = 0

    # Transfer coefficient for heat and moisture
    CH::T = 0

    # Surface roughness length (m)
    z0::T = 0

    # Surface thermal conductivity (W/m/K)
    ksurf::T = 0

    # Surface layer temperature (K)
    Ts1::T = 0

    # Surface layer thickness (m)
    Dz1::T = 0

    # Snow sublimation rate (kg/m^2/s)
    Esnow::T = 0

    # Heat flux into surface (W/m^2)
    Gsurf::T = 0

    # Sensible heat flux (W/m^2)
    Hsurf::T = 0

    # Latent heat flux (W/m^2)
    LEsrf::T = 0

    # Surface melt rate (kg/m^2/s)
    Melt::T = 0

    # Net radiation (W/m^2)
    Rnet::T = 0

    # Heat flux into soil (W/m^2)
    Gsoil::T = 0

    # Temporary soil variables
    soil_a::Vector{T} = fill(0, Nsoil)
    soil_b::Vector{T} = fill(0, Nsoil)
    soil_c::Vector{T} = fill(0, Nsoil)
    soil_Gs::Vector{T} = fill(0, Nsoil)
    soil_rhs::Vector{T} = fill(0, Nsoil)
    soil_dTs::Vector{T} = fill(0, Nsoil)
    soil_gamma::Vector{T} = fill(0, Nsoil)

    # Temporary snow variables
    snow_csnow::Vector{T} = fill(0, Nsmax)
    snow_E::Vector{T} = fill(0, Nsmax)
    snow_U::Vector{T} = fill(0, Nsmax)
    snow_Gs::Vector{T} = fill(0, Nsmax)
    snow_dTs::Vector{T} = fill(0, Nsmax)
    snow_a::Vector{T} = fill(0, Nsmax)
    snow_b::Vector{T} = fill(0, Nsmax)
    snow_c::Vector{T} = fill(0, Nsmax)
    snow_rhs::Vector{T} = fill(0, Nsmax)
    snow_gamma::Vector{T} = fill(0, Nsmax)

end


text:
function run_life_view!(results::Array{Float64,3},environment::Array{Float64,3},
    grn::NamedTuple,diff_adhesion::Array{Float64,2},diffusion_constants::Array{Float64,1},
    Ax::Array{Float64,2},Ay::Array{Float64,2},
    dc::Array{Float64,1},cell_config::CellularConfig,exp_config::ExperimentalConfig,grn_config::FGRNConfig)

    target_cell_areas = Dict{Float64,Float64}()
    n_div = Dict{Float64,Integer}()

    all_cells::Array{Float64,1} = [1.]
    n_div[1.] = 0
    target_cell_areas[1.] = 2*0.25*cell_config.cell_size_upper_limit

    if exp_config.est_strength == "evolve"
        affinity_1 = calculate_affinity(grn.rec_prom[:,:,1],grn_config.substrate_protein)
        affinity_2 = calculate_affinity(grn.rec_prom[:,:,2],grn_config.substrate_protein)
        affinity_3 = calculate_affinity(grn.rec_prom[:,:,3],grn_config.substrate_protein)
        affinity_4 = calculate_affinity(grn.rec_prom[:,:,4],grn_config.substrate_protein)
    else
        affinity_1 = 1 - exp_config.est_strength
        affinity_2 = 1 - exp_config.est_strength
        affinity_3 = 1 - exp_config.est_strength
        affinity_4 = 1 - exp_config.est_strength
    end

    left_segment = Int.(floor(cell_config.grid_size/3))
    middle_segment = Int.(floor(2*cell_config.grid_size/3))

    est_setup = ESTSetup(affinity_1,affinity_2,affinity_3,affinity_4,left_segment,middle_segment)

    for step in 1:size(results,3) # development

        results[:,:,step] = environment[:,:,2]

        environment_update2!(grn,environment,all_cells,diff_adhesion,target_cell_areas,cell_config,exp_config)

        if size(all_cells,1) == 0
            break
        end

        cell_mitosis_and_apoptosis!(environment,all_cells,n_div,cell_config)

        protein_dynamics!(environment,est_setup,grn,all_cells,target_cell_areas,diffusion_constants,Ax,Ay,dc,
            exp_config,cell_config,grn_config)

        for cell in all_cells
            if n_div[cell] > cell_config.max_split
                target_cell_areas[cell] = 0
            end
        end

    end
end

function run_life_view!(results::Array{Float64,4},environment::Array{Float64,3},
    grn::NamedTuple,diff_adhesion::Array{Float64,2}, diffusion_constants::Array{Float64,1},
    Ax::Array{Float64,2},Ay::Array{Float64,2},dc::Array{Float64,1},cell_config::CellularConfig,exp_config::ExperimentalConfig,grn_config::FGRNConfig)

    target_cell_areas = Dict{Float64,Float64}()
    n_div = Dict{Float64,Integer}()

    all_cells::Array{Float64,1} = [1.]
    n_div[1.] = 0
    target_cell_areas[1.] = 2*0.25*cell_config.cell_size_upper_limit

    if exp_config.est_strength == "evolve"
        affinity_1 = calculate_affinity(grn.rec_prom[:,:,1],grn_config.substrate_protein)
        affinity_2 = calculate_affinity(grn.rec_prom[:,:,2],grn_config.substrate_protein)
        affinity_3 = calculate_affinity(grn.rec_prom[:,:,3],grn_config.substrate_protein)
        affinity_4 = calculate_affinity(grn.rec_prom[:,:,4],grn_config.substrate_protein)
    else
        affinity_1 = 1 - exp_config.est_strength
        affinity_2 = 1 - exp_config.est_strength
        affinity_3 = 1 - exp_config.est_strength
        affinity_4 = 1 - exp_config.est_strength
    end

    left_segment = Int.(floor(cell_config.grid_size/3))
    middle_segment = Int.(floor(2*cell_config.grid_size/3))

    est_setup = ESTSetup(affinity_1,affinity_2,affinity_3,affinity_4,left_segment,middle_segment)

    for step in 1:size(results,4)

        results[:,:,:,step] = environment[:,:,:]

        environment_update2!(grn,environment,all_cells,diff_adhesion,target_cell_areas,cell_config,exp_config)

        if size(all_cells,1) == 0
            break
        end

        cell_mitosis_and_apoptosis!(environment,all_cells,n_div,cell_config)

        protein_dynamics!(environment,est_setup,grn,all_cells,target_cell_areas,diffusion_constants,Ax,Ay,dc,
            exp_config,cell_config,grn_config)

        for cell in all_cells
            if n_div[cell] > cell_config.max_split
                target_cell_areas[cell] = 0
            end
        end

    end
end

function run_life_view!(damage_times::Vector{Int64},damage_percentage::Float64,results::Array{Float64,3},environment::Array{Float64,3},
    grn::NamedTuple,diff_adhesion::Array{Float64,2}, diffusion_constants::Array{Float64,1},
    Ax::Array{Float64,2},Ay::Array{Float64,2},dc::Array{Float64,1},
    cell_config::CellularConfig,exp_config::ExperimentalConfig,grn_config::FGRNConfig)

    target_cell_areas = Dict{Float64,Float64}()
    n_div = Dict{Float64,Integer}()

    all_cells::Array{Float64,1} = [1.]
    n_div[1.] = 0
    target_cell_areas[1.] = 2*0.25*cell_config.cell_size_upper_limit

    if exp_config.est_strength == "evolve"
        affinity_1 = calculate_affinity(grn.rec_prom[:,:,1],grn_config.substrate_protein)
        affinity_2 = calculate_affinity(grn.rec_prom[:,:,2],grn_config.substrate_protein)
        affinity_3 = calculate_affinity(grn.rec_prom[:,:,3],grn_config.substrate_protein)
        affinity_4 = calculate_affinity(grn.rec_prom[:,:,4],grn_config.substrate_protein)
    else
        affinity_1 = 1 - exp_config.est_strength
        affinity_2 = 1 - exp_config.est_strength
        affinity_3 = 1 - exp_config.est_strength
        affinity_4 = 1 - exp_config.est_strength
    end

    left_segment = Int.(floor(cell_config.grid_size/3))
    middle_segment = Int.(floor(2*cell_config.grid_size/3))

    est_setup = ESTSetup(affinity_1,affinity_2,affinity_3,affinity_4,left_segment,middle_segment)

    for step in 1:size(results,3)

        if step in damage_times
            d_int = Int64.(floor(damage_percentage*size(all_cells,1)))
            damage_pixels = sample(all_cells,d_int,replace = false)
            for dp in damage_pixels
                environment[environment[:,:,1] .== dp,:] .= 0
            end

        end

        results[:,:,step] = environment[:,:,2]

        environment_update2!(grn,environment,all_cells,diff_adhesion,target_cell_areas,cell_config,exp_config)

        if size(all_cells,1) == 0
            break
        end

        cell_mitosis_and_apoptosis!(environment,all_cells,n_div,cell_config)

        protein_dynamics!(environment,est_setup,grn,all_cells,target_cell_areas,diffusion_constants,Ax,Ay,dc,
            exp_config,cell_config,grn_config)

        for cell in all_cells
            if n_div[cell] > cell_config.max_split
                target_cell_areas[cell] = 0
            end
        end

    end
end

function run_life_view!(damage_times::Vector{Int64},damage_percentage::Float64,energies::Array{Float64,1},
    results::Array{Float64,3},environment::Array{Float64,3},
    grn::NamedTuple,diff_adhesion::Array{Float64,2}, diffusion_constants::Array{Float64,1},
    Ax::Array{Float64,2},Ay::Array{Float64,2},dc::Array{Float64,1},
    cell_config::CellularConfig,exp_config::ExperimentalConfig,grn_config::FGRNConfig)

    target_cell_areas = Dict{Float64,Float64}()
    n_div = Dict{Float64,Integer}()

    all_cells::Array{Float64,1} = [1.]
    n_div[1.] = 0
    target_cell_areas[1.] = 2*0.25*cell_config.cell_size_upper_limit

    if exp_config.est_strength == "evolve"
        affinity_1 = calculate_affinity(grn.rec_prom[:,:,1],grn_config.substrate_protein)
        affinity_2 = calculate_affinity(grn.rec_prom[:,:,2],grn_config.substrate_protein)
        affinity_3 = calculate_affinity(grn.rec_prom[:,:,3],grn_config.substrate_protein)
        affinity_4 = calculate_affinity(grn.rec_prom[:,:,4],grn_config.substrate_protein)
    else
        affinity_1 = 1 - exp_config.est_strength
        affinity_2 = 1 - exp_config.est_strength
        affinity_3 = 1 - exp_config.est_strength
        affinity_4 = 1 - exp_config.est_strength
    end

    left_segment = Int.(floor(cell_config.grid_size/3))
    middle_segment = Int.(floor(2*cell_config.grid_size/3))

    est_setup = ESTSetup(affinity_1,affinity_2,affinity_3,affinity_4,left_segment,middle_segment)

    for step in 1:size(results,3)

        if step in damage_times
            d_int = Int64.(floor(damage_percentage*size(all_cells,1)))
            damage_pixels = sample(all_cells,d_int,replace = false)
            for dp in damage_pixels
                environment[environment[:,:,1] .== dp,:] .= 0
            end

        end

        results[:,:,step] = environment[:,:,2]

        energies[step] = measure_energy(grn,environment,all_cells,diff_adhesion,target_cell_areas,cell_config,exp_config)

        environment_update2!(grn,environment,all_cells,diff_adhesion,target_cell_areas,cell_config,exp_config)

        if size(all_cells,1) == 0
            break
        end

        cell_mitosis_and_apoptosis!(environment,all_cells,n_div,cell_config)

        protein_dynamics!(environment,est_setup,grn,all_cells,target_cell_areas,diffusion_constants,Ax,Ay,dc,
            exp_config,cell_config,grn_config)

        for cell in all_cells
            if n_div[cell] > cell_config.max_split
                target_cell_areas[cell] = 0
            end
        end

    end
end

function run_life_fitness!(time_steps::Int64,environment::Array{Float64,3},
    all_cells::Array{Float64,1},target_cell_areas::Dict{Float64,Float64},n_div::Dict{Float64,Integer},
    grn::NamedTuple,diff_adhesion::Array{Float64,2}, diffusion_constants::Array{Float64,1},
    Ax::Array{Float64,2},Ay::Array{Float64,2},dc::Array{Float64,1},cell_config::CellularConfig,
    exp_config::ExperimentalConfig,grn_config::FGRNConfig)

    if exp_config.est_strength == "evolve"
        affinity_1 = calculate_affinity(grn.rec_prom[:,:,1],grn_config.substrate_protein)
        affinity_2 = calculate_affinity(grn.rec_prom[:,:,2],grn_config.substrate_protein)
        affinity_3 = calculate_affinity(grn.rec_prom[:,:,3],grn_config.substrate_protein)
        affinity_4 = calculate_affinity(grn.rec_prom[:,:,4],grn_config.substrate_protein)
    else
        affinity_1 = 1 - exp_config.est_strength
        affinity_2 = 1 - exp_config.est_strength
        affinity_3 = 1 - exp_config.est_strength
        affinity_4 = 1 - exp_config.est_strength
    end

    left_segment = Int.(floor(cell_config.grid_size/3))
    middle_segment = Int.(floor(2*cell_config.grid_size/3))

    est_setup = ESTSetup(affinity_1,affinity_2,affinity_3,affinity_4,left_segment,middle_segment)

    for step in 1:time_steps

        environment_update2!(grn,environment,all_cells,diff_adhesion,target_cell_areas,cell_config,exp_config)

        if size(all_cells,1) == 0
            break
        end

        cell_mitosis_and_apoptosis!(environment,all_cells,n_div,cell_config)

        protein_dynamics!(environment,est_setup,grn,all_cells,target_cell_areas,diffusion_constants,Ax,Ay,dc,
            exp_config,cell_config,grn_config)

        for cell in all_cells
            if n_div[cell] > cell_config.max_split
                target_cell_areas[cell] = 0
            end
        end

    end
end

max_stars_repo_name: harrybooth/FractalDevelopment

max_issues_repo_name: harrybooth/FractalDevelopment

max_forks_repo_name: harrybooth/FractalDevelopment


text:
const ut3 = UpperTriangular(eye(3))
const ut2 = UpperTriangular(eye(2))

@test_throws DimensionMismatch MixedModels.inject!(ut3,ut2)
@test MixedModels.inject!(UpperTriangular(zeros(3,3)),ut3) == eye(3)

speye41 = sparse(vcat(1:4, 1), vcat(1:4, 2), ones(5))
@test nonzeros(MixedModels.inject!(speye41, speye4)) == [1., 0, 1, 1, 1]


text:
module FromFile

export @from

using Requires

track(mod, path) = nothing

track_modules() = isinteractive()

function __init__()
    @require Revise="295af30f-e4ad-537b-8983-00126c2a3abe" track(mod, path) = Revise.track(mod, path)
end

macro from(path::String, ex::Expr)
    esc(from_m(__module__, __source__, path, ex))
end

macro from(path::String)
    esc(from_m(__module__, __source__, path, Expr(:block)))
end

function from_m(m::Module, s::LineNumberNode, path::String, root_ex::Expr)
    import_exs = if root_ex.head === :block
        filter(ex -> !(ex isa LineNumberNode), root_ex.args)
    else
        [root_ex]
    end

    all(ex -> ex.head === :using || ex.head === :import, import_exs) || error("expected using/import statement")

    root = Base.moduleroot(m)
    basepath = dirname(String(s.file))

    # file path should always be relative to the
    # module loads it, unless specified as absolute
    # path or the module is created interactively
    if !isabspath(path) && basepath != ""
        path = joinpath(basepath, path)
    end
    path = abspath(path)

    if root === Main
        file_module_sym = Symbol(path)
    else
        file_module_sym = Symbol(relpath(path, pathof(root)))
    end

    if isdefined(root, file_module_sym)
        file_module = getfield(root, file_module_sym)
    else
        file_module = Base.eval(root, :(module $(file_module_sym); include($path); end))

        # In interactive sessions, track generated module using Revise.jl if Revise has been loaded
        track_modules() && track(file_module, path)
    end

    return Expr(:block, map(import_exs) do ex
        loading = Expr(ex.head)

        for each in ex.args
            each isa Expr || continue

            if each.head === :(:) # using/import A: a, b, c
                each.args[1].args[1] === :(.) && error("cannot load relative module from file")
                push!(loading.args, Expr(:(:), Expr(:., fullname(file_module)..., each.args[1].args...), each.args[2:end]...) )
            elseif each.head === :(.) # using/import A, B.C
                each.args[1] === :(.) && error("cannot load relative module from file")
                push!(loading.args, Expr(:., fullname(file_module)..., each.args...))
            else
                error("invalid syntax $ex")
            end
        end
        return loading
    end...)
end

end


text:
# Make this work
# 1. Need a good example
# 2. Plot both preds in plotscalar
# 3. Isssue is this is whatif we get into a satisfying region?

# May be more than one gamma per constraint

function softeqgamma(γ0 = 1.0)
  γ = Ref(γ0)
  function op(x, y, k = Omega.globalkernel())
    SoftBool(k(Omega.d(x, y) * (γ.x)))
  end
  (op = lift(op), γ = γ)
end

# Moving Averages
"Exponential moving average"
function ema(α)

end

# Callbacks

struct EMA
  vals::Array{Float64}
end

function Base.push!(ema::EMA, val)
  push!(ema.vals, val)
end

average(ema::EMA) = Statistics.mean(ema.vals)

"Callback that updates `γs` to minimize the `var([E[pred] for pred in preds])`"
function updatescales(preds, γs, windowsize)
  @pre length(γs) == length(preds) + 1
  weights = Array{Float64}(undef, length(preds), length(windowsize))
  emas = [EMA() for i = 1:length(preds)]
  function updateγs(data, IterEnd)
    predvals = applymany(preds, ω)
    for pred in predvals

      push()
    end
    diffs = averages[2:end] .- averages[1]

  end
end


text:
function ElevationExtremumDetector(arg0::TopocentricFrame)
    return ElevationExtremumDetector((TopocentricFrame,), arg0)
end

function ElevationExtremumDetector(arg0::jdouble, arg1::jdouble, arg2::TopocentricFrame)
    return ElevationExtremumDetector((jdouble, jdouble, TopocentricFrame), arg0, arg1, arg2)
end

function g(obj::ElevationExtremumDetector, arg0::SpacecraftState)
    return jcall(obj, "g", jdouble, (SpacecraftState,), arg0)
end

function get_elevation(obj::ElevationExtremumDetector, arg0::SpacecraftState)
    return jcall(obj, "getElevation", jdouble, (SpacecraftState,), arg0)
end

function get_topocentric_frame(obj::ElevationExtremumDetector)
    return jcall(obj, "getTopocentricFrame", TopocentricFrame, ())
end


text:
include("./exportData.jl")
using .exportData, DelimitedFiles

inds = [78, 88, 100, 106, 114]

output = Array{Any,2}(undef,(1002,length(inds)+1));
for (n,ind) in enumerate(inds)
	t, data, BCL = getExpData(ind; tInds=1:1000);
	if n==1
		output[1,:] .= vcat("#ind", inds)
		output[2,1] = "#BCLs"
		output[3:end,1] .= t.-t[1]
	end
	output[2,n+1] = round(BCL; digits=2)
	output[3:end,n+1] .= data[:]
end

open("./fitting_data.txt", "w") do io
	writedlm(io, output)
end

output = Array{Any,2}(undef,(length(inds)+1,18))
for (n,ind) in enumerate(inds)
	P = readdlm("./pde/bbo/$(ind)_params.txt")
	if n==1
		output[n,1] = "#ind"
		output[n,2:end] .= P[1,:]
	end
	output[n+1,1] = ind
	output[n+1,2:end] .= round.(P[end,:]; digits=2)
end

open("./parameters.txt", "w") do io
	writedlm(io, output)
end


text:
function LGCPobs(field, data::AbstractVector{Tuple{Int, Int}}, A)
    # The first sum is over the data points and the second is over the grid.
    logA = log(A)
    s = sum((field[i] + logA)*n for (i, n) in data)
    #s - data[:gridcellsize]*sum(exp, field)
    s - A*sum(exp, field)
end

function PPobs(intensity, data::AbstractVector{Tuple{Int, Int}}, A)
    # The first sum is over the data points and the second is over the grid.
    logA = log(A)
    s = sum((log(intensity[i]) + logA)*n for (i, n) in data)
    #s - data[:gridcellsize]*sum(exp, field)
    s - A*sum(intensity)
end

function Distributions.logpdf(d::Union{LogGaussianCoxProcess, PoissonProcess}, di::Discretization, x)
    logpdf_(d, di, x)[1]
end

function logpdf_(d, di::Discretization, x::PointPattern)
    logpdf_(d, di, cellcounts(Vector{Tuple{Int,Int}}, di, x))
end
function logpdf_(d::PoissonProcess, di::Discretization, x::AbstractVector{Tuple{Int, Int}})
    issquare(di) || throw(ArgumentError("Only square discretizations supported."))
    d.window == di.window && size(d.intensity) == di.size || throw(ArgumentError("Incompatible discretization."))
    scale = step(di.gridx)
    l1 = PPobs(d.intensity, x, scale^2)
    l2 = sum(((cell, count),) -> logabsgamma(count+1)[1], x)
    l1 - l2
end
function logpdf_(d::LogGaussianCoxProcess, di::Discretization, x::AbstractVector{Tuple{Int, Int}})
    issquare(di) || throw(ArgumentError("Only square discretizations supported."))
    d.range > 0 || return -Inf
    d.sd > 0 || return -Inf
    nu = d.nu
    sd = d.sd
    scale = step(di.gridx)
    range = d.range / scale # range in cells
    fixedfield = d.mean
    prec = 1/sd^2
    (l, _, pz), random_field_mode, H = poisson_gp(x, fixedfield, prec, range, nu, scale^2, di.size)
    lr = -Inf
    if (l != -Inf)
        l1 = LGCPobs(fixedfield .+ random_field_mode, x, scale^2)
        l2 = sum(((cell, count),) -> logabsgamma(count+1)[1], x)
        lr = pz + l1 - l2
        l = l + l1 - l2
    end
    l, random_field_mode, H, lr
end

function logpdf__(d, di::Discretization, x::PointPattern, l, random_field_mode, H)
    y = cellcounts(Vector{Tuple{Int,Int}}, di, x)
    logpdf__(d, di, y, l, random_field_mode, H)
end
function logpdf__(d, di::Discretization, x::AbstractVector{Tuple{Int, Int}}, l, random_field_mode, H)
    issquare(di) || throw(ArgumentError("Only square discretizations supported."))
    d.range > 0 || return -Inf
    d.sd > 0 || return -Inf
    scale = step(di.gridx)
    fixedfield = d.mean
    if (l != -Inf)
        l1 = LGCPobs(fixedfield .+ random_field_mode, x, scale^2)
        l2 = sum(((cell, count),) -> logabsgamma(count+1)[1], x)
        l = l + l1 - l2
    end
    l, random_field_mode, H
end

function Distributions.mean(m::LogGaussianCoxProcess, d::Discretization)
    size(m.mean) == d.size || error("Discretizations don't agree")
    sd = m.sd
    fixedfield = m.mean
    step(d.gridx) * step(d.gridy) * @.(exp(fixedfield + 0.5*sd^2))
end
function Base.rand(m::LogGaussianCoxProcess, d::Discretization)
    # pts1, pts2 = @. LinRange(d.windowsize / d.gridsize, d.windowsize, d.gridsize)
    # GaussianRandomFields uses a different parametrization for range.
    cov = GaussianRandomFields.CovarianceFunction(2, GaussianRandomFields.Matern(m.range/4, m.nu, σ=m.sd))
    grf = GaussianRandomFields.GaussianRandomField(cov, GaussianRandomFields.CirculantEmbedding(), d.gridx, d.gridy)
    g = GaussianRandomFields.sample(grf)
    intensity = @. exp(m.mean + g)
    rand(PoissonProcess(intensity, d.window), d)
end

function Base.rand(x::PoissonProcess, d::Discretization)#; maxn::Int=typemax(Int))
    size(x.intensity) == d.size || error("Discretizations don't agree")
    x.window == d.window || error("Windows don't agree")
    sx = step(d.gridx)
    sy = step(d.gridy)
    rpois(λ) = rand(Distributions.Poisson(λ * sx*sy))
    counts = rpois.(x.intensity)
    n = sum(counts)
    # if n > maxn
    #     throw(ErrorException)
    # end
    xy = Vector{NTuple{2,Float64}}(undef, n)
    index = 1
    for i in 1:d.size[1], j in 1:d.size[2]
        for _ in 1:counts[i, j]
            xy[index] = i*sx - rand()*sx, j*sy - rand()*sy
            index += 1
        end
    end
    PointPattern(xy, d.window)
end


text:
using Parameters
using ModelParameters

Base.@kwdef mutable struct StreamNode{A} <: NetworkNode{A}
    @network_node

    # https://wiki.ewater.org.au/display/SD41/IHACRES-CMD+-+SRG
    d::A = Param(200.0, bounds=(10.0, 550.0))  # flow threshold
    d2::A = Param(2.0, bounds=(10.0, 500.0))   # flow threshold2
    e::A = Param(1.0, bounds=(0.1, 1.5))  # temperature to PET conversion factor
    f::A = Param(0.8, bounds=(0.01, 3.0))  # plant stress threshold factor (multiplicative factor of d)
    a::A = Param(200.0, bounds=(0.1, 100.0))
    b::A = Param(200.0, bounds=(0.001, 1.0))
    storage_coef::A = Param(2.9, bounds=(0.2, 10.0))
    alpha::A = Param(0.1, bounds=(1e-5, 1 - 1/10^9))

    level_params::Array{A} = [
        Param(-0.01, bounds=(-0.01, -10.0))
        Param(0.8, bounds=(0.0, 1.5))
        Param(4.5, bounds=(0.0, 20.0))
        Param(5.0, bounds=(1.0, 10.0))
        Param(0.35, bounds=(0.0, 1.0))
        Param(1.41, bounds=(-2.0, 2.0))
        Param(-1.45, bounds=(0.0, -2.5))
        Param(6.75, bounds=(0.0, 10.0))
        Param(150.0, bounds=(50.0, 200.0))
    ]

    storage::Array{Float64} = [100.0]
    quickflow::Array{Float64} = [100.0]
    slowflow::Array{Float64} = [100.0]
    outflow::Array{Float64} = []
    effective_rainfall::Array{Float64} = []
    et::Array{Float64} = []
    inflow::Array{Float64} = []
    level::Array{Float64} = []
end

function StreamNode(node_id::String, area::Float64, d::Float64, d2::Float64, e::Float64, f::Float64,
                    a::Float64, b::Float64, s_coef::Float64, alpha::Float64,
                    store::Float64, quick::Float64, slow::Float64)
    return StreamNode(
        node_id=node_id,
        area=area,
        d=d,
        d2=d2,
        e=e,
        f=f,
        a=a,
        b=b,
        storage_coef=s_coef,
        alpha=alpha,
        storage=[store],
        quickflow=[quick],
        slowflow=[slow]
        )
end

function update_state(s_node::StreamNode, storage::Float64, e_rainfall::Float64, et::Float64, qflow_store::Float64, sflow_store::Float64, outflow::Float64, level::Float64)
    push!(s_node.storage, storage)
    push!(s_node.effective_rainfall, e_rainfall)
    push!(s_node.et, et)
    push!(s_node.outflow, outflow)

    push!(s_node.quickflow, qflow_store)
    push!(s_node.slowflow, sflow_store)
    push!(s_node.level, level)
end

"""
Run node to calculate outflow and update state.

Parameters
----------
timestep: int, time step
rain: float, rainfall
evap: float, evapotranspiration
extractions: float, irrigation and other water extractions
gw_exchange: float, flux in ML - positive is contribution to stream, negative is infiltration
loss: float,

Returns
----------
float, outflow from node
"""
function run_node!(s_node::StreamNode,
                  rain::Float64,
                  evap::Float64,
                  inflow::Float64,
                  ext::Float64,
                  gw_exchange::Float64=0.0,
                  loss::Float64=0.0)::Tuple{Float64, Float64}
    arr_len = length(s_node.storage)
    current_store = s_node.storage[arr_len]

    interim_results = [0.0, 0.0, 0.0]
    @ccall ihacres.calc_ft_interim(interim_results::Ptr{Cdouble},
                                   current_store::Cdouble,
                                   rain::Cdouble,
                                   s_node.d::Cdouble,
                                   s_node.d2::Cdouble,
                                   s_node.alpha::Cdouble)::Cvoid

    (mf, e_rainfall, recharge) = interim_results

    et::Float64 = @ccall ihacres.calc_ET(
        s_node.e::Cdouble,
        evap::Cdouble,
        mf::Cdouble,
        s_node.f::Cdouble,
        s_node.d::Cdouble
    )::Cdouble

    cmd::Float64 = @ccall ihacres.calc_cmd(
        current_store::Cdouble,
        rain::Cdouble,
        s_node.d::Cdouble,
        s_node.d2::Cdouble,
        s_node.alpha::Cdouble
    )::Cdouble

    # var inflow = 0.0
    # for nid in s_node.prev_node:
    #     inflow += s_node.prev_node[nid].run(timestep, rain_evap, ext)
    # # End for
    push!(s_node.inflow, inflow)

    flow_results = [0.0, 0.0, 0.0]
    @ccall ihacres.calc_ft_flows(
        flow_results::Ptr{Cdouble},
        s_node.quickflow[arr_len]::Cdouble,
        s_node.slowflow[arr_len]::Cdouble,
        e_rainfall::Cdouble,
        recharge::Cdouble,
        s_node.area::Cdouble,
        s_node.a::Cdouble,
        s_node.b::Cdouble,
        loss::Cdouble
    )::Cvoid

    (quick_store, slow_store, outflow) = flow_results

    # if self.next_node:  # and ('dam' not in self.next_node.node_type):
    #     cmd, outflow = routing(cmd, s_node.storage_coef, inflow, outflow, ext, gamma=gw_exchange)
    # else:
    #     outflow = calc_outflow(outflow, ext)
    # # End if
    routing_res = [0.0, 0.0]
    @ccall ihacres.routing(
        routing_res::Ptr{Cdouble},
        cmd::Cdouble,
        s_node.storage_coef::Cdouble,
        inflow::Cdouble,
        outflow::Cdouble,
        ext::Cdouble,
        gw_exchange::Cdouble)::Cvoid

    (cmd, outflow) = routing_res

    level::Float64 = @ccall ihacres.calc_ft_level(
        outflow::Cdouble,
        s_node.level_params::Ptr{Cdouble}
    )::Cdouble

    update_state(s_node, cmd, e_rainfall, et, quick_store, slow_store, outflow, level)

    return outflow, level
end

"""
"""
function update_params!(node::StreamNode, area::Float64, d::Float64, d2::Float64, e::Float64, f::Float64,
                        a::Float64, b::Float64, s_coef::Float64, alpha::Float64,
                        p1::Float64, p2::Float64, p3::Float64, p4::Float64, p5::Float64, p6::Float64, p7::Float64, p8::Float64, CTF::Float64)
    node.area = area
    node.d = d
    node.d2 = d2
    node.e = e
    node.f = f
    node.a = a
    node.b = b
    node.storage_coef = s_coef
    node.alpha = alpha
    node.level_params = [p1, p2, p3, p4, p5, p6, p7, p8, CTF]

    return nothing
end

function reset!(s_node::StreamNode)
    s_node.storage = [s_node.storage[1]]

    s_node.quickflow = [s_node.quickflow[1]]
    s_node.slowflow = [s_node.slowflow[1]]
    s_node.outflow = []
    s_node.effective_rainfall = []
    s_node.et = []
    s_node.inflow = []
end

max_stars_repo_name: ConnectedSystems/camgraph

max_issues_repo_name: ConnectedSystems/camgraph

max_forks_repo_name: ConnectedSystems/camgraph


text:
using SEM
using Test

@test true


text:
#the statmech Vumps
#it made sense to seperate both vumpses as
# - leading_boundary primarily works on MPSMultiline
# - they search for different eigenvalues
# - ham vumps should use Lanczos, this has to use arnoldi
# - this vumps updates entire collumns (state[:,i]); incompatible with InfiniteMPS
# - (a)c-prime takes a different number of arguments
# - it's very litle duplicate code, but together it'd be a bit more convoluted (primarily because of the indexing way)

"
    leading_boundary(state,opp,alg,pars=params(state,ham))

    approximate the leading eigenvector for opp
"
function leading_boundary(state::InfiniteMPS,H,alg,pars=params(state,H))
    (st,pr,de) = leading_boundary(convert(MPSMultiline,state),H,alg,pars)
    return convert(InfiniteMPS,st),pr,de
end

function leading_boundary(state::MPSMultiline, H,alg::Vumps,pars = params(state,H))
    galerkin  = 1+alg.tol_galerkin
    iter       = 1

    newAs = similar(state.AL)

    while true
        eigalg = Arnoldi(tol=alg.tol_galerkin/10)

        acjobs = map(1:size(state,2)) do col
            @Threads.spawn eigsolve(RecursiveVec(state.AC[:,col]), 1, :LM, eigalg) do x
                tasks = map(1:length(x)) do row
                    @Threads.spawn ac_prime(x[row], row,col, state, pars)
                end

                RecursiveVec(circshift(fetch.(tasks),1))
            end
        end

        cjobs = map(1:size(state,2)) do col
            @Threads.spawn eigsolve(RecursiveVec(state.CR[:,col]), 1, :LM, eigalg) do x
                tasks = map(1:length(x)) do row
                    @Threads.spawn c_prime(x[row], row,col, state, pars)
                end

                RecursiveVec(circshift(fetch.(tasks),1))
            end
        end

        for col in 1:size(state,2)
            (e,vac,ch) = fetch(acjobs[col])
            (e,vc,ch) = fetch(cjobs[col])

            for row in 1:size(state,1)
                QAc,_ = TensorKit.leftorth!(vac[1][row], alg=TensorKit.QRpos())
                Qc,_  = TensorKit.leftorth!(vc[1][row], alg=TensorKit.QRpos())
                newAs[row,col] = QAc*adjoint(Qc)
            end

        end

        state = MPSMultiline(newAs; leftgauged=true,tol = alg.tol_gauge, maxiter = alg.orthmaxiter)
        galerkin = calc_galerkin(state, pars)
        alg.verbose && println("vumps @iteration $(iter) galerkin = $(galerkin)")

        if galerkin <= alg.tol_galerkin || iter>=alg.maxiter
            iter>=alg.maxiter && println("vumps didn't converge $(galerkin)")
            return state, pars, galerkin
        end

        (state,pars) = alg.finalize(iter,state,H,pars);

        iter += 1
    end
end

calc_galerkin(state::MPSMultiline, pars) = maximum([norm(leftnull(state.AC[row+1,col])'*ac_prime(state.AC[row,col], row,col, state, pars)) for (row,col) in Iterators.product(1:size(state,1),1:size(state,2))][:])


text:
module IteratedMaps

export LogisticMap, HenonMap, DuffingMap, TentMap

# package code goes here
include("base.jl")
include("maps.jl")

end # module

max_stars_repo_name: wkearn/IteratedMaps

max_issues_repo_name: wkearn/IteratedMaps

max_forks_repo_name: wkearn/IteratedMaps


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENCE in the project root.
# ------------------------------------------------------------------

module GslibIO

using Meshes
using Tables
using GeoStatsBase

using FileIO
using Printf
using DelimitedFiles

include("extended.jl")
include("legacy.jl")

end


text:
# ------------ Call sites ------------ #

@inline function (ctx::SimulateContext)(::typeof(cond),
                                        addr::T,
                                        c::Function,
                                        c_args::Tuple,
                                        a::Function,
                                        a_args::Tuple,
                                        b::Function,
                                        b_args::Tuple) where {T <: Address, K, J}
    visit!(ctx, addr)
    ret, cl = simulate(c, c_args...)
    if ret
        branch_ret, branch_cl = simulate(a, a_args...)
        br_tr = BranchTrace(cl, branch_cl)
        add_call!(ctx, addr, ConditionalBranchCallSite(br_tr, get_score(branch_cl) + get_score(cl), c, c_args, ret, a, a_args, branch_ret))
    else
        branch_ret, branch_cl = simulate(b, b_args...)
        br_tr = BranchTrace(cl, branch_cl)
        add_call!(ctx, addr, ConditionalBranchCallSite(br_tr, get_score(branch_cl) + get_score(cl), c, c_args, ret, b, b_args, ret))
    end
end


text:
using BinDeps
using Conda

Conda.add("biopython")

@BinDeps.setup

deps = [
	mafft = library_dependency("mafft")
	cmake = library_dependency("cmake")
]

provides(AptGet, Dict{Any,Any}("mafft" => mafft, "cmake" => cmake))