text:
module IrisToolsF

    include("iris_variables.jl")
    include("iris_types.jl ")
    include("iris_functions.jl")

end

max_stars_repo_path: Section 2/iris_toolsf.jl

max_stars_repo_name: lytemar/Julia-for-Data-Science-Video

max_issues_repo_path: Section 2/iris_toolsf.jl

max_issues_repo_name: lytemar/Julia-for-Data-Science-Video

max_forks_repo_path: Section 2/iris_toolsf.jl

max_forks_repo_name: lytemar/Julia-for-Data-Science-Video


text:
A = [1 2 3; -4 -5 -6]

A[2,3]

length(A)

size(A,1)

size(A,2)

sumA = 0
for i = 1:size(A,1)
    for j = 1:size(A,2)
        sumA += A[i,j]
    end
end

println("Sum of the elements in A = ", sumA)

sumA = 0
for i = 1:size(A,1), j = 1:size(A,2)
    sumA += A[i,j]
end

println("Sum of the elements in A = ", sumA)

sumAsquared = 0
for a in A
    sumAsquared += a^2
end

print("Sum of the squares of the elements in A = ", sumAsquared)

C = ones(5,3)               # 5-by-4 matrix with all elements == 1

ones(5)      # 1-d array (can be used as a column vector)

ones(5,1)    # 2-d array, column vector

ones(1,5)    # 2-d array, row vector

D = [A; C]                     # Vertical concatenation - arrays must have same # columns
E = [zeros(5) C]               # Horizontal concatenation - arrays must have same # rows
F = [A ones(2,2); zeros(1,5)]  # General vertical-horizontal concatenation - all dimensions must match

G = [1:5 ones(Int64,5) 101:105]

B = A.^2 .- 3A           # Apply function elementwise

B = @. A^2 - 3A          # Easier syntax - entire expression evaluated elementwise

A[1,:]             # Row 1

A[:,1]             # Column 1

A[:,2:3]           # Columns 2-3

A[:,[1,3]] .= 0          # Set columns 1 and 3 to zero
A

A[:,2:3] = 2A[:,1:2]     # Set columns 2 and 3 to twice of columns 1 and 2
A

ext: jl


text:
using ProgressMeter
using Test
using LinearAlgebra,Statistics,StatsBase,Distributions
using QuadGK
using Plots,NamedColors ; theme(:dark) #; plotlyjs()
using FileIO
using BenchmarkTools
using FFTW
using Random
Random.seed!(0)

push!(LOAD_PATH, abspath(@__DIR__,".."))

using  HawkesSimulator; const global H = HawkesSimulator

##
function plotvs(x::AbstractArray{<:Real},y::AbstractArray{<:Real})
  x,y=x[:],y[:]
  @info """
  The max differences between the two are $(extrema(x .-y ))
  """
  plt=plot()
  scatter!(plt,x,y;leg=false,ratio=1,color=:white)
  lm=xlims()
  plot!(plt,identity,range(lm...;length=3);linestyle=:dash,color=:yellow)
  return plt
end

##

nneus = 2
tauker = 0.5
trace_ker = H.Trace(tauker,nneus)
trace_boh = H.Trace(123.0,nneus)

popstate = H.PopulationStateExpKernel(nneus,trace_ker,trace_boh)
myweights = [0.31 -0.3
            0.9  -0.15]
myinputs = [5.0 , 5.0]
rates_analytic  = inv(I-myweights)*myinputs

connection = H.ConnectionExpKernel(myweights,trace_ker)
population = H.PopulationExpKernel(popstate,connection,myinputs)

n_spikes = 10_000
recorder = H.RecFullTrain(n_spikes,1)
network = H.RecurrentNetworkExpKernel(population,recorder)

function simulate!(network,num_spikes)
  t_now = 0.0
  H.reset!(network) # clear spike trains etc
  for _ in 1:num_spikes
    t_now = H.dynamics_step!(t_now,network)
  end
  return t_now
end

Tmax = simulate!(network,n_spikes)
rates = H.numerical_rates(recorder,nneus,Tmax)

##
error("shtop!")

##

##

function generate_poisson_train(rate::R,Tend::R;Tstart=0.0) where R
  # preallocate for efficiency (even if it's already super fast)
  nmax = round(Integer,rate*Tend*1.5)
  ret = Vector{Float64}(undef,nmax)
  k = 0
  t_spike_now = 0.0
  while t_spike_now <= Tend
    k+=1
    t_spike_now += -log(rand())/rate # add delta t
    ret[k] = t_spike_now
  end
  deleteat!(ret,k+1:nmax) # remove unnecessary elements
  if !iszero(Tstart)
    ret .+= Tstart # shift by start time, if needed
  end
  return ret
end

function generate_perturbed_train1(rate::R,dissimilarity::R,
    reference_train::Vector{R}) where R
  @assert 0.0 <= dissimilarity <= 1.0
  new_train = copy(reference_train)
  nspikes = length(reference_train)
  nchange = round(Integer,dissimilarity*nspikes) # how many spikes to remove?
  if iszero(nchange)
    return new_train
  end
  # select spikes to delete, and delete them
  idx_change = sample(1:nspikes,nchange;replace=false,ordered=true)
  deleteat!(new_train,idx_change)
  rate_add = rate*dissimilarity
  Tend = nspikes/rate
  # add replacement spikes
  new_spikes = generate_poisson_train(rate_add,Tend)
  new_train = sort!(vcat(new_train,new_spikes))
  return new_train
end

function quickrate(train::Vector{R},Tend::R) where R
  k = findfirst(>=(Tend),train)
  if isnothing(k)
    error("please set Tend correctly!")
  end
  return (k-1)/Tend
end

##


text:
# Julia code
# Copyright (c) 2022 the SLICOTMath.jl developers
# Portions extracted from SLICOT-Reference distribution:
# Copyright (c) 2002-2020 NICONET e.V.

function run_mc01xd(datfile, io=stdout)
    NIN = 5
    NOUT = 6
    NEV = 3
    LDWORK = 42
    ZERO = 0.0e0
    EVR = Array{Float64,1}(undef, NEV)
    EVI = Array{Float64,1}(undef, NEV)
    EVQ = Array{Float64,1}(undef, NEV)
    DWORK = Array{Float64,1}(undef, LDWORK)
    RT = Array{Float64,1}(undef, 2)
    RTS = Array{ComplexF64,1}(undef, 3)
    f = open(datfile,"r")
    readline(f)
    vs = split(readline(f))
    ALPHA = parse(Float64, replace(vs[1],'D'=>'E'))
    BETA = parse(Float64, replace(vs[2],'D'=>'E'))
    GAMMA = parse(Float64, replace(vs[3],'D'=>'E'))
    DELTA = parse(Float64, replace(vs[4],'D'=>'E'))
# interp call 1

    INFO = SLICOT.mc01xd!(ALPHA, BETA, GAMMA, DELTA, EVR, EVI, EVQ)
    @test INFO == 0
    INFO == 0 || return

if ( INFO!=0 )
else
# interp output 1
    println(io,"EVR:")
    _nr = 3
    show(io,"text/plain",EVR[1:_nr])
    println(io,)

# interp output 2
    println(io,"EVI:")
    _nr = 3
    show(io,"text/plain",EVI[1:_nr])
    println(io,)

# interp output 3
    println(io,"EVQ:")
    _nr = 3
    show(io,"text/plain",EVQ[1:_nr])
    println(io,)

if ( true )
    RTS = (EVR + im * EVI) ./ EVQ
# interp output 4
    println(io,"RTS:")
    _nr = 3
    show(io,"text/plain",RTS[1:_nr])
    println(io,)

end # if
end # if
    close(f)
end # run_X()


text:
function adder(x, y)
    a = x + y
    return a
end

max_stars_repo_name: wigging/julia-computing

max_issues_repo_name: wigging/julia-computing

max_forks_repo_name: wigging/julia-computing


text:
import AutomotiveDrivingModels.observe!

mutable struct MOBILDriver <: LaneChangeModel

    dir::Int
    #rec::SceneRecord
    mlon::LaneFollowingDriver
    safe_decel::Float64 # safe deceleration (positive value)
    politeness::Float64 # politeness factor (suggested p ∈ [0.2,0.5])
    advantage_threshold::Float64 # Δaₜₕ
    max_horizon::Float64

    function MOBILDriver(
        timestep::Float64;
        #rec::SceneRecord=SceneRecord(2,timestep),
        mlon::LaneFollowingDriver=IDMDriver(),
        safe_decel::Float64=2.0, # [m/s²]
        politeness::Float64=0.35,
        advantage_threshold::Float64=0.1,
        max_horizon::Float64=35.0,
        )

        retval = new()
        retval.dir = DIR_MIDDLE
        #retval.rec = rec
        retval.mlon = mlon
        retval.safe_decel = safe_decel
        retval.politeness = politeness
        retval.advantage_threshold = advantage_threshold
        retval.max_horizon = max_horizon
        retval
    end
end
get_name(::MOBILDriver) = "MOBILDriver"
function set_desired_speed!(model::MOBILDriver, v_des::Float64)
    set_desired_speed!(model.mlon, v_des)
    model
end
function AutomotiveDrivingModels.observe!(model::MOBILDriver, scene::Scene, roadway::Roadway, egoid::Int)
    observe_helper!(model,scene,roadway,egoid)
end
function AutomotiveDrivingModels.observe!(model::MOBILDriver, scene::Frame{Entity{VehicleState, BicycleModel, Int}}, roadway::Roadway, egoid::Int)
    observe_helper!(model,scene,roadway,egoid)
end

