text:
# Tudor Berariu, 2015

module Miners

using Scenarios;

export Gatherers;

# -----------------------------------------------------------------------------

#= Global constants =#

const HEIGHT           = 12;
const WIDTH            = 12;
const AGENTS_NO        = 8;
const WAREHOUSES_NO    = 1;
const MINES_NO         = 1;
const DISTANCE         = 6;
const RANGE            = 15;
const AGENT_RANGE      = 2.5;

# -----------------------------------------------------------------------------

#= Type aliases =#

typealias CellInfo    Uint8;
typealias Signal      Uint8;
typealias AgentState  Array{Uint8, 1};
typealias Action      Uint8;

# -----------------------------------------------------------------------------

#= CellInfo =#

#=
    +---+---+---+---+---+---+---+---+
    |   |   |   |   |   |   |   |N/S|
    +---+---+---+---+---+---+---+---+
=#
const NOTHING   = convert(CellInfo, 0);
const SOMETHING = convert(CellInfo, 1);

isNothing(cell::CellInfo)   = cell == 0;
isSomething(cell::CellInfo) = cell > 0;

#=
    Entity info
    +---+---+---+---+---+---+---+---+
    |   |   |   |   | W | B | A | 1 |
    +---+---+---+---+---+---+---+---+
=#
const AGENT    = SOMETHING | convert(CellInfo, 1<<1);
const BUILDING = SOMETHING | convert(CellInfo, 1<<2);
const WALL     = SOMETHING | convert(CellInfo, 1<<3);

isAgent(cell::CellInfo)    = (cell & AGENT) == AGENT;
isBuilding(cell::CellInfo) = (cell & BUILDING) == BUILDING;
isWall(cell::CellInfo)     = (cell & WALL) == WALL;

#=
    Agent info
    +---+---+---+---+---+---+---+---+
    |   |   | E | L | 0 | 0 | 1 | 1 |
    +---+---+---+---+---+---+---+---+
=#

const LOADED = SOMETHING | AGENT | convert(CellInfo, 1<<4)
const EMPTY  = SOMETHING | AGENT | convert(CellInfo, 1<<5)

isLoaded(agent::CellInfo) = (agent & LOADED) == LOADED;
isEmpty(agent::CellInfo)  = (agent & EMPTY) == EMPTY;

#=
    Building info
    +---+---+---+---+---+---+---+---+
    |   |   | M | W | 0 | 1 | 0 | 1 |
    +---+---+---+---+---+---+---+---+
=#

const WAREHOUSE = SOMETHING | BUILDING | convert(CellInfo, 1<<4);
const MINE      = SOMETHING | BUILDING | convert(CellInfo, 1<<5);

isWarehouse(building::CellInfo) = (building & WAREHOUSE) == WAREHOUSE;
isMine(building::CellInfo)      = (building & MINE     ) == MINE;

const CellInfoToChar = Dict{CellInfo, Char}(EMPTY     => char('\u26c4'),
                                            LOADED    => char('\u26c7'),
                                            WAREHOUSE => char('\u26ea'),
                                            MINE      => char('\u233e'),
                                            NOTHING   => ' ',
                                            WALL      => char('\u2588'));

# -----------------------------------------------------------------------------

#=
    Actions
=#

const DO_NOTHING     = convert(Action, 0) ;              # 00000000

const DO_MOVE        = convert(Action, 1<<0);            # ----0001
const DO_LOAD        = convert(Action, 1<<1);            # 00000010
const DO_UNLOAD      = convert(Action, 1<<2);            # 00000100
const DO_RANDOM_MOVE = convert(Action, 1<<3);            # 00001000

const UP             = DO_MOVE | convert(Action, 1<<4);  # 00010001
const DOWN           = DO_MOVE | convert(Action, 1<<5);  # 00100001
const LEFT           = DO_MOVE | convert(Action, 1<<6);  # 01000001
const RIGHT          = DO_MOVE | convert(Action, 1<<7);  # 10000001

const UP_LEFT        = UP   | LEFT ;                     # 01010001
const UP_RIGHT       = UP   | RIGHT;                     # 10010001
const DOWN_LEFT      = DOWN | LEFT ;                     # 01100001
const DOWN_RIGHT     = DOWN | RIGHT;                     # 10100001

const moveActions = [UP, DOWN, LEFT, RIGHT,
                     UP_LEFT, UP_RIGHT, DOWN_LEFT, DOWN_RIGHT];

const validActions = [DO_NOTHING,
                      DO_LOAD, DO_UNLOAD, #DO_RANDOM_MOVE,
                      UP, DOWN, LEFT, RIGHT,
                      UP_LEFT, UP_RIGHT, DOWN_LEFT, DOWN_RIGHT];

isMoveAction(a::Action) = (a & DO_MOVE) == DO_MOVE;

δy(a::Action) = (a & UP)   > DO_MOVE ? -1 : (a & DOWN)  > DO_MOVE ? 1 : 0;
δx(a::Action) = (a & LEFT) > DO_MOVE ? -1 : (a & RIGHT) > DO_MOVE ? 1 : 0;

nextCell(row::Integer, column::Integer, action::Action) =
    min(HEIGHT, max(1, row    + δy(action))),
    min(WIDTH,  max(1, column + δx(action)));

const ActionToChar = Dict{Action, Char}(DO_NOTHING     => char('\u2716'),
                                        DO_LOAD        => char('\u27f0'),
                                        DO_UNLOAD      => char('\u27f1'),
                                        DO_RANDOM_MOVE => char('?'),
                                        UP             => char('\u21d1'),
                                        DOWN           => char('\u21d3'),
                                        LEFT           => char('\u21d0'),
                                        RIGHT          => char('\u21d2'),
                                        UP_LEFT        => char('\u21d6'),
                                        UP_RIGHT       => char('\u21d7'),
                                        DOWN_LEFT      => char('\u21d9'),
                                        DOWN_RIGHT     => char('\u21d8'))

# -----------------------------------------------------------------------------

#= Signal =#

#=
    +---+---+---+---+---+---+---+---+
    | MINE SIGNAL   | DEPOSIT SIGNAL|
    +---+---+---+---+---+---+---+---+
=#

const NO_SIGNAL = convert(Signal, 0);
const NORTH     = convert(Signal, 1<<0);
const SOUTH     = convert(Signal, 1<<1);
const EAST      = convert(Signal, 1<<2);
const WEST      = convert(Signal, 1<<3);

const R4      = convert(Signal, 1 | 1<<1 | 1<<2 | 1<<3);

signalInfo(s::Signal, oldS::Signal) =
    ( s== NO_SIGNAL) ? (oldS & R4) : ((s<<4) | s);

const SignalToChar = Dict{Signal, Char}(NO_SIGNAL    => char(' '),
                                        WEST         => char('\u2190'),
                                        NORTH        => char('\u2191'),
                                        EAST         => char('\u2192'),
                                        SOUTH        => char('\u2193'),
                                        NORTH | WEST => char('\u2196'),
                                        NORTH | EAST => char('\u2197'),
                                        SOUTH | WEST => char('\u2199'),
                                        SOUTH | EAST => char('\u2198'));

# -----------------------------------------------------------------------------

#= Cartesian stuff =#

const ROW    = 1;
const COLUMN = 2;
const RC = [ROW, COLUMN];

sqEuclid(r1, c1, r2, c2) = (r1 - r2) * (r1 - r2) + (c1 - c2) * (c1 - c2);

# ---------------------------------------

type GlobalState
    board::Array{CellInfo, 2}
    warehouses::Array{Int64, 2}
    mines::Array{Int64, 2}
    agents::Array{Int64, 2}
    wSignal::Array{Signal, 2}
    mSignal::Array{Signal, 2}
    agentsMap::Array{CellInfo, 2}
end

# ---------------------------------------

function placeBuildings()
    assert(WAREHOUSES_NO == MINES_NO);

    warehouses = zeros(Int64, 2, WAREHOUSES_NO);
    mines = zeros(Int64, 2, MINES_NO);

    notOver = true;
    idx = 0;

    const MAX_ATTEMPTS = 10;

    while idx < WAREHOUSES_NO
        attemptNo = 0;
        found = false
        while attemptNo < MAX_ATTEMPTS && found == false
            attemptNo += 1;
            rw = rand(2:HEIGHT-1);
            cw = rand(2:WIDTH-1);
            isGood = true;
            for i = 1:idx
                if maximum(abs(warehouses[RC,i] - [rw, cw])) <=  DISTANCE ||
                    maximum(abs(mines[RC,i] - [rw, cw])) <= DISTANCE
                    isGood = false;
                    break;
                end
            end
            if isGood == false
                continue
            end
            p = filter(crd -> (crd[1] > 1) && (crd[1] < HEIGHT) &&
                       (crd[2] > 1) && (crd[2] < WIDTH) &&
                       (max(abs(crd[1]-rw), abs(crd[2]-cw)) == DISTANCE),
                       [(r, c) for r in rw-DISTANCE:rw+DISTANCE,
                        c in cw-DISTANCE:cw+DISTANCE]);
            if length(p) == 0
                continue
            end
            crd = rand(p);
            rm = crd[1];
            cm = crd[2];
            for i = 1:idx
                if maximum(abs(warehouses[RC,i] - [rm, cm])) <=  DISTANCE ||
                    maximum(abs(mines[RC,i] - [rm, cm])) <= DISTANCE
                    isGood = false;
                    break;
                end
            end
            if isGood == false
                continue;
            end
            found = true
            idx = idx + 1;
            warehouses[RC, idx] = [rw, cw];
            mines[RC, idx] = [rm, cm];
        end
        if found == false
            warehouses = zeros(Int64, 2, WAREHOUSES_NO);
            mines = zeros(Int64, 2, MINES_NO);
            idx = 0;
        end
    end

    return warehouses, mines;
end

function computeSignal(board::Array{CellInfo,2},
                       warehouses::Array{Int64,2}, mines::Array{Int64,2})
    mSignal = fill(NO_SIGNAL, size(board));
    wSignal = fill(NO_SIGNAL, size(board));
    for r in 1:HEIGHT
        for c in 1:WIDTH
            #-- Warehouses
            bestDistance = Inf
            wIdx = 0
            for w in 1:WAREHOUSES_NO
                d = sqEuclid(warehouses[ROW, w], warehouses[COLUMN, w], r, c);
                if d < bestDistance
                    bestDistance = d;
                    wIdx = w;
                end
            end
            if bestDistance <= RANGE * RANGE
                if r < warehouses[ROW, wIdx]
                    wSignal[r,c] = wSignal[r,c] | SOUTH;
                elseif r > warehouses[ROW, wIdx]
                    wSignal[r,c] = wSignal[r,c] | NORTH;
                end
                if c < warehouses[COLUMN, wIdx]
                    wSignal[r,c] = wSignal[r,c] | EAST;
                elseif c > warehouses[COLUMN, wIdx]
                    wSignal[r,c] = wSignal[r,c] | WEST;
                end
            end
            #-- Mines
            bestDistance = Inf;
            mIdx = 0;
            for m in 1:MINES_NO
                d = sqEuclid(mines[ROW, m], mines[COLUMN, m], r, c);
                if d < bestDistance
                    bestDistance = d;
                    mIdx = m;
                end
            end
            if bestDistance <= RANGE * RANGE
                if r < mines[ROW, mIdx]
                    mSignal[r,c] = mSignal[r,c] | SOUTH;
                elseif r > mines[ROW, mIdx]
                    mSignal[r,c] = mSignal[r,c] | NORTH;
                end
                if c < mines[COLUMN, mIdx]
                    mSignal[r,c] = mSignal[r,c] | EAST;
                elseif c > mines[COLUMN, mIdx]
                    mSignal[r,c] = mSignal[r,c] | WEST;
                end
            end
        end
    end
    return wSignal, mSignal;
end

function initialState()
    board = fill(NOTHING, HEIGHT, WIDTH);
    # Put warehouses and mines
    warehouses, mines = placeBuildings();
    for w in 1:WAREHOUSES_NO
        board[warehouses[ROW, w], warehouses[COLUMN, w]] = WAREHOUSE;
    end
    for m in 1:MINES_NO
        board[mines[ROW, m], mines[COLUMN, m]] = MINE;
    end
    # Put agents
    agents = zeros(Int64, 2, AGENTS_NO);
    agentsMap = zeros(Int64, HEIGHT, WIDTH);
    for a in 1:AGENTS_NO
        agents[ROW, a]    = rand(1:HEIGHT)
        agents[COLUMN, a] = rand(1:WIDTH)
        while (board[agents[ROW, a], agents[COLUMN, a]] != NOTHING)
            agents[ROW, a]    = rand(1:HEIGHT)
            agents[COLUMN, a] = rand(1:WIDTH)
        end
        board[agents[ROW, a], agents[COLUMN, a]] = EMPTY;
        agentsMap[agents[ROW, a], agents[COLUMN, a]] = a;
    end
    # Compute signals
    wSignal, mSignal = computeSignal(board, warehouses, mines);
    # Pack everything
    return GlobalState(board, warehouses, mines, agents, wSignal, mSignal,
                       agentsMap);
end

function perceiveMap(gs::GlobalState,
                     a::Int64,
                     oldState::AgentState = [NO_SIGNAL, NO_SIGNAL])
    const diameter = trunc(Int64, AGENT_RANGE) * 2 + 1;
    const mid = trunc(Int64, AGENT_RANGE) + 1;
    const row = gs.agents[ROW,a];
    const col = gs.agents[COLUMN,a];

    state = fill(NOTHING, diameter, diameter);

    for c in 1:diameter, r in 1:diameter
        if sqEuclid(r, c, mid, mid) <= AGENT_RANGE * AGENT_RANGE
            state[r, c] =
                ((1 <= col-mid+c <= WIDTH) && (1 <= row-mid+r <= HEIGHT)) ?
            gs.board[row-mid+r,col-mid+c] : WALL
        end
    end

    const wSignal = signalInfo(gs.wSignal[gs.agents[ROW,a],gs.agents[COLUMN,a]],
                               oldState[1]);
    const mSignal = signalInfo(gs.mSignal[gs.agents[ROW,a],gs.agents[COLUMN,a]],
                               oldState[2]);
    return vcat(wSignal, mSignal, state[:]);
end

