text:
const OUDouble = OU+OU
@testset "typestructure" begin
    @testset "AdditiveTimeSeriesModel" begin
        @test AdditiveTimeSeriesModel(OU(1,1),OU(1,1)) == AdditiveTimeSeriesModel{OU,OU,1,Float64}(ones(4))
        @test_throws MethodError AdditiveTimeSeriesModel(OU(1,1),CorrelatedOU(1,1,1/2))
        @test OUDouble == AdditiveTimeSeriesModel{OU,OU,1,Float64}
        @test_throws MethodError OU+CorrelatedOU
        @test OU(1,1)+OU(1,1) == AdditiveTimeSeriesModel{OU,OU,1,Float64}(ones(4))
    end
    @testset "ndims" begin
        @test ndims(OU) == 1
        @test ndims(Matern{10,55}) == 10
        @test ndims(CorrelatedOU(1,1,1/2)) == 2
    end
    @testset "nlowertriangle_dimension" begin
        @test nlowertriangle_dimension(OU) == 1
        @test nlowertriangle_dimension(Matern{10,55}) == 55
        @test nlowertriangle_dimension(CorrelatedOU(1,1,1/2)) == 3
    end
    @testset "npars" begin
        @test npars(OU) == 2
        @test npars(CorrelatedOU(1,1,1/2)) == 3
        @test npars(OUDouble) == 4
    end
    @testset "nlowertriangle_parameter" begin
        @test nlowertriangle_parameter(OU) == 3
        @test nlowertriangle_parameter(CorrelatedOU(1,1,1/2)) == 6
        @test nlowertriangle_parameter(OUDouble) == 10
    end
    @testset "indexLT" begin
        @test indexLT(2,3,3) == 5
        @test indexLT(4,2,5) == 8
    end
    @testset "nalias" begin
        @test nalias(OU(1,1)) isa Int
    end
    @testset "minbins" begin
        @test minbins(OU) isa Int
    end
end


text:
"""
    rand(::Type{DQMC}, m::Model, nslices::Int)

Draw random configuration.
"""
Base.rand(::Type{DQMC}, m::Model, nslices) = throw(MethodError(rand, (DQMC, m, nslices)))

"""
    nflavors(model)

Returns the number of activer fermion flavors of a given Quantum Monte Carlo
model.

The size of the hopping matrix, the interaction matrix, the greens matrix and
any derived quantity should have the
`(length(lattice(model)) * nflavors(model), length(lattice(model)) * nflavors(model)`.
"""
nflavors(::T) where {T <: Model} = error("nflavors() not implemented for $T")

"""
    hopping_matrix(mc::DQMC, m::Model)

Calculates the hopping matrix \$T_{i\\sigma, j\\sigma '}\$ where \$i, j\$ are
site indices and \$\\sigma , \\sigma '\$ are flavor indices (e.g. spin indices).
The hopping matrix should also contain potential chemical potential terms on the
diagonal.

A matrix element is the hopping amplitude for a hopping process: \$j,\\sigma '
\\rightarrow i,\\sigma\$.

Regarding the order of indices, if `T[i, σ, j, σ']` is your desired 4D hopping
array, then `reshape(T, (n_sites * n_flavors, :))` is the hopping matrix.
"""
hopping_matrix(mc::DQMC, m::Model) = throw(MethodError(hopping_matrix, (mc, m)))

"""
    interaction_matrix_exp!(
        mc::DQMC, m::Model, result::AbstractArray, conf, slice::Int, power::Float64=1.
    ) -> nothing

Calculate the interaction matrix exponential
`expV = exp(- power * delta_tau * V(slice))` and store it in `result::Matrix`.
Potential chemical potential terms should be part of the `hopping_matrix` and
not the interaction.

This is a performance critical method and one might consider efficient in-place
(in `result`) construction.
"""
function interaction_matrix_exp!(mc::DQMC, m::Model, result::AbstractArray, conf, slice, power=1.0)
    throw(MethodError(interaction_matrix_exp!, (mc, m, result, conf, slice, power)))
end

"""
    propose_local(mc::DQMC, m::Model, i::Int, conf) -> detratio, ΔE_boson, passthrough

Propose a local move for lattice site `i` of current configuration `conf` .
Returns the Green's function determinant ratio, the boson energy difference
`ΔE_boson = E_boson_new - E_boson`,
and potentially additional local move information `passthrough` (will be
forwarded to `accept_local!`).

See also [`accept_local!`](@ref).
"""
function propose_local(mc::DQMC, m::Model, i, slice, conf)
    throw(MethodError(propose_local, (mc, m, i, slice, conf)))
end

"""
    accept_local(mc::DQMC, m::Model, i::Int, slice::Int, conf, detratio, ΔE_boson, passthrough)

Accept a local move for site `i` at imaginary time slice `slice` of current
configuration `conf`. Arguments `detratio`, `ΔE_boson` and `passthrough`
correspond to output of `propose_local` for that local move.

See also [`propose_local`](@ref).
"""
function accept_local!(mc::DQMC, m::Model, i, slice, conf, detratio, ΔE_boson, passthrough)
    throw(MethodError(accept_local!, (mc, m, i, slice, conf, detratio, ΔE_boson, passthrough)))
end

ext: jl


text:
using XXZNumerics
using Test, Statistics, Random, LinearAlgebra

@testset "XXZNumerics.jl" begin
    # Write your tests here.
    include("geometry.jl")

    include("interaction.jl")

    include("spin_half.jl")

    include("entropy.jl")

    include("hamiltonian.jl")

    include("thermalization.jl")

end


text:
skipnan(x) = Iterators.filter(!isnan, x)

function get_dims(perspective::Symbol)
    dims = if perspective === :edgeon
        [1, 3]
    elseif perspective === :sideon
        [2, 3]
    elseif perspective === :faceon
        [1, 2]
    else
        error("`perspective` is $perspective, needs to be one of `:edgeon`, `:sideon`, `:faceon`")
    end
end

function LinearAlgebra.cross(a::AbstractMatrix, b::AbstractMatrix)
    @assert size(a) == size(b)
    @assert size(a, 1) == 3

    c = Array{typeof(a[1] * b[1])}(undef, 3, size(a, 2))
    @views @. c[1, :] = a[2, :] * b[3, :] - a[3, :] * b[2, :]
    @views @. c[2, :] = a[3, :] * b[1, :] - a[1, :] * b[3, :]
    @views @. c[3, :] = a[1, :] * b[2, :] - a[2, :] * b[1, :]

    return c
end


text:
"compose a markdown docstring based on a Doxygen XML element"
function build_docstring(fnode)
    io = IOBuffer()

    # code block with function definition
    print(io, "    ", text(fnode, "name"), "(")
    nspace = position(io)
    params = LibExpat.find(fnode, "param")
    if isempty(params)
        print(io, ")")
    end
    for param in params
        param !== first(params) && print(io, " "^nspace) # align
        islast = param === last(params)
        print(io, text(param, "{type}"))
        # declname not always present (with void)
        lastchar = islast ? ")" : ",\n"
        if !isempty(LibExpat.find(param, "{declname}"))
            print(io, " ", text(param, "{declname}"), lastchar)
        else
            print(io, lastchar)
        end
    end
    println(io, " -> ", text(fnode, "type"))

    # brief description
    if !isempty(LibExpat.find(fnode, "briefdescription")) &&
        !isempty(strip(LibExpat.find(fnode, "briefdescription#string")))
        println(io, "\n", text(fnode, "briefdescription"))
    end

    # parameters
    params = LibExpat.find(fnode, "detaileddescription/para/parameterlist/parameteritem")
    if !isempty(params)
        println(io, "\n### Parameters")
        for param in params
            print(io, "* **", text(param, "parameternamelist"), "**: ")
            println(io, text(param, "parameterdescription"))
        end
    end

    # returns
    return_elems = xpath(fnode, "detaileddescription/para/simplesect[@kind='return']")
    if !isempty(return_elems)
        println(io, "\n### Returns")
        println(io, text(return_elems[1], "para"))
    end

    String(take!(io))
end

# Change GENERATE_XML to YES in the Doxyfile from the GDAL SVN
# then run combine.xslt to create this XML file.
xmlfile = joinpath(workdir, "doxygen.xml")

xmlstring = readstring(xmlfile)
et = xp_parse(xmlstring)
functionnode(et, fname) = et["/doxygen/compounddef/sectiondef/memberdef[name='$fname']"]
text(node::LibExpat.ETree, el::AbstractString) = strip(LibExpat.find(node, "$el#string"))


text:
function now() return round(Int64, time() * 1000) end

function primes!(numbers)
    n = length(numbers)
    hn = n / 2
    i = 1
    while (i += 1) < hn
        if numbers[i] == 0
            continue
        end
        x = numbers[i]
        j = i
        while (j += x) <= n
            numbers[j] = 0
        end
    end
end

function benchmark(n)
    numbers = collect(1:2:n)
    ms = now()
    primes!(numbers)
    ms = now() - ms
    out = numbers[numbers .!= 0]
    println("range: ", length(numbers))
    println("primes: ", length(out))
    println("ms: ", ms, "/sec: ", ms / 1000)
end

benchmark(100000000)

max_stars_repo_name: OliverKovacs/primes

max_issues_repo_name: OliverKovacs/primes

max_forks_repo_name: OliverKovacs/primes


text:
#
# constants.jl --
#
# Constant definitions for Julia wrapper of X11 library.
#
#------------------------------------------------------------------------------
#
# Copyright (C) 2017, Éric Thiébaut.
#
# This file is part of Xlib.jl which is licensed under the MIT "Expat" License.
#

const X_PROTOCOL          = Cint(11) # current protocol version
const X_PROTOCOL_REVISION = Cint(0)  # current minor version

###############################################################################
# RESERVED RESOURCE AND CONSTANT DEFINITIONS
#------------------------------------------------------------------------------

const None            = Clong(0) # universal null resource or null atom
const ParentRelative  = Clong(1) # background pixmap in CreateWindow
const CopyFromParent  = Clong(0) # border pixmap in CreateWindow
const PointerWindow   = Clong(0) # destination window in SendEvent
const InputFocus      = Clong(1) # destination window in SendEvent
const PointerRoot     = Clong(1) # focus window in SetInputFocus
const AnyPropertyType = Clong(0) # special Atom, passed to GetProperty
const AnyKey          = Clong(0) # special Key Code, passed to GrabKey
const AnyButton       = Clong(0) # special Button Code, passed to GrabButton
const AllTemporary    = Clong(0) # special Resource ID passed to KillClient
const CurrentTime     = Clong(0) # special Time
const NoSymbol        = Clong(0) # special KeySym

###############################################################################
# EVENT DEFINITIONS
#------------------------------------------------------------------------------

# Input Event Masks. Used as event-mask window attribute and as arguments to
# Grab requests.  Not to be confused with event names.

const NoEventMask              = Clong(0)
const KeyPressMask             = Clong(1)<<0
const KeyReleaseMask           = Clong(1)<<1
const ButtonPressMask          = Clong(1)<<2
const ButtonReleaseMask        = Clong(1)<<3
const EnterWindowMask          = Clong(1)<<4
const LeaveWindowMask          = Clong(1)<<5
const PointerMotionMask        = Clong(1)<<6
const PointerMotionHintMask    = Clong(1)<<7
const Button1MotionMask        = Clong(1)<<8
const Button2MotionMask        = Clong(1)<<9
const Button3MotionMask        = Clong(1)<<10
const Button4MotionMask        = Clong(1)<<11
const Button5MotionMask        = Clong(1)<<12
const ButtonMotionMask         = Clong(1)<<13
const KeymapStateMask          = Clong(1)<<14
const ExposureMask             = Clong(1)<<15
const VisibilityChangeMask     = Clong(1)<<16
const StructureNotifyMask      = Clong(1)<<17
const ResizeRedirectMask       = Clong(1)<<18
const SubstructureNotifyMask   = Clong(1)<<19
const SubstructureRedirectMask = Clong(1)<<20
const FocusChangeMask          = Clong(1)<<21
const PropertyChangeMask       = Clong(1)<<22
const ColormapChangeMask       = Clong(1)<<23
const OwnerGrabButtonMask      = Clong(1)<<24

# Event names.  Used in "type" field in XEvent structures.  Not to be confused
# with event masks above.  They start from 2 because 0 and 1 are reserved in
# the protocol for errors and replies.

const KeyPress         = Cint( 2)
const KeyRelease       = Cint( 3)
const ButtonPress      = Cint( 4)
const ButtonRelease    = Cint( 5)
const MotionNotify     = Cint( 6)
const EnterNotify      = Cint( 7)
const LeaveNotify      = Cint( 8)
const FocusIn          = Cint( 9)
const FocusOut         = Cint(10)
const KeymapNotify     = Cint(11)
const Expose           = Cint(12)
const GraphicsExpose   = Cint(13)
const NoExpose         = Cint(14)
const VisibilityNotify = Cint(15)
const CreateNotify     = Cint(16)
const DestroyNotify    = Cint(17)
const UnmapNotify      = Cint(18)
const MapNotify        = Cint(19)
const MapRequest       = Cint(20)
const ReparentNotify   = Cint(21)
const ConfigureNotify  = Cint(22)
const ConfigureRequest = Cint(23)
const GravityNotify    = Cint(24)
const ResizeRequest    = Cint(25)
const CirculateNotify  = Cint(26)
const CirculateRequest = Cint(27)
const PropertyNotify   = Cint(28)
const SelectionClear   = Cint(29)
const SelectionRequest = Cint(30)
const SelectionNotify  = Cint(31)
const ColormapNotify   = Cint(32)
const ClientMessage    = Cint(33)
const MappingNotify    = Cint(34)
const GenericEvent     = Cint(35)
const LASTEvent        = Cint(36)     # must be bigger than any event

# Key masks. Used as modifiers to GrabButton and GrabKey, results of
# QueryPointer, state in various key-, mouse-, and button-related events.

const ShiftMask   = Cint(1)<<0
const LockMask    = Cint(1)<<1
const ControlMask = Cint(1)<<2
const Mod1Mask    = Cint(1)<<3
const Mod2Mask    = Cint(1)<<4
const Mod3Mask    = Cint(1)<<5
const Mod4Mask    = Cint(1)<<6
const Mod5Mask    = Cint(1)<<7

# Modifier names.  Used to build a SetModifierMapping request or to read a
# GetModifierMapping request.  These correspond to the masks defined above.

