text:
# License for this file: MIT (expat)
# Copyright 2017-2018, DLR Institute of System Dynamics and Control
#
# This file is part of module
#   Modia3D.Composition (Modia3D/Composition/_module.jl)
#

"""
    get_Object3DsWithoutParent!(assembly,frames) - Returns all frames within assembly, that do not have a parent frame

The function traverses assembly and stores all frames without a parent frame in vector frames.
"""
function get_Object3DsWithoutParent!(assembly::Modia3D.AbstractAssembly, frames::AbstractVector)::NOTHING
    assemblyType = typeof(assembly)

    for i = 1:fieldcount(assemblyType)
          field     = getfield(assembly, fieldname(assemblyType,i))
          fieldType = typeof(field)
          if fieldType == Object3D
             if hasNoParent(field)
                push!(frames,field)
             end
          elseif fieldType <: Modia3D.AbstractAssembly
             get_Object3DsWithoutParent!(field, frames)
          end
    end
    return nothing
end

"""
    world = get_WorldObject3D(assembly) - Returns the assembly world frame

The function traverses assembly and collects all frames that do not have a parent frame.
If exactly one such frame is present, it is returned.
Otherwise, an error is triggered
"""
function get_WorldObject3D(assembly::Modia3D.AbstractAssembly)::Object3D
   frames = Object3D[]
   get_Object3DsWithoutParent!(assembly, frames)
   if length(frames) == 1
      return frames[1]
   elseif length(frames) < 1
      error("\nError message from Modia3D.SimulationModel!(..):\n",
            typeof(assembly), " has no reference frame.")
   else
      frameNames = String[]
      for frame in frames
         push!(frameNames, ModiaMath.fullName(frame))
      end
      error("\nError message from Modia3D.SimulationModel!(..):\n",
            typeof(assembly), " has the following frames that do not have a parent frame:\n",
            frameNames, "\n",
            "This is not allowed (exactly one frame must not have a parent frame).")
   end
end

"""
    realScalar = getFlangeVariable(assembly::Modia3D.AbstractAssembly, name::String)

Return the reference of a Revolute or Prismatic flange variable, given the top-level `assembly`
and the full path `name` of the flange variable. If one of the following flange variable names
is provided, a reference to the respective RealScalar is returned:
```
flange_b.phi
flange_b.s
flange_b.tau
flange_b.f
der(<full-name>.flange_b.phi)
der(<full-name>.flange_b.s)
```

If one of the following flange names is provided, `nothing` is returned:
```
flange_a.phi
flange_a.s
flange_a.tau
flange_a.f
der(<full-name>.flange_a.phi)
der(<full-name>.flange_a.s)
```
"""
function getFlangeVariable(assembly::Modia3D.AbstractAssembly, name::String)::Union{ModiaMath.RealScalar,Nothing}
    # Handle der(..)
    lenName = length(name)
    if lenName >= 4 && SubString(name,1,4) == "der("
        if SubString(name,lenName,lenName) != ")" || lenName < 6
            error("Wrong name = \"", name, "\"")
        end
        name2 = SubString(name,5,lenName-1)
        der   = true
    else
        name2 = name
        der   = false
    end

    field = assembly
    while true
        #println("... before. name2 = ", name2)
        i = findnext(".", name2, 1)
        #println("... after.")
        if i==nothing || length(name2) <= i[1] || i[1] <= 1
            error("Name not found: \"", name, "\"")
        end
        fieldName = Symbol(SubString(name2,1,i[1]-1))
        name2     = SubString(name2,i[1]+1)
        if !isdefined(field, fieldName)
            error("Name not found: \"", name, "\"")
        end
        field = getfield(field, fieldName)

        # Handle the names "flange_a" and "flange_b"
        #println("... field name = ", fieldName, ", typeof(field) = ", typeof(field))
        if name2=="flange_a.phi"
            if typeof(field) == Modia3D.Composition.TreeJointRevolute
                return nothing
            else
                error("Name = \"", name, "\" is no flange variable of a Revolute joint")
            end
        elseif name2=="flange_a.tau"
            if typeof(field) == Modia3D.Composition.TreeJointRevolute
                return nothing
            else
                error("Name = \"", name, "\" is no flange variable of a Revolute joint")
            end
        elseif name2=="flange_a.s"
            if typeof(field) == Modia3D.Composition.TreeJointPrismatic
                return nothing
            else
                error("Name = \"", name, "\" is no flange variable of a Prismatic joint")
            end
        elseif name2=="flange_a.f"
            if typeof(field) == Modia3D.Composition.TreeJointPrismatic
                return nothing
            else
                error("Name = \"", name, "\" is no flange variable of a Prismatic joint")
            end
        elseif name2=="flange_b.phi"
            if typeof(field) == Modia3D.Composition.TreeJointRevolute
                return der ? getfield(field,:w) : getfield(field,:phi)
            else
                error("Name = \"", name, "\" is no flange variable of a Revolute joint")
            end
        elseif name2=="flange_b.tau"
            if typeof(field) == Modia3D.Composition.TreeJointRevolute
                if der
                    error("Name = \"", name, "\" is no flange variable of a Revolute joint")
                end
                return getfield(field,:tau)
            else
                error("Name = \"", name, "\" is no flange variable of a Revolute joint")
            end
        elseif name2=="flange_b.s"
            if typeof(field) == Modia3D.Composition.TreeJointPrismatic
                return der ? getfield(field,:v) : getfield(field,:s)
            else
                error("Name = \"", name, "\" is no flange variable of a Prismatic joint")
            end
        elseif name2=="flange_b.f"
            if typeof(field) == Modia3D.Composition.TreeJointPrismatic
                if der
                    error("Name = \"", name, "\" is no flange variable of a Prismatic joint")
                end
                return getfield(field,:f)
            else
                error("Name = \"", name, "\" is no flange variable of a Prismatic joint")
            end
        end
    end
end

"""
    model = Modia3D.Model(assembly::Modia3D.AbstractAssembly;
                          analysis::ModiaMath.AnalysisType=ModiaMath.DynamicAnalysis,
                          potentialNames::Vector{String} = fill("",0),
                          flowNames::Vector{String}      = fill("",0))

Generate a `Modia3D.Model` from an `assembly` generated with macro [`Modia3D.@assembly`](@ref)
and the type of `analysis` to be carried out on the `assembly`.
"""
struct Model
    name::String
    nz::Int   # Number of event indicators
    assembly::Modia3D.AbstractAssembly
    var::ModiaMath.ModelVariables
    analysis::ModiaMath.AnalysisType
    x_start::Vector{Float64}
    x_fixed::Vector{Bool}
    x_nominal::Vector{Float64}
    is_constraint::Vector{Bool}

    # For connection with Modia
    potentials::Vector{Union{ModiaMath.RealScalar,Nothing}}
    flows::Vector{Union{ModiaMath.RealScalar,Nothing}}

    function Model(assembly::Modia3D.AbstractAssembly;
                   analysis::ModiaMath.AnalysisType=ModiaMath.DynamicAnalysis,
                   potentialNames::Vector{String} = fill("",0),
                   flowNames::Vector{String}      = fill("",0))

        name = Modia3D.trailingPartOfName( string( typeof(assembly) ) )
        world = get_WorldObject3D(assembly)
        assembly._internal.referenceObject3D = world
        # println("\n... world reference frame: ", ModiaMath.fullName(world))

        # Determine Variables in assembly
        var = ModiaMath.ModelVariables(assembly, analysis=analysis)
        #ModiaMath.print_ModelVariables(var)

        # Set initial values for x
        x_start   = zeros(var.nx)
        x_fixed   = fill(false,var.nx)
        x_nominal = fill(1.0  ,var.nx)
        ModiaMath.copy_start_to_x!(var, x_start, x_fixed, x_nominal)
        # println("... x0 = ", x_start)

        # Last nfc equations are the constraint equations
        is_constraint = fill(false, var.nx)
        for i = (var.nx-var.nfc+1):var.nx
            is_constraint[i] = true
        end

        # Construct Scene(..) object
        so = assembly._internal.sceneOptions
        sceneOptions::SceneOptions = typeof(so) == NOTHING ? SceneOptions() : so
        scene = Scene(sceneOptions)
        scene.analysis = analysis
        assembly._internal.scene = scene
        if !scene.options.enableContactDetection
            scene.collide = false
        end

        # Build tree for optimized structure or standard structure
        # collision handling is only available for optimized structure
        nz = 0
        if scene.options.useOptimizedStructure
            build_superObjs!(scene, world)
            if !scene.options.enableContactDetection
               scene.collide = false
            end

            if scene.collide
                initializeContactDetection!(world, scene)
                if scene.collide
                    nz = scene.options.contactDetection.contactPairs.nz
                    append!(scene.allVisuElements, world.contactVisuObj1)
                    append!(scene.allVisuElements, world.contactVisuObj2)
                    append!(scene.allVisuElements, world.supportVisuObj1A)
                    append!(scene.allVisuElements, world.supportVisuObj1B)
                    append!(scene.allVisuElements, world.supportVisuObj1C)
                    append!(scene.allVisuElements, world.supportVisuObj2A)
                    append!(scene.allVisuElements, world.supportVisuObj2B)
                    append!(scene.allVisuElements, world.supportVisuObj2C)
                end
            end
            initializeMassComputation!(scene)
        else
            build_tree!(scene, world)
            if scene.options.enableContactDetection
                error("Collision handling is only possible with the optimized structure. Please set useOptimizedStructure = true in Modia3D.SceneOptions.")
            end
        end

        # Initialize connections between signals and frames, joints, ...
        build_SignalObject3DConnections!(assembly)
        scene.initAnalysis = true

        # Generate potentials and flows
        potentials = Vector{Union{ModiaMath.RealScalar,Nothing}}(nothing, length(potentialNames))
        flows      = Vector{Union{ModiaMath.RealScalar,Nothing}}(nothing, length(flowNames))
        for i = 1:length(potentialNames)
            potentials[i] = getFlangeVariable(assembly, potentialNames[i])
            #println("... potentials[", i, "] = ", potentials[i]==nothing ? "nothing" : ModiaMath.instanceName(potentials[i]))
        end
        for i = 1:length(flowNames)
            flows[i] = getFlangeVariable(assembly, flowNames[i])
            #println("... flows[", i, "] = ", flows[i]==nothing ? "nothing" : ModiaMath.instanceName(flows[i]))
        end

        # Generate Model object
        new(name, nz, assembly, var, analysis, x_start, x_fixed, x_nominal, is_constraint, potentials, flows)
   end
end

"""
    simModel = SimulationModel(assembly::Modia3D.AbstractAssembly;
                               analysis::ModiaMath.dAnalysisType=ModiaMath.DynamicAnalysis,
                               startTime = 0.0, stopTime  = 1.0, tolerance = 1e-4,
                               interval  = (stopTime-startTime)/500.0,
                               maxStepSize = NaN, maxNumberOfSteps=missing)

Generate a `simulationModel` from an `assembly` generated with macro [`Modia3D.@assembly`](@ref)
and the type of `analysis` to be carried out on the `assembly`.
Additionally, default `startTime`, `stopTime`, `tolerance`, `interval`,
`maxStepSize`, `maxNumberOfSteps`, for the simulation
engine are defined. These values should be adapted so that assembly-specific, meaningful
defaults are provided.
"""
struct SimulationModel <: ModiaMath.AbstractSimulationModel
    modelName::String
    simulationState::ModiaMath.SimulationState
    var::ModiaMath.ModelVariables
    model::Model

    function SimulationModel(assembly::Modia3D.AbstractAssembly;
                             analysis::ModiaMath.AnalysisType=ModiaMath.DynamicAnalysis,
                             startTime = 0.0,
                             stopTime  = 1.0,
                             tolerance = 1e-4,
                             interval  = (stopTime-startTime)/500.0,
                             maxStepSize=NaN,
                             maxNumberOfSteps=missing,
                             hev = 1e-8,
                             scaleConstraintsAtEvents::Bool = true)
        model           = Model(assembly; analysis=analysis)
        simulationState = ModiaMath.SimulationState(model.name, getModelResidues!, model.x_start, ModiaMath.Variables.getVariableName;
                                x_fixed          = model.x_fixed,
                                x_nominal        = model.x_nominal,
                                is_constraint    = model.is_constraint,
                                nz               = model.nz,
                                defaultStartTime = startTime,
                                defaultStopTime  = stopTime,
                                defaultTolerance = tolerance,
                                defaultInterval  = interval,
                                defaultMaxStepSize = maxStepSize,
                                defaultMaxNumberOfSteps = maxNumberOfSteps,
                                getResultNames   = Composition.getResultNames,
                                getResult        = ModiaMath.Variables.getResult,
                                storeResult!     = ModiaMath.Variables.storeVariables!,
                                hev = hev,
                                scaleConstraintsAtEvents = scaleConstraintsAtEvents)

#=
      scene = assembly._internal.scene
      if useOptimizedStructure
          println("\n... SimulationModel: treeAccVelo begin")
          tree = scene.treeAccVelo
      else
          println("\n... SimulationModel: tree begin")
          tree = scene.tree
      end

      for obj in tree
         if hasParent(obj)
            println(ModiaMath.instanceName(obj), " (parent = ", ModiaMath.instanceName(obj.parent), ")")
            println("    r_rel = ", obj.r_rel, ", R_rel = ", obj.R_rel)
         else
            println(ModiaMath.instanceName(obj), " (no parent)")
         end
         if hasJoint(obj)
            println("    joint type = ", typeof(obj.joint))
         end
         if objectHasMass(obj)
            println("    m = ", obj.massProperties.m, " rCM = ", obj.massProperties.rCM, " I = ", obj.massProperties.I)
         end
      end
      println("...  end\n\n")

      println("... allVisuElements:")
      for obj in scene.allVisuElements
          println(ModiaMath.instanceName(obj))
      end
      println("... end allVisuElements")
