text:
# WARNING: This file was generated from parser.rl using ragel. Do not edit!
const samparser_start  = 70
const samparser_first_final  = 70
const samparser_error  = 0
const samparser_en_main  = 70
const _samparser_nfa_targs = Int8[0, 0 , ]
const _samparser_nfa_offsets = Int8[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 , ]
const _samparser_nfa_push_actions = Int8[0, 0 , ]
const _samparser_nfa_pop_trans = Int8[0, 0 , ]
Ragel.@generate_read!_function(
"samparser",
SAMReader,
SAMRecord,
begin
	if p == pe
		@goto _test_eof

	end
	if cs  == 70
		@goto st_case_70
	elseif cs  == 0
		@goto st_case_0
	elseif cs  == 1
		@goto st_case_1
	elseif cs  == 2
		@goto st_case_2
	elseif cs  == 3
		@goto st_case_3
	elseif cs  == 4
		@goto st_case_4
	elseif cs  == 5
		@goto st_case_5
	elseif cs  == 6
		@goto st_case_6
	elseif cs  == 7
		@goto st_case_7
	elseif cs  == 8
		@goto st_case_8
	elseif cs  == 9
		@goto st_case_9
	elseif cs  == 10
		@goto st_case_10
	elseif cs  == 11
		@goto st_case_11
	elseif cs  == 12
		@goto st_case_12
	elseif cs  == 13
		@goto st_case_13
	elseif cs  == 14
		@goto st_case_14
	elseif cs  == 15
		@goto st_case_15
	elseif cs  == 16
		@goto st_case_16
	elseif cs  == 17
		@goto st_case_17
	elseif cs  == 18
		@goto st_case_18
	elseif cs  == 19
		@goto st_case_19
	elseif cs  == 20
		@goto st_case_20
	elseif cs  == 21
		@goto st_case_21
	elseif cs  == 22
		@goto st_case_22
	elseif cs  == 23
		@goto st_case_23
	elseif cs  == 24
		@goto st_case_24
	elseif cs  == 25
		@goto st_case_25
	elseif cs  == 26
		@goto st_case_26
	elseif cs  == 27
		@goto st_case_27
	elseif cs  == 28
		@goto st_case_28
	elseif cs  == 29
		@goto st_case_29
	elseif cs  == 71
		@goto st_case_71
	elseif cs  == 30
		@goto st_case_30
	elseif cs  == 31
		@goto st_case_31
	elseif cs  == 32
		@goto st_case_32
	elseif cs  == 33
		@goto st_case_33
	elseif cs  == 34
		@goto st_case_34
	elseif cs  == 35
		@goto st_case_35
	elseif cs  == 36
		@goto st_case_36
	elseif cs  == 37
		@goto st_case_37
	elseif cs  == 38
		@goto st_case_38
	elseif cs  == 39
		@goto st_case_39
	elseif cs  == 40
		@goto st_case_40
	elseif cs  == 41
		@goto st_case_41
	elseif cs  == 42
		@goto st_case_42
	elseif cs  == 43
		@goto st_case_43
	elseif cs  == 44
		@goto st_case_44
	elseif cs  == 45
		@goto st_case_45
	elseif cs  == 46
		@goto st_case_46
	elseif cs  == 47
		@goto st_case_47
	elseif cs  == 48
		@goto st_case_48
	elseif cs  == 49
		@goto st_case_49
	elseif cs  == 50
		@goto st_case_50
	elseif cs  == 51
		@goto st_case_51
	elseif cs  == 52
		@goto st_case_52
	elseif cs  == 53
		@goto st_case_53
	elseif cs  == 54
		@goto st_case_54
	elseif cs  == 55
		@goto st_case_55
	elseif cs  == 56
		@goto st_case_56
	elseif cs  == 57
		@goto st_case_57
	elseif cs  == 58
		@goto st_case_58
	elseif cs  == 59
		@goto st_case_59
	elseif cs  == 60
		@goto st_case_60
	elseif cs  == 61
		@goto st_case_61
	elseif cs  == 62
		@goto st_case_62
	elseif cs  == 63
		@goto st_case_63
	elseif cs  == 64
		@goto st_case_64
	elseif cs  == 65
		@goto st_case_65
	elseif cs  == 66
		@goto st_case_66
	elseif cs  == 72
		@goto st_case_72
	elseif cs  == 67
		@goto st_case_67
	elseif cs  == 68
		@goto st_case_68
	elseif cs  == 69
		@goto st_case_69
	end
	@goto st_out
	@label st_case_70
	if (data[1+(p )])== 64
		@goto ctr84

	end
	if 33 <= (data[1+(p )])&& (data[1+(p )])<= 126
		@goto ctr83

	end
	@goto st0
	@label st_case_0
	@label st0
	cs = 0
	@goto _out
	@label ctr83
	Ragel.@anchor!
	@goto st1
	@label ctr85
	# need to anchor here for the next record
	Ragel.@anchor!
	Ragel.@yield 1
	Ragel.@anchor!
	@goto st1
	@label ctr86
	line = Ragel.@ascii_from_anchor!()
	tag = line[2:3]
	if tag == "HD"
		if haskey(input.header, "HD")
			error("found multiple @HD lines")
		end
		input.header[tag] = parse_keyvals(line[5:end-1])
else
	if !haskey(input.header, tag)
		input.header[tag] = []
	end
	if tag == "CO"
		push!(input.header[tag], line[5:end-1])
else
	push!(input.header[tag], parse_keyvals(line[5:end-1]))
end
		end
@goto yield
Ragel.@anchor!
@goto st1
@label st1
p+= 1
if p == pe
	@goto _test_eof1

end
@label st_case_1
if (data[1+(p )])== 9
	@goto ctr0

end
if (data[1+(p )])> 63
	if 	65 <= (data[1+(p )])&& (data[1+(p )])<= 126
	@goto st1

end

elseif 	(data[1+(p )])>= 33
	@goto st1

end
@goto st0
@label ctr0
Ragel.@copy_from_anchor!(output.name)
@goto st2
@label st2
p+= 1
if p == pe
	@goto _test_eof2

end
@label st_case_2
if 48 <= (data[1+(p )])&& (data[1+(p )])<= 57
	@goto ctr3

end
@goto st0
@label ctr3
Ragel.@anchor!
@goto st3
@label st3
p+= 1
if p == pe
	@goto _test_eof3

end
@label st_case_3
if (data[1+(p )])== 9
	@goto ctr4

end
if 48 <= (data[1+(p )])&& (data[1+(p )])<= 57
	@goto st3

end
@goto st0
@label ctr4
output.flag = Ragel.@int64_from_anchor!
@goto st4
@label st4
p+= 1
if p == pe
	@goto _test_eof4

end
@label st_case_4
if 33 <= (data[1+(p )])&& (data[1+(p )])<= 126
	@goto ctr6

end
@goto st0
@label ctr6
Ragel.@anchor!
@goto st5
@label st5
p+= 1
if p == pe
	@goto _test_eof5

end
@label st_case_5
if (data[1+(p )])== 9
	@goto ctr7

end
if 33 <= (data[1+(p )])&& (data[1+(p )])<= 126
	@goto st5

end
@goto st0
@label ctr7
Ragel.@copy_from_anchor!(output.refname)
@goto st6
@label st6
p+= 1
if p == pe
	@goto _test_eof6

end
@label st_case_6
if 48 <= (data[1+(p )])&& (data[1+(p )])<= 57
	@goto ctr9

end
@goto st0
@label ctr9
Ragel.@anchor!
@goto st7
@label st7
p+= 1
if p == pe
	@goto _test_eof7

end
@label st_case_7
if (data[1+(p )])== 9
	@goto ctr10

end
if 48 <= (data[1+(p )])&& (data[1+(p )])<= 57
	@goto st7

end
@goto st0
@label ctr10
output.pos = Ragel.@int64_from_anchor!
@goto st8
@label st8
p+= 1
if p == pe
	@goto _test_eof8

end
@label st_case_8
if 48 <= (data[1+(p )])&& (data[1+(p )])<= 57
	@goto ctr12

end
@goto st0
@label ctr12
Ragel.@anchor!
@goto st9
@label st9
p+= 1
if p == pe
	@goto _test_eof9

end
@label st_case_9
if (data[1+(p )])== 9
	@goto ctr13

end
if 48 <= (data[1+(p )])&& (data[1+(p )])<= 57
	@goto st9

end
@goto st0
@label ctr13
output.mapq = Ragel.@int64_from_anchor!
@goto st10
@label st10
p+= 1
if p == pe
	@goto _test_eof10

end
@label st_case_10
if (data[1+(p )])== 42
	@goto ctr15

end
if 48 <= (data[1+(p )])&& (data[1+(p )])<= 57
	@goto ctr16

end
@goto st0
@label ctr15
Ragel.@anchor!
@goto st11
@label st11
p+= 1
if p == pe
	@goto _test_eof11

end
@label st_case_11
if (data[1+(p )])== 9
	@goto ctr17

end
@goto st0
@label ctr17
Ragel.@copy_from_anchor!(output.cigar)
@goto st12
@label st12
p+= 1
if p == pe
	@goto _test_eof12

end
@label st_case_12
if 33 <= (data[1+(p )])&& (data[1+(p )])<= 126
	@goto ctr18

end
@goto st0
@label ctr18
Ragel.@anchor!
@goto st13
@label st13
p+= 1
if p == pe
	@goto _test_eof13

end
@label st_case_13
if (data[1+(p )])== 9
	@goto ctr19

end
if 33 <= (data[1+(p )])&& (data[1+(p )])<= 126
	@goto st13

end
@goto st0
@label ctr19
Ragel.@copy_from_anchor!(output.next_refname)
@goto st14
@label st14
p+= 1
if p == pe
	@goto _test_eof14

end
@label st_case_14
if 48 <= (data[1+(p )])&& (data[1+(p )])<= 57
	@goto ctr21

end
@goto st0
@label ctr21
Ragel.@anchor!
@goto st15
@label st15
p+= 1
if p == pe
	@goto _test_eof15

end
@label st_case_15
if (data[1+(p )])== 9
	@goto ctr22

end
if 48 <= (data[1+(p )])&& (data[1+(p )])<= 57
	@goto st15

end
@goto st0
@label ctr22
output.next_pos = Ragel.@int64_from_anchor!
@goto st16
@label st16
p+= 1
if p == pe
	@goto _test_eof16

end
@label st_case_16
if (data[1+(p )])== 45
	@goto ctr24

end
if 48 <= (data[1+(p )])&& (data[1+(p )])<= 57
	@goto ctr25

end
@goto st0
@label ctr24
Ragel.@anchor!
@goto st17
@label st17
p+= 1
if p == pe
	@goto _test_eof17

end
@label st_case_17
if 48 <= (data[1+(p )])&& (data[1+(p )])<= 57
	@goto st18

end
@goto st0
@label ctr25
Ragel.@anchor!
@goto st18
@label st18
p+= 1
if p == pe
	@goto _test_eof18

end
@label st_case_18
if (data[1+(p )])== 9
	@goto ctr27

end
if 48 <= (data[1+(p )])&& (data[1+(p )])<= 57
	@goto st18

end
@goto st0
@label ctr27
output.tlen = Ragel.@int64_from_anchor!
@goto st19
@label st19
p+= 1
if p == pe
	@goto _test_eof19

end
@label st_case_19
if (data[1+(p )]) == 42
	@goto ctr28
elseif (data[1+(p )]) == 46
	@goto ctr29
elseif (data[1+(p )]) == 61
	@goto ctr29
end
if (data[1+(p )])> 90
	if 97 <= (data[1+(p )])&& (data[1+(p )])<= 122
		@goto ctr29

	end

elseif (data[1+(p )])>= 65
	@goto ctr29

end
@goto st0
@label ctr28
Ragel.@anchor!
@goto st20
@label st20
p+= 1
if p == pe
	@goto _test_eof20

end
@label st_case_20
if (data[1+(p )])== 9
	@goto ctr30

end
@goto st0
@label ctr30
Ragel.@copy_from_anchor!(output.seq)
@goto st21
@label st21
p+= 1
if p == pe
	@goto _test_eof21

end
@label st_case_21
if 33 <= (data[1+(p )])&& (data[1+(p )])<= 126
	@goto ctr31

end
@goto st0
@label ctr31
Ragel.@anchor!
@goto st22
@label st22
p+= 1
if p == pe
	@goto _test_eof22

end
@label st_case_22
if (data[1+(p )]) == 9
	@goto ctr32
elseif (data[1+(p )]) == 10
	@goto ctr33
end
if 33 <= (data[1+(p )])&& (data[1+(p )])<= 126
	@goto st22

end
@goto st0
@label ctr32
Ragel.@copy_from_anchor!(output.qual)
empty!(output.optional_fields)
@goto st23
@label ctr46
optfieldstr = Ragel.@ascii_from_anchor!()
tag = optfieldstr[1:2]
typ = optfieldstr[4]
if typ == 'A'
	value = optfieldstr[6]
elseif typ == 'Z'
	value = optfieldstr[6:end]
elseif 	typ == 'H'
	value = parse_hexbytearray(optfieldstr[6:end])
elseif 	typ == 'B'
	eltyp = auxtype[UInt8[optfieldstr[6]]]
	value = [parse(eltyp, x) for x in split(optfieldstr[8:end], ',')]
else
			value = parse(auxtype[UInt8(typ)], optfieldstr[6:end])
end
		output.optional_fields[tag] = value
	@goto st23
@label st23
p+= 1
if 	p == pe
	@goto _test_eof23

end
@label st_case_23
if (data[1+(p )])> 90
	if 97 <= (data[1+(p )])&& (data[1+(p )])<= 122
		@goto ctr35

	end

elseif (data[1+(p )])>= 65
	@goto ctr35

end
@goto st0
@label ctr35
Ragel.@anchor!
@goto st24
@label st24
p+= 1
if p == pe
	@goto _test_eof24

end
@label st_case_24
if (data[1+(p )])< 65
	if 48 <= (data[1+(p )])&& (data[1+(p )])<= 57
		@goto st25

	end

elseif (data[1+(p )])> 90
	if 97 <= (data[1+(p )])&& (data[1+(p )])<= 122
		@goto st25

	end

else
	@goto st25

end
@goto st0
@label st25
p+= 1
if p == pe
	@goto _test_eof25

end
@label st_case_25
if (data[1+(p )])== 58
	@goto st26

end
@goto st0
@label st26
p+= 1
if p == pe
	@goto _test_eof26

end
@label st_case_26
if (data[1+(p )]) == 65
	@goto st27
elseif (data[1+(p )]) == 66
	@goto st30
elseif (data[1+(p )]) == 67
	@goto st40
elseif (data[1+(p )]) == 72
	@goto st44
elseif (data[1+(p )]) == 73
	@goto st40
elseif (data[1+(p )]) == 83
	@goto st40
elseif (data[1+(p )]) == 90
	@goto st47
elseif (data[1+(p )]) == 99
	@goto st40
elseif (data[1+(p )]) == 102
	@goto st50
elseif (data[1+(p )]) == 105
	@goto st40
elseif (data[1+(p )]) == 115
	@goto st40
end
@goto st0
@label st27
p+= 1
if p == pe
	@goto _test_eof27

end
@label st_case_27
if (data[1+(p )])== 58
	@goto st28

end
@goto st0
@label st28
p+= 1
if p == pe
	@goto _test_eof28

end
@label st_case_28
if 33 <= (data[1+(p )])&& (data[1+(p )])<= 126
	@goto st29

end
@goto st0
@label st29
p+= 1
if p == pe
	@goto _test_eof29

end
@label st_case_29
if (data[1+(p )]) == 9
	@goto ctr46
elseif (data[1+(p )]) == 10
	@goto ctr47
end
@goto st0
@label ctr33
Ragel.@copy_from_anchor!(output.qual)
empty!(output.optional_fields)
input.state.linenum += 1
@goto st71
@label ctr47
optfieldstr = Ragel.@ascii_from_anchor!()
tag = optfieldstr[1:2]
typ = optfieldstr[4]
if typ == 'A'
	value = optfieldstr[6]
elseif typ == 'Z'
	value = optfieldstr[6:end]
elseif 	typ == 'H'
	value = parse_hexbytearray(optfieldstr[6:end])
elseif 	typ == 'B'
	eltyp = auxtype[UInt8[optfieldstr[6]]]
	value = [parse(eltyp, x) for x in split(optfieldstr[8:end], ',')]
else
			value = parse(auxtype[UInt8(typ)], optfieldstr[6:end])
end
		output.optional_fields[tag] = value
	input.state.linenum += 1
	@goto st71
@label st71
p+= 1
if 	p == pe
	@goto _test_eof71

end
@label st_case_71
if (data[1+(p )])> 63
	if 65 <= (data[1+(p )])&& (data[1+(p )])<= 126
		@goto ctr85

	end

elseif (data[1+(p )])>= 33
	@goto ctr85

end
@goto st0
@label st30
p+= 1
if p == pe
	@goto _test_eof30

end
@label st_case_30
if (data[1+(p )])== 58
	@goto st31

end
@goto st0
@label st31
p+= 1
if p == pe
	@goto _test_eof31

end
@label st_case_31
if (data[1+(p )])== 44
	@goto st32

end
@goto st0
@label st32
p+= 1
if p == pe
	@goto _test_eof32

end
@label st_case_32
if (data[1+(p )]) == 43
	@goto st33
elseif (data[1+(p )]) == 45
	@goto st33
elseif (data[1+(p )]) == 46
	@goto st34
end
if 48 <= (data[1+(p )])&& (data[1+(p )])<= 57
	@goto st39

end
@goto st0
@label st33
p+= 1
if p == pe
	@goto _test_eof33

end
@label st_case_33
if (data[1+(p )])== 46
	@goto st34

end
if 48 <= (data[1+(p )])&& (data[1+(p )])<= 57
	@goto st39

end
@goto st0
@label st34
p+= 1
if p == pe
	@goto _test_eof34

end
@label st_case_34
if 48 <= (data[1+(p )])&& (data[1+(p )])<= 57
	@goto st35

end
@goto st0
@label st35
p+= 1
if p == pe
	@goto _test_eof35

end
@label st_case_35
if (data[1+(p )]) == 9
	@goto ctr46
elseif (data[1+(p )]) == 10
	@goto ctr47
elseif (data[1+(p )]) == 44
	@goto st32
elseif (data[1+(p )]) == 69
	@goto st36
elseif (data[1+(p )]) == 101
	@goto st36
end
if 48 <= (data[1+(p )])&& (data[1+(p )])<= 57
	@goto st35

end
@goto st0
@label st36
p+= 1
if p == pe
	@goto _test_eof36

end
@label st_case_36
if (data[1+(p )]) == 43
	@goto st37
elseif (data[1+(p )]) == 45
	@goto st37
end
if 48 <= (data[1+(p )])&& (data[1+(p )])<= 57
	@goto st38

end
@goto st0
@label st37
p+= 1
if p == pe
	@goto _test_eof37

end
@label st_case_37
if 48 <= (data[1+(p )])&& (data[1+(p )])<= 57
	@goto st38

end
@goto st0
@label st38
p+= 1
if p == pe
	@goto _test_eof38

end
@label st_case_38
if (data[1+(p )]) == 9
	@goto ctr46
elseif (data[1+(p )]) == 10
	@goto ctr47
elseif (data[1+(p )]) == 44
	@goto st32
end
if 48 <= (data[1+(p )])&& (data[1+(p )])<= 57
	@goto st38

end
@goto st0
@label st39
p+= 1
if p == pe
	@goto _test_eof39

end
@label st_case_39
if (data[1+(p )]) == 9
	@goto ctr46
elseif (data[1+(p )]) == 10
	@goto ctr47
elseif (data[1+(p )]) == 44
	@goto st32
elseif (data[1+(p )]) == 46
	@goto st34
elseif (data[1+(p )]) == 69
	@goto st36
elseif (data[1+(p )]) == 101
	@goto st36
end
if 48 <= (data[1+(p )])&& (data[1+(p )])<= 57
	@goto st39

end
@goto st0
@label st40
p+= 1
if p == pe
	@goto _test_eof40

end
@label st_case_40
if (data[1+(p )])== 58
	@goto st41

end
@goto st0
@label st41
p+= 1
if p == pe
	@goto _test_eof41

end
@label st_case_41
if (data[1+(p )]) == 43
	@goto st42
elseif (data[1+(p )]) == 45
	@goto st42
end
if 48 <= (data[1+(p )])&& (data[1+(p )])<= 57
	@goto st43

end
@goto st0
@label st42
p+= 1
if p == pe
	@goto _test_eof42

end
@label st_case_42
if 48 <= (data[1+(p )])&& (data[1+(p )])<= 57
	@goto st43

end
@goto st0
@label st43
p+= 1
if p == pe
	@goto _test_eof43

end
@label st_case_43
if (data[1+(p )]) == 9
	@goto ctr46
elseif (data[1+(p )]) == 10
	@goto ctr47
end
if 48 <= (data[1+(p )])&& (data[1+(p )])<= 57
	@goto st43

end
@goto st0
@label st44
p+= 1
if p == pe
	@goto _test_eof44

end
@label st_case_44
if (data[1+(p )])== 58
	@goto st45

end
@goto st0
@label st45
p+= 1
if p == pe
	@goto _test_eof45

end
@label st_case_45
if (data[1+(p )])> 57
	if 65 <= (data[1+(p )])&& (data[1+(p )])<= 70
		@goto st46

	end

elseif (data[1+(p )])>= 48
	@goto st46

end
@goto st0
@label st46
p+= 1
if p == pe
	@goto _test_eof46

end
@label st_case_46
if (data[1+(p )]) == 9
	@goto ctr46
elseif (data[1+(p )]) == 10
	@goto ctr47
end
if (data[1+(p )])> 57
	if 65 <= (data[1+(p )])&& (data[1+(p )])<= 70
		@goto st46

	end

elseif (data[1+(p )])>= 48
	@goto st46

end
@goto st0
@label st47
p+= 1
if p == pe
	@goto _test_eof47

end
@label st_case_47
if (data[1+(p )])== 58
	@goto st48

end
@goto st0
@label st48
p+= 1
if p == pe
	@goto _test_eof48

end
@label st_case_48
if 32 <= (data[1+(p )])&& (data[1+(p )])<= 126
	@goto st49

end
@goto st0
@label st49
p+= 1
if p == pe
	@goto _test_eof49

end
@label st_case_49
if (data[1+(p )]) == 9
	@goto ctr46
elseif (data[1+(p )]) == 10
	@goto ctr47
end
if 32 <= (data[1+(p )])&& (data[1+(p )])<= 126
	@goto st49

end
@goto st0
@label st50
p+= 1
if p == pe
	@goto _test_eof50

end
@label st_case_50
if (data[1+(p )])== 58
	@goto st51

end
@goto st0
@label st51
p+= 1
if p == pe
	@goto _test_eof51

end
@label st_case_51
if (data[1+(p )]) == 43
	@goto st52
elseif (data[1+(p )]) == 45
	@goto st52
elseif (data[1+(p )]) == 46
	@goto st53
end
if 48 <= (data[1+(p )])&& (data[1+(p )])<= 57
	@goto st55

end
@goto st0
@label st52
p+= 1
if p == pe
	@goto _test_eof52

end
@label st_case_52
if (data[1+(p )])== 46
	@goto st53

end
if 48 <= (data[1+(p )])&& (data[1+(p )])<= 57
	@goto st55

end
@goto st0
@label st53
p+= 1
if p == pe
	@goto _test_eof53

end
@label st_case_53
if 48 <= (data[1+(p )])&& (data[1+(p )])<= 57
	@goto st54

end
@goto st0
@label st54
p+= 1
if p == pe
	@goto _test_eof54

end
@label st_case_54
if (data[1+(p )]) == 9
	@goto ctr46
elseif (data[1+(p )]) == 10
	@goto ctr47
elseif (data[1+(p )]) == 69
	@goto st41
elseif (data[1+(p )]) == 101
	@goto st41
end
if 48 <= (data[1+(p )])&& (data[1+(p )])<= 57
	@goto st54

end
@goto st0
@label st55
p+= 1
if p == pe
	@goto _test_eof55

end
@label st_case_55
if (data[1+(p )]) == 9
	@goto ctr46
elseif (data[1+(p )]) == 10
	@goto ctr47
elseif (data[1+(p )]) == 46
	@goto st53
elseif (data[1+(p )]) == 69
	@goto st41
elseif (data[1+(p )]) == 101
	@goto st41
end
if 48 <= (data[1+(p )])&& (data[1+(p )])<= 57
	@goto st55

end
@goto st0
@label ctr29
Ragel.@anchor!
@goto st56
@label st56
p+= 1
if p == pe
	@goto _test_eof56

end
@label st_case_56
if (data[1+(p )]) == 9
	@goto ctr30
elseif (data[1+(p )]) == 46
	@goto st56
elseif (data[1+(p )]) == 61
	@goto st56
end
if (data[1+(p )])> 90
	if 97 <= (data[1+(p )])&& (data[1+(p )])<= 122
		@goto st56

	end

elseif (data[1+(p )])>= 65
	@goto st56

end
@goto st0
@label ctr16
Ragel.@anchor!
@goto st57
@label st57
p+= 1
if p == pe
	@goto _test_eof57

end
@label st_case_57
if (data[1+(p )]) == 61
	@goto st58
elseif (data[1+(p )]) == 68
	@goto st58
elseif (data[1+(p )]) == 80
	@goto st58
elseif (data[1+(p )]) == 83
	@goto st58
elseif (data[1+(p )]) == 88
	@goto st58
end
if (data[1+(p )])< 72
	if 48 <= (data[1+(p )])&& (data[1+(p )])<= 57
		@goto st57

	end

elseif (data[1+(p )])> 73
	if 77 <= (data[1+(p )])&& (data[1+(p )])<= 78
		@goto st58

	end

else
	@goto st58

end
@goto st0
@label st58
p+= 1
if p == pe
	@goto _test_eof58

end
@label st_case_58
if (data[1+(p )])== 9
	@goto ctr17

end
if 48 <= (data[1+(p )])&& (data[1+(p )])<= 57
	@goto st57

end
@goto st0
@label ctr84
Ragel.@anchor!
@goto st59
@label ctr87
line = Ragel.@ascii_from_anchor!()
tag = line[2:3]
if tag == "HD"
	if haskey(input.header, "HD")
		error("found multiple @HD lines")
	end
	input.header[tag] = parse_keyvals(line[5:end-1])
else
			if !haskey(input.header, tag)
		input.header[tag] = []
	end
	if tag == "CO"
		push!(input.header[tag], line[5:end-1])
else
	push!(input.header[tag], parse_keyvals(line[5:end-1]))
end
		end
@goto yield
Ragel.@anchor!
@goto st59
@label st59
p+= 1
if p == pe
	@goto _test_eof59

end
@label st_case_59
if (data[1+(p )])== 67
	@goto st67

end
if (data[1+(p )])> 90
	if 	97 <= (data[1+(p )])&& (data[1+(p )])<= 122
	@goto st60

end

elseif 	(data[1+(p )])>= 65
	@goto st60

end
@goto st0
@label st60
p+= 1
if p == pe
	@goto _test_eof60

end
@label st_case_60
if (data[1+(p )])< 65
	if 48 <= (data[1+(p )])&& (data[1+(p )])<= 57
		@goto st61

	end

elseif (data[1+(p )])> 90
	if 97 <= (data[1+(p )])&& (data[1+(p )])<= 122
		@goto st61

	end

else
	@goto st61

end
@goto st0
@label st61
p+= 1
if p == pe
	@goto _test_eof61

end
@label st_case_61
if (data[1+(p )])== 9
	@goto st62

end
@goto st0
@label st62
p+= 1
if p == pe
	@goto _test_eof62

end
@label st_case_62
if (data[1+(p )])> 90
	if 97 <= (data[1+(p )])&& (data[1+(p )])<= 122
		@goto st63

	end

elseif (data[1+(p )])>= 65
	@goto st63

end
@goto st0
@label st63
p+= 1
if p == pe
	@goto _test_eof63

end
@label st_case_63
if (data[1+(p )])< 65
	if 48 <= (data[1+(p )])&& (data[1+(p )])<= 57
		@goto st64

	end

elseif (data[1+(p )])> 90
	if 97 <= (data[1+(p )])&& (data[1+(p )])<= 122
		@goto st64

	end

else
	@goto st64

end
@goto st0
@label st64
p+= 1
if p == pe
	@goto _test_eof64

end
@label st_case_64
if (data[1+(p )])== 58
	@goto st65

end
@goto st0
@label st65
p+= 1
if p == pe
	@goto _test_eof65

end
@label st_case_65
if 32 <= (data[1+(p )])&& (data[1+(p )])<= 126
	@goto st66

end
@goto st0
@label st66
p+= 1
if p == pe
	@goto _test_eof66

end
@label st_case_66
if (data[1+(p )]) == 9
	@goto st62
elseif (data[1+(p )]) == 10
	@goto ctr80
end
if 32 <= (data[1+(p )])&& (data[1+(p )])<= 126
	@goto st66

end
@goto st0
@label ctr80
input.state.linenum += 1
@goto st72
@label st72
p+= 1
if p == pe
	@goto _test_eof72

end
@label st_case_72
if (data[1+(p )])== 64
	@goto ctr87

end
if 33 <= (data[1+(p )])&& (data[1+(p )])<= 126
	@goto ctr86

end
@goto st0
@label st67
p+= 1
if p == pe
	@goto _test_eof67

end
@label st_case_67
if (data[1+(p )])== 79
	@goto st68

end
if (data[1+(p )])< 65
	if 48 <= (data[1+(p )])&& (data[1+(p )])<= 57
		@goto st61

	end

elseif (data[1+(p )])> 90
	if 97 <= (data[1+(p )])&& (data[1+(p )])<= 122
		@goto st61

	end

else
	@goto st61

end
@goto st0
@label st68
p+= 1
if p == pe
	@goto _test_eof68

end
@label st_case_68
if (data[1+(p )])== 9
	@goto st69

end
@goto st0
@label st69
p+= 1
if p == pe
	@goto _test_eof69

end
@label st_case_69
if (data[1+(p )])== 10
	@goto ctr80

end
@goto st69
@label st_out
@label _test_eof1
cs = 1
@goto _test_eof
@label _test_eof2
cs = 2
@goto _test_eof
@label _test_eof3
cs = 3
@goto _test_eof
@label _test_eof4
cs = 4
@goto _test_eof
@label _test_eof5
cs = 5
@goto _test_eof
@label _test_eof6
cs = 6
@goto _test_eof
@label _test_eof7
cs = 7
@goto _test_eof
@label _test_eof8
cs = 8
@goto _test_eof
@label _test_eof9
cs = 9
@goto _test_eof
@label _test_eof10
cs = 10
@goto _test_eof
@label _test_eof11
cs = 11
@goto _test_eof
@label _test_eof12
cs = 12
@goto _test_eof
@label _test_eof13
cs = 13
@goto _test_eof
@label _test_eof14
cs = 14
@goto _test_eof
@label _test_eof15
cs = 15
@goto _test_eof
@label _test_eof16
cs = 16
@goto _test_eof
@label _test_eof17
cs = 17
@goto _test_eof
@label _test_eof18
cs = 18
@goto _test_eof
@label _test_eof19
cs = 19
@goto _test_eof
@label _test_eof20
cs = 20
@goto _test_eof
@label _test_eof21
cs = 21
@goto _test_eof
@label _test_eof22
cs = 22
@goto _test_eof
@label _test_eof23
cs = 23
@goto _test_eof
@label _test_eof24
cs = 24
@goto _test_eof
@label _test_eof25
cs = 25
@goto _test_eof
@label _test_eof26
cs = 26
@goto _test_eof
@label _test_eof27
cs = 27
@goto _test_eof
@label _test_eof28
cs = 28
@goto _test_eof
@label _test_eof29
cs = 29
@goto _test_eof
@label _test_eof71
cs = 71
@goto _test_eof
@label _test_eof30
cs = 30
@goto _test_eof
@label _test_eof31
cs = 31
@goto _test_eof
@label _test_eof32
cs = 32
@goto _test_eof
@label _test_eof33
cs = 33
@goto _test_eof
@label _test_eof34
cs = 34
@goto _test_eof
@label _test_eof35
cs = 35
@goto _test_eof
@label _test_eof36
cs = 36
@goto _test_eof
@label _test_eof37
cs = 37
@goto _test_eof
@label _test_eof38
cs = 38
@goto _test_eof
@label _test_eof39
cs = 39
@goto _test_eof
@label _test_eof40
cs = 40
@goto _test_eof
@label _test_eof41
cs = 41
@goto _test_eof
@label _test_eof42
cs = 42
@goto _test_eof
@label _test_eof43
cs = 43
@goto _test_eof
@label _test_eof44
cs = 44
@goto _test_eof
@label _test_eof45
cs = 45
@goto _test_eof
@label _test_eof46
cs = 46
@goto _test_eof
@label _test_eof47
cs = 47
@goto _test_eof
@label _test_eof48
cs = 48
@goto _test_eof
@label _test_eof49
cs = 49
@goto _test_eof
@label _test_eof50
cs = 50
@goto _test_eof
@label _test_eof51
cs = 51
@goto _test_eof
@label _test_eof52
cs = 52
@goto _test_eof
@label _test_eof53
cs = 53
@goto _test_eof
@label _test_eof54
cs = 54
@goto _test_eof
@label _test_eof55
cs = 55
@goto _test_eof
@label _test_eof56
cs = 56
@goto _test_eof
@label _test_eof57
cs = 57
@goto _test_eof
@label _test_eof58
cs = 58
@goto _test_eof
@label _test_eof59
cs = 59
@goto _test_eof
@label _test_eof60
cs = 60
@goto _test_eof
@label _test_eof61
cs = 61
@goto _test_eof
@label _test_eof62
cs = 62
@goto _test_eof
@label _test_eof63
cs = 63
@goto _test_eof
@label _test_eof64
cs = 64
@goto _test_eof
@label _test_eof65
cs = 65
@goto _test_eof
@label _test_eof66
cs = 66
@goto _test_eof
@label _test_eof72
cs = 72
@goto _test_eof
@label _test_eof67
cs = 67
@goto _test_eof
@label _test_eof68
cs = 68
@goto _test_eof
@label _test_eof69
cs = 69
@goto _test_eof
@label _test_eof
if p == eof
	if cs  == 72
		line = Ragel.@ascii_from_anchor!()
		tag = line[2:3]
		if tag == "HD"
			if haskey(input.header, "HD")
				error("found multiple @HD lines")
			end
			input.header[tag] = parse_keyvals(line[5:end-1])
	else
		if !haskey(input.header, tag)
			input.header[tag] = []
		end
		if tag == "CO"
			push!(input.header[tag], line[5:end-1])
	else
		push!(input.header[tag], parse_keyvals(line[5:end-1]))
end
end
				@goto yield
			elseif cs  == 71
# need to anchor here for the next record
				Ragel.@anchor!
				Ragel.@yield 0
			end

end
@label _out
end)