function observe_helper!(model::MOBILDriver, scene::Union{Scene,Frame{Entity{VehicleState, BicycleModel, Int}}}, roadway::Roadway, egoid::Int)
    #rec = model.rec
    #update!(rec, scene)

    #vehicle_index = findfirst(rec[0], egoid)
    vehicle_index = egoid
    veh_ego = scene[vehicle_index]
    v = veh_ego.state.v
    egostate_M = veh_ego.state

    #left_lane_exists = convert(Float64, get(N_LANE_LEFT, rec, roadway, vehicle_index)) > 0
    #right_lane_exists = convert(Float64, get(N_LANE_RIGHT, rec, roadway, vehicle_index)) > 0
    lane = roadway[scene[egoid].state.posF.roadind.tag]
    left_lane_exists = (n_lanes_left(lane, roadway) > 0 )
    right_lane_exists = (n_lanes_right(lane, roadway) > 0 )
    fore_M = get_neighbor_fore_along_lane(scene, vehicle_index, roadway, VehicleTargetPointFront(), VehicleTargetPointRear(), VehicleTargetPointFront(),max_distance_fore=model.max_horizon)
    rear_M = get_neighbor_rear_along_lane(scene, vehicle_index, roadway, VehicleTargetPointFront(), VehicleTargetPointFront(), VehicleTargetPointRear(),max_distance_rear=model.max_horizon)

    # accel if we do not make a lane change
    accel_M_orig = rand(observe!(reset_hidden_state!(model.mlon), scene, roadway, egoid)).a
    model.dir = DIR_MIDDLE

    advantage_threshold = model.advantage_threshold
    if right_lane_exists

        rear_R = get_neighbor_rear_along_right_lane(scene, vehicle_index, roadway, VehicleTargetPointFront(), VehicleTargetPointFront(), VehicleTargetPointRear(),max_distance_rear=model.max_horizon)

        # candidate position after lane change is over
        footpoint = get_footpoint(convert2vehicle(veh_ego))
        lane = roadway[veh_ego.state.posF.roadind.tag]
        lane_R = roadway[LaneTag(lane.tag.segment, lane.tag.lane - 1)]
        roadproj = proj(footpoint, lane_R, roadway)
        frenet_R = Frenet(RoadIndex(roadproj), roadway)
        egostate_R = VehicleState(frenet_R, roadway, veh_ego.state.v)

        Δaccel_n = 0.0
        passes_safety_criterion = true
        if rear_R.ind != 0
            id = scene[rear_R.ind].id
            accel_n_orig = rand(observe!(reset_hidden_state!(model.mlon), scene, roadway, id)).a
            veh_ego = Entity(veh_ego, egostate_R)
            scene[vehicle_index] = veh_ego
            accel_n_test = rand(observe!(reset_hidden_state!(model.mlon), scene, roadway, id)).a

            # body = inertial2body(get_rear(convert2vehicle(scene[rear_R.ind])), get_front(convert2vehicle(veh_ego))) # project target to be relative to ego
            body = inertial2body(get_rear(convert2vehicle(veh_ego)), get_front(convert2vehicle(scene[rear_R.ind])))
            s_gap = body.x

            veh_ego = Entity(veh_ego, egostate_M)
            scene[vehicle_index] = veh_ego
            passes_safety_criterion = accel_n_test ≥ -model.safe_decel && s_gap ≥ 0
            Δaccel_n = accel_n_test - accel_n_orig
        end

        if passes_safety_criterion

            Δaccel_o = 0.0
            if rear_M.ind != 0
                id = scene[rear_M.ind].id
                accel_o_orig = rand(observe!(reset_hidden_state!(model.mlon), scene, roadway, id)).a
                veh_ego = Entity(veh_ego, egostate_R)
                scene[vehicle_index] = veh_ego
                accel_o_test = rand(observe!(reset_hidden_state!(model.mlon), scene, roadway, id)).a
                veh_ego = Entity(veh_ego, egostate_M)
                scene[vehicle_index] = veh_ego
                Δaccel_o = accel_o_test - accel_o_orig
            end

            veh_ego = Entity(veh_ego, egostate_R)
            scene[vehicle_index] = veh_ego
            accel_M_test = rand(observe!(reset_hidden_state!(model.mlon), scene, roadway, egoid)).a
            veh_ego = Entity(veh_ego, egostate_M)
            scene[vehicle_index] = veh_ego
            Δaccel_M = accel_M_test - accel_M_orig

            Δaₜₕ = Δaccel_M + model.politeness*(Δaccel_n + Δaccel_o)
            if Δaₜₕ > advantage_threshold
                model.dir = DIR_RIGHT
                advantage_threshold = Δaₜₕ
            end
        end
    end

    if left_lane_exists
        rear_L = get_neighbor_rear_along_left_lane(scene, vehicle_index, roadway, VehicleTargetPointFront(), VehicleTargetPointFront(), VehicleTargetPointRear(),max_distance_rear=model.max_horizon)

        # candidate position after lane change is over
        footpoint = get_footpoint(convert2vehicle(veh_ego))
        lane = roadway[veh_ego.state.posF.roadind.tag]
        lane_L = roadway[LaneTag(lane.tag.segment, lane.tag.lane + 1)]
        roadproj = proj(footpoint, lane_L, roadway)
        frenet_L = Frenet(RoadIndex(roadproj), roadway)
        egostate_L = VehicleState(frenet_L, roadway, veh_ego.state.v)

        Δaccel_n = 0.0
        passes_safety_criterion = true
        if rear_L.ind != 0
            id = scene[rear_L.ind].id
            accel_n_orig = rand(observe!(reset_hidden_state!(model.mlon), scene, roadway, id)).a
            veh_ego = Entity(veh_ego, egostate_L)
            scene[vehicle_index] = veh_ego
            accel_n_test = rand(observe!(reset_hidden_state!(model.mlon), scene, roadway, id)).a
            # body = inertial2body(get_rear(convert2vehicle(scene[rear_L.ind])), get_front(convert2vehicle(veh_ego))) # project target to be relative to ego
            body = inertial2body(get_rear(convert2vehicle(veh_ego)), get_front(convert2vehicle(scene[rear_L.ind])))
            s_gap = body.x
            veh_ego = Entity(veh_ego, egostate_M)
            scene[vehicle_index] = veh_ego
            passes_safety_criterion = accel_n_test ≥ -model.safe_decel && s_gap ≥ 0
            Δaccel_n = accel_n_test - accel_n_orig
        end

        if passes_safety_criterion

            Δaccel_o = 0.0
            if rear_M.ind != 0
                id = scene[rear_M.ind].id
                accel_o_orig = rand(observe!(reset_hidden_state!(model.mlon), scene, roadway, id)).a
                veh_ego = Entity(veh_ego, egostate_L)
                scene[vehicle_index] = veh_ego
                accel_o_test = rand(observe!(reset_hidden_state!(model.mlon), scene, roadway, id)).a
                veh_ego = Entity(veh_ego, egostate_M)
                scene[vehicle_index] = veh_ego
                Δaccel_o = accel_o_test - accel_o_orig
            end

            veh_ego = Entity(veh_ego, egostate_L)
            scene[vehicle_index] = veh_ego
            accel_M_test = rand(observe!(reset_hidden_state!(model.mlon), scene, roadway, egoid)).a
            veh_ego = Entity(veh_ego, egostate_M)
            scene[vehicle_index] = veh_ego
            Δaccel_M = accel_M_test - accel_M_orig

            Δaₜₕ = Δaccel_M + model.politeness*(Δaccel_n + Δaccel_o)
            if Δaₜₕ > advantage_threshold
                model.dir = DIR_LEFT
                advantage_threshold = Δaₜₕ
            end
        end
    end

    # chack if front vehicle is doing the same action
    fore = get_neighbor_fore_along_lane(scene, vehicle_index, roadway, VehicleTargetPointFront(), VehicleTargetPointRear(), VehicleTargetPointFront(),max_distance_fore=model.max_horizon)
    if fore.ind != 0
        if model.dir == DIR_LEFT && scene[fore.ind].state.posF.t > 0.0 && sin(scene[fore.ind].state.posF.ϕ) > 0.0
            model.dir == DIR_MIDDLE
        elseif model.dir == DIR_RIGHT && scene[fore.ind].state.posF.t < 0.0 && sin(scene[fore.ind].state.posF.ϕ) < 0.0
            model.dir == DIR_MIDDLE
        end
    end

    model
end

Base.rand(model::MOBILDriver) = LaneChangeChoice(model.dir)


text:
using MyGraph
using DataStructures

G = DictGraph(true)
B = Array{OrderedDict, 1}()
ctr = 0

addEdge!(G,0,2)
addEdge!(G,0,4)
addEdge!(G,0,3)
addEdge!(G,1,1)
addEdge!(G,1,4)
addEdge!(G,1,2)
addEdge!(G,2,4)
addEdge!(G,3,4)
addEdge!(G,3,5)
addEdge!(G,4,5)
addEdge!(G,5,1)
addEdge!(G,5,6)

for v in vertices(G)                        #para cada vertice em G
    push!(B, BFS(G,v))                      #inserir em B o BSF a partir do vertive V
end

for D in B
    local R = Array{Int, 1}()               #Inicializa R
    R = [0 for _ = 1:length(G.relation)]    #Vetor R do tamando da quantidade de verices no grafo.
    for k in keys(D)
        R[k + 1] = D[k][2]
    end
    println("$ctr => $R")
    global ctr += 1
end


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

"""
    polymat(xs, d)

Return the matrix of monomials for the iterator `xs`, i.e.
for each item `x = (x₁, x₂,…, xₙ)` in `xs`, evaluate
the monomial terms of the expansion `(x₁ + x₂ + ⋯ + xₙ)ᵈ`
for a given degree `d`.

The resulting matrix has a number of rows that is equal
to the number of items in the iterator `xs`. The number
of columns is a function of the degree. For `d=0`, a
single column of ones is returned that corresponds to
the constant term `x₁⁰⋅x₂⁰⋅⋯⋅xₙ⁰` for all items in `xs`.
"""
function polymat(xs, d)
  x  = first(xs)
  n  = length(x)
  es = Iterators.flatten(multiexponents(n, d) for d in 0:d)
  ps = [[prod(x.^e) for x in xs] for e in es]
  reduce(hcat, ps)
end

"""
    trend(data, vars; degree=1)

Return the deterministic spatial trend for the variables `vars`
in the spatial `data`. Approximate the trend with a polynomial
of given `degree`.

## References

* Menafoglio, A., Secchi, P. 2013. [A Universal Kriging predictor
  for spatially dependent functional data of a Hilbert Space]
  (https://doi.org/10.1214/13-EJS843)
"""
function trend(data, vars::AbstractVector{Symbol}; degree=1)
  𝒯 = values(data)
  𝒟 = domain(data)

  # build polynomial drift terms
  coords(𝒟, i) = coordinates(centroid(𝒟, i))
  xs = (coords(𝒟, i) for i in 1:nelements(𝒟))
  F  = polymat(xs, degree)

  # eqs 25 and 26 in Menafoglio, A., Secchi, P. 2013.
  ms = map(vars) do var
    z  = Tables.getcolumn(𝒯, var)
    a  = (F'F \ F') * z
    F * a
  end

  ctor  = Tables.materializer(𝒯)
  means = ctor((; zip(vars, ms)...))

  georef(means, 𝒟)
end

trend(data, var::Symbol; kwargs...) =
  trend(data, [var]; kwargs...)

function trend(data; kwargs...)
  𝒯 = values(data)
  s = Tables.schema(𝒯)
  vars = collect(s.names)
  trend(data, vars; kwargs...)
end


text:
#=
# Variable struct and related functions
=#

"""
    Variable

Represents a value that is not defined by an independent function of coordinates.
Note that this does not have to be an unknown that is solved for. It can hold
values that are dependent on other variables such as transformed variables, or
any other value that the user has control over.

This is also useful for setting up data arrays that match those of unknowns.
The values for a variable v can be directly accessed as `v.values` and has
dimensions [C, N] where C is the number of components (SCALAR=1, etc.) and N is
the number of nodes or cells.

This should be built with the `variable` function.
"""
mutable struct Variable
    symbol::Symbol          # symbol used in expressions referring to this variable
    symvar::Array{Basic,1}  # Array of symbolic layer variable symbols(Basic symbols)
    index::Int              # index in the Finch list of variables
    type::String            # constants for SCALAR, VECTOR, etc.
    location::String        # constant for NODAL or CELL
    discretization::String  # constant for CG, DG, FV
    values::Array{Float64}  # a C x N array, C is number of components (SCALAR=1, etc.)

    # These are mainly for ARRAY type vars. Should there be a separate struct for those?
    indexer                  # indexer object or array of them
    total_components::Int   # Number of components
    #transposed::Bool        # If true, the values are actually trasposed -> N x C array (better for large C?)

    # not currently used, may be removed
    dependson               # a list of variables that this depends on
    ready::Bool             # Is this variable's values ready? Can dependent variables use it?
end

"""
    VariableTransform

Two sets of variables and a function that transforms one into the other.
This is built with the `variableTransform` function.
"""
struct VariableTransform
    from::Union{Variable, Array{Variable}}  # Transformed from this
    to::Union{Variable, Array{Variable}}    # to this
    func                                    # using this function
    # This function takes one number(or array matching from)
    # and returns one number(or and array matching to)
    # NOTE: if "from" is an array, "to" must also be an array (size can be different)
end

# For printing, write the variable symbol only
Base.show(io::IO, x::Variable) = print(io, string(x.symbol));

Base.:+(x::Finch.Variable, y::Finch.Variable) = x.values + y.values;
Base.:-(x::Finch.Variable, y::Finch.Variable) = x.values - y.values;
Base.:*(x::Finch.Variable, y::Finch.Variable) = x.values .* y.values;
Base.:/(x::Finch.Variable, y::Finch.Variable) = x.values ./ y.values;
Base.:^(x::Finch.Variable, y::Number) = x.values .^ y;

# Performs the transform on all of the variable values.
# Writes the results directly into "to".
function transform_variable_values(xform)
    is_array = false;
    if typeof(xform.from) <: Array
        from_len = length(xform.from);
        N = length(xform.from[1].values);
        is_array = true;
    else
        from_len = 1;
        N = length(xform.from.values);
    end
    if typeof(xform.to) <: Array
        to_len = length(xform.to);
    else
        to_len = 1;
    end

    if is_array
        input = zeros(from_len);
        output = zeros(to_len);
        for i=1:N
            for j=1:from_len
                input[j] = xform.from[j].values[i];
            end
            output = xform.func(input);
            for j=1:to_len
                xform.to[j].values[i] = output[j];
            end
        end
    else
        for i=1:N
            xform.to.values[i] = xform.func(xform.from.values[i]);
        end
    end
    return nothing;
end

max_stars_repo_name: paralab/Finch

max_issues_repo_name: paralab/Finch

max_forks_repo_name: paralab/Finch


text:
function create_geotiff(img_filename, outfile=replace(img_filename, ".png" => ".tif");
                        demrsc=nothing, demrscfile=nothing)

    if isnothing(demrsc)
        demrsc = isnothing(demrscfile) ? Sario.load(Sario.find_rsc_file(".")) : Sario.load(demrscfile)
    end
    # gdalwarp -s_srs EPSG:4326 -t_srs EPSG:3857 out1.tif out2.tif
    # ullr means upper left, lower right (lon, lat) points
    north, south, east, west = nsew(demrsc)
    cmd = `gdal_translate -of GTiff -a_nodata 0 -a_srs EPSG:4326 -a_ullr $west $north $east $south $img_filename $outfile`
    println("Running:")
    println(cmd)
    run(cmd)
end

"""return tuple of (north, south, east, west) from demrsc data"""
function nsew(demrsc)
    @unpack x_first, y_first, x_step, y_step, width, file_length = demrsc
    north = y_first
    west = x_first
    east = west + width * x_step
    south = north + file_length * y_step
    return north, south, east, west
end


