text:
module LKJ

"""
A Julia package implementing the Lewandowski, Kurowicka and Joe (LKJ) probability distribution on the space of correlation matrices.
The probability density function of a d-dimensional LKJ-distribution (''d \\geq 2'') with parameter ''\\boldsymbol{\\eta} \\geq 1'' is:
```math
f(\\mathbf{A}; \\boldsymbol{\\eta}) = c_d |det(A)|^{\\eta-1}
```
where ''\\mathbf{A}'' is ''d \\times d'' correlation matrix and ''c_d'' is a normalization constant depending on the dimension 'd'.
"""
LKJ

using Distributions
using SpecialFunctions: lgamma
using LinearAlgebra
using ArgCheck

import Base: length, size
import Distributions: rand, logpdf, params

export LKJcorr, LKJcorrChol
export length, size, params, rand, logpdf, transform_free_to_chol, transform_chol_to_free, log_jacobian_det_free_to_chol
export link, invlink

## LKJ distribution for correlation matrices
struct LKJcorr <: ContinuousMatrixDistribution
    d::Int
    eta::Float64

    function LKJcorr(d::Int, eta::Float64)
        @argcheck (d >= 2) & (eta >= 1)
        new(d, eta)
    end
end

## LKJ distribution for cholesky factors of correlation matrices
struct LKJcorrChol <: ContinuousMatrixDistribution
    d::Int
    eta::Float64

    function LKJcorrChol(d::Int, eta::Float64)
        @argcheck (d >= 2) & (eta >= 1)
        new(d, eta)
    end
end

params(distr::Union{LKJcorr, LKJcorrChol}) = (distr.d, distr.eta)
length(distr::Union{LKJcorr, LKJcorrChol}) = distr.d * distr.d
size(distr::Union{LKJcorr, LKJcorrChol}) = (distr.d, distr.d)

function rand(distr::Union{LKJcorr, LKJcorrChol}, n::Int)
    res=[rand(distr)]
    for i in 2:n
        push!(res,rand(distr))
    end
    res
end

rand(distr::LKJcorr) = _rand(distr.d, distr.eta)
rand(distr::LKJcorrChol) = cholesky!(_rand(distr.d, distr.eta)).L

# function _rand(d::Int, eta::Float64)
#     beta = eta + (d-2)/2
#     u = rand(Beta(beta, beta), 1)[1]
#
#     stdn = Normal(0, 1)
#
#     R = ones(Float64, d, d)
#     R[1,2] = 2*u - 1
#     R[2,1] = R[1,2]
#
#     for k in 2:(d-1)
#         beta -= 0.5
#
#         ## sample point uniformly from (k-1)-sphere
#         u = rand(stdn, k)
#         u ./= sqrt(sum(u.^2))
#
#         w = sqrt(rand(Beta(k/2, beta), 1)[1]) .* u
#         z = cholesky(R[1:k, 1:k]).U' * w
#
#         R[1:k, k+1] = z
#         R[k+1, 1:k] = z'
#     end
#     R
# end

function _rand(d::Int, eta::Float64)
    beta = eta + (d-2)/2
    u = rand(Beta(beta, beta))

    stdn = Normal(0, 1)

    R = ones(Float64, d, d)
    R[1,2] = 2*u - 1
    R[2,1] = R[1,2]

    for k in 2:(d-1)
        beta -= 0.5

        ## sample point uniformly from (k-1)-sphere
        u = rand(stdn, k)
        u ./= sqrt(sum(u.^2))

        w = sqrt(rand(Beta(k/2, beta))) .* u
        z = cholesky(@view R[1:k, 1:k]).L * w

        R[1:k, k+1] .= z
        R[k+1, 1:k] .= z
    end
    Hermitian(R)
end

function iscorr(A::AbstractMatrix{T}) where T<:Real
    all(diag(A) .≈ one(T)) && isposdef(A)
end

function logpdf(distr::LKJcorr, x::AbstractMatrix{T}; norm=false) where T<:Real
    !(iscorr(x)) && return -T(Inf)

    n, m = size(x)

    if n != distr.d
        throw(ArgumentError("Number of rows/columns ($n) does not match dimension of distribution ($(distr.d))"))
    end

    c = norm ? log_normalizing_const(distr.d, distr.eta) : 0

    (distr.eta - 1)*log(det(x)) + c
end