# header only

const samheaderparser_start  = 12
const samheaderparser_first_final  = 12
const samheaderparser_error  = 0
const samheaderparser_en_main  = 12
const _samheaderparser_nfa_targs = Int8[0, 0 , ]
const _samheaderparser_nfa_offsets = Int8[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 , ]
const _samheaderparser_nfa_push_actions = Int8[0, 0 , ]
const _samheaderparser_nfa_pop_trans = Int8[0, 0 , ]
function parse_samheader(data)
header = SAMHeader()

p = 0
pe = eof = endof(data)
anchor = 0

cs = convert(Int, samheaderparser_start )
if p == pe
@goto _test_eof

end
if cs  == 12
@goto st_case_12
elseif cs  == 0
@goto st_case_0
elseif cs  == 1
@goto st_case_1
elseif cs  == 2
@goto st_case_2
elseif cs  == 3
@goto st_case_3
elseif cs  == 4
@goto st_case_4
elseif cs  == 5
@goto st_case_5
elseif cs  == 6
@goto st_case_6
elseif cs  == 7
@goto st_case_7
elseif cs  == 8
@goto st_case_8
elseif cs  == 13
@goto st_case_13
elseif cs  == 9
@goto st_case_9
elseif cs  == 10
@goto st_case_10
elseif cs  == 11
@goto st_case_11
end
@goto st_out
@label st_case_12
if (data[1+(p )])== 64
@goto ctr12

