text:
using DataFrames
using StanSample, Random

#Random.seed!(123)

cd(@__DIR__)

LBA_stan = "
functions{

     real lba_pdf(real t, real b, real A, real v, real s){
          //PDF of the LBA model

          real b_A_tv_ts;
          real b_tv_ts;
          real term_1;
          real term_2;
          real term_3;
          real term_4;
          real pdf;

          b_A_tv_ts = (b - A - t*v)/(t*s);
          b_tv_ts = (b - t*v)/(t*s);
          term_1 = v*Phi(b_A_tv_ts);
          term_2 = s*exp(normal_lpdf(b_A_tv_ts|0,1));
          term_3 = v*Phi(b_tv_ts);
          term_4 = s*exp(normal_lpdf(b_tv_ts|0,1));
          pdf = (1/A)*(-term_1 + term_2 + term_3 - term_4);

          return pdf;
     }

     real lba_cdf(real t, real b, real A, real v, real s){
          //CDF of the LBA model

          real b_A_tv;
          real b_tv;
          real ts;
          real term_1;
          real term_2;
          real term_3;
          real term_4;
          real cdf;

          b_A_tv = b - A - t*v;
          b_tv = b - t*v;
          ts = t*s;
          term_1 = b_A_tv/A * Phi(b_A_tv/ts);
          term_2 = b_tv/A   * Phi(b_tv/ts);
          term_3 = ts/A     * exp(normal_lpdf(b_A_tv/ts|0,1));
          term_4 = ts/A     * exp(normal_lpdf(b_tv/ts|0,1));
          cdf = 1 + term_1 - term_2 + term_3 - term_4;

          return cdf;

     }

     real lba_lpdf(matrix RT, real k, real A, vector v, real s, real tau){

          real t;
          real b;
          real cdf;
          real pdf;
          vector[rows(RT)] prob;
          real out;
          real prob_neg;

          b = A + k;
          for (i in 1:rows(RT)){
               t = RT[i,1] - tau;
               if(t > 0){
                    cdf = 1;

                    for(j in 1:num_elements(v)){
                         if(RT[i,2] == j){
                              pdf = lba_pdf(t, b, A, v[j], s);
                         }else{
                              cdf = (1-lba_cdf(t, b, A, v[j], s)) * cdf;
                         }
                    }
                    prob_neg = 1;
                    for(j in 1:num_elements(v)){
                         prob_neg = Phi(-v[j]/s) * prob_neg;
                    }
                    prob[i] = pdf*cdf;
                    prob[i] = prob[i]/(1-prob_neg);
                    if(prob[i] < 1e-10){
                         prob[i] = 1e-10;
                    }

               }else{
                    prob[i] = 1e-10;
               }
          }
          out = sum(log(prob));
          return out;
     }
}

data{
     int N;
     int Nc;
     vector[N] rt;
     vector[N] choice;
}

parameters {
     real<lower=0> k;
     real<lower=0> A;
     real<lower=0> tau;
     vector<lower=0>[Nc] v;
}

model {
     real s;
     matrix[N,2] RT;
     s=1;
     RT[:,1] = rt;
     RT[:,2] = choice;
     k ~ normal(.5,1)T[0,];
     A ~ normal(.5,1)T[0,];
     tau ~ normal(.5,.5)T[0,];
     for(n in 1:Nc){
          v[n] ~ normal(2,1)T[0,];
     }
     RT ~ lba(k,A,v,s,tau);
}
";

LBA_data = (Nc = 3, N = 200, rt = [
  0.5172563185035218, 0.8128454208564759, 0.5287927016034656, 0.6052793579150275,
  0.723151190095471, 0.7460198463976369, 0.45431036513273687, 0.9057048332956246,
  0.5706689553572384, 0.8084971308318065, 0.619517451443632, 0.5983122383151855,
  0.6435242973178197, 1.1321032066600665, 0.6647108084245479, 0.49330145136481923,
  0.6149854662216923,0.6511213339085504,0.5568235912567999,0.6044452252118269,
  0.49566438865268736,0.6279681980301477,0.8117378786097601,0.6172740150854072,
  0.4948989304604639,0.7675369226331317,0.46353371786271824,0.645651483256734,
  0.8963288513778418,0.7146600910671628,0.5408767499701641,0.9650536609042444,
  0.5578421275764064,0.5795453896321787,0.6057684089639942,1.5631874975130597,
  0.629808141744541,1.022421056465441,0.6996055605296965,0.5106721416624329,
  0.5917700094202267,0.7091020907029699,0.6119760415323193,0.7014004524524569,
  0.6011453541927583,0.6164693504868568,0.5663668177526784,0.8931165405119605,
  0.5635466158342433,0.5220944220791173,0.5495697677181935,0.462466932393734,
  0.5415455854737724,0.6369646691761369,0.6809963375645367,0.552897656359747,
  0.6188821565306044,0.5030067555354056,0.6057070496540626,0.6072234852359264,
  0.5030390039084438,0.7605180052353402,0.8753965811982178,0.7100046668672728,
  0.6786320746373478,0.5204538842008222,0.5344689749468778,0.5434745621403249,
  0.5981603048774546,0.7328587072762123,0.644715791453861,0.6243597629188848,
  1.3855200658432651,0.48084093691859586,0.7689218001385749,0.7695036016681377,
  1.0869514748343296,0.7959213856975516,0.5348614356661059,0.4936584374011437,
  0.7096991230294177,0.5243555890435926,0.6320511613857858,0.5769494290780035,
  0.7307312216352477,0.5158226425292626,0.6857721952574352,0.5311504860506873,
  0.6877066575282882,0.5670345547148583,0.7249837316204284,0.7040416637968604,
  0.5454428582470049,0.5977294133097858,0.48001836194920766,0.5152010316305199,
  0.5819798987992713,0.788609952938505,0.5176348868183556,0.6683149041712193,
  0.6896216443937007,0.554724513161415,0.6300558955376948,0.6229862854781085,
  0.54458864060411,0.6312182280659482,0.5528622430791532,0.5677220820816027,
  0.6027445523827388,0.5023157798139705,0.6977735437611517,0.809475166191995,
  0.5704580867030262,0.5046427158534664,0.5255592626852288,0.6143857329367953,
  0.6178402270963708,0.5631464929269593,0.5965224748140201,0.7298831748338575,
  0.6946700553045693,0.4820532443228644,0.6736483636790365,0.6737106523311083,
  0.9232974483628132,0.7354470000550888,0.6726116531295407,0.5883287005892515,
  0.6160263479549937,0.5573678777576276,0.7625656142155866,0.5292714975831675,
  0.7246075325766295,0.7708827544279314,0.5336499540403102,0.8051382314453717,
  0.6442494256307568,0.9334582322639116,0.9461083952852346,0.6196330941672238,
  0.5586124687600272,0.5566129176765593,0.7904878901502628,0.7479393473413253,
  0.6574586477816032,0.4944981578764378,0.6510666622321069,0.6166906369729925,
  0.5000436387722222,0.4752189615622799,0.5959538562155668,0.6938059053886498,
  0.6438274327143417,0.5589433865365045,0.5846831738823322,0.6314707748237262,
  0.5140271877113509,0.49346493970846383,0.6515051208839449,0.8294415751644921,
  2.0321926771647276,0.49412525362775855,0.5303283741433148,0.8523027635055367,
  0.5414768731824661,0.4695826111771622,0.48833888190741115,0.7167178014989934,
  1.1408781345289776,0.7140074385674222,0.6272348029864163,0.48052613023089863,
  0.6281778536199709,0.6012504548489284,0.5462307192768121,0.6239689876094862,
  0.8432350862861097,0.7694134803246382,0.5391600922718652,0.63566692258107,
  0.871898665375137,0.5966805648137725,0.5756731591390922,0.9990694566686539,
  0.5242361805544791,0.4813390113561128,0.5665089258206126,0.5094476254864568,
  0.553923456556689,0.6379083181936438,0.5265558116761913,0.6318610211301654,
  0.4840139865984098,0.567837742224435,1.1478087390790883,0.6157830093722181,
  0.4597544454463623,0.7692541965561321,0.5552938695617227,0.5893517515699693],
  choice = [1,2,3,3,3,2,1,3,3,3,2,2,1,3,3,3,1,3,1,3,3,2,3,3,3,3,3,2,1,2,3,2,2,3,2,3,2,3,3,3,2,3,2,2,1,
  3,2,3,3,3,1,2,2,3,3,2,2,3,3,2,1,2,1,3,1,3,3,2,2,2,2,2,2,3,1,3,1,1,2,2,3,3,3,3,3,3,2,3,3,3,1,2,3,2,2,
  3,2,1,3,3,3,3,3,1,3,3,3,3,2,3,2,2,3,2,3,2,2,3,3,2,3,3,3,2,2,3,3,3,2,2,1,3,1,3,1,3,3,1,1,3,1,3,1,1,2,
  2,1,2,3,2,3,2,2,1,1,3,2,3,2,3,3,1,1,3,2,2,3,3,2,2,1,2,3,1,1,3,2,3,2,3,3,3,3,2,2,3,3,3,3,3,3,2,2,2,2,2,3,3,2,3]
);

# This run tests passing a data file name as data in the stan_sample() call

sm = SampleModel("LBA", LBA_stan;
  method = StanSample.Sample(adapt = StanSample.Adapt(delta = 0.95)));

@time rc = stan_sample(sm; data="LBA.R")

if success(rc)
  samples = read_samples(sm, :dataframe)

  # Show the same output in DataFrame format
  sdf = StanSample.read_summary(sm)
  display(sdf)
end


text:
using .LinearSpace

struct LongDelays{T} end
struct Threading{T} end

import Base: transpose, +, -, *

abstract type AbstractOperator end

#@linearspace AbstractOperator{T} T

"""
*Atomic operator*: one that assemblechunk can deal with
"""
abstract type Operator <: AbstractOperator end

mutable struct TransposedOperator <: Operator
    op::Operator
end

scalartype(op::TransposedOperator) = scalartype(op.op)
defaultquadstrat(op::TransposedOperator, tfs, bfs) = defaultquadstrat(op.op, tfs, bfs)

mutable struct LinearCombinationOfOperators{T} <: AbstractOperator
    coeffs::Vector{T}
    ops::Vector
end

function scalartype(op::LinearCombinationOfOperators{T}) where T
    W = promote_type(T, scalartype(op.ops[1]))
    for i in 2:length(op.ops)
        W = promote_type(W, scalartype(op.ops[i]))
    end
    W
end

function derive(a::LinearCombinationOfOperators)
    coeffs = copy(a.coeffs)
    ops = [derive(op) for op in a.ops]
    return LinearCombinationOfOperators(coeffs, ops)
end

function +(a::LinearCombinationOfOperators, b::AbstractOperator)
    LinearCombinationOfOperators(
        [a.coeffs;[1.0]],
        [a.ops;[b]]
    )
end

function +(a::LinearCombinationOfOperators, b::LinearCombinationOfOperators)
    LinearCombinationOfOperators(
        [a.coeffs; b.coeffs],
        [a.ops; b.ops]
    )
end

+(a::AbstractOperator, b::LinearCombinationOfOperators) = b + a
+(a::AbstractOperator, b::Number) = a + (b * Identity())
+(a::Number, b::AbstractOperator) = b + a

*(a::Number, b::AbstractOperator) = LinearCombinationOfOperators([a], [b])
*(a::Number, b::LinearCombinationOfOperators) = LinearCombinationOfOperators(a * b.coeffs, b.ops)
-(a::AbstractOperator, b::AbstractOperator) = a + (-1.0) * b
-(a::AbstractOperator) = (-1.0) * a

function +(a::AbstractOperator, b::AbstractOperator)
    LinearCombinationOfOperators(
        [1.0, 1.0],
        [a, b]
    )
end

transpose(op::TransposedOperator) = op.op
transpose(op::Operator) = TransposedOperator(op)

defaultquadstrat(lc::LinearCombinationOfOperators, tfs, bfs) =
    [defaultquadstrat(op,tfs,bfs) for op in lc.ops]

function assemble(operator::AbstractOperator, test_functions, trial_functions;
    storage_policy = Val{:bandedstorage},
    long_delays_policy = LongDelays{:compress},
    threading = Threading{:multi},
    quadstrat=defaultquadstrat(operator, test_functions, trial_functions))

    Z, store = allocatestorage(operator, test_functions, trial_functions,
        storage_policy, long_delays_policy)
    assemble!(operator, test_functions, trial_functions, store, threading; quadstrat)
    return Z()
end

function assemblerow(operator::AbstractOperator, test_functions, trial_functions,
    storage_policy = Val{:bandedstorage},
    long_delays_policy = LongDelays{:ignore};
    quadstrat=defaultquadstrat(operator, test_functions, trial_functions))

    Z, store = allocatestorage(operator, test_functions, trial_functions,
        storage_policy, long_delays_policy)
    assemblerow!(operator, test_functions, trial_functions, store; quadstrat)

    Z()
end

function assemblecol(operator::AbstractOperator, test_functions, trial_functions,
    storage_policy = Val{:bandestorage},
    long_delays_policy = LongDelays{:ignore};
    quadstrat=defaultquadstrat(operator, test_functions, trial_functions))

    Z, store = allocatestorage(operator, test_functions, trial_functions,
        storage_policy, long_delays_policy)
    assemblecol!(operator, test_functions, trial_functions, store; quadstrat)

    Z()
end

function allocatestorage(operator::AbstractOperator, test_functions, trial_functions,
    storage_trait, longdelays_trait)

    T = promote_type(
        scalartype(operator)       ,
        scalartype(test_functions) ,
        scalartype(trial_functions),
    )
    Z = Matrix{T}(undef,
        numfunctions(test_functions),
        numfunctions(trial_functions),
    )
    fill!(Z, 0)
    store(v,m,n) = (Z[m,n] += v)
    return ()->Z, store
end

function allocatestorage(operator::LinearCombinationOfOperators,
        test_functions::SpaceTimeBasis, trial_functions::SpaceTimeBasis,
        storage_policy::Type{Val{:bandedstorage}},
        long_delays_policy::Type{LongDelays{:ignore}})

    # TODO: remove this ugly, ugly patch
    return allocatestorage(operator.ops[end], test_functions, trial_functions,
        storage_policy, long_delays_policy)
end