function logpdf(distr::LKJcorrChol, x::LowerTriangular; norm=false)
    !(iscorr(x*x')) && return -Inf

    n, m = size(x)

    if n != m
        throw(ArgumentError("x is not square"))
    end

    if n != distr.d
        throw(ArgumentError("Number of rows/columns ($n) does not match dimension of distribution ($(distr.d))"))
    end

    c = norm ? log_normalizing_const(distr.d, distr.eta) : 0

    sum((distr.d .- (1:distr.d) .+ 2*distr.eta .- 2).*log.(diag(x))) + c
end

function log_normalizing_const(d::Int, eta::T) where T <: Real
    c = zero(T)

    if eta == one(T)
        for k in 1:floor(Int, (d-1)/2)
            c += lgamma(2*k)
        end

        if iseven(d)
            c += d*(d-2)/4 * log(pi) + (3*d^2 - 4*d) / 4*log(2) + d*lgamma(d/2) - (d-1)*lgamma(d)
        else
            c += (d^2 - 1)/4 * log(pi) - (d-1)^2 / 4*log(2) - (d-1)*lgamma((d+1)/2)
        end
    else
        c = (2*eta + d - 3)*log(2) + 2*lgamma(eta + d/2 - 1) - lgamma(2*eta + d - 2) - (d-2)*lgamma(eta + (d-1)/2)

        for k in 2:(d-1)
            c += k/2 * log(pi) + lgamma(eta + (d-1-k)/2)
        end
    end
    c
end

## Stan Manual (version 2.17) Section 35.11. Cholesky Factors of Correlation Matrices
## see https://github.com/stan-dev/math/blob/develop/stan/math/prim/mat/fun/cholesky_corr_constrain.hpp
## y = vector of unconstrained values
function transform_free_to_chol(y::AbstractVector{T}, K) where T
    z = tanh.(y)
    w = zeros(T, K, K)
    w[1,1] = 1.0

    k = 1
    for i in 2:K
        w[i,1] = z[k]
        k += 1
        sum_sqs = w[i,1]^2
        for j in 2:(i-1)
            if sum_sqs > 1.0
                w[i,j] = 0.0
            else
                w[i,j] = z[k] * sqrt(1.0 - sum_sqs)
            end
            k += 1
            sum_sqs += w[i,j]^2
        end

        if sum_sqs > 1.0
            w[i,i] = 0.0
        else
            w[i,i] = sqrt(1.0 - sum_sqs)
        end
    end

    LowerTriangular(w)
end

function invlink(distr::LKJcorrChol,y::AbstractArray)
    K=distr.d
    k_choose_2=K * (K-1) / 2
    length(y) != k_choose_2 && throw("distribution size is not compatible with vector size")
    transform_free_to_chol(y, K)
end

## see https://github.com/stan-dev/math/blob/develop/stan/math/prim/mat/fun/cholesky_corr_free.hpp
## w = Cholesky factor of correlation matrix
function transform_chol_to_free(w::LowerTriangular)
    K = size(w,1)
    y = zeros(eltype(w), Int(K*(K-1)/2))

    k = 1
    for i in 2:K
        y[k] = atanh(w[i,1])
        k += 1
        sum_sqs = w[i,1]^2

        for j in 2:(i-1)
            y[k] = atanh(w[i,j] / sqrt(1.0 - sum_sqs))
            k += 1
            sum_sqs += w[i,j]^2
        end
    end
    y
end

link(distr::LKJcorrChol,x::AbstractArray)=transform_chol_to_free(LowerTriangular(x))

## K = (1 + sqrt(1 + 8*length(y))) / 2

"""
    log_jacobian_det_free_to_chol(y, K)

Log absolute determinant of Jacobian of transformation from unconstrained vector (length K*(K-1)/2) to cholesky factor.

# Examples
```jldoctest
julia> log_jacobian_det_free_to_chol([1.0, 2.0, 3.0], 3)
-9.461226912195217
```
"""
function log_jacobian_det_free_to_chol(y, K)
    w = transform_free_to_chol(y, K)
    log_jacobian_det_free_to_chol(y, K, w)
end

function log_jacobian_det_free_to_chol(y, K, w::LowerTriangular)
    res = -2*sum(log.(cosh.(y)))

    for i in 2:K
        sum_sq = 0.0
        for j in 1:(i-1)
            ##res += 0.5*log(1 - sum(w[i,1:(j-1)].^2))
            res += 0.5*log(1 - sum_sq)
            sum_sq += w[i,j]^2
        end
    end
    res
end

end


text:
#=
  Copyright (c) 2021
  https://github.com/Ismael-VC/Varvara.jl/blob/main/CONTRIBUTORS.md

  Permission to use, copy, modify, and distribute this software for any
  purpose with or without fee is hereby granted, provided that the above
  copyright notice and this permission notice appear in all copies.

  THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
  WITH REGARD TO THIS SOFTWARE.
=#

module UxnASM

using OffsetArrays: OffsetMatrix
using ..Uxn: ByteVector

const TRIM = 0x0100
const LENGTH = 0x10000
const LABELS = 512
const MACROS = 256
const ByteMatrix{UInt8} = OffsetMatrix{UInt8,Matrix{UInt8}}

ByteMatrix(dim₁::UInt8, dim₂::UInt8)::ByteMatrix{UInt8} = OffsetMatrix(zeros(UInt8, dim₁, dim₂), 0:dim₁-1, 0:dim₂-1)

struct Macro
	name::String  # ByteVector
  items::ByteMatrix
	len::UInt8
end

typedef struct {
	char name[64];
	Uint16 addr, refs;
} Label;

typedef struct {
	Uint8 data[LENGTH];
	Uint16 ptr, length, llen, mlen;
	Label labels[LABELS];
	Macro macros[MACROS];
} Program;

Program p;
static Uint16 addr = 0;

/* clang-format off */

static char ops[][4] = {
	"LIT", "INC", "POP", "DUP", "NIP", "SWP", "OVR", "ROT",
	"EQU", "NEQ", "GTH", "LTH", "JMP", "JCN", "JSR", "STH",
	"LDZ", "STZ", "LDR", "STR", "LDA", "STA", "DEI", "DEO",
	"ADD", "SUB", "MUL", "DIV", "AND", "ORA", "EOR", "SFT"
};

static int   cpos(char *s, char a){ int i = 0; char c; while((c = s[i++])) if(c == a) return i; return -1; }
static int   scmp(char *a, char *b, int len) { int i = 0; while(a[i] == b[i]) if(!a[i] || ++i >= len) return 1; return 0; } /* string compare */
static int   sihx(char *s) { int i = 0; char c; while((c = s[i++])) if(!(c >= '0' && c <= '9') && !(c >= 'a' && c <= 'f')) return 0; return i > 1; } /* string is hexadecimal */
static int   shex(char *s) { int n = 0, i = 0; char c; while((c = s[i++])) if(c >= '0' && c <= '9') n = n * 16 + (c - '0'); else if(c >= 'a' && c <= 'f') n = n * 16 + 10 + (c - 'a'); return n; } /* string to num */
static int   slen(char *s) { int i = 0; while(s[i]) ++i; return i; } /* string length */
static char *scpy(char *src, char *dst, int len) { int i = 0; while((dst[i] = src[i]) && i < len - 2) i++; dst[i + 1] = '\0'; return dst; } /* string copy */
static char *scat(char *dst, const char *src) { char *ptr = dst + slen(dst); while(*src) *ptr++ = *src++; *ptr = '\0'; return dst; } /* string cat */

#pragma mark - Helpers

/* clang-format on */

#pragma mark - I/O

static Macro *
findmacro(char *name)
{
	int i;
	for(i = 0; i < p.mlen; ++i)
		if(scmp(p.macros[i].name, name, 64))
			return &p.macros[i];
	return NULL;
}

static Label *
findlabel(char *name)
{
	int i;
	for(i = 0; i < p.llen; ++i)
		if(scmp(p.labels[i].name, name, 64))
			return &p.labels[i];
	return NULL;
}

static Uint8
findopcode(char *s)
{
	int i;
	for(i = 0; i < 0x20; ++i) {
		int m = 0;
		if(!scmp(ops[i], s, 3))
			continue;
		if(!i) i |= (1 << 7); /* force keep for LIT */
		while(s[3 + m]) {
			if(s[3 + m] == '2')
				i |= (1 << 5); /* mode: short */
			else if(s[3 + m] == 'r')
				i |= (1 << 6); /* mode: return */
			else if(s[3 + m] == 'k')
				i |= (1 << 7); /* mode: keep */
			else
				return 0; /* failed to match */
			m++;
		}
		return i;
	}
	return 0;
}

static void
pushbyte(Uint8 b, int lit)
{
	if(lit) pushbyte(findopcode("LIT"), 0);
	p.data[p.ptr++] = b;
	p.length = p.ptr;
}

static void
pushshort(Uint16 s, int lit)
{
	if(lit) pushbyte(findopcode("LIT2"), 0);
	pushbyte((s >> 8) & 0xff, 0);
	pushbyte(s & 0xff, 0);
}

static void
pushword(char *s)
{
	int i = 0;
	char c;
	while((c = s[i++])) pushbyte(c, 0);
}

static char *
sublabel(char *src, char *scope, char *name)
{
	return scat(scat(scpy(scope, src, 64), "/"), name);
}

#pragma mark - Parser

static int
error(char *name, char *msg)
{
	fprintf(stderr, "%s: %s\n", name, msg);
	return 0;
}

static int
makemacro(char *name, FILE *f)
{
	Macro *m;
	char word[64];
	if(findmacro(name))
		return error("Macro duplicate", name);
	if(sihx(name) && slen(name) % 2 == 0)
		return error("Macro name is hex number", name);
	if(findopcode(name) || scmp(name, "BRK", 4) || !slen(name) || scmp(name, "include", 8))
		return error("Macro name is invalid", name);
	if(p.mlen == MACROS)
		return error("Too many macros", name);
	m = &p.macros[p.mlen++];
	scpy(name, m->name, 64);
	while(fscanf(f, "%63s", word) == 1) {
		if(word[0] == '{') continue;
		if(word[0] == '}') break;
		if(m->len > 64)
			return error("Macro too large", name);
		scpy(word, m->items[m->len++], 64);
	}
	return 1;
}

static int
makelabel(char *name)
{
	Label *l;
	if(findlabel(name))
		return error("Label duplicate", name);
	if(sihx(name) && slen(name) % 2 == 0)
		return error("Label name is hex number", name);
	if(findopcode(name) || scmp(name, "BRK", 4) || !slen(name))
		return error("Label name is invalid", name);
	if(p.llen == LABELS)
		return error("Too many labels", name);
	l = &p.labels[p.llen++];
	l->addr = addr;
	l->refs = 0;
	scpy(name, l->name, 64);
	return 1;
}

static int
addref(Label *l, Uint8 rel)
{
	if(rel) {
		int pos = cpos(l->name, '/');
		if(pos != -1) {
			char parent[64];
			Label *rl = findlabel(scpy(l->name, parent, pos));
			if(rl)
				++rl->refs;
		}
	}
	return ++l->refs;
}

static int
skipblock(char *w, int *cap, char a, char b)
{
	if(w[0] == b) {
		*cap = 0;
		return 1;
	}
	if(w[0] == a) *cap = 1;
	if(*cap) return 1;
	return 0;
}

static int
walktoken(char *w)
{
	Macro *m;
	if(findopcode(w) || scmp(w, "BRK", 4))
		return 1;
	switch(w[0]) {
	case '[': return 0;
	case ']': return 0;
	case '\'': return 1;
	case '.': return 2; /* zero-page: LIT addr-lb */
	case ',': return 2; /* relative:  LIT addr-rel */
	case ':': return 2; /* absolute:      addr-hb addr-lb */
	case ';': return 3; /* absolute:  LIT addr-hb addr-lb */
	case '$': return shex(w + 1);
	case '#': return slen(w + 1) == 4 ? 3 : 2;
	case '"': return slen(w + 1);
	}
	if((m = findmacro(w))) {
		int i, res = 0;
		for(i = 0; i < m->len; ++i)
			res += walktoken(m->items[i]);
		return res;
	}
	if(sihx(w) && slen(w) == 2)
		return 1;
	else if(sihx(w) && slen(w) == 4)
		return 2;
	return error("Invalid token", w);
}

static int
parsetoken(char *w)
{
	Label *l;
	Macro *m;
	if(w[0] == '.' && (l = findlabel(w + 1))) { /* zero-page */
		if(l->addr > 0xff)
			return error("Address is not in zero page", w);
		pushbyte(l->addr, 1);
		return addref(l, 1);
	} else if(w[0] == ',' && (l = findlabel(w + 1))) { /* relative */
		int off = l->addr - p.ptr - 3;
		if(off < -126 || off > 126)
			return error("Address is too far", w);
		pushbyte((Sint8)off, 1);
		return addref(l, 0);
	} else if(w[0] == ':' && (l = findlabel(w + 1))) { /* raw */
		pushshort(l->addr, 0);
		return addref(l, 1);
	} else if(w[0] == ';' && (l = findlabel(w + 1))) { /* absolute */
		pushshort(l->addr, 1);
		return addref(l, 1);
	} else if(findopcode(w) || scmp(w, "BRK", 4)) { /* opcode */
		pushbyte(findopcode(w), 0);
		return 1;
	} else if(w[0] == '"') { /* string */
		pushword(w + 1);
		return 1;
	} else if(w[0] == '\'') { /* char */
		pushbyte((Uint8)w[1], 0);
		return 1;
	} else if(w[0] == '#') { /* immediate */
		if(slen(w + 1) == 1)
			pushbyte((Uint8)w[1], 1);
		if(sihx(w + 1) && slen(w + 1) == 2)
			pushbyte(shex(w + 1), 1);
		else if(sihx(w + 1) && slen(w + 1) == 4)
			pushshort(shex(w + 1), 1);
		else
			return error("Invalid hexadecimal literal", w);
		return 1;
	} else if(sihx(w)) { /* raw */
		if(slen(w) == 2)
			pushbyte(shex(w), 0);
		else if(slen(w) == 4)
			pushshort(shex(w), 0);
		else
			return error("Invalid hexadecimal value", w);
		return 1;
	} else if((m = findmacro(w))) {
		int i;
		for(i = 0; i < m->len; ++i)
			if(!parsetoken(m->items[i]))
				return error("Invalid macro", m->name);
		return 1;
	}
	return error("Invalid token", w);
}

static int
doinclude(FILE *f, int (*pass)(FILE *))
{
	char word[64];
	FILE *finc;
	int ret;
	if(fscanf(f, "%63s", word) != 1)
		return error("End of input", "include");
	if(!(finc = fopen(word, "r")))
		return error("Include failed to open", word);
	ret = pass(finc);
	fclose(finc);
	return ret;
}

static int
pass1(FILE *f)
{
	int ccmnt = 0;
	char w[64], scope[64], subw[64];
	while(fscanf(f, "%63s", w) == 1) {
		if(skipblock(w, &ccmnt, '(', ')')) continue;
		if(slen(w) >= 63)
			return error("Pass 1 - Invalid token", w);
		if(w[0] == '|') {
			if(!sihx(w + 1))
				return error("Pass 1 - Invalid padding", w);
			addr = shex(w + 1);
		} else if(w[0] == '%') {
			if(!makemacro(w + 1, f))
				return error("Pass 1 - Invalid macro", w);
		} else if(w[0] == '@') {
			if(!makelabel(w + 1))
				return error("Pass 1 - Invalid label", w);
			scpy(w + 1, scope, 64);
		} else if(w[0] == '&') {
			if(!makelabel(sublabel(subw, scope, w + 1)))
				return error("Pass 1 - Invalid sublabel", w);
		} else if(scmp(w, "include", 8)) {
			if(!doinclude(f, pass1))
				return 0;
		} else if(sihx(w))
			addr += slen(w) / 2;
		else
			addr += walktoken(w);
	}
	rewind(f);
	return 1;
}

static int
pass2(FILE *f)
{
	int ccmnt = 0, cmacr = 0;
	char w[64], scope[64], subw[64];
	while(fscanf(f, "%63s", w) == 1) {
		if(w[0] == '%') continue;
		if(w[0] == '&') continue;
		if(w[0] == '[') continue;
		if(w[0] == ']') continue;
		if(skipblock(w, &ccmnt, '(', ')')) continue;
		if(skipblock(w, &cmacr, '{', '}')) continue;
		if(w[0] == '|') {
			if(p.length && (Uint16)shex(w + 1) < p.ptr)
				return error("Pass 2 - Memory overwrite", w);
			p.ptr = shex(w + 1);
			continue;
		} else if(w[0] == '$') {
			if(p.length && (Uint16)(p.ptr + shex(w + 1)) < p.ptr)
				return error("Pass 2 - Memory overwrite", w);
			p.ptr += shex(w + 1);
			continue;
		} else if(w[0] == '@') {
			scpy(w + 1, scope, 64);
			continue;
		} else if(scmp(w, "include", 8)) {
			if(!doinclude(f, pass2))
				return 0;
			continue;
		}
		if(w[1] == '&' && (w[0] == '.' || w[0] == ',' || w[0] == ';' || w[0] == ':'))
			scpy(sublabel(subw, scope, w + 2), w + 1, 64);
		if(!parsetoken(w))
			return error("Pass 2 - Unknown label", w);
	}
	return 1;
}

static void
cleanup(char *filename)
{
	int i;
	for(i = 0; i < p.llen; ++i)
		if(p.labels[i].name[0] >= 'A' && p.labels[i].name[0] <= 'Z')
			continue; /* Ignore capitalized labels(devices) */
		else if(!p.labels[i].refs)
			fprintf(stderr, "--- Unused label: %s\n", p.labels[i].name);
	fprintf(stderr, "Assembled %s in %.2fkb(%.2f%% used), %d labels, %d macros.\n", filename, (p.length - TRIM) / 1024.0, p.length / 652.80, p.llen, p.mlen);
}

int
main(int argc, char *argv[])
{
	FILE *src, *dst;
	if(argc < 3)
		return !error("usage", "input.tal output.rom");
	if(!(src = fopen(argv[1], "r")))
		return !error("Invalid Input", argv[1]);
	if(!pass1(src) || !pass2(src))
		return !error("Assembly", "Failed to assemble rom.");
	if(!(dst = fopen(argv[2], "wb")))
		return !error("Invalid Output", argv[2]);
	fwrite(p.data + TRIM, p.length - TRIM, 1, dst);
	cleanup(argv[2]);
	return 0;
}
=#

end  # module

ext: jl


text:
using JuliaWeBWorK
using Test

@testset "JuliaWeBWorK.jl" begin
    # Write your tests here.
end


text:
function repl_startdebugger_request(conn, params::String, crashreporting_pipename)
    hideprompt() do
        debug_pipename = params
        try
            @debug "Trying to connect to debug adapter."
            socket = Sockets.connect(debug_pipename)
            try
                DebugAdapter.startdebug(socket, (err, bt)->global_err_handler(err, bt, crashreporting_pipename, "Debugger"))
            finally
                close(socket)
            end
        catch err
            global_err_handler(err, catch_backtrace(), crashreporting_pipename, "Debugger")
        end
    end
end

function remove_lln!(ex::Expr)
    for i in length(ex.args):-1:1
        if ex.args[i] isa LineNumberNode
            deleteat!(ex.args, i)
        elseif ex.args[i] isa Expr
            remove_lln!(ex.args[i])
        end
    end
end

macro enter(command)
    remove_lln!(command)
    :(JSONRPC.send_notification(conn_endpoint[], "debugger/enter", $(string(command))))
end

macro run(command)
    remove_lln!(command)
    :(JSONRPC.send_notification(conn_endpoint[], "debugger/run", $(string(command))))
end

max_stars_repo_name: mortenpi/julia-vscode

max_issues_repo_name: mortenpi/julia-vscode

max_forks_repo_name: mortenpi/julia-vscode


text:
module Metalhead

using Flux
using Flux: outputsize, Zygote
using Functors
using BSON
using Artifacts, LazyArtifacts

import Functors

# Models
include("utilities.jl")
include("alexnet.jl")
include("vgg.jl")
include("resnet.jl")
include("googlenet.jl")
include("inception.jl")
include("squeezenet.jl")
include("densenet.jl")

export  AlexNet,
        VGG, VGG11, VGG13, VGG16, VGG19,
        ResNet, ResNet18, ResNet34, ResNet50, ResNet101, ResNet152,
        GoogLeNet, Inception3, SqueezeNet,
        DenseNet, DenseNet121, DenseNet161, DenseNet169, DenseNet201

# use Flux._big_show to pretty print large models
for T in (:AlexNet, :VGG, :ResNet, :GoogLeNet, :Inception3, :SqueezeNet, :DenseNet)
  @eval Base.show(io::IO, ::MIME"text/plain", model::$T) = _maybe_big_show(io, model)
end

end # module


text:
# This file was generated, do not modify it. # hide
fit!(classif, rows=train)
ŷ = predict_mode(classif, rows=test)
accuracy(ŷ, y[test]) |> r3

max_issues_repo_name: alan-turing-institute/MLJTutorials

max_forks_repo_name: alan-turing-institute/MLJTutorials


text:
module TestPursesBasic

using Test

using Purses

@testset "Basic functionality" begin
    x = 1.5
    purse = Purse(x)
    @test Purses.value(purse) === Purses.value(x) === x
    @test Purses.cache(purse) === ()
    @test_throws BoundsError Purses.cache(purse, 1)
    @test Purses.cache_signature(typeof(purse)) === Tuple{}

    x = 1.5
    purse = Purse(x, inv, sqrt, cbrt)
    @test Purses.value(purse) === x
    @test Purses.cache(purse) === (inv(x), sqrt(x), cbrt(x))
    @test Purses.cache_signature(typeof(purse)) === Tuple{typeof(inv), typeof(sqrt), typeof(cbrt)}
    @test inv(purse) === inv(x)
    @test sqrt(purse) === sqrt(x)
    @test cbrt(purse) === cbrt(x)
    @test sin(purse) === sin(x)

    x = 1.5
    y = 2.0
    purse = Purse(x, inv, sqrt => y, cbrt)
    @test Purses.value(purse) == x
    @test sqrt(purse) != sqrt(x)
    @test sqrt(purse) == y
end

@testset "Conversion" begin
    x = 1.5
    y = 1.0
    purse = Purse(x)
    @test convert(typeof(purse), purse) === purse
    @test convert(Any, purse) === purse
    @test convert(Float64, purse) === x
    @test convert(Purse{Float32}, purse) === Purse(Float32(x))
    @test convert(Purse{typeof(x),Tuple{typeof(+)}}, purse) === convert(Purse{typeof(x),Tuple{typeof(+)}}, x) === Purse(x, +)

    @test Purse{Float64,Tuple{typeof(-)},Tuple{Float64}}(x) === Purse(x, -)
    @test_throws InexactError Purse{Float64,Tuple{typeof(-)},Tuple{Int64}}(x)
    @test Purses.cache(Purse{Float64,Tuple{typeof(-)},Tuple{Int64}}(y), 1) isa Int64

    purses = Purse{Float64,Tuple{typeof(-)},Tuple{Float64}}[1.0, 2.0, 8.5]
    @test purses == map(x -> Purse(x, -), [1.0, 2.0, 8.5])
end

@testset "Function registration" begin
    my_inc(x) = x + 1

    @test Purses.register!(my_inc) == (my_inc,)

    x = 1.5
    purse = Purse(x, my_inc)
    @test Purses.cache_signature(typeof(purse)) === Tuple{typeof(my_inc)}
    @test my_inc(purse) === my_inc(x)

    @test Purses.register!(my_inc; force=true) == (my_inc,)
    @test my_inc(purse) === my_inc(x)

    x = 1.5
    y = 2.0
    purse = Purse(x, my_inc => y)
    @test my_inc(purse) !== my_inc(x)
    @test my_inc(purse) === y

    my_cos(x) = cos(x)
    purse_cos = Purse(x, cos)
    purse_my_cos = Purse(x, my_cos)
    @test my_cos(purse_cos) == my_cos(purse_my_cos) == cos(purse_cos) == cos(purse_my_cos)
end

end # module


text:
@ghdef mutable struct Gist
    url::Union{HTTP.URI, Nothing}
    forks_url::Union{HTTP.URI, Nothing}
    commits_url::Union{HTTP.URI, Nothing}
    id::Union{String, Nothing}
    description::Union{String, Nothing}
    public::Union{Bool, Nothing}
    owner::Union{Owner, Nothing}
    user::Union{Owner, Nothing}
    truncated::Union{Bool, Nothing}
    comments::Union{Int, Nothing}
    comments_url::Union{HTTP.URI, Nothing}
    html_url::Union{HTTP.URI, Nothing}
    git_pull_url::Union{HTTP.URI, Nothing}
    git_push_url::Union{HTTP.URI, Nothing}
    created_at::Union{Dates.DateTime, Nothing}
    updated_at::Union{Dates.DateTime, Nothing}
    forks::Union{Vector{Gist}, Nothing}
    files::Union{Dict, Nothing}
    history::Union{Vector{Dict}, Nothing}
end

Gist(id::AbstractString) = Gist(Dict("id" => id))

namefield(gist::Gist) = gist.id

###############
# API Methods #
###############

# creating #
#----------#

@api_default gist(api::GitHubAPI, gist_obj::Gist; options...) = gist(api::GitHubAPI, name(gist_obj); options...)

@api_default function gist(api::GitHubAPI, gist_obj, sha = ""; options...)
    !isempty(sha) && (sha = "/" * sha)
    result = gh_get_json(api, "/gists/$(name(gist_obj))$sha"; options...)
    g = Gist(result)
end

@api_default function gists(api::GitHubAPI, owner; options...)
    results, page_data = gh_get_paged_json(api, "/users/$(name(owner))/gists"; options...)
    map(Gist, results), page_data
end

@api_default function gists(api::GitHubAPI; options...)
    results, page_data = gh_get_paged_json(api, "/gists/public"; options...)
    return map(Gist, results), page_data
end

# modifying #
#-----------#

@api_default create_gist(api::GitHubAPI; options...) = Gist(gh_post_json(api, "/gists"; options...))
@api_default edit_gist(api::GitHubAPI, gist; options...) = Gist(gh_patch_json(api, "/gists/$(name(gist))"; options...))
@api_default delete_gist(api::GitHubAPI, gist; options...) = gh_delete(api, "/gists/$(name(gist))"; options...)

# stars #
#------#

@api_default star_gist(api::GitHubAPI, gist; options...) = gh_put(api, "/gists/$(name(gist))/star"; options...)
@api_default unstar_gist(api::GitHubAPI, gist; options...) = gh_delete(api, "/gists/$(name(gist))/star"; options...)

@api_default function starred_gists(api::GitHubAPI; options...)
    results, page_data = gh_get_paged_json(api, "/gists/starred"; options...)
    return map(Gist, results), page_data
end

# forks #
#-------#

@api_default create_gist_fork(api::GitHubAPI, gist::Gist; options...) = Gist(gh_post_json(api, "/gists/$(name(gist))/forks"; options...))

@api_default function gist_forks(api::GitHubAPI, gist; options...)
    results, page_data = gh_get_paged_json(api, "/gists/$(name(gist))/forks"; options...)
    return map(Gist, results), page_data
end


text:
include("motion.jl")
include("estimate.jl")
include("correspondence.jl")
include("measurement.jl")
include("sparsification.jl")


text:
export VCFIterator, getVersion, getFilename, getSampleNames, getINFOProperties, getFORMATProperties, getContigs, getFilters, getReference, next1, eof, close, reset
import Base.eof, Base.close, Base.reset
import GZip

#Main class for VCFiterator
#type VCFIterator
mutable struct VCFIterator
    version::AbstractString
    filename::AbstractString
    samples::Array{AbstractString, 1}
    infoTypes::Dict{AbstractString, Any}
    infoFlags::Array{AbstractString, 1}
    formatTypes::Dict{AbstractString, Any}
    contigs::Dict{AbstractString, Int64}
    filters::Dict{AbstractString, AbstractString}
    reference::AbstractString
    file
end

#Stores INFO and FORMAT section properties
#type InfoField
mutable struct InfoField
    Number::AbstractString
    Type::AbstractString
    Description::AbstractString
end

#Standard get methods
getVersion(vc::VCFIterator) = vc.version
getFilename(vc::VCFIterator) = vc.filename
getSampleNames(vc::VCFIterator) = vc.samples
getINFOProperties(vc::VCFIterator) = vc.infoTypes
getFORMATProperties(vc::VCFIterator) = vc.formatTypes
getContigs(vc::VCFIterator) = vc.contigs
getFilters(vc::VCFIterator) = vc.filters
getReference(vc::VCFIterator) = vc.reference

#Gets next variant.
function next1(vc::VCFIterator)
    if !eof(vc.file)
        line = readline(vc.file)
        if length(line)<2
            return false
        end
        return generateVariant(vc, line)
    else
        return false
    end
end

#Checks if it is eof.
function eof(vc::VCFIterator)
    eof(vc.file)
end

#Closes the parser
function close(vc::VCFIterator)
    close(vc.file)
end

#resets the parser to the first variant of the VCF file.
function reset(vc::VCFIterator)
    close(vc.file)
    vc.file = open(vc.filename)
    f2=splitext(filename)
    if (uppercase(f2[2])==".GZ")
        close(vc.file)
        f=GZip.open(vc.filename)
    end

    while !eof(vc.file)
        line = chomp(readline(vc.file))
        if line[1:2] != "##"
            break
        end
    end
end

#Parses info/format line
function readInfoLine!(typesdict::Dict{AbstractString, Any}, line::AbstractString, flaglist)
    #parse the line
    entries = split(line[2:end-1], ',')
    Id = split(entries[1], '=')[2]
    Number = split(entries[2], '=')[2]
    Type = split(entries[3], '=')[2]
    Description = split(entries[4], '=')[2]
    temp = InfoField(Number, Type, Description[2:end-1])

    #add to dict
    typesdict[Id] = temp

    if Type == "Flag"
        push!(flaglist, Id)
    end
end

#Parses contig line
function readContigLine!(contigs::Dict{AbstractString, Int64}, line::AbstractString)
    entries = split(line[2:end-1], ',')
    Id = split(entries[1], '=')[2]
    #println(entries[2])
    Length = parse(Int32,split(entries[2], '=')[2])
    contigs[Id] = Length
end

#Parses filter line
function readFilterLine!(filters::Dict{AbstractString, AbstractString}, line::AbstractString)
    entries = split(line[2:end-1], ',')
    Id = split(entries[1], '=')[2]
    Description = split(entries[2], '=')[2]
    filters[Id] = Description[2:end-1]
end

#Constructor for VCFIterator
function VCFIterator(filename::AbstractString)
    f=open(filename)
    f2=splitext(filename)
    if (uppercase(f2[2])==".GZ")
        close(f)
        f=GZip.open(filename)
    end
    vc = VCFIterator("VCFv4.1", filename, [], Dict{AbstractString, Any}(), [], Dict{AbstractString, Any}(), Dict{AbstractString, Int64}(), Dict{AbstractString, AbstractString}(), "", f)
    parseHeader(vc)

    #return at the end
    vc
end

#Parses through header lines
function parseHeader(vc::VCFIterator)
    while !eof(vc.file)
        line = chomp(readline(vc.file))
        if line[1:2] != "##"
            samples = split(line, "\t")[10:end]
            vc.samples = samples
            break
        end

        #read the info line
        if length(line) >= 6 && line[1:6]=="##INFO"
            line = line[8:end]
            readInfoLine!(vc.infoTypes, line, vc.infoFlags)
        elseif length(line) >= 8 && line[1:8]=="##FORMAT"
            line = line[10:end]
            readInfoLine!(vc.formatTypes, line, AbstractString[])
        elseif length(line) >= 8 && line[1:8]=="##contig"
            line = line[10:end]
            readContigLine!(vc.contigs, line)
        elseif length(line) >= 8 && line[1:8]=="##FILTER"
            line = line[10:end]
            readFilterLine!(vc.filters, line)
        elseif length(line) >= 11 && line[1:11]=="##reference"
            line = line[13:end]
            vc.reference = line
        end

    end
end

#Helper function that convertss x to float except when x is a period.
function floatField(x)
    if x == "."
        return "."
    else
        return parse(Float64,x)
    end
end

#Helper function that convertss x to int except when x is a period.
function intField(x)
    if x == "."
        return "."
    else
        #println(x)
        return parse(Int32,x)
    end
end

#Constructor for variant given a line of vcf file
function generateVariant(vc::VCFIterator, line::AbstractString)
    #remove newline and split by tabs.
    line = split(chomp(line), "\t")

    #Assigining variables
    chr = convert(AbstractString, line[1])
    pos = intField(line[2])
    id_ = split(line[3])
    ref = convert(AbstractString, line[4])
    alt = split(line[5],",")
    if line[6] == "."
        qual = "."
    else
        qual = float(line[6])
    end
    filter = split(line[7],";")

    #creating INFO field
    information = Dict{AbstractString, Any}()
    tempINFO = split(line[8],";")

    #instantiate all flags to be false
    if tempINFO[1] != "."
        for i in 1:length(vc.infoFlags)
            information[vc.infoFlags[i]] = false
        end
        for i in 1:length(tempINFO)
            entry = tempINFO[i]
            id = split(entry,'=')[1]
            if ! haskey(vc.infoTypes, id)
                continue
            end
            #check if it is flag or not
            contentType = (vc.infoTypes[id]).Type

            if contentType != "Flag"
                contents = split(entry,'=')[2]
                #parse through contents and convert types
                contents = split(contents, ",")

                if contentType == "Integer"
                    contents = map(intField, contents)
                elseif contentType == "Float"
                    contents = map(floatField, contents)
                elseif contentType != "String" && contentType != "Char"
                    println("INVALID TYPE",":\t",entry)
                end

                information[id] = contents
            else
                #set variable to true if the flag exists
                information[id] = true
            end
        end
    end

    #creating FORMAT fields if they exist
    if length(line)>8
        format = Dict{AbstractString, Dict{AbstractString, Any}}()
        formatEntries = split(line[9],":")
        for i in 10:length(line)
            tempFORMAT = split(line[i],":")
            sample = vc.samples[i-9]
            tempFORMATdict = Dict{AbstractString, Any}()

            for j in 1:length(tempFORMAT)
                id = formatEntries[j]
                contents = split(tempFORMAT[j], ",")
                if ! haskey(vc.formatTypes, id)
                    continue
                end

                contentType = vc.formatTypes[id].Type

                if contentType == "Integer"
                    contents = map(intField, contents)
                elseif contentType == "Float"
                    contents = map(floatField, contents)
                elseif contentType != "String" && contentType != "Char"
	            println("INVALID TYPE",":\t",entry)
                end

                tempFORMATdict[id] = contents
            end
            format[sample] = tempFORMATdict
        end
    else
        format = Dict{AbstractString,Dict{AbstractString,Any}}()
    end

    #construct a variant
    Variant(chr, pos, id_, ref, alt, qual, filter, information, format)
end


text:
"""
Problem 1:

Given a path given by an input array of commands (direction + integer),
compute the final position relative to the start.

Problem 2:
Same, but changing slightly the rules.

"""

using DelimitedFiles

# Read input
function readInput() :: Array{Tuple{String,Int}}
  data = readdlm("./inputs/day2.txt", ' ', String)
  data = map(x -> (x[1],parse(Int,x[2])), data[i,:] for i in 1:size(data)[1] )
  return data
end

# Converts direction from string to cartesian index.
function toIndex(s::String)
  if s == "up"
    return CartesianIndex(0,-1)
  elseif s == "down"
    return CartesianIndex(0,1)
  elseif s == "forward"
    return CartesianIndex(1,0)
  else
    return CartesianIndex(-1,0)
  end
end

# Follows a path given by an array of steps (distance + direction). Return final position.
function followPath(path:: Array{Tuple{String,Int}}) :: CartesianIndex{2}
  pos = CartesianIndex(0,0)
  for step in path
    pos += step[2] * toIndex(step[1])
  end
  return pos
end

# Modified version for the second part of the problem.
function followPath2(path:: Array{Tuple{String,Int}}) :: CartesianIndex{2}
  aim = 0
  pos = CartesianIndex(0,0)
  for step in path
    dir = toIndex(step[1])
    aim += step[2] * dir[2]
    pos += CartesianIndex(step[2] * dir[1] ,aim * step[2] * dir[1])
  end
  return pos
end

path = readInput()

pos = followPath(path)
println("Problem 1:")
println(pos[1]*pos[2])

pos = followPath2(path)
println("Problem 2:")
println(pos[1]*pos[2])


text:
# RUN: julia -e "import Brutus; Brutus.lit(:emit_optimized)" --startup-file=no %s 2>&1 | FileCheck %s

function haspi(x::Union{Int64, Float64})
    if x isa Int64
        return x + 1
    end
end
emit(haspi, Union{Int64, Float64})

# CHECK: Core.MethodMatch(Tuple{typeof(Main.Main.haspi), Union{Float64, Int64}}, svec(), haspi(x::Union{Float64, Int64}) in Main.Main at /home/mccoy/Dev/brutus/test/Codegen/canonicalize/haspi.jl:3, true)after translating to MLIR in JLIR dialect:module  {
# CHECK-NEXT:   func nested @"Tuple{typeof(Main.haspi), Union{Float64, Int64}}"(%arg0: !jlir<"typeof(Main.haspi)">, %arg1: !jlir<"Union{Float64, Int64}">) -> !jlir<"Union{Nothing, Int64}"> attributes {llvm.emit_c_interface} {
# CHECK-NEXT:     "jlir.goto"()[^bb1] : () -> ()
# CHECK-NEXT:   ^bb1:  // pred: ^bb0
# CHECK-NEXT:     %0 = "jlir.constant"() {value = #jlir.isa} : () -> !jlir<"typeof(isa)">
# CHECK-NEXT:     %1 = "jlir.constant"() {value = #jlir.Int64} : () -> !jlir.DataType
# CHECK-NEXT:     %2 = "jlir.call"(%0, %arg1, %1) : (!jlir<"typeof(isa)">, !jlir<"Union{Float64, Int64}">, !jlir.DataType) -> !jlir.Bool
# CHECK-NEXT:     "jlir.gotoifnot"(%2)[^bb3, ^bb2] {operand_segment_sizes = dense<[1, 0, 0]> : vector<3xi32>} : (!jlir.Bool) -> ()
# CHECK-NEXT:   ^bb2:  // pred: ^bb1
# CHECK-NEXT:     %3 = "jlir.pi"(%arg1) : (!jlir<"Union{Float64, Int64}">) -> !jlir.Int64
# CHECK-NEXT:     %4 = "jlir.constant"() {value = #jlir<"#<intrinsic #2 add_int>">} : () -> !jlir<"typeof(Core.IntrinsicFunction)">
# CHECK-NEXT:     %5 = "jlir.constant"() {value = #jlir<"1">} : () -> !jlir.Int64
# CHECK-NEXT:     %6 = "jlir.call"(%4, %3, %5) : (!jlir<"typeof(Core.IntrinsicFunction)">, !jlir.Int64, !jlir.Int64) -> !jlir.Int64
# CHECK-NEXT:     %7 = "jlir.pi"(%6) : (!jlir.Int64) -> !jlir<"Union{Nothing, Int64}">
# CHECK-NEXT:     "jlir.return"(%7) : (!jlir<"Union{Nothing, Int64}">) -> ()
# CHECK-NEXT:   ^bb3:  // pred: ^bb1
# CHECK-NEXT:     %8 = "jlir.constant"() {value = #jlir.nothing} : () -> !jlir.Nothing
# CHECK-NEXT:     %9 = "jlir.pi"(%8) : (!jlir.Nothing) -> !jlir<"Union{Nothing, Int64}">
# CHECK-NEXT:     "jlir.return"(%9) : (!jlir<"Union{Nothing, Int64}">) -> ()
# CHECK-NEXT:   }
# CHECK-NEXT: }

# CHECK:   func nested @"Tuple{typeof(Main.haspi), Union{Float64, Int64}}"(%arg0: !jlir<"typeof(Main.haspi)">, %arg1: !jlir<"Union{Float64, Int64}">) -> !jlir<"Union{Nothing, Int64}"> attributes {llvm.emit_c_interface} {
# CHECK-NEXT:     "jlir.goto"()[^bb1] : () -> ()
# CHECK-NEXT:   ^bb1:  // pred: ^bb0
# CHECK-NEXT:     %0 = "jlir.constant"() {value = #jlir.Int64} : () -> !jlir.DataType
# CHECK-NEXT:     %1 = "jlir.isa"(%arg1, %0) : (!jlir<"Union{Float64, Int64}">, !jlir.DataType) -> !jlir.Bool
# CHECK-NEXT:     "jlir.gotoifnot"(%1)[^bb3, ^bb2] {operand_segment_sizes = dense<[1, 0, 0]> : vector<3xi32>} : (!jlir.Bool) -> ()
# CHECK-NEXT:   ^bb2:  // pred: ^bb1
# CHECK-NEXT:     %2 = "jlir.pi"(%arg1) : (!jlir<"Union{Float64, Int64}">) -> !jlir.Int64
# CHECK-NEXT:     %3 = "jlir.constant"() {value = #jlir<"1">} : () -> !jlir.Int64
# CHECK-NEXT:     %4 = "jlir.add_int"(%2, %3) : (!jlir.Int64, !jlir.Int64) -> !jlir.Int64
# CHECK-NEXT:     %5 = "jlir.pi"(%4) : (!jlir.Int64) -> !jlir<"Union{Nothing, Int64}">
# CHECK-NEXT:     "jlir.return"(%5) : (!jlir<"Union{Nothing, Int64}">) -> ()
# CHECK-NEXT:   ^bb3:  // pred: ^bb1
# CHECK-NEXT:     %6 = "jlir.constant"() {value = #jlir.nothing} : () -> !jlir.Nothing
# CHECK-NEXT:     %7 = "jlir.pi"(%6) : (!jlir.Nothing) -> !jlir<"Union{Nothing, Int64}">
# CHECK-NEXT:     "jlir.return"(%7) : (!jlir<"Union{Nothing, Int64}">) -> ()
# CHECK-NEXT:   }
# CHECK-NEXT: }

max_stars_repo_name: JuliaLabs/brutus

max_issues_repo_name: JuliaLabs/brutus

max_forks_repo_name: JuliaLabs/brutus


text:
# Autogenerated wrapper script for Binutils_jll for x86_64-linux-gnu-target_libc+musl-target_os+linux-target_arch+x86_64
export bindir

using Zlib_jll
JLLWrappers.@generate_wrapper_header("Binutils")
JLLWrappers.@declare_file_product(bindir)
function __init__()
    JLLWrappers.@generate_init_header(Zlib_jll)
    JLLWrappers.@init_file_product(
        bindir,
        "bin",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
#= Bool functions =#
isqnode(ex) = false
isqnode(ex::QuoteNode) = true
isqnodewithdim(ex) = false
isqnodewithdim(ex::Expr) = ex.head == :call && length(ex.args) == 3 &&
    ex.args[1] == :| && isqnode(ex.args[2])
isintindexproduct(ex) = false
isintindexproduct(ex::Expr) = ex.head == :call && length(ex.args) == 3 &&
    ex.args[1] == :* && ex.args[2:3] .|> typeof |> Set == (Int,QuoteNode) |> Set

isrhspairedindex(ex) = false
isrhspairedindex(ex::Expr) = ex.head == :tuple &&
    all(ex.args .|> x -> isqnodewithdim(x) || isqnode(x))
isrhsindex(ex) = isqnode(ex) || isrhspairedindex(ex) || isintindexproduct(ex)
isrhstensor(ex) = false
isrhstensor(ex::Expr) = ex.head == :ref && all(ex.args[2:end] .|> isrhsindex)
isrhssimpletensor(ex) = false
isrhssimpletensor(ex::Expr) = ex.head == :ref && all(ex.args[2:end] .|> isqnode)

islhspairedindex(ex) = false
islhspairedindex(ex::Expr) = ex.head == :tuple && all(ex.args .|> isqnode)
islhsindex(ex) = isqnode(ex) || islhspairedindex(ex)
islhstensor(ex) = false
islhstensor(ex::Expr) = ex.head == :vect && all(ex.args .|> islhsindex)
islhssimpletensor(ex) = false
islhssimpletensor(ex::Expr) =  ex.head == :vect && all(ex.args .|> isqnode)

istensorproduct(ex) = false
istensorproduct(ex::Expr) = ex.head == :call && ex.args[1] == :* &&
    all(ex.args[2:end] .|> isrhstensor)
issimpletensorproduct(ex) = false
issimpletensorproduct(ex::Expr) = ex.head == :call && ex.args[1] == :* &&
    all(ex.args[2:end] .|> isrhssimpletensor)

#= elementary functions to parse the Expr =#
function indint(ex::Expr)
    if typeof(ex.args[2]) == QuoteNode
        ex.args[2],ex.args[3]
    else
        ex.args[3],ex.args[2]
    end
end

function dupinds(indslis)
    nodup = QuoteNode[]
    dup = QuoteNode[]
    for inds in indslis
        for i in inds
            if i in dup
                error(i," appears theree times")
            elseif i in nodup
                filter!(x->x!=i,nodup)
                push!(dup,i)
            else
                push!(nodup,i)
            end
        end
    end
    nodup,dup
end

# main steps =#
function headlhsrhs(ex::Expr)
    if ex.head == :call && ex.args[1] == :(<=) && length(ex.args) == 3 &&
        islhstensor(ex.args[2]) && istensorproduct(ex.args[3])
        ex.args[1],ex.args[2],ex.args[3]
    elseif ex.head == :call && ex.args[1] == :(=>) && length(ex.args) == 3 &&
        islhstensor(ex.args[3]) && istensorproduct(ex.args[2])
        ex.args[1],ex.args[3],ex.args[2]
    else
        error("cannot parse the Expr")
    end
end

function bonddimdict(ex::Expr)
    resdict = Dict{QuoteNode,Union{Int,Symbol,Expr,Nothing}}()
    pairedinddict = Dict{Array{QuoteNode,1},Expr}()
    for t in ex.args[2:end]
        posreshaped = posoriginal = 1
        for ind in t.args[2:end]
            if isqnode(ind)
                resdict[ind] = :(size($(t.args[1]),$posoriginal))
                posreshaped += 1
                posoriginal += 1
            elseif isrhspairedindex(ind)
                pairedind = QuoteNode[]
                for i in ind.args
                    if typeof(i) == Expr
                        resdict[i.args[2]] = i.args[3]
                        push!(pairedind,i.args[2])
                    elseif typeof(i) == QuoteNode
                        if !haskey(resdict,i)
                            resdict[i] = nothing
                        end
                        push!(pairedind,i)
                    end
                end
                pairedinddict[pairedind] = :(size($(t.args[1]),$posoriginal))
                posreshaped += length(ind.args)
                posoriginal += 1
            elseif isintindexproduct(ind)
                indname,posshift = indint(ind)
                resdict[indname] =
                    :(size($(t.args[1]))[$posoriginal:$posoriginal+$posshift-1]|>prod)
                posreshaped += 1
                posoriginal += posshift
            end
        end
    end
    for i in filter(x->resdict[x]==nothing,keys(resdict))
        for j in filter(x->(i in x),keys(pairedinddict))
            if count(x->resdict[x]==nothing,j) == 1
                denom =
                Expr(:call,:*,map(x->resdict[x],filter(x->resdict[x]!=nothing,j))...)
                resdict[i] = Expr(:call,:div,pairedinddict[j],denom)
                break
            end
            error("cannot predict the dimension of ",i)
        end
    end
    resdict
end

function tosimpletensor!(ex,bddict::Dict{QuoteNode,<:Any})
    ex.args[2:end] .= map(x->_tosimpletensor(x,bddict),ex.args[2:end])
end
function _tosimpletensor(ex,bddict::Dict{QuoteNode,<:Any})
    if isrhssimpletensor(ex)
        return ex
    end
    tname = Expr(:call,:reshape,ex.args[1])
    exout = Expr(:ref,tname)
    for i in ex.args[2:end]
        if isqnode(i)
            push!(tname.args,bddict[i])
            push!(exout.args,i)
        elseif isintindexproduct(i)
            push!(tname.args,bddict[indint(i)[1]])
            push!(exout.args,indint(i)[1])
        elseif isrhspairedindex(i)
            for j in i.args
                if typeof(j) == QuoteNode
                    push!(tname.args,bddict[j])
                    push!(exout.args,j)
                else
                    push!(tname.args,bddict[j.args[2]])
                    push!(exout.args,j.args[2])
                end
            end
        end
    end
    exout
end

function taketrace!(ex::Expr)
    ex.args[2:end] .= map(_taketrace,ex.args[2:end])
end
function _taketrace(ex::Expr)
    tname,tind = ex.args[1],ex.args[2:end]
    if length(dupinds([tind])[2]) == 0
        ex
    else
        newtind = dupinds([tind])[1]
        Expr(:ref,:($tracefunc($tname,$tind,$newtind)),newtind...)
    end
end

function makepairwised(ex::Expr,ord::Union{Array{QuoteNode,1},Array{Nothing,1}})
    if length(ex.args) == 3
        return ex
    end
    nodup,dup = dupinds([i.args[2:end] for i in ex.args[2:end]])
    contracttensor = Int[]
    contractind = nothing
    for i in ord
        if i in dup
            contractind = i
            break
        end
    end
    if contractind != nothing
        append!(contracttensor,
            filter(x->(contractind in ex.args[x].args[2:end]),2:length(ex.args)))
    else
        push!(contracttensor,2,3)
    end
    newind,dst = dupinds([i.args[2:end] for i in ex.args[contracttensor]])
    exout = Expr(:call,:*,
        map(i->ex.args[i],filter(x->!(x in contracttensor),2:length(ex.args)))...)
    push!(exout.args,Expr(:ref,Expr(:call,:*,ex.args[contracttensor]...),newind...))
    makepairwised(exout,ord)
end

function parsetensorproduct(ex::Expr)
    if istensorproduct(ex.args[1])
        l = parsetensorproduct(ex.args[1].args[2])
        r = parsetensorproduct(ex.args[1].args[3])
        :($contractfunc($l,$(ex.args[1].args[2].args[2:end]|>Tuple),
            $r,$(ex.args[1].args[3].args[2:end]|>Tuple),$(ex.args[2:end]|>Tuple)))
    else
        :($(ex.args[1]))
    end
end

function toindreshape(indslis,bddict::Dict{QuoteNode,<:Any})
    resultindex = []
    resultshape = []
    for i in indslis
        if isqnode(i)
            push!(resultindex,i)
            push!(resultshape,bddict[i])
        else
            append!(resultindex,i.args)
            push!(resultshape,Expr(:call,:*,map(x->bddict[x],i.args)...))
        end
    end
    resultindex,resultshape
end

#= definition of macro =#
function tensorproductmain(ex::Expr,ord::Union{Array{QuoteNode,1},Array{Nothing,1}})
    head,lhs,rhs = headlhsrhs(ex)
    bdims = bonddimdict(rhs)
    tosimpletensor!(rhs,bdims)
    taketrace!(rhs)
    rhs = makepairwised(rhs,ord)
    lhsind,lhsreshape = toindreshape(lhs.args,bdims)
    rhs = Expr(:ref,rhs,lhsind...)
    rhs = parsetensorproduct(rhs)
    if length(lhs.args) != length(lhsind)
        rhs = Expr(:call,:reshape,rhs,lhsreshape...)
    end
    rhs
end

"""
    @tensorfunc(optex, block)
    @tensorfunc(block)

    TODO: write docs
"""
macro tensorfunc(ord::Expr,ex::Expr)
    esc(tensorproductmain(ex,order(ex,ord)))
end

macro tensorfunc(ex::Expr)
    esc(tensorproductmain(ex,order(ex)))
end


text:
### A Pluto.jl notebook ###
# v0.17.3

using Markdown
using InteractiveUtils

# ╔═╡ 9be09c2c-84fc-11ec-34d5-032693be6678
md"""

## Guides

content for this page is coming soon ...
"""

# ╔═╡ 2c809ee8-9d3e-441e-b351-eb3f8f11e875
md"""
* [Working with Genie apps](/docs/guides/Working-With-Genie-Apps.html)
* [Using Genie Plugins](/docs/guides/Genie-Plugins.html)
* [Developing an API backend](/docs/guides/Simple-API-backend.html)
* [Using Genie in an interactive environment](/docs/guides/Interactive-environment.html)
"""

# ╔═╡ Cell order:
# ╠═9be09c2c-84fc-11ec-34d5-032693be6678
# ╟─2c809ee8-9d3e-441e-b351-eb3f8f11e875


text:
using Test
using TestSetExtensions
using LinearAlgebra
using Qaintessent
using Qaintmodels

function isunitary(cb::Vector{CircuitGate})
    sparse_matrix(cb) * sparse_matrix(Base.adjoint(cb)) ≈ I
end

@testset ExtendedTestSet "quantum Fourier transform" begin
    for N in 1:4
        c = qft_circuit(N)

        Fmat = [exp(2*π*im*j*k/2^N)/sqrt(2^N) for k in 0:(2^N-1), j in 0:(2^N-1)]

        @test Qaintessent.sparse_matrix(c) ≈ Fmat
        @test isunitary(c)
    end
end

@testset ExtendedTestSet "toffoli circuit test" begin
    N = 3
    cgc = toffoli_circuit(1, (3, 2))

    # reference
    toffoli = circuit_gate(1, X, 3, 2)

    @test Qaintessent.sparse_matrix(toffoli) ≈ Qaintessent.sparse_matrix(cgc)
end

@testset ExtendedTestSet "vbe adder test" begin
    for N in 1:4
        M = 3N + 1
        adder = vbe_adder_circuit(N)
        ψ = fill(0.0 + 0.0im, 2^M)
        a = rand(0:2^(N-1))
        b = rand(0:2^(N-1))

        index = b << N + a
        ψ[index+1] = 1.0
        ψ = apply(ψ, adder)
        answer = ((findall(x -> x == 1, ψ)[1] - 1) % (2^2N)) >> N
        @test answer == (a + b) % (2^N)
    end
end

@testset ExtendedTestSet "qcla out of place adder test" begin
    for N in 1:4
        anc = N - count_ones(N) - floor(Int, log(N))
        M = 3N + anc + 1

        cgc = qcla_out_adder_circuit(N)
        ψ = fill(0.0 + 0.0 * im, 2^M)
        a = rand(0:2^(N-1))
        b = rand(0:2^(N-1))

        index = b << N + a

        ψ[index+1] = 1.0

        ψ = apply(ψ, cgc)
        answer = (findall(x -> x == 1, ψ)[1] - 1) >> 2N
        @test answer == a + b
    end
end

@testset ExtendedTestSet "qcla in place adder test" begin
    for N in 1:4
        n = 1
        anc = 0
        while 2^n < N
            anc += N ÷ 2^n -1
            n += 1
        end
        M = 3N + anc
        cgc = qcla_inplace_adder_circuit(N)
        a = rand(0:2^(N-1))
        b = rand(0:2^(N-1))

        index = b << N + a
        ψ = fill(0.0 + 0.0 * im, 2^M)

        ψ[index+1] = 1.0

        ψ = apply(ψ, cgc)
        answer = ((findall(x -> x == 1, ψ)[1] - 1) >> N)
        @test answer == a + b
    end
end


text:
using ITensors
using ITensorVisualizationBase
using Test

@testset "Basic tests without any backend" begin
  N = 10
  s(n) = Index([QN("Sz", 0) => 1, QN("Sz", 1) => 1]; tags="S=1/2,Site,n=$n")
  l(n) = Index([QN("Sz", 0) => 10, QN("Sz", 1) => 10]; tags="Link,l=$n")
  h(n) = Index([QN("Sz", 0) => 5, QN("Sz", 1) => 5]; tags="ham,Link,l=$n")
  s⃗ = [s(n) for n in 1:N]
  l⃗ = [l(n) for n in 1:(N - 1)]
  h⃗ = [h(n) for n in 1:(N - 1)]

  # Add some more indices between two of the tensors
  x = Index([QN("Sz", 0) => 2]; tags="X")
  y = Index([QN("Sz", 0) => 2]; tags="Y")

  n = 2
  ψn1n2 = randomITensor(l⃗[n - 1], s⃗[n], s⃗[n + 1], l⃗[n + 1], dag(x), dag(y))
  hn1 = randomITensor(dag(h⃗[n - 1]), s⃗[n]', dag(s⃗[n]), h⃗[n], x, y)
  hn2 = randomITensor(dag(h⃗[n]), s⃗[n + 1]', dag(s⃗[n + 1]), h⃗[n + 1])
  ELn0 = randomITensor(l⃗[n - 1]', h⃗[n - 1], dag(l⃗[n - 1]))
  ERn2 = randomITensor(l⃗[n + 1]', dag(h⃗[n + 1]), dag(l⃗[n + 1]))

  tn = [ELn0, ψn1n2, hn1, hn2, ERn2]

  R = @visualize ELn0 * ψn1n2 * hn1 * hn2 * ERn2
  R1 = @visualize ELn0 * ψn1n2 * hn1
  R2 = @visualize R1 * hn2 * ERn2 vertex_labels = ["T1", "T2", "T3"]
  tn2 = @visualize tn
  T = @visualize ELn0

  @test R ≈ ELn0 * ψn1n2 * hn1 * hn2 * ERn2
  @test R1 ≈ ELn0 * ψn1n2 * hn1
  @test R2 ≈ ELn0 * ψn1n2 * hn1 * hn2 * ERn2
  @test all(tn .== tn2)
  @test T == ELn0

  R = @visualize figR ELn0 * ψn1n2 * hn1 * hn2 * ERn2
  R1 = @visualize figR1 ELn0 * ψn1n2 * hn1
  R2 = @visualize figR2 R1 * hn2 * ERn2 vertex_labels = ["T1", "T2", "T3"]
  T = @visualize figT T

  fig_tn = @visualize_noeval tn

  @test isnothing(figR)
  @test isnothing(figR1)
  @test isnothing(figR2)
  @test isnothing(fig_tn)
  @test isnothing(figT)
end


text:
export Cora

"""
    Cora

The Cora citation network dataset from Ref. [1].
Nodes represent documents and edges represent citation links.
Each node has a predefined feature with 1433 dimensions.
The dataset is designed for the node classification task.
The task is to predict the category of certain paper.
The dataset is retrieved from Ref. [2].

## Statistics

- Nodes: 2708
- Edges: 10556
- Number of Classes: 7
- Label split:
    - Train:  140
    - Val:    500
    - Test:  1000

The split is the one used in the original paper [1] and
doesn't consider all nodes.

## Interface

- [`Cora.dataset`](@ref)

## References

[1]: [Deep Gaussian Embedding of Graphs: Unsupervised Inductive Learning via Ranking](https://arxiv.org/abs/1707.03815)
[2]: [Planetoid](https://github.com/kimiyoung/planetoid
"""
module Cora

using DataDeps
using ..MLDatasets: datafile, read_planetoid_data
using DelimitedFiles: readdlm

using PyCall

const DEPNAME = "Cora"
# LINK = "https://github.com/shchur/gnn-benchmark/raw/master/data/npz"
# LINK = "https://github.com/abojchevski/graph2gauss/raw/master/data/"
const LINK = "https://github.com/kimiyoung/planetoid/raw/master/data"
const DOCS = "https://github.com/kimiyoung/planetoid"
const DATA = "ind.cora." .* ["x", "y", "tx", "allx", "ty", "ally", "graph", "test.index"]

function __init__()
    register(DataDep(
        DEPNAME,
        """
        Dataset: The $DEPNAME dataset.
        Website: $DOCS
        """,
        map(x -> "$LINK/$x", DATA),
        "81de017067dc045ebdb8ffd5c0e69a209973ffdb1fe2d5b434e94d3614f3f5c7",  # if checksum omitted, will be generated by DataDeps
        # post_fetch_method = unpack
    ))
end

"""
    dataset(; dir=nothing, reverse_edges=true)

Retrieve the Cora dataset. The output is a named tuple with fields
```juliarepl
julia> keys(Cora.dataset())
(:node_features, :node_labels, :adjacency_list, :train_indices, :val_indices, :test_indices, :num_classes, :num_nodes, :num_edges, :directed)
```
In particular, `adjacency_list` is a vector of vector,
where `adjacency_list[i]` will contain the neighbors of node `i`
through outgoing edges.

If `reverse_edges=true`, the graph will contain
the reverse of each edge and the graph will be undirected.

See also [`Cora`](@ref).

## Usage Examples

```juliarepl
using MLDatasets: Cora

data = Cora.dataset()
train_labels = data.node_labels[data.train_indices]
```
"""
dataset(; dir=nothing, reverse_edges=true) =
    read_planetoid_data(DEPNAME, dir=dir, reverse_edges=reverse_edges)

end #module


text:
"""
# plot_models

Multiple regression coefficient plot for multiple models.

$(SIGNATURES)

### Required arguments
```julia
* `models`                             : Vector of `SampleModel`s to compare
* `pars`                               : Vector of parameters to include in comparison
```
### Optional arguments
```julia
* `fig=""`                             : File to store plot
* `title=""`                           : Title for plot
```

### Return values
```julia
* `(s, f)`                             : (particles, plot)
```

"""
function plot_models(models::Vector{SampleModel}, pars::Vector{Symbol};
  fig="", title="")

  mnames = [models[i].name for i in 1:length(models)]
  for i in 1:length(mnames)
    if length(mnames[i]) > 9
      mnames[i] = mnames[i][1:9]
    end
  end

  s = Vector{NamedTuple}(undef, length(models))
  for (mindx, mdl) in enumerate(models)
    df = read_samples(mdl, :dataframe)
    m, l, u = estimparam(df)
    d = Dict{Symbol, NamedTuple}()
    for (indx, par) in enumerate(names(df))
        d[Symbol(par)] = (mean=m[indx], lower=l[indx], upper=u[indx])
    end
    s[mindx] =   (; d...)
  end

  plot_models(s, mnames, pars; fig, title)
end

export
  plot_models


text:
function conic_form!(x::Variable, unique_conic_forms::UniqueConicForms=UniqueConicForms())
    if !has_conic_form(unique_conic_forms, x)
        add_to_id_to_variables!(unique_conic_forms, x)
        if vexity(x) == ConstVexity()
            # do exactly what we would for a constant
            objective = ConicObj()
            objective[objectid(:constant)] = (vec([real(x.value);]),vec([imag(x.value);]))
            cache_conic_form!(unique_conic_forms, x, objective)
        else
            objective = ConicObj()
            vec_size = length(x)

            objective[x.id_hash] = (real_conic_form(x), imag_conic_form(x))
            objective[objectid(:constant)] = (spzeros(vec_size, 1), spzeros(vec_size, 1))
            # placeholder values in unique constraints prevent infinite recursion depth
            cache_conic_form!(unique_conic_forms, x, objective)
            if !(x.sign == NoSign() || x.sign == ComplexSign())
                conic_form!(x.sign, x, unique_conic_forms)
            end
            for set in x.sets
                conic_form!(set, x, unique_conic_forms)
            end
        end
    end
    return get_conic_form(unique_conic_forms, x)
end


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

module GitTools

using ..Pkg
import LibGit2
using Printf

Base.@kwdef mutable struct MiniProgressBar
    max::Float64 = 1.0
    header::String = ""
    color::Symbol = :white
    width::Int = 40
    current::Float64 = 0.0
    prev::Float64 = 0.0
    has_shown::Bool = false
    time_shown::Float64 = 0.0
end

const NONINTERACTIVE_TIME_GRANULARITY = Ref(2.0)
const PROGRESS_BAR_PERCENTAGE_GRANULARITY = Ref(0.1)

function showprogress(io::IO, p::MiniProgressBar)
    perc = p.current / p.max * 100
    prev_perc = p.prev / p.max * 100
    # Bail early if we are not updating the progress bar,
    # Saves printing to the terminal
    if p.has_shown && !((perc - prev_perc) > PROGRESS_BAR_PERCENTAGE_GRANULARITY[])
        return
    end
    if !isinteractive()
        t = time()
        if p.has_shown && (t - p.time_shown) < NONINTERACTIVE_TIME_GRANULARITY[]
            return
        end
        p.time_shown = t
    end
    p.prev = p.current
    p.has_shown = true
    n_filled = ceil(Int, p.width * perc / 100)
    n_left = p.width - n_filled
    print(io, "    ")
    printstyled(io, p.header, color=p.color, bold=true)
    print(io, " [")
    print(io, "="^n_filled, ">")
    print(io, " "^n_left, "]  ", )
    @printf io "%2.1f %%" perc
    print(io, "\r")
end

function transfer_progress(progress::Ptr{LibGit2.TransferProgress}, p::Any)
    progress = unsafe_load(progress)
    @assert haskey(p, :transfer_progress)
    bar = p[:transfer_progress]
    @assert typeof(bar) == MiniProgressBar
    if progress.total_deltas != 0
        bar.header = "Resolving Deltas:"
        bar.max = progress.total_deltas
        bar.current = progress.indexed_deltas
    else
        bar.max = progress.total_objects
        bar.current = progress.received_objects
    end
    showprogress(stdout, bar)
    return Cint(0)
end

const GITHUB_REGEX =
    r"^(?:git@|git://|https://(?:[\w\.\+\-]+@)?)github.com[:/](([^/].+)/(.+?))(?:\.git)?$"i
const GIT_PROTOCOL = Ref{Union{String, Nothing}}(nothing)

setprotocol!(proto::Union{Nothing, AbstractString}=nothing) = GIT_PROTOCOL[] = proto

# TODO: extend this to more urls
function normalize_url(url::AbstractString)
    m = match(GITHUB_REGEX, url)
    if m === nothing || GIT_PROTOCOL[] === nothing
        url
    elseif GIT_PROTOCOL[] == "ssh"
        "ssh://git@github.com/$(m.captures[1]).git"
    else
        "$(GIT_PROTOCOL[])://github.com/$(m.captures[1]).git"
    end
end

function clone(url, source_path; header=nothing, kwargs...)
    @assert !isdir(source_path) || isempty(readdir(source_path))
    Pkg.Types.printpkgstyle(stdout, :Cloning, header == nothing ? "git-repo `$url`" : header)
    transfer_payload = MiniProgressBar(header = "Fetching:", color = Base.info_color())
    callbacks = LibGit2.Callbacks(
        :transfer_progress => (
            @cfunction(transfer_progress, Cint, (Ptr{LibGit2.TransferProgress}, Any)),
            transfer_payload,
        )
    )
    url = normalize_url(url)
    print(stdout, "\e[?25l") # disable cursor
    try
        return LibGit2.clone(url, source_path; callbacks=callbacks, kwargs...)
    catch err
        rm(source_path; force=true, recursive=true)
        err isa LibGit2.GitError || rethrow()
        if (err.class == LibGit2.Error.Net && err.code == LibGit2.Error.EINVALIDSPEC) ||
           (err.class == LibGit2.Error.Repository && err.code == LibGit2.Error.ENOTFOUND)
            Pkg.Types.pkgerror("Git repository not found at '$(url)'")
        else
            Pkg.Types.pkgerror("failed to clone from $(url), error: $err")
        end
    finally
        print(stdout, "\033[2K") # clear line
        print(stdout, "\e[?25h") # put back cursor
    end
end

function fetch(repo::LibGit2.GitRepo, remoteurl=nothing; header=nothing, kwargs...)
    if remoteurl === nothing
        remoteurl = LibGit2.with(LibGit2.get(LibGit2.GitRemote, repo, "origin")) do remote
            LibGit2.url(remote)
        end
    end
    remoteurl = normalize_url(remoteurl)
    Pkg.Types.printpkgstyle(stdout, :Updating, header == nothing ? "git-repo `$remoteurl`" : header)
    transfer_payload = MiniProgressBar(header = "Fetching:", color = Base.info_color())
    callbacks = LibGit2.Callbacks(
        :transfer_progress => (
            @cfunction(transfer_progress, Cint, (Ptr{LibGit2.TransferProgress}, Any)),
            transfer_payload,
        )
    )
    print(stdout, "\e[?25l") # disable cursor
    try
        return LibGit2.fetch(repo; remoteurl=remoteurl, callbacks=callbacks, kwargs...)
    catch err
        err isa LibGit2.GitError || rethrow()
        if (err.class == LibGit2.Error.Repository && err.code == LibGit2.Error.ERROR)
            Pkg.Types.pkgerror("Git repository not found at '$(remoteurl)'")
        else
            Pkg.Types.pkgerror("failed to fetch from $(remoteurl), error: $err")
        end
    finally
        print(stdout, "\033[2K") # clear line
        print(stdout, "\e[?25h") # put back cursor
    end
end

end # module


text:
struct BIDSDirectory <: AbstractPath
end

struct ParticipantPath <: AbstractPath
end

function participantlabel end

function sessionlabel end

function tasklabel end

anatdata
funcdata
dwidata
fmapdata
behdata
megdata
eegdata
ieegdata
channelsdata


text:
using Statistics, Plots

#------------------------------------------------------------------------------
struct OPT
	xMin::Float64
	xMax::Float64
	xIni::Vector{Float64}
	VolMax::Float64
	Tol::Float64
	MaxIter::Int
	OCMove::Float64
	OCEta::Float64
end

#------------------------------------------------------------------------------
function InitOPT(xIni::Vector{Float64}, VolMax::Float64, Tol::Float64, MaxIter::Int,
	xMin::Float64=1e-4*mean(xIni), xMax::Float64=1e+4*mean(xIni),
	OCMove::Float64=xMax, OCEta::Float64=0.5)::OPT

	OPT(xMin, xMax, xIni, VolMax, Tol, MaxIter, OCMove, OCEta)
end

#------------------------------------------------------------------------------
function TrussTop(FEA::FEM, Opt::OPT)::Tuple{Vector{Vector{Float64}}, Vector{Float64}, FEM}

	Iter, Tol, Change = 0, Opt.Tol, Opt.Tol + 1
	x = Opt.xIni; xTemp = similar(x)

	∂f∂x = zeros(FEA.NElem)
	∂g∂x = zeros(FEA.NElem)
	U = zeros(2FEA.NNode)

	xHist = Vector{Float64}[]; sizehint!(xHist, Opt.MaxIter + 1)
	fHist = Float64[]; sizehint!(fHist, Opt.MaxIter + 1)

	while (Iter < Opt.MaxIter) && (Change > Tol)
		Iter += 1

		# Compute cost functionals and analysis sensitivities
		f = ObjectiveFcn!(∂f∂x, U, x, FEA)
		g = ConstraintFcn!(∂g∂x, x, FEA, Opt)

		# Store design vars. and obj. function history
		push!(xHist, x)
		push!(fHist, f)

		# Update design variable and analysis parameters
		Change = UpdateScheme!(x, xTemp, ∂f∂x, g, ∂g∂x, Opt)
	end

	# Store last design var. and obj. function value
	f = ObjectiveFcn!(∂f∂x, U, x, FEA)
	push!(xHist, x); push!(fHist, f)

	# Return arguments
	xHist, fHist, FEA
end

#------------------------------------------------------------------------------
function ObjectiveFcn!(∂f∂x::Vector{Float64}, U::Vector{Float64}, x::Vector{Float64}, FEA::FEM)::Float64

	FEAnalysis!(U, x, FEA)

	@fastmath @inbounds @. FEA.Sparse.k = -U[FEA.Sparse.i] * FEA.Sparse.k0 * U[FEA.Sparse.j]
	@fastmath cumsum!(FEA.Sparse.k, FEA.Sparse.k)
	∂f∂x[1] = FEA.Sparse.k[FEA.CumSumElemNDof2[1]]
	@fastmath ∂f∂x[2:end] = diff(FEA.Sparse.k[FEA.CumSumElemNDof2])

	@fastmath f = FEA.Force ⋅ U
end

# function ObjectiveFcn!(∂f∂x::Vector{Float64}, U::Vector{Float64}, x::Vector{Float64}, FEA::FEM)::Float64

# 	FEAnalysis!(U, x, FEA)

# 	curr = 0.; prev = 0.; idx = 0
# 	for i in eachindex(FEA.Sparse.i)
# 		curr -= U[FEA.Sparse.i[i]] * FEA.Sparse.k0[i] * U[FEA.Sparse.j[i]]
# 		if (i % FEA.ElemNDof[idx+1]^2) == 0
# 			∂f∂x[idx+=1] = curr - prev
# 			prev = curr
# 		end
# 	end

# 	f = FEA.Force ⋅ U
# end

#------------------------------------------------------------------------------
function ConstraintFcn!(∂g∂x::Vector{Float64}, x::Vector{Float64}, FEA::FEM, Opt::OPT)

	@fastmath ∂g∂x .= FEA.ElemLen
	@fastmath g = FEA.ElemLen ⋅ x - Opt.VolMax
end

#------------------------------------------------------------------------------
function UpdateScheme!(x::Vector{Float64}, xTemp::Vector{Float64}, ∂f∂x::Vector{Float64},
	g::Float64, ∂g∂x::Vector{Float64}, Opt::OPT)::Float64

	xMin, xMax, move, η = Opt.xMin, Opt.xMax, Opt.OCMove, Opt.OCEta

	@fastmath @. ∂f∂x = min(∂f∂x, 0.)
	@fastmath @. xTemp = -∂f∂x / ∂g∂x
	@fastmath λ₁, λ₂ = 0., 1.2 * maximum(xTemp)

	@fastmath @. xTemp = x
	while (λ₂ - λ₁) > (1e-10 * (1 + λ₂))
	  λₘ = (λ₁ + λ₂) / 2.

	  @fastmath @. x = xMin + (xTemp - xMin) * ((-∂f∂x / ∂g∂x / λₘ) ^ η)
	  @fastmath @. x = min(x, xTemp + move, xMax)
	  @fastmath @. x = max(x, xTemp - move, xMin)

	  @fastmath if (g + ∂g∂x ⋅ (x - xTemp))>0 λ₁ = λₘ else λ₂ = λₘ end
	end

	@fastmath @. xTemp = abs(x - xTemp) / (1. + xTemp)
	@fastmath Δ = maximum(xTemp)
end

# function UpdateScheme!(x::Vector{Float64}, xTemp::Vector{Float64}, ∂f∂x::Vector{Float64},
# 	g::Float64, ∂g∂x::Vector{Float64}, Opt::OPT)::Float64

# 	xMin, xMax, move, η = Opt.xMin, Opt.xMax, Opt.OCMove, Opt.OCEta

# 	λ₁, λ₂, = 0., -Inf
# 	for i in eachindex(x)
# 		∂f∂x[i] = min(∂f∂x[i], 0.)
# 		λ₂ = max(λ₂, -∂f∂x[i] / ∂g∂x[i])
# 		xTemp[i] = x[i]
# 	end
# 	λ₂ *= 1.2

# 	while (λ₂ - λ₁) > (1e-10 * (1 + λ₂))
# 		c, λₘ = g, (λ₁ + λ₂) / 2.

# 		for i in eachindex(x)
# 			x[i] = xMin + (xTemp[i] - xMin) * ((-∂f∂x[i] / ∂g∂x[i] / λₘ) ^ η)
# 			x[i] = min(x[i], xTemp[i] + move, xMax)
# 			x[i] = max(x[i], xTemp[i] - move, xMin)

# 			c += ∂g∂x[i] * (x[i] - xTemp[i])
# 		end
# 		if c>0 λ₁ = λₘ else λ₂ = λₘ end
# 	end

# 	Δ = -Inf
# 	for i in eachindex(x)
# 		Δ = max(Δ, abs(x[i] - xTemp[i]) / (1 + xTemp[i]))
# 	end
# 	return Δ
# end

#------------------------------------------------------------------------------
function TrussTopView(xH::Vector{Vector{Float64}}, fH::Vector{Float64}, fem::FEM, Filter::Float64)

	# Apply visualization filter
	x = xH[end]; MaxArea = maximum(x);
	inactive = findall(x .< Filter * MaxArea)
	active = setdiff(1:length(x), inactive);
	xEff = copy(x); xEff[inactive] .= 0

	# Summary results
	@printf("Number of bars: %d of %d\n", length(active), length(x));
	@printf("Objective function value: %g\n", fH[end]);
	@printf("Number of iterations: %d\n", length(fH) - 1);

	# Calc stress
	U = zeros(2fem.NNode); FEAnalysis!(U, x, fem)
	σ = zeros(fem.NElem); Stress!(σ, U, fem)

	# Final topology
	FinalTopology = plot()
	# println(active)
	for e in active
		X = vec(fem.Node[fem.Elem[e,:],1])
		Y = vec(fem.Node[fem.Elem[e,:],2])
		clr = "blue"; if σ[e]<0; clr="red"; end
		global FinalTopology = plot!(X,Y,color=clr,linewidth=xEff[e]/MaxArea*5,legend=false, title="Final topology",
		aspect_ratio=:equal, grid=false, axis=false, ticks=false)
	end
	ActiveNodes = unique(fem.Elem[active,:])
	X = vec(fem.Node[ActiveNodes,1])
	Y = vec(fem.Node[ActiveNodes,2])
	FinalTopology = scatter!(X,Y,color="white",markersize=5, legend=false)
	display(FinalTopology)
end

max_stars_repo_name: ricardoaf/jutop

max_issues_repo_name: ricardoaf/jutop

max_forks_repo_name: ricardoaf/jutop