end
@goto st0
@label st_case_0
@label st0
cs = 0
@goto _out
@label ctr12
anchor = p + 1
@goto st1
@label ctr13
line = String(data[anchor:p])
tag = line[2:3]
if tag == "HD"
if haskey(header, "HD")
		error("found multiple @HD lines")
		end
header[tag] = parse_keyvals(line[5:end-1])
else
if !haskey(header, tag)
header[tag] = []
end
if tag == "CO"
push!(header[tag], line[5:end-1])
else
push!(header[tag], parse_keyvals(line[5:end-1]))
end
end
anchor = p + 1
@goto st1
@label st1
p+= 1
if p == pe
@goto _test_eof1

end
@label st_case_1
if (data[1+(p )])== 67
@goto st9

end
if (data[1+(p )])> 90
if 97 <= (data[1+(p )])&& (data[1+(p )])<= 122
@goto st2

end

elseif (data[1+(p )])>= 65
@goto st2

end
@goto st0
@label st2
p+= 1
if p == pe
@goto _test_eof2

end
@label st_case_2
if (data[1+(p )])< 65
if 48 <= (data[1+(p )])&& (data[1+(p )])<= 57
@goto st3

end

elseif (data[1+(p )])> 90
if 97 <= (data[1+(p )])&& (data[1+(p )])<= 122
@goto st3

end

else
@goto st3

end
@goto st0
@label st3
p+= 1
if p == pe
@goto _test_eof3

end
@label st_case_3
if (data[1+(p )])== 9
@goto st4

end
@goto st0
@label st4
p+= 1
if p == pe
@goto _test_eof4

end
@label st_case_4
if (data[1+(p )])> 90
if 97 <= (data[1+(p )])&& (data[1+(p )])<= 122
@goto st5

end

elseif (data[1+(p )])>= 65
@goto st5

end
@goto st0
@label st5
p+= 1
if p == pe
@goto _test_eof5

end
@label st_case_5
if (data[1+(p )])< 65
if 48 <= (data[1+(p )])&& (data[1+(p )])<= 57
@goto st6

end

elseif (data[1+(p )])> 90
if 97 <= (data[1+(p )])&& (data[1+(p )])<= 122
@goto st6

end

else
@goto st6

end
@goto st0
@label st6
p+= 1
if p == pe
@goto _test_eof6

end
@label st_case_6
if (data[1+(p )])== 58
@goto st7

end
@goto st0
@label st7
p+= 1
if p == pe
@goto _test_eof7

end
@label st_case_7
if 32 <= (data[1+(p )])&& (data[1+(p )])<= 126
@goto st8

end
@goto st0
@label st8
p+= 1
if p == pe
@goto _test_eof8

end
@label st_case_8
if (data[1+(p )]) == 9
@goto st4
elseif (data[1+(p )]) == 10
@goto st13
end
if 32 <= (data[1+(p )])&& (data[1+(p )])<= 126
@goto st8

end
@goto st0
@label st13
p+= 1
if p == pe
@goto _test_eof13

end
@label st_case_13
if (data[1+(p )])== 64
@goto ctr13

end
@goto st0
@label st9
p+= 1
if p == pe
@goto _test_eof9

end
@label st_case_9
if (data[1+(p )])== 79
@goto st10

end
if (data[1+(p )])< 65
if 48 <= (data[1+(p )])&& (data[1+(p )])<= 57
@goto st3

end

elseif (data[1+(p )])> 90
if 97 <= (data[1+(p )])&& (data[1+(p )])<= 122
@goto st3

end

else
@goto st3

end
@goto st0
@label st10
p+= 1
if p == pe
@goto _test_eof10

end
@label st_case_10
if (data[1+(p )])== 9
@goto st11

end
@goto st0
@label st11
p+= 1
if p == pe
@goto _test_eof11

end
@label st_case_11
if (data[1+(p )])== 10
@goto st13

end
@goto st11
@label st_out
@label _test_eof1
cs = 1
@goto _test_eof
@label _test_eof2
cs = 2
@goto _test_eof
@label _test_eof3
cs = 3
@goto _test_eof
@label _test_eof4
cs = 4
@goto _test_eof
@label _test_eof5
cs = 5
@goto _test_eof
@label _test_eof6
cs = 6
@goto _test_eof
@label _test_eof7
cs = 7
@goto _test_eof
@label _test_eof8
cs = 8
@goto _test_eof
@label _test_eof13
cs = 13
@goto _test_eof
@label _test_eof9
cs = 9
@goto _test_eof
@label _test_eof10
cs = 10
@goto _test_eof
@label _test_eof11
cs = 11
@goto _test_eof
@label _test_eof
if p == eof
if cs  == 13
line = String(data[anchor:p])
tag = line[2:3]
if tag == "HD"
if haskey(header, "HD")
error("found multiple @HD lines")
end
header[tag] = parse_keyvals(line[5:end-1])
else
if !haskey(header, tag)
header[tag] = []
end
if tag == "CO"
push!(header[tag], line[5:end-1])
else
push!(header[tag], parse_keyvals(line[5:end-1]))
end
end
end

end
@label _out
@assert cs >= samheaderparser_first_final

return header
end


text:
using PyCall
using PyPlot
@pyimport matplotlib.animation as anim

#surf(disc.disc_data.flow.xi,t,yout[:,1:N])

plt = PyPlot
# First set up the figure, the axis, and the plot element we want to animate
fig = plt.figure()
#ax = plt.axes(xlim=(0, 0.5), ylim=(0.00323439, 0.00323441))
ax = plt.axes(xlim=(0, 0.5), ylim=(0.000, 0.005))
global lines = [ax[:plot]([], [], lw=2)[1];ax[:plot]([], [], lw=2)[1]]

# initialization function: plot the background of each frame
function init()
    global lines
	for line = lines
		line[:set_data]([], [])
	end
	plt.legend(["Nonlinear", "Linear"])
    return (lines,None)
end

# animation function.  This is called sequentially
function animate(i)
	k = 1
	global lines
	for line = lines
		x = [xi]
		if k == 1
			y = [yout[i*4,N+1:end-1][:]]
		else
			y = [youtlin[i*4,N+1:end-1][:]]
		end
		line[:set_data](x, y)
		#line = plt.plot(x,y)
		k += 1
	end
    return (lines,None)
end