=#
       new(model.name, simulationState, model.var, model)
   end
end

ModiaMath.print_ModelVariables(model::Model) = ModiaMath.print_ModelVariables(model.var)
print_ModelVariables(          model::Model) = ModiaMath.print_ModelVariables(model.var)

ModiaMath.print_ModelVariables(simModel::SimulationModel) = ModiaMath.print_ModelVariables(simModel.model.var)
print_ModelVariables(          simModel::SimulationModel) = ModiaMath.print_ModelVariables(simModel.model.var)

getResultNames(simModel::SimulationModel) = simModel.model.var.result_names

function initializeFlowVariables(scene::Scene)
  flowVarInput   = scene.flowVarInput
  flowVarOutput  = scene.flowVarOutput
  @assert(length(flowVarInput) == length(flowVarOutput))
  for i=eachindex(flowVarInput)
    flowVarInput[i].value  = 0.0
    flowVarOutput[i].value = 0.0
  end
end

function setFlowVariables(scene::Scene)
  flowVarInput   = scene.flowVarInput
  flowVarOutput  = scene.flowVarOutput
  @assert(length(flowVarInput) == length(flowVarOutput))
  for i=eachindex(flowVarInput)
    flowVarInput[i].value  += flowVarOutput[i].value
  end
end

function setPotentialVariables(scene::Scene)
  potentialVarInput   = scene.potentialVarInput
  potentialVarOutput  = scene.potentialVarOutput
  @assert(length(potentialVarInput) == length(potentialVarOutput))
  for i=eachindex(potentialVarInput)
    potentialVarInput[i].value = potentialVarOutput[i].value
  end
end

function computationForcesAndTorques(scene::Scene, sim::ModiaMath.SimulationState)
    for forceTorque in scene.uniqueForceTorques
        computeTorque(forceTorque, sim)
    end
end

function computationSignals(scene::Scene, sim::ModiaMath.SimulationState)
    for signal in scene.uniqueSignals
        computeSignal(signal, sim)
    end
end

# Return a table of actual variable and residue values from nonlinear solver in case of error
function ModiaMath.getVariableAndResidueValues(simModel::SimulationModel)
   var = simModel.model.var
   v_table = ModiaMath.get_variableValueTable(var)
   r_table = ModiaMath.get_residueValueTable(var, simModel.simulationState.derxev)
   return (v_table, r_table)
end

# Only temporarily here. Shall be moved to ModiaMath
"""
    copy_x_to_variables!(x::Vector{Float64}, vars::ModelVariables)

Copy `x` of the integrator interface to the model variables `vars`.
"""
function copy_x_to_variables!(x::Vector{Float64}, m::ModiaMath.ModelVariables)::Nothing
    @assert(length(x) == m.nx)

    for v in m.x_var
        #println("... typeof(", fullName(v), ") = ", typeof(v), ", isimmutable(v) = ", isimmutable(v))
        if ModiaMath.isScalar(v)
            v.value = x[ v.ivar ]
        elseif isimmutable(v.value)
            # v is an immutable array (e.g. SVector{3,Float64})
            v.value = x[ v.ivar:v.ivar + length(v.value) - 1 ]
        else
            vv = v.value
            for j in 1:length(vv)
                vv[j] = x[ v.ivar + j - 1 ]
            end
        end
    end
    return nothing
end

"""
    copy_derx_to_variables!(time::Float64, derx::Vector{Float64}, vars::ModelVariables)

Copy `time` and `derx` of the integrator interface to the model variables `vars`.
"""
function copy_derx_to_variables!(time::Float64, derx::Vector{Float64}, m::ModiaMath.ModelVariables)::Nothing
    @assert(length(derx) == m.nx)

    m.var[1].value = time
    for v in m.derx_var
        if ModiaMath.isScalar(v)
            v.value = derx[ v.ivar ]
        elseif isimmutable(v.value)
            # v is an immutable array (e.g. SVector{3,Float64})
            v.value = derx[ v.ivar:v.ivar + length(v.value) - 1 ]
        else
            vv = v.value
            for j in 1:length(vv)
                vv[j] = derx[ v.ivar + j - 1 ]
            end
        end
    end

    return nothing
end

"""
    modelf1!(model, x, potentials)

Copy the potential variables defined for `model::Modia3D.Model` from `x` to vector `potentials`,
"""
function model_f1!(model::Model, x::Vector{Float64}, potentials::Vector{Float64})::Nothing
    @assert( length(model.potentials) == length(potentials) )

    # Copy (time,x) to model variables
    copy_x_to_variables!(x, model.var)

    # Copy selected potential variables from model to output argument potentials
    for i = 1:length(potentials)
        potentials[i] = model.potentials[i] == nothing ? 0.0 : model.potentials[i].value
    end
    return nothing
end

const str_DUMMY = "dummyDistance(nothing,nothing)"
const zeroVector = Vector{Float64}()

getModelResidues!(simModel::SimulationModel, time::Float64, _x::Vector{Float64}, _derx::Vector{Float64}, _r::Vector{Float64}, _w::Vector{Float64}) =
    model_f2!(simModel.model, simModel.simulationState, time, _x, _derx, zeroVector, _r, _w)

function model_f2!(m::Model, sim::ModiaMath.SimulationState,
                   time::Float64, _x::Vector{Float64}, _derx::Vector{Float64}, _flows::Vector{Float64},
                   _r::Vector{Float64}, _w::Vector{Float64})
   # println("... time = ", time, ", x = ", _x, ", derx = ", _derx)
   world = m.assembly._internal.referenceObject3D
   scene = m.assembly._internal.scene
   if scene.options.useOptimizedStructure
      tree          = scene.treeAccVelo
   else
      tree          = scene.tree
   end
   cutJoints        = scene.cutJoints
   var              = m.var
   analysis         = m.analysis

   # Handle initialization and termination of model
   if ModiaMath.isInitial(sim)
      # println("... isInitial = true")
      if scene.collide
         initializeContactDetection!(world, scene)
         append!(scene.allVisuElements, world.contactVisuObj1)
         append!(scene.allVisuElements, world.contactVisuObj2)

         append!(scene.allVisuElements, world.supportVisuObj1A)
         append!(scene.allVisuElements, world.supportVisuObj1B)
         append!(scene.allVisuElements, world.supportVisuObj1C)
         append!(scene.allVisuElements, world.supportVisuObj2A)
         append!(scene.allVisuElements, world.supportVisuObj2B)
         append!(scene.allVisuElements, world.supportVisuObj2C)
      end
      if scene.visualize
         initializeVisualization(Modia3D.renderer[1], scene.allVisuElements)
      end
   end

   if ModiaMath.isTerminal(sim)
      # println("... isTerminal = true")
      if scene.visualize
         # println("... visualization closing")
         closeVisualization(Modia3D.renderer[1])
      end
      if scene.collide
         closeContactDetection!(m.assembly)
      end
      return
   end

   # Check input arguments
   @assert(length(_x)    == var.nx)
   @assert(length(_derx) == var.nx)
   @assert(length(_r)    == var.nx)
   @assert(length(_w)    == 0)
   @assert(length(_flows) == length(m.flows))

   storeResults = ModiaMath.isStoreResult(sim) && (scene.visualize || var.nwr > 0 || var.nwc > 0 )

   # Initialize residues _r with zero
   for i in eachindex(_r)
      _r[i] = 0.0
   end

   # Copy x and derx to variables
   if length(m.potentials) == 0
       ModiaMath.copy_x_and_derx_to_variables!(time, _x, _derx, var)
   else
       copy_derx_to_variables!(time, _derx, var)
   end

   # Compute signals
   initializeFlowVariables(scene)
   computationSignals(scene, sim)
   setPotentialVariables(scene)
   setFlowVariables(scene)

   # Copy flows to model variables
   for i = 1:length(_flows)
       if m.flows[i] != nothing
          m.flows[i].value = _flows[i]
       end
   end

   # Initialize force/torque of world-frame
   world.dynamics.f = ModiaMath.ZeroVector3D
   world.dynamics.t = ModiaMath.ZeroVector3D

   # Compute positions, velocities, accelerations in a forward recursion

   open("log.txt", "a") do file

   #println("\n... Forward recursion")
   for obj in tree
      # Compute kinematics
      computeKinematics!(obj.joint, obj, analysis, time)

      if typeof(obj.visualizationFrame) != NOTHING
         obj.visualizationFrame.r_abs = obj.r_abs
         obj.visualizationFrame.R_abs = obj.R_abs
      end

      # Compute forces and torques
      initializeFlowVariables(scene)
      setPotentialVariables(scene)
      computationForcesAndTorques(scene,sim)
      setFlowVariables(scene)

      # Initialize forces and torques
      dynamics::Object3Ddynamics = obj.dynamics

      if analysis != ModiaMath.KinematicAnalysis
         hasMassProperties = false
         if scene.options.useOptimizedStructure
            if objectHasMass(obj)
               massProperties    = obj.massProperties
               hasMassProperties = true
            end
         else
            if dataHasMass(obj)
               massProperties    = obj.data.massProperties
               hasMassProperties = true
            end
         end

         if hasMassProperties
            # Compute inertia forces / torques
            mass = massProperties.m
            rCM  = massProperties.rCM
            I    = massProperties.I
            w    = dynamics.w
            #println("w = ",w)
            grav = gravityAcceleration(scene.options.gravityField, obj.r_abs)
            # println("grav = ", grav)
            if rCM === ModiaMath.ZeroVector3D
               dynamics.f = -mass*( obj.R_abs*(dynamics.a0  - grav) )
               dynamics.t = -(I*dynamics.z + cross(w, I*w))
               # println("(dynamics.a0  - grav) = ", (dynamics.a0  - grav) )
               #println("dynamics.t = ", dynamics.t)
               #println(" ")
            else
               dynamics.f = -mass*( obj.R_abs*(dynamics.a0 - grav) +
                                  cross(dynamics.z, rCM) + cross(w, cross(w, rCM)))
               dynamics.t = -(I*dynamics.z + cross(w, I*w)) + cross(rCM, dynamics.f)
            end

         else
            dynamics.f = ModiaMath.ZeroVector3D
            dynamics.t = ModiaMath.ZeroVector3D
         end
      end
      #println("    ", ModiaMath.instanceName(obj), ": f = ", dynamics.f, ", t = ", dynamics.t)

      #if ModiaMath.isInitial(sim)
      #   println(ModiaMath.instanceName(obj), ": f = ", dynamics.f, ", t = ", dynamics.t)
      #end
   end # for

   # Compute contact forces/torques
   if scene.collide
      ch = scene.options.contactDetection
      chpairs = ch.contactPairs

      # Compute signed distances of all contact shapes during zero-crossing computation
      setComputationFlag(ch)
      if ModiaMath.isEvent(sim)    # with Event
         selectContactPairsWithEvent!(sim, ch, world) #sim#
      elseif ModiaMath.isZeroCrossing(sim) # no Event
         selectContactPairsNoEvent!(sim, ch, world) #sim
      else
         getDistances!(ch, world)
      end

      # Handle zero crossing event
      simh = sim.eventHandler
      for (pairID, pair) in ch.contactDict
        obj1 = pair.obj1
        obj2 = pair.obj2
        rContact      = (pair.contactPoint1 + pair.contactPoint2)/2.0
        contactNormal = pair.contactNormal

        #simh = sim.eventHandler
        #println( "time = ", sim.time, ": ", ModiaMath.instanceName(obj1), " ", ModiaMath.instanceName(obj2),
        #         " i = ", i, ", initial = ", simh.initial, ", event = ", simh.event, " change = ", chpairs.changeDirection[i] )
        if ModiaMath.isEvent(sim)
            # Include contact pair material into collision pair
            if haskey(ch.lastContactDict, pairID)
                # use material (reference) from previous event
                pair.contactPairMaterial = ch.lastContactDict[pairID].contactPairMaterial    # improve later (should avoid to inquire pairID twice)
            else
                # determine contact pair material
                pair.contactPairMaterial = contactStart(obj1, obj2, rContact, contactNormal,
                                                         scene.options.elasticContactReductionFactor)
                simh.restart = max(simh.restart, ModiaMath.Restart)
                simh.newEventIteration = false
                if ModiaMath.isLogEvents(simh.logger)
                    name1 = ModiaMath.instanceName(obj1)
                    name2 = ModiaMath.instanceName(obj2)
                    n     = contactNormal
                    println("        distance(", name1, ",", name2, ") = ", pair.distanceWithHysteresis, " became < 0")
                    Printf.@printf("            contact normal = [%.3g,%.3g,%.3g], contact position = [%.3g,%.3g,%.3g], c_res=%.3g, d_res=%.3g\n",
                                   n[1],n[2],n[3],rContact[1],rContact[2],rContact[3], pair.contactPairMaterial.c_res, pair.contactPairMaterial.d_res)
                end
            end
        end

        # println("length(ch.dictCommunicate) ", length(ch.dictCommunicate) )
        (f1,f2,t1,t2) = responseCalculation(pair.contactPairMaterial, obj1, obj2, rContact, contactNormal,
                                             pair.distanceWithHysteresis, time, file)

        # Transform forces/torques in local part frames
        obj1.dynamics.f += obj1.R_abs*f1
        obj1.dynamics.t += obj1.R_abs*t1
        obj2.dynamics.f += obj2.R_abs*f2
        obj2.dynamics.t += obj2.R_abs*t2
#=
               if time > 0.785 && time < 0.787  && String(ModiaMath.instanceName(obj1)) != "table.box1"
                  println("obj1= \"", ModiaMath.instanceName(obj1), "\" obj2 = ", ModiaMath.instanceName(obj2), " f = ", obj1.dynamics.f, " t = ", obj1.dynamics.t, " rContact = ", rContact, " ctNormal[i] ", ModiaMath.Vector3D(chpairs.contactNormal[i]) ," time = ", time)
               end
