text:
 using Pkg
 Pkg.instantiate()
 Pkg.activate(".")
 using Psychometrics
 using Distributions
 using LinearAlgebra
 using Dates
 using Random
 using SharedArrays
 using JLD2
 using StatsBase
 using DelimitedFiles
 using FileIO
reps=100
folder="test/simulation/batch50 NT400 Dgain"
items = load(string(folder,"/true values/true_items.jld2"), "items")

a_s = zeros(I_field,reps);
b_s = zeros(I_field,reps);
for rep = 1:reps
    # if isfile(string("test/online-calibration/rep_", rep, "_items_est.jld2"))
        # @load string("test/online-calibration/rep_", rep, "_items_est.jld2") items_est
        # a = map(i -> i.parameters.a, items_est)
        # a_s[:, rep] = copy(a)
        # b = map(i -> i.parameters.b, items_est)
        # b_s[:, rep] = copy(b)
    # else
        @load string(folder,"/rep_", rep, "_items_est_field.jld2") items_est_field
        a = map(i -> i.parameters.a, items_est_field)
        a_s[:, rep] = copy(a)
        b = map(i -> i.parameters.b, items_est_field)
        b_s[:, rep] = copy(b)
    #end

end

a_true = map(i -> i.parameters.a, items[(I_operational + 1) : I_total]);
a_true = hcat([a_true for rep = 1:reps]...);
b_true = map(i -> i.parameters.b, items[(I_operational + 1) : I_total]);
b_true = hcat([b_true for rep=1:reps]...);
diff_a = (a_s - a_true);
diff_b = (b_s - b_true);

using CSV
using DataFrames
a_true = map(i -> i.parameters.a, items[(I_operational + 1) : I_total]);
b_true = map(i -> i.parameters.b, items[(I_operational + 1) : I_total]);

RMSE_df = DataFrame(a_true = Vector{Float64}(a_true), b_true = Vector{Float64}(b_true), a_RMSE=[sqrt(sum(r.^2)/reps) for r in eachrow(diff_a)],b_RMSE=[sqrt(sum(r.^2)/reps) for r in eachrow(diff_b)]);
CSV.write(
    string(folder,"/RMSE.csv"),
    RMSE_df
)

CSV.write(
    string(folder,"/BIAS.csv"),
    DataFrame(a_BIAS=[sum(r.^2)/reps for r in eachrow(diff_a)],
    b_BIAS=[sum(r)/reps for r in eachrow(diff_b)])
)
using Plots
p = plot(layout=(5,5));
p_a = plot(ylim=(0,maximum(RMSE_df.a_RMSE)+0.15), xticks = (round.(unique(b_true), digits=2), string.(round.(unique(b_true), digits=2))), xlabel= "b", legend= :best);
global i=1
for a in unique(a_true)
    p_a= plot!(p_a, RMSE_df[RMSE_df.a_true .== a,:].b_true, RMSE_df[RMSE_df.a_true .== a,:].a_RMSE, label="a = " .*string.(round(a, digits=2)));
    i+=1
end
plot(p_a)
savefig(p_a,string(folder,"/a_RMSE.pdf"))
p_b = plot(ylim=(0,0.25), xticks = (round.(unique(a_true), digits=2), string.(round.(unique(a_true), digits=2))), xlabel= "a", legend= :best);
global i=1
for b in unique(b_true)
    p_b= plot!(p_b, RMSE_df[RMSE_df.b_true .== b,:].a_true, RMSE_df[RMSE_df.b_true .== b,:].b_RMSE, label="b = " .*string.(round(b, digits=2)));
    i+=1
end
plot(p_b)
savefig(p_b,string(folder,"/b_RMSE.pdf"))

println(sqrt(sum(diff_a.^2)/I_field/reps))
println(sqrt(sum(diff_b.^2)/I_field/reps))


text:
function setTtlMode(rhd::FPGA,mode)

    value = 0;
    value += mode[1] ? 1 : 0;
    value += mode[2] ? 2 : 0;
    value += mode[3] ? 4 : 0;
    value += mode[4] ? 8 : 0;
    value += mode[5] ? 16 : 0;
    value += mode[6] ? 32 : 0;
    value += mode[7] ? 64 : 0;
    value += mode[8] ? 128 : 0;

    SetWireInValue(rhd,WireInTtlOut,value,0x000000ff)
    UpdateWireIns(rhd)
end

function manual_trigger(rhd,trigger,triggerOn)
    #Trigger 0-7
    SetWireInValue(rhd,WireInManualTriggers, (triggerOn ? 1 : 0) << trigger,1<<trigger);
    UpdateWireIns(rhd)
end

clearTtlOut(rhd::FPGA)=(SetWireInValue(rhd,WireInTtlOut, 0x0000);UpdateWireIns(rhd))

function getTtlIn(rhd::FPGA,ttlInArray)

    UpdateWireOuts(rhd)
    ttlIn=GetWireOutValue(rhd,WireOutTtlIn)
    for i=1:16
        ttlInArray[i] = 0
        if (ttlIn & (1 << (i-1))) > 0
            ttlInArray[i] = 1
        end
    end

    ttlInArray
end

function setTtlOut(rhd::FPGA,ttlOutArray)

    ttlOut=Int32(0)
    for i=1:16
        if ttlOutArray[i]>0
            ttlOut += (1<< (i-1))
        end
    end

    SetWireInValue(rhd,WireInTtlOut,ttlOut)
    UpdateWireIns(rhd)
    nothing
end

function sendTimePulse(fpga::FPGA,val::Bool)

    if val==true
        fpga.ttloutput += (1 << (16-1))
    else
        fpga.ttloutput -= (1 << (16-1))
    end

    SetWireInValue(fpga,WireInTtlOut,fpga.ttloutput)
    UpdateWireIns(fpga)
    nothing
end

#=
Fast settle function allows TTL inputs to black amplifier channels
Can be used to prevent amplifiers from saturation (such as during stimulation)
=#
function enableExternalFastSettle(rhd::FPGA,enable)

    SetWireInValue(rhd,WireInMultiUse, (enable ? 1 : 0))
    UpdateWireIns(rhd)
    if (OPEN_EPHYS)
        ActivateTriggerIn(rhd,TrigInExtFastSettle,0)
    else
        ActivateTriggerIn(rhd,TrigInConfig,6)
    end

    nothing
end

function setExternalFastSettleChannel(rhd::FPGA,channel)

    #error checking goes here

    SetWireInValue(rhd,WireInMultiUse,channel)
    UpdateWireIns(rhd)
    if (OPEN_EPHYS)
        ActivateTriggerIn(rhd,TrigInExtFastSettle,1)
    else
        ActivateTriggerIn(rhd,TrigInConfig,7)
    end
end

ext: jl


text:
using Base.Test
using MathProgComplex

@testset "WB2 real formulation, fixed phase, order 2" begin
    sols = SortedDict(0.976 => (2, 905.76, 905.76),
                    0.983 => (2, 905.73, 903.12),
                    0.989 => (2, 905.73, 900.84),
                    0.996 => (2, 905.73, 898.17),
                    1.002 => (2, 905.73, 895.86),
                    1.009 => (2, 905.73, 893.16),
                    1.015 => (2, 905.73, 890.82),
                    1.022 => (3, 905.73, 888.08),
                    1.028 => (3, 905.73, 885.71),
                    1.035 => (2, 882.97, 882.97))

    testfolder = joinpath(Pkg.dir("MathProgComplex"), "Mosek_runs", "tests")
    @testset "v2max=$v2max, rmeqs=$rmeqs, addball=$addball" for (v2max, (dmax, obj_fixedphase, obj_rankrel)) in sols, rmeqs in Set([true, false]), addball in Set([true, false])

        # info("Working on WB2, no eqs, v2max=$v2max, d=$d")
        problem = buildPOP_WB2(v2max=v2max, rmeqs=rmeqs, setnetworkphase=true, addball=addball)

        logpath = joinpath(testfolder, "WB2_v2max_$(v2max)_rankrel")
        ispath(logpath) && rm(logpath, recursive=true); mkpath(logpath)
        println("Saving file at $logpath")

        relax_ctx = set_relaxation(problem; hierarchykind=:Real,
                                            d = dmax)

        primobj, dualobj = run_hierarchy(problem, relax_ctx, logpath, save_pbs=true);
        @show primobj, dualobj, obj_fixedphase
        if (0.983 ≤ v2max ≤ 1.028)
            @test_broken primobj ≈ obj_fixedphase atol=1e-2
            @test_broken dualobj ≈ primobj atol=mosek_optgap*min(abs(primobj), abs(dualobj))
        else
            @test primobj ≈ obj_fixedphase atol=1e-2
            @test dualobj ≈ primobj atol=mosek_optgap*min(abs(primobj), abs(dualobj))
        end
    end