# call the animator.  blit=True means only re-draw the parts that have changed.
myanim = anim.FuncAnimation(fig, animate, init_func=init,
                               frames=2500, interval=1);

#myanim[:save]("PyPlots-sinplot.mp4", extra_args=["-vcodec", "libx264", "-pix_fmt", "yuv420p"])

ext: jl

max_stars_repo_name: flavioluiz/PortHamiltonian

max_issues_repo_name: flavioluiz/PortHamiltonian

max_forks_repo_name: flavioluiz/PortHamiltonian


text:
# Function to compute log(sum(exp.(x))) using [Sebastian Nowozin: Streaming Log-sum-exp Computation.](http://www.nowozin.net/sebastian/blog/streaming-log-sum-exp-computation.html)

_f1(x1, x2) = x1 == x2 ? zero(x1 - x2) : -abs(x1 - x2)
_f2(x1, x2, a) = x1 > x2 ? oftype(a, x1) : oftype(a, x2)
_f3(x1, x2, r, a) = x1 > x2 ? (r + one(r)) * exp(a) : r + exp(a)
_f4(x1, x2, r1, r2, a) = x1 > x2 ? r1 + (r2 + one(r2)) * exp(a) : r2 + (r1 + one(r1)) * exp(a)

function _logsumexp_onepass_op(x1::AbstractVector, x2::AbstractVector)
    a = _f1.(x1,x2)
    xmax = _f2.(x1,x2,a)
    r = exp.(a)
    return xmax, r
end

function _logsumexp_onepass_op(x, (xmax, r)::Tuple)
    a = _f1.(x, xmax)
    _xmax = _f2.(x,xmax,a)
    _r = _f3.(x, xmax, r, a)
    return _xmax, _r
end
_logsumexp_onepass_op(xmax_r::Tuple, x) = _logsumexp_onepass_op(x, xmax_r)

function _logsumexp_onepass_op((xmax1, r1)::Tuple, (xmax2, r2)::Tuple)
    a = _f1.(xmax1, xmax2)
    xmax = _f2.(xmax1, xmax2)
    r = _f4.(xmax1, xmax2, r1, r2, a)
    return xmax, r
end


text:
"""
Time stepping solver

"""
function step!(
    KS::T1,
    uq::UQ1D,
    faceL::Interface1D1F,
    cell::ControlVolume1D1F,
    faceR::Interface1D1F,
    dt,
    RES,
    AVG,
    coll = :bgk::Symbol,
) where {T1<:AbstractSolverSet} # 1D1F1V

    if uq.method == "galerkin"

        #--- update conservative flow variables: step 1 ---#
        # w^n
        w_old = deepcopy(cell.w)
        prim_old = deepcopy(cell.prim)

        # flux -> w^{n+1}
        @. cell.w += (faceL.fw - faceR.fw) / cell.dx
        cell.prim .= uq_conserve_prim(cell.w, KS.gas.γ, uq)

        # locate variables on random quadrature points
        wRan = chaos_ran(cell.w, 2, uq)
        primRan = chaos_ran(cell.prim, 2, uq)

        # pure collocation settings
        #wRan =
        #    chaos_ran(cell.w, 2, uq) .+
        #    (chaos_ran(faceL.fw, 2, uq) .- chaos_ran(faceR.fw, 2, uq)) ./ cell.dx
        #primRan = uq_conserve_prim(wRan, KS.gas.γ, uq)

        # temperature protection
        if min(minimum(primRan[1, :]), minimum(primRan[end, :])) < 0
            @warn "negative temperature update"
            #wRan = chaos_ran(w_old, 2, uq)
            #primRan = chaos_ran(prim_old, 2, uq)

            @info "filtering"
            #=
            Lnorm = ones(uq.nr+1)
            for i in eachindex(Lnorm)
                Lnorm[i] = sum(uq.op.quad.weights .* 2 .* abs.(evaluate(i-1, uq.op.quad.nodes, uq.op)))
                #Lnorm[i] = exp(sum(uq.op.quad.weights .* abs.(evaluate(i-1, uq.op.quad.nodes, uq.op))))
            end
            =#
            iter = 1
            while min(minimum(primRan[1, :]), minimum(primRan[end, :])) < 0
                filter!(cell.w, 2, iter * dt, maximum(KS.ib.wL .- KS.ib.wR), uq.op, :l2)
                #filter!(cell.w, Lnorm, 1e-3, 2, :adapt)
                wRan = chaos_ran(cell.w, 2, uq)
                primRan = chaos_ran(cell.prim, 2, uq)

                iter += 1
                iter > 4 && break
            end
        end

        cell.w .= ran_chaos(wRan, 2, uq)
        cell.prim .= ran_chaos(primRan, 2, uq)

        #--- update particle distribution function ---#
        # flux -> f^{n+1}
        #@. cell.f += (faceL.ff - faceR.ff) / cell.dx

        fRan =
            chaos_ran(cell.f, 2, uq) .+
            (chaos_ran(faceL.ff, 2, uq) .- chaos_ran(faceR.ff, 2, uq)) ./ cell.dx

        # source -> f^{n+1}
        tau = uq_vhs_collision_time(primRan, KS.gas.μᵣ, KS.gas.ω, uq)

        gRan = zeros(KS.vSpace.nu, uq.op.quad.Nquad)
        for j in axes(gRan, 2)
            gRan[:, j] .= KitBase.maxwellian(KS.vSpace.u, primRan[:, j])
        end

        # BGK term
        for j in axes(fRan, 2)
            @. fRan[:, j] = (fRan[:, j] + dt / tau[j] * gRan[:, j]) / (1.0 + dt / tau[j])
        end

        cell.f .= ran_chaos(fRan, 2, uq)

        #--- record residuals ---#
        @. RES += (w_old[:, 1] - cell.w[:, 1])^2
        @. AVG += abs(cell.w[:, 1])

    elseif uq.method == "collocation"

        #--- update conservative flow variables: step 1 ---#
        # w^n
        w_old = deepcopy(cell.w)
        prim_old = deepcopy(cell.prim)

        # flux -> w^{n+1}
        @. cell.w += (faceL.fw - faceR.fw) / cell.dx
        cell.prim .= uq_conserve_prim(cell.w, KS.gas.γ, uq)

        # temperature protection
        if minimum(cell.prim[end, :]) < 0
            @warn "negative temperature update"
            cell.w .= w_old
            cell.prim .= prim_old
        end

        #--- update particle distribution function ---#
        # flux -> f^{n+1}
        @. cell.f += (faceL.ff - faceR.ff) / cell.dx

        # source -> f^{n+1}
        tau = uq_vhs_collision_time(cell.prim, KS.gas.μᵣ, KS.gas.ω, uq)

        g = zeros(KS.vSpace.nu, uq.op.quad.Nquad)
        for j in axes(g, 2)
            g[:, j] .= maxwellian(KS.vSpace.u, cell.prim[:, j])
        end

        # BGK term
        for j in axes(cell.f, 2)
            @. cell.f[:, j] = (cell.f[:, j] + dt / tau[j] * g[:, j]) / (1.0 + dt / tau[j])
        end

        #--- record residuals ---#
        @. RES += (w_old[:, 1] - cell.w[:, 1])^2
        @. AVG += abs(cell.w[:, 1])

    end

end