const ShiftMapIndex   = Cint(0)
const LockMapIndex    = Cint(1)
const ControlMapIndex = Cint(2)
const Mod1MapIndex    = Cint(3)
const Mod2MapIndex    = Cint(4)
const Mod3MapIndex    = Cint(5)
const Mod4MapIndex    = Cint(6)
const Mod5MapIndex    = Cint(7)

# Button masks.  Used in same manner as Key masks above. Not to be confused
# with button names below.

const Button1Mask = Cint(1)<<8
const Button2Mask = Cint(1)<<9
const Button3Mask = Cint(1)<<10
const Button4Mask = Cint(1)<<11
const Button5Mask = Cint(1)<<12

const AnyModifier = Cint(1)<<15  # used in GrabButton, GrabK

# Button names. Used as arguments to GrabButton and as detail in ButtonPress
# and ButtonRelease events.  Not to be confused with button masks above.  Note
# that 0 is already defined above as "AnyButton".

const Button1 = Cint(1)
const Button2 = Cint(2)
const Button3 = Cint(3)
const Button4 = Cint(4)
const Button5 = Cint(5)

# Notify modes

const NotifyNormal       = Cint(0)
const NotifyGrab         = Cint(1)
const NotifyUngrab       = Cint(2)
const NotifyWhileGrabbed = Cint(3)

const NotifyHint         = Cint(1) # for MotionNotify events

# Notify detail

const NotifyAncestor         = Cint(0)
const NotifyVirtual          = Cint(1)
const NotifyInferior         = Cint(2)
const NotifyNonlinear        = Cint(3)
const NotifyNonlinearVirtual = Cint(4)
const NotifyPointer          = Cint(5)
const NotifyPointerRoot      = Cint(6)
const NotifyDetailNone       = Cint(7)

# Visibility notify

const VisibilityUnobscured        = Cint(0)
const VisibilityPartiallyObscured = Cint(1)
const VisibilityFullyObscured     = Cint(2)

# Circulation request

const PlaceOnTop    = Cint(0)
const PlaceOnBottom = Cint(1)

# protocol families

const FamilyInternet  = Cint(0) # IPv4
const FamilyDECnet    = Cint(1)
const FamilyChaos     = Cint(2)
const FamilyInternet6 = Cint(6) # IPv6

# authentication families not tied to a specific protocol

const FamilyServerInterpreted = Cint(5)

# Property notification

const PropertyNewValue  = Cint(0)
const PropertyDelete    = Cint(1)

# Color Map notification

const ColormapUninstalled = Cint(0)
const ColormapInstalled   = Cint(1)

# GrabPointer, GrabButton, GrabKeyboard, GrabKey Modes

const GrabModeSync  = Cint(0)
const GrabModeAsync = Cint(1)

# GrabPointer, GrabKeyboard reply status

const GrabSuccess     = Cint(0)
const AlreadyGrabbed  = Cint(1)
const GrabInvalidTime = Cint(2)
const GrabNotViewable = Cint(3)
const GrabFrozen      = Cint(4)

# AllowEvents modes

const AsyncPointer   = Cint(0)
const SyncPointer    = Cint(1)
const ReplayPointer  = Cint(2)
const AsyncKeyboard  = Cint(3)
const SyncKeyboard   = Cint(4)
const ReplayKeyboard = Cint(5)
const AsyncBoth      = Cint(6)
const SyncBoth       = Cint(7)

# Used in SetInputFocus, GetInputFocus

const RevertToNone        = Cint(None)
const RevertToPointerRoot = Cint(PointerRoot)
const RevertToParent      = Cint(2)

###############################################################################
# ERROR CODES
#------------------------------------------------------------------------------

const Success           = Cint(0)  # everything's okay
const BadRequest        = Cint(1)  # bad request code
const BadValue          = Cint(2)  # int parameter out of range
const BadWindow         = Cint(3)  # parameter not a Window
const BadPixmap         = Cint(4)  # parameter not a Pixmap
const BadAtom           = Cint(5)  # parameter not an Atom
const BadCursor         = Cint(6)  # parameter not a Cursor
const BadFont           = Cint(7)  # parameter not a Font
const BadMatch          = Cint(8)  # parameter mismatch
const BadDrawable       = Cint(9)  # parameter not a Pixmap or Window
const BadAccess         = Cint(10) # depending on context:
                                   # - key/button already grabbed
                                   # - attempt to free an illegal
                                   #   cmap entry
                                   # - attempt to store into a read-only
                                   #   color map entry.
                                   # - attempt to modify the access control
                                   #   list from other than the local host.
const BadAlloc          = Cint(11) # insufficient resources
const BadColor          = Cint(12) # no such colormap
const BadGC             = Cint(13) # parameter not a GC
const BadIDChoice       = Cint(14) # choice not in range or already used
const BadName           = Cint(15) # font or color name doesn't exist
const BadLength         = Cint(16) # Request length incorrect
const BadImplementation = Cint(17) # server is defective

const FirstExtensionError = Cint(128)
const LastExtensionError  = Cint(255)

###############################################################################
# WINDOW DEFINITIONS
#------------------------------------------------------------------------------

# Window classes used by CreateWindow
# Note that CopyFromParent is already defined as 0 above

const InputOutput = Cint(1)
const InputOnly   = Cint(2)

# Window attributes for CreateWindow and ChangeWindowAttributes

const CWBackPixmap       = Clong(1)<<0
const CWBackPixel        = Clong(1)<<1
const CWBorderPixmap     = Clong(1)<<2
const CWBorderPixel      = Clong(1)<<3
const CWBitGravity       = Clong(1)<<4
const CWWinGravity       = Clong(1)<<5
const CWBackingStore     = Clong(1)<<6
const CWBackingPlanes    = Clong(1)<<7
const CWBackingPixel     = Clong(1)<<8
const CWOverrideRedirect = Clong(1)<<9
const CWSaveUnder        = Clong(1)<<10
const CWEventMask        = Clong(1)<<11
const CWDontPropagate    = Clong(1)<<12
const CWColormap         = Clong(1)<<13
const CWCursor           = Clong(1)<<14

# ConfigureWindow structure

const CWX           = Cint(1)<<0
const CWY           = Cint(1)<<1
const CWWidth       = Cint(1)<<2
const CWHeight      = Cint(1)<<3
const CWBorderWidth = Cint(1)<<4
const CWSibling     = Cint(1)<<5
const CWStackMode   = Cint(1)<<6

# Bit Gravity

const ForgetGravity    = Cint(0)
const NorthWestGravity = Cint(1)
const NorthGravity     = Cint(2)
const NorthEastGravity = Cint(3)
const WestGravity      = Cint(4)
const CenterGravity    = Cint(5)
const EastGravity      = Cint(6)
const SouthWestGravity = Cint(7)
const SouthGravity     = Cint(8)
const SouthEastGravity = Cint(9)
const StaticGravity    = Cint(10)

# Window gravity + bit gravity above

const UnmapGravity = Cint(0)

# Used in CreateWindow for backing-store hint

const NotUseful  = Cint(0)
const WhenMapped = Cint(1)
const Always     = Cint(2)

# Used in GetWindowAttributes reply

const IsUnmapped   = Cint(0)
const IsUnviewable = Cint(1)
const IsViewable   = Cint(2)

# Used in ChangeSaveSet

const SetModeInsert = Cint(0)
const SetModeDelete = Cint(1)

# Used in ChangeCloseDownMode

const DestroyAll      = Cint(0)
const RetainPermanent = Cint(1)
const RetainTemporary = Cint(2)

# Window stacking method (in configureWindow)

const Above    = Cint(0)
const Below    = Cint(1)
const TopIf    = Cint(2)
const BottomIf = Cint(3)
const Opposite = Cint(4)

# Circulation direction

const RaiseLowest  = Cint(0)
const LowerHighest = Cint(1)

# Property modes

const PropModeReplace = Cint(0)
const PropModePrepend = Cint(1)
const PropModeAppend  = Cint(2)

###############################################################################
# GRAPHICS DEFINITIONS
#------------------------------------------------------------------------------

# Graphics functions, as in GC.alu

const GXclear        = Cint(0x0) # 0
const GXand          = Cint(0x1) # src AND dst
const GXandReverse   = Cint(0x2) # src AND NOT dst
const GXcopy         = Cint(0x3) # src
const GXandInverted  = Cint(0x4) # NOT src AND dst
const GXnoop         = Cint(0x5) # dst
const GXxor          = Cint(0x6) # src XOR dst
const GXor           = Cint(0x7) # src OR dst
const GXnor          = Cint(0x8) # NOT src AND NOT dst
const GXequiv        = Cint(0x9) # NOT src XOR dst
const GXinvert       = Cint(0xa) # NOT dst
const GXorReverse    = Cint(0xb) # src OR NOT dst
const GXcopyInverted = Cint(0xc) # NOT src
const GXorInverted   = Cint(0xd) # NOT src OR dst
const GXnand         = Cint(0xe) # NOT src OR NOT dst
const GXset          = Cint(0xf) # 1

# LineStyle

const LineSolid      = Cint(0)
const LineOnOffDash  = Cint(1)
const LineDoubleDash = Cint(2)

# capStyle

const CapNotLast    = Cint(0)
const CapButt       = Cint(1)
const CapRound      = Cint(2)
const CapProjecting = Cint(3)

# joinStyle

const JoinMiter = Cint(0)
const JoinRound = Cint(1)
const JoinBevel = Cint(2)

# fillStyle

const FillSolid          = Cint(0)
const FillTiled          = Cint(1)
const FillStippled       = Cint(2)
const FillOpaqueStippled = Cint(3)

# fillRule

const EvenOddRule = Cint(0)
const WindingRule = Cint(1)

# subwindow mode

const ClipByChildren   = Cint(0)
const IncludeInferiors = Cint(1)

# SetClipRectangles ordering

const Unsorted = Cint(0)
const YSorted  = Cint(1)
const YXSorted = Cint(2)
const YXBanded = Cint(3)

# CoordinateMode for drawing routines

const CoordModeOrigin   = Cint(0) # relative to the origin
const CoordModePrevious = Cint(1) # relative to previous point

# Polygon shapes

const _Complex  = Cint(0) # paths may intersect
const Nonconvex = Cint(1) # no paths intersect, but not convex
const Convex    = Cint(2) # wholly convex

# Arc modes for PolyFillArc

const ArcChord    = Cint(0) # join endpoints of arc
const ArcPieSlice = Cint(1) # join endpoints to center of arc

# GC components: masks used in CreateGC, CopyGC, ChangeGC, OR'ed into
# GC.stateChanges

const GCFunction          = Clong(1)<<0
const GCPlaneMask         = Clong(1)<<1
const GCForeground        = Clong(1)<<2
const GCBackground        = Clong(1)<<3
const GCLineWidth         = Clong(1)<<4
const GCLineStyle         = Clong(1)<<5
const GCCapStyle          = Clong(1)<<6
const GCJoinStyle         = Clong(1)<<7
const GCFillStyle         = Clong(1)<<8
const GCFillRule          = Clong(1)<<9
const GCTile              = Clong(1)<<10
const GCStipple           = Clong(1)<<11
const GCTileStipXOrigin   = Clong(1)<<12
const GCTileStipYOrigin   = Clong(1)<<13
const GCFont              = Clong(1)<<14
const GCSubwindowMode     = Clong(1)<<15
const GCGraphicsExposures = Clong(1)<<16
const GCClipXOrigin       = Clong(1)<<17
const GCClipYOrigin       = Clong(1)<<18
const GCClipMask          = Clong(1)<<19
const GCDashOffset        = Clong(1)<<20
const GCDashList          = Clong(1)<<21
const GCArcMode           = Clong(1)<<22

const GCLastBit           = Cint(22)

###############################################################################
# FONTS
#------------------------------------------------------------------------------

# used in QueryFont -- draw direction

const FontLeftToRight = Cint(0)
const FontRightToLeft = Cint(1)

const FontChange      = Cint(255)

###############################################################################
#  IMAGING
#------------------------------------------------------------------------------

# ImageFormat -- PutImage, GetImage

const XYBitmap = Cint(0) # depth 1, XYFormat
const XYPixmap = Cint(1) # depth == drawable depth
const ZPixmap  = Cint(2) # depth == drawable depth

###############################################################################
#  COLOR MAP STUFF
#------------------------------------------------------------------------------

# For CreateColormap

const AllocNone = Cint(0)      # create map with no entries
const AllocAll  = Cint(1)      # allocate entire map writeable

# Flags used in StoreNamedColor, StoreColors

const DoRed   = Cint(1)<<0
const DoGreen = Cint(1)<<1
const DoBlue  = Cint(1)<<2

###############################################################################
# CURSOR STUFF
#------------------------------------------------------------------------------

# QueryBestSize Class

const CursorShape  = Cint(0) # largest size that can be displayed
const TileShape    = Cint(1) # size tiled fastest
const StippleShape = Cint(2) # size stippled fastest

###############################################################################
# KEYBOARD/POINTER STUFF
#------------------------------------------------------------------------------

const AutoRepeatModeOff     = Cint(0)
const AutoRepeatModeOn      = Cint(1)
const AutoRepeatModeDefault = Cint(2)

const LedModeOff            = Cint(0)
const LedModeOn             = Cint(1)

# masks for ChangeKeyboardControl

const KBKeyClickPercent     = Clong(1)<<0
const KBBellPercent         = Clong(1)<<1
const KBBellPitch           = Clong(1)<<2
const KBBellDuration        = Clong(1)<<3
const KBLed                 = Clong(1)<<4
const KBLedMode             = Clong(1)<<5
const KBKey                 = Clong(1)<<6
const KBAutoRepeatMode      = Clong(1)<<7

const MappingSuccess        = Cint(0)
const MappingBusy           = Cint(1)
const MappingFailed         = Cint(2)

const MappingModifier       = Cint(0)
const MappingKeyboard       = Cint(1)
const MappingPointer        = Cint(2)

###############################################################################
# SCREEN SAVER STUFF
#------------------------------------------------------------------------------

const DontPreferBlanking    = Cint(0)
const PreferBlanking        = Cint(1)
const DefaultBlanking       = Cint(2)

const DisableScreenSaver    = Cint(0)
const DisableScreenInterval = Cint(0)

const DontAllowExposures    = Cint(0)
const AllowExposures        = Cint(1)
const DefaultExposures      = Cint(2)

# for ForceScreenSaver

const ScreenSaverReset      = Cint(0)
const ScreenSaverActive     = Cint(1)