function allocatestorage(operator::LinearCombinationOfOperators,
        test_functions::SpaceTimeBasis, trial_functions::SpaceTimeBasis,
        storage_policy::Type{Val{S}},
        long_delays_policy::Type{LongDelays{L}}) where {L,S}

    # TODO: remove this ugly, ugly patch
    return allocatestorage(operator.ops[end], test_functions, trial_functions,
        storage_policy, long_delays_policy)
end

function assemble!(operator::Operator, test_functions::Space, trial_functions::Space, store,
    threading::Type{Threading{:multi}} = Threading{:multi};
    quadstrat=defaultquadstrat(operator, test_functions, trial_functions))

    @info "Multi-threaded assembly:"

    P = Threads.nthreads()
    numchunks = P
    @assert numchunks >= 1
    splits = [round(Int,s) for s in range(0, stop=numfunctions(test_functions), length=numchunks+1)]

    Threads.@threads for i in 1:P
        lo, hi = splits[i]+1, splits[i+1]
        lo <= hi || continue
        test_functions_p = subset(test_functions, lo:hi)
        store1 = (v,m,n) -> store(v,lo+m-1,n)
        assemblechunk!(operator, test_functions_p, trial_functions, store1; quadstrat)
    end

end

function assemble!(operator::Operator, test_functions::Space, trial_functions::Space, store,
    threading::Type{Threading{:single}};
    quadstrat=defaultquadstrat(operator, test_functions, trial_functions))

    @info "Single-threaded assembly"

    assemblechunk!(operator, test_functions, trial_functions, store; quadstrat)
end

function assemble!(op::TransposedOperator, tfs::Space, bfs::Space, store;
    quadstrat=defaultquadstrat(op, tfs, bfs))

    store1(v,m,n) = store(v,n,m)
    assemble!(op.op, bfs, tfs, store1; quadstrat)
end

function assemble!(op::LinearCombinationOfOperators, tfs::AbstractSpace, bfs::AbstractSpace,
    store, threading = Threading{:multi};
    quadstrat=defaultquadstrat(op, tfs, bfs))

    for (a,A,qs) in zip(op.coeffs, op.ops, quadstrat)
        store1(v,m,n) = store(a*v,m,n)
        assemble!(A, tfs, bfs, store1; quadstrat=qs)
    end
end

# Support for direct product spaces
function assemble!(op::Operator, tfs::DirectProductSpace, bfs::Space, store, threading = Threading{:multi};
    quadstrat=defaultquadstrat(op, tfs[1], bfs))
    I = Int[0]
    for s in tfs.factors push!(I, last(I) + numfunctions(s)) end
    for (i,s) in enumerate(tfs.factors)
        store1(v,m,n) = store(v,m + I[i], n)
        assemble!(op, s, bfs, store1; quadstrat)
    end
end

function assemble!(op::Operator, tfs::Space, bfs::DirectProductSpace, store, threading=Threading{:multi};
    quadstrat=defaultquadstrat(op, tfs, bfs[1]))
    J = Int[0]
    for s in bfs.factors push!(J, last(J) + numfunctions(s)) end
    for (j,s) in enumerate(bfs.factors)
        store1(v,m,n) = store(v,m,n + J[j])
        assemble!(op, tfs, s, store1; quadstrat)
    end
end

function assemble!(op::Operator, tfs::DirectProductSpace, bfs::DirectProductSpace, store, threading=Threading{:multi};
    quadstrat=defaultquadstrat(op, tfs[1], bfs[1]))
    I = Int[0]
    for s in tfs.factors push!(I, last(I) + numfunctions(s)) end
    for (i,s) in enumerate(tfs.factors)
        store1(v,m,n) = store(v,m + I[i],n)
        assemble!(op, s, bfs, store1; quadstrat)
    end
end

ext: jl


text:
__precompile__()

module PDFIO

include("Common.jl") #Module Common
include("Cos.jl")    #Module Cos
include("PD.jl")     #Module PD

using .Common
export CDTextString, CDDate, CDRect, getUTCTime

using .PD
export  PDDoc,
            pdDocOpen,
            pdDocClose,
            pdDocGetCatalog,
            pdDocGetNamesDict,
            pdDocGetInfo,
            pdDocGetCosDoc,
            pdDocGetPage,
            pdDocGetPageCount,
            pdDocGetPageRange,
            pdDocHasPageLabels,
            pdDocGetPageLabel,
            pdDocGetOutline,
            pdDocHasSignature,
            pdDocValidateSignatures,
        PDPage,
            pdPageGetContents,
            pdPageIsEmpty,
            pdPageGetCosObject,
            pdPageGetContentObjects,
            pdPageGetMediaBox,
            pdPageGetFonts,
            pdPageGetCropBox,
            pdPageExtractText,
            pdPageGetPageNumber,
        PDPageObject,
            PDPageObjectGroup,
                PDPage_BeginGroup, PDPage_EndGroup,
            PDPageTextObject,
            PDPageMarkedContent,
            PDPageElement,
            PDPageTextRun,
            PDPageInlineImage,
        PDFont,
            pdFontIsBold,
            pdFontIsItalic,
            pdFontIsFixedW,
            pdFontIsAllCap,
            pdFontIsSmallCap,
        PDOutline,
            PDDestination,
            PDOutlineItem,
                pdOutlineItemGetAttr

using .Cos
export  CosDoc,
            cosDocOpen,
            cosDocClose,
            cosDocGetRoot,
            cosDocGetObject,
            cosDocIsEncrypted,
        CosObject,
            CosNull, CosDict, CosString, CosArray, CosStream,
            CosIndirectObjectRef, set!,
        CosBoolean,
            CosTrue, CosFalse,
        CosNumeric,
            CosFloat, CosInt,
        CosName,
            @cn_str,
        CosTreeNode,
            createTreeNode

end # module


text:
"""
Utility functions for Probabilistic Programming - 3
Wouter Kouw
Last update: 07-01-2020
"""

import MCMCChains: Chains

function one_hot(index::Number, length_vector::Integer; nonzero::Bool=false)
    "Map integer to one-hot encoding"

    # Initialize encoded vector
    encoded = zeros(length_vector,)

    # Set i-th index to 1
    encoded[index] = 1.0

    # Check for strictly positive vector
    if nonzero

        # Add epsilon
        encoded .+= 1e-8

        # Re-normalize to 1
        encoded ./= sum(encoded)
    end

    return encoded
end

function mean_chain(chain::Chains, param::Symbol)
    "Quick function to extract sample average from MCMCChains object"
    return vec(convert(Array{Float64,2}, mean(chain[param].value.data, dims=[1, 3])[:,:,1]))
end

function median_chain(chain::Chains, param::Symbol)
    "Quick function to extract median from MCMCChains object"
    return vec(convert(Array{Float64,2}, median(chain[param].value.data, dims=[1, 3])[:,:,1]))
end;

max_stars_repo_name: lisakoolen/BMLIP

max_issues_repo_name: lisakoolen/BMLIP

max_forks_repo_name: lisakoolen/BMLIP


text:
using ACE, BenchmarkTools

pool = ACE.ArrayCache{Float64}()
pool32 = ACE.ArrayCache{Float32}()
genpool = ACE.GenArrayCache()

##

function dosomething!(_A, x)
   A = parent(_A)
   @inbounds begin
      A[1] = 1; A[2] = x
      for n = 3:length(A)
         A[n] = 2 * x * A[n-1] - A[n-2]
      end
   end
   return _A
end

function acquire_N(pool, N, len)
   x = rand()
   for n = 1:N
      x -= 1e-12
      A = ACE.acquire!(pool, Float64, len)
      dosomething!(A, x)
      ACE.release!(A)
   end
   return nothing
end

##

N = 1000
len = 100
acquire_N(pool, N, len)
acquire_N(genpool, N, len)
acquire_N(pool32, N, len)

@btime acquire_N($pool, $N, $len)
@btime acquire_N($genpool, $N, $len)
@btime acquire_N($pool32, $N, $len)

##

A = ACE.acquire!(pool, Float64, len)
gA = ACE.acquire!(genpool, Float64, len)
aA = ACE.acquire!(pool32, Float64, len)
x = rand()

@btime dosomething!($A, $x)
@btime dosomething!($gA, $x)
@btime dosomething!($aA, $x)

##

# @profview acquire_N(genpool, 100_000_000, len)

##
nothing
##

# ##

# function runn(x, N, len)
#    @assert length(x) >= len
#    fill!(x, 0)
#    for n = 1:N, i = 1:len
#       x[i] += i/n
#    end
#    return nothing
# end

# x = rand(1000)
# @btime runn($x, 1000, 1000)

# ##

# function runn2(cache, N, len)
#    x = ACE.acquire!(cache, len)
#    fill!(x.A, 0)
#    for n = 1:N, i = 1:len
#       x[i] += i/n
#    end
#    ACE.release!(x)
#    return nothing
# end

# @btime runn2($pool, 1000, 1000)

# ##

# N = 1000

# for len in [1000, 300, 100, 30, 10]
#    N = 1 # (1000 ÷ len) * 2
#    println("len = $len")
#    x = zeros(len)
#    print(" preallocate: ")
#    @btime runn($x, $N, $len)
#    print("      cached: ")
#    @btime runn2($pool, $N, $len)
#    print("garbage coll: ")
#    @btime runn(zeros($len), $N, $len)
# end


text:
#=
Copyright 2017 INSIGNEO Institute for in silico Medicine

Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at

   http://www.apache.org/licenses/LICENSE-2.0

Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.
=#

# The anastomosis is solved as a reversed bifurcation. Refer to [bifurcations](bifurcations.html)
# documentation for further details.

#
function solveAnastomosis(v1 :: Vessel, v2 :: Vessel, v3 :: Vessel)

  #Unknowns vector
  U = [v1.u[end],
       v2.u[end],
       v3.u[ 1 ],
       sqrt(sqrt(v1.A[end])),
       sqrt(sqrt(v2.A[end])),
       sqrt(sqrt(v3.A[ 1 ]))]

  #Parameters vector
  k1 = sqrt(0.5*3*v1.gamma[end])
  k2 = sqrt(0.5*3*v2.gamma[end])
  k3 = sqrt(0.5*3*v3.gamma[ 1 ])
  k = [k1, k2, k3]

  W = calculateWstarAn(U, k)
  J = calculateJacobianAn(v1, v2, v3, U, k)
  F = calculateFofUAn(v1, v2, v3, U, k, W)

  #Newton-Raphson
  nr_toll_U = 1.e-5
  nr_toll_F = 1.e-5

  while true
    dU = J\(-F)
    # U_new = U + 0.01*dU
    U_new = U + dU

    if any(isnan(dot(F,F)))
      println(F)
      @printf "error at anastomosis with vessels %s, %s, and %s \n" v1.label v2.label v3.label
      break
    end

    u_ok = 0
    f_ok = 0
    for i in 1:length(dU)
      if abs(dU[i]) <= nr_toll_U || abs(F[i]) <= nr_toll_F
        u_ok += 1
        f_ok += 1
      end
    end

    if u_ok == length(dU) || f_ok == length(dU)
      U = U_new
      break
    else
      U = U_new
      W = calculateWstarBif(U, k)
      F = calculateFofUBif(v1, v2, v3, U, k, W)
    end
  end

  #Update vessel quantities
  v1.u[end] = U[1]
  v2.u[end] = U[2]
  v3.u[ 1 ] = U[3]

  v1.A[end] = U[4]*U[4]*U[4]*U[4]
  v2.A[end] = U[5]*U[5]*U[5]*U[5]
  v3.A[ 1 ] = U[6]*U[6]*U[6]*U[6]

  v1.Q[end] = v1.u[end]*v1.A[end]
  v2.Q[end] = v2.u[end]*v2.A[end]
  v3.Q[ 1 ] = v3.u[ 1 ]*v3.A[ 1 ]

  v1.P[end] = pressure(v1.A[end], v1.A0[end], v1.beta[end], v1.Pext)
  v2.P[end] = pressure(v2.A[end], v2.A0[end], v2.beta[end], v2.Pext)
  v3.P[ 1 ] = pressure(v3.A[ 1 ], v3.A0[ 1 ], v3.beta[ 1 ], v3.Pext)

  v1.c[end] = waveSpeed(v1.A[end], v1.gamma[end])
  v2.c[end] = waveSpeed(v2.A[end], v2.gamma[end])
  v3.c[ 1 ] = waveSpeed(v3.A[ 1 ], v3.gamma[ 1 ])

end

#
function calculateWstarAn(U :: Array, k :: Array)

  W1 = U[1] + 4*k[1]*U[4]
  W2 = U[2] + 4*k[2]*U[5]
  W3 = U[3] - 4*k[3]*U[6]

  return [W1, W2, W3]
end

#
function calculateFofUAn(v1 :: Vessel, v2 :: Vessel, v3 :: Vessel,
                           U :: Array,   k :: Array,   W :: Array)

  f1 = U[1] + 4*k[1]*U[4] - W[1]

  f2 = U[2] + 4*k[2]*U[5] - W[2]

  f3 = U[3] - 4*k[3]*U[6] - W[3]

  f4 = U[1]*(U[4]*U[4]*U[4]*U[4]) + U[2]*(U[5]*U[5]*U[5]*U[5]) - U[3]*(U[6]*U[6]*U[6]*U[6])

  f5 = v1.beta[end]*( U[4]^2/sqrt(v1.A0[end])  - 1  ) -
     ( v3.beta[ 1 ]*((U[6])^2/sqrt(v3.A0[1]) - 1) )

  f6 = v2.beta[end]*( U[5]^2/sqrt(v2.A0[end])  - 1  ) -
     ( v3.beta[ 1 ]*((U[6]^2)/sqrt(v3.A0[1]) - 1) )

  return [f1, f2, f3, f4, f5, f6]
end

#
function calculateJacobianAn(v1 :: Vessel, v2 :: Vessel, v3 :: Vessel,
                               U :: Array,   k :: Array)

  J = eye(6)

  J[1,4] =  4*k[1]
  J[2,5] =  4*k[2]
  J[3,6] = -4*k[3]

  J[4,1] =  U[4]*U[4]*U[4]*U[4]
  J[4,2] =  U[5]*U[5]*U[5]*U[5]
  J[4,3] = -U[6]*U[6]*U[6]*U[6]
  J[4,4] =  4*U[1]*(U[4]*U[4]*U[4])
  J[4,5] =  4*U[2]*(U[5]*U[5]*U[5])
  J[4,6] = -4*U[3]*(U[6]*U[6]*U[6])

  J[5,4] =  2*v1.beta[end]*U[4]/sqrt(v1.A0[end])
  J[5,5] =  0.
  J[5,6] = -2*v3.beta[ 1 ]*U[6]/sqrt(v3.A0[ 1 ])

  J[6,5] =  2*v2.beta[end]*U[5]/sqrt(v2.A0[end])
  J[6,6] = -2*v3.beta[ 1 ]*U[6]/sqrt(v3.A0[ 1 ])

  return J