function step!(
    KS::SolverSet,
    uq::AbstractUQ,
    faceL::Interface1D4F,
    cell::ControlVolume1D4F,
    faceR::Interface1D4F,
    dt::AbstractFloat,
    RES::Array{<:AbstractFloat,2},
    AVG::Array{<:AbstractFloat,2},
    coll = :bgk::Symbol,
    isMHD = true::Bool,
) #// mixture 1D4F1V

    if uq.method == "galerkin"

        #--- update conservative flow variables: step 1 ---#
        # w^n
        w_old = deepcopy(cell.w)
        prim_old = deepcopy(cell.prim)

        # flux -> w^{n+1}
        #@. cell.w += (faceL.fw - faceR.fw) / cell.dx
        #cell.prim .= get_primitive(cell.w, KS.gas.γ, uq)

        # locate variables on random quadrature points
        #wRan = get_ran_array(cell.w, 2, uq)
        #primRan = get_ran_array(cell.prim, 2, uq)

        # pure collocation settings
        wRan =
            chaos_ran(cell.w, 2, uq) .+
            (chaos_ran(faceL.fw, 2, uq) .- chaos_ran(faceR.fw, 2, uq)) ./ cell.dx
        primRan = uq_conserve_prim(wRan, KS.gas.γ, uq)

        # temperature protection
        if min(minimum(primRan[5, :, 1]), minimum(primRan[5, :, 2])) < 0
            println("warning: temperature update is negative")
            wRan = chaos_ran(w_old, 2, uq)
            primRan = chaos_ran(prim_old, 2, uq)
        end

        #=
        # source -> w^{n+1}
        # DifferentialEquations.jl
        tau = get_tau(cell.prim, KS.gas.mi, KS.gas.ni, KS.gas.me, KS.gas.ne, KS.gas.Kn[1], uq)
        for j in axes(wRan, 2)
        prob = ODEProblem( mixture_source,
                    vcat(wRan[1:5,j,1], wRan[1:5,j,2]),
                    dt,
                    (tau[1], tau[2], KS.gas.mi, KS.gas.ni, KS.gas.me, KS.gas.ne, KS.gas.Kn[1], KS.gas.γ) )
        sol = solve(prob, Rosenbrock23())

        wRan[1:5,j,1] .= sol[end][1:5]
        wRan[1:5,j,2] .= sol[end][6:10]
        for k=1:2
        primRan[:,j,k] .= KitBase.conserve_prim(wRan[:,j,k], KS.gas.γ)
        end
        end

        # explicit
        tau = get_tau(cell.prim, KS.gas.mi, KS.gas.ni, KS.gas.me, KS.gas.ne, KS.gas.Kn[1], uq)
        mprim = get_mixprim(cell.prim, tau, KS.gas.mi, KS.gas.ni, KS.gas.me, KS.gas.ne, KS.gas.Kn[1], uq)
        mw = get_conserved(mprim, KS.gas.γ, uq)
        for k=1:2
        cell.w[:,:,k] .+= (mw[:,:,k] .- w_old[:,:,k]) .* dt ./ tau[k]
        end
        cell.prim .= get_primitive(cell.w, KS.gas.γ, uq)

        wRan .= get_ran_array(cell.w, 2, uq);
        primRan .= get_ran_array(cell.prim, 2, uq);
        =#

        #--- update electromagnetic variables ---#
        # flux -> E^{n+1} & B^{n+1}
        #@. cell.E[1,:] -= dt * (faceL.femR[1,:] + faceR.femL[1,:]) / cell.dx
        #@. cell.E[2,:] -= dt * (faceL.femR[2,:] + faceR.femL[2,:]) / cell.dx
        #@. cell.E[3,:] -= dt * (faceL.femR[3,:] + faceR.femL[3,:]) / cell.dx
        #@. cell.B[1,:] -= dt * (faceL.femR[4,:] + faceR.femL[4,:]) / cell.dx
        #@. cell.B[2,:] -= dt * (faceL.femR[5,:] + faceR.femL[5,:]) / cell.dx
        #@. cell.B[3,:] -= dt * (faceL.femR[6,:] + faceR.femL[6,:]) / cell.dx
        @. cell.ϕ -= dt * (faceL.femR[7, :] + faceR.femL[7, :]) / cell.dx
        @. cell.ψ -= dt * (faceL.femR[8, :] + faceR.femL[8, :]) / cell.dx

        ERan = chaos_ran(cell.E, 2, uq)
        ERan[1, :] .-=
            dt .* (
                evaluatePCE(faceL.femR[1, :], uq.op.quad.nodes, uq.op) .+
                evaluatePCE(faceR.femL[1, :], uq.op.quad.nodes, uq.op)
            ) ./ cell.dx
        ERan[2, :] .-=
            dt .* (
                evaluatePCE(faceL.femR[2, :], uq.op.quad.nodes, uq.op) .+
                evaluatePCE(faceR.femL[2, :], uq.op.quad.nodes, uq.op)
            ) ./ cell.dx
        ERan[3, :] .-=
            dt .* (
                evaluatePCE(faceL.femR[3, :], uq.op.quad.nodes, uq.op) .+
                evaluatePCE(faceR.femL[3, :], uq.op.quad.nodes, uq.op)
            ) ./ cell.dx
        BRan = chaos_ran(cell.B, 2, uq)
        BRan[1, :] .-=
            dt .* (
                evaluatePCE(faceL.femR[4, :], uq.op.quad.nodes, uq.op) .+
                evaluatePCE(faceR.femL[4, :], uq.op.quad.nodes, uq.op)
            ) ./ cell.dx
        BRan[2, :] .-=
            dt .* (
                evaluatePCE(faceL.femR[5, :], uq.op.quad.nodes, uq.op) .+
                evaluatePCE(faceR.femL[5, :], uq.op.quad.nodes, uq.op)
            ) ./ cell.dx
        BRan[3, :] .-=
            dt .* (
                evaluatePCE(faceL.femR[6, :], uq.op.quad.nodes, uq.op) .+
                evaluatePCE(faceR.femL[6, :], uq.op.quad.nodes, uq.op)
            ) ./ cell.dx

        # source -> ϕ
        #@. cell.ϕ += dt * (cell.w[1,:,1] / KS.gas.mi - cell.w[1,:,2] / KS.gas.me) / (KS.gas.lD^2 * KS.gas.rL)

        # source -> U^{n+1}, E^{n+1} and B^{n+1}
        mr = KS.gas.mi / KS.gas.me

        #ERan = get_ran_array(cell.E, 2, uq)
        #BRan = get_ran_array(cell.B, 2, uq)

        xRan = zeros(9, uq.op.quad.Nquad)
        for j in axes(xRan, 2)
            A, b = em_coefficients(
                primRan[:, j, :],
                ERan[:, j],
                BRan[:, j],
                mr,
                KS.gas.lD,
                KS.gas.rL,
                dt,
            )
            xRan[:, j] .= A \ b
        end

        lorenzRan = zeros(3, uq.op.quad.Nquad, 2)
        for j in axes(lorenzRan, 2)
            lorenzRan[1, j, 1] =
                0.5 * (
                    xRan[1, j] + ERan[1, j] + (primRan[3, j, 1] + xRan[5, j]) * BRan[3, j] -
                    (primRan[4, j, 1] + xRan[6, j]) * BRan[2, j]
                ) / KS.gas.rL
            lorenzRan[2, j, 1] =
                0.5 * (
                    xRan[2, j] + ERan[2, j] + (primRan[4, j, 1] + xRan[6, j]) * BRan[1, j] -
                    (primRan[2, j, 1] + xRan[4, j]) * BRan[3, j]
                ) / KS.gas.rL
            lorenzRan[3, j, 1] =
                0.5 * (
                    xRan[3, j] + ERan[3, j] + (primRan[2, j, 1] + xRan[4, j]) * BRan[2, j] -
                    (primRan[3, j, 1] + xRan[5, j]) * BRan[1, j]
                ) / KS.gas.rL
            lorenzRan[1, j, 2] =
                -0.5 *
                (
                    xRan[1, j] + ERan[1, j] + (primRan[3, j, 2] + xRan[8, j]) * BRan[3, j] -
                    (primRan[4, j, 2] + xRan[9, j]) * BRan[2, j]
                ) *
                mr / KS.gas.rL
            lorenzRan[2, j, 2] =
                -0.5 *
                (
                    xRan[2, j] + ERan[2, j] + (primRan[4, j, 2] + xRan[9, j]) * BRan[1, j] -
                    (primRan[2, j, 2] + xRan[7, j]) * BRan[3, j]
                ) *
                mr / KS.gas.rL
            lorenzRan[3, j, 2] =
                -0.5 *
                (
                    xRan[3, j] + ERan[3, j] + (primRan[2, j, 2] + xRan[7, j]) * BRan[2, j] -
                    (primRan[3, j, 2] + xRan[8, j]) * BRan[1, j]
                ) *
                mr / KS.gas.rL
        end

        ERan[1, :] .= xRan[1, :]
        ERan[2, :] .= xRan[2, :]
        ERan[3, :] .= xRan[3, :]

        #--- update conservative flow variables: step 2 ---#
        primRan[2, :, 1] .= xRan[4, :]
        primRan[3, :, 1] .= xRan[5, :]
        primRan[4, :, 1] .= xRan[6, :]
        primRan[2, :, 2] .= xRan[7, :]
        primRan[3, :, 2] .= xRan[8, :]
        primRan[4, :, 2] .= xRan[9, :]

        for j in axes(wRan, 2)
            wRan[:, j, :] .= KitBase.mixture_prim_conserve(primRan[:, j, :], KS.gas.γ)
        end

        cell.w .= ran_chaos(wRan, 2, uq)
        cell.prim .= ran_chaos(primRan, 2, uq)
        cell.E .= ran_chaos(ERan, 2, uq)
        cell.lorenz .= ran_chaos(lorenzRan, 2, uq)
        cell.B .= ran_chaos(BRan, 2, uq)

        #--- update particle distribution function ---#
        # flux -> f^{n+1}
        #@. cell.h0 += (faceL.fh0 - faceR.fh0) / cell.dx
        #@. cell.h1 += (faceL.fh1 - faceR.fh1) / cell.dx
        #@. cell.h2 += (faceL.fh2 - faceR.fh2) / cell.dx
        #@. cell.h3 += (faceL.fh3 - faceR.fh3) / cell.dx

        #h0Ran = get_ran_array(cell.h0, 2, uq)
        #h1Ran = get_ran_array(cell.h1, 2, uq)
        #h2Ran = get_ran_array(cell.h2, 2, uq)
        #h3Ran = get_ran_array(cell.h3, 2, uq)

        h0Ran =
            chaos_ran(cell.h0, 2, uq) .+
            (chaos_ran(faceL.fh0, 2, uq) .- chaos_ran(faceR.fh0, 2, uq)) ./ cell.dx
        h1Ran =
            chaos_ran(cell.h1, 2, uq) .+
            (chaos_ran(faceL.fh1, 2, uq) .- chaos_ran(faceR.fh1, 2, uq)) ./ cell.dx
        h2Ran =
            chaos_ran(cell.h2, 2, uq) .+
            (chaos_ran(faceL.fh2, 2, uq) .- chaos_ran(faceR.fh2, 2, uq)) ./ cell.dx
        h3Ran =
            chaos_ran(cell.h3, 2, uq) .+
            (chaos_ran(faceL.fh3, 2, uq) .- chaos_ran(faceR.fh3, 2, uq)) ./ cell.dx

        # force -> f^{n+1} : step 1
        for j in axes(h0Ran, 2)
            _h0 = @view h0Ran[:, j, :]
            _h1 = @view h1Ran[:, j, :]
            _h2 = @view h2Ran[:, j, :]
            _h3 = @view h3Ran[:, j, :]

            shift_pdf!(_h0, lorenzRan[1, j, :], KS.vSpace.du[1, :], dt)
            shift_pdf!(_h1, lorenzRan[1, j, :], KS.vSpace.du[1, :], dt)
            shift_pdf!(_h2, lorenzRan[1, j, :], KS.vSpace.du[1, :], dt)
            shift_pdf!(_h3, lorenzRan[1, j, :], KS.vSpace.du[1, :], dt)
        end

        # force -> f^{n+1} : step 2
        for k in axes(h1Ran, 3), j in axes(h1Ran, 2)
            @. h3Ran[:, j, k] +=
                2.0 * dt * lorenzRan[2, j, k] * h1Ran[:, j, k] +
                (dt * lorenzRan[2, j, k])^2 * h0Ran[:, j, k] +
                2.0 * dt * lorenzRan[3, j, k] * h2Ran[:, j, k] +
                (dt * lorenzRan[3, j, k])^2 * h0Ran[:, j, k]
            @. h2Ran[:, j, k] += dt * lorenzRan[3, j, k] * h0Ran[:, j, k]
            @. h1Ran[:, j, k] += dt * lorenzRan[2, j, k] * h0Ran[:, j, k]
        end

        # source -> f^{n+1}
        #tau = get_tau(cell.prim, KS.gas.mi, KS.gas.ni, KS.gas.me, KS.gas.ne, KS.gas.Kn[1], uq)
        tau = get_tau(primRan, KS.gas.mi, KS.gas.ni, KS.gas.me, KS.gas.ne, KS.gas.Kn[1])

        # interspecies interaction
        for j in axes(primRan, 2)
            primRan[:, j, :] .= KitBase.aap_hs_prim(
                primRan[:, j, :],
                tau,
                KS.gas.mi,
                KS.gas.ni,
                KS.gas.me,
                KS.gas.ne,
                KS.gas.Kn[1],
            )
        end

        gRan = zeros(KS.vSpace.nu, uq.op.quad.Nquad, 2)
        for k in axes(gRan, 3)
            for j in axes(gRan, 2)
                gRan[:, j, k] .= KitBase.maxwellian(KS.vSpace.u[:, k], primRan[:, j, k])
            end
        end

        # BGK term
        for j in axes(h0Ran, 2)
            Mu, Mv, Mw, MuL, MuR = KitBase.mixture_gauss_moments(primRan[:, j, :], KS.gas.K)
            for k in axes(h0Ran, 3)
                @. h0Ran[:, j, k] =
                    (h0Ran[:, j, k] + dt / tau[k] * gRan[:, j, k]) / (1.0 + dt / tau[k])
                @. h1Ran[:, j, k] =
                    (h1Ran[:, j, k] + dt / tau[k] * Mv[1, k] * gRan[:, j, k]) /
                    (1.0 + dt / tau[k])
                @. h2Ran[:, j, k] =
                    (h2Ran[:, j, k] + dt / tau[k] * Mw[1, k] * gRan[:, j, k]) /
                    (1.0 + dt / tau[k])
                @. h3Ran[:, j, k] =
                    (h3Ran[:, j, k] + dt / tau[k] * (Mv[2, k] + Mw[2, k]) * gRan[:, j, k]) /
                    (1.0 + dt / tau[k])
            end
        end

        cell.h0 .= ran_chaos(h0Ran, 2, uq)
        cell.h1 .= ran_chaos(h1Ran, 2, uq)
        cell.h2 .= ran_chaos(h2Ran, 2, uq)
        cell.h3 .= ran_chaos(h3Ran, 2, uq)

        #--- record residuals ---#
        @. RES += (w_old[:, 1, :] - cell.w[:, 1, :])^2
        @. AVG += abs(cell.w[:, 1, :])
        #=
        #-- filter ---#
        λ = 0.00001
        for k in 1:2, i in 1:5
        filter!(cell.w[i,:,k], λ)
        filter!(cell.prim[i,:,k], λ)
        end
        for k in 1:2, i in axes(cell.h0, 1)
        filter!(cell.h0[i,:,k], λ)
        filter!(cell.h1[i,:,k], λ)
        filter!(cell.h2[i,:,k], λ)
        filter!(cell.h3[i,:,k], λ)
        end
        for k in 1:2, i in 1:3
        filter!(cell.lorenz[i,:,k], λ)
        end
        for i in 1:3
        filter!(cell.E[i,:], λ)
        filter!(cell.B[i,:], λ)
        end
        filter!(cell.ϕ, λ)
        filter!(cell.ψ, λ)
        =#

    elseif uq.method == "collocation"

        #--- update conservative flow variables: step 1 ---#
        # w^n
        w_old = deepcopy(cell.w)
        prim_old = deepcopy(cell.prim)

        # flux -> w^{n+1}
        @. cell.w += (faceL.fw - faceR.fw) / cell.dx
        for j in axes(cell.prim, 2)
            cell.prim[:, j, :] .= mixture_conserve_prim(cell.w[:, j, :], KS.gas.γ)
        end

        # temperature protection
        if min(minimum(cell.prim[5, :, 1]), minimum(cell.prim[5, :, 2])) < 0
            @warn "negative temperature update"
            cell.w .= w_old
            cell.prim .= prim_old
        end

        #=
        # source -> w^{n+1}
        # DifferentialEquations.jl
        tau = get_tau(cell.prim, KS.gas.mi, KS.gas.ni, KS.gas.me, KS.gas.ne, KS.gas.Kn[1])
        for j in axes(wRan, 2)
        prob = ODEProblem( mixture_source,
            vcat(cell.w[1:5,j,1], cell.w[1:5,j,2]),
            dt,
            (tau[1], tau[2], KS.gas.mi, KS.gas.ni, KS.gas.me, KS.gas.ne, KS.gas.Kn[1], KS.gas.γ) )
        sol = solve(prob, Rosenbrock23())

        cell.w[1:5,j,1] .= sol[end][1:5]
        cell.w[1:5,j,2] .= sol[end][6:10]
        for k=1:2
        cell.prim[:,j,k] .= KitBase.conserve_prim(cell.w[:,j,k], KS.gas.γ)
        end
        end
        =#
        #=
        # explicit
        tau = get_tau(cell.prim, KS.gas.mi, KS.gas.ni, KS.gas.me, KS.gas.ne, KS.gas.Kn[1])
        mprim = get_mixprim(cell.prim, tau, KS.gas.mi, KS.gas.ni, KS.gas.me, KS.gas.ne, KS.gas.Kn[1])
        mw = get_conserved(mprim, KS.gas.γ)
        for k=1:2
        cell.w[:,:,k] .+= (mw[:,:,k] .- cell.w[:,:,k]) .* dt ./ tau[k]
        end
        cell.prim .= get_primitive(cell.w, KS.gas.γ)
        =#

        #--- update electromagnetic variables ---#
        # flux -> E^{n+1} & B^{n+1}
        @. cell.E[1, :] -= dt * (faceL.femR[1, :] + faceR.femL[1, :]) / cell.dx
        @. cell.E[2, :] -= dt * (faceL.femR[2, :] + faceR.femL[2, :]) / cell.dx
        @. cell.E[3, :] -= dt * (faceL.femR[3, :] + faceR.femL[3, :]) / cell.dx
        @. cell.B[1, :] -= dt * (faceL.femR[4, :] + faceR.femL[4, :]) / cell.dx
        @. cell.B[2, :] -= dt * (faceL.femR[5, :] + faceR.femL[5, :]) / cell.dx
        @. cell.B[3, :] -= dt * (faceL.femR[6, :] + faceR.femL[6, :]) / cell.dx
        @. cell.ϕ -= dt * (faceL.femR[7, :] + faceR.femL[7, :]) / cell.dx
        @. cell.ψ -= dt * (faceL.femR[8, :] + faceR.femL[8, :]) / cell.dx

        # source -> ϕ
        #@. cell.ϕ += dt * (cell.w[1,:,1] / KS.gas.mi - cell.w[1,:,2] / KS.gas.me) / (KS.gas.lD^2 * KS.gas.rL)

        # source -> U^{n+1}, E^{n+1} and B^{n+1}
        mr = KS.gas.mi / KS.gas.me

        x = zeros(9, uq.op.quad.Nquad)
        for j in axes(x, 2)
            A, b = em_coefficients(
                cell.prim[:, j, :],
                cell.E[:, j],
                cell.B[:, j],
                mr,
                KS.gas.lD,
                KS.gas.rL,
                dt,
            )
            x[:, j] .= A \ b
        end

        #--- calculate lorenz force ---#
        for j in axes(cell.lorenz, 2)
            cell.lorenz[1, j, 1] =
                0.5 * (
                    x[1, j] + cell.E[1, j] + (cell.prim[3, j, 1] + x[5, j]) * cell.B[3, j] -
                    (cell.prim[4, j, 1] + x[6, j]) * cell.B[2, j]
                ) / KS.gas.rL
            cell.lorenz[2, j, 1] =
                0.5 * (
                    x[2, j] + cell.E[2, j] + (cell.prim[4, j, 1] + x[6, j]) * cell.B[1, j] -
                    (cell.prim[2, j, 1] + x[4, j]) * cell.B[3, j]
                ) / KS.gas.rL
            cell.lorenz[3, j, 1] =
                0.5 * (
                    x[3, j] + cell.E[3, j] + (cell.prim[2, j, 1] + x[4, j]) * cell.B[2, j] -
                    (cell.prim[3, j, 1] + x[5, j]) * cell.B[1, j]
                ) / KS.gas.rL
            cell.lorenz[1, j, 2] =
                -0.5 *
                (
                    x[1, j] + cell.E[1, j] + (cell.prim[3, j, 2] + x[8, j]) * cell.B[3, j] -
                    (cell.prim[4, j, 2] + x[9, j]) * cell.B[2, j]
                ) *
                mr / KS.gas.rL
            cell.lorenz[2, j, 2] =
                -0.5 *
                (
                    x[2, j] + cell.E[2, j] + (cell.prim[4, j, 2] + x[9, j]) * cell.B[1, j] -
                    (cell.prim[2, j, 2] + x[7, j]) * cell.B[3, j]
                ) *
                mr / KS.gas.rL
            cell.lorenz[3, j, 2] =
                -0.5 *
                (
                    x[3, j] + cell.E[3, j] + (cell.prim[2, j, 2] + x[7, j]) * cell.B[2, j] -
                    (cell.prim[3, j, 2] + x[8, j]) * cell.B[1, j]
                ) *
                mr / KS.gas.rL
        end

        cell.E[1, :] .= x[1, :]
        cell.E[2, :] .= x[2, :]
        cell.E[3, :] .= x[3, :]

        #--- update conservative flow variables: step 2 ---#
        cell.prim[2, :, 1] .= x[4, :]
        cell.prim[3, :, 1] .= x[5, :]
        cell.prim[4, :, 1] .= x[6, :]
        cell.prim[2, :, 2] .= x[7, :]
        cell.prim[3, :, 2] .= x[8, :]
        cell.prim[4, :, 2] .= x[9, :]

        for j in axes(cell.w, 2)
            cell.w[:, j, :] .= KitBase.mixture_prim_conserve(cell.prim[:, j, :], KS.gas.γ)
        end

        #--- update particle distribution function ---#
        # flux -> f^{n+1}
        @. cell.h0 += (faceL.fh0 - faceR.fh0) / cell.dx
        @. cell.h1 += (faceL.fh1 - faceR.fh1) / cell.dx
        @. cell.h2 += (faceL.fh2 - faceR.fh2) / cell.dx
        @. cell.h3 += (faceL.fh3 - faceR.fh3) / cell.dx

        # force -> f^{n+1} : step 1
        for k in axes(cell.h0, 3)
            for j in axes(cell.h0, 2)
                _h0 = @view cell.h0[:, j, k]
                _h1 = @view cell.h1[:, j, k]
                _h2 = @view cell.h2[:, j, k]
                _h3 = @view cell.h3[:, j, k]

                shift_pdf!(_h0, cell.lorenz[1, j, k], KS.vSpace.du[1, k], dt)
                shift_pdf!(_h1, cell.lorenz[1, j, k], KS.vSpace.du[1, k], dt)
                shift_pdf!(_h2, cell.lorenz[1, j, k], KS.vSpace.du[1, k], dt)
                shift_pdf!(_h3, cell.lorenz[1, j, k], KS.vSpace.du[1, k], dt)
            end
        end

        # force -> f^{n+1} : step 2
        for k in axes(cell.h1, 3), j in axes(cell.h1, 2)
            @. cell.h3[:, j, k] +=
                2.0 * dt * cell.lorenz[2, j, k] * cell.h1[:, j, k] +
                (dt * cell.lorenz[2, j, k])^2 * cell.h0[:, j, k] +
                2.0 * dt * cell.lorenz[3, j, k] * cell.h2[:, j, k] +
                (dt * cell.lorenz[3, j, k])^2 * cell.h0[:, j, k]
            @. cell.h2[:, j, k] += dt * cell.lorenz[3, j, k] * cell.h0[:, j, k]
            @. cell.h1[:, j, k] += dt * cell.lorenz[2, j, k] * cell.h0[:, j, k]
        end

        # source -> f^{n+1}
        tau = uq_aap_hs_collision_time(
            cell.prim,
            KS.gas.mi,
            KS.gas.ni,
            KS.gas.me,
            KS.gas.ne,
            KS.gas.Kn[1],
            uq,
        )

        # interspecies interaction
        prim = deepcopy(cell.prim)
        #for j in axes(prim, 2)
        #    prim[:, j, :] .= KitBase.aap_hs_prim(
        #        cell.prim[:, j, :],
        #        tau,
        #        KS.gas.mi,
        #        KS.gas.ni,
        #        KS.gas.me,
        #        KS.gas.ne,
        #        KS.gas.Kn[1],
        #    )
        #end

        g = zeros(KS.vSpace.nu, uq.op.quad.Nquad, 2)
        for j in axes(g, 2)
            g[:, j, :] .= mixture_maxwellian(KS.vSpace.u, prim[:, j, :])
        end

        # BGK term
        for j in axes(cell.h0, 2)
            Mu, Mv, Mw, MuL, MuR = mixture_gauss_moments(prim[:, j, :], KS.gas.K)
            for k in axes(cell.h0, 3)
                @. cell.h0[:, j, k] =
                    (cell.h0[:, j, k] + dt / tau[k] * g[:, j, k]) / (1.0 + dt / tau[k])
                @. cell.h1[:, j, k] =
                    (cell.h1[:, j, k] + dt / tau[k] * Mv[1, k] * g[:, j, k]) /
                    (1.0 + dt / tau[k])
                @. cell.h2[:, j, k] =
                    (cell.h2[:, j, k] + dt / tau[k] * Mw[1, k] * g[:, j, k]) /
                    (1.0 + dt / tau[k])
                @. cell.h3[:, j, k] =
                    (cell.h3[:, j, k] + dt / tau[k] * (Mv[2, k] + Mw[2, k]) * g[:, j, k]) /
                    (1.0 + dt / tau[k])
            end
        end

        #--- record residuals ---#
        @. RES += (w_old[:, 1, :] - cell.w[:, 1, :])^2
        @. AVG += abs(cell.w[:, 1, :])

    end