=#
      end

      if ModiaMath.isEvent(sim)
            # Should be possible to make this more efficient
            for (pairID, pair) in ch.lastContactDict
                if !haskey(ch.contactDict, pairID)
                    contactEnd(pair.contactPairMaterial, pair.obj1, pair.obj2)
                    simh.restart = max(simh.restart, ModiaMath.Restart)
                    simh.newEventIteration = false
                    if !ModiaMath.isLogEvents(simh.logger)
                        break
                    end
                    name1 = ModiaMath.instanceName(pair.obj1)
                    name2 = ModiaMath.instanceName(pair.obj2)
                    println("        distance(", name1, ",", name2, ")  became > 0")
                end
            end

            # Save contactDict in lastContactDict
            empty!(ch.lastContactDict)
            for (pairID, pair) in ch.contactDict
                push!(ch.lastContactDict, pairID => pair)
            end
      end
   end

   # For all cut-joints, set forces/torques at the frames and compute cut-joint residues
   for joint in cutJoints
      computePositionResidues!(joint, time)
      if analysis == ModiaMath.DynamicAnalysis
         computeVelocityResidues!(joint, time)
      end
      if analysis != ModiaMath.KinematicAnalysis
         computeCutForcesAndToques!(joint, time)
      end
   end

   #println("\n... Backward recursion")

   # Compute forces/torques and residues in a backward recursion
   if analysis != ModiaMath.KinematicAnalysis
      for i = length(tree):-1:1
         obj = tree[i]
         computeForceTorqueAndResidue!(obj.joint, obj, analysis, time)
      end

      #if ModiaMath.isInitial(sim)
      #   for i = length(tree):-1:1
      #      obj = tree[i]
      #      println(ModiaMath.instanceName(obj), ": f = ", obj.dynamics.f, ", t = ", obj.dynamics.t)
      #   end
      #end
   end

   # Visualize at a communication point

   if scene.visualize && storeResults
      # Compute positions of frames that are only used for visualization
      if scene.options.useOptimizedStructure
          for obj in scene.treeVisu
              parent = obj.parent
              obj.r_abs = obj.r_rel ≡ ModiaMath.ZeroVector3D ? parent.r_abs : parent.r_abs + parent.R_abs'*obj.r_rel
              obj.R_abs = obj.R_rel ≡ ModiaMath.NullRotation ? parent.R_abs : obj.R_rel*parent.R_abs
              if typeof(obj.visualizationFrame) != NOTHING
                 obj.visualizationFrame.r_abs = obj.r_abs
                 obj.visualizationFrame.R_abs = obj.R_abs
              end
          end
      end
      visualize!(Modia3D.renderer[1], time)
   end

   # Copy variables to residues
   ModiaMath.copy_variables_to_residue!(var, _x, _derx, _r)

   #println("... x = ", _x, ", derx = ", _derx, ", residue = ", _r)

end   # from open
   return nothing
end


text:
include("../imports.jl")

@testset "TabularRegression" begin
    df = DataFrame(A = 1:4, B = ["M", "F", "F", "M"], C = 10:13)
    td = TableDataset(df)
    targets = [rand(2) for _ in 1:4]
    method = TabularRegression(2, td; catcols=(:B,), contcols=(:A,))
    testencoding(method.encodings, method.blocks)
    DLPipelines.checkmethod_core(method)
    @test_nowarn methodlossfn(method)
    @test_nowarn methodmodel(method)

    @testset "`encodeinput`" begin
        row = mockblock(method.blocks[1])

        xtrain = encodeinput(method, Training(), row)
        @test length(xtrain[1]) == length(method.blocks[1].catcols)
        @test length(xtrain[2]) == length(method.blocks[1].contcols)

        @test eltype(xtrain[1]) <: Number
    end

    @testset "`encodetarget`" begin
        target = 11
        y = encodetarget(method, Training(), target)
        @test target == y
    end

    @test_nowarn method = TabularRegression(2, td)

    FastAI.test_method_show(method, ShowText(Base.DevNull()))
end

ext: jl


text:
cd(@__DIR__)
using AdvancedPS, Random, Parameters, Distributions

priors = (
    θ=(Beta(1, 1),),
)

bounds = ((0,1),)

N = 10
k = rand(Binomial(N, .5))
data = (N=N,k=k)

function loglike(θ, data)
    @unpack N,k = data
    n_sim = 10^4
    counter(_) = rand(Binomial(N, θ)) == k ? (return 1) : (return 0)
    cnt = mapreduce(counter, +, 1:n_sim)
    return log(cnt/n_sim)
end

# loglike(θ, data) = logpdf(Binomial(data.N, θ), data.k)

loglike(θ) = loglike(θ..., data)

model = DEModel(priors=priors, model=loglike)

de = DE(bounds=bounds, burnin=1000, priors=priors, σ=.01)
n_iter = 2000
@elapsed chains = sample(model, de, MCMCThreads(), n_iter, progress=true)
println(chains)


text:
push!(LOAD_PATH, "/home/ian/Grad/Research/Disks/DiskJockey/movie/")
# push!(LOAD_PATH, "/n/home07/iczekala/DiskJockey/")
# push!(LOAD_PATH, "/pool/scout0/DiskJockey/")

using ArgParse

s = ArgParseSettings()
@add_arg_table s begin
    # "--opt1"
    # help = "an option with an argument"
    "--run_index", "-r"
    help = "Output run index"
    default = 0
    arg_type = Int
end

parsed_args = parse_args(ARGS, s)
run_index = parsed_args["run_index"]

# Determine the proper vmin and vmax scaling for the dust images, channel maps,
# and transforms of channel maps
using consonance

using DiskJockey.image
using DiskJockey.visibilities

import YAML
config = YAML.load(open("config.yaml"))

# Store a global max value, and then go through each image, see if it beats
# this value, if so, record the new one, then close image to keep the memory
# low.

imgdir = moviedir * @sprintf("g%02d/", run_index)

fname = imgdir * @sprintf("gimage%04d.out", 1)
im = imread(fname)

dpc = config["parameters"]["dpc"][1]
skim = imToSky(im, dpc)
nchan = length(skim.lams)
println("$nchan channels")

global vmax_g = maximum(skim.data)
global vmin_g = minimum(skim.data)

vis_fft = transform(skim, 1)

global vmax_v = maximum(abs(vis_fft.VV))
global vmin_v = minimum(abs(vis_fft.VV))

# Load all of the images into memory.
# Determine which image*.out files are in this directory.
imgfunc = x -> contains(x, "gimage") && contains(x, ".out")
imglist = filter(imgfunc, readdir(imgdir))

nimg = length(imglist)

# Go through each one, determine the maximum and minimum, and if it beats the
# current value, replace it
for fname in imglist
    im = imread(imgdir * fname)
    skim = imToSky(im, dpc)
    min,max = extrema(skim.data)
    if max > vmax_g
        vmax_g = max
    end
    if min < vmin_g
        vmin_g = min
    end
    for i=1:nchan
        vis_fft = transform(skim, i)
        min,max = extrema(abs(vis_fft.VV))
        if max > vmax_v
            vmax_v = max
        end
        if min < vmin_v
            vmin_v = min
        end
    end

end

println("Gas minimum and maximum $vmin_g, $vmax_g")
println("Visibility minimum and maximum $vmin_v, $vmax_v")

max_stars_repo_name: elnjensen/DiskJockey

max_issues_repo_name: elnjensen/DiskJockey

max_forks_repo_name: elnjensen/DiskJockey


text:
export elbo, em