end


text:
#=---------
REGBANK A
----------=#
module regbankA
#=----------------------------
Initialize registers to zero
-----------------------------=#
 s0 = 0
 s1 = 0
 s2 = 0
 s3 = 0
 s4 = 0
 s5 = 0
 s6 = 0
 s7 = 0
 s8 = 0
 s9 = 0
 sA = 0
 sB = 0
 sC = 0
 sD = 0
 sE = 0
 sF = 0

#=-------------------------------------
Text parsing stuff (you can ignore :])
--------------------------------------=#
registers = Dict("S0" => s0,"S1" => s1,"S2" => s2,"S3" => s3,"S4" => s4,"S5" => s5,"S6" => s6,
"S7" => s7,"S8" => s8,"S9" => s9,"SA" => sA,"SB" => sB,"SC" => sC,"SD" => sD,"SE" => sE,"SF" => sF)
#=----------------------------------=#
#=----Getters and setters allow access to registers from other modules----=#
function get(sX)
return registers[sX]
end

function set(sX,value)
  registers[sX] = value
end

end

#=--------------
REGBANK B
---------------=#
module regbankB #Registers in bank A
#=----------------------------
Initialize registers to zero
-----------------------------=#
s0 = 0
s1 = 0
s2 = 0
s3 = 0
s4 = 0
s5 = 0
s6 = 0
s7 = 0
s8 = 0
s9 = 0
sA = 0
sB = 0
sC = 0
sD = 0
sE = 0
sF = 0

#=-------------------------------------
Text parsing stuff (you can ignore :])
--------------------------------------=#
registers = Dict("S0" => s0,"S1" => s1,"S2" => s2,"S3" => s3,"S4" => s4,"S5" => s5,"S6" => s6,
"S7" => s7,"S8" => s8,"S9" => s9,"SA" => sA,"SB" => sB,"SC" => sC,"SD" => sD,"SE" => sE,"SF" => sF)
#=----------------------------------=#

function get(sX)
  return registers[sX]
end

function set(sX,value)
  registers[sX] = value
end

end

max_stars_repo_name: kinoroy/KCPSM-Simulator

max_issues_repo_name: kinoroy/KCPSM-Simulator

max_forks_repo_name: kinoroy/KCPSM-Simulator


text:
using BinDeps
using Compat

@BinDeps.setup

libopenni2 = library_dependency("libOpenNI2", aliases=["libOpenNI2"])

provides(AptGet, "libopenni2-dev", libopenni2)

#=
@osx_only begin
    using Homebrew
    Homebrew.add("homebrew/science/openni2")
    provides(Homebrew.HB, "openni2", openni2)
end
=#

@BinDeps.install Dict(:libOpenNI2 => :libOpenNI2)


text:
# From https://gist.github.com/mcabbott/80ac43cca3bee8f57809155a5240519f
function _repeat(x::AbstractArray, counts::Integer...)
    N = max(ndims(x), length(counts))
    size_y = ntuple(d -> size(x,d) * get(counts, d, 1), N)
    size_x2 = ntuple(d -> isodd(d) ? size(x, 1+d÷2) : 1, 2*N)

    ## version without mutation
    # ignores = ntuple(d -> reshape(Base.OneTo(counts[d]), ntuple(_->1, 2d-1)..., :), length(counts))
    # y = reshape(broadcast(first∘tuple, reshape(x, size_x2), ignores...), size_y)

    # ## version with mutation
    size_y2 = ntuple(d -> isodd(d) ? size(x, 1+d÷2) : get(counts, d÷2, 1), 2*N)
    y = similar(x, size_y)
    reshape(y, size_y2) .= reshape(x, size_x2)
    y
end

function ChainRulesCore.rrule(::typeof(_repeat), x::AbstractArray, counts::Integer...)
    size_x = size(x)
    function repeat_pullback_1(dy_raw)
        dy = unthunk(dy_raw)
        size2ndims = ntuple(d -> isodd(d) ? get(size_x, 1+d÷2, 1) : get(counts, d÷2, 1), 2*ndims(dy))
        reduced = sum(reshape(dy, size2ndims); dims = ntuple(d -> 2d, ndims(dy)))
        return (NoTangent(), reshape(reduced, size_x), map(_->NoTangent(), counts)...)
    end
    return _repeat(x, counts...), repeat_pullback_1
end


text:
#
# Advent of Code
# Day 13: A Maze of Twisty Little Cubicles
#

println("Advent of Code")
println("Day 13: A Maze of Twisty Little Cubicles")
println("")

import Base.Collections: PriorityQueue, enqueue!, dequeue!, peek
import Base.Order.Reverse

DEBUG = false

# test case
# SALT = 10
# (STARTX, STARTY) = (2,2)
# (GOALX, GOALY) = (8,5)
# MAZE_SIZE = 10

# part 1
SALT = 1362
(STARTX, STARTY) = (2,2)
(GOALX, GOALY) = (32,40)
MAZE_SIZE = 50

DIRS_X = [1, 0, -1, 0]
DIRS_Y = [0, 1, 0, -1]

MAZE = Dict([0=>".", 1=>"#", 2=>"O", 3=>"S", 4=>"F"])

function generateMaze(width::Int, height::Int)
  newMaze = Array{Int}(width, height)

  for y in 0:(height - 1)
    for x in 0:(width - 1)
      seed = x*x + 3*x + 2*x*y + y + y*y + SALT
      if isodd(sum([parse(string(item)) for item in bin(seed)]))
        # odd? wall.
        newMaze[x+1,y+1] = 1
      else
        # even? open.
        newMaze[x+1,y+1] = 0
      end
    end
  end
  return newMaze
end

function printMaze(maze::Array{Int})
  (rows, columns) = size(maze)
  for y in 1:rows
    row = ""
    for x in 1:columns
      row *= MAZE[maze[x,y]]
    end
    println(row)
  end
end

type Node
  x::Int
  y::Int
  distance::Float64
  priority::Float64
end

function updatePriority(node, goalx, goaly)
  node.priority = node.distance + sqrt((goalx - node.x)^2 + (goaly - node.y)^2) * 10.
end