end

#--- mixture 1D3F ---#
function step!(
    KS::SolverSet,
    uq::AbstractUQ,
    faceL::Interface1D3F,
    cell::ControlVolume1D3F,
    faceR::Interface1D3F,
    dt::AbstractFloat,
    RES::Array{<:AbstractFloat,2},
    AVG::Array{<:AbstractFloat,2},
    coll = :bgk::Symbol,
    isMHD = true::Bool,
)

    if uq.method == "galerkin"

        @assert size(cell.w, 2) == uq.nr + 1

        #--- update conservative flow variables: step 1 ---#
        # w^n
        w_old = deepcopy(cell.w)
        prim_old = deepcopy(cell.prim)

        # flux -> w^{n+1}
        @. cell.w += (faceL.fw - faceR.fw) / cell.dx
        cell.prim .= get_primitive(cell.w, KS.gas.γ, uq)

        # locate variables on random quadrature points
        wRan = chaos_ran(cell.w, 2, uq)
        primRan = chaos_ran(cell.prim, 2, uq)

        # temperature protection
        if min(minimum(primRan[5, :, 1]), minimum(primRan[5, :, 2])) < 0
            @warn "negative temperature update"
            wRan = chaos_ran(w_old, 2, uq)
            primRan = chaos_ran(prim_old, 2, uq)
        end

        if isMHD == false
            tauRan = uq_aap_hs_collision_time(
                primRan,
                KS.gas.mi,
                KS.gas.ni,
                KS.gas.me,
                KS.gas.ne,
                KS.gas.Kn[1],
                uq,
            )
            mprimRan = uq_aap_hs_prim(
                primRan,
                tauRan,
                KS.gas.mi,
                KS.gas.ni,
                KS.gas.me,
                KS.gas.ne,
                KS.gas.Kn[1],
                uq,
            )
            mwRan = uq_prim_conserve(mprimRan, KS.gas.γ, uq)
            for k in axes(wRan, 3)
                @. wRan[:, :, k] += (mwRan[:, :, k] - wRan[:, :, k]) * dt / tauRan[k]
            end
            primRan .= uq_conserve_prim(wRan, KS.gas.γ, uq)
        end

        #--- update electromagnetic variables ---#
        # flux -> E^{n+1} & B^{n+1}
        @. cell.E[1, :] -= dt * (faceL.femR[1, :] + faceR.femL[1, :]) / cell.dx
        @. cell.E[2, :] -= dt * (faceL.femR[2, :] + faceR.femL[2, :]) / cell.dx
        @. cell.E[3, :] -= dt * (faceL.femR[3, :] + faceR.femL[3, :]) / cell.dx
        @. cell.B[1, :] -= dt * (faceL.femR[4, :] + faceR.femL[4, :]) / cell.dx
        @. cell.B[2, :] -= dt * (faceL.femR[5, :] + faceR.femL[5, :]) / cell.dx
        @. cell.B[3, :] -= dt * (faceL.femR[6, :] + faceR.femL[6, :]) / cell.dx
        @. cell.ϕ -= dt * (faceL.femR[7, :] + faceR.femL[7, :]) / cell.dx
        @. cell.ψ -= dt * (faceL.femR[8, :] + faceR.femL[8, :]) / cell.dx

        ERan = chaos_ran(cell.E, 2, uq)
        BRan = chaos_ran(cell.B, 2, uq)

        # source -> ϕ
        #@. cell.ϕ += dt * (cell.w[1,:,1] / KS.gas.mi - cell.w[1,:,2] / KS.gas.me) / (KS.gas.lD^2 * KS.gas.rL)

        # source -> U^{n+1}, E^{n+1} and B^{n+1}
        mr = KS.gas.mi / KS.gas.me

        xRan = zeros(9, uq.op.quad.Nquad)
        for j in axes(xRan, 2)
            A, b = em_coefficients(
                primRan[:, j, :],
                ERan[:, j],
                BRan[:, j],
                mr,
                KS.gas.lD,
                KS.gas.rL,
                dt,
            )
            xRan[:, j] .= A \ b
        end

        lorenzRan = zeros(3, uq.op.quad.Nquad, 2)
        for j in axes(lorenzRan, 2)
            lorenzRan[1, j, 1] =
                0.5 * (
                    xRan[1, j] + ERan[1, j] + (primRan[3, j, 1] + xRan[5, j]) * BRan[3, j] -
                    (primRan[4, j, 1] + xRan[6, j]) * BRan[2, j]
                ) / KS.gas.rL
            lorenzRan[2, j, 1] =
                0.5 * (
                    xRan[2, j] + ERan[2, j] + (primRan[4, j, 1] + xRan[6, j]) * BRan[1, j] -
                    (primRan[2, j, 1] + xRan[4, j]) * BRan[3, j]
                ) / KS.gas.rL
            lorenzRan[3, j, 1] =
                0.5 * (
                    xRan[3, j] + ERan[3, j] + (primRan[2, j, 1] + xRan[4, j]) * BRan[2, j] -
                    (primRan[3, j, 1] + xRan[5, j]) * BRan[1, j]
                ) / KS.gas.rL
            lorenzRan[1, j, 2] =
                -0.5 *
                (
                    xRan[1, j] + ERan[1, j] + (primRan[3, j, 2] + xRan[8, j]) * BRan[3, j] -
                    (primRan[4, j, 2] + xRan[9, j]) * BRan[2, j]
                ) *
                mr / KS.gas.rL
            lorenzRan[2, j, 2] =
                -0.5 *
                (
                    xRan[2, j] + ERan[2, j] + (primRan[4, j, 2] + xRan[9, j]) * BRan[1, j] -
                    (primRan[2, j, 2] + xRan[7, j]) * BRan[3, j]
                ) *
                mr / KS.gas.rL
            lorenzRan[3, j, 2] =
                -0.5 *
                (
                    xRan[3, j] + ERan[3, j] + (primRan[2, j, 2] + xRan[7, j]) * BRan[2, j] -
                    (primRan[3, j, 2] + xRan[8, j]) * BRan[1, j]
                ) *
                mr / KS.gas.rL
        end

        ERan[1, :] .= xRan[1, :]
        ERan[2, :] .= xRan[2, :]
        ERan[3, :] .= xRan[3, :]

        #--- update conservative flow variables: step 2 ---#
        primRan[2, :, 1] .= xRan[4, :]
        primRan[3, :, 1] .= xRan[5, :]
        primRan[4, :, 1] .= xRan[6, :]
        primRan[2, :, 2] .= xRan[7, :]
        primRan[3, :, 2] .= xRan[8, :]
        primRan[4, :, 2] .= xRan[9, :]

        for j in axes(wRan, 2)
            wRan[:, j, :] .= KitBase.mixture_prim_conserve(primRan[:, j, :], KS.gas.γ)
        end

        cell.w .= ran_chaos(wRan, 2, uq)
        cell.prim .= ran_chaos(primRan, 2, uq)
        cell.E .= ran_chaos(ERan, 2, uq)
        cell.lorenz .= ran_chaos(lorenzRan, 2, uq)
        cell.B .= ran_chaos(BRan, 2, uq)

        #--- update particle distribution function ---#
        # flux -> f^{n+1}
        @. cell.h0 += (faceL.fh0 - faceR.fh0) / cell.dx
        @. cell.h1 += (faceL.fh1 - faceR.fh1) / cell.dx
        @. cell.h2 += (faceL.fh2 - faceR.fh2) / cell.dx

        h0Ran = chaos_ran(cell.h0, 3, uq)
        h1Ran = chaos_ran(cell.h1, 3, uq)
        h2Ran = chaos_ran(cell.h2, 3, uq)

        # force -> f^{n+1} : step 1
        for k in axes(h0Ran, 4)
            for j in axes(h0Ran, 3)
                for i in axes(h0Ran, 2)
                    _h0 = @view h0Ran[:, i, j, k]
                    _h1 = @view h1Ran[:, i, j, k]
                    _h2 = @view h2Ran[:, i, j, k]

                    shift_pdf!(_h0, cell.lorenz[1, j, k], KS.vSpace.du[1, i, k], dt)
                    shift_pdf!(_h1, cell.lorenz[1, j, k], KS.vSpace.du[1, i, k], dt)
                    shift_pdf!(_h2, cell.lorenz[1, j, k], KS.vSpace.du[1, i, k], dt)
                end
            end
        end

        for k in axes(h0Ran, 4)
            for j in axes(h0Ran, 3)
                for i in axes(h0Ran, 1)
                    _h0 = @view h0Ran[i, :, j, k]
                    _h1 = @view h1Ran[i, :, j, k]
                    _h2 = @view h2Ran[i, :, j, k]

                    shift_pdf!(_h0, lorenzRan[2, j, k], KS.vSpace.dv[i, 1, k], dt)
                    shift_pdf!(_h1, lorenzRan[2, j, k], KS.vSpace.dv[i, 1, k], dt)
                    shift_pdf!(_h2, lorenzRan[2, j, k], KS.vSpace.dv[i, 1, k], dt)
                end
            end
        end

        # force -> f^{n+1} : step 2
        for k in axes(h1Ran, 4), j in axes(h1Ran, 3)
            @. h2Ran[:, :, j, k] +=
                2.0 * dt * lorenzRan[3, j, k] * h1Ran[:, :, j, k] +
                (dt * lorenzRan[3, j, k])^2 * h0Ran[:, :, j, k]
            @. h1Ran[:, :, j, k] += dt * lorenzRan[3, j, k] * h0Ran[:, :, j, k]
        end

        # source -> f^{n+1}
        tau = uq_aap_hs_collision_time(
            primRan,
            KS.gas.mi,
            KS.gas.ni,
            KS.gas.me,
            KS.gas.ne,
            KS.gas.Kn[1],
            uq,
        )

        if isMHD == false
            # interspecies interaction
            for j in axes(primRan, 2)
                primRan[:, j, :] .= KitBase.aap_hs_prim(
                    primRan[:, j, :],
                    tau,
                    KS.gas.mi,
                    KS.gas.ni,
                    KS.gas.me,
                    KS.gas.ne,
                    KS.gas.Kn[1],
                )
            end
        end

        H0Ran = zeros(KS.vSpace.nu, KS.vSpace.nv, uq.op.quad.Nquad, 2)
        H1Ran = similar(H0Ran)
        H2Ran = similar(H0Ran)
        for k in axes(H0Ran, 4), j in axes(H0Ran, 3)
            H0Ran[:, :, j, k] .=
                maxwellian(KS.vSpace.u[:, :, k], KS.vSpace.v[:, :, k], primRan[:, j, k])
            @. H1Ran[:, :, j, k] = H0Ran[:, :, j, k] * primRan[4, j, k]
            @. H2Ran[:, :, j, k] =
                H0Ran[:, :, j, k] * (primRan[4, j, k]^2 + 1.0 / (2.0 * primRan[5, j, k]))
        end

        # BGK term
        for k in axes(h0Ran, 4)
            for j in axes(h0Ran, 3)
                @. h0Ran[:, :, j, k] =
                    (h0Ran[:, :, j, k] + dt / tau[k] * H0Ran[:, :, j, k]) /
                    (1.0 + dt / tau[k])
                @. h1Ran[:, :, j, k] =
                    (h1Ran[:, :, j, k] + dt / tau[k] * H1Ran[:, :, j, k]) /
                    (1.0 + dt / tau[k])
                @. h2Ran[:, j, :, k] =
                    (h2Ran[:, :, j, k] + dt / tau[k] * H2Ran[:, :, j, k]) /
                    (1.0 + dt / tau[k])
            end
        end

        cell.h0 .= ran_chaos(h0Ran, 3, uq)
        cell.h1 .= ran_chaos(h1Ran, 3, uq)
        cell.h2 .= ran_chaos(h2Ran, 3, uq)

        #--- record residuals ---#
        @. RES += (w_old[:, 1, :] - cell.w[:, 1, :])^2
        @. AVG += abs(cell.w[:, 1, :])

    elseif uq.method == "collocation"

        w_old = deepcopy(cell.w)
        prim_old = deepcopy(cell.prim)

        # flux -> w^{n+1}
        @. cell.w += (faceL.fw - faceR.fw) / cell.dx
        cell.prim = uq_conserve_prim(cell.w, KS.gas.γ, uq)

        # temperature protection
        if min(minimum(cell.prim[5, :, 1]), minimum(cell.prim[5, :, 2])) < 0
            @warn "negative temperature update"
            cell.w .= w_old
            cell.prim .= prim_old
        end

        if isMHD == false
            # explicit
            tau = uq_aap_hs_collision_time(
                cell.prim,
                KS.gas.mi,
                KS.gas.ni,
                KS.gas.me,
                KS.gas.ne,
                KS.gas.Kn[1],
                uq,
            )
            mprim = uq_aap_hs_prim(
                cell.prim,
                tau,
                KS.gas.mi,
                KS.gas.ni,
                KS.gas.me,
                KS.gas.ne,
                KS.gas.Kn[1],
                uq,
            )
            mw = uq_prim_conserve(mprim, KS.gas.γ, uq)
            for k = 1:2
                @. cell.w[:, :, k] += (mw[:, :, k] - cell.w[:, :, k]) * dt / tau[k]
            end
            cell.prim .= uq_conserve_prim(cell.w, KS.gas.γ, uq)
        end

        #--- update electromagnetic variables ---#
        # flux -> E^{n+1} & B^{n+1}
        @. cell.E[1, :] -= dt * (faceL.femR[1, :] + faceR.femL[1, :]) / cell.dx
        @. cell.E[2, :] -= dt * (faceL.femR[2, :] + faceR.femL[2, :]) / cell.dx
        @. cell.E[3, :] -= dt * (faceL.femR[3, :] + faceR.femL[3, :]) / cell.dx
        @. cell.B[1, :] -= dt * (faceL.femR[4, :] + faceR.femL[4, :]) / cell.dx
        @. cell.B[2, :] -= dt * (faceL.femR[5, :] + faceR.femL[5, :]) / cell.dx
        @. cell.B[3, :] -= dt * (faceL.femR[6, :] + faceR.femL[6, :]) / cell.dx
        @. cell.ϕ -= dt * (faceL.femR[7, :] + faceR.femL[7, :]) / cell.dx
        @. cell.ψ -= dt * (faceL.femR[8, :] + faceR.femL[8, :]) / cell.dx

        # source -> ϕ
        #@. cell.ϕ += dt * (cell.w[1,:,1] / KS.gas.mi - cell.w[1,:,2] / KS.gas.me) / (KS.gas.lD^2 * KS.gas.rL)

        # source -> U^{n+1}, E^{n+1} and B^{n+1}
        mr = KS.gas.mi / KS.gas.me

        x = zeros(9, uq.op.quad.Nquad)
        for j in axes(x, 2)
            A, b = em_coefficients(
                cell.prim[:, j, :],
                cell.E[:, j],
                cell.B[:, j],
                mr,
                KS.gas.lD,
                KS.gas.rL,
                dt,
            )
            x[:, j] .= A \ b
        end

        #--- calculate lorenz force ---#
        for j in axes(cell.lorenz, 2)
            cell.lorenz[1, j, 1] =
                0.5 * (
                    x[1, j] + cell.E[1, j] + (cell.prim[3, j, 1] + x[5, j]) * cell.B[3, j] -
                    (cell.prim[4, j, 1] + x[6, j]) * cell.B[2, j]
                ) / KS.gas.rL
            cell.lorenz[2, j, 1] =
                0.5 * (
                    x[2, j] + cell.E[2, j] + (cell.prim[4, j, 1] + x[6, j]) * cell.B[1, j] -
                    (cell.prim[2, j, 1] + x[4, j]) * cell.B[3, j]
                ) / KS.gas.rL
            cell.lorenz[3, j, 1] =
                0.5 * (
                    x[3, j] + cell.E[3, j] + (cell.prim[2, j, 1] + x[4, j]) * cell.B[2, j] -
                    (cell.prim[3, j, 1] + x[5, j]) * cell.B[1, j]
                ) / KS.gas.rL
            cell.lorenz[1, j, 2] =
                -0.5 *
                (
                    x[1, j] + cell.E[1, j] + (cell.prim[3, j, 2] + x[8, j]) * cell.B[3, j] -
                    (cell.prim[4, j, 2] + x[9, j]) * cell.B[2, j]
                ) *
                mr / KS.gas.rL
            cell.lorenz[2, j, 2] =
                -0.5 *
                (
                    x[2, j] + cell.E[2, j] + (cell.prim[4, j, 2] + x[9, j]) * cell.B[1, j] -
                    (cell.prim[2, j, 2] + x[7, j]) * cell.B[3, j]
                ) *
                mr / KS.gas.rL
            cell.lorenz[3, j, 2] =
                -0.5 *
                (
                    x[3, j] + cell.E[3, j] + (cell.prim[2, j, 2] + x[7, j]) * cell.B[2, j] -
                    (cell.prim[3, j, 2] + x[8, j]) * cell.B[1, j]
                ) *
                mr / KS.gas.rL
        end

        cell.E[1, :] .= x[1, :]
        cell.E[2, :] .= x[2, :]
        cell.E[3, :] .= x[3, :]

        #--- update conservative flow variables: step 2 ---#
        cell.prim[2, :, 1] .= x[4, :]
        cell.prim[3, :, 1] .= x[5, :]
        cell.prim[4, :, 1] .= x[6, :]
        cell.prim[2, :, 2] .= x[7, :]
        cell.prim[3, :, 2] .= x[8, :]
        cell.prim[4, :, 2] .= x[9, :]

        for j in axes(cell.w, 2)
            cell.w[:, j, :] .= KitBase.mixture_prim_conserve(cell.prim[:, j, :], KS.gas.γ)
        end

        #--- update particle distribution function ---#
        # flux -> f^{n+1}
        @. cell.h0 += (faceL.fh0 - faceR.fh0) / cell.dx
        @. cell.h1 += (faceL.fh1 - faceR.fh1) / cell.dx
        @. cell.h2 += (faceL.fh2 - faceR.fh2) / cell.dx

        # force -> f^{n+1} : step 1
        for k in axes(cell.h0, 4)
            for j in axes(cell.h0, 3)
                for i in axes(cell.h0, 2)
                    _h0 = @view cell.h0[:, i, j, k]
                    _h1 = @view cell.h1[:, i, j, k]
                    _h2 = @view cell.h2[:, i, j, k]

                    shift_pdf!(_h0, cell.lorenz[1, j, k], KS.vSpace.du[1, i, k], dt)
                    shift_pdf!(_h1, cell.lorenz[1, j, k], KS.vSpace.du[1, i, k], dt)
                    shift_pdf!(_h2, cell.lorenz[1, j, k], KS.vSpace.du[1, i, k], dt)
                end
            end
        end

        for k in axes(cell.h0, 4)
            for j in axes(cell.h0, 3)
                for i in axes(cell.h0, 1)
                    _h0 = @view cell.h0[i, :, j, k]
                    _h1 = @view cell.h1[i, :, j, k]
                    _h2 = @view cell.h2[i, :, j, k]

                    shift_pdf!(_h0, cell.lorenz[2, j, k], KS.vSpace.dv[i, 1, k], dt)
                    shift_pdf!(_h1, cell.lorenz[2, j, k], KS.vSpace.dv[i, 1, k], dt)
                    shift_pdf!(_h2, cell.lorenz[2, j, k], KS.vSpace.dv[i, 1, k], dt)
                end
            end
        end

        # force -> f^{n+1} : step 2
        for k in axes(cell.h1, 4), j in axes(cell.h1, 3)
            @. cell.h2[:, :, j, k] +=
                2.0 * dt * cell.lorenz[3, j, k] * cell.h1[:, :, j, k] +
                (dt * cell.lorenz[3, j, k])^2 * cell.h0[:, :, j, k]
            @. cell.h1[:, :, j, k] += dt * cell.lorenz[3, j, k] * cell.h0[:, :, j, k]
        end

        # source -> f^{n+1}
        tau = uq_aap_hs_collision_time(
            cell.prim,
            KS.gas.mi,
            KS.gas.ni,
            KS.gas.me,
            KS.gas.ne,
            KS.gas.Kn[1],
            uq,
        )

        # interspecies interaction
        prim = deepcopy(cell.prim)
        if isMHD == false
            for j in axes(prim, 2)
                prim[:, j, :] .= KitBase.aap_hs_prim(
                    cell.prim[:, j, :],
                    tau,
                    KS.gas.mi,
                    KS.gas.ni,
                    KS.gas.me,
                    KS.gas.ne,
                    KS.gas.Kn[1],
                )
            end
        end

        H0 = zeros(KS.vSpace.nu, KS.vSpace.nv, uq.op.quad.Nquad, 2)
        H1 = similar(H0)
        H2 = similar(H0)
        for k in axes(H0, 4), j in axes(H0, 3)
            H0[:, :, j, k] .=
                maxwellian(KS.vSpace.u[:, :, k], KS.vSpace.v[:, :, k], prim[:, j, k])
            @. H1[:, :, j, k] = H0[:, :, j, k] * prim[4, j, k]
            @. H2[:, :, j, k] =
                H0[:, :, j, k] * (prim[4, j, k]^2 + 1.0 / (2.0 * prim[5, j, k]))
        end

        # BGK term
        for k in axes(cell.h0, 4), j in axes(cell.h0, 3)
            @. cell.h0[:, :, j, k] =
                (cell.h0[:, :, j, k] + dt / tau[k] * H0[:, :, j, k]) / (1.0 + dt / tau[k])
            @. cell.h1[:, :, j, k] =
                (cell.h1[:, :, j, k] + dt / tau[k] * H1[:, :, j, k]) / (1.0 + dt / tau[k])
            @. cell.h2[:, :, j, k] =
                (cell.h2[:, :, j, k] + dt / tau[k] * H2[:, :, j, k]) / (1.0 + dt / tau[k])
        end

        #--- record residuals ---#
        @. RES += (w_old[:, 1, :] - cell.w[:, 1, :])^2
        @. AVG += abs(cell.w[:, 1, :])

    end