end


text:
# build the mock dynamic library
run(`make`)

isdefined(:DL_LOAD_PATH) || (DL_LOAD_PATH = ASCIIString[])
# put the mock MultiNest library on DL_LOAD_PATH
push!(DL_LOAD_PATH, dirname(@__FILE__()))

using MultiNest
using Base.Test

# write your own tests here
@test 1 == 1

# done with tests, clean up
run(`make clean`)


text:
function regulateperiphery!(system::CVSystem,n::Int64,terms::Vector{Int64})
    for i = 1:length(terms)
        if system.branches.group[terms[i]] == "lower"
            system.branches.term[terms[i]].R[2] = system.cns.R2L[n+2];
            system.branches.term[terms[i]].R[3] = system.cns.R3L[n+2];
            system.branches.term[terms[i]].C[4] = system.cns.C4L[n+2];
            system.branches.term[terms[i]].C[5] = system.cns.C5L[n+2];
            system.branches.term[terms[i]].V0[4] = system.cns.V4L[n+2];
            system.branches.term[terms[i]].V0[5] = system.cns.V5L[n+2];
        else
            system.branches.term[terms[i]].R[2] = system.cns.R2U[n+2];
            system.branches.term[terms[i]].R[3] = system.cns.R3U[n+2];
            system.branches.term[terms[i]].C[4] = system.cns.C4U[n+2];
            system.branches.term[terms[i]].C[5] = system.cns.C5U[n+2];
            system.branches.term[terms[i]].V0[4] = system.cns.V4U[n+2];
            system.branches.term[terms[i]].V0[5] = system.cns.V5U[n+2];
        end
    end
end

max_stars_repo_name: dcanuto/CardioModeling

max_issues_repo_name: dcanuto/CardioModeling