###############################################################################
# HOSTS AND CONNECTIONS
#------------------------------------------------------------------------------

# for ChangeHosts

const HostInsert = Cint(0)
const HostDelete = Cint(1)

# for ChangeAccessControl

const EnableAccess  = Cint(1)
const DisableAccess = Cint(0)

# Display classes used in opening the connection
# Note that the statically allocated ones are even numbered and the
# dynamically changeable ones are odd numbered

const StaticGray  = Cint(0)
const GrayScale   = Cint(1)
const StaticColor = Cint(2)
const PseudoColor = Cint(3)
const TrueColor   = Cint(4)
const DirectColor = Cint(5)

# Byte order  used in imageByteOrder and bitmapBitOrder

const LSBFirst = Cint(0)
const MSBFirst = Cint(1)

###############################################################################
# Constants from `Xlib.h`

const XlibSpecificationRelease        = Cint(6)
const X_HAVE_UTF8_STRING              = Cint(1)

const True                            = Cint(1)
const False                           = Cint(0)

const QueuedAlready                   = Cint(0)
const QueuedAfterReading              = Cint(1)
const QueuedAfterFlush                = Cint(2)

const AllPlanes                       = ~zero(Culong)

const XNRequiredCharSet               = "requiredCharSet"
const XNQueryOrientation              = "queryOrientation"
const XNBaseFontName                  = "baseFontName"
const XNOMAutomatic                   = "omAutomatic"
const XNMissingCharSet                = "missingCharSet"
const XNDefaultString                 = "defaultString"
const XNOrientation                   = "orientation"
const XNDirectionalDependentDrawing   = "directionalDependentDrawing"
const XNContextualDrawing             = "contextualDrawing"
const XNFontInfo                      = "fontInfo"

const XIMPreeditArea                  = Clong(0x0001)
const XIMPreeditCallbacks             = Clong(0x0002)
const XIMPreeditPosition              = Clong(0x0004)
const XIMPreeditNothing               = Clong(0x0008)
const XIMPreeditNone                  = Clong(0x0010)
const XIMStatusArea                   = Clong(0x0100)
const XIMStatusCallbacks              = Clong(0x0200)
const XIMStatusNothing                = Clong(0x0400)
const XIMStatusNone                   = Clong(0x0800)

const XNVaNestedList                  = "XNVaNestedList"
const XNQueryInputStyle               = "queryInputStyle"
const XNClientWindow                  = "clientWindow"
const XNInputStyle                    = "inputStyle"
const XNFocusWindow                   = "focusWindow"
const XNResourceName                  = "resourceName"
const XNResourceClass                 = "resourceClass"
const XNGeometryCallback              = "geometryCallback"
const XNDestroyCallback               = "destroyCallback"
const XNFilterEvents                  = "filterEvents"
const XNPreeditStartCallback          = "preeditStartCallback"
const XNPreeditDoneCallback           = "preeditDoneCallback"
const XNPreeditDrawCallback           = "preeditDrawCallback"
const XNPreeditCaretCallback          = "preeditCaretCallback"
const XNPreeditStateNotifyCallback    = "preeditStateNotifyCallback"
const XNPreeditAttributes             = "preeditAttributes"
const XNStatusStartCallback           = "statusStartCallback"
const XNStatusDoneCallback            = "statusDoneCallback"
const XNStatusDrawCallback            = "statusDrawCallback"
const XNStatusAttributes              = "statusAttributes"
const XNArea                          = "area"
const XNAreaNeeded                    = "areaNeeded"
const XNSpotLocation                  = "spotLocation"
const XNColormap                      = "colorMap"
const XNStdColormap                   = "stdColorMap"
const XNForeground                    = "foreground"
const XNBackground                    = "background"
const XNBackgroundPixmap              = "backgroundPixmap"
const XNFontSet                       = "fontSet"
const XNLineSpace                     = "lineSpace"
const XNCursor                        = "cursor"
const XNQueryIMValuesList             = "queryIMValuesList"
const XNQueryICValuesList             = "queryICValuesList"
const XNVisiblePosition               = "visiblePosition"
const XNR6PreeditCallback             = "r6PreeditCallback"
const XNStringConversionCallback      = "stringConversionCallback"
const XNStringConversion              = "stringConversion"
const XNResetState                    = "resetState"
const XNHotKey                        = "hotKey"
const XNHotKeyState                   = "hotKeyState"
const XNPreeditState                  = "preeditState"
const XNSeparatorofNestedList         = "separatorofNestedList"

const XBufferOverflow                 = Cint(-1)
const XLookupNone                     = Cint(1)
const XLookupChars                    = Cint(2)
const XLookupKeySym                   = Cint(3)
const XLookupBoth                     = Cint(4)

const XIMReverse                      = Clong(1)
const XIMUnderline                    = Clong(1)<<1
const XIMHighlight                    = Clong(1)<<2
const XIMPrimary                      = Clong(1)<<5
const XIMSecondary                    = Clong(1)<<6
const XIMTertiary                     = Clong(1)<<7
const XIMVisibleToForward             = Clong(1)<<8
const XIMVisibleToBackword            = Clong(1)<<9
const XIMVisibleToCenter              = Clong(1)<<10
const XIMPreeditUnKnown               = Clong(0)
const XIMPreeditEnable                = Clong(1)
const XIMPreeditDisable               = Clong(1)<<1
const XIMInitialState                 = Clong(1)
const XIMPreserveState                = Clong(1)<<1

const XIMStringConversionLeftEdge     = Cint(0x00000001)
const XIMStringConversionRightEdge    = Cint(0x00000002)
const XIMStringConversionTopEdge      = Cint(0x00000004)
const XIMStringConversionBottomEdge   = Cint(0x00000008)
const XIMStringConversionConcealed    = Cint(0x00000010)
const XIMStringConversionWrapped      = Cint(0x00000020)
const XIMStringConversionBuffer       = Cint(0x0001)
const XIMStringConversionLine         = Cint(0x0002)
const XIMStringConversionWord         = Cint(0x0003)
const XIMStringConversionChar         = Cint(0x0004)
const XIMStringConversionSubstitution = Cint(0x0001)
const XIMStringConversionRetrieval    = Cint(0x0002)

const XIMHotKeyStateON                = Clong(0x0001)
const XIMHotKeyStateOFF               = Clong(0x0002)


text:
mutable struct raw_t <: LCMType
    utime::Int64
    length::Int32
    data::Vector{UInt8}
end

@lcmtypesetup(raw_t,
    data => (length,)
)


text:
#   Author: V. Vitale
#   Feb 2022

include("abstractTN.jl")

# MPS A-matrix is a 3-index tensor, A[i,s,j]
#    s
#    |
# i -A- j
#
# [s] acts on the local Hilbert space
# [i,j] act on the virtual bonds

mutable struct MPS <: AbstractTN
  data::Dict
  N::Int
end

MPS() = MPS(Dict(), 0)

function MPS_dot(A::MPS,B::MPS)
    E=ones(1,1)
    for i in 1:A.N
        @tensor temp[:] := E[-1,1]*A.data[i][1,-2,-3]
        #@tensor E[:] := temp[1,2,-1] * conj( B.data[i][1,2,-2] )
        @tensor E[:] := temp[1,2,-2] * conj( B.data[i][1,2,-1] )
    end
    return E[1]
end

Base.:*(A::MPS,B::MPS)=MPS_dot(A,B)

function Normalize!(A::MPS)
    norm = MPS_dot(A,A)
    for i in 1:A.N
        A.data[i] = A.data[i]/norm^(1/(2*A.N))
    end
end

function Initialize!(A::MPS,d::Int,chi::Int,N::Int)
    temp=Dict()
    temp[1] = im*rand(1,d,chi)
    for i in 2:N-1
        temp[i]= im*rand(chi,d,chi)
    end
    temp[N] = im*rand(chi,d,1)
    A.N=N
    A.data=temp
end

function Initialize!(s::String,A::MPS,N::Int)
    if s=="Brydges_Closed"
        chi=1
        d=2
        temp=Dict()
        temp[1] = im*zeros(1,d,chi)
        for i in 2:N-1
            temp[i]= im*zeros(chi,d,chi)
        end
        temp[N] = im*zeros(chi,d,1)

        temp[1][1,:,1] = [0 1]
        for i in 2:2:N-1
            temp[i][1,:,1]= [1 0]
            temp[i+1][1,:,1]= [0 1]
        end
        temp[N][1,:,1] = [1 0]
        A.N=N
        A.data=temp
    elseif s=="Brydges_Open"
        chi=1
        d=4
        p=0.008326
        temp=Dict()
        temp[1] = im*zeros(1,d,chi)
        for i in 2:N-1
            temp[i]= im*zeros(chi,d,chi)
        end
        temp[N] = im*zeros(chi,d,1)

        temp[1][1,:,1] = reshape([p 0 ; 0 1-p],d)
        for i in 2:2:N-1
            temp[i][1,:,1]= reshape([ 1 0 ; 0 0],d)
            temp[i+1][1,:,1]= reshape([p 0 ; 0 1-p],d)
        end
        temp[N][1,:,1] = reshape([ 1 0 ; 0 0],d)
        A.N=N
        A.data=temp
    end
end

function right_normalize!(A::MPS)
    for i in A.N:-1:1
        sA = size(A.data[i])
        U, S, V = svd(reshape(A.data[i],sA[1], sA[2]*sA[3]), full=false)
        V=V'
        S /= norm(S)
        A.data[i] = reshape(V,(:, sA[2], sA[3]))
        if i>1
            S=diagm(S)
            @tensor A.data[i-1][:] := A.data[i-1][-1,-2,2] * U[ 2,3 ] * S[ 3,-3 ]
        end
    end
end

function right_orthogonalize!(A::MPS)
    for i in A.N:-1:1
        sA = size(A.data[i])
        U, S, V = svd(reshape(A.data[i],sA[1], sA[2]*sA[3]), full=false)
        V=V'
        A.data[i] = reshape(V,(:, sA[2], sA[3]))
        if i>1
            S=diagm(S)
            @tensor A.data[i-1][:] := A.data[i-1][-1,-2,2] * U[ 2,3 ] * S[ 3,-3 ]
        end
    end
end

function left_normalize!(A::MPS)
    for i in 1:A.N
        sA = size(A.data[i])
        U,S,V = svd(reshape(A.data[i],(sA[1]*sA[2],sA[3])),full=false)
        S /= norm(S)
        V=V'
        A.data[i] = reshape( U,( sA[1], sA[2], :))
        if i<A.N
            S=diagm(S)
            @tensor A.data[i+1][:] := S[-1,1 ] * V[ 1,2 ] * A.data[i+1][2,-2,-3]
        end
    end
end

function left_orthogonalize!(A::MPS)
    for i in 1:A.N
        sA = size(A.data[i])
        U,S,V = svd(reshape(A.data[i],(sA[1]*sA[2],sA[3])),full=false)
        V=V'
        A.data[i] = reshape( U,( sA[1], sA[2], :))
        if i<A.N
            S=diagm(S)
            @tensor A.data[i+1][:] := S[-1,1 ] * V[ 1,2 ] * A.data[i+1][2,-2,-3]
        end
    end
end

function move_orthogonality_center!(A::MPS,b::Int)
    right_orthogonalize!(A)
    for i in 1:b
        sA = size(A.data[i])
        U,S,V = svd(reshape(A.data[i],(sA[1]*sA[2],sA[3])),full=false)
        #S /= norm(S)
        V=V'
        A.data[i] = reshape( U,( sA[1], sA[2], :))
        if i<A.N
            S=diagm(S)
            @tensor A.data[i+1][:] := S[-1,1 ] * V[ 1,2 ] * A.data[i+1][2,-2,-3]
        end
    end
end

function calc_entropy(A::MPS)
    M=copy(A)
    Sent = zeros(M.N)
    for i in 1:M.N
        sM = size(M.data[i])
        U,S,V = svd(reshape(M.data[i],(sM[1]*sM[2],sM[3])),full=false)
        #S /= norm(S)
        V=V'
        M.data[i] = reshape( U,( sM[1], sM[2], :))
        if i<M.N
            @tensor M.data[i+1][:] := diagm(S)[-1,1 ] * V[ 1,2 ] * M.data[i+1][2,-2,-3]
        end
        Sent[i] = sum(-dot(S.^2,log.(S.^2)))
    end

    return Sent
end

function calc_Renyi2(A::MPS)
    M=copy(A)
    Sent = zeros(M.N)
    for i in 1:M.N
        sM = size(M.data[i])
        U,S,V = svd(reshape(M.data[i],(sM[1]*sM[2],sM[3])),full=false)
        #S /= norm(S)
        V=V'
        M.data[i] = reshape( U,( sM[1], sM[2], :))
        if i<M.N
            @tensor M.data[i+1][:] := diagm(S)[-1,1 ] * V[ 1,2 ] * M.data[i+1][2,-2,-3]
        end
        Sent[i] = -log(sum(S.^4))
    end

    return Sent
end

function calc_purity(A::MPS)
    M=copy(A)
    Sent = zeros(M.N)
    for i in 1:M.N
        sM = size(M.data[i])
        U,S,V = svd(reshape(M.data[i],(sM[1]*sM[2],sM[3])),full=false)
        #S /= norm(S)
        V=V'
        M.data[i] = reshape( U,( sM[1], sM[2], :))
        if i<M.N
            @tensor M.data[i+1][:] := diagm(S)[-1,1 ] * V[ 1,2 ] * M.data[i+1][2,-2,-3]
        end
        Sent[i] = sum(S.^4)
    end

    return Sent
end

function calc_trace(A::MPS)
    M=copy(A)
    Sent = zeros(M.N)
    for i in 1:M.N
        sM = size(M.data[i])
        U,S,V = svd(reshape(M.data[i],(sM[1]*sM[2],sM[3])),full=false)
        #S /= norm(S)
        V=V'
        M.data[i] = reshape( U,( sM[1], sM[2], :))
        if i<M.N
            @tensor M.data[i+1][:] := diagm(S)[-1,1 ] * V[ 1,2 ] * M.data[i+1][2,-2,-3]
        end
        Sent[i] = sum(S.^2)
    end

    return Sent
end

