text:
module Test_fdigenspec

using FaultDetectionTools
using DescriptorSystems
using LinearAlgebra
using Polynomials
using Random
using Test

# Testing examples for fdigenspec
println("Test_fdigenspec")
@testset "fdigenspec & fdichkspec" begin

p = 4; mu = 0; md = 0; mf = 4; n = 0;
sysf = fdimodset(rss(n,p,mf),faults = 1:mf);
@time S_weak = fdigenspec(sysf, atol = 1.e-7, FDtol = 1.e-5)
@test size(S_weak,1) == 15

@time S_strong = fdigenspec(sysf, atol = 1.e-7, FDtol = 1.e-5, FDGainTol = 1.e-3, FDfreq=[0,1])
@test S_weak == S_strong

@time rdims, orders, leastorders = fdichkspec(sysf, S_weak; atol = 1.e-7, FDtol = 1.e-5)
@test all(rdims .> 0) && all(orders .== 0)
@time rdims, orders, leastorders = fdichkspec(sysf, S_weak; atol = 1.e-7, FDtol = 1.e-5, minimal = true)
@test all(rdims .> 0) && all(orders .== 0) && all(leastorders .== 0)

## Example without control and disturbance inputs
p = 4; mu = 0; md = 0; mf = 4; n = 2;
sysf = fdimodset(rss(n,p,mf),faults = 1:mf);
@time S_weak = fdigenspec(sysf, atol = 1.e-7, FDtol = 1.e-5)
@test size(S_weak,1) == 15

@time S_strong = fdigenspec(sysf, atol = 1.e-7, FDtol = 1.e-5, FDGainTol = 1.e-3, FDfreq=[0,1])
@test S_weak == S_strong

@time rdims, orders, leastorders = fdichkspec(sysf, S_weak; atol = 1.e-7, FDtol = 1.e-5)
@test all(rdims .> 0) && all(orders .== n)

@time rdims, orders, leastorders = fdichkspec(sysf, S_weak; atol = 1.e-7, FDtol = 1.e-5, minimal = true)
@test all(rdims .> 0) && all(orders .== n) && all(leastorders .<= n)

## Yuan et al. IJC (1997)  weak synthesis with pole assignment and
#  enforcing least global order of 6!
p = 3; mu = 1; md = 0; mf = 8; n = 4;
A = [
-1 1 0 0
1 -2 1 0
0 1 -2 1
0 0 1 -2  ];
Bu = [1 0 0 0]';
Bf = [
1 0 0 0 1 0 0 0
0 1 0 0 -1 1 0 0
0 0 1 0 0 -1 1 0
0 0 0 1 0 0 -1 1];
C = [
 1 0 0 0
0 0 1 0
0 0 0 1];
Du = zeros(p,mu); Df = zeros(p,mf);

sys = dss(A,[Bu Bf],C,[Du Df]);
sysf = fdimodset(sys, c = 1:mu, d = mu .+ (1:md) , f = (mu+md) .+ (1:mf));

@time S_weak = fdigenspec(sysf, atol = 1.e-7, FDtol = 1.e-5)
@test size(S_weak,1) == 18 &&
sort(binmat2dec(S_weak)) == sort(binmat2dec([
0   0   0   1   0   0   1   1
0   1   1   0   1   1   1   0
0   1   1   1   1   1   0   1
0   1   1   1   1   1   1   1
1   0   1   0   1   1   1   0
1   0   1   1   1   1   0   1
1   0   1   1   1   1   1   1
1   1   0   0   1   1   0   0
1   1   0   1   1   1   1   1
1   1   1   0   0   1   1   0
1   1   1   0   1   0   1   0
1   1   1   0   1   1   1   0
1   1   1   1   0   1   0   1
1   1   1   1   0   1   1   1
1   1   1   1   1   0   0   1
1   1   1   1   1   0   1   1
1   1   1   1   1   1   0   1
1   1   1   1   1   1   1   1
]))

@time S_strong = fdigenspec(sysf, atol = 1.e-7, FDtol = 1.e-5, FDGainTol = 1.e-3, FDfreq=[0,1])
@test size(S_strong,1) == 12 &&
sort(binmat2dec(S_strong)) == sort(binmat2dec([
      0   0   0   1   0   0   1   1
      0   1   1   0   1   1   1   0
      0   1   1   1   1   1   0   1
      0   1   1   1   1   1   1   1
      1   0   1   0   1   1   1   0
      1   0   1   1   1   1   0   1
      1   0   1   1   1   1   1   1
      1   1   0   0   1   1   0   0
      1   1   0   1   1   1   1   1
      1   1   1   0   1   1   1   0
      1   1   1   1   1   1   0   1
      1   1   1   1   1   1   1   1
   ]))

@time rdims, orders, leastorders = fdichkspec(sysf, S_weak; atol = 1.e-7, FDtol = 1.e-5, minimal = true)
@test all(rdims .> 0) && all(orders .<= n)

@time rdims, orders, leastorders = fdichkspec(sysf, S_weak; FDfreq=[0,1], atol = 1.e-7, FDtol = 1.e-5, minimal = true)
@test sort(binmat2dec(S_strong)) == sort(binmat2dec(S_weak[rdims .> 0,:]))

# some timing results
@time S_weak = fdigenspec(sysf, atol = 1.e-7, FDtol = 1.e-5);
Stest = dec2binmat(Vector(1:2^mf-1));
@time rdims, orders, leastorders = fdichkspec(sysf, Stest; atol = 1.e-7, FDtol = 1.e-5, minimal = true);
@test size(S_weak,1) == 18 && sort(binmat2dec(S_weak)) == sort(binmat2dec(Stest[rdims .> 0,:]))

# Example 5.4m - Solution of an EFDP using EFDSYN
# Uses DescriptorSystems.jl v1.0 (or later) FaultDetectionTools.jl V0.1 (or later)

# define s as an improper transfer function
s = rtf('s');
# define Gu(s) and Gd(s)
Gu = [(s+1)/(s-2); (s+2)/(s-3)];     # enter Gu(s)
Gd = [(s-1)/(s+2); 0];               # enter Gd(s)
Gf = [(s+1)/(s-2) 0; (s+2)/(s-3) 1]; # enter Gf(s)
p = 2; mu = 1; md = 1; mf = 2;       # set dimensions

# setup the synthesis model with faults with Gf(s) = [ Gu(s) [0;1]]
sysf = fdimodset(dss([Gu Gd]),c =1,d = 2,f = 1,fs = 2);

@time S_weak = fdigenspec(sysf, atol = 1.e-7, FDtol = 1.e-5)
@test S_weak == Bool[1 1]

@time S_strong = fdigenspec(sysf, atol = 1.e-7, FDtol = 1.e-5, FDGainTol = 1.e-3, FDfreq=[0,1])
@test S_strong == Bool[1 1]

Stest = dec2binmat(Vector(1:2^mf-1));

@time rdims, orders, leastorders = fdichkspec(sysf, Stest; atol = 1.e-7, FDtol = 1.e-5, minimal = true)
@test rdims == [0, 0, 1] && orders == [-1, -1, 1] && leastorders == [-1, -1, 1]

## Model Niemann 1998, Optim. Appl. Meth.
n = 5; mu = 2; md = 1; mf = 2; m = mu+md+mf; p = 5;
a = [
-0.0782 0.2939 0.0220 0.0208 -0.0291
0.0077 -0.0278 0.0015 -0.0015 0.0026
-1.3251 5.32608 -0.5263 0.2214 -0.4777
1.0809 -4.4452 0.3770 -0.4631 0.4032
2.1532 -8.6386 0.7811 -0.5745 0.7816]*1000;
bd = [
    0.0443
   -0.0042
    0.7910
   -0.6598
   -1.2881]*1000;
bu = [
    -0.007197 0.003005
0.003473 0.000284
1.218502 -0.032866
1.322502 0.020147
-0.082265 0.024388 ];
bf = bu;
c = eye(n); dd = zeros(n,md); du = zeros(n,mu); df = zeros(n,mf);

sys = dss(a,[bu bd bf],c,[du dd df]);
sysf = fdimodset(sys, c = 1:mu, d = mu .+ (1:md) , f = (mu+md) .+ (1:mf));

@time S_weak = fdigenspec(sysf, atol = 1.e-7, FDtol = 1.e-5)
@test S_weak == Bool[1 1; 0 1; 1 0]

@time S_strong = fdigenspec(sysf, atol = 1.e-7, FDtol = 1.e-5, FDGainTol = 1.e-3, FDfreq=[0,1])
@test S_strong == Bool[1 1; 0 1; 1 0]

Stest = dec2binmat(Vector(1:2^mf-1));

@time rdims, orders, leastorders = fdichkspec(sysf, Stest; atol = 1.e-7, FDtol = 1.e-5, minimal = true)
@test rdims == [3, 3, 4] && orders == [3, 3, 4] && leastorders == [1, 1, 1]

## Example 1, Xiong & Saif 2000, Int. J. Robust Nonlinear Control
n = 4; mu = 2; md = 2; mf = 2; m = mu+md+mf; p = 4;
a = [ -9.9477 -0.7476 0.2632 -5.0337
      52.1659 2.7452 5.5532 -24.4221
      26.0922 2.6361 -4.1975 -19.2774
      0.0 0.0 1.0 0.0];
bu = [ 0.4422 0.1761
       3.5446 -7.5922
        -5.52 4.49
        0.0 0.0];