function neighbours(gs::GlobalState, ag::Int64)
    const row = gs.agents[ROW,ag];
    const col = gs.agents[COLUMN,ag];
    const mid = trunc(Int64, AGENT_RANGE) + 1;
    const diameter = trunc(Int64, AGENT_RANGE) * 2 + 1;

    nbs = Array(Int64, 0);

    for c in 1:diameter, r in 1:diameter
        if sqEuclid(r, c, mid, mid) <= AGENT_RANGE * AGENT_RANGE
            if ((1 <= col-mid+c <= WIDTH) && (1 <= row-mid+r <= HEIGHT))
                cellInfo = gs.agentsMap[row-mid+r,col-mid+c];
                if cellInfo > 0 && cellInfo != ag
                    push!(nbs, gs.agentsMap[row-mid+r, col-mid+c]);
                end
            end
        end
    end
    return nbs;
end

function sortAgents(gs::GlobalState, actions::Array{Action, 1})
    sortedAgents = zeros(Int64, AGENTS_NO);
    crtIdx = 1;

    isBlocked = fill(true, AGENTS_NO);
    waitsOn = zeros(Int64, AGENTS_NO);

    # First put in the finall list all agents that do not wait on others to move
    for i in 1:AGENTS_NO
        if isMoveAction(actions[i])
            const crtRow = gs.agents[ROW,i];
            const crtColumn = gs.agents[COLUMN,i];
            nextRow, nextColumn = nextCell(crtRow, crtColumn, actions[i]);
            if (((nextRow != crtRow) || (nextColumn != crtColumn)) &&
                isAgent(gs.board[nextRow, nextColumn]) &&
                isBlocked[gs.agentsMap[nextRow, nextColumn]])
                waitsOn[i] = gs.agentsMap[nextRow, nextColumn];
                continue;
            end
        end
        sortedAgents[crtIdx] = i;
        crtIdx = crtIdx + 1;
        isBlocked[i] = false;
    end

    # Put agents that wait for agents already put in the final list
    foundOne = true;
    while foundOne && crtIdx <= AGENTS_NO
        foundOne = false;
        for i in find(isBlocked)
            if isBlocked[waitsOn[i]] == false
                foundOne = true
                sortedAgents[crtIdx] = i;
                crtIdx = crtIdx + 1;
                isBlocked[i] = false;
            end
        end
    end
    return vcat(sortedAgents[1:(crtIdx-1)], find(isBlocked));
end

function doLoad!(gs::GlobalState, ag::Int64)
    row = gs.agents[ROW, ag];
    col = gs.agents[COLUMN, ag];
    if isEmpty(gs.board[row, col])
        for m in 1:MINES_NO
            if sqEuclid(gs.mines[ROW,m], gs.mines[COLUMN,m], row, col) <= 2
                gs.board[row,col] = LOADED;
                return
            end
        end
    end
    nothing
end

function doUnload!(gs::GlobalState, ag::Int64)
    row = gs.agents[ROW, ag];
    col = gs.agents[COLUMN, ag];
    if isLoaded(gs.board[row, col])
        for w in 1:WAREHOUSES_NO
            if sqEuclid(gs.warehouses[ROW,w], gs.warehouses[COLUMN,w],
                        row, col) <= 2
                gs.board[row, col] = EMPTY;
                return w;
            end
        end
    end
    return 0;
end

function doRandomMove!(gs::GlobalState, ag::Int64)
    doMove!(gs, ag, moveActions[rand(1:end)])
    nothing
end

function doMove!(gs::GlobalState, ag::Int64, action::Action)
    const row = gs.agents[ROW, ag];
    const column = gs.agents[COLUMN, ag];
    const newRow, newColumn = nextCell(row, column, action);
    if gs.board[newRow, newColumn] == NOTHING
        gs.board[newRow, newColumn] = gs.board[row, column];
        gs.board[row, column] = NOTHING;
        gs.agents[ROW, ag] = newRow;
        gs.agents[COLUMN, ag] = newColumn;
        gs.agentsMap[row, column] = NOTHING;
        gs.agentsMap[newRow, newColumn] = ag;
    end
    nothing
end

reward(dropCount::Int64) = 2.5 ^ dropCount - 1;

function sizeHint()
    x = zero(Int64)
    const mid = trunc(Int64, AGENT_RANGE) + 1;
    const diameter = trunc(Int64, AGENT_RANGE) * 2 + 1;
    const cells = [(sqEuclid(r,c,mid,mid) <= AGENT_RANGE * AGENT_RANGE) ? 1 : 0
                   for r in 1:diameter, c in 1:diameter];
    const cellsNo = sum(cells)-1;
    for i = 1:AGENTS_NO-1
        l = one(Int64);
        for j = 0:i
            l *= 2 * (cellsNo - j);
        end
        x += l;
    end
    x *= 16 * 4 * (cellsNo - AGENTS_NO + 1);
    x += cellsNo * 16 * 4;
    println(x);
    return x;
end

function doActions!(gs::GlobalState, actions::Array{Action, 1})
    rewards        = zeros(Float64, AGENTS_NO);
    rewarders      = zeros(Int64, AGENTS_NO);
    warehouseDrops = zeros(Int64,   WAREHOUSES_NO);

    for ag in sortAgents(gs, actions)
        if actions[ag] != DO_NOTHING
            if actions[ag] == DO_LOAD
                doLoad!(gs, ag);
            elseif actions[ag] == DO_UNLOAD
                wIdx = doUnload!(gs, ag);
                if wIdx > 0
                    warehouseDrops[wIdx] = warehouseDrops[wIdx] + 1;
                    rewarders[ag] = wIdx;
                end
            elseif actions[ag] == DO_RANDOM_MOVE
                doRandomMove!(gs, ag);
            else
                doMove!(gs, ag, actions[ag]);
            end
        end
    end

    for ag in 1:AGENTS_NO
        if rewarders[ag] > 0
            if warehouseDrops[rewarders[ag]] > 0
                rewards[ag] = reward(warehouseDrops[rewarders[ag]]);
            end
        end
    end

    return rewards, warehouseDrops;
end

function printState(state::AgentState)
    wSignal = state[1];
    crtWSignal = wSignal>>4;
    oldWSignal = (wSignal<<4)>>4;
    mSignal = state[2];
    crtMSignal = mSignal>>4;
    oldWSignal = (mSignal<<4)>>4;
    println("W: $(SignalToChar[crtWSignal]) - $(SignalToChar[crtWSignal])");
    println("M: $(SignalToChar[crtMSignal]) - $(SignalToChar[crtMSignal])");
    l = int(sqrt(length(state[3:end])));
    printBoard(reshape(state[3:end], l, l));
end

function printSignalMap(signalMap::Array{Signal, 2})
    for col in 1:(size(signalMap,2)+2)
        print(char('\u2588'));
    end
    println();
    for row in 1:size(signalMap,1)
        print(char('\u2588'));
        for col in 1:size(signalMap,2)
            print(SignalToChar[signalMap[row, col]]);
        end
        println(char('\u2588'));
    end
    for col in 1:(size(signalMap,2)+2)
        print(char('\u2588'));
    end
    println();
end

function printBoard(board::Array{CellInfo, 2})
    for col in 1:(size(board,2)+2)
        print(char('\u2588'));
    end
    println();
    for row in 1:size(board,1)
        print(char('\u2588'));
        for col in 1:size(board,2)
            print(CellInfoToChar[board[row, col]]);
        end
        println(char('\u2588'));
    end
    for col in 1:(size(board,2)+2)
        print(char('\u2588'));
    end
    println();
end

Gatherers = Scenario{AgentState, Action}(AGENTS_NO, WAREHOUSES_NO,
                                         initialState, perceiveMap,
                                         validActions, doActions!,
                                         neighbours);

end

max_stars_repo_name: tudor-berariu/phd

max_issues_repo_name: tudor-berariu/phd

max_forks_repo_name: tudor-berariu/phd


text:
"""
Update inversion variable with the initial model.
At the same time, update the lower and upper bounds on the inversion variable.
"""
function initialize!(pa::Param)

	randn!(pa.mx.last_x)  # reset last_x

	# use mod_initial as starting model according to parameterization
	Models.Seismic_get!(pa.mx.x, pa.mod_initial, pa.parameterization)

	# bounds
	Seismic_xbound!(pa.mx.lower_x, pa.mx.upper_x, pa)

end

# finalize xfwi or wfwi on pa
function finalize!(pa::Param, xminimizer, attrib_mod)
	# update modm
	x_to_modm!(pa, xminimizer)

	# update calculated data using the minimizer
	F!(pa, nothing, attrib_mod)

	# put minimizer into modi
	Models.Seismic_reparameterize!(pa.modi, xminimizer, pa.parameterization)

	# update initial model using minimizer, so that running xfwi! next time will start with updated model
	copyto!(pa.mod_initial, pa.modi)
end

"""
Full Waveform Inversion using `Optim` and `Ipopt` packages.
This method updates `pa.modm` and `pa.dcal`. More details about the
optional parameters can be found in the documentation of the `Optim`
package.
The gradiet is computed using the adjoint state method.
`pa.modi` is used as initial model if non-zero.

# Arguments

* `pa::Param` : inversion object (updated inside method)
* `obj::Union{LS,LS_prior}` : which objective function?
  * `=LS()`
  * `=LS_prior([1.0, 0.5])`

# Optional Arguments

* `optim_options` : see Optim.jl package for optimization options
* `optim_scheme=LBFGS()` : see Optim.jl for other options
* `bounded_flag=true` : use box constraints, see Optim.jl

# Outputs

* updated `modm` in the input `Param`
* returns the result of optimization as an Optim.jl object
  * `=:migr_finite_difference` same as above but *not* using adjoint state method; time consuming; only for testing, TODO: implement autodiff here
"""
function xfwi!(pa::Param, obj::Union{LS,LS_prior}, attrib_mod;
	   optim_scheme=LBFGS(),
	   optim_options=Optim.Options(show_trace=true,
	   store_trace=true,
	   extended_trace=false,
	   outer_iterations=2,
	   iterations=5,
	   f_tol=1e-5,
	   g_tol=1e-8,
	   x_tol=1e-5, ),
	   bounded_flag=false, solver=nothing,
	   ipopt_options=[["max_iter", 5]])

	global to
	reset_timer!(to)

	println("updating modm and modi...")
	println("> xfwi: number of inversion variables:\t", xfwi_ninv(pa))

	initialize!(pa)

	f(x) = ζfunc(x, pa.mx.last_x,  pa, obj, attrib_mod)
	g!(storage, x) = ζgrad!(storage, x, pa.mx.last_x, pa, obj, attrib_mod)
	if(!bounded_flag)
		"""
		Unbounded LBFGS inversion, only for testing
		"""
		@timeit to "xfwi!" begin
			res = optimize(f, g!, pa.mx.x, optim_scheme, optim_options)
		end
	else
		"""
		Bounded LBFGS inversion
		"""
		if (solver == :ipopt)
			function eval_f(x)
			    return ζfunc(x, pa.mx.last_x,  pa, obj, attrib_mod)
			end

			function eval_grad_f(x, grad_f)
			    ζgrad!(grad_f, x, pa.mx.last_x, pa, obj, attrib_mod)
			end

			function void_g(x, g)

			end

			function void_g_jac(x, mode, rows, cols, values)

			end

			prob = createProblem(size(pa.mx.x)[1], pa.mx.lower_x, pa.mx.upper_x, 0,
				Array{Float64}(undef,0), Array{Float64}(undef,0), 0, 0,
					eval_f, void_g, eval_grad_f, void_g_jac, nothing)

			addOption(prob, "hessian_approximation", "limited-memory")

			if !(ipopt_options === nothing)
			    for i in 1:size(ipopt_options)[1]
				addOption(prob, ipopt_options[i][1], ipopt_options[i][2])
			    end
			end

			prob.x = pa.mx.x

			@timeit to "xfwi!" begin
				res = solveProblem(prob)
			end
		else
			@timeit to "xfwi!" begin
				res = optimize(f, g!, pa.mx.lower_x, pa.mx.upper_x,pa.mx.x, Fminbox(optim_scheme), optim_options)
			end
                end
	end
	println(to)

	# print some stuff after optimization...
	if (solver == :ipopt)
		pa.verbose && println(ApplicationReturnStatus[res])
	else
		pa.verbose && println(res)
	end

	# finalize optimization, using the optimizer...
	if solver == :ipopt
		finalize!(pa, prob.x, attrib_mod)
	else
		finalize!(pa, Optim.minimizer(res), attrib_mod)
	end

	return res
end

"""
Return gradient at the first iteration, i.e., a migration image
"""
function xfwi!(pa::Param, obj::Migr, attrib_mod)

	global to
	reset_timer!(to)

	initialize!(pa)

	g1=pa.mx.gm[1]
	@timeit to "xfwi!" begin
		grad!(g1, pa.mx.x, pa.mx.last_x,  pa, attrib_mod)
	end
	println(to)

	# convert gradient vector to model
	visualize_gx!(pa.gmodm, pa.modm, pa.gmodi, pa.modi, g1, pa)

	return pa.gmodi, g1
end

"""
Return gradient at the first iteration, i.e., a migration image, without using
adjoint state method.
"""
function xfwi!(pa::Param, obj::Migr_fd, attrib_mod)
	global to
	reset_timer!(to)

	initialize!(pa)
	gx=pa.mx.gm[1]

	println("> number of functional evaluations:\t", xfwi_ninv(pa))

	f(x) = ζfunc(x, pa.mx.last_x,  pa, LS(), attrib_mod)

	gx = Calculus.gradient(f,pa.mx.x) # allocating new gradient vector here

	println(to)

	# convert gradient vector to model
	visualize_gx!(pa.gmodm, pa.modm, pa.gmodi, pa.modi, gx, pa)

	return pa.gmodi, gx
end

		#=
		Harvest the Optim result to plot functional and gradient --- to be updated later
		if(extended_trace)
			# convert gradient vector to model
			gmodi = [Models.Seismic_zeros(pa.modi.mgrid) for itr=1:Optim.iterations(res)]
			gmodm = [Models.Seismic_zeros(pa.modm.mgrid) for itr=1:Optim.iterations(res)]
			modi = [Models.Seismic_zeros(pa.modi.mgrid) for itr=1:Optim.iterations(res)]
			modm = [Models.Seismic_zeros(pa.modm.mgrid) for itr=1:Optim.iterations(res)]
			for itr=1:Optim.iterations(res)
				# update modm and modi
				Seismic_x!(modm[itr], modi[itr], Optim.x_trace(res)[itr], pa, -1)

				# update gmodm and gmodi
				Seismic_gx!(gmodm[itr],modm[itr],gmodi[itr],modi[itr],Optim.trace(res)[itr].metadata["g(x)"],pa,-1)
			end

		else
			f = Optim.minimum(res)
		end
		=#