function to_dm_MPS(A::MPS)
    M=MPS()
    M.N=A.N
    for i in 1:M.N
        sA=size(A.data[i])
        @tensor M.data[i][:] := A.data[i][-1,-3,-5]*conj(A.data[i][-2,-4,-6])#*conj(A.data[i][-2,-4,-6])
        M.data[i]= reshape(M.data[i],(sA[1]*sA[1],sA[2]*sA[2],sA[3]*sA[3]))
    end
    return M
end

"
function rdm_list(A::MPS,r::Array)
    #move_orthogonality_center!(A,1)
    sA=size(A.data[1])
    @tensor rd_rho[:] := A.data[1][-1,-3,-2]
    rd_rho=reshape(rd_rho,(sA[1],sA[3],isqrt(sA[2]),isqrt(sA[2])))
    for j in 2:A.N
        sA=size(A.data[j])
        @tensor M[:] := A.data[j][-1,-3,-2]
        M=reshape(M,(sA[1],sA[3],isqrt(sA[2]),isqrt(sA[2])))
        if j ∈ r
            st=size(rd_rho)
            sM=size(M)
            @tensor rd_rho[:] :=rd_rho[-1,1,-3,-5]*M[1,-2,-4,-6]
            rd_rho=reshape(rd_rho,(st[1],sM[2],st[3]*sM[3],st[4]*sM[4]))
        elseif j ∉ r
            @tensor rd_rho[:] :=rd_rho[-1,1,-3,-4]*M[1,-2,2,2]
        end
    end
    @tensor rd_rho[:] :=rd_rho[1,1,-1,-2]
    return rd_rho
end
"

function rdm_rho(A::MPS,r::Array)
    #move_orthogonality_center!(A,1)
    sA=size(A.data[r[1]])
    rd_rho=reshape(A.data[r[1]],(sA[1],isqrt(sA[2]),isqrt(sA[2]),sA[3]))
    for j in r[2:end]
        sA=size(A.data[j])
        M=reshape(A.data[j],(sA[1],isqrt(sA[2]),isqrt(sA[2]),sA[3]))
        st=size(rd_rho)
        sM=size(M)
        @tensor rd_rho[:] :=rd_rho[-1,-2,-4,1]*M[1,-3,-5,-6]
        rd_rho=reshape(rd_rho,(st[1],st[2]*sM[2],st[3]*sM[3],sM[4]))
    end
    for j in r[1]-1:-1:1
        sA=size(A.data[j])
        M=reshape(A.data[j],(sA[1],isqrt(sA[2]),isqrt(sA[2]),sA[3]))
        st=size(rd_rho)
        sM=size(M)
        @tensor rd_rho[:] :=M[-1,2,2,1]*rd_rho[1,-2,-3,-4]
    end
    for j in r[end]+1:A.N
        sA=size(A.data[j])
        M=reshape(A.data[j],(sA[1],isqrt(sA[2]),isqrt(sA[2]),sA[3]))
        st=size(rd_rho)
        sM=size(M)
        @tensor rd_rho[:] :=rd_rho[-1,-2,-3,1]*M[1,2,2,-4]
    end
    @tensor rd_rho[:] :=rd_rho[1,-1,-2,1]
    return rd_rho
end

function rdm_from_state(A::MPS,r::Array)
    L=ones(1,1)
    for j in 1:r[1]-1
        @tensor L[:] :=L[1,2]*A.data[j][1,3,-3]*conj(A.data[j])[2,3,-4]
    end
    R=ones(1,1)
    for j in A.N:-1:r[end]+1
        @tensor R[:] :=A.data[j][-1,3,1]*conj(A.data[j])[-2,3,2]*R[1,2]
    end

    if length(r)>1
        @tensor L[:] := L[1,2]*A.data[r[1]][1,-1,-3]*conj(A.data[r[1]])[2,-2,-4]
        @tensor R[:] := A.data[r[end]][-1,-3,1]*conj(A.data[r[end]])[-2,-4,2]*R[1,2];
        j=1
        while r[1]+j<r[end]-j
            sA=size(A.data[r[1]+j])
            sL=size(L)
            @tensor L[:] := L[-1,-3,1,2]*A.data[r[1]+j][1,-2,-5]*conj(A.data[r[1]+j])[2,-4,-6]
            L=reshape(L,(sL[1]*sA[2],sL[2]*sA[2],sA[3],sA[3]))
            sA=size(A.data[r[end]-j])
            sR=size(R)
            @tensor R[:] := A.data[r[end]-j][-1,-3,1]*conj(A.data[r[end]-j])[-2,-5,2]*R[1,2,-4,-6]
            R=reshape(R,(sA[1],sA[1],sR[3]*sA[2],sR[4]*sA[2]))
            j+=1
        end
        if mod(length(r),2)==0
            @tensor rd_rho[:] := L[-1,-3,1,2]*R[1,2,-2,-4]
        else
            sA=size(A.data[r[1]+j])
            sL=size(L)
            @tensor L[:] := L[-1,-3,1,2]*A.data[r[1]+j][1,-2,-5]*conj(A.data[r[1]+j])[2,-4,-6]
            L=reshape(L,(sL[1]*sA[2],sL[2]*sA[2],sA[3],sA[3]))
            @tensor rd_rho[:] := L[-1,-3,1,2]*R[1,2,-2,-4]
        end
        return reshape(rd_rho,(2^length(r),2^length(r)))
    else
        @tensor rd_rho[:] := L[1,2]*A.data[r[1]][1,-1,3]*conj(A.data[r[1]])[2,-2,4]*R[3,4]
        return reshape(rd_rho,(2,2))
    end
end

function rdm_from_dm(A::MPS,r::Array)
    L=ones(1)
    for j in 1:r[1]-1
        sA=size(A.data[j])
        M=reshape(A.data[j],(sA[1],isqrt(sA[2]),isqrt(sA[2]),sA[3]))
        @tensor L[:] :=L[1]*M[1,2,2,-1]
    end
    R=ones(1)
    for j in A.N:-1:r[end]+1
        sA=size(A.data[j])
        M=reshape(A.data[j],(sA[1],isqrt(sA[2]),isqrt(sA[2]),sA[3]))
        @tensor R[:] :=M[-1,2,2,1]*R[1]
    end

    if length(r)>1
        sA=size(A.data[r[1]])
        M=reshape(A.data[r[1]],(sA[1],isqrt(sA[2]),isqrt(sA[2]),sA[3]))

        @tensor L[:] := L[1]*M[1,-1,-2,-3]
        sA=size(A.data[r[end]])
        M=reshape(A.data[r[end]],(sA[1],isqrt(sA[2]),isqrt(sA[2]),sA[3]))
        @tensor R[:] := M[-1,-2,-3,1]*R[1];
        j=1
        while r[1]+j<r[end]-j
            sA=size(A.data[r[1]+j])
            M=reshape(A.data[r[1]+j],(sA[1],isqrt(sA[2]),isqrt(sA[2]),sA[3]))
            sL=size(L)
            @tensor L[:] := L[-1,-3,1]*M[1,-2,-4,-5]
            L=reshape(L,(sL[1]*isqrt(sA[2]),sL[2]*isqrt(sA[2]),sL[3]))

            sA=size(A.data[r[end]-j])
            M=reshape(A.data[r[end]-j],(sA[1],isqrt(sA[2]),isqrt(sA[2]),sA[3]))
            sR=size(R)

            @tensor R[:] := M[-1,-3,-5,1]*R[1,-2,-4]
            R=reshape(R,(sA[1],sR[2]*isqrt(sA[2]),sR[3]*isqrt(sA[2])))
            j+=1
        end
        if mod(length(r),2)==0
            @tensor rd_rho[:] := L[-1,-3,1]*R[1,-2,-4]
        else
            sA=size(A.data[r[1]+j])
            M=reshape(A.data[r[1]+j],(sA[1],isqrt(sA[2]),isqrt(sA[2]),sA[3]))
            sL=size(L)
            @tensor L[:] := L[-1,-3,1]*M[1,-2,-4,-5]
            L=reshape(L,(sL[1]*isqrt(sA[2]),sL[2]*isqrt(sA[2]),sL[3]))

            @tensor rd_rho[:] :=L[-1,-3,1]*R[1,-2,-4]
        end
        return reshape(rd_rho,(2^length(r),2^length(r)))
    else
        sA=size(A.data[r[1]])
        M=reshape(A.data[r[1]],(sA[1],isqrt(sA[2]),isqrt(sA[2]),sA[3]))

        @tensor L[:] := L[1]*M[1,-1,-2,4]*R[4]
        return reshape(rd_rho,(2,2))
    end

end

max_stars_repo_name: aanniikkeettr/Tensor-Networks

max_issues_repo_name: aanniikkeettr/Tensor-Networks

max_forks_repo_name: aanniikkeettr/Tensor-Networks


text:
###############################################################################
#
# Wave length parameter set
#
###############################################################################
"""
    mutable struct WaveLengths{FT}

Struct for pre-set wave length parameters

# Fields
$(TYPEDFIELDS)
"""
Base.@kwdef mutable struct WaveLengths{FT}
    # Wave length (WL) boundaries
    "Minimal WL for PAR `[nm]`"
    minwlPAR::FT = FT(400.0)
    "Maximal WL for PAR `[nm]`"
    maxwlPAR::FT = FT(700.0)
    "Maximal WL for NIR `[nm]`"
    maxwlNIR::FT = FT(2500.0)
    "Minimal WL for SIF excitation `[nm]`"
    minwle  ::FT = FT(400.0)
    "Maximal WL for SIF excitation `[nm]`"
    maxwle  ::FT = FT(750.0)
    "Minimal WL for SIF emission/fluorescence `[nm]`"
    minwlf  ::FT = FT(640.0)
    "Maximal WL for SIF emission/fluorescence `[nm]` "
    maxwlf  ::FT = FT(850.0)

    # Wave length lists
    "Standard wave length `[nm]`"
    sWL::Array{FT,1} = [collect(FT(400.0):FT(10.0):FT( 650.1));
                        collect(FT(655.0):FT( 5.0):FT( 770.1));
                        collect(FT(780.0):FT(25.0):FT(2400.1))]
    "Differential wavelength"
    dWL::Array{FT,1} = diff(sWL)

    "Leaf optical parameter set"
    optis::LeafOpticals = LeafOpticals{FT}()
    "Wave length `[nm]`"
    WL  ::Array{FT,1}  = optis.lambda

    "Index of WLE in WL"
    iWLE::Array{Int,1} = findall( (WL .>= minwle) .& (WL .<= maxwle) )
    "Index of WLF in WL"
    iWLF::Array{Int,1} = findall( (WL .>= minwlf) .& (WL .<= maxwlf) )
    "index of wlPAR in WL"
    iPAR::Array{Int,1} = findall( (WL .>= minwlPAR) .& (WL .<= maxwlPAR) )
    "index of wlNIR in WL"
    iNIR::Array{Int,1} = findall( (WL .>= maxwlPAR) .& (WL .<= maxwlNIR) )

    "excitation wave length `[nm]`"
    WLE     ::Array{FT,1} = WL[iWLE]
    "Fluorescence wave length `[nm]`"
    WLF     ::Array{FT,1} = WL[iWLF]
    "Wave length for PAR"
    WL_iPAR ::Array{FT,1} = WL[iPAR];
    "Differential wave length for PAR"
    dWL_iPAR::Array{FT,1} = dWL[iPAR];
    "Differential wave length for iWLE"
    dWL_iWLE::Array{FT,1} = dWL[iWLE];

    # local storage of dimension information
    "Length of WL_iPAR"
    nPAR::Int = length(iPAR)
    "Length of WL"
    nWL ::Int = length(WL)
    "length of WLE"
    nWLE::Int = length(iWLE)
    "length of WLF"
    nWLF::Int = length(iWLF)
end

max_stars_repo_name: climate-machine/Land


text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:hydrogen
#     text_representation:
#       extension: .jl
#       format_name: hydrogen
#       format_version: '1.3'
#       jupytext_version: 1.11.2
#   kernelspec:
#     display_name: Julia 1.6.2
#     language: julia
#     name: julia-1.6
# ---

# %%
using Distributions
using StatsPlots

n = 10
L = 10^5
μ = 1.0
dist = Exponential(μ)
X̄ = [mean(rand(dist, n)) for _ in 1:L]
histogram(X̄; norm=true, alpha=0.3, label="sample means")
plot!(Gamma(n, μ/n); label="Gamma(n, μ/n)", lw=1.5)
title!("Samples: Exponential(μ = $μ),  n = $n")

# %%
using Distributions
using Random
using Roots

function ci_exp(n, x̄; α = 0.05)
    t0 = log(quantile(Gamma(n, x̄/n), 1 - α/2))
    L = find_zero(t -> cdf(Gamma(n, exp(t)/n), x̄) - (1 - α/2), t0) |> exp
    t0 = log(quantile(Gamma(n, x̄/n), α/2))
    R = find_zero(t -> cdf(Gamma(n, exp(t)/n), x̄) - α/2, t0) |> exp
    (L, R)
end

is_in(ci, x) = ci[1] ≤ x ≤ ci[2]

function prob_true_val_is_in_ci_exp(; n = 10, L = 10^5, μ = 1.0, α = 0.05)
    dist = Exponential(μ)
    tmp = zeros(n)
    prob = count(_ -> is_in(ci_exp(n, mean(rand!(dist, tmp)); α), μ), 1:L)/L
end

for n in (10, 20, 40, 100)
    @eval @show prob_true_val_is_in_ci_exp(n = $n)
end

# %%

max_stars_repo_path: 0021/confidence interval of the exponential dist.jl

max_stars_repo_name: genkuroki/public

max_issues_repo_name: genkuroki/public

max_forks_repo_name: genkuroki/public


text:
const FONTS = RelocatableFolders.@path joinpath(@__DIR__, "..", "..", "assets", "fonts")

function full_fontpath(fontname::AbstractString)
    isfile(fontname) && return fontname
    return joinpath(FONTS, fontname)
end

const _cached_fonts = Dict{String, FTFont}()

"""
    load_font(str)

Load a font. `str` must be either the full path to the font file, or a font in
the package font folder.

A font at a given location is cached for further use.
"""
function load_font(str)
    path = full_fontpath(str)
    get!(_cached_fonts, path) do
        FTFont(path)
    end
end

# Loading the font directly here lead to FreeTypeAbstraction to fail with error
# code 35, because handles to fonts are C pointer that cannot be fully
# serialized at compile time
const _default_fonts = Dict(
    :regular => joinpath("ComputerModern", "cmr10.ttf"),
    :italic => joinpath("ComputerModern", "cmmi10.ttf"),
    :bold => joinpath("ComputerModern", "cmb10.ttf"),
    :bolditalic => joinpath("ComputerModern", "cmmib10.ttf"),
    :math => joinpath("ComputerModern", "cmmi10.ttf")
)