function elbo(ks::KalmanSmoother,z,θ,build)
    T = length(z)

    mθ = build(θ)

    μ = mθ.μ
    Σ0 = mθ.Σ0
    F = mθ.F
    Q = mθ.Q
    Θ = mθ.Θ
    Qp = Θ*Q*Θ'
    H = mθ.H
    R = mθ.R

    S11 = ks.S11
    S00 = ks.S00
    S10 = ks.S10

    Syy = ks.Syy

    Sμ = ks.P[1] + (ks.x[1]-μ)*(ks.x[1]-μ)'

    l1 = logdet(Σ0) + tr(Σ0\Sμ)
    l2 = T * logdet(Q) + tr(pinv(Qp)*(S11 - F*S10' - S10*F' + F*S00*F'))
    l3 = T * logdet(R) + tr(pinv(R)*Syy)

    l1 + l2 + l3
end

function elbo(m::LinearStateSpaceModel,z,θ,build)
    kf = kalman_filter(m,z)
    ks = kalman_smoother(kf,m)

    elbo(ks,z,θ,build)
end

# This is not ready for prime time
function elbo_gradient(ks::KalmanSmoother,z,θ)
    T = length(z)

    mθ = SSTFI(θ)
    μ = mθ.μ
    Σ0 = mθ.Σ0
    F = mθ.F
    Q = mθ.Q
    Θ = mθ.Θ
    Qp = Θ*Q*Θ'
    H = mθ.H
    R = mθ.R

    dm = Filtering._gradient(θ)
    dF = dm.dF
    dQ = dm.dQ
    dH = dm.dH
    dR = dm.dR
    dΘ = dm.dΘ
    dμ = dm.dμ
    dΣ0 = dm.dΣ0

    S11 = ks.S11
    S00 = ks.S00
    S10 = ks.S10

    Syy = ks.Syy

    Sμ = ks.P[1] + (ks.x[1]-μ)*(ks.x[1]-μ)'

    # l1
    l1 = logdet(Σ0) + tr(Σ0\Sμ)
    dl1 = zeros(length(θ0))

    # l2
    Θp = pinv(Θ)
    Sf =(S11 - F*S10' - S10*F' + F*S00*F')
    l2 = T*logdet(Q) + tr(pinv(Qp)*Sf)
    dl2 = [T*tr(Q\dQ[:,:,i]) +
        tr(-Θp'*(Q\dQ[:,:,i])*(Q\Θp)*Sf +
            pinv(Qp)*(-dF[:,:,i]*S10' - S10*dF[:,:,i]' + dF[:,:,i]*S00*F' + F*S00*dF[:,:,i]'))
        for i in 1:length(θ)]

    # l3
    l3 = T*logdet(R) + tr(R\Syy)
    dl3 = [T*tr(R\dR[:,:,i]) + tr(-(R\dR[:,:,i])*(R\Syy)) for i in 1:length(θ)]

    dl1 + dl2 + dl3
end

function em_step(build,θ,z)
    opt = optimize(t->elbo(build(θ),z,t,build),θ)
    Optim.minimizer(opt)
end

function em_step2(build,θ,z;optim_options...)
    m = build(θ)
    kf = kalman_filter(m,z)
    ks = kalman_smoother(kf,m)
    opt = optimize(t->elbo(ks,z,t,build),θ,Optim.Options(;optim_options...))
    Optim.minimizer(opt)
end

struct EMResults
    θs
    ls
end

function em(build,θ0,z;tol=1e-4,max_iterations=1000,verbose=false,optim_options...)
    m0 = build(θ0)
    kf0 = kalman_filter(m0,z)
    ks0 = kalman_smoother(kf0,m0)
    ls = zeros(max_iterations)
    ls[1] = filter_likelihood(kf0)
    if verbose
            println("Iteration: ",1)
            println("Log likelihood: ",ls[1])
    end

    θs = zeros(length(θ0),max_iterations)
    θs[:,1] = θ0[:]
    s = 1
    for i in 2:max_iterations
        θs[:,i] = em_step2(build,θs[:,i-1],z;optim_options...)
        ls[i] = filter_likelihood(kalman_filter(build(θs[:,i]),z))
        if verbose
            println("Iteration: ",i)
            println("Log likelihood: ",ls[i])
        end
        if ls[i] > ls[i-1]
            @warn "Expectation maximization increased the log likelihood"
            break
        end
        s+=1
        if (ls[i-1] - ls[i])<tol
            break
        end
    end
    EMResults(θs[:,1:s],ls[1:s])
end


text:
# Piracy that used to be requires, allowing Zyogote.jl to be specialized for SciML

function ∇tmap(cx, f, args...)
    ys_and_backs = SciMLBase.tmap((args...) -> Zygote._pullback(cx, f, args...), args...)
    if isempty(ys_and_backs)
        ys_and_backs, _ -> (NoTangent(), NoTangent())
    else
        ys, backs = Zygote.unzip(ys_and_backs)
        function ∇tmap_internal(Δ)
            Δf_and_args_zipped = SciMLBase.tmap((f, δ) -> f(δ), backs, Δ)
            Δf_and_args = Zygote.unzip(Δf_and_args_zipped)
            Δf = reduce(Zygote.accum, Δf_and_args[1])
            (Δf, Δf_and_args[2:end]...)
        end
        ys, ∇tmap_internal
    end
end

function ∇responsible_map(cx, f, args...)
    ys_and_backs = SciMLBase.responsible_map((args...) -> Zygote._pullback(cx, f, args...), args...)
    if isempty(ys_and_backs)
        ys_and_backs, _ -> (NoTangent(), NoTangent())
    else
        ys, backs = Zygote.unzip(ys_and_backs)
        ys, function ∇responsible_map_internal(Δ)
            # Apply pullbacks in reverse order. Needed for correctness if `f` is stateful.
            Δf_and_args_zipped = SciMLBase.responsible_map((f, δ) -> f(δ), Zygote._tryreverse(SciMLBase.responsible_map, backs, Δ)...)
            Δf_and_args = Zygote.unzip(Zygote._tryreverse(SciMLBase.responsible_map, Δf_and_args_zipped))
            Δf = reduce(Zygote.accum, Δf_and_args[1])
            (Δf, Δf_and_args[2:end]...)
        end
    end
end

ZygoteRules.@adjoint function SciMLBase.tmap(f, args::Union{AbstractArray,Tuple}...)
    ∇tmap(__context__, f, args...)
end

ZygoteRules.@adjoint function SciMLBase.responsible_map(f, args::Union{AbstractArray,Tuple}...)
    ∇responsible_map(__context__, f, args...)
end


text:
using LatticeDiracOperators
using Gaugefields
using Test
using LinearAlgebra

@testset "LatticeDiracOperators.jl" begin
    @testset "Wilson HMC" begin
        println("Wilson HMC")
        include("wilsonhmc.jl")
        @test true
    end

    @testset "Staggered HMC" begin
        println("Staggered HMC")
        include("hmc.jl")
        @test true
    end

    @testset "Basic operations" begin
        include("basic.jl")
    # Write your tests here.
        @test true
    end
end


text:
using NNlib: kernel_size,  output_size, channels_in, channels_out
using NNlib: ∇conv_data!, ∇conv_filter!, conv!

function QMaskedConv(filter_sz::NTuple{N,Integer},
                     channels=1=>1;
                     bias_first=false) where N
    padding = Int[]
    for f=filter_sz
        p = bias_first ? f : f-1
        push!(padding, p)
        push!(padding, 0)
    end
    padding = Tuple(padding)

    return Conv(filter_sz, channels, logℒ, pad=padding)
end

struct ConvCache{Tc,To,Td,Tb,Cd,Co,Cfo,Cdo,Te}
    σ::Tc
    out::To
    δℒℒ::Td

    b_rshpd::Tb
    cdims::Cd
    conv_out::Co
    conv_filter_out::Cfo
    conv_data_out::Cdo

    out2::Te
    valid::Bool
end

cache(l::Conv, in_T, in_sz) = begin
    if length(in_sz) == length(size(l.weight))-2
        in_sz   = (in_sz...,1,1)
        @info "did"
    end
    in_rspd = similar(l.weight, in_T, in_sz)

    b_rshpd = reshape(l.bias, map(_->1, l.stride)..., :, 1)
    cdims = DenseConvDims(in_sz, size(l.weight);
                          stride=l.stride, padding=l.pad,
                          dilation=l.dilation)
    xT = in_T
    wT = eltype(l.weight)
    conv_T = promote_type(xT, wT)
    conv_out = similar(l.weight, conv_T,
                    NNlib.output_size(cdims)...,
                    NNlib.channels_out(cdims), in_sz[end])

    out = similar(conv_out, promote_type(conv_T, eltype(l.bias)))

    δℒℒ = similar(out)
    #
    conv_filter_out = similar(δℒℒ, kernel_size(cdims)...,
                        channels_in(cdims), channels_out(cdims))
    conv_data_out   = similar(δℒℒ, NNlib.input_size(cdims)...,
                        channels_in(cdims),
                            size(δℒℒ)[end])

    ConvCache(in_rspd, out, δℒℒ,
              b_rshpd, cdims, conv_out, conv_filter_out, conv_data_out,
              0, false)
end

layer_out_type_size(l::Conv, in_T, in_sz) = begin
    if length(in_sz) == length(size(l.weight))-2
        in_sz   = (in_sz...,1,1)
        @info "did"
    end
    cdims = DenseConvDims(in_sz, size(l.weight);
                          stride=l.stride, padding=l.pad,
                          dilation=l.dilation)
    conv_T = promote_type(in_T, eltype(l.weight))
    return promote_type(conv_T, eltype(l.bias)), (NNlib.output_size(cdims)..., NNlib.channels_out(cdims)..., in_sz[end])
end

function (l::Conv)(c::ConvCache, x)
    σ    = l.σ
    b    = c.b_rshpd
    cdims = c.cdims

    x = isnothing(c.σ) ? x : copyto!(c.σ, x)

    conv!(c.conv_out, x, l.weight, cdims)
    c.conv_out .+= b
    c.out .= σ.(c.conv_out)
    return c.out
end

function backprop(∇, l::Conv, c::ConvCache, δℒ)
    δℒℒ = c.δℒℒ

    # Compute the actual sensitivity
    copyto!(δℒℒ, δℒ)
    if l.σ == identity
        δℒℒ .= c.conv_out
    elseif l.σ == logℒ
        δℒℒ .*= ∂logℒ.(c.conv_out)
    else
        throw("notknown $(l.σ)")
    end

    conv_filter_out = ∇conv_filter!(c.conv_filter_out, c.σ, δℒℒ, c.cdims)
    conv_data_out   = ∇conv_data!(c.conv_data_out, δℒℒ, l.weight, c.cdims)

    ∇.weight  .= c.conv_filter_out
    ∇.bias    .= sum(δℒℒ) # sum(δℒℒ, dims=1)
    return c.conv_data_out
end


text:
push!(LOAD_PATH, "../src/")

using Documenter, HSM

makedocs(
    modules = [HSM],
    sitename=  "HSM.jl",
    authors = "Andrew Lin",
    pages = [
        "Home" => "index.md",
        "reference.md"
    ]
)
deploydocs(
    repo = "github.com/AndrewWasHere/HSM.jl.git",
)


text:
### Monte Carlo routines for
### the 2D Ising model

import StatsBase: countmap
import StaticArrays: SVector

"""
    randomConfiguration(L)

Returns a LxL Matrix filled randomly with ±1.
"""
function randomConfiguration(L::Int)
    return rand(Int8[-1,1], L,L)
end

"""
    frustratedConfiguration(L)

Returns a configuration where all bonds are frustrated.
"""
function frustratedConfiguration(L::Int)
    return Int8[(-1)^(i+j) for i=1:L,j=1:L]
end

"""
    H(state, J=1.,h=0.)

Energy of Ising configuration for a given coupling `J` and external field `h`.
"""
function H(state::Matrix{Int8},J=1.,h=0.)
    L = size(state,1)
    s = h*(sum(state))

    @inbounds for i in 1:L, j in 1:L
        s += -J*state[i,j]*(state[i%L+1,j]+state[i,j%L+1])
    end
    return s
end

"""
    dH(state,pos,J=1.,h=0.)

Energy difference between the current state `state`
and one spin flipped at position `(i,j)` given coupling and field strength.
"""
function dH(state::Matrix{Int8},i,j,J=1.,h=0.)
    L = size(state,1)
    @inbounds return 2*J*state[i,j]*( state[i%L+1,j]+state[i,j%L+1]+
        state[i==1 ? L : i-1,j]+state[i,j==1 ? L : j-1] ) + 2*h*state[i,j]
end

"""
    m(state)

Magnetization per spin
"""
m(state) = Float64(sum(state)/length(state))

@inline neighbors(i,j,L) = CartesianIndex.([(i%L+1,j),(i,j%L+1),(i==1 ? L : i-1,j),(i,j==1 ? L : j-1)])
@inline function neighbors(i,L)
    r = (i-1)%L
    c = div(i-1, L)
    # @info r,c
    SVector{4, Int32}(1 + (r+1)%L + c*L, 1 + ((r-1)>=0 ? r-1 : L-1) + c*L, 1 + r + L*((c+1)%L), 1 + r + L*((c-1)>=0 ? c-1 : L-1))
end
@inline function neighbors_leftup(i,L)
    r = (i-1)%L
    c = div(i-1, L)
    # @info r,c
    SVector{2, Int32}(1 + ((r-1)>=0 ? r-1 : L-1) + c*L, 1 + r + L*((c-1)>=0 ? c-1 : L-1))
end

"""
    cluster_sizes(state)

A cluster is a connected region of aligned spins. Retrieve the sizes of all of them.
"""
function cluster_sizes(state)
    # Check for each element if they belong to the same cluster
    # as their left or upper neighbor.
    cluster = zeros(Int32, size(state))
    max_cluster = 0
    pointer = 1
    new_cluster = true
    while pointer <= length(state)
        nn = neighbors_leftup(pointer, size(state, 1))
        new_cluster = true
        for n in nn
            if cluster[n]!=0 && state[pointer]==state[n]
                cluster[pointer] = cluster[n]
                new_cluster = false
            end
        end
        if new_cluster
            cluster[pointer] = max_cluster += 1
        end
        pointer += 1
    end

    return countmap(reshape(cluster, length(cluster))) |> values |> collect
end

"""
    metropolis_step!(state,beta,h)

Perform one step of the Metropolis algorithm.
`state` is mutated.
"""
function metropolis_step!(state::Matrix{Int8},beta,h)
    i = rand(1:size(state,1))
    j = rand(1:size(state,2))
    dh = dH(state,i,j,1.,h)
    if dh <= 0 || rand()<exp(-beta*dh)
        state[i,j] *= -1
    end
    return nothing
end

"""
    metropolis_sweep!(state,n,beta,h)

Perform n Metropolis steps.
"""
function metropolis_sweep!(state,n,beta,h)
    for _ in 1:n
        metropolis_step!(state,beta,h)
    end
end

"""
    init(L,beta,h,sweep)

Setup a configuration with `LxL` spins and perform an intial thermal sweep of
`sweep` timesteps.
"""
function init(L,beta,h,sweep)
    state = fill(Int8(-1), L, L)
    # Initial sweep to get into the steady state
    metropolis_sweep!(state,sweep,beta,h)
    return state
end

"""
    run_metropolis(L, beta, h;Tmax=1, sweep=0, sample_interval=1)

Sets up a random state and runs the Metropolis algorithm for a given set of parameters.
Samples in defined intervals along the Markov-Chain.

An initial thermal sweep to go to equilibrium may be specified.

Returns an array of averaged observables: [E, E^2, m, m^2, m^4] with total
energy E and magnetisation per spin m.

# Arguments
- L::Integer:            Linear system size
- beta::Float:            Inverse temperature
- h::Float:               External field
- Tmax::Integer:          Number of steps
- sweep::Integer:         Length of the initial sweep
- sample_interval::Integer: sample interval

# Example
```
julia> run_metropolis(50, 0., 0.;Tmax=50*10^3*50^2,sample_interval=10*50^2,sweep=10^3*50^2)
5-element Array{Float64,1}:
  -39.9664      # <E>
 1753.5168      # <E^2>
    0.010382079999999908 # <m>
    0.00011572633600000021 # <m^2>
    1.6924325969920386e-8 # <m^4>
```
"""
function run_metropolis(L::Int, beta,h;Tmax::Int=1,sweep::Int=0,sample_interval::Int=1)
    ## Initialise a random state
    state = init(L,beta,h,sweep)
    return _run_metropolis!(state,beta,h,Tmax=Tmax,sample_interval=sample_interval)
end

"""
    _run_metropolis!(state::Matrix{Int8},beta,h;Tmax::Int=1,sample_interval::Int=1)

Analogous to `run_metropolis` but mutates an existing state. Called by `run_metropolis`.
"""
function _run_metropolis!(state::Matrix{Int8},beta,h;Tmax::Int=1,sample_interval::Int=1)

    ## Define a matrix to record the observables.
    ## One entry for each observable, e.g E and m.
    ## Preallocating the matrix gives much better performance than
    ## constructing it on the fly.
    observables = zeros(Float64, 5)

    k = 0 #counts the number of samples

    t = 0 #simulation steps
    e = 0.
    mag = 0.

    @inbounds begin
        while(t<Tmax)
            ## Take the defined no. of steps before
            ## recording a measurement
            metropolis_sweep!(state,sample_interval,beta,h)

            ## Record observables
            e = H(state,1.,h)
            mag = m(state)
            observables[1] += e
            observables[2] += e^2
            observables[3] += mag
            observables[4] += mag^2
            observables[5] += mag^4

            ## increment counters
            k+=1
            t+=sample_interval
        end
    end
    ## Return sample means
    return observables/k
end

"""
    metropolis_timeseries(L, β, Tmax; sample_interval=L^2, sweep=1000)

Follows along a Markov chain in time.
Initialises a random configuration, sweeps it for `sweep*L^2` timsteps, and records
magnetetisation and energy every `sample_interval` until `Tmax`.

Returns a tuple of these time series.
"""
function metropolis_timeseries(L::Int, beta, Tmax; sample_interval=L^2, sweep=1000)
    state = init(L,beta,0.,sweep*L^2)
    #time series magnetization
    tsm = Vector{Float64}(undef, div(Tmax,sample_interval))
    #time series energy
    tse = Vector{Float64}(undef, div(Tmax,sample_interval))

    t=0
    k=1
    mag0 = m(state)
    while(t<Tmax)
        metropolis_sweep!(state,sample_interval,beta,0.)
        tsm[k] = m(state)
        tse[k] = H(state)
        t+=sample_interval
        k+=1
    end
    return tsm, tse
end

### Wolff algorithm ###
### --------------- ###

"""
    cluster!(cluster_state, state, i,j, p)

Builds a cluster around position `(i,j)` with acceptance rate `p`.

`cluster_state` is a boolean matrix of the same size as the system. It is set to `true` whenever
a site belongs to the cluster.
"""
function cluster!(cluster_state, state, i,j, p)
    @inbounds begin L = size(state, 1)
        s = state[i,j]
        cluster_state[i,j] = true
        for neighbor in neighbors(i,j,L)
            if state[neighbor] == s && !cluster_state[neighbor] && rand()<p
                cluster_state[neighbor] = true
                cluster!(cluster_state, state, Tuple(neighbor)..., p)
            end
        end
    end
end

"""
    wolff_step!(state, cluster_state, beta)

Build a cluster and flip it.
"""
function wolff_step!(state, cluster_state, beta, h)
    @inbounds begin
        L = size(state,1)
        i,j = rand(1:L,2)
        fill!(cluster_state, false)
        cluster!(cluster_state, state, i,j, 1.0-exp(-2.0*beta))
        state[cluster_state] *= -1
    end
    return nothing
end

function wolff_sweep!(state, cluster, n, args...)
    for _ in 1:n
        wolff_step!(state, cluster, args...)
    end
end

"""
    _run_wolff!(state,cluster,beta,h;Tmax=1,sample_interval=1)

Run the Wolff cluster algorithm on a given state. You need to provide a cluster matrix, e.g
`similar(state,Bool)`.

See also: [`_run_metropolis`]
"""
function _run_wolff!(state::Matrix{Int8},cluster,beta,h;Tmax::Int=1,sample_interval::Int=1)

    ## Define a matrix in which to record the observables.
    ## One row for each observable, e.g E and m.
    ## Preallocating the matrix gives much better performance than
    ## constructing it on the fly.
    observables = zeros(Float64,5)

    k = 0 #counts the number of samples

    # pick the middle lattice site for the spin correlation

    t=0 #number of simulation steps
    r=0 #repeat counter
    e=0.
    mag=0.
    @inbounds begin
        while(t<Tmax)
            ## Do the defined no. of steps before
            ## taking a measurement
            for _ in 1:sample_interval
                wolff_step!(state, cluster, beta, h)
            end
            ## Record observables
            e = H(state,1.,0.)
            mag = m(state)|>abs
            observables[1] += e
            observables[2] += e^2
            observables[3] += mag
            observables[4] += mag^2
            observables[5] += mag^4

            ## increment counters
            k+=1
            t+=sample_interval
        end
    end
    ## Return statistics about the observables
    return observables/k
end

"""
    run_wolff(L, beta,h;Tmax=1,sweep=0,sample_interval=1)

See also: [`run_metropolis`]
"""
function run_wolff(L::Int, beta,h;Tmax::Int=1,sweep::Int=0,sample_interval::Int=1)
    ## Initialise a random state
    function init(L,beta)
        state = randomConfiguration(L)
        cluster = zeros(Bool,L,L)
        # Initial sweep to get into the steady state
        for _ in 1:sweep
            wolff_step!(state, cluster, beta, h)
        end
        return state,cluster
    end

    state,cluster = init(L,beta)
    return _run_wolff!(state,cluster,beta,0.,Tmax=Tmax,sample_interval=sample_interval)
end

max_stars_repo_name: skleinbo/Ising

max_issues_repo_name: Yaren-Github/Ising

max_forks_repo_name: Yaren-Github/Ising


text:
abstract type MatchingMethod end
struct ΛTAME_GramMatching <: MatchingMethod end
struct ΛTAME_rankOneMatching <: MatchingMethod end

struct ΛTAME_Return{T} <: returnType
	matchScore::Union{Int,Vector{Int}}
	motifCounts::NTuple{2,Union{Int,Vector{Int}}}
	matching::Union{Dict{Int,Int},Vector{Int}}
	rank1Matchingindices::Union{Nothing,Tuple{Int,Int}}
	embedding::Tuple{Matrix{T},Matrix{T}}
	profile::Union{Nothing,Dict{String,Array{Float64,1}}}
end

struct ΛTAME_MultiMotif_Return{T} <: returnType
#    motifDistributions::NTuple(2,Vector{T}) # ttv({A,B},1,modes=k-1)
	aggregateMotifScore::Int
	matchScore::Union{Int,Vector{Int}}
	motifCounts::NTuple{2,Union{Int,Vector{Int}}}
	matching::Union{Dict{Int,Int},Vector{Int}}
	rank1Matchingindices::Union{Nothing,Tuple{Int,Int}}
	embedding::Tuple{Matrix{T},Matrix{T}}
	profile::Union{Nothing,Dict{String,Array{Float64,1}}}
end

#=------------------------------------------------------------------------------
              Routines for searching over alpha/beta parameters
------------------------------------------------------------------------------=#

function ΛTAME_param_search(A,B;#duck type tensor inputs
                            iter::Int = 15,tol::Float64=1e-6,
							alphas::Array{F,1}=[.5,1.0],
							betas::Array{F,1} =[1000.0,100.0,10.0,1.0,0.0,0.1,0.01,0.001],
							matchingMethod::MatchingMethod=ΛTAME_rankOneMatching(),kwargs...) where {F <: AbstractFloat}

	motifA = size(A.indices,2)
	motifB = size(B.indices,2)

	max_triangle_match = min(motifA ,motifB)
    #max_triangle_match = min(size(A.indices,1),size(B.indices,1))
    total_triangles = size(A.indices,1) + size(B.indices,1)

    U = Array{Float64,2}(undef,A.n,iter)
    V = Array{Float64,2}(undef,B.n,iter)
	best_U = Array{Float64,2}(undef,A.n,iter)
	best_V = Array{Float64,2}(undef,B.n,iter)

    best_TAME_PP_tris = -1
    best_i  = -1
	best_j = -1
	best_matching = Dict{Int,Int}()

    for α in alphas
        for beta in betas

			U,V = ΛTAME(A,B,beta,iter,tol,α)
			if typeof(matchingMethod) === ΛTAME_rankOneMatching
				search_tris, i, j, matching = search_Krylov_space(A,B,U,V)

			elseif typeof(matchingMethod) === ΛTAME_GramMatching
				search_tris, _,matching = TAME_score(A,B,U*V';kwargs...)
				i = -1
				j = -1
			else
				throw(TypeError("typeof MatchingMethod must be either ΛTAME_rankOneMatching or ΛTAME_GramMatching"))
			end
			println("α:$(α) -- β:$(beta) finished -- tri_match:$search_tris -- max_tris $(max_triangle_match) -- best tri_match:$best_TAME_PP_tris")
            if search_tris > best_TAME_PP_tris
                best_TAME_PP_tris = search_tris
                best_i = i
				best_j = j
				best_U = U
				best_V = V
				best_matching = matching
            end
        end
    end
	println("best i:$best_i -- best j:$best_j")
	if typeof(matchingMethod) === ΛTAME_rankOneMatching
		return ΛTAME_Return(best_TAME_PP_tris,(motifA,motifB),
						    best_matching,(best_i, best_j),(best_U, best_V),nothing)
	else
		return ΛTAME_Return(best_TAME_PP_tris,(motifA,motifB),
							best_matching,nothing,(best_U, best_V),nothing)
	end

end

function ΛTAME_param_search(A_tensors::Array{SymTensorUnweighted{S},1},B_tensors::Array{SymTensorUnweighted{S},1};
							iter::Int = 15,tol::Float64=1e-6,
							alphas::Array{F,1}=[.5,1.0],
							betas::Array{F,1} =[1000.0,100.0,10.0,1.0,0.0,0.1,0.01,0.001],
							matchingMethod::MatchingMethod=ΛTAME_rankOneMatching(),kwargs...) where {F <: AbstractFloat, S <: Motif}

	motifCountA = [size(A.indices,2) for A in A_tensors]
	motifCountB = [size(B.indices,2) for B in B_tensors]
	max_motif_match = [min(A_motifCount,B_motifCount) for (A_motifCount,B_motifCount) in zip(motifCountA,motifCountB)]
	total_triangles = [size(A.indices,2) + size(B.indices,2) for (A,B) in zip(A_tensors,B_tensors)]
	#total_triangles = size(A.indices,1) + size(B.indices,1)

	m = maximum([A.n for A in A_tensors])
	n = maximum([B.n for B in B_tensors])

	U = Array{Float64,2}(undef,m,iter)
	V = Array{Float64,2}(undef,n,iter)
	best_U = Array{Float64,2}(undef,m,iter)
	best_V = Array{Float64,2}(undef,n,iter)

	best_matching_score = -1
	best_i  = -1
	best_j = -1
	best_matched_motifs::Array{Int,1} = []
	best_matching = Dict{Int,Int}()

	for α in alphas
		for beta in betas

			U,V = ΛTAME(A_tensors,B_tensors,beta,iter,tol,α)

			if typeof(matchingMethod) === ΛTAME_rankOneMatching
				matching_score, matched_motifs, i, j, matching = search_Krylov_space(A_tensors,B_tensors,U,V)
			elseif typeof(matchingMethod) === ΛTAME_GramMatching
				matching_score, matched_motifs, matching = TAME_score(A_tensors,B_tensors,U*V';kwargs...)
				i = -1
				j = -1
			else
				throw(TypeError("typeof MatchingMethod must be either ΛTAME_rankOneMatching or ΛTAME_GramMatching"))
			end

			if matching_score > best_matching_score
				best_matching_score = matching_score
				best_matched_motifs = matched_motifs
				best_i = i
				best_j = j
				best_U = U
				best_V = V
				best_matching = matching
			end
			println("α:$(α) -- β:$(beta) finished -- motif_match:$matched_motifs -- max_motifs:$max_motif_match -- matching_score:$matching_score -- best score:$best_matching_score")

		end
	end

	println("best i:$best_i -- best j:$best_j")
	if typeof(matchingMethod) === ΛTAME_rankOneMatching
		return ΛTAME_MultiMotif_Return(best_matching_score,best_matched_motifs, (motifCountA,motifCountB),
					best_matching,(best_i, best_j),(best_U, best_V),nothing)
	else
		return ΛTAME_MultiMotif_Return(best_matching_score,best_matched_motifs, (motifCountA,motifCountB),
					best_matching,nothing,(best_U, best_V),nothing)
	end

end

function ΛTAME_param_search_profiled(A,B;
	                                 iter::Int = 15,tol::Float64=1e-6, alphas::Array{F,1}=[.5,1.0],
									 betas::Array{F,1} =[1000.0,100.0,10.0,1.0,0.0,0.1,0.01,0.001],
									 matchingMethod::MatchingMethod=ΛTAME_rankOneMatching(),kwargs...) where {F <: AbstractFloat}

	motifA = size(A.indices,2)
	motifB = size(B.indices,2)

	max_motif_match = min(motifA ,motifB)
    best_matched_motifs = -1
    best_i  = -1
	best_j = -1
	best_matching = Dict{Int,Int}()

	m = A.n
	n = B.n

	results = Dict(
		"TAME_timings" => Array{Float64,1}(undef,length(alphas)*length(betas)),
		"Matching Timings"=> Array{Float64,1}(undef,length(alphas)*length(betas)),
		"Scoring Timings"=> Array{Float64,1}(undef,length(alphas)*length(betas)),
	)
	exp_index = 1

    U = Array{Float64,2}(undef,m,iter)
    V = Array{Float64,2}(undef,n,iter)
	best_U = Array{Float64,2}(undef,m,iter)
	best_V = Array{Float64,2}(undef,n,iter)

    for α in alphas
        for beta in betas

			((U,V),runtime) = @timed ΛTAME(A,B,beta,iter,tol,α)
			results["TAME_timings"][exp_index] = runtime

			if typeof(matchingMethod) === ΛTAME_rankOneMatching
				((matched_motifs, i, j, matching,scoring_time),runtime) = @timed search_Krylov_space(A,B,U,V;returnScoringTimings=returnTimings(),kwargs...)
				results["Matching Timings"][exp_index] = runtime - scoring_time
				results["Scoring Timings"][exp_index] = scoring_time
			elseif typeof(matchingMethod) === ΛTAME_GramMatching
				(matched_motifs, gaped_motifs, matching, matching_time, scoring_time)= TAME_score(A,B,U*V';return_timings=returnTimings(),kwargs...)
				i = -1
				j = -1
				best_U = copy(U)
				best_V = copy(V)
				results["Matching Timings"][exp_index] = matching_time
				results["Scoring Timings"][exp_index] = scoring_time
			else
				throw(TypeError("typeof MatchingMethod must be either ΛTAME_rankOneMatching or ΛTAME_GramMatching"))
			end

			exp_index +=1

			if matched_motifs > best_matched_motifs
				best_matched_motifs = matched_motifs
				best_i = i
				best_j = j
				best_U = U
				best_V = V
				best_matching = matching
			end

			println("α:$(α) -- β:$(beta) finished -- motif_match:$matched_motifs -- max_motif $max_motif_match -- best motif_match: $best_matched_motifs")
        end
    end

	println("best i:$best_i -- best j:$best_j")
	if typeof(matchingMethod) === ΛTAME_rankOneMatching
		return ΛTAME_Return(best_matched_motifs,(motifA ,motifB),
						    best_matching,(best_i, best_j),(best_U, best_V),results)
	else
		return ΛTAME_Return(best_matched_motifs,(motifA ,motifB),
							best_matching,nothing,(best_U, best_V),results)
	end

end

function ΛTAME_param_search_profiled(A_tensors::Array{SymTensorUnweighted{S},1},B_tensors::Array{SymTensorUnweighted{S},1};
							iter::Int = 15,tol::Float64=1e-6,
							alphas::Array{F,1}=[.5,1.0],
							betas::Array{F,1} =[1000.0,100.0,10.0,1.0,0.0,0.1,0.01,0.001],
							matchingMethod::MatchingMethod=ΛTAME_rankOneMatching(),kwargs...) where {F <: AbstractFloat, S <: Motif}

	motifCountA = [size(A.indices,2) for A in A_tensors]
	motifCountB = [size(B.indices,2) for B in B_tensors]
	max_motif_match = [min(A_motifCount,B_motifCount) for (A_motifCount,B_motifCount) in zip(motifCountA,motifCountB)]

	m = maximum([A.n for A in A_tensors])
	n = maximum([B.n for B in B_tensors])

	U = Array{Float64,2}(undef,m,iter)
	V = Array{Float64,2}(undef,n,iter)
	best_U = Array{Float64,2}(undef,m,iter)
	best_V = Array{Float64,2}(undef,n,iter)

	profiling = Dict(
		"TAME_timings" => Array{Float64,1}(undef,length(alphas)*length(betas)),
		"Matching Timings"=> Array{Float64,1}(undef,length(alphas)*length(betas)),
		"Scoring Timings"=> Array{Float64,1}(undef,length(alphas)*length(betas)),
	)
	exp_index = 1

	best_matching_score = -1
	best_i  = -1
	best_j = -1
	best_matched_motifs::Array{Int,1} = []
	best_matching = Dict{Int,Int}()

	for α in alphas
		for beta in betas

			(U,V),ΛTAME_runtime = @timed ΛTAME(A_tensors,B_tensors,beta,iter,tol,α)
			profiling["TAME_timings"][exp_index] = ΛTAME_runtime

			if typeof(matchingMethod) === ΛTAME_rankOneMatching
				((matching_score, matched_motifs, i, j, matching,scoring_time), full_runtime) = @timed search_Krylov_space(A_tensors,B_tensors,U,V;returnScoringTimings=returnTimings())
				profiling["Matching Timings"][exp_index] = full_runtime - scoring_time
				profiling["Scoring Timings"][exp_index] = scoring_time
			elseif typeof(matchingMethod) === ΛTAME_GramMatching
				((matching_score, matched_motifs, matching, matching_time, scoring_time),matching_time) = @timed TAME_score(A_tensors,B_tensors,U*V';return_timings=returnTimings(),kwargs...)
				i = -1
				j = -1
				profiling["Matching Timings"][exp_index] = matching_time
				profiling["Scoring Timings"][exp_index] = scoring_time
			else
				throw(TypeError("typeof MatchingMethod must be either ΛTAME_rankOneMatching or ΛTAME_GramMatching"))
			end

			exp_index +=1

			if matching_score > best_matching_score
				best_matching_score = matching_score
				best_matched_motifs = matched_motifs
				best_i = i
				best_j = j
				best_U = U
				best_V = V
				best_matching = matching
			end
			println("α:$(α) -- β:$(beta) finished -- motif_match:$matched_motifs -- max_motifs:$max_motif_match -- matching_score:$matching_score -- best score:$best_matching_score")

		end
	end

	println("best i:$best_i -- best j:$best_j")
	if typeof(matchingMethod) === ΛTAME_rankOneMatching
		return ΛTAME_MultiMotif_Return(best_matching_score,best_matched_motifs, (motifCountA,motifCountB),
					best_matching,(best_i, best_j),(best_U, best_V),profiling)
	else
		return ΛTAME_MultiMotif_Return(best_matching_score,best_matched_motifs, (motifCountA,motifCountB),
					best_matching,nothing,(best_U, best_V),profiling)
	end

end

#=------------------------------------------------------------------------------
             		    Spectral Relaxation Routines
------------------------------------------------------------------------------=#

"""------------------------------------------------------------------------------
  The LambdaTAME method. Method starts with a uniform prior and computes the
  power method for each of the tensors passed it. These iterates are stored and
  the best rank one matching is picked out of all the quadratic pairs of vectors
  computed. This portion of the algorith computes the contractions, the full
  procedure can be found in the associated '_param_search' algorithms.

  Inputs
  ------
  * A,B - (ThirdOrderSymTensor):
	The tensors to align against one another.
  * β - (Float):
	The shift to use on the iterations.
  * α -(float):
	The mixing parameter for combining in the starting iterates.
  * 'max_iter' - (Int):
	The maximum number of iterations to run.
  * tol - (Float):
	The tolerence to solve the algorithm to. Computes the tolerance by measuring
	the absolute value of the difference between the computed eigenvalues.
  * 'update_user' - (Int):
	Specifies how frequently output messages should be printed to the user.
	Default is -1 which means no output messages are printed, else if the value
	is k, then the kth iterate will print out an update.
  Output
  ------
  * U, V - (Array{Float,2})
    Returns the computed contractions started with a uniform prior.
------------------------------------------------------------------------------"""
function ΛTAME(A::ThirdOrderSymTensor, B::ThirdOrderSymTensor, β::Float64,
               max_iter::Int,tol::Float64,α::Float64;update_user=-1)

    U = zeros(A.n,max_iter+1)
    V = zeros(B.n,max_iter+1) #store initial in first column

    U[:,1] = ones(A.n)
    U[:,1] /=norm(U[:,1])

    V[:,1] = ones(B.n)
    V[:,1] /=norm(V[:,1])

    sqrt_β = β^(.5)

    lambda = Inf
    i = 1

    while true

        U[:,i+1] = tensor_vector_contraction(A,U[:,i])
        V[:,i+1] = tensor_vector_contraction(B,V[:,i])

        lambda_A = (U[:,i+1]'*U[:,i])
        lambda_B = (V[:,i+1]'*V[:,i])
        new_lambda = lambda_A*lambda_B

        if β != 0.0
            U[:,i+1] .+= sqrt_β*U[:,i+1]
            V[:,i+1] .+= sqrt_β*V[:,i+1]
        end

        if α != 1.0
            U[:,i+1] = α*U[:,i+1] + (1 -α)*U[:,1]
            V[:,i+1] = α*V[:,i+1] + (1 -α)*V[:,1]
        end

        U[:,i+1] ./= norm(U[:,i+1])
        V[:,i+1] ./= norm(V[:,i+1])

		if update_user != -1 && i % update_user == 0
            println("iteration $(i)    λ_A: $(lambda_A) -- λ_B: $(lambda_B) -- newλ: $(new_lambda)")
		end

        if abs(new_lambda - lambda) < tol || i >= max_iter
            return U[:,1:i], V[:,1:i]
        else
            lambda = new_lambda
            i += 1
        end

    end

end

function ΛTAME(A::SymTensorUnweighted{S}, B::SymTensorUnweighted{S}, β::Float64,
	max_iter::Int,tol::Float64,α::Float64;update_user=-1) where {S <: Motif}

	U = zeros(A.n,max_iter+1)
	V = zeros(B.n,max_iter+1) #store initial in first column

	U[:,1] = ones(A.n)
	U[:,1] /=norm(U[:,1])

	V[:,1] = ones(B.n)
	V[:,1] /=norm(V[:,1])

	sqrt_β = β^(.5)

	lambda = Inf
	i = 1

	A_buf = zeros(A.n)
	B_buf = zeros(B.n)

	while true

		DistributedTensorConstruction.contraction!(A,U[:,i],A_buf)
		DistributedTensorConstruction.contraction!(B,V[:,i],B_buf)

		U[:,i+1] .= A_buf
		V[:,i+1] .= B_buf

		lambda_A = (U[:,i+1]'*U[:,i])
		lambda_B = (V[:,i+1]'*V[:,i])
		new_lambda = lambda_A*lambda_B

		if β != 0.0
			U[:,i+1] .+= sqrt_β*U[:,i+1]
			V[:,i+1] .+= sqrt_β*V[:,i+1]
		end

		if α != 1.0
			U[:,i+1] = α*U[:,i+1] + (1 -α)*U[:,1]
			V[:,i+1] = α*V[:,i+1] + (1 -α)*V[:,1]
		end

		U[:,i+1] ./= norm(U[:,i+1])
		V[:,i+1] ./= norm(V[:,i+1])

		if update_user != -1 && i % update_user == 0
			println("iteration $(i)    λ_A: $(lambda_A) -- λ_B: $(lambda_B) -- newλ: $(new_lambda)")
		end

		if abs(new_lambda - lambda) < tol || i >= max_iter
			return U[:,1:i], V[:,1:i]
		else
			lambda = new_lambda
			A_buf .= 0.0
			B_buf .= 0.0
			i += 1
		end

	end

end

function ΛTAME(A_tensors::Array{SymTensorUnweighted{S},1}, B_tensors::Array{SymTensorUnweighted{S},1}, β::Float64,
	max_iter::Int,tol::Float64,α::Float64;update_user=-1) where {S <: Motif}

	@assert length(A_tensors) == length(B_tensors)
	for i =1:length(A_tensors)
		@assert A_tensors[i].order == B_tensors[i].order
	end

	m = maximum([tensor.n for tensor in A_tensors])
	n = maximum([tensor.n for tensor in B_tensors])

	U = zeros(m,max_iter+1)
	V = zeros(n,max_iter+1) #store initial in first column

	U[:,1] = ones(m)
	U[:,1] /=norm(U[:,1])

	V[:,1] = ones(n)
	V[:,1] /=norm(V[:,1])

	sqrt_β = β^(.5)

	lambda = Inf
	i = 1

	A_buf = zeros(m)
	B_buf = zeros(n)

	while true

		#DistributedTensorConstruction.embedded_contraction!(A_tensors,U[:,i],A_buf)
		#DistributedTensorConstruction.embedded_contraction!(B_tensors,V[:,i],B_buf)
		embedded_contraction!(A_tensors,U[:,i],A_buf)
		embedded_contraction!(B_tensors,V[:,i],B_buf)

		#contraction_divide_out!(A_tensors,U[:,i],A_buf)
		#contraction_divide_out!(B_tensors,V[:,i],B_buf)

		U[:,i+1] .= A_buf
		V[:,i+1] .= B_buf

		lambda_A = (U[:,i+1]'*U[:,i])
		lambda_B = (V[:,i+1]'*V[:,i])
		new_lambda = lambda_A*lambda_B

		if β != 0.0
			U[:,i+1] .+= sqrt_β*U[:,i+1]
			V[:,i+1] .+= sqrt_β*V[:,i+1]
		end

		if α != 1.0
			U[:,i+1] = α*U[:,i+1] + (1 -α)*U[:,1]
			V[:,i+1] = α*V[:,i+1] + (1 -α)*V[:,1]
		end

		U[:,i+1] ./= norm(U[:,i+1])
		V[:,i+1] ./= norm(V[:,i+1])

		if update_user != -1 && i % update_user == 0
			println("iteration $(i)    λ_A: $(lambda_A) -- λ_B: $(lambda_B) -- newλ: $(new_lambda)")
		end

		if abs(new_lambda - lambda) < tol || i >= max_iter
			return U[:,1:i], V[:,1:i]
		else
			lambda = new_lambda
			A_buf .= 0.0
			B_buf .= 0.0
			i += 1
		end

	end

end

max_stars_repo_name: chuckcol/LambdaTAME

max_issues_repo_name: chuckcol/LambdaTAME

max_forks_repo_name: chuckcol/LambdaTAME


text:
# Additional postprocessors for @planar and @plansor

# since temporaries were taken out in preprocessing by _decompose_planar_contractions, they
# are not identified by the parsing step of TensorOperations, so we have to manually fix this
# Step 1: we have to find the new name that TO.tensorify assigned to these temporaries
# since it parses `tmp[] := a[] * b[]` as `newtmp = similar...; tmp = contract!(.... , newtmp, ...)`
function _update_temporaries(ex, temporaries)
    if ex isa Expr && ex.head == :(=)
        lhs = ex.args[1]
        i = findfirst(==(lhs), temporaries)
        if i !== nothing
            rhs = ex.args[2]
            if rhs isa Expr && rhs.head == :call && rhs.args[1] == :add!
                newname = rhs.args[6]
                temporaries[i] = newname
            elseif rhs isa Expr && rhs.head == :call && rhs.args[1] == :contract!
                newname = rhs.args[8]
                temporaries[i] = newname
            elseif rhs isa Expr && rhs.head == :call && rhs.args[1] == :trace!
                newname = rhs.args[6]
                temporaries[i] = newname
            else
                @error "lhs = $lhs , rhs = $rhs"
            end
        end
    elseif ex isa Expr
        for a in ex.args
            _update_temporaries(a, temporaries)
        end
    end
    return ex
end

# Step 2: we find `newtmp = similar_from_...` and replace with `newtmp = cached_similar_from...`
function _annotate_temporaries(ex, temporaries)
    if ex isa Expr && ex.head == :(=)
        lhs = ex.args[1]
        i = findfirst(==(lhs), temporaries)
        if i !== nothing
            rhs = ex.args[2]
            if !(rhs isa Expr && rhs.head == :call && rhs.args[1] == :similar_from_indices)
                @error "lhs = $lhs , rhs = $rhs"
            end
            newrhs = Expr(:call, :cached_similar_from_indices,
                            QuoteNode(lhs), rhs.args[2:end]...)
            return Expr(:(=), lhs, newrhs)
        end
    elseif ex isa Expr
        return Expr(ex.head, [_annotate_temporaries(a, temporaries) for a in ex.args]...)
    end
    return ex
end

# add correct modules (`GlobalRef`) to various functions
const _TOFUNCTIONS = (:similar_from_indices, :cached_similar_from_indices,
                        :scalar, :IndexError)

function _add_modules(ex::Expr)
    if ex.head == :call && ex.args[1] in _TOFUNCTIONS
        return Expr(ex.head, GlobalRef(TensorOperations, ex.args[1]),
                        (_add_modules(ex.args[i]) for i in 2:length(ex.args))...)
    elseif ex.head == :call && ex.args[1] == :add!
        @assert ex.args[4] == :(:N)
        argind = [2,3,5,6,7,8]
        return Expr(ex.head, GlobalRef(TensorKit, Symbol(:planar_add!)),
                        (_add_modules(ex.args[i]) for i in argind)...)
    elseif ex.head == :call && ex.args[1] == :trace!
        @assert ex.args[4] == :(:N)
        argind = [2,3,5,6,7,8,9,10]
        return Expr(ex.head, GlobalRef(TensorKit, Symbol(:planar_trace!)),
                        (_add_modules(ex.args[i]) for i in argind)...)
    elseif ex.head == :call && ex.args[1] == :contract!
        @assert ex.args[4] == :(:N) && ex.args[6] == :(:N)
        argind = vcat([2,3,5], 7:length(ex.args))
        return Expr(ex.head, GlobalRef(TensorKit, Symbol(:planar_contract!)),
                        (_add_modules(ex.args[i]) for i in argind)...)
    else
        return Expr(ex.head, (_add_modules(e) for e in ex.args)...)
    end
end
_add_modules(ex) = ex


text:
targetSystem = Sys.WORD_SIZE;
println("Sys.WORD_SIZE indicates whether the target system is 32-bit or 64-bit. ['$targetSystem']");

integerType = typeof(1);
println("integerType is ['$integerType']");

#Unsigned integers are input and output using the 0x prefix and
#hexadecimal (base 16) digits 0-9a-f (the capitalized digits A-F also work for input).
unsignedValue = 0x1;
unsignedValuesType = typeof(unsignedValue);
println("unsignedValue(0x1) has value: ['$unsignedValue'] and is of type $unsignedValuesType");

#over-flow behavior
maxInt = typemax(Int64);
println("Max integer value: [$maxInt]");

overflowResult = maxInt + 1; #overflow by 1
println("Max integer after an overflow by 1: [$overflowResult]");

println("floating-point literals:");
println("------------------------");

value = 1.0;
println("1.0 -> $value");

value = 1.;
println("1. -> $value");

value = 0.5;
println("0.5 -> $value");

value = .5;
println(".5 -> $value");

value = -1.23;
println("-1.23 -> $value");

value = 1e10;
println("1e10 -> $value");

value = 2.5e-4;
println("2.5e-4 -> $value");

println(typeof(value));

println();
println("Hexadecimal floating-point literals, only have Float64 values, with p preceding the base-2 exponent");

hexValue = 0x1p0;
println("0x1p0 -> $hexValue");

hexValue = 0x1.8p3;
println("0x1.8p3 -> $hexValue");

hexValue = 0x.4p-1;
println("0x.4p-1 -> $hexValue");

println();
println("Underscores can be used as a digit seperator (10,000 would be 10_000)");
seperatedDigit = 10_000;
println("10_0000 -> $seperatedDigit");

println();
println("Floating-Point zero");

areZerosEqual = 0.0 == -0.0;
println("0.0 == -0.0 -> $areZerosEqual");
println("0.0 -> " * bitstring(0.0));
println("-0.0 -> " * bitstring(-0.0));

println();
println("Numeric Literal Coefficients");
x = 3;

answer = 2x^2 - 3x + 1;
println("2x^2 - 3x + 1 = " * string(answer));

answer = 1.5x^2 - .5x + 1;
println("1.5x^2 - .5x + 1 = " * string(answer));

answer = 2(x-1)^2 - 3(x -1) + 1;
println("2(x-1)^2 - 3(x -1) + 1 = " * string(answer));

println();
println("Literal zero and one");
println("zero(Float32) -> " * string(zero(Float32)));
println("zero(1.0) -> " * string(zero(1.0)));
println("one(Int64) -> " * string(one(Int64)));
println("one(BigFloat) -> " * string(one(BigFloat)));

max_stars_repo_name: CubicleJockey/julia-playground

max_issues_repo_name: CubicleJockey/julia-playground

max_forks_repo_name: CubicleJockey/julia-playground


text:
using GAIO

# logistic map
μ = 4.0
f(x) = μ.*x.*(1.0.-x)

center, radius = 0.5, 0.5
P = BoxPartition(Box(center, radius), 256)
F = BoxMap(f, P; no_of_points=400)

T = TransferOperator(F, P[:])
(λ, ev) = eigs(T)

plot(abs ∘ ev[1])


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

const EMPTY_VECTOR = Vector{Any}()

mutable struct InferenceResult
    linfo::MethodInstance
    args::Vector{Any}
    result # ::Type, or InferenceState if WIP
    src::Union{CodeInfo, Nothing} # if inferred copy is available
    function InferenceResult(linfo::MethodInstance)
        if isdefined(linfo, :inferred_const)
            result = Const(linfo.inferred_const)
        else
            result = linfo.rettype
        end
        return new(linfo, EMPTY_VECTOR, result, nothing)
    end
end

function get_argtypes(result::InferenceResult)
    result.args === EMPTY_VECTOR || return result.args # already cached
    linfo = result.linfo
    toplevel = !isa(linfo.def, Method)
    atypes::SimpleVector = unwrap_unionall(linfo.specTypes).parameters
    nargs::Int = toplevel ? 0 : linfo.def.nargs
    args = Vector{Any}(undef, nargs)
    if !toplevel && linfo.def.isva
        if linfo.specTypes == Tuple
            if nargs > 1
                atypes = svec(Any[ Any for i = 1:(nargs - 1) ]..., Tuple.parameters[1])
            end
            vararg_type = Tuple
        else
            vararg_type = rewrap(tupleparam_tail(atypes, nargs), linfo.specTypes)
        end
        args[nargs] = vararg_type
        nargs -= 1
    end
    laty = length(atypes)
    if laty > 0
        if laty > nargs
            laty = nargs
        end
        local lastatype
        atail = laty
        for i = 1:laty
            atyp = atypes[i]
            if i == laty && isvarargtype(atyp)
                atyp = unwrapva(atyp)
                atail -= 1
            end
            while isa(atyp, TypeVar)
                atyp = atyp.ub
            end
            if isa(atyp, DataType) && isdefined(atyp, :instance)
                # replace singleton types with their equivalent Const object
                atyp = Const(atyp.instance)
            elseif isconstType(atyp)
                atyp = Const(atyp.parameters[1])
            else
                atyp = rewrap_unionall(atyp, linfo.specTypes)
            end
            i == laty && (lastatype = atyp)
            args[i] = atyp
        end
        for i = (atail + 1):nargs
            args[i] = lastatype
        end
    else
        @assert nargs == 0 "invalid specialization of method" # wrong number of arguments
    end
    result.args = args
    return args
end

function cache_lookup(code::MethodInstance, argtypes::Vector{Any}, cache::Vector{InferenceResult})
    method = code.def::Method
    nargs::Int = method.nargs
    method.isva && (nargs -= 1)
    for cache_code in cache
        # try to search cache first
        cache_args = cache_code.args
        if cache_code.linfo === code && length(cache_args) >= nargs
            cache_match = true
            # verify that the trailing args (va) aren't Const
            for i in (nargs + 1):length(cache_args)
                if isa(cache_args[i], Const)
                    cache_match = false
                    break
                end
            end
            cache_match || continue
            for i in 1:nargs
                a = argtypes[i]
                ca = cache_args[i]
                # verify that all Const argument types match between the call and cache
                if (isa(a, Const) || isa(ca, Const)) && !(a === ca)
                    cache_match = false
                    break
                end
            end
            cache_match || continue
            return cache_code
        end
    end
    return nothing
end

max_stars_repo_name: djsegal/julia-fork

max_issues_repo_name: djsegal/julia-fork

max_forks_repo_name: djsegal/julia-fork


text:
module CytofImg

import Cytof5
using Plots
pyplot()
using StatPlots
using Distributions # std, mean, quantile
using RCall # TODO: GET RID OF THIS DEPENDENCY!

startupMsg = """
To scale font size globally for plots (for presentations mostly), do this:
```
julia> Plots.scalefontsizes(2) # scales all texts by factor of 2

# To reset
julia> Plots.reset_defaults()
```
"""
println(startupMsg)

function getPosterior(sym::Symbol, monitor)
  return [ m[sym] for m in monitor]
end

"""
turn on legend with `legend=true`

Note that xlabel / ylabel can be supplied by either:

```
julia> xlabel!("my xlab")
# OR
julia> plotZ(Z, xlabel="my xlab")
```

"""
function plotZ(Z::Matrix{T}; kw...) where {T <: Number}
  # FIXME: doesn't work for Z having either 1 row or 1 column!
  J, K = size(Z)
  if J < 2 || K < 2
    println("heatmap is not supported for single row / column matrices")
    return
  end

  img = heatmap(Z, c=:Greys, legend=:none, border=true, bordercolor=:lightgrey; kw...)
  if J > 1
    hline!((2:J) .- .5, c=:lightgrey)
  end
  if K > 1
    vline!((2:K) .- .5, c=:lightgrey)
  end

  return img
end

function plotY_heatmap(y::Matrix{T}; clim=(-5, 5), c=:pu_or) where {T <: Number}
  img = heatmap(y, c=c, clim=clim, background_color_inside=:black, bordercolor=:transparent);
  return img
end

function plotYZ(y::Matrix{V}, Z::Matrix{T}; ycolor=:balance, clim=(-5, 5), heights=[.7, .3], kw...) where {T <: Number, V <: Number}
  # TODO: change proportions
  heatY = plotY_heatmap(y, c=ycolor, clim=clim)
  heatZ = plotZ(Matrix(Z'))
  img = plot(heatY, heatZ, layout=grid(2, 1, heights=heights); kw...);
  return img
end

function annotateHeatmap(mat::Matrix{T}; digits=2, textsize=6) where T
  rows, cols = size(mat)
  for c in 1:cols
    for r in 1:rows
      annotate!([(float(c), float(r),
                  text("$(round(mat[r,c], digits=digits))", textsize))])
    end
  end
end

function plotPost(x::Vector{T}; a::Float64=0.05, q_digits::Int=3, sd_digits::Int=3, add=false,
                  trace::Bool=true, useDensity::Bool=true, showAccRate::Bool=true,
                  parent=1, offset=1, traceFont=font(16), kw...) where {T <: Number}
  denColor = (0, .99, :steelblue)
  if useDensity
    if add
      img = StatPlots.density!(x, label="", grid=false, yaxis=false,
                               linecolor=:transparent,
                               fill=denColor, bordercolor=:transparent; kw...)
    else
      img = StatPlots.density(x, label="", grid=false, bordercolor=:transparent, yaxis=false,
                              fill=denColor, linecolor=:transparent; kw...)
    end
  else
    if add
      img = histogram!(x, label="", normed=true, grid=false,
                       bordercolor=:white, yaxis=false,
                       linecolor=:transparent, c=:steelblue; kw...)
    else
      img = histogram(x, label="", normed=true, grid=false,
                      bordercolor=:white, yaxis=false,
                      linecolor=:transparent, c=:steelblue; kw...)
    end
  end

  if trace
    accRate = length(unique(x)) / length(x)
		traceTitle = showAccRate ? "acc: $(Int(round(accRate*100)))%" : ""
    #                                x    y   W   H
    plot!(x, inset_subplots=(parent, bbox(0, .1, .3, .3, :top, :right)), subplot=parent+offset,
          axis=false, legend=false, title=traceTitle,
          c=:grey70, grid=false, titlefont=traceFont, linewidths=.5)

  end

  xmean = mean(x)
  xsd = std(x)
  q = quantile(x, [a/2, 1-a/2])

  xticks!(round.([q; xmean], digits=q_digits), axiscolor=:transparent; kw...)
  #vline!(q, c=:red, linewidths=1, line=:dot, label="95% CI", legend=:best; kw...)
  #vline!([xmean], c=:red, linewidths=2, legend=:best,
  #       label="Mean: $(round(xmean, digits=q_digits))"; kw...)
  #hline!([0], label="SD: $(round(xsd, digits=sd_digits))",
  #       linewidth=0, bgcolor_legend=:transparent; kw...)
  vline!(q, c=:red, label="", linewidths=1, line=:dot; kw...)
  vline!([xmean], c=:red, label="", linewidths=1; kw...)
  hline!([0], linewidth=0, label="", bgcolor_legend=:transparent; kw...)

  # TODO:
  # Option for traceplot

  return img
end

function plotPosts(X::Matrix{T}; a::Float64=0.05, q_digits::Int=3,
                   sd_digits::Int=3, cor_digits::Int=3, titles=:none,
                   detail=false, hist2d::Bool=false, trace::Bool=true,
                   traceFont=font(16), kw...) where {T <: Number}
  N = size(X, 2)

  if titles == :none
    titles = fill("", N)
  end

  img = plot(;layout=(N, N))
  counter = 0
  offset = 0
  for r in 1:N
    for c in 1:N
      counter += 1

      if r == c
        offset += 1
        plotPost(X[:, r], a=a, q_digits=q_digits, sd_digits=sd_digits,
                 subplot=counter, add=true, title=titles[c], trace=trace,
                 parent=counter, offset=N*N-counter+offset,
                 traceFont=traceFont; kw...)
      elseif r < c
        if hist2d
          if detail
            histogram2d!(X[:, c], X[:,r], subplot=counter, colorbar=:none, grid=false)
          else
            histogram2d!(X[:, c], X[:,r], subplot=counter, colorbar=:none, grid=false,
                         axis=false)
          end
        else
          plot!(X[:, c], X[:,r], subplot=counter, c=:grey70, grid=false,
                axis=false, legend=false, linewidths=.5)
        end
      else # r > c
        corXrc = cor(X[:, r], X[:, c])
        annotate!([(0.5, 0.5, text("r=$(round(corXrc, digits=cor_digits))",
                                   Int(ceil(abs(corXrc) * 5))+5, :center))],
                  subplot=counter, axis=false, grid=false)
      end
    end
  end

  return img
end

function postProbMiss(b0, b1, i::Int;
                      y::Vector{Float64}=collect(-10:.1:10),
                      credibility::Float64=.95)

  N, I = size(b0)
  @assert size(b0) == size(b1)

  len_y = length(y)
  alpha = 1 - credibility
  p_lower = alpha / 2
  p_upper = 1 - alpha / 2

  pmiss = hcat([Cytof5.Model.prob_miss.(yi, b0[:,i], b1[:,i]) for yi in y]...)
  pmiss_mean = vec(mean(pmiss, dims=1))
  pmiss_lower = [ quantile(pmiss[:, col], p_lower) for col in 1:len_y ]
  pmiss_upper = [ quantile(pmiss[:, col], p_upper) for col in 1:len_y ]

  return (pmiss_mean, pmiss_lower, pmiss_upper, y)
end

function pairwiseAlloc(Z, W, i)
  J, K = size(Z)
  A = zeros(J, J)
  for r in 1:J
    for c in 1:J
      for k in 1:K
        if Z[r, k] == 1 && Z[c, k] == 1
          A[r, c] += W[i, k]
          A[c, r] += W[i, k]
        end
      end
    end
  end

  return A
end # pairwiseAlloc

function estimate_ZWi_index(monitor, i)
  As = [pairwiseAlloc(m[:Z], m[:W], i) for m in monitor]

  Amean = mean(As)
  mse = [ mean((A - Amean) .^ 2) for A in As]

  return argmin(mse)
end

function countmap(x::Vector{T}) where {T <: Number}
  out = Dict{Int, T}()
  for xi in x
    if xi in keys(out)
      out[xi] += 1
    else
      out[xi] = 1
    end
  end
  return out
end

function reorder_lam_i(lam_i::Vector{T}, W_i::Vector{V}) where {T <: Number, V <: Number}
  N = length(lam_i)
  new_lam_i = Vector{Int}(undef, N)
  K = length(W_i)
  ord = sortperm(W_i)

  for k in 1:K
    new_lam_i[lam_i .== ord[k]] .= k
  end

  return new_lam_i
end

# TODO: Implement the threshold
function yZ_inspect(monitor, y::Vector{Matrix{Float64}},
                    i::Int; thresh=0.7, ycolor=:balance,
                    addLines::Bool=true, marker_names = missing,
                    clim=(-5,5), W_digits=1, kw...)
  idx_best = estimate_ZWi_index(monitor, i)
  Zi = monitor[idx_best][:Z]
  Wi = monitor[idx_best][:W][i, :]
  lami = monitor[idx_best][:lam][i]
  new_lami = reorder_lam_i(lami, Wi)
  ord = sortperm(new_lami)
  K = length(Wi)
  J = size(y[i], 2)

  W_order = sortperm(Wi)

  cumProb = 0.0; k = K
  for w in reverse(Wi[W_order])
    cumProb += w
    if cumProb > thresh
      break
    end
    k -= 1
  end
  W_order = W_order[k:end]

  plotYZ(y[i][ord, :], Zi[:, W_order], clim=clim, ycolor=ycolor)

  # TODO: cell types on left, percentage on right
  W_val = map(x->"$(x)%", round.(Wi[W_order] * 100, digits=W_digits))
  #plot!(yticks=(1:K, W_order), subplot=2)
  plot!(yticks=(1:length(W_order), W_val), subplot=2, ymirror=true)

  # add marker names if provided
  if !ismissing(marker_names)
    plot!(xticks=(1:J, marker_names), subplot=1, rotation=90)
  end

  if addLines
    groups = cumsum(sort(collect(values(countmap(new_lami)))))
    hline!(groups, c=:yellow, subplot=1, linewidths=1, legend=false)
  end
end

function postSummary(x::Matrix{T}; alpha::Float64=.05) where {T <: Number}
  K = size(x, 2)
  xmean = mean(x, dims=1)
  xsd = std(x, dims=1)
  xqLower = [ quantile(x[:, k], alpha / 2) for k in 1:K ]
  xqLower = reshape(xqLower, 1, K)
  xqUpper = [ quantile(x[:, k], 1- alpha / 2) for k in 1:K ]
  xqUpper = reshape(xqUpper, 1, K)
  xmedian = median(x, dims=1)

  summaryMatrix = Matrix([xmean; xsd; xqLower; xmedian; xqUpper]')
  println("mean   sd  lower  upper  median")
  return summaryMatrix
end

# TODO: GET RID OF THIS!
R"""
ari <- function(x, y) {
  x = as.vector(x)
  y = as.vector(y)
  n = length(x)
  stopifnot(n == length(y))
  tab = table(x, y)
  index = sum(choose(tab, 2))
  rowSumsChoose2.sum = sum(choose(rowSums(tab), 2))
  colSumsChoose2.sum = sum(choose(colSums(tab), 2))
  expected.index = rowSumsChoose2.sum * colSumsChoose2.sum/choose(n,
      2)
  max.index = mean(c(rowSumsChoose2.sum, colSumsChoose2.sum))
  (index - expected.index)/(max.index - expected.index)
}
"""
ari = R"ari"

function confusion(x::Vector{T}, y::Vector{T}) where {T <: Number}
  @assert length(y) == length(x)
  n = length(x)
  confusionMat = zeros(Int, n, n)
  for i in 1:n
    for j in 1:n
      if x[i] == y[j]
        confusionMat[i, j] += 1
      end
    end
  end
  return confusionMat
end

#=
function ari(x::Vector{T}, y::Vector{T}) where {T <: Number}
  @assert length(y) == length(x)
  n = length(x)
  tab = confusion(x, y)
  # TODO: https://en.wikipedia.org/wiki/Rand_index
end
=#

end


text:
type Event
    ptr::Ptr{Void}

    function Event(ptr::Ptr{Void})
        event = new(ptr)
        finalizer(event, destroy)
        event
    end
end

immutable KeyEvent
    etype::Cint
    key_code::Cint
    alt::Cint
    control::Cint
    shift::Cint
    system::Cint
end

immutable TextEvent
    etype::Cint
    unicode::UInt32
end

immutable MouseMoveEvent
    etype::Cint
    x::Cint
    y::Cint
end

immutable MouseButtonEvent
    etype::Cint
    button::Cint
    x::Cint
    y::Cint
end

immutable MouseWheelEvent
    etype::Cint
    delta::Cint
    x::Cint
    y::Cint
end

immutable JoystickMoveEvent
    etype::Cint
    joystick_id::UInt32
    axis::JoystickAxis
    position::Cfloat
end

immutable JoystickButtonEvent
    etype::Cint
    joystick_id::UInt32
    button::UInt32
end

immutable JoystickConnectEvent
    etype::Cint
    joystick_id::UInt32;
end

immutable SizeEvent
    etype::Cint
    width::UInt32
    height::UInt32
end

baremodule EventType
    const CLOSED = 0
    const RESIZED = 1
    const LOST_FOCUS = 2
    const GAINED_FOCUS = 3
    const TEXT_ENTERED = 4
    const KEY_PRESSED = 5
    const KEY_RELEASED = 6
    const MOUSE_WHEEL_MOVED = 7
    const MOUSE_BUTTON_PRESSED = 8
    const MOUSE_BUTTON_RELEASED = 9
    const MOUSE_MOVED = 10
    const MOUSE_ENTERED = 11
    const MOUSE_LEFT = 12
    const JOYSTICK_BUTTON_PRESSED = 13
    const JOYSTICK_BUTTON_RELEASED = 14
    const JOYSTICK_MOVED = 15
    const JOYSTICK_CONNECTED = 16
    const JOYSTICK_DISCONNECTED = 17
    const TOUCH_BEGAN = 18
    const TOUCH_MOVED = 19
    const TOUCH_ENDED = 20
    const SENSOR_CHANGED = 21
end

function Event()
    Event(ccall((:new_sjEvent, "libjuliasfml"), Ptr{Void}, ()))
end

function destroy(event::Event)
    ccall((:sjEvent_destroy, "libjuliasfml"), Void, (Ptr{Void},), event.ptr)
end

function get_type(event::Event)
    return ccall((:sjEvent_eventType, "libjuliasfml"), Int32, (Ptr{Void},), event.ptr)
end

function get_size(event::Event)
    return ccall((:sjEvent_eventSize, "libjuliasfml"), SizeEvent, (Ptr{Void},), event.ptr)
end

function get_key(event::Event)
    return ccall((:sjEvent_eventKey, "libjuliasfml"), KeyEvent, (Ptr{Void},), event.ptr)
end

function get_text(event::Event)
    return ccall((:sjEvent_eventText, "libjuliasfml"), TextEvent, (Ptr{Void},), event.ptr)
end

function get_mousebutton(event::Event)
    return ccall((:sjEvent_eventMouseButton, "libjuliasfml"), MouseButtonEvent, (Ptr{Void},), event.ptr)
end

function get_mousemove(event::Event)
    return ccall((:sjEvent_eventMouseMove, "libjuliasfml"), MouseMoveEvent, (Ptr{Void},), event.ptr)
end

function get_mousewheel(event::Event)
    return ccall((:sjEvent_eventMouseWheel, "libjuliasfml"), MouseWheelEvent, (Ptr{Void},), event.ptr)
end


text:
module SpelledOut

using DecFP: Dec64
using Formatting: format

export spelled_out, Spelled_out, Spelled_Out, SPELLED_OUT

include("en.jl")
include("es.jl")

"""
```julia
function spelled_out(
    number::Number;
    lang::Symbol = <locale language>,
    dict::Symbol = :modern
)
```

Spells out a number in words, in lowercase, given a specified language.  The default language is the one used on your system.  The `dict` keyword argument only applies to some languages; see "Supported Languages" in the docs for more information.

---

### Examples

```julia
julia> spelled_out(12, lang = :en)
"twelve"

julia> spelled_out(112, lang = :en)
"one hundred twelve"

julia> spelled_out(112, lang = :en_GB)
"one hundred and twelve"

julia> spelled_out(112, lang = :en, dict = :european)
"one hundred twelve"

julia> spelled_out(112, lang = :es)
"ciento doce"
```
"""
function spelled_out(
    number::Number;
    lang::Symbol = Symbol(ENV["LANG"][1:(findfirst(==('.'), ENV["LANG"]) - 1)]), # Symbol(first(split(ENV["LANG"], '.')))
    dict::Symbol = :modern
)

    if lang ∈ (:en, :en_US)
        return spelled_out_en(number, british = false, dict = dict)
    elseif lang ∈ (:en_UK, :en_GB, :en_NZ, :en_AU)
        return spelled_out_en(number, british = true, dict = dict)
    elseif lang ∈ (:es,)
        return spelled_out_es(number,dict)
    end

    throw(error("We do not support $lang yet.  Please make an issue and someone might be able to help you, or feel free to contribute."))

    return nothing
end

end # end module


text:
# Load a common workspace and solve a lasso problem
using MAT
vars = matread("../compare/vars_lasso.mat")

A = vars["A"]
b = vars["b"]
tau = vars["tau"]


text:
function GeoMagneticModelLoader()
    return GeoMagneticModelLoader(())
end

function get_models(obj::GeoMagneticModelLoader)
    return jcall(obj, "getModels", Collection, ())
end

function load_data(obj::GeoMagneticModelLoader, arg0::InputStream, arg1::JString)
    return jcall(obj, "loadData", void, (InputStream, JString), arg0, arg1)
end

function still_accepts_data(obj::GeoMagneticModelLoader)
    return jcall(obj, "stillAcceptsData", jboolean, ())
end


text:
"""
`SavingsAccount` is a kind of `Account` that generates interest daily.
We will use Delegation pattern to gain existing `Account` functionality.
"""
struct SavingsAccount
    acct::Account
    interest_rate::Float64

    SavingsAccount(account_number, balance, date_opened, interest_rate) = new(
        Account(account_number, balance, date_opened),
        interest_rate
    )
end

# Forward assessors
account_number(sa::SavingsAccount) = account_number(sa.acct)
balance(sa::SavingsAccount) = balance(sa.acct)
date_opened(sa::SavingsAccount) = date_opened(sa.acct)

deposit!(sa::SavingsAccount, amount::Real) = deposit!(sa.acct, amount)
withdraw!(sa::SavingsAccount, amount::Real) = withdraw!(sa.acct, amount)
transfer!(sa1::SavingsAccount, sa2::SavingsAccount, amount::Real) = transfer!(
    sa1.acct, sa2.acct, amount)

# new accessor
interest_rate(sa::SavingsAccount) = sa.interest_rate

# new behavior
function accrue_daily_interest!(sa::SavingsAccount)
    interest = balance(sa.acct) * interest_rate(sa) / 365
    deposit!(sa.acct, interest)
end


text:
using ABBAj, Test, DelimitedFiles
time_series = load_sample();

# test 1
symbols, model = fit_transform(time_series,20, 0.1);
r_time_series = inverse_transform(symbols, model, time_series[1]);

write("symbols1", String(symbols));
symbols_backup = read("symbols1", String);
@test symbols_backup == String(symbols)

writedlm("r_time_series1", r_time_series)
r_time_series_backup = readdlm("r_time_series1", '\t', Float64, '\n');
@test r_time_series_backup == hcat(r_time_series)

# test 2
symbols, model = fit_transform(time_series,30, 0.1);
r_time_series = inverse_transform(symbols, model, time_series[1]);

write("symbols2", String(symbols));
symbols_backup = read("symbols2", String);
@test symbols_backup == String(symbols)

writedlm("r_time_series2", r_time_series)
r_time_series_backup = readdlm("r_time_series2", '\t', Float64, '\n');
@test r_time_series_backup == hcat(r_time_series)

# test 3
symbols, model = fit_transform(time_series,30, 0.5);
r_time_series = inverse_transform(symbols, model, time_series[1]);

write("symbols3", String(symbols));
symbols_backup = read("symbols3", String);
@test symbols_backup == String(symbols)

writedlm("r_time_series3", r_time_series)
r_time_series_backup = readdlm("r_time_series3", '\t', Float64, '\n');
@test r_time_series_backup == hcat(r_time_series)

# test 4
symbols, model = fit_transform(time_series,30, 0.1, 5);
r_time_series = inverse_transform(symbols, model, time_series[1]);

write("symbols4", String(symbols));
symbols_backup = read("symbols4", String);
@test symbols_backup == String(symbols)

writedlm("r_time_series4", r_time_series)
r_time_series_backup = readdlm("r_time_series4", '\t', Float64, '\n');
@test r_time_series_backup == hcat(r_time_series)

max_stars_repo_name: nla-group/Julia-ABBA

max_issues_repo_name: nla-group/Julia-ABBA

max_forks_repo_name: nla-group/Julia-ABBA


text:
using CauViz
using Test