ext: jl


text:
module SocatTimeSeries

using Dates

export SeriesEntry, tostring, makeminuteseries

struct SeriesEntry
  timeindex::Int64
  fco2::Float64

  function SeriesEntry(date::DateTime, fco2::Float64)
    return new(secondofday(date), fco2)
  end

  function SeriesEntry(index::Int64, fco2::Float64)
    return new(index, fco2)
  end
end

function secondofday(date)
  return hour(date) * 3600 + minute(date) * 60 + second(date)
end

function tostring(series::Vector{SeriesEntry})::String
  json = IOBuffer()

  for entry in series
    print(json, "$(entry.timeindex),$(entry.fco2);")
  end

  result = String(take!(json))
  close(json)
  return result
end

function makeminuteseries(secondseries::Vector{SeriesEntry})

  minuteseries = SeriesEntry[]

  sum = 0.0
  count = 0
  currentminute = -1

  for secondentry in secondseries

    entryminute = floor(Int64, secondentry.timeindex / 60)
    if entryminute != currentminute
      if count > 0
        push!(minuteseries, SeriesEntry(currentminute, sum / count))
      end

      sum = 0
      count = 0
      currentminute = entryminute
    end

    sum += secondentry.fco2
    count += 1
  end

  if count > 0
    push!(minuteseries, SeriesEntry(currentminute, sum / count))
  end

  return minuteseries
end

end


text:
using NumericalMethodsforEngineers, IterativeSolvers

a = [10. 1. -5.; -20. 3. 20.; 5. 3. 5.]
b = [1., 2., 6.]

c = a \ b

# p_2_13 Fortran version

x0 = [1., 1., 1.]
tol = 1.0e-5
limit = 1000

# bicgstab defined in src/ch02/bicgstab.jl
iters, converged, x = bicgstab(a, b, x0; tol=tol, limit=limit, trace=true)
println("\nSolution Vector: $x")
println("Converged: $(iters < limit)")
println("Number of iters: $(iters)")

@test round(x, 9) == round(c, 9)


text:
function vol_given_values_and_times(vals::Vector, times::Vector, asset::Symbol)
    duration = maximum(times) - minimum(times)
    returns = simple_differencing(vals[2:end], vals[1:(end-1)])
    return sqrt(sum(returns .^ 2)/duration)
end

"""
    default_num_grids(ts::SortedDataFrame)

This gives a default number of intervals to divide a series of ticks over for the purposes of estimating returns and volatility.
### Inputs
* `ts` - The tick data.
### Returns
* An integer for the number of intervals.
"""
function default_num_grids(ts::SortedDataFrame)
    min_ticks = minimum(map( a -> length(ts.groupingrows[a]) , collect(keys(ts.groupingrows)) ))
    return Int(max(floor(min_ticks / 100), 3))
end

"""
    two_scales_volatility(vals::Vector, times::Vector, asset::Symbol, num_grids::Real)

Calculates volatility with the two scales method of Zhang, Mykland, Ait-Sahalia 2005. The amount of time for the grid spacing is by default this is a tenth of the total duration
by default. If this doesn't make sense for your use of it then choose a spacing at which you expect the effect of microstructure noise will be small.
### Inputs
* `vals` - The prices at each instant in time.
* `times` - The times corresponding to each element in `vals`.
* `asset` - The name of the asset.
* `num_grids` - Number of grids used in order in two scales estimation.
### Returns
* A scalar for the estimated volatility of the asset.
* A scalar for the estimated microstructure noise variance.

    two_scales_volatility(ts::SortedDataFrame, assets::Vector{Symbol} = get_assets(ts);
                          num_grids::Real = default_num_grids(ts))

Calculates volatility with the two scales method of Zhang, Mykland, Ait-Sahalia 2005. The amount of time for the grid spacing is by default this is a tenth of the total duration
by default. If this doesn't make sense for your use of it then choose a spacing at which you expect the effect of microstructure noise will be small.
### Inputs
* `ts` - The tick data.
* `assets` - The assets you want to estimate volatilities for.
* `num_grids` - Number of grids used in order in two scales estimation.
### Returns
* A  `Dict` with estimated volatilities for each asset.
* A  `Dict` with estimated microstructure noise variances for each asset.

### References
Zhang L, Mykland PA, Aït-Sahalia Y (2005). "A Tale of Two Time Scales: Determining Integrated Volatility with Noisy High-Frequency Data." Journal of the American Statistical Association, 100(472), 1394–1411. ISSN 01621459. doi:10.1198/016214505000000169.
"""
function two_scales_volatility(vals::Vector, times::Vector, asset::Symbol, num_grids::Real)
    dura  = maximum(times) - minimum(times)
    if (dura < eps()) | (length(vals) < 10) return NaN, NaN end
    num_grids = Int(  min( floor(length(times)/4),   max( 2,  floor(num_grids)  ) ))
    avg_vol   = mean(map(i ->  vol_given_values_and_times(vals[i:num_grids:end], times[i:num_grids:end], asset), 1:num_grids ))
    all_vol   = vol_given_values_and_times(vals, times, asset)

    pure_vol  = ((1- 1/num_grids)^(-1)) * ( avg_vol - (1/num_grids) *all_vol   )
    noise     = (1/(2*length(vals))) * ((all_vol^2)*dura - (pure_vol^2)*dura )
    return pure_vol, noise
end

function two_scales_volatility(ts::SortedDataFrame, assets::Vector{Symbol} = get_assets(ts);
                               num_grids::Real = default_num_grids(ts))
    vols = Dict{Symbol,eltype(ts.df[:,ts.value])}()
    micro_noise_var = Dict{Symbol,eltype(ts.df[:,ts.value])}()
    for a in assets
        vals = ts.df[ts.groupingrows[a],ts.value]
        times = ts.df[ts.groupingrows[a],ts.time]
        if length(vals) < 10
            @warn "There was not enough data to use the two_scales_volatility method. Returning NaN."
            vols[a] = NaN
            micro_noise_var[a] = NaN
        else
            pure_vol, noise = two_scales_volatility(vals, times, a, num_grids)
            vols[a] = pure_vol
            micro_noise_var[a] = max(0.0,noise)
        end
    end
    return vols, micro_noise_var
end


text:
my_f(x,y) = 2x + y

my_ff(x,y) = 2x + y + 1

max_stars_repo_name: aivodji/Yg

max_issues_repo_name: aivodji/Yg

max_forks_repo_name: aivodji/Yg


text:
using CommutativeRings: len, GFImpl

@testset "iterate" begin
    @test length(collect(ZZ/7)) == 7
    @test length(collect(GF(3,5))) == 3^5
    @test length(collect(Monic((ZZ/2)[:x], 5))) == 2^5
    @test eltype(Monic((ZZ/2)[:x], 5)) == (ZZ/2)[:x]
end

@testset "ofindex $T" for (T, r) in ((Int, 0), (UInt, 0), (ZZ{Int8}, 0), (ZZ/77, 77), (GF(2,5),2^5), (QQ{Int}, 0))
    @test iszero(ofindex(0, T))
    @test isone(ofindex(1, T))
    @test ofindex(10, T) isa T
    @test len(T) == r
end

@testset "ofindex $T degree $n" for (T, n) in ((GF(2)[:z], 3),)
    @test ofindex(0, T, n) == monom(T, n)
    @test ofindex(1, T, n) == monom(T, n) + 1
    @test ofindex(10, T, n) isa T
end

@testset "ofindex Quotient" begin
    P = (ZZ/5)[:x]
    x = monom(P)
    Q = P / ( x^3 + x + 1)
    n = length(Q)
    @test length(Q) == 125
    @test sort(ofindex.(0:n-1, Ref(Q))) == sort(collect(Q))
end

@testset "factors" begin
    @test length(factors(12)) == 6
    @test factors(12) |> collect |> sort == [1, 2, 3, 4, 6, 12]
end


text:
using Dates
using DataStructures

"""
    UncertaintyNDistribution

describes the uncertainty distribution of an injection (load or renewable)
We consider `ts` as the target timepoint at which the realisation will have place,
and `ech` is the observation point.
The farthest we are from ts, the greater is the uncertainty: we assume an increase  by a ratio of `time_factor` for each hour.
i.e.: factor_l = max{ 0, ( ts - ech )/1h } * time_factor
adjusted_sigma = factor_l * sigma
random_injection = Normal( mu,  adjusted_sigma)

# Arguments
    - `min_value::Float64` : minimum value that can be assigned to the injection
    - `max_value::Float64` : maxiimum value that can be assigned to the injection
    - `mu::Float64` : average injection value
    - `sigma::Float64` : base sigma value to be used as a standard deviation value for the injection uncertainty distribution
    - `time_factor::Float64` : per hour increase ratio of the uncertainty's standard deviation
"""
struct UncertaintyNDistribution
    id::String
    min_value::Float64
    max_value::Float64
    mu::Float64
    sigma::Float64
    time_factor::Float64
end

UncertaintiesDistribution = SortedDict{String, UncertaintyNDistribution}

function add_uncertainty_distribution!(uncertainties_distribution::UncertaintiesDistribution,
                        id::String,
                        min_value::Float64, max_value::Float64,
                        mu::Float64, sigma::Float64,
                        time_factor::Float64=1.)
    uncertainties_distribution[id] = UncertaintyNDistribution(id, min_value, max_value,
                                                            mu, sigma,
                                                            time_factor)
end

# uncertainties[ech][nodal_injection_name][ts][scenario_name]
InjectionUncertainties = SortedDict{Dates.DateTime,
                            SortedDict{String, Float64}
                            }
UncertaintiesAtEch = SortedDict{String, InjectionUncertainties}
Uncertainties = SortedDict{Dates.DateTime, UncertaintiesAtEch}

function add_uncertainty!(uncertainties_at_ech::UncertaintiesAtEch, nodal_injection_name::String, ts::Dates.DateTime, scenario_name::String, val::Float64)
    get!(uncertainties_at_ech, nodal_injection_name, SortedDict{Dates.DateTime, SortedDict{String, Float64}}())
    get!(uncertainties_at_ech[nodal_injection_name], ts, SortedDict{String, Float64}())

    uncertainties_at_ech[nodal_injection_name][ts][scenario_name] = val
    return uncertainties_at_ech
end
function add_uncertainty!(uncertainties::Uncertainties, ech::Dates.DateTime, nodal_injection_name::String, ts::Dates.DateTime, scenario_name::String, val::Float64)
    uncertainties_at_ech = get!(uncertainties, ech, SortedDict{String, SortedDict{Dates.DateTime, SortedDict{String, Float64}}}())
    add_uncertainty!(uncertainties_at_ech, nodal_injection_name, ts, scenario_name, val)
    return uncertainties
end

function get_scenarios(uncertainties_at_ech::UncertaintiesAtEch)::Vector{String}
    for (_, uncertainties_by_ts) in uncertainties_at_ech
        for (_, by_scenario) in uncertainties_by_ts
            return collect(keys(by_scenario))
        end
    end
    return Vector{String}()
end

function get_scenarios(uncertainties::Uncertainties)::Vector{String}
    for (_, uncertainties_at_ech) in uncertainties
        return get_scenarios(uncertainties_at_ech)
    end
    return Vector{String}()
end

function get_target_timepoints(uncertainties::Uncertainties)::Vector{Dates.DateTime}
    for (_, uncertainties_at_ech) in uncertainties
        for (_, by_ts) in uncertainties_at_ech
            return collect(keys(by_ts))
        end
    end
    return Vector{Dates.DateTime}()
end

function get_horizon_timepoints(uncertainties::Uncertainties)::Vector{Dates.DateTime}
    return collect(keys(uncertainties))
end

function get_uncertainties(uncertainties::Uncertainties, ech::Dates.DateTime)::UncertaintiesAtEch
    return uncertainties[ech]
end

function get_uncertainties(uncertainties_at_ech::UncertaintiesAtEch, injection_name::String)::InjectionUncertainties
    return uncertainties_at_ech[injection_name]
end
function get_uncertainties(uncertainties::Uncertainties, ech::Dates.DateTime, injection_name::String)::InjectionUncertainties
    uncertainties_at_ech = get_uncertainties(uncertainties, ech)
    return get_uncertainties(uncertainties_at_ech, injection_name)
end

function get_uncertainties(injection_uncertainties::InjectionUncertainties, ts::Dates.DateTime)
    return injection_uncertainties[ts]
end
function get_uncertainties(uncertainties_at_ech::UncertaintiesAtEch, injection_name::String, ts::Dates.DateTime)
    injection_uncertainties = get_uncertainties(uncertainties_at_ech, injection_name)
    return get_uncertainties(injection_uncertainties, ts)
end
function get_uncertainties(uncertainties::Uncertainties, ech::Dates.DateTime, injection_name::String, ts::Dates.DateTime)
    uncertainties_at_ech = get_uncertainties(uncertainties, ech)
    return get_uncertainties(uncertainties_at_ech, injection_name, ts)
end

function get_uncertainties(injection_uncertainties::InjectionUncertainties, ts::Dates.DateTime, scenario::String)
    return get_uncertainties(injection_uncertainties, ts)[scenario]
end
function get_uncertainties(uncertainties_at_ech::UncertaintiesAtEch, injection_name::String, ts::Dates.DateTime, scenario::String)
    injection_uncertainties = get_uncertainties(uncertainties_at_ech, injection_name)
    return get_uncertainties(injection_uncertainties, ts, scenario)
end
function get_uncertainties(uncertainties::Uncertainties, ech::Dates.DateTime, injection_name::String, ts::Dates.DateTime, scenario::String)
    uncertainties_at_ech = get_uncertainties(uncertainties, ech)
    return get_uncertainties(uncertainties_at_ech, injection_name, ts, scenario)
end

function sum_uncertainties(uncertainties_at_ech::UncertaintiesAtEch,
    ids::Vector{String},
    ts::Dates.DateTime, s::String)
    sum = 0.
    for id in ids
        sum += get_uncertainties(uncertainties_at_ech, id, ts)[s]
    end
    return sum
end

function compute_prod(uncertainties_at_ech::UncertaintiesAtEch,
    limitables_gens, ts::Dates.DateTime, s::String)
    generators_ids = map(Networks.get_id, limitables_gens)
    return sum_uncertainties(uncertainties_at_ech, generators_ids, ts, s)
end
function compute_prod(uncertainties_at_ech::UncertaintiesAtEch,
    network::Networks.Network, ts::Dates.DateTime, s::String)
    limitables_list = Networks.get_generators_of_type(network, Networks.LIMITABLE)
    return compute_prod(uncertainties_at_ech, limitables_list, ts, s)