const _default_font_mapping = Dict(
    :text => :regular,
    :delimiter => :regular,
    :digit => :regular,
    :function => :regular,
    :punctuation => :regular,
    :symbol => :math,
    :char => :italic
)

const _default_font_modifiers = Dict(
    :rm => Dict(:bolditalic => :bold, :italic => :regular),
    :it => Dict(:bold => :bolditalic, :regular => :italic),
    :bf => Dict(:italic => :bolditalic, :regular => :bold)
)

"""
    FontFamily([fonts, font_mapping, font_modifiers, slant_angle])

A set of font for LaTeX rendering.

# Fields
  - `font_mapping` a dict mapping the different character types (`:digit`,
    `:function`, `:punctuation`, `:symbol`, `:variable`) to a font identifier.
    Default to `MathTeXEngine._default_font_mapping`
  - `fonts` a dict mapping font identifier to a font path. Default to
    `MathTeXEngine._default_fonts` which represents the NewComputerModern font.
  - `font_modifiers` a dict of dict, one entry per font command supported in the
    font set. Each entry is a dict that maps a font identifier to another.
    Default to `MathTeXEngine._default_font_modifiers`.
  - `slant_angle` the angle by which the italic fonts are slanted, in degree
"""
struct FontFamily
    fonts::Dict{Symbol, String}
    font_mapping::Dict{Symbol, Symbol}
    font_modifiers::Dict{Symbol, Dict{Symbol, Symbol}}
    special_chars::Dict{Char, Tuple{String, Char}}
    slant_angle::Float64
    thickness::Float64
end

FontFamily(fonts) = FontFamily(
    fonts,
    _default_font_mapping,
    _default_font_modifiers,
    _symbol_to_computer_modern,
    15,
    0.0375)
FontFamily() = FontFamily(_default_fonts)

"""
    get_font([font_family=FontFamily()], fontstyle)

Get the FTFont object representing a font in the given font family. When called
with a single argument uses the default font family.
"""
function get_font(font_family::FontFamily, fontstyle::Symbol)
    return load_font(font_family.fonts[fontstyle])
end
get_font(fontstyle::Symbol) = get_font(FontFamily(), fontstyle)

"""
    get_fontpath([font_family::FontFamily], fontstyle)

Similar to `get_font` but return the path of the font instead of the FTFont
object.
"""
get_fontpath(font_family::FontFamily, fontstyle::Symbol) = full_fontpath(font_family.fonts[fontstyle])
get_fontpath(fontstyle::Symbol) = get_fontpath(FontFamily(), fontstyle)

function is_slanted(font_family, char_type)
    font_id = font_family.font_mapping[char_type]
    return font_id == :italic
end

slant_angle(font_family) = font_family.slant_angle * π / 180

# Few helper functions
"""
    thickness(font::FontFamily)

The thickness of the underline for the given font set.
"""
thickness(font_family::FontFamily) = font_family.thickness

"""
    xheight(font::FTFont)

The height of the letter x in the given font, i.e. the height of the letters
without neither ascender nor descender.
"""
xheight(font::FTFont) = inkheight(TeXChar('x', font))

"""
    xheight(font::FontFamily)

The height of the letter x in the given font family, i.e. the height of the letters
without neither ascender nor descender.
"""
xheight(font_family) = xheight(get_font(font_family, :regular))


text:
"""
    lon,lat,data = GeoDatasets.landseamask(;resolution='l',grid=5)

Load the land-sea-lake raster from basemap: 0 is ocean, 1 is land and 2 is lake.
`grid` is the resolution in arc minutes and should be either 1.25, 2.5, 5 or 10.
The parameter `resolution` should be either `'c'`,`'l'`,`'i'`,`'h'` or `'f'`
(standing for crude, low, intermediate, high and full resolution)

The data is downloaded from [basemap](https://github.com/matplotlib/basemap/tree/master/lib/mpl_toolkits/basemap_data).

The data is originally from [GMT](http://gmt.soest.hawaii.edu/) and distributed under the terms of the
[GPL licences](https://github.com/matplotlib/basemap/blob/a551c2361314670dc7d95143190d3d48ba750d20/LICENSE_data).

```julia
using GeoDatasets
lon,lat,data = GeoDatasets.landseamask(;resolution='c',grid=5)
```

"""
function landseamask(;resolution='l',grid=5)
    if !(grid ∈ Set([1.25,2.5,5.,10.]))
        error("grid should be either 1.25, 2.5, 5 or 10")
    end
    if !(resolution ∈ Set(['c','l','i','h','f',"c","l","i","h","f"]))
        error("resolution should be either 'c','l','i','h' or 'f'")
    end

    url = "https://raw.githubusercontent.com/matplotlib/basemap/v1.2.2rel/lib/mpl_toolkits/basemap/data/lsmask_$(grid)min_$(resolution).bin"
    @debug "url for landseamask: $url"

    # download remote file if it is not yet available
    @RemoteFile(lsmask, url)
    download(lsmask)

    # resolution in arc seconds
    res = round(Int,grid * 60)

    lon = (-180*3600 + res/2:res:180*3600 - res/2) / 3600
    lat = (-90*3600 + res/2:res:90*3600 - res/2) / 3600

    data = zeros(UInt8,length(lon)*length(lat));
    read!(GzipDecompressorStream(open(path(lsmask))),data);

    return lon,lat,reshape(data,length(lon),length(lat))
end

"""
    LandSeaMask

Encapsulates the `landseamask` data, and provides the `is_land` function.
"""
struct LandSeaMask{Tlon, Tlat, Tdata}
    lon::Tlon
    lat::Tlat
    data::Tdata
end

function LandSeaMask(; resolution='l', grid=5)
    return LandSeaMask(GeoDatasets.landseamask(; resolution=resolution, grid=grid)...)
end

"""
    is_land(mask::LandSeaMask, lat::Real, lon::Real)

`true` if the nearest lat-lon in the `mask` to `lat`-`lon` is land. Nearest is calculated
by flooring the lat and lon to the nearest point in the `mask` and return the value of that
entry.
"""
is_land(mask::LandSeaMask, lat::Real, lon::Real) = nearest_point(mask, lat, lon) == 1
is_land(mask::LandSeaMask, (lat, lon)::Tuple{Real, Real}) = is_land(mask, lat, lon)

function nearest_point(mask::LandSeaMask, lat::Real, lon::Real)
    return mask.data[nearest_point_idx(mask.lon, lon), nearest_point_idx(mask.lat, lat)]
end

function nearest_point_idx(xs::StepRangeLen, x::Real)
    return convert(Int, floor((x - first(xs)) / (last(xs) - first(xs)) * length(xs)))
end

max_stars_repo_name: JuliaGeo/GeoDatasets

max_issues_repo_name: JuliaGeo/GeoDatasets

max_forks_repo_name: JuliaGeo/GeoDatasets


text:
# blob.jl
import WaveletScattering: ScatteredBlob
# bank.jl
import WaveletScattering: Bank1D
# inputlayer.jl
import WaveletScattering: InputLayer, InputLayerState
# fourierlayer.jl
import WaveletScattering: FourierLayer, FourierLayerState
# morlet1d.jl
import WaveletScattering: Spec1D
# path.jl
import WaveletScattering: Path, PathKey
# node.jl
import WaveletScattering: RealFourierNode, InvComplexFourierNode
# waveletlayer.jl
import WaveletScattering: WaveletLayer
# invfourierlayer.jl
import WaveletScattering: InvFourierLayer
# pointwise.jl
import WaveletScattering: Modulus, PointwiseLayer, PointwiseLayerState
# poolinglayer.jl
import WaveletScattering: PoolingLayer
# scatteringlayer.jl
import WaveletScattering: ScatteringLayer

J = 8
data = zeros(Float32, 2^J, 2)
data[1] = 1.0f0
bank1 = Bank1D(Spec1D(log2_size=J, n_filters_per_octave=1))
bank2 = Bank1D(Spec1D(log2_size=J, n_filters_per_octave=1))

backend = Mocha.CPUBackend()
signal = InputLayer(
    data = data,
    tops = [:signal],
    symbols = [:time, :chunk])

fourier = FourierLayer(
    bottoms = [:signal],
    pathkeys = [PathKey(:time)],
    tops = [:fourier])

wavelets = WaveletLayer(
    bank = bank1,
    bottoms = [:fourier],
    tops = [:wavelets])

invfourier = InvFourierLayer(
    bottoms = [:wavelets],
    pathkeys = [PathKey(:time)],
    tops = [:invfourier])

modulus = PointwiseLayer(
    bottoms = [:invfourier],
    tops = [:modulus],
    ρ = Modulus())

fourier2 = FourierLayer(
    bottoms = [:modulus],
    name = "fourier2",
    pathkeys = [PathKey(:time)],
    tops = [:fourier2])

scattering = ScatteringLayer(
    bank = bank2,
    bottoms = [:fourier2],
    name = "scattering",
    sibling_mask_factor = 1.0,
    previous_bank = bank2,
    tops = [:scattering])

layers = Mocha.Layer[
    signal,
    fourier,
    wavelets,
    invfourier,
    modulus,
    fourier2]

Mocha.init(backend)
net = Mocha.Net("network", backend, layers)

@test isa(net, Mocha.Net{Mocha.CPUBackend})


text:
using FINUFFT
using Test

@testset "Error handling" begin

    xj = zeros(10)
    cj = complex(zeros(10))
    iflag = 1
    tol = 1e-15
    ms = 10

    nufft1d1(xj, cj, iflag, tol, ms) # This should pass

    @info("Testing error handling")

    # Tolerance too small
    try
        nufft1d1(xj, cj, iflag, 1e-100, ms)
    catch e
        @test e.errno == FINUFFT.ERR_EPS_TOO_SMALL
    end

    # Allocate too much
    opts = finufft_default_opts()
    opts.spread_kerevalmeth = 0
    opts.upsampfac = maxintfloat(typeof(opts.upsampfac))
    try
        nufft1d1(xj, cj, iflag, tol, ms, opts)
    catch e
        @test e.errno == FINUFFT.ERR_MAXNALLOC
    end

    # Too small upsampfac
    opts = finufft_default_opts()
    opts.spread_kerevalmeth = 0
    opts.upsampfac = 0
    try
        nufft1d1(xj, cj, iflag, tol, ms, opts)
    catch e
        @test e.errno == FINUFFT.ERR_UPSAMPFAC_TOO_SMALL
    end

    # Wrong beta
    opts = finufft_default_opts()
    opts.upsampfac = 0
    try
        nufft1d1(xj, cj, iflag, tol, ms, opts)
    catch e
        @test e.errno == FINUFFT.HORNER_WRONG_BETA
    end

end


text:
## R code 4.27
m4.1 <- quap(
    alist(
        height ~ dnorm( mu , sigma ) ,
        mu ~ dnorm( 178 , 20 ) ,
        sigma ~ dunif( 0 , 50 )
      )
)

## R code 4.31
m4.2 <- quap(
    alist(
        height ~ dnorm( mu , sigma ) ,
        mu ~ dnorm( 178 , 0.1 ) ,
        sigma ~ dunif( 0 , 50 )
    ) , data=d2 )

## R code 4.42
m4.3a <- quap(
    alist(
        height ~ dnorm( mu , sigma ) ,
        mu <- a + b*( weight - xbar ) ,
        a ~ dnorm( 178 , 20 ) ,
        b ~ dlnorm( 0 , 1 ) ,
        sigma ~ dunif( 0 , 50 )
    ) , data=d2 )

## R code 4.43
m4.3b <- quap(
    alist(
        height ~ dnorm( mu , sigma ) ,
        mu <- a + exp(log_b)*( weight - xbar ),
        a ~ dnorm( 178 , 20 ) ,
        log_b ~ dnorm( 0 , 1 ) ,
        sigma ~ dunif( 0 , 50 )
    ) ,
    data=d2 )

## R code 4.48
N <- 10
dN <- d2[ 1:N , ]
mN <- quap(
    alist(
        height ~ dnorm( mu , sigma ) ,
        mu <- a + b*( weight - mean(weight) ) ,
        a ~ dnorm( 178 , 20 ) ,
        b ~ dlnorm( 0 , 1 ) ,
        sigma ~ dunif( 0 , 50 )
    ) , data=dN )

## R code 4.65
m4.5 <- quap(
    alist(
        height ~ dnorm( mu , sigma ) ,
        mu <- a + b1*weight_s + b2*weight_s2 ,
        a ~ dnorm( 178 , 20 ) ,
        b1 ~ dlnorm( 0 , 1 ) ,
        b2 ~ dnorm( 0 , 1 ) ,
        sigma ~ dunif( 0 , 50 )
    ) ,
    data=d )

## R code 4.69
m4.6 <- quap(
    alist(
        height ~ dnorm( mu , sigma ) ,
        mu <- a + b1*weight_s + b2*weight_s2 + b3*weight_s3 ,
        a ~ dnorm( 178 , 20 ) ,
        b1 ~ dlnorm( 0 , 1 ) ,
        b2 ~ dnorm( 0 , 10 ) ,
        b3 ~ dnorm( 0 , 10 ) ,
        sigma ~ dunif( 0 , 50 )
    ) ,
    data=d )

## R code 4.76
m4.7 <- quap(
    alist(
        T ~ dnorm( mu , sigma ) ,
        mu <- a + B %*% w ,
        a ~ dnorm(6,10),
        w ~ dnorm(0,1),
        sigma ~ dexp(1)
    ),
    data=list( T=d2$temp , B=B ) ,
    start=list( w=rep( 0 , ncol(B) ) ) )

## R code 4.79
m4.7alt <- quap(
    alist(
        T ~ dnorm( mu , sigma ) ,
        mu <- a + sapply( 1:1124 , function(i) sum( B[i,]*w ) ) ,
        a ~ dnorm(6,10),
        w ~ dnorm(0,1),
        sigma ~ dexp(1)
    ),
    data=list( T=d2$temp , B=B ) ,
    start=list( w=rep( 0 , ncol(B) ) ) )


text:
module NaiveBayes

using Distributions
using HDF5
using KernelDensity
using Interpolations
using LinearAlgebra
using StatsBase
import StatsBase: fit, predict