text:
#=
This file is auto-generated. Do not edit.
=#
"""
    mutable struct SEXS <: AVR
        Ta_Tb::Float64
        Tb::Float64
        K::Float64
        Te::Float64
        V_lim::NamedTuple{(:min, :max), Tuple{Float64, Float64}}
        ext::Dict{String, Any}
        states::Vector{Symbol}
        n_states::Int
        states_types::Vector{StateTypes}
        internal::InfrastructureSystemsInternal
    end

Parameters of Simplified Excitation System Model - SEXS in PSSE

# Arguments
- `Ta_Tb::Float64`: Ratio of lead and lag time constants, validation range: `(0, nothing)`
- `Tb::Float64`: Lag time constant, validation range: `(eps(), nothing)`
- `K::Float64`: Gain, validation range: `(0, nothing)`
- `Te::Float64`: Field circuit time constant in s, validation range: `(0, nothing)`
- `V_lim::NamedTuple{(:min, :max), Tuple{Float64, Float64}}`: Field voltage limits
- `ext::Dict{String, Any}`
- `states::Vector{Symbol}`: The states are:	Vf: Voltage field,	Vr: Lead-lag state
- `n_states::Int`: SEXS has 2 states
- `states_types::Vector{StateTypes}`: SEXS has 2 differential states
- `internal::InfrastructureSystemsInternal`: power system internal reference, do not modify
"""
mutable struct SEXS <: AVR
    "Ratio of lead and lag time constants"
    Ta_Tb::Float64
    "Lag time constant"
    Tb::Float64
    "Gain"
    K::Float64
    "Field circuit time constant in s"
    Te::Float64
    "Field voltage limits"
    V_lim::NamedTuple{(:min, :max), Tuple{Float64, Float64}}
    ext::Dict{String, Any}
    "The states are:	Vf: Voltage field,	Vr: Lead-lag state"
    states::Vector{Symbol}
    "SEXS has 2 states"
    n_states::Int
    "SEXS has 2 differential states"
    states_types::Vector{StateTypes}
    "power system internal reference, do not modify"
    internal::InfrastructureSystemsInternal
end

function SEXS(Ta_Tb, Tb, K, Te, V_lim, ext=Dict{String, Any}(), )
    SEXS(Ta_Tb, Tb, K, Te, V_lim, ext, [:Vf, :Vr], 2, [StateTypes.Differential, StateTypes.Differential], InfrastructureSystemsInternal(), )
end

function SEXS(; Ta_Tb, Tb, K, Te, V_lim, ext=Dict{String, Any}(), states=[:Vf, :Vr], n_states=2, states_types=[StateTypes.Differential, StateTypes.Differential], internal=InfrastructureSystemsInternal(), )
    SEXS(Ta_Tb, Tb, K, Te, V_lim, ext, states, n_states, states_types, internal, )
end

# Constructor for demo purposes; non-functional.
function SEXS(::Nothing)
    SEXS(;
        Ta_Tb=0,
        Tb=0,
        K=0,
        Te=0,
        V_lim=(min=0.0, max=0.0),
        ext=Dict{String, Any}(),
    )
end

"""Get [`SEXS`](@ref) `Ta_Tb`."""
get_Ta_Tb(value::SEXS) = value.Ta_Tb
"""Get [`SEXS`](@ref) `Tb`."""
get_Tb(value::SEXS) = value.Tb
"""Get [`SEXS`](@ref) `K`."""
get_K(value::SEXS) = value.K
"""Get [`SEXS`](@ref) `Te`."""
get_Te(value::SEXS) = value.Te
"""Get [`SEXS`](@ref) `V_lim`."""
get_V_lim(value::SEXS) = value.V_lim
"""Get [`SEXS`](@ref) `ext`."""
get_ext(value::SEXS) = value.ext
"""Get [`SEXS`](@ref) `states`."""
get_states(value::SEXS) = value.states
"""Get [`SEXS`](@ref) `n_states`."""
get_n_states(value::SEXS) = value.n_states
"""Get [`SEXS`](@ref) `states_types`."""
get_states_types(value::SEXS) = value.states_types
"""Get [`SEXS`](@ref) `internal`."""
get_internal(value::SEXS) = value.internal

"""Set [`SEXS`](@ref) `Ta_Tb`."""
set_Ta_Tb!(value::SEXS, val) = value.Ta_Tb = val
"""Set [`SEXS`](@ref) `Tb`."""
set_Tb!(value::SEXS, val) = value.Tb = val
"""Set [`SEXS`](@ref) `K`."""
set_K!(value::SEXS, val) = value.K = val
"""Set [`SEXS`](@ref) `Te`."""
set_Te!(value::SEXS, val) = value.Te = val
"""Set [`SEXS`](@ref) `V_lim`."""
set_V_lim!(value::SEXS, val) = value.V_lim = val
"""Set [`SEXS`](@ref) `ext`."""
set_ext!(value::SEXS, val) = value.ext = val
"""Set [`SEXS`](@ref) `states_types`."""
set_states_types!(value::SEXS, val) = value.states_types = val


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

module REPLTests
    include("repl.jl")
end
module REPLCompletionsTest
    include("replcompletions.jl")
end
module LineEditTest
    include("lineedit.jl")
end
module DocviewTest
    include("docview.jl")
end
module TerminalMenusTest
    include("TerminalMenus/runtests.jl")
end

max_stars_repo_name: vanillajonathan/julia

max_issues_repo_name: vanillajonathan/julia

max_forks_repo_name: vanillajonathan/julia


text:
const flag_list = Any[ :!, :assert,
    :cat, :glue, :lazy, :nolazy,
    :strided, :julienne, :batched,
    :named, :nameddims, :namedarray, :axis ]

#= TODO

# this gives an error -- change to sizeinfer! which updates dict, and call that?
# or perhaps just check at :staticslice & throw an error
@pretty @cast A[k]{i,j} := B[i,(j,k)]  k:length(C)

# I thought this was broken for a bit, but now looks OK, add to tests?
@pretty @cast A[i,j] := B[i][-j]

# named...
B = rand(2,3);
@cast A[i,j] := B[j,i] named
@reduce D[i] := sum(j)  A[i,j] named

=#

struct SizeDict
    dict::Dict{Any, Any}
    checks::Vector{Any}
    seen::Vector{Any}
    topex::Vector{Any}
    rightnames::Vector{Any}
end

SizeDict() = SizeDict(Dict{Any, Any}(), Any[], Any[], Any[], Any[])

#="""
    flat, getafix, putsz, negated = parse!(sdict, A, outer, inner)

Use this for `A[outer...][inner...]`.
* `flat` is a list of symbols, naming the indices, in order.
* `negated` is the subset which had - in front.
* `getafix` is the thing for making views of given input or target array, : in non-fixed directions.
* `putsz` is a tuple of `sz[n]` and products, which you want for reshaping to the very final array;
  the numbers `n` refer to position in `flat`.
