text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "UCX"
version = v"1.10.0"
tag = v"1.10.0-rc3"

# Collection of sources required to complete build
sources = [
    ArchiveSource("https://github.com/openucx/ucx/releases/download/v$(tag)/ucx-$(version).tar.gz",
                  "cb9dcc6b53ab141a2dd4cd121e32a35a1537fd73d44164f41f5da956e6c099ec"),
    DirectorySource("./bundled"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/ucx-*

# Apply all our patches
if [ -d $WORKSPACE/srcdir/patches ]; then
for f in $WORKSPACE/srcdir/patches/*.patch; do
    echo "Applying patch ${f}"
    atomic_patch -p1 ${f}
done
fi

update_configure_scripts --reconf

FLAGS=()
FLAGS+=(--prefix=${prefix})
FLAGS+=(--build=${MACHTYPE})
FLAGS+=(--host=${target})
FLAGS+=(--disable-debug)
FLAGS+=(--disable-assertions)
FLAGS+=(--disable-params-check)
FLAGS+=(--disable-static)
FLAGS+=(--disable-profiling)
FLAGS+=(--enable-shared)
FLAGS+=(--enable-mt)
FLAGS+=(--enable-frame-pointer)
FLAGS+=(--enable-cma)
FLAGS+=(--with-rdmacm=${prefix})

if [[ "${target}" != *aarch64* ]]; then
    FLAGS+=(--with-cuda=${prefix}/cuda)
fi

if [[ "${target}" == *x86_64* ]]; then
    FLAGS+=(--with-rocm=${prefix})
fi

./configure ${FLAGS[@]}

# For a bug in `src/uct/sm/cma/Makefile` that I did't have the time to look
# into, we have to build with `V=1`
make -j${nproc} V=1
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    Platform("x86_64", "linux"; libc="glibc"),
    Platform("aarch64", "linux"; libc="glibc"),
    Platform("powerpc64le", "linux"; libc="glibc"),
]

# The products that we will ensure are always built
products = [
    LibraryProduct("libuct", :libuct),
    LibraryProduct("libucm", :libucm),
    LibraryProduct("libucp", :libucp),
    LibraryProduct("libucs", :libucs),
    ExecutableProduct("ucx_info", :ucx_info),
    ExecutableProduct("ucx_perftest", :ucx_perftest),
    ExecutableProduct("ucx_read_profile", :ucx_read_profile),
]

# Dependencies that must be installed before this package can be built
# - librdmacm -> provided through rdma-core, need glibc 2.15
# - libibcm   -> legacy libibverbs
# - knem  -> kernel module
# - xpmem -> kernel module
# - CUDA -> Figure out how version dependent we are
#   - gdrcopy -> kernel module
# - ROCM -> TODO

cuda_version = v"11.2.0"
rocm_version = v"3.7.0"

dependencies = [
    Dependency(PackageSpec(name="CompilerSupportLibraries_jll", uuid="e66e0078-7015-5450-92f7-15fbd957f2ae")),
    Dependency(PackageSpec(name="NUMA_jll", uuid="7f51dc2b-bb24-59f8-b771-bb1490e4195d")),
    Dependency(PackageSpec(name="rdma_core_jll", uuid="69dc3629-5c98-505f-8bcd-225213cebe70")),
    BuildDependency(PackageSpec(name="CUDA_full_jll", version=cuda_version)),
    BuildDependency(PackageSpec(name="hsa_rocr_jll", version=rocm_version))
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies;
               preferred_gcc_version=v"5")

max_stars_repo_name: jebej/Yggdrasil

max_issues_repo_name: jebej/Yggdrasil

max_forks_repo_name: jebej/Yggdrasil


text:
# types and functions describing the topology of the reference element

const TopoIdxType = Int

"""
### ODLCommonTools.ElementTopology

  This type describes the topology of the reference element.  For now, it only
  contains face information, but eventually will needed more info.

  The way to think about element topology is that an element is defined by its
  vertices, numbered 1 through n.  All other parts of the element's topology
  (edges, faces, etc.) can be defined in terms of the vertices.  It is then
  possible to construct the relationship between the higher dimensional
  topological entities directly (ie. edges to faces, etc.), although this is
  not yet implemented.

  Terminology:
    numFaces: the number of faces an element has.  A triangle has 3 faces (the
              edges) and a tet has 4 faces
    numEdges: The number of edges

  TODO: consider making the indices Cints or UInt8s?

  Fields:
    face_verts: Tdim x numFaces array holding the indices of the vertices of
                each face
    edge_verts: The 2 x numEdges array holding the indices of the vertices of
                each face.  Note that this array specifies the orientation of
                the edge, ie. the edge is directed from the  first vertex to
                the second vertex.
                In 2D this is the same as face_verts.

    face_edges: numEdgesPerFace x numFaces array containins, for each face,
                the index of the edges on the face in the element list of
                edges.

  Constructor:
    face_verts: the indices of the vertices that define each face
    edge_verts: the indices of the vertices that define each edge.
                If there are not know, thie argument can be omitted, however
                the face_edges field will not be usable
    topo2:      Only needed in 3D when edge_verts are supplied.  Describes the
                topology of a face of the element
"""

struct ElementTopology{Tdim}
  face_verts::Array{TopoIdxType, 2}
  edge_verts::Array{TopoIdxType, 2}
  face_edges::Array{TopoIdxType, 2}
  face_edges_flipped::Array{Bool, 2}

  function ElementTopology{Tdim}(
                           face_verts::AbstractArray{I1, 2},
                           edge_verts::AbstractArray{I2, 2}=zeros(Int, 0,0);
                           topo2::ElementTopology=ElementTopology{1}() ) where {Tdim, I1 <: Integer, I2 <: Integer}

    # do sanity checks

    # check all vertices are within range
    for i=1:length(face_verts)
      @assert face_verts[i] > 0
      @assert face_verts[i] <= Tdim + 1
    end

    # check faces are distinct
    for i=1:size(face_verts, 2)
      curr_face = sort(face_verts[:, i])
      for j=(i+1):size(face_verts, 2)
        @assert sort(face_verts[:, j]) != curr_face
      end
    end

    if Tdim == 2
      # face_verts == edge_verts implies the edge ordering v1 -> v2, v2 -> v3,
      # v3 -> v1
      # this is not strictly necessary (some people might reverse the direction
      # of the 3rd edge) but seems good enough for now
      edge_verts = face_verts
      face_edges = Array{TopoIdxType}(1, 3)
      face_edges_flipped = Array{Bool}(1, 3)

      for i=1:3
        face_edges[1, i] = i
        face_edges_flipped[1, i] = false
      end
    elseif Tdim == 3
      if size(edge_verts, 1) == 0  # edge verts not known
        face_edges = Array{TopoIdxType}(0, 0)
        face_edges_flipped = Array{Bool}(0, 0)
      else
        face_edges, face_edges_flipped = construct_face_edges(topo2, face_verts, edge_verts)
      end
    end

    return new(face_verts, edge_verts, face_edges, face_edges_flipped)

  end  # end function

  function ElementTopology{Tdim}() where Tdim
    face_verts = zeros(TopoIdxType, 0, 0)
    edge_verts = zeros(face_verts)
    face_edges = zeros(face_verts)
    face_edges_flipped = Array{Bool}(0, 0)
    return new(face_verts, edge_verts, face_edges, face_edges_flipped)
  end
end

"""
  Default constructor that uses Pumi topology
"""
function ElementTopology3()
  face_verts = [1 1 2 1; 2 2 3 3; 3 4 4 4]
  return ElementTopology{3}(face_verts)
end

"""
  Constructs a dummy 2d ElementTopology
"""
function ElementTopology2()
  face_verts = [1 2 3; 2 3 1]
  return ElementTopology{2}(face_verts)
end

#------------------------------------------------------------------------------
# functions to assist the constructor

"""
  This function uses a brute force lookup approach to compute the
  face_edges array.  Throws an exception if the algorithm cannot construct
  the array.

  Inputs:
    topo2: an ElementTopology2 to related the vertices of each face to their
           edges
    face_verts:
    edge_verts:

  Outputs:

"""
function construct_face_edges(topo2::ElementTopology{2},
                              face_verts::AbstractMatrix,
                              edge_verts::AbstractMatrix)
  dim = size(face_verts, 1)
  numFacesPerElement = size(face_verts, 2)
  numEdgesPerElement = size(edge_verts, 2)
  numEdgesPerFace = size(topo2.face_verts, 2)

  # index of the edge of each face in edge_verts
  face_edges = zeros(TopoIdxType, numEdgesPerFace, numFacesPerElement)

  # is the edge of the face oriented the same as the edge of the tet
  face_edges_flipped = Array{Bool}(numEdgesPerFace, numFacesPerElement)

  for i=1:numFacesPerElement
    for j=1:numEdgesPerFace  # loop over edges on this face
      v1 = face_verts[topo2.face_verts[1, j], i]
      v2 = face_verts[topo2.face_verts[2, j], i]

      edge_verts_k = [v1, v2]

      # search edge list
      is_flipped = false
      matched_edgenum = 0
      for k=1:numEdgesPerElement
        matches_edge, is_flipped = compare_edge(edge_verts_k, edge_verts[:, k])
        if matches_edge
          matched_edgenum = k; break
        end
      end

      @assert matched_edgenum != 0  # we should have found the edge somewhere

      face_edges[j, i] = matched_edgenum
      face_edges_flipped[j, i] = is_flipped

    end
  end

  return face_edges, face_edges_flipped
end

function compare_edge(edges1::AbstractVector, edges2::AbstractVector)

  matches = false
  is_flipped = false  # value doesn't matter if matches is false
  if edges1[1] == edges2[1] && edges1[2] == edges2[2]
    matches = true
    is_flipped = false
  elseif edges1[1] == edges2[2] && edges1[2] == edges2[1]
    matches = true
    is_flipped = true
  else
    matches = false
    is_flipped = false  # arbitrary
  end

  return matches, is_flipped
end

ext: jl


text:
"""
    self_efficacy2scores(echo::DataFrame)::DataFrame

Return the scores for the echo questionnaire, which is about self-efficacy.
None of the items appear to be reversed, so this method just returns the sum.

This function assumes questions [1, 14] which holds for all questionnaires in from 2018 to the time of writing (2021).
"""
function self_efficacy2scores(echo::DataFrame)::DataFrame
    df = disallowmissing(echo)
    questions = [Symbol("v$i") for i in 1:14]
    cols = [df[:, col] for col in questions]
    df.score = [sum(t) for t in zip(cols...)]
    select!(df, Not(questions))
end


text:
#== # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
#
#   Pre-defined formats.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # ==#

export tf_unicode, tf_ascii_dots, tf_ascii_rounded, tf_borderless, tf_compact,
       tf_markdown, tf_mysql, tf_simple, tf_unicode_rounded, tf_matrix

const tf_unicode = TextFormat()

const tf_ascii_dots = TextFormat(
    up_right_corner     = '.',
    up_left_corner      = '.',
    bottom_left_corner  = ':',
    bottom_right_corner = ':',
    up_intersection     = '.',
    left_intersection   = ':',
    right_intersection  = ':',
    middle_intersection = ':',
    bottom_intersection = ':',
    column              = ':',
    row                 = '.'
)

const tf_ascii_rounded = TextFormat(
    up_right_corner     = '.',
    up_left_corner      = '.',
    bottom_left_corner  = ''',
    bottom_right_corner = ''',
    up_intersection     = '.',
    left_intersection   = ':',
    right_intersection  = ':',
    middle_intersection = '+',
    bottom_intersection = ''',
    column              = '|',
    row                 = '-'
)

const tf_borderless = TextFormat(
    up_right_corner     = ' ',
    up_left_corner      = ' ',
    bottom_left_corner  = ' ',
    bottom_right_corner = ' ',
    up_intersection     = ' ',
    left_intersection   = ' ',
    right_intersection  = ' ',
    middle_intersection = ' ',
    bottom_intersection = ' ',
    column              = ' ',
    row                 = ' ',
    hlines              = [:header]
)

const tf_compact = TextFormat(
    up_right_corner     = ' ',
    up_left_corner      = ' ',
    bottom_left_corner  = ' ',
    bottom_right_corner = ' ',
    up_intersection     = ' ',
    left_intersection   = ' ',
    right_intersection  = ' ',
    middle_intersection = ' ',
    bottom_intersection  = ' ',
    column              = ' ',
    row                 = '-'
   )

const tf_markdown = TextFormat(
    left_intersection   = '|',
    right_intersection  = '|',
    middle_intersection = '|',
    column              = '|',
    row                 = '-',
    hlines              = [:header]
)

const tf_mysql = TextFormat(
    up_right_corner     = '+',
    up_left_corner      = '+',
    bottom_left_corner  = '+',
    bottom_right_corner = '+',
    up_intersection     = '+',
    left_intersection   = '+',
    right_intersection  = '+',
    middle_intersection = '+',
    bottom_intersection = '+',
    column              = '|',
    row                 = '-'
   )

const tf_simple = TextFormat(
    up_right_corner     = '=',
    up_left_corner      = '=',
    bottom_left_corner  = '=',
    bottom_right_corner = '=',
    up_intersection     = ' ',
    left_intersection   = '=',
    right_intersection  = '=',
    middle_intersection = ' ',
    bottom_intersection  = ' ',
    column              = ' ',
    row                 = '='
)

const tf_unicode_rounded = TextFormat(
    up_right_corner     = '╮',
    up_left_corner      = '╭',
    bottom_left_corner  = '╰',
    bottom_right_corner = '╯',
    up_intersection     = '┬',
    left_intersection   = '├',
    right_intersection  = '┤',
    middle_intersection = '┼',
    bottom_intersection = '┴',
    column              = '│',
    row                 = '─'
)

const tf_matrix = TextFormat(
    left_intersection   = '│',
    right_intersection  = '│',
    row                 = ' ',
    vlines              = [:begin,:end],
    hlines              = [:begin,:end]
)


text:
using Test, Infiltrator, REPL

foo = 3
function f(x)
    x = 2
    y = [1, 2, 3]
    @infiltrate
    return x .+ y
end

function g(x)
    x *= 12
    @infiltrate x > 12
    return x
end

if Sys.isunix() && VERSION >= v"1.1.0"
    using TerminalRegressionTests

    function run_terminal_test(func, result, commands, validation)
        # TerminalRegressionTests.automated_test(joinpath(@__DIR__, validation), commands) do emuterm
        TerminalRegressionTests.create_automated_test(joinpath(@__DIR__, validation), commands) do emuterm
            repl = REPL.LineEditREPL(emuterm, true)
            repl.interface = REPL.setup_interface(repl)
            repl.specialdisplay = REPL.REPLDisplay(repl)

            Infiltrator.TEST_TERMINAL_REF[] = repl.t
            Infiltrator.TEST_NOSTACK[] = true
            Infiltrator.TEST_REPL_REF[] = repl
            @test func() == result
            Infiltrator.TEST_TERMINAL_REF[] = nothing
            Infiltrator.TEST_NOSTACK[] = false
            Infiltrator.TEST_REPL_REF[] = nothing
        end
    end

    run_terminal_test(() -> f(3), [3, 4, 5],
                      ["?\n", "@trace\n", "@locals\n", "x.*y\n", "foo\n", "0//0\n", "@stop\n", "\x4"],
                      "Julia_f_$(VERSION.major).$(VERSION.minor).multiout")

    @test f(3) == [3, 4, 5] # `@stop`ped `@infiltrate` should not open a prompt

    Infiltrator.clear_stop()

    run_terminal_test(() -> f(3), [3, 4, 5],
                      ["@locals\n", "\x4"],
                      "Julia_f2_$(VERSION.major).$(VERSION.minor).multiout")

    @test g(1) == 12 # conditional @infiltrate should not open a prompt

    run_terminal_test(() -> g(2), 24,
                      ["?\n", "@trace\n", "@locals\n", "x\n", "\x4"],
                      "Julia_g_$(VERSION.major).$(VERSION.minor).multiout")
else
    @warn "Skipping UI tests on non unix systems"
end


text:
index_changed = joinpath(@__DIR__, "fixtures", "index_changed.html")
cp(joinpath(@__DIR__, "fixtures", "index.html"), index_changed, force=true)
OutdatedWarning.add_old_docs_notice(index_changed)
output = read(index_changed, String)

@test occursin("data-outdated-warner", output)
# The following is sensitive to how Gumbo.jl writes HTML, but let's hope that doesn't change often:
@test replace(output, "\r\n" => "\n") == replace(read(joinpath(@__DIR__, "fixtures", "index_after.html"), String), "\r\n" => "\n")

rm(index_changed)

transient_path = joinpath(@__DIR__, "fixtures", "transient")
cp(joinpath(@__DIR__, "fixtures", "pre"), transient_path, force=true)
OutdatedWarning.generate(transient_path)

for (root, _, files) in walkdir(transient_path)
    for file in files
        content = read(joinpath(root, file), String)
        expected = read(joinpath(replace(root, "transient" => "post"), file), String)
        @test replace(content, "\r\n" => "\n") == replace(expected, "\r\n" => "\n")
    end
end

rm(transient_path, recursive=true)


text:
#
# Correctness Tests
#

using Base.Test
using SimilarityMetrics

my_tests = ["test/cosine.jl",
            "test/expl1.jl",
            "test/expl2.jl",
            "test/jaccard.jl",
            "test/matching.jl",
            "test/overlap.jl",
            "test/pearson.jl"]

@printf "Running tests:\n"

for my_test in my_tests
    @printf " * %s\n" my_test
    include(my_test)
end


text:
## PART 2

# Define the error functions

function σf1(v̄̂,σ̂)
    ddvs = zeros(5);
    ddvs[1] = 1 / v̄̂[2]^2;
    ddvs[2] = v̄̂[1]^2 / v̄̂[2]^4;
    return sqrt( dot( ddvs, σ̂.^2) );
end

function σf2(v̄̂,σ̂)
    ddvs = zeros(5);
    ddvs[3] = exp(v̄̂[3] - v̄̂[4])^2;
    ddvs[4] = exp(v̄̂[3] - v̄̂[4])^2;
    return sqrt( dot( ddvs, σ̂.^2) );
end

function σf3(v̄̂,σ̂)
    lv5s = log(v̄̂[5])^2;
    ddvs = zeros(5);
    ddvs[1] = lv5s / v̄̂[2]^2;
    ddvs[2] = lv5s * v̄̂[1]^2 / v̄̂[2]^4;
    ddvs[3] = lv5s / v̄̂[4]^2;
    ddvs[4] = lv5s * v̄̂[3]^2 / v̄̂[4]^4;
    ddvs[5] = ( ( v̄̂[1]/v̄̂[2] + v̄̂[3]/v̄̂[4] ) / v̄̂[5] )^2;
    return sqrt( dot( ddvs, σ̂.^2) );
end

println("Calculating σ̂f1, σ̂f2, σ̂f3 for sample sizes N1, N2");
σ̂f1N1 = σf1(v̄̂,σ̂1);
σ̂f2N1 = σf2(v̄̂,σ̂1);
σ̂f3N1 = σf3(v̄̂,σ̂1);

σ̂f1N2 = σf1(v̄̂,σ̂2);
σ̂f2N2 = σf2(v̄̂,σ̂2);
σ̂f3N2 = σf3(v̄̂,σ̂2);

println("\n\tNaive standard deviations for sample size 1,000");
println("\tσ̂f1N1 = $σ̂f1N1");
println("\tσ̂f2N1 = $σ̂f2N1");
println("\tσ̂f3N1 = $σ̂f3N1");

println("\n\tNaive standard deviations for sample size 10,000");
println("\tσ̂f1N2 = $σ̂f1N2");
println("\tσ̂f2N2 = $σ̂f2N2");
println("\tσ̂f3N2 = $σ̂f3N2");

println("\n\tTrue standard deviations for sample size 1,000");
println("\tσ̂truef1N1 = $σ̂truef1N1");
println("\tσ̂truef2N1 = $σ̂truef2N1");
println("\tσ̂truef3N1 = $σ̂truef3N1");

println("\n\tTrue standard deviations for sample size 10,000");
println("\tσ̂truef1N2 = $σ̂truef1N2");
println("\tσ̂truef2N2 = $σ̂truef2N2");
println("\tσ̂truef3N2 = $σ̂truef3N2");


text:
"""
    main_mode1
A main file to launch PSCOPF in MODE 1.

Parameters:
    instance_path : path to input data directory
    output_path : path to output directory
"""

using Dates

root_path = dirname(@__DIR__)
push!(LOAD_PATH, root_path);
cd(root_path)
include(joinpath(root_path, "src", "PSCOPF.jl"));

#########################
# INPUT & PARAMS
#########################

# instance_path is the path to the input data directory containing :
# pscopf_units : the description of the units
# pscopf_init : the description of the initial state of units (just before the first interest timepoit T)
# pscopf_limits : the flow capacity of each branch
# pscopf_ptdf : the ptdf coefficients per (branch, bus_id)
# pscopf_uncertainties : the nodal injections (for each bus and each limitable)
instance_path = ( length(ARGS) > 0 ? ARGS[1] :
                    joinpath(@__DIR__, "..", "data", "2buses_usecase", "instance") )

# output_path is the path where output files will be write_commitment_schedule
#NOTE: all files in output_path, except those starting with pscopf_, will be deleted
output_path = length(ARGS) > 1 ? ARGS[2] : instance_path

#########################
# EXECUTION
#########################

# load Data
network = PSCOPF.Data.pscopfdata2network(instance_path)
uncertainties = PSCOPF.PSCOPFio.read_uncertainties(instance_path)
generators_init_state = PSCOPF.PSCOPFio.read_initial_state(instance_path)

# Launch mode 1
mode = PSCOPF.PSCOPF_MODE_1

ts1 = Dates.DateTime("2015-01-01T11:00:00") #11h
TS = PSCOPF.create_target_timepoints(ts1) #T: 11h, 11h15, 11h30, 11h45
ECH = PSCOPF.generate_ech(network, TS, mode) #ech: -4h, -1h, -30mins, -15mins, 0h

sequence = PSCOPF.generate_sequence(network, TS, ECH, mode)

PSCOPF.rm_non_prefixed(output_path, "pscopf_")
exec_context = PSCOPF.PSCOPFContext(network, TS, mode,
                                    generators_init_state,
                                    uncertainties, nothing,
                                    output_path)

PSCOPF.run!(exec_context, sequence)

max_stars_repo_name: klorel/pscopf

max_issues_repo_name: klorel/pscopf

max_forks_repo_name: klorel/pscopf


text:
@testset "$TEST $G" begin

@test length(collect(edges(G()))) == 0

ga = G(10,20; seed=1)
gb = G(10,20; seed=1)

@test length(collect(edges(ga))) == 20
@test collect(edges(ga)) == collect(edges(gb))
for e in edges(ga)
    @test has_edge(ga, e)
end
@test length(unique(collect(edges(ga)))) == 20

ga = erdos_renyi(10,20, G; seed=1)
gb = erdos_renyi(10,20, DG; seed=1)

# TODO implement length(edges(g)) somehow
# @test length(collect(edges(ga))) == length(edges(ga)) == ne(ga)
# @test length(collect(edges(gb))) == length(edges(gb)) == ne(gb)
@test length(collect(edges(ga))) == ne(ga)
@test length(collect(edges(gb))) == ne(gb)

for e in edges(ga)
    @test has_edge(ga, e)
end
@test length(unique(collect(edges(ga)))) == 20

ga = G(10)
add_edge!(ga, 3, 2)
add_edge!(ga, 3, 10)
add_edge!(ga, 5, 10)
add_edge!(ga, 10, 3)

eit = edges(ga)
# TODO?
# @test Base.IteratorSize(typeof(eit)) == Base.HasLength()

i = 0
for e in eit
    @test src(e) <= dst(e)
    i += 1
end
@test i == ne(ga)

@test [e for e in eit] == [E(2, 3), E(3, 10), E(5,10)]

ga = DG(10)
add_edge!(ga, 3, 2)
add_edge!(ga, 3, 10)
add_edge!(ga, 5, 10)
add_edge!(ga, 10, 3)

eit = edges(ga)

@test [e for e in eit] == [E(3, 2), E(3, 10), E(5,10), E(10,3)]
@test collect(eit) == [e for e in eit]
# g = CompleteGraph(10)
# @test collect(edges(g,1:3)) == collect(edges(g,[1:3;]))
# @test length(collect(edges(g, 1:3))) == 3
# @test length(collect(edges(g, 1:0))) == 0
#
#
# g = CompleteDiGraph(10)
# @test collect(edges(g,1:3)) == collect(edges(g,[1:3;]))
# @test length(collect(edges(g, 1:3))) == 6
# @test length(collect(edges(g, 1:0))) == 0

end # testset


text:
module SA

##export simulated_annealing, random_search, grid_search

using Distributions
using SharedArrays

function grid_search(M, lower, upper, objective::Function)
    max_obj = objective(upper)
    best = lower
    K = length(M)

    ranges = map(i -> lower[i]:(upper[i]-lower[i])/M[i]:upper[i], 1:K)

    @inbounds for el in Base.product(ranges...)
        new_obj = objective(el)
        if new_obj > max_obj
            best, max_obj = el, new_obj
        end
    end

    best, max_obj
end

function simulated_annealing(M::Int, init::T,
                             kernel::Function,
                             objective::Function,
                             temperature::Function = t -> 1/log(t)) where T
    x = Vector{T}(undef,M+1)
    x[1] = init

    old_obj = objective(x[1])
    max_obj = old_obj
    accept = 0
    best_x = x[1]

    obj = Vector{T}(undef,M+1)
    obj[1] = old_obj

    @inbounds for t in 1:M
        proposal = kernel(x[t], temperature(t+1))
        new_obj = objective(proposal)

        if log(rand()) <= (new_obj - old_obj) / temperature(t+1)
            x[t+1] = proposal
            old_obj = new_obj

            if max_obj < new_obj
                max_obj = new_obj
                best_x = proposal
                accept += 1
            end
        else
            x[t+1] = x[t]
        end

        obj[t+1] = old_obj
    end

    max_obj, best_x, accept/(M+1)
    ##max_obj, best_x, obj[end], x[end], accept/(M+1)
end

function random_search(M::Int, init::T, kernel::Function, objective::Function) where T
    max_obj = objective(init)
    x::T = init
    accept = 0

    for t in 1:M
        proposal = kernel(x)
        new_obj = objective(proposal)

        if new_obj > max_obj
            x = proposal
            max_obj = new_obj
            accept += 1
        end
    end

    max_obj, x, accept/(M+1)
end

function test_kernel(x, temp)
    rand(Normal(0, 1), 2)
end

function test_kernel2(x, temp)
    rand(Uniform(-4,4), 2)
end

function test_obj(x)
    pdf(MultivariateNormal([1.0, -0.5], eye(2)), x)
end

function rosenbrock(x)
    -((1 - x[1])^2 + 100*(x[2] - x[1]^2)^2)
end

function himmelblau(x)
  -((x[1]^2 + x[2] - 11)^2 + (x[1] + x[2]^2 - 7)^2)
end

end


text:
# Benchmark the calculation of a distance in a PDB file
# The distance is the closest distance between any atoms of residues 50 and 60
#   of chain A in 1AKE

using Bio.Structure

pdb_filepath = "pdbs/1AKE.pdb"
struc = read(pdb_filepath, PDB)

# Run to JIT compile
distance(struc['A'][50], struc['A'][60])

elapsed = @elapsed distance(struc['A'][50], struc['A'][60])

println(elapsed)

max_stars_repo_name: chemlove/pdb-benchmarks

max_issues_repo_name: chemlove/pdb-benchmarks

max_forks_repo_name: chemlove/pdb-benchmarks


text:
 @testset "Direct Ellipse Fit" begin
    @info "Test: Direct Ellipse Fit"

    @testset "Numerical" begin
        truth = SVector(4.053658735425687e-6,
                        -1.3487952962436518e-7,
                        1.8534714216409155e-5,
                        -0.001965777847001004,
                        -0.009242570343232327,
                         0.9999553541288334)

        𝛉 = fit_ellipse(data, DirectEllipseFit())
        @test all(𝛉 .≈ truth)
    end

end


text:
module PigeonViz

using Distributed
using RobotOS
using StaticArrays
using Interpolations
using Contour
using JLD2
using RobotOS
@rosimport std_msgs.msg: ColorRGBA
@rosimport geometry_msgs.msg: Point
@rosimport visualization_msgs.msg: Marker
@rosimport safe_traffic_weaving.msg: XYThV
@rosimport auto_messages.msg: from_autobox
rostypegen(@__MODULE__)
import .std_msgs.msg: ColorRGBA
import .geometry_msgs.msg: Point
import .visualization_msgs.msg: Marker
import .safe_traffic_weaving.msg: XYThV
import .auto_messages.msg: from_autobox
import StaticArrays: SUnitRange
import Interpolations: GriddedInterpolation

@inline mod2piF(x::T) where {T} = mod(x, 2*T(pi))
@inline function adiff(x, y)
    d = mod2piF(x - y)
    ifelse(d <= π, d, d - 2*oftype(d, π))
end

### HJICache
struct BicycleState{T} <: FieldVector{6,T}
    E::T    # world frame "x" position of CM
    N::T    # world frame "y" position of CM
    ψ::T    # world frame heading of vehicle
    Ux::T   # body frame longitudinal speed
    Uy::T   # body frame lateral speed
    r::T    # yaw rate (dψ/dt)
end
struct SimpleCarState{T} <: FieldVector{4,T}
    E::T    # world frame "x" position of CM
    N::T    # world frame "y" position of CM
    ψ::T    # world frame heading of vehicle
    V::T    # speed of vehicle (velocity assumed to be in the heading direction)
end
struct HJIRelativeState{T} <: FieldVector{7,T}
    ΔE::T
    ΔN::T
    Δψ::T
    Ux::T
    Uy::T
    V::T
    r::T
end
function HJIRelativeState(us::BicycleState, them::SimpleCarState)
    cψ, sψ = sincos(us.ψ)
    ΔE, ΔN = @SMatrix([cψ sψ; -sψ cψ])*SVector(them.E - us.E, them.N - us.N)
    HJIRelativeState(ΔE, ΔN, adiff(them.ψ, us.ψ), us.Ux, us.Uy, them.V, us.r)
end
struct HJICache
    grid_knots::NTuple{7,Vector{Float32}}
    V ::GriddedInterpolation{Float32,7,Float32,Gridded{Linear},NTuple{7,Vector{Float32}}}
    ∇V::GriddedInterpolation{SVector{7,Float32},7,SVector{7,Float32},Gridded{Linear},NTuple{7,Vector{Float32}}}
end
function HJICache(fname::String)
    @load fname grid_knots V_raw ∇V_raw
    V  = interpolate(Float32, Float32, grid_knots, V_raw, Gridded(Linear()))
    ∇V = interpolate(Float32, SVector{7,Float32}, grid_knots, Array(reinterpret(SVector{7,Float32}, ∇V_raw)), Gridded(Linear()))
    HJICache(grid_knots, V, ∇V)
end
function Base.getindex(cache::HJICache, x::HJIRelativeState{T}) where {T}
    if all(cache.grid_knots[i][1] <= x[i] <= cache.grid_knots[i][end] for i in 1:length(cache.grid_knots))
        (V=cache.V(x[1], x[2], x[3], x[4], x[5], x[6], x[7]), ∇V=cache.∇V(x[1], x[2], x[3], x[4], x[5], x[6], x[7]))    # avoid splatting penalty
    else
        (V=T(Inf), ∇V=zeros(SVector{7,T}))
    end
end

const cache = HJICache(joinpath(@__DIR__, "../deps/BicycleCAvoid.jld2"))

### /HJI_values
function initialize_HJI_values_marker!(m)
    m.header.frame_id = "x1"
    m.type = Marker[:TRIANGLE_LIST]
    m.pose.orientation.w = 1.0
    m.scale.x = 1.0
    m.scale.y = 1.0
    m.scale.z = 1.0
    m.frame_locked = true
    X = cache.grid_knots[1]
    Y = cache.grid_knots[2]
    for i in 1:length(X)-1, j in 1:length(Y)-1
        push!(m.points, Point(X[i]  , Y[j]  , -0.1)); push!(m.colors, ColorRGBA(1.0, 1.0, 1.0, 1.0))
        push!(m.points, Point(X[i+1], Y[j]  , -0.1)); push!(m.colors, ColorRGBA(1.0, 1.0, 1.0, 1.0))
        push!(m.points, Point(X[i+1], Y[j+1], -0.1)); push!(m.colors, ColorRGBA(1.0, 1.0, 1.0, 1.0))
        push!(m.points, Point(X[i]  , Y[j]  , -0.1)); push!(m.colors, ColorRGBA(1.0, 1.0, 1.0, 1.0))
        push!(m.points, Point(X[i]  , Y[j+1], -0.1)); push!(m.colors, ColorRGBA(1.0, 1.0, 1.0, 1.0))
        push!(m.points, Point(X[i+1], Y[j+1], -0.1)); push!(m.colors, ColorRGBA(1.0, 1.0, 1.0, 1.0))
    end
end
function update_HJI_values_marker!(m, q)
    X = cache.grid_knots[1]
    Y = cache.grid_knots[2]
    ct = 1
    for i in 1:length(X)-1, j in 1:length(Y)-1
        for (x,y) in ((X[i]  , Y[j]  ),
                      (X[i+1], Y[j]  ),
                      (X[i+1], Y[j+1]),
                      (X[i]  , Y[j]  ),
                      (X[i]  , Y[j+1]),
                      (X[i+1], Y[j+1]))
            q = HJIRelativeState{Float32}(x, y, q[3], q[4], q[5], q[6], q[7])
            color = m.colors[ct]
            color.r, color.g, color.b = value_to_RGB(cache[q].V)
            ct += 1
        end
    end
end
function value_to_RGB(V, V_lo=-3.0, V_hi=20.0, C_lo = SVector(1.0, 0.5, 0.0), C_hi = SVector(0.0, 0.5, 1.0))
    x = clamp(ifelse(V < 0, 0.5*(V_lo - V)/V_lo, 0.5 + 0.5*V/V_hi), 0.0, 1.0)
    (1 - x)*C_lo + x*C_hi
end

const HJI_values_marker = Marker()
initialize_HJI_values_marker!(HJI_values_marker)

### /HJI_contour
function initialize_HJI_contour_marker!(m)
    m.header.frame_id = "x1"
    m.type = Marker[:LINE_STRIP]
    m.pose.orientation.w = 1.0
    m.scale.x = 0.1
    m.scale.y = 0.1
    m.scale.z = 0.1
    m.frame_locked = true
    m.color = ColorRGBA(1.0, 1.0, 1.0, 1.0)
end
function update_HJI_contour_marker!(m, q)
    X = cache.grid_knots[1]
    Y = cache.grid_knots[2]
    c = contour(X, Y, [cache[HJIRelativeState{Float32}(x, y, q[3], q[4], q[5], q[6], q[7])].V for x in X, y in Y], Float32(0))
    if isempty(c.lines[1].vertices)
        resize!(m.points, 0)
    else
        m.points = [Point(x, y, -0.05) for (x,y) in c.lines[1].vertices]
    end
end

const HJI_contour_marker = Marker()
initialize_HJI_contour_marker!(HJI_contour_marker)

### /xbox_car/xythv
const other_car_state = fill(zeros(SimpleCarState{Float64}))
function other_car_callback(msg::XYThV)
    other_car_state[] = SimpleCarState(msg.x, msg.y, msg.th - pi/2, msg.v)
end

### /from_autobox
function from_autobox_callback(msg::from_autobox, HJI_values_pub, HJI_contour_pub)
    try
        current_state = BicycleState(msg.E_m, msg.N_m, msg.psi_rad, msg.ux_mps, msg.uy_mps, msg.r_radps)
        relative_state = HJIRelativeState(current_state, other_car_state[])
        update_HJI_values_marker!(HJI_values_marker, relative_state)
        update_HJI_contour_marker!(HJI_contour_marker, relative_state)
        publish(HJI_values_pub, HJI_values_marker)
        publish(HJI_contour_pub, HJI_contour_marker)
    catch err
        RobotOS.logwarn("Error: $err\n$(stacktrace(catch_backtrace()))")
    end
end

### ROS node init
function start_ROS_node()
    init_node("pigeon_viz", anonymous=false)
    other_car = RobotOS.get_param("human", "/xbox_car")
    HJI_values_pub = Publisher{Marker}("/HJI_values", queue_size=1)
    HJI_contour_pub = Publisher{Marker}("/HJI_contour", queue_size=1)
    Subscriber{from_autobox}("/from_autobox", from_autobox_callback, (HJI_values_pub, HJI_contour_pub), queue_size=1)
    Subscriber{XYThV}("$(other_car)/xythv", other_car_callback, queue_size=1)
    @spawn spin()
end

end

max_stars_repo_name: Michelle-NYX/CollisionAvoidance

max_issues_repo_name: Michelle-NYX/CollisionAvoidance

max_forks_repo_name: Michelle-NYX/CollisionAvoidance


text:
module FM

#
# Single cell dynamics of the fuel moisture model including the Kalman predict/update mechanisms.
#
#

type FMModel
    k::Int
    m_ext::Vector{Float64}
    P::Matrix{Float64}
    Tk::Vector{Float64}
    Trk::Float64
    r0::Float64
    rk::Float64
    S::Float64
    latlon::Tuple

    function FMModel(latlon, k, m0, P0, Tk)
        # construct initial state
        m_ext = zeros(2*k+3)
        m_ext[1:k] = m0

        # initialize with defaults
        new(k, m_ext, copy(P0), Tk * 3600, 14.0 * 3600, 5.0, 8.0, 2.5, latlon)
    end

end

function advance_model(c::FMModel, Ed::Float64, Ew::Float64,
                       r::Float64, dt::Float64, Q::Matrix{Float64})

    # acquire local variables
    Tk = c.Tk
    k = c.k
    m_ext = c.m_ext
    m = m_ext[1:k]
    dTk = m_ext[k+1:2*k]
    dE, dS, dTrk = m_ext[2*k+1:2*k+3]

    # add assimilation deltas to moisture equilibria
    Ed = max(0, Ed + dE)
    Ew = max(0, Ew + dE)

    # compute equilibria
    equi = copy(m)
    rlag = zeros(k)
    model_ids = zeros(Int32, k)

    # switch between models, first branch is rain model, rest are for no-rain conditions
    if r > c.r0
        equi[:] = c.S + dS
        model_ids[:] = 3
        rlag[:] = 1.0 / (c.Trk + dTrk) * (1.0 - exp(- (r - c.r0) / c.rk))
    else
        model_ids[:] = 4
        for i in [1:k]
            # drying model
            if equi[i] > Ed
                model_ids[i] = 1
                equi[i] = Ed
            # wetting model
            elseif equi[i] < Ew
                model_ids[i] = 2
                equi[i] = Ew
            end
        end
        rlag[:] = 1.0 ./ (Tk + dTk)
    end

    # select integration method according to actual change
    # in each fuel
    change = dt * rlag
    m_new = zeros(k)
    for i in 1:k
        if change[i] < 0.01
            m_new[i] = m[i] + (equi[i] - m[i]) * (1.0 - exp(-change[i]))
        else
            m_new[i] = m[i] + (equi[i] - m[i]) * change[i] * (1.0 - 0.5 * change[i])
        end
    end

    # compute Jacobian and update local covariance matrix
    J = zeros((2*k+3,2*k+3))
    for i in 1:k

        if change[i] < 0.01
            J[i,i] = model_ids[i] < 4 ? exp(-change[i]) : 1.0
            dmi_dchng = (equi[i] - m[i]) * exp(-change[i])
            dmi_dequi = (1.0 - exp(-change[i]))
        else
            J[i,i] = model_ids[i] < 4 ? 1.0 - change[i] * (1.0 - 0.5 * change[i]) : 1.0
            dmi_dchng = (equi[i] - m[i]) * (1.0 - change[i])
            dmi_dequi = change[i] * (1.0 - 0.5 * change[i])
        end

        if r < c.r0 # drying/wetting model or dead zone
            if model_ids[i] < 4 # if not dead zone
                J[i,k+i] = dmi_dchng * (-dt) * (Tk[i] + dTk[i])^(-2)
                J[i,2*k+1] = dmi_dequi
            end
        else # rain model
            J[i,2*k+2] = dmi_dequi
            J[i,2*k+3] = dmi_dchng * dt * (exp(-(r - c.r0) / c.rk) - 1.0) * (c.Trk + dTrk)^(-2)
        end

        J[k+i,k+i] = 1.0

    end
    J[2*k+1,2*k+1] = 1.0
    J[2*k+2,2*k+2] = 1.0
    J[2*k+3,2*k+3] = 1.0

    # update covariance (using computed Jacobian) and state
    c.P = J*c.P*J' + Q
    c.m_ext[1:k] = m_new
end

function kalman_update(c::FMModel, O::Vector{Float64},
                       V::Matrix{Float64}, fuel_types::Vector{Int})
    P = c.P
    Nf = size(fuel_types,1)
    H = zeros((Nf,2*c.k+3))

    # construct observation operator for the given fuel types
    ndx = 1
    for i in fuel_types
        H[ndx,i] = 1.0
        ndx += 1
    end

    # Kalman update matrices
    K = P * H' / (H * P * H' + V)

    # update the state
    c.m_ext = c.m_ext + K * (O - c.m_ext[fuel_types])
    c.P -= K * H * P

    # return the kalman update
    return K[fuel_types]

end

end


text:
module VortexHelperVortexSwitchGate
using ..Ahorn, Maple

function gateFinalizer(entity)
    x, y = Ahorn.position(entity)

    width = Int(get(entity.data, "width", 8))
    height = Int(get(entity.data, "height", 8))

    entity.data["nodes"] = [(x + width, y)]
end

@mapdef Entity "VortexHelper/VortexSwitchGate" VortexSwitchGate(
			x::Integer, y::Integer,
			width::Integer=8, height::Integer=8,
			sprite::String="block",
			behavior::String="crush",
			persistent::Bool=false,
			crushDuration::Number=0.75,
			nodes::Array{Tuple{Integer, Integer}, 1}=Tuple{Integer, Integer}[])

const textures = ["block", "mirror", "temple", "stars"]
const behaviors = ["crush", "shatter"]

const placements = Ahorn.PlacementDict(
    "Custom Switch Gate ($(uppercasefirst(behavior))) (Vortex Helper)" => Ahorn.EntityPlacement(
        VortexSwitchGate,
		"rectangle",
		Dict{String, Any}(
			"behavior" => behavior
        ),
        gateFinalizer
    ) for behavior in behaviors
)

Ahorn.editingOptions(entity::VortexSwitchGate) = Dict{String, Any}(
    "sprite" => textures,
	"behavior" => behaviors
)

Ahorn.nodeLimits(entity::VortexSwitchGate) = 1, 1
Ahorn.minimumSize(entity::VortexSwitchGate) = 16, 16
Ahorn.resizable(entity::VortexSwitchGate) = true, true

function Ahorn.selection(entity::VortexSwitchGate)
	behavior = get(entity.data, "behavior", "crush")
    x, y = Ahorn.position(entity)
    stopX, stopY = Int.(entity.data["nodes"][1])

    width = Int(get(entity.data, "width", 8))
    height = Int(get(entity.data, "height", 8))

    return (behavior == "crush") ? [Ahorn.Rectangle(x, y, width, height), Ahorn.Rectangle(stopX, stopY, width, height)] : Ahorn.Rectangle(x, y, width, height)
end

iconResource = "objects/switchgate/icon00"

function renderGateSwitch(ctx::Ahorn.Cairo.CairoContext, x::Number, y::Number, width::Number, height::Number, sprite::String)
    iconSprite = Ahorn.getSprite(iconResource, "Gameplay")

    tilesWidth = div(width, 8)
    tilesHeight = div(height, 8)

    frame = "objects/switchgate/$sprite"

    for i in 2:tilesWidth - 1
        Ahorn.drawImage(ctx, frame, x + (i - 1) * 8, y, 8, 0, 8, 8)
        Ahorn.drawImage(ctx, frame, x + (i - 1) * 8, y + height - 8, 8, 16, 8, 8)
    end

    for i in 2:tilesHeight - 1
        Ahorn.drawImage(ctx, frame, x, y + (i - 1) * 8, 0, 8, 8, 8)
        Ahorn.drawImage(ctx, frame, x + width - 8, y + (i - 1) * 8, 16, 8, 8, 8)
    end

    for i in 2:tilesWidth - 1, j in 2:tilesHeight - 1
        Ahorn.drawImage(ctx, frame, x + (i - 1) * 8, y + (j - 1) * 8, 8, 8, 8, 8)
    end

    Ahorn.drawImage(ctx, frame, x, y, 0, 0, 8, 8)
    Ahorn.drawImage(ctx, frame, x + width - 8, y, 16, 0, 8, 8)
    Ahorn.drawImage(ctx, frame, x, y + height - 8, 0, 16, 8, 8)
    Ahorn.drawImage(ctx, frame, x + width - 8, y + height - 8, 16, 16, 8, 8)

    Ahorn.drawImage(ctx, iconSprite, x + div(width - iconSprite.width, 2), y + div(height - iconSprite.height, 2))
end

function Ahorn.renderSelectedAbs(ctx::Ahorn.Cairo.CairoContext, entity::VortexSwitchGate, room::Maple.Room)
	behavior = get(entity.data, "behavior", "crush")

	if behavior == "crush"
		sprite = get(entity.data, "sprite", "block")
		startX, startY = Int(entity.data["x"]), Int(entity.data["y"])
		stopX, stopY = Int.(entity.data["nodes"][1])

		width = Int(get(entity.data, "width", 32))
		height = Int(get(entity.data, "height", 32))

		renderGateSwitch(ctx, stopX, stopY, width, height, sprite)
		Ahorn.drawArrow(ctx, startX + width / 2, startY + height / 2, stopX + width / 2, stopY + height / 2, Ahorn.colors.selection_selected_fc, headLength=6)
	end
end

function Ahorn.renderAbs(ctx::Ahorn.Cairo.CairoContext, entity::VortexSwitchGate, room::Maple.Room)
    sprite = get(entity.data, "sprite", "block")

    x = Int(get(entity.data, "x", 0))
    y = Int(get(entity.data, "y", 0))

    width = Int(get(entity.data, "width", 32))
    height = Int(get(entity.data, "height", 32))

    renderGateSwitch(ctx, x, y, width, height, sprite)
end

end

max_stars_repo_name: catapillie/VortexHelper

max_issues_repo_name: catapillie/VortexHelper

max_forks_repo_name: catapillie/VortexHelper


text:
using RCall

#install.packages("rgexf", dependencies=TRUE)
R"library(igraph);library(rgexf)"
R"g = gexf.to.igraph(read.gexf('MCMPOSTER.gexf'))"
R"weights = E(g)$weight"
R"weights = (weights-min(weights))"
R"E(g)$weight= weights/max(weights)"

 R"""undirected = function(x){
     if (length(x)==1) {  return(x[1])}
     else if (length(x)>2) {stop('Too many variables, use igraphs simplify function')}
     else {return(abs(x[1]-x[2]))}
         }"""#makes graph undirected

R"g = simplify(g, edge.attr.comb='sum')"
R"g = as.undirected(g, mode = c('collapse'),  edge.attr.comb = undirected)"

R"write_graph(g, 'mcl.input', 'ncol')"

cmd =  `/Users/$(ENV["USER"])/local/bin/mcl mcl.input --abc -o mcl.out -I 1.4`
run(   cmd  )

f = open("mcl.input")
s = readstring(f)
close(f)

using PyCall,DataFrames,RCall
unshift!(PyVector(pyimport("sys")["path"]), "")
@pyimport ncdata

filename = "./dunmcm.nc"

data = ncdata.get(filename,grp=0)
specs = names!(DataFrame(data["spec"]),[Symbol(i)for i in data["sc"]])
specs = specs[:,[Symbol(i) for i in filter(x->!ismatch(r"EMISS|DUMMY",x),data["sc"])]]

t = 216

R"v = V(g)$name"
@rget v

vv= [specs[Symbol(i)][t] for i in v ]
@rput vv
R"v = log10(vv)"
R"v= (v+abs(min(v)))"
R"v=1e-10 + v/max(v)"
R"V(g)$conc = v"

#=
for i in v
ii = "z"*i
R"g=g+vertex($(ii), conc=0)"
R"g = g + edge($(ii), $(i), weight =1000000)"
end
=#

R"""g1.gexf <- igraph.to.gexf(g)
f <- file('mcmspheres.gexf')
writeLines(g1.gexf$graph, con = f)
close(f)
"""

data = split(s,"\n")
pop!(data)

'''
data = [push!(split(i," "),"Borneo MCM") for i in data]
open("reactionedges.json", "w") do f
        write(f, replace(JSON.json(data),r"\"(\d+\.\d+)\"",s"\1"))
     end
    '''

     #=
n = 100 # length(names(specs))

label = R"a = cluster_optimal(g,weights=weights)"

f = open("spinglass.out","w")
for i in 1:length(label)
R"str = paste(unlist(a[$i]),collapse='\t')"
@rget str
if str != ""
write(f,str)
write(f,"\n")
end
end
close(f)

println("spinglass")
=#
print("fini")

max_stars_repo_name: wolfiex/RopaJL

max_issues_repo_name: wolfiex/RopaJL

max_forks_repo_name: wolfiex/RopaJL


text:
# Functions that find special events, like e.g. lyrics, or tracknames
# are contained here

const NOTRACKNAME = "No track name found"

"""
    trackname(track::MIDI.MIDITrack)

Return the name of the given `track` as a string,
by finding the `TrackNameEvent`.

If no such event exists, `"No track name found"` is returned.
"""
function trackname(track::MIDI.MIDITrack)
    for event in track.events
        if event isa TrackNameEvent
            return event.text
        end
    end
    return NOTRACKNAME
end

"`tracknames(m::MIDIFile) = trackname.(m.tracks)`"
tracknames(m::MIDIFile) = trackname.(m.tracks)

"""
    addtrackname!(track::MIDI.MIDITrack, name::String)

Add a name to the given `track` by attaching the
`TrackNameEvent` to the start of the `track`.
"""
function addtrackname!(track::MIDI.MIDITrack, name::String)
    trackname = TrackNameEvent(0, name)

    # Remove existing name
    for (i, event) in enumerate(track.events)
        if event isa TrackNameEvent
            deleteat!(track.events, i)
            break
        end
    end

    addevent!(track, 0, trackname)
end

"""
    findtextevents(eventtype, track)
Find all text events specifield by `eventtype` in the `track`.
The `eventtype` can be `TextEvent, LyricEvent, MarkerEvent, which will find the
appropriate meta events.

For convenience, this function does not return the events themselves.
Instead, it returns three vectors: the first is the strings of the events,
the second is the indices of the events in the `track` and the third is
the absolute position of the events (since start of `track`).

*Notice* - common music score editors like e.g. MuseScore, GuitarPro, etc., do
not export the lyrics and text information when exporting midi files.

*Notice* - Cubase can read the marker events and MuseScore can read the lyrics
events. We haven't seen any editor that can read the text events, so far.
"""
function findtextevents(eventtype, track)
    events = String[]
    idxs = Int[]
    for (i, event) ∈ enumerate(track.events)
        if event isa eventtype
            # lol without copy the events are lost!
            push!(events, String(copy(event.data)))
            push!(idxs, i)
        end
    end
    length(idxs) == 0 && error("Found no events of such type")
    abspos = get_abs_pos(idxs, track)
    return events, idxs, abspos
end


text:
type TwCoordinate
  lon::Float32
  lat::Float32
end

function TwCoordinate(d::Dict)
  TwCoordinate(
    get(d, "coordinates", [0])[1],
    get(d, "coordinates", [0,0])[2]
  )
end

max_stars_repo_name: nick-paul/JuliaTwitter

max_issues_repo_name: nick-paul/JuliaTwitter

max_forks_repo_name: nick-paul/JuliaTwitter


text:
for fn in (:yearmonthday, :yearmonth, :monthday, :dayofmonth,
           :dayofweek, :isleapyear, :daysinmonth, :daysinyear, :dayofyear, :dayname, :dayabbr,
           :dayofweekofmonth, :daysofweekinmonth, :monthname, :monthabbr,
           :quarterofyear, :dayofquarter)
    @eval $fn(td::TimeDate) = $fn(Date(td))
    @eval $fn(tdz::TimeDateZone) = $fn(Date(tdz))
end

for fn in (:firstdayofweek, :lastdayofweek, :firstdayofmonth, :lastdayofmonth,
           :firstdayofyear, :lastdayofyear, :firstdayofquarter, :lastdayofquarter)
    @eval $fn(td::TimeDate) = $fn(Date(td))
    @eval $fn(tdz::TimeDateZone) = TimeDateZone($fn(Date(tdz)),in_zone(tdz))
end

function tonext(td::TimeDate, dow::Int64)
   fast_time = Microsecond(td) + Nanosecond(td)
   dtm = DateTime(td)
   dtm = tonext(dtm, dow)
   tmdt = TimeDate(dtm)
   if !isempty(fast_time)
       tmdt += fast_time
   end
   return tmdt
end
tonext(td::TimeDate, dow::Int32) = tonext(td, Int64(dow))

function toprev(td::TimeDate, dow::Int64)
   fast_time = Microsecond(td) + Nanosecond(td)
   dtm = DateTime(td)
   dtm = toprev(dtm, dow)
   tmdt = TimeDate(dtm)
   if !isempty(fast_time)
       tmdt += fast_time
   end
   return tmdt
end
toprev(td::TimeDate, dow::Int32) = toprev(td, Int64(dow))

function tonext(fn::Function, td::TimeDate)
   fast_time = Microsecond(td) + Nanosecond(td)
   dtm = DateTime(td)
   dtm = tonext(fn, dtm)
   tmdt = TimeDate(dtm)
   if !isempty(fast_time)
       tmdt += fast_time
   end
   return tmdt
end

function toprev(fn::Function, td::TimeDate)
   fast_time = Microsecond(td) + Nanosecond(td)
   dtm = DateTime(td)
   dtm = toprev(fn, dtm)
   tmdt = TimeDate(dtm)
   if !isempty(fast_time)
       tmdt += fast_time
   end
   return tmdt
end

canonical(x::ZonedDateTime) = x
canonical(x::TimeDateZone) = x


text:
struct AcceleratedArray{T, N, A <: AbstractArray{T, N}, I <: AbstractIndex} <: AbstractArray{T, N}
	parent::A
	index::I
end

AcceleratedVector{T, A <: AbstractVector{T}, I <: AbstractIndex} = AcceleratedArray{T, 1, A, I}
AcceleratedMatrix{T, A <: AbstractMatrix{T}, I <: AbstractIndex} = AcceleratedArray{T, 2, A, I}

Base.parent(a::AcceleratedArray) = a.parent
Base.IndexStyle(::Type{AcceleratedArray{<:Any, <:Any, A, <:Any}}) where A = IndexStyle(A)
Base.axes(a::AcceleratedArray) = axes(parent(a))
Base.size(a::AcceleratedArray) = size(parent(a))
@propagate_inbounds Base.getindex(a::AcceleratedArray, i::Int) = getindex(parent(a), i)
@propagate_inbounds function Base.getindex(a::AcceleratedArray{<:Any, N}, inds::Vararg{Int, N}) where {N}
	getindex(parent(a), inds...)
end

function Base.summary(io::IO, a::AcceleratedArray)
	print(io, summary(parent(a)), " + ", summary(a.index))
end

# Disable mutation for now... simplifies the acceleration indices


text:
using ADCME
using PyCall
using LinearAlgebra
using PyPlot
using Random
using FwiFlow
# include("../ops.jl")
# Random.seed!(233)

function saturation(s0,dporodt,pt,perm,poro,qw,qo,muw,muo,sref,dt,h)
    saturation_ = load_op_and_grad("./build/libSaturation","saturation")
    s0,dporodt,pt,perm,poro,qw,qo,muw,muo,sref,dt,h = convert_to_tensor([s0,dporodt,pt,perm,poro,qw,qo,muw,muo,sref,dt,h], [Float64,Float64,Float64,Float64,Float64,Float64,Float64,Float64,Float64,Float64,Float64,Float64])
    saturation_(s0,dporodt,pt,perm,poro,qw,qo,muw,muo,sref,dt,h)
end

# TODO: specify your input parameters
m = 100
n = 10
h = 0.01
x = zeros(n, m)
y = zeros(n, m)
for i = 1:m
    for j = 1:n
        x[j, i] = h*i
        y[j, i] = h*j
    end
end
t = 3.0
s0 = @. (x^2 + y^2)/(1+x^2+y^2) * exp(-t)
dporodt = exp(t) * zeros(n, m)
pt = @. (x^2+y^2)
perm = rand(n, m)
poro = exp(t) * ones(n, m)
qw = ones(n, m)
qo = ones(n, m)
muw = 2.0
muo = 3.0
sref = s0
dt = 0.01
h = 0.1
u = sat_op2(s0,dporodt,pt,perm,poro,qw,qo,muw,muo,sref,dt,h)
u3 = sat_op(s0,pt,perm,poro,qw,qo,muw, muo,sref,dt,h)

# u3 = sat_op(s0,dporodt,pt,perm,poro,qw,qo,muw,muo,sref,dt,h)
sess = Session(); init(sess)
# @show run(sess, u)

@show run(sess, u3-u)

# uncomment it for testing gradients
error()

# TODO: change your test parameter to `m`
#       in the case of `multiple=true`, you also need to specify which component you are testings
# gradient check -- v
# s0 qw qo
function scalar_function(m)
    # return sum(saturation(s0,dporodt,pt,perm,poro,qw,qo,muw,muo,sref,dt,h)^2)
    return sum(sat_op2(m,dporodt,pt,perm,poro,qw,qo,muw,muo,sref,dt,h)^2)
    # return sum(sat_op(m,pt,perm,poro,qw,qo,muw,muo,sref,dt,h)^2)
end

# TODO: change `m_` and `v_` to appropriate values
m_ = constant(0.7*rand(n, m))
v_ = rand(n, m)

# m_ = constant(s0)
# v_ = rand(n, m)

y_ = scalar_function(m_)
dy_ = gradients(y_, m_)
ms_ = Array{Any}(undef, 5)
ys_ = Array{Any}(undef, 5)
s_ = Array{Any}(undef, 5)
w_ = Array{Any}(undef, 5)
gs_ =  @. 1 / 10^(1:5)

for i = 1:5
    g_ = gs_[i]
    ms_[i] = m_ + g_*v_
    ys_[i] = scalar_function(ms_[i])
    s_[i] = ys_[i] - y_
    w_[i] = s_[i] - g_*sum(v_.*dy_)
end

sess = Session(); init(sess)
sval_ = run(sess, s_)
wval_ = run(sess, w_)
close("all")
loglog(gs_, abs.(sval_), "*-", label="finite difference")
loglog(gs_, abs.(wval_), "+-", label="automatic differentiation")
loglog(gs_, gs_.^2 * 0.5*abs(wval_[1])/gs_[1]^2, "--",label="\$\\mathcal{O}(\\gamma^2)\$")
loglog(gs_, gs_ * 0.5*abs(sval_[1])/gs_[1], "--",label="\$\\mathcal{O}(\\gamma)\$")

plt.gca().invert_xaxis()
legend()
xlabel("\$\\gamma\$")
ylabel("Error")


text:
function main(::Robot)
  a = []
  summa = 0
  side = Ost
  while !isborder(r, Nord)
  movess!(r, side, a)
    if side == Ost
      side = West
    else
      side = Ost
    end
  end
  movess!(r, side, a)
  k = length(a)
  while length(a) > 0
    summa += pop!(a)
  end
  println(summa/k)

end

function movess!(::Robot, side, a)
  while !isborder(r, side)
    if ismarker(r)
      push!(a, temperature(r))
    end
    move!(r,side)
  end
  if !isborder(r, Nord)
    move!(r, Nord)
  end
end

max_stars_repo_name: VladislavBukin/mirea-progs

max_issues_repo_name: VladislavBukin/mirea-progs

max_forks_repo_name: VladislavBukin/mirea-progs


text:
#=
NOTE currently we do not restrict the sparsity pattern to be chordal here (at the
cost of not being able to obtain "closed form" hess sqrt and inv hess oracles)

TODO
- improve efficiency of hessian calculations using structure
- use inbounds
- improve efficiency of frontal matrix scattering etc
- maybe allow passing more options to CHOLMOD
=#

import SuiteSparse.CHOLMOD

"""
$(TYPEDEF)

CHOLMOD sparse Cholesky-based implementation for the sparse positive semidefinite
cone [`PosSemidefTriSparse`](@ref). Note only BLAS floating point types are
supported.
"""
struct PSDSparseCholmod <: PSDSparseImpl end

mutable struct PSDSparseCholmodCache{T <: BlasReal, R <: RealOrComplex{T}} <:
    PSDSparseCache{T, R}
    sparse_point
    sparse_point_map
    symb_mat
    supers
    super_map
    parents
    ancestors
    num_cols
    num_rows
    J_rows
    L_idxs
    F_blocks
    L_blocks
    S_blocks
    L_pr_blocks
    S_pr_blocks
    L_pr_pr_blocks
    map_blocks
    temp_blocks
    rel_idxs
    PSDSparseCholmodCache{T, R}() where {T <: BlasReal, R <: RealOrComplex{T}} =
        new{T, R}()
end

function setup_extra_data!(
    cone::PosSemidefTriSparse{PSDSparseCholmod, T, R},
    ) where {T <: BlasReal, R <: RealOrComplex{T}}
    cone.cache = cache = PSDSparseCholmodCache{T, R}()

    # setup symbolic factorization
    side = cone.side
    dim_R = length(cone.row_idxs)

    fake_point = [R(l) for l in 1:dim_R]
    sparse_point = cache.sparse_point = CHOLMOD.Sparse(Hermitian(
        sparse(cone.row_idxs, cone.col_idxs, fake_point, side, side), :L))
    sparse_point_map = cache.sparse_point_map = zeros(Int, dim_R)
    sx_ptr = Ptr{R}(unsafe_load(sparse_point.ptr).x)
    for l in 1:dim_R
        sparse_point_map[Int(unsafe_load(sx_ptr, l))] = l
    end
    CHOLMOD.@cholmod_param supernodal = 2 begin
        cache.symb_mat = CHOLMOD.symbolic(sparse_point, postorder=true)
    end
    symb_mat = cache.symb_mat

    f = unsafe_load(pointer(symb_mat))
    @assert f.n == cone.side
    @assert f.is_super != 0

    num_super = Int(f.nsuper)
    supers = cache.supers = [unsafe_load(f.super, i) + 1 for i in 1:num_super]
    push!(supers, side + 1)
    fs = [unsafe_load(f.s, i) + 1 for i in 1:Int(f.ssize)]
    fpi = [unsafe_load(f.pi, i) + 1 for i in 1:num_super]
    push!(fpi, length(fs) + 1)
    fpx = [unsafe_load(f.px, i) + 1 for i in 1:num_super]

    # construct super_map
    # super_map[k] = s if column k is in supernode s
    super_map = cache.super_map = zeros(Int, side)
    for s in 1:num_super, k in supers[s]:(supers[s + 1] - 1)
        super_map[k] = s
    end

    # construct supernode tree
    parents = cache.parents = zeros(Int, num_super)
    for k in 1:num_super
        nn = supers[k + 1] - supers[k] # n cols
        nj = fpi[k + 1] - fpi[k] # n rows
        na = nj - nn # n rows below tri
        @assert nn >= 0
        @assert nj >= 0
        @assert na >= 0

        if !iszero(na)
            kparloc = fs[fpi[k] + nn]
            kparent = super_map[kparloc]
            @assert supers[kparent] <= kparloc < supers[kparent + 1]
            parents[k] = kparent
        end
    end
    ancestors = cache.ancestors = [Int[] for k in 1:num_super]
    for k in 1:num_super
        curr = k
        while true
            push!(ancestors[k], curr)
            curr = parents[curr]
            if curr == 0
                break
            end
        end
    end

    num_cols = cache.num_cols = Vector{Int}(undef, num_super)
    num_rows = cache.num_rows = Vector{Int}(undef, num_super)
    L_idxs = cache.L_idxs = Vector{UnitRange{Int}}(undef, num_super)
    J_rows = cache.J_rows = Vector{Vector}(undef, num_super)
    F_blocks = cache.F_blocks = Vector{Matrix{R}}(undef, num_super)
    L_blocks = cache.L_blocks = Vector{Matrix{R}}(undef, num_super)
    L_pr_blocks = cache.L_pr_blocks = Vector{Matrix{R}}(undef, num_super)
    L_pr_pr_blocks = cache.L_pr_pr_blocks = Vector{Matrix{R}}(undef, num_super)
    S_blocks = cache.S_blocks = Vector{Matrix{R}}(undef, num_super)
    S_pr_blocks = cache.S_pr_blocks = Vector{Matrix{R}}(undef, num_super)
    temp_blocks = cache.temp_blocks = Vector{Matrix{R}}(undef, num_super)
    rel_idxs = cache.rel_idxs = [Tuple{Int, Int}[] for k in 1:num_super]
    for k in reverse(1:num_super)
        num_col = num_cols[k] = supers[k + 1] - supers[k]
        num_row = num_rows[k] = fpi[k + 1] - fpi[k] # n rows
        @assert 0 < num_col <= num_row

        J_k = J_rows[k] = fs[fpi[k]:(fpi[k + 1] - 1)]
        @assert length(J_k) == num_row

        L_idxs[k] = fpx[k] - 1 .+ (1:(num_row * num_col))

        F_blocks[k] = zeros(R, num_row, num_row)
        num_below = num_row - num_col
        L_blocks[k] = zeros(R, num_row, num_col)
        L_pr_blocks[k] = zeros(R, num_row, num_col)
        L_pr_pr_blocks[k] = zeros(R, num_row, num_col)
        S_blocks[k] = zeros(R, num_below, num_below)
        S_pr_blocks[k] = zeros(R, num_below, num_below)
        temp_blocks[k] = zeros(R, num_row, num_col)

        if num_row > num_col
            # setup relative indices in frontal matrix
            rel_idx = rel_idxs[k]
            k_par = parents[k]
            I_k = J_k[(num_col + 1):end]
            for (idx_i, i) in enumerate(I_k), (idx_j, j) in
                enumerate(J_rows[k_par])
                if i == j
                    push!(rel_idx, (idx_i, idx_j))
                end
            end
        end
    end

    iperm = invperm(symb_mat.p)
    map_blocks = cache.map_blocks = Vector{Tuple{Int, Int, Int, Bool, Bool}}(
        undef, dim_R)
    for i in 1:dim_R
        row = iperm[cone.row_idxs[i]]
        col = iperm[cone.col_idxs[i]]
        if row < col
            (row, col) = (col, row)
            swapped = true
        else
            swapped = false
        end
        super = super_map[col]
        J = J_rows[super]
        row_idx = findfirst(r -> r == row, J)
        col_idx = col - supers[super] + 1
        map_blocks[i] = (super, row_idx, col_idx, row != col, swapped)
    end

    return
end

function update_feas(
    cone::PosSemidefTriSparse{PSDSparseCholmod, T, R},
    ) where {T <: BlasReal, R <: RealOrComplex{T}}
    @assert !cone.feas_updated
    point = cone.point
    cache = cone.cache
    sparse_point = cache.sparse_point

    # svec scale point and copy directly to CHOLDMOD Sparse data structure
    sx_ptr = Ptr{R}(unsafe_load(sparse_point.ptr).x)
    if cone.is_complex
        idx = 1
        @inbounds for (p_idx, s_idx) in enumerate(cache.sparse_point_map)
            if cone.row_idxs[p_idx] == cone.col_idxs[p_idx]
                p = Complex(point[idx])
                idx += 1
            else
                p = Complex(point[idx], point[idx + 1]) / cone.rt2
                idx += 2
            end
            unsafe_store!(sx_ptr, p, s_idx)
        end
    else
        @inbounds for (p_idx, s_idx) in enumerate(cache.sparse_point_map)
            p = point[p_idx]
            if cone.row_idxs[p_idx] != cone.col_idxs[p_idx]
                p /= cone.rt2
            end
            unsafe_store!(sx_ptr, p, s_idx)
        end
    end

    # update numeric factorization
    CHOLMOD.@cholmod_param supernodal = 2 begin
        cone.is_feas = isposdef(CHOLMOD.cholesky!(cache.symb_mat, sparse_point;
            check = false))
    end

    cone.feas_updated = true
    return cone.is_feas
end

function update_grad(
    cone::PosSemidefTriSparse{PSDSparseCholmod, T, R},
    ) where {T <: BlasReal, R <: RealOrComplex{T}}
    @assert cone.is_feas
    cache = cone.cache
    num_cols = cache.num_cols
    num_rows = cache.num_rows

    # update L blocks from CHOLMOD numerical factorization
    lx_ptr = Ptr{R}(unsafe_load(pointer(cache.symb_mat)).x)
    @inbounds for k in 1:length(num_cols)
        L_block = cache.L_blocks[k]
        for (l, lx_idx) in enumerate(cache.L_idxs[k])
            L_block[l] = unsafe_load(lx_ptr, lx_idx)
        end
    end

    # build inv blocks
    @inbounds for k in reverse(1:length(num_cols))
        num_col = num_cols[k]
        num_row = num_rows[k]
        idxs_n = 1:num_col
        F_block = cache.F_blocks[k]
        L_block = cache.L_blocks[k]

        @views L_n = L_block[idxs_n, :]
        @views F_nn = F_block[idxs_n, idxs_n]
        copyto!(F_nn, L_n)
        LAPACK.potri!('L', F_nn)

        if num_row > num_col
            idxs_a = (num_col + 1):num_row
            @views F_aa = Hermitian(F_block[idxs_a, idxs_a], :L)
            @views F_an = F_block[idxs_a, idxs_n]
            @views L_a = rdiv!(L_block[idxs_a, :], LowerTriangular(L_n))

            F_aa.data .= 0
            F_par = Hermitian(cache.F_blocks[cache.parents[k]], :L)
            rel_idx = cache.rel_idxs[k]
            for (i, j) in rel_idx, (i2, j2) in rel_idx # TODO only lower tri
                F_aa.data[i, i2] = F_par[j, j2]
            end

            mul!(F_an, F_aa, L_a, -1, false)
            mul!(F_nn, F_an', L_a, -1, true)

            copyto!(cache.S_blocks[k], F_aa.data) # use in Hessian calculations
        end

        @views cache.temp_blocks[k] = F_block[:, idxs_n]
    end

    smat_to_svec_sparse!(cone.grad, cache.temp_blocks, cone)
    cone.grad .*= -1

    cone.grad_updated = true
    return cone.grad
end

# for each column in identity, get hess prod to build explicit hess
function update_hess(cone::PosSemidefTriSparse{PSDSparseCholmod})
    @assert cone.grad_updated
    isdefined(cone, :hess) || alloc_hess!(cone)
    cache = cone.cache
    temp_blocks = cache.temp_blocks
    rt2 = cone.rt2
    invrt2 = inv(rt2)
    H = cone.hess.data

    H_idx_j = 1
    @inbounds for (j, (super_j, row_idx_j, col_idx_j, scal_j, swapped_j)) in
        enumerate(cache.map_blocks)
        for H_block in temp_blocks
            H_block .= 0
        end

        if cone.is_complex
            # real part j
            temp_blocks[super_j][row_idx_j, col_idx_j] = (scal_j ? invrt2 : 1)
            _hess_step1(cone, cache.ancestors[super_j])
            _hess_step2(cone, cache.ancestors[super_j], false)
            out_blocks = _hess_step3(cone)
            H_idx_i = 1
            for i in 1:j
                (super_i, row_idx_i, col_idx_i, scal_i, swapped_i) =
                    cache.map_blocks[i]
                HijR = out_blocks[super_i][row_idx_i, col_idx_i]
                if scal_i
                    HijR *= rt2
                end
                H[H_idx_i, H_idx_j] = real(HijR) # real part i
                H_idx_i += 1
                if row_idx_i != col_idx_i
                    # complex part i
                    H[H_idx_i, H_idx_j] = swapped_i ? -imag(HijR) : imag(HijR)
                    H_idx_i += 1
                end
            end
            H_idx_j += 1

            if row_idx_j != col_idx_j
                # complex part j
                for H_block in temp_blocks
                    H_block .= 0
                end
                temp_blocks[super_j][row_idx_j, col_idx_j] =
                    (scal_j ? invrt2 : 1) * im
                _hess_step1(cone, cache.ancestors[super_j])
                _hess_step2(cone, cache.ancestors[super_j], false)
                out_blocks = _hess_step3(cone)
                H_idx_i = 1
                for i in 1:j
                    (super_i, row_idx_i, col_idx_i, scal_i, swapped_i) =
                        cache.map_blocks[i]
                    HijR = out_blocks[super_i][row_idx_i, col_idx_i]
                    if scal_i
                        HijR *= rt2
                    end
                    # real part i
                    H[H_idx_i, H_idx_j] = swapped_j ? -real(HijR) : real(HijR)
                    H_idx_i += 1
                    if row_idx_i != col_idx_i
                        # complex part i
                        H[H_idx_i, H_idx_j] = xor(swapped_i, swapped_j) ?
                            -imag(HijR) : imag(HijR)
                        H_idx_i += 1
                    end
                end
                H_idx_j += 1
            end
        else
            temp_blocks[super_j][row_idx_j, col_idx_j] = (scal_j ? invrt2 : 1)
            _hess_step1(cone, cache.ancestors[super_j])
            _hess_step2(cone, cache.ancestors[super_j], false)
            out_blocks = _hess_step3(cone)

            for i in 1:j
                (super_i, row_idx_i, col_idx_i, scal_i, swapped_i) =
                    cache.map_blocks[i]
                H[i, j] = out_blocks[super_i][row_idx_i, col_idx_i]
                if scal_i
                    H[i, j] *= rt2
                end
            end
        end
    end

    cone.hess_updated = true
    return cone.hess
end

function hess_prod_slow!(
    prod::AbstractVecOrMat,
    arr::AbstractVecOrMat,
    cone::PosSemidefTriSparse,
    )
    cone.use_hess_prod_slow_updated || update_use_hess_prod_slow(cone)
    @assert cone.hess_updated
    cone.use_hess_prod_slow || return hess_prod!(prod, arr, cone)
    @assert cone.grad_updated
    cache = cone.cache
    temp_blocks = cache.temp_blocks

    @inbounds for i in 1:size(arr, 2)
        @views svec_to_smat_sparse!(temp_blocks, arr[:, i], cone)
        _hess_step1(cone, eachindex(cache.num_cols))
        _hess_step2(cone, eachindex(cache.num_cols), false)
        _hess_step3(cone)
        @views smat_to_svec_sparse!(prod[:, i], temp_blocks, cone)
    end

    return prod
end

function _hess_step1(
    cone::PosSemidefTriSparse{PSDSparseCholmod},
    supernode_list::AbstractVector{Int},
    )
    @assert cone.grad_updated
    cache = cone.cache
    temp_blocks = cache.temp_blocks

    @inbounds for k in supernode_list
        num_col = cache.num_cols[k]
        num_row = cache.num_rows[k]
        idxs_n = 1:num_col
        idxs_a = (num_col + 1):num_row
        F_block = cache.F_blocks[k]
        temp_block = temp_blocks[k]

        @views F_block[idxs_a, idxs_a] .= 0
        @views F_block[:, idxs_n] .= temp_block
    end

    @inbounds for k in supernode_list
        num_col = cache.num_cols[k]
        num_row = cache.num_rows[k]
        idxs_n = 1:num_col
        F_block = cache.F_blocks[k]
        temp_block = temp_blocks[k]

        if num_row > num_col
            outer_L_prod(cone, k)
            idxs_a = (num_col + 1):num_row
            @views F_aa = Hermitian(F_block[idxs_a, idxs_a], :L)

            F_par = cache.F_blocks[cache.parents[k]]
            rel_idx = cache.rel_idxs[k]
            for (i, j) in rel_idx, (i2, j2) in rel_idx # TODO only lower tri
                F_par[j, j2] += F_aa[i, i2]
            end
        end

        @views copyto!(temp_block, F_block[:, idxs_n])
    end

    return temp_blocks
end

function _hess_step2(
    cone::PosSemidefTriSparse{PSDSparseCholmod},
    supernode_list::AbstractVector{Int},
    save_L_pr::Bool,
    )
    cache = cone.cache
    temp_blocks = cache.temp_blocks

    @inbounds for k in supernode_list
        num_col = cache.num_cols[k]
        num_row = cache.num_rows[k]
        idxs_n = 1:num_col
        temp_block = temp_blocks[k]
        L_pr_block = cache.L_pr_blocks[k]

        @views temp_block_n = temp_block[idxs_n, :]
        copytri!(temp_block_n, 'L', true)
        @views L_n = LowerTriangular(cache.L_blocks[k][idxs_n, :])
        if save_L_pr
            @views copyto!(L_pr_block[idxs_n, :], temp_block_n)
        end
        ldiv!(L_n, temp_block_n)
        ldiv!(L_n', temp_block_n)
        rdiv!(temp_block, L_n')
        rdiv!(temp_block, L_n)

        if num_row > num_col
            idxs_a = (num_col + 1):num_row
            @views temp_block_a = temp_block[idxs_a, :]
            if save_L_pr
                @views copyto!(L_pr_block[idxs_a, :], temp_block_a)
            end
            @views F_an = cache.F_blocks[k][idxs_a, idxs_n]
            mul!(F_an, Hermitian(cache.S_blocks[k], :L), temp_block_a)
            copyto!(temp_block_a, F_an)
        end
    end

    return temp_blocks
end

function _hess_step3(
    cone::PosSemidefTriSparse{PSDSparseCholmod},
    )
    cache = cone.cache
    temp_blocks = cache.temp_blocks

    @inbounds for k in reverse(1:length(cache.num_cols))
        num_col = cache.num_cols[k]
        num_row = cache.num_rows[k]
        idxs_n = 1:num_col
        F_block = cache.F_blocks[k]
        temp_block = temp_blocks[k]
        @. @views F_block[:, idxs_n] = temp_block
        if num_row > num_col
            idxs_a = (num_col + 1):num_row
            @views L_a = cache.L_blocks[k][idxs_a, :]
            @views F_aa = Hermitian(F_block[idxs_a, idxs_a], :L)
            @views F_an = F_block[idxs_a, idxs_n]
            @views F_nn = Hermitian(F_block[idxs_n, idxs_n], :L)
            F_aa.data .= 0
            F_par = Hermitian(cache.F_blocks[cache.parents[k]], :L)
            rel_idx = cache.rel_idxs[k]
            for (i, j) in rel_idx, (i2, j2) in rel_idx # TODO only lower tri
                F_aa.data[i, i2] = F_par[j, j2]
            end
            mul!(F_nn.data, F_an', L_a, -1, true)
            mul!(F_an, F_aa, L_a, -1, true)
            mul!(F_nn.data, L_a', F_an, -1, true)
        end
        @views copyto!(temp_block, F_block[:, idxs_n])
    end

    return temp_blocks
end

function outer_L_prod(
    cone::PosSemidefTriSparse{PSDSparseCholmod},
    k::Int,
    )
    cache = cone.cache

    num_col = cache.num_cols[k]
    num_row = cache.num_rows[k]
    idxs_n = 1:num_col
    F_block = cache.F_blocks[k]
    idxs_a = (num_col + 1):num_row
    @views L_a = cache.L_blocks[k][idxs_a, :]
    @views F_aa = Hermitian(F_block[idxs_a, idxs_a], :L)
    @views F_an = F_block[idxs_a, idxs_n]
    @views F_nn = Hermitian(F_block[idxs_n, idxs_n], :L)
    mul!(F_aa.data, L_a, F_an', -1, true)
    mul!(F_an, L_a, F_nn, -1, true)
    mul!(F_aa.data, F_an, L_a', -1, true)

    return F_block
end

function dder3(
    cone::PosSemidefTriSparse{PSDSparseCholmod},
    dir::AbstractVector,
    )
    @assert cone.grad_updated
    cache = cone.cache
    temp_blocks = cache.temp_blocks

    @views svec_to_smat_sparse!(temp_blocks, dir, cone)
    _hess_step1(cone, eachindex(cache.num_cols))
    _hess_step2(cone, eachindex(cache.num_cols), true)
    _hess_step3(cone)
    @inbounds for k in eachindex(cache.num_cols)
        idxs_a = (cache.num_cols[k] + 1):cache.num_rows[k]
        F_block = cache.F_blocks[k]
        @. @views cache.S_pr_blocks[k] = F_block[idxs_a, idxs_a]
        F_block .= 0
    end

    @inbounds for k in eachindex(cache.num_cols)
        num_col = cache.num_cols[k]
        num_row = cache.num_rows[k]
        idxs_n = 1:num_col
        idxs_a = (num_col + 1):num_row
        F_block = cache.F_blocks[k]
        @views L_n = LowerTriangular(cache.L_blocks[k][idxs_n, :])
        @views L_pr = cache.L_pr_blocks[k][idxs_a, :]
        @views temp_block_a = temp_blocks[k][idxs_a, :]

        if num_row > num_col
            outer_L_prod(cone, k)
            @views F_aa = F_block[idxs_a, idxs_a]
            mul!(temp_block_a, L_pr, L_n)
            mul!(F_aa, temp_block_a, temp_block_a', -2, true)
            F_par = cache.F_blocks[cache.parents[k]]
            rel_idx = cache.rel_idxs[k]
            for (i, j) in rel_idx, (i2, j2) in rel_idx
                F_par[j, j2] += F_aa[i, i2]
            end
        end

        # transform block from linearized factorization into
        # block for linearized inverse
        S = cache.S_blocks[k]
        L_pr_pr = cache.L_pr_pr_blocks[k]
        @views D_pr = cache.L_pr_blocks[k][idxs_n, :]
        @views temp_block_n = temp_blocks[k][idxs_n, :]
        @views copyto!(L_pr_pr, F_block[:, idxs_n])
        @views L_pr_pr_a = L_pr_pr[idxs_a, :]
        mul!(L_pr_pr_a, L_pr, D_pr, -2, true)
        rdiv!(L_pr_pr_a, L_n')
        rdiv!(L_pr_pr_a, L_n)
        @. @views temp_block_n = -L_pr_pr[idxs_n, :]
        ldiv!(L_n, D_pr)
        mul!(temp_block_n, D_pr', D_pr, 2, true)
        ldiv!(L_n, temp_block_n)
        ldiv!(L_n', temp_block_n)
        rdiv!(temp_block_n, L_n')
        rdiv!(temp_block_n, L_n)
        mul!(temp_block_a, S, L_pr)
        mul!(temp_block_n, L_pr', temp_block_a, 2, true)
        mul!(temp_block_a, cache.S_pr_blocks[k], L_pr, 2, false)
        mul!(temp_block_a, S, L_pr_pr_a, -1, true)
    end

    _hess_step3(cone)
    smat_to_svec_sparse!(cone.dder3, cache.temp_blocks, cone)
    cone.dder3 ./= 2

    return cone.dder3
end

function svec_to_smat_sparse!(
    blocks::Vector{Matrix{T}},
    vec::AbstractVector{T},
    cone::PosSemidefTriSparse{PSDSparseCholmod, T, T},
    ) where {T <: BlasReal}
    @assert length(vec) == length(cone.row_idxs)
    for b in blocks
        b .= 0
    end
    @inbounds for (i, (super, row_idx, col_idx, scal, swapped)) in enumerate(cone.cache.map_blocks)
        vec_i = vec[i]
        if scal
            vec_i /= cone.rt2
        end
        blocks[super][row_idx, col_idx] = vec_i
    end
    return blocks
end

function smat_to_svec_sparse!(
    vec::AbstractVector{T},
    blocks::Vector{Matrix{T}},
    cone::PosSemidefTriSparse{PSDSparseCholmod, T, T},
    ) where {T <: BlasReal}
    @assert length(vec) == length(cone.row_idxs)
    @inbounds for (i, (super, row_idx, col_idx, scal, swapped)) in enumerate(cone.cache.map_blocks)
        vec_i = blocks[super][row_idx, col_idx]
        if scal
            vec_i *= cone.rt2
        end
        vec[i] = vec_i
    end
    return vec
end

function svec_to_smat_sparse!(
    blocks::Vector{Matrix{Complex{T}}},
    vec::AbstractVector{T},
    cone::PosSemidefTriSparse{PSDSparseCholmod, T, Complex{T}},
    ) where {T <: BlasReal}
    for b in blocks
        b .= 0
    end
    idx = 1
    @inbounds for (i, (super, row_idx, col_idx, scal, swapped)) in enumerate(cone.cache.map_blocks)
        if scal
            blocks[super][row_idx, col_idx] = Complex(vec[idx], (swapped ? -vec[idx + 1] : vec[idx + 1])) / cone.rt2
            idx += 2
        else
            blocks[super][row_idx, col_idx] = vec[idx]
            idx += 1
        end
    end
    @assert idx == length(vec) + 1
    return blocks
end

function smat_to_svec_sparse!(
    vec::AbstractVector{T},
    blocks::Vector{Matrix{Complex{T}}},
    cone::PosSemidefTriSparse{PSDSparseCholmod, T, Complex{T}},
    ) where {T <: BlasReal}
    idx = 1
    @inbounds for (i, (super, row_idx, col_idx, scal, swapped)) in enumerate(cone.cache.map_blocks)
        vec_R = blocks[super][row_idx, col_idx]
        if scal
            vec_R *= cone.rt2
            vec[idx] = real(vec_R)
            vec[idx + 1] = swapped ? -imag(vec_R) : imag(vec_R)
            idx += 2
        else
            vec[idx] = real(vec_R)
            idx += 1
        end
    end
    @assert idx == length(vec) + 1
    return vec
end


text:
export ProductMeasure

using MappedArrays
using Base: @propagate_inbounds

struct ProductMeasure{T} <: AbstractMeasure
    data :: T

    ProductMeasure(μs...) = new{typeof(μs)}(μs)
    ProductMeasure(μs) = new{typeof(μs)}(μs)
end

Base.size(μ::ProductMeasure) = size(μ.data)

function Base.show(io::IO, μ::ProductMeasure{NamedTuple{N,T}}) where {N,T}
    io = IOContext(io, :compact => true)
    print(io, "Product(",μ.data, ")")
end

function Base.show(io::IO, μ::ProductMeasure{T}) where {T <: Tuple}
    io = IOContext(io, :compact => true)
    print(io, join(string.(μ.data), " ⊗ "))
end

function Base.show_unquoted(io::IO, μ::ProductMeasure, indent::Int, prec::Int)
    if Base.operator_precedence(:*) ≤ prec
        print(io, "(")
        show(io, μ)
        print(io, ")")
    else
        show(io, μ)
    end
    return nothing
end

# ProductMeasure(m::NTuple{N, Measure{X}}) where {N,X} = ProductMeasure(m...)

Base.length(m::ProductMeasure{T}) where {T} = length(m.data)

function Base.:*(μ::ProductMeasure{Tuple{}}, ν::N) where {X, N <: AbstractMeasure}
    ProductMeasure((ν,))
end

function Base.:*(μ::ProductMeasure{X}, ν::ProductMeasure{Y}) where {X,Y}
    data = (μ.data..., ν.data...)
    ProductMeasure(data...)
end

function Base.:*(μ, ν::ProductMeasure{Y}) where {Y}
    data = (μ, ν.data...)
    ProductMeasure(data...)
end

function Base.:*(μ::ProductMeasure{X}, ν::N) where {X, N <: AbstractMeasure}
    data = (μ.data..., ν)
    ProductMeasure(data...)
end

function Base.:*(μ::M, ν::N) where {M <: AbstractMeasure, N <: AbstractMeasure}
    data = (μ, ν)
    ProductMeasure(data...)
end

using Tullio

function logdensity(d::ProductMeasure, x::NTuple)
    mapreduce(logdensity, +, d.data, x)
end

function logdensity(d::ProductMeasure{A}, x::NTuple) where {A <: AbstractArray}
    mapreduce(logdensity, +, d.data, x)
end

function logdensity(d::ProductMeasure{A}, x) where {T, A<:AbstractArray}
    data = d.data
    @boundscheck size(data) == size(x) || throw(BoundsError)
    s = 0.0
    Δs(i) = logdensity(data[i], x[i])
    for i in eachindex(x)
        s += Δs(i)
    end
    return s
end

# @propagate_inbounds function MeasureTheory.logdensity(d::ProductMeasure{A}, x) where{T, A<:AbstractArray{T,1}}
#     data = d.data
#     @boundscheck size(data) == size(x) || throw(BoundsError)
#     @tullio s = logdensity(data[i], x[i])
#     s
# end

# @propagate_inbounds function MeasureTheory.logdensity(d::ProductMeasure{A}, x) where{T, A<:AbstractArray{T,2}}
#     data = d.data
#     @boundscheck size(data) == size(x) || throw(BoundsError)
#     @tullio s = @inbounds logdensity(data[i,j], x[i,j])
#     s
# end

# @propagate_inbounds function MeasureTheory.logdensity(d::ProductMeasure{A}, x) where{T, A<:AbstractArray{T,3}}
#     data = d.data
#     @boundscheck size(data) == size(x) || throw(BoundsError)
#     @tullio s = @inbounds logdensity(data[i,j,k], x[i,j,k])
#     s
# end

# @propagate_inbounds function MeasureTheory.logdensity(d::ProductMeasure{A}, x) where{T, A<:AbstractArray{T,4}}
#     data = d.data
#     @boundscheck size(data) == size(x) || throw(BoundsError)
#     @tullio s = @inbounds logdensity(data[i,j,k,l], x[i,j,k,l])
#     s
# end

# @propagate_inbounds function MeasureTheory.logdensity(d::ProductMeasure{A}, x) where{T, A<:AbstractArray{T,5}}
#     data = d.data
#     @boundscheck size(data) == size(x) || throw(BoundsError)
#     @tullio s = @inbounds logdensity(data[i,j,k,l,m], x[i,j,k,l,m])
#     s
# end

export rand!
using Random: rand!, GLOBAL_RNG, AbstractRNG

@propagate_inbounds function Random.rand!(rng::AbstractRNG, d::ProductMeasure, x::AbstractArray)
    @boundscheck size(d.data) == size(x) || throw(BoundsError)

    @inbounds for j in eachindex(x)
        x[j] = rand(rng, eltype(x), d.data[j])
    end
    x
end

# function Base.rand(rng::AbstractRNG, d::ProductMeasure)
#     return rand(rng, sampletype(d), d)
# end

# function Base.rand(T::Type, d::ProductMeasure)
#     return rand(Random.GLOBAL_RNG, T, d)
# end

# function Base.rand(d::ProductMeasure)
#     T = sampletype(d)
#     return rand(Random.GLOBAL_RNG, T, d)
# end

function sampletype(d::ProductMeasure{A}) where {T,N,A <: AbstractArray{T,N}}
    S = @inbounds sampletype(d.data[1])
    Array{S, N}
end

function sampletype(d::ProductMeasure{<: Tuple})
    Tuple{sampletype.(d.data)...}
end

# TODO: Pull weights outside
basemeasure(μ::ProductMeasure) = ProductMeasure(basemeasure.(μ.data))

# function logdensity(μ::ProductMeasure{Aμ}, x::Ax) where {Aμ <: MappedArray, Ax <: AbstractArray}
#     μ.data
# end

export marginals

marginals(d::ProductMeasure) = d.data

function testvalue(d::ProductMeasure{NamedTuple{N,T}}) where {N,T}
    vals = values(d.data)
    NamedTuple{N}(testvalue.(vals))
end


text:
### Views
to_index(::BlockIndexRange) = throw(ArgumentError("BlockIndexRange must be converted by to_indices(...)"))

@inline to_indices(A, inds, I::Tuple{BlockIndexRange{1,R}, Vararg{Any}}) where R =
    (unblock(A, inds, I), to_indices(A, _maybetail(inds), tail(I))...)

# splat out higher dimensional blocks
# this mimics view of a CartesianIndex
@inline to_indices(A, inds, I::Tuple{BlockIndexRange, Vararg{Any}}) =
    to_indices(A, inds, (BlockRange.(Block.(I[1].block.n), tuple.(I[1].indices))..., tail(I)...))

# In 0.7, we need to override to_indices to avoid calling linearindices
@inline to_indices(A, I::Tuple{BlockIndexRange, Vararg{Any}}) =
    to_indices(A, axes(A), I)

if VERSION >= v"1.2-"  # See also `reindex` definitions in views.jl
    reindex(idxs::Tuple{BlockSlice{<:BlockRange}, Vararg{Any}},
            subidxs::Tuple{BlockSlice{<:BlockIndexRange}, Vararg{Any}}) =
        (@_propagate_inbounds_meta; (BlockSlice(BlockIndexRange(Block(idxs[1].block.indices[1][Int(subidxs[1].block.block)]),
                                                                subidxs[1].block.indices),
                                                idxs[1].indices[subidxs[1].indices]),
                                    reindex(tail(idxs), tail(subidxs))...))
else  # if VERSION >= v"1.2-"
    reindex(V, idxs::Tuple{BlockSlice{<:BlockRange}, Vararg{Any}},
            subidxs::Tuple{BlockSlice{<:BlockIndexRange}, Vararg{Any}}) =
        (@_propagate_inbounds_meta; (BlockSlice(BlockIndexRange(Block(idxs[1].block.indices[1][Int(subidxs[1].block.block)]),
                                                                subidxs[1].block.indices),
                                                idxs[1].indices[subidxs[1].indices]),
                                        reindex(V, tail(idxs), tail(subidxs))...))
end  # if VERSION >= v"1.2-"

# _splatmap taken from Base:
_splatmap(f, ::Tuple{}) = ()
_splatmap(f, t::Tuple) = (f(t[1])..., _splatmap(f, tail(t))...)

# De-reference blocks before creating a view to avoid taking `global2blockindex`
# path in `AbstractBlockStyle` broadcasting.
@inline function Base.unsafe_view(
        A::BlockArray{<:Any, N},
        I::Vararg{BlockSlice{<:BlockIndexRange{1}}, N}) where {N}
    @_propagate_inbounds_meta
    B = A[map(x -> x.block.block, I)...]
    return view(B, _splatmap(x -> x.block.indices, I)...)
end

@inline function Base.unsafe_view(
        A::PseudoBlockArray{<:Any, N},
        I::Vararg{BlockSlice{<:BlockIndexRange{1}}, N}) where {N}
    @_propagate_inbounds_meta
    return view(A.blocks, map(x -> x.indices, I)...)
end

@inline function Base.unsafe_view(
        A::ReshapedArray{<:Any, N, <:AbstractBlockArray{<:Any, M}},
        I::Vararg{BlockSlice{<:BlockIndexRange{1}}, N}) where {N, M}
    @_propagate_inbounds_meta
    # Note: assuming that I[M+1:end] are verified to be singletons
    return reshape(view(A.parent, I[1:M]...), Val(N))
end

@inline function Base.unsafe_view(
        A::Array{<:Any, N},
        I::Vararg{BlockSlice{<:BlockIndexRange{1}}, N}) where {N}
    @_propagate_inbounds_meta
    return view(A, map(x -> x.indices, I)...)
end

"""
    unblock(block_sizes, inds, I)

Returns the indices associated with a block as a `BlockSlice`.
"""
function unblock(A, inds, I)
    B = first(I)
    if length(inds) == 0
        # Allow `ones(2)[Block(1)[1:1], Block(1)[1:1]]` which is
        # similar to `ones(2)[1:1, 1:1]`.
        BlockSlice(B,Base.OneTo(1))
    else
        BlockSlice(B,inds[1][B])
    end
end

to_index(::Block) = throw(ArgumentError("Block must be converted by to_indices(...)"))
to_index(::BlockIndex) = throw(ArgumentError("BlockIndex must be converted by to_indices(...)"))
to_index(::BlockRange) = throw(ArgumentError("BlockRange must be converted by to_indices(...)"))

@inline to_indices(A, inds, I::Tuple{Block{1}, Vararg{Any}}) =
    (unblock(A, inds, I), to_indices(A, _maybetail(inds), tail(I))...)

# splat out higher dimensional blocks
# this mimics view of a CartesianIndex
@inline to_indices(A, inds, I::Tuple{Block, Vararg{Any}}) =
    to_indices(A, inds, (Block.(I[1].n)..., tail(I)...))

@inline to_indices(A, inds, I::Tuple{BlockRange{1,R}, Vararg{Any}}) where R =
    (unblock(A, inds, I), to_indices(A, _maybetail(inds), tail(I))...)

# splat out higher dimensional blocks
# this mimics view of a CartesianIndex
@inline to_indices(A, inds, I::Tuple{BlockRange, Vararg{Any}}) =
    to_indices(A, inds, (BlockRange.(tuple.(I[1].indices))..., tail(I)...))

# In 0.7, we need to override to_indices to avoid calling linearindices
@inline to_indices(A, I::Tuple{Block, Vararg{Any}}) =
    to_indices(A, axes(A), I)

@inline to_indices(A, I::Tuple{BlockRange, Vararg{Any}}) =
    to_indices(A, axes(A), I)

# The first argument for `reindex` is removed as of
# https://github.com/JuliaLang/julia/pull/30789 in Julia `Base`.  So,
# we define 2-arg `reindex` for Julia 1.2 and later.
if VERSION >= v"1.2-"

# BlockSlices map the blocks and the indices
# this is loosely based on Slice reindex in subarray.jl
reindex(idxs::Tuple{BlockSlice{<:BlockRange}, Vararg{Any}},
        subidxs::Tuple{BlockSlice{<:BlockRange}, Vararg{Any}}) =
    (@_propagate_inbounds_meta; (BlockSlice(BlockRange(idxs[1].block.indices[1][Int.(subidxs[1].block)]),
                                            idxs[1].indices[subidxs[1].indices]),
                                 reindex(tail(idxs), tail(subidxs))...))

reindex(idxs::Tuple{BlockSlice{BlockRange{1,Tuple{UnitRange{Int}}}}, Vararg{Any}},
        subidxs::Tuple{BlockSlice{Block{1,Int}}, Vararg{Any}}) =
    (@_propagate_inbounds_meta; (BlockSlice(Block(idxs[1].block.indices[1][Int(subidxs[1].block)]),
                                            idxs[1].indices[subidxs[1].indices]),
                                 reindex(tail(idxs), tail(subidxs))...))

function reindex(idxs::Tuple{BlockSlice{Block{1,Int}}, Vararg{Any}},
        subidxs::Tuple{BlockSlice{Block{1,Int}}, Vararg{Any}})
    (idxs[1], reindex(tail(idxs), tail(subidxs))...)
end

else  # if VERSION >= v"1.2-"

reindex(V, idxs::Tuple{BlockSlice{<:BlockRange}, Vararg{Any}},
        subidxs::Tuple{BlockSlice{<:BlockRange}, Vararg{Any}}) =
    (@_propagate_inbounds_meta; (BlockSlice(BlockRange(idxs[1].block.indices[1][Int.(subidxs[1].block)]),
                                            idxs[1].indices[subidxs[1].indices]),
                                    reindex(V, tail(idxs), tail(subidxs))...))

reindex(V, idxs::Tuple{BlockSlice{BlockRange{1,Tuple{UnitRange{Int}}}}, Vararg{Any}},
        subidxs::Tuple{BlockSlice{Block{1,Int}}, Vararg{Any}}) =
    (@_propagate_inbounds_meta; (BlockSlice(Block(idxs[1].block.indices[1][Int(subidxs[1].block)]),
                                            idxs[1].indices[subidxs[1].indices]),
                                    reindex(V, tail(idxs), tail(subidxs))...))

function reindex(V, idxs::Tuple{BlockSlice{Block{1,Int}}, Vararg{Any}},
        subidxs::Tuple{BlockSlice{Block{1,Int}}, Vararg{Any}})
    subidxs[1].block == Block(1) || throw(BoundsError(V, subidxs[1].block))
    (idxs[1], reindex(V, tail(idxs), tail(subidxs))...)
end

end  # if VERSION >= v"1.2-"

#################
# support for pointers
#################

const BlockOrRangeIndex = Union{RangeIndex, BlockSlice}

## BlockSlice1 is a convenience for views
const BlockSlice1 = BlockSlice{Block{1,Int},UnitRange{Int}}

block(A::BlockSlice) = block(A.block)
block(A::Block) = A

getblock(A::SubArray{<:Any,N,<:Any,NTuple{N,BlockSlice1}}) where N =
    getblock(parent(A), Int.(block.(parentindices(A)))...)
getblock(A::Adjoint, k::Int, j::Int) = getblock(parent(A), j, k)'
getblock(A::Transpose, k::Int, j::Int) = transpose(getblock(parent(A), j, k))

strides(A::SubArray{<:Any,N,<:BlockArray,NTuple{N,BlockSlice1}}) where N =
    strides(getblock(A))

for Adj in (:Transpose, :Adjoint)
    @eval strides(A::SubArray{<:Any,N,<:$Adj{<:Any,<:BlockArray},NTuple{N,BlockSlice1}}) where N =
        strides(getblock(A))
end

unsafe_convert(::Type{Ptr{T}}, V::SubArray{T, N, <:Any, <:NTuple{N, BlockSlice1}}) where {T,N} =
    unsafe_convert(Ptr{T}, getblock(parent(V), Int.(Block.(parentindices(V)))...))

unsafe_convert(::Type{Ptr{T}}, V::SubArray{T,N,PseudoBlockArray{T,N,AT},<:Tuple{Vararg{BlockOrRangeIndex}}}) where {T,N,AT} =
    unsafe_convert(Ptr{T}, V.parent) + (Base.first_index(V)-1)*sizeof(T)


text:
struct SerialExecution{T <: AbstractFloat,
                       A <: AbstractVector,
                       S <: LQSolver} <: AbstractExecution
    subproblems::Vector{SubProblem{T,A,S}}

    function SerialExecution(::Type{T}, ::Type{A}, ::Type{S}) where {T <: AbstractFloat, A <: AbstractVector, S <: LQSolver}
        return new{T,A,S}(Vector{SubProblem{T,A,S}}())
    end
end

function init_subproblems!(ph::AbstractProgressiveHedgingSolver, subsolver::QPSolver, execution::SerialExecution)
    for i = 1:ph.nscenarios
        push!(execution.subproblems,SubProblem(WS(ph.stochasticprogram, scenario(ph.stochasticprogram,i); solver = subsolver),
                                               i,
                                               probability(ph.stochasticprogram,i),
                                               decision_length(ph.stochasticprogram),
                                               subsolver))
    end
    update_iterate!(ph)
    return ph
end

function iterate!(ph::AbstractProgressiveHedgingSolver, ::SerialExecution)
    return iterate_nominal!(ph)
end

function init_workers!(::AbstractProgressiveHedgingSolver, ::SerialExecution)
    return nothing
end

function close_workers!(::AbstractProgressiveHedgingSolver, ::SerialExecution)
    return nothing
end

function resolve_subproblems!(ph::AbstractProgressiveHedgingSolver, execution::SerialExecution{T,A}) where {T <: AbstractFloat, A <: AbstractVector}
    Qs = A(undef, length(execution.subproblems))
    # Update subproblems
    reformulate_subproblems!(execution.subproblems, ph.ξ, penalty(ph))
    # Solve sub problems
    for (i, subproblem) ∈ enumerate(execution.subproblems)
        Qs[i] = subproblem()
    end
    # Return current objective value
    return sum(Qs)
end

function update_iterate!(ph::AbstractProgressiveHedgingSolver, execution::SerialExecution)
    # Update the estimate
    ξ_prev = copy(ph.ξ)
    ph.ξ[:] = sum([subproblem.π*subproblem.x for subproblem in execution.subproblems])
    # Update δ₁
    ph.data.δ₁ = norm(ph.ξ-ξ_prev, 2)^2
    return nothing
end

function update_subproblems!(ph::AbstractProgressiveHedgingSolver, execution::SerialExecution)
    # Update dual prices
    update_subproblems!(execution.subproblems, ph.ξ, penalty(ph))
    return nothing
end

function update_dual_gap!(ph::AbstractProgressiveHedgingSolver, execution::SerialExecution)
    # Update δ₂
    ph.data.δ₂ = sum([s.π*norm(s.x-ph.ξ,2)^2 for s in execution.subproblems])
    return nothing
end

function calculate_objective_value(ph::AbstractProgressiveHedgingSolver, execution::SerialExecution)
    return sum([get_objective_value(s) for s in execution.subproblems])
end

function fill_submodels!(ph::AbstractProgressiveHedgingSolver, scenarioproblems::StochasticPrograms.ScenarioProblems, execution::SerialExecution)
    for (i, submodel) in enumerate(scenarioproblems.problems)
        fill_submodel!(submodel, execution.subproblems[i])
    end
end

function fill_submodels!(ph::AbstractProgressiveHedgingSolver, scenarioproblems::StochasticPrograms.DScenarioProblems, execution::SerialExecution)
    j = 0
    @sync begin
        for w in workers()
            n = remotecall_fetch((sp)->length(fetch(sp).problems), w, scenarioproblems[w-1])
            for i in 1:n
                k = i+j
                @async remotecall_fetch((sp,i,x,μ,λ) -> fill_submodel!(fetch(sp).problems[i],x,μ,λ),
                                        w,
                                        scenarioproblems[w-1],
                                        i,
                                        get_solution(execution.subproblems[k])...)
            end
            j += n
        end
    end
end
# API
# ------------------------------------------------------------
struct Serial <: Execution end

function (execution::Serial)(::Type{T}, ::Type{A}, ::Type{S}) where {T <: AbstractFloat, A <: AbstractVector, S <: LQSolver}
    return SerialExecution(T,A,S)
end

function str(::Serial)
    return "serial execution"
end