end


text:
"""
    Utility functions for MCMC
"""
############################################################################################
#Import
include("tuning/tuning.jl")
include("algorithm.jl")
include("postprocessing.jl")
############################################################################################
# Export


text:
export OffPolicy

using MacroTools: @forward

"""
    OffPolicy(π_target::P, π_behavior::B) -> OffPolicy{P,B}
"""
Base.@kwdef struct OffPolicy{P,B} <: AbstractPolicy
    π_target::P
    π_behavior::B
end

(π::OffPolicy)(env) = π.π_behavior(env)

@forward OffPolicy.π_behavior RLBase.get_priority, RLBase.get_prob


text:
module SfMatter

using Unitful

import PhysicalConstants.CODATA2018: k_B, g_n, R
import ..SfGeometry: volume

export Material, density, yield_strength, bouyancy, compression_energy

struct Material
	density::Unitful.Density
	yield_strength::Union{Missing, Unitful.Pressure}
end

Material(density::Unitful.Density) = Material(density, missing)

# empty definition for overriding
mass(::Nothing) = nothing

density(x) = mass(x) / volume(x) |> u"kg/m^3"
density(x::Unitful.Density) = x
density(m::Material) = m.density

yield_strength(m::Material) = m.yield_strength

"""
	bouyancy(ρ_body::Unitful.Density, ρ_fluid::Unitful.Density; g = g_n)

Acceleration due to bouyancy effects experienced by an object with density `ρ_body` in a fluid with density `ρ_fluid`, undergoing acceleration `g` (which defaults to Earth gravity).
"""
bouyancy(ρ_body::Unitful.Density, ρ_fluid::Unitful.Density; g = g_n) = g * ρ_fluid / ρ_body