bd = [ 0 0; 0 1 ; 1 0; 0 0 ];
bf = [bu[:,1] zeros(n,1)];
c = [ 1 0 0 0; 0 1 0 0; 0 0 1 0; 0 1 1 1];
dd = zeros(n,md); du = zeros(n,mu); df = [zeros(n,mf-1) [0 0 1 0]'];

sys = dss(a,[bu bd bf],c,[du dd df]);
sysf = fdimodset(sys, c = 1:mu, d = mu .+ (1:md) , f = (mu+md) .+ (1:mf));

@time S_weak = fdigenspec(sysf, atol = 1.e-7, FDtol = 1.e-5)
@test S_weak == Bool[1 1; 0 1; 1 0]

@time S_strong = fdigenspec(sysf, atol = 1.e-7, FDtol = 1.e-5, FDGainTol = 1.e-3, FDfreq=[0,1])
@test S_strong == Bool[1 1; 0 1; 1 0]

Stest = dec2binmat(Vector(1:2^mf-1));

@time rdims, orders, leastorders = fdichkspec(sysf, Stest; atol = 1.e-7, FDtol = 1.e-5, minimal = true)
@test rdims == [1, 1, 2] && orders == [1, 2, 2] && leastorders == [1, 2, 1]

## Example 2, Xiong & Saif 2000, Int. J. Robust Nonlinear Control
n = 5; mu = 3; md = 2; mf = 2; m = mu+md+mf; p = 3;
a = [
0 0 -0.0034 0 0
0 -0.041 0.0013 0 0
0 0 -1.1471 0 0
0 0 -0.0036 0 0
0 0.094 0.0057 0 -0.051];
bu = [
    -1 0 0
     0 0 0
     0 0 0.948
     0.916 -1 0
    -0.598 0 0];
bd = [
     0 1
     0.062 -0.132
     0 -7.189
     0 0
     0 0 ];
bf = bu[:,1:2];
c = [
1 0 0 0 0
0 0 1 0 0
0 0 0 1 0];
dd = zeros(p,md); du = zeros(p,mu); df = zeros(p,mu-1);

sys = dss(a,[bu bd bf],c,[du dd df]);

sysf = fdimodset(sys, c = 1:mu, d = mu .+ (1:md) , f = (mu+md) .+ (1:mf));

@time S_weak = fdigenspec(sysf, atol = 1.e-7, FDtol = 1.e-5)
@test S_weak == Bool[1 1; 0 1; 1 0]

@time S_strong = fdigenspec(sysf, atol = 1.e-7, FDtol = 1.e-5, FDGainTol = 1.e-3, FDfreq=[0,1])
@test S_strong == Bool[1 1; 0 1; 1 0]

Stest = dec2binmat(Vector(1:2^mf-1));

@time rdims, orders, leastorders = fdichkspec(sysf, Stest; atol = 1.e-7, FDtol = 1.e-5, minimal = true)
@test rdims == [1, 1, 2] && orders == [1, 1, 2] && leastorders == [1, 1, 1]

## Linearized Boeing-747 (Varga,2007)

A = [
-0.4861 0.000317 -0.5588 0 -2.04e-6
0 -0.0199 3.0796 -9.8048 8.98e-5
1.0053 -0.0021 -0.5211 0 9.30e-6
1 0 0 0 0
0 0 -92.6 92.6 0];

Bu = [
-0.1455 -0.1455 -0.1494 -0.1494 -1.2860 0.0013 0.0035 0.0035 0.0013
0 0 0 0 -0.3122 0.1999 0.1999 0.1999 0.1999
-0.0071 -0.0071 -0.0074 -0.0074 -0.0676 -0.0004 -0.0004 -0.0004 -0.0004
0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0];

C = [
0 0 1 0 0
0 -0.0199 3.0796 -9.8048 8.98e-5
0 0 0 1 0
1 0 0 0 0
0 0 -92.6 92.6 0
0 0 0 0 1 ];

Du = [
 0 0 0 0 0 0 0 0 0
0 0 0 0 -0.3122 0.1999 0.1999 0.1999 0.1999
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  ];

Bf = Bu[:,1:5]; Df = Du[:,1:5];

p = 6; n = 5; mu = 9; mf = 5; md = 0;

sys = dss(A,[Bu Bf],C,[Du Df]);

sysf = fdimodset(sys, c = 1:mu, d = mu .+ (1:md) , f = (mu+md) .+ (1:mf));

@time S_weak = fdigenspec(sysf, atol = 1.e-7, FDtol = 1.e-5)
@test sort(binmat2dec(S_weak)) == sort(binmat2dec([
      0   0   0   0   1
      0   0   1   1   0
      0   0   1   1   1
      1   1   0   0   0
      1   1   0   0   1
      1   1   1   1   0
      1   1   1   1   1
      ]))

@time S_strong = fdigenspec(sysf, atol = 1.e-8, FDtol = 1.e-5, FDGainTol = 0.00001, FDfreq=[0,1])
@test sort(binmat2dec(S_strong)) == sort(binmat2dec([
    1  1  1  1  1
    0  0  1  1  1
    1  1  0  0  1
    1  1  1  1  0
         ]))

Stest = dec2binmat(Vector(1:2^mf-1));

@time rdims, orders, leastorders = fdichkspec(sysf, Stest; atol = 1.e-8, FDtol = 0.00001, minimal = true)
@test rdims[rdims .> 0] == [4, 4, 5, 4, 5, 5, 6] && orders[rdims .> 0] == [3, 4, 4, 4, 4, 5, 5] &&
leastorders[rdims .> 0] == [1, 2, 1, 2, 1, 1, 1]

@time rdims, orders, leastorders = fdichkspec(sysf, Stest; atol = 1.e-8, FDtol = 0.00001, FDGainTol = 0.00001, FDfreq=[0,1], minimal = true)
@test rdims[rdims .> 0] == [5, 5, 5, 6] && orders[rdims .> 0] == [4, 4, 5, 5] &&
leastorders[rdims .> 0] == [1, 1, 1, 1]

@time rdims, orders, leastorders = fdichkspec(sysf, S_weak; atol = 1.e-8, FDtol = 0.000001, FDGainTol = 0.00001, FDfreq=0)
@test rdims[rdims .> 0] == [6, 5, 5, 5] && orders[rdims .> 0] == [5, 4, 4, 5]

@time rdims, orders, leastorders = fdichkspec(sysf, S_weak; atol = 1.e-8, FDtol = 0.000001, FDGainTol = 0.00001, FDfreq=0, minimal = true)
@test rdims[rdims .> 0] == [6, 5, 5, 5] && orders[rdims .> 0] == [5, 4, 4, 5] && leastorders[rdims .> 0] == [1, 1, 1, 1]

S_weak[sortperm(binmat2dec(S_weak)),:]   # reordering according to increasing decimal values of rows
S_weak[sortperm(binmat2dec(S_weak[sum(S_weak,dims=2)[:],:])),:]  # reording according to increasing numbers of zeros in rows

# Lecture 6 example ()
s = rtf('s'); # define the Laplace variable s
Gu = [1/s; 1/s] # enter Gu(s)
Gd = [0; s/(s+3)]; # enter Gd(s)
Gf = [(s+1)/(s+2); 1/(s+2)]; # enter Gd(s)
# build state space model of [Gu(s) Gd(s) Gf(s)] and set input groups
sys = dss([Gu Gd Gf],minimal=true)
sysf = fdimodset(dss([Gu Gd Gf],minimal=true),c =1,d = 2,f = 3)

@time S_weak = fdigenspec(sysf, atol = 1.e-7, FDtol = 1.e-5)
@test S_weak ==  dec2binmat(1)

@time S_strong = fdigenspec(sysf, atol = 1.e-7, FDtol = 1.e-5, FDGainTol = 0.00001, FDfreq=[0,1])
@test isempty(S_strong)

@time rdims, ords, lords = fdichkspec(sysf,trues(1,1), FDfreq = 0)
@test rdims[1] == 0

s = rtf('s')
sys = dss([1/(s^2+1) 1/s]);
@time S_weak = fdigenspec(sys, atol = 1.e-7, FDtol = 1.e-5)
@test S_weak ==  Bool[1 1]

@time S_strong = fdigenspec(sys, atol = 1.e-7, FDtol = 1.e-5, FDGainTol = 0.00001, FDfreq=[0,1])
@test isempty(S_strong)

@time rdims, ords, lords = fdichkspec(sys, S_weak, FDfreq=[1])
@test rdims[1] == 0

s = rtf('s')
sys = dss([1/s/(s^2+1) 1/s/(s^2+1)/(s+2)]);
@time S_weak = fdigenspec(sys, atol = 1.e-7, FDtol = 1.e-5)
@test S_weak ==  Bool[1 1]

@time S_strong = fdigenspec(sys, atol = 1.e-7, FDtol = 1.e-5, FDGainTol = 0.00001, FDfreq=[0,1])
@test S_weak == S_strong

@time rdims, ords, lords = fdichkspec(sys, S_weak, FDfreq=[0,1])
@test rdims[1] == 1

## Example with solvable strong synthesis (frequencies are also zeros)
#
s = rtf('s');
# define 2x2 Gf(s) without zeros
Gf = [s/(s+1) 0; (s^2+4)/(s+3)^2 1]; # enter Gf(s)

@time sysf = fdimodset(dss(Gf; minimal = true, atol = 1.e-7),faults = 1:2);
S_weak = fdigenspec(sysf)
S_strong = fdigenspec(sysf, FDfreq = [0], FDtol = 1.e-7, FDGainTol = 1.e-3)
r, ord, lord = fdichkspec(sysf, S_weak; FDfreq = [0], FDGainTol = 1.e-3)
@test S_strong == S_weak[r.> 0,:]

S_strong = fdigenspec(sysf, FDfreq = [2], FDtol = 1.e-7, FDGainTol = 1.e-3)
r, ord, lord = fdichkspec(sysf, S_weak; FDfreq = [2], FDGainTol = 1.e-3)
@test S_strong == S_weak[r.> 0,:]

S_strong = fdigenspec(sysf, FDfreq = [0,2], FDtol = 1.e-7, FDGainTol = 1.e-3)
r, ord, lord = fdichkspec(sysf, S_weak; FDfreq = [0,2], FDGainTol = 1.e-3)
@test S_strong == S_weak[r.> 0,:]

end # test fdigenspec

end # module

max_stars_repo_name: andreasvarga/FaultDetectionTools

max_issues_repo_name: andreasvarga/FaultDetectionTools

max_forks_repo_name: andreasvarga/FaultDetectionTools


text:
module FlowNodeTest

using Test
using ReactiveMP
using ReactiveMP: getL, getα, getβ, getκ, getλ

include("../../test_helpers.jl")

using .ReactiveMPTestingHelpers

@testset "FlowNode" begin

    # test for basics in the flow.jl file
    @testset "Basics" begin

        # tests for node creation
        @testset "Creation" begin

            # create example Flow node
            node = make_node(Flow, meta = 1)

            # check whether Flow node creation has been a success
            @test functionalform(node)          === Flow
            @test sdtype(node)                  === Deterministic()
            @test name.(interfaces(node))       === (:out, :in)
            @test factorisation(node)           === ((1, 2), )
            @test localmarginalnames(node)      === (:out_in, )
            @test metadata(node)                === 1

        end

        # tests for meta data
        @testset "MetaData" begin

            @testset "Default" begin

                # create example meta data
                meta = FlowMeta(compile(FlowModel(2, (AdditiveCouplingLayer(PlanarFlow()),))))

                # check whether the Flow node has a default metadata structure
                @test_throws ErrorException ReactiveMP.default_meta(Flow)

                # check whether the getmodel function works
                @test typeof(getmodel(meta)) <: ReactiveMP.AbstractCompiledFlowModel
                @test getmodel(meta) == meta.model
                @test getapproximation(meta) == meta.approximation
                @test getapproximation(meta) == Linearization()
            end

            @testset "Linearization" begin

                # create example meta data
                meta = FlowMeta(compile(FlowModel(2, (AdditiveCouplingLayer(PlanarFlow()),))), Linearization())

                # check whether the Flow node has a default metadata structure
                @test_throws ErrorException ReactiveMP.default_meta(Flow)

                # check whether the getmodel function works
                @test typeof(getmodel(meta)) <: ReactiveMP.AbstractCompiledFlowModel
                @test getmodel(meta) == meta.model
                @test getapproximation(meta) == meta.approximation
                @test typeof(getapproximation(meta)) == Linearization

            end

            @testset "Unscented" begin

                # create example meta data
                meta = FlowMeta(compile(FlowModel(2, (AdditiveCouplingLayer(PlanarFlow()),))), Unscented(3))

                # check whether the Flow node has a default metadata structure
                @test_throws ErrorException ReactiveMP.default_meta(Flow)

                # check whether the getmodel function works
                @test typeof(getmodel(meta)) <: ReactiveMP.AbstractCompiledFlowModel
                @test getmodel(meta) == meta.model
                @test getapproximation(meta) == meta.approximation
                @test typeof(getapproximation(meta)) == Unscented
                @test getL(getapproximation(meta))   == 3
                @test getα(getapproximation(meta))   == 1e-3
                @test getβ(getapproximation(meta))   == 2.0
                @test getκ(getapproximation(meta))   == 0.0
                @test getλ(getapproximation(meta))   == 3e-6 - 3

            end

        end

    end

    # test for model specifics in the flow_models folder
    @testset "Models" begin

        include("flow_models/test_flow_model.jl")

    end

    # test for layer specifics in the layers folder
    @testset "Layers" begin

        include("layers/test_additive_coupling_layer.jl")
        include("layers/test_input_layer.jl")
        include("layers/test_permutation_layer.jl")

    end

    # test for coupling flow specifics in the coupling_flows folder
    @testset "Coupling flows" begin

        include("coupling_flows/test_planar_flow.jl")
        include("coupling_flows/test_radial_flow.jl")

    end

end

end

ext: jl


text:
# syntax: proto3
using ProtoBuf
import ProtoBuf.meta

const DataClass = (;[
    Symbol("DATA_CLASS_UNKNOWN") => Int32(0),
    Symbol("DATA_CLASS_SCALAR") => Int32(1),
    Symbol("DATA_CLASS_TENSOR") => Int32(2),
    Symbol("DATA_CLASS_BLOB_SEQUENCE") => Int32(3),
]...)

mutable struct SummaryDescription <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function SummaryDescription(; kwargs...)
        obj = new(meta(SummaryDescription), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end # mutable struct SummaryDescription
const __meta_SummaryDescription = Ref{ProtoMeta}()
function meta(::Type{SummaryDescription})
    ProtoBuf.metalock() do
        if !isassigned(__meta_SummaryDescription)
            __meta_SummaryDescription[] = target = ProtoMeta(SummaryDescription)
            allflds = Pair{Symbol,Union{Type,String}}[:type_hint => AbstractString]
            meta(target, SummaryDescription, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_SummaryDescription[]
    end
end
function Base.getproperty(obj::SummaryDescription, name::Symbol)
    if name === :type_hint
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    else
        getfield(obj, name)
    end
end

mutable struct HistogramProto <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function HistogramProto(; kwargs...)
        obj = new(meta(HistogramProto), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end # mutable struct HistogramProto
const __meta_HistogramProto = Ref{ProtoMeta}()
function meta(::Type{HistogramProto})
    ProtoBuf.metalock() do
        if !isassigned(__meta_HistogramProto)
            __meta_HistogramProto[] = target = ProtoMeta(HistogramProto)
            pack = Symbol[:bucket_limit,:bucket]
            allflds = Pair{Symbol,Union{Type,String}}[:min => Float64, :max => Float64, :num => Float64, :sum => Float64, :sum_squares => Float64, :bucket_limit => Base.Vector{Float64}, :bucket => Base.Vector{Float64}]
            meta(target, HistogramProto, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, pack, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_HistogramProto[]
    end
end
function Base.getproperty(obj::HistogramProto, name::Symbol)
    if name === :min
        return (obj.__protobuf_jl_internal_values[name])::Float64
    elseif name === :max
        return (obj.__protobuf_jl_internal_values[name])::Float64
    elseif name === :num
        return (obj.__protobuf_jl_internal_values[name])::Float64
    elseif name === :sum
        return (obj.__protobuf_jl_internal_values[name])::Float64
    elseif name === :sum_squares
        return (obj.__protobuf_jl_internal_values[name])::Float64
    elseif name === :bucket_limit
        return (obj.__protobuf_jl_internal_values[name])::Base.Vector{Float64}
    elseif name === :bucket
        return (obj.__protobuf_jl_internal_values[name])::Base.Vector{Float64}
    else
        getfield(obj, name)
    end
end

mutable struct SummaryMetadata_PluginData <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function SummaryMetadata_PluginData(; kwargs...)
        obj = new(meta(SummaryMetadata_PluginData), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end # mutable struct SummaryMetadata_PluginData
const __meta_SummaryMetadata_PluginData = Ref{ProtoMeta}()
function meta(::Type{SummaryMetadata_PluginData})
    ProtoBuf.metalock() do
        if !isassigned(__meta_SummaryMetadata_PluginData)
            __meta_SummaryMetadata_PluginData[] = target = ProtoMeta(SummaryMetadata_PluginData)
            allflds = Pair{Symbol,Union{Type,String}}[:plugin_name => AbstractString, :content => Vector{UInt8}]
            meta(target, SummaryMetadata_PluginData, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_SummaryMetadata_PluginData[]
    end
end
function Base.getproperty(obj::SummaryMetadata_PluginData, name::Symbol)
    if name === :plugin_name
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    elseif name === :content
        return (obj.__protobuf_jl_internal_values[name])::Vector{UInt8}
    else
        getfield(obj, name)
    end
end

mutable struct SummaryMetadata <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function SummaryMetadata(; kwargs...)
        obj = new(meta(SummaryMetadata), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end # mutable struct SummaryMetadata
const __meta_SummaryMetadata = Ref{ProtoMeta}()
function meta(::Type{SummaryMetadata})
    ProtoBuf.metalock() do
        if !isassigned(__meta_SummaryMetadata)
            __meta_SummaryMetadata[] = target = ProtoMeta(SummaryMetadata)
            allflds = Pair{Symbol,Union{Type,String}}[:plugin_data => SummaryMetadata_PluginData, :display_name => AbstractString, :summary_description => AbstractString, :data_class => Int32]
            meta(target, SummaryMetadata, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_SummaryMetadata[]
    end
end
function Base.getproperty(obj::SummaryMetadata, name::Symbol)
    if name === :plugin_data
        return (obj.__protobuf_jl_internal_values[name])::SummaryMetadata_PluginData
    elseif name === :display_name
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    elseif name === :summary_description
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    elseif name === :data_class
        return (obj.__protobuf_jl_internal_values[name])::Int32
    else
        getfield(obj, name)
    end
end

mutable struct Summary_Image <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function Summary_Image(; kwargs...)
        obj = new(meta(Summary_Image), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end # mutable struct Summary_Image
const __meta_Summary_Image = Ref{ProtoMeta}()
function meta(::Type{Summary_Image})
    ProtoBuf.metalock() do
        if !isassigned(__meta_Summary_Image)
            __meta_Summary_Image[] = target = ProtoMeta(Summary_Image)
            allflds = Pair{Symbol,Union{Type,String}}[:height => Int32, :width => Int32, :colorspace => Int32, :encoded_image_string => Vector{UInt8}]
            meta(target, Summary_Image, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_Summary_Image[]
    end
end
function Base.getproperty(obj::Summary_Image, name::Symbol)
    if name === :height
        return (obj.__protobuf_jl_internal_values[name])::Int32
    elseif name === :width
        return (obj.__protobuf_jl_internal_values[name])::Int32
    elseif name === :colorspace
        return (obj.__protobuf_jl_internal_values[name])::Int32
    elseif name === :encoded_image_string
        return (obj.__protobuf_jl_internal_values[name])::Vector{UInt8}
    else
        getfield(obj, name)
    end
end

mutable struct Summary_Audio <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function Summary_Audio(; kwargs...)
        obj = new(meta(Summary_Audio), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end # mutable struct Summary_Audio
const __meta_Summary_Audio = Ref{ProtoMeta}()
function meta(::Type{Summary_Audio})
    ProtoBuf.metalock() do
        if !isassigned(__meta_Summary_Audio)
            __meta_Summary_Audio[] = target = ProtoMeta(Summary_Audio)
            allflds = Pair{Symbol,Union{Type,String}}[:sample_rate => Float32, :num_channels => Int64, :length_frames => Int64, :encoded_audio_string => Vector{UInt8}, :content_type => AbstractString]
            meta(target, Summary_Audio, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_Summary_Audio[]
    end
end
function Base.getproperty(obj::Summary_Audio, name::Symbol)
    if name === :sample_rate
        return (obj.__protobuf_jl_internal_values[name])::Float32
    elseif name === :num_channels
        return (obj.__protobuf_jl_internal_values[name])::Int64
    elseif name === :length_frames
        return (obj.__protobuf_jl_internal_values[name])::Int64
    elseif name === :encoded_audio_string
        return (obj.__protobuf_jl_internal_values[name])::Vector{UInt8}
    elseif name === :content_type
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    else
        getfield(obj, name)
    end
end

mutable struct Summary_Value <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function Summary_Value(; kwargs...)
        obj = new(meta(Summary_Value), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end # mutable struct Summary_Value
const __meta_Summary_Value = Ref{ProtoMeta}()
function meta(::Type{Summary_Value})
    ProtoBuf.metalock() do
        if !isassigned(__meta_Summary_Value)
            __meta_Summary_Value[] = target = ProtoMeta(Summary_Value)
            fnum = Int[7,1,9,2,3,4,5,6,8]
            allflds = Pair{Symbol,Union{Type,String}}[:node_name => AbstractString, :tag => AbstractString, :metadata => SummaryMetadata, :simple_value => Float32, :obsolete_old_style_histogram => Vector{UInt8}, :image => Summary_Image, :histo => HistogramProto, :audio => Summary_Audio, :tensor => TensorProto]
            oneofs = Int[0,0,0,1,1,1,1,1,1]
            oneof_names = Symbol[Symbol("value")]
            meta(target, Summary_Value, allflds, ProtoBuf.DEF_REQ, fnum, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, oneofs, oneof_names)
        end
        __meta_Summary_Value[]
    end
end
function Base.getproperty(obj::Summary_Value, name::Symbol)
    if name === :node_name
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    elseif name === :tag
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    elseif name === :metadata
        return (obj.__protobuf_jl_internal_values[name])::SummaryMetadata
    elseif name === :simple_value
        return (obj.__protobuf_jl_internal_values[name])::Float32
    elseif name === :obsolete_old_style_histogram
        return (obj.__protobuf_jl_internal_values[name])::Vector{UInt8}
    elseif name === :image
        return (obj.__protobuf_jl_internal_values[name])::Summary_Image
    elseif name === :histo
        return (obj.__protobuf_jl_internal_values[name])::HistogramProto
    elseif name === :audio
        return (obj.__protobuf_jl_internal_values[name])::Summary_Audio
    elseif name === :tensor
        return (obj.__protobuf_jl_internal_values[name])::TensorProto
    else
        getfield(obj, name)
    end
end

mutable struct Summary <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function Summary(; kwargs...)
        obj = new(meta(Summary), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end # mutable struct Summary
const __meta_Summary = Ref{ProtoMeta}()
function meta(::Type{Summary})
    ProtoBuf.metalock() do
        if !isassigned(__meta_Summary)
            __meta_Summary[] = target = ProtoMeta(Summary)
            allflds = Pair{Symbol,Union{Type,String}}[:value => Base.Vector{Summary_Value}]
            meta(target, Summary, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_Summary[]
    end
end
function Base.getproperty(obj::Summary, name::Symbol)
    if name === :value
        return (obj.__protobuf_jl_internal_values[name])::Base.Vector{Summary_Value}
    else
        getfield(obj, name)
    end
end

export DataClass, SummaryDescription, HistogramProto, SummaryMetadata_PluginData, SummaryMetadata, Summary_Image, Summary_Audio, Summary_Value, Summary


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

struct MCMCSampleID
    chainid::Int32
    chaincycle::Int32
    stepno::Int64
    sampletype::Int64
end

struct MCMCSampleIDVector{
    VInt32<:AbstractVector{Int32},
    VInt64<:AbstractVector{Int64}
} <: BATDataVector{MCMCSampleID}
    chainid::VInt32
    chaincycle::VInt32
    stepno::VInt64
    sampletype::VInt64
end

export MCMCSampleIDVector

MCMCSampleIDVector() = MCMCSampleIDVector(Vector{Int32}(), Vector{Int32}(), Vector{Int64}(), Vector{Int64}())

MCMCSampleIDVector(chain::MCMCIterator) = MCMCSampleIDVector()

Base.size(xs::MCMCSampleIDVector) = size(xs.chainid)

Base.@propagate_inbounds Base.getindex(xs::MCMCSampleIDVector, i::Integer)  =
    MCMCSampleID(xs.chainid[i], xs.chaincycle[i], xs.stepno[i], xs.sampletype[i])

Base.@propagate_inbounds function Base.setindex!(xs::MCMCSampleIDVector, x::MCMCSampleID, i::Integer)
    xs.chainid[i] = x.chainid
    xs.chaincycle[i] = x.chaincycle
    xs.stepno[i] = x.stepno
    xs.sampletype[i] = x.sampletype
    xs
end

Base.IndexStyle(xs::MCMCSampleIDVector) = IndexStyle(xs.chainid)

function Base.push!(xs::MCMCSampleIDVector, x::MCMCSampleID)
    push!(xs.chainid, x.chainid)
    push!(xs.chaincycle, x.chaincycle)
    push!(xs.stepno, x.stepno)
    push!(xs.sampletype, x.sampletype)
    xs
end

function Base.append!(A::MCMCSampleIDVector, B::MCMCSampleIDVector)
    append!(A.chainid, B.chainid)
    append!(A.chaincycle, B.chaincycle)
    append!(A.stepno, B.stepno)
    append!(A.sampletype, B.sampletype)
    A
end

function Base.append!(A::MCMCSampleIDVector, B::AbstractVector{<:MCMCSampleID})
    for x in B
        push!(A, x)
    end
    A
end

Base.@propagate_inbounds function Base.view(A::MCMCSampleIDVector, idxs)
    MCMCSampleIDVector(
        view(A.chainid, idxs),
        view(A.chaincycle, idxs),
        view(A.stepno, idxs),
        view(A.sampletype, idxs)
    )
end

Base.convert(::Type{AbstractMCMCCallback}, x::MCMCSampleIDVector) = MCMCAppendCallback(x)

MCMCAppendCallback(x::MCMCSampleIDVector, nonzero_weights::Bool = true) =
    MCMCAppendCallback(x, 1, get_sample_ids!, nonzero_weights)

function read_fom_hdf5(input, ::Type{MCMCSampleIDVector})
    MCMCSampleIDVector(
        input["chainid"][:],
        input["chaincycle"][:],
        input["stepno"][:],
        input["sampletype"][:],
    )
end

function write_to_hdf5(output, sampleids::MCMCSampleIDVector)
    output["chainid"] = sampleids.chainid
    output["chaincycle"] = sampleids.chaincycle
    output["stepno"] = sampleids.stepno
    output["sampletype"] = sampleids.sampletype
    nothing
end


text:
n = 10 # rand(1:10)
A = matrixdepot("randcorr", n)
@test issymmetric(A) # symmetric

# positive semidefinite
for eig in eigvals(A)
    @test eig >= 0
end

# 1s on the diagonal
@test diag(A) ≈ ones(Float64, n)
println("'randcorr' passed test...")


text:
module EQmap

    # Modules
    using Plots, DataFrames

    # Methods
    include("gen_map.jl")

end


text:
using Pkg
Pkg.activate(".")
Pkg.build()
Pkg.test(; coverage=true)


text:
@testset "Measure" begin
    @testset "AABox" begin
        for T in Floats
            @test measure(setup_AABox2(T)) ≈ 1
            @test measure(setup_AABox3(T)) ≈ 1
        end
    end

    @testset "LineSegment" begin
        for T in Floats
            @test measure(setup_LineSegment2(T)) ≈ sqrt(T(2))
            @test measure(setup_LineSegment3(T)) ≈ sqrt(T(3))
        end
    end

    @testset "QuadraticSegment" begin
        for T in Floats
            q₂ = QuadraticSegment(Point{2,T}(0,0), Point{2,T}(1,0), Point{2,T}(1//2, 0))
            @test measure(q₂) ≈ 1
            q₂ = setup_QuadraticSegment2(T)
            @test abs(measure(q₂) - 1.4789428575445974) < 1e-6
            q₃ = QuadraticSegment(Point{3,T}(0,0,0),
                                  Point{3,T}(1,0,0),
                                  Point{3,T}(1//2,0,0))
            @test measure(q₃) ≈ 1
            q₃ = setup_QuadraticSegment3(T)
            @test abs(measure(q₃) - 1.7978527887818835) < 1e-6
        end
    end

    @testset "Polygon" begin
        for T in Floats
            @test measure(setup_Triangle2(T)) ≈ 1//2
            @test measure(setup_Triangle3(T)) ≈ 1/√(T(2))
            @test abs(measure(setup_Quadrilateral2(T)) - 1) < 1e-6
            @test abs(measure(setup_Quadrilateral3(T)) - 1.28078927557) < 1e-6
        end
    end

    @testset "QuadraticPolygon" begin
        for T in Floats
            @test abs(measure(setup_QuadraticTriangle2(T)) - 0.4333333333) < 1e-6
            @test abs(measure(setup_QuadraticTriangle3(T)) - 0.9058066937) < 1e-6
            @test abs(measure(setup_QuadraticQuadrilateral2(T)) - 0.86666666) < 1e-6
            @test abs(measure(setup_QuadraticQuadrilateral3(T)) - 1.39710509) < 1e-6
        end
    end
end

max_stars_repo_name: khurrumsaleem/MOCNeutronTransport

max_issues_repo_name: khurrumsaleem/MOCNeutronTransport

max_forks_repo_name: khurrumsaleem/MOCNeutronTransport


text:
function load_vfa_dicom(dir)
    dicoms = load_dicom_dir(dir)
    num_images = length(dicoms)
    unique_flip_angles = unique(lookup.(dicoms, "Flip Angle"))
    num_flip_angles = length(unique_flip_angles)
    num_slices = Int(num_images / num_flip_angles)

    dummy_image = lookup(dicoms[1], "Pixel Data")
    image_size = size(dummy_image)
    signal = zeros(num_images, image_size...)
    flip_angles = zeros(num_images)
    for dicom in dicoms
        instance = lookup(dicom, "Instance Number")
        signal[:, :, instance] = lookup(dicom, "Pixel Data")
        flip_angles[instance] = lookup(dicom, "Flip Angle")
    end
    signal = reshape(signal, (image_size..., num_slices, num_flip_angles))
    flip_angles = reshape(flip_angles, (num_slices, num_flip_angles))[1, :]
    @. flip_angles = deg2rad(flip_angles)
    TR = lookup(dicoms[1], "Repetition Time")
    return (; signal, angles = flip_angles, TR)
end

function load_dce_dicom(dir; num_slices::Integer)
    dicom_files = load_dicom_dir(dir)
    num_images = length(dicom_files)
    num_timepoints = Int(num_images / num_slices)

    dummy_image = lookup(dicom_files[1], "Pixel Data")
    image_size = size(dummy_image)
    signal = zeros(image_size..., num_images)
    timepoints = zeros(num_images)
    for dicom in dicom_files
        instance = lookup(dicom, "Instance Number")
        signal[:, :, instance] = lookup(dicom, "Pixel Data")
        timepoints[instance] = lookup(dicom, "Trigger Time")
    end
    signal = reshape(signal, (image_size..., num_slices, num_timepoints))
    timepoints = reshape(timepoints, (num_slices, num_timepoints))[1, :] ./ 1000 ./ 60
    TR = lookup(dicom_files[1], "Repetition Time")
    flip_angle = deg2rad(lookup(dicom_files[1], "Flip Angle"))
    return (; signal, t = vec(timepoints), TR, angle = flip_angle)
end

function load_dicom_dir(dir)
    dicom_files = list_dicom_in_dir(dir)
    num_files = length(dicom_files)
    dummy_data = dcm_parse(dicom_files[1])
    dicom_data = Vector{typeof(dummy_data)}(undef, num_files)
    for (index, file) in enumerate(dicom_files)
        dicom_data[index] = dcm_parse(file)
    end
    return dicom_data
end

function list_dicom_in_dir(dir)
    files = readdir(dir)
    dicom_files = joinpath.(dir, files[is_dicom.(files)])
    return dicom_files
end

is_dicom(file::AbstractString) = splitext(file)[2] == ".dcm"


text:
function _versioncheck(version::VersionNumber, valid_versions)
    supported = version ∈ valid_versions
    supported || error(LOGGER, ArgumentError(
        string(
            "Unsupported version ($version). ",
            "Expected a value between ($valid_versions)."
        )
    ))
end

function _versioncheck(version::String, valid_versions)
    return _versioncheck(VersionNumber(version), valid_versions)
end

# Cache of the versioninfo and image, so we don't compute these every time.
const _CACHE = Dict(
    :MANIFEST => "",
    :PROJECT => "",
    :PKGS => Dict{String, VersionNumber}(),
    :IMAGE => "",
)

# Make sure _CACHE access is thread-safe
const _CACHE_LOCK = ReentrantLock()

function _pkgs()
    lock(_CACHE_LOCK) do
        if isempty(_CACHE[:PKGS])
            for (pkg, ver) in Pkg.installed()
                # BSON can't handle Void types
                if ver !== nothing
                    global _CACHE[:PKGS][pkg] = ver
                end
            end
        end

        return _CACHE[:PKGS]
    end
end

function _env()
    lock(_CACHE_LOCK) do
        if isempty(_CACHE[:PROJECT]) && isempty(_CACHE[:MANIFEST])
            _CACHE[:PROJECT] = read(Base.active_project(), String)
            _CACHE[:MANIFEST] = read(
                joinpath(dirname(Base.active_project()), "Manifest.toml"),
                String
            )
        end

        return (_CACHE[:PROJECT], _CACHE[:MANIFEST])
    end
end

function _image()
    if isempty(_CACHE[:IMAGE]) && haskey(ENV, "JLSO_IMAGE")
        return ENV["JLSO_IMAGE"]
    end

    return _CACHE[:IMAGE]
end

function Pkg.activate(jlso::JLSOFile, path=tempdir(); kwargs...)
    mkpath(path)
    open(io -> Pkg.TOML.print(io, jlso.project), joinpath(path, "Project.toml"), "w")
    open(io -> Pkg.TOML.print(io, jlso.manifest), joinpath(path, "Manifest.toml"), "w")
    Pkg.activate(path; kwargs...)
end

function Pkg.activate(f::Function, jlso::JLSOFile, path=tempdir(); kwargs...)
    curr_env = dirname(Base.active_project())
    try
        Pkg.activate(jlso, path; kwargs...)
        f()
    finally
        Pkg.activate(curr_env)
    end
end

max_stars_repo_name: invenia/JLSO

max_issues_repo_name: invenia/JLSO

max_forks_repo_name: invenia/JLSO


text:
module Utility

using Rematch
using StatsBase

# internal modules
using ..Intervals, ..Nodes

import ..Graphs:
    reverse_complement, reverse_complement!,
    SNPMap, InsMap, DelMap

# exports
export random_id, log
export Alignment
export cigar, uncigar
export hamming_align

export contiguous_trues
export make_consensus, alignment_alleles

export write_fasta, read_fasta, name
export read_paf

Maybe{T} = Union{Nothing,T}

# ------------------------------------------------------------------------
# random functions

# random string of fixed length
"""
	random_id(;len=10, alphabet=UInt8[])

Generate a random string of length `len` drawn from letters in `alphabet`.
"""
function random_id(;len=10, alphabet=UInt8[])
    if length(alphabet) == 0
        alphabet = ['A','B','C','D','E','F','G','H','I','J','K','L','M',
                    'N','O','P','Q','R','S','T','U','V','W','X','Y','Z']
    end
    return join(sample(alphabet, len))
end

# helper functions w/ common functionality
"""
	write_fasta(io::IO, name, seq)

Output a single FASTA record with sequence `seq` and name `name` to IO stream `io`.
"""
function write_fasta(io::IO, name, seq)
    write(io::IO, '>', name, '\n')
    write(io::IO, columns(seq), '\n')
end

# ------------------------------------------------------------------------
# banded alignment data structure

"""
	struct Score <: AbstractArray{Float64,2}
		rows::Int
		cols::Int
		band::NamedTuple{(:lower, :upper)}
		data::Array{Float64}
		offset::Array{Int}
		starts::Array{Int}
		stops::Array{Int}
	end

Store information about a banded pairwise alignment.
"""
struct Score <: AbstractArray{Float64,2}
    rows::Int
    cols::Int
    band::NamedTuple{(:lower, :upper)}
    data::Array{Float64}
    offset::Array{Int}
    starts::Array{Int}
    stops::Array{Int}
end

function Score(rows, cols; band=(lower=Inf,upper=Inf))
    num_elts = [(r == 0) ? 0 : min(r+band.upper, cols) - max(r-band.lower, 1) + 1 for r in 0:rows]
    offset   = round.(Int, cumsum(num_elts))
    starts   = round.([max(r-band.lower, 1) for r in 1:rows])
    stops    = round.([min(r+band.upper, cols) for r in 1:rows])
    return Score(rows, cols, band, zeros(offset[end]), offset, starts, stops)
end

# data stored in row-major form
# @inline index(S::Score, i, j) = (j-S.starts[i]+1) + S.offset[i]
@inline index(S::Score, i, j) = j + S.cols*(i-1)
@inline data(S::Score, i, j)  = S.data[index(S,i,j)]

# --------------------------------
# overloads of base operators

Base.size(S::Score)        = (S.rows, S.cols)
Base.ndims(T::Type{Score}) = 2

Base.BroadcastStyle(::Type{Score})              = Base.Broadcast.ArrayStyle{Score}()
Base.show(io::IO, S::Score)                     = Base.show(io, reshape(S.data, S.rows, S.cols))
Base.show(io::IO, ::MIME"text/plain", S::Score) = Base.show(io, reshape(S.data, S.rows, S.cols))
Base.display(S::Score)                          = Base.display(reshape(S.data, S.rows, S.cols))

function Base.lastindex(S::Score, d)
    if d == 1
        return S.rows
    elseif d == 2
        return S.cols
    else
        error("incorrect number of dimensions")
    end
end

# TODO: better (less-explicit) indexing
#       would be nice to keep as abstract ranges if we could, e.g. 1:2:10
#       right now we allocate memory every time
@inline Base.getindex(S::Score, inds...) = Base.getindex(S.data, index(S,r,c))
    # TODO: make this faster!
    # close(out)
    #       causes huge slow down
    # rows, cols = inds
    # ι = [index(S,r,c) for r in rows for c in cols]
    # if length(ι) == 1
    #     # if cols < S.starts[rows] || cols > S.stops[rows]
    #     #     println("out of bound access: (", rows, ",", cols, ")" )
    #     #     return -Inf
    #     # else
    #     return Base.getindex(S.data, ι[1])
    #     # end
    # else
    #     return Base.getindex(S.data, ι)
    # end
# end

function Base.setindex!(S::Score, X, inds...)
    i, j = inds
    Base.setindex!(S.data, X, index(S,i,j))
end

# ------------------------------------------------------------------------
# cigar functions

# alignment -> cigar
"""
	cigar(seq₁::Array{UInt8}, seq₂::Array{UInt8})

Given two sequences, `seq₁` and `seq₂`, perform a pairwise banded alignment and return the cigar string of alignment.
"""
function cigar(seq₁::Array{UInt8}, seq₂::Array{UInt8})
    if length(seq₁) != length(seq₂)
        error("not an alignment")
    end

    aln = IOBuffer()
    M, I, D = 0, 0, 0
    for (c₁, c₂) in zip(seq₁, seq₂)
        @match (Char(c₁), Char(c₂)) begin
            ('-','-') => error("both columns are gaps")
            ('-', _ ) => begin
                if I > 0
                    write(aln, "$(I)I")
                    I = 0
                elseif M > 0
                    write(aln, "$(M)M")
                    M = 0
                end
                D += 1
            end
            ( _ ,'-') => begin
                if D > 0
                    write(aln, "$(D)D")
                    D = 0
                elseif M > 0
                    write(aln, "$(M)M")
                    M = 0
                end
                I += 1
            end
            ( _ , _ ) => begin
                if D > 0
                    write(aln, "$(D)D")
                    D = 0
                elseif I > 0
                    write(aln, "$(I)I")
                    I = 0
                end
                M += 1
            end
        end
    end

    if I > 0
        write(aln, "$(I)I")
        I = 0
    elseif M > 0
        write(aln, "$(M)M")
        M = 0
    elseif D > 0
        write(aln, "$(D)D")
        D = 0
    end

    return String(take!(aln))
end

uncigar(cg::T) where T <: AbstractArray{Tuple{Int,Char}} = cg

"""
	uncigar(cg::String)

Return an interator over intervals of alignment defined by cigar string `cg`.
"""
function uncigar(cg::String)
    chan = Channel{Tuple{Int, Char}}(0)
    @async begin
        i₁, i₂ = 1, 1
        while i₁ <= length(cg)
            while isdigit(cg[i₂])
                i₂ += 1
            end
            put!(chan, (parse(Int,cg[i₁:i₂-1]),cg[i₂]))
            i₁ = i₂ + 1
            i₂ = i₁
        end
        close(chan)
    end

    return chan
end

# ------------------------------------------------------------------------
# alignment functions

# paf alignment pair
"""
	mutable struct Hit
		name::String
		length::Int
		start::Int
		stop::Int
		seq::Maybe{Array{UInt8,1}}
	end

Hit is one side of a pairwise alignment between homologous sequences.
"""
mutable struct Hit
    name::String
    length::Int
    start::Int
    stop::Int
    seq::Maybe{Array{UInt8,1}}
end

"""
	mutable struct Alignment{T <: Union{String,Nothing,Array{Tuple{Int,Char}}}}
		qry::Hit
		ref::Hit
		matches::Int
		length::Int
		quality::Int
		orientation::Bool
		cigar::T
		divergence::Union{Float64,Nothing}
		align::Union{Float64,Nothing}
	end

Alignment is a pairwise homologous alignment between two sequences.
"""
mutable struct Alignment{T <: Union{String,Nothing,Array{Tuple{Int,Char}}}}
    qry::Hit
    ref::Hit
    matches::Int
    length::Int
    quality::Int
    orientation::Bool
    cigar::T
    divergence::Union{Float64,Nothing}
    align::Union{Float64,Nothing}
end

# dynamic programming
"""
	align(seq₁::Array{UInt8}, seq₂::Array{UInt8}, cost::Score)

Perform a pairwise alignment using Needleman-Wunsch style dynamic programming between
`seq₁` and `seq₂` given `cost`. The cost is defined by the Score structure.
"""
function align(seq₁::Array{UInt8}, seq₂::Array{UInt8}, cost)
    flip = if length(seq₁) > length(seq₂)
        seq₁, seq₂ = seq₂, seq₁
        true
    else
        false
    end

    L₁, L₂ = length(seq₁)+1, length(seq₂)+1

    # initialize matrices
    score = (
         M = zeros(L₁,L₂), #Score(L₁,L₂,band=cost.band),
         I = zeros(L₁,L₂), #Score(L₁,L₂,band=cost.band),
         D = zeros(L₁,L₂), #Score(L₁,L₂,band=cost.band),
    )

    # NOTE: upper and lower could be flipped
    bound(x,d) = isinf(x) ? size(score.M,d) : min(size(score.M,d),x)
    for j in 1:bound(cost.band.upper,2)
        score.M[1,j] = cost.gap(j-1)
    end
    for i in 1:bound(cost.band.lower,1)
        score.M[i,1] = cost.gap(i-1)
    end

    score.I[1,1:bound(cost.band.upper,2)] .= -Inf
    score.D[1:bound(cost.band.lower,1),1] .= -Inf

    # fill in bulk
    for i in 2:L₁
        lb = round(Int, max(i-cost.band.lower, 2))
        ub = round(Int, min(i+cost.band.upper, L₂))
        for j in lb:ub
            score.I[i,j] = max(
               score.M[i-1,j]+cost.open,
               score.I[i-1,j]+cost.extend,
            )
            score.D[i,j] = max(
               score.M[i,j-1]+cost.open,
               score.D[i,j-1]+cost.extend,
            )
            score.M[i,j] = max(
               score.M[i-1,j-1]+cost.match(seq₁[i-1],seq₂[j-1]),
               score.I[i,j],
               score.D[i,j],
            )
        end
    end

    # backtrack
    a₁, a₂ = IOBuffer(), IOBuffer()
    i, j = size(score.M)
    while i > 1 && j > 1
        if score.M[i,j] == score.D[i,j]
            k = 1
            while j > k && ((score.M[i,j-k] + cost.gap(k)) != score.M[i,j])
                k += 1
            end
            write(a₁,'-'^k)
            write(a₂,seq₂[j-1:-1:j-k])
            j -= k
        elseif score.M[i,j] == score.I[i,j]
            k = 1
            while i > k && ((score.M[i-k,j] + cost.gap(k)) != score.M[i,j])
                k += 1
            end
            write(a₁,seq₁[i-1:-1:i-k])
            write(a₂,'-'^k)
            i -= k
        else
            write(a₁,seq₁[i-1])
            write(a₂,seq₂[j-1])
            i -= 1
            j -= 1
        end
    end

    if i > 1
        write(a₁,seq₁[i-1:-1:1])
        write(a₂,'-'^(i-1))
    elseif j > 1
        write(a₂,seq₂[j-1:-1:1])
        write(a₁,'-'^(j-1))
    elseif j > 1 && i > 1
        error("invalid backtraversal")
    end

    b₁ = reverse(take!(a₁))
    b₂ = reverse(take!(a₂))
    return flip ? (b₂, b₁) : (b₁, b₂)
end

"""
	hamming_align(qry::Array{UInt8,1}, ref::Array{UInt8,1})

Perform a simple alignment of `qry` to `ref` by minimizing hamming distance.
Useful for fast, approximate alignments of small sequences.
"""
function hamming_align(qry::Array{UInt8,1}, ref::Array{UInt8,1})
    matches = [
        let
            len = min(length(ref)-x+1, length(qry))
            sum(qry[1:len] .== ref[x:x+len-1])
        end
    for x ∈ 1:length(ref) ]

    return argmax(matches)
end

include("static/watson-crick.jl")
# XXX: This is the major allocator for us.
#      Will require serious thought as to how to not create so much garbage
"""
	reverse_complement(seq::Array{UInt8})

Return a newly allocated sequence array that is the reverse complement of `seq`.
"""
reverse_complement(seq::Array{UInt8}) = UInt8[wcpair[nuc+1] for nuc in reverse(seq)]

"""
	reverse_complement!(hit::Hit)

Reverse complement the qry of Hit in place.
"""
function reverse_complement!(hit::Hit)
    if hit.seq !== nothing
        hit.seq = reverse_complement(hit.seq)
    end

    start, stop = hit.start, hit.stop

    hit.start = 1 + (hit.length - stop)
    hit.stop  = hit.length - (start - 1)
end

function reverse_complement!(aln::Alignment)
    reverse_complement!(aln.qry)
    aln.orientation = ~aln.orientation
end

# ------------------------------------------------------------------------
# alignment modification

"""
	cost = (
		open   = -6.0,
		extend = -1.0,
		band   = (
			lower = Inf,
			upper = Inf,
		),
		gap    = k -> k == 0 ? 0 : cost.open + cost.extend*(k-1),
		match  = (c₁, c₂) -> 6.0*(c₁ == c₂) - 3.0,
	)

`cost` are the default dynamic alignment parameters used.
"""
const cost = (
    open   = -6.0,
    extend = -1.0,
    band   = (
        lower = Inf,
        upper = Inf,
    ),
    gap    = k -> k == 0 ? 0 : cost.open + cost.extend*(k-1),
    match  = (c₁, c₂) -> 6.0*(c₁ == c₂) - 3.0,
)
# TODO: come up with a better function name
"""
	enforce_cutoff!(a::Alignment, χ)

Ensure that the alignment `a` does not have insertion or deletion segments larger than `χ`.
Return the list of segments created by parsing the alignment such that all segments are larger than `χ`.
"""
function enforce_cutoff!(a::Alignment, χ)
    δqₗ, δqᵣ = a.qry.start-1, a.qry.length - a.qry.stop
    δrₗ, δrᵣ = a.ref.start-1, a.ref.length - a.ref.stop

    s₁ = a.orientation ? a.qry.seq : reverse_complement(a.qry.seq)
    s₂ = a.ref.seq

    # left side of match
    if (0 < δqₗ ≤ χ) && (δrₗ == 0 || δrₗ > χ)
        a.qry.start = 1
        pushfirst!(a.cigar, (δqₗ, 'I'))
    elseif (0 < δrₗ ≤ χ) && (δqₗ == 0 || δqₗ > χ)
        a.ref.start = 1
        pushfirst!(a.cigar, (δrₗ, 'D'))
    elseif (0 < δqₗ ≤ χ) && (0 < δrₗ <= χ)
        a₁, a₂ = align(s₁[1:δqₗ], s₂[1:δrₗ], cost)
        cg     = collect(uncigar(cigar(a₁, a₂)))

        a.qry.start = 1
        a.ref.start = 1

        prepend!(a.cigar, cg)
        a.length += length(a₁)
    end

    # right side of match
    if (0 < δqᵣ ≤ χ) && (δrᵣ == 0 || δrᵣ > χ)
        a.qry.stop  = a.qry.length
        push!(a.cigar, (δqᵣ, 'I'))
    elseif (0 < δrᵣ ≤ χ) && (δqᵣ == 0 || δqᵣ > χ)
        a.ref.stop  = a.ref.length
        push!(a.cigar, (δrᵣ, 'D'))
    elseif (0 < δqᵣ ≤ χ) && (δrᵣ <= χ)
        a₁, a₂ = align(s₁[end-δqᵣ+1:end], s₂[end-δrᵣ+1:end], cost)
        cg     = collect(uncigar(cigar(a₁, a₂)))

        a.qry.stop = a.qry.length
        a.ref.stop = a.ref.length

        append!(a.cigar, cg)
        a.length += length(a₁)
    end
end

# ------------------------------------------------------------------------
# insertion -> intervals

# ------------------------------------------------------------------------
# io functions

log(msg...) = println(stderr, msg...)

# fasta sequence record
"""
	struct Record
		seq::Array{UInt8}
		name::String
		meta::String
	end

A record obtained when parsing a single entry of a FASTA file.
"""
struct Record
    seq  :: Array{UInt8}
    name :: String
    meta :: String
end

name(r::Record) = isempty(r.meta) ? r.name : r.name * " " * r.meta

NL = '\n'
Base.show(io::IO, rec::Record) = print(io, ">$(rec.name) $(rec.meta)$(NL)$(String(rec.seq[1:40]))...$(String(rec.seq[end-40:end]))")

"""
	read_fasta(io::IO)

Parse a FASTA file from IO stream `io`.
Return an iterator over all records.
"""
function read_fasta(io::IO)
    chan = Channel{Record}(0)
    @async begin
        buf=IOBuffer()
        line=readline(io)
        while !isempty(line) && line[1] == '>'
            words      = split(line[2:end])
            name, meta = words[1], join(words[2:end], " ")

            line=readline(io)

            while !isempty(line) && line[1] != '>'
                write(buf,rstrip(line))
                line=readline(io)
            end
            put!(chan, Record(take!(buf), name, meta))
        end

        close(buf)
        close(chan)
    end

    return chan
end

function Base.show(io::IO, h::Hit)
    print(io, "$(h.name)[$(h.length)]: ($(h.start),$(h.stop))")
end

function Base.show(io::IO, a::Alignment)
    print(io, "qry: $(a.qry)", '\t')
    print(io, "ref: $(a.ref)", '\t')
    print(io, "polarity: $(a.orientation)")
end

"""
	read_paf(io::IO)

Parse a PAF file from IO stream `io`.
Return an iterator over all pairwise alignments.
"""
function read_paf(io::IO)
    chan = Channel{Alignment}(0)

    int(x)   = parse(Int,x)
    float(x) = parse(Float64,x)
    last(x)  = split(x,':')[end]

    @async begin
        for row in eachline(io)
            elt = split(strip(row))
            cg = nothing
            dv = nothing
            as = nothing
            for x in elt[13:end]
                if startswith(x, "cg:")
                    cg = last(x)
                elseif startswith(x, "de:f")
                    dv = float(last(x))
                elseif startswith(x, "AS:i")
                    as = int(last(x))
                end
            end

            # XXX: important: shift PAF entry to be 1 indexed and right-inclusive to match julia
            put!(chan, Alignment(
                        Hit(
                            elt[1],
                            int(elt[2]),   # length
                            int(elt[3])+1, # start
                            int(elt[4]),   # stop
                            nothing,
                        ),
                        Hit(
                            elt[6],
                            int(elt[7]),   # length
                            int(elt[8])+1, # start
                            int(elt[9]),   # stop
                            nothing,
                        ),
                        int(elt[10]),
                        int(elt[11]),
                        int(elt[12]),
                        elt[5] == "+",
                        cg,
                        dv,
                        as)
            )
        end
        close(chan)
    end

    return chan
end

# ------------------------------------------------------------------------
# string processing functions

"""
	columns(s; nc=80)

Partition string `s` into an array of strings such that no string is longer than `nc` characters.
"""
function columns(s; nc=80)
    nr   = ceil(Int64, length(s)/nc)
    l(i) = 1+(nc*(i-1))
    r(i) = min(nc*i, length(s))
    rows = [String(s[l(i):r(i)]) for i in 1:nr]
    return join(rows,'\n')
end

function test()
    # println(">testing fasta parse...")
    # open("data/test.fna") do io
    #     for record in read_fasta(io)
    #         println(String(record.seq))
    #         println(record.seq)
    #     end
    # end
    # println(">done!")

    # println(">testing paf parse")
    # open("data/test.paf") do io
    #     for aln in read_paf(io)
    #         println(aln)
    #     end
    # end
    # println("done!")

    println(">testing cigar serialization...")
    s₁ = Vector{UInt8}("A-TCGT-GTCA-TAGC")
    s₂ = Vector{UInt8}("AGG-GTCGTCAGT-GC")
    cg = cigar(s₁, s₂)
    println("-->", cg)
    println(">done!")

    cost = (
        open   = -0.75,
        extend = -0.5,
        band   = (
            lower = Inf,
            upper = Inf,
        ),
        gap    = k -> k == 0 ? 0 : cost.open + cost.extend*(k-1),
        match  = (c₁, c₂) -> 2.0*(c₁ == c₂) - 1.0,
    )

    seq = (N) -> Vector{UInt8}(random_id(;len=N, alphabet=['A','C','G','T']))

    s = [ (seq(100), seq(100)) for i in 1:100 ]
    # println("1: ", String(a₁))
    # println("2: ", String(a₂))
end

function contiguous_trues(x)
    intervals = Interval{Int}[]

    l      = 1
    inside = false

    for (r,σ) ∈ enumerate(x)
        if σ
            if !inside
                l = r
                inside = true
            end

            continue
        end

        if inside
            inside = false
            push!(intervals, Interval(l, r))
        end
    end

    # check for case that extends to right boundary
    if inside
        push!(intervals, Interval(l, length(x)+1))
    end

    return IntervalSet(intervals)
end

# ------------------------------------------------------------------------
# alignment processing functions

make_consensus(alignment) = [mode(view(alignment,i,:)) for i in 1:size(alignment,1)]

function alignment_alleles(ref, aln, nodes)
    isdiff = aln .!= ref
    refdel = ref .== UInt8('-')
    alndel = aln .== UInt8('-')

    δ = (
        snp =   isdiff .& .~refdel .& .~alndel,
        del = .~refdel .&   alndel,
        ins =   refdel .& .~alndel,
    )

    coord   = cumsum(.!refdel)
    refgaps = contiguous_trues(refdel)
    gaps    = Dict{Int,Int}(coord[gap.lo] => length(gap) for gap in refgaps)

    mutate = Dict{Node,SNPMap}(
            node => SNPMap(
                   coord[l] => aln[l,i]
                for l in findall(δ.snp[:,i])
            )
        for (i,node) in enumerate(nodes)
    )

    delete = Dict{Node,DelMap}(
            node => DelMap(
                      del.lo => length(del)
                for del in contiguous_trues(δ.del[.~refdel,i])
             )
        for (i,node) in enumerate(nodes)
    )

    Δ(I) = (R = containing(refgaps, I)) == nothing ? 0 : I.lo - R.lo
    insert = Dict{Node,InsMap}(
            node => InsMap(
                      (coord[ins.lo],Δ(ins)) => aln[ins,i]
                for ins in contiguous_trues(δ.ins[:,i])
             )
        for (i,node) in enumerate(nodes)
    )

    sequence = ref[.~refdel]

    return (
        gaps     = gaps,
        mutate   = mutate,
        delete   = delete,
        insert   = insert,
        sequence = sequence,
    )
end

end

max_stars_repo_name: neherlab/pangraph

max_issues_repo_name: neherlab/pangraph

max_forks_repo_name: neherlab/pangraph


text:
function oscillatory_basis_exp(omega, n)
    x = IdentityFunction()
    w1 = exp(im*omega*x)
    w2 = exp(-im*omega*x)
    if iseven(n)
        w1 * ChebyshevT(n>>1) ⊕ w2 * ChebyshevT(n>>1)
    else
        w1 * ChebyshevT(n>>1+1) ⊕ w2 * ChebyshevT(n>>1)
    end
end

function oscillatory_basis_cos(omega, n)
    x = IdentityFunction()
    w1 = cos(omega*x)
    w2 = sin(omega*x)
    if iseven(n)
        w1 * ChebyshevT(n>>1) ⊕ w2 * ChebyshevT(n>>1)
    else
        w1 * ChebyshevT(n>>1+1) ⊕ w2 * ChebyshevT(n>>1)
    end
end

function log_basis(n)
    b1 = Chebyhev(n>>1) → 0..1
    b2 = Chebyhev(n>>1+1) → 0..1
    iseven(n) ? b1 ⊕ Log() * b1 : b2 ⊕ Log() * b1
end

function exp_basis(n)
    b1 = Chebyhev(n>>1) → 0..1
    b2 = Chebyhev(n>>1+1) → 0..1
    iseven(n) ? b1 ⊕ Exp() * b1 : b2 ⊕ Exp() * b1
end


text:
module NeXLMatrixCorrection

using Reexport
using Requires

@reexport using NeXLCore  # Re-export so NeXLCore automatically is available with being explicitly loaded

# Abstract class for a base matrix correction algorithm (like XPP or ???)
include("matrixcorrection.jl")
export MatrixCorrection # Abstract class for the algorithm implementing the Z- & A-terms
export Z # Atomic number correction
export A # Absorption correction
export ZA # Combined ZA (ϕ(ρz)) correction (more efficient than Z(...)*A(...))
export ϕ # ϕ(ρz) function
export ϕabs # ϕ(ρz) function (absorbed)
export F  # Generated intensity function
export Fχ # Absorbed intensity function
export matrixcorrection # Factory method for MatrixCorrection algorithms (XPP and NullCorrection)
export continuumcorrection # Factory method for creating matrix correction algorithms for continuum correction
export correctcontinuum # The method to calculate the continuum correction
export χ # Calculates reduced mass-absorption coefficient

# Abstract class for a base fluorescence correction
include("fluorescencecorrection.jl")
export FluorescenceCorrection # Abstract class for the algorithm implementing the F-term
export NullFluorescence
# export F # Compute the flurorescence correction
export fluorescencecorrection # Factory method for FluorescenceCorrection algorithms (Reed and NullFluorescence)

# Abstract class for a base coating correction algorithm
include("coating.jl")
export CoatingCorrection
export NullCoating  # 100% transmission (no correction)
export Coating      # A basic multi-layer coating for ultra-thin coatings
export carboncoating # Build a carbon coating
export coatingcorrection # Factory method for CoatingCorrection algorithms (Null or Coating)

# Pulls together the MatrixCorrection, FluorescenceCorrection and CoatingCorrection into a single object
include("zafcorrection.jl")
export ZAFCorrection
export zafcorrection # Factory method for ZAFCorrection
export ZAFc # Combined correction factor ZAF + coating
export beamEnergy

# Use massthickness to estimate the thickness of an ultra-thin coating on a substrate.
export massthickness # Estimate the mass-thickness of a coating on a substrate
export coatingasfilm # Calculates the coating thickness and creates a Film object to represent it.
export kcoating # Calculate the measured k for an ultra-thin coating on a substrate

# Wraps multiple ZAFCorrection objects to calculate the full matrix correction for one or more characteristic x-rays
# This addresses the need to
include("multizaf.jl")
export MultiZAF # Represents a multiline ZAF correction
export gZAFc # Combined correction factor (generation + ZAF + coating)
export k # Computed k-ratio
export aspure # The k-ratio for this X-ray relative to a pure element

# Implements Pouchou & Pichoir's XPP  ϕ(ρz) model
include("xpp.jl")
export XPP # <: MatrixCorrection
export xpp # Build an XPP correction

include("citzaf.jl")
export CitZAF # <: MatrixCorrection

include("riveros.jl")
export Riveros1993 # <: MatrixCorrection

# Implements Pouchou & Pichoir's XPP  ϕ(ρz) model with full uncertainty calculation
include("xppu.jl")
export JzLabel
export Ju
export StepMJZbarb, StepDPT, StepQlaOoS

# Implements Merlet's XPhi algorithm
include("xphi.jl")
export XPhi

# Implements the null matrix correction
include("nullmatrixcorrection.jl")
export NullCorrection

# Implements Reed's 1991 fluorescence correction
include("reed.jl")
export ReedFluorescence
export reedFluorescence # Construct a structure encapsulating the Reed fluorescence correction model

# For picking which characteristic lines to use to perform the matrix correction
include("kratioopt.jl")
export KRatioOptimizer # Abstract class
export SimpleKRatioOptimizer # A very simple implmentation of KRatioOptimizer
export optimizeks # The method required of KRatioOptimizer

include("unmeasuredelement.jl")
# Abstract class for calculating elements by difference or stoichiometry or ???
export UnmeasuredElementRule
export isunmeasured # Does this UnmeasuredElementRule handle this element?
export NullUnmeasuredRule # Don't do anything..
export ElementByStoichiometry # Any element by stoichiometry
export OByStoichiometry # O by stoichiometry
export ElementByFiat # Just say it..
export MultiUnmeasuredElementRule # Combine rules sequentially

# Performs iteration to estimate the composition from measured k-ratios
include("iterate.jl")
# Abstract class for updating between iteration steps
export UpdateRule # A rule implementing update(...)
export NaiveUpdateRule # Simple iteration
export WegsteinUpdateRule # Wegstein iteration (gradient)
export RecordingUpdateRule

export ConvergenceTest # Test for convergence using converged(...)
export RMSBelowTolerance, AllBelowTolerance, IsApproximate # Difference implementations of ConvergenceTest

export Iteration # Defines the iteration procedure
export IterationResult # The output from quantify(...)

export quantify # Perform the iteration on KRatio(s) or FilterFitResult
export estimatecoating # Use a measured k-ratio to estimate the coating thickness on a material

include("aspure.jl")
#export aspure

include("supportedthinfilms.jl")
export SupportedThinFilms
export nlayers
export χs
export outer, inner

#include("xfilm.jl")
#export XFilm

include("helpers.jl")
export zaf # Generates a table of ZAF correction factors

include("defaultstandards.jl")
export getstandards # Get a list of suggested standards for an element.

function __init__()
    @require Gadfly = "c91e804a-d5a3-530f-b6f0-dfbca275c004" include("gadflysupport.jl")
end

end

max_stars_repo_name: NicholasWMRitchie/NeXLMatrixCorrection

max_issues_repo_name: NicholasWMRitchie/NeXLMatrixCorrection

max_forks_repo_name: NicholasWMRitchie/NeXLMatrixCorrection


text:
using JSON

# Return the long image SHA256 identifier using various ways of referencing images
function full_image_sha(image::AbstractString)
    json = JSON.parse(read(`docker inspect $image`, String))
    return last(split(json[1]["Id"], ':'))
end

const DOCKER_SPAWN_REGEX = r"^Spawning container: [0-9a-z]{12}$"

@testset "DockerManager" begin
    # A docker image name which is expected to not exist on the local system.
    mock_image = "x5cn2a7hzq4k"

    @testset "constructors" begin
        @testset "inner" begin
            mgr = DockerManager(2, mock_image, Minute(10), ip"1.0.0.0", ip"2.0.0.0")

            # Validate that no additional fields were added without the tests being updated
            @test fieldcount(DockerManager) == 5

            @test mgr.num_workers == 2
            @test mgr.image == mock_image
            @test mgr.timeout == Minute(10)
            @test mgr.min_ip == ip"1.0.0.0"
            @test mgr.max_ip == ip"2.0.0.0"

            @test launch_timeout(mgr) == Minute(10)
            @test desired_workers(mgr) == (2, 2)
        end

        @testset "keywords" begin
            mgr = DockerManager(
                3;
                image=mock_image,
                timeout=Minute(12),
                min_ip=ip"3.0.0.0",
                max_ip=ip"4.0.0.0",
            )

            @test mgr.num_workers == 3
            @test mgr.image == mock_image
            @test mgr.timeout == Minute(12)
            @test mgr.min_ip == ip"3.0.0.0"
            @test mgr.max_ip == ip"4.0.0.0"
        end

        @testset "defaults" begin
            patch = @patch AWSClusterManagers.image_id() = mock_image

            apply(patch) do
                mgr = DockerManager(0)

                @test mgr.num_workers == 0
                @test mgr.image == mock_image
                @test mgr.timeout == AWSClusterManagers.DOCKER_TIMEOUT
                @test mgr.min_ip == ip"0.0.0.0"
                @test mgr.max_ip == ip"255.255.255.255"
            end
        end

        @testset "num workers" begin
            # Define keywords which are required to avoid mocking
            kwargs = Dict(:image => mock_image)

            @test_throws ArgumentError DockerManager(-1; kwargs...)
            @test desired_workers(DockerManager(0; kwargs...)) == (0, 0)
            @test desired_workers(DockerManager(2; kwargs...)) == (2, 2)

            # Note: DockerManager does not support ranges
            @test_throws MethodError DockerManager(1, 2; kwargs...)
            @test_throws MethodError DockerManager(3:4; kwargs...)
        end
    end

    @testset "equality" begin
        patch = @patch AWSClusterManagers.image_id() = mock_image

        apply(patch) do
            @test DockerManager(3) == DockerManager(3)
        end
    end

    @testset "addprocs" begin
        @testset "success" begin
            patch = @patch function read(cmd::AbstractCmd, ::Type{String})
                # Extract original `override_command` using image position:
                # "docker run [OPTIONS] IMAGE [COMMAND] [ARG...]"
                i = findfirst(isequal(mock_image), collect(cmd))
                override_cmd = Cmd(cmd[(i + 1):end])
                @async run(override_cmd)
                return "000000000001"
            end

            apply(patch) do
                # Get an initial list of processes
                init_procs = procs()
                # Add a single AWSBatchManager worker
                added_procs = @test_log LOGGER "notice" DOCKER_SPAWN_REGEX begin
                    addprocs(DockerManager(1, mock_image))
                end
                # Check that the workers are available
                @test length(added_procs) == 1
                @test procs() == vcat(init_procs, added_procs)
                # Remove the added workers
                rmprocs(added_procs; waitfor=5.0)
                # Double check that rmprocs worked
                @test init_procs == procs()
            end
        end

        @testset "worker timeout" begin
            patch = @patch function read(cmd::AbstractCmd, ::Type{String})
                # Avoiding spawning a worker process
                return "000000000002"
            end

            @test_throws TaskFailedException apply(patch) do
                @test_log LOGGER "notice" DOCKER_SPAWN_REGEX begin
                    addprocs(DockerManager(1, mock_image, Second(1)))
                end
            end
        end
    end
end


text:
import Distributed

@inline function Distributed.launch(manager::ExistingProcessManager,
                                    params::Dict,
                                    launched::Array,
                                    launch_ntfy::Condition)
    while !isempty(manager.wconfigs)
        wconfig = pop!(manager.wconfigs)
        push!(launched, wconfig)
        notify(launch_ntfy)
    end
    return nothing
end


text:
module TestPlots
using CustomTest
using Base.Test
using ControlSystems, Plots
using VisualRegressionTests, ControlExamplePlots

default(show=false)

funcs, refs, eps = getexamples()
# Make it easier to pass tests on different systems
# Set to a factor 2 of common errors
eps = [0.15, 0.015, 0.1, 0.01, 0.01, 0.02, 0.01, 0.15, 0.15, 0.01]
res = genplots(funcs, refs, eps=eps, popup=false)

##Explicit enumeration for simpler debugging

#"bode.png"
@test  res[1] |> success
#"nyquist.png"
@test  res[2] |> success
#"sigma.png"
@test  res[3] |> success
#"nichols.png"
@test  res[4] |> success
#"step.png"
@test  res[5] |> success
#"impulse.png"
@test  res[6] |> success
#"lsim.png"
@test  res[7] |> success
#"margin.png"
@test  res[8] |> success
#"gangoffour.png"
@test  res[9] |> success
#"pzmap.png"
@test  res[10] |> success

end


text:
module ICapeThermoMaterial

    import CapeOpen: MaterialObject
    import ..ECapeExceptions: ECapeNoImpl, ECapeUnknown
    export clearallprops!,copyfrommaterial!,creatematerial!,getoverallprop!
    export getoveralltpfraction!,getpresentphases!,getsinglephaseprop!
    export gettpfraction!,gettwophaseprop!,setoverallprop,setpresentphases
    export setsinglephaseprop,settwophaseprop

    """
      Remove all stored Physical Property values.
    """
    function  clearallprops!(
        #= [in] =# this::MaterialObject)

		end

    """
      Creates a Material Object with the s ration as the current Material Object
    """
    function  copyfrommaterial!(
        #= [in] =# this::MaterialObject,
        #= [in] =# source::MaterialObject)
		end

    """
      Copies all the stored non-constant Physical  Properties (which have been set using the
      SetSinglePhaseProp SetTwoPhaseProp or SetOverallProp) from the  source Material Object to the
      current instance Material Object.
    """
    function  creatematerial!(
        #= [retval][out] =# materialObject::MaterialObject)
		end

    """
      Retrieves non- nt Physical Property values for the overall mixture.
    """
    function  getoverallprop!(
        #= [in] =# this::MaterialObject,
        #= [in] =# property::AbstractString,
        #= [in] =# basis::AbstractString,
        #= [out][in] =# results::Vector{Float64})
		end

    """
      Retrieves tempe pressure an on for the overall mixture.
    """
    function  getoveralltpfraction!(
        #= [in] =# this::MaterialObject,
        #= [out][in] =# temperature!::Float64,
        #= [out][in] =# pressure!::Float64,
        #= [out][in] =# composition!::Vector{Float64})
		end

    """
      Returns Phase l for the Phaes that are currently present in the Material Object.
    """
    function  getpresentphases!(
        #= [in] =# this::MaterialObject,
        #= [out][in] =# phaseLabels::Vector{AbstractString},
        #= [out][in] =# phaseStatus!::Vector{Int})
		end

    """
      Retrieves single non-constan Physical Property values for a mixture.
    """
    function  getsinglephaseprop!(
        #= [in] =# this::MaterialObject,
        #= [in] =# property::AbstractString,
        #= [in] =# phaseLabel::AbstractString,
        #= [in] =# basis::AbstractString,
        #= [out][in] =# results::Vector{Float64})
		end

    """
      Retrieves temperature, pressure and composition for a Phase.
    """
    function  gettpfraction!(
        #= [in] =# this::MaterialObject,
        #= [in] =# phaseLabel::AbstractString,
        #= [out][in] =# temperature!::Float64,
        #= [out][in] =# pressure!::Float64,
        #= [out][in] =# composition!::Vector{Float64})
		end

    """
      Retrieves two-p on-constant P perty values for a mixture.
    """
    function  gettwophaseprop!(
        #= [in] =# this::MaterialObject,
        #= [in] =# property::AbstractString,
        #= [in] =# phaseLabels::Vector{AbstractString},
        #= [in] =# basis::AbstractString,
        #= [out][in] =# results::Vector{Float64})
		end

    """
      Sets non-constant property values f ll mixture.
    """
    function  setoverallprop(
        #= [in] =# this::MaterialObject,
        #= [in] =# property::AbstractString,
        #= [in] =# basis::AbstractString,
        #= [in] =# values::Vector{Float64})
		end

    """
      Allows the PME or the Property Package to specify the list of Phases that are currently present.
    """
    function  setpresentphases(
        #= [in] =# this::MaterialObject,
        #= [in] =# phaseLabels::Vector{AbstractString} ,
        #= [in] =# phaseStatus::Vector{Int})
		end

    """
      Sets single-phas -constant pro for a mixture.
    """
    function  setsinglephaseprop(
        #= [in] =# this::MaterialObject,
        #= [in] =# property::AbstractString,
        #= [in] =# phaseLabel::AbstractString,
        #= [in] =# basis::AbstractString,
        #= [in] =# values::Vector{Float64})
		end

    """
      Sets two-phase nstant property values for a mixture.
    """
    function  settwophaseprop(
        #= [in] =# this::MaterialObject,
        #= [in] =# property::AbstractString,
        #= [in] =# phaseLabels::Vector{AbstractString},
        #= [in] =# basis::AbstractString,
        #= [in] =# values::Vector{Float64})
		end

end


text:
"""
    function get_metadata_urls(
        time_periods::Vector{Tuple{Int64, Int64}}
    )::Vector{String}

Creates an array of URLs for the metadata files

Parameters
* time_periods: Vector of time periods (year, quarter) to get metadata files Vector{Tuple{year, quarter}}

Returns: Vector of metadata urls
"""
function get_metadata_urls(time_periods::Vector{Tuple{Int64, Int64}})::Vector{String}

    # get urls for all time get_time_periods
    edgar_prefix = "https://www.sec.gov/Archives/"
    urls = [edgar_prefix * "edgar/full-index/"*string(i[1])*"/QTR"*string(i[2])*"/master.zip" for i in time_periods]

    return urls
end

"""
    function download_metadata(url::String; dest::String,
        skip_file=false::Bool, verbose=false::Bool
    )

Download filing metadata CSV file

Parameters
* url: URL where metadata file is hosted
* dest: Destination folder
* skip_file: If true, existing files will be skipped
* verbose: Print out log
"""
function download_metadata(url::String; dest::String, skip_file=false::Bool, verbose=false::Bool)

    # get full file path for download
    full_file = split(url, "/")[end-2] * "-" * split(url, "/")[end-1] * ".tsv"
    full_file = joinpath(dest, full_file)
    if verbose; println(full_file); end

    #TODO: unique temp files, so we can async download metadata
    temp_file = "main.idx"
    temp_zip = "main.zip"

    # check if we skip the download
    if isfile(full_file) & skip_file
        if verbose; println("Skipping " * full_file); end
        return
    end

    # download, unzip file
    HTTP.download(url, temp_zip)
    zarchive = ZipFile.Reader(temp_zip)
    for f in zarchive.files
        @assert f.name == "master.idx"
        out = open(temp_file,"w")
        write(out,read(f,String))
        close(out)
    end
    close(zarchive)
    rm(temp_zip)

    # import unziped file for cleaning
    f = open(temp_file, "r")
    metadata = readlines(f)[10:end] # skip fluff at top
    close(f)
    rm(temp_file)

    # save metadata file
    f = open(full_file, "w")
    for line in metadata
        if occursin("|", line) # skip "----------" line
            write(f, line * "\n")
        end
    end
    close(f)

    return
end

"""
    function download_metadata_files(start_year::Int64, end_year=nothing::Union{Int64, Nothing};
        quarters=[1, 2, 3, 4]::Vector{Int64},
        skip_file=false::Bool,
        dest="../metadata/"::String,
        verbose=false::Bool,
        download_rate=10::Int
    )

Download all metadata files over a time range

Parameters
* start_year: first year in range
* end_year: last year in range
* quarters: Quarters of the year to download files from [1,2,3,4]
* skip_file: If true, existing files will be skipped
* verbose: Print out log
"""
function download_metadata_files(start_year::Int64, end_year=nothing::Union{Int64, Nothing};
                        quarters=[1, 2, 3, 4]::Vector{Int64},
                        skip_file=false::Bool,
                        dest="../metadata/"::String,
                        verbose=false::Bool)

    # verify download_rate is valid (less than 10 requests per second, more than 0)
    if download_rate > 10
        download_rate = 10
        println("download_rate of more than 10 per second(", download_rate, ") is not valid. download_rate has been set to 10/second.")
    elseif download_rate < 1
        download_rate = 1
        println("download_rate of less than 1 per second(", download_rate, ") is not valid. download_rate has been set to 1/second.")
    end

    # create download folder if needed
    println("Metadata Destination:  " * dest)
    if !isdir(dest)
        mkdir(dest)
    end

    # get current year, quarter to prevent errors trying to get future data
    current_date = Dates.now()
    current_year = Dates.year(current_date)
    current_quarter = Dates.quarterofyear(current_date)

    # set end year to current year if no year is specified
    if end_year == nothing
        end_year = current_year
    end

    # get an array of dates to download metadata
    years = collect(start_year:end_year)
    time_periods = [(y, q) for y in years for q in quarters if (q <= current_quarter || y < current_year) && (q > 2 || y > 1993)]

    # get download urls
    urls = get_metadata_urls(time_periods)

    # download metadata files
    @showprogress 1 "Downloading Metadata..."  for idx in eachindex(urls)
        download_metadata(urls[idx]; dest=dest,
                    skip_file=skip_file, verbose=verbose)
    end

    return
end


text:
# const PSDGaussian{T} = Gaussian{Vector{T}, PSDMatrix{T}}
# const PSDGaussianList{T} = StructArray{PSDGaussian{T}}
const SRGaussian{T,S,M} = Gaussian{Vector{T},SRMatrix{T,S,M}}
const SRGaussianList{T,S,M} = StructArray{SRGaussian{T,S,M}}

copy(P::Gaussian) = Gaussian(copy(P.μ), copy(P.Σ))
similar(P::Gaussian) = Gaussian(similar(P.μ), similar(P.Σ))
function Base.copy!(dst::Gaussian, src::Gaussian)
    copy!(dst.μ, src.μ)
    copy!(dst.Σ, src.Σ)
    return dst
end
RecursiveArrayTools.recursivecopy(P::Gaussian) = copy(P)
show(io::IO, g::Gaussian) = print(io, "Gaussian($(g.μ), $(g.Σ))")
show(io::IO, ::MIME"text/plain", g::Gaussian{T,S}) where {T,S} =
    print(io, "Gaussian{$T,$S}($(g.μ), $(g.Σ))")
size(g::Gaussian) = size(g.μ)
ndims(g::Gaussian) = ndims(g.μ)

Base.:*(M, g::SRGaussian) = Gaussian(M * g.μ, X_A_Xt(g.Σ, M))
# GaussianDistributions.whiten(Σ::SRMatrix, z) = Σ.L\z

function mul!(g_out::SRGaussian, M::AbstractMatrix, g_in::SRGaussian)
    mul!(g_out.μ, M, g_in.μ)
    X_A_Xt!(g_out.Σ, g_in.Σ, M)
    return g_out
end

var(p::SRGaussian{T}) where {T} = diag(p.Σ)
std(p::SRGaussian{T}) where {T} = sqrt.(var(p))
mean(s::SRGaussianList{T}) where {T} = s.μ
var(s::SRGaussianList{T}) where {T} = diag.(s.Σ)
std(s::SRGaussianList{T}) where {T} = map(std, s)


text:
export getUjProjMatrix
export getVectorMassMatrix, getDifferentialOperators,GetLinearElasticityOperator,GetLinearElasticityOperatorFullStrain

function ddxCN(n::Int64,h) ## From nodes to cells
	D = (1/h)*ddx(n);
	return D
end

function ddxNC(n::Int64,h) ## From cells to nodes
	#A = 1/h*spdiagm((-ones(n),ones(n)),[-1,0],n+1,n)
	I, J, V = SparseArrays.spdiagm_internal(-1 => fill(-1.0,n), 0 => fill(1.0,n));
	D = sparse(I, J, (1/h)*V,n+1,n);
	## This is for boundary conditions: need to figure this out
	D[1,1] = 1.0/h;
	D[end,end] = -1.0/h;
	return D
end

function getDifferentialOperators(M::RegularMesh,optype = 1)
	n = M.n;
	h = M.h;
	if length(n)==3
		# Face sizes
		nf1 = prod(n + [1; 0; 0])
		nf2 = prod(n + [0; 1; 0])
		nf3 = prod(n + [0; 0; 1])
		nf  = [nf1; nf2; nf3];

		# Notation Dij = derivative of component j in direction i
		tmp = ddxCN(n[1],h[1])
		D11 = kron(speye(n[3]),kron(speye(n[2]),tmp))

		tmp = ddxNC(n[1],h[1])
		D12 = kron(speye(n[3]),kron(speye(n[2]+1),tmp))

		tmp = ddxNC(n[1],h[1])
		D13 = kron(speye(n[3]+1),kron(speye(n[2]),tmp))

		tmp = ddxNC(n[2],h[2])
		D21 = kron(speye(n[3]),kron(tmp,speye(n[1]+1)))

		tmp = ddxCN(n[2],h[2])
		D22 = kron(speye(n[3]),kron(tmp,speye(n[1])))

		tmp = ddxNC(n[2],h[2])
		D23 = kron(speye(n[3]+1),kron(tmp,speye(n[1])))

		tmp = ddxNC(n[3],h[3])
		D31 = kron(tmp,kron(speye(n[2]),speye(n[1]+1)))

		tmp = ddxNC(n[3],h[3])
		D32 = kron(tmp,kron(speye(n[2]+1),speye(n[1])))

		tmp = ddxCN(n[3],h[3])
		D33 = kron(tmp,kron(speye(n[2]),speye(n[1])))

		# Tensor sizes
		t = [size(D11,1); size(D12,1); size(D13,1);
			size(D21,1); size(D22,1); size(D23,1);
			size(D31,1); size(D32,1); size(D33,1);]

		Div = [D11 D22 D33]
		if optype==1
			vectorDiv = [Div; spzeros(t[2]+t[3]+t[4],sum(nf)); Div; spzeros(t[6]+t[7]+t[8],sum(nf)); Div];
			vectorGrad  = 0.5*[blockdiag(D11,D12,D13);blockdiag(D21,D22,D23);blockdiag(D31,D32,D33)]
			vectorGradT = 0.5*blockdiag([D11;D21;D31],[D12;D22;D32],[D13;D23;D33])
			D11 = 0; D12 = 0; D13 = 0; D22 = 0; D21 = 0; D23 = 0; D31 = 0; D32 = 0; D33 = 0;
			vectorGrad  = (vectorGrad  + vectorGradT);
		else
			vectorGrad = blockdiag([D11;D21;D31],[D12;D22;D32],[D13;D23;D33]);
			vectorDiv = 0;
			vectorGradT = 0;
		end
	else
		nf1 = prod(n + [1; 0]);
		nf2 = prod(n + [0; 1]);
		nf  = [nf1; nf2];

		# Notation Dij = derivative of component j in direction i
		tmp = ddxCN(n[1],h[1]);
		D11 = kron(speye(n[2]),tmp);

		tmp = ddxNC(n[1],h[1]);
		D12 = kron(speye(n[2]+1),tmp);

		tmp = ddxNC(n[2],h[2])
		D21 = kron(tmp,speye(n[1]+1))

		tmp = ddxCN(n[2],h[2])
		D22 = kron(tmp,speye(n[1]))

		# Tensor sizes
		t = [size(D11,1); size(D12,1); size(D21,1); size(D22,1);]

		Div = [D11 D22]

		if optype==1
			vectorDiv = [Div; spzeros(t[2]+t[3],sum(nf)); Div; ];
			vectorGrad  = [blockdiag(D11,D12);blockdiag(D21,D22)]
			vectorGradT = blockdiag([D11;D21],[D12;D22])
			vectorGrad  = 0.5*(vectorGrad  + vectorGradT)
		else
			# vectorGrad = [[D11;D21],[D12;D22]];
			vectorGrad = blockdiag([D11;D21],[D12;D22]); ## diagGrad
			vectorDiv = 0;
			vectorGradT = 0;
		end
	end
	return vectorGrad, vectorDiv, Div,nf,vectorGradT
end

function getTensorMassMatrix(M::RegularMesh, mu::Vector; saveAvMat::Bool = true)

	if M.dim==3
		n = M.n;
		# (Ae1, Ae2, Ae3) = getEdgeAveragingMatrices(n)
		# mu    = vec(mu);
		# Ae3mu = Ae3*mu;
		# Ae1mu = Ae1*mu;
		# Ae2mu = Ae2*mu;

		Ae = getEdgeAverageMatrix(M; saveMat = saveAvMat, avN2C = avN2C_Nearest);
		AeMu = Ae'*mu;
		Ae1mu = AeMu[1:M.nf[1]];
		Ae2mu = AeMu[(M.nf[1]+1):(M.nf[1]+M.nf[2])];
		Ae3mu = AeMu[(M.nf[1]+M.nf[2]+1):end];
		m = [mu;Ae3mu;Ae2mu;Ae3mu;mu;Ae1mu;Ae2mu;Ae1mu;mu];
		Mass = spdiagm(m);
		len = [0;length(mu)+length(Ae3mu)+length(Ae2mu);length(Ae3mu)+length(mu)+length(Ae1mu);length(Ae2mu)+length(Ae1mu)+length(mu)];
		len[3] = len[3]+len[2];
		len[4] = len[4] + len[3];
	else
		# (Ae1, Ae2) = getEdgeAveragingMatrices(n)
		# mu    = vec(mu);
		# Ae1mu = Ae1*mu;
		# Ae2mu = Ae2*mu;
		# m = [mu;Ae1mu;Ae2mu;mu];
		# Mass = spdiagm(m);
		# len = [0;length(mu)+length(Ae1mu);length(Ae2mu)+length(mu)];
		# len[3] = len[3]+len[2];

		# In 2D the strain tensor (when u is on faces) is defined on the nodes.
		An = getNodalAverageMatrix(M; saveMat = saveAvMat, avN2C = avN2C_Nearest)
		AnMu = An'*mu;
		Mass = spdiagm([mu;AnMu;AnMu;mu]);
		len_nodes = (M.n[1]+1)*(M.n[2]+1);
		len = [0;length(mu)+len_nodes;len_nodes+length(mu)];
		len[3] = len[3]+len[2];
	end
	return Mass,len
end

function getUjProjMatrix(n,j)
	if length(n)==3
		nf1 = prod(n + [1; 0; 0])
		nf2 = prod(n + [0; 1; 0])
		nf3 = prod(n + [0; 0; 1])
		nf  = [0;nf1;nf2+nf1;nf3+nf1+nf2];
		Pj = speye(nf1+nf2+nf3);
		Pj = Pj[:,nf[j]+1:nf[j+1]];
	else
		nf1 = prod(n + [1; 0])
		nf2 = prod(n + [0; 1])
		nf  = [0;nf1;nf2+nf1];
		Pj = speye(nf1+nf2);
		Pj = Pj[:,nf[j]+1:nf[j+1]];
	end
	return Pj;
end

function GetLinearElasticityOperator(M::RegularMesh, mu::Vector,lambda::Vector)
	vecGRAD,~,Div,nf, = getDifferentialOperators(M,2);
	Mu     = getTensorMassMatrix(M,mu,saveAvMat=false)[1];
	LambdaMuCells = spdiagm(lambda[:] + mu[:]);
	Rho	   = getFaceMassMatrix(M,1e-5*mu, saveMat = false, avN2C = avN2C_Nearest);
	G = vecGRAD'*Mu*vecGRAD

	H = Div'*LambdaMuCells*Div + vecGRAD'*Mu*vecGRAD + Rho;
	return H;
end

function GetLinearElasticityOperatorFullStrain(M::RegularMesh, mu::Vector,lambda::Vector)
	vecGRAD,vecDiv,Div,nf, = getDifferentialOperators(M,1);
	Mu     = getTensorMassMatrix(M,mu,saveAvMat=false)[1];
	Lambda = getTensorMassMatrix(M,lambda,saveAvMat=false)[1];
	Rho	   = getFaceMassMatrix(M,1e-5*mu, saveMat = false, avN2C = avN2C_Nearest);
	H = vecGRAD'*(Lambda*vecDiv + 2.0*Mu*vecGRAD) + Rho;
	return H;
end


text:
### Saving results

function saveresults(number_zones, state_matrix, transition_matrix, path_to_results, data_set, C)

    # prepare the resulting location of cars, the share of parking and driving cars in each zone
    driving_cars = zeros(number_zones, T);
    parking_cars = zeros(number_zones, T);
    for i=1:C
        for t=1:T
            index = state_matrix[i,t];
            driving_activity = transition_matrix[i,t,1];
            parking_cars[index,t] += 1;
            if(driving_activity == 1)
                driving_cars[index,t] += 1;
            end
        end
    end

    # normalize with the total number of cars
    parking_cars = parking_cars/C;

    # calculating circadian rythm of driving activity
    traffic_resultsmatrix = sum(driving_cars,dims=1);
    min_sampled = minimum(traffic_resultsmatrix);
    max_sampled = maximum(traffic_resultsmatrix);
    for i=1:24
        traffic_resultsmatrix[i] = (traffic_resultsmatrix[i]-min_sampled)/(max_sampled - min_sampled);
    end

    # convert to DataFrame for saving as CSV files
    parking_cars = convert(DataFrame, parking_cars);
    traffic_resultsmatrix = convert(DataFrame, traffic_resultsmatrix);

    # create the header for the CSV file and save the results
    header_vector = fill("", T); # creates an empty string vector in the desired dimension
    for t = 1:T
        header_vector[t] = string("t = ", t, "h");
    end
    results_string_parking = string(path_to_results,"/results_parkingdensities_",data_set);
    results_string_driving = string(path_to_results,"/results_trafficactivity_",data_set);
    CSV.write(results_string_parking, parking_cars, header=header_vector);
    CSV.write(results_string_driving, traffic_resultsmatrix, header=header_vector);

end

max_stars_repo_name: ArsamAryandoust/CarParkingMaps

max_issues_repo_name: ArsamAryandoust/CarParkingMaps

max_forks_repo_name: ArsamAryandoust/CarParkingMaps


text:
"""
$(TYPEDEF)
"""
struct LinearSolution{T,N,uType,R,A,C} <: AbstractLinearSolution{T,N}
  u::uType
  resid::R
  alg::A
  retcode::Symbol
  iters::Int
  cache::C
end

function build_linear_solution(alg,u,resid,cache;
                               retcode = :Default,
                               iters = 0)
  T = eltype(eltype(u))
  N = length((size(u)...,))
  LinearSolution{T,N,typeof(u),typeof(resid),typeof(alg),typeof(cache)}(u,resid,alg,retcode,iters,cache)
end

struct QuadratureSolution{T,N,uType,R,P,A,C} <: AbstractQuadratureSolution{T,N}
  u::uType
  resid::R
  prob::P
  alg::A
  retcode::Symbol
  chi::C
end

function build_solution(prob::AbstractQuadratureProblem,
                        alg,u,resid;calculate_error = true,
                        chi = nothing,
                        retcode = :Default, kwargs...)

  T = eltype(eltype(u))
  N = length((size(u)...,))

  QuadratureSolution{T,N,typeof(u),typeof(resid),
                     typeof(prob),typeof(alg),typeof(chi)}(
                     u,resid,prob,alg,retcode,chi)
end


text:
# single line comment
#= multi
   line
   comments
   https://learnxinyminutes.com/docs/julia/
=#

# -- types and operators --

3
3.2
2+1im
2//3 # rational

div(5,2) # 2
2^2 # 4
6%5 # 1

true
false

s = "string"
println(s[1]) # s, julia indexes from 1
println("expression in string $(s[3:4])")

@printf "%s\n" "hello"
println("println()")
println(bits(Int64(12345)))

λ = 1 # unicode characters are recognized

# -- variables and collections --

a = Int64[] # empty array
b = [4, 5, 6]
c = [4; 5; 6]
println(b[1], b[end])

mata = [1 2; 3 4] # 2x2 array
push!(a, 1) # add element to last position
push!(a, 2)
append!(a, b) # append b elements to a
pop!(b)
push!(b, 6)
shift!(a) # all elements move left and pop the first element
unshift!(a, 7) # all elements go right and put the element to head

arr = [5,4,6]
sort(arr) # will not change arr
sort!(arr) # arr will be changed

some_var = 5
try
	var_not_defined
	a[0]
	a[end+1]
catch e
	println(e)
end

a = [1:5;] # ~ Array(1:5)
a[1:3] # slicing
a[2:end] # slicing

in(1, [1,2,3]) # true
1 in [1,2,3] # true
length(a)

tup = (1,2,3) # tuples are immutable
x,y,z = tup # unpack tuples
(1,) != (1) # one-element tuple

d = Dict("one"=> 1, "two"=> 2)
d["one"]
keys(d)
values(d)
in(("one"=> 1), d) # true
in(("two"=> 3), d) # false
haskey(d, "one")
# d["four"] # ERROR
get(d, "one", 4) # 1
get(d, "four", 4) # 4

Set()
seta = Set([1,1,2,2,3,3]) # Set supports push!(), in() and so on ...
println(seta)
# intersect() union() setdiff()

# -- control flow --
if 1 < 2
	println("msg")
elseif 3 < 4
	println("ridiculous")
else
	println("else")
end

for fruit in ["apple", "peach", "banana"]
	println("$fruit is a kind of fruit")
end
for mapping in Dict("one"=> 1, "two"=> 2)
# for (k,v) in Dict("one"=> 1, "two"=> 2) # also works
	println("$(mapping[1]) => $(mapping[2])")
end

j = 1
while j < 4
	j+=j
end

try
	error("example error")
catch e
	println("caught $e")
end

# -- functions --
function add(x, y)
	println("adding $x and $y")
	x+y
end
add(5, 6)
f_add(x, y) = x+y # compact assignment of functions
f_foo(x, y) = x+y, x-y # returns a tuple

function va_args(args...)
	return args
end
va_args(1, 2, 3)

add([5,6]...) # equiv to add(5,6), the ... is called a splat.

function f_positional(a, b, x=5, y=6) # positional arguments
	return "$a $b and $x $y"
end
function f_keyword(; k1=4, name2="hello") # keyword argument
	return Dict("k1"=> k1, "name2"=> name2)
end

function create_adder(x) # Julia has first class functions
	adder = function (y)
		return x+y
	end
	return adder
end

println( (x -> x>2)(3) ) # true; this is stabby lambda syntax

map(+, [1,2,3,4]) # built-in higher order functions
filter(x -> x>5, [3,4,5,6,7])

println([(x -> x+1)(i) for i in [1, 2, 3]]) # list comprehensions

# -- types --
typeof(5) # Int64
type MyType
	real::Float64
	imag::Float64
	tag::AbstractString
	stub # same as stub::Any
end
aaa = MyType(0.1, 0.2, "test", 0)
bbb = typeof(aaa)(0, 0, "bbb", 0)

abstract Cat # just a name and node in the type hierarchy
println(subtypes(AbstractString))
println(super(Int64))

type Lion <: Cat # <: is the subtype operator
	name_color
	roar::AbstractString
end
Lion(roar::AbstractString) = Lion("gree", roar) # overloading constructor
type Panther <: Cat
	eye_color
	Panther() = new(green) # will only have this constructor
end

# -- multiple dispatch --
function meow(animal::Lion)
	animal.roar
end
function meow(animal::Panther)
	"grrr"
end
issubtype(Lion, Cat)

# -- under the hood
sqare_area(l) = l*l
code_native(sqare_area, (Int32,))
code_native(sqare_area, (Float64,))

# -- quick guide done

max_stars_repo_name: cdluminate/MyNotes

max_issues_repo_name: cdluminate/MyNotes

max_forks_repo_name: cdluminate/MyNotes