function findShortestPath(startx::Int, starty::Int, endx::Int, endy::Int, maze::Array{Int})
  #
  # Implements the A* path-finding algorithm, based on example found at
  # http://code.activestate.com/recipes/577519-a-star-shortest-path-algorithm/
  #
  mazeWidth, mazeHeight = size(maze)

  evaluatedNodes = zeros(Float64, (mazeWidth, mazeHeight)) # "closed nodes map"
  potentialNodes = zeros(Float64, (mazeWidth, mazeHeight)) # "open nodes map"
  directions = ones(Int64, (mazeWidth, mazeHeight)) # directions to end

  priorityQueues = Vector{PriorityQueue}()                     # reverse priority queues because
  push!(priorityQueues, PriorityQueue(Node, Float64, Reverse)) # we want to take the highest
  push!(priorityQueues, PriorityQueue(Node, Float64, Reverse)) # priority, not the lowest

  # initialization
  priorityQueueIndex = 1

  startNode = Node(startx, starty, 0., 0.)
  updatePriority(startNode, endx, endy)
  enqueue!(priorityQueues[priorityQueueIndex], startNode, startNode.priority)
  potentialNodes[startx, starty] = startNode.priority

  # A* search
  while length(priorityQueues[priorityQueueIndex]) > 0
    # get next node with highest priority
    nextNode = deepcopy(dequeue!(priorityQueues[priorityQueueIndex]))
    (x, y) = (nextNode.x, nextNode.y)
    potentialNodes[x, y] = 0.
    evaluatedNodes[x, y] = 1.

    # are we there yet?
    if x == endx && y == endy
      # done, work backwards from finish to start
      shortestPath = 0
      maze[x, y] = 4 # finish symbol
      while !(x == startx && y == starty)
        direction = directions[x,y]
        shortestPath += 1
        x -= DIRS_X[direction]
        y -= DIRS_Y[direction]
        maze[x,y] = 2 # path symbol
      end
      maze[x,y] = 3 # start symbol
      return (maze, shortestPath)
    end

    # where to move now?
    for direction in 1:4
      # try each direction
      newx = x + DIRS_X[direction]
      newy = y + DIRS_Y[direction]

      # eliminate disallowed places:
      #    outside the maze
      #    wall
      #    already evaluated

      if !(newx < 1 || newx > mazeWidth || newy < 1 || newy > mazeHeight || maze[newx, newy] == 1 || evaluatedNodes[newx, newy] == 1)
        # make a new node
        newNode = Node(newx, newy, nextNode.distance + 10., nextNode.priority)
        updatePriority(newNode, endx, endy)

        if DEBUG println("New node at ", newx, ", ", newy, " with priority ", newNode.priority) end

        if potentialNodes[newx, newy] == 0.
          # if this wasn't a potential node before, assign the potential w/ new priority
          potentialNodes[newx, newy] = newNode.priority

          # put the node in the priority queue
          enqueue!(priorityQueues[priorityQueueIndex], newNode, newNode.priority)

          # indicate how we got here
          directions[newx, newy] = direction

        elseif potentialNodes[newx, newy] > newNode.priority
          # if we're putting a higher priority node here, swap for it
          potentialNodes[newx, newy] = newNode.priority

          # how did we get here
          directions[newx, newy] = direction

          # swap the old node for this higher-priority node
          while !(peek(priorityQueues[priorityQueueIndex]).x == newx && peek(priorityQueues[priorityQueueIndex]).y == newy)
            tempNode = dequeue!(priorityQueues[priorityQueueIndex])
            enqueue!(priorityQueues[3-priorityQueueIndex], tempNode)
          end

          tempNode = dequeue!(priorityQueues[priorityQueueIndex]) # remove the node

          if length(priorityQueues[priorityQueueIndex]) > length(priorityQueues[1-priorityQueueIndex])
            priorityQueueIndex = 3 -  priorityQueueIndex
          end

          while length(priorityQueues[priorityQueueIndex]) > 0
            tempNode = dequeue!(priorityQueues[priorityQueueIndex])
            enqueue!(priorityQueues[3-priorityQueueIndex], tempNode)
          end

          priorityQueueIndex = 3 -  priorityQueueIndex

          # put the better node on
          enqueue!(priorityQueues[priorityQueueIndex], newNode)
        end
      end
    end
  end

  # didn't find solution

  return (maze, 99999)

end

maze = generateMaze(MAZE_SIZE, MAZE_SIZE)

(maze, shortestPath) = findShortestPath(STARTX, STARTY, GOALX, GOALY, maze)
println("The shortest path is $shortestPath.")
printMaze(maze)

# part 2 is trying to find all the places we can reach within a distance
# of 50 steps

# Just brute force it.

maze = generateMaze(MAZE_SIZE, MAZE_SIZE)

nPaths = 0

for currentx in 1:50
  for currenty in 1:50
    if maze[currentx, currenty] == 0
      validMaze = deepcopy(maze)
      (maze, shortestPath) = findShortestPath(STARTX, STARTY, currentx, currenty, maze)
      if !(shortestPath <= 50)
        maze = validMaze
      else
        nPaths += 1
      end
    end
  end
end

# Now count non-zero and non-one entries

nReached = 0
for currentx in 1:50
  for currenty in 1:50
    if maze[currentx, currenty] > 1
        nReached += 1
    end
  end
end

println("")
println("$nReached distinct locations can be reached within 50 steps, by $nPaths paths.")

printMaze(maze)


text:
using Documenter, Genie
using App, AppServer, Authentication, Authorization
using Cache, Commands, Configuration, Cookies, Channels
using DatabaseSeeding, Error, FileTemplates, Generator
using Helpers, Inflector, Input, Logger, Macros, Migration
using Renderer, REPL, Router, Sessions, Tester, Toolbox, Util
using FileCacheAdapter, FileSessionAdapter, Encryption

push!(LOAD_PATH,  "../../src",
                  "../../src/cache_adapters",
                  "../../src/session_adapters")

makedocs()


text:
module NaiveGAflux

using Base: release
using Reexport
@reexport using NaiveNASflux
using NaiveNASflux: FluxDense, FluxConv, FluxConvolutional, FluxNoParLayer, FluxParNorm, FluxRnn, FluxBatchNorm
using NaiveNASflux: nograd, layertype
using NaiveNASlib.Advanced, NaiveNASlib.Extend
import Flux
using Flux: Dense, Conv, ConvTranspose, DepthwiseConv, CrossCor, LayerNorm, BatchNorm, InstanceNorm, GroupNorm,
            MaxPool, MeanPool, Dropout, AlphaDropout, GlobalMaxPool, GlobalMeanPool, cpu, gpu, WeightDecay,
            SamePad, params
using Random
using Logging
using Statistics
import CUDA
import MemPool
import IterTools
import Functors
using Functors: @functor, functor, fmap
using Printf

using Setfield

# For temporary storage of program state for pause/resume type of operations
using Serialization

const rng_default = MersenneTwister(abs(rand(Int)))
const modeldir = "models"

# Fitness
export fitness, AbstractFitness, LogFitness, GpuFitness, AccuracyFitness, TrainThenFitness, TrainAccuracyFitness, MapFitness, EwmaFitness, TimeFitness, SizeFitness, AggFitness

# Candidate
export evolvemodel, AbstractCandidate, CandidateModel, CandidateOptModel, FittedCandidate, model, opt, lossfun

# Evolution
export evolve, AbstractEvolution, NoOpEvolution, AfterEvolution, EliteSelection, SusSelection, TournamentSelection, CombinedEvolution, EvolutionChain, PairCandidates, ShuffleCandidates, EvolveCandidates

# Population
export Population, generation

# misc types
export Probability, MutationShield, ApplyIf, RemoveIfSingleInput, PersistentArray, ShieldedOpt

# Iterators. These should preferably come from somewhere else, but I haven't found anything which fits the bill w.r.t repeatability over subsets
export RepeatPartitionIterator, SeedIterator, MapIterator, GpuIterator, BatchIterator, ShuffleIterator, TimedIterator, TimedIteratorStop, StatefulGenerationIter

# Persistence
export persist

# Vertex selection types
export AbstractVertexSelection, AllVertices, FilterMutationAllowed

# mutation types
export AbstractMutation, MutationProbability, WeightedMutationProbability, HighUtilityMutationProbability, LowUtilityMutationProbability, MutationChain, RecordMutation, LogMutation, MutationFilter, PostMutation, VertexMutation, NoutMutation, AddVertexMutation, RemoveVertexMutation, AddEdgeMutation, RemoveEdgeMutation, KernelSizeMutation, KernelSizeMutation2D, ActivationFunctionMutation, PostMutation, OptimizerMutation, LearningRateMutation, AddOptimizerMutation