max_forks_repo_name: dcanuto/CardioModeling


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:5] <= 1)

                     @variable(m, -36.7700820821577 <= q <= 37.6713142442229)

                     add_NL_constraint(m, :(log(1 + exp(0.9266287961121966 + 0.9630341445616355*log(1 + exp(0.9080159995418975 + 0.10769169268185319*log(1 + exp(-0.02399021036854254 + -0.4943608904262171*$(x[1]) + -0.3028702711871105*$(x[2]) + 0.7730701583831108*$(x[3]) + -0.5656579689281509*$(x[4]) + 0.08835909373115314*$(x[5]))) + 0.2597058016512084*log(1 + exp(-0.6908874409241781 + -0.9493883109743697*$(x[1]) + 0.1034335229123764*$(x[2]) + 0.3901179347519834*$(x[3]) + -0.16073609758033847*$(x[4]) + 0.052390589356133166*$(x[5]))) + 0.8597063403720888*log(1 + exp(0.7510327333371474 + 0.48933101371416265*$(x[1]) + 0.5299590321169751*$(x[2]) + 0.05546767568522437*$(x[3]) + -0.4663274814822995*$(x[4]) + -0.6567285643500202*$(x[5]))) + 0.17950756837409587*log(1 + exp(-0.27157848416188113 + 0.560003468778171*$(x[1]) + 0.2933473130688298*$(x[2]) + 0.8808820816540863*$(x[3]) + 0.6726368959629845*$(x[4]) + -0.19953521813356856*$(x[5]))))) + 0.8197021439185934*log(1 + exp(-0.9264374759333807 + -0.7028105280156569*log(1 + exp(-0.02399021036854254 + -0.4943608904262171*$(x[1]) + -0.3028702711871105*$(x[2]) + 0.7730701583831108*$(x[3]) + -0.5656579689281509*$(x[4]) + 0.08835909373115314*$(x[5]))) + 0.6656745963467641*log(1 + exp(-0.6908874409241781 + -0.9493883109743697*$(x[1]) + 0.1034335229123764*$(x[2]) + 0.3901179347519834*$(x[3]) + -0.16073609758033847*$(x[4]) + 0.052390589356133166*$(x[5]))) + -0.9099166873549671*log(1 + exp(0.7510327333371474 + 0.48933101371416265*$(x[1]) + 0.5299590321169751*$(x[2]) + 0.05546767568522437*$(x[3]) + -0.4663274814822995*$(x[4]) + -0.6567285643500202*$(x[5]))) + 0.9082646170966706*log(1 + exp(-0.27157848416188113 + 0.560003468778171*$(x[1]) + 0.2933473130688298*$(x[2]) + 0.8808820816540863*$(x[3]) + 0.6726368959629845*$(x[4]) + -0.19953521813356856*$(x[5]))))) + -0.8564198456770362*log(1 + exp(-0.7118497063916944 + 0.36586202274128254*log(1 + exp(-0.02399021036854254 + -0.4943608904262171*$(x[1]) + -0.3028702711871105*$(x[2]) + 0.7730701583831108*$(x[3]) + -0.5656579689281509*$(x[4]) + 0.08835909373115314*$(x[5]))) + -0.9129427463398758*log(1 + exp(-0.6908874409241781 + -0.9493883109743697*$(x[1]) + 0.1034335229123764*$(x[2]) + 0.3901179347519834*$(x[3]) + -0.16073609758033847*$(x[4]) + 0.052390589356133166*$(x[5]))) + -0.5191517354247042*log(1 + exp(0.7510327333371474 + 0.48933101371416265*$(x[1]) + 0.5299590321169751*$(x[2]) + 0.05546767568522437*$(x[3]) + -0.4663274814822995*$(x[4]) + -0.6567285643500202*$(x[5]))) + 0.42477634545139953*log(1 + exp(-0.27157848416188113 + 0.560003468778171*$(x[1]) + 0.2933473130688298*$(x[2]) + 0.8808820816540863*$(x[3]) + 0.6726368959629845*$(x[4]) + -0.19953521813356856*$(x[5]))))) + 0.7868467121279954*log(1 + exp(0.8429853377330745 + -0.6377098680254001*log(1 + exp(-0.02399021036854254 + -0.4943608904262171*$(x[1]) + -0.3028702711871105*$(x[2]) + 0.7730701583831108*$(x[3]) + -0.5656579689281509*$(x[4]) + 0.08835909373115314*$(x[5]))) + 0.04468698484435718*log(1 + exp(-0.6908874409241781 + -0.9493883109743697*$(x[1]) + 0.1034335229123764*$(x[2]) + 0.3901179347519834*$(x[3]) + -0.16073609758033847*$(x[4]) + 0.052390589356133166*$(x[5]))) + 0.3607403638969684*log(1 + exp(0.7510327333371474 + 0.48933101371416265*$(x[1]) + 0.5299590321169751*$(x[2]) + 0.05546767568522437*$(x[3]) + -0.4663274814822995*$(x[4]) + -0.6567285643500202*$(x[5]))) + 0.4004266402265344*log(1 + exp(-0.27157848416188113 + 0.560003468778171*$(x[1]) + 0.2933473130688298*$(x[2]) + 0.8808820816540863*$(x[3]) + 0.6726368959629845*$(x[4]) + -0.19953521813356856*$(x[5]))))))) + log(1 + exp(-0.9209290943615045 + -0.04560727661330999*log(1 + exp(0.9080159995418975 + 0.10769169268185319*log(1 + exp(-0.02399021036854254 + -0.4943608904262171*$(x[1]) + -0.3028702711871105*$(x[2]) + 0.7730701583831108*$(x[3]) + -0.5656579689281509*$(x[4]) + 0.08835909373115314*$(x[5]))) + 0.2597058016512084*log(1 + exp(-0.6908874409241781 + -0.9493883109743697*$(x[1]) + 0.1034335229123764*$(x[2]) + 0.3901179347519834*$(x[3]) + -0.16073609758033847*$(x[4]) + 0.052390589356133166*$(x[5]))) + 0.8597063403720888*log(1 + exp(0.7510327333371474 + 0.48933101371416265*$(x[1]) + 0.5299590321169751*$(x[2]) + 0.05546767568522437*$(x[3]) + -0.4663274814822995*$(x[4]) + -0.6567285643500202*$(x[5]))) + 0.17950756837409587*log(1 + exp(-0.27157848416188113 + 0.560003468778171*$(x[1]) + 0.2933473130688298*$(x[2]) + 0.8808820816540863*$(x[3]) + 0.6726368959629845*$(x[4]) + -0.19953521813356856*$(x[5]))))) + -0.04187997574232316*log(1 + exp(-0.9264374759333807 + -0.7028105280156569*log(1 + exp(-0.02399021036854254 + -0.4943608904262171*$(x[1]) + -0.3028702711871105*$(x[2]) + 0.7730701583831108*$(x[3]) + -0.5656579689281509*$(x[4]) + 0.08835909373115314*$(x[5]))) + 0.6656745963467641*log(1 + exp(-0.6908874409241781 + -0.9493883109743697*$(x[1]) + 0.1034335229123764*$(x[2]) + 0.3901179347519834*$(x[3]) + -0.16073609758033847*$(x[4]) + 0.052390589356133166*$(x[5]))) + -0.9099166873549671*log(1 + exp(0.7510327333371474 + 0.48933101371416265*$(x[1]) + 0.5299590321169751*$(x[2]) + 0.05546767568522437*$(x[3]) + -0.4663274814822995*$(x[4]) + -0.6567285643500202*$(x[5]))) + 0.9082646170966706*log(1 + exp(-0.27157848416188113 + 0.560003468778171*$(x[1]) + 0.2933473130688298*$(x[2]) + 0.8808820816540863*$(x[3]) + 0.6726368959629845*$(x[4]) + -0.19953521813356856*$(x[5]))))) + 0.009897273687461627*log(1 + exp(-0.7118497063916944 + 0.36586202274128254*log(1 + exp(-0.02399021036854254 + -0.4943608904262171*$(x[1]) + -0.3028702711871105*$(x[2]) + 0.7730701583831108*$(x[3]) + -0.5656579689281509*$(x[4]) + 0.08835909373115314*$(x[5]))) + -0.9129427463398758*log(1 + exp(-0.6908874409241781 + -0.9493883109743697*$(x[1]) + 0.1034335229123764*$(x[2]) + 0.3901179347519834*$(x[3]) + -0.16073609758033847*$(x[4]) + 0.052390589356133166*$(x[5]))) + -0.5191517354247042*log(1 + exp(0.7510327333371474 + 0.48933101371416265*$(x[1]) + 0.5299590321169751*$(x[2]) + 0.05546767568522437*$(x[3]) + -0.4663274814822995*$(x[4]) + -0.6567285643500202*$(x[5]))) + 0.42477634545139953*log(1 + exp(-0.27157848416188113 + 0.560003468778171*$(x[1]) + 0.2933473130688298*$(x[2]) + 0.8808820816540863*$(x[3]) + 0.6726368959629845*$(x[4]) + -0.19953521813356856*$(x[5]))))) + -0.4138607422831795*log(1 + exp(0.8429853377330745 + -0.6377098680254001*log(1 + exp(-0.02399021036854254 + -0.4943608904262171*$(x[1]) + -0.3028702711871105*$(x[2]) + 0.7730701583831108*$(x[3]) + -0.5656579689281509*$(x[4]) + 0.08835909373115314*$(x[5]))) + 0.04468698484435718*log(1 + exp(-0.6908874409241781 + -0.9493883109743697*$(x[1]) + 0.1034335229123764*$(x[2]) + 0.3901179347519834*$(x[3]) + -0.16073609758033847*$(x[4]) + 0.052390589356133166*$(x[5]))) + 0.3607403638969684*log(1 + exp(0.7510327333371474 + 0.48933101371416265*$(x[1]) + 0.5299590321169751*$(x[2]) + 0.05546767568522437*$(x[3]) + -0.4663274814822995*$(x[4]) + -0.6567285643500202*$(x[5]))) + 0.4004266402265344*log(1 + exp(-0.27157848416188113 + 0.560003468778171*$(x[1]) + 0.2933473130688298*$(x[2]) + 0.8808820816540863*$(x[3]) + 0.6726368959629845*$(x[4]) + -0.19953521813356856*$(x[5]))))))) + log(1 + exp(0.41055131347137674 + 0.5074173875637427*log(1 + exp(0.9080159995418975 + 0.10769169268185319*log(1 + exp(-0.02399021036854254 + -0.4943608904262171*$(x[1]) + -0.3028702711871105*$(x[2]) + 0.7730701583831108*$(x[3]) + -0.5656579689281509*$(x[4]) + 0.08835909373115314*$(x[5]))) + 0.2597058016512084*log(1 + exp(-0.6908874409241781 + -0.9493883109743697*$(x[1]) + 0.1034335229123764*$(x[2]) + 0.3901179347519834*$(x[3]) + -0.16073609758033847*$(x[4]) + 0.052390589356133166*$(x[5]))) + 0.8597063403720888*log(1 + exp(0.7510327333371474 + 0.48933101371416265*$(x[1]) + 0.5299590321169751*$(x[2]) + 0.05546767568522437*$(x[3]) + -0.4663274814822995*$(x[4]) + -0.6567285643500202*$(x[5]))) + 0.17950756837409587*log(1 + exp(-0.27157848416188113 + 0.560003468778171*$(x[1]) + 0.2933473130688298*$(x[2]) + 0.8808820816540863*$(x[3]) + 0.6726368959629845*$(x[4]) + -0.19953521813356856*$(x[5]))))) + 0.7865169798440914*log(1 + exp(-0.9264374759333807 + -0.7028105280156569*log(1 + exp(-0.02399021036854254 + -0.4943608904262171*$(x[1]) + -0.3028702711871105*$(x[2]) + 0.7730701583831108*$(x[3]) + -0.5656579689281509*$(x[4]) + 0.08835909373115314*$(x[5]))) + 0.6656745963467641*log(1 + exp(-0.6908874409241781 + -0.9493883109743697*$(x[1]) + 0.1034335229123764*$(x[2]) + 0.3901179347519834*$(x[3]) + -0.16073609758033847*$(x[4]) + 0.052390589356133166*$(x[5]))) + -0.9099166873549671*log(1 + exp(0.7510327333371474 + 0.48933101371416265*$(x[1]) + 0.5299590321169751*$(x[2]) + 0.05546767568522437*$(x[3]) + -0.4663274814822995*$(x[4]) + -0.6567285643500202*$(x[5]))) + 0.9082646170966706*log(1 + exp(-0.27157848416188113 + 0.560003468778171*$(x[1]) + 0.2933473130688298*$(x[2]) + 0.8808820816540863*$(x[3]) + 0.6726368959629845*$(x[4]) + -0.19953521813356856*$(x[5]))))) + -0.6322955384333278*log(1 + exp(-0.7118497063916944 + 0.36586202274128254*log(1 + exp(-0.02399021036854254 + -0.4943608904262171*$(x[1]) + -0.3028702711871105*$(x[2]) + 0.7730701583831108*$(x[3]) + -0.5656579689281509*$(x[4]) + 0.08835909373115314*$(x[5]))) + -0.9129427463398758*log(1 + exp(-0.6908874409241781 + -0.9493883109743697*$(x[1]) + 0.1034335229123764*$(x[2]) + 0.3901179347519834*$(x[3]) + -0.16073609758033847*$(x[4]) + 0.052390589356133166*$(x[5]))) + -0.5191517354247042*log(1 + exp(0.7510327333371474 + 0.48933101371416265*$(x[1]) + 0.5299590321169751*$(x[2]) + 0.05546767568522437*$(x[3]) + -0.4663274814822995*$(x[4]) + -0.6567285643500202*$(x[5]))) + 0.42477634545139953*log(1 + exp(-0.27157848416188113 + 0.560003468778171*$(x[1]) + 0.2933473130688298*$(x[2]) + 0.8808820816540863*$(x[3]) + 0.6726368959629845*$(x[4]) + -0.19953521813356856*$(x[5]))))) + 0.03598168048021799*log(1 + exp(0.8429853377330745 + -0.6377098680254001*log(1 + exp(-0.02399021036854254 + -0.4943608904262171*$(x[1]) + -0.3028702711871105*$(x[2]) + 0.7730701583831108*$(x[3]) + -0.5656579689281509*$(x[4]) + 0.08835909373115314*$(x[5]))) + 0.04468698484435718*log(1 + exp(-0.6908874409241781 + -0.9493883109743697*$(x[1]) + 0.1034335229123764*$(x[2]) + 0.3901179347519834*$(x[3]) + -0.16073609758033847*$(x[4]) + 0.052390589356133166*$(x[5]))) + 0.3607403638969684*log(1 + exp(0.7510327333371474 + 0.48933101371416265*$(x[1]) + 0.5299590321169751*$(x[2]) + 0.05546767568522437*$(x[3]) + -0.4663274814822995*$(x[4]) + -0.6567285643500202*$(x[5]))) + 0.4004266402265344*log(1 + exp(-0.27157848416188113 + 0.560003468778171*$(x[1]) + 0.2933473130688298*$(x[2]) + 0.8808820816540863*$(x[3]) + 0.6726368959629845*$(x[4]) + -0.19953521813356856*$(x[5]))))))) + log(1 + exp(-0.7045693885767101 + 0.765322888364711*log(1 + exp(0.9080159995418975 + 0.10769169268185319*log(1 + exp(-0.02399021036854254 + -0.4943608904262171*$(x[1]) + -0.3028702711871105*$(x[2]) + 0.7730701583831108*$(x[3]) + -0.5656579689281509*$(x[4]) + 0.08835909373115314*$(x[5]))) + 0.2597058016512084*log(1 + exp(-0.6908874409241781 + -0.9493883109743697*$(x[1]) + 0.1034335229123764*$(x[2]) + 0.3901179347519834*$(x[3]) + -0.16073609758033847*$(x[4]) + 0.052390589356133166*$(x[5]))) + 0.8597063403720888*log(1 + exp(0.7510327333371474 + 0.48933101371416265*$(x[1]) + 0.5299590321169751*$(x[2]) + 0.05546767568522437*$(x[3]) + -0.4663274814822995*$(x[4]) + -0.6567285643500202*$(x[5]))) + 0.17950756837409587*log(1 + exp(-0.27157848416188113 + 0.560003468778171*$(x[1]) + 0.2933473130688298*$(x[2]) + 0.8808820816540863*$(x[3]) + 0.6726368959629845*$(x[4]) + -0.19953521813356856*$(x[5]))))) + 0.9087169409161477*log(1 + exp(-0.9264374759333807 + -0.7028105280156569*log(1 + exp(-0.02399021036854254 + -0.4943608904262171*$(x[1]) + -0.3028702711871105*$(x[2]) + 0.7730701583831108*$(x[3]) + -0.5656579689281509*$(x[4]) + 0.08835909373115314*$(x[5]))) + 0.6656745963467641*log(1 + exp(-0.6908874409241781 + -0.9493883109743697*$(x[1]) + 0.1034335229123764*$(x[2]) + 0.3901179347519834*$(x[3]) + -0.16073609758033847*$(x[4]) + 0.052390589356133166*$(x[5]))) + -0.9099166873549671*log(1 + exp(0.7510327333371474 + 0.48933101371416265*$(x[1]) + 0.5299590321169751*$(x[2]) + 0.05546767568522437*$(x[3]) + -0.4663274814822995*$(x[4]) + -0.6567285643500202*$(x[5]))) + 0.9082646170966706*log(1 + exp(-0.27157848416188113 + 0.560003468778171*$(x[1]) + 0.2933473130688298*$(x[2]) + 0.8808820816540863*$(x[3]) + 0.6726368959629845*$(x[4]) + -0.19953521813356856*$(x[5]))))) + -0.29202574365958966*log(1 + exp(-0.7118497063916944 + 0.36586202274128254*log(1 + exp(-0.02399021036854254 + -0.4943608904262171*$(x[1]) + -0.3028702711871105*$(x[2]) + 0.7730701583831108*$(x[3]) + -0.5656579689281509*$(x[4]) + 0.08835909373115314*$(x[5]))) + -0.9129427463398758*log(1 + exp(-0.6908874409241781 + -0.9493883109743697*$(x[1]) + 0.1034335229123764*$(x[2]) + 0.3901179347519834*$(x[3]) + -0.16073609758033847*$(x[4]) + 0.052390589356133166*$(x[5]))) + -0.5191517354247042*log(1 + exp(0.7510327333371474 + 0.48933101371416265*$(x[1]) + 0.5299590321169751*$(x[2]) + 0.05546767568522437*$(x[3]) + -0.4663274814822995*$(x[4]) + -0.6567285643500202*$(x[5]))) + 0.42477634545139953*log(1 + exp(-0.27157848416188113 + 0.560003468778171*$(x[1]) + 0.2933473130688298*$(x[2]) + 0.8808820816540863*$(x[3]) + 0.6726368959629845*$(x[4]) + -0.19953521813356856*$(x[5]))))) + -0.14636407330293322*log(1 + exp(0.8429853377330745 + -0.6377098680254001*log(1 + exp(-0.02399021036854254 + -0.4943608904262171*$(x[1]) + -0.3028702711871105*$(x[2]) + 0.7730701583831108*$(x[3]) + -0.5656579689281509*$(x[4]) + 0.08835909373115314*$(x[5]))) + 0.04468698484435718*log(1 + exp(-0.6908874409241781 + -0.9493883109743697*$(x[1]) + 0.1034335229123764*$(x[2]) + 0.3901179347519834*$(x[3]) + -0.16073609758033847*$(x[4]) + 0.052390589356133166*$(x[5]))) + 0.3607403638969684*log(1 + exp(0.7510327333371474 + 0.48933101371416265*$(x[1]) + 0.5299590321169751*$(x[2]) + 0.05546767568522437*$(x[3]) + -0.4663274814822995*$(x[4]) + -0.6567285643500202*$(x[5]))) + 0.4004266402265344*log(1 + exp(-0.27157848416188113 + 0.560003468778171*$(x[1]) + 0.2933473130688298*$(x[2]) + 0.8808820816540863*$(x[3]) + 0.6726368959629845*$(x[4]) + -0.19953521813356856*$(x[5]))))))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
function utility(consumption::Float64,new_debt::Bool,reserve::Bool,tozero_l::Bool)

    #Put here Your Favorite Utility Function
    check1 = ((new_debt) | tozero_l)
    check2 = reserve
    check3 = consumption > 0
    ( check1 & check2 & check3) ? ( return( (consumption ^ (1 - σ)) / (1 - σ) ) ) : ( return(-Inf) )

end


text:
using Mimi                                          #Make sure to call Mimi again

@defcomp emissions begin
    regions     = Index()                           #The regions index must be specified for each component

    E           = Variable(index=[time, regions])   #Total greenhouse gas emissions
    E_Global    = Variable(index=[time])            #Global emissions (sum of regional emissions)
    sigma       = Parameter(index=[time, regions])  #Emissions output ratio
    YGROSS      = Parameter(index=[time, regions])  #Gross output - Note that YGROSS is now a parameter
end

function timestep(state::emissions, t::Int)
    v = state.Variables
    p = state.Parameters
    d = state.Dimensions

    #Define an eqation for E
    for r in d.regions
        v.E[t,r] = p.YGROSS[t,r] * p.sigma[t,r]
    end

    #Define an equation for E_Global
    for r in d.regions
        v.E_Global[t] = sum(v.E[t,:])
    end
end


text:
function adaptCrossover(p_cr::Vector{Float64}, M::Vector{Float64})
    p_best = p_cr[argmin(M)]

    for i = 1:length(p_cr)
        if M[i] > 0.3
            pn = abs(p_best .+ 0.3 * randn())
            if pn > 1.0
                pn = 1.0
            end

            if pn < 0.0
                pn = p_best
            end
            p_cr[i] = pn
        end
    end

    return p_cr
end

function resizePop!(P::Array, N_new::Int, K::Int)
    N = length(P)

    if N == N_new
        return P
    end

    f = zeros(N)
    for i = 1:N
        f[i] = P[i].f
    end

    ids = sortperm(P, lt = is_better)[1:N_new]
    return P[ids]
end


text:
@testset "tempscoll" begin

    n = 10
    ϵ = 0.1
    rng = MersenneTwister(42)

    p = SquareParticles(rng, n, ϵ)

    q = [[0.5       , 0.5       ],
         [0.20006862, 0.38139954],
         [0.32267135, 0.72687046],
         [0.68044674, 0.78095695],
         [0.77473372, 0.51989917],
         [0.4949439 , 0.23201753],
         [0.74204115, 0.28248482],
         [0.63101615, 0.3650337 ],
         [0.5846235 , 0.64552143],
         [0.45656816, 0.76644607]]

    v = [[ 0.,   1.],
         [ 1.,   0.],
         [-1.,   0.],
         [ 0.,  -1.],
         [ 1.,   0.],
         [ 0.,   1.],
         [ 0.,   1.],
         [-1.,   0.],
         [-1.,   0.],
         [ 0.,  -1.]]

    for i in 1:n
        p.q[i] = q[i]
        p.v[i] = v[i]
    end

    @test Uchiyama.compute_dt(p, 4, 1, 1) ≈ 0.130701845
    @test Uchiyama.compute_dt(p, 9, 1, 1) ≈ 0.015072465
    @test Uchiyama.compute_dt(p, 4, 2, 1) ≈ 0.339967765
    @test Uchiyama.compute_dt(p, 6, 2, 1) ≈ 0.122128645
    @test Uchiyama.compute_dt(p, 8, 2, 1) ≈ 0.123656685
    @test Uchiyama.compute_dt(p, 7, 3, 3) ≈ 0.412507920
    @test Uchiyama.compute_dt(p, 6, 4, 1) ≈ 0.267221131
    @test Uchiyama.compute_dt(p, 7, 4, 1) ≈ 0.180033270
    @test Uchiyama.compute_dt(p, 8, 5, 2) ≈ 0.405573952
    @test Uchiyama.compute_dt(p, 9, 5, 2) ≈ 0.367756020
    @test Uchiyama.compute_dt(p, 8, 6, 1) ≈ 0.034544210
    @test Uchiyama.compute_dt(p, 8, 7, 7) ≈ 0.885761939

end


text:
using PFKernels
using PFKernels.PowerSystem
using PFKernels.ForwardDiff

AT = PFKernels.AT

# Matlab-like residual function
function residual(V, Ybus, Sbus, pv, pq)
    # form mismatch vector
    mis = V .* conj(Ybus * V) - Sbus
    # form residual vector
    F = [real(mis[pv])
        real(mis[pq])
        imag(mis[pq]) ]
    return F
end

# Driver for dual numbers
function residual!(F, v_m, v_a,
                                ybus_re, ybus_im,
                                pinj, qinj, pv, pq, nbus)
    npv = length(pv)
    npq = length(pq)
    t1s{N} = ForwardDiff.Dual{Nothing,Float64, N} where N
    FT = t1s{2}
    V = Vector
    adF = V{FT}(undef, length(F))
    adv_m = V{FT}(undef, length(v_m))
    adv_a = V{FT}(undef, length(v_a))
    adybus_re_nzval = V{FT}(undef, length(ybus_re.nzval))
    adybus_re_colptr = V{Int64}(undef, length(ybus_re.colptr))
    adybus_re_rowval = V{Int64}(undef, length(ybus_re.rowval))
    adybus_im_nzval = V{FT}(undef, length(ybus_im.nzval))
    adybus_im_colptr = V{Int64}(undef, length(ybus_im.colptr))
    adybus_im_rowval = V{Int64}(undef, length(ybus_im.rowval))
    adpinj = V{FT}(undef, length(pinj))
    adqinj = V{FT}(undef, length(qinj))
    adpv = V{FT}(undef, length(pv))
    adpq = V{FT}(undef, length(pq))

    adF .= F
    adv_m .= v_m
    adv_a .= v_a
    adybus_re_nzval .= ybus_re.nzval
    adybus_im_nzval .= ybus_im.nzval
    adpinj .= pinj
    adqinj .= qinj
    adpv .= pv
    adpq .= pq

    dF = AT(adF)
    dv_m = AT(adv_m)
    dv_a = AT(adv_a)
    dybus_re_nzval = AT(adybus_re_nzval)
    dybus_re_colptr = AT(ybus_re.colptr)
    dybus_re_rowval = AT(ybus_re.rowval)
    dybus_im_nzval = AT(adybus_im_nzval)
    dybus_im_colptr = AT(ybus_im.colptr)
    dybus_im_rowval = AT(ybus_im.rowval)
    dpinj = AT(adpinj)
    dqinj = AT(adqinj)
    dpv = AT(pv)
    dpq = AT(pq)

    PFKernels.kernel!(PFKernels.backend, dF, dv_m, dv_a,
                    dybus_re_nzval, dybus_re_colptr, dybus_re_rowval,
                    dybus_im_nzval, dybus_im_colptr, dybus_im_rowval,
                    dpinj, dqinj, dpv, dpq, nbus)
    F .= ForwardDiff.value.(dF)
end

const PS = PowerSystem

# Inject GPU kernels
PFKernels.include(joinpath(dirname(@__FILE__), "kernels.jl"))

local_case = "case14.raw"
V, Ybus, Sbus, Sload, ref, pv, pq = PFKernels.loaddata(local_case)
npv = size(pv, 1)
npq = size(pq, 1)
F = zeros(Float64, npv + 2*npq)
F♯ = residual(V, Ybus, Sbus, pv, pq)
Vm = abs.(V)
Va = angle.(V)
ybus_re, ybus_im = PFKernels.Spmat{Vector{Int}, Vector{Float64}}(Ybus)
pbus = real(Sbus)
qbus = imag(Sbus)
residual!(F, Vm, Va,
    ybus_re, ybus_im,
    pbus, qbus, pv, pq, length(V))
@assert(F ≈ F♯)


text:
# Unit tests for implementations/nbodysimulator/nbodysimulator_base.jl

@testset "implementations/nbodysimulator/nbodysimulator_base.jl" begin
    # TODO
end


text:
module ObjFileBase

export ObjectHandle, SectionRef, SymbolRef, debugsections

export printfield, printentry, printfield_with_color, deref,
    sectionaddress, sectionoffset, sectionsize, sectionname,
    load_strtab, readmeta

import Base: read, seek, readbytes, position, show

########## ObjFileBase.jl - Basic shared functionality for all object files ####
#
#   This package provides basic functionality that is shared among object files.
#   It was written with ELF/MachO/COFF in mind, but should be easily adaptable
#   to other object files as well.
#
#   The family of Object file implementations have a somewhat different focus than
#   a lot of other ones. While certainly usable to write a linker or compiler,
#   they are more intended for interactive exploration and verification and some
#   of the design choices reflect this. Nevertheless, this does not mean that
#   performance should be unnecessarily compromised.
#
# # # Basic Concepts
#
# # ObjectHandle
#
#   Provides an abstraction over a specific instance of the object file.
#   Currently many implementations put the IO object to be read from in this
#   handle, which has proven convenient for interactive usage.
#
# # SectionIterator, SymbolIterator, etc.
#
#   Provides support for iterating over Symbols, Sections, Relocations, etc.
#   These should have a reference to the original ObjectHandle and read the
#   desired information on demand. Rather than returning the read datastructure
#   directly, it should return in in an appropriate *Ref object.
#
# # SetionRef, SymbolRef
#
#   Convenience objects (mostly for interactive usage) that contain a reference
#   to the original object handle as well as any information that may be needed
#   to get further information about the Referenced object that is not
#   explicitly encoded in the data structure itself. Examples may include:
#
#       - Offset in the file (especially for variable length records
#            where only the header is part of the datastructure)
#       - Symbol/Section index
#       - Etc.
#
################################################################################

# # # Tools

macro mustimplement(sig)
    fname = sig.args[1]
    arg1 = sig.args[2]
    if isa(arg1,Expr)
        arg1 = arg1.args[1]
    end
    :($(esc(sig)) = error(typeof($(esc(arg1))),
                          " must implement ", $(Expr(:quote,fname))))
end

################################# Interfaces ###################################

# # # ObjectHandle

abstract ObjectHandle

# External facing interface
readmeta{T<:ObjectHandle}(io::IO, ::Type{T}) =
    error("$T must register a readmeta method")

@mustimplement readheader(oh::ObjectHandle)

# Querying general properties
@mustimplement endianness(oh::ObjectHandle)

# ObjectHandle IO interface
@mustimplement seek(oh::ObjectHandle, args...)
@mustimplement position(oh::ObjectHandle)

##
#  These are parameterized on the type of Object Handle. An imaginary Foo
#  fileformat might declare:
#
#   module Foo
#       immutable FooHandle; ... end
#       immutable Section <: ObjFileBase.Section{FooHandle}
#       ...; end
#   end
##

abstract SectionRef{T<:ObjectHandle}
abstract Section{T<:ObjectHandle}

# The size of the actual data contained in the section. This should exclude any
# padding mandated by the file format e.g. due to alignment rules
@mustimplement sectionsize(section::Section)

# The offset of the section in the file
@mustimplement sectionoffset(section::Section)

# The address of the section in virtual memory
@mustimplement sectionaddress(section::Section)

# The name of the section
@mustimplement sectionname(section::SectionRef)

# Retrieving the actual section datastructure
@mustimplement deref(section::SectionRef)

# Retrieving the object handle
@mustimplement handle(section::SectionRef)

abstract SymbolRef{T<:ObjectHandle}
abstract SymtabEntry{T<:ObjectHandle}

sectionsize(x::SectionRef) = sectionsize(deref(x))
sectionaddress(x::SectionRef) = sectionaddress(deref(x))
sectionoffset(x::SectionRef) = sectionoffset(deref(x))

handleT{T}(::Union(Type{SectionRef{T}}, Type{Section{T}}, Type{SymbolRef{T}},
    Type{SymtabEntry{T}})) = T

abstract StrTab

function load_strtab
end
@mustimplement strtab_lookup(s::StrTab, offset)

################################# Utilities ####################################

typealias SectionOrRef{T} Union(Section{T},SectionRef{T})

sectionsize(sect::SectionRef) = sectionsize(deref(sect))
sectionoffset(sect::SectionRef) = sectionoffset(deref(sect))

seek{T<:ObjectHandle,S}(oh::T, section::Section{S}) =
    (@assert T <: S; seek(oh,sectionoffset(section)))

seek(section::SectionRef) = seek(handle(section), deref(section))

function readbytes{T<:ObjectHandle,S}(oh::T,sec::Section{S})
    @assert T <: S
    seek(oh, sec)
    readbytes(oh, sectionsize(sec))
end
readbytes(sec::SectionRef) = readbytes(handle(sec),deref(sec))

typealias Maybe{T} Union(T,Nothing)

# # # Higher level debug info support
immutable DebugSections{T<:ObjectHandle, Sect}
    oh::T
    debug_abbrev::Maybe{Sect}
    debug_aranges::Maybe{Sect}
    debug_frame::Maybe{Sect}
    debug_info::Maybe{Sect}
    debug_line::Maybe{Sect}
    debug_loc::Maybe{Sect}
    debug_macinfo::Maybe{Sect}
    debug_pubnames::Maybe{Sect}
    debug_ranges::Maybe{Sect}
    debug_str::Maybe{Sect}
    debug_types::Maybe{Sect}

end

function DebugSections{T}(oh::T; debug_abbrev = nothing, debug_aranges = nothing,
    debug_frame = nothing, debug_info = nothing, debug_line = nothing,
    debug_macinfo = nothing, debug_pubnames = nothing, debug_loc= nothing,
    debug_ranges = nothing, debug_str = nothing, debug_types = nothing)
    DebugSections(oh, debug_abbrev, debug_aranges, debug_frame, debug_info,
        debug_line, debug_loc, debug_macinfo, debug_pubnames, debug_ranges,
        debug_str, debug_types)
end

function DebugSections{T<:ObjectHandle}(oh::T, sections::Dict)
    DebugSections(oh,
        debug_abbrev = get(sections, "debug_abbrev", nothing),
        debug_aranges = get(sections, "debug_aranges", nothing),
        debug_frame = get(sections, "debug_frame", nothing),
        debug_info = get(sections, "debug_info", nothing),
        debug_line = get(sections, "debug_line", nothing),
        debug_loc = get(sections, "debug_loc", nothing),
        debug_macinfo = get(sections, "debug_macinfo", nothing),
        debug_pubnames = get(sections, "debug_pubnames", nothing),
        debug_ranges = get(sections, "debug_ranges", nothing),
        debug_str = get(sections, "debug_str", nothing),
        debug_types = get(sections, "debug_types", nothing))
end

function show(io::IO, dsect::DebugSections)
    println(io, "Debug Sections for $(dsect.oh)")
    println(io,"========================= debug_abbrev =========================")
    println(io,dsect.debug_abbrev)
    println(io,"======================== debug_aranges =========================")
    println(io,dsect.debug_aranges)
    println(io,"========================= debug_frame ==========================")
    println(io,dsect.debug_frame)
    println(io,"========================= debug_info ===========================")
    println(io,dsect.debug_info)
    println(io,"========================= debug_line ===========================")
    println(io,dsect.debug_line)
    println(io,"========================== debug_loc ===========================")
    println(io,dsect.debug_loc)
    println(io,"======================== debug_macinfo =========================")
    println(io,dsect.debug_macinfo)
    println(io,"======================= debug_pubnames =========================")
    println(io,dsect.debug_pubnames)
    println(io,"======================== debug_ranges ==========================")
    println(io,dsect.debug_ranges)
    println(io,"=========================== debug_str ==========================")
    println(io,dsect.debug_str)
    println(io,"========================= debug_types ==========================")
    println(io,dsect.debug_types)
end

@mustimplement debugsections(oh::ObjectHandle)

export findindexbyname, findcubyname

function findindexbyname
end

function findcubyname
end

# Utils
# JIT Utils

function is_jit_section(s::Section)
    sectionaddress(s) > 0x100000
end
is_jit_section(s::SectionRef) = is_jit_section(deref(s))

@mustimplement replace_sections_from_memory(h::ObjectHandle, buffer)

# Printing utils
function printfield(io::IO,string,fieldlength; align = :right)
    (align == :right) && print(io," "^max(fieldlength-length(string),0))
    print(io,string)
    (align == :left) && print(io," "^max(fieldlength-length(string),0))
end
function printfield_with_color(color,io::IO,string,fieldlength; align = :right)
    (align == :right) && print(io," "^max(fieldlength-length(string),0))
    print_with_color(color,io,string)
    (align == :left) && print(io," "^max(fieldlength-length(string),0))
end
printentry(io::IO,header,values...) = (printfield(io,header,21);println(io," ",values...))

# User facing interfaces

immutable MagicMismatch;
    message
end

function readmeta(io::IO)
    ts = subtypes(ObjFileBase.ObjectHandle)
    pos = position(io)
    for T in ts
        seek(io,pos)
        try
            return readmeta(io, T)
        catch e
            if !isa(e,MagicMismatch)
                rethrow(e)
            end
        end
    end
    error("""
        Object file is not any of $(join(ts, ", "))!
        To force one object file use readmeta(io,T).
        If the format you want is not listed, make sure
        the appropriate pacakge is loaded before calling
        this function.
        """)
end

readmeta(file::String) = readmeta(open(file,"r"))

# Others

function getSectionLoadAddress
end

end # module

max_stars_repo_name: ambujone/DIDebug

max_issues_repo_name: ambujone/DIDebug

max_forks_repo_name: ambujone/DIDebug


text:
getrange(t, n) = range(minimum(t); stop = maximum(t), length = n)
itp(x) = LinearInterpolation(1:length(x), x)
interp(itp, tᵢ) = itp()

"""
Linearly interpolates two vectors x and y on a linear grid consisting of `nsteps`.
"""
function interp(x::Vector, y::Vector, nsteps::Int)
    # Interpolate
    itp = interpolate((x,), y, Gridded(Linear()))
    # Interpolate at the given resolution
    x_fills = LinRange(minimum(x), maximum(x), nsteps)
    y_fills = itp(x_fills)

    return collect(x_fills), y_fills
end

"""
Linearly interpolates two vector x and y on a linear grid consisting of `nsteps`.
"""
function interp(x, y, range::LinRange)
    itp = interpolate((x,), y, Gridded(Linear()))

    # Interpolate at the given resolution
    return itp(range)
end

"""
    interp!(ȳ::Vector, itp)

Interpolate using the pre-computed interpolation instance `itp` into
the pre-allocated vector `ȳ`.
"""
function interp!(ȳ::Vector, itp)
    y_fills .= itp(x_fills)

    return collect(x_fills), y_fills
end


text:
"""
    FEMMDeforLinearESNICEModule

Formulation for the small displacement, small strain deformation
model for Nodally-Integrated Continuum Elements (NICE).

The approximation is  originally from Dohrmann et al IJNME 47 (2000).
The formulation was subsequently developed in Krysl, P. and Zhu, B.
Locking-free continuum displacement finite elements with nodal
integration, International Journal for Numerical Methods in Engineering,
76,7,1020-1043,2008.

The stabilization scheme comes from papers on energy-sampling stabilization
for mean-strain elements (Krysl).
"""
module FEMMDeforLinearESNICEModule

__precompile__(true)

using FinEtools.FTypesModule: FInt, FFlt, FCplxFlt, FFltVec, FIntVec, FFltMat, FIntMat, FMat, FVec, FDataDict
import FinEtools.FENodeSetModule: FENodeSet
import FinEtools.FESetModule: AbstractFESet, FESetH8, FESetT4, manifdim, nodesperelem, gradN!
import FinEtools.IntegDomainModule: IntegDomain, integrationdata, Jacobianvolume
import FinEtoolsDeforLinear.FEMMDeforLinearBaseModule: AbstractFEMMDeforLinear
import FinEtoolsDeforLinear.DeforModelRedModule: AbstractDeforModelRed, DeforModelRed3D
import FinEtoolsDeforLinear.MatDeforLinearElasticModule: AbstractMatDeforLinearElastic, tangentmoduli!, update!, thermalstrain!
import FinEtoolsDeforLinear.MatDeforElastIsoModule: MatDeforElastIso
import FinEtools.FieldModule: ndofs, gatherdofnums!, gatherfixedvalues_asvec!, gathervalues_asvec!, gathervalues_asmat!
import FinEtools.NodalFieldModule: NodalField, nnodes
import FinEtools.CSysModule: CSys, updatecsmat!
import FinEtools.FENodeToFEMapModule: FENodeToFEMap
import FinEtoolsDeforLinear.DeforModelRedModule: nstressstrain, nthermstrain, Blmat!, divmat, vgradmat
import FinEtools.AssemblyModule: AbstractSysvecAssembler, AbstractSysmatAssembler, SysmatAssemblerSparseSymm, startassembly!, assemble!, makematrix!, makevector!, SysvecAssembler
using FinEtools.MatrixUtilityModule: add_btdb_ut_only!, complete_lt!, loc!, jac!, locjac!, adjugate3!
import FinEtoolsDeforLinear.FEMMDeforLinearBaseModule: stiffness, nzebcloadsstiffness, mass, thermalstrainloads, inspectintegpoints
import FinEtools.FEMMBaseModule: associategeometry!
import FinEtoolsDeforLinear.MatDeforModule: rotstressvec!
import FinEtools.MatModule: massdensity
import LinearAlgebra: mul!, Transpose, UpperTriangular, eigvals, det
At_mul_B!(C, A, B) = mul!(C, Transpose(A), B)
A_mul_B!(C, A, B) = mul!(C, A, B)
import LinearAlgebra: norm, qr, diag, dot, cond, I, cross
import Statistics: mean
import StatsBase: geomean

"""
    AbstractFEMMDeforLinearESNICE <: AbstractFEMMDeforLinear

Abstract FEMM type for Nodally Integrated Continuum Elements (NICE) with energy-sampling stabilization.
"""
abstract type AbstractFEMMDeforLinearESNICE <: AbstractFEMMDeforLinear end

# Tetrahedron
# The coefficient set below was obtained by fitting the ratio of energies
# true/approximate for the finite element model of six tetrahedra arranged
# into a rectangular block and subject to pure bending
# Fitting for a small aspect-ratio range (1.0 to 10)
_T4_stabilization_parameters = (2.101588423297799,  1.311321055432958)

mutable struct _NodalBasisFunctionGradients
    gradN::FFltMat
    patchconn::FIntVec
    Vpatch::FFlt
end

function _make_stabilization_material(material::M) where {M}
    ns = fieldnames(typeof(material))
    E = 0.0; nu = 0.0
    if :E in ns
        E = material.E
        if material.nu < 0.3
            nu = material.nu
        else
            nu = 0.3 + (material.nu - 0.3) / 2.0
        end
    else
        if :E1 in ns
            E = mean([material.E1, material.E2, material.E3])
            nu = min(material.nu12, material.nu13, material.nu23)
        else
            error("No clues on how to construct the stabilization material")
        end
    end
    return  MatDeforElastIso(material.mr, 0.0, E, nu, 0.0)
end

"""
    FEMMDeforLinearESNICET4{MR<:AbstractDeforModelRed, S<:FESetT4, F<:Function, M<:AbstractMatDeforLinearElastic} <: AbstractFEMMDeforLinearESNICE

FEMM type for Nodally Integrated Continuum Elements (NICE) based on the 4-node tetrahedron.
"""
mutable struct FEMMDeforLinearESNICET4{MR<:AbstractDeforModelRed, S<:FESetT4, F<:Function, M<:AbstractMatDeforLinearElastic} <: AbstractFEMMDeforLinearESNICE
    mr::Type{MR}
    integdomain::IntegDomain{S, F} # geometry data
    mcsys::CSys # updater of the material orientation matrix
    material::M # material object
    stabilization_material::MatDeforElastIso
    nodalbasisfunctiongrad::Vector{_NodalBasisFunctionGradients}
    ephis::Vector{FFlt}
    nphis::Vector{FFlt}
end

"""
    FEMMDeforLinearESNICEH8{MR<:AbstractDeforModelRed, S<:FESetH8, F<:Function, M<:AbstractMatDeforLinearElastic} <: AbstractFEMMDeforLinearESNICE

FEMM type for Nodally Integrated Continuum Elements (NICE) based on the 8-node
hexahedron.
"""
mutable struct FEMMDeforLinearESNICEH8{MR<:AbstractDeforModelRed, S<:FESetH8, F<:Function, M<:AbstractMatDeforLinearElastic} <: AbstractFEMMDeforLinearESNICE
    mr::Type{MR}
    integdomain::IntegDomain{S, F} # geometry data
    mcsys::CSys # updater of the material orientation matrix
    material::M # material object
    stabilization_material::MatDeforElastIso
    nodalbasisfunctiongrad::Vector{_NodalBasisFunctionGradients}
    ephis::Vector{FFlt}
    nphis::Vector{FFlt}
end

function FEMMDeforLinearESNICET4(mr::Type{MR}, integdomain::IntegDomain{S, F}, mcsys::CSys, material::M) where {MR<:AbstractDeforModelRed,  S<:FESetT4, F<:Function, M<:AbstractMatDeforLinearElastic}
    @assert mr == material.mr "Model reduction is mismatched"
    @assert (mr == DeforModelRed3D) "3D model required"
    stabilization_material = _make_stabilization_material(material)
    return FEMMDeforLinearESNICET4(mr, integdomain, mcsys, material, stabilization_material, _NodalBasisFunctionGradients[], fill(zero(FFlt), 1), fill(zero(FFlt), 1))
end

function FEMMDeforLinearESNICET4(mr::Type{MR}, integdomain::IntegDomain{S, F}, material::M) where {MR<:AbstractDeforModelRed,  S<:FESetT4, F<:Function, M<:AbstractMatDeforLinearElastic}
    @assert mr == material.mr "Model reduction is mismatched"
    @assert (mr == DeforModelRed3D) "3D model required"
    stabilization_material = _make_stabilization_material(material)
    return FEMMDeforLinearESNICET4(mr, integdomain, CSys(manifdim(integdomain.fes)), material, stabilization_material, _NodalBasisFunctionGradients[], fill(zero(FFlt), 1), fill(zero(FFlt), 1))
end

function FEMMDeforLinearESNICEH8(mr::Type{MR}, integdomain::IntegDomain{S, F}, mcsys::CSys, material::M) where {MR<:AbstractDeforModelRed,  S<:FESetH8, F<:Function, M<:AbstractMatDeforLinearElastic}
    @assert mr == material.mr "Model reduction is mismatched"
    @assert (mr == DeforModelRed3D) "3D model required"
    stabilization_material = _make_stabilization_material(material)
    return FEMMDeforLinearESNICEH8(mr, integdomain, mcsys, material, stabilization_material, _NodalBasisFunctionGradients[], fill(zero(FFlt), 1), fill(zero(FFlt), 1))
end

function FEMMDeforLinearESNICEH8(mr::Type{MR}, integdomain::IntegDomain{S, F}, material::M) where {MR<:AbstractDeforModelRed,  S<:FESetH8, F<:Function, M<:AbstractMatDeforLinearElastic}
    @assert mr == material.mr "Model reduction is mismatched"
    @assert (mr == DeforModelRed3D) "3D model required"
    stabilization_material = _make_stabilization_material(material)
    return FEMMDeforLinearESNICEH8(mr, integdomain, CSys(manifdim(integdomain.fes)), material, stabilization_material, _NodalBasisFunctionGradients[], fill(zero(FFlt), 1), fill(zero(FFlt), 1))
end

function _buffers1(self::AbstractFEMMDeforLinearESNICE, geom::NodalField)
    fes = self.integdomain.fes
    nne = nodesperelem(fes); # number of nodes for element
    sdim = ndofs(geom);            # number of space dimensions
    mdim = manifdim(fes); # manifold dimension of the element
    # Prepare buffers
    loc = fill(zero(FFlt), 1, sdim); # quadrature point location -- buffer
    J = fill(zero(FFlt), sdim, mdim); # Jacobian matrix -- buffer
    adjJ = fill(zero(FFlt), sdim, mdim); # Jacobian matrix -- buffer
    csmatTJ = fill(zero(FFlt), mdim, mdim); # intermediate result -- buffer
    gradN = fill(zero(FFlt), nne, mdim);
    xl = fill(zero(FFlt), nne, mdim);
    return loc, J, adjJ, csmatTJ, gradN, xl
end

function _buffers2(self::AbstractFEMMDeforLinearESNICE, geom::NodalField, u::NodalField, npts::FInt)
    fes = self.integdomain.fes
    ndn = ndofs(u); # number of degrees of freedom per node
    nne = nodesperelem(fes); # number of nodes for element
    sdim = ndofs(geom);            # number of space dimensions
    mdim = manifdim(fes); # manifold dimension of the element
    nstrs = nstressstrain(self.mr);  # number of stresses
    elmatdim = ndn*nne;             # dimension of the element matrix
    # Prepare buffers
    ecoords = fill(zero(FFlt), nne, ndofs(geom)); # array of Element coordinates
    dofnums = zeros(FInt, elmatdim); # degree of freedom array -- buffer
    loc = fill(zero(FFlt), 1, sdim); # quadrature point location -- buffer
    J = fill(zero(FFlt), sdim, mdim); # Jacobian matrix -- buffer
    csmatTJ = fill(zero(FFlt), mdim, mdim); # intermediate result -- buffer
    Jac = fill(zero(FFlt), npts);
    D = fill(zero(FFlt), nstrs, nstrs); # material stiffness matrix -- buffer
    Dstab = fill(zero(FFlt), nstrs, nstrs); # material stiffness matrix -- buffer
    elmat = fill(zero(FFlt), elmatdim, elmatdim);      # element matrix -- buffer
    B = fill(zero(FFlt), nstrs, elmatdim); # strain-displacement matrix -- buffer
    return ecoords, dofnums, loc, J, csmatTJ, Jac, D, Dstab, elmat, B
end

function _buffers3(self::AbstractFEMMDeforLinearESNICE, geom::NodalField, u::NodalField)
    fes = self.integdomain.fes
    ndn = ndofs(u); # number of degrees of freedom per node
    nne = nodesperelem(fes); # number of nodes for element
    sdim = ndofs(geom);            # number of space dimensions
    mdim = manifdim(fes); # manifold dimension of the element
    nstrs = nstressstrain(self.mr);  # number of stresses
    elmatdim = ndn*nne;             # dimension of the element matrix
    # Prepare buffers
    elmat = fill(zero(FFlt), elmatdim, elmatdim);      # element matrix -- buffer
    dofnums = zeros(FInt, elmatdim); # degree of freedom array -- buffer
    B = fill(zero(FFlt), nstrs, elmatdim); # strain-displacement matrix -- buffer
    DB = fill(zero(FFlt), nstrs, elmatdim); # strain-displacement matrix -- buffer
    elvecfix = fill(zero(FFlt), elmatdim); # vector of prescribed displ. -- buffer
    elvec = fill(zero(FFlt), elmatdim); # element vector -- buffer
    gradN = fill(zero(FFlt), nne, mdim); # intermediate result -- buffer
    return dofnums, B, DB, elmat, elvec, elvecfix, gradN
end

function _patchconn(fes, gl, thisnn)
    # Generate patch connectivity for a given node (thisnn)
    # from the connectivities of the finite elements attached to it.
    return vcat(collect(setdiff(Set([i for j=1:length(gl) for i in fes.conn[gl[j]]]), thisnn)), [thisnn])
end

function _computenodalbfungrads(self, geom)
    # # Compute the nodal basis function gradients.
    # # Return the cell array of structures with attributes
    # %      bfun_gradients{nix}.Nspd= basis function gradient matrix
    # #        bfun_gradients{nix}.Vpatch= nodal patch volume
    # #        bfun_gradients{nix}.patchconn= nodal patch connectivity

    fes = self.integdomain.fes
    npts,  Ns,  gradNparams,  w,  pc = integrationdata(self.integdomain);
    loc, J, adjJ, csmatTJ, gradN, xl = _buffers1(self, geom)

    # Get the inverse map from finite element nodes to geometric cells
    fen2fe = FENodeToFEMap(fes.conn, nnodes(geom));
    # Initialize the nodal gradients, nodal patch, and patch connectivity
    bfungrads = fill(_NodalBasisFunctionGradients(fill(0.0, 0, 0), fill(0, 0), 0.0), nnodes(geom));
    # Now loop over all finite element nodes in the map
    lnmap = fill(0, length(fen2fe.map)); # Local node map: buffer to speed up operations
    for nix = 1:length(fen2fe.map)
        gl = fen2fe.map[nix];
        thisnn = nix; # We are at this node
        if !isempty(gl) # This node has an element patch in this block
            # establish local numbering of all nodes of the patch @ node thisnn
            p = _patchconn(fes, gl, thisnn);
            np = length(p);
            lnmap[p] .= 1:np;# now store the local numbers
            c = reshape(geom.values[thisnn, :], 1, ndofs(geom))
            updatecsmat!(self.mcsys, c, J, 0);
            gradNavg = fill(0.0, np, ndofs(geom));# preallocate strain-displacement matrix
            Vpatch = 0.0;
            for k = 1:length(gl)
                i = gl[k]
                kconn = collect(fes.conn[i]);
                pci = findfirst(cx -> cx == thisnn, kconn);# at which node in the element are we with this quadrature point?
                @assert 1 <= pci <= nodesperelem(fes)
                # centered coordinates of nodes in the material coordinate system
                for cn = 1:length(kconn)
                    xl[cn, :] = (reshape(geom.values[kconn[cn], :], 1, ndofs(geom)) - c) * self.mcsys.csmat
                end
                jac!(J, xl, gradNparams[pci])
                At_mul_B!(csmatTJ, self.mcsys.csmat, J); # local Jacobian matrix
                Jac = Jacobianvolume(self.integdomain, J, c, fes.conn[i], Ns[pci]);
                Vpatch += Jac * w[pci];
                sgradN = gradNparams[pci] * adjugate3!(adjJ, J);
                gradNavg[lnmap[kconn],:] += (w[pci] .* sgradN);
            end
            @assert Vpatch != 0
            gradNavg ./= Vpatch;
            bfungrads[nix] = _NodalBasisFunctionGradients(gradNavg, p, Vpatch);
            lnmap[p] .= 0; # Restore the buffer to pristine condition
        end
    end
    self.nodalbasisfunctiongrad = bfungrads
    return self
end

function _tetaspectratiovol(X)
    edge1 = vec(X[2, :] - X[1, :])
    edge2 = vec(X[3, :] - X[1, :])
    edge3 = vec(X[4, :] - X[1, :])
    edge4 = vec(X[3, :] - X[2, :])
    edge5 = vec(X[4, :] - X[3, :])
    edge6 = vec(X[4, :] - X[2, :])
    V = dot(edge3, cross(edge1, edge2))
    A1 = norm(cross(edge1, edge2)) # This is twice the area of the triangle
    A2 = norm(cross(edge2, edge3))
    A3 = norm(cross(edge3, edge1))
    A4 = norm(cross(edge4, edge6))
    #h1, h2, h3, h4 = V/A1, V/A2, V/A3, V/A4
    L1, L2, L3, L4, L5, L6 = norm(edge1), norm(edge2), norm(edge3), norm(edge4), norm(edge5), norm(edge6)
    # the heights and the edge lengths will now be used to derive a composite
    # index: aspect ratio index. If this cannot be computed without generating
    # not-a-number or infinity, this number is assumed to be better represented
    # with 1.0 (perfect ratio), so that the element shape is then governed by
    # the other aspect ratio values.
    f = maximum
    f([L1, L2, L4]) != 0.0 && A1 != 0.0 ? ar1 = V/A1/f([L1, L2, L4]) : ar1 = 1.0
    f([L3, L2, L5]) != 0.0 && A2 != 0.0 ? ar2 = V/A2/f([L3, L2, L5]) : ar2 = 1.0
    f([L1, L3, L6]) != 0.0 && A3 != 0.0 ? ar3 = V/A3/f([L1, L3, L6]) : ar3 = 1.0
    f([L6, L5, L4]) != 0.0 && A4 != 0.0 ? ar4 = V/A4/f([L6, L5, L4]) : ar4 = 1.0
    return ar1, ar2, ar3, ar4, V/6
end

"""
    associategeometry!(self::F,  geom::NodalField{FFlt}) where {F<:FEMMDeforLinearESNICET4}

Associate geometry field with the FEMM.

Compute the  correction factors to account for  the shape of the  elements.
"""
function associategeometry!(self::F,  geom::NodalField{FFlt}; stabilization_parameters = _T4_stabilization_parameters) where {F<:FEMMDeforLinearESNICET4}
    (a, b) = stabilization_parameters
    fes = self.integdomain.fes
    self.ephis = fill(zero(FFlt), count(fes))
    evols = fill(zero(FFlt), count(fes))
    self.nphis = fill(zero(FFlt), nnodes(geom))
    nvols = fill(zero(FFlt), nnodes(geom))
    for i = 1:count(fes) # Loop over elements
        ar1, ar2, ar3, ar4, V = _tetaspectratiovol(geom.values[collect(fes.conn[i]), :])
        evols[i] = V;
        # If the aspect ratios are not reasonable, such as when the element is a
        # sliver or inverted, we turn off the stabilization for the element
        # by setting its stabilization factor to zero.
        if min(ar1, ar2, ar3, ar4) <= 0
            self.ephis[i] = 0.0
        else
            self.ephis[i] = (1.0 / (b * min(ar1, ar2, ar3, ar4) ^a) + 1.0) ^(-1)
        end
        # Accumulate: the stabilization factor at the node is the weighted mean
        # of the stabilization factors of the elements at that node
        for k = 1:nodesperelem(fes)
            nvols[fes.conn[i][k]] += evols[i]
            self.nphis[fes.conn[i][k]] += self.ephis[i] * evols[i]
        end
    end # Loop over elements
    @assert any(isnan.(self.ephis)) == false
    # Now scale the values at the nodes with the nodal volumes
    for k = 1:length(nvols)
        if nvols[k] != 0.0
            self.nphis[k] /= nvols[k]
        end
    end
    @assert any(isnan.(self.nphis)) == false
    # Now calculate the nodal basis function gradients
    return _computenodalbfungrads(self, geom)
end

"""
    associategeometry!(self::F,  geom::NodalField{FFlt}) where {F<:FEMMDeforLinearESNICEH8}

Associate geometry field with the FEMM.

Compute the  correction factors to account for  the shape of the  elements.
"""
function associategeometry!(self::F,  geom::NodalField{FFlt}) where {F<:FEMMDeforLinearESNICEH8}
    fes = self.integdomain.fes
    self.ephis = fill(zero(FFlt), count(fes))
    evols = fill(zero(FFlt), count(fes))
    self.nphis = fill(zero(FFlt), nnodes(geom))
    nvols = fill(zero(FFlt), nnodes(geom))
    npts, Ns, gradNparams, w, pc = integrationdata(self.integdomain);

    for i = 1:count(fes) # Loop over elements
        X = geom.values[collect(fes.conn[i]), :]
        V = 0.0;
        for j = 1:npts
            J = X' * gradNparams[j]
            Jac = det(J)
            @assert Jac > 0 "Nonpositive Jacobian"
            V = V + Jac * w[j]
            hs = sum(J .* J, dims=1)
            phi = 2*(1.0 + self.stabilization_material.nu)*minimum(hs)/maximum(hs)
            self.ephis[i] = max(self.ephis[i], phi/(1+phi))
        end
        evols[i] = V;
        # Accumulate: the stabilization factor at the node is the weighted mean of the stabilization factors of the elements at that node
        for k = 1:nodesperelem(fes)
            nvols[fes.conn[i][k]] += evols[i]
            self.nphis[fes.conn[i][k]] += self.ephis[i] * evols[i]
        end
    end # Loop over elements
    # Now scale the values at the nodes with the nodal volumes
    for k = 1:length(nvols)
        self.nphis[k] /= nvols[k]
    end
    # Now calculate the nodal basis function gradients
    return _computenodalbfungrads(self, geom)
end

"""
    stiffness(self::AbstractFEMMDeforLinearESNICE, assembler::A,
      geom::NodalField{FFlt},
      u::NodalField{T}) where {A<:AbstractSysmatAssembler, T<:Number}

Compute and assemble  stiffness matrix.
"""
function stiffness(self::AbstractFEMMDeforLinearESNICE, assembler::A, geom::NodalField{FFlt}, u::NodalField{T}) where {A<:AbstractSysmatAssembler, T<:Number}
    fes = self.integdomain.fes
    npts,  Ns,  gradNparams,  w,  pc = integrationdata(self.integdomain);
    ecoords, dofnums, loc, J, csmatTJ, Jac, D, Dstab = _buffers2(self, geom, u, npts)
    tangentmoduli!(self.material, D, 0.0, 0.0, loc, 0)
    tangentmoduli!(self.stabilization_material, Dstab, 0.0, 0.0, loc, 0)
    elmatsizeguess = 4*nodesperelem(fes)*ndofs(u)
    startassembly!(assembler, elmatsizeguess, elmatsizeguess, nnodes(u), u.nfreedofs, u.nfreedofs);
    for nix = 1:length(self.nodalbasisfunctiongrad)
        gradN = self.nodalbasisfunctiongrad[nix].gradN
        patchconn = self.nodalbasisfunctiongrad[nix].patchconn
        Vpatch = self.nodalbasisfunctiongrad[nix].Vpatch
        c = reshape(geom.values[nix, :], 1, ndofs(geom))
        updatecsmat!(self.mcsys, c, J, 0);
        nd = length(patchconn) * ndofs(u)
        Bnodal = fill(0.0, size(D, 1), nd)
        Blmat!(self.mr, Bnodal, Ns[1], gradN, c, self.mcsys.csmat);
        elmat = fill(0.0, nd, nd) # Can we SPEED it UP?
        DB = fill(0.0, size(D, 1), nd)
        add_btdb_ut_only!(elmat, Bnodal, Vpatch, D, DB)
        add_btdb_ut_only!(elmat, Bnodal, -self.nphis[nix]*Vpatch, Dstab, DB)
        complete_lt!(elmat)
        dofnums = fill(0, nd)
        gatherdofnums!(u, dofnums, patchconn); # retrieve degrees of freedom
        assemble!(assembler, elmat, dofnums, dofnums); # assemble symmetric matrix
    end # Loop over elements
    Kn = makematrix!(assembler);
    dofnums, B, DB, elmat, elvec, elvecfix, gradN = _buffers3(self, geom, u)
    # OPTIMIZATION: switch to a single-point quadrature rule here
    startassembly!(assembler, nodesperelem(fes)*ndofs(u), nodesperelem(fes)*ndofs(u), count(fes), u.nfreedofs, u.nfreedofs);
    for i = 1:count(fes) # Loop over elements
        gathervalues_asmat!(geom, ecoords, fes.conn[i]);
        fill!(elmat,  0.0); # Initialize element matrix
        for j = 1:npts # Loop over quadrature points
            locjac!(loc, J, ecoords, Ns[j], gradNparams[j])
            Jac = Jacobianvolume(self.integdomain, J, loc, fes.conn[i], Ns[j]);
            # Do the  following only if the element is well shaped and the
            # stabilization factor is positive; if the element is so distorted
            # that its Jacobian is non-positive, skip the following step.
            if self.ephis[i] > 0  && Jac != 0.0
                updatecsmat!(self.mcsys, loc, J, fes.label[i]);
                At_mul_B!(csmatTJ, self.mcsys.csmat, J); # local Jacobian matrix
                gradN!(fes, gradN, gradNparams[j], csmatTJ);
                Blmat!(self.mr, B, Ns[j], gradN, loc, self.mcsys.csmat);
                add_btdb_ut_only!(elmat, B, self.ephis[i]*Jac*w[j], Dstab, DB)
            end
        end # Loop over quadrature points
        complete_lt!(elmat)
        gatherdofnums!(u, dofnums, fes.conn[i]); # retrieve degrees of freedom
        assemble!(assembler, elmat, dofnums, dofnums); # assemble symmetric matrix
    end # Loop over elements
    return makematrix!(assembler) + Kn;
end

function stiffness(self::AbstractFEMMDeforLinearESNICE, geom::NodalField{FFlt},  u::NodalField{T}) where {T<:Number}
    assembler = SysmatAssemblerSparseSymm();
    return stiffness(self, assembler, geom, u);
end

"""
nzebcloadsstiffness(self::AbstractFEMMDeforLinear,  assembler::A,
  geom::NodalField{FFlt},
  u::NodalField{T}) where {A<:AbstractSysvecAssembler, T<:Number}

Compute load vector for nonzero EBC for fixed displacement.
"""
function nzebcloadsstiffness(self::AbstractFEMMDeforLinearESNICE,  assembler::A, geom::NodalField{FFlt}, u::NodalField{T}) where {A<:AbstractSysvecAssembler, T<:Number}
    error("Not implemented yet")
end

function nzebcloadsstiffness(self::AbstractFEMMDeforLinearESNICE, geom::NodalField{FFlt}, u::NodalField{T}) where {T<:Number}
    assembler = SysvecAssembler()
    return  nzebcloadsstiffness(self, assembler, geom, u);
end

"""
   inspectintegpoints(self::AbstractFEMMDeforLinearESNICE, geom::NodalField{FFlt},  u::NodalField{T}, dT::NodalField{FFlt}, felist::FIntVec, inspector::F,  idat, quantity=:Cauchy; context...) where {T<:Number, F<:Function}

Inspect integration point quantities.

# Arguments
- `geom` - reference geometry field
- `u` - displacement field
- `dT` - temperature difference field
- `felist` - indexes of the finite elements that are to be inspected:
    The fes to be included are: `fes[felist]`.
- `context`    - structure: see the update!() method of the material.
- `inspector` - functionwith the signature
       idat = inspector(idat, j, conn, x, out, loc);
  where
   `idat` - a structure or an array that the inspector may
          use to maintain some state,  for instance minimum or maximum of
          stress, `j` is the element number, `conn` is the element connectivity,
          `out` is the output of the update!() method,  `loc` is the location
          of the integration point in the *reference* configuration.

# Return
The updated inspector data is returned.
"""
function inspectintegpoints(self::AbstractFEMMDeforLinearESNICE, geom::NodalField{FFlt},  u::NodalField{T}, dT::NodalField{FFlt}, felist::FIntVec, inspector::F,  idat, quantity=:Cauchy; context...) where {T<:Number, F<:Function}
    fes = self.integdomain.fes
    npts,  Ns,  gradNparams,  w,  pc = integrationdata(self.integdomain);
    ecoords, dofnums, loc, J, csmatTJ, Jac, D, Dstab = _buffers2(self, geom, u, npts)
    # Sort out  the output requirements
    outputcsys = self.mcsys; # default: report the stresses in the material coord system
    for apair in pairs(context)
    	sy, val = apair
    	if sy == :outputcsys
    		outputcsys = val
    	end
    end
    t= 0.0
    dt = 0.0
    qpdT = 0.0; # node temperature increment
    qpstrain = fill(zero(FFlt), nstressstrain(self.mr), 1); # total strain -- buffer
    qpthstrain = fill(zero(FFlt), nthermstrain(self.mr)); # thermal strain -- buffer
    qpstress = fill(zero(FFlt), nstressstrain(self.mr)); # stress -- buffer
    out1 = fill(zero(FFlt), nstressstrain(self.mr)); # stress -- buffer
    out =  fill(zero(FFlt), nstressstrain(self.mr));# output -- buffer
    outtot =  fill(zero(FFlt), nstressstrain(self.mr));# output -- buffer
    xe = fill(0.0, nodesperelem(fes), ndofs(geom))
    eue = fill(zero(T), nodesperelem(fes) * ndofs(u))
    gradN = fill(zero(FFlt), nodesperelem(fes), ndofs(u));
    B = fill(zero(FFlt), nstressstrain(self.mr), nodesperelem(fes) * ndofs(u));
    # Loop over  all the elements and all the quadrature points within them
    for ilist = 1:length(felist) # Loop over elements
    	i = felist[ilist];
    	gathervalues_asmat!(geom, ecoords, fes.conn[i]);# retrieve element coords
    	for nix = fes.conn[i] # For all nodes connected by this element
    		nodalgradN = self.nodalbasisfunctiongrad[nix].gradN
    		patchconn = self.nodalbasisfunctiongrad[nix].patchconn
    		Vpatch = self.nodalbasisfunctiongrad[nix].Vpatch
    		ue = fill(zero(T), length(patchconn) * ndofs(u))
    		gathervalues_asvec!(u, ue, patchconn);# retrieve element displacements
    		loc = reshape(geom.values[nix, :], 1, ndofs(geom))
    		updatecsmat!(self.mcsys, loc, J, nix);
    		nd = length(patchconn) * ndofs(u)
    		Bnodal = fill(0.0, size(D, 1), nd)
    		Blmat!(self.mr, Bnodal, Ns[1], nodalgradN, loc, self.mcsys.csmat);
    		updatecsmat!(outputcsys, loc, J, nix); # Update output coordinate system
    		# Quadrature point quantities
    		A_mul_B!(qpstrain, Bnodal, ue); # strain in material coordinates
    		qpdT = dT.values[nix] # Quadrature point temperature increment
    		thermalstrain!(self.material, qpthstrain, qpdT)
    		# Real Material: update the state and return the output
    		out = update!(self.material, qpstress, out, vec(qpstrain), qpthstrain, t, dt, loc, nix, quantity)
    		copyto!(outtot, out)
    		# Stabilization Material: update the state and return the output
    		out = update!(self.stabilization_material, qpstress, out, vec(qpstrain), qpthstrain, t, dt, loc, nix, quantity)
    		outtot .+= -self.nphis[nix].*out
    		pci = findfirst(cx -> cx == nix, fes.conn[i]);# at which node are we?
    		locjac!(loc, J, geom.values, fes.conn[i], Ns[pci], gradNparams[pci])
    		updatecsmat!(self.mcsys, loc, J, fes.label[i]);
    		At_mul_B!(csmatTJ, self.mcsys.csmat, J); # local Jacobian matrix
    		gradN!(fes, gradN, gradNparams[pci], csmatTJ);
    		Blmat!(self.mr, B, Ns[pci], gradN, loc, self.mcsys.csmat);
    		gathervalues_asvec!(u, eue, fes.conn[i]);# retrieve element displacements
    		A_mul_B!(qpstrain, B, eue); # strain in material coordinates
    		out = update!(self.stabilization_material, qpstress, out, vec(qpstrain), qpthstrain, t, dt, loc, nix, quantity)
    		outtot .+= +self.ephis[i].*out
    		out, outtot = outtot, out # swap in the total output
    		if (quantity == :Cauchy)   # Transform stress tensor,  if that is "out"
    		    (length(out1) >= length(out)) || (out1 = zeros(length(out)))
    		    rotstressvec!(self.mr, out1, out, transpose(self.mcsys.csmat))# To global coord sys
    		    rotstressvec!(self.mr, out, out1, outputcsys.csmat)# To output coord sys
    		end
    		# Call the inspector
    		idat = inspector(idat, i, fes.conn[i], ecoords, out, loc);
    	end # Loop over nodes
    end # Loop over elements
    return idat
end

function inspectintegpoints(self::AbstractFEMMDeforLinearESNICE, geom::NodalField{FFlt},  u::NodalField{T}, felist::FIntVec, inspector::F, idat, quantity=:Cauchy; context...) where {T<:Number, F<:Function}
    dT = NodalField(fill(zero(FFlt), nnodes(geom), 1)) # zero difference in temperature
    return inspectintegpoints(self, geom, u, dT, felist, inspector, idat, quantity; context...);
end

"""
    infsup_gh(self::AbstractFEMMDeforLinearESNICE, assembler::A, geom::NodalField{FFlt}, u::NodalField{T}) where {A<:AbstractSysmatAssembler, T<:Number}

Compute the matrix to produce the norm of the divergence of the displacement.

This matrix is used in the numerical infsup test (Klaus-Jurgen Bathe, The
inf-sup condition and its evaluation for mixed finite element methods,
Computers and Structures 79 (2001) 243-252.)

!!! note
This computation has not been optimized in any way. It can be expected to be
inefficient.
"""
function infsup_gh(self::AbstractFEMMDeforLinearESNICE, assembler::A, geom::NodalField{FFlt}, u::NodalField{T}) where {A<:AbstractSysmatAssembler, T<:Number}
	fes = self.integdomain.fes
	npts,  Ns,  gradNparams,  w,  pc = integrationdata(self.integdomain);
	elmatsizeguess = 4*nodesperelem(fes)*ndofs(u)
	startassembly!(assembler, elmatsizeguess, elmatsizeguess, nnodes(u) + count(fes), u.nfreedofs, u.nfreedofs);
	for nix = 1:length(self.nodalbasisfunctiongrad)
	    gradN = self.nodalbasisfunctiongrad[nix].gradN
	    patchconn = self.nodalbasisfunctiongrad[nix].patchconn
	    Vpatch = self.nodalbasisfunctiongrad[nix].Vpatch
	    c = reshape(geom.values[nix, :], 1, ndofs(geom))
	    nd = length(patchconn) * ndofs(u)
	    divm = divmat(self.mr, Ns[1], gradN, c);
	    elmat = (transpose(divm) * divm) * Vpatch
	    dofnums = fill(0, nd)
	    gatherdofnums!(u, dofnums, patchconn); # retrieve degrees of freedom
	    assemble!(assembler, elmat, dofnums, dofnums); # assemble symmetric matrix
	end # Loop over nodes
	return makematrix!(assembler);
end

function infsup_gh(self::AbstractFEMMDeforLinearESNICE, geom::NodalField{FFlt},  u::NodalField{T}) where {T<:Number}
    assembler = SysmatAssemblerSparseSymm();
    return infsup_gh(self, assembler, geom, u);
end

"""
    infsup_sh(self::AbstractFEMMDeforLinearESNICE, assembler::A, geom::NodalField{FFlt}, u::NodalField{T}) where {A<:AbstractSysmatAssembler, T<:Number}

Compute the matrix to produce the seminorm of the displacement (square root of
the sum of the squares of the derivatives of the components of displacement).

This matrix is used in the numerical infsup test (Klaus-Jurgen Bathe, The
inf-sup condition and its evaluation for mixed finite element methods,
Computers and Structures 79 (2001) 243-252.)

!!! note
This computation has not been optimized in any way. It can be expected to be
inefficient.
"""
function infsup_sh(self::AbstractFEMMDeforLinearESNICE, assembler::A, geom::NodalField{FFlt}, u::NodalField{T}) where {A<:AbstractSysmatAssembler, T<:Number}
	fes = self.integdomain.fes
	npts,  Ns,  gradNparams,  w,  pc = integrationdata(self.integdomain);
	elmatsizeguess = 4*nodesperelem(fes)*ndofs(u)
	startassembly!(assembler, elmatsizeguess, elmatsizeguess, nnodes(u) + count(fes), u.nfreedofs, u.nfreedofs);
	for nix = 1:length(self.nodalbasisfunctiongrad)
	    gradN = self.nodalbasisfunctiongrad[nix].gradN
	    patchconn = self.nodalbasisfunctiongrad[nix].patchconn
	    Vpatch = self.nodalbasisfunctiongrad[nix].Vpatch
	    c = reshape(geom.values[nix, :], 1, ndofs(geom))
	    nd = length(patchconn) * ndofs(u)
	    vgradm = vgradmat(self.mr, Ns[1], gradN, c);
	    elmat = (transpose(vgradm) * vgradm) * Vpatch
	    dofnums = fill(0, nd)
	    gatherdofnums!(u, dofnums, patchconn); # retrieve degrees of freedom
	    assemble!(assembler, elmat, dofnums, dofnums); # assemble symmetric matrix
	end # Loop over nodes
	return makematrix!(assembler);
end

function infsup_sh(self::AbstractFEMMDeforLinearESNICE, geom::NodalField{FFlt},  u::NodalField{T}) where {T<:Number}
    assembler = SysmatAssemblerSparseSymm();
    return infsup_sh(self, assembler, geom, u);
end

end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

export latex

function wrapblock(f, io, env)
    println(io, "\\begin{", env, "}")
    f()
    println(io, "\\end{", env, "}")
end

function wrapinline(f, io, cmd)
    print(io, "\\", cmd, "{")
    f()
    print(io, "}")
end

# Block elements

latex(io::IO, md::MD) = latex(io, md.content)

function latex(io::IO, content::Vector)
    for c in content
        latex(io, c)
    end
end

function latex(io::IO, header::Header{l}) where l
    tag = l < 4 ? "sub"^(l-1) * "section" : "sub"^(l-4) * "paragraph"
    wrapinline(io, tag) do
        latexinline(io, header.text)
    end
    println(io)
end

function latex(io::IO, code::Code)
    wrapblock(io, "verbatim") do
        # TODO latex escape
        println(io, code.code)
    end
end

function latexinline(io::IO, code::Code)
    wrapinline(io, "texttt") do
        print(io, latexesc(code.code))
    end
end

function latex(io::IO, md::Paragraph)
    for md in md.content
        latexinline(io, md)
    end
    println(io)
    println(io)
end

function latex(io::IO, md::BlockQuote)
    wrapblock(io, "quote") do
        latex(io, md.content)
    end
end

function latex(io::IO, f::Footnote)
    print(io, "\\footnotetext[", f.id, "]{")
    latex(io, f.text)
    println(io, "}")
end

function latex(io::IO, md::Admonition)
    wrapblock(io, "quote") do
        wrapinline(io, "textbf") do
            print(io, md.category)
        end
        println(io, "\n\n", md.title, "\n")
        latex(io, md.content)
    end
end

function latex(io::IO, md::List)
    # `\begin{itemize}` is used here for both ordered and unordered lists since providing
    # custom starting numbers for enumerated lists is simpler to do by manually assigning
    # each number to `\item` ourselves rather than using `\setcounter{enumi}{<start>}`.
    #
    # For an ordered list starting at 5 the following will be generated:
    #
    # \begin{itemize}
    #   \item[5. ] ...
    #   \item[6. ] ...
    #   ...
    # \end{itemize}
    #
    pad = ndigits(md.ordered + length(md.items)) + 2
    fmt = n -> (isordered(md) ? "[$(rpad("$(n + md.ordered - 1).", pad))]" : "")
    wrapblock(io, "itemize") do
        for (n, item) in enumerate(md.items)
            print(io, "\\item$(fmt(n)) ")
            latex(io, item)
            n < length(md.items) && println(io)
        end
    end
end

function show(io::IO, ::MIME"text/latex", md::HorizontalRule)
    println(io, "\\rule{\\textwidth}{1pt}")
end

# Inline elements

function latexinline(io::IO, md::Vector)
    for c in md
        latexinline(io, c)
    end
end

function latexinline(io::IO, md::AbstractString)
    latexesc(io, md)
end

function latexinline(io::IO, md::Bold)
    wrapinline(io, "textbf") do
        latexinline(io, md.text)
    end
end

function latexinline(io::IO, md::Italic)
    wrapinline(io, "emph") do
        latexinline(io, md.text)
    end
end

function latexinline(io::IO, md::Image)
    wrapblock(io, "figure") do
        println(io, "\\centering")
        wrapinline(io, "includegraphics") do
            print(io, md.url)
        end
        println(io)
        wrapinline(io, "caption") do
            latexinline(io, md.alt)
        end
        println(io)
    end
end

latexinline(io::IO, f::Footnote) = print(io, "\\footnotemark[", f.id, "]")

function latexinline(io::IO, md::Link)
    wrapinline(io, "href") do
        print(io, md.url)
    end
    print(io, "{")
    latexinline(io, md.text)
    print(io, "}")
end

const _latexescape_chars = Dict{Char, AbstractString}(
   '~'=>"{\\sim}", '^'=>"\\^{}", '\\'=>"{\\textbackslash}")
for ch in "&%\$#_{}"
    _latexescape_chars[ch] = "\\$ch"
end

function latexesc(io, s::AbstractString)
    for ch in s
        print(io, get(_latexescape_chars, ch, ch))
    end
end

latex(md) = sprint(latex, md)
latexinline(md) = sprint(latexinline, md)
latexesc(s) = sprint(latexesc, s)

show(io::IO, ::MIME"text/latex", md::MD) = latex(io, md)

max_stars_repo_name: caodi/julia

max_issues_repo_name: caodi/julia

max_forks_repo_name: caodi/julia


text:
using Turing

GKernel(var) = (x) -> Normal(x, sqrt.(var))

varname(s::Symbol) = nothing, s
function varname(expr::Expr)
    # Initialize an expression block to store the code for creating uid
    local sym
    @assert expr.head == :ref "expr needs to be an indexing expression, e.g. :(x[1])"
    indexing_ex = Expr(:block)
    # Add the initialization statement for uid
    push!(indexing_ex.args, quote indexing_list = [] end)
    # Initialize a local container for parsing and add the expr to it
    to_eval = []; pushfirst!(to_eval, expr)
    # Parse the expression and creating the code for creating uid
    find_head = false
    while length(to_eval) > 0
        evaling = popfirst!(to_eval)   # get the current expression to deal with
        if isa(evaling, Expr) && evaling.head == :ref && ~find_head
            # Add all the indexing arguments to the left
            pushfirst!(to_eval, "[", insdelim(evaling.args[2:end])..., "]")
            # Add first argument depending on its type
            # If it is an expression, it means it's a nested array calling
            # Otherwise it's the symbol for the calling
            if isa(evaling.args[1], Expr)
                pushfirst!(to_eval, evaling.args[1])
            else
                # push!(indexing_ex.args, quote pushfirst!(indexing_list, $(string(evaling.args[1]))) end)
                push!(indexing_ex.args, quote sym = Symbol($(string(evaling.args[1]))) end) # store symbol in runtime
                find_head = true
                sym = evaling.args[1] # store symbol in compilation time
            end
        else
            # Evaluting the concrete value of the indexing variable
            push!(indexing_ex.args, quote push!(indexing_list, string($evaling)) end)
        end
    end
    push!(indexing_ex.args, quote indexing = reduce(*, indexing_list) end)
    return indexing_ex, sym
end

genvn(sym::Symbol) = VarName(gensym(), sym, "", 1)
function genvn(expr::Expr)
    ex, sym = varname(expr)
    VarName(gensym(), sym, eval(ex), 1)
end

function randr(vi::Turing.VarInfo,
               vn::Turing.VarName,
               dist::Distribution,
               spl::Turing.Sampler,
               count::Bool = false)
    if ~haskey(vi, vn)
        r = rand(dist)
        Turing.push!(vi, vn, r, dist, spl)
        return r
    elseif is_flagged(vi, vn, "del")
        unset_flag!(vi, vn, "del")
        r = rand(dist)
        Turing.RandomVariables.setval!(vi, Turing.vectorize(dist, r), vn)
        return r
    else
        if count Turing.checkindex(vn, vi, spl) end
        Turing.updategid!(vi, vn, spl)
        return vi[vn]
    end
end

function insdelim(c, deli=",")
	return reduce((e, res) -> append!(e, [res, deli]), c; init = [])[1:end-1]
end