"""
	bouyancy(ρ_body::Unitful.Density, ρ_fluid::Unitful.Density; g = g_n)

Acceleration due to bouyancy effects experienced by an object made of `body` in a fluid composed of `fluid`, undergoing acceleration `g` (which defaults to Earth gravity).
"""
bouyancy(body::Material, fluid::Material; g = g_n) = g * density(fluid) / density(body)

"""
    bouyancy(v_body::Unitful.Volume, ρ_fluid::Unitful.Density; g = g_n)

Bouyancy force felt by an object of volume `v_body` in a fluid of density `ρ_fluid`, undergoing acceleration `g` (which defaults to Earth gravity).
"""
bouyancy(v_body::Unitful.Volume, ρ_fluid::Unitful.Density; g = g_n) = g * ρ_fluid * v_body

"""
    bouyancy(v_body::Unitful.Volume, fluid; g = g_n)

Bouyancy force felt by an object of volume `v_body` in a fluid composed of `fluid`, undergoing acceleration `g` (which defaults to Earth gravity).
"""
bouyancy(v_body::Unitful.Volume, fluid; g = g_n) = g * density(fluid) * v_body

"""
    compression_energy(n, T, v_initial, v_final)

Work done to isothermally compress `n` moles of a substance from a volume of `v_initial` to `v_final`, maintaining a temperature of `T`.

If `v_final` is larger than `v_initial`, the result will be negative.

Ref: https://physics.stackexchange.com/a/37676/
"""
compression_energy(n, T, v_initial, v_final) = n * R * T * log(v_initial/ v_final)

end

max_stars_repo_name: icantbelieveitsnotphysics/SfPhysics

max_issues_repo_name: icantbelieveitsnotphysics/SfPhysics

max_forks_repo_name: icantbelieveitsnotphysics/SfPhysics


text:
# ------------------------------------------------------------------------------------------
# # Miscellaneous Algorithms
#
# Topics:
# * Lyapunov exponent from numerical data
# * Broomhead-King coordinates
# * Finding fixed points of any map of any order
# * Detecting chaos with GALI
# * Docstrings
#
#
# *warning: things will get really fast here*
#
# ---
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# # Lyapunov exponent from numerical data
# It is possible to compute a maximum Lyapunov exponent from numerical data, stored either
# as a `Dataset` or as a `Reconstruction`. This is done with the function:
#
# ```julia
# E = numericallyapunov(R::AbstractDataset, ks; kwargs...)
# ```
# It returns `E = [E(k) for k ∈ ks]`, where `E(k)` is the average logarithmic distance
# between states of a `neighborhood` for `k` steps (`k` must be integer).
#
# **For a more detailed description, please see the [documentation
# page](https://juliadynamics.github.io/DynamicalSystems.jl/latest/chaos/nlts/#numerical-
# lyapunov-exponent).**
#
#
# ------------------------------------------------------------------------------------------

using DynamicalSystems

# ------------------------------------------------------------------------------------------
# Let's compute the lyapunov exponent for reference and compare our numerically calculated
# lyapunov exponent to it!
# ------------------------------------------------------------------------------------------

ds = Systems.henon()
λ = lyapunov(ds, 100000)

# ------------------------------------------------------------------------------------------
# Below we create a trajectory, produce a `Reconstruction`, and then pass that
# `Reconstruction` to `numericallyapunov`:
# ------------------------------------------------------------------------------------------

data = trajectory(ds, 1000000)
x = data[:, 1] #fake measurements for the win!

R = Reconstruction(x, 2, 2)
ks = 1:20
E = numericallyapunov(R, ks)

# ------------------------------------------------------------------------------------------
# If the data experiences exponential separation of nearby states, then it must hold that
#
# $$
# E(k) \approx \lambda\times(k \Delta \! t ) + E(0)
# $$
#
# where $\Delta \! t$ is the time difference between successive data points (which is 1 for
# discrete systems).
#
# *This means that, $\lambda$, the slope of the plot E vs. k, is our numerically determined
# lyapunov exponent.*
# ------------------------------------------------------------------------------------------

using PyPlot
figure(figsize=(6,4))
plot(ks, E .- E[1])
xlabel("k"); ylabel("E");

# ------------------------------------------------------------------------------------------
# Above we see that the scaling behavior holds very nicely!
#
# We can use linear_regions to compute the slope automatically:
# ------------------------------------------------------------------------------------------

λ_numeric = linear_region(ks, E)[2]

println("λ - λ_numeric = $(abs(λ - λ_numeric))")

# ------------------------------------------------------------------------------------------
# `numericallyapunov` also takes a number of keyword arguments. You can decide:
#
# 1. Which states of the data should be used for reference states (by default *all* states
# are used).
# 2. How this logarithmic distance should be computed.
# 3. What type of `neighborhood` to use: fixed mass or fixed size.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# # Broomhead-King coordinates
#
# This alternative/improvement of the traditional delay coordinates `Reconstruction` can be
# a very powerful tool. An example where it shines is on noisy data where that noise creates
# the effect of superficial dimensions.
#
# To show this, we'll use `Distributions` to create random numbers.
# ------------------------------------------------------------------------------------------

using Distributions

ds = Systems.gissinger()
data = trajectory(ds, 1000.0, dt = 0.05)
x = data[:, 1]

L = length(x)
distrib = Normal(0, 0.1)
s = x .+ rand(distrib, L)

U, S = broomhead_king(s, 40)

figure(figsize= (10,6))
subplot(1,2,1)
plot(U[:, 1], U[:, 2])
title("Broomhead-King of s")

subplot(1,2,2)
R = Reconstruction(s, 2, 30)
plot(columns(R)...; color = "C3")
title("2D Reconstruction of s")

tight_layout();

# ------------------------------------------------------------------------------------------
# # Finding fixed points of maps
# Finding unstable (or stable) periodic orbits of a discrete mapping analytically rapidly
# becomes impossible for higher orders of fixed points. Fortunately there is an algorithm
# from Schmelcher & Diakonos for doing so numerically.
#
# The function that performs the algorithm is:
# ```julia
# periodicorbits(ds::DiscreteDynamicalSystem, o, ics, args..)
# ```
# where `o` is the order of fixed points to find and `ics` are the initial conditions to
# start searching from.
#
# This function can take a lot extra optional and keyword arguments and thus it is good
# practice to read the [documentation page](https://juliadynamics.github.io/DynamicalSystems
# .jl/latest/chaos/periodicity/#detecting-stable-and-unstable-periodic-orbits-of-maps) on
# it!!!
#
# I will now compute fixed points of order 3 & 4 of the standard map:
# ------------------------------------------------------------------------------------------

ds = Systems.standardmap()

xs = linspace(0, 2π, 21); ys = copy(xs)
ics = [SVector{2}(x,y) for x in xs for y in ys] # initial conditions container

FP3 = periodicorbits(ds, 3, ics) # order 3 periodic points

FP4 = periodicorbits(ds, 4, ics) # order 4 periodic points

figure(figsize = (6,4))

plot([s[1] for s in FP3], [s[2] for s in FP3], marker = "o", color = "C0", ls = "None")
plot([s[1] for s in FP4], [s[2] for s in FP4], marker = "D", color = "C1", ls = "None");

# ------------------------------------------------------------------------------------------
# Weeeeeeell this isn't really helpful as it doesn't show the phase-space of the standard
# map...
#
# Let's start by creating the phasespace plot and then overlaying these points!
# ------------------------------------------------------------------------------------------

iters = 500
dataset = trajectory(ds, iters)
for x in xs
    for y in ys
        append!(dataset, trajectory(ds, iters, SVector{2}(x, y)))
    end
end

# plot phasespace
figure(figsize = (6,4))
plot(dataset[:, 1], dataset[:, 2], ls = "None", ms= 0.8, marker = ".", color = "black", alpha = 0.5)
xlim(xs[1], xs[end])
ylim(ys[1], ys[end])

# plot fixed points
plot([s[1] for s in FP3], [s[2] for s in FP3], marker = "o", color = "C0", ls = "None")
plot([s[1] for s in FP4], [s[2] for s in FP4], marker = "D", color = "C1", ls = "None");

# ------------------------------------------------------------------------------------------
# * The function did **not** find all fixed points of order 3 or 4.
# * This happened because we did not set the other arguments of `periodicorbits`, and
# therefore a random value was chosen for them.
#
# For more, I am pointing you to the [documentation example](https://juliadynamics.github.io
# /DynamicalSystems.jl/latest/chaos/periodicity/#standard-map-example), where the following
# picture is computed:
#
# <img src="standardmap_fp.png" alt="Fixed points of the standard map" style="width: 800px;"
# align="left"/>
#
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# # Detecting Chaos using `gali`
# * In (for example) Hamiltonian systems, it is quite often the case that regular orbits
# (i.e. periodic) co-exist with chaotic orbits.
# * Being able to detect and distinguish chaotic from regular behavior is crucial in the
# study of dynamical systems.
# * Lyapunov exponents *can* be used for this task but they are not efficient.
#
# * One of the newest methods to do this is "GALI", generalized alignment index method.
#
# * GALI uses the fact that deviation vectors (that live on tangent space) tend to align for
# chaotic motion, while they stay "not-aligned" for regular motion.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ---
# Here I will only show a simple example of the standard map, which is a system with **mixed
# phase space**.
#
# For *regular* orbits, `gali` stays approximately constant. For *chaotic* orbits, `gali`
# decays exponentially.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# We'll initialize with a random chaotic initial condition,
# ------------------------------------------------------------------------------------------

sm = Systems.standardmap(;k = 1.0)

# ------------------------------------------------------------------------------------------
# get a trajectory,
# ------------------------------------------------------------------------------------------

tr = trajectory(sm, 10000)

# ------------------------------------------------------------------------------------------
# and then plot it:
# ------------------------------------------------------------------------------------------

using PyPlot; figure(figsize=(6,4))
plot(tr[:, 1], tr[:, 2], ls = "None", marker = "o", ms = 0.2);

# ------------------------------------------------------------------------------------------
# To compute `gali`, I use the following call signature:
# ```
# gali(ds::DynamicalSystem, tfinal, k::Int; kwargs...) -> g, t
# ```
# where `k` is the order of GALI I want and `tfinal` specifies for how much time to evolve
# the system and the deviation vectors.
#
# Let's compute the gali_2 for this trajectory:
# ------------------------------------------------------------------------------------------

g, t = gali(sm, 1000, 2)
g

# ------------------------------------------------------------------------------------------
# * `g` quickly decays to zero for a chaotic orbit!
# ------------------------------------------------------------------------------------------

# Initialize with random regular initial condition:
sm = Systems.standardmap([π + 0.01*rand(), 0.01*rand()];k = 1.0)
# Get trajectory:
tr = trajectory(sm, 10000)
# Plot it:
using PyPlot; figure(figsize=(4,2))
plot(tr[:, 1], tr[:, 2], ls = "None", marker = "o", ms = 0.2);

# gali with regular initial condition:
g, t = gali(sm, 1000, 2)
g

# ------------------------------------------------------------------------------------------
# * `g` stays constant for a regular orbit!
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# # Docstrings
# ------------------------------------------------------------------------------------------

?numericallyapunov

?broomhead_king

?periodicorbits

?lambdamatrix

?gali

max_stars_repo_name: grenkoca/JuliaTutorials

max_issues_repo_name: grenkoca/JuliaTutorials

max_forks_repo_name: grenkoca/JuliaTutorials


text:
module Tweets

using App, SearchLight, Channels

export Tweet

mutable struct Tweet <: AbstractModel
  ### internals
  _table_name::String
  _id::String

  ### fields
  id::Nullable{SearchLight.DbId}
  tweet_id::String
  text::String
  user_name::String
  screen_name::String
  avatar_url::String
  favorite_count::Int
  retweet_count::Int
  polarity::Int
  subjective::Bool

  ### extra values
  search_hash::String

  ### callbacks
  after_save::Function

  ### constructor
  Tweet(;
    id = Nullable{SearchLight.DbId}(),
    tweet_id = "",
    text = "",
    user_name = "",
    screen_name = "",
    avatar_url = "",
    favorite_count = 0,
    retweet_count = 0,
    polarity = -100,
    subjective = false,

    search_hash = "",

    after_save = (m::Tweet) -> Channels.broadcast(["tweets-" * m.search_hash], "tweet", SearchLight.to_dict(m))

  ) = new("tweets", "id", id, tweet_id, text, user_name, screen_name, avatar_url,
          favorite_count, retweet_count, polarity, subjective, search_hash, after_save)
end

const EMPTY_DATA = ([], [0, 0, 0], [0, 0], "[]", "[]")

function tweets_stats(tweets::Vector{Tweet}) :: Tuple{Vector{Int},Vector{Int}}
  favorites_count = Int[]
  retweets_count = Int[]

  for tweet in tweets
    push!(favorites_count, tweet.favorite_count)
    push!(retweets_count, tweet.retweet_count)
  end

  favorites_count, retweets_count
end

function tweets_aggregates(filters) :: Tuple{Vector{Int},Vector{Int}}
  result = SearchLight.query_raw("
    SELECT
      SUM((polarity = -1)::int) AS negative,
      SUM((polarity = 0)::int) AS neutral,
      SUM((polarity = 1)::int) AS positive,
      SUM((subjective = TRUE)::int) AS subjective,
      SUM((subjective = FALSE)::int) AS objective
    FROM tweets
    WHERE $(filters[1]) $(filters[2] |> string )") |> first

  map(result) do (r)
    r == Union{} ? 0 : r
  end

  [result[1], result[2], result[3]], [result[4], result[5]]
end

end

max_stars_repo_name: essenciary/genie-demo-tweet-stats

max_issues_repo_name: essenciary/genie-demo-tweet-stats

max_forks_repo_name: essenciary/genie-demo-tweet-stats