# Crossover types
export AbstractCrossover, VertexCrossover, CrossoverSwap, OptimizerCrossover, LearningRateCrossover

# architecture spaces
export AbstractArchSpace, LoggingArchSpace, VertexSpace, NoOpArchSpace, ArchSpace, ConditionalArchSpace, RepeatArchSpace, ArchSpaceChain, ForkArchSpace, ResidualArchSpace, FunctionSpace, GlobalPoolSpace

#  Other search space types
export BaseLayerSpace, AbstractParSpace, SingletonParSpace, Singleton2DParSpace, ParSpace, ParSpace1D, ParSpace2D, CoupledParSpace, NamedLayerSpace, LoggingLayerSpace, DenseSpace, ConvSpace, BatchNormSpace, PoolSpace, LayerVertexConf, Shielded, ConcConf

#weight inits
export AbstractWeightInit, DefaultWeightInit, IdentityWeightInit, PartialIdentityWeightInit, ZeroWeightInit

# functions
export nparams

# Pre-built programs for fitting data
export AutoFlux

# Visulization
export PlotFitness, ScatterPop, ScatterOpt, MultiPlot, CbAll

include("util.jl")
include("shape.jl")
include("archspace.jl")
include("mutation.jl")
include("crossover.jl")
include("candidate.jl")
include("fitness.jl")
include("evolve.jl")
include("population.jl")
include("iterators.jl")
include("app/AutoFlux.jl")
include("visualize/callbacks.jl")

end # module


text:
# Copyright (c) 2021, Eric Sabo
# All rights reserved.
#
# This source code is licensed under the BSD-style license found in the
# LICENSE file in the root directory of this source tree.

using AbstractAlgebra
using Nemo

import LinearAlgebra: tr

function ⊕(A::T, B::T) where T <: Union{fq_nmod_mat, gfp_mat}
    base_ring(A) == base_ring(B) || error("Matrices must be over the same base ring in directsum.")

    return vcat(hcat(A, zero_matrix(base_ring(B), size(A, 1), size(B, 2))),
        hcat(zero_matrix(base_ring(A), size(B, 1), size(A, 2)), B))
end

directsum(A::T, B::T) where T <: Union{fq_nmod_mat, gfp_mat} = A ⊕ B
⊗(A::T, B::T) where T <: Union{fq_nmod_mat, gfp_mat} = kronecker_product(A, B)
kron(A::T, B::T) where T <: Union{fq_nmod_mat, gfp_mat} = kronecker_product(A, B)
tensorproduct(A::T, B::T) where T <: Union{fq_nmod_mat, gfp_mat} = kronecker_product(A, B)
kroneckerproduct(A::T, B::T) where T <: Union{fq_nmod_mat, gfp_mat} = kronecker_product(A, B)
# nrows(A::T) where T = size(A, 1)
# ncols(A::T) where T = size(A, 2)

# I think we should avoid length checking here and return it for entire matrix if given
# Hammingweight(v::T) where T <: Union{fq_nmod_mat, gfp_mat, Vector{S}} where S <: Integer = count(i->(i != 0), v)
function Hammingweight(v::T) where T <: Union{fq_nmod_mat, gfp_mat, Vector{S}} where S <: Integer
    count = 0
    for i in 1:length(v)
        if !iszero(v[i])
            count += 1
        end
    end
    return count
end
weight(v::T) where T <: Union{fq_nmod_mat, gfp_mat, Vector{S}} where S <: Integer = Hammingweight(v)
wt(v::T) where T <: Union{fq_nmod_mat, gfp_mat, Vector{S}} where S <: Integer = Hammingweight(v)
Hammingdistance(u::T, v::T) where T <: Union{fq_nmod_mat, gfp_mat, Vector{S}} where S <: Integer = Hammingweight(u .- v)
distance(u::T, v::T) where T <: Union{fq_nmod_mat, gfp_mat, Vector{S}} where S <: Integer = Hammingweight(u .- v)
dist(u::T, v::T) where T <: Union{fq_nmod_mat, gfp_mat, Vector{S}} where S <: Integer = Hammingweight(u .- v)

function tr(x::fq_nmod, K::FqNmodFiniteField, verify::Bool=false)
    L = parent(x)
    if verify
        # shouldn't need Int casting here but just in case...
        Int64(characteristic(L)) == Int64(characteristic(K)) || error("The given field is not a subfield of the base ring of the element.")
        degree(L) % degree(K) == 0 || error("The given field is not a subfield of the base ring of the element.")
    end
    n = div(degree(L), degree(K))
    return sum([x^(q^i) for i in 0:(n - 1)])
end

function _expandelement(x::fq_nmod, K::FqNmodFiniteField, basis::Vector{fq_nmod}, verify::Bool=false)
    return [tr(x * i, K, verify) for i in basis]
end

function _expandrow(row::Vector{fq_nmod}, K::FqNmodFiniteField, basis::Vector{fq_nmod}, verify::Bool=false)
    return hcat([_expandelement(row[i], K, basis, verify)]...)
end

function expandmatrix(M::fq_nmod_mat, K::FqNmodFiniteField, basis::Vector{fq_nmod})
    L = base_ring(M)
    L != K || return M
    # shouldn't need Int casting here but just in case...
    Int64(characteristic(L)) == Int64(characteristic(K)) || error("The given field is not a subfield of the base ring of the element.")
    degree(L) % degree(K) == 0 || error("The given field is not a subfield of the base ring of the element.")
    n = div(degree(L), degree(K))
    n == length(basis) || error("Provided basis is of incorrect size for the given field and subfield.")
    # should really check if it is a basis
    return vcat([_expandrow(M[r, :], K, basis) for r in 1:size(M, 1)])
end

function symplecticinnerproduct(u::fq_nmod_mat, v::fq_nmod_mat)
    (size(u, 1) == 1 || size(u, 2) == 1) || error("First argument of symplectic inner product is not a vector: dims = $(size(u, 1)).")
    (size(v, 1) == 1 || size(v, 2) == 1) || error("Second argument of symplectic inner product is not a vector: dims = $(size(v, 1)).")
    length(u) == length(v) || error("Vectors must be the same length in symplectic inner product.")
    iseven(length(u)) || error("Vectors must have even length in symplectic inner product.")
    base_ring(u) == base_ring(v) || error("Vectors must be over the same field in symplectic inner product.")
    ncols = div(length(u), 2)
    return sum([u[i + ncols] * v[i] - v[i + ncols] * u[i] for i in 1:ncols])
end
# SIP(u::fq_nmod_mat, v::fq_nmod_mat) = symplecticinnerproduct(u, v)

# since changed method here, need to add error checks on dimensions
function aresymplecticorthogonal(A::fq_nmod_mat, B::fq_nmod_mat)
    AEuc = hcat(A[:, div(size(A, 2), 2) + 1:end], -A[:, 1:div(size(A, 2), 2)])
    iszero(AEuc * transpose(B)) || return false
    return true
end

# function traceinnerproduct(u::fq_nmod_mat, v::fq_nmod_mat)
#
# end
# TIP(u::fq_nmod_mat, v::fq_nmod_mat) = traceinnerproduct(u, v)

function Hermitianinnerproduct(u::fq_nmod_mat, v::fq_nmod_mat)
    (size(u, 1) == 1 || size(u, 2) == 1) || error("First argument of Hermitian inner product is not a vector: dims = $(size(u, 1)).")
    (size(v, 1) == 1 || size(v, 2) == 1) || error("Second argument of Hermitian inner product is not a vector: dims = $(size(v, 1)).")
    length(u) == length(v) || error("Vectors must be the same length in Hermitian inner product.")
    base_ring(u) == base_ring(v) || error("Vectors must be over the same field in Hermitian inner product.")
    q2 = order(base_ring(u))
    issquare(q2) || error("The Hermitian inner product is only defined over quadratic field extensions.")
    q = Int64(sqrt(q2))
    return sum([u[i] * v[i]^q for i in 1:length(u)])