end
function compute_prod(uncertainties::Uncertainties,
    network::Networks.Network, ech::Dates.DateTime, ts::Dates.DateTime, s::String)
    uncertainties_at_ech = get_uncertainties(uncertainties, ech)
    return compute_prod(uncertainties_at_ech, network, ts, s)
end

function compute_load(uncertainties_at_ech::UncertaintiesAtEch,
    buses_list, ts::Dates.DateTime, s::String)
    buses_ids = map(Networks.get_id, buses_list)
    return sum_uncertainties(uncertainties_at_ech, buses_ids, ts, s)
end
function compute_load(uncertainties_at_ech::UncertaintiesAtEch,
    network::Networks.Network, ts::Dates.DateTime, s::String)
    return compute_load(uncertainties_at_ech, Networks.get_buses(network), ts, s)
end
function compute_load(uncertainties::Uncertainties,
    network::Networks.Network, ech::Dates.DateTime, ts::Dates.DateTime, s::String)
    uncertainties_at_ech = get_uncertainties(uncertainties, ech)
    return compute_load(uncertainties_at_ech, network, ts, s)
end

max_stars_repo_name: klorel/pscopf

max_issues_repo_name: klorel/pscopf

max_forks_repo_name: klorel/pscopf


text:
module Poisson2D
export calculate_potential!

function calculate_potential!()
end

end


text:
module StatsbombPitch

using Plots

"""
plot_pitch()

Plots a football pitch according to Statsbomb dimensions

Arguments:
line_color = color of the line
line_width = width of the lines

"""
function plot_pitch(;line_color, line_width)

    x1 = [0, 120, 120, 0, 0]; y1 = [0, 0, 80, 80, 0] ## pitch outlines
    x2 = [0, 18, 18, 0, 0]; y2 = [18, 18, 62, 62, 18] ## left d-box
    x3 = [102, 120, 120, 102, 102]; y3 = [18, 18, 62, 62, 18] ## right d-box
    x4 = [0, 6, 6, 0, 0]; y4 = [30, 30, 50, 50, 30] ## left 6-yard box
    x5 = [114, 120, 120, 114, 114]; y5 = [30, 30, 50, 50, 30] ## right 6-yard box

    ##goals
    x6 = [0,0]; y6 = [36, 44]
    x7 = [120, 120]; y7=y6

    ##halfway line and penalty spots
    x8 = [60,60]; y8 = [0,80]

    p = plot(x1, y1, color=line_color, lw=line_width)
    plot!(p, x2, y2, color=line_color, lw=line_width)
    plot!(p, x3, y3, color=line_color, lw=line_width)
    plot!(p, x4, y4, color=line_color, lw=line_width)
    plot!(p, x5, y5, color=line_color, lw=line_width)
    plot!(p, x6, y6, color=line_color, lw=3*line_width)
    plot!(p, x7, y7, color=line_color, lw=3*line_width)
    plot!(p, x8, y8, color=line_color, lw=line_width)
    scatter!(p, [12, 60, 108], [40, 40, 40], color=line_color, markersize=2)

    xs, ys = Plots.unzip(Plots.partialcircle(0, 2π, 100, 4))

    plot!(p, xs[1:25], ys[1:25], color=line_color, lw=line_width)
    plot!(p, xs[25:50].+120, ys[25:50], color=line_color, lw=line_width)
    plot!(p, xs[50:75].+120, ys[50:75].+80, color=line_color, lw=line_width)
    plot!(p, xs[75:100], ys[75:100].+80, color=line_color, lw=line_width)

    xs, ys = Plots.unzip(Plots.partialcircle(0, 2π, 100, 8))
    plot!(p, xs.+60, ys.+40, linecolor = line_color, fillalpha = 0.01, aspect_ratio = 1, lw=line_width)
    plot!(p, legend=false, grid=false, border=:none, xlims=(-5, 125), ylims=(-3, 83))

    ys.+=40; xs.+=12
    plot!(p, xs[xs.>18], ys[xs.>18], color=line_color, lw=line_width)

    xs.+=96
    plot!(p, xs[xs.<102], ys[xs.<102], color=line_color, lw=line_width)
    plot!(p, yflip=true)

    p

end

end