export  NBModel,
        MultinomialNB,
        GaussianNB,
        KernelNB,
        HybridNB,
        fit,
        predict,
        predict_proba,
        predict_logprobs,
        restructure_matrix,
        to_matrix,
        write_model,
        load_model,
        get_feature_names,
        train

include("nbtypes.jl")
include("common.jl")
include("hybrid.jl")
include("gaussian.jl")
include("multinomial.jl")

end


text:
mutable struct Controller4 <: Controller
    Δt::Real
    g::Real
    M::Real

    c::Real
    k::Real
    Fs::Real

    lqr::LQR
    F::Function
    control!::Function

    function Controller4(F,Δt,g,M,c,k,Fs,lqr)
        new(Δt,g,M,c,k,Fs,lqr,F,control!)
    end
end

function control!(mechanism::Mechanism{T,Nn,Nb},controller::Controller4,kstep) where {T,Nn,Nb}
    Δt = controller.Δt
    g = controller.g
    M = controller.M
    c = controller.c
    k = controller.k
    Fs = controller.Fs
    lqr = controller.lqr
    F = controller.F
    cart = mechanism.bodies[1]
    pend1 = mechanism.bodies[2]
    pend2 = mechanism.bodies[3]

    zdes = F(zeros(6),I,kstep*Δt)
    ydes = zdes[1]
    xd1 = lqr.xd[1]
    xd2 = lqr.xd[2]
    xd3 = lqr.xd[3]
    lqr.xd[1] = [xd1[1];ydes;xd1[3]]
    lqr.xd[2] = [xd2[1];ydes;xd2[3]]
    lqr.xd[3] = [xd3[1];ydes;xd3[3]]

    Δz = zeros(T,Nb*12)
    qvm = QuatVecMap()
    for (id,body) in pairs(mechanism.bodies)
        colx = (id-1)*12+1:(id-1)*12+3
        colv = (id-1)*12+4:(id-1)*12+6
        colq = (id-1)*12+7:(id-1)*12+9
        colω = (id-1)*12+10:(id-1)*12+12

        state = body.state
        Δz[colx] = state.xsol[2]-lqr.xd[id]
        Δz[colv] = state.vsol[2]-lqr.vd[id]
        Δz[colq] = rotation_error(state.qsol[2],lqr.qd[id],qvm)
        Δz[colω] = state.ωsol[2]-lqr.ωd[id]
    end

    v = cart.state.vc[2]
    ω1 = pend1.state.ωc[1]
    ω2 = pend2.state.ωc[1]

    Fcart = -(lqr.K[1][1]*Δz)[1]
    if abs(Fcart)<Fs
        Fcart = sign(Fcart)*(Fs+0.1)
    end
    Fcart, Fpend1, Fpend2 = friction(Fcart, 0, 0, v, ω1, ω2, k, c, Fs)

    setForce!(mechanism, geteqconstraint(mechanism, 4), [Fcart])
    setForce!(mechanism, geteqconstraint(mechanism, 5), [Fpend1])
    setForce!(mechanism, geteqconstraint(mechanism, 6), [Fpend2])

    return
end