* `sdict::SizeDict` ends up with `sdict.dict[:i] = ` an expr for length of this index.
  To avoid saving these (if you don't have A yet) set `A = nothing`.

    parse!(sdict, A, outer, reduce_inner; allowranges=true)

This is for `A[outer...] = sum(inner...) ...` LHS, which are allowed to have sum(a:2, ...) ranges.
The allowranges flag *also* disables the use of `size(first(A), 2)` etc,
to disable `size(A,2)` set `A = nothing` again.

    parse!(sdict, nothing, [], rex; allowranges=true, flags=vector)

For dimensions & annotations.
Will look for flags because it was given something to push them into.
You should now do this first, so that `sdict` is most likely to have these (neater) entries.
They are added using savesize!(sdict,...) which now puts later sizes into list of checks.
"""=#
function parse!(sdict::SizeDict, A, outer, inner; allowranges=false, flags=nothing)

    flat = Any[]
    getafix = Any[]
    putsz = Any[]
    negated = Any[]

    ### INNER dimensions come first in flat
    # also used for options / reduced dimensions
    for (din, ind) in enumerate(inner)

        # look for flags only in rex case
        if flags != nothing && ind in flag_list # global const
            push!(flags, ind)
            continue
        end

        # in rex or sum(a, b:3, c), some indices may have explicit ranges
        if allowranges

            if @capture(ind , i_:ilength_ )
                savesize!(sdict, i, ilength)
                if flags != nothing
                    push!(flags, :assert) # turn on size checks if you can
                else
                    # for sum(i:3, j) can't see flags, so can't do that, sorry
                end
            else
                i = ind    # replace nothing from @capture
                # i = stripminus!(negated, ind) # now i isa Symbol # this gives errors
            end
            # TODO should I have stripminus! here? @pretty @cast A[i,j] := B[i][-j]
            push!(flat, i) # all need to go into flat!

        # for simple slicing, always want iflat, but size(first(A),d) only if we know A.
        else
            i = stripminus!(negated, ind) # now i isa Symbol
            push!(flat, i)

            if !isnothing(A) # only save these if we can access A.
                savesize!(sdict, i, :( size(first($A), $din) ) )
            end
        end
    end

    ### OUTER dimensions never allow flags or ranges,
    # but do allow tuples etc, and fixed indices.
    for (dout, ind) in enumerate(outer)
        if isconstant(ind)
            # if ind == :_  ...         # treat _ almost like 1, now handled by needview! / rview
            if ind isa Expr             # constant slice with $c
                @assert ind.head == :($)
                ind = ind.args[1]
            end
            push!(getafix, ind)         # for view(A, getafix...) if RHS / in-place
            push!(putsz, 1)             # or if driving the other direction, make size=1.

        else
            push!(getafix, Colon())
            ind = stripminus!(negated, ind) # this looks inside for -, and returns Symbol or Vector{Symbol}
            if !isnothing(A)
                savesize!(sdict, ind, :( size($A, $dout) ) )
            end

            push!(putsz, szwrap(ind) ) # this may be product  sz_i * sz_j  or just one
            push_or_append!(flat, ind)
        end
    end

    if !isnothing(A) && length(outer)>0
        N = length(outer)
        str = "expected a $N-tensor $A[" * join(outer, ", ") * "]"
        push!(sdict.checks, :(TensorCast.@assert_ ndims($A)==$N $str) )
    end

    return flat, getafix, putsz, negated
end

isconstant(i::Int) = true
isconstant(i::Symbol) = i == :_
isconstant(ex::Expr) = ex.head == :($)

function findcheck(i::Symbol, ind::Vector, where=nothing)
    d = findfirst(isequal(i), ind)
    d isa Nothing && throw(MacroError("can't find index $i", where))
    return d
end
findcheck(i::Expr, ind::Vector, where=nothing) =
    throw(MacroError("don't know what to do with index $i, expected a single symbol", where))

function savesize!(store::SizeDict, ind, long)
    if !haskey(store.dict, ind) # then we don't yet have a size for this index (or tuple)
        store.dict[ind] = long
    else                        # we have seen it before,
        if isa(ind, Symbol)     # but list of checks is only for actual indices
            str = "range of index $ind must agree"
            push!(store.checks, :(TensorCast.@assert_ $(store.dict[ind]) == $long $str) )
        end
    end
end

push_or_append!(list, ind::Symbol) = push!(list, ind)
push_or_append!(list, inds::Vector) = append!(list, inds)

szwrap(i::Symbol) = Symbol(:sz_,i)
szwrap(ijk::Vector) = :( TensorCast.star($([ Symbol(:sz_,i) for i in ijk ]...)) )

#="""
    stripminus(negated, i)

For any index, or vector of indices, or vector of indices containing tuples / backslash combos,
this pushes every letter with a minus in front into negated list.

It returns an index, or a tidied-up vector of indices, in which tuples etc have become vectors.
"""=#
stripminus!(negated::Vector, ind::Symbol) = ind  # get one Symol & you're done

stripminus!(negated::Vector, ind::Int) = throw(MacroError("can't handle fixed index $ind here"))

function stripminus!(negated::Vector, ind::Expr)
    # minus one index
    if @capture(ind, -i_Symbol)
        push!(negated, i)       # add that to the negated list
        return i                # and go home.
    elseif @capture(ind, (i_Symbol)' )
        return Symbol(i, '′')   # normalise i' to be i′
        # return Symbol(i, "_prime") # allow i' but keep distinct?

    # tuple notation
    elseif @capture(ind, -(ijk__,) )
        ijkplus = stripminus!(negated, ijk)
        append!(negated, ijkplus)

    elseif @capture(ind, (ijk__,) )

    # backslash notation -- if you want more than four \ then I'm impressed! TODO make recursive
    elseif @capture(ind, i_\j_\k_\l_\m )
        ijk = Any[i,j,k,l,m]
    elseif @capture(ind, i_\j_\k_\l_ )
        ijk = Any[i,j,k,l]
    elseif @capture(ind, i_\j_\k_ )
        ijk = Any[i,j,k]
    elseif @capture(ind, i_\j_ )
        ijk = Any[i,j]

    elseif @capture(ind, i_⊗j_⊗k_⊗l_⊗m )
        ijk = Any[i,j,k,l,m]
    elseif @capture(ind, i_⊗j_⊗k_⊗l_ )
        ijk = Any[i,j,k,l]
    elseif @capture(ind, i_⊗j_⊗k_ )
        ijk = Any[i,j,k]
    elseif @capture(ind, i_⊗j_ )
        ijk = Any[i,j]

    else
        throw(MacroError("stuck on $ind, doesn't look like a valid index, or set of indices"))
    end

    # now we have a vector ijk of symbols, or expressions :(-i), no more
    return Any[ stripminus!(negated, i) for i in ijk ]
end

#="""
    oddunique(negated)

Returns a list in which anything repeated evenly many times has been removed, then `unique`.
"""=#
function oddunique(negated)  # -1 * -1 = +1
    set = Set{Symbol}()
    for i in negated
        if i in set
            pop!(set, i)
        else
            push!(set, i)
        end
    end
    Any[ i for i in set ]
end

#="""
    checkrepeats(flat)

Throws an error if there are repeated indices.
Also a good place to check whether flat is flat... e.g. A[i][j\\k] will arrive here
"""=#
function checkrepeats(flat, msg="", where=nothing)
    once = Set{Symbol}()
    twice = Set{Symbol}()
    for i in flat
        if i != nothing
            isa(i, Symbol) || throw(MacroError("can't handle index $i, expected a single symbol" * msg, where))

            if i in once
                push!(twice, i)
            else
                push!(once, i)
            end
        end
    end
    if length(twice) > 0
        str = join(twice, ", ")
        throw(MacroError("index $str repeated" * msg, where))
    end
    nothing
end

#="""
    sizeinfer(store, icannon, where, leaveone=true)

This is the point of SizeDict.
The goal is to produce a canonical vector of sizes, corresponding to vector of symbols icannon.
If these sizes are known, easy!

But for unknown ones, we do a second pass, looking for entries in sizedict like [:i, :j]
which come from tuples of indices, for which we know the product of their dimensions.
"""=#
function sizeinfer(store::SizeDict, icanon::Vector, where=nothing, leaveone = true)
    sizes = Any[ (:) for i in icanon ]

    # first pass looks for single indices whose lengths are known directly
    for pair in store.dict
        if isa(pair.first, Symbol)
            d = findcheck(pair.first, icanon, where)
            sizes[d] = pair.second  # so write it into the answer
        end
    end

    if leaveone && count(isequal(:), sizes) < 2
        V && count(isequal(:), sizes)==1 &&  @info "sizeinfer -- leaving one : in sizes" repr(sizes)
        return sizes
    end
    V && @info "sizeinfer -- first pass" sizes

    # second pass looks for tuples (whose product-length is known) where exactly one entry has unknown length
    for pair in store.dict
        if isa(pair.first, Vector)
            known = haskey.(Ref(store.dict), pair.first)

            if sum(.!known) == 1    # bingo! now work out its size:
                num = pair.second
                denfacts = [ store.dict[i] for i in pair.first[known] ]
                if length(denfacts) > 1
                    den = :( *($(denfacts...)) )
                else
                    den = :( $(denfacts[1]) )
                end
                rat = :( $num ÷ $den )

                d = findcheck(first(pair.first[.!known]), icanon, where)
                sizes[d] = rat      # save that size

                str = "inferring range of $(icanon[d]) from range of $(join(pair.first, " ⊗ "))"
                push!(store.checks, :( TensorCast.@assert_ rem($num, $den)==0 $str) )
            end
        end
    end
    V && @info "sizeinfer -- second pass" sizes

    unknown = Any[ i for (d,i) in enumerate(icanon) if sizes[d] == (:) ]
    str = join(unknown, ", ")
    length(unknown) <= 1 || throw(MacroError("unable to infer ranges for indices $str", where))

    return sizes
end

"""
    star(x,y,...)

Like `*` but intended for multiplying sizes, and understands that `:` is a wildcard.
"""
star(x,y) = *(x,y)
star(::Colon,y) = Colon()
star(x,::Colon) = Colon()
star(x,y,zs...) = star(star(x,y), zs...)

#="""
    needview!([:, 3, :])   # true, need view(A, :,3,:)
    needview!([:, :_, :])  # false, can use rview(A, :,1,:)

Mutates the given vector, replacing symbol `:_` with `1`.
If the vector contains only colons & underscores, then the result is suitable for use with `rview`,
but if not, we need a real view, so it returns `true`.
"""=#
function needview!(getafix::Vector)
    out = false
    for i=1:length(getafix)
        if getafix[i] == :_
            getafix[i] = 1
        elseif getafix[i] isa Int || getafix[i] isa Symbol
            out = true
        elseif getafix[i] isa Colon
        else
            error("this should never happen, getafix[i] = ",getafix[i])
        end
    end
    out
end

"""
    @assert_ cond str

Like `@assert`, but prints both the given string and the condition.
"""
macro assert_(ex, str)
    strex = Main.Base.string(ex)
    msg = str * ": " * strex
    return esc(:($ex || throw(DimensionMismatch($msg))))
end

function unparse(str::String, exs...)
    @capture(exs[1], left_ = right_ ) && return string(str, " ", left, " = ", right, "  ", join(exs[2:end],"  "))
    @capture(exs[1], left_ := right_ ) && return string(str, " ", left, " := ", right, "  ", join(exs[2:end],"  "))
    @capture(exs[1], left_ |= right_ ) && return string(str, " ", left, " |= ", right, "  ", join(exs[2:end],"  "))
    @capture(exs[1], left_ == right_ ) && return string(str, " ", left, " == ", right, "  ", join(exs[2:end],"  "))
    @capture(exs[1], left_ -> right_ ) && return string(str, " ", left, " -> ", right, "  ", join(exs[2:end],"  "))
    @capture(exs[1], left_ => right_ ) && return string(str, " ", left, " => ", right, "  ", join(exs[2:end],"  "))
    return string(exs)
end

function m_error(str::String, where=nothing)
    if where == nothing
        @error str
    else
        @error str  input=where.str _module=where.mod  _line=where.src.line  _file=string(where.src.file)
    end
end

struct MacroError <: Exception
    msg::String
    where::Union{Nothing, NamedTuple}
    MacroError(msg::String, where=nothing) = new(msg, where)
end

function Base.showerror(io::IO, err::MacroError)
    print(io, err.msg)
    if err.where != nothing
        printstyled(io, " \n    ", err.where.str; color = :normal)# :blue)
    end
end


text:
using Test
using MathOptInterface
const MOI = MathOptInterface
const MOIU = MOI.Utilities

function structeq(a::T, b::T) where {T}
    return all(f -> getfield(a, f) == getfield(b, f), fieldnames(T))
end

@testset "parsefunction" begin
    @test structeq(MOIU.parsefunction(:x), MOIU.ParsedSingleVariable(:x))
    @test structeq(
        MOIU.parsefunction(:([x, y, z])),
        MOIU.ParsedVectorOfVariables([:x, :y, :z]),
    )

    @test structeq(
        MOIU.parsefunction(:(x + y + 2.0)),
        MOIU.ParsedScalarAffineFunction(
            MOIU.ParsedScalarAffineTerm.([1.0, 1.0], [:x, :y]),
            2.0,
        ),
    )
    @test structeq(
        MOIU.parsefunction(:(x + -3y + 2.0)),
        MOIU.ParsedScalarAffineFunction(
            MOIU.ParsedScalarAffineTerm.([1.0, -3.0], [:x, :y]),
            2.0,
        ),
    )
    @test structeq(
        MOIU.parsefunction(:(2 * x * y + y + 1.0)),
        MOIU.ParsedScalarQuadraticFunction(
            MOIU.ParsedScalarAffineTerm.([1.0], [:y]),
            MOIU.ParsedScalarQuadraticTerm.([2.0], [:x], [:y]),
            1.0,
        ),
    )

    err = ErrorException("Expected `+`, got `-`.")
    @test_throws err MOIU.parsefunction(:(x - y))

    @test structeq(
        MOIU.parsefunction(:([x, 2x + y + 5.0])),
        MOIU.ParsedVectorAffineFunction(
            MOIU.ParsedVectorAffineTerm.(
                [1, 2, 2],
                MOIU.ParsedScalarAffineTerm.([1.0, 2.0, 1.0], [:x, :x, :y]),
            ),
            [0.0, 5.0],
        ),
    )
    @test structeq(
        MOIU.parsefunction(:([x, 2x + y + 5.0, 1 * x * x])),
        MOIU.ParsedVectorQuadraticFunction(
            MOIU.ParsedVectorAffineTerm.(
                [1, 2, 2],
                MOIU.ParsedScalarAffineTerm.([1.0, 2.0, 1.0], [:x, :x, :y]),
            ),
            MOIU.ParsedVectorQuadraticTerm.(
                [3],
                MOIU.ParsedScalarQuadraticTerm.([2.0], [:x], [:x]),
            ),
            [0.0, 5.0, 0.0],
        ),
    )
end

@testset "separatelabel" begin
    @test MOIU.separatelabel(:(variables:x)) == (:variables, :x)
    @test MOIU.separatelabel(:(variables:x, y)) == (:variables, :((x, y)))
    @test MOIU.separatelabel(:(minobjective:x+y)) == (:minobjective, :(x + y))
    @test MOIU.separatelabel(:(con1:2x <= 1)) == (:con1, :(2x <= 1))
    @test MOIU.separatelabel(:(con1:[x, y] in S)) == (:con1, :([x, y] in S))
end

@testset "loadfromstring" begin
    @testset "one variable" begin
        s = """
        variables: x
        bound: x >= 1.0
        """
        model = MOIU.Model{Float64}()
        x = MOI.add_variable(model)
        MOI.set(model, MOI.VariableName(), x, "x")
        bound = MOI.add_constraint(
            model,
            MOI.SingleVariable(x),
            MOI.GreaterThan(1.0),
        )
        MOI.set(model, MOI.ConstraintName(), bound, "bound")

        model2 = MOIU.Model{Float64}()
        MOIU.loadfromstring!(model2, s)
        MOIU.test_models_equal(model, model2, ["x"], ["bound"])
    end

    @testset "linear constraints" begin
        s = """
        variables: x, y
        linear1: x + y >= 1.0
        linear2: x + y <= 1.0
        linear3: x + y == 1.0
        """
        model = MOIU.Model{Float64}()
        x = MOI.add_variable(model)
        y = MOI.add_variable(model)
        MOI.set(model, MOI.VariableName(), x, "x")
        MOI.set(model, MOI.VariableName(), y, "y")
        linear1 = MOI.add_constraint(
            model,
            MOI.ScalarAffineFunction(
                MOI.ScalarAffineTerm.([1.0, 1.0], [x, y]),
                0.0,
            ),
            MOI.GreaterThan(1.0),
        )
        linear2 = MOI.add_constraint(
            model,
            MOI.ScalarAffineFunction(
                MOI.ScalarAffineTerm.([1.0, 1.0], [x, y]),
                0.0,
            ),
            MOI.LessThan(1.0),
        )
        linear3 = MOI.add_constraint(
            model,
            MOI.ScalarAffineFunction(
                MOI.ScalarAffineTerm.([1.0, 1.0], [x, y]),
                0.0,
            ),
            MOI.EqualTo(1.0),
        )
        MOI.set(model, MOI.ConstraintName(), linear1, "linear1")
        MOI.set(model, MOI.ConstraintName(), linear2, "linear2")
        MOI.set(model, MOI.ConstraintName(), linear3, "linear3")

        model2 = MOIU.Model{Float64}()
        MOIU.loadfromstring!(model2, s)
        MOIU.test_models_equal(
            model,
            model2,
            ["x", "y"],
            ["linear1", "linear2", "linear3"],
        )
    end

    @testset "minimization: linear objective" begin
        s = """
        variables: x, y
        minobjective: x + -2y + 1.0
        """
        model = MOIU.Model{Float64}()
        x = MOI.add_variable(model)
        y = MOI.add_variable(model)
        MOI.set(model, MOI.VariableName(), x, "x")
        MOI.set(model, MOI.VariableName(), y, "y")
        MOI.set(
            model,
            MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
            MOI.ScalarAffineFunction(
                MOI.ScalarAffineTerm.([1.0, -2.0], [x, y]),
                1.0,
            ),
        )
        MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)

        model2 = MOIU.Model{Float64}()
        MOIU.loadfromstring!(model2, s)
        MOIU.test_models_equal(model, model2, ["x", "y"], String[])
    end

    @testset "maximization: linear objective" begin
        s = """
        variables: x, y
        maxobjective: x + -2y + 1.0
        """
        model = MOIU.Model{Float64}()
        x = MOI.add_variable(model)
        y = MOI.add_variable(model)
        MOI.set(model, MOI.VariableName(), x, "x")
        MOI.set(model, MOI.VariableName(), y, "y")
        MOI.set(
            model,
            MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
            MOI.ScalarAffineFunction(
                MOI.ScalarAffineTerm.([1.0, -2.0], [x, y]),
                1.0,
            ),
        )
        MOI.set(model, MOI.ObjectiveSense(), MOI.MAX_SENSE)

        model2 = MOIU.Model{Float64}()
        MOIU.loadfromstring!(model2, s)
        MOIU.test_models_equal(model, model2, ["x", "y"], String[])
    end

    @testset "SOC constraints" begin
        s = """
        variables: x, y, z
        varsoc: [x,y,z] in SecondOrderCone(3)
        affsoc: [2x,y+1,-1*z] in SecondOrderCone(3)
        affsoc2: [1.0,2.0,3.0] in SecondOrderCone(3)
        """
        model = MOIU.Model{Float64}()
        x = MOI.add_variable(model)
        y = MOI.add_variable(model)
        z = MOI.add_variable(model)
        MOI.set(model, MOI.VariableName(), x, "x")
        MOI.set(model, MOI.VariableName(), y, "y")
        MOI.set(model, MOI.VariableName(), z, "z")
        varsoc = MOI.add_constraint(
            model,
            MOI.VectorOfVariables([x, y, z]),
            MOI.SecondOrderCone(3),
        )
        affsoc = MOI.add_constraint(
            model,
            MOI.VectorAffineFunction(
                MOI.VectorAffineTerm.(
                    [1, 2, 3],
                    MOI.ScalarAffineTerm.([2.0, 1.0, -1.0], [x, y, z]),
                ),
                [0.0, 1.0, 0.0],
            ),
            MOI.SecondOrderCone(3),
        )
        affsoc2 = MOI.add_constraint(
            model,
            MOI.VectorAffineFunction(
                MOI.VectorAffineTerm{Float64}[],
                [1.0, 2.0, 3.0],
            ),
            MOI.SecondOrderCone(3),
        )
        MOI.set(model, MOI.ConstraintName(), varsoc, "varsoc")
        MOI.set(model, MOI.ConstraintName(), affsoc, "affsoc")
        MOI.set(model, MOI.ConstraintName(), affsoc2, "affsoc2")

        model2 = MOIU.Model{Float64}()
        MOIU.loadfromstring!(model2, s)
        MOIU.test_models_equal(
            model,
            model2,
            ["x", "y", "z"],
            ["varsoc", "affsoc", "affsoc2"],
        )
    end

    @testset "Invalid variable name" begin
        s = """
        variables: x
        bound: y >= 1.0
        """
        model = MOIU.Model{Float64}()
        err = ErrorException("Invalid variable name y.")
        @test_throws err MOIU.loadfromstring!(model, s)
    end
end


text:
module Mixtape

using MacroTools: isexpr, @capture, rmlines, unblock, postwalk
using CodeInfoTools
using Core.Compiler
using Core.Compiler: MethodInstance, NativeInterpreter, CodeInfo, CodeInstance, WorldView,
                     OptimizationState, Const, widenconst, MethodResultPure, CallMeta
import Core.Compiler.abstract_call_known
using LLVM
using LLVM.Interop
using LLVM_full_jll

# Cache
using Core.Compiler: WorldView

# AbstractInterpreter
import Core.Compiler: AbstractInterpreter, InferenceResult, InferenceParams, InferenceState,
                      OptimizationParams
import Core.Compiler: get_world_counter, get_inference_cache, code_cache, lock_mi_inference,
                      unlock_mi_inference
import Core.Compiler: retrieve_code_info, validate_code_in_debug_mode

# Optimizer
import Core.Compiler: optimize
import Core.Compiler: CodeInfo, convert_to_ircode, copy_exprargs, slot2reg, compact!,
                      coverage_enabled, adce_pass!
import Core.Compiler: ssa_inlining_pass!, getfield_elim_pass!, type_lift_pass!, verify_ir,
                      verify_linetable

# JIT
using GPUCompiler: GPUCompiler, CompilerJob
import GPUCompiler: FunctionSpec, AbstractCompilerTarget, AbstractCompilerParams

import Base: show, push!

resolve(x) = x
resolve(gr::GlobalRef) = getproperty(gr.mod, gr.name)
resolve(c::Core.Const) = c.val

# Exports.
export CompilationContext, allow, transform, optimize!, show_after_inference,
       show_after_optimization, debug, @ctx, @intrinsic

export widen_invokes

include("cache.jl")
include("world.jl")
include("utility_transforms.jl")
include("interpreter.jl")
include("codegen.jl")
include("optimize.jl")
include("llvmopt.jl")
include("interface.jl")
include("init.jl")

end # module


text:
# Copyright 2022 RelationalAI, Inc.
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License

# Execute the given query string.
# Example:
# $ julia --proj=. examples/exec_async.jl "nhd-test-1" "nhd-s" "2+2"
# Transaction is done...
# JSON3.Object{Vector{UInt8}, Vector{UInt64}} with 2 entries:
#   :id    => "261f5d59-f1b4-f778-4c13-a7993871c972"
#   :state => "CREATED"

import RAI
using RAI: Context, HTTPError, exec_async, load_config, show_result, get_transaction

include("parseargs.jl")

function run(database, engine, source; profile)
    conf = load_config(; profile = profile)
    ctx = Context(conf)
    txn = exec_async(ctx, database, engine, source)
    println("Transaction is created...")
    display(txn)
    println()
end

function main()
    args = parseargs(
        "database", Dict(:help => "database name", :required => true),
        "engine", Dict(:help => "engine name", :required => true),
        "command", Dict(:help => "rel source string"),
        ["--file", "-f"], Dict(:help => "rel source file"),
        "--readonly", Dict(:help => "readonly query (default: false)", :action => "store_true"),
        "--profile", Dict(:help => "config profile (default: default)"))
    try
        source = nothing
        if !isnothing(args.command)
            source = args.command
        elseif !isnothing(args.file)
            source = open(args.file, "r")
        end
        isnothing(source) && return  # nothing to execute
        run(args.database, args.engine, source; profile = args.profile)
    catch e
        e isa HTTPError ? show(e) : rethrow()
    end
end

main()

max_stars_repo_name: RelationalAI/rai-sdk-julia

max_issues_repo_name: RelationalAI/rai-sdk-julia

max_forks_repo_name: RelationalAI/rai-sdk-julia


text:
using AbstractStorage, AzSessions, AzStorage, Dates, JSON, Random, Serialization, Test

credentials = JSON.parse(ENV["AZURE_CREDENTIALS"])
AzSessions.write_manifest(;client_id=credentials["clientId"], client_secret=credentials["clientSecret"], tenant=credentials["tenantId"])

session = AzSession(;protocal=AzClientCredentials, client_id=credentials["clientId"], client_secret=credentials["clientSecret"], resource="https://storage.azure.com/")

storageaccount = ENV["STORAGE_ACCOUNT"]
@info "storageaccount=$storageaccount"

for container in containers(;storageaccount=storageaccount,session=session)
    rm(AzContainer(container;storageaccount=storageaccount,session=session))
end
@info "sleeping for 60 seconds to ensure Azure clean-up from any previous run"
sleep(60)

@testset "Error codes" begin
    @test unsafe_load(cglobal((:N_HTTP_RETRY_CODES, AzStorage.libAzStorage), Cint)) == 2
    x = unsafe_load(cglobal((:HTTP_RETRY_CODES, AzStorage.libAzStorage), Ptr{Clong}))
    y = unsafe_wrap(Array, x, (2,); own=false)
    @test y == [500,503]

    @test unsafe_load(cglobal((:N_CURL_RETRY_CODES, AzStorage.libAzStorage), Cint)) == 4
    x = unsafe_load(cglobal((:CURL_RETRY_CODES, AzStorage.libAzStorage), Ptr{Clong}))
    y = unsafe_wrap(Array, x, (4,); own=false)
    @test y == [7,28,55,56]
end

@testset "Containers, equivalent" begin
    x = AzContainer("foo/bar"; storageaccount="baz", nthreads=10)
    y = AzContainer("foo/bar"; storageaccount="baz", nthreads=11)
    z = AzContainer("foo/bar"; storageaccount="fiz", nthreads=11)
    w = AzContainer("foo/fiz"; storageaccount="baz", nthreads=11)
    @test x == y
    @test x != z
    @test x != w
end

@testset "Containers, list" begin
    sleep(1)
    r = lowercase(randstring(MersenneTwister(millisecond(now())+0)))
    c = AzContainer("foo-$r-a", storageaccount=storageaccount, session=session)
    mkpath(c)

    write(c, "bar", "one")
    write(c, "baz", "two")

    l = readdir(c)

    @test "bar" ∈ l
    @test "baz" ∈ l

    containers(storageaccount=c.storageaccount, session=c.session)
    @test isdir(c)
    rm(c)
    @test !isdir(c)
end

@testset "Containers, prefix, list" begin
    sleep(1)
    r = lowercase(randstring(MersenneTwister(millisecond(now())+1)))
    c = AzContainer("foo-$r-b", prefix="prefix", storageaccount=storageaccount, session=session)
    mkpath(c)

    write(c, "bar", "one")
    write(c, "baz", "two")

    l = readdir(c)

    @test "bar" ∈ l
    @test "baz" ∈ l

    l = readdir(c, filterlist=false)
    @test "prefix/bar" ∈ l
    @test "prefix/baz" ∈ l

    containers(storageaccount=c.storageaccount, session=c.session)
    @test isdir(c)
    rm(c)
    @test !isdir(c)
end

@testset "Containers, prefix, list, alt construction" begin
    sleep(1)
    r = lowercase(randstring(MersenneTwister(millisecond(now())+2)))
    c = AzContainer("foo-$r-c/prefix", storageaccount=storageaccount, session=session)
    mkpath(c)

    write(c, "bar", "one")
    write(c, "baz", "two")

    l = readdir(c)

    @test "bar" ∈ l
    @test "baz" ∈ l

    l = readdir(c, filterlist=false)
    @test "prefix/bar" ∈ l
    @test "prefix/baz" ∈ l

    containers(storageaccount=c.storageaccount, session=c.session)
    @test isdir(c)
    rm(c)
    @test !isdir(c)
end

@testset "Containers, dirname" begin
    sleep(1)
    r = lowercase(randstring(MersenneTwister(millisecond(now())+3)))
    c = AzContainer("foo-$r-d", storageaccount=storageaccount, session=session)
    @test dirname(c) == "foo-$r-d"
end

@testset "Containers, prefix, dirname" begin
    sleep(1)
    r = lowercase(randstring(MersenneTwister(millisecond(now())+4)))
    c = AzContainer("foo-$r-e", prefix="prefix", storageaccount=storageaccount, session=session)
    @test dirname(c) == "foo-$r-e/prefix"
end

@testset "Containers, size, prefix=$prefix" for prefix in ("", "prefix")
    sleep(1)
    suffix = prefix == "" ? "-foo" : "-bar"
    r = lowercase(randstring(MersenneTwister(millisecond(now())+5)))
    c = AzContainer("foo-$r-f$suffix", prefix=prefix, storageaccount=storageaccount, session=session)
    mkpath(c)

    write(c, "bar", rand(UInt8,10))
    write(c, "baz", rand(UInt8,11))
    @test filesize(c, "bar") == 10
    @test filesize(c, "baz") == 11

    rm(c)
end

@testset "Containers, bytes, nthreads=$nthreads, prefix=$prefix" for nthreads in (1, 2), prefix in ("","prefix")
    sleep(1)
    suffix = prefix == "" ? "-foo" : "-bar"
    r = lowercase(randstring(MersenneTwister(millisecond(now())+6)))
    c = AzContainer("foo-$r-g$suffix", prefix=prefix, storageaccount=storageaccount, session=session, nthreads=nthreads)
    mkpath(c)

    N = 10
    if nthreads == 2
        N = round(Int, AzStorage._MINBYTES_PER_BLOCK * nthreads * 5 * (1+rand()) / 8)
        nblks = AzStorage.nblocks(nthreads, N*8)
        bytes_per_block = div(8*N, nblks) + (rem(8*N, nblks) == 0 ? 0 : 1)
        @info "bytes_per_block=$bytes_per_block, AzStorage._MAXBYTES_PER_BLOCK=$(AzStorage._MAXBYTES_PER_BLOCK)"
    end

    x = rand(N)
    write(c, "bar", x)
    y = read!(c, "bar", Vector{Float64}(undef, N))
    @test x ≈ y
    rm(c)
end

@testset "Containers, bytes, nested folder, prefix=$prefix" for prefix in ("","prefix")
    sleep(1)
    r = lowercase(randstring(MersenneTwister(millisecond(now())+7)))
    suffix = prefix == "" ? "-foo" : "-bar"
    c = AzContainer("foo-$r-e$suffix", prefix=prefix, storageaccount=storageaccount, session=session)
    mkpath(c)

    x = rand(10)
    write(c, "bar/baz", x)
    y = read!(c, "bar/baz", Vector{Float64}(undef, 10))
    @test x ≈ y
    rm(c)
end

@testset "Containers, string, prefix=$prefix" for prefix in ("", "prefix")
    sleep(1)
    r = lowercase(randstring(MersenneTwister(millisecond(now())+8)))
    suffix = prefix == "" ? "-foo" : "-bar"
    c = AzContainer("foo-$r-f$suffix", prefix=prefix, storageaccount=storageaccount, session=session)
    mkpath(c)

    write(c, "bar", "hello world\n")
    @test read(c, "bar", String) == "hello world\n"
    rm(c)
end

@testset "Containers, rm, prefix=$prefix" for prefix in ("prefix","")
    sleep(1)
    r = lowercase(randstring(MersenneTwister(millisecond(now())+9)))
    suffix = prefix == "" ? "-foo" : "-bar"
    c = AzContainer("foo-$r-g$suffix", prefix=prefix, storageaccount=storageaccount, session=session)
    mkpath(c)
    write(c, "bar", "bar\n")
    write(c, "foo", "foo\n")
    rm(c, "bar")
    @test readdir(c) == ["foo"]
    rm(c, "baz")
    @test readdir(c) == ["foo"]
    rm(c)
end

@testset "Containers, isfile, prefix=$prefix" for prefix in ("","prefix")
    sleep(1)
    r = lowercase(randstring(MersenneTwister(millisecond(now())+10)))
    suffix = prefix == "" ? "-foo" : "-bar"
    c = AzContainer("foo-$r-h$suffix", prefix=prefix, storageaccount=storageaccount, session=session)
    mkpath(c)
    write(c, "bar", "bar\n")
    @test isfile(c, "bar")
    write(c, "bar/baz", "bar\n")
    @test isfile(c, "bar/baz")
    @test !isfile(c, "notanobject")
    rm(c)
end

@testset "Containers, cp, prefix=$prefix" for prefix in ("", "prefix")
    sleep(1)
    r = lowercase(randstring(MersenneTwister(millisecond(now())+11)))
    suffix = prefix == "" ? "-foo" : "-bar"
    src = AzContainer("foo-$r-i$suffix", storageaccount=storageaccount, session=session)
    r = lowercase(randstring(MersenneTwister(millisecond(now())+12)))
    dst = AzContainer("foo-$r-j$suffix", storageaccount=storageaccount, session=session)

    mkpath(src)
    write(src, "bar", "one")
    write(src, "baz", "two")

    cp(src, dst)

    @test read(dst, "bar", String) == "one"
    @test read(dst, "baz", String) == "two"

    rm(src)
    rm(dst)
end

@testset "Containers, json" begin
    sleep(1)
    r = lowercase(randstring(MersenneTwister(millisecond(now())+13)))
    c = AzContainer("foo-$r-k", storageaccount=storageaccount, session=session, nthreads=2, nretry=10)
    mkpath(c)
    _c = Container(AzContainer, JSON.parse(json(c)), c.session)
    @test _c.storageaccount == storageaccount
    @test _c.containername == "foo-$r-k"
    @test _c.nretry == 10
    @test _c.nthreads == 2
    rm(c)
end

@testset "Containers, serialization" begin
    sleep(1)
    r = lowercase(randstring(MersenneTwister(millisecond(now())+14)))
    c = AzContainer("foo-$r-k", storageaccount=storageaccount, session=session, nthreads=2, nretry=10)
    mkpath(c)
    x = (a=rand(10), b=rand(10))
    serialize(c, "bar", x)
    _x = deserialize(c, "bar")
    @test x.a ≈ _x.a
    @test x.b ≈ _x.b
    rm(c)
end

@testset "Object, bytes" begin
    sleep(1)
    r = lowercase(randstring(MersenneTwister(millisecond(now())+15)))
    c = AzContainer("foo-$r-k", storageaccount=storageaccount, session=session, nthreads=2, nretry=10)
    io = open(c, "bar")
    x = rand(10)
    write(io, x)
    _x = read!(io, zeros(10))
    rm(c)
    @test x ≈ _x
end

@testset "Object, bytes" begin
    sleep(1)
    r = lowercase(randstring(MersenneTwister(millisecond(now())+16)))
    c = AzContainer("foo-$r-k", storageaccount=storageaccount, session=session, nthreads=2, nretry=10)
    io = open(c, "bar")
    write(io, "hello")
    x = read(io, String)
    rm(c)
    @test x == "hello"
end

@testset "Object, isfile" begin
    sleep(1)
    r = lowercase(randstring(MersenneTwister(millisecond(now())+17)))
    c = AzContainer("foo-$r-k", storageaccount=storageaccount, session=session, nthreads=2, nretry=10)
    io = open(c, "bar")
    write(io, "hello")
    @test isfile(io)
    rm(c)
end

@testset "Object, rm" begin
    sleep(1)
    r = lowercase(randstring(MersenneTwister(millisecond(now())+18)))
    c = AzContainer("foo-$r-k", storageaccount=storageaccount, session=session, nthreads=2, nretry=10)
    io = open(c, "bar")
    write(io, "hello")
    @test isfile(io)
    rm(io)
    @test !isfile(io)
    rm(c)
end

@testset "Object, joinpath" begin
    sleep(1)
    r = lowercase(randstring(MersenneTwister(millisecond(now())+19)))
    c = AzContainer("foo-$r-k", storageaccount=storageaccount, session=session, nthreads=2, nretry=10)
    io = joinpath(c, "bar", "baz")
    write(io, "hello")
    @test read(io, String) == "hello"
    rm(c)
end

@testset "Object, serialization" begin
    sleep(1)
    r = lowercase(randstring(MersenneTwister(millisecond(now())+20)))
    c = AzContainer("foo-$r-k", storageaccount=storageaccount, session=session, nthreads=2, nretry=10)
    io = open(c, "bar")
    x = (a=rand(10), b=rand(10))
    serialize(io, x)
    _x = deserialize(io)
    @test x.a ≈ _x.a
    @test x.b ≈ _x.b
end

@testset "Object, number of blocks calculation" begin
    nthreads = 16
    @test AzStorage.nblocks(nthreads, nthreads*AzStorage._MAXBYTES_PER_BLOCK) == nthreads
    @test AzStorage.nblocks(nthreads, nthreads*AzStorage._MINBYTES_PER_BLOCK) == nthreads
    @test AzStorage.nblocks(nthreads, div(nthreads*AzStorage._MINBYTES_PER_BLOCK,2)) == div(nthreads, 2)
    @test AzStorage.nblocks(nthreads, 2*nthreads*AzStorage._MAXBYTES_PER_BLOCK) == 2*nthreads
    @test_throws ErrorException AzStorage.nblocks(nthreads, 2*AzStorage._MAXBYTES_PER_BLOCK*AzStorage._MAXBLOCKS_PER_BLOB)
end

# this failed because of a bug in the block-list when using exactly 10 blocks
# Anusha found the failing example.
@testset "Anusha's example" begin
    sleep(1)
    r = lowercase(randstring(MersenneTwister(millisecond(now())+21)))
    c = AzContainer("foo-$r-l", storageaccount=storageaccount, session=session, nthreads=2, nretry=10)
    mkpath(c)
    x = rand(2801,13821)
    write(c, "bar", x)
    _x = read!(c, "bar", Array{Float64,2}(undef,2801,13821))
    @test x ≈ _x
    rm(c)
end

@testset "writedlm and readdlm" begin
    sleep(1)
    a = rand(1000,1000)
    r = lowercase(randstring(MersenneTwister(millisecond(now())+22)))
    c = AzContainer("foo-$r-m", storageaccount=storageaccount, session=session, nthreads=2, nretry=10)
    io = open(c, "bar")
    writedlm(io,a)
    _a = readdlm(io)
    @test _a ≈ a
    rm(c)
end

@testset "Containers, bytes, SubArray" begin
    sleep(1)
    a = rand(10,20,3)
    r = lowercase(randstring(MersenneTwister(millisecond(now())+23)))
    c = AzContainer("foo-$r-n", storageaccount=storageaccount, session=session, nthreads=2, nretry=10)
    mkpath(c)
    for i = 1:3
        write(c, "bar$i", @view a[:,:,i])
    end
    _a = zeros(10,20,3)
    for i = 1:3
        read!(c, "bar$i", @view _a[:,:,i])
    end
    @test a ≈ _a
    rm(c)
end

@testset "Containers, bytes, non-contiguous throws" begin
    sleep(1)
    a = rand(10,20)
    r = lowercase(randstring(MersenneTwister(millisecond(now())+24)))
    c = AzContainer("foo-$r-o", storageaccount=storageaccount, session=session, nthreads=2, nretry=10)
    mkpath(c)
    @test_throws ErrorException write(c, "bar", @view a[1:2:end,1:2:end])
    rm(c)
end

@testset "Container, minimal dictionary" begin
    sleep(1)
    r = lowercase(randstring(MersenneTwister(millisecond(now())+25)))
    c = AzContainer("foo-$r-o", storageaccount=storageaccount, session=session, nthreads=2, nretry=10)
    _c = minimaldict(c)
    @test _c["storageaccount"] == storageaccount
    @test _c["containername"] == "foo-$r-o"
    @test _c["prefix"] == ""
    @test length(_c) == 3
end


text:
# H = JΣⱼ (½ S⁺ⱼS⁻ⱼ₊₁ + ½ S⁻ⱼS⁺ⱼ₊₁) + J₂Σⱼ (½ S⁺ⱼZⱼ₊₁S⁻ⱼ₊₂ + ½ S⁻ⱼZⱼ₊₁S⁺ⱼ₊₂) - h Σⱼ Sⱼᶻ
function ITensors.OpSum(::Model"xx_extended", n1, n2; J=1.0, J₂=1.0, h=0.0)
  opsum = OpSum()
  if J != 0.0
    opsum += 0.25 * J, "S+", n1, "S-", n2
    opsum += 0.25 * J, "S-", n1, "S+", n2
    opsum += 0.25 * J, "S+", n2, "S-", n2 + 1
    opsum += 0.25 * J, "S-", n2, "S+", n2 + 1
  end
  if J₂ != 0
    opsum += 0.5 * J₂, "S+", n1, "Sz", n2, "S-", n2 + 1
    opsum += 0.5 * J₂, "S-", n1, "Sz", n2, "S+", n2 + 1
  end
  if h != 0
    opsum += -h / 3, "Sz", n1
    opsum += -h / 3, "Sz", n2
    opsum += -h / 3, "Sz", n2 + 1
  end
  return opsum
end

# H = JΣⱼ (½ S⁺ⱼS⁻ⱼ₊₁ + ½ S⁻ⱼS⁺ⱼ₊₁) + J₂Σⱼ (½ S⁺ⱼS⁻ⱼ₊₂ + ½ S⁻ⱼS⁺ⱼ₊₂) - h Σⱼ Sⱼᶻ
function ITensors.MPO(::Model"xx_extended", s; J=1.0, J₂=1.0, h=0.0)
  N = length(s)
  os = OpSum()
  if h != 0
    for j in 1:N
      os += -h, "Sz", j
    end
  end
  if J != 0
    for j in 1:(N - 1)
      os += 0.5 * J, "S+", j, "S-", j + 1
      os += 0.5 * J, "S-", j, "S+", j + 1
    end
  end
  if J₂ != 0
    for j in 1:(N - 2)
      os += 0.5 * J₂, "S+", j, "Sz", j + 1, "S-", j + 2
      os += 0.5 * J₂, "S-", j, "Sz", j + 1, "S+", j + 2
    end
  end
  return splitblocks(linkinds, MPO(os, s))
end

nrange(::Model"xx_extended") = 3

function ITensorInfiniteMPS.reference(
  ::Model"xx_extended", ::Observable"energy"; N=1000, J=1.0, J₂=1.0, h=0.0, filling=0.5
)
  mat = zeros(N, N)
  for j in 1:(N - 1)
    mat[j, j + 1] = -J
    mat[j + 1, j] = -J
  end
  for j in 1:(N - 2)
    mat[j, j + 2] = -J₂
    mat[j + 2, j] = -J₂
  end
  temp = sort(eigvals(mat))
  return sum(temp[1:round(Int64, filling * N)]) / N / 2 - h * (1 / 2 + filling)
end

function ITensorInfiniteMPS.reference(
  ::Model"xx_extended", ::Observable"energy"; N=1000, J=1.0, J₂=1.0, h=0.0, filling=0.5
)
  λ(x) = -J / 2 * cos(x * 2 * pi / N) - J₂ / 2 * cos(x * 4 * pi / N)
  temp = sort(λ.(0:(N - 1)))
  return sum(temp[1:round(Int64, filling * N)]) / N - h * (2 * filling - 1)
end


text:
################################################################################
#                          TEST: DCM <=> Quaternions
################################################################################

for k = 1:samples
    # Conversion test
    # ===============

    # Create a random quaternion.
    q1 = Quaternion(randn(), randn(), randn(), rand())
    q1 = q1/norm(q1)

    # Make sure that the real part is positive.
    (q1.q0 < 0) && (q1 *= -1)

    # Create a random DCM using Euler Angles.
    dcm2 = angle2dcm(randn(), randn(), randn(), rot_seq_array[rand(1:9)])

    # Convert the quaternion in DCM and the DCM in quaternion.
    q2   = dcm2quat(dcm2)
    dcm1 = quat2dcm(q1)

    # Do the inverse conversion and compare.
    error_dcm = dcm2 - quat2dcm(q2)
    @test norm(error_dcm) < 1e-10

    error_q = q1 - dcm2quat(dcm1)
    @test norm(error_q) < 1e-10

    # Multiplication test
    # ===================

    # Sample a 3x1 vector.
    v = randn(3)

    # Compute the rotated vector using the quaternions.
    v_rot_q = vect(inv(q1*q2)*v*(q1*q2))

    # Compute the rotated vector using the DCMs.
    v_rot_dcm = dcm2*dcm1*v

    # Compare.
    @test norm(v_rot_q-v_rot_dcm) < 1e-10
end


text:
function integralweight!(Wim::Matrix{Float64},
    X::Matrix{Float64}, Y::AbstractArray{Bool, 1},
    groupindex::IntegerVector, gammaM::Vector{Float64},
    wi::Vector{Float64}, ghw::Vector{Float64},
    llN::Vector{Float64}, llN2::Vector{Float64},
    xb::Vector{Float64}, N::Int, J::Int,
    n::Int, C::Int, ngh::Int)
    #A_mul_B!(xb, X, betas)
    ll = 0.
    for jcom in 1:C
        for ix in 1:ngh
            wtmp = log(ghw[ix])+log(wi[jcom])
            ixM = ix+ngh*(jcom-1)
            fill!(llN, gammaM[ixM])
            Yeppp.add!(llN, llN, xb)
            negateiftrue!(llN, Y)
            # Yeppp.exp!(llN, llN)
            # log1p!(llN)
            log1pexp!(llN, llN, llN2, N)

            for i in 1:n
                Wim[i, ixM] = wtmp
            end
            for i in 1:N
                @inbounds Wim[groupindex[i], ixM] -= llN[i]
            end
        end
    end
    for i in 1:n
        ll += logsumexp(Wim[i,:])
    end
    for i in 1:n
        u = maximum(Wim[i, :])
        for jcol in 1:C*ngh
            @inbounds Wim[i, jcol] = Wim[i, jcol] - u
        end
    end
    Yeppp.exp!(Wim, Wim)
    for i in 1:n
        u = sum(Wim[i, :])
        for jcol in 1:C*ngh
            @inbounds Wim[i, jcol] = Wim[i, jcol] / u
        end
    end
    return ll - n*logπ/2
end

function updateθ!(wi::Vector{Float64}, mu::Vector{Float64},
    sigmas::Vector{Float64}, X::Matrix{Float64},
    Y::AbstractArray{Bool, 1}, groupindex::IntegerVector,
    gammaM::Vector{Float64}, Wim::Matrix{Float64},
    Wm::Matrix{Float64}, sn::Vector{Float64},
    an::Real, N::Int, J::Int, n::Int, C::Int, ngh::Int)

    # A_mul_B!(xb, X, betas)
    mean!(Wm, Wim)
    divide!(Wm, Wm, sum(Wm))
    for kcom in 1:C
        ind = (1+ngh*(kcom-1)):ngh*kcom
        wi[kcom] = sum(Wm[ind])
        mu[kcom] = wsum(gammaM[ind], Wm[ind]) / wi[kcom]
        sigmas[kcom] = sqrt((wsum((gammaM[ind] .- mu[kcom]).^2, Wm[ind]) + 2 * an * sn[kcom]^2/n) / (wi[kcom] + 2 * an/n))
    end

end
function updateβ!(β::Vector{Float64}, X::Matrix{Float64},
    Y::AbstractArray{Bool, 1}, groupindex::IntegerVector,
    minStepFac::Real, betadevtol::Real,
    XWX::Matrix{Float64}, XWY::Vector{Float64},
    Xscratch::Matrix{Float64},
    gammaM::Vector{Float64}, Wim::Matrix{Float64},
    lln::Vector{Float64}, llN::Vector{Float64},
    llN2::Vector{Float64}, llN3::Vector{Float64}, xb::Vector{Float64},
    N::Int, J::Int, n::Int, C::Int, ngh::Int, Qmaxiteration::Int)
    M = C*ngh
    dev0 = negdeviance(β, X, Y, groupindex,
    gammaM, Wim, lln, llN, llN2, xb, N, J, n, M)

    for iterbeta in 1:Qmaxiteration
        deltabeta!(XWY, XWX, X, Y, groupindex, β, Xscratch, gammaM, Wim,
        llN, llN2, llN3, xb, N, J, n, M)
        if sum(abs2, XWY) < J*1e-16
            break
        end
        f = 1.
        dev = negdeviance(β .+ f .* XWY, X, Y, groupindex,
        gammaM, Wim, lln, llN, llN2, xb, N, J, n, M)
        while dev - dev0 > abs(dev)*betadevtol*0.001 #dev < dev0
            f /= 2.0
            f > minStepFac || error("step-halving failed at beta = $(β), deltabeta=$(XWY), dev=$(dev), dev0=$dev0, f=$f")
            dev = negdeviance(β .+ f .* XWY, X, Y, groupindex,
            gammaM, Wim, lln, llN, llN2, xb, N, J, n, M)
        end
        for j in 1:J
            β[j] += f * XWY[j]
        end

        if (dev0 - dev)/dev < 0.001*betadevtol || dev == 0.0
            break
        end
        dev0 = dev
    end

end
function deltabeta!(XWY::Vector{Float64},
    XWX::Matrix{Float64}, X::Matrix{Float64},
    Y::AbstractArray{Bool, 1}, groupindex::IntegerVector,
    β::Vector{Float64}, Xscratch::Matrix{Float64},
    gammaM::Vector{Float64}, Wim::Matrix{Float64},
    llN::Vector{Float64},llN2::Vector{Float64},
    llN3::Vector{Float64}, xb::Vector{Float64},
    N::Int, J::Int, n::Int, M::Int)
    LinearAlgebra.mul!(xb, X, β)
    fill!(XWX, 0.)
    fill!(XWY, 0.)
    fill!(llN2, 0.0)
    fill!(llN3, 0.0)

    for jcol in 1:M
        fill!(llN, gammaM[jcol])
        Yeppp.add!(llN, llN, xb)
        logistic!(llN, llN, N)

        @inbounds for i in 1:N
            Wim_im = Wim[groupindex[i],jcol]
            llN2[i] += Wim_im * llN[i]*(1.0 -llN[i]) #working weight
            llN3[i] += -llN[i]*Wim_im #working response without denominator
            if Y[i]
                llN3[i] += Wim_im
            end
        end
    end
    #relocate!(llN3, Wim[:, jcol], groupindex, N)
    copy!(Xscratch, X)
    LinearAlgebra.BLAS.gemv!('T', 1.0, Xscratch, llN3, 1.0, XWY)
    LinearAlgebra.mul!(Xscratch, Diagonal(llN2), Xscratch)
    LinearAlgebra.BLAS.gemm!('T', 'N', 1.0, Xscratch, X, 1.0, XWX)

    ldiv!(cholesky!(Hermitian(XWX, :U)), XWY)
end
#actual deviance, no longer negative
function negdeviance(beta2::Vector{Float64}, X::Matrix{Float64},
    Y::AbstractArray{Bool, 1}, groupindex::IntegerVector,
    gammaM::Vector{Float64}, Wim::Matrix{Float64},
    lln::Vector{Float64}, llN::Vector{Float64},
    llN2::Vector{Float64}, xb::Vector{Float64},
    N::Int, J::Int, n::Int, M::Int)
    dev = 0.
    LinearAlgebra.mul!(xb, X, beta2)
    for jcol in 1:M
        fill!(llN, gammaM[jcol])
        Yeppp.add!(llN, llN, xb)
        negateiftrue!(llN, Y)
        log1pexp!(llN, llN, llN2, N)
        fill!(lln, 0.0)
        for i in 1:N
            @inbounds lln[groupindex[i]] += llN[i]
        end
        dev += wsum(lln, Wim[:, jcol])
    end
    dev
end

function latentgmm(X::Matrix{Float64},
    Y::AbstractArray{Bool, 1}, groupindex::IntegerVector,
    ncomponent::Int, β_init::Vector{Float64},
    wi_init::Vector{Float64}, mu_init::Vector{Float64},
    sigmas_init::Vector{Float64};
    maxiteration::Int=100, tol::Real=.005,
    ngh::Int=100, ghx::Vector=zeros(ngh), ghw::Vector=zeros(ngh),
    sn::Vector{Float64}=sigmas_init, an::Float64=1.0/maximum(groupindex),
    debuginfo::Bool=false, Qmaxiteration::Int=30,
    whichtosplit::Int=1, tau::Real=.5, taufixed::Bool=false,
    mu_lb::Vector=fill(-Inf, ncomponent),
    mu_ub::Vector=fill(Inf, ncomponent),
    Wim::Matrix{Float64}=zeros(maximum(groupindex), ncomponent*ngh),
    Wm::Matrix{Float64}=zeros(1, ncomponent*ngh),
    lln::Vector{Float64}=zeros(maximum(groupindex)),
    llN::Vector{Float64}=zeros(length(Y)),
    llN2::Vector{Float64}=zeros(length(Y)),
    llN3::Vector{Float64}=zeros(length(Y)),
    Xscratch::Matrix{Float64}=copy(X),
    xb::Vector{Float64}=zeros(length(Y)),
    gammaM::Vector{Float64}=zeros( ncomponent*ngh),
    dotest::Bool=false, epsilon::Real=1e-6,
    updatebeta::Bool=true, pl::Bool=false, ptau::Bool=false, bn::Real=3.0 / maximum(groupindex))

    # initialize theta
    length(wi_init) == length(mu_init) == length(sigmas_init) == ncomponent || error("The length of initial values should be $ncomponent")
    N,J=size(X)
    length(β_init) == J || error("Initial values of fixed efffect coefficients should have same dimension as X")
    n = maximum(groupindex)
    M = ncomponent*ngh

    wi = copy(wi_init)
    mu = copy(mu_init)
    sigmas = copy(sigmas_init)
    β = copy(β_init)

    wi_old = ones(ncomponent)./ncomponent
    mu_old = zeros(ncomponent)
    sigmas_old = ones(ncomponent)
    beta_old = randn(J)
    if !taufixed || (ghx[1] == 0.0)
        ghx, ghw = gausshermite(ngh)
    end
    XWX = zeros(J, J)
    XWY = zeros(J)

    ll0=-Inf
    ll = 0.
    alreadystable = maxiteration <= 3
    for iter_em in 1:maxiteration
        for ix in 1:ngh, jcom in 1:ncomponent
            ixM = ix+ngh*(jcom-1)
            gammaM[ixM] = ghx[ix]*sigmas[jcom]*sqrt(2)+mu[jcom]
        end

        copy!(wi_old, wi)
        copy!(mu_old, mu)
        copy!(sigmas_old, sigmas)

        LinearAlgebra.mul!(xb, X, β)
        ll=integralweight!(Wim, X, Y, groupindex, gammaM, wi, ghw, llN, llN2, xb, N, J, n, ncomponent, ngh)
        lldiff = ll - ll0
        ll0 = ll

        if dotest
            if (lldiff < epsilon) && (iter_em > 3)
                break
            end
        end
        if debuginfo
            println("At $(iter_em)th iteration:")
        end
        if updatebeta && (mod1(iter_em, 3) == 1 || alreadystable)
            copy!(beta_old, β)
            updateβ!(β, X, Y, groupindex, .001, .001,
            XWX, XWY, Xscratch, gammaM, Wim, lln, llN, llN2, llN3,
            xb, N, J, n, ncomponent, ngh, Qmaxiteration)
            if debuginfo
                println("beta=", β)
            end
        end
        updateθ!(wi, mu, sigmas, X, Y, groupindex,
        gammaM, Wim, Wm, sn, an, N, J, n, ncomponent, ngh)

        for kcom in 1:ncomponent
            if wi[kcom] < bn
                wi[kcom] = bn
            end
        end
        tmpsum = sum(wi)
        for kcom in 1:ncomponent
            wi[kcom] = wi[kcom] / tmpsum
        end

        if taufixed
            Yeppp.max!(mu, mu, mu_lb)
            Yeppp.min!(mu, mu, mu_ub)
            wi_tmp = wi[whichtosplit]+wi[whichtosplit+1]
            wi[whichtosplit] = wi_tmp*tau
            wi[whichtosplit+1] = wi_tmp*(1-tau)
        end
        if any(wi .< 1e-8)
            @warn("Some elements of $wi are too small. Consider another starting value or reduce the number of components. Give up.")
            break
        end
        if debuginfo
            println("wi=$wi")
            println("mu=$mu")
            println("sigma=$sigmas")
            println("ll=",ll)
        end

        if !dotest
            if stopRule(vcat(β, wi, mu, sigmas), vcat(beta_old, wi_old, mu_old, sigmas_old), tol=tol)
                if debuginfo
                    println("latentgmm converged at $(iter_em)th iteration")
                end
                break
            end
        end
        if (iter_em == maxiteration) && (maxiteration > 50) && !taufixed && !dotest
            @warn("latentgmm not converge! $(iter_em), $(taufixed),
            $(ll), $(lldiff), $(wi), $(mu), $(sigmas), $(β)")
            println("latentgmm not converge! $(iter_em), $(taufixed),
            $(ll), $(lldiff), $(wi), $(mu), $(sigmas), $(β)")
        end
    end

    return(wi, mu, sigmas, β, ll)
end

function latentgmm(X::Matrix{Float64},
    Y::AbstractArray{Bool, 1}, groupindex::IntegerVector,
    ncomponent::Int; opts...)

    if ncomponent > 1
        wi_init, mu_init, sigmas_init, betas_init, ml_tmp = LatentGaussianMixtureModel.latentgmm(X, Y, groupindex, 1, ones(size(X)[2]), [1.0], [0.], [1.])
        gamma_init = LatentGaussianMixtureModel.predictgamma(X, Y, groupindex, wi_init, mu_init, sigmas_init, betas_init);
        mingamma = minimum(gamma_init) - 3*std(gamma_init)
        maxgamma = maximum(gamma_init) + 3*std(gamma_init)
        wi_init, mu_init, sigmas_init, ml_tmp = LatentGaussianMixtureModel.gmm(gamma_init, ncomponent)
    else
        betas_init = ones(size(X)[2])
        wi_init=[1.0;]
        mu_init=[0.0;]
        sigmas_init=[1.0;]
        mingamma = -Inf
        maxgamma = Inf
    end
    return LatentGaussianMixtureModel.latentgmm(X, Y, groupindex, ncomponent, betas_init, wi_init, mu_init, sigmas_init; mu_lb=ones(ncomponent).* mingamma, mu_ub=ones(ncomponent).* maxgamma, opts...)
end

function latentgmmrepeat(X::Matrix{Float64},
    Y::AbstractArray{Bool, 1}, groupindex::IntegerVector,
    ncomponent::Int; opts...)

    wi_init, mu_init, sigmas_init, betas_init, ml_tmp = LatentGaussianMixtureModel.latentgmm(X, Y, groupindex, 1, ones(size(X)[2]), [1.0], [0.], [1.])
    gamma_init = LatentGaussianMixtureModel.predictgamma(X, Y, groupindex, wi_init, mu_init, sigmas_init, betas_init);
    mingamma = minimum(gamma_init) - 3*std(gamma_init)
    maxgamma = maximum(gamma_init) + 3*std(gamma_init)

    wi_init, mu_init, sigmas_init, ml_tmp = LatentGaussianMixtureModel.gmm(gamma_init, ncomponent)
    return latentgmmrepeat(X, Y,
           groupindex, ncomponent, betas_init, wi_init,
           ones(ncomponent).*mingamma, ones(ncomponent).*maxgamma,
           0.25 .* sigmas_init, 2 .* sigmas_init; opts...)
end

"""
The interior step for `EMtest`
"""
function latentgmmrepeat(X::Matrix{Float64},
    Y::AbstractArray{Bool, 1}, groupindex::IntegerVector, C::Int,
    betas0::Vector{Float64}, wi_C1::Vector{Float64},
    mu_lb::Vector{Float64}, mu_ub::Vector{Float64},
    sigmas_lb::Vector{Float64}, sigmas_ub::Vector{Float64};
    taufixed::Bool=false, whichtosplit::Int64=1, tau::Float64=0.5,
    ntrials::Int=25, ngh::Int=100,
    sn::Vector{Float64}=sigmas_ub ./ 2, an=1/maximum(groupindex),
    debuginfo::Bool=false,
    gammaM::Vector = zeros(ngh*C),
    Wim::Matrix = zeros(maximum(groupindex), ngh*C),
    llN::Vector=zeros(length(Y)),
    llN2::Vector = zeros(length(Y)),
    llN3::Vector{Float64}=zeros(length(Y)),
    Xscratch::Matrix{Float64}=copy(X),
    xb::Vector=zeros(length(Y)), tol::Real=.005,
    pl::Bool=false, ptau::Bool=false, bn::Real=3.0 /maximum(groupindex),optional_EM_iteration::Bool=false, inparallel::Bool=false)

    n = maximum(groupindex)
    tau = min(tau, 1-tau)
    ghx, ghw = gausshermite(ngh)

    if inparallel
        wi = SharedArray(repeat(wi_C1, 1, 4*ntrials))
        mu = SharedArray(zeros(C, 4*ntrials))
        sigmas = SharedArray(ones(C, 4*ntrials))
        betas = SharedArray(repeat(betas0, 1, 4*ntrials))
        ml = SharedArray(-Inf .* ones(4*ntrials))
        @sync @distributed for i in 1:4*ntrials
            wi[:, i] = rand(Dirichlet(C, 1.0))
            mu[:, i] = rand(C) .* (mu_ub .- mu_lb) .+ mu_lb
            sigmas[:, i] = rand(C) .* (sigmas_ub .- sigmas_lb) .+ sigmas_lb

            wi[:, i], mu[:, i], sigmas[:, i], betas[:, i], ml[i] =
                 latentgmm(X, Y, groupindex, C, betas0,
                 wi[:, i], mu[:, i], sigmas[:, i],
                 whichtosplit=whichtosplit, tau=tau,
                 ghx=ghx, ghw=ghw, mu_lb=mu_lb, mu_ub=mu_ub,
                 maxiteration=16, sn=sn, an=an,
                 gammaM = gammaM, Wim=Wim,
                 llN=llN, llN2=llN2, llN3=llN3,
                 Xscratch=Xscratch, xb=xb,
                 Qmaxiteration=30, taufixed=taufixed, ngh=ngh,
                 dotest=false, tol=tol, bn=bn)
        end
        mlperm = sortperm(ml)
        @sync @distributed for j in 1:ntrials
            i = mlperm[4*ntrials+1 - j] # start from largest ml
            wi[:, i], mu[:, i], sigmas[:, i], betas[:, i], ml[i] =
                latentgmm(X, Y, groupindex, C, betas[:, i],
                wi[:, i], mu[:, i], sigmas[:, i],
                whichtosplit=whichtosplit, tau=tau, ghx=ghx, ghw=ghw,
                mu_lb=mu_lb,mu_ub=mu_ub, maxiteration=2000,
                sn=sn, an=an, gammaM = gammaM,
                Wim=Wim, llN=llN, llN2=llN2, llN3=llN3,
                Xscratch=Xscratch, xb=xb,
                Qmaxiteration=30, taufixed=taufixed, ngh=ngh,
                dotest=false, tol=tol, bn=bn)
        end
    else
        wi = repeat(wi_C1, 1, 4*ntrials)
        mu = zeros(C, 4*ntrials)
        sigmas = ones(C, 4*ntrials)
        betas = repeat(betas0, 1, 4*ntrials)
        ml = -Inf .* ones(4*ntrials)
        for i in 1:4*ntrials
            wi[:, i] = rand(Dirichlet(C, 1.0))
            mu[:, i] = rand(C) .* (mu_ub .- mu_lb) .+ mu_lb
            sigmas[:, i] = rand(C) .* (sigmas_ub .- sigmas_lb) .+ sigmas_lb

            wi[:, i], mu[:, i], sigmas[:, i], betas[:, i], ml[i] =
                 latentgmm(X, Y, groupindex, C, betas0,
                 wi[:, i], mu[:, i], sigmas[:, i],
                 whichtosplit=whichtosplit, tau=tau,
                 ghx=ghx, ghw=ghw, mu_lb=mu_lb, mu_ub=mu_ub,
                 maxiteration=16, sn=sn, an=an,
                 gammaM = gammaM, Wim=Wim,
                 llN=llN, llN2=llN2, llN3=llN3,
                 Xscratch=Xscratch, xb=xb,
                 Qmaxiteration=30, taufixed=taufixed, ngh=ngh,
                 dotest=false, tol=tol, bn=bn)
        end
        mlperm = sortperm(ml)
        for j in 1:ntrials
            i = mlperm[4*ntrials+1 - j] # start from largest ml
            wi[:, i], mu[:, i], sigmas[:, i], betas[:, i], ml[i] =
                latentgmm(X, Y, groupindex, C, betas[:, i],
                wi[:, i], mu[:, i], sigmas[:, i],
                whichtosplit=whichtosplit, tau=tau, ghx=ghx, ghw=ghw,
                mu_lb=mu_lb,mu_ub=mu_ub, maxiteration=2000,
                sn=sn, an=an, gammaM = gammaM,
                Wim=Wim, llN=llN, llN2=llN2, llN3=llN3,
                Xscratch=Xscratch, xb=xb,
                Qmaxiteration=30, taufixed=taufixed, ngh=ngh,
                dotest=false, tol=tol, bn=bn)
        end
    end

    mlmax, imax = findmax(ml[mlperm[(3*ntrials+1):4*ntrials]])
    imax = mlperm[3*ntrials+imax]

    if optional_EM_iteration
        re=latentgmm(X, Y, groupindex, C,
            betas[:, imax], wi[:, imax], mu[:, imax], sigmas[:, imax],
             maxiteration=2, an=an, sn=sn, debuginfo=false, ngh=ngh,
             mu_lb=mu_lb, mu_ub=mu_ub,
             tol=0., pl=pl, ptau=ptau, whichtosplit=whichtosplit, bn=bn)
        debuginfo && println("Trial:", re)
        return(re)
    else
        return(wi[:,imax], mu[:,imax], sigmas[:,imax], betas[:,imax], ml[imax])
    end
end

function EMtest(X::Matrix{Float64},
    Y::AbstractArray{Bool, 1}, groupindex::IntegerVector,
    C0::Int; vtau::Vector{Float64}=[.5;],
    ntrials::Int=25, ngh::Int=100, debuginfo::Bool=false,
    ctauparallel::Bool=true, tol::Real=0.001, bn::Real=3.0 / maximum(groupindex),  optional_EM_iteration::Bool=false, inparallel::Bool=false)

    C1 = C0 + 1
    n = maximum(groupindex)
    M = ngh * C1
    an1 = 1/n
    N,J=size(X)
    llN = zeros(N)
    llN2 = zeros(N)
    llN3 = zeros(N)
    Xscratch = copy(X)
    xb = zeros(N)

    #gamma_init, betas_init, sigmas_tmp = maxposterior(X, Y, groupindex)
    wi_init, mu_init, sigmas_init, betas_init, ml_C0 =
    latentgmm(X, Y, groupindex, 1, ones(J),
    [1.0], [0.], [1.], maxiteration=100, an=an1, sn=ones(C0), tol=.005)
    gamma_init = predictgamma(X, Y, groupindex,
        wi_init, mu_init, sigmas_init, betas_init)
    wi_init, mu_init, sigmas_init, ml_tmp = gmm(gamma_init, C0)
    mingamma = minimum(gamma_init) - 3*std(gamma_init)
    maxgamma = maximum(gamma_init) + 3*std(gamma_init)

    wi_init, mu_init, sigmas_init, betas_init, ml_C0 = latentgmmrepeat(X, Y,
       groupindex, C0, betas_init, wi_init,
       ones(C0).*mingamma, ones(C0).*maxgamma,
       0.25 .* sigmas_init, 2 .* sigmas_init,
       taufixed=false,
       ntrials=ntrials, ngh=ngh,
       sn=std(gamma_init).*ones(C0), an=an1,
       debuginfo=debuginfo,
       llN=llN, llN2=llN2, xb=xb, tol=tol,
       pl=false, ptau=false, bn=bn, inparallel=inparallel)

    # wi_init, mu_init, sigmas_init, betas_init, ml_C0 =
    #     latentgmm(X, Y, groupindex, C0, betas_init, wi_init, mu_init,
    #     sigmas_init, maxiteration=2000, an=an1,
    #     sn=std(gamma_init).*ones(C0), ngh=100, dotest=false, tol=.001,
    #     Qmaxiteration=5, pl=false, ptau=false)

    if C0 > 1
        trand=asymptoticdistribution(X, Y, groupindex, wi_init, mu_init, sigmas_init, betas_init)
    end

    if debuginfo
        println("ml_C0=", ml_C0)
    end

    or = sortperm(mu_init)
    wi0 = wi_init[or]
    mu0 = mu_init[or]
    sigmas0 = sigmas_init[or]
    betas0 = betas_init
    an = decidepenalty(wi0, mu0, sigmas0, n)

    N,J=size(X)
    gammaM = zeros(ngh*C1)
    Wim = zeros(n, ngh*C1)
    lr = 0.0

    lrv = zeros(length(vtau), C0)
    kl_tmp=zeros(length(vtau), C0)
    ll_tmp = zeros(length(vtau), C0)
    for whichtosplit in 1:C0, i in 1:length(vtau)

         #whichtosplit = mod1(irun, C0)
         #i = cld(irun, C0)
         ind = [1:whichtosplit; whichtosplit:C0;]
         if C1==2
             mu_lb = mingamma .* ones(2)
             mu_ub = maxgamma .* ones(2)
         elseif C1>2
             mu_lb = [mingamma; (mu0[1:(C0-1)] .+ mu0[2:C0])./2;]
             mu_ub = [(mu0[1:(C0-1)] .+ mu0[2:C0])./2; maxgamma;]
             mu_lb = mu_lb[ind]
             mu_ub = mu_ub[ind]
         end
         sigmas_lb = 0.25 .* sigmas0[ind]
         sigmas_ub = 2 .* sigmas0[ind]

         wi_C1 = wi0[ind]
         wi_C1[whichtosplit] = wi_C1[whichtosplit]*vtau[i]
         wi_C1[whichtosplit+1] = wi_C1[whichtosplit+1]*(1-vtau[i])
         mm_tmp = latentgmmrepeat(X, Y,
            groupindex, C1, betas0, wi_C1,
            mu_lb, mu_ub, sigmas_lb, sigmas_ub,
            taufixed=true, whichtosplit=whichtosplit, tau=vtau[i],
            ntrials=ntrials, ngh=ngh,
            sn=sigmas0[ind], an=an,
            debuginfo=debuginfo, gammaM = gammaM, Wim=Wim,
            llN=llN, llN2=llN2, xb=xb, tol=tol,
            pl=false, ptau=false, bn=bn, optional_EM_iteration=optional_EM_iteration, inparallel=inparallel)
         lrv[i, whichtosplit] = mm_tmp[5]
        if debuginfo
            show(mm_tmp)
            println(whichtosplit, " ", vtau[i], "->",
            lrv[i, whichtosplit])
        end
    end
    lr = maximum(lrv)

    if debuginfo
        println("lr=", lr)
    end
    Tvalue = 2*(lr - ml_C0)
    if C0 == 1
        pvalue = 1 - cdf(Chisq(2), Tvalue)
    else
        pvalue = mean(trand .> Tvalue)
    end
    if debuginfo
        println("The Test Statistic is ", Tvalue)
        println("The p vlaue of the test is ", pvalue)
    end
    return(Tvalue, pvalue, ml_C0, lr)
end