text:
@userplot AbundancePlot
@recipe function f(plt::AbundancePlot; topabund=10, srt=collect(1:nsites(plt.args[1])))
    abun = plt.args[1]
    typeof(abun) <: AbstractComMatrix || error("AbundancePlot not defined for $(typeof(abun))")

    topabund = min(topabund, nfeatures(abun))
    2 <= topabund < 12 || error("n must be between 2 and 12")

    top = filterabund(abun, topabund)

    rows = featurenames(top)

    yflip := true
    bar_position := :stack
    label := rows
    GroupedBar((1:nsamples(top), Matrix(occurrences(top)[:,srt]')))
end

"""
    AnnotationBar(::Array{<:Color,1})

Array of colors for plotting a grid. Use `plot(::AnnotationBar; kwargs...)`
"""
struct AnnotationBar
    colors::Array{<:Color,1}
end

"""
    annotationbar(colors::AbstractArray{<:Color,1})
    annotationbar(labels::AbstractArray{<:AbstractString,1})
    annotationbar(labels::AbstractArray{<:AbstractString,1}, colormap::Dict{<:AbstractString,<:Color})
    annotationbar(labels::AbstractArray{<:AbstractString,1}, colors::AbstractArray{<:Color,1})

Create an `AnnotationBar` from a vector of colors or a vector of `labels` (`String`s).
- A vector of `Color`s will be used directly to create an `AnnotationBar`
- For an array of `String`s, each unique value will be assigned a random color from
  `Colors.color_names`.
- A vector of labels may be passed with a `Dict` with keys for each `label` mapping
  to a `Color` value.
- A vector of `labels` and `colors` may be used.
    - If the lengths of these vectors are the same,
      the `colors` will be used to create an AnnotationBar
    - If the length of `labels` is longer than the length of `colors`,
      each unique label will be assigned to one of the `colors`
"""
annotationbar(colors::AbstractArray{<:Color,1}) = AnnotationBar(colors)

function annotationbar(labels::AbstractArray{<:AbstractString,1}, colormap::Dict{<:AbstractString,<:Color})
    colorkeys = keys(colormap)
    all(label-> in(label, colorkeys), labels) || throw(ArgumentError("Colormap missing keys $(setdiff(labels, colorkeys))"))
    AnnotationBar([colormap[label] for label in labels])
end

function annotationbar(labels::AbstractArray{<:AbstractString,1}, colors::AbstractArray{<:Color,1})
    if length(labels) == length(colors)
        return annotationbar(colors)
    end
    ulabels = unique(labels)
    length(colors) < length(ulabels) && throw(ArgumentError("Must have at least 1 color for each unique label"))
    colormap = Dict(l => colors[i] for (i, l) in enumerate(ulabels))
    annotationbar(labels, colormap)
end

# if no colors are passed, give a unique random color to each label
function annotationbar(labels::AbstractArray{<:AbstractString,1})
    ulabels = unique(labels)
    colors = rand(keys(Colors.color_names), length(ulabels))
    colormap = Dict(l => parse(Colorant, colors[i]) for (i, l) in enumerate(ulabels))
    annotationbar(labels, colormap)
end

@recipe function f(bar::AnnotationBar)
    xs = Int[]
    for i in 1:length(bar.colors)
        append!(xs, [0,0,1,1,0] .+ (i-1))
    end
    xs = reshape(xs, 5, length(bar.colors))
    ys = hcat([[0,1,1,0,0] for _ in bar.colors]...)

    fc = reshape(bar.colors, 1, length(bar.colors))

    seriestype := :path
    fill := (0,1)
    fillcolor := fc
    legend := false
    color --> :black
    ticks := false
    framestyle := false
    xaxis --> false
    yaxis --> false
    xs, ys
end

# From Michael K. Borregaard (posted to slack 2018-05-18)
# Usage:
#
# x = randn(100)
# y = randn(100) + 2
# y[y.<0] = 0
#
# zeroyplot(x,y, yscale = :log10, color = :red, markersize = 8)
# @userplot ZeroYPlot
# @recipe function f(h::ZeroYPlot)
#            length(h.args) != 2 && error("zeroyplot only defined for x and y input arguments")
#            x, y = h.args
#            val0 = y .== 0
#
#            layout := @layout([a
#                               b{0.05h}])
#
#            markeralpha --> 0.4
#            seriestype --> :scatter
#
#            @series begin
#                primary := false
#                subplot := 2
#                yscale := :identity
#                ylim := (-1,1)
#                yticks := ([0],[0])
#                grid := false
#                x[val0], y[val0]
#            end
#
#            subplot := 1
#            x[.!val0], y[.!val0]
# end


text:
using Documenter, MixFit

makedocs(sitename="MixFit Documentation", modules = [MixFit])
deploydocs(repo = "github.com/the-sushi/MixFit.jl.git")


text:
export ldlt_symm

"""
Higams' ldlt_symm translated from matlab. Performs  a so called
 BKK  Bounded Bunch Kaufman factorization of A0, that means
 ||L|| is bounded and bounded away from zero.
"""
function ldlt_symm(A0::Array{T,2}, piv::Char = 'r') where {T}

    #LDLT_SYMM  Block LDL^T factorization for a symmetric indefinite matrix.
    #     Given a Hermitian matrix A,
    #     [L, D, P, RHO, NCOMP] = LDLT_SYMM(A, PIV) computes a permutation P,
    #     a unit lower triangular L, and a real block diagonal D
    #     with 1x1 and 2x2 diagonal blocks, such that  P*A*P' = L*D*L'.
    #     PIV controls the pivoting strategy:
    #       PIV = 'p': partial pivoting (Bunch and Kaufman),
    #       PIV = 'r': rook pivoting (Ashcraft, Grimes and Lewis).
    #     The default is partial pivoting.
    #     RHO is the growth factor.
    #     For PIV = 'r' only, NCOMP is the total number of comparisons.

    #     References:
    #     J. R. Bunch and L. Kaufman, Some stable methods for calculating
    #        inertia and solving symmetric linear systems, Math. Comp.,
    #        31(137):163-179, 1977.
    #     C. Ashcraft, R. G. Grimes and J. G. Lewis, Accurate symmetric
    #        indefinite linear equation solvers. SIAM J. Matrix Anal. Appl.,
    #        20(2):513-561, 1998.
    #     N. J. Higham, Accuracy and Stability of Numerical Algorithms,
    #        Second edition, Society for Industrial and Applied Mathematics,
    #        Philadelphia, PA, 2002; chap. 11.

    #    This routine does not exploit symmetry and is not designed to be
    #     efficient.
    #
    #    Adapted from N. Higham Matlab code.
    #    JPD april 23 2015

    # Since array contents are mutable and modified,
    # copy the input matrix to keep it unchanged outside the
    # scope of the function

    A = copy(A0)

    # minimal checks for conforming inputs
    isequal(triu(A)', tril(A)) || error("Must supply Hermitian matrix.")
    piv in ['p', 'r'] || error("Pivoting must be \"'p'\" or \" 'r' \".")

    n, = size(A)

    k = 1
    # D = eye(n,n);
    # L = eye(n,n);
    global D = Matrix{T}(1.0I, n, n)
    global L = Matrix{T}(1.0I, n, n)
    if n == 1
        D = A
    end

    global pp = collect(1:n)

    maxA = norm(A, Inf)
    global ρ = maxA

    global ncomp = 0
    s = 1

    α = T((1 + sqrt(17)) / 8)
    while k < n
        (λ, vr) = findmax(abs.(A[k+1:n, k]))
        r = vr[1] + k
        if λ > 0
            swap = false
            if abs.(A[k, k]) >= α * λ
                s = 1
            else
                if piv == 'p'
                    temp = A[k:n, r]
                    temp[r-k+1] = 0
                    σ = norm(temp, Inf)
                    if α * λ^2 <= abs.(A[k, k]) * σ
                        s = 1
                    elseif abs.(A[r, r]) >= α * σ
                        swap = true
                        m1 = k
                        m2 = r
                        s = 1
                    else
                        swap = true
                        m1 = k + 1
                        m2 = r
                        s = 2
                    end
                    if swap
                        A[[m1, m2], :] = A[[m2, m1], :]
                        L[[m1, m2], :] = L[[m2, m1], :]
                        A[:, [m1, m2]] = A[:, [m2, m1]]
                        L[:, [m1, m2]] = L[:, [m2, m1]]

                        pp[[m1, m2]] = pp[[m2, m1]]
                    end
                elseif piv == 'r'
                    j = k
                    pivot = false
                    λ_j = λ
                    while ~pivot
                        (temp1, vr) = findmax(abs.(A[k:n, j]))
                        global ncomp = ncomp + n - k
                        r = vr[1] + k - 1
                        temp = A[k:n, r]
                        temp[r-k+1] = 0.0
                        λᵣ = max(maximum(temp), -minimum(temp))
                        global ncomp = ncomp + n - k
                        if α * λᵣ <= abs.(A[r, r])
                            pivot = true
                            s = 1
                            A[k, :], A[r, :] = A[r, :], A[k, :]
                            L[k, :], L[r, :] = L[r, :], L[k, :]
                            A[:, k], A[:, r] = A[:, r], A[:, k]
                            L[:, k], L[:, r] = L[:, r], L[:, k]
                            pp[k], pp[r] = pp[r], pp[k]
                        elseif λ_j == λᵣ
                            pivot = true
                            s = 2
                            A[k, :], A[j, :] = A[j, :], A[k, :]
                            L[k, :], L[j, :] = L[j, :], L[k, :]
                            A[:, k], A[:, j] = A[:, j], A[:, k]
                            L[:, k], L[:, j] = L[:, j], L[:, k]
                            pp[k], pp[j] = pp[j], pp[k]
                            k1 = k + 1
                            A[k1, :], A[r, :] = A[r, :], A[k1, :]
                            L[k1, :], L[r, :] = L[r, :], L[k1, :]
                            A[:, k1], A[:, r] = A[:, r], A[:, k1]
                            L[:, k1], L[:, r] = L[:, r], L[:, k1]
                            pp[k1], pp[r] = pp[r], pp[k1]
                        else
                            j = r
                            s = 1
                            λ_j = λᵣ
                        end
                    end
                end
            end

            if s == 1

                D[k, k] = A[k, k]
                A[k+1:n, k] = A[k+1:n, k] / A[k, k]
                L[k+1:n, k] = A[k+1:n, k]
                i = k+1:n
                A[i, i] = A[i, i] - A[i, k:k] * A[k:k, i]
                A[i, i] = 0.5 * (A[i, i] + A[i, i]')

            elseif s == 2

                E = A[k:k+1, k:k+1]
                D[k:k+1, k:k+1] = E
                i = k+2:n
                C = A[i, k:k+1]
                temp = C / E
                L[i, k:k+1] = temp
                A[i, k+2:n] = A[i, k+2:n] - temp * C'
                A[i, i] = 0.5 * (A[i, i] + A[i, i]')
            end

            # Ensure diagonal real (see LINPACK User's Guide, p. 5.17).
            for i = k+s:n
                A[i, i] = real(A[i, i])
            end

            if k + s <= n
                val, = findmax(abs.(A[k+s:n, k+s:n]))
                global ρ = max(ρ, val)
            end

        else  # Nothing to do.
            s = 1
            D[k, k] = A[k, k]
        end

        k = k + s

        if k == n
            D[n, n] = A[n, n]
            break
        end

    end
    #P=eye(n)
    #P = P[pp,:]
    return L, D, pp, ρ, ncomp
end


text:
baremodule Monoids
    using ..Types
end
MonoidUnion = Union{AbstractMonoid, libgb.GrB_Monoid}
function monoidnames(name)
    simple = splitconstant(name)[2]
    containername = Symbol(simple, "_MONOID_T")
    exportedname = Symbol(simple * "_MONOID")
    return containername, exportedname
end
function createmonoids()
    builtins = ["GrB_MIN", "GrB_MAX", "GrB_PLUS", "GrB_TIMES", "GxB_ANY", "GrB_LOR", "GrB_LAND", "GrB_LXOR", "GrB_LXNOR", "GxB_EQ", "GxB_BOR", "GxB_BAND", "GxB_BXOR", "GxB_BXNOR"]
    for name ∈ builtins
        containername, exportedname = monoidnames(name)
        structquote = quote
            struct $containername <: AbstractMonoid
                pointers::Dict{DataType, libgb.GrB_Monoid}
                name::String
                $containername() = new(Dict{DataType, libgb.GrB_Monoid}(), $name)
            end
        end
        @eval(Types, $structquote)
        constquote = quote
            const $exportedname = Types.$containername()
            export $exportedname
        end
        @eval(Monoids, $constquote)
    end
end

function load(monoid::AbstractMonoid)
    booleans = ["GxB_ANY", "GrB_LOR", "GrB_LAND", "GrB_LXOR", "GrB_LXNOR", "GxB_EQ"]
    integers = ["GrB_MIN", "GrB_MAX", "GrB_PLUS", "GrB_TIMES", "GxB_ANY"]
    unsignedintegers = ["GrB_MIN", "GrB_MAX", "GrB_PLUS", "GrB_TIMES", "GxB_ANY", "GxB_BOR", "GxB_BAND", "GxB_BXOR", "GxB_BXNOR"]
    floats = ["GrB_MIN", "GrB_MAX", "GrB_PLUS", "GrB_TIMES", "GxB_ANY"]
    name = monoid.name

    if name ∈ booleans
        constname = name * ((isGxB(name) ? "_BOOL_MONOID" : "_MONOID_BOOL"))
        monoid.pointers[Bool] = load_global(constname)
    end

    if name ∈ integers
        monoid.pointers[Int8] = load_global(name * (isGxB(name) ? "_INT8_MONOID" : "_MONOID_INT8"))
        monoid.pointers[Int16] = load_global(name * (isGxB(name) ? "_INT16_MONOID" : "_MONOID_INT16"))
        monoid.pointers[Int32] = load_global(name * (isGxB(name) ? "_INT32_MONOID" : "_MONOID_INT32"))
        monoid.pointers[Int64] = load_global(name * (isGxB(name) ? "_INT64_MONOID" : "_MONOID_INT64"))
    end

    if name ∈ unsignedintegers
        monoid.pointers[UInt8] = load_global(name * (isGxB(name) ? "_UINT8_MONOID" : "_MONOID_UINT8"))
        monoid.pointers[UInt16] = load_global(name * (isGxB(name) ? "_UINT16_MONOID" : "_MONOID_UINT16"))
        monoid.pointers[UInt32] = load_global(name * (isGxB(name) ? "_UINT32_MONOID" : "_MONOID_UINT32"))
        monoid.pointers[UInt64] = load_global(name * (isGxB(name) ? "_UINT64_MONOID" : "_MONOID_UINT64"))
    end

    if name ∈ floats
        monoid.pointers[Float32] = load_global(name * (isGxB(name) ? "_FP32_MONOID" : "_MONOID_FP32"))
        monoid.pointers[Float64] = load_global(name * (isGxB(name) ? "_FP64_MONOID" : "_MONOID_FP64"))
    end
end

operator(monoid::MonoidUnion) = libgb.GxB_Monoid_operator(monoid)
xtype(monoid::MonoidUnion) = xtype(operator(monoid))
ytype(monoid::MonoidUnion) = ytype(operator(monoid))
ztype(monoid::MonoidUnion) = ztype(operator(monoid))


text:
#------------------------------
# Convenience function for in benchmarking
#------------------------------

# define convenience function to generate data for simpler benchmark loop
function generate_all(n, p, q, q_info; seed::Int=12)
    # set seed for reproducibility
    Random.seed!(seed)
    # choose simulation setting for baseline variables
    baseline_simulation = "groupsonly" # alternatives: "trueparams" or "groupsonly"
    # generate time dependent variables
    xs, tvals, group1, group2 = generate_xs(n, p, true_u0, sol_group1, sol_group2, maxntps=1, σ_var=0.1, σ_ind=0.1); # default vals: t_start=1.5, t_end=10, maxntps=10, dt=0.1, σ_var=0.1, σ_ind=0.5
    # generate baseline variables
    if baseline_simulation == "trueparams"
        x_baseline = generate_baseline(n, q, q_info, group1, true_odeparams_group1, true_odeparams_group2); # defaul vals: σ_info=0.1, σ_noise=0.1
    elseif baseline_simulation == "groupsonly"
        x_baseline = generate_baseline(n, q, q_info, group1,  σ_info=0.5, σ_noise=0.5); # default vals: σ_info=1, σ_noise=1
    else
        error("Please select either 'trueparams' or 'groupsonly' as mode for simulation of baseline variables")
    end
    return xs, x_baseline, tvals, group1, group2
end

function run_benchmark(trainingdata, m)
    L = loss_wrapper(m)
    ps = getparams(m)
    opt = ADAM(0.0005)
    for epoch in 1:35
        Flux.train!(L, ps, trainingdata, opt)
    end
end

using Printf
function prettymemory(b)
    if b < 1024
        return string(b, " bytes")
    elseif b < 1024^2
        value, units = b / 1024, "KiB"
    elseif b < 1024^3
        value, units = b / 1024^2, "MiB"
    else
        value, units = b / 1024^3, "GiB"
    end
    return string(@sprintf("%.2f", value), " ", units)
end


text:
include("rmstest.jl")
"""
Runs all files whose names start with "Test" in the current directory
or any directories below it
"""
function runtests(path=".")
    for (root,dirs,files) in walkdir(path)
        for dir in dirs
            runtests(joinpath(root,dir))
        end
        for file in files
            if startswith(file,"Test")
                println(joinpath(root,file))
                include(joinpath(root,file))
            end
        end
    end
end

runtests()

@test true


text:
using Genie, Genie.Router, Genie.Renderer.Json

Genie.config.run_as_server = true
Genie.config.cors_allowed_origins = ["*"]

route("/random", method=POST) do
  dim = parse(Int, get(params(), :dim, "2"))
  num = parse(Int, get(params(), :num, "3"))

  (:random => rand(dim,num)) |> json
end

up(; open_browser = false)


text:
function remove!(a, item)
    deleteat!(a, findall(x->x==item, a))
end

function token_login(sock::TCPSocket)
    csrf_token = repr(uuid4().value)
    println(sock, "CSRF#> ", csrf_token)
    print(sock, "CSRF#> ")
    str = readline(sock)
    if str == csrf_token
        return true
    end
    return false
end

function local_repl(tcp_sock::TCPSocket)
    server = listen(6969)
    chn = Vector{String}()
    m = ""
    while true
        sock::TCPSocket = accept(server)
        @async begin
            try
                write(sock,"Connected to bot repl\r\n")
                if !(token_login(sock))
                    close(sock)
                end
                while true
                    write(sock, "#> ")
                    str = readline(sock)
                    try
                        repl_commands(tcp_sock, str, chn)
                    catch ex
                        println(sock, "Exception on: ", str)
                        println(sock, "with $ex")
                        close(sock)
                    end
                    if str == "ls"
                        if chn != ""
                            println(sock, "Channels joined:")
                            for (i, ch) in enumerate(chn)
                                println(sock, i, ": ", ch)
                            end
                        else
                            println(sock, "No channel joined")
                        end
                    end
                end
            catch y
                println(stderr, "Caught exception: $y")
                close(sock)
            end
        end
    end
end

function repl_commands(tcp_sock::TCPSocket, str::String, chn::Vector{String})
    if ( m = match(r"(?<cmd>.*) :(?<index>\w+) (?<msg>.*)", str)) !== nothing
        if occursin("say", m[:cmd])
            println(str)
            msg = m[:msg]
            ch_index = parse(Int32, m[:index])
            Utils.irc_send(tcp_sock, chn[ch_index], String(msg))
        end
    elseif  ( m = match(r"(.*) (.*)", str) ) !== nothing
        if m.captures[1] == "join"
            temp_chn = m.captures[2]
            push!(chn, temp_chn)
            err = Utils.irc_join(tcp_sock, String(temp_chn))
            println(stdout, "irc_join: ", err)
        elseif m.captures[1] == "part"
            if m.captures[2] == "all"
                for ch in chn
                    Utils.irc_send(tcp_sock, "PART #$ch")
                end
                chn = String[]
            else
                old_chn = m.captures[2]
                remove!(chn, old_chn)
                Utils.irc_send(tcp_sock, "PART #$old_chn")
            end
        elseif occursin("say", m.captures[1])
            msg = split(str, "say", limit=2)[2]
            Utils.irc_send(tcp_sock, String(chn[1]), String(msg))
        end
    end

    return chn
end

max_stars_repo_name: Iaggelis/bald-bot

max_issues_repo_name: Iaggelis/bald-bot

max_forks_repo_name: Iaggelis/bald-bot


text:
module mbas001
using FinEtools
using Test
function test()
   fens,fes = Q4block(1.3, 3.1,3,2); # Mesh
   @test manifdim(fes) == 2
   bfes = meshboundary(fes)
   @test manifdim(bfes) == 1
   true
end
end
using .mbas001
mbas001.test()

module mbas002
using FinEtools
using Test
function test()
   fens,fes = H8block(1.3, 3.1, 2.7, 3, 2, 4); # Mesh
   @test manifdim(fes) == 3
   bfes = meshboundary(fes)
   @test manifdim(bfes) == 2
   true
end
end
using .mbas002
mbas002.test()

module mbas003
using FinEtools
using Test
function test()
   fens,fes = T4block(1.3, 3.1, 2.7, 3, 2, 4); # Mesh
   @test manifdim(fes) == 3
   bfes = meshboundary(fes)
   @test manifdim(bfes) == 2
   true
end
end
using .mbas003
mbas003.test()

module mbas004
using FinEtools
using Test
function test()
   fens,fes = Q4block(1.3, 3.1,3,2); # Mesh
   @test nodesperelem(fes) == 4
   bfes = meshboundary(fes)
   @test nodesperelem(bfes) == 2
   true
end
end
using .mbas004
mbas004.test()

module mbas005
using FinEtools
using Test
function test()
   fens,fes = H8block(1.3, 3.1, 2.7, 3, 2, 4); # Mesh
   @test nodesperelem(fes) == 8
      bfes = meshboundary(fes)
      @test nodesperelem(bfes) == 4
   true
end
end
using .mbas005
mbas005.test()

module mbas006
using FinEtools
using Test
function test()
   fens,fes = T4block(1.3, 3.1, 2.7, 3, 2, 4); # Mesh
   @test nodesperelem(fes) == 4
      bfes = meshboundary(fes)
      @test nodesperelem(bfes) == 3
   true
end
end
using .mbas006
mbas006.test()

module mbas007
using FinEtools
using Test
function test()
   fens,fes = T10block(1.3, 3.1, 2.7, 3, 2, 4); # Mesh
   @test nodesperelem(fes) == 10
      bfes = meshboundary(fes)
      @test nodesperelem(bfes) == 6
   true
end
end
using .mbas007
mbas007.test()

module mbas008
using FinEtools
using Test
function test()
 fens,fes = L2block(1.37, 3); # Mesh
 @test nodesperelem(fes) == 2
 bfes = meshboundary(fes)
 @test nodesperelem(bfes) == 1
 true
end
end
using .mbas008
mbas008.test()

module mbas009
using FinEtools
using Test
function test()
   fens,fes = Q8block(1.3, 3.1,3,2); # Mesh
   @test nodesperelem(fes) == 8
   bfes = meshboundary(fes)
   @test nodesperelem(bfes) == 3
   true
end
end
using .mbas009
mbas009.test()

module mbas010
using FinEtools
using FinEtools.AlgoBaseModule: fieldnorm
using Test
function test()
   fens,fes = Q8block(1.3, 3.1,13,12); # Mesh
   geom = NodalField(fens.xyz)
   femm = FEMMBase(IntegDomain(fes, GaussRule(2, 3)))
   regions = [FDataDict("femm"=>femm)]
   # x = NodalField(reshape(fens.xyz[:, 1], count(fens), 1))
   # targetfields = [x]
   c = NodalField(fill(1.0, count(fens), 1))
   targetfields = [c]
   modeldata = FDataDict("fens"=>fens, "regions"=>regions, "targetfields"=>targetfields, "geom"=>geom, "elementsize"=>3.1 / 12)
   @test abs(fieldnorm(modeldata) - sqrt(1.3 * 3.1)) <= 1.0e-5
   true
end
end
using .mbas010
mbas010.test()

module mbascs01
using FinEtools
using FinEtools.AlgoBaseModule: fieldnorm
using LinearAlgebra
using Test
function test()
    center = [0.0, 0.0, 0.0]
    ez = [0.0, 0.0, 1.0]
    v = rand(3)
    XYZ = reshape([0.2, 0.3, 0.4], 1, 3)
    tangents = rand(3, 3)
    rfcsmat = [0.5547001962252291 -0.8320502943378437 0.0; 0.8320502943378436 0.5547001962252293 0.0; 0.0 0.0 1.0]
    function compute!(csmatout::FFltMat, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt)
        # Cylindrical coordinate system
        xyz = XYZ[:] .- center
        xyz[3] = 0.0
        csmatout[:, 1] = xyz / norm(xyz)
        cross3!(v, ez, vec(csmatout[:, 1]))
        csmatout[:, 2] = v / norm(v)
        csmatout[:, 3] = ez
        return csmatout
    end
    csys = CSys(3, 3, compute!)
    updatecsmat!(csys, XYZ, tangents, 0)
    @test norm(csys.csmat -  rfcsmat) / norm(rfcsmat) <= 1.0e-6
    true
end
end
using .mbascs01
mbascs01.test()

module mbascs02
using FinEtools
using FinEtools.AlgoBaseModule: fieldnorm
using LinearAlgebra
using Test
function test()
    center = [0.0, 0.0, 0.0]
    ez = [0.0, 0.0, 1.0]
    v = rand(3)
    XYZ = reshape([0.2, 0.3, 0.4], 1, 3)
    tangents = rand(3, 3)
    rfcsmat = [0.5547001962252291 -0.8320502943378437 0.0; 0.8320502943378436 0.5547001962252293 0.0; 0.0 0.0 1.0]
    # function compute!(csmatout::FFltMat, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt)
    #     # Cylindrical coordinate system
    #     xyz = XYZ[:] .- center
    #     xyz[3] = 0.0
    #     csmatout[:, 1] = xyz / norm(xyz)
    #     cross3!(v, ez, vec(csmatout[:, 1]))
    #     csmatout[:, 2] = v / norm(v)
    #     csmatout[:, 3] = ez
    #     return csmatout
    # end
    csys = CSys(rfcsmat)
    updatecsmat!(csys, XYZ, tangents, 0)
    @test norm(csys.csmat -  rfcsmat) / norm(rfcsmat) <= 1.0e-6
    true
end
end
using .mbascs02
mbascs02.test()

module mbascs03
using FinEtools
using FinEtools.AlgoBaseModule: fieldnorm
using LinearAlgebra
using Test
function test()
    center = [0.0, 0.0, 0.0]
    ez = [0.0, 0.0, 1.0]
    v = rand(3)
    XYZ = reshape([0.2, 0.3, 0.4], 1, 3)
    tangents = rand(3, 3)
    rfcsmat = Matrix(1.0 * I, 3, 3)
    # function compute!(csmatout::FFltMat, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt)
    #     # Cylindrical coordinate system
    #     xyz = XYZ[:] .- center
    #     xyz[3] = 0.0
    #     csmatout[:, 1] = xyz / norm(xyz)
    #     cross3!(v, ez, vec(csmatout[:, 1]))
    #     csmatout[:, 2] = v / norm(v)
    #     csmatout[:, 3] = ez
    #     return csmatout
    # end
    csys = CSys(3)
    updatecsmat!(csys, XYZ, tangents, 0)
    @test norm(csys.csmat -  rfcsmat) / norm(rfcsmat) <= 1.0e-6
    true
end
end
using .mbascs03
mbascs03.test()

module mbascs04
using FinEtools
using FinEtools.AlgoBaseModule: fieldnorm
using LinearAlgebra
using Test
function test()
    XYZ = reshape([0.2, 0.3, 0.4], 1, 3)
    rfcsmat = Matrix(1.0 * I, 3, 3)
    tangents = rand(3, 3)
    csys = CSys(3, 3)
    updatecsmat!(csys, XYZ, tangents, 0)
    @test norm(csys.csmat -  rfcsmat) / norm(rfcsmat) <= 1.0e-6
    true
end
end
using .mbascs04
mbascs04.test()

module mbascs05
using FinEtools
using FinEtools.AlgoBaseModule: fieldnorm
using LinearAlgebra
using Test
function test()
    XYZ = reshape([0.2, 0.3, 0.4], 1, 3)
    rfcsmat = reshape([0.37139067635410367; 0.5570860145311555; 0.7427813527082073], 3, 1)
    tangents = reshape([0.2, 0.3, 0.4], 3, 1)
    csys = CSys(3, 1)
    updatecsmat!(csys, XYZ, tangents, 0)
    # @show csys.csmat
    @test norm(csys.csmat -  rfcsmat) / norm(rfcsmat) <= 1.0e-6
    true
end
end
using .mbascs05
mbascs05.test()

module mbascs06
using FinEtools
using FinEtools.AlgoBaseModule: fieldnorm
using LinearAlgebra
using Test
function test()
    XYZ = reshape([0.2, 0.3, 0.4], 1, 3)
    rfcsmat = reshape([0.37139067635410367; 0.5570860145311555; 0.7427813527082073], 3, 1)
    tangents = reshape([0.2 0.0; 0.0 0.5; 0.4 0.2], 3, 2)
    csys = CSys(3, 2)
    updatecsmat!(csys, XYZ, tangents, 0)
    # @show csys.csmat
    n1 = cross(vec(tangents[:, 1]), vec(tangents[:, 2]))
    n1 = n1 / norm(n1)
    # @show n1
    n2 = cross(vec(csys.csmat[:, 1]), vec(csys.csmat[:, 2]))
    n2 = n2 / norm(n2)
    # @show n2
    @test norm(n1 -  n2) / norm(n2) <= 1.0e-6
    true
end
end
using .mbascs06
mbascs06.test()

module mbascs06Alpha
using FinEtools
using FinEtools.MatrixUtilityModule: complete_lt!, add_mggt_ut_only!, add_gkgt_ut_only!, add_btdb_ut_only!, add_btsigma!, add_nnt_ut_only!
using LinearAlgebra
using Test
function test()
    let
        N = 8
        Ke::FFltMat, gradN::FFltMat, mult::FFlt = fill(0.0, N, N), rand(N, 2), 3.0
        add_mggt_ut_only!(Ke::FFltMat, gradN::FFltMat, mult::FFlt)
        complete_lt!(Ke::FFltMat)
        @test norm(Ke - (mult .* (gradN * gradN'))) / norm(Ke) <= 1.0e-9
    end

    let
        N = 8
        Ke::FFltMat, gradN::FFltMat, Jac_w::FFlt = fill(0.0, N, N), rand(N, 3), 0.33
        kappa_bar::FFltMat = rand(3, 3)
        kappa_bar = kappa_bar + kappa_bar'
        kappa_bargradNT::FFltMat = fill(0.0, 3, N)
        add_gkgt_ut_only!(Ke::FFltMat, gradN::FFltMat, Jac_w::FFlt,
            kappa_bar::FFltMat, kappa_bargradNT::FFltMat)
        complete_lt!(Ke::FFltMat)
        @test norm(Ke - (Jac_w .* (gradN * kappa_bar * gradN'))) / norm(Ke) <= 1.0e-9
    end

    let
        N = 12
        Ke::FFltMat, B::FFltMat, Jac_w::FFlt = fill(0.0, N, N), rand(3, N), 0.33
        D::FFltMat = rand(3, 3)
        D = D + D'
        DB::FFltMat = fill(0.0, 3, N)
        add_btdb_ut_only!(Ke::FFltMat, B::FFltMat, Jac_w::FFlt, D::FFltMat, DB::FFltMat)
        complete_lt!(Ke::FFltMat)
        @test norm(Ke - (Jac_w .* (B' * D * B))) / norm(Ke) <= 1.0e-9
    end

    let
        N = 16
        Fe::FFltVec, B::FFltMat, coefficient::FFlt, sigma::FFltVec = fill(0.0, N), rand(3, N), 0.33, rand(3)
        D::FFltMat = rand(3, 3)
        D = D + D'
        add_btsigma!(Fe::FFltVec, B::FFltMat, coefficient::FFlt, sigma::FFltVec)
        @test norm(Fe - (coefficient * (B' * sigma))) / norm(Fe) <= 1.0e-9
    end

    let
        M = 21
        Ke::FFltMat, Nn::FFltMat, Jac_w_coeff::FFlt = fill(0.0, M, M), rand(M, 1), 0.533
        add_nnt_ut_only!(Ke::FMat{FFlt}, Nn::FFltMat, Jac_w_coeff::FFlt)
        complete_lt!(Ke::FFltMat)
        @test norm(Ke - ((Nn*Nn')*(Jac_w_coeff))) / norm(Ke) <= 1.0e-9
    end

    true
end
end
using .mbascs06Alpha
mbascs06Alpha.test()

module msymmx
using FinEtools
using FinEtools.MatrixUtilityModule: symmetrize!
using LinearAlgebra
using Test
function test()
    a = rand(6, 6)
    b = deepcopy(a)
    symmetrize!(a)
    @test norm(a - (b + b')/2) / norm(a) <= 1.0e-9
    true
end
end
using .msymmx
msymmx.test()

module mgausr111
using FinEtools
using FinEtools.MeshExportModule: VTK
using LinearAlgebra
using Test
function test()
  fens,fes = H8spheren(3.1, 5)
  @test (count(fens), count(fes)) == (175, 108)
  # File = "mesh.vtk"
  # VTK.vtkexportmesh(File, fens, fes)
  geom  =  NodalField(fens.xyz)
  results = []
  for order in 1:10
      femm  =  FEMMBase(IntegDomain(fes, GaussRule(3, order)))
      V = integratefunction(femm, geom, (x) ->  1.0)
      push!(results, V)
  end
  @test norm(results - [15.06847962861544, 15.137585913524491, 15.137585913524477, 15.137585913524472, 15.137585913524491, 15.13758591352459, 15.137585913524601, 15.137585913524363, 15.137585913524278, 15.137585913524612]) <= 1.0e-9
  true
end
end
using .mgausr111
mgausr111.test()

module mtrapr111
using FinEtools
using FinEtools.MeshExportModule: VTK
using LinearAlgebra
using Test
function test()
  fens,fes = H8spheren(3.1, 5)
  @test (count(fens), count(fes)) == (175, 108)
  # File = "mesh.vtk"
  # VTK.vtkexportmesh(File, fens, fes)
  geom  =  NodalField(fens.xyz)
  femm  =  FEMMBase(IntegDomain(fes, TrapezoidalRule(3)))
  V = integratefunction(femm, geom, (x) ->  1.0)
  @test abs(V - 15.27579848334253) / V <= eps(V)
  true
end
end
using .mtrapr111
mtrapr111.test()

module mtrapr112
using FinEtools
using FinEtools.MeshExportModule: VTK
using LinearAlgebra
using Test
function test()
  rin::FFlt, rex::FFlt, nr::FInt, nc::FInt, Angl::FFlt, orientation::Symbol = 100.0, 200.0, 3, 6, pi/3, :a
  fens, fes = T3annulus(rin::FFlt, rex::FFlt, nr::FInt, nc::FInt, Angl::FFlt, orientation::
      Symbol)
  # File = "mesh.vtk"
  # VTK.vtkexportmesh(File, fens, fes)
  geom  =  NodalField(fens.xyz)
  femm  =  FEMMBase(IntegDomain(fes, TrapezoidalRule(2)))
  V1 = integratefunction(femm, geom, (x) ->  1.0)
  #
  femm  =  FEMMBase(IntegDomain(fes, GaussRule(2, 2)))
  V2 = integratefunction(femm, geom, (x) ->  1.0)
  @test abs(V1 - V2) / V1 <= 1000 * eps(V1)
  true
end
end
using .mtrapr112
mtrapr112.test()

module mfld1a1
using FinEtools
using LinearAlgebra
using Test

cleandchi() = deepcopy(FinEtools.NodalFieldModule.NodalField{Float64}(
    [0.0 0.0 0.0 -0.002084946198827584 0.0 0.0; 6.938893903907228e-18 1.0842021724855044e-19 -0.14078343955580946 -0.001474279595600101 -2.710505431213761e-20 1.0842021724855044e-19; 2.7755575615628914e-17 0.0 -0.19909784957735863 -8.605854744103691e-19 0.0 -5.421010862427522e-20; 2.7755575615628914e-17 1.3552527156068805e-20 -0.14078343955580952 0.001474279595600101 0.0 0.0; 0.0 0.0 0.0 0.0020849461988275853 0.0 0.0],
    [0 0 0 1 0 0; 2 3 4 5 6 7; 8 9 10 11 12 13; 14 15 16 17 18 19; 0 0 0 20 0 0],
    Bool[1 1 1 0 1 1; 0 0 0 0 0 0; 0 0 0 0 0 0; 0 0 0 0 0 0; 1 1 1 0 1 1],
    [0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0], 20)
)
cleandchiv() = deepcopy([-0.002084946198827584, 6.938893903907228e-18, 1.0842021724855044e-19, -0.14078343955580946, -0.001474279595600101, -2.710505431213761e-20, 1.0842021724855044e-19, 2.7755575615628914e-17, 0.0, -0.19909784957735863, -8.605854744103691e-19, 0.0, -5.421010862427522e-20, 2.7755575615628914e-17, 1.3552527156068805e-20, -0.14078343955580952, 0.001474279595600101, 0.0, 0.0, 0.0020849461988275853])
function test()
    dchi = cleandchi()
    v = gathersysvec(dchi)
    @test norm(v - cleandchiv()) / norm(cleandchiv()) <= 1.0e-6
    v[:] .= 0.0
    gathersysvec!(dchi, v)
    @test norm(v - cleandchiv()) / norm(cleandchiv()) <= 1.0e-6
    elv = fill(0.0, 12)
    gathervalues_asvec!(dchi, elv, [1, 2])
    @test norm(elv - [0.0, 0.0, 0.0, -0.002084946198827584, 0.0, 0.0, 6.938893903907228e-18, 1.0842021724855044e-19, -0.14078343955580946, -0.001474279595600101, -2.710505431213761e-20, 1.0842021724855044e-19]) / norm(elv) <= 1.0e-6
    gatherfixedvalues_asvec!(dchi, elv, [1, 2])
    norm(elv) <= 1.0e-15
    @test FinEtools.FieldModule.anyfixedvaluenz(dchi, [1, 2]) == false
    edn = fill(0, 12)
    gatherdofnums!(dchi, edn, [1, 2])
    @test norm(edn - vec([0 0 0 1 0 0; 2 3 4 5 6 7]')) == 0
    dchi = cleandchi()
    numberdofs!(dchi)
    @test norm(dchi.dofnums - cleandchi().dofnums) == 0

    dchi = cleandchi()
    setebc!(dchi, 1, true, 1, -0.013)
    @test dchi.is_fixed[1, 1] == true
    @test dchi.fixed_values[1, 1] == -0.013
    setebc!(dchi, 5, true, 4, 0.61)
    @test dchi.is_fixed[5, 4] == true
    @test dchi.fixed_values[5, 4] == 0.61

    dchi = cleandchi()
    setebc!(dchi, [1, 5], true, 2, [-0.013, 2.0])
    @test dchi.is_fixed[1, 2] == true
    @test dchi.fixed_values[1, 2] == -0.013
    @test dchi.is_fixed[5, 2] == true
    @test dchi.fixed_values[5, 2] == 2.0

    dchi = cleandchi()
    setebc!(dchi, [1, 5], true, 2, -10.0)
    @test dchi.is_fixed[1, 2] == true
    @test dchi.fixed_values[1, 2] == -10.0
    @test dchi.is_fixed[5, 2] == true
    @test dchi.fixed_values[5, 2] == -10.0

    dchi = cleandchi()
    setebc!(dchi, [1, 5], true, 2)
    @test dchi.is_fixed[1, 2] == true
    @test dchi.fixed_values[1, 2] == 0.0
    @test dchi.is_fixed[5, 2] == true
    @test dchi.fixed_values[5, 2] == 0.0

    dchi = cleandchi()
    setebc!(dchi, [1, 5], 2, [-10.0, 10.0])
    @test dchi.is_fixed[1, 2] == true
    @test dchi.fixed_values[1, 2] == -10.0
    @test dchi.is_fixed[5, 2] == true
    @test dchi.fixed_values[5, 2] == 10.0

    dchi = cleandchi()
    setebc!(dchi, [1, 5], 2)
    @test dchi.is_fixed[1, 2] == true
    @test dchi.fixed_values[1, 2] == 0.0
    @test dchi.is_fixed[5, 2] == true
    @test dchi.fixed_values[5, 2] == 0.0

    dchi = cleandchi()
    setebc!(dchi, [1, 5], true, [2, 3], -10.0)
    @test dchi.is_fixed[1, 2] == true
    @test dchi.fixed_values[1, 2] == -10.0
    @test dchi.is_fixed[5, 3] == true
    @test dchi.fixed_values[5, 3] == -10.0

    dchi = cleandchi()
    setebc!(dchi, [1, 5], true, [2, 3])
    @test dchi.is_fixed[1, 2] == true
    @test dchi.fixed_values[1, 2] == 0.0
    @test dchi.is_fixed[5, 3] == true
    @test dchi.fixed_values[5, 3] == 0.0

    dchi = cleandchi()
    setebc!(dchi, [1, 5])
    for i in [1, 5]
        @test any([!dchi.is_fixed[i, idx] for idx in 1:ndofs(dchi)]) == false
        @test any([dchi.fixed_values[i, idx] != 0.0 for idx in 1:ndofs(dchi)]) == false
    end

    dchi = cleandchi()
    setebc!(dchi, 3)
    for i in [3]
        @test any([!dchi.is_fixed[i, idx] for idx in 1:ndofs(dchi)]) == false
        @test any([dchi.fixed_values[i, idx] != 0.0 for idx in 1:ndofs(dchi)]) == false
    end

    dchi = cleandchi()
    setebc!(dchi)
    for i in 1:nents(dchi)
        @test any([dchi.is_fixed[i, idx] for idx in 1:ndofs(dchi)]) == false
        @test any([dchi.fixed_values[i, idx] != 0.0 for idx in 1:ndofs(dchi)]) == false
    end

    dchi = cleandchi()
    v = gathersysvec(dchi)
    dchi.values[:] .= 0.0
    scattersysvec!(dchi, v)
    v = gathersysvec(dchi)
    @test norm(v - cleandchiv()) / norm(cleandchiv()) <= 1.0e-6
    v[:] .= 0.0
    gathersysvec!(dchi, v)
    @test norm(v - cleandchiv()) / norm(cleandchiv()) <= 1.0e-6

    dchi = cleandchi()
    v = gathersysvec(dchi)
    dchi.values[:] .= 0.0
    incrscattersysvec!(dchi, v)
    incrscattersysvec!(dchi, v)
    incrscattersysvec!(dchi, v)
    v = gathersysvec(dchi)
    @test norm(v - 3 .* cleandchiv()) / norm(cleandchiv()) <= 1.0e-6

    dchi = cleandchi()
    setebc!(dchi)
    numberdofs!(dchi)
    dofnums, prescribedvalues = prescribeddofs(cleandchi(), dchi)
    @test norm(dofnums - [1, 2, 3, 5, 6, 25, 26, 27, 29, 30]) == 0
    @test norm(prescribedvalues - [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]) <= 1.0e-9
    true
end
end
using .mfld1a1
mfld1a1.test()

module mvass1
using FinEtools
using LinearAlgebra
using Test
function test()
    N = 30
    v = fill(0.0, N)
    a = SysvecAssembler()
    startassembly!(a,  N)
    vec, dofnums = rand(3), [1, 7, 2]
    assemble!(a, vec, dofnums)
    for (i, p) in zip(dofnums, vec)
        v[i] += p
    end
    vec, dofnums = rand(7), [29, 15, 1, 7, 3, 6, 2]
    assemble!(a, vec, dofnums)
    for (i, p) in zip(dofnums, vec)
        v[i] += p
    end
    w = makevector!(a)
    @test norm(v-w) / norm(w) <= 1.0e-9
end
end
using .mvass1
mvass1.test()

module mvass2
using FinEtools
using LinearAlgebra
using Test
function test()
    N = 30
    v = fill(0.0, N)
    a = SysvecAssembler()
    startassembly!(a,  N)
    vec, dofnums = rand(3), [1, 7, 2]
    assemble!(a, vec, dofnums)
    for (i, p) in zip(dofnums, vec)
        v[i] += p
    end
    vec, dofnums = rand(7), [29, 0, 1, 7, 3, 0, 2]
    assemble!(a, vec, dofnums)
    for (i, p) in zip(dofnums, vec)
        if i > 0
            v[i] += p
        end
    end
    w = makevector!(a)
    @test norm(v-w) / norm(w) <= 1.0e-9
end
end
using .mvass2
mvass2.test()

module mmhrzass1
using FinEtools
using LinearAlgebra
using Test
function test()
    N = 30
    elem_mat_dim = 8
    elem_mat_nmatrices = 3
    elmat = rand(elem_mat_dim, elem_mat_dim)
    elmat = (elmat + elmat')
    Mref = fill(0.0, N, N)
    add!(Mref, e, d) = begin
        em2 = sum(sum(e, dims = 1));
        dem2 = sum(diag(e));
        ffactor = em2/dem2
        for i in 1:length(d)
            Mref[d[i], d[i]] += ffactor*e[i, i]
        end
        Mref
    end
    a = SysmatAssemblerSparseHRZLumpingSymm()
    startassembly!(a, elem_mat_dim, 0, elem_mat_nmatrices, N, 0)
    dofnums = [10, 29, 15, 1, 7, 3, 6, 2]
    assemble!(a, elmat,  dofnums, dofnums)
    add!(Mref, elmat, dofnums)
    dofnums = [15, 1, 7, 11, 29, 3, 8, 4]
    assemble!(a, elmat,  dofnums, dofnums)
    add!(Mref, elmat, dofnums)
    dofnums = reshape([25, 1, 17, 13, 30, 13, 18, 14], elem_mat_dim, 1)
    assemble!(a, elmat,  dofnums, dofnums)
    add!(Mref, elmat, dofnums)
    M = makematrix!(a)
    s = 0.0
    for i in 1:N
        s += abs(M[i, i]-diag(Mref)[i])
    end
    @test s / norm(elmat) <= 1.0e-9
end
end
using .mmhrzass1
mmhrzass1.test()

module mgram1
using FinEtools
using Test
function test()
    W = 1.1;
    L = 12.;
    t =  0.32;
    nl, nt, nw = 2, 3, 4;

    fens,fes  = H8block(L,W,t, nl,nw,nt)
    geom  =  NodalField(fens.xyz)
    psi  =  NodalField(fill(1.0, count(fens), 1))
    numberdofs!(psi)

    femm  =  FEMMBase(IntegDomain(fes, GaussRule(3, 2)))
    G = innerproduct(femm, geom, psi)
    v = gathersysvec(psi)

    @test abs(v' * G * v - 4.224) / 4.224 <=  1.0e-5
    true
end
end
using .mgram1
mgram1.test()

module mdistributedl1
using FinEtools
using Test
function test()
    W = 1.1;
    L = 12.;
    t = 4.32;
    nl, nt, nw = 2, 3, 4;

    fens,fes  = H8block(L,W,t, nl,nw,nt)
    geom  =  NodalField(fens.xyz)
    psi  =  NodalField(fill(1.0, count(fens), 1))
    numberdofs!(psi)

    femm  =  FEMMBase(IntegDomain(fes, GaussRule(3, 2)))
    fi = ForceIntensity([11.0])
    F = distribloads(femm, geom, psi, fi, 3)
    @test abs(sum(F) - (*).(L,W,t, 11.0)) / 667 <=  1.0e-5
    fi = ForceIntensity(11.0)
    F = distribloads(femm, geom, psi, fi, 3)
    @test abs(sum(F) - (*).(L,W,t, 11.0)) / 667 <=  1.0e-5
true
end
end
using .mdistributedl1
mdistributedl1.test()

module m1djac1
using FinEtools
using Test
function test()
    Lx=1900.0;# length of the box, millimeters
    Ly=800.0; # length of the box, millimeters
    th = 21.0

    fens,fes = Q4block(Lx,Ly,3,2); # Mesh
    bfes = meshboundary(fes)
    bfemm  =  FEMMBase(IntegDomain(bfes, GaussRule(1, 2), th))

    geom  =  NodalField(fens.xyz)
    circumarea = integratefunction(bfemm, geom, (x) ->  1.0, 2)
    @test abs(circumarea - 2*(Lx+Ly)*th) <= 1.0e-9
true
end
end
using .m1djac1
m1djac1.test()

module m1djac2
using FinEtools
using Test
function test()
    Lx=1900.0;# length of the box, millimeters
    Ly=800.0; # length of the box, millimeters
    th = 21.0

    fens,fes = Q4block(Lx,Ly,3,2); # Mesh
    bfes = meshboundary(fes)
    bfemm  =  FEMMBase(IntegDomain(bfes, GaussRule(1, 2)))

    geom  =  NodalField(fens.xyz)
    circumference = integratefunction(bfemm, geom, (x) ->  1.0, 1)
    @test abs(circumference - 2*(Lx+Ly)) <= 1.0e-9
true
end
end
using .m1djac2
m1djac2.test()

module m1djac3
using FinEtools
using Test
function test()
    Lx=1900.0;# length of the box, millimeters
    Ly=800.0; # length of the box, millimeters
    th = 21.0

    fens,fes = Q4block(Lx,Ly,3,2); # Mesh
    bfes = meshboundary(fes)
    bfemm  =  FEMMBase(IntegDomain(bfes, GaussRule(1, 2), 7.0))

    geom  =  NodalField(fens.xyz)
    circumference = integratefunction(bfemm, geom, (x) ->  1.0, 3)
    @test abs(circumference/7 - 2*(Lx+Ly)) <= 1.0e-9
true
end
end
using .m1djac3
m1djac3.test()

module m2djac3
using FinEtools
using Test
function test()
    Lx=1900.0;# length of the box, millimeters
    Ly=800.0; # length of the box, millimeters
    th = 7.0

    fens,fes = Q4block(Lx,Ly,3,2); # Mesh
    bfes = meshboundary(fes)
    femm  =  FEMMBase(IntegDomain(fes, GaussRule(2, 2), th))

    geom  =  NodalField(fens.xyz)
    volume = integratefunction(femm, geom, (x) ->  1.0, 3)
    @test abs(volume - (Lx*Ly)*th) <= 1.0e-7
true
end
end
using .m2djac3
m2djac3.test()

module mgrad2mm
using FinEtools
using FinEtools.FESetModule: gradN!
using LinearAlgebra
using Test
function test()
  rho=1.21*1e-9;# mass density
  c =345.0*1000;# millimeters per second
  bulk= c^2*rho;
  Lx=1900.0;# length of the box, millimeters
  Ly=800.0; # length of the box, millimeters

  fens,fes = Q4block(Lx,Ly,3,2); # Mesh
  # show(fes.conn)
  gradN, gradNparams, redJ = fill(0.0, nodesperelem(fes), 2), bfundpar(fes, [0.57, 0.57]), Matrix(1.0*I, 2, 2)

  gradN!(fes, gradN::FFltMat, gradNparams::FFltMat, redJ::FFltMat)
  @test norm(gradN - [-0.10750000000000001 -0.10750000000000001; 0.10750000000000001 -0.39249999999999996; 0.39249999999999996 0.39249999999999996; -0.39249999999999996 0.10750000000000001]) <= 1.0e-6
end
end
using .mgrad2mm
mgrad2mm.test()

 module mgrad3mm
 using FinEtools
 using FinEtools.FESetModule: gradN!
 using LinearAlgebra
 using Test
 function test()
   rho=1.21*1e-9;# mass density
   c =345.0*1000;# millimeters per second
   bulk= c^2*rho;
   Lx=1900.0;# length of the box, millimeters
   Ly=800.0; # length of the box, millimeters
   Lz=800.0; # length of the box, millimeters

   fens,fes = H8block(Lx,Ly,Lz, 3,2, 4); # Mesh
   # show(fes.conn)
   gradN, gradNparams, redJ = fill(0.0, nodesperelem(fes), 3), bfundpar(fes, [0.57, 0.57, -0.57]), Matrix(1.0*I, 3, 3)

   gradN!(fes, gradN::FFltMat, gradNparams::FFltMat, redJ::FFltMat)
   # @show gradN
   @test norm(gradN - [-0.0843875 -0.0843875 -0.023112500000000005; 0.0843875 -0.30811249999999996 -0.0843875; 0.30811249999999996 0.30811249999999996 -0.30811249999999996; -0.30811249999999996 0.0843875 -0.0843875; -0.023112500000000005 -0.023112500000000005 0.023112500000000005; 0.023112500000000005 -0.0843875 0.0843875; 0.0843875 0.0843875 0.30811249999999996; -0.0843875 0.023112500000000005 0.0843875]) <= 1.0e-6
end
end
using .mgrad3mm
mgrad3mm.test()

module mgffip1
using FinEtools
using LinearAlgebra
using FinEtools.MeshExportModule: VTK
using FinEtools.AlgoBaseModule: evalconvergencestudy
import FinEtools.FEMMBaseModule: inspectintegpoints
using  FinEtools.MatrixUtilityModule: locjac!
using Test

Lx=1900.0;# length of the box, millimeters
Ly=800.0; # length of the box, millimeters
nx, ny = 10, 9
th = 7.0

f(x) = cos(0.93 * pi * x[1]/Lx) + sin(1.7 * pi * x[2]/Ly)
g(x) = -cos(0.93 * pi * x[1]/Lx) + 2*sin(1.7 * pi * x[2]/Ly)

function test()

    fens,fes = Q4block(Lx,Ly,nx, ny); # Mesh
    femm  =  FEMMBase(IntegDomain(fes, GaussRule(2, 2), th))
    psi  =  GeneralField(fill(0.0, count(fens), 2))
    for i in 1:count(fens)
        psi.values[i, 1] = f(fens.xyz[i, :])
        psi.values[i, 2] = g(fens.xyz[i, :])
    end
    numberdofs!(psi)
    @test psi.nfreedofs == 220

    psi  =  GeneralField(fill(0.0, count(fens)))
    for i in 1:count(fens)
        psi.values[i, 1] = f(fens.xyz[i, :])
    end
    numberdofs!(psi)
    @test psi.nfreedofs == 110
true
end
end
using .mgffip1
mgffip1.test()

module msurfn1
using FinEtools
using LinearAlgebra
using FinEtools.MeshExportModule: VTK
using Test

function test()
        ndimensions = 2
        tangents, fe_label = reshape([1.0; 1.0], 2, 1), 0
        n = SurfaceNormal(ndimensions::FInt)
        normal =  updatenormal!(n, [0.0 0.0 0.0], tangents::FFltMat, fe_label::FInt)
        @test norm(normal-[0.7071067811865475, -0.7071067811865475]) <= 1.0e-5

        ndimensions = 3
        tangents, fe_label = reshape([1.0 0.0; 1.0 0.0; 0.0 1.0], 3, 2), 0
        n = SurfaceNormal(ndimensions::FInt)
        normal =  updatenormal!(n, [0.0 0.0 0.0], tangents::FFltMat, fe_label::FInt)
        @test norm(normal-[0.7071067811865475, -0.7071067811865475, 0.0]) <= 1.0e-5
        true
    end
end
using .msurfn1
msurfn1.test()

module mnsimpr1
using FinEtools
using FinEtools.MeshExportModule
using Test
function test()
    xs = collect(linearspace(0.0, pi / 2, 5))
    fens, fes = L2blockx(xs)
    @test count(fes) == 4

    geom  =  NodalField(fens.xyz)
    femm  =  FEMMBase(IntegDomain(fes, NodalSimplexRule(3)))
    V1 = integratefunction(femm, geom, (x) ->  1.0)
    femm  =  FEMMBase(IntegDomain(fes, SimplexRule(3, 4)))
    V2 = integratefunction(femm, geom, (x) ->  1.0)
    @test abs(V1 - V2) / V1 <= 1.0e-5
    true
end
end
using .mnsimpr1
mnsimpr1.test()

module mnsimpr2
using FinEtools
using FinEtools.MeshExportModule
using Test
function test()
    xs = collect(linearspace(0.0, pi / 2, 5))
    ys = collect(linearspace(0.0, 1.0, 6).^2)
    fens, fes = T3blockx(xs, ys, :a)
    for i = 1:count(fens)
        a, y = fens.xyz[i,:]
        fens.xyz[i,1] = sin(a) * (y + 0.5)
        fens.xyz[i,2] = cos(a) * (y + 0.5)
    end
    @test count(fes) == 40

    geom  =  NodalField(fens.xyz)
    femm  =  FEMMBase(IntegDomain(fes, NodalSimplexRule(3)))
    V1 = integratefunction(femm, geom, (x) ->  1.0)
    femm  =  FEMMBase(IntegDomain(fes, SimplexRule(3, 4)))
    V2 = integratefunction(femm, geom, (x) ->  1.0)
    @test abs(V1 - V2) / V1 <= 1.0e-5
    true
end
end
using .mnsimpr2
mnsimpr2.test()

module mnsimpr3
using FinEtools
using FinEtools.MeshExportModule
using Test
function test()
    xs = collect(linearspace(0.0, pi / 2, 5))
    ys = collect(linearspace(0.0, 1.0, 6).^2)
    zs = collect(linearspace(0.0, 1.0, 3))
    fens, fes = T4blockx(xs, ys, zs, :a)
    for i = 1:count(fens)
        a, y, z = fens.xyz[i,:]
        fens.xyz[i,1] = sin(a) * (y + 0.5)
        fens.xyz[i,2] = cos(a) * (y + 0.5)
        fens.xyz[i,3] = z
    end
    @test count(fes) == 240

    geom  =  NodalField(fens.xyz)
    femm  =  FEMMBase(IntegDomain(fes, NodalSimplexRule(3)))
    V1 = integratefunction(femm, geom, (x) ->  1.0)
    femm  =  FEMMBase(IntegDomain(fes, SimplexRule(3, 4)))
    V2 = integratefunction(femm, geom, (x) ->  1.0)
    @test abs(V1 - V2) / V1 <= 1.0e-5
    true
end
end
using .mnsimpr3
mnsimpr3.test()

module mntpr1
using FinEtools
using FinEtools.MeshExportModule
using Test
function test()
    xs = collect(linearspace(0.0, pi / 2, 5))
    fens, fes = L2blockx(xs)
    @test count(fes) == 4

    geom  =  NodalField(fens.xyz)
    femm  =  FEMMBase(IntegDomain(fes, NodalTensorProductRule(1)))
    V1 = integratefunction(femm, geom, (x) ->  1.0)
    femm  =  FEMMBase(IntegDomain(fes, TrapezoidalRule(1)))
    V2 = integratefunction(femm, geom, (x) ->  1.0)
    @test abs(V1 - V2) / V1 <= 1.0e-5
    true
end
end
using .mntpr1
mntpr1.test()

module mntpr2
using FinEtools
using FinEtools.MeshExportModule
using Test
function test()
    xs = collect(linearspace(0.0, pi / 2, 5))
    ys = collect(linearspace(0.0, 1.0, 6).^2)
    fens, fes = Q4blockx(xs, ys)
    for i = 1:count(fens)
        a, y = fens.xyz[i,:]
        fens.xyz[i,1] = sin(a) * (y + 0.5)
        fens.xyz[i,2] = cos(a) * (y + 0.5)
    end
    @test count(fes) == 20

    geom  =  NodalField(fens.xyz)
    femm  =  FEMMBase(IntegDomain(fes, NodalTensorProductRule(2)))
    V1 = integratefunction(femm, geom, (x) ->  1.0)
    femm  =  FEMMBase(IntegDomain(fes, TrapezoidalRule(2)))
    V2 = integratefunction(femm, geom, (x) ->  1.0)
    @test abs(V1 - V2) / V1 <= 1.0e-5
    true
end
end
using .mntpr2
mntpr2.test()

module mntpr3
using FinEtools
using FinEtools.MeshExportModule
using Test
function test()
    xs = collect(linearspace(0.0, pi / 2, 5))
    ys = collect(linearspace(0.0, 1.0, 6).^2)
    zs = collect(linearspace(0.0, 1.0, 3))
    fens, fes = H8blockx(xs, ys, zs)
    for i = 1:count(fens)
        a, y, z = fens.xyz[i,:]
        fens.xyz[i,1] = sin(a) * (y + 0.5)
        fens.xyz[i,2] = cos(a) * (y + 0.5)
        fens.xyz[i,3] = z
    end
    @test count(fes) == 40

    geom  =  NodalField(fens.xyz)
    femm  =  FEMMBase(IntegDomain(fes, NodalTensorProductRule(3)))
    V1 = integratefunction(femm, geom, (x) ->  1.0)
    femm  =  FEMMBase(IntegDomain(fes, TrapezoidalRule(3)))
    V2 = integratefunction(femm, geom, (x) ->  1.0)
    @test abs(V1 - V2) / V1 <= 1.0e-5
    true
end
end
using .mntpr3
mntpr3.test()

module mr3m1
using FinEtools
using LinearAlgebra
using Test
function test()
    a = rand(3)
    R = rotmat3(a)
    R1 = fill(0.0, 3, 3)
    rotmat3!(R1, a)
    @test norm(R - R1) <= 1.0e-9
    @test norm(R*R1' - I) <= 1.0e-9
end
end
using .mr3m1
mr3m1.test()

module MMATD1
using FinEtools
using LinearAlgebra
using Test
mutable struct Mat<:AbstractMat
    mass_density:: Float64
end

function test()
    m = Mat(133.0)
    @test massdensity(m) == 133.0
end
end
using .MMATD1
MMATD1.test()

module mmassembl2
using FinEtools
using Test
import LinearAlgebra: norm, cholesky
function test()
    a = SysmatAssemblerSparse(0.0)
    startassembly!(a, 5, 5, 3, 7, 7)
    m = [0.24406   0.599773    0.833404  0.0420141
    0.786024  0.00206713  0.995379  0.780298
    0.845816  0.198459    0.355149  0.224996]
    assemble!(a, m, [1 7 5], [5 0 1 4])
    m = [0.146618  0.53471   0.614342    0.737833
 0.479719  0.41354   0.00760941  0.836455
 0.254868  0.476189  0.460794    0.00919633
 0.159064  0.261821  0.317078    0.77646
 0.643538  0.429817  0.59788     0.958909]
    assemble!(a, m, [2 3 1 0 5], [6 7 3 4])
    A = makematrix!(a)
    # @show Matrix(A)
    @test abs.(maximum([0.833404 0.0 0.460794 0.05121043 0.24406 0.254868 0.476189;
0.0 0.0 0.614342 0.737833 0.0 0.146618 0.53471; 0.0 0.0 0.00760941 0.836455 0.0 0.479719 0.41354; 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.355149 0.0 0.59788 1.183905 0.845816 0.643538 0.429817; 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.995379 0.0 0.0 0.780298 0.786024 0.0 0.0]  - A)) < 1.0e-5
    # @test abs(maximum(T_i)-1380.5883006341187) < 1.0e-3
end
end
using .mmassembl2
mmassembl2.test()

module mbas010x1
using FinEtools
using FinEtools.AlgoBaseModule: fieldnorm
using Test
function test()
    a, b, na, nb = 1.3, 3.1, 3, 2
    fens,fes = Q8block(a, b, na, nb); # Mesh
    geom = NodalField(fens.xyz)
    femm = FEMMBase(IntegDomain(fes, GaussRule(2, 3), true))

    V = integratefunction(femm, geom, (x) ->  1.0, 3)
    @test abs(V - pi * a^2 * b) / V <= eps(V)
    true
end
end
using .mbas010x1
mbas010x1.test()

module mbas010x2
using FinEtools
using FinEtools.AlgoBaseModule: fieldnorm
using Test
function test()
    a, b, na, nb = 1.3, 3.1, 3, 2
    fens,fes = Q8block(a, b, na, nb); # Mesh
    geom = NodalField(fens.xyz)
    bfes  =   meshboundary(fes);
    femm = FEMMBase(IntegDomain(bfes, GaussRule(1, 3), true))

    S = integratefunction(femm, geom, (x) ->  1.0, 2)
    @test abs(S - (2 * pi * a * b + 2 * pi * a^2)) / S <= eps(S)
    true
end
end
using .mbas010x2
mbas010x2.test()

module mbas010x3
using FinEtools
using FinEtools.AlgoBaseModule: fieldnorm
using Test
function test()
    a, b, na, nb = 1.3, 3.1, 3, 2
    fens,fes = T3block(a, b, na, nb); # Mesh
    geom = NodalField(fens.xyz)
    femm = FEMMBase(IntegDomain(fes, NodalSimplexRule(2), true))

    V = integratefunction(femm, geom, (x) ->  1.0, 3)
    @test abs(V - pi * a^2 * b) / V <= eps(V)
    true
end
end
using .mbas010x3
mbas010x3.test()

module mbas010x4
using FinEtools
using Test
function test()
    a, b, na, nb = 1.3, 3.1, 3, 2
    fens,fes = T3block(a, b, na, nb); # Mesh
    geom = NodalField(fens.xyz)
    bfes  =   meshboundary(fes);
    femm = FEMMBase(IntegDomain(bfes, NodalSimplexRule(1), true))

    S = integratefunction(femm, geom, (x) ->  1.0, 2)
    @test abs(S - (2 * pi * a * b + 2 * pi * a^2)) / S <= eps(S)
    true
end
end
using .mbas010x4
mbas010x4.test()

module massnomat1
using  SparseArrays: sparse, spzeros
using FinEtools
using FinEtools.MeshImportModule
using FinEtools.MeshExportModule.MESH
using Test
function test()
    a = SysmatAssemblerSparse(0.0, true)
    startassembly!(a, 5, 5, 3, 7, 7)
    m = [0.24406   0.599773    0.833404  0.0420141
        0.786024  0.00206713  0.995379  0.780298
        0.845816  0.198459    0.355149  0.224996]
    assemble!(a, m, [1 7 5], [5 2 1 4])
    m = [0.146618  0.53471   0.614342    0.737833
         0.479719  0.41354   0.00760941  0.836455
         0.254868  0.476189  0.460794    0.00919633
         0.159064  0.261821  0.317078    0.77646
         0.643538  0.429817  0.59788     0.958909]
    assemble!(a, m, [2 3 1 7 5], [6 7 3 4])
    A = makematrix!(a)
    @test A == spzeros(7, 7)
    a.nomatrixresult = false
    A = makematrix!(a)
    @test A[1, 1] ≈ 0.833404
    @test A[5, 1] ≈ 0.355149
    @test A[7, 6] ≈ 0.159064
    @test A[3, 7] ≈ 0.41354
    @test A[7, 7] ≈ 0.261821
    true
end
end
using .massnomat1
massnomat1.test()


text:
file = open(f->read(f, String), "day_13.txt")

ids = [x.captures[1] for x in eachmatch(r"(\d+)",file)]
ids = parse.(Int64,ids)

(start,ids) = (ids[1], ids[2:end])

calulate_delay(x) = x -((start) % x )

index_min = argmin(calulate_delay.(ids))

print("awnser1 : $(ids[index_min] * calulate_delay(ids[index_min])) ")

#we can solve question 2 as a linear algebra problem,
#for example: The input 2,3 can be solve as:
#start_time + 2X = c
#start_time + 3X = c -1

#There is problably a way to find the basis of the problem and resolve it;
#I dont know now, lets just learn more about linear algebra, RIP

a = [2 1;1 0]
b = [4;3]
a\b