end
# HIP(u::fq_nmod_mat, v::fq_nmod_mat) = Hermitianinnerproduct(u, v)

function Hermitianconjugatematrix(A::fq_nmod_mat)
    B = copy(A)
    q2 = order(base_ring(A))
    issquare(q2) || error("The Hermitian conjugate is only defined over quadratic field extensions.")
    q = Int64(sqrt(q2))
    return B .^ q
end

function printstringarray(A::Vector{String}, withoutIs=false)
    for a in A
        if !withoutIs
            println(a)
        else
            for i in a
                if i == 'I'
                    print(' ')
                else
                    print(i)
                end
            end
            print('\n')
        end
    end
end
printchararray(A::Vector{Vector{Char}}, withoutIs=false) = printstringarray(setchartostringarray(A), withoutIs)
printsymplecticarray(A::Vector{Vector{T}}, withoutIs=false) where T <: Integer = printstringarray(setsymplectictostringarray(A), withoutIs)

function entropy(x::Real)
    x != 0 || return 0
    (0 < x <= 1 - 1 / q) || error("Number should be in the range [0, 1 - 1/order(field)].")
    F = parent(x)
    q = order(F)
    return x * (log(q, q - 1) - log(q, x)) - (1 - x) * log(q, 1 - x)
end

function FpmattoJulia(M::fq_nmod_mat)
    degree(base_ring(M)) == 1 || error("Cannot promote higher order elements to the integers.")
    # Fp = [i for i in 0:Int64(characteristic(base_ring(M)))]
    A = zeros(Int64, size(M))
    for r in 1:size(M, 1)
        for c in 1:size(M, 2)
            # A[r, c] = Fp[findfirst(x->x==M[r, c], Fp)]
            A[r, c] = coeff(M[r, c], 0)
        end
    end
    return A
end