function run4!(str; vis::Bool=true)
    open("Files/E4.jl",create=true,write=true) do file
        write(file,"
        Δt = 0.01

        # Parameters from Table 1.2
        g = 9.81
        m1 = 0.135
        m2 = 0.123
        l1 = 0.140
        l2 = 0.135
        L1 = 0.280
        L2 = 0.270
        # I1
        # I2
        M = 4.5
        n1 = 0.026
        n2 = 0.010
        c = 0.002
        k = 0.3
        Fs = 8

        # Joint axes
        ex = [1.0;0.0;0.0]
        ey = [0.0;1.0;0.0]

        p1 = [0.0;0.0;l1] # joint connection point
        p2 = [0.0;0.0;l2] # joint connection point

        # Desired orientation
        θ1 = 0
        θ2 = 0

        # Links
        origin = Origin{Float64}()
        cart = Box(0.1, 0.5, 0.1, M)
        pend1 = Box(0.055, 0.055/2, L1, m1, color = RGBA(1,0,0))
        pend2 = Box(0.055, 0.055/2, L2, m2)

        # Constraints
        joint1 = EqualityConstraint(Prismatic(origin, cart, ey))
        joint2 = EqualityConstraint(Revolute(cart, pend1, ex; p2 = -p1))
        joint3 = EqualityConstraint(Revolute(pend1, pend2, ex; p1 = p1, p2 = -p2))

        links = [cart;pend1;pend2]
        constraints = [joint1;joint2;joint3]

        mech = Mechanism(origin, links, constraints, Δt = Δt)

        xinit = 0.0
        theta1init = pi
        theta2init = 0
        Q = diagm(ones(6))
        R = 0.1
        F(z,K,t) = K*zeros(6)

        ### Begin Student Input

        ")
        write(file, str)
        write(file, "

        ### End Student Input

        setPosition!(origin,cart,Δx = [0;xinit;0])
        setPosition!(cart,pend1,p2 = -p1, Δq = QuatRotation(RotX(theta1init)))
        setPosition!(pend1,pend2,p1 = p1, p2 = -p2, Δq = QuatRotation(RotX(theta2init)))

        xd = [[[0;0.0;0.0]]; [p1]; [2*p1+p2]]

        Q1 = diagm([0;Q[1,1];0; 0;Q[2,2];0; 0;0;0; 0;0;0])
        Q2 = diagm([0;0;0; 0;0;0; Q[3,3]*4;0;0; Q[4,4];0;0])
        Q3 = diagm([0;0;0; 0;0;0; Q[5,5]*4;0;0; Q[6,6];0;0])
        Q = [[Q1];[Q2];[Q3]]
        R = [ones(1,1)*R[1]]

        lqr = LQR(mech, getid.(links), [getid(constraints[1])], Q, R, Inf, xd=xd)

        controller = Controller4(F,Δt,g,M,c,k,Fs,lqr)

        steps = Base.OneTo(Int(20/Δt))
        storage = Storage{Float64}(steps,3)

        try
            simulate!(mech,storage,controller,record = true)
        catch
            println(\"Unstable behavior\")
        end")

        if vis
            write(file, "
            visualize(mech,storage)")
        end
    end

    try
        include("Files/E4.jl")
    catch
        println("Error. Bad code.")
    end

    return
end

max_stars_repo_name: janbruedigam/ACRLab

max_issues_repo_name: janbruedigam/ACRLab

max_forks_repo_name: janbruedigam/ACRLab


text:
using Plots
using Random
using StatsBase

#=
 所属ID 0      家庭のみ
        1-40   社会組織
        41-60  学校

 地域ID 1-10   住宅地
        11-15  オフィス街
=#

NumberOfAgent   = 10000
NumberOfCompany = 40

mutable struct Agent
    infulState       :: UInt8
    org_id           :: UInt32
    count_of_contact :: UInt32
end

struct Family
    region     :: UInt16
    member_ids :: Array{UInt32,1}
end

struct Company
    region     :: UInt16
    member_ids :: Array{UInt32,1}
    table      :: Array{UInt32,2}
end

struct School
    member_ids :: Array{UInt32,1}
    table      :: Array{UInt32,2}
end

struct Vehicle
    passenger::Array{UInt32, 2} # 4x200 = 800
end

function initInhabitans(NumberOfAgent::Int64)

    inhabitants = Array{Agent,1}()
    families = Array{Family,1}()

    sizehint!(inhabitants, NumberOfAgent)

    parent_org_ids = zeros(UInt32,5000)

    # 所属IDを作成
    for i in 1:length(parent_org_ids)-1000
        parent_org_ids[i] = rem(i,40)+1 # 1-40
    end

    shuffle!(parent_org_ids)
    let
        hist = zeros(UInt32,5000)

        for i in 1:length(parent_org_ids)
            hist[parent_org_ids[i]+1]=hist[parent_org_ids[i]+1]+1
        end
        for i in 1:41
            j = i-1
            histj=hist[i]
            #println("$j,$histj")
        end
        #println(hist)
    end

    # 地域は1-10
    j=1
    for i::Int64 in 1:div(NumberOfAgent,4)

        region = div(i-1,250)+1
        # parent1
        push!(inhabitants, Agent(0,parent_org_ids[j],0))
        j += 1
        # parent2
        push!(inhabitants, Agent(0,parent_org_ids[j],0))
        j += 1

        # child 1
        push!(inhabitants, Agent(0,41+(region-1)*2,0))

        # child 2
        push!(inhabitants, Agent(0,41+(region-1)*2+1,0))

        push!(families,
              Family(region, [(i-1)*4+1,(i-1)*4+2,(i-1)*4+3,(i-1)*4+4]))
    end

    println("num of family: $(length(families))")

    #=
    num_family_in_region = zeros(Int64,10);
    for i in 1:length(families)
        println("family $i -> region: $(families[i].region)")
        num_family_in_region[families[i].region] += 1
    end

    for i in 1:length(num_family_in_region)
        println("region $i -> num of family : $(num_family_in_region[i])")
    end

    num_family_in_region = zeros(Int64,10)
    for i in 1:length(families)
        num_family_in_region[families[i].region] += 4
    end

    num_students_in_school = zeros(Int64,20)
    for i in 1:length(inhabitants)
        if 41 <= inhabitants[i].org_id <= 60
            num_students_in_school[inhabitants[i].org_id-40]+=1
        end
    end

    for i in 1:length(num_students_in_school)
        println("school $i -> num of student : $(num_students_in_school[i])")
    end
    =#

    return inhabitants, families;
end

function initCompany(inhabitans)

    companies = Array{Company,1}()

    for i in 1:NumberOfCompany
        push!(companies,Company(div(i-1,8)+11,[],zeros(UInt32, 10, 10)))
    end

    for i in 1:length(inhabitans)
        if 0 < inhabitans[i].org_id <= 40
            push!(companies[inhabitans[i].org_id].member_ids,i)
        end
    end

    # 席配置
    for i in 1:length(companies)
        #println("company $i -> region $(companies[i].region)")
        l = 1
        for j=1:10, k=1:10
            companies[i].table[j,k] = companies[i].member_ids[l]
            l = l + 1
        end
    end

    return companies
end

function initSchool(inhabitants)
    schools = Array{School,1}()

    for i in 1:20
        push!(schools,School([],zeros(UInt32,10,25)))
    end

    for i in 1:length(inhabitants)
        if 41 <= inhabitants[i].org_id <= 60
            schoolid = inhabitants[i].org_id-40
            push!(schools[schoolid].member_ids,i)
        end
    end

    for i in 1:length(schools)
        m = copy(schools[i].member_ids)
        shuffle!(m)
        for j=1:10,k=1:25
            schools[i].table[j,k] = m[(j-1)*25+k];
        end
    end

    return schools;
end

function ride_vehicles(companies::Array{Company,1},inhabitants::Array{Agent,1})
    vehicles = Array{Vehicle,1}()
    to_region_passengers = Array{Array{UInt32,1},1}()
    #num_of_p_in_region = zeros(UInt32, 5)

    for i in 11:15
        push!(to_region_passengers,[])
    end

    for i in 11:15, j in 1:length(companies)
        if companies[j].region == i
            #println("company $j -> $(companies[j].region)")
            #num_of_p_in_region[i-10] += length(companies[j].member_ids)
            to_region_passengers[i-10] = append!(to_region_passengers[i-10],companies[j].member_ids)
        end
    end
    #=
    for i in 1:length(to_region_passengers)
        println("region $i -> $(length(to_region_passengers[i]))")
    end
    =#
    for i in 1:5
        shuffle!(to_region_passengers[i])
        push!(vehicles,Vehicle(reshape(to_region_passengers[i],4,200)))
    end
    return vehicles
end

function create_familycell_array(families,inhabitans)
    familycell_array = zeros(UInt8,100,100)
    id = 1
    for i=1:50,j=1:50
        #println("($i,$j)")
        familycell_array[2*(i-1)+1,2*(j-1)+1] = inhabitans[families[id].member_ids[1]].infulState
        familycell_array[2*(i-1)+2,2*(j-1)+1] = inhabitans[families[id].member_ids[2]].infulState
        familycell_array[2*(i-1)+1,2*(j-1)+2] = inhabitans[families[id].member_ids[3]].infulState
        familycell_array[2*(i-1)+2,2*(j-1)+2] = inhabitans[families[id].member_ids[4]].infulState
        id = id + 1
    end
    return familycell_array
end

function create_company_plots(companies::Array{Company,1},inhabitants::Array{Agent,1})

    company_plots = Array{Plots.Plot,1}()

    for i in 1:length(companies)
        cell = zeros(UInt8,10,10)
        for j=1:10, k=1:10
            cell[j,k] = inhabitants[companies[i].table[j,k]].infulState
        end
        push!(company_plots,heatmap(cell,
                                    c=cgrad(:matter,3),
                                    legend=false,
                                    xlims=(1,10),
                                    title="company $i"))
    end

    println("plots = $(length(company_plots))")

    return plot(company_plots[1],company_plots[2], company_plots[3], company_plots[4], company_plots[5],
                company_plots[6],company_plots[7], company_plots[8], company_plots[9], company_plots[10],
                company_plots[11],company_plots[12], company_plots[13], company_plots[14], company_plots[15],
                company_plots[16],company_plots[17], company_plots[18], company_plots[19], company_plots[20],
                company_plots[21],company_plots[22], company_plots[23], company_plots[24], company_plots[25],
                company_plots[26],company_plots[27], company_plots[28], company_plots[29], company_plots[30],
                company_plots[31],company_plots[32], company_plots[33], company_plots[34], company_plots[35],
                company_plots[36],company_plots[37], company_plots[38], company_plots[39], company_plots[40],
                layout=(5,8))
end

function create_school_plots(schools::Array{School,1}, inhabitants)
    school_plots = Array{Plots.Plot,1}()
    for i in 1:length(schools)
        cell = zeros(UInt8,10,25)
        for j=1:10, k=1:25
            cell[j,k] = inhabitants[schools[i].table[j,k]].infulState
        end
        push!(school_plots,heatmap(cell,
                                   c=cgrad(:matter,3),
                                   legend=false,
                                   title="school $i"))
    end
    return plot(school_plots[1],school_plots[2],school_plots[3],school_plots[4],
                school_plots[5],school_plots[6],school_plots[7],school_plots[8],
                school_plots[9],school_plots[10],school_plots[11],school_plots[12],
                school_plots[13],school_plots[14],school_plots[15],school_plots[16],
                school_plots[17],school_plots[18],school_plots[19],school_plots[20],
                layout=(5,4))
end

function create_vehilcle_plots(vehilcles::Array{Vehicle,1}, inhabitants::Array{Agent,1})

    vehilcle_plots = Array{Plots.Plot,1}()
    for i in 1:length(vehilcles)
        cell = zeros(UInt8,4,200)
        for j=1:4, k=1:200
            cell[j,k] = inhabitants[vehilcles[i].passenger[j,k]].infulState
        end
        push!(vehilcle_plots,heatmap(cell,
                                     c=cgrad(:matter,3),
                                     legend=false,
                                     title="vehilcle $i"))
    end
    println(length(vehilcle_plots))
    return plot(vehilcle_plots[1],
                vehilcle_plots[2],
                vehilcle_plots[3],
                vehilcle_plots[4],
                vehilcle_plots[5],
                layout = (5,1))
end

inhabitants,families = initInhabitans(NumberOfAgent)

inhabitants[1].infulState = 1
inhabitants[2].infulState = 1
inhabitants[3].infulState = 1
inhabitants[4].infulState = 1

#println(familycell_array)

companies = initCompany(inhabitants)
schools = initSchool(inhabitants)

for d in 1:1
    # 3時間が通勤時間
    # 8時間が会社
    # 学校を11時間
    # 家を13時間
    vehilcles = ride_vehicles(companies, inhabitants)
end

family_plot = heatmap(create_familycell_array(families,inhabitants),
                      c=cgrad(:matter,3,categorical=true),
                      legend=false,
                      gridstyle=:solid,
                      title="family")

company_plot  = create_company_plots(companies, inhabitants)

school_plot   = create_school_plots(schools, inhabitants)

vehilcle_plot = create_vehilcle_plots(vehilcles, inhabitants)

plot(family_plot,
     school_plot,
     company_plot,
     vehilcle_plot,
     layout=@layout[a c; b; d],legend=false,size=(2600,2200))

#=
anim = Animation()

for i in 1:10
    graph_plot = plot(t->sinpi(t+i/5), range(0, 2, length=100), label="sin(t)")
    plot!(graph_plot,t->cospi(t+i/5), range(0, 2, length=100), label="cos(t)")
    frame(anim, graph_plot)
end

gif(anim, "plot.gif", fps=15)
=#

max_stars_repo_name: seisyuu-hantatsushi/julia

max_issues_repo_name: seisyuu-hantatsushi/julia

max_forks_repo_name: seisyuu-hantatsushi/julia


text:
module TestEmptyProxy

include("preamble.jl")

==ₜₗ(_, _) = false
==ₜₗ(x::T, y::T) where {T} = length(x) == length(y)

@testset begin
    @test collect(Empty(Vector)) == Union{}[]
    @test collect(Empty(Dict)) == Union{}[]
    @test collect(Empty(Set)) == Union{}[]
end

@testset "similar" begin
    @test similar(Empty(Vector)) ==ₜ Union{}[]
    @test similar(Empty(Vector{Union{}})) ==ₜ Union{}[]
    @test similar(Empty(Vector{Int})) ==ₜ Int[]
    @test similar(Empty(Vector), Int) ==ₜ Int[]
    @test similar(Empty(Vector{Union{}}), Int) ==ₜ Int[]
    @test similar(Empty(Vector{Int}), Symbol) ==ₜ Symbol[]

    @test similar(Empty(Vector), 3) ==ₜₗ resize!(Union{}[], 3)
    @test similar(Empty(Vector{Union{}}), 3) ==ₜₗ resize!(Union{}[], 3)
    @test similar(Empty(Vector{Int}), 3) ==ₜₗ resize!(Int[], 3)
    @test similar(Empty(Vector), Int, 3) ==ₜₗ resize!(Int[], 3)
    @test similar(Empty(Vector{Union{}}), Int, 3) ==ₜₗ resize!(Int[], 3)
    @test similar(Empty(Vector{Int}), Symbol, 3) ==ₜₗ resize!(Symbol[], 3)

    @test similar(Empty(Vector), (3,)) ==ₜₗ resize!(Union{}[], 3)
    @test similar(Empty(Vector{Union{}}), (3,)) ==ₜₗ resize!(Union{}[], 3)
    @test similar(Empty(Vector{Int}), (3,)) ==ₜₗ resize!(Int[], 3)
    @test similar(Empty(Vector), Int, (3,)) ==ₜₗ resize!(Int[], 3)
    @test similar(Empty(Vector{Union{}}), Int, (3,)) ==ₜₗ resize!(Int[], 3)
    @test similar(Empty(Vector{Int}), Symbol, (3,)) ==ₜₗ resize!(Symbol[], 3)
end

@testset "Empty(Dict)" begin
    @test !haskey(Empty(Dict), :key)
end

end  # module


text:
module FEMMPointGroundedSpringModule

using LinearAlgebra: norm, Transpose, mul!
using FinEtools
using FinEtools.IntegDomainModule: IntegDomain
import FinEtoolsDeforLinear.MatDeforElastIsoModule: MatDeforElastIso
using ..FESetCorotBeamModule: FESetL2CorotBeam, local_frame_and_def!, local_mass!, local_stiffness!, natural_forces!, local_geometric_stiffness!, local_forces!, MASS_TYPE_CONSISTENT_WITH_ROTATION_INERTIA

"""
    FEMMPointGroundedSpring{S<:AbstractFESet}

Type for finite element modeling machine for grounded springs.
"""
mutable struct FEMMPointGroundedSpring{S<:AbstractFESet} <: AbstractFEMM
    integdomain::IntegDomain{S} # integration domain data
    _stiffmatrix::FFltMat
end

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

Compute the stiffness matrix
"""
function stiffness(self::FEMMPointGroundedSpring, assembler::ASS, geom0::NodalField{FFlt}, u1::NodalField{T}, Rfield1::NodalField{T}, dchi::NodalField{TI}) where {ASS<:AbstractSysmatAssembler, T<:Number, TI<:Number}
    fes = self.integdomain.fes
    dofnums = zeros(FInt, 1, 6);
    startassembly!(assembler, 6, 6, count(fes), dchi.nfreedofs, dchi.nfreedofs);
    for i = 1:count(fes) # Loop over elements
        gatherdofnums!(dchi, dofnums, fes.conn[i]); # degrees of freedom
        assemble!(assembler, self._stiffmatrix, dofnums, dofnums);
    end # Loop over elements
    return makematrix!(assembler);
end

function stiffness(self::FEMMPointGroundedSpring, geom0::NodalField{FFlt}, u1::NodalField{T}, Rfield1::NodalField{T}, dchi::NodalField{TI}) where {T<:Number, TI<:Number}
    assembler = SysmatAssemblerSparseSymm();
    return stiffness(self, assembler, geom0, u1, Rfield1, dchi);
end

end # module


text:
function Initi_model_struct(H::T,Hc::T,N::T,Wf::T,B::T,αf::T,β₁::T,β₂::T,
    Nₛ::T,rₘ::T,a_Jmax::T,b_Jmax::T,i::T,Kₓₗ₀::T) where {T<:Float64}
    cons = Constants()
    env = EnvironmentStruct()
    kinetic = PhotoKineticRates()
    photo = PhotoPar(kinetic,env.Tₐ)
    hydPar = HydraulicsPar(;Kₓₗ₀=Kₓₗ₀,i=i)
    treepar = TreePar(;αf=αf,β₁=β₁,β₂=β₂,rₘ=rₘ,Nₛ=Nₛ,a_Jmax=a_Jmax,b_Jmax=b_Jmax)

    Hs = H-Hc
    As = Wf/treepar.αf
    Ww = treepar.ρw*As*(treepar.β₁*H+treepar.β₂*Hs)
    treesize = TreeSize(Wf,Ww,H,Hs,As,B,N)

    model = CCPHStruct(cons,env,treepar,treesize,photo,hydPar)

    return model,kinetic
end

function Initi_model_struct(j::Integer,data::RoData,
    αf::T,β₁::T,β₂::T,Nₛ::T,rₘ::T,a_Jmax::T,b_Jmax::T,i::T,Kₓₗ₀::T) where {T<:Float64}
    data_ind = Find_data_ind(j)
    H = data.H[data_ind]
    N = data.N[data_ind]
    Wf = data.Wf[data_ind]
    B = data.B[data_ind]
    Hc = data.Hc
    model,kinetic = Initi_model_struct(H,Hc,N,Wf,B,αf,β₁,β₂,Nₛ,rₘ,a_Jmax,b_Jmax,i,Kₓₗ₀)

    return model,kinetic
end

function Run_RO_C_F_CCPH(par::Array{Float64,1},weatherts_F::WeatherTS,
    weatherts_C::WeatherTS,data_F::RoData,data_C::RoData;sim_steps::Integer=83)

    αf,β₁,β₂,b_Jmax = 460.0,1.27,-0.27,0.0
    Nₛ,rₘ,a_Jmax,Kₓₗ₀,Nₛ_C,i = par[1:6]

    #Fertilized stand
    GPP_model_F = zeros(sim_steps+1)
    EC_model_F = zeros(sim_steps+1)
    for j = 1:sim_steps+1

        model,kinetic = Initi_model_struct(j,data_F,αf,β₁,β₂,Nₛ,rₘ,a_Jmax,b_Jmax,i,Kₓₗ₀)

        GPP_model_F[j],EC_model_F[j] = CalcModelOutput!(j,model,weatherts_F,kinetic)
    end

    #Control stand
    GPP_model_C = zeros(sim_steps+1)
    EC_model_C = zeros(sim_steps+1)
    for j = 1:sim_steps+1

        model,kinetic = Initi_model_struct(j,data_C,αf,β₁,β₂,Nₛ_C,rₘ,a_Jmax,b_Jmax,i,Kₓₗ₀)

        GPP_model_C[j],EC_model_C[j] = CalcModelOutput!(j,model,weatherts_C,kinetic)
    end

    return GPP_model_F,GPP_model_C,EC_model_F,EC_model_C
end

function calc_logP_term(ydata::T,ymodel::T,a::T,b::T) where {T<:Float64}
    e = abs(ydata-ymodel)
    σ = a+b*ymodel
    σ>0.0||error("Negative variance")
    return e/σ+log(σ)
end

function Calc_logP(
    GPP_model_F::Array{T,1},Ec_model_F::Array{T,1},data_F::RoData,
    GPP_model_C::Array{T,1},Ec_model_C::Array{T,1},data_C::RoData,
    a_GPP::T,b_GPP::T,a_Ec::T,b_Ec::T) where {T<:Float64}

    logP = 0.0
    for i = 1:length(data_F.GPP)
        logP += calc_logP_term(data_F.GPP[i],GPP_model_F[i],a_GPP,b_GPP)
        logP += calc_logP_term(data_C.GPP[i],GPP_model_C[i],a_GPP,b_GPP)

        logP += calc_logP_term(data_F.E_C[i],Ec_model_F[i],a_Ec,b_Ec)
        logP += calc_logP_term(data_C.E_C[i],Ec_model_C[i],a_Ec,b_Ec)
    end

    return  -logP
end

function Check_para_RO_C_F_CCPH(par::Array{Float64,1})
    #Nₛ,rₘ,a_Jmax,Kₓₗ₀,Nₛ_C,i,X0,τ,Smax,τ_C,a_GPP,b_GPP,a_Ec,b_Ec
    lo = [0.001,10.0,0.01,0.001,0.001,0.1,-18.0,1.0,10.0,1.0,0.0001,0.0001,0.0001,0.0001]

    up = [0.1,40.0,1.0,0.1,0.1,6.0,-0.9,15.0,25.0,15.0,5.0,3.0,5.0,3.0]

    return any(par.<lo)||any(par.>up)
end

function Post_distri_RO_C_F_CCPH(par::Array{Float64,1},RO_data::RO_raw_data)
    if Check_para_RO_C_F_CCPH(par)
        return post = -Inf
    else
        try
            X0,τ,Smax,τ_C = par[7:10]
            weatherts_F,GPP_data_F = create_weather_struct_RO(RO_data;stand_type="Fertilized",X0=X0,τ=τ,Smax=Smax)
            weatherts_C,GPP_data_C = create_weather_struct_RO(RO_data;stand_type="Control",X0=X0,τ=τ_C,Smax=Smax)

            data_F,data_C = Create_RoData_C_F(GPP_data_F,GPP_data_C,weatherts_F,weatherts_C)

            GPP_model_F,GPP_model_C,EC_model_F,EC_model_C = Run_RO_C_F_CCPH(par,weatherts_F,weatherts_C,data_F,data_C)

            a_GPP,b_GPP,a_Ec,b_Ec = par[11:14]

            logP = Calc_logP(
            GPP_model_F,EC_model_F,data_F,
            GPP_model_C,EC_model_C,data_C,
            a_GPP,b_GPP,a_Ec,b_Ec)

            return post = logP
        catch err
            println("Parameters Error: ", err)
            return post = -Inf
        end
    end
end

function test_adaptive_rwm()

    file_name = "adaptive_rwm_RO_C_GPP_Ec_20220228"
    n_chains = 4
    n_iter = 30000

    RO_data = Load_RO_data()

    #=
    par_guess = [0.011976551772617215, 25.42466609569836, 0.5009242298864488, 0.023503724920524185,
    0.010452800544803813,1.0,-1.312709470268466, 3.3784814235281933, 17.510523050796454, 9.566445194647766,
    0.237,0.0845,0.0688,0.146]
    =#

    par_guess = [0.004320402400586008, 27.297822676817912, 0.029549759752040473, 0.0693215277873528,
    0.0023441457252274013, 2.2501166086478843, -6.776410663775599, 13.19517957686665,
    20.362175060472445, 2.6466352227229257, 0.4645185292099437, 0.08043063207519931,
    0.14985616067841498, 0.1500009201257627]

    @show Post_distri_RO_C_F_CCPH(par_guess,RO_data)

    out_mat_tot = Array{Matrix{Float64},1}(undef,n_chains)

    Threads.@threads for j = 1:n_chains
        # Run n_iter iterations of the Adaptive Metropolis method:
        out = AdaptiveMCMC.adaptive_rwm(par_guess,
        x::Array{Float64,1}->Post_distri_RO_C_F_CCPH(x::Array{Float64,1},RO_data),
        n_iter; algorithm=:ram)

        out_mat_tot[j] = transpose(out.X)
    end

    try
        save("./output/"*file_name*".jld","samples",out_mat_tot)

        @show typeof(out_mat_tot)
        @show size(out_mat_tot)

        n_features = size(out_mat_tot[1])[2]

        histogram(out_mat_tot[1],layout=n_features,legend=false)
        savefig("./plots/"*file_name*".svg")

        plot(;layout=n_features,legends=false)
        for i = 1:n_features
            for j=1:n_chains
                plot!(out_mat_tot[j][:,i],subplot=i)
            end
        end
        savefig("./plots/"*file_name*"_trace_plot.svg")

        # Calculate '95% credible intervals':
        para_stat = mapslices(x->"$(mean(x)) ± $(1.96std(x))", out_mat_tot[1], dims=1)
        for line in para_stat
            println(line)
        end
    catch err
        println("Error: ",err)
    end
end

function find_max_adaptive_rwm()

    file_name = "adaptive_rwm_RO_C_GPP_Ec_20220228"
    out_mat_tot = load("./output/"*file_name*".jld","samples")
    out_mat = out_mat_tot[1]
    n_features = size(out_mat)[2]
    n_samples = size(out_mat)[1]

    RO_data = Load_RO_data()

    #=
    log_p_vec = Array{Float64,1}(undef,n_samples)

    Threads.@threads for j = 1:n_samples
        println("Done: $(j)/$(n_samples)")
        para_val = out_mat[j,:]
        log_p_vec[j] = Post_distri_RO_C_F_CCPH(para_val,RO_data)
    end

    save("./output/"*file_name*"_log_p.jld","log_p",log_p_vec)
    =#
    log_p_vec = load("./output/"*file_name*"_log_p.jld","log_p")

    @show min_ind = argmax(log_p_vec)
    @show log_p_vec[min_ind]
    @show par = reshape(mean(out_mat, dims=1),n_features)
    @show par = out_mat[min_ind,:]

    @show Post_distri_RO_C_F_CCPH(par,RO_data)

    X0,τ,Smax,τ_C = par[7:10]
    weatherts_F,GPP_data_F = create_weather_struct_RO(RO_data;stand_type="Fertilized",X0=X0,τ=τ,Smax=Smax)
    weatherts_C,GPP_data_C = create_weather_struct_RO(RO_data;stand_type="Control",X0=X0,τ=τ_C,Smax=Smax)

    data_F,data_C = Create_RoData_C_F(GPP_data_F,GPP_data_C,weatherts_F,weatherts_C)

    GPP_model_F,GPP_model_C,EC_model_F,EC_model_C = Run_RO_C_F_CCPH(par,weatherts_F,weatherts_C,data_F,data_C)

    a_GPP,b_GPP,a_Ec,b_Ec = par[11:14]

    @show GPP_R²_F =  calcR²(data_F.GPP,GPP_model_F)
    @show GPP_cor_F = cor(data_F.GPP,GPP_model_F)
    @show gₛ_R²_F =  calcR²(data_F.E_C,EC_model_F)
    @show gₛ_cor_F = cor(data_F.E_C,EC_model_F)

    @show GPP_R²_C =  calcR²(data_C.GPP,GPP_model_C)
    @show GPP_cor_C = cor(data_C.GPP,GPP_model_C)
    @show gₛ_R²_C =  calcR²(data_C.E_C,EC_model_C)
    @show gₛ_cor_C = cor(data_C.E_C,EC_model_C)

    println("Fertilized")
    clac_GPP_R²_annual(data_C.GPP,GPP_model_F,weatherts_F.date)
    println("Control")
    clac_GPP_R²_annual(data_C.GPP,GPP_model_C,weatherts_C.date)

    σ_GPP_F = a_GPP.+GPP_model_F*b_GPP
    σ_GPP_C = a_GPP.+GPP_model_C*b_GPP

    σ_EC_F = a_Ec.+b_Ec*EC_model_F
    σ_EC_C = a_Ec.+b_Ec*EC_model_C

    c = -log(0.025*2) #Use for calculating 95% credible intervals

    plot(weatherts_F.date,data_F.GPP,label="Data",ylabel="GPP")
    plot!(weatherts_F.date,GPP_model_F+c*σ_GPP_F)
    plot!(weatherts_F.date,GPP_model_F-c*σ_GPP_F)
    pl1 = plot!(weatherts_F.date,GPP_model_F,label="Model")
    plot(weatherts_C.date,data_C.GPP,label="Data",ylabel="GPP")
    plot!(weatherts_C.date,GPP_model_C+c*σ_GPP_C)
    plot!(weatherts_C.date,GPP_model_C-c*σ_GPP_C)
    pl2 = plot!(weatherts_C.date,GPP_model_C,label="Model")

    plot(weatherts_F.date,data_F.E_C,label="Data",ylabel="E_C")
    plot!(weatherts_F.date,EC_model_F+c*σ_EC_F)
    plot!(weatherts_F.date,EC_model_F-c*σ_EC_F)
    pl3 = plot!(weatherts_F.date,EC_model_F,label="Model")
    plot(weatherts_C.date,data_C.E_C,label="Data",ylabel="E_C")
    plot!(weatherts_C.date,EC_model_C+c*σ_EC_C)
    plot!(weatherts_C.date,EC_model_C-c*σ_EC_C)
    pl4 = plot!(weatherts_C.date,EC_model_C,label="Model")

    plot(pl1,pl3,pl2,pl4,layout=(2,2),legends=false)
    savefig("./plots/"*file_name*"_result.svg")

    # Calculate '95% credible intervals':
    para_stat = mapslices(x->"$(mean(x)) ± $(1.96std(x))", out_mat, dims=1)
    for line in para_stat
        println(line)
    end
end

function run_simple_trait_model_C_F_tuning_RO_ts()
    #Init Ro data
    RO_data = Load_RO_data()

    nsamples = 5000
    nchains = 12

    par_guess = [0.0028160038624083067, 27.33276184313149,
    0.028643648082103606, 0.09336572808750652,
     0.0019727630599060625, 2.2083338895077254,
      -6.817991596443347, 13.184687442117964,
       20.387232607829905, 2.6564297126274803,
        0.434540406001315, 0.06910534451564838,
         0.13628204854679268, 0.1383012394158192]

    x_current = [par_guess for i = 1:nchains]
    P_current = [Post_distri_RO_C_F_CCPH(par_guess,RO_data) for i = 1:nchains]

    println(P_current)

    q_vec = [0.04,24.0,0.033,0.01,0.04,1.0,4.0,7.0,16.0,7.0,
    0.237,0.0845,0.0688,0.146]/40.0

    metropolis_mcmc(x_current,P_current,
    x::Array{Float64,1}->Post_distri_RO_C_F_CCPH(x::Array{Float64,1},RO_data),q_vec,"adaptive_rwm_RO_C_GPP_Ec_20220301";
    n_chains=nchains,n_samples=nsamples,burn_in=7000,sample_freq=5)
end

run_simple_trait_model_C_F_tuning_RO_ts()
#test_adaptive_rwm()
#find_max_adaptive_rwm()


text:
using DashBootstrapComponents, DashHtmlComponents

tooltip = html_div([
    html_p([
        "I wonder what ",
        html_span(
            "floccinaucinihilipilification",
            id="tooltip-target",
            style=Dict("textDecoration" => "underline", "cursor" => "pointer"),
        ),
        " means?",
    ]),
    dbc_tooltip(
        "Noun: rare, " * "the action or habit of estimating something as worthless.",
        target="tooltip-target",
    ),
]);

max_stars_repo_name: harisbal/dash-bootstrap-components

max_issues_repo_name: harisbal/dash-bootstrap-components

max_forks_repo_name: harisbal/dash-bootstrap-components


text:
immutable AlignmentScoring
	match::Int
	mismatch::Int
	gap_start::Int
	gap::Int
end

immutable FOGSAABranch
	m::Int
	x::Int
	y::Int
end

function FOGSAA(X, Y, s::AlignmentScoring)
	const XY = 0x1
	const xY = 0x2
	const Xy = 0x3

	g0 = endof(X) - endof(Y)
	gM = g0 == 0? 0 : s.gap_start

	if g0 < 0
		gM += -g0 * s.gap
		m0 = endof(X)
	else
		gM += g0 * s.gap
		m0 = endof(Y)
	end

	gm = gM + m0 * s.mismatch
	gM += m0 * s.match

	alignments = zeros(UInt8, endof(X) + 1, endof(Y) + 1)
	scores = fill(gm, endof(X) + 1, endof(Y) + 1)
	scores[1, 1] = 0

	Mi = mi = gM - gm + 1
	queue = Vector{Vector{FOGSAABranch}}(Mi)
	queue[Mi] = [FOGSAABranch(1, 1, 1)]

	function enqueue(M, m, x, y, alignment)
		if alignment & 1 != 0
			x = -x
		end

		if alignment & 2 != 0
			y = -y
		end

		branch = FOGSAABranch(m, x, y)

		i = M - gm + 1

		if i > Mi
			Mi = i
		end

		if i < mi || mi == 0
			mi = i
		end

		if isdefined(queue, i)
			j = 1

			while j <= endof(queue[i]) && queue[i][j].m <= m
				j += 1
			end

			insert!(queue[i], j, branch)
		else
			queue[i] = [branch]
		end
	end

	function dequeue()
		x = queue[Mi][end].x
		y = queue[Mi][end].y

		pop!(queue[Mi])

		alignment = 0x0

		if x < 0
			x = -x
			alignment |= 1
		end

		if y < 0
			y = -y
			alignment |= 2
		end

		while !isdefined(queue, Mi) || endof(queue[Mi]) == 0
			if Mi == mi
				Mi = mi = 0
				break
			end

			Mi -= 1
		end

		x, y, alignment
	end

	ox = endof(X) + 1
	oy = endof(Y) + 1

	function print_alignment()
		X! = X[ox:end]
		Y! = Y[oy:end]

		x = ox
		y = oy

		while x != 1 && y != 1
			if alignments[x, y] == xY
				y -= 1

				X! = string("-", X!)
				Y! = string(Y[y], Y!)
			elseif alignments[x, y] == Xy
				x -= 1

				X! = string(X[x], X!)
				Y! = string("-", Y!)
			else
				x -= 1
				y -= 1

				X! = string(X[x], X!)
				Y! = string(Y[y], Y!)
			end
		end

		if x != 1
			t = x - 1
			X! = string(X[1:t], X!)
			Y! = string(repeat("-", t), Y!)
		elseif y != 1
			t = y - 1
			Y! = string(Y[1:t], Y!)
			X! = string(repeat("-", t), X!)
		end

		r = endof(X!) - endof(Y!)

		if r < 0
			X! = string(X!, repeat("-", -r))
		else
			Y! = string(Y!, repeat("-", r))
		end

		println(X!)
		println(Y!)
	end

	while Mi - 1 + gm > scores[ox, oy]
	let
		x, y, alignment = dequeue()

		while x <= endof(X) && y <= endof(Y)
			sorted = 0x0
			score = [(scores[x, y] +
				  (X[x] == Y[y]? s.match : s.mismatch)),
				 (scores[x, y] + s.gap +
				  (alignment == xY? 0 : s.gap_start)),
				 (scores[x, y] + s.gap +
				  (alignment == Xy? 0 : s.gap_start))]
			M = Vector{Int}(3)
			m = Vector{Int}(3)

			if score[XY] > scores[x + 1, y + 1]
				scores[x + 1, y + 1] = M[XY] = score[XY]
			let
				rx = endof(X) - x
				ry = endof(Y) - y
				g = rx - ry
				r = ry

				if g != 0
					M[XY] += s.gap_start

					if g < 0
						g = -g
						r = rx
					end
				end

				M[XY] += g * s.gap
				m[XY] = M[XY] + r * s.mismatch
				M[XY] += r * s.match

				if M[XY] > scores[ox, oy]
					sorted = XY
				end
			end
			end

			if score[xY] > scores[x, y + 1]
				scores[x, y + 1] = M[xY] = score[xY]
			let
				rx = endof(X) + 1 - x
				ry = endof(Y) - y
				g = ry - rx
				r = rx

				if g < 0
					g = -g
					r = ry
					M[xY] += s.gap_start
				end

				M[xY] += g * s.gap
				m[xY] = M[xY] + r * s.mismatch
				M[xY] += r * s.match

				if M[xY] > scores[ox, oy]
					if sorted == 0x0 ||
					  (M[xY] >= M[XY] &&
					   (M[xY] > M[XY] ||
					    m[xY] > m[XY]))
						sorted <<= 2
						sorted |= xY
					else
						sorted |= xY << 2
					end
				end
			end
			end

			if score[Xy] > scores[x + 1, y]
				scores[x + 1, y] = M[Xy] = score[Xy]
			let
				rx = endof(X) - x
				ry = endof(Y) + 1 - y
				g = rx - ry
				r = ry

				if g < 0
					g = -g
					r = rx
					M[Xy] += s.gap_start
				end

				M[Xy] += g * s.gap
				m[Xy] = M[Xy] + r * s.mismatch
				M[Xy] += r * s.match

				if M[Xy] > scores[ox, oy]
					t = sorted >> 2

					if t == 0 ||
					  (M[Xy] >= M[t] &&
					   (M[Xy] > M[t] ||
					    m[Xy] > m[t]))
						u = sorted & 0b11

						if u == 0 ||
						  (M[Xy] >= M[u] &&
						   (M[Xy] > M[u] ||
						    m[Xy] > m[u]))
							sorted <<= 2
							sorted |= Xy
						else
							sorted = u
							t = t << 4 | Xy << 2
							sorted |= t
						end
					else
						sorted |= Xy << 4
					end
				end
			end
			end

			alignment = sorted & 0b11

			if alignment == 0x0
				@goto prune
			end

			a = (sorted & 0b1100) >> 2

			if a != 0 && M[a] > m[alignment]
				enqueue(M[a], m[a], a == xY? x : x + 1,
					a == Xy? y : y + 1, a)

				b = sorted >> 4

				if b != 0 && M[b] > m[a]
					enqueue(M[b], m[b], b == xY? x : x + 1,
						b == Xy? y : y + 1, b)
				end
			end

			x = alignment == xY? x : x + 1
			y = alignment == Xy? y : y + 1
			alignments[x, y] = alignment
		end

		ox = x
		oy = y
@label prune
	end
	end

	print_alignment()
end

max_stars_repo_name: kjenova/FOGSAA

max_issues_repo_name: kjenova/FOGSAA

max_forks_repo_name: kjenova/FOGSAA