function pseudoinverse(M::fq_nmod_mat, verify::Bool=true)
    # let this fail elsewhere if not actually over a quadratic extension
    if degree(base_ring(M)) != 1
        M = deepcopy(quadratictosymplectic(M)')
    else
        M = deepcopy(M')
    end

    nr, nc = size(M)
    MS = MatrixSpace(base_ring(M), nr, nr)
    _, E = rref(hcat(M, MS(1)))
    E = E[:, (nc + 1):end]
    pinv = E[1:nc, :]
    dual = E[nc + 1:nr, :]

    if verify
        _, Mrref = rref(M)
        MScols = MatrixSpace(base_ring(M), nc, nc)
        E * M == Mrref || error("Pseudoinverse calculation failed (transformation incorrect).")
        Mrref[1:nc, 1:nc] == MScols(1) || error("Pseudoinverse calculation failed (failed to get I).")
        iszero(Mrref[nc + 1:nr, :]) || error("Pseudoinverse calculation failed (failed to get zero).")
        pinv * M == MScols(1) || error("Pseudoinverse calculation failed (eq 1).")
        M' * pinv' == MScols(1) || error("Pseudoinverse calculation failed (eq 2).")
        iszero(M' * dual') || error("Failed to correctly compute dual (rhs).")
        iszero(dual * M) || error("Failed to correctly compute dual (lhs).")
    end
    return pinv
end

function istriorthogonal(G::fq_nmod_mat, verbose::Bool=false)
    Int(order(base_ring(G))) == 2 || error("Triothogonality is only defined over 𝔽_2.")
    nr, nc = size(G)
    for r1 in 1:nr
        for r2 in 1:nr
            @views g1 = G[r1, :]
            @views g2 = G[r2, :]
            @views if !iszero(sum([g1[1, i] * g2[1, i] for i in 1:nc]))
                verbose && println("Rows $r1 and $r2 are not orthogonal.")
                return false
            end
        end
    end

    for r1 in 1:nr
        for r2 in 1:nr
            for r3 in 1:nr
                @views g1 = G[r1, :]
                @views g2 = G[r2, :]
                @views g3 = G[r3, :]
                @views if !iszero(sum([g1[1, i] * g2[1, i] * g3[1, i] for i in 1:nc]))
                    verbose && println("Rows $r1, $r2, and $r3 are not orthogonal.")
                    return false
                end
            end
        end
    end
    return true
end

function istriorthogonal(G::Matrix{Int})
    nr, nc = size(G)
    for r1 in 1:nr
        for r2 in 1:nr
            @views g1 = G[r1, :]
            @views g2 = G[r2, :]
            @views if !iszero(sum([g1[1, i] * g2[1, i] for i in 1:nc]) % 2)
                return false
            end
        end
    end

    for r1 in 1:nr
        for r2 in 1:nr
            for r3 in 1:nr
                @views g1 = G[r1, :]
                @views g2 = G[r2, :]
                @views g3 = G[r3, :]
                @views if !iszero(sum([g1[1, i] * g2[1, i] * g3[1, i] for i in 1:nc]) % 2)
                    return false
                end
            end
        end
    end
    return true
end

max_stars_repo_name: esabo/CodingTheory

max_issues_repo_name: esabo/CodingTheory

max_forks_repo_name: esabo/CodingTheory


text:
@doc doc"""
Convenience methods for accessing CLV matrices.

Methods `CLV.C`, `CLV.R` and `CLV.G` are the accessor to the matrices
``C``, ``R`` and ``G`` in

```math
G_{n+k} C_{n+k} D_{k,n} = M_{k,n} G_n C_n = G_{n+k} R_{k,n} C_n
```

(Eq. 32, Ginelli et al., 2013).
"""
module CLV

"""``M_n`` (the cocycle)"""
function M end

"""``G_n`` (`Q` from the QR decomposition of ``M_{k,n-k} G_{n-k}``)"""
function G end

"""``R_n`` (`R` from the QR decomposition of ``M_{k,n-k} G_{n-k}``)"""
function R end

"""``R_{n-k}``"""
function R_prev end

"""``C_n``"""
function C end

"""``D_n``"""
function D end

using ..CovariantVectors: ForwardDynamics, BackwardRelaxer, BackwardDynamics,
    BackwardDynamicsWithD, BackwardPass
import ...LyapunovExponents: tangent_propagate

R_prev(fitr::ForwardDynamics) = fitr.le_solver.R
R_prev(bitr::BackwardPass) = bitr.R_history[end-bitr.i-1]
R(bitr::BackwardPass)      = bitr.R_history[end-bitr.i]
R_next(bitr::BackwardPass) = bitr.R_history[end-bitr.i+1]

const R₋ = R_prev
const R₊ = R_next

G(fitr::ForwardDynamics) = fitr.le_solver.tangent_state
C(bitr::Union{BackwardRelaxer, BackwardDynamics}) = bitr.C

D(bitr::BackwardDynamicsWithD) = Diagonal(bitr.D_diag)
M(fitr::ForwardDynamics) = tangent_propagate(fitr.le_solver, I)

end


text:
"""
2-D Gradient

    gradient2(a::Array, dx::Float64=1, dy::Float=1)
    :param a: matrix
    :param dx: stepsize
    :param dy: stepsize

    :return dxdu: derivatives along x
    :return dydv: derivatives along y
"""
function gradient2(a::Array, dx::Float64=1., dy::Float64=1.)
    m,n = size(a);
    dxdu = zeros(m,n);
    dydv = zeros(m,n);

    for i=1:m
        dxdu[i,:] = gradient(vec(a[i,:]), dx)
    end

    for i=1:n
        dydv[:,i] = gradient(a[:,i], dy)
    end

    return dxdu, dydv
end

"""
Creates Rectangular Grid in 2-D space

    meshgrid(a::LinSpace,b::LinSpace)
"""
function meshgrid(a::LinSpace,b::LinSpace)
    grid_a = [i for i in a, j in b]';
    grid_b = [j for i in a, j in b]';

    return grid_a, grid_b
end

"""
Linear interpolation

    approx(xd, yd, xi)
    :param xd: x samples
    :param yd: response samples
    :param xi: new x samples
"""
function approx(xd, yd, xi)
    nd = length(xd);
    ni = length(xi);

    yi = zeros(ni);
    for i in 1:ni
        if (xi[i] <= xd[1])
            t = (xi[i]-xd[1]) / (xd[2] - xd[1]);
            yi[i] = (1.0 - t) * yd[1] + t * yd[2];
        elseif (xd[nd] <= xi[i])
            t = (xi[i] - xd[nd-1]) / (xd[nd] - xd[nd-1]);
            yi[i] = (1.0 - t) * yd[nd-1] + t * yd[nd];
        else
            for k in 2:nd
                if (xd[k-1] <= xi[i] && xi[i] <= xd[k])
                    t = (xi[i] - xd[k-1]) / (xd[k] - xd[k-1]);
                    yi[i] = (1.0 - t) * yd[k-1] + t * yd[k];
                    break
                end
            end
        end
    end

    return yi
end

"""
Trapezoidal Integration

    trapz(x, y, dim=1)
    :param x: vector of time samples
    :param y: array of response samples
    :param dim: dimension along which to integrate
"""
function trapz(x::Array{Float64, 1}, y::Array{Float64}, dim::Integer=1)
    perm = [dim:max(ndims(y),dim); 1:dim-1];
    y = permutedims(y, perm);
    if ndims(y) == 1
        m = 1;
    else
        m = size(y,1);
    end

    if m == 1
        M = length(y);
        out = sum(diff(x).*(y[1:M-1] + y[2:M])/2.);
    else
        out = diff(x).' * (y[1:m-1,:] + y[2:m,:])/2.;
        siz = size(y); siz = collect(siz); siz[1] = 1;
        out = reshape(out, tuple(siz...));
        out = ipermutedims(out, perm);
        ind = find(collect(size(out)).==1);
        out = squeeze(out,ind[1]);
        if length(out) == 1;
            out = out[1];
        end
    end

    return out
end

"""
Cumulative Trapezoidal Integration

    cumtrapz(x, y, dim=1)
    :param x: vector describing time samples
    :param y: array describing response
    :param dim: dimension to integrate over
"""
function cumtrapz(x::Array{Float64, 1}, y::Array{Float64}, dim::Integer=1)
    perm = [dim:max(length(size(y)),dim); 1:dim-1];
    y = permutedims(y, perm);
    if ndims(y) == 1
        n = 1;
        m = length(y);
    else
        m, n = size(y);
    end

    if n == 1
        dt = diff(x)/2.0;
        z = [0; cumsum(dt.*(y[1:(m-1)] + y[2:m]))];
    else
        dt = repmat(diff(x)/2.0,1,n);
        z = [zeros(1,n); cumsum(dt.*(y[1:(m-1), :] + y[2:m, :]),1)];
        z = ipermutedims(z, perm);
    end

    return z
end

"""
Cumulative Trapezoidal Integration using midpoint

    cumtrapzmid(x, y, c)
    :param x: time samples
    :param y: resposne samples
    :param c: midpoint
"""
function cumtrapzmid(x, y, c)
    a = length(x);
    mid = round(Integer, a/2);

    # case < mid
    fn = zeros(a);
    tmpx = x[(mid-1):-1:1];
    tmpy = y[(mid-1):-1:1];
    tmp = c + cumtrapz(tmpx, tmpy);
    fn[1:(mid-1)] = reverse(tmp);

    # case >= mid
    fn[mid:a] = c + cumtrapz(x[mid:a],y[mid:a]);

    return fn

end

"""
Multivariate Normal random number generation

    mvnrand(mu, C, n)
    :param mu: mean vector
    :param C: covariance matrix
    :param n: number of samples
"""
function mvnrand(mu, C, n)
    tmp = cholfact(C, :U, Val{true});
    R = tmp[:U];
    R = Array(R);
    R = R[:, tmp.piv];
    retval = randn(n, size(R,1)) * R;
    retval += transpose(repmat(mu, 1, n));
    return transpose(retval)
end

"""
Linear interpolation when response contains flat regions

    interp1_flat(x, y, xx)
    :param x: time samples
    :param y: response samples
    :param xx: new time samples
"""
function interp1_flat(x,y,xx)
    flat = find(diff(x).<=0);
    n = length(flat);

    if n==0
        tmp = interpolate((x,), y, Gridded(Linear()))
        yy = tmp[xx];
    else
        yy = zeros(size(xx));
        i1 = 1;
        if flat[1] == 1
            i2 = 1;
            j = xx.==x[i2];
            yy[j] = minimum(y[i2:i2+1]);
        else
            i2 = flat[1];
            j = (xx.>=x[i1]) & (xx.<=x[i2]);
            tmp = interpolate((x[i1:i2],), y[i1:i2], Gridded(Linear()))
            yy[j] = tmp[xx[j]];
            i1 = copy(i2);
        end
        for k = 2:n
            i2 = flat[k];
            if i2 > i1+1
                j = (xx.>=x[i1]) & (xx.<=x[i2]);
                yi = interpolate((x[i1+1:i2],), y[i1+1:i2], Gridded(Linear()))
                yy[j] = tmp[xx[j]];
            end
            j = xx.==x[i2];
            yy[j] = minimum(y[i2:i2+1]);
            i1 = copy(i2);
        end
        i2 = length(x);
        j = (xx.>=x[i1]) & (xx.<=x[i2]);
        if i1+1 == i2
            yy[j] = y[i2];
        else
            tmp = interpolate((x[i1+1:i2],), y[i1+1:i2], Gridded(Linear()))
            yy[j] = tmp[xx[j]];
        end
    end

    return yy
end


text:
using Test,Printf,LinearAlgebra

import Distributions
using MatrixCompletion

# import Random,Distributions
# import MatrixCompletion.Losses
# import MatrixCompletion.Losses:train,SGD,train_logistic
# import MatrixCompletion.Convex.ADMM:complete
# import MatrixCompletion.Utilities.RandomMatrices:rand
# import MatrixCompletion.Utilities.Sampling:sample,BernoulliModel
# import MatrixCompletion.Utilities.Indexing:construct_index_tracker,construct_type_matrix
# import MatrixCompletion.Utilities.Indexing:Bernoulli,Gaussian,Poisson,Gamma,NegativeBinomial,Missing,DIST_FLAGS
# import MatrixCompletion.Convex.ADMM:MatrixCompletionModel,
#        predict,
#        _get_column_distribution_type,
#        _get_complete_type_matrix,
#        _predict_bernoulli,
#        _predict_poisson,
#        _predict_gamma,
#        _predict_negative_binomial,
# _predict_gaussian

# using MatrixCompletion.Utilities

# function test_train_logistic(;size = 1000,ρ = 0.1,γ = 0.2,maxIter = 20)

#     y = Int.(Random.bitrand(size));
#     mle_x = train(Losses.Logistic(), Random.rand(size), y, zeros(size), ρ, iter = maxIter, γ = γ);
#     @test sum((Int.(sign.(mle_x)) .+ 1) / 2 .== y)
#     size > 0.99
# end

# test_train_logistic()           #

# function test_train_logistic_optimized(f;size = 1000,ρ = 0.1,γ = 0.2,maxIter = 20)
#     y = Int.(Random.bitrand(size));
#     mle_x = train_logistic(Random.rand(size), y, zeros(size), ρ, iter = maxIter, γ = γ);
#     @test sum((Int.(sign.(mle_x)) .+ 1) / 2 .== y) / size > 0.99
# end

# function test_train_logistic_sgd(size = 3000 * 9000, ρ = 0.1, γ = 0.2, ep = 0)
#     y = Int.(Random.bitrand(size));
#     mle_x = train(SGD(), Losses.Logistic(), Random.rand(size), y, zeros(size), ρ, epoch = ep, γ = γ, num_of_batch = 20);
#     print(sum((Int.(sign.(mle_x)) .+ 1) / 2 .== y) / size);
#     @test sum((Int.(sign.(mle_x)) .+ 1) / 2 .== y) / size > 0.99
# end

# function unit_test_train_poisson(;sz = 500,ρ = 0,step_size = 0.1,maxIter = 100)
#     # sz = 1000
#     # max_iter = 200
#     # γ = 0.1
#     # ρ = 0
#     y = Random.rand(Distributions.Poisson(10), sz) * 1.0
#     mle_x = train(Losses.Poisson(), Random.rand(sz), y, zeros(sz), ρ, iter = max_iter, γ = step_size);
#     recoveredX = round.(exp.(mle_x));
#     errRate = sum(abs.(recoveredX .- y) .> 1) / sz;
# #    @test 1 - errRate >= 0.99
#     return 1- errRate;
# end

# function POISSON_SMALL_TEST_SET_LOOSE()
#     @test unit_test_train_poisson(sz=1000,ρ=0,step_size=0.1,maxIter=200) > 0.9
#     @test unit_test_train_poisson(sz=3000,ρ=0,step_size=0.1,maxIter=200) > 0.9
#     @test unit_test_train_poisson(sz=5000,ρ=0,step_size=0.1,maxIter=200) > 0.9
#     @test unit_test_train_poisson(sz=1000,ρ=0.1,step_size=0.1,maxIter=200) > 0.9
#     @test unit_test_train_poisson(sz=3000,ρ=0.1,step_size=0.1,maxIter=200) > 0.9
#     @test unit_test_train_poisson(sz=5000,ρ=0.1,step_size=0.1,maxIter=200) > 0.9
# end

# function test_train_gamma(;size = 500,ρ = 0.05,γ = 0.1,maxIter = 100)
# end

# function test_train_negative_binomial(;size = 500,ρ = 0.05,γ = 0.1,maxIter = 100)
# end

# function test_complete_type_matrix()
#     test1_input = Array{Union{Missing,DIST_FLAGS},2}(undef, 10, 10)
#     test1_input[:,1:2] .= Ref(Gaussian);
#     test1_input[:,3:4] .= Ref(Bernoulli);
#     test1_input[:,5:6] .= Ref(Gamma);
#     test1_input[:,7:8] .= Ref(Poisson);
#     test1_input[:,9:10] .= Ref(NegativeBinomial);
#     test1_expect = deepcopy(test1_input);
#     test1_input[diagind(test1_input)] .= missing;
#     test1_output = _get_complete_type_matrix(test1_input);
#     @test test1_output == test1_expect;
#     test2_input = Array{DIST_FLAGS,2}(undef, 10, 10)
#     test2_input[:,1:2] .= Ref(Gaussian);
#     test2_input[:,3:4] .= Ref(Bernoulli);
#     test2_input[:,5:6] .= Ref(Gamma);
#     test2_input[:,7:8] .= Ref(Poisson);
#     test2_input[:,9:10] .= Ref(NegativeBinomial);
#     test2_expect = deepcopy(test2_input);
#     test2_output = _get_complete_type_matrix(test2_input);
#     @test test2_output == test2_expect;
# end

# function test_predict_bernoulli()
#     size = 500;ρ = 0.1;γ = 0.1;maxIter = 500;
#     y = Int.(Random.bitrand(size));
#     mle_x = train(Losses.Logistic(), Random.rand(size), y, zeros(size), ρ, iter = maxIter, γ = γ);
#     @test sum(_predict_bernoulli(mle_x) .== y) / size > 0.99
# end

# test_predict_bernoulli()

using MatrixCompletion

# function unit_test_admm(;input_matrix,                  = rand(Distributions.Gaussian(0,1),500,500,3),
#                         distribution_type_matrix        = provide()
#                         sampling_model                  = BernoulliModel(),
#                         sampling_rate                   = 0.8,
#                         max_iter_inner_gradient_descent = 3,
#                         max_iter_admm                   = 200,
#                         stop_tol_admm                   = 1e-5,
#                         debug_mode                      = false,
#                         use_auto_diff                   = true,
#                         λ                               = 5e-1,
#                         μ                               = 5e-4,
#                         σ                               = 0.3,
#                         τ                               = 1.618)

# end

function accuracyImputedBinaryPart(;truth::Array{Float64,2} = nothing, completedMatrix::Array{Float64,2} = nothing)
    typeM = construct_type_matrix(truth);
#     binaryColumns = find(x->x==BINARY,typeM[1,:]);
    binaryColumns = findall(x->x == Bernoulli, typeM[1,:]);
    imputedBinaryPart  = (sign.(completedMatrix[:,binaryColumns]) .+ 1) / 2
    return sum(Int.(truth[:,binaryColumns] .== imputedBinaryPart)) / (length(binaryColumns) * size(truth)[1]);
end

function accuracyImputedContinuousPart(;truth::Array{Float64,2} = nothing,completedMatrix::Array{Float64,2} = nothing)
    typeM = construct_type_matrix(truth);
    continuousColumns = findall(x->x == Gaussian, typeM[1,:])
    imputedContinuousPart = completedMatrix[:,continuousColumns];
    imputedContinuousPart - truth[:,continuousColumns]
    return norm(imputedContinuousPart - truth[:,continuousColumns])^2 / norm(truth[:,continuousColumns])^2
end

function test_admm_with_autodiff_smallinput(;gd_iter = 3,dbg = false)
    admm_test_matrix1 = rand([(Distributions.Bernoulli(0.7), 100 => 50, 3),(Distributions.Gaussian(3, 1), 100 => 50, 3)])
    admm_test_matrix_missing1 = sample(BernoulliModel(), x = admm_test_matrix1, rate = 0.8)
    @time admm_test_matrix_output_1 = complete(A = admm_test_matrix_missing1, maxiter = 200, use_autodiff = true, gd_iter = gd_iter, debug_mode = dbg);
    gaussian_acc = accuracyImputedContinuousPart(truth = admm_test_matrix1, completedMatrix = admm_test_matrix_output_1)
    bernoulli_acc = accuracyImputedBinaryPart(truth = admm_test_matrix1, completedMatrix = admm_test_matrix_output_1)
    @printf("gaussian acc: %f\n", gaussian_acc)
    @printf("bernoulli acc: %f\n",bernoulli_acc)
end

function test_admm_without_autodiff_smallinput(;gd_iter = 3,dbg = false)
    admm_test_matrix1 = rand([(Distributions.Bernoulli(0.7), 100 => 50, 3),(Distributions.Gaussian(3, 1), 100 => 50, 3)])
    admm_test_matrix_missing1 = sample(BernoulliModel(), x = admm_test_matrix1, rate = 0.8)
    @time admm_test_matrix_output_1 = complete(A = admm_test_matrix_missing1, maxiter = 200, use_autodiff = false, gd_iter = gd_iter, debug_mode = dbg)
    gaussian_acc = accuracyImputedContinuousPart(truth = admm_test_matrix1, completedMatrix = admm_test_matrix_output_1)
    bernoulli_acc = accuracyImputedBinaryPart(truth = admm_test_matrix1, completedMatrix = admm_test_matrix_output_1)
    @printf("gaussian acc: %f\n", gaussian_acc)
    @printf("bernoulli acc: %f\n",bernoulli_acc)
end

function test_admm_without_autodiff_largeinput(;gd_iter = 3,dbg = false)
    admm_test_matrix1 = rand([(Distributions.Bernoulli(0.7), 6000 => 3000, 10),(Distributions.Gaussian(3, 1), 6000 => 3000, 10)])
    admm_test_matrix_missing1 = sample(BernoulliModel(), x = admm_test_matrix1, rate = 0.8)
    @time admm_test_matrix_output_1 = complete(A = admm_test_matrix_missing1, maxiter = 200, use_autodiff = false, gd_iter = gd_iter, debug_mode = dbg)
    gaussian_acc = accuracyImputedContinuousPart(truth = admm_test_matrix1, completedMatrix = admm_test_matrix_output_1)
    bernoulli_acc = accuracyImputedBinaryPart(truth = admm_test_matrix1, completedMatrix = admm_test_matrix_output_1)
end

#
#
#
# admm_test_matrix2 = rand([(Distributions.Bernoulli(0.7),5000=>2500,100),(Distributions.Gaussian(3,1),5000=>2500,100)])
# admm_test_matrix_missing2 = sample(BernoulliModel(),x = admm_test_matrix2,rate = 0.8)
#
# a = construct_type_matrix(admm_test_matrix2)
# construct_index_trakcer(input_type_matrix=a)
# admm_test_matrix_output_2 = complete(A = admm_test_matrix_missing2)
# accuracyImputedContinuousPart(truth=admm_test_matrix2,completedMatrix = admm_test_matrix_output_2)
# accuracyImputedBinaryPart(truth=admm_test_matrix2,completedMatrix = admm_test_matrix_output_2)
#
#
#
# #
# # _get_column_distribution_type(type_mat[:,1])
# # _get_column_distribution_type(type_mat[:,51])
# # type_mat[:,90] .= Ref(_get_column_distribution_type(type_mat[:,51]))
#
#
# complete_type_mat = _get_complete_type_matrix(type_mat)
# complete_type_mat[:,90]
# predict(x=ot,obs=type_mat)
#
#
# MatrixCompletionModel(observed=admm_test_matrix_missing1,completed=ot,type_matrix=type)
#
#
# test_train_logistic()
# test_train_poisson()
# test_complete_type_matrix()


text:
"""
   find_worst_telluric_in_each_chunk( chunk_list_timeseries, array_of_spectra )
Returns a vector of the worst telluric (i.e., smallest value) within each chunk at any observations.

Warning: Assumes a tellurics value in metadata for each spectra, such as is provided by EXPRES.
"""
function find_worst_telluric_in_each_chunk end

"""
   `make_clean_line_list_from_tellurics_expres(line_list, expres_data; Δv_to_avoid_tellurics )`
Returns a new line list that excludes lines with telluric contamination.
Inputs:
- line_list:  Dataframe containing field lambda
- expres_data:  Array of spectra
- Δv_to_avoid_tellurics:  in m/s
Outputs:
- line_list_without_tellurics:   DataFrame with fields, lambda, weight, lambda_lo, and lambda_hi.
Warning: Currently, assumes a tellurics value in metadata for each spectra, such as is provided by EXPRES.
"""
function make_clean_line_list_from_tellurics_expres(line_list::DataFrame, expres_data::DT; Δv_to_avoid_tellurics::Real = default_Δv_to_avoid_tellurics,
            v_center_to_avoid_tellurics::Real = 0.0
               ) where { T1<:Real, A1<:AbstractArray{T1}, T2<:Real, A2<:AbstractArray{T2}, T3<:Real, A3<:AbstractArray{T3}, IT<:EXPRES.AnyEXPRES, ST<:Spectra2DBasic{T1,T2,T3,A1,A2,A3,IT}, DT<:AbstractArray{ST,1} }
   @assert 1000 <= Δv_to_avoid_tellurics <= 50000
   line_list_to_search_for_tellurics = copy(line_list)
   # TODO: Check sign convention for v_center_to_avoid_tellurics
   line_list_to_search_for_tellurics.lambda_lo = line_list_to_search_for_tellurics.lambda./calc_doppler_factor(Δv_to_avoid_tellurics).*calc_doppler_factor(v_center_to_avoid_tellurics)
   line_list_to_search_for_tellurics.lambda_hi = line_list_to_search_for_tellurics.lambda.*calc_doppler_factor(Δv_to_avoid_tellurics).*calc_doppler_factor(v_center_to_avoid_tellurics)
   chunk_list_timeseries = RvSpectML.make_chunk_list_timeseries(expres_data,line_list_to_search_for_tellurics)
   line_list_to_search_for_tellurics.min_telluric_model_all_obs = EXPRES.find_worst_telluric_in_each_chunk( chunk_list_timeseries, expres_data)
   line_list_no_tellurics_df = line_list_to_search_for_tellurics |> @filter(_.min_telluric_model_all_obs == 1.0) |> DataFrame
end

max_stars_repo_name: RvSpectML/EchelleCCFs

max_issues_repo_name: RvSpectML/EchelleCCFs

max_forks_repo_name: RvSpectML/EchelleCCFs


text:
export project_bounds!

function project_bounds!(x::Vector{TF},LB::TF,UB::TF) where {TF<:Real}
  """
  Computes the projection of x onto the set of constraints LB <= x <= UB
  x is a vector but the bounds are scalars. Uses some manual tricks like @inbounds and multi threading, attempt to be faster than native julia: x. = max.(LB,min.(x,UB))
  """
  Threads.@threads for j in (1:length(x))
    @inbounds x[j] = max(LB,min(x[j],UB))
  end
    return x
end

function project_bounds!(x::Vector{TF},LB::Vector{TF},UB::Vector{TF}) where {TF<:Real}
  """
  Computes the projection of x onto the set of constraints LB <= x <= UB
  x is a vector and the bounds are vectors as well. Uses some manual tricks like @inbounds and multi threading, attempt to be faster than native julia: x. = max.(LB,min.(x,UB))
  """
@inbounds Threads.@threads for j in (1:length(x))
    #x[i] = max.(LB,min.(x[i],UB))
  @inbounds  x[j] = min(x[j],UB[j])
  @inbounds  x[j] = max(LB[j],x[j])
  end
    return x
end

function project_bounds!(x::Array{TF,2},LB::Vector{TF},UB::Vector{TF},mode::Tuple{String,String}) where {TF<:Real}
  """
  Computes the projection of x onto the set of constraints LB <= x <= UB
  x is a matrix and the bounds are per row or column.
  """

  if mode[2] == "x"
    Threads.@threads for i=1:size(x,2)
      @inbounds x[:,i].=min.(max.(x[:,i],LB),UB)
    end
  elseif mode[2] == "z"
    Threads.@threads for i=1:size(x,1)
      @inbounds x[i,:].=min.(max.(x[i,:],LB),UB)
    end
  end
  x = vec(x)
return x
end

function project_bounds!(x::Array{TF,3},LB::Vector{TF},UB::Vector{TF},mode::Tuple{String,String}) where {TF<:Real}
  """
  Computes the projection of x onto the set of constraints LB <= x <= UB
  x is a 3D array and the bounds are per fiber (x/y/z).
  """

if mode[1] == "fiber"
  if mode[2] == "x"
    Threads.@threads for i=1:size(x,2)
      for j=1:size(x,3)
        @inbounds x[:,i,j].=min.(max.(x[:,i,j],LB),UB)
      end
    end
  elseif mode[2] == "y"
    Threads.@threads for i=1:size(x,1)
      for j=1:size(x,3)
        @inbounds x[i,:,j].=min.(max.(x[i,:,j],LB),UB)
      end
    end
  elseif mode[2] == "z"
    Threads.@threads for i=1:size(x,1)
      for j=1:size(x,2)
        @inbounds x[i,j,:].=min.(max.(x[i,j,:],LB),UB)
      end
    end
  end
elseif mode[1] == "slice"
    error("bound constraints per slice of a tensor currently not implemented, yet...")
end
  x = vec(x)
return x
end

max_stars_repo_name: slimgroup/SetIntersectionProjection

max_issues_repo_name: slimgroup/SetIntersectionProjection

max_forks_repo_name: slimgroup/SetIntersectionProjection